第二章弹性力学问题有限元方法的一般原理和表达格式.docx

第二章弹性力学问题有限元方法的一般原理和表达格式.docx

《第二章弹性力学问题有限元方法的一般原理和表达格式.docx》由会员分享，可在线阅读，更多相关《第二章弹性力学问题有限元方法的一般原理和表达格式.docx（39页珍藏版）》请在冰豆网上搜索。

第二章弹性力学问题有限元方法的一般原理和表达格式

第二章弹性力学问题有限元方法的一般原理和表达格式

2.1引言

本章将讨论通过弹性力学变分原理建立弹性力学问题有限元法列式的基本步骤。

最小位能原理的未知场变量是位移，以结点位移为基本未知量，并以最小位能为基础建立的有限单元位移元。

它是有限元方法中应用最普遍的单元。

对于一个力学或物理问题，在建立其数学模型以后，用有限元方法对它进行分析的首要步骤是选择单元形式。

平面问题三结点三角形单元是有限元方法最早采用，而且至今仍经常采用的单元形式。

我们将以此作为典型，讨论如何应用广义坐标建立单元位移模式与位移插值函数，以及如何根据最小位能原理建立有限元求解方程的原理、方法与步骤，并进而导出弹性力学问题有限元方法的一般列式。

2.2弹性力学平面问题的有限元列式

2.2.1单元位移模式及插值函数

典型的三结点三角形单元结点编码为每个结点有两个位移分量，如图2.2所示。

每个结点的位移可用位移矢量勺表示，即

］亿丿‘加）

每个单元有6个结点位移分量（称为6个自由度），于是单元结点的位移向量可表示为

■为单元结点位務列阵。

1.单元的位移模式和广义坐标

在有限元方法中单元的位移模式，是指在单元内位移的插值函数，其一般形式采用多项式作为近似函数，因为多项式运算简单，并且随着项数的增多，可以逼近任何一段光滑的函数曲线。

假设3结点三角形单元位移模式选取一次多项式

卩\+02兀+戸3『

（2.2.1）

v=pA+p5x+pby

16

它的矩阵形式是

u=（/）p（2.2.2）

其中

LvJL°卩」

0=[1xy]

0=口02AAAAJ

由于三个结点也在单元内，满足位移模式，于是得

%=屈+02兀+03儿

itj=0i+02勺+03儿（223）

%=01+02几+03九上式是关于久02仏的线性方程组。

久0异3是待定常数，也称为广义坐标。

它可由（2.2.3）式求岀。

上式的系数行列式是

1m片

D=1=2A（2.2.4）

1几几

上式中当i,j,m的编号顺序与坐标顺序一致时（此处为逆时针），A值为正，其大小为三角形面积，因此为了方便单元的编号一般与坐标的转向一致。

兀=2（血件+如仃+»%）

（2.2.5）

同理，y也有三个线性方程组匕=04+0妙+06儿

b=04+05厂+亿儿

f=04+05心+06X”

由上面方程组可求得0'05，几，即

V76力%Xmym

*■（"必fj+5%）

V,

=右@M+昭+bmvm）（2.2.6）

V.

土（們+仍+讣）

在（2.2.5）和

（2.2.6）式中

=Xjyj心儿

=9•儿一佔

=yj-ym（ij响（227）

1X：

上式下标（i.jjn）轮换，可得u'尸Cj及％饥，5o

2.位移插值函数

将求得的广义坐标人,02,…,06代入（221）式中，则位移函数可表示为结点位

移的函数，即

u=Niui+NjUj+N,nu,n

V=Nivi+NjVj+Nnlvm（2.2.8）

其中

Ni=丄（心+bix+ciy）（iJjfJ）（229）

2A

称为单元的插值函数或形函数，这里它是龙,y的一次函数，其中

G，Z?

i，Cj,及d”?

—,C”r是常数，由表达式可知，它完全由单元的大小和

方位确定，一旦单元确定了，这些常数也完全确定。

（2.2.9）中的单元面积可表示为

A=£D=£（©+aj+勺”）=|（叽-bjC）（2.2.10）

（2.2.8）式的矩阵形式是

「叮

Nj

Vr

Nj

知

丿,“_

=[^,叫IN,,]

久_

=[[7V]f[N].[N]m}xe=[N]ae（2211）

[N]称为插值函数矩阵或形函数矩阵。

插值函数具有如下性质：

（1）在结点上插值函数的值有

1当/=「

N）（x.，儿）=吊=|0当）知⑺几也）（2212）

也就是说在i结点上R=l,在j,m点上弘=0。

这一点与位移插值函数的表达式一致。

（2）在单元中任一点各插值函数之和应等于1,即

+N,+Nm=1（2.2.3）

以上结论可以由下面得到解释，若单元发生刚体位移，如x方向有刚体位移％,显然单元内任何一点的位移都应该为5,当然结点位移也为心，即ut=iij=»由（2.2.8）式得

U=Ng+NjUj+Nmum=（Nj+Nj+Nm）uQ=w（）

于是有M+y+M”=i。

若插值函数不满足此要求，在不能反映单元的刚体位移，用以求解必然得不到正确的结果。

单元内1的性质称为规一性。

（3）对于现在的单元，插值函数是线性的，在单元内部及单元的边界上位移也是线性的，可由结点上的位移值唯一确定。

由于相邻单元公共结点的位移是相等的，因此保证了相邻单元公共边界上位移的连续性。

3.应变矩阵和应力矩阵

⑴应变

确定了单元位移后，可以很方便地利用儿何方程和物理方程求得单元的应变和应力。

在（1.4.21）式的儿何方程中，位移用（2.2.11）式代入，得到单元应变为

理瓷Ji

0

=hl,［禺［叫”=瞰（2.2」4）

B称为应变矩阵，L是平面的微分算子。

应变矩阵的分块矩阵［B］,是

竺0

dx

（I,j.m）（2.2.15）

0竺dydN.dN.

dydx

对（2.2.9）式求导得

（2.2.16）

dM1fcNi1

=—b：

9=—c；dx2A62A

代入（2.2.15）式得

b,0

（/,j.ni）（2.2.17）

3结点三角形单元的应变矩阵为

B=[[B]r.[禺-0

b、0°Cj5勺

5

b；

（2.2.18）

式中»ct,bjq及乞,c,”由单元结点坐标确定，当单元确定了，这些常数也完全确定，因此B市常量矩阵。

而应变£=Ba",当单元结点位移确定后，单元内的应变都是常量，因此3结点三角形单元为常应变单元。

在应变变化较大的区域，单元划分应适当加密，否则不能反映应变的真实变化而导致较大误差。

⑵应力

单元应力可以根据物理方程求得，而对于平面问题

E'

6=匚討乙+心、）

E'

^■=—7

1—V

E9

Tx-=2（1+胪＞'

（平面应力）

（平面应力）

（平面应变）

（平面应变）

因此，对于平面应力问题，弹性矩阵D为

0

1-v2

0

1-v

由（2214）可知，c=Bae,从而可得

a=Ds=DBa（=Sae（2.2.19）

上式中S写成分块矩阵的形式，则

S=DB=D^B]t[B]j[S].

（2.2.20）

S称为应力矩阵。

写成分块矩阵的形式

c>°

L222d

对于平面应变问题：

Et—占一，—,对于具体计算时只要对不同的l-v21-v

赋值，就得到不同的平面问题。

从而可以到，三大物理参量，都可以用单元结点位移向量表示：

U

=N2

=Bae

=DBac=Sae

由于N,B,S都是已知的矩阵，只要求得则单元内的位移、应变和应力就可以就得，问题是：

如何求结点位移向量

2.2.2利用最小位能原理建立有限元方程

最小位能原理的泛函总位能口卩的表达式如（1.4.51）

在平面问题中最小位能原理的总位能与（1.4.51）式的形式一致。

现改为矩阵形式：

（2.2.24）

其中，t是二维体的厚度；f是体积力；T是表面力。

积分域也作了相应的改变。

对于离散模型，系统位能是各单元位能的利将（2211）及（2214）式代入（2.2.24）

式，即得到离散模型的总位能为

np=Snp=工巴丄占

BlDBtdxdyae^

-工LLNrftdxd^~刃”*NFds）（2.2.5）

将结构总位能的各项矩阵表达成各个单元总位能的各对应项矩阵之和，隐含着要求单元各项矩阵的阶数和结构各项矩阵的结束（自山度数）相同。

Ke=\BTDBtdxdyP；=JNTftdxdy

P；=\NFdS

Js；

（2.2.28）

Ke和Pe称为单元的刚度矩阵和单元等效结点载荷列阵。

于是

口厂工丐工/p（2.2.25b）

eee

由于每个单元的结点不同，相应的结点位移不同，为了统一起见，引入一个单元结点自山度和整体结构结点自山度的转换矩阵G,从而将单元结点位移列阵/用整体结构结点位移列阵&表示，即

=Ga（2.2.26）

其中

v,u2v2••-vr.…//„vjy

12・・・2/-12/・・•2m-12m・・・2y-l2j…2n

6x2/;

10・-00

01---00

00---00

00---00

00---10

00---01

00…0

00---0

10---0

01---0

00---0

00---0

于是（2.2.26）〜（2.2.28）式代入

（2.2.25b）式,

则离散形式的总位能可表示为

（2.2.27）

（2.2.29）

Yl^a1-^（GrKeG）a-a1^（G'Pc）

2QQ

令

K=^GTKeG

P=/GV（2.2.30）

K和P分别称为整体结构刚度矩阵和结构结点载荷列阵，于是总位能可表示为

np=-alKa-a1P（2.2.31）

卩2

山于离散形式的总位能的未知量是结构的结点位移根据变分原理，泛函口‘,取

驻值的条件是它的一次变分为零，况Ip=0,即

dnn

化=—^8a=0卩oa

由于加的任意性，则欲使上式成立，只有

dnn

—=0（2.2.32）

加

上式表明：

总位能对任一自由度分量求导都为零，若结点数为n,总的自由度数

为2n（平面问题），即得2n个方程。

其矩阵形式为

Ka

（2.2.33）

Ill（2.230）式可以看出结构整体刚度矩阵K和结构结点载荷列阵P都是山单元

刚度矩阵Ke和单元等效结点载荷列阵/集合而成。

2.2.3单元刚度矩阵1.单元刚度矩阵

III（2.2.28）式定义的单元刚度矩阵，山于应变矩阵B对于3结点三角形单

元是常量矩阵，因此有

K爲

K爲

K爲

K爲

K爲

Ke=BTDBtA=

（2.2.34）

将弹性矩阵D和应变矩阵B代入上式，任意分块矩阵可表示成

K“=BrDBstA=细]—，"

..1—v叽+—crc,

；1一"[

VCrbs+—hrCs

J1一"z

1—"f2c工、+

2

（2.2.35）

显然

（Kj=K”

（2.2.37）

2・单元刚度矩阵的力学意义和性质

为了进一步理解单元刚度矩阵的物理意义，同样可以利用最小位能原理建立一个单元的有限元方程

Keae=Pe+Fe（2.2.38）

是单元等效结点载荷（体力和面力引起的等效结点力），Fe是其他单元对该

单元的作用力,

P和Fe统称为单元结点力。

叫

・・・

+

—

V

P八

・・・

・・・

・・・

・・•

kJ

pmy

（2.2.40）

这是单元结点平衡方程，每个结点在x和y方向上各有一个平衡方程，3个结点

共有6个平衡方程。

令//,.=!

而其余结点位移皆等于零，由上式可得

e

Pi、

心心

=

+

F亦

（2.2.41）

上式表明：

单元刚度矩阵第一列的元素的物理意义是当i结点在x方向有一单位

位移色=1,其他结点位移都为零时，在单元各结点位移方向上施加结点力的大小。

当然，单元在这些结点力作用下应处于平衡，因此在x和y方向上结点力之和应为零，即

在x方向Kh+K3i+K,]=0

在y方向++=0（2.2.42）

单元刚度矩阵的性质：

（1）对称性

这一性质在迦辽金形式的加权余量法和里兹法已加以说明，在3结点三角形单元的刚度矩阵的具体表达式（2.2.37）已加以证明。

实际上对于各种形式的单元都普遍具有这种对称性。

（2）奇异性

III（2.2.42）可以见到，由于单元刚度矩阵的每一列都代表结构的平衡，各列元素之和等于零。

因此单元刚度矩阵是奇异的。

因此有限元方程组的解不唯一，而是包含了一系列的刚体位移，为消除结构的刚体位移，结构应有足够的约束。

（3）主元恒大于零

单元刚度矩阵对角线上的元素称为主元，由（2.235）及（2236）式可见他们是恒正的。

K"的意义是要使i结点在x方向的位移©=1,在x方向所施加的结点力。

显然恒为正。

例2.1给定3结点三角形如图2.3所示。

边长“=2,b=4,c=5,厚度/=1,材料

常数E=2x\05MPa,v=0.2o计算单元的刚度矩阵，并验证矩阵的对称性，奇

异性和主元恒正。

解：

按（2.2.7）式计算，得到

q=

=12,

=

0,G

3=0

2

~4上2=

4,b

3=0

C]=

_3,

6=

■

0,c

3=3

'b

0

b

0b

0_

1

1

j

m

B=

2A

0

Cj

0

Jo

J

Si

b.

Cj

bj5

h.n_

19.6

7.2

-16

-4.8

—3.6

-2.4

7.2

15.4

-2.4

-6.4

-4.8

-9

Ke

=R1DRtA=

8681

-16

-2.4

0

0

0

2.4

-4.8

—640

6.4

+4.8

+0

-3.6

-4.80

4.8

3.6

0

-2.4

-9

2.4

0

0

9

2.2.4单元等效结点载荷列阵

111（2.2.28）式得到单元等效结掉载荷是pc=P；十P,

体积力的等效结点载荷:

P；=^NTftdxdy

其中

P；=

p;

7/

tdxdy

Pe

L讥

[Jy]

L—

-

「笊

r

A

P"

p-_

=\M

Jnr1f

fy_

tdxdy

O加）

面积力的等效结点载荷:

P；=\N『TtdS

P鳥f[T1

P；=lx=fMtdS

L^iLrvJ

常见的载荷如下：

1.均质等候单元的自重

单元的单位体积重量为恣，坐标方向如图2.4所示。

按照（2.2.28）式，应有

P；=

tdxdv

✓

（2.2.45）

P—

=i

p

0

0'

■0'

-pg-

tdxdy

一

0

-

_o_

nNipgtdxdy一

~~P8^

o,jm

[N打

其中，每个结点的等效结点载荷为

（2.2.46）

自重的等效结点载荷是

Pep=_1^m[01010l]7（2.2.47）

2•均布压力

设在ij边上有一压力q,如图2.5所示。

ij边长为1,与x轴的夹角为Q。

压力

q在x和y方向的分量么,亦分别为

qx=qs\na

qy=-qcosa=j（x.-x,）作用在单元边界上的面力为

（2.2.48）

在边界上取局部坐标s,山面积坐标可知，沿ij边插值函数可以写作

1—CC

Ni=一,N，Nm=0（2.2.49）

1/7/

将（2.2.28）代入（2.2.48）及（2249）式代入（2.2.28）,得

&=\NqJds=J（1-\qjds=Lq（y.-y.）

兀=A（勺-九）

P；x=[Njqjds=qxtds=如（y,--儿）

Piy=尹（勺-“）

X一儿Xj—Xj

0o]7

（2.2.50）

3・x方向均布力（如图2.6）

町010.0Of（2.2.51）

2

4.x方向三角形分布载荷

载荷作用在ij边，如图2.7所示。

这时边界面积力写成局部坐标s的函数,即

（i-yk

0

则单元等效结点载荷为

（2.2.52）

2.2.5结构刚度矩阵和结构结点载荷列阵的集成

（2.2.30）式给出了结构刚度矩阵和结构结点载荷列阵山单元刚度矩阵和单

元等效结点载荷列阵集成的表达式。

1.单元刚度矩阵的转换

刚度矩阵的转换表示为

1•…i・・・j・・・m・・・n

iTo•…o•…o•…o・・・o

io…K；…K：

j…K：

m…0

（2.2.53）

GTKeG=-=j0…K；…K；…Kl.n…0

0…K爲

mi

其中n为结构结点总数为单元编号。

在上式的矩阵中除标明的9个子块

外，其它位置上的元素均为零。

i,j,m均为双行和双列。

单元矩阵的这种变换起到两个作用：

（1）将单元刚度矩阵扩大到与整个结构刚度矩阵同阶，便于矩阵相加；

（2）将单元刚度矩阵中各子矩阵按照单元结点的实际编码安放在扩大的矩阵中。

即对号入座。

2.单元等效结点载荷列阵的转换

单元等效结点载荷列阵的转换表示为

（2.2.54）

i,j,m均为双行，单元等效结点载荷包括体积力和面积力等的等效结点载荷。

山

（2.2.54）可见单元等效结点载荷列阵的转换是将单元结点载荷扩大到与结构载荷列阵同价，并将单元结点载荷按结点自由度顺序入位。

即对号入座。

3.结构刚度矩阵和结构载荷列阵的集成

举例：

4.结构刚度矩阵的性质和特点

（1）对称性

（2）奇异性；

（3）稀疏性

（4）非零元素在矩阵对角线附近呈带状分布。

半带宽：

3=2（D+1）,D为相邻结点的最大差值2.2.6引入位移边界条件

对于求解位移场问题，山于最后所得的有限元方程的系数矩阵为奇异阵，所表示的物理含义：

满足方程的位移包含刚体位移。

为了求解一个具体结构的受力行为，应有足够的约束，以消除结构的刚体位移，从而消除刚度矩阵的奇异性。

1.直接代入法

在方程组（2233）中将已知结点位移的自由度消去，得到一组修正方程，用以

解其他待定的结点位移。

其原理是按结点位移已知和代定重新组合方程

K』

'pt:

（2.2.55）

其中乙为待定结点位移，乙为已知结点位移；山刚度矩阵的对称性可知

Kha=K：

b。

由（2.2.55）式展开，得

Kaaaa+Kahah=Pa（2256）

曲于勺已知，最后的求解方程可写成

KZ*=P（2.2.57）

其中

卍=K“,a=aa,P'=P-Kahah（2.2.58）

若总体结点位移为n个，其中已知结点位移m个，则得到一组求解n-m个待定结点位移的修正方程组，K”为n-m阶方阵。

修正方程组的意义是在原来方程中,只保留与未知结点位移相应的n-m个方程，并将方程中左端的已知位移和相应刚度系数的乘积移至方程右端作为载荷的修正项。

这种方法：

组成的新的方程阶数降低了，但结点位移的顺序已被破坏，这给编程带来一定的麻烦。

2.对角元素改1法（主元置1法）

当给定位移是零位移时，可以在系数矩阵K中与零结点位移相对应的行和列中，将主对角线元素改为1,其它元素改为0,并在载荷列阵中将与零结点位移相对应的元素改为0。

例如：

勺=0

心

■J

心

a,

厶

■

P1

0

0

-1

0

=

0

（2.2.59）

心

K”2

…0:

如

3.对角元素乘大数法（主元赋大值）

当有结点位移为给定位移勺=勾，对笫j个方程的方程的对角线元素乘以大

值1O10,并将什•用^10,0代替，则

氐…

g…

■A■

心

■■

••

■•

心T

心

a,

■

Pl

••

心•-

/C.JO10…

心lop

0

©•-

K叫

%.

_P"-

（2.2.60）

则笫j个方程变为

（伯+Kj2a2+…+K力•10⑴勺+5=/C,10,（）a.

由于^10,（）»Kr项比其它项要大得多，于是

KAO^a.^心汕勾

总是能保证

对于多个给定位移时则按其所在位置进行修正。

2.2.7线性代数方程组的求解及应力计算

有限元求解方程在引入位移边界条件，消除了K的奇异性后，就可以求解了。

具体解法将在笫六章专门讲解。

2.3广义坐标有限元的一般格式

2.3.1广义坐标有限元位移模式的选择和插值函数的构造

常用的二维单元有三角形或矩形，常用的三维单元有四面体、五面体或平行六面体。

同样形状的单元还可有不同的单元结点数，如二维三角形3结点单元、三角形6结点单元、三角形10结点单元，因此单元种类较多。

图2.11中列举了一些常用的单元形式。

如何选择合适的单元进行计算，涉及求解问题的类型、对精度的要求以及经济性多方面的因素。

当然，不同的单元位移模式也不相同。

本节讨论选择广义坐标有限元位移模式的一般原则，以及建立位移插值函数的一般步骤。

1.选择广义坐标有限元位移模式的一般原则

（1）广义坐标是山结点场变量确定的，因此它的个数应与结点自山度数相等。

如3结点三角形单元有六个自由度，广义坐标（0,）应取6个，因此两个方向的位移各取三项多项式。

对4结点的矩形单元广义坐标数为8,位移函数可取四项多项式作为近似函数。

（2）选择多项式时，位移模式必须完备，即包含可能有的变形。

对于C。

类问题,就是位移模式应有常数项和坐标的一次项。

位移模式中的常数项和一次项反映了单元刚体位移和常应变特性。

当单元数趋于无穷时，单元缩小趋于一点，此时单元应变应趋于常应变。

（3）多项式的选取应山低阶到高阶，尽量选取完全多项式以提高单元精度。

（4）多项式的选取不应偏惠于某一坐标，即x或y应等同出现。

一般可根据巴斯卡三角形（或杨辉三角形）选取。

具体见表2.1

2.建立广义坐标有限元位移插值函数的一般步骤

i）假设位移模式；

ii）将各结点坐标代入，得到关于广义坐标的线性方程组，从而求出广义坐标；

iii）将广义坐标几回代入一般位移模式中得到山单元结点位