2.5 弹性力学有限元求解

力场基本方程为

其中，M为质量矩阵；C为阻尼矩阵；K为刚度矩阵；u为节点位移矩阵。

静态问题：

2.5.1 基于最小势能原理的变分法求解

2.5.1.1 变分法推导控制方程

以简支梁为例，利用变分法推导控制方程。两端简支梁受一横向载荷（见图2-9），其弹性势能是由梁弯曲变形引起的应变能，外力势能是由梁的扰度引起横向载荷做功产生势能。设梁的扰度为v（x），外力势能为Π₁，弹性势能为Π₂，则

弹性势能：

外力势能：

总势能为Π=Π₁+Π₂，即

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

将式（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处的扰度为v_i，转角为θ_i，则每个单元的位移插值函数可表示为

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

式中

将式（2-36）代入式（2-25）有

由式（2-38）可知，求得每个单元的δΠ_e，代入式（2-38）就可以求得每个节点的v_i和θ_i（i=1，2，n，n+1），从而得到位移函数v_e（x）。

1. 求单元位移函数

图2-11所示为任一单元（e）变形后的位移曲线v_e（x），其在i，j节点处所示的位移与转角均为正方向。现设单元位移函数v_e（x）为

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

代入v_e（x）= a₀+a₁x + a₂x²+ a₃x³和v′_e（x）= a₁+2a₂x +3a₃x²，得

解得：

代入式（2-39），按v_i，θ_i，v_j，θ_j的顺序加以整理得

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

式中

令

N =（N₁N₂N₃N₄），δ^e=（v_i θ_i v_j θ_j）^T，得

2. 形状函数

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

3. 求Π_e（单元总势能）

由式（2-43）得

将式（2-43）写成如下形式：

将式（2-45）、式（2-46）代入式（2-37），有

4. 求δΠ_e

因为N_i是x的已知函数，所以δΠ是由节点位移的变分而引起的。故式（2-47）的变分为

将提到积分号外

对式（2-49）第1项进行积分

将代入式（2-50）并积分得

对式（2-47）第2项进行积分

将N_i代入式（2-52）并积分得

不考虑梁单元轴向位移时，式（2-51）是梁单元的刚度矩阵K^e，而式（2-54）中大括号{}内的第1项，是梁单元由节点位移v_i、θ_i、v_i、θ_j引起的节点力V_i、M_i、V_j、M_j。式（2-53）是承受均布载荷q的两端固定梁的固端反力，由上向下依次为V_0i、M_0i、V_0j、M_0j，将式（2-54）写成：

由δΠ=0求解节点位移，将每个单元的δΠ_e式（2-55）代入式（2-38），得

由式（2-54）可知，与δ^e中节点位移的排列顺序是一样的，若将式（2-56）各项按梁的节点顺序排列，注意到δδ=（δv₁δθ₁……δv_n+1δθ_n+1）的任意性，则由式（2-56）可得

式（2-57）中KZ是由每个单元刚度矩阵K^e集合而成（整体刚度矩阵），FZ₀是由每个单元的固端反力集合而成，若将该矩阵前加上“-”号，则它就是等效节点载荷。

2.5.2 二维问题的有限元求解

1. 三角形单元的有限元求解

（1）位移函数 图2-13为任一三角单元，其三个节点的局部码为1、2、3，以逆时针为序。其节点坐标为（x₁，y₁）、（x₂，y₂）、（x₃，y₃），其节点位移为（u₁，v₁）、（u₂，v₂）、（u₃，v₃）。该单元内任一点（x，y）处的位移函数为

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

若节点位移和节点坐标值均已知，则由式（2-59）可解出α₁、α₂、α₃，由式（2-60）可解出α₄、α₅、α₆，即

式中

合并相同位移量整理式（2-58）得

令

则

令

则式（2-65）变为

式（2-66）是以节点位移表示的三角单元的位移函数，其中N₁、N₂、N₃是形状函数。其特点是N_i在本节点的值是1，在其他节点的值为零，图2-14为N₁为1的情况。

（2）单元刚度矩阵 由节点位移求应变，将式（2-63）代入式（2-5）有

写成矩阵形式为

令

则式（2-68）变为

B称为应变矩阵，该应变为常应变，即在单元内各点应变均为一个常数，这是由于所设位移函数是线性函数的缘故。

由节点位移求应力，将式（2-70）代入式（2-10）有

令

则式（2-71）变为

S称为应力矩阵，因为D、B均为常数，所以在单元内σ为常数，故称三角形为常应力单元。

由节点位移求节点力—单元刚度矩阵。

图2-15a为平面内的任一三角单元，设平面受力后节点产生位移u₁、v₁、u₂、v₂、u₃、v₃（其内部各点的位移由位移函数决定），同时产生节点力U₁、V₁、U₂、V₂、U₃、V₃（节点位移与节点力用一箭头表示），而其内力为σ，即σ=DBδ^e。现给该单元一虚位移（见图2-15b），此时3个节点将发生虚位移δu₁、δv₁、δu₂、δv₂、δu₃、δv₃，内部将产生虚应变。

式中

δ δ^e =（δu₁δv₁δu₂δv₂δu₃δv₃）^T

该三角单元的外力虚功为

式中

F^e =（U₁V₁U₂V₂U₃V₃）^T

一般弹性体的内力虚功为

W_int =∫_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是对称的，所以K^e也是对称的。对不同的问题，K^e中的B和D是不同的。一般情况下，B为函数矩阵，需经积分运算。对于平面问题的简单三角形单元，B是常数矩阵。式（2-80）虽然是由平面问题简单三角形推导而得，但具有普遍性，是位移法有限元分析中普通适用的单元刚度阵表达式。

2. 矩形单元有限元求解

（1）位移函数 图2-16是长为2a、宽为2b的矩形单元，局部坐标节点码如图所示。因其有8个节点位移，故其位移函数可设为

式中

上面的N_i（i=1，2，3，4）是根据形状函数的特点（在i节点N_i=1，在其他节点N_i=0）而得到的。

（2）单元刚度矩阵 将式（2-82）代入式（2-5）

式中

δ^e =（u₁v₁u₂v₂u₃v₃u₄v₄）^T

由式（2-72）得

由式（2-80）得

式中

1）如果域内有不同方位的矩形单元，须将它们的单元刚度矩阵K^e都变换到整体坐标方向才能构成总刚度矩阵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项节点位移，以列阵表示为

{δ}^e={u₁ v₁ w₁ u₂ v₂ w₂ u₃ v₃ w₃ u₄ v₄ w₄}^T

{δ}^e称为单元节点位移，这种单元具有12个自由度。单元变形时，单元内各点也有沿x、y、z方向的位移u、v、w，一般应为坐标x、y、z的函数。对于这种简单的四面体单元，其内部位移可假设为坐标的线性函数，即

u=a₁+a₂x+a₃y+a₄z

v=a₅+a₆x+a₇y+a₈z

w=a₉+a₁₀x+a₁₁y+a₁₂z

上式含有12个a参数，可以由单元的12项节点位移确定，将4个节点的坐标值代入上式中的u式，对1、2、3、4这4个节点分别有：

u₁=a₁+a₂x₁+a₃y₁+a₄z₁

u₂=a₁+a₂x₂+a₃y₂+a₄z₂

u₃=a₁+a₂x₃+a₃y₃+a₄z₃

u₄=a₁+a₂x₄+a₃y₄+a₄z₄

式中，V由计算下列行列式得到：

V是四面体的体积，方程式（2-87）中的系数α_i、β_i、γ_i和δ_i（i=1，2，3，4）由下式计算：

在方程式（2-87）中，以v_i或w_i项代替所有的u_i，可以得到v或w的表达式。u的位移表达式和v及w的表达式相似，可以等同地写成以形状函数和未知节点位移表示的展开形式：

式中，形状函数为

2. 单元刚度矩阵

将式（2-88）代入几何关系式（2-6），经过微分运算，可以得到单元内应变为

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

式中，下标后的逗号表示相对后面的变量取微分，把式中的下标i分别替换为1、2、3、4，就可以得到子矩阵B₁、B₂、B₃、B₄。

将方程式（2-89）代入方程式（2-92），经过微分运算，B₁表示为

同理得到B₂、B₃和B₄为

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

这里，V^e为单元体积，由于简单四面体单元为常应变单元，故积分有：

[K]^e= [B]^T[D][B]V^e

3. 不同节点数单元的位移函数

（1）8节点六面体单元8节点六面体单元如图2-18所示，每个节点沿其坐标方向x、y、z共有3个平移自由度。

以节点位移和形状函数表示的单元位移函数为

（2）10节点四面体单元10节点四面体单元如图2-19所示，是在三维线性四面体单元基础上建的一种高阶单元。与4节点的四面体单元相比，10节点四面体单元更适于精度要求较高、边界为曲线时的模型。

单元的位移函数可表示为

u = u_I（2S₁-1）S₁+u_J（2S₂-1）S₂+ u_K（2S₃-1）S₃+u_L（2S₄-1）S₄+

4×（u_MS₁S₂+ u_NS₂S₃+ u_OS₁S₃+ u_PS₁S₄+ u_QS₂S₄+u_RS₃S₄）

v = v_I（2S₁-1）S₁+ v_J（2S₂-1）S₂+ v_K（2S₃-1）S₃+ v_L（2S₄-1）S₄+

4×（v_MS₁S₂+ v_NS₂S₃+ v_OS₁S₃+ v_PS₁S₄+ v_QS₂S₄+ v_RS₃S₄）

w = w_I（2S₁-1）S₁+w_J（2S₂-1）S₂+ w_K（2S₃-1）S₃+w_L（2S₄-1）S₄+

4×（w_MS₁S₂+w_NS₂S₃+ w_OS₁S₃+ w_PS₁S₄+w_QS₂S₄+ w_RS₃S₄）

（3）20节点六面体单元20节点四面体单元如图2-20所示，它是在8节点六面体单元的基础上建立的一种高阶单元。与8节点的六面体单元相比，20节点四面体单元更适于精度要求高、边界为曲线模型。

单元的位移函数可表示为

位移v和w分量的位移函数，与u向位移分量相似。

不同节点的单元刚度矩阵推导与4节点四面体推导类似，这里不再赘述。

4. 载荷分配

三维弹性体内如受有均布的体积力（如重力）作用，对于简单的四面体单元，则可以逐个单元计算出其整个单元的全部体积力，再平均分配到四个结点上，即每个节点分配到1/4的单元体积力。

如果单元的某个表面作用有均布的面积力（如气体压力），也可将此面上的全部面积力平均分配到相连的三个结点上，即每个结点分配到三角面上面积力总和的1/3。

如果体积力、面积力不是均匀的，应按做功相等的原则等效分配，即

其中，{Q_V}^e、{Q_S}^e为e单元内分布体积力和分布面积力分配到单元结点的载荷，[N]为形状函数矩阵，{q}和{p}分别为单位体积力和面积力，V^e、S^e则为受有分布力的单元体积和面积。

2.5.4 弹性力学有限元求解方法

1. 求解基本流程

1）建立研究对象的近似模型；

2）将研究对象分割成有限数量的单元；

3）对每个单元提出一个近似解；

4）将所有单元按标准方法组合成一个与原有系统近似的系统；

5）用数值方法求解这个近似系统；

6）计算结果处理与结果验证。

2. 主要步骤

1）列出以位移为待解未知量的有限元的平衡方程。

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

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