第三章：拉格朗日视角（2）：有限元仿真

有限元方法(FEM)是物理仿真中最重要的数值方法之一，广泛应用于固体力学、流体力学和传热等领域。本章将深入探讨有限元方法的理论基础，从弱形式推导开始，逐步介绍材料模型、网格类型和高级特性。通过学习本章内容，读者将掌握使用有限元方法进行弹性体仿真的核心技术，并了解拓扑优化等前沿应用。

3.1 弱形式与拉格朗日有限元入门

有限元方法的核心思想是将连续的偏微分方程问题转化为离散的线性代数问题。这个转化过程的关键是弱形式(weak form)，它不仅降低了对解的光滑性要求，还自然地引入了边界条件的处理方式。

3.1.1 强形式vs弱形式

在物理仿真中，我们经常需要求解偏微分方程(PDE)。以泊松方程为例，它描述了稳态热传导、静电势等物理现象：

强形式： $$\nabla \cdot \nabla u = f \quad \text{in } \Omega$$ $$u = g \quad \text{on } \partial\Omega_D \quad \text{(Dirichlet边界)}$$ $$\nabla u \cdot n = h \quad \text{on } \partial\Omega_N \quad \text{(Neumann边界)}$$ 其中：

• $\Omega$ 是求解域
• $\partial\Omega = \partial\Omega_D \cup \partial\Omega_N$ 是边界
• $n$ 是外法向量
• $f$ 是源项（如热源密度）
• $g$ 是已知的边界值（如温度）
• $h$ 是已知的边界通量（如热流）

强形式要求解在每个点都满足微分方程，这意味着：

1. 解必须二阶可微（$C^2$连续）
2. 在不连续点（如材料界面）难以处理
3. 数值离散化不直观

弱形式推导： 引入测试函数（也称权函数）$w \in V$，其中 $V = \{w : w \in H^1(\Omega), w = 0 \text{ on } \partial\Omega_D\}$。

将强形式两边同时乘以测试函数并在域上积分： $$\int_\Omega w(\nabla \cdot \nabla u) \, d\Omega = \int_\Omega wf \, d\Omega$$ 这就是弱形式的起点。注意到左边包含了$u$的二阶导数，我们将通过分部积分降低导数阶数。

3.1.2 分部积分与散度定理

分部积分是将强形式转化为弱形式的关键工具。回顾一维的分部积分公式： $$\int_a^b u'v \, dx = [uv]_a^b - \int_a^b uv' \, dx$$ 在多维情况下，我们使用散度定理（也称高斯定理）： $$\int_\Omega \nabla \cdot \mathbf{F} \, d\Omega = \int_{\partial\Omega} \mathbf{F} \cdot n \, d\Gamma$$ 详细推导过程： 从 $\int_\Omega w(\nabla \cdot \nabla u) \, d\Omega$ 开始，注意到 $\nabla \cdot \nabla u = \nabla \cdot (\nabla u)$。

使用恒等式： $$\nabla \cdot (w\nabla u) = \nabla w \cdot \nabla u + w(\nabla \cdot \nabla u)$$ 重排得： $$w(\nabla \cdot \nabla u) = \nabla \cdot (w\nabla u) - \nabla w \cdot \nabla u$$ 代入积分： $$\int_\Omega w(\nabla \cdot \nabla u) \, d\Omega = \int_\Omega \nabla \cdot (w\nabla u) \, d\Omega - \int_\Omega \nabla w \cdot \nabla u \, d\Omega$$ 应用散度定理于第一项： $$\int_\Omega w(\nabla \cdot \nabla u) \, d\Omega = \int_{\partial\Omega} w \nabla u \cdot n \, d\Gamma - \int_\Omega \nabla w \cdot \nabla u \, d\Omega$$ 处理边界积分： 边界积分需要根据边界类型分别处理： $$\int_{\partial\Omega} w \nabla u \cdot n \, d\Gamma = \int_{\partial\Omega_D} w \nabla u \cdot n \, d\Gamma + \int_{\partial\Omega_N} w \nabla u \cdot n \, d\Gamma$$ 由于测试函数在Dirichlet边界上为零（$w = 0$ on $\partial\Omega_D$），第一项消失。 在Neumann边界上，我们有 $\nabla u \cdot n = h$，因此： $$\int_{\partial\Omega} w \nabla u \cdot n \, d\Gamma = \int_{\partial\Omega_N} wh \, d\Gamma$$ 弱形式的最终形式： $$\int_\Omega \nabla w \cdot \nabla u \, d\Omega = \int_\Omega wf \, d\Omega + \int_{\partial\Omega_N} wh \, d\Gamma$$ 这就是泊松方程的弱形式，也称为变分形式。注意到：

1. 只需要$u$和$w$的一阶导数（$H^1$空间）
2. Neumann边界条件自然地出现在弱形式中（自然边界条件）
3. Dirichlet边界条件通过测试函数空间强制满足（本质边界条件）

3.1.3 Galerkin方法

Galerkin方法是将无限维函数空间问题转化为有限维线性代数问题的一般框架。有限元方法是Galerkin方法的一个特例，其特点是使用具有局部支撑的基函数。

Galerkin离散化的步骤：

1. 构造有限维子空间： 选择有限维试函数空间 $V_h \subset V$： $$V_h = \text{span}\{\phi_1, \phi_2, ..., \phi_n\}$$ 其中$h$表示网格尺寸，$\phi_i$是基函数（也称形函数）。

2. 解的离散表示： 将近似解$u_h$表示为基函数的线性组合： $$u_h(x) = \sum_{j=1}^n u_j \phi_j(x)$$ 其中$u_j$是待求的系数（节点值）。

3. Galerkin投影： 要求残差与所有测试函数正交。选择测试函数 $w_h = \phi_i$（$i = 1, 2, ..., n$），代入弱形式： $$\int_\Omega \nabla \phi_i \cdot \nabla u_h \, d\Omega = \int_\Omega \phi_i f \, d\Omega + \int_{\partial\Omega_N} \phi_i h \, d\Gamma$$

4. 线性方程组： 将$u_h$的表达式代入： $$\sum_{j=1}^n u_j \int_\Omega \nabla \phi_i \cdot \nabla \phi_j \, d\Omega = \int_\Omega \phi_i f \, d\Omega + \int_{\partial\Omega_N} \phi_i h \, d\Gamma$$ 定义：

• 刚度矩阵：$K_{ij} = \int_\Omega \nabla \phi_i \cdot \nabla \phi_j \, d\Omega$
• 载荷向量：$f_i = \int_\Omega \phi_i f \, d\Omega + \int_{\partial\Omega_N} \phi_i h \, d\Gamma$

得到线性系统： $$\mathbf{K}\mathbf{u} = \mathbf{f}$$ Galerkin方法的数学性质：

• 最佳近似性：Galerkin解是真解在$V_h$中的最佳近似（在能量范数意义下）
• 对称性：如果原问题是自伴的，则刚度矩阵对称
• 正定性：对于椭圆型问题，刚度矩阵正定

3.1.4 试函数与测试函数

有限元方法的一个关键特征是使用具有局部支撑的分片多项式基函数。这些基函数的选择直接影响了方法的精度、稳定性和计算效率。

基函数的要求：

1. 完备性：能够精确表示常数和线性函数（对于二阶PDE）
2. 连续性：满足所需的连续性要求（$C^0$对于二阶PDE）
3. 局部支撑：每个基函数只在少数单元上非零
4. 插值性：$\phi_i(x_j) = \delta_{ij}$（Kronecker delta）

一维基函数：

1. 线性基函数（帽子函数）： $$\phi_i(x) = \begin{cases} \frac{x - x_{i-1}}{x_i - x_{i-1}} & x_{i-1} \leq x \leq x_i \\ \frac{x_{i+1} - x}{x_{i+1} - x_i} & x_i \leq x \leq x_{i+1} \\ 0 & \text{otherwise} \end{cases}$$ 性质：

• 在节点$x_i$处值为1，在其他节点处为0
• 支撑域为$[x_{i-1}, x_{i+1}]$
• 导数为分片常数

2. 二次基函数： 需要在单元中点引入额外节点。对于标准单元$[-1, 1]$： $$N_1(\xi) = \frac{\xi(\xi-1)}{2}, \quad N_2(\xi) = 1-\xi^2, \quad N_3(\xi) = \frac{\xi(\xi+1)}{2}$$

3. Hermite基函数： 同时插值函数值和导数值，用于需要$C^1$连续的问题： $$H_i^0(x) = 1 - 3\xi^2 + 2\xi^3, \quad H_i^1(x) = \xi - 2\xi^2 + \xi^3$$ 其中$\xi = (x-x_i)/(x_{i+1}-x_i)$。

二维基函数：

1. 三角形单元的线性基函数： 使用重心坐标$(L_1, L_2, L_3)$： $$\phi_i = L_i, \quad i = 1, 2, 3$$ 其中$L_i$是第$i$个重心坐标，满足$L_1 + L_2 + L_3 = 1$。

2. 四边形单元的双线性基函数： 在参考单元$[-1,1]^2$上： $$N_i(\xi, \eta) = \frac{1}{4}(1 + \xi_i\xi)(1 + \eta_i\eta)$$ 其中$(\xi_i, \eta_i)$是角点坐标。

p型与h型细化：

• h型细化：减小单元尺寸$h$，保持多项式阶数$p$固定
• p型细化：增加多项式阶数$p$，保持网格固定
• hp型细化：同时调整$h$和$p$，达到指数收敛率

基函数的数值性质：

• 条件数：高阶基函数可能导致刚度矩阵条件数增大
• 稀疏性：局部支撑保证刚度矩阵稀疏
• 正交性：某些基函数（如Legendre多项式）具有正交性，改善数值性质

3.2 变形与弹性基础

连续介质力学是有限元仿真的理论基础。本节将介绍描述材料变形的基本概念，包括变形梯度、应变度量和应力张量，以及它们之间的关系。这些概念对于理解和实现非线性有限元至关重要。

3.2.1 变形梯度张量F

变形梯度张量$\mathbf{F}$是描述材料局部变形的最基本量，它建立了参考构型和当前构型之间的联系。

定义与几何意义： 考虑一个连续体从参考构型（初始/未变形状态）到当前构型（变形后状态）的映射： $$\mathbf{x} = \boldsymbol{\chi}(\mathbf{X}, t)$$ 其中：

• $\mathbf{X} \in \Omega_0$：参考构型中的材料点位置
• $\mathbf{x} \in \Omega$：当前构型中的材料点位置
• $\boldsymbol{\chi}$：变形映射
• $t$：时间参数

变形梯度定义为： $$\mathbf{F} = \frac{\partial \mathbf{x}}{\partial \mathbf{X}} = \frac{\partial \boldsymbol{\chi}}{\partial \mathbf{X}}$$ 在指标记号下： $$F_{ij} = \frac{\partial x_i}{\partial X_j}$$ 与位移的关系： 引入位移场 $\mathbf{u}(\mathbf{X}, t) = \mathbf{x} - \mathbf{X}$，则： $$\mathbf{F} = \mathbf{I} + \frac{\partial \mathbf{u}}{\partial \mathbf{X}} = \mathbf{I} + \nabla_0 \mathbf{u}$$ 其中$\nabla_0$表示相对于参考坐标的梯度算子。

物理解释： 变形梯度描述了材料线元的变形：

• 参考构型中的线元：$d\mathbf{X}$
• 当前构型中的线元：$d\mathbf{x} = \mathbf{F} d\mathbf{X}$

变形梯度的重要性质：

1. 可逆性条件： $$J = \det(\mathbf{F}) > 0$$ $J$称为Jacobian或体积比，表示局部体积变化：

• $J > 1$：体积膨胀
• $J < 1$：体积压缩
• $J = 1$：体积不变（不可压缩）
• $J \leq 0$：材料穿透或翻转（非物理）

2. 极分解： 任何可逆的变形梯度可以唯一分解为： $$\mathbf{F} = \mathbf{R}\mathbf{U} = \mathbf{V}\mathbf{R}$$ 其中：

• $\mathbf{R}$：旋转张量（正交张量）
• $\mathbf{U}$：右拉伸张量（对称正定）
• $\mathbf{V}$：左拉伸张量（对称正定）

3. 客观性（框架不变性）： 在刚体运动下，变形梯度变换为： $$\mathbf{F}^* = \mathbf{Q}\mathbf{F}$$ 其中$\mathbf{Q}$是旋转张量。这要求本构关系必须满足客观性。

4. 特殊变形模式： - 刚体平移：$\mathbf{F} = \mathbf{I}$ - 刚体旋转：$\mathbf{F} = \mathbf{R}$（正交张量） - 单轴拉伸：$\mathbf{F} = \text{diag}(\lambda, 1/\sqrt{\lambda}, 1/\sqrt{\lambda})$（不可压缩情况） - 简单剪切：$\mathbf{F} = \mathbf{I} + \gamma \mathbf{e}_1 \otimes \mathbf{e}_2$

3.2.2 格林应变与柯西应变

应变是描述变形程度的无量纲量。由于有限变形下应变的定义不唯一，我们需要了解各种应变度量及其适用场景。

为什么需要应变度量： 变形梯度$\mathbf{F}$包含了刚体旋转，不是真正的应变度量。我们需要从$\mathbf{F}$中提取出纯变形信息。

右柯西-格林变形张量： $$\mathbf{C} = \mathbf{F}^T \mathbf{F}$$ 性质：

• 对称张量：$\mathbf{C} = \mathbf{C}^T$
• 客观性：在刚体旋转下不变
• 主不变量：$I_1 = \text{tr}(\mathbf{C})$，$I_2 = \frac{1}{2}[(\text{tr}\mathbf{C})^2 - \text{tr}(\mathbf{C}^2)]$，$I_3 = \det(\mathbf{C}) = J^2$

格林-拉格朗日应变张量： $$\mathbf{E} = \frac{1}{2}(\mathbf{C} - \mathbf{I}) = \frac{1}{2}(\mathbf{F}^T \mathbf{F} - \mathbf{I})$$ 展开形式： $$\mathbf{E} = \frac{1}{2}(\nabla_0 \mathbf{u} + (\nabla_0 \mathbf{u})^T + (\nabla_0 \mathbf{u})^T \nabla_0 \mathbf{u})$$ 在小变形下（$||\nabla_0 \mathbf{u}|| \ll 1$），忽略二次项： $$\mathbf{E} \approx \boldsymbol{\epsilon} = \frac{1}{2}(\nabla_0 \mathbf{u} + (\nabla_0 \mathbf{u})^T)$$ 这就是工程应变（无穷小应变）。

其他常用应变度量：

1. 左柯西-格林变形张量： $$\mathbf{B} = \mathbf{F}\mathbf{F}^T$$ 在当前构型中度量变形。

2. Hencky应变（对数应变）： $$\mathbf{E}_H = \frac{1}{2}\ln(\mathbf{C}) = \ln(\mathbf{U})$$ 其中$\mathbf{U}$来自极分解$\mathbf{F} = \mathbf{R}\mathbf{U}$。

优点：

• 对于单轴拉伸，$E_H = \ln(\lambda)$，具有可加性
• 压缩和拉伸对称：$E_H(-\lambda) = -E_H(\lambda)$

3. Almansi应变： $$\mathbf{e} = \frac{1}{2}(\mathbf{I} - \mathbf{B}^{-1})$$ 在当前构型中定义的应变度量。

主拉伸与主方向： 对$\mathbf{C}$进行特征值分解： $$\mathbf{C} = \sum_{i=1}^3 \lambda_i^2 \mathbf{N}_i \otimes \mathbf{N}_i$$ 其中：

• $\lambda_i$：主拉伸比（$\mathbf{F}$的奇异值）
• $\mathbf{N}_i$：拉格朗日主方向

体积与等容变形分解： 将变形梯度分解为体积和等容部分： $$\mathbf{F} = J^{1/3}\bar{\mathbf{F}}$$ 其中$\bar{\mathbf{F}} = J^{-1/3}\mathbf{F}$是等容变形梯度，满足$\det(\bar{\mathbf{F}}) = 1$。

相应的等容右柯西-格林张量： $$\bar{\mathbf{C}} = \bar{\mathbf{F}}^T\bar{\mathbf{F}} = J^{-2/3}\mathbf{C}$$

3.2.3 应力张量

应力描述了材料内部的力分布。在有限变形下，有多种应力度量，每种都有其物理意义和应用场景。

应力的物理定义： 考虑材料内部一个面元，其法向量为$\mathbf{n}$，作用在该面元上的力为$\mathbf{t}$（牵引力），则应力张量定义了它们之间的关系。

柯西应力张量（真实应力）： $$\mathbf{t} = \boldsymbol{\sigma} \mathbf{n}$$ 其中$\mathbf{t}$和$\mathbf{n}$都在当前构型中定义。柯西应力的物理意义：

• $\sigma_{ij}$：作用在法向为$\mathbf{e}_j$的面上，方向为$\mathbf{e}_i$的应力分量
• 对称张量：$\boldsymbol{\sigma} = \boldsymbol{\sigma}^T$（角动量守恒）
• 主应力：$\boldsymbol{\sigma}$的特征值

第一Piola-Kirchhoff应力（PK1，名义应力）： 将当前构型的力与参考构型的面积关联： $$\mathbf{t}_0 = \mathbf{P} \mathbf{N}$$ 其中$\mathbf{N}$是参考构型中的法向量。

与柯西应力的关系（Piola变换）： $$\mathbf{P} = J \boldsymbol{\sigma} \mathbf{F}^{-T}$$ 或反过来： $$\boldsymbol{\sigma} = \frac{1}{J} \mathbf{P} \mathbf{F}^T$$ 性质：

• 非对称张量（一般情况下）
• 功共轭于变形梯度：$\dot{W} = \mathbf{P} : \dot{\mathbf{F}}$
• 从应变能导出：$\mathbf{P} = \frac{\partial \psi}{\partial \mathbf{F}}$

第二Piola-Kirchhoff应力（PK2）： 完全在参考构型中定义的对称应力度量： $$\mathbf{S} = \mathbf{F}^{-1} \mathbf{P} = J \mathbf{F}^{-1} \boldsymbol{\sigma} \mathbf{F}^{-T}$$ 性质：

• 对称张量：$\mathbf{S} = \mathbf{S}^T$
• 功共轭于格林应变：$\dot{W} = \mathbf{S} : \dot{\mathbf{E}}$
• 从应变能导出：$\mathbf{S} = 2\frac{\partial \psi}{\partial \mathbf{C}} = \frac{\partial \psi}{\partial \mathbf{E}}$

Kirchhoff应力： $$\boldsymbol{\tau} = J \boldsymbol{\sigma}$$ 在不可压缩材料中特别有用，因为$J = 1$时$\boldsymbol{\tau} = \boldsymbol{\sigma}$。

应力张量之间的转换关系总结： $$\begin{aligned} \boldsymbol{\sigma} &= \frac{1}{J} \mathbf{P} \mathbf{F}^T = \frac{1}{J} \mathbf{F} \mathbf{S} \mathbf{F}^T \\ \mathbf{P} &= J \boldsymbol{\sigma} \mathbf{F}^{-T} = \mathbf{F} \mathbf{S} \\ \mathbf{S} &= J \mathbf{F}^{-1} \boldsymbol{\sigma} \mathbf{F}^{-T} = \mathbf{F}^{-1} \mathbf{P} \end{aligned}$$ 应力功率与共轭对： 不同的应力-应变对计算的功率相同： $$\dot{W} = \boldsymbol{\sigma} : \mathbf{D} = \mathbf{P} : \dot{\mathbf{F}} = \mathbf{S} : \dot{\mathbf{E}} = \frac{1}{2}\mathbf{S} : \dot{\mathbf{C}}$$ 其中$\mathbf{D} = \frac{1}{2}(\mathbf{L} + \mathbf{L}^T)$是变形率张量，$\mathbf{L} = \dot{\mathbf{F}}\mathbf{F}^{-1}$是速度梯度。

3.2.4 本构关系与材料参数

本构关系描述了应力和应变之间的关系，是材料模型的核心。不同的材料需要不同的本构模型，从简单的线弹性到复杂的粘弹塑性。

线弹性本构关系： 对于各向同性线弹性材料，在小应变下： $$\boldsymbol{\sigma} = \mathbf{C} : \boldsymbol{\epsilon}$$ 其中$\mathbf{C}$是四阶弹性张量。对于各向同性材料： $$\sigma_{ij} = \lambda \epsilon_{kk} \delta_{ij} + 2\mu \epsilon_{ij}$$ 或用矩阵形式（Voigt记号）： $$\begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{23} \\ \sigma_{13} \\ \sigma_{12} \end{bmatrix} = \begin{bmatrix} \lambda + 2\mu & \lambda & \lambda & 0 & 0 & 0 \\ \lambda & \lambda + 2\mu & \lambda & 0 & 0 & 0 \\ \lambda & \lambda & \lambda + 2\mu & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu \end{bmatrix} \begin{bmatrix} \epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{33} \\ 2\epsilon_{23} \\ 2\epsilon_{13} \\ 2\epsilon_{12} \end{bmatrix}$$ 材料参数及其物理意义：

1. 杨氏模量（Young's modulus）E： - 定义：单轴拉伸时应力与应变的比值 - 物理意义：材料的刚度 - 单位：Pa（帕斯卡） - 典型值：钢材 200 GPa，橡胶 0.01-0.1 GPa

2. 泊松比（Poisson's ratio）ν： - 定义：横向应变与纵向应变的负比值 - 物理意义：材料的体积变化倾向 - 范围：$\nu \in [0, 0.5)$（热力学稳定性要求） - 特殊值：

   • $\nu = 0$：软木（无横向变形）
   • $\nu \approx 0.5$：橡胶（近似不可压缩）
   • $\nu \approx 0.3$：大多数金属

3. Lamé参数： $$\lambda = \frac{E\nu}{(1+\nu)(1-2\nu)}, \quad \mu = \frac{E}{2(1+\nu)}$$

• $\lambda$：Lamé第一参数（无直接物理意义）
• $\mu$：剪切模量，也记作$G$

4. 体积模量（Bulk modulus）K： $$K = \frac{E}{3(1-2\nu)} = \lambda + \frac{2\mu}{3}$$

• 物理意义：抵抗体积变化的能力
• 当$\nu \to 0.5$时，$K \to \infty$（不可压缩）

参数之间的转换关系： 给定任意两个独立参数，可以计算其他参数： $$\begin{aligned} E &= \frac{9KG}{3K + G} = \frac{\mu(3\lambda + 2\mu)}{\lambda + \mu} \\ \nu &= \frac{3K - 2G}{2(3K + G)} = \frac{\lambda}{2(\lambda + \mu)} \\ K &= \frac{\lambda + 2\mu}{3} = \frac{E}{3(1-2\nu)} \\ G &= \mu = \frac{E}{2(1+\nu)} = \frac{3K(1-2\nu)}{2(1+\nu)} \end{aligned}$$ 材料稳定性要求： 热力学第二定律要求应变能为正定，这导致：

• $E > 0$（拉伸刚度为正）
• $G > 0$（剪切刚度为正）
• $K > 0$（体积刚度为正）
• $-1 < \nu < 0.5$（实际材料通常$\nu > 0$）

各向异性材料： 对于正交各向异性材料，需要9个独立参数：

• 3个杨氏模量：$E_1, E_2, E_3$
• 3个剪切模量：$G_{12}, G_{23}, G_{31}$
• 3个泊松比：$\nu_{12}, \nu_{23}, \nu_{31}$

对于一般各向异性材料，弹性张量$\mathbf{C}$有21个独立分量。

3.3 超弹性材料模型

3.3.1 应变能密度函数

超弹性材料通过应变能密度函数$\psi(\mathbf{F})$完全描述其力学行为： $$\mathbf{P} = \frac{\partial \psi}{\partial \mathbf{F}}$$ 应变能密度函数必须满足：

• 客观性：$\psi(\mathbf{QF}) = \psi(\mathbf{F})$ 对所有旋转矩阵$\mathbf{Q}$
• 材料对称性：反映材料的各向同性或各向异性
• 凸性：保证材料稳定性

3.3.2 Neo-Hookean模型

Neo-Hookean模型是最简单的超弹性模型，适用于橡胶类材料： $$\psi(\mathbf{F}) = \frac{\mu}{2}(\text{tr}(\mathbf{F}^T\mathbf{F}) - 3) - \mu \ln(J) + \frac{\lambda}{2}\ln^2(J)$$ 对应的第一Piola-Kirchhoff应力： $$\mathbf{P}(\mathbf{F}) = \mu(\mathbf{F} - \mathbf{F}^{-T}) + \lambda \ln(J) \mathbf{F}^{-T}$$ Neo-Hookean模型的特点：

• 在小变形下退化为线弹性
• 可以处理大变形
• 数学形式简单，计算效率高

3.3.3 Corotated模型

Corotated模型通过极分解处理大旋转： $$\mathbf{F} = \mathbf{R}\mathbf{S}$$ 其中$\mathbf{R}$是旋转部分，$\mathbf{S}$是拉伸部分。

应变能密度函数： $$\psi(\mathbf{F}) = \mu \sum_i (\sigma_i - 1)^2 + \frac{\lambda}{2}(J - 1)^2$$ 其中$\sigma_i$是$\mathbf{F}$的奇异值。

第一Piola-Kirchhoff应力： $$\mathbf{P}(\mathbf{F}) = 2\mu(\mathbf{F} - \mathbf{R}) + \lambda(J - 1)J\mathbf{F}^{-T}$$

3.3.4 Mooney-Rivlin模型

Mooney-Rivlin模型用于更精确地描述橡胶材料： $$\psi = C_{10}(I_1 - 3) + C_{01}(I_2 - 3) + \frac{1}{D_1}(J - 1)^2$$ 其中：

• $I_1 = \text{tr}(\mathbf{C})$：第一不变量
• $I_2 = \frac{1}{2}[(\text{tr}(\mathbf{C}))^2 - \text{tr}(\mathbf{C}^2)]$：第二不变量
• $C_{10}, C_{01}, D_1$：材料参数

3.4 基于六面体网格的拉格朗日有限元

3.4.1 等参单元与形函数

六面体单元使用三线性形函数，在参考坐标$(\xi, \eta, \zeta) \in [-1,1]^3$中： $$N_i(\xi, \eta, \zeta) = \frac{1}{8}(1 + \xi_i\xi)(1 + \eta_i\eta)(1 + \zeta_i\zeta)$$ 其中$(\xi_i, \eta_i, \zeta_i)$是第$i$个节点的参考坐标。

几何映射（等参映射）： $$\mathbf{x} = \sum_{i=1}^8 N_i(\xi, \eta, \zeta) \mathbf{x}_i$$

3.4.2 数值积分

使用高斯积分计算单元矩阵： $$\mathbf{K}_e = \int_{\Omega_e} \mathbf{B}^T \mathbf{D} \mathbf{B} \, d\Omega \approx \sum_{g=1}^{n_g} w_g \mathbf{B}^T(\xi_g) \mathbf{D} \mathbf{B}(\xi_g) |\mathbf{J}(\xi_g)|$$ 其中：

• $n_g$：高斯点数（通常$2×2×2$）
• $w_g$：高斯权重
• $\mathbf{J}$：雅可比矩阵

3.4.3 单元刚度矩阵组装

每个六面体单元贡献一个$24×24$的单元刚度矩阵（3D情况，每节点3个自由度）： $$\mathbf{K} = \bigcup_{e=1}^{n_e} \mathbf{K}_e$$ 组装过程需要：

1. 计算每个单元的刚度矩阵
2. 根据节点编号映射到全局矩阵
3. 处理边界条件

3.4.4 边界条件处理

Dirichlet边界条件（位移约束）：

• 直接设置：$u_i = \bar{u}_i$
• 罚函数法：添加大刚度$\alpha(u_i - \bar{u}_i)^2$
• 拉格朗日乘子法：引入约束力

Neumann边界条件（力边界）： $$\mathbf{f}_i = \int_{\Gamma_N} N_i \mathbf{t} \, d\Gamma$$ 其中$\mathbf{t}$是表面力。

3.5 基于四面体网格的拉格朗日有限元

3.5.1 线性四面体单元

对于线性四面体，变形梯度在单元内为常数： $$\mathbf{F} = \mathbf{D}_s \mathbf{D}_m^{-1}$$ 其中：

• $\mathbf{D}_s = [\mathbf{x}_1 - \mathbf{x}_4, \mathbf{x}_2 - \mathbf{x}_4, \mathbf{x}_3 - \mathbf{x}_4]$：变形边矩阵
• $\mathbf{D}_m = [\mathbf{X}_1 - \mathbf{X}_4, \mathbf{X}_2 - \mathbf{X}_4, \mathbf{X}_3 - \mathbf{X}_4]$：参考边矩阵

3.5.2 常应变假设

线性四面体单元假设应变在单元内均匀分布：

• 优点：实现简单，适合大变形
• 缺点：需要更多单元才能准确捕捉弯曲变形
• 可能过于刚硬（过度约束）

3.5.3 体积锁定与缓解策略

当材料接近不可压缩（$\nu \to 0.5$）时，出现体积锁定现象：

原因：体积模量$K \to \infty$，数值误差被放大

缓解方法：

1. 混合有限元：分别插值位移和压力
2. 降阶积分：使用较少的积分点
3. F-bar方法：修正变形梯度的体积部分
4. 增强应变元：添加额外的应变模式

3.5.4 网格生成与质量控制

高质量的四面体网格需要满足：

• 最小角度：避免过于尖锐的角度（建议>15°）
• 纵横比：控制最长边与最短边的比值（建议<10）
• 体积：避免退化的单元（体积接近0）

常用网格生成方法：

• Delaunay三角剖分
• 前沿推进法
• 八叉树细分

3.6 可逆(Invertible)有限元

3.6.1 单元翻转问题

当$\det(\mathbf{F}) \leq 0$时，单元发生翻转：

• 物理上不合理（材料穿透自身）
• 数值上导致发散（应力无穷大）
• 仿真崩溃

翻转通常发生在：

• 大变形情况
• 碰撞处理
• 不当的时间步长

3.6.2 SVD分解与修复

使用奇异值分解(SVD)检测和修复翻转： $$\mathbf{F} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^T$$ 修复步骤：

1. 计算SVD
2. 检查$\det(\mathbf{U})$和$\det(\mathbf{V})$的符号
3. 如果为负，翻转相应的列
4. 确保所有奇异值为正

3.6.3 能量势垒方法

在应变能中添加势垒项防止翻转： $$\psi_{barrier}(J) = \begin{cases} 0 & J > \epsilon \\ \alpha(J - \epsilon)^2 & J \leq \epsilon \end{cases}$$ 其中$\epsilon$是小的正数（如0.01）。

3.6.4 稳定性保证

确保数值稳定的技术：

1. 投影方法：将$\mathbf{F}$投影到可行域
2. 线搜索：在Newton迭代中使用线搜索
3. 自适应时间步：检测到接近翻转时减小时间步
4. 隐式积分：提高大变形下的稳定性

3.7 拓扑优化

3.7.1 最小柔度问题

拓扑优化的标准形式： $$\begin{aligned} \min_{\boldsymbol{\rho}} \quad & c(\boldsymbol{\rho}) = \mathbf{u}^T \mathbf{K}(\boldsymbol{\rho}) \mathbf{u} \\ \text{s.t.} \quad & \mathbf{K}(\boldsymbol{\rho})\mathbf{u} = \mathbf{f} \\ & \sum_{e=1}^{n_e} v_e \rho_e \leq V_{max} \\ & 0 < \rho_{min} \leq \rho_e \leq 1 \end{aligned}$$ 其中：

• $\rho_e$：单元$e$的密度
• $c$：结构柔度（compliance）
• $V_{max}$：材料体积约束

3.7.2 SIMP方法

固体各向同性材料惩罚(SIMP)方法： $$E_e(\rho_e) = \rho_e^p E_0$$ 其中：

• $p$：惩罚参数（通常$p=3$）
• $E_0$：实体材料的杨氏模量

惩罚中间密度，推动设计向0-1分布。

3.7.3 敏感度分析

目标函数对设计变量的导数： $$\frac{\partial c}{\partial \rho_e} = -p\rho_e^{p-1} \mathbf{u}_e^T \mathbf{k}_0 \mathbf{u}_e$$ 其中$\mathbf{k}_0$是密度为1时的单元刚度矩阵。

敏感度滤波（避免棋盘格）： $$\frac{\widetilde{\partial c}}{\partial \rho_e} = \frac{\sum_{i \in N_e} H_{ei} \rho_i \frac{\partial c}{\partial \rho_i}}{\rho_e \sum_{i \in N_e} H_{ei}}$$

3.7.4 优化准则法(OC)

基于最优性条件的更新公式： $$\rho_e^{new} = \begin{cases} \max(\rho_{min}, \rho_e - m) & \text{if } \rho_e B_e^{\eta} \leq \max(\rho_{min}, \rho_e - m) \\ \min(1, \rho_e + m) & \text{if } \rho_e B_e^{\eta} \geq \min(1, \rho_e + m) \\ \rho_e B_e^{\eta} & \text{otherwise} \end{cases}$$ 其中：

• $B_e = -\frac{\partial c/\partial \rho_e}{\lambda \frac{\partial V/\partial \rho_e}}$
• $\lambda$：拉格朗日乘子（通过二分法确定）
• $m$：移动限制（如0.2）
• $\eta$：阻尼系数（如0.5）

3.8 高级Taichi特性（1）

3.8.1 可微编程基础

Taichi支持自动微分，使物理仿真可微：

@ti.data_oriented
class ElasticBody:
    def __init__(self):
        self.x = ti.Vector.field(3, dtype=ti.f32, shape=n_verts, needs_grad=True)
        self.v = ti.Vector.field(3, dtype=ti.f32, shape=n_verts)
        self.f = ti.Vector.field(3, dtype=ti.f32, shape=n_verts)
        self.U = ti.field(dtype=ti.f32, shape=(), needs_grad=True)  # 势能

3.8.2 自动微分

Taichi提供前向和反向自动微分：

with ti.Tape(loss):
    forward_simulation()
# 梯度自动计算在 x.grad 中

关键特性：

• 支持复杂控制流
• 高效的梯度累积
• 与并行计算兼容

3.8.3 梯度计算与优化

利用自动微分进行参数优化：

@ti.kernel
def compute_loss():
    # 计算目标函数
    loss[None] = 0.0
    for i in range(n_verts):
        loss[None] += 0.5 * (x[i] - target[i]).norm_sqr()

# 优化循环
for iter in range(max_iters):
    with ti.Tape(loss):
        simulate()
        compute_loss()

    # 更新参数
    lr = 0.01
    for i in range(n_params):
        params[i] -= lr * params.grad[i]

3.8.4 物理过程的端到端优化

自动微分使得端到端优化成为可能：

应用场景：

1. 材料参数识别：从观测数据反推材料属性
2. 形状优化：优化结构形状以满足性能要求
3. 控制优化：设计最优控制序列
4. 逆向设计：从目标性能设计材料分布

优势：

• 无需手动推导梯度
• 支持复杂的物理模型
• 与神经网络无缝集成

本章小结

本章深入介绍了有限元方法的理论基础和实现技术：

核心概念：

• 弱形式将PDE转化为积分方程，降低对解的光滑性要求
• 变形梯度$\mathbf{F}$描述了材料的局部变形
• 超弹性材料通过应变能密度函数定义本构关系
• 不同网格类型（六面体、四面体）各有优劣

关键公式：

• 弱形式：$\int_\Omega \nabla w \cdot \nabla u \, d\Omega = \int_\Omega wf \, d\Omega$
• Neo-Hookean应力：$\mathbf{P} = \mu(\mathbf{F} - \mathbf{F}^{-T}) + \lambda \ln(J) \mathbf{F}^{-T}$
• SIMP插值：$E(\rho) = \rho^p E_0$

实践要点：

• 选择合适的材料模型和单元类型
• 处理数值不稳定性（翻转、锁定）
• 利用自动微分进行优化设计

练习题

基础题

习题3.1 推导一维杆的弱形式。给定强形式： $$-\frac{d}{dx}\left(EA\frac{du}{dx}\right) = f(x), \quad 0 < x < L$$ 其中$E$是杨氏模量，$A$是横截面积，$f(x)$是分布载荷。

习题3.2 给定二维变形梯度： $$\mathbf{F} = \begin{bmatrix} 1.2 & 0.3 \\ 0.1 & 0.8 \end{bmatrix}$$ 计算：(a) 体积比$J$ (b) 右柯西-格林张量$\mathbf{C}$ (c) 格林应变$\mathbf{E}$

习题3.3 对于Neo-Hookean材料，证明在小变形极限下（$\mathbf{F} = \mathbf{I} + \boldsymbol{\epsilon}$，$||\boldsymbol{\epsilon}|| \ll 1$），第一Piola-Kirchhoff应力退化为线弹性。

挑战题

习题3.4 设计一个算法检测四面体单元是否接近翻转，并提出一个修复策略。考虑数值稳定性和计算效率。

习题3.5 推导SIMP方法中体积约束的拉格朗日乘子$\lambda$的更新公式。使用KKT条件。

习题3.6 考虑一个2D悬臂梁拓扑优化问题。梁长$L=2$，高$H=1$，左端固定，右端中点受向下的集中力$F=1$。材料体积限制为30%。设计一个完整的优化流程，包括敏感度滤波和收敛准则。

习题3.7 实现一个简单的可微有限元求解器，用于优化悬臂梁的材料分布以最小化末端位移。使用自动微分计算梯度。

for iter in range(100):
    with ti.Tape(loss):
        assemble_K()
        solve_system()
        compute_loss()

    # 梯度下降
    for i, j in density:
        density[i,j] -= lr * density.grad[i,j]
        density[i,j] = clamp(density[i,j], 0.1, 1.0)

习题3.8 分析并比较六面体单元和四面体单元在以下情况下的性能：(a) 弯曲主导的变形 (b) 近似不可压缩材料 (c) 大旋转变形。讨论各自的优缺点和适用场景。

常见陷阱与错误 (Gotchas)

1. 单元翻转：大变形时忘记检查$\det(\mathbf{F}) > 0$，导致仿真爆炸
2. 体积锁定：对近似不可压缩材料使用完全积分，导致过度刚硬
3. 数值积分阶数：积分阶数过低导致秩缺陷，过高导致锁定
4. 材料参数：泊松比接近0.5时数值不稳定，需要特殊处理
5. 网格质量：忽视网格质量检查，退化单元导致矩阵奇异
6. 边界条件：Dirichlet和Neumann边界条件混淆，导致错误结果
7. 应力度量：混淆不同应力张量（PK1、PK2、Cauchy），导致错误的本构实现
8. 优化收敛：拓扑优化中忽略敏感度滤波，出现棋盘格现象

最佳实践检查清单

有限元实现

• [ ] 验证单个单元的刚度矩阵对称性
• [ ] 检查组装后的全局矩阵是否满足预期的性质（对称、正定）
• [ ] 实现patch test验证收敛性
• [ ] 监控能量守恒（对于保守系统）
• [ ] 检查单元质量指标

材料模型

• [ ] 验证小变形极限下退化为线弹性
• [ ] 检查应力-应变曲线的物理合理性
• [ ] 测试各向同性材料的旋转不变性
• [ ] 验证能量函数的凸性

数值稳定性

• [ ] 实现单元翻转检测和修复
• [ ] 对不可压缩材料使用适当的公式
• [ ] 选择合适的时间积分方案
• [ ] 实现自适应时间步长控制

优化设计

• [ ] 使用体积约束确保可制造性
• [ ] 应用敏感度滤波避免数值不稳定
• [ ] 设置合理的收敛准则
• [ ] 验证梯度计算的正确性（有限差分检验）

性能优化

• [ ] 利用单元矩阵的对称性减少计算
• [ ] 使用稀疏矩阵存储和求解器
• [ ] 实现并行化（单元级并行）
• [ ] 考虑使用自适应网格细化