考虑接触区域空间分布的岩体裂隙渗透特性数值模拟研究

Numerical Study on the Permeability Characteristics of Rock Fracture Considering the Spatial Distribution of Contact Area

摘要

为研究岩体裂隙接触区域空间分布对其渗透特性的影响，采用随机参数法构建理想接触的裂隙结构，并利用接触率和缺项对接触区域空间分布状态进行表征。基于格子Boltzmann方法，采用半反弹格式描述裂隙渗流过程中流固耦合作用机制，建立模拟理想接触岩体裂隙渗流演化机制的数值计算模型；结合接触裂隙渗流模型和网格精度验证该模型的有效性和计算精度，并讨论接触率、接触空间分布和接触几何形状等因素对岩体裂隙渗透特性的影响。结果表明：岩体裂隙接触率越大，其过流空间越小，流体绕流路径越长，使渗流流速变缓，从而降低了其渗透率；截面平均流速、截面平均压降与裂隙过流面积成负相关关系。当裂隙接触区域空间分布较为离散时，其对流体的阻碍作用加剧，渗流流线的曲折度增大，流体流速减小，裂隙的过流能力降低；在接触率为11.77%的情况下，缺项从1.42增加至4.76时，渗透率降低50.3%；当椭圆柱接触的纵横轴比较小时，其几何形貌呈扁平细长状，接触区域的奇异性增大，对流体流动的阻碍作用更加显著，纵横轴比从1.0减小到0.2，裂隙内平均流速从1.49×10^-4降至1.08×10^-4 m/s。研究成果能为岩体接触裂隙渗透特性的定量评价提供重要理论支撑。

Abstract

Objective Accurately understanding the permeability characteristics of rock contact fractures is crucial for the safety evaluation of deep underground engineering. A numerical model for simulating the ideal contact fracture seepage process is proposed based on the lattice Boltzmann method to analyze the evolution of permeability in fractures under localized contact areas, considering the coupling effect between the fluid and contact areas. Methods This study applied a random parameter method to construct an idealized contact fracture structure, and the contact ratio and lacunarity were selected to characterize the spatial distribution of contact areas. Based on the lattice Boltzmann method, a second-order accurate half-bounce-back scheme was utilized to describe the coupling mechanism between the fluid and localized contact areas, and a numerical model was proposed to simulate the seepage process in ideal contacted rock fractures. The effectiveness and computational accuracy of the model were validated using two classic examples. Streamline tortuosity and permeability were introduced to characterize the flow morphology, and the effects of fracture contact ratio, spatial distribution of contact areas, and contact geometry on fracture permeability were investigated to analyze the impact of contact areas on fluid flow in rock fractures. Results and Discussions Under the same driving pressure, as the contact ratio within the rock fracture increased, the fluid flow space became more compressed, the seepage path became more tortuous, the flow velocity decreased, and the fracture permeability was reduced. When the contact ratios were 3.95%, 8.28%, 11.77%, 15.56%, and 20.37%, the average seepage velocities were 1.82×10^-4, 1.54×10^-4, 1.43×10^-4, 1.38×10^-4, and 1.25×10^-4 m/s, respectively. As the contact ratio increased from 3.95% to 20.37%, the fracture permeability decreased from 4.54×10^-8 m² to 3.14×10^-8 m². The spatial distribution of contact areas significantly altered the permeability characteristics of the fractures. At the same contact ratio, a more concentrated distribution of contact areas, indicating a higher lacunarity, resulted in shorter average flow paths around obstacles and a more uniform fluid velocity distribution. The average seepage velocity decreased with decreasing lacunarity. When the lacunarity Λ values were 4.76, 3.21, 2.30, 1.70, and 1.42, the corresponding fracture permeabilities were 4.78×10^-8, 4.15×10^-8, 3.66×10^-8, 3.47×10^-8, and 3.18×10^-8 m², respectively. As the lacunarity increased from 1.42 to 1.70, 2.30, 3.21, and 4.76, the fracture permeability increased by 9.1%, 15.1%, 30.5%, and 50.3%, respectively. The shape of the contact area was also an important factor affecting the permeability characteristics. At a contact ratio C=11.77%, with the contact center position and the inclination angle of the elliptical major axis unchanged, an increase in the aspect ratio of the contact area resulted in a smoother geometric shape, a smaller blockage area in the flow direction, enhanced flow capacity, and a higher average seepage velocity, which collectively led to increased fracture permeability. When the aspect ratios were 0.2, 0.4, 0.6, 0.8, and 1.0, the corresponding average seepage velocities were 1.08×10^-4, 1.38×10^-4, 1.46×10^-4, 1.48×10^-4, and 1.49×10^-4 m/s, respectively, and the permeabilities were 2.71×10^-8, 3.46×10^-8, 3.65×10^-8, 3.70×10^-8, and 3.71×10^-8 m², respectively. Conclusions This study developed a numerical model to simulate seepage evolution in ideally contacted rock fractures, effectively capturing fluid flow patterns within contact fractures and revealing the influence of local contact on fracture permeability. The main conclusions are as follows: as the contact ratio of a rock fracture increases, the available flow space decreases, leading to a longer fluid flow path, which slows the seepage velocity and reduces permeability. In addition, the flow space within the fracture is negatively correlated with both the average flow velocity and the average pressure drop across the cross-section. When the spatial distribution of contact areas is more dispersed, obstruction to fluid flow intensifies, increasing the tortuosity of the seepage path, reducing fluid velocity, and weakening the fracture flow capacity. When the contact ratio is 11.77%, permeability decreases by 50.3% as the missing item increases from 1.42 to 4.76. When the aspect ratio of elliptical contact areas is small, the geometric shape becomes flat and slender, increasing the singularity of the contact areas and significantly obstructing fluid flow. These results provide essential theoretical support for the quantitative evaluation of permeability characteristics in rock contact fractures.

Graphical abstract

关键词

岩体裂隙 / 渗透特性 / 接触率 / 缺项 / 格子Boltzmann方法

Key words

rock fracture / permeability characteristics / contact ratio / lacunarity / lattice Boltzmann method

本刊网刊

深部地下工程中，页岩气开采^[1]、地下水资源开发^[2]、核废料处理^[3]及CO₂地质封存^[4]等，均涉及岩体裂隙的渗流问题。由于受到沉积历史、应力作用等影响，裂隙内部会因存在填充物或岩体碎块而产生接触区域，接触区域的空间分布表现出强随机性，这使得裂隙内的流体运动状态非常复杂，理论模型计算准确性大大降低。为此，研究接触区域空间分布对岩体裂隙渗流特性的影响，对于地下水资源开采、地下能源开发等具有重要意义。

在岩体裂隙渗流特性演化方面，学者开展了大量的研究，并在试验研究、理论分析和数值模拟方面取得了丰硕的成果。在试验研究方面，Chen等^[5]考虑裂隙几何特征影响开展了岩体裂隙渗流试验，并建立了渗透率与裂隙几何特征间的关系。Wang等^[6]基于剪切流动可视化试验装置，研究了剪切过程中接触面积和裂隙开度等因素对粗糙裂缝非线性渗流特性的影响。Li等^[7]根据岩体裂隙剪切渗流试验，提出了评价接触面积和表面粗糙度对岩体裂隙渗流影响的经验公式。Li等^[8]通过裂隙渗流试验，分析了法向应力作用下分形维数和接触率对非线性渗流行为的影响。由于试验条件的限制，目前只针对特定的试验环境研究宏观渗流问题，无法精确获取裂隙壁面接触面积、接触空间分布及接触几何形状等因素引起渗流流速、渗流压降等流动特性的实时演化。在理论分析方面，Walsh^[9]基于Maxwell有效介质方法研究了裂隙中单圆柱状障碍物对渗流流动的影响，并提出了接触率与水力开度的关系。Zimmerman^[10‒11]等研究了圆柱状、椭圆柱状和不规则粗糙形状接触裂隙的渗透特性，并讨论了流体在粗糙岩体裂隙中的流动行为。然而，接触区域在裂隙内的空间分布具有较大的随机性，使得其渗流形态非常复杂。基于理想假设条件建立的理论模型难以反映裂隙内流体的真实流动状况。数值模拟方法因其具有计算成本低、计算效率高、能够模拟实际环境且可视化效果好等优点，在计算岩体裂隙渗流方面得到广泛应用。Wei等^[12]基于有限单元法考虑两个圆柱状接触对裂隙渗流特性的影响，研究了接触区域角度、间距和半径与裂隙渗流流速、局部压强间的联系。Koyama等^[13]将接触区域处理为0孔径单元，采用有限单元法研究了接触裂隙产生的绕流效应。Xiong等^[14]基于有限体积法对三维糙壁裂隙进行了渗流场的数值计算，认为接触区域会增加速度分布的复杂性。基于连续介质理论的传统数值计算方法（有限单元法、有限体积法等），虽然计算过程方便快捷，但无法精准处理裂隙渗流过程中的流固耦合作用，对于裂隙结构复杂性与流体之间相互运动的反映效果较差^[15]。

格子Boltzmann方法，能够反映宏观物理量与微观分子运动间的联系^[16]，具有物理背景清晰^[17]、计算方法简便^[18]和复杂边界便于处理^[19]等优点。在考虑粗糙壁面几何形貌影响下模拟粗糙岩体裂隙渗流方面得到了广泛应用^[20‒21]，目前的研究成果主要揭示裂隙壁面粗糙程度与非线性渗流行为间的联系，尚无法诠释裂隙接触区域大小、空间分布形态及接触几何形状等因素与渗透特性间的关联性。

鉴于此，利用随机参数法构建随机分布的三维理想接触裂隙结构，基于格子Boltzmann方法，采用半反弹格式实现渗流流体与裂隙壁面、接触区域间的流固耦合作用，建立模拟接触裂隙渗流过程的数值计算模型，并讨论裂隙接触率、接触空间分布及接触几何形状等因素对裂隙渗透特性的影响。

1 物理模型

假定岩体裂隙渗流处于层流状态，且流体为不可压缩的牛顿流体，则流体在裂隙内的运动满足：

1）动量守恒方程

\rho \frac{\partial \mathit{u}}{\partial t}+\rho (\mathit{u}\nabla )\mathit{u}=-\nabla P+\nabla (\nabla {\upsilon }_{\mathrm{p}}\mathit{u})（1）

式中，ρ为密度， u 为宏观速度矢量，t为时间，P为压力，{\upsilon }_{\mathrm{p}}为流体的运动黏度，\nabla为哈密尔顿算子。

2）质量守恒方程

\frac{\partial \rho }{\partial t}+\nabla (\rho \mathit{u})=\mathrm{0}（2）

流体在裂隙中流动时，接触区域的阻碍作用会改变其运动路径，从而影响整体的渗透特性。在计算过程中，将岩体假定为刚体，不考虑其变形及接触区域的运动，假设流体与裂隙上下壁面及接触区域之间为无滑移边界条件。

2 格子Boltzmann模型

2.1 格子Boltzmann方程

基于D3Q19模型^[22]（三维空间， e₀～ e₁₈为19个离散速度矢量）进行计算，如图1所示。

单松弛的格子Boltzmann方程（不含外力项）可表示为：

\begin{array}{c}{f}_{\alpha }(\mathit{r}+{\mathit{e}}_{\alpha }\mathrm{\Delta }t,t+\mathrm{\Delta }t)=f_{\alpha }(\mathit{r},t)+\frac{\mathrm{1}}{{\tau }_{\mathrm{f}}}[f_{\alpha }^{\mathrm{e}\mathrm{q}}(\mathit{r},t)-f_{\alpha }(\mathit{r},t)],\\ \alpha =\mathrm{0,1},\cdots ,\mathrm{18}\end{array}（3）

式中：f_α ( r,t)、f_α^eq( r,t)分别为空间位置 r 处t时刻α方向的粒子分布函数和平衡态分布函数，Δt为时间步长，τ_f为无量纲松弛时间， e_α 为α方向上的速度矢量。

对于D3Q19模型，19个方向的速度矢量集合 e 为：

\mathit{e}=c\left[\begin{array}{ccccccccccccccccccc}\mathrm{0}& \mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& -\mathrm{1}& \mathrm{1}& -\mathrm{1}& \mathrm{1}& -\mathrm{1}& \mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{1}& -\mathrm{1}& -\mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& -\mathrm{1}& \mathrm{1}& -\mathrm{1}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& -\mathrm{1}& -\mathrm{1}& \mathrm{1}& \mathrm{1}& -\mathrm{1}& -\mathrm{1}& \mathrm{1}\end{array}\right]（4）

式中，c为格子速度，c=Δx/Δt，其中，Δx为格子步长。

平衡态分布函数f_α^eq( r,t)可表示为：

f_{\alpha }^{\mathrm{e}\mathrm{q}}(\mathit{r},t)=\rho {\omega }_{\alpha }\left[\mathrm{1}+\frac{{\mathit{e}}_{\alpha }\mathit{u}}{c_{\mathrm{s}}^{\mathrm{2}}}+\frac{({\mathit{e}}_{\alpha }{\mathit{u})}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{s}}^{\mathrm{4}}}-\frac{{\mathit{u}}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{s}}^{\mathrm{2}}}\right]（5）

式中：c_s为格子声速，c{}_{\mathrm{s}}^{\mathrm{2}}=c²/3；ω_α 为α方向上的权重系数，其取值为：

{\omega }_{\alpha }=\left\{\begin{array}{l}\mathrm{1}/\mathrm{3},\alpha =\mathrm{0};\\ \mathrm{1}/\mathrm{18},\alpha =\mathrm{1,2},\dots ,\mathrm{6};\\ \mathrm{1}/\mathrm{36},\alpha =\mathrm{7,8},\dots ,\mathrm{18}\end{array}\right.（6）

将格子Boltzmann方程用Chapman‒Enskog展开，可推导出流体的动量守恒及质量守恒方程。同时，得到流体密度、速度、压力、运动黏度为：

\left\{\begin{array}{l}\rho ={\displaystyle \sum _{\alpha =\mathrm{0}}^{\mathrm{18}}}{f}_{\alpha },\\ \mathit{u}=\frac{\mathrm{1}}{\rho }{\displaystyle \sum _{\alpha =\mathrm{0}}^{\mathrm{18}}}{f}_{\alpha }{\mathit{e}}_{\alpha },\\ P=\rho c_{\mathrm{s}}^{\mathrm{2}},\\ {\upsilon }_{\mathrm{p}}=c_{\mathrm{s}}^{\mathrm{2}}\left({\tau }_{\mathrm{f}}-\frac{\mathrm{1}}{\mathrm{2}}\right)\mathrm{\Delta }t\end{array}\right.（7）

2.2 边界处理

采用Zhang等^[23]提出用具有2阶精度的半反弹格式来处理渗流流体与固体壁面、接触区域间的相互作用，其演化方程为：

f_{-\mathit{i}}({\mathit{x}}_{\mathrm{l}},t+\mathrm{\Delta }t)=f_{\mathit{i}}({\mathit{x}}_{\mathrm{l}},t)（8）

式中， x₁为与固体壁面相邻的流体节点， i 为流体节点指向固体壁面方向，- i 表示其相反方向。

对于压力驱动的裂隙出入口边界则采用Guo等^[24]提出的非平衡外推格式：将边界节点 x_s处的粒子分布函数分为平衡态与非平衡态，其平衡态部分可由式（5）求得，而非平衡态部分则由相邻流体节点 x_sf的分布函数代替：

f_{\alpha }({\mathit{x}}_{\mathrm{s}},t)=f_{\alpha }^{\mathrm{e}\mathrm{q}}({\mathit{x}}_{\mathrm{s}},t)+[f_{\alpha }({\mathit{x}}_{\mathrm{s}\mathrm{f}},t)-f_{\alpha }^{\mathrm{e}\mathrm{q}}({\mathit{x}}_{\mathrm{s}\mathrm{f}},t)]（9）

2.3 单位转换

对于物理单位与格子单位之间的转换，可通过定义转换因子得到。首先，确定长度转换因子T_l 、运动黏度转换因子T_υ 和质量转换因子T_m ：

\left\{\begin{array}{l}{T}_l=\frac{l_{\mathrm{P}}}{l_{\mathrm{L}}},\\ T_{\mathrm{\upsilon }}=\frac{{\upsilon }_{\mathrm{P}}}{{\upsilon }_{\mathrm{L}}},\\ T_m=\frac{m_{\mathrm{P}}}{m_{\mathrm{L}}}\end{array}\right.（10）

式中，l_P、m_P为物理单位的长度和质量，l_L、υ_L、m_L为对应格子单位的长度、运动黏度和质量。

确定上述3个转换因子后，可以通过单位间的量纲分析，进而确定其他物理量的转换因子，如表1所示。

3 三维理想接触裂隙的生成及表征

岩体裂隙接触区域的空间分布具有强烈随机性，其对流体运动行为的影响非常复杂。为了突出岩体裂隙接触区域大小、空间分布状态及接触区域形状等因素对其渗透特性的影响，将上下裂隙面简化为光滑壁面，接触区域的形状统一设为圆柱状或椭圆柱状。

3.1 三维理想接触裂隙的生成

为了模拟岩体裂隙接触区域的随机分布，采用随机参数法构建了三维理想接触裂隙，其实施步骤如下：

1）根据接触区域数量，采用随机函数在计算域内生成接触中心坐标(x_i,y_i )，其中，x_i 、y_i 分别服从区间[x_min,x_max]、[y_min,y_max]随机分布，x_min、x_max、y_min、y_max根据计算域的范围确定。

2）对于圆柱状接触，通过随机函数在生成接触半径范围[D_min,D_max]内，生成接触半径D_i；同理，对于椭圆柱状接触，则需采用随机函数确定横轴a_i 、纵轴b_i 及横轴倾角θ_i。

3）根据接触中心坐标(x_i,y_i )与圆柱或椭圆柱几何形状，控制两个接触间不能产生相交区域，然后将接触面在垂直方向延伸生成空间接触柱。

4）重复步骤1）～3），直至接触区域满足设置要求，即完成三维理想接触裂隙构建。

3.2 三维理想接触裂隙的表征

为了衡量裂隙接触程度，定义接触率C为：

C=\frac{C_{\mathrm{c}}}{C_{\mathrm{a}}}\times \mathrm{100}\mathrm{\%}（11）

式中，C_c为接触面积，C_a为裂隙面面积。

缺项是表征裂隙接触区域空间分布状态的重要参数^[25]。缺项值越大，表明接触区域分布越集中；反之则分布越均匀。对于绝对均匀分布的物体，其缺项值为1。

本文基于滑移盒计数法来计算裂隙接触的缺项值，其计算过程如下：

1）确定滑移盒尺寸。选择一个适当的、边长为r的正方形盒子用于计算缺项，r_{\mathrm{e}}\le r\ll \mathrm{m}\mathrm{i}\mathrm{n}(N_x,N_y)^[26]，其中，r_e为单元尺度，N_x 和N_y 为计算域在x和y方向的长和宽。

2）移动盒子。将选定的盒子以固定步长在计算域上平移。遍历计算域的所有位置，且盒子之间无重叠。

3）统计盒计数。对于每个滑动的盒子位置，统计盒子内的接触点数量。

4）缺项计算。根据盒子内统计接触点的数量，计算缺项为：

Q(M,r)=\frac{n(M,r)}{N(r)}（12）

Z_Q^{(q)}(r)=\sum _M{M}^{q}Q(M,r)（13）

\Lambda =\frac{Z_Q^{(\mathrm{2})}(r)}{{\left[Z_Q^{(\mathrm{1})}(r)\right]}^{\mathrm{2}}}（14）

式（12）～（14）中：Q(M,r)为概率密度函数，描述边长为r的盒子中具有M个接触单元的盒子数量n(M,r)占总滑移盒次数N(r)的比例；Z_Q^{(q)}(r)为统计动差函数，q为正整数系数；Λ为缺项，定义为当q=2时的统计动差函数与当q=1时的统计动差函数的平方之比。

3.3 三维理想接触裂隙实例

将计算域的长宽高设置为N_x ×N_y ×N_z =200×200×20的网格，不同接触率模型计算参数如表2所示。根据表2的参数，利用所提出的随机参数法生成了接触率分别为8.28%、11.77%和15.56%的圆柱状接触裂隙模型，图2为不同接触率的接触裂隙截面，图3为接触率C=11.77%的三维接触裂隙。

不同空间分布模型计算参数如表3所示。为了构建不同空间分布状态的接触裂隙，根据表3的设计参数生成了接触率为11.77%，缺项Λ分别为3.21、2.30和1.70的3种裂隙结构，不同空间分布的接触裂隙截面如图4所示。

4 模型验证

4.1 接触裂隙渗流模型验证

为了验证所提格子Boltzmann模型在处理接触裂隙渗流方面的有效性，分别建立了单个圆柱状和椭圆柱状接触的理想裂隙结构，其模型如图5所示。采用的格子单位计算模型为N_x ×N_y ×N_z =300×300×50格子，上下边界设置一层固体区域，采用半反弹边界处理固液界面，故裂隙机械开度h₀=47，出入口两端采用压差驱动。

由文献[9‒10]可知，单个圆柱状和椭圆柱状接触裂隙的水力开度h_圆柱状、h_椭圆柱状的解析解分别为：

h_{\mathrm{圆柱}\mathrm{状}}^{\mathrm{3}}=h_{\mathrm{0}}^{\mathrm{3}}\left(\frac{\mathrm{1}-C}{\mathrm{1}+C}\right)（15）

h_{\mathrm{椭圆}\mathrm{柱状}}^{\mathrm{3}}=h_{\mathrm{0}}^{\mathrm{3}}\left(\frac{\mathrm{1}-\beta C}{\mathrm{1}+\beta C}\right)（16）

式（15）、（16）中：h₀为机械开度；β为与椭圆纵横轴之比n相关的参数，β=(1+n)²/4n。

在计算域中心位置构建了接触率为5%、10%、15%、20%、25%的圆柱状接触及接触率为10%，纵横轴比n为1.0、1.3、2.0、3.0、4.0的椭圆柱状接触裂隙。同时，基于格子Boltzmann方法，通过流量计算了渗流裂隙的水力开度，圆形与椭圆接触裂隙水力开度的数值解与解析解对比如图6所示。由图6可知，无论是圆柱状接触还是椭圆柱状接触裂隙，本文计算得到的水力开度数值解与解析解^[9‒10]均具有非常高的吻合度，这充分证明了本文计算模型在处理接触裂隙渗流方面的正确性。

4.2 网格精度验证

为了验证网格精度对计算结果的影响，针对光滑平板Poiseuille流进行了数值计算。计算区域的长宽高L×W×H为10 mm×10 mm×1 mm，格子步长Δx分别设为0.200、0.100、0.050和0.033 mm。驱动压差ΔP=0.01 Pa，流体运动黏度υ=1.0×10^-6 m²/s。

表4为不同网格精度计算结果。由表4可知：当Δx=0.100和0.200 mm时，计算网格数较少，计算耗时较短，但计算精度偏低；当Δx=0.033 mm时，计算精度有了极大的提升，但计算耗时又显著增长；当Δx=0.050 mm时，计算耗时相对较短，且具有较高的计算精度，故在后续的数值计算中选择Δx=0.050 mm。

5 分析讨论

为研究岩体裂隙接触区域对其渗流特性的影响，考虑接触率、接触空间分布及接触几何形状等因素影响，构建理想接触裂隙模型，基于格子Boltzmann方法，研究岩体接触裂隙的渗流特性。计算模型L×W×H=10 mm×10 mm×1 mm，计算网格为200 mm×200 mm×20 mm，h₀=0.85 mm，流体的运动黏度υ=1.0×10^-6 m²/s。研究重点是空间分布对裂隙渗透特性的影响，故将驱动压差设置为较低水平，避免复杂流态的产生。计算模型如图7所示。

为了便于分析岩体裂隙中接触区域对流体运动行为的影响，引入曲折度η²表征流线的弯曲程度^[27]：

{\eta }^{\mathrm{2}}={\left(\frac{L_{\mathrm{a}}}{L}\right)}^{\mathrm{2}}（17）

式中，L_a为流线长度，为了减少流线长度的计算误差，针对整个渗流场选用90条以上的贯通流线，并将其平均值作为流线长度。

将渗透率K作为评价裂隙的渗透特性变化的依据，其表达式为：

K=\frac{Q\rho {\upsilon }_{\mathrm{p}}L}{\mathrm{\Delta }PA}（18）

式中，Q为裂隙截面流量，A为截面面积。

5.1 接触率对裂隙渗透特性的影响

为讨论接触率C对裂隙渗透特性的影响，计算C为3.95%、8.28%、11.77%、15.56%和20.37%的理想接触裂隙，z=0.5 mm截面渗流流线分布如图8所示。图8中，u_x 为流场内x方向上的流速。接触区域的存在会改变流体的运动路径，从而影响其渗流场的速度分布。在接触区域附近，流体受到阻碍作用，原有的流动形态破坏，流速的大小、方向及分布均发生变化，并产生了局部阻力^[28]。流体越靠近接触区域，其流速越小，在远离接触壁面处会形成流速高速区，与Wei等^[12]观察到的现象一致。对比分析发现，在相同驱动压力作用下，岩体裂隙内的接触率越高，流体的过流空间受到压缩，渗流流线更加曲折，渗流流速越小，流体的运动形态更加复杂。当C为3.95%、8.28%、11.77%、15.56%和20.37%时，流体平均渗流流速分别为1.82×10^-4、1.54×10^-4、1.43×10^-4、1.38×10^-4、1.25×10^-4 m/s。

为进一步分析接触率对流体运动形态的影响，建立了流线曲折度η²与接触率C的关系，如图9所示。

由图9可知，流线曲折度η²随着接触率C的增加而增大，且两者间近似呈线性关系。在低接触率C=3.95%情况下，流体渗流较为顺畅，流线较平滑，曲折度较低（图8（a））。随着接触率的增加，流体绕流路径逐渐增长，流线变得愈发曲折，曲折度显著增加（图8（e））。

图10为岩体裂隙渗透率随着接触率的变化趋势。由图10可知，当接触率为3.95%时，裂隙渗透率最大为4.54×10^-8 m²。随着接触率的增大，裂隙渗透率的降低趋势偏离线性关系，这说明即使在理想接触情况下，裂隙接触率并不是影响渗透率的唯一因素。随着接触率增大，流线曲折度增大，流体流动路径变长，不同流层间的摩擦阻力增大，流体在裂隙中的整体流动速度变缓，进而降低裂隙的渗透性。

图11为不同接触率情况下裂隙过流面积A_cs与截面平均流速u_avg的关系。

由图11可知，裂隙过流面积与截面平均流速呈负相关关系，过流面积A_cs的减少对水流有加速作用，且这种作用随接触率C的增大更为显著。原因是接触区域占据了裂隙内原有的过流空间，使得过流空间压缩，被接触区域阻挡的流体涌入接触区域之间的空隙，使其承受更多的流体通过而导致流速变大。

图12为不同接触率情况下裂隙过流面积A_cs与截面平均压降ΔP_avg的关系。由图12可知，裂隙截面平均压降随过流面积的增大而减小。这主要是由于裂隙内压力势能与流体动能间存在相互转化关系，裂隙空隙内渗流流速增加的动力来自压力势能的减小。随着裂隙接触率的增大，压降主要发生在过流面积较小的截面，此时接触区域面积较大，速度梯度较大，黏滞力作用效果显著，压降增大明显。

5.2 接触空间分布对裂隙渗透特性的影响

为讨论接触区域空间分布对接触裂隙渗流的影响，对C=11.77%和缺项Λ分别为4.76、3.21、2.30、1.70和1.42的裂隙结构进行渗流流动的数值计算，不同空间分布情况下z=0.5 mm截面的渗流流线分布如图13所示。由图13可知，在相同接触率和驱动压力作用下，裂隙接触区域分布越集中（缺项越大），平均绕流流线越短，流体流速分布越均匀，高流速区域的数量少但面积大，且裂隙平均渗流流速随缺项的减小而减小。当缺项Λ为4.76、3.21、2.30、1.70及1.42时，流体平均渗流流速分别为1.91×10^-4、1.66×10^-4、1.47×10^-4、1.39×10^-4、1.27×10^-4 m/s。

图14为不同缺项裂隙出口处的流速分布。由图14可知，裂隙接触区域的空间分布改变了其出口处的流速分布。当接触区域较为集中（缺项较大）时，裂隙出口的流速大且分布均匀。随着缺项的减小，出口流速的波动现象加剧。这是因为接触率相同的情况下，接触区域的空间分布不同，其与流体的接触面积也有所差异。当缺项Λ为4.76、2.30及1.42时，两者接触面积分别为10.34、32.71、46.46 mm²。接触面积越大，其对流体的阻碍作用越显著，出口处的平均流速也越小。当接触较为均匀（缺项较小）时，在接触区域的影响下，裂隙内形成了多条优势渗流通道，导致出口处流速分布出现波动。

为进一步分析接触空间分布对流体渗流流态的影响，建立了流线曲折度η²与缺项Λ的关系，如图15所示。由图15可知，随着缺项Λ的增大，流线曲折度η²呈逐渐减小的趋势。当缺项较大时，接触区域集中分布，只有靠近接触区域的流体出现明显的绕流现象，其余流线仍较为顺直，流线的曲折度较低（图13（a））。随着缺项的减小，接触区域分布逐渐趋于均匀化，阻碍流体运动的接触区域数量增多，流体与固体的接触面积增加，流线曲折度也逐渐增大（图13（e））。而流线曲折度越大，其水力路径越长，这将导致流体间黏滞阻力的作用路径延长，整体流动阻力增大。同时，流线曲折度的增大，意味着流速大小、方向的变化加剧，导致局部阻力增大^[29]，并进一步影响其渗透特性。

图16为岩体裂隙渗透率随缺项的变化趋势。由图16可知，在接触率相同的情况下，缺项越大，岩体裂隙内的渗透率越高。当缺项从1.42增加至1.70、2.30、3.21及4.76时，裂隙渗透率增幅分别为9.1%、15.1%、30.5%、50.3%。原因是较大的缺项产生较小的流线曲折度，从而导致流体间和流固之间的摩擦阻力减小，产生较小的局部水头损失，进而使裂隙的渗透率增大。文献[30]的计算结果也得到了类似的演化规律，这表明接触区域的空间分布是影响裂隙渗透率的重要因素。

5.3 接触几何形状对裂隙渗透特性的影响

接触区域的几何形状会改变裂隙内流体流动的路径，从而影响其渗透特性。为此，当C=11.77%时，接触中心位置和椭圆横轴倾角不变的情况下，通过改变椭圆柱纵横轴比值n构建了不同接触几何形状的裂隙，并开展了其渗流流动的数值计算，不同纵横轴比情况下z=0.5 mm截面渗流流线分布如图17所示。

在接触率、接触空间位置及驱动压力相同的情况下，接触区域的n越大，其几何形状越圆滑，沿来流方向上接触区域的阻流面积越小，裂隙的过流能力越强，平均渗流流速也越大。当n为0.2、0.4、0.6、0.8和1.0时，裂隙平均渗流速度分别为1.08×10^-4、1.38×10^-4、1.46×10^-4、1.48×10^-4、1.49×10^-4 m/s。

裂隙渗流的流线曲折度η²与纵横轴比n的关系如图18所示。

由图18可知，η²随n的增加而减小，且变化趋势逐渐趋于平缓。当n较小（n=0.2）时，接触区域呈细长状，在来流方向上对过水断面的压缩作用显著，形成较窄的过流通道，流线渗流路径变长，流线的曲折度增加，进而导致流体运动的阻碍增加（图17（a））。随着n的增大，接触区域趋于圆滑，不同倾角情况下阻流面积的奇异性减小，其对流体的阻碍作用变小，渗流流线的曲折度也相应地减小。

岩体裂隙的渗透率K与纵横轴比n的关系如图19所示。

由图19可知，在接触率、接触空间分布位置相同的情况下，n越大，裂隙渗透率越高，且随n的增大，渗透率的变化趋于平缓。椭圆横纵比越大，接触区域形状越圆滑，单个接触区域的阻流奇异性越小，导致整体裂隙渗流的阻流程度减弱，裂隙渗透率增加。

6 结 论

基于格子Boltzmann方法，采用半反弹格式处理流固耦合作用边界，建立了模拟理想接触裂隙渗流流动过程的数值计算模型，并讨论了接触率、接触空间分布及接触几何形状等因素对接触裂隙渗流特性的影响，得到了以下结论：

1）基于所提出的数值计算模型，针对单个圆柱状和椭圆柱状的理想接触裂隙进行了渗流流动的数值计算，水力开度的计算结果与解析解具有较高的吻合度，说明该模型在解决岩体接触裂隙渗流问题的有效性。

2）裂隙接触率越大，其内部的过流空间越小，渗流流态越复杂，流线曲折度也越大，从而导致裂隙渗流流动速度变缓，渗透率减小。同时，截面平均流速、截面平均压降与裂隙过流面积成负相关关系。过流面积的减少对水流有加速作用，对压降有增强作用。

3）在裂隙接触率相同的情况下，接触区域分布越均匀，阻碍流体流动的接触区域数量越多，流体与固体的接触面积增加，其渗流流线的曲折度增大，进而导致局部阻力增加，渗流流速减小，裂隙渗透率降低。

4）裂隙内椭圆柱状接触的纵横轴比越小，其对流体渗流阻碍作用的奇异性越大，沿来流方向上阻流面积越大，导致渗流流线曲折度与局部阻力均增大，进而降低裂隙的渗透率。

参考文献

[1] Wang Shihao, Zhang Yanbin, Wu Haiyi,et al.A kinetic mo-del for multicomponent gas transport in shale gas reserv-oirs and its applications[J].Physics of Fluids,2022,34(8):082002. doi:10.1063/5.0101272

[2] Berkowitz B.Characterizing flow and transport in fractu-red geological media:A review[J].Advances in Water Resources,2002,25(8/9/10/11/12):861‒884. doi:10.1016/s0309-1708(02)00042-8

[3] Tsang C F, Bernier F, Davies C.Geohydromechanical proc-esses in the Excavation Damaged Zone in crystalline rock,rock salt,and indurated and plastic clays—In the context of radioactive waste disposal[J].International Journal of Rock Mechanics and Mining Sciences,2005,42(1):109‒125. doi:10.1016/j.ijrmms.2004.08.003

[4] Zhao Yan, Yang Liu, Xi Ruru,et al.CO₂‒H₂O two-phase displacement characteristics of low permeability core using nuclear magnetic resonance and magnetic resonance imaging techniques[J].Rock and Soil Mechanics,2023,44(6):1636‒1644.

[5] 赵艳,杨柳,奚茹茹,等.基于核磁共振和磁共振成像的低渗透岩芯CO₂‒H₂O两相驱替特征研究[J].岩土力学,2023,44(6):1636‒1644.

[6] Chen Yuedu, Liang Weiguo, Lian Haojie,et al.Experimental study on the effect of fracture geometric characteristics on the permeability in deformable rough-walled fractures[J].International Journal of Rock Mechanics and Mining Sciences,2017,98:121‒140. doi:10.1016/j.ijrmms.2017.07.003

[7] Wang Changsheng, Liu Richeng, Jiang Yujing,et al.Effect of shear-induced contact area and aperture variations on nonlinear flow behaviors in fractal rock fractures[J].Journal of Rock Mechanics and Geotechnical Engineering,2023,15(2):309‒322. doi:10.1016/j.jrmge.2022.04.014

[8] Li Bo, Jiang Yujing, Koyama T,et al.Experimental study of the hydro-mechanical behavior of rock joints using a parallel-plate model containing contact areas and artificial fractures[J].International Journal of Rock Mechanics and Mining Sciences,2008,45(3):362‒375. doi:10.1016/j.ijrmms.2007.06.004

[9] Li Man, Liu Xianshan, Li Yu,et al.Effect of contact areas on seepage behavior in rough fractures under normal stress[J].International Journal of Geomechanics,2022,22(4):04022019. doi:10.1061/(asce)gm.1943-5622.0002330

[10] Walsh J B.Effect of pore pressure and confining pressure on fracture permeability[J].International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abst-racts,1981,18(5):429‒435. doi:10.1016/0148-9062(81)90006-1

[11] Zimmerman R W, Chen Diwen, Cook N G W.The effect of contact area on the permeability of fractures[J].Journal of Hydrology,1992,139(1/2/3/4):79‒96. doi:10.1016/0022-1694(92)90196-3

[12] Zimmerman R W, Bodvarsson G S.Hydraulic conductivity of rock fractures[J].Transport in Porous Media,1996,23(1):1‒30. doi:10.1007/bf00145263

[13] Wei Xianfa, Ma Haichun, Qian Jiazhong,et al.Study on the geometric characteristics effect of contact area on fracture seepage[J].Physics of Fluids,2023,35:016603. doi:10.1063/5.0131145

[14] Koyama T, Li B, Jiang Y,et al.Numerical modelling of fluid flow tests in a rock fracture with a special algorithm for contact areas[J].Computers and Geotechnics,2009,36(1/2):291‒303. doi:10.1016/j.compgeo.2008.02.010

[15] Xiong Feng, Jiang Qinghui, Ye Zuyang,et al.Nonlinear flow behavior through rough-walled rock fractures:The effect of contact area[J].Computers and Geotechnics,2018,102:179‒195. doi:10.1016/j.compgeo.2018.06.006

[16] Tan Yunliang, Yin Yanchun, Teng Guirong,et al.Simulation research of gas seepage based on Lattice Boltzmann method[J].Journal of China Coal Society,2014,39(8):1446‒1454.

[17] 谭云亮,尹延春,滕桂荣,等.基于Lattice Boltzmann方法的瓦斯渗流模拟研究[J].煤炭学报,2014,39(8):1446‒1454.

[18] Jin Lei, Zeng Yawu, Cheng Tao,et al.Seepage characteristics of soil-rock mixture based on lattice Boltzmann met-hod[J].Chinese Journal of Geotechnical Engineering,2022,44(4):669‒677.

[19] 金磊,曾亚武,程涛,等.基于格子Boltzm-ann方法的土石混合体的渗流特性研究[J].岩土工程学报,2022,44(4):669‒677.

[20] Huang Mengmeng, Zhang Buyang, Lu Yao,et al.LBM ba-sed simulation analyses of the scattered droplet deposition in printed OLEDs[J].Advanced Engineering Sciences,2023,55(6):54‒65.

[21] 黄萌萌,张不扬,鲁瑶,等.基于格子Boltzmann的喷印OLED散点墨滴沉积仿真分析[J].工程科学与技术,2023,55(6):54‒65.

[22] Ma Jiuchen, Cui Afeng, Linhai Lyu,et al.Flow characteristics of suspended particles in aquifer using coupled LBM‒DEM numerical model[J].Advanced Engineering Scien-ces,2023,55(5):149‒160.

[23] 马玖辰,崔阿凤,吕林海,等.基于LBM‒DEM耦合计算模型的含水层内悬浮颗粒运动特性[J].工程科学与技术,2023,55(5):149‒160.

[24] Pu Dan, Li Miao, Shen Linfang,et al.The effects of channel width on particle sedimentation in fluids using a coupled lattice Boltzmann-discrete element model[J].Physics of Fl-uids,2023,35(5):053307. doi:10.1063/5.0147826

[25] Wang Min, Chen Yifeng, Ma Guowei,et al.Influence of surface roughness on nonlinear flow behaviors in 3D self-affine rough fractures:Lattice Boltzmann simulations[J].Advances in Water Resources,2016,96:373‒388. doi:10.1016/j.advwatres.2016.08.006

[26] Ma Guowei, Ma Chunlei, Chen Yun.An investigation of nolinear flow behaviour along rough-walled fractures co-nsidering the effects of fractal dimensions and contact areas[J].Journal of Natural Gas Science and Engineering,2022,104:104675. doi:10.1016/j.jngse.2022.104675

[27] Qian Y H, D'Humières D, Lallemand P.Lattice BGK models for navier-stokes equation[J].Europhysics Letters(EPL),1992,17(6):479‒484. doi:10.1209/0295-5075/17/6/001

[28] Zhang Ting, Shi Baochang, Guo Zhaoli,et al.General bo-unce-back scheme for concentration boundary condition in the lattice-Boltzmann method[J].Physical Review E,2012,85:016701. doi:10.1103/physreve.88.029903

[29] Guo Zhaoli, Zheng Chuguang, Shi Baochang.Non-equilib-rium extrapolation method for velocity and pressure bou-ndary conditions in the lattice Boltzmann method[J].Chinese Physics,2002,11(4):366‒374. doi:10.1088/1009-1963/11/4/310

[30] Xia Yuxuan, Cai Jianchao, Perfect E,et al.Fractal dimension,lacunarity and succolarity analyses on CT images of reservoir rocks for permeability prediction[J].Journal of Hydrology,2019,579:124198. doi:10.1016/j.jhydrol.2019.124198

[31] Allain C, Cloitre M.Characterizing the lacunarity of random and deterministic fractal sets[J].Physical Review A,1991,44(6):3552‒3558. doi:10.1103/physreva.44.3552

[32] Brown S, Caprihan A, Hardy R.Experimental observation of fluid flow channels in a single fracture[J].Journal of Geophysical Research:Solid Earth,1998,103(B3):5125‒5132. doi:10.1029/97jb03542

[33] Finnemore E J, Franzini J B.Fluid Mechanics with Engineering Applications[M].New York:McGraw‒Hill Companies Inc,2002.

[34] Ju Yang, Zhang Qingang, Yang Yongming,et al.An experimental investigation on the mechanism of fluid flow th-rough single rough fracture of rock[J].Science China Tec-hnological Sciences,2013,56(8):2070‒2080. doi:10.1007/s11431-013-5274-6

[35] Pan Pengzhi, Feng Xiating, Xu Dingping,et al.Modelling fluid flow through a single fracture with different contacts using cellular automata[J].Computers and Geotechnics,2011,38(8):959‒969. doi:10.1016/j.compgeo.2011.07.002

基金资助

国家自然科学基金项目(42167022)

国家自然科学基金项目(42067043)

国家自然科学基金项目(41931294)