高温后预制裂纹花岗岩损伤特性研究

2023-08-23王伏春黄守国

岩土工程技术 2023年4期

王伏春黄守国黄聪

（中南勘察基础工程有限公司，湖北武汉 430081）

0 引言

许多深部地下工程需要研究高温对岩石类材料力学与变形性能的影响，例如高放射性核废料的地质处置、增强型地热系统、地下煤气化及火灾后构筑物的修复重建工程等[1-7]。了解和掌握高温作用后岩石力学性质演变规律及损伤破坏特性，对围岩稳定性研究、地热热储层评估、岩石材料修复等具有重要意义。因此，高温作用后岩石的力学特性是岩石力学研究领域的关键课题之一。

花岗岩力学性质较稳定，渗透率较低，是许多新兴地下工程的理想介质。通过开展室内试验对高温后花岗岩物理力学特性已有广泛研究，包括单轴抗压强度、抗拉强度、弹性模量、纵波速度、体积、质量损失、孔隙度、渗透率、峰值应变、泊松比、导热系数、孔隙孔径变化规律、断裂韧性、微裂纹的扩展与发育、宏观破坏模式等[8-14]，并取得了丰硕的研究成果。然而，室内试验不能直观地显示热致微裂纹的产生及相互作用机理，且针对同一样品不可重复。近年来，随着计算机技术在岩土工程领域的应用与发展，学者们开始通过数值模拟的方法对热损伤岩石进行研究。Zhao 等[15]利用实验数据建立了颗粒离散元模型，在此基础上，揭示了北山花岗岩在压缩过程中出现负泊松比的原因。Sun 等[16]基于颗粒流数值模拟方法，对Lac du Bonnet 花岗岩热致微观和宏观裂纹过程进行了模拟，揭示了温度对花岗岩力学性能和声发射特性的影响。周喻等[17]利用颗粒流理论和PFC 程序，基于矩张量理论，建立了细观尺度上岩石声发射模拟方法，并应用该方法揭示了岩石的损伤破坏机制。宿辉等[18]同样利用颗粒流离散元程序，建立了片麻花岗岩热固耦合数值模型，并对不同温度处理后的岩石损伤演化规律进行了分析。Li 等[19]通过引入Knuth-Durstenfeld shuffle 算法，建立了考虑矿物组成的岩石非均质数值模型，模拟的岩石热损伤结果与实验结果一致，证明了该模型的可靠性。

节理作为一种重要的地质结构，普遍存在于花岗岩岩体中，对花岗岩的裂纹扩展和破坏模式有重要影响。然而，目前有关高温后含裂纹花岗岩损伤特性的研究较少，高温后预制裂纹花岗岩损伤演化机理尚不明确。因此，根据已发表文章中花岗岩实验数据，建立了颗粒流数值模型，模拟了不同温度作用后含裂纹花岗岩单轴压缩实验。在此基础上，对高温后预制裂纹花岗岩损伤特性进行了研究。

1 颗粒流数值模型的建立与验证

1.1 花岗岩模型构建

花岗岩数值模型的部分参数取自Huang 等[20]的实验研究成果。建立的模型中石英含量为23.8%，钠长石含量为52.5%，白云母含量为23.7%。模型尺寸高100 mm，宽50 mm，孔隙率为0.3%，包含大约1.1×104个颗粒，颗粒半径0.25～1.8 mm，服从均匀分布，颗粒密度为2730 kg/m3。预制裂纹位于试件中心，长度20 mm，宽度2 mm，与水平方向夹角为45°。由于PFC 中的颗粒不能进一步分割，所以预制裂纹的表面会出现局部粗糙现象。花岗岩试验样品与建立完成的数值模型如图1 所示。

图1 花岗岩试验样品与数值模型

1.2 力学参数的选取

PFC2D 中的平行黏结接触模型可以较好地应用于岩石模型[21]。因此，采用平行黏结模型进行建模，借鉴Shi 等[22]提出的标定方法对模型微观参数校准，数值模拟单轴压缩曲线与实际试验曲线对比见图2，可以看出模拟结果与实验结果匹配度较高。同时，花岗岩样品与数值模拟的宏观破坏模式也具有良好的一致性。校准后模型的细观物理力学参数见表1。

表1 模型细观物理力学参数

图2 试验结果与数值模拟结果

1.3 热学参数的选取

在PFC2D 热力学模块中，颗粒温度的变化会引起其半径的改变，以此产生热致应变，如式（1）所示

式中：ΔR为颗粒半径的变化量；α为颗粒的线性热膨胀系数；R为颗粒半径；ΔT为温度变化量。

温度产生的热致应变会使颗粒间胶结键的法向力发生变化，当胶结力增加到一定程度时，即达到模型所设定的拉伸或剪切强度，胶结键随即发生断裂，模型产生热致裂纹，具体公式为

式中：ΔFn为胶结法向力改变量；kn为胶结的法向刚度；A为胶结的横截面面积；L为胶结长度。

为了准确反映温度对花岗岩矿物组份的影响，需要对模型的各个矿物组份进行热导率、线性热膨胀系数、比热等热物性参数赋值。根据文献[23]的研究成果，花岗岩模型的热物性参数取值如表2 所示。

表2 模型热物性参数

在PFC 热力学模块中，颗粒之间的热通道具有热阻这一属性，不能直接赋值热导率参数，需要将热导率与热阻进行计算与转化，得到相应的热阻值赋给热通道，具体公式为

式中：η为热阻；k为热导率；n为孔隙率；V(b)为颗粒圆盘的体积；l(p)为热通道的长度。

2 研究方案

Isaka 等[24]研究表明，当温度梯度小于0.3 ℃/mm 时，可忽略由温度梯度引起的热应力，而5 ℃/min 的加温速率对应温度梯度约为0.27 ℃/mm。因此，本研究将模型的升温速率设置为5 ℃/min。此外，考虑到温度超过573 ℃石英会发生α-β 相变，由此产生体积膨胀，为了更好地模拟这一实际情况，参考已有研究经验[25]，当加热温度大于570 ℃时，将代表石英矿物的颗粒半径增大1.0046 倍。具体的模拟方案如下：

（1）通过颗粒分组、伺服消除内应力、赋值胶结参数、预制裂纹等步骤生成花岗岩模型。

（2）模型采用四周边界同时加温的方式，加温示意图如图3 所示，模型初始温度设置为20 ℃，温度边界随时间不断变化，变化速率为5 ℃/min，当温度边界达到预定温度水平后（设置了5 个温度水平：20℃，200 ℃，400 ℃，600 ℃，800 ℃），恒温至模型温度场均匀。加热过程中同时进行力学计算，循环模型直至平衡比极限达到0.0001。

图3 升温过程示意图

（3）为避免快速冷却对模型造成影响，温度边界的降温速率同样采用5 ℃/min，降至室温后再次恒温至温度场均匀。

（4）根据已有研究[26]，为模拟准静态的单轴加载过程，采用0.04 mm/min 的速率对模型进行单轴压缩加载，直至模型发生破坏，加载过程中时刻记录轴向应力及轴向应变的变化情况。

3 模拟结果与分析

3.1 力学性质

花岗岩模型单轴压缩的轴向应力-应变曲线见图4，大致可以分为压密阶段、弹性变形阶段、裂纹萌生阶段、裂纹扩展阶段及峰后变形阶段5 个阶段。随着温度升高，应力-应变曲线在初始变形阶段的非线性特性逐渐增强。但是数值模型没有考虑常温下花岗岩内部的原生裂纹，所以20 ℃的单轴压缩曲线没有体现出压密阶段。此外，随着温度升高，峰后仍呈现出脆性破坏特征，但应变软化现象趋于显著。

图4 不同温度水平下的应力-应变曲线

图5 是不同温度热处理后试件模拟结果的峰值应力、峰值应变及弹性模量的变化情况。从图中可以看出，峰值应力随热处理温度的增加呈现递减趋势，常温下的峰值应力为83 MPa，此后4 个温度水平下的峰值应力分别降低至80.5 MPa、76.1 MPa、64.6 MPa 和30.3 MPa。同样的，弹性模量随热处理温度的增加呈现下降趋势，与常温相比（34.6 GPa），此后4 个温度水平下的弹性模量分别降低了15.90%、20.52%、35.26%、63.58%。相反，峰值应变随热处理温度的增加呈现递增趋势，常温下的峰值应变为0.26%，此后四个温度水平下的峰值应变分别为0.275%、0.279%、0.286%和0.475%。

图5 不同温度水平下峰值应力、峰值应变及弹性模量变化

3.2 破坏模式

当颗粒之间所受作用力超过颗粒间的胶结强度就会发生断裂，断裂后通过调用Fracture 文件，可以实现在胶结断裂处插入一条短线（即一个DFN）来模拟微裂纹的产生。微裂纹逐渐累计，最终形成宏观破坏，即岩石的最终破坏是内部损伤累计的结果。

图6 显示了不同轴向应变水平下，热处理后预制裂纹花岗岩从微裂纹产生到试件完全失稳破坏的完整过程。图中绿色表示拉伸破坏，蓝色表示拉剪破坏，红色表示压剪破坏。从图中可以看出，所有温度水平下，拉伸破坏都占据主导位置。热处理温度为20 ℃、200 ℃、400 ℃和600 ℃时，加载初期，首先产生两个主要宏观裂纹，起裂位置位于预制裂纹的两端，裂纹扩展角度与预制裂纹角度呈大致垂直方向，随着轴向应变的增大，裂纹扩展长度不断增加。当加载到一定水平时，宏观裂纹数目增加，但仍主要起裂于预制裂纹的两个尖端，模型边缘部分也出现一些局部破坏，破碎程度增加，最终模型失去强度，形成宏观破坏。通过裂纹演化过程可以发现，该温度区间内，单轴压缩过程中裂纹扩展与预制裂纹的存在有关，主要宏观裂纹均起裂于预制裂纹的两端，与预制裂纹呈不同角度扩展。说明预制裂纹花岗岩的裂纹扩展和破坏模式受预制裂纹控制，与热处理温度无关，这与已有研究结果一致[20]。

图6 花岗岩数值模型的破坏模式

但是值得注意的是，800 ℃处理后的预制裂纹花岗岩试件裂纹扩展模式有所不同。如图6（e）所示，加载初期，裂纹并不是沿着预制裂纹两端起裂，而是非均匀的分布在试件内部。随着轴向应变增加，这些裂纹进一步扩展，并在其他位置产生新的裂纹，最终形成无规则的宏观破坏，破碎程度显著增加。此外，可以注意到，该温度水平下由单轴压缩产生的裂纹数量相比于200 ℃、400 ℃和600 ℃明显减少。分析上述现象的原因，800 ℃高温处理致使试件产生更多的热致微裂纹，且微裂纹数目和尺寸不断累积，导致热处理过程中试件就产生了许多宏观裂纹，这些裂纹严重影响了模型结构的完整性，从而热损伤主导了单轴压缩过程中预制裂纹花岗岩的宏观破坏模式，并导致单轴压缩过程产生裂纹数目较少，这与已有研究的结果是一致的[7]，即随着热处理温度增加，由单轴压缩过程产生的裂纹数目减少。

3.3 预制裂纹花岗岩热损伤特征分析

平行黏结模型不仅能传递颗粒间的力与力矩，还能体现热膨胀现象，可以更加准确地反映岩石的热力耦合行为。模型在加热及缓慢冷却至室温后，分析其内部应力的变化、微观损伤的产生、微裂纹数量及分布等情况，可以更好地理解试件的宏观力学性质及破坏特征。

图7 是模型试件加温前、高温状态和冷却后的应力分布与微裂纹分布图。其中红色表示拉伸应力，黑色表示压缩应力，蓝色表示微裂纹，力链的粗细代表力的大小。从图中看出，微裂纹的产生位置与应力集中区位置有较高的重叠率，表明高温导致花岗岩内部产生不均匀热应力，当热应力达到其胶结强度后就会产生破坏。

如图7（a）、图7（c）、图7（e）和图7（g）所示，加热前模型内部应力分布比较均匀，压应力远大于拉应力。当加热到高温状态时，模型内部应力分布不均匀性增加，出现较多的应力集中区（较粗力链形成的局部区域），这些应力集中区随机分布在模型内部。随着热处理温度升高，这些应力集中区会首先萌生微裂纹。为了更直观地研究热处理温度与热致微裂纹数量之间的关系，统计了热致微裂纹数量随加热温度的变化情况。如图8 所示，在400 ℃之前只有极少量的微裂纹产生（400 ℃时产生了6 个拉裂纹和1 个剪裂纹）。随着热处理温度增加，热致微裂纹的数量显著增加，且相互连通形成各大的裂纹，600 ℃时产生了741 个拉裂纹，124 个剪裂纹，800 ℃时产生1840 个拉裂纹，1360 个剪裂纹，且当热处理温度超过550 ℃以后，微裂纹数量迅速增加。升温过程中，微裂纹一方面是由于矿物的膨胀系数不同造成的，这也是本研究数值模型产生热致裂纹的主要原因；另一方面与高温导致的结合水蒸发及矿物组份化学变化等密切相关，这是本研究中200 ℃和400 ℃温度处理下产生热致裂纹数目较少的原因之一，但这超出了本文的研究范围。

如图7（b）、图7（d）、图7（f）和图7（h）所示，模型冷却降温后，应力集中区消失，应力分布变得相对均匀，并且拉应力区主要出现在模型的边缘部分，这与以往研究结果是一致的[27]。此外值得注意的是，200 ℃和400 ℃温度处理后，在试件预制裂纹的两端重新形成了两个明显的压应力集中区，这也是单轴压缩过程中，裂纹仍沿着预制裂纹两端起裂的原因，见图6（b）、图6（c）。600 ℃温度处理后，虽然预制裂纹两端没有出现应力集中现象，但是热处理产生的热致裂纹数目较少，单轴压缩过程中，破坏仍然从预制裂纹两端开始并扩展至整个试件，见图6（d）。相比之下，800 ℃温度处理后，热致裂纹显著增加，微裂纹之间相互连接贯通，形成更大区域的损伤，甚至宏观裂纹。因此，此时热致裂纹对轴压作用下试件破坏模式的影响会超过预制裂纹，见图6（e）。

对试件热处理后的微裂纹进行定量分析，结果如图9 所示。图中黑色表示剪裂纹，红色表示拉裂纹。从图中可知，微裂纹数量随着热处理温度的升高而显著增加。热处理温度为200 ℃时没有热致微裂纹产生，400 ℃时仅产生了8 个热致微裂纹。600 ℃时，热致裂纹数量开始急剧上升，试件产生了471 个剪切裂纹与1200 个拉伸裂纹，这可能与本研究在一定程度上考虑石英相变有关（573 ℃）。热处理温度为800 ℃时，试件产生了2193 个剪切裂纹与2204个拉伸裂纹。

综上所述，预制裂纹花岗岩数值模型在热处理温度为20 ℃、200 ℃和400 ℃时，高温对试件造成的损伤是非常有限的，而且高温处理导致试件预制裂纹两端形成应力集中，因此轴向荷载下，试件始终沿着预制裂纹两端开始破坏。热处理温度为600 ℃时，虽然热致裂纹数目显著增加，但是相比于预制裂纹，热损伤对试件影响较小，在轴压作用下，预制裂纹两端仍然会最先产生破坏。而热处理温度为800 ℃时，试件热损伤已经到达一定程度，甚至出现宏观破坏，试件强度显著降低，此时在轴向压力下，试件会沿着热损伤部位继续破坏，且破碎程度明显增加。