陳 成， 周正明， 張先偉

(中國科學(xué)院 武漢巖土力學(xué)研究所 巖土力學(xué)與工程國家重點實驗室，武漢 430071)

交通荷載作用下軟土地基的長期服役性能是軟土工程中的一個重要問題。 長期循環(huán)荷載作用下地基土的沉降變形隨時間逐漸累加，殘余變形受諸多因素影響，如循環(huán)荷載水平，頻率以及地基土的排水條件等[1]。

目前，計算循環(huán)荷載作用下土體累積殘余變形的模型可以分為兩種。一種是建立較為復(fù)雜的本構(gòu)模型，其大致可分為兩類，一類是黏彈性模型，如Hardin-Drnevich 模型[2]和Iwan模型[3]等；一類是彈塑性模型，如運動硬化嵌套面本構(gòu)模型[4]以及邊界面模型[5]等。黏彈性模型數(shù)學(xué)表述相對簡單便于程序化，且計算量較小，在工程實踐中應(yīng)用較多。但此類模型不能反映應(yīng)變軟化、應(yīng)力路徑的影響。彈塑性模型由于其理論嚴(yán)密，可較真實地反映土體在振動荷載作用下每一個循環(huán)內(nèi)的變形特征以及考慮應(yīng)變軟化、應(yīng)力路徑等。但模型計算中一般采用傳統(tǒng)的小步長積分方法，因而會產(chǎn)生巨大的計算量，計算效率方面有待進(jìn)一步的研究。另外一種計算方法是首先在試驗或?qū)崪y資料基礎(chǔ)上建立土體累積應(yīng)變與初始固結(jié)狀態(tài)、靜應(yīng)力狀態(tài)、動應(yīng)力及循環(huán)次數(shù)的關(guān)系的擬合曲線，再結(jié)合分層總和法等方法來預(yù)測土體在動荷載作用下的沉降。如Monismith等[6]建立的累積塑性應(yīng)變與循環(huán)次數(shù)間的指數(shù)模型，以及在其基礎(chǔ)上考慮靜強度參數(shù)、初始靜偏應(yīng)力等因素建立的修正模型[7]。黃茂松等[8-9]基于臨界狀態(tài)理論，提出了“動偏應(yīng)力水平”、“相對偏應(yīng)力水平”概念，歸一化地考慮了圍壓、靜偏應(yīng)力以及動應(yīng)力等因素的影響，建立了累積塑性應(yīng)變與累積孔壓經(jīng)驗?zāi)Ｐ?。這類方法的好處是無需建立復(fù)雜的數(shù)學(xué)模型，模型參數(shù)較少且計算過程簡單，工程實用性較強。

不同于地震荷載，交通荷載引起的地層應(yīng)力幅值較小，但振動持續(xù)時間長，累計振動次數(shù)非常高，其引起的微振動不像地震那樣強烈，一般不會造成地面建筑物的直接坍塌和破壞。但是，這種振動是長期存在和反復(fù)發(fā)生的，會導(dǎo)致地基沉降或不均勻下沉，影響地面建筑物的安全和使用壽命。因此，對于長期循環(huán)荷載作用下土體累積變形的計算，關(guān)鍵不是計算每一次循環(huán)荷載作用過程中的變形情況，而是計算長期累積變形量。為此，本文基于前期開展的泥炭質(zhì)土不排水循環(huán)三軸壓縮試驗結(jié)果[10]，采用分離型非線性蠕變有限元法對泥炭質(zhì)土循環(huán)累積殘余變形簡化計算方法進(jìn)行了研究。根據(jù)長期循環(huán)荷載作用下泥炭質(zhì)土累積殘余變形發(fā)展規(guī)律與靜荷載下的蠕變趨勢相類似的特點，將土體的循環(huán)累積殘余變形等效為蠕變，基于靜力學(xué)蠕變有限元求解方法，建立土體循環(huán)累積殘余變形的有限元計算模型。最后，通過對泥炭質(zhì)土室內(nèi)不排水循環(huán)三軸壓縮試驗結(jié)果以及地鐵荷載下泥炭質(zhì)土地基的長期變形的數(shù)值模擬，驗證了該方法的有效性。

1 循環(huán)累積變形簡化計算方法

1.1 基本假定

目前，靜力學(xué)蠕變模型(如過應(yīng)力模型、非穩(wěn)態(tài)流動面模型等)中通常將總應(yīng)變拆分為可恢復(fù)的彈性應(yīng)變和不可恢復(fù)的黏塑性應(yīng)變。其中，將率相關(guān)的蠕變變形與率無關(guān)的塑性變形耦合在一起視為不可恢復(fù)的應(yīng)變，雖然一定程度上簡化了模型在算法實現(xiàn)上的難度，但物理本質(zhì)上有待商榷。如Dafalias[11]通過試驗研究認(rèn)為將土體不可恢復(fù)應(yīng)變拆分為率相關(guān)部分和率無關(guān)部分更加合理。鑒于此，本研究中總應(yīng)變采用分離型拆分方式，即

(1)

(2)

1.2 彈塑性變形

本研究中，采用臨界狀態(tài)理論計算與時間無關(guān)的彈塑性變形。對于彈性變形，其剛度常數(shù)為

(3)

對于塑性變形，采用遵循修正劍橋模型的屈服準(zhǔn)則以及相關(guān)聯(lián)流動法則，塑性勢函數(shù)與屈服面函數(shù)相同，即

(4)

1.3 循環(huán)累積變形

根據(jù)長期循環(huán)荷載作用下土體累積殘余變形發(fā)展規(guī)律與靜荷載下土體蠕變趨勢相類似的特點[12-14]，將土體的循環(huán)累積殘余變形等效為靜荷載下的蠕變。在這一等效方法基礎(chǔ)上，通過引入蠕變勢函數(shù)來描述土體的三維循環(huán)累積變形性質(zhì)，即

(5)

式中：φ為蠕變乘子；f為蠕變勢函數(shù)。可以看出，在蠕變乘子和蠕變勢函數(shù)已知的條件下就可以確定土體的循環(huán)累積應(yīng)變速率張量，進(jìn)而可以確定每一增量步的循環(huán)累積變形。

對于蠕變勢函數(shù)的確定，可取的方法是采用或改進(jìn)經(jīng)典彈塑性模型的勢函數(shù)。目前在黏彈塑性本構(gòu)模型研究中，普遍將蠕變勢函數(shù)設(shè)置為與塑性勢函數(shù)相同的形式[15-16]。本文同樣采用這一方法，蠕變勢函數(shù)取為

(6)

(7)

圖1 不排水蠕變試驗應(yīng)力路徑示意圖Fig.1 Stress path in undrained triaxial creep test

上述蠕變勢函數(shù)確定后，各循環(huán)累積應(yīng)變速率的方向也就隨之確定了。若再能確定應(yīng)變速率的大小，即蠕變乘子，則循環(huán)累積應(yīng)變速率分量的求解將得以解決。一般而言，蠕變乘子的確定方法有多種，如王者超等[17-18]借用等效應(yīng)力和等效應(yīng)變概念，建立了等效蠕變速率與等效蠕變應(yīng)力間的關(guān)系，此時蠕變乘子即等于等效蠕變速率的大小。Borja等[19]分別提出了兩種方法確定蠕變乘子的大小，一種是基于體積蠕變公式計算體積蠕變率，而后通過式(7)確定出蠕變乘子大小，進(jìn)而反求出偏蠕變速率。另一種方法恰好相反，先基于不排水偏蠕變公式計算偏蠕變速率，在利用式(7)確定出蠕變乘子大小，進(jìn)而反求出體積蠕變速率。李興照等[20]提出的模型中采用了第一種方法，并取得了較好的模擬結(jié)果。結(jié)合前期開展的泥炭土不排水三軸循環(huán)加載試驗，同時考慮到長期循環(huán)荷載作用下土中存在一定的累積超孔隙水壓力導(dǎo)致很難建立準(zhǔn)確描述土體累積塑性體積應(yīng)變的經(jīng)驗公式，本文采用上述Borja等提出的第二種方法，具體方法如下。首先，利用式(5)將偏應(yīng)變速率張量表述為

(8)

(9)

(10)

此時，有

(11)

由式(9)和式(11)可求出蠕變乘子φ，即

(12)

結(jié)合式(5)和式(12)，則土體的三維循環(huán)累積應(yīng)變速率張量可表述為

(13)

1.4 硬化規(guī)律

(14)

式中：θ=(1+e)/(λ-k)，λ為e-lnp空間中壓縮線斜率。

對于不排水循環(huán)三軸試驗，可將動偏應(yīng)力qd等價為靜力不排水三軸蠕變試驗中的偏應(yīng)力q，其應(yīng)力路徑可采用圖1 中OAB表示，其中OA段表示第一次循環(huán)中的加載過程，并將這一過程視為率無關(guān)的彈塑性變形過程，而AB段表示率相關(guān)的循環(huán)加載過程。土體狀態(tài)由O點變化到A點時，土體發(fā)生塑性硬化，前期固結(jié)壓力由pcO發(fā)展到pcA；隨后，當(dāng)土體狀態(tài)由A點變化到B點時，這一過程中平均有效應(yīng)力逐漸減小，應(yīng)力狀態(tài)點均位于屈服面內(nèi)部，土體只發(fā)生時間硬化，前期固結(jié)壓力由pcA發(fā)展到pcB。

2 數(shù)值實現(xiàn)

與經(jīng)典彈塑性問題不同，率相關(guān)問題的數(shù)值分析需對時間項進(jìn)行離散化， 即在離散的時刻上進(jìn)行各種力學(xué)參量的計算。

2.1 應(yīng)力積分算法

率相關(guān)問題中，土體的應(yīng)力應(yīng)變關(guān)系可表述為

(15)

(16)

(17)

對于這一典型的彈塑性問題，本文采用隱式圖形返回算法進(jìn)行積分點的本構(gòu)積分[21]，并給出一致性切線剛度張量。主要思路為：

(18)

步驟2利用試探應(yīng)力速率張量，并結(jié)合前一增量步結(jié)束時的應(yīng)力應(yīng)變狀態(tài)，檢查土體是否發(fā)生屈服。若F≤0，土體沒有進(jìn)入屈服狀態(tài)，試探應(yīng)力即為最終應(yīng)力，一致性切線剛度張量就是彈性剛度張量。返回步驟1，進(jìn)行下一步長的計算；若F>0，表明土體發(fā)生屈服，進(jìn)入下一步驟；

步驟3根據(jù)隱式圖形返回算法，通過推導(dǎo)可得

(19)

(20)

式中：Δt為增量步的時間增量；dφ為增量塑性乘子，可通過塑性一致性條件求取，即將式(19)～(21)代入式(4)，有

(22)

隨后，利用 Newton-Raphson 迭代方法求出dφ大小。據(jù)此，可更新應(yīng)力張量和一致性切線剛度張量。返回步驟1，進(jìn)行下一步長的計算。

2.2 程序?qū)崿F(xiàn)

借助大型商業(yè)軟件Abaqus提供的用戶子程序接口，將上述應(yīng)力積分算法編寫成UMAT子程序。UMAT子程序主要任務(wù)是根據(jù)Abaqus主程序?qū)氲膽?yīng)變增量確定出應(yīng)力增量，并更新相應(yīng)的狀態(tài)變量，同時給出材料的雅克比矩陣。基本流程為：首先，計算出時間增量步Δt內(nèi)的循環(huán)累積應(yīng)變，并將循環(huán)累積應(yīng)變從總應(yīng)變增量中去除，將率相關(guān)問題轉(zhuǎn)化為率無關(guān)問題；隨后，利用隱式應(yīng)力積分算法確定積分點應(yīng)力狀態(tài)，并更新彈、塑性應(yīng)變等狀態(tài)變量；最后，將更新后的應(yīng)力狀態(tài)下的一致性切線剛度矩陣賦值給雅克比矩陣數(shù)組。

3 方法驗證

為了驗證上述基于靜力學(xué)蠕變理論的土體循環(huán)累積變形簡化計算方法的適用性與穩(wěn)定性，本文對昆明泥炭質(zhì)土不排水循環(huán)三軸試驗結(jié)果以及地鐵循環(huán)荷載作用下泥炭質(zhì)土的長期變形進(jìn)行數(shù)值模擬。

3.1 泥炭質(zhì)土累積塑性應(yīng)變速率

泥炭(質(zhì))土是由古生物殘體在缺少空氣的條件下經(jīng)過復(fù)雜生物化學(xué)過程以及數(shù)萬年的地質(zhì)作用形成的特殊類軟土，是昆明市區(qū)廣泛分布的軟土地層，具有高有機質(zhì)含量、大孔隙比、高天然含水量、低承載力、高壓縮性等較差的工程地質(zhì)特性。作者曾對昆明泥炭質(zhì)土開展了系統(tǒng)的循環(huán)三軸壓縮試驗，其中不同圍壓下循環(huán)荷載引起的累積塑性應(yīng)變與循環(huán)周次曲線如圖2所示。圖中ηd為循環(huán)應(yīng)力比，定義為循環(huán)偏應(yīng)力幅值與土樣固結(jié)壓力之比。基于雙對數(shù)坐標(biāo)下累積應(yīng)變速率與循環(huán)周次間的線性關(guān)系，陳成等研究中建立了試樣軸向累積塑性應(yīng)變與動荷載循環(huán)周次間的經(jīng)驗關(guān)系，即

(23)

式中:a，b,m和n1為模型參數(shù)；N為循環(huán)周次；qd為循環(huán)偏應(yīng)力幅值；qf為靜破壞偏應(yīng)力；qs為循環(huán)偏應(yīng)力基準(zhǔn)值，用于描述非雙向正弦模式的單向循環(huán)加載。它與天然地基土中非等向固結(jié)引起的靜偏應(yīng)力不同，前者屬于循環(huán)荷載的一部分。為了描述靜偏應(yīng)力的影響效應(yīng)，借鑒Chai等的模型，將式(23)進(jìn)一步擴展為

(24)

式中：qs0為靜偏應(yīng)力，反映了土體的固結(jié)偏應(yīng)力狀態(tài)；n2為模型參數(shù)，Chai等的模型中建議n2=1。

此外，對于靜破壞偏應(yīng)力qf，有

(25)

將式(24)中的循環(huán)荷載作用次數(shù)看作時間度量單位，求導(dǎo)后可得出軸向循環(huán)累積應(yīng)變速率，將其代入式(13)即可利用上述UMAT子程序進(jìn)行計算。實際計算中，循環(huán)荷載引起的累積變形計算分為三步。第一步是地基土初始應(yīng)力的確定，包括地基土的圍壓水平、固結(jié)偏應(yīng)力狀態(tài)以及靜抗剪強度等；第二步是只施加一次循環(huán)荷載，用于獲取地基土中循環(huán)應(yīng)力幅值的分布情況。在前兩步的基礎(chǔ)上，第三步進(jìn)行地基土累積變形的計算。

3.2 單元試驗?zāi)M

有限元分析中，采用 1/4 土樣尺寸(直徑39.1 mm、高度80 mm)的幾何模型，單元類型為二階20節(jié)點三維實體單元(C3D20)。按照荷載控制方式進(jìn)行加載。首先，建立初始的自平衡應(yīng)力狀態(tài)；然后在試樣頂部按靜力方式施加與動偏應(yīng)力相等的荷載，最后進(jìn)行循環(huán)累積變形計算。

模型計算參數(shù)包括率無關(guān)的修正劍橋模型參數(shù)以及循環(huán)累積應(yīng)變經(jīng)驗?zāi)Ｐ椭械臄M合參數(shù)。針對昆明泥炭質(zhì)土，模型計算參數(shù)如表1所示(修正劍橋模型參數(shù)λ和κ通過一維壓縮試驗，參數(shù)M通過不排水三軸壓縮試驗確定)。

表1 模型參數(shù)匯總 Tab.1 Summary of model parameters

數(shù)值計算結(jié)果與實測結(jié)果對比同樣見圖2?？梢钥闯觯煌瑖鷫汉蛣悠珣?yīng)力水平組合試驗條件下，有限元計算結(jié)果與實測值相對吻合，表明本文提出的基于靜力學(xué)蠕變理論的土體循環(huán)累積變形簡化計算方法可以很好地描述循環(huán)偏應(yīng)力水平以及圍壓水平對泥炭質(zhì)土試樣循環(huán)累積變形特性的影響規(guī)律。

圖2 循環(huán)累積應(yīng)變與循環(huán)次數(shù)關(guān)系Fig.2 Accumulated strain versus number of cycles

3.3 地鐵荷載下泥炭質(zhì)土的長期變形模擬

采用平面應(yīng)變模型模擬地鐵循環(huán)荷載作用下泥炭質(zhì)土地基的長期變形。考慮隧道結(jié)構(gòu)對稱性，取一半隧道以及土體進(jìn)行計算。計算幾何模型和網(wǎng)格劃分如圖3所示，包括地基土，隧道襯砌以及泥漿等代層。其中，隧道埋深8.8 m，隧洞直徑6.2 m，襯砌和等代層的厚度分別為0.35 m和0.1 m。模型水平方向長度為5倍的洞徑，即31 m。豎向為30 m，包含三層地基土，分別是7 m厚黏土層，13 m厚的泥炭質(zhì)土層以及10 m厚的粉砂下臥層。有限元模型共包括1 655個4節(jié)點完全積分平面應(yīng)變單元(CPE4)。土體的邊界條件為左右兩邊水平方向位移約束，底部為豎向位移約束。同時在泥炭質(zhì)土層與粉砂層間設(shè)置透水邊界條件。計算過程中，只考慮泥炭質(zhì)土層的循環(huán)累積變形特性，計算參數(shù)如表1所示，其中e0取為4.68。黏土層采用修正劍橋模型，計算參數(shù)分別為：λ=0.173，κ=0.034，μ=0.33，e0=1.40以及M=1.22。粉砂層采用摩爾庫倫模型，計算參數(shù)分別為：壓縮模量E=13.5 MPa，μ=0.29，c=5 kPa以及φ=28°。

整個計算過程共分為4個步驟，即：①地應(yīng)力平衡，確定地基土初始圍壓、靜偏應(yīng)力水平以及土體的靜力強度等；②隧道開挖與支護；③施加一次地鐵運行引起的動荷載，獲取地基土中循環(huán)應(yīng)力幅值；④地鐵荷載下地基土累積變形計算。對于地鐵運行引起的動荷載，按行駛速度80 km/h計取為69 kPa，作用于寬度為1.44 m的軌道上[25]。

圖3 有限元網(wǎng)格劃分Fig.3 Finite element mesh

圖4 給出了泥炭質(zhì)土層中循環(huán)偏應(yīng)力幅值的分布云圖?？梢钥闯觯淼赖撞康鼗了艿降难h(huán)應(yīng)力幅值最大。隨與隧道距離的增加，循環(huán)應(yīng)力逐漸減小。

隧道上方以及遠(yuǎn)場地基土所受到的循環(huán)應(yīng)力均趨于0。圖5給出了地鐵運行2 000 000引起的地基土豎向位移分布(不包括隧道開挖引起的豎向位移)。圖6給出了隧道底部以及隧道正上方地表的沉降歷時曲線。其中，地鐵運行按年平均20 萬次計?？梢钥闯觯淼赖撞康某两盗看笥谒淼郎戏降乇淼某两盗?。隨著運行時間的增長，兩者的沉降量均逐漸增加，且增長速率逐漸減小。在運行2年、4年、6年、8年、10年情況下，隧道正上方地表的沉降量分別為1.14 cm，1.73 cm，2.19 cm，2.58 cm,2.92 cm。此外，列車運營荷載引起隧道結(jié)構(gòu)發(fā)生下沉，進(jìn)而使隧道上部土體產(chǎn)生變形，包括豎向沉降和朝隧道結(jié)構(gòu)方向的水平向位移，淺層土體朝隧道結(jié)構(gòu)方向的水平向位移一定程度上阻礙了其沉降的趨勢，因而導(dǎo)致隧道底部沉降略大于隧道正上方地表沉降。高靜連[26]也得出了類似的計算成果。

圖4 循環(huán)應(yīng)力分布Fig.4 Distribution of cyclic stress

圖5 地鐵運營引起的沉降分布Fig.5 Distribution of settlements induced by subway operation

圖6 沉降時間曲線Fig.6 Settlements versus time

4 結(jié) 論

基于長期循環(huán)荷載作用下土體累積殘余變形發(fā)展規(guī)律與靜荷載下土體蠕變變形規(guī)律相似的特點，提出了長期循環(huán)荷載作用下土體累積變形簡化計算方法，得出如下結(jié)論：

(1) 基于臨界狀態(tài)理論，將土體循環(huán)累積殘余變形等效為靜荷載下的蠕變，同時將循環(huán)荷載作用次數(shù)看作時間度量單位，通過引入與修正建橋模型屈服面一致的蠕變勢，結(jié)合不排水條件下土體偏循環(huán)累積應(yīng)變經(jīng)驗?zāi)Ｐ?，建立了可描述土體三維循環(huán)累積變形的簡化方法。

(2) 分項考慮率相關(guān)循環(huán)累積應(yīng)變與率無關(guān)彈塑性應(yīng)變，并利用隱式圖形返回算法進(jìn)行本構(gòu)積分。在此基礎(chǔ)上，基于Abaqus二次開發(fā)接口，編制了UMAT子程序。

(3) 采用提出的循環(huán)累積變形簡化計算方法對昆明泥炭質(zhì)土不排水循環(huán)三軸試驗結(jié)果以及地鐵荷載下泥炭質(zhì)土地基的長期變形特性進(jìn)行了數(shù)值模擬，驗證了該方法的適用性與穩(wěn)定性。

循環(huán)荷載作用下地基土的變形增長與所受臨界動應(yīng)力水平、循環(huán)速率頻率等因素密切相關(guān)，本文后續(xù)研究將進(jìn)一步探討這些因素對昆明泥炭質(zhì)土的循環(huán)累積變形的影響，并將其量化擴展至式(24)中，進(jìn)而使本文提出的計算方法在實際應(yīng)用中可反映這些因素的影響。