非飽和參數(shù)交互作用的邊坡穩(wěn)定性分析

郭金鑫， 史 勇

(內(nèi)蒙古工業(yè)大學(xué)土木工程學(xué)院，呼和浩特 010051)

降雨入滲是造成自然邊坡以及工程邊坡失穩(wěn)的關(guān)鍵因素。構(gòu)成邊坡的土體結(jié)構(gòu)是非飽和土?xí)r，降雨入滲過程的滲流場與應(yīng)力場存在流固耦合作用，雨水在土體孔隙中做滲流運動，產(chǎn)生的滲透體積力會影響土體應(yīng)力場；土體由于受力變形，土體之間的孔隙發(fā)生變化，影響土體的滲透系數(shù)，進(jìn)而影響水在土體中的滲流[1-2]，土體抗剪強度降低并伴隨地下水位升高導(dǎo)致的孔隙壓力增大、穩(wěn)定地下水位以上區(qū)域呈暫態(tài)飽和狀態(tài)，這對邊坡整體的穩(wěn)定性造成影響[3]。因此，研究降雨入滲對邊坡穩(wěn)定性的影響具有重要意義。

中外學(xué)者對降雨入滲導(dǎo)致的邊坡失穩(wěn)開展了大量的研究，Wang等[4]建立了地下水滲流對排土場邊坡穩(wěn)定性影響的安全系數(shù)，且隨著地下水位的升高，邊坡穩(wěn)定系數(shù)逐漸減小，直至達(dá)到穩(wěn)定值。王述紅等[5]運用有限元軟件模擬了降雨條件下的邊坡滲流和應(yīng)力場分布，通過強度折減法對關(guān)鍵位移點和塑性應(yīng)變進(jìn)行了失穩(wěn)判斷。常金源等[6]基于Green-Ampt模型構(gòu)建了降雨入滲條件下淺層滑坡的概念模型，實現(xiàn)了降雨前不同地下水位狀態(tài)的邊坡安全系數(shù)與降雨時間的函數(shù)表征。王曼等[7]基于強度折減法實現(xiàn)了邊坡失穩(wěn)滑移面和滑移位置的計算，并與極限平衡法進(jìn)行了滑移特征點的位移和塑性區(qū)域?qū)Ρ?。宋亞萌等[8]利用Geostudio進(jìn)行了邊坡在不同非飽和參數(shù)作用下的失穩(wěn)分析，并基于灰色關(guān)聯(lián)度方法實現(xiàn)了降雨入滲過程的孔壓及安全系數(shù)變化規(guī)律研究。蔡欣育等[9]對基于不同類型降雨下的邊坡失效特征，利用飽和-非飽和滲流原理和有限元軟件對降雨類型及不同降雨強度下的邊坡滲流穩(wěn)定性進(jìn)行了模擬，確定了不同特征區(qū)的孔壓、安全系數(shù)變化。

上述研究在邊坡安全系數(shù)求解、考慮非飽和參數(shù)的單獨變化時邊坡穩(wěn)定性的影響方面取得了有益的成果[6,10-12]。然而，滲流場和應(yīng)力場耦合作用的本質(zhì)決定了非飽和參數(shù)的相互影響，因此，考慮滲流-應(yīng)力耦合作用、分析降雨入滲條件下非飽和參數(shù)對邊坡穩(wěn)定特性的直接作用和參數(shù)之間的耦合作用，將會進(jìn)一步完善已有研究成果的有效性。

基于此，通過降雨入滲條件下的瞬態(tài)滲流和應(yīng)力變化，確定邊坡孔壓和飽和度分布特征，然后，基于方差分析方法和強度折減法計算表征邊坡安全系數(shù)的位移值和應(yīng)力變化特征，最后基于瞬態(tài)降雨入滲條件，確定影響邊坡穩(wěn)定的非飽參數(shù)的滲透系數(shù)、降雨強度、黏聚力及內(nèi)摩擦角對其穩(wěn)定性的直接作用規(guī)律和參數(shù)間耦合特征。

1 非飽和流固耦合與方差理論

1.1 流固耦合機理

降雨入滲是典型的非飽和流固耦合過程，土體的非完全密實特性使得入滲水沿孔隙發(fā)生滲流運動，并產(chǎn)生滲流水壓且以滲透體積力的方式形成土體結(jié)構(gòu)的應(yīng)變場，進(jìn)而形成位移和體應(yīng)變的變化，降雨入滲影響土體的滲透系數(shù)和滲流場。由水力學(xué)基本原理可知[13]，連續(xù)多孔巖土體介質(zhì)的滲透體積力與水力梯度成正比例函數(shù)關(guān)系，可表示為

(1)

(2)

式中：γw為土體結(jié)構(gòu)和流體性質(zhì)決定的常數(shù)；h為降雨總水頭；f為降雨入滲的滲流體積力；fx、fy及fz分別為f在x、y、z向的分量；θ1、θ2及θ3分別為滲流體積力與分量值間的夾角；Jx、Jy及Jz分別為體單元在x、y、z向上的水力坡降。

在巖土工程的數(shù)值分析中，運用有限元法對應(yīng)力場求解時，單元的滲透體積力可通過式(3)、式(4)轉(zhuǎn)化為等效結(jié)點荷載。

(3)

(4)

式中：Ω為滲透體積；N為坡體單元形函數(shù)矩陣；F及ΔF分別為滲流體積力與增量值產(chǎn)生的等效節(jié)點力。

在進(jìn)行滲流應(yīng)力耦合分析時，采用孔隙率與滲透系數(shù)之間的經(jīng)驗關(guān)系式：

(5)

式(5)中：n0、n分別為邊坡土體初始狀態(tài)及降雨入滲后的孔隙率；k、k(δij)分別為入滲前、入滲后的滲透系數(shù)。

1.2 滲流-應(yīng)力場耦合模型

滲透體積力等外荷載引起土體應(yīng)力場的變化和土體的固結(jié)變形，其孔隙比和孔隙率發(fā)生變化，進(jìn)而影響其滲透系數(shù)，最終引起滲流場的變化[14]。

應(yīng)力場作用下的二維穩(wěn)定滲流場可表征為

(6)

式(6)中：[k(δij)]為節(jié)點形成的單元的滲流系數(shù)矩陣；H為總水頭；H1為Γ1的已知函數(shù)；n2、n3分別為Γ2、Γ3的法向；Γ1為水頭邊界；Γ2為流量邊界；Γ3為滲流域邊界；q為單位滲水量。

材料滲透系數(shù)及飽和度與基質(zhì)吸力的關(guān)系可表示為

Kw=awKws/[aw+(bwua-bwub)cw]

(7)

式(7)中：Kw為當(dāng)前滲透系數(shù)；Kws為土體飽和滲透系數(shù)；ua、ub分別為邊坡土體中的氣壓及水壓，其中，ua=0，ub=0；aw、bw及cw為材料系數(shù)。

Sr=Si+as(Sn-Si)/[as+(bsua-bsuw)cs]

(8)

式(8)中：Sr、Si分別為飽和度、殘余飽和度，其中，Si=0.08；Sn為最大飽和度，Sn=1。

滲流場作用下，應(yīng)力場的函數(shù)表征模型為

(9)

式(9)中：Ti為邊界Su上的面力分布；K為坡體剛度矩陣；ΔF為外載荷導(dǎo)致的節(jié)點載荷增量；ΔFs為滲流體積力引起的節(jié)點載荷增量；Δδ為位移變化；δ、δ0分別為節(jié)點位移和節(jié)點初始位移；Su、Sδ分別為位移和應(yīng)力邊界。

由式(8)、式(9)可知，滲流-應(yīng)力耦合的矩陣函數(shù)表達(dá)式為

(10)

式(10)中:f為降雨滲流場的水量分布。

利用直接耦合法將應(yīng)力場的節(jié)點位移和滲流場的孔隙水壓力作為待求解量，通過ABAQUS軟件的滲流—應(yīng)力耦合單元進(jìn)行迭代求解。

1.3 邊坡穩(wěn)定影響因素分析的基礎(chǔ)

土體飽和滲透系數(shù)、降雨強度、材料屬性的摩爾庫倫參數(shù)及摩擦角是決定土體邊坡穩(wěn)定的關(guān)鍵內(nèi)、外部因素。為分析以上4種參數(shù)對邊坡穩(wěn)定性的影響，需要對其進(jìn)行交互作用分析。利用輸入?yún)?shù)的變化方差對輸出結(jié)果變化方差的作用量，最終實現(xiàn)輸入?yún)?shù)及參數(shù)間的相互作用下對輸出參數(shù)的量化程度[15]。以邊坡模型在降雨入滲下的孔壓、飽和度、安全系數(shù)和有效應(yīng)力為輸出參數(shù)，以土體飽和滲透系數(shù)、降雨強度、材料屬性的摩爾庫倫參數(shù)及摩擦角為輸入?yún)?shù)進(jìn)行全局量化分析。參數(shù)定義為模型y=y(x1,x2, …,xn)的輸入?yún)?shù)，輸入?yún)?shù)i對輸出參數(shù)的作用關(guān)系可定義為

di(xi)=

(11)

式(11)中：di(xi)表示參數(shù)xi變化對輸出值di的影響程度，即輸出值對參數(shù)變化的敏感程度，其中，xi為輸入?yún)?shù)(土體飽和滲透系數(shù)、降雨強度、材料屬性的摩爾庫倫參數(shù)及摩擦角)；y為輸出參數(shù)(孔壓、飽和度、安全系數(shù)和有效應(yīng)力)；Δ為同一輸入?yún)?shù)變化前后的差值。

建立關(guān)于k個輸入?yún)?shù)的模型，i=1,2,…，k，輸入?yún)?shù)在p個級別上取值。通過多次抽樣代替重復(fù)取樣，并利用多次取樣計算確定的輸出參數(shù)的均值表征輸入?yún)?shù)對輸出的直接作用，利用統(tǒng)計量標(biāo)準(zhǔn)差表征輸入?yún)?shù)間的交互作用。

2 邊坡穩(wěn)定性有限元建模

2.1 邊坡穩(wěn)定性模型及邊界條件

對某工程實地邊坡選取坡面，利用ABAQUS軟件建立二維平面模型。有限元網(wǎng)格模型及幾何尺寸如圖1所示。

坡體土質(zhì)的彈性模量是10 MPa，泊松比0.3，土體的飽和滲透系數(shù)Kws為0.018 m/h。初始孔隙比為1.0。勻質(zhì)土坡的坡高為25 m，坡高比為1∶1.5。圖1中AB和EG為不透水邊界，DC和CB為降雨入滲邊界，施加的降雨強度初始值為20 mm/h。坡度體的網(wǎng)格劃分單元為多孔介質(zhì)流固耦合單元CPE4P。

根據(jù)模型輸入?yún)?shù)土體飽和滲透系數(shù)、降雨強度、材料屬性的摩爾庫倫參數(shù)及摩擦角的初始值，不飽和參數(shù)Kws為0.018，降雨強度為20 mm/h，黏聚力c為15 kPa，內(nèi)摩擦角φ為30°，4個參數(shù)的設(shè)定變化范圍分別為±10%、±10%、±10%和±5%?；贛orris量化分析方法和統(tǒng)計學(xué)原理可知4項參數(shù)有15種組合方式。如表1所示。利用不同參數(shù)組合及分布區(qū)間取值，進(jìn)行邊坡穩(wěn)定模型計算，可確定4項輸入?yún)?shù)對坡體穩(wěn)定特征的作用值。

圖1 坡體有限元模型及尺寸示意Fig.1 Finite element model and dimension of slope

表1 輸入?yún)?shù)及波動區(qū)間

2.2 邊坡穩(wěn)定性模型求解

降雨入滲過程中，采用折減法對土體邊坡穩(wěn)定進(jìn)行分析時，即計算庫倫參數(shù)和摩擦角變化導(dǎo)致的單元應(yīng)力大小，且與屈服應(yīng)力進(jìn)行比較，當(dāng)大于屈服應(yīng)力時單元失穩(wěn)，土體達(dá)到極限破壞形成連續(xù)滑動面[16-17]。土體的摩爾庫倫參數(shù)及摩擦角與折減后抗剪強度函數(shù)關(guān)系可表達(dá)為

(12)

式(12)中：Fr為強度折減系數(shù)；c、φ分別為折減前的邊坡土體黏聚力和內(nèi)摩擦角；cm、φm分別為折減后的邊坡土體黏聚力和內(nèi)摩擦角。

基于強度折減法進(jìn)行數(shù)值計算時，當(dāng)不飽和參數(shù)Kws=0.018，φ=30°，c=16.5,q=20 mm/h時，邊坡穩(wěn)定性特征如圖2所示。

圖2 非飽和參數(shù)作用下的邊坡穩(wěn)定性仿真結(jié)果Fig.2 Simulation results of slope stability under unsaturated parameters

由圖2可知，當(dāng)不飽和參數(shù)Kws=0.018，φ=30°，c=16.5,q=20 mm/h時，土體邊坡的水平和豎直方向位移隨抗剪切強度值減小而增大。當(dāng)邊坡達(dá)到失穩(wěn)狀態(tài)時，被選特征單元如坡腳和坡頂特征點的位移發(fā)生突變。孔壓的分布呈線性特征，坡體的底部孔壓為110.6 MPa，頂部為122.7 MPa。由圖2(b)可知，飽和度的分布呈線性特征，孔壓小于0的區(qū)域是非飽和土體區(qū)，其飽和度由下往上迅速減小，孔壓為負(fù)值，最大飽和度為1.08；孔壓大于0的區(qū)域是飽和土體區(qū)，其最大飽和度為0.08。由圖2(c)可知，滑移面通過坡腳點呈圓弧狀分布，在邊坡頂部滑動位移值最大，沿圓弧分布線呈遞減趨勢變化。由圖2(d)可知，塑性變形區(qū)域出現(xiàn)在坡腳位置，并具有由坡面向里逐漸增加的趨勢，最大值位于邊坡的最底部，符合應(yīng)力分布的規(guī)律。

3 邊坡失穩(wěn)參數(shù)作用分析

3.1 單參數(shù)作用下的穩(wěn)定性分析

基于強度折減法可確定邊坡穩(wěn)定的安全系數(shù)，初始狀態(tài)強度折減系數(shù)為1。將由土體抗剪強度參數(shù)可計算獲得初始狀態(tài)的單元位移和應(yīng)力。當(dāng)土體飽和滲透系數(shù)Kms在其取值區(qū)間內(nèi)均勻分布時(參數(shù)取值同2.1節(jié))，孔壓、飽和系數(shù)、特征點位移及應(yīng)力分布如圖3所示。

圖3 飽和滲透系數(shù)對邊坡穩(wěn)定性影響Fig.3 Influence of saturated permeability coefficient on slope stability

由圖3可知，頂部孔壓隨著土體飽和滲透系數(shù)的增大而降低，土體飽和強度隨滲透系數(shù)的增大而增大，這與土體強度本質(zhì)特征相符合。特征點的總位移和等效塑性應(yīng)變隨土體飽和系數(shù)呈降低趨勢，且突變的特征點位移位置處等效塑性應(yīng)變也較大，達(dá)到0.012 4。

圖4 降雨強度對邊坡穩(wěn)定性影響Fig.4 Influence of rainfall intensity on slope stability

由圖4可知，在4項參數(shù)中僅有降雨強度發(fā)生變化時(參數(shù)取值同2.1節(jié))，頂部孔壓隨著降雨強度增大而降低，這與圖2顯示的土體孔隙的吸力隨飽和度的增大而降低特征一致。圖4中，特征點的總位移和等效塑性應(yīng)變隨土體飽和系數(shù)呈上升趨勢，在降雨強度最大時特征點位移和等效塑性應(yīng)變最大，分別達(dá)到11.55 mm和0.018%。

圖5為黏聚力對邊坡穩(wěn)定性影響，即此時其他3項參數(shù)為恒定值。由圖5可知，土體抗剪強度參數(shù)的摩爾庫倫值增大，孔壓與水位飽和度恒定，這是由于隨著黏聚力增大，土體抗剪強度增加，降雨入滲條件下土體抗下滑能力增強，特征點位移降低到8.5 mm，特征單元塑性變形減小，最小為0.005%，此時邊坡穩(wěn)定性增強。

圖5 黏聚力對邊坡穩(wěn)定性影響Fig.5 Influence of cohesion on slope stability

由圖6可知，土體抗剪強度參數(shù)的內(nèi)摩擦角增大，受降雨強度恒定值的作用，孔壓與水位飽和度保持恒定。這是由于隨著內(nèi)摩擦角的增大，土體的抗剪強度增加，土體結(jié)構(gòu)單元間的抗滑移能力增強，特征體單元間的位移降低至8.2 mm，塑性變形程度減小。

圖6 內(nèi)摩擦角的對邊坡穩(wěn)定性影響Fig.6 Effect of internal friction angle on slope stability

3.2 參數(shù)交互作用下的穩(wěn)定性分析

非飽和滲流狀態(tài)下，決定流固耦合作用的外部降水量的土體抗剪強度參數(shù)不是單獨作用，參數(shù)間存在交互影響并使表征坡體穩(wěn)定特征的孔壓、飽和度、位移及塑性應(yīng)變發(fā)生變化。依據(jù)Morris確定的輸入?yún)?shù)組合進(jìn)行有限元仿真分析可確定15組參數(shù)作用下，4種穩(wěn)定性表征參數(shù)的變化特點。

由圖7(a)可知，15組參數(shù)中的前4組(1、2、3、4)是4項參數(shù)的單獨作用，橫坐標(biāo)是參數(shù)組合的標(biāo)準(zhǔn)差，表示參數(shù)間的耦合交互作用?？v坐標(biāo)是參數(shù)組合的均值，表示參數(shù)對邊坡穩(wěn)定性孔壓特征的直接作用。參數(shù)組合2的均值及標(biāo)準(zhǔn)差>參數(shù)組合1>參數(shù)組合3>參數(shù)組合4，即降雨強度對邊坡穩(wěn)定性孔壓特征的影響最為明顯。參數(shù)組6、7、8、9的標(biāo)準(zhǔn)差>參數(shù)組合5的均值和標(biāo)準(zhǔn)差，即前者的參數(shù)組合形式對邊坡穩(wěn)定性的孔壓特征有正向作用，參數(shù)組合4對孔壓特征具有負(fù)向作用，與參數(shù)組合1形成組合5后二者之間的交互耦合作用降低了對邊坡穩(wěn)定性孔壓特征的影響。

由圖7(b)可知，參數(shù)組合2的均值及標(biāo)準(zhǔn)差>參數(shù)組合1>參數(shù)組合3>參數(shù)組合4，即降雨強度對邊坡穩(wěn)定性孔壓特征的影響最為明顯，且降雨強度與其他3項參數(shù)的耦合交互關(guān)系較弱。

雙參數(shù)組合的直接作用效果顯示，參數(shù)組合由高到低排序為參數(shù)組合6>參數(shù)組合7>參數(shù)組合5>參數(shù)組合10>參數(shù)組合8>參數(shù)組合9，其中最小的兩組參數(shù)組合8和參數(shù)組合9均含有共同的降雨強度。結(jié)合單參數(shù)分析時的直接作用效果可知，降雨強度對飽和度的影響是負(fù)作用。由標(biāo)準(zhǔn)差可知參數(shù)間的耦合作用由高到低分別是參數(shù)組合8>參數(shù)組合9>參數(shù)組合5>參數(shù)組合10>參數(shù)組合6>參數(shù)組合7。

由圖7(c)可知，均值代表的直接作用效果中，參數(shù)組合2>參數(shù)組合4>參數(shù)組合3>參數(shù)組合1，標(biāo)準(zhǔn)差由大至小為參數(shù)組合4>參數(shù)組合3>參數(shù)組合2>參數(shù)組合1，即飽和系數(shù)與其他3項之間的耦合作用最小，內(nèi)摩擦角與其他3項的耦合作用對邊坡穩(wěn)定性位移特征影響最大。

由雙參數(shù)組合的均值分布可知，參數(shù)組合10>參數(shù)組合6>參數(shù)組合5>參數(shù)組合8>參數(shù)組合9>參數(shù)組合7。由單參數(shù)的均值可知，飽和滲透系數(shù)與內(nèi)摩擦角之間為負(fù)向耦合作用，因此導(dǎo)致二者參數(shù)組合7對位移變化的直接作用和耦合作用降低。多參數(shù)作用時，降雨強度、黏聚力及摩擦角之間是正耦合作用，因此參數(shù)組合13均值及標(biāo)準(zhǔn)差較大，三者與飽和滲透系數(shù)共同作用后，參數(shù)組合15對位移變化的直接作用效果和耦合作用效果明顯降低，即飽和滲透系數(shù)對邊坡穩(wěn)定的特征點位移變化是負(fù)作用。

由圖7(d)可知，均值表征的直接作用效果上，參數(shù)組合3>參數(shù)組合2>參數(shù)組合4>參數(shù)組合1，標(biāo)準(zhǔn)差依次為參數(shù)組合3>參數(shù)組合4>參數(shù)組合1>參數(shù)組合2，綜合均值與標(biāo)準(zhǔn)差系數(shù)可知，飽和滲透系數(shù)對邊坡穩(wěn)定的等效塑性應(yīng)變影響最弱。

圖7 參數(shù)交互作用對孔壓、飽和度、位移及等效塑性 應(yīng)變的影響Fig.7 Influence of parameter interaction on pore pressure, saturation, displacement and equivalent plastic strain

由雙參數(shù)組合的均值可知，參數(shù)組合6>參數(shù)組合10>參數(shù)組合5>參數(shù)組合8>參數(shù)組合9>參數(shù)組合9，標(biāo)準(zhǔn)差由高到低排序為參數(shù)組合10>參數(shù)組合6>參數(shù)組合8>參數(shù)組合9>參數(shù)組合7>參數(shù)組合5，綜合均值與標(biāo)準(zhǔn)差可知，飽和滲透系數(shù)、內(nèi)摩擦角及黏聚力組合對邊坡穩(wěn)定的塑性應(yīng)變特征作用較其他雙參數(shù)組合較為顯著，且黏聚力與二者之間呈正耦合作用。降雨強度與其他3項參數(shù)之間呈負(fù)向作用。

4項參數(shù)組合15作用時，對邊坡穩(wěn)定的塑性應(yīng)變特征、孔壓變化特征、飽和度特征及特征單元位移作用結(jié)果顯示，參數(shù)組合15的直接作用和耦合作用均處于較小水平，這表明四項參數(shù)間的共同波動時，交互耦合作用不是單項參數(shù)的作用效果代數(shù)堆疊。

4 結(jié)論

(1)單參數(shù)作用下，邊坡穩(wěn)定性特征孔壓、飽和度、位移及塑性應(yīng)變隨著飽和滲透系數(shù)的增大而降低，穩(wěn)定性特征隨降雨強度的增大而變大。黏聚力和內(nèi)摩擦角的增大對飽和度及孔壓作用不顯著，但邊坡穩(wěn)定特征增強。

(2)參數(shù)交互作用下，降雨強度和滲透系數(shù)對邊坡穩(wěn)定性孔壓及飽和度特征的直接影響和間接影響較內(nèi)摩擦角和黏聚力更為明顯。降雨強度對特征單元位移直接影響最大，內(nèi)摩擦角與其他三組參數(shù)的耦合作用最顯著。黏聚力對邊坡穩(wěn)定的直接影響最大，且對塑性變形作用時與其他參數(shù)的耦合最為緊密。

(3)土體飽和滲透系數(shù)、黏聚力和內(nèi)摩擦角及外部降雨強度對邊坡穩(wěn)定的孔壓、飽和度、位移及塑性應(yīng)變的直接作用和耦合作用均處于較小水平，4項參數(shù)間的共同變化時，交互耦合作用不等于單項參數(shù)的作用效果代數(shù)堆疊。