基于分形理论的非达西渗流参数研究

冯雪冬，张毅，由晗杨，李文英，马艳飞

(1. 山东理工大学 资源与环境工程学院, 山东 淄博 255049；2.山东理工大学 农业工程与食品科学学院, 山东 淄博 255049)

石油在开采、炼制、贮运及使用过程中往往伴随着严重的环境问题，其中地下水污染是环境问题之一[1-2]。研究地下水污染，需要从石油污染物渗流方程的角度出发，达西定律作为多孔介质中饱和渗流的基本方程，描述了流体渗流速度与压力梯度之间的线性关系，目前被广泛应用于地下水、石油工程等领域[3]。然而，理论研究与室内试验发现，当流体黏度或渗流速度增大时，渗流速度与压力梯度之间呈非线性关系，达西定律不再成立。因此，对非达西渗流的研究成为当前渗流领域的研究热点。

从现有研究成果来看，描述非达西渗流的方程有很多，Forchheimer方程因其具有明确的物理意义，被广泛应用于非达西渗流的计算中[4]。但目前国内对Forchheimer方程中渗流参数的研究主要集中在数值模拟方面[5-6]，多孔介质的粒径构成及孔隙结构等特征与非达西渗流参数的关系研究却相对较少。同时，由于多孔介质的颗粒和孔隙具有不规则性、自相似性和非线性等特点[7]，使得多孔介质孔隙结构与非达西渗流参数之间的关系很难用传统线性分析方法进行定量化描述。

分形理论作为定量描述多孔介质的新方法，其分形维数值既可以准确反映多孔介质的粒径构成，又能真实描述多孔介质的孔隙结构[8]。因此，本文利用分形理论，探讨了多孔介质分形特征与Forchheimer型非达西渗流参数之间的关系，以期为石油工程领域的非达西渗流研究提供一定的理论依据。

1 材料与方法

1.1 试验材料

试验所用多孔介质来源于山东省东营市广饶县和淄博市临淄区的农田，按0～20 cm，20～40 cm和40～60 cm深度由上向下分层采集各土层的土壤样品，每层重复取样3次，其中广饶县0～20 cm，20～40 cm和40～60 cm的土壤编号依次为GR1、GR2和GR3，临淄区0～20 cm，20～40 cm和40～60 cm的土壤编号依次为LZ1、LZ2和LZ3；所采多孔介质样品带回实验室自然风干、除杂，过20目筛后，储存在密封容器中；采用Mastersizer 2000激光粒度分析仪测定多孔介质样品的颗粒组成，并按照国际制分级标准来划分土壤质地类型。

试验所用流体为0#柴油，购自中国石化齐鲁石化股份有限公司，密度为0.82 g·cm-3，动力黏度为3.11 mPa·s，油水界面张力为37.30 mN·m-1，表面张力为23.90 mN·m-1。

1.2 渗流试验

试验采用自制的一维土柱渗流试验装置，如图1所示。为保证试样的均匀性，试验时按标准化分层方法进行装样，每次2 cm厚，装样高度为20 cm；进油端采用蠕动泵进样，并且在供油管上设有分流球阀，以控制进油端油压；土柱出口端连接恒定液面高度的油箱，以保持土柱出流端液面高度恒定；土柱进出油端分别装有测压管，用于测定进出油端的油压；柴油的渗流速度由量筒和秒表进行测定。

图1 试验装置示意图Fig.1 The schematic diagram of the experiment equipment

1.3 渗流理论方程

1856年法国工程师达西通过大量试验研究，发现了渗流速度与压力梯度之间的线性关系[9]，达西渗流方程可表示为：

(1)

式中：J为压力梯度，MPa·m-1；ρ为流体密度，kg·L-1；g为重力加速度，m·s-2；K为渗透系数，m·s-1；v为渗流速度，m·s-1。

Forchheimer型非达西渗流方程与达西渗流方程相比，其主要区别在于增加了一个二次项，意义上则表示压力梯度的损失不再完全由黏滞力所决定，而是由黏滞力和惯性力共同决定[6]，其具体表达形式为：

(2)

式中：β为非达西渗流因子，m-1；Kf为非达西渗透系数，m·s-1。

1.4 分形维数计算方法

本试验用粒径的重量分布对多孔介质分形维数进行计算[10-11]：

(3)

式中：M(r

对该式两边取对数，以lgR，lg (M(r

D=3-b

(4)

2 结果与分析

2.1 压力梯度与渗流速度的关系

调节分流球阀，测定不同压力梯度下柴油的渗流速度，将试验结果绘制成图，如图2所示。由图2可知，多孔介质中柴油的渗流曲线均呈类抛物线，表明渗流速度与压力梯度不再满足线性关系，很可能发生了非达西渗流。为进一步证实上述分析结果的正确性，采用Forchheimer方程对试验结果进行拟合，拟合结果见表1。由表1可知，Forchheimer方程的拟合效果较好，相关系数均大于0.96，说明柴油在渗流过程中确实发生了非达西渗流。

表1 渗流速度与压力梯度的曲线拟合方程及相关系数Tab.1 The fitting equations of seepage velocity and pressure gradient

2.2 多孔介质的分形特征

图3为lgR与lg (M(r

与广饶地区相比，临淄地区多孔介质的分形维数值相对较大，但两地区分形维数值的变化趋势也存在相似之处：随着土壤深度增加，分形维数值也相应增大，在40～60 cm土层都达到最大值。从表2中还可以发现，多孔介质的分形维数值与粉粒、粘粒含量的变化趋势基本一致，即多孔介质中粉粒、粘粒含量越高，土壤颗粒分形维数值就越大。

2.3 分形维数与粒径含量的关系

为进一步探讨分形维数值与粒径含量的关系，对粒径含量与分形维数值进行回归分析，拟合结果如图4所示。由回归分析结果可以看出，分形维数值与砂粒含量之间的相关性达到极显著水平(P<0.01)，相关关系为负相关(R=-0.985)；分形维数值与粉粒含量之间呈极显著正相关关系(P<0.01，R=0.970)，与粘粒含量也呈极显著正相关关系(P<0.01，R=0.995)。分形维数值与多孔介质粒径含量的相关性较高，这表明分形维数值能够在一定程度上反映多孔介质的粒径组成，即分形维数值的变化与粉粒、粘粒含量密切相关，分形维数值越大，粉粒、粘粒含量越高，质地越细，微小孔隙越多，结构越复杂，这与张世熔等[14]的研究结果一致。

表2 不同多孔介质的粒径组成及分形维数Tab.2 Particle size composition and fractal dimension of different porous media

(a)分形维数与砂粒含量的关系 (b)分形维数与粉粒含量的关系 (c)分形维数与黏粒含量的关系图4 分形维数与粒径含量的关系Fig.4 Relationship between fractal dimension and particle size content of soils

2.4 分形维数与非达西渗流参数的关系

从Forchheimer方程表达式可以看出，非达西渗透系数和非达西渗流因子分别代表着渗流过程中粘滞力对柴油的影响和惯性力对柴油的影响。当多孔介质孔隙较小、流速较低时，黏滞力占主导地位，惯性力可忽略，这时压力梯度与渗流速度满足线性关系，服从达西定律；当多孔介质孔隙增大、流速加快时，惯性力的影响逐渐增大，压力梯度与渗流速度的线性关系发生偏移，即发生非达西渗流。因此，Forchheimer方程渗流参数的变化与多孔介质性质(尤其是孔隙结构)密切相关。

本研究中分形维数值能反映多孔介质的粒径构成，而粒径构成又影响着多孔介质的孔隙结构，因此通过建立分形维数值与非达西渗流参数之间的回归方程，探讨它们之间的相关关系，可以为进一步研究Forchheimer型非达西渗流参数特征提供一定的研究思路和方法。以非达西渗流参数，分形维数值为横、纵坐标，进行曲线拟合，回归方程及相关系数如图5、图6所示。

由图5可以看出，分形维数值与非达西渗透系数之间呈线性负相关关系，且二者的相关程度较高。从微观角度来看，分形维数值越大，多孔介质中粉粒、粘粒含量越高，多孔介质越紧实，微小孔隙越多，从而使得多孔介质的渗透性变差。与非达西渗透系数的变化趋势相反，非达西渗流因子随分形维数值的增大而增大(图6)，反之，当分形维数值变小，多孔介质颗粒变粗，压力梯度与渗流速度之间的非线性关系越不明显。

综上所述，多孔介质的粒径大小、孔隙结构等在很大程度上决定着非达西渗流参数的大小。分形维数值较直观地反映了多孔介质的颗粒组成，通过建立分形维数值与非达西渗流参数之间的相关关系，从而可以定量地描述Forchheimer型非达西渗流的参数特征。

3 结论

1) 柴油在多孔介质中的渗流为非达西渗流，Forchheimer方程能较为准确地描述压力梯度与渗流速度之间的关系。

2) 多孔介质具有统计意义上的自相似性，即分形性，其中分形维数值与砂粒含量呈极显著负相关，

图5 分形维数与非达西渗透系数的关系Fig.5 Relationship between fractal dimension and non-Darcy permeability coefficient

图6 分形维数与非达西渗流因子的关系Fig.6 Relationship between fractal dimension and turbulence factor

与粉粒、粘粒含量呈极显著正相关。

3) 分形维数值与非达西渗透系数之间呈线性负相关关系，与非达西渗流因子之间呈线性正相关关系，因此，分形维数可以作为描述非达西渗流参数变化趋势的重要指标。