劉寧寧，時(shí) 勇，孫華林

（山東省海河淮河小清河流域水利管理服務(wù)中心，濟(jì)南 250100）

0 引 言

洪水演進(jìn)是指根據(jù)河段上斷面的洪水過程推求河段下斷面未來的洪水過程，其理論依據(jù)是洪水波在河道中的傳播規(guī)律。天然河道中的洪水波屬于淺水長波，它的運(yùn)動規(guī)律可以采用圣維南方程組進(jìn)行描述，但圣維南方程組屬于雙曲型擬線性的偏微分方程，直到目前還不能給出這類方程的解析解，只能用數(shù)值方法近似求解，不僅計(jì)算復(fù)雜，而且需要詳細(xì)的河道地形資料。因此，常在實(shí)踐中采用近似的計(jì)算方法，即采用圣維南方程組的簡化形式。圣維南方程組的簡化形式多采用運(yùn)動波方程和擴(kuò)散波方程。經(jīng)研究發(fā)現(xiàn)，擴(kuò)散波不考慮慣性項(xiàng)但考慮附加比降項(xiàng)，能反映回水影響，能用于平原河網(wǎng)地區(qū)，其適用范圍相較運(yùn)動波更廣［1-3］。

Lattice Boltzmann 方法是20世紀(jì)80年代末期提出的一種模擬流體流動［4］的數(shù)值計(jì)算方法，是建立在介觀動理學(xué)模型基礎(chǔ)上，以粒子速度分布函數(shù)為研究對象，建立微觀粒子與宏觀量之間的聯(lián)系。它打破了傳統(tǒng)數(shù)值計(jì)算方法的理論框架，不再受連續(xù)介質(zhì)條件的束縛，并在求解非線性偏微分方程組時(shí)，將非線性偏微分方程組轉(zhuǎn)化為Lattice Boltzmann 演進(jìn)方程，從而變?yōu)楹唵蔚木€性方程，使得計(jì)算效率得到大大提高。這些特點(diǎn)促使Lattice Boltzmann 方法在模擬大規(guī)模流動和求解非線性方程中得到了廣泛的應(yīng)用［5-8］。在水力學(xué)方面：程永光等［9］采用Lattice Boltzmann 方法模擬了一維、二維明渠非恒定流，并模擬了實(shí)際工程中二維圓形潰壩波的傳播；張東輝等［10］通過采用Lattice Boltzmann 方法求解Richards 方程，對土壤下滲過程進(jìn)行了模擬；吳家陽［11］提出了基于Lattice Boltzmann 方法的河庫瞬變流模擬方法，并與GPU 結(jié)合進(jìn)行并行計(jì)算，使得計(jì)算效率大大提升。

在水文學(xué)中，所謂“物質(zhì)特性”的傳輸，表面上是流速、流量、濃度的改變，但其物理本質(zhì)卻是質(zhì)量、動量、能量的傳輸，從微觀角度來理解，實(shí)質(zhì)上就是組成宏觀量的粒子的碰撞和遷移，因此，Lattice Boltzmann 方法應(yīng)用于水文學(xué)中，具有較強(qiáng)的物理基礎(chǔ)。擴(kuò)散波方程是修正的Burgers 方程參數(shù)取特定值時(shí)的特例，是其在水文學(xué)中的應(yīng)用，近年來，Lattice Boltzmann 方法在求解Burgers 方程方面有較多的研究［12-14］。本文將Lattice Boltzmann方法運(yùn)用到擴(kuò)散波方程中，通過算例進(jìn)行對比分析。

1 擴(kuò)散波方程

河道中的洪水波運(yùn)動可以采用擴(kuò)散波方程進(jìn)行描述，當(dāng)河段下邊界為不受回水頂托影響的自由下邊界，初始流量為零，且不考慮旁側(cè)入流時(shí)，可用下列定解問題描述：

式中：Q為流量；Cd為擴(kuò)散波的波速；μ為擴(kuò)散波的擴(kuò)散系數(shù)；x為河長；t為時(shí)間；I（t）為河段上斷面的入流過程；L為河段長度。

2 Lattice Boltzmann 模型

2.1 Lattice Boltzmann方程

Lattice Boltzmann 方法是將空間劃分為均勻的網(wǎng)格，并在每個(gè)節(jié)點(diǎn)上設(shè)置粒子分布函數(shù)，根據(jù)粒子分布函數(shù)在網(wǎng)格上的運(yùn)動特性，碰撞項(xiàng)采用BGK 近似［15］，得到Lattice BGK方程：

式中：fα為流體粒子的分布函數(shù)；Ωα為碰撞項(xiàng)；τ為弛豫時(shí)間；為平衡態(tài)分布函數(shù)；α為方向，與選取的離散速度模型有關(guān)。

Lattice Boltzmann 方法認(rèn)為在Δt時(shí)間內(nèi)粒子能以恒定的速度e→沿α方向從一個(gè)節(jié)點(diǎn)遷移到下一個(gè)相鄰節(jié)點(diǎn)，并在這個(gè)節(jié)點(diǎn)上與原有粒子發(fā)生碰撞，取代原有分布函數(shù)，其碰撞和遷移過程分別由式（3）和（4）確定：

碰撞過程：

在碰撞過程中，碰撞項(xiàng)需要滿足質(zhì)量和動量守恒條件：

2.2 擴(kuò)散波方程的Lattice Boltzmann 解法

構(gòu)造Lattice Boltzmann 模型的關(guān)鍵在于選擇合適的平衡態(tài)分布函數(shù)，而平衡態(tài)分布函數(shù)的具體表達(dá)式又與離散速度模型的構(gòu)造有關(guān)［16］。Qian 等人提出的DnQb模型是最具有代表性的離散速度模型，其中，Dn表示空間維數(shù)，n可取1、2、3；Qb表示離散的速度數(shù)，b可取3、5、7、9、15、19 等［6］。本文嘗試采用D1Q5速度模型，以期提高計(jì)算精度。D1Q5 為一維五速模型，即在每個(gè)網(wǎng)格點(diǎn)上存在5 類不同速度的粒子，粒子速度=c[0，-1，-2，1，2］，α=0，1，2，3，4，c為網(wǎng)格步長Δx與時(shí)間步長Δt之比，即：c=ΔxΔt。同時(shí)，采用對空間坐標(biāo)不進(jìn)行多尺度處理，而對時(shí)間坐標(biāo)引入五個(gè)時(shí)間尺度的多尺度展開方式［10］，得到：

式中：ε為Knudsen 數(shù)，表示粒子平均自由路徑與宏觀流動中特征長度的比值。

并將Chapman-Enskog展開式應(yīng)用于分布函數(shù)，得到：

對式（2）的第一項(xiàng)進(jìn)行二階Taylor展開，并將式（6）和（7）代入，最終可得到五個(gè)不同時(shí)間尺度的離散Boltzmann方程：

由于在Chapman-Enskog 展開中假設(shè)局部范圍內(nèi)已接近平衡態(tài)，故對式（7）兩邊求和，得到：

為建立D1Q5 模型的平衡態(tài)分布函數(shù)與擴(kuò)散波方程中流量這個(gè)宏觀量之間的聯(lián)系，并檢驗(yàn)多尺度展開方式選取的正確性，需將五個(gè)不同時(shí)間尺度的離散Boltzmann 方程通過求解各階矩的方式還原到式（1）。具體步驟詳見文獻(xiàn)［17］可得：

將式(9)+ε×式(10)+ε2×式(11)+ε3×式(12)五個(gè)方向求和得：

與宏觀方程式（1）相比，當(dāng)k=時(shí)，式（15）只增加了一項(xiàng)o(ε4)，可以通過τ來控制。這樣就得到了具有四階截?cái)嗾`差的擴(kuò)散波方程。D1Q5 模型對應(yīng)的擴(kuò)散波方程的平衡態(tài)分布函數(shù)由式（14）求解，得到：

根據(jù)宏觀初始流量，由式（16）計(jì)算各節(jié)點(diǎn)平衡態(tài)分布函數(shù)，通過粒子的碰撞和遷移過程，更新每個(gè)節(jié)點(diǎn)的分布函數(shù)，并通過式（14）即可得到不同時(shí)刻各節(jié)點(diǎn)的宏觀流量。

3 應(yīng)用實(shí)例

此次采用兩個(gè)算例，分別取自文獻(xiàn)［18］和［19］。算例一：長江上游龍街—喬家河段，河段長247 km，平均河床坡度0.001 12。1967年6月18日至20日，在河段上游和下游兩端的兩個(gè)站點(diǎn)測量了一次洪水事件。根據(jù)龍街和巧家兩水文站的有關(guān)資料得到波速Cd=5.92 m/s，擴(kuò)散系數(shù)μ=19 456 m2/s。算例二：湖南資水柘溪水電站下游江南—夫溪河段，河段長20.1 km，1967年5月13日至14日，發(fā)生一次洪水過程。根據(jù)江南和夫溪兩水文站的有關(guān)資料得到波速Cd=2.792 m/s，擴(kuò)散系數(shù)μ=8 376 m2/s。

3.1 模擬結(jié)果

采用Muskingum 法、Lattice Boltzmann 方法的D1Q5 模型以及解析解法對所給算例進(jìn)行求解。算例一：對于Muskingum法，為了保證長河段條件下的線性條件，需將長河段分成4 個(gè)子河段即N=4，作分段連續(xù)演算，子河段xe=0.45，Ke=Δt=3 h。對于解析解法，為保證計(jì)算精度，劃分矩形條塊時(shí)，取計(jì)算時(shí)段Δt=3 h。對于Lattice Boltzmann 方法的D1Q5 模型，取弛豫時(shí)間τ=1.5，空間步長Δx=1 000 m，時(shí)間步長Δt=10 s。算例二：對于Muskingum 法，將長河段分成4 個(gè)子河段即N=4，作分段連續(xù)演算，子河段xe=0，Ke=Δt=0.5 h。對于解析解法，取計(jì)算時(shí)段Δt=0.5 h。對于Lattice Boltzmann 方法的D1Q5 模型，取弛豫時(shí)間τ=1.5，空間步長Δx=402 m，時(shí)間步長Δt=5 s。上述3 種解法所得計(jì)算結(jié)果及精度分析見表1和圖1。

由表1和圖1可以看出，這3種方法計(jì)算得到的洪峰相對誤差、峰現(xiàn)時(shí)間、確定性系數(shù)基本一致，這表明Lattice Boltzmann方法可以應(yīng)用于河道洪水波的演算。

圖1 3種方法計(jì)算與實(shí)測流量過程比較Fig.1 The comparisons of observed and computed hydrograph

表1 3種方法計(jì)算精度統(tǒng)計(jì)表Tab.1 The comparisons of observed and computed values

3.2 影響因素分析

Lattice Boltzmann 方法的計(jì)算精度主要受到弛豫時(shí)間、步長選取、速度模型、邊界格式等因素的制約［10］。在水文預(yù)報(bào)中，判斷預(yù)報(bào)精度的指標(biāo)有4 個(gè)：①洪峰誤差；②徑流量誤差；③峰現(xiàn)時(shí)間；④確定性系數(shù)。洪峰、徑流量的誤差在20%以內(nèi)，峰現(xiàn)時(shí)間在一個(gè)計(jì)算時(shí)長之內(nèi)，確定性系數(shù)在0.7 以上，即認(rèn)為是合格的預(yù)報(bào)，對于誤差的分析，4 個(gè)指標(biāo)應(yīng)綜合考慮。經(jīng)計(jì)算，步長和弛豫時(shí)間的選取對峰現(xiàn)時(shí)間沒有影響，徑流量誤差可以由確定性系數(shù)反應(yīng)，故以洪峰誤差和確定性系數(shù)作為分析對象，將Lattice Boltzmann方法的計(jì)算結(jié)果與實(shí)測流量進(jìn)行比較。

3.2.1 步長對計(jì)算精度的影響

首先分析空間步長對計(jì)算精度的影響。選擇弛豫時(shí)間τ=2，時(shí)間步長Δt=5s，綜合考慮河長、計(jì)算效率、計(jì)算穩(wěn)定性等因素，對于算例一，分別取空間步長Δx=500 m、Δx=1 000 m、Δx=2 470 m、Δx=4 940 m、Δx=12 350 m 進(jìn)行計(jì)算；對于算例二，分別取空間步長Δx=402 m、Δx=804 m、Δx=1 005 m、Δx=1 340 m、Δx=2 010 m進(jìn)行計(jì)算，計(jì)算結(jié)果詳見表2。

表2 不同空間步長對計(jì)算精度的影響Tab.2 The influence of different space step on calculation accuracy

其次分析時(shí)間步長對計(jì)算精度的影響。選擇弛豫時(shí)間τ=2。綜合考慮河長、計(jì)算效率、計(jì)算穩(wěn)定性等因素，對于算例一，在空間步長Δx=2 470 m 的情況下，分別取時(shí)間步長Δt=1 s、Δt=2 s、Δt=5 s、Δt=10 s、Δt=20 s進(jìn)行計(jì)算；對于算例二，在空間步長Δx=1 005 m 的情況下，分別取時(shí)間步長Δt=1 s、Δt=2 s、Δt=5 s、Δt=10 s、Δt=20 s進(jìn)行計(jì)算，計(jì)算結(jié)果詳見表3。

由表2～表3可知，時(shí)間步長比空間步長對計(jì)算精度的影響要小一些。這是因?yàn)?，空間步長的選取直接影響到相同河段長情況下，網(wǎng)格劃分的多少，進(jìn)而影響粒子的碰撞遷移次數(shù)。因此，在采用Lattice Boltzmann 方法進(jìn)行洪水演算時(shí)，應(yīng)結(jié)合上下游斷面長度，選取合適的空間步長。同時(shí)，在計(jì)算效率允許的情況下，也應(yīng)綜合考慮粒子速度（Δx/Δt）的取值，因?yàn)榱Ｗ铀俣龋é/Δt）取值太小，會導(dǎo)致計(jì)算不穩(wěn)定。

表3 不同時(shí)間步長對計(jì)算精度的影響Tab.3 The influence of different time step on calculation accuracy

3.2.2 弛豫時(shí)間對計(jì)算精度的影響

弛豫時(shí)間τ是一個(gè)無綱量參數(shù)，表示粒子分布函數(shù)趨于平衡態(tài)的時(shí)間。在理論上，為了滿足Lattice Boltzmann 方程的線性穩(wěn)定性，弛豫時(shí)間τ的取值應(yīng)大于0.5［20］。但從計(jì)算精度的角度出發(fā)，對于擴(kuò)散波方程的Lattice Boltzmann 方法D1Q5 速度模型，應(yīng)分析弛豫時(shí)間τ在哪個(gè)取值范圍，可以使計(jì)算誤差最小。為討論弛豫時(shí)間τ的取值對計(jì)算精度的影響，應(yīng)在特定的空間和時(shí)間步長情況下，選取不同的弛豫時(shí)間τ進(jìn)行計(jì)算分析。對于算例一，選取空間步長Δx=1 000 m，時(shí)間步長Δt=5 s；對于算例二，選取空間步長Δx=402 m，時(shí)間步長Δt=5 s。計(jì)算結(jié)果詳見表4。

表4 不同弛豫時(shí)間對計(jì)算精度的影響Tab.4 The influence of different relaxation time on calculation accuracy

由表4可知，弛豫時(shí)間τ的取值過大，會導(dǎo)致誤差增大，取值過小，在特定的粒子速度（Δx/Δt）下，會導(dǎo)致計(jì)算不穩(wěn)定。因此對于洪水演算，建議弛豫時(shí)間τ的取值取在［1.5，3］范圍內(nèi)為宜。

4 結(jié) 論

本文選用擴(kuò)散波方程來描述河道洪水波運(yùn)動，采用Lattice Boltzmann 方法的D1Q5 模型來求解該偏微分方程，并給出了詳細(xì)的多尺度處理和分布函數(shù)確定步驟。采用兩個(gè)算例來說明Lattice Boltzmann 方法的應(yīng)用，同時(shí)采用Muskingum 法和解析解法來驗(yàn)證Lattice Boltzmann 方法的準(zhǔn)確性，證明了Lattice Boltzmann 方法可以很好地求解擴(kuò)散波方程和預(yù)測洪水過程。此外，分析了步長和弛豫時(shí)間的選取對計(jì)算精度的影響，結(jié)果表明，空間步長的取值比時(shí)間步長的取值對精度的影響大，建議弛豫時(shí)間τ的取值取在［1.5，3］范圍內(nèi)為宜。