下载掌阅APP，畅读海量书库
立即打开
在有限元法中经常会用到弹性力学的基本方程和与之等效的变分原理，现将它们连同相应的矩阵表达形式和张量表达形式引述于后。关于它们的详细推导可从弹性力学的有关教材中查到。
1.弹性力学基本方程的矩阵形式
在载荷作用下，弹性体内任意一点的应力状态可由6个应力分量 σ x ， σ y ， σ z ， τ xy ， τ yz ， τ zx 来表示。其中， σ x ， σ y ， σ z 为正应力； τ xy ， τ yz ， τ zx 为剪应力。应力分量的正负号规定如下：如果某一个面的外法线方向与坐标轴的正方向一致，那么这个面上的应力分量就以沿坐标轴正方向为正，沿坐标轴负方向为负；相反，如果某一个面的外法线方向与坐标轴的负方向一致，那么这个面上的应力分量就以沿坐标轴负方向为正，沿坐标轴正方向为负。
应力分量的矩阵形式称为应力列阵或应力向量，即
在载荷作用下，弹性体还将产生位移和变形，即弹性体位置的移动和形状的改变。
弹性体内任意一点的位移都可由沿直角坐标轴方向的3个位移分量 u ， v ， w 来表示。它的矩阵形式为
u 称作位移列阵或位移向量。
弹性体内任意一点的应变都可以由6个应变分量 ε x ， ε y ， ε z ， γ xy ， γ yz ， γ zx 来表示。其中， ε x ， ε y ， ε z 为正应变； γ xy ， γ yz ， γ zx 为剪应变。应变的正负号与应力的正负号相对应，即正应变以伸长为正，缩短为负；剪应变以两个沿坐标轴正方向的线段组成的角度变小为正，变大为负。
应变的矩阵形式为
ε 称作应变列阵或应变向量。
对于三维问题，弹性力学基本方程可写成如下几种形式。
1）平衡方程
在 V 域内，弹性体内任意一点沿 x 轴、 y 轴、 z 轴的平衡方程为
式中，
为单位体积的体积力沿
x
轴、
y
轴、
z
轴的分量，且有
τ
xy
=
τ
yx
，
τ
yz
=
τ
zy
，
τ
zx
=
τ
xz
。
平衡方程的矩阵形式为
式中， A 是微分算子，即
是体积力向量，
。
2）几何方程—应变-位移关系
在微小位移和微小变形情况下，略去位移导数的高次幂，则应变列阵和位移列阵间的几何关系为
几何方程的矩阵形式为
式中， L 为微分算子，即
3）物理方程—应力-应变关系
弹性力学中应力-应变之间的转换关系也称弹性关系。对于各向同性的线弹性材料，应力通过应变的表达式可用矩阵形式表示为
式中
D 称为弹性矩阵。它完全取决于弹性体材料的弹性模量 E 和泊松比 v 。
也可以采用拉梅（Lamé）常数 G 和 λ 表征弹性体的弹性，它们和 E 、 v 的关系为
G 也称剪切弹性模量。观察式（2-12）可得
物理方程中的弹性矩阵 D 亦可表示为
物理方程的另一种表示形式是
式中， C 是柔度矩阵。 C = D -1 ，它和弹性矩阵是互逆关系。
弹性体的全部边界为
S
，对于一部分边界已知外力
，称为力的边界条件，这部分边界用
S
σ
表示；对于另一部分边界已知位移
，称为几何边界条件或位移边界条件，这部分边界用
S
u
表示。这两部分边界构成弹性体的全部边界，即
4）力的边界条件
在边界上弹性体单位面积的内力为
T
x
，
T
y
，
T
z
，已知在边界
S
σ
上弹性体单位面积上作用的面积力为
，根据平衡应有
设边界外法线的方向余弦为 n x ， n y ， n z ，则边界上弹性体的内力为
式（2-17）的矩阵形式为
式中
5）位移边界条件
在边界上弹性体的位移为
u
，
v
，
w
，已知在边界
S
u
上弹性体的位移为
，根据平衡应有
其矩阵形式为
以上是三维弹性力学问题中的一组基本方程和边界条件。同样，对于平面问题、轴对称问题等也可以得到类似的方程和边界条件。
把弹性力学的基本方程记作一般形式，即
并有 S σ + S u = S ， S 为弹性体全部边界。
各类弹性力学问题有关矩阵符号如表2-1所示。板与壳的基本方程将在本书有关章节中给出。
表2-1 各类弹性力学问题有关矩阵符号
续表
6）弹性体的应变能和余能
单位体积的应变能（应变能密度）为
应变能是正定函数，只有当弹性体内所有点都没有应变时（ ε = 0 ），应变能才为零。
单位体积的余能（余能密度）为
余能也是正定函数。在线性弹性力学中弹性体的应变能等于余能。
2.弹性力学基本方程的张量形式
弹性力学基本方程亦可用笛卡儿张量符号来表示，使用附标求和的约定可以得到十分简练的方程表达形式。
在直角坐标系
x
1
，
x
2
，
x
3
中，应力张量和应变张量都是对称的二阶张量，分别用
σ
和
ε
表示，且有
σ
ij
=
σ
ji
和
ε
ij
=
ε
ji
。其他位移张量、体积力张量、面积力张量等都是一阶张量，用
等表示。下面将分别给出弹性力学基本方程及边界条件的张量形式和张量形式的展开式。
1）平衡方程
式中，下标 j 表示对独立坐标 x j 求偏导数； σ ij ， j 项的下标中 j 出现两次，表示该项在该指标的取值范围（1，2，3）内遍历求和。该重复指标称为哑指标。式（2-25）的展开式为
与式（2-4）比较可见，当 x 1 ， x 2 ， x 3 是笛卡儿坐标系时有
2）几何方程
其展开式为
与式（2-6）、式（2-7）比较可见，当 x 1 ， x 2 ， x 3 是笛卡儿坐标系时有
3）物理方程
广义胡克定律假定每个应力分量与各个应变分量成比例。广义胡克定律可以用张量符号表示为
式中，81个比例常数 D ijkl 称为弹性常数，是四阶张量 D 的分量。由于应力张量是对称张量，因此 D ijkl 的前2个指标具有对称性，由于应变张量也是对称张量，因此 D ijkl 的后2个指标也具有对称性，即
当变形过程是绝热或等温过程时，还有
考虑了上述对称性后，对于最一般的线弹性材料，即在不同方向具有不同弹性性质的材料，81个弹性常数中有21个是独立的。对于各向同性的线弹性材料，独立的弹性常数只有2个，即拉梅常数 G 和 λ 或弹性模量 E 和泊松比 ν ，此时弹性张量可以简化为
此时广义胡克定律可以表示为
式中
式（2-34）的展开式为
上面两式中拉梅常数 G 和 λ 与弹性模量 E 和泊松比 ν 的关系如式（2-12）所示。
物理方程的另一种形式为
式中， C ijkl 是柔度张量 C 的分量，它和刚度张量 D 的分量 D ijkl 是互逆关系。
3.力的边界条件
式中
式中， n j 是边界外法线 n 的三个方向余弦。
将式（2-39）代入式（2-38）后，它的展开式为
4.位移边界条件
6.应变能和余能
单位体积应变能，即
单位体积余能，即
3节点三角形单元是有限元法中最早提出的并且至今仍广泛应用的单元，由于三角形单元对复杂边界有较强的适应能力，因此很容易将一个二维域离散成有限个三角形单元。在边界上以若干段直线近似原来的曲线边界，随着单元增多，这种拟合将趋于精确。本书在讨论如何应用有限元法分析各类具体问题时，将以平面问题3节点三角形单元为例来阐明弹性力学平面问题有限元分析的表达格式和一般步骤。
1.单元位移模式及插值函数的构造
典型的3节点三角形单元节点编码为（ i ， j ， m ），以逆时针方向编码为正向，每个节点都有2个位移分量，即
每个单元都有6个节点位移，也就是6个节点自由度，即
1）单元的位移模式和广义坐标
在有限元法中单元的位移模式（有时称为位移函数）一般采用多项式作为近似函数，因为多项式运算简便，并且随着项数的增多可以逼近任何一段光滑的函数曲线。多项式的选取应由低次到高次。
3节点三角形单元位移模式选取一次多项式，即
它的矩阵形式为
式中
φ 称为位移模式，它表示位移作为坐标（ x ， y ）的函数中包含的项次。单元内的位移是坐标（ x ， y ）的线性函数； β 1 ～ β 6 是待定系数，称为广义坐标。6个广义坐标可由单元的6个节点位移来表示。在式（2-46）中代入节点 i 的坐标（ x i ， y i ）可得到节点 i 在 x 轴方向的位移 u i ，同理可得 u j 和 u m 。它们表示为
解式（2-51）可以得到广义坐标由节点位移表示的表达式。式（2-51）的系数行列式是
式中， A 是三角形单元的面积。
广义坐标 β 1 ～ β 3 为
同理，利用3个节点在 y 轴方向的位移，即式（2-47），可求得
在式（2-53）和式（2-54）中
式中，（ i ， j ， m ）表示下标轮换，如 i → j ， j → m ， m → i 。以下同此。
2）位移插值函数
将求得的广义坐标 β 1 ～ β 6 代入式（2-46）和式（2-47），可将位移函数表示成节点位移的函数，即
式中
N i ， N j ， N m 称为单元的插值函数或形函数，对于当前情况，它是坐标（ x ， y ）的一次函数。 a j ， b j ， c j ， a m ， b m 是常数，取决于3节点三角形单元的坐标。
式（2-58）中的三角形单元的面积 A 可用式（2-55）的系数表示为
式（2-56）和式（2-57）的矩阵形式为
式中， N 称为插值函数矩阵或形函数矩阵； a e 称为单元节点位移列阵。
插值函数具有如下性质。
在节点上插值函数的值有
即有 N i （ x i ， y i ）=1， N i （ x j ， y j ）= N i （ x m ， y m ）=0。也就是说在节点 i 上 N i =1，在节点 j 、节点 m 上 N i =0。由式（2-56）和式（2-57）可见，当 x = x i ， y = y i ，即在节点 i 上，应有 u = u i ，因此必然要求 N i =1， N j = N m =0。其他两个形函数也具有同样的性质。此性质称为Kronecker Delta性质。
在单元中任意一点各插值函数之和应等于1，即
因为若单元发生刚体位移，如 z 轴方向有刚体位移 u 0 ，则单元内（包括节点上）到处应有位移 u 0 ，即 u i = u j = u m = u 0 ，又由式（2-56）得到
因此，必然要求 N i + N j + N m =1。若插值函数不满足此要求，则不能反映单元的刚体位移，用来求解必然得不到正确结果。单元的各个节点位移插值函数之和等于1的性质称为归一性。
对于现在的单元，插值函数是线性的，单元内部及单元边界上的位移也是线性的，可由节点上的位移值唯一地确定。相邻单元公共节点的位移是相等的，这保证了相邻单元在公共边界上位移的连续性。
3）应变矩阵和应力矩阵
确定了单元位移后，可以很方便地利用几何方程和物理方程求得单元的应变和应力。
式（2-22）的几何方程中的位移用式（2-60）代入，得到单元应变为
式中， B 为应变矩阵； L 为平面问题的微分算子（见表2-1）。
应变矩阵 B 的分块子矩阵是
对式（2-58）求导可得
将式（2-66）代入式（2-65）可得
3节点三角形单元的应变矩阵是
式中， b i ， b j ， b m ， c i ， c j ， c m 由式（2-55）确定，它们是单元形状参数。当单元的节点坐标确定后，这些参数都是常量［与坐标变量（ x ， y ）无关］，因此 B 是常量阵。当单元的节点位移 a e 确定后，由 B 转换求得的单元应变都是常量，也就是说在载荷作用下单元中各点的 ε x 、 ε y 及 γ xy 相同。因此，3节点三角形单元称为常应变单元。在应变梯度较大（就是应力梯度较大）的部位，单元划分应适当密集，否则不能反映应变的真实变化面，从而导致较大误差。
单元应力可以根据物理方程求得，即在式（2-22）的物理方程中代入式（2-64）可得
式中， S 为应力矩阵。将平面应力或平面应变的弹性矩阵（见表2.1）及式（2-68）代入式（2-70），可以得到计算平面应力或平面应变问题的单元应力矩阵。 S 的分块矩阵为
式中， E 0 和 v 0 为材料常数。
对于平面应力问题有
对于平面应变问题有
与应变矩阵 B 相同，应力矩阵 S 也是常量阵，即3节点三角形单元中各点的应力是相同的。在很多情况下，不单独定义应力矩阵 S ，而直接用 DB 进行应力计算。
2.利用最小位能原理建立有限元方程
最小位能原理的泛函总位能 Π p 的表达式在平面问题中的矩阵表达形式为
式中， t 为二维体厚度； f 为作用在二维体内的体积力； T 为作用在二维体边界上的面积力。
对于离散模型，系统位能是各单元位能的和，将式（2-60）及式（2-64）代入式（2-74），即可得到离散模型的总位能为
将结构总位能的各项矩阵表达成各个单元总位能各项矩阵之和，要求单元各项矩阵的阶数（单元的节点自由度数）和结构各项矩阵的阶数（结构的节点自由度数）相同。为此需要引入单元节点自由度和结构节点自由度的转换矩阵 G ，从而将单元节点位移列阵 a e 用结构节点位移列阵 a 表示，即
式中
式中， n 为结构的节点数。
令
式中， K e 和 P e 分别为单元刚度矩阵和单元等效节点载荷列阵。将式（2-76）～式（2-78）代入式（2-75），则可得离散形式的泛函总位能为
并令
式中， K 和 P 分别为结构整体刚度矩阵和结构节点载荷列阵。这样一来，式（2-79）就可以表示为
由于离散形式的泛函总位能 Π p 的未知变量是结构的节点位移 a ，根据变分原理，泛函总位能 Π p 取值的条件是它的一次变分为零， δΠ p =0，即
这样就得到有限元的求解方程为
式中， K 和 P 由式（2-80）求得。由式（2-80）可以看出，结构整体刚度矩阵 K 和结构节点载荷列阵 P 分别是由单元刚度矩阵 K e 和单元等效节点载荷列阵 P e 集合而成的。
以上表述的是基于弹性力学最小位能原理形成有限元求解方程的一般原理。在具体计算时涉及单元刚度矩阵的形成、单元等效节点载荷列阵的形成，以及分别集合单元刚度矩阵和单元等效节点载荷列阵形成的结构刚度矩阵和结构等效节点载荷列阵，还有给定位移边界条件的引入等问题。这些将在后面进行讨论。
1.单元刚度矩阵的形成
由于对于3节点三角形单元，应变矩阵 B 是常量阵，因此根据式（2-78）定义的单元刚度矩阵为
将弹性矩阵 D 和应变矩阵 B 代入式（2-84）后，它的任一分块矩阵可表示为
式中
由式（2-85）可以得到
由此可知，单元刚度矩阵是对称矩阵。
2.单元刚度矩阵的力学意义和性质
为了进一步理解单元刚度矩阵的力学意义，可以利用最小位能原理建立一个单元的求解方程，从而得到
式中， P e 为单元等效节点载荷； F e 为其他相邻单元对该单元的作用力； P e 和 F e 统称节点力； a e ， P e ， F e 依次表示为
式（2-88）的展开式为
这是单元节点平衡方程，每个节点在 x 轴方向和 y 轴方向上各有一个平衡方程，3个节点共有6个平衡方程。方程左端是用单元节点位移表示的单元节点内力，方程右端是单元节点力（外载荷和相邻单元的作用力之和）。
令
由式（2-92）可得
在力学分析中常用的求解方法有两种：隐式（Implicit）和显式（Explicit）。
隐式有限元与显式有限元最大的区别在于是否迭代，所有物理量是否在同一时刻获得。通过迭代求解平衡方程（位移、速度和加速度），不管是用隐式方法，还是用显式方法（前向或后向欧拉求解方法）求解本构方程（应力和应变），都叫作隐式有限元；用显式方法—时间积分求解本构方程叫作显式有限元。
显式求解和隐式求解的递推公式一个用显函数表示，一个用隐函数表示。
隐式方法适合处理时间跨度大、需要集合整体刚度矩阵的问题，或者对线性方程组求解，但存在收敛性的问题，在时间选择方面需要考虑步长，步长太大误差偏大，步长太小浪费计算资源，耗时过长。
虽然隐式方法在某些情况下在计算效率方面具有优势，但在动力学分析中，需要取很短的时间步长来捕捉瞬时响应。对此显式动力学方法更加合适，由于该方法在时间上分别对单个单元进行连续分析，因此时间步骤可以取得非常小，容易编写程序进行并行分布式计算。因此，在动力学分析中主要采用显式方法，特别是对于响应时间很短的问题。