基于中心组合设计的颗粒流平直节理模型宏-细观参数相关性研究

2021-04-16

中南大学学报（自然科学版） 2021年3期

(中南大学土木工程学院，湖南长沙，410075)

在矿产资源开采、边坡、隧道、石油开采等岩体工程中，天然岩石材料存在的细观缺陷(如裂纹或孔隙)极大地影响着岩石的宏观力学性能。在外载作用下，裂纹或孔隙的演变发展形成宏观缺陷会引起岩体工程灾害隐患，因此，研究岩石细观缺陷发展至宏观破坏过程及其机理，对岩体工程的强度计算、安全评估和防灾减灾等具有重要理论和实际意义。由于现有试验测试手段的局限性和岩石内部缺陷的复杂性，人们难以在岩体工程中实时监测到岩石细观损伤积累至宏观破坏的全过程，故多借助于数值模拟方法。目前，研究岩石裂纹扩展的数值模拟方法主要有离散元法、有限元法、边界元法等，其中，基于非连续介质力学理论的离散元法在模拟岩石裂隙扩展时具有独特的优势，因为该方法既能避免选取复杂的本构关系和裂尖奇异性，又能够从细观层面上模拟岩石宏观力学特性并揭示其破坏机理。离散元颗粒流软件PFC(particle flow code)在岩石类材料的数值模拟中应用较广泛[1]。PFC 软件是将岩石离散为刚性颗粒，通过定义颗粒接触本构关系来描述颗粒之间的相互作用，以此模拟岩石的宏观力学特性。

目前，PFC 颗粒黏结模型大多为接触黏结(linear contact bond)、平行黏结(parallel bond)和平直节理模型(flat jointed model)。在接触黏结模型中，黏结断裂时颗粒接触刚度仍然存在，这与岩石材料不相符[2]；在平行黏结模型中，模拟得到的岩石抗拉强度与抗压强度的比值偏大，与岩石材料的力学性能也不相符[3-4]，通常采用簇平行和新的胶接模型来解决该问题[5-6]。在平直节理模型中[7]，圆形颗粒被构造成多边形从而抑制了颗粒的旋转，使得该模型计算得到的单轴抗拉强度和抗压强度之比减小，更符合岩石力学特性。采用颗粒流PFC2D软件模拟岩石的宏观力学特性时，首先需要确定颗粒黏结模型的细观参数。目前，采用较多的是试错法，如文献[8-9]基于试错法标定了平行黏结的细观参数；文献[10]基于试错法标定了平直节理接触模型的细观参数，并与平行黏结模型的细观参数进行对比，分析了平直节理模型的优势，但未建立宏-细观参数之间的关系式。由于试错法具有一定的盲目性，细观参数较多且标定过程较繁琐，从而导致试验次数很多，因此，研究者寻求更有效的细观参数定量确定方法，如：文献[11]采用单因素和多因素分析方法，分析了平行黏结模型的细观参数对宏观参数的影响，建立了平行黏结模型的宏-细观参数关系式；文献[12]采用BP神经网络法，建立了簇平行黏结模型的宏-细观参数的非线性模型，但BP 神经网络需要大量的随机组合样本才能具有较高的精度；文献[13]采用球面对称设计法建立了平行黏结宏-细观参数的关系式，但考虑的宏观参数有限；文献[14-15]通过对平直节理模型细观参数进行正交分析，得到了细观参数对宏观参数的影响规律，并利用回归分析建立了宏-细观参数关系式，但因试验次数较少导致精度较低。总之，目前关于PFC2D平直节理接触模型的细观参数定量确定方法偏少。为提高精度、减少试验次数，需要寻求一种更合适的方法来建立PFC2D平直节理接触模型的宏-细观参数关系式。中心组合设计(又称星点设计，central composite design)[16-17]是一种多因素五水平的试验设计方法，相对于正交试验和均匀化试验，它在提高精度的同时不需要大量增加实验次数，特别适用于多因素、多水平试验情形。

本文作者基于PFC2D平直节理模型，通过控制变量法得到细观参数的取值范围和对宏观参数的影响规律，采用中心组合设计和回归分析，建立岩石PFC2D宏-细观参数关系式以定量确定细观参数，模拟计算岩石在单轴和双轴压缩下完整岩石试件的应力-应变曲线和含多裂纹岩石试件的断裂轨迹，并与试验结果进行对比，以验证所建立的岩石PFC2D宏-细观参数关系式的有效性。

1 PFC平直节理模型

平直节理模型是针对硬岩提出的一种模型，如图1所示。该模型将圆形颗粒构造成多边形，从而限制了颗粒的旋转，致使模型单轴抗压强度和抗拉强度的比值增大，更加符合岩石类材料的力学性能。平直节理模型采用扁平面将2 个颗粒黏结，黏结面由片段的单元组成。每个单元具有黏结和未黏结2种模式，当黏结界面上的荷载单元达到黏接强度时，黏结单元便会破坏成非黏结单元，其本构关系如图2所示。表1所示为平直节理模型的细观参数，包括颗粒细观参数(6个)和平直节理参数(7个)。

图1 PFC2D平直节理模型[7]Fig.1 Flat joint model of PFC2D[7]

图2 PFC2D平直节理模型的本构模型[7]Fig.2 Constitutive relation of PFC2D flat joint model[7]

2 PFC细观参数确定方法

2.1 细观参数取值范围的确定

为了模拟岩石类材料宏观力学特性，首先需要确定PFC2D平直节理模型的细观参数。由于平直节理模型的细观参数较多(共13个)，为了减少标定次数，根据前人的研究成果[13]，进行以下假设：

表1 PFC2D平直节理模型的细观参数[7]Table 1 Microscopic parameters of PFC2D flat joint model[7]

1)颗粒接触模量Eb等于平直节理接触模量Ec；

2)颗粒刚度比(kn/ks)b等于平直节理接触刚度比kn/ks；

3)最大、最小颗粒半径之比Rmax/Rmin=1.66；

4)颗粒密度ρ=2 300 kg/m3；

5)平直节理半径系数-λ=1；

6)交界面段数N=4；

7)平直节理抗拉强度小于抗剪强度。

为了确定PFC2D细观参数，需要分析细观参数对岩石材料宏观力学参数(含弹性模量E、泊松比ν、抗压强度σf、抗拉强度σt、黏聚力C2、摩擦角φ)的影响，以此确定细观参数的范围。针对剩余的7个细观参数采用控制变量法分析，选取适当的初始细观参数，即Ec=40 GPa，kn/ks=1，C1=30 MPa，σc=3 MPa，μb=0.1，μ=0.1，R=0.2 mm[18]，在改变单一细观参数的情况下开展岩石材料的单轴压缩(E，ν，σf)、单轴拉伸(σt)和双轴压缩(C2，φ)数值模拟试验，计算得到的单一细观参数对宏观力学参数的影响结果如图3～9 所示，其中横坐标为PFC2D细观参数，纵坐标为归一化处理后的岩石宏观参数，即细观参数改变后通过模拟计算得到的宏观参数与由初始细观参数模拟计算得到的宏观参数之比。

从图3可知：接触模量Ec对宏观弹性模量E影响最大，对其他宏观参数几乎没有影响(在图中表现为其他参数曲线基本重合)。由图4可知：刚度比kn/ks主要影响泊松比ν，随着kn/ks的增加，宏观弹性模量E、抗拉强度σt略减少(在图中表现为E，σt和C2曲线基本重合)，其他宏观参数基本保持不变。由图5可知：平直节理黏聚力C1主要影响σf和黏聚力C2，对泊松比ν和摩擦角φ的影响较小。由图6可知：平直节理抗拉强度σc主要影响σt，随着σc增加，σf略增加。由图7和图8可知：颗粒摩擦因素μb和平直节理内摩擦因数μ主要影响摩擦角φ，且μb对φ的影响更显著；随着μ增加，σf略微增加，C2略微减少，对其他宏观参数影响较小(在图中表现为E、σt曲线基本重合)。由于μb和μ对宏观参数的影响具有相同的规律，可取μ=μb[7]。由图9可知，颗粒半径R对宏观参数影响的幅度较小，这与文献[19]中的结论一致。考虑到计算效率和精度，本文选取最小半径Rmin=0.3 mm。因此，待定的独立细观参数只有5个：Ec，kn/ks，C1，σc和μb。

图3 细观参数Ec对宏观参数的影响Fig.3 Effect of microscopic parameters Ec on macroscopic parameters

图4 细观参数kn/ks对宏观参数的影响Fig.4 Effect of microscopic parameters kn/ks on macroscopic parameters

图5 细观参数C1对宏观参数的影响Fig.5 Effect of microscopic parameters C1 on macroscopic parameters

图6 细观参数σc对宏观参数的影响Fig.6 Effect of microscopic parameters σc on macroscopic parameters

图7 细观参数μb对宏观参数的影响Fig.7 Effect of microscopic parameters μb on macroscopic parameters

图8 细观参数μ对宏观参数的影响Fig.8 Effect of microscopic parameters μ on macroscopic parameters

图9 细观参数R对宏观参数的影响Fig.9 Effect of microscopic parameters R on macroscopic parameters

由上述分析可知，细观参数Ec，kn/ks，σc和μb与宏观参数关联度最大的分别是E，v，σt和φ，可由对应的宏观力学参数范围来确定该细观参数范围。细观参数C1与宏观参数σf和C2关联度均较大，考虑到岩石类材料的黏聚力C2差异较大，变化范围难以确定，故C1的取值范围由σf的取值范围来确定。对于大多数岩石类材料，宏观力学参数的大致范围为[20-21]：E=30～120 GPa，v=0.10～0.35，σf=50～300 MPa，σt=3～25 MPa，φ=30°～65°。故由图3～9 可大致确定细观参数的取值范围为：Ec=35～110 GPa，kn/ks=1.6～5.3，C1=25～150 MPa，σc=5.3～30 MPa，μb=0.1～0.7。

2.2 基于中心组合设计样本的细观参数计算公式

中心组合设计又称星点设计，是一种多因素多水平的试验设计方法，它是在二水平析因设计方法的基础之上加上极值点和中心点组成的。可见，中心组合设计相对于正交、均匀设计具有较高的实验精度，并且也不过多地增加实验次数，更适合于多水平多因素的试验设计。为了定量确定PFC2D细观参数，采用中心组合设计法建立细观参数与宏观参数之间的关系式。考虑到待定的细观参数有5个(Ec，kn/ks，C1，σc和μb)，利用Design-Expert 软件生成5 因素5 水平的中心组合设计代码表(表2)，X1，X2，X3，X4和X5分别代表Ec，kn/ks，C1，σc和μb。每个元素均有5 个水平代码：-2.38，-1，0，+1 和+2.38，其中，代码-2.38，0 和+2.38对应的值分别为该因素范围的最小值、平均值、最大值，代码-1和+1对应的值可通过“任意2个水平代码之差的比值与该水平代码对应值之差的比值相等”原则计算求得。以Ec为例，其范围为35～110 GPa，则Ec的水平代码-2.38，0和+2.38所对应的值分别为35，72.5和110。设其水平代码-1和+1对应的值分别为x和y，根据上述原则，有

经求解得到x=56.74，y=88.26。同理，可求得其他因素水平代码对应的值，见表3。将表3中各因素水平代码对应的值代入表2，可得到中心组合设计试验表，如表4所示。

根据上述控制变量法分析(图3～9)、中心组合设计试验表(表4)可知：E受Ec和kn/ks的影响，v主要受kn/ks的影响，σf受C1，σc和μb的影响，σt受σc和kn/ks的影响，C2主要受C1和μb的影响，φ主要受C1和μb的影响。对6个数值模型宏观力学性能(E，v，σf，σt，C2，φ)和5 个细观参数(Ec，kn/ks，C1，σc，μb)，利用SPSS软件进行非线性回归，得到模型宏观力学参数与PFC2D细观参数之间的关系式：

表2 中心组合设计代码表Table 2 Central composite design code table

表3 细观参数的中心组合设计代码值Table 3 CCD code value of microscopic parameters

其中，R2为相关系数。可见，相关系数R2均较高，表明拟合精度较高，能够准确地反映岩石宏-细观参数的关联性。若已知岩石的宏观力学参数，则由式(2)可计算出PFC2D细观参数，其中Ec，kn/ks和σc由E，v和σt的关系式联立求解，C1和μb通过相关系数R2较高的σf和φ关系式求得。

3 试验验证

3.1 岩石单轴和双轴压缩试验

为了验证宏-细观参数关系式(2)，选取产自于云南的红砂岩开展试验研究。表5所示为本试验测得的红砂岩以及其他4 种岩石材料[22-23]的宏观力学参数。在TRW-30001 轴伺服试验机上开展岩石的单轴(直径×高度为50 mm×100 mm的圆柱体试件)、双轴(长×宽×高为100 mm×100 mm×40 mm 的立方体试件)压缩试验，其中双轴压缩中侧向压力为6 MPa，加载速率为500 N/s。在试验过程中自动记录应力-应变曲线。

根据不同岩石的宏观力学参数(表5)和模型宏观力学参数与PFC2D细观参数之间的关系式，可分别计算得到相应的PFC2D细观参数(Ec，kn/ks，C1，σc和μb)，如表6所示。采用该细观参数，分别开展不同岩石的单轴压缩(E，ν，σf)、单轴拉伸(σt)和双轴压缩(C2，φ)模拟计算，得到的岩石宏观力学参数(E，v，σf，σt，C2和φ)见表7，并与试验结果进行对比。从表7可见：各参数的实验结果与模拟结果相对误差均小于10%，模拟精度较高，从而验证了采用中心组合设计建立的模型宏观力学参数与PFC2D细观参数之间的关系式是合理、有效的。

表4 中心组合设计试验表Table 4 Central composite design table

表5 不同岩石的宏观力学参数Table 5 Macroscopic mechanical parameters of different rocks

表6 不同岩石的PFC2D细观参数计算结果Table 6 Calculated results of PFC2D microscopic parameters of different rocks

图10所示为红砂岩在单轴、双轴压缩下数值模拟计算得到的应力-应变曲线，并与试验结果进行对比。从图10可见：在单轴、双轴压缩下，实测的红砂岩应力-应变曲线(OABCD)可分为压密阶段(OA)、弹性变形阶段(AB)、损伤阶段(BC)和破坏阶段(CD)。模拟结果(OB′C′D曲线)中没有出现试验曲线的压密阶段，这是因为PFC2D建模时颗粒为连续体，而实际岩石存在有天然的孔隙、裂纹等缺陷，在压缩时会有微孔、微裂缝等闭合压密阶段(OA)。计算得到的红砂岩单轴压缩弹性模量、抗压强度分别为11.38 GPa 和59.45 MPa，实测的单轴弹性模量、抗压强度分别为10.69 GPa 和57.20 MPa，相对误差分别为6.45%和3.93%；计算得到的双轴压缩弹性模量和抗压强度分别为10.83 GPa和78.20 MPa，实测的双轴弹性模量和抗压强度分别为11.82 GPa 和82.40 MPa，相对误差分别为8.38%和5.10%。可见，弹性模量和抗压强度模拟结果(未考虑初始缺陷)均略高于实验结果，双轴压缩比单轴压缩的抗压强度更高。忽略初始压密阶段，单、双轴压缩试验曲线与模拟曲线基本重合，弹性模量、抗压强度均较吻合，从而验证了采用宏-细观参数关系式(2)计算PFC2D 细观参数方法的合理性和有效性。

表7 不同岩石的宏观力学参数(模拟与试验结果对比)Table 7 Macroscopic mechanical parameters of different rocks(Comparison of simulation and test results)

图10 红砂岩单、双轴压缩应力-应变曲线Fig.10 Stress-strain curves of sandstone

3.2 岩石多裂纹断裂试验

图11 含3条平行裂纹的红砂岩试件模型Fig.11 Model of sandstone specimen with three parallel cracks

为了进一步验证宏-细观参数关系式(2)，对含3条平行裂纹(裂纹长度为20 mm，宽度为1 mm)的红砂岩试件进行断裂试验。采用与上述压缩试验相同的红砂岩，宏观力学参数、细观参数与表5和表6所示的相同。试件是长×宽×高为100 mm×100 mm×20 mm的立方体，预制3条等长平行裂纹AB，CD和EF，裂纹长度为20 mm，宽度为1 mm，倾角α=45°，间距为15 mm，如图11所示。在DNS100 伺服加载系统上进行单轴压缩断裂试验，最大荷载为1 000 kN，以位移加载速率为0.1 mm/s进行加载，采用摄像机记录试件断裂破坏过程。

图12和图13所示分别为单轴压缩下红砂岩试件细观裂纹演变和宏观裂纹扩展轨迹的数值模拟结果，其中，黑色、红色、白色分别代表颗粒、细观裂纹、宏观裂纹。从图12和图13可见：在初始加载时，3 条裂纹的6 个裂纹尖端都会产生初始损伤而出现微裂纹(图12(a))，此时没有宏观裂纹(图13(a))；随着进一步加载，各裂纹尖端的损伤不断积累会产生越来越多的微裂纹并集聚(图12(b))，由于3条裂纹之间的损伤区相互作用，裂纹之间的微裂纹首先贯通而形成了最初的2 条宏观裂纹BC和DE(图13(b))。当外载不断增加到破坏荷载时，上部裂纹AB的尖端A和下部裂纹EF的尖端F附近微裂纹急剧增多(图12(c))，最后形成了2条沿轴向加载方向的宏观主裂纹(尖端A向上延伸、尖端F向下延伸)以及衍生的次生裂纹，破坏轨迹如图13(c)所示。

图12 单轴压缩下含3条平行裂纹的红砂岩试件细观裂纹演变过程(模拟结果)Fig.12 Micro-crack evolutions of sandstone specimen with three parallel cracks under uniaxial compression(simulation results)

图13 单轴压缩下含3条平行裂纹的红砂岩试件宏观裂纹扩展轨迹(模拟结果)Fig.13 Macro-crack trajectories of sandstone specimen with three parallel cracks under uniaxial compression(simulation results)

图14所示为摄像机记录的单轴压缩下红砂岩试件裂纹扩展轨迹试验结果。从图14可见：在加载初期，各裂纹尖端会有局部损伤产生微裂纹，但没有出现宏观裂纹(图14(a))；随着外载增加，3条裂纹之间的微裂纹聚集、演变、贯通而产生了2条宏观裂纹；当外载不断增加到破坏荷载时，裂纹尖端A和尖端F的微裂纹不断发育，最后形成了沿轴向方向延伸的2 条主裂纹(尖端A向上延伸、尖端F向下延伸)和衍生的次生裂纹，破坏轨迹如图14(c)所示，其中，试件的上部碎片是试验机卸载时能量释放冲压所致。

比较图13和图14可见：除了最后破坏阶段的次生裂纹有细微差别外(这是实际岩石材料的非均质性所致)，单轴压缩下含3 条平行裂纹的红砂岩试件裂纹扩展轨迹试验结果与PFC2D模拟结果较吻合，从而进一步验证了采用宏-细观参数关系式(2)计算PFC2D细观参数方法的合理性和可行性。

图14 单轴压缩下含3条平行裂纹的红砂岩试件裂纹扩展轨迹(试验结果)Fig.14 Propagation trajectories of sandstone specimen with three parallel cracks under uniaxial compression(test results)