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

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

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

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

1、增量迭代法求解公式

為了追蹤平衡路徑，必須以增量方式多次求解非線(xiàn)性方程組。 非線(xiàn)性解通過(guò)每個(gè)增量的線(xiàn)性響應(yīng)來(lái)近似。 線(xiàn)性近似在線(xiàn)性化解和實(shí)際非線(xiàn)性解之間產(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)δ表示每次迭代的增量。

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

以上是增量步的解釋?zhuān)酉聛?lái)對(duì)迭代步進(jìn)行講解。

考慮到結(jié)構(gòu)在第 i-1 步處于平衡狀態(tài)，希望在第 i 步達(dá)到平衡。 為了消除源自線(xiàn)性化增量的殘余力，通過(guò)在每個(gè)增量步驟中執(zhí)行Newton-Raphson之類(lèi)的校正迭代來(lái)完成的，直到滿(mǎn)足收斂標(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 迭代算法的切線(xiàn)剛度矩陣在所有迭代步都會(huì)進(jìn)行更新，考慮大型系統(tǒng)時(shí)矩陣更新會(huì)耗費(fèi)較多計(jì)算資源，因此可使用修正Newton-Raphson 算法，切線(xiàn)矩陣僅在每個(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ì)稱(chēng)矩陣，所以在做大型計(jì)算時(shí)存儲(chǔ)和計(jì)算上效率都會(huì)明顯降低，為了克服這個(gè)問(wèn)題，Batoz & Dhatt (1979) 提出了一種技術(shù)來(lái)克服這個(gè)問(wèn)題，方法是將系統(tǒng)分解為兩個(gè)使用原始切線(xiàn)剛度矩陣的系統(tǒng)，因此原始系統(tǒng)的帶狀和對(duì)稱(chēng)特性保持不變，具體如下所示

（35）

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

（36）

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

（37）

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

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

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

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

在每個(gè)增量步，首先執(zhí)行預(yù)測(cè)階段的迭代。 目的是使用上一步得到的切線(xiàn)剛度矩陣進(jìn)行單次線(xiàn)性分析來(lái)計(jì)算預(yù)測(cè)解。 具體來(lái)講，這一階段，對(duì)于第i個(gè)增量步，首先要進(jìn)行初始切線(xiàn)剛度矩陣的計(jì)算，直接基于上一步結(jié)束時(shí)獲得的結(jié)果（即節(jié)點(diǎn)位移 和內(nèi)力對(duì)應(yīng)的剛度矩陣）。 然后根據(jù)方程（35）通過(guò)線(xiàn)性分析計(jì)算位移的切線(xiàn)增量 。 該階段的殘余位移增量為零，因?yàn)楹雎粤藖?lái)自上一步的殘余力。 位移的預(yù)測(cè)增量， ，可根據(jù)方程(36)獲得。 但僅是位移的切線(xiàn)增量乘以載荷比的增量。 載荷比增量是用約束方程（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ù)非線(xiàn)性程度調(diào)整預(yù)測(cè)階段的荷載比增量是一個(gè)優(yōu)秀的非線(xiàn)性算法所應(yīng)具備的功能。 即當(dāng)結(jié)構(gòu)響應(yīng)基本程線(xiàn)性時(shí)，提高增量，當(dāng)非線(xiàn)性較大時(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è)階段荷載比增量的方法又可以分為三類(lèi)：

荷載增量法；

外力功增量法；

弧長(zhǎng)增量法。

我們重點(diǎn)介紹第三種弧長(zhǎng)增量法，因?yàn)榛￠L(zhǎng)法的優(yōu)勢(shì)在于在于可以追蹤平衡路徑，準(zhǔn)確捕捉snap-through和snap-back現(xiàn)象，如圖14所示。 為了數(shù)值算法上增量步大小一致性，校正階段也采用同類(lèi)型的弧長(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í)。 如果滿(mǎn)足收斂，則預(yù)測(cè)解被認(rèn)為平衡狀態(tài)，可致性下一增量步，否則開(kāi)始第一次校正迭代。

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

載荷比和位移的迭代增量取決于約束方程（37）定義的超曲面。 如果執(zhí)行的迭代不僅涉及位移修正，還涉及載荷比的修正，則稱(chēng)為連續(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)滿(mǎn)足平衡要求，如下圖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)法也分多種，如線(xiàn)性弧長(zhǎng)法和恒定弧長(zhǎng)法。 這里僅介紹恒定弧長(zhǎng)法，即由位移和載荷增量的范數(shù)定義的弧長(zhǎng)增量 ΔL 在整

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

（44）

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

（45）

作者：SimPC博士

審核編輯：湯梓紅