有限元分析程序设计

本文格式为Word版，下载可任意编辑——有限元分析程序设计结构有限元分析程序设计

绪论

§0.1开设―有限元程序设计‖课程的意义和目的§0.2课程特点§0.3课程安排§0.4课程要求§0.5基本方法复习

$0.1意义和目的

1.有限元数值分析技术本身要求工程设计研究人员把握

1).有限元数值分析技术的完善标志着现代计算力学的真正成熟和实用化,已在各种力学中得到了广泛的应用。譬如：，已杨为工程结构分析中最得以收敛的技术手段，现代功用大致有：

a).现代结构论证。对结构设计从内力，位移等方面进行优劣评定，从而进

行结构优化设计。

b)可取代部份试验，局部试验+有限元分析，是现代工程设计研究方法的一大

特点。

c)结构的各种功能分析（疲乏断裂，可靠性分析等）都以有限元分析工具作为

核心的计算工具。

2).有限元数值分析本身包括着理论+技术实现（本身功用所绝定的）

有限元数值分析本身包括着泛函理论+分片插值函数+程序设计2.有限元分析的技术实现（近十佘年的事）更依靠于计算机程序设计

有限元分析的技术取得的巨大的成就，从某种意义上说，得益于计算机硬件技术的发展和程序设计技术的发展，这两者的依靠性在当代表现得更加突出。（如可视化技术）3.从学习的角度，不仅要学习理论，而且要从程序设计设计角度对这些理论的技术实现有

一个深入的了解，应当致力于把握这些技术实现能力，从而开发它，发展它。（理论本身还有待于进一步完美相应的程序设计必需去开发）

4.程序设计不仅是实现有限元数值分析的工具和桥梁，而且在以下诸方面也有意义：

1).精通基本概念，深化理论认识；

2).锻炼实际工程分析，实际动手的能力；3).获得以后工作中必备的工具。（作业+老师给元素库）

目的：通过陈述有限元程序设计的技术与技巧，便能达到自编自读的能力。§0.2课程特点

总描述：理论+算法+数据结构（程序设计的意义）

理论：有限元算法，构造，步骤，解的等外性，收敛性，稳定性，误差分析

算法；指求解过程的技术方法，含两方面的含义；a.有限元数值分析算法，b,与数据结构有关的算法（总刚稀疏存贮，提取，节点优化编号等）

数据结构：指各向量矩阵存贮管理与实现，辅助管理结构（指针，数据记录等）具体特点：

理论性强：能量泛函理论+有限元构造算法+数据结构构造算法内容繁杂：理论方法+技术方法+技术技巧

技巧性强：排序，管理结构（指针生成，整型运算等）

§0.3课程安排

①.单元刚度矩阵及元素设计（单元刚阵算法，杆梁平面分析，板弯非协调元等）②.总刚的形式及程序设计（单刚提前准备，技术繁杂）

③.l边界条件及程序设计（等效荷载计算，位移边界条件置入，多工况的对称性）④.总刚线性方程组求解（LDLT分解，分块算法，子结构算法，波前法）⑤．单元应力计算+应力处理与改善。

⑥.数据处理（数据分类，压缩存贮，节点排序方法）

⑦,变带宽存贮的优化理论，图论的理论，有限元的图结构，存贮管理复核）⑧.有限元议程全稀疏管理与求解策略。

说明：仅线性部份，复材，接触，弹塑性等不包括，基本部份。

实践性作业安排：

1.作业：总的结构管理程序+子功能模块的编程，一个题的计算实践2.送有限元元素库。§0.4课程要求1.先行要求

2.作业要求（计算机编程得出正确结果）3.课程校核要求§0.5基本方法复习

0.5.1结构化程序设计方法0.5.2有限元分析方法回想0.5.3Fortran语言回想0.5.4结构化程序设计

1.基本结构：构成一个问题从输入到求解输出的基本程序形式：

Imput——→Process——→Output(输入)（处理）（输出）

三种基本形式：

a).树形结构（顺序执行结构）

Imput——→p1→……….——→pn→Output(多用于程序结构call1,call2……..)b).t选择结构（条件执行结构）YInputBeginNP1OutputP2内部算法实现:IF…..THEN;SWHICH,CASEc).循环结构(反复执行结构)特点:结构特征简单明白,易读,易调试`.尽量少用GOTO语句2.整体结构(算法语言系统结构)积木式(Fortran):每个设定的功能分析团体的一个模块,每个模块又称作整体结构的素材,主模块象积木一样堆积.语体不联系,但有通讯方法沟通模块间信息,各种模块有各自的特征语体,mainprogame……,subroutinesub….

嵌套式(Pascall):主模块与子模块相互嵌套,各模块的特征,语体一致proceduremainprocefuressub1……..

proceduresub2……….Endsub1………endsub2…….Endmain

函数式(c语言):主要特点是功能模块作为库函数调用,需用时在库内调用,每一个函数有表征语句,这种语体接近自然思维,而且对系统资源的调配应用更完善.面向对象的程序设计:实施过程的可视化+控制性

3.结构化程序设计方法

a).Top—Down(自上面下),系统性强,选择性强.

b).CriticalComponentFirst(关键部份优先),先抓主要矛盾,分清重,缓,急.c).独立调试,总体联调,(软件设计的社会化作业).4.程序设计要点

a).自觉有意识地设计一个良好的程序结构,做到:易读,易懂,易管理,易修改,易发展.b).做到规律明了,说明完整.

c).要有工艺设计概念有框图,有步骤.5).结构化程序设计原则

a).尽可能通用性好（适应各种规模的复题，？的扩大依据程序设计指标而定）b).整体精炼，明了；避免GOTO。c).省机时，省存贮，计算精度高，（算法上下功夫，要理论分析加技巧）d).输入数据少，格式简单。

e).输出结果简明，忌讳打印过多（与具体调试过程不一样）。f).易读易维护，易发展。§0.5.2.有限无方法求解过程回想一.力学模型的分级管理

有限无程序对力学模型的数据按

一级：结构级（有点广义，不仅指具体结构，也指模型题目的规模）

`二级：单元级三级：节点级2.基本关系离散化由结构单元节点（集合）3.描述参数A).节点描述参数(1).节点位置(总体坐标系下的坐标).

(2).节点局部坐标(按节点的约束方向制定的特别坐标系x＇，y＇，z＇,v如斜支撑)(3).节点的性质(自由,固定,指定位移,附属其它节点).(4).节点力:(Fx,Fy,Fz,Mx,My,Mz)(5).节点位移:(u,v,w,θx,θy,θz)B.单元描述参数

(1).材料特性参数不清E,G,γ→[D]

(2).节点的几何刚度参数(即面积A,板厚H,梁抗弯模量I)(3).单元的局部坐标.(用于应力分析等,如图形曲面)(4).单元的节点编号

(5).单元的几何矩阵营(节点变形与应力关系矩阵)(6).单元刚度矩阵[K]

(7).单元的应力,应变向量,(有限元分析多用向量,而不用矩阵(张量))

??11?????22??????11?12?13??33????21?22?23????12????31?32?33???????13?????23???11??????22???11????1?33????21??12??2????1?31?13??????2?23?1?122?221?3221??132?1??23?2?1??332??结构描述参数

单元总数,节点总数,单元娄型总数,结构材料种娄数,节点自由度数(控制题目规模)二.基本公式系统

1.单元刚度计算公式

[K]e??[B]T[D][B]dvVc2.单元刚阵组合[K]=ΣATKA3.单元节点荷载计算

{P}e??[N]T{P}dsSe4.节点荷载组装:

?P???APe5.位移约束关系:

????d???d???6.总刚方程解:

?K??????P???????K?e?P?7.应变计算:

?????B????8.应力计算:

?????D??????D??B????9.支撑反力计算:

Ri?Kiidi三.有限元分析的模块组织.输入原始参数?计算总刚规模形成单元刚阵形成总刚方程坐标变换单刚向总刚投放输入边界条件(对称条件)形成单元荷载阵形成各荷载工况的节点荷载阵向总节点荷载阵投放总刚分解回代求出位移及输出计算应变、应力

调整单元位移列阵调整几何、弹性矩阵

四．结构分析的原始输入数据

1．题目规模`节点数目：NNP单元数目：NE2．节点数据单元人坐标：XE（NNP，3）3．单元数据单元节点编号：ME（NE，3）、ME（NE，2）材料特性：E、N单元几何参数：I、RI（惯矩）

4．荷载数据外荷载作用点，坐标及大小：PA（NNP，1）§0.5.3Fortran语言回想

1．子模块（子程序）subroutine

a.特点：独立性强，只要输入输出接口，象一个黑匣子，与外界无关。b.作用：完成一个独立的功能（求应力，矩阵分解，投放等）

c.格式：subroutinefunction(ip1,ip2,rp1,rp2,io1,io2,ro1,ro2).(其中ip1,ip2,rp1,rp2,是输入形参，io1,io2.rp1,rp2是输出形参)2．数据传递形式

1）.COMMON公共块语句传递，（公共块的内容不能作为形参）

a.公共块分为无名公共块和有名公共块

b.公共块的参数不能作为子程序的参数出现，

c.公共块名一致，其内容在不同公共块中可以标志符不同（但其长度应一致）d.通用原始数据放入公共块（作为实参错误率大）e.尽可能不放数值，安眠组一般可作成可调长度f.格式Common/comm/…….Subroutinefun()

Common/comm/……..2).形参————实参对应a.实参不能开拓存贮单元，子程序内定义语句中的形参数组由主程序定义，在子程序中仅形式定义（即仅说明是数组，因而大小无所谓）b.格式：DimensionRP(1000),RO(1000)…….

Callsub1(RP,RO)………END

DIMENTIONIBANK

SUBROUTINESUB1(RP,RO,NE)DIMENTIONRP(1),RO(1),SP(50)

DIMENTIONRP(NE,1),RO(NE,1)(形参的动态定义，实参不能)3）.数组长度自动调整方法。PROGRAMMAIN

INPLICITREAL*8(A—H,OZ)CHARACTAR*20TRCOMMON/COMM/….

DIMENTIONIBANK(),RBANK(),IP1(),IP2()IP1(1)=…IP(N)=…..IP2(1)=……IP2(N)=…..CALLSUB1(IBANK(IP1(1)),IBANK(IP1(N)),RBANK(IP2((1)),….)

…..END

SUBROUTINESUB1(II1,IO2,….RI1,….RO1….NE)DIMENTIONII1(1),IO2(NE,1),RI1(1),RO1(1)

*.公共块分无名与有名公共块：公共块参数不能子程序形参中出现。公共块一致，其内容在不同程序中可以标识符不同。（但长度要一致）*.公共块使用原则：

1)通用原始数据代入公共块，作为实参出错率大2)尽可能不放数组，数组要作成可调的（自动可调）。

数组长度自动调整方法。

§1.单元刚阵几元素设计

主讲内容：工程常用元素的单元刚阵计算与编程。§1.1杆、梁的单元刚阵§1.2平面问题的三角形、矩形单元刚阵§1.3板弯问题的三角形、矩形单元刚阵§1.4非、拟协调之介绍及分片检验方法

§1.0本讲内容序

§1.0.1单元刚阵在总刚中的地位和作用

由有限元分析过程可以知道，最主要的两步骤：

1)进行单元的特性分析，建立单元级的刚度方程（平衡方程），得到单元的刚

度矩阵；

[K]e{?}e?{R}e

注：[K]e中各元素的位置与{?}e和{R}e要有严格的对应关系

2)进行结构整体分析，集合所有单元的刚度方程建立全结构的刚度方程，设得

到全结构的刚度矩阵。即

[K]??[K]e??ATKeA

A为布尔矩阵

总刚[K]e由单元刚阵按节点编号顺序投放而成，故单元刚阵是整个有限元方

法的第一步；是总装刚阵的原材料。

解释：单元节点平衡；结构总刚对号和迭加的原理§1.0.2单元刚阵形式

1.基本公式：

[K]e??[B]T[D][B]dV

V（1.1）

[B]成为几何矩阵，作用是把应变与位移联系起来反映了单元的形变特性与

性质。固而是坐标的函数，规模与节点数目（自由度数目）有关，不同类型的单元有不同的[B]。（依据设定的位移形函形式与节点位移的性质）。[B]对常应元是常数。

只取决于材料特性E，与坐标无关，对于平面应力问题：[D]为弹性矩阵，?，

??1?0?E???1[D]?0?（1.2）21???1???00??2??对于平面应变问题，将上式E?2.形成方法

1)直接赋值。有了积分出来的[K]e各元素显式，即可逐项（元素）直接赋值。在赋值时，有数字变化规律的应按特点赋值，以缩短程序长度，

提高执行效率。

2)矩阵乘法。用矩阵乘法进行赋值。

§1.1杆、梁的单元刚阵

§1.1.1杆元

对于各种变截面等轴力杆其单元刚阵显示：

y

?E;??代入上式即得。

1??(1??)(1??)EAeq?1?1?（1.3）[K]???L??11?ex’uj(Pj)y’Aeq为等效截面面积。

对于不同的截面变化状况可由材料力学ui(Pi)方法计算得到。

*对于等截面杆，可设定形状函数，利用公式（1.1）计算，即由：

x???Ni?i

i?12

x'N1?1?,

Lx'N2?

L

???dNidx'?i??1?1??11????L??2?(1.4)

[K]e?BTDBV?BTEBAL?EA?1?1???L??11?(1.5)

*对于变截面杆，采用假设形状函数的方法不妥（误差大），因此，常用材料

力学推倒方法：即公式：

pdxPLdxd??????(等轴力)（1.6）

EAE0A相对于1点：

∵???1??2∴P1?

?

P2?L0E(?1??2)?1dxA?L0??1?E?1?1???1??2?dxA同理，相对于2点：

?于是得：

L0?1?E??11????1??2?dxA[K]e?E?1?1???11?1?dx?A?令：Aeq?L（1.6）

0?L0L故得：（1.3）式1dxA

思考题：变截不均的弯轴力杆的刚阵形式如何？

对于满足杆元截面积随坐标x'单调变化的函数，Aeq计算见P28～30

*形成方法可采用直接赋值语句完成见P30FAKE1§1.1.2梁元

应用公式（1.1）进行直接积分，得到显式，然后用上赋值语句编程。1.梁元形及刚阵显式

工程上使用的梁元一般为复合受力状态，包括轴向、扭转及各平面内的弯曲，因此对于一般的三维空间中的梁元节点位移包括：线性位移u、v、

?及转角位移?x、?y、?z；节点力位包括：三个方向的力Px、Py、Pz和三个方向的弯矩。

y

z

故梁元每个节点有6个自由度，即：

y’x’z’x?u,v,?;?iiixi,?yi,?zi;uj,vj,?j;?xj,?yj,?zj?

??Pxi,Pyi,Pzi;Mxi,Myi,Mzi;Pxj,Pyj,Pzj;Mxj,Myj,Mzj?1)对于轴向拉伸与扭转可采用线性形函描述

?uxi????u??c0?c1x??N10N20???xi?即：?????????u??c?cx0N0N3?12??xj??x??2????xj??xxN1?1?N2?

LL?uxi??du????dx?1??1010???xi??BU??应变定义?????u?d?x?0?101L???xj?????dx????xj??0??EAD???0GIx??(1.7)

（1.8）

Ix???2dA

A截面极惯矩（1.9）

?K?e??0BTDBdx

L注意：这里的体积分、面积分为常数包括在弹性阵内

?EAL??GIx???L???GIx??L??K?e1?EA?L?GIx?L???EA??L?GI??xL?EAL(1.10)

2)在xy平面内的弯曲，不计剪切，节点位移为?vi?zivj?zj?；节点力为

?PyiMziPyjMzj?形的选取是依据弯曲理论（材料力学）的基本假设的

两直法线公设所确定即：

EIzd4uydx4?0

（1.11）

进行积分，可得三次挠度曲线

v?c0?c1x?c2x2?c3x3

（1.12）

代入节点位移系统可得到形状函数

?vi?????zi?v??v1v2v3v4???

?vj????zj??xxv1?1?3()2?2()3LLx2x3v2?x?2?LL2

x2x3v3?3()?2()LLx2x3v4???LL2（1.13）

（1.14）

*应变定义：梁的曲率定义为广义定义应变，梁的弯矩定义为广义应力，即：

d2v??2

dx??Mz（1.15）

按梁的挠度方程，可有广义的应力应变关系

??E??(1??M)EIz(1.16)

其中

E?EIz

?vi?????zi?B4????vj????zj??代入（1.13式），得：

??B?u???B1B2B3（1.17）

612?xL2L346B2???2xLL

612B3?2?3xLL26B4???2xLLB1??L（1.18）

由式（1.1），于是：

?K?e??0BTEIzBdx?K?e2（1.19）

对?12??称?EIz?6L4L2??3?L??12?6L12?22?6L2L?6L4L??（1.20）

3)

??在yz平面内的弯曲，节点位移与力系统为：

?P形函推导过程类似。刚阵显式：

e?K?3i?yi?j?yj?MziPzjziMzj?

对?12??EIy?6L4L2称??3?

???12?6L12L?22?6L2L?6L4L??（1.21）

4)综合上述三种单独刚度，可以得空间梁元素在局部坐标系下复合状态的刚阵

?K?e???K?ie

i显示见P42～43间，即书公式（2－24），（注意节点位移、节点力的排序）。

2.形成梁元刚阵的子程序设计

方法：直接赋值。见P43～44FAKE4

§1.2平面问题的三角元、矩形单元刚阵

§1.2.1平面应力三角元。（常应变元）1.公式。

??1?0?E???1?D??0?21???1????00?2??(1.22)

?bi?B??1?0?2???cibi?yi?ymci?xm?xj

??0cibibj0cj0cjbjbm0cm0??cm?bm??(1.23)

y’其中：在板平面内

Qmbj?ym?yj

转换

bm?yi?yj1

3

PiTjR

1(bicj?bjci)22

x’2.形成B阵的轮换循环赋值方法。

I=1,3

I1=I+1-3*(I/3)I2=I+2-3*((I+1)/3)KI=ME(K,I1)KJ=ME(K,I2)

B(1,2*I-1)=XE(KI,2)-XE(KJ,2)B的第1行1,3,5之元素B(2,2*I)=XE(KJ,1)-XE(KI,1)B的第2行2,4,6之元素B(3,2*I-1)=B(2,2*I)PP34B(3,2*I)=B(1,2*I-1)未完当平面元处于空间位置时，（即节点坐标为三维数据）。平面内的几何参数需要计算，计算公式：1）PQ长度：y’dPQ?(XQ?XP)2?(YQ?YP)2?(ZQ?ZP)2QmTRiPj2）PQ的方向余弦：lPQ?(XQ?XP)dPQmPQ?(YQ?YP)dPQnPQ?(ZQ?ZP)dPQ3）PT长度（按线段投影公式）x’dQT?lPQ(XR?XP)?mPQ(YR?YP)?nPQ(ZR?ZP)

4）TR长度：

222dTR?dPR?dPT?(XR?XP)?(YR?YP)?(ZR?ZP)?dPT

5）QT长度：

dQT?dPQ?dTP

6）TR的方向余弦

lTR?(XR?XT)dTR?XR?(XP?dPT?lPQ)dTRmTR?(YR?YT)dTRnTR?(ZR?ZT)dTR???Y??ZRR?(YP?dPT?mPQ)dTR?(ZP?dPT?nPQTR?)?d?7）?面积

S??1dPQ?dTR28）刚度阵中的参数

eB?E24?(1?u)9）在?B?中各长度的取值

x31?0

y31?dPQ

记号规则：

x32?dTRx21?dTR

y32?dTQy21?dPT

xij?xi?xjyij?yi?yj

利用上面的子程序，可直接赋值形成局部坐标下的?B?阵。

3.局部系下刚阵计算与编程。

?K?

e?T??EC??BDBdv?BDBV?BDB

4(1?u)2?TT（1.24）

其中C为板厚

利用矩阵乘法可直接运算得到?K?间见P35

eP43m

§1.2.2平面应力四节点单元刚度1．单元平面内局部下形状函数。

节点位移和节点力系统：

i12j???e??ui?P?e??PxiviPyiujPxjvjumPyjvmPxmupPymvpPxpPyp

??形状函数：

y1x(1?)(1?)4aby1xvj?(1?)(1?)4ab

y1xvm?(1?)(1?)4aby1xvp?(1?)(1?)4abvi?0?v1?y?v1?x?v2?x0?v2?y0?v2?y?v2?x（1.25）

??v1???x?B???0???v?1???y?v3?x0?v3?y0?v3?y?v3?x?v4?x0?v4?y?0???v4??y??v4???x??（1.26）

0b?y0b?y0?(b?y)??(b?y)???B??1?0?(a?x)0?(a?x)0a?xa?x?4ab??b?ya?xb?ya?x?(b?y)???(a?x)?(b?y)?(a?x)?（1.27）

显见得?B?不再是常数阵，说明应变在平面内线性变化。2.矩阵平面应力刚阵形成

公式：?Ke????B??D??B?ds

Ts由于不再是常数阵，只能手工积分每一元素，设得到刚阵显式，然后用赋值方法编程。

?K?e见公式PP38-39

（2-19）

刚度形成子程序FAKE3见PP39-40

思考题：变轴力、变截面杆单元刚阵的形成方法。（用柔度方法）

作业：（1）阅读《有限元》PP27-60内容

（2）有一空间梁元，

截面尺寸：

E?210000?Pa

??0.3?25

长度1m.

编程计算单元刚度矩阵。?50

3．程序框图:

单元节点行循环，IR=1,NPF单元节点列循环IS=1，NPFKO(IR)

Ｋdd=Ｋgd=Ｋgg=

???[Bjjkd]T[D][Bd]|J|HiHjHk]T[D][Bd]|J|HiHjHk

???[Bjjkjjkg???[BT][D][Bg]|J|HiHjHkg②、Guass消元法消去中间节点，得聚缩刚阵；

~[Kdd]??d????d?

~~[Kgd]??d??[Kgg]?g?g??????

????[K?1~][Kgd]??d?gg③、程序系统不再介绍。

§2总体刚阵的形成及程序设计

从有限元的计算过程可知，形成单刚、集合总刚、臵边界条件、求等效载荷列阵后，才能得到位移，应变及应力，课程正是按此工作进行，最终达到我们的目的。

§2.1总刚的特征及其组装时所需要效能的问题1总刚特性及投放过程

总体结构刚阵的形成是各单元刚阵按节点序排列的集合矩阵（对号入座）。①、对称性：依据Betti互易原理（讲：在结构体的两点1、2作用单位位移

11位移，（１）先1点作用,再2点作用?1?K11?K22?K12

2211(2)先2点作用,再1点作用?2?K22?K11?K21

22?1??2?K12?K21

②、带状性：（刚阵元素较集中在对角线的位臵）一个节点自由度对应总刚的

一行（一个节点不可能和结构的所有节点相连），这一行总是和相关的节点刚阵元素相关，