基于修正牛頓-拉普森迭代的數(shù)值流形法

張亞軍，莫思陽，張友良

(海南大學(xué)土木建筑工程學(xué)院，海南 ?？?570228)

1 引言

數(shù)值流形法[1，2]是石根華博士以數(shù)學(xué)流形為基礎(chǔ)，利用有限覆蓋技術(shù)建立起來的一種新型數(shù)值方法。該方法可以統(tǒng)一處理連續(xù)問題和非連續(xù)問題，分析小變形開裂和非連續(xù)大位移狀態(tài)并模擬運(yùn)動(dòng)全過程。有限元法(FEM)和非連續(xù)變形分析法(DDA)都是數(shù)值流形法的特例。該方法提出后，已廣泛應(yīng)用于許多巖土工程實(shí)際領(lǐng)域，如裂隙擴(kuò)展[3-6]、洞室開挖[7-9]、滲流[10-12]等。

數(shù)值流形法最初開發(fā)建立在線彈性假設(shè)基礎(chǔ)上。聶治豹等[13]將數(shù)值流形法與邊界元法進(jìn)行耦合，提出了一種具有邊界積分形式的數(shù)值流形法來進(jìn)行二維彈性靜力問題研究。徐棟棟等[14]在物理覆蓋位移函數(shù)中增加一種新的局部位移函數(shù)，在不增加自由度的基礎(chǔ)上克服了線性問題，并通過均質(zhì)邊坡彈性力學(xué)算例進(jìn)行了有效性驗(yàn)證。駱少明等[15]采用規(guī)則的六面體數(shù)學(xué)網(wǎng)格建立三維數(shù)值流形有限覆蓋系統(tǒng)，對(duì)三維數(shù)值流形法進(jìn)行了理論推導(dǎo)等研究，建立了線彈性變形三維數(shù)值流形分析的理論體系。上述數(shù)值流形方法多基于線性本構(gòu)關(guān)系，而實(shí)際巖土體結(jié)構(gòu)具有材料和幾何雙重非線性的特性需要進(jìn)行彈塑性分析。目前對(duì)于非線性問題(如彈塑性問題等)的研究還比較少[16，17]，因此，發(fā)展非線性數(shù)值流形法是必要的。由于牛頓-拉普森法多用于非線性求解中，同時(shí)也可以和其它非線性求解方法聯(lián)合使用。王軍祥等[18]將Newton-Raphson法與Arc-Length法相結(jié)合，用來求解針對(duì)巖石彈塑性階段發(fā)生的應(yīng)變軟化等問題。賈碩等[19]將Woodbury公式引入局部非線性分析中，提高了矩陣求逆的效率，在此基礎(chǔ)上利用N-R法和其它方法求解具體算例鋼框架和薄鋼板，并分析了相應(yīng)求解方法對(duì)應(yīng)的計(jì)算誤差、時(shí)間復(fù)雜度以及收斂快慢程度等方面。但是將Newton-Raphson迭代法應(yīng)用于數(shù)值流形非線性問題求解的方面尚未發(fā)現(xiàn)相關(guān)研究工作。

針對(duì)上述問題，本文嘗試將牛頓-拉普森迭代法和數(shù)值流形法相結(jié)合來分析材料非線性問題，提出了基于修正牛頓-拉普森迭代的數(shù)值流形法。在迭代求解過程中只需在初始迭代步建立剛度矩陣，此后迭代步沿用，可降低求解大型工程問題過程的耗時(shí)，提升了非線性求解的效率。

2 NMM覆蓋位移函數(shù)

整體位移函數(shù)由定義在各個(gè)物理覆蓋上的覆蓋位移函數(shù)通過權(quán)函數(shù)連接形成。定義在物理覆蓋Ui上的覆蓋函數(shù)ui(x，y)和vi(x，y)一般級(jí)數(shù)形式為

(1)

其中，fij(x，y)為多項(xiàng)式的基本級(jí)數(shù)；di，2j-1，di，2j為物理覆蓋的自由度，是級(jí)數(shù)每項(xiàng)系數(shù)對(duì)應(yīng)的未知量。覆蓋函數(shù)階數(shù)可取0階到高階，對(duì)應(yīng)于0，1，2階，m分別為1，3，6。

覆蓋函數(shù)為0階時(shí)，其常量函數(shù)形式為

(2)

物理覆蓋的每一個(gè)交集定義為流形單元，流形單元是物理覆蓋的再剖分。若流形單元e是q個(gè)物理覆蓋Ue(1)，Ue(2)，…，Ue(q)的交集，每個(gè)物理覆蓋上的權(quán)函數(shù)為we(r)(x，y)，則單元e的整體位移函數(shù)為

(3)

將式(1)代入式(3)右側(cè)，可得到以變形矩陣Tij(x，y)表達(dá)的整體位移函數(shù)為

(4)

(5)

(6)

其中，{De(r)}為物理覆蓋的自由度，為2m階向量。

3 修正牛頓-拉普森迭代方法

考慮非線性時(shí)，材料的應(yīng)力-應(yīng)變關(guān)系不再滿足胡克定律，其剛度矩陣不是一個(gè)常數(shù)，而與應(yīng)變和位移值有關(guān)。此時(shí)，結(jié)構(gòu)剛度方程是一個(gè)非線性方程組

[K]{u}={Fext}

(7)

其中，[K]為剛度矩陣；{u}為節(jié)點(diǎn)位移向量；{Fext}為節(jié)點(diǎn)荷載向量。

求解非線性方程組的方法包括增量法和迭代法。其中，增量法又分為始點(diǎn)增量法和中點(diǎn)增量法。

增量法的原理是把計(jì)算過程分為幾個(gè)荷載增量步，增量步可以等分，也可以不等分。假定每個(gè)荷載增量步內(nèi)的方程是線性的，在該增量步內(nèi)的剛度矩陣[K]為常數(shù)，不同增量步的剛度矩陣[K]可以不同。每加載一個(gè)荷載增量{ΔFext}，得到該增量所對(duì)應(yīng)的位移增量{Δui}，將所有荷載增量步的位移增量進(jìn)行疊加，得到總位移{u}。

而利用迭代法求解方程組時(shí)，不需劃分荷載步，一次施加全部荷載，然后逐步調(diào)整位移值，使得式(7)的變形式(8)得到滿足。

[K]{u}-{Fext}=0

(8)

牛頓-拉普森迭代法是一種在實(shí)數(shù)域和復(fù)數(shù)域上求解方程的方法，它迫使解在每個(gè)荷載增量末端達(dá)到一定容許范圍內(nèi)的平衡收斂，可有效解決荷載增量法中因誤差積累導(dǎo)致解發(fā)散的問題。由于具有收斂性好，形成方程數(shù)量少，計(jì)算速度快的特點(diǎn)，牛頓-拉普森迭代法被廣泛應(yīng)用于非線性方程組的求解中。

流形單元?jiǎng)偠染仃嚍?/p>

[Ke(r)e(s)]=?[Be(r)]T[Depc][Be(s)]

(9)

其中，[Be(r)]，[Be(s)]均為流形元應(yīng)變矩陣；[Depc]為一致切線模量矩陣；采用四邊形數(shù)學(xué)網(wǎng)格時(shí)，r，s=1，2，3，4。則

(10)

修正牛頓法在傳統(tǒng)牛頓法基礎(chǔ)上做了改進(jìn)。如圖1所示，若采用傳統(tǒng)的牛頓-拉普森迭代法，剛度矩陣應(yīng)在每個(gè)迭代步中根據(jù)非一致切線模量矩陣[Dep]建立，傳統(tǒng)牛頓法的缺點(diǎn)在于剛度矩陣需要在每個(gè)迭代步中重新形成并求逆，在求解大型問題時(shí)尤其費(fèi)時(shí)，且剛度矩陣求逆時(shí)可能會(huì)出現(xiàn)奇異或病態(tài)的問題。修正牛頓-拉普森迭代方法則可以克服以上困難，只需要在初始迭代步時(shí)建立剛度矩陣并求逆，此后的迭代步即可沿用初始迭代步的剛度矩陣，如圖2所示。相比于牛頓法，修正牛頓法雖然降低了迭代的收斂速度，卻大大節(jié)省了計(jì)算時(shí)間。本文采用修正牛頓-拉普森迭代法進(jìn)行分析計(jì)算。

修正的牛頓-拉普森迭代可以表達(dá)為下列方程

(11)

{ui+1}={ui}+{Δui}

(12)

使用修正牛頓-拉普森迭代法分析非線性問題的流程為：

1)確定初始位移{ui}。修正牛頓-拉普森迭代是一種梯度迭代算法，初始值的確定對(duì)算法的收斂速度有重要影響。在初始迭代步，{u0}={0}；

(13)

{Δσi}={σi}-{σi-1}

(14)

其中，{Δσi}為當(dāng)前迭代步單元應(yīng)力{σi}與上一迭代步單元應(yīng)力{σi-1}之間的差值。當(dāng)前迭代步和上一迭代步的應(yīng)力計(jì)算均與選取的本構(gòu)模型有關(guān)，具體過程在第3節(jié)論述。

5)利用式(12)計(jì)算下一迭代步更新后的位移向量{ui+1}，同時(shí)作為下一迭代步的初始位移值；

6)重復(fù)過程2)～5)，直至不平衡力小于一定范圍時(shí)，認(rèn)為結(jié)果收斂，迭代結(jié)束。

4 摩爾-庫侖模型在NMM中的應(yīng)用

土是一種很復(fù)雜的材料，土的本構(gòu)關(guān)系受到諸多因素的影響，如密度、含水量、應(yīng)力歷史等。已提出的土的本構(gòu)模型包括線彈性模型、非線性模型、彈塑性模型等[20，21]。

摩爾-庫侖模型屬于彈塑性模型，主要用于描述巖土介質(zhì)的抗剪破壞行為，適用于模擬一般巖土體的力學(xué)行為，如邊坡穩(wěn)定和地下開挖等。

本構(gòu)模型不同，迭代步內(nèi)的應(yīng)力算法則相應(yīng)不同。以摩爾-庫侖模型為例，其應(yīng)力算法為隱式返回映射算法[22]。具體計(jì)算流程如下：

1)預(yù)測單元彈性應(yīng)力并求主應(yīng)力

根據(jù)已知的單元B矩陣和位移值求得單元應(yīng)變，進(jìn)一步按彈性方法求出單元應(yīng)力作為預(yù)測應(yīng)力值，在主應(yīng)力方向分解后得到主應(yīng)力。

2)檢查應(yīng)力是否進(jìn)入塑性狀態(tài)

若滿足式(15)，說明流形元尚處于彈性狀態(tài)，則可將第1)步所求預(yù)測單元彈性應(yīng)力作為實(shí)際應(yīng)力值。若不滿足式(15)，則進(jìn)入第3)步。

(15)

(16)

3)返回映射。

存在以下三種返回映射類型，按順序依次進(jìn)行：①返回至主平面；②返回至棱邊(右或左)；③返回至頂點(diǎn)。在執(zhí)行當(dāng)前返回映射類型后根據(jù)式(17) 檢查應(yīng)力結(jié)果的有效性，若有效，則直接進(jìn)入第4)步，若檢查為無效，則執(zhí)行下一類型的映射。

σ1≥σ2≥σ3

(17)

4)組裝應(yīng)力張量并更新彈性應(yīng)變。

5 算例

該算例引用自文獻(xiàn)[23]，計(jì)算模型如圖3所示。計(jì)算模型約束條件為：邊坡底邊AB受雙向固定約束；垂直側(cè)邊AF和BC受x向固定約束，y向自由。采用四邊形數(shù)學(xué)網(wǎng)格對(duì)其進(jìn)行網(wǎng)格劃分，則形成如圖4所示的流形元覆蓋系統(tǒng)，共141個(gè)物理覆蓋，114個(gè)流形單元。彈性模量E為100000000 N/m2，泊松比ν為0.35，重度γ為20000 N/m3。摩爾庫侖模型各參數(shù)見表1。

表1 材料參數(shù)

根據(jù)有限元求解理論，剪脹角為5度下安全系數(shù)為1.09。根據(jù)上述提出的求解方法，求解時(shí)采用0階覆蓋函數(shù)，利用修正牛頓-拉普森迭代法對(duì)該邊坡算例進(jìn)行非線性分析，并將程序得到的計(jì)算結(jié)果與有限元理論值進(jìn)行比較，對(duì)比結(jié)果如表2所示。

表2 umax計(jì)算值與理論值

從上述結(jié)果可以看出，計(jì)算結(jié)果與算例有限元理論值比較吻合，說明利用修正牛頓-拉普森迭代法求解材料非線性問題是可行的。

6 結(jié)論

為了推進(jìn)數(shù)值流形法在材料非線性方面的應(yīng)用，據(jù)此將牛頓-拉普森引入到數(shù)值流形法中，提出一種可以求解非線性的修正牛頓-拉普森迭代的數(shù)值流形法。首先采用摩爾-庫侖本構(gòu)模型，求解材料非線性問題，內(nèi)外力差滿足收斂容許值即收斂，以力為基礎(chǔ)的收斂提供了收斂的絕對(duì)量度；其次牛頓-拉普森迭代迫使計(jì)算結(jié)果在每個(gè)計(jì)算步的末端達(dá)到平衡，相比于常見的荷載增量法，可避免荷載增量過程中誤差積累而導(dǎo)致結(jié)果失去平衡的問題；最后由于每個(gè)計(jì)算步中的彈性矩陣假設(shè)為常數(shù)矩陣，則該方法可求解任意階數(shù)的數(shù)值流形問題。通過算例驗(yàn)證了該方法的可行性，為數(shù)值流形法在非線性領(lǐng)域的應(yīng)用奠定理論基礎(chǔ)。