蘇軍偉， 王 樂， 呂 高， 劉媛箐

(1.西安交通大學(xué) 人居環(huán)境與建筑工程學(xué)院， 陜西 西安710049； 2. 西安石油大學(xué) 機(jī)械工程學(xué)院， 陜西 西安710065)

1 研究背景

地下水是地球數(shù)量豐度分布廣泛的淡水資源，對(duì)于人類生產(chǎn)生活均有著重要意義[1]?，F(xiàn)階段，研究者多采用有限元方法、有限差分和觀測(cè)統(tǒng)計(jì)等方法[2-6]，從宏觀角度解析地下水流動(dòng)及人類活動(dòng)對(duì)地下水質(zhì)量及污染物運(yùn)移的影響。鄭燦政等[2]采用MODFLOW軟件基于有限差分方法，研究了地下結(jié)構(gòu)物對(duì)地下水流場(chǎng)的影響，發(fā)現(xiàn)地下結(jié)構(gòu)物不僅影響地下水流場(chǎng)的穩(wěn)定時(shí)間，而且使得迎水面水位上升，背水面水位下降。吳樂等[7]采用有限差分方法對(duì)北京西山地區(qū)含水層系統(tǒng)進(jìn)行數(shù)值模擬，分析不同開采條件下含水層系統(tǒng)響應(yīng)特征，發(fā)現(xiàn)維持現(xiàn)狀開采至2030年，巖溶水水位下降22 m，第四系承壓水水位下降約28 m。代鋒剛等[8]以山西長(zhǎng)治盆地潞安礦區(qū)為例, 通過野外調(diào)查、相似材料模擬實(shí)驗(yàn)和數(shù)值模擬手段, 分析了群礦開采驅(qū)動(dòng)下含水層結(jié)構(gòu)破壞對(duì)地下水流場(chǎng)的影響。雖然現(xiàn)有研究已經(jīng)取得了豐碩成果，但目前的宏觀地下水?dāng)?shù)值模擬僅重視模型研究，不重視具體水文地質(zhì)條件的研究，影響了宏觀數(shù)值模擬結(jié)果的準(zhǔn)確性[9]。此外，傳統(tǒng)的達(dá)西尺度分析無法準(zhǔn)確描述復(fù)雜結(jié)構(gòu)中的滲流規(guī)律以及微納尺度流動(dòng)中的非達(dá)西效應(yīng)[10]。

從微觀角度上來說，地下水流動(dòng)于巖石孔隙之中，通過真實(shí)巖石孔隙中地下水的流動(dòng)分析能夠?yàn)楹暧^地下水流動(dòng)提供進(jìn)一步的認(rèn)識(shí)分析。隨著圖像學(xué)和計(jì)算流體力學(xué)的發(fā)展，通過對(duì)巖石進(jìn)行CT掃描獲得物理模型[11-12]，并進(jìn)一步對(duì)該模型進(jìn)行網(wǎng)格劃分得到真實(shí)的巖石孔隙計(jì)算域，結(jié)合計(jì)算流體力學(xué)，使得地下水流動(dòng)及污染物遷移在真實(shí)巖石孔隙內(nèi)的研究成為可能。近年來，基于真實(shí)巖心的數(shù)值模擬技術(shù)發(fā)展迅速，其中，哈佛大學(xué)Datta等[13]在微尺度下研究了三維多孔結(jié)構(gòu)內(nèi)的流體流動(dòng)行為，發(fā)現(xiàn)復(fù)雜孔隙結(jié)構(gòu)內(nèi)流體流動(dòng)并不是隨機(jī)的，且沿流動(dòng)方向的流動(dòng)速度呈指數(shù)分布。斯坦福大學(xué)的Roman等[14]采用微型PIV技術(shù)測(cè)量了孔隙內(nèi)的單相流動(dòng)速度，且采用數(shù)值模擬驗(yàn)證了實(shí)驗(yàn)結(jié)果，并對(duì)不相融兩相流動(dòng)進(jìn)行了數(shù)值模擬研究。Yang等[15]采用多種數(shù)值方法模擬了流體在巖石孔隙內(nèi)的流動(dòng)行為，從宏觀變量和微觀變量?jī)蓚€(gè)方面對(duì)不同數(shù)值方法的模擬結(jié)果進(jìn)行比較。清華大學(xué)王沫然等[10]對(duì)地下深層巖石微納米孔隙內(nèi)氣體滲流進(jìn)行了跨尺度混合模擬，為頁(yè)巖氣勘探開發(fā)提供了相關(guān)理論支持。Su等[16-17]基于歐拉-拉格朗日架構(gòu)下，開發(fā)了顆粒追蹤的精確算法，對(duì)真實(shí)巖石中的驅(qū)油過程進(jìn)行了多相數(shù)值模擬，揭示了真實(shí)巖石中多相流動(dòng)機(jī)制及流動(dòng)行為。以上研究表明了基于真實(shí)巖心的流體流動(dòng)行為研究逐漸展開，且由單相不斷擴(kuò)展至多相，研究結(jié)論也進(jìn)一步解釋了宏觀尺度研究中觀察到的現(xiàn)象[18-20]。

本文從微觀角度出發(fā)，通過對(duì)真實(shí)巖心進(jìn)行掃描，建立真實(shí)微觀巖石孔隙結(jié)構(gòu)模型，并基于此對(duì)地下水在巖石孔隙內(nèi)的流動(dòng)進(jìn)行數(shù)值模擬，探討不同流量及不同出入口位置對(duì)地下水流動(dòng)的動(dòng)力學(xué)行為影響，以期完善地下水流動(dòng)的微尺度數(shù)值模擬方法，為地下水的流動(dòng)行為提供微尺度解釋。

2 數(shù)學(xué)模型

根據(jù)連續(xù)介質(zhì)假設(shè)，只要流動(dòng)特征尺度遠(yuǎn)大于流體微元尺度，可采用N-S方程描述流動(dòng)，從目前實(shí)驗(yàn)測(cè)量看出，連續(xù)假設(shè)可以應(yīng)用于10 nm以上的流體流動(dòng)[21]，因此對(duì)于微米級(jí)真實(shí)孔隙尺度，連續(xù)介質(zhì)假設(shè)仍然適用。

質(zhì)量方程：

(1)

動(dòng)量方程：

(2)

式中：ρ為密度;u為速度;p為壓強(qiáng);ν為黏度。壓力和速度耦合采用PISO算法[22]。上述數(shù)學(xué)模型通過有限體積法進(jìn)行離散。

3 物理模型及計(jì)算方法

本研究中通過巖石掃描技術(shù)結(jié)合數(shù)字圖像處理技術(shù)得到真實(shí)巖石孔隙結(jié)構(gòu)的物理模型，二維真實(shí)孔隙尺度巖心的物理模型如圖1所示。該模型長(zhǎng)、寬分別為0.010 3和0.008 53 m，且模型內(nèi)分布著多條微米級(jí)復(fù)雜、曲折、相互貫通的孔隙通道。地下水從左側(cè)區(qū)域向右側(cè)區(qū)域流動(dòng)，為更好地說明入口和出口位置，圖1(a)中標(biāo)注的數(shù)字1、2、3、4、5、6、7、8分別對(duì)應(yīng)入口1、入口2、入口3、入口4、出口5、出口6、出口7及出口8。

采用非結(jié)構(gòu)化網(wǎng)格對(duì)該物理模型進(jìn)行網(wǎng)格劃分，圖1(b)中給出了圖1(a)中方框區(qū)域?qū)?yīng)網(wǎng)格的細(xì)部劃分，顯示了在微米級(jí)孔隙通道內(nèi)，網(wǎng)格分布均勻密布。經(jīng)網(wǎng)格獨(dú)立性檢驗(yàn)，本文選擇總量為221 742的網(wǎng)格進(jìn)行計(jì)算。

對(duì)不同的注入速度及出入口位置的算例設(shè)置如表1所示，算例a、b和c分別探討相同入口位置不同入口速度下孔隙尺度地下水流動(dòng)行為；算例b、d及e分別探討不同入口位置相同流量下孔隙尺度內(nèi)地下水流動(dòng)行為；算例f為多入口及出口的算例便于與以上算例的對(duì)比分析，其入口流量與算例b一致。水的密度為1 000 kg/m3，運(yùn)動(dòng)黏度設(shè)置為0.000 001 m2/s。

采用OpenFOAM開源流體力學(xué)軟件中的icoFoam求解器，設(shè)定時(shí)間步長(zhǎng)為10-6秒，最大庫(kù)朗數(shù)為0.1，各個(gè)物理量殘差設(shè)置為10-6。

4 結(jié)果分析與討論

圖2為整個(gè)巖石孔隙通道內(nèi)的速度分布。由圖2可以看出，不同區(qū)域的地下水流動(dòng)速度差異較大，尤其是在四個(gè)端角附近，流動(dòng)速度受邊界阻滯的影響，流動(dòng)較為緩慢，且在相鄰?fù)ǖ纼?nèi)受孔道大小及兩側(cè)壓差的影響，流動(dòng)速度也存在較大變化；此外，在單向封閉孔隙通道——盲道(通道僅存在一個(gè)出入口)內(nèi)地下水流動(dòng)近乎停滯，以上現(xiàn)象說明真實(shí)巖石孔隙通道內(nèi)的地下水流動(dòng)具有復(fù)雜性和非均勻性。對(duì)比圖2(a)、2(b)及2(c)，發(fā)現(xiàn)隨著入口速度的增大，巖石中部孔隙內(nèi)的速度不斷增大，而左側(cè)上部區(qū)域速度變化差異較小，這主要是由于左側(cè)上部存在著多處不連通孔隙且受邊界條件影響，共同形成了流動(dòng)死區(qū)。為了進(jìn)一步說明不連通孔隙對(duì)地下水流動(dòng)的影響，對(duì)圖2(a)中方框區(qū)域進(jìn)行了放大分析，詳見圖4敘述。此外，對(duì)比圖2(b)、2(d)及2(e)，發(fā)現(xiàn)隨著入口位置的改變，地下水流動(dòng)的速度分布發(fā)生明顯變化。當(dāng)入口位置位于中下部時(shí)(圖2(d))，對(duì)比圖2(b)，發(fā)現(xiàn)巖石孔隙下部區(qū)域的流動(dòng)情況得到改善，下部水流速度明顯增高，僅在巖石上部左側(cè)和上部右側(cè)孔隙存在著明顯的速度低值區(qū)；當(dāng)入口位置位于中上部時(shí)(圖2(e))，對(duì)比圖2(b)，發(fā)現(xiàn)左側(cè)上部巖石孔隙的地下水流動(dòng)速度明顯增加，且這一變化也體現(xiàn)在中部及中下部區(qū)域，僅下部水流速度沒有明顯變化。這意味著在流量相同的情況下，改變?nèi)肟谖恢媚軌蚋纳莆⒂^巖石內(nèi)的地下水流動(dòng)情況，且無論是改變?nèi)肟谒俣然蛘{(diào)整入口位置，單向封閉孔隙通道內(nèi)的流動(dòng)情況都無法得到改善。需要說明的是，本算例中出口位置并沒有變化。當(dāng)出口位置上移或下移時(shí)，受出口位置的影響，在中部及右側(cè)區(qū)域的地下水流動(dòng)方向?qū)㈦S之變化，并增大出口周邊區(qū)域的地下水流動(dòng)速度。圖2(f)給出了多入口、多出口下的巖石孔隙結(jié)構(gòu)內(nèi)水流速度分布，對(duì)比同一流量下算例b，發(fā)現(xiàn)整個(gè)巖石孔隙區(qū)域中部的水流速度明顯較低，而巖石右上部區(qū)域速度有著一定的升高。

圖1 二維真實(shí)巖石孔隙結(jié)構(gòu)模型及網(wǎng)格劃分

參數(shù)算例a算例b算例c算例d算例e算例f入口速度/(m·s-1)0.0010.0020.0030.003010.007360.0009入口位置入口1入口1入口1入口2入口3入口1、2、3、4出口位置出口6出口6出口6出口6出口6出口5、6、7、8

圖3進(jìn)一步給出了孔隙結(jié)構(gòu)內(nèi)的y方向速度分布，通過設(shè)置兩個(gè)速度區(qū)間區(qū)分了y正方向和負(fù)方向的水流流動(dòng)行為，間接反映了孔隙結(jié)構(gòu)內(nèi)水的流動(dòng)方向。發(fā)現(xiàn)所有算例中，在圖3(a)方框區(qū)域內(nèi)的y方向速度分布極為接近，差異主要體現(xiàn)在巖石孔隙結(jié)構(gòu)的左側(cè)區(qū)域。其中隨著流量的增大，對(duì)比圖3(a)～3(c)，發(fā)現(xiàn)y方向速度分布沒有明顯差異，這也意味著入口速度的不斷增大并沒有明顯改變水流在孔隙結(jié)構(gòu)內(nèi)的y方向流動(dòng)。而入口位置的改變對(duì)巖石孔隙結(jié)構(gòu)內(nèi)速度分布影響明顯，對(duì)比圖3(b)、3(d)及3(e)，發(fā)現(xiàn)左側(cè)區(qū)域受入口位置改變，速度分布存在明顯差異，意味著在這些區(qū)域水的流動(dòng)方向變化劇烈，但受影響的區(qū)域明顯有限，主要集中于入口周邊區(qū)域。如圖3(f)所示，多入口及多出口下，y方向速度分布與圖3(a)～3(c)中的速度分布較為近似，這可能是由于入口1和出口6為孔隙結(jié)構(gòu)面積最大的兩個(gè)出入口，導(dǎo)致了水流主要經(jīng)由這兩個(gè)口流入和流出，這也顯示了孔隙結(jié)構(gòu)內(nèi)水流方向往往受面積大的入口和出口影響。

圖4所示為各算例部分孔隙結(jié)構(gòu)速度分布(圖2(a)中方框區(qū)域放大部分)，發(fā)現(xiàn)無論入口速度大小、入口位置及出入口數(shù)量的改變，位于圖中左側(cè)中部區(qū)域的單向封閉孔隙通道內(nèi)(圖4(c)方框2區(qū)域)的地下水流動(dòng)速度始終極低；通道內(nèi)徑越小，相對(duì)的地下水流速度越高；通道內(nèi)的地下水流動(dòng)速度高值位于通道中心，受壁面剪切影響貼近壁面處地下水流動(dòng)速度較小。此外，在圖4(a)、4(b)以及4(c)中部區(qū)域(圖4(e)方框3區(qū)域)，隨著入口地下水注入速度增大，孔隙通道內(nèi)流動(dòng)速度并沒有進(jìn)一步增大，其可能的原因是連接該孔隙的上下兩側(cè)通道內(nèi)地下水流動(dòng)速度較為接近，沒有形成明顯的壓差，且入口速度的增大也沒有改變這一現(xiàn)象，導(dǎo)致該區(qū)域速度較低；上部區(qū)域(圖4(c)方框1區(qū)域)由于孔隙通道較細(xì)，速度變化較為明顯，速度持續(xù)增大。

注：圖(a)～(f)分別對(duì)應(yīng)算例a～f

注：圖(a)～(f)分別對(duì)應(yīng)算例a～f

對(duì)比相同流量下不同入口位置時(shí)的算例，如圖4(b)、4(d)及4(e)所示，發(fā)現(xiàn)地下水流動(dòng)速度分布區(qū)域差異主要體現(xiàn)在圖4(e)中的方框3和方框4區(qū)域。當(dāng)入口位置由中部移至下部時(shí)，對(duì)比圖4(b)及4(d)，發(fā)現(xiàn)方框3區(qū)域的地下水流動(dòng)速度增高，同時(shí)方框4區(qū)域的地下水流動(dòng)速度也有一定程度升高；當(dāng)入口位置由中部移至上部時(shí)，對(duì)比圖4(b)及4(e)，發(fā)現(xiàn)方框3區(qū)域的地下水流動(dòng)速度明顯增大，而方框4區(qū)域的速度明顯減小。如圖4(f)所示，當(dāng)存在多入口及多出口時(shí)，在方框3和方框4區(qū)域，對(duì)比圖4(b)，孔隙通道中部水流速度有著少許增大，其他區(qū)域速度分布沒有明顯差異。

上述討論分析了孔隙通道內(nèi)地下水流動(dòng)的速度分布，但仍然無法精確描述地下水流動(dòng)方向等行為，為此圖5給出了部分孔隙通道內(nèi)的地下水流線。對(duì)比圖5(a)、5(b)及5(c)，發(fā)現(xiàn)隨著入口速度的增大，孔隙結(jié)構(gòu)內(nèi)的流線沒有發(fā)生明顯的差異，意味著水流方向一致；旋渦主要分布于單向封閉孔隙內(nèi)(對(duì)應(yīng)圖5(d)方框4及方框3區(qū)域)及上下連通通道內(nèi)(對(duì)應(yīng)圖5(d)方框2區(qū)域)，這導(dǎo)致了地下水不易流入也難以流出。需要說明的是上下連通的孔隙通道內(nèi)部旋渦的形成可能是由于流速較為接近，上下沒有明顯的壓差時(shí)，該孔隙通道內(nèi)形成了兩個(gè)旋渦，影響了地下水向上方和下方通道的流動(dòng)。

當(dāng)入口位置發(fā)生改變時(shí)，對(duì)比圖5(b)、5(d)及5(e)，發(fā)現(xiàn)地下水流動(dòng)方向發(fā)生了較大變化。尤其是當(dāng)入口位置由中部下移，對(duì)比圖5(b)與5(d)，在方框1及方框2區(qū)域均發(fā)生了明顯的流線改變，其中圖5(d)中方框1區(qū)域地下水由底部向上部流動(dòng)，而在方框2區(qū)域中原有旋渦消失，地下水由下向上流動(dòng)；當(dāng)入口位置由中部上移(如圖5(e)所示)，同樣方框2區(qū)域內(nèi)的旋渦消失，此外在方框3區(qū)域內(nèi)地下水流線發(fā)生轉(zhuǎn)向，流動(dòng)方向?yàn)橛傻撞肯蛏喜苛鲃?dòng)。值得一提的是入口位置的改變并不影響單向封閉通道內(nèi)的水流流動(dòng)形態(tài)。當(dāng)存在多入口及多出口時(shí)，如圖5(f)所示，其流線與圖5(d)接近，而與圖5(b)中的方框2和方框1區(qū)域有著明顯差異?？傊?，無論改變?nèi)肟谒俣?、入口位置或出入口?shù)量，巖石孔隙結(jié)構(gòu)內(nèi)均存在著流動(dòng)相似區(qū)域，而流動(dòng)方向差異較為明顯的區(qū)域集中在縱向孔隙通道內(nèi)。

圖6給出了中軸線位置處的x方向速度，其中中軸線位置為x=0.004 45 m。從圖6中發(fā)現(xiàn)，速度經(jīng)歷著增大、減小、中斷的變化過程，這也意味著在每個(gè)孔道內(nèi)，貼近孔道兩側(cè)部位受孔壁的影響速度較低，孔道中心速度較高，這與圖4中的現(xiàn)象一致。此外，速度的中斷也意味著真實(shí)巖石內(nèi)孔道間是不連續(xù)的，這也符合巖石的真實(shí)地質(zhì)構(gòu)造，而現(xiàn)有研究采用多孔介質(zhì)進(jìn)行假設(shè)時(shí)會(huì)忽略這一現(xiàn)象。當(dāng)入口速度不斷增大時(shí)，對(duì)比圖6(a)、6(b)及6(c)，發(fā)現(xiàn)各個(gè)孔道內(nèi)速度值均不斷增大且不同孔隙通道內(nèi)的速度峰值排布序列沒有發(fā)生變化，僅從流動(dòng)速度上來說，在靠近y=0.000 2 m以及y=0.003 5 m附近的孔道內(nèi)存在兩處地下水流動(dòng)速度的較高值，這意味著這兩個(gè)孔道為地下水流動(dòng)的優(yōu)勢(shì)通道。當(dāng)改變?nèi)肟谖恢脮r(shí)，對(duì)比圖6(b)、6(d)及6(e)，發(fā)現(xiàn)不同孔隙通道內(nèi)的速度峰值排布序列發(fā)生變化，尤其是圖6(d)中原有的優(yōu)勢(shì)通道僅在y=0.000 2 m處，而在圖6(e)中y=0.000 2 m、y=0.003 5 m、及y=0.005 2 m處均存在著較高的速度值，這意味著改變?nèi)肟谖恢每赡軙?huì)帶來地下水流動(dòng)優(yōu)勢(shì)通道的改變。當(dāng)存在多個(gè)入口和出口時(shí)，如圖6(f)所示，在y=0.000 2 m孔道內(nèi)地下水流動(dòng)速度最大；對(duì)比圖6(b)發(fā)現(xiàn)，在y=0.003 5 m處速度減小，但不同孔隙通道內(nèi)的速度峰值排布序列較為一致。

注：圖(a)～(f)分別對(duì)應(yīng)算例a～f

注：圖(a)～(f)分別對(duì)應(yīng)算例a～f

圖7進(jìn)一步給出了不同速度權(quán)重區(qū)間下的占比，其中單個(gè)網(wǎng)格內(nèi)速度權(quán)重的計(jì)算公式為：速度權(quán)重=網(wǎng)格點(diǎn)速度/孔隙孔道內(nèi)最大速度，根據(jù)速度權(quán)重的大小歸屬至0～1中的10個(gè)等分區(qū)間內(nèi)得到不同速度權(quán)重區(qū)間的占比。從圖7中可以發(fā)現(xiàn)，算例a、b、c的數(shù)據(jù)基本一致，巖石孔隙內(nèi)速度權(quán)重區(qū)間為0～0.1的占比最高，均在0.75以上，意味著地下水流動(dòng)低速區(qū)占據(jù)了整個(gè)巖石的絕大部分；速度權(quán)重區(qū)間為0.1～0.2的區(qū)域同樣占比較高，但明顯低于速度權(quán)重0～0.1，上述現(xiàn)象意味著低速區(qū)在地下水孔隙通道流動(dòng)速度中占據(jù)主導(dǎo)。對(duì)比算例a、算例b及算例c，發(fā)現(xiàn)改變?nèi)肟谒俣鹊拇笮?，在不同速度?quán)重區(qū)間下，未發(fā)生明顯變化，均為低速區(qū)占據(jù)主要權(quán)重區(qū)間。對(duì)比算例b、算例d及算例e，發(fā)現(xiàn)改變?nèi)肟谖恢媚軌蛞欢ǔ潭雀淖兯俣葯?quán)重區(qū)間占比，尤其是算例e中，速度權(quán)重區(qū)間為0～0.1的占比相比于其他兩個(gè)算例增高了約0.06，意味著低速區(qū)的范圍進(jìn)一步擴(kuò)大，這對(duì)于地下水的污染物運(yùn)移不利。而多入口多出口下，在速度權(quán)重區(qū)間為0～0.1的占比最低為0.75，明顯低于算例b的占比，這意味著增加入口和出口數(shù)量能夠改善孔隙結(jié)構(gòu)內(nèi)的水流整體流動(dòng)情況。

圖7 各算例巖石孔隙水流不同速度權(quán)重區(qū)間占比

5 結(jié) 論

本研究針對(duì)巖石孔隙內(nèi)地下水流動(dòng)行為進(jìn)行微觀數(shù)值模擬研究，探討了不同出入口位置及入口速度對(duì)地下水流動(dòng)行為的影響，主要結(jié)論如下：

(1)真實(shí)巖石孔隙內(nèi)地下水流動(dòng)過程非常復(fù)雜，地下水流速呈現(xiàn)非均一性。隨著入口速度的不斷增大，巖石孔隙內(nèi)不同的速度權(quán)重基本沒有差別，但改變著巖石孔隙內(nèi)的速度大小。其中，單孔隙內(nèi)成為地下水流動(dòng)死區(qū)，不受外界入口位置及入口速度大小的影響。

(2)入口位置的變化會(huì)改變地下水在巖石孔隙中流動(dòng)的優(yōu)勢(shì)通道數(shù)量和位置，且能夠改變孔隙通道內(nèi)的流動(dòng)方向，尤其是能夠影響孔隙結(jié)構(gòu)內(nèi)低速區(qū)的分布權(quán)重。

(3)多入口多出口下，巖石中部速度降低，地下水流動(dòng)方向受面積較大的出口和入口影響，不同孔隙通道內(nèi)速度排布序列沒有明顯變化，速度權(quán)重0～0.1的占比最小。