賀香蘭，周佳慶，魏 凱，王 敏

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

0 引 言

裂隙是巖體滲流的主要通道，控制著巖體的滲流特性，隨著高壩樞紐工程和抽水蓄能電站的興建，巖體可能承受高滲壓、高水力梯度作用，此時節(jié)理裂隙等結(jié)構(gòu)面易張開擴(kuò)展，甚至發(fā)生水力劈裂，巖體滲流往往表現(xiàn)出非達(dá)西特征，特別是對于裂隙充分發(fā)育的破碎巖體，滲流特性更為復(fù)雜。黃先伍等[1]、李順才等[2]在破碎巖石構(gòu)成的多孔介質(zhì)滲流試驗的基礎(chǔ)上，研究了非達(dá)西滲流特性與孔隙率的關(guān)系；Thauvin和Mohanty[3]利用數(shù)值模擬的方法從細(xì)觀結(jié)構(gòu)的角度研究了多孔介質(zhì)非達(dá)西滲流的成因；胡大偉等[4]研究了單軸壓縮破壞的大理巖峰后非線性滲流特征演化規(guī)律。

對于非達(dá)西流，傳統(tǒng)的達(dá)西滲流理論已不能準(zhǔn)確表征巖體的滲透特性，因此滲流流態(tài)的判斷和分析理論的選擇是研究介質(zhì)滲流的關(guān)鍵。為此，國內(nèi)外學(xué)者對達(dá)西流與非達(dá)西流的臨界條件開展了大量研究，主要采用雷諾數(shù)Re和福希海默數(shù)F0作為臨界參數(shù)[5]。孟如真等[6]、李健等[7]分別通過現(xiàn)場壓水試驗和室內(nèi)滲流試驗，得到了量級為10～100的臨界雷諾數(shù)；Zimmerman等[8]通過對單裂隙砂巖進(jìn)行滲流試驗和數(shù)值模擬，發(fā)現(xiàn)Re大于20時，滲流滿足Forchheimer非達(dá)西方程；Javadi等[9]開展了粗糙裂隙面滲流試驗，得到的臨界雷諾數(shù)介于0.001～25；同時還有文獻(xiàn)提出了高達(dá)2 300的臨界雷諾數(shù)[10]。另一方面，Zeng和Grigg等[11]根據(jù)福希海默數(shù)F0的定義及其物理意義，建議取F0=0.11作為達(dá)西流和非達(dá)西流的臨界值；隨后，Macini等[12]通過研究天然砂粒的滲流試驗，得到值為0.4左右的臨界福希海默數(shù)，該值對應(yīng)28%的非達(dá)西效應(yīng)；文獻(xiàn)[13]在進(jìn)行數(shù)值模擬時采用了F0=0.01作為臨界值。通過眾多研究可以發(fā)現(xiàn)，雷諾數(shù)Re和福希海默數(shù)F0臨界值的取值范圍變化較大，且目前尚未得到臨界值的變化規(guī)律，因此采用參數(shù)Re和F0判斷滲流流態(tài)存在一定的主觀性和不確定性。

綜上所述，盡管滲流的非達(dá)西性得到了大量的研究，但涉及破碎巖體的非達(dá)西滲流的研究成果相對較少，并且達(dá)西流和非達(dá)西流的臨界問題仍未得到有效地解決?；诖?，本文通過開展破碎花崗巖在不同圍壓條件下的高水力梯度非達(dá)西滲流試驗，研究破碎巖體復(fù)雜裂隙介質(zhì)的滲流特性，并在試驗結(jié)果基礎(chǔ)上，提出一種簡單直觀方便的非達(dá)西流判定方法。該方法可將臨界值的不確定性問題轉(zhuǎn)化為實際工程目標(biāo)需求的容許誤差，具有較好的實用價值。

1 試驗系統(tǒng)與步驟

1.1 試驗系統(tǒng)介紹

破碎巖體非達(dá)西滲流試驗在武漢大學(xué)與法國里爾科技大學(xué)(USTL)聯(lián)合研發(fā)的Triaxial Cell三軸試驗系統(tǒng)上進(jìn)行，見圖1。該系統(tǒng)由圍壓、軸壓和孔隙水壓3套獨(dú)立加載部分構(gòu)成，具有應(yīng)力、應(yīng)變、位移和流量4種控制加載方式，可施加最大圍壓60 MPa，最大軸壓375 MPa，最大滲透水壓力50 MPa。其中滲透水壓力由升級后的高精度液壓泵PHMP 50-1000型施加，該流量泵可容納100 mL體積水，流量范圍為0.000 41～66.8 mL/min。系統(tǒng)出口端接入測量精度為0.01 g的SA2202S-CW 型高精度自動電子天平，以測量出口水流流量。

1.2 試驗步驟

試驗巖樣取自廣東省陽江抽水蓄能電站高壓隧洞區(qū)，為中?；◢弾r，經(jīng)取芯機(jī)、剛性切割機(jī)和打磨機(jī)加工后制備得到Φ50 mm×100 mm的圓柱狀標(biāo)準(zhǔn)巖石試樣，見圖2。破碎巖體滲流試驗過程如下。

(1)將飽和的標(biāo)準(zhǔn)巖石試樣在預(yù)定圍壓條件下施加軸向壓力，直至試樣破壞，得到破碎巖樣。

(2)以位移加載模式控制巖樣保持在破碎初始狀態(tài)，并以0.1 MPa的滲透水壓持續(xù)進(jìn)行12 h的過流，保證試樣充分飽和。

(3)在設(shè)定的圍壓條件下，逐級施加滲透水壓(見圖3)，測量并記錄每級水壓穩(wěn)定后的出口水流流量，完成該圍壓下的滲流試驗后，改變圍壓，重復(fù)上述過程，直至圍壓達(dá)到30 MPa。

圖1 非達(dá)西滲流試驗裝置Fig.1 Sketch of non-Darcy flow test equipment

圖2 試驗前后花崗巖試樣Fig.2 Granite samples before and after test

圖3 滲透水壓加載過程曲線Fig.3 Curves of the applied water pressure

2 試驗結(jié)果與分析

試驗采用陽江破碎花崗巖標(biāo)準(zhǔn)巖樣，每個試樣開展了不少于10級圍壓下的不同水力梯度的滲流試驗。試驗后的3個破碎巖樣見圖2(b)。由圖2(b)可見，試樣破碎程度較高，縱向裂隙發(fā)育，部分貫通試樣，同時存在少量橫向裂紋。本文給出了試樣1～3巖石試樣在不同圍壓條件下兩端水力梯度▽P與滲流流量Q之間的試驗曲線，見圖4。由圖4可見，▽P～Q關(guān)系曲線皆呈下凸?fàn)?，表明隨著水力梯度的增大，滲流流量的增長幅度逐漸減小。 換言之，隨著流量的增大，增加同等的流量所需消耗的水頭也越來越大。另外，對于同一試樣，圍壓增大可能引起巖樣壓縮，裂隙壓密，過水通道減小，導(dǎo)水性降低，因此對于相同的流量，圍壓越大，所對應(yīng)的水力梯度也越大，故在▽P～Q關(guān)系曲線圖中表現(xiàn)為圍壓越大，曲線越偏離橫軸，偏離程度與圍壓大小的變化基本保持一致。但也存在例外，例如試樣3，在圍壓σ3達(dá)到25.0 MPa后，繼續(xù)增大σ3對▽P～Q曲線的偏離作用略顯微弱，因為此時裂隙的壓密程度接近最大，巖樣透水性基本不再變化。

圖4 試樣在不同圍壓下的▽P～Q關(guān)系曲線Fig.4 Variation of pressure gradient ▽P with flow rate Q under various confining stresses

圖4所示的▽P～Q非線性關(guān)系主要由水流慣性阻力作用的逐漸增強(qiáng)產(chǎn)生[14]，可采用Forchheimer(1901)提出的二項式方程來描述其滲流規(guī)律：

-▽P=AQ+BQ2

(2)

式中：▽P為壓力梯度，Pa/m；Q為滲流流量，m3/s；A為線性項系數(shù)，(Pa·s)/m4；B為非線性項系數(shù)，(Pa·s2)/m7；Ah為過流面積，m2；μ為流體的動力黏度，Pa·s；ρ為密度，kg/m3；k為過流介質(zhì)的固有滲透率，m2；β為非達(dá)西系數(shù)，m-1，又被稱為紊流因子、非達(dá)西滲流慣性系數(shù)，是研究滲流非線性的重要參數(shù)[15-17]。

由式(1)可知，F(xiàn)orchheimer方程。在形式上僅僅是對達(dá)西方程增加了一個二次項(慣性項)，意義上則表示非達(dá)西滲流的能量損失不再僅由黏滯阻力控制，而是由黏滯阻力和慣性阻力2者共同決定[18]，換言之，滲流非達(dá)西效應(yīng)取決于慣性阻力作用的相對大小。

圖5 Forchheimer方程系數(shù)A和B與圍壓的關(guān)系曲線Fig.5 Variation of coefficients A and B with confining stress

采用Forchheimer非達(dá)西滲流方程對▽P～Q數(shù)據(jù)進(jìn)行擬合，得到各曲線的擬合系數(shù)見表1。由表1可知，F(xiàn)orchheimer方程對非線性滲流的擬合效果很好，相關(guān)系數(shù)R2值均介于0.983 6和0.999 9之間，接近于1。為了更為直觀地分析，圖5給出了擬合系數(shù)A、B與圍壓的關(guān)系曲線。由圖5可知，線性項系數(shù)A隨圍壓嚴(yán)格單調(diào)遞增，而非線性項系數(shù)B則隨著圍壓的增大呈總體上升趨勢。由式(2)可知，系數(shù)A反映的是破碎巖體的過流能力，圍壓增大，巖樣內(nèi)部的裂隙被壓密，過流能力減弱，導(dǎo)致系數(shù)A增大；系數(shù)B表征流動的慣性效應(yīng)，破碎巖樣的裂隙分布復(fù)雜，形態(tài)各異，在圍壓增大的過程中，縱向裂隙閉合，橫向裂隙擴(kuò)展，2者對滲流的非線性具有相反的作用機(jī)制，且作用強(qiáng)度具有波動性，但總體上縱向裂隙作用強(qiáng)于橫向裂隙作用，因而系數(shù)B在局部存在起伏波動，但總體上增大。

由Forchheimer方程的式(1)和式(2)可知，當(dāng)忽略滲流流體的物理特性變化時，線性項系數(shù)A取決于巖石固有滲透率k，且與k成反比，非線性項系數(shù)B則由非達(dá)西系數(shù)β確定，并與之成正比。因此圖5所示的系數(shù)A、B的變化關(guān)系在一定程度上反映了非達(dá)西系數(shù)β與固有滲透率k之間存在的相關(guān)關(guān)系；另一方面，在裂隙和多孔介質(zhì)的滲流研究中，非達(dá)西系數(shù)β與固有滲透率k之間廣泛采用冪函數(shù)關(guān)系描述[16,17,19,20]?；诒疚牡臐B流試驗結(jié)果，計算出3個破碎花崗巖巖樣在各圍壓下的固有滲透率k與非達(dá)西系數(shù)β，并將2者關(guān)系曲線繪制于圖6中，采用冪函數(shù)關(guān)系式進(jìn)行擬合，相關(guān)系數(shù)R2=0.916 7，擬合程度較高，表達(dá)式如下：

β=6.962×10-13k-1.709

(3)

3 非達(dá)西滲流評估

如前所述，雷諾數(shù)Re和福希海默數(shù)F0是劃分達(dá)西流和非達(dá)西流常用的指標(biāo)，但是確定臨界值卻十分困難，因此本文在試驗結(jié)果分析得到的非達(dá)西系數(shù)與固有滲透率關(guān)系式的基礎(chǔ)上，建立另一類簡單直觀的非達(dá)西滲流評估方法。

表1 Forchheimer方程擬合系數(shù)Tab.1 The best-fitted coefficients of Forchheimer Equation

圖6 非達(dá)西系數(shù)β與固有滲透率k關(guān)系曲線Fig.6 Variation of non-Darcy coefficient β with intrinsic permeability k

3.1 水力梯度比ΦP

假設(shè)破碎巖體的滲流流量Q已知，則由式(1)計算得到的非達(dá)西水力梯度|-▽P|F與由達(dá)西定律計算所得的達(dá)西水力梯度|-▽P|D2者的比值ΦP滿足[14]：

(4)

上式表明非達(dá)西與達(dá)西水力梯度的比值ΦP由福希海默數(shù)F0決定，因此參數(shù)ΦP也具有表征滲流非達(dá)西程度的意義。將式(2)代入式(4)，并考慮式(3)所得的β～k冪函數(shù)關(guān)系，則有：

(5)

可知，ΦP是固有滲透率k和滲流流速v的函數(shù)。

根據(jù)式(4)和式(5)，圖7(a)和圖7(b)分別做出了ΦP～F0關(guān)系曲線和不同ΦP值下的v～k關(guān)系曲線。由圖7可知，ΦP大于1，表明滲流慣性效應(yīng)的存在將引起額外的壓力損失，因此產(chǎn)生與達(dá)西流相等的流量需要更大的水力梯度。在評估滲流的非線性時，只要確定了水力梯度比ΦP的值，就可從圖7(b)找到對應(yīng)的臨界v～k曲線，該曲線上部區(qū)域的滲流為非達(dá)西流，下部區(qū)域滲流的非線性則忽略不計，可采用達(dá)西定律進(jìn)行滲流分析。雖然忽略非達(dá)西效應(yīng)的滲流計算結(jié)果必然存在偏差，但只要在容許范圍內(nèi)則影響不大。假設(shè)某工程要求水力梯度的容許相對誤差為10%，即對應(yīng)臨界水力梯度比ΦP=1.11(該值對應(yīng)Zeng等[11]提出的達(dá)西流與非達(dá)西流的臨界值F0=0.11)，圖7(b)中ΦP=1.11對應(yīng)v～k曲線即是相應(yīng)的臨界曲線。

圖7 水力梯度比值ΦP評估圖Fig.7 Evaluation of pressure gradient ratio ΦP

3.2 流量比ΦQ

假設(shè)破碎巖體的水力梯度|-▽P|已知，則由式(2)解得的非達(dá)西流量QF為[14]：

(6)

定義流量比ΦQ為非達(dá)西流量QF與由達(dá)西定律計算得到的達(dá)西流量QD的比值。為求得ΦQ，聯(lián)立式(2)和達(dá)西方程，則：

AQD=AQF+BQ2F

(7)

故：

(8)

聯(lián)立式(2)、式(3)和式(8)，得：

(9)

顯然，ΦQ是固有滲透率k和水力梯度|-▽P|的函數(shù)。

同樣地，根據(jù)式(8)和式(9)分別得出了ΦQ～F0關(guān)系曲線和不同ΦQ值下的|-▽P|～k關(guān)系曲線，見圖8。由圖8(a)可知，流量比ΦQ的值為0～1，因為在相同的水力梯度作用下，非達(dá)西流動的慣性項不可忽略，致使非達(dá)西流量小于達(dá)西流量。在評估滲流的非線性時，只要流量比ΦQ的值被確定后，就可以從圖8(b)找到對應(yīng)的臨界|-▽P|～k曲線，該曲線下部區(qū)域的滲流可忽略非線性，采用達(dá)西定律進(jìn)行滲流分析，其他區(qū)域的滲流作為非達(dá)西流處理。若流量的容許相對誤差定為20%，則臨界流量比ΦQ=0.83，圖8(b)中ΦQ=0.83的|-▽P|～k曲線即是相應(yīng)的臨界曲線。

圖8 流量比值ΦQ評估圖Fig.8 Evaluation of flow rate ratio ΦQ

與雷諾數(shù)Re和福希海默數(shù)F0相比，水力梯度比ΦP和流量比ΦQ2個評估參數(shù)的物理意義更為清晰直觀，通過限制非達(dá)西水力梯度或者非達(dá)西流量相對于達(dá)西流情況的偏離程度來判定滲流的非線性，對工程應(yīng)用具有重要的意義。例如油氣開發(fā)更關(guān)注開采量，因此可選擇表示流量偏離的參數(shù)ΦQ來考慮滲流非線性，其臨界值大小由工程實踐中可容許的流量偏差決定；當(dāng)分析與巖體滲透壓力有關(guān)的穩(wěn)定性問題時，選擇水力梯度比ΦP來判定滲流非線性更為合適。但值得指出的是，式(5)和式(9)是在式(3)的基礎(chǔ)上推導(dǎo)得到的，冪函數(shù)關(guān)系式(3)的擬合系數(shù)值與巖體介質(zhì)的滲流特性有關(guān)，因此對于不同的透水介質(zhì)，式(5)和式(9)中各系數(shù)值也會有所差異。

4 結(jié) 語

本文對3個陽江花崗巖破碎后，開展了不同圍壓條件下的高水力梯度滲流試驗，并對之進(jìn)行了分析討論，在試驗結(jié)果的基礎(chǔ)上提出了一種簡單直觀的非達(dá)西滲流評估方法。具體得到如下結(jié)論。

(1)試驗結(jié)果表明，破碎巖體中的滲流在高水力梯度下表現(xiàn)出顯著的非線性，采用Forchheimer方程擬合▽P～Q關(guān)系的相關(guān)系數(shù)為0.983 6～0.999 9。

(2)隨著圍壓的上升，與巖樣過水能力有關(guān)的線性項系數(shù)A嚴(yán)格單調(diào)遞增，反映滲流慣性阻力作用的非線性系數(shù)B呈總體增大趨勢；非達(dá)西系數(shù)β與固有滲透率k具有相關(guān)性較好的冪函數(shù)關(guān)系。

(3)評估滲流的非線性時，采用水力梯度比ΦP和流量比ΦQ更加簡單直觀，且物理意義清晰，可根據(jù)工程具體情況，通過限制非達(dá)西流與達(dá)西流情況下2者的水力梯度比或者流量比來確定滲流非線性的臨界條件，決定是否采用達(dá)西定律進(jìn)行滲流分析。

□

[1] 黃先伍, 唐 平, 繆協(xié)興, 等. 破碎砂巖滲透特性與孔隙率關(guān)系的試驗研究[J]. 巖土力學(xué), 2005,26(9):1 385-1 388.

[2] 李順才, 繆協(xié)興, 陳占清, 等. 承壓破碎巖石非 Darcy 滲流的滲透特性試驗研究[J]. 工程力學(xué), 2008,25(4):85-92.

[3] Thauvin F, Mohanty K K. Network modeling of non-Darcy flow through porous media[J]. Transport in Porous Media, 1998,31(1):19-37.

[4] 胡大偉, 周 輝, 謝守益, 等. 峰后大理巖非線性滲流特征及機(jī)制研究[J]. 巖石力學(xué)與工程學(xué)報, 2009,28(3):451-458.

[5] Zhou J Q, Hu S H, Fang S, et al. Nonlinear flow behavior at low Reynolds numbers through rough-walled fractures subjected to normal compressive loading[J]. International Journal of Rock Mechanics and Mining Sciences, 2015,80:202-218.

[6] 孟如真, 胡少華, 陳益峰, 等. 高滲壓條件下基于非達(dá)西流的裂隙巖體滲透特性研究[J]. 巖石力學(xué)與工程學(xué)報, 2014,33(9):1 756-1 764.

[7] 李 健, 黃冠華, 文 章, 等. 兩種不同粒徑石英砂中非達(dá)西流動的實驗研究[J]. 水利學(xué)報, 2008,39(6):726-732.

[8] Zimmerman R W, Al-yaarubi A, Pain C C, et al. Non-linear regimes of fluid flow in rock fractures[J]. International Journal of Rock Mechanics and Mining Sciences, 2004,41:163-169.

[9] Javadi M, Sharifzadeh M, Shahriar K, et al. Critical Reynolds number for nonlinear flow through rough-walled fractures: The role of shear processes[J]. Water Resources Research, 2014,50(2):1 789-1 804.

[10] 周創(chuàng)兵, 陳益峰, 姜清輝, 等. 復(fù)雜巖體多場廣義耦合分析導(dǎo)論[M]. 北京: 中國水利水電出版社, 2008.

[11] Zeng Z, Grigg R. A criterion for non-Darcy flow in porous media[J]. Transport in Porous Media, 2006,63(1):57-69.

[12] Macici P, Mesini E, Viola R. Laboratory measurements of non-Darcy flow coefficients in natural and artificial unconsolidated porous media[J]. Journal of Petroleum Science and Engineering, 2011,77(3):365-374.

[13] Liu R, Li B, Jiang Y. Critical hydraulic gradient for nonlinear flow through rock fracture networks: the roles of aperture, surface roughness, and number of intersections[J]. Advances in Water Resources, 2016, 88:53-65.

[14] Chen Y F, Zhou J Q, Hu S H, et al. Evaluation of Forchheimer equation coefficients for non-Darcy flow in deformable rough-walled fractures[J]. Journal of Hydrology, 2015,529:993-1 006.

[15] Firdaouss M, Guermond J L, Quéré L P. Nonlinear corrections to Darcy's law at low Reynolds numbers[J]. Journal of Fluid Mechanics, 1997,343:331-350.

[16] 許 凱, 雷學(xué)文, 孟慶山, 等. 非達(dá)西滲流慣性系數(shù)研究[J]. 巖石力學(xué)與工程學(xué)報, 2012,31(1):164-170.

[17] Geertsma J. Estimating the coefficient of inertial resistance in fluid flow through porous media[J]. Society of Petroleum Engineers Journal, 1974,14(5):445-450.

[18] Zhou J Q, Hu S H, Chen Y F, et al. The friction factor in the forchheimer equation for rock fractures[J]. Rock Mechanics and Rock Engineering, 2016,49(8):3 055-3 068.

[19] Li D, Engler T W. Literature review on correlations of the non-Darcy coefficient[C]∥SPE Permian Basin Oil and Gas Recovery Conference. Society of Petroleum Engineers, 2001.

[20] 王 媛, 秦 峰, 夏志皓, 等. 深埋隧洞涌水預(yù)測非達(dá)西流模型及數(shù)值模擬[J]. 巖石力學(xué)與工程學(xué)報, 2012,31(9):1 862-1 868.