单轴压缩条件下岩石细观力学参数特性数值研究

艾 涛 朱珍德 田 源 朱 端

(1. 河海大学 岩土力学与堤坝工程教育部重点实验室， 南京 210098; 2. 河海大学 江苏省岩土工程技术工程研究所， 南京 210098)

0 引 言

工程岩体内部蕴含着大量的节理、裂隙，在外部荷载的扰动下岩体内裂隙将会逐渐扩展及贯通，从而引起整个岩体的破坏．因此，在细观层次对岩石进行研究可以很好地揭示岩石在外部荷载作用下细观损伤断裂过程的破坏机理．计算机技术的发展和数值模拟方法的成熟，给复杂岩石结构变形与破坏细观分析提供了有效的方法，而它计算速度快、精度高等优点也使得其在解决土木工程问题上得到了广泛的应用．

T.X.Bai等[1]利用有限元软件FRANC对裂隙岩体中影响裂纹扩展规律的几个重要因素如弹性模量、预制裂纹深度等进行了数值研究．朱珍德等[2]利用扫描电镜在岩石的细观结构参数量化方面进行了大量的研究并编制了参数提取的程序．唐春安等开发了岩石破坏过程分析系统RFPA(Rock Failure Process Analysis)[3]，在此系统基础上进行了岩石三轴破裂过程数值模拟[4]．魏玉寒等[5]基于扩展有限元法研究了三维含裂隙岩体中裂隙倾角对于表面裂纹扩展的影响．另外还有一些典型的数值模型，如格构(Lattice)模型[6]、颗粒(Granular)模型[7]、梁一颗粒模型[8]等等，这些模型引入介质的非均匀性，对研究非均匀性带来的材料破裂过程的复杂性颇有帮助．现存的岩石受荷破坏理论中，传统的宏观经典弹塑性理论力学模型忽略了岩石细观结构的非均匀因素，而细观损伤统计模型虽然考虑了非均匀的影响，但依然不足以表达整个破坏过程的细观损伤复杂性．要深入研究岩石的破坏机理需要将宏、细观损伤破坏问题结合起来，建立多尺度的力学模型．在非均匀性研究的基础上，通过考虑宏、细观结合的数值分析方法将有助于更深入地了解岩石受荷破坏机理．

本文考虑岩石介质的非均匀性，结合细观统计损伤模型的思想，采用Weibull分布来描述离散后岩石细观单元的力学性质参数，并应用Monte-Carlo方法编制程序进行单元力学参数的赋值．在此基础上，通过ABAQUS有限元软件的二次开发实现对单轴压缩荷载作用下岩石细观单元的力学参数的影响的模拟研究．

1 有限元法数值模型

1.1 ABAQUS软件平台以及用户子程序接口

ABAQUS作为当前国际上最为通用的有限元计算分析软件之一[9]，其非线性力学分析功能具有独特的优势，从而在土木、水利、机械等领域得到了广泛应用．尽管如此，国内岩土工程领域比较通用的一些本构模型，现有的ABAQUS版本却并不具备，这无疑限制了ABAQUS在国内岩土工程领域的应用与推广．但是ABAQUS提供了二次开发的接口，通过其强大的分析求解平台[10]，可使困难的分析简单化，使复杂的过程层次化．

ABAQUS为了方便用户针对自己研究领域开发本构模型，采用了FORTRAN语言接口，提供了若干用户子程序(User Subroutines)和在编程时可以调用的实用工具(Utility Routines)[11-12]．常用的接口程序有：用于定义不同物理场的自定义场变量子程序USDFLD、用于定义材料力学行为的用户材料子程序UMAT、用于定义不同单元模式的用户单元子程序UEL等，对于材料、单元、算法等各方面均有涉及．对于细观单元模型，即可利用用户子程序USDFLD将表征单元损伤程度值定义为场变量，而场变量的值依赖于积分点上应力、应变状态，在分析步的每一增量步中每一积分点都将调用USDFLD子程序进行不断的迭代更新．

1.2 单元力学性质的赋值

在前人的研究基础上，采用Weibull分布来描述离散后岩石细观单元的力学性质参数，其分布密度函数为：

(1)

式中，u为细观单元体的力学性质参数(强度、弹性模量等)；u0为与细观单元力学性质的平均值有关的参数，其数值并非该平均值；m为分布函数的形状参数，其物理意义反映了岩石介质的均匀性，即均质度．

基于Monte-Carlo方法和统计描述相结合的方法，通过编制程序实现对单元力学性质(弹性模量和强度)的初始化赋值．以弹性模量的赋值进行说明，设所有单元的弹性模量的平均值为E0，则基于式(1)弹性模量weibull分布函数的积分为：

(2)

其反函数为：

(3)

1.3 数值模型的建立

数值试样采用尺寸为50 mm×100 mm二维模型，如图1所示，试样加载方式为以刚性体控制的位移加载模式，在位移加载的分析步中，采用100个固定的增量步．数值模型采用CPS单元(平面应力单元)，考虑到计算精度和计算资源的协调问题，进一步选用CPS4R，即细网格划分的二维平面线性减缩积分单元进行计算，其常用于大应变分析模拟．

图1 岩石单轴压缩数值模型

建模过程中先运用ABAQUS前处理软件对岩石细观损伤数值模型的几何尺寸、单元划分、边界条件等建立初始模型；然后根据数值模型文件(INP文件)的特点，通过编制程序，生成分别包括单元的信息、单元集合的信息、材料参数的信息的文本文件，并用参数嵌入INP文件中；接着利用Fortran语言编制USDFLD子程序和相关程序进行各细观单元力学性质的赋值、修改，通过调用Utility Routines中GETVRM函数，就可以利用强度准则进行判断，对单元的力学参数进行折减以进行下一增量步的计算；最后保存模型以待运行分析．

2 数值模型的力学参数研究

数值模型中涉及到单元的大小以及细观单元的力学参数包括：弹性模量和强度的均值E0、fc 0和均质度m，压拉比χ，残余强度系数λ．对于细观单元的力学参数，岩石破坏细观试验研究提供了破坏的细观机理，但由于目前岩石细观结构力学性质的试验研究还很少，细观单元的力学参数还很难通过细观试验获取，只有通过大量的数值模型研究来对细观力学参数进行标定．如没有特别说明，数值模拟中材料泊松比μ和摩擦角θ分别采用0.25和30°．设置不同单元数目进行数值破坏形态模拟，图2为采用1 mm×1 mm的单元尺寸(5 000个单元)下的数值试样破坏形态图，从破坏形态上看，该单元划分可以满足数值模拟要求．

图2 5 000个单元数目下数值试样的破坏形态

2.1 均质度m对宏观的影响

为得出m对宏观应力应变曲线的影响，建立均质度分别为1.5、3、6的模型，具体参数见表1，结果如图3所示．由图3可见，数值试样的宏观响应对于细观单元整体均质度的变化很敏感，随着m的增大，试样的脆性明显随之增大，宏观的弹性模量和强度也相应提高．m越大，数值试样越均匀，宏观力学参数越接近细观参数，随着m的增大，曲线的峰值应变也随之增大．

表1 不同单元尺寸的数值试样

图3 不同m下的应力应变曲线

2.2 宏细观弹性模量和强度的关系

为了得出一定均质度下，数值试样的细观单元的E0，fc 0与宏观弹性模量和宏观强度的关系，建立的数值试样以及宏观弹性模量和强度的结果见表2．

表2 不同细观弹性模量和强度下数值试样及其结果

对m=1.5和3下的数值试验结果与细观单元的强度和弹性模量的拟合结果如图4～5所示，可见细观单元强度、弹性模量和宏观强度、弹性模量的拟合具有很好的线性关系，可以得到：

(4)

(5)

式(4)(5)中，fc 0和E0为Weibull分布下细观单元强度和弹性模量的均值，fc和E为数值试样宏观的强度和弹性模量．

图4 细观单元强度与宏观强度的拟合曲线

图5 细观单元弹性模量和宏观弹性模量的拟合曲线

2.3 压拉比以及残余强度系数对宏观的影响

1)压拉比的影响

岩石细观单元的的压拉比χ反映了岩石的抗压和抗拉能力的差异，不同压拉比系数下数值试样的应力应变曲线如图6所示．所有试样的应力应变曲线都经历了线性变形、非线性变形、不断劣化和残余变形4个阶段．随着χ的降低，数值试样的峰值强度逐渐提高，所有曲线都显示相同线性段，残余强度在χ>5时，随着χ的降低有所提高但不明显，χ=3时残余强度提高很明显．对于相同的单元压缩强度fc 0=300 MPa，不同的压拉比对应不同的单元受拉强度．

图6 不同压拉比数值试样应力应变曲线

对于本文中摩擦角为30°，当χ=3时，如果单元发生受拉破坏，由于σ3≤ft=-fc 0/3⟹-3σ3≥fc 0⟹σ1-3σ3≥fc 0，可见，当χ≤3时，单元如果满足最大拉应力准则则必然满足摩尔库伦准则，即单元在受剪破坏一定发生在受拉破坏之前．之所以随着压拉比系数的降低，数值试样的宏观强度逐渐增加，是因为压拉比越小，细观单元的剪切破坏越占主导，受拉破坏的单元数目越少．

2)残余强度系数的影响

为研究细观单元的残余强度系数λ对宏观应力力学性质的影响，分别建立残余强度系数为0.1，0.3，0.5的数值试件．当压拉比为5时，3个数值试样的应力应变曲线如图7所示．

图7 不同细观残余强度下的应力应变曲线

由图可见，随着细观单元残余强度系数的提高，试样的峰值强度有一定的增加，试样的延性有一定的增强，数值试样宏观的破坏后残余强度并没有增加，几乎保持不变，略有降低．这主要是因为残余强度系数只是对受剪破坏的单元起作用，受拉破坏的单元几乎没有承载力，所以压拉比越小，单元强度残余系数的影响越大．虽然在峰值应力附近，承载力有一定的提高，但如果受剪破坏的单元受力状态最终变为拉伸破坏时，单元将不再具有承载力．因此，试样最后宏观残余强度都很小，变化不大，甚至随着λ的提高，试样最终破坏后的宏观残余强度略有降低．

3 岩石单轴压缩过程的数值模拟

通过对数值试验中细观单元参数的研究可知，数值试样细观参数的选取取决于宏、细观的试验．为研究岩石在单轴荷载下的破坏过程，结合武沂泉[13]的细观压缩SEM试验及其宏观的应力应变曲线建立对应的数值试样．

武沂泉[13]利用采自四川锦屏二级水电站工程现场的大理岩作为试样进行单轴加载．使用SEM对试件采用分层式扫描拍摄以采集试件在各级荷载下的SEM图像试件在各级荷载下SEM图像，并通过图像处理提取了裂纹信息，得到了细观裂纹随应力加载的变化规律及应力应变曲线．由宏观应力应变曲线得出的宏观的峰值强度为77.46 MPa，峰值应变为0.003 4，弹性模量为30 571 MPa，根据公式(4)或(5)可得到细观单元的强度fc 0=902.8 MPa和弹性模量E0=50 818.7 MPa．综合考虑其它参数的影响并通过试算，最终单元的细观参数情况确定见表3．

表3 压缩数值试样的参数表

不同应力状态下单轴压缩试验的SEM图片及数值试样破坏情况的试样破坏过程及形态如图8所示．可见，当加载到应力20 MPa左右时，试件表面变化不是很明显，由于局部初始损伤和细观不均匀性的存在，只有少量的局部黑点处出现拉应力集中，对应SEM图片中的白色发亮处，这些微裂纹方向均与轴向荷载方向一致或成小夹角，且在宽度上都很细．继续加载达到60 MPa，此时由扫描图片可以看出微裂纹数目和宽度都有所增加，裂纹扩展方向依然与轴向荷载几乎一致，0～60 MPa这一过程右面数值试样的破裂情况处于微破裂均匀随机发展阶段，此阶段细观单元的破坏是随机出现并比较均匀分布在试样中．当应力超过60 MPa，数值试样的破裂情况进入微破裂非均匀发展阶段，单元的破坏由随机均匀的破坏转为在某些单元附近发生贯通，并逐渐向集中的某几个局部区域发展，破坏单元的数目持续增多，在对应的左边扫描图像中裂纹的宽度急剧增大，裂纹进入不稳定扩展阶段(图8(c))．应力应变曲线的峰值后下降段对应数值试样的局部及宏观破裂发展阶段，微破裂的局部化趋于显著，单元的破坏数目迅速增多，破裂单元迅速贯通，当达到峰值时，数值试样完全破坏，形成破坏的剪切带(图8(d))．

图8 不同应力状态下SEM图片及数值试样破坏情况

由以上岩石单轴压缩的真实试验与数值模拟的结果对照可见，两者在应力应变曲线以及破坏形态和机理上都存在一定的相关性，数值模拟的结果在一定程度上可以合理地反映岩石单轴荷载作用下应力应变曲线以及从细观破坏直至宏观破坏的整个过程．

单轴压缩下数值试验得到的应力应变全曲线以及对应的SEM试验试样曲线如图9所示，可见数值试样和试验试样的应力应变曲线大体上是很接近的，说明数值试样的建立是合理的，也验证了细观参数对宏观试样的影响．从两者的曲线对比来看，数值试样的应力应变曲线反映了单轴压缩荷载作用下岩石的非线性变形阶段、弱化、以及残余强度的特征．

图9 数值试样与试验试样应力应变全曲线

4 结 论

本文基于岩石介质的非均匀性理论以及细观力学的基本思想，对细观单元各主要参数的影响以及荷载作用下的破裂过程和机理进行了分析．主要得出以下结论：

1)随着单元均质度的增大，试样的脆性增加，数目越多，越能反映试样的破坏形态．

2)细观单元和宏观试样的弹性模量、强度之间存在很好的线性关系．

3)随着细观单元的压拉比的降低、数值试样的宏观强度呈现增大趋势，而细观单元的残余强度系数对数值试样的宏观的残余强度影响甚微，随着其提高，对应力应变曲线硬化段有一定的提高．

4)岩石单轴压缩数值模型结果表明数值试样和试验试样的应力应变曲线和破坏情况大体上是很接近的，说明数值试样的建立是合理的．