第二章非线性代数方程组的解法-Read

在非线性力学中，有多种类型的非线性问题，如材料非线性、几何非线性、接触非线性等。无论是哪一类非线性问题，经过有限元离散后,它们都归结为求解一个非线性代数方程组：

ψ1(12n)0ψ2(12n)0ψn(12n)0其中1,2,,n是未知量，ψ1,ψ2,,ψn是1,2,,n的非线性函数，现引用矢量记号

δ[12n]T ψ[ψ1ψ2ψn]T

上述方程组可表示为

ψ(δ)0

还可以将它改写为

ψ(δ)F(δ)RK(δ)δR0

其元素kij是矢量δ的函数，R为已知矢量。在位移有限元中, δK(δ)是一个nn的矩阵，

代表未知的结点位移，F(δ)是等效结点力，R为等效结点荷载，方程ψ(δ)0表示结点

的平衡方程。

在线弹性有限元中，线性代数方程组

KδR0

可以毫无困难地求解，但对非线性方程组ψ(δ)0则不行。一般来说，难以求得其精确解，

通常采用数值解法，把非线性问题转化为一系列线性问题。为了使这一系列线性解收敛于非线性解，曾经有过许多方法，但这些解法都有一定的局限性。某一解法对某一类非线性问题有效，但对另一类问题可能不合适。因而，根据问题性质正确选用求解方法成为非线性有限元的一个极重要的问题。本章将介绍有限元分析中常见的各种求解非线性方程组的数值方法。

2.1 迭代法

前面已经提到，目前求解非线性方程组的方法一般为线性化方法。若对总荷载进行线性化处理，则称为迭代法。 2.1.1直接迭代法

对非线性方程组

K(δ)δR0 (2-1)

0设其初始的近似解为δδ，由此确定近似的K矩阵

K0K(δ0)

根据式〈2-1〉可得出改进的近似解

δ1(K0)1R

重复这一过程，以第i次近似解求出第i＋1次近似解的迭代公式为

7

KiK(δi)δ直到

i1(K)Ri1 (2-2)

δiδi1δi (2-3)

变得充分小，即近似解收敛时，终止迭代。

在迭代过程中,得到的近似解一般不会满足〈2-1〉式,即

ψ(δi)K(δi)δiR0

ψ(δ)作为对平衡偏离的一种度量，称为失衡力。

图2-1 F~为凸曲线

图2-2 F~为凹曲线

对于一个单变量问题的非线性方程，直接迭代法的计算过程如图2-1和图2-2所示，它们分别给出F~为凸和凹曲线时的迭代过程。可以看出K()就是过曲线上点(,F())

8

与原点的割线斜率。对于单变量问题，这一迭代过程是收敛的，但对多自由度情况，由于未知量通过矩阵K耦合，迭代过程可能不收敛。

2.1.2 Newton—Raphson方法

Newton—Raphson方法是求解非线性方程组

ψ(δ)F(δ)R0 (2-4)

的一个著名方法，简称Newton法。以下将介绍这种方法。

设(δ)为具有一阶导数的连续函数，δδ是方程（2-4）的第i次近似解。若

iψiψ(δi)F(δi)R0

希望能找到一个更好的、方程（2-4）的近似解为

δδi1δiδi (2-5) ii将（2-5）代入（2-4），并在δδ附近按一阶Taylor级数展开，则ψ(δ)在δ处的线性近

似公式为

ψi1ψi(其中

ψii)δ δ(ψiψ)()δδi δδ1ψ()2ψ1ψ2ψn δn引入记号

iKTKT(δi)(ψi) δ假定δi1为真实解，则由

iψ(δi1)ψ(δiδi)ψiKTδi0

解出修正量δ为

i1ii1δi(KT)ψ(KT)(RFi) (2-6)

i由于这样确定的δ仅考虑了Taylor级数的线性项，因而按式（2-6）和（2-5）求出的新解仍然是近似解。这样，Newton法的迭代公式可归纳为

i1ii1δi(KT)ψ(KT)(RFi)iKT(iδi1通过点(ψiF)()iδδδiδi (2-7)

对于单变量的非线性问题，其迭代过程见图2-3和2-4，可以看出KT()是F~曲线上

F())的切线斜率

9

Newton法的收敛性是好的，但对某些非线性问题，如理想塑性和塑性软化问题，在迭代过程中KT可能是奇异或病态的，于是KT的求逆就会出现困难。为此，可引入一个阻尼因子，使矩阵KTI或者成为非奇异的，或者使它的病态减弱。这儿I为nn阶的单位矩阵。的作用是改变矩阵KT主对角线元素不占优的情况。当变大时，收敛速度变

ii

iii

图2-3

图2-4

慢，当→0时，收敛速度最快。引入后，将用下式代替（2-6）

ii

iδi(KTiI)1ψi (2-8)

2.1.3修正的Newton-Raphson法

采用直接迭代法和Newton法求解非线性方程组时，在迭代过程的每一步都需要重新计算K。如将Newton法迭代公式中的

i0KT改用初始矩阵KTKT(δ)，就成了修正的Newton-Raphson法

i

T0K

R （简称修正的Newton法）。此时，仅第一步迭代需要完全求解一个线性方程组，并将三角分解后的KT存贮起来，以后的每一步迭代都采用公式

0



0

01iδi(KT)ψ (2-9) 图2－5

i这样，只需按式(2-9)右端的ψ进行回代即可。

修正Newton法的每一步迭代所用的计算时间较少，但迭代的收敛速度降低。为了提高收敛速度，可引入过量修正因子w。在按（2-9）式求出δ之后，采用下式计算新解

iiδi1δiwiδi (2-10)

10

wi为大于1的正数。可以采用一维搜索的方法确定wi，此时将δi看作n维空间中的搜索方向，希望在这一方向上找到一个更好的近似值，即使不能得到精确解（使ψ(δ)0的解）,

i但可通过选择w使ψ(δ)在搜索方向上的分量为零，即

(δi)Tψ(δiwiδi)0 (2-11)

i这是一个关于w的单变量非线性方程。

0

在应用修正的Newton法时，还可以在每经过若干次迭代后再重新计算一个新的KT，

也可达到提高收敛速度的目的。

2.1.4拟Newton法

前面所谈的Newton法，每次迭代后需要重新计算一个新的矩阵KT，而修正的Newton法保持KT不变。拟Newton法的主要特点是每次迭代后用一个简单的方法修正K，K的修正要满足以下的拟牛顿方程

0

Ki1(δi1δi)ψ(δi1)ψ(δi) (2-12)

对于单变量情况，上式中的Ki1是导数(ψδ)i的近似表达式，实际上就是割线劲度矩阵。由图2-6可知

0(K0)10(K0)1(RF0)

100

010 (K)1100FF111(K1)1(RF1)

┈┈ 图2-6 拟Newton法

(Ki11ii1i (2-13) )i1ii显然Ki1就是相应于ii1i与ii1i的割线劲度矩阵。但实际上对于

i多维情况，无法由（2-13）式求出K。.我们可仿照位移的迭代公式来建立劲度矩阵逆矩阵

的迭代公式：

(Ki1)1(Ki)1(Ki)1 (2-14)

ii那么只要由δ和ψ求出(K)，就可以确定(Ki1i11)。.这儿，修正矩阵(Ki)1的

秩m≥1，通常取m=1或2。对于秩为m的NN阶矩阵，总可以将它表示为ABT的形式，A和B均为Nm阶矩阵。得到(K(1) 秩1算法

修正矩阵(K)表示为

i1)后，再由它求出δi

δi(Ki1)1ψi (2-15)

i11(Ki)1ABT (2-16)

其中A和B均为N×1阶向量。将（2-16）式代入（2-14）后，再将（2-14）式代入（2-15）式可得

11

ABTψiδi(Ki)1ψi

若Bψ0，则

TiA[δi(Ki)1ψi]/BTψi (2-17)

将（2-17）式代入（2-16）得

(Ki)1[δi(Ki)1ψi]BT (2-18)

式中

1BTψi0TiTi1当ψ当ψii0时0时 (2-19)

若取B(ψ)(K)，由（2-18）和（2-19）式得

(ψi)T(Ki)1i （当ψ0时） (2-20) (K)[δ(K)ψ]iTi1i(ψ)(K)ψi1ii1i根据上式和式（2-14）求出的(Kii1ii11) 是不对称的，因而式（2-20）是非对称秩1算法。

若取Bδ(K)ψ，由（2-18）和（2-19）式可得

[δi(Ki)1ψi]Tii1iδ(K)ψ （当时） (K)[δ(K)ψ]ii1iTi[δ(K)ψ]ψi1ii1i(2-21)

可以看出，只要初始逆矩阵(K)是对称的，那么按式（2-21）和（2-14）求出的(K总是对称矩阵。所以式（2-21）是对称秩1算法。

(2) 秩2算法

一个NN阶的秩2矩阵，总可以表示为

01i11)

(K)A1i1B1TTA2TA1B1TA2B2 (2-22)

B2式中A1、A2、B1和B2均为N×1维向量。将上式代入（2-14），再代入（2-15）得

TA1B1TψiA2B2ψiδi(Ki)1ψi (2-23)

为满足（2-23）式，可取

δi(Ki)1ψi A2 A1TiTiB1ψB2ψ代回（2-22）式得

T(Ki)11δiB1T2(Ki)1ψiB2 (2-24)

其中

11B1Tψi0

当ψ当ψii0时0时 (2-25)

12

1T2B2ψi0为了使它具有更普遍的意义，考虑作以下的变换

当ψ当ψii0时0时 (2-26)

TTB1TB2 B2T BTB1ψiB2ψiT1则显然有B1ψiB2ψi1 (2-27) 于是式（2-24）成为

TT(Ki)1δiB1(Ki)1ψiB2 (2-28)

引入参数，将B1和B2取为如下的组合形式

TTTT(δi)TiTi1B[1(ψ)(K)ψ](ψ)(K) iTi(δ)ψT1iTi1i（(δ)ψ0） (2-29)

iTi(ψi)T(Ki)1iTB[1(δ)ψ](δ) iTi1i(ψ)(K)ψT2iTi（(ψ)(K)ψ0） (2-30)

显然，这样选择的B1和B2满足（2-27）式。现将（2-29）和（2-30）代入（2-28）式得

TTiTi1i(K)i1δi(δi)T(Ki)1ψi(ψi)T(Ki)1(δi)Tψi(ψi)T(Ki)1ψiiTi1ii1iiTi1δi(δi)TiTi(K)ψ(ψ)(K) (2-31) [(ψ)(K)ψ(δ)ψiTiiTi1i(δ)ψ(ψ)(K)ψδi(ψi)T(Ki)1(Ki)1ψi(δi)T]可以看出，只要初始逆矩阵(K)是对称的，那么按（2-14）和（2-31）式得到的(K总是对称矩阵，因而（2-31）式是对称的秩2算法。

如令式（2-31）中的=0，便得到DFP〈Davidon-Fletcher-Powell〉公式

01i11)δi(δi)T(Ki)1ψi(ψi)T(Ki)1 (2-32) (K)(δi)Tψi(ψi)T(Ki)1ψii1如令((δ)ψ)，可得 BFS（Broyden-Fletcher-shanno）公式

iTi1 13

(K)ii1(ψi)T(Ki)1ψiδi(δi)T(1)iTiiTi(δ)ψ(δ)ψiTi1i1iiTδ(ψ)(K)(K)ψ(δ)(δi)Tψi (2-33)

式（2-20）、（2-21）、（2-32）和（2-33）中的任一个与式（2-14）、（2-15）联立，便构成解非线注方程组的拟Newton法。实践表明，BFS的秩2算法是目前最成功的算法之一，它具有较好的数值稳定性。

从以上迭代法的计算过程可以看出，随着迭代次数的增加，失衡力ψKδR逐渐减小，并趋于平衡。可见迭代法就是用总荷载作用下不平衡的线性解去逼*衡的非线性解，迭代的过程就是消除失衡力的过程。对于不同的迭代方法，这一过程的快慢也就是收敛速度是不同的。一般来说，Newton法最快，拟Newton法次之，修正的Newton法最慢。理论上可证明，Newton法的收敛速度为二次，修正的Newton法收敛速度只有一次，BFS秩2拟Newton法的收敛速度介于一次和二次之间。不过，各种方法的效率不仅与收敛速度有关，还与每一步迭代所需的计算量有关。Newton法每一步的计算量最大，拟Newton法次之，修正Newton法最小。因而对于某个具体问题，往往需要进行数值实验，才能判断哪个方法较高。一般来说，不同问题可选用不同方法，究竟用哪一种？与所研究问题的性质、计算规模及容许误差等因素有关。

iii2.2 增量法

在用线性方法求解非线性方程组时，若对荷载增量进行线性化处理，则称增量法。它的基本思想是将荷载分成许多小的荷载部分（增量），每次施加一个荷载增量。此时，假定方程是线性的，劲度矩阵K为常矩阵。对不同级别的荷载增量，K变化的。这样，对每级增量求出位移增量δ，对它累加，就可得到总位移。实际上就是以一系列的线性问题代替了非线性问题。 2.2.1 Euler法

设R为总荷载，引入参数一一荷载因子，令R=R，则非线性方程组成为

ψ(δ,)F(δ)RF(δ)R0 (2-34) 问题成为对一个任意给定的（≥0），求δδ()。现设δ是对应于的解，而δδ是对应于的解，则有

ψ(δ,)ψ(δδ,)0 (2-35)

对上式按Taylor级数展开，

ψ(δδ,)ψ(δ,)略去高次项，令KT(δ,)ψψδ δψψR和式（2-34），并注意，可得

δ1δKTR (2-36)

这就是增量法的基本公式，现设

0012M1

将分成M个增量

14

mmm1此时，第m级荷载增量为

m1Mm1 (2-37)

RmRmRm1mR (2-38)

迭代公式成为

11δmKT(δm1m1)RmKT,m1Rmδmδm1δm0，R00，δ00。m一般可取等分值。根据位移增量δm，可求

初始值可取0εmσm出应变增量

和应力增量

，则

(2-39)

εmεm1εm (2-40)

σmσm1σm (2-41)

KT,0为初始切线劲度矩阵，KT,m是对应于第m级荷载开始时应力状态σm的切线劲度矩阵，

式（2-39）是基本的增量法，又称Euler法，对一维问题，其求解过程如图2-7所示。

（α） （b）

图2-7

从图2-7可以看出，每步计算都会引起偏差，使折线偏离曲线，解答产生漂移，随着求解步数的增加，由于偏差的积累使最后的解答离开真解较远，从而降低了计算精度，为此须对这一方法做些改进。 2.2.2修正的Euler法

将由Euler法第m级荷载增量求得的δm作为中间结果，记为δm，它与前一级结果δm1加权平均为

'δmδm1(1)δm (2-42)

式中为加权系数，由δm确定KT,m，并代替式（2-39）中的KT,m1则得

'1δmKT,mRmδmδm1δm上式就是修正Euler法的基本公式，实际计算步骤为

(1) 由荷载增量

(2-43)

Rm(1)(RmRm1) (2-44)

15

按下式计算中间位移δm

1δmKT,m1Rmδmδm1δm(2) 计算相应于δm的劲度矩阵KT,m

(3) 施加全部荷载增量Rm，按式（2-43）计算δm 单变量的修正的Euler求解过程见图2-8。

当＝0.5时，修正的Euler法就是Runge-Kutta法。

(2-45)

2.3 混合法

如对同一非线性方程组混合使用增量法和迭代法，则称为混合法或逐步迭代法。一般在总体上采用Euler增量法，而在同一级荷载增量内，采用迭代法。 2.3.1 Euier--Newton法

在增量步内采用Newton迭代法。现以δm和δm分别表示第m级荷载增量时δ的初值和终值，以Rm表示第m级增量时R的终值，则由式（2-7）得第m增量步的迭代公式

0δmδm1

0RmmR

ii1ii1iδm(KT,m1)(mRFm)(Kt,m1)(mRψm)δi1mδδimim (2-46)

逐步迭代过程可由图2-9所示的一维问题清楚地看出。

如果每一增量步内只迭代一次，此时

1δmδm 0δmδm

则对第m增量步有

010δm(KT,m1)(mRψm)δmδm1δm

(2-47)

这实际上是对Euler法所产生的，与真解偏差的修正，因而称为自修正方法。

16

当前一增量步的计算结果精确时，式（2-34）成立，则自修正方法回到Euler法。 2.3.2 Euler-修正Newton法

在每一增量步内，采用修正的Newton法，取其初始的切线劲度矩阵作为这一增量步内不变的劲度矩阵，则由于

i0KT,m1KT,m1

所以，在第m增量步内

0δmδm1 RmmR

i01i01iδm(KT,m1)(mRFm)(Kt,m1)(mRψm)δ

i1mδδimim (2-48)

按同样的思路，还可给出Euler-拟Newton法。

2.4 迭代收敛准则及增量法的步长选择

在采用迭代法或混合法求解非线性方程组时，必须给出迭代的收敛准则，否则就无法终止迭代计算。收敛准则取得不合适，会使计算结果不精确或多费机时。目前，在非线性有限元计算中常用的迭代收敛准则有以下三种

1． 位移准则

δidδi (2-49)

式中符号

表示范数，它的定义是

δ(δTδ)12 (2-50) δ(δTδ)12 (2-51)

d为位移收敛容差，是事先指定的一个很小的正数，通常取

0.1%≤d≤5%

当材料硬化明显时，位移增量的微小变化将引起失衡力的很大偏差。还有，当相邻两次迭代得到的位移增量范数之比跳动较大时，会使一个应当能收敛的问题判定为不收敛。对于这两种情况，不能使用位移收敛准则。

2． 失衡力准则

ψ(δi)qR (2-52)

q是失衡力收敛容差。当失衡力很小时，可认为式（2-1）满足，逼近了真解。

当材料表现出明显软化时，或材料接近理想塑性时，失衡力的微小变化将引起位移增量的很大偏差，此时不能采用失衡力准则。

3．能量准则

这是一种比较好的收敛准则，因为它同时考虑了位移增量和失衡力。能量准则是把每次迭代后的内能增量与初始内能增量相比较。内能增量指失衡力在位移糟量上所做的功。这一准则可表示为

(δi)Tψie((δ0)Tψ0) (2-53)

e是能量收敛容差。

在用增量法或混合法求解非线性方程组时，需要合理选择荷载增量的步长。步长太大使计算不收敛，步长太短则增加了计算时间。

17

由于增量法是一种线性化方法，应当根据问题的非线性程度来选择合理的步长。一般来说，随着非线性程度的增大，步长要减小。

P. G. Bergan 等提出根据不同荷载阶段的结构劲度，计算增量步长的方法。初始（线弹性）劲度可用下式度量

R1TR1 (2-54) STδ1R1*1第m增量步的劲度用式（2-55）度量

TRmRm （m＝1, 2, …） (2-55) STδmRm*m第m增量步的劲度参数则表示为

*SmSmS1* (2-56)

由于

R11R RmmR

由式（2-54）～（2-56）可得

Smmδ1TR1δRTm (2-57)

这样就可按下面的递推公式选择荷载增量的步长

mm1S (m≥3) (2-58)

Sm2Sm1式中S是劲度参数的变化值，事先给出，可在0.05~0.2范围内选取。如果给定进入非线性状态后的第一个荷载增量因子2，就可按上式确定3，4，……。

18