在ABAQUS中開 發實現Mortara界面循環本構模型

任 宇，陳仁朋，陳云敏

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

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

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

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

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

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

1 Mortara界面循環模型

1.1 Mortara界面循環模型

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

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

1.2 本構模型的數學表達式

1.2.1 屈服函數

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

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

圖1 界面循環模型的屈服面Fig. 1 Yield surfaces of the interface model

1.2.2 映射法則

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

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

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

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

1.2.3 循環流動法則

界面模型在第i周循環下流動法通過下式表示

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

式中：bχ為模型參數．

Mortara界面循環模型是一種高度非線性的界面彈塑性模型，能夠較為精確地對循環加載下結構-土界面的復雜力學行為進行描述，特別是能夠很好地對循環下界面土體的剪縮性狀進行模擬．該模型是一個非常有應用價值和發展前途的巖土本構模型．

2 界面子程序開發

2.1 UINTER子程序

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

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

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

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

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

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

每個次階計算結束時，需要對應力及硬化參數計算結果進行檢查，判斷其是否滿足屈服準則，如果不滿足即可以轉入屈服面偏移修正計算流程中，只有滿足了屈服準則后才轉入下一次階計算流程．詳細的積分步驟在文獻[10]中給出．

2.3 子程序結構

基于本構模型公式和應力積分算法，參照ABAQUS用戶界面子程序的接口規范，進行了UINTER的程序設計及代碼編寫．代碼由 Fortran語言編寫．圖2中給UINTER的計算流程．文獻[10]中給出了具體的程序代碼．

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

3 數值驗證

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

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

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

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

數值計算與解析解計算結果對比情況如圖 4-6所示．圖 4-6中分別給出了循環剪切下的剪應力 τ與剪切位移w、法向應力σn與剪切位移w、剪應力τ與法向應力σn間的變化規律．從圖中可以看出，本文數值模擬結果與文獻解析解間有著很好的一致性，驗證了本子程序對界面循環加載性狀的模擬能力，說明本文UINTER成功實現了Mortara界面循環模型的開發．

圖4.循環剪切試驗中～wτ關系的解析解與數值解對比Fig. 4 Comparison between relationships of ～wτ obtained from analytical and numerical solution

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

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

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

4 結 語

Mortara界面循環模型是一個優秀的界面本構模型，能夠很好地對循環荷載作用下土-結構界面復雜的力學特性進行描述，包括應力-應變滯回特性、法向體積剪縮、剪切強度弱化等．通過本文研究，得到以下結論：

基于ABAQUS用戶子程序開發平臺UINTER，實現了對Mortara界面循環模型的二次開發．算例結果顯示，自開發的UINTER子程序具有很好的運算穩定性及精度，能夠模擬出土-結構界面的循環加載性狀．

利用 ABAQUS強大的非線性有限求解能力及優秀的前后處理界面，結合本文所開發的 Mortara界面循環模型的UINTER子程序，有望實現對循環荷載作用下土與結構相互作用問題數值分析．

子程序的成功開發極大擴展了該界面本構模型在巖土研究及設計中的應用，同時也為其他界面本構模型在ABAQUS中的開發提供了借鑒和參考．

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] 任宇. 長期豎向循環荷載作用下樁的累積變形特性試驗及理論研究[D]. 杭州: 浙江大學, 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.