張學(xué)羿, 竇 智(河海大學(xué)地球科學(xué)與工程學(xué)院, 江蘇 南京 210098)

地下水是可供人類利用的重要資源,而隨著社會經(jīng)濟持續(xù)發(fā)展和人口快速增長,工業(yè)污水排放、化肥農(nóng)藥施用、垃圾填埋等導(dǎo)致了大量污染物進入地下水,地下水污染現(xiàn)已成為全世界共同關(guān)注的問題之一。持續(xù)泄露的可溶性污染物在地下水流運動的長期作用下可運移數(shù)公里甚至數(shù)百公里遠(yuǎn),嚴(yán)重影響了地下水環(huán)境。研究治理地下水污染,首先要了解此類污染的過程和特征,對其運移過程中水流場和濃度場變化的研究也十分關(guān)鍵[1]。孫家強等[2]以黏土模擬低滲透包氣帶環(huán)境,探究了混溶流體在包氣帶中的遷移規(guī)律；左自波等[3]模擬了土體非飽和帶中污染物在降雨條件下的遷移過程；余期沖等[4]研究了死端孔隙對溶質(zhì)運移的影響。但以上研究中對于孔隙微觀結(jié)構(gòu)影響的探討不多。

黏土多孔介質(zhì)結(jié)構(gòu)十分復(fù)雜,并影響著土體的宏觀特征(孔隙度、毛細(xì)壓力、相對滲透率等),進而影響了溶質(zhì)的運移。對于黏土微觀孔隙結(jié)構(gòu)研究的主要問題在于多孔介質(zhì)的描述問題,即如何采用數(shù)學(xué)的方法對微觀孔隙結(jié)構(gòu)進行定量描述。大多傳統(tǒng)研究中采用只能反映“整體”性能的連續(xù)介質(zhì)模型,故不能反映其微觀特征[5]。對于微觀研究,首先需構(gòu)建多孔介質(zhì)的微觀結(jié)構(gòu),Tacher L等[6]采用一種離散的減少距離的方法生成圓形、橢圓形的兩相粒狀多孔介質(zhì)；Pilotti等[7]采用球體沉降法生成三維多孔介質(zhì)；這兩類方法生成的孔隙結(jié)構(gòu)比較簡單,雖然也能反映一定的特征,但是和實際形態(tài)仍有差距。Dong Hu等[8]、劉向君等[9]、黃家國等[10]等學(xué)者采用掃描成像技術(shù)分析微觀孔隙結(jié)構(gòu),但成本高昂。Wang M等[11]結(jié)合多孔介質(zhì)構(gòu)造的孔隙生長模型和叢生理論,提出了四參數(shù)隨機生長方法(Quartet Structure Generation Set, QSGS)；周瀟等[12]利用該方法探討了滲流速度與孔隙通道大小和土顆粒大小之間的關(guān)系。與傳統(tǒng)的圓形(或橢圓形)顆粒構(gòu)造的黏土微觀結(jié)構(gòu)相比,使用隨機生長四參數(shù)生成法構(gòu)造的黏土微觀結(jié)構(gòu)與真實情況更接近。

黏土多孔介質(zhì)內(nèi)部通道形態(tài)各異,且分布不規(guī)律,在滲流過程中常出現(xiàn)一些沿著優(yōu)先途徑的集中流動,統(tǒng)稱為優(yōu)先流[13]。本文利用隨機生長四參數(shù)法構(gòu)建多孔介質(zhì)模型,通過Stokes Equation (Stokes方程) 與Advection-Diffusion Equation (對流擴散方程,ADE) 的耦合計算,模擬黏土多孔介質(zhì)中可溶性污染物的運移過程,旨在探討?zhàn)ね廖⒂^結(jié)構(gòu)對可溶性污染物運移的影響。

1 控制方程及問題的建立

1.1 黏土微觀孔隙結(jié)構(gòu)的生成

為獲得黏土多孔介質(zhì)的微觀孔隙結(jié)構(gòu),采用了隨機生長四參數(shù)生成法(QSGS)生成多孔介質(zhì)的孔隙輪廓,并用AutoCAD建立可供有限元分析軟件使用的模型。該方法的構(gòu)造過程接近于多孔介質(zhì)的生成過程,與傳統(tǒng)的圓形或橢圓形顆粒構(gòu)成的多孔介質(zhì)相比,其固相顆粒形態(tài)不規(guī)則且大小具有隨機性,生成的多孔介質(zhì)孔隙包含連通孔隙和非連通孔隙。模型的構(gòu)建主要步驟如下：

(1)利用Matlab實現(xiàn)Wang M等[11]提出的隨機生長四參數(shù)生成法的程序編譯,生成孔隙度n為0.7的多孔介質(zhì)輪廓,保存固相顆粒的輪廓坐標(biāo)數(shù)據(jù)：

①確定一定大小的構(gòu)造區(qū)域(本文采用199×199的節(jié)點范圍),初始狀態(tài)該范圍內(nèi)全為孔隙。隨機分布初始固體顆粒,分布概率Pa應(yīng)小于最終多孔介質(zhì)固相所占比例。游歷構(gòu)造區(qū)域內(nèi)的每個節(jié)點并在[0,1]內(nèi)取隨機數(shù),隨機數(shù)小于或等于Pa的節(jié)點設(shè)置為初始固體顆粒。

②在初始固體顆粒位置按不同方向的生長概率Pi生成固相顆粒：取8個生長方向(i=1,2,…,8),其中包括4個主要生長方向(i=1,2,3,4)和4個次生長方向(i=5,6,7,8),生長概率不都相同且P1～4>P5～8,本文中P1～4=0.1、P5～8=0.025,可獲得各向異性的多孔介質(zhì)。各方向相鄰節(jié)點在[0,1]內(nèi)取隨機數(shù), 隨機數(shù)小于或等于Pi的節(jié)點生長為固體顆粒。

③反復(fù)執(zhí)行步驟②,直至孔隙率達(dá)到設(shè)定值n。生成固相顆粒輪廓圖并保存輪廓坐標(biāo)數(shù)據(jù)。

(2)運行AutoCAD中的VBA宏程序,將固相顆粒的輪廓坐標(biāo)數(shù)據(jù)自動繪制成大小為19.9 mm×19.9 mm的輪廓模型,并添加四周邊界后,保存為可供有限元分析軟件COMSOL使用的DXF文件。

(3)將所得的DXF文件導(dǎo)入COMSOL,模型不同區(qū)域設(shè)置為固相材料或液相材料。

1.2 控制方程

微觀結(jié)構(gòu)黏土的可溶性污染物運移問題是水流場和濃度場的兩場耦合問題,因此,分別對水流場和濃度場的控制方程進行闡述。為了更真實地模擬黏土微觀結(jié)構(gòu)中水流場的狀態(tài),采用Stokes方程：

-ρgφ+μu=0

(1)

(2)

式中：ρ——水密度/(kg·m-3)；

g——重力加速度/(m·s-2)；

φ——水頭高度/m；

μ——水的動力黏滯系數(shù)/(kg·m-1·s-1)；

u——水流速度矢量/(m·s-1)。

通過上式,可解出黏土多孔介質(zhì)中的非均勻速度場,將非均勻速度場代入對流-擴散方程：

(3)

式中：u——由式(1)計算得到的二維水流場速度矢量；

c——可溶性污染物的濃度/(mol·L-3)；

Df——分子擴散系數(shù)/(m2·s-1)。

本文假設(shè)Df=1.0×10-9m2/s。需要注意的是,式(3)表示本文考慮的可溶性污染物具有理想特征,即不考慮污染物的化學(xué)反應(yīng)特征,且在顆粒表面無吸附。

1.3 問題的建立和求解

本文旨在利用數(shù)值模擬的方法探討各向異性黏土多孔介質(zhì)對可溶性污染物運移的影響。使用上述DXF格式的圖像文件在COMSOL中建立孔隙度n為0.7的多孔介質(zhì)幾何模型(圖1)。其水流場設(shè)置為：左邊和右邊的進口和出口邊界為定壓力邊界,顆粒邊界和上下邊界為無滑移邊界,水流在5 Pa的壓力下流動。

圖1 孔隙度n=0.7的多孔介質(zhì)模型Fig.1 Porous media model (n=0.7)

COMSOL是一款基于多物理場耦合的有限元計算分析軟件,可以用來求解線性、非線性問題及與時間有關(guān)的穩(wěn)態(tài)、瞬態(tài)問題。本文通過Stokes Equation (Stokes方程) 與Advection-Diffusion Equation (對流擴散方程,ADE) 的耦合,利用有限元的方法對計算區(qū)域進行網(wǎng)格剖分細(xì)化、求解,輸出速度場和隨時間變化的濃度場數(shù)據(jù)及圖形。

2 黏土微觀結(jié)構(gòu)中水流場分析

黏土微觀結(jié)構(gòu)中的水流場是研究可溶性污染物運移的基礎(chǔ)。水流場的控制方程由式(1)和式(2)構(gòu)成。為真實反映多孔介質(zhì)地下水水流狀態(tài),由進口和出口5 Pa壓力差控制的平均水流場速度為0.005 m/s,相應(yīng)的水流場分布如圖2所示。

圖2 水流場分布(單位：m·s-1)Fig.2 Flow field distribution (unit: m·s-1)

圖3顯示了孔隙度n=0.7的黏土微觀孔隙通道中在5 Pa的壓力下水流場的分布情況。其分布十分復(fù)雜,受孔隙連通情況、孔隙形態(tài)、顆粒分布和大小的影響明顯,最快流速達(dá)到0.0478 m/s,而最慢流速接近于0。從圖3選取的微觀孔隙通道的流速等值線圖可以發(fā)現(xiàn),孔隙中的流速基本滿足速度由中心區(qū)域向四周降低的特征。而此通道隙寬較窄且連通性好,共與四個有水流流動的通道相連接,出現(xiàn)了流速峰值。

計算區(qū)域存在類似于土體滲流中沿著優(yōu)先途徑集中流動的優(yōu)先流現(xiàn)象,優(yōu)先流與水流方向大致平行。幾條高速流在某些通道處匯聚成更大的優(yōu)先流,又在某些通道處分散開,導(dǎo)致了水流場分布的不均勻性,而這種優(yōu)先流分布具有隨機性。

圖3 微觀孔隙通道內(nèi)水流場分布(單位：m·s-1)Fig.3 Distribution of water flow in micro-pore channel (unit：m·s-1)

為研究孔隙內(nèi)流速分布特征,選取x=5 mm和x=15 mm兩個剖面上孔隙流速。由圖4可知,接近顆粒邊界的節(jié)點速度幾乎為零,這是由于將顆粒表面設(shè)置成無滑移邊界條件。每條截面上均有多個速度峰值,結(jié)合速度曲線和截線上孔隙通道特征,我們發(fā)現(xiàn)速度峰值往往出現(xiàn)在孔隙連通情況較好的位置,如x=15 mm剖面上的最大速度其所在孔隙位置上四個方向上均有水流流動,并且孔隙通道隙寬較??；而低速區(qū)域一般孔隙連通情況差或位于死端孔隙。正是因為介質(zhì)中顆粒和孔隙分布的隨機性,造成水流速度分布的差異而導(dǎo)致了污染物運移在介質(zhì)中的各向異性。

圖4 剖面上水流分布Fig.4 The distribution of water on the profile

用相同的方法構(gòu)建孔隙度n=0.6的黏土多孔介質(zhì)模型,分別對不同孔隙度模型的兩個截面上的流速做均勻性分析。本文采用相對標(biāo)準(zhǔn)偏差CV(CV值與速度均勻性呈負(fù)相關(guān))、均勻性指數(shù)γv和克里斯琴森均勻系數(shù)CU三個指標(biāo)進行評價[14]。

表1 截面速度均勻性指數(shù)

通過對比以上三個指標(biāo)可發(fā)現(xiàn),其中兩模型不同位置的CV值均超過了1,這表明剖面速度分布十分不均勻；同時兩剖面位置的γv較小,CU接近于0,結(jié)果也表明了均勻性差的特征；而孔隙度變化對水流場均勻性的影響不明顯。我們認(rèn)為這是剖面上孔隙和顆粒的隨機分布所導(dǎo)致的結(jié)果,即多孔介質(zhì)顆粒形態(tài)、分布情況和孔隙連通性、大小共同影響了水流場分布特征。

3 可溶性污染物的濃度場穿透曲線分析

為了更真實地反應(yīng)可溶性染污物在黏土微觀結(jié)構(gòu)中地下水流作用下污染和運移的過程,本文采用脈沖式注入可溶性污染物,濃度C0=1 mol/L。計算總時長為30 s,其中0～10 s為可溶性污染物從左邊進口邊界注入時段,隨后左邊進口邊界改注純水。

3.1 初始水頭壓力對可溶性污染物運移過程的影響分析

為了探究不同水頭壓力條件下可溶性污染物溶質(zhì)運移的情況,將初始壓力分別設(shè)置為5 Pa和10 Pa,在其他條件不變的情況下,進行計算模擬。圖5～6分別顯示了可溶性污染物在初始壓力在P=5 Pa和10 Pa下的運移過程。

圖5 當(dāng)水頭壓力為5 Pa時污染物運移過程(n=0.7)Fig.5 Contaminant transport process (5 Pa,n=0.7)

圖6 當(dāng)水頭壓力為10 Pa時污染物運移過程(n=0.7)Fig.6 Contaminant transport process (10 Pa,n=0.7)

通過對比分析可以看出孔隙的分布、顆粒形狀以及孔隙度對于可溶性污染物的運移情況有較大影響?？扇苄晕廴疚飻U散是一個首先占據(jù)大孔隙通道,然后向狹小孔隙擴散的過程：大孔隙通道內(nèi)水流流速快、流通性好,注入的可溶性污染物在高速水流作用下充滿這類孔隙通道；而大孔隙通道周邊狹小孔隙中的水流流速小,對流作用不明顯,污染物在分子擴散作用下向小孔隙內(nèi)緩慢遷移。在同樣的時間內(nèi),高壓力條件下污染物能夠擴散得更遠(yuǎn),且侵入微小孔隙的能力更強。

同時,其擴散與水流流速也有關(guān)聯(lián),在顆粒表面水流速度較小處,可溶性污染物侵入顆粒表面的能力遠(yuǎn)遠(yuǎn)小于其在水流方向上的運移能力。污染物前緣鋒面的運動方向與孔隙形態(tài)、分布有很大關(guān)聯(lián),鋒面總體上沿著水流方向在連通較好的孔隙中移動,且污染物的運移主方向與水流的優(yōu)勢流方向一致。

在運移過程中,可以發(fā)現(xiàn)一些區(qū)域的濃度變化滯后于濃度中心的運移：在注入可溶性污染物一段時間后,大部分研究區(qū)域為高濃度的污染物,但也有零星分布的區(qū)域內(nèi)污染物濃度較低或為零；而改注純水后,當(dāng)大部分研究區(qū)域濃度降為零時這些區(qū)域污染物濃度仍然較高。這是由于在流動區(qū)域存在可流動區(qū)域(Mobile region)和靜止區(qū)域(Immobile region),靜止區(qū)往往是一些微小空隙或死端孔隙。如圖7所示,在靜止區(qū)域內(nèi)水流幾乎不流動,主要是在分子擴散作用主導(dǎo)下的溶質(zhì)運移過程,因此,在運移過程后期靜止區(qū)域仍殘留高濃度的污染物。此類污染會在黏土介質(zhì)中殘留很長時間,且不會隨著水流的流動而消失。

圖7 不同流動區(qū)域水流場及濃度場分布Fig.7 Flow field and concentration field distribution in different flow areas

3.2 孔隙度對可溶性污染物運移過程的影響分析

使用同樣的方法構(gòu)建n=0.6的黏土多孔介質(zhì)模型,在其他條件不變的情況下進行模擬計算。為保證注入的污染物基本穿過右出口界面,將計算總時長延長至50 s,污染物注入時長不變。對比發(fā)現(xiàn),初始壓力差異的影響與孔隙度n=0.7的模型表現(xiàn)一致,體現(xiàn)了水流場對于污染物遷移的作用。

通過對比相同壓力條件下不同孔隙度模型中的污染物運移過程,可以發(fā)現(xiàn)污染物前緣也沿著水流優(yōu)勢流方向移動,而可溶性污染物沿著大孔隙通道運移的過程更突出。當(dāng)污染物抵達(dá)右邊界時,與孔隙度n=0.7的模型內(nèi)濃度場情況不同,仍有較多區(qū)域保持低濃度狀態(tài),小孔隙度介質(zhì)中不同空隙區(qū)域的濃度差異更明顯,污染物擴散至這些區(qū)域需要更長的時間,滯后現(xiàn)象更明顯。對于孔隙度大的介質(zhì)模型,污染物前緣鋒面抵達(dá)右邊界耗時短,污染物更容易充滿整個介質(zhì)。

圖8 當(dāng)水頭壓力為5 Pa時污染物運移過程(n=0.6)Fig.8 Contaminant transport process (5 Pa,n=0.6)

圖9 當(dāng)水頭壓力為10 Pa時污染物運移過程(n=0.6)Fig.9 Contaminant transport process (10 Pa,n=0.6)

3.3 穿透曲線分析

圖10顯示了在不同壓力邊界條件,黏土微觀結(jié)構(gòu)中可溶性污染物運移的穿透曲線。在本文中,對于不同孔隙度模型下不同初始壓力的算例,穿透曲線的計算位置均選擇在多孔介質(zhì)模型右出水口處。在COMSOL中輸出時長為50 s(步長0.05 s)的無量綱濃度變化情況,并以時間為橫坐標(biāo)繪制穿透曲線,其中無量綱濃度的表達(dá)式為：

(4)

式中：Ly——出口邊界y方向的長度/m；

u——由式(1)計算得到的出口邊界水流速度/(m·s-1)；

C——可溶性污染物的濃度/(mol·L-3)；

C0——污染物初始濃度/(mol·L-3)。

圖10 不同初始壓力下黏土微觀結(jié)構(gòu)中可溶性污染物運移穿透曲線Fig.10 Breakthrough curves of transport of soluble contaminants under different initial pressures

四條曲線反映出污染物的到達(dá)時間存在明顯的先后差異,曲線峰值依次右移且不斷降低,穿透曲線的斜率逐漸減小。對比不同孔隙度的穿透曲線,我們發(fā)現(xiàn)隨著孔隙度降低穿透曲線峰值隨時間軸向右移動且降低,達(dá)到峰值后下降速度減緩,這表明污染物在小孔隙度介質(zhì)中遷移得更慢且同時間內(nèi)殘留污染物更多；在同孔隙度中高初始壓力條件下穿透曲線的斜率更大更陡峭,在達(dá)到峰值后下降更快,這是由于對流作用是影響污染物運移的主要因素,在更大的初始壓力作用下,水流速度更快。可以發(fā)現(xiàn)當(dāng)初始壓力更大時,穿透曲線的峰值濃度更接近于1,也表明了相對于低壓力情況下運移后期在微觀結(jié)構(gòu)中殘留的污染物較少,這一點通過對比圖5(f)和圖6(f)的運移圖像也能說明。

同時,所得的穿透曲線呈現(xiàn)不對稱分布,污染物的流出存在滯后效應(yīng),即在曲線達(dá)到峰值后污染物濃度的下降趨勢小于達(dá)到前的上升趨勢且時間變長,不同情況下曲線均表現(xiàn)出一定拖尾現(xiàn)象,在低壓力、小孔隙度條件下表現(xiàn)最明顯。本文中沒有考慮溶質(zhì)化學(xué)反應(yīng)和吸附現(xiàn)象,造成該現(xiàn)象的原因我們認(rèn)為是介質(zhì)孔隙形態(tài)、顆粒的不均勻分布造成的；孔隙大小和隨機分布情況影響著在該孔隙內(nèi)水流場的分布和濃度場的變化,狹小孔隙和連通性較差的孔隙通道往往在運移過程后期仍殘留較大濃度的污染物,此類孔隙越多,一定時間內(nèi)則殘留的污染物越多且水流帶出污染物的效果越緩慢,使得穿透曲線的拖尾現(xiàn)象越明顯。

4 結(jié)論

(1)由于黏土微觀結(jié)構(gòu)中顆粒與孔隙分布具有隨機性,顆粒表面形狀復(fù)雜,其水流場分布具有很強的隨機性,水流流速分布十分不均勻。

(2)水流場中存在的優(yōu)先流現(xiàn)象,與孔隙大小及其連通情況密切相關(guān),速度峰值往往出現(xiàn)在連通性較好且孔隙寬度較窄的位置。

(3)可溶性污染物的運移過程與水流場密切相關(guān),水流場中的靜止區(qū)域由于水流幾乎不流動而受分子擴散作用主導(dǎo),往往容易殘留高濃度污染物,是污染后期的主要污染殘留位置。

(4)微觀結(jié)構(gòu)中的孔隙隨機分布(如一些微小空隙或死端孔隙)是穿透曲線出現(xiàn)“托尾”現(xiàn)象的主要原因,其“托尾”現(xiàn)象隨初始壓力的增大而減弱。