虚功原理¶

基于应力与守恒定律中导出的平衡方程和边界条件，推导作为连续体力学边值问题弱形式的虚功原理。有限元法的离散化以该弱形式为出发点。本章给出当前构形形式（用 Cauchy 应力和 Almansi 应变的线性部分表示）以及参考构形形式（用第二 Piola-Kirchhoff 应力和 Green-Lagrange 应变表示），在说明二者等价后，再确认其向小变形情况的归结。

平衡方程与边界条件¶

设作用于连续体的物体力（单位质量）为 \(\boldsymbol{g}\)，考虑在当前构形中占据区域 \(\Omega\) 的物体。边界 \(\Gamma\) 分为位移规定为 \(\bar{\boldsymbol{u}}\) 的几何边界 \(\Gamma_B\)，以及表面力规定为 \(\bar{\boldsymbol{t}}\) 的力学边界 \(\Gamma_t\)，并满足 \(\Gamma = \Gamma_B \cup \Gamma_t\)、\(\Gamma_B \cap \Gamma_t = \emptyset\)。对于静力问题，将应力与守恒定律所示动量守恒定律中的惯性项去掉，得到

\[ \nabla_x \cdot \boldsymbol{\sigma} + \rho \boldsymbol{g} = \boldsymbol{0} \quad \text{在} \ \Omega \]

作为平衡方程。边界条件为

\[ \boldsymbol{\sigma} \boldsymbol{n} = \bar{\boldsymbol{t}} \quad \text{在} \ \Gamma_t \]

\[ \boldsymbol{u} = \bar{\boldsymbol{u}} \quad \text{在} \ \Gamma_B \]

由上式给出。以下将虚功原理导出为平衡方程和力学边界条件 \(\boldsymbol{\sigma} \boldsymbol{n} = \bar{\boldsymbol{t}}\) 的弱形式。几何边界条件 \(\boldsymbol{u} = \bar{\boldsymbol{u}}\) 通过测试函数的选择纳入。

当前构形形式的弱形式¶

在弱形式中，分别定义未知位移的容许空间和测试函数空间为

\[ \mathcal{U} = \{ \boldsymbol{u} \in [H^1(\Omega)]^d \mid \boldsymbol{u} = \bar{\boldsymbol{u}} \ \text{在} \ \Gamma_B \} \]

\[ \mathcal{V} = \{ \delta \boldsymbol{u} \in [H^1(\Omega)]^d \mid \delta \boldsymbol{u} = \boldsymbol{0} \ \text{在} \ \Gamma_B \} \]

其中 \(d\) 为空间维数，\(H^1(\Omega)\) 是函数本身及其一阶弱导数均平方可积的 Sobolev 空间，\(\delta\) 为变分符号。在当前构形表示中，\(\Omega\) 是变形后的构形；在实际数值求解中，会将其拉回到参考构形或已知的中间构形进行处理。

将平衡方程乘以权函数 \(\delta \boldsymbol{u} \in \mathcal{V}\)，并应用 Gauss 散度定理和力学边界条件，可得当前构形中的虚功原理如下。

\[ \int_{\Omega} \boldsymbol{\sigma} : \delta \boldsymbol{A}_{(L)}\, dv = \int_{\Gamma_t} \delta \boldsymbol{u}^T \bar{\boldsymbol{t}}\, d\Gamma + \int_{\Omega} \delta \boldsymbol{u}^T \rho \boldsymbol{g}\, dv \]

这里，\(\boldsymbol{A}_{(L)}\) 是 Almansi 应变张量的线性部分，定义为

\[ \boldsymbol{A}_{(L)} = \frac{1}{2}\left( \nabla_x \boldsymbol{u} + (\nabla_x \boldsymbol{u})^T \right), \qquad A_{(L)ij} = \frac{1}{2}\left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \]

其变分为 \(\delta \boldsymbol{A}_{(L)} = \tfrac{1}{2}(\nabla_x \delta \boldsymbol{u} + (\nabla_x \delta \boldsymbol{u})^T)\)。也就是说，求 \(\boldsymbol{u} \in \mathcal{U}\)，使其对任意 \(\delta \boldsymbol{u} \in \mathcal{V}\) 都满足虚功方程。左端为内力虚功，右端为规定表面力和物体力产生的外力虚功。

由于该式写在变形后的（当前构形）区域上，实际求解时会重新选择初始构形 \(\Omega_0\)（参考构形）或已知的中间构形作为参考构形，改写为增量形式后再求解。关于参考构形的具体选择（Total Lagrange / Updated Lagrange）以及增量分解，参见增量分析框架。

初始构形形式的弱形式¶

考虑在参考构形中占据区域 \(\Omega_0\) 的物体，并将其边界 \(\Gamma_0\) 分为 \(\Gamma_{0B} \cup \Gamma_{0t}\)。将当前构形表示拉回到参考构形后，应力与应变的共轭对变为第二 Piola-Kirchhoff 应力 \(\boldsymbol{S}\) 与 Green-Lagrange 应变 \(\boldsymbol{E}\)。此时，初始构形中的虚功原理为

\[ \int_{\Omega_0} \boldsymbol{S} : \delta \boldsymbol{E}\, dV = \int_{\Gamma_{0t}} \delta \boldsymbol{u}^T \bar{\boldsymbol{t}}\, d\Gamma_0 + \int_{\Omega_0} \delta \boldsymbol{u}^T \rho_0 \boldsymbol{g}\, dV \]

其中 \(\rho_0\) 是参考构形中的质量密度，根据质量守恒定律 \(\rho_0 = J\rho\)，它与当前构形中的物体力表示等价。

当前构形形式与初始构形形式的等价性¶

两种表示中的内力虚功通过变形梯度 \(\boldsymbol{F}\) 和体积比 \(J = \det \boldsymbol{F}\) 的变换相互一致，即

\[ \int_{\Omega_0} \boldsymbol{S} : \delta \boldsymbol{E}\, dV = \int_{\Omega} \boldsymbol{\sigma} : \delta \boldsymbol{A}_{(L)}\, dv \]

外力项通过质量守恒定律和表面力变换也彼此等价。因此，当前构形的虚功方程与初始构形的虚功方程，是同一原理在不同构形中的表示。以参考构形为参照的求解方法对应 Total Lagrange 法，以当前构形（前一收敛构形）为参照的求解方法对应 Updated Lagrange 法。

向小变形的归结¶

在小变形假设 \(\boldsymbol{F} \approx \boldsymbol{I}\)、\(J \approx 1\) 下，当前构形与参考构形的区别消失，第二 PK 应力与 Cauchy 应力一致（\(\boldsymbol{S} \to \boldsymbol{\sigma}\)），Green-Lagrange 应变和 Almansi 应变的线性部分都归结为小应变 \(\boldsymbol{\varepsilon}\)。

\[ \boldsymbol{\varepsilon} = \nabla_S \boldsymbol{u} = \frac{1}{2}\left( \nabla \boldsymbol{u} + (\nabla \boldsymbol{u})^T \right), \qquad \varepsilon_{ij} = \frac{1}{2}\left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \]

此时，虚功原理归结为由 Cauchy 应力 \(\boldsymbol{\sigma}\) 和小应变 \(\boldsymbol{\varepsilon}\) 表示的弱形式

\[ \int_{\Omega} \boldsymbol{\sigma} : \delta \boldsymbol{\varepsilon}\, dV = \int_{\Gamma_t} \delta \boldsymbol{u}^T \bar{\boldsymbol{t}}\, d\Gamma + \int_{\Omega} \delta \boldsymbol{u}^T \rho \boldsymbol{g}\, dV \]

\[ \delta \boldsymbol{u} = \boldsymbol{0} \quad \text{在} \ \Gamma_B \]

这就是小变形、线性弹性静力分析中直接用于离散化的弱形式（线性弹性静力分析（入门与附录）以该形式为出发点，说明从单元刚度 \(\boldsymbol{K}^e\) 的构造到总体方程 \(\boldsymbol{K}\boldsymbol{U} = \boldsymbol{F}\) 的组装）。

代入线性弹性本构关系 \(\boldsymbol{\sigma} = \boldsymbol{\mathsf{C}} : \boldsymbol{\varepsilon}\)，并用 Voigt 记法写成 \(\hat{\sigma} = D\, \hat{\varepsilon}\) 后，弱形式变为

\[ \int_{\Omega} \delta \hat{\varepsilon}^T D\, \hat{\varepsilon}\, dV = \int_{\Gamma_t} \delta \boldsymbol{u}^T \bar{\boldsymbol{t}}\, d\Gamma + \int_{\Omega} \delta \boldsymbol{u}^T \rho \boldsymbol{g}\, dV \]

即上述形式。