基于CT 三维重建与逆向工程技术的煤体数字模型的建立

王 刚 ，杨鑫祥，张孝强，薛 娇，李文鑫

（1.山东科技大学 矿山灾害预防控制-省部共建国家重点试验室培育基地，山东 青岛 266590；2.山东科技大学 矿业与安全工程学院，山东 青岛 266590）

1 引 言

煤是一种多孔介质，它的微观孔隙结构对煤的宏观物理化学特性有着重要的影响，比如煤的孔隙尺度、比表面积以及连通性会影响煤的渗透率、孔隙度、吸水性等。研究人员对煤的认识过去多局限在宏观层面上，由于煤层气开采、井下水灾治理等工程问题均涉及煤的微观结构，国内外对煤的微观结构研究日益增多。由于多数的煤深埋于地下，通过观察法和实地探测法进行研究难度较大，研究人员多借助试验方法研究煤体的微观结构。常规试验方法主要有压汞法[1]、低温氮气吸附法[2]、切片法[3]和CT 法。CT 法与其他方法相比，不仅可以获得表征煤体孔隙结构的参数[4]，还可以在对煤体破坏尽可能小的情况下，定量表征煤体的微观孔隙结构，建立可视化的物理模型。

在对煤的微观研究中，CT 技术主要用3D 表征煤体的孔隙裂隙，比如煤层的割理及孔径间距[5]、构造煤渗流孔[6]、煤体内矿物的赋存特征[7]等。Simons 等[8]使用CT 重构技术与彩色图像分析技术对煤体进行了定量表征，表征的精度约为53 µm。Coles 等[9]利用 LBM 方法模拟研究了孔隙介质多相流，对砂岩表征的精度达到了30 μm。硬件的改善逐渐提高了3D 表征的效果，Van Geet 等[10]设计开发了微焦点CT 装置，最小观测精度可达0.35 μm，赵阳升等[11]使用的μCT225kVFCB 型高精度显微CT 系统可分辨1～2 μm 的孔隙和1 µm 宽的裂缝。其他研究人员通过CT 研究了煤体的受载破碎过程，Nie 等[12]通过设计开发的CT 系统研究了细观尺度下的煤体破裂过程。Li 等[13]通过CT 技术研究了煤体结构对煤体物理性质的影响，随压力升高将煤体结构的变化划分了5个阶段。

此外，对于煤体瓦斯的吸附与运移过程，Karacan 等[14]通过定量CT 技术进行了相应的研究，He 等[15]则通过煤样单轴压缩后可视化的孔隙及细观裂隙，进一步研究了煤与瓦斯耦合作用机制。研究人员也通过CT 技术对煤体的热破裂规律进行了探索，分析了次烟煤热干燥时的煤体孔隙空间的变化[16]、定量评价了煤热解时各影响因素的影响程度[17]、对瘦煤内部的热破裂过程进行了观测和分析，并得出热破裂发生的阈值温度[18]。其他研究手段的加入，弥补了CT 研究的不足，Cai 等[19]通过CT 重构技术结合声发射和超声技术研究了煤体在三轴加载时煤体破坏产生的裂隙及其对煤体渗透系数的影响，Yao 等[20]对比了压汞法（MIP）、恒速压汞法（CMP）、低场NMR（LENMR）与显微CT（µCT）在表征煤体孔隙结构上的差异。

CT 三维重建技术对煤体的微观结构的研究已日益成熟，但研究的内容多是使用重构后的模型对微观结构进行表征，进行定量分析与计算。基于重构后的三维模型进行数值模拟并进一步探讨其内部渗流规律的研究较少。一直以来，如何将CT 扫描数据用于数值模拟是摆在众多相关研究人员面前的一个难题。随着数字图形图像处理技术的发展，出现了一些商业三维数据可视化及建模软件，比如Avizo[21-22]、Simpleware[23-24]等，它们可以通过简单地操作完成CT 数据到三维网格的构建，其高度集成性与自动化，一方面为研究者提供了方便，节省研究时间，另一方面却造成整个转换过程的不可控，使得研究人员处于被动局面。为了更加灵活建立微观结构数字模型，由工业造型[25-28]得到启发，提出了逆向工程优化技术。

本文以高精度CT 系统扫描得到的CT 图片为基础，通过自行编写Matlab 三维重建程序，结合逆向工程技术，提出了一种将CT 三维数据转换为数字模型的方法，进而建立了煤体有限元模型并进行了煤体瓦斯渗流模拟。

2 CT 三维重建与逆向工程优化

2.1 CT 扫描试验

使用太原理工大学的µCT2225KvFCB 高精度CT 系统，对大柳塔煤矿的长焰煤原煤样品进行了CT扫描，煤质具有低灰、低硫、低磷和中高发热量的特点，属高挥发分长焰煤。在CT 扫描试验前，要对原始煤样进行前期加工处理，为了避免对煤样造成破坏，产生人为裂隙，将煤样制作成截面直径＜ 3 mm 的煤柱。进行CT 扫描试验时，煤样的放大倍数为100 倍，可观测到的最小孔隙直径为1.94 μm。所选试验设备及煤样尺寸如图1 所示。

2.2 煤体三维重建

三维重建技术就是将二维平面数据经过一定处理形成三维数据并进行显示。三维重建技术早期用在医学上，进行解剖显示及假肢设计等。随着CT技术和图形处理技术的发展，三维重建技术开始更多的应用在工业上，形成了以几何描述为目的的表面重构技术。三维重建的方法主要分为体绘制法和面绘制法[31-33]。使用Matlab 语言，编制了基于体绘制法的三维重建软件。

阈值分割是决定重构后的模型能否精确描述实际物理模型的关键一步，对于煤体的三维重建，随着CT 扫描环境及扫描试件的不同，会生成具有不同颜色特征的CT 图像。针对CT 试验扫描得到的CT 图像，结合扫描试件孔隙结构特征，对比多种阈值分割方法，最终采用了OSTU 阈值选择算法。重构经验表明，使用OSTU 算法重构得到的孔隙率偏大。同时，考虑到煤体重构后的表面比较粗糙，后续的逆向工程优化过程的光滑操作会造成重构后模型的体积损失，因此，在选取合理阈值时，应将后续优化的体积损失也考虑在内。

使用μCT225kVFCB 型CT 试验系统得到的CT图片的大小为2 041×2 041 像素，前期的加工处理会对煤样产生影响，为了能够反映煤体的微观结构并且兼顾数据量过大带来的重构计算压力，对CT图片进行了裁剪操作，选用了4个面积为60×60 像素的区域，对煤体内的孔隙进行了三维重建。重构后的三维模型大小为60×60×60 像素，实际大小为0.12 mm×0.12 mm×0.12 mm。重构完成后，输出stl格式的三维数据文件，将体素模型转换为基于三角面片的表面模型（见图2）。

图2 CT 切片及重构区域Fig.2 CT slice and the reconstruction areas

2.3 逆向工程优化

为了使得三维重建得到的模型能够导入主流的数值模拟软件进行更加深入的研究，需要将重构后模型进行优化处理，矢量化表面模型。当今主流数值模拟软件如Ansys、COMSOL和Abaqus 等均支持以NURBS 曲面数据为基础的CAD 数字模型[34]，逆向工程技术即提供了一种CT三维数据到CAD 数字模型的转换方法。

苏轼《守岁》诗说：“儿童强不睡，相守夜欢哗。”[28]128 旧时守岁要有火。 屋内炉火越旺越好，以辟邪除恶，也预示来年红火兴旺。 民国《晋祠志》记录除夕时“庙祝洁除内外，更贴楹联。 至夕陈设祭品，烧炬以待五鼓迎神，闾巷亦然。”[29]38 山西煤炭资源丰富，有架火塔、垒旺火的习俗。 午夜时鞭炮齐鸣，家长将旺火点燃，顿时火势熊熊，“旺气冲天”。 男女老少齐来观赏，儿童雀跃，欢欢喜喜。 陕北神木等地的旺火也很出名。

一般情况下，工业产品设计过程是一个“从无到有”的过程，可称为“正向工程”。逆向工程则是一个“从有到无”的过程。就是根据已经存在的产品模型，获取产品的设计数据，进行二次设计和生产的过程。CT 三维重建技术也逐渐成为一种获取产品设计数据的方法。逆向工程数据处理流程如图3 所示。

图3 逆向工程优化流程图Fig.3 Flowchart of optimizing reverse engineering

通常将坐标化的CT 三维数据称为点云数据，对由点云数据建立的三角形网格逐步优化得到平滑曲面，细化曲面得到子曲面，用NURBS 曲面逐一拟合子曲面，得到保持一定连续性的曲面样条，生成以NURBS 曲面表示的CAD 模型，完成CT 三维数据到CAD 数字模型转换。

使用Geomagic 逆向工程技术软件对重构后的模型进行了优化处理，首先要读取三维重建后得到的stl 文件（图4(a)），对可能存在的拓扑错误进行修复并光滑表面[35]（图4(b)），然后提取三维模型的轮廓线（图4(c)），以轮廓线为边界，将整个三维模型划分为多个曲面片（图4(d)），随后对曲面片进一步划分形成格栅（图4(e)），最后拟合生成NURBS曲面（图4(f)），形成CAD 数字模型，完成整个优化工作。

图4 逆向工程优化处理Fig.4 Reverse engineering optimization

常用的逆向工程软件多数具有数据查询及统计分析算法，这极大地方便了研究人员对三维数据的信息获取。例如煤体孔隙模型的体积与表面积均可以通过统计分析算法获取。利用逆向工程软件Geomagic 获取了煤体孔隙模型优化处理前后的孔隙度，如表1 所示。

表1 优化处理前后孔隙度变化值Table 1 Changes of porosity values before and after optimization

优化处理前后的孔隙度偏差值保持在3%左右，因此，在三维重建选取阈值时，应考虑到3%孔隙度损失。

3 煤体瓦斯渗流模拟

煤体瓦斯渗流即指瓦斯在煤体内的流动。煤体瓦斯的渗透系数是指瓦斯通过煤体的难易程度，它直接影响煤体中瓦斯的赋存量、运移量和涌出强度，同时也是瓦斯抽放难易程度及煤与瓦斯突出强度的重要参数。因此，对煤体渗透系数的研究对预防瓦斯动力灾害、煤层气开采及提高瓦斯抽放效率有着重要意义。

3.1 煤体有限元模型建立

将逆向工程优化后生成的CAD 数字模型导入到ICEM CFD 软件中进行网格划分。ICEM CFD 可以输出多种网格供CFX、Fluent、Abaqus、Nastran等软件使用，生成的网格具有一定的通用性，且其对三维CFD 网格的支持较好。由于煤体孔隙结构不规则，初始网格划分采用全局自动生成，对边缘网格进行细化处理，删除部分低质量网格，对局部网格进行特殊处理，逐步提高网格质量。将图2 中A区域CT 三维数据转化为CAD 数字模型，进行网格划分后的结果如图5 所示。

图5 煤体数字模型的网格划分Fig.5 Meshing of digital coal model

图6 边界条件设置示意图Fig.6 Schematic diagram of boundary conditions

3.2 结果分析

对建立好的模型进行解算求解，得到A、B、C、D 4个区域在X、Y、Z 3个方向的煤体孔隙压力分布，以及流场的速度分布。为节省篇幅，只列出了A、B、C、D 4个区域在X 方向上的孔隙压力分布和在X、Y 方向上的流场的速度分布，如图7 所示。

图7 煤体孔隙压力云图与瓦斯流线图Fig.7 Pressure nephograms of coal pore and gas streamline

在模拟结果的后处理时，为了方便对比，设置了统一的压力和速度标尺。

通过达西方程计算其渗透系数。

式中：Q为体积流量（m3·s-1）；A为流过面积（m2）；K为渗透系数（m·s-1）；Δh为水头差（m）；Δl为两端距离（m）。其中出入口边界质量流量可通过后处理的数据提取直接获得，结合甲烷密度可得体积流量。将质量流量与甲烷密度代入式（1）则得到式（2）。

式中：G为质量流量（kg·s-1）；g为重力加速度（m·s-2）；ΔP为两端压力差（Pa）。

根据建立的有限元模型，在渗透系数的每次计算中，g、Δl和ΔP 的值不变。渗透系数取决于G与A 的值，计算得到的渗透系数及其他参数的值如表2 所示。

根据模拟结果，沿着瓦斯流动的方向，煤体孔隙受到的压力逐渐降低，而瓦斯流动速度逐渐变大。当孔隙通道突然收缩时，瓦斯流动速度突然升高，孔隙压力则先降低后升高。孔隙压力的剧变体现在颜色上则是其发生的显著性变化，为了能够直观地观察到这种变化，可更改颜色标尺的范围，图7 中A 区域的局部放大图即是更改颜色标尺后的孔隙压力图。

对于其他孔隙通道变化平缓区域，压力变化相对较小，颜色变化不明显。正是由于压力在孔隙通道的急剧变化，从而导致了孔隙压力的层状分布。

表2 煤体不同渗流方向的渗透系数Table 2 Hydraulic conductivities of coal at different seepage directions

渗透系数作为煤体渗流的宏观参数，在微观尺度上，结构的差异会使得渗透系数相差较大，但其差异仍然保持在同一数量级。渗透系数也可以通过流线的丰满程度较直观地体现出来，相比A、B、D 3个模型在Y 方向的流线，C 模型的Y 方向流线更加丰满，而其渗透系数的值也明显高于其他3个模型。

4个区域的煤体在不同渗流方向上渗透系数的差异表明，微观尺度下的煤体渗透系数呈现各向异性。

4 结论与展望

（1）拥有复杂细节的精细煤体孔隙模型提高了后期有限元模型网格划分的难度，孔隙结构的精细表征与数值模拟的精确性之间存在着矛盾。二者平衡后的结果就是在尽量保持真实孔隙结构的前提下，对模型进行优化，将数值模拟的计算误差控制在可接受的范围内。正是由于允许误差的存在，逆向工程优化才变得可行。

（2）结合逆向工程技术，可以将煤体CT 三维数据转化煤体CAD 数字模型，并直接获取模型的孔隙体积、表面积以及孔隙度。转换得到的模型可被多种有限元模拟软件支持，间接拓宽了煤体CT三维数据的应用领域。CAD 数字模型的存在，使得网格划分过程可控，不再依赖专用的CT 数据处理软件，划分网格在精度上得以提高，CFD 模拟结果更加真实。

（3）建立的煤体有限元模型能够简化表征煤体微观孔隙结构，形象直观地描述瓦斯渗流速度场和压力场的空间分布，并且可被用于计算煤体的渗透系数。煤体孔隙介质和流体的属性共同决定了渗透系数的大小，煤体实际的渗透系数有待于进一步的试验验证。

（4）本文仅仅利用煤体有限元模型进行了渗流场的定性分析，并且假定渗流规律服从达西定律，如何从定量角度分析煤体内部渗流场及其流动规律有待于进一步研究。EDEM 等新型离散元软件也提供了对CAD 三维模型的支持，这使得由CT 三维数据转换得到的CAD 数字模型也可以进行离散元分析。CAD 数字模型可通过3D 打印技术制作成物理模型。对于使用3D 打印生成的物理模型进行渗流试验以及通过研究3D 打印材料的力学特性，基于相似材料模拟的方法，对3D 打印生成的物理模型进行加卸载试验等都值得进一步探索。

[1]LI H Y,OGAWA Y,SHIMADA S.Mechanism of methane flow through sheared coals and its role on methane recovery[J].Fuel,2003,82:1271-1279.

[2]降文萍,宋孝忠,钟玲文.基于低温液氮试验的不同煤体结构煤的孔隙特征及其对瓦斯突出影响[J].煤炭学报,2011,36(4):609-614.JIANG Wen-ping,SONG Xiao-zhong,ZHONG Ling-wen.Research on the pore properties of different coal body structure coals and the effects on gas outburst based on the low-temperature nitrogen adsorption method[J].Journal of China Coal Society,2011,36(4):609-614.

[3]张慧.煤孔隙的成因类型及其研究[J].煤炭学报,2001,26(1):40-44.ZHANG Hui.Genetical type of pores in coal reservoir and its research significance[J].Journal of China Coal Society,2001,26(1):40-44.

[4]宫伟力,安里千,赵海燕,等.基于图像描述的煤岩裂隙CT 图像多尺度特征[J].岩土力学,2010,31(2):371-376,381.GONG Wei-li,AN Li-qian,ZHAO Hai-yan,et al.Multiple scale characterization of CT image for coal rock fractures based on image description[J].Rock and Soil Mechanics,2010,31(2):371-376,381.

[5]MAZUMDER S,WOLF K,ELEWAUT K,et al.Application of X-ray computed tomography for analyzing cleat spacing and cleat aperture in coal samples[J].International Journal of Coal Geology,2006,68:205-222.

[6]宋晓夏,唐跃刚,李伟,等.基于显微的构造煤渗流孔精细表征[J].煤炭学报,2013,38(3):435-440.SONG Xiao-xia,TANG Yue-gang,LI Wei,et al.Advanced characterization of seepage pores in deformed coals based on micro-CT[J].Journal of China Coal Society,2013,38(3):435-440.

[7]GOLAB A,WARD C R,PERMANA A,et al.High-resolution three-dimensional imaging of coal using microfocus X-ray computed tomography,with special reference to modes of mineral occurrence[J].International Journal of Coal Geology,2013,113:97-108.

[8]SIMONS F J,VERHELST F,SWENNEN R.Quantitative characterization of coal by means of micro focal X-ray computed microtomography(CMT)and color image analysis(CIA)[J].International Journal of Coal Geology,1997,34:69-88.

[9]COLES M E,HAZLETT R D,SPANNE P,et al.Pore level imaging of fluid transport using synchrotron X-ray microtomography[J].Journal of Petroleum Science and Engineering,1998,19:55-63.

[10]VAN G M,SWENNEN R,DAVID P.Quantitative coal characterisation by means of microfocus X-ray computer tomography,colour image analysis and back-scattered scanning electron microscopy[J].International Journal of Coal Geology,2001,46:11-25.

[11]赵阳升,孟巧荣,康天合,等.显微CT 试验技术与花岗岩热破裂特征的细观研究[J].岩石力学与工程学报,2008,27(1):28-34.ZHAO Yang-sheng,MENG Qiao-rong,KANG Tian-he,el al.Micro-CT experimental technology and mesoinvestigation on thermal fracturing characteristics of granite[J].Chinese Journal of Rock Mechanics and Engineering,2008,27(1):28-34.

[12]NIE B S,HE X Q,LI X C.Meso-structures evolution rules of coal fracture with the computerized tomography scanning method[J].Engineering Failure Analysis,2014,41:81-88.

[13]LI S,TANG D Z,XU H,et al.Advanced characterization of physical properties of coals with different coal structures by nuclear magnetic resonance and X-ray computed tomography[J].Computers &Geosciences,2012,48:220-227.

[14]KARACAN C O,OKANDAN E.Adsorption and gas transport in coal microstructure:Investigation and evaluation by quantitative X-ray CT imaging[J].Fuel,2001,80:509-520.

[15]HE M C,WANG C G,FENG J L,et al.Experimental investigations on gas desorption and transport in stressed coal under isothermal conditions[J].International Journal of Coal Geology,2010,83:377-386.

[16]MATHEWS J P,PONE J D N,MITCHELL G D,et al.High-resolution X-ray computed tomography observations of the thermal drying of lump-sized subbituminous coal[J].Fuel Processing Technology,2011,92:58-64.

[17]NAUDÉ G,HOFFMAN J,THERON S J,et al.The use of X-ray computed tomography in the characterisation of coal and associated char reductants[J].Minerals Engineering,2013,52:143-154.

[18]于艳梅,胡耀青,梁卫国,等.瘦煤热破裂规律显微CT 试验[J].煤炭学报,2010,35(10):1696-1700.YU Yan-mei,HU Yao-qing,LIANG Wei-guo,et al.Micro-CT experimental research of lean coal thermal cracking laws[J].Journal of China Coal Society,2010,35(10):1696-1700.

[19]CAI Y D,LIU D M,MATHEWS J P,et al.Permeability evolution in fractured coal Combining triaxial confinement with X-ray computed tomography,acoustic emission and ultrasonic techniques[J].International Journal of Coal Geology,2014,122:91-104.

[20]YAO Y B,LIU D M.Comparison of low-field NMR and mercury intrusion porosimetry in characterizing pore size distributions of coals[J].Fuel,2012,92:152-158.

[21]TAKHAR P S,MAIER D E,CAMPANELLA O H,et al.Hybrid mixture theory based moisture transport and stress development in corn kernels during drying:Validation and simulation results[J].Journal of Food Engineering,2011,106:275-282.

[22]刘向君,朱洪林,梁利喜.基于微CT 技术的砂岩数字岩石物理试验[J].地球物理学报,2014,57(4):1133-1140.LIU Xiang-jun,ZHU Hong-lin,LIANG Li-xi.Digital rock physics of sandstone based on micro-CT technology[J].Chinese Journal of Geophysics,2014,57(4):1133-1140.

[23]YANG B H,WU A X,WANG C L,et al.Three-dimensional simulation of pore scale fluid flow in granular ore media with realistic geometry[J].Transactions of Nonferrous Metals,2012,22:3081-3086.

[24]BODLA K K,GARIMELLA S V,MURTHY J Y.3D reconstruction and design of porous media from thin sections[J].International Journal of Heat and Mass Transfer,2014,73:250-264.

[25]GÁLVEZ A,IGLESIAS A.Particle swarm optimization for non-uniform rational B-spline surface reconstruction from clouds of 3D data points[J].Information Sciences,2012,192:174-192.

[26]GÁLVEZ A,IGLESIAS A.A new iterative mutually coupled hybrid GA-PSO approach for curve fitting in manufacturing[J].Applied Soft Computing,2013,13:1491-1504.

[27]WANG J,GU D X,YU Z Y,et al.A framework for 3D model reconstruction in reverse engineering[J].Computers &Industrial Engineering,2012,63:1189-1200.

[28]BÉNiÈRE R,SUBSOL G,GESQUIÈRE G,et al.A comprehensive process of reverse engineering from 3D meshes to CAD models[J].Computer-Aided Design,2013,45:1382-1394.

[29]刘中华,徐素国,胡耀青,等.钙芒硝盐岩溶蚀试验研究[J].岩石力学与工程学报,2010,29(S2):3616-3621.LIU Zhong-hua,XU Su-guo,HU Yao-qing,et al.Experimental study of corrosion of glauberite salt rock[J].Chinese Journal of Rock Mechanics and Engineering,2010,29(S2):3616-3621.

[30]肖晓春,徐军,潘一山,等.功率超声影响的煤中甲烷气促解规律试验研究[J].岩石力学与工程学报,2013,32(1):65-71.XIAO Xiao-chun,XU Jun,PAN Yi-shan,et al.Experimental study of promoting methane desorption law in coal considering power ultrasound effect[J].Chinese Journal of Rock Mechanics and Engineering,2013,32(1):65-71.

[31]ZHU X L,AI S G,LU X F,et al.Collapse models of aluminum foam sandwiches under static three-point bending based on 3D geometrical reconstruction[J].Computational Materials Science,2014,85:38-45.

[32]ZHU X L,AI S G,LU X F,et al.Thermal conductivity of closed-cell aluminum foam based on the 3D geometrical reconstruction[J].International Journal of Heat and Mass Transfer,2014,72:242-249.

[33]YANG B H,WU A X,MIAO X X,et al.3D characterization and analysis of pore structure of packed ore particle beds based on computed tomography images[J].Transactions of Nonferrous Metals,2014,24:833-838.

[34]YOO D J.Three-dimensional surface reconstruction of human bone using a B-spline based interpolation approach[J].Computer-Aided Design,2011,43:934-947.

[35]李长圣,张丹,王宏宪,等.基于CT 扫描的土石混合体三维数值网格的建立[J].岩土力学,2014,35(9):2731-2736.LI Chang-sheng,ZHANG Dan,WANG Hong-xian,et al.3D mesh generation for soil-rock mixture based on CT scanning[J].Rock and Soil Mechanics,2014,35(9):2731-2736.