巖體裂隙網(wǎng)絡(luò)滲流變水溫影響分析

徐維生 ，周創(chuàng)兵

（1. 武漢大學(xué) 水資源與水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430072；2. 武漢大學(xué) 水利水電學(xué)院，武漢 430072；3. 長(zhǎng)江科學(xué)院 水利部巖土力學(xué)與工程重點(diǎn)實(shí)驗(yàn)室，武漢 430010）

1 引 言

滲流問(wèn)題是巖土工程中常見(jiàn)的問(wèn)題，當(dāng)巖體裂隙較稀疏，巖塊滲透系數(shù)非常小時(shí)，應(yīng)采用巖體離散裂隙網(wǎng)絡(luò)模型。裂隙網(wǎng)絡(luò)滲流模型是以線性立方定律為基礎(chǔ)發(fā)展起來(lái)的，然而巖體裂隙滲流本構(gòu)關(guān)系并不一定全是線性的[1－2]。以往對(duì)于這類(lèi)滲流本構(gòu)非線性問(wèn)題大多按線性簡(jiǎn)化，當(dāng)滲流偏離線性滲流較遠(yuǎn)時(shí)按線性簡(jiǎn)化會(huì)帶來(lái)較大誤差[3－5]。

目前，考慮水流溫度影響的主要有連續(xù)介質(zhì)（等效連續(xù)）[6－12]和裂隙網(wǎng)絡(luò)介質(zhì)模型[13－15]。等效連續(xù)介質(zhì)模型主要適合于裂隙分布較密集、集中，表征單元體較小的情況；而裂隙網(wǎng)絡(luò)介質(zhì)模型主要適合于滲流由主要的大稀疏裂隙控制的滲流[16]。

關(guān)于裂隙變水溫流動(dòng)，目前對(duì)于單裂隙內(nèi)變水溫影響大都將裂隙內(nèi)水溫按平均溫度處理，這樣單裂隙內(nèi)水頭分布為直線，也即裂隙水流與水力坡降成線性關(guān)系。然而裂隙滲流本構(gòu)關(guān)系與水流溫度相關(guān)，單裂隙變水溫滲流裂隙水頭分布為曲線，也即裂隙變溫水流與水力坡降成非線性關(guān)系[2]。本文從水溫對(duì)裂隙水流滲透系數(shù)影響出發(fā)，以水流黏滯系數(shù)與溫度關(guān)系試驗(yàn)的經(jīng)驗(yàn)公式為基礎(chǔ)建立模型，推導(dǎo)相關(guān)公式，分析變水溫對(duì)單裂隙及裂隙網(wǎng)絡(luò)滲流的影響。

2 單裂隙變溫水流滲流場(chǎng)分析

2.1 立方定律

對(duì)于不考慮水流溫度場(chǎng)變化的單裂隙水流，水流滿足立方定律：

式中：q為單寬流量(m3/(s·m))；v為滲流速度(m/s)；b為裂隙寬度(m)；k為滲透系數(shù)(m/s)；υ為運(yùn)動(dòng)黏滯系數(shù)（m2/s）；J為水力坡降。

2.2 單裂隙變水溫水流滲流場(chǎng)分析

2.2.1 單裂隙變水溫水流滲流場(chǎng)近似解析解

考慮水流溫度場(chǎng)變化的裂隙水流，可以得出水流在裂隙有限小長(zhǎng)度內(nèi)也滿足立方定律：

式中：qfin、vfin、bfin、Kfin、Jfin、υfin分別為有限小長(zhǎng)度內(nèi)裂隙單寬流量(m3/(s·m))、裂隙流速(m/s)、裂隙隙寬(m)、裂隙滲透系數(shù)(m/s)、裂隙水力坡降、裂隙水流運(yùn)動(dòng)黏滯系數(shù)（m2/s）。

根據(jù)立方定律可知，裂隙滲透系數(shù)與水流運(yùn)動(dòng)黏滯系數(shù)成反比，而水流運(yùn)動(dòng)黏滯系數(shù)又是溫度的函數(shù)[1]。目前，廣泛采用以下經(jīng)驗(yàn)公式：

式中：t為水溫（℃）。

（1）單裂隙線性變水溫分布規(guī)律

水流溫度在單裂隙內(nèi)成線性分布，且裂隙兩端點(diǎn)溫度不相等，則有

式中：x1、x2為點(diǎn)1、2的坐標(biāo)值；t1、t2為點(diǎn)1、2的水流溫度值；x為長(zhǎng)度（m）。

（2）單裂隙線性變水溫流量公式

將水頭表示成溫度場(chǎng)的函數(shù)，q為裂隙的單寬流量，裂隙滲透系數(shù)k是溫度t的函數(shù)，運(yùn)動(dòng)黏滯系數(shù)υ(t)是溫度值得函數(shù)，即

式中：H為水頭值(m)。

式中：H1、H2為點(diǎn) 1、2的水頭值；k(t(x))為滲透系數(shù)，是水流溫度t的函數(shù)，t又是坐標(biāo)x的函數(shù)。

令單裂隙變溫水流等效運(yùn)動(dòng)黏滯系數(shù)為υv，計(jì)算公式為

求流量的關(guān)鍵在于求式（11）中積分項(xiàng)，查積分公式可得

對(duì)應(yīng)積分公式可得，a=0.000 221，b=0.033，c=1，b2＞4ac。

（3）單裂隙線性變水溫水流水頭公式

裂隙中水溫值為t時(shí)的水頭為H，可推得

（4）單裂隙線性變水溫水流水頭分布規(guī)律

2.2.2 單裂隙變水溫水流滲流場(chǎng)分析的有限元解

對(duì)于考慮溫度場(chǎng)影響的單裂隙線性立方水流，通過(guò)第n段的單寬流量為qn，則有

式中：Jn為第n段水力坡降；bn為第n段平均隙寬；υn為裂隙第n段平均運(yùn)動(dòng)黏滯系數(shù)；Hn-1、Hn分別為第n-1、第n點(diǎn)水頭；qn為第n段單寬流量；Δl 為裂隙段長(zhǎng)度。

將式（20）中所有等式左右兩邊相加及將式前n（ n = 1,2,...,N - 1）項(xiàng)等式相加，可得

聯(lián)立式（21）、（22）求解，可得線性立方流第n裂隙段上水頭為

式中：Hn為節(jié)點(diǎn)n水頭。

求解式（21）裂隙單寬流量q為

2.3 單裂隙滲流場(chǎng)水溫影響分析

以一100 m長(zhǎng)單裂隙為研究對(duì)象，分工況對(duì)裂隙水流進(jìn)行分析，裂隙示意圖見(jiàn)圖1。分兩種工況，工況表見(jiàn)表1。

圖1 裂隙示意圖Fig.1 Sketch map of fracture

表1 計(jì)算工況Table 1 Computing conditions

計(jì)算結(jié)果見(jiàn)表2，由表2結(jié)果可得圖2。

表2 計(jì)算結(jié)果Table 2 Computed Results

圖2 單裂隙變水溫水流水頭值Fig.2 Water head in single fracture with variable temperature fluid

由表2、圖2可以看出：?jiǎn)瘟严秲?nèi)水流由高溫流向低溫時(shí)，水頭分布曲線為凸曲線，此時(shí)按線性滲流簡(jiǎn)化滲流水頭整體偏小；單裂隙內(nèi)水流由低溫流向高溫時(shí)，水頭分布曲線為凹曲線，此時(shí)按線性滲流簡(jiǎn)化滲流水頭整體偏大。當(dāng)裂隙內(nèi)節(jié)點(diǎn)多到一定程度就非常接近解析解。解析解和子結(jié)構(gòu)解吻合較好，解析解與子結(jié)構(gòu)解結(jié)果非常接近，有些節(jié)點(diǎn)數(shù)值相等，而與不考慮變溫作用的水頭值相差較大。

3 裂隙網(wǎng)絡(luò)滲流溫度影響子結(jié)構(gòu)法

裂隙內(nèi)水溫不同，裂隙內(nèi)不同點(diǎn)的水流黏滯系數(shù)不同，則水流滲透系數(shù)亦不同。因此，當(dāng)裂隙較長(zhǎng)、溫度變化較大時(shí)，有必要分析裂隙內(nèi)部水頭分布趨勢(shì)。

（1）裂隙網(wǎng)絡(luò)變水溫滲流近似解析解法

式中：TLJ為水力傳導(dǎo)矩陣；T LJ =diag(Tlj1,Tlj2,???,T ljm)m×m，A稱(chēng)為裂隙網(wǎng)絡(luò)的n×m階銜接矩陣。該矩陣中的元素aij：當(dāng)j線元不銜接i節(jié)點(diǎn)，取0；當(dāng)j線元銜接i節(jié)點(diǎn)，指向離開(kāi)i節(jié)點(diǎn)，取-1，指向i節(jié)點(diǎn)取1。

式中：Bj為裂隙網(wǎng)絡(luò)中第j條裂隙的等效水力隙寬；υvj為第j條裂隙的水流等效運(yùn)動(dòng)黏滯系數(shù)，裂隙兩端點(diǎn)溫度不相等時(shí)計(jì)算方法見(jiàn)式（9），相等時(shí)見(jiàn)式（4）。求得裂隙節(jié)點(diǎn)水頭后，可根據(jù)單裂隙滲流近似解析法求解裂隙流量和裂隙內(nèi)任意一點(diǎn)水頭值，具體求法見(jiàn)式（23）、（24）。

（2）裂隙網(wǎng)絡(luò)變水溫滲流子結(jié)構(gòu)法

利用裂隙等效水力隙寬公式和裂隙單寬流量公式先求得裂隙網(wǎng)絡(luò)中的各裂隙節(jié)點(diǎn)水頭，再根據(jù)式（23）計(jì)算每條裂隙段的任一子單元水頭和流速。裂隙節(jié)點(diǎn)水頭公式為

式中：TLD為等效水力傳導(dǎo)矩陣；T LD=diag(Tld1,Tld2,? ??,T ldm)m×m，且有

式中：bjd為裂隙網(wǎng)絡(luò)中第 j條裂隙第d段（子單元）的平均隙寬；υjd為裂隙網(wǎng)絡(luò)中第j條裂隙第d段（子單元）的平均運(yùn)動(dòng)黏滯系數(shù)。

求得裂隙的節(jié)點(diǎn)水頭之后，根據(jù)式（23）、（24）分別計(jì)算裂隙的等分點(diǎn)（子單元端點(diǎn)）水頭和單寬流量，然后求裂隙子單元流速。

4 二維巖體裂隙網(wǎng)絡(luò)變溫滲流分析

選擇一200 m×200 m研究區(qū)域，區(qū)域內(nèi)有5條裂隙見(jiàn)圖3，分別為裂隙BH（裂隙1），裂隙HE（裂隙2），裂隙HF（裂隙3），裂隙HD（裂隙4），裂隙HG（裂隙5），AC、EF為定水頭邊界，CE、AF為零流量邊界。裂隙2的10等分長(zhǎng)度為L(zhǎng)。

圖3 簡(jiǎn)單裂隙網(wǎng)絡(luò)Fig.3 Simple fracture network

分7種工況進(jìn)行滲流計(jì)算，具體見(jiàn)表3。

表3 工況表Table 3 Conditions table

表4為工況1、2條件下裂隙2、3內(nèi)的水頭分布值。

表4 工況1、2 裂隙2、3內(nèi)水頭分布Table 4 Water head distribution in fracture No. 1 and No. 2 in condition No. 1 and No. 2

從表4、圖4和圖5計(jì)算的工況1、2條件下的裂隙2、3的水頭分布可以看出：?jiǎn)瘟严秲?nèi)，高水溫地方水力坡降小，低水溫地方水力坡降大。對(duì)于裂隙2，水從下往上流，是高溫流向低溫，裂隙2的水頭分布曲線為凸曲線；對(duì)于裂隙 3，水從上往下流，是低溫流向高溫，裂隙3的水頭分布曲線為凹曲線。裂隙水流整體溫度值越高，裂隙各節(jié)點(diǎn)整體水頭值越大。驗(yàn)證了文中前部分解析解是合理的。

表5為各工況下裂隙2、3流量值。由表可以看出，裂隙存在偏流效應(yīng)，即節(jié)點(diǎn)水流有偏向水流溫度高的區(qū)域流動(dòng)的趨勢(shì)；溫度梯度越大，流量相差越大，水頭值相差越大；反之溫度梯度越小，則流量相差越小，水頭值相差越小。

表5 各工況裂隙2、3流量表Table 5 Flows of fracture No.2 and No.3 under various conditions

5 結(jié) 論

（1）單裂隙變水溫水流內(nèi)水流水頭與水力坡降成非線性關(guān)系，當(dāng)水流由高溫向低溫流動(dòng)時(shí)，水頭分布曲線為凸曲線，此時(shí)按線性滲流簡(jiǎn)化滲流水頭值整體偏??；當(dāng)由低溫向高溫流動(dòng)時(shí)，水頭分布曲線為凹曲線，此時(shí)按線性滲流簡(jiǎn)化滲流水頭值整體偏大。

（2）單裂隙內(nèi)，高水溫處水力坡降小，低水溫處水力坡降大；裂隙網(wǎng)絡(luò)內(nèi)裂隙平均水溫越高，流速越快；裂隙網(wǎng)絡(luò)內(nèi)存在與裂隙寬度相似的溫度偏流效應(yīng)，即節(jié)點(diǎn)水流有偏向水流溫度高的流動(dòng)趨勢(shì)；溫度梯度越大，流量、水頭值相差越大；反之溫度梯度越小，則流量、水頭值相差越小。

[1]柴軍瑞, 徐維生. 大壩工程滲流非線性問(wèn)題[M]. 北京:中國(guó)水利水電出版社, 2010.

[2]徐維生. 高壩工程滲流力學(xué)中的幾類(lèi)非線性問(wèn)題研究[博士學(xué)位論文D]. 西安: 西安理工大學(xué), 2011.

[3]徐維生, 柴軍瑞, 陳興周, 等. 巖體裂隙網(wǎng)絡(luò)非線性非立方滲流研究與應(yīng)用[J]. 巖土力學(xué), 2009, 30(增刊 1):53－57.XU Wei-sheng, CHAI Jun-rui, CHEN Xing-zhou, et al.Study of nonlinear non-cubic seepage in fracture network of rock mass and its application[J]. Rock and Soil Mechanics, 2009, 30(Supp. 1): 53－57.

[4]徐維生, 柴軍瑞, 陳興周, 等. 巖體粗糙裂隙網(wǎng)絡(luò)滲流分析的子單元法[J]. 水動(dòng)力學(xué)研究與進(jìn)展(A輯), 2010,25(2): 175－182.XU Wei-sheng, CHAI Jun-rui, CHEN Xing-zhou, et al.The sub-element method on the analysis for the seepage in the coarse fractures network of rock mass[J]Chinese Journal of Hydrodynamics (Seris A), 2010, 25(2): 175－182.

[5]徐維生, 柴軍瑞, 陳興周. 含自由面巖體裂隙網(wǎng)絡(luò)非線性滲流分析的傳導(dǎo)矩陣調(diào)整法[J]. 水動(dòng)力學(xué)研究與進(jìn)展(A輯), 2010, 25(4): 484－482.XU Wei-sheng, CHAI Jun-rui, CHEN Xing-zhou. The conductivity matrix reforming method about non-linear seepage with free surface in fracture network of rock mass[J]. Chinese Journal of Hydrodynamics (Seris A),2010, 25(4): 484－492.

[6]楊更社, 周春華, 田應(yīng)國(guó). 寒區(qū)軟巖隧道的水熱耦合數(shù)值模擬與分析[J]. 巖土力學(xué), 2006, 26(8): 2258－1262.YANG Geng-she, ZHOU Chun-hua, TIAN Ying-guo.Numerical simulation and analysis of moisture-heat coupling for soft rock tunnel in cold regions[J]. Rock and Soil Mechanics, 2006, 26(8): 2258－1262.

[7]張學(xué)富, 喻文兵, 劉志強(qiáng). 寒區(qū)隧道滲流場(chǎng)和溫度場(chǎng)耦合問(wèn)題的三維非線性分析[J]. 巖土工程學(xué)報(bào), 2006,28(9): 1095－1100.ZHANG Xue-fu, YU Wen-bing, LIU Zhi-qiang.Three-dimensional nonlinear analysis for coupled problem of seepage field and temperature field of cold regions tunnels[J]. Chinese Journal of Geotechnical Engineering, 2006, 28(9): 1095－1100.

[8]盛金昌. 多孔介質(zhì)流-固-熱三場(chǎng)全耦合數(shù)學(xué)模型及數(shù)值模擬[J]. 巖石力學(xué)與工程學(xué)報(bào), 2006, 25(增刊 1):3030－3033.SHENG Jin-chang. Fully coupled thermo-hydromechanical model of saturated porous media and numerical modelling[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(Supp.1): 3030－3033.

[9]陳波, 李寧, 糕瑞花. 多孔介質(zhì)的變形場(chǎng). 滲流場(chǎng)-溫度場(chǎng)耦合有限元分析[J]. 巖石力學(xué)與工程學(xué)報(bào), 2001,20(4): 467－472.CHEN Bo, LI Ning, GAO Rui-hua. FEM analysis on fully coupled thermo-hydro-mechanic behavior of porous media[J]. Chinese Journal of Rock Mechanics and Engineering, 2001, 20(4): 467－472.

[10]GUO R P. Thermo-hydro-mechanical modelling of the buffer/container experiment[J]. Engineering Geology,2011, 122(3－4): 303－315.

[11]ENGELER I, FRANSSEN H J H, MULLER R, et al. The importance of coupled modelling of variably saturated groundwater flow-heat transport for assessing river-aquifer interactions[J]. Journal of Hydrology, 2011,397(3－4): 295－305.

[12]FREI S, FLECKENSTEIN J H, KOLLET S J, et al.Patterns and dynamics of river-aquifer exchange with variably-saturated flow using a fully-coupled model[J].Journal of Hydrology, 375(3－4): 383－393.

[13]王如賓, 柴軍瑞. 非連續(xù)裂隙巖體滲流場(chǎng)與溫度場(chǎng)耦合分析[J]. 工程勘察, 2006, (6): 18－21.WANG Ru-bin, CHAI Jun-rui. Analysis of coupling seepage field and temperature field in discrete fracture rock mass[J]. Journal of Geotechnical Investigation &Surveying, 2006, (6): 18－21.

[14]ZHANG Y J, YANG C S. Coupled thermo-hydromechanical-migratory model for dual-porosity medium and numerical analysis[J]. Journal of Central South University of Technology, 2011, 18(4): 1256－1262.

[15]ZHANG Y J, ZHANG W Q. 3D thermo-hydromechanical-migratory coupling model and FEM analyses for dual-porosity medium[J]. Science China—Technological Sciences, 2010, 53(8): 2172－2182.

[16]仵彥卿. 巖土水力學(xué)[M]. 北京: 科學(xué)出版社, 2009: 6－7.