基于流固耦合的强震大型滑坡水力激发效应研究

2023-08-31周先强

水科学与工程技术 2023年4期

周先强

（福建省水投勘测设计有限公司，福州 350001）

为研究地震作用下滑坡的动力响应规律，本文通过有限元法实现动力和渗流的耦合，建立滑坡数值模型。通过前人试验研究验证模型参数的正确性，对是否考虑地下水工况下的滑坡加速度和位移响应规律进行计算，研究孔隙水压力的分布情况，并对软弱层带的有效应力和孔压响应进行分析。

1 流固耦合模型

1.1 基本理论

本文基于完全非线性的有限差分法，通过有限元法实现动力和渗流的耦合进行分析。边坡土体模型选用Finn 模型，在动力作用下，相对于摩尔-库伦模型而言，Finn 模型增加动态孔隙水压力程序，能够模拟材料孔隙水压力的积累及超孔隙水压力的形成［5］。

本文在进行非完全流固耦合分析时考虑到力学时标与扩散时间之间的联系和流体的扰动属性，对流固刚度比进行合理选择。在地震作用下，相对总体计算过程而言其力学计算时间较短，因此强震作用下可不考虑边坡渗流的作用，进行不排水计算分析。当模型整体失衡时，在耦合计算时需要流固刚度比。本文研究的滑坡类型属于刚性骨架，即流固刚度比小于1，固体刚度远大于流体刚度，因此可以对流体与固体的完全耦合进行忽略。

1.2 模型尺寸和监测点布置

本文通过有限元法建立滑坡数值模型，该模型长度设置为4500 m，倾角18°，高度为2400 m，其中，滑坡的软弱层带的厚度约5 m。根据相关地质勘察报告，该滑坡主要土体为白云岩，因此，将白云岩作为滑坡数值模型的整体。通过孔隙介质体来对软弱层带进行简化，该孔隙介质体仅需要设定角砾岩带参数，对下部地层的水文条件进行忽略。而软弱层带与地下水位线之间的坡体区域认为均达到饱和状态。根据实际滑坡情况，并参考相关文献，设置模型底部为固定端，模型两侧为法向固定，模型上部为自由端，地震波由模型底部进行输入。此外，为在流固耦合分析过程中保证计算精度的前提下兼顾计算效率，在网格划分时对软弱层带及坡面位置进行局部加密，共划分出34551 个10 节点的有限元网格。

为了更好地观测地震作用下滑坡的动力响应规律和变形特性，在模型中总共布置有15 个监测点进行监测。其中，监测点J1～J6 共6 个监测点设置在沿软弱层带分布，相邻监测点间距设置为500 m，对J1、J3 和J6 监测点位置每个位置处再布置有3 个垂直监测点，垂直间距为20 m，即K1～K9。此外，将观测点分为3 个区，以便更好地观测软弱层带的孔隙水压力变化情况。

1.3 模型参数设定与验证

模型上下硬层的泥质砂岩本构模型为弹性模型，弹性模量为5.8×104MPa。渗流计算过程中的材料颗粒通常认为是不可压缩的，并假设坡体为各向同性材料模型。通过下式可计算出材料的剪切模量和体积模量：

式中G 为泥岩的剪切模量；K 为泥岩的体积模量。

依据相关模拟方案的参数设定［6］和实际工程经验，本文采用的主要材料泥质砂岩和软弱层带的模型输入参数如表1。

表1 软弱层带和白云岩数值模型输入参数

通过利用表1 数据进行动三轴试验的模拟，围压和循环偏应力均选择为100 kPa，与室内试验条件完全相同，将模拟结果与前人试验结果［7］进行对比，图1 展示数值模拟与室内试验的孔压比曲线对比情况。从图1 可看出，本文数值模拟结果与前人试验研究结果趋势和数值都较为接近，随着动力持续时间的增加，并未出现明显差异，进而证明本文在参数选取和本构模型选择方面合理、正确。

图1 动三轴试验孔压比曲线

利用上述模型对地震荷载施加前的竖向应力和孔压进行计算，计算结果显示，竖向应力分布主要受高度影响，主要分布在滑坡底部，未出现明显波动，在坡脚处存在少量的应力集中现象，软弱层带的应力约在13～17 MPa 之间。由于软弱层带下方设置为不排水条件，因此不存在孔压，而在软弱层带上方，孔压亦随水位以下的深度分布。

1.4 加载方案制定

地震作用下，地震波的吸收是通过模型四周的自由场边界来完成的。局部阻尼系数取0.015，局部阻尼可消耗模型中的地震能量。模型的自振频率为0.25 Hz，在地震作用下可能会发生共振。

动荷载施加位置位于模型底部，输入荷载为加速度时程。可通过下式对滑坡滑动方向地震波进行换算：

式中aH（t）为滑坡滑动方向的水平加速度；aEW（t）为清平波东西加速度分量；aEW（t）为清平波南北加速度分量。

为了对地震作用下模型进行优化，减小计算时间提高计算效率，截取天然波主震部分为模型输入地震波，即13～33 s 段的地震波进行输入分析。

2 结果与讨论

2.1 滑坡动力响应

无水位工况下滑坡加速度响应情况如图2。图2（a）为垂直波作用下的响应曲线，图2（b）为水平波作用下的响应曲线。从图2可看出，相较于上硬层，滑坡下硬层在地震作用下的加速度响应更为强烈，变化幅值更大，最大幅度在动力时间15 s 左右。对比图2（a）和图2（b），同一动力持续时间和位置处，水平波造成的坡体加速度变化相较于垂直波更甚。

图2 无水位工况下加速度响应曲线

有水位工况下滑坡加速度响应情况如图3。图3（a）为垂直波作用下的响应曲线，图3（b）为水平波作用下的响应曲线。从图3可看出，在两种地震波作用下，其上硬层的加速度振幅和峰值均大于下硬层，这一点与无水位工况下的现象恰好相反。对比图3 和图4 发现，水位对滑坡加速度响应的主要影响在上硬层，有水位工况下的上硬层峰值相较于无水位工况提高3 倍之多，可见地下水对坡体动力响应影响程度之巨。此外，地下水的存在，使得水平波和垂直波之间的差异减小，滑坡对垂直波更为敏感。

图3 有水位工况下加速度响应曲线

图4 无水位工况下瞬时位移相应曲线

图4 和图5展示有无水位工况下滑坡在两种方向地震波作用下的瞬时位移情况。从图4 可看出，无论是垂直波还是水平波，下硬层的位移变化幅度和峰值均大于上硬层。对比图4(a)和图4(b)发现，水平波造成的滑坡瞬时位移相较于垂直波而言更大。从图5 中可看出，亦呈现出上文规律，即存在地下水时，上硬层对地震波的响应更为剧烈，上硬层的位移变化幅值和峰值均大于下硬层。此外，无地下水存在时，滑坡位移对水平波更为敏感，而存在地下水时，滑坡对垂直波更为敏感。

图5 有水位工况下瞬时位移相应曲线

2.2 超空隙水压力积累与激发

从上节可看出，在有水位工况下，滑坡对垂直地震波更为敏感，因此，选用垂直地震波作为输入波，对孔隙水压力分布情况进行分析。

孔压主要分布在软弱层带附近位置，最大孔压在B 区域监测点J3 附近，相较于无震工况下的孔压分布，最大孔压位置并未发生改变。随着动力持续时间的增加，孔压持续增大，最大增加到约30 MPa，相较于为施加地震荷载作用前的4 MPa，地震无疑大大增加滑坡孔压的大小。

图6 展示地震荷载作用下软弱层带孔压响应随动力持续时间的变化情况。从图6 可看出，上下硬层的应力较为接近，且随动力持续时间的增加变化趋势亦未出现明显差异。

图6 软弱层带孔压响应

随着动力持续时间的增加，软弱层孔压不断波动，较大的两次波动发生在4 s 和12 s 时刻，在4 s左右时，由于振动张拉导致孔压出现下降，而在12 s时刻附近，由于振动冲压导致孔压出现上升，此时软弱层孔压达到峰值，接近于上下硬层的孔压。

2.3 有效应力

为进一步研究地震作用下的滑坡动力响应规律，更为真实地模拟地震动环境，进行双向地震动模拟。

图7 展示了双向地震动作用下软弱层带有效应力的相应情况。从图7 可看出，在地震荷载施加初期，软弱层带有效应力数值较小且起伏不大，表现较为稳定，出现缓慢降低趋势，在7.5 s 之后开始发生较大波动，主要表现在中前段监测点J5 位置处，中后段即监测点J1 位置处在12 s 时刻左右发生较大突变，有效应力大小出现较大增加。在动力施加时段内，有效应力峰值最大的为软弱层带中前段即J5 监测点位置处。在J1 位置处孔隙水压力大于上覆应力，因此数值表现为负，而软弱层带中后段有效应力绝对值最小，有效应力较低。此外，从数值上看，滑坡中后段有效应力为负，超孔压会引起方向向上的内力，对滑坡造成托顶效应，有可能导致滑坡灾害发生，存在一定安全隐患。