杨鸿发 许向宁 杨华阳 蹇代君

摘 要： 为了提高有限元法在解决滑坡较大变形时网格畸变的计算效率和精度，以九寨沟区域滑坡处地层土壤材料参数为基础，采用SPH法建立边坡模型，通过改变粒子密度和瑞利阻尼参数，研究了九寨沟区域滑坡的缩尺斜面模型滑坡面发生状况，比较了不同参数下滑动面产生过程中的差异。结果表明，在应用SPH法分析地震边坡变形的过程中，平滑距离的调整易使计算精度发生变化，从而产生一定程度的误差，而平滑距离也会影响边坡的最终变形形态;在不同的SPH粒子密度和不同的瑞利阻尼参数下，滑动面的产生过程也呈现出显著差异。SPH法有效改善了有限元处理大变形时网格失真等问题，研究结果可为地质灾害易发区的滑坡灾害防治提供一定参考。

关键词： 岩土力学;九寨沟;滑坡;SPH;粒子密度;瑞利阻尼

中图分类号： P315   文献标识码：  A

doi：  10.7535/hbgykj.2020yx06007

Simulation study of Jiuzhaigou landslide based on SPH method

YANG Hongfa 1， XU Xiangning 1， YANG Huayang 1， JIAN Daijun 2

（1. 405 Geological Team， Sichuan Provincial Bureau of Geology and Mineral Exploration and Development， Chengdu， Sichuan  611830， China; 2. Jiuzhaigou Scenic Spot Administration， Aba Tibetan and Qiang Autonomous Prefecture， Sichuan  623400， China）

Abstract：

In order to improve the calculation efficiency and accuracy of the finite element method in solving mesh distortion of large landslide deformation， based on the soil materials parameters of the landslides in the Jiuzhaigou area， the SPH method was utilized to establish the slope model. By changing the particle density and Rayleigh damping parameters， the landslide surface occurrence of the scaled slope model of Jiuzhaigou lanolslides was studied and the differences in the sliding surface generation process under different parameters were analyzed. The results show that in the analysis of seismic slope deformation by SPH method， the calculation accuracy is easy to change when adjusting the smoothing distance， which results in a certain  degree  of error， and the final deformation of the slope can be affected by the smoothing distance. With the different SPH particle densities and Rayleigh damping parameters， the generation process of the sliding surface also shows significant differences. SPH method can improve the mesh distortion of large deformation by finite element method， and the results of this study can provide a reference for landslide diaster prevention and control in other geological disaster prone areas.

Keywords：

geotechnical mechanics; Jiuzhaigou; landslide; SPH; particle density; Rayleigh damping

山區地震发生时一般伴随不同程度的崩塌、滑坡和泥石流等灾害。2017年8月8日，四川省阿坝藏族羌族自治州九寨沟县发生了7.0级地震，震中位于九寨沟县漳扎镇世界自然遗产九寨沟风景名胜区内，此次地震震源深度为20 km，最大烈度为9度，详细的九寨沟地震烈度图如图1所示  [1] 。此次地震共造成四川、甘肃2省8个县受灾，25人死亡，525人受伤，6人失联，176 492人受灾， 73 671 间房屋不同程度受损。多个自然遗产点（如火花海、树正群海、五花海、熊猫海）遭到不同程度的破坏，尤其是地震活动引起景区内大量山体出现不同程度的崩滑，地形，地貌，以及大面积森林植被、耕地和草地遭到严重损毁。10年内2次强烈地震破坏的叠加作用，致使九寨沟地质灾害隐患更加突出，对九寨沟景区景观的生态安全构成了巨大威胁  [2-3] 。

针对地震诱发的崩塌、滑坡，常用的研究方法有数值分析法、人工智能法、监测法等  [4-5] 。通过数值分析法预测并分析地震引起的滑坡崩塌具有重大意义。有限元法是常用于滑坡解析的数值分析法  [6-8] ，虽然有限元法在预测边坡稳定性方面具有良好的效果，但采用有限元法解析九寨沟边坡大变形时出现了网格失真、精度大幅降低等问题，对于边坡的最终变形状态和滑动破坏范围的预测还不够充分。

SPH（smoothed particle hydrodynmics）法将连续体看作多数粒子的集合，属于粒子分析法，是一种适合处理大变形问题的解析方法，主要应用于流体分析领域。作为对压缩性流体的分析方法，SPH法最初被用于宇宙物理学领域  [9] 。由于不需要网格并具有适合大变形解析等优点，学者们开始将其运用于普通流体的解析，后来在一般弹塑性体的大变形解析中也加入了SPH法  [10-12] 。BUI等  [13] 将Drucker-Prager准则导入SPH法解析中，解析结果与砂土模型实验结果高度相似。NAILI等  [14] 将液化后的地基视为黏性流体并使用SPH法解析，成功分析出流动地基对建筑物的地下构造物的荷载作用。

因此，在本研究中利用SPH法代替有限元法对地震时的边坡状态进行了分析验证，并提出了瑞利阻尼的导入方法。瑞利阻尼是一种广泛采用的正交阻尼模型，与频率相关。在进行土层地震反应时有较好的解耦性能。

通过对九寨沟滑坡发生处的地层土壤材料建立边坡模型并进行解析，確认了粒子密度（SPH法中特有的参数）以及瑞利阻尼在边坡崩塌过程中的影响。

1 SPH法的理论基础

1.1 基本理论

SPH法是以连续体作为研究对象，使用多数粒子集合的近似计算方法。一般进行核近似和粒子近似2种近似计算方式。核近似如式（1）所示。

SPH粒子和核函数的详细关系可参考图2，求解物理量 f时相当于求出离中央的黑色粒子一定距离内的粒子（图中的深色粒子）所具有的值的加权平均值，在浅色粒子的位置 ，核函数的值为零则不含在计算中。

平滑距离 是常用参数，取值过小易导致SPH近似计算的粒子数减少，从而导致计算精度的降低。为了解决计算误差问题，使用CHEN等  [15] 提出的CSPM（corrective smoothed particle method）法，将式（2）和式（3）用式（7）和式（8）代替：

1.3 弹塑性本构方程及应力增量

在 SPH法弹塑性分析中，可以使用有限元法中求应力增量的弹塑性矩阵法。形变速度张量和速度梯度由式（13）—式（15）求得：

1.4 瑞利阻尼的使用

在地震分析中经常会用到瑞利阻尼，有限元法中的瑞利阻尼矩阵[ C ]利用系数αR及βR表示 如下：

此时，式（17）中的刚性矩阵[ K ]在分析对象时相当于切线刚度矩阵。

2 模型和解析条件

2.1 模型建立

图3为典型地震滑坡高分影像特征，九寨沟滑坡发生处地层土壤材料的物理性质见表1。图4中的圆弧线为使用表1数据进行的经典斜面模型实验所得到的滑动面，九寨沟景区熊猫海沿岸发生的滑坡按等比例缩小的模型实验如图5所示，实曲线即为模型实验中产生的滑动面。笔者采用4种不同粒子排布密度的模型进行解析，如图6所示。SPH粒子用规则格状排布，粒子 间隔定义为dp，粒子总数定义为np。边界条件为底面固定，侧面只在水平向固定，粒子可以上下方向自由运动，平滑距离h为2.6dp。形变速度和加速度由上述各式分别计算得到。瑞利阻尼在一次模式和二次模式中使用相同的衰减常数。为了考 察不同衰减对分析结果产生的影响，使用瑞利阻尼分别为2%，5%，10%时进行对比分析。

2.2 边坡自重作用分析

为重现边坡的崩塌过程，对模型进行了自重作用分析。解析中因为弹塑性本构方程的使用，在震动时有可能发生粒子塑性化。图7为重力加速度逐渐增加时的分析。为了减少自重作用分析所需的时间，采用BUI等  [16] 提出的衰减项代替运动方程：

式中 ξ为衰减强度参数，一般为0.002～0.005，本文采用ξ=0.005。

图8表示了坡顶角垂直方向速度随时间发生的变化。在粒子密度不同的模型中，衰减项均发挥了充分的效果，可以确认在分析结束时粒子不会产生震动。而图9表示了自重作用分析结束时沿垂直方向的轴应力 σ  yy 的分布。所有模型的结果都很相似，对于初始应力状态，粒子密度没有太大的差异。

2.3 粒子密度对地震解析结果的影响

水平向输入的地震加速度如图10所示，地震作用下解析得到的边坡变形情况如图11所示，在震动开始后经过10 s， dp =0.5 m时，最大剪切变形出现的地方和经典模型试验结果非常相似。 dp =1.0 m和2.0 m时，最大剪切变形出现的地方与经典斜面模型试验相比相对靠下。 dp =5.0 m时没有出现明显的滑动线。如图11 b）所示，开始后经过30 s震动结束时， dp 为0.5，1.0，2.0 m时，震后得到的结果形状几乎相同。 dp 为10， 2.0 m时，剪切变形超过60%的区域基本一致。但在 dp= 0.5 m时，在较深的位置也出现了剪切形变超过60%的滑动线并延伸到了顶部。另一方面， dp =5.0 m时， 坡趾处出现剪切变形超过60%的区域但没有延伸到顶部，与其他情况有明显差别。

从解析结果可以看出，利用SPH法对边坡崩塌的分析中，根据设定的不同粒子密度，其结果也会受到一定影响。

dp 为1.0 m和2.0 m时，得到了相似的结果。当粒子密度更大（ dp =0.5 m）时，出现了更多的滑动线，由于分析条件的差异，可以得出完全不同的崩塌状态。经典斜面模型实验结果是边坡在图4和图5所示圆弧线处滑落，而解析结果更倾向于整个坡面崩塌。

2.4 瑞利阻尼的作用對地震结果的影响

图12为瑞利阻尼分别为2%，5%，10%时的解析结果。震后10 s和震后30 s最大剪切变形出现的地方相对一致。瑞利阻尼为2%，5%，10%时，剪切形变超过60%的区域全部延伸到了顶部。随着瑞利阻尼从10%下降到2%，剪切形变超过60%的部分在逐渐增多，可以清晰地看到至少产生2个滑动面。从图12可以观察到，施加的衰减越小，边坡的剪切变形越大，并且区域分布也不相同。从解析结果看，滑动面的下部情况与经典斜面模型实验相对吻合，滑动面上部最大剪切变形部分仍有待分析。综上所述，由于瑞利阻尼可导致不同斜面的崩溃，需要进一步对使用的多个参数进行研究。

3 结 语

针对九寨沟区域分布的滑坡进行了分析，采用相对于传统有限元法更为精确的SPH法建模并与室内斜面模型实验结果进行了对比，研究了不同粒子密度以及瑞利阻尼作用对于边坡滑动的影响。

1） 模型解析时，通过采用不同的SPH粒子密度和不同的瑞利阻尼，滑动面的产生过程随之呈现出显著区别。粒子密度为1.0 m和2.0 m时，计算结果与模型实验较为相似。解析时施加的瑞利阻尼越小，边坡的剪切变形越大。

2）采用SPH法对地震边坡分析时，调整平滑距离易使计算精度发生变化，从而产生一定程度的误差，而平滑距离会显著影响边坡的最终变形，因此计算时需要充分考虑使用合理的平滑距离。

3）利用SPH法对九寨沟发生的真实滑坡进行了缩小模型的再现分析，使用单纯的粒子解析避免大变形失真问题，结果可为九寨沟等地质灾害易发区域的滑坡灾害防治提供重要参考。

虽然本研究对粒子密度和瑞利阻尼进行了一定的参数分析，但是SPH法模拟结果与斜面模型实验产生的滑动面上部未能较好地吻合，原因可能在于各类参数的适用性范围还存在一定的差别。在今后的研究中还需要对粒子密度、瑞利阻尼等各种影响参数的适用范围作进一步分析。

参考文献/References：

[1]  徐锡伟，陈桂华，王启欣，等.九寨沟地震发震断层属性及青藏高原东南缘现今应变状态讨论[J].地球物理学报，2017，60（10）：4018-4026.

XU Xiwei， CHEN Guihua， WANG Qixin，et al. Discussion on seismogenic structure of Jiuzhaigou earthquake and its implication for current strain state in the southeastern Qinghai-Tibet Plateau[J]. Chinese Journal of Geophysics， 2017，60（10）：4018-4026.

[2]  房立华， 吴建平， 苏金蓉， 等. 四川九寨沟Ms7.0地震主震及其余震序列精定位[J]. 科学通报，2018，63（7）：649-662.

FANG Lihua， WU Jianping， SU Jinrong，et al. Relocation of mainshock and aftershock sequence of the Ms7.0 Sichuan  Jiuzhaigou  earthquake[J]. Chinese Science Bulletin，2018，63（7）：649-662.

[3]  李勇， 邵崇建， 李芃宇， 等. 九寨沟Ms 7.0级地震的左旋走滑作用与动力机制[J]. 成都理工大学学报（自然科学版）， 2017， 44（6）：641-648.

LI Yong，SHAO Chongjian，LI Pengyu，et al. Left-lateral strike-slip effect and dynamic mechanism of the Jiuzhaigou  Ms 7.0  earthquake in the eastern margin of Tibetan Plateau， China [J]. Journal of Chengdu University of Technology （Science and Technology Edition），2017， 44（6）：641-648.

[4]  黎玺克.遗传算法优化BP神经网络的岩质边坡稳定性预测[J].河北工业科技，2020，37（3）：164-169.

LI Xike. Prediction of rock slope stability based on BP neural network optimized by genetic algorithm[J]. Hebei Journal of Industrial Science and Technology， 2020，37（3）：164-169.

[5]  于远亮，靳静.公路高边坡稳定性远程监控技术研究[J].河北工业科技，2014，31（3）：239-242.

YU Yuanliang， JIN Jing. Study on remote monitoring technology about highway slope stability[J]. Hebei Journal of Industrial Science and Technology， 2014，31（3）：239-242.

[6]  WATKINSON I M， HALL R. Impact of communal irrigation on the 2018 Palu earthquake-triggered landslides[J]. Nature Geoscience， 2019， 12（11）： 940-945.

[7]  FRANCI A， DE-POUPLANA I， CELIGUETA M A， et al. Application of the particle finite element method （PFEM） to the numerical simulation of landslides[C]// International Conference of Natural Hazards and Infrastructure （ICONHIC）. Athens： National Technical University of Athens， 2019：  1-10.

[8]  WANG  Liang， ZHANG Xue， ZANIBONI F， et al. Mathematical optimization problems for particle finite element analysis applied to 2D landslide modeling[J]. Mathematical Geosciences， 2019.

doi： 10.1007/s11004-019-09837-1.

[9]  MONAGHAN J J. Simulating free surface flows with SPH[J]. Journal of Computational Physics， 1994， 110（2）： 399-406.

[10]  HU Wei， GUO Guannan， HU Xiaozhe， et al. A consistent spatially adaptive smoothed particle hydrodynamics method for fluid-structure interactions[J]. Computer Methods in  Applied  Mechanics and Engineering， 2019， 347： 402-424.

[11]  BAYAREH  M， NOURBAKHSH A， ROUZBAHANI F，  et al . Simulation of sand particles flow using a weakly compressible smoothed particle hydrodynamics method （WCSPH）[J]. Annales de Chimie： Science des Materiaux， 2019， 43（1）： 43-51.

[12]  ONO  Y， OKAMOTO R， KOHNO M， et al.3D simulation of the aratozawa landslides base on smoothed particle hydrodynamics method[J]. Journal of Japan Society of Civil Engineers， 2017， 73（4）： 346-356.

[13]  BUI  H H， FUKAGAWA R， SAKO K， et al. Lagrangian meshfree particles method （SPH） for large deformation and failure flows of geomaterial using elastic-plastic soil constitutive model[J]. International Journal for Numerical and Analytical Methods in Geomechanics， 2008， 32（12）： 1537-1570.

[14]  NAILI M， MTSUSHIMA T， YAMADA Y.A 2D smoothed particle hydrodynamics method for liquefaction induced lateral spreading analysis [J]. Journal of Applied Mechanics， 2015， 8：591-599.

[15]  CHEN J K， BERAUN J E， JIH C J. Completeness of corrective smoothed particle method for linear elastodynamics[J]. Computational Mechanics， 1999， 24（4）： 273-285.

[16]  BUI H， FUKAGAWA R， SAKO K. A study of the matter of SPH application to saturated soil problems[C]// Proceedings of the 5th International SPHERIC Workshop.Manchester：  University  of Manchester， 2015： 354-361.