非線性方程組可以表示為：
$\pmb \psi(\pmb a)=\pmb P(\pmb a)-\pmb Q=0　　　?。ǎ保? mathimg=$
在位移為基本未知量的有限元分析中， $\pmb a$ 是結(jié)點(diǎn)位移向量， $\pmb Q$ 是結(jié)點(diǎn)載荷向量。對(duì)于線性方程組 $\pmb K \pmb a=\pmb Q$ ，由于 $\pmb K$ 是常數(shù)矩陣，可以沒(méi)有困難直接求解。對(duì)于非線性方程組需要迭代求解。

1 直接迭代法

只適用于與變形歷史無(wú)關(guān)的非線性問(wèn)題，例如非線性彈性問(wèn)題，利用形變理論分析的彈塑性問(wèn)題，穩(wěn)態(tài)蠕變問(wèn)題等。對(duì)于依賴于變形歷史的非線性問(wèn)題，應(yīng)力需要由應(yīng)變所經(jīng)歷的路徑?jīng)Q定，直接迭代法不適用。例如加載路徑不斷變化或涉及卸載和反復(fù)加載等彈塑性問(wèn)題。此時(shí)要利用增量理論。
可以指出，當(dāng) $P(a)-a$ 是凸曲線（如圖所示），其中 $a$ 是標(biāo)量，即系統(tǒng)為單自由度情形，通常解收斂。

3.PNG

1.1直接迭代法的求解步驟

(1)方程可以改寫為 $\pmb K(\pmb a) \pmb a=\pmb Q$
選擇一個(gè)初始的試探解 $\pmb a＝\pmb a^{(0)}$ 代入上式可以得到初始的 $\pmb K^{(0)}$ ，接著可以得到第一步迭代的位移解 $\pmb a^{(1)}=\pmb (K^{(0)})^{-1} \pmb Q$ ，由此類推得到第 $n$ 次迭代后的近似解 $\pmb a^{(n)}＝\pmb (K^{(n-1)})^{-1} \pmb Q$ ，一直到誤差的某種范數(shù)小于某個(gè)規(guī)定的容許小量 $er$ ，即 $||e||=||\pmb a^{(n)}-\pmb a^{(n-1)}||≤er$ ，上述迭代終止。

1.2常系數(shù)矩陣的直接迭代法

為避免每次迭代對(duì)剛度矩陣 $\pmb K$ 進(jìn)行求逆計(jì)算，一般可以將剛度矩陣設(shè)定為常數(shù)進(jìn)行迭代過(guò)程，單自由度系統(tǒng)下的常剛度直接迭代法如圖所示

4.PNG

利用下式求

a^{(1)}

的修正量

\Delta a^{(1)}

\Delta a^{(1)}=\pmb (K^{(0)})^{-1}(\pmb Q-(\pmb K^{(1)}a^{(1)}) )　　　?。ǎ玻? mathimg=

由此可以得到

a^{(2)} = a^{(1)}+\Delta a^{(1)}

同理可以得到

\Delta a^{(n-1)} = (\pmb K^{(0)})^{-1}(\pmb Q-(\pmb K^{(n-1)} a^{(n-1)}) )

a^{(n)} = a^{(n-1)}+\Delta a^{(n-1)}

非線性材料重新形成

\pmb K^{(n-1)}

的工作量遠(yuǎn)小于對(duì)其進(jìn)行分解求逆的工作量（用于迭代的剛度矩陣和用于計(jì)算內(nèi)力的剛度矩陣注意分開(kāi)）

2 Newton-Raphson方法

一般情況下，（１）不能被精確地滿足，即 $\pmb \psi(a^{(n)}) \neq 0$ ，為了得到進(jìn)一步的近似解 $a^{(n+1)}$ （假設(shè) $a^{(n+1)}$ 為某單自由度方向上的位移），將 $\pmb \psi(a^{(n+1)})$ 表示成在 $a^{(n)}$ 附近的僅保留線性項(xiàng)的Taylor展開(kāi)式
$\pmb \psi(a^{(n+1)})=\pmb \psi(a^{(n)})+(\frac{d\pmb\psi}{da})_{n}\Delta a^{(n)}=0　?。ǎ常? mathimg=$
且有 $a^{(n+1)} = a^{(n)}+\Delta a^{(n)}$
式中 $d\psi/da$ 是切線矩陣(tangent matrix)，即
$\frac{d\pmb\psi}{da}=\frac{d\pmb P}{da}=\pmb K_T(a)$
$\Delta a^{(n)} = (\pmb K_T^{(n)})^{-1}(\pmb Q-\pmb P^{(n)})$
$\pmb K_T^{(n)}=\pmb K_T(a^{(n)})$ $\pmb P^{(n)}=\pmb P(a^{(n)})$
其迭代過(guò)程如圖所示

5.PNG

為了克服N-R方法對(duì)于每次迭代需要形成并求逆一個(gè)新的切線矩陣所帶來(lái)的麻煩，可以仿照直接迭代法采用一種修正的N-R方法。即

\pmb K_T^{(n)}=\pmb K_T^{(0)})

如圖所示

6.PNG

3 增量法

核心思想是將載荷分解成若干增量步，即 $\pmb Q_0,\pmb Q_1,\pmb Q_2,\pmb Q_3,\pmb Q_4......$ ，相應(yīng)的位移也分為同樣的步數(shù)，即 $\pmb a_0,\pmb a_1,\pmb a_2,\pmb a_3,\pmb a_4......$ ，每?jī)刹街g的增長(zhǎng)量稱之為增量。增量解法的一般做法是假設(shè)第 $m$ 步的載荷 $\pmb Q_m$ 和相應(yīng)的位移 $a_m$ 為已知，然后將載荷增加為 $\pmb Q_{m+1}(=\pmb Q_{m}+\Delta \pmb Q_{m})$ ，如果每一步的載荷增量 $\Delta \pmb Q_{m}$ 足夠小，則解的收斂性可以保證。
使用增量法的（１）式改寫為如下形式（位移均以單自由度為例）
$\pmb \psi(a)=\pmb P(a)- \lambda \pmb Q_0=0$ ,其中 $\lambda$ 是用來(lái)表示載荷變化的參數(shù)，將上式對(duì) $\lambda$ 求導(dǎo)，可以得到
$\frac{d\pmb P}{da}\frac{da}{d\lambda}-\pmb Q_0=\pmb K_T \frac{da}{d\lambda}=0$
$\frac{da}{d\lambda}=\pmb K_T^{-1}(a)\pmb Q_0　　　?。ǎ矗? mathimg=$
其中 $\pmb K_T$ 為定義的切線剛度矩陣。
（４）式是一典型的常微分方程組問(wèn)題，下面介紹的是有限元中對(duì)這類方程組的求解方法
1.歐拉方法(單自由度為例)
$a_{m+1}-a_m=\Delta a_m=\pmb K_T^{-1}(a_m)\pmb Q_0 \Delta \lambda_m= (\pmb K_T)^{-1}_{m} \Delta \pmb Q_m　　　?。ǎ担? mathimg=$
顯然為了滿足精度要求， $\Delta \lambda_m$ 必須是足夠小的量。使用Runge-Kutte方法可以改善精度
此方法會(huì)導(dǎo)致解的漂移（與真實(shí)曲線上的解產(chǎn)生較大的誤差）
為了克服每一增量步解的誤差可能導(dǎo)致的解的漂移，可以將（５）式改寫為
$a_{m+1}-a_m=\Delta a_m=\pmb K_T^{-1}(a_m)(\pmb Q_{m+1}-\pmb P(a_m))　　　（６）$
這里解釋一下 $(\pmb Q_{m+1}-\pmb P(a_m))=(\pmb Q_{m+1}-\pmb Q_{m}+\pmb Q_{m}-\pmb P(a_m))$
此方法稱為考慮平衡矯正的歐拉增量方法，即將前增量步的載荷和內(nèi)力的不平衡量合并到當(dāng)前增量步求解，一定程度上避免了解的漂移。

7.PNG

2.N-R方法
為了改進(jìn)歐拉法的精度，現(xiàn)在更多采用N-R方法，如果采用N-R方法，是在每一增量步內(nèi)進(jìn)行迭代
則對(duì)于 $\lambda$ 的 $m+1$ 次增量步的第 $n+1$ 次迭代可以表示為
$\psi_{m+1}^{(n+1)}=\pmb P(a_{m+1}^{(n+1)})-\pmb Q_{m+1}=\pmb P(a_{m+1}^{(n)})-\pmb Q_{m+1}+(\pmb K_T^n)_{m+1}\Delta a_m^n=0　　?。ǎ叮? mathimg=$
表示成前述的N-R形式為 $\Delta a_m^n=(\pmb K_T^n)_{m+1}^{-1}(\pmb Q_{m+1}-\pmb P(a_{m+1}^{(n+1)}))$
采用mN-R迭代
$(\pmb K_T^n)_{m+1}=(\pmb K_T^0)_{m+1}=\pmb K_T(a_m)$

8.PNG

9.PNG

4 加速收斂的方法(Aitken加速法)

利用mN-R方法求解非線性方程組時(shí)，可以避免每次迭代重新形成和求逆切線矩陣，但是降低了收斂速度，特別是 $P-a$ 曲線突然趨于平坦的情況
有Aitken加速的迭代和無(wú)Aitken加速的迭代如圖所示

10.PNG

11.PNG

核心原理是將通過(guò)初始切線剛度

\pmb K_T^0

得到的第m+1個(gè)增量步的第n個(gè)迭代解

\Delta a_{m}^{n}

利用局部割線剛度轉(zhuǎn)換出新的迭代解

\Delta \widetilde a_{m}^{n}

，以

n=0

說(shuō)明此轉(zhuǎn)換過(guò)程

K_s \Delta a_{m}^{0}=(K_T)_m(\Delta a_{m}^{0}-\Delta a_{m}^{1})

可以改寫為

\frac{(K_T)_m}{K_s}=\frac{\Delta a_{m}^{0}}{(\Delta a_{m}^{0}-\Delta a_{m}^{1})}=\alpha^{(1)}

K_s \Delta \widetilde a_{m}^{1}=(K_T)_m\Delta a_{m}^{1}

\Delta \widetilde a_{m}^{1}=\frac{(K_T)_m}{K_s}\Delta a_{m}^{1}=\alpha^{(1)}\Delta a_{m}^{1}

\alpha^{(1)}

被稱作加速因子
Aitken加速收斂方法是每隔一次迭代進(jìn)行一次加速（為了保證精度的同時(shí)提高收斂速度）。
推廣到

N

個(gè)自由度系統(tǒng)，Aitken方法可以表示為

\pmb a_{m+1}^{(n+1)}=\pmb a_{m+1}^{n}+\pmb \alpha^{(n)}\Delta \pmb a_{m}^{n}

其中

\pmb \alpha^{(n)}

是對(duì)角矩陣，其元素

\alpha_i^{(n)}=\left\{\begin{matrix}1 \\ \frac{\Delta a_{i,m}^{n-1}}{\Delta a_{i,m}^{n-1}-\Delta a_{i,m}^{n}}\end{matrix}\right.

5載荷增量步長(zhǎng)的自動(dòng)選擇

在研究載荷增量步長(zhǎng)自動(dòng)選擇的方法時(shí)，首先是假設(shè)載荷的分布模式是給定的，變化的只是他的幅值，在此情況下，外載荷可以表示為：
$^t \overline {\pmb F} ={^t}P \overline{\pmb F_0}$
等效結(jié)點(diǎn)載荷向量可以表示為：
$^t \pmb Q_l ={^t}P \pmb Q_0$
$P$ 是載荷幅值，載荷分布實(shí)際是 $P(t)$ 的分布
$^{t+\Delta t}P = {^t}P+\Delta P$
求解上述載荷分布方程依然是一個(gè)N-R迭代過(guò)程，設(shè) $^{t+\Delta t}P^{(0)}={^t}P$
$^{t+\Delta t}P^{(n+1)} = ^{t+\Delta t}P^{(n)}+\Delta P^{(n)}$
$^{t+\Delta t}P^{(n+1)}$ 是經(jīng)過(guò)n+1次迭代修正后得到的 $^{t+\Delta t}P$ 數(shù)值
下面是幾種常用的自動(dòng)選擇載荷步長(zhǎng)的方法

5.1規(guī)定”本步剛度參數(shù)“的變化量以控制載荷增量

此方法對(duì)于計(jì)算結(jié)構(gòu)的極限載荷很有效，利用本步剛度參數(shù)可以使步長(zhǎng)調(diào)整的比較合理，并可以減少總的增量步數(shù)，適合于計(jì)算由理想彈塑性材料組成結(jié)構(gòu)的極限載荷情況。

5.1.1 ”本步剛度參數(shù)“的概念

第i增量步結(jié)構(gòu)的總體剛度可用下式度量
$S_p^{(i)*}=\frac{\Delta \pmb Q^{(i)T}\Delta \pmb Q^{(i)}}{\Delta \pmb a^{(i)T} \Delta \pmb Q^{(i)}}$
初始結(jié)構(gòu)總體剛度的度量是
$S_p^{(0)*}=\frac{\pmb Q^{(e)T}\pmb Q^{(e)}}{\pmb a_e^{T} \pmb Q_e}$
其中 $\pmb Q_e$ 和 $\pmb a_e$ 是載荷向量和按彈性分析得到的位移向量。

5.1.2 增量步長(zhǎng)的自動(dòng)選擇

(1) 利用”本步剛度參數(shù)“的規(guī)定變化量自動(dòng)選擇增量步長(zhǎng)。
$S_p^{(i)}= \frac{S_p^{(i)*}}{S_p^{(0)*}}$
$S_p^{(i)}$ 稱為第i步的”本步剛度參數(shù)“，它代表結(jié)構(gòu)本身的剛度性質(zhì)，與載荷增量大小無(wú)關(guān)，當(dāng)結(jié)構(gòu)處于完全彈性時(shí)， $S_p^{(i)}=1$ ，隨著塑性區(qū)擴(kuò)大，結(jié)構(gòu)變軟， $S_p^{(i)}$ 逐漸減小，到達(dá)極限載荷時(shí)， $S_p^{(i)}=0$
對(duì)于比例加載的情況：
$\pmb Q_e=p_e\pmb Q_0$ 　　　　　 $\Delta \pmb Q^{(i)}=\Delta p_i\pmb Q_0$
$p_e$ 是極限載荷參數(shù)， $\pmb Q_0$ 是 $p=1$ 時(shí)的結(jié)點(diǎn)載荷向量
利用本步剛度參數(shù)可以使步長(zhǎng)調(diào)整得比較合理，并可以減少總的增量步
(2) 增量步長(zhǎng)的自動(dòng)選擇
利用(1)得到的”本步剛度參數(shù)“的規(guī)定變化量自動(dòng)選擇增量步長(zhǎng)，具體的算法步驟如下
1、求彈性極限載荷參數(shù) $p_e$
先施加任意載荷 $p \pmb Q_0$ ，假定結(jié)構(gòu)為完全彈性求解，求出結(jié)構(gòu)的最大等效應(yīng)力 $\overline \sigma_{max}$ ，則有
$p_e=\frac{\sigma_{s0}}{\overline \sigma_{max}}$
$\sigma_{s0}$ 是材料的初始屈服應(yīng)力
2、給定第一步載荷參數(shù)增量 $\Delta p_1=p_e/N,$ $N$ 的值可以事先給定，用 $p_1=p_e+\Delta p_1$ 求解第一增量步
3、給定第二及以后各增量步的剛度參數(shù)變化的預(yù)測(cè)值 $\Delta \widetilde S_p$ ，其大小決定步長(zhǎng)的大小，例如可在0.05~0.2之間選擇，并給定剛度的最小允許值 $S_p^{min}(S_p^{min}\neq 0)$ ，則每步載荷參數(shù)增量為
$\Delta p_i=\Delta p_{i-1} \frac{min(\Delta \widetilde S_p,S_p^{(i-1)}-S_p^{min})}{\lvert S_p^{(i-2)}-S_p^{(i-1)}\rvert}　　　　（i=2,3,......）$
然后用 $p_i=p_{i-1}+\Delta p_i$ 求解第 $i$ 增量步
然后可以通過(guò)(1)中的公式計(jì)算出”本步剛度參數(shù)“

12.PNG

從圖中可見(jiàn)

\Delta p_i

是在給定預(yù)測(cè)值

\Delta \widetilde S_p

情況下，利用

\Delta p_{i-1}/{\lvert S_p^{(i-2)}-S_p^{(i-1)}\rvert}

線性外推得到的

(y=kx)

，為了得到更好的適應(yīng)復(fù)雜的曲線，可以使用二次外推。

5.2規(guī)定某個(gè)結(jié)點(diǎn)的位移增量以確定載荷增量

以 $mN-R$ 迭代為例，它可以表示為
$^\tau \pmb K_{eq}\Delta \pmb a^{(n)}=\Delta \pmb Q^{(n)}　　　(n=0,1,2,...)$
在載荷增量步長(zhǎng)自動(dòng)控制的求解方法中， $\Delta \pmb Q^{(n)}$ 可以表示成
$\Delta \pmb Q^{(n)}=^{t+\Delta t} \pmb Q_l^{(n+1)}-^{t+\Delta t} \pmb Q_i^{(n)}　　　?。ǎ罚? mathimg=$
其中 $^{t+\Delta t} \pmb Q_l^{(n+1)}=^{t+\Delta t} \pmb Q_l^{(n)}+\Delta \pmb Q_l^{(n)}=(^{t+\Delta t}p^{(n)}+\Delta p^{(n)})\pmb Q_0$
$^{t+\Delta t} \pmb Q_i^{(n)}=\sum_{e} \int_{V_e}\pmb B^T {^{t+\Delta t}\pmb \sigma^{(n)}}dV$
以上兩式中:
$^{t+\Delta t}p^{(0)}={^t}p$ 　　　　　　　 $^{t+\Delta t}\pmb \sigma^{(0)}={^t}\pmb \sigma$
因此（７）可以改寫為 $\Delta \pmb Q^{(n)}=\Delta \pmb Q_u^{(n)}+\Delta \pmb Q_l^{(n)}=\Delta \pmb Q_u^{(n)}+\Delta p^{(n)}\pmb Q_0$
其中 $\Delta \pmb Q_u^{(n)}={^{t+\Delta t}}p^{(n)}\pmb Q_0-{^{t+\Delta t}} \pmb Q_i^{(n)}$
$\Delta p^{(n)}$ 是在第 $n$ 次迭代中由某個(gè)規(guī)定的約束條件來(lái)確定的載荷因子增量 $\Delta p$ 的第 $n$ 次修正量。在現(xiàn)在的方法中，這個(gè)條件就是某個(gè)結(jié)點(diǎn)的位移增量的大小，例如規(guī)定 $\Delta \pmb a$ 中的最大分量 $\Delta a_g$ 是給定的，此條件可以表示為
$\Delta a_g^{(n)}=\pmb b^T \Delta \pmb a^{(n)}=\left\{\begin{matrix}\Delta l（n=0）\\0 （n=1,2,...）\end{matrix}\right.（８）$
其中 $\Delta l$ 是 $\Delta a_g$ 的規(guī)定值， $\pmb b$ 是除第g個(gè)元素為1，其余元素為零的向量，具體迭代算法步驟如下：
(1)計(jì)算對(duì)于節(jié)點(diǎn)載荷模式 $\pmb Q_0$ 的位移模式 $\pmb a_0$ ，即
$\pmb a_0={^\tau} \pmb K_{eq}^{-1}\pmb Q_0$
(2)計(jì)算對(duì)于不平衡結(jié)點(diǎn)力向量 $\Delta \pmb Q_u^{(n)}$ 的位移增量修正值 $\Delta \pmb a_u^{(n)}$ 和 $n$ 次迭代后位移增量修正值的全量 $\Delta \pmb a^{(n)}$ ，即
$\Delta \pmb a_u^{(n)}={^\tau} \pmb K_{eq}^{-1}\Delta \pmb Q_u^{(n)}$
$\Delta \pmb a^{(n)}=\Delta \pmb a_u^{(n)}+\Delta p^{(n)}\pmb a_0$
其中 $\Delta p^{(n)}$ 是待定載荷因子增量的修正值
(3)利用條件(8)確定 $\Delta p^{(n)}$ ，由于
$\Delta a_g^{(n)}=\pmb b^T \Delta \pmb a^{(n)}=\pmb b^T \Delta \pmb a_u^{(n)}+\Delta p^{(n)}\pmb b^T\pmb a_0=\left\{\begin{matrix}\Delta l（n=0）\\0 （n=1,2,...）\end{matrix}\right.（８）$
從上式可以得到
$\Delta p^{(n)}=\frac{\Delta a_g^{(n)}-\pmb b^T \Delta \pmb a_u^{(n)}}{\pmb b^T\pmb a_0}$
這樣就確定了第 $n$ 次迭代后的載荷因子增量
載荷因子增量求出之后，可以知道每一步的外載荷，從而使用 $mN-R$ 法迭代得到位移，進(jìn)而求出應(yīng)變，應(yīng)力和當(dāng)前增量步的內(nèi)力等物理量。并檢驗(yàn)收斂準(zhǔn)則是否滿足，如未滿足，回到步驟二進(jìn)行新的一次迭代，直至收斂準(zhǔn)則滿足為止。關(guān)于每一個(gè)增量步某個(gè)指定結(jié)點(diǎn)位移增量 $\Delta l$ 本身的選擇，通常的方法是第1個(gè)增量步可以由某個(gè)給定的載荷因子增量 $\Delta p_1$ (例如令 $\Delta p_1=p_e/N$ ，其中p_e是彈性極限載荷因子， $N$ 可取5~10)，通過(guò)求解得到 $\Delta l_1$ ，以后各增量步的 $\Delta l_i$ 可由下式確定，即
$\Delta l_i=\sqrt{\frac{r_0}{r_{i-1}}}\Delta l_{i-1}$
其中 $r_0=4,5,6$