0.   引言

颗粒流模型（PFC）的细观力学参数是通过“试错法”进行间接获取的，但该方法常常耗费科研人员大量的时间和精力。国内外学者对岩石颗粒流模型的宏细观参数的相关性及取值方法究展开了一系列研究。专家们对宏-细观参数关系的研究深度从定性关系研究（Huang^[1]）到定量关系研究（Yang等^[2]，丛宇等^[3]），再到少部分的理论关系（刘新荣等^[4]、陈亚东等^[5]）研究。目前较多的研究为细观参数的优化及算法的改进研究：一开始Yoon^[6]利用中心复合设计法及优化参数来获得合适的细观参数，随后学者们^[7-8]陆续的提出了众多改进的算法（改进的粒子群优化（PSO）标定方法等）来获得细观参数。虽然这些方法最终都获得了所需的细观参数，但并不能从机理上解释细观参数是如何作用的，同时理论方面的研究也较少，因此一套系统的关于宏-细观参数关系的理论研究是有重要意义的。

本文以接触黏结模型（简称CBM）为研究对象，首先筛选出有代表性的四颗粒-二粒径单元体，随后理论推导获得代表性单元体的宏细观参数解析解。再通过试验设计（简称DOE）定性研究岩石材料宏观力学参数的响应。随后基于大量数值试验，获得宏细观参数基于解析方程的拟合关系式。最终基于宏细观参数的内在联系建立半解析查表法。

1.   颗粒流材料宏-细观参数的解析关系

1.1   代表性均质模型的筛选

筛选一个能够用于构成均质颗粒流模型的代表性单元体是研究颗粒流黏结材料宏观-细观参数理论关系的前提。石崇等^[9]对等粒径的颗粒流模型进行了理论研究，并用几何特征角θ（如图 3所示）及配位数N来描述等径颗粒流模型的细观几何特征，其中几何特征角θ为产生接触的两个颗粒的圆心连线与水平线的夹角，配位数N表示模型内每个颗粒上的平均接触个数。本文选取颗粒数目为四的单元体进行分析，并设置两种粒径以研究粒径比对单元体力学特性的影响。该单元体的几何特征可通过粒径比R₂/R₁与几何特征角θ进行定量描述。随着R₂/R₁与θ的变化，颗粒的相对位置及接触黏结点的数目亦发生变化。对于用于构成均质颗粒流模型的基本单元体，其粒径比与几何特征角需满足下述条件：

（1）颗粒边界不重叠。

（2）均质模型中，小颗粒之间不应发生大量接触。

如图 1所示，根据条件（1），（2）可绘制出a，b，c，d，e，f 6条曲线，图 1中区域Ⅰ、Ⅱ、Ⅲ满足条件（1），（2）。标记为1~10的单元体为可用于构成均质模型的单元体。本文以单元体5为例进行分析研究。如图 2所示，由单元体5构成的颗粒流模型为均质模型，该模型的力学特性可通过分析模型的基本单元即单元体5获得，故将单元体5称为代表性单元体。假设均质模型水平排列a个基本单元，竖向排列b个基本单元，则均质模型的细观组构特征如表 1所示。由表 1可知，均质模型的配位数N≈5，与本文的数值试验模型（N=3.5~5.5）相符。

图  1  二粒径-四颗粒单元体的几何形态

Figure  1.  Geometric morphology of units with four granules

下载: 全尺寸图片 幻灯片

图  2  代表性单元体

Figure  2.  Representative unit

下载: 全尺寸图片 幻灯片

表  1  均质模型的几何参数

Table  1.  Geometric parameters of homogeneous model

Rmax/Rmin	θ	颗粒配位数N
		小颗粒	大颗粒	平均值
${\text{(1, }}\sqrt 2 + 1)$	(30°, 45°)	4	6	$5 - {{\left( {5a + 5b} \right)} \mathord{\left/ {\vphantom {{\left( {5a + 5b} \right)} {\left( {2ab + a + b} \right)}}} \right. } {\left( {2ab + a + b} \right)}}$

下载: 导出CSV 

| 显示表格

1.2   均质模型的细观破坏模式

（1）直接拉伸状态下

以直接拉伸状态为例详细推导基本单元的内力分布规律及变形特性。在沿y轴方向的拉伸条件下均质模型及基本单元的受力情况如图 3所示。

图  3  拉伸条件下均质模型及其基本单元的受力情况

Figure  3.  Force analysis of homogeneous model and its basic unit under tensile conditions

下载: 全尺寸图片 幻灯片

平衡方程：

如图 4（a）所示，假设在基本单元内，接触点C_i（i=1~5）中存在的内力为：法向力f_ni（如绿色箭头所示，以拉为正），切向力f_si（如红色箭头所示，以顺时针旋转为正）。分别对O_A，O_B，O_C进行受力分析，得到6组互相独立的力平衡方程，如式（1）所示。

图  4  拉伸条件下基本单元的内力分析及变形特性

Figure  4.  Analysis of internal forces and deformation characteristics of basic unit under tensile conditions

下载: 全尺寸图片 幻灯片

变形协调方程：

如图 4（b）所示，接触点C_i（i=1~5）的法、切向位移为${\varDelta _{ni}}$，${\varDelta _{si}}$ (i=1~5)。令O_D的圆心位移为零，则O_A，O_B，O_C圆心位置产生的x轴向、y轴向位移分别为${\boldsymbol{\varDelta} _i}$ (i=1~6)。本文中，${\boldsymbol{\varDelta} _{{\text{n}}i}}$，${\boldsymbol{\varDelta} _{{\text{s}}i}}$及${\boldsymbol{\varDelta} _i}$为矢量，${\boldsymbol{\varDelta} _{{\text{n}}i}}$与${\boldsymbol{\varDelta} _{{\text{s}}i}}$的正方向与法、切向接触力的一致，${\boldsymbol{\varDelta} _i}$的正方向与x-y坐标系的一致。颗粒圆心位移${\varDelta _{}}_1$~${\varDelta _{}}_6$与接触点法、切向位移${\varDelta _{{\text{n}}i}}$，${\varDelta _{{\text{s}}i}}$(i=1~5)的关系如下所示。

物理方程：

CBM接触点处颗粒-颗粒的法向力-法向应力、切向力-切向应力关系见公式（3）。其中A=2Rt，2D模型中t=1，R=min(R₁, R₂)，R₁，R₂分别为基本单元中O_A，O_B的半径。

颗粒-颗粒的法向力-法向位移、切向力-切向位移关系见公式（4）。其中Δ_n，Δ_s为接触点上的法向，切向变形量，k_n，k_s为接触点上的法向，切向刚度。

颗粒-颗粒的刚度与弹性模量的关系如公式（5）所示。其中，L= R₁+R₂= R_min+ R_max。

联立公式（1），（2），（4）可得，外应力为${\sigma _{\text{y}}}$=F_yi/ (2cosθ(R_min+R_max))时，各接触点的接触力及位移如公式（6）所示，基本单元的内力分布规律及变形特性，如图 5（a）所示。

图  5  拉伸条件下基本单元的内力分布规律及变形特性

Figure  5.  Distribution of internal forces and deformation characteristics of unit under tensile conditions

下载: 全尺寸图片 幻灯片

其中，$m = 1 + {\sin ^2}\theta + {\cos ^2}\theta {k_{\text{n}}}/{k_s}$，$n = (1 + {\sin ^2}\theta )/$ ${k_{\text{n}}}/{k_s} + {\cos ^2}\theta$。

同理可得，在沿x方向的直接拉伸作用下，外应力${\sigma _x}$=${F_{xi}}$/(2sinθ(R_min+R_max))，各接触点的接触力及位移如公式（7）所示，基本单元的内力分布规律及变形特性如图 5（b）所示。

根据式（6），（7）可得模型的宏观抗拉强度${\sigma _{{\text{t}}x}}$或${\sigma _{{\text{t}}y}}$与接触法向黏结强度${\sigma _{{\text{cn}}}}$、接触切向黏结强度${\tau _{{\text{cs}}}}$的关系，详见表 2。由表 2可知，接触强度比${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$为决定破坏模式的重要因素。在下文中式子1/(1+R_max/ R_min)均简写为RR，刚度比k_n/k_s均简写为kk。

表  2  颗粒流材料的破坏模式

Table  2.  Failure modes of PFC materials

序号	判据	τcs/σcn < 判据			τcs/σcn > 判据
		破坏模式	宏观强度与细观参数关系式		破坏模式	宏观强度与细观参数关系式
1	$\frac{{{\text{2}}\tan \theta }}{{kk + 1}}$	拉致剪裂	${\sigma _{{\text{t}}x}} = {\tau _{{\text{cs}}}}{\text{RR}}\frac{m}{{{{\sin }^2}\theta }}$		拉致拉裂	${\sigma _{{\text{t}}x}} = {\sigma _{{\text{cn}}}}{\text{RR}}\left( {2\cot \theta + \frac{{4\tan \theta }}{{kk + 1}}} \right)$
2	$\frac{{{R_{{\text{max}}}}}}{{2{R_{{\text{min}}}}}}\frac{{\cos \theta }}{{m - 1}}$	拉致剪裂	${\sigma _{{\text{t}}y}} = {\tau _{{\text{cs}}}}{\text{RR}}\frac{{2m}}{{{{\cos }^2}\theta }}$		拉致拉裂	${\sigma _{{\text{t}}y}} = {\sigma _{{\text{cn}}}}(1 - {\text{RR}})\left( {\frac{1}{{\cos \theta }} + \frac{1}{{{{\cos }^3}\theta (kk + {{\tan }^2}\theta )}}} \right)$
3	$\frac{{{R_{{\text{max}}}}}}{{{R_{{\text{min}}}}}}\frac{1}{{\cos \theta (kk - 1)}}$	压致剪裂	${f_{{\text{c}}x}} = {\tau _{{\text{cs}}}}{\text{RR}}(2 + {\cos ^2}\theta (kk + 1))$		压致拉裂	${f_{{\text{c}}x}} = {\sigma _{{\text{cn}}}}(1 - {\text{RR}})\left( {\frac{{\cos \theta }}{{{{\sin }^2}\theta }} + \frac{2}{{\cos \theta {{\sin }^2}\theta (kk - 1)}}} \right)$
4	—	压致剪裂	${f_{{\text{c}}y}} = {\tau _{{\text{cs}}}}{\text{RR}}(4{\tan ^2}\theta + 2(kk + 1))$		—	—
5	$\frac{{2{R_{{\text{min}}}}}}{{{R_{{\text{max}}}}}}\frac{{{\text{(}}1 + {{\sin }^2}\theta {\text{)}}kk}}{{\cos \theta }}$	剪致剪裂	${\tau _{\max }} = {\tau _{{\text{cs}}}}(1 - {\text{RR}})\frac{n}{{\cos \theta }}$		剪致拉裂	${\tau _{\max }} = 2{\sigma _{{\text{cn}}}}{\text{RR}}\frac{n}{{n - {{\cos }^2}\theta }}$

下载: 导出CSV 

| 显示表格

（2）单轴压缩状态下

式（6），（7）。沿y轴方向、x轴方向的压应力为${\sigma _y}$=F_yi/(2(R_min+R_max)cosθ)，${\sigma _x}$=F_xi/(2(R_min+R_max)sinθ)，基本单元的内力分布规律及变形特性如图 6所示。模型的单轴抗压强度f_cx或f_cy与${\tau _{{\text{cs}}}}$，${\sigma _{{\text{cn}}}}$的关系详见表 2。

图  6  单轴压缩条件下基本单元的内力分布规律及变形特性

Figure  6.  Distribution of internal forces and deformation characteristics of unit under uniaxial compression conditions

下载: 全尺寸图片 幻灯片

单轴压缩状态下，各接触点的接触力及位移同公

（3）简单剪切状态下

简单剪切状态下，剪应力τ_xy=τ_yx =F_yxi/(2(R_max+ R_min)cosθ)，基本单元的内力分布规律及变形特性如图 7所示，各接触点的接触力及位移如公式（8）所示。根据公式（8），纯剪条件下模型的最大剪应力${\tau _{{\text{max}}}}$与接触切向黏结强度${\tau _{{\text{cs}}}}$、接触法向黏结强度${\sigma _{{\text{cn}}}}$的关系详见表 2。

图  7  简单剪切条件下基本单元的内力分布规律及变形特性

Figure  7.  Distribution of internal forces and deformation characteristics of unit under simple shear conditions

下载: 全尺寸图片 幻灯片

1.3   均质模型宏-细观参数的解析关系

（1）抗拉强度

在方案A、B中，${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$的范围为1.67~50，k_n/k_s > 1，计算可得：${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$大于判据1，2（详见表 2），既宏观抗拉强度仅与法向接触强度有关。根据表 2，宏观抗拉强度${\sigma _{\text{t}}}$与细观参数的关系可简化为式（9），（10）。式中，a_D，b_D，c_D，d_D为与几何特征角有关的系数。

学者们往往是通过巴西劈裂试验来获得类岩石材料的抗拉强度，结合课题组的已有研究^[10]，巴西抗拉强度${\sigma _{\text{t}}}$的表达式如下。K_B为统计获得的拟合系数。

（2）单轴抗压强度

根据表 2，单轴抗压强度f_c与细观参数的关系如式（12），（13）所示。其中，a_U，b_U，c_U，d_U为与几何特征角有关的系数。

（3）弹性模量

由式（3），（5）~（7）可得，宏观弹性模量E_x，E_y与细观参数的关系：

根据上述关系式，弹性模量的表达式可归纳为

式中：e_U，f_U，g_U为与几何特征角有关的系数。

（4）泊松比

由公式（3），（5）~（7）可得，泊松比${\nu _{xy}}$，${\nu _{yx}}$与细观参数的关系：

根据上述关系式，泊松比的表达式可归纳为

式中：h_U为与几何特征角有关的系数。

（5）剪切强度

由上述推导可得基本单元在具有y轴向压应力及剪切应力状态下各接触点的应力分布。随着剪应力的增加，接触黏结点C₁，C₃先发生拉伸破坏，剪切强度${\tau _{\text{f}}}$与接触法向黏结强度${\sigma _{{\text{cn}}}}$，法向应力$\sigma$的关系如下：

式中：$a = 2\sin \theta mn/(mkk(1 + {\sin ^2}\theta ) + {\sin ^2}\theta )$；$b =$ $\sin \theta \cos \theta n/(mkk(1 + {\sin ^2}\theta ) + {\sin ^2}\theta )$。

在接触破坏后，细观摩擦系数发挥作用，此时，剪切强度${\tau _{\text{f}}}$与细观参数的关系如下：

式中：$a=\mathrm{sin}\theta \mathrm{cos}\theta n/(mn-{\mathrm{cos}}^{2}\theta m+{\mathrm{sin}}^{2}\theta n)\text{ }$；$b = 2\mathrm{sin}\theta \mathrm{cos}\theta (1+kk)/(mn-{\mathrm{cos}}^{2}\theta m+{\mathrm{sin}}^{2}\theta n)-\mathrm{cot}\theta ；c= 2\mathrm{sin}\theta \mathrm{cos}\theta (1+kk)/(mn-{\mathrm{cos}}^{2}\theta m+{\mathrm{sin}}^{2}\theta n)。$

根据上述推导，可得基本单元在具有y轴向压应力${\sigma _3}$（最小主应力）及x轴向压应力${\sigma _1}$（最大主应力）下各接触点的应力分布。随着${\sigma _1}$的增加，接触点C₁~C₄发生剪切破坏，应力${\sigma _1}$，${\sigma _3}$与${\tau _{{\text{cs}}}}$的关系如下：

基于莫尔库仑准则可知，

根据式（22），（23），可得黏聚力c：

在y轴向压应力为${\sigma _1}$及x轴向压应力为${\sigma _3}$时，黏聚力c与细观参数的关系同式（24）。

由上述分析可知，黏聚力与细观接触强度呈线性相关，根据式（20），（24），可知黏聚力c与细观参数的关系，如下式所示：

值得注意的是，式（20），（22）中未考虑接触摩摩擦系数μ对宏观摩擦系数tan(φ)的影响。将θ=30°~45°，kk=3.5代入式（25），可得a_n=1.15~1.51，a_s=4.6~6.33。

根据式（21），可得宏观摩擦系数tan(φ)与细观参数的关系，如式（26）所示。代入θ=30°~45°，kk=3.5，得b_S=0.122~0.144，c_S=-1.518~-0.77，d_S= 0.214~0.227。

2.   数值试验及结果

2.1   基于Plackrtt-Burman设计法的敏感性分析

采用Plackrtt-Burman设计方法，基于巴西试验、单轴压缩试验及直剪试验，对细观参数进行敏感性分析。细观参数的较小值、较大值参照以往学者^{[2, 6]}的研究选取。随后对宏观响应进行因子分析，置信水平取95%（a=0.05），作出标准化效应直方图，如图 8所示，图中跨越参考线的条形项在统计意义上显著。

图  8  标准化效应的Pareto图

Figure  8.  Pareto charts of standardization effects

下载: 全尺寸图片 幻灯片

2.2   数值模拟方案

依据细观参数的敏感性分析结果，设计9组方案进行数值模拟。试验变量的取值范围参照以往学者^{[2, 6]}的研究进行选取，试验常量的取值参考傅晏^[11]模拟微风化砂岩的细观参数进行选取（见表 3）。

表  3  细观参数取值表

Table  3.  Values of microscopic parameters

细观参数		试验常量取值
几何参数	最小粒径Rmin/mm	0.2
	粒径比Rmax/Rmin	1.67
细观强度参数	接触法向黏结强度σcn/MPa	20
	接触切向黏结强度τcs/MPa	20/50
	颗粒摩擦系数μ	1.7
细观刚度参数	刚度比kn/ks	3.5
	颗粒接触模量Ec/GPa	10

下载: 导出CSV 

| 显示表格

方案A：对56组（不同粒径）模型进行5次直接拉伸试验，${\sigma _{{\text{cn}}}}$分别为1，5，15，20及30 MPa，k_n/k_s均为3.5，共280次。方案B：对21组（不同粒径）模型进行12次直接拉伸试验，${\sigma _{{\text{cn}}}}$分别为1，5，15 MPa，k_n/k_s分别为1，2，3.5及6，共252次；方案C：对模型（R_min=0.2，R_max/R_min=1.67）进行单轴压缩试验，${\sigma _{{\text{cn}}}}$，${\tau _{{\text{cs}}}}$分别为5，10，20，30，50，70 MPa，共30次。方案D：对36组（不同粒径）模型进行单轴压缩试验，k_n/k_s为1.5~11，强度比${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$为0.5，1及2，共392次。方案E：对3组（不同粒径）模型进行单轴压缩试验，μ为0.1，0.9，1.7，2.5及5，共18次。方案F：对8组（不同粒径）模型进行单轴压缩试验，E_c为1，5，10，15及20，共40次。方案G：对模型（R_min=0.2，R_max/R_min=1.67）进行直剪试验，${\sigma _{{\text{cn}}}}$，${\tau _{{\text{cs}}}}$分别为5，10，20，30，50，70 MPa，共90次。方案H：对的模型（R_min=0.2，R_max/R_min=1.67）进行直剪试验，μ为0.1，1，1.7，3，5，7及10，共24次。方案Ⅰ：对模型（R_min=0.2，R_max/R_min=1.67）进行直剪试验，k_n/k_s为1.5，2，2.5，3，3.5，4，5，6，9，共24次。

本文直接拉伸、单轴压缩及直接剪切试验的数值模型的尺寸分别为50 mm×100 mm、50 mm×50 mm及50 mm×50 mm，部分模型如图 9所示。

图  9  数值试件模型

Figure  9.  Numerical models

下载: 全尺寸图片 幻灯片

2.3   结果分析

（1）抗拉强度

基于式（9），（11）对抗拉强度（巴西劈裂^[10]、方案A的结果）与${\sigma _{{\text{cn}}}}$进行线性回归分析，获得拟合系数K_B，K_D，详见表 4。由表 4可知，细观参数一定时，K_D约为K_B的两倍；K_D，K_B随着R_max/R_min的增长而增长。K_B的相关系数R²主要在0.99~1.0，K_D的相关系数R²均大于0.999，由此可知宏观抗拉强度与法向黏结强度呈显著的线性相关。

表  4  拟合系数K_B与K_D

Table  4.  Fitting coefficients K_B and K_D

系数	Rmin/ mm	Rmax/Rmin
		1	2	3	4	5	6	7	8
KB	0.1	0.609	0.768	0.997	1.090	1.086	1.368	2.015	2.568
	0.15	0.589	0.861	1.137	1.271	1.303	1.317	1.616	2.044
	0.2	0.602	0.817	0.987	1.270	1.583	1.949	2.559	1.983
	0.25	0.737	0.970	0.953	1.213	1.734	2.099	1.926	1.870
	0.3	0.677	0.779	1.120	1.386	1.612	1.692	2.198	2.054
	0.35	0.633	0.839	1.183	1.348	1.988	1.478	1.877	2.008
	0.4	0.541	1.064	1.155	1.247	1.948	1.795	2.155	2.698
系数	Rmin/	Rmax/Rmin
	mm	1.2	1.67	2	3	4	5	6	7
KD	0.1	1.320	1.414	1.677	2.175	2.216	2.960	3.443	3.894
	0.15	1.234	1.444	1.614	2.009	2.099	2.698	3.260	3.751
	0.2	1.308	1.425	1.622	1.971	2.450	2.630	3.204	—
	0.3	1.479	1.479	1.605	2.067	2.330	—	—	—
	0.4	1.186	1.425	1.527	1.844	—	—	—	—
	0.5	1.267	1.475	1.542	—	—	—	—	—

下载: 导出CSV 

| 显示表格

将方案B的结果代入式（9），（10），可得系数a_D，b_D，c_D，d_D的拟合值，详见表 5。将θ=30°~45°代入上述公式，可得各系数的理论值范围，详见表 6。由表 5，6可知，a_D，b_D的拟合值较理论值偏小，c_D拟合值略小于理论值，d_D拟合值与理论值吻合。由此可知，对于抗拉强度而言，式（10）更具代表性。

表  5  方案B的拟合结果

Table  5.  Fitting results of Plan B

Rmin/ mm	拟合结果	Rmax/Rmin				拟合结果	Rmax/Rmin
		1.67	2	3	4		1.67	2	3	4
0.1	aD	1.3	1.47	1.83	—	cD	0.78	0.73	0.62	—
	bD	0.65	0.85	1.07	—	dD	1.54	1.01	0.62	—
	R2	0.961	0.86	0.85	—	R2	0.92	0.93	0.72	—
0.15	aD	1.32	1.39	1.78	1.83	cD	—	—	0.61	—
	bD	0.48	0.62	0.97	1.26	dD	—	—	0.57	—
	R2	0.89	0.74	0.98	0.99	R2	—	—	0.94	—
0.2	aD	1.27	1.41	—	2.16	cD	0.76	0.70	—	—
	bD	0.67	0.76	—	1.19	dD	1.14	0.67	—	—
	R2	0.99	0.92	—	1	R2	0.99	1.00	—	—
0.3	aD	—	1.41	1.77	2.1	cD	—	0.71	0.61	0.55
	bD	—	0.89	1.03	1.17	dD	—	0.77	0.55	0.26
	R2	—	0.92	0.9	0.99	R2	—	0.97	0.81	1.00
0.4	aD	1.26	1.36	1.61	—	cD	0.76	0.68	—	—
	bD	0.71	0.65	0.99	—	dD	1.05	1.06	—	—
	R2	0.99	0.92	0.99	—	R2	0.99	0.88	—	—
0.5	aD	1.34	1.32	—	—	cD	0.77	0.70	0.57	—
	bD	0.66	0.92	—	—	dD	0.92	0.89	0.42	—
	R2	0.94	0.96	—	—	R2	0.87	0.96	0.99	—

下载: 导出CSV 

| 显示表格

表  6  拟合系数的理论值与模拟值

Table  6.  Theoretical and simulated ranges of fitting coefficients

试验类型	系数	理论值范围	拟合值范围
直接拉伸试验	aD	2~3.46	1.26~2.16
	bD	2.31~4	0.48~1.26
	cD	1.15~1.41	0.55~0.78
	dD	0.33~1	0.26~1.54
单轴压缩试验	aU	1.41~3.46	1.11~1.93
	bU	5.67~9.24	1.01~2.53
	cU	1.33~4	1.74~7.17
	dU	1~3	0.08~0.77
	eU	0.33~1.04	0.25~0.86
	fU	0.77~2	2.02~4.5
	gU	0.33~1	1.95~3.45
	hU	1~3	1.698~2.21

下载: 导出CSV 

| 显示表格

（2）单轴抗压强度

接触强度对单轴抗压强度f_c的影响如图 10所示。由图 10可知，f_c与接触强度较小值min(${\sigma _{{\text{cn}}}}$，${\tau _{{\text{cs}}}}$)呈较好的线性相关。该结论与表 2中判据3的结果相符。

图  10  单轴抗压强度f_c与细观接触强度的关系（方案C）

Figure  10.  Relationship between f_c and contact bond strength (Plan C)

下载: 全尺寸图片 幻灯片

以${\sigma _{{\text{cn}}}}$为x，以${\tau _{{\text{cs}}}}$为y，以f_c为z，基于式z=ax+by，对方案C结果进行拟合，拟合结果见表 7。由该表可知，${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$ < 1时，f_c主要受${\tau _{{\text{cs}}}}$影响，可认为在${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$ < 0.5时，将不再受${\sigma _{{\text{cn}}}}$的影响；${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$ > 1时，f_c主要受${\sigma _{{\text{cn}}}}$影响，可认为${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$ > 2后，将不再受${\tau _{{\text{cs}}}}$的影响。

表  7  单轴抗压强度的拟合结果（方案C）

Table  7.  Fitting results of uniaxial compressive strength (Plan C)

τcs/σcn	拟合结果
	a	b	R2
0.25~0.428	0.0398	1.7570	0.999
0.5~0.714	0.1392	1.4726	0.999
1	a+b=1.388		0.997
1.4~2	1.591	-0.0740	0.998
2.33~5	1.549	-0.0318	0.999
6~14	1.419	-0.0010	1

下载: 导出CSV 

| 显示表格

k_n/k_s对f_c的影响如图 11所示。由图 11可知，随着k_n/k_s的增大，f_c均为先增大后减小，但各强度比转折点的位置不同，且强度比越大，转折点出现越早。该结论与理论结果相符。由表 2可知，随着k_n/k_s的增加（${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$不变），判据3减小，试件从压致剪裂转变为压致拉裂。压致剪裂的试件与k_n/k_s呈线性正相关（见式（13）），压致拉裂的试件与k_n/k_s呈负相关（见式（12））。

图  11  单轴抗压强度vs刚度比k_n/k_s（方案D）

Figure  11.  Uniaxial compressive strengths vs. stiffness ratios (Plan D)

下载: 全尺寸图片 幻灯片

根据式（12），（13）对部分模拟结果进行拟合，获得系数a_U，b_U，c_U，d_U，如表 8所示，系数理论值与模拟值的范围见表 6。由表 6可知，a_U，c_U的模拟值与理论值基本吻合，b_U，d_U的模拟值偏小于理论值。这说明，本文选取的模型，在用于分析k_n/k_s对f_c的影响时，其代表性不足。

表  8  方案D单轴抗压强度拟合结果

Table  8.  Fitting results of uniaxial compressive strength (Plan D)

τcs/σcn	Rmax/Rmin	Rmin/mm	式（13）
			cU	dU	R2
0.5	1.67	0.2	2.949	0.331	0.95
	2	0.15	3.924	0.165	0.88
	2	0.3	3.377	0.251	0.98
	3	0.1	4.333	0.331	0.99
	4	0.15	4.696	0.589	0.91
	5	0.15	7.170	0.267	0.87
1	1.2	0.1	2.435	0.081	0.99
	1.2	0.4	1.736	0.342	0.99
	1.67	0.3	2.358	0.270	0.98
	1.67	0.5	2.735	0.098	0.99
	2	0.1	3.907	0.283	0.99
	2	0.3	2.302	0.395	0.86
	2	0.4	2.689	0.244	0.98
	3	0.1	3.464	0.200	0.98
	3	0.4	2.941	0.358	0.98
	4	0.1	4.628	0.350	0.93
	4	0.2	5.166	0.147	0.98
	5	0.1	5.526	0.140	1.00
	5	0.15	4.638	0.617	0.91
	6	0.15	6.239	0.330	1.00
	6	0.2	4.615	0.767	0.99
	7	0.15	7.654	0.485	0.99
τcs/σcn	Rmax/Rmin	Rmin/mm	式（12）
			aU	bU	R2
2	1.67	0.2	1.932	1.007	0.99
	2	0.15	1.723	1.077	0.71
	2	0.3	1.461	2.190	0.72
	3	0.1	1.210	2.186	0.83
	3	0.15	1.482	1.047	0.99
	3	0.2	1.116	2.301	1.00
	3	0.3	1.133	2.529	0.99

下载: 导出CSV 

| 显示表格

摩擦系数μ及接触模量E_c对峰值应力f_c的影响如图 12所示。由图 12（a）可知，随着μ的增大，f_c呈现先增大后平缓的变化趋势，其拐点大致在μ=2处；由图 12（b）可知，随着E_c的增大，f_c呈现先减小后平缓的变化趋势，其拐点大致在E_c=10 GPa处。

图  12  单轴抗压强度与摩擦系数及接触模量的关系(方案E，F)

Figure  12.  Relationship among f_c, E_c and μ (Plan E, F)

下载: 全尺寸图片 幻灯片

（3）弹性模量

接触模量E_c对杨氏模量E及泊松比$\nu$的影响如图 13所示。由图 13（a）可知，E与E_c呈显著的线性相关；由图 13（b）可知，随着E_c的增大，$\nu$呈现先增大后平缓的变化趋势，其拐点大致在E_c=10 GPa处。

图  13  接触模量E_c与杨氏模量E、泊松比$\nu$的的关系（方案F）

Figure  13.  Relationship among E_c, E and ν (Plan F)

下载: 全尺寸图片 幻灯片

基于式（16），对方案F的结果进行线性拟合，得到拟合参数K_E，详见表 9。再基于式（16）对方案D的结果进行拟合，得到拟合系数e_U，f_U，g_U，详见表 10。由表 6可知，e_U的模拟值与理论值吻合的较好，f_U的模拟值比理论值偏大，而g_U的模拟值与理论值则具有较大差异。

表  9  拟合参数K_E及其拟合优度R²

Table  9.  Fitting coefficient K_E and its goodness of fit R²

颗粒几何参数			拟合参数及R2
Rmax/Rmin	Rmin/mm		KE	R2
1.67	0.1		1.5596	0.9945
	0.2		1.5052	0.9938
2	0.15		1.6553	0.9947
	0.3		1.632	0.9939
3	0.1		2.272	0.99
	0.3		2.0795	0.9921
4	0.15		2.5885	0.9938
	0.2		2.4767	0.9889

下载: 导出CSV 

| 显示表格

表  10  方案D中弹性模量的拟合系数及其拟合优度R²

Table  10.  Fitting coefficient of elastic modulus and its goodness of fit R²(Plan D)

颗粒几何参数			拟合参数及R2
Rmax/Rmin	Rmin/mm		eU	fU	gU	R2
1.2	0.15		0.339	2.201	2.616	0.999
	0.2		0.294	2.565	3.086	0.998
	0.3		0.295	2.579	3.174	0.997
	0.4		0.295	2.629	3.446	0.999
	0.5		0.314	2.026	2.507	0.981
	1		0.254	2.634	3.360	0.996
1.67	0.1		0.436	2.271	2.414	0.998
	0.15		0.352	2.815	3.063	0.996
	0.2		0.385	2.411	2.573	0.992
	0.3		0.329	2.931	3.214	0.992
2	0.1		0.460	2.567	2.374	0.998
	0.15		0.423	2.617	2.508	0.997
	0.2		0.422	2.453	2.274	0.998
	0.4		0.365	3.090	3.207	0.991
3	0.1		0.541	3.387	2.518	0.996
	0.15		0.581	2.770	2.184	0.973
	0.2		0.466	3.776	2.827	0.996
	0.3		0.557	2.865	2.176	0.981
4	0.1		0.672	3.854	2.186	0.996
	0.15		0.641	4.311	2.758	0.993
	0.2		0.701	3.387	1.955	0.993
	0.3		0.662	3.610	2.270	0.988
5	0.1		0.864	4.030	2.428	0.953
	0.15		0.788	4.547	2.394	0.989

下载: 导出CSV 

| 显示表格

（4）泊松比

基于式（19）对方案D的结果进行拟合，得到拟合系数h_U，结果见表 11，其相关系数R²的范围为0.86~0.99。由表 11可知，颗粒粒径对h_U的影响较小。同时由表 6可知，h_U的模拟值均在理论值范围内。

表  11  方案D中泊松比的拟合系数h_U及其R²（方案D）

Table  11.  Fitting coefficient h_U and its goodness of fit R²(Plan D)

Rmin/mm	强度比	Rmax/Rmin
		1.2	1.67	2	3	4	5	6
0.1	0.5	—	—	—	1.960	—	—	—
	1	2.127	2.210	2.017	1.916	1.906	1.980	2.050
	2	—	2.161	2.125	1.987	2.015	2.125	—
0.15	0.5	—	—	2.039	—	2.054	2.027	—
	1	1.982	1.914	1.960	1.94	1.962	1.981	1.933
	2	—	2.026	1.987	1.983	2.096	2.052	—
0.2	0.5	—	1.987	—	—	—	—	—
	1	1.933	1.896	1.871	2.019	1.766	1.804	2.104
	2	—	1.948	1.983	2.093	1.964	1.910	—
0.3	0.5	—	—	2.055	—	—	—	—
	1	1.901	1.851	1.916	1.858	1.883	—	—
	2	—	2.001	1.938	1.953	1.935	—	—
0.4	0.5	—	—	—	—	—	—	—
	1	1.822	1.941	1.904	1.698	—	—	—
	2	—	2.016	2.072	—	—	—	—
0.5	0.5	—	—	—	—	—	—	—
	1	1.788	1.851	1.874	—	—	—	—
	2	—	1.909	1.991	—	—	—	—

下载: 导出CSV 

| 显示表格

（5）黏聚力

接触强度对对黏聚力c的影响如图 14所示。由图 14可知，c与接触强度得较小值min(${\sigma _{{\text{cn}}}}$, ${\tau _{{\text{cs}}}}$)呈较好的线性相关。以${\sigma _{{\text{cn}}}}$为x，以${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$为y，以c为z，基于式z=x(a+by)进行拟合，结果见表 12。由表 12可知，${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$ < 0.7时，c的主要影响因素为${\tau _{{\text{cs}}}}$，${\sigma _{{\text{cn}}}}$的影响极小；${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$ > 1.4时，c的主要影响因素为${\sigma _{{\text{cn}}}}$，${\tau _{{\text{cs}}}}$的影响极小。最后，根据方案Ⅰ的结果可知，当k_s/k_n=2~6时，k_s/k_n对黏聚力c的影响较小。

图  14  黏聚力与细观强度的关系(方案G)

Figure  14.  Relationship between c and contact bond strecngth

下载: 全尺寸图片 幻灯片

表  12  系数a，b

Table  12.  Fitting coefficients a and b

切法向强度比	拟合公式：c =σcn×(a+b×τcs/σcn)
	a	b	R2
0.25~0.428	0.013	0.551	0.997
0.5~0.714	0.005	0.516	0.991
1	a+b=0.414		0.998
1.4~2	0.474	-0.035	0.999
2.33~5	0.412	0.004	0.991
6~14	0.426	-0.002	0.996

下载: 导出CSV 

| 显示表格

（6）内摩擦角

细观摩擦系数μ对宏观摩擦系数tanφ的影响如图 15所示。由图 15可知，tanφ随着μ的增大呈先增大后平缓的变化趋势；平均粒径较大的模型的tanφ略大于平均粒径较小的模型。基于式（26），对模拟结果进行拟合，如表 13所示。最后，根据方案Ⅰ的结果可知，当k_s/k_n=2~6时，k_s/k_n对tanφ的影响较小。

图  15  宏观摩擦系数tanφ vs接触摩擦系数μ(方案H)

Figure  15.  tanφ vs μ (Plan H)

下载: 全尺寸图片 幻灯片

表  13  系数b_S，c_S，d_S的拟合结果(方案H)

Table  13.  Fitting results of coefficients b_S, c_S and d_S (Plan H)

颗粒几何参数			拟合系数及R2
Rmin /mm	Rmax/Rmin		bS	cS	dS	R2
0.15	1.67		0.6482	0.8027	2.1002	0.639
0.2	2		0.6644	0.1074	0.4247	0.723
0.3	3		0.8122	0.9226	2.3140	0.968

下载: 导出CSV 

| 显示表格

3.   半解析查表法（SLT法）

半解析查表法（semi-analytic lookup table，简称SLT）为基于CBM模型的颗粒流材料细观参数的标定方法。本方法基于巴西劈裂试验、单轴压缩试验及直剪试验的试验结果，对类岩石材料的细观参数进行标定。具体步骤如图 16所示。

图  16  SLT法的具体步骤

Figure  16.  Specific steps of SLT method

下载: 全尺寸图片 幻灯片

4.   应用

本文采用某滑坡的中风化砂岩来验证SLT法，该砂岩的室内试验力学特性如表 14所示。根据SLT法：

表  14  砂岩的力学特性及细观参数的选取

Table  14.  Mechanical properties of sandstone and selection of microscopic parameters

	细观参数取值								巴西试验			单轴压缩试验									三轴试验强度/MPa
	Rmin/ mm	Rmax/ Rmin	σcn/ MPa	ks/kn	Ec/ GPa	τcs/ MPa	µ		σt/MPa			fc/MPa			E/GPa			$\nu$			围压1 MPa			围压6 MPa
									值	误差		值	误差		值	误差		值	误差		值	误差		值	误差
砂岩	—	—	—	—	—	—	—		5.99	—		64	—		9.195	—		0.226	—		71.89	—		95.45	—
1号	0.2	2	34.896	2.172	14.165	69.792	1.7		4.77	-20.4		50.69	-20.8		9.005	-2.1		0.239	+5.8		—			—
2号	0.2	2	51.235	2.172	14.165	51.235	1.7		6.96	+16.2		60.35	-5.7		9.005	-2.1		0.239	+5.8		—			—
3号	0. 2	2	45	2.172	14.165	90	1.7		6.58	+9.8		67.31	+5.2		8.972	-2.4		0.234	+3.5		72.55	+0.9		91.18	-4.5

下载: 导出CSV 

| 显示表格

（1）颗粒粒径取R_min=0.2mm，R_max/R_min=2。

（2）根据表 11，得${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$为1，2时，系数h_U分别为1.871，1.983，在${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$未确定时，h_U取平均值1.927，将$\nu$=0.226代入式（19），得kk=2.172。

（3）根据表 10，可得e_U=0.422，f_U=2.453，g_U=2.274，将RR=1/3，kk=2.172代入式（16），得E_c=14.165 GPa。

（4）基于直剪试验获得的φ值与三轴试验获的有差距，无法对μ直接进行标定，暂定μ=1.7。

（5）参考表 8，当${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$=2，估算a_U=1.635、b_U=1.448，代入式（12），得${\sigma _{{\text{cn}}}}$=34.896，则${\tau _{{\text{cs}}}}$=69.792；当${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$=1，估算c_U=3.095、d_U=0.257，代入式（13），得${\tau _{{\text{cs}}}}$=51.235，则${\sigma _{{\text{cn}}}}$=51.235。

由此，确定了两组细观参数，为表 1中的1，2号参数，模拟结果如表 14所示，表 14中误差的单位为%。结果显示，模量与泊松比的拟效果好，误差在6%以内。发现，1号参数中（${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$=2），颗粒流材料的压拉强度比f_c/${\sigma _{\text{t}}}$与试的一致。而在前述理论中，f_c由${\tau _{{\text{cs}}}}$，${\sigma _{{\text{cn}}}}$决定，${\sigma _{\text{t}}}$由${\sigma _{{\text{cn}}}}$决定，故建议在3号参数中保持${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$=2。1号参数中，${\sigma _{\text{t}}}$，f_c均偏小20%，表明${\sigma _{{\text{cn}}}}$，${\tau _{{\text{cs}}}}$均偏小；2号参数中，${\sigma _{\text{t}}}$偏大，f_c较好，表明${\sigma _{{\text{cn}}}}$偏大，既34.896 < ${\sigma _{{\text{cn}}}}$ < 51.235。此时，在3号参数中，基于式（11），调整${\sigma _{{\text{cn}}}}$为45 MPa。由前述分析可知，${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$=2时，f_c由${\sigma _{{\text{cn}}}}$决定。此时，若将${\sigma _{{\text{cn}}}}$由51.235 MPa调整至45 MPa，在${\tau _{{\text{cs}}}}$不变（51.235 MPa）时，f_c必下降，既更加偏小。故，此时${\tau _{{\text{cs}}}}$应大于51.235 MPa。结合上述建议保持${\tau _{{\text{cs}}}}$/${\sigma _{{\text{cn}}}}$=2，则取${\tau _{{\text{cs}}}}$=90 MPa，并进行三轴试验以验证。结果表明，3号参数中，巴西抗拉强度模拟值偏大，其余参数模拟效果较好。

5.   结论

本文基于理论推导与数值模拟的结果，建立了适用于颗粒流CBM模型的半解析查表法。结果表明：

（1）代表性单元体与实际颗粒流模型具有相同的宏、细观参数变化规律。

（2）各系数的理论值与拟合值基本相符，但部分与刚度比相关的系数呈现较大差异。

（3）半解析查表法包含两部分内容，一是文中具体的参数标定步骤；二是在参数微调过程中，利用宏-细观参数的内在关联判断调整方向，这些内在关联可通过理论分析的结果获得。

（4）半解析查表法提出了明确的参数标定顺序及标定方式。应用结果表明：半解析查表法能够快速且较为准确的获得岩石材料的细观参数值。