基于正交响应面法的砂土细观参数标定

摘要：颗粒离散元的模拟精度主要取决于所标定细观参数的准确性，为提升细观参数的标定精度和效率，本文首先采用抗转动线性接触模型模拟砂土三轴压缩试验；再利用正交试验中的方差分析法筛选出宏观参数的控制性细观参数；然后建立宏观参数与控制性细观参数的响应面方程，依据室内三轴试验确定的砂土宏观参数，代入所建的响应面方程求解得到控制性细观参数；最后结合三轴试验结果验证所标定细观参数的准确性。结果表明：宏观参数包括弹性模量、泊松比、峰值摩擦角和剪胀角，其中弹性模量的控制性细观参数为有效模量和摩擦系数，泊松比的控制性细观参数为刚度比与摩擦系数，峰值摩擦角和剪胀角的控制性细观参数都为摩擦系数；通过三轴试验得出宏观参数取值，弹性模量为83.83 MPa、泊松比为0.45、峰值摩擦角为34.47°、剪胀角为8.93°；进一步标定细观参数取值，有效模量为153.35 MPa，刚度比为2.16，摩擦系数为0.45；利用标定的细观参数进行三轴试验模拟，结果与室内试验结果在低围压（100 kPa）下峰值强度相差7.2%，在高围压（300、500 kPa）下误差都控制在了±15.0%以内，证明了所提出的标定方法具有准确性和可靠性。

关键词：砂土；离散元法；细观参数；正交试验；响应面法；宏观参数；弹性模量；泊松比

doi：10.13278/j.cnki.jjuese.20230132

中图分类号：TU4

文献标志码：A

Supported by the Natural Science Foundation of Hebei Province （E2020201017），the Scientific Research Fund of Institute of Engineering Mechanics， China Earthquake Administration （2019EEEVL0202） and the Science and Technology Research Project of Higher Education Institutions in Hebei Province （ZD2020157）

Calibration of Sand Mesoscopic Parameters Based on Orthogonal-Response Surface Method

Liu Hongshuai"" Zhang Dongtao1，2

1. College of Civil Engineering and Architecture， Hebei University，Baoding 071000，Hebei，China

2. Institute of Geotechnical Engineering， Hebei University，Baoding 071000，Hebei，China

Abstract： The simulation accuracy of the particle discrete element method （DEM） mainly depends on the accuracy of calibrated mesoscopic parameters. To improve the calibration accuracy and efficiency of mesoscopic parameters， firstly，the rolling-resistance linear contact model was used to simulate sand triaxial compression tests. Then， the variance analysis method in orthogonal experiments is used to select the controlling mesoscopic parameters （CMP）， and response surface equations are established to relate the macro parameters to the CMP. By substituting the macro parameters determined from triaxial tests into the response surface equations， the values of the CMP are solved. Finally， the accuracy of the calibrated micro parameters is verified by comparing the simulation results with the triaxial test results. The results show that the macro parameters include elastic modulus， Poisson’s ratio， dilatancy angle， and peak friction angle. The CMP for the elastic modulus are effective modulus and friction coefficient， while the controlling micro parameters for Poisson’s ratio are stiffness ratio and friction coefficient. The CMP of the peak friction angle and dilatancy angle are the friction coefficient. The elastic modulus， Poisson’s ratio， peak friction angle and dilatancy angle from the triaxial tests， are 83.83 MPa， 0.45， 34.47° and 8.93°， respectively. The corresponding mesoscopic parameters are determined as follows： Effective modulus of 153.35 MPa， stiffness ratio of 2.16， and friction coefficient of 0.45. Comparing simulation results from the calibrated mesoscopic parameters with the triaxial tests， the difference in peak strength is 7.2% at a low level of confining pressure （100 kPa）， and the errors are within 15.0% at higher level of confining pressures （300 and 500 kPa）.

Key words： sand;" discrete element method; mesoscopic parameter; orthogonal test; response surface method； marco parameter； elastic modulus； Poisson’s ratio

0 引言

砂土作为浅表地层的天然散粒体或建筑材料，广泛应用于各种工程建设中。砂土颗粒间基本无黏结，在地震荷载、降雨和高压差水头等作用下，可能会发生液化、渗透变形、滑坡等非连续大变形破坏，造成大量工程破坏、人员伤亡等严重后果[14]。离散元法基于牛顿第二定律模拟砂土颗粒的运动，引入接触模型处理颗粒单元之间的相互作用[5]，通过颗粒接触行为，得出材料内部受力变形的演化过程，为深入理解砂土液化、冲刷、管涌和应变局部化的力学机理提供了有效方法。离散元的模拟精度主要取决于细观参数标定的准确性（即细观参数能否真实地反映砂土颗粒间的接触力学行为）。因此，如何快速、准确地标定砂土细观模型参数对于研究砂土力学行为及其灾变机理具有重要的理论价值和工程应用意义。

目前，离散元细观模型参数的标定方法主要分为两类：一为“试错法”。即通过大量试算不断修改模型参数，待与室内试验结果相匹配后结束[610]。该方法计算量大、随机性高、耗时长、标定精度较低，且多数研究主要基于偏应力轴向应变关系曲线直接得到最终的细观参数，但未明确参数标定的详细过程。二为“优化法”。可以分为单一法和复合法。其中单一法常用的有正交试验法

和PlackettBurman（PB）试验设计法。此类方法是

将细观参数作为自变量，宏观参数看作因变量，设计不同自变量在不同水平下的正交组合。先对每一组合进行单元试验模拟（如直接剪切试验、三轴压缩试验等），分析各细观参数对宏观参数的影响规律；再通过比较各细观参数对宏观参数影响的相对大小，筛选出各宏观参数的控制性细观参数，从而确定调整哪些细观参数能快速拟合室内试验结果[1115]。虽然正交试验法和PB设计法明确了细观参数对宏观参数的影响规律，在很大程度上缩小了试算范围，但是并不能直接确定控制性细观参数的取值，仍需缩小范围进一步试算。为此，部分学者在单一法的基础上，结合等值线法[16]、随机森林法[17]等提出了复合优化法。

相较于正交试验法和PB设计法，复合优化法在明确宏细观参数之间的敏感性后不需继续试错，而是通过建立数学关系式，直接求出细观参数的具体值，这不仅降低了计算成本，且模拟精度也可以得到保障。

本文基于离散元软件PFC3D，选取抗转动线性接触模型，采用正交试验确定控制性细观参数，运用响应面法构建各宏观参数与控制性细观参数的预测方程，提出一种新的砂土细观参数标定方法——正交响应面复合优化法，并利用三轴试验结果验证所提出方法的有效性。

1 砂土离散元建模及试验指标

1.1 接触模型

砂土离散元模型若考虑颗粒不规则形状，则需用颗粒簇来模拟，会造成巨大的计算代价。已有研究结果[1819]表明，只要能合理考虑颗粒间的转动阻力，球形颗粒亦能很好地模拟不规则颗粒的力学行为；因此，本文选取抗转动线性接触模型模拟砂土的力学行为。

抗转动线性接触模型包括法向、切向和抗转动三种接触类型（图1），其共同点是都包含弹簧及与之并联的阻尼器，弹簧和阻尼器分别反映颗粒接触面的力学行为和能量耗散特性。

法向接触同时还包含分离器，用于模拟法向只受压不受拉的力学行为。通常假定法向压力与颗粒间的法相位移重叠量成正比，若忽略阻尼力，法向力Fn计算公式为：

Fn=KnUn，Un≥0；Fn=0，Unlt;0。（1）

式中：Kn为法向刚度；Un为颗粒间法向重叠量。

切向接触与法向接触不同的是滑片代替了法向接触中的分离器，用于模拟接触点超过剪切强度后的塑性滑动。切向力Fs的计算公式为：

Fs，max=μFn； Fs=Fs0+KsΔUs。（2）

式中：Fs，max为切向力阈值，超过该值后切向力不再增大；μ为摩擦系数；Fs0为计算开始前的切向力；Ks为切向刚度；ΔUs为颗粒间切向位移重叠量的变化量。

抗转动接触模型由抗转动弹簧和阻尼器组成，用于等效弥补圆形颗粒不能模拟真实颗粒的抗转动能力。抗转动力矩Mr的计算公式为：

Mr，max=μrR-Fn；Mr=Mr0+KrΔθr；1R-=1RA+1RB。（3）

式中：Mr，max为抗转动力矩阈值，超过该值后抗转动力矩不再增大；μr为转动阻力系数；R-为等效半径；Mr0为计算开始前的抗转动力矩；Kr为扭转刚度，与接触点的切向刚度有关；Δθr为相对转角的变化量；RA和RB分别为相接触两个球的半径。

1.2 离散元模拟

离散元模拟分为数值制样、固结、加载三个阶段。1）数值制样：以砂土室内三轴压缩试验为研究对象，试样高100 mm，直径49.9 mm，所使用砂土的级配如图2所示（粒径范围为0.1～ 0.5 mm）。考虑到若按实际粒径大小生成试验颗粒，数目过于庞大，计算代价难以接受；为减小计算量，本文采用与土工试验粒径相当、缩小数值试样尺寸的方法。根据相关研究成果，只要满足试验的容器直径（本文为49.9 mm）与最大粒径的比值大于5[2021]，缩小试样直径或扩大颗粒直径所得到的结果即能满足精度要求。本文试验颗粒的最大粒径为0.5 mm，故将土样尺寸缩小为高15 mm、直径7.5 mm，可满足室内三轴试验高径比为2的条件。采用刚性墙模拟三轴橡皮膜，上、下底部各生成一个无限大的正方形墙体模拟加载板。按照试验级配生成9 913个颗粒，颗粒密度设置为2 650 kg/m3，细观接触模型为抗转动

线性模型。在生成试样过程中，有效模量和刚度比取自正交试验设计方案，先将颗粒间、颗粒与墙体的摩擦系数和抗转动系数均设为0，以便生成均匀的数值试样；之后控制墙体应力为10 kPa，改变颗粒间的摩擦系数和抗转动系数，用于生成不同密实度的试样。2）固结：设置固结应力分别为100、300、500 kPa进行双向固结，待达到指定围压后关闭上下墙伺服。3）加载：对固结完成后的试样进行加载，按5％的轴向应变率加载，待轴向应变达到20%后完成加载。在加载过程中记录试验的偏应力、轴向应变和体积应变。为模拟理想的三轴应力状态，整个数值试验忽略了墙体和颗粒间的摩擦力和颗粒重力。

1.3 宏观力学参数确定方法

本文选择工程应用最为广泛的弹性模量E50、泊松比υ、剪胀角和峰值摩擦角θ为试验宏观参数。相关参数的选取规则如图3所示。

各宏观力学参数计算公式如下。

1）弹性模量：E50=σ50εa，50。（4）

2）泊松比：υ=-εr，50εa，50。（5）

3）剪胀角：=arcsinεp，v-2εp，a+εp，v。（6）

4）峰值摩擦角：θ=arcsin3Mp6+Mp。（7）

其中：

Mp=qfpf[22]；（8）

qf=σ1f—σ3f；（9）

pf=σ1f+2σ3f3。（10）

式中：σ50为峰值强度取一半时对应的偏应力；

εa，50为弹性模量所对应的轴向应变；

εr，50为弹性模量所对应的径向应变；

εp，v为峰值强度对应的体积应变；

εp，a为峰值强度对应的轴向应变；

Mp为峰值强度比；σ1f和σ3f分别为极限状态下第一和第三主应力。

2 细观参数正交试验及其方差分析

每个细观参数都会对宏观参数有不同程度的影响，同时，各细观参数之间的交互效应对宏观力学参数的影响也不可忽略。为此，本节采用正交试验中的方差分析法对细观参数进行敏感性分析，给出各宏观参数的控制性细观参数，为建立宏观参数和细观参数的非线性关系提供依据。

2.1 细观参数取值范围初步确定

本文采用的细观本构模型为抗转动线性接触模型，主要包括有效模量、摩擦系数、抗转动系数、刚度比四个细观参数。为减小搜索范围，本论文收集分析了上述细观参数的经验范围值[2324]：有效模量一般在64～300 MPa之间，在满足弹性球体刚度模型的基础上刚度比应该在1.0～1.5之间。参照已有的研究结果，本次正交试验设计取值范围见表1。

2.2 正交试验方案及结果

本节在前述基础上设计了四因素五水平的正交试验，共25组，具体试验方案如表1所示。四因素表示不同的细观参数，五水平表示细观参数的不同取值大小。根据正交试验设计进行三轴压缩试验模拟，结合所给出的宏观参数确定方法，计算出四因素五水平所对应的弹性模量、泊松比、峰值摩擦角和剪胀角。正交试验汇总结果列于表2，可得出各宏观参数范围：弹性模量为24～281 MPa；泊松比为0.35～0.53；峰值摩擦角为25.3°～39.2°；剪胀角为5.3°～16.3°。

2.3 方差分析

方差分析能够同时分析多个细观参数对宏观参数影响的显著性，并能按照影响的重要程度进行排序，从而确定控制宏观力学参数的细观参数（控制性细观参数）。F统计量是方差分析法的重要组成部分，表示组间方差与组内方差的比值。 F统计量通常用伴随概率p来表示，本文取检验显著水平因子α=0.05，当plt;α时认为细观参数对宏观参数影响显著（p值越小表示越显著），否则认为影响不显著，可忽略不计。表2中各细观参数对宏观参数影响的正交试验方差统计结果见图4。分析图4可以得出：弹性模量的控制性细观参数分别为有效模量和摩擦系数，其余两个参数的影响可忽略不计；泊松比的控制性细观参数分别为摩擦系数和刚度比，其余两个细观参数可忽略不计；峰值摩擦角和剪胀角的控制性细观参数只有摩擦系数，其余三个参数可忽略不计。

3 响应面法

响应面法是一种多元统计分析方法，为用于建立因变量和自变量之间非线性关系的数学模型。因此本节采用响应面法建立宏观参数和控制性细观参数间的预测方程。

运用Design-Expert软件，对控制性细观参数和宏观参数进行中心复合设计（CCD）试验，采用二阶模型来近似表示响应。为便于表示，将细观参数有效模量、摩擦系数、抗转动系数和刚度比分别表示为x1、x2、x3和x4，宏观参数弹性模量、泊松比、峰值摩擦角和剪胀角分别表示为y1、y2、y3和y4。

3.1 弹性模量响应曲面

由前述的方差分析结果可知，抗转动系数和刚度比对弹性模量的影响可忽略，因此将这两个参数取值分别设为0.5和1.5。设计13组CCD试验，计算出弹性模量对控制性细观参数（有效模量和摩擦系数）的响应值，根据试验设计模拟三轴压缩试验，结果列于表3。采用二次回归法拟合得到宏观参数弹性模量与细观参数有效模量、摩擦系数间的响应面方程：

y1=109.28742+0.61084x1-338.21678x2-0.22381x1x2-0.0001624x21+258.44835x22。

（11）

运用表3的数据对式（11）进行检验，所得误差均在±10%以内，可见响应面方程具有较高的精度，可满足工程要求。

依据式（11）计算得到了宏观参数弹性模量与其控制性细观参数有效模量、摩擦系数间的响应曲面。图5a为弹性模量的等值线图，可以看出当有效模量一定时，随着摩擦系数增大，弹性模量逐渐减小；当摩擦系数一定时，弹性模量随有效模量的增大而增大。

由图5b可知：弹性模量和有效模量存在较强的正相关性，当有效模量增大时，弹性模量也随着增大；而弹性模量有随摩擦系数增大而减小的趋势，尤其是在摩擦系数小于0.50的情况下，弹性模量的下降速度要比摩擦系数大于0.50时更大，其对弹性模量的影响比大于0.50时更为显著；此外，有效模量方向的变化梯度要明显大于摩擦系数方向，这表明弹性模量的响应值对有效模量更加敏感。该特征在方差分析中的p值大小上也有体现。

3.2 泊松比响应曲面

泊松比的控制性细观参数为摩擦系数和刚度比，非控制性细观参数有效模量取150 MPa，抗转动系数取0.5，设计13组CCD试验计算泊松比对刚度比和摩擦系数的响应值，模拟三轴压缩试验，结果列于表4。

采用二次回归法拟合得到宏观泊松比与细观刚度比、摩擦系数间的响应面方程：

y2=0.2593+0.17943x2+0.078588x4+0.030833x2x4-0.15404x22-0.012173x24。

（12）

用表4的数据对式（12）进行检验，所得误差均在±10%以内，可见响应面方程具有较高的精度，满足工程要求。

依据式（12）的响应面方程计算可得宏观参数泊松比与细观参数摩擦系数、刚度比的响应曲面，见图6。

图6a为泊松比的等值线图，可见：当刚度比一定时，随着摩擦系数增大，泊松比逐渐增大；当摩擦系数一定时，泊松比随刚度比的增大而增大。

如图6b所示：泊松比与摩擦系数和刚度比都呈现较强的正相关性，当摩擦系数小于0.50时，曲线更加陡峭，表明摩擦系数较小时对泊松比的影响较大；刚度比也有类似的规律，当刚度比小于2.00时，泊松比的变化梯度比刚度比大于2.00时更大；此外刚度比方向上的变化梯度要明显大于摩擦系数方向的，说明泊松比对刚度比的变化更加敏感。

3.3 峰值摩擦角、剪胀角与摩擦系数拟合

由于影响峰值摩擦角和剪胀角的控制性细观参数只有摩擦系数，因此只需对摩擦系数取不同值进行重复试验，得出结果后进行多项式拟合即可。将摩擦系数设为0.20、0.40、0.60、0.08和1.00。进行三轴压缩试验模拟得到峰值摩擦角分别为29.82°、33.93°、36.17°，37.41°和38.15°，剪胀角分别为7.57°、8.94°、9.06°、8.07°和7.20°。通过二次曲线拟合得到响应面方程：

y3=25.19576+26.8282x2-13.80355x22;

（13）

y4=5.856+11.17357x2-9.98214x22。

（14）

依据式（13）与式（14）响应面曲线方程分别计算得到峰值摩擦角、剪胀角与摩擦系数的拟合曲线图（图7、8）。由图7可见，峰值摩擦角与摩擦系数呈正相关性，且摩擦系数在0.60以下时斜率更大，可见在摩擦系数较小时，峰值摩擦角对摩擦系数更加敏感。由图8可知，剪胀角随摩擦系数的增大呈现

先增大后减小的趋势，在摩擦系数为0.55时剪胀角

达到峰值。这可能是由于：剪胀角是描述试样剪胀规律的参数，在较小的摩擦系数下，颗粒间的摩擦阻力较小，试样受到剪切时更容易发生滑移，导致剪胀现象增强，剪胀角增大；随着摩擦系数的增大，颗粒间摩擦阻力增加，导致剪胀现象减弱，剪胀角增长速度减小；当摩擦系数进一步增大时，颗粒间摩擦阻力达到一定程度，导致颗粒紧密堆积且发生互锁现象使土体体积变化受到限制，剪胀角开始减小。

4 参数标定及试验验证

首先依据室内试验获取宏观力学参数，将其代入响应面预测方程，求解得到控制性细观参数；再结合已有数值模拟经验，综合确定控制性细观参数取值；然后根据标定出的细观参数进行三轴压缩试验模拟，得出宏观应力应变曲线，并与试验曲线进行对比分析；最后完成细观参数标定的验证工作。

将围压设为100 kPa对福建标准砂进行室内三轴试验，所选福建砂的物理指标如表5所示，图9为室内三轴试验相应的偏应力轴向应变和体积应变轴向应变曲线。

依据图9可分别确定弹性模量、泊松比、峰值摩擦角和剪胀角，结果如表6所示。细观参数中，摩擦系数作为单因素变量分别与峰值摩擦角（式（13））、剪胀角（式（14））成对应关系，代入相关公式计算得到的摩擦系数分别为0.45、0.63，二者差别显著。考虑到砂土为摩擦材料，因此本次标定所使用的摩擦系数取0.45。由式（11）和式（12）可知，摩擦系数确定后，弹性模量与有效模量、泊松比与刚度比成对应关系。综合表6的弹性模量、泊松比和式（11）（12），便可得到有效模量和刚度比。由前述可知，抗转动系数为所有宏观参数的非控制性细观参数，本文采用Gu等[24]建议值，抗转动系数取为0.5。细观参数标定的最终结果列于表7。

根据标定的细观参数建立100、300、500 kPa三种围压条件下的三轴压缩试验模拟，给出应力应变曲线，并与试验结果进行对比，以此来验证所标定细观参数的可靠性。对比模拟结果和试验结果（图10）可看出，两者的总体趋势基本一致，能够较合理地预测应变软化现象。与标定围压相一致的情况下（100 kPa），拟合效果最好，预测峰值强度比试验值低7.2%；随围压的增大，即距标定围压的差别越大，预测峰值强度的误差越大（围压为300 kPa时，峰值预测误差12.1%；围压为500 kPa时，峰值预测误差13.4%），但峰值预测误差均小于15.0%。由此可知，当使用所标定出的细观参数进行不同于标定围压的试验时，其预测的强度峰值误差会变大。因此若模拟围压变化较大时，为提高模拟精度，建议采用模拟围压上下限的试验数据联合标定来保证结果的可靠性。

5 结论与建议

1）弹性模量的控制性细观参数为有效模量和摩擦系数，泊松比的控制性细观参数为摩擦系数和刚度比，峰值摩擦角和剪胀角的控制性细观参数为摩擦系数。

2）所提出的正交响应面细观参数标定法所标定出的细观参数能够较合理地反映砂土的应力应变特性。与标定围压相一致的应力情况下 （100 kPa），预测峰值强度比试验值低7.2%。

3）当使用所标定出的细观参数进行不同于标定围压的试验时，其预测的峰值强度误差会变大，围压为300 kPa时，峰值预测误差12.1%；围压为500 kPa时，峰值预测误差13.4%。

4）若模拟试验围压变化较大时，为提高模拟精度，建议采用模拟围压上下限的试验数据联合标定来保证结果的可靠性。

参考文献（References）：

[1] He S Y，Wang， X L， Fan H B， et al. The Study on Loess Liquefaction in China： A Systematic Review[J]. Natural Hazards， 2020， 103（2）： 16391669.

[2] 苏德辰，杨战兵，孙爱萍，等. 山西云岗石窟侏罗系地震液化砂岩柱的发现及其大地构造意义[J]. 地质学报，2019， 93（8）： 18141830.

Su Dechen， Yang Zhanbing， Sun Aiping， et al. Discovery of the Jurassic Seismic Liquefaction Sandstone Pipes in the Yungang Grottoes， Shanxi Province and Its Tectonic Significance[J]. Acta Geologica Sinica， 2019， 93（8）：" 18141830.

[3] 刘红帅，宋东松，张东涛，等. 干砂小应变动力学参数的弯曲元试验[J].吉林大学学报（地球科学版），2023，53（6）：18261834.

Liu Hongshuai， Song Dongsong， Zhang Dongtao， et al. Small Strain Dynamic Parameters of Dry Sand by Bender Element[J]. Journal of Jilin University （Earth Science Edition）， 2023，53（6）：18261834.

[4] 刘鑫，苗雪青，黄良，等. 5种排水条件下饱和重塑黄土三轴剪切特性[J].吉林大学学报（地球科学版）， 2023， 53（5）：14991509.

Liu Xin，Miao Xueqing，Huang Liang，et al. Triaxial Shear Behavior of Saturated Remolded Loess Subjected to Five Drainage Conditions[J]. Journal of Jilin University （Earth Science Edition）， 2023， 53（5）：14991509.

[5] Cundall P A， Strack O D L. A Discrete Numerical Model for Granular Assemblies[J]. Geotechnic， 1979， 29（1）： 4765.

[6] 邓志颖，范祥，何忠明，等. L形边裂隙单圆孔砂岩单轴压缩力学特性[J]. 吉林大学学报（地球科学版）， 2023， 53（2）： 541554.

Deng Zhiying， Fan Xiang， He Zhongming， et al. Uniaxial Compressive Mechanical Properties of Sandstone with Single Circular Hole and L-Shaped Side Fissures [J]. Journal of Jilin University （Earth Science Edition）， 2023， 53（2）： 541554.

[7] 孙佳政，傅翼，仇雅诗，等. RB模式下砂土非极限主动土压力的离散元模拟与理论研究[J]. 岩土力学，2023， 44（2）： 603614.

Sun Jiazheng， Fu Yi， Qiu Yashi， et al. Discrete Element Simulation and Theoretical Study on Non-Limit Active Earth Pressure of Sand Under RB Mode[J]. Rock and Soil Mechanics， 2023， 44（2）：" 603614.

[8] 叶斌，宋思聪，倪雪倩. 制样方法对砂土液化力学性质影响的离散元模拟[J]. 同济大学学报（自然科学版），2022， 50（7）： 9981008.

Ye Bin， Song Sicong， Ni Xueqian. Effects of Sample Preparation Method on Sand Liquefaction Mechanical Properties with Discrete Element Simulation[J]. Journal of Tongji University （Natural Science）， 2022， 50（7）： 9981008.

[9] 杨文剑，刘怀忠，谢红强，等. 基于直剪试验的堆积碎石类材料细观参数标定方法及模拟验证[J]. 工程科学与技术，2022， 54（3）： 4654.

Yang Wenjian， Liu Huaizhong， Xie Hongqiang， et al. Mesoscopic Parameter Calibration Method of Accumulated Debris Materials Based on Direct Shear Test and Simulation Verification[J]. Advanced Engineering Sciences， 2022， 54（3）： 4654.

[10] 蒋明镜，石安宁，刘俊，等. 结构性砂土力学特性三维离散元分析[J]. 岩土工程学报，2019， 41（增刊2）： 14.

Jiang Mingjing， Shi Anning， Liu Jun， et al. Three-Dimensional Distinct Element Analysis of Mechanical Properties of Structured Sands[J]. Chinese Journal of Geotechnical Engineering， 2019， 41（Sup.2）： 14.

[11] 许江波，曹宝花，余洋林，等. 基于PFC3D的黄土三轴试验细观参数敏感性分析[J]. 工程地质学报，2021， 29（5）： 13421353.

Xu Jiangbo， Cao Baohua， Yu Yanglin， et al. Sensitivity Analysis of Meso-Parameters in Loess Triaxial Test Based on PFC3D[J]. Journal of Engineering Geology， 2021， 29（5）：" 13421353.

[12] 张宝玉，张昌锁，王晨龙，等. PFC2D平节理模型细观参数标定方法[J]. 计算力学学报，2021， 38（5）： 665673.

Zhang Baoyu， Zhang Changsuo， Wang Chenlong， et al. Calibration Method of Meso-Parameters of PFC2D Flat-Joint Model[J]. Chinese Journal of Computational Mechanics， 2021， 38（5）： 665673.

[13] 陈坤，廖良健，朱瑞，等. 煤系土离散元模型宏细观参数标定及其敏感性分析[J]. 南昌工程学院学报，2023， 42（1）： 3944.

Chen Kun， Liao Liangjian， Zhu Rui， et al. Calibration and Sensitivity Analysis of Macro and Meso Parameters of Discrete Element Model for Coal Measure Soil[J]. Journal of Nanchang Institute of Technology， 2023， 42（1）： 3944.

[14] 崔熙灿，张凌凯，王建祥. 高堆石坝砂砾石料的细观参数反演及三轴试验模拟[J]. 农业工程学报，2022， 38（4）： 113122.

Cui Xican， Zhang Lingkai， Wang Jianxiang， et al. Inversion of Meso Parameters and Triaxial Test Simulation of the Gravel Materials for High Rockfill Dam[J]. Transactions of the Chinese Society of Agricultural Engineering， 2022， 38（4）： 113122.

[15]

汪洋，孟祥林，高进，等. 基于PFC2d的平行粘结模型宏细观参数标定研究[J]. 四川建筑，2022， 42（6）： 137140.

Wang Yang， Meng Xianglin， Gao Jin， et al. Research on Calibration of Macro and Micro Parameters of Linear Parallel Bond Model Based on PFC2d[J]. Sichuan Architecture，2022， 42（6）： 137140.

[16] 王晋伟，迟世春，邵晓泉，等. 正交–等值线法在堆石料细观参数标定中的应用[J]. 岩土工程学报，2020， 42（10）： 18671875.

Wang Jinwei， Chi Shichun， Shao Xiaoquan， et al. Application of Orthogonal-Contour Method in Calibration of Microscopic Parameters of Rockfill Materials[J]. Chinese Journal of Geotechnical Engineering， 2020， 42（10）： 18671875.

[17] 马少坤，黄骁，韦榕宽，等. 考虑滚动阻抗的线性接触模型离散元宏细观参数敏感性研究[J]. 中国安全生产科学技术，2021， 17（6）： 104110.

Ma Shaokun， Huang Xiao， Wei Rongkuan， et al. Study on Sensitivity of Macro and Micro Parameters of Discrete Element in Linear Contact Model Considering Rolling Resistance[J]. Journal of Safety Science and Technology， 2021， 17（6）： 104110.

[18] Wensrich C M， Katterfeld A. Rolling Friction as a Technique for Modelling Particle Shape in DEM[J]. Powder Technology， 2012，217： 409417.

[19] Jiang M J， Shen Z F， Wang J F. A Novel Three-Dimensional Contact Model for Granulates Incorporating Rolling and Twisting Resistances[J]. Computers and Geotechnics， 2015， 65： 147163.

[20] 刘海涛，程晓辉. 粗粒土尺寸效应的离散元分析[J]. 岩土力学，2009， 30（增刊1）： 287292.

Liu Haitao， Cheng Xiaohui. Discrete Element Analysis for Size Effects of Coarse-Grained Soils[J]. Rock and Soil Mechanics， 2009， 30（Sup.1）： 287292.

[21] 尹小涛，郑亚娜，马双科. 基于颗粒流数值试验的岩土材料内尺度比研究[J]. 岩土力学，2011， 32（4）： 12111215.

Yin Xiaotao， Zheng Yana， Ma Shuangke. Study of Inner Scale Ratio of Rock and Soil Material Based on Numerical Test of Particle Flow Code[J]. Rock and Soil Mechanics， 2011， 32（4）： 12111215.

[22] 童森杰，黄茂松，时振昊，等. 基于初始状态参数的砂土峰值摩擦角和极限强度[J]. 岩石力学与工程学报，2022， 41（2）： 389398.

Tong Senjie， Huang Maosong， Shi Zhenhao， et al. Peak Frictional Angle and Ultimate Drainage Strength of Sand Based on Initial State Parameter[J]. Chinese Journal of Rock Mechanics and Engineering， 2022， 41（2）： 389398.

[23] Goldenberg C， Goldhirsch I. Friction Enhances Elasticity in Granular Solids[J]. Nature， 2005， 435： 188191.

[24] Gu X， Zhang J， Huang X. DEM Analysis of Monotonic and Cyclic Behaviors of Sand Based on Critical State Soil Mechanics Framework[J]. Computers and Geotechnics， 2020， 128： 103787.