1. 引言
在边坡治理中,为了维护自然景观,常常通过构建植物生长所需土壤环境和种植植物来对裸露岩质边坡进行生态修复。但由于护坡土体与岩石的弹性模型相差过大[1] [2]难以形成一个统一整体,因此护坡土体与岩石之间的接触面成为了薄弱滑面。降雨作用下,土体易沿着土体与岩石接触面发生滑动失稳,对安全形成严重的潜在威胁。因此研究暴雨作用下岩质边坡生态护坡体稳定性具有重要的意义。
岩质边坡植物护坡土体稳定性由护坡土体本身的力学特性和护坡土体与基岩接触面剪切力学特性共同决定。汪优、DU P、LI D J分别从接触面粗糙度[3]、剪切方式[4]、土体类型[5]等方面进行了实验研究,认为土体–岩石接触面剪应力–剪切位移曲线呈现出应变硬化型和应变软化型两种不同类型。界面抗剪强度大于土体抗剪强度时,接触面剪切变形以界面附近土体剪切变形为主[3],剪应力–剪切位移曲线表现出应变软化特性。植物还能改变土体的水力特性,研究表明,植物根系能影响土体的渗透系数[6]、基质吸力[7]、保水性[8]、物理化学特性[8]和剪胀性[9]。植物根系对土体渗透系数的影响存在两种不同的观点:一部分学者研究后认为根系生长过程中对土体的挤密作用会降低土体的孔隙率,土体的渗透系数也随之降低[10];另一部分学者认为植物根系的生长、腐烂会在土体内部形成优势渗流通道,增大土体的渗透系数[11]。然而,在暴雨工况下,综合考虑根土复合体的力学和水力特性的护坡土体稳定性的研究仍较为匮乏。
2. 模型选取、开发与适用性验证
2.1. 护坡土体本构模型选取
护坡土体本质上为根土复合体,选用Duncan-Chang模型作为根土复合体的本构模型。Duncan-Chang模型为非线性弹性本构模型,该模型采用切向弹性模量和切向泊松比作为土体的弹性模量和泊松比[12]。土体的切线弹性模量和切线泊松比的计算公式分别如式(1)和式(2)所示。
(1)
(2)
Et——切线弹性模量(kPa);
vt——土体的切线泊松比;
——最大主应力(kPa);
——最小主应力(kPa);
pa——大气压(kPa);
Rf——破坏比;
c、
——土体的粘聚力(kPa)和内摩擦角(˚);
2.2. 排版规范的完整性
护坡土体–基岩接触面相互作用力学模型选取
选取Guo等[13]提出的接触面剪切损伤模型作为根土复合体与风化基岩接触面的相互作用力学模型。该模型为从曲线拟合角度出发建立的模型,参数易得,模拟结果较为准确,且能同时模拟接触面剪切过程中的剪切硬化和剪切软化特征。剪切损伤模型的表达式如式(3)所示。
(3)
式中:
——作用在剪切面上的剪应力(kPa);
s——接触面上的剪切位移(mm);
G——接触面上剪切刚度(kPa/mm);
——残余强度(kPa);
——屈服强度(kPa);
F0、m——待求解的变量。
2.3. 模型开发与适用性验证
利用ABAQUS中预留的UMAT子程序接口对Duncan-Chang模型进行开发,同时利用FRIC子程序接口对接触面剪切损伤模型进行开发。采用Fortran语言进行编程,定义根土复合体的本构方程以及土岩接触面的接触方程,将其作为子程序调用到ABAQUS中。子程序开发完毕后,并进行了三轴剪切实验,通过将数值模拟结果与实验数据对比来进行验证(见图1),结果显示平均误差只有4.9%,因此认为使用所选用的Duncan-Chang模型和接触面剪切损伤模型是合适的。
Figure 1. Comparison of numerical simulation results with experimental results
图1. 数值模拟结果与试验结果对比
3. 护坡土体稳定性数值模拟研究
3.1. 研究工况设置
护坡土体稳定性受多种因素的影响,本文选取边坡坡率、植物类型和基岩风化程度作为研究因素。边坡坡率选取1:0.75、1:1.0和1:1.25。对于植物类型,设置无植物、狗牙根草和香根草三种研究工况;对于岩石风化类型,设置未风化、弱风化和强风化三种研究工况。
3.2. 模型建立与参数选取
方型框架护坡是岩质边坡锚杆框架护坡工程中常用的框架类型,如图2所示。模型建立时选取单个框架单元内的护坡土体进行分析,护坡土体尺寸为3 m × 3 m × 0.1 m,数值模拟模型如图3所示。
Figure 2. Framework structure ecological slope protection
图2. 框架植被护坡
Figure 3. Numerical simulation model
图3. 数值模拟模型
岩质边坡护坡土体稳定性数值模拟研究中的岩体部分采用线性弹性本构模型,见表1。护坡土体的本构模型采用第2节中开发的Duncan-Chang模型,见表2。根土复合体–基岩接触面相互作用力学模型采用剪切损伤模型,见表3。
Table 1. Rock mass model parameters [14]
表1. 岩体模型参数[14]
风化程度 |
抗压强度(MPa) |
切线弹性模量(GPa) |
切线泊松比 |
未风化 |
150 |
45 |
0.20 |
弱风化 |
75 |
24 |
0.31 |
强风化 |
16 |
9 |
0.40 |
Table 2. Soil model parameters [15] [16]
表2. 土体模型参数[15] [16]
土体类型 |
K |
n |
Rf |
c (kPa) |
Φ (˚) |
F |
G |
D |
无植物 |
59.977 |
1.759 |
0.8547 |
12.295 |
17.54 |
0.0859 |
0.2643 |
2.5037 |
狗牙根 |
79.864 |
2.043 |
0.8733 |
15.64 |
18.58 |
0.0899 |
0.2642 |
2.4939 |
香根草 |
97.731 |
1.983 |
0.8780 |
18.024 |
19.60 |
0.0871 |
0.2120 |
3.1500 |
Table 3. Root-soil composite-bedrock contact surface model parameters [17]
表3. 根土复合体–基岩接触面的模型参数[17]
计算工况 |
E (MPa) |
|
c (kPa) |
φ (˚) |
无植物 |
10 |
0.3 |
9.03 |
18.17 |
狗牙根草 |
11.23 |
19.25 |
香根草 |
19.73 |
26.65 |
3.3. 数值模拟结果
3.3.1. 开挖坡率对护坡土体的影响
边坡坡率分别为1:0.75、1:1.0和1:1.25时,未风化基岩 + 无植物工况的护坡土体位移云图如图4,底部框架的位移限制使得护坡土体位移值分布呈现出从坡底土体到坡顶土体逐渐增大的趋势。当边坡坡率从1:0.75分别降低至1:1.0和1:1.25时,未风化基岩 + 无植物工况的护坡土体最大位移从8.11 mm分别降低至6.95 mm和6.08 mm,下降比例分别为14.3%和25%,表明开挖坡率的降低能有效地减少边坡位移。
(a) (b)
(c)
Figure 4. Displacement cloud map of unweathered bedrock and non-vegetated slope soil, (a) slope ratio 1:0.75 condition, (b) slope ratio 1:1.0 condition, (c) slope ratio 1:1.25 condition
图4. 未风化基岩 + 无植物工况的护坡土体位移云图,(a) 坡率1:0.75工况,(b) 坡率1:1.0工况,(c) 坡率1:1.25工况
3.3.2. 植物类型对护坡土体的影响
以未风化基岩 + 坡率1:1.0的工况为例,无植物、种植狗牙根草和种植香根草的护坡土体位移云图如图5。当护坡土体无植物种植时,最大位移为6.95 mm;当种植狗牙根草时,最大位移为2.75 mm,相比无植物工况,位移下降了60.4%;当种植香根草时,最大位移为2.11 mm,相比无植物工况,位移下降了69.6%。表明植物根系的加筋作用可以大幅提高护坡土体的抗剪强度,显著减低变形。
(a) (b)
(c)
Figure 5. Displacement cloud map of unweathered bedrock and slope ratio 1:1.0 condition slope soil, (a) non-vegetated, (b) planted with Bermuda grass, (c) planted with vetiver grass
图5. 未风化基岩 + 坡率1:1.0工况的护坡土体位移云图,(a) 无植物,(b) 种植狗牙根草,(c) 种植香根草
3.3.3. 基岩风化程度对护坡土体的影响
以种植香根草 + 坡率1:1.0的工况为例,护坡土体位移云图见图6。当基岩未风化时护坡土体的最大位移为2.11 mm,弱风化时最大位移为2.03 mm,强风化时最大位移为1.96 mm。风化程度越高反而位移越小,原因为基岩风化之后会有很多裂隙,种植植物后植物根系伸展进入裂隙中,增强了护坡土体与基岩接触面的剪切强度,进而降低变形。
(a) (b)
(c)
Figure 6. Displacement cloud map of slope soil planted with vetiver grass and slope ratio 1:1.0 condition, (a) unweathered bedrock, (b) weakly weathered bedrock, (c) strongly weathered bedrock
图6. 种植香根草 + 坡率1:1.0工况的护坡土体位移云图,(a) 未风化基岩,(b) 弱风化基岩,(c) 强风化基岩
4. 降雨作用下护坡土体稳定性数值模拟研究
4.1. 研究工况设置
植物类型改变土体抗剪强度、根土复合体–风化基岩接触面剪切力学特性和根土复合体水力特性,从而影响降雨作用下护坡土体的稳定性。因此,选取植物类型和降雨强度两个研究因素,建立数值模拟模型分析降雨作用下护坡土体稳定性,研究工况与研究水平如表4所示。
Table 4. Research factors and research levels
表4. 研究因素与研究水平
研究因素 |
研究水平 |
植物类型 |
无植物 |
狗牙根草 |
香根草 |
降雨强度 |
暴雨 |
特大暴雨 |
4.2. 模型建立与参数选取
4.2.1. 模型建立
数值模拟模型与图3相同,选取边坡坡率为1:1.00,对单个框架内的护坡土体进行建模分析。框架对护坡土体的影响简化于模型边界条件中,同时将护坡土体边界设置为不透水边界。模型共设置两个分析步:在第一个分析步中对模型施加重力荷载,进行地应力平衡;在第二个分析步中施加降雨边界条件,降雨时长为24 h。
4.2.2. 参数选取
1) 土体饱和度–基质吸力关系
程鹏等[17]利用现场监测方法获得了无植物、种植狗牙根草工况和种植香根草工况的土体水土特征曲线及非饱和渗透系数,计算参数见表5。
Table 5. Calculation parameters
表5. 计算参数
计算工况 |
E (MPa) |
|
c (kPa) |
φ (˚) |
饱和渗透系数
|
无植物 |
10 |
0.3 |
9.03 |
18.17 |
5.00 × 10−6 |
狗牙根草 |
11.23 |
19.25 |
6.77 × 10−7 |
香根草 |
19.73 |
26.65 |
1.02 × 10−5 |
2) 岩体渗透特性
胡云进[18]将裂隙岩体等效为连续介质体,建立了其饱和非饱和渗流数学模型。本模型中将岩体等效为连续介质体,风化岩体渗透参数参考该文献中的风化基岩试验数据,其渗透系数为1.67 × 10−7 m/s。
3) 降雨强度
暴雨工况的降雨量为75 mm/24 h,特大暴雨工况的降雨量取深圳市2022年“5.12”特大暴雨事件的24小时降雨量,为267.1 mm。降雨雨型为单峰型,降雨时长为24 h。两种工况的小时降雨量分别如图7所示。
(a)
(b)
Figure 7. Rainfall conditions, (a) hourly rainfall (rainstorm condition), (b) hourly rainfall (extreme rainstorm condition)
图7. 降雨工况,(a) 小时降雨量(暴雨工况),(b) 小时降雨量(特大暴雨工况)
4.3. 数值模拟结果
4.3.1. 无植物工况计算结果
施加24 h暴雨工况后,护坡土体位移云图如图8(a)所示。降雨作用下,近底部框架处护坡土体位移出现较大增长,土体的最大位移值为5.4 mm。中部和顶部土体的位移变化量较小。与施加暴雨前(护坡土体的最大位移绝对值为4.389 mm)相比,护坡土体的最大位移绝对值增加1.02 mm,增幅约为23.24%。施加24 h降雨后,护坡土体位移矢量图如图8(b)所示。与初始状态相比,底部土体位移向临空侧方向发展,中部和顶部土体位移方向仍为平行坡面向下。施加24 h降雨后,无植物工况土体内部塑性区如图8(c)所示。降雨初期,底部表层土体最先出现塑性区,随降雨时长增大,塑性区范围逐渐扩大,施加24 h降雨后土体塑性区扩展至土–岩接触面。
(a) (b)
(c)
Figure 8. Calculation results of slope soil under non-vegetated condition in rainstorm condition (t = 24 h; unit: m), (a) displacement cloud map, (b) displacement vector map, (c) plastic zone cloud map
图8. 暴雨工况下无植物工况护坡土体计算结果(t = 24 h;单位:m),(a) 位移云图,(b) 位移矢量图,(c) 土体塑性区云图
特大暴雨工况下,未种植植物时,护坡土体在t = 7.57 h时计算不收敛,护坡土体发生变形破坏。此时,护坡土体位移云图如图9(a)所示。特大暴雨作用下,框架底部土体位移迅速发展,土体的最大位移为17.64 mm,护坡土体的最大位移绝对值是施加降雨前的4.02倍。变形破坏时的护坡土体位移矢量图和塑性区云图分别如图9(b)、图9(c)所示。受特大暴雨影响,底部表层土体出现塑性区,且随降雨时长增大,护坡土体塑性区逐渐扩展并发展至土–岩接触面。土体内部塑性区与土–岩接触面共同形成薄弱滑带,塑性区内部的土体被侧向挤出,土体位移方向指向临空侧。同时,底部土体侧向挤出导致中部和上部土体失去竖向支撑,中部和上部土体岩土–岩接触面向下滑动。
(a) (b)
(c)
Figure 9. Calculation results of slope soil under non-vegetated condition in extreme rainstorm condition (t = 24 h; unit: m), (a) displacement cloud map, (b) displacement vector map, (c) plastic zone cloud map
图9. 特大暴雨工况下无植物工况护坡土体计算结果(t = 24 h;单位:m),(a) 位移云图,(b) 位移矢量图,(c) 土体塑性区云图
4.3.2. 狗牙根草工况计算结果
施加24 h暴雨工况后,护坡土体位移云图和位移矢量图如图10(a)和图10(b)所示。降雨作用下,狗牙根草工况护坡土体位移分布与初始状态相似,呈现出沿坡底到坡顶土体位移逐渐增大的规律。另外,降雨入渗导致护坡土体内的负孔隙水压力压力迅速减小,土体有效应力减小,因此土体位移较初始状态减小。施加降雨后护坡土体的最大位移绝对值为1.908 mm,与初始状态相比,土体的最大位移绝对值减小1.149 mm,降幅为37.59%。狗牙根草工况护坡土体塑性应变云图如图10(c)所示。降雨作用下,底部表层土体内出现小面积塑性区,塑性区并未贯通至根土复合体–风化基岩接触面,塑性区面积较无植物相比明显变小。
(a) (b)
(c)
Figure 10. Calculation results of slope soil under Bermuda grass condition in rainstorm condition (t = 24 h; unit: m), (a) displacement cloud map, (b) displacement vector map, (c) plastic zone cloud map
图10. 暴雨工况下种植狗牙根草工况护坡土体计算结果(t = 24 h;单位:m),(a) 位移云图,(b) 位移矢量图,(c) 土体塑性区云图
特大暴雨工况下,种植狗牙根草时护坡土体在t = 7.73 h时计算不收敛,护坡土体发生变形破坏。此时,护坡土体位移云图如图11(a)所示。特大暴雨作用下,框架底部表层土体位移迅速发展,土体的最大位移为10.97 mm,土体的最大位移绝对值是施加降雨前的3.59倍。种植有狗牙根草的护坡土体变形破坏时的土体的位移矢量图与塑性区云图分别如图11(b)和图11(c)所示。与无植物工况类型,特大降雨作用下底部表层土体出现塑性区,且塑性区随降雨时长的增大逐渐扩展至土–岩接触面。但是,种植有狗牙根草的护坡土体的塑性区范围较未种植植物时小,尤其是接触面处的塑性区的范围远小于无植物工况。由于塑性区的出现与发展,底部土体位移向临空侧方向发展,上部土体沿着接触面向下滑动。
(a) (b)
(c)
Figure 11. Calculation results of slope soil under Bermuda grass condition in extreme rainstorm condition (t = 24 h; unit: m), (a) displacement cloud map, (b) displacement vector map, (c) plastic zone cloud map
图11. 特大暴雨工况下种植狗牙根草工况护坡土体计算结果(t = 24 h;单位:m),(a) 位移云图,(b) 位移矢量图,(c) 土体塑性区云图
4.3.3. 香根草草工况计算结果
施加24 h暴雨工况后,护坡土体位移云图如图12(a)所示。降雨作用下,底部1 m范围内土体的位移增长显著,土体的最大位移出现在距底部框架0.7 m左右的区域,最大位移值为8.88 mm。同时,降雨作用下顶部土体的负孔隙水压力迅速减小,土体有效应力降低,因此顶部土体的位移较施加降雨前减小。与施加降雨前的初始状态相比,护坡土体的最大位移绝对值增长了6.265 mm,增幅约为239.76%。护坡土体位移矢量图如图12(b)所示。距底部框架1 m范围内的土体的位移方向发生偏转,由平行坡面向下偏转为近乎垂直坡面方向。中部和顶部土体的位移方向与初始状态一致,为平行坡面向下。护坡土体塑性应变云图如图12(c)所示。虽然施加24 h降雨后,底部土体的位移出现显著增长,位移方向发生偏转,但是土体内部并未出现贯通至根土复合体–基岩接触面的塑性区,仅坡脚表层极小范围内的土体内部出现塑性区。
(a) (b)
(c)
Figure 12. Calculation results of slope soil under vetiver grass condition in rainstorm condition (t = 24 h; Unit: m), (a) displacement cloud map, (b) displacement vector map, (c) plastic zone cloud map
图12. 暴雨工况下种植香根草工况护坡土体计算结果(t = 24 h;单位:m),(a) 位移云图,(b) 位移矢量图,(c) 土体塑性区云图
特大暴雨工况下,种植香根草时护坡土体在t = 7.87 h时计算不收敛,护坡土体发生变形破坏。此时,护坡土体位移云图如图13(a)所示。特大暴雨作用下,框架底部表层土体位移迅速发展,土体的最大位移为20.14 mm,护坡土体的最大位移量绝对值是施加降雨前的7.71倍。种植有香根草的护坡土体发生变形破坏时的位移矢量图如图13(b)。与施加降雨前相比,底部2 m范围内的土体的位移方向均发生偏转,土体位移方向由初始的平行坡面向下偏转为垂直坡面方向。种植有香根草的护坡土体发生变形破坏时的塑性区云图如图13(c)所示。与无植物工况和狗牙根草工况不同,种植有香根草时的土–岩接触面抗剪强度大幅提高,护坡土体不易沿着接触面发生滑动。随降雨时长增加,近底部框架的护坡土体表层出现圆弧形的贯通塑性区,护坡土体变形破坏方式以表层土体局部滑动为主。
(a) (b)
(c)
Figure 13. Calculation results of slope soil under vetiver grass condition in extreme rainstorm condition (t = 24 h; unit: m), (a) displacement cloud map, (b) displacement vector map, (c) plastic zone cloud map
图13. 特大暴雨工况下种植香根草工况护坡土体计算结果(t = 24 h;单位:m),(a) 位移云图,(b) 位移矢量图,(c) 土体塑性区云图
4.4. 护坡土体变形破坏机理分析
降雨作用下,不同植物类型工况中护坡土体向临空侧方向的位移量出现显著增长的时间与土体中出现正孔隙水压力的时间呈现出高度相关性。简而言之,土体中正孔隙水压力的出现与发展是降雨作用下护坡土体发生变形破坏的主要因素。降雨作用下,土体中出现正孔隙水压力后,孔隙水压力减弱了土颗粒之间的咬合作用,降低了相邻土颗粒接触面上的正应力,使得土体的抗剪强度降低,土体中更容易出现变形破坏。
与无植物工况相比,狗牙根草根系沿水平生长导致土体的渗透系数减小,降雨无法及时、快速入渗,表层土体孔隙水压力迅速增大但深层土体的孔隙水压力增长较小,特大暴雨作用下土体变形较无植物工况小。种植香根草时,香根草根系在土体中易形成优势渗流通道,土体的渗透系数较无植物工况大。虽然香根草根系增大了土体与接触面的抗剪强度,但降雨作用下土体孔隙水压力增长迅速,土体变形较无植物工况大。
5. 结论
以岩质边坡植物护坡土体稳定性为研究对象,使用ABAQUS数值模拟方法分析了岩质边坡植物护坡土体变形破坏机理,研究了植物类型、土体类型、基岩风化程度、边坡开挖坡率、降雨作用等因素对岩质边坡护坡土体稳定性的影响。本文取得了以下研究成果:
1) 开发了非线性弹性本构模型与非线性接触面力学模型。使用接触面剪切试验数据对非线性接触面力学模型进行验证,平均误差仅为4.9%。开发的本构模型和接触面力学模型实现了对含根系土体力学特性的模拟,为岩质边坡植物护坡土体稳定性分析提供了研究手段。
2) 使用开发的本构模型和接触面模型建模分析了边坡开挖坡率、植物类型和基岩风化程度对岩质边坡植物护坡土体稳定性的影响。开挖坡率增大,边坡最大位移也随之增大。不考虑降雨时,种植植物后均能够显著减低土体的滑动变形,变形量随基岩风化程度的增加而减小。考虑降雨时,以须根系为主的狗牙根草能降低变形量,以直根系为主的香根草则会增大变形量。
3) 对于所在地区降雨量少、岩石风化程度较大的岩质边坡,建议选用以直根系为主的香根草作为护坡植物,直根系加筋作用明显,有利于岩质边坡护坡土体的稳定性。边坡所在区域降雨量较大时,宜采用以须根系为主的狗牙根草作为护坡植物,须根系能降低土体渗透系数,减少降雨入渗,提高护坡土体的稳定性。正孔隙水压力是导致岩质边坡植物护坡土体变形破坏的主要因素。因此,在降雨量较大的地区进行岩质边坡植被修护时,可采用在护坡土体表面增设有效的排水措施来增大排水能力,降低土体中孔隙水压力的增长速率,提高护坡土体的稳定性。
基金项目
中国铁路广州局集团有限公司科技研究开发计划(2022K057-Z);国家重点研发计划项目(2021YFC3001002)。
NOTES
*通讯作者。