巖石孔隙內地下水流動的微觀數值模擬研究

2019-11-26 06:20:14蘇軍偉劉媛箐

水資源與水工程學報 2019年5期

關鍵詞：區域

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

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

1 研究背景

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

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

本文從微觀角度出發，通過對真實巖心進行掃描，建立真實微觀巖石孔隙結構模型，并基于此對地下水在巖石孔隙內的流動進行數值模擬，探討不同流量及不同出入口位置對地下水流動的動力學行為影響，以期完善地下水流動的微尺度數值模擬方法，為地下水的流動行為提供微尺度解釋。

2 數學模型

根據連續介質假設，只要流動特征尺度遠大于流體微元尺度，可采用N-S方程描述流動，從目前實驗測量看出，連續假設可以應用于10 nm以上的流體流動[21]，因此對于微米級真實孔隙尺度，連續介質假設仍然適用。

質量方程：

(1)

動量方程：

(2)

式中：ρ為密度;u為速度;p為壓強;ν為黏度。壓力和速度耦合采用PISO算法[22]。上述數學模型通過有限體積法進行離散。

3 物理模型及計算方法

本研究中通過巖石掃描技術結合數字圖像處理技術得到真實巖石孔隙結構的物理模型，二維真實孔隙尺度巖心的物理模型如圖1所示。該模型長、寬分別為0.010 3和0.008 53 m，且模型內分布著多條微米級復雜、曲折、相互貫通的孔隙通道。地下水從左側區域向右側區域流動，為更好地說明入口和出口位置，圖1(a)中標注的數字1、2、3、4、5、6、7、8分別對應入口1、入口2、入口3、入口4、出口5、出口6、出口7及出口8。

采用非結構化網格對該物理模型進行網格劃分，圖1(b)中給出了圖1(a)中方框區域對應網格的細部劃分，顯示了在微米級孔隙通道內，網格分布均勻密布。經網格獨立性檢驗，本文選擇總量為221 742的網格進行計算。

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

采用OpenFOAM開源流體力學軟件中的icoFoam求解器，設定時間步長為10-6秒，最大庫朗數為0.1，各個物理量殘差設置為10-6。

4 結果分析與討論

圖2為整個巖石孔隙通道內的速度分布。由圖2可以看出，不同區域的地下水流動速度差異較大，尤其是在四個端角附近，流動速度受邊界阻滯的影響，流動較為緩慢，且在相鄰通道內受孔道大小及兩側壓差的影響，流動速度也存在較大變化；此外，在單向封閉孔隙通道——盲道(通道僅存在一個出入口)內地下水流動近乎停滯，以上現象說明真實巖石孔隙通道內的地下水流動具有復雜性和非均勻性。對比圖2(a)、2(b)及2(c)，發現隨著入口速度的增大，巖石中部孔隙內的速度不斷增大，而左側上部區域速度變化差異較小，這主要是由于左側上部存在著多處不連通孔隙且受邊界條件影響，共同形成了流動死區。為了進一步說明不連通孔隙對地下水流動的影響，對圖2(a)中方框區域進行了放大分析，詳見圖4敘述。此外，對比圖2(b)、2(d)及2(e)，發現隨著入口位置的改變，地下水流動的速度分布發生明顯變化。當入口位置位于中下部時(圖2(d))，對比圖2(b)，發現巖石孔隙下部區域的流動情況得到改善，下部水流速度明顯增高，僅在巖石上部左側和上部右側孔隙存在著明顯的速度低值區；當入口位置位于中上部時(圖2(e))，對比圖2(b)，發現左側上部巖石孔隙的地下水流動速度明顯增加，且這一變化也體現在中部及中下部區域，僅下部水流速度沒有明顯變化。這意味著在流量相同的情況下，改變入口位置能夠改善微觀巖石內的地下水流動情況，且無論是改變入口速度或調整入口位置，單向封閉孔隙通道內的流動情況都無法得到改善。需要說明的是，本算例中出口位置并沒有變化。當出口位置上移或下移時，受出口位置的影響，在中部及右側區域的地下水流動方向將隨之變化，并增大出口周邊區域的地下水流動速度。圖2(f)給出了多入口、多出口下的巖石孔隙結構內水流速度分布，對比同一流量下算例b，發現整個巖石孔隙區域中部的水流速度明顯較低，而巖石右上部區域速度有著一定的升高。

圖1 二維真實巖石孔隙結構模型及網格劃分

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

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

注：圖(a)～(f)分別對應算例a～f

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

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

當入口位置發生改變時，對比圖5(b)、5(d)及5(e)，發現地下水流動方向發生了較大變化。尤其是當入口位置由中部下移，對比圖5(b)與5(d)，在方框1及方框2區域均發生了明顯的流線改變，其中圖5(d)中方框1區域地下水由底部向上部流動，而在方框2區域中原有旋渦消失，地下水由下向上流動；當入口位置由中部上移(如圖5(e)所示)，同樣方框2區域內的旋渦消失，此外在方框3區域內地下水流線發生轉向，流動方向為由底部向上部流動。值得一提的是入口位置的改變并不影響單向封閉通道內的水流流動形態。當存在多入口及多出口時，如圖5(f)所示，其流線與圖5(d)接近，而與圖5(b)中的方框2和方框1區域有著明顯差異。總之，無論改變入口速度、入口位置或出入口數量，巖石孔隙結構內均存在著流動相似區域，而流動方向差異較為明顯的區域集中在縱向孔隙通道內。

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

注：圖(a)～(f)分別對應算例a～f

圖7進一步給出了不同速度權重區間下的占比，其中單個網格內速度權重的計算公式為：速度權重=網格點速度/孔隙孔道內最大速度，根據速度權重的大小歸屬至0～1中的10個等分區間內得到不同速度權重區間的占比。從圖7中可以發現，算例a、b、c的數據基本一致，巖石孔隙內速度權重區間為0～0.1的占比最高，均在0.75以上，意味著地下水流動低速區占據了整個巖石的絕大部分；速度權重區間為0.1～0.2的區域同樣占比較高，但明顯低于速度權重0～0.1，上述現象意味著低速區在地下水孔隙通道流動速度中占據主導。對比算例a、算例b及算例c，發現改變入口速度的大小，在不同速度權重區間下，未發生明顯變化，均為低速區占據主要權重區間。對比算例b、算例d及算例e，發現改變入口位置能夠一定程度改變速度權重區間占比，尤其是算例e中，速度權重區間為0～0.1的占比相比于其他兩個算例增高了約0.06，意味著低速區的范圍進一步擴大，這對于地下水的污染物運移不利。而多入口多出口下，在速度權重區間為0～0.1的占比最低為0.75，明顯低于算例b的占比，這意味著增加入口和出口數量能夠改善孔隙結構內的水流整體流動情況。

圖7 各算例巖石孔隙水流不同速度權重區間占比

5 結論

本研究針對巖石孔隙內地下水流動行為進行微觀數值模擬研究，探討了不同出入口位置及入口速度對地下水流動行為的影響，主要結論如下：

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

(2)入口位置的變化會改變地下水在巖石孔隙中流動的優勢通道數量和位置，且能夠改變孔隙通道內的流動方向，尤其是能夠影響孔隙結構內低速區的分布權重。