程鵬達 朱心廣,? 馮 春,2) 王曉亮

?(中國科學(xué)院力學(xué)研究所,北京 100190)

?(中國科學(xué)院大學(xué)工程科學(xué)學(xué)院,北京 100049)

??(北京理工大學(xué),北京 100081)

引言

水體被認為是由水、溶解性物質(zhì)、再懸浮性物質(zhì)、水生生物和底泥組成的自然綜合體[1].污染物質(zhì)進入水體后會沉積到底泥中并逐漸富集,形成具有一定厚度的含有各種污染物的沉積物層,使得底泥成為污染物質(zhì)的蓄積庫.當(dāng)水力條件比較復(fù)雜時,沉積物會在一定的水流作用下發(fā)生沖刷及再懸浮等過程,使大量的污染物被重新釋放出來,造成水體的二次污染.當(dāng)外源性污染逐步得到控制后,污染底泥的內(nèi)源釋放效應(yīng)越來越明顯.例如,上海市環(huán)境科學(xué)院對淀山湖污染負荷輸入進行了實測,發(fā)現(xiàn)除上、下游來水負荷之外,底泥污染物釋放已與大氣干濕沉降、地表徑流各占1/3 左右[2-3].污染物在水環(huán)境中的遷移轉(zhuǎn)化、歸宿和泥沙運動密切相關(guān),水動力學(xué)條件和泥沙性質(zhì)共同影響著污染物的遷移轉(zhuǎn)化過程[4].一方面,水體中再懸浮泥沙的運動、輸移直接影響著污染物時空分布;另一方面,污染物隨泥沙沉積到底床中,沉積底泥成為潛在的重要污染“源”或“匯”,而在一定的水動力擾動和環(huán)境條件下會發(fā)生“源”、“匯”機制的轉(zhuǎn)換[5-7].底泥再懸浮污染物釋放過程中,首先水動力學(xué)條件增強了泥?水界面附近污染物的擴散和混合能力,然后泥?水界面水流的剪應(yīng)力和紊動強度會導(dǎo)致受污染底泥顆粒的再懸浮,進而引起污染物向上覆水體大量釋放[8].因此,底泥再懸浮引起上覆水體流場特性變化,是正確理解和完整描述非吸附性污染物在水環(huán)境中的遷移轉(zhuǎn)化關(guān)鍵.

以往底泥污染物釋放研究主要集中在湖泊等相對靜止或流速較低的水體當(dāng)中,將底泥視為多孔介質(zhì),系統(tǒng)分析和研究了污染物的滲流場?濃度場耦合輸運.朱廣偉等[9]在野外調(diào)查和實驗室試驗分析中觀測到太湖沉積物?水界面處物質(zhì)交換主要發(fā)生在沉積物表層5～10 cm 范圍內(nèi),水動力作用(風(fēng)浪、湖流) 導(dǎo)致的表層底泥再懸浮將使孔隙水中的營養(yǎng)鹽得以大量釋放進入上覆水體.王道增等[10-12]采用水槽實驗的方法研究了底泥污染物在動態(tài)水流條件下(不同水深、流速)的釋放擴散規(guī)律及其對上覆水體濁度、COD,BOD 等水質(zhì)指標(biāo)的影響.而對于底泥再懸浮,傳統(tǒng)泥沙再懸浮分析大多集中于泥沙和水流相互作用[13-14].底泥再懸浮后污染物釋放研究以實驗觀測為主,水動力學(xué)釋放機理研究較少.Fetters 等[15]基于伊利諾伊州德普湖和緬因州樸次茅斯海軍造船廠的底泥再懸浮實驗觀察,指出沉積物短期(< 4 h)再懸浮,引起金屬污染物向上覆水體釋放,釋放量與水動力學(xué)條件、氧化還原反應(yīng)、生化反應(yīng)有關(guān).Matisoff等[16]基于加拿大溫尼伯湖底泥實驗研究,指出湖泊藻類水華大部分與懸浮顆粒有關(guān),營養(yǎng)鹽在再懸浮后的遷移和歸宿是很復(fù)雜的.大多數(shù)(95%～99%)懸浮物質(zhì)來自超過23 年累計的高達7 cm 的底部沉積物,該層沉積物在幾十年內(nèi)仍然是一個重要的、活躍的內(nèi)部營養(yǎng)鹽負荷源.Lepage 等[17]基于2012 年6 月法國羅納河上游3 座大壩的防淤沖洗觀測數(shù)據(jù),指出沖洗期間的平均懸浮顆粒物濃度比2011 年—2016 年記錄的洪水事件期間平均高出6～8 倍,但污染物濃度低于基流或洪水期間,這可能是因為沖洗過程中污染物僅來自庫區(qū)再懸浮底泥,而洪水期污染物來源更為廣泛.隨底泥再懸浮,固相懸浮顆粒濃度的增加比較直觀可見,而顆粒懸浮后污染物釋放的過程易被忽略了.因此,底泥再懸浮污染物釋放總量容易被低估或夸大.為了解水體底泥再懸浮污染物釋放的物理過程和影響因素,本文建立上覆水體?底泥?污染物的耦合力學(xué)模型,研究不同水動力學(xué)條件下,流場速度、泥沙顆粒體積分數(shù)、污染物濃度、湍動能以及時間等參數(shù)之間的關(guān)系,得到流場特性與污染物釋放之間的關(guān)系,并得到對流和湍流擴散兩種作用對上覆水體污染物濃度的影響.

1 控制方程

再懸浮泥沙的主要成分是比表面積較大的細顆粒物,是污染物的主要的輸運載體[18].前期實驗測量得到受污染底泥屬于粉沙和黏土混合的細顆粒泥沙,小于5 μm 的土顆粒占8%以上,污染底泥的中值粒徑D50在25～50 μm 左右,但表層浮泥和淤泥的D50都小于30μm.本文中選擇D50=30μm 為底泥起動的代表粒徑,由于底泥粒徑較小,在本文上覆水體?底泥?污染物模型中再懸浮底泥被視為一種懸浮液.基于大量實驗數(shù)據(jù)分析[19-20],該類混合細顆粒的懸浮液黏度可以表示為顆粒體積分數(shù)的函數(shù).因此,可以認為該類懸浮液(水?顆粒) 是具有宏觀性質(zhì)(例如密度和黏度) 的單一流動連續(xù)體.模型假設(shè)如下:(1)每相的密度大致恒定;(2)共享相同的壓力;(3)與宏觀流動的時間尺度相比,顆粒弛豫時間較短.

再懸浮液的密度和黏度分別由下式給出

式中,ρ1,ρ2和μ1,μ2分別是流體和顆粒的密度和黏度,φ 是顆粒的體積分數(shù).

考慮到再懸浮液可以作為連續(xù)介質(zhì)處理,流場由不可壓縮質(zhì)量和動量守恒方程控制.

式中,u為速度,ρ 為再懸浮液密度,t為時間,p為壓力,τ為黏性應(yīng)力和雷諾應(yīng)力總和,g是重力加速度.τ=μγ,其中μ 是再懸浮液黏度,γ=?u+?uT.f為流體?顆粒的阻力,與顆粒直徑和雷諾數(shù)相關(guān).此外,本文流場計算時,采用標(biāo)準k–ε 湍流模型[21],模型中包括一些經(jīng)驗參數(shù),這些參數(shù)的值是通過適用于各種湍流的數(shù)據(jù)多次迭代得出的,其中Cμ=0.09,σk=0.10,σε=1.3,C1ε=1.44,C2ε=1.92.

底泥顆粒的輸運方程由下式?jīng)Q定

式中,Nφ是粒子的總擴散通量,受兩種不同的機制的影響,第一是由粒子碰撞頻率梯度引起的通量Nc,Nc=?Kca2φ?(γφ);第二是由懸浮黏度梯度引起的通量Nμ,Nμ=?Kμa2(γφ2/μ)?μ.其中Kc和Kμ是經(jīng)驗擴散系數(shù),a是顆粒半徑,γ 是局部剪切速率.Nc驅(qū)動顆粒向剪切速率和顆粒濃度梯度的相反方向運動;Nμ驅(qū)動顆粒向懸浮黏度梯度的相反方向運動.本文主要考慮平衡浮力的細顆粒懸浮,從而忽略了布朗運動和重力沉降引起的附加擴散通量.

污染物濃度的輸運方程由下式?jīng)Q定

式中,c是相對濃度,為無量綱參數(shù);D為擴散系數(shù);R為單位時間單位體積空間內(nèi)因化學(xué)反應(yīng)生成或消失量,為源和匯項.

底泥顆粒再懸浮液的黏度通常寫為顆粒體積分數(shù)的函數(shù)[23-25].本文采用MPQ 模型[22-23]

式中,φm是固相顆粒最大填充濃度,取值為0.62.

2 數(shù)值模型

2.1 幾何模型

前期底泥再懸浮污染物釋放實驗采用循環(huán)水槽,如圖1 所示.水槽寬度為0.25 m,沙床長度為4 m,水深為0.10 m,底泥厚度為0.08 m.

圖1 循環(huán)水槽示意圖Fig.1 Diagrammatic sketch of circulating water channel

數(shù)值模擬中使用的幾何模型與實驗尺寸一致如圖2 所示,水槽長度方向為X方向,寬度方向為Y方向,深度方向為Z方向.選取水槽上覆水體中間剖面作為研究對象.

圖2 幾何模型Fig.2 Geometric model

2.2 邊界條件和初始條件

流場方程計算中,水槽壁面為無滑移邊界

流場方程計算中,頂部表面為對稱邊界,法向速度為0 和剪切應(yīng)力為0,是Dirichlet 條件和Neumann條件的組合,其中n為邊界的外法線方向

顆粒和濃度輸運方程計算中,水槽壁面和頂面為無通量邊界條件

左右邊界在流場和輸運方程計算中均采用周期性邊界條件,即等速度,等污染物濃度,等顆粒體積分數(shù),相等的湍流動能和相等的湍流耗散率.

初始速度設(shè)定為0,為靜水狀態(tài),初始壓力與水深相關(guān),其中p0為大氣壓力

污染底泥顆粒粒徑D50=0.03 mm,含水量為57.5%,上覆水體中顆粒體積分數(shù)為0.污染底泥中污染物相對濃度為1,上覆水體中相對濃度設(shè)定為0.

上覆水和底泥的密度ρ1和ρ2分別為997 kg/m3和2650 kg/m3.上覆水的黏度設(shè)定為1.0×10?3Pa·s.

經(jīng)過網(wǎng)格獨立性測試之后,數(shù)值計算采用82 800個單元的結(jié)構(gòu)化網(wǎng)格,并在泥?水界面和底邊界處加密,泥?水界面加密層數(shù)為16 層,底邊界加密層數(shù)為8 層,如圖3 所示.

圖3 界面網(wǎng)格示意圖Fig.3 Diagrammatic sketch of interface grid

2.3 數(shù)值方法

本文采用有限元方法[26-28]求解流場和污染物濃度場.首先利用Galerkin 有限元法建立連續(xù)模型方程在空間上的離散形式.該方法是基于微分方程的弱形式,將動量方程(2) 乘以權(quán)函數(shù)w∈H1(?)3(H1(?)3∈H1(?),H1(?)是Hilbertian-Sobolev 空間)并在域? 上積分得到[29],如下

式中,速度場的解u∈H1(?)3,權(quán)函數(shù)w∈H1(?)3.

基于P1D-P2 單元方法[29],速度場求解采用不連續(xù)的分段線性形函數(shù)(P1D形函數(shù)),壓力求解采用連續(xù)的分段二次形函數(shù)(P2 形函數(shù)),該單元方法具有理想的Ladyzhenskaya-Babuska-Brezzi (LBB) 穩(wěn)定性.速度場的形函數(shù)和解都是在子空間H1(?)3∈H1(?)上獨立計算的,對于任意一個子空間(單元),速度場的形函數(shù)和解可以由該空間上的形函數(shù)ψ(P1D形函數(shù))的線性組合表示

式中,N是單元節(jié)點數(shù),wj是節(jié)點j處權(quán)函數(shù)的值,uk是節(jié)點k處的速度.ψj,ψk分別是該單元節(jié)點j,k處的形函數(shù),且僅在j,k處有值,其他節(jié)點上為0.

壓力p由連續(xù)分段二次形函數(shù)(P2 形函數(shù))表示

式中,pl是節(jié)點l處的壓力值,Ψl分別是節(jié)點l處的形函數(shù),且僅在l處有值,其他節(jié)點上為0.

由此,可以給出一個N×N的線性方程組離散動量方程,并求解速度場u和壓力p.

式中,M,A,K和C分別是質(zhì)量、對流、應(yīng)力和梯度矩陣.Fleft是方程左側(cè)的阻力項矩陣.b和fleft分別是浮力和方程右邊的阻力項.

同樣,連續(xù)方程也可以通過有限元方法離散如下

式中,Ci為散度算子,r是應(yīng)用Dirichlet 邊界條件的曲面積分項.

方程離散后,數(shù)值計算采用壓力投影法(pressure projection method)[31-32].該方法根據(jù)Helmholtz-Hodge分解定理,引入一個中間速度場和壓力場,對速度和壓力解耦后分步計算,步驟如下:

(1) 使用最新測試解計算近似速度場,un+1=θutent+(1 ?θ)un,θ=0.5.

(2)求解壓力Poisson 方程,得到測試壓力場ptent.

(3)求解離散后的動量方程得到速度場u?.該速度場u?并不能滿足連續(xù)方程,從而得到修正壓力?p.

(4)用?p修正u?得到速度場測試解utent.

(5)將測試解utent代入流場作為下一次迭代的最新測試解.

當(dāng)Picard 迭代達到極限或收斂,基于此時的utent和ptent可得到最終的un+1和pn+1.重復(fù)上述求解方法,直到達到所需的求解時間或達到穩(wěn)定狀態(tài).

3 模型驗證

本文通過計算不同流速條件下的上覆水中非吸附性污染物濃度變化,并比較計算結(jié)果與水槽實驗結(jié)果,驗證數(shù)值模型的準確性.實驗時,緩慢加水調(diào)節(jié)至實驗設(shè)計的流速和水深,通過分析上覆水體中代表性物質(zhì)的相對濃度C/C0,得到底泥再懸浮中物質(zhì)釋放的特征.模型驗證中,計算4 組不同流速條件下,非吸附污染物濃度隨時間變化,如圖4 所示,4 種流速0.13 m/s,0.23 m/s,0.35 m/s 和0.50 m/s 分別對應(yīng)計算結(jié)果1,2,3 和4.由于再懸浮過程中污染物釋放在短時間內(nèi)完成,受制于實驗條件,底泥再懸浮污染物濃度測量工況有限.本文選取流速為0.15 m/s 時,上覆水體中污染物濃度隨時間變化,如圖4 所示.同時選取靜水條件下,上覆水體污染物濃度隨時間變化.可以看出,受到水動力學(xué)條件影響,上覆水體中污染物濃度在短時間內(nèi)(< 2 min)達到平衡,數(shù)值模擬結(jié)果和實驗測量結(jié)果基本吻合.相對于靜水條件下,污染物從底泥孔隙水中通過分子擴散向上覆水體釋放,底泥再懸浮污染物釋放量在短時間內(nèi)可達靜水時的5 倍.

圖4 數(shù)值計算結(jié)果和實驗結(jié)果比較Fig.4 Comparison of numerical results and experimental results

4 結(jié)果和討論

底泥顆粒隨流場的運動非常復(fù)雜,具有很大的隨機性.為了定量分析,選取幾何模型中上覆水體中心剖面作為研究對象,引入平均速度、相對平均濃度和平均湍流動能的概念如下

式中,C為相對平均濃度,為無量綱值,U為平均速度(m/s),K為平均湍動能(m2/s2).

選取幾何模型中上覆水體中心剖面作為研究對象,底泥起動過程中上覆水體平均速度隨時間的變化如圖5 所示.底泥起動初期,不同工況的水體平均速度均很小,隨時間延長,水體平均速度迅速增大并在短時間內(nèi)達到穩(wěn)定.不同工況下,穩(wěn)定后流速值如表1 所示.

圖5 底泥再懸浮后速度隨時間變化Fig.5 The variation of the velocity with time after sediment resuspended

表1 不同工況下平均流速Table 1 Average velocity under different cases

選取Case3 和Case4 數(shù)據(jù),分析不同時間底泥顆粒體積分數(shù)φ 在水深方向的分布,如圖6.不同流速條件下,φ 的垂向分布均隨時間快速平衡,并均在2 min 內(nèi)達到同一個穩(wěn)定值,即在短時間內(nèi)上覆水體中φ 達到均勻,這與實驗現(xiàn)象接近.同一時刻,流速越快上覆水體中φ 越均勻.對于細顆粒底泥D50=30 μm,其空間結(jié)構(gòu)分布均勻,顆粒向上的再懸浮和向下的沉降運動之間的平衡導(dǎo)致了顆粒體積分數(shù)分布的平衡,當(dāng)水動力條件不變時,這種平衡不會打破.

圖6 不同時間上覆水體顆粒體積分數(shù)垂向分布Fig.6 The Vertical distribution of particle volume fraction in overlying water at different time

一般情況下,污染物在水體中擴散受到對流作用、分子擴散、湍流作用、吸附解吸等眾多物理化學(xué)機制的影響.選取Case1,Case3,Case4 和Case5 的數(shù)據(jù),分析不同流速時,污染物平均相對濃度C、平均湍動能K與時間的關(guān)系,如圖 7所示.在不同流速下,上覆水體中C均在短時間內(nèi)(< 3 min)達到平衡.這一現(xiàn)象和φ 短時間內(nèi)在上覆水體達到平衡存在對應(yīng)關(guān)系.不同流速下,隨底泥起動上覆水中K均迅速增大,流速越快,K增大越快,C達到平衡時間越短.從Case1 到Case5,流速提高了約11.6 倍,K的峰值也增加超過11 倍.當(dāng)水體中φ達到平衡后,K隨之趨于平衡,C也趨于平衡.在底泥起動過程中,K和水體中φ 變化密切相關(guān),進入水體的細顆粒底泥改變了上覆水的流場結(jié)構(gòu),增大了上覆水的K.K的增大降低了上覆水體中C達到平衡的時間.迅速進入上覆水體的細顆粒底泥,改變了上覆水體流場特性,進而影響到污染物的釋放.對于動水條件下底泥中非吸附性污染物的釋放,水動力學(xué)條件是其擴散的關(guān)鍵.

圖7 C–K–t 關(guān)系曲線Fig.7 The curve of C–K–t

進一步研究不同流速條件下,對流作用和湍流作用對污染物釋放的貢獻.通過分析污染物相對平均濃度隨流速和時間變化,得到污染物擴散總通量NT,對流擴散通量NC,湍流擴散通量ND和流場特性(雷諾數(shù)Re) 之間的關(guān)系,如圖8.對于非吸附性污染物,NT隨流速線性增大,也即隨著雷諾數(shù)線性增大.當(dāng)雷諾數(shù)較小時(Re< 9 000),污染物擴散過程中,對流和湍流擴散貢獻基本一致,湍流擴散貢獻略大于對流擴散的貢獻.當(dāng)雷諾數(shù)增大后(9 000 15 000),湍流擴散貢獻迅速下降,污染物擴散主要由對流作用主導(dǎo).可以認為隨著流速增大,底泥再懸浮非吸附性污染物釋放逐漸由對流擴散主導(dǎo).通常情況下,富含污染物的細顆粒底泥總是存在于流速較低的河流、湖庫等水域中,受到自然或人為擾動后,細顆粒底泥很容易懸浮,水體的流場特性隨之改變,此時湍流擴散和對流擴散貢獻接近,湍流擴散量不可輕易忽略.

圖8 污染物釋放通量與雷諾數(shù)關(guān)系曲線Fig.8 The curve of release flux and Reynolds number

5 結(jié)論

底泥再懸浮污染物釋放過程是上覆水體?底泥?污染物耦合的釋放過程.底泥起動后復(fù)雜的流場特性是底泥再懸浮釋放污染物的主要影響因素.細顆粒底泥易受到水流作用的影響,由于其空間結(jié)構(gòu)相對均勻,沉降速度較小,再懸浮時水深方向的泥沙顆粒在很短時間內(nèi)(< 2 min)就能達到均勻,并保持穩(wěn)定.迅速進入上覆水體的底泥顆粒,改變了上覆水體流動特性,進而影響到污染物的釋放.平均湍動能在底泥細顆粒泥沙起動后達到峰值,隨垂向泥沙顆粒分布均勻而降低,并隨時間逐漸達到穩(wěn)定.流速增加,湍動能隨之增加,上覆水體中污染物濃度達到平衡的時間也越短.對流作用和湍流擴散作用均對上覆水體污染物濃度達到平衡的過程有影響.對于非吸附介質(zhì),當(dāng)雷諾數(shù)較小時,污染物擴散過程中,對流和湍流擴散貢獻基本一致,湍流擴散貢獻略大于對流擴散的貢獻.隨雷諾數(shù)增大,對流擴散貢獻逐漸大于湍流擴散的貢獻.當(dāng)雷諾數(shù)較大時,湍流擴散貢獻迅速下降,污染物擴散主要由對流作用主導(dǎo).通常情況下,富含污染物的細顆粒底泥總是存在于流速較低的河流、湖庫等水域中,受到自然或人為擾動后,細顆粒底泥很容易懸浮,此時湍流擴散和對流擴散對污染物釋放的貢獻接近,湍流擴散量不可輕易忽略.