張學(xué)羿, 竇 智(河海大學(xué)地球科學(xué)與工程學(xué)院, 江蘇 南京 210098)
地下水是可供人類(lèi)利用的重要資源,而隨著社會(huì)經(jīng)濟(jì)持續(xù)發(fā)展和人口快速增長(zhǎng),工業(yè)污水排放、化肥農(nóng)藥施用、垃圾填埋等導(dǎo)致了大量污染物進(jìn)入地下水,地下水污染現(xiàn)已成為全世界共同關(guān)注的問(wèn)題之一。持續(xù)泄露的可溶性污染物在地下水流運(yùn)動(dòng)的長(zhǎng)期作用下可運(yùn)移數(shù)公里甚至數(shù)百公里遠(yuǎn),嚴(yán)重影響了地下水環(huán)境。研究治理地下水污染,首先要了解此類(lèi)污染的過(guò)程和特征,對(duì)其運(yùn)移過(guò)程中水流場(chǎng)和濃度場(chǎng)變化的研究也十分關(guān)鍵[1]。孫家強(qiáng)等[2]以黏土模擬低滲透包氣帶環(huán)境,探究了混溶流體在包氣帶中的遷移規(guī)律;左自波等[3]模擬了土體非飽和帶中污染物在降雨條件下的遷移過(guò)程;余期沖等[4]研究了死端孔隙對(duì)溶質(zhì)運(yùn)移的影響。但以上研究中對(duì)于孔隙微觀結(jié)構(gòu)影響的探討不多。
黏土多孔介質(zhì)結(jié)構(gòu)十分復(fù)雜,并影響著土體的宏觀特征(孔隙度、毛細(xì)壓力、相對(duì)滲透率等),進(jìn)而影響了溶質(zhì)的運(yùn)移。對(duì)于黏土微觀孔隙結(jié)構(gòu)研究的主要問(wèn)題在于多孔介質(zhì)的描述問(wèn)題,即如何采用數(shù)學(xué)的方法對(duì)微觀孔隙結(jié)構(gòu)進(jìn)行定量描述。大多傳統(tǒng)研究中采用只能反映“整體”性能的連續(xù)介質(zhì)模型,故不能反映其微觀特征[5]。對(duì)于微觀研究,首先需構(gòu)建多孔介質(zhì)的微觀結(jié)構(gòu),Tacher L等[6]采用一種離散的減少距離的方法生成圓形、橢圓形的兩相粒狀多孔介質(zhì);Pilotti等[7]采用球體沉降法生成三維多孔介質(zhì);這兩類(lèi)方法生成的孔隙結(jié)構(gòu)比較簡(jiǎn)單,雖然也能反映一定的特征,但是和實(shí)際形態(tài)仍有差距。Dong Hu等[8]、劉向君等[9]、黃家國(guó)等[10]等學(xué)者采用掃描成像技術(shù)分析微觀孔隙結(jié)構(gòu),但成本高昂。Wang M等[11]結(jié)合多孔介質(zhì)構(gòu)造的孔隙生長(zhǎng)模型和叢生理論,提出了四參數(shù)隨機(jī)生長(zhǎng)方法(Quartet Structure Generation Set, QSGS);周瀟等[12]利用該方法探討了滲流速度與孔隙通道大小和土顆粒大小之間的關(guān)系。與傳統(tǒng)的圓形(或橢圓形)顆粒構(gòu)造的黏土微觀結(jié)構(gòu)相比,使用隨機(jī)生長(zhǎng)四參數(shù)生成法構(gòu)造的黏土微觀結(jié)構(gòu)與真實(shí)情況更接近。
黏土多孔介質(zhì)內(nèi)部通道形態(tài)各異,且分布不規(guī)律,在滲流過(guò)程中常出現(xiàn)一些沿著優(yōu)先途徑的集中流動(dòng),統(tǒng)稱(chēng)為優(yōu)先流[13]。本文利用隨機(jī)生長(zhǎng)四參數(shù)法構(gòu)建多孔介質(zhì)模型,通過(guò)Stokes Equation (Stokes方程) 與Advection-Diffusion Equation (對(duì)流擴(kuò)散方程,ADE) 的耦合計(jì)算,模擬黏土多孔介質(zhì)中可溶性污染物的運(yùn)移過(guò)程,旨在探討?zhàn)ね廖⒂^結(jié)構(gòu)對(duì)可溶性污染物運(yùn)移的影響。
為獲得黏土多孔介質(zhì)的微觀孔隙結(jié)構(gòu),采用了隨機(jī)生長(zhǎng)四參數(shù)生成法(QSGS)生成多孔介質(zhì)的孔隙輪廓,并用AutoCAD建立可供有限元分析軟件使用的模型。該方法的構(gòu)造過(guò)程接近于多孔介質(zhì)的生成過(guò)程,與傳統(tǒng)的圓形或橢圓形顆粒構(gòu)成的多孔介質(zhì)相比,其固相顆粒形態(tài)不規(guī)則且大小具有隨機(jī)性,生成的多孔介質(zhì)孔隙包含連通孔隙和非連通孔隙。模型的構(gòu)建主要步驟如下:
(1)利用Matlab實(shí)現(xiàn)Wang M等[11]提出的隨機(jī)生長(zhǎng)四參數(shù)生成法的程序編譯,生成孔隙度n為0.7的多孔介質(zhì)輪廓,保存固相顆粒的輪廓坐標(biāo)數(shù)據(jù):
①確定一定大小的構(gòu)造區(qū)域(本文采用199×199的節(jié)點(diǎn)范圍),初始狀態(tài)該范圍內(nèi)全為孔隙。隨機(jī)分布初始固體顆粒,分布概率Pa應(yīng)小于最終多孔介質(zhì)固相所占比例。游歷構(gòu)造區(qū)域內(nèi)的每個(gè)節(jié)點(diǎn)并在[0,1]內(nèi)取隨機(jī)數(shù),隨機(jī)數(shù)小于或等于Pa的節(jié)點(diǎn)設(shè)置為初始固體顆粒。
②在初始固體顆粒位置按不同方向的生長(zhǎng)概率Pi生成固相顆粒:取8個(gè)生長(zhǎng)方向(i=1,2,…,8),其中包括4個(gè)主要生長(zhǎng)方向(i=1,2,3,4)和4個(gè)次生長(zhǎng)方向(i=5,6,7,8),生長(zhǎng)概率不都相同且P1~4>P5~8,本文中P1~4=0.1、P5~8=0.025,可獲得各向異性的多孔介質(zhì)。各方向相鄰節(jié)點(diǎn)在[0,1]內(nèi)取隨機(jī)數(shù), 隨機(jī)數(shù)小于或等于Pi的節(jié)點(diǎn)生長(zhǎng)為固體顆粒。
③反復(fù)執(zhí)行步驟②,直至孔隙率達(dá)到設(shè)定值n。生成固相顆粒輪廓圖并保存輪廓坐標(biāo)數(shù)據(jù)。
(2)運(yùn)行AutoCAD中的VBA宏程序,將固相顆粒的輪廓坐標(biāo)數(shù)據(jù)自動(dòng)繪制成大小為19.9 mm×19.9 mm的輪廓模型,并添加四周邊界后,保存為可供有限元分析軟件COMSOL使用的DXF文件。
(3)將所得的DXF文件導(dǎo)入COMSOL,模型不同區(qū)域設(shè)置為固相材料或液相材料。
微觀結(jié)構(gòu)黏土的可溶性污染物運(yùn)移問(wèn)題是水流場(chǎng)和濃度場(chǎng)的兩場(chǎng)耦合問(wèn)題,因此,分別對(duì)水流場(chǎng)和濃度場(chǎng)的控制方程進(jìn)行闡述。為了更真實(shí)地模擬黏土微觀結(jié)構(gòu)中水流場(chǎng)的狀態(tài),采用Stokes方程:
-ρgφ+μu=0
(1)
(2)
式中:ρ——水密度/(kg·m-3);
g——重力加速度/(m·s-2);
φ——水頭高度/m;
μ——水的動(dòng)力黏滯系數(shù)/(kg·m-1·s-1);
u——水流速度矢量/(m·s-1)。
通過(guò)上式,可解出黏土多孔介質(zhì)中的非均勻速度場(chǎng),將非均勻速度場(chǎng)代入對(duì)流-擴(kuò)散方程:
(3)
式中:u——由式(1)計(jì)算得到的二維水流場(chǎng)速度矢量;
c——可溶性污染物的濃度/(mol·L-3);
Df——分子擴(kuò)散系數(shù)/(m2·s-1)。
本文假設(shè)Df=1.0×10-9m2/s。需要注意的是,式(3)表示本文考慮的可溶性污染物具有理想特征,即不考慮污染物的化學(xué)反應(yīng)特征,且在顆粒表面無(wú)吸附。
本文旨在利用數(shù)值模擬的方法探討各向異性黏土多孔介質(zhì)對(duì)可溶性污染物運(yùn)移的影響。使用上述DXF格式的圖像文件在COMSOL中建立孔隙度n為0.7的多孔介質(zhì)幾何模型(圖1)。其水流場(chǎng)設(shè)置為:左邊和右邊的進(jìn)口和出口邊界為定壓力邊界,顆粒邊界和上下邊界為無(wú)滑移邊界,水流在5 Pa的壓力下流動(dòng)。
圖1 孔隙度n=0.7的多孔介質(zhì)模型Fig.1 Porous media model (n=0.7)
COMSOL是一款基于多物理場(chǎng)耦合的有限元計(jì)算分析軟件,可以用來(lái)求解線性、非線性問(wèn)題及與時(shí)間有關(guān)的穩(wěn)態(tài)、瞬態(tài)問(wèn)題。本文通過(guò)Stokes Equation (Stokes方程) 與Advection-Diffusion Equation (對(duì)流擴(kuò)散方程,ADE) 的耦合,利用有限元的方法對(duì)計(jì)算區(qū)域進(jìn)行網(wǎng)格剖分細(xì)化、求解,輸出速度場(chǎng)和隨時(shí)間變化的濃度場(chǎng)數(shù)據(jù)及圖形。
黏土微觀結(jié)構(gòu)中的水流場(chǎng)是研究可溶性污染物運(yùn)移的基礎(chǔ)。水流場(chǎng)的控制方程由式(1)和式(2)構(gòu)成。為真實(shí)反映多孔介質(zhì)地下水水流狀態(tài),由進(jìn)口和出口5 Pa壓力差控制的平均水流場(chǎng)速度為0.005 m/s,相應(yīng)的水流場(chǎng)分布如圖2所示。
圖2 水流場(chǎng)分布(單位:m·s-1)Fig.2 Flow field distribution (unit: m·s-1)
圖3顯示了孔隙度n=0.7的黏土微觀孔隙通道中在5 Pa的壓力下水流場(chǎng)的分布情況。其分布十分復(fù)雜,受孔隙連通情況、孔隙形態(tài)、顆粒分布和大小的影響明顯,最快流速達(dá)到0.0478 m/s,而最慢流速接近于0。從圖3選取的微觀孔隙通道的流速等值線圖可以發(fā)現(xiàn),孔隙中的流速基本滿足速度由中心區(qū)域向四周降低的特征。而此通道隙寬較窄且連通性好,共與四個(gè)有水流流動(dòng)的通道相連接,出現(xiàn)了流速峰值。
計(jì)算區(qū)域存在類(lèi)似于土體滲流中沿著優(yōu)先途徑集中流動(dòng)的優(yōu)先流現(xiàn)象,優(yōu)先流與水流方向大致平行。幾條高速流在某些通道處匯聚成更大的優(yōu)先流,又在某些通道處分散開(kāi),導(dǎo)致了水流場(chǎng)分布的不均勻性,而這種優(yōu)先流分布具有隨機(jī)性。
圖3 微觀孔隙通道內(nèi)水流場(chǎng)分布(單位: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兩個(gè)剖面上孔隙流速。由圖4可知,接近顆粒邊界的節(jié)點(diǎn)速度幾乎為零,這是由于將顆粒表面設(shè)置成無(wú)滑移邊界條件。每條截面上均有多個(gè)速度峰值,結(jié)合速度曲線和截線上孔隙通道特征,我們發(fā)現(xiàn)速度峰值往往出現(xiàn)在孔隙連通情況較好的位置,如x=15 mm剖面上的最大速度其所在孔隙位置上四個(gè)方向上均有水流流動(dòng),并且孔隙通道隙寬較??;而低速區(qū)域一般孔隙連通情況差或位于死端孔隙。正是因?yàn)榻橘|(zhì)中顆粒和孔隙分布的隨機(jī)性,造成水流速度分布的差異而導(dǎo)致了污染物運(yùn)移在介質(zhì)中的各向異性。
圖4 剖面上水流分布Fig.4 The distribution of water on the profile
用相同的方法構(gòu)建孔隙度n=0.6的黏土多孔介質(zhì)模型,分別對(duì)不同孔隙度模型的兩個(gè)截面上的流速做均勻性分析。本文采用相對(duì)標(biāo)準(zhǔn)偏差CV(CV值與速度均勻性呈負(fù)相關(guān))、均勻性指數(shù)γv和克里斯琴森均勻系數(shù)CU三個(gè)指標(biāo)進(jìn)行評(píng)價(jià)[14]。
表1 截面速度均勻性指數(shù)
通過(guò)對(duì)比以上三個(gè)指標(biāo)可發(fā)現(xiàn),其中兩模型不同位置的CV值均超過(guò)了1,這表明剖面速度分布十分不均勻;同時(shí)兩剖面位置的γv較小,CU接近于0,結(jié)果也表明了均勻性差的特征;而孔隙度變化對(duì)水流場(chǎng)均勻性的影響不明顯。我們認(rèn)為這是剖面上孔隙和顆粒的隨機(jī)分布所導(dǎo)致的結(jié)果,即多孔介質(zhì)顆粒形態(tài)、分布情況和孔隙連通性、大小共同影響了水流場(chǎng)分布特征。
為了更真實(shí)地反應(yīng)可溶性染污物在黏土微觀結(jié)構(gòu)中地下水流作用下污染和運(yùn)移的過(guò)程,本文采用脈沖式注入可溶性污染物,濃度C0=1 mol/L。計(jì)算總時(shí)長(zhǎng)為30 s,其中0~10 s為可溶性污染物從左邊進(jìn)口邊界注入時(shí)段,隨后左邊進(jìn)口邊界改注純水。
為了探究不同水頭壓力條件下可溶性污染物溶質(zhì)運(yùn)移的情況,將初始?jí)毫Ψ謩e設(shè)置為5 Pa和10 Pa,在其他條件不變的情況下,進(jìn)行計(jì)算模擬。圖5~6分別顯示了可溶性污染物在初始?jí)毫υ赑=5 Pa和10 Pa下的運(yùn)移過(guò)程。
圖5 當(dāng)水頭壓力為5 Pa時(shí)污染物運(yùn)移過(guò)程(n=0.7)Fig.5 Contaminant transport process (5 Pa,n=0.7)
圖6 當(dāng)水頭壓力為10 Pa時(shí)污染物運(yùn)移過(guò)程(n=0.7)Fig.6 Contaminant transport process (10 Pa,n=0.7)
通過(guò)對(duì)比分析可以看出孔隙的分布、顆粒形狀以及孔隙度對(duì)于可溶性污染物的運(yùn)移情況有較大影響??扇苄晕廴疚飻U(kuò)散是一個(gè)首先占據(jù)大孔隙通道,然后向狹小孔隙擴(kuò)散的過(guò)程:大孔隙通道內(nèi)水流流速快、流通性好,注入的可溶性污染物在高速水流作用下充滿這類(lèi)孔隙通道;而大孔隙通道周邊狹小孔隙中的水流流速小,對(duì)流作用不明顯,污染物在分子擴(kuò)散作用下向小孔隙內(nèi)緩慢遷移。在同樣的時(shí)間內(nèi),高壓力條件下污染物能夠擴(kuò)散得更遠(yuǎn),且侵入微小孔隙的能力更強(qiáng)。
同時(shí),其擴(kuò)散與水流流速也有關(guān)聯(lián),在顆粒表面水流速度較小處,可溶性污染物侵入顆粒表面的能力遠(yuǎn)遠(yuǎn)小于其在水流方向上的運(yùn)移能力。污染物前緣鋒面的運(yùn)動(dòng)方向與孔隙形態(tài)、分布有很大關(guān)聯(lián),鋒面總體上沿著水流方向在連通較好的孔隙中移動(dòng),且污染物的運(yùn)移主方向與水流的優(yōu)勢(shì)流方向一致。
在運(yùn)移過(guò)程中,可以發(fā)現(xiàn)一些區(qū)域的濃度變化滯后于濃度中心的運(yùn)移:在注入可溶性污染物一段時(shí)間后,大部分研究區(qū)域?yàn)楦邼舛鹊奈廴疚?但也有零星分布的區(qū)域內(nèi)污染物濃度較低或?yàn)榱悖欢淖⒓兯?當(dāng)大部分研究區(qū)域濃度降為零時(shí)這些區(qū)域污染物濃度仍然較高。這是由于在流動(dòng)區(qū)域存在可流動(dòng)區(qū)域(Mobile region)和靜止區(qū)域(Immobile region),靜止區(qū)往往是一些微小空隙或死端孔隙。如圖7所示,在靜止區(qū)域內(nèi)水流幾乎不流動(dòng),主要是在分子擴(kuò)散作用主導(dǎo)下的溶質(zhì)運(yùn)移過(guò)程,因此,在運(yùn)移過(guò)程后期靜止區(qū)域仍殘留高濃度的污染物。此類(lèi)污染會(huì)在黏土介質(zhì)中殘留很長(zhǎng)時(shí)間,且不會(huì)隨著水流的流動(dòng)而消失。
圖7 不同流動(dòng)區(qū)域水流場(chǎng)及濃度場(chǎng)分布Fig.7 Flow field and concentration field distribution in different flow areas
使用同樣的方法構(gòu)建n=0.6的黏土多孔介質(zhì)模型,在其他條件不變的情況下進(jìn)行模擬計(jì)算。為保證注入的污染物基本穿過(guò)右出口界面,將計(jì)算總時(shí)長(zhǎng)延長(zhǎng)至50 s,污染物注入時(shí)長(zhǎng)不變。對(duì)比發(fā)現(xiàn),初始?jí)毫Σ町惖挠绊懪c孔隙度n=0.7的模型表現(xiàn)一致,體現(xiàn)了水流場(chǎng)對(duì)于污染物遷移的作用。
通過(guò)對(duì)比相同壓力條件下不同孔隙度模型中的污染物運(yùn)移過(guò)程,可以發(fā)現(xiàn)污染物前緣也沿著水流優(yōu)勢(shì)流方向移動(dòng),而可溶性污染物沿著大孔隙通道運(yùn)移的過(guò)程更突出。當(dāng)污染物抵達(dá)右邊界時(shí),與孔隙度n=0.7的模型內(nèi)濃度場(chǎng)情況不同,仍有較多區(qū)域保持低濃度狀態(tài),小孔隙度介質(zhì)中不同空隙區(qū)域的濃度差異更明顯,污染物擴(kuò)散至這些區(qū)域需要更長(zhǎng)的時(shí)間,滯后現(xiàn)象更明顯。對(duì)于孔隙度大的介質(zhì)模型,污染物前緣鋒面抵達(dá)右邊界耗時(shí)短,污染物更容易充滿整個(gè)介質(zhì)。
圖8 當(dāng)水頭壓力為5 Pa時(shí)污染物運(yùn)移過(guò)程(n=0.6)Fig.8 Contaminant transport process (5 Pa,n=0.6)
圖9 當(dāng)水頭壓力為10 Pa時(shí)污染物運(yùn)移過(guò)程(n=0.6)Fig.9 Contaminant transport process (10 Pa,n=0.6)
圖10顯示了在不同壓力邊界條件,黏土微觀結(jié)構(gòu)中可溶性污染物運(yùn)移的穿透曲線。在本文中,對(duì)于不同孔隙度模型下不同初始?jí)毫Φ乃憷?穿透曲線的計(jì)算位置均選擇在多孔介質(zhì)模型右出水口處。在COMSOL中輸出時(shí)長(zhǎng)為50 s(步長(zhǎng)0.05 s)的無(wú)量綱濃度變化情況,并以時(shí)間為橫坐標(biāo)繪制穿透曲線,其中無(wú)量綱濃度的表達(dá)式為:
(4)
式中:Ly——出口邊界y方向的長(zhǎng)度/m;
u——由式(1)計(jì)算得到的出口邊界水流速度/(m·s-1);
C——可溶性污染物的濃度/(mol·L-3);
C0——污染物初始濃度/(mol·L-3)。
圖10 不同初始?jí)毫ο吗ね廖⒂^結(jié)構(gòu)中可溶性污染物運(yùn)移穿透曲線Fig.10 Breakthrough curves of transport of soluble contaminants under different initial pressures
四條曲線反映出污染物的到達(dá)時(shí)間存在明顯的先后差異,曲線峰值依次右移且不斷降低,穿透曲線的斜率逐漸減小。對(duì)比不同孔隙度的穿透曲線,我們發(fā)現(xiàn)隨著孔隙度降低穿透曲線峰值隨時(shí)間軸向右移動(dòng)且降低,達(dá)到峰值后下降速度減緩,這表明污染物在小孔隙度介質(zhì)中遷移得更慢且同時(shí)間內(nèi)殘留污染物更多;在同孔隙度中高初始?jí)毫l件下穿透曲線的斜率更大更陡峭,在達(dá)到峰值后下降更快,這是由于對(duì)流作用是影響污染物運(yùn)移的主要因素,在更大的初始?jí)毫ψ饔孟?水流速度更快。可以發(fā)現(xiàn)當(dāng)初始?jí)毫Ω髸r(shí),穿透曲線的峰值濃度更接近于1,也表明了相對(duì)于低壓力情況下運(yùn)移后期在微觀結(jié)構(gòu)中殘留的污染物較少,這一點(diǎn)通過(guò)對(duì)比圖5(f)和圖6(f)的運(yùn)移圖像也能說(shuō)明。
同時(shí),所得的穿透曲線呈現(xiàn)不對(duì)稱(chēng)分布,污染物的流出存在滯后效應(yīng),即在曲線達(dá)到峰值后污染物濃度的下降趨勢(shì)小于達(dá)到前的上升趨勢(shì)且時(shí)間變長(zhǎng),不同情況下曲線均表現(xiàn)出一定拖尾現(xiàn)象,在低壓力、小孔隙度條件下表現(xiàn)最明顯。本文中沒(méi)有考慮溶質(zhì)化學(xué)反應(yīng)和吸附現(xiàn)象,造成該現(xiàn)象的原因我們認(rèn)為是介質(zhì)孔隙形態(tài)、顆粒的不均勻分布造成的;孔隙大小和隨機(jī)分布情況影響著在該孔隙內(nèi)水流場(chǎng)的分布和濃度場(chǎng)的變化,狹小孔隙和連通性較差的孔隙通道往往在運(yùn)移過(guò)程后期仍殘留較大濃度的污染物,此類(lèi)孔隙越多,一定時(shí)間內(nèi)則殘留的污染物越多且水流帶出污染物的效果越緩慢,使得穿透曲線的拖尾現(xiàn)象越明顯。
(1)由于黏土微觀結(jié)構(gòu)中顆粒與孔隙分布具有隨機(jī)性,顆粒表面形狀復(fù)雜,其水流場(chǎng)分布具有很強(qiáng)的隨機(jī)性,水流流速分布十分不均勻。
(2)水流場(chǎng)中存在的優(yōu)先流現(xiàn)象,與孔隙大小及其連通情況密切相關(guān),速度峰值往往出現(xiàn)在連通性較好且孔隙寬度較窄的位置。
(3)可溶性污染物的運(yùn)移過(guò)程與水流場(chǎng)密切相關(guān),水流場(chǎng)中的靜止區(qū)域由于水流幾乎不流動(dòng)而受分子擴(kuò)散作用主導(dǎo),往往容易殘留高濃度污染物,是污染后期的主要污染殘留位置。
(4)微觀結(jié)構(gòu)中的孔隙隨機(jī)分布(如一些微小空隙或死端孔隙)是穿透曲線出現(xiàn)“托尾”現(xiàn)象的主要原因,其“托尾”現(xiàn)象隨初始?jí)毫Φ脑龃蠖鴾p弱。