土壤渗透^[1]、低渗煤岩增透^[2]、盐岩储气^[3]及花岗岩封存高放核废料^[4]均涉及液体或气体渗透，即流体通过低渗多孔介质的流动规律与相互作用。空间杂乱分布的孔隙-裂隙相互连通构成流体流经的通道，而孔隙网络的空间定量描述与精确重构是准确判断渗透特性的几何基础。对于煤与瓦斯共采工程，低渗煤岩孔隙尺寸甚至达到纳米尺度，直接对探视精细结构的实验设备精度与准度带来了巨大挑战^[5–6]，而低渗理论研究仍是解决问题的窗口。根据Carman-Kozeny经验公式^[7]，多孔介质渗透率与毛细管迂曲度近似呈反比关系，迂曲度的重要性不言而喻，但鲜有基于低渗煤岩开展迂曲度的定量研究。迂曲度反映流体流经毛细孔道历程的复杂程度，研究结果证实孔隙度与迂曲度存在正向关系^[8–9]。迂曲度与渗透率都是孔隙结构的几何相关参量，是有效连通孔隙度的直接反映，因此采用渗透率半经验公式间接测试迂曲度也是主要方法，即采用理论与实验结合的方法可有效避免直接几何测量方法的误差。Purcell^[10]较早提出用汞开展毛管压力与渗透率测试。由于研究对象集中于渗透率，忽略了迂曲度的影响，导致渗透率误差较大^[11–13]。Boundreau^[14]对有限颗粒构成的迂曲度开展模型研究，Sen等^[15]认为孔隙结构分布具有自相似特性，即迂曲度可利用分形理论研究。另外迂曲度带来的流体非线性扰动也是热点^[16–17]，可见迂曲度的定量表述利于开展相关方向研究。早期对迂曲度的认识仅是对渗透率的误差修正，目前仍有文献资料直接对迂曲度给出估值，缺乏测量依据，易引起困惑。随着认识深入，迂曲度被定义为反映孔隙结构流程的弯曲程度，才禀赋了其内在物理意义。国内也有文献提出直接用孔隙度计算迂曲度^[18]，认为迂曲度曲线具有分维特征^[19]，提出了分形毛细管力模型^[20]。

低渗煤岩孔隙通道的弯折程度影响流体的能量损耗，大流阻导致流体通过困难，定量上反映渗透率低。迂曲度刻画弯折程度为定量解析渗透率变化提供了载体，但分形截面的粗糙度在迂曲度方向上的积分效应却未引起足够的重视。基于此，本文侧重在纵向分形迂曲度结构与横向分形截面结构开展低渗理论研究，并建立气-液驱替的分形运动方程。

1 基于Hagen-Poiseuille方程的迂曲度推导

低渗煤岩作为多孔介质，孔隙一般在微米-纳米尺度，微米孔可定义为毛细管孔。孔道既包括盲孔等非有效通道，也包括连通的有效通道。孔隙空间分布十分复杂，截面形状与尺寸不同位置迥异，具有随机性。甚至，在水力等外力作用下还存在着开闭状态。理论上，完备解决基于跨尺度精细结构的定量描述仍存在挑战，分形理论在解决无序结构的自相似问题上表现了极强的适应性^[15]。建立煤岩的分形毛细管模型，如图1所示。

圆柱煤岩纵截面任取矩形截面，考虑一维渗流。沿 $x$ 方向分布一系列不同截面形状的毛细管孔道，假设 $x{\textit{z}}$ 切面上孔道均匀分布，且截面形状相同，流体在孔道中均匀流动。设单元体横截面面积为 $A$ ，长度为 $L$ ，单位面积分布毛细管根数为 $n$ ，每根毛细管弯曲度为 $\tau $ ，毛细管截面为圆形，半径为 $r$ ，单根毛细管孔隙体积为 $V$ 。对应的煤岩基质也假设为毛细管实体，设单位面积分布根数为 $N$ 。

迂曲度反映了孔道在空间分布的复杂程度，存在多种定义^[21]。设 ${L_{\rm e}}$ 为流动路径实际长度，定义迂曲度：

$\tau = {\left( {\frac{{{L_{\rm{e}}}}}{L}} \right)^2}$

(1)

还有贝尔定义 $\tau = {\left( {{{{L_{\rm{e}}}}/L}} \right)^{ - 2}}$ ^[22]；科林斯定义： $\tau = {{{L_{\rm{e}}}}/L}$ ^[23]。

根据Hagen-Poiseuille公式^[24]，截面孔隙通道流量 $q$ ：

$q = nA \cdot \frac{{\Delta p{\text{π}} {r^4}}}{{8\mu \tau L}}$

(2)

式中， $\Delta p$ 为压差， $\mu $ 为黏滞系数。

根据达西定律：

$q = \frac{{KA\Delta p}}{{\mu \tau L}}$

(3)

式中， $K$ 为渗透率。定义孔隙度^[25]：

$\phi = \frac{{An{\text{π}} {r^2}L\tau }}{{AL}} = n{\text{π}} {r^2}\tau $

(4)

联立式（2）、（3）与（4），得：

$\tau = \frac{{\phi {r^2}}}{{8K}}$

(5)

若截面均为圆形孔道，定义水力学半径平均值 ${R_{\rm{h}}}$ ：

${R_{\rm{h}}} = \frac{S}{{\textit{Z}}} = \frac{{{\text{π}} {r^2}}}{{2{\text{π}} r}} = \frac{r}{2}$

(6)

式中， $S$ 为流体流过总截面积， ${\textit{Z}}$ 为总的润湿周长。

根据经验，存在如下关系^[26]：

${R_{\rm{h}}} = \frac{{\phi {D_{\rm{p}}}}}{{6\left( {1 - \phi } \right)}}$

(7)

式中， ${D_{\rm{p}}}$ 为基质等效半径。

结合式（5）、（6）与（7），迂曲度：

$\tau = \frac{{{\phi ^3}D_{\rm{p}}^2}}{{72K\left( {1 - \phi } \right)}}$

(8)

系统由 $N$ 个毛细管实体基质构成，则：

$\frac{{N{\text{π}} D_{\rm{p}}^2}}{4} = \left( {1 - \phi } \right)AL$

(9)

则迂曲度可定义为：

$\tau = \frac{{{\phi ^3}AL}}{{18KN{\text{π}} }}$

(10)

另一种方法是引入毛细管压力公式，不同于Hagen-Poiseuille公式，毛细管压力公式局限于致密岩石、微尺度以下孔道渗流模型。

单根毛管压力^[23]：

$p = \frac{{2\sigma \cos \;\theta }}{r}$

(11)

式中， $\sigma $ 为界面张力， $\theta $ 为接触角。

结合式（2），得：

$q = \frac{{V\Delta p}}{{2\mu }}{\left( {\frac{{\sigma \cos \;\theta }}{{\tau Lp}}} \right)^2}$

(12)

联立达西定律式（3），得：

$K = \frac{{nV}}{{2AL}}{\left( {\frac{{\sigma \cos \;\theta }}{{\tau p}}} \right)^2}$

(13)

孔道一般含有一定的水分或气体，考虑饱和度 $s$ ，则式（13）更改为：

$K = \frac{\phi }{2}{\left( {\frac{{\sigma \cos \;\theta }}{\tau }} \right)^2}\int_0^1 {\frac{{{\rm d}s}}{{{p^2}}}} $

(14)

则迂曲度重新定义为：

$\tau = \sigma \;\cos \;\theta \sqrt {\frac{\phi }{{2K}}\int_0^1 {\frac{{{\rm d}s}}{{{p^2}}}} } $

(15)

饱和度与孔隙体积紧密相关，定义为单根毛细管充满液体的饱和度，即针对单根毛细管取1，针对整体是单根毛细管体积与孔隙度的比值。如果假设毛细管压力 $p$ 恒定，则：

$\tau = \sigma \cos \;\theta \sqrt {\frac{\phi }{{2pK}}} $

(16)

式（15）与（16）即为致密煤岩孔隙迂曲度的毛细管模型算法。结合压汞实验测毛管力，即可测定迂曲度。

2 低渗煤岩分形孔道迂曲度计算

基于科林斯定义，Wheatcraft和Tyler^[26]建立了多孔介质迂曲流管特性的分形标度关系：

${L_{\rm{e}}} = \lambda _{\min }^{1 - {D_{\rm{T}}}}{L^{{D_{\rm{T}}}}}$

(17)

则：

$\tau = \frac{{{L_{\rm{e}}}}}{L} = {\left( {\frac{L}{{{\lambda _{\min }}}}} \right)^{{D_{\rm{T}}} - 1}}$

(18)

式中， ${\lambda _{\min }}$ 为毛细管孔径的最小值， ${D_{\rm{T}}}$ 为迂曲度分形维数，对于线性孔道渗流， $1 \le {D_{\rm{T}}} \le 2$ 。

联立式（15）与（18），得：

${D_{\rm{T}}} = 1 + \frac{{\ln \left( {\sigma \cos \;\theta \sqrt {\displaystyle\frac{\phi }{{2K}}\int_0^1 {\displaystyle\frac{{{\rm d}s}}{{{p^2}}}} } } \right)}}{{\ln \left( {L/r} \right)}}$

(19)

式（19）说明利用压汞实验可以确定迂曲孔道分维。定义分形影响系数：

$\overline \lambda = {\left( {\frac{L}{r}} \right)^{{D_{\rm{T}}} - 1}}$

(20)

如图2（a）为长径比系数和分形影响系数的关系，当迂曲度分形维数 ${D_{\rm{T}}}$ 确定时，分形影响系数 $\overline \lambda $ 与长径比近似成线性关系。图2（b）为迂曲度分维与分形影响系数的关系，当长径比确定时，分形影响系数 $\overline \lambda $ 与迂曲度分形维数 ${D_{\rm{T}}}$ 呈非线性关系， ${D_{\rm{T}}}$ 与 $L/r$ 越大影响越明显，说明 ${D_{\rm{T}}}$ 需要长程累积效应，短程情况下孔道横断面尺寸与流程在一个数量级上，圆形孔道阻滞特性并不明显；另一方面，应考虑迂曲度在流程上的非线性积分。

煤岩取自平顶山矿区，埋深约1 000 m，加工成几何尺寸： ${\varPhi }$ 25 mm $ \times $ $H$ 30 mm，开展高压压汞实验，共3个样品。实验前利用工业CT对孔隙扫描，精度为30 μm，难见孔隙分布，可见煤岩孔隙尺寸十分微小。设计压力从0.007到200 MPa，对应孔隙尺寸为2.5 μm到4 nm（样品1与样品2）、25 μm到4 nm（样品3）。连续性加压，进汞速度控制为7 MPa/min，设备精度为0.01%。

图3为3个煤样的压汞曲线图，进汞及回汞曲线差异明显，初始压汞部分存在陡峭部分，即压力极速增长，说明汞未有效压入，对应的孔隙部分存在缺失，即煤岩十分致密。实验测定样品1： ${\phi _1} = 5.8\% $ ； ${K_1} = 0.003 \times 10^{-15}$ m/s；样品2： ${\phi _2} = 5.6\% $ ； ${K_2} = 0.007\times $ 10^–15 m/s；样品3： ${\phi _3} = 4.1\% $ ； ${K_3} = 0.001\times 10^{-15}$ m/s；取 $\sigma = $ 0.48 N/m； $\theta = $ 140°（对于汞，均为定值），利用式（15）结合图4数据积分，得迂曲度： ${\tau _1} = 2.25,$ ${\tau _2} = 1.27,{\tau _3} = 2.78$ ，平均迂曲度 $\overline \tau = 2.1$ 。结合式（19）与孔隙半径分布，可获得3个样品的迂曲分维 ${D_{{\rm{T}}1}}{\rm{ = }}1.071$ ， ${D_{{\rm{T}}2}}{\rm{ = }}1.021$ 和 ${D_{{\rm{T}}3}}{\rm{ = }}1.103$ 。定义迂曲系数 $\tau {D_{\rm{T}}}$ 。

图4为渗透率与迂曲系数的关系，呈直线关系。迂曲系数可以很好的反映渗透率变化，是一种几何结构参数，包括了迂曲度与迂曲分维的影响。

3 低渗煤岩分形截面分维理论推导

由于毛细管微-纳米的结构“半径”，毛细管横截面润湿周长分形特性常被忽略，甚至同迂曲度分形概念混淆。假设毛细管截面为正 $n$ 边形，外接圆半径为 $r$ 。定义 $n$ 边形的边长为 $a$ ，中心到边的垂直距离为 $h$ ，则：

$a = 2r\sin \frac{{\text{π}} }{n}$

(21)

$h = r\cos \frac{{\text{π}} }{n}$

(22)

根据Carman-Kozeny方程^[7]，圆形毛细管截面渗透率 ${K_{\rm{c}}}$ ：

${K_{\rm{c}}} = \frac{{{\phi ^3}}}{{2{\tau ^2}{S^2}}}$

(23)

则正 $n$ 边形与圆形截面渗透率关系为：

${K_n} = \frac{{\sin \displaystyle\frac{{\text{π}} }{n}}}{{\displaystyle\frac{{\text{π}} }{n}}} \cdot {\cos ^3}\displaystyle\frac{{\text{π}} }{n} \cdot {K_{\rm c}}$

(24)

当 $n$ 趋近于正无穷时，存在：

$\mathop {\lim }\limits_{n \to + \infty } {K_n} = {K_{\rm{c}}}$

(25)

即正 $n$ 边形毛细管渗透率等同于外接圆半径为 $r$ 圆截面毛细管模型的渗透率。

图5为基于正六边形的5次分形截面迭代，考虑正 $n$ 边形毛细管截面的自相似分形，如以正六边形为母元。第1次迭代：1）将母元每条边等分三份；2）去掉中间三分之一边长；3）以三分之一边长为子元边长，重生一个正六边形。依次循环第2、3、4次迭代，其他参数同上。

相应地4次迭代后，可以建立正 $n$ 边形分形毛细管渗透率与圆形毛细管截面渗透率 ${K_{\rm{c}}}$ 的关系（表1），比值系数仅与迭代次数和母元边数有关系。

根据Washburn方程知，气-液两相流运动界面方程同时考虑截面积 $S$ 与界面周长 $C$ 的影响。

对于分形结构，量测的面积存在近似关系（表2）：

${S_i} = {S^0}{\varepsilon _i}^{1 - {D_{\rm s}}}$

(26)

式中： ${S_i}$ 为第 $i$ 次迭代对应的截面积，对于严格自相似分形， ${S_0} = {S^0}$ ； ${\varepsilon _i}$ 为对应面积的码尺； ${D_{\rm s}}$ 为截面积的分维。

量测的周长存在近似关系（表3）：

${C_i} = {C^0}{\delta _i}^{1 - {D_{\rm C}}}$

(27)

式中： ${C_i}$ 为第 $i$ 次迭代对应的周长，对于严格自相似分形， ${C_0} = {C^0}$ ； ${\delta _i}$ 为对应周长的码尺； ${D_{\rm C}}$ 为周长的分维。

根据体积覆盖法^[27]，可以建立对应 ${\varepsilon _i}$ 、 ${\delta _i}$ 码尺与覆盖格子数 ${N_{\varepsilon i}}$ 与 ${N_{\delta i}}$ （ $i = 0,1,2,3 \cdots $ ）的关系，定义最新生成分形子结构，即最小正 $n$ 变形的面积 ${\varepsilon _i}$ ：

${\varepsilon _i} = \frac{{3n}}{2}a_i^2\tan \frac{{2{\text{π}} }}{n},i = 0,1,2,3, \cdots $

(28)

定义最新生成分形子结构的最小边长 ${\delta _i}$ ：

${\delta _i} = {a_i},i = 0,1,2,3, \cdots $

(29)

则截面积的分维 ${D_{\rm s}}$ ：

${D_{\rm s}} = - \frac{{\lg \left( {{N_{\varepsilon ,i}}/{N_{\varepsilon ,0}}} \right)}}{{\lg \left( {{\varepsilon _i}/{\varepsilon _0}} \right)}}$

(30)

则周长的分维 ${D_{\rm C}}$ ：

${D_{\rm C}} = - \frac{{\lg \left( {{N_{\delta ,i}}/{N_{\delta ,0}}} \right)}}{{\lg \left( {{\delta _i}/{\delta _0}} \right)}}$

(31)

4 气液两相流分形运动方程

图6为气液两相驱替模型。

如图6所示，基于毛细管结构建立气液两相驱替模型，取迂曲方向上一小直段 $\Delta L$ ，包含气液两相及界面，气体为 ${{\rm{N}}_2}$ ，液体为 ${{\rm{H}}_{\rm{2}}}{\rm{O}}$ 。某一个时刻，液体长度为 $x$ ；下一时刻，液体长度为 $x + \Delta x$ ， $\Delta x$ 为负的增量。忽略气体质量，对液体建立Washburn方程^[28–29]，主要包括三方面影响：驱替压力 ${F_{\rm p}}$ 、毛细管力 ${F_{\rm c}}$ 与壁面阻力 ${F_{\rm f}}$ ，考虑分形截面影响，则分别为：

${F_{\rm p}} = \Delta p \cdot {S^0}{\varepsilon _i}^{1 - {D_s}}$

(32)

式中， $\Delta p$ 为驱动力， $\Delta p = {p_{\rm in}} - {p_{\rm out}}$ 。

${F_{\rm c}} = {p_{\rm c}} \cdot {S^0}{\varepsilon _i}^{1 - {D_{\rm s}}}$

(33)

式中， ${p_{\rm c}}$ 为毛细管阻力。

${F_ {\rm f} }= \tau \left( {x + {\rm{d}}x} \right) \cdot {C^0}{\delta _i}^{1 - {D_{\rm C}}}$

(34)

式中， $\tau $ 为粗糙壁面的平均黏性切应力， $\tau = \left( {{{4\eta }/r}} \right){v_{\rm m}}$ ， $\eta $ 为水的黏度， ${v_{\rm m}}$ 为水的平均流速， ${v_{\rm m}} = {{{\rm{d}}x}/{{\rm d}t}}$ 。

根据牛顿第二定律：

$ma = {F_{\rm p}} + {F_{\rm c}} + {F_{\rm f}}$

(35)

将式（32）、（33）与（34）代入式（35），即：

$\begin{aligned}[b] \rho \cdot {S^0}{\varepsilon _i}^{1 - {D_{\rm s}}} \cdot x{\rm{d}}x = & \left( {\Delta p - {p_{\rm c}}} \right){S^0}{\varepsilon _i}^{1 - {D_{\rm s}}} - \\& {C^0}{\delta _i}^{1 - {D_{\rm C}}} \cdot \frac{{4\eta }}{r}\frac{{{\rm{d}}x}}{{{\rm{d}}t}}\left( {x + {\rm{d}}x} \right)\end{aligned} $

(36)

忽略惯性项与高阶黏性切应力，则

${C^0}{\delta _i}^{1 - {D_{\rm C}}} \cdot \frac{{4\eta }}{r}\frac{{{\rm{d}}x}}{{{\rm{d}}t}}x = \left( {\Delta p - {p_{\rm c}}} \right){S^0}{\varepsilon _i}^{1 - {D_{\rm s}}}$

(37)

此方程即为Washburn方程分形形式，简称分形Washburn方程。对方程积分，可以得到界面位置方程：

${x^2} = \frac{{{S^0}{\varepsilon _i}^{1 - {D_{\rm s}}}}}{{{C^0}{\delta _i}^{1 - {D_{\rm C}}}}} \cdot \left( {\Delta p - {p_{\rm c}}} \right)\frac{r}{{2\eta }}t$

(38)

式中， ${{{S^0}{\varepsilon _i}^{1 - {D_{\rm s}}}}/{{C^0}{\delta _i}^{1 - {D_{\rm C}}}}}$ 定义为分形截面粗糙度系数。

同时可以得到界面运动 $x$ 处平均速度 ${v_{\rm m}}$ 的方程：

${v_{\rm m}} = \frac{{{S^0}{\varepsilon _i}^{1 - {D_{\rm s}}}}}{{{C^0}{\delta _i}^{1 - {D_{\rm C}}}}}\left( {\Delta p - {p_{\rm c}}} \right)\frac{r}{{4\eta x}}$

(39)

5 分形运动方程应用分析

利用Nano CT对煤岩进行纳米-微米尺度分析，样品近似圆柱形，尺寸约65 μm。重构后截取60 μm立方体，孔隙结构十分复杂（图7）。图8为基于纳米与微米尺度的孔径与孔隙体积关系，并对内部孔隙结构进行定量分析，孔径在0.1～20 μm均有分布，主要集中2～20 μm，可见在微米尺度上孔周表面仍有可能存在更细微的结构。

进一步基于低场NMR开展两相N₂–H₂O驱替实验研究，驱替压路径设计如图9所示。低渗煤岩样品尺寸为 ${\varPhi }$ 25 mm $ \times $ $H$ 50 mm。渗透压 $\Delta p$ 设计路径为1、3、5、6与7 MPa，驱替时长分别为1、2、3、6与4 h。

根据孔径分布，选取18 μm为边长的正六边形为母元进行分形结构分析，则对应的码尺分别为18、6、2、0.67与0.22 μm，基本可以覆盖主要的孔径分布，根据式（30）与（31），可计算得到面积分维 ${D_{\rm S}} =$ $ 1.092$ （图10（a））与周长分维 ${D_{\rm C}} =1.601$ （图10（b））。

针对母元，毛细管力 ${p_{\rm c}}$ 相比驱替压 $\Delta p$ 可以忽略不计，对于常温水： $\eta = 0.001\;{\rm{Pa}} \cdot {\rm{s}}$ 。根据式（46）可以建立驱替时间与界面运动距离的关系（图11（a））；联合式（46）与（47）可以建立界面速度与驱替时间的关系（图11（b））。

界面运动距离随着时间的增长，呈现非线性增长，随着驱替压的提高，增长幅度趋于缓慢；对于速度，随着时间增长快速衰减（图11（a））。

在1、3、5、6与7 MPa下，启动速度分别为111.4、192.9、249.1、272.9与294.7 μm/s，可见启动速度随着驱替压的增加也在增加。但速度随距离的衰减极快（图11（b）），不同驱替压下界面速度与界面距离单调递减，遵循相同的衰减路径。初始1 h后的衰减速度，1 MPa速度最大，7 MPa速度最小；对应运动距离，1 MPa距离最小，7 MPa距离最大。压差提供了长程驱替的能量，压差大，驱替距离远，速度衰减也最快。初始速度衰减快的原因主要为 $\left( {{S^0}{\varepsilon _i}^{1 - {D_{\rm s}}}} \right)/\left( {{C^0}{\delta _i}^{1 - {D_{\rm C}}}} \right)$ ，即分形结构的初始阻滞作用。

图12为5 MPa作用下，不同迭代次数分形结构界面运动距离、速度与驱替时间的关系。界面运动距离，母元结构中移动最远：62.13 cm；1次迭代后距离衰减至26.05 cm、衰减58.1%；2次迭代后距离为10.92 cm、衰减82.4%；3次迭代后距离为4.58 cm、衰减92.6%；4次迭代为3.33 cm、衰减94.6%。可见，分形结构不同于光滑圆形结构，随着壁面粗糙程度增加，阻滞作用越趋于明显。对于液体，壁面粗糙度是不可忽略的因素。速度降落也很明显，初始速度，母元结构最大；第四次迭代分形结构速度最小。因此，壁面粗糙度同时对驱替液体界面运动距离和速度存在阻滞作用，这就需要更高的驱替压抵抗新增的阻滞。

基于压汞试验并计算平均迂曲度 $\overline \tau = 2.1$ ，样品高为5 cm，则流动路径实际长度 ${L_{\rm{e}}}{\rm{ = }}10.5\;{\rm{cm}}$ ，根据图13，实验驱替下累计界面运动距离可知，对应的启动驱替压约为3 MPa。结合NMR实验过程驱水信号分布（图14），从初始饱和至1 MPa驱替压，驱替时长1 h，效果不是很明显。而3 MPa下驱替2 h后，具有明显的效果。可见分形运动方程可以有效地反映实验情况。随着驱替压的增加，更复杂分形结构孔或者更小直径孔内的水分陆续被驱替。当驱替压为7 MPa时，累计驱替距离47.71 cm，约是 ${L_{\rm{e}}}$ 的5倍，相当于从低压到高压增长过程中，煤岩被成功驱替了5次，成功驱替意味着水分从煤岩一端到另一端。由此可见，为实现成功驱替，需要保证界面运动距离持续增长，驱替压的持续提高可实现这一过程。但从驱替效果看，经历3、5与6 MPa后，大部分水分已经被驱排，7 MPa下仍可有效驱排复杂分形结构与微-纳米孔内水分。实际操作中，需要根据具体需要确定适当的驱替压，因为在高渗透压下长时驱替仅是极小的含水量。综上，基于分形运动方程利于科学定量确定驱替压范围，对应的基于分形截面模型计算驱替距离与速度，有效地符合实验结果。

6 结论与讨论

1）建立了分形毛细管模型，推导了基于Hagen–Poiseuille方程的低渗煤岩迂曲度两种数学表达式。第一种为理论公式，属传统定义；第二种结合压汞试验，简化了迂曲度的计算。进一步引入孔径-长度分形标度关系，获得了迂曲孔道的分维表达式。

2）开展了低渗煤岩压汞实验，获得了煤的迂曲度与迂曲分维。定义迂曲系数为迂曲度与迂曲分维的乘积，发现渗透率与迂曲系数呈线性关系。

3）建立了低渗煤岩分形截面模型，基于Carman–Kozeny方程计算渗透率验证了分形截面的可靠性。建立了截面积与周长的分维表达式，同时基于Washburn方程建立了界面运动位置方程与运动速度方程。

4）基于Nano–CT重构了低渗煤岩纳米-微米孔隙结构，获得了孔径与孔隙体积分布，结合理论模型，计算了面积分维与周长分维。并结合分形运动方程建立了运动距离、速度同时间的关系；速度同距离的关系；不同分形结构的驱替规律。结果表明：系数 $\left( {{S^0}{\varepsilon _i}^{1 - {D_{\rm s}}}} \right)/\left( {{C^0}{\delta _i}^{1 - {D_{\rm C}}}} \right)$ 对运动距离与速度均存在明显的阻滞作用， $\left( {{S^0}{\varepsilon _i}^{1 - {D_{\rm s}}}} \right)/\left( {{C^0}{\delta _i}^{1 - {D_{\rm C}}}} \right)$ 主要反映了面积分维与周长分维的共同作用。

5）基于NMR实验开展了两相N₂-H₂O驱替实验，结合理论模型获得了实验驱替下界面运动距离与速度分布。通过驱水信号分布证实了分形运动方程的有效性。

综上，在孔隙度相同条件下，煤岩低渗机制应考虑毛细孔道纵向与横向的影响。纵向主要为迂曲系数；横向主要为分形截面粗糙度系数的影响。而迂曲系数与分形截面粗糙度系数在分析压汞实验、Nano-CT重构实验与NMR实验结果时也是分开考虑的，同时考虑两个权重的影响仍有待研究。需要说明的是：

1）式（10）建立在Hagen-Poiseuille公式与达西定律的基础上，是理论求解，但模型要求存在并行分形毛细管。这点在实际验证时又带来了新的不确定性，即截面毛细管分布密度的确定问题，优点是式（10）并不要求模型严格建立在微米尺度。式（16）是式（15）的简化，是在式（10）建立原理基础上，结合压汞试验，引入饱和度巧妙避免了毛细管分布密度问题。但这种方法测定的是所有分形毛细管的综合效应，即迂曲度概念是建立在分形毛细管串联的基础上。

2）毛细管孔道，上述迂曲度定义既无法理清单条长程孔道迂曲度差异；更是忽略了迂曲度禀赋界定渗透性含义。对于分叉孔道，仍有必要基于能量原理^[30]深入研究分叉孔道迂曲度的基本关系。