粗糙裂隙各向異性對非達西滲流特性的影響

馬 亞 楠，楊 志 兵，熊 小 鋒，陳 益 峰

(1.武漢大學 水資源與水電工程科學國家重點實驗室，湖北 武漢 430072； 2.武漢大學 水工巖石力學教育部重點實驗室，湖北 武漢 430072)

0 引 言

在高壩樞紐、高水頭抽水蓄能電站、煤礦開采和深埋隧洞等工程活動中，受滲透水壓和裂隙形貌的影響，巖體滲流具有強烈的非線性，發(fā)生滲漏、突透水和滲透破壞的風險急劇增大[1-4]，因此非線性滲流對于工程建設和運行穩(wěn)定有著重要影響。巖體滲流主要是通過裂隙結構面進行的。受地質沉積環(huán)境、溫度變化和地質構造作用等因素的影響，天然裂隙面存在各向異性特征[5]，使得裂隙滲流模擬及預測存在不確定性。研究各向異性粗糙裂隙中的非線性滲流規(guī)律，對于預測宏觀尺度裂隙巖體滲流，預防隧洞涌水、內水外滲、煤礦突水和瓦斯突出等事故具有重要理論意義和參考價值。

近年來已有許多學者開展了關于裂隙巖體滲流的研究[6-7]。基于平行板假設的立方定律被廣泛用于描述裂隙介質中的滲流[8]。以往的研究針對天然裂隙粗糙特性的影響，提出了一系列修正的立方定律。偏離線性規(guī)律的非達西滲流現(xiàn)象在數(shù)值模擬、室內外試驗和工程應用中廣泛存在[9-10]。研究發(fā)現(xiàn)裂隙滲流的非線性受到法向應力[11]、剪切錯動[12]、接觸面積[13]等因素的影響，其本質上是裂隙幾何形狀和開度的變化對裂隙滲流的控制作用。非達西滲流細觀機制和宏觀表征的研究得到了普遍關注。在細觀機制方面，理論分析結果顯示當慣性效應占主導作用時滲流出現(xiàn)非達西特征[14]；此外，數(shù)值模擬結果表面渦旋區(qū)的產(chǎn)生和發(fā)展是非達西滲流發(fā)生的主要原因[15]。在宏觀表征方面，Izbash公式和Forchheimer公式常用于描述非達西滲流的壓力梯度與流量關系。其中Forchheimer公式物理意義明確，參數(shù)取值穩(wěn)定性強，被廣泛用于裂隙介質非線性流的表征中。目前主要通過建立裂隙形貌與非達西滲流公式參數(shù)的關系來對非線性系數(shù)進行量化，如采用峰值起伏度[16]、水文彎曲度和表面分形冪律關系[17]或提出新的參數(shù)[18]來表征裂隙形貌，建立起非線性系數(shù)表征模型。但這些研究沒有考慮裂隙的各向異性及滲流方向性的影響。

開度場的各向異性可顯著影響裂隙滲流規(guī)律。剪切作用是比較常見的導致粗糙裂隙產(chǎn)生各向異性特征的方式。在各向異性粗糙裂隙中，裂隙滲透特性受各向異性程度、相對粗糙度等因素的共同影響[19-20]；在剪切粗糙裂隙中，垂直于剪切方向的滲透性往往大于平行于剪切方向的滲透性[21]。剪切裂隙各向異性程度與剪切位移、平均開度、粗糙度、法向應力等因素密切相關[22-24]；同時，剪切過程中也伴隨著開度場的改變和通道流的形成與變化。目前對于各向異性粗糙裂隙滲流規(guī)律的研究成果多為定性認識，且對滲流的各向異性與開度、粗糙度等裂隙幾何形貌特征的關系研究相對較少，缺乏對各向異性裂隙非線性滲流特性的量化表征，因此亟需深入研究。

本文采用地質統(tǒng)計隨機模型生成表面符合分形特征的不同各向異性程度的粗糙裂隙，通過直接數(shù)值求解Navier-Stokes方程進行模擬，分析粗糙裂隙的非線性滲流特征，闡明裂隙的各向異性程度、裂隙開度、粗糙度等因素對非線性滲流的影響規(guī)律，提出考慮裂隙各向異性和粗糙度的非線性滲流參數(shù)表征模型，并推導出相應的臨界雷諾數(shù)表達式，為裂隙網(wǎng)絡及宏觀尺度巖體滲流研究提供理論依據(jù)。

1 研究方法

1.1 粗糙裂隙生成

本文采用Brown[25]提出的數(shù)學模型來生成隨機粗糙裂隙。在該模型中，通過分形維數(shù)D和標準差控制裂隙上、下表面的粗糙度，通過設置不匹配長度來控制上下面的不匹配程度及開度的空間分布。本文選取的分形維數(shù)為2.5，不匹配波數(shù)為8(不匹配長度=裂隙長度/不匹配波數(shù))。開度均值和標準差的選取參照文獻[17]試驗中的數(shù)值。選擇平均開度為0.5，1.0，1.5 mm；開度標準差為0.1，0.2，0.3。采用Brown方法生成的裂隙為隨機裂隙，其中開度為0.5 mm和1.0 mm的裂隙是通過開度為1.5 mm的裂隙平移上表面得到的。

根據(jù)Brown模型，為考慮裂隙各向異性，將裂隙表面的功率譜密度G(k)表達為式(1)[26]：

G(k)∝[(kx/a)2+(ky/b)2]-(3.5-D)

(1)

式中：kx，ky為x，y方向上的波數(shù)，D為裂隙表面的分形維數(shù)，參數(shù)a和b表征裂隙的各向異性程度。

如圖1所示，當b/a=1時，裂隙開度場為各向同性；當b/a<1時，裂隙開度場為各向異性，空間關聯(lián)性或連通性在y方向上優(yōu)于x方向，即y方向為開度場各向異性優(yōu)勢方向；當b/a>1時，x方向為優(yōu)勢方向。本文選取y方向為優(yōu)勢方向，各向異性程度表征參數(shù)b/a取值1，1/2，1/4用于研究裂隙表面的各向異性對滲流的影響。

本文中裂隙面的分辨率為512×512，生成裂隙面的大小為128 mm×128 mm，x-y方向網(wǎng)格大小為0.25 mm，截取其中120 mm×120 mm用于滲流計算。生成的隨機粗糙裂隙的編號和幾何參數(shù)如表1所列。圖1為F1H3A1、F1H3A2、F1H3A3的開度分布圖。

圖1 不同各向異性程度b/a的開度分布Fig.1 Aperture field of different anisotropy factors b/a

表1 生成裂隙幾何參數(shù)

1.2 控制方程

對于單個粗糙裂隙中不可壓縮牛頓流體的定常流動，根據(jù)動量和質量守恒定律，流動過程可用Navier-Stokes方程和連續(xù)性方程來描述[27]。

ρ(u·?u)-μ?2u=-?P

(2)

?·u=0

(3)

式中：u為速度矢量，m/s；P為壓力，Pa；μ為動力黏度，Pa·s；ρ為流體密度，kg/m3；等式(2)中，ρ(u·?u)和μ?2u分別代表流動過程中的慣性力和黏性力，式(3)是代表質量守恒的連續(xù)性方程。

本次模擬中入口處設定為流量邊界條件，出口處設定為0壓力邊界條件，其余設置為無滑移邊界。z方向劃分網(wǎng)格7～10層。模型網(wǎng)格數(shù)量(2～4)×106。通過有限元軟件求解Navier-Stokes方程[15]得到不同裂隙的流場。分別將流動方向設置為沿x方向和y方向，即垂直于優(yōu)勢方向和平行于優(yōu)勢方向，來研究裂隙滲流的各向異性。圖2為沿x方向和沿y方向的兩種流動情形的邊界條件。

圖2 流動邊界條件Fig.2 Flow boundary conditions

Forchheimer方程常被用來宏觀表征非線性流，具體如下：

-?P=AQ+BQ2

(4)

(5)

(6)

式中：Q為流量，m3/s；A為描述黏性力能量損耗的系數(shù)，kg/(m5·s)；B為描述慣性力能量損耗的系數(shù)，kg/m8；w為裂隙垂直于流動方向上的寬度，m，在本文中為0.12 m；eh為等效水力開度，m；β為非達西系數(shù)或慣性系數(shù)，1/m，取決于裂隙幾何結構；bD為無量綱化的系數(shù)。

本文中設eh為僅由幾何形狀控制的變量，沿用Fourar等[28]和Nowamooz等[29]的定義，通過線性階段的系數(shù)A反算得到。

在裂隙滲流中，常用雷諾數(shù)Re和Forchheimer數(shù)Fo來表征流態(tài)的變化，定義為如下形式：

(7)

(8)

Re可以反映慣性力和黏性力的相對強弱，F(xiàn)o可以反映流動偏離線性的程度。目前常使用非線性程度因子α來區(qū)分達西流和非達西流。

(9)

當α大于臨界值時，可以認為水流已經(jīng)發(fā)展為非達西流態(tài)，參照Zeng & Grigg[30]等的取值，臨界值取0.1。

模擬的雷諾數(shù)范圍為0.01～300。滲流數(shù)值模擬參數(shù)取值如表2所列。

表2 滲流數(shù)值模擬模型參數(shù)

1.3 裂隙粗糙度及各向異性的表征

目前裂隙粗糙程度表征方法可歸納為三大類：經(jīng)驗方法(JRC)、統(tǒng)計參數(shù)(平均值、均方根、方差等)和分形維數(shù)表征法[31]。本文采用一階導數(shù)均方根Z2來描述裂隙的粗糙度[32]。

(10)

式中：n為采樣點數(shù)，zi為第i個點的z坐標，Δd為采樣點間隔。

對于三維曲面，可以在流動方向上取多條線計算Z2值，取其平均值作為該曲面的Z2值。根據(jù)生成裂隙表面的精度，沿x方向取512條剖面線，計算平均值作為垂直優(yōu)勢方向的Z2值，Z2，⊥。沿y方向取剖面線平均值為平行優(yōu)勢方向的Z2值，Z2，∥。

如表3所列，對于隨機生成的裂隙，在b/a不變的情況下，其Z2值隨均方差值增大而增大，但Z2，⊥/Z2，∥的值維持穩(wěn)定，說明Z2可以較好地表示裂隙的粗糙度及各向異性。

隨機生成多組不同各向異性程度的粗糙裂隙并計算其Z2，⊥/Z2，∥值，發(fā)現(xiàn)Z2，⊥/Z2，∥與裂隙的b/a值有很強的相關關系，擬合相關度大于0.99。Yin等的研究表明，Z2與分形維數(shù)和標準差之間也存在很好的相關性[33]。因此，可以用Z2來表征裂隙粗糙度及各向異性。

表3 不同開度均方差隨機裂隙的一階導數(shù)均方根

圖3 Z2，⊥/Z2，∥與各向異性表征指標b/a的關系Fig.3 Relationship between Z2，⊥/Z2，∥ and anisotropy factor b/a

1.4 各向異性粗糙裂隙非線性系數(shù)表征模型

Chen等提出了雙參數(shù)模型[34]，用峰值粗糙度和等效水力開度對式(6)中的無量綱非達西系數(shù)進行表征。根據(jù)不同粗糙度裂隙在法向應力下的流量-壓力梯度實驗數(shù)據(jù)，借助Levenberg-Marquardt優(yōu)化算法，得到模型擬合系數(shù)的值：m=0.022，n=0.666。

(11)

(12)

式中：ξ為裂隙的峰值粗糙度，m。

對于各向異性裂隙，不同滲流方向的ξ值相同，因此上述表征模型不能很好地體現(xiàn)裂隙各向異性對非線性滲流的影響。本文參考上述雙參數(shù)模型，用Z2和eh作為參數(shù)來表征非達西系數(shù)β，具體形式如下：

(13)

式中：κ，ψ，γ為擬合系數(shù)。

本文選取的表達式為經(jīng)驗表達式，不是嚴格的理論公式，選取該形式基于以下原因：① 形式上與雙參數(shù)公式相似，用Z2代替ξ來表征裂隙形貌對于非線性滲流的影響。② 非達西系數(shù)與eh近似呈冪函數(shù)關系，在多個研究中已經(jīng)得到證明[35]。③ 對于光滑的平板模型情況，有Z2=0，從而得到β=0。公式可以退化為光滑平板模型[8]。

2 數(shù)值模擬結果和分析

2.1 壓力梯度-流量關系

根據(jù)數(shù)值模擬得到的結果，用Forchheimer公式對裂隙滲流的壓力梯度和流量曲線進行擬合。如圖4所示，不同粗糙度、機械開度和各向異性程度下裂隙滲流的?P～Q曲線均與Forchheimer公式擬合良好(R2>0.99)，曲線呈現(xiàn)非線性特征，說明在選擇的進口流量范圍內，裂隙均已發(fā)展為非線性流動。

圖4 壓力梯度-流量關系Fig.4 Relationship between pressure gradient and flow rate for different fractures

從圖4(a)、(b)可以看出，隨著粗糙度增加和開度的減小(即相對粗糙度的增加)，?P～Q曲線的非線性程度增強。從圖4(c)可以看出，對于b/a為1的裂隙，垂直于優(yōu)勢方向(x方向)與平行于優(yōu)勢方向(y方向)的?P～Q曲線基本一致。隨b/a的減小，開度場各向異性增強，不同方向的?P～Q曲線差別變大，且可以看出，垂直于優(yōu)勢方向的非線性程度更強。

2.2 表觀滲透性

采用表觀水力開度來反映裂隙滲流的表觀滲透性，以此研究滲流的非達西特性。

(14)

式中：Ta為表觀導水系數(shù)，m4；ea為表觀水力開度，m。

圖5給出了不同各向異性程度的裂隙水力開度ea和機械開度em之比隨雷諾數(shù)的變化。可以看出，在雷諾數(shù)很小時(Re<1)，ea/em幾乎保持不變，此時滲流處于達西狀態(tài)，滲透率接近固有滲透率。隨著雷諾數(shù)的不斷增大，ea/em不斷減小，裂隙滲透性不斷降低。從圖5可以看出，對于b/a為1的粗糙裂隙，垂直于優(yōu)勢方向與平行于優(yōu)勢方向的ea/em值基本一致，說明兩個方向的滲透性大致相等。b/a<1時，不同方向的ea/em值差異較大，且平行于優(yōu)勢方向滲透性更強。

圖5 水力開度與機械開度比值隨雷諾數(shù)變化Fig.5 Relationship between ratio of hydraulic aperture to mean aperture and Reynolds number for fractures of different anisotropy factors

為了量化導水系數(shù)的各向異性，定義η為導水系數(shù)比：

(15)

式中：Ta，x，Ta，y為垂直于優(yōu)勢方向與平行于優(yōu)勢方向的表觀導水系數(shù)。

圖6為不同粗糙度、不同開度、不同各向異性程度的裂隙導水系數(shù)比隨雷諾數(shù)的變化。從圖中可以看出，在雷諾數(shù)較小時，η值基本不變。隨著雷諾數(shù)增大，η值不斷增大，隨著水流慣性效應的增大，流動狀態(tài)更加復雜，平行于流動方向的導水系數(shù)變化較小，因此η值不斷增大。從圖6(a)和(b)可以看出，隨相對粗糙度的增大，裂隙內空間結構更加復雜，導水系數(shù)差異更大。從圖6(c)可以看出，當b/a為1時，η值接近于1，說明垂直于優(yōu)勢方向與平行于優(yōu)勢方向的導水性在線性階段和非線性階段都基本相等，裂隙滲流具有各向同性。當b/a從1/2減小到1/4時，導水性各向異性程度增大。

圖6 導水系數(shù)比η隨雷諾數(shù)變化Fig.6 Relationship between ratio of hydraulic conductivity and Reynolds number for different fractures

2.3 臨界雷諾數(shù)

非線性程度因子取0.1時可以得到流動的臨界雷諾數(shù)。圖7給出了臨界雷諾數(shù)隨著各向異性程度的變化。從圖7中可以看出，在b/a為1的裂隙中，垂直于優(yōu)勢方向與平行于優(yōu)勢方向的臨界雷諾數(shù)基本一致。隨著b/a的減小，不同方向的臨界雷諾數(shù)差異變大，且平行于優(yōu)勢方向的臨界雷諾數(shù)大于垂直優(yōu)勢方向的臨界雷諾數(shù)，即更不容易發(fā)生非線性流動。

圖7 非達西程度因子α隨雷諾數(shù)變化Fig.7 Relationship between non-Darcian flow factor and Reynolds number for fractures of different anisotropy factors

3 各向異性非線性系數(shù)表征模型驗證

3.1 擬合系數(shù)的確定

根據(jù)數(shù)值模擬得到不同裂隙的?P～Q數(shù)據(jù)，采用Forchheimer公式擬合，得到相應的系數(shù)A、B，根據(jù)式(5)、(6)反算得到eh和β。借助Levenberg-Marquardt優(yōu)化算法，根據(jù)式(13)對模擬值進行回歸分析，得到擬合系數(shù)κ=0.006 13，ψ=1.443 58，γ=1.558 9，相關系數(shù)R2=0.968，表明該模型可以很好表征各向異性粗糙裂隙滲流非線性系數(shù)。

注：圖中數(shù)據(jù)點代表數(shù)值模擬值，曲面為模型擬合結果。圖8 參數(shù)模型回歸分析結果Fig.8 Regression analysis to obtain model parameters

將擬合結果代入式(13)、(6)中，可以得到各向異性粗糙裂隙非線性滲流的參數(shù)化模型：

(16)

(17)

圖9為各向異性模型和雙參數(shù)模型對粗糙裂隙F3H3A3的擬合結果。對于平行優(yōu)勢方向(y方向)的滲流，雙參數(shù)模型和各向異性模型均擬合較好。對于垂直優(yōu)勢方向(x方向)的滲流，各向異性模型擬合較好，雙參數(shù)模型有較大偏差。同時，可以觀察到不同方向的雙參數(shù)模型擬合結果較為接近，這是因為不同方向的峰值粗糙度相同，且等效水力開度相差不大。當裂隙形貌和流動方向確定時，各向異性模型能更好地擬合裂隙巖體中非線性滲流特征。

圖9 F3H3A3裂隙不同方向上壓力梯度-流量關系Fig.9 Relationship between pressure gradient and flow rate of different directions in Fracture F3H3A3

3.2 試驗結果對比

Rong等通過巴西劈裂得到不同粗糙度裂隙，并開展了滲流試驗[18]。圖10為試驗數(shù)據(jù)與本文模型擬合結果，可以看出，該模型可以較好地預測裂隙的非線性系數(shù)，表明本文提出的各向異性裂隙非線性滲流參數(shù)表征模型對于試驗結果仍然有很好的適用性。

注：圖中數(shù)據(jù)點為試驗測量值，曲面為模型擬合結果。圖10 試驗數(shù)據(jù)與模型擬合結果Fig.10 Experimental results and model fitting results

3.3 臨界雷諾數(shù)模型

根據(jù)非線性程度因子α的表達式(9)和非線性項系數(shù)B的參數(shù)化表征公式(17)，可以推導出適用于各向異性粗糙裂隙滲流的臨界雷諾數(shù)計算公式

(18)

計算得到的臨界雷諾數(shù)在15～360之間，數(shù)值模擬得到的臨界雷諾數(shù)10～200之間，兩者范圍基本一致，說明使用式(18)計算臨界雷諾數(shù)結果是可靠的。

3.4 表征模型影響因素討論

本文提出的各向異性非線性系數(shù)表征模型適用于不同粗糙度、開度和各向異性程度的裂隙。當裂隙表面分形維數(shù)變化時，因為分形維數(shù)與Z2相關性較強[33]，理論上表征模型的形式可保持不變，但擬合系數(shù)可能會有所變化，因此還需要進一步研究。

此外，本文采用考慮各向異性的隨機數(shù)學模型生成粗糙裂隙，因而開度場的關聯(lián)長度具有各向異性，且只需要由一個不匹配波數(shù)(空間頻率)來控制。當不匹配波數(shù)減小時，各個方向關聯(lián)長度均相應地增大，但放大裂隙的尺寸以后(因此需要在更大尺度上進行滲流模擬計算)，開度分布仍與原來相似，該模型將仍適用。

4 結 論

本文基于地質統(tǒng)計理論和隨機模型生成了不同粗糙度、開度和各向異性程度的粗糙裂隙，開展了三維數(shù)值模擬研究，重點探究了粗糙裂隙各向異性對非線性滲流特征的影響，提出并驗證了考慮各向異性的裂隙非線性滲流參數(shù)表征模型，主要結論如下：

(1) 與沿平行于優(yōu)勢方向的流動相比，垂直于優(yōu)勢方向的滲透性更弱，且更易發(fā)展為非線性流動。用導水系數(shù)比η表征表觀導水度的各向異性，其隨相對粗糙度和開度場各向異性的增大而增大。在線性階段η值基本恒定，隨著雷諾數(shù)的增加，流動發(fā)展為非線性，η逐漸增大。

(2) 本文提出以Z2和eh為參數(shù)的粗糙裂隙非線性系數(shù)表征模型，可用于各向異性粗糙裂隙，因而與文獻中已有的類似模型相比更具有通用性。通過對比發(fā)現(xiàn)模型與試驗結果有較好的一致性。根據(jù)提出的模型，推導出臨界雷諾數(shù)的經(jīng)驗計算公式，公式計算結果與模擬結果大致相符。