矩形狹縫裂隙水流運動特征數值模擬研究

張帥領,張春艷,江承陽,楊亞東

(1.中國電建集團河南省電力勘測設計院有限公司,河南 鄭州 450007;2.華北水利水電大學 地球科學與工程學院,河南 鄭州 450046)

0 引言

裂隙作為巖溶含水系統的重要組成部分,對巖溶水流運動起到重要作用[1]。裂隙水流運動規律的研究對巖溶水流運動規律、溶質運移規律以及裂隙巖體滲流理論的發展起到重要作用[2-4]。裂隙水流運動一般可以用納維-斯托克斯方程(簡稱NS方程)來描述[5]。但是由于NS方程解的復雜性,限制了方程的實際應用。眾多研究者按照水流運動實際情況,對NS方程進行了不同程度的簡化,第一層將NS方程化為斯托克斯(Stokes)方程;第二層將Stokes方程簡化為雷諾方程;第三層將雷諾方程簡化為立方定律(CL),由于立方定律的簡潔性,得到了廣泛的引用。文章以具有四個裂隙面的矩形狹縫裂隙為研究對象(如圖1),探討立方定律在矩形狹縫裂隙中的適用性以及矩形狹縫裂隙水流運動特性。

圖1 矩形狹縫平行裂隙示意圖

NS方程化為Stokes方程的簡化條件為:慣性力可以忽略,即流速較小,粘性力占主導地位,雷諾數Re較小。研究者主要從雷諾數出發,研究了Stokes方程的適用條件。Sharp[6],Iwai[7]以及Schrauf[8]利用野外裂隙樣本通過室內試驗分析得到,當雷諾數大于1到10的某個數值時,慣性力不可忽略。Zimmerman[9]提出當Re(b)/∧b(波長)小于8時,慣性力對滲流量修正系數的影響最高達10%。與NS方程相比,Stokes方程有了很大程度的簡化,但是三維Stokes方程的解仍非常復雜。

對LCL進行簡化,簡化條件為兩裂隙面光滑且平行分布(即uy=uz=0),則三維的流體流動可概化為一維流動,設流體沿x軸流動,則沿x軸方向的裂隙滲流量為:

Q=Cb3J

(1)

式中:Q為裂隙間的流量,m3·s-1;b為裂隙開度,m;J為水力坡度,無量綱;C為常數,1 s-1。

方程(1)由于流量Q與裂隙開度b的三次方呈正比,故稱為立方定律(Cubic Law,CL)(Boussinesq[16])。對于直線流,常數C可表達為:

(2)

聯合方程(1)和(2),可得裂隙滲流量方程為:

(3)

式中:W為裂隙寬度,m;L為裂隙長度,m;ρ為流體密度,kg·m-3;g為重力加速度,m·s-2;μ為動力黏滯系數,N·s·m-2。

將流量方程Q=Wbv以及水力坡度與水頭損失之間的關系J=hf/L代入式(3),得到立方定律的另一種表達形式:

(4)

式中:hf為水流通過裂隙時產生的水頭損失,m;v為水流流速,m·s-1。

式(4)表明,恒溫條件下,水流流經光滑裂隙時的水頭損失與水流流速之間呈線性關系。

由于立方定律的簡潔性,立方定律得到了廣泛應用[17],但由于立方定律的局限性,國內外許多學者主要從裂隙開度[18-19]、裂隙壁面粗糙性[20]、節理粗糙度系數修正[21]、面積接觸率[22]等方面進行了研究,并基于各自的研究對象提出了立方定律的修正形式[23]。綜合國內外研究成果,絕大多數研究者注重于開放裂隙中立方定律的應用,在實際野外地區,完全開放的裂隙較少存在,矩形狹縫裂隙廣泛存在,針對這一問題,以過水斷面為矩形的狹縫裂隙(如圖1所示)為研究對象,通過數值模擬,對比NS的解,對立方定律的有效性進行研究,提出以極限流速及極限雷諾數來衡量立方定律是否適用,提出適用于矩形狹縫裂隙的修正立方定律,并對裂隙流的水力特性進行了研究,該研究不僅對以往的理論研究進行一定的補充,而且為巖溶含水系統水流運動特征的研究提供理論支撐。

1 矩形狹縫裂隙數值模型的建立

1.1 控制方程

非穩定流條件下,對于不可壓縮流體,NS方程[24]為:

(5)

式中:f為體積力,N;p為壓力,Pa;t為時間,s。

穩定流情況下,對于不可壓縮流體,NS方程為:

ρ(v)v=μ2v-p

(6)

對于兩光滑平板間的流體,其運動符合泊肅葉定律,設裂隙兩側壁分布在z=±b/2處,那么流速分布表達式[25]為:

(7)

1.2 邊界條件

對不同尺寸的裂隙進行二維和三維模擬,不同模擬情景下,邊界條件及裂隙尺寸的設置見表1。

表1 不同模擬情景下邊界條件、裂隙尺寸設置

1.3 模型驗證(S0模擬結果)

對尺寸為W=2 cm,b=0.840 0 mm,L=20 cm的裂隙中的水流運動進行數值模擬,數值模型設置溫度T=13℃,水的密度ρ=999.4 kg/m3,動力粘度μ=0.001 2 Pa·s,重力加速度g=9.8 m/s2,邊界條件的設置見表1中模擬情景S0。

水力坡度試驗測量值與數值模擬值對比如圖2所示。經計算,試驗測量值與數值模擬值相對誤差((數值模擬值-試驗測量值)/試驗測量值)在10%以內。

圖2 試驗測量值與數值模擬值對比圖

2 模型結果與分析

2.1 立方定律適用性驗證(S1模擬結果)

數值模擬得到7個不同裂隙寬度(W=0.5,1,2,3,4,5,6 cm),6個不同裂隙開度(b=0.05,0.06,0.07,0.08,0.09,0.1 cm),共42個不同尺寸的裂隙水流模擬結果。以數值模擬所得水力坡度與立方定律計算所得水力坡度的比值(J-NS/J-CL)為研究對象(J-NS與J-CL對比圖見附件1),來確定立方定律的適用范圍?，F將不同裂隙寬度(W),不同裂隙開度(b)條件下,J-NS/J-CL隨流速變化作對比分析,如圖3所示(以W=1,2 cm為例,其余結果見附件2)。結果表明,流速較小時,J-NS與J-CL重合度較高,表明流速較小時,可用立方定律來描述裂隙水流,然而隨著流速的增大,J-CL與J-NS之間的差值越來越大,且J-CL均小于J-NS,表明流速較大時,立方定律計算值與實際值偏離越大,立方定律不可以再用來描述裂隙水流運動。

圖3 不同裂隙寬度及開度條件下,水力坡度NS值與CL值之比J-NS/J-CL隨流速變化

結果表明,在本次研究中(裂隙寬度分布于0.5～6 cm之間,裂隙開度分布于0.05～0.1 cm之間),J-NS/J-CL分布于0.9～2之間,J-CL與J-NS的差值隨流速的增大而增大,隨裂隙開度的增大而增大,隨裂隙寬度的增大而減小,因為流速越大,慣性力的作用越不可忽略。為確定立方定律的適用范圍,定義相對偏差為((CL計算值-NS方程模擬值)/NS方程模擬值),相對偏差達到10%時對應的流速稱為極限流速(vlim),對應的雷諾數稱為極限雷諾數(Relim)。將不同裂隙寬度條件下,極限流速以及極限雷諾數隨裂隙開度的變化作對比分析,分別如圖4(a)(b)所示?？傮w來講,極限流速或極限雷諾數隨裂隙寬度的增大而增大,隨裂隙開度的增大而減小,說明裂隙寬度越大,開度越小,立方定律越適用。本次研究中,極限流速的最大值為30.08 cm/s,極限雷諾數的最大值為75.2。

(a)極限流速 (b)極限雷諾數

2.2 修正立方定律的提出

(8)

式中的n值可由S1模擬結果擬合確定(NS方程模擬值與修正立方定律計算值對比結果見附件3)。不同裂隙寬度條件下,n值隨裂隙開度的變化如圖5所示。結果表明,裂隙開度越大,n值越大;總體來講,裂隙寬度越大,n值越小,本次研究中,n值分布于1.4～1.9。

圖5 不同裂隙寬度條件下,n值隨裂隙開度的變化

2.3 矩形狹縫平行裂隙滲透系數分析

根據修正立方定律和達西定律,矩形狹縫平行裂隙滲透系數可以用下式計算:

(9)

分析矩形狹縫裂隙滲透系數隨裂隙寬度、裂隙開度以及裂隙寬度與裂隙開度比值的變化(分別如圖6(a)(b)(c)所示),同一裂隙開度條件下,滲透系數隨裂隙寬度的增大而增大,但增大幅度相對較小。隨裂隙開度的增大而增大,增大幅度相對較大。滲透系數隨裂隙寬度與裂隙開度比值的增大而增大,通過曲線擬合得到二者之間呈指數關系。本次研究中,矩形狹縫裂隙滲透系數分布于10cm·s-1-70 cm·s-1。

圖6 滲透系數隨裂隙寬度、裂隙開度以及裂隙寬度與裂隙開度比值的變化圖

2.4 流速分布剖面分析

2.4.1 裂隙開度對流速分布剖面影響分析(S2模擬結果)

不同裂隙開度條件下,X=10 cm處流速分布如圖7所示(以b=0.05,0.1 cm為例,其余結果見附件4)。從圖中可以看出,二維裂隙面流速剖面形狀變化與理想泊肅葉流體流速變化形狀基本一致,數值上有所差異。流速的最大值位于裂隙中心線處。流速最小值位于裂隙壁處且為零。分析不同開度條件下,流速模擬值與泊肅葉值之差(絕對偏差)及相對偏差沿質點位置(Z軸)的變化(如圖8所示),絕對偏差隨裂隙開度的增大而增大,且離裂隙中心線(X軸)越近,偏差越大。相對偏差沿質點位置(Z軸)呈波動變化,變幅很小;相對偏差隨裂隙開度的增大而增大,表明裂隙開度越大,流速分布越不符合標準拋物線分布。

圖7 不同裂隙開度條件下,流速剖面分布

圖8 不同開度條件下,流速模擬值與泊肅葉值的絕對偏差及相對偏差沿質點位置(Z軸)的變化圖

2.4.2 流速大小對流速分布影響分析(S3模擬結果)

不同流速條件下,流速分布剖面如圖9所示。從圖中可以看出,無論流速大小,模擬所得流速分布大致符合拋物線分布,流速最大值位于裂隙中心線處(即X軸上),離裂隙壁越近,流速越小,在裂隙壁處流速為零。分析不同流速條件下,流速模擬值與泊肅葉值之差(絕對偏差)及相對偏差變化(如圖10所示),流速越大,絕對偏差越大,在裂隙中心線處最大,離裂隙壁越近,偏差越小。而相對偏差隨流速的變化幅度相對較小,表現為集中型。

圖9 不同流速條件下,流速剖面分布對比(NS值與泊肅葉值對比圖)

圖10 不同流速條件下,流速模擬值與泊肅葉值的絕對偏差及相對偏差沿質點位置(Z軸)的變化圖

2.4.3 裂隙寬度對流速分布影響分析(S4模擬結果)

取XY面中心線處Z=b/2,沿X=L/2質點流速為研究對象,不同裂隙寬度條件下,流速分布如圖11所示,從圖中可以看出流速值沿Y軸波動較小,在裂隙壁處驟減為0,隨著裂隙寬度的增大,流速波動幅度變小,裂隙寬度大于等于3 cm(對應的裂隙寬度與開度之比W/b≥40)的四條曲線重合度較高,表明隨著裂隙寬度及裂隙寬度與開度之比的增大,裂隙寬度的大小對流速分布曲線的影響減弱。

圖11 不同裂隙寬度條件下,XY面中心線處Z=b/2,沿X=L/2質點流速分布圖

3 結語

以矩形狹縫平行單裂隙為研究對象,對不同尺寸的矩形狹縫平行裂隙進行NS方程模擬計算,并將模擬結果與CL計算結果進行對比,首先對立方定律在矩形狹縫平行裂隙中的適用性進行驗證,發現流速較小時,CL值與NS值吻合較好,流速越大,立方定律越不適用,進而提出極限流速和極限雷諾數的概念來確定立方定律在矩形狹縫裂隙中的適用范圍,在本次研究中,極限流速的最大值為30.08 cm/s,極限雷諾數的最大值為75.2;基于修正立方定律和達西定律,提出矩形狹縫裂隙中立方定律的修正形式,式中的參數n值分布與1.4～1.9之間,矩形狹縫裂隙滲透系數隨裂隙寬度的變化幅度較小,隨裂隙開度的增大而增大,隨裂隙寬度與裂隙開度比值的增大而增大。分析流速剖面分布隨裂隙開度、裂隙寬度以及流速的變化規律,結果表明裂隙開度越大,流速越大,流速分布越不符合標準拋物線分布;裂隙寬度越大,對流速分布曲線的影響越弱。