颗粒离散元模拟边坡动力响应的研究进展

常文斌，王 平，2，于一帆，王会娟

(1.中国地震局兰州地震研究所 黄土地震工程重点实验室， 甘肃 兰州 730070；2.西安理工大学 岩土工程研究所， 陕西 西安 710048)

边坡是一个复杂多变的系统，其稳定性受多种条件的影响，如自然条件、边坡结构、人类活动等，这些影响因素大多是随机的、可变的，其不确定性增加了边坡稳定性的研究难度。

边坡在动荷载作用下的稳定性分析较为复杂，它包含了工程地震学和岩体力学等多学科的交叉以及诸多技术问题。动荷载包括地震荷载、风荷载等自然荷载及机械荷载、爆炸荷载等人为荷载。与其他动荷载不同的是，地震荷载并非一种直接的力学作用，而是地面发生位移时坡体结构在惯性作用下产生了一定的运动差，导致边坡产生严重的变形。同时，地震荷载的频谱特性及传播方向对坡体的破坏均有很大的影响，加之需要考虑岩土体的动本构关系及动强度，导致边坡动力稳定性问题变得更加复杂。然而随着计算机科学技术的发展，数值模拟技术已成为边坡动力分析领域不可或缺的手段之一。当下主要的数值方法包括有限差分法、有限元法、离散元法[1]及边界元法等。现阶段，有限差分法及有限元法是边坡动力分析常用的两种方法，由于这两种方法将坡体视为连续介质，因而在动力分析方面有着独特的优势。但将岩土体视为连续介质则不能真实的体现出介质材料的流动变形特征，造成了有限差分法及有限元法在模拟边坡变形时的局限性。

1971年美国学者 Cundall[2]提出了离散单元法。随后Cundall等[3-4]首次把颗粒流法定为一种特殊的离散单元法。在此基础上，美国Itasca 公司开发出了二维颗粒流软件PFC2D(Particle Flow Code in 2Dimensions, PFC2D)[5]和三维颗粒流软件PFC3D(Particle Flow Code in 3Dimensions, PFC3D)。PFC中，球和簇的运动遵守牛顿第二运动定律，墙体的运动通过用户来指定。颗粒之间的接触可以选用不同相互作用定律的接触模型来施加，在受力体系下，通过颗粒间内力和位移的不断更新来模拟介质体的变形[6]。颗粒流方法很好的解决了常规数值模拟方法在分析大变形岩土工程破坏时的局限性，能够更好的揭示微观及细观介质材料的破坏机理，因此备受广大学者青睐。

颗粒流离散元法解决了岩土工程大变形及介质的微观破坏问题，可更方便的处理非连续介质力学问题，但是对于边坡的动力响应研究由于启动较晚，目前仍有诸多问题需要解决。本文总结归纳了颗粒流进行边坡动力响应研究时，阻尼参数标定、动荷载输入方式和动力分析参量这三个方面的研究现状及存在的问题。分析过程中基于地震动力作用的复杂及特殊性，着重对离散元颗粒流方法在边坡地震动力响应方面存在的关键问题进行了分析和探讨。

1 颗粒流中阻尼参数标定

1.1 PFC程序能量耗散机制

PFC程序中，单元之间的能量是通过摩擦滑动来进行耗散的。由于单元之间的滑动摩擦并非主动发生，因此在所需的计算时步内，仅靠摩擦滑动系统很难达到一个稳定的状态。鉴于此，PFC中提供了3种阻尼耗能机制：局部阻尼、组合阻尼、黏性阻尼。其中，局部阻尼通过对每一个颗粒单元施加一个与其速度方向相反的阻尼力。通过局部阻尼理论推导部分可知，施加于颗粒上的阻尼力于颗粒的不平衡力成正比，两者之间的比值即为局部阻尼系数。PFC中局部阻尼力的施加则通过设置局部阻尼系数来实现，默认值为0.7。组合阻尼同时考虑了颗粒的速度与不平衡力，但其能量耗散率较低，一般不被考虑使用。黏性阻尼施加于PFC模型中每一个接触上，包括法向阻尼和切向阻尼。黏性阻尼的方向与接触处的相对速度方向相反，其定义通过指定黏性阻尼比参数来实现。

许多学者对PFC中阻尼参数的选取做了相关研究[7-9]，结果表明：在静力情况下，PFC中采用局部阻尼以吸收单元之间的动能，其阻尼取值仅影响系统达到稳定解的时间。而在动力问题中，应力波在岩土体的传播过程中，阻尼取值直接影响着岩土体的振动，从而直接影响结算结果，因此必须确保阻尼取值的科学性。当下PFC在岩土工程中的模拟研究多集中于颗粒细观参数等对模拟结果的影响，阻尼系数取值的重要性被多数学者忽略，无形中影响了模拟结果的准确性。

1.2 颗粒阻尼参数的选取

颗粒流模拟过程中，阻尼参数的选取往往存在主观经验性及随意性。为解决这一问题，钟文镇等[7]采用声频采样技术对钢球自由落体物理模型中的阻尼系数进行了标定，得出的阻尼系数(0.5左右)与实际较为相符。同时，通过颗粒在容器中随机堆积的模拟试验验证了所标定系数的可靠性。杨冰[8]在PFC2D中采取不同的局部阻尼系数及黏性阻尼系数进行了颗粒的自由落体试验，计算不同阻尼系数所对应的能量衰减率，拟合出了颗粒自由落体时能量衰减率与阻尼系数的关系，验证了颗粒材料通过设置相应的阻尼系数也可用于工程的模拟。周国庆等[9]通过钢球自由下落和悬臂梁的剪弯受力模拟试验，明确了颗粒流中局部阻尼系数和黏性阻尼系数的适用范围。结果表明静态问题可采用局部阻尼减少计算时间，动态问题需根据不同情况选取局部阻尼或黏性阻尼，并通过钢球的自由落体试验确定了钢球的法向黏性阻尼比为0.16。

由此可见，PFC中阻尼参数的选取对动力问题计算结果影响较大，在进行边坡动力分析模拟工作时更需要注意。以地震作用为例，地震作用下，颗粒集合体边坡发生滑动后，多数颗粒为单向运动，此时阻尼类型的选取尤为关键。如果对颗粒设置局部阻尼，由于局部阻尼是对颗粒速度相反方向施加一个阻尼力，只要颗粒存在速度就会产生作用。这会造成此后计算的每一时步都会产生与颗粒运动相反的阻尼力，导致颗粒的运动速度急剧减小，能量耗散迅速，从而影响最终的模拟结果。因此局部阻尼不适合地震作用下的颗粒流模型，部分学者在研究过程中很容易忽视这一点，对地震边坡模型采用局部阻尼且系数过大。相比之下，黏性阻尼比较适合地震边坡模型，它相当于在颗粒间每一个接触点上施加法向和切向阻尼器，阻尼力由系统的有限质量及自振频率等决定，并在计算循环时添加到接触力上。因此在颗粒流边坡动力分析过程中，根据不同的动力荷载及研究问题，选取合适的阻尼类型是确保模拟计算成功的前提，值得重视。

2 颗粒流中动荷载的输入方式

在离散元颗粒流中，岩土体介质材料通过颗粒及其之间的接触特性来表示。应力波在弹性状态的连续性介质中的传播过程，同样可以在颗粒流中进行模拟。

2.1 动边界条件

在颗粒流方法中，应力波可以通过对颗粒或者墙体施加来模拟动荷载。无论应力波施加于颗粒还是墙体，均要提及一个问题：应力波传播的动边界条件。

对于数值模型一般采用的单面边坡，在动力问题中，模型边界条件对动力分析结果有较大影响。分析模型的范围越大，分析结果就越贴近事实，但大的模型由于颗粒数目较多会导致计算量变大，这时候需要设置一定的人工边界。全局人工边界与局部人工边界是人工边界中最常用的两类边界。其中全局人工边界是一种较为精确的人工边界，其边界参数与相邻介质的特性相联系，藕联性较强，对于颗粒离散元这种没有控制方程的计算方法并不适用。局部人工边界主要有黏性边界[10]、黏弹性边界[11]、吸收层边界[12]与透射边界[13]。黏性边界、黏弹性边界以及透射边界在有限元及有限差分法中有大量的应用，以往用离散元颗粒流进行边界处理时大多参考这些数值方法的边界处理方法[14-16]。

颗粒流中动边界存在多种形式，其主要目的都在于吸收掉入射波的动能，防止其在边界处发生反射，从而达到模拟无限介质的效果。吸收层边界通过在模型外施加一层具有耗能作用的缓冲区吸收透射波(对边界吸收层颗粒施加一定的阻尼系数)，该吸收区又叫做阻尼层[12]。该边界最大的优点就是边界的几何形状可以任意设置，适用于任意入射角、各种频率的入射波。但是这种对模型边界颗粒施加阻尼系数的高阻尼吸波方法并非完美。石崇等[17-18]在对颗粒流中应力波的传播进行研究时，对黏滞透射边界和高阻尼吸收边界两种边界的合理性进行了验证。结果表明应力波传播过程中高阻尼边界容易导致边界附近的应力峰值明显降低，有悖于自然界中应力波的传播规律。而黏滞透射边界情况下虽然应力波传播距离略有延迟，但应力波的振动规律较为接近。周兴涛等[19-20]在二维颗粒流软件中分别建立了黏性、黏弹性及自由边界等人工动力模型，探讨了边界颗粒分布模式对边界作用的影响，验证了黏弹性边界的正确性。当下颗粒流边坡动力分析中，黏滞透射边界多被学者使用，其模拟结果也与振动台及离心机试验结果较为相符。

当下连续性介质方法与非连续介质方法耦合计算方法已被广泛使用，包括离散元与有限差分法[21-23]，离散元与有限元等[24-28]。不同数值方法的耦合取长补短，以PFC与FLAC耦合为例，李祥龙[29-30]采用颗粒离散元来模拟边坡的潜在破坏区域，用FLAC模拟边坡其他部位，进行了边坡的地震动力响应模拟。这种方法不仅解决了颗粒流模拟大型边坡计算量过大的问题，同时避免了离散元中动边界条件设置不成熟的问题。该方法对于离散元与其他数值方法耦合模拟边坡动力响应的研究提供了新的思路。

2.2 应力波的作用方式

应力波的作用方式包括应力波的作用方向(不同入射方向)和作用性质,同时应力波的频率强度、持时及类型对数值模拟结果均存在一定的影响。通常研究中施加的爆破荷载[31-33]、机械荷载[34-35]等为动力荷载在模拟时应力波的作用方式相对明确，影响较小。但对于类似地震波的不可控荷载，其作用方式则存在较多的争议。

目前，在边坡地震颗粒流值模拟分析中，地震动普遍采用一致性输入。实际中由于远近场的作用，地震动在经过不同路径传播到边坡底部时，边坡底部不同位置接收到的地震动是存在一定差异的，所以当下一致性输入地震动的方法是不完全正确的。由于目前边坡地震动力响应数值分析仍集中于有限元及有限差分法等方法，地震动作用方式的相关研究也都是基于这些方法展开。许多学者基于地震波的不同入射方向及不同震相的性质做了较多的研究[36-37]。国内以孙进忠教授为代表的团队在这方面做了相当多的研究工作[38-40]。结果表明不同的地震波作用方式对于边坡的动力响应的确存在着一定的影响。这些研究成果对于在颗粒流方法中进行应力波输入时有着很大的启发。

即便如此，采用地震台站记录的典型强震地震动数据对模型直接输入或进行简单调整后作为地震动输入(颗粒流中对颗粒或墙体施加速度)，虽然输入的是实际地震动，但由于岩土体参数、传播路径、震源参数等的差异，所得分析结果也只能进行边坡动力破坏机理的研究，并不能反映实际的动力响应情况[41]。因此，颗粒流中应力波(尤其地震波)的作用方式这是一个相当复杂问题，如何利用震源和岩土介质特性来确定适合于数值边坡动力分析的地震动类型，其入射方向、幅值、频谱和持时该如何确定，仍有待于更深层次的研究。

3 颗粒流中边坡动力分析方式

3.1 颗粒流中边坡动力分析对象

动荷载作用下，边坡动力响应分析的对象主要是速度、位移、加速度、应力等。在颗粒离散元边坡动力分析中，现有的文献大多都是对振动台模型试验、有限元、有限差分法等研究方法与研究结果的重复，多数学者关注点仍局限于边坡不同部位的速度、位移及加速度变化情况[42-46]，在此不再赘述。诚然，这些参量可以在一定程度上体现边坡的动力响应，但没有充分体现颗粒离散元的优越性。例如颗粒流中，介质的破裂能够定义岩土体损伤、力链可以用来考察能量传输等。颗粒流边坡模拟的过程中，可以对颗粒间的不同类型的接触力变化及坡体内裂纹的分布情况进行监测[47-48]，坡体中不同类型的裂隙(拉裂隙、剪裂隙)展布及裂隙数目可以清楚的展示出来。Wu等[49]采用PFC2D对含裂缝黄土边坡的地震作用下破坏过程进行模拟，边坡破坏过程中裂纹分布与振动台模拟结果类似。同时，颗粒流中每个计算步都会累积计算球、墙和接触的能量，诸如动能、摩擦能、应变能等能量；测量圆方法可以对小范围内颗粒的应力、孔隙率、应变率等进行监测。丁建辉[50]和苏燕等[51]采用二维颗粒流软件对地震降雨型滑坡形成机理进行了研究，通过对地震及降雨作用下坡体内的速度场、位移场、应力场、孔隙率、颗粒接触数及滑动比进行分析，揭示了地震降雨滑坡的形成机制。当下，将颗粒流的这些细观分析优势与边坡动力响应结合起来，而不是仅仅局限于速度、位移等传统的参数分析，有助于从不同角度分析边坡的破坏机理。

3.2 颗粒流中边坡加速度响应相关问题

在颗粒流中进行边坡动力响应分析时，由于不能直接监测坡体内测点处颗粒的加速度，目前多数学者采用的做法是将记录到的速度时程曲线进行微分，对得到的加速度时程曲线进行分析[15-16]。由于颗粒离散元的特殊性，在动荷载作用的过程中，边坡某一监测点的颗粒在小范围内处于一个不断运动的状态，这种运动状态受到其周围颗粒的碰撞和应力波的共同影响。因此在以往的研究中，通过微分速度时程曲线得出颗粒的加速度时程曲线与原始应力波曲线存在较大差异的现象多有发生。那么在颗粒流中基于这种方法所得加速度动力响应的可靠性就值得商榷。如何有效地解决或避免当下加速度响应分析存在偏差的这一问题，是颗粒流边坡加速度响应分析工作的关键。刘一飞[12]采用Fish语言编程实现对颗粒加速度的定义，即遍历颗粒的在某一瞬间的不平衡力，通过不平衡力除以颗粒质量得到颗粒的加速度，但其研究结果表明通过不平衡为计算的加速度分布有很大的随机性和离散性，其峰值加速度(PGA)并不准确，参考意义并不是很大。笔者建议尝试在颗粒流建模过程中，首先尽量设置颗粒的粒径较小，以减小颗粒在小范围内的运动。其次，在保证某一监测点处局部颗粒的总体动力响应情况一致的前提下，通过对监测点小范围内的颗粒的监测数据求取平均值，以减少偶然性误差。

在此基础上，应力波频谱分析在动力响应中的应用较少，并不是很成熟。地震频谱可以定量的揭示地震动的动力特性，以不同频率分量的表现来表现地震动及其对结构的作用。傅里叶频谱[52]和反应频谱[53]是地震波频谱特性分析时常用的频谱，其中，傅里叶频谱反应了地震波在坡体内部传播过程中自身频率各分量的变化情况，反应频谱表征了地震动对上部自由弹性结构的影响特征。将频谱特征响应引入颗粒流边坡动力分析中，结合颗粒流边坡的宏观破坏机制，分析结果对于边坡的防护治理及坡体上部临坡建筑物的抗震设计意义重大。然而，颗粒流中边坡频谱响应特征的分析需要基于监测点处加速度时程曲线提取方法的成熟，需要进一步的探索。

4 结论及展望

当前离散元方法已成为研究边坡变形破坏的主要手段之一，开展在颗粒离散元基础上的边坡动力响应研究，对于现有的边坡动力分析研究成果是一个极大的补充。本文从阻尼参数标定、动荷载输入方式与动力分析三个方面，对离散元颗粒流在边坡动力响应领域的研究进展进行了总结与分析。目前值得注意的关键问题如下：

(1)颗粒流边坡动力破坏模拟过程中阻尼类型的选取及阻尼系数的标定应引起重视。离散元模型中，组合阻尼由于能量耗散率较低，一般不建议使用；局部阻尼适用于静力作用下吸收颗粒单元之间的动能减少系统达到稳定的时间，而在动力作用下对能量的耗散较大，因此不适应动力分析时使用；黏性阻尼的阻尼力在振动过程中是不断循环添加至单元上的，因此适用于动力响应分析时使用。因此在模拟过程中，应根据不同情况选取对应的阻尼类型及阻尼系数，以提高计算的准确性。

(2)如何结合岩土体材料特性及震源机制等对边坡模型合理的输入对应的地震波。在模拟边坡动力破坏的过程中，地震动参数的选取往往至关重要。不同三要素(幅值、频谱、持时)的地震动数据对模型计算结果的影响不可忽视；其次地震动的远近场输入、斜入射、垂直入射等输入方式对模拟结果也存在一定的影响。如何根据边坡的实际场地条件及震源机制选取合理的地震动输入方式仍需进一步的研究。

(3)离散元边坡加速度响应分析中颗粒加速度的获取方式。颗粒离散元的离散特性使其在加速度响应分析时无法直接导出颗粒的加速度值。目前两种常用的方法均存在一定的误差，不平衡力转化法得到的结果随机性及离散性较大，速度时程微分法得到的颗粒速度时程曲线受多个因素影响。因此更准确的获取模型破坏过程中颗粒的加速度值的方法有待进一步研究。

最后，颗粒离散元模拟边坡动力响应的研究工作仍有很大的发展空间。首先，将离散元颗粒流可模拟岩土体大变形及细观破坏的优势与动力分析相结合。例如离散元模拟过程中可检测坡体内裂纹数量、颗粒间接触的破坏数量、接触力强度的变化情况等。将颗粒流的这些细观分析优势与边坡动力响应结合起来，而不是仅仅局限于传统的参量的分析，有助于从不同角度分析边坡的破坏机理。其次，颗粒离散元与其他数值方法耦合模拟边坡动力响应的相关研究已初露头角，如颗粒离散元与有限元、有限差分法、光滑粒子流方法等的耦合可以使边坡动力响应的数值方法面向更多的研究，也为边坡的稳定性分析提供了新的途径。