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

(1.武漢大學(xué) 水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072； 2.武漢大學(xué) 水工巖石力學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072)

0 引 言

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

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

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

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

1 研究方法

1.1 粗糙裂隙生成

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

根據(jù)Brown模型，為考慮裂隙各向異性，將裂隙表面的功率譜密度G(k)表達(dá)為式(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所示，當(dāng)b/a=1時(shí)，裂隙開度場為各向同性；當(dāng)b/a<1時(shí)，裂隙開度場為各向異性，空間關(guān)聯(lián)性或連通性在y方向上優(yōu)于x方向，即y方向?yàn)殚_度場各向異性優(yōu)勢方向；當(dāng)b/a>1時(shí)，x方向?yàn)閮?yōu)勢方向。本文選取y方向?yàn)閮?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用于滲流計(jì)算。生成的隨機(jī)粗糙裂隙的編號和幾何參數(shù)如表1所列。圖1為F1H3A1、F1H3A2、F1H3A3的開度分布圖。

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

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

1.2 控制方程

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

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

(2)

?·u=0

(3)

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

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

圖2 流動(dòng)邊界條件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為裂隙垂直于流動(dòng)方向上的寬度，m，在本文中為0.12 m；eh為等效水力開度，m；β為非達(dá)西系數(shù)或慣性系數(shù)，1/m，取決于裂隙幾何結(jié)構(gòu)；bD為無量綱化的系數(shù)。

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

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

(7)

(8)

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

(9)

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

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

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

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

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

(10)

式中：n為采樣點(diǎn)數(shù)，zi為第i個(gè)點(diǎn)的z坐標(biāo)，Δd為采樣點(diǎn)間隔。

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

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

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

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

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

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

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

(11)

(12)

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

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

(13)

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

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

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

2.1 壓力梯度-流量關(guān)系

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

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

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

2.2 表觀滲透性

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

(14)

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

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

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

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

(15)

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

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

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

2.3 臨界雷諾數(shù)

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

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

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

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

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

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

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

(16)

(17)

圖9為各向異性模型和雙參數(shù)模型對粗糙裂隙F3H3A3的擬合結(jié)果。對于平行優(yōu)勢方向(y方向)的滲流，雙參數(shù)模型和各向異性模型均擬合較好。對于垂直優(yōu)勢方向(x方向)的滲流，各向異性模型擬合較好，雙參數(shù)模型有較大偏差。同時(shí)，可以觀察到不同方向的雙參數(shù)模型擬合結(jié)果較為接近，這是因?yàn)椴煌较虻姆逯荡植诙认嗤业刃﹂_度相差不大。當(dāng)裂隙形貌和流動(dòng)方向確定時(shí)，各向異性模型能更好地?cái)M合裂隙巖體中非線性滲流特征。

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

3.2 試驗(yàn)結(jié)果對比

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

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

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

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

(18)

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

3.4 表征模型影響因素討論

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

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

4 結(jié) 論

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

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

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