高　瑜，葉　咸，夏　強

?

基于等效連續(xù)介質(zhì)模型的單裂隙滲流數(shù)值模擬研究

高瑜1,3，葉咸2,3，夏強3

(1.云南省地質(zhì)環(huán)境監(jiān)測院，云南 昆明 650216；2.云南省交通規(guī)劃設計研究院陸地交通氣象災害防治技術國家工程實驗室，云南 昆明 650041；3.成都理工大學地質(zhì)災害防治與地質(zhì)環(huán)境保護國家重點實驗室，四川 成都 610059)

基巖裂隙介質(zhì)的空間變異性導致裂隙滲流情況十分復雜。在等效連續(xù)介質(zhì)模型的基礎上，建立數(shù)值分析模型研究了裂隙寬度和滲透系數(shù)對單裂隙滲流的影響。利用單裂隙立方定律對裂隙流量進行解析計算，并在單一變量研究原則下，利用MODFLOW創(chuàng)建了光滑平行板單裂隙介質(zhì)的數(shù)值模型。據(jù)等效理論，將單裂隙寬度放大10倍的條件下，解析解計算結(jié)果比模擬結(jié)果高出6個數(shù)量級；不同裂隙寬度條件下，數(shù)值模擬結(jié)果與解析計算結(jié)果都表現(xiàn)出通過單個裂隙的水流量隨著裂隙的寬度的增加而增加；不同滲透系數(shù)條件下，穩(wěn)定流與非穩(wěn)定流條件下得到的流量都隨滲透系數(shù)的減小而變小。

等效連續(xù)介質(zhì)；單裂隙；裂隙寬度；滲透系數(shù)；數(shù)值模擬

巖體由于應力作用、風化作用等往往存在裂隙，這使得巖體的滲流性質(zhì)與一般的多孔介質(zhì)有很大差別，具有如各向異性、非均勻性等復雜性與特殊性。單裂隙面是巖體復雜裂隙網(wǎng)絡構(gòu)成的基本單元，其幾何特征如裂隙的方向、寬度、充填性和粗糙性等等，控制著巖體的滲透機制[1]。因此，只有首要研究單個裂隙面的滲透特征，才能為合理地預測工程巖體中復雜的滲流狀態(tài)提供可靠依據(jù)[2]。

目前，單裂隙滲流問題的研究成果多基于理論解析或者現(xiàn)場試驗，存在一定的局限性，因而將理論研究與數(shù)值模擬結(jié)果有利的結(jié)合，開展單個光滑裂隙的水流特征研究顯得十分必要。而國內(nèi)大部分數(shù)值模擬研究都是基于等效連續(xù)介質(zhì)開展的，如隧道涌水量的計算、巖溶管道水流量[3]等。模擬結(jié)果與實際情況之間存在著怎樣的誤差，誤差有多大，其中有哪些重要的影響因素，關于這些方面的研究則寥寥無幾。本文在等效連續(xù)介質(zhì)模型的基礎上，利用MODFLOW進行單裂隙滲流的數(shù)值分析，對于巖體復雜裂隙網(wǎng)絡的滲流特征研究具有顯著的意義。

1　單裂隙滲流概念模型及其解析解

巖體裂隙在分布和發(fā)育上表現(xiàn)為不均勻性和各向異性。與之相對應的，裂隙介質(zhì)的地下水滲流顯然也存在非均質(zhì)各向異性的特點。裂隙介質(zhì)存在氣態(tài)、固態(tài)、液態(tài)三相特征，氣態(tài)主要是空氣，并且大多存在于非飽和帶中；固態(tài)稱為固體骨架；液態(tài)是地下水、水與其他物質(zhì)的混合物或是其他流體(如石油等)。地下水有在裂隙巖體中的存在形式較多，一般有薄膜水、吸附水、毛管力和重力水等，本文主要研究的是重力水的運動。由于巖體中的裂隙系統(tǒng)十分復雜，由多條平行或者交叉的單裂隙組成，而裂隙巖體的透水性比巖體基質(zhì)部分的透水性大很多，故研究將巖石本身的滲透性忽略不計，可以視為流體僅在裂隙中流動[4]。

光滑平行板模型是單裂隙最基本的概念模型(如圖1)，假設裂隙由兩塊一定間距且表面光滑的平行板構(gòu)成。該概念模型將復雜的單裂隙滲流模型簡單化，利用單裂隙的立方定律可以很容易的計算和預測出兩平行板之間的穩(wěn)定滲流量。由于光滑平行板模型構(gòu)成簡單，在預測裂隙中的滲流時，平行板模型和立方定律是極易讓人接受的[5]。本文在光滑平行板模型的基礎上展開研究。

圖1　光滑平行板裂隙模型

然而，光滑平行板模型是一種較為理想的假設情況，實際上巖體裂隙并非平面或平行的空間，且?guī)r體裂隙面的粗糙度和隙寬具有空間變異性。因此立方定律在粗糙裂隙中的滲流的預測中是不合適的，使得我們對于實際工程問題中涉及到裂隙巖體的溶質(zhì)運移問題和水流運動問題難度增大。大多數(shù)情況下，要定量描述單裂隙的滲流問題需利用線性達西定律或基于達西定律的立方定律(簡稱LCL)。1856年，達西定律出現(xiàn)以后，被廣泛應用在空隙介質(zhì)和裂隙介質(zhì)中，達西認為裂隙介質(zhì)中的流體平均流速與水力梯度成線性關系，就服從達西定律。但是，許多研究成果都表明當雷諾數(shù)較大時，該線性關系就不再成立。假設裂隙內(nèi)速度場是穩(wěn)定、一維、一致的，線性的，構(gòu)成裂隙的板面是平行且光滑的，同時是忽略裂隙寬度所造成的邊界效應，則立方定律就演變成Navier-Stokes(NS)方程的簡單解?；谝陨霞僭O，立方定律得出了滲流量與裂隙開度的立方成比例的結(jié)論，相關實驗已經(jīng)驗證了該成果的可靠[6]。

Zimmerman 和Bodvarsson(1996年)[7]認為，當雷諾數(shù)小于1 時的流動，慣性力的影響將會變小，它導致滲流場非線性地依賴于水力梯度。粗糙裂隙中滲流的完全描述需要用到三維的NS方程(包含連續(xù)性方程)，NS方程僅被用來模擬具有理想幾何形狀的裂隙中的滲流[8]。

水在裂隙內(nèi)流動又分為層流和紊流，層流的運動規(guī)律和紊流是同樣的。當水流損失和與流速呈線性關系時，就稱為達西定律[9]。達西定律公式：

(1)

式中：V為水流的平均速度；J為水力坡度；K為滲透系數(shù)。

對于等寬度的光滑裂隙中的線性流運動，裂隙內(nèi)平均流速為：

(2)

式中：g為重力加速度；b為裂隙的寬度；vw為水的運動粘滯系數(shù)。并得裂隙的滲透公式：

(3)

則，光滑裂隙的單寬流量計算公式為：

(4)

光滑裂隙的總流量計算公式為：

Q=qL

天際線層次面積比例與評價結(jié)果間并不能視為單純的線性相關。同時，由于天際線樣本及受訪樣本數(shù)量的限制，當前無法確定具體的臨界點比例值。但為了從統(tǒng)計學的角度，對天際線層次面積比例與評價結(jié)果間的相關性進行驗證與說明，以下以0.442的比例值為大致臨界點，將“天際線層次面積比例”評價結(jié)果散點圖劃分為a、b兩段，分別視為兩段大致呈現(xiàn)出線性相關的數(shù)據(jù)集，便于后文相關系數(shù)的計算。

(5)

其中，水的運動粘滯系數(shù)：

(6)

式中：u為水的動力粘滯系數(shù)，ρ為水的密度。

2　MODFLOW數(shù)值模擬試驗

2.1數(shù)值模型的建立

通過查閱相關資料，比對多組數(shù)據(jù)，選取與模擬條件最接近的一組數(shù)據(jù)進行模型建立。模型為一個長20 m，寬20 m，高20 m的立方體。本次研究不考慮裂隙的傾斜角度，故將裂隙設置為水平裂隙，裂隙寬度b?；灸Ｐ椭衎=0.2 m。非穩(wěn)定流模擬時，模擬時間總長設定為10 a，即3 650 d，模型無降雨入滲補給。根據(jù)有限差分法的基本思路，在對模型離散化處理時，每個單元體在水平方向上的長度為0.5 m，因此模型被剖分成為40列，40行的網(wǎng)格。垂直方向上，將模型分成7層，1、2層和6、7層的滲透系數(shù)kx，ky相同，3、4、5層的滲透系數(shù)kx，ky相同，巖體和裂隙各向異性，kz取值均設置為kx或者ky的百分之一，即垂向滲透系數(shù)比水平方向的小。其中3、4、5層即為裂隙，1、2、6、7層為巖體。在基礎模型中，在穩(wěn)定流的基礎上，假設流經(jīng)裂隙的水流不變，將左右兩邊均設置為定水頭邊界。如圖2，將所有層位(第1、2、3、4、5、6、7層)的左側(cè)水頭高度設置為20 m，右側(cè)水頭高度設置為18 m。

圖2　邊界條件示意圖

2.2基本模型數(shù)值解與解析解

在穩(wěn)定流條件下運行模型，水均衡中in=37.059 m3/day，out=37.059 m3/day，平均值為37.059 m3/day，即為4.29×10-4m3/s。這個值即為通過裂隙的地下水流量。

本文根據(jù)單個裂隙的定義，建立一個用于基本研究的幾何模型，模型參數(shù)如表1所示，根據(jù)公式6，算得水的運動粘滯系數(shù)為Vm=1×10-6m2/s2。通過光滑裂隙的單寬流量為653.3 m3/s。則總流量Q為20×653.3 =13 066.7 m3/s。

表1　基本模型參數(shù)表

由上可知，基本模型的公式求解結(jié)果與模擬結(jié)果相差較大，考慮到解析算法和模型中的水流為兩種物理機制的水流，可以進行幾何近似調(diào)整，把實際的裂隙寬度放大10倍，再放到等效連續(xù)介質(zhì)模型中。所以，下文中解析算法的結(jié)果均為模型裂隙寬度的0.1倍計算所得。

2.3數(shù)值試驗方案的設定

本研究主要是通過設計一個單裂隙巖體的模型，對不同裂隙寬度和不同巖體和裂隙的滲透系數(shù)比值條件下的裂隙巖體中水流情況進行模擬，并初步探討裂隙寬度、滲透系數(shù)比值對單裂隙的滲流是否有影響，如果有影響，那么裂隙寬度和滲透系數(shù)比值分別是如何對單裂隙的滲流產(chǎn)生影響的。由立方定律可知，流經(jīng)平行光滑單裂隙的單寬流量只受到裂隙寬度的影響，與滲透系數(shù)無關。所以此處，為了使試驗具有可信度，在穩(wěn)定流條件下，分別對不同裂隙寬度和不同滲透系數(shù)比值情況下做了五組數(shù)值模擬，另加五組在非穩(wěn)定條件下不同滲透系數(shù)比值的數(shù)值模擬。非穩(wěn)定流條件下，需要增加參數(shù)貯水系數(shù)，裂隙的貯水率值參考巖溶地區(qū)的值，而除裂隙的巖體可以考慮用砂巖的。據(jù)《地下水流數(shù)值模擬技術要求》，設定巖體的貯水系數(shù)為3e-3，裂隙的貯水系數(shù)為5e-5[10]。以下所有模擬均建立在上一節(jié)單裂隙基本模型的基礎上。數(shù)值模擬試驗方案如表2所示。

表2　數(shù)值試驗方案表

3　數(shù)值模擬結(jié)果對比分析

3.1不同隙寬b模擬結(jié)果對比分析

通過數(shù)值模型模擬計算和解析計算，得到不同隙寬下的模擬結(jié)果和解析計算結(jié)果，見表3及圖3。如表3所示。隨著裂隙寬度的增加，流量也在增加。在模擬結(jié)果中，裂隙寬度與流量呈線性關系；在解析結(jié)果中，裂隙寬度與流量呈指數(shù)關系，如圖3。模擬結(jié)果與解析解計算結(jié)果對比分析顯示，計算結(jié)果與模擬的結(jié)果相差很大，解析解計算結(jié)果比模擬結(jié)果高出6個數(shù)量級。

表3　不同裂隙寬度流量對比表

圖3　不同裂隙寬度的流量計算結(jié)果曲線

3.2不同滲透系數(shù)比值模擬結(jié)果對比分析

3.2.1穩(wěn)定流模擬

從立方定律可知，流量的變化與滲透系數(shù)無關。但是從表6的模擬結(jié)果來看，隨著滲透系數(shù)的減小，流量變小。表4為隙寬保持在0.2 m時，改變滲透系數(shù)，其他參數(shù)設置按照數(shù)值試驗方案表(表2)模擬所得結(jié)果與解析計算結(jié)果對比表。模擬的裂隙流量隨著裂隙滲透系數(shù)的減小和巖體滲透系數(shù)和裂隙滲透系數(shù)比值的增大而減小，且呈拋物線形式如圖4所示，與改變隙寬的相比，流量曲線由直線變成拋物線。但是，計算的裂隙流量保持不變。

表4　不同滲透系數(shù)比值模擬流量表(穩(wěn)定流)

3.2.2非穩(wěn)定流

在利用Visual Modflow模擬非穩(wěn)定流問題時，模擬計算可以隨時停止、暫停，考慮到模型運行問題引起的誤差，現(xiàn)增加五組不同滲透系數(shù)下的非穩(wěn)定流數(shù)值模擬，模擬3 650 d的時間，將時間劃分為十個步長，輸出時間分別為140、309、511、754、1 046、1 396、1 816、2 319、2 924、3 650 days。

非穩(wěn)定流模擬條件下，模擬結(jié)果如表6所示，隨著滲透系數(shù)的減小，流經(jīng)單裂隙的流量變小。從數(shù)值上看，與穩(wěn)定流相比差別不大。模擬流量與滲透系數(shù)的關系如圖5所示，呈線性關系。

將模型第三層的速度參數(shù)導出，求平均值，得到在改變滲透系數(shù)比值的條件下的流速如表5。模擬流速隨著裂隙滲透系數(shù)的增大而增大，且呈線性關系，但是解析計算流速與滲透系數(shù)的大小無關，保持不變。

圖4　不同滲透系數(shù)流量模擬結(jié)果(穩(wěn)定流)

編號裂隙滲透系數(shù)kx,ky(m/s)模擬流速(m/s)計算流速(m/s)11.00E-036.67E-0421.00E-046.67E-0531.00E-056.67E-0641.00E-066.67E-0751.00E-076.67E-083.27E+01

表6　不同滲透系數(shù)比值模擬結(jié)果表(非穩(wěn)定流)

圖5　不同滲透系數(shù)流量模擬結(jié)果(非穩(wěn)定流)

將模型第三層最后時間步長的速度參數(shù)導出，求平均值，得到在改變滲透系數(shù)的條件下的流速。解析計算的流速不隨裂隙滲透系數(shù)的改變而改變，模擬流速隨著裂隙滲透系數(shù)的變小而變小，且變化情況與穩(wěn)定流的一致，即巖體滲透系數(shù)與裂隙滲透系數(shù)的比值越小，模擬所得裂隙中的水流速度越小。

4　結(jié)語

在單一變量研究原則的指導下，本文利用MODFLOW創(chuàng)建了光滑平行板單裂隙介質(zhì)的數(shù)值模型，分別對寬度和滲透系數(shù)做了五組模擬，共15種數(shù)值模型。通過模擬分析了不同寬度和滲透系數(shù)比值下單裂隙水流滲透的情況，得出以下結(jié)論：

(1)分析不同裂隙寬度條件下的流量，數(shù)值模擬結(jié)果中裂隙寬度與滲流量呈一次函數(shù)關系，解析計算結(jié)果中裂隙寬度與滲流量呈指數(shù)關系，均表現(xiàn)出隨著裂隙的寬度的增加，通過單個裂隙的滲流量增加。

(2)分析不同滲透系數(shù)比值條件下的流量，無論是穩(wěn)定流還是非穩(wěn)定流條件下得到的流量都隨滲透系數(shù)的減小而變小。穩(wěn)定流條件下，模擬結(jié)果中裂隙滲透系數(shù)與流量呈二次函數(shù)關系，非穩(wěn)定流條件下，模擬結(jié)果中裂隙滲透系數(shù)與流量呈一次函數(shù)關系。穩(wěn)定流和非穩(wěn)定流條件下的模擬流速都與裂隙系數(shù)呈線性關系，且隨裂隙滲透系數(shù)的增加而增加。

(3)以等效理論，將單裂隙寬度放大10倍的條件下，解析計算結(jié)果比模擬結(jié)果高出6個數(shù)量級。結(jié)果相差較大，說明利用等效連續(xù)介質(zhì)模型的模擬方法是存在問題的，現(xiàn)今對于裂隙水流的研究很多模擬都是基于等效連續(xù)介質(zhì)的模擬，這種方法的模擬結(jié)果值得懷疑，這將下一步研究工作的重點。

[1]艾萬民.巖體卸荷應力場-滲流場耦合試驗研究及數(shù)值分析[D].重慶:重慶大學.2005.

[2]宋曉晨,徐衛(wèi)亞.裂隙巖體滲流概念模型研究[J].巖土力學.2004,25(2):226-232.

[3]葉咸.昆明盆地北東側(cè)巖溶水系統(tǒng)探析[D].成都:成都理工大學.2014.

[4]黃勇.多尺度裂隙介質(zhì)中的水流和溶質(zhì)運移隨機模擬研究[D].江蘇:河海大學.2004.

[5]趙劍.公路隧道防排水機理及噴射防水試驗裝置的研究[D].重慶:重慶交通大學.2007.

[6]宋曉晨.裂隙巖體滲流非連續(xù)介質(zhì)數(shù)值模型研究及工程應用[D].江蘇:河海大學.2004.

[7]Zimmerman R W,Bodvarsson G S.Hydraulic conductivityof rock fractures[J].Transport in Porous Media,1996,23:1-30.

[8]Mourezenko V V,Thovert J -F,Adler P M.Permeabilityof a single fracture; validity of the Reynolds equation[J].Journal of Physics Ⅱ France,1995,(5):465-482.

[9]周志芳,王錦國.裂隙介質(zhì)水動力學[M].中國水利水電出版社.2004,1:1-302.

[10]徐軍祥,邢立亭.濟南泉域巖溶水數(shù)值預報與供水保泉對策[J].地質(zhì)調(diào)查與研究.2008,31(3):210-225.

2016-05-03

高瑜(1991-)，女，云南昆明人，助理工程師，主要從事水文地質(zhì)研究。

P641.73

A

1004-1184(2016)05-0040-04