基于库水位涨落与降雨耦合作用下新田崩塌体变形分析

刘东泽， 谭 燕， 卢应发

(湖北工业大学土木建筑与环境学院， 武汉 430068)

2020年11月1日，水利部、国家发展改革委公布，三峡工程完成整体竣工验收全部程序，标志着三峡工程正式完工，自三峡工程进行大规模蓄水以来，在防洪与发电等方面发挥了巨大作用，但也带来了一系列问题。其中三峡库区的滑坡与崩塌问题最为严重，其形成原因也十分复杂，大量学者对库区的滑坡、崩塌形成机理进行研究，肖捷夫等[1-2]、江强强等[3]、朱元甲等[4]通过建立大型滑坡实验对滑坡进行分析，认为水位的快速降落与降雨是影响滑坡变形的重要因素，闫国强等[5]、刘新荣等[6]、周剑等[7]对典型滑坡成因进行分析，认为强降雨与库水位升降会使滑坡产生较大变形，王英珍等[8]、夏军强等[9]对黄河、荆江河崩岸建立理论模型，分析多种因素对河岸的影响，还有大量学者利用数值模拟对滑坡的形成因素进行分析[10-16]。

显然，在诸多因素中，库水位变化与降雨对三峡库区的滑坡崩塌形成有重要影响，因此现建立库水位降落联合降雨的耦合数值分析模型，以新田崩塌体为对象，建立二维、三维计算模型，分析崩塌体在不同环境下的变形规律，并建立数值预测模型，预测最危险情况下的变形特征。

1 新田崩塌体概况

新田崩塌体区位于武隆县江口镇，地处乌江右岸，崩塌体前缘为乌江主河道，崩塌体全貌见图1。崩塌体地貌形态为一斜坡，坡体地形起伏较大，崩塌体地形上陡中缓下陡。据野外调查测绘，崩塌体后缘受地形特征及岩土界面控制，呈圈椅状展开，崩塌体高程275 m，纵长250 m，横宽400 m，平均厚度15 m，主变形方向157°。崩塌体两侧边界的划定主要以地形特征变化及变形特征综合确定，西侧边界延伸走向170°，以土坎变形及地形特征、冲沟为界，其外围未见变形现象，而崩塌体区内存在变形。斜坡中上部堆积有厚度较大的崩坡积物，前缘受乌江水不断冲刷形成高陡的临空面，为崩塌体变形提供了空间。并且，崩塌体物质结构较为松散，水流易于下渗，下伏为相对隔水的灰岩、页岩，地下水沿上部崩坡积物与灰岩页岩接触面汇集，大大增加了其接触面的变形性能。

图1 新田崩塌体全貌

2 数值分析模型

2.1 计算参数

根据勘察资料，选取滑坡体的主剖面为计算剖面(图2)。滑体的变形模量、凝聚力和摩擦角，类比周边崩塌体的参数取值，综合确定新田崩塌体数值计算参数具体见表1。

表1 新田崩塌体有限元计算物理力学参数取值表

图2 新田崩塌体地质剖面图

2.2 边界条件

由于崩塌体前缘被库水淹没，则模型的水头边界为浸没的崩塌前缘，崩塌体表面处取降雨引起入渗的流量边界，当降雨强度大于坡面岩土体的入渗速度时，取岩土体的入渗速度值作为边界流量值；当降雨强度小于坡面岩土体的入渗速度时，取降雨强度值作为边界流量值，取降雨强度计算。模型底面和两侧为自由渗流边界，在此由于基岩的渗透性很小，可认为是不透水边界。各边界条件的表达如图3所示。

v为入渗速度，f为产流

2.3 计算模型

2.3.1 二维计算模型

根据新田崩塌体的地质条件和地形地貌特征，利用GEO-SLOPE有限元软件对新田崩塌体的主剖面进行稳定性计算，采用四边形单元进行有限元网格剖分，节点数为3 806，单元数为3 775，网格图如图4所示。

图4 二维计算网格模型图

2.3.2 三维计算模型

为进一步分析崩塌体的变形特征，结合二维计算结果进行三维建模，选取新田崩塌体三维数值计算模型的范围为：沿乌江水流方向为720 m，垂直乌江水流方向为720 m，模型底面高程为50 m。计算域包含崩塌体、与基岩接触带和基岩，整个计算域剖分了38 291个六面体单元，共计43 200个节点，计算模型中的断层采用实体单元模拟，三维计算模型与网格如图5所示。

图5 三维计算网格模型图

2.4 工况设计

三峡工程建成后，水位呈周期性变化，库区水位升级变化过程如图6所示。三峡库区在汛期来临前，水位会降至145 m，在遇到洪水时，水位会达到162 m，洪水过后水位又会迅速下降，在蓄水时，水位能达到175 m，然后又会缓慢下降。

图6 三峡水库运行水位曲线图

武隆县库区地处长江支流乌江两岸，乌江库水位的升降参考长江调度，依据武隆区降雨资料统计分析得到：重现期50 年一遇的汛期日降雨量q全=63.371 mm/d，非汛期日降雨量q枯=45.484 mm/d；重现期20 年一遇的汛期日降雨量q全=59.065 mm/d，非汛期日降雨量q枯=40.571 mm/d。在工况设计时，根据水位调动与降雨资料，设置9种工况对崩塌体的稳定性进行分析，工况组合见表2，同时为分析崩塌体破坏机理，设置5种工况进行三维变形分析，详细组合见表3。

表2 二维计算工况及荷载组合

表3 三维计算工况及荷载组合

3 崩塌体二维稳定性分析

3.1 稳定性计算结果

用Morgenstern-Prince法对新田崩塌体主剖面进行稳定性计算，当蓄水时，对崩塌体的稳定性无影响，即工况1条件下的稳定系数为1.094，工况5条件下的稳定系数为1.044。

计算中考虑饱和-非饱和渗流场与应力场耦合的有限元法，在不同的外部环境下得到新田崩塌体的稳定性系数，水位涨落耦合不同强度降雨得到耦合变化曲线如图7所示。

图7 不同工况下稳定系数变化图

3.2 崩塌体稳定性分析

在水位缓慢下降过程中，由图7(a)可知，在仅有水位变化的条件下，稳定系数变化不大，但在水位下降时遭遇3 d的持续降雨，稳定性系数开始急剧降低，主要原因是库水位降落与降雨引起崩塌体渗流场的改变，同时在降雨早期阶段，降雨强度对稳定性影响并不明显。

在水位降落速度较快时，根据图7(b)可知，降雨强度越大，早期对稳定性影响较大，但随着时间的推移，不同强度降雨对稳定性影响逐渐减小，当在50 年一遇的暴雨时，崩塌体在第7 天稳定性系数就降至1.042，而在降雨强度稍弱的20 年一遇暴雨时，则要达到第9 天。

根据图7(c)可知，在145 m上升至175 m过程中，仅有库水位上涨时，稳定性基本不变，表明水位上涨对崩塌体影响不大，在遭遇暴雨时，早期稳定性出现了短暂的上升，随后稳定性开始降低，总体而言，稳定性降低较小。

联合分析，得到以下结论：工况6为新田崩塌体最危险工况，即水位快速下降与50年一遇暴雨同时出现的情况对崩塌体最为不利，同时在水位上升时可以短暂提高崩塌体稳定性。

4 崩塌体三维稳定性分析

4.1 孔压计算结果及分析

在工况1(水库蓄水至175 m水位)条件下得到孔压计算结果，由图8(a)所示的孔压分布图可知：崩塌体175 m水位以下区域的孔压受水位的直接影响。

结合工况2(175 m水位缓降到145 m)与工况4(坝前162 m水位骤降到145 m水位)孔压计算结果[图8(b)、图8(d)]，分析可知：在没有降雨影响时，崩塌体的孔隙水压力主要受水位下降的影响，在整个崩塌体中，主要受影响的区域为涉水的滑坡前缘，随着水位的降低，崩塌体前缘的地下渗流场发生较大变化，地下水位也随之下降，崩塌体前缘的孔隙水压力消散缓慢，表现为应力滞后现象，崩塌体的中后部区域的孔隙水压力变化不大，可以忽略不计。

图8 不同条件下孔压分布图

在水位下降时遭遇50年一遇降雨的情况下，得到工况3(175 m水位缓降到145 m+50年一遇降雨)与工况5(坝前162 m水位骤降到145 m+50年一遇降雨)计算结果[图8(c)、图8(e)]，结果表明：崩塌体前缘地下水分布受库水位变化和降雨联合作用影响，崩塌体中后部的水压力受降雨影响，但对于崩塌整体而言，孔隙水压力变化不大。

综合图8所示5种工况下孔隙水压力云图可知，水位下降对涉水库岸的孔隙水压力影响较大，在不涉水区域影响不大，降雨主要影响区域为崩塌表层。

4.2 应力计算结果及分析

水库蓄水至175 m水位后大部分区域拉、压应力分布变化较小。根据工况1条件下的主应力场分布云图(图9) 可知，拉应力(第一主应力)主要分布在崩塌体前缘和崩塌体后缘；压应力(第三主应力)主要分布在崩塌体上部、中部和下部(除崩塌体前缘和后缘)的区域，其绝对值的最大值为616 kPa。

图9 工况1条件下应力场分布图

在没有强降雨的条件下，选取工况2与工况4第三主应力云图[图10(a)、图10(c)]，水位由175 m缓降到145 m时，压应力的分布均有所改变，压应力值的绝对值均有所增大，压应力(第三主应力)最大值的绝对值由616 kPa增大为625 kPa；压应力(第三主应力)主要分布在崩塌体上部、中部和下部。库水位由坝前162 m骤降到145 m时，压应力场整体变化较大，压应力(第三主应力)的绝对值最大值为627 kPa；压应力(第三主应力)主要分布在崩塌体上部、中部和下部。

结合工况3与工况5的应力云图[图10(b)、图10(d)]，在库水位下降时遭遇50 年一遇降雨时，压应力的分布均有所改变，压应力值的绝对值均有所减小，库水位由175 m缓降到145 m又遇到50 年一遇降雨时，压应力(第三主应力)最大值的绝对值由625.6 kPa减小为625.2 kPa。库水位由坝前162 m骤降到145 m又遇到50年一遇降雨时，压应力(第三主应力)最大值的绝对值由627 kPa减小为625 kPa。

图10 不同工况下第三应力场分布图

结合孔隙水压力的变化可知，水位变化与降雨作用会导致渗流场的改变，而渗流场与应力场会相互影响，因此在不同工况下应力场的改变主要是渗流场的改变。

4.3 位移计算结果及分析

图11为工况1条件下水平、垂直位移分布图，在蓄水状态下崩塌体下部发生较大水平位移，最大位移值为0.005 m，崩塌体中上部水平位移相对较小，主要原因是崩塌体前缘涉水导致崩塌体下部产生较大变形。在崩塌体上部发生较大垂直位移，方向向下，最大位移值为0.032 m，崩塌体中下部垂直位移相对较小。这主要是由于崩塌体中上部较为陡峭，使得崩塌体在自重作用下产生较大的向下垂直位移，崩塌体下部较为平缓，向下垂直位移较小。

图11 工况1条件下位移分布图

对比工况2～5的水平位移分布图(图12)，四个工况下崩塌体中上部均发生较大水平位移，在没有降雨影响时，水位降落产生的水平位移均比蓄水条件下的位移大，即静止水位下产生的水平位移最小，主要原因是水位下降改变渗流场，使其产生向外的渗透力，导致崩塌体的滑动力增加，由此表现出更大的水平位移，但水位下降速率对水平位移影响不大。在强降雨条件下，崩塌体的渗透力进一步增大，滑动力也随着暴雨的持续而增大，导致工况3、5的位移比工况2、4所得水平位移较大，其中工况5产生的水平位移最大，工况5在崩塌体中上部发生较大水平位移，量值为0.014～0.019 m，最大位移值为0.019 m；崩塌体下部水平位移较小，量值为0～0.014 m。

图12 不同工况下水平位移场分布图

对比工况2～5的垂直位移分布图(图13)，崩塌体上部和中部发生较大垂直位移，方向向下，崩塌体下部垂直位移较小，其中库水位下降叠加强降雨条件下产生的垂直位移最为明显，在工况5条件下崩塌体位移最大，主要是由于库水下降引起崩塌体向下滑移，同时在自重和降雨作用下，使崩塌体产生向下的垂直位移增大。上部最大垂直位移值为0.077 m。在没有降雨影响时，对比工况2、4垂直位移场分布，在工况2崩塌体上部和中部垂直位移量值为0.033～0.042 m，最大位移值为0.043 m，在工况4上部和中部垂直位移量值为0.032～0.048 m，最大位移值为0.048 m，结果表明水位下降速率对崩塌体中上部的垂直位移有影响。

图13 不同工况下垂直位移场分布图

4.4 等效塑性应变分布计算结果及分析

取5个工况下主剖面的等效塑性应变分布计算结果如图14所示，在蓄水状态下，等效塑性应变较大值集中分布在与基岩接触带的中部和上部，等效塑性应变最大值为0.007。在强降雨的影响下，工况3等效塑性应变较大值分布范围较工况2有所增大。等效塑性应变较大值由0.011增大为0.013。工况5等效塑性应变分布范围较工况4有所增大。等效塑性应变较大值由0.011增大为0.015。

综合工况1～工况5的数值计算分析结果，新田崩塌体在工况5(162 m水位骤降到145 m+50 年一遇降雨)条件下，库水位下降与强降雨的影响造成滑体内渗流场发生很大的变化，使坡体内的渗透力(降雨入渗和坡体内外水位差引起)发生很大的改变，进而引起位移发生很大的变化；相比较其他4个工况，新田崩塌体在工况5(162 m水位骤降到145 m+50年一遇降雨)条件下变形增量最大，变形破坏主要发生在崩塌体后部和前部，后部主要发生大范围垂直位移。因此，工况5(162 m水位骤降到145 m+50 年一遇降雨)为对新田崩塌体稳定性最不利工况。

5 数值预测模型

由数值计算分析结果可知库水位由162 m水位骤降至145 m水位遭遇50 年一遇暴雨为对新田崩塌体稳定性最不利工况。提取崩塌体主滑剖面上地表位移监测点XT01、XT02和XT03(图15)在最不利工况条件下的位移变化及滑带代表性位置点的等效塑性剪应变等信息(图16)，通过这些信息的数据拟合，建立崩塌体变形破坏的数值预测模型。

图15 监测点在三维模型中的位置

图16 各监测点位移变化趋势

为了获得崩塌体等效塑性应变的预测模型，在崩塌体与基岩接触处分别选取的A、B、C三个代表点进行分析，A、B、C三点位置如图17所示，等效塑性应变变化趋势如图18所示。

图17 代表点位置分布图

图18 等效塑性应变变化趋势

综合监测位移趋势与等效塑性应变变化趋势可知工况5(库水位由162 m水位骤降至145 m水位+50 年一遇暴雨)条件下，处于崩塌体中后部的地表变形监测点XT03和XT02的位移变化比处于崩塌体前部的地表变形监测点XT01明显；由于崩塌体较陡，所有监测点的垂直位移要比水平位移变化明显。

同时崩塌体与基岩接触处的前中后部三个监测点A点、B点以及C点的等效塑性应变均大于0，说明崩塌体可能发生整体的破坏；C点和B点的等效塑性应变值大于A点，说明崩塌体中后部更容易破坏。

6 结论

(1)根据二维稳定性分析可知，新田崩塌体最危险工况为水位从162 m水位快速降至145 m，同时遭遇50 年一遇暴雨，在该工况条件下，崩塌体的稳定系数为1.042，崩塌体处于欠稳定状态。

(2)在库水位上升的过程中，整体稳定性系数相较于水位下降过程有所提升，即新田崩塌体库水位上升比水位下降时安全。

(3)根据三维变形破坏分析可知，新田崩塌体在工况5(162 m水位骤降到145 m+50 年一遇降雨)条件下，变形最大和等效塑性剪应变增量最大，工况5为对新田崩塌体稳定性最不利的工况。

(4)根据三维变形破坏分析可知，在工况5(162 m水位骤降到145 m+50 年一遇降雨)条件下，新田崩塌体在崩塌体中后部和前缘分别出现局部的拉剪和压剪破坏，但未发生整体破坏。其中，崩塌体中后部滑带较前缘滑带陡，破坏的范围较前缘大。

(5)综合二维、三维计算分析可知，库水位改变对新田崩塌体整体稳定性影响不大，但会导致涉水区域的局部破坏，强降雨对整体稳定性影响较大，在后期监测中，要重点关注暴雨天气。