跳转至

第六章 有限元 (Finite Elements)

§1 作者

本章是Ch 1—5理论与实验的"数值实现"。Humphrey的策略："3种弱形式 (Galerkin/虚功/最小势能) → 1个4-1 Lagrange 单元 → 1个轴对称算例 (axisymmetric) → 1个膜膨胀算例"。全章是Ch 7、8、10所有有限元计算的"母版"。Humphrey明确说"商业有限元软件 (Abaqus, ANSYS) 不适合软组织 (生物力学市场份额小)，所以大多数重要贡献仍由自定义代码 完成"。

§2 内容概述

分七节： - 6.1 基本方程：参考形式 (Eq. 6.2—6.7) + 3种弱形式：加权残差 (Eq. 6.8—6.9)、虚功 (Eq. 6.10—6.15)、最小势能 (Eq. 6.16—6.18)。 - 6.2 插值、积分、求解器：Hermite/Lagrange 插值 (Eq. 6.20—6.34) + isoparametric 映射 (Eq. 6.36—6.39) + Gauss 求积 (Eq. 6.40—6.43) + 边界条件装配 (Eq. 6.44—6.45) + 矩阵分解 (Eq. 6.46—6.48) + Newton-Raphson (Eq. 6.49—6.51)。 - 6.3 一个示范性公式：不可压缩 4-1 单元 (4 nodes for u + 1 node for p)，完整推导 \(\Delta\Pi\) 的矩阵组装 (Eq. 6.79—6.93)。 - 6.4 膜膨胀：Fried 1982 经典解 (轴对称 Mooney-Rivlin 膜) — Ch 7血管壁/动脉瘤膜的"原始版"。 - 6.5 反有限元 (Inverse FE)：用测量的位移反推本构参数。Ch 5、7数据反演的桥梁。 - 6.6 习题 (15题)。 - 6.7 参考文献 (38条)。

§3 核心方程与概念

3.1 强形式 → 弱形式 (Eq. 6.2—6.18)

(1) 参考形式平衡 (Eq. 6.2, 6.6) $\(\nabla_0\cdot\mathbf{P} + \rho_0\mathbf{b} = \mathbf{0} \quad \text{in } \Omega_0\)$ 其中 \(\mathbf{P} = \mathbf{S}\cdot\mathbf{F}^T\) (Ch 3.43)，\(\mathbf{S} = 2\partial W/\partial\mathbf{C} - p\mathbf{C}^{-1}\) (Ch 3.74)。

(2) 加权残差法 (Eq. 6.8—6.9) 设试函数 \(\mathbf{v}(\mathbf{X})\) 满足位移边界条件 (kinematically admissible)： $\(\int_{\Omega_0} \mathbf{R}_\Omega(\mathbf{v})\cdot\mathbf{w}\, dV + \int_{\partial\Omega_0} \mathbf{R}_{\partial\Omega}(\mathbf{v})\cdot\mathbf{w}\, dA = 0\)$ 其中 \(\mathbf{R}_\Omega = \nabla_0\cdot(\mathbf{S}\cdot\mathbf{F}^T) + \rho_0\mathbf{b}\)。Galerkin法 = \(\mathbf{w} = \delta\mathbf{v}\) 是试函数的导数 → 系数矩阵对称 (Ch 6.1.1末段)。

(3) 虚功原理 (Eq. 6.10—6.15) Holzapfel 2000 形式： $\(\int_{\partial(\Omega_0)} \mathbf{T}^{(N)}\cdot\delta\mathbf{u}\, dA = \int_{\partial\Omega_{0T}} \bar{\mathbf{T}}^{(N)}\cdot\delta\mathbf{u}\, dA\)$ 推导得： $\(\int_{\Omega_0} \mathbf{P}:\nabla_0(\delta\mathbf{u})\, dV = \int_{\Omega_0} \rho_0\mathbf{b}\cdot\delta\mathbf{u}\, dV + \int_{\partial\Omega_{0T}} \bar{\mathbf{T}}^{(N)}\cdot\delta\mathbf{u}\, dA\)$ 关键 (Eq. 6.15)：\(\int\mathbf{P}:\nabla_0(\delta\mathbf{u})dV = \int\mathbf{S}:\delta\mathbf{E}\, dV\)（等价于势能形式）。

(4) 最小总势能 (Eq. 6.16—6.18) $\(\Pi = \int_{\Omega_0} W(\mathbf{E})\, dV - \int_{\Omega_0} \rho_0\mathbf{b}\cdot\mathbf{u}\, dV - \int_{\partial\Omega_{0T}} \bar{\mathbf{T}}^{(N)}\cdot\mathbf{u}\, dA\)$ \(\delta\Pi = 0\) 给出平衡。仅对超弹 (保守) 适用。与虚功相比，"虚功更一般，不要求材料守恒" (6.1.3末段)。

3.2 插值 (Eq. 6.20—6.39)

(5) Hermite 插值 (Eq. 6.20—6.22) \(v(\xi) = \sum_{j=1}^{n}\sum_{s=1}^{N+1} {}^{(N)}H_j^s(\xi) \frac{d^s v_j}{d\xi^s}\)，关键性质 \({}^{(N)}H_j^s(\xi_k) = \delta_{jk}\delta_{sr}\)。Ch 10 双三次 Hermite 是 Hunter-McCulloch 1998 心肌有限元的标准选择。

(6) Lagrange 插值 (Eq. 6.23—6.27) \(N=0\) 阶 Hermite = Lagrange： $\(L^j(\xi) = \prod_{k\neq j} \frac{\xi - \xi_k}{\xi_j - \xi_k}, \quad v(\xi) = \sum_{j=1}^{n} L^j(\xi) v_j\)$ Ch 5 marker tracking (Eq. 5.9) 就是 Lagrange 插值在 2D 的应用。

(7) Isoparametric 映射 (Eq. 6.36) \(X = \sum L^j(\xi) X^j\)，\(u = \sum L^j(\xi) u^j\)（位移和坐标用相同插值）。作用：把任意几何形状映射到 \(\xi \in [-1, 1]\) 标准化区间 → Gauss 积分可重用。

(8) Gauss 求积 (Eq. 6.40—6.43) \(n\) 阶 Gauss 求积 = 精确积分 \(2n-1\) 阶多项式。Table 6.1 给出 \(m=1\) 到 \(m=4\) 的 Gauss 点和权重。m=2 即可精确积分 3 阶多项式。Eq. 6.42：2 个 Gauss 点 \(\pm 1/\sqrt{3}\) 积分 1D。

3.3 求解器 (Eq. 6.44—6.51)

(9) 边界条件装配 (Eq. 6.44—6.45) Dirichlet (\(u_2 = \bar{u}_2\))：置 \(K_{21} = K_{23} = 0, K_{22} = 1, Q_2 = \bar{u}_2\)。保留对称性的"分块"形式 (Eq. 6.45)。

(10) 矩阵分解 (Eq. 6.46—6.48) LDL 分解：\([K] = [L][D][L]^T\)。前向替换求 \(\{g\} = [L]^{-1}\{Q\}\)，后向替换求 \(\{q\}\)。比 Gauss 消去数值稳定。

(11) Newton-Raphson (Eq. 6.49—6.51) \(q^{i+1} = q^i - (\partial g/\partial q)^{-1} g(q^i)\)。\(\partial g/\partial q\) = 切线刚度。Full NR 每步更新；Modified NR 不更新。非线性越强 → 更新越频繁。

3.4 不可压缩有限元 (Eq. 6.52—6.93)

(12) 罚方法 (Penalty, Eq. 6.52) $\(W = W(\mathbf{C}) - \frac{1}{2}k(III_C - 1)^2\)$ \(k \to \infty\) 强制 \(III_C = 1\)。优点：不需要额外的 \(p\) 自由度；缺点：病态刚度矩阵。

(13) 4-1 单元 (Eq. 6.72—6.79) 4 节点双线性插值 (位移) + 1 节点常数 Lagrange 乘子 \(p\)。9 个自由度/单元。问题：constant \(p\) 在某些边界条件下出现 "checkerboarding" 棋盘式伪应力模式。解决：Brezzi-Babuska inf-sup 条件 (Glowinski-LeTallec 1982)。

(14) 增量公式 (Eq. 6.54—6.64) \(\Delta\Pi = \Pi^{t+\Delta t} - \Pi^t\) + Taylor展开 \(W + \frac{1}{2}\partial^2 W/\partial E^2 \cdot \Delta E^2\)。Eq. 6.64给出 \(\Delta\Pi\) 的线性 + 二次项。

(15) 轴对称有限元 (Eq. 6.66—6.93) \(r = R + u_R(R,Z), \theta = \Theta, z = Z + u_Z(R,Z)\)。\(\mathbf{C}\) 和 \(\mathbf{E}\) 的分量在柱坐标下 (Eq. 6.68—6.69)。Eq. 6.79 给 \(\Delta\varepsilon = [G]\{\Delta u^A\}\) 矩阵形式。最终 (Eq. 6.93)： $\(\begin{pmatrix} [K] & [L] \\ [L]^T & 0 \end{pmatrix} \begin{Bmatrix} \Delta q \\ \Delta p \end{Bmatrix} = \begin{Bmatrix} Q \\ -\frac{1}{2}\int(III_C - 1)dV \end{Bmatrix}\)$ 9x9 单元方程，半正定（\(p\) 对应的零行）→ 三重分解 (triple factoring)。

3.5 膜膨胀 (6.4节)

(16) Fried 1982 膜膨胀 轴对称 + Mooney-Rivlin → 数值解。Ch 7血管膨胀/动脉瘤膨胀是这个公式的扩展。

3.6 反有限元 (6.5节)

(17) Inverse FE 已知测量位移 \(u^{exp}(\mathbf{X})\)，反推本构参数。算法：最小化 \(\|u^{FE}(\mathbf{X}) - u^{exp}(\mathbf{X})\|^2\)。Ch 7血管/动脉瘤中实测 marker 跟踪 → 反演材料参数的标准方法。

§4 关键算法或建模方法

  1. 3种弱形式等价性 (Galerkin = 虚功 = 最小势能)：在超弹条件下三者等价于同一代数系统。Ch 6 核心trick。
  2. Newton-Raphson 增量求解：每步 \(\Delta q\) 由 \([\partial g/\partial q] \Delta q = -g\) 求。Ch 7、8、10 所有非线性问题的标准。
  3. 不可压缩约束的3种实现 (Lagrange / Penalty / Augmented Lagrangian)：各有所长。Ch 7血管用 Lagrange + 4-1 单元。
  4. Isoparametric 映射 (Eq. 6.36)：把任意几何标准化到 \(\xi \in [-1, 1]\) → Gauss 积分可重用 → 编程简单。
  5. 4-1 单元 + Brezzi-Babuska inf-sup 检验：确保 \(p\) 离散化合理，避免伪模式。Ch 7血管壁的"沿径向 8 单元"是这套。
  6. Cauchy应力后处理 (Eq. 6.71)：先在 \(\kappa_0\) 算 \(\mathbf{S}\)，再用 \(t_{ij} = \det F \cdot F_{iA} F_{jB} S_{AB}\) 转到 \(\kappa_t\)。
  7. Inverse FE (6.5节)：测 \(\to\) 算材料参数的反问题。是 Ch 5 实验数据接口 → Ch 7 血管本构反演的桥梁。

§5 关键结论

  1. 有限元 = 强形式 (ODE/PDE) → 弱形式 (积分) → 离散代数三步。6.1节末段总结：从总势能出发最自然。
  2. 3种弱形式在超弹条件下等价：虚功 + Galerkin + 最小势能 → 同一 \([\mathbf{K}]\{\mathbf{q}\} = \{\mathbf{Q}\}\)。Ch 7 所有血管有限元都用此。
  3. 不可压缩约束 = 病态刚度：Lagrange 乘子 + 4-1 单元半正定；Penalty 病态但简单；Augmented Lagrangian 折中。Ch 7 血管实际用 augmented Lagrangian (Srinivasan-Perucchio 1994)。
  4. 位移增量法 (Eq. 6.64) 是大变形问题的标准：每步小变形近似 → 累积。Ch 6.3.3轴对称是 Ch 7 血管管壁、Ch 8 动脉瘤、Ch 10 心肌的"母公式"。
  5. 4-1 单元"checkerboarding"是经典陷阱：常数 \(p\) 离散化常出问题 → 必须用 inf-sup 检验。Ch 7、10 论文常见此陷阱。
  6. 后处理求 Cauchy 应力 (Eq. 6.71) 而不是直接 \(\sigma = J^{-1}\mathbf{F}\mathbf{S}\mathbf{F}^T\)：Humphrey强调先算 \(\mathbf{S}\)，再转到 \(\mathbf{t}\) —— 因为 \(\mathbf{S}\) 定义在 \(\kappa_0\)（不变），更易算。
  7. 商业软件 (Abaqus, ANSYS) 不适合软组织：6.0节末段警告 → "自定义代码仍是主流"。Ch 7、10 实际论文多用 ABAQUS + UMAT/HYPERELAS 用户子程序。

§6 挑战和开放性问题

  1. 不可压缩的"常数 \(p\) 棋盘模式"：Babuska inf-sup 检验提供了判据，但实际网格设计时仍常出伪模式——Ch 7 血管 3 层结构 + 每层不同本构时尤其难。
  2. Large strain 增量法的"漂移"问题：每步假设小应变 → 大应变累积 100+ 步后可能偏离真解。线搜索 + 弧长法 (Riks) 是 Ch 8 动脉瘤破裂分析的关键，但 Humphrey 6.2.5 节只提 Newton-Raphson，未讲弧长法。
  3. 膜 vs 实体单元的选择：6.4节膜（Fried 1982）只用于"无限薄"血管壁；Ch 7 实际血管 \(h/R \sim 0.1\) → 实体单元更安全。6.4节没讲"何时膜假设失效"。
  4. Inverse FE 的不适定性：已知 \(\mathbf{u}^{exp}\) 反推 \(\mathbf{W}\) 是ill-posed (多个 \(W\) 可能给出相同 \(\mathbf{u}\))。6.5节没给正则化策略。
  5. 纤维增强本构的有限元实现：Ch 7、10 血管/心肌是横观各向同性 (4 个不变量 \(I_C, II_C, IV_C, V_C\))，4-1 单元直接套用 (Eq. 6.91) 时 \(S_{AB}\) 需重新推导。Humphrey 把这留给读者 (Eq. 6.93末段)。
  6. 本构梯度 (constitutive gradient) 有限元：\(\mathbf{K} = \partial^2 W / \partial \mathbf{C}^2\) 在大变形下可能不正定——可能导致 Newton-Raphson 步"不下降"。Ch 6 没说这种情况怎么办。
  7. 商业软组织有限元软件的现状：6.0节说"市场小"——但 2026 年实际看，Ansys、FEBio、ABAQUS 都有专门的生物力学包。Humphrey 2002 年的预言并未完全应验。

§7 个人反思与批判性分析

Humphrey的有限元章的"轴对称 + 4-1 单元 + 增量法"组合，是心血管力学研究的"事实标准" —— Ch 7 血管、Ch 8 动脉瘤、Ch 10 心肌的所有数值工作都从这套出发。但 Humphrey 也明确说"商业软件不适用"——这一论断在 2026 年已部分失效：FEBio, ABAQUS + UMAT, ANSYS Mechanical + 用户子程序已是工业/学术标准。

值得自己复现的推导： - Eq. 6.64 (增量势能)：从 \(\Pi\) Taylor展开到 \(\Delta\Pi\) 线性+二次项 → 写出 \(\delta\Delta\Pi = 0\) 的方程 (Eq. 6.91)。 - Eq. 6.79 (\(\Delta\varepsilon = [G]\{\Delta u^A\}\))：4 节点双线性形函数的 Jacobian 链式法则 → 这步是"几何非线性"的核心。 - Eq. 6.93 (单元刚度矩阵组装)：把 \(\Delta\Pi\) 中所有项写成 \(\{q\}^T [K] \{q\} + \{Q\}^T \{q\}\) 形式。 - Gauss 积分的实现：Table 6.1 的 \(m=2\) 即可。 - Newton-Raphson 收敛准则：相对残差 \(< 10^{-6}\)，相对位移增量 \(< 10^{-4}\)。

与作者的潜在对话： - 6.0节末段"商业软件不适用"——但 2026 年 FEBio 已成生物力学事实标准 (Maas, Ateshian 等)。您对 2026 年的研究者推荐FEniCS / FEBio / ABAQUS + UMAT 中的哪一套？ - 6.2.5节只提 Newton-Raphson。Ch 8 动脉瘤破裂问题常需弧长法 (Riks 1979) 处理 limit-point instability。为什么不加一小节"弧长法"？ - 6.3.1节提到"4-1 单元不能在纯位移边界条件下使用" (Oden-Carey 1984) — 但 Ch 7 血管都是纯位移边界！Ch 7 是否改用 9-3 单元？ - 6.5节反有限元只 1 页——Ch 7 血管本构反演的核心工具，建议加 5-10 页（含 Levenberg-Marquardt + Sobolev 正则化）。

总结性批评： - 6.3节 4-1 单元推导占了 ~25 页，但 4-1 单元实际上已过时——现代 Ch 7 血管有限元多用 Q1Q0 (B-bar) 或 Q1P0。6.3节应至少给 Q1Q0 引用。 - 6.4节膜膨胀 用了 Fried 1982 经典解——但 Ch 7 实际血管是3层 (intima, media, adventitia) 复合结构，每层不同本构。6.4节没给层合板公式。 - 6.5节反有限元 1 页——Humphrey自己实验室的Inverse FE 工作 (Humphrey, Yin 1986 等) 显然涉及更多内容，但没在6.5节展开。 - 6.6节习题 15题多数是代数推导 (Eq. 6.79的具体矩阵)，缺"实际血管问题的有限元建模"练习。 - 6.0节末段"自定义代码仍是主流"——但 2026 年看，FEBio 已被广泛使用 (生物力学专用的开源 FEA)。应更新对软件生态的判断。

§8 重要参考文献

  • [X1] Oden JT (1972) Finite Elements of Nonlinear Continua. McGraw-Hill. — Ch 6 最主要追随的"圣典"，非线性有限元理论源头。
  • [X2] Bathe K-J (1996) Finite Element Procedures. Prentice-Hall. — Newton-Raphson、弧长法、不可压缩处理的现代标准参考。
  • [X3] Cook RD, Malkus DS, Plesha ME (1989) Concepts and Applications of Finite Element Analysis. Wiley. — 工科有限元入门标准教材。
  • [X4] Zienkiewicz OC, Taylor RL (1991) The Finite Element Method. McGraw-Hill. — 4卷本有限元百科。
  • [X5] Bonet J, Wood RD (1997) Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge. — 几何非线性 + 有限元结合的现代教材。
  • [X6] Crisfield MA (1997) Non-Linear Finite Element Analysis of Solids and Structures. Wiley. — 弧长法等高级非线性技术。
  • [X7] Holzapfel GA (2000) Nonlinear Solid Mechanics: A Continuum Approach for Engineering. Wiley. — 6.1.2节虚功形式 (Eq. 6.10) 来源。
  • [X8] Tong P, Rossettos JN (1977) Finite-Element Method: Basic Technique and Implementation. MIT Press. — 6.2节 Hermite/Lagrange 插值的标准参考。
  • [X9] Hunter PJ, McCulloch AD, ter Keurs HEDJ (1997) Modelling the mechanical properties of cardiac muscle. Prog Biophys Mol Biol 69:289-331. — 双三次 Hermite 在心脏有限元中的应用。
  • [X10] Glowinski R, LeTallec P (1982, 1984) Numerical Analysis of Variational Inequalities. North-Holland. — Lagrange 乘子 + Brezzi-Babuska 条件。
  • [X11] Malkus DS, Hughes TJR (1978) Mixed finite element methods - reduced and selective integration techniques. Comp Methods Appl Mech Eng 15:63-81. — Reduced integration 处理不可压缩。
  • [X12] Simo JC, Taylor RL (1982) Penalty function formulations for incompressible non-linear analyses. Comp Methods Appl Mech Eng — Penalty 方法。
  • [X13] Srinivasan S, Perucchio R (1994) Finite element analysis of the embryonic heart. — Augmented Lagrangian 在心脏早期形态的应用。
  • [X14] Costa KD, Hunter PJ, Rogers JM, Guccione JM, Waldman LK, McCulloch AD (1996) A three-dimensional finite element method for large elastic deformations in ventricular myocardium. J Biomech Eng 118:452-463. — Isochoric 插值。
  • [X15] Fried I (1982) Finite element computation of large rubber membrane deformations. Int J Numer Methods Eng 18:653-660. — 6.4节膜膨胀的来源。
  • [X16] Oden JT, Key JE (1970) Numerical analysis of finite axisymmetric deformations of incompressible elastic solids. — 4节点三角形轴对称有限元最早工作。
  • [X17] Pian THH, Tong P (1970) Finite element methods in continuum mechanics. Adv Appl Mech 12:1-93. — 增量法的奠基。
  • [X18] Brezzi F, Babuska I (1972-1974) inf-sup 条件的基础工作。
  • [X19] Maas SA, Ateshian GA, Weiss JA (2017) FEBio: History and overview. — 2026 年的开源生物力学 FEA 软件，应在 2002 年版本上加一节"现代 FEA 软件生态" (FEBio, FEBio, BioEng, FEBio, FEBio + UMAT).
  • [X20] Vawter DL (1980) An incremental finite element and large deformation analysis of a pulmonary artery. — Humphrey 实验室早期有限元工作。