考虑颗粒形状的离散元法在岩土工程中的应用进展

张世民, 陈万里, 刘嘉英, 孙宏磊

(1.浙江工业大学土木工程学院, 杭州 310014; 2.浙大城市学院工程学院, 杭州 310015; 3.城市基础设施智能化浙江省工程研究中心, 杭州 310015)

岩土工程中往往涉及大量的散粒体材料,散粒体物理力学性质复杂,以往常用的宏观连续体力学理论体系分析其应力应变特征,但这种传统方法忽略了散粒体的离散性,对散粒体内部作用机理认识不足,较难描述应变软化、应力路径相关性等特征,致理论与实际相离。为探究岩土类散粒体力学特征的内部作用机理,节约研究成本,提高研究效率,离散元法(discrete element method, DEM)[1]由Cundall等[2]于1979年提出,并逐渐在岩土数值分析中得到了应用。

离散元法的基本思想为将整个介质材料看作一系列离散的单元,每一单元颗粒均可以独立运动并遵从牛顿运动定律,对每一单元求解其运动方程以进一步迭代计算。通常而言,构成散粒体的颗粒,如岩土工程中的粗粒土、砂土等,具有不规则的几何形状,因此在离散元模拟中,也应考虑单元颗粒的几何形状特性。

早在20世纪30年代,就有学者认识到颗粒形状会对土体力学特性产生一定程度的影响[3-4],土体颗粒由于形状所带来的摩擦咬合,表面凹凸所产生的嵌锁作用都是提供土体抗剪强度不可忽略的重要因素。真实工况中,由于介质材料的几何形状的不规则性,现场试验较难精确还原材料形状参数,许多学者提出了颗粒形状的离散元模拟方法:①以简单形状或形状参数为分类依据的典型形状模拟法[5-7];②基于三维扫描重现颗粒轮廓的实物扫描法[8-9]。采用这些方法进行离散元试验,试验结果均表明了颗粒形状对散粒体材料力学特性产生影响[10-12]。然而大多数的研究都受限于模拟形状方法的局限性,还有待进一步完善。

基于目前岩土工程计算模拟领域离散元法应用的现状,现细致介绍离散元法的模拟方式及模拟成果,分析各模拟颗粒形状的方法的优劣并对各方法的实际应用性进行评价,最后对考虑颗粒形状离散元模拟方法进行展望。

1 离散元方法应用现状

1.1 细观机制分析

采用离散元法进行散粒体物理力学性质模拟,很大一部分是用来分析颗粒材料复杂力学响应的细观机制。目前常对配位数、微观摩擦角、组构张量等进行分析,进而建立宏微观联系。对散粒体材料力学性质的机制分析目前还处于定性阶段,部分学者提出的定量关系的确定方法具有一定的参考意义,为本构模型的建立提供了思路。杨升等[13]模拟直剪试验分析了颗粒初始配位数与内摩擦角的关系以及不同颗粒试样的体积变化规律,从细观上描述了不同形状颗粒的空间活跃程度。张玲等[14]对碎石和黏土颗粒进行数值三轴试验分析了筋材开孔率、筋材抗拉刚度等因素对界面摩擦特性的影响。徐小敏等[15]模拟三轴试验探究了初始杨氏模量和初始剪切模量与颗粒刚度比之间的关系,颗粒刚度比和初始泊松比之间的关系。刘勇等[16]根据粗粒土室内试验结果和数值试验对参数进行标定,研究了颗粒间摩擦系数、细观连接强度以及孔隙率对宏观力学特性的影响。王怡舒等[17]进行了不同滑动摩擦系数和滚动摩擦系数组合下的常规三轴模拟,研究了颗粒接触摩擦对颗粒宏细观力学特性的影响规律。王舒永等[18]利用三维扫描技术建立了不同含石量的土石混合体离散元模型,从微观角度解释了土石混合土剪切过程中一系列的变形并建立了宏细观参数之间的线性方程。路德春等[19]确定了组构张量与应变张量的关系并用隐式方程的方式表达了出来,其中定量关系的确定方法为建立大尺度的本构提供了一种潜在思路。张家伟等[20]通过离散元法生成不同数量和角度的裂缝岩,定量分析了裂隙数量和岩石损伤变量之间的关系。闫洪超等[21]基于Hertz接触理论模拟砂土直剪试验推导出具有普适性的砂土宏细观力学关系的模型,为定量构建土的本构关系奠定了基础。

1.2 复杂应力路径模拟

散体类岩土材料在实际中常经历复杂应力路径作用,常规物理试验无法还原许多客观存在的三维复杂应力路径,离散元数值试验方法能够简便地实现对复杂应力路径的模拟,揭示复杂应力路径下散粒体力学性质的演化规律。薛龙等[22]通过对3个主应力大小和方向的任意控制实现复杂应力路径,显示了数值试验方法在深入研究三维复杂应力路径下粒状介质力学响应方面所拥有的能力和优势。刘嘉英等[23-24]进行真三轴应力路径模拟试验分析了加载过程中宏微观力学参数的演化过程以及组构张量与应力张量的联系。李涛等[25]为探究结构性土粒间胶结的关键特征,建立了结构性土离散元模型开展多种应力路径下的离散元数值分析,再现了结构性土体的主要宏观力学性质。郝晓平[26]研究了砂岩试样在不同应力路径下宏观应力-应变响应及细观接触力的分布规律发现试样破坏时的峰值应力及体积应变与围压有强相关性。

1.3 多尺度耦合分析

为研究环境及地质灾害、第三方施工行为等复杂载荷扰动下基坑、管道、桩体等宏观变形与散粒体颗粒的细观力学行为,需要考虑基坑、管道、水体等与周围土体的耦合作用,鉴于离散元法在土体分层、压缩、运移、流变等大变形领域时有其独特的优势,可采用离散元与其他数值方法耦合这样不仅提高了计算效率,还提升了模拟精度,进而能够更好地应用于岩土工程之中。徐涛龙等[27]利用SPH-FEM算法和PFC-FLAC耦合分析技术对特定管土耦合场景成功实践,丰富了土中多物理场耦合领域。高涛等[28]采用离散元-有限元差分法的跨尺度耦合进行桩和周围土体接触过程中的稳定性分析,详细解释了下沉过程中周围土的细观变形、应力分布和桩体自身的变化及原因。禹海涛等[29]采用有限差分法(finite difference method,DEM-FDM)方法研究隧道开挖诱发近断层错动的物理机制。杨江坤等[30]基于DEM-FDM耦合方法构建了三维分离式霍普金森压杆(split hopkinson pressure bar,SHPB)冲击数值模型为细观角度理解岩石动态冲击破坏提供了新的思路。李伟一等[31]利用CFD-DEM研究间断级配砂土的渗流现象对土体抗剪强度的影响,为从本源上理解砂土渗透侵蚀做出了重要贡献。熊书春等[32]基于CFD-DEM提出了优化算法使用粘结颗粒模型提升了计算效率还解决了三维情况下非球形大颗粒在流场中的运动模拟难题。徐汪豪等[33]使用刚体动力学(rigid body dynamics, RBD)-离散元(DEM)耦合的数值模拟方法,实现了在数值模拟中的滚刀动态被动转动,进一步贴近工程实际。

2 考虑颗粒形状的离散元方法

自然界中的颗粒形状并非绝对的圆球形,因此在颗粒材料离散元数值分析中,已有关于颗粒形状的研究[5-9]。从目前已有的模拟颗粒形状方法来看大致可以分为两大类:一类是模拟颗粒的典型形状,一类是模拟颗粒的真实形状。

2.1 模拟典型颗粒形状

模拟典型颗粒形状是目前的主流模拟方法,其方法为设置特殊形状颗粒替换圆形颗粒以达到模拟典型形状的目的,起初只能较为笼统地表达颗粒的形状,后发展为利用形状参数对颗粒形状进行量化,对颗粒形状量化再与圆形颗粒作对照能简便、高效地探究颗粒形状与宏细观参数之间的内部机理,适用于机制分析。John等[5]首先尝试采用数值模拟方法模拟不同的颗粒形状但仅局限于笼统得表达颗粒形状如较圆、不规则形状等,随后Wadell[3]提出以球形度n(表征椭圆颗粒、多面体颗粒的形貌参数,形貌上越接近球的颗粒,其球形度越接近于1)。来表示颗粒形状与圆球形的偏离程度。在此基础上刘清秉等[34]提出以球形度为关键形状特征参数生成试样。定量分析了球形度对土体力学特性的影响。

宋二祥[7]参考实际工程中堆石料的几何特性将堆石料颗粒形状大致分为4类,利用粘合法将半径相同的球颗粒粘结成密排六方结构并按照单颗粒的精细程度及计算效率分别生成不同球度的试样。分别采用13、23、33颗粒簇作为球形度0.878、0.617和0.446的颗粒模型并保证各组试样的圆度和表面粗糙度相同。颗粒之间采用平行黏结模型无需修正密度。在参考已有研究的基础上采用最大与最小颗粒粒径路径之比为1.25,试样边长与颗粒直径之比在5～9。对不同颗粒形状的试样进行常规三轴试验分析了峰值摩擦角、平均接触力等微观变量的变化规律。

周光军等[35]研究砾石含量以及砾石长宽比对砂-砾石混合物宏观力学特性的影响。颗粒形状的量化方法经历了一系列的发展逐渐变得成熟,但量化参数却缺乏工程统一的应用标准。有不少学者采用长细比[36]、磨圆度[37]、轴向系数[38]、扁平度[39]等参数对颗粒形状进行表征取得了不少研究成果但不同形状量化参数适用的工况缺乏统一的标准,导致同一工况下不同学者研究不同的颗粒形状参数,彼此之间的联系与参考价值有所影响长此以往不利于此领域的发展。

王家全等[40]在保证生成的Clump颗粒体积不变的前提下,在三个大小相同的圆形颗粒中间塞入一个可改变直径D的圆形颗粒控制生成的Cump颗粒的球形度n。图1为Clump颗粒生成示意图,图2为球形度示意图。林呈祥等[41]运用类似方法对月壤颗粒进行量化统计后研究颗粒凹凸度对月壤抗剪强度指标的影响。魏婕等[42]通过Clump命令预生成几种形状的颗粒在0～360°方位角内随机分布来模拟不同形状的颗粒和圆形颗粒之间的区别。杨升等[13]通过Clump命令生成椭圆形颗粒和哑铃型颗粒来分析不同颗粒的配位数和体积变化及空间转动特性也取得了较好的结果。利用Clump方法模拟颗粒典型形状的方法在定量生成所需量化参数具有便捷性,在研究颗粒形状对土体力学特性的影响上直观得出结论,便于机理研究,但实际工程中参数标定困难、计算量大、误差不可控等原因都是实际工程中应用所必须克服的难题。

图1 Clump颗粒生成示意图[40]Fig.1 The schematic diagram of clump generation[40]

图2 Clump颗粒示意图[40]Fig.2 The schematic diagram of particle shape[40]

2.2 模拟颗粒真实形状

颗粒真实形状模拟能够真实反映散粒体实际情况,模拟颗粒真实形状首先需基于扫描设备对颗粒外轮廓进行扫描,其次利用圆形颗粒填充外轮廓和利用数学方程再现颗粒形貌来达到模拟真实形状的目的,这是目前模拟颗粒真实形状的主流方法。

2.2.1 颗粒真实形状扫描

纯利用数学方程或颗粒变形对颗粒真实形状建模十分困难,所以常借助成像设备来对颗粒真实形状重构。二维真实形状模拟常用到显微镜或高清照相机并通过图像处理、数学分析方法计算颗粒真实形状二维形状参数。陈海洋等[43]利用光学显微镜获取钙质砂的二维几何形状成功对钙质砂的几何形态进行描述,研究发现钙质砂颗粒形状具有明显分形特征,与颗粒粒径密切相关。为进一步提高重构形状的准确性,近年来X射线计算机断层扫描技术(computed tomography, CT)成为研究颗粒三维形态特征的主要手段。当X射线穿透物体后,其强度会随物质成分、密度、光源强度等变化,而CT扫描过程则是基于X射线的变化来生成物体的三维影像,将试件置于X射线下旋转360°完成CT扫描过程[44]。早期的三维扫描技术主要用于研究三维剪切破坏过程中颗粒的运动以及接触行为[42]。由于所获得的CT影像体元像素具有锯齿状的特性,近期高精度CT扫描技术也逐渐应用于扫描颗粒形态特征和微观结构[8-9],通过后期图像处理技术可精确计算体积、表面积等三位尺寸[45]。但是边界搜索算法具有较大的误差,进而使得砂粒表面积和表面曲率等形态参数的精确测定成为一个巨大的挑战[46]。崔博等[47]运用三维激光扫描技术提出砾石形状的评定标准,通过三轴压缩模拟试验分析砾石土的力学性质。徐文杰等[48]根据现场试验反演得到土石混合体参数,为离散元法在实际岩土工程中的应用提供了新的思路。

2.2.2 真实颗粒形状建模

Matsushim等[49]和Katagiri等[50]在一定区域内设置数个大小不同的圆形颗粒,让颗粒之间互相运动变形,最后对稳定下来的各个圆形颗粒黏结成为一个整体,这样就完成了对实际形状的模拟。但是由于起始圆的位置具有随机性,所以试验会有一定误差,同时无法量化生成的颗粒形状,总体而言此模拟方法效率较低,工程应用困难。Wang等[51]在此方法的基础上在基本圆相邻较多的区域用一个更大的基本圆来代替相邻的基本圆,提高了模拟效率,此种方法不仅可以模拟平面形状还可以对空间立面进行模拟。薛明华等[52]采用Clump模块生成真实天然砂样模型进行直剪试验模拟揭示颗粒形状对砂性土宏细观力学性质的影响。为了优化计算Ferellec等[53]设置了一系列参数来优化计算效率,但是由于不同环境和材料的优化参数差异较大导致此种方法较难推广应用。Shamsi[54]采用此方法来研究晶粒形状对于颗粒土力学行为的影响。虽然保证了模拟精度但参数选取较为繁琐。李皋等[55]运用填充法对砾岩真实形状建模并开展单轴抗压模拟试验,分析了颗粒形态和颗粒粒径大小对砾岩力学性质的影响,更好地反映了砾岩的力学性质。袁斌等[56]基于等效椭球体改进棱度指数的计算公式可定量生成量化后的颗粒形状,虽参数调整较为繁琐但已经实现了实际中的应用,应用性大幅提升,在应用碎石铁道方面有较大潜力。

徐文杰等[48]基于三维扫描设备利用无嵌套方式将相邻球颗粒互相嵌套在一起来,在块体三维几何空间范围内均匀填充一定比例的直径在1D～1.5D(D为标准圆球颗粒的直径)的颗粒,以块体三维几何网络作为固定边界,采用Yade进行计算并不断调整球体颗粒半径直至所有颗粒的接触力在规定范围之内。基于原位土体直剪试验结果反演得到离散元细观参数,通过数值模拟与现场试验的剪切应力关系曲线验证了此建模方法的可行性和合理性。所构建的块石三维形态数据库极大地优化了土石混合体模拟效率,也为在真实工况下模拟土石混合体颗粒真实形状打下了结实的基础。

边学成等[57]建立计算机成像装置利用三个互相垂直的高清摄像机获取颗粒外轮廓,进行道砟颗粒材料直剪试验。由于依赖三个互相垂直的高清摄像机获取颗粒外轮廓,成像设备的精度会对模拟精度的产生的较大影响,此外仅靠三个互相垂直的成像设备会丢失颗粒表面特性,但误差在允许范围内。Reid等[58]运用LS-DEM提供的技术模拟三轴试验,这种模拟试验的预测水平较高,不仅能够预测剪切带的位置,还能预测剪切带的方向、厚度和倾角。

部分学者运用数字方程进行真实颗粒形状建模,此种方法模拟精度十分灵活,计算效率的优化较为成熟,采用数学方法来模拟颗粒形状无疑为颗粒形状的模拟提供了一种新的思路。Bowman等[59]首次使用数学傅里叶级数描述符来表征二维空间中颗粒的外轮廓,然而,真实颗粒的微观形态远比二维颗粒的外轮廓要复杂得多。Mollon等[60]将二维傅里叶扩展至三维利用三个真实的二维正交平面来构建颗粒的三维表面。然而,由于人工选择原始剖面的三个正交截面,会丢失颗粒的表面特性尤其是十分不规则的颗粒。

Shi等[61]开发的球谐函数分析程序,用以表征颗粒的微观形态,并以离散的方式确定介观效应。此方法可根据所需的精细程度来选择球谐函数的在再现阶数。苏栋等[62]获取颗粒表面的三维几何坐标,再通过球谐函数进行三维重构实现了对骨料颗粒的精确重构和表征分析。Cui等[63]通过球谐函数生成不规则块石,用以研究土石混合体的力学性质。进行土石混合体的崩滑流试验模拟结果与实验测试结果吻合较好。利用球谐函数重构可对粗、细骨料进行模拟在土石混合体上的应用效果较好,计算精度灵活,计算效率优化较为成熟。李景哲等[64]基于高斯差分小波提出了重构骨料形状的方法与球谐重构法相比可有效避免振铃效应导致的重构误差增大的现象。王嗣强等[65]发展了基于水平集方法的任意形态接触算法用以准确计算单元间的接触方向和重叠量,最大程度保留颗粒外轮廓几何特性。

利用初始圆球颗粒运动和变形的方法由于初始圆颗粒的位置和大小具有随机性,需反复调试,变形参数难以参考所以导致此种方法较难推广应用。由于扫描设备成像具有锯齿状特性所以对精度要求较高并且模拟的过程十分耗时,虽有学者提出优化模拟效率的参数,但对于不同的环境和材料还需另外试验选取新的优化参数导致应用性不强。从另一个角度,由于扫描设备的问题这反过来限制了试样的尺寸,总的来说利用三维扫描设备模拟颗粒真实形状的方法还需进一步发展,其中如何平衡模拟精度以及计算效率是其中的关键。

3 存在问题

尽管采用各类方法对不同颗粒形状的岩土类散粒体已有一定的研究,但在考虑颗粒形状进行离散元模拟及分析时,还存在以下问题。

(1)在工程适用性方面,目前缺乏统一的工程应用颗粒形状分类标准,建议在充分认识颗粒形状参数对土力学特性影响的基础上构建出一套行之有效的颗粒形状分类标准。在此基础上进行大规模离散元数值分析,以期合理描述工程问题,采用离散元方法为工程问题提出合理的解决方案。

(2)进行颗粒材料真实形状模拟时,成像设备的扫描精度要求较高,且利用数学函数模拟真实形状对于使用者本身有门槛要求,不便于推广应用。运用数字方程来模拟实际形状的方法对内凹型的颗粒模拟效果不理想。

(3)参数标定一直是离散元模拟试验的一大难点,无论是试错法还是现场试验法参数标定还是十分繁琐。

(4)离散元法模拟真实形状主要集中于块石、砾石、卵石、土石混合体等材料,但真实形状千万种,单一种形状就费时费力,实际工程中一个部位的模拟就需几十甚至上百万个颗粒单元,一次模拟就需数百个小时,若是提升计算效率模拟精度又无法有效保证,因此需要一种兼具效率与精度的广泛适用性方法。

4 结论与展望

离散元法在岩土颗粒材料数值分析中已得到了一定的应用,尤其在机制分析、复杂应力路径模拟和多尺度耦合研究方面。颗粒的形貌特征对土体力学性质产生影响已经得到了学者们的普遍认同,并取得一定的研究成果。

(1)离散元方法在岩土数值分析主要集中于细观机制分析、复杂应力路径模拟、多尺度耦合等方面,其中细观机制分析为离散元岩土工程应用的热点,但大部分还停留在定性分析阶段,不少学者提出的定量分析方法具有一定的参考意义。

(2)模拟颗粒典型形状由于模拟门槛较低,可便捷生成量化后的颗粒形状有适用于机理研究。近年来模拟真实形状已可应用小规模工程,但参数选取困难需、误差不可避免等都是实际工程应用的难点。在提升计算效率方面可以从大数据神经网络出发通过对算法的优化从根本上提升计算效率。

(3)散粒体的加卸载过程在岩土工程中无处不在,岩土颗粒材料在加卸载过程中的力学响应具有一定的特殊性,目前对于散粒体加卸载的研究主要集中于应力路径,研究颗粒形状对加卸载过程的影响仍然较少。在后续的研究中可以对这方面进行深入分析。