改进的SPH边界处理方法与土体大变形模拟

，

(1.中交第二航务工程局有限公司，武汉 430040；2.中交公路长大桥建设国家工程研究中心有限公司，北京 250073;3.上海交通大学 船舶海洋与建筑工程学院，上海 200240)

1 引 言

人工堆积的尾矿坝和粉砂质边坡等土工构筑物极容易发生滑移失稳等灾害，往往会造成极大的人身和财产损失。准确模拟这类土体大变形问题，对工程设计和灾害防御具有重大意义。土力学中常见的数值模拟方法如FEM(有限元法)和DEM(离散元法)在求解这类问题时都存在缺陷。FEM常因网格束缚导致计算结果不准确；DEM则由于选取计算参数困难，较难得出精确的数值模拟结果[1]。光滑粒子流体动力学SPH(Smoothed Particle Hydrodynamics)方法是一种拉格朗日型无网格粒子方法，起初主要用于求解天体物理学问题[2,3]。SPH方法将问题域离散成有限个携带物质信息的粒子，每个粒子都在周围粒子影响下遵守守恒定律运动，十分适用于求解大变形问题。目前SPH已广泛应用于流体力学及固体力学等各个领域，成功模拟了爆炸、穿甲及溃坝等大变形问题。近年来，越来越多的学者开始应用SPH方法模拟土体大变形问题。Bui等[4]将SPH方法应用在土体大变形问题中，采用Drucker-Prager屈服模型，模拟了不同内摩擦角的土坡滑移现象。黄雨等[5]基于等效牛顿黏度系数的概念，采用Bingham模型描述土体大变形，模拟了土体流动。Chen等[6]将各向异性的Drucker-Prager模型与修正Kondner和Zelasko模型应用到SPH中，用以分析土体的非线性特性及边坡的动力反应和滑移特征。Wang等[7,8]基于混合物原理，建立了两相流模型，模拟了射流冲刷和水下滑坡问题。

在SPH方法中，粒子的场函数由支持域内其他粒子场函数的加权平均而得，其中权函数是具有偶函数性质的光滑函数。因此，若想获得较高的计算精度，需支持域内具有足够多的粒子，且粒子整体呈对称分布。但这两条性质在边界附近都很难满足[9]。该边界附近粒子缺失的问题，需要通过边界处理方法来解决。常用的边界处理方法包括排斥力法和虚粒子法，排斥力法易受初始扰动的影响[10]，而虚粒子法则难以处理复杂边界[11]，且存在边界零粒子层问题[12]。近年来，一种新兴的边界处理方法正引起学者们的注意。Ferrand等[13]提出了能够应用于任意形状边界的边界处理方法，即统一半解析壁面边界条件USAW(unified semi-analytical wall boundary conditions)处理方法，成功模拟了溃坝流问题。Leroy等[14]改进了USAW边界处理方法，并将其与ISPH方法结合。Mayrhofer等[15]应用USAW方法成功模拟了三维溃坝流问题。Cercos-Pita等[16]提供了一套基于GPU并行计算的USAW-SPH开源代码。文献[17]基于USAW方法，提出了两种不规则粒子分布，更好地模拟了复杂边界溃坝问题。USAW边界处理方法在流体动力学上的成功，给土体大变形问题的边界处理提供了新思路。本文首次应用USAW边界处理方法模拟土体大变形问题。

以往USAW方法中，研究者一般生成规则的初始粒子分布，其中内部粒子质量相同，而边界粒子质量取决于该点在边界处所形成的角度[13]。但在这种方法中，边界粒子的质量会影响问题域的模拟，并且会出现边界零粒子层，影响模拟精度。结合USAW基本理论，本文提出了边界粒子无体积和无质量的做法。

2 土体大变形理论模型

土体的控制方程由质量守恒和动量守恒方程构成：

(1，2)

定义土体为理想弹塑性材料，且服从Drucker-Prager屈服准则。在该模型中，屈服函数可表示为

(3)

(4，5)

本文算例是大变形问题，选取与刚体旋转无关的Jaumann应力率,

(6)

(7)

弹塑性材料可将材料的应力应变关系分解为弹性部分和塑性部分。弹性部分采用胡克定律求解，塑性部分则遵守塑性流动法则，最后获得本构方程[7]：

(8)

(9)

(10)

式中Ψ为塑性势函数中的剪胀角，参考文献[4]取值为0，本文选用的塑性势函数g为

(11)

3 SPH理论和边界处理方法

在SPH方法中，将问题域离散成多个粒子，每个粒子具有体积和质量等属性，并携带应力、速度和位置等信息，能够在力的作用下运动。粒子a的场函数fa(如压强或速度)及其导函数fa，均可由其支持域内各粒子的场函数插值近似得到，

(12)

Wa b

(13)

式中Vb为粒子b的体积，P为支持域内所有粒子点的集合，本文符号较多，为避免混淆，将其列入表1。在SPH方法中，目标粒子的场函数通过周围粒子对应函数的加权平均求得，见式(12)。式(13)是SPH方法最基础的导函数近似式，该式将导函数的加权平均转换成了核函数导数的加权平均，方便了导函数的求解。Wa b=W(xa b,h)为光滑核函数，是SPH方法加权近似的权函数，具有紧支性和归一性，其大小随粒子a与粒子b之间相对距离的增大而减小；xa b=xa-xb，xa和xb分别指粒子a和b的位置矢量，h为光滑长度。本文选取Wendland五阶核函数：

(14)

式中Ra b=xa b/h；αd为保证光滑核函数归一性的量，在二维算例中，其值为7/(4πh2)。更详细的SPH理论推导可参考文献[7]。

表1 符号含义

Tab.1 Symbol definition

符号 含义a,b粒子(含边界粒子)e仅边界粒子s边界上两个粒子之间的线段P支持域所有粒子的集合Ω整个计算域F内部粒子的集合(除边界)S边界上所有线段的集合

严格来讲，式(12,13)只适用于远离边界处的粒子，当粒子靠近边界时，由于边界的截断，粒子的支持域内将缺少部分粒子。为了弥补边界粒子缺失，USAW边界处理方法在插值近似时引入了修正因子(renormalization factor)：

(15)

式中γa为粒子a的修正因子，定义为

(16)

式中Ω为粒子a的支持域，W(xa-x,h)为粒子a的光滑函数。当粒子a远离壁面时，支持域完整，γa=1；靠近壁面时，由于边界的截断，支持域破缺，导致γa<1。γa可采用时间步进法[13]或解析方法[17]求解，本文采用后者。在计算函数导数时，与式(13)不同，USAW边界处理方法保留了边界积分项，

(17)

(18)

式中l为线段s的长度。

边界粒子的密度可通过SPH插值近似得到，

(19)

(20)

而粒子a处速度的偏导，可由式(21)离散得到。

(21)

式(21)采用常用的SPH变化技巧[7]，具有一定的反对称性，符合作用力与反作用力原理。粒子a上其他场函数的偏导数可做类似处理。

传统的SPH方法大多假定边界粒子与内部粒子属性一致，使用相同的方法来求解密度和压强；在USAW方法中，则采取不同的方法。

内部粒子的密度，通过修正后的密度求和法[13]来求解：

(22)

式中 上标n指的是时间步。式(22)是运用传统的密度求和法，将两个相邻时间步的密度求解方程相减得到[13]。同时，为了减少数值振荡，本文在方程中添加了δ-SPH 技术[18]，即式中最后一项，其中，δ为控制修正大小的一个值，本文取0.1。

为了减少数值振荡，除上文介绍过的δ-SPH 技术外，本文还采用了Marrone等[19]提出的人工粘性技术。其不仅能够将动能转变为热能，同时也提供了冲击的耗散，避免了粒子互相靠近时的穿透现象。在模拟过程中发现，固体粒子在滑坡过程中容易出现聚集成块的现象，即所谓的张力不稳定现象。为了减小这种张力不稳定，本文采取了由Monaghan[20]提出的人工应力法。在数值模拟中，有时应力会跃出屈服面外，这是不允许的。本文采用张裂处理(Tension crack treatment)和比例拉回(scaling back procedure)的办法将屈服面外的应力映射回屈服面上[4]，并采用二阶精度的蛙跳法(Leap -Frog)进行时间步进。具体可参见文献 [7,8]。

图1 USAW中的边界定义

Fig.1 Boundary definition in USAW method

关于SPH方法的粒子分布，一般先将问题域划分为若干个相邻的矩形体积单元，再取矩形单元的中心位置为粒子位置，粒子的体积等于对应单元的体积。除问题域外，还需布置若干粒子来模拟边界线，这些粒子是位于问题域外的虚粒子。虚粒子虽相对于边界静止，但同样参与计算，如图2(a)所示。图中黑色实心点为内部粒子，其体积之和为模拟的问题域，灰色粒子点为边界粒子，虚线为问题域的边界。该方法由于边界粒子无法与壁面贴合，较难模拟复杂壁面形状。

USAW边界处理方法可以处理任意形状边界，部分原因是该方法直接假定边界粒子的连线作为壁面，如图2(b)所示。图中边界粒子仅在问题域内部分上色，是为了反应传统USAW方法中边界粒子的质量分布[13]，即边界处粒子质量一般是内部粒子的一半，直角拐点处粒子质量是内部粒子质量的1/4。在这种分布方法中，边界粒子体积占据了一定内部问题域的体积，即相当于在边界处固定了一排静止不动的内部粒子，从而导致问题域小于需模拟的尺寸，是不合理的。

本文发现USAW边界处理方法中，边界粒子可以不具有体积和质量。USAW方法在求解导函数时，将积分项分成了体积积分项和边界上的面积分项，在边界的面积分项中，仅边界粒子的密度参与了计算，见式(17)。同时在边界处密度和压强的求解式(25,26)中，边界粒子的质量或体积都不参与计算。故在边界处采用这种无体积和无质量的边界粒子十分合理。

图2 有质量边界粒子的分布

Fig.2 Distribution of boundary particles with mass in different methods

本文提出边界粒子无体积和无质量的处理方法，如图3所示，其中空心粒子为边界粒子。USAW边界处理方法将场函数导数分为体积分和面积分的累加，而这种边界粒子无质量的分布方法，也对应地将粒子区分为有体积和无体积两部分，具有更高的精度，也更符合USAW边界处理方法的本质。使用无体积和无质量的边界粒子，既能保证边界粒子与壁面完美贴合，正确模拟复杂形状，又能准确地求解导函数中的面积分。

4 数值计算结果分析

模拟沙土滑移，基本模型如图4所示，在一个矩形槽左下角，存在一个矩形沙柱，初始时刻沙粒静止，当快速拔除右边的挡沙板之后，沙柱在重力作用下垮塌。其中灰色矩形是初始时刻的沙土，hi和di分别为初始高和宽，高宽比为k；虚线是滑坡终止时的沙面形状，其中h∞和d∞分别为滑坡结束时沙土的最大高度和最大宽度，δd为滑坡终止时滑移的最远距离。

本文模拟不同初始高宽比k下的滑坡运动，分析初始高宽比对滑移最远距离和最终沙粒高度的影响，并与Lube等[21]测量的实验值进行比较，验证算例的可靠性。在SPH模拟中，采用与实验所用的粗石英砂相同的参数，即内摩擦角为31°，粘聚力为0。长宽比k分别为0.5，1和3，初始长度di恒为0.2 m，各算例初始粒子间距均为0.002 m，光滑长度均取初始粒子间距的1.2倍，3个模型的总粒子数分别为5350，10650和31950，计算时间步长为5×10-5s，计算结果如图5所示。

图3 无质量边界粒子分布

Fig.3 Distribution of boundary particles without mass

图4 溃沙模型和基本符号

Fig.4 Landslide model and basic symbols

3个算例的滑移最远距离d∞分别为0.35 m，0.51 m和1.13 m，最终沙粒高度h∞分别为0.1 m，0.19 m和0.29 m。Lube等[21]分析了滑坡的实验数据，总结出滑移距离δd与初始宽度di的无量纲比值ζ(k)可表示为

(23)

式中c1=1.6,c2=2.2；当k位于1.8～2.8之间时，没有准确函数可表达。滑移后左端高度h∞与初始宽度di比值ξ(k)可表示为

(24)

数值计算值和实验值的比较结果列入表2，相对误差计算方法为实验值与模拟值之差的绝对值除以实验值。由表2可知，模拟结果与实验值的ζ(k)和ξ(k)相对误差较小，在5%左右，属工程误差范围之内，故可认为本文方法能够较好地模拟溃沙过程。

图5 不同高宽比下的溃沙终止形态

Fig.5 Landslide simulations with different aspect ratios

从图6(a)可以看出，内摩擦角增大时，滑移距离减小。从图6(b)可以看出，黏聚力增大时，滑移距离减小；黏聚力继续增大时，沙体的表面将变得粗糙，沙柱右上角在滑动过程中在沙面上形成一个尖点，滑动终止时也无法消失，且黏聚力越大，最后保留的尖点也越大。

SPH方法模拟运动学问题时，常出现边界零粒子层现象[12]，在运动介质的舌尖与边界之间，经常会出现一个大约为1个粒子大小的空隙。该零粒子层会影响边界层内的数值模拟精度，而且粒子数越少，单个粒子所占体积越大，零粒子层就越厚。本文采用三种边界处理方法，模拟了图5(a)的溃沙问题，并观察其舌尖与边界之间是否存在零粒子层现象，如图7所示。

表2 模拟值与实验值之间的相对误差

Tab.2 Relative error between simulated and experimental values

Kζ(k)模拟值ζ(k)实验值相对误差/%ξ(k)模拟值ξ(k)实验值相对误差%0.50.750.86.250.50.5011.551.63.130.951534.654.61.11.351.556.4

图6 不同参数下的模拟结果

Fig.6 Simulation results under different soil parameters

可以看出，使用虚粒子法模拟时，粒子舌尖与边界之间明显存在零粒子层现象；在传统的USAW边界处理方法中，依然可以观察到边界零粒子层的存在，但较虚粒子法稍薄，主要原因是其边界粒子质量仅为内部粒子质量的1/2；图7(c)显示，本文提出的边界粒子无质量的USAW方法有效消除了边界零粒子层现象，而且其舌尖位置更加接近实验位置。

探索滑坡冲击楔形体的运动过程，并分析滑坡对楔形体的冲击力。如图8(a)所示，在问题域的左下角，存在矩形沙土，沙土高为0.6 m，宽为0.2 m，底部距原点0.3 m处存在楔形体凸起，楔形体呈等边直角三角形形状，斜边长为0.18 m，方向朝下。SPH模拟的算例共由8400个粒子组成，其余参数均与上文保持一致。

滑动前0.8 s内的沙体轮廓如图8(a)所示。滑坡开始0.14 s后，滑坡才抵达楔形体底部，随后越来越多沙土覆盖楔形体左侧，楔形体左侧压力相应急剧上升；t=0.22 s时，沙土刚好全部覆盖楔形体；此后，部分沙土冲过楔形体，形成一个勾形舌尖；t=0.34 s时，舌尖落入地面，与楔形体形成一个空穴，并随后填满；然后沙土继续向前滑行，速度不断减慢，约进行到0.8 s时，停止滑行，最远距离达到1.1 m，整个过程中有少数粒子脱离舌尖，以较大速度抛物前行。

图7 边界零粒子层

Fig.7 Zero particle layer near boundary

图8 楔形体溃沙轮廓图

Fig.8 Profiles of landslide on a wedge

沙土冲击建筑物所造成的压力是工程设计关注的重点。该算例中，楔形体左侧压力随时间变化如图9所示。可以看出，滑坡开始0.14 s后，沙土接触楔形体并对楔形体施压，楔形体左侧压力对应急剧上升；0.22 s时，沙土刚好全部覆盖楔形体，楔形体左侧所受压力达到最大值，约为15.5 kN；此后，沙土冲过楔形体，对应的楔形体左侧压力急剧减小，在0.4 s时，沙粒与楔形体右侧壁面相互作用，从而导致左侧壁面压强形成一个小波峰；然后，楔形体左侧沙土减少，且沙土流过楔形体的速度逐渐降低，左侧压强缓慢减小，最终静止时压力值约为3 kN。边界压力的强烈振荡一直是限制SPH方法发展的重要阻碍，本文方法模拟的边界压力振荡较小。

图9 楔形体左侧压力随时间变化情况

Fig.9 Time evolution of the pressure on the left ide of the wedge

5 结 论

SPH方法是拉格朗日型无网格粒子法，十分适用于处理大变形问题。本文应用SPH方法，结合改进后的USAW边界处理方法，较精确地模拟了滑坡问题，具体结论如下。

(1) 探讨并改进了USAW边界处理方法，首次使用梯形法求解γa s，并将传统的SPH数值技术如 δ-SPH技术、人工粘性以及人工应力与USAW边界处理方法相结合，改善了USAW方法的计算效果。

(2) 提出了边界粒子无质量的处理方法，该方法适用于模拟复杂边界问题，并提高了USAW边界处理方法的精度，还能避免边界零粒子层现象。

(3) 首次应用USAW边界处理方法模拟滑坡问题，证明了该方法模拟土体大变形问题的可行性，分析了内摩擦角和黏聚力等土体物理特性参数对滑坡过程的影响，并计算模拟了滑坡冲击楔形体过程，研究了楔形体所受压力的变化规律。

本文将改进的SPH方法应用于滑坡问题计算，通过比较计算结果与实验结果，验证了方法的有效性。但对于土力学中大量的其他问题并没有涉及，这将是今后的研究内容之一。

:

[1] 潘建平，曾庆筠.SPH方法在土工大变形分析中的应用研究进展[J].安全与环境学报，2015，15(5):144-150.(PAN Jian-ping，ZENG Qing-yun.Research review over the advances of the application of SPH method to the analysis of the serious geotechnical deformation[J].JournalofSafetyandEnvironment，2015，15(5):144-150.(in Chinese))

[2] Lucy L B.A numerical approach to the testing of the fission hypothesis[J].TheAstronomicalJournal，1977，82:1013-1024.

[3] Gingold R A，Monaghan J J.Smoothed particle hydrodynamics:theory and application to non-spherical stars[J].MonthlyNoticesoftheRoyalAstronomicalSociety，1977，181(3):375-389.

[4] 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].InternationalJournalforNumericalandAnalyticalMethodsinGeomechanics，2008，32(12):1537-1570.

[5] 黄 雨，郝 亮，谢 攀，等.土体流动大变形的 SPH 数值模拟[J].岩土工程学报，2009，31(10):1520-1524.(HUANG Yu，HAO Liang，XIE Pan，et al.Numerical simulation of large deformation of soil flow based on SPH method[J].ChineseJournalofGeotechnicalEngineering，2009，31(10):1520-1524.(in Chinese))

[6] Chen W，Qiu T.Simulation of earthquake -induced slope deformation using SPH method[J].InternationalJournalforNumericalandAnalyticalMe-thodsinGeomechanics，2014，38(3):297-330.

[7] Wang C，Wang Y，Peng C，et al.Smoothed particle hydrodynamics simulation of water-soil mixture flows[J].JournalofHydraulicEngineering，2016，142(10):04016032.

[8] Wang C，Wang Y，Peng C，et al.Two -fluid smoothed particle hydrodynamics simulation of submerged granular column collapse[J].MechanicsResearchCommunications，2017，79:15-23.

[9] 周 杰，徐胜利.SPH方法的运动边界条件研究[J].计算力学学报，2016，33(3):412-417.(ZHOU Jie，XU Sheng-li.Research on the moving boundary condition in SPH methods[J].ChineseJournalofComputationalMechanics，2016，33(3):412-417. (in Chinese))

[10] 上官子柠，马庆位，韩端锋.SPH 边界处理中两种排斥力模型的对比分析[J].船舶力学，2014，18(1):37-44.(SHANGGUAN Zi-ning，MA Qing-wei，HAN Duan-feng.Comparisons of two repulsive models for boundary treatment in SPH [J].JournalofShipMechanics，2014，18(1):37-44.(in Chinese))

[11] 艾孜海尔·哈力克，热合买提江·依明，开依沙尔·热合曼，等.不同质量粒子分布 SPH 方法及其应用[J].振动与冲击，2015，34(22):62-67.(AZHAR Halik，RAHMATJAN Imin，KAYSAR Rahman，et al.SPH method with different mass particle distributions and its application[J].JournalofVibrationandShock，2015，34(22):62-67.(in Chinese))

[12] Ni X Y，Feng W B.Numerical simulation of wave overtopping based on DualSPHysics [J].AppliedMechanicsandMaterials，2013，405:1463-1471.

[13] Ferrand M，Laurence D R，Rogers B D，et al.Unified semi-analytical wall boundary conditions for inviscid，laminar or turbulent flows in the meshless SPH method[J].InternationalJournalforNumericalMethodsinFluids，2013，71(4):446-472.

[14] Leroy A，Violeau D，Ferrand M，et al.Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[J].JournalofComputationalPhysics，2014，261:106-129.

[15] Mayrhofer A，Ferrand M，Kassiotis C，et al.Unified semi-analytical wall boundary conditions in SPH:analytical extension to 3-D[J].NumericalAlgorithms，2015，68(1):15-34.

[16] Cercos-Pita J L.AQUAgpusph，a new free 3D SPH solver accelerated with OpenCL[J].ComputerPhy-sicsCommunications，2015，192:295-312.

[17] 骆 钊，汪 淳.统一半解析边界处理SPH方法与不规则粒子分布[J].水动力学研究与进展A辑，2017，32(2):189-197.(LUO Zhao，WANG Chun.Unified semi-analytical wall boundary treatment in SPH and irregular particle distribution[J].ChineseJournalofHydrodynamics，2017，32(2):189-197.(in Chinese))

[18] Leroy A，Violeau D，Ferrand M，et al.Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[J].JournalofComputationalPhysics，2014，261:106-129.

[19] Marrone S，Antuono M，Colagrossi A，et al.δ-SPH model for simulating violent impact flows[J].ComputerMethodsinAppliedMechanicsandEnginee-ring，2011，200(13):1526-1542.

[20] Monaghan J J.SPH without a tensile instability[J].JournalofComputationalPhysics，2000，159(2):290-311.

[21] Lube G，Huppert H E，Sparks R S J，et al.Collapses of two -dimensional granular columns[J].PhysicalReviewE，2005，72(4):041301.