考慮雙重沉積模式的多孔介質(zhì)中顆粒三維遷移模型研究

劉曉明 ，涂樹杰 ，陽(yáng)棟 ，楊澤曦 ，劉濤

（1.湖南大學(xué) 土木工程學(xué)院，湖南 長(zhǎng)沙 410082；2.中國(guó)建筑第五工程局有限公司，湖南 長(zhǎng)沙 410004；3.中交公路規(guī)劃設(shè)計(jì)院有限公司，北京 100088）

多孔介質(zhì)中懸浮顆粒的遷移規(guī)律研究在地下水回灌［1］、核廢料處置［2］、反濾層設(shè)計(jì)［3］、地下污染物遷移［4］等領(lǐng)域具有重要意義，是當(dāng)前環(huán)境巖土工程領(lǐng)域一個(gè)重要的研究方向［5］.

國(guó)內(nèi)外學(xué)者對(duì)多孔介質(zhì)中顆粒遷移過(guò)程及機(jī)理的研究已開(kāi)展了較多的工作.Ahfir 等［6］提出了一種描述脈沖注入顆粒物在多孔介質(zhì)中的遷移及沉積特性數(shù)學(xué)模型.Bedrikovetsky 等［7］提出了考慮彌散效應(yīng)及顆粒物尺寸阻滯作用的深層顆粒物遷移數(shù)學(xué)模型，并基于土柱試驗(yàn)與傳統(tǒng)模型進(jìn)行了對(duì)比.陳云敏等［8］建立了污染物在層狀土中的一維擴(kuò)散模型，通過(guò)分離變量法得到了解析解，闡明了污染物遷移擊穿防污屏障的內(nèi)在機(jī)理.試驗(yàn)研究方面，陳星欣等［9］通過(guò)室內(nèi)土柱試驗(yàn)探討了懸浮顆粒的濃度對(duì)其在飽和多孔介質(zhì)中遷移和沉積特性的影響.劉泉聲等［10］通過(guò)模擬試驗(yàn)研究了顆粒粒徑對(duì)多孔介質(zhì)中懸浮顆粒遷移-沉積特性的影響，并根據(jù)粒徑比不同將懸浮顆粒在多孔介質(zhì)中的遷移-沉積類型劃分為濾餅過(guò)濾型、遷移-沉積型、自由遷移型三種.

經(jīng)典顆粒遷移模型認(rèn)為多孔介質(zhì)中的懸浮顆粒沉積是吸附作用造成的.Sakthivadivel 等［11］發(fā)現(xiàn)，當(dāng)懸浮顆粒粒徑與多孔介質(zhì)粒徑比值大于0.05 時(shí)，存在顯著的篩濾作用（即懸浮顆粒粒徑大于多孔介質(zhì)的孔隙吼道尺寸，形成的篩濾沉積）.而Bradford等［12］研究表明，懸浮顆粒中位粒徑與多孔介質(zhì)粒徑比值大于0.005的低粒徑比時(shí)，篩濾作用依然發(fā)揮著重要的作用.因此，在低粒徑比時(shí)，同時(shí)考慮顆粒篩濾作用和吸附作用的預(yù)測(cè)結(jié)果與試驗(yàn)結(jié)果更為接近.處理因篩濾作用而產(chǎn)生的額外沉積機(jī)制的一種方法是假設(shè)存在雙重沉積模式，即篩濾效應(yīng)沉積與吸附效應(yīng)沉積.對(duì)于雙重沉積模式的研究，Simoni 等［13］通過(guò)將顆粒劃分為兩個(gè)子群，并賦予不同的碰撞效率值，成功描述了沉積顆粒濃度隨深度非指數(shù)分布現(xiàn)象.Tufenkji 等［14］提出了考慮“快”和“慢”顆粒沉積共同影響的雙沉積模式，提高了顆粒遷移沉積的預(yù)測(cè)精度.Katzourakis 等［15］通過(guò)拉普拉斯變換方法求得了雙重沉積模式下顆粒遷移一維解析解，更好地?cái)M合了現(xiàn)有的實(shí)驗(yàn)數(shù)據(jù).由于考慮雙重沉積后，問(wèn)題求解變得復(fù)雜，這些研究都是在一維條件下進(jìn)行的，不符合顆粒遷移的空間特性.

本文采用雙重沉積模式下的沉積動(dòng)力學(xué)方程，考慮附著顆粒的再釋放和大顆粒的篩濾效應(yīng)對(duì)顆粒遷移的經(jīng)典模型進(jìn)行修正，建立可同時(shí)考慮篩濾效應(yīng)沉積與吸附效應(yīng)沉積的多孔介質(zhì)中顆粒三維遷移模型；然后通過(guò)拉普拉斯變換和傅里葉變換及其逆變換求得顆粒三維遷移模型的通解，得到點(diǎn)源和圓形面源注入條件下半無(wú)限空間顆粒遷移的解析表達(dá)式，并通過(guò)解的退化及對(duì)突破曲線的參數(shù)反演驗(yàn)證了解的正確性；最后，基于導(dǎo)得的解析解，對(duì)圓形面源恒定濃度注入情況下的水動(dòng)力彌散系數(shù)、篩濾系數(shù)、吸附系數(shù)、釋放系數(shù)等參數(shù)對(duì)顆粒遷移過(guò)程的影響機(jī)理進(jìn)行了詳細(xì)的分析.

1 顆粒遷移三維控制方程及通解

1.1 控制方程建立

考慮飽和均勻流場(chǎng)中的對(duì)流、三維水動(dòng)力彌散和顆粒沉積效應(yīng)，多孔介質(zhì)液相中懸浮物顆粒運(yùn)輸?shù)馁|(zhì)量平衡方程為：

式中：C是液相中的懸浮物濃度，量綱為ML-3；t為時(shí)間，量綱為T；v為孔隙流體平均流速，量綱為L(zhǎng)T-1；x為平行于流動(dòng)方向的空間坐標(biāo)，量綱為L(zhǎng)；y，z為垂直于流動(dòng)方向的空間坐標(biāo)，量綱為L(zhǎng)；Dx，Dy，Dz分別為x，y，z方向上的水動(dòng)力彌散系數(shù)，量綱為L(zhǎng)2T-1；n為孔隙率，量綱為1；C*為單位體積內(nèi)沉積的顆?？傎|(zhì)量，量綱為ML-3.C與C*均為變量t、x、y、z的函數(shù).

影響多孔介質(zhì)中顆粒沉積的機(jī)制主要有兩種，即大顆粒的篩濾作用和小顆粒的吸附作用［11，16］.圖1是兩種顆粒篩濾作用和吸附作用的示意圖.

圖1 篩濾作用與吸附作用示意圖Fig.1 Schematic diagram of adsorption and sieve filtration

單位體積內(nèi)沉積的顆?？傎|(zhì)量C*與篩濾作用沉積項(xiàng)C1*和吸附作用沉積項(xiàng)C2*的關(guān)系為：

其中，篩濾作用沉積項(xiàng)C1*被認(rèn)為是不可逆的，而吸附作用沉積項(xiàng)C2*是可逆的.顆粒沉積動(dòng)力學(xué)方程可以描述為：

式中：kstr為篩濾系數(shù)，量綱為T-1，與粒徑比呈指數(shù)相關(guān)［17］；kat=3(1 -n)vαη/2d50為顆粒吸附系數(shù)，量綱為T-1，其中，d50為多孔介質(zhì)的平均粒徑，α為顆粒吸附效率，η為收集效率，其計(jì)算公式見(jiàn)文獻(xiàn)［18］；kde為顆粒釋放系數(shù)，量綱為T-1.

式（1）右側(cè)第一項(xiàng)為對(duì)流項(xiàng)，體現(xiàn)了顆粒遷移的滲流作用.式（1）右側(cè)最后一項(xiàng)為顆粒沉積項(xiàng)，由沉積動(dòng)力學(xué)方程式（2）～（4）共同控制，本文的顆粒遷移模型是由式（1）～（4）聯(lián)合建立的，因此考慮了滲流和顆粒遷移/沉積的耦合作用.

假設(shè)初始時(shí)多孔介質(zhì)中無(wú)顆粒，顆粒從半無(wú)限空間表面注入，相應(yīng)初始條件和邊界條件可表示為：

初始條件式（5）（6）表示多孔介質(zhì)中既無(wú)沉積顆粒，也無(wú)懸浮顆粒.邊界條件式（7）表示在半無(wú)限空間表面注入強(qiáng)度隨時(shí)間變化的顆粒注入源.邊界條件式（8）～（10）為半無(wú)限域的理想情況.

需要指出的是，考慮極限情景下細(xì)顆?？赡軣o(wú)法穿透多孔介質(zhì)，本文提出的顆粒遷移模型對(duì)顆粒的粒徑應(yīng)該有限定要求，參照文獻(xiàn)［19-20］，限定細(xì)顆粒與多孔介質(zhì)的粒徑比小于0.1.當(dāng)細(xì)顆粒粒徑大于100 μm 時(shí)，顆粒的重力沉降較為明顯，而本文模型中只考慮了對(duì)流與水動(dòng)力彌散作用，因此應(yīng)用本文模型時(shí)，細(xì)顆粒粒徑應(yīng)小于100 μm.

1.2 變換域的解

首先，通過(guò)拉普拉斯變換求得滿足初始條件式（5）下方程式（3）的解為：

同理可求得方程（4）的解為：

將式（11）（12）對(duì)t求導(dǎo)后代入式（1）得到液相中懸浮物濃度表達(dá)式為：

式（13）需采用積分變換方法求解.對(duì)t和x進(jìn)行拉普拉斯變換，變換變量為s和ω，對(duì)y和z進(jìn)行傅里葉變換，變換變量為α和β.利用式（6）（7）可得變換域上的解為：

1.3 時(shí)空域的解

時(shí)空域中懸浮物濃度是通過(guò)對(duì)拉普拉斯域和傅里葉域的解進(jìn)行反演得到的.首先對(duì)ω求拉普拉斯逆變換，使用拉普拉斯變換的位移性質(zhì)、卷積定理和拉普拉斯變換表得到［21］：

將式（15）兩邊同乘2B?exp[(A-B)x]并運(yùn)用邊界條件（8）可以得到：

式中：L-1為拉普拉斯逆變換算子為任意函數(shù)f0(t)的拉普拉斯變換；a為任意常數(shù)；I0為零階第一類修正貝塞爾函數(shù).

假設(shè)：

式中：I1為一階第一類修正貝塞爾函數(shù).

對(duì)exp(-Et)進(jìn)行兩次傅里葉逆變換，得到：

最后，根據(jù)傅里葉變換的卷積定理，由式（17）（24）（27）得到考慮雙重沉積條件下懸浮顆粒濃度的三維解析解：

2 半無(wú)限空間表面注入顆粒問(wèn)題解析解

將邊界條件式（7）的顆粒注入函數(shù)g(t，y，z)可以表示為：

式中：W(y，z)為顆粒注入源的幾何形態(tài)；G(t)為顆粒的注入方式.

2.1 點(diǎn)源注入顆粒

當(dāng)顆粒注入源為點(diǎn)源時(shí)，顆粒注入源的幾何形態(tài)可表示為：

式中：δ(?)為Dirac delta 函數(shù)；y0和z0為點(diǎn)源的坐標(biāo).

將式（31）（32）代入式（28），利用Dirac delta 函數(shù)的篩選性質(zhì)可得點(diǎn)源注入情況下顆粒遷移的解析解：

2.2 面源注入顆粒

假定半無(wú)限多孔介質(zhì)表面局部區(qū)域D內(nèi)，在t=0時(shí)刻連續(xù)注入濃度隨時(shí)間t變化G(t)的顆粒.那么多孔介質(zhì)內(nèi)部顆粒物的遷移過(guò)程相當(dāng)于無(wú)數(shù)個(gè)強(qiáng)度為G(t)dydz點(diǎn)源注入情形下顆粒物遷移過(guò)程在空間上的疊加.于是，多孔介質(zhì)內(nèi)顆粒的濃度可由式（33）在空間域上通過(guò)積分求得，此即源函數(shù)法的基本思想［23］.因此

考慮半無(wú)限體表面圓形區(qū)域連續(xù)注入濃度隨時(shí)間變化的顆粒（圓心為O，半徑為a）的情形.假定垂直水流方向水動(dòng)力彌散系數(shù)相等（即Dy=Dz=Dr），而平行水流方向水動(dòng)力彌散系數(shù)為Dx，此時(shí)顆粒濃度場(chǎng)可采用圓柱坐標(biāo)（x，r）來(lái)表示，式（35）可寫成：

將式（36）對(duì)變量φ進(jìn)行積分，并注意貝塞爾函數(shù)的性質(zhì)，整理后可得：

3 顆粒遷移解析解驗(yàn)證

3.1 解的一維退化

當(dāng)顆粒注入源半徑a=∞時(shí)，式（38）方便地退化成一維情形.

式中：tp為注入時(shí)間；Ω(t，x，r)為式（40）中G(τ)=C0的特殊情況.

對(duì)于瞬時(shí)注入顆粒情況，有：

式中：I=m/(nvS)為面源強(qiáng)度，量綱為MTL-3；m為注入顆粒質(zhì)量，量綱為M；S為面源面積，量綱為L(zhǎng)2；t'為顆粒注入的時(shí)刻，量綱為T.

利用Dirac delta 函數(shù)的篩選性質(zhì)，式（40）可以轉(zhuǎn)化為：

令式（43）中t'=0，kstr=0 和kde=0，即可得到基于經(jīng)典對(duì)流彌散模型的顆粒遷移解析解［24-25］.

3.2 驗(yàn)證1：恒定濃度注入

本節(jié)采用式（40）（41）對(duì)Syngouna 等［26］的試驗(yàn)進(jìn)行分析，并與單重沉積模式顆粒遷移解析解進(jìn)行對(duì)比，以檢驗(yàn)雙重沉積模式顆粒遷移解析解的正確性.Syngouna 等［26］進(jìn)行了兩種流速（v=1.21 cm/min 和0.74 cm/min）下兩種黏土顆粒（高嶺石（KGa-1b）和蒙脫石（STx-1b））在填充玻璃珠的水平柱中的運(yùn)輸試驗(yàn).柱直徑D為2.5 cm，長(zhǎng)L為30 cm，孔隙率n為0.42，其他相關(guān)參數(shù)見(jiàn)表1.

表1 恒定濃度注入時(shí)顆粒運(yùn)輸參數(shù)Tab.1 Particle transport parameters during constant concentration injection

試驗(yàn)數(shù)據(jù)的分析是通過(guò)PEST（Parameter Estima?tion）自動(dòng)率定套件程序進(jìn)行的.PEST 套件是基于GML（Gauss-Marquardt-Levenberg）算法開(kāi)發(fā)的獨(dú)立于模型參數(shù)估計(jì)和不確定性分析的綜合軟件，具有逆海森方法和最速下降法的優(yōu)點(diǎn)，可以通過(guò)較少的模型運(yùn)行次數(shù)，得到最優(yōu)的參數(shù)結(jié)果［27］.采用PEST程序進(jìn)行參數(shù)率定共分為三個(gè)步驟：

第一步為確定目標(biāo)函數(shù).本文采用單目標(biāo)進(jìn)行參數(shù)優(yōu)化，以懸浮物理論流出濃度與實(shí)測(cè)流出濃度差的平方和最小為優(yōu)化目標(biāo)，目標(biāo)函數(shù)公式為：

式中：f為目標(biāo)函數(shù)；i為時(shí)段序號(hào)；m為總時(shí)段數(shù)；Cmod，i為第i時(shí)段的懸浮物理論流出濃度；Cobs，i表示第i時(shí)段的懸浮物實(shí)測(cè)流出濃度.

第二步為選擇參數(shù)范圍，為需要率定的參數(shù)選擇合理的范圍.需要率定的參數(shù)為Dx，kstr，kat，kde4個(gè).

第三步為輸入懸浮物實(shí)測(cè)流出濃度數(shù)據(jù)并執(zhí)行迭代，輸出率定的參數(shù)結(jié)果以及懸浮物理論流出濃度.

試驗(yàn)的突破曲線和單、雙重沉積模式顆粒遷移解析解的計(jì)算結(jié)果如圖2 所示（x=L=30 cm），相關(guān)參數(shù)見(jiàn)表1.對(duì)于Syngouna 等［26］考慮的兩種流速情況，雙重沉積模式解析解成功地?cái)M合了KGa-1b 和STx-1b的突破曲線（x=L=30 cm），而相同條件下單重沉積模式下的解對(duì)該組試驗(yàn)數(shù)據(jù)的效果不理想，說(shuō)明雙重沉積模式下顆粒遷移解析解效果更優(yōu).

圖2 恒定濃度注入時(shí)顆粒遷移突破曲線Fig.2 Particle migration breakthrough curve at constant concentration injection

3.3 驗(yàn)證2：瞬時(shí)注入

本節(jié)采用瞬時(shí)注入顆粒遷移解析式（43）對(duì)Bai等［20］的試驗(yàn)進(jìn)行分析，并與單重沉積模式顆粒遷移解析解進(jìn)行對(duì)比，以檢驗(yàn)雙重沉積模式顆粒遷移解析解的性能.試驗(yàn)數(shù)據(jù)的分析方法與3.2 節(jié)一致.Bai 等［20］進(jìn)行了兩種流速（v=0.384 cm/s 和0.576 cm/s）下兩種粒徑的球形二氧化硅顆粒（中值粒徑D50=1.02 μm 和47 μm）在石英砂（D50=2.25 mm）柱中的運(yùn)輸試驗(yàn).柱直徑D為7 cm，長(zhǎng)L為30 cm，孔隙率為0.45.

定義相對(duì)濃度CR=CVP/m，其中VP為柱的孔隙體積.試驗(yàn)的突破曲線和單、雙重沉積模式顆粒遷移解析解的計(jì)算結(jié)果如圖3 所示（x=L=30 cm），相關(guān)參數(shù)見(jiàn)表2.由圖3 可知，解析式（43）可以成功模擬瞬時(shí)注入時(shí)顆粒遷移的突破曲線，且在顆粒粒徑較大（D50=47 μm）時(shí)的性能比單重沉積模式下顆粒遷移解析解性能更好.通過(guò)分析可知，當(dāng)懸浮物粒徑為1.02 μm 時(shí)，粒徑比遠(yuǎn)小于0.005，顆粒物沉積由吸附作用主導(dǎo)；當(dāng)懸浮物粒徑為47 μm 時(shí)，粒徑比大于0.005，篩濾作用在顆粒物沉積過(guò)程中發(fā)揮著重要的作用，此時(shí)采用雙重沉積模式的顆粒遷移模型效果更好.

表2 瞬時(shí)注入時(shí)顆粒運(yùn)輸參數(shù)Tab.2 Particle transport parameters during instantaneous injection

圖3 瞬時(shí)注入時(shí)顆粒遷移突破曲線Fig.3 Particle migration breakthrough curve during instantaneous injection

4 圓形面源注入時(shí)顆粒遷移規(guī)律

本節(jié)以圓形面源注入時(shí)的雙重沉積模式顆粒遷移解析解［式（38）］來(lái)研究多孔介質(zhì)中顆粒三維遷移規(guī)律.選擇了恒定濃度注入這一典型的顆粒注入方式進(jìn)行分析.

4.1 恒定濃度注入時(shí)顆粒遷移規(guī)律

參考文獻(xiàn)［28］～［30］擬定計(jì)算參數(shù)：徑向水動(dòng)力彌散系數(shù)Dr=0.1 cm2/s，軸向水動(dòng)力彌散系數(shù)Dx=1 cm2/min，滲流速度v=0.2 cm/min，篩濾系數(shù)kstr=0.02/min，顆粒吸附系數(shù)kat=0.01/min，顆粒釋放系數(shù)kde=0.001/min.除特別說(shuō)明外，本文算例均采用前述數(shù)值進(jìn)行計(jì)算.

圖4 為半無(wú)限多孔材料表面存在恒定濃度注入時(shí)懸浮顆粒物濃度演化（tp=200 min，a=5 cm）.引入無(wú)量綱濃度C/C0，由圖4 可知，C/C0在深度方向和水平方向隨著時(shí)間逐漸變化，由于深度方向存在滲流作用，顆粒物在深度方向的遷移速度明顯快于水平方向.當(dāng)時(shí)間較大（如tp=160 min）時(shí)，顆粒遷移過(guò)程逐漸穩(wěn)定，顆粒物濃度在空間上的分布保持不變.

圖4 恒定濃度注入時(shí)懸浮顆粒物濃度演化（tp=200 min，a=5 cm）Fig.4 Concentration evolution of suspended particles at constant concentration injection（tp=200 min，a=5 cm）

圖5 給出了恒定濃度注入時(shí)顆粒遷移過(guò)程（tp=100 min，a=5 cm）.由圖5（a）可知，在顆粒注入期間（ttp時(shí)），顆粒物濃度隨深度增加依次達(dá)到峰值并快速下降.注入的顆粒物在對(duì)流和彌散作用下逐漸向多孔介質(zhì)的深度和水平方向遷移并沉積.因此，在一定深度下，多孔介質(zhì)中的顆粒物濃度會(huì)逐漸增加.深度越大，顆粒物濃度達(dá)到峰值的時(shí)間就越長(zhǎng)，峰值濃度就越小.

圖5 恒定濃度注入時(shí)顆粒遷移過(guò)程（tp=100 min，a=5 cm）Fig.5 Particle migration process during constant concentration injection（tp=100 min，a=5 cm）

圖5（b）為恒定濃度注入時(shí)（tp=100 min，a=5 cm）多孔介質(zhì)中顆粒物濃度沿深度的分布.在顆粒注入期間（即t=20 min→40 min→80 min 時(shí)），多孔介質(zhì)中顆粒物濃度隨時(shí)間增大而增大，隨深度增加而減小并趨于0.在顆粒停止注入后（即t=120 min→140 min→180 min 時(shí)），多孔介質(zhì)表面顆粒濃度迅速下降，濃度峰值隨時(shí)間增加而減小，峰值濃度所處的深度隨時(shí)間增加而增加.

4.2 參數(shù)敏感性分析

為了進(jìn)一步研究恒定濃度注入時(shí)顆粒遷移模型對(duì)各模型參數(shù)的敏感性，分別對(duì)參數(shù)Dx、kstr、kat、kde取不同值時(shí)的懸浮物突破曲線進(jìn)行了計(jì)算（x=10 cm，tp=100 min），參數(shù)的敏感性分析結(jié)果見(jiàn)圖6.

圖6 恒定濃度注入時(shí)顆粒遷移敏感性分析（x=10 cm，tp=100 min）Fig.6 Sensitivity analysis of particle migration at constant concentration injection（x=10 cm，tp=100 min）

由圖6 可知，懸浮物突破曲線存在峰值，顆粒濃度峰值隨著參數(shù)Dx和kde值的增加而增加，隨著kstr和kat值的增加而減小.通過(guò)分析可知，水動(dòng)力彌散系數(shù)Dx值越大，顆粒彌散效應(yīng)越顯著，突破時(shí)間越快，而停止注入后（t>tp）的濃度越小.此外，kstr和kat越大，顆粒物在固體基質(zhì)上的沉積越多，從而液相中顆粒峰值濃度越??；kde值越大，沉積顆粒從固體基質(zhì)中脫離越多，從而導(dǎo)致液相中峰值濃度越大.

5 結(jié)論

1）建立了考慮篩濾效應(yīng)與吸附效應(yīng)的雙重沉積模式顆粒三維遷移模型，利用拉普拉斯變換、傅里葉變換及其變換反演，導(dǎo)出了在一維滲流和三維彌散效應(yīng)條件下顆粒遷移問(wèn)題通解.給出了點(diǎn)源注入及面源注入顆粒時(shí)顆粒物濃度時(shí)空分布的解析表達(dá)式，并通過(guò)解的退化及對(duì)突破曲線的參數(shù)反演驗(yàn)證了解的正確性.

2）圓形面源恒定濃度注入情況下，在顆粒注入期間，多孔介質(zhì)中顆粒物濃度隨時(shí)間增大而增大，隨深度增加而減小并趨于0.在顆粒停止注入后，多孔介質(zhì)中顆粒物濃度隨時(shí)間增大而減小，在深度上存在濃度峰值，濃度峰值所處的深度隨時(shí)間增加而增加.

3）對(duì)面源注入情況下水動(dòng)力彌散系數(shù)、沉積參數(shù)等的影響機(jī)理分析表明，水動(dòng)力彌散效應(yīng)加速了顆粒物的遷移，使顆粒突破時(shí)間變快，峰值濃度變高.篩濾系數(shù)kstr和顆粒吸附系數(shù)kat越大，顆粒釋放系數(shù)kde越小，顆粒物在固體基質(zhì)上的沉積越多，懸浮顆粒峰值濃度越小.