玄武岩三维细观孔隙模型重构与直接拉伸数值试验

工程科学学报，第41卷，第8期：997-1006,2019年8月 Chinese Journal of Engineering,Vol.41,No.8:997-1006,August 2019 DOI:10.13374/j.issn2095-9389.2019.08.005;http://journals.ustb.edu.cn 玄武岩三维细观孔隙模型重构与直接拉伸数值试验 郎颗娴12)，梁正召12)，董卓12) 1)大连理工大学海岸和近海工程国家重点实验室，大连1160242)大连理工大学岩石破裂与失稳研究所，大连116024 ☒通信作者，E-mail:LiangZZ@dlut.cdu.cm 摘要由于岩石材料的不透明性和多孔隙特性，通过传统的物理试验或数值模拟很难真实体现其内部三维细观结构.本文 基于CT扫描技术、边缘检测算法、滤波算法、三维点阵映射与重构算法，构建了可以表征玄武岩试样内部孔隙结构的三维细 观非均匀数值模型.结合并行计算进行直接拉伸数值试验，研究了内部孔隙结构特征对试样破坏机制及抗拉强度的影响.研 究结果表明：加载初期在试样孔隙处产生初始裂纹，随着荷载的增加初始裂纹逐渐沿横向扩展最终形成宏观拉伸破坏裂纹， 并且孔隙含量和分布位置对试样拉伸断裂的位置具有重要影响.随着孔隙率增高，试样破坏过程中的声发射数目和能量逐渐 减小.拉伸破坏模式呈现脆性破坏特征，同时孔隙的存在削弱了试样的抗拉强度 关键词玄武岩：孔隙；CT扫描：数字图像处理：有限元 分类号TG142.71 Three-dimensional microscopic model reconstruction of basalt and numerical direct tension tests LANG Ying-xian'),LIANG Zheng-zhao DONG Zhuo2) 1)State Key Laboratory of Coastal and Offshore Engineering.Dalian University of Technology,Dalian 116024,China 2)Institute of Rock Instability Seismicity Research,Dalian University of Technology,Dalian 116024,China XCorresponding author,E-mail:Liangzz@dlut.edu.cn ABSTRACT The presence of discontinuities and randomly distributed pores in basalt specimens greatly affects their engineering properties,such as the failure mechanism and strength.Therefore,investigating the mechanical and fracture behaviors of basalt affect- ed by the pre-existing defects is important for underground engineering,mining engineering,foundation engineering,and rock breaking and blasting.Laboratory tests have been widely used to research the failure mechanism of rocks under different conditions.However,it is difficult to clearly show the internal or spatial crack evolution during rock failure process in laboratory tests.Recently,X-ray comput- erized tomography (CT)and numerical tests have been used to detect the internal microstructures of rock specimens and to study their failure mechanism and strength.In addition,tensile strength is an important mechanical property of rock material.The direct tensile test is theoretically the simplest and most effective method for understanding the tensile behavior of rock.However,it is difficult to car- ry out in practical condition,because the sample processing and test procedures are complicated,also the experimental process of each sample cannot be repeated and has limited results.Due to the opacity of rocks,it is difficult to examine the three-dimensional internal structures of rocks through traditional physical and numerical experiments.In the present research,a 3D numerical method was pro- posed for simulating porous rock failure based on CT technology,the edge detection algorithm,filtering algorithm,and 3D matrix map- ping method.Direct tensile tests were carried out based on the parallel finite element method to study the effect of the porosity and pore distribution on the failure mechanism and tensile strength.The results indicate that initial cracks at the beginning of loading usually oc- cur in pores,and then with the raising of load the initial cracks propagate along the direction perpendicular to the loading direction and 收稿日期：2018-07-18 基金项目：国家重点研发计划资助项目(2016YFB0201000):国家自然科学基金资助项目(51779031,51878190)

工程科学学报,第 41 卷,第 8 期:997鄄鄄1006,2019 年 8 月 Chinese Journal of Engineering, Vol. 41, No. 8: 997鄄鄄1006, August 2019 DOI: 10. 13374 / j. issn2095鄄鄄9389. 2019. 08. 005; http: / / journals. ustb. edu. cn 玄武岩三维细观孔隙模型重构与直接拉伸数值试验 郎颖娴1,2) , 梁正召1,2) 苣 , 董 卓1,2) 1) 大连理工大学海岸和近海工程国家重点实验室, 大连 116024 2) 大连理工大学岩石破裂与失稳研究所, 大连 116024 苣通信作者, E鄄mail: LiangZZ@ dlut. edu. cn 摘 要 由于岩石材料的不透明性和多孔隙特性,通过传统的物理试验或数值模拟很难真实体现其内部三维细观结构. 本文 基于 CT 扫描技术、边缘检测算法、滤波算法、三维点阵映射与重构算法,构建了可以表征玄武岩试样内部孔隙结构的三维细 观非均匀数值模型. 结合并行计算进行直接拉伸数值试验,研究了内部孔隙结构特征对试样破坏机制及抗拉强度的影响. 研 究结果表明:加载初期在试样孔隙处产生初始裂纹,随着荷载的增加初始裂纹逐渐沿横向扩展最终形成宏观拉伸破坏裂纹, 并且孔隙含量和分布位置对试样拉伸断裂的位置具有重要影响. 随着孔隙率增高,试样破坏过程中的声发射数目和能量逐渐 减小. 拉伸破坏模式呈现脆性破坏特征,同时孔隙的存在削弱了试样的抗拉强度. 关键词 玄武岩; 孔隙; CT 扫描; 数字图像处理; 有限元 分类号 TG142郾 71 收稿日期: 2018鄄鄄07鄄鄄18 基金项目: 国家重点研发计划资助项目(2016YFB0201000);国家自然科学基金资助项目(51779031,51878190) Three鄄dimensional microscopic model reconstruction of basalt and numerical direct tension tests LANG Ying鄄xian 1,2) , LIANG Zheng鄄zhao 1,2) 苣 , DONG Zhuo 1,2) 1) State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China 2) Institute of Rock Instability & Seismicity Research, Dalian University of Technology, Dalian 116024, China 苣Corresponding author, E鄄mail: LiangZZ@ dlut. edu. cn ABSTRACT The presence of discontinuities and randomly distributed pores in basalt specimens greatly affects their engineering properties, such as the failure mechanism and strength. Therefore, investigating the mechanical and fracture behaviors of basalt affect鄄 ed by the pre鄄existing defects is important for underground engineering, mining engineering, foundation engineering, and rock breaking and blasting. Laboratory tests have been widely used to research the failure mechanism of rocks under different conditions. However, it is difficult to clearly show the internal or spatial crack evolution during rock failure process in laboratory tests. Recently, X鄄ray comput鄄 erized tomography (CT) and numerical tests have been used to detect the internal microstructures of rock specimens and to study their failure mechanism and strength. In addition, tensile strength is an important mechanical property of rock material. The direct tensile test is theoretically the simplest and most effective method for understanding the tensile behavior of rock. However, it is difficult to car鄄 ry out in practical condition, because the sample processing and test procedures are complicated, also the experimental process of each sample cannot be repeated and has limited results. Due to the opacity of rocks, it is difficult to examine the three鄄dimensional internal structures of rocks through traditional physical and numerical experiments. In the present research, a 3D numerical method was pro鄄 posed for simulating porous rock failure based on CT technology, the edge detection algorithm, filtering algorithm, and 3D matrix map鄄 ping method. Direct tensile tests were carried out based on the parallel finite element method to study the effect of the porosity and pore distribution on the failure mechanism and tensile strength. The results indicate that initial cracks at the beginning of loading usually oc鄄 cur in pores, and then with the raising of load the initial cracks propagate along the direction perpendicular to the loading direction and

·998· 工程科学学报，第41卷，第8期 eventually form macroscopic tensile cracks.The porosity and pore distribution have significant influences on the position of macroscopic tensile cracks.The acoustic emission (AE)event numbers and the accumulative AE energy are gradually decreased as the porosity in- creased.In addition,the brittle failure primarily determines the tensile failure mode and the presence of pores weakens the tensile strength of basalt samples. KEY WORDS basalt:pore:CT;digital image;finite element method 在自然界中岩石往往承受压缩荷载作用，因此 近年来CT扫描技术和数值试验方法的广泛应 研究者经常关注岩石在受压状态下的物理力学性质 用使得能够突破实验室客观实验条件的限制，非常 和破坏机制.但是岩石作为一种典型的脆性材料， 理想地实现各种加载条件，并且能针对试样内部孔 拉伸破坏是其破坏的基本模式之一，因此准确测试 隙结构对其力学性质和破裂过程的影响展开研究. 岩石在拉伸荷载作用下的力学性质并分析其产生的 Yu等8]利用CT扫描技术结合有限元方法研究了 破坏模式对边坡和隧道等地下工程具有重要的工程 节理岩体在压缩荷载条件下的力学破坏特性.J山 意义 等[]根据天然砂岩孔隙的统计分布和特征函数，利 目前，确定岩石在拉伸受力状态下性质的室内 用随机数算法和FLAC3D软件生成三维岩石孔隙数 物理试验方法主要有巴西盘间接劈裂试验和直接拉 值模型，研究了孔隙结构参数对试样力学性质的影 伸试验2-3).张少华等[)分别利用弧压劈裂、角压 响，但其研究方法为间接劈裂法.同时Wag等0] 劈裂和直接拉伸三种试验方法测试了不同类型岩石 采用蒙特卡罗方法生成随机骨料结构，并结合有限 抗拉强度，结果表明间接劈裂法得到的抗拉强度与 元方法模拟了混凝土材料中的裂纹扩展过程.但是 实际值相比偏小，同时指出在条件允许情况下，应尽 利用随机算法生成的数值模型不能准确的表征岩石 量采用直接拉伸方法测定岩石抗拉强度.窦庆峰 试样内部孔隙结构.于庆磊2]利用RFPA2D软件基 等[)对不同类型岩石试样进行抗拉强度试验，结果 于数字图像方法研究了细观结构对花岗岩断裂破坏 同样证明直接拉伸法较间接劈裂法优越.然而在直 特性的影响，其研究结果表明将CT扫描技术和RF- 接拉伸试验中存在试样加工困难且两端固定及试验 PA数值分析软件相结合可以得到良好的效果. 装置和试验过程复杂等问题[6，尽管Nova和Za- 因此，在上述研究的基础之上，本文利用CT扫 ninetti)]以及Okubo和Fukuit8]针对上述问题采用 描和数字图像处理技术，考虑岩石材料非均匀性，利 了各种不同的方法来提高试验结果的准确性，但是 用RFPA3D软件构建可以表征玄武岩内部孔隙结构 室内物理试验方法是一个花费高的、耗时的过程，且 的精细三维数值模型.结合有限元并行计算方法进 只能得到有限结果，同时室内试验往往需要借助高 行直接拉伸试验，深入研究玄武岩内部孔隙结构特 速摄影机和声发射技术研究试样的破坏过程9-]. 性对试样拉伸破坏模式、声发射特性和抗拉强度的 由于玄武岩内部含有各种不规则、跨尺度的孔 影响，为进一步了解玄武岩的拉伸破坏机制提供参考 隙和天然裂缝等缺陷，具有不均匀性和不可见性的 1构建数值模型 特点，这些缺陷的存在对玄武岩试样的宏微观破坏 模式和物理力学性质具有显著影响[].目前研 1.1试样与CT扫描处理 究者对玄武岩内部孔隙结构特征进行了一系列的研 本研究选取的玄武岩试样来自黑龙江省五大连 究.Bubeck等I),Schaefer等Iw以及Al-Harthi 池市，大约在地质年代的更新世晚期、第四纪早期， 等[]利用试验方法研究了孔隙率和孔隙形状特征 如图1所示.通过观察试样可以看出：试样1孔隙 对玄武岩试样工程性质和力学强度的影响，但是上 尺寸相对较大且分布稀疏：试样2表面孔隙含量较 述研究均是在压缩荷载条件下进行的.Canaker 少且尺寸较小：试样3孔隙较多且表面具有较大尺 等16)基于神经网络方法研究了玄武岩试样的抗压 寸孔隙：试样4孔隙尺寸较小但分布密集：试样5孔 和抗拉强度，但是其采用巴西劈裂法测定抗拉强度. 隙尺寸较小并呈“线状”分布.为获取玄武岩内部 同时Heap等[)利用数值试验方法研究了压缩荷载 孔隙结构信息，利用太原理工大学与中国工程物理 条件下孔隙率和气孔形状对火山岩脆性强度的影 研究院应用电子学研究所共同研制的 响，但是其数值模型是二维的且孔隙是随机生成的. μCT225 kVFCB型高精度显微CT试验系统对5个试 上述研究均未讨论试样内部孔隙结构对其最终破坏 样进行扫描.经统计分析得知本文选取的五个试样 形态的影响. 孔隙孔径主要集中在0~3mm范围内，而该cT扫

工程科学学报,第 41 卷,第 8 期 eventually form macroscopic tensile cracks. The porosity and pore distribution have significant influences on the position of macroscopic tensile cracks. The acoustic emission (AE) event numbers and the accumulative AE energy are gradually decreased as the porosity in鄄 creased. In addition, the brittle failure primarily determines the tensile failure mode and the presence of pores weakens the tensile strength of basalt samples. KEY WORDS basalt; pore; CT; digital image; finite element method 在自然界中岩石往往承受压缩荷载作用,因此 研究者经常关注岩石在受压状态下的物理力学性质 和破坏机制. 但是岩石作为一种典型的脆性材料, 拉伸破坏是其破坏的基本模式之一,因此准确测试 岩石在拉伸荷载作用下的力学性质并分析其产生的 破坏模式对边坡和隧道等地下工程具有重要的工程 意义[1] . 目前,确定岩石在拉伸受力状态下性质的室内 物理试验方法主要有巴西盘间接劈裂试验和直接拉 伸试验[2鄄鄄3] . 张少华等[4] 分别利用弧压劈裂、角压 劈裂和直接拉伸三种试验方法测试了不同类型岩石 抗拉强度,结果表明间接劈裂法得到的抗拉强度与 实际值相比偏小,同时指出在条件允许情况下,应尽 量采用直接拉伸方法测定岩石抗拉强度. 窦庆峰 等[5]对不同类型岩石试样进行抗拉强度试验,结果 同样证明直接拉伸法较间接劈裂法优越. 然而在直 接拉伸试验中存在试样加工困难且两端固定及试验 装置和试验过程复杂等问题[6] ,尽管 Nova 和 Za鄄 ninetti [7]以及 Okubo 和 Fukui [8] 针对上述问题采用 了各种不同的方法来提高试验结果的准确性,但是 室内物理试验方法是一个花费高的、耗时的过程,且 只能得到有限结果,同时室内试验往往需要借助高 速摄影机和声发射技术研究试样的破坏过程[9鄄鄄10] . 由于玄武岩内部含有各种不规则、跨尺度的孔 隙和天然裂缝等缺陷,具有不均匀性和不可见性的 特点,这些缺陷的存在对玄武岩试样的宏微观破坏 模式和物理力学性质具有显著影响[11鄄鄄12] . 目前研 究者对玄武岩内部孔隙结构特征进行了一系列的研 究. Bubeck 等[13] , Schaefer 等[14] 以 及 Al鄄Harthi 等[15]利用试验方法研究了孔隙率和孔隙形状特征 对玄武岩试样工程性质和力学强度的影响,但是上 述研究均是在压缩荷载条件下进行的. 覶anakc覦 等[16]基于神经网络方法研究了玄武岩试样的抗压 和抗拉强度,但是其采用巴西劈裂法测定抗拉强度. 同时 Heap 等[17]利用数值试验方法研究了压缩荷载 条件下孔隙率和气孔形状对火山岩脆性强度的影 响,但是其数值模型是二维的且孔隙是随机生成的. 上述研究均未讨论试样内部孔隙结构对其最终破坏 形态的影响. 近年来 CT 扫描技术和数值试验方法的广泛应 用使得能够突破实验室客观实验条件的限制,非常 理想地实现各种加载条件,并且能针对试样内部孔 隙结构对其力学性质和破裂过程的影响展开研究. Yu 等[18]利用 CT 扫描技术结合有限元方法研究了 节理岩体在压缩荷载条件下的力学破坏特性. Ju 等[19]根据天然砂岩孔隙的统计分布和特征函数,利 用随机数算法和 FLAC 3D软件生成三维岩石孔隙数 值模型,研究了孔隙结构参数对试样力学性质的影 响,但其研究方法为间接劈裂法. 同时 Wang 等[20] 采用蒙特卡罗方法生成随机骨料结构,并结合有限 元方法模拟了混凝土材料中的裂纹扩展过程. 但是 利用随机算法生成的数值模型不能准确的表征岩石 试样内部孔隙结构. 于庆磊[21] 利用 RFPA 2D软件基 于数字图像方法研究了细观结构对花岗岩断裂破坏 特性的影响,其研究结果表明将 CT 扫描技术和 RF鄄 PA 数值分析软件相结合可以得到良好的效果. 因此,在上述研究的基础之上,本文利用 CT 扫 描和数字图像处理技术,考虑岩石材料非均匀性,利 用 RFPA 3D软件构建可以表征玄武岩内部孔隙结构 的精细三维数值模型. 结合有限元并行计算方法进 行直接拉伸试验,深入研究玄武岩内部孔隙结构特 性对试样拉伸破坏模式、声发射特性和抗拉强度的 影响,为进一步了解玄武岩的拉伸破坏机制提供参考. 1 构建数值模型 1郾 1 试样与 CT 扫描处理 本研究选取的玄武岩试样来自黑龙江省五大连 池市,大约在地质年代的更新世晚期、第四纪早期, 如图 1 所示. 通过观察试样可以看出:试样 1 孔隙 尺寸相对较大且分布稀疏;试样 2 表面孔隙含量较 少且尺寸较小;试样 3 孔隙较多且表面具有较大尺 寸孔隙;试样 4 孔隙尺寸较小但分布密集;试样 5 孔 隙尺寸较小并呈 “线状冶分布. 为获取玄武岩内部 孔隙结构信息,利用太原理工大学与中国工程物理 研 究 院 应 用 电 子 学 研 究 所 共 同 研 制 的 滋CT225kVFCB 型高精度显微 CT 试验系统对5 个试 样进行扫描. 经统计分析得知本文选取的五个试样 孔隙孔径主要集中在 0 ~ 3 mm 范围内,而该 CT 扫 ·998·

郎颖娴等：玄武岩三维细观孔隙模型重构与直接拉伸数值试验 ·999· 描系统可以分辨微米级的裂隙和孔隙[]，因此满足 建模所要求的识别精度 试样1 试样2 试样3 试样4 试样5 图1玄武岩试样 Fig.1 Basalt samples 将试样沿试样高度方向自上而下扫描1500层， 图像，图中各像素点灰度值的变化范围为0~255 每层扫描间隔约为66.7μm,成像大小为2048像素 (黑色的灰度值为0，白色的灰度值为255)，不同灰 ×2048像素，扫描单元尺寸为70μm×70um.若试 度值代表不同的岩石结构类型.利用数字图像处理 样中孔隙尺寸小于CT扫描识别孔隙尺寸，则在CT 技术，可获取玄武岩试样内部孔隙结构信息，作为构 扫描图像中无法显示.图2(a)为试样1的CT扫描 建三维细观孔隙结构数值模型的基础. 图2CT扫描图像预处理.(a)CT图像：(b)二值化处理：(c)5×5滤波窗口：(d)10×10滤波窗口 Fig.2 CT images of sample 1 after processing:(a)CT images;(b)binarization processing;(c)5 x5 filtering window;(d)10 x10 filtering win- dow 1.2CT扫描图像预处理 能够显著区分岩石基质和孔隙介质，就认为阈值的 直接拉伸试验所需试样一般可分为以下三类： 选取是合理的.又因试样的CT扫描试验条件相同， 圆柱形、长方体、哑铃型(“狗骨形”)，本文选取长方 所得的CT扫描图像中岩石基质和孔隙介质的灰度 体试样进行直接拉伸试验.首先利用数学分析工具 值变化范围基本一致，因此可选取同一阈值进行二 MATLAB自主编写计算程序，对CT扫描图像进行 值化处理 剪裁、二值化与中值滤波处理，然后利用RFPA3D有 250 限元软件构建玄武岩三维细观孔隙结构数值计算 200 模型. 由于CT扫描试验获取的扫描图像外围存在黑 150 边，因此通过剪裁处理将黑边裁去（如图2(a)所 100 示).然后对图2(a)进行二值化处理以获取CT图 50 孔隙 --T=121 像中的孔隙结构信息.在图2(a)中随机选取一条 255075100125150175200225250 同时穿过岩石基质与孔隙介质的扫描线，统计该扫 像素点 描线上各个像素点的灰度值变化如图3所示，以此 图3扫描线上灰度变化曲线 来确定二值化图像分割阈值T.由图3可知，当灰度 Fig.3 Gray values along the scanned line processing 值大于所选阈值时表示岩石基质，当灰度值小于所 选阈值时表示孔隙介质及噪点，因此只要所选阈值 本文根据公式(1)对图2(a)进行二值化处理

郎颖娴等: 玄武岩三维细观孔隙模型重构与直接拉伸数值试验 描系统可以分辨微米级的裂隙和孔隙[22] ,因此满足 建模所要求的识别精度. 图 1 玄武岩试样 Fig. 1 Basalt samples 将试样沿试样高度方向自上而下扫描 1500 层, 每层扫描间隔约为 66郾 7 滋m,成像大小为 2048 像素 伊 2048 像素,扫描单元尺寸为 70 滋m 伊 70 滋m. 若试 样中孔隙尺寸小于 CT 扫描识别孔隙尺寸,则在 CT 扫描图像中无法显示. 图 2(a)为试样 1 的 CT 扫描 图像,图中各像素点灰度值的变化范围为 0 ~ 255 (黑色的灰度值为 0,白色的灰度值为 255),不同灰 度值代表不同的岩石结构类型. 利用数字图像处理 技术,可获取玄武岩试样内部孔隙结构信息,作为构 建三维细观孔隙结构数值模型的基础. 图 2 CT 扫描图像预处理. (a) CT 图像; (b) 二值化处理; (c) 5 伊 5 滤波窗口; (d) 10 伊 10 滤波窗口 Fig. 2 CT images of sample 1 after processing: (a) CT images; (b) binarization processing; (c) 5 伊 5 filtering window; (d) 10 伊 10 filtering win鄄 dow 1郾 2 CT 扫描图像预处理 直接拉伸试验所需试样一般可分为以下三类: 圆柱形、长方体、哑铃型(“狗骨形冶),本文选取长方 体试样进行直接拉伸试验. 首先利用数学分析工具 MATLAB 自主编写计算程序,对 CT 扫描图像进行 剪裁、二值化与中值滤波处理,然后利用 RFPA 3D有 限元软件构建玄武岩三维细观孔隙结构数值计算 模型. 由于 CT 扫描试验获取的扫描图像外围存在黑 边,因此通过剪裁处理将黑边裁去(如图 2 ( a) 所 示). 然后对图 2(a) 进行二值化处理以获取 CT 图 像中的孔隙结构信息. 在图 2( a)中随机选取一条 同时穿过岩石基质与孔隙介质的扫描线,统计该扫 描线上各个像素点的灰度值变化如图 3 所示,以此 来确定二值化图像分割阈值 T. 由图 3 可知,当灰度 值大于所选阈值时表示岩石基质,当灰度值小于所 选阈值时表示孔隙介质及噪点,因此只要所选阈值 能够显著区分岩石基质和孔隙介质,就认为阈值的 选取是合理的. 又因试样的 CT 扫描试验条件相同, 所得的 CT 扫描图像中岩石基质和孔隙介质的灰度 值变化范围基本一致,因此可选取同一阈值进行二 值化处理. 图 3 扫描线上灰度变化曲线 Fig. 3 Gray values along the scanned line processing 本文根据公式(1)对图 2(a)进行二值化处理. ·999·

·1000· 工程科学学报，第41卷，第8期 其原理如下：选取阈值T,若图像中的像素点灰度值 机数算法生成试样内部孔隙结构的传统方法，本文 大于或等于阈值T,则其灰度值被重新定义为1：若 建模方法更能真实表征岩石内部孔隙结构.根据上 像素点的灰度值小于阈值T,则其灰度值被重新定 述三维模型重构原理，构建如图5所示的直接拉伸 义为0.通过二值化处理，利用岩石基质与孔隙在 数值模型，通过与试样1的实体模型对比可知，利用 CT扫描图像中灰度值的差异来表征玄武岩试样内 该方法构建的数值模型可以精细表征玄武岩内部孔 部孔隙结构，如图2(b)所示 隙结构特征 1,fx,y)≥T f(x,)={0,fx,)<T (1) 式中，f(x,y)为像素点(x,y)的初始灰度值，f(x,y) 放大 为二值化处理后像素点(x,y)的灰度值. 通过图2(b)可知，二值化处理后的图像中含有 较多噪点，这些噪点的存在会对构建数值模型时孔 隙的识别产生影响，因此需要进一步对图像进行 中值滤波处理.中值滤波法是一种非线性平滑技 术，利用该方法可将图像中各个像素点的灰度值 设置为该点邻域窗口内的所有像素点灰度值的中 值.本文利用对比法确定中值滤波邻域窗口的大 小.图2(c)和(d)分别为5×5、10×10邻域窗口 处理的效果图，经对比可知采用10×10邻域窗口 图4细观图像转化为矢量化结构 进行中值滤波处理可以清晰的识别岩石基质和孔 Fig.4 Single-layer grid model based on a structural characterization 隙介质信息.按上述数字图像处理过程分别对试 样1~5的CT扫描图像进行处理，从而得到构建 直接拉伸数值模型所需的岩石基质和孔隙结构 信息 1.3三维模型重构 本文利用RFPA3D有限元软件构建可以表征玄 武岩三维细观孔隙结构的直接拉伸数值模型[23-2] 首先将数字图像处理后的图片中的信息转换为建模 所需的矢量化数据.如图4所示，数字图像是由正 方形像素点组成的.在三维空间中，若认为图像具 有一定的厚度t,则可将每一个像素点看作一个三维 有限元网格.将图片按照扫描次序进行排列叠加， 图5试样1实体模型与数值模型 当定义的每张图片的厚度1足够小时，其构建的三 Fig.5 Physical and numerical models of sample I 维模型与实际岩体中的孔隙尺寸误差可忽略不计. 本文模型尺寸50mm×50mm×100mm,模型规 将各个像素点的角点坐标转换为相应的矢量空间物 理位置（其中每个像素点的厚度为t,边长为a),并 模可达270万单元，单元尺寸约为0.5mm×0.5mm 根据像素点灰度值不同将其归类为不同岩石结构， ×0.36mm.直接拉伸数值模型岩石基质细观单元 并赋予相应的材料参数. 基本力学参数见表1，孔隙介质由空气基元[23-2]组 成，采用位移加载分布施加，位移加载量为每步 由于CT扫描精度的局限性，当孔隙尺寸小于 扫描精度时，在CT扫描图像中无法显示：其次，CT 0.001mm. 扫描图像中的噪点影响孔隙结构的识别：第三，计算 表1模型参数 能力的局限性，使得所构建数值模型可以表征的最 Table 1 Mechanical parameters of the numerical model 小孔隙尺寸即为最小单元尺寸.尽管上述三种影响 弹性模量/ 摩擦角/ 抗压强度/ 材料 泊松比 GPa () MPa 因素无法避免，但是本文数值模型中无法表征的微 玄武岩 哈 0.26 中 92 小孔隙对数值模拟结果影响不大，且相对于利用随

工程科学学报,第 41 卷,第 8 期 其原理如下:选取阈值 T,若图像中的像素点灰度值 大于或等于阈值 T,则其灰度值被重新定义为 1;若 像素点的灰度值小于阈值 T,则其灰度值被重新定 义为 0. 通过二值化处理,利用岩石基质与孔隙在 CT 扫描图像中灰度值的差异来表征玄武岩试样内 部孔隙结构,如图 2(b)所示. f忆(x,y) = 1, f(x,y)逸T 0, f(x,y) < { T (1) 式中,f(x,y)为像素点(x,y)的初始灰度值,f忆(x,y) 为二值化处理后像素点(x,y)的灰度值. 通过图 2(b)可知,二值化处理后的图像中含有 较多噪点,这些噪点的存在会对构建数值模型时孔 隙的识别产生影响,因此需要进一步对图像进行 中值滤波处理. 中值滤波法是一种非线性平滑技 术,利用该方法可将图像中各个像素点的灰度值 设置为该点邻域窗口内的所有像素点灰度值的中 值. 本文利用对比法确定中值滤波邻域窗口的大 小. 图 2( c)和( d)分别为 5 伊 5、10 伊 10 邻域窗口 处理的效果图,经对比可知采用 10 伊 10 邻域窗口 进行中值滤波处理可以清晰的识别岩石基质和孔 隙介质信息. 按上述数字图像处理过程分别对试 样 1 ~ 5 的 CT 扫描图像进行处理,从而得到构建 直接拉伸数值模型所需的岩石基质和孔隙结构 信息. 1郾 3 三维模型重构 本文利用 RFPA 3D有限元软件构建可以表征玄 武岩三维细观孔隙结构的直接拉伸数值模型[23鄄鄄24] . 首先将数字图像处理后的图片中的信息转换为建模 所需的矢量化数据. 如图 4 所示,数字图像是由正 方形像素点组成的. 在三维空间中,若认为图像具 有一定的厚度 t,则可将每一个像素点看作一个三维 有限元网格. 将图片按照扫描次序进行排列叠加, 当定义的每张图片的厚度 t 足够小时,其构建的三 维模型与实际岩体中的孔隙尺寸误差可忽略不计. 将各个像素点的角点坐标转换为相应的矢量空间物 理位置(其中每个像素点的厚度为 t,边长为 a),并 根据像素点灰度值不同将其归类为不同岩石结构, 并赋予相应的材料参数. 由于 CT 扫描精度的局限性,当孔隙尺寸小于 扫描精度时,在 CT 扫描图像中无法显示;其次,CT 扫描图像中的噪点影响孔隙结构的识别;第三,计算 能力的局限性,使得所构建数值模型可以表征的最 小孔隙尺寸即为最小单元尺寸. 尽管上述三种影响 因素无法避免,但是本文数值模型中无法表征的微 小孔隙对数值模拟结果影响不大,且相对于利用随 机数算法生成试样内部孔隙结构的传统方法,本文 建模方法更能真实表征岩石内部孔隙结构. 根据上 述三维模型重构原理,构建如图 5 所示的直接拉伸 数值模型,通过与试样 1 的实体模型对比可知,利用 该方法构建的数值模型可以精细表征玄武岩内部孔 隙结构特征. 图 4 细观图像转化为矢量化结构 Fig. 4 Single鄄layer grid model based on a structural characterization 图 5 试样 1 实体模型与数值模型 Fig. 5 Physical and numerical models of sample 1 本文模型尺寸 50 mm 伊 50 mm 伊 100 mm,模型规 模可达 270 万单元,单元尺寸约为 0郾 5 mm 伊 0郾 5 mm 伊 0郾 36 mm. 直接拉伸数值模型岩石基质细观单元 基本力学参数见表 1,孔隙介质由空气基元[23鄄鄄24] 组 成,采用位移加载分布施加,位移加载量为每步 0郾 001 mm. 表 1 模型参数 Table 1 Mechanical parameters of the numerical model 材料 弹性模量/ GPa 泊松比 摩擦角/ (毅) 抗压强度/ MPa 玄武岩 26 0郾 26 38 92 ·1000·

郎颖娴等：玄武岩三维细观孔隙模型重构与直接拉伸数值试验 ·1001· 表2模型计算孔隙率 2数值试样“计算孔隙率” Table 2 Porosities calculated by the numerical model 模型序号 计算孔隙率/% 在上述表征玄武岩三维精细孔隙结构数值模型 试样1 7.78 的基础上，利用统计学方法计算各个模型的计算孔 试样2 8.52 隙率，具体过程如下：通过图2(d)可知，预处理后的 试样3 25.75 CT扫描图像内的黑色像素点代表孔隙，通过统计试 试样4 13.95 试样5 10.43 样预处理CT扫描图像的黑色像素数目N,利用公式 (2)和公式(3)可计算出孔隙体积与模型的计算孔 隙率，计算结果如表2所示.从表2可知试样1孔 3 数值试验结果与分析 隙率最小，试样3孔隙率最大 3.1试样表面裂纹分析 V=N.V (2) 图6为5个试样在直接拉伸试验不同加载步 P.=方×100% (3) (Step)条件下产生拉伸破坏的最终裂纹形态，其中 同一试样的左右两图分别表示在不同视角条件下的 式中，P,为计算孔隙率，V为孔隙体积，N为试样黑色 破坏形态.通过分析可知：在孔隙含量较高的试样3 像素个数，V,为模型单元体积，V。为直接拉伸数值试 中部形成了一条相接的、贯通整个试样的横向宏观 样体积 拉伸破坏裂纹.而在其余试样中，主要产生相互啮 Z方向位移mm Z方向位移mm Z方向位移m Z方向位移ma -3037x101 -5.922x102■ -5.q22x102 -3.037×10 -2564x103 -5.011×102 5.011×102 -2564×10 -4099x102 4.099x102 -2.092x10-3 -2.092x101 -3.188×10 -3.188×10- -1.619x10 -1619x10 -2.277×103 -2277×12 -1.147×10-2 -1.147x10 -1366x10-2 -6.743x103 -6743x1 -1366×10- -4.545x10 4.545x103 -2.018×103 -2018x103 4568×10 4.568x103 2708x103 2.708x103 1368x10 1368×10- 7.433x10 7433x103 2279x103 2.279x10-2 1.216×102 1.216x10-■ Z方向位移mm Z方向位移 Z方向位移mm Z方向位移/mm -3.672×10-2 -3672×102 -3.436×102 -3.436×102 -3.010x10 -3010x10- -2.820x10 -2.820x10 -2348×10-2 -2348×10 -2205×10 -2205×10 -1.686x10 -1.686x10 -1.598×10 -1.598×10 1024x102 -1024×102 -9.735×10 -9.735×10 -3.618x10 -3.618×10 -3.578x10 -3.578x10 3.001×10-1 3001×103 2.578×10 2.578×10 9.621x10 9.621×10 8.734×10 8.734x10 1.24×10- 1.624×10 1.489x10-3 1.489x10- X2.286x10 2286x10 2.105x10 2.105×10 (d) Z方向位移mm Z方向位移mm -3467x101■ -3467×102■ -2.789x102 -2789x102 -2112x10 -2112x10 -1.435x10 -1.35×10 -7.572x10 -7.572x10- -7.977x10- -7.9T7×10 5.976x10 5.976x10 1.275x103 1.275x10 1.952x102 L.952×10 2.630x102 2630×102 y 图6玄武岩试样不同视角下拉伸断裂形态.(a)试样1(Step=18):(b)试样2(Stcp=12):(c)试样3(Step=15):(d)试样4(Step= 15):(e)试样5(Step=15) Fig.6 Crack morphologies of basalt specimens under direct tensile stress:(a)sample 1 (Step=18);(b)sample 2 (Step =12);(c)sample 3 (Step=15);(d)sample 4 (Step=15);(e)sample 5 (Step=15)

郎颖娴等: 玄武岩三维细观孔隙模型重构与直接拉伸数值试验 2 数值试样“计算孔隙率冶 在上述表征玄武岩三维精细孔隙结构数值模型 的基础上,利用统计学方法计算各个模型的计算孔 隙率,具体过程如下:通过图 2(d)可知,预处理后的 CT 扫描图像内的黑色像素点代表孔隙,通过统计试 样预处理 CT 扫描图像的黑色像素数目 N,利用公式 (2)和公式(3)可计算出孔隙体积与模型的计算孔 隙率,计算结果如表 2 所示. 从表 2 可知试样 1 孔 隙率最小,试样 3 孔隙率最大. V = N·V1 (2) 籽v = V V0 伊 100% (3) 式中,籽v为计算孔隙率,V 为孔隙体积,N 为试样黑色 像素个数,V1为模型单元体积,V0为直接拉伸数值试 图 6 玄武岩试样不同视角下拉伸断裂形态. (a) 试样1 (Step = 18); (b) 试样2 (Step = 12); (c) 试样3 (Step = 15); (d) 试样4 (Step = 15); (e) 试样 5 (Step = 15) Fig. 6 Crack morphologies of basalt specimens under direct tensile stress: (a) sample 1 (Step = 18); ( b) sample 2 ( Step = 12); ( c) sample 3 (Step = 15); (d) sample 4 (Step = 15); (e) sample 5 (Step = 15) 样体积. 表 2 模型计算孔隙率 Table 2 Porosities calculated by the numerical model 模型序号 计算孔隙率/ % 试样 1 7郾 78 试样 2 8郾 52 试样 3 25郾 75 试样 4 13郾 95 试样 5 10郾 43 3 数值试验结果与分析 3郾 1 试样表面裂纹分析 图 6 为 5 个试样在直接拉伸试验不同加载步 (Step)条件下产生拉伸破坏的最终裂纹形态,其中 同一试样的左右两图分别表示在不同视角条件下的 破坏形态. 通过分析可知:在孔隙含量较高的试样 3 中部形成了一条相接的、贯通整个试样的横向宏观 拉伸破坏裂纹. 而在其余试样中,主要产生相互啮 ·1001·