考虑断裂韧度的PFC3D细观参数标定新方法

谭攀，饶秋华，李卓，张权，易威

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

在矿产开采、隧道开挖、水利运输等岩体工程中，岩石的非连续性(如裂纹、孔洞)会极大地影响岩石的力学性能。在外载作用下，岩石中裂纹容易起裂、扩展、贯通而导致断裂，从而引起工程安全隐患问题，因此，研究岩石破坏机理对岩体工程的安全评估、防灾减灾等具有重大的理论和实际意义。由于理论分析三维多裂纹扩展过程较复杂，实验测试也难以揭示三维多裂纹扩展的细观机理，因此，数值模拟成为一种分析含裂隙岩石的力学性能和破坏机理的新途径。目前，数值模拟方法主要有连续介质方法(有限元、扩展有限元、边界元、流行元等)和非连续介质方法(块状离散元、颗粒流离散元等)。传统的连续介质方法(如有限元)[1]由于存在网格的重划分和裂纹尖端的奇异性等问题，在处理裂隙起裂和延伸时具有一定的局限性。非连续介质方法(如颗粒流离散元)[2]将系统理想化为一系列分离单元(单元在接触点处黏结)，以单元间接触的破坏直接表达损伤，无需对裂纹扩展进行假设，能自然地模拟裂纹的延伸问题，其中颗粒流软件PFC 应用较广[2−4]。PFC 软件是基于牛顿第二定律、定义颗粒间的接触本构关系来模拟材料的力学行为。常见的黏结模型主要有3 种：接触黏结模型(contact bonded model，CBM)、平行黏结模型(parallel bond model，PBM)和平直节理模型(flat-joint model，FJM)。由于CBM和PBM 在模拟岩石力学特性时[2,5−6]，存在岩石的单轴抗压强度σf与抗拉强度σt的比值σf/σt(3～7)远远低于常值(10～20)、内摩擦角偏低、强度包络线仅为直线等问题，采用簇平行模型和新的胶接模型(修改本构)[6−7]可部分解决上述问题或仅能解决二维问题，而FJM 通过构造多边形抑制颗粒的旋转，从而克服了CBM和PBM的不足，能更加真实地模拟岩石的力学行为。在PFC 模拟过程中，首先需要确定细观参数，通常采用试错法、神经网络法、试验设计法(Plackett-Burman 设计、正交设计、球面正交设计等)等方法，其中应用最多的是试错法。王刚等[8−10]采用试错法标定了PFC2D 的PBM和FJM 细观参数，但试错法标定繁琐，具有盲目性。周喻等[11−12]采用BP 神经网络和优化的BP网络对PBM 细观参数进行了反演，其与试错法相比能快速标定细观参数，但需要大量样本或较小的细观参数范围才能得到较高精度，且无法得到宏−细观参数关系式。为了得到PFC 宏−细观参数关系式，YOON[13]采用PB设计法分析了CBM细观参数对宏观参数的敏感性，建立了CBM 的宏−细观参数表达式；彭霞等[14]用球面正交设计建立了PBM 的宏−细观参数表达式；CHEN[15]用正交设计法建立了FJM的宏−细观参数表达式，但他们考虑的宏观参数有限。SUN 等[16]用PB 设计法分析了FJM 细观参数对宏观参数的敏感性，筛选了主因素，用中心组合设计法(CCD)建立了宏−细观参数表达式，但只考虑了强度参数，未考虑断裂参数。而裂纹的扩展与断裂韧度密切相关，有少量文献考虑了I型断裂韧度KIC，如王勇[17]用支持向量机模型对PBM的KIC进行了标定，但未涉及II型断裂韧度KIIC的标定，且没有宏细观参数表达式。目前，关于PFC3D 细观参数定量确定方法较少，急需一种同时考虑强度和断裂参数的PFC3D宏−细观参数关系式的方法。为此，本文作者基于PFC3D 平直节理模型，同时考虑强度参数和岩石I型拉伸断裂韧度KIC和II型剪切断裂韧度KIIC，根据其宏观力学性能和单因素法确定细观参数范围；采用PB+CCD 试验设计相结合的方法，提出一种综合考虑强度参数和断裂韧度的岩石PFC3D 细观参数标定的新方法，建立岩石PFC3D 宏−细观参数关系式，对单轴和三轴压缩下完整岩石试件的应力−应变曲线和含裂纹岩石试件的断裂轨迹进行模拟计算，并与实验结果进行比较，从而验证该方法的合理性和有效性。

1 PFC3D平直节理模型

PFC3D平直节理模型(flat-joint model，FJM)是由POTYONDY[18]提出的一种新的接触模型，它的基本原理如图1所示。物体由n个多边形颗粒组成，当2 个颗粒表面的最小距离小于g0(初始表面间隙)时，自动形成平直节理接触。颗粒之间的接触界面由片段的单元组成(图1(a))，每个单元有黏结和未黏结2种形式。当黏结单元上的法向应力超过平直节理抗拉强度σc时，黏结单元会发生拉伸破坏而成为非黏结单元，此时，单元的剪切强度τc由库仑准则确定：

图1 PFC3D平直节理接触模型[18]Fig.1 Flat-jointed contact model of PFC3D[18]

式中：cb和φb分别为平直节理黏聚力和内摩擦角；为单元上的法向应力。当黏结单元上的切向应力超过抗剪强度τc时，黏结单元会发生剪切破坏而成为非黏结单元，此时，单元的残余剪切强度为

式中：μc平直节理摩擦因数。传统的接触黏结模型(CBM)和平行黏结模型(PBM)均采用圆形颗粒，且g0=0(无间隙接触)，颗粒破坏后仍可发生自由旋转，这与实际情况不符。平直节理模型(FJM)因采用多边形颗粒，且g0可以大于0(即认为颗粒之间允许间隙接触)，增大了颗粒接触数与颗粒数的比值，能有效抑制黏结破坏后颗粒的旋转，从而提高了颗粒间的互锁效应，因此，模拟得到的抗压强度σf与抗拉强度σt之比σf/σt更接近于实验结果，能更好地反映材料真实的力学性能。

PFC3D平直节理模型(FJM)的细观参数见表1，包括颗粒参数和节理参数共15 个。进行如下假设[2,18−19]：1)Eb=Ec，(kn/ks)b=kn/ks，μb=uc；2)λ=1；3)Rmax/Rmin=1.66；4)ρ=2 690 kg/m3；5)Nr=1，Nα=3。简化后的8个细观参数为Ec，kn/ks，σc，cb，φb，gr，μc和Rmin。

表1 PFC3D平直节理接触模型细观参数Table 1 Micro-parameters of PFC3D flat-joint model

2 考虑断裂韧度的PFC3D 细观参数定量确定方法

基于PFC3D 简化后的8 个细观参数，首先，采用单因素分析法获得细观参数对宏观参数的影响规律，依据宏观力学参数范围确定细观参数范围；然后，采用PB+CCD 试验设计相结合的方法建立岩石PFC3D宏−细观参数关系式；最后，采用遗传算法得到PFC3D细观参数的优化解。

2.1 考虑断裂韧度的PFC3D细观参数范围

以岩石为例，其宏观力学参数范围一般为[20−23]：单轴抗压强度σf=70～250 MPa，抗拉强度σt=4～18 MPa，弹性模量E=30～80 GPa，泊松比υ=0.1～0.4，I 型断裂韧度KIC=0.5～2.0 MPa·m1/2，II 型断裂韧度KIIC=1.5～6.0 MPa·m1/2。选取初始的8个细观参数为：Ec=50 GPa，kn/ks=2.0，σc=10 MPa，cb=120 MPa，φb=10°，gr=0.3，μc=0.3，Rmin=1.0 mm。采用单因素分析，每次只改变1个细观参数，其他7个细观参数保持不变，分别模拟计算：单轴压缩试件(图2(a)，直径×高为50 mm×100 mm)的σf，E和υ；直接拉伸试件(图2(b)，直径×高为50 mm×100 mm)的σt；带裂纹中心圆孔圆盘压缩试件(图2(c)，直径×高为76 mm×30 mm，裂缝长度为12 mm，中心孔直径15 为mm)的KIC；单边平面压剪试件(图2(d)，长×宽×高为50 mm×50 mm×50 mm，裂缝长度30 mm)的KIIC共6 个宏观参数。其中，模拟计算得到的直接拉伸试验抗拉强度σt比巴西劈裂试验计算的σt精度更高[15,19]。带裂纹中心圆孔圆盘(HCCD)模拟计算得到的I 型断裂韧度KIC比人字形切槽巴西圆盘试样(CCNBD)的KIC更接近国际实验值[23−24]。

图2 数值模拟试件图Fig.2 Figures of numerical simulation specimen

平直节理内摩擦角φb主要影响岩石内摩擦角φ和完整岩石HB (Hoek-Brown)强度准则参数mi[19]，而安装间隙比gr对应着岩石的配位数(coordination number,CN)[26]，因此，本文通过φ和mi共同确定φb，通过CN确定gr。φb对内摩擦角φ和H-B强度参数mi的影响如图3所示。从图3可见φ和mi都随φb的增加显著增加。一般岩石的φ为30°～55°，mi为5～20[25]，综合确定φb的范围为10°～30°。gr对配位数CN 的影响如图4所示。从图4可见：当gr为0，0.2，0.4，0.5，0.6 时，对应的CN 分别为5.00，8.08，9.32，9.64和9.65；当gr大于0.5 时，CN 趋于稳定，配位数(CN)的范围为6～10[26]，最终确定gr的范围为0.2～0.5。

图3 细观参数φb对内摩擦角φ和H-B强度参数mi的影响Fig.3 Effect of micro-parameters φb on φ and mi

图4 细观参数gr对配位数的影响Fig.4 Effect of micro-parameters gr on coordination number

图5所示为PFC3D 模拟计算得到的细观参数对宏观参数的影响，其中细观参数均进行了归一化处理。从图5可见：平直节理黏聚力cb和安装间隙比gr对单轴抗压强度σf影响较大(图5(a))；由于gr范围已确定，所以，由σf的范围推得cb的范围为35～220 MPa；cb和gr间存在交互作用，需要对cb的范围进行调整，本文将cb的范围调整为50～120 MPa。此外，平直节理抗拉强度σc对抗拉强度σt的影响显著且呈正相关(图5(b))，由σt的范围可推得σc的范围为4～22 MPa。平直节理刚度比kn/ks对泊松比υ影响显著且呈正相关(图5(c))，由υ的范围推得kn/ks的范围为1.5～2.8。平直节理弹性模量Ec对弹性模量E的影响最大且呈正相关(图5d)，由E的范围推得Ec的范围为35～90 GPa。平直节理摩擦因数μc对岩石宏观力学参数基本无影响(图5(a)～5(f))，一般取0.5[2]。为了分析μc与其余细观参数是否具有交互影响，本文选取μc的范围为0.3～0.7。

从图5(e)可知：平直节理抗拉强度σc对I 型断裂韧度KIC影响最大，且KIC随σc增加而增加。由于σc对抗拉强度σt影响显著(图5(b))，因此，σc对KIC和σt的影响也最大。文献[15，27]中，σc的范围只通过σt求得，忽略了KIC的影响，难以准确模拟裂隙岩石裂纹拉伸扩展过程，ZHANG[28]也表明KIC与σt是相关的。本文同时考虑KIC和σt对σc的影响来确定σc的范围，由KIC的范围推得σc的范围为6～27 MPa，最终结合σt对σc的影响共同确定σc的范围为6～22 MPa。

从图5(f)可知：平直节理黏聚力cb、安装间隙比gr和平直节理刚度比kn/ks都对II 型断裂韧度KIIC影响较大；KIIC随cb和gr增大而增大，随kn/ks增大而减小，但KIIC对cb的影响最显著，且gr和kn/ks的范围已确定，所以，由KIIC的范围推得cb的范围为20～190 MPa。根据文献[28−29]，cb的范围仅通过σf求得，忽略了KIIC的影响，难以准确模拟裂隙岩石裂纹剪切扩展过程。本文同时考虑cb对KIIC和σf的影响来确定cb的范围，结合σb对KIIC的影响最终共同确定σb的范围为50～120 MPa。

由图5(a)～5(d)可知：当Rmin较小时(0.7～1.3 mm)，Rmin对宏观强度参数(σf，σt，E和υ)影响较小；Rmin与断裂韧度KIC和KIIC整体上呈正相关关系(图5(e)和图5(f))，这与POTYONDY等[2]得出的结论一致，因此，可依据KIC和KIIC确定Rmin的范围。综合考虑Rmin影响KIC和KIIC的最大范围为1.1～1.2 mm，最小范围为0.8～1.1 mm，最终确定Rmin的取值范围为0.8～1.2 mm。根据ASTM[30]规范要求，当岩石试件最小直径尺寸L与平均颗粒尺寸Daver之比L/Daver>10时，离散颗粒的尺寸对完整岩石结果(宏观强度参数)影响不大。本文模拟计算所采取的岩石试件最小直径L=50 mm，最大颗粒半径Rmax=1.66Rmin[2]，颗粒平均半径Daver=(2Rmax+2Rmin)/2=2.128～3.192 mm，其对应的L/Daver范围为15.7～23.5，超过10，因此，颗粒尺寸对宏观强度参数影响不大。但由图5(e)和5(f)可知颗粒粒径对断裂参数有较大影响，且呈正相关关系。

图5 细观参数对宏观参数的影响Fig.5 Effect of micro-parameters on macro-parameters

综上所述，细观参数的取值范围为：Ec=35～90 GPa，kn/ks=1.5～2.8，σc=6～22 MPa，cb=50～120 MPa，φb=10°～30°，gr=0.2～0.5，μc=0.3～0.7，Rmin=0.8～1.2 mm。

2.2 基于PB+CCD 方法的PFC3D 宏−细观参数关系式

每个细观参数对宏观参数都有一定影响，细观参数间可能还存在耦合作用。PB 设计是一种较好的筛选方法，它能筛选出对宏观参数影响较大的细观参数，但它得到的宏−细观参数关系式呈线性关系，不能考虑细观参数间的耦合项。因此，在PB设计的基础上，针对筛选出的主要细观参数进行中心组合设计(CCD)，通过考虑主要细观参数间的耦合作用，建立更加合理的宏−细观参数非线性关系式。

2.2.1 基于PB 设计的PFC3D 宏−细观参数线性关系式

考虑到待定的细观参数有8个，利用数据统计分析软件MinTab[31]生成8 因素2 水平的PB 设计代码表(表2)，其中，各因素(即细观参数)有2个水平代码即−1和1，分别代表该细观参数范围的最小值和最大值。以Ec为例，其范围为35～90 GPa，则水平代码−1和1 对应的值分别为35 GPa和90 GPa，其他细观参数同理可得。PB 设计试验表如表3所示。根据PB 设计试验表，选取每行中的8 个细观参数进行PFC3D 模拟计算，数值模拟结果如表4所示。从表4可见模拟得到的宏观力学参数范围如下：σf为65.107～243.045 MPa，E为24.763～97.764 GPa，υ为 0.079～0.517，σt为 3.848～21.07 MPa，KIC为0.471～2.583 MPa·m1/2，KIIC为1.920～7.920 MPa·m1/2，除了第4组的细观参数极端组合导致宏观参数偏高外，其他组宏观参数与岩石的宏观力学参数范围基本吻合(σf=70～250 MPa，E=30～80 GPa，υ=0.1～0.4，σt=4～18 MPa，KIC=0.5～2.0 MPa·m1/2，KIIC=1.5～6.0 MPa·m1/2)，说明本文选取的细观参数范围是合理的。

表2 PB设计代码表(细观参数代码值)Table 2 PB design code table(code values of micro-parameters)

表3 PB设计试验表(细观参数实际值)Table 3 PB design test table(actual values of micro-parameters)

表4 PFC3D模拟结果(宏观参数)Table 4 PFC3D simulation results(Macro-parameters)

利用MinTab 软件，对6 个宏观力学参数(σf，E，υ，σt，KIC和KIIC)和8个细观参数(Ec，kn/ks，cb，σc，φb，gr，μc和Rmin)进行线性回归(回归系数见表5)，得到的宏−细观参数线性关系式如下(其中细观参数为代码值)：

表5所示为PB 设计结果分析，其中P表征细观参数对宏观参数的影响程度。P越小，表示细观参数影响越显著。以P=0.05 为界，P≤0.05 说明影响显著，否则影响不显著。从表5可见：影响σf最大的3 个细观参数依次为cb，gr和σc，σf都会随着cb，gr和σc增加而增加(正相关)；影响σt最大的3个细观参数依次为σc，gr和kn/ks，且σt与σc和gr呈正相关，σt与kn/ks呈负相关；影响E最大的3 个细观参数依次为Ec，gr和kn/ks，且E与Ec和gr呈正相关，E与kn/ks呈负相关；影响υ最大的3 个细观参数依次为kn/ks，σc和cb，且υ与kn/ks(呈正相关)，因为kn/ks越大，相当于切向刚度比越小，更容易发生侧向变形，因此，泊松比υ会更大；影响KIC最大的3个细观参数依次为σc，cb和gr，KIC都会随着cb，gr和σc的增大而增大(正相关)，且拉伸断裂韧度KIC与细观抗拉强度σc关联度更大；影响KIIC最大的3 个细观参数依次为cb，σc和gr，KIIC与cb，gr和σc呈正相关，剪切断裂韧度KIIC与细观黏聚力cb关联度更大。影响各宏观参数最大的3个细观参数见表6(按影响程度从左到右依次降序排列)。

表5 PB设计结果分析Table 5 Analysis of PB design results

表6 影响各宏观参数最大的3个细观参数Table 6 Three micro-parameters with the greatest influence in each macro-parameter

2.2.2 基于CCD 的PFC3D 宏−细观参数非线性关系式

MinTab 生成的3 因素5 水平CCD 设计代码表见表7，其中，列代表试验次数，行代表细观参数。CCD中因素(即细观参数)有5个水平，分别为−1.682，−1，0，1和1.682，其水平−1和1 代表的实际值分别与PB 设计的水平−1和1 代表的实际值一致，其他水平可以通过等比例求得。以Ec为例，细观参数范围为35～90 GPa，35 GPa 对应水平−1，90 GPa对应水平+1。设水平−1.682和水平1.682代码对应的实际值分别为x和y，则通过比值相等求得

得x=16.25 GPa，y=108.75 GPa，其他细观参数实际值同理可得(见表8)。将表8中的水平代码实际值代入CCD 设计代码表(表7)中进行数值模拟，模拟结果见表9，其中X1，X2和X3分别表示对各自宏观因素影响最大的3个细观参数。例如，对于σf，X1，X2和X3分别为cb，gr和σc；对于σt，X1，X2和X3分别为σc，gr和kn/ks。

表7 CCD设计代码表Table 7 CCD design code table

表8 细观参数实际值Table 8 Actual value of micro-parameters

表9 中心组合设计模拟结果Table 9 Simulation results of CCD

选用MinTab 软件，考虑主要细观参数本身的一次项和二次项、主要细观参数间的耦合影响项(二次项)，对表9中数据进行非线性回归分析，结果如表10所示。从表10可见：主要细观参数P≤0.05，与PB 设计的筛选结果一致；在抗压强度σf的影响因素中，cb·gr的P=0<<0.05(显著交互影响)，与3 个主要细观参数的影响程度一致(P=0)。由于cb·gr项的系数为正，其交互作用会显著增大σf，因此，在确定cb的范围时考虑了cb和gr的交互作用。从表4可见：模拟计算得到σf范围为64.107～243.045 MPa，与一般岩石的范围50～250 MPa 相符，表明考虑交互作用得到的cb范围是合理的；在影响泊松比υ的因素中，虽然(kn/ks)·σc的P=0.008，但远小于主要细观参数kn/ks和σc的影响程度，因此，在确定kn/ks的范围时未考虑其交互作用，同理，在弹性模量E的因素中，虽然Ec·(kn/ks)和Ec·gr的交互作用分别会使E减少和增大，但它们的影响程度P远小于主要细观参数Ec，kn/ks和gr的影响程度。因此，在确定Ec的范围时，可对Ec·(kn/ks)和Ec·gr的交互作用忽略不计。

回归分析得到的宏−细观非线性关系式如式(5)所示，各项系数见表10，其中为了避免与式(3)混淆，式(5)中的宏观参数用上标“′”表示，同样，细观参数也为代码值。

表10 中心组合设计结果分析Table 10 Result analysis of CCD

2.3 基于PB+CCD 方法的PFC3D 细观参数的优化解

采用PB 方法得到宏−细观线性关系式，采用CCD法得到宏−细观非线性关系式。为了求得与宏观参数(σf，E，υ，σt，KIC和KIIC)对应的细观参数优化解，建立如下优化模型的细观参数。

1)自变量。每个细观参数都有范围，在PB设计中，每个细观参数的代码值范围为[−1，1]，在CCD 中每个细观参数的代码值范围为[−1.682，1.682]。为了获得更接近实验值的解，选取较大的代码值范围[−1.682，1.682]，即

2)目标函数。以8 个细观参数(Ec，kn/ks，σc，cb，φb，gr，μc和Rmin)作为自变量，采用精度较高的宏−细观参数非线性关系式(式(5))计算宏观力学参数(σf，E，υ，σt，KIC和KIIC)，优化细观参数以实现模拟得到的宏观力学参数与实验值差距最小。为了保证目标函数的权重一致，需要对目标函数归一化处理：

式中：σ′f，E′，υ′，σ′t，K′IC和K′IIC为非线性关系式宏观值；σ*f，E*，υ*，σ*t，和为实验值。

3)约束条件。以8 个细观参数(Ec，kn/ks，σc，cb，φb，gr，μc和Rmin)作为自变量，采用宏−细观参数线性关系式(式(3))计算宏观力学参数(σf，E，υ，σt，KIC和KIIC)，模拟得到的宏观力学参数值与实验值相等(等式约束条件)。此外，根据岩石的抗压强度与抗拉强度之比一般在10～20之内[19]，由式(3)和式(5)计算得到的σ′f/σ′t/和σf/σt需满足不等式约束条件。等式约束条件为

不等式约束条件为

通过上述建立优化模型，采用遗传算法求得满足材料宏观力学参数的细观参数优化解(代码值)。

3 实验验证

基于上述建立的PFC3D细观参数优化模型(式(6)～(9))，通过遗传算法计算得到PFC3D 细观参数代码值，并根据等比例原则(参考式(4))转化为细观参数实际值。在此基础上，开展完整岩石压缩实验和裂隙岩石断裂试验的PFC3D 模拟计算，并将计算得到的宏观力学参数、应力−应变曲线和断裂轨迹与实验结果进行对比验证。

3.1 完整岩石单轴压缩实验

3.1.1 宏观力学参数

选取大理岩圆柱体试件(直径×高为50 mm×100 mm)的单轴压缩实验[17]开展PFC3D 模拟计算。大理岩的宏观力学参数为：σf=138.8 MPa，σt=8.391 MPa，E=41.5 GPa，υ=0.23，KIC=1.057 MPa·m1/2。II 型断裂韧度KIIC在文献[17]中未提供，KIIC可根据脆性岩石KIIC/KIC=2～4 来确定[32]，KIIC=2.114～4.228 MPa·m1/2。根据大理岩宏观力学参数和细观参数优化模型(式(6)～(9))，通过遗传算法可计算得到PFC3D 细观参数代码值解集，但这些解都集中在最优点附近。为了更好地与实验值匹配，本文选取细观参数(其表现的宏观性能参数σf，σt，E，υ，KIC和KIIC与实验值最接近)，并将其转化为实际值(见表11)进行数值模拟，得到的宏观力学参数如表12所示。从表12可见模拟计算值与实验值之间的相对误差均小于10%，精度较高，从而验证了基于PB+CCD 方法建立的考虑断裂韧度的PFC3D宏−细观关系式是合理的、有效的。

表11 大理岩PFC3D细观参数Table 11 PFC3D micro-parameters of marble

3.1.2 应力−应变曲线

图6所示为模拟计算得到的完整岩石单轴压缩应力−应变(σ−ε曲线)以及实验曲线。从图6可见：模拟的应力−应变曲线(O′B′C′)可分为弹性段(OB′)、损伤段(B′C′)和破坏段(C′D′)，但没有出现类似实验曲线(OABCD)的初始压密过程(OA)。这是因为在PFC3D 建模过程中，颗粒生成后紧密接触，不存在压密阶段。而模拟与实验相比，在压缩实验中，天然岩石内部存在的微孔和裂隙等缺陷会发生闭合而产生不可逆变形(OA段)。两者直线段斜率(弹性模量E)和应力最高点(抗压强度σf)十分接近(见表12)，不考虑压密阶段时，将实验曲线向原点水平左移后与模拟曲线基本重合，两者吻合度较高，进一步验证了基于PB+CCD方法建立的PFC3D宏−细观关系式的可靠性。

图6 完整大理岩试件单轴压缩应力−应变曲线Fig.6 Stress−strain curves of intact rocks under uniaxial compression

表12 大理岩宏观力学参数的PFC3D模拟值与实验值比较Table 12 Comparison of PFC3D simulation value and experimental value of of marble macro-mechanical parameters

3.1.3 破坏模式

图7所示为大理岩试件单轴压缩破坏的PFC3D数值模拟结果。为了更好地观察破坏模式，采用颗粒间的黏结状态表示试件破坏模式。图7中，灰色、褐色和黑色分别代表未破坏(黏结)、拉伸破坏和剪切破坏。从图7可见：试件以拉伸破坏(褐色)为主，剪切(黑色)较少，且破坏面大致沿着轴向呈60°～70°夹角，这与实验结果(图8)较吻合，再次验证了本文所建立的PFC3D 细观参数标定方法的有效性。

图7 大理岩试件单轴压缩数值模拟结果Fig.7 Numerical simulation results of marble specimen under uniaxial compression

图8 大理岩试件单轴压缩实验结果Fig.8 Experimental results of marble specimen under uniaxial compression

3.2 完整岩石三轴压缩实验

选取大理岩圆柱体试件(直径×高为50 mm×100 mm)三轴压缩实验[17]开展PFC3D模拟计算，其宏观力学参数和PFC3D 细观参数(见表11和表12)均与单轴压缩实验的相同。图9所示为大理岩在不同围压σ3下的三轴压缩实验结果与PFC3D 模拟计算结果(围压通过wall 伺服施加)。从图9可见：模拟得到的三轴压缩应力−应变曲线(图9(a))与单轴压缩应力−应变曲线类似，分为弹性、损伤和破坏3个阶段，均未出现实验曲线(见图6)的初始压密阶段，且与水平左移后的实验曲线基本重合。所不同的是，模拟曲线没有经过原点，是因为模拟试验是先施加围压，然后施加轴压，当轴压等于围压时才开始记录应力与应变，所以，模拟曲线的起点横坐标为0。此外，比较模拟曲线和实验曲线可知：随着围压增加，抗压强度增大，但弹性模量变化较小。

图9 大理岩三轴压缩实验结果与PFC3D模拟计算结果对比Fig.9 Comparsion of experimental results of marble under triaxial compression and PFC3D simulation

图9(b)所示为模拟计算得到的峰值强度σ1和围压σ3、基于库仑准则得到的C和φ及基于Hoek-Brown(HB)准则拟合得到的非线性强度包络线和相应实验结果，模拟结果和实验结果基本重合，模拟精度较高。

3.3 裂隙岩石断裂实验

3.3.1 宏观力学参数

为了进一步验证考虑断裂参数的PFC3D宏−细观参数非线性关系式的合理性，选取含平行双裂纹的大理岩试件单轴压缩断裂实验[32−33]开展PFC3D模拟计算。大理岩试件长×宽×高为152 mm×76 mm×32 mm，其中2条等长的平行斜裂纹长度为13 mm，宽度为1.3 mm，间距为13 mm，斜裂纹与水平面倾角为30°(见图11)。大理岩的宏观力学参数为：σf=88.3 MPa，σt=5.1 MPa，E=49 GPa，υ=0.19，KIC=0.65 MPa·m1/2。II 型断裂韧度KIIC在文献[33]中未提供，KIIC可根据脆性岩石KIIC/KIC=2～4 来确定[32]，KIIC=1.3～2.6 MPa·m1/2。

为了验证在标定PFC3D 细观参数时考虑断裂参数的必要性，根据大理岩宏观力学参数和细观参数优化模型(式(6)～式(9))以及遗传算法，同时计算考虑断裂韧度和不考虑断裂韧度的最优的PFC3D 细观参数代码值，分别为[−0.386 0，−0.206 0，−0.980 0，−0.899 0，−0.456 0，0.199 0，0.097 9，0.743 0]，[−0.298 0，0.047 0，−0.916 0，−0.906 0，0.331 0，−0.032 0，−0.732 0，0.738 0]。其中，考虑断裂韧度时，其细观参数通过综合考虑材料强度参数(σf，σt，E和υ)及断裂韧度(KIC和KIIC)(即将式(5)作为目标函数)求得；不考虑断裂韧度时，其细观参数通过只考虑材料强度参数(σf，σt，E和υ)(即只将式(5)的前4 项作为目标函数)求得。然后，将该细观参数值转化为细观参数实际值，如表13所示。

表13 大理岩PFC3D细观参数Table 13 PFC3D micro-parameters of marble

表14所示为采用PFC3D 细观参数实际值模拟计算得到的宏观力学参数。从表14可见：通过考虑和不考虑断裂韧度确定的细观参数模拟的宏观力学参数与实验值之间的相对误差均小于10%，KIIC的相对误差也小于10%，表明该两组细观参数均合理、有效。

表14 大理岩力学参数的模拟计算值与实验值比较Table 14 Comparison of simulated and experimental values of mechanical parameters in marble

3.3.2 裂纹扩展过程

图10所示为考虑断裂韧度的大理岩双裂纹试件单轴压缩裂纹扩展过程的PFC3D数值模拟结果。从图10可见：双裂纹的外尖端B′和C′及裂纹中部H′和I′点发生拉伸起裂(图10(b)和图10(f))；随着进一步加载，从裂纹尖端B′和C′点起裂的裂纹沿着加载方向扩展(尖端B′向上延伸，尖端C′向下延伸，与原裂纹呈约30°的方向扩展)，受到试件两端压密阻碍；同时，B′D′点之间产生1条次生裂纹(处于萌生阶段)。由于2 条预制裂纹的干涉影响，从H′点和I′点起裂的裂纹沿着垂直于原裂纹的方向扩展(图10(c)和图10(g))；当荷载不断增加时，试件除预制裂纹以外的位置也产生拉伸破坏，但预制裂纹附件拉伸破坏较集中(图10(d))。最终，B′D′间产生的次生裂纹逐渐贯通，从H′和I′点起裂的2条裂纹沿着垂直原裂纹的方向，分别在双裂纹的内尖端D′和A′点贯通，从B′和C′点起裂的2 条裂纹(拉伸裂纹为主)分别沿轴向贯通至上、下表面E′和F′点(图10(e))。

图10 大理岩双裂纹试件单轴压缩裂纹扩展过程模拟结果(考虑断裂韧度)Fig.10 Simulation results of crack propagation trajectory of double fractured marble under uniaxial compression(considering fracture toughness)

图11所示为大理岩双裂纹试件单轴压缩实验结果[34]，其中T表示拉伸裂纹，箭头表示裂纹扩展方向。虽然大理岩试件中的初始裂缝与加载方向呈60°角，预制裂纹面上承受压−剪复合应力作用，但在裂纹尖端附近的任意其他方向上仍受到拉应力的作用。从数值模拟结果(图10(b)和图10(c))也可看出：2条预制裂纹从外尖端起裂后，沿着与原裂纹呈约30°的方向扩展，且破坏模式以拉伸破坏为主(褐色)，因此，破坏机理为I 型断裂；在单轴压缩作用下，双裂纹扩展过程是从双裂纹的外尖端B和C点与裂纹中部的H和I点同时起裂，随后，从B和C点起裂后的裂纹分别沿着轴向加载方向向上、向下扩展，同时产生1 条次生裂纹BG沿着轴向加载方向向下扩展(图11(b)和图11(c))；从H和I点起裂的裂纹逐渐向双裂纹内部扩展，直到与双裂纹的内尖端D和A点贯通，最终的破坏轨迹为2条沿轴向加载方向(BE和CF)以及2 条垂直于原裂纹方向(HD和IA)的宏观主裂纹以及衍生的次生裂纹BG(图11(b))。可见，除了次生裂纹略有细微差别外(实际岩石材料的非均质性所致)，实验结果与PFC3D 模拟计算结果较吻合，从而验证了综合考虑强度和断裂参数的PFC3D宏−细观参数标定方法的有效性。

图12所示为不考虑断裂韧度的大理岩双裂纹试件单轴压缩裂纹扩展过程的PFC3D 数值模拟结果。将图12所示结果与图10考虑断裂韧度的结果及图11所示实验结果进行对比，验证考虑断裂韧度的必要性。从图11和图12可见：在加载初期，裂纹也是沿着预制裂纹外尖端B″和C″点以及预制裂纹中部的H″和I″点拉伸起裂(图12(b)和12(f))；随着进一步加载，从裂纹尖端B″和C″点起裂的裂纹，基本沿着轴向扩展，同时B″D″点之间产生1条次生裂纹。但从H″和I″点始起裂的裂纹沿着轴向扩展(图12(c)和图12(g)并在中部逐渐贯通，这与考虑断裂韧度的模拟结果和实验结果不一致；随着荷载不断增加，试件除预制裂纹以外的位置也逐渐产生大量拉伸破坏，伴随预制裂纹尖端少量剪切破坏(图12(d))。最终B″D″点之间的次生裂纹贯通，H″和I″点始起裂的裂纹沿着轴向在中部贯通，从B″和C″点起裂的裂纹(拉伸裂纹为主)分别贯通至上、下表面E″和F″点(图12(e))。对比图12(e)和图10(e)，考虑断裂韧度的裂纹扩展轨迹B′E′和C′F′要比不考虑断裂韧度的裂纹扩展轨迹B″E″和C″F″更聚集。

图11 大理岩双裂纹试件单轴压缩实验结果Fig.11 Experimental results of Double Crack marble specimen under uniaxial compression

图12 大理岩双裂隙单轴压缩裂纹扩展轨迹模拟结果(不考虑断裂韧度)Fig.12 Simulation results of crack propagation trajectory of double fractured marble under uniaxial compression(without considering fracture toughness)

综上可知，考虑和不考虑断裂韧度的双裂隙大理岩单轴压缩的扩展轨迹在2条预制裂纹起裂后之间的扩展模型存在明显差异。当不考虑断裂韧度时，2条原裂纹内尖端附近起裂后直接贯通，只形成1 条分支裂纹；而考虑断裂韧度时，2 条原裂纹内尖端附近起裂后，各自扩展形成2 条分支裂纹，这与实验结果相吻合。因为大理岩双裂纹试件的破坏机理为I型断裂，虽然表13中考虑和不考虑断裂韧度的大理岩PFC3D 细观参数差异不大，但考虑了断裂韧度时，计算得到的I型断裂韧度KIC和抗拉强度σt更接近于实验值，因此，考虑断裂韧度的裂纹扩展轨迹模拟结果要比不考虑断裂韧度的模拟结果与实验结果更吻合。

4 结论

1)基于PB+CCD 试验设计方法和遗传算法，提出了综合考虑强度和断裂参数(I型和II型断裂韧度)的PFC3D 宏−细观参数标定新方法，建立了PFC3D平直节理模型的宏−细观参数表达式，并利用遗传算法求得了PFC3D 细观参数的优化解。该方法可推广至其他接触模型(接触黏结模型、平行黏结模型)，为PFC3D定量分析工程材料多裂纹扩展的细观机理提供参考。

2)基于PB设计计算得到的PFC3D细观参数对宏观参数的影响规律为：平直节理黏聚力cb对单轴抗压强度σf和II 型断裂韧度KIIC的影响最大，平直节理抗拉强度σc对宏观抗拉强度σt和I 型断裂韧度KIC的影响最大，平直节理弹性模量Ec和刚度比kn/ks分别对弹性模量E和泊松比υ的影响最大，且cb与σf和KIIC呈正相关关系，σc与σt和KIC呈正相关关系，Ec与E呈正相关关系，kn/ks与υ呈正相关关系。

3)基于CCD方法计算得到的PFC3D细观参数交互作用对宏观参数的影响规律为：平直节理黏聚力cb和安装间隙比gr的交互作用对单轴抗压强度σf有显著影响，且呈正相关；平直节理刚度比kn/ks和抗拉强度σc的交互作用对泊松比υ有显著影响，且呈负相关；Ec与kn/ks之间、Ec与gr之间存在交互作用，Ec·(kn/ks)与E呈负相关关系，Ec·gr与E呈正相关关系。

4)基于PFC3D宏−细观参数表达式定量确定的细观参数，模拟计算了完整岩石的单轴和三轴压缩应力−应变曲线、宏观力学参数和破坏模式，与试验结果较吻合。

5)通过对比PFC3D 细观参数标定时考虑和不考虑断裂韧度下模拟计算得到的大理岩双裂纹试件单轴压缩断裂轨迹，考虑断裂韧度的试件断裂轨迹与试验结果更吻合。