郭德偉

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

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

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

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

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

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

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

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

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

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

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

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

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

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

對上式在B點展開寫為

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

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

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

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

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

圖1 回映算法的計算模式

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

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

(1)當?shù)趎個應變增量計算結(jié)束后，某積分點的應力狀態(tài)為σ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)開始迭代計算。當?shù)降?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，則跳出迭代循環(huán)。TOLER是預先設定的計算精度。此時的應力已更新至屈服面上σn+1(C點)(可能有誤差，但誤差在精度容許的范圍之內(nèi))。對塑性應變、彈性應變和求解狀態(tài)變量進行更新:

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

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

3 Druker－Prager模型

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

式中:

4 UMAT二次開發(fā)

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

ABAQUS主程序調(diào)用UMAT子程序進行求解的過程如下:

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

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

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

5 Druker－Prager模型驗證

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

圖3 數(shù)值模擬分析模型

圖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用戶材料子程序開發(fā)及應用［D］．武漢:華中科技大學，2005