柴學(xué)銳，孫冬梅，S. Semprich

(1.天津大學(xué) 水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072；2.格拉茨技術(shù)大學(xué) 土力學(xué)及基礎(chǔ)工程研究所，奧地利 格拉茨 A-8010)

坡外部水位升高，岸坡土體內(nèi)的地下水位也隨之發(fā)生改變，原來(lái)被氣占據(jù)的孔隙空間逐漸被水所補(bǔ)充，在這種“水驅(qū)替氣”的過(guò)程中，孔隙壓力的變化將影響岸坡的穩(wěn)定性。因此，研究水位上升作用下岸坡的滲流場(chǎng)變化情況對(duì)進(jìn)一步研究岸坡穩(wěn)定性變化規(guī)律至關(guān)重要[1-2]。

水位上升作用下的岸坡是一種常見(jiàn)的飽和-非飽和土坡。在以往的飽和-非飽和滲流問(wèn)題中，一般僅研究水相而忽略氣相流動(dòng)，盡管這種僅考慮水相的方法簡(jiǎn)單易行并且已經(jīng)積累了一定經(jīng)驗(yàn)，但得出的結(jié)論與實(shí)際工況相比有很大出入。因此，為了更好地模擬非飽和區(qū)流體的運(yùn)移規(guī)律，國(guó)內(nèi)外有學(xué)者認(rèn)為采用水-氣二相流模型分析飽和-非飽和滲流問(wèn)題更符合實(shí)際[3-7]。

基于水-氣二相流模型，分析了水位上升過(guò)程中，飽和-非飽和土坡內(nèi)滲流場(chǎng)的變化規(guī)律；采用Bishop簡(jiǎn)化法[8]推導(dǎo)出同時(shí)考慮孔隙氣壓力、毛細(xì)壓力和坡外水壓力等貢獻(xiàn)的土坡安全系數(shù)計(jì)算公式，并計(jì)算了不同時(shí)刻邊坡的安全系數(shù)；分析安全系數(shù)的變化規(guī)律及各個(gè)力對(duì)安全系數(shù)的貢獻(xiàn)，為進(jìn)一步研究水位上升作用下土坡的失穩(wěn)機(jī)理提供理論依據(jù)。

1 水-氣二相流模型

近年來(lái)，多孔介質(zhì)的多相流理論研究和應(yīng)用越來(lái)越受到學(xué)者重視[9]。對(duì)于巖土工程滲流問(wèn)題而言，通常認(rèn)為滲流系統(tǒng)絕對(duì)恒溫，土體骨架不變形，處于平衡狀態(tài)的水相(l)和氣相(g)在多孔介質(zhì)中流動(dòng)，且滿足Darcy定律。

1.1 基本控制方程

當(dāng)非飽和土體中流體的擴(kuò)散作用忽略不計(jì)時(shí)，飽和-非飽和水-氣二相流模型的基本控制方程為

(1)

1.2 本構(gòu)方程

1.2.1 毛細(xì)壓力-飽和度關(guān)系模型

近年來(lái)，許多學(xué)者均提出了毛細(xì)壓力和飽和度的關(guān)系模型，如Brooks and Corey(BC)模型[10]和Van Genuchten(VG)模型[11]等，其中BC模型的進(jìn)氣值不具備連續(xù)性，而Van Genuchten(VG)模型不僅其表達(dá)形式簡(jiǎn)單，且模型的關(guān)系曲線中進(jìn)氣值和殘余水飽和度處過(guò)度平滑，因此在水-氣二相流模型中Van Genuchten(VG)模型得到廣泛應(yīng)用

(2)

式中：pc為水相和氣相交界面處的毛細(xì)壓力，等于孔隙水壓力和孔隙氣壓力的差值；P0為進(jìn)氣值；參數(shù)Sle為有效水飽和度，Sle=(Sl-Slr)/(Sls-Slr)，其中Sl為水飽和度，Slr為殘余水飽和度，Sls為飽和水飽和度；λ、ζ為模型的擬合參數(shù)，其值與土壤自身性質(zhì)有關(guān)；Pmax為最大毛細(xì)壓力。

1.2.2 相對(duì)滲透率和飽和度關(guān)系

土體孔隙中同時(shí)存在氣相流體和水相流體，一種流體流動(dòng)將會(huì)干擾另一種流體的流動(dòng)，當(dāng)水飽和度變化時(shí)，水相和氣相之間的交界面也將隨之變化，水流通道的彎曲度將會(huì)改變，導(dǎo)致相對(duì)滲透率發(fā)生改變。1976年，Mualem提出了相對(duì)滲透率和飽和度的關(guān)系表達(dá)式，將反應(yīng)毛細(xì)壓力和飽和度關(guān)系的VG模型表達(dá)式帶入相對(duì)滲透率和飽和度的關(guān)系表達(dá)式中，可得到TOUGH2中相對(duì)滲透率和飽和度關(guān)系的改進(jìn)模型，即VGM模型[12]

(3)

式中：krl為水相相對(duì)滲透率；krg為氣相相對(duì)滲透率；τ為迂曲因子，其取值大小與土體孔隙彎曲度有關(guān)，這里取0.5。

1.3 邊界條件處理方法

邊界條件實(shí)際上是指研究區(qū)域邊界所處的條件。在數(shù)學(xué)模型中，通過(guò)在滲流場(chǎng)邊界施加一薄層虛擬單元來(lái)實(shí)現(xiàn)各種邊界條件的數(shù)值模擬。其中，第一類邊界條件(Dirichlet條件)和第二類邊界條件(Neumann條件)為水-氣二相流模型根據(jù)液相和氣相劃分的邊界條件。

1.3.1 第一類邊界條件

1.3.2 第二類邊界條件

第二類邊界條件主要用來(lái)描述系統(tǒng)與外界的流量交換情況，且第二類條件可以是常量，也可以是隨時(shí)間變化的量。將上升引起的加載作用等效為邊界單元上有一個(gè)流入的流量，該值可通過(guò)式(4)求解

(4)

式中：m為源匯項(xiàng)，流入為正；b為擬合系數(shù)，擬合系數(shù)的取值與土壤的性質(zhì)有關(guān)；V為邊界單元的體積；C為邊界單元的壓縮性；ΔP為水壓力的變化量；Δt為水壓力變化所需的時(shí)間。以往的研究通常將水位變化過(guò)程等效為一系列特征水位，這種處理方法不能很好地模擬邊界上相態(tài)的轉(zhuǎn)變，而將水位上升過(guò)程等效為邊界單元上有一個(gè)流入的流量，可以實(shí)現(xiàn)水位的連續(xù)變化和邊界上相態(tài)的轉(zhuǎn)變。

1.4 初始條件

初始條件即為t=0時(shí)刻土體中各個(gè)狀態(tài)變量的分布情況。求解土體穩(wěn)定滲流問(wèn)題時(shí)，只需給出相應(yīng)的邊界條件即可，不需給定初始條件，但是求解土體非穩(wěn)定滲流問(wèn)題時(shí)，除了給出邊界條件之外，還需給定初始條件。對(duì)于坡外水位上升作用下的非穩(wěn)定滲流問(wèn)題，迭代求解的時(shí)間取決于給定的初始條件的合理性，當(dāng)初始條件接近真實(shí)值時(shí)，迭代收斂速度會(huì)加快，所需要的模擬時(shí)間會(huì)減少，因此初始條件的給出應(yīng)該以接近實(shí)際物理意義為準(zhǔn)則。

2 非飽和土坡穩(wěn)定性分析方法

2.1 非飽和土坡抗剪強(qiáng)度公式

Fredlund等[13]認(rèn)為非飽和土的抗剪強(qiáng)度由凈法向應(yīng)力、基質(zhì)吸力和有效黏聚力組成，由此提出了非飽和土的抗剪強(qiáng)度公式

τ=c′+(σn-pg)tanΦ′+(pg-pw)tanΦb

(5)

圖1 Bishop簡(jiǎn)化法下的土條作用力Fig.1 Forces acting on a slice in simplified Bishop′s Method

式中：τ為非飽和土的抗剪強(qiáng)度；c′為有效黏聚力；σn為總法向應(yīng)力；pg=pg-patm為孔隙氣壓力；pw=pw-patm為孔隙水壓力；Φ′表示有效內(nèi)摩擦角；Φb表示抗剪強(qiáng)度隨基質(zhì)吸力增加的曲線傾角；(σn-pg)為凈法向應(yīng)力；(pg-pw)為基質(zhì)吸力，其值等于負(fù)的毛細(xì)壓力[14]。

2.2 非飽和土坡安全系數(shù)計(jì)算公式

Bishop簡(jiǎn)化法常用來(lái)計(jì)算某一危險(xiǎn)圓弧滑動(dòng)面上的安全系數(shù)。圖1為一均質(zhì)非飽和土坡，土坡中任一土條i上的作用力如圖所示：bi為土條i的寬度；αi為土條i的底面與水平方向的夾角；βi為土坡側(cè)面與水平方向的夾角；Wi為土條i的重力；Ni為土條i底部的正法向應(yīng)力；Ti為土條i底部的切向力；hi為土條i的高度；pi為作用在土條i上的坡外水壓力。

根據(jù)土條i在垂直方向上受力平衡

(6)

根據(jù)圓心力矩平衡

(7)

式中：di為坡外水壓力pi到圓心的垂直距離，R為半徑。

根據(jù)極限平衡條件，土條i底面的抗剪力Ti可表示為

(8)

則非飽和土坡圓弧滑動(dòng)面上的安全系數(shù)F可表示為

(9)

3 算例分析

3.1 算例說(shuō)明

圖2 邊坡幾何形狀剖面圖(單位：m)Fig.2 Cross section of soil slope

圖2為一個(gè)各向同性的勻質(zhì)土坡，該土坡長(zhǎng)58.0 m，寬30.0 m，左側(cè)坡高25.0 m，右側(cè)坡高10.0 m，坡面傾角40°，坡內(nèi)坡外的水位水平，且水位高程均為10.0 m。為了保證模擬結(jié)果的精確程度，對(duì)土坡進(jìn)行網(wǎng)格剖分，若網(wǎng)格尺寸剖分過(guò)小，則計(jì)算時(shí)間較長(zhǎng)且計(jì)算結(jié)果精度并未有明顯提升，若網(wǎng)格尺寸剖分過(guò)大，將導(dǎo)致模擬結(jié)果精確程度降低，因此，根據(jù)筆者已有經(jīng)驗(yàn)將網(wǎng)格剖分成水平方向上從0.5 m到1.0 m不等，垂直方向的網(wǎng)格尺寸均為0.5 m。在土體處于完全飽和的狀態(tài)下，采用Bishop簡(jiǎn)化方法計(jì)算得到最危險(xiǎn)圓弧滑動(dòng)面的圓心坐標(biāo)為(43.0，31.1)，半徑為21.1 m，且其中一條半徑與水平方向垂直。在該土坡上設(shè)置一個(gè)追蹤面A-A，該追蹤面處于土坡的中間偏左位置，其水平坐標(biāo)是25.0 m。

表1 土體特性參數(shù)Tab.1 Parameters of different soils

整個(gè)模擬過(guò)程中，設(shè)定研究區(qū)處于恒溫系統(tǒng)中，表1給出土體的特性參數(shù)取值。

3.2 滲流場(chǎng)分析

初始穩(wěn)態(tài)的邊界條件包括地表(如圖2中ab邊所示)和坡外水位以上的右側(cè)邊界(如圖2中bc邊所示)為大氣邊界，坡外水位以下(如圖2中cd和de邊所示)為水相邊界，區(qū)域左側(cè)(如圖2中af邊所示)和底部(如圖2中ef邊所示)為不透水邊界，在此邊界條件下，運(yùn)行模型直至毛細(xì)壓力與重力平衡。水位上升過(guò)程中，坡外水位在52 h內(nèi)從10.0 m快速上升至21.2 m，上升速度為0.598×10-4m/s，該過(guò)程中，地表面和坡外水位以上的右側(cè)邊界仍為大氣邊界，區(qū)域左側(cè)和底部仍為不透水邊界，但原來(lái)在水位以上的右側(cè)氣相邊界隨著水位上升逐漸被水淹沒(méi)，轉(zhuǎn)變?yōu)樗噙吔纭?/p>

模擬過(guò)程中，在邊界單元上施加源匯項(xiàng)，以實(shí)現(xiàn)由于水位上升引起的加載作用導(dǎo)致的邊界單元上的水壓力變化和相態(tài)的轉(zhuǎn)變。初始的源匯項(xiàng)值可以通過(guò)式(4)計(jì)算得到，經(jīng)過(guò)試算，得出水位上升過(guò)程的源匯項(xiàng)值為2.70×10-44kg/s。在邊界單元上施加源匯項(xiàng)處理水位上升引起的加載作用可以實(shí)現(xiàn)邊界單元上水壓力的連續(xù)變化，且比以往將水位上升過(guò)程簡(jiǎn)化為一系列不連續(xù)的特征水位更符合實(shí)際。

在穩(wěn)定滲流狀態(tài)下，各個(gè)物理量的分布狀態(tài)如圖3所示，非飽區(qū)的孔隙氣壓力水頭為0，孔隙水壓力水頭小于0，因此毛細(xì)壓力水頭小于0，在飽和區(qū)，孔隙氣壓力水頭等于孔隙水壓力水頭，毛細(xì)壓力水頭為0。觀察圖4，隨著坡外水位上升，坡外水逐漸流入土坡內(nèi)部，原來(lái)被氣體填充的土壤孔隙逐漸被水填充，水飽和度逐漸增加，飽和區(qū)范圍逐漸上移，且毛細(xì)壓力水頭，孔隙水壓力水頭和孔隙氣壓力水頭均逐漸增加，其中接近土坡表面處的各物理量的變化更為顯著。由圖5表示的水位穩(wěn)定后滲流場(chǎng)各個(gè)物理量分布情況可知，飽和區(qū)范圍繼續(xù)上移，非飽和區(qū)的水飽和度繼續(xù)增加，毛細(xì)壓力水頭、孔隙水壓力水頭和孔隙氣壓力水頭也進(jìn)一步增加。

3-a 水飽和度3-b 毛細(xì)壓力3-c 孔隙水壓力3-d 孔隙氣壓力圖3 初始穩(wěn)態(tài)滲流場(chǎng)狀態(tài)分布Fig.3 Distributions of water saturation, capillary pressure and pore pressures in initial steady state

4-a 水飽和度4-b 毛細(xì)壓力4-c 孔隙水壓力4-d 孔隙氣壓力圖4 52 h滲流場(chǎng)狀態(tài)分布Fig.4 Distributions of water saturation, capillary pressure and pore pressures in 52 h

圖6為追蹤面A-A在初始穩(wěn)態(tài)，第52小時(shí)和第200小時(shí)的孔隙氣壓力、孔隙水壓力、毛細(xì)壓力和水飽和度分布情況。觀察圖6，水位上升以前，非飽和土體內(nèi)的孔隙氣壓力為0，但是由于毛細(xì)壓力的存在，在飽和土和非飽和土的接觸面上孔隙氣壓力有一個(gè)微小的波動(dòng)。水位上升過(guò)程中，隨著坡外水體流入土坡內(nèi)部，孔隙氣壓力的大小逐漸增加，水位穩(wěn)定后，坡外水體繼續(xù)流入土坡內(nèi)部，土體內(nèi)的孔隙氣壓力繼續(xù)增加。初始狀態(tài)下，土體內(nèi)的孔隙水壓力水頭隨著高程的降低呈現(xiàn)線性增加的趨勢(shì)，在水位上升過(guò)程中和水位穩(wěn)定后，隨著坡外水體流入土坡內(nèi)部，孔隙水壓力水頭均在逐漸增加。毛細(xì)壓力為孔隙水壓力與孔隙氣壓力的差值，且在水位上升作用下，毛細(xì)壓力值逐漸增加。水位上升前，土坡的水飽和度隨高程的降低逐漸升高，在水位上升過(guò)程中和水位穩(wěn)定后，由于坡外水流入土坡內(nèi)部，土坡的水飽和度繼續(xù)增加。

圖6 初始穩(wěn)態(tài),第52小時(shí)和第200小時(shí)追蹤面A-A上的孔隙氣壓力、孔隙水壓力、毛細(xì)壓力和水飽和度分布情況Fig.6 Simulation results of pore-air pressure, pore-water pressure, capillary pressure and water saturation at section A-A at initial steady state, 52 h and 200 h

3.3 飽和-非飽和土坡穩(wěn)定性分析

圖7為土坡安全系數(shù)隨時(shí)間變化過(guò)程，同時(shí)考慮孔隙氣壓力、毛細(xì)壓力、坡外水壓力和土體重力貢獻(xiàn)的安全系數(shù)的初始值為1.54，隨著水位的上升，安全系數(shù)逐漸增加，當(dāng)水位升至21.2 m時(shí)，安全系增加至最大值2.67，水位穩(wěn)定后，安全系數(shù)開(kāi)始逐漸下降。不考慮孔隙氣壓力貢獻(xiàn)的安全系數(shù)明顯大于考慮孔隙氣壓力和毛細(xì)壓力貢獻(xiàn)的安全系數(shù)，其值從初始時(shí)的1.54上升到2.86后開(kāi)始逐漸下降。不考慮毛細(xì)壓力貢獻(xiàn)的安全系數(shù)在三者中最小，其值從初始的1.06上升至2.21后開(kāi)始逐漸下降。水位上升過(guò)程中，對(duì)土坡穩(wěn)定性起主要作用的且有利于土坡穩(wěn)定的坡外水壓力逐漸增加，使得安全系數(shù)逐漸增加；當(dāng)水位穩(wěn)定后，坡外水壓力保持不變，而不利于土坡穩(wěn)定的孔隙氣壓力值逐漸增加使得安全系數(shù)逐漸降低。

圖7 危險(xiǎn)圓弧滑動(dòng)面上安全系數(shù) 隨時(shí)間變化曲線Fig.7 Time evolution of the safety factors on the slip surface

圖8為水位上升過(guò)程中的各個(gè)力對(duì)安全系數(shù)的貢獻(xiàn)及作用在危險(xiǎn)圓弧滑動(dòng)面上的各個(gè)力隨時(shí)間的變化情況。由于土體自身重力在水位上升初期為1 640.40 kN/m，到模擬結(jié)束后為1 665.00 kN/m，上升了24.60 kN/m，較其自身而言，增加的重力可以忽略不計(jì)，因此土體自身重力對(duì)安全系數(shù)的貢獻(xiàn)幾乎不變，將不做詳細(xì)討論。

圖8-a為水位上升作用下毛細(xì)壓力對(duì)安全系數(shù)的貢獻(xiàn)隨時(shí)間變化的曲線及作用在危險(xiǎn)圓弧滑動(dòng)面上的總毛細(xì)壓力隨時(shí)間的變化規(guī)律。觀察圖8-a，作用在危險(xiǎn)圓弧滑動(dòng)面上總的毛細(xì)壓力值小于0，隨著水位上升毛細(xì)壓力值逐漸下降，當(dāng)水位穩(wěn)定后，毛細(xì)壓力值緩慢上升，因此毛細(xì)壓力對(duì)安全系數(shù)的貢獻(xiàn)先逐漸增加后逐漸減小，所以，安全系數(shù)圖像中不考慮毛細(xì)壓力貢獻(xiàn)的安全系數(shù)與考慮孔隙氣壓力和毛細(xì)壓力貢獻(xiàn)的安全系數(shù)的差值先逐漸增加后逐漸減小。

圖8-b為水位上升作用下孔隙氣壓力對(duì)安全系數(shù)的貢獻(xiàn)隨時(shí)間變化的曲線及作用在危險(xiǎn)圓弧滑動(dòng)面上總的孔隙氣壓力隨時(shí)間的變化規(guī)律。如圖8-b所示，水位上升過(guò)程中作用在危險(xiǎn)圓弧滑動(dòng)面上總的孔隙氣壓力先快速增加，即從0 kN/m增加至805 kN/m，由于孔隙氣壓力不利于邊坡穩(wěn)定，總的孔隙氣壓力增加導(dǎo)致孔隙氣壓力對(duì)安全系數(shù)的不利貢獻(xiàn)從0上升至-0.16，當(dāng)水位穩(wěn)定后，總的孔隙氣壓力繼續(xù)由805 kN/m逐漸增加至1 068 kN/m，因此安全系數(shù)圖像中，不考慮孔隙氣壓力貢獻(xiàn)的安全系數(shù)與考慮孔隙氣壓力和毛細(xì)壓力貢獻(xiàn)的安全系數(shù)的差值先快速增加后緩慢增加。

圖8-c為水位上升作用下坡外水壓力對(duì)安全系數(shù)的貢獻(xiàn)隨時(shí)間變化的曲線及作用在危險(xiǎn)圓弧滑動(dòng)面上總的坡外水壓力隨時(shí)間的變化規(guī)律。水位上升過(guò)程中坡外水壓力從0 kN/m快速上升至960 kN/m，使得坡外水壓力對(duì)安全系數(shù)的貢獻(xiàn)從0快速增加至1.56，當(dāng)水位穩(wěn)定以后，由于坡外水位保持不變，因此作用在危險(xiǎn)圓弧滑動(dòng)面上總的坡外水壓力值始終為960 kN/m，坡外水壓力對(duì)安全系數(shù)的貢獻(xiàn)也始終穩(wěn)定在1.56。

8-a 毛細(xì)壓力8-b 孔隙氣壓力8-c 坡外水壓力圖8 滑動(dòng)面上各個(gè)力對(duì)安全系數(shù)的貢獻(xiàn)及各個(gè)力隨時(shí)間的變化過(guò)程Fig.8 Contributions to the safety factors from different kinds of pressures and different kinds of pressures acting on the slip surface

4 結(jié)論

通過(guò)建立水-氣二相流模型，研究了水位上升作用下土坡滲流場(chǎng)變化規(guī)律，并進(jìn)一步分析了土坡穩(wěn)定性變化規(guī)律，得到以下結(jié)論：

(1)將水位上升引起的飽和-非飽和滲流問(wèn)題作為水-氣二相流問(wèn)題研究更符合實(shí)際物理意義。從物理現(xiàn)象上來(lái)看，水位上升是土壤中“水驅(qū)替氣”的過(guò)程，因此，隨著水位上升，孔隙氣壓力、孔隙水壓力、毛細(xì)壓力和水飽和度均逐漸增加。

(2)岸坡的安全系數(shù)隨著水位的上升逐漸增加，水位穩(wěn)定后，安全系數(shù)逐漸下降，且考慮孔隙氣壓力和毛細(xì)壓力貢獻(xiàn)的安全系數(shù)大于不考慮毛細(xì)壓力貢獻(xiàn)的安全系數(shù)，但小于不考慮孔隙氣壓力貢獻(xiàn)的安全系數(shù)。因此，毛細(xì)壓力的存在有利于土坡穩(wěn)定，而孔隙氣壓力的存在不利于土坡穩(wěn)定。

(3)水位上升作用下，毛細(xì)壓力逐漸減小，水位穩(wěn)定后，毛細(xì)壓力逐漸增加，因此毛細(xì)壓力對(duì)安全系數(shù)的貢獻(xiàn)先逐漸增加后逐漸減?。凰簧仙^(guò)程中，孔隙氣壓力先快速增加，當(dāng)水位穩(wěn)定后，孔隙氣壓力緩慢增加，因此孔隙氣壓力對(duì)安全系數(shù)的不利貢獻(xiàn)先快速增加后緩慢增加；坡外水壓力逐漸增加后穩(wěn)定不變，因此，坡外水壓力對(duì)安全系數(shù)的貢獻(xiàn)先逐漸增加后保持不變。