单排水孔轴对称填充流场解析解

张 帅，叶 飞*，胡 超，郑 双，符文熹

(1.四川大学 水利水电学院，四川 成都 610065；2.四川大学 水力学与山区河流保护国家重点实验室，四川 成都 610065)

岩体是由裂隙、溶隙、填充介质和岩石基质组成的不连续介质。由于岩石基质自身的渗透特性较低，分布于天然岩体中的裂隙、溶隙、填充介质起着主要导水作用。在碳酸岩层中，地表水通过裂隙、溶隙进入岩层储水溶腔，补给地下水体。工程中，为了降低岩体内部渗透压力，通常会布置排水孔进行释压排水，然而随着时间推移，碎屑物质和部分化学结晶物会堆积在排水孔内部，影响排水孔泄水功能。开展部分充填排水孔的流场理论研究是解决坝基渗漏、隧道涌（突）水和油气开采等问题的关键。

前人通过数值试验和物理试验对岩体内部排水裂隙、孔洞流体运动进行大量的研究。在数值研究方面，邵建立等通过COMSOL数值模拟软件对岩体孔隙-裂隙双重介质渗流过程进行了数值模拟，得到了不同形状路基下的速度场和压力场。王志良等通过格子玻尔兹曼方法研究了单裂隙面渗流的流场特征，发现渗流特征主要受裂隙表面形态影响，且粗糙裂隙渗流规律与立方定律误差较大。黄震等对裂隙非饱和渗流进行简化，利用COMSOL Multiphysics分析了裂隙对渗流的影响，并基于溶洞-裂隙系统建立岩溶区裂隙导水突水模型。薛娈鸾将岩体内部裂隙和排水孔视为高孔隙率的填充物，基于复合单元法提出渗流与应力耦合的算法。

在试验研究方面，李琛亮等通过研制双重介质渗流水力特性试验系统，研究双重介质的水力性态和渗流机制，较好地解决了空隙介质、裂隙介质和边界条件的模拟。速宝玉和詹美礼通过试验研究填充裂隙的渗流特性，发现充填裂隙渗透率主要与充填介质颗粒组成、颗粒粒径、孔隙率和裂隙宽度有关。王鹏飞等通过3D打印技术制作不同粗糙度和裂隙宽度的岩体，通过渗透试验研究了不同围压条件下，粗糙度和裂隙宽度变化对贯通充填裂隙渗流特性的影响。沈振中等研制了光滑岩溶管道与裂隙交叉的渗流试验模型，并在不同水力情况下，通过48组不同岩溶排水管道和裂隙交叉的渗流水力试验得出岩溶管道和裂隙交叉时排水管道汇流量的数学表达式。许欢等通过试验分析边坡水平排水孔的淤堵因素，发现淤堵的主要影响因素是渗流路径、倾角、和泥沙絮凝等原因。

此外，填充介质与自由流体间的界面边界条件是自由流与渗流耦合的关键，大致分为3类：速度滑移界面边界条件，连续界面边界条件，应力跳跃界面边界条件。Beavers和Joseph应用Darcy-Stokes耦合模型推导了速度滑移界面条件下单边充填裂隙速度分布的解析解，且通过试验验证了推导的理论结果。舒付军等在连续边界条件下推导出部分充填平板排水裂隙岩体等效渗透系数表达式，并通过相应试验验证解析解的正确性。李琪等在应力跳跃边界条件下推导出非对称充填平板通道的流速分布解析解，并以前人的研究作为特例验证所推导结果的合理性。

综上所述，目前的岩体排水渗流理论研究大多是基于排水裂隙，而排水孔的渗流大多集中在数值模拟和模型试验上，对于轴对称填充排水孔的渗流场的理论研究较少，尚缺乏完善的含填充排水孔流场流速分布解析解。因此，有必要对轴对称填充排水孔进行研究，确定其流场流速分布解析解。本文采用Navier-Stokes方程控制岩体贯通排水孔中纯水的运动，采用Darcy-Brinkman方程控制轴对称填充介质中渗流水的运动；在应力跳跃界面边界条件下推导出岩体内部部分填充排水孔的流场流速分布解析解，并与数值解对比，验证解析解的合理性。本文结论对揭示含填充排水孔渗流影响因素、完善和发展排水孔渗流理论具有重要意义，同时为分析岩体内部含填充排水孔的导水特性提供参考。

1 理论描述

岩体含填充排水孔如图1所示。图1中，

R

为孔内半径，

R

为模型总半径，

L

为模型

x

方向的长度。填充介质是被视为孔隙率为

n

渗透率为

K

的各向同性介质，本文假设水是不可压缩的牛顿流体，并且在

x

方向是充分发展的层流，忽略进出口的边界的影响。

图 1 岩体内部轴对称填充排水孔模型Fig. 1 Model of axisymmetric filled drainage holes in the rock mass

1.1 排水孔纯水流控制方程

不可压缩纯水流在排水孔中满足连续性方程：

式中，

v

、

v

和

v

分别为排水孔中纯水在

x

、

y

和

z

方向的实际流速，m/s。

排水孔中纯水的控制方程（Navier-Stokes方程）

在

x

方向可以简化为：

式中：

f

为

x

方向上质量力，kg·m/s；ρ为水的密度，kg/m；

p

为

x

方向上的压强，kg/ms；υ为水的运动黏度，m/s。由于孔中纯水是沿

x

方向流动的层流，

y

和

z

方向上的流速为0，所以∂

v

/∂

y

=∂

v

/∂

z

=0，将其代入式（1），可知∂

v

/∂

x

=0，且

v

在

x

方向上是一个常数，因此∂

v

/∂

x

=0。由于水是恒定流且质量力只包含重力项，所以∂

v

/∂

t

=0，

f

=0。将以上条件代入式（2），其简化结果如下：

式中，µ为水的动力黏度，kg/ms，µ= υρ。

由式（3）可知，∂

p

/∂

x

与

x

无关，并且-∂

p

/∂

x

=Δ

P

/

L

，由于管道内部是轴对称的，式（3）中括号里面的项可以在柱坐标中表示为：

将式（4）代入式（3），岩体内部排水孔中纯水的控制方程可见式（5）：

1.2 填充介质渗流控制方程

填充介质中渗流水控制方程满足连续性方程（式（6））和Darcy-Brinkman方程（式（7））：

式中：

v

、

v

和

v

分别为渗流水在

x

、

y

和

z

方向上的真实流速，m/s；

n

为填充介质孔隙率；

K

为填充介质渗透率，m。

渗流控制方程简化过程与纯水流控制方程相似，填充介质中渗流水控制方程最终简化结果如式（8）所示：

1.3 流场分布解析解

为了便于求解，引入以下无量纲参数：

式中：ξ为相对半径；

M

为黏度比；µ为填充介质有效黏度，kg/ms；

Da

为达西数；

S

为填充介质形状参数；

U

为无量纲速度；

G

为压强在

x

方向变化速率，

G

=-∂

p

/∂

x

=Δ

P

/

L

。

将无量纲参数（式（9））代入式（5）和（8），控制方程无量纲形式如下：

式中，

U

和

U

为纯水流和渗流的无量纲速度。通过求解式（10）和（11），得到

U

和

U

的解析解：

式中，

A

、

A

、

A

和

A

为待求系数，

I

和

K

为修正贝塞尔函数。

如图1所示，岩体内部部分填充排水孔满足以下边界条件：

1）在孔中心（ξ=ξ=0），纯水流流速达到最大值：

2）在填充介质和纯水交界面（ξ=ξ=γ=

R

/

R

），满足流速连续剪应力跳跃边界条件：

3）在填充介质和不透水岩石界面（ξ=ξ=1），流速等于0：

将以上边界条件（式（14）～（17））代入式（12）和（13），可求出待求系数

A

、

A

、

A

和

A

的理论表达式，具体结果可见式（18）～（21）：

式（19）～（21）中：

此外，对式（12）和（13）进行积分还可求得排水孔中纯水流和渗流的流量表达式：

式中，

Q

为排水孔中纯水流的无量纲流量，

Q

为排水孔中填充介质中渗流的无量纲流量。

另外，纯水流和填充介质界面处水流拖曳力可由Newton内摩擦定律求得：

式中，τ为界面拖曳力，Pa。

2 有效性检验

通过4阶Runge-Kutta方法，结合打靶法，求得岩体内部部分充填排水孔流场流速分布数值解（步长为0.05），数值解和解析解的对比结果见图2。在排水孔纯水流流动的自由流体区，两者结果吻合良好。然而，在水渗流的填充介质区，靠近交界面处的数值解略大于解析解。

图2 流速分布数值解和解析解对比Fig. 2 Comparison of velocity distribution between analytical solution and numerical solution

误差产生的原因是数值解忽视界面应力跳跃导致的（β = 0）。在实际问题中存在应力跳跃的现象，在界面处颗粒的缝隙间（图3中A点），流体间的应力是连续的。但在颗粒表面上（图3中B点），流体间应力会发生跳跃。由于目前应力跳跃系数并没有系统的求解方式，Ochoa-Tapia和Whitaker通过对试验结果拟合得出β的范围在-1.0～1.5，在讨论部分将分析β对模型的影响。总体上本模型的流速分布解析解与数值解相差较小，由此证明模型解析解是有效的。

图3 界面应力跳跃Fig. 3 Stress jump on the interface

此外，在相同工况下，将本文解析解计算结果与Poulikakos和Kazmierczak模型计算结果对比，见图4和5。结果发现Poulikakos和Kazmierczak的解析解是本模型当应力跳跃系数β = 0的一个特例。

图4 Da变化情况下的流速分布Fig. 4 Velocity distribution for representative values of Da

图5 γ变化情况下的流速分布Fig. 5 Velocity distribution for representative values ofγ

3 分析与讨论

根据流速分布推导结果可知，流场流速分布主要受4个典型参数的控制：黏度比（

M

）、应力跳跃系数（β）、达西数（

Da

）和空隙相对开度（γ）。

3.1 黏度比M对流速分布U的影响

Givler和Altobelli通过试验发现

M

的范围在1.0～7.5。Al-Azmi通过分析Vafai和Kim提出的模型，发现在速度和应力连续的界面边界条件下，

M

的变化对流场流速变化影响很小。很明显，图6和7中的速度随着

M

的增加而减小，这是因为

M

的增加会使填充介质的孔隙率和渗透率降低，导致颗粒对流体的阻力增大。此外，当β=1.5时，

M

的变化对流速分布的影响很小，这与Al-Azmi的结论相同（图6）；然而，当β=-0.8时，

M

的变化对流速分布的影响很大（图7）。

图6 M变化对流速分布的影响(β=1.5)Fig. 6 Efect of M on velocity distribution(β=1.5)

图7 M变化对流速分布的影响(β=-0.8)Fig. 7 Effect of M on velocity distribution(β=-0.8)

图8为

M

和β的变化对界面流速的影响。β越小，

M

的变化对界面速度的影响越大；同时，当β<0时，

M

的变化会导致界面速度明显的变化，说明

M

对实际流体速度分布的影响不可忽视。

图8 M变化对界面流速的影响Fig. 8 Effect of M on the interface velocity

3.2 应力跳跃系数β对流速分布U的影响

在纯水流与填充介质的交界面处应力是连续的，可取β=0。然而，界面处应力并不总是连续的，Ochoa-Tapia和Whitaker的研究表明β的范围在-1.0～1.5。

图9为β变化对速度分布的影响。从图9可以看出，随着β的减小，流场速度明显增大。当β从1.5减小到-0.8时，速度的增长率逐渐增大，且β < 0时速度变化更为显著。Ochoa-Tapia和Whitaker发现β与填充介质的渗透率呈负相关。随着β的增加，填充介质对流体的阻力增大，流场速度减小，这与Kuznetsov的结论相似。

图9 β变化对流速分布的影响Fig. 9 Effect of β on velocity distribution

3.3 达西数Da对流速分布U的影响

Da

对流速分布的影响见图10，流场流速随着

Da

的增加而增加。这是因为

Da

与填充介质的渗透率成正比，

Da

的增加会使颗粒对流体的阻力减小，能量损失减小。当

Da

<10时，填充介质中的渗流流速非常小，填充介质几乎可视为不透水，此时Darcy-Brinkman公式不再适用，排水孔中纯水流满足泊肃叶定律；当

Da

=1时，填充介质可视为全透水，整个模型中流体可全视为纯水流。

图10 Da变化对流速分布的影响Fig. 10 Effect of Da on velocity distribution

从上述分析可知，β和

Da

都间接反映了填充介质的渗透率，图11中进一步分析了β和

Da

对界面流速的影响。从图11可知，界面流速与

Da

呈正相关，与β呈负相关。然而，在β=-0.8时存在异常点，此时

Da

=1时的界面速度低于

Da

=0.1时的界面速度，这与

Da

与界面速度呈正相关的结论并不完全吻合。这是因为β的取值范围是Ochoa-Tapia和Whitaker通过Beavers和Joseph的试验数据所得到的拟合结果，β的拟合值范围与实际结果有一定的偏差。因此，当β=-0.8时，拟合结果的偏差可能会导致界面速度分析结果有轻微的差异，且只有当β<-0.5和

Da

>10时才会出现轻微的差异。

图11 Da变化对界面流速的影响Fig. 11 Effect of Da on the interface velocity

3.4 相对开度γ对流速分布U的影响

图12为γ变化对流速分布的影响，且γ越大，流场速度越大。γ变化对流速分布的影响分为两个方面：一方面，根据质量守恒定律，流量相同的情况下，γ增大，过流面积增大，理论上相应的流速应该减小；另一方面，当孔隙开度增大时，多孔介质阻力会减小，理论上流速会增大。

图12 γ变化对流速分布的影响Fig. 12 Effect of γ on velocity distribution

图13和14进一步研究了不同β和

Da

值下γ对界面速度的影响。

图13 β变化对界面流速的影响Fig. 13 Effect of β on the interface velocity

图13中，界面速度随γ的增大先增大后减小。此外，界面速度的转折点发生在不同的γ值下。当β>-0.8时，对应的转折点是γ=0.8；而当β=-0.8时，对应的转折点是γ=0.7。图14中，当

Da

=1时，排水孔内流体可以视为纯水流，填充介质对流体的阻力可以忽略，界面速度随着γ的增加而减小。当

Da

=0.1时，速度先增大后减小，与图13中界面速度变化情况相似；当

Da

<0.1时，界面速度随γ的增加而增加，这与图12中界面速度的变化相似。

图14 Da变化对流速分布的影响Fig. 14 Effect of Da on the interface velocity

4 本文模型与传统泊肃叶定律对比

传统泊肃叶定律仅考虑水流在孔洞区域中流动，忽视了填充介质区域的渗流。当填充介质渗透率较大时，填充区域的渗水量不可忽视。

本文同时考虑非填充区域的纯水流和填充区域的渗流。由图15可知，在解释含充填排水孔的水力特性时，本文理论模型比传统的理论模型更合理。

图15 本模型与传统泊肃叶定律对比Fig. 15 Comparison between present model and Poiseuille’s law

5 结 论

联立Navier-Stokes方程和Darcy-Brinkman方程，在应力跳跃边界条件下推导出了自由流-渗流相耦合的岩体内部轴对称填充排水孔流场流速解析解。分析典型参数对流速分布解析解的影响，具体结果如下：

1）流场流速与黏度比（

M

）呈负相关。当应力跳跃系数β>0时，

M

对0流速分布的影响很小；然而当β<0时，

M

对流速分布的影响十分显著，且不可忽略。2）流场流速与应力跳跃系数（β）呈负相关。β增加，填充介质对流体的阻力增加，这会导致流场流速减小。当β<0，典型参数（

M

、

Da

和γ）的变化对流速分布的影响更显著。3）流场流速与达西数（

Da

）呈正相关。当

Da

<10时，可认为填充介质区域不透水，Darcy-Brinkman方程不再适用。当

Da

=1时，填充介质对流体的阻力可忽略不计，通道内的速度分布与泊肃叶流相似。4）流场流速与空隙相对开度（γ）呈正相关。然而，在不同的

Da

和β值下，界面速度并不总是与γ呈正相关。随着γ的逐渐增加，当

Da

=1时，界面速度逐渐减小；当

Da

=0.1时，界面速度先增大后减小；当

Da

<0.1时，界面速度逐渐增大。

5）本模型同时考虑孔洞纯水流和填充介质中渗流，比传统泊肃叶定律更合理。模型理论解可为分析部分填充排水孔泄流能力提供参考。