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

高　瑜，葉　咸，夏　強

?

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

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

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

基巖裂隙介質(zhì)的空間變異性導(dǎo)致裂隙滲流情況十分復(fù)雜。在等效連續(xù)介質(zhì)模型的基礎(chǔ)上，建立數(shù)值分析模型研究了裂隙寬度和滲透系數(shù)對單裂隙滲流的影響。利用單裂隙立方定律對裂隙流量進(jìn)行解析計算，并在單一變量研究原則下，利用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ù)值模擬

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

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

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

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

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

圖1　光滑平行板裂隙模型

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

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

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

(1)

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

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

(2)

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

(3)

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

(4)

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

Q=qL

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

(5)

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

(6)

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

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

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

通過查閱相關(guān)資料，比對多組數(shù)據(jù)，選取與模擬條件最接近的一組數(shù)據(jù)進(jìn)行模型建立。模型為一個長20 m，寬20 m，高20 m的立方體。本次研究不考慮裂隙的傾斜角度，故將裂隙設(shè)置為水平裂隙，裂隙寬度b?；灸Ｐ椭衎=0.2 m。非穩(wěn)定流模擬時，模擬時間總長設(shè)定為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取值均設(shè)置為kx或者ky的百分之一，即垂向滲透系數(shù)比水平方向的小。其中3、4、5層即為裂隙，1、2、6、7層為巖體。在基礎(chǔ)模型中，在穩(wěn)定流的基礎(chǔ)上，假設(shè)流經(jīng)裂隙的水流不變，將左右兩邊均設(shè)置為定水頭邊界。如圖2，將所有層位(第1、2、3、4、5、6、7層)的左側(cè)水頭高度設(shè)置為20 m，右側(cè)水頭高度設(shè)置為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é)果相差較大，考慮到解析算法和模型中的水流為兩種物理機(jī)制的水流，可以進(jìn)行幾何近似調(diào)整，把實際的裂隙寬度放大10倍，再放到等效連續(xù)介質(zhì)模型中。所以，下文中解析算法的結(jié)果均為模型裂隙寬度的0.1倍計算所得。

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

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

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

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

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

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

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

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

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

3.2.1穩(wěn)定流模擬

從立方定律可知，流量的變化與滲透系數(shù)無關(guān)。但是從表6的模擬結(jié)果來看，隨著滲透系數(shù)的減小，流量變小。表4為隙寬保持在0.2 m時，改變滲透系數(shù)，其他參數(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ù)的關(guān)系如圖5所示，呈線性關(guān)系。

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

圖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ù)導(dǎo)出，求平均值，得到在改變滲透系數(shù)的條件下的流速。解析計算的流速不隨裂隙滲透系數(shù)的改變而改變，模擬流速隨著裂隙滲透系數(shù)的變小而變小，且變化情況與穩(wěn)定流的一致，即巖體滲透系數(shù)與裂隙滲透系數(shù)的比值越小，模擬所得裂隙中的水流速度越小。

4　結(jié)語

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

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

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

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

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

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

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

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

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

[6]宋曉晨.裂隙巖體滲流非連續(xù)介質(zhì)數(shù)值模型研究及工程應(yīng)用[D].江蘇:河海大學(xué).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ì)水動力學(xué)[M].中國水利水電出版社.2004,1:1-302.

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

2016-05-03

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

P641.73

A

1004-1184(2016)05-0040-04