Mohr-Coulomb準(zhǔn)則主應(yīng)力回映算法及在地基問題的應(yīng)用

許夢飛，姜諳男

(大連海事大學(xué) 道路與橋梁工程研究所，大連 116026)

1 引言

M-C(Mohr-Coulomb)屈服準(zhǔn)則是應(yīng)用最廣泛和應(yīng)用時間最長的巖土屈服準(zhǔn)則[1]。M-C準(zhǔn)則在π平面上的投影是不規(guī)則的六角形，導(dǎo)致其屈服函數(shù)和塑性勢函數(shù)的導(dǎo)數(shù)存在奇異點，給數(shù)值計算帶來了困難[2]。為了簡化計算，國際流行的大型有限元軟件ANSYS和ABAQUS以及美國MSC公司的MARC和NASTRAN等多通過采用 D -P(Druker-Prager)準(zhǔn)則近似替代M-C準(zhǔn)則[3]，或?qū)-C屈服面進行角點光滑化計算[4，5]。

D -P 系列屈服準(zhǔn)則包括，(1) M-C外角點外接圓準(zhǔn)則(D -P1)，(2) M-C內(nèi)角點外接圓準(zhǔn)則(D -P2),(3) M-C內(nèi)切圓準(zhǔn)則(D -P3)，(4) M-C等面積圓準(zhǔn)則(D -P4)，(5) M-C匹配 D -P 圓準(zhǔn)則(D -P5)，由于其在π平面上的圖形是一系列圓，因此為數(shù)值計算提供了極大便利[6]。

角點光滑法是指將屈服面尖角磨圓，當(dāng)羅德角超過一定范圍時，將屈服函數(shù)和屈服勢函數(shù)的導(dǎo)數(shù)設(shè)置為0[7]。D -P替代法和角點光滑法都從本質(zhì)上改變了屈服函數(shù)表達式，勢必導(dǎo)致計算結(jié)果出現(xiàn)偏差。

針對此類屈服函數(shù)奇異點問題，文獻[8，9]最早提出了六維應(yīng)力空間中的返回映射算法及本構(gòu)矩陣的建立方法。為了便于問題的求解，減少應(yīng)力空間維度(6維減少到3維)，使回映策略在幾何上更加直觀，文獻[10-13]實現(xiàn)了主應(yīng)力空間中Tresca準(zhǔn)則和H-B(Hoek-Brown)準(zhǔn)則的返回映射法。

本文論述了主應(yīng)力空間中Mohr-Coulomb屈服準(zhǔn)則的流動法則，將奇異點處的流動向量表示為多個屈服面上流動向量的線性組合，解決了數(shù)值求解的難題；給出了應(yīng)力更新過程中回映區(qū)域的判定方法、塑性因子的求解方程和應(yīng)力更新方程，推導(dǎo)了一致切線模量的表達式。在此基礎(chǔ)上,利用C++語言編制了有限元求解程序，通過對經(jīng)典巖土地基問題的求解，驗證了程序的有效性和精確性。

2 主應(yīng)力空間的隱式積分算法

2.1 多屈服面流動法則

Mohr-Coulomb準(zhǔn)則在主應(yīng)力空間中的多屈服面表達式為(拉正壓負)

(1)

(2)

Ni=?Ψi/?σ

(3)

假設(shè)主應(yīng)力位于σ1≥σ2≥σ3區(qū)域，在不失一般性的情況下，只對π平面上1/6區(qū)域進行討論，有四種可能的塑性流動情況，如圖2所示。

(1) 當(dāng)應(yīng)力點在屈服函數(shù)面上時，只有一個塑性因子不為0。

(4)

(1-sinΨ)e3?e3

(5)

(2) 當(dāng)應(yīng)力點位于右棱線上，有兩個塑性因子不為0。

(6)

Nb=N6=(1+sinΨ)e1?e1-(1-sinΨ)e2?e2

(7)

式中N6為與塑性勢函數(shù)Ψ6正交的流動向量。

(3) 當(dāng)應(yīng)力點位于左棱線上，塑性應(yīng)變表達式

圖1 分段線性軟化(硬化)函數(shù)曲線

Fig.1 Piecewise linear hardening(softening) function curves

與式(6)相同，此時

Nb=N2=(1+sinΨ)e2?e2-(1-sinΨ)e3?e3

(8)

式中N2為與Ψ2正交的流動向量。

(4) 當(dāng)應(yīng)力點位于MC準(zhǔn)則的尖點處，6個塑性因子均不為0。

(9)

(10)

當(dāng)應(yīng)力點在屈服平面、棱線和尖點處時，其表達式分別為

(11)

2.2 M-C準(zhǔn)則主應(yīng)力空間應(yīng)力回映算法

對于彈塑性材料，其本構(gòu)關(guān)系取決于加載歷史，具有非線性和非單值對應(yīng)等特點。因此，選取合適的積分算法，在應(yīng)變路徑上對本構(gòu)方程進行積分對于彈塑性材料的有限元求解過程顯得尤為重要。

在求解率相關(guān)的積分算法中，通常采用歐拉算法將本構(gòu)方程離散為多個時間步，由tn時間步中所得的應(yīng)力σn、應(yīng)變εn和內(nèi)變量αn對tn + 1時刻的應(yīng)力狀態(tài)進行唯一求解。然而，求解過程中應(yīng)力和內(nèi)變量的更新值容易出現(xiàn)不滿足屈服函數(shù)的情況，導(dǎo)致計算結(jié)果不精確。

應(yīng)力回映算法具有彈性預(yù)測和塑性修正兩個過程，在每個迭代步中對更新變量進行塑性修正，使其滿足屈服函數(shù)，有效解決了上述問題。

圖2 多屈服面流動法則

Fig.2 Flow rule of multisurface plastic potential function

M-C準(zhǔn)則的主應(yīng)力空間應(yīng)力回映算法步驟如下。

步驟1彈性預(yù)測

已知tn時刻的狀態(tài)變量以及當(dāng)前迭代步的應(yīng)變增量Δε，彈性預(yù)測狀態(tài)如下，

(12)

(13)

則當(dāng)前計算步仍在彈性區(qū)，對所有變量進行更新：

(14)

如果式(13)不成立，則需進行塑性修正，將試算應(yīng)力映射回屈服面上。

步驟2塑性修正

(1) 假設(shè)更新應(yīng)力位于屈服函數(shù)主平面上，則由塑性流動向量式(3)可知，需對塑性因子Δγ進行求解,假設(shè)初值：

(15)

此時相應(yīng)的屈服函數(shù)為

(16)

建立New-Raphson迭代式對Δγ進行求解：

σ1∶=σ2∶=σ3∶=pn + 1

(17)

進行收斂性判斷：

(18)

(19)

圖3 應(yīng)力回映算法過程

Fig.3 Stress mapping progress

當(dāng)|f|≤tol(tol為容許誤差，本文取值為 1e -5)時，計算完成，更新應(yīng)力狀態(tài)：

(20)

(2) 對更新后的應(yīng)力進行判定，如果σ1≥σ2≥σ3，則應(yīng)力更新成立，對應(yīng)力、應(yīng)變和內(nèi)變量進行更新，退出當(dāng)前迭代步；否則假設(shè)返回應(yīng)力位于棱線上。

在M-C準(zhǔn)則的任意棱線上有偏應(yīng)力si=sj，其偏應(yīng)力更新方程為

(21)

T1=1-sinΨ，T2=-2，T3=1+sinΨ

(22)

(1-sinΨ)σ1-2σ2+(1+sinΨ)σ3

(23)

則由圖4的幾何關(guān)系可知，當(dāng)S>0時，更新應(yīng)力應(yīng)在右棱線上；S<0時，更新應(yīng)力在左棱線。

當(dāng)返回應(yīng)力位于右棱線時，根據(jù)塑性流動法則式(7，8)，對Δγa和Δγb進行求解。

設(shè)定迭代計算初始值：

(24)

在右棱線處應(yīng)滿足屈服方程f1=f6=0，即

(25)

(26)

(27)

建立Δγa和Δγb的Newton-Raphson迭代式子：

(28)

(29)

圖4 棱線回映區(qū)域判斷

Fig.4 Selection of edge return

式中d為殘差矩陣，表達式為

(30)

式中

(31)

(32)

進行收斂性判斷：

(33)

(34)

當(dāng)滿足|f1|+|f6|

(35)

當(dāng)返回左棱線時，將式(27,32)進行重寫，其他表達式不變：

(36)

(37)

此時的應(yīng)力更新表達式為

(38)

(3) 對更新后的應(yīng)力重新進行判斷，如果滿足σ1≥σ2≥σ3，返回右(左)棱線成立，否則返回尖點。

(39)

從圖5可以看出，此時有pn + 1=p，即

(40)

設(shè)置迭代初始值：

(41)

由式(41)建立殘差方程：

(42)

(43)

計算殘差，進行收斂性檢驗：

(44)

(45)

(46)

如果|r|≤tol，則計算收斂，對應(yīng)力進行更新：

σ1∶=σ2∶=σ3∶=pn + 1

(47)

2.3 一致切線模量表達式

一致切線模量能夠在求解過程中保持牛頓迭代法的平方收斂速度，提高計算效率，避免材料在塑性狀態(tài)時產(chǎn)生偽加載和卸載問題[14,15]。本質(zhì)上看，彈塑性一致切線模量需要求得應(yīng)力張量對應(yīng)變張量的微分表達式，即

(i,j=1,2,3) (48)

(49)

當(dāng)返回應(yīng)力位于主屈服面上時，對式(20)求導(dǎo)可得

(50)

圖5 回映至尖點

Fig.5 Return to apex

根據(jù)一致性條件，對屈服函數(shù)進行求導(dǎo)可得d Δγ。

(51)

式中a的表達式同式(31)。

當(dāng)返回應(yīng)力位于右(左)棱線時，求解過程同上，對屈服函數(shù)f1和f6(或f1和f2)和式(36)(或式(39))進行求導(dǎo)，與線性關(guān)系式進行聯(lián)立即可求得一致切線模量表達式。

當(dāng)返回應(yīng)力位于尖點時，對式(39)進行求導(dǎo)得

(52)

由一致性條件對式(40)進行求導(dǎo)可得

(53)

由式(53)可得

(54)

將式(54)代入式(52)求得一致性切線模量表達式為

(55)

3 主應(yīng)力回映算法有限元程序編制

本程序采用面向?qū)ο蟮腃++語言編制而成，并在VS2010環(huán)境下調(diào)試通過。

非線性巖土材料的有限元求解通常采用增量迭代法，將總荷載f分為多個荷載增量Δf逐級施加進行求解，非線性有限元增量方程為

KΔdi=Δfi

(56)

(57)

式中K為剛度矩陣,Δdi為位移增量,f0為荷載初始值,d0為位移初始值。

采用Newton-Raphson迭代法對式(56)進行求解，具體步驟如下。

(1) 設(shè)置初始位移d0、應(yīng)變ε0和應(yīng)力σ0。

(58)

(59)

(60)

(6) 利用式(61)求解位移增量Δdi

(61)

4 巖土地基問題的求解

圖6 有限元計算流程

Fig.6 FEM program flow char

當(dāng)荷載p=40 kN/m2時，地基內(nèi)的位移和應(yīng)力云圖如圖9所示。

圖7 地基承載力計算模型

Fig.7 Calculation model of foundation bearing capacity

圖8 不同荷載作用下塑性區(qū)演化過程

Fig.8 Evolution progress of plastic zone

傳統(tǒng)有限元計算軟件在處理Mohr-Coulomb尖點問題時，通常采用基于廣義米塞斯屈服準(zhǔn)則的DP系列準(zhǔn)則或角點光滑法等手段進行近似計算。圖10為DP系列準(zhǔn)則在π平面上的圖形。

圖9 應(yīng)力和位移云圖

Fig.9 Contour of stress and displacement

圖10 DP系列屈服準(zhǔn)則在π平面上的圖形

Fig.10 Yield surface on the deviatoric plane of DP criterion

為了驗證本文所編程序的正確性，分別采用本文方法、等效DP準(zhǔn)則法和角點光滑法對巖土地基問題進行求解。對于一個承受均勻垂直荷載的半無限大、無重量地基，其承載力的Prandtl解為

(62)

qu=(π+2)c

(63)

由文獻[17]可知，對于平面應(yīng)變問題，在關(guān)聯(lián)流動法則下采用Mohr-Coulomb內(nèi)切圓(D -P3)或在非關(guān)聯(lián)流動法則下采用Mohr-Coulomb匹配應(yīng)力圓準(zhǔn)則(D -P5)所得結(jié)果最為精確。因此,本文利用文獻[18]的有限元程序，只選用D -P系列中的D -P3和D -P5準(zhǔn)則與角點光滑法進行對比計算，結(jié)果列入表1。

由表1可知，屈服準(zhǔn)則和內(nèi)摩擦角φ對計算結(jié)果精度均有較大影響；關(guān)聯(lián)流動法則下，本文算法最大計算誤差為0.65%，采用角點光滑法計算所得最大誤差為2.01%，D -P3準(zhǔn)則下計算最大誤差為2.30%；非關(guān)聯(lián)流動法則下，本文算法最大計算誤差為0.51%，角點光滑法計算所得最大誤差為2.64%，D -P5準(zhǔn)則下計算最大誤差為1.84%。從計算結(jié)果對比可以看出，M-C準(zhǔn)則的主應(yīng)力回映算法不需要對屈服方程進行簡化處理，因此其計算結(jié)果具有較高的精度。

表1 不同計算方法下地基承載力計算結(jié)果(單位：kN/m2)

Tab.1 Foundation bearing capacity using different method (unit:kN/m2)

計算方法φ/(°)關(guān)聯(lián)流動法則非關(guān)聯(lián)流動法則(Ψ=φ/2)10151015Prandtl(M1)41.58角點光滑法(M2)42.42DP3/DP5(M3)42.49本文(M4)41.75(M2-M1)/M10.0201(M3-M1)/M10.0219(M4-M1)/M10.004054.71-55.6042.5955.9742.2455.0741.740.01630.02430.02300.01590.00650.0038-56.1555.7254.990.02640.01840.0051

圖11 不同硬化規(guī)律下地基P -S曲線

Fig.11 P -S curve of foundation under different hardening rule

5 結(jié) 論

本文提出了Mohr-Coulomb屈服準(zhǔn)則在主應(yīng)力空間中的應(yīng)力回映算法，解決了其在棱線和尖點處的數(shù)值奇異點問題；利用C++語言編制了有限元計算程序，對巖土地基承載力進行求解，計算結(jié)果同兩種近似方法相比更接近理論解，驗證了程序的精確性和有效性；通過硬化參數(shù)的設(shè)置，所編程序還能夠反映硬化(軟化)巖土材料總體承載力的變化規(guī)律，為工程設(shè)計提供了有利工具。