陳 成， 周正明， 張先偉

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

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

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

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

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

1.1 基本假定

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

(1)

(2)

1.2 彈塑性變形

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

(3)

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

(4)

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

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

(5)

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

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

(6)

(7)

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

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

(8)

(9)

(10)

此時，有

(11)

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

(12)

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

(13)

1.4 硬化規(guī)律

(14)

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

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

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

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

2.1 應力積分算法

率相關問題中，土體的應力應變關系可表述為

(15)

(16)

(17)

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

(18)

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

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

(19)

(20)

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

(22)

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

2.2 程序實現(xiàn)

借助大型商業(yè)軟件Abaqus提供的用戶子程序接口，將上述應力積分算法編寫成UMAT子程序。UMAT子程序主要任務是根據(jù)Abaqus主程序導入的應變增量確定出應力增量，并更新相應的狀態(tài)變量，同時給出材料的雅克比矩陣?；玖鞒虨椋菏紫龋嬎愠鰰r間增量步Δt內的循環(huán)累積應變，并將循環(huán)累積應變從總應變增量中去除，將率相關問題轉化為率無關問題；隨后，利用隱式應力積分算法確定積分點應力狀態(tài)，并更新彈、塑性應變等狀態(tài)變量；最后，將更新后的應力狀態(tài)下的一致性切線剛度矩陣賦值給雅克比矩陣數(shù)組。

3 方法驗證

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

3.1 泥炭質土累積塑性應變速率

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

(23)

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

(24)

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

此外，對于靜破壞偏應力qf，有

(25)

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

3.2 單元試驗模擬

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

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

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

數(shù)值計算結果與實測結果對比同樣見圖2。可以看出，不同圍壓和動偏應力水平組合試驗條件下，有限元計算結果與實測值相對吻合，表明本文提出的基于靜力學蠕變理論的土體循環(huán)累積變形簡化計算方法可以很好地描述循環(huán)偏應力水平以及圍壓水平對泥炭質土試樣循環(huán)累積變形特性的影響規(guī)律。

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

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

采用平面應變模型模擬地鐵循環(huá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厚的泥炭質土層以及10 m厚的粉砂下臥層。有限元模型共包括1 655個4節(jié)點完全積分平面應變單元(CPE4)。土體的邊界條件為左右兩邊水平方向位移約束，底部為豎向位移約束。同時在泥炭質土層與粉砂層間設置透水邊界條件。計算過程中，只考慮泥炭質土層的循環(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個步驟，即：①地應力平衡，確定地基土初始圍壓、靜偏應力水平以及土體的靜力強度等；②隧道開挖與支護；③施加一次地鐵運行引起的動荷載，獲取地基土中循環(huán)應力幅值；④地鐵荷載下地基土累積變形計算。對于地鐵運行引起的動荷載，按行駛速度80 km/h計取為69 kPa，作用于寬度為1.44 m的軌道上[25]。

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

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

隧道上方以及遠場地基土所受到的循環(huán)應力均趨于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。此外，列車運營荷載引起隧道結構發(fā)生下沉，進而使隧道上部土體產生變形，包括豎向沉降和朝隧道結構方向的水平向位移，淺層土體朝隧道結構方向的水平向位移一定程度上阻礙了其沉降的趨勢，因而導致隧道底部沉降略大于隧道正上方地表沉降。高靜連[26]也得出了類似的計算成果。

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

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

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

4 結 論

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

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

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

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

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