基于格子玻爾茲曼方法的植被水流數(shù)值模擬研究

趙 芳，白鳳朋，馮振鵬，姜家良，劉 標(biāo)

(1.中冶南方城市建設(shè)工程技術(shù)有限公司，湖北 武漢 430077；2.長江水資源保護(hù)科學(xué)研究所，湖北 武漢 430051)

水生植被是河道生態(tài)系統(tǒng)的重要組成部分[1]，植被的存在增加了河床的阻力，改變了河道中的水流結(jié)構(gòu)，使河道水流流速重新分布，水流結(jié)構(gòu)的調(diào)整改變了河流生物棲息地的特征，會(huì)進(jìn)一步影響河流生態(tài)系統(tǒng)的演化。植被與水體的交互作用機(jī)理一直是環(huán)境水力學(xué)領(lǐng)域的研究熱點(diǎn)，水生植被修復(fù)也是河流生態(tài)修復(fù)的重要組成部分[2-3]。

近二十年來，格子Boltzmann方法吸引了眾多的研究學(xué)者的目光，開始有研究學(xué)者將這一新的方法引入計(jì)算水力學(xué)領(lǐng)域[4-5]。在植被水流研究中，數(shù)值模型可以彌補(bǔ)室內(nèi)水槽和野外試驗(yàn)的缺點(diǎn)，較為準(zhǔn)確地獲得整個(gè)計(jì)算區(qū)域的水動(dòng)力特征。Jiménez等將植被的影響概化為拖曳力，采用二維Lattice Boltzmann方法模擬剛性植被水流[6]；Gac將植被當(dāng)作固體邊界，基于大渦模擬技術(shù)建立了模擬植被水流的三維格子Boltzmann數(shù)學(xué)模型[7]。

本文基于格子Boltzmann方法建立了擬植被影響的二維淺水流動(dòng)數(shù)學(xué)模型，并與基于Godunov型有限體積方法建立的二維淺水?dāng)?shù)學(xué)模型的計(jì)算結(jié)果進(jìn)行了對(duì)比，測試和對(duì)比結(jié)果表明，建立的格子Boltzmann數(shù)值模型可以成功模擬植被影響的水流流場調(diào)整再平衡過程。

1 控制方程

采取張量表達(dá)方式，水深平均的二維淺水方程組形式如下：

(1)

(2)

式中，h—水深；i、j—笛卡爾坐標(biāo)系空間指標(biāo)，并且采用了愛因斯坦求和約定;ui—i方向上的流速：i=x，ui=u；i=y，ui=v;vf—運(yùn)動(dòng)粘性系數(shù)；Fi—源項(xiàng)，包括地形項(xiàng)、摩擦阻力項(xiàng)以及植被引起的阻力項(xiàng)：

(3)

2 剛性植被處理方法

對(duì)于生長在河床上的剛性(挺水)植被，往往占有一定的過水面積，其對(duì)水流結(jié)構(gòu)有顯著的影響，流速分布形式如圖1所示。

圖1 剛性植被垂向流速分布

常用的基于拖曳力系數(shù)的拖曳力繞流阻力表達(dá)式為：

(4)

在x、y方向上的剛性植被引起的拖曳力表達(dá)式分別為：

(5)

式中，Cd—拖曳力系數(shù)，λ—單位體積水體中植被的擋水面積，UC—植被層的水深平均流速，uc、vc—x、y方向上的植被層的水深平均流速。

在非淹沒狀態(tài)下U=UC；在淹沒狀態(tài)下，采用Stone and Shen(2002)提出的計(jì)算公式：

(6)

式中，η—流速校正系數(shù)，約為1.0；hv—?jiǎng)傂灾仓甑母叨龋籙—整個(gè)水體的水深平均流速。

單位體積水體中植被的擋水面積λ的計(jì)算公式為：

(7)

式中，αv—植被形狀系數(shù)，當(dāng)植株為圓柱體時(shí)，αv為1.0；Vd—在植被區(qū)域植被體積占整個(gè)水體的百分比；dv—植株圓柱體直徑。

3 格子Boltzmann方法

格子Boltzmann模型如下：

fα(x+eαΔt，t+Δt)=fα(x，t)-

(8)

選擇D2Q9格子模型，1—8粒子沿著特定的格子鏈方向以特定的速度遷移到各自相鄰的節(jié)點(diǎn)；0粒子速度為0，留在中心節(jié)點(diǎn)。D2Q9模型粒子速度矢量eα如圖2所示。

圖2 D2Q9格子類型

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

(9)

(10)

宏觀變量水深h、流速ui與粒子分布函數(shù)的關(guān)系如下：

(11)

4 模型的驗(yàn)證與應(yīng)用

4.1 植被覆蓋一側(cè)的矩形河道

試驗(yàn)在一傾斜循環(huán)式水槽中進(jìn)行，水槽全長25.5m，寬1m，高1m，底坡為0.0005，水槽無植被時(shí)曼寧系數(shù)是0.01。數(shù)值計(jì)算2個(gè)試驗(yàn)工況，分別記為工況1和工況2。工況1：水深為0.0457m，平均流速為0.32m/s，植被間距0.028m，植被密度c為0.0023；工況2：水深為0.0428m，平均流速為0.276m/s，植被間距0.02m，植被密度c為0.0044。數(shù)值模型相關(guān)計(jì)算參數(shù)：Δx=Δy=0.02m，Δt=0.005s，τ=0.55，Smagorinsky常數(shù)Cs=0.2，Cd=1.5，二次流附加阻力影響系數(shù)k=0.2。

將格子Boltzmann數(shù)學(xué)模型的計(jì)算結(jié)果同采用有限體積方法計(jì)算結(jié)果進(jìn)行對(duì)比，如圖3—4所示。格子Boltzmann數(shù)學(xué)模型的計(jì)算結(jié)果與試驗(yàn)測量結(jié)果吻合較好，模型可以成功模擬植被對(duì)水流結(jié)構(gòu)的調(diào)整作用，例如：植被區(qū)流速減小，自由水流區(qū)流速變大，植被區(qū)與自由水流區(qū)存在著較大的流速梯度等。格子Boltzmann數(shù)學(xué)模型的計(jì)算結(jié)果同采用有限體積方法的計(jì)算結(jié)果相差無幾。

圖3 工況1斷面平均流速計(jì)算值與實(shí)測值對(duì)比

圖4 工況2斷面平均流速計(jì)算值與實(shí)測值對(duì)比

4.2 植被對(duì)撫河故道河道水流運(yùn)動(dòng)影響

撫河故道位于撫河下游西岸南昌縣黃馬鄉(xiāng)、三江鎮(zhèn)境內(nèi)，是撫河防洪非工程措施的重要組成部分，起于撫河下游左岸支流箭江口堵口處，止于崗前大壩，是贛撫平原西總干渠的一段，全長約18.0km。撫河故道上游來水主要由2部分組成：一是自焦石攔河大壩引入撫河水，經(jīng)西總干渠流入撫河故道；二是當(dāng)撫河干流箭江口以上河段出現(xiàn)較大洪水時(shí)，箭江分洪閘執(zhí)行分洪功能，分泄部分撫河洪水進(jìn)入撫河故道。植被主要分布在河道上游灘地和河心洲。正常流量下植被區(qū)域往往不過水；當(dāng)撫河故道執(zhí)行分洪功能時(shí)，水流會(huì)流入植被區(qū)，對(duì)河道的行洪能力產(chǎn)生一定的影響。

本文主要研究撫河故道分洪時(shí)植被對(duì)水流的影響，對(duì)比有植被與無植被水位和流速的變化。計(jì)算區(qū)域選擇植被分布集中的上游河段，約為7km長，植被分布情況如圖5所示，每個(gè)區(qū)域植被的種類和相關(guān)參數(shù)見表1。

圖5 計(jì)算區(qū)域植被分布圖

表1 計(jì)算區(qū)域植被參數(shù)

上游來水量汛期實(shí)測為65m3/s，下游出口水位為實(shí)測21.76m，水流處在恒定流時(shí)考慮箭江分洪閘從撫河分洪。根據(jù)箭江分洪閘管理規(guī)定，設(shè)定分洪流量為200m3/s，并且為持續(xù)分洪。在分洪時(shí)有植被與無植被工況下的水位和流速變化情況如圖6—7所示。由圖6可以看出，受到邊灘植被和箭江分洪閘至主河道之間植被的阻水影響，河道上游水位壅高，水位抬升明顯，水位最大增幅達(dá)到0.6m，增加了上游洪水的威脅；下游河道沒有植被分布，同時(shí)受到出口邊界的影響，2種工況的水位相差無幾；水位抬升作用自上游至下游逐漸減弱。圖7表明有植被的主河道的流速小于無植被工況，流速增加的地方主要位于河道邊灘。這是由于邊灘上的植被增加了河床的粗糙度，引起水位上升和河道過水面積增加，在無植被工況下不能過水的部分邊灘開始過水。

圖6 水位變化圖(有植被工況減去無植被工況)

圖7 流速變化圖(有植被工況減去無植被工況)

5 結(jié)論

本文以淺水動(dòng)力學(xué)和基于格子Boltzmann方法理論知識(shí)為基礎(chǔ)，建立了模擬淺水流動(dòng)的二維格子Boltzmann數(shù)值模型。運(yùn)用室內(nèi)水槽試驗(yàn)對(duì)模型進(jìn)行了驗(yàn)證，并與基于Godunov型有限體積方法建立的二維淺水?dāng)?shù)學(xué)模型的計(jì)算結(jié)果進(jìn)行了對(duì)比，驗(yàn)證結(jié)果表明建立的格子Boltzmann數(shù)值模型可以成功模擬植被水流流場分布和植被引起的水流調(diào)整再穩(wěn)定過程，利用建立的模型揭示了植被對(duì)撫河故道行洪時(shí)水位、流速的影響規(guī)律。本文在如何描述柔性植被對(duì)水流的影響、模型并行計(jì)算等方面存在不足。