有限元 2-弹性力学平面问题有限单元法(2.1三角形单元,2.2几个问题的讨论)

第2章 弹性力学平面问题有限单元法

2.1 三角形单元(triangular Element)

三角形单元是有限元分析中的常见单元形式之一，它的优点是:

①对边界形状的适应性较好，②单刚形式及其推导比较简单,故首先介绍之。 一、结点位移和结点力列阵

设右图为从某一结构中取出的一典型三角形单元。

在平面应力问题中，单元的每个结点上有沿x 、y 两个方向的力和位移，单元的结点位移列阵规定为： 相应结点力列阵为: (式2-1-1)

二、单元位移函数和形状函数

前已述及，有限单元法是一种近似方法，在单元分析中，首先要求假定（构

造）一组在单元内有定义的位移函数作为近似计算的基础。即以结点位移为已知量，假定一个能表示单元内部（包括边界）任意点位移变化规律的函数。 构造位移函数的方法是：以结点(i,j,m)为定点。以位移(u i ,v i ,…u m v m )为定点上的函数值，利用普通的函数插值法构造出一个单元位移函数。 在平面应力问题中，有u,v 两个方向的位移,若假定单元位移函数是线性的,则可表示成:

(,)123

u u x y x y ααα==++

546(,)v v x y x y ααα==++ (2-1-2)a

式中的6个待定常数α1 ,…, α6 可由已知的6个结点位移分量(3个结点的坐标)

{}???

??

?????=?????

????

?????????????=m j i m e

d d d d m j j i v u v u v u i {}

i

i j j m X Y X (2-1-1)Y X Y i

e

j m m F F F F ??

??

????

????

??==??????????????????

确定。将3个结点坐标(x i,y i ),(x j,y j ),(x m,y m )代入上式得如下两组线性方程:

1

23i i i u x y ααα

=++

123

j j j u x y ααα=++ (a)

1

23m m m u x y ααα

=++

和

54

6i i i v x y ααα

=++

546

j j j v x y ααα=++ (b)

54

6m m m v x y ααα

=++

利用线性代数中解方程组的克来姆法则,由(a)可解出待定常数1α 、2α 、

3α :

11A A

α=

22

A A

α

=

33

A A

α

=

式中行列式:

1

i i i j j j m m m u x y A u x y u x y =

2111i i j j m m

u y A u y u y =

3111i i j j

m m

x u A x u x u =

2111i i j j m m

A

x y A x y x y ==

A 为△ijm 的面积,只要A 不为0,则可由上式解出:

11()2m m i i

j j a u a u a u A α=++ 21()2m m i i

j j bu b u b u A α=++ （C ）

31()2m m

i i j j c u c u c u A α=

++

式中:

m m i j j a x y x y =- m m j i i a x y x y =- m i j j i a x y x y =-

m i j b y y =- m j i b y y =- m i j b y y =- （d ）

m i j c x x =- m j i c x x =- m j i c x x =-

为了书写方便，可将上式记为: m m i j i a x y x y =- m i

j b

y y =- (,,)i j m m i j

c x x =-

(,,)i j m 表示按顺序调换下标,即代表采用i,j,m 作轮换的方式便可得到(d)

式。 将(c)式代入2-1-2中,整理后可得:

(,)(,)(,)m m i i j j

u N x y u N x y u N x y u =++

同理: (,)(,)(,)m m i i j j

v N x y v N x y v N x y v =++ （2-1-2）b 式中: 1()2i

i i i N a b x c y A

=++ (,,)i j m (2-1-3)

将三角形单元的位移函数用矩阵表示：

或：

4)a -1-(2

v u i v u v u m 0 j 0 i 0 0 m 0 j 0 ),(m j i i ?????

?

?????????????????????

???=??????????=m N N N N N N v u f i y x 4)b

-1-(2 }]{[}{e d N f v u =???

???????=

三、单元的应变和应力

１、应变──几何矩阵[B]

由弹性力学知,弹性力学平面问题的几何方程： x

u x

ε?=?； y v y ε?=? ； xy u v y x ε??=+??

用矩阵表示

00x y xy x u v y y x εεε?

?

??????????????????????????

????????????

???????

?

???=?????

或， {}{}H f ε????

= (2-1-5)

[H]称为微分符矩阵,又称为微分算子,“[H]{f}”实际上不是一般的矩阵乘,

可以称为微分符矩阵[H]作用在{f}上,其作用规律符合矩阵乘积规定，实际上是按[H]对{f}求导。

将2-1-4式的{f}=[N]{d}代入: {ε

}=[H][N]{d}=[B]{d} 2-1-6

式中: (2-1-7)

0000

00000000012m i j m i j m i j m i j m m i

i j j x

N N N B H N y N N N y x b b b c c c c b c b c b A ????????????

????????????????

?

???

?????????

?

?

????????

???

??

???==?????=

称为几何矩阵，对于上述三角形单元，[B]是常量矩阵,因此常把这种三角形单

元称为常应变单元

２、应力矩阵[S]

由弹性力学知,由应力求应变的物理方程是:

()1x x y E

σμσε=-

()1y y x E

σμσε=-

2(1)xy xy xy

E

G

τμγτ+==

由上式解出应力,得到由应变表示应力的物理方程: ()21x y x E εμεσμ

=+-

()21y x y E εμεσμ

=+-

2121xy

xy E μγτμ-

=?- 用矩阵表示:

{}

{}

x y xy D σσσετ????

????

??????????

??==

21010(1)00

21D E μ

μμμ

?

?

??????????

???????

?

=-- 2-1-8

称为弹性矩阵。

将2-1-6式{ε}=[B]{d}代入上式得: {}{}{}

{}

e

e

D D B d S d σε????????????????=== 2-1-9 式中:

221111*********m m i i j j m m i i j j m m i i j j S D B b

c b c b c E b c b c b c A c b c b c b μμμμμμμμ

μμμμμ???????

????????????

???

?

?

?

???

???

?

????????

?

==-------

2-1-10

上式称为应力矩阵,是结点位移与应力之间的关系矩阵,在上述三角形单元中它也是一个常量矩阵(常应力单元)。

四、单元刚度矩阵

有了几何矩阵[B]和弹性矩阵[D]后,我们便可将其代入在§1-3中推导出的单元刚度矩阵的一般表达式:

T

v

K B D B dv ??????????

??????

=?

对于平面问题，积分

dz t =?，是单元的厚度，并假定t 在单元内不变化(常数)，

所以三角形单元的单刚:

T K t B D B dxdy ????????????????

=??

积分式内的两个矩阵都是常量,矩阵乘后,积分得三角形单元的 单元刚度矩阵：

ii ij im ji jj jm mm mi mj K K K K K K K K K K ????

????????????

= 2-1-11

子块:

21122114(1)22r s r s r s r s rs r s r s r s r s b b c c b c c b Et K A c b b c c c b b μμμμμμμ????

????

??

???

?

--++=

---++ (r,s=i,j,m)

上式常称为各向同性常应变单刚(stiffness matrix for isotropic constant strain triangle in plane stress)

上面的[K]是在图1及2-1-1式的位移排列顺序和(2-1-2)a 下导得的,对此作些改变,便会获得不同的单元刚度矩阵形式,如将位移列阵的排列改为: {d}=[ui uj um vi vj vm ]T 。

（可自行推导,此种单元刚度矩阵的显式可参见《Finite Element Analysis Fundamentals 》）

作业3：

1．求形函数N i (x,y)在三角形形心（x c ,y c ）上的函数值。

2．设图中i 点有水平位移u i =1，试由单元刚度方程写出 各链杆的反力；并证明各水平、竖向反力之和为0。

3．求图示单元1、2的单元刚度矩阵和应力矩阵。

（结论：将单元逆时针转动180度，则单刚无影响而应力矩阵反号）

2.2 三角形单元中几个问题的讨论

一、形函数的物理意义 由(2-1-4)a

{

}

{}

000(,)

000

i i e

m i j j m j i j m m u v N N N u u f x y N d

v v N N N u v ?????????

???????????????????

????

???????

??

???????????

==

可以看出,当u ι =1,其它5个位移分量为0时,u(x,y)=Ni(x,y),或当v i =1其余为0时 v (x,y)=N i (x,y)(如图所示)。

故形函数Ni(x,y)表示当结点i 发生单位位移时,在单元内部产生的位移分布状态,函数Nj ,Nm 亦具有类似的性质,因此,Ni ,Nj ,Nm 称为位移的形状(态)函数,简称形函数,[N]即称形函数矩阵（shape function ）。

二、形函数的几何意义

在单元分析中，很重要的一步是构造位移函数，得出以形函数和结点位移乘积表示的单元位移场。

{}

{}

000000

i i e m i

j j m j i j m m u v N N N u f N d v N N N u v ????

???????????

???=?

?????

?

?

??

???

???????????

= 其中形函数: 12i i i i N a b x c y A ?? ??

?=++ (),,i j m 而

m i j b y y =-

x m i j c x =-

将其代回:

j

m m j i y x y x a -=

1

()()21121m m m m i j i j j j j m m

N x y x y y y x x x y A x y x y A

x y ??

????

=-+-+-= 若从单元内任意点p(x,y)向各顶点引连线,将其分成三个小三角形,则上式中的行列式恰好是小三角形pjm 面积的2倍

即: 1211221i x y A

A i N x y i j j A A A

x y m m

=== 同理可得:

j j A N A = m m A N A

=

由此可得如下结论(三角形单元的几何性质)

1. 任意一点形函数之和等于1,（ Ni +Nj +Nm =1）

2. 形函数为≥0,且 ≤1的值, 0≤(Ni ,Nj ,Nm )≤1

3. 顶点坐标上的形函数值:

当(x,y)坐标取在i 点时, Ni =1, Nj =Nm =0 当(x,y)坐标取在j 点时, Ni =1, Nj =1, Nm =0 当(x,y)坐标取在m 点时

, Ni =Nj =0, Nm =1

因此:Ni ,Nj ,Nm 又称为面积坐标

面积坐标的概念在讲叙六结点三角形单元时将得到应用。

三、有限元的收敛性

解的收敛性也可理解为一个问题的解的精度，较粗的分，影响一个实际问题的解的精度可分成三个方面：

实际物理问题→①理想化力学模型→②有限元求解方法（解法）→③数字截断。

此处仅讨论解法。绪论中提到，有限元作为一种数值方法可以认为是李兹法（弹性力学解）的一种特殊形式，不同之处在于有限元法的形函数（在弹力称试探函数）是定义于单元（子域）而不是全域。李兹法的收敛条件是要求试函数具有完全性和连续性。那么它在有限元法中又是如何具体体现的？ 可从两方面：

严格的数学论证（用变分原理，可参考王勖成书pp60-63）； 物理方面，也就是单元模型问题（三角形单元的位移摸式）。

１、位移模式的收敛性条件 ① 完备性

一个有限元解的必要条件是要求该单元的位移函数必须能表示刚体位移和常应变状态。 例如对于上述平面问题的三角形单元，其位移函数应能保证单元能够在不产生应变的前提下，在其平面内的任意方向上平移和转动，如图示悬臂深梁即可解释上述之理由。

现考察上述三角形单元位移函数对刚体位移的描述。

设单元产生刚体位移，原点的水平位移U 0

；竖向位移V0；转角φ0， 则任

意点位移：

00(,)(,)u x y u y v x y v x

=-=+00φφ

显然它包含在我们所设的位移模式中。

对于常应变问题，可以这样理解：当有限元网络不断细分，每个单元的尺度都

趋于很小的极限尺寸时，每个单元的应变即应趋于常数，结构内部任何复杂变化都可以被近似。由于三角形单元本身就是常应变单元，故它自然满足。

②连续性

又称为协调性要求。是指变形后相邻单元在公共边界处的位移（变形）是连续的，即相邻单元之间既不发生“间隙”又不“重叠”。

从能量原理上讲，如果单元在交界面上的位移不连续，将在交界面引起无限大的应变，而我们在建立能量方程（总势能表达式）或泛函数公式时没有考虑这种情况，因此，有限元解就不可能（难于）收敛于真正解。

对于上述三角形单元，由于位移函数是线性的，即单元中的一条直线变形后仍为直线，而相邻单元在两个公共点上的位移又应是相等的（位移谐调），所以条件②得到满足。从而保证了变形后相邻单元在公共边界处的位移（变形）是连续的

条件２得不到满足时，不能收敛到精确解，但不等于它不收敛，对有些位移函数，尽管它不能满足条件②，但当单元大小划分得恰当时，仍能获得很好的近似结果，这便非协调元问题。

因此常把条件１称为必要条件，条件２称为充分条件。

２、位移法解的下限性质

用能量变分法可以证明，有限元位移法得到的位移解，总体上不大于真正解。即解具有下限性质。

位移解下限性质可以这样理解：所划分单元是原来连续体的一部分，具有无限多个自由度，在假定了单元的位移函数后，自由度限制为只有以结点位移表示的有限自由度，即位移函数对单元的变形进行了约束和限制，使单元的刚度较实际加强了，因此连续体刚度随之加强，故求得的近似解总体上（而不是每一点）将小于精度解。

又由此可知，杆件结构原本只在结点约束，且位移模式与假设者相同（精确方法：力法，位移法），故有限元法对于杆件结构所求得的是精确解。

３、求解误差问题

结构内力分析程序的有效性和解题速度等，在很大程度上都取决于方程组的求解方法。线性方程组的解法很多，限于篇幅这里将不作具体介绍。读者可参阅有关书籍和借鉴已有程序中的方法。

引起求解过程中出现误差的原因很多，也很复杂。除了程序设计本身的错误等问题外，其中还有许多属于数学中的误差理论问题。但不管怎样，一般来说，遇到下列几种情况，方程组的解就有可能出现问题（“病态”）１、总刚度矩阵中元素间的数量级相差很远，且无一定规律；

２、消去过程中出现严重的有效数位损失或数值很小的元素等；

３、矩阵的行列式值相对来说很小，或其中某些行（列）接近线性相关。为了尽量克服方程组“病态”，保证解的正确性，通常采取的办法有：

１、改变算法，寻求一种较好的求解方法：

２、提高数学的有效位，如采用双精度（甚至双双精度）变量等；

３、尽量避免使矩阵元素间数量级相差太大；

４、尽量减小带宽。一般来说，带宽较小的矩阵比带宽较大的矩阵误差会更小；

５、尽量将位移值较大的位移分量排在位移列阵的前面，比如在高层结构计算简图中，结点编号采用从上到下的顺序等。

五、单刚的性质

１、对称性（位移互等定理）

２、奇异性（即，不能根据结点力反求结点位移）

３、分块性质，按结点划分子块

六、荷载等效变换

在有限元分析中，都须把作用在单元内部或边界上的荷载统统按静力等效原则转换到结点上，变成结点荷载代入方程组后,方可求解结点位移（如图）。

静力等效的原则是转换后的结点荷载与原荷载在任意虚位移上所的虚功相等。下面利用虚位移原理推导其一般公式。

荷载变换的一般公式

设在单元内部或边界上的任意点M(x,y,z)处作用有荷载P(如图)。其外荷载列阵用{P}表示,单元的等效结力列阵用{R}表示。

设由于该单元发生某一虚位移场{f}*

而引起的结点虚位移为{d}* ,

}*

荷载作用点M的相应虚位移为{f

M

由静力等效原则有:

}*T{P}

{d}*T{R}={f

M

因虚位移{f}* 应符合位移模式规律:{f}*=[N]{d}*

将其代入得:

{d}*T {R}={d}*T [N M ]T {P} 由此得:

{R}=[N M ]T {P} (2-2-1)

* [N M ]—N 在M 点的形函数值。

如果将{P}看成是作用在单元上的任意广义力(体积力、集中力、弹性力)，则荷载等效变换的一般表达式为:

{R}=∫V [N]T {R}dv (2-2-2)

下面具体分析三角形单元上常见荷载的转换

1、集中力

设在单元ij 边界上的d 点处作用一沿x 方向的集中力P,其作用点到i,j 的距离分别为li ,lj (如图) 由式2-2-1得

0000000000ix i i iy i jx j j jy j m mx m my m R N N R N R N N P P R N N R N R N ??????

??????????????????

??????????????????

?????

?????????

??????????

????????????

??????????

???

???==

根据形函数Ni(i,j,m)的几何性质,它们在d 点上的函数值分别为

j i i l A N A l ==; ; 0m N = 于是得:

{}0000T

j i P R l l l ??????

= (2-2-3)

2、分布力

设在单元ij 边界上的作用有沿x 方向的三角形分布力, 取微段ds, 视其为集中力

{}

0q t sds dp l ??????

????????

=

l

l N i

j =

根据集中力计算公式 (2-2-3)，可得dp 产生的等效结点力:

{}

20000T tqs dR ds s l s l

??????=-

于是积分可得:

{}{}

202300132000000l l s l s tqs tql R dR ds ds l ????????????????????????

????????????

????

????

??????

????

-===??

(2-2-4)

如果是沿x 方向作用均布载, 则有

{}[]T tql

R 0001012

= (20101008更正) 由上知, 作用在ij 边界上沿x 方向的力, 等效变换后只存在i,j 结点x 方向的分量,其余四个分量为0,如果外荷载方向是任意的,则可首先将其分解成x,y 方向，然后再作变换。

3、体积力

设在单元内作用有均匀分布的体积力,单位面积上的体积力分量： {}x y W W t W ????

?

???

?

?

= 由式2-2-1得:

{}

0000i i x y A m m N N W R tdxdy W N N ????

?

???

????????????

??????????

=?? 由于t,Wx ,Wy 在单元内为常量,可提到积分号外, 而 ?Ni(x,y)dxdy=A/3 (i,j,m) 故

{},

,,,,3T

x y x y x y At R W W W W W W ??

?

?

??

= 2-2-4

从上面三种荷载的等效变换结果可以看出，它们均符合刚体力学中力的分解、合成原则，其原因是我们所取三角形单元的位移函数是线性函数。换句话说，荷载等效变换的结果与单元位移函数（形函数）有关。

作业4：

１、求图示三角形单元的等效结点荷载

已知：t=1.5cm, Wx =0 Wy =-0.003kg/cm3

P=6kg q=2kg/cm a=20cm

部分答案:

Rjx =-15.0111kg, Rjy =-8.9265kg。

２、用有限元法计算图示结构的各支座反力。

E=constant, μ=0.15

参考答案:X2 =0.4039→，Y2 =0

（分成两个三角形单元）