1. 引言
岩石是一种我们在日常生活中经常接触到的天然的多孔材料，其内部包含着大量不规则、跨尺度的孔隙，这些孔隙直接影响着岩石的宏观物理、力学和化学性质 [1] [2] [3] 。由于各种岩石的组分和结构各异，形成的年代不同，其中还有许多裂隙，致使其力学性质相差很大。在地质学中，无论从已知外力和边界条件去推测地块变形，或从地块的变形去反推地块所受的边界条件和外力，都需了解当时当地的岩石力学性质。探讨研究孔隙结构与岩石宏观物理力学性质之间的内在关系，对于解决石油、地质、采矿、冶金、土木和水利工程中的实际问题具有十分重要的意义。
由于岩石中天然孔隙形态复杂、分布无序、数量庞大，难以定量而准确地描述其结构特征，大量研究表明，岩石孔隙分布具有显著的统计分形特性，分形理论是刻画岩石孔隙复杂性和不规则性的有效方法。马新仿等 [4] [5] 研究了低渗透岩石的孔隙结构，建立了毛管压力和孔隙大小概率密度分布的分形几何模型，计算了孔隙结构的分维数和孔径概率密度分布；J. Lai和G. W. Wang [6] 利用高压压汞试验测试了Bashijiqike致密含气砂岩的孔径分布特征，发现砂岩孔隙系统主要由晶界大孔和晶内微孔组成，孔隙分形维数可用来评价微孔孔隙结构和储层非均匀性；A. Giri等 [7] 针对岩石孔隙度相同而孔径分布不同的特点，利用RBBD模型重构不同孔隙率的沉积岩三维结构，并进行孔隙分布的多重分形分析，基于此提出了一种控制岩石输运性质的结构参数；Z. Y. Zhang和A. Weller [8] 用毛管压力曲线和核磁共振研究了24个岩样的孔隙结构，通过3种方法计算了岩石孔隙表面分维和体积分维，并用于渗透率预测。这些成果极大深化了人们对岩石内孔隙孔径分布规律的认识。
存在于大自然的岩石大多都有着或多或少的孔洞，孔洞的大小及分布规律也有很大的区别，这些孔洞的存在会改变岩石弹性模量等力学性质，弱化岩体结构的稳定性，对岩石有着巨大的影响，在这之前已经有大量专家人员对此做了许多相关的研究，也有了很多成果。但是在研究不同分形维数情况下孔洞孔径分布对岩石力学性质等的研究并不太多见。本文在阅读各类文献和资料基础上，通过RFPA软件建立一个二维的长方形的岩石模型，然后通过随机数生成器进行随机数的产生，用生成的随机数来表示我们要在岩石模型中挖去的孔隙的位置，最后在单轴压缩下对岩石的以下内容进行研究和探讨：
1) 孔隙率对岩石应力–应变曲线、抗压强度、弹性模量、破坏方式的影响。
2) 孔隙孔径分布分形维数对岩石应力–应变曲线、抗压强度、弹性模量、破坏方式的影响。
2. 数值试验概述
2.1. 材料参数
选择采用花岗岩试件材料属性其材料属性见表1。
2.2. 孔洞半径
根据孔隙率和分形维数的不同总计需要建立9个模型，各个模型下各圆形小孔半径见表2~表4。孔隙率分成三种：1) 孔隙面积占整个模型的2%，即孔隙总面积为100mm × 50mm × 2% = 100 mm2。2) 岩石孔隙面积占整个模型的5%，即孔隙总面积为100 mm × 50 mm × 5% = 250 mm2。3) 岩石孔隙面积占整个数值模型的10%，即孔隙总面积为100 mm × 50 mm × 10% = 500 mm2。根据资料花岗岩的孔隙孔径分
Table 1. System resulting data of standard experiment
表1. 标准试验系统结果数据
Table 2. The radius of each circular hole when the fractal dimension is 2
表2. 分形维数为2时各圆形小孔半径
Table 3. The radius of each circular hole when the fractal dimension is 2.5
表3. 分形维数为2.5时各圆形小孔半径
布分形维数在本次数值试验中分别选取2、2.5和3。所有的孔隙都是由各种半径不同的圆形小孔来表示，在这次试验中选取五种不同大小半径的圆形小孔。
式中，D是模型孔洞孔径分布的分形维数，r表示圆形小孔的半径，N表示相对应的半径大小的圆形小孔的个数。
2.3. 数值试验方法
运用RFPA软件对所研究岩石的破坏过程进行数值模拟 [9] [10] ，建立二维长方形模型：宽50 mm，高100 mm，划分为250 × 500总计125,000个网格数，使用Linux随机数生成器来确定圆形小孔在模型上的具体位置，建立了含有孔隙的岩石的二维模型，具体见图1~图9。采用位移加载方法来施加载荷，方式是单轴压缩，每一步的载荷增量为0.005 mm，求解时使用直接法，西塔值为1，时间步长为1，最大迭代次数为300，迭代精度为1e−010。
Figure 1. The 2D model of the pore when the porosity is 2% and the fractal dimension is 2
图1. 孔隙率为2%，分形维数为2时孔隙的二维模型
Table 4. The radius of each circular hole when the fractal dimension is 3
表4. 分形维数为3时各圆形小孔半径
Figure 2. The 2D model of the pore when the porosity is 2% and the fractal dimension is 2.5
图2. 孔隙率为2%，分形维数为2.5时孔隙的二维模型
Figure 3. The 2D model of the pore when the porosity is 2% and the fractal dimension is 3
图3. 孔隙率为2%，分形维数为3时孔隙的二维模型
3. 数值试验结果分析
3.1. 破坏过程分析
通过对不同孔隙率和不同分形维数的岩石试件破坏过程的分析，可以得到以下规律：
Figure 4. The 2D model of the pore when the porosity is 5% and the fractal dimension is 2
图4. 孔隙率为5%，分形维数为2时孔隙的二维模型
Figure 5. The 2D model of the pore when the porosity is 5% and the fractal dimension is 2.5
图5. 孔隙率为5%，分形维数为2.5时孔隙的二维模型
1) 在相同分形维数情况下，岩石试样随着孔隙率的增加破坏的速度越来越快，即孔隙率越大完全破坏时所需加载的步数越少。
2) 几乎所有模型都是先贯通最大的两个孔洞，并最终大致沿着两个孔洞圆心连成的直线方向与上下边界贯通，只有在孔隙率 = 2，分形维数D = 3的时候破环方式不同，可能是由于孔隙率太小并且分形维数有点大导致各个圆孔的大小差别不是太明显引起的原因。
Figure 6. The 2D model of the pore when the porosity is 5% and the fractal dimension is 3
图6. 孔隙率为5%，分形维数为3时孔隙的二维模型
Figure 7. The 2D model of the pore when the porosity is 10% and the fractal dimension is 2
图7. 孔隙率为10%，分形维数为2时孔隙的二维模型
3) 所有模型的破坏方式都是拉伸破坏和剪切破坏相结合。
4) 在同一孔隙率下，分形维数D越大在完全破坏时明显的裂纹越多。
Figure 8. The 2D model of the pore when the porosity is 10% and the fractal dimension is 2.5
图8. 孔隙率为10%，分形维数为2.5时孔隙的二维模型
Figure 9. The 2D model of the pore when the porosity is 10% and the fractal dimension is 3
图9. 孔隙率为10%，分形维数为3时孔隙的二维模型
3.2. 孔隙率对岩石应力–应变曲线的影响分析
在相同分形维数情况下不同孔隙率岩石试件模型的应力–应变曲线比较图见图10~图12。
Figure 10. Stress-strain diagram when fractal dimension D = 2
图10. 分形维数D = 2时应力–应变图
Figure 11. Stress-strain diagram when fractal dimension D = 2.5
图11. 分形维数D = 2.5时应力–应变图
可得结论：
1) 在同一种分形维数情况下，峰值应力随着孔隙率的增大而减小，孔隙率为2%时峰值应力最大，孔隙率为10%时最小，且差距比较大。
2) 上面三组岩石模型在达到峰值应力前的大部分时间都呈线性变化，处于弹性变形阶段；且应力增长的速度随着孔隙率的变化而改变，孔隙率越小应力增长速度越快。
3) 在同一种分形维数情况下，孔隙率越大达到峰值应力所需的加载步数越少，在本次实验中，孔隙率为10%时所需的步数最少，孔隙率为2%时所需的步数最多。
3.3. 孔隙分形维数对应力–应变曲线的影响分析
在相同孔隙率下不同分形维数岩石试件模型的应力–应变曲线比较图见图13~图15。
可得结论：
1) 在同一孔隙率情况下，峰值应力随着分形维数D的增加而增加，在分形维数D为2时峰值应力最小，在分形维数D为3时最大，但是差距不是太大。
2) 上面三组岩石模型在达到峰值应力前的大部分时间都呈线性变化，处于弹性变形阶段；且通过对曲线斜率的计算知道应力增加的速度随着分形维数D的增加都加快，但变化不明显。
Figure 12. Stress-strain diagram when fractal dimension D = 3
图12. 分形维数D = 3时应力–应变图
Figure 13. Stress-strain diagram when porosity = 2%
图13. 孔隙率 = 2%时应力–应变图
Figure 14. Stress-strain diagram when porosity = 5%
图14. 孔隙率 = 5%时应力–应变图
Figure 15. Stress-strain diagram when porosity = 10%
图15. 孔隙率 = 10%时应力–应变图
3.4. 孔隙率和孔隙分形维数对抗压强度的影响分析
根据数值模拟实验结果得到不同孔隙率和不同分形维数的模型的抗压强度数据，进而得到孔隙率与抗压强度的关系图和分形维数与抗压强度的关系图，见图16，图17。
可得结论：
1) 在同一分形维数情况下，抗压强度的大小随着孔隙率的变化而发生变化，孔隙率越小抗压强度越大，在本次实验中，孔隙率为2%时抗压强度最大，孔隙率为10%时抗压强度最小，且相互之间差距比较大。
2) 在相同孔隙率的情况下，抗压强度的大小随着分形维数的变化也发生变化，分形维数D越大时抗压强度越大，在这次实验中，分形维数D等于3时抗压强度最大，分形维数D等于2时抗压强度最小，
Figure 16. The relationship between porosity and compressive strength
图16. 孔隙率与抗压强度的关系
Figure 17. The relationship between fractal dimension and compressive strength
图17. 分形维数与抗压强度的关系
但彼此之间差距不明显。
3) 抗压强度的差值随着孔隙率间差距的变大而变小。
4) 在孔隙率为2%时，分形维数D等于2.5和分形维数D等于3情况下抗压强度特别接近，在孔隙率为10%时，分形维数D等于2和分形维数D等于2.5情况下抗压强度接近，在孔隙率等于5%时抗压强度在不同分形维数情况下有一定的差距。
3.5. 孔隙率和孔隙分形维数对弹性模量的影响
根据之前的不同情况下的应力应变曲线计算得到了不同孔隙率和不同分形维数时岩石模型的弹性模量E，进而得到孔隙率与弹性模量的关系图和分形维数与弹性模量的关系图，见图18，图19。
Figure 18. The relationship between porosity and elastic modulus
图18. 孔隙率与弹性模量的关系
Figure 19. The relationship between fractal dimension and elastic modulus
图19. 分形维数与弹性模量的关系
可得结论：
1) 在同一分形维数情况下，弹性模量的大小随着孔隙率的变化而发生变化，孔隙率越小弹性模量越大，在本次实验中，孔隙率为2%时弹性模量最大，孔隙率为10%时弹性模量最小，且相互之间差距比较大。
2) 在相同孔隙率的情况下，弹性模量的大小随着分形维数的变化也发生变化，分形维数D越大时弹性模量越大，在这次实验中，分形维数D等于3时弹性模量最大，分形维数D等于2时抗压强度最小，但相互之间差距不太明显。
4. 结论与展望
1) 在同一分形维数情况下，岩石试样随着孔隙率的增加破坏的速度越来越快,即孔隙率越大完全破坏时所需加载的步数越少；在弹性阶段应力增长的速度随着孔隙率的减小而变快，但是达到峰值应力所需的加载步数却随着孔隙率的减小而增多；孔隙率越小，岩石的抗压强度和弹性模量就越大，且差距比较明显。
2) 在同一孔隙率情况下，分形维数D越大岩石模型在完全破坏时明显的裂纹就越多；在弹性阶段应力增加的速度随着分形维数D的增加都变快，但不明显；分形维数D越大，岩石模型的抗压强度和弹性模量也越大，但差距较小。
3) 本文中全部选用摩尔库伦破坏准则，而没有将其与其它的Hoek-Brown准则、Double-shear准则和Single-shear准则进行对比，探讨各种准则的适用范围。
基金项目
江苏省大学生创新训练计划(201710290019Y)；国家自然科学基金(51604260)；中国博士后科学基金特别资助项目(2016T90519)；江苏省博士后科研资助计划项目(1501085B)。
NOTES
*通讯作者。