張 帥，葉 飛*，胡 超，鄭 雙，符文熹

(1.四川大學(xué) 水利水電學(xué)院，四川 成都 610065；2.四川大學(xué) 水力學(xué)與山區(qū)河流保護(hù)國家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610065)

巖體是由裂隙、溶隙、填充介質(zhì)和巖石基質(zhì)組成的不連續(xù)介質(zhì)。由于巖石基質(zhì)自身的滲透特性較低，分布于天然巖體中的裂隙、溶隙、填充介質(zhì)起著主要導(dǎo)水作用。在碳酸巖層中，地表水通過裂隙、溶隙進(jìn)入巖層儲(chǔ)水溶腔，補(bǔ)給地下水體。工程中，為了降低巖體內(nèi)部滲透壓力，通常會(huì)布置排水孔進(jìn)行釋壓排水，然而隨著時(shí)間推移，碎屑物質(zhì)和部分化學(xué)結(jié)晶物會(huì)堆積在排水孔內(nèi)部，影響排水孔泄水功能。開展部分充填排水孔的流場理論研究是解決壩基滲漏、隧道涌（突）水和油氣開采等問題的關(guān)鍵。

前人通過數(shù)值試驗(yàn)和物理試驗(yàn)對(duì)巖體內(nèi)部排水裂隙、孔洞流體運(yùn)動(dòng)進(jìn)行大量的研究。在數(shù)值研究方面，邵建立等通過COMSOL數(shù)值模擬軟件對(duì)巖體孔隙-裂隙雙重介質(zhì)滲流過程進(jìn)行了數(shù)值模擬，得到了不同形狀路基下的速度場和壓力場。王志良等通過格子玻爾茲曼方法研究了單裂隙面滲流的流場特征，發(fā)現(xiàn)滲流特征主要受裂隙表面形態(tài)影響，且粗糙裂隙滲流規(guī)律與立方定律誤差較大。黃震等對(duì)裂隙非飽和滲流進(jìn)行簡化，利用COMSOL Multiphysics分析了裂隙對(duì)滲流的影響，并基于溶洞-裂隙系統(tǒng)建立巖溶區(qū)裂隙導(dǎo)水突水模型。薛孌鸞將巖體內(nèi)部裂隙和排水孔視為高孔隙率的填充物，基于復(fù)合單元法提出滲流與應(yīng)力耦合的算法。

在試驗(yàn)研究方面，李琛亮等通過研制雙重介質(zhì)滲流水力特性試驗(yàn)系統(tǒng)，研究雙重介質(zhì)的水力性態(tài)和滲流機(jī)制，較好地解決了空隙介質(zhì)、裂隙介質(zhì)和邊界條件的模擬。速寶玉和詹美禮通過試驗(yàn)研究填充裂隙的滲流特性，發(fā)現(xiàn)充填裂隙滲透率主要與充填介質(zhì)顆粒組成、顆粒粒徑、孔隙率和裂隙寬度有關(guān)。王鵬飛等通過3D打印技術(shù)制作不同粗糙度和裂隙寬度的巖體，通過滲透試驗(yàn)研究了不同圍壓條件下，粗糙度和裂隙寬度變化對(duì)貫通充填裂隙滲流特性的影響。沈振中等研制了光滑巖溶管道與裂隙交叉的滲流試驗(yàn)?zāi)Ｐ?，并在不同水力情況下，通過48組不同巖溶排水管道和裂隙交叉的滲流水力試驗(yàn)得出巖溶管道和裂隙交叉時(shí)排水管道匯流量的數(shù)學(xué)表達(dá)式。許歡等通過試驗(yàn)分析邊坡水平排水孔的淤堵因素，發(fā)現(xiàn)淤堵的主要影響因素是滲流路徑、傾角、和泥沙絮凝等原因。

此外，填充介質(zhì)與自由流體間的界面邊界條件是自由流與滲流耦合的關(guān)鍵，大致分為3類：速度滑移界面邊界條件，連續(xù)界面邊界條件，應(yīng)力跳躍界面邊界條件。Beavers和Joseph應(yīng)用Darcy-Stokes耦合模型推導(dǎo)了速度滑移界面條件下單邊充填裂隙速度分布的解析解，且通過試驗(yàn)驗(yàn)證了推導(dǎo)的理論結(jié)果。舒付軍等在連續(xù)邊界條件下推導(dǎo)出部分充填平板排水裂隙巖體等效滲透系數(shù)表達(dá)式，并通過相應(yīng)試驗(yàn)驗(yàn)證解析解的正確性。李琪等在應(yīng)力跳躍邊界條件下推導(dǎo)出非對(duì)稱充填平板通道的流速分布解析解，并以前人的研究作為特例驗(yàn)證所推導(dǎo)結(jié)果的合理性。

綜上所述，目前的巖體排水滲流理論研究大多是基于排水裂隙，而排水孔的滲流大多集中在數(shù)值模擬和模型試驗(yàn)上，對(duì)于軸對(duì)稱填充排水孔的滲流場的理論研究較少，尚缺乏完善的含填充排水孔流場流速分布解析解。因此，有必要對(duì)軸對(duì)稱填充排水孔進(jìn)行研究，確定其流場流速分布解析解。本文采用Navier-Stokes方程控制巖體貫通排水孔中純水的運(yùn)動(dòng)，采用Darcy-Brinkman方程控制軸對(duì)稱填充介質(zhì)中滲流水的運(yùn)動(dòng)；在應(yīng)力跳躍界面邊界條件下推導(dǎo)出巖體內(nèi)部部分填充排水孔的流場流速分布解析解，并與數(shù)值解對(duì)比，驗(yàn)證解析解的合理性。本文結(jié)論對(duì)揭示含填充排水孔滲流影響因素、完善和發(fā)展排水孔滲流理論具有重要意義，同時(shí)為分析巖體內(nèi)部含填充排水孔的導(dǎo)水特性提供參考。

1 理論描述

巖體含填充排水孔如圖1所示。圖1中，

R

為孔內(nèi)半徑，

R

為模型總半徑，

L

為模型

x

方向的長度。填充介質(zhì)是被視為孔隙率為

n

滲透率為

K

的各向同性介質(zhì)，本文假設(shè)水是不可壓縮的牛頓流體，并且在

x

方向是充分發(fā)展的層流，忽略進(jìn)出口的邊界的影響。

圖 1 巖體內(nèi)部軸對(duì)稱填充排水孔模型Fig. 1 Model of axisymmetric filled drainage holes in the rock mass

1.1 排水孔純水流控制方程

不可壓縮純水流在排水孔中滿足連續(xù)性方程：

式中，

v

、

v

和

v

分別為排水孔中純水在

x

、

y

和

z

方向的實(shí)際流速，m/s。

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

在

x

方向可以簡化為：

式中：

f

為

x

方向上質(zhì)量力，kg·m/s；ρ為水的密度，kg/m；

p

為

x

方向上的壓強(qiáng)，kg/ms；υ為水的運(yùn)動(dòng)黏度，m/s。由于孔中純水是沿

x

方向流動(dòng)的層流，

y

和

z

方向上的流速為0，所以?

v

/?

y

=?

v

/?

z

=0，將其代入式（1），可知?

v

/?

x

=0，且

v

在

x

方向上是一個(gè)常數(shù)，因此?

v

/?

x

=0。由于水是恒定流且質(zhì)量力只包含重力項(xiàng)，所以?

v

/?

t

=0，

f

=0。將以上條件代入式（2），其簡化結(jié)果如下：

式中，μ為水的動(dòng)力黏度，kg/ms，μ= υρ。

由式（3）可知，?

p

/?

x

與

x

無關(guān)，并且-?

p

/?

x

=Δ

P

/

L

，由于管道內(nèi)部是軸對(duì)稱的，式（3）中括號(hào)里面的項(xiàng)可以在柱坐標(biāo)中表示為：

將式（4）代入式（3），巖體內(nèi)部排水孔中純水的控制方程可見式（5）：

1.2 填充介質(zhì)滲流控制方程

填充介質(zhì)中滲流水控制方程滿足連續(xù)性方程（式（6））和Darcy-Brinkman方程（式（7））：

式中：

v

、

v

和

v

分別為滲流水在

x

、

y

和

z

方向上的真實(shí)流速，m/s；

n

為填充介質(zhì)孔隙率；

K

為填充介質(zhì)滲透率，m。

滲流控制方程簡化過程與純水流控制方程相似，填充介質(zhì)中滲流水控制方程最終簡化結(jié)果如式（8）所示：

1.3 流場分布解析解

為了便于求解，引入以下無量綱參數(shù)：

式中：ξ為相對(duì)半徑；

M

為黏度比；μ為填充介質(zhì)有效黏度，kg/ms；

Da

為達(dá)西數(shù)；

S

為填充介質(zhì)形狀參數(shù)；

U

為無量綱速度；

G

為壓強(qiáng)在

x

方向變化速率，

G

=-?

p

/?

x

=Δ

P

/

L

。

將無量綱參數(shù)（式（9））代入式（5）和（8），控制方程無量綱形式如下：

式中，

U

和

U

為純水流和滲流的無量綱速度。通過求解式（10）和（11），得到

U

和

U

的解析解：

式中，

A

、

A

、

A

和

A

為待求系數(shù)，

I

和

K

為修正貝塞爾函數(shù)。

如圖1所示，巖體內(nèi)部部分填充排水孔滿足以下邊界條件：

1）在孔中心（ξ=ξ=0），純水流流速達(dá)到最大值：

2）在填充介質(zhì)和純水交界面（ξ=ξ=γ=

R

/

R

），滿足流速連續(xù)剪應(yīng)力跳躍邊界條件：

3）在填充介質(zhì)和不透水巖石界面（ξ=ξ=1），流速等于0：

將以上邊界條件（式（14）～（17））代入式（12）和（13），可求出待求系數(shù)

A

、

A

、

A

和

A

的理論表達(dá)式，具體結(jié)果可見式（18）～（21）：

式（19）～（21）中：

此外，對(duì)式（12）和（13）進(jìn)行積分還可求得排水孔中純水流和滲流的流量表達(dá)式：

式中，

Q

為排水孔中純水流的無量綱流量，

Q

為排水孔中填充介質(zhì)中滲流的無量綱流量。

另外，純水流和填充介質(zhì)界面處水流拖曳力可由Newton內(nèi)摩擦定律求得：

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

2 有效性檢驗(yàn)

通過4階Runge-Kutta方法，結(jié)合打靶法，求得巖體內(nèi)部部分充填排水孔流場流速分布數(shù)值解（步長為0.05），數(shù)值解和解析解的對(duì)比結(jié)果見圖2。在排水孔純水流流動(dòng)的自由流體區(qū)，兩者結(jié)果吻合良好。然而，在水滲流的填充介質(zhì)區(qū)，靠近交界面處的數(shù)值解略大于解析解。

圖2 流速分布數(shù)值解和解析解對(duì)比Fig. 2 Comparison of velocity distribution between analytical solution and numerical solution

誤差產(chǎn)生的原因是數(shù)值解忽視界面應(yīng)力跳躍導(dǎo)致的（β = 0）。在實(shí)際問題中存在應(yīng)力跳躍的現(xiàn)象，在界面處顆粒的縫隙間（圖3中A點(diǎn)），流體間的應(yīng)力是連續(xù)的。但在顆粒表面上（圖3中B點(diǎn)），流體間應(yīng)力會(huì)發(fā)生跳躍。由于目前應(yīng)力跳躍系數(shù)并沒有系統(tǒng)的求解方式，Ochoa-Tapia和Whitaker通過對(duì)試驗(yàn)結(jié)果擬合得出β的范圍在-1.0～1.5，在討論部分將分析β對(duì)模型的影響。總體上本模型的流速分布解析解與數(shù)值解相差較小，由此證明模型解析解是有效的。

圖3 界面應(yīng)力跳躍Fig. 3 Stress jump on the interface

此外，在相同工況下，將本文解析解計(jì)算結(jié)果與Poulikakos和Kazmierczak模型計(jì)算結(jié)果對(duì)比，見圖4和5。結(jié)果發(fā)現(xiàn)Poulikakos和Kazmierczak的解析解是本模型當(dāng)應(yīng)力跳躍系數(shù)β = 0的一個(gè)特例。

圖4 Da變化情況下的流速分布Fig. 4 Velocity distribution for representative values of Da

圖5 γ變化情況下的流速分布Fig. 5 Velocity distribution for representative values ofγ

3 分析與討論

根據(jù)流速分布推導(dǎo)結(jié)果可知，流場流速分布主要受4個(gè)典型參數(shù)的控制：黏度比（

M

）、應(yīng)力跳躍系數(shù)（β）、達(dá)西數(shù)（

Da

）和空隙相對(duì)開度（γ）。

3.1 黏度比M對(duì)流速分布U的影響

Givler和Altobelli通過試驗(yàn)發(fā)現(xiàn)

M

的范圍在1.0～7.5。Al-Azmi通過分析Vafai和Kim提出的模型，發(fā)現(xiàn)在速度和應(yīng)力連續(xù)的界面邊界條件下，

M

的變化對(duì)流場流速變化影響很小。很明顯，圖6和7中的速度隨著

M

的增加而減小，這是因?yàn)?p>M

的增加會(huì)使填充介質(zhì)的孔隙率和滲透率降低，導(dǎo)致顆粒對(duì)流體的阻力增大。此外，當(dāng)β=1.5時(shí)，

M

的變化對(duì)流速分布的影響很小，這與Al-Azmi的結(jié)論相同（圖6）；然而，當(dāng)β=-0.8時(shí)，

M

的變化對(duì)流速分布的影響很大（圖7）。

圖6 M變化對(duì)流速分布的影響(β=1.5)Fig. 6 Efect of M on velocity distribution(β=1.5)

圖7 M變化對(duì)流速分布的影響(β=-0.8)Fig. 7 Effect of M on velocity distribution(β=-0.8)

圖8為

M

和β的變化對(duì)界面流速的影響。β越小，

M

的變化對(duì)界面速度的影響越大；同時(shí)，當(dāng)β<0時(shí)，

M

的變化會(huì)導(dǎo)致界面速度明顯的變化，說明

M

對(duì)實(shí)際流體速度分布的影響不可忽視。

圖8 M變化對(duì)界面流速的影響Fig. 8 Effect of M on the interface velocity

3.2 應(yīng)力跳躍系數(shù)β對(duì)流速分布U的影響

在純水流與填充介質(zhì)的交界面處應(yīng)力是連續(xù)的，可取β=0。然而，界面處應(yīng)力并不總是連續(xù)的，Ochoa-Tapia和Whitaker的研究表明β的范圍在-1.0～1.5。

圖9為β變化對(duì)速度分布的影響。從圖9可以看出，隨著β的減小，流場速度明顯增大。當(dāng)β從1.5減小到-0.8時(shí)，速度的增長率逐漸增大，且β < 0時(shí)速度變化更為顯著。Ochoa-Tapia和Whitaker發(fā)現(xiàn)β與填充介質(zhì)的滲透率呈負(fù)相關(guān)。隨著β的增加，填充介質(zhì)對(duì)流體的阻力增大，流場速度減小，這與Kuznetsov的結(jié)論相似。

圖9 β變化對(duì)流速分布的影響Fig. 9 Effect of β on velocity distribution

3.3 達(dá)西數(shù)Da對(duì)流速分布U的影響

Da

對(duì)流速分布的影響見圖10，流場流速隨著

Da

的增加而增加。這是因?yàn)?p>Da

與填充介質(zhì)的滲透率成正比，

Da

的增加會(huì)使顆粒對(duì)流體的阻力減小，能量損失減小。當(dāng)

Da

<10時(shí)，填充介質(zhì)中的滲流流速非常小，填充介質(zhì)幾乎可視為不透水，此時(shí)Darcy-Brinkman公式不再適用，排水孔中純水流滿足泊肅葉定律；當(dāng)

Da

=1時(shí)，填充介質(zhì)可視為全透水，整個(gè)模型中流體可全視為純水流。

圖10 Da變化對(duì)流速分布的影響Fig. 10 Effect of Da on velocity distribution

從上述分析可知，β和

Da

都間接反映了填充介質(zhì)的滲透率，圖11中進(jìn)一步分析了β和

Da

對(duì)界面流速的影響。從圖11可知，界面流速與

Da

呈正相關(guān)，與β呈負(fù)相關(guān)。然而，在β=-0.8時(shí)存在異常點(diǎn)，此時(shí)

Da

=1時(shí)的界面速度低于

Da

=0.1時(shí)的界面速度，這與

Da

與界面速度呈正相關(guān)的結(jié)論并不完全吻合。這是因?yàn)棣碌娜≈捣秶荗choa-Tapia和Whitaker通過Beavers和Joseph的試驗(yàn)數(shù)據(jù)所得到的擬合結(jié)果，β的擬合值范圍與實(shí)際結(jié)果有一定的偏差。因此，當(dāng)β=-0.8時(shí)，擬合結(jié)果的偏差可能會(huì)導(dǎo)致界面速度分析結(jié)果有輕微的差異，且只有當(dāng)β<-0.5和

Da

>10時(shí)才會(huì)出現(xiàn)輕微的差異。

圖11 Da變化對(duì)界面流速的影響Fig. 11 Effect of Da on the interface velocity

3.4 相對(duì)開度γ對(duì)流速分布U的影響

圖12為γ變化對(duì)流速分布的影響，且γ越大，流場速度越大。γ變化對(duì)流速分布的影響分為兩個(gè)方面：一方面，根據(jù)質(zhì)量守恒定律，流量相同的情況下，γ增大，過流面積增大，理論上相應(yīng)的流速應(yīng)該減小；另一方面，當(dāng)孔隙開度增大時(shí)，多孔介質(zhì)阻力會(huì)減小，理論上流速會(huì)增大。

圖12 γ變化對(duì)流速分布的影響Fig. 12 Effect of γ on velocity distribution

圖13和14進(jìn)一步研究了不同β和

Da

值下γ對(duì)界面速度的影響。

圖13 β變化對(duì)界面流速的影響Fig. 13 Effect of β on the interface velocity

圖13中，界面速度隨γ的增大先增大后減小。此外，界面速度的轉(zhuǎn)折點(diǎn)發(fā)生在不同的γ值下。當(dāng)β>-0.8時(shí)，對(duì)應(yīng)的轉(zhuǎn)折點(diǎn)是γ=0.8；而當(dāng)β=-0.8時(shí)，對(duì)應(yīng)的轉(zhuǎn)折點(diǎn)是γ=0.7。圖14中，當(dāng)

Da

=1時(shí)，排水孔內(nèi)流體可以視為純水流，填充介質(zhì)對(duì)流體的阻力可以忽略，界面速度隨著γ的增加而減小。當(dāng)

Da

=0.1時(shí)，速度先增大后減小，與圖13中界面速度變化情況相似；當(dāng)

Da

<0.1時(shí)，界面速度隨γ的增加而增加，這與圖12中界面速度的變化相似。

圖14 Da變化對(duì)流速分布的影響Fig. 14 Effect of Da on the interface velocity

4 本文模型與傳統(tǒng)泊肅葉定律對(duì)比

傳統(tǒng)泊肅葉定律僅考慮水流在孔洞區(qū)域中流動(dòng)，忽視了填充介質(zhì)區(qū)域的滲流。當(dāng)填充介質(zhì)滲透率較大時(shí)，填充區(qū)域的滲水量不可忽視。

本文同時(shí)考慮非填充區(qū)域的純水流和填充區(qū)域的滲流。由圖15可知，在解釋含充填排水孔的水力特性時(shí)，本文理論模型比傳統(tǒng)的理論模型更合理。

圖15 本模型與傳統(tǒng)泊肅葉定律對(duì)比Fig. 15 Comparison between present model and Poiseuille’s law

5 結(jié) 論

聯(lián)立Navier-Stokes方程和Darcy-Brinkman方程，在應(yīng)力跳躍邊界條件下推導(dǎo)出了自由流-滲流相耦合的巖體內(nèi)部軸對(duì)稱填充排水孔流場流速解析解。分析典型參數(shù)對(duì)流速分布解析解的影響，具體結(jié)果如下：

1）流場流速與黏度比（

M

）呈負(fù)相關(guān)。當(dāng)應(yīng)力跳躍系數(shù)β>0時(shí)，

M

對(duì)0流速分布的影響很??；然而當(dāng)β<0時(shí)，

M

對(duì)流速分布的影響十分顯著，且不可忽略。2）流場流速與應(yīng)力跳躍系數(shù)（β）呈負(fù)相關(guān)。β增加，填充介質(zhì)對(duì)流體的阻力增加，這會(huì)導(dǎo)致流場流速減小。當(dāng)β<0，典型參數(shù)（

M

、

Da

和γ）的變化對(duì)流速分布的影響更顯著。3）流場流速與達(dá)西數(shù)（

Da

）呈正相關(guān)。當(dāng)

Da

<10時(shí)，可認(rèn)為填充介質(zhì)區(qū)域不透水，Darcy-Brinkman方程不再適用。當(dāng)

Da

=1時(shí)，填充介質(zhì)對(duì)流體的阻力可忽略不計(jì)，通道內(nèi)的速度分布與泊肅葉流相似。4）流場流速與空隙相對(duì)開度（γ）呈正相關(guān)。然而，在不同的

Da

和β值下，界面速度并不總是與γ呈正相關(guān)。隨著γ的逐漸增加，當(dāng)

Da

=1時(shí)，界面速度逐漸減??；當(dāng)

Da

=0.1時(shí)，界面速度先增大后減?。划?dāng)

Da

<0.1時(shí)，界面速度逐漸增大。

5）本模型同時(shí)考慮孔洞純水流和填充介質(zhì)中滲流，比傳統(tǒng)泊肅葉定律更合理。模型理論解可為分析部分填充排水孔泄流能力提供參考。