考虑壁面浸润性的光滑岩体微裂隙渗流特性数值模拟研究

申林方，李腾风，王志良，李 泽，王鹏宇

(昆明理工大学建筑工程学院，昆明，云南 650500)

岩体在长期的自然风化、地质构造变动及人类活动等作用下，形成了大量平均隙宽仅为0.2 mm的微裂隙[1]。这些微裂隙构成了地下水运移的主要通道，影响着岩体的整体渗透特性。而岩体微裂隙的渗流问题广泛存在于实际工程中[2-3]，在油气开采时人们希望提高岩体渗透率以增加开采率；而在水利、地下空间开发、核废弃物处置等工程建设中，又希望降低其渗透性增加结构安全。因此研究岩体微裂隙的渗流机制，具有非常重要的理论和工程应用价值。

由于岩体微裂隙的隙宽微小，固液界面表面积与液体体积之比较大，使得固液壁面间的浸润性显著影响流体的运动特性。与宏观裂隙相比，微裂隙具有渗透性弱、持水性强、导水性差、隐蔽性强等特点。同时，岩石通常由多种矿物组成，而不同矿物成分的浸润性存在着显著差异，如：文象花岗岩中石英的水接触角约为144.0°、萤石的水接触角约为74.6°等[4]，从而进一步加剧了岩体微裂隙渗流的复杂性。因此，有必要考虑壁面浸润性影响对岩体微裂隙的渗流机制展开研究。

在微尺度下，流体受到裂隙壁面浸润性的影响，导致宏观条件下的无滑移边界不再适用[5]。针对于微尺度渗流特性研究，国内外学者在实验研究、理论研究及数值模拟等方面，已进行了大量的研究工作，并取得了较好的进展。在实验研究方面：Holt 等[6]通过实验研究发现当水通过微尺度碳纳米管时，其流速比Poiseuille 流理论解高3 个数量级。Song 等[7]采用实验方法研究了超疏水表面毫米级气-水界面的滑移速度。霍素斌等[8]基于水在超疏水和超亲水微通道内的渗流实验，指出超疏水微通道内的流动阻力显著降低，最大可达25%。总体而言，采用实验方法研究虽然可以直观得到微尺度下的渗流规律，但其花费高、耗时长，且实验结果往往受控于实验方法、仪器精度等[9]。在理论研究方面：Navier 最早注意到固体壁面附近的滑移现象，并认为壁面处流体的滑移速度与其剪切速率成正比[10]。Spikes 等[11]通过引入滑移模型，建立了新的描述牛顿流体运动的数学方程。Yang 等[12]研究了流体在疏水壁面矩形微通道内的流动，并提出了一种改进方法用以确定疏水平行板微通道的滑移系数。Gennes[13]根据剪应力与摩擦阻力间的关系，推导出疏水壁面滑移长度的计算公式。对于微通道渗流的理论研究，通常其基本假定过于理想化，所得结论只能用于定性分析，与实际情况尚存在一定的偏差。随着计算机技术的快速发展，数值计算方法在研究微通道内的流体流动方面得到了广泛应用。Martell 等[14]基于直接数值模拟方法(DNS)，研究了超疏水表面在微通道中的减阻性能，并讨论了不同粗糙度影响下的滑移速度、壁面剪应力及雷诺应力等。刘赵淼和逄燕[15]模拟了液体在圆形和梯形截面微通道内的流动，并分析了不同压力差、微通道尺寸和表面粗糙度等因素对流动摩擦系数的影响。杨大勇和刘莹[16]基于有限单元法，研究了微通道表面粗糙度对其电渗流特性的影响。然而，常规的宏观数值计算方法，无法真实的反映微通道固液壁面间的作用机制，考虑壁面浸润性影响从微观的角度揭示其渗流特征。

基于分子动理论发展起来的Boltzmann 方法(LBM)，是一种宏观离散、微观连续的介观数值计算方法。它通过考虑不同相不同组分间分子作用力的影响，描述流体与流体以及流体与固体之间的微观作用机制，能够弥补宏观数值计算方法的缺陷，从微观角度出发解释宏观物理现象，因此在微通道多相流体界面动力学的研究中得到了广泛应用。付宇航等[17]基于伪势格子Boltzmann 方法模拟了倾斜壁面浸润性梯度驱动液滴的运动过程。黄桥高等[18]采用格子Boltzmann 方法研究了壁面浸润性对疏水表面滑移流动及减阻特性的影响规律。Genty 等[19]基于格子Boltzmann 方法，研究了深部非饱和泥岩微裂隙内部的有效扩散系数。然而考虑壁面浸润性影响的多相流问题分布于不同的研究领域，现有成果相对分散且各自研究的侧重点也有所不同，尚不足以全面揭示岩体微裂隙的渗流机制。鉴于此，本文基于Shan-Chen伪势模型的格子Boltzmann 方法，考虑壁面浸润性的影响建立了光滑岩体微裂隙渗流的数值模型，结合两个经典算例验证了计算模型的有效性，并讨论了裂隙壁面浸润性、裂隙隙宽、压力梯度及流体黏滞性等因素对其渗流特性的影响。

1 壁面滑移边界

当流体在裂隙中流动时，经典的宏观流体力学认为壁面附近流体速度为零，即无滑移边界。随着裂隙隙宽的减小，固液壁面接触面积与流体体积之比随之增大，导致壁面与流体间的相互作用不可忽略，从而促使与壁面接触的部分流体产生切向速度，发生相对运动，称之为滑移边界[20]。

图1 分别为无滑移、黏滞层、正滑移边界的示意图。其中，us为壁面处的流体流速，在无滑移和黏滞层边界中us=0。当壁面强疏水时，壁面与附近流体之间存在排斥力，此时产生正滑移。反之，当壁面强亲水时，壁面处流体在吸引力作用下，易形成“类固体层”，即黏滞层。

图1 滑移边界类型Fig. 1 Slip boundary type

2 格子Boltzmann 模型

2.1 格子Boltzmann 方程

2.2 边界处理

2.3 单位转换

本文将所有计算参数进行无量纲化处理，并确保无量纲化前后的流动准则数一致，将无量纲参数作为联系的桥梁，实现物理单位与格子单位之间的单位转换。对于本文的渗流问题，引入如下无量纲参数：

2.4 收敛判断依据

2.5 算法验证

2.5.1 蒸汽中的悬浮液滴模拟

基于Shan-Chen 伪势模型的格子Boltzmann 方法，模拟了悬浮在蒸汽中液滴的演化过程。计算区域网格数为128× 128 格子，计算参数采用格子单位，假定初始状态下蒸汽(密度ρg为0.032)中存在一个椭圆形液滴(密度ρl为0.279)，气液相粒子间的相互作用强度G =1/3，液滴的长轴与短轴之比为1.25，如图4 所示，计算模型四周均采用周期性边界。在气液两相分子间作用力的影响下，液滴逐渐由椭圆形演化为圆形，并趋于稳定。将稳定后蒸汽与液滴密度平均值的位置作为气液相界面，并用于计算液滴半径。

根据Laplace 定律，稳定后液滴的内外压差与其半径成反比，可表示为[24]：

式中： pl为液滴内压力； pg为蒸汽中的压力；γgl为气液表面张力；r 为稳定后的液滴半径。

图4 悬浮在蒸汽中的液滴Fig. 4 A suspended liquid droplet in vapor

本文模拟了长轴长度从12 增至40 格子数的情况下，椭圆形液滴的演化过程，并得到了稳定后液滴内外压差与其半径间的关系，如图5 所示。由图可知，圆形液滴内外压差与其半径的倒数呈线性关系，线性拟合曲线关系式为：

拟合曲线的相关系数为 R2=0.9986，这与Laplace 定律相符，也充分证明了计算程序的正确性。

图5 液滴内外压差与半径间的关系Fig. 5 Relationship between pressure difference and droplet radius

2.5.2 壁面接触角模拟

为了研究固体壁面的浸润性，假定液态水在固体壁面形成液滴，基于伪势Shan-Chen 模型的格子Boltzmann 方法，模拟了不同浸润性壁面的接触角。计算模型采用200× 150 的网格系统，计算参数均采用格子单位。初始状态下，在固体壁面的中央设置一个半径r=30 的半圆形液滴，其液相密度为ρl=0.279，液滴周围的气相密度为ρg=0.032，两相粒子间的相互作用强度G =1/3，计算模型四周均采用周期性边界。通过改变虚拟壁面密度ρs从0.06 增至0.18，使不同相间的分子作用力发生变化，从而改变了固体壁面的浸润性。液滴在分子间作用力的影响下，其形态逐渐发展演化，稳定后得到不同浸润性壁面的接触角θ。经计算得到，不同浸润性固体壁面接触角与虚拟壁面密度的对应关系如表2 所示，其典型接触角的几何形态如图6 所示。在本文后续的分析讨论中，以此为基础，通过设置不同的虚拟壁面密度，选取相应的接触角。

表2 接触角与壁面密度对应关系Table 2 Relationship between contact angle and wall density

图6 不同固体壁面接触角的几何形态Fig. 6 Geometric shapes of different solid wall contact angles

3 数值分析研究

为了研究光滑岩体微裂隙的渗流特性，基于Shan-Chen 伪势模型的格子Boltzmann 方法，考虑壁面浸润性影响建立了相应的数值模型。计算模型如图7 所示，微裂隙的L×W=1.6 mm×0.2 mm，划分成800×100 的网格，在初始状态下裂隙中充满液体，且液体的运动黏度ν=1.0×10-6m2/s，裂隙渗流的驱动力为重力，系统中计算参数物理单位与格子单位间的转换，如表1 所示。基于建立的微裂隙渗流数值模型，研究了裂隙壁面浸润性、裂隙隙宽、压力梯度、流体黏滞性等因素对微裂隙渗流特性的影响，并将计算结果与无滑移边界的Poiseuille 流进行了对比分析。

图7 计算模型Fig. 7 Computational model

在压力驱动下，无滑移边界Poiseuille 流平均流速的计算公式为[28]：

3.1 壁面浸润性对微裂隙渗流特性的影响

为了研究壁面浸润性对岩体微裂隙渗流特性的影响，针对裂隙隙宽W=0.2 mm 不同接触角的裂隙方案，进行了渗流场的数值计算，并将计算结果与无滑移边界的解析解进行了对比，其相对误差 δr如表3 所示。由表3 可知，对于疏水壁面，由于其对附近流体的排斥作用而产生加速效果，使得平均渗流流速大于无滑移边界，且随着裂隙壁面接触角的增加(疏水程度增强)，其相对误差逐渐增大。而亲水壁面由于壁面附近流体受到吸引而产生黏滞效应，使得其平均渗流流速小于无滑移边界，且随着接触角的减小(亲水程度增强)，相对误差逐渐增大。但从总体变化趋势看，疏水壁面对微裂隙渗流的影响要比亲水壁面显著。

表3 不同接触角微裂隙的平均渗流流速Table 3 Average velocity in a micro-fracture with different contact angles

3.2 隙宽对微裂隙渗流特性的影响

岩体微裂隙隙宽变化改变了固液界面表面积与流体体积之比，从而进一步影响到壁面浸润性与流体渗流间的关系。为此，针对重力作用下接触角分别为57°、84°、114°和147°的四种裂隙，研究了裂隙宽度从0.08 mm 增加至0.32 mm 时，岩体微裂隙的渗流规律，其平均渗流流速与无滑移边界相对误差的变化规律，如图8 所示。由图可知，在壁面浸润性相同的情况下，随着裂隙隙宽的增加，其渗流流速与无滑移边界的相对误差逐渐减小。说明壁面浸润性对微裂隙渗流特性的影响，随隙宽的增加而逐渐减弱。此外，壁面的亲/疏水性愈强，微裂隙隙宽对其渗流特性的影响愈发显著。接触角为147°、114°、84°及57°的壁面，微裂隙隙宽从0.32 mm 降至0.08 mm 的过程中，其平均流速与无滑移边界的相对误差分别增加了7.42%、2.07%、1.18%及3.19%。

图8 平均流速相对误差随微裂隙隙宽的变化趋势Fig. 8 Relative error variation of average velocity with micro-fracture width

3.3 压力梯度对微裂隙渗流特性的影响

压力梯度以及流体与固体壁面间的作用力，共同影响着岩体微裂隙的渗流特性。为此，保持裂隙隙宽W=0.2 mm 不变，针对接触角为147°、114°、84°及57°的壁面，将压力梯度G*从0.6ρg 逐渐增加至1.8ρg，研究了压力梯度对微裂隙渗流特性的影响。裂隙平均流速与压力梯度的对应关系如图9 所示，由图可知，裂隙的平均渗流流速随压力梯度的增加而增大，且两者呈线性关系，这一规律与无滑移边界Poiseuille 流的理论公式相一致。同时，裂隙壁面疏水能力越强，其直线斜率越大，说明强疏水壁面对微裂隙渗流具有显著的促进作用。对于亲水或弱疏水壁面，其斜率变化并不显著。

图9 微裂隙平均流速与压力梯度的关系Fig. 9 Relationship between average velocity of microfracture and pressure gradient

3.4 流体黏滞性对微裂隙渗流特性的影响

岩体微裂隙中流体的黏滞性，随着温度、溶质成分及浓度等因素的不同而有所差异。为此，在裂隙宽度W=0.2 mm 及压力梯度G*=ρg不变的情况下，分析了流体运动黏度从0.6×10-6m2/s 增加至1.6×10-6m2/s 的过程中，壁面接触角为57°、84°、114°和147°时，裂隙渗流平均流速的变化规律。微裂隙平均渗流流速与流体运动黏度间的关系，如图10 所示。由图可知，随着流体运动黏度的增加，流动阻力增大，微裂隙渗流的平均流速逐渐减小，两者呈反比例关系，这与Poiseuille 流的解析式(26)相符。当流体的运动黏度相同时，裂隙渗流的平均流速随接触角的增加而增大，且对于强疏水壁面其影响尤为显著。

图10 微裂隙平均流速与流体运动黏度的关系Fig. 10 Relationship between average velocity of microfracture and kinematic viscosity of fluid

4 结论

基于Shan -Chen 伪势模型的格子Boltzmann 方法，考虑壁面浸润性的影响建立了光滑岩体微裂隙渗流的数值模型，结合蒸汽中悬浮液滴和壁面接触角的模拟，证明了该模型的准确性和有效性。最后，考虑岩体壁面浸润性、裂隙隙宽、压力梯度及流体黏滞性等因素影响，研究了岩体微裂隙的渗流特性，得出以下结论。

(1)疏水壁面对附近流体的排斥作用，使得微裂隙的平均渗流流速要大于无滑移边界，且随着疏水程度的增强，两者间的相对误差逐渐增大；而亲水壁面对流体的黏滞作用使其产生减速效果。总体而言，疏水壁面对微裂隙渗流的影响比亲水壁面更加显著。

(2)壁面浸润性对光滑岩体微裂隙渗流特性的影响，随着隙宽的减小而逐渐增强，且壁面的亲/疏水性越强，微裂隙隙宽对其渗流特性的影响也越显著。

(3)岩体微裂隙的平均渗流流速随压力梯度的增加而增大，且两者呈线性关系。同时，由于强疏水壁面对微裂隙渗流具有显著的促进作用，导致其斜率较大。

(4)随着流体运动黏度的增加，流动阻力增大，导致微裂隙渗流的平均流速逐渐减小，流体的运动黏度与渗流平均流速呈反比例关系。