工程地质学报  2018, Vol. 26 Issue (6): 1463-1472   (3773 KB)    
Article Options
  • PDF (3773 KB)
  • Full Text HTML
  • Abstract
  • Figures
  • References
  • History
  • 收稿日期:2017-09-18
  • 收到修改稿日期:2018-03-19
  • 扩展功能
    把本文推荐给朋友
    加入引用管理器
    Email Alert
    文章反馈
    浏览反馈信息
    本文作者相关文章
    杨婷婷
    杨永森
    邱流潮

    引用本文  

    杨婷婷, 杨永森, 邱流潮. 2018. 基于物质点法的土体流动大变形过程数值模拟[J]. 工程地质学报, 26(6): 1463-1472. doi: 10.13544/j.cnki.jeg.2017-453.
    YANG Tingting, YANG Yongsen, QIU Liuchao. 2018. MPM based numerical simulation of large deformation process of soil flow[J]. Journal of Engineering Geology, 26(6): 1463-1472. doi: 10.13544/j.cnki.jeg.2017-453.

    基于物质点法的土体流动大变形过程数值模拟
    杨婷婷, 杨永森, 邱流潮    
    中国农业大学水利与土木工程学院 北京 100083
    摘要:土体滑坡作为一种自然地质灾害,受自然因素和人类活动的影响在我国时有发生,给周围居民的生命和财产安全带来了很大威胁,日益受到人们的广泛关注。滑坡防治也逐渐成为工程研究的热点之一。土体本质上是一种具有复杂组成结构的颗粒材料堆积体,通过对颗粒流动的模拟可以深入理解自然界中的土体流动现象,如滑坡、泥石流等,进而预测灾害破坏范围及改进相应工程防护措施。但由于土体流动是一个涉及大变形及大位移的复杂流动过程,传统的基于网格的有限元法(FEM)由于网格畸变,并不适合这类问题的研究。本文采用物质点法(MPM)模拟土体流动大变形问题。作为一种源于particle-in-cell(PIC)法的无网格法,兼具欧拉法和拉格朗日法的优点,因而,物质点法在处理大变形问题上具有独特的优势。目前,国内外利用物质点法模拟边坡滑动问题已有不少研究,但对相关参数进行敏感性分析的较少。本文基于物质点法模拟了黏性土体及无黏性土体流动大变形问题,并进行了参数敏感性分析,包括土体材料的内摩擦角、黏聚力、高宽比、底面坡度对土体滑动距离的影响规律。本文计算中采用Drucker-Prager(DP)弹塑性本构模型描述土的非线性特性。研究结果表明:(1)基于物质点法得到的土体的流动形态、滑动距离以及自然休止角等数值模拟结果均与文献中的实验结果基本吻合,验证了物质点法模拟土体大变形力学行为的精度及有效性;(2)随着内摩擦角、黏聚力的增大,滑动距离相应减小;(3)坡度对边坡稳定的影响是显著的,随着底面坡度的增大,滑动距离相应增大;(4)当土柱高宽比较小时,与滑动距离呈线性增长关系。其中,内摩擦角和黏聚力反映了土体的抗剪切性能,因此通过工程措施提高土体的抗剪能力可以降低土体滑坡带来的危害。研究结果为探索土体滑动破坏规律,降低滑动破坏范围提供了可靠的参考。
    关键词土体流动    物质点法    大变形    滑动距离    无网格法    
    MPM BASED NUMERICAL SIMULATION OF LARGE DEFORMATION PROCESS OF SOIL FLOW
    YANG Tingting, YANG Yongsen, QIU Liuchao    
    College of Water Resources & Civil Engineering, China Agricultural University, Beijing 100083
    Abstract: As a natural geologic disaster, soil landslide is affected by natural factors and human activities. It brings a great threat to the safety of the surrounding residents' life and property, attracting more worries from people. Landslide prevention and control are also one of the hot spots in engineering research. Soil is essentially a kind of deposits consists of granular material with complicated composition and structure. By means of numerical simulation of granular material flows, we can well understand the phenomena of soil flows in the nature, such as landslide and debris flows. It is therefore possible to predict the landslide hazard zone and to improve the design of engineering protection measurements. However, the finite element method(FEM) is not suitable for simulating soil flow with large deformation and large displacement due to the finite element method is sensitive to mesh distortion. In this paper, we simulate the soil flows using the material point method(MPM). MPM is a mesh-free method and originating from the particle-in-cell(PIC)method. It combines the Eulerian description and Lagrangian description and has distinct advantages in solving the large deformation and large displacement problem. Currently, there have been many studies on the simulation of slope sliding using MPM, but less attentions has paid to sensitivity analysis on relevant parameters. In this paper, the large deformations of the cohesive soil and the non-cohesive soil slopes due to gravity are numerically investigated with MPM.The influence of the internal friction angle, the cohesive force, the aspect ratio and bottom surface gradient on the landslide run-out are analysed. In this MPM simulation, the elasto-palstic constitutive models based on the Drucker-Prager(DP)yield criterion is used for modeling nonlinear characters of soil flows. The Drucker-Prager yield criterion is a pressure-dependent model for determining whether a material has failed or undergone plastic yielding. The yielding surface of the Drucker-Prager criterion may be considered depending on the internal friction angle of the material and its cohesion. The simulated results show the following featurs. (a) There are well agreements in the flow pattern, the sliding distance and natural angle of repose between the simulated results and the corresponding experimental results, which validates the ability of MPM to modeling soil flows. (b) The sliding distance decreases with the increase of the internal friction angle and cohesive force. (c) The steepness of the soil slope has significant influence on its stability. The sliding distance increases with the increase of the steepness of the soil slope. (d) For a comparatively small aspect ratio, the sliding distance increases linearly with aspect ratio. It is noteworthy that the internal friction angle and cohesive force of soil reflect its shear resistance. Therefore, the damage due to landslide can be reduced by improving the shear resistance of soil through engineering measurements. The results of the simulation provide a reliable reference to explore the law of the soil sliding hazardous behavior and reduce the sliding damage range.
    Key words: Soil flow    Material point method    Large deformation    Run-out    Mesh-free method    

    0 引言

    颗粒材料在人为扰动或自然因素(如降雨、地震等)干扰下,受重力作用会产生失稳流动的现象,这种现象在自然界中意味着自然灾害的发生,如泥石流、雪崩、滑坡等,会危害人民生命和公共财产安全。一直以来,土体颗粒流动导致的滑坡问题由于其危害的严重性而受到广泛的关注(曹文等,2017王畯才等,2017)。为了研究颗粒材料失稳流动现象,目前常用的数值模拟手段包括有限元法(FEM)、有限体积法(FVM)等基于网格的方法以及离散元法(DEM)、物质点法(MPM)等无网格法。由于无网格法采用无拓扑连接关系的点来离散材料区域,并基于这些点建立场变量插值,构造的形函数随插值节点的移动而变化,不需要划分网格,所以在处理大变形问题的时候,可以彻底或部分地消除对网格的依赖,避免网格畸变,一定程度上提高计算的精度和效率,更适合用来研究颗粒材料流动问题。离散元法可以更真实地反映岩土结构的细观结构特点,而且采用显式差分法计算动力平衡方程,可以很方便地对非线性或大变形问题进行计算,因此在岩土工程中应用十分广泛。但是离散元法也有一定的局限性,主要是颗粒细观参数与相应宏观力学参数间的标定理论还不够完善,同时计算量受颗粒数目的限制不能过大,在实际工程的大尺寸模拟中还无法应用。物质点法(MPM)由Sulsky et al. (1994)首次提出,采用拉格朗日法对材料进行离散,由离散后的质点来代表材料并携带所有的物质信息如质量、速度、应力、应变、位移等;采用欧拉法对空间进行背景网格划分,在网格节点上求解动量方程和计算空间导数,兼具拉格朗日法和欧拉法的优点,可以避免对流项的干扰和背景网格的畸变。对于大变形及破碎断裂问题的处理有着显著的优势,可用于研究高速碰撞、爆炸、新型材料的模拟、地质灾害等问题。本文拟采用物质点法对颗粒流动问题进行数值模拟。

    自从物质点法提出后,被广泛用于各个领域。在岩土工程方面,Bardenhagen et al.(2000, 2001)利用物质点法建立了颗粒性材料模型,考虑了颗粒内部变形及颗粒间相互作用,之后又对接触算法进行了优化,模拟了压力在理想的颗粒性材料中的传播过程。Cummins et al. (2002)提出了物质点隐式解法,更好地描述了颗粒间的碰撞问题。这些从材料细观结构入手的研究为岩土问题的持续探索与讨论奠定了基础。此后,Daphalapurkar et al. (2007)基于广义插值的物质点法研究了动态裂纹的扩展问题,Johansson et al. (2007)分析了断层面引起的土体变形问题,Dong et al. (2015)介绍了物质点法在离岸岩土工程上的应用。相比之下,国内虽然起步较晚,但发展迅速。张洪武等(2009)模拟了饱和多孔介质与固体间的相互作用,分析了饱和土体的动力学响应问题,黄鹏(2010)基于OpenMP开展物质点法并行计算,研究了金属及岩土的动力冲击和侵彻问题,张雄等(2016)以刚性容器内液体晃动为例研究了流固耦合问题,王双等(2016)史卜涛等(2016)将物质点强度折减法应用于各类边坡工程问题中,从而实现了边坡稳定性评价。

    关于物质点模拟边坡滑动问题,目前国内外已有不少研究,利用物质点法研究了边坡失效、泥石流等问题(Woo,2009Andersen et al., 2010)。Bandara et al.(2016)考虑了非饱和边坡在雨水入渗情况下的滑动问题,一些学者利用物质点法分析了高边坡、人工堆填土边坡的滑动破坏问题(孙玉进等,2015张巍等,2017)等,但对相关参数的影响分析不多。本文基于物质点法对无黏性颗粒柱塌陷及黏性土体变形问题进行模拟,并进行参数敏感性分析,从而找出不同参数对破坏程度的影响规律,有助于更好地理解土体滑动破坏行为,为相关工程防护提供参考。

    1 物质点法
    1.1 运动方程及其积分弱形式

    当采用更新拉格朗日格式描述物体运动过程时,若不考虑热量交换,连续体运动的动量方程如下:

    $ \frac{{\partial {\sigma _{ij}}}}{{\partial {x_j}}} + \rho {b_i} = \rho {{\ddot u}_i}\left( {i,j = 1,2,3} \right) $ (1)

    其中,σ为Cauchy应力;ρ为当前密度;x为空间坐标;b为作用于物体单位质量上的体积力,$\ddot u$为物体的加速度。

    边界条件为:

    $ \begin{array}{*{20}{c}} {\left( {{n_j}{\sigma _{ij}}} \right)\left| {_{{\mathit{\Gamma }_t}}} \right. = {{\bar t}_i}}\\ {{v_i}\left| {_{{\mathit{\Gamma }_u}}} \right. = {{\bar v}_i}} \end{array} $ (2)

    其中,ΓtΓu分别为给定面力边界和给定位移边界。

    从数学的角度来看,动量方程作为偏微分方程,要满足在求解域Ω内处处可解的条件,直接求解是十分困难的。这里应用加权余量法,建立了和原微分方程及定解条件等价的弱形式,只要求动量方程在某种平均意义下满足,在此基础上再进行求解。当引入比应力σijs=σij/ρ和比边界面力$\bar t_i^s = {\bar t_i}/\rho $后,动量方程(1)和给定面力边界条件(2)的等效积分形式为:

    $ \int\limits_\Omega {\rho {{\ddot u}_i}\delta {u_i}{\rm{d}}V} + \int\limits_\Omega {\rho \sigma _{ij}^s\delta {u_{i,j}}{\rm{d}}V} - \int\limits_\Omega {\rho {b_i}\delta {u_i}{\rm{d}}V} - \int\limits_{{\Gamma _t}} {\rho \bar t_i^s\delta {u_{i,j}}{\rm{d}}A} = 0 $ (3)

    1.2 物质点离散

    物质点法将计算域Ω离散为一系列物质点,并将连续体的质量集中于这些物质点上,物质点携带了速度、应力、应变等所有物质信息,可以代替原来的整体,此时连续体的密度可表示为:

    $ \rho \left( {{x_i}} \right) = \sum\limits_{p = 1}^{{n_p}} {{m_p}\delta \left( {{x_i} - {x_{ip}}} \right)} $ (4)

    其中,np为质点总数;mp为质点P的质量;δ为Dirac Delta函数;xip为质点p的坐标。

    一般情况下,在计算过程中物质点质量不会改变,质点也不会凭空增加或减少,自动满足质量守恒定律。若将式(4)代入式(3),可将积分形式的虚功方程转化为物质点信息求和的形式(张雄等,2013):

    $ \begin{array}{l} \sum\limits_{p = 1}^{{n_p}} {{m_p}{{\ddot u}_{ip}}\delta {u_{ip}}} + \sum\limits_{p = 1}^{{n_p}} {{m_p}\sigma _{ijp}^s\delta {u_{ip,j}}} - \sum\limits_{p = 1}^{{n_p}} {{m_p}{b_{ip}}\delta {u_{ip}}} \\ - \sum\limits_{p = 1}^{{n_p}} {{m_p}\bar t_{ip}^s{h^{ - 1}}\delta {u_{ip}}} = 0 \end{array} $ (5)

    式中,h是为了将式(3)左端最后一项边界积分转化为体积分而引入的假想边界层厚度。

    物质点法求解动量方程时,在时间步起始与终止时,利用欧拉法划分网格来集成动量方程和更新位置及其他物理量,这样可以避免网格畸变。在每个时间步内,采用拉格朗日格式,物质点与背景网格完全固连,两者一起运动,因此可以通过建立在背景网格节点上的有限元插值形函数NI(xi)来实现信息在质点与背景网格之间的相互映射,从而避免欧拉对流项的处理。

    这里背景网格采用八节点六面体单元,则背景网格节点I的形函数为:

    $ {N_l} = \frac{1}{8}\left( {1 + \xi {\xi _I}} \right)\left( {1 + \eta {\eta _I}} \right)\left( {1 + \xi {\xi _I}} \right),I = 1,2, \cdots ,8 $ (6)

    其中,ξIηIζI为节点I的自然坐标,其取值分别为±1。

    对于物理量β(这里代表位移及其导数、虚位移),物质点与背景网格节点的插值关系为:

    $ {\beta _l} = \sum\limits_{P = 1}^{{n_p}} {{\beta _p}{N_I}\left( {{x_p}} \right)} $ (7)

    其中,NI(xp)为节点I的形函数在质点p处的值,用带有下标I的量表示网格节点的变量,带有下标p的量表示质点携带的变量。

    将式(7)代入式(5)中,得到背景网格节点的运动方程:

    $ {{\dot p}_{iI}} = f_{iI}^{{\mathop{\rm int}} } + f_{iI}^{{\rm{ext}}},{x_I} \notin {\mathit{\Gamma }_u} $ (8a)

    $ f_{iI}^{{\mathop{\rm int}} } = - \sum\limits_{p = 1}^{{n_p}} {{N_{IP,j}}{\sigma _{ijp}}\frac{{{m_p}}}{{{\rho _p}}}} $ (8b)

    $ f_{iI}^{{\rm{ext}}} = \sum\limits_{P = 1}^{{n_p}} {{m_p}{N_{Ip}}{b_{ip}}} + \sum\limits_{p = 1}^{{n_p}} {{N_{Ip}}{{\bar t}_{ip}}{h^{ - 1}}\frac{{{m_p}}}{{{\rho _p}}}} $ (8c)

    其中,PiI为第I个网格节点在i方向的动量;fiIint为节点内力;fiIext为节点外力;σijp为质点p的应力。

    对于式(8)可以采用逐步积分法求解,如中心差分法、Newmark法和Wilson法等(Nairn,2003)。中心差分法作为显式格式可以由当前时刻状态量直接求解计算下一刻系统状态量,其优点是不需要求解联立方程组,使计算过程大大简化。但是显式方法是条件稳定算法,只有时间划分在临界步长范围内才能保证稳定,一般多用于冲击爆炸等需要关注应力波传播过程的高频响应问题中。对于隐式求解方法,需要求解同时涉及系统的当前时刻状态量和下一时刻状态量的方程来确定下一时刻状态量,需要对方程组进行迭代,尽管单步内计算量增加,但是隐式方法一般是无条件稳定算法,因此可以加大时间步长而不影响最终结果的稳定性,适用于求解长时间低频响应问题(张雄等,2013)。本文采用显式方法求解。

    1.3 计算流程

    通过背景网格形函数可以将积分结果映射回质点以更新质点的位置和速度,即分别利用背景网格的加速度场和速度场更新质点速度和位置。根据应力更新可以在每个时间步开始或结束时进行,可以把物质点法分为3种求解格式:在时间步开始时进行应力更新的格式USF(Update Stress First)、在时间步结束时进行应力更新的格式USL(Update Stress Last)、针对USL改进后的格式MUSL(Modified Stress Last)。其中,USL数值计算不稳定,且具有较强的数值耗散,MUSL在USL的基础上做了改进,但总体能量依旧趋于缓慢耗散,而USF具有较好的能量守恒特性(Bardenhagen,2002Nairn,2003)。因此,本文计算采用USF格式,具体步骤如下:

    (1) 利用式(7)将物质点的质量和动量映射到背景网格上,得到网格节点上的质量和动量,在施加边界条件后,得到网格节点的速度viIn-1/2

    $ v_{iI}^{n - 1/2} = \frac{{p_{iI}^{n - 1/2}}}{{m_I^n}} = \frac{1}{{m_I^n}}\sum\limits_{P = 1}^{{n_p}} {{m_p}v_{ip}^{n - 1/2}N_{Ip}^n} $ (9)

    (2) 根据式(9)来计算各物质点的应变增量Δεijpn-1/2和旋量增量ΔΩijpn-1/2

    $ \Delta \varepsilon _{ijp}^{n - 1/2} = \Delta t/2\sum\limits_{I = 1}^8 {\left( {N_{Ip,j}^nv_{iI}^{n - 1/2} + N_{Ip,i}^nv_{jI}^{n - 1/2}} \right)} $ (10a)

    $ \Delta \Omega _{ijp}^{n - 1/2} = \Delta t/2\sum\limits_{I = 1}^8 {\left( {N_{Ip,j}^nv_{iI}^{n - 1/2} - N_{Ip,i}^nv_{jI}^{n - 1/2}} \right)} $ (10b)

    (3) 应用本构积分法更新各物质点的应力和密度,其中在本构方程中采用焦曼应力率来消除刚体转动对应力的影响。

    (4) 根据式(8)计算节点内力fiIint,n和外力fiIext,n,并在背景网格上积分动量方程:

    $ P_{iI}^{n + 1/2} = P_{iI}^{n - 1/2} + \left( {f_{iI}^{{\mathop{\rm int}} ,n} + f_{iI}^{{\rm{ext}},n}} \right) \cdot \Delta t $ (11)

    (5) 将式(11)计算结果利用背景网格节点与物质点间的插值关系映射回物质点,更新质点位置、速度、加速度。

    (6) 重置背景网格,从第一步开始新时间步的计算。

    1.4 材料的本构关系

    岩土可看作颗粒体堆积或胶结而成的多相体材料,其抗压强度远大于抗拉强度,受剪时颗粒会膨胀,实验表明,其屈服极限值不是一个常数,而是与该平面上的正应力σn有关,此时金属材料的Tresca和Mises屈服准则不再适用于这类材料,因为这两种模型都与静水压力无关(余同希,1989庄茁,2002)。

    Mohr-Coulomb(MC)屈服准则可以看作是Coulomb摩擦定律在连续体和多轴应力-应变状态中的拓展,目前在岩土工程中得到了广泛的应用。其屈服准则可以表示为:

    $ {\tau _n} = c - {\sigma _n}\tan \varphi $ (12)

    其中,c为黏聚力;φ为内摩擦角;σ为材料某点应力;n为某一斜截面外法向单位向量;τnσn分别为斜面正应力、剪应力(图 1)。

    图 1 MC模型临界屈服状态 Fig. 1 The MC model critical yield state

    当满足条件σ3σ2σ1时,可确定应力Mohr圆与式(12)所表示的直线相切时,对应的屈服函数为:

    $ f\left( \sigma \right) = {\sigma _1} - {\sigma _3} + \left( {{\sigma _1} + {\sigma _3}} \right)\sin \varphi - 2c\cos \varphi = 0 $ (13)

    Drucker-Prager(DP)模型是常用的弹塑性本构模型,该模型的屈服面与压力有关,还考虑了剪胀性和非关联塑性的影响。它所对应的屈服函数为:

    $ f\left( \sigma \right) = \sqrt {{J_2}} + \alpha {I_1} - k = 0 $ (14)

    其中,I1=σii为应力第一不变量,J2=sijsij/2为应力第二不变量;αk为材料参数。两模型间参数αkcφ的关系为:

    $ \alpha = \frac{{2\sin \varphi }}{{\sqrt 3 \left( {3 \pm \sin \varphi } \right)}} $ (15)

    $ k = \frac{{6c\cos \varphi }}{{\sqrt 3 \left( {3 \pm \sin \varphi } \right)}} $ (16)

    对于干颗粒材料,可能处于堆积静止、倾斜流动、剧烈晃动等不同状态,其颗粒行为表现可以在固体、液体、气体间转化,所以如何确定其本构模型,一直存有争议。其中,屈服准则是描述颗粒流的一个主要依据。

    当剪应力小于临界剪应力,颗粒处于静止状态,可被视为固体;当剪应力大于临界剪应力时,流动发生且与剪切率密切相关。Jop et al. (2006)指出,描述颗粒流的3D屈服准则的形式与DP模型类似。相对于MC准则,DP模型可能更符合颗粒流的流动特性。

    此外,DP模型所对应屈服面更为光滑(图 2),可看作是MC准则在π平面上去除奇异点的光滑近似处理结果,解决了因MC屈服面不光滑给数值计算带来的困难。因此我们采用DP模型作为材料的本构关系。

    图 2 π平面上的屈服面 Fig. 2 Shape of yield criteria on the π-plane

    2 数值计算与分析

    在这一部分,我们采用物质点法模拟土体大变形及破坏问题,数值计算基于张雄等(2013)开发的MPM3D程序。首先,利用物质点法模拟实验室中颗粒柱塌陷的过程,该过程类似于无黏性土坡在重力作用下的滑动问题,通过比较实验与模拟结果,可以验证物质点法适用于模拟无黏土体失稳问题。其次,将物质点法应用于黏性土体塌陷问题中,模拟结果与文献中的SPH法模拟结果作对比来讨论物质点法模拟黏性土体失稳问题的准确性。

    2.1 无黏性土体大变形过程模拟

    2010年,Mangeney等进行了一个简易的颗粒柱崩塌失稳实验,实验操作流程如图 3所示,它包含一个宽20 cm、高14 cm的深槽,一个可控制的闸门,和一个可以使砂砾自由流动的平面。实验材料为直径600~800 μm的干的玻璃微珠,来模拟自然界中无黏性的岩土崩塌过程。整个实验过程使用高速摄影机记录下来,用于分析具有自由面的颗粒柱在不同时刻的流动形态以及最终的滑动距离,来估计实际工程中滑坡土方量、判断滑坡堆积形态及预判滑坡是否可能威胁临近建筑物。其中,具体实验的细节和试验方法可以参考文献。

    图 3 颗粒流实验的操作流程图(Mangeney et al., 2010Farin et al., 2014) Fig. 3 Operation process of granular experiments

    为了能够与该实验结果进行对比,数值模拟中的模型参数及尺寸与实验保持一致。模型为高14 cm,宽20 cm的棱柱体,遵从DP屈服准则,具体材料参数见表 1。模型底部为固定边界来模拟坚硬基岩,约束左侧及前后两侧法线方向的位移来模拟平面应变情况,上部为自由边界来模拟颗粒自由流动。物质点离散时,采用背景网格为1 cm,质点间距0.5 cm,即在每一个网格单元里有8个物质质点。计算时采用USF显式计算格式。

    表 1 颗粒材料参数 Table 1 Granular material parameters

    当底面坡度θ为0时,通过比较不同时刻实验与物质点法模拟的颗粒流在垂直方向的自由面轮廓图(图 4),我们可以很好的观察在水平面上颗粒流动的形态,颗粒不断向右侧移动,在经历了大变形后最终静止形成了自然休止角。速度云图表明流动主要集中在自由表面附近的浅层子区域内,且滑动带前缘位置的颗粒水平位移最大。当闸门开启后,颗粒的速度快速增加,t=0.18 s时,颗粒最大速度接近1 m · s-1,随后以一段稳定速度流动,接着速度逐渐减小,t=0.5 s时,颗粒最大速度 < 0.2 m · s-1,但还是存在颗粒缓慢流动的现象,不断对最终的流动形态产生轻微的改变,直至t=1.02 s时颗粒达到最终平衡状态。

    图 4 颗粒塌陷时自由表面对比图(a),(b),(c)及速度分布图(d),(e),(f) Fig. 4 Free surface contrast(left) and plastic strain distribution(right) of the collapsing sand column a. t=0.18 s; d. t=0.18 s; b. t=0.30 s; e. t=0.30 s; c. t=0.50 s; f. t=0.50 s

    模拟结果与试验吻合较好,计算的滑动距离与自然休止角均与实验保持基本一致,说明利用物质点法模拟无黏性土体材料的大变形是合适的。在流动前期,模拟的颗粒柱流动发展较快一些,而在流动快结束时,两条轮廓线基本重合。其中,在前期两者间的差异可能是数值计算中没有考虑在实验中闸门的开启过程而导致的。左侧土体的厚度有轻微的下降,而实验中厚度保持不变;相比于实验,前端流动的形态更圆润一点

    2.2 黏性土体大变形过程模拟

    2008年,Bui et al.(2008)为了证明SPH法能够很好地模拟黏性与非黏性土体大变形过程,处理好土颗粒间的拉伸问题,建立了一个高2 m、宽4 m的矩形土体模型,模拟其在重力作用下的变形过程。在模拟黏性土体时,他们将土体黏聚力参数由0增加至5 kPa。

    用物质点法模拟黏性土体大变形及破坏过程时,应注意建立的模型及参数同上述保持一致,模拟其在重力作用下的破坏过程。其中,土材料参数见表 2,将土体黏聚力参数增加至5 kPa以模拟黏性土壤,采用D-P模型屈服准则,按照内接圆等效计算D-P模型参数。该模型采用5000个质点进行离散,背景网格单元尺寸为0.08 m,相邻物质点间距为0.04 m,其他计算过程及属性设置同2.1。

    表 2 土壤材料参数 Table 2 Soil material parameters

    为了验证模拟结果是否准确,我们将基于物质点法的计算结果与SPH法模拟结果(Bui et al., 2008)进行了比较(图 5)。变形过程根据体系动能的变化可分为加速滑动、减速滑动以及稳定阶段3个阶段(孙玉进等,2015),其中彩色云图表示塑性应变的分布。从图中可以清楚的观察到,随着破坏过程的进行,塑性区在不断扩大,且逐渐形成明显的剪切带,直至土颗粒在开始塌陷2 s后停止移动,变形处于平衡状态。经过对比,物质点法和SPH法均能有效模拟该滑坡失稳过程,而且这两种算法在不同时刻的边坡外形和剪切带能基本吻合。在这一部分模拟中,土体表现更像是黏性土壤,而非颗粒性材料。由于黏性土体颗粒之间存在着拉应力,故而变形程度较小,随着黏性的增大,土体甚至可以保持站立状态,不发生变形。

    图 5 SPH法(左)与物质点法(右)在不同时刻的变形对比图 Fig. 5 Deformation contrast at different times with SPH method(left) and MPM(right)

    前述两个数值模拟结果表明,采用D-P模型及其材料参数是基本准确的,物质点法能够模拟土体大变形力学行为,无论是黏土或无黏土。此外,物质点法计算效率高,可以避免有限元中的网格畸变以及SPH法计算中的拉伸失稳定问题。

    2.3 参数敏感性分析

    由于土体失稳破坏运动特性直接受其材料的特征参数影响,为了寻找力学参数与土体滑动距离之间的关系,本文采用物质点法研究内摩擦角、黏聚力、坡度、高宽比对土体失稳破坏后滑动距离的影响,通过改变相应材料参数进行多组数值模拟来实现这个过程。在研究高宽比对滑动距离的影响时,由于物体的尺寸不同,不能保证初始时刻的滑动距离(物体宽度)相同,因此在这里研究的滑动距离统一都是相对滑动距离Lr=(L2-L1)/L1,其中,L1表示颗粒柱初始宽度,L2表示对应时刻(t=0.5 s或1 s)的滑动距离;模型建立过程参考2.1,除需要进行敏感分析的参数外,其他参数与2.1保持一致,图 6~ 图 8是得到的模拟结果。

    图 6 内摩擦角、黏聚力与滑动距离间的关系 Fig. 6 The relationship between angle of friction, cohesive strength and sliding distance

    图 7 底面坡度与滑动距离间的关系 Fig. 7 The relationship between the bottom slope and sliding distance

    图 8 高宽比与滑动距离间的关系(t=0.5 s) Fig. 8 The relationship between aspect ratio and sliding distance

    图 6a的总体趋势来看,在同一时刻,随着内摩擦角的增大,滑动距离相应减小。当内摩擦角增大或减小到一定程度时,滑动距离倾向于稳定,即内摩擦角对滑动距离的影响减弱;观察图 6b可知,当黏聚力增大时,对应的滑动距离减小,且变化速率由快变慢,当黏聚力增加到一定程度,滑动距离趋于0,即基本不发生变形;由图 7图 8可知,坡度与高宽比对边坡稳定的影响是显著的,随着坡度变得越来越陡,滑动距离也会相应增大,增速由快变慢;从图 8中我们可以看到,在一定范围内,高宽比与相对滑动距离间呈现出线性增长关系。研究表明,当高宽比较小时(一般 < 1.7)时,与滑动距离间呈线性关系,颗粒柱上表面存在一个内部静止区与外部流动区,这和田成武(2014)中的理论分析结果是一致的。

    3 结论

    本文基于物质点法对有黏及无黏土体流动大变形问题进行数值模拟,首先验证了物质点法适用于模拟土体大变形力学行为。其次对相关参数进行了影响分析,包括内摩擦角、黏聚力、高宽比、坡度。研究表明:随着内摩擦角、黏聚力的增大,滑动距离相应减小,变化速率由快变慢;随着高宽比、平面坡度的增大,土体滑动距离也随之增加,滑动破坏范围变大。高宽比与相对滑动距离在一定范围内呈现线性增长关系。其中内摩擦角和黏聚力反映了土体的抗剪切性能,因此可以采取一些工程措施来提高土体的抗剪能力,如排水固结,加筋锚固等,从而减少或降低土体滑坡带来的危害。

    土体作为一种不连续、非均质、各向异性的复杂多孔介质,涉及到固、气、液三相的组成与交互,在这里我们将其简化为单一的固体颗粒性材料,忽略了其他因素的影响,因此模拟结果不能完全符合实际情况,水分、土壤孔隙结构等对土体滑动的影响还有待今后进一步研究。

    参考文献
    Andersen S, Andersen L. 2010. Modelling of landslides with the material point method[J]. Computational Geosciences, 14(1): 137~147.
    Bandara S, Ferrari A, Lalou L. 2016. Modelling landslides in unsaturated slopes subjected to rainfall infiltration using material point method[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 40(9): 1358~1380. DOI:10.1002/nag.v40.9
    Bardenhagen S G, Brackbill J U, Sulsky D. 2000. The material-point method for granular materials[J]. Computer Methods in Applied Mechanics and Engineering, 187(3-4): 529~541. DOI:10.1016/S0045-7825(99)00338-2
    Bardenhagen S G, Guilkey J E, Roessig K M, et al. 2001. An improved contact algorithm for the material point method and application to stress propagation in granular material[J]. CMES-Computer Modeling in Engineering & Sciences, 2(4): 509~522.
    Bardenhagen S G. 2002. Energy conservation error in the material point method for solid mechanics[J]. Journal of Computational Physics, 180(1): 383~403. DOI:10.1006/jcph.2002.7103
    Bui H H, Fukagawa R, Sako K, et al. 2008. 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, 32(12): 1537~1570. DOI:10.1002/nag.v32:12
    Cao W, Li W C, Tang B, et al. 2017. PFC study on building of 2D and 3D landslides models[J]. Journal of Engineering Geology, 25(2): 455~462.
    Cummins S J, Brackbill J U. 2002. An implicit particle-in-cell method for granular materials[J]. Journal of Computational Physics, 180(2): 506~548. DOI:10.1006/jcph.2002.7101
    Daphalapurkar N P, Lu H B, Coker D, et al. 2007. Simulation of dynamic crack growth using the generalized interpolation material point(GIMP)method[J]. International Journal of Fracture, 143(1): 79~102. DOI:10.1007/s10704-007-9051-z
    Dong Y, Ma J, Wang D, et al. 2015. Assessment of applicability of the material point method in offshore geotechnical engineering[C]//Computer Methods and Recent Advances in Geomechanics. kyoto: [s.n.]: 117-122.
    Farin M, Mangeney A, Roche O. 2014. Fundamental changes of granular flow dynamics, deposition, and erosion processes at high slope angles:Insights from laboratory experiments[J]. Journal of Geophysical Research-Earth Surface, 119(3): 504~532. DOI:10.1002/2013JF002750
    Huang P. 2010. Material point method for metal and soil impact dynamics problems[D]. Beijing: Tsinghua University.
    Ionescu L R, Mangeney A, Bouchut F, et al. 2015. Viscoplastic modelling of granular column collapse with pressure dependent rheology[J]. Journal of Non-Newtonian Fluid Mechanics, 219: 1~18. DOI:10.1016/j.jnnfm.2015.02.006
    Johansson J, Konagai K. 2007. Fault induced permanent ground deformations:Experimental verification of wet and dry soil, numerical findings' relation to field observations of tunnel damage and implications for design[J]. Soil Dynamics and Earthquake Engineering, 27(10): 938~956. DOI:10.1016/j.soildyn.2007.01.007
    Jop P, Forterre Y, Pouliquen O. 2006. A constitutive law for dense granular flows[J]. Natrue, 441(7094): 727~730. DOI:10.1038/nature04801
    Mangeney A, Roche O, Hungr O, et al. 2010. Erosion and mobility in granular collapse over sloping beds[J]. Journal of Geophysical Research-Earth Surface, 115: F03040.
    Nairn J A. 2003. Material point method calculations with explicit cracks[J]. Computer Modeling in Engineering & Sciences, 4(6): 649~663.
    Shi B T, Zhang Y, Zhang W. 2016. Strength reduction material point method for slope stability[J]. Chinese Journal of Geotechnical Engineering, 38(9): 1678~1684.
    Sulsky D, Chen Z, Schreyer H L. 1994. A particle method for history-dependent materials[J]. Computer Methods in Applied Mechanics and Engineering, 118(1): 179~196.
    Sun Y J, Song E X. 2015. Simulation of large-displacement landslide by material point method[J]. Chinese Journal of Geotechnical Engineering, 37(7): 1218~1225.
    Tian W C. 2014. The application of the material point method in geotechnical engineering problems[D]. Hangzhou: Zhejiang University.
    Wang J C, Lu K L, Zhu D Y. 2017. Indoor modeling test for regularity of deposit position of landslide-debris avalanches[J]. Journal of Engineering Geology, 25(6): 1509~1517.
    Wang S, Li X C, Shi L, et al. 2016. Material point strength reduction method and its application to slope engineering[J]. Rock and Soil Mechanics, 37(9): 2672~2678.
    Woo K S. 2009. Numerical simulation of landslides and debris flows using an enhanced material point method[D]. Washington D. C: University of Washington.
    Yu T X. 1989. Plastic mechanics[M]. Beijng: Higher Education Press.
    Zhang H W, Wang K P, Chen Z. 2009. Material point method for dynamic analysis of saturated porous media(Ⅰ):Coupling material point method[J]. Chinese Journal of Geotechnical Engineering, 31(10): 1505~1511.
    Zhang H W, Wang K P, Chen Z. 2009. Material point method for dynamic analysis of saturated porous media(Ⅱ):Dynamic contact analysis between saturated porous media and solid bodies[J]. Chinese Journal of Geotechnical Engineering, 31(11): 1672~1679.
    Zhang H W, Wang K P. 2010. Material point method for dynamic analysis of saturated porous media(Ⅲ):Two-phase material point method[J]. Chinese Journal of Geotechnical Engineering, 32(4): 507~513.
    Zhang W, Shi B T, Shi B, et al. 2017. Material point method for run-out process simulation of soil landslides and application[J]. Journal of Engineering Geology, 25(3): 815~823.
    Zhang X, Lian Y P, Liu Y, et al. 2013. Material point method[M]. Beijng: Tsinghua University Publishing House: 38~52.
    Zhang X, Zhang F. 2016. Fluid structure interaction incompressible material point method and its applications in sloshing problem[J]. Chinese Journal of Computational Mechanics, 33(4): 582~587.
    Zhang X, Wang T S. 2007. Computational dynamics[M]. Beijng: Tsinghua University Publishing House: 147~169.
    Zhuang Z. 2002. Nonlinear finite elements for continua and structures[M]. Beijng: Tsinghua University Publishing House.
    曹文, 李维朝, 唐斌, 等. 2017. PFC滑坡模拟二、三维建模方法研究[J]. 工程地质学报, 25(2): 455~462.
    黄鹏. 2010.金属及岩土冲击动力学问题的物质点法研究[D].北京: 清华大学. http://cdmd.cnki.com.cn/Article/CDMD-10003-1011280524.htm
    史卜涛, 张云, 张巍. 2016. 边坡稳定性分析的物质点强度折减法[J]. 岩土工程学报, 38(9): 1678~1684.
    孙玉进, 宋二祥. 2015. 大位移滑坡形态的物质点法模拟[J]. 岩土工程学报, 37(7): 1218~1225.
    田武成. 2014.物质点法在岩土工程问题中的应用[D].杭州: 浙江大学. http://www.cnki.com.cn/Article/CJFDTOTAL-JXJC201706200.htm
    王畯才, 卢坤林, 朱大勇. 2017. 基于室内模型试验的滑坡碎屑流堆积分布规律研究[J]. 工程地质学报, 25(6): 1509~1517.
    王双, 李小春, 石露, 等. 2016. 物质点强度折减法及其在边坡中的应用[J]. 岩土力学, 37(9): 2672~2678.
    余同希. 1989. 塑性力学[M]. 北京: 高等教育出版社.
    张巍, 史卜涛, 施斌, 等. 2017. 土质滑坡运动全过程物质点法模拟及其应用[J]. 工程地质学报, 25(3): 815~823.
    张雄, 王天舒. 2007. 计算动力学[M]. 北京: 清华大学出版社: 147~169.
    张雄, 廉艳平, 刘岩, 等. 2013. 物质点法[M]. .
    张雄, 张帆. 2016. 流固耦合不可压物质点法及其在晃动问题中的应用[J]. 计算力学学报, 33(4): 582~587.
    张洪武, 王鲲鹏, 陈震. 2009. 基于物质点方法饱和多孔介质动力学模拟(Ⅰ):耦合物质点方法[J]. 岩土工程学报, 31(10): 1505~1511. DOI:10.3321/j.issn:1000-4548.2009.10.005
    张洪武, 王鲲鹏, 陈震. 2009. 基于物质点方法饱和多孔介质动力学模拟(Ⅱ):饱和多孔介质与固体间动力接触分析[J]. 岩土工程学报, 31(11): 1672~1679. DOI:10.3321/j.issn:1000-4548.2009.11.005
    张洪武, 王鲲鹏. 2010. 基于物质点方法饱和多孔介质动力学模拟(Ⅲ):两相物质点方法[J]. 岩土工程学报, 32(4): 507~513.
    庄茁. 2002. 连续体和结构的非线性有限元[M]. 北京: 清华大学出版社.
    Andersen S, Andersen L. 2010. Modelling of landslides with the material point method[J]. Computational Geosciences, 14(1): 137~147.
    Bandara S, Ferrari A, Lalou L. 2016. Modelling landslides in unsaturated slopes subjected to rainfall infiltration using material point method[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 40(9): 1358~1380. DOI:10.1002/nag.v40.9
    Bardenhagen S G, Brackbill J U, Sulsky D. 2000. The material-point method for granular materials[J]. Computer Methods in Applied Mechanics and Engineering, 187(3-4): 529~541. DOI:10.1016/S0045-7825(99)00338-2
    Bardenhagen S G, Guilkey J E, Roessig K M, et al. 2001. An improved contact algorithm for the material point method and application to stress propagation in granular material[J]. CMES-Computer Modeling in Engineering & Sciences, 2(4): 509~522.
    Bardenhagen S G. 2002. Energy conservation error in the material point method for solid mechanics[J]. Journal of Computational Physics, 180(1): 383~403. DOI:10.1006/jcph.2002.7103
    Bui H H, Fukagawa R, Sako K, et al. 2008. 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, 32(12): 1537~1570. DOI:10.1002/nag.v32:12
    Cao W, Li W C, Tang B, et al. 2017. PFC study on building of 2D and 3D landslides models[J]. Journal of Engineering Geology, 25(2): 455~462.
    Cummins S J, Brackbill J U. 2002. An implicit particle-in-cell method for granular materials[J]. Journal of Computational Physics, 180(2): 506~548. DOI:10.1006/jcph.2002.7101
    Daphalapurkar N P, Lu H B, Coker D, et al. 2007. Simulation of dynamic crack growth using the generalized interpolation material point(GIMP)method[J]. International Journal of Fracture, 143(1): 79~102. DOI:10.1007/s10704-007-9051-z
    Dong Y, Ma J, Wang D, et al. 2015. Assessment of applicability of the material point method in offshore geotechnical engineering[C]//Computer Methods and Recent Advances in Geomechanics. kyoto: [s.n.]: 117-122.
    Farin M, Mangeney A, Roche O. 2014. Fundamental changes of granular flow dynamics, deposition, and erosion processes at high slope angles:Insights from laboratory experiments[J]. Journal of Geophysical Research-Earth Surface, 119(3): 504~532. DOI:10.1002/2013JF002750
    Huang P. 2010. Material point method for metal and soil impact dynamics problems[D]. Beijing: Tsinghua University.
    Ionescu L R, Mangeney A, Bouchut F, et al. 2015. Viscoplastic modelling of granular column collapse with pressure dependent rheology[J]. Journal of Non-Newtonian Fluid Mechanics, 219: 1~18. DOI:10.1016/j.jnnfm.2015.02.006
    Johansson J, Konagai K. 2007. Fault induced permanent ground deformations:Experimental verification of wet and dry soil, numerical findings' relation to field observations of tunnel damage and implications for design[J]. Soil Dynamics and Earthquake Engineering, 27(10): 938~956. DOI:10.1016/j.soildyn.2007.01.007
    Jop P, Forterre Y, Pouliquen O. 2006. A constitutive law for dense granular flows[J]. Natrue, 441(7094): 727~730. DOI:10.1038/nature04801
    Mangeney A, Roche O, Hungr O, et al. 2010. Erosion and mobility in granular collapse over sloping beds[J]. Journal of Geophysical Research-Earth Surface, 115: F03040.
    Nairn J A. 2003. Material point method calculations with explicit cracks[J]. Computer Modeling in Engineering & Sciences, 4(6): 649~663.
    Shi B T, Zhang Y, Zhang W. 2016. Strength reduction material point method for slope stability[J]. Chinese Journal of Geotechnical Engineering, 38(9): 1678~1684.
    Sulsky D, Chen Z, Schreyer H L. 1994. A particle method for history-dependent materials[J]. Computer Methods in Applied Mechanics and Engineering, 118(1): 179~196.
    Sun Y J, Song E X. 2015. Simulation of large-displacement landslide by material point method[J]. Chinese Journal of Geotechnical Engineering, 37(7): 1218~1225.
    Tian W C. 2014. The application of the material point method in geotechnical engineering problems[D]. Hangzhou: Zhejiang University.
    Wang J C, Lu K L, Zhu D Y. 2017. Indoor modeling test for regularity of deposit position of landslide-debris avalanches[J]. Journal of Engineering Geology, 25(6): 1509~1517.
    Wang S, Li X C, Shi L, et al. 2016. Material point strength reduction method and its application to slope engineering[J]. Rock and Soil Mechanics, 37(9): 2672~2678.
    Woo K S. 2009. Numerical simulation of landslides and debris flows using an enhanced material point method[D]. Washington D. C: University of Washington.
    Yu T X. 1989. Plastic mechanics[M]. Beijng: Higher Education Press.
    Zhang H W, Wang K P, Chen Z. 2009. Material point method for dynamic analysis of saturated porous media(Ⅰ):Coupling material point method[J]. Chinese Journal of Geotechnical Engineering, 31(10): 1505~1511.
    Zhang H W, Wang K P, Chen Z. 2009. Material point method for dynamic analysis of saturated porous media(Ⅱ):Dynamic contact analysis between saturated porous media and solid bodies[J]. Chinese Journal of Geotechnical Engineering, 31(11): 1672~1679.
    Zhang H W, Wang K P. 2010. Material point method for dynamic analysis of saturated porous media(Ⅲ):Two-phase material point method[J]. Chinese Journal of Geotechnical Engineering, 32(4): 507~513.
    Zhang W, Shi B T, Shi B, et al. 2017. Material point method for run-out process simulation of soil landslides and application[J]. Journal of Engineering Geology, 25(3): 815~823.
    Zhang X, Lian Y P, Liu Y, et al. 2013. Material point method[M]. Beijng: Tsinghua University Publishing House: 38~52.
    Zhang X, Zhang F. 2016. Fluid structure interaction incompressible material point method and its applications in sloshing problem[J]. Chinese Journal of Computational Mechanics, 33(4): 582~587.
    Zhang X, Wang T S. 2007. Computational dynamics[M]. Beijng: Tsinghua University Publishing House: 147~169.
    Zhuang Z. 2002. Nonlinear finite elements for continua and structures[M]. Beijng: Tsinghua University Publishing House.
    曹文, 李维朝, 唐斌, 等. 2017. PFC滑坡模拟二、三维建模方法研究[J]. 工程地质学报, 25(2): 455~462.
    黄鹏. 2010.金属及岩土冲击动力学问题的物质点法研究[D].北京: 清华大学. http://cdmd.cnki.com.cn/Article/CDMD-10003-1011280524.htm
    史卜涛, 张云, 张巍. 2016. 边坡稳定性分析的物质点强度折减法[J]. 岩土工程学报, 38(9): 1678~1684.
    孙玉进, 宋二祥. 2015. 大位移滑坡形态的物质点法模拟[J]. 岩土工程学报, 37(7): 1218~1225.
    田武成. 2014.物质点法在岩土工程问题中的应用[D].杭州: 浙江大学. http://www.cnki.com.cn/Article/CJFDTOTAL-JXJC201706200.htm
    王畯才, 卢坤林, 朱大勇. 2017. 基于室内模型试验的滑坡碎屑流堆积分布规律研究[J]. 工程地质学报, 25(6): 1509~1517.
    王双, 李小春, 石露, 等. 2016. 物质点强度折减法及其在边坡中的应用[J]. 岩土力学, 37(9): 2672~2678.
    余同希. 1989. 塑性力学[M]. 北京: 高等教育出版社.
    张巍, 史卜涛, 施斌, 等. 2017. 土质滑坡运动全过程物质点法模拟及其应用[J]. 工程地质学报, 25(3): 815~823.
    张雄, 王天舒. 2007. 计算动力学[M]. 北京: 清华大学出版社: 147~169.
    张雄, 廉艳平, 刘岩, 等. 2013. 物质点法[M]. .
    张雄, 张帆. 2016. 流固耦合不可压物质点法及其在晃动问题中的应用[J]. 计算力学学报, 33(4): 582~587.
    张洪武, 王鲲鹏, 陈震. 2009. 基于物质点方法饱和多孔介质动力学模拟(Ⅰ):耦合物质点方法[J]. 岩土工程学报, 31(10): 1505~1511. DOI:10.3321/j.issn:1000-4548.2009.10.005
    张洪武, 王鲲鹏, 陈震. 2009. 基于物质点方法饱和多孔介质动力学模拟(Ⅱ):饱和多孔介质与固体间动力接触分析[J]. 岩土工程学报, 31(11): 1672~1679. DOI:10.3321/j.issn:1000-4548.2009.11.005
    张洪武, 王鲲鹏. 2010. 基于物质点方法饱和多孔介质动力学模拟(Ⅲ):两相物质点方法[J]. 岩土工程学报, 32(4): 507~513.
    庄茁. 2002. 连续体和结构的非线性有限元[M]. 北京: 清华大学出版社.