非线性有限元方法及实例分析梁军河海大学水利水电工程学院，南京（210098）摘 要：对在地下工程稳定性分析中常用的非线性方程组的求解方法进行研究，讨论了非线性计算的迭代收敛准则，并利用非线性有限元方法分析了一个钢棒单轴拉伸的实例。
关键词：非线性有限元，方程组求解，实例分析1引 言有限单元法已成为一种强有力的数值解法来解决工程中遇到的大量问题，其应用范围从固体到流体，从静力到动力，从力学问题到非力学问题。
有限元的线性分析已经设计工具被广泛采用。
但对于绝大多数水利工程中遇到的实际问题如地下洞室等，将其作为非线性问题加以考虑更符合实际情况。
根据产生非线性的原因，非线性问题主要有3种类型[1]：1．材料非线性问题（简称材料非线性或物理非线性） 2．几何非线性问题3．接触非线性问题（简称接触非线性或边界非线性）2 非线性方程组的求解在非线性力学中，无论是哪一类非线性问题，经过有限元离散后，它们都归结为求解一个非线性代数方程组[2]：()()()00021212211=……==n n n n δδδψδδδψδδδψΛΛΛ （1.1）其中n δδδ,,,21Λ是未知量，n ψψψ,,,21Λ是n δδδ,,,21Λ的非线性函数，引用矢量记号[]T n δδδδΛ21＝ （1.2） []T n ψψψψΛ21= （1.3）上述方程组（1.1）可表示为()0=δψ （1.4）可以将它改写为()()()0=−≡−≡R K R F δδδδψ （1.5）其中()δK 是一个的矩阵，其元素是矢量的函数，n n ×ijk R 为已知矢量。
在位移有限元中，δ代表未知的结点位移，()δF 是等效结点力，R 为等效结点荷载，方程()0=δψ表示结点平衡方程。
在线弹性有限元中，线性方程组0＝－R K δ （1.6）可以毫无困难地求解，但对线性方程组()0=δψ则不行。
一般来说，难以求得其精确解，通常采用数值解法，把非线性问题转化为一系列线性问题。
为了使这一系列线性解收敛于非线性解，曾经有过许多方法，但这些解法都有一定的局限性。
某解法对某一类非线性问题有效，但对另一类问题可能不合适。
因而，根据问题性质正确选用求解方法成为非线性有限元的一个极重要的问题。
常见的求解非线性方程组的数值方法有迭代法、增量法和混合法[3][4][5]。
2.1. 迭代法在每次迭代过程中都施加全部荷载，但逐步修改位移和应变，使之满足非线性的应力－应变关系。
迭代法还分为直接迭代法、Newton-Raphson 方法、修正的Newton-Raphson 方法和拟Newton 法。
①直接迭代法 对非线性方程组()0＝－R K δδ （1.7）设其初始的近似解为，由此确定近似的0δδ＝K 矩阵()00δK K ＝ （1.8）可得改进的近似解()R K 101−＝δ （1.9）重复这一过程，以第i 次近似解求出第1+i 次近似解的迭代公式为()i i K K δ＝ （1.10）()R K i i 11−+＝δ （1.11）直到变得充分小，即近似收敛时，终止迭代。
i i iδδδ−=Δ+1对于单变量问题，这一迭代过程是收敛的，单对于多自由度情况，由于未知量通过矩阵K 耦合，迭代过程可能不收敛。
在迭代过程中，得到的近似解一般不会满足()0＝－R K δδ即()()0≠−≡R K i i i δδδψ （1.12）()δψ作为对平衡偏离的一种度量，称为失衡力。
② Newton-Raphson 方法Newton-Raphson 方法是求解非线性方程组()()()0=−≡−≡R K R F δδδδψ （1.13） 的一个著名方法，简称Newton 法。
设()δψ为具有一阶导数的连续函数，是方程（1.13）的第i 次近似解。
若iδδ＝()()0≠−≡=R F i i i δδψψ （1.14）希望能找到一个更好的、方程（1.13）的近似解为i i i δδδδΔ+=+1＝ （1.14）将（1.15）代入（1.13），并在附近按一阶Tayor 级数展开，则iδδ＝()δψ在iδ处的线性近似公式为iiii δδψψψΔ⎟⎠⎞⎜⎝⎛∂∂+=+1（1.16）引入记号()iiT i TK K ⎟⎠⎞⎜⎝⎛∂∂≡=δψδ Newton 法的迭代公式可归纳为()()(ii TiiTi F R K K −=−=Δ−−11ψδ) （1.17）iii TF K ⎟⎠⎞⎜⎝⎛∂∂=⎟⎠⎞⎜⎝⎛∂∂=δδψ （1.18）i i i δδδΔ+=+1 （1.19）Newton 法的收敛性是比较好的，但对于某些非线性问题，如理想塑性和塑性软化问题，在迭代过程中可能事奇异或病态的，于是的求逆就会出现困难。
为此，可引入一个阻尼因子T K T K η，使矩阵或者成为非奇异的，或者使它的病态减弱。
I K i i Tη+③修正的Newton-Raphson 方法上述两种方法求解非线性方程组时，在迭代过程中的每一步都需要重新计算。
如将Newton 法迭代公式中的改用初始矩阵iT K i T K ()0δT T K K =，就成为修正的Newton-Raphson 法。
此时，仅第一步迭代需要完全求解线性方程组，并将三角分解后的存储起来，以后的每一步迭代都采用公式T K ()i T i K ψδ1−−=Δ （1.20）修正的Newton 法的每一步迭代所用的计算时间较少，但迭代的收敛速度降低。
④拟Newton 法拟Newton 法的主要特点是每次迭代后用一个简单的方法修正K ，K 的修正要满足以下的拟牛顿方程()()()i i i i i K δψδψδδ−=−+++111 （1.21）仿照位移的迭代公式来建立劲度拟矩阵的迭代公式：()()()1111−−−+Δ+=i i i K K K （1.22） 只要由和求出。
iδΔiψΔ()1−Δi K 综上所述，迭代法就是用总荷载作用下不平衡的线性解去逼近平衡的非线性解，迭代的过程就是消除失衡力的过程。
对于所采用的迭代方法，Newton 法收敛最快，修正的Newton 法最慢。
但对于某个具体问题，往往需要进行数值实验，才能判断哪个方法较好。
2.2．增量法增量法与迭代法的最大不同就在于用线性方法求解非线性方程组时，对荷载增量进行线性化处理。
基本思想是将荷载分成许多小的荷载部分（增量），每次施加一个荷载增量。
此时，假定方程是线性的，劲度矩阵K 为常矩阵，对每级增量求出位移增量δΔ，对它累加，就可以得到总位移。
增量法分为Euler 法和修正的Euler 法。
①Euler 法设R 为总荷载，引入荷载因子参数λ，令R R λ=，非线性 方程组称为()()()0＝－＝－＝，R F R F λδδλδψ （1.23）迭代公式成为()mm m T m R R K Δ=Δ−−−111λδδ （1.24）②修正的Euler 法将由Euler 法第级荷载增量求得的m m δ作为中间结果，记为m δ′，它与前一级结果1−m δ加权平均为()'11mm m δθθδδθ−+=−− （1.25）式中θ为加权系数，由θδ−m 确定可得修正的Euler 法的基本公式θ−m T K ,即mm T m R K Δ=Δ−−1,θδm m m δδδΔ+=−1 （1.26）2.3．混合法对于一个实际的非线性问题往往用一种非线性解法是得不到满意得结果的，由此混合使用增量法和迭代法，则称之为混合法或逐步迭代法。
一般在总体上采用Euler 法，而在同一级荷载增量内，采用迭代法。
3 收敛准则在采用迭代法和混合法求解非线性方程组时，必须给出迭代的收敛准则，否则就无法终止计算。
目前，在非线性有限元计算中常用的迭代收敛准则有位移准则、失衡力准则和能量准则。
3.1．位移准则id i δαδ≤Δ （1.27）式中符号表示范数，d α为位移收敛容许误差，是事先指定的一个很小的正数，通常取%5%1.0≤≤d α3.2．失衡力准则()Rq i αδψ≤ （1.28）qα是失衡力容许误差。
当材料表现处明显的软化时，或材料接近理想塑性时，失衡力的微小变化将引起位移增量的很大误差，此时不能采用失衡力误差。
3.3．能量准则能量准则同时考虑位移增量和失衡力，是把每次迭代后的内能增量与初始内能增量相比较。
内能增量指失衡力在位移增量上所做的功。
能量准则可表示为()()()00ψδαψδTe iTi Δ≤Δ （1.29）e α是能量收敛容许误差。
4 非线性分析实例4.1．问题的提出一直径为10mm ，长为100mm 的钢棒，将其沿轴向拉伸20mm ，求其最后的变形情况和应力分布。
弹性模量，泊松比3200E EX =3.0=NUXY 。
塑性时的应力－应变关系如下表1表1 塑性时的应力－应变关系数据点 1 2 3 4 5 6 7 应变 0.002 0.003 0.004 0.005 0.006 0.008 0.01 应力 400 416.3 429 438 445 456.8 465.3 数据点 8 9 10 11 12 13 14 应变 0.015 0.02 0.03 0.05 0.1 0.2 0.3 应力 480.8 491.8 507.4 527.5 555.7 586 603由于模型和载荷都是轴对称的，因此建模时可以只取钢棒轴向截面的1/4，用平面轴对称单元来实现。
4.2．结果分析经过计算得出以下成果：1. 应力－应变曲线如图1所示。
2. 钢棒单轴拉伸等效应力分布如图2。
3. 钢棒X方向应力分布等值线如图3所示。
4. 钢棒Y方向应力分布等值线如图4所示。
分析钢棒受力后的结果，可以看出：（1）对于钢棒的单轴拉伸，由于模型和载荷都是轴对称的，因此建模时可以只取钢棒轴向截面的1/4，用平面轴对称单元来实现。
而且为了产生颈缩现象，在建模时钢棒端部和中部截面作了5％的误差，使其诱导出颈缩现象。
（2）从钢棒单轴拉伸等效应力图上可以看出，最大主应力出现在受拉一侧的拉伸受力处，其值为603.4Mpa；最小主应力出现在钢棒的中部，其值为38.37Mpa。
对于钢棒X方向应力分布可以看出，最大和最小应力的位置发生了变化，出现在钢棒中下部，靠近受拉一侧。
图1 应力－应变曲线图2 钢棒单轴拉伸等效应力分布图3 钢棒X方向应力分布等值线图4 钢棒Y方向应力分布等值线5小结本文主要对在地下工程稳定性分析中常用的非线性方程组的求解方法进行研究，列出了一些基本的求解方法，讨论了非线性计算的迭代收敛准则。
通过一个钢棒单轴拉伸的材料非线性分析实例来说明非线性分析的过程。
参考文献[1] 任青文.非线性有限单元法（上，下）[M].河海大学.2003[2] 吕和祥，蒋和洋.非线性有限元[M].北京：化学工业出版社.1992，10[3] 韦博成，万方焕，朱宏图译.非线性回归分析及应用[M].北京：中国统计出版社.1997[4] Shi G.H. & Goodman,R.E., Underground support design using block theory to determine key block bolts requirements, Proc. Symp. on Rock Mech , in Design of Tunnels, 1983:81-105[5] Shou KJ.A three-dimensional hybrid boundary element method for non-linear analysis of a weak plane near an underground excavation.Tunnelling Underground Space Techno 2000，16(2):215-226The Non-linear Finite Element Method and the ExampleAnalyzeLiang JunWater Conservancy and Engineering College of Hohai University, Nanjing (210098)AbstractTo the commonly used nonlinear simultaneous equation solution method conducts the research in the underground project stability analysis, discussed the nonlinear computation iterative restraining criterion, and has analyzed the example using the non-linear finite element method which a steel rod single axle stretches.Keywords: The non-linear finite element, The system of equations solves, Examples analyze作者简介：梁军（1981—），男，江苏盐城人，在读硕士研究生，从事水工结构方面的研究工作。
