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

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

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

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

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

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

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

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

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

圖1 模型計算步驟

1.1 D2Q16多速度淺水模型

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

(1)

(2)

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

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

(3)

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

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

(4)

弗勞德數則變為:

(5)

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

(6)

式(1)變為如下形式:

(7)

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

(8)

粒子i的平衡分布函數的表達式為:

(9)

式中:ωα為權重。

ωα權重值如下:

(10)

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

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

圖2 D2Q16格子模型

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

h=∑fα

(11)

(12)

1.2 D2Q5對流擴散模型

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

(13)

eα表達形式為:

(14)

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

ex、ey分別表示為:

ex=Δx/Δt

(15)

ey=Δy/Δt

(16)

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

平衡分布函數的表達式為:

(17)

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

(18)

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

圖3 D2Q5格子模型

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

(19)

2 應用案例

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

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

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

表1 植被區間隔分布研究各工況參數

圖5 植被區間隔布置情況數值模擬模型圖

(20)

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

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

3 結果與討論

3.1 間隔植被覆蓋情況流場模擬結果

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

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

圖6 縱斷面流速沿程變化

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

圖7 橫斷面流速分布

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

(21)

(22)

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

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

圖8 橫向流速場總體分布

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

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

圖9 橫向流速沿程變化

圖10 橫向流速最大值與植被間隔的關系

(23)

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

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

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

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

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

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

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

圖14 污染物濃度衰減與植被間隔的關系

其關系可以表達為:

(24)

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

4 結 論

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

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

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

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