Thermal cracking mechanism of granite during heating and cooling processes
-
摘要: 由于高温试验设备的限制，在实验室中研究岩石真实热破裂规律通常需要依据冷却后试样的细微观结构观测结果来进行分析推断，无法获得热破裂的实时动态演化过程。为此，基于考虑温度和裂纹滑移效应的节理本构关系，构建花岗岩热力耦合UDEC晶粒模型，以深入探究加热和冷却过程中花岗岩的实时热破裂行为。研究发现，在加热状态下，花岗岩热诱导微裂纹在75℃左右开始萌生，随着温度的升高，微裂纹数量在α→β石英相变温度附近快速达到峰值，但微裂纹密度在600℃降温到25℃的冷却过程中并没有明显变化。虽然冷却效应引起的裂纹数量变化可以忽略不计，但会导致裂纹开度的增加或减少。在加热过程中，微裂纹的萌生主要是由于相邻晶粒的热膨胀不同导致局部应力累积造成的。而石英高温相变引起的细微观结构变化可以增强不同晶粒之间的相互作用，导致晶粒尺度上的压缩和剪切运动加剧，使得热应力裂纹继续变形拓展。在冷却过程中，由于加热过程中热破裂导致的细微观应力释放和冷却效应导致的不同矿物晶体收缩，使得微裂纹的数量虽然不变，但其形态却能发生较为明显的变化，从而使得冷却后花岗岩的宏观应力-应变行为受到很大影响。基于离散元数值模拟对花岗岩热力耦合试验结果作出了细微观阐释，揭示了加热和冷却过程中花岗岩的实时热破裂机制，进一步加深了对高温岩石热力耦合规律的理解。Abstract: Due to the limitations in high-temperature test equipment, the studies on the real thermal cracking of rocks in laboratory typically involve inverse analysis based on microstructure observations of cooled specimens. The real-time cracking evolution at high temperatures cannot be obtained through this method. Therefore, in this study, a thermo-mechanical coupled UDEC grain-based model for granite is established based on the modified joint constitutive law considering temperature and crack slip effects so as to investigate the real-time thermal cracking behavior of granite during heating and cooling processes. It is found that the thermally induced microcracking in granite begins to occur at around 75℃ under heating conditions. The number of microcracks rapidly increase near the α→β quartz phase transition temperature, but the microcrack density does not change significantly during the cooling process. Although the change in the crack number caused by the cooling effects is negligible, it can lead to an increase or decrease in the crack opening. During the heating process, the initiation of microcracks is mainly formed by the local stress accumulation due to different thermal expansions of the adjacent grains. The microstructure changes caused by quartz transition can enhance the interaction between different grains, leading to the increasing compression and shear motion on the grain level. This results in thermal-induced cracks continuing to deform and develop. During the cooling process, the local microscopic stress release due to the thermal cracking during heating and the shrinkage of different mineral crystals due to the cooling effects make the number of microcracks hardly change, but their morphological characteristics can change more significantly. This greatly affects the macroscopic stress-strain behaviors of granite after cooling. The findings of thermo-mechanical coupling tests on granite based on the discrete element numerical simulations are interpreted in a micro-meso-scale manner, revealing the real-time thermal cracking mechanism of granite during heating and cooling, further promoting the understanding of the thermo-mechanical coupling of high-temperature rocks.
-
Keywords:
- granite /
- grain-based model /
- real-time high temperature /
- thermal cracking
-
0. 引言
地热能源开发对于中国的绿色可持续发展和国家能源安全具有重要意义。在地热井建设过程中，需要在高温花岗岩地层中进行钻进，传统机械破岩方法面临效率低、钻头寿命短、钻速低等问题，使得钻井作业费用占总成本预算的比例可超过60%[1]。为了克服传统破岩钻进方式的不足，高温射流、微波、激光等常被用作物理辅助手段，通过高温效应使岩石产生劣化，进而加快破岩效率。这种劣化通过即时热量变化在岩石中诱发实时破裂（即实时高温热破裂），如果不能对实时热破裂机制有着较好的理解，就无法有效地发展和利用这些新型辅助破岩方式，甚至在未来工程实践中引起井壁围岩破裂失稳等灾难性后果。为此，研究花岗岩的实时高温热破裂机制对提高地热资源开发效率、促进高温辅助破岩技术的发展具有重要的意义。
目前，花岗岩在热处理后室温下的力学行为已经得到了广泛研究，但对实时高温下花岗岩热破裂规律的理解还不够充分[2-3]。近年来，越来越多的研究者在实时温度下（real-time temperature, 简称RT）进行了花岗岩热力耦合特性测试。例如，Yin等[4]对高温处理后（after heat treatment, 简称AT）和800℃实时高温下的花岗岩样品进行了单轴压缩试验。他们发现，与AT花岗岩相比，RT花岗岩延性转变的临界温度更低，峰值应力和弹性模量更小，而峰值应变和损伤更大。Kumari等[5]对花岗岩进行了实时三轴试验，试验围压为10~90 MPa，峰值温度为25~300℃。他们发现在选定的温度范围内，温度升高可导致花岗岩的强度升高和剪切参数降低。Wang等[6]对RT和AT花岗岩样品进行了单轴压缩试验，并比较了不同热处理条件下试样的应力-应变行为，发现RT试样的峰值应力和弹性模量低于AT试样，而峰值应变则表现出相反的趋势。Ma等[2]开发了一个实时高温真三轴系统，用来测试实时高温（最高400℃）下温度和水平应力对花岗岩力学行为的影响，他们发现RT对内聚力的影响比AT更大，并把这种现象归因于水平应力和实时高温的耦合效应。
RT和AT花岗岩的不同力学行为通常被解释为不同热载荷（如冷却过程）影响细微观热破裂行为进而导致试样内部细微观结构变化[7]。然而，由于当前设备无法对花岗岩实时热破裂进行直接观察，尚不清楚这些热应力裂纹的具体变形规律。目前进行岩石热破裂试验的主要困难是当前光学显微设备不能在高温条件下使用，通常只能对冷却至室温后的样品进行损伤观测。然而，高温花岗岩的裂纹网络是加热裂纹和冷却裂纹共同组成的，通过图像分析等手段无法区分裂纹是由加热引起的还是冷却引起的[8]。因此，只对冷却后试样进行观测无法获得热应力裂纹真正的实时高温演化过程。此外，由于冷却温度梯度变化、不同矿物热胀冷缩特性等的影响[9-10]，冷却过程可能会显著改变加热过程中形成的微观结构[11-12]。例如，Browning等[8]推测，在冷却过程中，岩浆岩中产生了更多的裂纹，并且根据加热和冷却过程中声发射信号的差异，发现这些裂缝的平均尺度更大。
综上所述，在热载荷作用下岩石内部的热破裂行为的动态可视化目前仍难以通过试验手段获得，只能根据冷却后试样的显微观察结果来进行定性推测。为了克服试验设备的不足，本研究将在试验数据的基础上采用离散元数值模拟手段对花岗岩的热破裂机制进行分析。多晶体花岗岩可以由离散UDEC块体模型表征，并通过块体接触相互连接和作用。UDEC块体可以分离，并自动检测新的接触，因此可以利用块体张拉和滑动来模拟矿物颗粒的热应力裂纹的产生和扩展[13]。鉴于此，本文将基于开发的花岗岩热力耦合UDEC晶粒模型来深入探究加热和冷却过程中花岗岩的实时热破裂规律及其对RT和AT试样宏观热力耦合特性的影响。
1. 数值模型建立
本节研究对象为德国的Eibenstock花岗岩（EG），室温下EG的密度为2.60 g/cm3，单轴抗压强度（UCS）为125.66 MPa[14]。作为典型的多晶体硬岩，其矿物成分主要为50%长石、44%石英和6%云母。这种非均质特征可以在UDEC中模拟为刚性或可变形块体单元的离散集合，从而可以在晶粒（或块体）边界张开处模拟裂缝的发育拓展[15]。此外，为模拟矿物晶体形态的随机性和岩石试样的非均质特性，在UDEC几何模型中利用其自带的Voronoi单元体生成功能，生成随机Voronoi单元，继而通过对这些Voronoi块体及其接触赋予相应的性质参数，就可以模拟不同的矿物组分（见图 1（a））。一旦局部应力超过接触的微观强度，代表不同矿物颗粒和接触的Voronoi块体边界就会张开或滑移，产生晶内或晶间裂纹。Mohr-Coulomb接触模型通常用于模拟块体接触的开裂行为[15]。在弹性范围内（即发生拉伸或剪切破坏之前），法向和剪切力增量（以压缩方向为正）根据以下公式计算[13]：
![]()
$$ \Delta {F^{\text{n}}} = - {k_{\text{n}}}\Delta {u^{\text{n}}}{A_{\text{c}}}\text{，} $$ (1) $$ \Delta F_i^{\text{s}} = - {k_{\text{s}}}\Delta u_i^{\text{s}}{A_{\text{c}}}。 $$ (2) 式中：$\Delta {F^{\text{n}}}$和$\Delta F_i^{\text{s}}$分别为法向力增量和剪切力矢量增量；${k_{\text{n}}}$和${k_{\text{s}}}$分别为法向和剪切刚度；$\Delta {u^n}$和$\Delta u_i^{\text{s}}$分别为法向位移增量和剪切位移矢量增量；${A_{\text{c}}}$为接触面积；i为直角坐标系中的分量。因此，接触处累积的拉伸法向力和剪切力矢量为
$$ {F^{\text{n}}} = {F^{\text{n}}} + \Delta {F^{\text{n}}} = {F^{\text{n}}} - {k_{\text{n}}}\Delta {u^{\text{n}}}{A_{\text{c}}}\text{，} $$ (3) $$ F_i^{\text{s}} = F_i^{\text{s}} + \Delta F_i^{\text{s}} = F_i^{\text{s}} - {k_{\text{s}}}\Delta u_i^{\text{s}}{A_{\text{c}}}。 $$ (4) 剪切力$F_{}^{\text{s}}$根据以下公式计算：
$$ {F^{\text{s}}} = {{\text{(}}F_i^{\text{s}}F_i^{\text{s}}{\text{)}}^{1/2}}。 $$ (5) 拉伸法向力和剪切力分别被限制为${T_{{\text{max}}}}$和${S_{{\text{max}}}}$：
$$ {T_{{\text{max}}}} = - {\sigma _{{\text{t0}}}}{A_{\text{c}}}\text{，} $$ (6) $$ {S_{{\text{max}}}} = {c_0}{A_{\text{c}}} + {F^{\text{n}}}{\text{tan}}{\varphi _{\text{0}}}。 $$ (7) 式中：${\sigma _{{\text{t0}}}}$为初始抗拉强度；${c_0}$为初始内聚力；${\varphi _{\text{0}}}$为初始内摩擦角。将内聚力和内摩擦角定义为抗剪强度参数。接触处的破坏（拉伸或滑动）是根据控制方程（式（6），（7））定义的。一旦检测到破坏，峰值强度参数（剪切和拉伸）被残余值所取代。在拉伸和剪切破坏开始时，接触弱化是通过内聚力和抗拉强度的快速损失来模拟的。残余拉伸法向力${T_{{\text{res}}}}$为
$$ {T_{{\text{res}}}} = - \sigma _{\text{t}}^{{\text{res}}}{A_{\text{c}}}。 $$ (8) 式中：$\sigma _{\text{t}}^{{\text{res}}}$为接触残余抗拉强度。一旦拉伸法向力超过${T_{{\text{max}}}}$，新的拉力修正为
$$ {F^{\text{n}}} \leqslant {T_{{\text{max}}}} \Rightarrow {F^{\text{n}}} = {T_{{\text{res}}}}。 $$ (9) 剪切应力不能超过残余强度，新的接触剪切力按以下方式修正：
$$ {S_{{\text{res}}}} = {c_{{\text{res}}}}{A_{\text{c}}} + {F^{\text{n}}}{\text{tan}}{\varphi _{{\text{res}}}}\text{，} $$ (10) $$ {F^{\text{s}}} \geqslant {S_{{\text{max}}}} \Rightarrow F_i^{\text{s}} = F_i^{\text{s}}\frac{{{S_{{\text{res}}}}}}{{{F^{\text{s}}}}}。 $$ (11) 式中：${S_{{\text{res}}}}$为剪切破坏后的最大剪切力；${c_{{\text{res}}}}$为残余黏聚力；${\varphi _{{\text{res}}}}$为残余内摩擦角。为简单起见，通常假设残余抗拉强度（$\sigma _{\text{t}}^{{\text{res}}}$）和黏聚力（${c_{{\text{res}}}}$）为零[13, 16]。这种抗拉强度和内聚力的瞬时损失通常应用于准脆性材料。
本研究中的晶粒模型采用了最新开发的考虑温度和裂纹滑移效应的节理本构关系[17]（见图 1（a））。对于屈服前的完整接触，考虑热效应的最大剪切力${S_{{\text{temp - max}}}}$和最大张拉法向力${T_{{\text{temp - max}}}}$分别为
$$ \begin{aligned} S_{\text {temp-max }} & =c_{\text {temp }} A_{\mathrm{c}}+F^{\mathrm{n}} \tan \varphi_{\text {temp }} \\ & =f_{\text {temp }-c} c_0 A_{\mathrm{c}}+F^{\mathrm{n}} \tan \left(f_{\text {temp }-\varphi} \varphi_0\right) \end{aligned} $$ (12) $$ {T_{{\text{temp - max}}}} = - {\sigma _{{\text{temp - }}t}}{A_{\text{c}}} = {f_{{\text{temp - }}t}}{\sigma _{{\text{t0}}}}{A_{\text{c}}}。 $$ (13) 式中：$ {c_{{\text{temp}}}} $，$ {\varphi _{{\text{temp}}}} $，${\sigma _{{\text{temp - t}}}}$分别为受温度影响的内聚力、内摩擦角和抗拉强度；$ {f_{{\text{temp - }}c}} $，$ {f_{{\text{temp}} - \varphi }} $，${f_{{\text{temp - }}t}}$分别为相应强度参数的温度依赖性函数。对于考虑温度和滑移效应的剪切破坏，其残余剪切参数按照以下方式进行修正：
$$ c_{{\text{temp - res}}}^{{\text{disp}}} = f_{{\text{res - }}c}^{{\text{disp}}}{c_{{\text{temp}}}} = f_{{\text{res - }}c}^{{\text{disp}}}{f_{{\text{temp - }}c}}{c_0}\text{，} $$ (14) $$ \varphi _{{\text{temp - res}}}^{{\text{disp}}} = f_{{\text{res}} - \varphi }^{{\text{disp}}}{\varphi _{{\text{temp}}}} = f_{{\text{res}} - \varphi }^{{\text{disp}}}{f_{{temp} - \varphi }}{\varphi _0}\text{，} $$ (15) $$ \begin{aligned} S_{\text {temp-res }}^{\text {disp }} & =c_{\text {temp-res }}^{\text {disp }} A_{\mathrm{c}}+F^{\mathrm{n}} \tan \varphi_{\text {temp-res }}^{\text {disp }} \\ & =\left(f_{\text {resc }}^{\text {disp }} f_{\text {temp-c }} c_0\right) A_{\mathrm{c}}+F^{\mathrm{n}} \tan \left(f_{\text {res }-\varphi}^{\text {disp }} f_{\text {temp- }-\varphi} \varphi_0\right) 。 \end{aligned} $$ (16) 式中：$c_{{\text{temp - res}}}^{{\text{disp}}}$和$\varphi _{{\text{temp - res}}}^{{\text{disp}}}$分别为考虑裂纹滑移和温度效应的残余内聚力和残余内摩擦角；$f_{{\text{res - }}c}^{{\text{disp}}}$和$f_{{\text{res}} - \varphi }^{{\text{disp}}}$分别为考虑温度和滑移效应的残余剪切函数；$S_{{\text{temp - res}}}^{{\text{disp}}}$为考虑裂纹滑移和温度效应的残余最大抗剪力。对于EG，考虑热效应和裂纹滑移位移的强弱化函数已通过实验室试验数据进行了反分析推算[17]。
在UDEC模型中，传热可以被模拟为各向同性瞬态热传导。热传导的基本定律为傅里叶定律，其中传热率与能量流动方向中的负温度梯度成正比[13, 18]：
$$ {Q_i} = - {k_{ij}}\frac{{\partial T}}{{\partial {x_j}}}。 $$ (17) 式中：${Q_i}$为i正方向上的热通量（W/m2）；${k_{ij}}$为导热张量（W/mK）；$\partial T/\partial {x_j}$为热流方向上的（负）温度梯度（K/m）；负号表示传导发生在温度降低的方向。任何物质的热容量是将该物质的温度升高1摄氏度或开尔文所需的热量，而物质的比热容则是单位质量该物质的热容量（J/kg/K）[19]。因此，温度变化的另一个基本定律可以写成
$$ \frac{{\partial T}}{{\partial t}} = \frac{{{Q_{{\text{net}}}}}}{{M \cdot {C_{\text{P}}}}}。 $$ (18) 式中：${Q_{{\text{net}}}}$为流入质量M（kg）的净热量（J/S）；${C_{\text{P}}}$为比热（J/kg/K）。对于二维传导，扩散方程可以通过结合方程（17），（18）获得，如下所示：
$$ \begin{aligned} \frac{\partial T}{\partial t} & =\frac{1}{\rho \cdot C_{\mathrm{p}}}\left[\frac{\partial Q_x}{\partial x}+\frac{\partial Q_y}{\partial y}\right] \\ & =\frac{1}{\rho \cdot C_{\mathrm{p}}}\left[k_x \frac{\partial^2 T}{\partial x^2}+k_y \frac{\partial^2 T}{\partial y^2}\right]。 \end{aligned} $$ (19) 式中：$ \rho $为质量密度。该方程控制数值模型中的热传导。对于相互接触的块体，热量可自由通过节理。对于已张开裂纹，为简化热传导过程，通过设置热公差范围来控制热量在不同块体间的流动[20]。在公差距离范围内的网格点被视为一个网格点，热流可以自由通过。
在模拟过程中，应力应变可随时与瞬态传热相耦合。由于准静态力学问题的能量变化通常可以忽略不计，因此不考虑应力引起的温度变化，而温度变化则会导致应力和应变的变化[13]，如下所示：
$$ \Delta {\sigma _{ij}} = - {\delta _{ij}}3K\alpha \Delta {\rm T},\Delta {\varepsilon _{ij}} = {\delta _{ij}}{\alpha _t}\Delta {\rm T}。 $$ (20) 式中：$\Delta {\sigma _{ij}}$为应力增量；$\Delta {\varepsilon _{ij}}$为应变增量；${\delta _{ij}}$为克罗内克函数（当i=j时，${\delta _{ij}} = 1$；当$ i \ne j $时，${\delta _{ij}} = 0$）；K为体积模量；$\alpha $为线性热膨胀系数；$\Delta {\rm T}$为温度变化量。在运行热模块的计算之前，需要对建好的模型（几何结构、材料模型和性质、边界条件、内部条件等）进行力学平衡求解，继而执行热力耦合计算过程[13]。
2. 模型校验
根据实验室试验[17]中采用的花岗岩试样尺寸，建立直径为50 mm、高度为110 mm的矩形模型（见图 2）。根据EG花岗岩的矿物含量百分比，随机定义Voronoi块体为不同的矿物组分，即50%长石、44%石英和6%云母。继而将不同矿物组分的热力性质（见表 1，2）赋值给相应矿物晶体及其接触。根据先前文献[12，21，22]中的数据，石英晶粒的线性热膨胀系数（CLTE）具有温度依赖性，并且由于可逆石英相变呈现出急剧热胀冷缩特性（见图 1（b））。研究发现α↔β石英相变温度具有约2.5℃的滞后[21]；因此在加热冷却期间α→β和β→α的石英相变温度分别设置为573℃，570℃。β→α相变的最大线性膨胀系数（750×10-6/℃）设置为略低于α→β相变的最大线性膨胀系数（800×10-6/℃）。以10℃/min的速率将设置好的花岗岩模型加热至400℃和600℃目标温度，并在目标温度下进行保温，直至整个模型的温度均匀分布（见图 2）。继而以-10℃/min的降温速率将模型缓慢冷却至室温，并对冷却后的试样进行单轴压缩模拟试验，如图 2所示。在单轴压缩模拟中，基于标准试样离散元模型的单轴压缩试验加载速度计算方程[20]，以0.021 m/s的恒定速度对试样的端部和底部进行加载，直至试样破裂失稳。在此过程中，通过试样两端的各9个均匀分布的监测点来实时收集应力-应变数据。如图 3所示，使用所开发模型模拟的应力-应变曲线和试样的最终破坏模式与实验室试验结果吻合良好，验证了数值模型的有效性。
性质(1) ρm/
(kg·m-3)νm Em/
GPaPct/
%性质
接触φm/
(°)cm/
MPaσtm/
MPakn/
(GPa·m-1)晶体 石英 2650 0.16 37.5 44 Q-Q(2)
M-M(2)30
1560
4013
7119802
51376长石 2620 0.19 22.5 50 F-F(2)
Q-F(2)25
27.550
5511
1274105
119802云母 3050 0.22 15 6 Q-M(2)
F-M(2)22.5
2050
4510
9119802
74105注：（1）微观密度($ {\rho _{\text{m}}} $)、微观杨氏模量($ {E_{\text{m}}} $)、微观泊松比($ {v_{\text{m}}} $)、微观内聚力($ {c_{\text{m}}} $)和微观摩擦角($ {\varphi _{\text{m}}} $)，数据源自文献[7]；（2）石英、云母、长石晶粒内的接触分别表示为Q-Q、M-M、F-F，晶粒间的接触为Q-F、Q-M、F-M，其强度为相邻矿物晶体的平均强度；接触法向刚度$ {k_{\text{n}}} = 5(k + {\text{4}}G/3)/\Delta {z_{{\text{min}}}} $，其中K和G分别是体积模量和剪切模量，$ \Delta {z_{{\text{min}}}} $是法线方向上相邻区域的最小宽度[13]；为求简化，晶粒接触切向刚度（$ {k_{\text{s}}} $）均设置为$ {k_{\text{n}}}/{\text{4}} $。 性质(1) αt(2) /
(10-6·K-1)k/
(W·m-1·K)Cp/
(J·kg-1·K)晶体 石英 13.40 7.70 750.0 长石 3.20 2.38 640.0 云母 12.81 2.01 769.5 注：（1）热性质$ {\alpha _{\text{t}}}({\text{CLTE}}) $，k（导热系数），$ {C}_{\text{p}}(比热) $，数据源自文献[7]，其中k和$ {C_{\text{p}}} $均设置为恒定值；（2）石英在冷却期间的CLTE遵循图 1（b）中的温度依赖性曲线[12, 21-22]。因为云母晶体仅占6%，并且没有经历可逆的相变，为求简化，其在冷却过程中的线性热膨胀系数被设定为加热过程中数值的1/2。 ![]()
图 2 花岗岩晶粒模型加热加载模拟试验方案Figure 2. Scheme of heating-loading simulation tests for granite grain-based model![]()
3. 加热冷却过程中花岗岩的热力学特性
图 4显示了热处理后实验室试样的p波速度以及加热冷却过程中数值模拟的热裂纹数量。在加热过程中，微裂纹数量随温度升高而逐渐增加，但在冷却过程中，该数量大致保持不变。模拟的加热和冷却裂纹数量之和与实验室的p波速度测量结果显示两者表征的热裂纹密度变化趋势具有高度一致性。模拟结果表明热处理后花岗岩中的裂纹数量主要是由加热应力引起的。在冷却过程中仅产生了少量的微裂纹，这是因为花岗岩的微破裂主要与热膨胀系数较小的矿物晶体（如长石）周围的累积张拉应力有关，一旦局部应力超过晶粒内或晶粒间的接触强度，就会产生新的微裂纹，进而使累积的局部集中应力得到释放，形成新的接触强度-应力平衡[20]。在冷却过程中，由于加热产生的随机热破裂遍布试样，使得绝大部分晶粒附近很难形成应力集中，因此很难产生新的微裂纹[20]。
![]()
图 4 加热冷却过程中微裂纹密度演化规律Figure 4. Evolution of simulated crack density during heating and cooling processes除了裂纹数量外，微裂纹类型和变形特征在加热-冷却过程中也发生了显著变化。由于矿物晶体的不断膨胀，尤其是石英在α→β相变时体积急剧增加，使得剪切裂纹数量在500℃至573℃的温度区间内持续增长（见图 4）。逐步活跃的剪切破裂行为表明，高温（尤其是石英相变温度）诱发的细微观结构变化可以增强不同晶粒之间的相互作用，导致晶粒尺度上的局部压缩、拉伸以及剪切运动加剧。图 5展示了实时高温和加热-冷却后花岗岩试样的剪切裂纹位移分布情况。在加热过程中，裂纹剪切位移随着温度的升高而不断增加。600℃试样的微裂纹剪切位移大多超过400℃试样的3倍以上。因此，加热过程中，在573℃ α→β石英相变温度之前，花岗岩中的细微观结构变化以新裂纹萌生和先前热应力裂纹的扩展为主。尽管微裂纹数量从600℃和400℃冷却至25℃的过程中没有明显变化（见图 4），但裂纹在冷却过程中会发生显著变形。在冷却期间，微裂纹的剪切位移既有增加又有减少（见图 5）。这意味着冷却过程中的裂纹活动比加热阶段更为复杂[3]，即加热阶段诱发热应力裂纹在冷却过程中会进一步滑动、收缩、加宽或闭合，同时也会产生少量的新裂纹。由于矿物晶体破裂释放了局部应力，冷却过程中新产生的微裂纹非常少。裂纹的剪切滑移或收缩变形主要是由于不同的晶粒变形和相互作用导致的。如果矿物晶粒在冷却过程中的体积显著减小，则相邻晶粒会受到变形晶粒形变位移的影响。因此，矿物晶粒在发生热破裂后仍可通过块体或碎片边缘的相互作用而产生移动。这些相互作用是随机而又无方向的，并且高度依赖于矿物的聚集模式，这就解释了为什么早期形成的裂纹可以同时发生滑动和收缩。
![]()
图 5 数值模拟中AT和RT试样应力应变曲线及其裂纹剪切位移分布规律Figure 5. Stress-strain curves and crack shear displacements of AT and RT samples in numerical simulations尽管冷却过程中的裂纹数量没有明显的变化，但微观结构的变化显著影响了热损伤花岗岩样品在单轴压缩下的宏观力学行为。对RT和AT样品进行单轴压缩试验的数值模拟，其应力-应变曲线如图 5所示。曲线的非线性压密阶段随着温度升高而增大，表明温度越高（比如600℃），样品的微裂纹密度越大。与加热后试样的应力-应变曲线相比，RT样品的斜率和峰值变小，这一特征也与大多数热损伤花岗岩的力学响应规律相一致[23]。
4. 花岗岩颗粒尺度热破裂特性
为了进一步验证EG花岗岩的热破裂规律是否也适用于其他花岗岩，本研究选用山东日照花岗岩的CT图像建立了细观热力耦合颗粒模型。日照花岗岩的单轴抗压强度（UCS）约为180 MPa，平均杨氏模量为31.7 GPa，其主要矿物组分为石英（37%）、斜长石（36.2%）、钾长石（20.6%）、云母（2.7%）、方解石（2.5%）和黏土矿物（1%）[24]。在实验室测试中，为避免过高的加热速度造成热冲击影响，花岗岩试样以2.5℃/min的速率缓慢加热至目标温度，并在目标温度下保持3 h以保证试样温度均匀分布，最后在加热炉内缓慢冷却至室温。冷却后的花岗岩样品通过Nano Voxel-2200系列高分辨率CT扫描系统进行扫描[23]。
4.1 基于CT图像的数值建模
利用数字图像处理（DIP）技术，在UDEC中复制花岗岩晶体颗粒的形态和分布特征。如图 6所示，首先使用图像处理程序ImageJ将CT图像转换为大小为100×100像素的8位灰度图像。8位灰度图像的亮度范围为[0，255]。每个像素具有256种可能亮度中的一种，即所谓的灰度值。然后利用ImageJ输出图像的像素数据矩阵并以ASCII形式保存。利用UDEC的程序读取功能，读取灰度值信息并赋值于相应位置的Voronoi块体（见图 6）。具有特定灰度值（即在定义的亮度范围内）的单元（块体）即为不同的矿物晶粒（见图 6）。对于日照花岗岩CT图像，阈值164和217可将灰度图像划分为与晶体形态吻合的图案。因此，灰色值在[0，164]、（164，217）和[217，255]范围内的块体分别被设定为长石、石英和云母。
![]()
图 6 基于CT图像生成晶粒模型步骤（在CT图像中，浅色区域表示高密度矿物，灰色区域表示低密度矿物；云母矿物显示为白色，长石显示为深灰色，石英矿物显示为浅灰色[24]）Figure 6. Scheme of generating numerical model from CT image (in the CT image, lighter color regions indicate high-density minerals, and gray color regions indicate low-density minerals; mica minerals are shown in white, feldspars in dark gray and quartz minerals in light gray[24])4.2 实验室和模拟结果
图 7展示了不同目标温度下加热冷却处理后花岗
![]()
岩试样截面的CT图像，其中，试样的孔隙和裂缝显示为黑色。400℃热处理后的花岗岩试样中，长石颗粒中仅产生了少量可见的微裂纹，而600℃热处理后的试样中微裂纹长度和宽度明显增大。数值模拟可以很好地再现热损伤花岗岩CT图像的热破裂模式和规律。为了更好地了解实时细微观热破裂机制，图 8（a）展示了不同加热和冷却阶段日照花岗岩模型中实时裂纹分布情况。微裂纹在75℃左右开始萌生，200℃时微裂纹数量逐渐增加，400℃时早期生成的微裂纹进一步扩展并伴有新的裂纹萌生。在600℃时微裂纹的数量达到最大值，但裂纹密度在整个冷却过程中没有明显变化。加热阶段诱发的热应力裂纹由于温度的升高而不断变宽，并在石英相变温度附近急剧增加。在冷却过程中加热阶段诱发的微裂纹产生较为明显的变形，尤其是收缩效应导致宽度略微减小，其中从600℃冷却到25℃裂纹宽度最大减少了3 μm左右（见图 8（b））。此外，裂纹的法向位移变化也说明了冷却过程中裂纹的张开或闭合等变形行为，且大多数裂纹的张开度随着材料的冷却收缩效应而减少（见图 8（b））。图 9（a）展示了加热冷却过程中实时最大主应力（正表示拉应力）的情况。由于石英具有较大的膨胀系数，其临近的长石晶粒中会产生拉应力来约束石英晶粒的膨胀。随着温度的不断升高，应力逐渐增大并进一步累积，产生更多的裂纹。由于局部的应力集中得到释放，冷却过程中几乎不产生新裂纹。因此裂纹的扩展或收缩主要由于晶粒尺寸水平上不同矿物的变形和相互作用所致。图 9（b）展示了不同温度下各矿物晶体的体积变化。长石晶粒在冷却过程中体积显著变小，相邻晶粒因此受到长石晶粒运动的影响，在块体边缘相互作用，接触破坏处的相邻矿物晶粒更易于发生移动，导致裂缝的扩大或收缩。这与EG花岗岩的试验模拟结果相一致。
![]()
图 8 加热冷却过程中花岗岩微裂纹的数量和法向位移Figure 8. Numbers and normal displacements of microcracks in granite during a heating-cooling cycle![]()
图 9 加热冷却过程中花岗岩的最大主应力（正表示拉伸）和体积应变Figure 9. Maximum principal stresses (positive for tension) and volumetric strains of microcracks in granite during a heating-cooling cycle5. 结论
本研究采用基于花岗岩热力耦合晶粒模型，深入探究了加热和冷却对高温花岗岩实时热破裂行为的影响机制。利用随机分布和图像处理技术，在UDEC模型中复制了花岗岩矿物晶体的形态特征，并从宏细观尺度对花岗岩的实时热破裂机制进行了揭示。
（1）在加热过程中，花岗岩的热应力裂纹的起始温度在75℃左右，新裂纹的发育和热应力裂纹的扩展是加热过程中花岗岩细微观结构的主要变化。其中，当温度超过400℃后，越来越多在较低温度下诱发的裂纹随着温度的升高逐渐加宽。由于α→β石英相变导致石英晶体颗粒体积的急剧变化，致使先前形成的微裂纹变形加剧。在冷却过程中，微裂纹数量没有明显变化，但裂纹宽度变化较为明显。冷却过程中的热破裂行为主要表现为加热阶段诱发裂纹的进一步加宽、收缩或闭合，以及冷却热梯度引起的少量新裂纹。
（2）加热和冷却过程热破裂机制的不同，主要和局部热应力的累积、释放以及颗粒尺度的材料变形有关。温度升高时，热膨胀系数较大的矿物晶粒体积增加趋势受到热膨胀系数较小的相邻晶粒的约束，诱发晶粒间接触应力累积，并在应力超过接触强度时产生微裂纹。温度下降时，由于加热阶段热破裂引发的集中应力释放，新的微裂纹很难形成。然而，冷却过程中的不同矿物晶体变形和相互作用仍然会导致裂纹张开度增加或减小。
总之，加热和冷却过程对热破裂活动都可产生显著影响，仅仅通过观察实验室中加热-冷却后岩石样品的裂纹分布规律，并不能准确地揭示实时高温下岩石的热破裂过程。数值模拟手段在揭示岩石实时热破裂机制方面发挥着不可或缺的作用。
-
![]()
![]()
图 2 花岗岩晶粒模型加热加载模拟试验方案
Figure 2. Scheme of heating-loading simulation tests for granite grain-based model
![]()
![]()
图 4 加热冷却过程中微裂纹密度演化规律
Figure 4. Evolution of simulated crack density during heating and cooling processes
![]()
图 5 数值模拟中AT和RT试样应力应变曲线及其裂纹剪切位移分布规律
Figure 5. Stress-strain curves and crack shear displacements of AT and RT samples in numerical simulations
![]()
图 6 基于CT图像生成晶粒模型步骤（在CT图像中，浅色区域表示高密度矿物，灰色区域表示低密度矿物；云母矿物显示为白色，长石显示为深灰色，石英矿物显示为浅灰色[24]）
Figure 6. Scheme of generating numerical model from CT image (in the CT image, lighter color regions indicate high-density minerals, and gray color regions indicate low-density minerals; mica minerals are shown in white, feldspars in dark gray and quartz minerals in light gray[24])
![]()
![]()
图 8 加热冷却过程中花岗岩微裂纹的数量和法向位移
Figure 8. Numbers and normal displacements of microcracks in granite during a heating-cooling cycle
![]()
图 9 加热冷却过程中花岗岩的最大主应力（正表示拉伸）和体积应变
Figure 9. Maximum principal stresses (positive for tension) and volumetric strains of microcracks in granite during a heating-cooling cycle
性质(1) ρm/
(kg·m-3)νm Em/
GPaPct/
%性质
接触φm/
(°)cm/
MPaσtm/
MPakn/
(GPa·m-1)晶体 石英 2650 0.16 37.5 44 Q-Q(2)
M-M(2)30
1560
4013
7119802
51376长石 2620 0.19 22.5 50 F-F(2)
Q-F(2)25
27.550
5511
1274105
119802云母 3050 0.22 15 6 Q-M(2)
F-M(2)22.5
2050
4510
9119802
74105注：（1）微观密度(ρm)、微观杨氏模量(Em)、微观泊松比(vm)、微观内聚力(cm)和微观摩擦角(φm)，数据源自文献[7]；（2）石英、云母、长石晶粒内的接触分别表示为Q-Q、M-M、F-F，晶粒间的接触为Q-F、Q-M、F-M，其强度为相邻矿物晶体的平均强度；接触法向刚度kn=5(k+4G/3)/Δzmin，其中K和G分别是体积模量和剪切模量，Δzmin是法线方向上相邻区域的最小宽度[13]；为求简化，晶粒接触切向刚度（ks）均设置为kn/4。 
下载: 导出CSV
下载: 导出CSV
-
[1] ROSSI E, KANT M A, BORKELOH O, et al. Experiments on rock-bit interaction during a combined thermo-mechanical drilling method[C]// 43rd Workshop on Geothermal Reservoir Engineering 2018. Proceedings of a Meeting Held 12-14 February 2018, Stanford, California, Curran Associates, Inc., 2019: 1874.
[2] MA X, WANG G, HU D, et al. Mechanical properties of granite under real-time high temperature and three-dimensional stress[J]. International Journal of Rock Mechanics and Mining Sciences, 2020, 136: 104521. doi: 10.1016/j.ijrmms.2020.104521
[3] YANG Z, YANG S Q, TIAN W L. Peridynamic simulation of fracture mechanical behaviour of granite specimen under real-time temperature and post-temperature treatments[J]. International Journal of Rock Mechanics and Mining Sciences, 2021, 138: 104573. doi: 10.1016/j.ijrmms.2020.104573
[4] YIN T B, SHU R H, LI X B, et al. Comparison of mechanical properties in high temperature and thermal treatment granite[J]. Transactions of Nonferrous Metals Society of China, 2016, 26(7): 1926-1937. doi: 10.1016/S1003-6326(16)64311-X
[5] KUMARI W G P, RANJITH P G, PERERA M S A, et al. Temperature-dependent mechanical behaviour of Australian Strathbogie granite with different cooling treatments[J]. Engineering Geology, 2017, 229: 31-44. doi: 10.1016/j.enggeo.2017.09.012
[6] WANG Z, HE A, SHI G, et al. Temperature effect on AE energy characteristics and damage mechanical behaviors of granite[J]. International Journal of Geomechanics, 2018, 18(3): 04017163. doi: 10.1061/(ASCE)GM.1943-5622.0001094
[7] WANG F, KONIETZKY H. Thermal cracking in granite during a heating-cooling cycle up to 1000℃: laboratory testing and real-time simulation[J]. Rock Mechanics and Rock Engineering, 2022, 55(3): 1411-1428. doi: 10.1007/s00603-021-02740-4
[8] BROWNING J, MEREDITH P, GUDMUNDSSON A. Cooling-dominated cracking in thermally stressed volcanic rocks[J]. Geophysical Research Letters, 2016, 43(16): 8417-8425. doi: 10.1002/2016GL070532
[9] NORDLUND E, ZHANG P, DINEVA S, et al. Impact of Fire on the Stability of Hard Rock Tunnels in Sweden[M]. Stockholm: Stiftelsen Bergteknisk Forskning-Befo, 2015.
[10] ZHANG B, TIAN H, DOU B, et al. Macroscopic and microscopic experimental research on granite properties after high-temperature and water-cooling cycles[J]. Geothermics, 2021, 93: 102079. doi: 10.1016/j.geothermics.2021.102079
[11] GLOVER P W J, BAUD P, DAROT M, et al. α/β phase transition in quartz monitored using acoustic emissions[J]. Geophysical Journal International, 1995, 120(3): 775-782. doi: 10.1111/j.1365-246X.1995.tb01852.x
[12] SIPPEL J, SIEGESMUND S, WEISS T, et al. Decay of natural stones caused by fire damage[J]. Geological Society, London, Special Publications, 2007, 271(1): 139-151. doi: 10.1144/GSL.SP.2007.271.01.15
[13] ITASCA. Universal Distinct Element Code[M]. Minneapolis: Itasca Consulting Group, Inc, 2018.
[14] WANG F. Thermal Induced Cracking of Granite: Laboratory Investigations and Numerical Simulations[D]. Germany: TU Bergakademie Freiberg, 2020.
[15] GHAZVINIAN E, DIEDERICHS M S, QUEY R. 3D random Voronoi grain-based models for simulation of brittle rock damage and fabric-guided micro-fracturing[J]. Journal of Rock Mechanics and Geotechnical Engineering, 2014, 6(6): 506-521. doi: 10.1016/j.jrmge.2014.09.001
[16] KAZERANI T, ZHAO J. Micromechanical parameters in bonded particle method for modelling of brittle material failure: micromechanical parameters in bonded particle method[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2010, 34(18): 1877-1895. doi: 10.1002/nag.884
[17] WANG F, KONIETZKY H, PANG R, et al. Grain-based discrete element modeling of thermo-mechanical response of granite under temperature[J]. Rock Mechanics and Rock Engineering, 2023, 56(7): 5009-5027. doi: 10.1007/s00603-023-03316-0
[18] FOURIER J B J, FREEMAN A. The Analytical Theory of Heat[M]. Cambridge: Cambridge University Press, 1878.
[19] WAPLES D W, WAPLES J S. A review and evaluation of specific heat capacities of rocks, minerals, and subsurface fluids. part 2: fluids and porous rocks[J]. Natural Resources Research, 2004, 13(2): 123-130. doi: 10.1023/B:NARR.0000032648.15016.49
[20] WANG F, KONIETZKY H, HERBST M, et al. Mechanical responses of grain-based models considering different crystallographic spatial distributions to simulate heterogeneous rocks under loading[J]. International Journal of Rock Mechanics and Mining Sciences, 2022, 151: 105036. doi: 10.1016/j.ijrmms.2022.105036
[21] SORRELL C A, ANDERSON H U, ACKERMANN R J. Thermal expansion and the high-low transformation in quartz. II: dilatometric studies[J]. Journal of Applied Crystallography, 1974, 7(5): 468-473. doi: 10.1107/S0021889874010223
[22] BALDO J B, DOS SANTOS W N. Phase transitions and their effects on the thermal diffusivity behaviour of some SiO2 polymorphs[J]. Cerâmica, 2002, 48: 172-177.
[23] WANG F, KONIETZKY H, HERBST M. Thermal effect of load platen stiffness during high-temperature rock- mechanical tests[J]. Computers and Geotechnics, 2020, 126: 103721. doi: 10.1016/j.compgeo.2020.103721
[24] FAN L F, GAO J W, DU X L, et al. Spatial gradient distributions of thermal shock-induced damage to granite[J]. Journal of Rock Mechanics and Geotechnical Engineering, 2020, 12(5): 917-926. doi: 10.1016/j.jrmge.2020.05.004
-
其他相关附件
计量
- 文章访问数: 535
- HTML全文浏览量: 78
- PDF下载量: 115
