SimPC博士：幾何非線性有限元基本原理及matlab編程

三、非線性系統(tǒng)求解算法

非線性系統(tǒng)求解難點(diǎn)在于，即在任何平衡構(gòu)型下，切線剛度矩陣的計(jì)算取決于結(jié)構(gòu)的變形幾何形狀和單元的內(nèi)力（見(jiàn)part1剛度矩陣推導(dǎo)）。 這些屬性是從節(jié)點(diǎn)位移獲得的，但節(jié)點(diǎn)位移是未知數(shù)。 因此，無(wú)法解析求解平衡方程組，需要運(yùn)用數(shù)值方法進(jìn)行處理。

求解非線性方程組追蹤平衡路徑可分為單步法和增量迭代法，本文重點(diǎn)介紹增量迭代法。 該方法將總載荷分解為一系列增量步，通過(guò)增量線性系統(tǒng)的直接或迭代求解獲得相應(yīng)載荷增量步的位移增量。 因此，通過(guò)每一步中的一個(gè)或幾個(gè)線性分析來(lái)近似施加載荷和節(jié)點(diǎn)位移之間的非線性關(guān)系。 求解非線性平衡系統(tǒng)的兩類增量方法：第一種是純?cè)隽炕騿尾椒椒?，其中使用單個(gè)剛度矩陣來(lái)表示每個(gè)分析步驟中的載荷位移關(guān)系。 第二種是增量迭代法，它執(zhí)行多次線性分析，在每一步迭代求解增量系統(tǒng)，以便在數(shù)值公差范圍內(nèi)收斂到平衡解。 本文主要針對(duì)增量迭代法進(jìn)行講解。

1、增量迭代法求解公式

為了追蹤平衡路徑，必須以增量方式多次求解非線性方程組。 非線性解通過(guò)每個(gè)增量的線性響應(yīng)來(lái)近似。 線性近似在線性化解和實(shí)際非線性解之間產(chǎn)生殘差。 該殘差通過(guò)通過(guò)迭代進(jìn)行收斂校正。 在后文中，增量步i-1對(duì)應(yīng)的構(gòu)型是最新獲得的平衡配置，而增量步 i 對(duì)應(yīng)的構(gòu)型是當(dāng)前未知的需要求解的構(gòu)型，符號(hào)Δ用來(lái)表示變量在每個(gè)增量步步的增量，符號(hào)δ表示每次迭代的增量。

切線剛度矩陣 K 在下式中定義，即內(nèi)力對(duì)上一步對(duì)應(yīng)的平衡構(gòu)型下的節(jié)點(diǎn)位移的導(dǎo)數(shù)。

?（26）

因此，內(nèi)力表式為

（27）

求解位移增量的增量系統(tǒng)表達(dá)式為，其中表示荷載比，表示第i增量步荷載比增量，與總荷載乘積即為荷載增量

（28）

梁?jiǎn)卧那芯€矩陣的推導(dǎo)在part1中通過(guò)虛擬位移原理、運(yùn)動(dòng)學(xué)描述（更新拉格朗日法或共旋法）和有限元方法來(lái)離散每個(gè)單元位移場(chǎng)來(lái)得到。 然而，由于內(nèi)力是位移的非線性函數(shù)，方程的線性化增量系統(tǒng)的解不滿足平衡。 即外力和內(nèi)力之間的不平衡，產(chǎn)生了殘余力矢量。 在每個(gè)增量中通常采用迭代過(guò)程，以通過(guò)消除殘余力來(lái)確保平衡。

以上是增量步的解釋，接下來(lái)對(duì)迭代步進(jìn)行講解。

考慮到結(jié)構(gòu)在第 i-1 步處于平衡狀態(tài)，希望在第 i 步達(dá)到平衡。 為了消除源自線性化增量的殘余力，通過(guò)在每個(gè)增量步驟中執(zhí)行Newton-Raphson之類的校正迭代來(lái)完成的，直到滿足收斂標(biāo)準(zhǔn)并建立新的平衡狀態(tài)。 在每次迭代中，由下標(biāo) j = 1,2,3... 表示，得到對(duì)應(yīng)迭代步的載荷比增量 δλij 和節(jié)點(diǎn)位移增量δUij。 這些迭代增量表示相應(yīng)步驟中載荷和位移值的校正。 因此，在第j次迭代之后，通過(guò)累加迭代增量來(lái)更新第i步的總增量，如下式所示：

（29）

基于上述迭代步對(duì)荷載和位移增量的修正，整個(gè)分析的載荷比和節(jié)點(diǎn)位移的總值更新如下：

（30）

在第 j 次迭代后得到殘余力矢量，由外部和內(nèi)部節(jié)點(diǎn)力的總值之間的差值給出：

（31）

基于上述殘余力矢量，可根據(jù)下式求得在第i步的第j次迭代中位移迭代增量

（32）

圖12

上述迭代算法最常用的是Newton-Raphson迭代算法，如圖12所示，標(biāo)準(zhǔn) Newton-Raphson 迭代算法的切線剛度矩陣在所有迭代步都會(huì)進(jìn)行更新，考慮大型系統(tǒng)時(shí)矩陣更新會(huì)耗費(fèi)較多計(jì)算資源，因此可使用修正Newton-Raphson 算法，切線矩陣僅在每個(gè)增量步驟的開(kāi)始進(jìn)行計(jì)算，并在所有后續(xù)迭代中保持不變，即， 對(duì)于 j > 1。 修正Newton-Raphson 算法具有較低的計(jì)算成本 每次迭代都比標(biāo)準(zhǔn)版本，但收斂通常較慢，因?yàn)樗ǔＰ枰诿總€(gè)增量步驟中進(jìn)行更多的迭代，但相較切向剛度矩陣更新耗費(fèi)的資源來(lái)說(shuō)還是值得的。 對(duì)于上述方程x求解來(lái)說(shuō)，因?yàn)榇嬖?img alt="poYBAGPbKqiAddh2AAABxf_93es931.png" src="https://file.elecfans.com/web2/M00/8C/39/poYBAGPbKqiAddh2AAABxf_93es931.png" />這個(gè)未知量，因此還需要一個(gè)附加方程來(lái)與之組成方程組進(jìn)行求解，這個(gè)附加方程的一般形式為

（33）

式中向量A和標(biāo)量 B 和 C 是常數(shù)，可以根據(jù)不同的求解方法采用不同的值。 上式x與組成的方程組寫(xiě)成矩陣形式為

（34）

由于上述等式左端的矩陣不再是對(duì)稱矩陣，所以在做大型計(jì)算時(shí)存儲(chǔ)和計(jì)算上效率都會(huì)明顯降低，為了克服這個(gè)問(wèn)題，Batoz & Dhatt (1979) 提出了一種技術(shù)來(lái)克服這個(gè)問(wèn)題，方法是將系統(tǒng)分解為兩個(gè)使用原始切線剛度矩陣的系統(tǒng)，因此原始系統(tǒng)的帶狀和對(duì)稱特性保持不變，具體如下所示

（35）

位移迭代增量的解是上述切線增量和殘差增量增量的線性組合：

（36）

其中迭代步j(luò)對(duì)應(yīng)的荷載比增量可通過(guò)下式求得

（37）

上述增量迭代法求解所用到的公式總結(jié)如下表所示

（1）增量迭代法求解步驟

上述增量迭代求解可分為兩個(gè)階段，分別為預(yù)測(cè)階段和校正階段。 預(yù)測(cè)階段相當(dāng)于第一次迭代，其余迭代屬于校正階段。

① 預(yù)測(cè)階段

在每個(gè)增量步，首先執(zhí)行預(yù)測(cè)階段的迭代。 目的是使用上一步得到的切線剛度矩陣進(jìn)行單次線性分析來(lái)計(jì)算預(yù)測(cè)解。 具體來(lái)講，這一階段，對(duì)于第i個(gè)增量步，首先要進(jìn)行初始切線剛度矩陣的計(jì)算，直接基于上一步結(jié)束時(shí)獲得的結(jié)果（即節(jié)點(diǎn)位移 和內(nèi)力對(duì)應(yīng)的剛度矩陣）。 然后根據(jù)方程（35）通過(guò)線性分析計(jì)算位移的切線增量 。 該階段的殘余位移增量為零，因?yàn)楹雎粤藖?lái)自上一步的殘余力。 位移的預(yù)測(cè)增量， ，可根據(jù)方程(36)獲得。 但僅是位移的切線增量乘以載荷比的增量。 載荷比增量是用約束方程（37）計(jì)算得到的，該約束方程定義了一個(gè)超曲面來(lái)限制增量迭代方法的校正解。 下圖13為單自由度系統(tǒng)的預(yù)測(cè)階段示意圖。

圖13

因此預(yù)測(cè)階段的求解的核心目的就是計(jì)算預(yù)測(cè)荷載比增量。 在某分析過(guò)程的第一個(gè)增量步的預(yù)測(cè)階段荷載比增量須人為指定，一般為最大荷載的10%~20%。 在接下來(lái)的迭代過(guò)程中，荷載比增量則由約束方程（37）所定義的迭代算法進(jìn)行計(jì)算，后文荷載比增量的計(jì)算方法本質(zhì)就是對(duì)約束方程（37）的定義。

預(yù)測(cè)階段荷載比增量對(duì)求解有很大的影響，直接與收斂性相關(guān)。 該階段，荷載比增量過(guò)大可能導(dǎo)致收斂問(wèn)題，荷載比增量過(guò)小耗費(fèi)計(jì)算資源，精度的提高也有限。 因此使得程序能夠自動(dòng)根據(jù)非線性程度調(diào)整預(yù)測(cè)階段的荷載比增量是一個(gè)優(yōu)秀的非線性算法所應(yīng)具備的功能。 即當(dāng)結(jié)構(gòu)響應(yīng)基本程線性時(shí)，提高增量，當(dāng)非線性較大時(shí)，減小增量。 而且當(dāng)求解至荷載極限點(diǎn)達(dá)到一個(gè)即將失穩(wěn)的平衡態(tài)時(shí)，所追蹤的位移增量對(duì)應(yīng)的荷載增量必須是負(fù)值才符合常理，因此算法還應(yīng)能夠判斷增量正負(fù)的選擇。

a、預(yù)測(cè)階段荷載比增量計(jì)算方法

荷載比增量的計(jì)算取決于約束方程（37）的定義，本節(jié)計(jì)算方法本質(zhì)即約束方程的定義。 兩種計(jì)算預(yù)測(cè)階段荷載比增量的方法，一種基于迭代次數(shù)的計(jì)算方法，另一種是根據(jù)系統(tǒng)剛度的計(jì)算方法。 本文重點(diǎn)針對(duì)第一種方法進(jìn)行介紹。 基于迭代次數(shù)的計(jì)算方法計(jì)算預(yù)測(cè)階段荷載比增量的方法又可以分為三類：

荷載增量法；

外力功增量法；

弧長(zhǎng)增量法。

我們重點(diǎn)介紹第三種弧長(zhǎng)增量法，因?yàn)榛￠L(zhǎng)法的優(yōu)勢(shì)在于在于可以追蹤平衡路徑，準(zhǔn)確捕捉snap-through和snap-back現(xiàn)象，如圖14所示。 為了數(shù)值算法上增量步大小一致性，校正階段也采用同類型的弧長(zhǎng)法。 弧長(zhǎng)法公式具體推導(dǎo)過(guò)程大家可參考相關(guān)文獻(xiàn)，這里只列出用圓柱弧長(zhǎng)法和球形弧長(zhǎng)法計(jì)算預(yù)測(cè)階段荷載比增量的最終公式

（38）

（39）

其中J為迭代調(diào)整因子，表達(dá)式如下

（40）

式中，Ni 和 Ni-1 分別是當(dāng)前增量步所需的迭代次數(shù)，以及前一增量步中實(shí)現(xiàn)收斂的迭代次數(shù)。 指數(shù)變量 η 通常在 0.5 到 1.0 的范圍內(nèi)，但通常采用較低的值。

圖14

b、預(yù)測(cè)階段荷載比增量符號(hào)計(jì)算方法

上文提到預(yù)測(cè)階段的求解的核心目的除了需要確定荷載比增量大小外，另一個(gè)重要的工作就是完成荷載比增量正負(fù)符號(hào)的確定。 確定荷載比增量符號(hào)最常見(jiàn)方法的是基于剛度參數(shù)的行為進(jìn)行確定，這些剛度參數(shù)常見(jiàn)的有 CSP（Current Stiffness Parameter）和 GSP（Generalized Stiffness Parameter）。 因?yàn)閰?shù)CSP會(huì)在位移極限點(diǎn)附近出現(xiàn)一些數(shù)值不穩(wěn)定性，所以本文主要介紹更為通用的GSP參數(shù)。 GSP的具體表達(dá)式為

?（41）

可以看出，GSP參數(shù)的分子是一個(gè)正數(shù)，分母由當(dāng)前和前一步位移向量的標(biāo)量積給出，該標(biāo)量積的符號(hào)可以反應(yīng)前一步和當(dāng)前步的增量是否在同一個(gè)“方向”上，如果同向則為正，如果反向則為負(fù)，如圖 15所示。 當(dāng)兩個(gè)增量步之間存在荷載極限點(diǎn)時(shí)，兩個(gè)位移向量的方向是不同的，因此GSP<0。 因此，每次 GSP為負(fù)值時(shí)，必須反轉(zhuǎn)預(yù)測(cè)階段的荷載比增量的符號(hào)。 增量符號(hào)可根據(jù)下式確定，其中初始增量步假定增量符號(hào)為正值。

（42）

圖15

（2）迭代矯正階段

增量迭代方法的校正階段旨在通過(guò)迭代循環(huán)消除由預(yù)測(cè)階段產(chǎn)生的殘余力來(lái)恢復(fù)結(jié)構(gòu)平衡。 該迭代循環(huán)從更新總荷載比和總節(jié)點(diǎn)位移 開(kāi)始，將預(yù)測(cè)增量（和 ）添加到上一步的結(jié)果（ 和 ）。 隨著幾何構(gòu)型的更新，相應(yīng)的內(nèi)力 被計(jì)算出來(lái)，殘余力可以通過(guò)外部和內(nèi)部節(jié)點(diǎn)力之間的差異來(lái)獲得。

檢查迭代是否收斂的方法，最常見(jiàn)的是基于殘余力、殘余位移或這些殘余結(jié)果產(chǎn)生的功。 本文采用的收斂標(biāo)準(zhǔn)是基于力的檢查，如下式所示

（43）

式中，分子分母分別為殘余力矢量和參考載荷矢量的歐幾里得范數(shù)，，該比率必須低于給定的容差ε，通常在 10-5 到 10-3 的數(shù)量級(jí)。 如果滿足收斂，則預(yù)測(cè)解被認(rèn)為平衡狀態(tài)，可致性下一增量步，否則開(kāi)始第一次校正迭代。

每次校正迭代的第一個(gè)過(guò)程是計(jì)算切線剛度矩陣，要根據(jù)最新構(gòu)型的節(jié)點(diǎn)位移和內(nèi)力來(lái)對(duì)剛度矩陣進(jìn)行更新。 如果采用修正Newton-Raphson 迭代方案，則跳過(guò)此步驟，并使用在預(yù)測(cè)階段的切線矩陣。 如方程(35)所示，位移的切線增量和殘余增量分別用參考載荷矢量和最后獲得的殘余力矢量計(jì)算。 然后，根據(jù)相應(yīng)的約束方程（37）計(jì)算載荷比的迭代增量。 最后，用方程(36)得到位移的迭代增量。。

載荷比和位移的迭代增量取決于約束方程（37）定義的超曲面。 如果執(zhí)行的迭代不僅涉及位移修正，還涉及載荷比的修正，則稱為連續(xù)法，因?yàn)樗梢愿櫝鰳O限點(diǎn)的平衡路徑，例如上小節(jié)提到的弧長(zhǎng)法。 在這種情況下，控制修正解的約束面在一個(gè)或多個(gè)點(diǎn)處與平衡路徑相交。 獲得修正解后，接下來(lái)的程序與檢查預(yù)測(cè)解的收斂性相同：更新載荷比和節(jié)點(diǎn)位移的總值，計(jì)算外力、內(nèi)力和殘余力，最后檢查當(dāng)前迭代解的收斂性。 圖16 給出了所描述過(guò)程的示意圖。

圖16

與預(yù)測(cè)階段的荷載比計(jì)算方法相對(duì)應(yīng)，迭代階段的荷載比增量計(jì)算方法同樣有荷載控制法、外力功控制法和弧長(zhǎng)控制法（還有其他多種方法，就不一一列舉）。 傳統(tǒng)的牛頓拉普森迭代法本質(zhì)就是荷載增量控制法，在每個(gè)增量步中使用固定量的荷載比增量即預(yù)測(cè)階段的荷載增量，并在每次迭代后保持不變。 執(zhí)行校正迭代僅通過(guò)位移校正來(lái)滿足平衡要求，如下圖17所示。

圖17

當(dāng)試圖解決荷載極限點(diǎn)問(wèn)題時(shí)，這種方法的存在明顯的缺陷。 一旦在預(yù)測(cè)階段定義了固定荷載比增量，如果在增量中出現(xiàn)極限點(diǎn)，就無(wú)法修改荷載矢量。 盡管可減小的載荷比增量使其緩慢地接近極限點(diǎn)，但由此產(chǎn)生的剛度矩陣接近奇異的性質(zhì)使得難以追蹤結(jié)構(gòu)的后極限狀態(tài)響應(yīng)。 圖 18 顯示了使用荷載增量控制法跟蹤具有snap-though行為的系統(tǒng)的平衡路徑時(shí)的典型結(jié)果。

圖18

本文重點(diǎn)介紹適應(yīng)度較強(qiáng)的弧長(zhǎng)控制法，弧長(zhǎng)法的優(yōu)勢(shì)在于在于可以追蹤平衡路徑，準(zhǔn)確捕捉snap-through和snap-back現(xiàn)象。 由于弧長(zhǎng)法也分多種，如線性弧長(zhǎng)法和恒定弧長(zhǎng)法。 這里僅介紹恒定弧長(zhǎng)法，即由位移和載荷增量的范數(shù)定義的弧長(zhǎng)增量 ΔL 在整

個(gè)迭代過(guò)程中保持不變。 推導(dǎo)過(guò)程省略，這里直接給出恒定圓柱弧長(zhǎng)法計(jì)算迭代步荷載比增量的公式

（44）

求解上式可得迭代步荷載比增量的顯式表達(dá)式

（45）

作者：SimPC博士

審核編輯：湯梓紅