王俊奇,李鵬宇,董 曄,穆孟婧

(華北電力大學(xué) 可再生能源學(xué)院,北京 102206)

10.3969/j.issn.1671-8798.2017.05.007

2017-04-19

水資源與水電工程科學(xué)國家重點實驗室開放基金項目(2016SGG03)

王俊奇(1971— ),男,山西省朔州人,副教授,博士,主要從事水工結(jié)構(gòu)工程研究。E-mail:jqwangmail@163.com。

巖體裂隙網(wǎng)絡(luò)管單元多參數(shù)模型滲透張量研究

王俊奇,李鵬宇,董 曄,穆孟婧

(華北電力大學(xué) 可再生能源學(xué)院,北京 102206)

根據(jù)現(xiàn)場實測的巖體裂隙產(chǎn)生狀況,采用圓盤裂隙假設(shè),用Monte-Carlo方法生成3組裂隙形成裂隙網(wǎng)絡(luò)?；谒诹严吨械臏喜哿餍螒B(tài)將三維裂隙網(wǎng)絡(luò)簡化成空間一維管單元模型,且認為管單元的等效管徑與隙寬呈線性關(guān)系。設(shè)置每組裂隙等效管徑與隙寬比值為待求參數(shù),利用遺傳算法反演得到,進而求取等效滲透張量,并與整體單參數(shù)模型獲得的等效滲透張量相比較。與工程實測滲透張量對比,結(jié)果表明多參數(shù)模型較單參數(shù)模型能夠更好地表征滲透張量。

多參數(shù)模型;三維裂隙網(wǎng)絡(luò);管單元;滲透張量

滲流是流體通過多孔介質(zhì)或裂隙介質(zhì)的流動,裂隙巖體中流體的滲流分析對人類生產(chǎn)生活非常重要,如石油儲存、核廢料地下貯存等方面[1],因此裂隙巖體的滲流問題一直受到各國學(xué)者的廣泛關(guān)注。研究發(fā)現(xiàn)巖體中結(jié)構(gòu)面具有不連續(xù)性且發(fā)育極不規(guī)則,裂隙網(wǎng)絡(luò)的幾何構(gòu)成非常復(fù)雜。為方便研究裂隙的滲透性,通過計算機技術(shù)生成離散裂隙網(wǎng)絡(luò)模型,對仿真模型滲透特性展開研究的較多[2-4]。在這些研究中,利用滲透張量來表征裂隙滲透特性的方法得到廣泛的應(yīng)用,通過張量的性質(zhì)可以將非連續(xù)裂隙巖體的滲流問題轉(zhuǎn)化為連續(xù)介質(zhì)模型求解滲流場的問題[5],從而簡化了計算,使得大規(guī)模工程問題的解決成為可能。目前,研究滲透張量主要有實驗法和數(shù)值法兩種方法。在實驗法確定滲透張量方面,周志芳等[6]基于振蕩式微水試驗原理,通過單孔微水振蕩試驗現(xiàn)場確定了巖體的滲透系數(shù)張量。楊金保等[7]對取自丹江口水庫庫區(qū)的輝綠巖剪切裂隙進行了不同圍壓和裂隙水壓下的滲透試驗,基于該試驗對滲透系數(shù)進行考慮三向應(yīng)力和裂隙水壓力變化的修正,得到新的滲透張量的表達式。李克克等[8]通過現(xiàn)場實測獲得裂隙結(jié)構(gòu)面參數(shù),采用離散元法進行滲流分析。在數(shù)值法確定滲透張量方面,王培濤等[9]應(yīng)用離散裂隙網(wǎng)絡(luò)模型,通過平面滲流分析方法研究了節(jié)理巖體不同幾何分布情況下的滲透張量特征。Baghbanan等[10]基于UDEC計算平臺,采用離散裂隙網(wǎng)絡(luò)模型對滲透張量進行了研究,并對滲透的各異性進行了分析。從文獻資料中可以看出,數(shù)值法主要是在二維模型的情況下對裂隙滲流進行分析計算,這與實際情況相差較遠。在三維模型的研究方面,李新強等[11]使用邊界元法來計算三維裂隙網(wǎng)絡(luò)的滲透張量,通過滲透張量來定量分析裂隙滲流。Lee等[12]根據(jù)裂隙的幾何分布性質(zhì)生成裂隙網(wǎng)絡(luò),在裂隙圓盤上劃分三角形單元網(wǎng)絡(luò),用數(shù)值方法來模擬裂隙中的復(fù)雜流動和二氧化碳的運移問題。Cacas等[13]認為水在裂隙內(nèi)成溝槽流形態(tài),提出將三維裂隙網(wǎng)絡(luò)簡化成一維管單元,為滲流的分析計算提供了新的思路。Wang[14]等利用壓水實驗對管單元模型結(jié)果進行校正,通過回歸分析得到滲透張量。于青春等[15]建立了非連續(xù)裂隙網(wǎng)絡(luò)管單元滲流模型并對模型進行了校正。王俊奇等[16]采用單參數(shù)模型,通過反演法確定單元管徑,利用有限單元法計算得到裂隙巖體的滲透張量。

由于管單元模型接近水流在裂隙間的流動形態(tài),而且能有效地減少計算量,因此本研究也采用管單元模型進行滲流分析;不同之處在于,以往的管單元模型中只有一個參數(shù)作為等效管徑和隙寬的比值,本研究將在三維裂隙網(wǎng)絡(luò)模型中生成多組裂隙網(wǎng)絡(luò),并且對各組裂隙進行標記,設(shè)置相應(yīng)的參數(shù)作為管徑和隙寬的比值,通過有限單元法研究多參數(shù)模型確定的滲透張量。

1 三維裂隙網(wǎng)絡(luò)模型的建立與滲流計算方法

在工程實踐中,巖體結(jié)構(gòu)面的分布表面上看起來是雜亂無序的,但是通過將其分布的特征數(shù)據(jù)進行匯總可以發(fā)現(xiàn),結(jié)構(gòu)面的分布符合一定的概率特征。通過野外的實地測量采樣,進行統(tǒng)計分析,建立巖體幾何參數(shù)的統(tǒng)計模型,利用計算機技術(shù)生成統(tǒng)計特性與實地測量相符的裂隙網(wǎng)絡(luò)模型[17]。

1.1 裂隙網(wǎng)絡(luò)模型的生成

經(jīng)過多年的發(fā)展,目前有多種巖體結(jié)構(gòu)面的幾何系統(tǒng)模型用來表征巖體結(jié)構(gòu)系統(tǒng)特性的特定組合。在眾多的模型中,Baecher模型因其簡便實用而被廣泛應(yīng)用于巖石力學(xué)的相關(guān)分析中[18]。本研究采用這種模型假設(shè)巖體中的裂隙形狀為圓形,通過實地調(diào)查確定裂隙幾何參數(shù)的分布形式,利用Monte-Carlo方法生成裂隙網(wǎng)絡(luò),并通過優(yōu)化方法剔除對導(dǎo)水無用的裂隙。在生成的過程中分別對每組裂隙中的每條裂隙進行編號。

1.2 一維管單元模型

平面滲流假定認為裂隙為有厚度的圓盤,水在整個裂隙圓盤上流動。進行滲流計算時需對裂隙圓盤進行網(wǎng)格劃分,裂隙數(shù)目較多時,有限元計算工作量非常地大,給研究工作帶來很多不便。本研究根據(jù)前人的溝槽流假設(shè),將水在裂隙中的滲流路徑簡化為管單元模型,即把相交裂隙的圓心與交線的中心作為管單元節(jié)點,兩個相鄰節(jié)點組成一個管單元,將裂隙之間水的流動簡化為管元通道之間的流動。根據(jù)管元所在裂隙平面的組號再對管元進行分類,如圖1所示。

圖1 裂隙網(wǎng)絡(luò)管單元Fig.1 Pipe element of fracture network

1.3 滲流的計算法

假設(shè)管單元的水流通道為圓形,那么通過水流通道的流量[18]為:

(1)

式(1)中:ΔH為兩個節(jié)點之間的水頭差;L為兩節(jié)點之間的距離;D為管單元管徑;ρ為水的密度;g為重力加速度;μ為水動力黏滯系數(shù)(20 ℃時,μ=0.01 g/(cm·s))。

根據(jù)式(1)可以得到任意的管單元m中流過兩個節(jié)點i,j的流量為:

(2)

式(2)中,節(jié)點流量q以流入為正,流出為負。用矩陣可表示為:

(3)

由式(3)集成為總體滲透方程:

[K]{H}={Q},

其中,[K]稱為總滲透矩陣，它由所有的單元滲透矩陣[Km]疊加得到。

2 管徑反演與滲透張量確定

2.1 參數(shù)的反演

由式(1)可知,進行滲流計算須知道每個管單元直徑D的大小。根據(jù)等效管徑與裂隙隙寬的正比關(guān)系,設(shè)第i組裂隙的等效管徑與隙寬的比值為Ni,i=1,2,3,…,n;n為研究域中裂隙組編號。

基于改進后的遺傳算法,以滲透張量的模擬值與實測值之間的誤差最小為目標函數(shù),通過反演的方法整體上得到每組裂隙的等效管徑與隙寬的比值Ni。由于Ni的改變可以影響滲透張量主值的大小,所以可以將目標函數(shù)表示為:

(4)

2.2 滲透張量的確定

節(jié)理巖體滲透性具有明顯的各向異性,即各個方向上的滲透系數(shù)變化明顯,目前大多采用滲透張量來表示。依據(jù)文獻[19]中提到的在單位體積立方體的空間中,在兩個相對面邊界上施加水頭邊界,對其進行數(shù)值仿真試驗,可得到各方向的流量qx、qy和qz,但僅靠這3個數(shù)值并不能有效確定裂隙的滲透張量,所以在正六面體試件的某一對相對面(邊界)施加等水頭邊界,其余各面根據(jù)其相鄰的情況設(shè)定相應(yīng)的變水頭邊界條件,分別改變3組相對面的邊界條件,如圖2所示的計算模型進行試驗。

圖2 計算模型的邊界條件Fig.2 Boundary conditions of computational model

模型邊界沿ii方向施加等水頭邊界,另兩個方向為變水頭邊界,可得到6個流量數(shù)據(jù);進行3組仿真可得到18個數(shù)據(jù),然后基于迭加原理確定二階滲透張量。根據(jù)相應(yīng)的達西定律:

Q=kAΔP。

(5)

式(5)中:Q為流量向量;k為滲透張量;ΔP為水力坡降;A為過流面積。

由相應(yīng)的流量分量,通過矩陣變換得到滲透張量矩陣:

(6)

3 算 例

3.1 工程數(shù)據(jù)參照

文獻[11]通過實測得到裂隙幾何參數(shù)的信息如表1所示,共測得3組裂隙。

表1 裂隙網(wǎng)絡(luò)參數(shù)統(tǒng)計Table 1 Statistical parameters of fracture network

文獻[11]利用邊界元法算得表1裂隙網(wǎng)絡(luò)的滲透張量,并通過現(xiàn)場壓水試驗進行校核得到滲透主值及主方向如式(7),作為后面管徑反演的目標滲透張量。

(7)

3.2 多參數(shù)模型計算

采用自編裂隙網(wǎng)絡(luò)有限元程序,將表1中信息輸入到程序中,對裂隙網(wǎng)絡(luò)進行模擬。程序生成100 m×100 m×100 m的立方體區(qū)域,在其中截取20 m×20 m×20 m的立方體區(qū)域作為研究域。利用Monte-Carlo方法生成三維裂隙網(wǎng)絡(luò)。在生成域中共生成結(jié)構(gòu)面200 255條,其中第1組結(jié)構(gòu)面129 006條,第2組結(jié)構(gòu)面38 224條,第3組結(jié)構(gòu)面33 025條。在研究域中,剔除對滲流無作用的裂隙后共生成結(jié)構(gòu)面911條,其中第1組結(jié)構(gòu)面為480條,第2組結(jié)構(gòu)面為250條,第3組結(jié)構(gòu)面為181條,將其簡化成管單元模型,則得到4 104個單元,節(jié)點數(shù)為3 033個,生成的裂隙網(wǎng)絡(luò)如圖3所示。

圖3 三維裂隙網(wǎng)絡(luò)生成Fig.3 Generation of 3-D fracture network

在數(shù)值計算過程中,設(shè)定水頭邊界H1=30 m,H2=10 m。由于模型中生成了3組裂隙網(wǎng)絡(luò),根據(jù)2節(jié)內(nèi)容所述,設(shè)置參數(shù)N1、N2和N3分別作為3組裂隙的等效管徑與隙寬的比值。以式(7)作為擬合的標準,通過優(yōu)化擬合,得到3個參數(shù):13.17、12.06及3.01。

通過計算得到多參數(shù)模型水頭分布如圖4所示,滲透張量(m/s)如下:

(8)

圖4 三組模擬水頭分布圖Fig.4 Three groups of simulated head distribution

通過滲透張量得到的滲透主值和滲透主方向如下:

(9)

與式(7)數(shù)據(jù)相比,誤差為:

(10)

式(10)中:ε為主值相對誤差;θ為主方向絕對誤差。

3.3 單參數(shù)模型計算

采用文獻[16]的單參數(shù)模型方法,對三維裂隙網(wǎng)絡(luò)滲流模型進行數(shù)值模擬。在計算的過程中設(shè)置一個參數(shù)N作為3組裂隙等效管徑與隙寬的比值。在表1所示信息的情況下,得到參數(shù)N的值為12.34,其滲透張量結(jié)果如下:

(11)

滲透主系數(shù)與滲透主方向如下:

(12)

與式(7)數(shù)據(jù)相比,其誤差為:

(13)

式(13)中:ε為主值相對誤差;θ為主方向絕對誤差。

與多參數(shù)模型相比,單參數(shù)模型滲透主值的相對誤差較大。

3.4 隙寬為定值時不同模型的計算

3.2和3.3兩節(jié)中隙寬的分布服從均值為0.000 1,標準差為0.000 1的對數(shù)正態(tài)分布,本節(jié)將裂隙的所有隙寬設(shè)置為定值0.000 1,其他幾何參數(shù)保持不變。研究域尺寸為20 m×20 m×20 m的情況下多參數(shù)模型的參數(shù)值為19.51、16.55和14.05,滲透張量為:

(14)

(15)

(16)

單參數(shù)模型的參數(shù)值為16.46,滲透張量為:

(17)

(18)

(19)

改變研究域的尺寸大小得到多參數(shù)模型的參數(shù)均值與單參數(shù)模型的參數(shù),如表2所示。

表2 模型參數(shù)值Table 2 Parameter values of the model

從式(16)和式(19)的計算結(jié)果可以看出,兩種模型在定隙寬的情況下,滲透主值和滲透主方向的誤差在工程允許范圍內(nèi)。

3.5 分 析

通過3.2和3.3節(jié)的計算結(jié)果可以看出,當(dāng)隙寬服從對數(shù)正態(tài)分布時,多參數(shù)模型主方向的絕對誤差總和為63.04,單參數(shù)模型的誤差和為63.1,兩者相當(dāng);在滲透主值誤差方面,多參數(shù)模型的誤差明顯小于單參數(shù)模型。在實際情況中,不同的裂隙組之間的管寬比存在差異,多參數(shù)模型更加與實際情況相符合,所以在計算結(jié)果上比單參數(shù)模型更為精確。從兩種模型的計算結(jié)果可以發(fā)現(xiàn),主方向上的差異比較明顯,這是因為在優(yōu)化的過程中,在使目標函數(shù)式(4)達到要求的條件下,對應(yīng)主方向與目標張量的主方向保持基本一致,即兩者夾角為銳角即可。

當(dāng)裂隙網(wǎng)絡(luò)的隙寬為定值時,單參數(shù)模型的計算結(jié)果優(yōu)于多參數(shù)模型,這是因為當(dāng)隙寬為恒定值時,每個單元的管徑與隙寬的比值接近于同一個數(shù),所以采用單參數(shù)模型能得到更精確的結(jié)果。同時,由表2還發(fā)現(xiàn)恒定隙寬情況下,多參數(shù)模型的參數(shù)均值與單參數(shù)模型參數(shù)值近似相等,且在研究域尺寸不同的情況下參數(shù)值在17.00左右。

4 結(jié) 論

本研究采用三維離散裂隙網(wǎng)絡(luò)模型,將復(fù)雜的裂隙網(wǎng)絡(luò)簡化為管單元模型,在計算方面設(shè)置多個參數(shù)得到滲透張量,從算例的計算結(jié)果可以看出:當(dāng)隙寬服從對數(shù)正態(tài)分布時多參數(shù)模型誤差較小,這種情況下適宜選擇多參數(shù)模型;當(dāng)隙寬均值為恒定大小時,適宜采用單參數(shù)模型進行計算分析。同時,從式(10)的計算結(jié)果可以發(fā)現(xiàn),多參數(shù)模型的計算結(jié)果與工程實測結(jié)果之間的誤差比較小,有著較高的計算精度。這表明在工程應(yīng)用上,使用多參數(shù)模型進行模擬可以提高滲透分析的準確性,因而具有良好的理論與實用價值。

[1] JAVADI M, SHARIFZADEH M, SHAHRIAR K. A new geometrical model for non-linear fluid flow through rough fractures[J].Journal of Hydrology,2010,389(1):18.

[2] 王者超,張振杰,李術(shù)才,等.基于離散裂隙網(wǎng)絡(luò)法的地下石油洞庫洞室間水封性評價[J].山東大學(xué)學(xué)報(工學(xué)版),2016,46(2):94.

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

[4] WITHERSPOON P, LONG J. Some recent development in understanding the hydrology of fractured rocks[C]//The 28th US Symposium on Rock Mechanics.Tucson:CRC Press,1987:422.

[5] ANSARI M A, DEODHAR A, KUMAR U S, et al. Isotope hydrogeological factors control transport of radon-222 in hard rock fractured aquifer of bangalore, karnataka[M].Berlin:Springer International Publishing,2016:259.

[6] 周志芳,莊超,戴云峰,等.單孔振蕩式微水試驗確定裂隙巖體各向異性滲透參數(shù)[J].巖石力學(xué)與工程學(xué)報,2015,34(2):271.

[7] 楊金保,李蕾,謝妮.考慮應(yīng)力狀態(tài)的裂隙巖體滲透系數(shù)張量計算[J].人民黃河,2013,35(3):133.

[8] 李克克,張楊松.基于實測結(jié)構(gòu)面的巖體滲流特性分析[J].勘察科學(xué)技術(shù),2016(1):1.

[9] 王培濤,楊天鴻,于慶磊,等.基于離散裂隙網(wǎng)絡(luò)模型的節(jié)理巖體滲透張量及特性分析[J].巖土力學(xué),2013,34(增刊2):448.

[10] BAGHBANAN A, JING L. Hydraulic properties of fractured rock masses with correlated fracture length and aperture[J].International Journal of Rock Mechanics and Mining Sciences,2007,44(5):704.

[11] 李新強,陳祖煜.三維裂隙網(wǎng)絡(luò)滲流計算的邊界元法及程序[J].中國水利水電科學(xué)研究院學(xué)報,2006,49(2):86.

[12] LEE I H, NI C F. Fracture-based modeling of complex flow and CO2migration in three dimensional fractured rocks[J].Computer and Geosciences,2015,81:64.

[13] CACAS M C, LEDOUX E, MARSILY G D, et al. Modeling fracture flow with a stochastic discrete fracture network : calibration and validation: 1. the flow model[J].Water Resources Research,1990,26(3):479.

[14] WANG M, KULATILAKE P H S W, UM J, et al. Estimation of REV size and three-dimensional hydraulic conductivity tensor for a fractured rock mass through a single well packer test and discrete fracture fluid flow modeling[J].International Journal of Rock Mechanics and Mining Sciences,2002,39(7):887.

[15] 于青春,武雄,大西有三.非連續(xù)裂隙網(wǎng)絡(luò)管狀滲流模型及其校正[J].巖石力學(xué)與工程學(xué)報,2006,25(7):1469.

[16] 王俊奇,岳瀟,李闖.三維裂隙網(wǎng)絡(luò)管單元模型確定巖體滲透張量[J].華北電力大學(xué)學(xué)報(自然科學(xué)版),2015,42(1):108.

[17] 張發(fā)明,汪小剛,賈志欣,等.三維結(jié)構(gòu)面連通率的隨機模擬計算[J].巖石力學(xué)與工程學(xué)報,2004,23(9):1487.

[18] 汪小剛,賈志欣,張發(fā)明,等.巖體結(jié)構(gòu)面網(wǎng)格模擬原理及其工程應(yīng)用[M].北京:水利水電出版社,2010:83.

[19] ZHANG X, SANDERSON D J. Numerical modeling and analysis of fluid flow and deformation of fractured rock masses[M].Amsterdam: Elsevier Science Ltd,2002:53.

Studyontheseepagetensorofmulti-parametermodelofpipeelementsinfracturedrockmass

WANG Junqi, LI Pengyu, DONG Ye, MU Mengjing

(Rene Wable Energy School, North China Electric Power University, Beijing 102206, China)

Based on the fracture occurrence of rock mass measured on the spot, the Monte-Carlo method is employed to generate three groups of fracture networks by referring to the discfracture hypothesis. In accordance with the hypothesis of trench flow in the fracture, the 3-D fracture network is simplified to a model consisting of one-dimensional pipe elements, and the equivalent diameter of the pipe elements is considered to be linear with the fracture width. The ratio of equivalent diameters to apertures of each fracture being set as the unknown parameter, the equivalent seepage tensor is obtained by inverting the genetic algorithm,which can be compared with the equivalent seepage tensor obtained by the whole single parameter model. In comparison with the field data of seepage tensor, the results show that the multi-parameter model can better characterize the seepage tensor than the single parameter model.

multi-parameter model; 3-D fracture network; pipe element; seepage tensor

TU452

A

1671-8798(2017)05-0358-08