在ABAQUS中開 發(fā)實現(xiàn)Mortara界面循環(huán)本構(gòu)模型

任 宇，陳仁朋，陳云敏

(1.中國電力工程顧問集團(tuán)華東電力設(shè)計院，上海 200063; 2. 浙江大學(xué)建筑工程學(xué)院，浙江 杭州 310058)

ABAQUS由美國的HKS公司研制開發(fā)，是目前世界公認(rèn)的最為先進(jìn)的大型通用非線性有限元分析軟件，在幾何、材料非線性和復(fù)雜接觸問題等方面的分析能力居世界領(lǐng)先水平，特別適用于巖土工程這種包含高度材料非線性及復(fù)雜界面接觸問題的求解．

在巖土工程數(shù)值分析中，結(jié)構(gòu)-土界面特性的模擬一直是其中的重點及難點．特別對于樁、錨桿和土釘?shù)冗@些主要通過界面將荷載傳遞至周圍土中的土工結(jié)構(gòu)來說，界面的力學(xué)特性往往決定了整個土工結(jié)構(gòu)的承載及變形性狀[1]．

ABAQUS中自帶的界面模型包括切向摩擦模型和法向接觸模型兩類，其中摩擦模型用以模擬接觸面間的剪應(yīng)力傳遞規(guī)律，一般為庫倫摩擦模型；法向接觸模型則用來模擬接觸面間法向應(yīng)力傳遞規(guī)律，包括“硬接觸”（hard contact）和“軟接觸”(softened contact)兩種模型[2]．目前ABAQUS中的界面模型在使用時是相互獨立的，較少考慮了切向和法向間的耦合作用，如現(xiàn)有界面模型中就沒有考慮剪切作用對于法向應(yīng)力-應(yīng)變關(guān)系的影響等．然而這與界面剪切試驗中所觀察到的現(xiàn)象并不吻合：剪切通常會引起界面土體體積的改變，進(jìn)而會影響到法向應(yīng)力的大小[1-3]．更為重要的是，現(xiàn)有模型僅是針對于靜載下的界面性狀的，軟件中還沒有能夠用于模擬界面循環(huán)加載性狀的本構(gòu)模型．這一缺憾無疑極大限制了 ABAQUS在土工結(jié)構(gòu)循環(huán)加載性狀數(shù)值分析中的應(yīng)用范圍．

為了彌補自帶界面模型的不足，ABAQUS基于其子程序擴(kuò)展平臺，提供了允許用戶自定義所需界面特性的用戶界面子程序UINTER(User subroutine to define surface interaction behavior for contact surfaces)．UINTER 能夠讓用戶定義并使用ABAQUS中所沒有的界面本構(gòu)模型．這一功能極大地拓展了軟件在土與結(jié)構(gòu)相互作用分析中的應(yīng)用范圍．

基于ABAQUS提供的UINTER子程序，完成了 Mortara界面循環(huán)模型(Mortara cyclic interface model)的開發(fā)工作．對界面循環(huán)剪切試驗進(jìn)行了建模分析，將 ABAQUS有限元分析結(jié)果與文獻(xiàn)給出的數(shù)值解結(jié)果進(jìn)行了對比．結(jié)果表明，UINTER的求解效率和精度令人滿意，有望成為 ABAQUS中模擬結(jié)構(gòu)-土界面循環(huán)加載性狀的強有力分析工具．

1 Mortara界面循環(huán)模型

1.1 Mortara界面循環(huán)模型

Ghionna 和 Mortara[4]基于經(jīng)典彈塑性理論建立了一個模擬界面靜載性狀的二維彈塑性本構(gòu)模型．模型的屈服面方程借鑒了Vermeer[5]的偏應(yīng)力屈服面的概念，在 τ－σn空間內(nèi)體現(xiàn)為是一條原點出發(fā)微彎的曲線．模型采用非關(guān)聯(lián)流動法則，模型勢函數(shù)則與傳統(tǒng)劍橋模型的相類似．

與傳統(tǒng)彈塑性模型類似，在上述界面靜力模型中，卸載及再加載過程中的應(yīng)力路徑一直處于屈服面內(nèi)，此時界面的應(yīng)力-應(yīng)變關(guān)系均表現(xiàn)為彈性，并不能反映出循環(huán)加載下界面土顆粒重排導(dǎo)致塑性變形累積的特征．為此，Mortara等[6]對上述界面靜力本構(gòu)模型進(jìn)行了發(fā)展，在原有屈服面內(nèi)引入了一個小的運動屈服面；當(dāng)應(yīng)力點落在運動屈服面內(nèi)時，界面僅產(chǎn)生彈性變形，如若應(yīng)力試圖穿過運動屈服面，就會引發(fā)彈塑性變形，同時會拖移運動屈服面在外邊界面內(nèi)發(fā)生移動．該模型可以較好地表征出界面在卸載及再加載下的塑性變形產(chǎn)生過程．

1.2 本構(gòu)模型的數(shù)學(xué)表達(dá)式

1.2.1 屈服函數(shù)

Mortara界面循環(huán)模型包含有外邊界面和內(nèi)部運動屈服面（簡稱內(nèi)面）兩個屈服面．模型的外邊界面和內(nèi)部屈服面的方程分別為

式中：α為硬化-軟化參數(shù)，它的取值決定了外邊界面的大小及位置；α0為移動屈服面的邊線與中心軸線間的夾角；為以局部參考坐標(biāo)系中的應(yīng)力值，該坐標(biāo)是將原坐標(biāo)的τ=0軸旋轉(zhuǎn)θ得到的．

圖1 界面循環(huán)模型的屈服面Fig. 1 Yield surfaces of the interface model

1.2.2 映射法則

當(dāng)應(yīng)力點在邊界面內(nèi)移動時，需要建立映射法則來確定塑性模量隨屈服面移動的變化情況． 外邊界面內(nèi)的塑性模量由下式計算：

上式與 Dafalias和 Popov[7]所提出的映射函數(shù)形式相類似．其中:rθ為內(nèi)外屈服面對應(yīng)邊界線間的角度；0rθ則為應(yīng)力路徑發(fā)生折返穿過彈性區(qū)再次達(dá)到內(nèi)屈服面上時的rθ．Hα為加載至外邊界面時所對應(yīng)的塑性模量．h則由下式計算

式中: n和hχ為模型參數(shù)；0s為一任意參考位移，用于將hχ無量綱化；max[]p u 為最大法向累積塑性變形；inh 則由下式計算

其中，γ和δ為模型參數(shù)．式(4)中Rαθ則由下式定義

1.2.3 循環(huán)流動法則

界面模型在第i周循環(huán)下流動法通過下式表示

式中：η為應(yīng)力比，η =τ/ σn；d為剪脹率．可以看出，模型采用的流動法則在η?d平面內(nèi)表現(xiàn)為線性，為第i周循環(huán)中流動法則在η軸上的截距．模型中通過與法向塑性累積變形的關(guān)系來考慮循環(huán)對于流動法則的影響

式中：bχ為模型參數(shù)．

Mortara界面循環(huán)模型是一種高度非線性的界面彈塑性模型，能夠較為精確地對循環(huán)加載下結(jié)構(gòu)-土界面的復(fù)雜力學(xué)行為進(jìn)行描述，特別是能夠很好地對循環(huán)下界面土體的剪縮性狀進(jìn)行模擬．該模型是一個非常有應(yīng)用價值和發(fā)展前途的巖土本構(gòu)模型．

2 界面子程序開發(fā)

2.1 UINTER子程序

為了方便及鼓勵用戶自行定義符合特定問題的模型及單元，ABAQUS提供了大量的用戶自定義子程序（user subroutine），UINTER是其中支持用戶自定義界面特性的子程序．UINTER功能非常強大，不但可用于定義接觸面間的力學(xué)行為，也可對界面間的熱力學(xué)性狀進(jìn)行定義．用戶可以將所需的界面本構(gòu)關(guān)系編寫入UINTER，并在有限元分析中進(jìn)行調(diào)用．

ABAQUS默認(rèn)采用Newton-Raphson迭代法求解非線性方程組，當(dāng)每一個增量步開始時，ABAQUS主程序都會在界面每個積分點上調(diào)用UINTER，傳入當(dāng)前狀態(tài)應(yīng)力、應(yīng)變及狀態(tài)變量等相關(guān)參數(shù)；UINTER則根據(jù)傳入?yún)?shù)，求解出主程序給出應(yīng)變增量下所相應(yīng)的應(yīng)力增量，并更新應(yīng)力及相關(guān)其他參數(shù)，并向主程序提高更新后的雅克比矩陣；主程序根據(jù)應(yīng)力增量計算出殘余荷載，并根據(jù)殘余荷載求得相應(yīng)誤差，如果不滿足指定誤差條件，ABAQUS將進(jìn)行迭代，直至達(dá)到收斂準(zhǔn)則，然后進(jìn)入下一增量步的計算[2]．

2.2 帶誤差控制的顯式積分算法

Newton-Raphson法中的關(guān)鍵步驟是確定每次迭代中的殘余荷載，為此必須計算本次迭代結(jié)束時的總應(yīng)力．通過對本構(gòu)方程積分得到應(yīng)力增量的方法叫做應(yīng)力積分法，包括顯式積分算法和隱式積分算法兩種．

傳統(tǒng)顯式歐拉積分法的計算精度較低，預(yù)測應(yīng)力往往會偏離實際屈服面，易于產(chǎn)生屈服面漂移．Potts等[8]提出了旨在消除屈服面偏移的修正算法，極大地提高了顯示積分算法的精度．Potts等[9]將帶誤差控制的顯式歐拉積分算法與隱式積分法進(jìn)行了比較，指出兩種算法均能得到令人滿意的積分結(jié)果，但是相較于隱式算法，顯式算法對于復(fù)雜本構(gòu)模型的適用性更好．對于一些復(fù)雜的本構(gòu)模型，往往難以推導(dǎo)出適合于隱式回歸計算的公式．因此，本文選取 Potts等提出的帶誤差控制的修正歐拉顯式積分法對本構(gòu)方程進(jìn)行積分計算．

顯式積分計算步驟主要包括兩步：①彈性預(yù)測；②彈塑性應(yīng)變積分．彈性預(yù)測過程中假設(shè)變形增量為純彈性，按照彈性本構(gòu)關(guān)系計算相應(yīng)的應(yīng)力增量，如果增量結(jié)束時的應(yīng)力仍在屈服面內(nèi)，則說明彈性預(yù)測得到的應(yīng)力增量是準(zhǔn)確的；如果超過屈服面，則說明應(yīng)變不全是彈性的，應(yīng)該對彈塑性部分進(jìn)行積分以獲得正確的應(yīng)力增量．由于循環(huán)模型中存在內(nèi)外兩個屈服面，包含了兩種勢函數(shù)及硬化法則，因此還需要判斷顯式積分過程中的應(yīng)力狀態(tài)是否超過了外邊界面，如果超過則需要采用外邊界面所對應(yīng)的勢函數(shù)及硬化法則計算相應(yīng)的彈塑性方程矩陣．

每個次階計算結(jié)束時，需要對應(yīng)力及硬化參數(shù)計算結(jié)果進(jìn)行檢查，判斷其是否滿足屈服準(zhǔn)則，如果不滿足即可以轉(zhuǎn)入屈服面偏移修正計算流程中，只有滿足了屈服準(zhǔn)則后才轉(zhuǎn)入下一次階計算流程．詳細(xì)的積分步驟在文獻(xiàn)[10]中給出．

2.3 子程序結(jié)構(gòu)

基于本構(gòu)模型公式和應(yīng)力積分算法，參照ABAQUS用戶界面子程序的接口規(guī)范，進(jìn)行了UINTER的程序設(shè)計及代碼編寫．代碼由 Fortran語言編寫．圖2中給UINTER的計算流程．文獻(xiàn)[10]中給出了具體的程序代碼．

圖2.UINTER計算流程圖Fig 2. Flow chart for UINTER calculation

3 數(shù)值驗證

為了對所開發(fā)的UINTER子程序進(jìn)行驗證，筆者將自定義 UINTER程序接入 ABAQUS/Standard中，對界面循環(huán)剪切試驗進(jìn)行了數(shù)值模擬，并將計算結(jié)果與文獻(xiàn)[6]提供數(shù)值計算結(jié)果進(jìn)行了對比．用于驗證的有限元模型如圖3所示．模型由兩部分組成，一個尺寸為1.0 m×1.0 m×1.0 m的三維彈性正方體，壓在一尺寸為3.0 m×3.0 m×1.0 m的剛體上．彈性體與剛體間設(shè)置接觸面對，接觸面間力學(xué)特性由自定義UINTER子程序定義．彈性體劃分為1個單元，剛體劃分為9個單元．彈性體的彈性模量設(shè)為10GPa．為保證剪切過程中彈性體水平方向的尺寸保持恒定，彈性體的泊松比ν=0．

圖3.界面循環(huán)剪切試驗有限元分析模型Fig. 3 FE model for cyclic interface shear test

Mortara界面循環(huán)模型的計算參數(shù)取自文獻(xiàn)[6]，見表1．在所對應(yīng)界面循環(huán)剪切試驗中，試驗土體為Toyouta砂，界面為鋁板，粗糙度為30 μm，相對密實度DR=85%，法向初始應(yīng)力0nδ=150 kPa．分析中選取常法向剛度條件下的剪切試驗進(jìn)行分析，法向剛度K=500 kPa/mm，w△=0.8mm，[wmean]=0.4 mm．計算開始時首先在彈性體的法向施加初始法向應(yīng)力0nδ，剪切過程采用位移控制方式施加，對彈性體施加水平x向的位移，以模擬土與結(jié)構(gòu)間的剪切過程．

表1.算例的Mortara模型參數(shù)Tab. 1 Parameters of the model

數(shù)值計算與解析解計算結(jié)果對比情況如圖 4-6所示．圖 4-6中分別給出了循環(huán)剪切下的剪應(yīng)力 τ與剪切位移w、法向應(yīng)力σn與剪切位移w、剪應(yīng)力τ與法向應(yīng)力σn間的變化規(guī)律．從圖中可以看出，本文數(shù)值模擬結(jié)果與文獻(xiàn)解析解間有著很好的一致性，驗證了本子程序?qū)缑嫜h(huán)加載性狀的模擬能力，說明本文UINTER成功實現(xiàn)了Mortara界面循環(huán)模型的開發(fā)．

圖4.循環(huán)剪切試驗中～wτ關(guān)系的解析解與數(shù)值解對比Fig. 4 Comparison between relationships of ～wτ obtained from analytical and numerical solution

圖5.循環(huán)剪切試驗中 σ n ～ w關(guān)系的解析解與數(shù)值解對比Fig. 5 Comparison between relationships of σ n ～w obtained from analytical and numerical solution

從圖4及圖5中的模擬結(jié)果可以清晰看出，循環(huán)剪切下的界面將會產(chǎn)生明顯的滯回特性，同時循環(huán)剪切下的界面土體會隨循環(huán)產(chǎn)生剪縮，由此導(dǎo)致界面法向減小，剪應(yīng)力弱化．隨著循環(huán)次數(shù)的增加，剪應(yīng)力及法向應(yīng)力會逐漸趨于穩(wěn)定，并在τ～σn呈現(xiàn)出典型的“蝶型”變化規(guī)律（圖6）．這些與界面循環(huán)剪切試驗中所觀測到的現(xiàn)象是一致的．這不但證明了Mortara界面循環(huán)模型對于復(fù)雜循環(huán)特性的精確描述能力，也驗證了本子程序運算的正確性．

圖6.循環(huán)剪切試驗中 τ ～ σ n關(guān)系的解析解與數(shù)值解對比Fig. 6 Comparison between relationships of τ～σn obtained from analytical and numerical solution

4 結(jié) 語

Mortara界面循環(huán)模型是一個優(yōu)秀的界面本構(gòu)模型，能夠很好地對循環(huán)荷載作用下土-結(jié)構(gòu)界面復(fù)雜的力學(xué)特性進(jìn)行描述，包括應(yīng)力-應(yīng)變滯回特性、法向體積剪縮、剪切強度弱化等．通過本文研究，得到以下結(jié)論：

基于ABAQUS用戶子程序開發(fā)平臺UINTER，實現(xiàn)了對Mortara界面循環(huán)模型的二次開發(fā)．算例結(jié)果顯示，自開發(fā)的UINTER子程序具有很好的運算穩(wěn)定性及精度，能夠模擬出土-結(jié)構(gòu)界面的循環(huán)加載性狀．

利用 ABAQUS強大的非線性有限求解能力及優(yōu)秀的前后處理界面，結(jié)合本文所開發(fā)的 Mortara界面循環(huán)模型的UINTER子程序，有望實現(xiàn)對循環(huán)荷載作用下土與結(jié)構(gòu)相互作用問題數(shù)值分析．

子程序的成功開發(fā)極大擴(kuò)展了該界面本構(gòu)模型在巖土研究及設(shè)計中的應(yīng)用，同時也為其他界面本構(gòu)模型在ABAQUS中的開發(fā)提供了借鑒和參考．

References

[1] MORTARA G, MANGIOLA A, GHIONNA V N. Cyclic shear stress degradation and post-cyclic behavior from sand -steel interface direct shear tests[J]. Canadian Geotechnical Engineeing Journal, 2007, 44(7): 739-752.

[2] ABAQUS user’s and theory manuals[M]. Rhode Island,Hibbitt: Karlsson & Sorensen Inc., 2004.

[3] FAKHARIAN K, EVGIN E. Cyclic simple-shear behaviour of sand-steel interfaces under constant normal stiffness condition[J]. Journal of Geotechnical and Geoenvironmental Engineering, ASCE, 1997, 123(12):1096-1105.

[4] GHIONNA V N, MORTARA G. An elastoplastic model for sand-structure interface behavior[J]. Géotechnique,2002, 52(1): 41-50.

[5] VERMEER P A. A double hardening model for sand[J].Géotechnique, 1978, 28(4): 413-433.

[6] MORTARA G, BOULON M, GHIONNA V N. A 2-D constitutive model for cyclic interface behavior[J].International Journal for Numerical Analytical Methods in Geomechanics, 2002, 26(11): 1071-1096.

[7] DAFALIAS Y F, POPOV E P. Plastic internal variables formalism of cyclic plasticity[J]. ASME, Journal of Applied Mechanics, 1976, 43(4): 645-651.

[8] POTTS D M, GENS A. A critical assessment of methods of correcting for drift from yield surface in elasto-plastic finite element analysis[J]. International Journal for Numerical and Analytical Method in Geomechanics,1985, 9: 149-159.

[9] POTTS D M, GANENDRA D. An evaluation of substepping and implicit stress point algorithms[J].Computer Method In Applied Mechanics Engineering,1994, 119(3/4): 341-354.

[10] 任宇. 長期豎向循環(huán)荷載作用下樁的累積變形特性試驗及理論研究[D]. 杭州: 浙江大學(xué), 2013.REN Yu. Model test and theoretical study on the long-term deformation of single pile to cyclic axial loading [D]. Hangzhou: Zhejiang University, 2013.