余 成，雷文武

(重慶交通大學(xué) 河海學(xué)院，重慶 400074)

?

粗糙單裂隙溶質(zhì)運移優(yōu)先通道的模擬研究

余 成，雷文武

(重慶交通大學(xué) 河海學(xué)院，重慶 400074)

為考察巖石粗糙裂隙內(nèi)不均勻分布的裂隙空間對溶質(zhì)運移過程的影響，采用基于LCL方法的數(shù)值模型對不同裂隙寬度下的流場和溶質(zhì)運移過程進行模擬。結(jié)果表明：隨著裂隙寬度減小，流場逐漸從均勻的平行板流過渡到強烈不均勻的溝槽流，流線逐漸集中于部分優(yōu)先溝槽通道；溶質(zhì)運移鋒面的推進也從基本呈直線狀平行推進，逐漸開始不規(guī)則化，出現(xiàn)指流現(xiàn)象；隨著寬度繼續(xù)減小，溝槽通道進一步發(fā)展，形成貫穿整個流場的溶質(zhì)運移優(yōu)先通道；優(yōu)先運移通道的溶質(zhì)運移速度遠(yuǎn)大于其他部位，出現(xiàn)某些部位溶質(zhì)運移受困而滯留的現(xiàn)象。因此，在裂隙寬度較小時，粗糙裂隙內(nèi)的優(yōu)先溝槽通道對裂隙流動和溶質(zhì)運移過程有著非常重要的影響。

水利工程；粗糙裂隙；溶質(zhì)運移；數(shù)值模擬；優(yōu)先溝槽通道

0 引 言

粗糙裂隙內(nèi)流體流動與溶質(zhì)運移過程在巖土工程、油氣工程及環(huán)境工程領(lǐng)域受到極大的重視。裂隙流及溶質(zhì)運移過程可能直接影響地下污染物的釋放和運移過程[1]。在核廢料地下填埋的處理工作中，地下巖體裂隙的水力滲透特征及其對溶質(zhì)運移過程的影響也受到關(guān)注[2]。裂隙在地下巖體中常常以裂隙網(wǎng)絡(luò)的形式存在，但截至目前，人們對流體和溶質(zhì)在單條粗糙裂隙內(nèi)運動過程的了解還不夠深入，對于單裂隙內(nèi)水流和溶質(zhì)運移過程的研究仍然是水文學(xué)和環(huán)境科學(xué)的一大熱點[3-5]。

其中的難點主要在于裂隙空間和流場分布的不規(guī)則性。裂隙表面具有典型的分形特征[6]，其隨機性決定了裂隙表面的粗糙性和復(fù)雜性。在裂隙寬度較大時，表面粗糙性對流動和溶質(zhì)運移過程的影響較小，但裂隙寬度減小時，表面粗糙性對流場的影響則顯著增大[7]。另外由于巖石碎裂時常常伴隨著相互錯動，使得裂隙內(nèi)部空間分布強烈不均。隨著裂隙寬度減小，各處寬度的相對差異逐漸增大。寬度較小的部位水力滲透性較差，不利于流體流動，而寬度較大的部位則可能形成局部優(yōu)先溝槽通道[8]。優(yōu)先溝槽通道一方面使裂隙呈現(xiàn)局部透水特征，另一方面也影響著流動及溶質(zhì)運移過程的規(guī)則性。裂隙面的復(fù)雜性決定了溝槽通道的復(fù)雜性，以及流動和溶質(zhì)運移過程的復(fù)雜性。至今人們對于粗糙裂隙內(nèi)流動和溶質(zhì)運移過程依然缺乏更加深入的研究。

筆者運用LCL模型對粗糙裂隙內(nèi)的流動和溶質(zhì)運移過程進行模擬研究。流動模型顯示，隨著裂隙寬度減小，流場逐漸從均勻的平行板流過渡到強烈不均勻的溝槽流，流線逐漸集中于部分優(yōu)先溝槽通道；溶質(zhì)運移過程受到溝槽通道較強的影響；裂隙寬度較大時，濃度運移鋒面基本呈直線狀平行推進，寬度減小時，濃度推進鋒面逐漸不規(guī)則化，出現(xiàn)指流現(xiàn)象，寬度繼續(xù)減小，溝槽通道進一步發(fā)展，指流現(xiàn)象也進一步演化，形成貫穿整個流場的溶質(zhì)運移優(yōu)先通道；優(yōu)先運移通道的溶質(zhì)運移速度遠(yuǎn)大于部分局部的溶質(zhì)運移速度，出現(xiàn)某些局部溶質(zhì)運移受困而遲滯的現(xiàn)象。

1 粗糙裂隙面的隨機生成

裂隙表面是具有分形特征的粗糙表面。實踐中用于生成分形地形表面的建模方法有分?jǐn)?shù)維布朗運動法、傅立葉濾波法、小波變換法、中點隨機位移法等。筆者采用比較流行的菱形-方形細(xì)分法(Diamond-Square)。其原理可以簡述為：首先對于給定的正方形區(qū)域，假設(shè)4個角點A,B,C,D賦值高程分別為hA，hB，hC，hD，隨后的過程分為Diamond步和 Square步。

1.1 Diamond步

給正方形的中心點賦高程值為，如式(1)：

(1)

式中：δ0為擾動隨機變量，通?？珊唵蔚夭捎镁鶆蚍植糢(ξ,-ξ)(表示在該方形區(qū)域內(nèi)高程擾動的上限值)。

1.2Square步

(2)

式中：δ0同前。

如此方形ABCD可被分割為4個小方形，對每一個小方形可使用以上兩步進行遞歸，但在第i層遞歸中，擾動變量應(yīng)修改為，如式(3)：

(3)

式中：k為分形表面粗糙度，k∈[0，1]，越小分形表面越粗糙。

2 裂隙流動模型

粗糙裂隙內(nèi)的流動可以使用復(fù)雜邊界下的Navier-Stokes方程描述，但直接Navier-Stokes流動模擬具有對網(wǎng)格劃分精度要求高、非線性微分方程求解困難等缺點。實踐中人們更多地使用局部立方定律(LCL)。LCL認(rèn)為裂隙的上下表面在局部范圍內(nèi)為平行板，并具有局部滲透系數(shù)K：

(4)

式中：b為局部裂隙寬度；v為流體運動黏度系數(shù)。

LCL可將定常三維裂隙流簡化為二維平面滲流模型，在低雷諾數(shù)和通常幾何條件下，LCL方法與直接Navier-Stokes模擬之間的差異一般在5%～10%范圍內(nèi)[9]，因此被人廣泛接受。

利用隨機生成的粗糙表面作為裂隙空間的上下表面，通過上下表面間的水平和垂向相對運動可生成不均勻分布的裂隙空間，如圖1。裂隙空間的分布與表面粗糙性有關(guān)外，還與裂隙表面間的相對位置有關(guān)。相對的水平、垂向以及傾斜、旋轉(zhuǎn)等都可能使裂隙空間發(fā)生不規(guī)則變化。筆者主要考察上下表面在固定水平相對位移(x，y方向)時，垂直(z方向)位移變化過程中的裂隙空間變化，以及相應(yīng)地流動和溶質(zhì)運移過程所受到的影響。

模型采用了的有限差分網(wǎng)格，每個網(wǎng)格寬度為0.1mm(圖1)。模型 y 方向的兩側(cè)邊界被設(shè)定為定水頭邊界，水頭差為常數(shù)0.1mm；模型x方向的兩側(cè)邊界被設(shè)定為隔水邊界。

圖1 文中模型所采用的粗糙裂隙寬度的空間分布

3 溶質(zhì)運移模型

溶質(zhì)運移模型依賴于流動模型，但溶質(zhì)運移除了有對流運動外，還有明顯的擴散和彌散運動。由于速度邊界層的影響，裂隙流體速度在寬度方向上呈不均勻分布，導(dǎo)致較為顯著的彌散效應(yīng)。對于平行板流動，溶質(zhì)彌散系數(shù)D由流體性質(zhì)和平行板間距決定，如式(5)[10-11]：

(5)

式中：D*為分子擴散系數(shù)；b為裂隙寬度；U為裂隙流的平均流速。

對于裂隙流，理論上D在空間上呈不均勻分布，但由于較難統(tǒng)計裂隙各處流量，且在b較小時彌散作用對溶質(zhì)運移過程的影響比對流速度較小，因此筆者近似地認(rèn)為D在空間上為均勻分布，其大小利用裂隙流平均流速和裂隙平均寬度計算。

為了考察定常流場狀態(tài)下溶質(zhì)的運移過程，模型中溶質(zhì)被認(rèn)為從0時刻開始以定濃度的方式釋放在模型的定水頭邊界。通過考察溶質(zhì)分布狀態(tài)以及溶質(zhì)運移前鋒面的推進過程，研究優(yōu)先溝槽通道對溶質(zhì)運移過程的影響。裂隙流動模型采用MODFLOW程序?qū)崿F(xiàn)，溶質(zhì)運移模型采用MT3DMS程序?qū)崿F(xiàn)。

4 結(jié)果與討論

4.1 導(dǎo)水系數(shù)變化

則本例中，當(dāng)α=-0.119時，依然可以很好地滿足立方定律(圖2)，說明在h=0時，由于水平錯動形成的溝槽通道的影響，可以認(rèn)為是一個寬度為0.119mm的附加裂隙寬度。

圖2 粗糙裂隙導(dǎo)水系數(shù)與測量平均隙寬和校正平均隙寬之間的關(guān)系

4.2 裂隙流場變化

裂隙上下表面垂向運動時，各處間的滲透系數(shù)值變化差異巨大。通常較寬部位的值變化程度較小，而較窄部位的值受影響較大。相應(yīng)地，裂隙流場形態(tài)、等水頭線及流線分布也必然受到影響。如圖3，裂隙寬度較大時，流場狀態(tài)比較接近于平行板流，流線與等水頭線主要呈平行直線分布；但在寬度減小時，滲透系數(shù)的不均勻性最終使得流線和水頭分布也越加曲折。

圖3 不同裂隙寬度下的等水頭線

隨著裂隙各處值相對差異變大，水頭梯度逐漸集中于裂隙通道比較狹窄的部位。表明隨著溝槽通道的逐漸形成，這些溝槽通道尤其是其中比較狹窄的部位，將逐漸成為裂隙水力滲透性的重要影響因素。對應(yīng)地，流線分布特征(圖4)也反映出在優(yōu)先溝槽通道逐漸形成的過程中流動逐漸向溝槽通道不斷集中的過程。

圖4 不同裂隙寬度下的流線分布

4.3 溶質(zhì)運移鋒面變化

流動過程的變化必然引起溶質(zhì)運移過程的變化。優(yōu)先溝槽通道也很大程度地影響了裂隙流的溶質(zhì)運移過程。裂隙寬度變化以及溝槽通道形成的過程中，裂隙流動速度的變化不僅整體上影響了溶質(zhì)運移的速度，同時，裂隙流不斷集中的過程也很大程度上影響了濃度的分布狀態(tài)。

圖5顯示了不同裂隙寬度下，在特征時刻的典型濃度分布狀態(tài)。從濃度鋒面分布看，在裂隙寬度很大時，裂隙粗糙性的影響很小，濃度鋒面形狀主要呈平行推進，與平行板間的溶質(zhì)運移過程極為類似。

圖5 不同裂隙寬度下特征時刻的溶質(zhì)濃度分布

隨著裂隙寬度減小，流動逐漸向優(yōu)先溝槽通道集中，此時濃度運移鋒面逐漸出現(xiàn)不均勻變化，鋒面上開始出現(xiàn)局部的指流現(xiàn)象，但濃度的整體推進過程依然保持相對一致。此時優(yōu)先溝槽通道尚處于形成階段，對溶質(zhì)運移過程只具有部分的影響，還沒有起到?jīng)Q定性作用。

隨著裂隙寬度繼續(xù)減小，流線進一步向溝槽通道集中，溶質(zhì)運移過程也進一步發(fā)生較大變化。濃度鋒面的不均勻程度進一步加劇，指流現(xiàn)象進一步演化，發(fā)展出非常明顯的優(yōu)先運移通道。由于流場的不均勻性，各處間的溶質(zhì)運移速度差異明顯。此時優(yōu)先溝槽通道對溶質(zhì)運移過程逐漸產(chǎn)生決定性的影響作用，濃度鋒面推進一致性差，溶質(zhì)主要沿著優(yōu)先溝槽通道運移。

在寬度減小過程中，部分局部滲透系數(shù)減小迅速，其溶質(zhì)運移速度遠(yuǎn)小于優(yōu)先溝槽通道的溶質(zhì)運移速度，因此呈現(xiàn)出在某些部位溶質(zhì)可能陷入局部的“困境”，在裂隙內(nèi)運移相當(dāng)緩慢從而長期滯留，從而使裂隙出口處的濃度穿透曲線呈現(xiàn)非常顯著的拖尾現(xiàn)象。這與多孔介質(zhì)中由于機械彌散或者吸附過程的影響而造成的濃度曲線拖尾現(xiàn)象類似。在研究上可以與多孔介質(zhì)中的溶質(zhì)運移過程進行類比。

5 結(jié) 語

筆者利用LCL方法對粗糙裂隙間的裂隙流動以及溶質(zhì)運移過程進行了模擬研究。裂隙上下表面通過截取菱形-方形細(xì)分法所獲得的粗糙表面生成，利用上下表面間的相對運動生成不規(guī)則分布的裂隙空間。裂隙流動模型顯示，在裂隙寬度逐漸減小過程中，裂隙流逐漸從平行板流過渡到具有強烈不均勻特征的溝槽流，其中逐漸形成的優(yōu)先溝槽通道對裂隙流場有著較大的影響，流線和水頭梯度都逐漸集中于這些優(yōu)先溝槽通道。

溶質(zhì)運移過程也較大程度地受到優(yōu)先溝槽通道的影響。在裂隙寬度較大時，濃度運移的前鋒面基本呈直線狀的平行推進，在寬度減小時，受到優(yōu)先溝槽通道的影響，濃度鋒面逐漸不規(guī)則化，出現(xiàn)局部指流現(xiàn)象。在裂隙寬度繼續(xù)減小時，優(yōu)先溝槽通道繼續(xù)進一步發(fā)展，成為影響溶質(zhì)運移過程的主要因素。指流也進一步發(fā)展和演化，形成貫穿整個流場的溶質(zhì)運移優(yōu)先通道。優(yōu)先運移通道的溶質(zhì)運移速度遠(yuǎn)大于某些局部的溶質(zhì)運移速度，因此可能出現(xiàn)某些局部溶質(zhì)被困而滯留的現(xiàn)象。

[1] 錢家忠,汪家權(quán),葛曉光,等.我國北方型裂隙巖溶水流及污染物運移數(shù)值模擬研究進展[J].水科學(xué)進展,2003,14(4):509-512. Qian Jiazhong,Wang Jiaquan,Ge Xiaoguang,et al.Advances in research for numerical simulation of contaminant transport and flow in North China type fracture-karst media [J].Advances in Water Science,2003,14(4):509-512.

[2] Seetharam S C,Perko J,Jacques D,et al.Influence of fracture networks on radionuclide transport from solidified waste forms [J/OL].Nuclear Engineering and Design,2014,270:162-175 [2014-02-18].http://www.sciencedirect.com/science/article/pii/S002954931400020X.

[3] Wu Yushu,Ye Ming,Sudicky E A.Fracture-flow-enhanced matrix diffusion in solute transport through fractured porous media [J/OL].Transport in Porous Media,2010,81(1):21-34 [2014-02-11].http://link.springer.com/article/10.1007%2Fs11242-009-9383-4.

[4] Qian J Z,Chen Z,Zhan H B,et al.Solute transport in a filled single fracture under non-Darcian flow [J/OL].International Journal ofRock Mechanics and Mining Sciences,2011,48(1):132-140[2014-02-21].http://www.sciencedirect.com/science/article/pii/S1365160910001796.

[5] 鞠楊,張欽剛,楊永明,等.巖體粗糙單裂隙流體滲流機制的實驗研究[J].中國科學(xué):技術(shù)科學(xué),2013,43(10):1144-1154. Ju Yang,Zhang Qingang,Yang Yongming,et al.An experimental investigation on the mechanism of fluid flow through single rough fracture of rock [J].Scientia Sinica Technologica,2013,43(10):1144-1154.

[6] Borodich F M.Fractals and fractal scaling in fracture mechanics [J/OL].International Journal of Fracture,1999,95(1/2/3/4):239-259 [2014-02- 15].http://link.springer.com/article/10.1023%2FA%3A1018660604078.

[7] Witherspoon P A,Wang J S Y,Iwai K,et al.Validity of cubic law for fluid flow in a deformable rock fracture [J/OL].Water Resources Research,1980,16(6):1016-1024 [2014-02-23].http://onlinelibrary.wiley.com/doi/10.1029/WR016i006p01016/abstract.

[8] Xiong Xiangbin,Li Bo,Jiang Yujing,et al.Experimental and numerical study of the geometrical and hydraulic characteristics of a single rock fracture during shear [J/OL].International Journal of Rock Mechanics and Mining Sciences,2011,48(8):1292-1302 [2014-02-25].http://www.sciencedirect.com/science/article/pii/S136516091100147X.

[9] Brush D J,Thomson N R.Fluid flow in synthetic rough-walled fractures:Navier-Stokes,Stokes,and local cubic law simulations [J/OL].Water Resources Research,2003,39(4):1085-1100 [2014-02-01].http:// onlinelibrary.wiley.com/doi/10.1029/2002WR001346/abstract.

[10] Aris R.On the dispersion of a solute in a fluid flowing through a tube [J/OL].Proceedings of the Royal Society:A,1956,235:67-77 [2014-02-18].http://rspa.royalsocietypublishing.org/content/235/1200/ 67.

[11] Taylor G.Dispersion of soluble matter in solvent flowing slowly through a tube [J/OL].Proceedings of the Royal Society:A,1953,219:186-203 [2014-02-10].http://rspa.royalsocietypublishing.org/content/ 219/1137/186.

Numerical Modeling of Preferential Flow and Transport Channel between Rough-Walled Fractures

Yu Cheng, Lei Wenwu

(School of River & Ocean Engineering, Chongqing Jiaotong University, Chongqing 400074, China)

In order to show the impact of non-uniform void space on the flow and transport processes between rough-walled fractures, LCL modeling was employed to simulate the process of the solute transport under the flow with different crack width. The result suggests that, along with the decrease of fracture apertures, the flow field tends to get irregular, and streamlines tend to concentrate at preferential flow channels. Solute concentration front looks parallel when the fracture aperture is still significant, but gets irregular and looks similar to finger flow while fracture aperture decreases. As apertures decrease, preferential channels tend to form preferential transport channels. Solute transport speed at preferential transport channels tends to be much higher, thus some of the solute at low permeable regions seems to be trapped, or retarded within fractures. Therefore it can be concluded that when fracture apertures are small, preferential channels have significant impact on flow and transport processes between rough fractures.

hydraulic engineering; rough-walled fractures; solute transport; numerical modeling; preferential flow channel

10.3969/j.issn.1674-0696.2015.04.17

2015-01-12；

2015-04-01

國家自然科學(xué)基金項目(51409028)；重慶交通大學(xué)科研啟動基金項目(2010310002)

余 成(1984—)，重慶人，講師，博士，主要從事地下滲流及溶質(zhì)運移方面的研究。E-mail：beimingyu@pku.edu.cn。

P 641.2

A

1674-0696(2015)04-091-04