跳转至

第四章:有限元方法(Finite Element Method)

阅读笔记

书籍信息:T. Christian Gasser, Vascular Biomechanics (2022), Springer

章节概述:本章系统介绍了线性与非线性有限元方法(FEM),涵盖空间离散化、形状函数、变分法、有限元方程组、约束问题、方程组求解方法,以及在血管生物力学中的实际应用案例。


4.1 引言

有限元方法是求解偏微分方程(PDE)最广泛使用的数值方法,也是血管生物力学问题计算仿真的首选方法。该方法将连续体划分为大量有限元,通过互联节点近似求解。时间步进或预测-校正方法用于计算节点独立变量,如固体力学中的位移u或流体力学中的速度v。

偏微分方程分类

工程问题的物理性质取决于控制偏微分方程的类型:

  • 双曲型PDE:描述波动传播等问题,初始条件的不连续性会在系统中传播,甚至在非线性问题中会局部化形成激波。
  • 椭圆型PDE:描述弹性体变形等问题,解是光滑的,即使初始条件不连续。边界数据影响整个域,边界尖角会导致解的奇异性。
  • 抛物型PDE:描述热传导等问题,介于双曲型和椭圆型之间。

4.2 空间离散化

4.2.1 形状函数(Shape Function)

有限元的节点存储独立问题变量(如结构力学中的位移ui),形状函数Ni(ξ)用于在有限元Ωe内插值这些变量:

\[u(\xi) = \sum_{i=1}^{nnpe} N_i(\xi) u_i \tag{4.1}\]

其中nnpe为有限元节点数,ξ为自然坐标系。

形状函数必须满足的三个性质: 1. 单位性(Property of unity):∑ Ni(ξ) = 1 2. δ性质(Delta property):Ni(ξj) = δij 3. 线性无关性(Linear independence)

4.2.1.1 一维问题形状函数

两节点线性单元(−1 ≤ ξ ≤ 1): $\(N_1(\xi) = (1-\xi)/2 ; \quad N_2(\xi) = (1+\xi)/2 \tag{4.3}\)$

三节点二次单元: $\(N_1(\xi) = (\xi^2-\xi)/2 ; \quad N_2(\xi) = (\xi^2+\xi)/2 ; \quad N_3(\xi) = 1-\xi^2 \tag{4.4}\)$

4.2.1.2 二维问题形状函数

四边形双线性形状函数(自然坐标−1 ≤ ξ1 ≤ 1, −1 ≤ ξ2 ≤ 1): $\(N_1(\xi_1,\xi_2) = \frac{1}{4}(1-\xi_1)(1-\xi_2)\)$ $\(N_2(\xi_1,\xi_2) = \frac{1}{4}(1+\xi_1)(1-\xi_2)\)$ $\(N_3(\xi_1,\xi_2) = \frac{1}{4}(1+\xi_1)(1+\xi_2)\)$ $\(N_4(\xi_1,\xi_2) = \frac{1}{4}(1-\xi_1)(1+\xi_2) \tag{4.5}\)$

三角形单元形状函数(面积坐标): $\(N_1 = \xi_1 ; \quad N_2 = \xi_2 ; \quad N_3 = 1-\xi_1-\xi_2\)$

4.2.1.3 三维问题形状函数

八节点六面体(三线性)单元: $\(N_1 = \frac{1}{8}(1-\xi_1)(1-\xi_2)(1-\xi_3)\)$ $\(N_2 = \frac{1}{8}(1+\xi_1)(1-\xi_2)(1-\xi_3)\)$ $\(N_3 = \frac{1}{8}(1+\xi_1)(1+\xi_2)(1-\xi_3)\)$ $\(N_4 = \frac{1}{8}(1-\xi_1)(1+\xi_2)(1-\xi_3)\)$ $\(N_5 = \frac{1}{8}(1-\xi_1)(1-\xi_2)(1+\xi_3)\)$ $\(N_6 = \frac{1}{8}(1+\xi_1)(1-\xi_2)(1+\xi_3)\)$ $\(N_7 = \frac{1}{8}(1+\xi_1)(1+\xi_2)(1+\xi_3)\)$ $\(N_8 = \frac{1}{8}(1-\xi_1)(1+\xi_2)(1+\xi_3) \tag{4.6}\)$

四节点四面体(线性)单元: $\(N_1 = \xi_1 ; \quad N_2 = \xi_2 ; \quad N_3 = \xi_3 ; \quad N_4 = 1-\xi_1-\xi_2-\xi_3\)$

4.2.2 梯度插值

对称空间梯度: $\(\text{grad}_s u(\xi) = \frac{1}{2}\sum_{i=1}^{nnpe} [u_i \otimes \text{grad}N_i(\xi) + \text{grad}N_i(\xi) \otimes u_i] \tag{4.6}\)$

Voigt记号下梯度插值: $\(\text{grad}_{vs} u = B h \tag{4.7}\)$

4.2.3 混合与杂交有限元

混合有限元方法引入额外变量(多于必要数量)来插值解,可避免不可压缩问题中的体积锁定(volume locking)现象。


4.3 变分法(Calculus of Variations)

变分法是将强形式问题描述转换为弱形式(或积分形式)的一般方法。通过将控制方程乘以测试函数( admissible variation),在物体空间域上积分,并嵌入Dirichlet边界条件来完成。

4.3.1 扩散边值问题

强形式(Poisson方程): $\(\text{div}(\text{grad}c) + \alpha = 0 \quad \text{in} \quad \Omega\)$ $\(c = \bar{c} \quad \text{on} \quad \partial\Omega_c \quad \text{(Dirichlet)}\)$ $\(\text{grad}c \cdot n = q \quad \text{on} \quad \partial\Omega_q \quad \text{(Neumann)} \tag{4.8}\)$

弱形式: $\(\int_\Omega \text{grad}\delta c \cdot \text{grad}c \, d\nu = \int_{\partial\Omega_q} \delta c \, q \, ds + \int_\Omega \delta c \, \alpha \, d\nu \tag{4.9}\)$

4.3.2 对流-扩散边值问题

强形式: $\(v \cdot \text{grad}c - \nu \text{div}(\text{grad}c) + \alpha = 0 \quad \text{in} \quad \Omega \tag{4.10}\)$

弱形式: $\(\int_\Omega \delta c \, v \cdot \text{grad}c \, d\nu + \int_\Omega \nu \, \text{grad}\delta c \cdot \text{grad}c \, d\nu = \int_{\partial\Omega_q} \delta c \, q \, ds + \int_\Omega \delta c \, \alpha \, d\nu \tag{4.11}\)$

4.3.3 线性固体力学边值问题

Cauchy应力σ(u) = C : ε(u),其中C为弹性张量,ε(u) = grads u为工程应变。

强形式(Cauchy动量方程): $\(\text{div}\sigma + b_f - \rho \ddot{u} - \eta \dot{u} = 0 \quad \text{in} \quad \Omega \tag{4.15}\)$

弱形式(虚功原理): $\(\int_{\partial\Omega_t} t \cdot \delta u \, ds + \int_\Omega \delta u \cdot (b_f - \rho \ddot{u} - \eta \dot{u}) \, d\nu - \int_\Omega \sigma : \text{grads}\delta u \, d\nu = 0 \tag{4.17}\)$

4.3.4 非线性固体力学边值问题

大变形、非线性本构模型和非恒定外力导致非线性问题。

内力的变分: $\(\delta \text{int} = \int_\Omega \delta e : \sigma \, d\nu \tag{4.18}\)$

方向导数(线性化): $\(D_u \delta \text{int} = \int_\Omega [\text{grads}\delta u : \text{grads}\Delta u \, \sigma + \text{grads}\delta u : C : \text{grads}\Delta u] \, d\nu \tag{4.19}\)$

几何贡献(初始应力刚度): $\(D_u \delta \text{int}_{geo} = \int_\Omega \text{grads}\delta u : \text{grads}\Delta u \, \sigma \, d\nu \tag{4.20}\)$

材料贡献(材料刚度): $\(D_u \delta \text{int}_{mat} = \int_\Omega \text{grads}\delta u : C : \text{grads}\Delta u \, d\nu \tag{4.21}\)$

4.3.5 不可压缩流动边值问题

Eulerian描述下,流动由以下方程组描述: $\(\text{div}v = 0 \quad \text{in} \quad \Omega \quad \text{(连续性)}\)$ $\(\text{div}\sigma + b_f - \rho \frac{Dv}{Dt} = 0 \quad \text{in} \quad \Omega \quad \text{(动量守恒)} \tag{4.24}\)$


4.4 有限元方程

4.4.1 扩散问题

离散弱形式: $\(\delta h_i \left[ \int_{\Omega_e} B_{ai} B_{aj} \, d\nu \, h_j - \int_{\Omega_e} \alpha N_i \, d\nu + \int_{\partial\Omega_e^q} q N_i \, ds \right] = 0 \tag{4.26}\)$

代数线性方程组: $\(D h - f = 0 \tag{4.27}\)$

4.4.2 对流-扩散问题

离散弱形式: $\((K + D)h - f = 0 \tag{4.29}\)$

其中K为对流矩阵(非对称),D为扩散矩阵(对称)。

4.4.3 线性固体力学问题

二阶微分方程组: $\(M\ddot{h} + C\dot{h} + Kh - f = 0 \tag{4.31}\)$

其中M为质量矩阵,C为阻尼矩阵,K为刚度矩阵。

线性桁架单元: - 质量矩阵:\(M = \frac{\rho Al}{6}\begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}\)(一致质量矩阵) - 集中质量矩阵:\(M_{lumped} = \frac{m}{2}\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}\) - 刚度矩阵:\(K = \frac{EA}{l}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\)

4.4.4 非线性固体力学问题

节点力向量: $\(f_i = \int_{\Omega_e} B_{ai} \sigma_a \, d\nu \tag{4.36}\)$

材料刚度矩阵: $\(K_{mat_{ij}} = \int_{\Omega_e} B_{ai} C_{ab} B_{bj} \, d\nu \tag{4.37}\)$

几何刚度矩阵(初始应力刚度): $\(K_{geo_{ij}} = \int_{\Omega_e} \frac{\partial N_{ai}}{\partial x_c} \sigma_{ab} \frac{\partial N_{cj}}{\partial x_b} \, d\nu \tag{4.38}\)$

压力边界条件(跟随载荷): $\(\delta \text{ext}_p = -p_0 \int_{\partial\Omega_t} \delta u \cdot n \, ds \tag{4.39}\)$

4.4.5 不可压缩流动问题

方程组: $\(M\dot{h} + Ah^2 + Kh - D^T q - f = 0 \tag{4.45}\)$

约束条件: $\(Dh = 0 \quad \text{(连续性/不可压缩性)}\)$

4.4.6 数值积分(高斯-勒让德积分)

体积/面积积分从物理域变换到父域: $\(\int_{\Omega_e} F(x) \, d\nu \approx \sum_{l=1}^{n_{int}} F(\xi^l) \text{det}J(\xi^l) w_l \tag{4.47}\)$

一维高斯-勒让德积分表

n_int 积分点坐标和权重 (ξl, wl)
1 (0.0, 2.0)
2 (±0.57735, 2.0)
3 (0.0, 0.88889); (±0.77460, 0.55556)
4 (±0.86113, 0.34786); (±0.33998, 0.65214)
5 (0.0, 0.56889); (±0.90618, 0.23693); (±0.53847, 0.47863)
6 (±0.66121, 0.36076); (±0.23862, 0.46791); (±0.93247, 0.17132)

高斯-勒让德积分的阶次为2n-1,n个积分点可精确积分n次多项式函数。


4.5 约束问题(Constrained Problems)

许多生物力学问题的本质变量不能"自由"发展,必须满足约束条件。血管问题中流体的不可压缩性即为典型约束。

4.5.1 罚函数法(Penalty Method)

将Cauchy应力分解为σ = σ + σ^vol(偏应力和体积应力),罚函数法通过Cvol >> C来近似不可压缩材料。

优点:易于实现,不增加额外自由度 缺点:Cvol相对越大,矩阵条件数越大,求解越困难

罚函数法将势能泛函P(u)增强,对违反约束C(u)=0的状态添加惩罚项。

4.5.2 Lagrange约束法

Lagrange势能: $\(\mathcal{L}(u,p) = \int_\Omega \Psi(u) \, d\nu + \int_\Omega p(J(u)-1) \, d\nu \tag{4.50}\)$

其中p为Lagrange乘子,可识别为静水压力p = trσ/ndim。

变分陈述给出两个方程,分别对应位移和压力场。

缺点:导致刚度矩阵对角元素为零(pivot elements),无法使用高斯消元法快速求解。

4.5.3 增广Lagrange法(Augmented-Lagrange Method)

增广Lagrange势能: $\(\mathcal{A_L}(u,p) = \mathcal{L}(u,p) - \frac{1}{2\kappa}\int_\Omega p^2 \, d\nu \tag{4.53}\)$

当κ→∞时趋近于标准Lagrange势能。

优点:避免罚函数法的病态条件数问题,同时避免Lagrange法的零pivot问题


4.6 整体化(Globalization)

有限元方程的组装涉及将各单元贡献分配到相应全局节点。

组装操作: $\(K = \bigcup_{e=1}^{n_e} A^e K^e \quad \text{和} \quad f = \bigcup_{e=1}^{n_e} A^e f^e\)$

全局系统: $\(Kh - f = 0 \tag{4.56}\)$

全局刚度矩阵K是稀疏的(sparsely populated),在许多情况下也是对称的。节点编号策略影响矩阵带宽,带宽越小求解越快。

边界条件处理:Dirichlet边界条件(如u1=0)通过移除对应行和列来实现。


4.7 稳定化(Stabilization)

数值稳定性问题可能源于离散化本身,与物理稳定性无关。

4.7.1 有限元刚度正定性

Lyapunov稳定性要求正定刚度矩阵K。

沙漏不稳定(Hourglass Instability):四边形或六面体单元在欠积分时可能出现,需采取稳定化措施。

4.7.2 对流-扩散问题的稳定化

问题:当Peclet数Pe = vh/(2ν)较大时,Galerkin有限元解会出现虚假节点间振荡。

解决方案

  1. 全上游稳定化(Full Upwind Stabilization): 添加人工扩散ν* = vh/2,代价是解过于耗散。

  2. Petrov-Galerkin有限元: 使用不同形状函数插值物理量和测试函数: $\(S_1 = (1+\xi)/2 - 3\beta(1-\xi^2)/4\)$ $\(S_2 = (1-\xi)/2 + 3\beta(1-\xi^2)/4 \tag{4.62}\)$

β控制上游增强程度。

  1. 多维问题推广: 人工扩散张量: $\(\nu^* = \nu^* \frac{v \otimes v}{|v|^2} \tag{4.63}\)$

  2. SUPG(Streamline Upwind Petrov-Galerkin)稳定化: $\(P(\delta c) = v \cdot \text{grad}\delta c \quad \text{和} \quad \tau = \nu^*/|v|^2 \tag{4.65}\)$


4.8 有限元方程组的求解

4.8.1 稀疏线性系统求解

4.8.1.1 直接求解法

LU分解:K = LU,L和U分别为下三角和上三角矩阵。

主元选择策略(避免数值不稳定): - 部分主元(Partial pivoting) - 完全主元(Complete pivoting) - Rook主元

LDU分解:K = LDU,D为非奇异对角矩阵。

4.8.1.2 迭代求解法

定点迭代: $\(Ah_{n+1} = Ahn + f - Kh_n \tag{4.70}\)$

对角预条件可加速收敛。迭代法适合超过约100k未知量的大型系统。

4.8.2 时间积分

显式方法:用tn时刻信息显式计算tn+1时刻的值 隐式方法:需要求解隐式方程组计算tn+1

Euler积分: - 前向Euler(显式):一阶精度 - 后向Euler(隐式):一阶精度

稳定性条件(CFL准则): $\(\Delta t < \alpha \min(h/c) \tag{4.78}\)$

h为最小单元尺寸,c为材料声速。

4.8.3 非线性 formulations

非线性向量方程: $\(g(h) - f(h) = 0 \tag{4.73}\)$

非线性来源:大变形、非线性本构、对流项、跟随载荷等。

4.8.4 增量公式

运动方程: $\(m(\ddot{h}) + d(\dot{h}) + g(h) - f(h) = 0 \tag{4.75}\)$

增量关系: $\(M(h)\Delta\ddot{h} + D(h)\Delta\dot{h} + K(h)\Delta h - K_f(h)\Delta h = 0 \tag{4.76}\)$

4.8.5 显式求解

由MΔö_{n+1} = ... 求解加速度,然后更新速度和位移。

优点:对角集中质量矩阵时避免求解线性方程组 缺点:受CFL条件限制,对非线性问题控制困难

4.8.6 隐式求解

4.8.6.1 Newton-Raphson方法

迭代格式: $\(h \leftarrow h - \delta h \quad \text{其中} \quad K^*(h)\delta h = r(h) \quad \text{且} \quad h^0 = 0 \tag{4.82}\)$

Newton-Raphson在解附近具有二阶收敛速度。

4.8.6.2 载荷增量法

引入载荷因子λ(0 ≤ λ ≤ 1),逐级增加求解。

4.8.6.3 Dirichlet vs Neumann边界条件

位移控制方法:可计算应变软化 载荷控制方法:在极限载荷处终止

4.8.6.4 弧长法(Arc-Length / Continuation Methods)

用于处理" snap-back"行为和极限点问题。

引入载荷因子λ作为额外自由度,广义变量h = [u λ]^T。

弧长约束: $\(\Delta u \cdot \Delta u + \eta^2 \Delta\lambda^2 f \cdot f - \Delta h^2 = r_\lambda \tag{4.89}\)$

伪位移控制是一种简化弧长方法。

4.8.6.5 不可压缩流动方程求解

单体式求解:整体系统联立求解

Chorin投影法(解耦方法): 1. 计算中间速度\(h^*\) 2. 由Poisson方程求解压力q 3. 更新速度h


4.9-4.13 案例研究

4.9 平面双轴试验

问题描述:矩形血管壁样本,四角通过20个锚点连接到致动器。

材料:不可压缩Yeoh材料,Ψ(C) = c1(I1-3) + c2(I1-3)^2 - c1 = 1.1 kPa, c2 = 15.5 kPa

边界条件: - 对称轴:u1 = u2 = 0 - 锚点:直接指定位移

结果: - 中心处:Cauchy应力σ11 = 131.85 kPa, σ22 = 75.36 kPa - 主拉伸:λ1 = 1.513, λ2 = 1.200 - 反力:F1 = 1.5249 N, F2 = 1.2011 N

4.10 圆柱 vessel 充气

问题描述:长L = 8.0 cm, 内半径Ri = 8.0 mm, 外半径Ro = 11.0 mm, 压力pi = 26.0 kPa

材料:不可压缩Yeoh材料,c1 = 18.1 kPa, c2 = 138.0 kPa

验证:膜模型分析,与FEM结果对比

4.11 鼓胀充气实验

问题描述:直径D = 4.0 cm, 厚度H = 1.5 mm圆形贴片,夹持边界,施加内压pi

材料:不可压缩neoHookean材料,Ψ(C) = G(I1-3)/2, G = 422.5 kPa

边界条件处理:避免严重网格畸变,使用线性弹性边界条件(弹簧刚度kbc = 1.0 GPa m⁻¹)

结果:存在极限点,需使用弧长法或伪位移控制求解

4.12 导动脉网络流动

问题描述:腹主动脉以下主要导动脉网络,Newtonian流体,η = 3.5 mPa s

边界条件: - 入口流量:q_inl = 30 ml s⁻¹ - 出口压力:p_out = 13.0 kPa

方法:Stokes流动,Hagen-Poiseuille定律确定血管段阻力

4.13 圆柱管内流动

问题描述:主动脉模型,长l = 30.0 cm, 直径d = 2.4 cm

流体:Newtonian血液,η = 3.5 mPa s, ρ = 1060 kg m⁻³

稳态分析

  • 入口速度:v_inl = 40.0 cm s⁻¹
  • 出口压力:p_out = 13.3 kPa
  • Reynolds数:Re = 2880

脉动分析

  • 入口速度:v_inl = v0 sin(2πt), v0 = 40.0 cm s⁻¹
  • 观察到流动 reversal

瞬态分析与WindKessel边界条件

  • 三元素WindKessel模型
  • R1 = 0.1 mmHg s ml⁻¹, R2 = 1.8 mmHg s ml⁻¹, C = 0.9 ml mmHg⁻¹

公式汇总表

编号 公式名称/内容 公式
4.1 形状函数插值 \(u(\xi) = \sum_{i=1}^{nnpe} N_i(\xi) u_i\)
4.3 1D线性形状函数 \(N_1 = (1-\xi)/2 ; N_2 = (1+\xi)/2\)
4.5 2D双线性形状函数 \(N_i(\xi_1,\xi_2) = \frac{1}{4}(1\pm\xi_1)(1\pm\xi_2)\)
4.6 三线性形状函数 \(N_1 = \frac{1}{8}(1-\xi_1)(1-\xi_2)(1-\xi_3)\)
4.7 梯度插值 \(\text{grad}_{vs} u = B h\)
4.9 扩散问题弱形式 \(\int \text{grad}\delta c \cdot \text{grad}c = \int \delta c q ds + \int \delta c \alpha d\nu\)
4.11 对流-扩散弱形式 \(\int \delta c v \cdot \text{grad}c + \int \nu \text{grad}\delta c \cdot \text{grad}c = ...\)
4.17 虚功原理 \(\int t \cdot \delta u ds + \int \delta u \cdot (b_f - \rho \ddot{u} - \eta \dot{u}) d\nu - \int \sigma : \text{grads}\delta u d\nu = 0\)
4.19 方向导数 \(D_u \delta \text{int} = \int [\text{grads}\delta u : \text{grads}\Delta u \, \sigma + \text{grads}\delta u : C : \text{grads}\Delta u] d\nu\)
4.27 扩散有限元方程 \(Dh - f = 0\)
4.29 对流-扩散有限元方程 \((K + D)h - f = 0\)
4.31 线性固体力学方程 \(M\ddot{h} + C\dot{h} + Kh - f = 0\)
4.36 节点力向量 \(f_i = \int_{\Omega_e} B_{ai} \sigma_a d\nu\)
4.37 材料刚度矩阵 \(K_{mat_{ij}} = \int_{\Omega_e} B_{ai} C_{ab} B_{bj} d\nu\)
4.38 几何刚度矩阵 \(K_{geo_{ij}} = \int_{\Omega_e} \frac{\partial N_{ai}}{\partial x_c} \sigma_{ab} \frac{\partial N_{cj}}{\partial x_b} d\nu\)
4.45 不可压缩流动方程 \(M\dot{h} + Ah^2 + Kh - D^T q - f = 0\)
4.47 数值积分 \(\int_{\Omega_e} F(x) d\nu \approx \sum_{l=1}^{n_{int}} F(\xi^l) \text{det}J(\xi^l) w_l\)
4.50 Lagrange势能 \(\mathcal{L}(u,p) = \int_\Omega \Psi(u) d\nu + \int_\Omega p(J(u)-1) d\nu\)
4.53 增广Lagrange势能 \(\mathcal{A_L}(u,p) = \mathcal{L}(u,p) - \frac{1}{2\kappa}\int_\Omega p^2 d\nu\)
4.56 全局系统 \(Kh - f = 0\)
4.68 线性代数系统 \(Kh - f = 0\)
4.72 刚性微分方程 \(dy/dt = -yt ; y(0) = 1\)
4.73 非线性向量方程 \(g(h) - f(h) = 0\)
4.75 运动方程 \(m(\ddot{h}) + d(\dot{h}) + g(h) - f(h) = 0\)
4.78 CFL稳定性条件 \(\Delta t < \alpha \min(h/c)\)
4.82 Newton-Raphson迭代 \(h \leftarrow h - \delta h\) 其中 \(K^*(h)\delta h = r(h)\)
4.89 弧长约束 \(\Delta u \cdot \Delta u + \eta^2 \Delta\lambda^2 f \cdot f - \Delta h^2 = r_\lambda\)
4.96 不可压缩Yeoh材料Cauchy应力 \(\sigma_1 = \alpha(\lambda_1^2 - \lambda_1^{-2}\lambda_2^{-2})\)

关键概念总结

  1. 有限元方法的核心:将连续问题离散化,通过节点插值近似解
  2. 形状函数:必须满足单位性、δ性质和线性无关性
  3. 变分法:将PDE强形式转换为弱形式,是FEM的理论基础
  4. 约束处理:罚函数法(简单但病态)、Lagrange法(精确但有零pivot)、增广Lagrange法(结合两者优点)
  5. 稳定化技术:对流占优问题的SUPG等方法
  6. 求解策略:显式(条件稳定)vs 隐式(无条件稳定但需迭代)
  7. Newton-Raphson:非线性问题的基本迭代方法
  8. 弧长法:处理snap-through/snap-back等极限点问题

阅读日期:2026年5月21日