2維水動力模型中入流邊界優(yōu)化處理方法

侯精明，汪 煜，張兆安，李丙堯，王俊琿，張大偉

(1.西安理工大學 西北旱區(qū)生態(tài)水利國家重點實驗室，陜西 西安 710048；2.河海大學 水文水資源學院，江蘇 南京 210024；3.中國水利水電科學研究院，北京 100038)

在全球氣候復雜多變和極端惡劣天氣頻發(fā)的大背景下，洪水災害給人類造成了巨大的經(jīng)濟損失，威脅著人民群眾的生命財產(chǎn)安全。為做好防洪減災和預警預報工作，2維水動力學模型作為洪水演進過程模擬的一種重要技術(shù)手段，可及時準確地為應急決策部門提供水災害信息，便于高效開展搶險救援工作。

在2維淺水流動模擬中，邊界條件處理是極其重要的部分，若邊界條件處理的不恰當，就會影響計算結(jié)果，有時甚至會阻礙計算的正常進行。所以對邊界條件進行合理有效的處理是數(shù)值模擬過程中必須面對的問題。除了因極端暴雨形成的流域及城市洪水和堤壩瞬間潰決形成的江河洪水外，對于大部分洪水演進模擬而言，若在上游入口處沒有輸入確定的水位流量關(guān)系，通常就只是給定一個入流流量過程，通過水動力模型計算來獲得水流沿河道向下游傳播時的運動狀態(tài)，但是斷面流量在入流邊界處如何合理地分配水力要素就成了研究重點。目前，入流邊界處理方法已有很多，最簡單的方法就是用斷面流量除以入口寬度去計算單寬流量，但這樣處理在入流邊界附近會出現(xiàn)高處過水、水位間斷等非物理現(xiàn)象，而且對入流邊界的寬度有一定要求。為此，于普兵提出按過水斷面面積分配各網(wǎng)格上的流量；宋利祥根據(jù)曼寧和謝才公式，以邊界上的水深和網(wǎng)格邊長計算權(quán)重來分配流量。這些方法都與入流邊界上的水深有很大關(guān)系，有學者認為邊界上的水深就等于相鄰網(wǎng)格上的平均水深值，但這樣處理對于初始為干地形的河道洪水模擬而言無法實現(xiàn)流量分配；還有的學者設(shè)置初始條件時在入流邊界網(wǎng)格上賦予一定的水深，以彌補這一缺陷。但這樣給定的水深是否合理還尚未可知，而且對于不同入口寬度也會帶來一定的問題。

針對上述入流邊界處出現(xiàn)的問題，本文采用基于非結(jié)構(gòu)網(wǎng)格有限體積法建立的2維水動力模型，對入流模塊做了部分改進，依據(jù)給定斷面流量計算出的水位來分配入流邊界各網(wǎng)格上的水力要素，使得入口附近水流流態(tài)更為合理，模型能更好地適用于洪水演進過程的模擬計算及應用。

1 2維水動力模型

1.1 控制方程

模型以2維淺水方程為控制方程，其矢量形式可表示為：

式（1）～（2）中：

h

為水深；

q

和

q

分別為

x

、

y

方向的單寬流量；

g

為重力加速度；

u

、

v

分別為

x

、

y

方向的流速；

f

和

g

分別為

x

、

y

方向的通量矢量；

S

為源項矢量；

z

為河床底高程；

C

為床面摩擦系數(shù)，

C

=

gn

/

h

，其中

n

為曼寧系數(shù)。

1.2 數(shù)值計算方法

該水動力模型使用三角形非結(jié)構(gòu)網(wǎng)格剖分計算域，采用基于Godunov格式的有限體積法對2維淺水方程進行數(shù)值求解。在網(wǎng)格單元

i

內(nèi)，其積分形式可寫為：

式中，

Ω

為控制體

i

的面積。

應用高斯定理，將式（3）中通量項的面積分轉(zhuǎn)化為邊界上的線積分：

式中，

Г

為控制體

i

的邊界，

n

為邊界

Г

所對應外法線方向的單位向量。相應界面通量

F

(

q

)

·n

可以表示為：

式中，

n

和

n

為邊界外法向單位向量

n

在

x

和

y

方向的分量。在三角形非結(jié)構(gòu)網(wǎng)格

i

中，各通量如圖1所示，

j

、

j

、

j

分別為

i

的相鄰網(wǎng)格單元，則單元

i

內(nèi)通量項線積分可寫為：

圖1 三角形網(wǎng)格上各邊通量Fig. 1 Flux of each side on triangular grids

式中，

k

為網(wǎng)格邊的編號，

l

為第

i

個網(wǎng)格單元上第

k

條邊的長度。通量項采用HLLC近似黎曼解算器求解，底坡源項和摩阻源項的具體計算方法見文獻[15–16]，本文不再贅述。此外，引入GPU加速技術(shù)以大幅提高模型計算效率。

2 邊界條件

水動力模型中常用的邊界條件有：入流邊界、自由出流邊界、固壁邊界、水位邊界等。本文采用虛擬單元邊界處理方法，如圖2所示，對模型中的邊界條件均采用HLLC近似黎曼求解器求解界面通量，邊界左側(cè)即內(nèi)部網(wǎng)格上信息為已知，需構(gòu)造邊界右側(cè)即虛擬網(wǎng)格上的信息包括地形高程、水位、流速等。模型中設(shè)置所有邊界左右兩側(cè)網(wǎng)格上地形高程相等，不同邊界右側(cè)網(wǎng)格上的水力要素則根據(jù)邊界性質(zhì)分別構(gòu)造。

圖2 邊界處理方法示意圖Fig. 2 Diagram of boundary processing method

2.1 固壁邊界條件

在固壁邊界上，采用無滑移邊界條件，法向和切向流速都為0，認為左右兩側(cè)水深相等，則右側(cè)網(wǎng)格上信息構(gòu)造為：

式中，下標R表示的是右側(cè)虛擬網(wǎng)格，下標L表示的是左側(cè)內(nèi)部網(wǎng)格，

u

=

un

+

vn

和

u

=-

un

+

vn

分別表示邊界上的法向流速和切向流速。

2.2 自由出流邊界條件

該邊界上的擾動不會對計算域內(nèi)的水流流態(tài)產(chǎn)生影響，故而邊界右側(cè)網(wǎng)格的水力要素值均采用左側(cè)網(wǎng)格的值：

2.3 水位邊界條件

通常在水位邊界上給定一個水位變化過程

z

=

z

(

t

)，可算出邊界右側(cè)每個網(wǎng)格上的水深

h

；而邊界左側(cè)網(wǎng)格上水深和流速為已知，根據(jù)1維特征線理論：

2.4 單寬流量邊界條件

3 入流邊界–斷面流量邊界條件

對于一般的河道洪水模擬來說，通常都是在上游入口處給一個斷面流量過程，但是如何合理的分配流量就成了一個問題。模型采用三角形非結(jié)構(gòu)網(wǎng)格剖分計算域，其優(yōu)勢在于能夠靈活處理復雜地形邊界，所以模型中設(shè)置為垂直于入流邊界入流。最簡單的處理辦法就是用流量除以入口寬度來計算單寬流量：

式中，

Q

為入口處斷面流量，

B

為入口寬度。

由于在入流邊界附近，該處理方式對水流流態(tài)有一定影響，故而本文對入流模塊做了部分改進，采用一種基于地形網(wǎng)格數(shù)據(jù)的縱斷面提取方法，確定河道縱比降，在入流邊界網(wǎng)格上以均勻流的方式入流，確定邊界上的水位，進而分配邊界各網(wǎng)格上的水力要素，最后耦合至模型通量部分進行計算。

3.1 河道縱斷面及縱比降確定

對于河流中的某一段來說，河底高程沿縱斷面變化趨勢大概一致；再考慮到實際地形獲取以及網(wǎng)格剖分時產(chǎn)生的誤差，導致局部河道縱比降存在較大的不確定性，故選用河道平均縱比降計算。

沿深泓線的剖面稱為河道縱斷面，能反映河床的沿程變化。在實際工作中，通常以河槽底部轉(zhuǎn)折點高程為縱坐標，以河長為橫坐標，繪制出河道縱斷面圖。本文引入一種基于地形網(wǎng)格數(shù)據(jù)的河道縱斷面提取方法，從入流邊界開始，首先，確定入流邊界處高程最低的網(wǎng)格，以此為起點和搜索中心，自動搜索所在網(wǎng)格附近的相鄰網(wǎng)格，找出起點網(wǎng)格下游處高程最低的網(wǎng)格編號；然后，以該網(wǎng)格為搜索中心，繼續(xù)重復上述搜索步驟，并依次向下游搜索，直到出流邊界為止，搜索完成，所有搜索到的最低高程網(wǎng)格中心點依次連線近似為深泓線的位置。以Toce河地形為例，如圖3所示，在這一段河道中，深泓線一共經(jīng)過534個網(wǎng)格，全長59.85 m；Toce河物理模型試驗的真實地形如圖4所示。通過對比圖3的深泓線與圖4試驗地形上深泓線可以看出，該深泓線的位置較為準確。

圖3 Toce河測點分布及深泓線示意圖Fig. 3 View of the gauging point and thalweg of Toce river

圖4 Toce河物理模型試驗地形[24]Fig. 4 Physical model topography of Toce river[24]

按照河道平均比降計算方法，確定河道縱比降，考慮到從上游入流，上游附近區(qū)域的局部地形對入口處流態(tài)有一定影響，故選取上游斷面代替原方法的下游斷面河底高程為基點做一斜線，使得斜線以下的面積與原河底線以下面積相等，該斜線坡度即為河道的平均縱比降，計算式為：

圖5為河底高程沿程變化及平均縱比降示意圖，其中，

J

為原方法計算出的比降，

J

為本文調(diào)整后計算出的比降。從圖5可以看出：對于Toce河中的某一段來說，河底高程變化趨勢確實基本一致，但在下游河道出口附近變化較為劇大。與整段河道變化趨勢相比，以上游斷面河底高程為基點的平均縱比降計算方法更為合理，由式（13）計算出Toce河該段平均比降約為1.32%。

圖5 Toce河道縱斷面及平均縱比降示意圖Fig. 5 Vertical profile and average gradient of Toce river

3.2 均勻流入流

在入流邊界各網(wǎng)格上以均勻流的方式入流，通過給定斷面流量計算出所對應的水位。認為每個網(wǎng)格上水力半徑

R

=

h

，

J

統(tǒng)一取平均比降，則每個網(wǎng)格上的流量為：

式中，

A

為入流邊各網(wǎng)格上過水面積，

C

為謝才系數(shù)。假設(shè)入流邊上的水位為

z

時，則：

總斷面流量可表示為：

3.3 入口流量分配

根據(jù)謝才公式和曼寧公式，流速近似與水深的2/3次方成正比，則在入流邊界第1至

m

個網(wǎng)格上：

式中，

l

為第

k

條入流邊的長度，可得：

則第

i

條邊的流量為：

則第

i

條邊的單寬流量為：

然后，結(jié)合單寬流量邊界條件迭代求解邊界右側(cè)網(wǎng)格上的水深和流速，進而代入HLLC近似黎曼求解器計算界面通量。

4 入流邊界不同處理方法討論

以Toce河物理模型試驗模擬為例，對入流邊界采用了4種不同的處理方法，見表1。比較各方法所對應的入口流態(tài)，確定一種合理的入流邊界處理方法。入口處局部網(wǎng)格劃分如圖6所示，流量過程如圖7所示。

圖7 Toce河入流流量過程Fig. 7 Flow discharges of Toce river

表1 入流邊界不同處理方法
Tab. 1 Different processing method of inflow boundary

方法入口寬度/m流量分配方法右側(cè)網(wǎng)格信息構(gòu)造M13.407式（12）式（10）+（11）M21.703式（12）式（10）+（11）M33.407式（16）+（21）式（10）+（11）M43.407式（16）+（21）式（15）+（11）

圖6 入口附近局部網(wǎng)格劃分示意圖Fig. 6 View of local mesh generation near the inlet

按照流量除以入口寬度計算單寬流量，方法M1：入口寬度

B

=3.407 m，在入口處水位不平順，出現(xiàn)了高處過水的現(xiàn)象，經(jīng)分析是由于入口寬度太寬的原因?qū)е碌模环椒∕2：縮短入流邊界寬度

B

=1.703 m，發(fā)現(xiàn)只在入流邊界處入流，水位在兩側(cè)非入流邊界處存在間斷，結(jié)果如圖8和9中入流邊界內(nèi)部網(wǎng)格上的水位和流速分布所示。而對于河流橫斷面來說，水位是連續(xù)的，這種高處過水和水位存在間斷的現(xiàn)象極不合理。

圖8 入流邊界上水位分布Fig. 8 Water level distribution on inflow boundary

圖9 入流邊界上流速分布Fig. 9 Velocity distribution on inflow boundary

對于不同的入流流量來說，其過水斷面的水面寬度也是不同的，這就需要根據(jù)不同流量合理地確定水面寬度，所以本文采用了一種新的水面寬度自適應確定方法。方法M3：在每個網(wǎng)格上以均勻流的方式入流，用二分法逼近求解對應流量下的限定水位，用該水位計算出的水深去分配入流邊界上每個網(wǎng)格上的流量（高出限定水位的網(wǎng)格上不分流量），以確保流量守恒。再結(jié)合單寬流量邊界條件迭代求解邊界右側(cè)網(wǎng)格上的水力要素，進而計算界面通量。方法M4：從方法M3已求得的入流邊界右側(cè)網(wǎng)格上的水位和單寬流量，不用迭代，就可直接通過式（11）計算流速，進而計算通量。

在入流峰值時刻

t

= 20 s時，采用不同處理方法所對應入流邊界內(nèi)部網(wǎng)格上的水位和流速分布如圖8、9所示，入流邊界附近局部區(qū)域的水深分布如圖10所示。

圖10 入流邊界附近局部區(qū)域水深分布Fig. 10 Water depth distribution in local area near the inlet

從圖8～10可以看出，當在模型中設(shè)置一個固定寬度的入流邊界時，對于不同的入流流量來說，若簡單的按照所有入流邊界網(wǎng)格上用斷面輸入流量除以入口寬度去計算單寬流量的方法處理，在入流邊界內(nèi)部網(wǎng)格及附近局部區(qū)域會出現(xiàn)像高處過水、水位間斷等不合理的非物理現(xiàn)象。經(jīng)過本文改進過的流量分配方法模擬結(jié)果則較為合理，而且方法M4與方法M3處理效果相似，對整體結(jié)果影響不大，避免了迭代求解過程，節(jié)省了計算時間。

5 模型應用

該物理模型真實體現(xiàn)了Toce河上游約5 km的復雜河道地形，由ENEL在意大利米蘭市按照1∶100的比尺建立，作為歐盟支持的CADAM項目的一個標準模型，常被用來檢驗水動力模型的精度和穩(wěn)定性。該物理模型尺寸約為50 m×11 m，地形如圖4所示。采用三角形網(wǎng)格剖分計算域，共生成了具有42 014個節(jié)點、82 895個單元的非結(jié)構(gòu)網(wǎng)格，平均單元面積0.004 27 m。曼寧系數(shù)

n

=0.016 2 s/m，模型初始為干地形，水深為0，左側(cè)為入流邊界，右側(cè)為自由出流邊界，其余均為固壁邊界。在河道內(nèi)布置了一系列的測點，本文選擇其中部分測點來檢驗模型精度，其測點位置如圖3所示。模擬中庫朗數(shù)CFL=0.8。圖11、12和13分別為入流時間

t

=30、45、60 s時，河道中洪水演進過程的水深、流速及流態(tài)分布。

圖11 Toce河洪水演進過程水深分布Fig. 11 Water depth distribution of Toce river

圖12 Toce河洪水演進過程流速分布Fig. 12 Velocity distribution of Toce river

圖13 Toce河洪水演進過程流態(tài)分布Fig. 13 Flow state distribution of Toce river

圖11～13是采用優(yōu)化方法處理后水深、流速和流態(tài)的計算結(jié)果，從圖11～13中可以看出：在Toce河物理模型潰壩試驗洪水演進過程模擬中，急緩流交替，流態(tài)復雜，極其考驗著模型的準確性和穩(wěn)定性。圖14為河道中不同位置測點的水位隨時間變化過程。通過與觀測數(shù)據(jù)的比較可以看出，水位的計算值與試驗測量值總體上吻合較好，變化趨勢基本一致；從測點P1到P26的洪水傳播時間為42 s，與試驗觀測值40 s極為接近，采用納什效率系數(shù)（NSE）來評價模型質(zhì)量，經(jīng)計算，河道中5個測點的納什效率系數(shù)范圍為0.82～0.95，說明該模型計算精度較高，驗證了該入流邊界處理方法的可行性和模型的準確性；而且在單機上僅用時86 s就完成了該算例的模擬計算，表明該模型可以高效準確的模擬洪水演進過程，對于洪水災害管理決策意義重大。

圖14 Toce河道中各測點水位變化過程Fig. 14 Water level evolution process of gauging points in Toce river

6 結(jié) 論

為了高效準確的模擬洪水演進過程，本文使用一種較為合理的入流邊界流量分配方法，采用基于非結(jié)構(gòu)網(wǎng)格有限體積法建立2維水動力模型，改進了模型中入流模塊，并引入GPU加速技術(shù)對模型算法進行處理，以提高計算效率，將其應用于Toce河物理模型試驗的模擬，結(jié)果表明：

本文采用一種基于地形網(wǎng)格數(shù)據(jù)的河道縱斷面提取方法以確定河道比降，在入流邊網(wǎng)格上通過均勻流的方式入流，確定水位，依據(jù)入流邊界各網(wǎng)格上的水深和邊長去計算權(quán)重以分配流量。從Toce河試驗模擬結(jié)果可以看出，在上游入口處，水流流態(tài)更為合理。

從河道中各測點水位的模擬值與觀測值的比較可以看出，水位變化趨勢基本一致，洪水到達及傳播時間較為接近，河道中各測點的納什效率系數(shù)在0.82～0.95之間，驗證了該方法的可行性和模型的準確性。

本文利用GPU加速技術(shù)對模型算法進行處理，在單機上用時86 s完成了具有82 895個計算單元的復雜地形上180 s入流洪水演進過程的模擬試驗，計算效率較高。

通過本文對入流模塊的改進與實現(xiàn)，該水動力模型可以成功的應用于流域洪水演進過程的高效準確模擬，為洪水災害管理決策提供一個更直接有效的工具，更好地指導防洪減災及預警預報工作。