珠三角典型軟土硬化土模型及其工程應用研究

王祥秋，楊柱，鄭土永

(1.佛山科學技術學院 交通與土木建筑學院，廣東 佛山 528000;2.北京市市政工程設計研究總院，北京 100037)

目前，有限元計算方法在巖土工程中得到了廣泛的應用，而本構(gòu)模型的選取及模型參數(shù)的確定是進行數(shù)值計算的關鍵。宋二祥等[1]，王海波等[2]、王衛(wèi)東等[3-4]將硬化土本構(gòu)模型應用于深基坑工程開挖數(shù)值模擬；龐小朝等[5]、楊蘭強等[6]通過室內(nèi)試驗獲得了巖土體硬化土本構(gòu)模型參數(shù)，分析了硬化土本構(gòu)模型的可行性； Sukkarak等[7]利用硬化土本構(gòu)模型分析了大壩地基沉降規(guī)律；溫科偉等[8]、沈璽[9]基于硬化土本構(gòu)模型對地鐵隧道的變形特性進行分析；謝東武等[10]通過室內(nèi)三軸試驗確定了上海軟土小應變硬化土模型參數(shù)，并基于土單元數(shù)值模擬對模型參數(shù)的敏感性進行分析。研究結(jié)果表明，硬化土彈塑性本構(gòu)模型能綜合考慮黏性土的剪脹性和中性加載，能區(qū)分加荷和卸荷特性，對復雜環(huán)境下巖土與地下工程開挖數(shù)值分析具有很好的適用性。而目前內(nèi)嵌硬化土本構(gòu)模型(HS模型)的巖土分析軟件只有Midas GTS、Plaxis和Zsoil等，其它巖土分析軟件如FLAC3D、ABAQUS等有限元分析平臺尚未內(nèi)嵌硬化土本構(gòu)模型，在一定程度上影響了相關軟件在巖土工程領域的廣泛使用。為此，國內(nèi)外學者針對相關軟件開展了HS模型的二次開發(fā)研究，如姜兆華等[11]、王春波等[12]基于VC++編程環(huán)境，利用FLAC3D提供的二次開發(fā)平臺編制HS模型的有限差分程序，實現(xiàn)了HS模型在FLAC3D的二次開發(fā)，滿足了巖土工程領域大變形問題計算的需要。而ABAQUS作為具有強大非線性問題求解功能的有限元分析軟件，目前尚未實現(xiàn)硬化土模型內(nèi)嵌功能，開發(fā)研制基于ABAQUS平臺的HS模型分析功能則可滿足巖土工程領域非線性小變形問題計算的需要。

為此,本文基于ABAQUS軟件強大的自定義材料模型開發(fā)功能，通過推導硬化土本構(gòu)關鍵算法，基于Fortran編程環(huán)境，利用ABAQUS提供的用戶材料UMAT子程序接口，實現(xiàn)HS本構(gòu)模型的二次開發(fā)。利用GDS多功能應力路徑三軸儀對原狀土進行固結(jié)排水剪切實驗、三軸固結(jié)排水加載-卸載-再加載實驗以及三聯(lián)固結(jié)實驗測定珠三角地區(qū)典型軟土HS模型參數(shù)，并通過實驗數(shù)據(jù)與數(shù)值計算結(jié)果對比分析及深基坑工程實例監(jiān)測結(jié)果與模擬結(jié)果比較，驗證HS本構(gòu)模型UMAT子程序的合理性及可靠性。

1 硬化土本構(gòu)模型

硬化本構(gòu)模型(即HS模型)由Schanz等[13]在VERMEER雙硬化模型基礎上提出，該模型采用 Mohr-Coulomb破壞準則。土體彈性階段采用雙剛度分別考慮土體加載與卸載影響，并假定土體剛度與應力水平相關，豎向應變ε1與偏應力q之間滿足雙曲線關系[4-5]，即

(1)

其中

(2)

假定土體塑性階段屈服面由剪切屈服面和壓縮屈服面兩部分組成。剪切屈服面(剪切屈服函數(shù)Fs)定義為

(3)

(4)

其中

(5)

若采用關聯(lián)流動法則,則剪切屈服面塑性勢函數(shù)可表示為

(6)

(7)

式中：ψm為機動剪脹角；φm為機動摩擦角；φcv為臨界摩擦角；ψ為土體固有剪脹角。

壓縮屈服面(壓縮屈服函數(shù)Fv)定義為

(8)

其中

(9)

(10)

若采用關聯(lián)流動法則,則壓縮屈服面塑性勢函數(shù)可表示為

(11)

2 二次開發(fā)關鍵算法

基于ABAQUS軟件的二次開發(fā)中最關鍵的問題是如何選擇本構(gòu)模型的積分算法。常用的積分算法主要包括隱式積分算法和顯式積分算法，其中，顯式積分算法又包括基本剛度法、中點剛度法以及SLOAN提出的帶誤差控制的修正向后EULER返回算法。本文主要采用基本剛度法對硬化土本構(gòu)模型進行二次開發(fā)，其基本步驟如下：

1)計算彈性預測應力。在第n步應力{σn}已知的情況下，先按虎克定律預測第n+1步的試探應力{σn+1,trial}。

{σn+1,trial}={σn}+De{Δεn+1}，

(12)

式中：De為彈性剛度矩陣；{Δεn+1}為應變增量矩陣。

2)屈服準則判斷。將試探應力{σn+1,trial}分別代入式(3)和式(8)。如果Fv({σn+1,trial}，{kv})<0，F(xiàn)s({σn+1,trial}，{ks})<0,則表明材料此時處于彈性階段，硬化參數(shù){kv}和{ks}保持不變，剛度矩陣為彈性剛度矩陣De，上述試探應力即為第n+1個增量步的應力；如果Fv({σn+1,trial}，{kv})≥0或Fs({σn+1,trial}，{ks})≥0,則表明材料處于單屈服階段，可能發(fā)生剪切屈服或者壓縮屈服，此時應力增量按式(13)確定，剛度矩陣按式(14)確定；如果Fv({σn+1,trial}，{kv})≥0且Fs({σn+1,trial}，{ks})≥0,則表明材料同時發(fā)生剪切屈服和壓縮屈服，此時應力增量按式(13)確定，剛度矩陣按式(29)確定。

3)求解增量本構(gòu)方程對應的彈塑性矩陣及應力增量。根據(jù)塑性理論，則有彈塑性本構(gòu)關系的增量本構(gòu)方程的一般形式為

{Δσn+1}=Dep{Δεn+1}，

(13)

式中Dep為彈塑性剛度矩陣。

對于只發(fā)生剪切屈服或者壓縮屈服，則有增量本構(gòu)關系的彈塑性剛度矩陣為

(14)

對于雙屈服面模型剛度矩陣，部分做法是先計算柔度矩陣，再通過求逆得到剛度矩陣，但是求逆往往需要耗費大量計算時間[14]。本文采用David等[15]人的方法，直接推導雙屈服面模型彈塑性剛度矩陣,這種方法同樣可以適用于三屈服面模型。

當同時發(fā)生剪切屈服以及壓縮屈服時，則總應變增量主要由彈性應變增量{Δεe}、剪切應變增量{Δεps}以及壓縮應變增量{Δεpv}組成，即

{Δε}={Δεe}+{Δεps}+{Δεpv}，

(15)

而根據(jù)彈塑性理論，總應力增量可以表示為

{Δσ}=De{Δεe}

(16)

或者

{Δεe}=De-1{Δσ}。

(17)

將式(15)代入式(16)可得

{Δσ}=De({Δε}-{Δεps}-{Δεpv})，

(18)

其中，剪切塑性應變和壓縮塑性應變分別與剪切屈服函數(shù)及壓縮屈服函數(shù)對應的塑性勢函數(shù)相關，且

(19)

式中：λs和λv為塑性比例因子。

將式(19)代入式(18)則有

(20)

又由一致性條件得

(21)

將式(20)代入式(21)得

{

(22)

其中

(23)

式(22)也可改寫為

(24)

其中：

(25)

(26)

(27)

聯(lián)立式(24)解得

(28)

代入式(20)得

(29)

其中：

Ω=LssLvv-LsvLvs，

(30)

(31)

(32)

4)對屈服函數(shù)和勢函數(shù)求一階導數(shù)。為了便于數(shù)值計算，根據(jù)彈塑性理論，將硬化土本構(gòu)模型的屈服函數(shù)和勢函數(shù)表示為應力不變量的形式，其中：I1為第一應力不變量；J2為第二偏應力不變量；J3為第三偏應力不變量；θ為應力洛德角。則剪切屈服函數(shù)可以表達為

(33)

壓縮屈服函數(shù)可以表達為

Fv=

(34)

與剪切屈服函數(shù)相對應的塑性勢函數(shù)為

(35)

屈服函數(shù)流動矢量可以表示為

(36)

式中

(37)

對于剪切屈服函數(shù)，式(36)中參量C1、C2、C3按下式計算：

(38)

對于壓縮屈服函數(shù)，式(36)中參量C1、C2、C3則按下式計算：

(39)

與屈服函數(shù)流動矢量計算過程一樣，塑性勢函數(shù)流動矢量的表達式如下：

D1{a1}+D2{a2}+D3{a3}。

(40)

對于剪切屈服函數(shù)對應的塑性勢函數(shù)，式中參量D1、D2、D3按下式計算：

(41)

5)更新應力、應變和硬化參數(shù)。

{σn+1}={σn}+{Δσn+1}。

(42)

當同時發(fā)生剪切應變和壓縮應變時，第n+1步結(jié)束時應變更新為

(43)

當僅發(fā)生剪切應變或壓縮應變時，第n+1步結(jié)束時應變更新為

(44)

或者

(45)

剪切屈服硬化參數(shù)增量為

(46)

第n+1步結(jié)束時，剪切屈服硬化參量更新為

(47)

對于壓縮屈服，塑性體應變增量為

(48)

壓縮屈服硬化參數(shù)增量為

(49)

第n+1步結(jié)束時壓縮屈服硬化參量更新為

pc,n+1=pc,n+Δpc,n+1=

(50)

3 子程序開發(fā)關鍵技術

根據(jù)上述硬化土本構(gòu)模型關鍵算法，基于ABAQUS提供的UMAT子程序接口,利用Fortran編程語言環(huán)境可開發(fā)HS模型的UMAT子程序。子程序開發(fā)中應注意如下關鍵技術問題:

1)子程序編寫應符合HS模型和彈塑性有限元計算特點，同時其輸入輸出格式及變量名應注意與標準ABAQUS程序一致。

2)ABAQUS分析軟件中，應力應變符號以拉為正、壓為負，主應力、主應變排序與巖土力學有關規(guī)定相反。

3)程序編寫過程可以輸入write(7，*),從而可在保存于文件夾內(nèi)的msg文件中輸出所關心的變量并加以考察。

4)UMAT開始計算時，應避免ABAQUS主程序傳入初始應力為零而導致計算過程出現(xiàn)極大值現(xiàn)象。

UMAT子程序調(diào)用與求解流程如圖1所示。

圖1 UMAT子程序開發(fā)與工作流程圖

1)在ABAQUS程序求解時，每一個增量加載步開始時，ABAQUS主程序都會在單元積分點上調(diào)用UMAT子程序，傳入當前狀態(tài)的總應力、應變增量和用戶自定義狀態(tài)變量等基本信息，同時傳入主程序計算得出的應變增量。

3)變量更新值通過接口返回主程序，雅可比矩陣將同單元應變矩陣運算形成單元剛度矩陣，進而獲得總體剛度矩陣；主程序根據(jù)當前荷載增量求解位移增量并進行平衡校核；如果不滿足用戶指定的誤差或者缺省值，ABAQUS將進行迭代，直到滿足收斂條件為止，然后進行下一增量步的求解。

4 工程實例

4.1 工程概況

珠三角某地鐵車站采用地下三層島式結(jié)構(gòu)，基坑全長211.4 m，標準段寬20.9 m，車站基坑開挖深度為24.22～26.72 m，采用1 000 mm厚的地下連續(xù)墻作為圍護結(jié)構(gòu)，地連墻嵌固深度為6 m?？紤]到車站周邊環(huán)境條件復雜，車站深基坑工程采用明挖法施工，與地鐵1號線換乘節(jié)點處采用暗挖法施工?；舆吘€西北側(cè)距離某廣場地下室邊線僅2.48 m。自上而下設四道支撐和角撐，分別位于-1.7 m、-7.3 m、-14.1 m、-19.8 m和-24.2 m。豎向標準段第一道和第三道采用砼支撐，第二道和第四道采用φ609鋼支撐，并分別施加400 kN和600 kN預應力，換乘節(jié)點段和端頭井采用四道砼支撐，基坑支護結(jié)構(gòu)平面布置如圖2所示。

圖2 車站基坑支護結(jié)構(gòu)平面圖

4.2 三維有限元分析模型

采用ABAQUS大型有限元分析軟件對車站深基坑工程力學特性進行分析。為消除邊界效應影響，根據(jù)實際基坑的開挖深度及平面尺寸建立的模型空間尺寸為350 m×200 m×60 m(長×寬×高)。土體采用實體8節(jié)點減縮單元模擬，單元類型為C3D8R。車站結(jié)構(gòu)及連續(xù)墻變形均采用實體8節(jié)點協(xié)調(diào)單元模擬，單元類型為C3D8I；砼支撐、鋼支撐及立柱采用梁單元(Beam)，單元類型為B31，計算模型如圖3所示。地連墻與土體的作用、 車站結(jié)構(gòu)與土體的作用及地連墻與車站結(jié)構(gòu)的作用均選擇tie連接，只傳遞拉力和壓力且不產(chǎn)生相對位移。對模型整體施加重力荷載并限制模型側(cè)向位移及底部豎向位移。

(a)基坑整體分析模型

4.3 計算參數(shù)

根據(jù)巖土工程勘察報告，利用GDS多功能應力路徑三軸儀對原狀土進行固結(jié)排水剪切實驗、三軸固結(jié)排水加載-卸載-再加載試驗以及利用三聯(lián)固結(jié)儀進行常規(guī)固結(jié)實驗測定各土層HS模型參數(shù)，具體參數(shù)見表1。

表1 車站深基坑土層主要物理力學性能參數(shù)

基坑支護結(jié)構(gòu)(含地下連續(xù)墻、內(nèi)支撐)采用線彈性模型進行模擬，地下連續(xù)墻、混凝土支撐均采用C35混凝土，考慮到施工因素影響，其彈性模量按80%折減取為24 GPa，泊松比取0.2；鋼管支撐和鋼圍檁彈性模量取為210 GPa，泊松比為0.3。為了客觀地模擬地鐵車站深基坑施工開挖力學特性，根據(jù)深基坑支護結(jié)構(gòu)設計方案以及實際開挖工況,結(jié)合車站換乘節(jié)點與3號線地鐵車站深基坑開挖先后次序安排，確定有限元模擬的施工開挖步，計算工況見表2。

表2 計算工況

4.4 模擬結(jié)果與監(jiān)測數(shù)據(jù)對比分析

為了對比分析基于HS模型的有限元數(shù)值結(jié)果與現(xiàn)場監(jiān)測成果的吻合程度，任選基坑長邊典型監(jiān)測點ZQT4和ZQT14作為分析對象，如圖4所示。

圖4 車站基坑監(jiān)測點平面圖

由現(xiàn)場監(jiān)測和有限元分析可得監(jiān)測點ZQT4和ZQT14處地下連續(xù)墻深部位移在基坑關鍵施工步即第4施工步(開挖第一層土體)、第6施工步(開挖第三層土體)及第9施工步(開挖第六層土體)時的變化曲線(如圖5所示)。

由圖5可知，在車站深基坑開挖過程中，地下連續(xù)墻深部位移監(jiān)測值與計算值的變形趨勢基本吻合。在基坑開挖初期(第4步)，兩者誤差較??；隨著基坑開挖深度不斷增加，地下連續(xù)墻深部位移拐點以上部分的監(jiān)測值與計算值仍然具有較高的吻合度， 當基坑開挖到底時(第9步)， 監(jiān)測點ZQT4處地下連續(xù)墻深部位移監(jiān)測最大值為24.5 mm,數(shù)值模擬深部位移最大值為22.3 mm；監(jiān)測點ZQT14處地下連續(xù)墻深部位移監(jiān)測最大值為23.8 mm,數(shù)值模擬深部位移最大值為21.4 mm，兩者最大誤差僅為11.2%。但在拐點以下部分，其監(jiān)測值與有限元模擬值的誤差逐漸增大，其原因可能與軟土深基坑工程隨著開挖深度不斷增大，基坑開挖時空效應愈來愈明顯，以及由于基坑施工過程未能及時施加鋼支撐等因素有關。而基于MC模型得出的墻體深部位移計算值較HS模型以及現(xiàn)場監(jiān)測值偏小，主要原因在于MC模型只定義了一個彈性模量，未能考慮土體加載和卸載模量的差異性，且無法考慮應力路徑的影響，導致基坑開挖產(chǎn)生較大的坑底回彈，從而減小了墻體的變形。

(a)特征點1(ZQT4)

5 結(jié)論

1)深基坑工程實例分析結(jié)果表明，基于顯式積分基本剛度法構(gòu)建的硬化土本構(gòu)模型關鍵算法以及基于ABAQUS有限元分析平臺開發(fā)研制的HS模型UMAT子程序是合理可行的。

2)三維有限元數(shù)值模擬結(jié)果表明，硬化土本構(gòu)模型(HS模型)能有效模擬土體的硬化特性，而莫爾-庫侖模型(MC模型)只能模擬土體理想彈塑性變形；因此，與莫爾-庫侖模型(MC模型)相比，HS模型能更好地模擬軟土非線性應力應變特性。

3)基于ABAQUS非線性有限元分析平臺研發(fā)的HS模型UMAT子程序，可推廣應用于軟土地下工程復雜施工力學性態(tài)的分析模擬。如能綜合考慮軟土深基坑工程時空效應等因素影響，將會進一步提高有限元模擬效果。