徑流-滲流作用下松散土坡非均勻顆粒臨界失穩(wěn)分析

劉 悅，袁星宇，張 帥，符文熹，李龍國

（1.四川大學水利水電學院，成都610065；2.四川大學水力學與山區(qū)河流保護國家重點實驗室，成都610065）

0 引 言

松散土坡淺表物質(zhì)具有結(jié)構(gòu)松散、膠結(jié)程度低、抗干擾能力弱、孔隙度大、透水性強等特性，斜坡顆粒極易在降雨、地震、人類活動等因素影響下發(fā)生失穩(wěn)破壞［1，2］。我國西南地區(qū)松散土坡分布廣泛，降水豐富，發(fā)生滑坡和碎屑流災害嚴重。因此，有必要開展坡面徑流與坡體滲流聯(lián)動作用下松散土坡顆粒沖刷啟動機理研究。

近年來，松散土坡穩(wěn)定性問題受到國內(nèi)外學者廣泛關注。Rijn 等［3］基于泥沙顆粒的受力特征，通過綜合考慮橫縱向底坡對泥沙起動的疊加效應，建立了任意坡面上泥沙顆粒處于臨界起動狀態(tài)時的拖曳力公式；謝立全等［4］在進行坡面泥沙受力分析時考慮了滲透力作用，發(fā)現(xiàn)向坡外的滲流將大幅降低岸坡泥沙起動條件；Kumar等［5］考慮到岸坡滲流的影響，提出無量綱的水流功率概念，并將其運用到滲流對坡面泥沙起動過程的描述中；Zamani 等［6］通過顆粒離散元法建立了三維土-基礎-結(jié)構(gòu)體系，對松散體在不同多向振動作用下的響應進行模擬；田海等［7］進行離心模型試驗，對降雨條件下邊坡土體強度的變化進行了分析，并得出結(jié)論降雨會降低松散堆積體邊坡穩(wěn)定性；毛雪松等［8］基于Midas/GTS 建立松散堆積體路基上邊坡穩(wěn)定數(shù)值模型，分析了各坡形因素對邊坡穩(wěn)定性的影響；張巖巖［9］采用一維地表徑流運動波方程和二維地下水滲流Richards 方程，將地表徑流和地下水滲流耦合，并結(jié)合非飽和土Bishop 有效應力及極限平衡理論，研究了降雨作用下滑坡穩(wěn)定性問題；羅浩等［10］考慮棄渣場坡體結(jié)構(gòu)特征，通過建立邊坡顆粒流模型對降雨工況下邊坡變形失穩(wěn)過程進行了模擬；雷小芹等［11］構(gòu)建描述非飽和堆積土中細顆粒侵蝕-運移-沉積全過程的滲流潛蝕模型，定量研究了細顆粒運移對邊坡降雨穩(wěn)定性的影響。

綜合查閱的文獻分析表明，有關坡面徑流和坡體滲流對松散土坡失穩(wěn)的聯(lián)合作用鮮有報道。本文從單顆粒的滑動和滾動角度進行分析，分別采用Navier-Stokes方程［12］和擴展Brinkman-Darcy方程［13，14］描述坡面徑流和坡體滲流，聯(lián)合連續(xù)性方程，構(gòu)建徑流-滲流作用下的分析模型；將徑流-滲流界面視為應力跳躍邊界，推求出流固界面流速表達式，并結(jié)合顆粒的滑動和滾動受力分析，得到考慮徑流-滲流聯(lián)動作用下單顆粒穩(wěn)定安全系數(shù)；最后討論了單顆粒滑動、滾動安全系數(shù)對各影響因子的敏感強弱，確定安全系數(shù)敏感參數(shù)。

1 徑流-滲流作用下非均勻顆粒滑動和滾動受力分析

1.1 顆粒滑動分析

松散土坡受外界因素（如降雨）影響，表層顆粒會發(fā)生運動。本文假定土坡顆粒為不可壓縮非均勻顆粒（等效粒徑為di），顆粒只沿土坡傾向方向運動（圖1）。

如圖1 所示，徑流-滲流作用下單個顆粒受到拖曳力FD、上舉力FL、水下重力W，以及由于顆粒相對暴露度產(chǎn)生的附加質(zhì)量力FM；FD、FL、W和FM的計算分別見式（1）~式（4）［15，16］。

圖1 松散土坡淺層顆粒受力情況Fig.1 Force of shallow particles on loose soil slope

式中：CD為拖曳力阻力系數(shù)；ρ為水的密度，kg/m3；αD為與拖曳力相應的面積系數(shù)，取值范圍0~1；ub為流固界面流速，m/s。

式中：CL為上舉力系數(shù)；αL為與上舉力相應的面積系數(shù)，取值范圍0~1；其余符號同前。

式中：αW為與水下重力相應的面積系數(shù)，取值范圍0~1；γs為泥沙顆粒的容重，N/m3；γw為水的容重，N/m3；其余符號同前。

式中：dmax、dmin、dm分別為顆粒級配中最大、最小和平均粒徑，m；αM為與附加質(zhì)量力相應的面積系數(shù)，取值范圍為0~1；di為研究顆粒等效粒徑，m；dj為研究對象周圍顆粒的粒徑，m；f（dj）為河床表面混合物以重量計算dj顆粒分布函數(shù)；ζ為與顆粒相對暴露度（床面粗化度）有關的系數(shù)，顆粒完全暴露時ζ=0，顆粒完全隱蔽時ζ=1。

當土坡顆粒發(fā)生滑動時，顆粒還會受到周圍其他顆粒的摩阻力Ff，F(xiàn)f計算見式（5）。

式中：φ為土體內(nèi)摩擦角；其余符號同前。

將式（3）、式（4）代入式（5）中，F(xiàn)f可改寫為式（6）。

在此，定義顆粒抗滑穩(wěn)定安全系數(shù)SFs為：

式中：Frf為抗滑力，N；Fsf為滑動力，N。

從圖1 可見，抗滑力僅為摩擦阻力Ff，即式（6）。滑動力包含水下重力沿x方向的分量Wx=Wsinθ、附加質(zhì)量力沿x方向的分量FMx=FMsinθ和拖曳力FD。于是，式（7）可具體寫出如下：

分析式（8）可知，顆粒抗滑穩(wěn)定安全系數(shù)主要受到流固界面流速ub、土坡坡度θ、土體內(nèi)摩擦角φ、顆粒等效粒徑di、顆粒容重γs的影響。

1.2 顆粒滾動分析

在此，定義顆粒抗?jié)L動穩(wěn)定安全系數(shù)SFr為：

式中：Frm為抗?jié)L動力矩，N?m；Fsm為滾動力矩，N?m。

從圖2 可見，當土坡顆粒發(fā)生滾動時，支持力FN=0，設顆粒繞O點發(fā)生滾動，且W、FD、FL、FM的力臂分別為lW、lD、lL、lM。于是，F(xiàn)rm=W lw+FMlM，F(xiàn)sm=FDlD+FLlL，相應地，式（9）可具體寫出如下。

圖2 顆粒滾動趨勢下受力及力臂情況Fig.2 Force and force arm of rolling particles

分析式（10）可知，顆粒抗?jié)L動穩(wěn)定安全系數(shù)主要受到流固界面流速ub、顆粒等效粒徑di、顆粒容重γs的影響。

2 徑流-滲流作用土坡流場分析

2.1 模型描述

根據(jù)2.1 節(jié)、2.2 節(jié)顆粒安全系數(shù)計算公式可知，徑流-滲流聯(lián)合作用下松散土坡非均勻顆粒是否穩(wěn)定與流固界面流速ub密切相關；安全系數(shù)隨著ub的增加而減小。為推求松散土坡流速分布、求解流固界面流速ub，本文將坡面徑流及坡體滲流進行耦合分析，如圖3 所示，土坡坡長為L，松散層厚度為b，孔隙率為n，滲透率為K，下層基巖假設為不透水層；松散土坡坡面徑流、坡體滲流均為二維平面流，且沿y方向的流速為0；坡面徑流水深為h，沿x方向的水流流速為ux；坡體滲流沿x方向的水流流速為vx。

圖3 松散土坡徑流-滲流耦合分析模型Fig.3 Runoff-Seepage coupling analysis model on loose soil slope

假設圖3 中水流是不可壓縮的Newton 流體，為充分發(fā)展層流。采用連續(xù)性方程［式（11）］和Navier-Stokes 方程［式（12）］描述松散土坡坡面徑流運動［12，17］。

式中：ux、uy分別為松散土坡坡面徑流沿x、y方向的流速，m/s；fx為沿x方向的質(zhì)量力，N/kg；P為沿x方向的壓強，Pa；υ為水的運動黏滯系數(shù)，m2/s；t為時間，s；其余符號同前。

已知水流沿y方向的流速為0，故有?uy/?y=0，將其代入式（11）中可得?ux/?x=0；x方向流速ux在y方向不發(fā)生變化，即?ux/?y=0；水流為恒定流，則?ux/?t=0；沿x方向fx=gsinθ，?P/?x=-ΔP/L。因此，式（12）可寫為：

式中：μ為水的動力黏滯系數(shù)，Pa?m，其中μ=υ ρ；L為土坡長度，m。

采用連續(xù)性方程［式（14）］和擴展Brinkman-Darcy 方程［式（15）］描述松散土坡坡體滲流運動［13，14，18］。

式中：vx、vy分別為松散土坡坡體滲流沿x、y方向的流速，m/s；n為土顆粒孔隙率；K為土坡滲透率，m2，k=K g/υ，其中k為滲透系數(shù)，m/s，g為重力加速度，m/s2。

同理，式（15）可簡化為：

2.2 界面流速

為便于計算，現(xiàn)對流速計算公式進行無量綱化處理。引入無量綱參數(shù)c0、c1、c2、c3、c4、V，令c0=b/D、c1=y/D、c2=1/n、c3=K/D2、c4=1/（c2c3）0.5、V=μ v/（ρ g D2tanθ），其中D=b+h，v為水流流速。則式（13）、式（16）可進一步簡化為：

式中：Vr為松散土坡坡面徑流的無量綱流速；Vs為松散土坡坡體滲流的無量綱流速。

求解式（17）、式（18），得無量綱流速表達式為：

式中：a1、a2、b1、b2為待求系數(shù)。

將無量綱參數(shù)V=μ v/（ρ g D2tanθ）代入式（19）和式（20）中，即可反求出有量綱流速：

松散土坡坡面徑流ux和坡體滲流vx滿足以下邊界條件［19，20］：

（1）在松散土坡坡面徑流上表面（y=b+h）處，流速ux達到最大，即滿足?ux/?y=0；

（2）在土坡表面（y=b）處，滿足速度連續(xù)且剪應力跳躍的邊界條件，即ux=vx、c2?vx/?c1-?ux/?c1=β vx/c30.5，其中β為應力跳躍系數(shù)，與多孔介質(zhì)特性、界面位置密切相關，其取值范圍為-1～1.5；

（3）在土坡底部（y=0）處，固體顆粒與表面流體不存在相對位移，即滿足ux=0。

將上述邊界條件代入式（21）、式（22）中，解得待定系數(shù)a1、a2、b1、b2，見式（23）～（26）。

式（24）～（26）中d1、d2、d3、d4、d5、d6分別見式（27）～（32）。

因此，流固界面流速ub（y=b）為：

分析式（33）可知，流固界面流速主要受到松散層厚度b、坡面徑流水深h、土坡坡度θ、應力跳躍系數(shù)β、松散土坡孔隙率n及滲透率K的影響。

3 顆粒穩(wěn)定性影響因子分析

根據(jù)式（8）、式（10）和式（33），松散土坡非均勻顆粒抗滑穩(wěn)定安全系數(shù)SFs的影響因素有松散層厚度b、坡面徑流水深h、土坡坡度θ、應力跳躍系數(shù)β、土體內(nèi)摩擦角φ、顆粒等效粒徑di、顆粒容重γs、松散土坡孔隙率n、滲透率K；抗?jié)L動穩(wěn)定安全系數(shù)SFr的影響因素有b、h、θ、β、n、K、di及γs。

為定量了解上述因素對松散土坡非均勻顆粒穩(wěn)定性的影響程度，對各因子進行敏感性分析。敏感系數(shù)Si計算見式（34），敏感系數(shù)越大，安全系數(shù)對該影響因子越敏感［21，22］；敏感分析相關參數(shù)基準值如表1所示，其中拖曳力阻力系數(shù)CD、上舉力系數(shù)CL、顆粒所受各力的面積系數(shù)α及力臂l取值參考文獻［23］、［24］。

表1 相關計算參數(shù)基準值［23，24］Tab.1 The base value of related calculation parameters

式中：Si為安全系數(shù)影響因子xi的敏感系數(shù)；|ΔSFi/SF0|為安全系數(shù)SFi的相對變化率；SF0為各參數(shù)基準值計算出的安全系數(shù)基準值；|Δxi/x0|為影響因子xi的相對變化率；x0為各參數(shù)基準值。

3.1 坡面徑流水深h敏感性

令坡面徑流水深h從0.6 mm 增加到1.8 mm，其余參數(shù)取值如表1 所示。根據(jù)不同的h值，求得與h對應的SFs、SFr值，徑流水深h與安全系數(shù)SF關系如圖4所示。

圖4 安全系數(shù)與敏感系數(shù)隨徑流深變化Fig.4 Variations of safety factors and sensitivity factors with runoff depth of slope surface

從圖4 可以看出，當h從0.6 mm 增加到1.8 mm 時，安全系數(shù)SFs從1.36 下降到0.47，SFr從4.60 下降到0.88，SFs、SFr隨著h的增加而降低；敏感系數(shù)Shs從1.13 降低到0.46，Shr從3.75 降低到0.39，Shs、Shr隨著h的增加而降低。當h一定時，顆粒更易發(fā)生滑動失穩(wěn)。

取h=1.0 mm 時的敏感系數(shù)為影響因子h的敏感系數(shù)，即Shs=0.90、Shr=1.51，明顯SFr相較于SFs對h更敏感。

3.2 坡度θ敏感性

令土坡坡度θ從10°增加到20°，其余參數(shù)取值如表1 所示，計算得與θ對應的SFs、SFr值。土坡坡度θ與安全系數(shù)SF關系如圖5所示。

從圖5 可以看出，當θ從10°增加到20°時，安全系數(shù)SFs從1.85下降到0.56，SFr從5.11下降到1.26，SFs、SFr隨著θ的增加而降低；敏感系數(shù)Sθs從3.94降低到0.96，Sθr從5.93降低到0.92，Sθs、Sθr隨著θ的增加而降低。當θ一定時，顆粒更易發(fā)生滑動失穩(wěn)。

取θ=15°時的敏感系數(shù)為影響因子θ的敏感系數(shù)，即Sθs=1.77、Sθr=2.04。

3.3 土體內(nèi)摩擦角φ敏感性

令顆粒內(nèi)摩擦角φ從25°增加到40°，其余參數(shù)取值如表1所示，計算得與φ對應的SFs、SFr值。內(nèi)摩擦角φ與安全系數(shù)SF關系如圖6所示。

圖6 安全系數(shù)與敏感系數(shù)隨土體內(nèi)摩擦角變化Fig.6 Variations of safety factors and sensitivity factors with soil internal friction angle

從圖6 可以看出，當φ從25°增加到40°時，安全系數(shù)SFs從0.78上升到1.40，SFr恒為2.26，SFs隨著φ的增加而增大；敏感系數(shù)Sφs從1.11升高到1.52，Sφs隨著φ的增加而增大。當φ一定時，顆粒更易發(fā)生滑動失穩(wěn)。

取φ=30°時的敏感系數(shù)為影響因子φ的敏感系數(shù)，即Sφs=1.21。

3.4 顆粒容重γs敏感性

令顆粒容重γs從16 kN/m3增加到30 kN/m3，其余參數(shù)取值如表1 所示，計算得與γs對應的SFs、SFr值。顆粒容重γs與安全系數(shù)SF關系如圖7所示。

圖7 安全系數(shù)與敏感系數(shù)隨顆粒容重變化Fig.7 Variations of safety factors and sensitivity factors with volumetric weight of particles

從圖7可以看出，安全系數(shù)與容重基本呈線性關系，當γs從16 kN/m3增加到30 kN/m3時，安全系數(shù)SFs從0.48 上升到1.13，SFr從0.90 上升到3.01，SFs、SFr隨著γs的增加而增大；敏感系數(shù)Sγs從1.88降低到0.76，隨著γs的增加而降低，Sγr恒為1.67。當γs一定時，顆粒更易發(fā)生滑動失穩(wěn)。

土壤顆粒容重的變化，會同時引起抗滑力和滑動力的變化。根據(jù)文獻［21］、［25］可知，邊坡抗滑穩(wěn)定安全系數(shù)一般隨著巖土體容重的增加而減小，與上述得出的結(jié)論相悖。造成這一結(jié)果的主要原因是，松散狀態(tài)下顆粒間的黏聚力c可忽略，即c=0，重力增加抗滑力增加的比例大于滑動力增加的比例，所以得出了抗滑安全系數(shù)隨著容重的增加而增大的結(jié)論。

取γs=25 kN/m3時的敏感系數(shù)為影響因子γs的敏感系數(shù)，即Sγs=1.00、Sγr=1.67，明顯SFr相較于SFs對γs更敏感。

3.5 顆粒等效粒徑di敏感性

令顆粒等效粒徑di從1 mm 增加到25 mm，其余參數(shù)取值如表1所示。根據(jù)不同的di值，求得與di對應的SFs、SFr值，等效粒徑di與安全系數(shù)SF關系如圖8所示。

圖8 安全系數(shù)與敏感系數(shù)隨顆粒等效粒徑變化Fig.8 Variations of safety factors and sensitivity factors with equivalent diameter of particles

從圖8可以看出，安全系數(shù)與等效粒徑基本呈線性關系，當di從1 mm增加到25 mm時，安全系數(shù)SFs從0.48上升到1.04，SFr從0.92 上升到2.61，SFs、SFr隨著di的增加而增大；敏感系數(shù)Sds從0.72 降低到0.33，隨著di的增加而降低，Sdr恒為0.62。當di一定時，顆粒更易發(fā)生滑動失穩(wěn)。

取di=20 mm 時的敏感系數(shù)為影響因子di的敏感系數(shù)，即Sds=0.37、Sdr=0.62。

3.6 土坡孔隙率n敏感性

令土坡孔隙率n從0.25 增加到0.4，其余參數(shù)取值如表1 所示，計算得與n對應的SFs、SFr值。土坡孔隙率n與安全系數(shù)SF關系如圖9所示。

從圖9 可以看出，當n從0.25 增加到0.40 時，安全系數(shù)SFs從1.11 下降到0.90，SFr從2.94 下降到2.02，SFs、SFr隨著n的增加而降低；敏感系數(shù)Sns從0.63 降低到0.44，Snr從1.36 降低到0.66，Sns、Snr隨著n的增加而降低。當n一定時，顆粒更易發(fā)生滑動失穩(wěn)。

圖9 安全系數(shù)與敏感系數(shù)隨土坡孔隙率變化Fig.9 Variations of safety factors and sensitivity factors with porosity of loose soil slope

取n=0.35 時的敏感系數(shù)為影響因子n的敏感系數(shù)，即Sns=0.49、Snr=0.81。

3.7 土坡滲透率K敏感性

令土坡滲透率K從0.03 mm2增加到0.08 mm2，其余參數(shù)取值如表1 所示，計算得與K對應的SFs、SFr值。土坡滲透率K與安全系數(shù)SF關系如圖10所示。

圖10 安全系數(shù)與敏感系數(shù)隨土坡滲透率變化Fig.10 Variations of safety factors and sensitivity factors with permeability of loose soil slope

從圖10 可以看出，當K從0.03 mm2增加到0.08 mm2時，安全系數(shù)SFs從1.15下降到0.48，SFr從3.15下降到0.90，SFs、SFr隨著K的增加而降低；敏感系數(shù)SKs從0.89 降低到0.32，SKr從1.90降低到0.28，SKs、SKr隨著K的增加而降低。當K一定時，顆粒更易發(fā)生滑動失穩(wěn)。

取K=3.93×10-3mm2時的敏感系數(shù)為影響因子K的敏感系數(shù)，即SKs=0.78、SKr=1.39，明顯SFr相較于SFs對K更敏感。

3.8 應力跳躍系數(shù)β敏感性

令應力跳躍系數(shù)β從-0.3 增加到0.7，其余參數(shù)取值如表1所示，計算得與β對應的SFs、SFr值。應力跳躍系數(shù)β與安全系數(shù)SF關系如圖11所示。

圖11 安全系數(shù)與敏感系數(shù)隨應力跳躍系數(shù)變化Fig.11 Variations of safety factors and sensitivity factors with stress jump coefficient

從圖11 可以看出，當β從-0.3 增加到0.7 時，安全系數(shù)SFs從1.22 下降到0.46，SFr從3.53 下降到0.87，SFs、SFr隨著β的增加而降低；敏感系數(shù)Sβs從0.06 升高到0.09，Sβr從0.15 降低到0.08，Sβs隨著β的增加而增大，Sβr隨著β的增加而降低。當β一定時，顆粒更易發(fā)生滑動失穩(wěn)。

取β=0.1 時的敏感系數(shù)為影響因子β的敏感系數(shù)，即Sβs=0.08、Sβr=0.13。

經(jīng)計算，松散層厚度b的敏感系數(shù)小于10-10，松散層厚度b變化對安全系數(shù)SF的影響過小，幾乎可忽略不計，故不再對其進行詳細的敏感性分析；安全系數(shù)影響因子的敏感系數(shù)值如表2所示。

由表2 可知，Sθs＞Sφs＞Sγs≥1＞Shs＞SKs≥0.5＞Sns＞Sds＞Sβs＞Sbs、Sθr＞Sγr＞Shr＞SKr≥1＞Snr＞Sdr≥0.5＞Sβr＞Sbr。土坡坡度θ、土體內(nèi)摩擦角φ、顆粒容重γs對顆粒抗滑安全系數(shù)的影響最敏感，為抗滑安全系數(shù)敏感參數(shù)；坡面徑流水深h、滲透率K影響次之，為次敏感參數(shù)；顆粒等效粒徑di、松散土坡孔隙率n、應力跳躍系數(shù)β、松散層厚度b為不敏感參數(shù)。土坡坡度θ、顆粒容重γs、坡面徑流水深h、滲透率K對顆粒抗?jié)L動安全系數(shù)的影響最敏感，為抗?jié)L動安全系數(shù)敏感參數(shù)；松散土坡孔隙率n、顆粒等效粒徑di、影響次之，為次敏感參數(shù)；應力跳躍系數(shù)β、松散層厚度b為不敏感參數(shù)。

表2 安全系數(shù)影響因子敏感系數(shù)Tab.2 The sensitivity coefficient of safety factor’s impact factors

4 結(jié) 論

（1）對松散土坡淺層非均勻顆粒進行受力分析，分別采用Navier-Stokes方程和擴展Brinkman-Darcy方程描述坡面徑流和坡體滲流運動，并結(jié)合連續(xù)性方程構(gòu)建徑流-滲流聯(lián)合作用數(shù)學模型；視徑流-滲流界面為應力跳躍邊界，推求出流固界面流速解析解和考慮徑流-滲流聯(lián)動作用的單顆粒穩(wěn)定安全系數(shù)。可知，顆粒穩(wěn)定安全系數(shù)主要受到坡面徑流水深、松散層厚度、土坡坡度、應力跳躍系數(shù)、土體內(nèi)摩擦角、顆粒等效粒徑及容重、松散土坡孔隙率及滲透率的影響。

（2）通過敏感性分析，討論了單顆粒滑動、滾動安全系數(shù)對各影響因子的敏感強弱。結(jié)果表明，土坡坡度對顆粒穩(wěn)定性的影響相對最大，敏感系數(shù)分別為1.77 和2.04；抗滑安全系數(shù)的敏感參數(shù)為土坡坡度、土體內(nèi)摩擦角和顆粒容重；抗滑安全系數(shù)的敏感參數(shù)為土坡坡度、顆粒容重、坡面徑流水深及土坡滲透率。□