非线性有限元课程报告

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 101-=δ。

重复这一过程，以第i 次近似解求出第1+i 次近似解的迭代公式为：()i K K iδ=()R K i i 11-+=δ 1-2直到Δi i iδδδ-=+1变得充分小，即近似解收敛时，终止迭代。

在迭代过程中，得到的近似解一般不会满足式1-1，即()()0≠-=R K ii 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δ处的线性近似公式为ii i ⎪⎭⎫ ⎝⎛∂∂+=+δψψψ1Δiδ。

其中i iδδδψδψ=⎪⎭⎫⎝⎛∂∂=⎪⎭⎫ ⎝⎛∂∂， []n Tn ψψψδδδδψ.....,......,2121⎥⎦⎤⎢⎣⎡∂∂∂∂∂∂=⎪⎭⎫ ⎝⎛∂∂。

引入记号：()ii T iT K K ⎪⎭⎫ ⎝⎛∂∂==δψδ。

假定1+i δ为真实解，则由()()01=∆+=∆+=+i i T i i i i K δψδδψδψ，解出修正量为：()()()i iTi iTi F R K K -=-=∆--11ψδ 1-5 由于这样确定的iδ∆仅考虑了Taylor 级数的线性项，因而按式1-4和1-5求出的新解仍然是近似解。

这样，Newton 法的迭代公式可归纳为：()()()i iTi iTi F R K K -=-=∆--11ψδ iiiT F K ⎪⎭⎫ ⎝⎛∂∂=⎪⎭⎫ ⎝⎛∂∂=δδψ+=+i i δδ1Δi δ 1-6对于单变量的非线性问题()iT K δ是δ-F 曲线上通过点()()δδF ,的切线斜率。

Newton 法的收敛性是好的，但对某些非线性问题，如理想塑性和塑性软化问题，在迭代过程中T K 可能是奇异或病态的，于是T K 的求逆就会出现困难。

为此，可引入一个阻尼因子η，使矩阵I K iiT η+或者成为非奇异的，或者使它的病态减弱。

这里I 为n n ⨯阶的单位矩阵。

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

当iη变大时，收敛速度变慢，当iη→0时，收敛速度最快。

引入iη后，将用下式代替式1-5()i i iT i I K ψηδ1-+-=∆ 1-71.3 修正的Newton-Raphson 法采用直接迭代法和Newton 法求解非线性方程组时，在迭代过程的每一步都需要重新计算iT K 。

如将Newton 法迭代公式中的iT K 改用初始矩阵()0δT T K K =，就成了修正的Newton-Raphson 法。

此时，仅第一步迭代需要完全求解一个线性方程组，并将三角分解后的0T K 存贮起来，以后的每一步迭代都采用公式：()i Ti K ψδ1--=∆ 1-8这样，只需按式1-8右端的iψ进行回代即可。

修正的Newton 法的每一步迭代所用的计算时间较少，但迭代的收敛速度降低。

为了提高收敛速度，可引入过量修正因子iω。

在按式1-8求出iδ∆之后，采用下式计算新解：+=+i i δδ1i ωΔi δ 1-9i ω为大于1的正数。

可以采用一维搜索的方法确定i ω，此时将i δ∆看作n 维空间中的搜索方向，希望在这一方向上找到一个更好的近似值，即使不能得到精确解（使()0=δψ的解），但可通过选择iω使()()0=∆∆+∆i i i iTi δδωδψδ 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 T1 1-13这就是增量法的基本公式，现设1......0210=<<<<=M λλλλ，将λ分成M 个增量：1--=∆m m m λλλ，11=∆∑=m Mm λ 1-14此时，第m 级荷载增量为：--∆=-=∆R R R R m m m m λ1 1-15迭代公式成为()m m T mm m Tm R K R K ∆=∆=∆------11,111,λλδδm m m δδδ∆+=-1 1-16初始值可取00=λ，00=R ，00=δ， m λ∆一般可取等分值。

根据位移增量m δ∆，可求出应变增量m ε∆和应力增量m σ∆，则：m m m εεε∆+=-1 1-17 m m m σσσ∆+=-1 1-180,T K 为初始切线劲度矩阵，m T K ,是对应于第m 级荷载开始时应力状态m σ的切线劲度矩阵，式1-16是基本的增量法，又称Euler 法，对一维问题，其求解过程如图1-2所示。

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

1.4.2 修正的Euler 法将由Euler 法第m 级荷载增量求得的m δ作为中间结果，记为,m δ，它与前一级结果1-m δ加权平均为：(),1m m m δθθδδθθ-+=-- 1-19式中θ为加权系数，由θδ-m 确定θ-m T K ,，并代替式1-16中的1,-m T K ，则得：m m T m R K ∆=∆--1,θδm m m δδδ∆+=-1 1-20上式就是修正Euler 法的基本公式，实际计算步骤为： (1) 由荷载增量 ()()11----=∆m m m R R R θθ 1-21按下式计算中间位移θδ-m ，θθδ----∆=∆m m T m R K 11,，θθδδδ---∆+=m m m 1(2) 计算相应于θδ-m 的劲度矩阵θ-m T K ,。

(3) 施加全部荷载增量m R ∆，按式1-20计算m δ。

2 材料非线性的基本理论在材料非线性问题中，物理方程中的应力—应变关系不再是线性的。

例如在结构中的裂纹尖端存在应力集中现象，当外载荷达到一定数值时该部位进入塑性，而此时结构中的其他部位还处于弹性阶段。

又如长期处于高温条件下工作的结构，将发生蠕变变形，即在荷载和应力保持不变的情况下，变形或应变仍随着时间的进展而继续增长，这也不是线弹性的物理方程所能描述的。

材料非线性问题可以分为两类。

一类是不依赖于时间的弹塑性问题，其特点是当荷载作用以后，材料变形立即发生，并且不再随时间变化。

另一类是依赖于时间的黏弹塑性问题，其特点是荷载作用后，材料不仅立即发生变形，而且变形随时间而继续变化，在荷载保持不变的条件下，由于材料黏性而继续增长的变形称之为蠕变。

本文重点讨论弹塑性问题的有限元分析。

2.1 塑性力学的基本法则 2.1.1 初始屈服条件此条件规定材料开始塑性变形的应力状态。

对于初始各向同性材料，在一般应力状态下开始进入塑性流动的条件是()00=ijF σ，式中ijσ表示应力张量分量，()ijF σ0的几何意义可以理解为九维应力空间的一个超曲面。

1）V . Mises 屈服条件()0321200=-=s ij ij ij s s F σσ其中0s σ是材料的初始屈服应力，ij m ij ij s δσσ-=是偏斜应力张量分量，()33221131σσσσ++=m 是平均正应力。

并且有以下关系0321==σij ij s s 。