在结构非线性分析中, 结构平衡路径曲线上存在两种类型的临界点: 极值点和分叉点, 二者往往分别对应于结构的极值点型失稳和分叉点型失 稳^[1]。当运用有限元分析方法研究结构和连续体的非线性性质时, 应该追踪平衡路径, 进行增量计算, 并识别极值点和分叉点, 这就需要使用延拓算法去求解非线性方程组^[2]。

非线性有限元分析中的荷载增量法, 采用的是荷载控制的 Newton 方法, 其中包含延拓方法的思想。在平衡路径的极限点(力的转向点)附近, 切线刚度矩阵趋于奇异, 计算分析的稳定性和收敛性相当差, 以至无法进行, 于是引入位移控制技术。在位移转向点附近, 位移控制方法亦失效, 为此 Sabir等^[3]采用荷载控制和位移控制相互转换的控制技术。然而, 数值稳定、计算效率更高的方法还是Riks^[4]和 Wempner^[5]建议的弧长法。该方法能够同时控制结构的变形与荷载, 使分析沿着结构的平衡路径顺利进行。经 Crisfield^[6]和 Ramm^[7]的改进, 弧长延拓方法已成为结构非线性分析中平衡路径追踪最主要的算法。

弧长延拓方法初期主要用于弹性结构的大变形和稳定性分析^[8-12], 近年来开始应用于软化材料三维连续体稳定性的有限元分析^[1,13-16]。弧长延拓算法的编程实现往往要根据原有非线性有限元程序的特点, 采用不同的 Newton 迭代格式。本文研究弧长延拓方法 Newton 迭代的标准格式和两种编程实现的实用格式, 并分析它们之间的关系。

1 有限元系统的非线性方程组

由虚功原理导出的有限元系统的平衡方程^[17]如下:

\(\psi (\sigma )=P(\sigma )-R=\mathbf{0}\ \ (1)\),

\(P(\sigma )=\sum{c_{e}^{T}\int_{e}{{{B}^{\text{T}}}}}\sigma \text{d}V\ \ (2)\),

\(R=\sum{c_{e}^{\text{T}}}\int_{e}{{{N}^{\text{T}}}}p\text{d}V+\sum{c_{e}^{\text{T}}}\int_{Se}{{{N}^{\text{T}}}}q\text{d}S。\ \ (3) \)

其中,\(P(\sigma )\)和\(R\)分别是物体应力矢量\(\sigma\)和载荷矢量p, q的等效节点力; N和\({{c}_{e}}\)分别是有限单元的形函数和选择矩阵; \(B\)是单元应变和节点位移之间的转换矩阵; \(\psi (\sigma )=\mathbf{0}\)表示物体应力和外载处于平衡状态; \(\psi (\sigma )\ne \mathbf{0}\)表示物体处于不平衡(失衡)状态, 将\(-\psi (\sigma )=R-P(\sigma )\)称为失衡力矢量。

设有限元系统的自由度总数为\(N\), 则系统节点位移矢量\(a\)是\(N\)维矢量:

\[a={{[{{a}_{1}},{{a}_{2}},...,{{a}_{N}}]}^{\text{T}}}。 \ \ (4) \]

等效节点力\(P(\sigma )\)和\(R\)也是\(N\)维矢量, 如果将式(1)和(2)中应力\(\sigma\)用应变\(\varepsilon \)表示, 进而用位移表示, 可得到以位移\(a\)为变量的平衡方程:

\[\psi (a)=P(a)-R=\mathbf{0}。\ \ (5)\]

等效节点力 P(σ)与应力 σ 之间是线性关系。在线弹性材料情况下, P(a)与 a 之间也是线性关系。对于弹塑性材料, 应力是应变的非线性函数, 也是位移的非线性函数。于是, \(P(a)\)是\(a\)的非线性函数, 也就是说, 式(5)是非线性方程组。

2 早期的延拓算法——载荷增量法

方程式(5)是以全量形式给出的, 通常需要用Newton迭代算法进行数值求解。Newton迭代具有局部收敛性的特征, 要求迭代给定的位移初值\({{a}^{0}}\)离真解较近, 才能保证迭代解序列收敛到真解。当选择接近真解的初值难以做到时, 往往采用下面的延拓方法^[2]来处理。

引入载荷参数λ, 并嵌入式(5)中, 得到含参数λ的平衡方程组:

\[\psi (a,\lambda )=P(a)-\lambda R=\mathbf{0}。\ \ (6)\]

通常参数λ是事先指定的, 并且逐渐增大, 将其离散为

\[\left\{ \begin{align} & 0={{\lambda }_{0}}<{{\lambda }_{1}}<{{\lambda }_{2}}<...<{{\lambda }_{m}}<{{\lambda }_{m+1}}<...<{{\lambda }_{M}}=1,\, \\ & \Delta \lambda =\Delta {{\lambda }_{m}}={{\lambda }_{m+1}}-{{\lambda }_{m}}\, \\ \end{align} \right.\ \ (7)\]

其中,m为施加的载荷步数, Δλ为载荷参数的 增量。

随着λ的增加, 逐步对每个增量△λ求解在λ_m+1下的节点位移a_m+1, 称为求解式(6)的载荷增量法。在每一增量步内, 采用Newton迭代格式, 其初值取为前一增量步的解a_m。只要增量步足够小, 就会满足局部收敛性条件。式(6)的Jacobi矩阵(即系统的切线刚度矩阵)和Newton迭代公式分别为

\[J=\frac{\partial \psi }{\partial a}=\sum{c_{e}^{\text{T}}}\int_{e}{{{B}^{\text{T}}}{{D}_{ep}}B}\text{d}V{{c}_{e}},\ \ (8)\]

\[\left\{ \begin{align} & \Delta {{a}^{0}}=0, \\ & \Delta {{a}^{n+1}}=\Delta {{a}^{n}}-{{J}^{-1}}\psi ({{a}_{m}}+\Delta {{a}^{n}})\ \\ \end{align} \right.\ \ (9)\]

.其中, D_ep是弹塑性本构矩阵。

对于 Jacobi 矩阵是非奇异(detJ ≠ 0)的情况, 在增量步Δλ内迭代可以收敛到真解 a_m+1。当 m = M, λ_M = 1时, 就得到原方程组式(5)的解。这种引入载荷因子λ, 用增量方法求解非线性方程组的思想, 就是延拓方法的思想, 载荷增量法是早期的(或经典的)延拓方法。.

3 弧长延拓算法

对于某些强非线性问题, 例如金属结构大变形引起曲屈和突跳以及材料软化引起岩土结构的不稳定^[1,8,18], 载荷参数λ不再是单调的, 式(7)不再成立。同时, 在载荷参数λ转向点附近, 式(6)的Jacobi 矩阵J可能是奇异或病态的。这种情况下, 采用载荷增量法求解式(6)将会失效。

弧长延拓算法^[4-5]不再将载荷参数\(\lambda \)视为事前已知和单调增加的参数, 而是将它视为与位移分量同样待定的未知变量。于是, 式(6)成为含\(N+1\)个未知变量(\(a\)和\(\lambda \))的\(N\)个方程。为此需要引进一个辅助参数\(s\)(通常是\(N+1\)维空间解曲线的弧长), 并增加一个与位移矢量 a、荷载参数\(\lambda \)和弧长\(s\)关联的约束方程

\[B(a,\ \lambda ,\ s)=0,\ \ (10)\]

才能对问题求解。于是, 问题变为求解(N + 1)个变量的(N + 1) 阶方程组:

\[{{\psi }^{*}}(a,\lambda ,s)=\left\{ \begin{matrix} \psi (a,\lambda ) \\ B(a,\lambda ,s) \\ \end{matrix} \right\}=\mathbf{0}。\ \ (11)\]

辅助参数 s 可选用 N + 1 维空间中解曲线 a(s) 和λ(s)的弧长。弧长微分的平方为

\({{(\text{d}s)}^{2}}=d{{a}^{\text{T}}}\text{d}a+{{(\text{d}\lambda )}^{2}}\), (12)

弧长参数s是单调上升的, 将其离散化为

\[\left\{ \begin{align} & 0\text{=}{{s}_{0}}<{{s}_{1}}<{{s}_{1}}<\ldots <{{s}_{m}}<{{s}_{m+1}}<\ldots <{{s}_{M}}\,, \\ & \Delta s=\Delta {{s}_{m}}={{s}_{m+1}}-{{s}_{m}}\, \\ \end{align} \right.\ \ (13)\]

然后, 按弧长增量\(\Delta s\)逐步地对方程组增量求解。

弧长延拓算法的核心问题是构建具体的约束方程式(10)。约束方程可以采用各种不同的形式。Riks^[19]和 Belytschko 等^[20]在每个弧长步内采用的 方程为

\[B(a,\lambda ,s)=\Delta {{a}^{\text{T}}}\Delta a+{{(\Delta \lambda )}^{2}}-{{(\Delta s)}^{2}}=0,\ \ (14)\]

其中, Δs 是在位移-荷载(N + 1)维空间中给定的弧长增量。

式(14)是式(12)的有限增量表述形式。这时, 扩展后的方程组式(11)的Jacobi矩阵为

\[{{J}^{*}}=\left[ \begin{matrix} J & -R \\ 2\Delta {{a}^{\text{T}}} & 2\Delta \lambda \\ \end{matrix} \right]。\ \ (15)\]

显然, 矩阵\({{J}^{*}}\)不总是正则的。在使用Newton迭代算法时, 迭代初值总是取\(\Delta {{a}^{0}}=\mathbf{0},\ \ \Delta {{\lambda }^{0}}=0\)(\(a_{m+1}^{0}=a_{m}^{{}},\ \lambda _{m+1}^{0}=\lambda _{m}^{{}}\)), 相应的J^*是奇异的。

Crisfield^[6]给出对位移的控制要求Δa^TΔa= (Δs)², 从而写出如式(10)的约束方程:

\(B(a,\lambda ,s)=\Delta {{a}^{\text{T}}}\Delta a-{{(\Delta s)}^{2}}=0\)。 (16)

式(16)是一个球形路径控制的约束方程。为避免方程组式(11)的Jacobi矩阵失去对称性, 需要先导出求解\(\Delta \lambda \)的二次代数方程, 求出Δλ后再求解Δa。在建立二次代数方程时, 方程式的系数含 J^-1。

武际可等^[8]通过位移-荷载参数(N+1)维空间内解曲线的单位切向矢量τ 来建立约束方程。首先, 构造一个(N+1)维矢量:

\[v(a,\lambda )={{[{{v}_{1}},{{v}_{2}},\ldots ,{{v}_{N}},{{v}_{N+1}}]}^{\text{T}}}, \ \ (17)\]

\[{{v}_{i}}={{(-1)}^{i}}\det \left[ \frac{\partial \psi }{\partial {{a}_{1}}},\text{ }\frac{\partial \psi }{\partial {{a}_{2}}},\text{ }\ldots ,\text{ }\frac{\partial \hat{\psi }}{\partial {{a}_{i}}},\text{ }\ldots ,\text{ }\frac{\partial \psi }{\partial {{a}_{N}}},\text{ }\frac{\partial \psi }{\partial \lambda } \right],\ \ (18)\]

式中, “^”表示划去该列。

矢量\(v(a,\lambda )\)是 N+1 维空间中解曲线 a(s)和λ(s)的切向矢量, 而解曲线的单位切向矢量为

\[\tau (a,\lambda )=\frac{v}{||v||},\ \ (19)\]

式中, ||v||代表矢量\(\mathbf{v}\)的模。

对第\(m\)个弧长步Δs, 约束方程为

\[B(a,\lambda ,s)={{\tau }^{\text{T}}}(s)\left\{ \begin{align} & \Delta a(s) \\ & \Delta \lambda (s) \\ \end{align} \right\}-\Delta s=0\, \ \ (20)\]

式中, \(\Delta a=a(s)-{{a}_{m}},\,\,\Delta \lambda =\lambda (s)-{{\lambda }_{m}},\,\,\Delta s={{s}_{m+1}}-{{s}_{m}}\)。

这时, 方程组式(11)的Jacobi矩阵可写为

\[{{J}^{*}}=\left[ \begin{matrix} \frac{\partial \psi }{\partial a} & \frac{\partial \psi }{\partial \lambda } \\ \frac{\partial B}{\partial a} & \frac{\partial B}{\partial \lambda } \\ \end{matrix} \right]=\left[ \begin{matrix} J & {} & -R \\ {} & {{\tau }^{\text{T}}} & {} \\\end{matrix} \right].\]\[=\left[ \begin{matrix} {} & {} & {} & -{{R}_{1}} \\ {} & {} & {} & -{{R}_{2}} \\ {} & J & {} & \vdots \\ {} & {} & {} & -{{R}_{N}} \\ {{\tau }_{1}} & \cdots & {{\tau }_{N}} & {{\tau }_{N+1}} \\\end{matrix} \right]。\ \ (21)\]

在临界点\(a={{a}_{cr}},\ \ \lambda ={{\lambda }_{cr}}\)附近, 矩阵 J 是奇异或病态的, 但是矩阵 J^* 是非奇异的(可证明 det J^*= ||v||>0)。因此, 平衡路径上当遇到转向点时, 可以毫不困难地求解。

本文根据式(20)给出的约束方程, 讨论弧长延拓算法 Newton 迭代的标准格式和有限元程序编程实现时可采用的实用格式。

4 迭代算法的标准格式

为了表述简洁紧凑, 本文采用新的变量:

\[y=\left[ \begin{align} & a \\ & \lambda \\ \end{align} \right]={{\left[ {{a}^{\text{T}}}\quad \lambda \right]}^{\text{T}}}。\ \ (22)\]

式(11)和Jacobi矩阵J^*可分别写为

\[{{\psi }^{*}}(y,s)=\left[ \begin{align} & \psi (y) \\ & B(y,s) \\ \end{align} \right]=\mathbf{0}, \ \ (23)\]

\[{{J}^{*}}(y,s)=\frac{\partial {{\psi }^{*}}}{\partial y}。\ \ (24)\]

按弧长步增量求解。设\(s={{s}_{m}}\)时的解\({{y}_{m}}\)已经得到, 在本增量步内需用Newton法求\(s={{s}_{m+1}}\)时的解\({{y}_{m\text{+}1}}\)。设已经得一个近似解\(y_{m+1}^{n}\),n表示迭代次数, \(m+1\)表示弧长增量步数。为了求得一个更好的解, 将式(11)左端在\(y_{m+1}^{n}\)附近 Taylor 展开, 并截至线性项, 可得

\[\mathbf{0}={{\psi }^{*}}(y,s)\approx {{\psi }^{*}}(y_{m+1}^{n})+{{J}^{*}}(y_{m+1}^{n})(y-y_{m+1}^{n})。\ \ (25)\]

其中, y 的改进解为\(y_{m+1}^{n\text{+}1}\), 迭代初值\(y_{m+1}^{0}\)取自前一弧长步的解y_m, 可得到Newton法标准的迭代公式:

\[\left\{ \begin{align} & y_{m+1}^{0}={{y}_{m}}\ , \\ & y_{m+1}^{n+1}=y_{m+1}^{n}-{{({{J}^{*}}(y_{m+1}^{n}))}^{-1}}{{\psi }^{*}}(y_{m+1}^{n})\ \\ \end{align} \right.\ \ (26)\]

如果在式(26)中, 用矩阵\({{J}^{*}}({{y}_{m}})\)代替\({{J}^{*}}(y_{m+1}^{n})\), 就得到简化Newton法的迭代公式:

\[\left\{ \begin{align} & y_{m+1}^{0}={{y}_{m}}\ , \\ & y_{m+1}^{n+1}=y_{m+1}^{n}-{{({{J}^{*}}({{y}_{m}}))}^{-1}}{{\psi }^{*}}(y_{m+1}^{n})\ \\ \end{align} \right.\ \ (27)\]

如果仅迭代一次, 得到解\(y_{m+1}^{1}\):

\[y_{m+1}^{1}={{y}_{m}}-{{({{J}^{*}}({{y}_{m}}))}^{-1}}{{\psi }^{*}}({{y}_{m}}),\ \ (28)\]

即是自修正Euler方法的解^[18]。

在使用 Newton 法迭代(式(26))和简单 Newton法的迭代(式(27))时, 要根据式(20)和(11)正确地给出\({{J}^{*}}(y_{m+1}^{n})\), \({{J}^{*}}({{y}_{m}})\)和\({{\psi }^{*}}(y_{m+1}^{n})\)的表达式。矩阵\({{J}^{*}}\)的左上角子矩阵 J 是有限元系统的切线刚度矩阵, 将其加边扩展为 J^*(J^*称为广义切线刚度矩阵)。类似地, \({{\psi }^{*}}(y_{m+1}^{n})\)称为广义失衡力矢量, 是 N + 1 维矢量, 前\(N\)个分量构成的矢量是通常意义下的失衡力\(\psi (y_{m+1}^{n})\), 第 N + 1 个分量则为失约弧长\(B(y_{m+1}^{n})\)。随着迭代步数 n 的增加, 广义失衡力矢量趋于零, \((y_{m+1}^{n+1}-y_{m+1}^{n})\to \mathbf{0}\), \(y_{m+1}^{n\text{+}1}\)则收敛到真解\(y_{m+1}^{{}}\)。

Newton法和简化Newton法的迭代公式也可以用增量形式给出。考虑到

\[\left\{ \begin{align} & y_{m+1}^{n}\text{=}y_{m+1}^{0}+\Delta {{y}^{n}}={{y}_{m}}+\Delta {{y}^{n}}\ , \\ & y_{m+1}^{n+1}=y_{m+1}^{0}+\Delta {{y}^{n+1}}={{y}_{m}}+\Delta {{y}^{n+1}}\ , \\ \end{align} \right. (29)\]

Newton法增量形式的迭代公式为

\[\left\{ \begin{align} & \Delta {{y}^{0}}\text{=}0\ , \\ & \Delta {{y}^{n\text{+}1}}\text{=}\Delta {{y}^{n}}-\text{(}{{J}^{*}})_{{{y}_{m}}+\Delta {{y}^{n}}}^{-1}{{\psi }^{*}}({{y}_{m}}+\Delta {{y}^{n}})\ \\ \end{align} \right.\ \ (30)\]

简化Newton法增量形式的迭代公式为

\[\left\{ \begin{align} & \Delta {{y}^{0}}\text{=}\mathbf{0}\ , \\ & \Delta {{y}^{n\text{+}1}}\text{=}\Delta {{y}^{n}}-\text{(}{{J}^{*}})_{{{y}_{m}}}^{-1}{{\psi }^{*}}({{y}_{m}}+\Delta {{y}^{n}})\ \\ \end{align} \right.\ \ (31)\]

根据式(6), 失衡力矢量\(\psi \)通常可写为

\[\psi ({{y}_{m}}+\Delta {{y}^{n}})\ \text{=}P({{a}_{m}}+\Delta {{a}^{n}})-({{\lambda }_{m}}+\Delta {{\lambda }^{n}})R \text{=}P({{\sigma }_{m}}+\Delta {{\sigma }^{n}})-({{\lambda }_{m}}+\Delta {{\lambda }^{n}})R \text{=}P({{\sigma }_{m}})+P(\Delta {{\sigma }^{n}})-{{\lambda }_{m}}R-\Delta {{\lambda }^{n}}R \text{=}\psi ({{y}_{m}})+\psi (\Delta {{y}^{n}}),\ \ (32a)\]

广义失衡力\({{\psi }^{*}}({{y}_{m}}+\Delta {{y}^{n}})\)可写为

\[{{\psi}^{*}}({{y}_{m}}+\Delta {{y}^{n}})=\left[ \begin{matrix} \psi ({{y}_{m}}+\Delta {{y}^{n}}) \\ B({{y}_{m}}+\Delta {{y}^{n}}) \\ \end{matrix} \right]=\left[ \begin{matrix} \psi ({{y}_{m}})+\psi (\Delta {{y}^{n}}) \\ {{\tau }^{\text{T}}}({{y}_{m}}+\Delta {{y}^{n}})\Delta {{y}^{n}}-\Delta s \\ \end{matrix} \right]。\ \ (32b)\]

\[\psi (\Delta {{y}^{n}})=P(\Delta {{\sigma }^{n}})-\Delta {{\lambda }^{n}}R =\sum{c_{e}^{\text{T}}\int\limits_{e}{{{B}^{\text{T}}}\Delta {{\sigma }^{n}}\text{d}V}}-\Delta {{\lambda }^{n}}R。\ \ (33)\]

式中, \(\Delta {{\sigma }^{n}}\)是按弹塑性材料本构关系, 由\(\Delta {{\varepsilon }^{n}}\)计算的应力增量矢量。

采用简化 Newton 法时, 迭代的每一步都是按广义切线刚度 J^*(y_m)重新分配广义失衡力\({{\psi }^{*}}({{y}_{m}}\)\(+\Delta {{y}^{n}})\), 以便恢复平衡和满足约束条件。这种方法称为失衡力迭代方法, 更确切地说是广义失衡力迭代方法。

参照式(27)和(29), 简化Newton法增量形式的迭代公式(30)可以改写为

\[{{J}^{*}}({{y}_{m}})(\Delta {{y}^{n+1}}-\Delta {{y}^{n}})+{{\psi }^{*}}({{y}_{m}}+\Delta {{y}^{n}})=0\ \ (34)\]

其中, J^*(y_m)是\(y={{y}_{m}}\)时的广义切线刚度矩阵, 它是(N+1)×(N+1)的常数矩阵:

\({{J}^{*}}({{y}_{m}})=\left[ \begin{matrix} J({{a}_{m}})-R \\ {{\tau }^{\text{T}}}({{y}_{m}}) \\ \end{matrix} \right]\)。 (35)

其中, 左上角的子矩阵\(J({{a}_{m}})\)是通常意义下的切线刚度矩阵, 为 N ×N 的常数矩阵,\(J({{a}_{m}})=\)\(\sum{c_{e}^{\text{T}}}\int\limits_{e}{{{B}^{\text{T}}}{{D}_{E}}B\text{d}V{{c}_{e}}},\)其中积分项中的 D_E 是在 y= y_m情况下的弹塑性本构矩阵\({{({{D}_{ep}})}_{m}}\), 它在当前弧长步迭代过程中保持不变, 具有弹性本构矩阵的特点(\(\Delta \sigma ={{D}_{E}}\Delta \varepsilon \)), 但它又不是材料的弹性矩阵, 故称为准弹性矩阵。这时有

\[\Delta \sigma _{E}^{n}={{D}_{E}}\Delta {{\varepsilon }^{n}}={{D}_{E}}B\Delta {{a}_{e}}={{D}_{E}}B{{c}_{e}}a,\ \ (36)\]

\[\left\{ \begin{align} & J({{a}_{m}})\Delta {{a}^{n}}=\sum{c_{e}^{\text{T}}}\int\limits_{e}{{{B}^{\text{T}}}\Delta \sigma _{E}^{n}\text{d}V=P(\Delta \sigma _{E}^{n})}, \\ & {{J}^{*}}({{y}_{m}})\Delta {{y}^{n}}=\left[ \begin{matrix} P(\Delta \sigma _{E}^{n})-\Delta {{\lambda }^{n}}R \\ {{\tau }^{\text{T}}}({{y}_{m}})\Delta {{y}^{n}} \\ \end{matrix} \right]\ \\ \end{align} \right.\ \ (37)\]

\({{\psi }^{*}}({{y}_{m}}+\Delta {{y}^{n}})=\left[ \begin{matrix} \psi ({{y}_{m}})+P(\Delta {{\sigma }^{n}})-\Delta {{\lambda }^{n}}R \\ {{\tau }^{\text{T}}}({{y}_{m}}+\Delta {{y}^{n}})\Delta {{y}^{n}}-\Delta s \\ \end{matrix} \right]\)。 (38)

根据式(36)和(37), 简化 Newton 法的方程组式(34)可写为

\[\text{ }\left[ \begin{matrix} J({{a}_{m}})\quad \ -R \\ \tau _{m}^{n} \\ \end{matrix} \right]\left[ \begin{matrix} \Delta {{a}^{n+1}} \\ \Delta {{\lambda }^{n+1}} \\ \end{matrix} \right] =\left[ \begin{matrix} -\psi ({{y}_{m}})+P(\Delta \sigma _{E}^{n}-\Delta {{\sigma }^{n}}) \\ ({{\tau }^{\text{T}}}({{y}_{m}})-{{\tau }^{\text{T}}}({{y}_{m}}+\Delta {{y}^{n}}))\Delta {{y}^{n}}+\Delta s \\ \end{matrix} \right],\ \ (39)\]

相应的迭代公式为

\[\left\{ \begin{align} & \Delta {{y}^{0}}=0, \\ & \Delta {{y}^{n+1}}={{({{J}^{*}}({{y}_{m}}))}^{-1}}\left[ \begin{matrix} -\psi ({{y}_{m}})+P(\Delta \sigma _{E}^{n}-\Delta {{\sigma }^{n}}) \\ {{\tau }^{\text{T}}}({{y}_{m}})\Delta {{y}^{n}}-{{\tau }^{\text{T}}}({{y}_{m}}+\Delta {{y}^{n}})\Delta {{y}^{n}}+\Delta s \\ \end{matrix} \right]\ , \\ \end{align} \right.\ \ (40)\]

式中,

\[P(\Delta \sigma _{E}^{n}-\Delta \sigma _{{}}^{n})=\ \sum{c_{e}^{\text{T}}}\int\limits_{e}{{{B}^{\text{T}}}}(\Delta \sigma _{E}^{n}-\Delta \sigma _{{}}^{n})\text{d}A。\ \ (41)\]

其中, Δσ_E, Δσ和(Δσ_E-Δσ)分别为准弹性应力增量、应力增量(按真实本构计算得到)和准塑性应力增量, P(Δσ_E)称为准弹性应力节点力, P(Δσ_E - Δσ) 称为准塑性应力节点力。

式(40)是由\(\Delta {{y}^{n}}\)直接计算\(\Delta {{y}^{n+1}}\)的简化 Newton法的迭代公式, 迭代过程直至在允许的精度范围内, 达到\(\Delta {{y}^{N+1}}\approx \Delta {{y}^{N}}\)时停止, 并按真实的本构关系, 由\(\Delta {{y}^{N}}\)计算得到最后的应力增量\(\Delta {{\sigma }^{N}}\)和应力\({{\sigma }_{m+1}}\)。在迭代过程中, \(\Delta \sigma _{E}^{n}\)和\(\Delta {{\sigma }^{n}}\)分别是按准弹性关系和按真实的材料本构关系由\(\Delta {{a}^{n}}\)计算得到的应力增量矢量。在迭代过程中, \(P(\Delta \sigma _{E}^{n}-\Delta \sigma _{{}}^{n})\)趋于一个常数矢量\(P(\Delta \sigma _{E}^{N}-\Delta \sigma _{{}}^{N})\), 此矢量是收敛解\(\Delta {{y}^{N}}\)对应的准弹性应力和真实应力之差(Δσ_E - Δσ)的等效节点力, 而(σ_E - σ)可以视为一种“初应力”。于是, 迭代过程相当于寻找一个合适的“初应力”场的过程, 一旦找到“初应力”场, 就按这个“初应力”场和实际荷载求解位移和应变(准线性弹性问题), 就是原来的非线性问题的解。因此, 按式(40)迭代求解的方法称为初应力法^[17]。由于同时考虑失约弧长的迭代, 因 而是广义的初应力方法。

对于 N = 1 (位移自由度为 1)和\({{\psi }_{m}}=0\)的情况, 简化 Newton 法失衡力迭代和初应力迭代的示意图见图1。

5 迭代算法的两种实用格式

武际可等^[8]建议的弧长算法迭代格式分为以下两个步骤。

1) Euler预报。

\[{{y}^{E}}={{y}_{m}}+\tau ({{y}_{m}})({{s}_{m+1}}-{{s}_{m}}), \ \ (42)\]

或

\[\Delta {{y}^{E}}=\tau ({{y}_{m}})\Delta s\text{ },\ \ (43)\]

式中, \({{y}_{m}}\)是弧长 s=s_m 的解, \(\Delta s={{s}_{m+1}}-{{s}_{m}}\)是本弧长步的弧长增量。

2) 迭代修正。

\[\left\{ \begin{align} & \Delta {{y}^{1}}=\Delta {{y}^{E}}, \\ & \Delta {{y}^{n+1}}=\Delta {{y}^{n}}+{{({{J}^{*}}({{y}_{m}}))}^{-1}}\left[ \begin{matrix} -\psi \mathbf{(}{{y}_{m}}+\text{ }\!\!\Delta\!\!\text{ }{{y}^{n}}) \\ 0 \\ \end{matrix} \right]\ , \\ \end{align} \right.\ \ (44)\]

式中, \(-\psi ({{y}_{m}}+\Delta {{y}^{n}})\)是第\(n\)次迭代步的失衡力矢量。失约弧长项取值为0, 表明

\[\tau _{m}^{\text{T}}\Delta {{y}^{E}}=\tau _{m}^{\text{T}}\Delta {{y}^{n}}=\tau _{m}^{\text{T}}\Delta {{y}^{n+1}}=\tau _{m}^{\text{T}}\Delta {{y}^{N}}=\Delta s。\ \ (45)\]

式(45)相当于设定约束条件为 B(Δy, Δs)=\(\tau _{m}^{\text{T}}\Delta y\)\(-\Delta s=0\)。

在每一迭代步, 将近似解增量\(\Delta {{y}^{n}}\)投影到\({{\tau }_{m}}\)方向, 可得到弧长增量。

本文算法采用初应力迭代格式, 将两步合并完成:

\[\left\{ \begin{align} & \Delta {{y}^{0}}=\mathbf{0}\ , \\ & \Delta {{y}^{n+1}}=({{J}^{*}})_{m}^{-1}\left[ \begin{matrix} -\psi ({{y}_{m}})+{{P}^{n}} \\ \Delta s \\\end{matrix} \right],\ \ n=0,\text{ }1,\text{ }2,\cdots \ , \\ \end{align} \right.\ \ (46)\]

式中,

\[{{P}^{n}}=P(\Delta \sigma _{E}^{n}-\Delta \sigma _{{}}^{n})=\sum{c_{e}^{\text{T}}\int\limits_{e}{{{B}^{\text{T}}}}(\Delta \sigma _{E}^{n}-\Delta \sigma _{{}}^{n})\text{d}A\ ;}\]

\(\Delta \sigma _{{}}^{n}\)是按照材料本构关系(\(\Delta \sigma =\int{{{D}_{ep}}}\text{d}\varepsilon \)), 由\(\Delta \varepsilon _{{}}^{n}\)计算得到的应力增量; \(\Delta \sigma _{E}^{n}\)是按照准弹性关系\((\Delta \sigma ={{D}_{E}}\Delta \varepsilon ={{({{D}_{ep}})}_{m}}\Delta \varepsilon ),\) 由\(\Delta \varepsilon \)计算出的应力增量; Pⁿ 为初应力等效节点力矢量。初值\(P_{{}}^{0}\)取为零, 得到的第一次近似解就是 Euler 预报解, Δy¹ = Δy^E。当失约弧长项取为\(\Delta s\)时, 有\(\tau _{m}^{\text{T}}\Delta {{y}^{1}}=\)\(\tau _{m}^{\text{T}}\Delta {{y}^{n}}\)\(=\tau _{m}^{\text{T}}\Delta {{y}^{n+1}}\) \(=\tau _{m}^{\text{T}}\Delta {{y}^{N}}=\Delta s\), 相当于约束条件为

\[\tau _{m}^{\text{T}}\Delta y-\Delta s=0。\ \ (47)\]

上述两种实用格式分别是简化 Newton 法的失衡力迭代格式和初应力迭代格式, 且将约束方程简化为\(\tau _{m}^{\text{T}}\Delta y-\Delta s=0\)。在迭代过程中, 这两种实用格式用\(\tau _{m}^{{}}\)代替\(\tau ({{y}_{m}}+\Delta {{y}^{n}})\), 虽然略去高阶小量, 但在弧长步内不再更新切向量, 从而简化了算法, 利于编程实现, 也能够减小计算量。这两种实用格式采用更简单的约束方程式(47)来替代式(20), 因而迭代的收敛性是有保障的。

6 弧长法在边坡稳定分析中的应用

方建瑞等^[16]将弧长算法应用于边坡的稳定性研究, 边坡材料使用Duncan-Chang模型, 这是一种稳定材料模型, 因而其边坡失稳研究采用的是失衡机制。

基于失稳机制研究边坡的稳定性, 需要建立平衡路径曲线。对于边坡这种多自由度有限元系统, 需要引入广义力和广义位移来描述平衡路径。

设边坡的等效节点荷载为 R, 引入荷载参数λ, 荷载做功的变分为

\[\text{ }\!\!\delta\!\!\text{ }W=\lambda {{R}^{T }}\text{ }\!\!\delta\!\!\text{ }a=\lambda \text{ }\!\!\delta\!\!\text{ (}{{R}^{T }}a\text{)}\ \ （48）\]

如果定义广义力为\(F=\lambda \), 则相应广义位移为\(U={{R}^{T }}a\)。平衡路径曲线\(F=F(U)\)能直观地分析边坡有限元系统的稳定性。

基于初应力迭代格式的弧长法, 编制有限元分析程序, 计算某边坡失稳机制下的安全储备系数(安全系数)K_cr。该边坡有限元网格见图 2, 边坡高度为 9.5 m, 坡角为\(64{}^\circ \), 边坡仅受重力荷载作用, 岩土重度为24 kN/m³, 弹性模量为28.7 GPa, 泊松比为 0.27。单元选用平面应变的三角形网格, 两侧边界上节点的水平位移被约束, 底面节点固定。

采用强度折减法分析该边坡的稳定性, 强度折减系数设为ζ。边坡材料服从 Drucker-Prager 屈服准则:

\[f=\sqrt{{{J}_{2}}}-\alpha (\kappa ){{I}_{1}}-k(\kappa )=0,\ \ (49)\]

式中, I₁ 和J₁分别是应力张量的第一不变量和偏应力张量的第二不变量; 材料强度参数\(k(\kappa )\)和\(\alpha (\kappa )\)是塑性内变量K的函数, 其峰值分别为 k = 61.50 kPa和α = 0.066。选用等效塑性应变为塑性内变量, 边坡材料的应变软化可通过下式的负幂次形式来考虑:

\[\frac{k(\kappa )}{k}=\frac{\alpha (\kappa )}{\alpha }=\frac{1}{4}(3{{\text{e}}^{-{{\kappa }^{2}}}}+1)。\ \ (50)\]

选取强度折减系数\(\zeta =1.00\), 弧长增量\(\Delta s=0.1\), 采用初应力格式的弧长法有限元程序进行增量分析。当荷载参数\(\lambda \)\(=\)\(1.44\)时, 边坡达到临界状态, 此时边坡部分的变形如图3所示(边坡坡顶 A 点的水平位移为-0.01 mm, 竖向位移为0.20 mm), Mises应力和静水应力分布如图4和5所示。分别取折减系数ζ = 1.25, 1.45, 1.50, 依次使用弧长延拓算法, 计算得到相应的平衡路径曲线, 如图 6 所示。曲线上极值点(即临界点)的纵坐标即是稳定性临界点的载荷参数λ_cr。当荷载参数\(\lambda \)=1 时, 对应于边坡上作用有全部重力载荷的情况。图 6 中, ζ = 1.45 对应平衡路径曲线极值点的纵坐标\({{\lambda }_{cr}}=1\), 即全部重力荷载作用下、材料强度折减系数ζ = 1.45 时, 该边坡处于稳定状态的临界点, 即该边坡的安全系数K_cr=1.45。

7 结语

材料软化和构件屈曲往往引起工程结构的不稳定, 对这种稳定性问题的研究需要追踪结构的平衡路径。当采用有限元方法进行数值计算时, 问题最终都归结为一非线性代数方程的求解。如果直接采用Newton迭代法求解非线性有限元方程组, 在平衡路径曲线上的临界点附近, 由于切线刚度矩阵趋于奇异而无法得到收敛解。弧长延拓算法能够很好地解决这一问题。

在非线性有限元程序中编程实现弧长延拓算法时, 需要较为详细的迭代格式。本文给出弧长延拓算法Newton迭代和简化Newton迭代的标准格式, 并讨论了失衡力迭代和初应力迭代两种实用格式及其与标准格式之间的关系。如果弧长增量步是一阶小量, 则这两种实用格式与标准格式在约束方程上的差是二阶小量。由于这两种实用格式在弧长步内不更新切向量, 减少了大量高阶矩阵行列式的计算, 能够减小一个弧长步内的计算量, 同时也能保证实用格式的收敛性。

The authors have declared that no competing interests exist.

---

参考文献

文献选项

• 原文顺序
• 文献年度倒序
• 文中引用次数倒序
• 被引期刊影响因子

[1]	殷有泉. 岩石力学与岩石工程的稳定性.北京: 北京大学出版社, 2011 [本文引用: 3]
[2]	李庆杨, 莫孜中, 祁力群. 非线性方程组的数值解法.北京: 科学出版社, 1987 [本文引用: 2]
[3]	Sabir A B, Lock A C. The application of finite elements to the large-deflection geometrically nonli-near behavior of cylindrical shells // Proceedings of International Conference on Variational Mechanics . Southampton, 1972: Session VII [本文引用: 1]
[4]	Riks E. The application of Newton’s method to the problem of elastic stability . Journal of Applied Mecha-nics, 1972, 39: 1060-1065 [本文引用: 2]
[5]	Wempner G A. Discrete approximation related to nonlinear theories of solids. International Journal of Solids and Structures, 1971, 7: 1581-1599 [本文引用: 2]
[6]	Crisfield M A. A fast incremental/iterative solution procedure that handles snap-through. Computers and Structures, 1981, 13: 55-62 [本文引用: 2]
[7]	Ramm E. Strategies for tracing the non-linear res-ponse near limit points // Wunderlich W, Stein E, Bathe K J. Non-Linear finite element analysis in structural mechanics. Berlin: Springer-Verlag, 1981: 63-89 [本文引用: 1]
[8]	武际可, 苏先樾. 弹性系统的稳定性.北京: 科学出版社, 1994 [本文引用: 4]
[9]	叶康生, 陆天天, 袁驷. 结构几何非线性分析中临界点的直接求解 . 工程力学, 2010, 27(10): 1-6
[10]	罗永峰, 滕锦光, 沈祖炎. 确定结构分支点及跟踪平衡路径的改进弧长法 . 同济大学学报, 1997, 25(5): 502-506
[11]	朱菊芬, 初晓婷. 一种改进的弧长法及在结构后屈曲分析中的应用 . 应用数学和力学, 2002, 23(9): 961-967
[12]	周凌远, 李乔, 李彤梅, 等. 改进弧长法求解屈曲问题 . 西南交通大学学报, 2011, 46(6): 922-925 [本文引用: 1]
[13]	段云岭. 非线性方程组的解法: 局部弧长法 . 力学学报, 1997, 29(1): 116-122 [本文引用: 1]
[14]	陈肖峰, 王建华, 崔海勇. 自适应有限元结合弧长法求解软化问题 . 上海交通大学学报, 2003, 37(12): 1916-1918
[15]	刘国明, 卓家寿, 夏颂佑. 求解非线性有限元方程的弧长法及在工程稳定分析中的应用 . 岩土力学, 1993, 14(4): 57-67
[16]	方建瑞, 李志高, 朱合华. 弧长法在边坡稳定非线性有限元分析中的应用 . 水利学报, 2006, 37(9): 1142-1146 [本文引用: 2]
[17]	殷有泉. 非线性有限元基础.北京: 北京大学出版社, 2007 [本文引用: 2]
[18]	殷有泉. 岩石类材料塑性力学.北京: 北京大学出版社, 2014 [本文引用: 2]
[19]	Riks E. An incremental approach to the solutionof snapping and buckling problems. International Journal of Solids Structures, 1979, 15: 529-551 [本文引用: 1]
[20]	Belytschko T, Liu W K, Brian M.连续体和结构的非线性有限元. 庄茁, 译. 北京: 清华大学出版社, 2002 [本文引用: 1]