基于格子Boltzmann的植被對河道的流場和污染物輸運的影響研究

黃德法,楊佳剛,歐陽競一,高 軒

(1.中國電建集團華東勘測設(shè)計研究院有限公司,浙江 杭州 311122;2.湖北省水利水電規(guī)劃勘測設(shè)計院有限公司,湖北 武漢 430070;3.四川省眉山市東坡區(qū)人民政府辦公室,四川 眉山 620010;4.西北農(nóng)林科技大學(xué) 水利與建筑工程學(xué)院,陜西 楊凌 712100)

近年來,河流開發(fā)的程度越來越大,造成了河流生態(tài)系統(tǒng)破壞,出現(xiàn)了越來越多的河流污染問題,河流的修復(fù)和保護(hù)成為了人們越來越關(guān)心的問題[1],而植被在河流生態(tài)治理中有十分積極的作用,其對水流的物理阻擋作用會改變水流阻力和水流結(jié)構(gòu)[2-3],減輕了水流對河岸的沖刷作用,而水流結(jié)構(gòu)的變化會影響河道中顆粒物、污染物、營養(yǎng)物質(zhì)等溶質(zhì)的輸運,從而對水質(zhì)和河床形態(tài)造成影響。但自然界河流中植被的分布方式往往十分復(fù)雜,研究復(fù)雜植被的分布方式對河道流場和污染物擴散輸移的影響非常必要。

一些學(xué)者通過在水渠中移栽植被或者利用人工模擬植被的實驗方法,研究了植被對明渠水流的影響[4-5]。但隨著計算機技術(shù)的不斷發(fā)展,數(shù)值模擬方法由于其方便、快捷以及能模擬較復(fù)雜工況的優(yōu)點,逐漸地成為了研究含植被水流以及污染物輸移擴散的重要方法[6-7]?；荻郲8]分別將植被的作用概化為繞流拖曳阻力和曼寧底部摩阻,建立了水深平均的二維水流數(shù)值模型。魯俊等[9]采用分布算法以及有限差分求解控制方程,對含植被水渠中水流的流速和切應(yīng)力分布進(jìn)行了數(shù)值模擬。白鳳朋等[10]基于Godunov型有限體積法,將植被的作用視為一個拖曳力,建立了一個模擬淺水植被水流的二維數(shù)學(xué)模型?；蔽男诺萚11]則將植被的作用等效為河床的附加阻力,推導(dǎo)出了植被區(qū)域的等效曼寧系數(shù)以模擬植被水流的流速分布。

以上求解淺水方程以及對流擴散方程的傳統(tǒng)數(shù)值模擬方法主要包括有限體積法、有限差分法、有限元法等,這些方法具有良好的模擬效果,但在效率或精度方面仍有待改進(jìn)。除了這些傳統(tǒng)數(shù)值模擬方法外,格子Boltzmann方法(Lattice Boltzmann Method,LBM)在近年來成為了計算流體力學(xué)領(lǐng)域中比較重要的方法之一。與傳統(tǒng)方法基于連續(xù)介質(zhì)假定描述宏觀流體運動不同,LBM方法是一種基于介觀模擬尺度的計算流體力學(xué)方法,具有介于微觀分子動力學(xué)模型和宏觀連續(xù)模型的特點[12-13]。其具有易編程、計算效率高、并行化計算和便于處理復(fù)雜邊界等優(yōu)勢。這些優(yōu)勢使得可以將格子Boltzmann方法作為處理淺水流動和對流擴散的一個重要工具[14-17],但目前相關(guān)研究相對較少。

本文使用基于離散Boltzmann方法的D2Q16多速度淺水模型和基于格子Boltzmann方法的D2Q5對流擴散模型來研究植被間隔覆蓋情況下植被的不同分布條件對流場和污染物輸運的影響。

1 河道流場及污染物輸運模型

使用基于離散Boltzmann方法的D2Q16多速度淺水模型和基于格子Boltzmann方法的D2Q5模型,計算含植被河道的流場和污染物濃度場分布以及研究復(fù)雜分布植被對河道流場和污染物輸運的影響。模型計算步驟如圖1。

圖1 模型計算步驟

1.1 D2Q16多速度淺水模型

單松弛因子的Boltzmann淺水控制方程為[18-19]:

(1)

(2)

式中:f為粒子的分布函數(shù)f(r,c,t);F為作用力加速度;feq為平衡分布函數(shù);τ為松弛因子;g為重力加速度;V表示水深方向上的平均速度;h為水深;分布函數(shù)f(r,c,t)的物理意義為t時刻以r=(x,y)為中心的單位面積內(nèi),以速度c=(cx,cy)為中心的單位粒子速度區(qū)間內(nèi)的總粒子數(shù)量。宏觀變量可以通過對分布函數(shù)在粒子速度區(qū)間內(nèi)積分得到。

當(dāng)使用離散玻爾茲曼模型時,為方便計算,采用無量綱變量表示,即:

(3)

式中:頂部帶小角符號的表示無量綱化后的變量;r為物理長度;h0為數(shù)值模擬中的特征水深;Q表示流量;v為運動粘性系數(shù);v0為參考運動粘性系數(shù);P0為參考壓強;P為壓強。

因此可以將式(1)和式(2)轉(zhuǎn)換為無量綱的形式,同時引用一個“Knudsen”數(shù)如下:

(4)

弗勞德數(shù)則變?yōu)?

(5)

引入無量綱弛豫時間τ′:

(6)

式(1)變?yōu)槿缦滦问?

(7)

通過Hermite展開和Gauss型積分對粒子速度區(qū)間進(jìn)行離散,得到離散Boltzmann多速度淺水模型的控制方程[18]:

(8)

粒子i的平衡分布函數(shù)的表達(dá)式為:

(9)

式中:ωα為權(quán)重。

ωα權(quán)重值如下:

(10)

式中:α為粒子速度方向的序號。

在模型中,D2Q16的格子模型如圖2所示。

圖2 D2Q16格子模型

通過對整個粒子速度空間上的分布函數(shù)進(jìn)行積分,可以獲得宏觀量水深和速度,表達(dá)式如下:

h=∑fα

(11)

(12)

1.2 D2Q5對流擴散模型

通過基于格子Boltzmann方法的D2Q5模型求解對流擴散方程以用于模擬污染物的擴散輸移[20-21]。單松弛因子的Boltzmann對流擴散模型的控制方程可以表示為[22]:

(13)

eα表達(dá)形式為:

(14)

式中:ex和ey為x和y方向上的格子速度。

ex、ey分別表示為:

ex=Δx/Δt

(15)

ey=Δy/Δt

(16)

式中:Δx和Δy分別為x和y方向上的格子大小。

平衡分布函數(shù)的表達(dá)式為:

(17)

式中:λx、λy分別為x、y方向上無量綱后的物理擴散系數(shù),表達(dá)式如式(18);c′為污染物濃度。

(18)

在模型中,D2Q5的格子模型如圖3所示。

圖3 D2Q5格子模型

沿深度方向上的平均濃度值c′可以通過式(19)獲得:

(19)

2 應(yīng)用案例

自然界河流中植被分布往往十分復(fù)雜,同時,在一些城市的河道中,布置著一些人工生態(tài)植被用于凈化水質(zhì)、恢復(fù)生態(tài),其分布方式有:沿河道兩側(cè)以不同密度分布、沿河道兩側(cè)間隔分布和在河道中呈斑狀分布等[17]。其中沿河道兩側(cè)間隔分布是常見的分布方式之一,如圖4所示。

圖4 河道兩岸植被間隔分布

參考圖4中自然河流以及人工布置生態(tài)植被的河道中兩側(cè)植被區(qū)間隔分布的情況,利用D2Q16多速度淺水模型和D2Q5對流擴散模型模擬研究植被區(qū)間隔布置的情況下水流的流速分布以及污染物擴散的濃度分布,并研究植被區(qū)間隔的距離對流場和污染物輸運的影響。假定一段長度為5.8 m、寬度為1.95 m的河道,沿程坡降為1‰,自由流動無植被區(qū)域的床面粗糙度曼寧系數(shù)為0.01,入口處控制流量為28.6 L/s。植被直徑為0.5 cm,密度設(shè)為124 株/m2,即植被區(qū)的粗糙度nv設(shè)為0.06。設(shè)置4種研究工況,如表1,植被區(qū)布置分布圖見圖5,橫斷面分為CS-1、CS-2、CS-3三個斷面,分別代表第一段植被出口斷面、第二段植被入口斷面和第二段植被出口斷面,縱斷面分為LS-1、LS-2、LS-3三個斷面,對應(yīng)植被區(qū)、過渡區(qū)和中心無植被區(qū),K-1、K-2、K-3分別表示各個區(qū)域內(nèi)的污染物投入點的起始位置,投入點的水深為平均水深。其中,植被粗糙度nv計算公式為[23]:

表1 植被區(qū)間隔分布研究各工況參數(shù)

圖5 植被區(qū)間隔布置情況數(shù)值模擬模型圖

(20)

式中:k表示含植被渠道中因植被作用所引起的一種阻力系數(shù),根據(jù)經(jīng)驗取值,對于部分植被覆蓋的矩形渠道,k取0～0.3;d表示植被區(qū)的植被密度占比;nb為河床的曼寧摩擦系數(shù);Cd表示植被拖拽阻力系數(shù),Cd的經(jīng)驗取值范圍為0.8～3.5[24];αv為形狀系數(shù);Dv為植被要素的直徑。

模型的參數(shù)設(shè)置上,采用0.01 m×0.01 m大小的網(wǎng)格,即在長度和寬度方向分別離散為581和196個網(wǎng)格節(jié)點,dt=0.005 s,特征水深h0=0.029 m,弛豫時間τ=1.5,無量綱化后的污染物擴散系數(shù)λ=0.016。上下邊壁采用周期邊界,入口和出口處使用非平衡外推邊界類型,水深通過插值一點得到,出口處控制水深使平均水深為4.4 cm,流速通過插值一點得到。

3 結(jié)果與討論

3.1 間隔植被覆蓋情況流場模擬結(jié)果

使用工況1—工況4模擬植被區(qū)間隔分布條件下的流速沿程變化規(guī)律、橫斷面流速分布規(guī)律以及橫向流速的分布和變化規(guī)律。

據(jù)圖6(a)展示的是植被區(qū)(LS-1)的流速沿程變化規(guī)律?？梢钥吹礁鱾€工況下,這個區(qū)域的沿程流速在進(jìn)入第一段植被后都會有下降趨勢,離開第一段植被后速度逐漸升高,進(jìn)入第二段植被后再次下降。工況1情況下,在離開第一段植被后速度回升不明顯,是因為該工況下兩個植被區(qū)間隔較小,水流在離開第一個植被段后馬上就會受到第二個植被段的影響。圖6(b)展示的是過渡區(qū)(LS-2)的流速沿程變化規(guī)律。在該區(qū)域,流速變化十分不均勻,轉(zhuǎn)點非常明顯,這些轉(zhuǎn)點即為植被段和無植被段的交界處,不同工況下該區(qū)域的沿程流速變化規(guī)律大致相同,只是速度變化的位置發(fā)生了改變。圖6(c)展示的是中心無植被區(qū)(LS-3)的流速沿程變化規(guī)律。各個工況下,該區(qū)域的沿程流速在進(jìn)入第一個植被段的斷面后都會有較大的升高趨勢,離開第一個植被段的斷面后速度升高趨勢逐漸平緩,進(jìn)入第二個植被段的斷面后趨勢再次上升。通過三個區(qū)域的結(jié)果可以看到,植被區(qū)的間隔距離對于三個區(qū)域內(nèi)流速的漲幅影響不大。此外,流速沿程增加,說明還沒有形成均勻流。

圖6 縱斷面流速沿程變化

圖7展示的分別是CS-1、CS-2和CS-3三個橫斷面的4種工況下的流速分布規(guī)律。

圖7 橫斷面流速分布

由于4種工況下植被布置的位置不同,每個工況的橫斷面設(shè)置有所差異,但水流流速分布均呈中間大兩邊小的趨勢,即植被區(qū)流速整體較小,中心無植被區(qū)流速較大,過渡區(qū)存在較大的流速梯度,選取受植被影響充分的CS-3斷面,橫斷面流速分布仍然可以表達(dá)為式(21)形式,但系數(shù)變?yōu)榕c植被區(qū)間隔D有關(guān)。

(21)

(22)

通過建立植被段末端斷面的植被區(qū)、過渡區(qū)和中心無植被區(qū)流速與植被區(qū)間隔關(guān)系的公式,可以得到植被段末端斷面不同植被間隔分布下的植被區(qū)和中心無植被區(qū)的流速,以及過渡區(qū)不同位置處的流速。

圖8展示的是4種工況下的橫向流速場的總體分布。當(dāng)水流在經(jīng)過植被段時,受植被的影響會產(chǎn)生方向指向河道中心的橫向流速,在兩個植被區(qū)中間間隔的部分,水流又會以較小的流速流向兩邊邊壁,工況1中由于兩個植被區(qū)的間隔很小,所以在間隔的區(qū)域只出現(xiàn)了一小部分區(qū)域的水流流向邊壁?？傮w來看,橫向流速基本上是以y=1 m為中心線對稱的,但是也有局部地方不對稱,原因是本文中采用算法的誤差造成的。

圖8 橫向流速場總體分布

圖9展示的是植被區(qū)(LS-1)和過渡區(qū)(LS-2)的橫向流速沿程變化。各工況下,橫向流速均有兩個波峰,兩個波峰所在的斷面大致為植被段所在斷面。對比不同工況,橫向流速的最大值還與植被區(qū)的間隔距離有關(guān),在植被區(qū)和過渡區(qū)中,該規(guī)律一致,即植被區(qū)間隔距離越大,前一個植被區(qū)植被導(dǎo)致的水流橫向流速越小,而后一個植被區(qū)植被導(dǎo)致的水流橫向流速越大。橫向流速的峰值大小見表2,可以由表中數(shù)據(jù)得到,植被區(qū)中橫向流速兩個波峰的峰值大小Vmax1和Vmax2與植被區(qū)間隔D基本呈線性關(guān)系,如圖10,可概括為式(23),過渡區(qū)中滿足的關(guān)系也基本一致。

表2 橫向流速峰值大小及其出現(xiàn)位置 單位:m/s

圖9 橫向流速沿程變化

圖10 橫向流速最大值與植被間隔的關(guān)系

(23)

3.2 間隔植被覆蓋情況污染物輸運過程模擬結(jié)果

對4種工況下河道的污染物濃度場分布進(jìn)行模擬。圖11展示的分別是4種工況下從K-1處釋放污染物后整體的濃度場分布情況,圖12展示的是4種工況下從K-2處釋放污染物后整體的濃度場分布情況。圖13為K-1點位處釋放污染物時該點位縱截面的污染物濃度值。K-3處釋放的污染物由于受到兩側(cè)植被的對稱影響,仍然呈沿程不斷減小的順直條狀分布,這里不再展示。在每個工況下,污染物在沒有進(jìn)入植被段時都以直線向下游輸移,遇到植被后會發(fā)生橫向移動,同一個工況下,在河道較為邊緣處(K-1)釋放污染物的橫向移動的發(fā)生位置要比在過渡區(qū)前方(K-2)釋放污染物時靠前。在離開植被段后,污染物又會沿著直線往下游輸移。污染物濃度沿程不斷減小。

圖11 K-1處釋放污染物的濃度場分布

圖12 K-2處釋放污染物的濃度場分布

圖13 K-1處釋放污染物時縱截面污染物濃度值

接下來分析植被區(qū)間隔距離對于污染物濃度最終衰減大小以及污染物最終偏移距離的影響。據(jù)表3和圖14,可以看到從K-1釋放污染物時,4種工況下出口斷面污染物濃度衰減情況與植被間隔D呈線性關(guān)系,從K-2釋放污染物時,4種工況下出口斷面污染物濃度衰減情況與植被間隔D呈二次關(guān)系。

表3 出口斷面污染物濃度衰減大小和橫向偏移距離

圖14 污染物濃度衰減與植被間隔的關(guān)系

其關(guān)系可以表達(dá)為:

(24)

從K-1釋放污染物時,在出口斷面濃度最大值處離濃度投放點的橫向偏移距離在4種工況下均為0.03 m,從K-2釋放污染物時,在出口斷面濃度最大值處離濃度投放點的橫向偏移距離在4種工況下則均為0.06 m。說明植被區(qū)分布的間隔距離對于污染物橫向偏移的距離基本沒有影響。

4 結(jié) 論

本文基于格子Boltzmann方法模擬計算了植被區(qū)間隔布置條件下河道的流場和污染物輸運濃度場分布情況,主要得到以下結(jié)論:

(1) 在流速沿程變化中,兩個植被區(qū)之間的間隔中流速會出現(xiàn)迅速回升,過渡區(qū)的水流流速變化十分不均勻,植被段和無植被段的交界處的轉(zhuǎn)點明顯。

(2) 植被區(qū)和中心區(qū)流速呈直線分布,過渡區(qū)流速與植被密度呈二次冪相關(guān),植被的存在會使流線出現(xiàn)明顯的偏轉(zhuǎn),出現(xiàn)橫向流速,橫向流速的最大值與植被密度呈線性正相關(guān),橫向流速中出現(xiàn)的兩個波峰峰值分別與植被區(qū)間隔距離呈負(fù)相關(guān)性和正相關(guān)性,可概括為線性關(guān)系。

(3) 污染物呈一個帶有兩段明顯橫向偏移的條狀分布,污染物衰減幅度與植被區(qū)間隔呈正相關(guān)性,在植被區(qū)為線性關(guān)系,在過渡區(qū)為二次關(guān)系。