河谷下切对岸坡卸荷影响的数值模拟研究

李万才，邓辉，苏航，丁中辉，陆泌锋，郝浩，李安润

(1.成都理工大学地质灾害防治与地质环境保护国家重点实验室，四川成都610059； 2.四川省地质工程勘察院集团有限公司，四川成都610072；3.中交隧道局第四工程有限公司，四川成都610031)

经由漫长的地质演化，河谷地貌最终得以形成。在该过程中，强烈的侵蚀和剥蚀作用促使了深切河谷的形成。不少学者在深切河谷地应力场方面做了许多研究：黄润秋[1]在分析总结西南地区大量工程实践的基础上，建立了卸荷条件下岩石高边坡发育的动力过程及三阶段演化模式；梁瑶[2]围绕锦屏水电站所处的特殊地质地貌特征，对河谷地区地应力场和边坡开挖稳定性进行了深入的研究，确定河谷的发育演化是影响河谷地区地应力场成因的重要因素并详细讨论了地应力的松弛卸荷对边坡稳定性的影响；刘亚群等[3]通过对南水北调西线工程多个坝址区地应力测试的结果进行统计分析，初步得出深切河谷区最大水平主应力量值随深度呈分段线性关系的结论；端木杰超[4]通过建立深切河谷应力场的弹塑性有限元模型，模拟了侵蚀基准面及表层风化对河谷应力场的影响以及不同构造应力模式对河谷应力场的影响；黄书岭等[5]通过将复杂地质环境下三维地应力场模型与洞室群精细开挖计算模型相耦合，提出考虑河谷演化规律的厂址区域地应力场量化精细数值模型的建立方法和分析思路；李超等[6]通过数值分析研究了不同坡角的山坡在重力场和构造应力条件下山体中在垂直向和水平向典型剖面中初始地应力分布的特点；李永松等[7]提出了深切河谷应力场非线性系统分析方法并实际运用到乌东德水电站坝址区岩体应力场的研究过程中；胡斌等[8]基于实测数据，利用统计分析手段分析了某深切河谷和岸坡中的地应力场的特征，得到了河谷区和岸坡中最大水平应力的方向，并分析了河谷区最大水平应力量值与埋深的关系。

相较平坦的地形而言，由于受到强烈的剥蚀、侵蚀作用以及沉积作用等因素的影响，河谷地区地应力状态发生一定的变化，从而在局部形成特殊的应力场。深切河谷地区的地形高低不平，呈现剧烈起伏的形态，可将内部应力场看作为随河谷不断演化而造成的原始应力场不断变化的局部应力场。由于后期受到地形起伏的影响、地表剥蚀作用、河流的侵蚀作用以及风化等作用的综合影响，先期形成地应力场不断调整并发生改变。从河谷的发育机理来看，河谷地形的演变可近似看作为“岩体的开挖”，这种开挖是指对应于河谷的剥蚀下切(不考虑河谷水文的变化，同时不考虑岸坡孔隙水压的变化)。

1 下坝址区地质条件

下坝址位于色汝小河下游的位置。左岸岸坡坡度大多在35～50°之间，局部位置坡度在60°以上，右岸岸坡在3 078 m高程以下的岸坡坡度稍微陡峭，大多在50～60°之间，局部位置在60～70°之间，而在3 078 m高程以上的斜坡坡度多在30～40°之间，河谷总体呈“V”形，河谷形貌影响着河谷地应力分布[9-11]。坝址一带左岸受色汝河切割影响，形成高陡山脊，坝区右岸临江边坡高度多在2 000 m之上。坝址区岩性主要为英安岩，左岸出现多旋回含角砾英安岩和凝灰岩，坝址左岸至河床位置分布有印支期(γδ)的花岗闪长斑岩。河床部位为冲积砂、卵砾石夹近缘崩塌堆积碎块石。经坝址区地质测绘，坝址区的主要构造有韧性长大裂隙、剪切带、小断层、平硐揭露裂隙及柱状节理等地质构造。其中，坝址区韧性剪切带的宽度约30～150 m，剪切带内具有明显的岩体片理化以及糜棱化现象。坝址区出露地层主要为英安岩，物理地质现象主要有风化堆积体、碎屑流、泥石流、卸荷、崩塌、蚀变、碎裂松动岩体等不良地质现象。本区主要发育有怒江、澜沧江、金沙江及其他支流，沿两岸常发育有对称或不对称的河谷阶地，阶地在峡谷地段不发育，宽谷段发育有五级阶地[12-13]，见表1、图1。根据相关资料显示，自进入峡谷期以来，澜沧江河谷下切速率很快，下切速率在1.9 ～4.4 mm/a之间。

表1 西藏曲孜卡澜沧江河谷阶地特征

1.河流相砂砾石层；2.砂岩；3.泥页岩；4.阶地编号(形成时代)图1 西藏曲孜卡澜沧江河流阶地实测剖面

2 数值模拟建立

2.1 模型的建立

建模建立前期，采用CAD和Surfer软件进行高程数据导出，然后导入到ANSYS中进行模型的搭建和网格划分，最终再生成FLAC 3D可读入模型进行计算。针对五次河谷的下切，在模型计算过程中采用“null”命令依次对各个阶地赋为空模型来模拟河谷的下切，现今河谷按照现今实际情况进行模拟。计算模型按弹塑性材料考虑，破坏准则采用莫尔-库伦强度准则。

根据下坝工程平面地质图、建模区工程地质条件并结合河谷下切深度表(河谷下切深度取平均值)建立地质模型，见图2。模型X方向(与河道近于垂直)长750 m，Y方向(与河道近于平行)长200 m，Z方向(垂直方向)高程界于2 800～3 240 m。在模型建立过程中，河谷现今状态建立时考虑了强弱卸荷带和韧性剪切带，在恢复河谷阶地时未考虑韧性剪切带。经网格化之后，整个模型共划分为21 467个节点和116 450个单元，模型分为9个体，其中体1代表第I级阶地，体2代表第II级阶地，体3代表第III级阶地，体4代表第IV级阶地，体5代表第V级阶地，体6代表强卸荷岩体，体7代表弱卸荷岩体，体8代表基岩，体9代表韧性剪切带岩体。

图2 网格化后模型

2.2 参数选取及边界条件的确定

模型区内一共有4种材料，在模型建立过程中，考虑到网格划分，在韧性剪切带内并未划分卸荷带，但对其岩体参数进行了一定的折减。模型采用的力学参数取值参照贵阳勘测设计研究院提供的相关力学实验结果，考虑模型简化后进行参数取值，最终取值见表2。模型计算所需的体积模量、剪切模量可由下列公式计算得到：

(1)

(2)

式中K——体积模量；G——剪切模量；E——弹性模量；μ——泊松比。

表2 岩体力学参数取值

考虑到河谷下切以及岩体结构面的影响，在河谷5次下切过程中，韧性剪切带按照给出建议值进行取值，其他块体的岩体物理力学参数均取未卸荷带参数进行一定的折减再带入进行计算；在计算现今状态时，按照现今实际情况的岩体参数进行选取。

模型4个侧面以及底面采用位移约束条件，模型顶面采用自由边界，模型为应力应变状态。一般说来，自然斜坡卸荷作用随侧向地应力增大而增大，而在构造地应力不大的地区，侧向地应力主要由岩体自重应力产生；在深切河谷下切过程中，自重应力变化是导致斜坡发生卸荷变形的主要原因[14-16]。因此，在本次数值模拟过程中，仅考虑自重应力，未考虑水平应力的影响，自重应力则通过在Z方向上施加重力加速度来实现。

3 应力场特征

为了更好地反映河谷下切过程中应力变化情况，选取Y=100 m处的纵剖面为对象，下面将从该剖面的变化情况来分析应力场的变化情况。在模拟河谷下切时，将对应的阶地设为空模型来模拟阶地的下切。

河谷下切导致了主应力方向的变化。坡体深部最大主应力呈近竖直向且垂直于最小主应力，当接近临空面时，其方向发生变化，近似与坡面平行。坡体深部最小主应力呈近水平向，当接近临空面时，其方向也发生变化，近似与坡面垂直。

剖面最大主应力云见图3，剖面最小主应力云见图4。由图3可知，最大主应力云有以下特征。

a) 随着河谷的持续下切，应力发生重新分布，坡体最大主应力逐渐降低；在坡体顶部、坡表和韧性剪切带附近最大主应力值较小，最大主应力最大值分布区域近似沿河道呈现为对称形式。

b) 在同一高程位置，由坡表向坡内延伸，最大主应力逐渐增大。

c) 最大主应力主要为压应力，在越靠近坡表的位置，应力水平就越低。

a) 第一次下切

b) 第二次下切

c) 第三次下切

d)第四次下切

e)第五次下切

f)现今河谷图3 最大主应力演化过程

a) 第一次下切图4 最小主应力演化过程

b) 第二次下切

c) 第三次下切

d)第四次下切

e)第五次下切

f)现今河谷续图4 最小主应力演化过程

由图4可知，最小主应力云特征如下：① 随着河谷持续下切，最小主应力逐渐降低，应力集中位置由河谷位置逐渐向河道两侧转移;②最小主应力以压应力为主，越靠近坡表，应力水平就越低,在坡表及坡顶位置(尤其是在坡顶和韧性剪切带附近)分布有少量的拉应力区。

总体来看，河谷下切造成了最大、最小主应力值的降低，应力随着河谷下切发生重分布，应力方向也产生调整；应力最小值主要出现在坡顶和韧性剪切带附近，韧性剪切带对斜坡应力影响较大；整体来看，在现今状态下，斜坡坡脚位置有一定的应力集中现象；在同一高程下，坡表应力比坡体内部应力要低， 河谷下切后，坡体内部应力逐渐调整，坡体整体应力都出现一定的降低，在坡表部分位置(比如坡顶、韧性剪切带附近)甚至出现拉应力。

随着河谷持续下切，河谷岩体逐渐发生卸荷，河谷边坡应力逐渐发生调整，总体上边坡应力呈现出逐渐降低的趋势。从坡表到坡体内部，卸荷作用逐渐降低，应力降低幅度也逐渐减小；从坡体顶部到河谷不同高程位置，坡体岩体均发生卸荷作用，但高程越高的位置，坡体卸荷作用越明显，表现为应力降低幅度更大，而在河谷位置，由于卸荷作用时间短，卸荷作用不是很充分，应力降低幅度相对而言要小很多。

4 位移场特征

在模拟计算时，上一阶地计算完成后，只是将速度清零，而位移保留到下一阶地的计算过程中，即位移变化可以进行累积。为便于分析位移变化情况，选取上述剖面进行位移分析。随着河谷下切，岸坡岩体产生向河谷方向的位移。岸坡岩体的下切，主要导致坡体发生向河谷方向的位移，因此模型结果只导出了X方向(即垂直河道方向)的位移。

各级阶地下切后的剖面X方向(即垂直河道方向)位移变化情况见图5。从图5中可以看出，位移场有以下特征：从河谷开始下切到现今状态，坡体位移呈现增长趋势，即随着河谷下切，X方向位移逐渐增大，从第一次下切的最大位移约3.6 mm到现今的最大位移7.2 mm；整体来看，坡表的位移要大于坡体内部位移；受韧性剪切带的影响，坡体位移最大范围出现在左右两岸的韧性剪切带范围，在河谷位置的位移相对很小；河谷左岸的位移表现为正值(即位移沿着X正方向产生)，河谷右岸的位移表现为负值(即位移沿X负方向产生)，两者均向河道方向发生变形。由于韧性剪切带岩体物理力学参数相对周围岩体取值要低，在重力作用下，上覆岩体作用在韧性剪切带上，韧性剪切带发生向坡脚方向的变形，类似于发生倾倒变形，韧性剪切带就会挤压下覆岩体。因此，在韧性剪切带附近会出现较大位移，而远离韧性剪切带的岩体，受到韧性剪切带变形所造成的位移逐渐减小，位移相对较小。韧性剪切带在左右两岸均倾向坡体，因此在上覆岩体重力等作用下产生两个相反方向的位移。

a) 第一次下切

b) 第二次下切

c) 第三次下切

d)第四次下切

e)第五次下切

f)现今河谷图5 X方向位移演化过程

5 塑性区特征

各级阶地下切后的剖面塑性破坏区的变化情况见图6，其中tension和shear分别代表的是拉张塑性状态和剪切塑性状态。另外，none表示没有塑性区出现；n表示当前处于塑性状态的单元；p表示过去处于塑性状态的单元。从图6中可以看出，剖面绝大多数范围没有塑性区的分布；随着河谷的下切，坡体塑性区范围逐渐扩大，塑性区主要集中在右岸韧性剪切带范围和右岸坡表局部区域，在左岸韧性剪切带也有少量的分布；随着阶地的下切，塑性区由韧性坡表逐渐向坡体内部发展，这种现象在韧性剪切带内更为明显。在左岸韧性剪切带内，塑性区的分布范围主要集中在低高程位置，可能是由于左岸韧性剪切带的厚度较大、分布的高程要低，上覆岩体自重造成的韧性剪切带产生变形，但未达到破坏所造成。在西部常见的岩质高边坡中，边坡顶部及地形突出部位处的最小主应力往往要比岩体的抗拉强度大，而近于与坡表相互平行的最大主应力基本起不到任何作用，因此在最小主应力的作用下，常常形成与坡面平行的拉裂缝。最小主应力在韧性剪切带和坡顶等部位表现为拉应力，因此在这些部位容易出现塑性区。

a) 第一次下切

b) 第二次下切

c) 第三次下切

d)第四次下切

f)现今河谷图6 塑性区变化

结合河谷下切的应力场、位移场以及塑性区的变化情况来看，随着河谷不断下切，坡体的卸荷作用逐渐增强，表现为塑性区的扩大以及应力的降低，与此同时位移也增加。由于坡体两岸存在着韧性剪切带，对坡体而言，其属于相对软岩带，随着河谷的下切，产生的位移最大值以及主要的塑性区分布带就分布在韧性剪切带周边。在河谷卸荷过程中，由于韧性剪切带的物理力学参数相对坝址区主要岩性英安岩来说取值很小，其在卸荷过程中的应力变化以及变形等较大，卸荷作用对韧性剪切带造成的影响较大，导致在此处的卸荷现象相对坝址区其他位置而言要明显一些。

6 结论

通过对研究区地质情况进行概化，对河谷分为了5次下切来进行数值模拟，得到以下结论。

a) 随河谷下切，边坡岩体产生卸荷作用，造成应力产生重新分布，最大、最小主应力逐渐降低；压应力在坡体中占绝大部分，而有少量拉应力存在于坡表部位；总体来看，在同一高程，坡体内部应力水平大于坡表部位。

b) 在河谷下切过程中，位移呈现随河谷演化逐渐增大的趋势，左右两岸岩体均产生向河谷方向的位移，位移最大值出现在韧性剪切带内及韧性剪切带附近的影响范围内；而在河谷位置，始终为位移较小值分布区域；总体而言，坡表到坡体内部所产生的位移是逐渐减小的。

c) 随河谷的逐渐演化，坡体塑性区分布范围明显比上一级阶地下切前要大，而且塑性区分布深度也在逐渐加深，说明卸荷深度在不断向内部发展；塑性区的分布范围主要集中在右岸韧性剪切带内和右岸坡表位置，在左岸韧性剪切带表层分布有少量塑性区；在坡体表层，岩体更多的是产生剪切塑性区。根据数值模拟结果来看，相较于坝址区左岸，韧性剪切带对坝址区右岸边坡岩体的卸荷影响较大。