1.2 單裂隙中水流及溶質(zhì)運移問題的研究進(jìn)展

1.2.1 單裂隙中水流問題研究進(jìn)展

1.2.1.1 單裂隙水流理論模型的發(fā)展

張開度遠(yuǎn)遠(yuǎn)小于其長寬尺寸的斷裂稱為裂隙。關(guān)于裂隙中水流問題的研究已經(jīng)開展了幾十年［21^，22］，這些研究始終集中在四個基本方面［23］：①建立裂隙概念模型；②建立模型解析解或數(shù)值解的方法；③對裂隙水力特征進(jìn)行描述；④運用隨機(jī)方法描述裂隙隙寬及水文地質(zhì)參數(shù)。早期研究單裂隙時，通常把其簡化成平直、光滑且隙寬處處相等，通常忽略裂隙基質(zhì)本身的滲透性，因此在整個裂隙面上的滲透性是各向同性。等溫條件下不可壓縮牛頓流體Navier-Stokes方程如下

式中 ρ——流體密度；

μ——液體動力黏滯系數(shù)；

u_i——流體在i方向上的速度分量；

P——流體壓力；

g_i——方向i上的重力加速度分量。

對于穩(wěn)定流，式（1.1）中左邊項為0，于是得到

式(1.2)也被稱為Stokes方程，在平行板模型中，該方程的解如下

式中 b——裂隙寬度。

從式(1.3)中可以看出速度剖面呈拋物線形態(tài)，如圖1.2所示。通常被稱為泊肅葉流(Poiseuille's flow)。在寬度上對其積分可得

對其整理即可得到著名的裂隙水流立方定律

式中 q——單寬滲流量；

J——水力坡降；

b——裂隙寬度；

g——重力加速度；

v——水流的運動黏滯系數(shù)。

運用立方定律求得的裂隙內(nèi)部流場為層流且水流速度處處相等，而在紊流階段或微裂隙中運用立方定律求得的結(jié)果往往不盡如人意［24］。

理想裂隙在自然界是不存在的，天然裂隙面均為粗糙裂隙，其隙寬是沿程變化的。通常上下裂隙面會有一定程度的接觸，而且其有效隙寬往往取決于隙面所受法向應(yīng)力的大小。天然裂隙中水流在流動過程中遇到許多阻礙，即使對于裂隙中的純粹泊肅葉流，其與裂隙面接觸的占裂隙體積10%的區(qū)域只傳遞約5%的流量，Djik和Berkowitz［25］運用核磁共振技術(shù)，獲得飽和粗糙裂隙中水流圖像，其中水流速度剖面呈類似拋物線，但不完全對稱，這種速度分布對其中水流的影響十分明顯，因此人們開始懷疑立方定律是否一直有效。Hakami和Larsson［26］對瑞典?sp?地區(qū)采集的天然花崗巖裂隙進(jìn)行了水力實驗和隙寬測量，結(jié)果表明測量得到的平均隙寬是通過立方定律得到的水力隙寬的約1.4倍。隨著研究的深入，人們逐漸認(rèn)識到把裂隙簡單地簡化為平直光滑的平板模型還具有很多其他方面的局限性。例如，裂隙中局部水流密度和溫度的改變都將對流場產(chǎn)生影響，而立方定律只考慮了裂隙的隙寬、尺寸和水的黏滯系數(shù)(動力黏滯系數(shù)或運動黏滯系數(shù))，雖然流體的黏滯系數(shù)和溫度、密度密切相關(guān)，但是由溫度/密度的不均導(dǎo)致的熱對流/密度流卻被忽略。除了溫度、密度之外，裂隙壓力也對其中的流場產(chǎn)生一定的影響。Witherspoon和Wang通過室內(nèi)實驗證明，當(dāng)平板裂隙壓力超過10MPa時，用立方定律得出的結(jié)果和實際情況相差甚遠(yuǎn)［27］。這說明只有在一定的壓力范圍內(nèi)才能運用立方定律進(jìn)行簡化運算。

考慮到立方定律的局限性，眾多學(xué)者對其進(jìn)行了改進(jìn)，該研究主要是為了揭示和定量分析傳統(tǒng)立方定律在實際運用中的偏離情況。而改進(jìn)的關(guān)鍵問題則在于隙寬的確定。一般來說，隙寬主要有三種定義［28］，分別是平均隙寬<b>、機(jī)械隙寬b_m和等效水力隙寬b_h。平均隙寬<b>是指裂隙隙寬函數(shù)b(x,y)的平均值；機(jī)械隙寬b_m為裂隙面發(fā)生的最大閉合變形量；等效水力隙寬b_h是為了應(yīng)用立方定律于實際裂隙而提出的概念，即是將實驗所得裂隙滲流量代入立方定律反求得到的裂隙寬度。對于光滑平行板裂隙，這三種隙寬值是相等的;而對于實際粗糙裂隙，它們通常是不等的。

Lomize通過一系列實驗研究了粗糙裂隙對水流的影響，在此基礎(chǔ)上引進(jìn)了裂隙粗糙度的概念［1］。Snow將此粗糙度概念運用于天然裂隙的水流模擬中，并對傳統(tǒng)立方定律進(jìn)行了相應(yīng)的評價［29］。Lomize和Louis在大量實驗基礎(chǔ)上提出的立方定律修正公式，見表1.1和表1.2。

表1.1和表1.2中為紊流時裂隙滲透系數(shù)，V_f為裂隙水流平均速度，?為裂隙絕對糙率。由于天然裂隙中粗糙顆粒分布不均且凸起高度差異較大，因此在實際情況中仍無法確定?的值。紊流時，水頭損失與流速呈非線性關(guān)系，可用V^m=-表示，式中m為紊流時的非線性指數(shù)，其變化范圍為1～2[2]。光滑裂隙中紊流公式為q^m==，由此可得

通過室內(nèi)實驗可以得到lgJ和lgq之間的關(guān)系直線，其斜率即為m。

Barton等首次提出將巖石節(jié)理粗糙系數(shù)JRC運用于裂隙水力隙寬的求解，提出水力隙寬、機(jī)械隙寬和JRC之間的關(guān)系式［31］

于是立方定律修正為

上式中的難點在于JRC的確定，雖然很多學(xué)者提出一系列JRC的測量方法［32^-35］，然而由于其尺度效應(yīng)［36］及粗糙度的各向異性等條件限制，使得測量準(zhǔn)確度仍有待提高。

Iwai通過大量單裂隙水流實驗發(fā)現(xiàn)裂隙的面積接觸率ω(裂隙面接觸面積與總面積之比)與水流規(guī)律存在一定的聯(lián)系［37］。據(jù)此得出的理想裂隙模型中，沿水流法向方向，裂隙開度連續(xù)變化；而在沿水流方向，裂隙開度不變。進(jìn)而Walsh［38］和周創(chuàng)兵等［39］分別推導(dǎo)出如下公式

式中 κ——經(jīng)驗系數(shù)，Walsh建議κ=1，周創(chuàng)兵等建議κ=0。

Amadei等則提出如下的修正公式［40］

式中 σ_b——隙寬均方差。

Nolte對等石英二長巖芯試樣進(jìn)行滲流實驗,3個試樣的成果整理發(fā)現(xiàn)，隙寬指數(shù)n遠(yuǎn)遠(yuǎn)大于3，隨隙寬增大分別為7.6、8.3和9.8［41］。張有天等采用計算機(jī)生成人工裂隙和有限元數(shù)值計算方法，經(jīng)分析認(rèn)為，單寬流量q與b_m也不是3 次方關(guān)系，n均大于3，并隨裂隙粗糙程度的增加而增大［42］。耿克勤根據(jù)人工、天然光滑和粗糙裂隙的實驗結(jié)果分析得到，對于小開度裂隙層流而言，1.7≤n≤3.0，裂隙面幾何形態(tài)越光滑n值越大;對于中開度過渡狀態(tài)，0.8≤n≤1.4;對于大開度裂隙，0.3≤n≤0.48［43］。在此基礎(chǔ)上，許光祥等進(jìn)行了歸納總結(jié)，提出了次立方定律和超立方定律的概念，并通過試驗進(jìn)行了驗證［44］。

Neuzil和Tracy模擬了由一系列不同隙寬的平行板模型組成的裂隙，這些平行板模型的隙寬分布符合對數(shù)正態(tài)分布。模擬結(jié)果顯示，此裂隙中的水流主要通過隙寬較大的路徑，裂隙中水流存在優(yōu)勢流［45］。

Tsang等利用電阻代替水流發(fā)現(xiàn)裂隙的曲折程度對其中的水流有著較為明顯的影響，而且隙寬越小，曲折度的影響越明顯。和普通的光滑平行板模型相比，裂隙曲折度和隙面粗糙度對水流產(chǎn)生的阻滯影響可以將水流速度降低3個數(shù)量級以上［46］。Brown詳細(xì)研究了裂隙粗糙度對水流的影響［47］，他利用分形方法生成的粗糙裂隙，模擬了其中的水流狀態(tài)。他認(rèn)為，在隙寬較小的情況下，裂隙的曲折度對其中的水流有著十分明顯的影響，這也進(jìn)一步證明了Tsang的結(jié)論。

Rasmason和Neretnieks歷時3年在瑞典一個廢棄的Stripa礦里進(jìn)行了現(xiàn)場實驗［48］，通常稱這個實驗為Stripa-3D實驗，它揭示了裂隙巖體中水流和運移的主要特征。Heath［48］和Bourke［49］也在康沃爾做了現(xiàn)場實驗，他們的研究發(fā)現(xiàn)，水流很大程度上在裂隙巖體中被隔離的管道流里面流動，并沒有大面積的水流發(fā)生，水流僅僅在5%~20%的裂隙面內(nèi)流動。在加拿大安大略的白堊紀(jì)石灰?guī)r和美國伊利諾伊州東北部的志留紀(jì)白云巖中，Novakowski［50］、Raven［51］和Shapiro［52］分別進(jìn)行了現(xiàn)場示蹤劑實驗，發(fā)現(xiàn)通過裂隙隙寬估計的水力傳導(dǎo)系數(shù)和溶質(zhì)的遷移規(guī)律與現(xiàn)場觀測的不一致，如果裂隙隙寬的變化很小，則與實際觀測吻合得較好。因此，上述的現(xiàn)場實驗提供給我們這樣一個事實：裂隙面的平行板模型假定并不適用，水流和溶質(zhì)運移只在裂隙形成的管道或溝槽中發(fā)生。

Pyrak［53］等在不同壓力條件下，將熔化的金屬代替水流注入裂隙中，冷卻后打開，發(fā)現(xiàn)液體在其中的流動路徑是曲折的，沿一定的路徑流動，具有優(yōu)勢流和溝槽流的特征。為了有效地模擬裂隙中的溝槽流，Y.W.Tsang和C.F.Tsang［46^，54］提出了一種新的概念模型。他們運用對數(shù)正態(tài)分布函數(shù)并在平均、差異和空間相關(guān)長度基礎(chǔ)上生成隙寬統(tǒng)計分布，進(jìn)行了水流和溶質(zhì)運移數(shù)值模擬，預(yù)測了其中的穿透曲線，預(yù)測的結(jié)果和Moreno等［55］提供的數(shù)據(jù)能很好地吻合。

1.2.1.2 單裂隙物理模型實驗

由于現(xiàn)場水文地質(zhì)條件以及巖體的復(fù)雜性，開展現(xiàn)場實驗工作需要大量的人力和物力，測試方法和測試裝置的不同，所得到的結(jié)果也不一定完全一致［56］。室內(nèi)實驗相對來說花費較少而且有較精確的測試方法，因此大部分實驗通過室內(nèi)的小尺度模型進(jìn)行。

研究者最初是利用兩塊光滑平板，通過控制兩塊平行板之間的距離來模擬不同隙寬的裂隙，這類裝置較為簡單，使用的平行板通常為有機(jī)玻璃板或玻璃板，國內(nèi)外眾多學(xué)者都做過類似的實驗［1^{，24，30］}，圖1.3為速寶玉等制作的光滑裂隙水流實驗裝置。在此基礎(chǔ)上進(jìn)一步發(fā)展，人們開始在平行板上粘貼粗糙顆粒來模仿粗糙裂隙，大多數(shù)此類實驗都是利用標(biāo)準(zhǔn)粒徑的砂粒來充當(dāng)粗糙顆粒［57^-59］，并通過粘貼不同粒徑的砂粒和控制平行板之間的距離來獲得不同的粗糙裂隙。也有研究者利用機(jī)械加工后的鋼板來模擬粗糙裂隙。許光祥等［44］為了驗證其理論，加工了包括光滑裂隙在內(nèi)的5種不同形態(tài)的裂隙進(jìn)行水流實驗，如圖1.4所示。由于粘貼的砂粒粒徑或機(jī)械加工出的粗糙顆粒都比較均勻，雖然獲得了粗糙的裂隙，但和實際情況仍相距甚遠(yuǎn)。于是研究者開始利用混凝土來獲得仿真裂隙［56^，60-62］并進(jìn)行水流及溶質(zhì)運移實驗(圖1.5)。這在一定程度上讓研究者獲得更為接近實際的實驗數(shù)據(jù)，但是混凝土和天然巖體在材料特性上存在較大差別，如透水性、表面可濕性等，都對其內(nèi)部水流及裂隙本身形態(tài)產(chǎn)生影響，于是一些學(xué)者開始利用巖芯中的天然裂隙［63］，或者利用天然巖樣制作人工裂隙，如圖1.6所示。也有直接在野外進(jìn)行較大尺度的單裂隙水流實驗［64］。在過去的幾十年中，研究人員逐漸開始運用諸如裂隙粗糙面剖面測量、低熔點金屬注入和樹脂鑄造等技術(shù)來獲得更加接近真實裂隙的模型來進(jìn)行相關(guān)的模擬［46^{，47，53，65］}，這些模型考慮了裂隙隙寬變化及其對水流影響的實際情況，對研究變隙寬裂隙中的水流運動提供了極大的幫助。

1.2.1.3 裂隙水流問題的數(shù)值模擬

對于某些復(fù)雜問題，如在尋求裂隙的滲透性質(zhì)與影響溶質(zhì)運移的控制因素的關(guān)系時，現(xiàn)場實驗和室內(nèi)實驗的測試手段都不完美，這就需要求助于數(shù)值模擬，對于這些復(fù)雜問題，數(shù)值模擬是一種比較方便且精度較高的方法。一些室內(nèi)實驗無法模擬的現(xiàn)象通過數(shù)值模擬得到了較好的解決［56］。

自20世紀(jì)60年代末期以來，關(guān)于裂隙中水流等相關(guān)問題的數(shù)值模擬已經(jīng)發(fā)展了40多年［66］。據(jù)統(tǒng)計，截至1994年，已經(jīng)出現(xiàn)至少30余種求解此類相關(guān)問題的專門軟件［67］。Streltsova-Adams［68］在求解混合含水層井流問題時描述了幾種求解裂隙水流的解析法。Elsworth［69］提出了求解具有特定幾何形狀層流或紊流問題的解析方法，由于對問題中水流的形狀有特定要求，因此其應(yīng)用受到了限制。Amadei和Illangasekare［70］運用積分變換得到矩形裂隙中流體勢能場的連續(xù)表達(dá)式。這個方法可以在裂隙隙寬和粗糙度各向異性的情況下使用，同時也可以為進(jìn)一步研究立方定律提供依據(jù)。由于積分變換的運用，不需要再對裂隙進(jìn)行離散，但是裂隙的幾何形態(tài)不能過于復(fù)雜，運用僅限于裂隙形態(tài)較為簡單的情況。

在求解裂隙中水流及傳質(zhì)問題時，通常用到差分和積分法來求解平衡方程。對于空間求導(dǎo)問題，積分法往往獲得比有限差分法(FDM)更為廣泛的應(yīng)用，原因之一就是積分法能更好地適應(yīng)不規(guī)則幾何區(qū)域。常用的積分法包括有限單元法(FEM)［70^{，71，72］}和邊界元法(BEM)。Elsworth［69^，73］曾運用BEM-FEM混合方法來模擬裂隙水流問題，同時積分有限差分法(IFDM)也獲得了廣泛的運用［74^-76］。積分型有限差分法通常使用Crank-Nicholson法逼近來處理時間導(dǎo)數(shù)。為了在高Péclet數(shù)條件下獲得穩(wěn)定解，Sudicky和McLaren［77］引進(jìn)了Laplace變換。

20世紀(jì)80年代，在格子氣自動機(jī)(Lattice Gas Automata)理論的基礎(chǔ)上發(fā)展出了一種新型流場模擬方法——格子波茨曼方法(Lattice Boltzmann Method)［78^-84］。與傳統(tǒng)的模擬方法相比，該方法具有規(guī)則簡單、復(fù)雜邊界易處理、能適應(yīng)大規(guī)模并行計算等優(yōu)點，并迅速成為模擬復(fù)雜流場的新工具［85^-92］。

早在20世紀(jì)六七十年代，研究人員就已經(jīng)開始利用諸如時間序列法、譜分析法和蒙特卡羅法來模擬裂隙隙寬概率分布，這些方法的共同點就是通過對隨機(jī)參數(shù)的控制得到所需要的介質(zhì)特點，如非均質(zhì)性和各向異性。20世紀(jì)90年代，隨著分形技術(shù)的充分發(fā)展，以及對裂隙面分形特點認(rèn)識的逐漸成熟，人們開始采用分形技術(shù)對裂隙面的粗糙度和裂隙側(cè)面進(jìn)行分析和模擬［93^-95］。眾多的研究表明隙寬滿足對數(shù)正態(tài)分布或者高斯分布［26^，96］，同時，隙寬分布的方差較大，說明了模擬隙寬非均質(zhì)性的重要性［56］。