隐式迭代方法

发表于 2023-09-02 更新于 2025-03-25 分类于有限元方法阅读次数：

一般的有限元平衡方程

平衡态的方程可以表示为: \[ 0 = \mathbf{f}^{\rm int}(\mathbf{d}^{n+1}, t^{n+1}) - \mathbf{f}^{\rm ext}(\mathbf{d}^{n+1}, t^{n+1}) \triangleq \mathbf{r}^{n+1}(\mathbf{d}^{n+1}) \] 式中, $\mathbf{f}^{\rm int}$ 和 $\mathbf{f}^{\rm ext}$ 分别表示 $t^{n+1}$ 时刻的内力与外力. 内力与外力是关于节点位移 $\mathbf{d}^{n+1}$ 和时刻 $t^{n+1}$ 的函数. $\mathbf{r}^{n+1}: \mathbb{R}^n \mapsto \mathbb{R}^n$ 表示在时刻 $t^{n+1}$ 关于节点 $\mathbf{d}$ 的非线性函数, 方程 \[ \mathbf{r}^{n+1}(\mathbf{d})=0 \tag{1} \] 的零点即为使得系统在 $t^{n+1}$ 平衡的节点位移.

为求解上述平衡方程, 我们使用 Newton 迭代法. 设迭代的初始位移为 $\mathbf{d}_{(0)}$, 并定义第 $i$ 个迭代步方程的残差为 \[ \mathbf{r}_{(i)}^{n+1} \triangleq \mathbf{r}^{n+1}(\mathbf{d}_{(i)}^{n+1}), \quad i=0,1,2,\ldots \] 在 $\mathbf{d} = \mathbf{d}_{(i)}$ 处展开残差函数 $\mathbf{r}^{n+1}$, 就得到 \[ \mathbf{r}^{n+1}(\mathbf{d}) = \mathbf{r}_{(i)}^{n+1} + \left. \frac{\partial \mathbf{r} ^{\scriptsize n+1}}{\partial \mathbf{d}} \right|_{\mathbf{d} = \mathbf{d}_{(i)}} (\mathbf{d} - \mathbf{d}_{(i)}) + \mathcal{O}(\mathbf{d} - \mathbf{d}_{(i)})^2 \] 忽略关于 $(\mathbf{d} - \mathbf{d}_{(i)})$ 的高阶项, 并将上式代入到待求解的方程 (1) 当中, 就得到关于 $\mathbf{d}$ 的线性近似方程: \[ \mathbf{r}_{(i)}^{n+1} + \left. \frac{\partial \mathbf{r} ^{\scriptsize n+1}}{\partial \mathbf{d}} \right|_{\mathbf{d} = \mathbf{d}_{(i)}} (\mathbf{d} - \mathbf{d}_{(i)}) = 0 \] 为了方便表示上述方程, 重新定义符号: \[ A_{(i)}^{n+1} \triangleq \left. \frac{\partial \mathbf{r} ^{\scriptsize n+1}}{\partial \mathbf{d}} \right|_{\mathbf{d} = \mathbf{d}_{(i)}}, \quad \Delta \mathbf{d}_{(i+1)} \triangleq \mathbf{d} - \mathbf{d}_{(i)} \] 得到 \[ \mathbf{r}_{(i)}^{n+1} + A_{(i)}^{n+1} \Delta \mathbf{d}_{(i+1)} = 0 \tag{2} \] 求解方程 (2) 得到 $\Delta \mathbf{d}_{(i+1)}$, 下一迭代步的初值更新为 \[ \mathbf{d}_{(i+1)} = \mathbf{d}_{(i)} + \Delta \mathbf{d}_{(i+1)} \] 为完成 Newton 迭代法, 还需要设置收敛准则, 在满足准则之后, 跳出 Newton 迭代法.

基于残差的收敛准则

\[ \| \mathbf{r}_{(i)}^{n+1} \| \leq \varepsilon \max \left( \|\mathbf{f}^{\rm int}(\mathbf{d}_{(i)}^{n+1}, t^{n+1})\| , \|\mathbf{f}^{\rm ext}(\mathbf{d}_{(i)}^{n+1}, t^{n+1})\| \right) \]

基于位移增量的收敛准则

\[ \|\Delta \mathbf{d}_{(i+1)}\| \leq \varepsilon \| \mathbf{d}_{(i+1)}^{n+1} \| \]

平衡问题的隐式求解算法陈述如下:

初始化: 节点位移 $\mathbf{d}^{0} := \mathbf{0}$; 单元应力 $\boldsymbol{\sigma}^{0}$; 增量步索引 $n := 0$; 分析步时刻 $t := 0$;
开始计算增量步 $n+1$, 对应增量时间为 $\Delta t^{n+1}$. 初始化 Newton 迭代法: 迭代步索引 $i := 0$, $\mathbf{d}_{(i)}^{n+1} := \mathbf{d}^{n}$, 最大迭代步数 MAXITER
1. 计算系统内力 $\mathbf{f}^{\rm int}(\mathbf{d}_{(i)}^{n+1}, t^{n+1})$ 与外力 $\mathbf{f}^{\rm ext}(\mathbf{d}_{(i)}^{n+1}, t^{n+1})$, 得到内外力之差 $\mathbf{r}_{(i)}^{n+1}$
2. 计算系统刚度矩阵 $A_{(i)}^{n+1}$
3. 求解线性方程组 (2), 得到 $\Delta \mathbf{d}_{(i+1)}$
4. 更新节点位移 $\mathbf{d}_{(i+1)}^{n+1} := \mathbf{d}_{(i)}^{n+1} + \Delta \mathbf{d}_{(i+1)}$; 迭代步索引 $i := i+1$
5. 检查收敛性, 若不收敛, 则回到迭代步骤 1
更新节点位移: $\mathbf{d}^{n+1} := \mathbf{d}_{(i+1)}^{n+1}$; 分析步时刻 $t := t + \Delta t^{n+1}$; 增量步索引 $n := n+1$;

[!NOTE]

不同的收敛准则, 检查收敛性的位置是不同的, 上述算法设定的是基于位移增量的收敛准则

线性问题求解

对于线性问题, 系统的刚度 $A_{(i)}^{n+1} $ 对所有的增量步与迭代步保持恒定, 记为 $K$. 因此, 式 (2) 可以改写为 \[ K \Delta \mathbf{d}_{(i+1)} = \mathbf{f}^{\rm ext}(\mathbf{d}_{(i)}^{n+1}, t^{n+1}) -\mathbf{f}^{\rm int}(\mathbf{d}_{(i)}^{n+1}, t^{n+1}) \] 线性问题不需要使用多个增量步, 只需要在一个增量步中进行迭代即可, 因此将上式中关于增量步的上标去掉: \[ K \Delta \mathbf{d}_{(i+1)} = \mathbf{f}^{\rm ext}(\mathbf{d}_{(i)}) -\mathbf{f}^{\rm int}(\mathbf{d}_{(i)}) \]