洪心怡,劉志強,胡春洋,熊勇林

(寧波大學(xué) 土木工程與地理環(huán)境學(xué)院,浙江 寧波 315211)

0 引言

降雨入滲是引發(fā)邊坡失穩(wěn)破壞的常見誘因[1-2]。而在自然界中,邊坡土體多處于非飽和狀態(tài)。因此,開展非飽和土坡在降雨條件下滲透變形破壞的研究,對指導(dǎo)邊坡防護工作具有重要意義。

目前,針對邊坡在降雨入滲下的破壞問題,主要通過實驗與數(shù)值計算這2類方法進行研究[3-10]。實驗大致可分為2類,即探究不同初始條件,如不同坡腳、初始含水量和降雨強度等與邊坡穩(wěn)定性之間的關(guān)系(如李煥強等[4]、王福恒等[5])和探究降雨入滲作用下邊坡性質(zhì)的變化和形成機制(如張磊等[3]、王維早等[6]),非飽和土坡受到雨水滲入后,其整個坡體內(nèi)部的變化很難通過實驗來予以詳細地展示。數(shù)值計算由于其具有適應(yīng)性強、精度高、可重復(fù)等優(yōu)點受到學(xué)者們的廣泛應(yīng)用。劉子振等[7]通過不同的數(shù)值計算方法求得臨界平衡狀態(tài)下滑體條塊的相互作用力系數(shù)和非飽和邊坡安全系數(shù),平揚等[8]對降雨入滲條件下的膨脹土邊坡穩(wěn)定性進行了分析。汪洋等[9]研究了小降雨入滲的非飽和土邊坡穩(wěn)定性計算方法,榮冠等[10]編寫了非飽和滲流程序分析降雨入滲機理和模擬方法,張社榮等[11]研究了多層巖質(zhì)邊坡在巖層傾角、邊坡坡角、結(jié)構(gòu)面間距不同時的破壞機制。石振明等[12]通過改進Green-Amp入滲模型和強度折減法來分析多層非飽和土邊坡的穩(wěn)定性計算。謝秀棟等[13]通過編制使用程序計算相應(yīng)的可靠指標(biāo),對邊坡工程的整體穩(wěn)定可靠度進行分析。陳勇等[14]根據(jù)非飽和土Barcelona模型進行了三維有限元數(shù)值模擬,得出非飽和土邊坡的應(yīng)力位移分布。宗振邦等[15]采用隨機場-貝葉斯方法校準(zhǔn)了關(guān)鍵巖土參數(shù)并對邊坡穩(wěn)定性的可靠度進行了分析。但是,以上研究對象皆為單一均質(zhì)土邊坡,對于多層坡體的數(shù)值分析則相對較少。受自然環(huán)境中降雨、巖石風(fēng)化、動植物活動等因素的影響,土體通常在同一斷面出現(xiàn)分層現(xiàn)象。不同土層之間,其力學(xué)參數(shù)以及飽和滲透系數(shù)均會發(fā)生變化。而且目前為止,大部分的研究成果都是基于極限平衡法或強度折減法來分析邊坡的穩(wěn)定性,沒有對其進行瞬時滲流-變形耦合分析。

因此,本研究基于多相混合體理論,推導(dǎo)得到土-水-氣三相耦合控制方程組[16-17],同時結(jié)合ZHANG等[18]提出的非飽和土彈塑性本構(gòu)模型,采用有限元與有限差分法進行空間時間離散化,最終得到土-水-氣三相耦合有限元控制方程,并利用Fortran90語言編制成計算程序。以KITAMURA等[19]開展的2層不同濕密度的非飽和砂土邊坡模型實驗為研究對象,使用計算程序?qū)Σ煌涤陱姸扰c坡角的4個工況進行數(shù)值模擬分析,驗證了文中所提方法數(shù)值分析方法的可靠性及實用性。

1 土-水-氣三耦合有限元控制方程

1.1 控制方程

土-水-氣三相耦合控制方程分為3個部分:第1部分是混合體平衡方程式;第2部分是土-水連續(xù)方程式;第3部分是土-氣連續(xù)方程式。3個方程式是基于土-水-氣各相的動量守恒定律和質(zhì)量守恒定律推導(dǎo)而得[16-17]:

1)平衡方程式:

(1)

2)土-水連續(xù)方程式:

(2)

3)土-氣連續(xù)方程式:

(3)

式中:σjj與εii分別為應(yīng)力張量和應(yīng)變張量;ρ為土體密度;bi為體積力張量;n為孔隙率;pw為水壓;pa為氣壓;Sr為飽和度;k為滲透系數(shù);K為體積模量;γ為重度; 上標(biāo)s、w、a分別為固液體和氣體。

1.2 控制方程空間/時間離散化

將上節(jié)的平衡方程式、土-水和土-氣連續(xù)方程式進行空間和時間離散化后,即可得到土-水-氣三相耦合有限元方程。平衡方程是基于有限元法進行空間離散,而連續(xù)方程是基于有限差分法而完成的。

首先將平衡方程進行空間離散化,基于虛位移原理,并進行一系列的計算,可得平衡方程增量的弱形式:

(4)

文中采用ZHANG等[18]提出的以Bishop有效應(yīng)力和飽和度作為狀態(tài)變量的非飽和土彈塑性本構(gòu)模型。該模型在飽和修正劍橋模型的基礎(chǔ)上,通過添加飽和度對土體孔隙比的影響規(guī)律而建立,模型可以自由連續(xù)地同時描述飽和土和非飽和土的力學(xué)性質(zhì)。其應(yīng)力增量有限元表達式如式(5):

Δσ′=DepBΔuN-ERFSΔSr

(5)

將式(5)代入到式(4)中,得到如式(6):

(6)

眾所周知,非飽和土的飽和度與其基質(zhì)吸力之間有一定的關(guān)系,文中所采用的土-水特征曲線為ZHANG等[18]提出的能夠同時考慮干、濕循環(huán)效應(yīng)的土-水特征曲線(soil water characteristic curve, SWCC)。基質(zhì)吸力增量與飽和度的關(guān)系如式(7):

(7)

(8)

(9)

將式(6)～式(9)進行時間離散化,并計算得到平衡方程時間和空間離散化后的方程式:

(10)

將土-水和土-氣連續(xù)方程式進行時間和空間離散化,得到如式(11)、式(12)所示:

(11)

(12)

(13)

2 邊坡模型的數(shù)值計算

2.1 計算模型及其初始條件的確定

KITAMURA等[19]為了探究降雨強度和雙層邊坡坡角對非飽和土雙層邊坡的破壞影響,設(shè)計了如下實驗方案:首先選取日本常見的Shirasu砂土,制作了2種不同尺寸的邊坡模型,模型是由2種不同濕密度的土組成。將右側(cè)的土體記為A層,左側(cè)的土體記為B層。2種邊坡模型的實際尺寸以及孔壓計布設(shè)具體位置如圖1所示,并對每種邊坡模型A層頂部分別進行降雨強度為50、100 mm/h的散水實驗。在實驗的過程中,一直保持降雨狀態(tài)直至邊坡發(fā)生破壞,其余實驗條件皆相同。表1列出了4種工況的實驗條件。

圖1 實驗?zāi)Ｐ偷某叽缂翱讐河嫷姆植嘉恢脠DFig. 1 Size of the experimental model and the distribution position of the pore pressure gauge

表1 各工況的實驗條件對比表Table 1 Comparison of experimental conditions of various working conditions mm/h

文中選此實驗進行數(shù)值分析。圖2反映了數(shù)值計算所使用的2個模型在二維平面內(nèi)有限元網(wǎng)格的劃分以及邊界條件。2種不同尺寸的計算模型都劃分成1 500個單元,模型的底部固定,右側(cè)橫向約束,其余各邊自由。模型的上部,左側(cè)以及下部為排氣排水面,右側(cè)為不排水不排氣面。表2列出了實驗所用A、B這2層砂土的物理特性參數(shù)。從圖3可知非飽和度狀態(tài)變量ρs的計算方法,即通過在參考應(yīng)力pr時非飽和土和飽和土的孔隙比之差求得。表3列出了計算時所用的非飽和土本構(gòu)模型的材料參數(shù)。計算模型的初始應(yīng)力場通過2方面確定,首先根據(jù)自重應(yīng)力進行計算,其次在邊坡模型的制作過程中,考慮到對砂土的擊實,故在自重應(yīng)力的基礎(chǔ)上額外增加4 kPa的應(yīng)力。

圖2 數(shù)值計算模型及邊界條件 圖3 非飽和土及飽和土的e-lnp曲線Fig. 2 Numerical calculation model and boundary conditions Fig. 3 The e-lnp curve of unsaturated soil and saturated soil

表2 砂土物理特性表Table 2 Physical properties of sand

表3 砂土的材料參數(shù)表Table 3 Material parameters of sand

圖4顯示了A、B這2層土樣土-水特征曲線的實驗數(shù)據(jù)與計算曲線。從圖中可知,計算曲線能較好地擬合實驗數(shù)據(jù)。計算曲線是基于ZHANG等[18]提出的SWCC模型而得,計算所用的土-水特性曲線模型參數(shù)列于表4。在邊坡模型進行摸擬降雨實驗前,先測其4個工況A、B土層土樣的初始負孔壓值。因為邊坡模型的尺寸較小,故可以將初始氣壓視為大氣壓力,其大小設(shè)定為0 kPa。為了方便計算,根據(jù)實驗測得的初始負孔壓值,將4個工況邊坡的每層土整體都取為同一個負孔壓值。并通過圖4的土-水特征曲線,確定計算的初始飽和度,見表5。

表4 A、B層砂土的土-水特征曲線參數(shù)表Table 4 Parameters of the soil-water characteristic curves of the A、B layer sand

圖4 A、B層土-水特征曲線Fig. 4 Soil-water characteristic curves of layer A、B

表5 計算初始負孔壓值和初始飽和度Table 5 Calculation of initial negative pore pressure and initial saturation

一般而言,非飽和土的透水系數(shù)與其飽和度相關(guān)。文中基于MUALEM[20]所提的公式來確定A、B這2土層土樣飽和度與透水系數(shù)的關(guān)系。從圖5可知,A、B層土樣的透水系數(shù)隨著飽和度的增大而增大。

圖5 飽和度與A、B土層透水系數(shù)的關(guān)系Fig. 5 Relationship between saturation and permeability coefficient of A、B soil layer

由于實驗所用的砂土透氣性較好,模型的尺寸不大以及模型的3個邊界都為透氣界面,因此在實驗過程中產(chǎn)生的空氣壓力很小,可忽略不計。在此文中不再過多探究氣壓對于邊坡雨水入滲破壞的影響。

2.2 數(shù)值計算結(jié)果與實驗對比分析

4種工況實驗和計算結(jié)果的對比見圖6。在初始階段,各點位的負孔隙水壓都趨于平穩(wěn)。隨著降雨的持續(xù)進行,各點位的負孔隙水壓依次在某段時間內(nèi)快速上升。實驗和計算曲線終止的時間點為邊坡破壞時刻。從圖6可知,4種工況的計算結(jié)果在總體上與實驗數(shù)據(jù)吻合,可以較為準(zhǔn)確地反映在不同降雨強度和坡角條件下,邊坡各點位負孔隙水壓隨時間的變化過程以及預(yù)測邊坡破壞的時間。但是,計算曲線最后達到穩(wěn)定時的負孔隙水壓值總是高于實驗結(jié)果。其產(chǎn)生的原因是,計算所用的土體透水系數(shù)相對偏小,使得雨水在滲入土體時滯留增多。

圖6 負孔隙水壓力的計算值與實驗值的對比Fig. 6 Comparison of calculated and experimental values of negative pore water pressure

2.3 邊坡水分遷移的計算結(jié)果分析

飽和度和吸力值具有一一對應(yīng)的關(guān)系,飽和度越高,則吸力值越小,見圖4。故可從吸力在各時間點的變化來反映邊坡的水分遷移。圖7為4種工況在50 min和坡體破壞時吸力的計算結(jié)果云圖。雨水是從A層土體的頂部流入,并向下滲透。在初始階段A層土體的透水系數(shù)高于B層,故在短時間內(nèi)(50 min)B層只有少量的土體受到了滲透。邊坡破壞時底部吸力值仍然大于上部,其原因在于邊坡底面為排氣排水面,水可從底面排出,邊坡底面不會出現(xiàn)大量積水。

圖7 各工況的吸力分布圖Fig. 7 Suction distribution diagram of various working conditions

2.4 雙層非飽和土邊坡破壞機理分析

降雨入滲對于非飽和土邊坡的影響主要在于:雨水滲入使得土體的重度增大,同時又使其基質(zhì)吸力降低甚至喪失。在雙重作用下,邊坡發(fā)生失穩(wěn)破壞。工況2的B層土的飽和度和其他幾個工況不同,因此不將其作為對比實驗進行分析。

2.4.1 邊坡破壞時間的分析

實驗和計算曲線終止的時間點為邊坡破壞的時刻,見圖6。3個工況破壞所需的時間從長到短依次為:工況1(250 min)、工況3(230 min)、工況4(116 min)。影響邊坡破壞時間的因素主要有2點:1)從圖1可知,模型1的坡角(59°)小于模型2(68.3°),從力學(xué)平衡分析,模型1邊坡的穩(wěn)定性高于模型2。2)從降雨條件分析,工況3的降雨強度(50 mm/h)小于工況4(100 mm/h),即相同的時間內(nèi),工況4的累計雨量大于工況3和工況1,從而加快邊坡的失穩(wěn)破壞。

2.4.2 邊坡塑性破壞帶形成過程的分析

圖8為邊坡破壞前后瞬時塑性剪切應(yīng)變分布的計算結(jié)果。從3個工況的云圖可知,在邊坡遭受破壞以后,都有一條明顯的塑性剪切帶貫穿整個邊坡。其所在的位置皆在兩土層交界面偏向A層一側(cè)。其原因是雨水滲入整個邊坡后,由于兩土層滲透性的差異,使得A層的滲透程度遠高于B層,即A層的吸力小于B層(如圖7所示),因此破壞帶主要在A層形成。從3個工況邊坡破壞前后的變化可知塑性破壞帶的演化過程,即塑性破壞從左側(cè)的坡趾附近開始,隨著降雨的進行逐漸蔓延至坡頂,形成一條滑動破壞帶,最終導(dǎo)致邊坡發(fā)生漸進性倒坍滑動破壞。

圖8 各工況邊坡模型破壞前后的塑性剪切應(yīng)變分布Fig. 8 Plastic shear strain distribution before and after failure of the slope model under various conditions

2.4.3 降雨強度對邊坡破壞影響分析

圖9為工況3和工況4邊坡破壞后的位移矢量圖。矢量的方向表示土體滑動方向。2個工況的左側(cè)土體整體向下滑動。可從圖9中很明顯地發(fā)現(xiàn)滑動帶的具體位置。工況4滑動的土體規(guī)模小于工況3。造成上述現(xiàn)象的原因是:從累計雨量來看,兩者大致相同(工況3降雨強度50 mm/h,歷時230 min;工況4降雨強度100 mm/h,歷時116 min)。從圖7可知,工況4中水從A層滲入B層的量相對較少,使得更多的水滯留在A層中,導(dǎo)致工況4的A層基質(zhì)吸力小于工況3,此現(xiàn)象也可在圖6中得到佐證。因此相比于工況3,其滑動帶的位置會發(fā)生上移。

圖9 破壞后的邊坡模型位移矢量分布Fig. 9 Displacement vector distribution of the slope model after failure

2.4.4 塑性剪切帶上單個單元的應(yīng)力路徑分析

圖10給出了工況3模型中某2個單元的應(yīng)力路徑。圖10(a)為選取單元的示意圖,選取臨界破壞面上C點和塑性帶外一點D。圖10(b)中橫軸表示平均有效應(yīng)力,縱軸表示廣義剪切應(yīng)力,隨著時間的推移,將獲得的C點和D點的平均有效應(yīng)力和剪切應(yīng)力大小在圖10(b)中繪出。從圖中可知,當(dāng)雨水滲入土體中后,土體內(nèi)部的孔壓逐漸上升,土骨架的有效應(yīng)力逐漸減小。2個單元的應(yīng)力路徑不斷向橫軸的負方向移動,與此同界破壞面上C點首先達到臨界狀態(tài)線,表明塑性破壞帶內(nèi)的單元已到達臨界狀態(tài)線,此時土體發(fā)生了破壞,當(dāng)C點達到臨界狀態(tài)線時,D點還未達到臨界狀態(tài)線,表明塑性帶外的單元沒有達到臨界破壞狀態(tài),并未產(chǎn)生明顯滑動。

圖10 塑性剪切帶上的單元位置及其應(yīng)力路徑Fig. 10 Element position and stress path on plastic shear band

3 結(jié)論

為研究不同降雨強度和坡度的雙層非飽和土邊坡的破壞機制,文中基于ZHANG等[18]提出的以飽和度和Bishop有效應(yīng)力為狀態(tài)變量的非飽和土本構(gòu)模型,將其導(dǎo)入土-水-氣三相耦合的有限元程序中,并對KITAMURA等[19]的室內(nèi)邊坡降雨破壞模型實驗進行數(shù)值模擬,得出了如下結(jié)論:

1)負孔隙水壓的計算結(jié)果較為準(zhǔn)確地擬合實驗數(shù)據(jù),表明文中所提的數(shù)值計算方法能在同一套參數(shù)下較好地模擬不同降雨強度和坡角的雙層非飽和土邊坡的破壞實驗。 計算結(jié)果得到了邊坡2個土層之間水分遷移的變化情況。開始階段,水由A層頂部滲入,隨后開始向A層土體的底部蔓延,水分由A層滲入B層較為緩慢。

2)降雨強度與坡角的不同會影響邊坡發(fā)生破壞的時間。坡角越小或降雨強度越小,邊坡越不易發(fā)生滑動破壞。計算結(jié)果揭示了非飽和土邊坡塑性破壞帶形成的過程。隨著降雨的進行,邊坡的塑性破壞首先在左側(cè)的坡趾附近形成,然后逐漸延伸至坡頂,形成一條貫穿整個邊坡的塑性剪切帶,最后邊坡左側(cè)的土體整體沿著滑動帶下滑。此外,還可得到降雨強度對邊坡滑動規(guī)模大小的影響,即在總降雨量大致相同時,低強度長歷時的降雨使邊坡滑動的土體體積更大。相比于雙層坡體,單層邊坡在A、B層交界面兩側(cè)的吸力變化相對平穩(wěn)。將B層置換成A層砂土后,整個坡體的抗?jié)B能力有了明顯的提升。

由于文中選取試驗為小尺寸的室內(nèi)實驗,模型的三邊皆為透氣透水邊界,且試驗所用的砂土透氣性較好。因此沒有探究氣壓對于坡體破壞的影響。若采用大尺寸的邊坡模型或者試驗土樣的透氣系數(shù)較小時,可能將產(chǎn)生較大的氣壓,今后,將會在這些方向上進行探究。