基于可变形圆化多边形离散单元法的块体直剪数值试验

Direct Shear Numerical Tests of Blocks Based on the Deformable Spheropolygon-based Discrete Element Method

摘要

本文采用可变形圆化多边形离散单元法模拟块体直剪数值试验，并引入Ⅰ/Ⅱ复合型断裂标准模型，创建一种模拟准脆性材料断裂的方法来模拟直剪过程中块体的破碎效应。对刚性块体、可破碎块体和可变形块体分别开展块体直剪数值试验，分析刚性块体圆化程度对剪切应力、剪胀剪缩特性及法向接触力各向异性系数的影响；通过可破碎块体的应力‒位移关系，对可破碎块体抗拉强度进行敏感性分析，研究抗拉强度对剪切演化规律的影响，对比刚性块体、可破碎块体和可变形块体的剪切特性。直剪试验结果表明：块体圆化程度对块体剪切力学特性影响较大，尤其是切应力，刚性块体的自锁效应随着块体圆化程度的增大而减小；当考虑块体的可变形特性时，随着圆化程度的减小，剪切应力和法向接触力的各向异性系数也会随之增大；相同垂直向荷载作用下，刚性块体、可破碎块体和可变形块体的剪切特性存在较大差异，刚性块体棱角之间存在夹具效应，切应力较大，而可变形块体和可破碎块体内部发生变形，切应力较小。模拟结果说明，基于刚体假定的块体直剪数值试验无法反映块体真实的剪切力学特性，而基于可变形圆化多边形离散单元法的模拟准脆性材料断裂法能够较好地模拟可破碎块体和可变形块体的剪切力学特性。

Abstract

Objective Block materials are widely used in various fields of civil engineering, and conducting an in-depth study of the mechanical properties of block materials is essential. Traditional physical tests provide only limited macroscopic strength and deformation characteristics, and emerging technology tests are costly and complex; therefore, numerical simulation, as an effective alternative method, compensates for these deficiencies. This study uses the deformable spheropolygon-based polygon discrete element method to simulate the numerical test of block direct shear, and the standard model of Ⅰ/Ⅱ mixed-mode fracture is introduced to establish a method for simulating the fracture of quasi-brittle materials, simulating the block-breaking behavior during the direct shear process. Methods Block direct shear numerical tests on rigid blocks, crushable blocks, and deformable blocks were conducted, respectively. First, the numerical model adopted for the direct shear test of the block was presented, along with the calibration of parameters during the simulation. Then, direct shear tests of rigid blocks were simulated using rounded polygons, and the effects of the degree of rounding of rigid blocks on shear stress, shear dilation, and contraction behavior, and the anisotropy coefficient of the normal contact force were analyzed. Then, DSDEM was utilized to simulate the direct shear test of the crushable block, the stress displacement relationship of the crushable block was investigated, and a sensitivity analysis of the tensile strength of the crushable block was conducted to examine the effect of tensile strength on the shear evolution law. Finally, the shear characteristics of rigid blocks, crushable blocks, and deformable blocks were compared. Results and Discussions In the direct shear tests of rigid blocks, under the same vertical load, the shear stress decreased with increasing circularization radius of the block. The smaller the circularization radius was, the more pronounced the block interlocking effect was, and the greater the resulting shear stress was. With increasing vertical load, the difference in shear stress among blocks with different circularization radii increased. Under different vertical loads, when the degree of rounding was small, the vertical displacement of the top plate first decreased and then increased, and this behavior was weakly affected by the vertical load; when the degree of rounding was large, block dilation was suppressed, and only compression occurred. In addition, as the vertical load increased, the vertical displacement of the top plate decreased; as the degree of rounding of the block increased, the block surface became smoother, rotation during the shear process became easier, and overturning of blocks around the contact area occurred. In the direct shear tests of crushable blocks, the crushing rate of nodal units increased sharply at the early stage and showed no significant increase during the middle and later stages. The shear stress of the blocks increased with increasing vertically oriented load, and the stress displacement curves exhibited a softening trend after the shear stress peak under lower vertically oriented loads. During shearing under a vertical load of 0.30 MPa, the number of fractures occurring in the nodal units increased significantly with increasing horizontal displacement and then gradually decreased after the shear stress reached its peak. When examining the shear mechanical properties of blocks with identical shapes but different tensile strengths under the same vertical load of 1.0 MPa, the fracture rate of the nodal units decreased as the tensile strength increased. When comparing the shear characteristics of rigid blocks, crushable blocks, and deformable blocks, the rigid block consistently exhibited the highest shear stress under the same vertical load. The rigid block can not deform under extrusion, resulting in a significant fixture effect between block corners. The crushable block fractured during the shear process, leading to a looser block arrangement and lower shear stress between blocks. The internal deformation of the deformable block reduced the fixture effect between block asperities, and because the Young's modulus of the deformable block was low, the internal jamming effect during direct shear was weak, resulting in low shear stress. Conclusions In the rigid block simulation, a smaller degree of block rounding resulted in a stronger block interlocking effect and higher shear stress. Under larger vertical loads, the degree of block rounding had a more pronounced effect on shear stress, while higher vertical loads inhibited the shear expansion of the specimen and reduced the rate of increase of the anisotropy coefficient of the normal contact force. In the simulation of crushable blocks, shear stress increased with increasing vertical load. Under lower vertical loads, the stress-displacement curve exhibited a strain-softening trend after the shear stress reached its peak, whereas under higher vertical loads, the stress-displacement curve transitioned from strain softening to strain hardening, and no distinct peak was observed. When the tensile strength was low, the degree of block crushing had a greater influence on the shear mechanical properties. In addition, block shear stress increased with increasing tensile strength, and the degree of strain hardening also increased. Under the same vertical load, the shear stress of rigid blocks was the highest, while the shear stresses of deformable and crushable blocks were lower. This result confirms that the direct shear test based on the rigid body assumption has limitations when simulating blocks and cannot accurately characterize the mechanical properties of blocks in such tests.

Graphical abstract

关键词

剪切试验 / 可破碎块体 / 可变形块体 / 圆化多边形 / 离散单元法

Key words

shear tests / breakable blocks / deformable blocks / spheropolygon / discrete element method

本刊网刊

块体材料被广泛应用于土木工程的各个领域，如大体积堆石坝的填筑、铁路路基道砟碎石的铺设、砂质土的基坑开挖等。因此，对块体材料的力学特性进行深入研究具有重要意义。在土力学领域，传统的物理试验，如直剪试验和三轴试验，仅能提供有限的宏观强度与变形特性，细观信息则是通过采用一系列新兴技术进行试验测试而获得，如扫描式电子显微镜^[1]、核磁共振技术^[2]和X射线CT扫描^[3]，但这些试验成本高且过程复杂，数值模拟作为一种有效的替代方法可以弥补此类不足^[4‒5]。

离散单元法（DEM）作为一种代表性方法，已被广泛应用于块体材料的数值试验^[6‒7]，试验结果表明，块体的级配和剪切强度受块体形状影响最为明显^[8]。离散单元法的单元类型较多，其中圆盘因接触检测效率最高^[9]而被广泛采用^[10‒11]。然而，圆盘与块体原始形状相差较大，数值模拟与实际试验结果吻合较差。此外，圆盘不能反映块体之间的互锁性。为更好地模拟块体的几何特性，采用非球形单元模拟块体材料，如椭圆^[12]、球簇模型^[13]、超二次曲面模型^[14]等。然而，非球形单元无法准确模拟四边形块体的几何特性，多边形单元能够更加真实地反映岩石、土壤等块体材料的几何形态，更适用于不规则块体的模拟^[15‒16]。

多边形单元模拟方法规避了椭圆对称性的限制，但缺乏统一的接触力计算模型，单元之间的接触形式比较复杂，如点‒点接触、点‒边接触和边‒边接触等，而且计算效率较低^[17]。因无法计算顶点之间的法向接触力，采用多边形单元模拟点‒点接触时会导致能量不守恒，使计算结果产生偏差。与传统的离散单元法相比，有限‒离散单元耦合分析方法（FDEM）^[18]巧妙地规避了单元之间接触形式的判断，建立了全新的接触计算模型。传统的离散单元法通过嵌入量计算接触力，而FDEM则是采用势函数算法求解接触力^[19]。但FDEM对单元形式有严格的限制，仅适用于三角形单元，形状复杂的单元需离散为若干个三角形，从而导致单元数量增加，计算效率下降^[20]。

颗粒离散元计算效率高，但不能准确地表征不规则块体外形；块体离散元可以较好地表示单元形状，但计算效率低。而圆化多边形离散单元法则兼有颗粒离散元和块体离散元的优点。圆化多边形^[21]由多边形和圆盘的闵可夫斯基和构成，可以准确地描述形状复杂的单元^[22]，且圆化多边形之间的接触近似于圆盘的接触，与块体离散元相比，计算更为简单高效。

圆化多边形离散单元法已应用在诸多领域^[23]，但由于块体假定为刚性，圆化多边形离散单元法不能反映块体的变形，也不能模拟断裂全过程，比如离散单元的变形，裂纹的产生、碎裂以及相互之间的接触。

可变形圆化多边形离散单元法（DSDEM）^[24]将有限单元法与圆化多边形离散单元法相结合，能够模拟块体材料的变形特性。在此基础上引入虚拟裂缝模型（FCM）描述裂缝扩展状态。因此，该方法能够模拟岩石材料的变形特性和裂缝扩展状态，预测断裂的开始、演化过程，以及断裂后离散块体的碰撞。

综合现有研究成果发现，基于刚体假设的块体剪切试验存在局限性，块体的破碎和变形对直剪试验有重要的影响。因此，本文基于可变形圆化多边形离散单元法，建立一种模拟准脆性材料断裂的方法，对块体进行直剪数值试验，研究在不同垂直向荷载作用下，块体圆化程度对刚性块体剪切特性的影响，分析可破碎块体的抗拉强度与其剪切特性的关系。此外，比较不同外力作用下刚性块体、可破碎块体和可变形块体的剪切特性。

1 块体断裂模型

以可变形圆化多边形离散单元法^[24]为基础，将块体离散为有限单元，在每对相邻有限单元之间插入0厚度节理单元，如图1所示。节理单元在有限单元基础上生成，根据剖分的有限元单元之间的公共边界确定节理单元的数量和位置，节理单元断裂后各单元之间的相互作用力则采用DSDEM计算。

断裂模型采用Hillerborg等^[25]提出的虚拟裂缝模型（FCM），节理单元断裂模型如图2所示。图2中，w_{f\mathrm{s}}为材料达到抗剪强度时的切向裂缝张开度，f_{\mathrm{s}}为材料抗剪强度，f_{\mathrm{u}}为剪应力残余值，{\sigma }_{\mathrm{s}}为切应力，{\delta }_{\mathrm{s}}为切向位移，{\sigma }_{\mathrm{n}}为正应力，G_{f\mathrm{I}}为Ⅰ型断裂能，G_{f\mathrm{Ⅱ}}为Ⅱ型断裂能，w_{\mathrm{G}\mathrm{Ⅰ}}为材料达到Ⅰ型断裂能时的法向裂缝张开度，w_{\mathrm{G}\mathrm{Ⅱ}}为材料达到Ⅱ型断裂能时的切向裂缝张开度,k_{\mathrm{n}}为法向罚刚度，f_{\mathrm{t}}为材料抗拉强度，w_{f\mathrm{t}}为材料达到抗拉强度时的法向裂缝张开度，k_{\mathrm{s}}为切向罚刚度，{\delta }_{\mathrm{n}}为法向量位移。

Ⅰ型裂缝中，节理单元的{\sigma }_{\mathrm{n}}为：

（1）

f_{\mathrm{s}}采用具有拉伸截断的摩尔‒库仑破坏准则来定义：

（2）

式中，c为材料黏聚力，\phi为材料内摩擦角。

在Ⅱ型断裂中，{\sigma }_{\mathrm{s}}和{\delta }_{\mathrm{s}}的关系为：

（3）

由于节理单元正应力和切应力的计算具有独立性，通常分别对其拉伸破坏和剪切破坏进行判断。但对于脆性材料，该方法不能准确合理地评估材料的断裂特性。针对以上问题，Hillerborg等^[25]建立了Ⅰ/Ⅱ复合型断裂标准（图2（c））。当断裂发生时，Ⅰ/Ⅱ复合型断裂的法向和切向位移会小于纯Ⅰ型或Ⅱ型断裂的法向和切向位移。有限单元间的最大相对位移被分解为法向位移和切向位移，当法向位移和切向位移满足式（4）条件时，节理单元会发生Ⅰ/Ⅱ复合型断裂破坏：

{\left(\frac{{\delta }_{\mathrm{n}}-w_{f\mathrm{t}}}{w_{\mathrm{G}\mathrm{Ⅰ}}}\right)}^{\mathrm{2}}+{\left(\frac{{\delta }_{\mathrm{s}}-w_{f\mathrm{s}}}{w_{\mathrm{G}\mathrm{Ⅱ}}}\right)}^{\mathrm{2}}\ge \mathrm{1}（4）

节理单元中的应力按照式（5）进行积分计算节点力f_{\mathrm{j}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}}^{}：

f_{\mathrm{j}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}}^{}=\int N^{\mathrm{T}}\sigma \mathrm{d}\Omega（5）

式中，N为节理单元的形函数，\sigma为节理单元应力，\Omega为节理单元的积分域。

2 块体直剪试验及剪切特性分析

2.1 数值模型及参数标定

数值模型中的块体材料采用圆化多边形单元模拟，剪切试验的DEM模型如图3所示。粒径（圆盘半径）D为4～20 mm，圆盘面积等于不规则形状块体的面积。刚性箱体中，试样被模拟为7 192个圆化多边形单元，粒径分布如图4所示。初始状况下，试样由被压缩的松散块体组成，块体具备各向同性。上部箱体在水平和竖直方向上被固定。通过对上部加载板施加恒定垂直向荷载{\sigma }_y，同时使下部箱体以恒定剪切速度v_x水平向移动，对试样进行剪切试验。为模拟可变形块体和可破碎块体的剪切过程，将圆化多边形单元进一步划分为若干有限单元采用DSDEM进行模拟，将上述7 192个块体离散为53 900个四边形有限元网格，块体的有限元离散化过程如图5所示。

为使箱体内部为均匀应力场，试样重力设置为0^[26]。试验中，块体密度\rho =\mathrm{2}\mathrm{000}\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}，与王一伟等^[27]的取值相同，摩擦系数和材料刚度按经验取值，并通过试算不断优化数值模拟中细观参数取值。数值模拟的主要材料力学参数为：杨氏模量E=\mathrm{20.0}\mathrm{G}\mathrm{P}\mathrm{a}，泊松比\upsilon =\mathrm{0.2}，密度\rho =\mathrm{2}\mathrm{000}\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}，法向刚度K_{\mathrm{n}}=\mathrm{0.2}\mathrm{G}\mathrm{P}\mathrm{a}，切向刚度K_{\mathrm{s}}=\mathrm{0.2}\mathrm{G}\mathrm{P}\mathrm{a}，摩擦系数\mu =\mathrm{0.35}，f_{\mathrm{t}}=\mathrm{3.2}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{P}\mathrm{a}，G_{f\mathrm{I}}=\mathrm{120}\mathrm{N}/\mathrm{m}，c=\mathrm{20}\mathrm{M}\mathrm{P}\mathrm{a}，\phi =\mathrm{45}°，G_{f\mathrm{Ⅱ}}=\mathrm{1}\mathrm{500}\mathrm{N}/\mathrm{m}。

2.2 刚性块体直剪试验及剪切特性分析

目前，模拟块体土的非规则形状主要采用圆盘，但圆盘较难模拟块体土的基本几何特性，模拟精度较低。故采用圆化多边形模拟块体土，通过改变圆化半径模拟块体受到的不同程度的圆化，分析块体圆化程度对材料剪切应力、剪胀剪缩特性及法向接触力各向异性的影响。块体圆化半径\alpha越大，表示块体圆化程度越高，未圆化块体的\alpha为0，使用基于势函数的离散单元法对未圆化多边形块体进行数值模拟^[20]。

图6为不同垂直向荷载{\sigma }_y（0.30、0.50、0.75、1.00 MPa）作用下，5种不同\alpha（0、0.1、0.4、0.5、0.9）的圆盘试样的剪切应力与水平向位移的关系。在相同{\sigma }_y作用下，剪切应力随\alpha的增大而减小。试验结果表明：\alpha越小，块体互锁效应越明显，剪切应力越大；随着{\sigma }_y增大，圆化半径不同的块体之间的剪切应力差值增大，说明在垂直向荷载较大的情况下，块体圆化程度对剪切应力的影响较大。

图7为不同{\sigma }_y下顶板垂直向位移和块体水平向位移的关系。由图7可知：当圆化程度较小时，顶板垂直向位移先减小再增大，受垂直向荷载影响小；当圆化程度较大时，块体膨胀受到抑制，只发生压缩；随着垂直向荷载的增加，顶板垂直向位移减小，块体更易发生剪切收缩；随着块体圆化程度增大，块体更加光滑，在剪切过程中更易发生自由旋转，围绕接触的块体发生翻转。因此，与垂直方向相比，块体在水平方向上更易产生位移，增大垂直向荷载相当于增加块体的法向约束，降低了试样的剪胀性。

对于块体材料，在加载负荷过程中，试样的宏观强度和变形与细观演化密切相关。Rothenburg等^[28]提出使用傅里叶函数f_{\mathrm{n}}(\theta )近似描述法向接触力的各向异性的演化规律，其数学表达式为：

f_{\mathrm{n}}(\theta )=f_{\mathrm{0}}[\mathrm{1}+a_{\mathrm{n}}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{2}(\theta -{\theta }_{\mathrm{n}})]（6）

式中，f_{\mathrm{0}}为整体的平均法向接触力，\theta为法向接触力的方向角，{\theta }_{\mathrm{n}}为法向接触力的各向异性的主方向角，a_{\mathrm{n}}为法向接触力的各向异性的系数。

不同垂直向荷载作用下刚性块体法向接触力的各向异性演化，如图8所示。由图8可知：相同外力作用下，圆化程度较小的块体在剪切过程中的法向接触力各向异性系数始终较大；承受较大垂直向荷载的块体，其法向接触力各向异性系数变化幅度较小。试验表明，块体的自锁效应随着块体圆化程度的增大而减小，同时较大的垂直向荷载可以抑制法向接触力各向异性系数的增加速度。

2.3 可破碎块体直剪试验及剪切特性分析

块体材料在应力集中作用下发生剪切破碎，改变原始级配导致地基不均匀沉降进而影响工程安全^[29]。目前研究多集中于块体土宏观力学特性，鲜少研究块体土破碎时的细观力学特性^[30]。

通过改变外力大小，在0.30、0.50、0.75和1.00 MPa 4种垂直向荷载作用下，对材料的抗拉强度进行敏感性分析，研究抗拉强度对剪切演化规律的影响。不同垂直向荷载作用下断裂节理单元的累积分布与水平向位移的关系如图9所示。

由图9可知，在直剪试验初期，节理单元的破碎率急剧增加。随着块体的破碎，试样处于相对疏松的状态，在直剪试验的中后期，块体破碎量没有显著增加。当水平向位移增加至50 mm时，与较小块体相比，较大块体可以接触更多块体，承受更大荷载，更易在直剪过程中发生破碎。0.30 MPa垂直向荷载下可破碎块体的变形示意图如图10所示。可破碎块体的切应力与水平向位移的关系如图11所示。块体的剪切应力随着垂直向荷载的增加而增加，在较低垂直向荷载作用下，应力‒位移曲线在剪切应力达到峰值后呈软化趋势。在较大的垂直向荷载作用下，应力‒位移曲线由应变软化转变为应变硬化，无明显峰值。

节理单元破碎数量与剪切应力的关系如图12所示。由图12可知，在垂直向荷载为0.30 MPa的剪切过程中，随着水平向位移的增大，节理单元发生断裂的数量显著增加，在剪切应力达到峰值后逐渐减少。

为研究块体破碎程度对材料剪切力学特性的影响，对形状相同但抗拉强度不同的块体进行了一系列剪切试验。试验采用3.2、4.5、5.5、6.5和8.5 MPa 5种抗拉强度。在垂直向荷载1.0 MPa作用下，对试样进行剪切试验。断裂节理单元的累积分布与水平向位移的关系如图13所示。由图13可知，在抗拉强度较高的情况下，块体破碎程度对材料剪切力学特性的影响可以忽略不计。随着抗拉强度的增加，节理单元断裂的速率降低；但当抗拉强度较低时，其影响显著，不可忽略。不同抗拉强度下的剪切应力与水平向位移的关系如图14所示。由图14可知，随着部分块体的破碎，块体之间的排列变得疏松，块体之间的接触挤压相对较弱，因此剪切应力增长缓慢。随着抗拉强度的增加，剪切应力增加，块体的应变硬化程度也增加。

2.4 剪切特性对比

采用离散元模拟不可破碎的刚性块体，采用DSDEM模拟可破碎块体和不可破碎的可变形块体，其中，不可破碎的可变形块体的杨氏模量E分别为10.0和15.0 GPa。比较了0.30、0.50、0.75、1.00 MPa这4种不同垂直向荷载作用下刚性块体、可破碎块体和可变形块体的剪切力学特性。

不同{\sigma }_y下块体的剪切应力与水平向位移的关系如图15所示。

由图15可知，在相同垂直向荷载作用下，刚性块体的切应力始终最大。刚性块体被挤压时不能变形，块体的棱角之间存在较大的夹具效应。可破碎块体在剪切过程中被击碎，块体排列松散，因此块体之间的切应力较小。可变形块体的内部变形降低了块体棱角之间的夹具效应，且由于可变形块体的杨氏模量较小，直剪过程中块体内部的卡阻效应较低，切应力较小。本算例的计算结果表明，刚性块体不能准确地表征剪切试验中块体的力学特性，刚体假设存在局限性。

3 结 论

基于可变形圆化多边形离散单元法，建立了一种模拟准脆性材料断裂的方法，对土体多边形块体进行了直剪试验模拟，从宏观和细观两个角度研究其剪切特性，得出以下结论：

1）刚性块体模拟中，块体圆化程度越小，块体互锁效应越强，剪切应力越高。在垂直向荷载较大的情况下，块体圆化程度对切应力影响较大，同时较大的垂直向荷载可以抑制试样的剪胀性和法向接触力各向异性系数的增加速度。

2）可破碎块体模拟中，块体切应力随着垂直向荷载的增加而增加，在较低垂直向荷载作用下，应力‒位移曲线在切应力达到峰值后呈软化趋势，在较高垂直向荷载作用下，应力‒位移曲线由应变软化转为应变硬化，曲线无明显峰值。抗拉强度较低时，块体破碎程度对剪切力学特性影响较大，块体切应力随着抗拉强度的增加而增加，应变硬化程度也随之增加。

3）在相同垂直向荷载作用下，刚性块体切应力最大，可变形块体与可破碎块体切应力较小，验证了基于刚体假设模拟块体的直剪试验存在局限性，不能够准确表征试验中块体的力学特性。

参考文献

[1] Tang Yiqun, Li Jinzhang.Test method and application for microstructures of undisturbed silty sand and sandy silt[J].Environmental Earth Sciences,2018,77(18):657. doi:10.1007/s12665-018-7847-y

[2] Stannarius R.Magnetic resonance imaging of granular materials[J].Review of Scientific Instruments,2017,88(5):051806. doi:10.1063/1.4983135

[3] Yin Yanzhou, Cui Yifei, Liu Dingzhu,et al.Study on microscopic process of fine particle migration in loose soil[J].Advanced Engineering Sciences,2019,51(4):21‒29. doi:10.15961/j.jsuese.201900253

[4] 殷延洲,崔一飞,刘定竺,等.松散土体中细颗粒运移的微观过程研究[J].工程科学与技术,2019,51(4):21‒29. doi:10.15961/j.jsuese.201900253

[5] Tong Lihong, Fu Li, Xu Changjie.Effects of particle rotation on shear strength of granular material system[J].Engineering Mechanics,2024,41(11):125‒133. doi:10.6052/j.issn.1000-4750.2022.09.0808

[6] 童立红,傅力,徐长节.颗粒旋转对颗粒材料系统抗剪强度影响研究[J].工程力学,2024,41(11):125‒133. doi:10.6052/j.issn.1000-4750.2022.09.0808

[7] Cheng Shufan, Zeng Yawu, Gao Rui,et al.Finite‒discrete element model for slip debonding of cement anchors in clay rock[J].Chinese Journal of Geotechnical Engineering,2023,45(12):2594‒2603.

[8] 程树范,曾亚武,高睿,等.黏土岩中水泥锚固体滑脱的有限‒离散元模型初探[J].岩土工程学报,2023,45(12):2594‒2603.

[9] Zhou Xiaowen, Xu Yanbin, Zhao Shiwei,et al.Influence of eccentricity on induced anisotropy of granular media[J].Journal of South China University of Technology(Natural Science Edition),2022,50(11):141‒154. doi:10.12141/j.issn.1000-565X.210747

[10] 周小文,许衍彬,赵仕威,等.偏心率对颗粒介质次生各向异性的影响[J].华南理工大学学报(自然科学版),2022,50(11):141‒154. doi:10.12141/j.issn.1000-565X.210747

[11] Yao Fuhai, Ma Gang, Guan Shaoheng,et al.Interfacial shearing behavior analysis of rockfill using FDEM simulation with irregularly shaped particles[J].International Journal of Geomechanics,2020,20(3):04019193. doi:10.1061/(ASCE)GM.1943-5622.0001590

[12] Gao Huan, Zhai Yue, Wang Tienan,et al.Compressive mechanical properties and strength prediction model of concrete‒granite combined body under active confining pressure[J].Rock and Soil Mechanics,2022,43(11):2983‒2992.

[13] 高欢,翟越,汪铁楠,等.主动围压作用下混凝土‒花岗岩组合体抗压力学特性与强度预测模型研究[J].岩土力学,2022,43(11):2983‒2992.

[14] Xu Kun, Zhou Wei, Ma Gang,et al.Review of particle br-eakage simulation based on DEM[J].Chinese Journal of Geotechnical Engineering,2018,40(5):880‒889. doi:10.11779/CJGE201805013

[15] 徐琨,周伟,马刚,等.基于离散元法的颗粒破碎模拟研究进展[J].岩土工程学报,2018,40(5):880‒889. doi:10.11779/CJGE201805013

[16] Ali U, Kikumoto M, Ciantia M,et al.Systematic effect of particle roundness/angularity on macro-and microscopic behavior of granular materials[J].Granular Matter,2023,25(3):51. doi:10.1007/s10035-023-01341-y

[17] An Huicong, Ouyang Chaojun, Zhou Shu.Dynamic process analysis of the Baige landslide by the combination of DEM and long-period seismic waves[J].Landslides,2021,18(5):1625‒1639. doi:10.1007/s10346-020-01595-0

[18] Pan Shiting, Li Kai, Zhang Chaohui,et al.Influence of coarse aggregate shape on chloride diffusivity in concrete by numerical modelling[J].Materials Review,2022,36(10):85‒93. doi:10.11896/cldb.21030145

[19] 潘诗婷,李凯,张超慧,等.粗骨料形状对混凝土氯离子扩散性能影响的数值模拟研究[J].材料导报,2022,36(10):85‒93. doi:10.11896/cldb.21030145

[20] Wang Yanlin, Guo Qi, Sun Shanshan,et al.Discrete element simulation and analysis of ice-inclined structure interaction[J].Journal of Shanghai Jiaotong University,2022,56(9):1168‒1175. doi:10.16183/j.cnki.jsjtu.2021.182

[21] 王延林,郭麒,孙珊珊,等.冰与斜面结构作用过程的离散元模拟与分析[J].上海交通大学学报,2022,56(9):1168‒1175. doi:10.16183/j.cnki.jsjtu.2021.182

[22] Wang Siqiang, Ji Shunying.Non-linear contact model for super-quadric element considering the equivalent radius of curvature[J].Chinese Journal of Theoretical and Appl-ied Mechanics,2018,50(5):1081‒1092. doi:10.6052/0459-1879-18-103

[23] 王嗣强,季顺迎.考虑等效曲率的超二次曲面单元非线性接触模型[J].力学学报,2018,50(5):1081‒1092. doi:10.6052/0459-1879-18-103

[24] Khazeni A, Mansourpour Z.Influence of non-spherical shape approximation on DEM simulation accuracy by multi-sp-here method[J].Powder Technology,2018,332:265‒278. doi:10.1016/j.powtec.2018.03.030

[25] Liu Lu, Ji Shunying.A new contact detection method for arbitrary dilated polyhedra with potential function in discrete element method[J].International Journal for Numerical Methods in Engineering,2020,121(24):5742‒5765. doi:10.1002/nme.6522

[26] Mousakhani M, Thoeni K, Fityus S,et al.A novel and robust contact detection algorithm for arbitrarily shaped bl-ocks[J].Computers and Geotechnics,2021,140:104430. doi:10.1016/j.compgeo.2021.104430

[27] Munjiza A, Bangash T, John N W M.The combined finite-discrete element method for structural failure and collapse[J].Engineering Fracture Mechanics,2004,71(4/5/6):469‒483. doi:10.1016/s0013-7944(03)00044-4

[28] Zhang Hong, Wu Wei, Guan Xiaofei,et al.A contact detection algorithm of general polyhedrons based on local convex decomposition[J].Chinese Journal of Rock Mechanics and Engineering,2018,37(12):2709‒2720. doi:10.13722/j.cnki.jrme.2018.0853

[29] 张洪,武威,关晓飞,等.基于局部凸分解的一般多面体接触搜索方法[J].岩石力学与工程学报,2018,37(12):2709‒2720. doi:10.13722/j.cnki.jrme.2018.0853

[30] Liu Xunnan, Zhao Lanhao, Mao Jia,et al.Discrete element method using three dimensional distance potential[J].Rock and Soil Mechanics,2018,39(7):2639‒2650. doi:10.16285/j.rsm.2017.1702

[31] 刘勋楠,赵兰浩,毛佳,等.三维距离势离散单元法[J].岩土力学,2018,39(7):2639‒2650. doi:10.16285/j.rsm.2017.1702

[32] Hopkins M A.Polyhedra faster than spheres?[J].Engineering Computations,2014,31(3):567‒583. doi:10.1108/ec-09-2012-0211

[33] Wu Jie, Liu Lu, Ji Shunying.Analysis of wave attenuation and ice breaking characteristics based on DEM‒SPH coupling method[J].Chinese Journal of Hydrodynamics,2023,38(2):176‒186.

[34] 吴捷,刘璐,季顺迎.基于DEM‒SPH耦合方法的冰区波浪衰减及冰层破碎特性分析[J].水动力学研究与进展A辑,2023,38(2):176‒186.

[35] Fu Yanghua, Wang Guojun, Huang Yating,et al.Study on failure modes of interaction between sea ice and vertical structures[J].The Ocean Engineering,2023,41(3):27‒36. doi:10.16483/j.issn.1005-9865.2023.03.003

[36] 付仰华,王国军,黄亚婷,等.海冰与直立结构相互作用的破坏模式研究[J].海洋工程,2023,41(3):27‒36. doi:10.16483/j.issn.1005-9865.2023.03.003

[37] Zhou Yangshiqi, Zhao Lanhao, Shao Linyu,et al.A deformable spheropolygon-based discrete element method[J].Rock and Soil Mechanics,2022,43(7):1961‒1968.

[38] 周洋诗琦,赵兰浩,邵琳玉,等.可变形圆化多边形离散单元法[J].岩土力学,2022,43(7):1961‒1968.

[39] Hillerborg A, Modéer M, Petersson P E.Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[J].Cement and Concrete Research,1976,6(6):773‒781. doi:10.1016/0008-8846(76)90007-7

[40] Zhao Shiwei, Zhou Xiaowen, Liu Wenhui,et al.DEM simulation of direct shear tests considering particle angularity effect[J].Rock and Soil Mechanics,2015,36(S1):602‒608. doi:JournalArticle/5b3c1a5dc095d70f00a6f3b5

[41] 赵仕威,周小文,刘文辉,等.考虑颗粒棱角影响的直剪试验的离散元模拟[J].岩土力学,2015,36(S1):602‒608. doi:JournalArticle/5b3c1a5dc095d70f00a6f3b5

[42] Wang Yiwei, Liu Run, Sun Ruohan,et al.Macroscopic and mesoscopic correlation of granular materials based on rolling resistance linear contact model[J].Rock and Soil Mechanics,2022,43(4):945‒956.

[43] 王一伟,刘润,孙若晗,等.基于抗转模型的颗粒材料宏‒细观关系研究[J].岩土力学,2022,43(4):945‒956.

[44] Rothenburg L, Bathurst R J.Analytical study of induced anisotropy in idealized granular materials[J].Géotechnique,1989,39(4):601‒614. doi:10.1680/geot.1989.39.4.601

[45] Zhao Tiancheng, Wei Yufeng, Xue Yuan,et al.Study on the initial structural effects of particle breakage in coarse-grained soils[J].Journal of Chengdu University of Technology(Science & Technology Edition),2024,51(3):489‒499. doi:10.3969/j.issn.1671-9727.2024.03.11

[46] 赵天丞,魏玉峰,薛元,等.粗粒土初始结构效应对颗粒破碎的影响研究[J].成都理工大学学报(自然科学版),2024,51(3):489‒499. doi:10.3969/j.issn.1671-9727.2024.03.11

[47] Meng Minqiang, Xiao Yang, Sun Zengchun,et al.Recent progress on crushing‒strength‒energy dissipation of coarse gra-nular soil and biocementation at contacts[J].Scientia Sinica(Technologica),2022,52(7):999‒1021.

[48] 孟敏强,肖杨,孙增春,等.粗粒料及粒间微生物胶结的破碎‒强度‒能量耗散研究进展[J].中国科学:技术科学,2022,52(7):999‒1021.

基金资助

国家重点研发计划项目(2022YFC3005402)

水利部重大科技项目（SKS‒2022108）

河海大学水安全与水科学协同创新中心