巖土彈塑性模型的回映算法及ABAQUS子程序開發

郭德偉

(中國土木工程集團有限公司，北京100038)

數值分析是巖土工程中最重要的分析手段之一，為巖土工程的設計和建設提供極大的計算支持，因此，愈來愈受到工程界的重視和應用。巖土工程彈塑性分析的難點之一是巖土材料非線性彈塑性本構模型在計算機上的程序化，這需要對本構方程進行積分，以便得到新的應力增量。在復雜的應力路徑加載條件下，很難給出材料(包括巖土材料)彈塑性應力增量的顯示積分解析式，而只能進行數值積分，在過去的數十年，在這方面有大量的文獻報道［1－8］。這些文獻報道的方法大致分為兩類:顯示積分方法［1－2］和隱式向后歐拉算法［3－8］。隱式向后歐拉算法最早由 Krieg 等人提出［3］，經過發展，目前已成為應用最廣的算法。

巖土類材料大多為遵循Mohr－Coulomb屈服準則，在主應力空間內其π平面(偏平面)內的屈服函數和塑性勢函數的跡線存在不光滑的角點，在這些點上，這兩個函數的導數存在不連續的情況，將會造成數值計算的角點奇異。在處理角點處導數不連續的問題時，Clausen 和 Damkilde［4－6］做了卓有成效的研究工作。其他更一般的方法則是對屈服函數和塑性勢函數跡線進行光滑處理［9，10］。

本文首先介紹了巖土理想彈塑性本構模型的一般形式，然后再從一般形式推導出彈塑性計算的迭代格式和計算方法。最后在ABAQUS有限元軟件平臺上對Druker－Prager理想彈塑性本構模型進行用戶子程序的二次開發。并將數值模擬結果同巖石試樣的單軸壓縮試驗結果作對比。所做的工作只是對UMAT子程序開發的初步學習。以后還將開展深入的學習和應用。

1 巖土彈塑性模型的一般形式

在不考慮材料彈塑性耦合變形的情況下，材料的彈塑性應變增量可以分解為彈性部分與塑性部分之和

在理想塑性條件下，當應力滿足屈服條件，發生塑性變形，屈服函數和一致性條件可以表示為

在式(1)和式(2)中，ε和σ表示工程應變和應力矢量，其分量展開寫為

材料的應力增量由彈性應變確定

D是材料的彈性剛度矩陣，是一個6×6的對稱正定矩陣。對于一個有限應變增量，對上式進行積分，可得

上式說明在彈塑性計算時，一個有限應力增量是彈性應力增量與塑性應力增量之差，圖1形象地說明了這一算法的核心思想:當積分點的應力處于A點σA時，增加一個總應變Δε，它與彈性剛度矩陣的乘積即為彈性應力增量Δσe，A點的應力狀態疊加上這一應力增量后得到彈性應力預測點B的應力值σB，它已超出實際的屈服面之外，而此時，實際的應力狀態處于屈服面上的C點，即σC，故

因此，彈塑性計算的關鍵是計算出塑性應力增量Δσp，然后將彈性預測應力拉回到屈服面上。

根據塑性位勢理論，材料的塑性應變垂直于塑性勢面，且

λ為非負的塑性因子，G=G(σ)為材料的塑性勢函數，當屈服函數F(σ)=G(σ)時，材料是關聯流動的。但對于巖土材料，一般地，F(σ)≠G(σ)，所以是非關聯流動的。將(7)式左乘彈性剛度矩陣，然后積分可得到塑性應力增量

對上式在B點展開寫為

另外式(4)中，應力增量還可以寫成彈塑性剛度矩陣與總應變增量的乘積形式

其中Dep是材料的彈塑性剛度矩陣

當材料為非關聯流動時，塑性剛度矩陣Dp是非對稱矩陣，所以必須使用非對稱線性方程的求解策略。

從式(9)，可以看出，求解塑性應力增量Δσp的關鍵是確定塑性因子增量Δλ的值。將(1)式中屈服方程的左邊在B點做一階泰勒級數展開，這樣屈服方程就可寫為

將式(12)、式(6)和式(9)聯立，即可得到

圖1 回映算法的計算模式

2 巖土彈塑性計算的迭代格式

通常情況下，對于一個確定的應變增量Δ ε，并不能直接從彈性預測應力σB計算出σC，而要經過多次迭代才能修正到真實的應力狀態。每迭代一次，應力、彈性應變和塑性應變就更新一次，見圖2。因此，巖土彈塑性計算的迭代格式可以做如下敘述:

(1)當第n個應變增量計算結束后，某積分點的應力狀態為σn，彈性應變和塑性應變分別為εen和εpn;

(2)此時，開始第n+1個應變增量Δ εn+1的計算。計算得到的彈性預測應力為σ(0)n+1(B點)，對Δλ(0)n+1賦初始值為0。此時εe(0)n+1=εen+Δ εn+1、εp(0)n+1=εpn。將 σ(0)n+1代入屈服方程中，若F(σ(0)n+1)＞0，則要進行第(3)步的迭代計算;

(3)開始迭代計算。當迭代到第 k次時，求解 F(σ(k)n+1)、(aF/aσ)(k)n+1、(aG/aσ)(k)n+1，由式(13)、式(9)和式(6)計算出第k次迭代完成后的Δλ(k)n+1、Δσp(k)n+1和σ(k)n+1。第k次迭

代完成后的塑性應變增量為

(4)每迭代一次，判斷一次屈服情況。假設迭代至第m次時，應力滿足≤TOLER，則跳出迭代循環。TOLER是預先設定的計算精度。此時的應力已更新至屈服面上σn+1(C點)(可能有誤差，但誤差在精度容許的范圍之內)。對塑性應變、彈性應變和求解狀態變量進行更新:

(5)計算彈塑性剛度矩陣:在更新后的應力狀態點σn+1處，計算(aF/aσ)n+1、(aG/aσ)n+1，將其代入式(11)中，即可計算出Dep。

圖2 塑性應力的迭代算法示意

3 Druker－Prager模型

使用靜水壓力p和廣義剪應力q表述的的D－P(Druker－Prager)屈服準則為設塑性勢函數與屈服函數具有相似的形式，寫為式(14)和式(15)中，α、β和d是材料的力學參數，稱為Druker－Prager準則的廣義內摩擦角，廣義剪脹角和廣義單軸抗壓強度。可以通過Mohr－Coulomb準則的黏聚力c、內摩擦角φ和剪脹角φ換算得到。Druker－Prager模型的屈服函數和塑性勢函數的一階導數為:

式中:

4 UMAT二次開發

ABAQUS提供了用FORTRAN語言編寫的子程序接口，供用戶二次開發之用。UMAT子程序用來定義材料的物理特性，在每一個增量步的末尾，應力和與求解相關的狀態變量需要更新。在UMAT中，通過迭代獲得真實應力以及彈性、塑性應變，從而更新硬化參量及體積模量與剪切模量等狀態變量，通過DDSDDE數組提供材料本構模型的Jacobian矩陣，即應力增量對應變增量的變化率。UMAT子程序主要包括以下幾個部分:子程序定義語句、參數說明、用戶定義的。局部變量說明、用戶編寫的主體語句、子程序返回和結束語句。

ABAQUS主程序調用UMAT子程序進行求解的過程如下:

(1)在荷載步n+1開始時，ABAQUS主程序為UMAT子程序提供前n個荷載步計算結束后的應力矢量tnσ、總應變增量tnΔε和時間增量Δt;

(2)UMAT得到ABAQUS主程序提供的數據之后，根據UMAT設定的算法計算出材料的應力應變狀態。計算完成后將應力應變數據和Jacobian矩陣返回給ABAQUS主程序。第n+1個荷載步在UMAT內的應力應變計算完成;

(3)ABAQUS主程序根據UMAT返回的Jacobian矩陣(在此應力應變狀態下的剛度矩陣)按照有限元理論理論形成第n+1個荷載步計算所需的總體剛度矩陣，然后進行力與位移的計算。

5 Druker－Prager模型驗證

開發了Druker－Prager理想塑性的UMAT本構程序，對室內長方體巖樣的單軸壓縮試驗進行有限元模擬。長方體巖樣高度100mm，上下底面邊長50mm，沿模型垂直方向施加位移載荷，所施加的豎向位移為－2mm。材料參數如下:彈性模量10 GPa、泊松比0.25、粘聚力17MPa、內摩擦角44°，計算時設置剪脹角為40°。圖3和圖4分別是數值模擬的分析模型和巖樣單軸情況下的應力應變曲線。

圖3 數值模擬分析模型

圖4 巖樣的軸向應力應變曲線

［1］ S．W．Sloan，A．J．Abbo，D．Sheng．Refined explicit integration of elastoplasticmodels with automatic error control［J］．Engineering Computations，24(2001):121 –154

［2］ T．W．So·owski，D．Gallipoli．Explicit stress integration with error controlforthe Barcelona Basic ModelPartI:Algorithms formulations［J］．Computers and Geotechnics，37(2010):59 – 67

［3］ R．Krieg，D．Krieg．Accuracies of numerical solutionmethods for the elastic－ perfectly plasticmodel［J］．ASME Journal of Pressure Vessel Technology，(1977):510 –515

［4］ J．Clausen，L．Damkilde，L．Andersen．Efficient return algorithms for associated plasticity withmultiple yield planes［J］．International Journal for Numerical Methods in Engineering，66(2006):1036–1059

［5］ J．Clausen，L．Damkilde，L．Andersen．An efficient return algorithm for non－associated plasticity with linear yield criteria in principal stress space［J］．Computers and Structures，85(2007):1795–1807

［6］ J．Clausen，L．Damkilde．An exact implementation of the Hoek –Brown criterion for elasto－plastic finite element calculations［J］．International Journal of Rock Mechanics and Mining Science，45(2008):831–847

［7］ X．Wang，L．B．Wang，L．M．Xu．Formulation of the returnmapping algorithm for elastoplastic soilmodels［J］．Computers and Geotechnics，31(2004):315 –338

［8］ J．Huang，D．V．Griffiths．Observations on Return Mapping Algorithms for Piecewise Linear Yield Criteria［J］．International Journal of Geomechanics．8(2008):253 –265

［9］ A．J．Abbo，S．W．Sloan．A smooth hyperbolic approximation to the Mohr– Coulomb yield criterion．Computers and Structures．54(1995):427–441

［10］ 楊曼娟．ABAQUS用戶材料子程序開發及應用［D］．武漢:華中科技大學，2005