地震是地壳在受力破裂后快速释放能量的过程中产生振动并传播地震波的一种自然现象，是常见的自然灾害之一。由于迄今为止在技术层面仍然难以准确地预测地震，为了降低地震造成的损失，评价地震可能诱发的灾害形式并进行合理的抗震设计是当前最为有效的应对方法。地震发生时，震源处释放的能量往往以地震波的形式传至地壳表面，大量震害调查结果表明，剪切波和以瑞利波（Rayleigh wave）为代表的面波更具破坏性，其传播形式如图1所示。在一些特定的地质条件下，地震波引起的岩土体剧烈运动会产生流滑大变形，一般认为，在坡度小于3°的缓倾地基中这类变形相对有限。然而，在2018年苏拉威西地震中，坡度仅有1°～3°的缓倾场地却出现了数百米的流滑变形[1,2]。大量振动台离心模型试验证明[3]，此时仅输入剪切波不足以产生如此规模的变形，Li等[4]发现瑞利波可能在缓倾场地的侧向变形中起重要作用，在抗震设计中有必要考虑瑞利波等面波的影响。

图1剪切波和瑞利波的传播形式 

当涉及大变形问题时，传统的有限元方法中网格往往会产生扭曲和畸变，并导致刚度矩阵和质量矩阵变得接近奇异。此时动力计算的时间步长必须降低到相当小的程度，线性方程组的解也较不稳定，不论是计算精度还是效率都难以接受。为了解决上述网格畸变问题，许多学者尝试了不同的改进方法来处理岩土工程中的大变形问题，比如光滑粒子流体动力学法（SPH）、粒子有限元法（PFEM)[5]、物质点法（MPM)[6]、无网格伽辽金法（EFG)[7]等。其中，SPH避免了繁琐的数值积分，且可以较好地反映自由表面的运动情况，因而非常适合处理流滑大变形问题。

SPH最早于1977年由Gingold等[8]和Lucy[9]分别独立提出，是应用最早的无网格数值方法之一，广泛应用于天体物理学、流体动力学、冲击动力学、岩土工程等领域。Maeda等[10]最早使用SPH研究了岩土体的大变形问题，但没有考虑塑性变形，仍存在诸多局限性。Bui等[11]采用了Drucker-Prager理想弹塑性模型，充分讨论了使用SPH模拟岩土体的塑性变形时所遇到的各类问题。在此基础上，国内外许多学者都使用SPH研究了岩土工程中的各类大变形问题[12,13,14,15,16,17,18]。受计算能力的制约，早期的SPH程序一般只能处理数千个粒子，难以准确模拟实际震害现象。但随着GPU加速和并行计算技术的发展，现如今SPH已经可以求解千万级自由度的数值模型[16,17,18]，足以解决常见的岩土体大变形问题。

综上所述，为了研究瑞利波对缓倾场地流滑大变形的影响，本文使用SPH进行了数值模拟。首先介绍了SPH大规模计算平台的构建方法，随后通过二维平面应变模型与基于2018年苏拉威西地震Balaroa滑坡构建的三维模型，计算了瑞利波与剪切波输入下缓倾场地的变形特征，对比分析了瑞利波与剪切波对此类变形的影响，揭示了瑞利波在缓倾场地流滑大变形中所起的重要作用。

1、SPH数值分析方法

1.1控制方程

在连续介质动力计算中，控制方程为如下连续性方程和平衡方程：

式中：u为位移矢量；x为位置矢量；σ为总应力张量；g为重力加速度矢量；f为体力矢量。对于大部分较复杂的动力计算问题，直接求解上述偏微分方程组的解析解是非常困难的。此时可以使用紧支试函数加权残量法将该方程组转换为弱形式，从而使用数值计算方法求解计算域Ω内u(x)的近似解uh(x)，即

式中的近似解uh(x)通常可以表示为一组基函数的线性组合：

式（3）中W k (x)称为权函数，式（5）中N i (x)称为试函数，iu为离散点处的位移矢量。通过选取不同的权函数、试函数以及数值积分方法，可以得到不同的无网格数值计算方法。当权函数采用配点法、试函数采用核近似、积分方法采用梯形公式时，得到的数值方法即为广义的SPH，此时特定节点xI处的物理量uI可以表示为

式中：wIJ=w(x I,x J,h)为核函数；为点xI与xJ之间的距离；h为选定的支撑域半径。本文选用如下形式的四次样条函数作为核函数：

将式（6）代入式（1)～(4），即可得到控制方程的SPH近似形式：

式中，RIJ为提高数值稳定性所施加的人工黏滞阻力，本文使用如下形式的层流黏性力：

式中，υ0为运动黏度，取固定值0.05 m2/s。

1.2本构模型

从物理状态而言，发生流滑变形的土体可视为介于固态和液态之间，因此简化地使用流变学模型进行数值模拟通常能够较为全面地反映其物理性质。由于土体发生流滑变形时密度可能会发生较大的变化，本文采用弱压缩性流体的状态方程来描述平均压应力p（以压缩为正）与质量密度ρ之间的关系，表达式为

式中：ρ0为参考状态的密度，为ρ=ρ0时的胀缩波波速；β为常数。对于弱压缩性流体，建议取β=0.7。偏应力张量s与偏应变张量e间的关系选用Herschel-Bulkley-Papanastasiou流变学模型（简记为HBP模型）进行描述[19]，即

式中：为广义剪应变率；τ0为初始屈服应力；μ为运动黏度；m和n为材料参数，对本构关系的影响如图2所示。可以发现，当m较大时（比如取m=100）可以认为τ0为土体静止状态下的屈服应力，当n=1时，材料为典型的Bingham流体。Han等[17]证明了这套本构模型可以较为准确地反映饱和土弱化后的运动特征，并且数值稳定性很好。

1.3边界条件

边界条件的处理方法以及边界附近粒子的计算精度问题一直以来都是SPH的应用难点，常用的方法有虚粒子边界条件[20]、排斥粒子边界条件[21]、速度无滑移边界条件[22]、动力粒子边界条件[23]等。

动力粒子边界条件（Dynamic Boundary Particle Conditions,DBP）即在计算域外再添加若干粒子来反映外边界，这类粒子可称为“边界粒子”。在求解内部粒子的运动方程时，边界粒子就像其他常规粒子一样参与计算，但其本身的运动状态是由边界条件给定的。与此同时，边界粒子提供的压应力p对密度ρ的变化非常敏感，从而提供较大的斥力，来避免内部粒子穿过边界粒子层。本文只需要在边界上输入地震波的位移时程，因此选用动力粒子边界条件。

图2参数m与n对本构关系的影响 

DBP边界条件物理概念清晰，程序设计简单，对于给定位移的边界条件可以得到较好的处理效果，但在使用时需要注意以下两个问题：(1)时间步长Δt不能取得太大，否则内部粒子在一个时步内运动距离较远，可能会穿过边界粒子层；(2)支撑域半径h不能取得太大，边界附近的内部粒子的核函数支撑域不能越过边界粒子层，否则支撑域两侧不同粒子的权重差异较大，会严重影响计算精度。

1.4大规模计算平台

使用SPH求解偏微分方程时，支撑域内粒子的数量对计算结果的影响很大，而且为了确保边界附近的计算精度，支撑域半径h也不能取的太大。因此，构建SPH模型时粒子的间距不能取得太大，Peng等[18]证明了在模拟地基变形时粒子间距应小于2 m，否则可能会出现较大的误差。因此，在使用SPH模拟实际工程中的大规模地基变形时，往往需要构建上千万级自由度的数值模型，需要搭建大规模计算平台。

本文使用GPU加速以及并行计算的技术来处理此类大规模计算问题。与CPU相比，虽然GPU不适合处理复杂的计算任务，但数据的吞吐量要大很多，特别适合处理大规模且反复运行的简单数学运算，与SPH的特征不谋而合。Crespo等[13]证明了在SPH模型计算中，GPU加速技术可以有效地提高计算效率，最多可以提升70余倍，并且非常适合并行计算。

本文使用Dual SPHysics开源平台[13]进行动力计算，该平台集成了HBP本构模型和DBP边界条件，且支持GPU并行计算。并行计算则通过北京超级云计算中心提供的GPU并行服务实现，最多可同时调用8块计算显卡，共计32 GB显存，可以在24 h内完成本文涉及的千万级自由度SPH动力计算。

2、瑞利波对流滑大变形的影响

2.1模型构建

为了研究瑞利波对地基流滑大变形的影响，本节构建了如图3所示的水平与缓倾地基的平面应变模型，并使用Dual SPHysics软件进行动力计算，对比剪切波与瑞利波输入下地基的变形特征。在二维模型中，上层覆盖着10 m厚的软弱土层，下层为较坚硬的土层，水平方向宽100 m，倾斜场地的倾斜角为2°。

正在加载图片

图3理想双层地基的平面应变SPH模型 

Fig.3 Idealized two-layer plane strain SPH ground model

平面应变模型中粒子初始间隔d=0.1 m，支撑域半径h=0.173 m；状态方程中的参考状态密度ρ0=1.0 g/cm3，胀缩波波速参数β=0.7；依照Han等[17]根据泥浆的水槽试验给出的建议值，两种土本构参数取值如表1所示。

表1软弱土层与坚硬土层的本构参数

模型的左侧边界（即x=0 m）上输入瑞利波，底部边界（即z=0 m）上输入剪切波，位移时程如图4所示。其中，剪切波沿x方向是均匀分布的，瑞利波沿z方向是非均匀分布的，在z=0 m处位移保持为0，其余深度处则如图4所示，且两种波的平均能量密度相等。根据瑞利波的解析解，地面附近瑞利波竖直方向的运动产生的影响相对较小[4]，简化起见，只输入了水平方向的运动分量。此外，在实际地震中，瑞利波和剪切波往往是同时存在的，因此本节还考虑了同时输入瑞利波和剪切波的工况（记为“复合波”）。

图4平面应变模型中输入的位移时程 

2.2计算结果

不同地震波输入下水平与倾斜场地中软弱土层的水平残余位移分布情况如图5所示，可以发现，虽然各种输入工况下倾斜场地整体会比水平场地多出约2m的水平残余位移，但总体的分布规律是类似的。剪切波输入下两种模型的整体变形幅度都较小，且对于水平场地残余位移基本对称，最大正向位移为8.2 m，最大负向位移为-9.4 m，分别发生在场地左侧地表与右侧地表；而对于倾斜场地，最大正向与负向位移分别为8.6,-8.0 m，相差不大。

图5不同地震波输入下平面应变场地的水平残余位移云图 

瑞利波和复合波输入下水平场地变形状态都是非对称的，变形集中在左侧，最大正向位移分别为49.1,53.2 m；倾斜场地中最大正向位移则分别为60.0,59.0m，变形规律与水平场地基本一致，且变形都集中在地震动输入侧20 m范围内，说明地震波在HBP黏塑性流体中衰减得较快，直接影响范围有限。

综上所述，在使用HBP流变学模型和DBP边界条件的SPH数值模型中，不论是水平场地还是缓倾角的倾斜场地，仅输入剪切波时都不足以产生足够大的永久变形，而输入瑞利波后场地中则会出现较大的侧向变形，说明抗震分析中有必要考虑瑞利波的影响。

3、Balaroa滑坡数值分析

3.1工程背景

2018年苏拉威西地震诱发的流滑大变形是近年来岩土地震工程领域的研究与讨论热点，尤其是Jono Oge、Sibalaya、Petobo和Balaroa四处滑坡，发生滑动处地面坡度仅有1°～3°，但滑动的最大位移却达到了数百米甚至上千米。其中，Balaroa滑坡由于发生在帕鲁市的居民区，造成的人员伤亡和经济损失远比其余3处更大。许多学者对该滑坡的成因给出了不同的解释：部分学者认为，当地的饱和砂土地基在剪切波作用下发生了液化，并在重力作用下发生了流滑大变形[1,2]，但大量振动台离心模型试验证明[3]，仅输入剪切波时，产生的剪应变不足0.15；也有许多学者认为，当地运河、灌溉渠以及地下水的影响加剧了流滑的发展[25,26,27]，但Balaroa滑坡并不受上述因素影响[28]。

Balaroa地区地震前后的卫星图像如图6所示，其中黄色多边形表示流滑变形较大的区域，根据当地的地质勘探结果[2]，主要为1～7 m厚的可液化饱和粉细砂，本构参数按表1的软弱土层取值，其余地表土则为不易弱化或液化的粉质黏土，本构参数按表1的坚硬土层取值。Li等[4]曾使用有限元法定性揭示了面波在该滑坡中所起的重要作用，但受制于变形规模，并不能完全再现流滑大变形现象。因此，本节以Balaroa滑坡为研究对象，使用SPH，通过对比当地在不同地震波输入下的变形现象，来解释流滑大变形的成因。

3.2基于GIS的建模方法

基于以上卫星图像和物探结果，可以结合地理信息系统（GIS）构建三维SPH数值计算模型，具体操作流程如下：

(1）根据地理空间数据云Aster GDEM高程数据库绘制地表面，参照卫星图像确定软弱土与坚硬土的水平分布范围，并根据当地物探结果确定两种土的竖直分布范围，通过三次样条曲线插值的方法确定分界面的形状，并使用Arc GIS软件构建三维几何模型，地基土的本构参数取值参照第2节表1。

图6 Balaroa滑坡的卫星图像  

(2）使用Dual SPHysics软件读取上述几何模型，将该模型离散化为约1400万个流体粒子，并添加约200万个边界粒子，从而构建如图7所示的SPH模型。该模型在北京超级云计算中心平台上使用GPU加速以及并行计算技术进行动力分析。

图7 Balaroa滑坡的SPH模型 

(3）根据印度尼西亚气象气候和地球物理局（BMKG）及日本国际协力机构（JICA）实测得到的帕卢市当地实测地震波加速度时程[28]，使用消除多项式趋势法进行基线修正[24]，并根据弹性柱面波的能量密度解析解分解体波和面波[4]，从而得到如图8所示的边界粒子的位移时程。根据第2节的计算结果，地震波在HBP黏塑性流体中衰减得较快，可以忽略地震波散射的影响，因此，可在模型的左侧边界和底部边界上分别输入对应于图8中瑞利波和剪切波的给定位移时程，无需施加用以处理散射波的人工边界。根据当地目击者的观察，流滑在地震中持续了约1 min[2]，因此，在20 s的动力计算结束后继续计算40 s，使变形能够充分发展。

图8三维模型中输入的位移时程

Fig.8 Time histories of input displacement of 3D model

3.3计算结果分析

不同地震波输入下三维SPH模型计算得到的残余水平位移云图如图9所示。可以发现，剪切波几乎不会引起明显的变形，绝大部分区域的残余位移均在10 m以下，说明单纯输入剪切波不会对场地变形产生明显的影响，与平面应变模型得到的结论基本一致。

与剪切波相比，瑞利波作用下土体产生了较为明显的侧向变形，最大位移达到近100 m，但永久位移大于50 m的区域（即图9中的红色与黄色区域）比较窄，沿y方向的宽度仅在50～100 m左右，与图6的卫星图像中观察到的滑动范围不完全吻合。

同时输入两种波时，虽然最大位移仍然在100 m左右，与瑞利波输入的工况接近，但发生较大规模永久变形的范围扩大了不少。复合波输入下，永久位移大于50 m的区域沿y方向的宽度最大可以达到300 m以上，且该区域的轮廓与卫星图像中滑裂面的形状和位置都比较接近，主要分布在滑坡的西北侧。当然，虽然复合波输入下该SPH模型出现了近百米的水平位移，但Balaroa滑坡实际观测到的最大位移约有300m[2]，说明瑞利波只是诱发该缓倾场地发生流滑大变形的成因之一，当地的水文地质条件以及工农业活动也对流滑破坏有重要影响。

为了更清晰地反映不同地震波输入下Balaroa滑动的破坏形式，在y=325 m处（见图6）沿铅直方向截取了一个截面，并绘制了静力变形以及不同输入工况下的动力残余变形形态与位移云图，如图10所示。

可以发现，剪切波输入下几乎不会产生残余变形，而在瑞利波和复合波输入下，都会产生最大近百米的残余水平位移。但纯瑞利波输入下当x>180 m后，残余位移就衰减至不足50 m；而复合波输入下直到x>500 m时，残余位移才下降至50 m以下，说明此时面波的影响范围远大于纯瑞利波输入的工况。瑞利波和复合波输入下残余变形在x<100 m时是逐渐增大的，这是因为x<50 m时地面几乎水平，随后坡度逐渐变陡，变形才会逐渐增加；之后随着面波影响的衰减，变形又会逐渐减小。

除了位移云图外，根据震后截面的形态也可以得到一些有用的结论。剪切波作用下倾斜角沿x方向的分布比较均匀，保持在2°～3°左右。在瑞利波输入下，截面的坡度会变得不再均匀，而在复合波输入下，变形后的截面呈现出明显的三段：当x<150 m时发生非均匀滑动，坡度较陡；当150 m≤x<630 m时发生整体流动，坡度平缓且比较均匀；当630 m≤x<830 m时，滑下的土堆积在滑动体底部，使地面抬高了6～8 m，坡度较陡，也是建筑物碰撞与破坏的主要区域。上述模拟结果与现场调查比较一致[2]。

综上所述，三维SPH模型的计算结果表明，对于缓倾场地，仅输入剪切波并不会产生明显的变形；瑞利波输入下则会在局部产生较大的变形，但影响范围有限；如果同时输入剪切波和瑞利波，不但产生的变形会更大，影响范围也要宽得多，远远不是剪切波和瑞利波输入下计算结果的直接叠加。

图9不同地震波输入下Balaroa滑坡的残余水平位移云图  

图1 0不同地震波输入下代表性截面的残余水平位移云图  

Fig.10 Nephogram of residual horizontal displacement of chosen profile under different seismic wave inputs

4、结论

为了研究瑞利波对缓倾场地流滑大变形的影响，并解释2018年苏拉威西地震Balaroa滑坡的成因，本文使用SPH构建了二维与三维数值计算模型，通过Dual SPHysics开源平台对比分析了剪切波与瑞利波作用下残余变形的特征，得到以下3点结论。

(1)SPH可以有效地模拟实际工程中的软弱地基大变形问题：在建模时通过GIS技术可以便捷地构建真实场地的三维模型，在计算中通过GPU加速技术和并行计算技术可以高效地处理大规模场地动力计算问题，在边界上使用DBP可以实现瑞利波的非一致输入，在内部使用HBP流变学模型可以模拟地基土弱化后的应力应变关系。

(2）对于简化的二维平面应变模型，不论是水平场地还是缓倾场地，剪切波输入下变形都比较小，而瑞利波输入下则会出现数十米的残余变形。

(3）对于模拟Balaroa滑坡的三维缓倾场地模型，剪切波输入下不会发生明显的变形，瑞利波输入下场地顶部附近会发生近百米的大变形，但范围较小。而当同时输入剪切波与瑞利波时，发生流滑大变形的范围较大，与实际震害现象比较吻合，证明了瑞利波等面波在边坡变形中可能发挥重要作用，在抗震分析中有必要考虑这类面波的影响。

当然，HBP流变学模型并不足以完整地描述饱和砂土地基液化前后的动力学行为，还需采用精度更高的本构模型与数值计算方法展开进一步的研究。本文主要针对Balaroa滑坡进行计算分析，由于瑞利波兼具行波效应和沿竖向分布不均匀的特征，为了进一步推广本文的研究结论，还需讨论可液化地基的场地规模对于变形特征的影响。

文章来源：李云屹,王睿,张建民.瑞利波作用下缓倾场地流滑大变形的SPH数值模拟[J].岩土工程学报,2023,45(07):1333-1340.