從一實(shí)例分析看弧長法與牛頓—拉普森法

魏鵬，李建光，李延強(qiáng)，陳超

?

從一實(shí)例分析看弧長法與牛頓—拉普森法

魏鵬，李建光，李延強(qiáng)，陳超

（西南林業(yè)大學(xué)土木工程學(xué)院，云南 昆明 650000）

在求解的非線性有限元方程中，牛頓—拉普森法和弧長法是兩類重要的方法。牛頓—拉普森迭代法只能跟蹤位移載荷曲線的上升段，但無法跟蹤極值點(diǎn)以后的位移—載荷路徑，而弧長法可以全程跟蹤位移—載荷路徑。對牛頓—拉普森法和弧長法的原理以及實(shí)施步驟進(jìn)行了回顧，然后通過MATLAB對一則本構(gòu)關(guān)系為非線性的算例使用牛頓—拉普森法和弧長法進(jìn)行了計(jì)算并與理論解做了比較。數(shù)值結(jié)果表明，弧長法能很好跟蹤全過程位移—載荷路徑，并取得較好的效果，而牛頓—拉普森法在跟蹤到極值點(diǎn)時(shí)產(chǎn)生發(fā)散而無法繼續(xù)跟蹤極值點(diǎn)以后的路徑。

牛頓—拉普森法；弧長法；非線性有限元；MATLAB

1 牛頓—拉普森

對于非線性問題得到的非線性方程組，一般可表示為：

（）=0. （1）

寫成平衡方程形式：

（）=. （2）

式（2）中：（）為切向剛度矩陣；為位移矢量；為施加的載荷向量。

平衡迭代的過程表示為：

{i+1}={i}[i]-1（－i）. （3）

式（3）中：[i]為方程的Jacobian矩陣（即剛度矩陣）；{}為結(jié)點(diǎn)力矢量；{i}為內(nèi)力矢量；{－i}為不平衡力矢量。

牛頓—拉普森迭代步驟如下：①基于i時(shí)的結(jié)構(gòu)構(gòu)型計(jì)算i和i；②計(jì)算不平衡力矢量{－i}；③由i和不平衡力矢量{－i}計(jì)算位移增量；④更新位移向量△i；⑤重復(fù)②到④的過程直至計(jì)算收斂為止。

2 弧長法

弧長法的約束方程為：

{△}T{△}+2△2{}T{}=△2. （4）

式（4）中：為載荷比例系數(shù)；△為載荷增量；△為弧長半徑。

根據(jù)的不同，可以分為不同類型的弧長法：=1，球面弧長法；=0，柱面弧長法；等于當(dāng)前剛度參數(shù)值，橢球面弧長法。

弧長法的平衡迭代方程為：

其中：

以柱面弧長法為例推導(dǎo)，此時(shí)約束方程為：

{△}T{△}=△2. （9）

由位移增量關(guān)系可得：

為求得，把（6）式代入（10）式得：

把（12）代入（9）式可得關(guān)于一元二次方程：

式（12）中系數(shù)分別為：

對于收斂準(zhǔn)則，一般采用不平衡力準(zhǔn)則[3]。表達(dá)式為：

3 算例分析比較及討論

本文的算例如下：

受拉桿施加的力=20 kN，桿長=50 cm，截面積=2 cm2，材料的應(yīng)力應(yīng)變關(guān)系為：

=0（1－/0）. （17）

式（17）中：0=0.002，0=21 000 kN/cm2。

分別用牛頓—拉普森法和弧長法計(jì)算位移—載荷曲線，并與理論解做比較。理論解通過把應(yīng)力—應(yīng)變曲線轉(zhuǎn)化為載荷—位移曲線得到。過程如下：

=0（1－/0）；=/.

由以上兩式，并代入相關(guān)參數(shù)可得：

=﹣8 4002+840. （18）

下面利用牛頓—拉普森迭代法和弧長法（分別如圖1和圖2所示），借助MATLAB對算例進(jìn)行了實(shí)現(xiàn)，具體如圖3、圖4所示。

圖1 牛頓—拉普森迭代法

圖3 牛頓—拉普森法與理論解

圖4 弧長法與理論解

由圖1可知，用牛頓—拉普森迭代法計(jì)算得到的載荷—位移曲線與理論解的上升段是非常吻合的，但是由于在極值點(diǎn)附近出現(xiàn)發(fā)散，而不能繼續(xù)跟蹤曲線的下降段，但在跟蹤載荷—位移曲線的上升段時(shí)還是有效的。由圖2可知弧長法計(jì)算所得的載荷—位移曲線與理論解的全過程的逼近程度都是很好的。弧長法不僅可以跟蹤曲線的上升段，還可以跟蹤曲線中的下降段。因此，弧長法可以更好地跟蹤載荷—位移曲線的特點(diǎn)。

［1］Carrera E.A study on arc-length-type methods and their operation failures illustrated by a simple model［J］.Computers & Structures，1994，50（2）：217-229.

［2］Zhou Z，Murray D W. An incremental solution technique for unstable equilibrium paths of shell structures［J］. Computers & Structures，1995，55（5）：749-759.

［3］劉國，卓家壽，夏頌佑.求解非線性有限元方程的弧長法及在工程穩(wěn)定分析中的應(yīng)用［J］.巖土力學(xué)，1993（4）：57-67.

［4］李元齊，沈祖炎.弧長控制類方法使用中若干問題的探討與改進(jìn)［J］.計(jì)算力學(xué)學(xué)報(bào)，1998，15（4）：414-422.

［5］Bellini P X，Chulya A.An improved automatic incremental algorithm for efficient solution of nonlinear finite element equations［J］.Computers & Structures，1987，26（1）：99-110.

2095－6835（2018）24－0005－02

TM402

A

10.15913/j.cnki.kjycx.2018.24.005

魏鵬（1991—），男，研究方向?yàn)榉蔷€性。

〔編輯：嚴(yán)麗琴〕