非线性有限元课程报告

1 非线性问题的求解方法

无论是哪一类非线性问题，经过有限元离散后,它们都归结为求解一个非线性代数方程组：

()0......,211=n δδδψ

()0......,212=n δδδψ

……

()0......,21=n n δδδψ

其中：1δ,2δ……n δ是未知量，1ψ,2ψ……n ψ是1δ,2δ……n δ的非线性函数，现引用矢量记号

[]T n δδδδ (21)

[]T n ψψψψ (21)

上述方程组可表示为： ()0=δψ

还可以将它改写为：

()()()0=-=-=R K R F δδδδψ

()δK 是一个n n ⨯的矩阵，其元素ij k 是矢量δ的函数，R 为已知矢量。在位移有限元中，δ代表未知的结点位移，()δF 是等效结点力，R 为等效结点荷载，方程()0=δψ表示结点的平衡方程。

在线弹性有限元中，线性代数方程组0=-R K δ可以毫无困难地求解，但对非线性方程组()0=δψ则

不行。一般来说，难以求得其精确解，通常采用数值解法，把非线性问题转化为一系列线性问题。为了使这一系列线性解收敛于非线性解，曾经有过许多方法，但这些解法都有一定的局限性。某一解法对某一类非线性问题有效，但对另一类问题可能不合适。因而，根据问题性质正确选用求解方法成为非线性有限元的一个极重要的问题。 1.1 直接迭代法

目前求解非线性方程组的方法一般为线性化方法。若对总荷载进行线性化处理，则称为迭代法。 对非线性方程组：

()0=-R K δδ 1-1

设其初始的近似解为0

δδ=，由此确定近似的K 矩阵()0

δ

K K =，根据式1-1可得出改进的近似解

()R K 1

01-=δ。重复这一过程，以第i 次近似解求出第1+i 次近似解的迭代公式为：

()i K K i

δ=

()R K i i 1

1

-+=δ 1-2

直到Δi i i

δδ

δ-=+1

变得充分小，即近似解收敛时，终止迭代。在迭代过程中，得到的近似解一般不会满足式1-1，即()()0≠-=R K i

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

对于一个单变量问题的非线性方程，δ-F 为凸曲线时，直接迭代法的计算过程如图1-1所示，可以看出()δK 就是过曲线上点()()δδF ,与原点的割线斜率。对于单变量问题，这一迭代过程是收敛的，但对多自由度情况，由于未知量通过矩阵耦合，迭代过程可能不收敛。同样可以得出δ-F 为凹曲线时的情况。

图1-1

δ-F 为凸曲线

1.2 Newton-Raphson 方法

Newton —Raphson 方法是求解非线性方程组 ()()0=-=R F δδψ 1-3

的一个著名方法，简称Newton 法。

设()δψ为具有一阶导数的连续函数，i

δ

δ=是方程1-3的第i 次近似解。若

()()0=-==R F i i i δδψψ，希望能找到一个更好的、方程1-3的近似解为

+==+i i δδδ1Δi δ 1-4 将式1-4代入式1-3，并在i δδ=附近按一阶Taylor 级数展开，则()δψ在i

δ处的线性近似公式为

i

i i ⎪⎭

⎫ ⎝⎛∂∂+=+δψψψ1

Δi

δ。

其中i i

δδδψδψ=⎪⎭⎫

⎝⎛∂∂=⎪⎭⎫ ⎝⎛∂∂， []n T

n ψψψδδδδψ.....,......,2121⎥⎦

⎤⎢⎣⎡∂∂∂∂∂∂=⎪⎭⎫ ⎝⎛∂∂。 引入记号：()

i

i T i

T K K ⎪⎭⎫ ⎝⎛∂∂==δψδ。

假定1

+i δ为真实解，则由()()01=∆+=∆+=+i i T i i i i K δψδδψδψ，解出修正量为：

()

()()i i

T

i i

T

i F R K K -=-=∆--1

1

ψδ 1-5 由于这样确定的i

δ∆仅考虑了Taylor 级数的线性项，因而按式1-4和1-5求出的新解仍然是近似解。这样，Newton 法的迭代公式可归纳为：

()()()i i

T

i i

T

i F R K K -=-=∆--1

1

ψδ i

i

i

T F K ⎪⎭

⎫ ⎝⎛∂∂=⎪⎭⎫ ⎝⎛∂∂=δδψ

+=+i i δδ1Δi δ 1-6

对于单变量的非线性问题()i

T K δ是δ-F 曲线上通过点()()δδF ,的切线斜率。Newton 法的收敛性是

好的，但对某些非线性问题，如理想塑性和塑性软化问题，在迭代过程中T K 可能是奇异或病态的，于是T K 的求逆就会出现困难。为此，可引入一个阻尼因子η，使矩阵I K i

i

T η+或者成为非奇异的，或者使它的病态减弱。这里I 为n n ⨯阶的单位矩阵。i

η的作用是改变矩阵i T K 主对角线元素不占优的情况。当i

η变大时，收敛速度变慢，当i

η→0时，收敛速度最快。引入i

η后，将用下式代替式1-5

()

i i i

T i I K ψηδ1

-+-=∆ 1-7

1.3 修正的Newton-Raphson 法

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

T K 。如将Newton 法迭代公式中的i

T K 改用初始矩阵()0

δ

T T K K =，就成了修正的Newton-Raphson 法。此时，仅第

一步迭代需要完全求解一个线性方程组，并将三角分解后的0

T K 存贮起来，以后的每一步迭代都采用公式：

()

i T

i K ψδ1

--=∆ 1-8

这样，只需按式1-8右端的i

ψ进行回代即可。修正的Newton 法的每一步迭代所用的计算时间较少，但迭代的收敛速度降低。为了提高收敛速度，可引入过量修正因子i

ω。在按式1-8求出i

δ∆之后，采用下式计算新解：

+=+i i δδ1i ωΔi δ 1-9

i ω为大于1的正数。可以采用一维搜索的方法确定i ω，此时将i δ∆看作n 维空间中的搜索方向，希

望在这一方向上找到一个更好的近似值，即使不能得到精确解（使()0=δψ的解），但可通过选择i

ω使

()()

0=∆∆+∆i i i i

T

i δδωδ

ψδ 1-10

这是一个关于i

ω的单变量非线性方程。在应用修正的Newton 法时，还可以在每经过若干次迭代后再重新计算一个新的，也可达到提高收敛速度的目的。 1.4 增量法

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

设R 为总荷载，引入参数λ荷载因子，令-

=R R λ则非线性方程组成为

()()()0,=-=-=-

R F R F λδδλδψ 1-11

问题成为对一个任意给定的λ（0≥λ），求()λδδ=。现设δ是对应于λ的解，而δδ∆+是对应于λλ∆+的解，则有：

()()0,,=∆+∆+=λλδδψλδψ 1-12

对上式按Taylor 级数展开：

()()......,,+∆∂∂+∆∂∂+

=λλ

ψδδψλδψλδψ 略去高次项，令()δψλδ∂∂=,T K ，并注意--=∂∂R λ

ψ

和式1-11，可得：

λδ∆=∆-

-R K T

1 1-13

这就是增量法的基本公式，现设1......0210=<<<<=M λλλλ，将λ分成M 个增量：

1--=∆m m m λλλ，11=∆∑=m M

m λ 1-14