振蕩潛流帶沉積層-水界面污染物輸運的研究1)

陳金峰 張金龍 楊文武 董宇紅

(上海大學力學與工程科學學院,上海市應用數學和力學研究所,上海 200072)

引言

在工業化和城市化快速發展過程中,污染廢水和生活用水在排放到江河湖海之前會經過無害化處理[1].當外部污染源得到有效控制后,沉積層成為更突出的潛在污染源.相關研究表明,我國水體潛流帶沉積物的污染率高達80.1%[2].當水流遇到擾動或水體的環境條件(如酸堿度、溶解氧含量等)發生改變時,沉積層中的各種污染物如重金屬離子會再次釋放到上覆水體中,造成二次污染[3-4].沉積層-水交界面(sediment-water interface,SWI)是沉積層與上覆水體間進行物理、化學、生態作用的主要區域之一.上覆水的流動特性和沉積層的物質和結構屬性很大程度上影響跨越SWI 的物質輸運過程[5].研究河口海岸流動和沉積層中污染物的輸運機理和傳質規律,以及研究潛流帶流場統計特性、相干結構和猝發事件與物質輸運的關系是科學界、工程界共同關注的課題,對于認識和保護水資源生態具有重要的科學意義.

在河口、湖泊、海岸帶等地區,潮汐、內波、風荷載、水泵作業等多種因素會產生振蕩流動現象[6].不同因素所產生的振蕩流中,其振幅和頻率差異性很大,導致物質輸送中的動力學過程更為復雜.因此,有必要研究振蕩流條件下的潛流交換對沉積層內外污染物釋放與輸運的影響.

以往研究表明,在低速平緩流動和沉積物不發生懸浮的情況下,SWI 處通常存在厚度以毫米計的擴散邊界層,分子擴散是控制這一區域物質垂向輸運的主要機制[7].但是當上覆流動復雜化,如湍流狀態或振蕩流時,對SWI 區域的混合和輸運過程需要加以深入研究[8-9].Spalart 等[10]數值研究了受擾動的層流和間歇湍流模式,并分析了壁面應力和流速分布.Reidenbach 等[11]采用大渦模擬方法研究了高Reynolds 數下的振蕩湍流邊界層及轉捩問題,進一步拓展了此類問題的研究.

隨后,一些學者將研究聚焦到流體混合和傳質方面,涉及到的溶質包括溶解氧、營養鹽等.Lorke等[12]在瑞士的阿爾卑納赫湖測量了周期為18 小時的振蕩流中溶解氧垂向分布情況.數據顯示大幅振蕩的流動導致了SWI 區域質量通量的不穩定性.Higashino 等[13]通過數值模擬發現溶解氧濃度和Sherwood 數的波動幅度與振蕩周期、振幅和Schmidt 數有關.Tian 等[14]研究了非恒定的驅動條件下高Schmidt 數傳質過程.研究發現,低頻高振幅的振蕩壓力梯度驅動力對流動和物質輸運有顯著影響.Wang 等[15]進一步研究了振蕩流下的氧濃度垂向分布、SWI 區域的湍流通量和沉積層耗氧量隨時間的變化規律.Thomas 等[16]研究發現表面波誘導的振蕩流動增強了底棲生物對養分的吸收速率.雖然通過實驗和實地測量總結了一些質量輸運經驗公式可以預測水體沉積層溶質的吸收釋放率[17].

以上研究主要關注潛流交換中“自上而下”的物質輸送過程,而累積于沉積層內的污染物質則是在受到擾動后“自下而上”的釋放和擴散過程.Cheng 等[18]通過水槽實驗測量了沉積層內的磷向水體中釋放過程中的濃度變化,分析了平均流速、周期和振幅等水動力因素對其的影響.注意到以往對潛流帶中污染溶質輸送的研究多著眼于其分布特點和污染水平[19],較少涉及復雜水動力條件和底床物理特性(如滲透性、孔隙率等).基于以上分析,本文運用大渦模擬方法對具有沉積層的振蕩槽道流和其中高Schmidt 數污染物傳質問題展開研究.通過求解三維不可壓縮黏性Navier-Stokes (N-S)方程和對流擴散方程來模擬沉積層內外的流動和傳質過程,同時采用修正的Darcy-Brinkman-Forcheimer 模型描述有鋅離子污染溶質的可滲透沉積層.研究聚焦于振蕩流和界面上下物質輸運過程的耦合作用,探究不同振蕩條件對動量輸運和污染物溶質輸運的影響,以及由此產生的沉積污染物的釋放和輸送規律與機制.

1 數理模型與數值方法

1.1 物理模型與控制方程

在研究含有沉積層的污染溶質輸運中,湍流運動特性與跨SWI 的質量通量之間依賴關系是一個重點,這對實驗測量有相當的難度,而模型化和數值模擬則具有優勢.本文物理模型如圖1 所示,x,y,z分別代表流向、法向和展向,物理模型的計算域為Lx×Ly×Lz=2πh×1.1h×1.5πh,對應的網格數目為Nx×Ny×Nz=128×196×128.沉積層深度為0.1h,該局部區域法向網格數為35.流向和展向為均勻網格,法向網格采用伸縮變換,適當增加壁面、自由面和SWI 處的網格分辨率,網格間距為 Δx+=25.26,Δz+=18.95,Δy+=0.43～5.61 .槽道中水層高度h取10 cm作為特征長度,SWI 處的摩擦速度uτ為特征速度,沉積層內初始濃度C0=2.34 mg/L 為特征濃度.由此,無量綱的三維控制方程為

圖1 底部含有沉積層上部為自由面的槽道振蕩流動模型示意圖Fig.1 The sketch of the oscillatory flow model of a channel with a sediment layer at the bottom and a free surface at the top

式中上劃線“-”表示濾波后的變量.流向x、法向y、展向z的速度分量分別為u,v和w或ui(i=1,2,3),C為鋅離子濃度.摩擦Reynolds 數定義為Reτ=uτh/υ,Schmidt 數Sc=ν/Dw,描述同時存在動量擴散及分子擴散的無量綱參數.ν是水的運動黏性系數,Dw是金屬鋅離子在水中的分子擴散系數.水中鋅離子可被視為被動標量,其分子擴散系數相對于運動黏度小兩個量級,是一種高Schmidt 數的傳質現象.p表示濾波后的壓力.振蕩驅動力Pgδik(k=1,3)與流向形成一個夾角A,簡稱振蕩角;T為振蕩周期,該驅動力表示為

τij和qj分別為亞格子尺度湍流應力項和湍流質量通量項,這里采用的是Germano 等[20]提出的動力學亞格子模型.假定亞格子尺度湍流應力與大尺度運動應變率有如下關系,渦黏性系數進一步假設為可以得到如下數學表達式

式中S cT=νT/DT是借鑒分子輸運模型定義的湍流Schmidt 數,方程式(6)和式(7)中的G和S cT將用動力學方法求出[14,21].

一般來說,沉積層具有復雜的幾何結構和廣泛的長度尺度特征,這些性質與多孔介質相符.Whitaker[22]提出在一個體積為V,半徑為r的球型控制體上平均N-S 方程,只對多孔介質中的大尺度流動進行建模從而降低對網格的要求.這樣的體積平均方程被稱為VANS(volume-averaged Navier-Stokes equations)方程.對于一般的物理量 φ ,〈φ〉s表示相體平均,〈φ〉f表示固有體平均,相體平均被定義為

其中Vf是控制體中流體所占的體積,固有體平均被定義為

兩者之間的關系為

于是得到在沉積層中的流體的控制方程如下

當上覆水為較高Reynolds 數流動,具有多孔滲透性的下墊沉積層間隙流速與水力梯度并非線性關系.因此本文采用修正的Darcy-Brinkman-Forcheimer 非線性模型來反映高滲透沉積層中的阻力影響[23],即

其中達西數Dai=Ki/h2表示滲透率Ki的影響,ε 為孔隙率,常數a=150和b=1.75 .

在計算中,沿槽道流向x和展向z取周期性邊界條件,下壁面y=0 處取無滑移無穿透速度邊界條件;在自由面y=1.1 處,不考慮變形的影響,取剪切力為零的條件,即

初始時刻,沉積層內污染物濃度為1,水中濃度為0;邊界條件為下壁面濃度恒定為1,上方自由面處取濃度的法向通量為0.速度和濃度場在SWI 界面處采用階躍界面條件[24-25],即

設定上述界面條件,從物理上合理給出了界面的動量傳遞關系,是將VANS 方程(控制多孔介質內部流動) 與N-S 方程(控制自由流動流體) 在SWI 處耦合所需的應力條件.

1.2 數值方法和驗證

本文的時間離散方法為Verzicco 等[26]在Kim 等[27]分裂步方法基礎上改進得到的混合格式[28].對流項的時間推進用三階Runge-Kunta 法,黏性項采用Crank-Nicholson 隱式格式,空間導數離散采用二階中心差分.控制方程的離散都是基于交錯網格[29].

為了驗證程序計算的可靠性和合理性,本文將采用大渦模擬,首先對帶自由面的高Reynolds 數槽道流這一典型物理模型進行驗證計算,并將結果與Moser 等[30]的直接數值模擬(direct numerical simulation,DNS)結果對比.為此,與該DNS 模擬的物理和計算參數取值一致,即Reτ=590,計算域為Lx×Ly×Lz=2πh×h×πh.本文采用的網格為128×196×128,前述DNS 網格數為 384×257×384 .圖2 是平均流向速度沿法向分布,壁面坐標定義為y+=yuτ/ν .

圖2 流向平均速度分布Fig.2 Profiles of the mean streamwise velocity

從圖中可以看出,平均速度與Moser 等[30]的DNS 結果吻合良好.在黏性底層,平均速度符合u=y+的線性律.當y+＞40,平均速度剖面符合對數律的典型分布,為u=2.5lny++5.5 .有研究表明,當沉積層中的滲透率愈來愈小時,流動特性與規范槽道流趨于一致(Ref.[31]),結果對比說明本文數值計算方法以及程序能夠準確合理地預測高Reynolds數壁面剪切湍流特性.

然后,考慮到本文著重研究具有底部沉積層的高Schmidt 數傳質問題,除了將上述的LES 模擬的流場統計結果與DNS 結果進行比較外,還進行了另一個數值驗證,以便與O'Connor 等[32]的含有沉積層水槽流動實驗數據進行比較,特別是考察對高Schmidt 數標量(溶解氧)輸運和其濃度分布的模擬準確性.表1 顯示了實驗的相關物理參數,包括沉積層厚度h、流體黏性系數、摩擦Reynolds 數、基準濃度、Schmidt 數.這里計算與實驗保持一致,取S c=373,Reτ=180 .達西數和孔隙率表示沉積層的滲透性,計算取值參照實驗[32],即Da=1×10-4,ε=0.7.計算域、網格數和數值格式同上.圖3 為模擬所得平均濃度在垂向的分布,同時給出了與實驗數據的比較.圖中SWI 表示水-沉積層交界面.將SWI 的位置平移到y=0.08,規范模擬算例與實驗之間的坐標對應關系.圖中可見計算的平均濃度分布與實驗結果吻合良好.根據上述比較和驗證,可以證實本文計算方法和程序能夠準確預測帶有沉積層的槽道流高Schmidt 數傳質問題.

表1 流場物理參數(h=10 cm)Table 1 Physical parameters of flow field (h=10 cm)

圖3 平均濃度在垂向的分布Fig.3 Profiles of the mean concentration along the vertical direction

2 結果與討論

2.1 振蕩流動下濃度場與流場的相關性

在底部為高滲透沉積層的三維槽道振蕩流物質輸運問題中,剪切流動主導其傳質過程.隨機運動的流體微團和相干運動的大尺度結構為共同載體促進了高Schmidt 數標量輸運行為.通過湍流各尺度之間,以及流動作用于交界面來實現其中動量、能量和物質的混合和交換.本文算例的模擬參數均設置在水溫35 °C,摩擦Reynolds 數Reτ=515,Schmidt數S c=109,孔隙率ε取0.45,達西數Da為1.0 ×10-4.

表2 給出不同振蕩角A和振蕩周期T下的物理參數和統計量,〈〉 表示統計量在10 個周期內的時均值.水層和沉積層中的體積平均濃度分別定義為

表2 不同算例的計算參數和湍流通量等統計量Table 2 Calculation parameters and statistics such as turbulent flux for different cases

〈C′v′〉SWI和 〈J〉SWI分別為交界面處的濃度通量和分子擴散通量.

首先分析了SWI 上覆振蕩流的湍流結構、速度、能量和濃度變化之間的相關性.圖4 為振蕩驅動力Pg與體積平均的湍動能Ek、流向速度Ub、展向速度Wb隨時間的變化曲線.Ek的表達式為

體積平均的流向速度Ub、展向速度Wb均由摩擦速度uτ無量綱化得到.當流向振蕩驅動力占優時,即A為15°,Ub顯著大于Wb.隨著振蕩周期增大,Ek,Ub,Wb的波動幅度逐漸增大,如圖4(a)和圖4(b)所示.當振蕩驅動力逐漸過渡到展向為主時,即A從45°再到75°時,湍動能Ek的幅值變化并不明顯,但Ub與Wb振蕩幅值相當(圖4(c)),并逐漸發展出以展向速度Wb大幅振蕩為特征的流動,如圖4(d)所示.

圖4 振蕩驅動力 Pg與體積平均的展向速度 Wb、流向速度 Ub、湍動能 Ek的關系Fig.4 Oscillation driving force versus volume-averaged spanwise velocity Wb,streamwise velocity Ub,turbulent kinetic energy Ek

圖5 進一步給出沉積層內體積平均的濃度Cs、流向速度Ub和湍動能Ek隨時間變化曲線.由于濃度場對流場的延遲響應,可以發現Cs與Ek和Ub之間存在明顯的相位差.Ek呈現增長與衰弱的交替變化,并作用于傳質過程.前者在增強過程中帶動Cs持續增高,反之亦然,表明振蕩中的湍流強度與高滲透沉積層中污染物的準周期擴散呈現明顯關聯性.當流向振蕩驅動力占優時(A= 15°),高頻振蕩的湍流運動抑制了沉積層污染物濃度Cs的波動性(圖5(a));而低頻振蕩使得Cs的波動幅值逐漸增大(圖5(b)).當展向振蕩驅動力占優時(A= 75°),雖然Ub和Ek的波動幅度明顯低于A= 15°的情況,但平均濃度的波動仍較為顯著(圖5(a)和圖5(b)),這表明低頻振蕩流與沉積層內污染物濃度波動的關聯性更強.

圖5 體積平均的流向速度、湍動能和沉積層內濃度隨時間變化曲線Fig.5 Volume-averaged streamwise velocity,turbulent kinetic energy,and concentration in sediment versus time

隨時間變化的振蕩驅動力不僅對平均速度場和濃度場產生重要的影響,而且在其作用下,中小尺度湍流脈動和濃度脈動的統計量也會交替地增強或減弱,呈現準周期的變化.圖6 給出了不同振蕩角和頻率下相位平均的法向湍流強度分布.可以看到法向速度脈動增強與減弱的程度和相位均有所不同.低頻振蕩對于流動有顯著影響,其法向脈動峰值更高,波動更強烈.

為了觀察振蕩驅動作用對瞬時流場結構的影響,圖7 給出了A=45°,T=12 算例在SWI 上方y+=103處一個周期內不同相位的瞬時法向速度的水平等值線圖.圖中黑色實線和黑色虛線分別是法向湍流脈動速度v′=2.0和v′=-2.0 的等值線.在相位t/T=5/8,6/8 時,從圖7(b)和圖7(c)中可以看到流場中比較豐富的正負斑團結構.而如圖7(a)和圖7(d)所示,在相位t/T=4/8,7/8 時,因法向脈動較弱而斑團結構稀少,這與圖6(c)所顯示的相位平均湍流強度變化相對應.結合標注相位點的圖7(e),可以觀察到在第2 和3 點處,即t/T=5/8,6/8 處湍動能達到高位,此時法向速度脈動也達到峰值,如圖7(b)和圖7(c)所示,隨后湍流強度開始減弱.法向湍流強度的這一周期性變化將會強化湍流混合,使得其質量通量在法向發生間歇性的變化,進而促進其中的物質輸送和交換.

圖6 相位平均的法向湍流強度Fig.6 Profiles of phase-averaged wall-normal turbulent intensities

相應地,圖8 給出了A=45°,T=12 時在不同相位的法向濃度通量瞬時分布(z方向視圖).法向湍流濃度通量分為兩部分,藍色部分(C′v′=0.01)為正向通量,棕色(C′v′=-0.01)為負向通量.通過圖8(b)和圖8(c)發現C′v′=0.01 的發生概率更大,說明正向濃度通量即向上輸運起主導作用.同時濃度通量在不同相位變化顯著,其演化狀態與法向湍流強度具有一致性.當t/T=4/8,5/8 時,法向速度脈動的增強伴隨著濃度通量的增大.當t/T=6/8,濃度通量大幅度地增加,隨后在t/T=7/8 相位點減弱.說明法向湍流脈動在污染物溶質從沉積層往上覆水釋放過程中發揮重要作用.綜合圖7 可以判斷在相位t/T=4/8,5/8 擴散子區成為更強的傳質限制層,但在相位t/T=6/8,7/8 部分含能渦能夠突破限制主導質量輸送.

圖7 A=45°,T=12算例在x-z平面 y+=103 處法向速度脈動的等值線圖(v′=±2.0)Fig.7 Contour plot of the instantaneous wall-normal velocities(v′=±2.0) parallel to the x-zplane (y+=103) for the A=45°,T=12case

圖8 算例 A=45°,T=12 在不同相位的法向濃度通量瞬時等值面分布圖(z方向視圖).藍色和棕色分別代表法向濃度通量值為C′v′=0.01和C′v′=-0.01Fig.8 Contour maps of the instantaneous wall-normal convective concentration flux in different phases for A=45°,T=12 (z-direction view).The blue and the brown iso-surface means C′v′=0.01 and C′v′=-0.01

2.2 振蕩角和振蕩周期對平均速度和濃度的影響

從上覆水流場到其下高滲透沉積層的高Schmidt數傳質行為受制于近SWI 界面處湍流的主要特性,包括黏性底層和擴散子區的流動特性.圖9(a)給出流向平均速度沿法向的分布.需要注意的是各圖中出現的豎直細虛線表示沉積層和水層交界面SWI(y=0.1或y+=52).當展向振蕩驅動力占優時,即A=75°,長周期振蕩導致水層流向速度明顯減小.而當流向振蕩驅動力占優時(A=15°),則流向速度均維持較高水平.

圖9(b)給出平均濃度沿法向分布.可以看出,當振蕩周期較大時,無論A=15°還是A=75°,沉積層中污染物釋放速率加快,使得上覆水中的濃度水平更高,表2 中的〈Cw〉也反映了相同的結果.原因在于SWI 處速度梯度增大,該處剪切應力和摩阻增大.說明振蕩驅動力主導下有利于突破擴散子層限制增強在沉積層上下的動量和質量交換,導致污染溶質更易于向上釋放和擴散.

圖9 (a)流向平均速度和(b)平均濃度沿法向分布Fig.9 (a) Mean streamwise velocity and (b) mean concentration distribution along the normal direction

進一步,采用象限分析法考察了振蕩過程中如湍流猝發事件在動量傳遞和污染物輸運過程中所起的作用.圖10 統計了雷諾應力的變化.其中圖10(a)和圖10(b)分別是第二象限的上拋運動(u′＜0,v′＞0)和第四象限的下掃運動(u′＞0,v′＜0).可以看出,當振蕩角減小或周期增大時,SWI 之上流動的上拋下掃行為隨之增強,且上拋強于下掃.結合前面圖7 和圖8 的分析,可知污染物上揚事件與猝發事件的關聯性.這再次表明了沉積層中污染物釋放速率加快,導致上覆水中高濃度水平的內在原因.

圖10 雷諾應力在第二、四象限的分布Fig.10 Distribution of Reynolds stress in the second and fourth quadrants

圖11(a)給出了不同振蕩條件下湍流誘導的濃度通量沿法向分布.可以發現,隨著周期增大,濃度通量增大.當周期較小時(T=3),濃度通量對振蕩角的變化不敏感.而當周期較大時,隨著振蕩角減小,法向通量也隨之增大.值得注意的是,SWI 上方的濃度通量并非線性遞減,而是會在y+=150 處出現再次增高現象,這源自于該處法向速度脈動強度最大,如圖11(b)的y=0.3 處.以較低頻率進行的流向或展向振蕩會顯著提高法向濃度通量,促進污染物向上覆水中遷移擴散.

圖11 (a)濃度通量和(b)法向湍流強度沿法向分布Fig.11 Distribution of (a) concentration flux and (b) turbulent intensities along normal direction

在自然水環境系統中,由多種不同尺度且相互作用的水動力機制共同驅動著SWI 區域的物質交換,這些機制包括分子擴散、對流擴散和湍流滲透,而有效擴散系數能夠綜合反映多種機制的共同作用[33].從表2 中可以得知,隨著振蕩周期增大,當流向驅動力占優時,分子擴散通量增大;而當展向驅動力占優時,分子擴散通量減小,說明其濃度輸運的增強主要來自于湍流作用.值得注意的是,當T=12時,由于流向速度出現了間歇性的回流(圖4(b)的流向速度Ub出現負值),導致SWI 上方產生較大的湍流通量.

當滲透率較高時,湍流通量是在質量輸運中占主導地位的,而隨著周期增大,SWI 處的湍流通量增大.在此定義有效擴散率為分子擴散通量與湍流通量之和與分子擴散通量的比值,表達式為

圖12 統計了不同振蕩角和振蕩周期下的有效擴散系數的分布情況,圖中橫坐標為振蕩角,縱坐標為SWI 處的有效擴散率.從整體趨勢來看,隨著振蕩角或振蕩周期增大,擴散率隨之增大.而當振蕩驅動力逐漸過渡到以展向為主時,較長的振蕩周期對擴散率的增幅影響更大.

圖12 不同振蕩角和振蕩周期與有效擴散系數的關系Fig.12 Relationship between different oscillation angles and oscillation periods and effective diffusivity

3 結論

本文通過大渦模擬研究了有振蕩激勵下沉積層內外湍流流動特性、相干結構、猝發現象、統計特性和污染物輸運特性.對不同振蕩角與振蕩周期下的平均湍流強度、污染物溶質濃度分布、通量變化和瞬態結構進行了分析,得到結論如下.

(1)周期性振蕩驅動作用導致了體積平均的流向、展向速度、湍動能以及沉積層內污染物濃度呈現與振蕩角和振蕩周期直接相關的準周期性變化,且與振蕩驅動力存在不同的相位差,表明速度場和濃度場對振蕩驅動力的變化有著不同程度的時間延遲響應.

(2)當展向振蕩驅動力強于流向振蕩時,隨著振蕩周期增加,平均湍動能以及沉積層內污染物濃度波動幅度增大.由于SWI 處流體的上拋強于下掃,進而以對流擴散方式作用于濃度標量.此處有效擴散率增大,污染物溶質向水層中輸運增強,表明展向振蕩流與污染物濃度輸運的波動的關聯性更強.

(3) SWI 界面的濃度通量變化與法向速度脈動具有明確的相關性.較低頻率和展向振蕩為主的雙重條件會強化湍流混合,顯著提高有效擴散率.研究還證實SWI 處污染物上揚事件與湍流猝發事件的一致性,表明猝發行為在上揚過程中有效推動其中的物質輸送和交換.