
2.5 弹性力学有限元求解
力场基本方程为

其中,M为质量矩阵;C为阻尼矩阵;K为刚度矩阵;u为节点位移矩阵。
静态问题:

2.5.1 基于最小势能原理的变分法求解
2.5.1.1 变分法推导控制方程
以简支梁为例,利用变分法推导控制方程。两端简支梁受一横向载荷(见图2-9),其弹性势能是由梁弯曲变形引起的应变能,外力势能是由梁的扰度引起横向载荷做功产生势能。设梁的扰度为v(x),外力势能为Π1,弹性势能为Π2,则

图2-9 两端简支梁
弹性势能:

外力势能:

总势能为Π=Π1+Π2,即

由最小势能原理可知,粱在外力作用下处于稳定平衡的条件是其总势能取极小值。

将式(2-31)第1项中的微分变分符号互调,并将该式代入式(2-25)有

将式(2-32)的第1项分步积分两次变为

式(2-33)中前两项给出了边界条件,第3项给出了控制方程。无论是两端固定梁、两端简支梁,还是一端固定,另一端自由的情况,式(2-33)前两项均为0,依据变分法得控制方程:

由式(2-33)可知,若能找到一个位移函数v(x),它既满足式(2-33)的前两项(边界条件),又满足最后一项(控制方程),则它就是问题的精确解。但是,由于选择一个精确位移函数v(x)很不容易,所以在实际应用中往往只让v(x)精确地满足式(2-33)中的部分项,而近似地满足另外的项。
2.5.1.2 有限元法
将简支梁分为n个单元,其单元号与节点号如图2-10a所示。若任一节点i处的扰度为vi,转角为θi,则每个单元的位移插值函数可表示为

图2-10 受弯简支梁的插值函数挠度曲线

各单元位移插值函数所连成的曲线作为梁的位移函数,如图2-10b所示。总势能沿梁长L的积分变成沿每个单元长度le积分之和,即

式中

将式(2-36)代入式(2-25)有

由式(2-38)可知,求得每个单元的δΠe,代入式(2-38)就可以求得每个节点的vi和θi(i=1,2,n,n+1),从而得到位移函数ve(x)。
1. 求单元位移函数
图2-11所示为任一单元(e)变形后的位移曲线ve(x),其在i,j节点处所示的位移与转角均为正方向。现设单元位移函数ve(x)为

图2-11 梁单元位移函数

式中,ai(i=0,1,2,3)为待定系数。若其两端的位移与转角为已知,则由边界条件可求出ai。边界条件为

代入ve(x)= a0+a1x + a2x2+ a3x3和v′e(x)= a1+2a2x +3a3x2,得
解得:
代入式(2-39),按vi,θi,vj,θj的顺序加以整理得

上式就是插值形式的位移函数,写成矩阵形式为

式中

令
N =(N1N2N3N4),δe=(vi θi vj θj)T,得

2. 形状函数
式(2-42)中的Ni(i=1,2,3,4)称为形状函数。在梁单元中,它表示一个两端固定梁只产生一个单位位移时的梁弯曲形状,如图2-12所示。

图2-12 梁单元的形状函数
3. 求Πe(单元总势能)
由式(2-43)得


将式(2-43)写成如下形式:

将式(2-45)、式(2-46)代入式(2-37),有

4. 求δΠe
因为Ni是x的已知函数,所以δΠ是由节点位移的变分而引起的。故式(2-47)的变分为

将提到积分号外

对式(2-49)第1项进行积分

将代入式(2-50)并积分得

对式(2-47)第2项进行积分

将Ni代入式(2-52)并积分得


不考虑梁单元轴向位移时,式(2-51)是梁单元的刚度矩阵Ke,而式(2-54)中大括号{}内的第1项,是梁单元由节点位移vi、θi、vi、θj引起的节点力Vi、Mi、Vj、Mj。式(2-53)是承受均布载荷q的两端固定梁的固端反力,由上向下依次为V0i、M0i、V0j、M0j,将式(2-54)写成:

由δΠ=0求解节点位移,将每个单元的δΠe式(2-55)代入式(2-38),得

由式(2-54)可知,与δe中节点位移的排列顺序是一样的,若将式(2-56)各项按梁的节点顺序排列,注意到δδ=(δv1δθ1……δvn+1δθn+1)的任意性,则由式(2-56)可得

式(2-57)中KZ是由每个单元刚度矩阵Ke集合而成(整体刚度矩阵),FZ0是由每个单元的固端反力集合而成,若将该矩阵前加上“-”号,则它就是等效节点载荷。
2.5.2 二维问题的有限元求解
1. 三角形单元的有限元求解
(1)位移函数 图2-13为任一三角单元,其三个节点的局部码为1、2、3,以逆时针为序。其节点坐标为(x1,y1)、(x2,y2)、(x3,y3),其节点位移为(u1,v1)、(u2,v2)、(u3,v3)。该单元内任一点(x,y)处的位移函数为

图2-13 三角形单元

将节点坐标值和位移值代入式(2-58)有

若节点位移和节点坐标值均已知,则由式(2-59)可解出α1、α2、α3,由式(2-60)可解出α4、α5、α6,即

式中


合并相同位移量整理式(2-58)得

令

则

令

则式(2-65)变为

式(2-66)是以节点位移表示的三角单元的位移函数,其中N1、N2、N3是形状函数。其特点是Ni在本节点的值是1,在其他节点的值为零,图2-14为N1为1的情况。

图2-14 N1描述的形状
(2)单元刚度矩阵 由节点位移求应变,将式(2-63)代入式(2-5)有

写成矩阵形式为

令

则式(2-68)变为

B称为应变矩阵,该应变为常应变,即在单元内各点应变均为一个常数,这是由于所设位移函数是线性函数的缘故。
由节点位移求应力,将式(2-70)代入式(2-10)有

令

则式(2-71)变为

S称为应力矩阵,因为D、B均为常数,所以在单元内σ为常数,故称三角形为常应力单元。
由节点位移求节点力—单元刚度矩阵。
图2-15a为平面内的任一三角单元,设平面受力后节点产生位移u1、v1、u2、v2、u3、v3(其内部各点的位移由位移函数决定),同时产生节点力U1、V1、U2、V2、U3、V3(节点位移与节点力用一箭头表示),而其内力为σ,即σ=DBδe。现给该单元一虚位移(见图2-15b),此时3个节点将发生虚位移δu1、δv1、δu2、δv2、δu3、δv3,内部将产生虚应变。

图2-15 三角单元的节点力、内力与它们对应的虚位移与虚应变

式中
δ δe =(δu1δv1δu2δv2δu3δv3)T
该三角单元的外力虚功为

式中
Fe =(U1V1U2V2U3V3)T
一般弹性体的内力虚功为
Wint =∫vδ εT σdv
式中,δεT为对应虚位移的虚应变,σ为原平衡力系引起的应力。
将式(2-71)、式(2-74)代入上式,则三角单元的内力虚功为

根据虚位移原理,式(2-75)和式(2-76)相等。

若单元厚度为h、面积为A,再将、δe提到积分外,式(2-77)变为

由于的任意性,故有下式成立:

式(2-79)为一矢量等式,是6个代数方程,表示出e单元的6个节点力与6个节点位移分量之间的相互关系。
令

则

式(2-81)是节点力矩阵表达式,式(2-80)是三角单元刚度矩阵,因为是在整体坐标下推导的,故无需再进行坐标变换。由于弹性系数矩阵D是对称的,所以Ke也是对称的。对不同的问题,Ke中的B和D是不同的。一般情况下,B为函数矩阵,需经积分运算。对于平面问题的简单三角形单元,B是常数矩阵。式(2-80)虽然是由平面问题简单三角形推导而得,但具有普遍性,是位移法有限元分析中普通适用的单元刚度阵表达式。
2. 矩形单元有限元求解
(1)位移函数 图2-16是长为2a、宽为2b的矩形单元,局部坐标节点码如图所示。因其有8个节点位移,故其位移函数可设为


图2-16 矩形单元
式中

上面的Ni(i=1,2,3,4)是根据形状函数的特点(在i节点Ni=1,在其他节点Ni=0)而得到的。
(2)单元刚度矩阵 将式(2-82)代入式(2-5)

式中

δe =(u1v1u2v2u3v3u4v4)T
由式(2-72)得

由式(2-80)得

式中

1)如果域内有不同方位的矩形单元,须将它们的单元刚度矩阵Ke都变换到整体坐标方向才能构成总刚度矩阵KZ。坐标变换矩阵由4个矢量变换矩阵构成,它是一个8×8阶的对角矩阵。
2)至于节点载荷列阵的计算,可利用虚功等效原则去计算。
2.5.3 三维问题的有限元求解
1. 四面体三维单元形状函数
实际物体是立体的,弹性体受力作用后,其内部各点将沿x、y、z三个坐标方向发生位移,以u、v、w表示各点沿x、y、z方向的位移:
u = u(x、y、z)
v = v(x、y、z)
w = w(x、y、z)
三维结构可以划分成很多小四面体,大量的小四面体拼合起来,可以逼近任意形状的实际结构。以四面体顶点为节点,可以构造最简单的三维体单元。
图2-17表示一简单四面体单元,其中4个节点的编号设定为1、2、3、4。单元变形时,各节点都有x、y、z的3项位移,单元有4个节点,共有12项节点位移,以列阵表示为

图2-17 四面体单元
{δ}e={u1 v1 w1 u2 v2 w2 u3 v3 w3 u4 v4 w4}T
{δ}e称为单元节点位移,这种单元具有12个自由度。单元变形时,单元内各点也有沿x、y、z方向的位移u、v、w,一般应为坐标x、y、z的函数。对于这种简单的四面体单元,其内部位移可假设为坐标的线性函数,即
u=a1+a2x+a3y+a4z
v=a5+a6x+a7y+a8z
w=a9+a10x+a11y+a12z
上式含有12个a参数,可以由单元的12项节点位移确定,将4个节点的坐标值代入上式中的u式,对1、2、3、4这4个节点分别有:
u1=a1+a2x1+a3y1+a4z1
u2=a1+a2x2+a3y2+a4z2
u3=a1+a2x3+a3y3+a4z3
u4=a1+a2x4+a3y4+a4z4

式中,V由计算下列行列式得到:

V是四面体的体积,方程式(2-87)中的系数αi、βi、γi和δi(i=1,2,3,4)由下式计算:

在方程式(2-87)中,以vi或wi项代替所有的ui,可以得到v或w的表达式。u的位移表达式和v及w的表达式相似,可以等同地写成以形状函数和未知节点位移表示的展开形式:

式中,形状函数为

2. 单元刚度矩阵
将式(2-88)代入几何关系式(2-6),经过微分运算,可以得到单元内应变为

其中,应变矩阵 [B]是形状函数矩阵经微分算子矩阵作用的结果,[B]中任一个子矩阵 [Bi]为

式中,下标后的逗号表示相对后面的变量取微分,把式中的下标i分别替换为1、2、3、4,就可以得到子矩阵B1、B2、B3、B4。
将方程式(2-89)代入方程式(2-92),经过微分运算,B1表示为

同理得到B2、B3和B4为

[Bi]的每项元素都是由节点坐标决定的常数,因而四面体单元内各点的应变相同,即是常应变单元。由方程式(2-10)可知,单元内应力也是常数。一般受力情况下,构成三维体的有限个大小不同的四面体内的应力并不是常值,用常应力单元来代替它是近似的,单元间的应力是不连续的。只有当单元划分得很小时,单元内的应力才接近于常数。将三向应力中D的表达式及式(2-94)代入式(2-80)即可得到单元刚度矩阵:

这里,Ve为单元体积,由于简单四面体单元为常应变单元,故积分有:
[K]e= [B]T[D][B]Ve
3. 不同节点数单元的位移函数
(1)8节点六面体单元8节点六面体单元如图2-18所示,每个节点沿其坐标方向x、y、z共有3个平移自由度。

图2-18 8节点六面体单元
以节点位移和形状函数表示的单元位移函数为

(2)10节点四面体单元10节点四面体单元如图2-19所示,是在三维线性四面体单元基础上建的一种高阶单元。与4节点的四面体单元相比,10节点四面体单元更适于精度要求较高、边界为曲线时的模型。
单元的位移函数可表示为
u = uI(2S1-1)S1+uJ(2S2-1)S2+ uK(2S3-1)S3+uL(2S4-1)S4+
4×(uMS1S2+ uNS2S3+ uOS1S3+ uPS1S4+ uQS2S4+uRS3S4)
v = vI(2S1-1)S1+ vJ(2S2-1)S2+ vK(2S3-1)S3+ vL(2S4-1)S4+
4×(vMS1S2+ vNS2S3+ vOS1S3+ vPS1S4+ vQS2S4+ vRS3S4)
w = wI(2S1-1)S1+wJ(2S2-1)S2+ wK(2S3-1)S3+wL(2S4-1)S4+
4×(wMS1S2+wNS2S3+ wOS1S3+ wPS1S4+wQS2S4+ wRS3S4)
(3)20节点六面体单元20节点四面体单元如图2-20所示,它是在8节点六面体单元的基础上建立的一种高阶单元。与8节点的六面体单元相比,20节点四面体单元更适于精度要求高、边界为曲线模型。
单元的位移函数可表示为



图2-19 10节点四面体单元

图2-20 20节点六面体单元
位移v和w分量的位移函数,与u向位移分量相似。
不同节点的单元刚度矩阵推导与4节点四面体推导类似,这里不再赘述。
4. 载荷分配
三维弹性体内如受有均布的体积力(如重力)作用,对于简单的四面体单元,则可以逐个单元计算出其整个单元的全部体积力,再平均分配到四个结点上,即每个节点分配到1/4的单元体积力。
如果单元的某个表面作用有均布的面积力(如气体压力),也可将此面上的全部面积力平均分配到相连的三个结点上,即每个结点分配到三角面上面积力总和的1/3。
如果体积力、面积力不是均匀的,应按做功相等的原则等效分配,即

其中,{QV}e、{QS}e为e单元内分布体积力和分布面积力分配到单元结点的载荷,[N]为形状函数矩阵,{q}和{p}分别为单位体积力和面积力,Ve、Se则为受有分布力的单元体积和面积。
2.5.4 弹性力学有限元求解方法
1. 求解基本流程
1)建立研究对象的近似模型;
2)将研究对象分割成有限数量的单元;
3)对每个单元提出一个近似解;
4)将所有单元按标准方法组合成一个与原有系统近似的系统;
5)用数值方法求解这个近似系统;
6)计算结果处理与结果验证。
2. 主要步骤
1)列出以位移为待解未知量的有限元的平衡方程。

2)运用变分最小势能等有限元法求位移,由平衡方程得

3)将本构方程代入求解。
