摘    要：硬化土模型在描述軟土和較硬土的變形特性上有較好的表現(xiàn)，文章結(jié)合有限元軟件ABAQUS中的UMAT二次開發(fā)平臺(tái)，編寫了硬化土本構(gòu)模型子程序，提高了硬化土模型的泛用性，并提出了通過NewtonRapson迭代、Runge-Kutta迭代等數(shù)值方法求解任意應(yīng)變?cè)隽繉?duì)應(yīng)的應(yīng)力增量，最后通過室內(nèi)三軸壓縮試驗(yàn)數(shù)據(jù)驗(yàn)證了程序的正確性和合理性。

關(guān)鍵詞：硬化土模型;應(yīng)力更新算法;ABAQUS;二次開發(fā);

隨著現(xiàn)代巖土工程的發(fā)展，工程建設(shè)中遇到的問題逐漸從簡(jiǎn)單的穩(wěn)定性分析轉(zhuǎn)變?yōu)檩^精細(xì)的變形分析，能否精準(zhǔn)地進(jìn)行變形分析通常取決于計(jì)算使用的本構(gòu)模型[1]。由于巖土體復(fù)雜，盡管目前已提出了上百種本構(gòu)模型，但大多數(shù)模型僅能反映特定土體在特定情況下的力學(xué)行為，因此存在一定的局限性。巖土工程常用的Mohr-Coulomb模型和Drucker-Prager模型為理想彈塑性本構(gòu)模型，MCC模型為硬化彈塑性模型，難以同時(shí)反映土體的剪切硬化和壓縮硬化，采用Mohr-Coulomb強(qiáng)度理論作為屈服準(zhǔn)則，從Vermeer雙硬化模型發(fā)展而來的硬化土（HS）本構(gòu)模型[2]作為一種雙屈服面硬化彈塑性本構(gòu)模型，在描述軟土和較硬土的變形特性上有較好的表現(xiàn)[3]。

目前，除了PLAXIS、ZSoil等少數(shù)有限元軟件已嵌入HS模型，其他軟件使用該本構(gòu)仍需自行開發(fā)編寫相關(guān)程序。ABAQUS軟件在求解巖土等非線性問題上有突出的優(yōu)勢(shì)，有能為用戶提供編寫自定義本構(gòu)模型的二次開發(fā)平臺(tái)。徐遠(yuǎn)杰等[4]將Duncan-Chang本構(gòu)模型成功編成了UMAT子程序，岑威鈞和朱岳明[5]推導(dǎo)了平面應(yīng)變條件下UMAT子程序所需的彈塑性剛度矩陣，為后續(xù)學(xué)者開發(fā)UMAT子程序提供了支撐，使許多本構(gòu)模型被廣泛應(yīng)用于巖土工程數(shù)值模擬中。因此，為有效地?cái)U(kuò)展HS模型的應(yīng)用范圍，可選用ABAQUS作為HS模型的開發(fā)平臺(tái)。

為實(shí)現(xiàn)HS模型在ABAQUS上的應(yīng)用，文章利用ABAQUS用戶自定義材料子程序（UMAT）編寫HS本構(gòu)模型，優(yōu)化其應(yīng)力更新算法，并將模型的計(jì)算結(jié)果與室內(nèi)三軸壓縮試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，驗(yàn)證該本構(gòu)模型的合理性及可靠性，還分析了采用其他本構(gòu)時(shí)數(shù)值模擬的結(jié)果。

1 硬化土模型

硬化土模型中土體的剛度與應(yīng)力相關(guān)[4]。土體處于彈性階段時(shí)，采用雙剛度分別模擬加載與卸載時(shí)的力學(xué)特性，如圖1所示。土體處于塑性階段時(shí)，采用各向同性硬化法則與非關(guān)聯(lián)的流動(dòng)法則反映土體的剪脹性，同時(shí)引入帽蓋屈服面反映土體的壓縮硬化，形成了雙硬化本構(gòu)模型[5]。

圖1 硬化土模型彈性階段應(yīng)力-應(yīng)變關(guān)系

1.1 剪切屈服面

HS模型的剪切屈服面可用如下公式表示：

式中：Fs為剪切屈服函數(shù)；qa為極限偏應(yīng)力；E50為對(duì)應(yīng)50%強(qiáng)度時(shí)的割線模量；q為剪應(yīng)力；Eur為卸載再加載模量；γ p為塑性剪切應(yīng)變。

式中：E50ref為對(duì)應(yīng)參考圍壓σref時(shí)的E50模量；σ3為第三主應(yīng)力；c為黏聚力；φ為摩擦角；m為與土體性質(zhì)有關(guān)的冪指數(shù)。

式中：Eurref為對(duì)應(yīng)參考圍壓σref時(shí)的Eur模量。

式中：σ1為第一主應(yīng)力；f為處于破壞時(shí)的狀態(tài)；Rf為破壞比。

1.2 塑性勢(shì)函數(shù)

HS模型剪切屈服面采用的是不相適應(yīng)的流動(dòng)規(guī)則，其剪切塑性勢(shì)函數(shù)如下：

式中：Qs為剪切塑性勢(shì)函數(shù)；ψm為機(jī)動(dòng)剪脹角，由于不允許負(fù)剪脹角的存在，當(dāng)ψm<0時(shí)，取0。

式中：φm為機(jī)動(dòng)摩擦角；φcv為臨界摩擦角。

式中：φ為土體固有剪脹角。

1.3 壓縮屈服面

HS模型中為體現(xiàn)土體壓縮硬化特性加入了帽蓋屈服面，其屈服函數(shù)如下：

式中：F為壓縮屈服函數(shù)；δ為土體計(jì)算參數(shù)；σ2為第二主應(yīng)力；M為摩擦常數(shù)；Pc為前期固結(jié)應(yīng)力。

壓縮屈服面采用的是相適應(yīng)的流動(dòng)規(guī)則，其硬化定律如下：

式中：dp為硬化參數(shù)增量；Eoedref為切線模量；dεvp為塑性應(yīng)變?cè)隽俊?/p>

2 應(yīng)力更新計(jì)算方法優(yōu)化

ABAQUS中編寫自定義本構(gòu)模型需要推導(dǎo)出在笛卡爾三維（或二維）坐標(biāo)系下單元積分點(diǎn)的雅可比矩陣DDSDDE，即求解出該積分點(diǎn)各應(yīng)力增量對(duì)應(yīng)變?cè)隽康淖兓?σ/?ε，并據(jù)此對(duì)該增量步的總應(yīng)力和狀態(tài)變量進(jìn)行更新、輸出，通過接口返回ABAQUS主程序[6]。

2.1 變形階段劃分與試探應(yīng)力

對(duì)任意傳入的應(yīng)變?cè)隽浚燃俣ㄆ淙繛閺椥詰?yīng)變，由廣義胡克定律可得材料的彈性剛度矩陣：

式中：[De]為彈性剛度矩陣；λ為拉梅常數(shù)；G為剪切模量。

式中：E為彈性模量；μ為泊松比。

式中：{?σn+1}為第n+1步的應(yīng)力增量張量；{?σen+1}為第n+1步的彈性應(yīng)力增量張量。

式中：為試探應(yīng)力；為第n步的應(yīng)力張量；{?σn+1}為第n+1步的應(yīng)力張量。

將代入式（1）和式（9），若，則增量步全程處于彈性階段；若，則增量步從彈性階段過渡至彈塑性階段；若，則增量步全程處于彈塑性階段。

2.2 應(yīng)力更新計(jì)算

在彈塑性階段，材料在從到的過程中必然會(huì)經(jīng)過屈服應(yīng)力狀態(tài)，在從到此段內(nèi)只發(fā)生彈性變形，僅須用彈性剛度矩陣[De]更新應(yīng)力即可。在到這段則發(fā)生彈塑性變形，須求得彈塑性剛度矩陣[Dep]進(jìn)行求解。因此，為求解整個(gè)增量步的應(yīng)力變化，須確定在屈服面上的位置，可設(shè)：

式中：r為比例因子；為第n+1步的屈服應(yīng)力張量。

則：

由在屈服面上可得：

式中：為關(guān)于判斷應(yīng)力的函數(shù)；F (σn+r·?σn+1)為關(guān)于應(yīng)力及參數(shù)r的函數(shù)；F(r)為關(guān)于r的函數(shù)。

F(r)形式較復(fù)雜，很難求得準(zhǔn)確的解析解，因此可用Newton-Rapson迭代法求出數(shù)值解，鑒于Newton-Rapson迭代法對(duì)初值的要求較高，可先采用二分法大致確定初值范圍后再進(jìn)行迭代計(jì)算，迭代方程如下：

式中：rk+1為第k+1個(gè)參數(shù)r;rk為第k個(gè)參數(shù)r;F(rk)為關(guān)于rk的函數(shù)；F'(rk)為關(guān)于rk的函數(shù)的導(dǎo)函數(shù)。

得到r后，可將分成彈性增量{?σen+1}和彈塑性增量{?σepn+1}兩部分，通過彈性剛度矩陣反算得到真實(shí)彈性應(yīng)變?cè)隽縶?εe}：

式中：{?ε n+1}為第n+1步的應(yīng)變?cè)隽俊?/p>

由上式可得彈塑性應(yīng)變?cè)隽縶?ε ep}。

彈塑性階段的應(yīng)力更新求解式如下：

式中：?εep為塑性應(yīng)變?cè)隽浚粸閷?duì)塑性應(yīng)變求積分。

2.3 應(yīng)力更新迭代算法

彈塑性階段時(shí)，材料的彈塑性剛度矩陣是關(guān)于應(yīng)力的函數(shù)：

式中：[Dep]為材料的彈塑性剛度矩陣；σ為應(yīng)力。

在積分時(shí)，[Dep]會(huì)隨著應(yīng)力的變化而變化，積分的求解將十分困難。在ABAQUS中一般使用顯式歐拉方程進(jìn)行應(yīng)力更新，針對(duì)這種情況，要將步長(zhǎng)縮小，以保證程序計(jì)算的收斂性，極大地降低了計(jì)算效率。使用Runge-Kutta迭代法對(duì)式（21）進(jìn)行數(shù)值求解，能減少程序迭代計(jì)算的次數(shù)。因此文章使用Runge-Kutta法求解該常微分方程，迭代方程如下：

式中：σi+1為第i+1步時(shí)的應(yīng)力；σi為第i步時(shí)的應(yīng)力；Diep為不同σi時(shí)對(duì)應(yīng)的彈塑性剛度矩陣；h2為迭代步長(zhǎng)，通過公式計(jì)算可得或按需選取。

模型的彈塑性剛度矩陣為一般形式：

式中：A為塑性硬化模量，是硬化參數(shù)的函數(shù)；F、Q分別為屈服函數(shù)和塑性勢(shì)函數(shù)，采用相適應(yīng)流動(dòng)規(guī)則時(shí)兩者相等；?Q為對(duì)塑性勢(shì)函數(shù)求偏導(dǎo)；?F為對(duì)屈服函數(shù)求偏導(dǎo)；?σ為對(duì)應(yīng)力求偏導(dǎo)；T為線性代數(shù)矩陣轉(zhuǎn)置符號(hào)。

3 HS模型在程序中的實(shí)現(xiàn)

在計(jì)算中，ABAQUS主程序會(huì)調(diào)用用戶編寫的UMAT子程序，將主程序各單元積分點(diǎn)上的應(yīng)力、應(yīng)變傳遞給子程序[7]。子程序通過本構(gòu)關(guān)系計(jì)算得到雅可比矩陣及迭代更新后的應(yīng)力、應(yīng)變和用戶設(shè)置的狀態(tài)變量，并將相應(yīng)數(shù)據(jù)回傳給ABAQUS主程序。計(jì)算過程中，每一增量步的各迭代步都需調(diào)用UMAT子程序，直至全部增量步計(jì)算完成[8]。

一般來說，ABAQUS中UMAT子程序的編寫需要采用Fixed Format，但通過修改ABAQUS配置文件“win86＿64.env”可將UMAT編寫格式修改為Free Format，簡(jiǎn)化UMAT子程序的編寫流程，修改后原UMAT接口的固定格式也需相應(yīng)修改[9],UMAT子程序接口固定格式示意圖如圖2所示。

圖2 UMAT子程序接口固定格式示意圖

采用模塊化程序設(shè)計(jì)理念編寫UMAT子程序，即將需要重復(fù)使用的代碼功能整合成獨(dú)立運(yùn)行的小模塊，當(dāng)有需要時(shí)再調(diào)用，顯著提高了程序代碼的可讀性及修改效率[10]。流程圖如圖3所示。

4 模型驗(yàn)證

4.1 算例簡(jiǎn)介

為驗(yàn)證編寫的UMAT子程序的準(zhǔn)確性和可靠性，在ABAQUS上使用UMAT子程序?qū)θS壓縮試驗(yàn)?zāi)Ｐ瓦M(jìn)行數(shù)值模擬，并與室內(nèi)三軸試驗(yàn)結(jié)果進(jìn)行對(duì)比。數(shù)值模型為直徑39.1 mm、高80 mm的標(biāo)準(zhǔn)三軸試樣。模型邊界條件如下：底面設(shè)置法向位移約束；側(cè)面根據(jù)不同工況施加不同的徑向壓力（125 kPa、225 kPa和300 kPa）；頂面施加分級(jí)荷載，先施加與圍壓相等的固結(jié)壓力模擬試樣的固結(jié)過程，再逐級(jí)加載至試樣破壞[11]。模擬材料為花崗巖殘積土，土體相應(yīng)的參數(shù)由三軸試驗(yàn)得到，具體數(shù)值如表1所示。

圖3 HS模型UMAT子程序流程圖

4.2 結(jié)果對(duì)比

將數(shù)值模型計(jì)算結(jié)果與室內(nèi)三軸試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，具體結(jié)果如圖4所示。從圖4可以看出，花崗巖殘積土在應(yīng)力軟化前的應(yīng)變-應(yīng)力關(guān)系類似于雙曲線，與HS模型彈性階段變化趨勢(shì)相吻合[12]。在不同圍壓尤其是125 kPa和225 kPa下，數(shù)值模型計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)均一致，充分表明了文章開發(fā)的HS模型子程序的可靠性和準(zhǔn)確性[13]。

圖4 不同圍壓HS本構(gòu)數(shù)值模擬結(jié)果與試驗(yàn)數(shù)據(jù)對(duì)比

表1 土體計(jì)算參數(shù)

為了對(duì)比HS模型與現(xiàn)有本構(gòu)模型，文章將開發(fā)的HS模型與ABAQUS中的Mohr-Coulomb模型進(jìn)行比較，將數(shù)值結(jié)果與試驗(yàn)實(shí)測(cè)數(shù)據(jù)進(jìn)行對(duì)比，結(jié)果如圖5和圖6所示。從圖5、圖6可以看出，盡管Mohr-Coulomb模型能在一定程度上表現(xiàn)土體的應(yīng)力-應(yīng)變行為，但與HS模型相比仍有較大差距，HS模型既能體現(xiàn)土體變形前段的雙曲線應(yīng)力-應(yīng)變關(guān)系，又能很好地反映土體剪切屈服平臺(tái)，這充分表明了在ABAQUS上進(jìn)行HS模型開發(fā)的必要性[14]。

圖5 不同本構(gòu)模型三軸試驗(yàn)數(shù)值模擬結(jié)果對(duì)比（σ2=200 k Pa) 

圖6 不同本構(gòu)模型三軸試驗(yàn)數(shù)值模擬結(jié)果對(duì)比（σ2=300 k Pa)  

5 結(jié)論

文章利用ABAQUS的二次開發(fā)平臺(tái)完成了硬化土模型的數(shù)值實(shí)現(xiàn)，并成功將其應(yīng)用于三軸試驗(yàn)相關(guān)的數(shù)值模擬計(jì)算，得到了以下成果：

(1）推導(dǎo)出HS模型在不同應(yīng)力狀態(tài)下的應(yīng)力更新計(jì)算公式，提出了對(duì)任意給定應(yīng)變?cè)隽烤赏ㄟ^Runge-Kutta迭代等數(shù)值方法求解相應(yīng)應(yīng)力增量的算法。

(2）通過三軸壓縮模擬實(shí)驗(yàn)驗(yàn)證了HS模型UAMT子程序開發(fā)的可靠性和必要性，并通過與其他本構(gòu)模型對(duì)比，證明了HS模型能很好地模擬土體的應(yīng)力應(yīng)變關(guān)系。

(3）基于Runge-Kutta算法及Newton-Rapson迭代法對(duì)UMAT進(jìn)行優(yōu)化的算法設(shè)計(jì)流程及思路可以為其他本構(gòu)模型的開發(fā)提供借鑒，進(jìn)一步促進(jìn)巖土工程數(shù)值分析技術(shù)的發(fā)展。

參考文獻(xiàn)

[1] 岑威鈞,陳司寧,鄧同春,等.土石料雙屈服面彈塑性模型的二次開發(fā)算法與應(yīng)用[J].西南交通大學(xué)學(xué)報(bào),2018,53(3):582-588.

[2] 王仁輝.硬化土模型在Open Sees中的實(shí)現(xiàn)與應(yīng)用[D].開封:河南大學(xué),2021.

[3] 宗露丹,徐中華,翁其平,等.小應(yīng)變本構(gòu)模型在超深大基坑分析中的應(yīng)用[J].地下空間與工程學(xué)報(bào),2019,15(S1):231-242.

[4] 徐遠(yuǎn)杰,王觀琪,李健,等.在ABAQUS中開發(fā)實(shí)現(xiàn)Duncan-Chang本構(gòu)模型[J].巖土力學(xué),2004(7):1032-1036.

[5] 岑威鈞,朱岳明.基于ABAQUS的土石料本構(gòu)模型二次開發(fā)及其應(yīng)用[J].水利水電科技進(jìn)展,2005(6):78-81.

[6] 王正振,龔維明,戴國(guó)亮,等.考慮位移影響的土壓力非線性計(jì)算[J].巖土工程學(xué)報(bào),2019,41(S2):244-248.

[7] 王祥秋,楊柱,鄭土永.珠三角典型軟土硬化土模型及其工程應(yīng)用研究[J].山東理工大學(xué)學(xué)報(bào)(自然科學(xué)版),2022,36(1):19-26,32.

[8] 董正方,過晴,王仁輝,等.硬化土模型在OpenSees中的實(shí)現(xiàn)[J].中國(guó)科技論文,2023,18(2):193-203.

[9] 林德周.小應(yīng)變土體硬化模型參數(shù)試驗(yàn)研究及工程應(yīng)用[D].杭州:浙江大學(xué),2023.

[10] 劉大維.基于GA-BP的小應(yīng)變硬化土本構(gòu)模型的參數(shù)反演及在基坑工程中的應(yīng)用研究[D].西安:長(zhǎng)安大學(xué),2023.

[11] 姜兆華.基坑開挖時(shí)鄰近既有隧道的力學(xué)響應(yīng)規(guī)律研究[D].重慶:重慶大學(xué),2014.

[12] 湯道飛,王長(zhǎng)虹.小應(yīng)變硬化土本構(gòu)模型在FLAC3D中的開發(fā)及應(yīng)用[J].上海大學(xué)學(xué)報(bào)(自然科學(xué)版),2023,29(3):549-561.

[13] 姜兆華,張永興.硬化土模型在FLAC3D中的二次開發(fā)[J].解放軍理工大學(xué)學(xué)報(bào)(自然科學(xué)版),2013,14(5):524-529.

[14] 鄭土永.基于HS本構(gòu)模型軟土地鐵換乘車站深基坑力學(xué)特性研究[D].佛山:佛山科學(xué)技術(shù)學(xué)院,2022.

文章來源：工程技術(shù)研究