基于多場(chǎng)耦合數(shù)值分析的膨脹土邊坡淺層膨脹變形破壞機(jī)制研究

丁金華，陳仁朋，童 軍，龔壁衛(wèi)

（1. 長(zhǎng)江科學(xué)院 水利部巖土力學(xué)與工程重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430010；2. 浙江大學(xué) 巖土工程研究所，浙江 杭州 310058）

1 引 言

膨脹土具有脹縮性、裂隙性、超固結(jié)性等多種不良工程地質(zhì)特性，導(dǎo)致其分布地區(qū)頻發(fā)各種工程問(wèn)題或地質(zhì)災(zāi)害。大量膨脹土地區(qū)邊坡破壞實(shí)例的勘察結(jié)果表明，膨脹土邊坡滑坡具有明顯的漸進(jìn)性、淺層性、季節(jié)性和反復(fù)性等特征。已開建的南水北調(diào)中線工程總干渠沿線經(jīng)過(guò)的膨脹土地區(qū)累計(jì)長(zhǎng)度達(dá)340 km，渠坡在開挖過(guò)程中發(fā)生在邊坡淺表層一定深度范圍內(nèi)的牽引式滑坡破壞現(xiàn)象非常普遍，給工程處理帶來(lái)了極大的困難[1－2]。

國(guó)內(nèi)外很多學(xué)者從干濕循環(huán)、降雨、裂隙、強(qiáng)度等多方面開展了有關(guān)膨脹土邊坡破壞機(jī)制的研究工作，然而很多工程實(shí)例分析表明，即便采用膨脹土殘余強(qiáng)度進(jìn)行極限平衡計(jì)算，其安全系數(shù)也大于1。鑒于極限平衡法僅通過(guò)土的抗剪強(qiáng)度單一指標(biāo)反映邊坡的力系平衡，而不考慮變形問(wèn)題，越來(lái)越多的研究人員認(rèn)識(shí)到，對(duì)于膨脹土邊坡這種淺層漸進(jìn)性破壞，極限平衡分析是無(wú)法反映其失穩(wěn)狀態(tài)和破壞機(jī)制的[2]。有學(xué)者在極限平衡分析的基礎(chǔ)上陸續(xù)提出了各種改進(jìn)方法，典型的如Fredlund基于非飽和土理論提出的考慮基質(zhì)吸力影響的改進(jìn)極限平衡法，在安全系數(shù)計(jì)算公式中反映正/負(fù)孔隙水壓力的作用[3]。袁俊平[4]從非飽和膨脹土的裂隙概化模型入手，分析了裂隙對(duì)邊坡穩(wěn)定的影響。沈珠江[5]采用非飽和土簡(jiǎn)化固結(jié)理論模擬膨脹土渠道邊坡的降雨入滲，以反映邊坡在入滲過(guò)程中吸力喪失、有效應(yīng)力降低和土體膨脹回彈及水平變位的全過(guò)程。劉華強(qiáng)等[6]指出膨脹土邊坡穩(wěn)定性分析應(yīng)綜合考慮浸水后吸力消失、強(qiáng)度降低以及裂隙連通性等因素。

數(shù)值分析是另外一種研究膨脹土邊坡穩(wěn)定的有效手段[7－13]，但也存在若干問(wèn)題，如常規(guī)的飽和土本構(gòu)模型無(wú)法考慮膨脹土的吸濕膨脹效應(yīng)，只能忽略其膨脹性，代之以降低土體強(qiáng)度來(lái)近似反映，但膨脹土的吸濕膨脹變形不同于一般黏土在荷載作用下由顆粒剪切滑移形成的宏觀變形[14]，且對(duì)于非飽和膨脹土，不考慮孔隙中水/氣相互作用所得的模擬結(jié)果是偏于不安全的[15－16]。1999年，Alonso等[17]提出非飽和膨脹模型BExM，引入反映微觀結(jié)構(gòu)彈塑性變形的屈服面，可以考慮不同應(yīng)力路徑等復(fù)雜邊界條件以及干濕循環(huán)形成的累積塑性變形效應(yīng)，是目前公認(rèn)理論基礎(chǔ)嚴(yán)謹(jǐn)、體系較為完整的非飽和本構(gòu)模型。國(guó)內(nèi)也有學(xué)者以非飽和土理論和復(fù)合體損傷理論為基礎(chǔ)，引入非飽和彈塑性損傷本構(gòu)模型，對(duì)膨脹土邊坡進(jìn)行多場(chǎng)耦合數(shù)值模擬[18]，但受制于非飽和本構(gòu)模型的復(fù)雜性和參數(shù)獲取困難，很難在工程實(shí)踐中加以應(yīng)用和驗(yàn)證。1993年，繆協(xié)興[19]根據(jù)溫度應(yīng)力場(chǎng)的概念，針對(duì)膨脹巖首次提出了一種分析巖體受水作用而產(chǎn)生的應(yīng)力應(yīng)變場(chǎng)的新理論—“濕度應(yīng)力場(chǎng)”理論。該理論認(rèn)為，由于膨脹土吸水后的體積膨脹軟化現(xiàn)象與材料升溫后的體積膨脹軟化類似，且溫度場(chǎng)熱傳導(dǎo)方程和濕度場(chǎng)水分?jǐn)U散方程完全相似，因此，可以采用已較成熟的溫度應(yīng)力場(chǎng)理論來(lái)反映由于水分?jǐn)U散而形成的濕度場(chǎng)對(duì)應(yīng)力場(chǎng)的耦合影響。濕度應(yīng)力場(chǎng)理論以溫度場(chǎng)為媒介，可實(shí)現(xiàn)濕度和變形的耦合效應(yīng)，采用當(dāng)前大多數(shù)數(shù)值軟件即可進(jìn)行分析，但該理論基于溫度擴(kuò)散方程和水分?jǐn)U散方程的相似性，實(shí)質(zhì)上仍只是針對(duì)飽和滲流場(chǎng)的分析，不能反映膨脹土邊坡由降雨引起的土體非飽和-飽和入滲過(guò)程，也就無(wú)法模擬邊坡的漸進(jìn)性破壞。但是，該理論中有關(guān)膨脹土吸水產(chǎn)生膨脹類似于溫度變化產(chǎn)生膨脹，引入膨脹系數(shù)反映土體膨脹變形的思路是有參考價(jià)值的。

本文借鑒溫度場(chǎng)中溫度膨脹系數(shù)的定義和數(shù)值處理方法，結(jié)合膨脹土室內(nèi)無(wú)荷/有荷膨脹試驗(yàn)，對(duì)膨脹土濕度膨脹系數(shù)進(jìn)行反演擬合，獲得濕度場(chǎng)與溫度場(chǎng)、濕度場(chǎng)與膨脹系數(shù)的轉(zhuǎn)換關(guān)系；以降雨作用下膨脹土邊坡破壞的大型靜力物理模型試驗(yàn)成果為基礎(chǔ)，采用FLAC中的非飽和二相流模型模擬降雨引起的邊坡濕度場(chǎng)分布變化狀態(tài)；引入濕度場(chǎng)-溫度場(chǎng)及濕度場(chǎng)-膨脹系數(shù)的線性關(guān)系，提出一種適用于膨脹土邊坡的濕度場(chǎng)-膨脹變形場(chǎng)-應(yīng)力場(chǎng)的多場(chǎng)耦合數(shù)值分析方法，并藉此深入探討壓實(shí)膨脹土邊坡淺層漸進(jìn)性失穩(wěn)的全過(guò)程和破壞機(jī)制。

2 降雨作用下膨脹土邊坡淺層破壞的物理模型試驗(yàn)

2.1 模型試驗(yàn)破壞過(guò)程

邊坡物理模型試驗(yàn)方法及成果見(jiàn)文獻(xiàn)[20]，本文僅做簡(jiǎn)要敘述。

模型土樣取自河北邯鄲南水北調(diào)中線渠段，其黏土礦物成分主要為蒙脫石，含量為38%～45%，< 0.005 mm的黏粒含量達(dá)43.9%～53.5%，自由膨脹率大于 120%，屬?gòu)?qiáng)膨脹土。邊坡模型試驗(yàn)控制制備土樣的初始含水率為20%，干密度為1.60 g/cm3（壓實(shí)度為0.95）。模型邊坡高2.5 m（包括地基厚度為0.5 m），坡比1:1.5。試驗(yàn)中觀測(cè)位移、含水率、土體應(yīng)力以及環(huán)境溫/濕度和降雨量等。降雨強(qiáng)度平均為0.53 mm/h。

觀察試驗(yàn)全過(guò)程可見(jiàn)，受降水作用，試驗(yàn)持續(xù)27 h后邊坡下部即發(fā)生橫向裂縫，裂縫下方土體明顯向坡腳處滑移堆積，說(shuō)明此時(shí)邊坡在坡腳附近已局部失穩(wěn)。隨降水持續(xù)，裂縫進(jìn)一步向兩側(cè)擴(kuò)展，從模型箱側(cè)面觀察窗可發(fā)現(xiàn)邊坡中下部位內(nèi)部一定深度處已產(chǎn)生多處局部拉裂隙，同時(shí)坡頂也產(chǎn)生貫穿性的張拉裂縫。試驗(yàn)341 h時(shí)，邊坡上部近坡肩部位出現(xiàn)第3條貫穿性橫向裂縫，到385 h，坡體下部原第一條裂縫處土體發(fā)生快速塌落，繼而邊坡上部土體在2 min內(nèi)相繼滑塌，發(fā)生整體破壞。

試驗(yàn)結(jié)束后對(duì)邊坡進(jìn)行開挖，根據(jù)在軸線剖面預(yù)先埋置的標(biāo)志砂芯發(fā)現(xiàn)（見(jiàn)圖1），邊坡淺層多處發(fā)生了滑動(dòng)變形，甚至明顯的剪切錯(cuò)動(dòng)。圖1中同時(shí)繪出推測(cè)得到的滑裂面分布，可見(jiàn)邊坡從最初的淺表層局部滑動(dòng)面開始向深部發(fā)展，最終導(dǎo)致大范圍的滑坡產(chǎn)生，其中埋深約0.1 m和埋深約0.3 m處有明顯的剪切滑動(dòng)面，埋深0.3 m處的最大剪切位移可達(dá)20 cm左右。

圖1 模型試驗(yàn)滑坡后軸線剖面及滑動(dòng)面概化Fig.1 Axis section and sliding surfaces

2.2 試驗(yàn)監(jiān)測(cè)成果

試驗(yàn)過(guò)程中，在邊坡上部靠近坡肩處鉆孔取樣測(cè)量不同埋深處的含水率。由圖2可見(jiàn)，淺層0.2 m內(nèi)土體含水率在近100 h內(nèi)很快增大到30%以上，并在后續(xù)試驗(yàn)過(guò)程中基本穩(wěn)定在40%以上。隨試驗(yàn)進(jìn)程，水分逐漸向下遷移，在153～238 h時(shí)，表層0.3 m以內(nèi)土體含水率基本穩(wěn)定在40%以上，埋深0.3～0.4 m范圍內(nèi)的土體含水率有明顯增大。埋深0.5 m處的土體含水率在試驗(yàn)266～338 h內(nèi)變化顯著。發(fā)生整體滑坡后，測(cè)量結(jié)果顯示滑體上緣埋深0.6 m以下的土體水分變化不大，但在滑坡中下部，入滲影響深度基本接近 0.8 m?？梢?jiàn)，降雨引起邊坡的含水率場(chǎng)在時(shí)空分布上存在明顯的不均勻性。

圖2 模型試驗(yàn)實(shí)測(cè)含水率分布Fig.2 Distribution of actual measured water content of model test

圖3為由LVDT測(cè)得的邊坡面部分典型測(cè)點(diǎn)位移變化曲線。結(jié)合試驗(yàn)過(guò)程中的觀察現(xiàn)象，可知邊坡吸濕后的膨脹變形呈現(xiàn)以下特點(diǎn)：（1）邊坡土體在降水作用下吸濕發(fā)生膨脹變形的過(guò)程與土體含水率變化趨勢(shì)基本一致，但各部位膨脹變形量在時(shí)空上呈現(xiàn)明顯的不均勻分布。邊坡中下部土體膨脹變形速率明顯超過(guò)其他部位，最大變形量發(fā)生在坡體下部約1/3處（測(cè)點(diǎn)L5處）。與圖1對(duì)比可見(jiàn)膨脹變形越大處土體的剪切變形也越大；（2）坡肩和坡頂在試驗(yàn)初期隨著土體含水率的增大，表現(xiàn)為膨脹變形，但一段時(shí)間之后，隨著淺表層土體水分迅速增大到飽和狀態(tài)及坡肩張拉裂隙的發(fā)生，這些部位的表面變形表現(xiàn)為沉降，顯示出邊坡整體向下滑移的趨勢(shì)；（3）坡體中部的水平位移觀測(cè)點(diǎn)在滑坡后的位移量已達(dá)到約30 cm，坡頂水平位移約10 cm。

試驗(yàn)結(jié)束后在滑帶附近取樣進(jìn)行了室內(nèi)物理力學(xué)性質(zhì)試驗(yàn)，與試驗(yàn)前相比，滑帶附近土的含水率達(dá)40%以上，飽和度近90%，強(qiáng)度參數(shù)指標(biāo)均明顯降低，見(jiàn)表 1。采用殘余強(qiáng)度進(jìn)行模型邊坡的極限平衡穩(wěn)定計(jì)算，得到邊坡安全系數(shù)為1.90，無(wú)法反映滑坡的實(shí)際狀態(tài)。

表1 模型試驗(yàn)滑坡前后強(qiáng)膨脹土的物理力學(xué)性質(zhì)參數(shù)Table 1 Physico-mechanical parameters of high expansive soil before and after model test

3 多場(chǎng)耦合數(shù)值分析方法

3.1 多場(chǎng)耦合分析思路和計(jì)算流程

實(shí)現(xiàn)膨脹土邊坡多場(chǎng)耦合分析的關(guān)鍵是如何建立合理、準(zhǔn)確的濕度場(chǎng)和膨脹變形場(chǎng)的關(guān)系。本文借鑒“濕度應(yīng)力場(chǎng)”的理論假定，采用FLAC熱力學(xué)模型對(duì)膨脹土的室內(nèi)無(wú)荷/有荷膨脹試驗(yàn)進(jìn)行熱-固耦合預(yù)分析，獲得濕度場(chǎng)與溫度場(chǎng)、濕度場(chǎng)與膨脹系數(shù)的轉(zhuǎn)換關(guān)系，以此反映膨脹土由于水分增加而引起的膨脹變形效應(yīng)。

對(duì)于邊坡由于降雨引起的非飽和-飽和滲流場(chǎng)變化，采用FLAC的二相流模型進(jìn)行分析，獲得不同降雨時(shí)段對(duì)應(yīng)的濕度場(chǎng)分布狀態(tài)。然后，利用熱-固耦合預(yù)分析得到的濕度場(chǎng)與溫度場(chǎng)、濕度場(chǎng)與膨脹系數(shù)的轉(zhuǎn)化關(guān)系，進(jìn)行邊坡多場(chǎng)耦合分析。為方便起見(jiàn)，數(shù)值計(jì)算均采用節(jié)點(diǎn)飽和度代替含水率來(lái)反映邊坡不同部位的濕度場(chǎng)。

由于目前FLAC版本只能進(jìn)行飽和土單相流的流-固耦合計(jì)算或熱-固耦合計(jì)算，未提供非飽和二相流-熱-固模塊的耦合計(jì)算功能。因此，本文利用FISH語(yǔ)言編制相應(yīng)程序，最終實(shí)現(xiàn)膨脹土邊坡的濕度場(chǎng)-膨脹變形場(chǎng)-應(yīng)力場(chǎng)的耦合分析。數(shù)值分析流程可描述如下：

① 以膨脹土邊坡物理模型試驗(yàn)為研究對(duì)象，同時(shí)建立邊坡的熱-固力學(xué)概化模型和二相流概化模型。二者的剖分網(wǎng)格保持一致，分別賦予其相應(yīng)的模型參數(shù)、邊界條件及初始條件。

② 先計(jì)算力學(xué)模型中的自重應(yīng)力場(chǎng)，然后通過(guò)FISH程序?qū)⑵鋫鬟f給二相流模型作為初始應(yīng)力場(chǎng)，進(jìn)行第一時(shí)段非飽和滲流計(jì)算。

③ 根據(jù)熱-固耦合預(yù)分析得到的飽和度-溫度關(guān)系和飽和度-膨脹系數(shù)關(guān)系，編制FISH程序?qū)⒎秋柡蜐B流計(jì)算得到的濕度場(chǎng)分別轉(zhuǎn)化為溫度場(chǎng)和膨脹系數(shù)值，賦予邊坡力學(xué)模型作為下一時(shí)步計(jì)算的初始條件，同時(shí)將非飽和應(yīng)力場(chǎng)也傳遞給熱-固力學(xué)模型，求解該時(shí)段熱學(xué)和力學(xué)平衡。

④ 將力學(xué)模型計(jì)算得到的應(yīng)力場(chǎng)返回給二相流模型，作為下一時(shí)段非飽和計(jì)算的初始應(yīng)力場(chǎng)，重復(fù)上述過(guò)程直到計(jì)算時(shí)段結(jié)束。

強(qiáng)膨脹土的本構(gòu)關(guān)系采用摩爾-庫(kù)侖模型（M-C模型），其參數(shù)根據(jù)室內(nèi)試驗(yàn)確定。很多工程實(shí)踐和室內(nèi)試驗(yàn)都表明，膨脹土的“吸濕軟化”效應(yīng)非常明顯，表現(xiàn)為膨脹土遇水作用后，一方面其變形模量隨含水率增加而降低、泊松比增大；另一方面，由于含水率增大而使吸力降低，導(dǎo)致強(qiáng)度參數(shù)降低。本文在多場(chǎng)耦合計(jì)算中同時(shí)考慮了膨脹土的“應(yīng)變軟化”效應(yīng)，即當(dāng)土體飽和度達(dá)到 95%以上時(shí)，根據(jù)模型試驗(yàn)后滑帶土的室內(nèi)力學(xué)試驗(yàn)結(jié)果，重新調(diào)整強(qiáng)膨脹土的強(qiáng)度參數(shù)。

3.2 膨脹系數(shù)的熱-固耦合預(yù)分析

熱-固耦合預(yù)分析的目的在于確定濕度場(chǎng)和溫度場(chǎng)以及濕度場(chǎng)和膨脹系數(shù)的轉(zhuǎn)換關(guān)系，因此，以室內(nèi)大尺寸（φ = 152 mm，H = 50 mm）無(wú)荷/有荷膨脹試驗(yàn)為對(duì)象和各種含水率（飽和度）和干密度條件下的無(wú)荷/有荷膨脹試驗(yàn)得到的膨脹量為目標(biāo)，建立 FLAC各向同性熱傳導(dǎo)-力學(xué)耦合模型，計(jì)算反演不同條件下土的膨脹系數(shù)。

熱傳導(dǎo)系數(shù)取0.1 w/m·K，比熱為0.5 kJ/ kg·℃。假定初始狀態(tài)下含水率為 27.5%（相應(yīng)飽和度為76%）對(duì)應(yīng)初始溫度T0= 0 ℃，飽和度100%時(shí)對(duì)應(yīng)終了溫度T1= 100 ℃。

溫度邊界：室內(nèi)膨脹試驗(yàn)由試樣底部進(jìn)水，數(shù)值計(jì)算中將底部邊界設(shè)置為溫度邊界（100 ℃）。位移邊界：模型左右邊界取為X方向位移約束，底部為X-Y雙向位移約束。荷載邊界：與無(wú)荷膨脹試驗(yàn)對(duì)應(yīng)時(shí)上邊界取為自由邊界，模擬有荷膨脹試驗(yàn)時(shí)上邊界取為定荷載邊界。

反演擬合發(fā)現(xiàn)，無(wú)荷情況下的膨脹系數(shù)約在10-5量級(jí)，數(shù)值模擬精度較低，計(jì)算穩(wěn)定性差，最大誤差接近12%左右，而有荷情況下的膨脹系數(shù)均在10-4量級(jí)，數(shù)值模擬精度大多可控制在2%以內(nèi)。表2為各種條件下的膨脹系數(shù)擬合結(jié)果。

表2 熱-固耦合預(yù)分析得到的膨脹系數(shù)λTable 2 Swelling coefficient λ

由表2可見(jiàn)，膨脹系數(shù)λ受含水率、干密度、荷載的綜合影響，表現(xiàn)在：（1）膨脹系數(shù)與干密度呈正相關(guān)，但干密度變化對(duì)其影響較??；（2）僅在低荷載條件下（6.25 kPa），含水率與膨脹系數(shù)近似符合線性關(guān)系（見(jiàn)圖4），但荷載增大后各因素對(duì)膨脹系數(shù)的影響呈現(xiàn)復(fù)雜的相互交叉作用。

圖4 含水率-膨脹系數(shù)關(guān)系（ρd = 1.39 g/cm3）Fig.4 Relationship between water content and swelling coefficient

考慮到模型試驗(yàn)中邊坡淺層破壞主要發(fā)生在埋深0.2～0.5 m以內(nèi)（相當(dāng)于荷載約在10 kPa以內(nèi)），因此，耦合計(jì)算中采用本次預(yù)分析得到的荷載為6.25 kPa對(duì)應(yīng)的線性飽和度Sr-溫度T轉(zhuǎn)換關(guān)系及飽和度 Sr-膨脹系數(shù) λ關(guān)系來(lái)反映濕度變化引起的膨脹變形，即

4 多場(chǎng)耦合數(shù)值計(jì)算成果分析

4.1 邊坡濕度場(chǎng)分析

4.1.1 非飽和滲流模型及計(jì)算參數(shù)

FLAC二相流模型（two-phase flow model）假定介質(zhì)孔隙被兩種流體（氣/液）完全填充，且可用達(dá)西定律來(lái)描述每種流體。由于液相流動(dòng)在孔隙中占主要地位，導(dǎo)致氣相壓力高于液相壓力，二者之差即為毛細(xì)壓力，是飽和度的函數(shù)。

傳導(dǎo)定律（液相以上/下標(biāo)w表示，氣相以上/下標(biāo)g表示）可表示為

根據(jù)物理模型試驗(yàn)建立的二相流數(shù)值概化模型如圖5所示。圖中，同時(shí)表示出與模型試驗(yàn)典型測(cè)點(diǎn)位置相近的結(jié)點(diǎn)，3條垂線分別表示后續(xù)分析中所指的邊坡上、中、下部典型位置。為方便與試驗(yàn)結(jié)果對(duì)比，非飽和滲流數(shù)值計(jì)算的各個(gè)時(shí)段和模型試驗(yàn)降雨周期同步。

圖5 “多場(chǎng)”耦合數(shù)值分析的概化模型Fig.5 Generalized model of multifield coupling numerical analysis

強(qiáng)膨脹土的初始飽和度為76%，殘余飽和度近似取為 40%。根據(jù)邯鄲強(qiáng)膨脹土的土-水特征試驗(yàn)結(jié)果擬合得到的VG模型參數(shù)見(jiàn)表3。

表3 非飽和滲流計(jì)算參數(shù)Table 3 Parameters of two-phase flow model

邊坡表面設(shè)定為流量邊界。不考慮入滲強(qiáng)度與土體滲透性的關(guān)系，即忽略非飽和土入滲能力是與飽和度及吸力有關(guān)的變量，假定降雨全部入滲，以模型試驗(yàn)平均降雨強(qiáng)度0.53 mm/h作為流量邊界的入滲強(qiáng)度。

4.1.2 濕度場(chǎng)分布

根據(jù)二相流計(jì)算結(jié)果，提取對(duì)應(yīng)于模型試驗(yàn)鉆孔取樣測(cè)含水率位置處結(jié)點(diǎn)（埋深0.1 m和0.4 m處）的飽和度，繪出其隨降雨歷時(shí)的變化過(guò)程如圖6所示，圖中同時(shí)繪出實(shí)測(cè)含水率變化過(guò)程曲線。從圖中可以看出，非飽和滲流數(shù)值計(jì)算得到的土體濕度分布及變化過(guò)程和模型試驗(yàn)實(shí)測(cè)結(jié)果基本一致，表明本文采用的二相流數(shù)值模型和參數(shù)可適用于膨脹土邊坡的非飽和滲流計(jì)算，具有較高的準(zhǔn)確性。

圖6 數(shù)值計(jì)算飽和度與模型試驗(yàn)實(shí)測(cè)含水率的對(duì)比Fig.6 Comparison between measured water contents and calculated saturation degrees

圖7為不同降雨時(shí)段對(duì)應(yīng)的邊坡濕度場(chǎng)分布云圖，圖8為以邊坡中部測(cè)點(diǎn)為例示邊坡內(nèi)部非飽和浸潤(rùn)峰面隨降雨變化周期的推進(jìn)過(guò)程（圖中深色空心線條對(duì)應(yīng)降雨時(shí)刻，其中紅色實(shí)心線條對(duì)應(yīng)降雨間歇期）。從圖中可以看出，邊坡濕度場(chǎng)的變化與降雨密切相關(guān)，降雨期內(nèi)淺表層土體的暫態(tài)飽和區(qū)范圍明顯增大，干/濕土之間的濕潤(rùn)峰面非常明顯。降雨持續(xù)時(shí)間越長(zhǎng)，形成的暫態(tài)飽和區(qū)深度越大。在降雨間歇期，表層土的飽和度反而略有降低，受深層土體吸力的作用，水分繼續(xù)向下入滲，水力梯度減緩，濕潤(rùn)鋒被均化。

圖7 不同時(shí)段的濕度場(chǎng)分布狀態(tài)Fig.7 Distribution of moisture fields at different times

4.2 應(yīng)力、應(yīng)變場(chǎng)分析

4.2.1 邊界條件和計(jì)算參數(shù)

耦合計(jì)算力學(xué)模型網(wǎng)格剖分與二相流模型完全相同，如圖5所示。左右邊界設(shè)定為X方向位移約束，底部邊界為X-Y雙向約束。土體采用摩爾-庫(kù)侖模型（M-C模型），不考慮土的剪脹性。根據(jù)室內(nèi)物理力學(xué)試驗(yàn)結(jié)果，概化確定計(jì)算參數(shù)見(jiàn)表 4。對(duì)飽和度達(dá)到95%以上的土體，考慮應(yīng)變軟化效應(yīng)，根據(jù)模型試驗(yàn)滑坡后滑帶土的力學(xué)試驗(yàn)結(jié)果，降低其M-C模型參數(shù)。

圖8 非飽和浸潤(rùn)峰面隨降雨歷時(shí)的演進(jìn)過(guò)程（邊坡中部）Fig.8 Unsaturated infiltration peak at different times

表4 數(shù)值計(jì)算采用的強(qiáng)膨脹土力學(xué)強(qiáng)度參數(shù)Table 4 Mechanical strength parameters of expansive soil for numerical analysis

根據(jù)多場(chǎng)耦合分析流程，需要根據(jù)式（1）的飽和度-溫度關(guān)系將邊坡初始含水率場(chǎng)轉(zhuǎn)化為初始溫度場(chǎng)，作為耦合力學(xué)模型計(jì)算的初始條件。

4.2.2 膨脹變形

為便于與模型試驗(yàn)LVDT實(shí)測(cè)得到的典型測(cè)點(diǎn)坡面法向膨脹量進(jìn)行對(duì)比，將數(shù)值計(jì)算得到的相應(yīng)節(jié)點(diǎn)水平/垂直位移根據(jù)坡比轉(zhuǎn)換為法向變形量，繪出其變化時(shí)程曲線如圖9所示。由此可見(jiàn)，耦合數(shù)值分析得到的邊坡膨脹變形量及其發(fā)展過(guò)程與模型試驗(yàn)測(cè)量結(jié)果基本一致。在試驗(yàn)初期（0～153 h），計(jì)算位移和實(shí)測(cè)位移的誤差小于8%，到降雨后期，誤差有所增大。分析原因是，（1）在于耦合計(jì)算的位移精度與非飽和滲流計(jì)算時(shí)步大小密切相關(guān)，誤差存在累積效應(yīng)；（2）模型試驗(yàn)中邊坡地基側(cè)邊界并非完全剛性約束，而數(shù)值分析采用完全位移約束；（3）將數(shù)值計(jì)算得到的位移分量轉(zhuǎn)化為法向總變形時(shí)，沒(méi)有考慮邊坡坡比隨變形發(fā)展而實(shí)際發(fā)生變化所帶來(lái)的計(jì)算誤差，這種誤差在坡腳處更為明顯。

圖9 典型部位膨脹變形數(shù)值分析與模型試驗(yàn)實(shí)測(cè)值對(duì)比Fig.9 Comparison of swelling deformations between measurements and simulations typical measurement points

圖10 膨脹土邊坡的剪應(yīng)力場(chǎng)τxy分布（單位：Pa）Fig.10 Shear stress field τxy of expansive soil slope(unit: Pa)

計(jì)算位移與實(shí)測(cè)位移規(guī)律的基本一致性，說(shuō)明本文提出的多場(chǎng)耦合數(shù)值分析方法可以準(zhǔn)確地模擬降雨作用下膨脹土邊坡性態(tài)的變化，可以用來(lái)進(jìn)一步分析并評(píng)價(jià)邊坡的失穩(wěn)破壞過(guò)程。

4.2.3 應(yīng)力場(chǎng)分布

計(jì)算表明，考慮濕度場(chǎng)及膨脹變形耦合后邊坡的應(yīng)力應(yīng)變場(chǎng)發(fā)生重分布。限于篇幅，圖10僅繪出典型時(shí)段的剪應(yīng)力分布云圖。從計(jì)算結(jié)果可以看出，（1）豎向自重應(yīng)力變化不大，整體分布狀態(tài)變化不明顯，僅在近坡面附近由于暫態(tài)飽和區(qū)的擴(kuò)展，而使局部范圍內(nèi)的豎向應(yīng)力略有所增大；（2）大主應(yīng)力方向發(fā)生偏轉(zhuǎn)，基本與坡面平行，且隨降雨時(shí)段發(fā)展，在非飽和-飽和交界區(qū)域形成大主應(yīng)力集中區(qū)；（3）水平應(yīng)力分布狀態(tài)明顯改變，淺表層飽和區(qū)的水平應(yīng)力迅速增大，并在邊坡頂部、坡腳附近形成沿非飽和-飽和交界處的應(yīng)力集中區(qū)，隨非飽和浸潤(rùn)峰的推進(jìn)，水平應(yīng)力集中區(qū)也進(jìn)一步擴(kuò)展至相互貫通，最后形成明顯的貫穿全坡的集中帶；（4）剪應(yīng)力集中區(qū)首先在邊坡中下部、坡腳附近區(qū)域形成，降雨期剪應(yīng)力不僅明顯增大（其最大值可增大至初始狀態(tài)的5倍以上），且在邊坡表面和埋深較深處形成2個(gè)剪應(yīng)力集中區(qū)，以邊坡中部典型部位為例，圖11為剪應(yīng)力隨埋深的分布。隨非飽和-飽和區(qū)的擴(kuò)展，剪應(yīng)力集中帶也向深部上部擴(kuò)大，但其峰值略有均化的趨勢(shì)。邊坡下部剪應(yīng)力明顯大于邊坡中部。降雨間歇期由于表層土飽和度降低，深部土體飽和度增大，使得剪應(yīng)力集中區(qū)向深部發(fā)展，但其最大值較之降雨期有所減小。

圖11 剪應(yīng)力沿深度的分布（邊坡中部）Fig.11 Distribution of shearing stress along depth(at middle of slope)

4.2.4 邊坡穩(wěn)定性

本文耦合計(jì)算時(shí)段只模擬到模型試驗(yàn)進(jìn)行341 h，對(duì)應(yīng)邊坡產(chǎn)生較大變形，尚處于臨界狀態(tài)時(shí)刻，因此對(duì)于邊坡穩(wěn)定性采用應(yīng)力水平來(lái)進(jìn)行初步判斷。圖12為若干典型時(shí)刻應(yīng)力水平小于1的區(qū)域分布。從圖中可以看出，降雨初期邊坡下部近坡腳附近首先出現(xiàn)應(yīng)力水平小于1的局部區(qū)域，然后逐漸向兩方延伸，到226 h以后，臨界區(qū)已逐漸擴(kuò)展到從坡腳到邊坡中部，坡面以下埋深約0.2 m內(nèi)范圍，這與模型試驗(yàn)觀測(cè)到的該時(shí)刻邊坡淺層出現(xiàn)大量拉裂隙的現(xiàn)象是非常吻合的。降雨間歇期內(nèi)邊坡表層的含水率略有所降低，應(yīng)力重新調(diào)整分布后，臨界區(qū)范圍有所減?。?41 h時(shí)），說(shuō)明邊坡的穩(wěn)定性尚處于動(dòng)態(tài)變化過(guò)程中。

圖12 耦合數(shù)值分析的邊坡應(yīng)力水平臨界區(qū)域分布Fig.12 Critical zone (stress level less than 1) of expansive soil slope based on multifield coupling numerical analysis

4.3 不同計(jì)算模式的對(duì)比

實(shí)際工況中非飽和膨脹土邊坡穩(wěn)定受多種因素交叉影響，是典型的多場(chǎng)耦合作用問(wèn)題。本文提出的濕度場(chǎng)-變形場(chǎng)-應(yīng)力場(chǎng)耦合分析方法重點(diǎn)考慮了由于濕度變化引起的膨脹變形和土的強(qiáng)度軟化效應(yīng)。通過(guò)前述的數(shù)值分析與模型試驗(yàn)實(shí)測(cè)結(jié)果對(duì)比，已論證了該方法對(duì)膨脹土邊坡穩(wěn)定分析是基本合理準(zhǔn)確的。

為了進(jìn)一步比較膨脹變形對(duì)膨脹土邊坡穩(wěn)定的重要性，本文進(jìn)行了另外兩種計(jì)算模式的敏感性分析：（1）不考慮土體應(yīng)變軟化，僅采用膨脹系數(shù)反映膨脹變形；（2）忽略膨脹變形，僅考慮土體應(yīng)變軟化效應(yīng)。由于膨脹變形依托熱力學(xué)模塊進(jìn)行耦合分析，因此該模式只采用二相流模型的流-固耦合方式進(jìn)行計(jì)算。

圖 13為不同計(jì)算模式下得到的邊坡典型測(cè)點(diǎn)位移。與實(shí)測(cè)結(jié)果相比可知，吸濕膨脹變形以及應(yīng)變軟化都對(duì)膨脹土邊坡穩(wěn)定分析的準(zhǔn)確性具有重要影響，但二者的影響程度和影響方式略有不同：（1）忽略膨脹土的吸濕膨脹變形，僅考慮土強(qiáng)度的降低，計(jì)算得到的邊坡變形很小，這與普通黏土邊坡在降雨作用下的力學(xué)響應(yīng)相似，完全無(wú)法模擬出膨脹土邊坡在降水作用下變形破壞的逐漸發(fā)展過(guò)程，最終導(dǎo)致大大高估邊坡穩(wěn)定性。因此，對(duì)于膨脹土邊坡這種淺層漸進(jìn)性的滑坡，必須考慮吸濕膨脹變形的影響。（2）在降雨初期水分對(duì)土體強(qiáng)度特性的影響還不明顯的情況下，不考慮土體強(qiáng)度的軟化效應(yīng)，僅采用濕度場(chǎng)耦合膨脹變形的計(jì)算方法可以基本反映膨脹土邊坡變形的發(fā)展趨勢(shì)，但隨著水的影響的進(jìn)一步加劇，土體軟化效應(yīng)的重要性逐漸明顯，不考慮土強(qiáng)度的軟化就很難如實(shí)反映邊坡的長(zhǎng)期變形和穩(wěn)定狀態(tài)。

圖13 不同計(jì)算模式得到的邊坡變形對(duì)比Fig.13 Comparison of swelling deformation of different analysis methods

5 膨脹變形控制下的膨脹土邊坡淺層失穩(wěn)破壞機(jī)制

通過(guò)膨脹土邊坡的大型靜力物理模型試驗(yàn)及濕度場(chǎng)-膨脹變形場(chǎng)-應(yīng)力場(chǎng)多場(chǎng)耦合數(shù)值分析，揭示和闡明了膨脹變形作用是影響膨脹土邊坡淺層穩(wěn)定的制約性因素，其漸進(jìn)性失穩(wěn)破壞機(jī)制可闡述如下：（1）膨脹土邊坡淺層失穩(wěn)破壞主要受膨脹變形控制，是膨脹土膨脹性和環(huán)境條件綜合作用的結(jié)果。即便是均質(zhì)的膨脹土邊坡（不考慮超固結(jié)性，裂隙，干濕循環(huán)等因素），降雨入滲就可能引起膨脹土邊坡的淺表層剪切破壞，且這種淺層失穩(wěn)形式具有明顯的牽引性和多重滑動(dòng)性。（2）水力邊界條件變化會(huì)引起邊坡含水率場(chǎng)發(fā)生不均勻分布，繼而引起土體膨脹變形，邊坡內(nèi)應(yīng)力應(yīng)變場(chǎng)重分布，在非飽和-飽和交界區(qū)域形成應(yīng)力集中區(qū)，產(chǎn)生局部剪切破壞，并逐漸向邊坡深部擴(kuò)展，最終形成多重剪切面，邊坡發(fā)生失穩(wěn)破壞。因此，吸水膨脹性是膨脹土邊坡發(fā)生淺層滑動(dòng)的最本質(zhì)原因，而降雨等因素是提供土體水分發(fā)生變化的外部條件。（3）膨脹土邊坡淺層一定范圍內(nèi)土體含水率場(chǎng)分布具有明顯的時(shí)空不均勻性。降雨導(dǎo)致邊坡淺層形成明顯的暫態(tài)飽和區(qū)，降雨持續(xù)時(shí)間越長(zhǎng)，暫態(tài)飽和區(qū)范圍越大。降雨間歇期，暫態(tài)飽和區(qū)土體的飽和度略有降低，孔隙水壓力消散，受深層土體吸力的作用，水分會(huì)逐漸向坡體深部和上部擴(kuò)展，但深部土體含水率的變化較為緩慢。（4）膨脹土邊坡的膨脹變形過(guò)程受土體含水率分布的控制，也具有明顯的時(shí)空不均勻性。不同部位膨脹變形的發(fā)生發(fā)展過(guò)程不同步，且膨脹量大小差異幅度也較大。含水率首先發(fā)生變化的邊坡中下部土體其膨脹變形量率先增大，導(dǎo)致非飽和-飽和交界處產(chǎn)生明顯的剪切變形，邊坡局部產(chǎn)生張拉裂隙，但一段時(shí)間內(nèi)仍可保持邊坡的整體穩(wěn)定。隨降雨持續(xù)，膨脹變形繼續(xù)發(fā)展，非飽和-飽和浸潤(rùn)區(qū)向邊坡深部擴(kuò)展，交界處的剪切變形進(jìn)一步增大。（5）邊坡濕度場(chǎng)變化導(dǎo)致應(yīng)力-應(yīng)變場(chǎng)發(fā)生重分布，表現(xiàn)為豎向自重應(yīng)力的整體分布狀態(tài)變化不大，但邊坡內(nèi)水平應(yīng)力及剪應(yīng)力會(huì)顯著增大，并在邊坡內(nèi)發(fā)生重分布，特別會(huì)導(dǎo)致邊坡坡腳附近首先形成沿飽和-非飽和交界處的應(yīng)力集中區(qū)，發(fā)生局部剪切滑動(dòng)，但一段時(shí)間內(nèi)仍可保持邊坡的整體穩(wěn)定。隨非飽和浸潤(rùn)峰的推進(jìn)，水平應(yīng)力集中區(qū)和剪應(yīng)力集中區(qū)也進(jìn)一步擴(kuò)展至相互貫通，最后形成明顯的貫穿全坡的集中帶，在坡頂/坡肩附近產(chǎn)生張拉裂隙，最終導(dǎo)致邊坡整體滑動(dòng)，形成多重剪切破壞面，具有明顯的牽引式滑動(dòng)特點(diǎn)。

6 結(jié) 論

（1）膨脹土邊坡淺層破壞具有明顯的牽引性和多重滑動(dòng)性，其失穩(wěn)主要受膨脹變形控制，是膨脹土膨脹性和環(huán)境條件綜合作用的結(jié)果，超固結(jié)性、裂隙、干濕循環(huán)等因素僅起到一定的促進(jìn)作用。

（2）對(duì)于膨脹土邊坡淺層漸進(jìn)性失穩(wěn)破壞，極限平衡分析是無(wú)法準(zhǔn)確反映其失穩(wěn)機(jī)制的。只有采用考慮濕度場(chǎng)-膨脹變形場(chǎng)-應(yīng)力場(chǎng)多場(chǎng)耦合作用的數(shù)值分析方法，才能合理準(zhǔn)確地模擬邊坡在水分入滲作用下的變形破壞過(guò)程。

（3）在傳統(tǒng)飽和土本構(gòu)模型基礎(chǔ)上，借鑒目前較為成熟的溫度應(yīng)力場(chǎng)數(shù)值處理方法，將非飽和滲流場(chǎng)和膨脹變形場(chǎng)及應(yīng)力場(chǎng)進(jìn)行了初步耦合計(jì)算，物理概念清晰、參數(shù)易于確定，是一種便于解決工程實(shí)踐需要的實(shí)用數(shù)值模型和計(jì)算方法。通過(guò)對(duì)模型試驗(yàn)的仿真模擬，證明這種多場(chǎng)耦合分析模型和方法的可信性和可行性，為深入理解膨脹土邊坡淺層失穩(wěn)破壞機(jī)制和穩(wěn)定性分析提供了理論依據(jù)和有效分析手段。

（4）多場(chǎng)耦合數(shù)值分析的關(guān)鍵在于確定濕度場(chǎng)和膨脹變形場(chǎng)之間的關(guān)系?；跍囟扰蛎浐蜐穸扰蛎浀南嗨菩?，借鑒溫度應(yīng)力場(chǎng)的數(shù)值處理方法，根據(jù)室內(nèi)無(wú)荷/有荷膨脹試驗(yàn)成果，反演得到不同含水率、干密度和荷載作用下的膨脹系數(shù)變化規(guī)律。在耦合數(shù)值計(jì)算中采用了低荷載條件下對(duì)應(yīng)的含水率-膨脹系數(shù)的線性關(guān)系來(lái)反映邊坡淺層濕度變化引起的膨脹系數(shù)的變化規(guī)律。

[1] 李青云, 程展林, 龔壁衛(wèi), 等. 南水北調(diào)中線膨脹土巖地段渠道破壞機(jī)制和處理技術(shù)研究[J]. 長(zhǎng)江科學(xué)院院報(bào), 2009, 26(11): 1－9.LI Qing-yun, CHENG Zhan-lin, GONG Bi-wei, et al.Failure mechanism and treatment technology of expansive soil slope of middle route project of South-North Water diversion[J]. Journal of Yangtze River Scientific Research Institute, 2009, 26(11): 1–9.

[2] 程展林, 李青云, 郭熙靈, 等. 膨脹土邊坡穩(wěn)定性研究[J]. 長(zhǎng)江科學(xué)院院報(bào), 2011, 28(10): 102－111.CHENG Zhan-lin, LI Qing-yun, GUO Xi-lin, et al. Study on the stability of expansive soil slope[J]. Journal of Yangtze River Scientific Research Institute, 2011,28(10): 102–111.

[3] FREDLUND D G, MORGENSTERN N R, WIDGER R A. The shear strength of unsaturated soils[J]. Canadian Geotechnical Journal, 1978, 15(3): 313－321.

[4] 袁俊平. 非飽和膨脹土的裂隙概化模型與邊坡穩(wěn)定研究[D]. 南京: 河海大學(xué), 2003.

[5] 沈珠江, 米占寬. 膨脹土渠道邊坡降雨入滲和變形耦合分析[J]. 水利水運(yùn)工程學(xué)報(bào), 2004, 26(3): 7－11.SHEN Zhu-jiang, MI Zhan-kuan. Coupled analyses of seepage and deformation of expansive soil slopes during rainfall[J]. Hydro-Science and Engineering, 2004, 26(3):7－11.

[6] 劉華強(qiáng), 殷宗澤. 膨脹土邊坡穩(wěn)定分析方法研究[J]. 巖土力學(xué), 2010, 31(5): 1545－1549.LIU Hua-qiang, YIN Zong-ze. Research on analytical method of stability of expansive soil slope[J]. Rock and Soil Mechanics, 2010, 31(5): 1545－1549.

[7] NG C W W, SHI Q. A numerical investigation of the stability of unsaturated soil slopes subjected to transient seepage[J]. Computers and Geotechnics, 1998, 22(1): 1－28.

[8] CHO S E, LEE S R. Instability of unsaturated soil slopes due to infiltration[J]. Computers and Geotechnics, 2001,28(3): 185－208.

[9] 陳善雄, 陳守義. 考慮降雨的非飽和土邊坡穩(wěn)定性分析方法[J]. 巖土力學(xué), 2001, 22(4): 447－450.CHEN Shan-xiong, CHEN Shou-yi. Analysis of stability of unsaturated soil slope due to permeation of rainwater[J]. Rock and Soil Mechanics, 2001, 22(4): 447－450.

[10] 姚海林, 鄭少河, 李文斌, 等. 降雨入滲對(duì)非飽和膨脹土邊坡穩(wěn)定性影響的參數(shù)研究[J]. 巖石力學(xué)與工程學(xué)報(bào), 2002, 21(7): 1034－1039.YAO Hai-lin, ZHENG Shao-he, LI Wen-bin, et al.Parametric study on the effect of rain infiltration on stability of unsaturated expansive soil slope[J]. Chinese Journal of Rock Mechanics and Engineering, 2002,21(7): 1034－1039.

[11] 劉義高, 周玉峰, 鄭健龍, 等. 增濕條件下膨脹土路塹邊坡穩(wěn)定性數(shù)值分析[J]. 巖土工程學(xué)報(bào), 2007, 29(12):1870－1875.LIU Yi-gao，ZHOU Yu-feng，ZHENG Jian-long, et al.Numerical analysis of cutting slope stability in expansive soils considering water-increased state[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(12):1870－1875.

[12] DYE H B, HOUSTON S L, WELFERT B D. Influence of unsaturated soil properties uncertainty on moisture flow modeling[J]. Geotechnical and Geological Engineering,2011, 29(2): 161－169.

[13] HOSSAIN J. Geohazard potential of rainfall induced slope failure on expansive clay[D]. Texa: The University of Texas, 2012.

[14] 陳鐵林, 周成, 沈珠江. 黏土變形特性分析[J]. 巖土力學(xué), 2004, 25(4): 522－528.CHEN Tie-lin, ZHOU Cheng, SHEN Zhu-jiang. Analysis of deformation characteristics of clay[J]. Rock and Soil Mechanics, 2004, 25(4): 522－528.

[15] 王鐵行, 李寧, 謝定義, 等. 土體水熱力耦合問(wèn)題研究意義、現(xiàn)狀及建議[J]. 巖土力學(xué), 2005, 26(3): 488－493.WANG Tie-hang, LI Ning, XIE Ding-yi, et al. Necessity and means in research on soil coupled heatmoisture-stress issues[J]. Rock and Soil Mechanics,2005, 26(3): 488－493.

[16] 詹良通. 非飽和膨脹土邊坡中土水相互作用機(jī)制[J].浙江大學(xué)學(xué)報(bào)(工學(xué)版), 2006, 40(3): 494－500.ZHAN Liang-tong. Study on soil-water interaction in unsaturated expansive soil slope[J]. Journal of Zhejiang University (Engineering Science), 2006, 40(3): 494－500.

[17] ALONSO E E, VAUNAT J, GENS A. Modelling the mechanical behaviour of expansive clays[J]. Engineering Geology, 1999, 54(1): 173－183.

[18] 盧再華, 陳正漢, 方祥位, 等. 非飽和膨脹土的結(jié)構(gòu)損傷模型及其在土坡多場(chǎng)耦合分析中的應(yīng)用[J]. 應(yīng)用數(shù)學(xué)和力學(xué), 2006, 27(7): 781－788.LU Zai-hua, CHEN Zheng-han, FANG Xiang-wei, et al.Structural damage model of unsaturated expansive soil and its application in multi-field couple analysis on expansive soil slope[J]. Applied Mathematics and Mechanics, 2006, 27(7): 781－788.

[19] 繆協(xié)興, 楊成永, 陳至達(dá). 膨脹巖體中的濕度應(yīng)力場(chǎng)理論[J]. 巖土力學(xué), 1993, 14(4): 49－55.MIAO Xie-xing, YANG Cheng-yong, CHEN Zhi-da.Humidity stress field theory in swelling rock mass[J].Rock and Soil Mechanics, 1993, 14(4): 49－55.

[20] CHENG Zhan-lin, DING Jin-hua, RAO Xi-bao, et al.Physical model tests of expansive soil slope[J]. Chinese Journal of Geotechnical Engineering, 2014, 36(4): 716－723.