刘 悦，袁星宇，张 帅，符文熹，李龙国

（1.四川大学水利水电学院，成都610065；2.四川大学水力学与山区河流保护国家重点实验室，成都610065）

0 引 言

松散土坡浅表物质具有结构松散、胶结程度低、抗干扰能力弱、孔隙度大、透水性强等特性，斜坡颗粒极易在降雨、地震、人类活动等因素影响下发生失稳破坏［1，2］。我国西南地区松散土坡分布广泛，降水丰富，发生滑坡和碎屑流灾害严重。因此，有必要开展坡面径流与坡体渗流联动作用下松散土坡颗粒冲刷启动机理研究。

近年来，松散土坡稳定性问题受到国内外学者广泛关注。Rijn 等［3］基于泥沙颗粒的受力特征，通过综合考虑横纵向底坡对泥沙起动的叠加效应，建立了任意坡面上泥沙颗粒处于临界起动状态时的拖曳力公式；谢立全等［4］在进行坡面泥沙受力分析时考虑了渗透力作用，发现向坡外的渗流将大幅降低岸坡泥沙起动条件；Kumar等［5］考虑到岸坡渗流的影响，提出无量纲的水流功率概念，并将其运用到渗流对坡面泥沙起动过程的描述中；Zamani 等［6］通过颗粒离散元法建立了三维土-基础-结构体系，对松散体在不同多向振动作用下的响应进行模拟；田海等［7］进行离心模型试验，对降雨条件下边坡土体强度的变化进行了分析，并得出结论降雨会降低松散堆积体边坡稳定性；毛雪松等［8］基于Midas/GTS 建立松散堆积体路基上边坡稳定数值模型，分析了各坡形因素对边坡稳定性的影响；张岩岩［9］采用一维地表径流运动波方程和二维地下水渗流Richards 方程，将地表径流和地下水渗流耦合，并结合非饱和土Bishop 有效应力及极限平衡理论，研究了降雨作用下滑坡稳定性问题；罗浩等［10］考虑弃渣场坡体结构特征，通过建立边坡颗粒流模型对降雨工况下边坡变形失稳过程进行了模拟；雷小芹等［11］构建描述非饱和堆积土中细颗粒侵蚀-运移-沉积全过程的渗流潜蚀模型，定量研究了细颗粒运移对边坡降雨稳定性的影响。

综合查阅的文献分析表明，有关坡面径流和坡体渗流对松散土坡失稳的联合作用鲜有报道。本文从单颗粒的滑动和滚动角度进行分析，分别采用Navier-Stokes方程［12］和扩展Brinkman-Darcy方程［13，14］描述坡面径流和坡体渗流，联合连续性方程，构建径流-渗流作用下的分析模型；将径流-渗流界面视为应力跳跃边界，推求出流固界面流速表达式，并结合颗粒的滑动和滚动受力分析，得到考虑径流-渗流联动作用下单颗粒稳定安全系数；最后讨论了单颗粒滑动、滚动安全系数对各影响因子的敏感强弱，确定安全系数敏感参数。

1 径流-渗流作用下非均匀颗粒滑动和滚动受力分析

1.1 颗粒滑动分析

松散土坡受外界因素（如降雨）影响，表层颗粒会发生运动。本文假定土坡颗粒为不可压缩非均匀颗粒（等效粒径为di），颗粒只沿土坡倾向方向运动（图1）。

如图1 所示，径流-渗流作用下单个颗粒受到拖曳力FD、上举力FL、水下重力W，以及由于颗粒相对暴露度产生的附加质量力FM；FD、FL、W和FM的计算分别见式（1）~式（4）［15，16］。

图1 松散土坡浅层颗粒受力情况Fig.1 Force of shallow particles on loose soil slope

式中：CD为拖曳力阻力系数；ρ为水的密度，kg/m3；αD为与拖曳力相应的面积系数，取值范围0~1；ub为流固界面流速，m/s。

式中：CL为上举力系数；αL为与上举力相应的面积系数，取值范围0~1；其余符号同前。

式中：αW为与水下重力相应的面积系数，取值范围0~1；γs为泥沙颗粒的容重，N/m3；γw为水的容重，N/m3；其余符号同前。

式中：dmax、dmin、dm分别为颗粒级配中最大、最小和平均粒径，m；αM为与附加质量力相应的面积系数，取值范围为0~1；di为研究颗粒等效粒径，m；dj为研究对象周围颗粒的粒径，m；f（dj）为河床表面混合物以重量计算dj颗粒分布函数；ζ为与颗粒相对暴露度（床面粗化度）有关的系数，颗粒完全暴露时ζ=0，颗粒完全隐蔽时ζ=1。

当土坡颗粒发生滑动时，颗粒还会受到周围其他颗粒的摩阻力Ff，Ff计算见式（5）。

式中：φ为土体内摩擦角；其余符号同前。

将式（3）、式（4）代入式（5）中，Ff可改写为式（6）。

在此，定义颗粒抗滑稳定安全系数SFs为：

式中：Frf为抗滑力，N；Fsf为滑动力，N。

从图1 可见，抗滑力仅为摩擦阻力Ff，即式（6）。滑动力包含水下重力沿x方向的分量Wx=Wsinθ、附加质量力沿x方向的分量FMx=FMsinθ和拖曳力FD。于是，式（7）可具体写出如下：

分析式（8）可知，颗粒抗滑稳定安全系数主要受到流固界面流速ub、土坡坡度θ、土体内摩擦角φ、颗粒等效粒径di、颗粒容重γs的影响。

1.2 颗粒滚动分析

在此，定义颗粒抗滚动稳定安全系数SFr为：

式中：Frm为抗滚动力矩，N∙m；Fsm为滚动力矩，N∙m。

从图2 可见，当土坡颗粒发生滚动时，支持力FN=0，设颗粒绕O点发生滚动，且W、FD、FL、FM的力臂分别为lW、lD、lL、lM。于是，Frm=W lw+FMlM，Fsm=FDlD+FLlL，相应地，式（9）可具体写出如下。

图2 颗粒滚动趋势下受力及力臂情况Fig.2 Force and force arm of rolling particles

分析式（10）可知，颗粒抗滚动稳定安全系数主要受到流固界面流速ub、颗粒等效粒径di、颗粒容重γs的影响。

2 径流-渗流作用土坡流场分析

2.1 模型描述

根据2.1 节、2.2 节颗粒安全系数计算公式可知，径流-渗流联合作用下松散土坡非均匀颗粒是否稳定与流固界面流速ub密切相关；安全系数随着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 流体，为充分发展层流。采用连续性方程［式（11）］和Navier-Stokes 方程［式（12）］描述松散土坡坡面径流运动［12，17］。

式中：ux、uy分别为松散土坡坡面径流沿x、y方向的流速，m/s；fx为沿x方向的质量力，N/kg；P为沿x方向的压强，Pa；υ为水的运动黏滞系数，m2/s；t为时间，s；其余符号同前。

已知水流沿y方向的流速为0，故有∂uy/∂y=0，将其代入式（11）中可得∂ux/∂x=0；x方向流速ux在y方向不发生变化，即∂ux/∂y=0；水流为恒定流，则∂ux/∂t=0；沿x方向fx=gsinθ，∂P/∂x=-ΔP/L。因此，式（12）可写为：

式中：μ为水的动力黏滞系数，Pa∙m，其中μ=υ ρ；L为土坡长度，m。

采用连续性方程［式（14）］和扩展Brinkman-Darcy 方程［式（15）］描述松散土坡坡体渗流运动［13，14，18］。

式中：vx、vy分别为松散土坡坡体渗流沿x、y方向的流速，m/s；n为土颗粒孔隙率；K为土坡渗透率，m2，k=K g/υ，其中k为渗透系数，m/s，g为重力加速度，m/s2。

同理，式（15）可简化为：

2.2 界面流速

为便于计算，现对流速计算公式进行无量纲化处理。引入无量纲参数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为待求系数。

将无量纲参数V=μ v/（ρ g D2tanθ）代入式（19）和式（20）中，即可反求出有量纲流速：

松散土坡坡面径流ux和坡体渗流vx满足以下边界条件［19，20］：

（1）在松散土坡坡面径流上表面（y=b+h）处，流速ux达到最大，即满足∂ux/∂y=0；

（2）在土坡表面（y=b）处，满足速度连续且剪应力跳跃的边界条件，即ux=vx、c2∂vx/∂c1-∂ux/∂c1=β vx/c30.5，其中β为应力跳跃系数，与多孔介质特性、界面位置密切相关，其取值范围为-1～1.5；

（3）在土坡底部（y=0）处，固体颗粒与表面流体不存在相对位移，即满足ux=0。

将上述边界条件代入式（21）、式（22）中，解得待定系数a1、a2、b1、b2，见式（23）～（26）。

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

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

分析式（33）可知，流固界面流速主要受到松散层厚度b、坡面径流水深h、土坡坡度θ、应力跳跃系数β、松散土坡孔隙率n及渗透率K的影响。

3 颗粒稳定性影响因子分析

根据式（8）、式（10）和式（33），松散土坡非均匀颗粒抗滑稳定安全系数SFs的影响因素有松散层厚度b、坡面径流水深h、土坡坡度θ、应力跳跃系数β、土体内摩擦角φ、颗粒等效粒径di、颗粒容重γs、松散土坡孔隙率n、渗透率K；抗滚动稳定安全系数SFr的影响因素有b、h、θ、β、n、K、di及γs。

为定量了解上述因素对松散土坡非均匀颗粒稳定性的影响程度，对各因子进行敏感性分析。敏感系数Si计算见式（34），敏感系数越大，安全系数对该影响因子越敏感［21，22］；敏感分析相关参数基准值如表1所示，其中拖曳力阻力系数CD、上举力系数CL、颗粒所受各力的面积系数α及力臂l取值参考文献［23］、［24］。

表1 相关计算参数基准值［23，24］Tab.1 The base value of related calculation parameters

式中：Si为安全系数影响因子xi的敏感系数；|ΔSFi/SF0|为安全系数SFi的相对变化率；SF0为各参数基准值计算出的安全系数基准值；|Δxi/x0|为影响因子xi的相对变化率；x0为各参数基准值。

3.1 坡面径流水深h敏感性

令坡面径流水深h从0.6 mm 增加到1.8 mm，其余参数取值如表1 所示。根据不同的h值，求得与h对应的SFs、SFr值，径流水深h与安全系数SF关系如图4所示。

图4 安全系数与敏感系数随径流深变化Fig.4 Variations of safety factors and sensitivity factors with runoff depth of slope surface

从图4 可以看出，当h从0.6 mm 增加到1.8 mm 时，安全系数SFs从1.36 下降到0.47，SFr从4.60 下降到0.88，SFs、SFr随着h的增加而降低；敏感系数Shs从1.13 降低到0.46，Shr从3.75 降低到0.39，Shs、Shr随着h的增加而降低。当h一定时，颗粒更易发生滑动失稳。

取h=1.0 mm 时的敏感系数为影响因子h的敏感系数，即Shs=0.90、Shr=1.51，明显SFr相较于SFs对h更敏感。

3.2 坡度θ敏感性

令土坡坡度θ从10°增加到20°，其余参数取值如表1 所示，计算得与θ对应的SFs、SFr值。土坡坡度θ与安全系数SF关系如图5所示。

从图5 可以看出，当θ从10°增加到20°时，安全系数SFs从1.85下降到0.56，SFr从5.11下降到1.26，SFs、SFr随着θ的增加而降低；敏感系数Sθs从3.94降低到0.96，Sθr从5.93降低到0.92，Sθs、Sθr随着θ的增加而降低。当θ一定时，颗粒更易发生滑动失稳。

取θ=15°时的敏感系数为影响因子θ的敏感系数，即Sθs=1.77、Sθr=2.04。

3.3 土体内摩擦角φ敏感性

令颗粒内摩擦角φ从25°增加到40°，其余参数取值如表1所示，计算得与φ对应的SFs、SFr值。内摩擦角φ与安全系数SF关系如图6所示。

图6 安全系数与敏感系数随土体内摩擦角变化Fig.6 Variations of safety factors and sensitivity factors with soil internal friction angle

从图6 可以看出，当φ从25°增加到40°时，安全系数SFs从0.78上升到1.40，SFr恒为2.26，SFs随着φ的增加而增大；敏感系数Sφs从1.11升高到1.52，Sφs随着φ的增加而增大。当φ一定时，颗粒更易发生滑动失稳。

取φ=30°时的敏感系数为影响因子φ的敏感系数，即Sφs=1.21。

3.4 颗粒容重γs敏感性

令颗粒容重γs从16 kN/m3增加到30 kN/m3，其余参数取值如表1 所示，计算得与γs对应的SFs、SFr值。颗粒容重γs与安全系数SF关系如图7所示。

图7 安全系数与敏感系数随颗粒容重变化Fig.7 Variations of safety factors and sensitivity factors with volumetric weight of particles

从图7可以看出，安全系数与容重基本呈线性关系，当γs从16 kN/m3增加到30 kN/m3时，安全系数SFs从0.48 上升到1.13，SFr从0.90 上升到3.01，SFs、SFr随着γs的增加而增大；敏感系数Sγs从1.88降低到0.76，随着γs的增加而降低，Sγr恒为1.67。当γs一定时，颗粒更易发生滑动失稳。

土壤颗粒容重的变化，会同时引起抗滑力和滑动力的变化。根据文献［21］、［25］可知，边坡抗滑稳定安全系数一般随着岩土体容重的增加而减小，与上述得出的结论相悖。造成这一结果的主要原因是，松散状态下颗粒间的黏聚力c可忽略，即c=0，重力增加抗滑力增加的比例大于滑动力增加的比例，所以得出了抗滑安全系数随着容重的增加而增大的结论。

取γs=25 kN/m3时的敏感系数为影响因子γs的敏感系数，即Sγs=1.00、Sγr=1.67，明显SFr相较于SFs对γs更敏感。

3.5 颗粒等效粒径di敏感性

令颗粒等效粒径di从1 mm 增加到25 mm，其余参数取值如表1所示。根据不同的di值，求得与di对应的SFs、SFr值，等效粒径di与安全系数SF关系如图8所示。

图8 安全系数与敏感系数随颗粒等效粒径变化Fig.8 Variations of safety factors and sensitivity factors with equivalent diameter of particles

从图8可以看出，安全系数与等效粒径基本呈线性关系，当di从1 mm增加到25 mm时，安全系数SFs从0.48上升到1.04，SFr从0.92 上升到2.61，SFs、SFr随着di的增加而增大；敏感系数Sds从0.72 降低到0.33，随着di的增加而降低，Sdr恒为0.62。当di一定时，颗粒更易发生滑动失稳。

取di=20 mm 时的敏感系数为影响因子di的敏感系数，即Sds=0.37、Sdr=0.62。

3.6 土坡孔隙率n敏感性

令土坡孔隙率n从0.25 增加到0.4，其余参数取值如表1 所示，计算得与n对应的SFs、SFr值。土坡孔隙率n与安全系数SF关系如图9所示。

从图9 可以看出，当n从0.25 增加到0.40 时，安全系数SFs从1.11 下降到0.90，SFr从2.94 下降到2.02，SFs、SFr随着n的增加而降低；敏感系数Sns从0.63 降低到0.44，Snr从1.36 降低到0.66，Sns、Snr随着n的增加而降低。当n一定时，颗粒更易发生滑动失稳。

图9 安全系数与敏感系数随土坡孔隙率变化Fig.9 Variations of safety factors and sensitivity factors with porosity of loose soil slope

取n=0.35 时的敏感系数为影响因子n的敏感系数，即Sns=0.49、Snr=0.81。

3.7 土坡渗透率K敏感性

令土坡渗透率K从0.03 mm2增加到0.08 mm2，其余参数取值如表1 所示，计算得与K对应的SFs、SFr值。土坡渗透率K与安全系数SF关系如图10所示。

图10 安全系数与敏感系数随土坡渗透率变化Fig.10 Variations of safety factors and sensitivity factors with permeability of loose soil slope

从图10 可以看出，当K从0.03 mm2增加到0.08 mm2时，安全系数SFs从1.15下降到0.48，SFr从3.15下降到0.90，SFs、SFr随着K的增加而降低；敏感系数SKs从0.89 降低到0.32，SKr从1.90降低到0.28，SKs、SKr随着K的增加而降低。当K一定时，颗粒更易发生滑动失稳。

取K=3.93×10-3mm2时的敏感系数为影响因子K的敏感系数，即SKs=0.78、SKr=1.39，明显SFr相较于SFs对K更敏感。

3.8 应力跳跃系数β敏感性

令应力跳跃系数β从-0.3 增加到0.7，其余参数取值如表1所示，计算得与β对应的SFs、SFr值。应力跳跃系数β与安全系数SF关系如图11所示。

图11 安全系数与敏感系数随应力跳跃系数变化Fig.11 Variations of safety factors and sensitivity factors with stress jump coefficient

从图11 可以看出，当β从-0.3 增加到0.7 时，安全系数SFs从1.22 下降到0.46，SFr从3.53 下降到0.87，SFs、SFr随着β的增加而降低；敏感系数Sβs从0.06 升高到0.09，Sβr从0.15 降低到0.08，Sβs随着β的增加而增大，Sβr随着β的增加而降低。当β一定时，颗粒更易发生滑动失稳。

取β=0.1 时的敏感系数为影响因子β的敏感系数，即Sβs=0.08、Sβr=0.13。

经计算，松散层厚度b的敏感系数小于10-10，松散层厚度b变化对安全系数SF的影响过小，几乎可忽略不计，故不再对其进行详细的敏感性分析；安全系数影响因子的敏感系数值如表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。土坡坡度θ、土体内摩擦角φ、颗粒容重γs对颗粒抗滑安全系数的影响最敏感，为抗滑安全系数敏感参数；坡面径流水深h、渗透率K影响次之，为次敏感参数；颗粒等效粒径di、松散土坡孔隙率n、应力跳跃系数β、松散层厚度b为不敏感参数。土坡坡度θ、颗粒容重γs、坡面径流水深h、渗透率K对颗粒抗滚动安全系数的影响最敏感，为抗滚动安全系数敏感参数；松散土坡孔隙率n、颗粒等效粒径di、影响次之，为次敏感参数；应力跳跃系数β、松散层厚度b为不敏感参数。

表2 安全系数影响因子敏感系数Tab.2 The sensitivity coefficient of safety factor’s impact factors

4 结 论

（1）对松散土坡浅层非均匀颗粒进行受力分析，分别采用Navier-Stokes方程和扩展Brinkman-Darcy方程描述坡面径流和坡体渗流运动，并结合连续性方程构建径流-渗流联合作用数学模型；视径流-渗流界面为应力跳跃边界，推求出流固界面流速解析解和考虑径流-渗流联动作用的单颗粒稳定安全系数。可知，颗粒稳定安全系数主要受到坡面径流水深、松散层厚度、土坡坡度、应力跳跃系数、土体内摩擦角、颗粒等效粒径及容重、松散土坡孔隙率及渗透率的影响。

（2）通过敏感性分析，讨论了单颗粒滑动、滚动安全系数对各影响因子的敏感强弱。结果表明，土坡坡度对颗粒稳定性的影响相对最大，敏感系数分别为1.77 和2.04；抗滑安全系数的敏感参数为土坡坡度、土体内摩擦角和颗粒容重；抗滑安全系数的敏感参数为土坡坡度、颗粒容重、坡面径流水深及土坡渗透率。□