0.   引言

随着社会经济的快速发展以及西部水电开发进程的加快,中国西南地区正在或即将建设一批调节性能好的高堆石坝。作为堆石坝的主要筑坝材料,堆石料具有非连续性、非线性、非均匀性、各向异性等复杂性质,这些特征很难从连续介质力学的角度来理解^[1]。离散单元法在描述散粒体材料特性方面具有独特优势,目前已广泛应用于岩土介质力学特性的研究^[2],而岩土介质行为的真实重现很大程度上取决于接触模型的选取与细观参数的确定^[3]。因此,快速、准确的确定堆石料离散元模型细观参数,对研究堆石料力学特性、确保堆石坝工程安全具有重要的意义。

目前堆石料离散元模型细观参数的标定方法主要有两类,第一类为试错法^[4-10],即通过不断的试算、调整确定最终细观参数,但该方法计算量大、偶然性高、耗时长,而且大多数研究主要基于堆石料偏应力-轴变曲线直接给出最终的细观参数,并未说明参数标定的具体过程。第二类为智能法^[11],即在大量细观参数敏感性分析的基础上,通过智能优化算法寻求最优的参数组合使数值模拟结果与室内试验结果最接近。李守巨等^[12]针对堆石料室内三轴压缩试验结果,利用神经网络对堆石料的力学模型参数进行反演,估计了堆石料模型参数。马刚^[13]对粒群算法（POS）进行改进并将该算法与径向基函数（RBF）神经网络相结合,对堆石料细观参数进行搜索寻优标定了细观参数。Li等^[14]根据堆石材料三轴压缩试验数据,基于响应面法反演标定了堆石材料细观参数。李守巨等^[15]基于室内试验结果,采用响应面法对堆石料颗粒流模型细观参数进行了反演分析。杨杰等^[16]基于量子遗传算法（QGA）和支持向量机（SVM）,建立了堆石料细观参数标定模型。智能算法较试错法有长足改进,也取得了丰富成果,但该方法需要大量的分析细观参数敏感性以及训练较多的样本数,计算量仍然较大。

针对堆石料细观参数标定计算量大的现状,本文将正交试验设计和等值线法相结合,提出并利用“正交-等值线”法快速标定了堆石料细观参数。首先采用多次正交试验优化了细观参数取值范围；然后在分析宏观响应对细观参数敏感性的基础上,抓住主要影响因素,将宏观指标与细观参数进行组合并确定了细观参数的标定顺序；最后利用各宏观指标的等值线图确定了堆石料细观参数,并与不同围压下的室内试验结果进行对比验证。

1.   堆石料离散元模型

颗粒离散单元中,颗粒之间的相互接触关系即为材料的接触本构特性^[17]。堆石料棱角较多,颗粒形状极不规则,颗粒间咬合力大,且易破碎。目前针对颗粒形状不规则材料,许多学者^[18-20]选用滚动接触模型考虑颗粒形状的影响,因此本文采用PFC软件内置的滚动阻力模型模拟堆石料。本次数值模拟对象为邵磊^[21]博士论文中花岗岩堆石料室内排水三轴剪切试验结果,相对密度2.68,制样干密度2.09 g/cm³,颗粒级配曲线如图1所示,试验结果见图2。本文以围压800 kPa下室内试验结果开展细观参数标定。

图  1  颗粒级配曲线

Figure  1.  Grain-size distribution curves

下载: 全尺寸图片 幻灯片

图  2  室内三轴试验结果^[21]

Figure  2.  Results of laboratory triaxial shear tests^[21]

下载: 全尺寸图片 幻灯片

由于堆石料含有较多的小颗粒,若模拟完整级配曲线,总的颗粒数目将会特别庞大,远远超出当前计算机的运算能力,目前通用做法是将级配进行截断,且截断处理后的数值试样孔隙率一般较大,难以达到室内试样真实孔隙率^[4-10]。文献[7～9]研究了不同截断粒径对试验结果精度的影响,张宜^[9]指出截断粒径取15 mm时,其对宏观强度和变形特性的影响可以接受。借鉴以上研究,本文模拟选取的截断粒径为10 mm,粒径小于10 mm的颗粒用粒径10 mm的颗粒等量替换（如图1所示）,DEM试样为三维模型,初始孔隙率0.34,尺寸为直径$\Phi$=300 mm,高度H=650 mm。

室内通常采用分层填筑,分层振捣击实的方法制备三轴试样,为使数值试样贴近原试样,本文数值试样分6层碾压制备。蒋明镜^[2]提出分层欠压法对颗粒进行压实来获得比较均匀的试样,该方法为了减小连续层压实能量的传递,每层压实后的平均孔隙率略大于整个试样的目标孔隙率,将这一原理应用于每一层,直到最后一层压实后,所有层的平均孔隙率等于目标孔隙率。借鉴该方法制备试样,首先根据截断后的级配、三轴试样总体积以及试样孔隙率求出堆石料各粒组体积,在层高3倍范围内生成各粒组1/6体积所对应的颗粒；然后使其在重力作用下沉降,沉积稳定后用墙体匀速将颗粒压到H/6高度,检查并调整该层颗粒的粒径使之满足指定的级配和孔隙率；采用相同方法依次制备第2层到第6层试样。制样过程示意图见图3。制样过程中颗粒与颗粒（或墙体）间的摩擦系数为0,制样完成后施加预固结围压30 kPa,摩擦系数恢复为初始值。

图  3  制样过程示意图

Figure  3.  Scheme of preparation process of specimens

下载: 全尺寸图片 幻灯片

颗粒密度取室内堆石料密度2680 kg/m³,颗粒-墙体间的法向刚度与颗粒间的法向刚度相同,颗粒与墙体间的摩擦系数、滚动阻力系数、切向刚度均取为0,这样消除加载墙体对试样的摩擦作用,因此施加的应力始终保持在墙的法线方向。通过伺服控制对试样进行等围压固结形成固结试样,轴向加载采用应变控制方式,上下加载板同时向中间压缩,加载速率设定为0.005 m/s,此时颗粒的惯性效应很小^[22],为准静态加载,剪切时通过伺服保证围压偏差在1%以内^[10]。

2.   细观参数范围优化

由于细观参数的取值范围没有统一规定,假设对文中花岗岩堆石料的性质知道较少,因此在标定细观参数之前,先采用多次正交设计方法缩小细观参数的取值范围,从而可以确定细观参数敏感性分析时的取值范围并减轻后续参数标定的工作量。具体遵循原则：初次正交试验时细观参数取值范围略大些,保证室内试验结果包含于数值模拟结果范围内（否则重新选取范围）,然后通过多次正交试验优化细观参数的取值范围,直至前后2次细观参数取值范围差距满足期望。

2.1   正交试验设计

（1）试验指标及试验因素

正交试验的试验指标即为标定细观参数时选取的堆石料宏观特性指标。宏观特性指标选取遵循原则：①选定的指标能方便且容易从试验中获得；②选定的指标需要能较好描述所标定曲线的形态。如图4所示,描述偏应力-轴变关系曲线的指标有割线模量$E_{50}$,峰值强度$S_{\text{P}}$,峰后强度$S_{\text{r}}$等,描述体变-轴变关系曲线的指标有泊松比$\nu$,最大压缩量${\epsilon }_{\text{vmax}}$、剪胀角$\psi$等,具体采用指标需根据标定曲线形态和关注重点确定。以选取体变-轴变关系曲线的宏观指标为例,若随着轴向应变的增加,体变基本一直以压缩段为主（没有或者相变不明显）,建议选用$\nu$,${\epsilon }_{\text{vmax}}$这2个指标；若随着轴向应变的增加存在明显相变,试样体积达到最大压缩量（体积达到最小）后又开始反向增加甚至最终表现为剪胀现象,建议选用$\nu$,${\epsilon }_{\text{vmax}}$和$\psi$这3个指标或者部分指标。参考文献[13, 18, 19]并结合本次标定的室内堆石料试验结果曲线,本文选取的宏观指标为割线模量$E_{50}$,最大压缩量${\epsilon }_{\text{vmax}}$,峰值强度$S_{\text{P}}$,剪胀角$\psi$。$E_{50}$为偏应力-轴向应变曲线上50%峰值偏应力点的割线斜率,$E_{50}$和$\psi$的计算公式如下：

图  4  堆石料常规三轴剪切试验典型轴变-应力-体变关系^[23]

Figure  4.  Stress-strain-volume relation of rockfill materials in conventional triaxial tests^[23]

下载: 全尺寸图片 幻灯片

式中,${\epsilon }_{100\%}$和${\epsilon }_{50\%}$分别为峰值偏应力点和50%峰值偏应力点对应的轴向应变。

正交试验的试验因素即为离散元模型细观参数。模型的细观参数有4个,即法向刚度$k_{\text{n}}$,切向刚度$k_{\text{s}}$,滑动摩擦系数$\mu$,滚动阻力系数${\mu }_{\text{r}}$。为了方便,用有效模量$E^{*}$和法向切向刚度比$k_{\text{r}}$代替法向刚度$k_{\text{n}}$和切向刚度$k_{\text{s}}$。转换关系如下^[20]：

式中,$A=\text{π}{R}^2$（三维情况）,$A=2R$（二维情况）。球体-墙体接触时,R为球体的半径,$L=R$；球体-球体接触时,$R=\min \{R_1,R_2\}$,$L=R_1+R_2$,其中R₁和R₂分别为接触两球体的半径。

（2）试验方案及结果

正交试验首先需要确定各因素的水平,即各细观参数的取值。针对无黏性颗粒材料,Rowe测定的颗粒间摩擦角范围为22°～30°^[1],李广信^[24]总结的滑动摩擦角范围为12.4°～37.6°。对于$k_{\text{r}}$,蒋明镜等^[25]根据散粒体弹性理论给出$k_{\text{r}}$在1.5附近,建议取1.5,也有一些学者进行离散元模拟时取其它值^[4-10]。${\mu }_{\text{r}}$表征颗粒形状的影响,参考相关研究^[18-20]暂定取值范围0～1。$E^{*}$的取值主要取决于试样实际变形模量,且受数值试样疏密程度的影响。参考以上研究,设计初次正交试验细观参数范围及水平,如表1所示,选取4个水平,采用正交设计表L₁₆（4⁵）,开展了16次试验,表2为初次正交试验方案及宏观指标结果。由于本次数值模拟没有考虑颗粒的破碎,所以表2中剪胀角的模拟值均大于室内试验值,但其它宏观指标的室内试验值包含于模拟值范围内,故初选的细观参数范围合适。

表  1  初次正交试验设计表

Table  1.  Design of initial orthogonal tests

因素	取值范围	因素水平
		1	2	3	4
${\mu }_{\text{r}}$	0.3~1.0	0.3	0.5	0.7	1.0
$\mu$	0.3~1.0	0.3	0.5	0.7	1.0
$E^{*}$/100 MPa	1.0~10.0	1.0	4.0	7.0	10
$k_{\text{r}}$	1.0~2.5	1.0	1.5	2	2.5

下载: 导出CSV 

| 显示表格

表  2  初次正交试验方案及宏观指标结果

Table  2.  Schemes of initial orthogonal tests of micro-parameters and results of macroscopic indexes

方案	细观参数（自变量）				宏观指标（因变量）
		${\mu }_{\text{r}}$	$\mu$	$E^{*}$/100 MPa	$k_{\text{r}}$	$E_{50}$/100 MPa	$S_{\text{p}}$/MPa	${\epsilon }_{\text{vmax}}$/%	$\psi$/(°)
数值模拟	1	0.3	0.3	1.0	1.0	0.476	1.779	2.672	9.098
	2	0.3	0.5	4.0	1.5	1.391	2.897	1.164	14.142
	3	0.3	0.7	7.0	2.0	2.106	3.559	0.749	18.761
	4	0.3	1.0	10	2.5	2.691	4.219	0.540	22.030
	5	0.5	0.3	7.0	2.5	1.489	2.030	0.608	6.629
	6	0.5	0.5	10	2.0	2.589	3.562	0.553	15.871
	7	0.5	0.7	1.0	1.5	0.493	3.844	4.428	16.315
	8	0.5	1.0	4.0	1.0	1.833	5.928	1.896	21.954
	9	0.7	0.3	10	1.5	2.010	2.204	0.487	7.630
	10	0.7	0.5	7.0	1.0	2.413	4.067	0.889	16.238
	11	0.7	0.7	4.0	2.5	1.242	4.425	1.290	19.828
	12	0.7	1.0	1.0	2.0	0.422	4.488	4.629	16.151
	13	1.0	0.3	4.0	2.0	1.064	2.066	0.975	9.015
	14	1.0	0.5	1.0	2.5	0.376	2.723	3.211	12.711
	15	1.0	0.7	10	1.0	3.589	6.275	0.824	22.760
	16	1.0	1.0	7.0	1.5	2.536	7.487	1.147	28.529
室内试验	—	—	—	—	—	1.625	3.1926	0.912	2.031

下载: 导出CSV 

| 显示表格

2.2   范围优化

（1）第一次优化

采用极差分析法对表2试验数据进行分析^[26],结果如表3所示。表3中$k_i$为表2任一列因素水平取i（本次正交试验i=1,2,3或4）时所得试验结果的算数平均值。正交试验确定的优方案^[26]指所做试验范围内,满足期望的各因素较优水平组合,若指标越大越好,则选取使指标结果大的水平,反之相反。由于宏观指标模拟值很难与室内试验值正好相等,故我们以模拟值最接近室内试验值为目标,将表3中$k_i$与室内试验值进行比较得到宏观指标的优方案,列于表4。

表  3  初次正交试验极差分析结果

Table  3.  Range analysis results of initial orthogonal tests

宏观指标		因素
		${\mu }_{\text{r}}$	$\mu$	$E^{*}$/100 MPa	$k_{\text{r}}$
$E_{50}$ /100 MPa	k1	1.666	1.260	0.442	2.078
	k2	1.601	1.692	1.383	1.608
	k3	1.522	1.857	2.136	1.545
	k4	1.891	1.870	2.720	1.449
	极差R	0.369	0.610	2.278	0.629
$S_{\text{p}}$ /MPa	k1	3.114	2.020	3.209	4.512
	k2	3.841	3.312	3.829	4.108
	k3	3.796	4.526	4.286	3.349
	k4	4.638	5.530	4.065	3.419
	极差R	1.524	3.510	1.077	1.163
${\epsilon }_{\text{vmax}}$ /%	k1	1.281	1.186	3.735	1.412
	k2	1.871	1.454	1.331	1.807
	k3	1.824	1.823	0.848	1.727
	k4	1.539	2.053	0.601	1.570
	极差R	0.590	0.867	3.134	0.395
$\psi$ /(°)	k1	16.008	8.093	13.569	17.513
	k2	15.192	14.741	16.235	16.654
	k3	14.962	19.416	17.539	14.949
	k4	18.254	22.166	17.073	15.299
	极差R	3.292	14.073	3.970	2.564

下载: 导出CSV 

| 显示表格

表  4  室内试验宏观指标值及初次正交试验结果优方案

Table  4.  Values of macro-indexes and optimal combination of micro-parameters

宏观指标	初次正交试验细观参数组合
E50/100 MPa	${\mu }_r=0.3$$\mu =0.5$ E*=400 MPa $k_{\text{r}}=1.5$
Sp/MPa	${\mu }_r=0.3$$\mu =0.5$ E*=100MPa $k_{\text{r}}=2.0$
${\epsilon }_{\text{vmax}}$/%	${\mu }_r=0.3$$\mu =0.3$ E*=700 MPa $k_{\text{r}}=1.0$
$\psi$/(°)	${\mu }_r=0.7$ $\mu =0.3$ E*=100 MPa $k_{\text{r}}=2.0$

下载: 导出CSV 

| 显示表格

对比表4中各宏观指标的优方案可知,各宏观指标对应的细观参数组合取值不同,而最终的细观参数需要同时兼顾所有宏观指标,因此以所有宏观指标同一细观参数取值为边界得到新的细观参数范围,这样完成了细观参数的范围优化,结果列于表5。从表5可以看出,第一次优化后各细观参数的范围较初始范围均有所缩小,其中滑动摩擦系数$\mu$的变化最明显。

表  5  细观参数范围优化结果

Table  5.  Results of micro-parameters range

细观参数	初始范围	第一次优化	第二次优化
${\mu }_{\text{r}}$	0.3~1.0	0.3~0.7	0.3~0.7
$\mu$	0.3~1.0	0.3~0.5	0.3~0.5
$E^{*}$/100 MPa	1.0~10.0	1.0~7.0	2.5~7.0
$k_{\text{r}}$	1.0~2.5	1.0~2.0	1.0~2.0

下载: 导出CSV 

| 显示表格

（2）第二次优化

在第一次优化基础上,为使结果更加精确,细观参数选取5个水平,采用正交设计表L₂₅（5⁶）开展了25次试验。表6为第二次正交试验方案及结果,类似第一次优化,求得第二次优化后的细观参数范围,列于表5。对比两次优化后的细观参数范围发现,除$E^{*}$有较小幅度的变化外,其余细观参数范围没有变化,说明此时细观参数的范围基本满足期望,所以不再进一步优化。

表  6  第二次正交试验方案及结果

Table  6.  Schemes of micro-parameters and results of macroscopic indexes in second orthogonal tests

方案	细观参数（自变量）				宏观指标（因变量）
	${\mu }_{\text{r}}$	$\mu$	$E^{*}$/100 MPa	$k_{\text{r}}$	$E_{50}$/100 MPa	$S_{\text{P}}$/MPa	${\epsilon }_{\text{vmax}}$/%	$\psi$/(°)
1	0.30	0.30	1.00	1.00	0.354	1.779	2.672	9.098
2	0.30	0.35	2.50	1.25	0.744	2.147	1.498	10.127
3	0.30	0.40	4.00	1.50	1.117	2.439	1.070	10.241
4	0.30	0.45	5.50	1.75	1.466	2.702	0.841	12.362
5	0.30	0.50	7.00	2.00	1.791	2.910	0.692	14.338
6	0.40	0.30	2.50	1.50	0.658	1.945	1.411	8.587
7	0.40	0.35	4.00	1.75	1.009	2.300	1.034	8.744
8	0.40	0.40	5.50	2.00	1.347	2.644	0.825	10.983
9	0.40	0.45	7.00	1.00	1.984	3.135	0.776	13.609
10	0.40	0.50	1.00	1.25	0.434	2.933	3.714	13.212
11	0.50	0.30	4.00	2.00	0.901	1.997	0.970	9.272
12	0.50	0.35	5.50	1.00	1.409	2.547	0.857	9.205
13	0.50	0.40	7.00	1.25	1.787	2.935	0.738	11.786
14	0.50	0.45	1.00	1.50	0.402	2.689	3.418	12.588
15	0.50	0.50	2.50	1.75	0.806	3.253	1.804	15.025
16	0.60	0.30	5.50	1.25	1.245	2.131	0.794	9.845
17	0.60	0.35	7.00	1.50	1.593	2.594	0.682	10.329
18	0.60	0.40	1.00	1.75	0.370	2.374	3.103	10.133
19	0.60	0.45	2.50	2.00	0.749	2.967	1.676	12.623
20	0.60	0.50	4.00	1.00	1.321	3.783	1.402	16.838
21	0.70	0.30	7.00	1.75	1.426	2.136	0.638	7.004
22	0.70	0.35	1.00	2.00	0.340	2.043	2.795	10.135
23	0.70	0.40	2.50	1.00	0.808	2.934	1.767	6.006
24	0.70	0.45	4.00	1.25	1.202	3.396	1.275	13.891
25	0.70	0.50	5.50	1.50	1.560	3.811	1.008	16.327

下载: 导出CSV 

| 显示表格

综上,基于正交试验“试验次数少、代表性较好”的优势,通过较少的试验次数就能缩小所需细观参数的取值范围,故多次正交试验是一种很好的优化细观参数取值范围的途径,特别是如果对需要标定材料的细观参数没有经验可依、知之较少时,此方法的优越性将更加明显。

3.   细观参数敏感性分析

相比较单一因素敏感性分析,正交试验可以突破控制变量法固定某些参数的局限分析细观参数的敏感性^[27]。正交试验极差分析结果中,极差代表各细观参数水平值的改变对宏观指标的影响程度,计算公式为$R=\max \left\{k_i\right\}-\min \left\{k_i\right\}$,极差越大说明细观参数水平值的改变对宏观指标影响越大。为进一步对比各细观参数对宏观指标影响程度的相对大小,此处采用相对极差,其计算公式为$R_i/max\left\{R_1,R_2,R_3,R_4\right\}$,其中$R_i$为某一细观参数水平值的改变对应的宏观指标极差。图5（a）～（d）为表6中各宏观指标的相对极差分析结果。

图  5  宏观指标的极差分析结果

Figure  5.  Range analysis results of macroscopic indexes

下载: 全尺寸图片 幻灯片

分析图5（a）～（d）可发现,各细观参数对宏观指标影响程度的大小排序为,$E_{50}$：$E^{*}>\mu >k_{\text{r}}>{\mu }_{\text{r}}$；$S_{\text{p}}$：$\mu >{\mu }_{\text{r}}>E^{*}>k_{\text{r}}$；${\epsilon }_{\text{vmax}}$：$E^{*}>\mu >k_{\text{r}}>{\mu }_{\text{r}}$；$\psi$：$\mu >E^{*}>{\mu }_{\text{r}}>k_{\text{r}}$,可见各细观参数对不同宏观指标的影响程度不完全相同。

综合对比图5（a）～（d）发现,$E_{50}$和${\epsilon }_{\text{vmax}}$受细观参数影响程度大小相似,均是$E^{*}$的影响相比较于其它细观参数的影响绝对占优,$\mu$的影响程度稍大于$k_{\text{r}}$,${\mu }_{\text{r}}$的影响最小。$S_{\text{P}}$和$\psi$受细观参数的影响程度大小相似,均是$\mu$的影响最大,${\mu }_{\text{r}}$和$E^{*}$的影响程度差不多,$k_{\text{r}}$的影响最小。

4.   细观参数标定

4.1   标定方法

通过敏感性分析可知,每一个宏观指标受多个细观参数不同程度的影响,且标定的细观参数需同时满足所有宏观指标,这导致细观参数的标定更加困难。因此本文抓住主要影响因素,忽略次要影响因素,并借鉴响应面法和等值线的思想^[28],基于已有的宏观指标数据点,通过插值可以获取该区域内其它任意数值的等值线,进而求取期望的细观参数。标定步骤如下：

（1）上节敏感性分析可知,在优化后的细观参数范围内,$E^{*}$对$E_{50}$和${\epsilon }_{\text{vmax}}$的影响绝对占优,$k_{\text{r}}$,${\mu }_{\text{r}}$,$\mu$的影响可以忽略。$E^{*}$,$k_{\text{r}}$,${\mu }_{\text{r}}$,$\mu$这4个细观参数对$S_{\text{P}}$和$\psi$均有一定影响。综合考虑,本文先基于$E_{50}$和${\epsilon }_{\text{vmax}}$的值标定细观参数$E^{*}$和$k_{\text{r}}$。选用表6中不同$E^{*}$和$k_{\text{r}}$组合下$E_{50}$和${\epsilon }_{\text{vmax}}$的数据点,绘制一系列$E_{50}$和${\epsilon }_{\text{vmax}}$等值线图,然后通过数值为室内试验值的$E_{50}$和${\epsilon }_{\text{vmax}}$等值线的交点确定$E^{*}$和$k_{\text{r}}$（详见4.2节）。

（2）固定已经确定的$E^{*}$和$k_{\text{r}}$不变,令${\mu }_{\text{r}}$和$\mu$在优化后的细观参数范围内分别取不同的值,开展一系列数值模拟。同样地方法,以$S_{\text{P}}$和$\psi$室内试验值等值线的交点确定细观参数${\mu }_{\text{r}}$和$\mu$（详见4.3节）。

4.2   确定有效模量$E^{*}$和法向切向刚度比$k_{\text{r}}$

图6为基于表6中的数据绘制的$E_{50}$（黑色线）和${\epsilon }_{\text{vmax}}$（灰色线）等值线图。从图中可以看出,$E_{50}$随着$E^{*}$的增大而增大,$k_{\text{r}}$对$E_{50}$的影响与$E^{*}$的大小有关,但整体而言$E_{50}$基本随着$k_{\text{r}}$的增大而减小；$E^{*}$变化范围内$E_{50}$的变化幅度大于$k_{\text{r}}$变化范围内$E_{50}$的变化幅度,说明$E_{50}$对$E^{*}$的敏感程度大于对$k_{\text{r}}$的敏感程度,这与前文敏感性分析结论一致。${\epsilon }_{\text{vmax}}$随着$E^{*}$的增大而减小,$k_{\text{r}}$对${\epsilon }_{\text{vmax}}$的影响同样与$E^{*}$的大小有关；对比$E^{*}$和$k_{\text{r}}$各自变化范围内${\epsilon }_{\text{vmax}}$的变化幅度可看出,${\epsilon }_{\text{vmax}}$对$E^{*}$的敏感程度大于对$k_{\text{r}}$的敏感程度,与前文敏感性分析结果相同。

图  6  割线模量$E_{50}$（黑色线,100 MPa）和最大压缩量${\epsilon }_{\text{vmax}}$（灰色线,%）等值线图

Figure  6.  Contours of secant modulus E₅₀ (black lines, 100 MPa) and maximum compression${\epsilon }_{\text{vmax}}$(gray lines, %)

下载: 全尺寸图片 幻灯片

$E_{50}$和${\epsilon }_{\text{vmax}}$相互借助来标定细观参数$E^{*}$和$k_{\text{r}}$。从图6中可见$E_{50}$和${\epsilon }_{\text{vmax}}$的等值线之间相互交叉,其中虚线分别为模拟值取室内试验值（$E_{50}=162.5\text{ MPa}$和${\epsilon }_{\text{vmax}}=0.912\%$）的等值线,两条虚线相交于点A,则点A对应的横纵坐标（$E^{*}=657.7\text{ MPa}$、$k_{\text{r}}=1.376$）即为同时满足$E_{50}$和${\epsilon }_{\text{vmax}}$的细观参数值,列于表7。

表  7  细观参数标定结果

Table  7.  The results of micro-parameters after calibration

细观参数	$E^{*}$/MPa	$k_{\text{r}}$	${\mu }_{\text{r}}$	$\mu$
标定值	657.7	1.376	0.6717	0.4177

下载: 导出CSV 

| 显示表格

4.3   确定滑动摩擦系数$\mu$和滚动阻力系数${\mu }_{\text{r}}$

固定已经确定的$E^{*}=657.7\text{ MPa}$和$k_{\text{r}}=1.376$不变,在优化后的$\mu$和${\mu }_{\text{r}}$细观参数范围内,令$\mu$分别取0.3,0.35,0.40,0.45,0.5,${\mu }_{\text{r}}$分别取0.3,0.4,0.5,0.6,0.7,重新开展了25次数值试验获取$S_{\text{P}}$和$\psi$的值并绘制其等值线图,如图7所示（黑色为$S_{\text{P}}$等值线,灰色为$\psi$等值线）。从图7中可见,$S_{\text{P}}$随着${\mu }_{\text{r}}$的增大而增大,并且$\mu$越大,$S_{\text{P}}$随${\mu }_{\text{r}}$的增大其增幅越大；同样地,$S_{\text{P}}$随着$\mu$的增大而增大,并且${\mu }_{\text{r}}$越大,$S_{\text{P}}$随$\mu$的增大其增幅越大；$\mu$变化范围内$S_{\text{P}}$的变化幅度大于${\mu }_{\text{r}}$变化范围内$S_{\text{P}}$的变化幅度,说明$S_{\text{P}}$对$\mu$的敏感程度大于对${\mu }_{\text{r}}$的敏感程度,与前文敏感性分析结论一致。$\psi$受${\mu }_{\text{r}}$的影响较小,而随着$\mu$的增大,$\psi$有较大增幅,同样可见$\psi$对$\mu$的敏感程度大于对${\mu }_{\text{r}}$的敏感程度。

图  7  峰值强度$S_{\text{P}}$（黑色线,MPa）和剪胀角ψ（灰色线,（°））等值线图

Figure  7.  Contours of peak strength S_p (black lines, MPa) and dilation angle ψ (gray lines, in degree)

下载: 全尺寸图片 幻灯片

$S_{\text{P}}$和$\psi$相互借助来标定细观参数$\mu$和${\mu }_{\text{r}}$。由于数值模拟未考虑颗粒的破碎,导致数值模拟获取的剪胀角均大于室内试验值,无法绘出剪胀角的室内试验值等值线,此处选取黑色虚线与剪胀角等值线交点中对应剪胀角最小的交点作为满足要求的点。等值线$S_{\text{P}}=3.1926\text{ MPa}$和$\psi =$12.53°的交点B对应的横纵坐标（${\mu }_{\text{r}}=0.6717$,$\mu =0.4177$）即为确定的细观参数,列于表7。当然增加剪胀角等值线的绘制数量和精度会影响点B的位置,但相差不大。

5.   细观参数验证

采用表7中确定的细观参数开展模拟并将结果与与室内三轴试验结果进行对比,如图8所示。从图中可见,偏应力-轴变关系模拟结果与室内试验结果吻合很好,模拟结果较好的预测了$E_{50}$和$S_{\text{P}}$,在轴向应变小于2%时,偏应力-轴变曲线基本重合,并且$S_{\text{P}}$的位置均大体位于轴向应变5%附近,说明文中标定的细观参数能够较好反映堆石料的偏应力-轴变关系。

图  8  室内试验与数值模拟曲线（围压800 kPa）

Figure  8.  Stress-strain curves of experimental and numerical tests under confining pressure of 800 kPa

下载: 全尺寸图片 幻灯片

对比图8中模拟结果和室内试验结果的体变-轴变关系曲线发现,${\epsilon }_{\text{vmax}}$基本相同,但$\psi$相差较大,其中室内试验表现为明显的体缩,这主要是因为室内试验堆石料在一定荷载下会发生颗粒破碎,而本次数值模拟未考虑颗粒的破碎,导致剪胀角数值模拟值较大,进而采用本文方法确定细观参数时是在满足峰值强度的前提下,基于剪胀角模拟结果最小值确定了细观参数,所以采用表7的细观参数获取得剪胀角模拟值与试验值相差较大。正因为数值模拟不考虑颗粒破碎时的体变-轴变曲线与室内试验结果存在较大差异,所以目前大部分学者标定堆石料细观参数时主要关注偏应力-轴变曲线结果的吻合^{[4-9, 12, 14-16]},而本文方法通过较少的计算量便可以达到相同的效果。

为了进一步验证标定结果,同样采用表7中的细观参数模拟了围压分别为400 kPa和1200 kPa下的三轴排水剪切试验。数值模拟与室内试验对比结果如图9所示。从图9中可以看出,无论400 kPa还是1200 kPa围压下,偏应力-轴变关系的模拟结果与室内试验结果基本吻合,在轴向应变小于2%时,偏应力-轴变曲线接近基本重合,这进一步说明文中标定的细观参数能较好的反映室内试验堆石料的偏应力-轴变关系。从体变-轴变关系曲线看,数值模拟结果与室内试验结果存在一定差异,其中,最大压缩量${\epsilon }_{\text{vmax}}$在围压400 kPa时相差不大,但围压取1200 kPa时,室内试验值大于数值模拟值；另外1200 kPa围压下数值模拟和室内试验最终体变差距大于围压400 kPa下的差距,这主要是因为室内试验围压越大颗粒破碎量越大,颗粒破碎对变形的影响占比也增大,而本次数值模拟没有考虑颗粒破碎,所以导致数值模拟的误差随着颗粒破碎程度的增大而增大。

图  9  室内试验与数值模拟曲线（围压400 kPa,1200 kPa）

Figure  9.  Stress-strain curves of experimental and numerical tests under confining pressures of 400 and 1200 kPa

下载: 全尺寸图片 幻灯片

6.   结论与建议

采用“正交-等值线”法快速标定了堆石料离散元模型细观参数,相比较于目前堆石料参数标定方法,该方法采用相对较少的计算量就能达到同样的效果,主要结论和建议如下：

（1）借助正交试验“试验次数少、代表性较好”的优势,通过较少的试验次数可以缩小细观参数取值范围,进而可以确定细观参数敏感性分析时细观参数的取值范围并减轻后续参数标定的工作量。如果对需要标定材料的细观参数取值范围没有经验可依,知之较少时,该方法的优越性将更加明显。

（2）细观参数的敏感性分析可知,在优化后细观参数范围内,各细观参数对不同宏观指标的影响程度不完全相同。$E_{50}$：$E^{*}>\mu >k_{\text{r}}>{\mu }_{\text{r}}$；$S_{\text{P}}$：$\mu >{\mu }_{\text{r}}>$ $E^{*}>k_{\text{r}}$；${\epsilon }_{\text{vmax}}$：$E^{*}>\mu >k_{\text{r}}>{\mu }_{\text{r}}$；$\psi$：$\mu >E^{*}>{\mu }_{\text{r}}>k_{\text{r}}$。其中$E^{*}$对$E_{50}$和${\epsilon }_{\text{vmax}}$的影响绝对占优,$k_{\text{r}}$,${\mu }_{\text{r}}$,$\mu$影响可以忽略。$E^{*}$,$k_{\text{r}}$,${\mu }_{\text{r}}$,$\mu$这4个细观参数对$S_{\text{P}}$和$\psi$均有一定影响。

（3）基于细观参数敏感性的分析,抓住主要影响因素,借助等值线法基于离散的数据点通过插值获取其范围内任意值的优势,提出了细观参数的标定方法。即首先根据试样的割线模量$E_{50}$和最大压缩量${\epsilon }_{\text{vmax}}$确定有效模量$E^{*}$和法向切向刚度比$k_{\text{r}}$,然后根据试样的峰值强度$S_{\text{P}}$和剪胀角$\psi$确定滚动阻力系数${\mu }_{\text{r}}$和滑动摩擦系数$\mu$。

（4）标定后的细观参数值为：$E^{*}=657.7\text{ MPa}$,$k_{\text{r}}=1.376$,${\mu }_{\text{r}}$=0.6717,$\mu =0.4177$。DEM模型采用该细观参数能够很好的再现不同围压下室内堆石料的偏应力-轴变关系。

本文数值模拟过程中由于没有考虑颗粒的破碎,导致最终标定的细观参数未能很好的再现堆石料体变-轴变曲线,但本文的重点是离散元模型细观参数标定方法的提出,后续将考虑颗粒破碎开展研究。另外,针对滚动阻力模型以外的其它DEM模型,将进一步研究本文方法的普适性。