引用本文
孙英峰, 赵毅鑫, 王欣, 彭磊, 孙强. 基于同步辐射装置定量表征煤孔隙结构非均质性和各向异性[J]. 石油勘探与开发, 2019,46(6): 1128-1137.
SUN Yingfeng, ZHAO Yixin, WANG Xin, PENG Lei, SUN Qiang. Synchrotron radiation facility-based quantitative evaluation of pore structure heterogeneity and anisotropy in coal[J]. Petroleum Exploration & Development, 2019,46(6): 1128-1137.  
Doi: 10.11698/PED.2019.06.10
Permissions
Copyright©2019, 石油勘探与开发》编辑部
石油勘探与开发》编辑部 所有
基于同步辐射装置定量表征煤孔隙结构非均质性和各向异性
孙英峰1,2,3, 赵毅鑫1,3, 王欣4, 彭磊1,3, 孙强5
1. 中国矿业大学（北京）能源与矿业学院,北京 100083
2. 中国矿业大学（北京）应急管理与安全工程学院,北京 100083
3. 共伴生能源精准开采北京市重点实验室,北京 100083
4. 中国矿业大学（北京）力学与建筑工程学院,北京 100083
5. 中国石油勘探开发研究院,北京 100083
联系作者简介：赵毅鑫（1977-）,男,河南洛阳人,博士,中国矿业大学（北京）能源与矿业学院教授,主要从事煤储集层微观结构表征及其力学性质方面的研究。地址：北京市海淀区学院路丁11号,中国矿业大学（北京）能源与矿业学院,邮政编码：100083。E-mail：zhaoyx@cumtb.edu.cn

第一作者简介：孙英峰（1989-）,男,山西晋中人,博士,中国矿业大学（北京）在站博士后,主要从事煤储集层微观结构表征及煤层气开发方面的研究。地址：北京市海淀区学院路丁11号,中国矿业大学（北京）能源与矿业学院,邮政编码：100083。E-mail：yingfengsun@163.com

收稿日期: 2019-04-25
基金项目: 国家自然科学基金（51861145403,51874312）; 中国博士后科学基金（2018M641526）
摘要

为了定量表征煤孔隙结构非均质性和各向异性,利用同步辐射SAXS（小角X射线散射）获得两种不同变质程度煤样的SAXS图像,通过图像数据处理,得到表面分形（ D1）和孔分形（ D2）,由孔分形（ D2）定量表征两种煤样孔隙结构的非均质性,忻州窑煤的孔分形维数为2.74,唐山煤的孔分形维数为1.69,表明忻州窑煤的孔隙结构非均质性要比唐山煤强。利用同步辐射nano-CT（纳米CT）获得煤的三维孔隙结构成像,选取煤样兴趣区域（ROI）并分成一定数量的子块,通过计算子块孔隙度的相对标准偏差的极值来定量表征孔隙结构的非均质性,忻州窑煤的孔隙结构非均质性值为3.21,唐山煤的孔隙结构非均质性值为2.71,表明忻州窑煤的孔隙结构非均质性同样强于唐山煤,即同步辐射SAXS和同步辐射nano-CT获得的孔隙结构非均质性一致。考虑到孔隙结构各向异性和渗流能力各向异性之间的对应关系,孔隙结构各向异性的定量表征通过运用LBM（格子Boltzmann方法）方法计算孔隙结构的渗透率张量来实现,应用渗透率张量特征值和特征向量来表征孔隙结构的各向异性,该方法得到的孔隙结构各向异性得到三维孔隙结构几何形貌的验证。图11表1参48

关键词: 同步辐射SAXS; 同步辐射nano-CT; 煤; 孔隙结构; 非均质性; 各向异性
中图分类号:TE357 文献标志码:A 文章编号:1000-0747(2019)06-1128-10
Synchrotron radiation facility-based quantitative evaluation of pore structure heterogeneity and anisotropy in coal
SUN Yingfeng1,2,3, ZHAO Yixin1,3, WANG Xin4, PENG Lei1,3, SUN Qiang5
1. School of Energy and Mining Engineering, China University of Mining and Technology (Beijing), Beijing 100083, China
2. School of Emergency Management and Safety Engineering, China University of Mining and Technology (Beijing), Beijing 100083, China
3. Beijing Key Laboratory for Precise Mining of Intergrown Energy and Resources, China University of Mining and Technology (Beijing), Beijing 100083, China
4. School of Mechanics and Civil Engineering, China University of Mining and Technology (Beijing), Beijing 100083, China
5. Research Institute of Petroleum Exploration & Development, PetroChina, Beijing 100083, China;
Abstract

In order to quantify coal pore structure heterogeneity and anisotropy, synchrotron radiation SAXS (Small Angle X-ray Scattering) was applied to obtain the SAXS images of two different rank coal samples. The surface fractal dimension ( D1) and pore fractal dimension ( D2) were obtained by processing the image data. The pore structure heterogeneity of two coal samples was quantified by pore fractal dimension ( D2). Pore fractal dimension of Xinzhouyao coal is 2.74 and pore fractal dimension of Tangshan coal is 1.69. As a result, the pore structure heterogeneity of Xinzhouyao coal is stronger than that of Tangshan coal. 3D pore structure imaging was achieved by synchrotron radiation nano-CT. The selected Region of Interest (ROI) of coal sample was divided into a certain number of subvolumes. Pore structure heterogeneity was quantified by calculating the limit of the relative standard deviation of each subvolume's porosity. The heterogeneity value of Xinzhouyao coal pore structure is 3.21 and the heterogeneity value of Tangshan coal pore structure is 2.71. As a result, the pore structure heterogeneity of Xinzhouyao coal is also stronger than that of Tangshan coal, namely, pore structure heterogeneity from synchrotron radiation SAXS and synchrotron radiation nano-CT is consistent. Considering the corresponding relationship between the pore structure anisotropy and the permeability anisotropy, the quantification of pore structure anisotropy was realized by computing the permeability tensor of pore structure using the Lattice Boltzmann method (LBM), and the pore structure anisotropy was characterized by the eigenvalues and eigenvectors of the permeability tensor. The pore structure anisotropy obtained by the method proposed in this paper was validated by the pore structure geometrical morphology.

Keyword: synchrotron radiation SAXS; synchrotron radiation nano-CT; coal; pore structure; heterogeneity; anisotropy
0 引言

煤层气（CBM）已成为世界天然气资源的重要组成部分[1]。全球可采煤层气总量估计为（1.4~8.5）× 1014 m3[2]。煤是具有复杂孔隙结构的多孔介质[3], 目前定量表征煤孔隙结构的方法主要是压汞法[4, 5]和氮气/二氧化碳低压吸附法（LPA）[6]。压汞法由于高压注入可能破坏煤的孔隙结构或引起孔隙结构的变形[4]。LPA测孔的实验结果需要选用合适的分析模型（如Barrett-Joyner- Halenda（BJH）、Dubinin-Astakhov（D-A）、Horvath- Kawazoe（H-K）模型等）进行数据解释[7, 8, 9, 10, 11], 但是每个模型都有其适用的孔隙形状和孔径范围, 而且煤的孔隙结构具有复杂的几何形貌和宽泛的孔径范围, 所以很难找到一个单一的模型来准确分析煤中完整的孔径分布（PSD）, 因此, 有必要探索新的定量表征煤孔隙结构的方法。

2010年以来, 利用小角散射（SAXS）和nano-CT表征煤孔隙结构受到越来越多的关注。新一代实验室nano-CT可以实现纳米级的分辨率[12]。然而, 实验室nano-CT没有单色光束, 也不能实现高X射线相干, 这可能影响成像质量, 如射束硬化[13]。由于常规SAXS信号强度较弱, 测量时间就要拉长, 给实验带来很大的困难。测量一个样品往往需要几个小时或几十个小时, 即使采用高功率X射线源, 也难以满足许多实验工作的需要。而同步辐射光源不但强度高且光斑尺寸好, 可以使用很长的准直系统（可达10 m）, 从而大大提高实验的分辨率和灵敏度, 大幅度节约实验时间（一般样品曝光时间仅需几十秒）, 简化繁琐的数据修正工作。由于同步辐射光的波长连续可调, 在实验中可选择适当的波长来消除SAXS中较难解决的多次散射问题[14]。因此, 同步辐射SAXS与常规SAXS相比, 具有突出的优点。

国际上应用同步辐射SAXS对煤的孔隙结构开展了一些研究, Okolo G N等[15]采用低温氮吸附、压汞法和同步辐射SAXS 3种技术研究了4种煤样, 指出SAXS能探测的孔径范围更广, 且包括闭孔和开孔, 因此其测定样品的比表面积和孔隙度最大。Radlinski A P等[16]借助SAXS和SANS（小角中子散射）测量了7组原煤和相应粉末状样品, 研究认为球团模型充分代表了煤样的微观孔隙结构。赵毅鑫等[17]采用TEM（透射电子显微镜）和同步辐射SAXS测量了6种不同变质程度的煤样, 发现孔隙孔径分布与煤阶无关, 但随镜质组含量的变化而变化。Luo L等[18]采用同步辐射SAXS研究了不同变质程度的超细煤粉, 得出孔隙表面分形维数随煤阶的增加而减小。王博文等[19]采用同步辐射SAXS研究表明, 低挥发分烟煤脱灰分后, 其孔隙总数和体积增加, 微孔占比增大, 孔表面分形维数和孔比表面积增大, 但平均孔径和最可几孔径减小。作者团队2018年以来, 应用同步辐射nano-CT研究表征煤孔隙结构的实验方法和图像处理方法, 探索基于连通孔隙的CFD（计算流体动力学）数值模拟方法, CFD数值模拟得到的渗透率与其他学者得到的煤的渗透率在同一数量级[20]。

定量表征煤孔隙结构非均质性和各向异性对于揭示气体在煤层中的赋存和运移规律具有重要意义[21, 22, 23, 24, 25, 26]。根据文献研究, 对于孔隙结构非均质性和各向异性的定量表征方法已经开展了大量的研究[27, 28, 29, 30, 31, 32, 33, 34, 35, 36], 但是鲜有基于同步辐射装置定量表征煤孔隙结构非均质性和各向异性的研究。同步辐射SAXS和同步辐射nano-CT都具有纳米级的分辨率, 同步辐射SAXS孔径分辨率可以达到1 nm, 同步辐射nano-CT的空间分辨率可以达到30 nm, 而且对孔隙几何形貌没有要求, 这为定量研究煤的孔隙结构提供了可能性。因此, 本文利用同步辐射SAXS和同步辐射nano-CT, 进行定量表征煤孔隙结构非均质性和各向异性研究, 以期探索出新的定量表征煤孔隙结构非均质性和各向异性的方法。

1 样品和实验方法
1.1 煤样

实验中两个不同变质程度煤样分别采自山西省忻州窑煤矿11＃煤层和河北省唐山煤矿的9＃煤层。煤的工业分析[37]结果见表1。设备采用MDMDY-300全自动密度仪, 通过氦气置换法来测量两个煤样的真密度, 忻州窑煤的真密度为1.350 2 g/cm3, 唐山煤的真密度为1.535 5 g/cm3。根据中国石油天然气行业标准 SY/T 5124— 2012[38], 使用Leitz MPV-3显微光度计对抛光样品进行镜质组反射率测量和显微组分定量分析。从镜质组平均最大反射率来看, 忻州窑煤为气煤, 唐山煤为肥煤, 忻州窑煤的镜质组含量更高, 两个煤样的惰质组含量非常接近。根据中国石油天然气行业标准SY/T 5163— 2010[39], 应用D/MAX 2500 X射线衍射仪对两组煤样进行了X射线衍射（XRD）分析。XRD分析结果显示, 唐山煤的黏土矿物含量更高（见表1）。

表1 煤样成分分析
1.2 同步辐射SAXS实验

由于常规的孔隙结构表征方法难以探测到煤中的闭孔, 为了全面（包括闭孔和开孔）表征煤中孔隙结构的非均质性, 本文探索利用同步辐射SAXS定量表征煤中孔隙结构非均质性的方法, 同步辐射SAXS实验利用北京同步辐射装置的1W2A小角散射站开展。

同步辐射SAXS光路示意图如图1所示。角分辨率为0.5× 10-3rad, 能量分辨率约为1× 10-3, 样品处光强可达到1× 1011cps, 比常规X光机的强度高1× 103~1× 104倍, 入射X射线波长为0.154 nm, 光斑尺寸为1.4 mm× 0.2 mm。在同步辐射SAXS实验中, 实验样品选用粒径为0.18~0.25 mm（60~80目）的煤颗粒, 用3M胶带将样品安装在样品仓上, 将样品仓安装至综合平台上, 启动同步辐射SAXS装置, 设置曝光时间和曝光次数, 记录光电二极管读数, 保存散射图像数据。

图1 同步辐射SAXS光路示意图

当X射线照射到试样上, 如果试样内部存在纳米尺寸（2~100 nm）的密度不均匀区时, 则会在入射X射线的周围2° ~5° 的小角度范围内出现散射X射线。一个电子在不同方向的散射强度由汤姆逊（Thomson）公式决定[40]：

${{I}_{e}}\left( \theta \right)\text{=}{{I}_{0}}\frac{{{e}^{4}}}{{{m}^{2}}{{c}^{4}}}\frac{1}{{{a}^{2}}}\frac{1+{{\cos }^{2}}2\theta }{2}$ （1）

$I\left( q \right)\text{=}4\pi {{\left( {{\rho }_{m}}-{{\rho }_{p}} \right)}^{\text{2}}}{{\phi }_{s}}\left( 1-{{\phi }_{s}} \right)V\int_{0}^{\infty }{{{r}^{2}}{{\gamma }_{0}}\left( r \right)\frac{\sin \left( qr \right)}{qr}}\text{d}r$ （2）

对于多粒子体系, 如果粒子之间的相互距离远远大于粒子本身的尺寸, 散射强度则可近似为单个粒子散射强度的简单加和。因此, 对于M个不相干涉的粒子体系, 其散射强度可认为近似满足Guinier定律[41]：

$I\left( q \right)={{I}_{e}}M{{n}^{2}}\exp \left( -\frac{{{q}^{2}}R_{g}^{2}}{3} \right)$ （3）

${{R}_{g}}\text{=}\sqrt{\frac{3}{5}}R$ （4）

$q=\frac{4\pi \sin \theta }{\lambda }$ （5）

由（3）— （5）式可得, 散射强度随散射角度的增加而减小。通过SAXS实验观测可得到q和I(q)的关系, 对其进行数据处理, 进而获得散射体相关结构信息, 如散射体的回转半径、粒径大小、形状、分形维数、界面层厚度等。

1.3 同步辐射nano-CT实验

为了对比验证同步辐射SAXS定量表征煤中孔隙结构非均质性方法的有效性, 同时考虑到同步辐射nano-CT在探测闭孔和获得孔隙结构三维几何形貌方面具有的优势, 本文探索利用同步辐射nano-CT定量表征煤中孔隙结构非均质性的方法, 同步辐射nano-CT实验利用北京同步辐射装置的4W1A-X射线成像实验站开展。

同步辐射nano-CT的能量范围为5~12 keV, 光斑尺寸为10 µ m× 10 µ m, 可以达到的空间分辨率优于30 nm, 对于该同步辐射nano-CT的详细介绍见文献[42]。对于同步辐射nano-CT成像, 样品颗粒的尺寸要求小于10 μ m, 因此, 首先粉碎煤样, 然后选择尺寸接近10 μ m的颗粒, 并通过图2c所示的仪器粘在大头针的针尖上。由于后续图像校准需要参考点, 通过图2d所示的仪器将直径为0.5 µ m的球形金颗粒安装在样品上, 用作后期图像校准的参考点。最后, 将粘有样品的大头针夹在如图2b所示的样品转台上。在同步辐射nano-CT成像过程中, 样品台以预设的角度步（旋转速度）不连续地旋转以获得不同角度的投影图像。

图2 同步辐射nano-CT成像用到的仪器照片

2 实验结果和讨论
2.1 煤孔隙结构非均质性定量表征

2.1.1 同步辐射SAXS定量表征煤孔隙结构非均质性

同步辐射SAXS实验所获得的两种不同变质程度煤样的SAXS散射图像如图3所示。

图3 煤的散射图像

为了方便数据分析, 通过FIT2D软件将图3所示的SAXS散射图像转为对应的散射曲线, 如图4所示。其中, 横坐标|q|为倒空间的散射矢量的模, 纵坐标为散射强度。

图4 煤的散射曲线

分形物体的X光小角散射强度在分形区可用幂律表示：

$I\left( q \right)={{I}_{0}}{{q}^{-\alpha }}$ （6）

α 为[0, 4], 取$\ln \left| I\left( q \right) \right|-\ln \left| q \right|$曲线线性区域的斜率为d, 则α =– d。当0< α < 3时, 散射体为质量分形或孔分形, 其分形维数${{D}_{2}}=\alpha $; 当3< α < 4时, 散射体为表面分形, 其分形维数${{D}_{1}}=6-\alpha $。因此, 若$\ln \left| I\left( q \right) \right|-\ln q$曲线存在线性区域, 则表明存在分形现象, 进而可通过α 大小判断分形为表面分形还是孔分形[43]。

图5为两个煤样的|q|和散射强度I对应的双对数曲线。通过在不同|q|值范围下, 对两个煤样的双对数曲线作切线的方法, 可得相应|q|值范围段煤样孔隙的分形维数。

图5 煤样分形求解对应的双对数曲线

在低|q|值区域（-3~-1 nm-1）曲线的分形特征为表面分形（D1）, 由此可分析孔隙表面的光滑度和平整度, 其分形维数描述孔隙表面光滑度和平整度, 其值越大, 则说明孔隙表面越不规则。高|q|值区域（-1~1 nm-1）曲线的分形特征表现为孔分形（D2）, 其分形维数可分析孔隙空间分布的非均质性, 其值越大, 说明孔隙空间分布的非均质性越强。由图5的表面分形（D1）可以看出, 忻州窑煤孔隙表面更不规则, 由图5的孔分形（D2）可以看出, 忻州窑煤孔隙空间分布的非均质性更强。

2.1.2 同步辐射nano-CT定量表征煤孔隙结构非均质性

同步辐射nano-CT获得煤的成像以后, 图像处理流程主要包括图像校准、三维重构、兴趣区域（ROI）选取、噪声滤除和图像分割, 如图6所示。图像处理具体方法见文献[20]。

图6 同步辐射nano-CT图像处理

图像分割以后, 对孔隙结构的连通性进行了分析。连通孔隙指ROI立方体从一个面到对面连续连通的孔隙, 两个煤样中的连通孔隙如图7a和7b所示。通过分析计算, 忻州窑煤连通孔隙的体积是1.71 µ m3, 占总孔隙体积2.47 µ m3的69%, 而唐山煤连通孔隙的体积是2.70 µ m3, 占总孔隙体积3.24 µ m3的83%, 这说明唐山煤有更多的连通孔隙。

图7 连通孔隙的几何形貌

煤样ROI中不同子块孔隙度之间的差异可以表征孔隙结构非均质性, 本文中煤样ROI尺寸均为200像素× 200像素× 200像素, 像素尺寸为0.014 59 μ m× 0.014 59 μ m× 0.014 59 μ m。煤样ROI中不同子块孔隙度之间的差异可以通过计算各个子块孔隙度的相对标准偏差来定量表征：

$H=\frac{1}{\phi }\sqrt{\frac{\sum{{{\left( {{\phi }_{i}}-\phi \right)}^{2}}}}{N-1}}$ （7）

如图8所示, 当煤样ROI在（X, Y, Z）方向上分别被等分成s段, 煤样ROI被分成s3个子块, 各个子块孔隙度的相对标准偏差H随分割段数s的变化而变化。所以有必要提出一个与子块数量无关的参数来表征煤的三维孔隙结构的非均质性。虽然曲线随着分割段数s的增加而上升, 但是曲线变得越来越平缓。为了获得一个常数来表征煤的三维孔隙结构的非均质性, 用下面的函数来拟合图8中的数据, 函数H（s）的极值是一个常数：

$H\left( s \right)=\frac{As}{1+Bs}$ （8）

$HV=\underset{x\to \infty }{\mathop{\lim }}\, \frac{As}{1+Bs}=\frac{A}{B}$ （9）

图8 子块孔隙度相对标准偏差与分割子块数量的关系

HV值可以表征煤样三维孔隙结构的非均质性, HV值越大, 孔隙结构的非均质性越强。从（9）式可以发现该方法得到的非均质性值消除了分割子块数量的影响。

忻州窑煤：A=0.210 47, B=0.065 62, R2=0.999 41, HV=3.21; 唐山煤：A=0.175 47, B=0.064 63, R2=0.994 02, HV=2.71。

如图8所示, 两个煤样的非均质性值可以很好地用（8）式拟合。计算结果表明, 忻州窑煤的HV值为3.21, 而唐山煤的HV值为2.71; 可见, 忻州窑煤比唐山煤孔隙结构表现出更高的非均质性, 也可以从图8中两条曲线的对比中得到相同的结论。

关于两个煤样孔隙结构的非均质性, 忻州窑煤孔隙结构的非均质性要比唐山煤强, 同步辐射SAXS和同步辐射nano-CT得到的结果一致, 这也验证了这两种表征煤孔隙结构非均质性方法的有效性。虽然孔隙结构非均质性分级标准的建立需要对大量煤样进行统计, 但是借助同步辐射SAXS和同步辐射nano-CT分别获得的孔分形（D2）和HV值为定量描述孔隙结构的非均质性提供了手段, 这两个指标可以实现对多个样品的孔隙结构非均质性进行排序, 也为孔隙结构非均质性分级标准（强-中-弱）的建立提供了可能。

根据文献研究, 煤孔隙结构的非均质性随镜质组含量的增加而增加, 随着惰质组含量的增加而降低[44], 与灰分含量呈现线性正相关关系[45]。忻州窑煤的镜质组含量高于唐山煤, 两个煤样的惰质组和灰分含量接近, 忻州窑煤孔隙结构的非均质性强于唐山煤, 这与前人的相关研究结果一致。

2.2 煤孔隙结构各向异性定量表征

考虑到孔隙结构各向异性和渗透能力各向异性之间的对应关系, 孔隙结构各向异性的定量表征可以通过运用LBM（格子Boltzmann方法）方法计算渗透率张量来实现, 孔隙结构的各向异性通过渗透率张量的特征值和特征向量来表征。渗透率张量的特征值反映了由相应特征向量确定的3个方向上的渗透率。通过特征值的对比分析可以找到孔隙结构渗透性最好的方向, 渗透率张量特征值中最大值对应的特征向量确定的方向上渗透性最好, 而且特征值中的最大值也代表了孔隙结构可以达到的最佳的渗透率, 为对比不同样品的渗透性提供了参考。

本文中用到的LBM方法采用单松弛时间Bhatnagar- Gross-Krook（BGK）模型[46], 本模型通过求解一个离散的Boltzmann方程来模拟格子节点中迁移和碰撞的流体颗粒分布的变化, 然后宏观流体流动用虚拟流体颗粒的整体运动来近似。在方向${{e}_{i}}$上的流体颗粒分布函数${{f}_{i}}$在每一个时间步通过以下的方程来更新：

${{f}_{i}}\left( x+{{e}_{i}}, t+\Delta t \right)-{{f}_{i}}\left( x, t \right)=-\frac{1}{\tau }\left[ {{f}_{i}}\left( x, t \right)-{{f}_{\text{eq, }i}}\left( x, t \right) \right]$ （10）

式中${{f}_{\text{eq}}}$是截面平衡分布, 定义为：

${{f}_{\text{eq}, i}}={{w}_{i}}\rho \left[ 1+\frac{3{{e}_{i}}\cdot v}{c_{\text{s}}^{2}}+\frac{9{{\left( {{e}_{i}}\cdot v \right)}^{2}}}{c_{\text{s}}^{4}}-\frac{v\cdot v}{2c_{\text{s}}^{2}} \right]$ （11）

$\tau $是与流体动力黏度$\upsilon $相关的参数, 在Chapman- Enskog方程中表示为：

$\upsilon =\Delta tc_{s}^{2}\left( \tau -\frac{1}{2} \right)$ （12）

$\rho $和$v$是宏观的密度和速度, 它们由颗粒分布的变化决定：

$\rho =\sum\limits_{i=1}^{N}{{{f}_{i}}}$ （13）

$v=\frac{1}{\rho }\sum\limits_{i=1}^{N}{{{f}_{i}}{{e}_{i}}}$ （14）

根据达西定律, 渗透率张量可以由下式计算[47]：

${{k}_{ij}}=-\frac{{{\mu }_{v}}}{{{({{\nabla }_{x}}p)}_{j}}}\frac{1}{{{V}_{\Omega }}}\int_{\Omega }{{{v}_{i}}\left( x \right)d\Omega }$ （15）

通过3个独立的模拟, 在3个流动方向（X, Y, Z）的两个相对面上施加压力梯度并在剩余4个面上施加无流动边界条件, 计算出渗透率张量的3个对角值。

在渗透率张量是对称的和正定的假设下, 达西定律可以扩展为：

$\left[ \begin{matrix} \frac{\partial p}{\partial {{x}_{2}}} & \frac{\partial p}{\partial {{x}_{3}}} & 0 \\ \frac{\partial p}{\partial {{x}_{1}}} & 0 & \frac{\partial p}{\partial {{x}_{3}}} \\ 0 & \frac{\partial p}{\partial {{x}_{1}}} & \frac{\partial p}{\partial {{x}_{2}}} \\ \end{matrix} \right]\left[ \begin{matrix} {{k}_{12}} \\ {{k}_{13}} \\ {{k}_{23}} \\ \end{matrix} \right]=\left[ \begin{matrix} -{{\mu }_{v}}{{v}_{1}}-{{k}_{11}}\frac{\partial p}{\partial {{x}_{1}}} \\ -{{\mu }_{v}}{{v}_{2}}-{{k}_{22}}\frac{\partial p}{\partial {{x}_{2}}} \\ -{{\mu }_{v}}{{v}_{3}}-{{k}_{33}}\frac{\partial p}{\partial {{x}_{3}}} \\ \end{matrix} \right]$ （16）

如（16）式所示, 如果得到${{v}_{1}}$、${{v}_{2}}$、${{v}_{3}}$, 渗透率张量的非对角值可以由渗透率张量已经得到的对角值计算。

为了获得${{v}_{1}}$、${{v}_{2}}$、${{v}_{3}}$, 进行了另一组模拟, 在这组模拟中, 压力梯度同时施加在3个正交方向上。图9展示了4组LBM模拟的流速流线图, 可以发现流线分布与压力梯度的方向一致。

图9 渗透率张量计算得到的流速流线图

为了分析三维孔隙结构的各向异性, 计算了两个煤样ROI的渗透率张量的特征值和特征向量：忻州窑煤的渗透率张量计算结果是：

$\left[ \begin{matrix} 0.68 & 2.02 & 4.23 \\ {} & 0.87 & 0.40 \\ {} & {} & 0.11 \\ \end{matrix} \right]\times {{10}^{-6}}\ \text{ }\!\!\mu\!\!\text{ }{{\text{m}}^{2}}$

特征值是k1=-4.10× 10-6 μ m2, k2=0.42× 10-6μ m2, k3=5.33× 10-6 μ m2, 相应特征向量是：

${{E}_{1}}=\left[ \begin{matrix} -0.70 \\ 0.23 \\ 0.68 \\ \end{matrix} \right]\text{ }{{E}_{2}}=\left[ \begin{matrix} 0.12 \\ -0.90 \\ 0.42 \\ \end{matrix} \right]\text{ }{{E}_{3}}=\left[ \begin{matrix} 0.71 \\ 0.37 \\ 0.60 \\ \end{matrix} \right]$

唐山煤的渗透率张量计算结果是：

$\left[ \begin{matrix} 0.93 & 6.10 & -0.57 \\ {} & 0.69 & -0.21 \\ {} & {} & 0.57 \\ \end{matrix} \right]\times {{10}^{-6}}\ \text{ }\!\!\mu\!\!\text{ }{{\text{m}}^{2}}$

采用同样的方法, 计算出渗透率张量的特征值是k1=-5.30× 10-6 μ m2, k2=0.54× 10-6μ m2, k3=6.69× 10-6 μ m2, 相应的特征向量是：

${{E}_{1}}=\left[ \begin{matrix} 0.70 \\ -0.71 \\ 0.04 \\\end{matrix} \right]\text{ }{{E}_{2}}=\left[ \begin{matrix} 0.03 \\ 0.09 \\ 1.00 \\\end{matrix} \right]\text{ }{{E}_{3}}=\left[ \begin{matrix} 0.71 \\ 0.70 \\ -0.09 \\\end{matrix} \right]$

可以看出唐山煤的渗透率张量的特征值比忻州窑煤对应的特征值大, 这和唐山煤有更多的连通孔隙的相符合。

孔隙结构各向异性分析的目的是找到渗透性最好的方向, 但是值得注意的是, 两组煤样的渗透率张量特征值都有负值, 表明在相应特征向量所确定的方向上存在U形连通孔隙。渗透率张量特征值中最大的值对应的特征向量所确定的方向上渗透性最好, 在该方向上连通孔隙的截面积也最大, 渗透率张量特征值中第2大的值对应的特征向量所确定的方向上渗透性第2好, 在该方向上连通孔隙的截面积也第2大, 在图7a和图7b中, 相应特征向量所确定的方向被标出, 发现渗透率的各项异性和孔隙结构几何相貌的各向异性之间有明显的相关性。特征值比较大的方向, 连通孔隙的横截面积也比较大。

2.3 表征煤孔隙结构不同方法的对比分析

为了体现同步辐射SAXS和同步辐射nano-CT方法的特点和优势, 本文将同步辐射nano-CT和同步辐射SAXS测得的孔径分布与低温氮气吸附解吸和核磁共振冻融测得的结果进行了对比（见图10）。低温氮气吸附解吸和核磁共振冻融实验和数据处理方法在作者之前发表的论文中做了详细介绍[48]。从图10a可以发现, 孔径小于100 nm时, 低温氮气吸附解吸和核磁共振冻融测得的孔体积比同步辐射nano-CT要高, 这可以由同步辐射nano-CT分辨率的局限性解释, 一些孔径小于同步辐射nano-CT分辨率的孔隙不能被同步辐射nano-CT探测到, 但是可以被低温氮气吸附解吸和核磁共振冻融探测到, 而当孔径大于100 nm时, 同步辐射nano-CT测得的孔体积比其他两种方法要大, 这是由于低温氮气吸附解吸和核磁共振冻融探测不到闭孔, 而同步辐射nano-CT可以探测到闭孔所导致的。同步辐射nano-CT在探测闭孔方面具有优势, 但是由于分辨率的限制, 在探测孔径较小的孔隙方面存在不足。

图10 不同方法测得的煤孔径分布对比

从图10b可以发现, 同步辐射SAXS和低温氮气吸附解吸的测试结果具有很好的一致性, 在0~40 nm和40~80 nm孔径范围内分别具有两个峰。在20~25 nm孔径范围内, 同步辐射SAXS测得的孔隙体积要大于低温氮气吸附解吸, 这可以由同步辐射SAXS可以测得闭孔, 而低温氮气吸附解吸测不到闭孔来解释。在其他范围内同步辐射SAXS测得的孔体积小于低温氮气吸附解吸, 出现这种现象的原因还有待进一步研究。在笔者开展的应用同步辐射SAXS和核磁共振冻融测试页岩的孔隙结构中, 也出现在部分区域同步辐射SAXS测得的孔体积大于核磁共振冻融, 而在其他区域小于核磁共振冻融的情况[41]。虽然同步辐射SAXS测得的孔径范围不能覆盖全部的孔径范围, 但是如图11所示, 根据核磁共振冻融测试的结果, 孔径小于80 nm的孔体积占孔径500 nm以下的孔体积的74%; 所以, 同步辐射SAXS所得的孔径范围对于分析煤的孔隙结构还是具有一定意义的。

图11 核磁共振冻融测得的累计孔体积

总之, 与常规的孔隙结构表征方法相比, 虽然同步辐射SAXS可探测的孔径范围较小, 同步辐射nano-CT难以探测到一些孔径小于其分辨率的孔隙, 但二者的优势在于可以探测到闭孔, 且同步辐射nano-CT可以获得孔隙结构的三维几何形貌。另外, 两种方法的结合使用可以互相弥补各自在探测孔径较大和较小的孔隙时存在的劣势。

3 结论

针对定量表征煤孔隙结构非均质性和各向异性的研究存在的不足, 探索了利用同步辐射SAXS和同步辐射nano-CT定量表征煤孔隙结构非均质性和各向异性的方法。在利用同步辐射SAXS定量表征煤孔隙结构非均质性中, 由孔分形D2定量表征两个煤样孔隙结构的非均质性, 忻州窑煤的孔分形维数为2.74, 唐山煤的孔分形维数为1.69。在利用同步辐射nano-CT定量表征煤孔隙结构非均质性中, 把煤样ROI分成一定数量的子块, 通过计算各个子块孔隙度的相对标准偏差的极值来定量表征孔隙结构的非均质性, 忻州窑煤的孔隙结构非均质性值为3.21, 唐山煤的孔隙结构非均质性值为2.71。两种手段获得的煤孔隙结构的非均质性具有一致性, 即忻州窑煤孔隙结构的非均质性强于唐山煤。在利用同步辐射nano-CT定量表征煤孔隙结构各向异性中, 考虑到孔隙结构各向异性和渗透能力各向异性之间的对应关系, 孔隙结构各向异性的定量表征通过运用LBM方法计算孔隙结构的渗透率张量来实现, 利用渗透率张量的特征值和特征向量来定量表征孔隙结构的各向异性, 孔隙结构的形貌验证了本文提出的定量表征煤孔隙结构各向异性方法的有效性。

致谢：对中国科学院高能物理研究所李志宏和黄万霞在同步辐射SAXS和同步辐射nano-CT实验中给予的帮助、北京市计算中心李娇娇在数值计算中提供的帮助表示衷心感谢。

符号注释：

a— — 试样到接收器的距离, m; A— — 拟合参数, 无因次; B— — 拟合参数, 无因次; c— — 光速, 299 792 458 m/s; cs— — 声速, 346 m/s; d— — 斜率, 无因次; D— — 孔径, nm; D1— — 表面分形维数, 无因次; D2— — 质量或孔分形维数, 无因次; e— — 电子电荷, 1.6× 10-19 C; ei— — 方向向量; E1, E2, E3— — 特征向量, 无因次; fi— — 流体颗粒分布函数; feq— — 截面平衡分布函数; H— — 各个子块孔隙度的相对标准偏差, 无因次; HV— — 煤样三维孔隙结构的非均质性值, 无因次; ie— — 常量, 7.9× 10-20 cm2; Ie— — 单个电子对光子的散射强度; Ie（θ ）— — 单个电子在不同方向对光子的散射强度; I（q）— — 散射强度; I0— — 入射X射线光通量, phs/（mm2· s）; k1, k2, k3— — 渗透率张量的特征值, μ m2; kij— — 渗透率张量, μ m2; m— — 电子质量, 9.109 56× 10-31kg; M— — 不相干涉的粒子体系数量, 个; n— — 一个粒子中的电子数目, 个; N— — 子块数量, 个; p— — 流体压力, Pa; q— — 散射矢量, nm-1; r— — 颗粒半径, nm; R— — 粒子半径, nm; Rg— — 粒子的回转半径, nm; s— — 煤样ROI在（X, Y, Z）方向上分别被等分成的段数, 无因次; t— — 时间, s; Δ t— — 时间增量, s; 2θ — — 散射角; v— — 宏观速度, m/s; $v{}_{i}(x)$— — 速度分布函数; V— — SAXS辐照的体积, mm3; Vp— — 孔隙体积, cm3/g; VΩ — — 积分区域的体积, m3; wi— — 权系数, 无因次; x— — 流体颗粒的位置, m; α — — 常量, 无因次; γ 0— — 相关函数; λ — — 波长, nm; μ v— — 流体的运动黏度, m2/s; ρ — — 宏观密度, kg3/m3; ρ m— — 基质电子密度, eA-3; ρ p— — 微孔电子密度, eA-3; τ — — 松弛时间, s; υ — — 流体动力黏度, Pa· s; ϕ s— — 样品的孔隙度, %; ϕ i— — 第i个子块的孔隙度, %, ϕ — — 煤样整体的孔隙度, %; Ω — — 积分区域; x— — x方向梯度算子, Pa/m。

参考文献
[1] EIA. Annual Energy Outlook[R/OL]. (2013-1-21)[2019-4-20]. https://www.eia.gov/pressroom/presentations/sieminski_01212013.pdf. [本文引用:1]
[2] DONG Z, HOLDITCH S A, AYERS W B, et al. Probabilistic estimate of global coalbed methane recoverable resources[J]. SPE Economics & Management, 2015, 7(4): 1-9. [本文引用:1]
[3] VAN KREVELEN D W. Coal-typology, chemistry, physics, constitution[M]. Amsterdam: Elsevier Science Publishing Company Inc, 1993. [本文引用:1]
[4] MAHAMUD M M, NOVO M F. The use of fractal analysis in the textural characterization of coals[J]. Fuel, 2008, 87(2): 222-231. [本文引用:2]
[5] ROOTARE H M, PRENZLOW C F. Surface areas from mercury porosimeter measurements[J]. Journal of Physical Chemistry, 2002, 71(8): 2733-2736. [本文引用:1]
[6] RODRIGUES C F, SOUSA M J L D. The measurement of coal porosity with different gases[J]. International Journal of Coal Geology, 2002, 48(3/4): 245-251. [本文引用:1]
[7] BARRETT E P, JOYNER L G, HALENDA P P. The determination of pore volume and area distributions in porous substances[J]. Journal of the American Chemical Society, 1951, 73(1): 373-380. [本文引用:1]
[8] TERZYK A P, GAUDEN P A, KOWALCZYK P. What kind of pore size distribution is assumed in the Dubinin-Astakhov adsorption isotherm equation?[J]. Carbon, 2002, 40(15): 2879-2886. [本文引用:1]
[9] DOMBROWSKI R J, LASTOSKIE C M, HYDUKE D R. The Horvath-Kawazoe method revisited[J]. Colloids & Surfaces A: Physicochemical & Engineering Aspects, 2001, 187(1): 23-39. [本文引用:1]
[10] SAITO A, FOLEY H C. Curvature and parametric sensitivity in models for adsorption in micropores[J]. Aiche Journal, 2010, 37(3): 429-436. [本文引用:1]
[11] SHAO X H, ZHANG X R, WANG W C. Comparison of density functional theory and molecular simulation methods for pore size distribution of mesoporous materials[J]. Acta Physico-chimica Sinica, 2003, 19(6): 538-542. [本文引用:1]
[12] SAKDINAWAT A, ATTWOOD D. Nanoscale X-ray imaging[J]. Nature Photonics, 2009, 4(12): 840-848. [本文引用:1]
[13] BOSSA N, CHAURAND P, VICENTE J, et al. Micro- and nano-X-ray computed - tomography: A step forward in the characterization of the pore network of a leached cement paste[J]. Cement & Concrete Research, 2015, 67: 138-147. [本文引用:1]
[14] 李志宏, 吴忠华. 我国同步辐射小角X光散射装置[J]. 光散射学报, 2003, 15(3): 213-216.
LI Zhihong, WU Zhonghua. Facility of Chinese synchrotron radiation small angle X- Ray scattering[J]. Chinese Journal of Light Scattering, 2003, 15(3): 213-216. [本文引用:1]
[15] OKOLO G N, EVERSON R C, NEOMAGUS H W J P, et al. Comparing the porosity and surface areas of coal as measured by gas adsorption, mercury intrusion and SAXS techniques[J]. Fuel, 2015, 141: 293-304. [本文引用:1]
[16] RADLINSKI A P, MASTALERZ M, HINDE A L, et al. Application of SAXS and SANS in evaluation of porosity, pore size distribution and surface area of coal[J]. International Journal of Coal Geology, 2004, 59(3): 245-271. [本文引用:1]
[17] ZHAO Y, LIU S, ELSWORTH D, et al. Pore structure characterization of coal by synchrotron small-angle X-ray scattering and transmission electron microscopy[J]. Energy & Fuels, 2014, 28(6): 3704-3711. [本文引用:1]
[18] LUO L, LIU J, ZHANG Y, et al. Application of small angle X-ray scattering in evaluation of pore structure of superfine pulverized coal/char[J]. Fuel, 2016, 185: 190-198. [本文引用:1]
[19] 王博文, 李伟. 基于小角X射线散射的灰分对煤孔隙结构的影响研究[J]. 煤矿安全, 2017, 48(1): 144-148.
WANG Bowen, LI Wei. Effect of ash content on coal pore structure based on small angle X-ray scattering[J]. Safety in Coal Mines, 2017, 48(1): 144-148. [本文引用:1]
[20] ZHAO Y X, SUN Y F, LIU S M, et al. Pore structure characterization of coal by synchrotron radiation nano-CT[J]. Fuel, 2018, 215: 102-110. [本文引用:2]
[21] GAMSON P, BEAMISH B, JOHNSON D. Coal microstructure and secondary mineralization: Their effect on methane recovery[J]. Geological Society of London Special Publications, 1996, 109(1): 165-179. [本文引用:1]
[22] GAMSON P, BEAMISH B, JOHNSON D. Coal microstructure and micropermeability and their effects on natural gas recovery[J]. Fuel, 1993, 72(1): 87-99. [本文引用:1]
[23] HARPALANI S, ZHAO X. Microstructure of coal and its influence on flow of gas[J]. Energy Sources, 2007, 13(2): 229-242. [本文引用:1]
[24] TANG X, JIANG Z, LI Z, et al. The effect of the variation in material composition on the heterogeneous pore structure of high-maturity shale of the Silurian Longmaxi formation in the southeastern Sichuan Basin, China[J]. Journal of Natural Gas Science & Engineering, 2015, 23: 464-473. [本文引用:1]
[25] WAN Y, PAN Z, TANG S, et al. An experimental investigation of diffusivity and porosity anisotropy of a Chinese gas shale[J]. Journal of Natural Gas Science & Engineering, 2015, 23: 70-79. [本文引用:1]
[26] HU Q, GAO X, GAO Z, et al. Pore Accessibility and Connectivity of Mineral and Kerogen Phases for Shales[R]. Denver, USA: Proceedings of the Unconventional Resources Technology Conference, 2014. [本文引用:1]
[27] WANG Z, JIN X, WANG X, et al. Pore-scale geometry effects on gas permeability in shale[J]. Journal of Natural Gas Science & Engineering, 2016, 34: 948-957. [本文引用:1]
[28] KNOWLES J, ARMATAS G, HUDSON M, et al. Pore anisotropy and microporosity in nanostructured mesoporous solids[J]. Langmuir the Acs Journal of Surfaces & Colloids, 2006, 22(1): 410-418. [本文引用:1]
[29] POMONIS P J, ARMATAS G S. A method for the estimation of pore anisotropy in porous solids[J]. Langmuir the Acs Journal of Surfaces & Colloids, 2004, 20(16): 6719. [本文引用:1]
[30] BULTREYS T, BOEVER W D, CNUDDE V. Imaging and image-based fluid transport modeling at the pore scale in geological materials: A practical introduction to the current state-of-the-art[J]. Earth-Science Reviews, 2016, 155: 93-128. [本文引用:1]
[31] RAMSTAD T, OREN P E, BAKKE S. Simulation of two phase flow in reservoir rocks using a lattice boltzmann method[J]. SPE Journal, 2010, 15(4): 917-927. [本文引用:1]
[32] ICARDI M, BOCCARDO G, MARCHISIO D L, et al. Pore-scale simulation of fluid flow and solute dispersion in three-dimensional porous media[J]. Physical Review E, 2014, 90(1): 13032. [本文引用:1]
[33] KOROTEEV D, DINARIEV O, EVSEEV N, et al. Direct hydrodynamic simulation of multiphase flow in porous rock[J]. Petrophysics, 2014, 55(4): 294-303. [本文引用:1]
[34] ASTORINO M, SAGREDO J B, QUARTERONI A. A modular lattice boltzmann solver for GPU computing processors[J]. Sema Journal, 2012, 59(1): 53-78. [本文引用:1]
[35] BLUNT M J. Flow in porous media-pore-network models and multiphase flow[J]. Current Opinion in Colloid & Interface Science, 2001, 6(3): 197-207. [本文引用:1]
[36] YOUSSEF S, ROSENBERG E, GLAND N, et al. High Resolution CT and Pore-Network Models To Assess Petrophysical Properties Of Homogeneous and Heterogeneous Carbonates[R]. Abu Dhabi, UAE: SPE/EAGE Reservoir Characterization and Simulation Conference, 2007. [本文引用:1]
[37] Astm International. Stand ard test methods for proximate analysis of the analysis sample of coal and coke by instrumental procedures: ASTM D5142[S]. West Conshohocken: Astm International, 2009. [本文引用:1]
[38] 国家能源局. 沉积岩中镜质体反射率测定方法: SY/T 5124—2012 [S]. 北京: 石油工业出版社, 2012.
National Energy Administration. Method of determing microscopically the reflectance of vitrinite in sedimentary: SY/T 5124—2012[S]. Beijing: Petroleum Industry Press, 2012. [本文引用:1]
[39] 国家能源局. 沉积岩中黏土矿物和常见非黏土矿物X射线衍射分析方法: SY/T 5163—2010[S]. 北京: 石油工业出版社, 2010.
National Energy Administration. Analysis method for clay minerals and ordinary non-clay minerals in sedimentary rocks by the X-ray diffraction: SY/T 5163—2010[S]. Beijing: Petroleum Industry Press, 2010. [本文引用:1]
[40] 李志宏. SAXS方法及其在胶体和介孔材料研究中的应用[D]. 太原: 中国科学院山西煤炭化学研究所, 2002.
LI Zhihong. SAXS method and its application in colloid and mesoporous material[D]. Taiyuan: Institute of Coal Chemistry Chinese Academy of Science, 2002. [本文引用:1]
[41] ZHAO Y X, PENG L, LIU S M, et al. Pore structure characterization of shales using synchrotron SAXS and NMR cryoporometry[J]. Marine and Petroleum Geology, 2019, 102: 116-125. [本文引用:2]
[42] WANG S H, ZHANG K, WANG Z L, et al. A user-friendly nano-CT image alignment and 3D reconstruction platform based on LabVIEW[J]. Chinese Physics C, 2015, 39(1): 018001. [本文引用:1]
[43] WIJNEN P W J G, BEELEN T P M, RUMMENS K P J, et al. Silica gel from water glass: A SAXS study of the formation and ageing of fractal aggregates[J]. Journal of Applied Crystallography, 1991, 24(5): 759-764. [本文引用:1]
[44] 李振涛. 煤储层孔裂隙演化及对煤层气微观流动的影响[D]. 北京: 中国地质大学(北京), 2018.
LI Zhentao. Evolution of pore-fractures of coal reservoir and its impact on CBM microcosmic flow[D]. Beijing: China University of Geosciences (Beijing), 2018. [本文引用:1]
[45] 姚艳斌, 刘大锰. 煤储层精细定量表征与综合评价模型[M]. 北京: 地质出版社, 2013.
YAO Yanbin, LIU Dameng. Advanced quantitative characterization and comprehensive evaluation model of coalbed methane reservoirs[M]. Beijing: Geological Publishing House, 2013. [本文引用:1]
[46] SUCCI S. The Lattice Boltzmann Equation[M]. Oxford: Oxford University Press, 2001. [本文引用:1]
[47] SUN W C, KUHN M R, RUDNICKI J W. A multiscale DEM-LBM analysis on permeability evolutions inside a dilatant shear band [J]. Acta Geotechnica, 2013, 8(5): 465-480. [本文引用:1]
[48] ZHAO Y X, SUN Y F, LIU S M, et al. Pore structure characterization of coal by NMR cryoporometry[J]. Fuel, 2017, 190: 359-369. [本文引用:1]