第3章：拉格朗日视角（下）

本章深入探讨拉格朗日视角下的高级数值方法，重点介绍有限元方法（FEM）的理论基础与实现技巧。我们将从变分原理出发，推导弱形式方程，并探讨不同网格类型的实现策略。通过本章学习，你将掌握构建稳健、高效的固体力学求解器所需的核心技术，并了解如何将这些方法应用于拓扑优化等前沿领域。

3.1 有限元方法的弱形式推导

有限元方法是连续介质力学数值求解的基石。与有限差分方法直接离散微分方程不同，FEM从能量原理出发，通过变分方法得到弱形式方程。这种方法的优势在于自然满足力的平衡，降低了对解的光滑性要求，并且能够灵活处理复杂几何和边界条件。

3.1.1 从强形式到弱形式

考虑弹性体的平衡方程（强形式）： $$-\nabla \cdot \boldsymbol{\sigma} = \mathbf{f} \quad \text{在 } \Omega \text{ 内}$$ 其中 $\boldsymbol{\sigma}$ 是Cauchy应力张量，$\mathbf{f}$ 是体力密度（如重力）。这个方程表达了每个物质点的力平衡：应力的散度加上体力等于零（静力学情况）。

张量记号与指标表示

为了更清晰地理解推导过程，我们同时给出张量和指标表示。平衡方程的分量形式为： $$-\frac{\partial \sigma_{ij}}{\partial x_j} = f_i, \quad i,j = 1,2,3$$ 这里采用Einstein求和约定：重复指标表示求和。应力张量的对称性 $\sigma_{ij} = \sigma_{ji}$ 源于角动量守恒。

边界 $\partial\Omega = \Gamma_u \cup \Gamma_t$ 分为两部分，满足 $\Gamma_u \cap \Gamma_t = \emptyset$：

• 位移边界条件（Dirichlet）：$\mathbf{u} = \overline{\mathbf{u}}$ 在 $\Gamma_u$ 上
• 力边界条件（Neumann）：$\boldsymbol{\sigma} \cdot \mathbf{n} = \overline{\mathbf{t}}$ 在 $\Gamma_t$ 上

其中 $\mathbf{n}$ 是外法向量，$\overline{\mathbf{t}}$ 是施加的表面力（牵引力）。力边界条件的分量形式： $$\sigma_{ij} n_j = \overline{t}_i$$

函数空间的选择

弱形式推导的核心思想：将强形式方程乘以一个任意的测试函数（也称虚位移）$\mathbf{v} \in \mathcal{V}$，其中函数空间： $$\mathcal{V} = \{\mathbf{v} : \mathbf{v} \in H^1(\Omega), \mathbf{v} = \mathbf{0} \text{ 在 } \Gamma_u \text{ 上}\}$$ 这里 $H^1(\Omega)$ 是Sobolev空间，定义为： $$H^1(\Omega) = \left\{\mathbf{v} : \int_\Omega |\mathbf{v}|^2 \, d\Omega < \infty, \int_\Omega |\nabla \mathbf{v}|^2 \, d\Omega < \infty \right\}$$ 选择这个空间的物理意义：

• 函数本身平方可积：保证位移场具有有限的动能
• 导数平方可积：保证应变能有限
• 在 $\Gamma_u$ 上为零：满足位移边界条件的齐次化

分部积分与弱形式

将平衡方程乘以测试函数并在域内积分： $$-\int_\Omega \mathbf{v} \cdot (\nabla \cdot \boldsymbol{\sigma}) \, d\Omega = \int_\Omega \mathbf{v} \cdot \mathbf{f} \, d\Omega$$ 左端的指标形式： $$-\int_\Omega v_i \frac{\partial \sigma_{ij}}{\partial x_j} \, d\Omega$$ 应用分部积分公式（格林公式的张量形式）： $$\int_\Omega \mathbf{v} \cdot (\nabla \cdot \boldsymbol{\sigma}) \, d\Omega = -\int_\Omega \nabla \mathbf{v} : \boldsymbol{\sigma} \, d\Omega + \int_{\partial\Omega} \mathbf{v} \cdot (\boldsymbol{\sigma} \cdot \mathbf{n}) \, d\Gamma$$ 推导细节（以分量形式）： $$\int_\Omega v_i \frac{\partial \sigma_{ij}}{\partial x_j} \, d\Omega = -\int_\Omega \frac{\partial v_i}{\partial x_j} \sigma_{ij} \, d\Omega + \int_{\partial\Omega} v_i \sigma_{ij} n_j \, d\Gamma$$ 这里使用了双点积（双重收缩）运算： $$\nabla \mathbf{v} : \boldsymbol{\sigma} = \frac{\partial v_i}{\partial x_j} \sigma_{ij} = \text{tr}(\nabla \mathbf{v}^T \boldsymbol{\sigma})$$

边界条件的处理

由于边界分为两部分 $\partial\Omega = \Gamma_u \cup \Gamma_t$，边界积分可以分解： $$\int_{\partial\Omega} \mathbf{v} \cdot (\boldsymbol{\sigma} \cdot \mathbf{n}) \, d\Gamma = \int_{\Gamma_u} \mathbf{v} \cdot (\boldsymbol{\sigma} \cdot \mathbf{n}) \, d\Gamma + \int_{\Gamma_t} \mathbf{v} \cdot (\boldsymbol{\sigma} \cdot \mathbf{n}) \, d\Gamma$$ 关键观察：

1. 在 $\Gamma_u$ 上：$\mathbf{v} = \mathbf{0}$（测试函数的定义），因此第一项为零
2. 在 $\Gamma_t$ 上：$\boldsymbol{\sigma} \cdot \mathbf{n} = \overline{\mathbf{t}}$（自然边界条件）

代入得到弱形式的最终表达： $$\int_\Omega \nabla \mathbf{v} : \boldsymbol{\sigma} \, d\Omega = \int_\Omega \mathbf{v} \cdot \mathbf{f} \, d\Omega + \int_{\Gamma_t} \mathbf{v} \cdot \overline{\mathbf{t}} \, d\Gamma$$

虚功原理的物理解释

这个方程可以解释为虚功原理：

• 左端：内力在虚应变 $\nabla \mathbf{v}$ 上做的虚功（内部虚功）
• 右端第一项：体力在虚位移 $\mathbf{v}$ 上做的虚功
• 右端第二项：表面力在虚位移 $\mathbf{v}$ 上做的虚功

虚功原理表明：对于平衡状态，内部虚功等于外部虚功。

弱形式的优势

1. 降低光滑性要求： - 强形式需要 $\boldsymbol{\sigma} \in C^1$（应力场一阶连续可微） - 弱形式只需要 $\mathbf{u} \in H^1$（位移场平方可积且导数平方可积）

2. 自然边界条件： - Neumann边界条件自动满足，无需显式施加 - 这是变分方法的本质特征

3. 能量守恒： - 弱形式自然保证了能量守恒 - 这对于数值方法的稳定性至关重要

4. 适合数值离散： - 为有限元等数值方法提供了理想的出发点 - 允许使用简单的分片多项式逼近

3.1.2 变分原理与能量最小化

弱形式等价于最小势能原理。这种等价性不仅提供了有限元方法的理论基础，还为数值算法的设计和分析提供了重要工具。

总势能泛函的构造

对于超弹性材料，系统的总势能泛函为： $$\Pi(\mathbf{u}) = \int_\Omega \Psi(\mathbf{F}) \, d\Omega - \int_\Omega \mathbf{u} \cdot \mathbf{f} \, d\Omega - \int_{\Gamma_t} \mathbf{u} \cdot \overline{\mathbf{t}} \, d\Gamma$$ 其中：

• $\Psi(\mathbf{F})$ 是应变能密度函数，依赖于变形梯度 $\mathbf{F} = \mathbf{I} + \nabla \mathbf{u}$
• 第一项是内能（弹性势能）$U = \int_\Omega \Psi \, d\Omega$
• 第二、三项是外力势能 $W = -\int_\Omega \mathbf{u} \cdot \mathbf{f} \, d\Omega - \int_{\Gamma_t} \mathbf{u} \cdot \overline{\mathbf{t}} \, d\Gamma$

总势能可以写作：$\Pi = U + W$

线弹性材料的应变能密度

对于线弹性材料，应变能密度函数为： $$\Psi = \frac{1}{2} \boldsymbol{\varepsilon} : \mathbb{C} : \boldsymbol{\varepsilon} = \frac{1}{2} \lambda (\text{tr}\boldsymbol{\varepsilon})^2 + \mu \boldsymbol{\varepsilon} : \boldsymbol{\varepsilon}$$ 其中：

• $\boldsymbol{\varepsilon} = \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T)$ 是线性（无穷小）应变张量
• $\mathbb{C}$ 是四阶弹性张量
• $\lambda$ 和 $\mu$ 是Lamé常数，与工程常数的关系： $$\lambda = \frac{E\nu}{(1+\nu)(1-2\nu)}, \quad \mu = \frac{E}{2(1+\nu)}$$ 指标形式的应变能密度： $$\Psi = \frac{1}{2} C_{ijkl} \varepsilon_{ij} \varepsilon_{kl}$$ 对于各向同性材料： $$C_{ijkl} = \lambda \delta_{ij} \delta_{kl} + \mu (\delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk})$$

驻点条件与第一变分

平衡位形对应于势能的驻点。设 $\mathbf{u}$ 是真实位移，$\delta \mathbf{u} \in \mathcal{V}$ 是任意容许变分（满足齐次边界条件），则： $$\delta \Pi = \frac{d}{d\epsilon}\Pi(\mathbf{u} + \epsilon \delta \mathbf{u})\bigg|_{\epsilon=0} = 0 \quad \forall \delta \mathbf{u} \in \mathcal{V}$$ 计算第一变分（Gâteaux导数）：

内能的变分： $$\delta U = \int_\Omega \frac{\partial \Psi}{\partial \boldsymbol{\varepsilon}} : \delta \boldsymbol{\varepsilon} \, d\Omega$$ 其中 $\delta \boldsymbol{\varepsilon} = \frac{1}{2}(\nabla \delta \mathbf{u} + \nabla \delta \mathbf{u}^T)$ 是应变的变分。

对于线弹性材料： $$\frac{\partial \Psi}{\partial \boldsymbol{\varepsilon}} = \mathbb{C} : \boldsymbol{\varepsilon} = \boldsymbol{\sigma}$$ 因此： $$\delta U = \int_\Omega \boldsymbol{\sigma} : \delta \boldsymbol{\varepsilon} \, d\Omega = \int_\Omega \boldsymbol{\sigma} : \nabla \delta \mathbf{u} \, d\Omega$$ 这里利用了应力张量的对称性：$\boldsymbol{\sigma} : \delta \boldsymbol{\varepsilon} = \boldsymbol{\sigma} : \nabla \delta \mathbf{u}$。

外力势能的变分： $$\delta W = -\int_\Omega \delta \mathbf{u} \cdot \mathbf{f} \, d\Omega - \int_{\Gamma_t} \delta \mathbf{u} \cdot \overline{\mathbf{t}} \, d\Gamma$$ 总变分： $$\delta \Pi = \int_\Omega \boldsymbol{\sigma} : \nabla \delta \mathbf{u} \, d\Omega - \int_\Omega \delta \mathbf{u} \cdot \mathbf{f} \, d\Omega - \int_{\Gamma_t} \delta \mathbf{u} \cdot \overline{\mathbf{t}} \, d\Gamma = 0$$ 这正是我们之前推导的弱形式！

二阶变分与稳定性

稳定平衡要求势能取极小值，即二阶变分为正： $$\delta^2 \Pi = \int_\Omega \delta \boldsymbol{\varepsilon} : \mathbb{C} : \delta \boldsymbol{\varepsilon} \, d\Omega > 0 \quad \forall \delta \mathbf{u} \neq \mathbf{0}$$ 这等价于弹性张量 $\mathbb{C}$ 的正定性。对于各向同性材料，正定性条件为： $$\mu > 0, \quad 3\lambda + 2\mu > 0$$ 或用工程常数表示： $$E > 0, \quad -1 < \nu < 0.5$$

非线性问题的能量泛函

对于有限变形（几何非线性），应变能密度依赖于变形梯度 $\mathbf{F} = \nabla \mathbf{X} + \nabla \mathbf{u}$： $$\Pi(\mathbf{u}) = \int_{\Omega_0} \Psi(\mathbf{F}) \, dV - \int_{\Omega_0} \mathbf{u} \cdot \mathbf{f}_0 \, dV - \int_{\Gamma_{t0}} \mathbf{u} \cdot \overline{\mathbf{t}}_0 \, dA$$ 注意积分是在参考构型 $\Omega_0$ 上进行的。常用的超弹性模型包括：

Neo-Hookean模型： $$\Psi = \frac{\mu}{2}(I_1 - 3) - \mu \ln J + \frac{\lambda}{2}(\ln J)^2$$ 其中 $I_1 = \text{tr}(\mathbf{F}^T\mathbf{F})$，$J = \det(\mathbf{F})$。

St. Venant-Kirchhoff模型： $$\Psi = \frac{\lambda}{2}(\text{tr}\mathbf{E})^2 + \mu \mathbf{E} : \mathbf{E}$$ 其中 $\mathbf{E} = \frac{1}{2}(\mathbf{F}^T\mathbf{F} - \mathbf{I})$ 是Green-Lagrange应变。

最小原理的唯一性

定理（最小势能原理）：在所有满足位移边界条件的容许位移场中，真实位移场使总势能取最小值。

证明要点：

1. 对于线弹性问题，势能泛函是严格凸的
2. 凸泛函的驻点唯一且为全局最小值
3. 存在性由Lax-Milgram定理保证

重要性质：

1. 线弹性问题：势能泛函是二次的，Hessian矩阵（刚度矩阵）正定，解唯一
2. 非线性问题：可能存在多个驻点（分叉、屈曲），需要稳定性分析
3. 路径无关性：保守系统中，势能只依赖于当前构型，与加载路径无关

互补能原理

除了最小势能原理，还存在基于应力的互补能原理：

余能泛函： $$\Pi^*(\boldsymbol{\sigma}) = \int_\Omega \Psi^*(\boldsymbol{\sigma}) \, d\Omega - \int_{\Gamma_u} (\boldsymbol{\sigma} \cdot \mathbf{n}) \cdot \overline{\mathbf{u}} \, d\Gamma$$ 其中 $\Psi^*$ 是余能密度，通过Legendre变换得到： $$\Psi^*(\boldsymbol{\sigma}) = \boldsymbol{\sigma} : \boldsymbol{\varepsilon} - \Psi(\boldsymbol{\varepsilon})$$ 对于线弹性材料： $$\Psi^* = \frac{1}{2} \boldsymbol{\sigma} : \mathbb{S} : \boldsymbol{\sigma}$$ 其中 $\mathbb{S} = \mathbb{C}^{-1}$ 是柔度张量。

误差估计与收敛性

最小势能原理为有限元误差估计提供了基础：

Céa引理：设 $\mathbf{u}$ 是精确解，$\mathbf{u}^h$ 是有限元解，则： $$|\mathbf{u} - \mathbf{u}^h|_E \leq C \inf_{\mathbf{v}^h \in \mathcal{V}^h} |\mathbf{u} - \mathbf{v}^h|_E$$ 其中 $|\cdot|_E$ 是能量范数： $$|\mathbf{v}|_E^2 = \int_\Omega \nabla \mathbf{v} : \mathbb{C} : \nabla \mathbf{v} \, d\Omega$$ 这表明有限元解是真实解在有限元空间中的最佳逼近（在能量范数意义下）。

3.1.3 离散化与形函数

有限元离散化的核心是用有限维空间逼近无限维函数空间。这个过程涉及网格生成、形函数构造、数值积分和矩阵组装等关键步骤。

网格离散化

将连续域 $\Omega$ 划分为 $n_e$ 个非重叠单元： $$\Omega = \bigcup_{e=1}^{n_e} \Omega_e, \quad \Omega_i \cap \Omega_j = \emptyset \text{ if } i \neq j$$ 网格需要满足的条件：

1. 完备性：所有单元覆盖整个域，无遗漏
2. 相容性：相邻单元在公共边界上的节点匹配
3. 正则性：单元形状不过度扭曲，满足最小角度要求

形函数的构造

在每个单元内，位移场用形函数插值表示： $$\mathbf{u}^h(\mathbf{x}) = \sum_{i=1}^{n} N_i(\mathbf{x}) \mathbf{u}_i$$ 形函数的基本性质：

1. 插值性（Kronecker性质）： $$N_i(\mathbf{x}_j) = \delta_{ij} = \begin{cases} 1 & \text{if } i = j \\ 0 & \text{if } i \neq j \end{cases}$$

2. 单位分解（完备性）： $$\sum_{i=1}^{n} N_i(\mathbf{x}) = 1 \quad \forall \mathbf{x} \in \Omega$$ 这保证了刚体运动的精确表示。

3. 紧支撑性（局部性）： $$\text{supp}(N_i) = \bigcup_{e \in \mathcal{E}_i} \Omega_e$$ 其中 $\mathcal{E}_i$ 是包含节点 $i$ 的单元集合。

4. 连续性要求： - $C^0$ 连续：位移场连续（基本要求） - 高阶连续性：某些问题（如板壳）需要 $C^1$ 或更高

单元形函数的具体构造

一维线性单元（2节点）： $$N_1(\xi) = \frac{1-\xi}{2}, \quad N_2(\xi) = \frac{1+\xi}{2}, \quad \xi \in [-1,1]$$ 二维双线性单元（4节点）： $$N_i(\xi,\eta) = \frac{1}{4}(1 + \xi_i\xi)(1 + \eta_i\eta)$$ 其中 $(\xi_i, \eta_i)$ 是节点 $i$ 在参考单元中的坐标。

拉格朗日多项式基： 对于 $p$ 阶单元，形函数通过拉格朗日插值构造： $$N_i(\xi) = \prod_{j=0,j\neq i}^{p} \frac{\xi - \xi_j}{\xi_i - \xi_j}$$ 分层基函数（Hierarchical basis）： $$N_0 = \frac{1-\xi}{2}, \quad N_1 = \frac{1+\xi}{2}, \quad N_k = \phi_k(\xi) \text{ for } k \geq 2$$ 其中 $\phi_k$ 是高阶多项式，通常选择Legendre多项式。

应变-位移关系的矩阵形式

位移场的离散表示： $$\mathbf{u}^h = \sum_{i=1}^{n} N_i \mathbf{u}_i = \mathbf{N} \mathbf{d}$$ 其中： $$\mathbf{N} = \begin{bmatrix} N_1 & 0 & 0 & N_2 & 0 & 0 & \cdots \\ 0 & N_1 & 0 & 0 & N_2 & 0 & \cdots \\ 0 & 0 & N_1 & 0 & 0 & N_2 & \cdots \end{bmatrix}$$ 应变场的离散表示： $$\boldsymbol{\varepsilon}^h = \mathbf{L} \mathbf{u}^h = \mathbf{L} \mathbf{N} \mathbf{d} = \mathbf{B} \mathbf{d}$$ 其中微分算子矩阵： $$\mathbf{L} = \begin{bmatrix} \frac{\partial}{\partial x} & 0 & 0 \\ 0 & \frac{\partial}{\partial y} & 0 \\ 0 & 0 & \frac{\partial}{\partial z} \\ \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0 \\ 0 & \frac{\partial}{\partial z} & \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} & 0 & \frac{\partial}{\partial x} \end{bmatrix}$$ 应变-位移矩阵 $\mathbf{B} = \mathbf{L} \mathbf{N}$ 的具体形式： $$\mathbf{B}_i = \begin{bmatrix} \frac{\partial N_i}{\partial x} & 0 & 0 \\ 0 & \frac{\partial N_i}{\partial y} & 0 \\ 0 & 0 & \frac{\partial N_i}{\partial z} \\ \frac{\partial N_i}{\partial y} & \frac{\partial N_i}{\partial x} & 0 \\ 0 & \frac{\partial N_i}{\partial z} & \frac{\partial N_i}{\partial y} \\ \frac{\partial N_i}{\partial z} & 0 & \frac{\partial N_i}{\partial x} \end{bmatrix}$$

Galerkin方法与刚度矩阵推导

将离散位移代入弱形式： $$\int_\Omega \nabla \mathbf{v}^h : \boldsymbol{\sigma}^h \, d\Omega = \int_\Omega \mathbf{v}^h \cdot \mathbf{f} \, d\Omega + \int_{\Gamma_t} \mathbf{v}^h \cdot \overline{\mathbf{t}} \, d\Gamma$$ Galerkin方法选择测试函数与试探函数来自同一空间： $$\mathbf{v}^h = \sum_{i=1}^{n} N_i \mathbf{v}_i$$ 由于 $\mathbf{v}_i$ 的任意性，得到离散方程组： $$\sum_{j=1}^{n} K_{ij} u_j = F_i, \quad i = 1,2,...,n$$ 其中刚度矩阵元素： $$K_{ij} = \int_\Omega \mathbf{B}_i^T \mathbf{D} \mathbf{B}_j \, d\Omega$$ 载荷向量元素： $$F_i = \int_\Omega N_i \mathbf{f} \, d\Omega + \int_{\Gamma_t} N_i \overline{\mathbf{t}} \, d\Gamma$$

单元矩阵的计算

在单元级别，刚度矩阵计算为： $$\mathbf{K}^e = \int_{\Omega_e} \mathbf{B}^T \mathbf{D} \mathbf{B} \, d\Omega$$ 对于等参单元，使用坐标变换： $$\mathbf{K}^e = \int_{-1}^{1} \int_{-1}^{1} \int_{-1}^{1} \mathbf{B}^T(\boldsymbol{\xi}) \mathbf{D} \mathbf{B}(\boldsymbol{\xi}) \det(\mathbf{J}(\boldsymbol{\xi})) \, d\xi d\eta d\zeta$$ 材料刚度矩阵 $\mathbf{D}$ 对于各向同性线弹性材料： $$\mathbf{D} = \frac{E}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1-\nu & \nu & \nu & 0 & 0 & 0 \\ \nu & 1-\nu & \nu & 0 & 0 & 0 \\ \nu & \nu & 1-\nu & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} \end{bmatrix}$$

矩阵组装过程

全局矩阵通过单元贡献组装： $$\mathbf{K} = \sum_{e=1}^{n_e} \mathbf{L}_e^T \mathbf{K}^e \mathbf{L}_e$$ 其中 $\mathbf{L}_e$ 是局部-全局自由度映射矩阵。实际实现中，通常使用索引数组：

对于每个单元 e:
    计算单元刚度矩阵 K_e
    获取单元节点的全局编号 [i1, i2, ..., in_en]
    对于局部索引 (a, b):
        全局索引 (I, J) = (3*i_a + α, 3*i_b + β)
        K[I,J] += K_e[3*a+α, 3*b+β]

稀疏性与存储格式

有限元刚度矩阵具有稀疏性：

• 非零元素比例：$O(1/n)$
• 带宽：依赖于节点编号
• 常用存储格式：
• CSR（Compressed Sparse Row）
• CSC（Compressed Sparse Column）
• Skyline（天际线）格式

数值积分考虑

高斯积分点数的选择：

• 完全积分：精确积分多项式阶数 $\geq 2p-1$（$p$ 是形函数阶数）
• 减缩积分：使用较少积分点，可能导致零能模式
• 选择性减缩：对不同项使用不同积分阶数

积分点数与精度的关系： | 单元类型 | 完全积分 | 减缩积分 | 选择性减缩 |

3.2 六面体与四面体网格的FEM实现

网格类型的选择直接影响求解精度、计算效率和实现复杂度。六面体网格在规则几何上表现优异，而四面体网格在处理复杂几何时更加灵活。理解不同单元类型的数学基础和实现细节，是构建高效FEM求解器的关键。

3.2.1 六面体单元

八节点六面体单元（Hex8） 是工程计算中最常用的单元类型之一。其优势在于：精度高、收敛性好、每个自由度的计算成本低。

参考单元与形函数

在参考坐标系 $(\xi, \eta, \zeta) \in [-1,1]^3$ 中，八个节点的位置为： $$\begin{align} &\text{节点1: } (-1, -1, -1), \quad \text{节点2: } (1, -1, -1) \\ &\text{节点3: } (1, 1, -1), \quad \text{节点4: } (-1, 1, -1) \\ &\text{节点5: } (-1, -1, 1), \quad \text{节点6: } (1, -1, 1) \\ &\text{节点7: } (1, 1, 1), \quad \text{节点8: } (-1, 1, 1) \end{align}$$ 三线性形函数定义为： $$N_i(\xi, \eta, \zeta) = \frac{1}{8}(1 + \xi_i\xi)(1 + \eta_i\eta)(1 + \zeta_i\zeta)$$ 形函数的导数（在参考坐标系中）： $$\begin{align} \frac{\partial N_i}{\partial \xi} &= \frac{1}{8}\xi_i(1 + \eta_i\eta)(1 + \zeta_i\zeta) \\ \frac{\partial N_i}{\partial \eta} &= \frac{1}{8}(1 + \xi_i\xi)\eta_i(1 + \zeta_i\zeta) \\ \frac{\partial N_i}{\partial \zeta} &= \frac{1}{8}(1 + \xi_i\xi)(1 + \eta_i\eta)\zeta_i \end{align}$$

等参变换与雅可比矩阵

物理坐标和参考坐标通过等参变换关联： $$\mathbf{x}(\xi, \eta, \zeta) = \sum_{i=1}^{8} N_i(\xi, \eta, \zeta) \mathbf{x}_i$$ 雅可比矩阵定义坐标变换的局部线性近似： $$\mathbf{J} = \frac{\partial \mathbf{x}}{\partial \boldsymbol{\xi}} = \begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} & \frac{\partial z}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} & \frac{\partial z}{\partial \eta} \\ \frac{\partial x}{\partial \zeta} & \frac{\partial y}{\partial \zeta} & \frac{\partial z}{\partial \zeta} \end{bmatrix} = \sum_{i=1}^{8} \mathbf{x}_i \otimes \nabla_{\boldsymbol{\xi}} N_i$$ 形函数在物理坐标中的导数通过链式法则计算： $$\nabla_{\mathbf{x}} N_i = \mathbf{J}^{-T} \nabla_{\boldsymbol{\xi}} N_i$$

数值积分

六面体单元通常使用 $2 \times 2 \times 2$ 高斯积分（完全积分）： $$\int_{\Omega_e} f(\mathbf{x}) \, d\Omega = \int_{-1}^{1} \int_{-1}^{1} \int_{-1}^{1} f(\mathbf{x}(\boldsymbol{\xi})) \det(\mathbf{J}) \, d\xi d\eta d\zeta \approx \sum_{g=1}^{8} w_g f(\mathbf{x}_g) \det(\mathbf{J}_g)$$ 高斯点位置：$\xi_g, \eta_g, \zeta_g = \pm 1/\sqrt{3}$，权重：$w_g = 1$。

单元质量指标：

• 雅可比行列式：$\det(\mathbf{J}) > 0$（保证单元不翻转）
• 长宽比：评估单元的伸长程度
• 扭曲度：评估单元面的平面性偏离

3.2.2 四面体单元

四节点四面体单元（Tet4） 是自动网格生成的首选，特别适合复杂几何。虽然精度低于六面体，但其简单性和灵活性使其在许多应用中不可或缺。

体积坐标系

四面体的体积坐标（重心坐标）$L_i$ 定义为： $$L_i = \frac{V_i}{V}, \quad i = 1,2,3,4$$ 其中 $V_i$ 是由点 $\mathbf{x}$ 和对面（去掉节点 $i$ 的三角形）构成的四面体体积，$V$ 是整个四面体体积。

体积计算公式： $$V = \frac{1}{6}\det\begin{bmatrix} x_1 & y_1 & z_1 & 1 \\ x_2 & y_2 & z_2 & 1 \\ x_3 & y_3 & z_3 & 1 \\ x_4 & y_4 & z_4 & 1 \end{bmatrix}$$ 线性形函数恰好等于体积坐标：$N_i = L_i$。

形函数导数

体积坐标的导数可以通过几何关系直接计算： $$\nabla L_i = -\frac{\mathbf{n}_i A_i}{3V}$$ 其中 $\mathbf{n}_i$ 是对面的外法向量，$A_i$ 是对面的面积。

对于线性四面体，应变是常数，因此 $\mathbf{B}$ 矩阵在整个单元内保持不变： $$\mathbf{B} = \frac{1}{6V} \begin{bmatrix} b_1 & 0 & 0 & b_2 & 0 & 0 & b_3 & 0 & 0 & b_4 & 0 & 0 \\ 0 & c_1 & 0 & 0 & c_2 & 0 & 0 & c_3 & 0 & 0 & c_4 & 0 \\ 0 & 0 & d_1 & 0 & 0 & d_2 & 0 & 0 & d_3 & 0 & 0 & d_4 \\ c_1 & b_1 & 0 & c_2 & b_2 & 0 & c_3 & b_3 & 0 & c_4 & b_4 & 0 \\ 0 & d_1 & c_1 & 0 & d_2 & c_2 & 0 & d_3 & c_3 & 0 & d_4 & c_4 \\ d_1 & 0 & b_1 & d_2 & 0 & b_2 & d_3 & 0 & b_3 & d_4 & 0 & b_4 \end{bmatrix}$$ 其中： $$b_i = (-1)^i(y_j - y_k)(z_k - z_l) + (y_k - y_l)(z_j - z_k)$$ $$c_i = (-1)^i(z_j - z_k)(x_k - x_l) + (z_k - z_l)(x_j - x_k)$$ $$d_i = (-1)^i(x_j - x_k)(y_k - y_l) + (x_k - x_l)(y_j - y_k)$$ 这里 $(i,j,k,l)$ 是 $(1,2,3,4)$ 的循环排列。

3.2.3 锁定现象与缓解策略

体积锁定是不可压缩或近不可压缩材料（$\nu \to 0.5$）的主要数值问题。当泊松比接近0.5时，体积模量 $K = E/[3(1-2\nu)]$ 趋于无穷，导致体积应变被过度约束。

锁定的数学机理

对于完全积分的低阶单元，体积约束条件： $$\text{tr}(\boldsymbol{\varepsilon}) = \nabla \cdot \mathbf{u} = 0$$ 在单元级别施加了过多的约束，导致位移场过度刚化。以4节点四边形为例，每个积分点施加一个约束，4个积分点共4个约束，而单元只有8个自由度，扣除刚体运动后只剩2个变形模式。

选择性减缩积分（B-bar方法）

将应变分解为体积和偏差部分： $$\boldsymbol{\varepsilon} = \frac{1}{3}\text{tr}(\boldsymbol{\varepsilon})\mathbf{I} + \text{dev}(\boldsymbol{\varepsilon})$$ 对应的应力分解： $$\boldsymbol{\sigma} = K \text{tr}(\boldsymbol{\varepsilon}) \mathbf{I} + 2\mu \text{dev}(\boldsymbol{\varepsilon})$$ B-bar方法的核心思想是对体积应变使用单元平均值： $$\overline{\mathbf{B}} = \mathbf{B}_{\text{dev}} + \frac{1}{3}\overline{\mathbf{B}_{\text{vol}}} \mathbf{m} \mathbf{m}^T$$ 其中 $\mathbf{m} = [1, 1, 1, 0, 0, 0]^T$，$\overline{\mathbf{B}_{\text{vol}}}$ 是体积部分的单元平均。

混合格式

引入压力 $p$ 作为独立变量，形成混合变分原理： $$\begin{align} \int_\Omega \text{dev}(\boldsymbol{\varepsilon}) : 2\mu \text{dev}(\boldsymbol{\varepsilon}) \, d\Omega &+ \int_\Omega p \nabla \cdot \mathbf{v} \, d\Omega = \int_\Omega \mathbf{v} \cdot \mathbf{f} \, d\Omega \\ \int_\Omega q \left(\nabla \cdot \mathbf{u} - \frac{p}{K}\right) \, d\Omega &= 0 \end{align}$$ 这导致鞍点问题： $$\begin{bmatrix} \mathbf{A} & \mathbf{B}^T \\ \mathbf{B} & -\mathbf{C} \end{bmatrix} \begin{bmatrix} \mathbf{U} \\ \mathbf{P} \end{bmatrix} = \begin{bmatrix} \mathbf{F} \\ \mathbf{0} \end{bmatrix}$$ 稳定性要求满足inf-sup条件（LBB条件）。

3.3 边界条件处理技巧

边界条件的正确施加是有限元分析成功的关键。不当的处理可能导致方程组奇异、收敛困难或精度损失。本节详细讨论各类边界条件的数学基础和实现技巧。

3.3.1 位移边界条件

位移边界条件（Dirichlet条件）要求某些自由度取预定值。从数学角度看，这是对解空间的约束，需要在保持系统对称性和稳定性的前提下施加。

直接置换法（行列消元法）

最直观的方法是直接修改线性系统。对于约束 $u_i = \overline{u}_i$：

1. 保持对称性的实现：

对于第i个约束自由度：

- 将K的第i行和第i列移到右端项
- 设置 K[i,i] = 1, K[i,j] = K[j,i] = 0 (j≠i)
- 修改载荷：F[j] = F[j] - K[j,i] * ū[i]
- 设置 F[i] = ū[i]

2. 矩阵形式： 原系统： $$\begin{bmatrix} K_{aa} & K_{ab} \\ K_{ba} & K_{bb} \end{bmatrix} \begin{bmatrix} U_a \\ U_b \end{bmatrix} = \begin{bmatrix} F_a \\ F_b \end{bmatrix}$$ 其中 $U_b = \overline{U}_b$ 是约束自由度。修改后： $$K_{aa} U_a = F_a - K_{ab} \overline{U}_b$$

罚函数法

罚函数法通过在能量泛函中添加惩罚项来近似满足约束： $$\Pi_{\alpha}(\mathbf{u}) = \Pi(\mathbf{u}) + \frac{\alpha}{2} \sum_i (u_i - \overline{u}_i)^2$$ 导致修改后的刚度矩阵和载荷向量： $$\tilde{K}_{ii} = K_{ii} + \alpha, \quad \tilde{F}_i = F_i + \alpha \overline{u}_i$$ 参数选择：

• $\alpha$ 太小：约束满足不精确
• $\alpha$ 太大：矩阵条件数恶化
• 经验选择：$\alpha = 10^3 \sim 10^6 \times \max(K_{ii})$

优势：

• 保持矩阵规模不变
• 实现简单，易于并行化
• 适合处理接触等非线性约束

拉格朗日乘子法

引入拉格朗日乘子 $\lambda$ 作为约束反力： $$L(\mathbf{u}, \boldsymbol{\lambda}) = \Pi(\mathbf{u}) + \sum_i \lambda_i (u_i - \overline{u}_i)$$ 导致扩展的鞍点系统： $$\begin{bmatrix} \mathbf{K} & \mathbf{C}^T \\ \mathbf{C} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{U} \\ \boldsymbol{\lambda} \end{bmatrix} = \begin{bmatrix} \mathbf{F} \\ \mathbf{g} \end{bmatrix}$$ 其中 $\mathbf{C}$ 是约束矩阵，$\mathbf{g}$ 是约束值向量。

数值考虑：

• 鞍点系统需要特殊求解器
• 零对角块可能导致数值不稳定
• 可使用增广拉格朗日法改善条件数

3.3.2 周期性边界条件

周期性边界条件在多尺度分析、代表性体积元（RVE）计算和晶体塑性模拟中广泛应用。

数学描述

对于周期性单胞，位移场分解为： $$\mathbf{u}(\mathbf{x}) = \overline{\boldsymbol{\varepsilon}} \cdot \mathbf{x} + \tilde{\mathbf{u}}(\mathbf{x})$$ 其中 $\overline{\boldsymbol{\varepsilon}}$ 是宏观应变，$\tilde{\mathbf{u}}$ 是周期性波动。

周期性条件要求： $$\tilde{\mathbf{u}}(\mathbf{x} + \mathbf{L}_i) = \tilde{\mathbf{u}}(\mathbf{x})$$ 对于相对的边界节点对 $(+,-)$： $$\mathbf{u}^+ - \mathbf{u}^- = \overline{\boldsymbol{\varepsilon}} \cdot (\mathbf{x}^+ - \mathbf{x}^-)$$

实现策略

1. 主从节点法： - 选择一侧为主节点，另一侧为从节点 - 从节点自由度表示为：$\mathbf{u}^- = \mathbf{u}^+ - \Delta\mathbf{u}$ - 在组装时累加从节点贡献到主节点

2. 约束方程法： 使用拉格朗日乘子或罚函数施加线性约束： $$\mathbf{C}_{\text{per}} \mathbf{U} = \mathbf{g}_{\text{per}}$$

3. 角点处理： - 角点同时属于多个面，需要特殊处理 - 通常选择一个角点作为参考，其他角点相对于它约束

宏观载荷施加

施加宏观应变 $\overline{\boldsymbol{\varepsilon}}$： $$\mathbf{F}_{\text{macro}} = \int_{\partial\Omega} \mathbf{t} \otimes \mathbf{x} \, dS$$ 通过Hill-Mandel条件保证能量一致性： $$\overline{\boldsymbol{\sigma}} : \overline{\boldsymbol{\varepsilon}} = \frac{1}{V} \int_\Omega \boldsymbol{\sigma} : \boldsymbol{\varepsilon} \, dV$$

3.3.3 接触边界条件

接触是典型的单边约束问题，涉及几何非线性和状态切换。

Hertz-Signorini-Moreau条件

接触条件的数学表述： $$\begin{cases} g_N \geq 0 & \text{（间隙非负）} \\ t_N \leq 0 & \text{（压力非正）} \\ g_N \cdot t_N = 0 & \text{（互补条件）} \end{cases}$$ 对于摩擦接触，还需满足Coulomb定律： $$\begin{cases} |\mathbf{t}_T| \leq \mu |t_N| & \text{（粘着）} \\ |\mathbf{t}_T| = \mu |t_N| \Rightarrow \exists \lambda > 0: \dot{\mathbf{g}}_T = -\lambda \mathbf{t}_T & \text{（滑动）} \end{cases}$$

数值方法

1. 罚函数法： $$t_N = \begin{cases} -\epsilon_N g_N & \text{if } g_N < 0 \\ 0 & \text{otherwise} \end{cases}$$ 切向采用弹性-滑动模型： $$\mathbf{t}_T = \begin{cases} -\epsilon_T \mathbf{g}_T & \text{if } |\epsilon_T \mathbf{g}_T| < \mu |t_N| \\ -\mu |t_N| \frac{\mathbf{g}_T}{|\mathbf{g}_T|} & \text{otherwise} \end{cases}$$

2. 增广拉格朗日法： 结合拉格朗日乘子和罚函数： $$t_N = \lambda_N - \epsilon_N g_N$$ 通过Uzawa算法更新乘子： $$\lambda_N^{k+1} = \max(0, \lambda_N^k - \epsilon_N g_N^k)$$

3. Active Set策略： - 预测接触状态（接触/分离） - 求解约束问题 - 检查违反的约束并更新活动集 - 迭代直到收敛

接触搜索

高效的接触检测对大规模问题至关重要：

1. 空间划分： - 包围盒层次（BVH） - 八叉树/四叉树 - 空间哈希

2. 主从面识别： - 基于法向的投影 - 最近点搜索 - 保守的穿透检测

3. 接触积分： - 节点-面接触：简单但可能不守恒 - 面-面接触：守恒但实现复杂 - Mortar方法：数学严格但计算昂贵

3.4 可逆有限元与隐式求解器

大变形问题需要特殊处理以保证数值稳定性。

3.4.1 可逆性约束

当 $\det(\mathbf{F}) \leq 0$ 时，标准本构模型失效。可逆有限元通过修改应变能函数保证物理合理性： $$\Psi_{\text{invertible}}(\mathbf{F}) = \begin{cases} \Psi(\mathbf{F}) & \text{if } \det(\mathbf{F}) > \delta \\ \Psi(\mathbf{F}_{\text{proj}}) + \text{barrier} & \text{otherwise} \end{cases}$$ 其中 $\mathbf{F}_{\text{proj}}$ 是投影到可逆空间的变形梯度。

3.4.2 隐式时间积分

对于动力学问题，隐式Newmark方法： $$\mathbf{u}_{n+1} = \mathbf{u}_n + \Delta t \mathbf{v}_n + \frac{\Delta t^2}{2}[(1-2\beta)\mathbf{a}_n + 2\beta\mathbf{a}_{n+1}]$$ $$\mathbf{v}_{n+1} = \mathbf{v}_n + \Delta t[(1-\gamma)\mathbf{a}_n + \gamma\mathbf{a}_{n+1}]$$ 导致非线性系统： $$\mathbf{R}(\mathbf{u}_{n+1}) = \mathbf{M}\mathbf{a}_{n+1} + \mathbf{f}_{\text{int}}(\mathbf{u}_{n+1}) - \mathbf{f}_{\text{ext}} = \mathbf{0}$$

3.4.3 Newton-Raphson求解

线性化残差： $$\mathbf{R}(\mathbf{u} + \Delta\mathbf{u}) \approx \mathbf{R}(\mathbf{u}) + \frac{\partial \mathbf{R}}{\partial \mathbf{u}}\Delta\mathbf{u} = \mathbf{0}$$ 切线刚度矩阵： $$\mathbf{K}_T = \frac{\partial \mathbf{R}}{\partial \mathbf{u}} = \frac{1}{\beta\Delta t^2}\mathbf{M} + \frac{\partial \mathbf{f}_{\text{int}}}{\partial \mathbf{u}}$$ 收敛准则通常基于残差范数和位移增量。

3.4.4 线性求解器选择

大规模系统的求解策略：

1. 直接求解器（如Cholesky分解）： - 适用于中等规模问题 - 内存需求 $O(n^{3/2})$（3D问题）

2. 迭代求解器（如共轭梯度法）： - 适用于大规模问题 - 需要良好的预条件子

常用预条件包括Jacobi、不完全Cholesky分解和代数多重网格。

3.5 拓扑优化应用

拓扑优化寻找材料的最优分布以满足设计目标。

3.5.1 密度法（SIMP）

固体各向同性材料惩罚（SIMP）方法引入密度变量 $\rho \in [0,1]$： $$E(\rho) = \rho^p E_0$$ 其中 $p > 1$ 是惩罚参数，推动设计向0-1分布收敛。

优化问题表述： $$\begin{aligned} \min_{\rho} \quad & c = \mathbf{U}^T \mathbf{K}(\rho) \mathbf{U} \\ \text{s.t.} \quad & \mathbf{K}(\rho) \mathbf{U} = \mathbf{F} \\ & \sum_e v_e \rho_e \leq V^* \\ & 0 < \rho_{\min} \leq \rho_e \leq 1 \end{aligned}$$

3.5.2 灵敏度分析

目标函数对密度的导数： $$\frac{\partial c}{\partial \rho_e} = -p\rho_e^{p-1} \mathbf{u}_e^T \mathbf{k}_0 \mathbf{u}_e$$ 其中 $\mathbf{u}_e$ 是单元位移向量。

3.5.3 过滤技术

防止棋盘格现象的密度过滤： $$\tilde{\rho}_e = \frac{\sum_i H_{ei} v_i \rho_i}{\sum_i H_{ei} v_i}$$ 其中 $H_{ei}$ 是基于距离的权重函数。

3.5.4 优化算法

常用方法包括：

1. 优化准则法（OC）：基于KKT条件的启发式更新
2. 移动渐近线法（MMA）：凸近似的序列规划
3. 水平集方法：隐式表示拓扑边界

本章小结

本章深入探讨了有限元方法的理论基础和实现技巧：

1. 弱形式推导：从强形式出发，通过变分原理得到弱形式，奠定FEM的数学基础
2. 网格策略：六面体网格精度高但生成困难，四面体网格灵活但需要更多单元
3. 边界条件：掌握位移、力和接触边界的处理方法对求解稳定性至关重要
4. 大变形分析：可逆有限元和隐式积分保证了非线性问题的鲁棒性
5. 拓扑优化：展示了FEM在设计优化中的强大应用

关键公式回顾：

• 弱形式：$\int_\Omega \nabla \mathbf{v} : \boldsymbol{\sigma} \, d\Omega = \int_\Omega \mathbf{v} \cdot \mathbf{f} \, d\Omega + \int_{\Gamma_t} \mathbf{v} \cdot \overline{\mathbf{t}} \, d\Gamma$
• 刚度矩阵：$K_{ij} = \int_\Omega \mathbf{B}_i^T \mathbf{D} \mathbf{B}_j \, d\Omega$
• Newton-Raphson：$\mathbf{K}_T \Delta\mathbf{u} = -\mathbf{R}$
• SIMP插值：$E(\rho) = \rho^p E_0$

练习题

基础题

3.1 推导一维杆单元的刚度矩阵。考虑长度为 $L$、横截面积为 $A$、弹性模量为 $E$ 的杆单元，使用线性形函数。

3.2 证明四面体单元的体积坐标满足： $$L_1 + L_2 + L_3 + L_4 = 1$$ $$\nabla L_1 + \nabla L_2 + \nabla L_3 + \nabla L_4 = \mathbf{0}$$

3.3 对于平面应力问题，推导选择性减缩积分如何缓解体积锁定。考虑4节点双线性单元。

挑战题

3.4 设计一个算法，自动检测六面体网格中的扭曲单元。考虑雅可比行列式的正负性和条件数。

3.5 推导增广拉格朗日方法处理接触问题的迭代格式。考虑法向接触with摩擦。

3.6 分析拓扑优化中continuation方法的作用。如何选择SIMP惩罚参数 $p$ 的演化策略？

3.7 设计一个自适应网格细化策略用于有限元分析。基于后验误差估计器。

3.8 探讨无网格法（如MLS）相比FEM的优劣。什么情况下应该选择无网格方法？

常见陷阱与错误

1. 网格质量问题

• 错误：忽视网格质量检查，导致负雅可比
• 后果：求解发散或得到非物理结果
• 解决：实施网格质量度量，拒绝劣质单元

2. 数值积分精度

• 错误：对所有项使用相同的积分阶数
• 后果：可能导致锁定或零能模式
• 解决：针对不同项选择合适的积分方案

3. 边界条件冲突

• 错误：在同一节点同时施加位移和力边界条件
• 后果：方程组奇异
• 解决：仔细检查边界条件的相容性

4. 大变形线性化

• 错误：在大变形问题中使用小变形假设
• 后果：严重的精度损失
• 解决：使用几何非线性公式和适当的应力度量

5. 迭代求解器发散

• 错误：不当的预条件子或收敛准则
• 后果：计算时间过长或无法收敛
• 解决：选择问题相关的预条件子，调整收敛容差

最佳实践检查清单

前处理阶段

• [ ] 网格质量检查（雅可比、长宽比、扭曲度）
• [ ] 材料参数物理合理性验证
• [ ] 边界条件完整性和相容性检查
• [ ] 单位制一致性确认

求解器配置

• [ ] 根据问题特征选择合适的单元类型
• [ ] 非线性问题的载荷步长合理划分
• [ ] 选择适当的收敛准则（力/位移/能量）
• [ ] 线性求解器和预条件子匹配问题规模

结果验证

• [ ] 能量守恒检查（适用时）
• [ ] 应力集中区域网格收敛性研究
• [ ] 与解析解或基准解对比（如有）
• [ ] 物理合理性判断（变形模式、应力分布）

性能优化

• [ ] 利用对称性减少计算域
• [ ] 合理的自由度编号减少带宽
• [ ] 并行化潜力评估
• [ ] 存储格式优化（CSR、CSC等）

拓扑优化特定

• [ ] 体积约束合理性
• [ ] 过滤半径与网格尺寸匹配
• [ ] 惩罚参数演化策略
• [ ] 灰度单元处理方案