刘晓明, 杨本亮, 申琦瑜

(湖南大学 土木工程学院，长沙 410082)

1 研究背景

岩石裂隙的渗透规律因水资源利用、地下空间开发、能源开采而受到广泛重视。自然界中岩石类型繁多，在地壳浅层、表层广泛分布着具有气孔构造的岩石，同其他种类岩石一样，气孔构造岩石裂隙的渗透性也主要受分布在其中的裂隙控制，探明其渗透特性是在气孔构造岩石地区开展市政、水利工程、油气开采的一项重要前期工作。

对于岩石裂隙渗透性的研究，Poiseuiile对均质恒定运动条件下的黏性不可压缩流体进行理论推导[1]，得到裂隙渗流流量与裂隙宽度的立方定理。Lomize[2]、Louis[3]通过大量的单裂隙渗流试验，验证了层流时光滑裂隙的立方定理。岩石裂隙在自然界广泛存在，但并不光滑。Barton等[4]和周创兵等[5]相继进行了粗糙裂隙渗透试验，发现粗糙裂隙面凸起度Δ与裂隙宽度的比值对裂隙过流能力影响显著，分别提出了相应的修正立方公式。胡少华等[6]建立了考虑裂隙几何形貌的临界雷诺数模型。鞠杨等[7]通过试验发现裂隙的渗透能力随裂隙粗糙度增加呈指数下降。

随着计算流体力学(Computational Fluid Dynamics，CFD)理论逐渐成熟，数值模拟成为研究裂隙渗流一种新手段。程勤波等[8]基于Navier-Stokes(N-S)方程建立了岩石裂隙入渗水流数值模型，结果表明N-S方程能精确地描述岩石裂隙的水流运动过程。徐维生等[9]统计现场工程中裂隙宽度的分布特征，根据相应特征基于数值生成仿粗糙裂隙，探明了多种渗流本构关系下的单宽流量公式。于可等[10]通过Fluent进行数值模拟，发现粗糙裂隙在高速渗流时，渗透率与粗糙度呈负相关，与裂隙宽度呈正相关；段慕白等[11]和Javadi等[12]运用CFD对裂隙渗透进行了数值模拟，得到了相似结论。赵强等[13]对单裂隙微凸体粗糙裂隙进行仿真模拟，结果表明微凸体粗糙裂隙模型中，同一裂隙渗流通道中的隙宽较大处流速缓慢，隙宽较小处流速较大。

综上所述，目前岩石裂隙渗流研究主要针对一般粗糙裂隙，而对裂隙面具有下凹气孔的岩石裂隙渗流研究较少。鉴于此，本文以气孔构造玄武岩为例，开展气孔构造岩石裂隙渗透试验，经试验数据拟合，得到了裂隙面构造影响因子λ与雷诺数Re的关系。并基于数值软件建立模拟气孔构造裂隙渗流模型，在印证室内试验所得规律合理的基础上，进一步探讨了裂隙面气孔的存在对裂隙渗流的影响机理。

2 气孔构造玄武岩裂隙渗透试验

2.1 气孔构造玄武岩气孔特征

本研究所用的气孔构造岩石为海南省海口市的玄武岩，属第四系上更新统道堂组火山岩，呈灰色、灰褐色的隐晶质结构，具有气孔、微孔状构造。经过实地调查后，选取如图1所示的4种气孔特征差异较大的玄武岩试块进行研究。

图1 气孔构造玄武岩试样Fig.1 Basalt samples of vesicular structure

摄影获取图1所示玄武岩试块图像后，采用常用绘图软件AutoCAD中光栅图像功能对玄武岩气孔图像进行描绘，将描绘后的图片导入到Photoshop中，统计气孔图像像素后进行气孔面积的计算，计算得到切割后玄武岩表面的气孔特征参数如表1所示。

表1 气孔构造玄武岩气孔特征参数Table 1 Vesicular characteristic parameters of basaltwith vesicular structure

2.2 试验设备和方法

研究所用渗透试验设备如图2所示，设备分进水箱、出水箱、渗透试样3部分，进水箱高1.3 m，箱体连接进水管、测压管；出水箱高0.6 m，试验水头的大小通过控制出水管的高度设定。测压管读数精度为1 mm，即本装置可测试的最小水头差为1 mm，对应水头梯度为0.005。

图2 试验设备Fig.2 Test equipment

带裂隙的岩石试样采用两块切割整齐的玄武岩夹持一定厚度的铜条来模拟，如图3(a)所示。试样尺寸为200 mm×120 mm×60 mm(长×宽×高)。这样裂隙面的长度为120 mm，因所选铜条宽10 mm，所以裂隙渗流截面长度W为100 mm。然后将试样用1 mm厚的乳胶膜和四块铝板包围，使用夹具夹持后，可将试块裂隙周边密封，夹持后的试样如图3(b)所示。

图3 试样组装Fig.3 Assembling of sample

夹持中需采用扭力扳手控制力度，确保每次施加在夹具上的力相同，从而控制裂隙宽度相同。将夹具加持后的试样组装至图2所示试验设备后进行渗透试验。在裂隙渗透试验中，等效水力缝宽的获取尤为重要，王缓等[14]对裂隙等效水力隙宽确定问题进行了探讨，提出了将试验所得裂隙渗流量代入立方定律反求得到等效水力缝宽。因此，本文通过光滑裂隙渗流的渗流量进行反算获得裂隙渗流等效水力缝宽，采用光滑有机玻璃板模拟裂隙，选用的3种铜条，实测厚度分别为0.41、1.01、1.92 mm，进行光滑裂隙渗透试验，获得试验流量Qg,在层流状态下，试验流量与立方公式计算流量Q相等，即Qg=Q。计算获得对应的裂隙等效水力宽度b分别为0.406、0.933、1.643 mm，裂隙等效水力宽度小于实际观测值，这与Yeoiw等[15]得出的水力等效隙宽比机械隙宽小的结论一致。

开展玄武岩裂隙渗流试验时，先将夹具夹持后的玄武岩试样组装到试验设备上，打开设备的进水阀门，通水后对试样进行全方位转动，以排出试样气孔中的空气，待无气泡排出后，开始渗透试验，记录进、出水箱两侧水头差和对应的稳定出水流量，计算得到水头梯度J和单位时间过流流量Qs。通过调节进出水管高度，实现变换水头，重复以上步骤，直至所有试验完成。本研究共测试了4组玄武岩岩石裂隙试样，每种试样对应3个裂隙宽度、平均9个水头差，共得到108个水头梯度J-单位时间过流流量Qs数据对。

2.3 试验结果及分析

立方定理公式为

(1)

式中：q为单位长度的裂隙计算流量(cm2/s)；ρ为液体密度(kg/m3)；g为重力加速度(m/s2)；μ为黏滞系数，与温度有关；J为水头梯度；b为等效裂隙宽度。

光滑裂隙理论渗流流量为

Q=qW。

(2)

式中试样裂隙渗流截面长度W为10 cm。

将不同等效裂隙宽度对应的试验结果分别绘制在图4中，图中实线由式(2)计算得出。

图4 不同等效裂隙宽度下的玄武岩裂隙Qs-J曲线Fig.4 Curves of hydraulic gradient and flow in basaltfracture of different equivalent fracture widths

从图4的试验结果可以发现，与光滑裂隙相比，当等效裂隙宽度不同，气孔构造的存在对裂隙渗流的影响不同：当等效裂隙宽度较小(b=0.406 mm)时，4个玄武岩试样的过流流量均大于立方公式计算流量，说明此时气孔构造对裂隙水流通过是有利的；当等效裂隙宽度较大(b=1.643 mm)时，4种玄武岩的过流流量均显著小于立方公式计算流量，这说明在等效裂隙宽度较大时，气孔对裂隙水流通过是不利的。且在等效裂隙宽度较小(b=0.406 mm)、较大(b=1.643 mm)条件下，气孔面积占比最大的3#试样裂隙渗流量较气孔面积占比最小的1#试样裂隙(较接近光滑裂隙)渗流量更偏离立方定理公式计算流量，可见，裂隙面的气孔对裂隙渗流的影响确实存在，而且与渗流条件有关。

为分析试验渗流流量变化规律，将试验结果转换形式，采用流体分析常用参数雷诺数Re描述水流特征，以揭示这种特殊构造裂隙的渗透规律。在单裂隙渗透试验中，流体雷诺数Re按下式计算，即

(3)

式中：A为裂缝开口的横截面积;Dh为水力直径(cm)，等于4A/p，p表示裂隙的周长(cm)，p=2(W+b)。

雷诺数Re揭示了惯性运动与黏性运动之间的比率，根据有效裂隙宽度计算获得。通常，较小的雷诺数表示层流，其中黏性力起主要作用，而较大的雷诺数表示湍流惯性力主导的流动。

引入裂隙面构造影响因子λ，其计算方法为

(4)

λ实际上是气孔构造裂隙渗流流量与相同条件下光滑裂隙渗流流量的比值，它是一个无量纲参数。根据单裂隙渗流试验数据，可分别依据式(3)和式(4)获得雷诺数Re和影响因子λ。将试验数据计算所得的λ与Re的关系绘制在图5中，同时为了便于观察规律，图中还绘制了λ=1和Re=200参考线。

图5 雷诺数Re与裂隙面构造影响因子λ的变化关系Fig.5 Variation of Reynolds number (Re) againstinfluence factor (λ) of fracture surface structure

由图5可见，当雷诺数相对较小(Re<200)时，气孔构造裂隙的渗流流量大于光滑裂缝流量；当雷诺数Re>200后，气孔构造裂隙渗流流量比光滑裂缝的裂缝流量小，且随着雷诺数Re的增大，裂隙面构造影响因子λ减小趋势变缓。根据图5中数据，可得回归公式：

λ=13Re-0.5。

(5)

式(5)表明在不同的雷诺数下，气孔构造对裂隙水流流动的影响程度并不同。随着雷诺数Re的增大，裂隙面气孔的存在，对裂隙水流的影响逐渐由利于水流通过转化为不利于水流通过。根据式(1)—式(5)，整理可得气孔构造裂隙渗流流量公式为

(6)

式中p为裂隙的周长。

由于气孔构造裂隙的宽度远小于长度，故裂隙的周长p近似等于裂隙长度W的两倍，对气孔构造裂隙渗流流量公式简化为

(7)

与立方公式表征光滑裂隙的渗透规律一样，式(7)也可表征气孔构造裂隙的渗透规律。

受试验条件所限，试验中最小水头梯度为0.115，最大水头梯度为1.45，完成水流雷诺数Re的范围6.54～773.11，故试验所得渗透规律适用范围受限；另外，气孔对裂隙渗流影响机理尚不清晰，因此对气孔构造玄武岩裂隙渗流进行 CFD 数值计算，以期完善气孔构造岩石裂隙渗流规律的研究。

3 模拟气孔构造渗流模型的建立

3.1 模拟气孔构造裂隙几何模型

CFD数值模拟基于Fluent软件进行。玄武岩裂隙面气孔分布极不规则，无法对气孔进行真实模拟，现场采集回的玄武岩试样表面气孔之间未发现相互连通。因此，假设模拟气孔构造裂隙渗流模型中同一裂隙面上的气孔之间无水力联系，将裂隙面的气孔等效为大小相同、分布均匀的柱状孔洞，依据表1所列玄武岩裂隙气孔特征参数建立模拟气孔构造裂隙几何模型，如图6所示。

图6 模拟气孔构造裂隙几何模型Fig.6 Geometry model of vesicular structure offracture

图6显示的裂隙几何模型由平板和圆柱体组成，平板为模型裂隙，圆柱体代表裂隙面的气孔结构，两者接触面连通。几何模型中平板与圆柱体的外表面为边界面，水流从平板短边的一端流入，向另一端流出。

3.2 湍流模型

对于模拟气孔构造裂隙渗流过程中的湍流运动，引入湍流模型进行计算。Fluent中内置多种湍流模型，本文选择常用的可实现k-e模型进行分析。可实现k-e模型为湍流黏性提供了计算公式，为耗散率增加了新的传输方程，能够更加准确地描述平板绕流、对旋转流动、逆压梯度的边界层流动、流动分离等渗流现象，且对于复杂的二次流可取得较好的计算效果，其湍动能k及耗散率e输运公式为：

(8)

(9)

其中：

C1=max[0.43,η/(η+5)],

(10)

η=Sk/ε。

(11)

式中：μi为流体运动黏滞系数；xi、xj均为方向坐标；μt为湍流黏滞系数；Gk表示因平均速度梯度引起的湍动能产生项；Gb为由浮力影响引起湍动能产生项；ε为耗散率；YM表示可压缩湍流脉动膨胀对总耗散的影响；C1ε、C2ε、C3ε为常数；σk和σε分别表示为湍动能及其耗散率湍流普朗特数，分别取1.0、1.2；S为常数项；C2、υ为常数项。

可实现k-e模型与标准k-e模型的相比，其黏性系数Cμ不是常数，通过式(12)计算得到。

(12)

其中，Sij=(∂uj/∂xi+∂ui/∂xj)/2。

Cμ是平均应变率与旋度的函数，在平衡边界层惯性底层可得到Cμ=0.09，与标准k-e模型中采用的常数相同。

3.3 模型校核

将建立的模拟气孔构造裂隙几何模型划分网格后导入Fluent软件中计算，模型输入的计算参数依据试验条件确定，水流参数设置中，测得试验过程水温为12 ℃，水流密度为999.5 kg/m3，黏滞系数通过计算可得为1.24810-3Pa·s，水流入口处的压强通过试验水头差计算获得。数值计算后可获得几何模型中各位置的流速，图7给出了b=0.406 mm、J=0.9时的光滑裂隙和气孔构造裂隙的流场流速云图(取几何模型中心截面)，可见两种裂隙的最大流速、流速分布均有显著差异。

图7 光滑裂隙和气孔构造裂隙水流速度云图Fig.7 Contours of flow velocity in smooth fractureand vesicular structured fracture

为确认数值模拟选取参数的合理性，首先对气孔构造裂隙模型(H=15 mm、d=20 mm、w=18.85%)在不同等效裂隙宽度(b=0.406、0.933、1.643 mm)、不同水头梯度J条件下进行模拟计算，将计算结果的过流流量Q与水头梯度J关系和试验结果一并绘制在图8中。由图8可见，数值模拟结果与试验结果接近，说明数值模拟选用参数是合理的，可进行进一步的模拟计算。

图8 过流流量与水头梯度关系的数值模拟与试验结果对比Fig.8 Comparison of flow rate versus head gradientcurves between numerical simulation andexperimental results

4 渗流场数值模拟结果分析

4.1 经验公式的适用范围

限于条件，本文试验仅开展了Re在6.54～773.11范围内的室内裂隙渗透试验，对于更大范围的雷诺数Re模拟通过数值计算完成。在模拟中，完成了Re=0.002 5～2 000 (不同裂隙宽度、水头梯度、气孔构造特征) 的计算，结果绘制在图9中。

图9 扩大范围后的雷诺数Re与裂隙面构造影响因子λ的变化关系Fig.9 Variation of Reynolds number (Re) againstinfluence factor (λ) with expanded range

由图9可见，Re>200后，数值模拟的结果普遍大于试验结果，更接近水平线λ=1，说明Re>200后数值模拟中气孔的存在对水流的影响没有试验室条件下显著。这是因为试验的玄武岩试块不仅有尺寸较大的气孔，而且裂隙表面还存在岩石本身的微观起伏，气孔代表了岩石的粗构造，这是数值计算可以模拟的。而岩石的微观起伏是细构造，数值模型无法模拟，故Re>200后当数值模拟计算流量会系统性地略大于玄武岩室内试验流量。

考虑数值模拟和试验材料的系统差异后，可以发现，Re-λ关系在Re=2～2000范围内，数值模拟结果展现的渗透规律与式(5)基本相符，因此室内试验和数值模拟取得的规律近似，说明经验公式是有效的，且试验公式的适用范围可推广到Re=2～2 000范围。

4.2 气孔构造对渗流影响的作用机理

选取裂隙面气孔构造参数为H=15 mm、d=20 mm、w=18.85%的裂隙几何模型在不同雷诺数、不同裂隙宽度条件下的模拟结果进行分析，以阐明气孔构造对裂隙渗流的影响机理。

雷诺数、裂隙宽度条件可分为如下4种。

(1)条件1:Re=50、b=0.406 mm。

(2)条件2:Re=600、b=0.406 mm。

(3)条件3:Re=50、b=1.643 mm。

(4)条件4:Re=600、b=1.643 mm。

对4种气孔构造裂隙渗流模型流场矢量图进行对比分析，如图10所示。

图10 气孔构造裂隙水流运动矢量图Fig.10 Vectors of flow movement in vesicularstructured fracture

将图10(a)与图10(b)对比、图10(c)与图10(d)对比可知，气孔构造对裂隙渗流影响的机理为：当水流雷诺数较低时，在气孔中，靠近气孔边壁面水流流速相对裂隙水流流速很小(v孔≈5%v裂隙)，即靠近气孔边壁面水流流动趋于静止，这样，水流与壁面不发生显著碰撞、摩擦，减少了水流能量消耗，这部分水实际起到隔离裂隙运动水流与裂隙壁面的碰撞、摩擦的作用，因此相比一般裂隙，气孔构造促进了裂隙渗流；当水流雷诺数较高时，气孔内的流速较大(v孔>30%v裂隙)，裂隙面气孔内的水流旋涡明显，水流产生的扰动深入到裂隙水流中，可使裂隙水流产生额外扰动，形成更大的能量损失，因此，气孔构造可导致裂隙渗流流量下降。

改变裂隙面气孔体积(H和d发生变化)，建立裂隙面气孔构造参数为H=5 mm、d=5 mm、w=18.85%的裂隙几何模型，将相同水头梯度(J=0.9)的两种裂隙模型进行对比研究(图11)，取裂隙等效水力宽度b=0.406 mm。

v=0.016 3 m/sv=0.015 5 m/s(a) H=15 mmd=20 mm(b) H=5 mmd=5 mm图11 不同裂隙面气孔体积条件下构造裂隙水流运动矢量图Fig.11 Vectors of flow movement in vesicular struc-tured fracture with different vesicular volumes

由图11可见，在裂隙等效水力宽度b=0.406 mm时，两种裂隙渗流都处于紊流状态。在水头梯度相同(J=0.9)的条件下，裂隙面气孔较大的构造裂隙渗流流速大于裂隙面气孔较小的构造裂隙渗流流速，这印证了3#试样裂隙渗流量较1#试样裂隙(较接近光滑裂隙)渗流量更偏离立方定理公式计算流量的结论。

5 结 论

通过选取代表性气孔构造玄武岩试样进行裂隙渗透试验，并建立模拟气孔构造裂隙渗流模型进行分析，得出主要结论如下：

(1)试验结果表明，当等效裂隙宽度较小时，气孔构造岩石裂隙的过流流量大于光滑裂隙渗流流量；当等效裂隙宽度较大时，气孔构造岩石裂隙的过流流量显著小于光滑裂隙渗流流量，说明气孔构造的存在对岩石裂隙渗流有显著影响。

(2)对试验数据进行转换分析后，回归得到了气孔构造裂隙的计算公式。该公式反映的基本规律经数值模拟验证，适用于雷诺数Re∈[2，2 000]。

(3)通过数值模拟，揭示了气孔构造的影响机理：当雷诺数较低(Re<200)时，气孔中的水流隔离了裂隙水流与裂隙壁面的碰撞、摩擦，因此可减少水流阻力，使得渗流流量增大；雷诺数较大(Re>200)时，气孔中的水流产生的扰动深入到裂隙水流中，可使裂隙水流产生额外扰动，形成更大的能量损失，导致渗流量下降。