砂土狀態相關本構模型在ABAQUS中的數值實現方法

郭浩天 孫增春 梁放 崔昊 汪成貴 肖楊

DOI：?10.11835/j.issn.2096-6717.2023.081

收稿日期：2023?04?22

基金項目：國家自然科學基金（52078085）；中央高校科研基本業務費（2452022121）

作者簡介：郭浩天（1998-?），男，主要從事土體基本力學特性研究，E-mail：guohaotian20@163.com。

通信作者：孫增春（通信作者），男，博士，E-mail：cqsunzcu@163.com。

Numerical implementation method of state-dependent constitutive model of sand in ABAQUS

GUO Haotian¹，?SUN Zengchun²，?LIANG Fang¹，?CUI Hao¹，?WANG Chenggui¹，?XIAO Yang¹

（1. School of Civil Engineering， Chongqing University， Chongqing 400045， P. R. China；?2. College of Water Resources and Architectural Engineering， Northwest A&F University， Yangling 712100， Shaanxi， P. R. China）

土體的剪脹、應變硬化和應變軟化等應力-應變特性與土體當前所處的狀態相關，狀態參數（當前孔隙比與臨界孔隙比的差值）被廣泛用于表征不同狀態下砂土的力學行為^[1]。隨后，不同學者相繼提出了狀態孔隙比參數和、狀態壓力參數，以及狀態孔隙比-壓力參數^[2]。考慮到這一特性，可將狀態參數引入廣義塑性力學或邊界面塑性力學中，建立一系列適用于砂土的狀態相關塑性本構模型。

巖土工程問題的數值分析中，土體的力學響應受本構模型的應力-應變關系控制。在有限元框架內進行數值建模時，首先要在應力點水平上研究精確、高效的本構積分算法。當將復雜、高級的彈塑性本構模型應用到有限元程序中時，應力狀態需要在每個高斯點上使用積分算法進行多次更新，積分算法的性能直接影響數值計算的精確性、收斂性和穩定性。數值積分算法主要分為完全隱式、半隱式和顯式積分算法3類。隱式算法通常需要進行局部迭代，更新后的應力自動滿足屈服條件，應力點由彈性轉變為塑性狀態時，不需要對應力路徑與屈服面的交點進行局部定位，具有計算精度高、全局收斂速度快的優點。但一致性切線模量以及塑性勢面二階導數的推導較為復雜。與隱式算法不同，半隱式算法在塑性修正階段用已知應力點塑性勢函數的梯度確定映射回歸的方向，避免了二階求導的同時具有二階收斂速度，降低了數值計算的難度。但是其精確性和收斂性遠低于隱式算法，且缺乏一致的線性化。顯式算法通過計算塑性勢函數或屈服面的梯度，結合彈塑性剛度矩陣來預測當前增量步結束時的應力增量。由于不需要進行局部迭代，與隱式和半隱式算法相比，顯式算法具有更高的計算效率。此外，通過控制應力和狀態變量的局部誤差，自動劃分子增量步的大小能夠顯著提高顯式算法的精確性和收斂性^[3]。

采用泰勒級數展開構造的Runge-Kutta方案是一種高精度的標準數值方法，在高斯點水平的積分本構關系上得到了廣泛的應用。對于任意高階（n階）Runge-Kutta方案，計算階段k處的應力增量可視為執行計算的加權和，用式（1）表示。

式中：和分別為第k階段的試算應力和相對應的彈塑性剛度矩陣；為第k階段的應變增量；為與Runge-Kutta方案相關的系數。

數值積分的相對誤差可通過高階項與低階項的差值確定，表示為

式中：和為與Runge-Kutta方案相關的系數。

狀態相關本構模型的彈塑性剛度矩陣表示為

式中：和分別為塑性流動方向和加載方向向量；為彈性剛度矩陣；為塑性模量。

將五階Runge-Kutta -Dormand-Prince（RKDP5）方案與自動誤差控制的顯式積分算法相結合，在ABAQUS中采用用戶自定義材料子程序（UMAT）對砂土狀態相關本構模型進行數值實現，模型的具體表達式和參數可參考相關文獻^[4-5]。圖1為顯式Runge-Kutta應力積分算法的流程圖。

圖2為數值模擬與試驗數據的對比，可以看出，通過UMAT子程序數值計算得到的應力-應變曲線與試驗數據基本吻合，驗證了帶自動誤差控制的Runge-Kutta算法在ABAQUS中開發砂土狀態相關本構模型具有可行性。在顯式積分算法中，應變增量步大小和誤差容許值STOL會影響積分算法的精度、效率和穩定性。圖3為RKDP5算法計算獲得的局部誤差隨應變增量步和誤差容許值STOL的變化規律。從圖3可以看出，應變增量步長對局部誤差的影響隨STOL值的減小逐漸降低。類似地，STOL對局部誤差的影響隨著應變增量步的增加也在逐漸減弱。對于STOL的典型值（10^-3～10^-4），不同應變增量步下的局部誤差在10^-3～10^-5范圍內。

地基的承載特性是驗證本構模型有效性和應力積分算法可靠性的典型問題^[6]。加載過程中，基礎下土體的應變增量值變化范圍非常大，這要求在所有高斯點上采用精確和穩定的應力積分方案。由于問題的對稱性，選取四分之一的地基進行模擬。在ABAQUS中有限元模型采用C3D8R單元（21 025個單元，23 400個節點），基礎頂部采用位移控制方式進行加載。圖4為三維方形基礎承載力分析結果。從圖4（a）可以觀察到，在豎向加載過程中，基礎下方的土體向下移動形成楔形體，基礎附近的土體逐漸向上隆起。塑性區出現在基礎的底部，并逐漸向地面延伸形成兩條明顯的剪切帶。模型計算得到的荷載隨著位移的增加逐漸增加，直到達到一個穩定值，如圖4（c）所示。通過上述綜合分析表明，自動誤差控制的顯式Runge-Kutta算法不僅可用于單元尺度試驗的分析，而且在多單元的邊值問題的數值分析中也是可行的，進一步驗證了砂土狀態相關本構模型在ABAQUS應用中的可行性和有效性。

參考文獻

[1] ?吉華， 孫逸飛. 堆石料的三維應力分數階本構模型[J]. 土木與環境工程學報（中英文）， 2022， 44（4）： 27-34.

JI H， SUN Y F. Three-dimensional stress-fractional constitutive model for rockfill [J]. Journal of Civil and Environmental Engineering， 2022， 44（4）： 27-34. （in Chinese）

[2] ?XIAO Y， LIU H L， CHEN Y M， et al. Bounding surface plasticity model incorporating the state pressure index for rockfill materials [J]. Journal of Engineering Mechanics， 2014， 140（11）： 04014087.

[3] ?SLOAN S W， ABBO A J， SHENG D C. Refined explicit integration of elastoplastic models with automatic error control [J]. Engineering Computations， 2001， 18（1/2）： 121-194.

[4] ?SUN Z C， CHU J A， XIAO Y. Formulation and implementation of an elastoplastic constitutive model for sand-fines mixtures [J]. International Journal for Numerical and Analytical Methods in Geomechanics， 2021， 45（18）： 2682-2708.

[5] ?孫增春. 溫度效應下砂-黏土混合物的力學特性與本構模型研究[D]. 重慶： 重慶大學， 2022.

SUN Z C. Study on mechanical properties and constitutive model of sand-clay mixture under temperature effect [D]. Chongqing： Chongqing University， 2022. （in Chinese）

[6] ?郭毓熙， 章懿濤， 方祥位， 等. 不同含水和密實狀態下珊瑚砂地基承載特性試驗研究[J]. 土木與環境工程學報（中英文）， 2023， 45（5）： 49-57.

GUO Y X， ZHANG Y T， FANG X W， et al. Experimental study on bearing characteristics of coral sand foundation under different water content states and relative compactnesses [J]. Journal of Civil and Environmental Engineering， 2023， 45（5）： 49-57. （in Chinese）

（編輯??胡英奎）