岩石亚临界裂纹扩展应力强度因子Williams级数计算及数值模拟

郝瑞卿，李江腾，曹平，廖峻，刘博

(中南大学 资源与安全工程学院， 湖南 长沙，410083)

摘要：引入 Williams级数，提出一种计算双扭试验中应力强度因子的新方法，并运用ANSYS数值模拟软件，对双扭试验进行数值模拟。通过理论分析和数值模拟2种方法，分别计算大理岩和花岗岩在亚临界裂纹扩展中的应力强度因子和。同时，将理论分析和数值模拟所得的应力强度因子和与试验所得的应力强度因子进行比较。研究结果表明：理论分析、数值模拟与试验这3种方法所得的应力强度因子均随裂纹的扩展而减小，与断裂力学理论相符；对于同一种岩石试样，理论分析所得的与试验所得的的拟合曲线相比斜率偏小，但偏差较小；模拟所得的与试验所得的的拟合曲线斜率基本一致，说明理论分析、数值模拟与试验过程中应力强度因子减小规律相同，能反映试验的基本过程，可见理论分析与数值模拟结果可靠，能为岩土工程稳定与时间的相关性研究提供依据。

关键词：应力强度因子；双扭试验；Williams级数；亚临界裂纹扩展；数值模拟

中图分类号：TU318                 文献标志码：A         文章编号：1672-7207(2014)11-3949-08

Williams series calculation and numerical simulation for stress intensity factor of rock subcritical crack growth

HAO Ruiqing, LI Jiangteng, CAO Ping, LIAO Jun, LIU Bo

(School of Resources and Safety Engineering, Central South University, Changsha 410083, China)

Abstract: A new method to calculate the stress intensity factor was proposed by using Williams series and some numerical simulations were made with ANSYS software. Though theory analysis and simulation calculation, the stress intensity factors were obtained during the subcritical crack growth for marble and granite. And the stress intensity factor obtained by theory analysis, obtained by numerical simulation were compared with the stress intensity factor obtained by double-torsion test. The results show that the stress intensity factors obtained from three methods decrease while the crack propagates. And to the same rock, the slope of fitting curves for is lower than that of the fitting curves for , but the deviation is small. and have almost the same decline rate, and the slope of fitting curves is basically identical, which shows that the reduce rule of stress intensity factor is identical to those of the theory analysis, simulation calculation and test, and so the results of theory analysis and simulation calculation are reliable, which can provide the basis for studying time-dependency of rock stability.

Key words: stress intensity factor; double-torsion test; Williams series; numerical simulation; subcritical crack growth

近年来，人类修建的岩土工程规模不断扩大，为保证这些工程的顺利进行以及工程结束后使用的稳定性、安全性，要求岩体工程能保持长期的稳定性。根据断裂力学原理，在岩体裂纹尖端的应力强度因子K_Ⅰ未达到其断裂韧度时，裂纹会以一种较稳定的速度向前扩展，这种岩石裂纹稳定的、准静态扩展方式称为岩石的亚临界扩展。而当亚临界扩展进行到一定程度时，通常岩石内裂纹将会快速扩展，进而发生断裂，导致岩土工程失稳，说明岩土工程的稳定性与岩石裂纹扩展的时间相关：因此，研究岩石亚临界裂纹扩展规律对岩石稳定的时间相依性研究具有重要意义。目前，已有许多学者对岩石亚临界裂纹扩展进行了试验研究，如肖洪天等^[1]对三峡船闸高边坡花岗岩进行了双扭试验，得到了花岗岩亚临界裂纹扩展速度与应力强度因子之间的关系；张雯等^[2]分析了岩石亚临界裂纹扩展的应力腐蚀机理，认为亚临界扩展是由拉应力和裂纹尖端物质与环境中的腐蚀介质发生化学反应使化学键断裂这2种机制联合作用的结果，并测试了金川矿区的3种岩样，得到了应力强度因子与裂纹扩展速度之间的关系。袁海平等^[3]对金川矿区的软弱复杂矿岩进行了双扭试验，得到矿岩的断裂韧度及应力强度因子与裂纹扩展速度之间的线性关系；陈瑜等^[4-7]对岩石亚临界裂纹扩展进行了实验研究，并得到裂纹扩展的停滞速度和门槛值。汪亦显等^[8]得到了含水率较大的岩石亚临界裂纹扩展中应力强度因子与裂纹扩展速度的关系。上述研究均是从双扭试验中应力强度因子与柔度和应变能释放率之间的关系推导出应力强度因子的求解公式。本文作者运用Williams级数，对双扭试验中应力强度因子的数值解进行推导，并在上述学者所做试验研究的基础上，应用ANSYS数值模拟软件对岩石亚临界裂纹扩展进行数值模拟，将3种方法所得应力强度因子进行比较和分析，以便为岩土工程稳定的时间相依性研究提供基本参考。

1  双扭试验的基本原理

双扭试验是测试亚临界裂纹扩展速度的一种直接方法，它能够直观地监测裂纹的扩展过程，加载及预裂方式简单有效，适用性广，并且在确定应力强度因子K_Ⅰ时不需知道裂纹的长度，因此，该方法成为研究亚临界裂纹扩展的有效手段。双扭试件最先应用于研究玻璃、陶瓷等脆性材料的断裂特性，Ciccotti等^[9-12]将其应用于岩石材料，研究其断裂韧度及亚临界裂纹扩展规律。

图1所示为双扭试件示意图。该试件为一矩形薄板，在板的下表面沿长度方向的对称线开有1条通槽，以使加载时裂纹沿该槽扩展。

双扭试件可以看作由2个弹性扭转杆组成。对于小变形，且宽度远大于试件厚度的双扭试件，裂纹尖端的应力强度因子K_I^[9-12]为

式中：P为作用于扭杆上的荷载；w_m为扭臂的长度；μ为泊松比；d为试件的厚度；d_n为试件厚度与通槽厚度之差；w为试件的宽度。

图1  双扭试件示意图

Fig. 1  Schematic diagrams of double torsion specimen

式(1)为计算应力强度因子所用公式，可知应力强度因子与荷载、试件尺寸和泊松比等因素有关，而不依赖于裂纹长度a。

2  双扭试验裂纹尖端应力强度因子Williams级数解

2.1  作用力的转换

由图1可知：虽然双扭试验中试件端部受到力偶作用，但试件裂纹扩展是由于裂纹面受到垂直于裂纹面的应力作用，裂纹扩展原理与Ⅰ型裂纹的扩展原理相同，故可以在力偶作用面内用1对能产生相同裂纹扩展效果的拉力代替力偶。

力偶作用平面示意图如图2所示。根据弹性力学基本理论^[13]，力偶作用面内最大切应力为

式中：M为作用于杆端截面力偶；β为与d和w相关的系数。

图2  力偶作用平面示意图

Fig. 2  Plane of couple diagram

则模型底面任一点切应力为

式中：x为距坐标原点的距离；为切应力的最大值。可得底面切应力的合为

则模型所受到产生Ⅰ型裂纹的合力为

通过上述分析可知：双扭试验中裂纹的扩展由上述合力产生，故可用该力替代力偶。转换后的分析模型见图3。图3中：P为转换后的合力；r和θ为极坐标；b为试件长度；c为试件长度与裂纹长度之差。

图3  力偶转换后分析模型图

Fig. 3  Model diagram after changing couple

2.2  Williams级数表示的应力函数

力偶转换后，原双扭试验转换为平面问题，故可采用Williams级数进行分析。由文献[14]知，Williams提出的无穷级数应力函数分为偶函数和奇函数2部分，其中偶函数相当于Ⅰ型裂纹里对称加载，故可设图3所示模型的应力函数为

式中：为极坐标；C_j为待定系数。引入1个无量纲系数X_j^[14]，其与应力函数中第j项待定系数C_j的关系为

；j=1, 2, …, ∞         

式中：b为试件长度。根据直角坐标和极坐标的转换关系，可求得：

对于任意直角坐标系Ont与原直角坐标系Oxy之间夹角为ω，根据坐标系间关系可知

其二阶偏导数为

由Airy应力分量表达式^[14]可知

当ω=0时，，即

根据断裂力学原理可知

将式(13)代入，引入Williams极数，应力强度因子为

从式(15)可以看出：要解得应力强度因子，必须求解无量纲系数X₁。

2.3  边界配置法解方程组

在边界配置法中，为使方程可以求解，取应力函数前2m项组成的有限级数形式的近似解代替精确解，故式(13)可改写为

由于试样左右对称，故取1半进行分析。根据理论力学中力的简化可知：作用于裂纹口的1对拉力可转换为作用于模型长边中点的大小相等、方向相同的力和这个力与平移距离相乘得到的力偶。简化后的受力情况如图4所示。

图4  配置点设置图

Fig. 4  Diagram of collocation points

对于被研究部分m个配置点的设置如图4所示，在右边界设置4个点，编号为1，2，3和4；左边界设置3个点，编号为m-2，m-1和m；在上边界设置m-7个点，编号为5，6，…，m-3。

在边界配置法中，针对配置点建立边界条件，每个配置点有2个边界条件，即

；i=1, 2, …, m       

；i=1, 2, …, m         

式中：为非裂纹体的应力函数，

a为边裂纹长度；J为矩形截面惯性矩。则应力强度因子^[14]为

式中：c为试件长度减去裂纹长度，即c=b-a；λ₁和λ₂为与a和b相关的系数。

2.4  算例分析

为研究该计算方法的可行性，选取大理岩和花岗岩试块为例，具体计算试件尺寸及力学参数见表1。

表1  试样尺寸与力学参数

Table 1  Sizes and mechanical Parameters of specimens

根据李江腾^[16]试验数据中不同测试点的受力，经过前述力的转换可求得此计算方法中模型所受的单位厚度的力、力矩及开裂过程中裂纹长度a。大理岩试件具体参数见表2，花岗岩试件具体参数见表3。

将上述数据代入式(20)可得各试件计算点的应力强度因子。根据李江腾^[15]所测得的结果，并利用式(1)，可得各试件亚临界扩展时裂纹尖端的应力强度因子，见表4。

从表4可以看出：此计算方法所得的应力强度因子随着裂纹长度的增加而减小，与试验测得的应力强度因子变化规律相同；同时，计算所得的与试验测得的比较相对较大，并且随着裂纹占试样长度比例的增大而差距增大。

表2  大理岩试件计算点受力P、力矩M及裂纹长度a

Table 2  Force, couple and crack length of marble

表3  花岗岩试件计算点受力P、力矩M及裂纹长度a

Table 3  Force, couple and crack length of granite

表4  试验和计算所得应力强度因子

Table 4  Stress intensity factor of test and calculation      MN·m^-3/2

3  双扭试验的数值模拟

在试验和理论推导的基础上，本文应用ANSYS数值模拟软件，分别对上述大理岩和花岗岩试件进行常位移松弛法双扭试验的数值模拟，通过ANSYS软件分析，分别计算大理岩和花岗岩在亚临界裂纹扩展中裂纹尖端应力强度因子。具体模拟试件尺寸及力学参数见表1。

3.1  模型的建立

使用ANSYS计算应力强度因子，常规单元在裂纹尖端存在奇异性，为提高计算精度，本次模拟采用奇异单元，即二次三角(或五面体)单元，消除应力强度因子奇异性，提高计算精度。

由双扭试验可知，模型对称于沿长度方向所开的通槽，根据对称性，取试件整体的1/2建立几何模型。图5所示为试件网格划分模型。

图5  试样网格划分模型

Fig. 5  Meshing model of sample

在裂纹尖端点处，采用奇异单元。图6所示为裂纹尖端处网格划分模型。

图6  裂纹尖端处网格划分

Fig. 6  Meshing at crack-tip

3.2  加载求解

模型加载方式为：通槽所在边界处施加对称位移约束，反映另一半未画出模型对本模型的影响；在试验加载点对模型进行z轴方向的位移加载，模拟试验时的恒定位移；在钢珠放置点进行全位移约束，即限制钢珠放置点x，y和z方向的位移及转动。

在求解应力强度因子前，定义沿裂纹面的路径，以裂纹尖端作为路径第1点，沿裂纹面定义2个附加点，以指出裂纹扩展方向。

图7所示为试样1-1 y轴方向的应力云图。从图7可以看出：在裂纹尖端点处有明显的应力集中，与断裂力学理论相符。表5和表6所示分别为大理岩、花岗岩根据双扭试验测得的应力强度因子K_Ⅰ以及数值模拟所得的应力强度因子。

图7  试样1-1 y轴方向的应力云图

Fig. 7  Stress of y axis for sample 1-1

4  理论计算、模拟与试验结果对比 分析

大理岩和花岗岩根据试验测得的应力强度因子、Williams级数法计算的应力强度因子和模拟所得的应力强度因子的变化趋势分别如图8和图9所示。

表5  大理岩应力强度因子

Table 5  Stress intensity factor of marble

表6  花岗岩应力强度因子

Table 6  Stress intensity factor of granite

图8  大理岩应力强度因子变化趋势图

Fig. 8  Trend of stress intensity factor for marble

图9  花岗岩应力强度因子变化趋势图

Fig. 9  Yrend of stress intensity factor for granite

从图8和图9可以看出：岩石亚临界裂纹扩展过程中，随着裂纹的扩展，应力逐渐松弛，裂纹尖端的应力强度因子逐渐减小。

比较Williams级数法计算得到的应力强度因子和双扭试验测得的应力强度因子可知：Williams级数法计算得到的相对偏高，笔者认为虽然计算的式(15)是精确的，但式中无量纲系数X₁是近似解，因此，给计算结果引入误差。并且计算得到的随加载力的变化较小，笔者认为相对于试验法计算应力强度因子的公式，Williams级数法的计算公式会受到裂纹长度a的影响。在试验过程中，裂纹长度a变化较小，故由此计算法得到的应力强度因子亦变化较小。比较双扭试验测得的K_I和数值模拟所得的，可以看出数值模拟所得的应力强度因子相对偏低。这是由于：1) 模拟时对模型进行适当简化，使计算结果受到影响；2) 模拟软件使用的计算方法为有限元法，是一种近似解法，故与精确解有所差异；3) 由图1(c)可以看出试验时试件上、下表面的裂纹并不同时扩展，同一时刻裂纹尖端的应力强度因子在上、下表面并不相等，而模拟所得为均值。但从模拟所得与试验所得的拟合曲线可以看出：对于同一种岩石试样，3组试样的拟合曲线斜率基本一致，说明模拟与试验过程中应力强度因子减小规律相同，模拟能反映出试验的基本过程。因此，数值模拟的结果是可靠的，能为岩土稳定的时间相依性研究提供依据。

5  结论

1) 采用Williams级数及ANSYS数值模拟2种方法分别对大理岩、花岗岩在双扭试验中亚临界裂纹扩展进行了理论计算和数值模拟，得到了裂纹扩展时裂纹尖端的应力强度因子，为岩土工程稳定的时间相依性研究提供了基本依据。

2) 随着裂纹的扩展，应力逐渐松弛，裂纹尖端的应力强度因子逐渐减小，此减小规律符合断裂力学的相关理论。

3) 岩石一般所处地质条件较复杂，周围环境对岩石的亚临界裂纹扩展影响如水对裂纹亚临界扩展的影响较大，这在本文Williams级数计算和模拟分析中并没有得到体现；又如在地层深处，温度对岩石亚临界裂纹扩展的影响亦不能忽视。这些因素的影响还有待于进一步研究。

参考文献：

[1] 肖洪天, 杨若琼, 周维垣. 三峡船闸花岗岩亚临界裂纹扩展试验研究[J]. 岩石力学与工程学报, 1999, 18(4): 447-450.

XIAO Hongtian, YANG Ruoqiong, ZHOU Weiyuan. Testing study of subcritical crack growth of granite at the Three Gorges shiplock[J]. Chinese Journal of Rock Mechanics and Engineering, 1999, 18(4): 447-450.

[2] 张雯, 曹平, 袁海平, 等. 岩石亚临界裂纹扩展的应力腐蚀[J]. 土工基础, 2009, 23(2): 64-67.

ZHANG Wen, CAO Ping, YUAN Haiping, et al. Testing study of stress corrosion of the rock subcritical clack growth[J]. Soil Eng and Foundation, 2009, 23(2): 64-67.

[3] 袁海平, 曹平, 周正义. 金川矿岩亚临界裂纹扩展试验研究[J]. 中南大学学报(自然科学版), 2006, 37(2): 381-384.

YUAN Haiping, CAO Ping, ZHOU Zhengyi. Testing study of ubcritical crack growth of ore in Jinchuan Mine[J]. Journal of Central South University (Science and Technology), 2006, 37(2): 381-384.

[4] 陈瑜, 曹平, 张科, 等. 岩石亚临界裂纹扩展机制及工程门槛值研究[J]. 武汉理工大学学报, 2010, 32(19): 60-64.

CHEN Yu, CAO Ping, ZHANG Ke, et al. Research on subcritical clack growth and engineering thresholds of rocks[J]. Journal of Wuhan University of Technology, 2010, 32(19): 60-64.

[5] 万琳辉, 曹平, 黄永恒, 等. 岩石亚临界裂纹扩展及其应用研究[J]. 科技导报, 2010, 28(9): 42-46.

WAN Linhui, CAO Ping, HUANG Yongheng, et al. Subcritical crack growth of rocks and its applications[J]. Science & Technology Review, 2010, 28(9): 42-46.

[6] 李江腾, 曹平, 袁海平. 岩石亚临界裂纹扩展试验及门槛值研究[J]. 岩土工程学报, 2006, 28(3): 415-418.

LI Jiangteng, CAO Ping, YUAN Haiping. Study on subcritical crack growth and thresholds of rocks[J]. Chinese Journal of Geotechnical Engineering, 2006, 28(3): 415-418.

[7] 赵延林, 曹平, 万文, 等. 岩石裂纹亚临界扩展实验与压剪流变断裂模型[J]. 中南大学学报(自然科学版), 2014, 45(1): 276-286.

ZHAO Yanlin, CAO Ping, WAN Wen, et al. Rock cracks subcritical propagation test and compression-shear rheological fracture model[J]. Journal of Central South University (Science and Technology), 2014, 45(1): 276-286.

[8] 汪亦显, 曹平, 陈瑜, 等. 膨胀性软岩水腐蚀损伤断裂力学效应实验研究[J]. 中南大学学报(自然科学版), 2011, 42(6): 1685-1691.

WANG Yixian, CAO Ping, CHEN Yu, et al. Testing study on damage and fracture mechanical effects for swelling soft rock during water corrosion[J]. Journal of Central South University (Science and Technology), 2011, 42(6): 1685-1691.

[9] Ciccotti M, Negri N, Gonzato G, et al. Practical application of an improved methodology for the double torsion load relaxation method[J]. International Journal of Rock Mechanics and Mining Sciences, 2001, 38(4): 569-576.

[10] Saadaoui M, Reynaud P, Fantozzi G, et al. Slow crack growth study of plaster using the double torsion method[J]. Ceramics International, 2000, 26(4): 435-439.

[11] Nara Y, Kaneko K. Study of subcritical crack growth in andesite using the Double Torsion test[J]. International Journal of Rock Mechanics and Mining Sciences, 2005, 42(4): 521-530.

[12] Ciccotti M, Gonzato G, Mulargia F. The double torsion loading configuration for fracture propagation: An improved methodology for the load-relaxation at constant displacement[J]. International Journal of Rock Mechanics and Mining Sciences, 2000, 37(7): 1103-1113.

[13] 徐芝纶. 弹性力学简明教程[M]. 北京: 高等教育出版社, 2005: 173-174.

XU Zhilun. A concise course in elasticity[M]. Beijing: Higher Education Press, 2005: 173-174.

[14] 袁懋昶. 断裂力学理论及其工程应用[M]. 重庆: 重庆大学出版社, 1989: 174-196.

YUAN Maochang. The theory of fracture mechanics and application in engineering[M]. Chongqing: Chongqing University Press, 1989: 174-196.

[15] 李江腾. 硬岩矿柱失稳及时间相依性研究[D]. 长沙: 中南大学资源与安全工程学院, 2005: 66-77.

LI Jiangteng. Study on instability and time dependence of pillar in hard rock mass[D]. Changsha: Central South University. School of Resources and Safety Engineering, 2005: 66-77.

(编辑  陈灿华)

收稿日期：2013-11-10；修回日期：2014-01-21

基金项目(Foundation item)：国家自然科学基金资助项目(51374246，51274249)；湖南省水利厅科技项目(2010) (Projects(51374246, 51274249) supported by National Natural Science Foundation of China; Project(2010) supported by Science Foundation of Water Resources Department of Hunan Province)

通信作者：李江腾(1964-)，女，湖南株洲人，博士，教授，从事固体力学及岩土工程的教学与研究；电话：13973113109；E-mail: ljtcsu@163.com