花岗岩是地壳中最常见的岩石之一，同时也是广泛使用的建筑材料，在宏观上表现出均质性，但在微观上因为矿物颗粒级配、成分比例以及微孔隙等的存在而呈现非均质性，从而在受热作用时容易产生微裂缝。微裂缝的逐步扩展、发育显然会影响岩石在热力作用下的宏观稳定性。

前人开展了大量的花岗岩受热破坏的研究，通过测量受热后花岗岩的弹性模量、力学强度、渗透率和地震波速度等参数，给出岩石热损伤的基本特征[1-11]。一般性的结论为，随着温度升高，花岗岩内部产生裂缝，而裂缝导致了渗透率的升高以及弹性模量、强度和地震波速度的降低。而且通过实验可观测到在300°C左右花岗岩物性参数有较明显变化，随后在石英发生相变的573°C发生更明显的变化[10-17]。这些实验研究侧重于花岗岩的均质化的热力学损伤和物性变化，并没有深入分析花岗岩内部非均质性的影响。

对微观尺度花岗岩微裂纹形态的观测方案包括：①二维高精度的显微观测（如SEM、偏光显微镜等），所观测的是某加载状态后某切面的静止结构[11,16,18]；②声发射记录可以提供三维裂缝动态发育过程观测，但是精度一般只能达到mm级别，无法确定裂缝与样品中矿物颗粒和先存孔隙之间的关系[17,19-22] ；③CT扫描提供高精度三维结构图像，分辨率一般为μm级别，单个CT扫描所提供的样品在某一状态的静止三维图像[23-26]。

原位实验与多次CT扫描结合可以获得岩石裂纹发展的系列动态图像，并可能识别裂缝与矿物颗粒及先存孔隙的关系。不过，目前可以识别花岗岩中原始孔隙的高精度观测并不多，高精度动态三维观测更是屈指可数[27-28]，而与高精度动态观测配套的动力学研究则更少。

热力作用下，花岗岩微观尺度非均匀结构导致其内部应力分布差异明显，这些内部应力差异无法观测，只能通过数值模拟方法获取。关于岩石样品内部结构与应力状态的数值模拟研究较多[25,29-33] ，广泛采用的有限单元方法可以考虑不同结构模型和本构，但基于连续介质力学的有限单元方法较难处理新生裂缝的形成扩展问题，离散元方法则在裂缝的形成扩展的模拟上具有得天独厚的优势。

根据矿物颗粒建模获得GBM（grain based model）并采用离散元进行模拟计算[34]，是矿物颗粒尺度非均质性特征研究的新技术。其中模型构建可以根据岩石结构的统计数据随机生成[35-36]，也可以根据SEM或者CT技术等获得的岩石结构真实图像实施[11,37-38]。但多数GBM建模仅考虑不同矿物颗粒，忽略花岗岩中先存的微小孔隙。Zhang等[38]尝试在模型中人为设置若干随机分布的先存裂缝，发现由此可以获得更接近实验曲线的应力应变关系，揭示了先存裂缝对岩石变形特征的影响。考虑天然微孔隙影响并与高精度动态CT观测对比的花岗岩热损伤数值模拟研究未曾见诸报道。

本文对一个Westerly花岗岩样品原位加热和多次扫描得到的四维（空间三维+时间演化）高精度CT图像[28]进行分析，获取其中矿物颗粒的分布以及孔隙、裂隙的扩展演化特征，并结合GBM离散元数值模拟，分析花岗岩非均质性对微破裂形成与演化的影响。

1.   原位动态CT观测与分析

1.1   图像数据及基本特征

本文所分析的Westerly花岗岩样品呈圆柱状，直径2.5 mm，高度8 mm。使用安装在芝加哥同步辐射光源（APS）CT扫描线站的加热装置，首先在室温下进行第1次扫描，2 min内快速加热到50°C进行第2次扫描；随后以1 °C/min的速率加热，设定每15°C为温度间隔，加热到设定温度后保持稳定10 min，然后在设定温度下进行CT扫描（耗时约25 min）；再循环加热、扫描。最高加热扫描温度为395°C，共获25组高分辨率三维岩石加热CT数据[28]。

受成像探测器的尺寸限制，仅扫描了样品一端（长度~2.7 mm），扫描分辨率为1.3 μm。每一个扫描步的数据均为2048张2048×2048体像素的.jpgf格式灰度图片，堆叠构成三维图像（图1（a））。

图  1  Westerly花岗岩初始状态CT扫描原始图像的正交切面显示（a）及图像切割后进行分析处理的圆柱状空间（b） ，其中右侧数字为原始数据体中水平切片编号 Figure  1.  (a) Original CT images of Westerly granite shown in orthogonal slices, and (b) the cropped cylinder in which horizontal slice numbers are labele

图像操作首先对图像进行切割。如图1（a） 所示的原始CT图像的顶、底部数百张切片需要去除；同时，假设用一个标准圆柱体对样品进行环状切割。切割后得到直径为1600体像素、高度为1280体像素的圆柱体（图1（b））。

选取中间3个水平切面（位置见图1（b））不同加热温度值时的图像进行观测对比（图2）。灰度图像中灰度值最小（接近黑色）的位置对应孔隙和裂隙。初始状态（25°C）时样品中就可以观察到较多先存孔隙，包含不同矿物边界（如图中红色箭头所示）和相同矿物内部（如图中黄色箭头所示）两类孔隙。215°C时，矿物边界及矿物内部孔隙较初始状态有明显增加，除了先存裂缝变宽、延长，局部位置出现新裂缝，如图中椭圆形所示位置。注意到一些新生裂缝出现于原始图像含较多微小孔隙或裂隙部位，如图中方框所示位置。后续升温过程中，大部分裂缝都沿着这些已形成的裂隙延伸扩展，可观察到裂隙的宽度有所增大，矿物边界及内部破碎加剧，395°C时裂缝明显更多。

图  2  原位加热Westerly花岗岩样品不同切片热裂纹扩展演化CT图像 Figure  2.  Three CT slices of in-situ heated Westerly granite showing the development of microcracks

1.2   图像分割与矿物识别

采用阈值分割方案对25个图像数据中的孔隙进行了分割，即将图像数据中灰度值最小、表现为近于黑色的所有体像素识别为孔隙。部分图像中存在的环状伪影（如图2（b）中395°C时圆心附近深色半圆环）有可能被分割为孔隙，所以对图像进行了细致的手动操作，以清除伪影造成的假孔隙，确保分割获得的孔隙准确可信。3个温度值对应的孔隙结构体渲染见图3，图中显示样品在初始状态即存在较多先存孔隙，且明显不均匀分布，随着温度的升高，内部孔隙逐步增加、裂缝增强。

图  3  Westerly花岗岩3个不同温度值时的孔隙结构体渲染图 Figure  3.  Volume rendering of the in-situ heated Westerly granite micropores and cracks at three different temperatures

图1和图2所示CT灰度图像中可见明显不同灰度的区域，这些不同灰度值皆因物质质量差异引起，即图像灰度分别对应不同矿物。Schrank等[28]对该组图像数据进行过矿物识别，本文参考其方案对初始扫描步的3张切片图像进行了不同矿物的分割，将孔隙以外的不同灰度分割为斜长石、微斜长石、石英和云母4种矿物。

对应图2初始状态的3个切片，4种矿物和孔隙均分割后的图像见图4。从图中可见，颗粒内部孔隙主要存在于斜长石和微斜长石矿物颗粒中，尤其在斜长石内部更多；石英矿物颗粒中可见少量裂隙，而云母矿物颗粒中几乎不存在孔隙。不同矿物颗粒边界之间存在大量明显先存裂缝，主要分布于石英−微斜长石以及石英−云母矿物边界，显然，不同切面上孔隙结构和各种矿物占比及分布差异均十分明显。

图  4  Westerly花岗岩CT扫描图像的矿物及孔隙分割结果 Figure  4.  Pore segmentation and different minerals in three slices of the granite sample

1.3   热裂纹特征

采用基于逾渗理论的微观结构定量分析程序CTSTA[39] 对不同温度下的孔隙进行分析。图像分割后由孔隙和固体构成的二值图像中，团簇定义为由共面的孔隙体像素构成的内部相互连通的结构，可以是相互交叉、相连的孔隙和裂隙所构成的大且复杂的形态，也可以是仅少量体像素构成的简单结构。CTSTA程序根据团簇定义实施团簇识别和标注后，输出数据体中每一个团簇的大小、形态和对孔隙度的贡献等参数，以及所有团簇的统计结果。连通性是孔隙结构的重要特征，当一个团簇将某方向两个相对的外边界面相连接时，定义该方向是连通的。由于圆柱状样品仅存在一对平行的上下边界，难以判别另两个方向的连通性，因此截取其中11303体积（最大内接正方体）进行了连通性的分析。

图5显示孔隙度、孔隙比表面积以及孔隙团簇数量均随温度的升高而总体增大，但是在65~95°C范围内，孔隙度和比表面积发生明显下降，团簇数量也有所下降。孔隙团簇数量在230~245°C大量降低，但这一时间段比表面积仅微小降低，孔隙度并未降低，只是增幅很小。结合CTSTA输出结果中不同大小团簇的数量和体积等参数变化，可以将花岗岩加热过程划分为以下阶段：

图  5  Westerly花岗岩样品孔隙度、孔隙比表面积（a） 及团簇数量 （b）随温度变化 Figure  5.  Porosity and specific surface area (a) and, number of pore clusters (b) versus temperature of the in-situ heated Westerly granite

①低温阶段（ < 65°C），孔隙度增加、最大团簇体积增加，小团簇（体积小于100体像素）数量几乎没有变化，意味着该阶段为空气升温使得孔隙膨胀，但并没有新的裂隙形成；②65~95°C之间，孔隙度、团簇数量和最大团簇体积、小团簇数量均减小，意味着矿物颗粒开始发生膨胀，造成孔隙体积减小，并造成了部分小孔隙的闭合；③95°C后，孔隙度、比表面积和团簇数量整体增加，最大团簇的体积和对孔隙度贡献率也增加，意味着裂缝的新生、扩展和相连。

230°C时孔隙实现3个方向连通；245°C团簇总数和小团簇数量均减小，同时比表面积略减，意味着大量独立的团簇合并到连通孔隙团簇中，但该温度新生的裂缝数量少于合并减少的裂缝数量；随后各参数均随温度上升，表明新裂缝生成速度大于裂缝合并速度，且随着温度升高新裂缝生成加速。

2.   GBM离散元数值模拟方法

为了模拟分析上述花岗岩热损伤过程，开展GBM离散元数值模拟。因高精度三维CT扫描数据量极大，离散元方法实施三维模型计算代价过高，本文采用二维模型进行计算，所获结果主要与图2和图5的裂缝发育规律对比。

根据图4中3个切面上4种矿物和孔隙的分布，构建了3个二维结构模型。综合考虑计算量和精度，将4×4像素大小（~5.2 μm）作为模型颗粒的平均粒径。将矿物晶体轮廓导入PFC2 D建立GBM模型；孔隙也以一种材料颗粒导入，但在矿物分组赋参前将其删除，因此模型中存在以空区表示的先存孔隙。构建完成的GBM与图4相似，即图4中不同颜色区域均被代表不同矿物的颗粒充填，孔隙部分则为空区。

模型中的颗粒假设为刚体，通过颗粒之间接触参数定义矿物晶粒和样品整体的力学特征，因此需要根据实验数据对颗粒参数进行标定。除了颗粒本身的力学参数，还需要定义颗粒之间的粘结参数，本文选用平行粘结模型。

前人开展过大量Westerly花岗岩相关实验研究。参考Petružálek等[40] 测定的Westerly花岗岩单轴压缩实验的应力应变曲线，本文通过数值模拟测试，在采用表1所列参数时，获得了与实验数据拟合很好的应力应变曲线，且模拟计算的破坏形态也与实验观测接近。热裂纹演化发展是一个热−力学耦合问题，还需要考虑热学参数。颗粒模型方法中，颗粒的热参数与宏观参数是一致的，本文采用与前人[28,41]实验和计算相似的热参数（表1）。

表  1  Westerly花岗岩试样矿物晶粒内颗粒微观热−力学参数 Table  1.  Thermal and mechanical parameters of different mineral particles for DEM

3.   模拟结果分析

3.1   天然（含孔隙）结构模型的热裂纹

采用表1中参数对图4中3个二维结构模型进行加热计算。结果显示，125°C时，模型内部仅产生极少量热裂纹，图中难以观测到；到215°C时各模型均出现了数十条微裂纹，与CT扫描图像所观察到的微裂缝扩展特征较为一致；300°C之后，Westerly花岗岩颗粒边界与颗粒内部裂纹数量急剧升高。

图6展示了3个温度时矿物颗粒、先存孔隙和新生孔隙的分布。显然，热裂纹主要沿着矿物颗粒边界以及先存孔隙进一步扩展，云母−石英边界上新生微裂纹最明显。

图  6  依据CT图像构建的真实模型热裂纹演化数值模拟结果（白色为先存孔隙，红色为新生热裂纹） Figure  6.  Simulated microcracks at different temperatures in 2D models constructed from CT sclies (white indicates pre-existing pores/cracks, red indicates newly formed pores/cracks)

模型中不同类型裂纹数量随温度的变化见图7。统计表明，拉伸型裂纹在全部破裂中占据主导地位；矿物晶界裂纹远多于晶内裂纹；不同矿物颗粒边界上，石英−云母边界新生裂缝多，其次为微斜长石−云母边界或微斜长石−石英边界。注意到3个不同结构模型的不同裂缝数量差异显著。

图  7  花岗岩不同切片GBM热裂纹数量特征 Figure  7.  Statistics of thermal cracks and temperatures in 3 GBMs of granite

3.2   无孔隙结构模型的热裂纹

为探讨模型中先存孔隙对热裂纹形成的影响，对图4结构模型进行孔隙充填，建立无先存孔隙的花岗岩颗粒模型。在相同的温度条件下，模拟计算了3个二维模型的热裂纹形成。结果显示热裂纹仍主要分布于矿物颗粒边界，热裂纹演化过程与图7相似（不另附图），但裂缝数量有所不同。

表2给出了考虑和忽略先存孔隙两种结构模型的Westerly花岗岩在不同温度时的热裂纹数量对比。数据表明，在125°C时，含先存孔隙的Westerly花岗岩模型已经出现极少量热裂纹；215°C时，含先存孔隙的花岗岩模型产生明显裂缝，而无孔隙的Westerly花岗岩内部几乎未产生任何裂纹，仅有1450切片在云母−微斜长石颗粒边界产生了一条裂纹；305°C时，含先存孔隙的Westerly花岗岩热裂纹数量远多于无孔隙的相同结构样品；但是到395°C时，热裂纹数量似乎与计算模型中是否包含先存孔隙无相关性。

表  2  有/无先存孔隙时不同温度下计算热裂纹数目对比 Table  2.  Number of thermal cracks at different temperatures in models with or without pre-existing pores

3.3   热应力及其演化

图8展示模型在加热过程中热应力的演化。其中，图8（a）~图8（c）分别为3个结构模型在395°C时的米塞斯应力分布，其余各幅图为1450切片局部（图8（c）中方框）的放大。

图  8  热应力分布及演化特征 Figure  8.  Thermal stress distribution and evolution

3个计算模型的米塞斯应力分布显示，石英矿物中的应力值明显高于其余3种矿物，并且在石英矿物狭窄连接部位应力值最高。

从表1可知，4种矿物成分中，石英矿物的热力学参数明显高于其他3种矿物，这势必造成石英矿物内部的高应力。类比Schrank等[28]的结果，热应力分布具有相似特征，但由于边界条件的不同，应力值大小无可比性。

图8 （d）~图8（f） 为1450切片模型局部区域分别在95°C、230°C和395°C时的应力分布。显然，随着温度升高，石英矿物的应力快速升高，并且应力升高部位逐步扩展到与其他矿物交界部位。由于离散元建模时不同颗粒之间的接触具有一定的随机性，同时模型中包含微小孔隙，高应力值位置在图中显示为较离散、且呈点状分布，但总体分布特征依然清晰。

前文图5 显示95°C时样品孔隙度最小，将之解释为矿物颗粒膨胀挤压孔隙空间，且该温度尚未产生裂缝。图8（d）显示尽管此时矿物晶体内的应力值还较低，与之对应的位移分布（图8 （g）），其显示出颗粒的膨胀确实已经造成颗粒从模型中心点向外缘的逐步扩展。

图8（h）和图8（i） 为1450切面模型中移除孔隙结构的计算结果，与图8（e）、图8（f） 对比发现，因天然样品中孔隙较小，有孔隙和无孔隙结构的模型总体应力分布差异不大；但是，无空隙模型局部区域的应力值高于有孔隙模型的相同部位，具体表现在：①椭圆所围的石英矿物区域内无孔隙结构模型近白色（应力值高于30 MPa）区域更多，②石英矿物的边缘和其相邻矿物的边缘部位应力值也更高，见图8（i） 中箭头所示部位与图8（f） 相应部位的对比。

4.   结论

通过对Westerly花岗岩加热过程高精度动态CT图像进行观测和定量分析，并开展热裂纹形成的GBM离散元二维数值模拟，二者关键特征具一致性，获得了微观尺度非均匀结构热损伤的基本特征。

CT观测发现原始状态的花岗岩中存在较多天然孔隙。受热伊始，气体膨胀导致孔隙增大，随后的矿物颗粒膨胀导致孔隙减小，温度~100°C起热裂纹开始形成，且随温度升高热裂缝形成加速、孤立孔隙相互合并、之后连通形成大的裂隙网络。裂缝发育主要沿矿物颗粒边界或以原始孔隙为起始点扩展，总体受岩石内部非均质结构的控制。

基于CT图像的GBM离散元数值模拟较好地再现热裂纹形成的动力学过程，其显示矿物晶粒和晶界力学性能决定了花岗岩的热损伤，微裂纹主要沿着矿物边界分布，且以张性破裂为主。受热后花岗岩中的石英矿物内产生更高的热应力，并逐步使与之相邻矿物的边缘应力也升高，造成石英矿物与其他矿物边界发生更多破裂。先存孔隙的存在使得岩石内部整体热应力值有所减小，同时显著降低了花岗岩加热后的起裂温度，其实含先存孔隙的花岗岩在200°C以下就已发生损伤。

未来可进一步开展三维数值模拟，以更好地再现花岗岩微裂缝扩展、连通的空间特征，并获得裂缝发育过程中岩石内部三维应力状态的演化规律。

1.中山大学地球科学与工程学院，广东 珠海519082

2.广东省地质过程与矿产资源探查重点实验室机构，广东 珠海519082

3.南方海洋科学与工程广东省实验室（珠海），广东 珠海519082

基金项目: 国家自然科学基金（基于岩石微观损伤的破裂失稳研究——从三轴实验同步动态CT观测到数值模拟（42072251））。