有限元方法求解弹性力学问题的基本原理

北京理工大学 | 李明健

有限元方法（Finite element method，FEM）是工程中应用最广泛的一种数值方法，起源于 Courant 的工作，正式开始于 1960 年 Clough 的工作。其基本原理是加权余量法和变分法，可用于固体、流体、传热、电磁等多个领域。

弹性力学中的有限元方法是从变分原理（通常是从最小势能原理）出发而建立的。

弹性力学基本方程的矩阵形式：

弹性体中，应力有 6 个独立分量，写成列阵形式：

\boldsymbol \sigma = [\sigma_x , \sigma_y, \sigma_z, \tau_{xy} ,\tau_{yz}, \tau_{xz}]^\text T

应变同样有 6 个分量，写为：

\boldsymbol \varepsilon = [\varepsilon_x , \varepsilon_y, \varepsilon_z, \gamma_{xy} ,\gamma_{yz}, \gamma_{xz}]^\text T

位移有 3 个方向的分量：

\boldsymbol u = [u , v, w]^\text T

外力（体积力）为：

\boldsymbol f = [f_x,f_y, f_z]^\text T

由于静力学平衡方程为：

\left\{ \begin{align} &\frac{\partial \sigma_x}{\partial x}+\frac{\partial \tau_{yx}}{\partial y}+\frac{\partial \tau_{zx}}{\partial z} + f_x=0\\ &\frac{\partial \tau_{xy}}{\partial x}+\frac{\partial \sigma_{y}}{\partial y}+\frac{\partial \tau_{zy}}{\partial z} + f_y=0\\ &\frac{\partial \tau_{xz}}{\partial x}+\frac{\partial \tau_{yz}}{\partial y}+\frac{\partial \sigma_{z}}{\partial z} + f_z=0 \end{align} \right.

写成矩阵形式为：

\boldsymbol A \boldsymbol \sigma + \boldsymbol f =0

其中：

\boldsymbol {A} = \begin{bmatrix} \dfrac{\partial}{\partial x} & 0 & 0 & \dfrac{\partial}{\partial y} & 0 & \dfrac{\partial}{\partial z} \\ 0 & \dfrac{\partial}{\partial y} & 0 & \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial z} & 0 \\ 0 & 0 & \dfrac{\partial}{\partial z} & 0 & \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial x} \end{bmatrix}

几何方程为：

\left\{ \begin{align} &\varepsilon_x = \frac{\partial u}{\partial x}\\ &\varepsilon_y = \frac{\partial v}{\partial y}\\ &\varepsilon_z = \frac{\partial w}{\partial z}\\ &\gamma_{xy} = \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\\ &\gamma_{yz} = \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y}\\ &\gamma_{zx} = \frac{\partial w}{\partial x} + \frac{\partial u}{\partial z} \end{align} \right.

写成矩阵形式为：

\boldsymbol \varepsilon = \boldsymbol L \boldsymbol u

其中：

\boldsymbol {L}=\left[\begin{array} {ccc} \dfrac{\partial}{\partial x} & 0 & 0 \\ 0 & \dfrac{\partial}{\partial y} & 0 \\ 0 & 0 & \dfrac{\partial}{\partial z} \\ \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial x} & 0 \\ 0 & \dfrac{\partial}{\partial z} & \dfrac{\partial}{\partial y} \\ \dfrac{\partial}{\partial z} & 0 & \dfrac{\partial}{\partial x} \end{array}\right]=\boldsymbol{A}^{\text{T}}

线弹性问题的本构关系（物理方程）为：

\left\{ \begin{align} &\sigma_x = \frac{E}{(1+\mu)(1-2\mu)} \left[ (1- \mu) \varepsilon_x + \mu ( \varepsilon_y+\varepsilon_z)\right] \\ &\sigma_y = \frac{E}{(1+\mu)(1-2\mu)} \left[ (1- \mu) \varepsilon_y + \mu ( \varepsilon_x+\varepsilon_z)\right] \\ &\sigma_z = \frac{E}{(1+\mu)(1-2\mu)} \left[ (1- \mu) \varepsilon_z + \mu ( \varepsilon_x+\varepsilon_y)\right] \\ &\tau_{xy} = \frac{E}{2(1+\mu)} \gamma_{xy} \\ &\tau_{yz} = \frac{E}{2(1+\mu)} \gamma_{yz} \\ &\tau_{zx} = \frac{E}{2(1+\mu)} \gamma_{zx} \end{align} \right.

写成矩阵形式为：

\boldsymbol \sigma = \boldsymbol D \boldsymbol \varepsilon

其中弹性矩阵为：

\boldsymbol D = \begin{bmatrix} \lambda+2G & \lambda & \lambda & 0 & 0 & 0 \\ \lambda & \lambda+2G & \lambda & 0 & 0 & 0 \\ \lambda & \lambda & \lambda+2G & 0 & 0 & 0 \\ 0 & 0 & 0 & G & 0 & 0\\ 0 & 0 & 0 & 0 & G & 0\\ 0 & 0 & 0 & 0 & 0 & G\\ \end{bmatrix}

$\lambda$ $G$ （又称为 剪切模量）为：

\lambda=\frac{E\mu}{(1+\mu)(1-2\mu)},\quad G=\frac{E}{2(1+\mu)}

也可以反过来写成：

\boldsymbol \varepsilon =\boldsymbol C \boldsymbol \sigma

$\boldsymbol C=\boldsymbol D^{-1}$ 为柔度矩阵。

位移边界条件为：

\boldsymbol u = \overline{ \boldsymbol u},\quad \text{on} \; S_1

力边界条件为：

\begin{array}{l} p_{x}=n_{x} \sigma_{x}+n_{y} \tau_{y x}+n_{z} \tau_{z x} \\ p_{y}=n_{x} \tau_{x y}+n_{y} \sigma_{y}+n_{z} \tau_{z y} \\ p_{z}=n_{x} \tau_{z x}+n_{y} \tau_{y z}+n_{z} \sigma_{z} \end{array}

矩阵形式：

\boldsymbol n \boldsymbol \sigma = \overline{ \boldsymbol p} ,\quad \text{on} \; S_2

其中边界面上的外法线为：

\boldsymbol{n}=\left[\begin{array}{cccccc} n_{x} & 0 & 0 & {n}_{y} & 0 & {n}_z \\ 0 & n_{y} & 0 & {n}_x & {n}_z & 0 \\ 0 & 0 & n_z & 0 & {n}_y & {n}_x \end{array}\right]

现在的问题，是要根据已知的力和位移边界条件，求解出弹性体内的位移场、应力场、应变场等。

有限元基本原理：

由于复杂情况下的连续场很难直接从数学上直接求解，因此将结构划分为有限个单元（一维问题划分线段，二维为三角形、四边形，三维为四面体、六面体等），然后在每一个单元上采用线性（或高阶）近似，这样，就可以用这些有限个单元的结果来逼近连续的无限自由度问题，即为有限元方法。

连续体划分为单元后，每个单元有若干个广义节点位移自由度，我们可以用节点上的位移来插值单元内任意一点的位移：

\boldsymbol u =\boldsymbol N \boldsymbol u_e = \sum_i^{n_p} N_i \boldsymbol u_i

$\boldsymbol u$ $\boldsymbol u_e$ 的问题。

$(u_i,v_i,w_i)$ ，则单元节点位移为：

\boldsymbol u_e = [u_1,v_1,w_1,u_2,v_2,w_2,...,u_8,v_8,w_8]^\text T

我们用这 8 个节点（共 24 个自由度）插值任意一点，其中形函数为：

N_i = \frac{1}{8}(1+\xi_i \xi)(1+\eta_i \eta )(1+ \zeta_i \zeta)

形函数需要满足“本点为一，它点为零”、以及任一点处“和为一”的性质。

根据前述的几何方程，可以得到单元内的应变为：

\boldsymbol \varepsilon = \boldsymbol L \boldsymbol u = \boldsymbol L \boldsymbol N \boldsymbol u_e = \boldsymbol B \boldsymbol u_e

$\boldsymbol B=\boldsymbol L \boldsymbol N$ 称为应变矩阵。

单元应力为：

\boldsymbol \sigma = \boldsymbol D \boldsymbol \varepsilon = \boldsymbol D \boldsymbol B \boldsymbol u_e

单元应变能密度为：

U_0 = \frac{1}{2} \boldsymbol \sigma \boldsymbol \varepsilon = \frac{1}{2} \boldsymbol \varepsilon^\text T \boldsymbol D \boldsymbol \varepsilon

单元应变能为：

U_e=\frac{1}{2} \int_\Omega \boldsymbol \varepsilon^\text T \boldsymbol D \boldsymbol \varepsilon \text d \Omega

单元的势能为变形能减外力功：

\begin{aligned} \Pi_e &= U_e-W_e \\ &= \frac{1}{2} \int_\Omega \boldsymbol \varepsilon^\text T \boldsymbol D \boldsymbol \varepsilon \text d \Omega - \int_\Omega \boldsymbol u^\text T \boldsymbol f \text d \Omega - \int_S \boldsymbol u^\text T \boldsymbol p \text d S \\ &= \frac{1}{2} \boldsymbol u_e^\text T \int_\Omega \boldsymbol B^\text T \boldsymbol D \boldsymbol B \text d \Omega \boldsymbol u_e - \boldsymbol u_e^\text T \int_\Omega \boldsymbol N^\text T \boldsymbol f \text d \Omega - \boldsymbol u_e^\text T \int_S \boldsymbol N^\text T \boldsymbol p \text d S \end{aligned}

其中，我们定义 单元刚度矩阵 为：

\boldsymbol k_e = \int_\Omega \boldsymbol B^\text T \boldsymbol D \boldsymbol B \text d \Omega

定义 单元等效节点载荷列阵 为：

\boldsymbol P_e = \int_\Omega \boldsymbol N^\text T \boldsymbol f \text d \Omega + \int_S \boldsymbol N^\text T \boldsymbol p \text d S

系统的总势能，是各个单元的势能之和：

\Pi = \sum \Pi_e = \boldsymbol u ^\text T \frac{1}{2} \sum \boldsymbol G^\text T \boldsymbol k^e \boldsymbol G \; \boldsymbol u - \boldsymbol u ^\text T \sum \boldsymbol G^\text T \boldsymbol P_e

$\boldsymbol G$ 是单元自由度编号对应总体编号的转换矩阵。

定义 总体刚度矩阵：

\boldsymbol K = \sum \boldsymbol G^\text T \boldsymbol k^e \boldsymbol G

定义 结构节点载荷列阵：

\boldsymbol P =\sum \boldsymbol G^\text T \boldsymbol P_e

则总势能就可以表示为：

\Pi = \frac{1}{2} \boldsymbol u ^\text T \boldsymbol K \boldsymbol u - \boldsymbol u ^\text T \boldsymbol P

根据最小势能原理，在所有可能的位移场中，真实位移场使系统总势能最小，即对势能的一阶变分为零：

\frac{\delta \Pi}{\delta \boldsymbol u} =0

即得到静力学有限元的求解方程：

\boldsymbol K \boldsymbol u= \boldsymbol P

动力学有限元需要考虑惯性力和阻尼力，写为：

\boldsymbol M \ddot{ \boldsymbol {u} } + \boldsymbol C \dot{\boldsymbol u} + \boldsymbol K \boldsymbol u= \boldsymbol P

以上。