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

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

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

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

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

1 2維水動(dòng)力模型

1.1 控制方程

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

式（1）～（2）中：

h

為水深；

q

和

q

分別為

x

、

y

方向的單寬流量；

g

為重力加速度；

u

、

v

分別為

x

、

y

方向的流速；

f

和

g

分別為

x

、

y

方向的通量矢量；

S

為源項(xiàng)矢量；

z

為河床底高程；

C

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

C

=

gn

/

h

，其中

n

為曼寧系數(shù)。

1.2 數(shù)值計(jì)算方法

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

i

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

式中，

Ω

為控制體

i

的面積。

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

式中，

Г

為控制體

i

的邊界，

n

為邊界

Г

所對(duì)應(yīng)外法線方向的單位向量。相應(yīng)界面通量

F

(

q

)

·n

可以表示為：

式中，

n

和

n

為邊界外法向單位向量

n

在

x

和

y

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

i

中，各通量如圖1所示，

j

、

j

、

j

分別為

i

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

i

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

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

式中，

k

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

l

為第

i

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

k

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

2 邊界條件

水動(dòng)力模型中常用的邊界條件有：入流邊界、自由出流邊界、固壁邊界、水位邊界等。本文采用虛擬單元邊界處理方法，如圖2所示，對(duì)模型中的邊界條件均采用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 固壁邊界條件

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

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

u

=

un

+

vn

和

u

=-

un

+

vn

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

2.2 自由出流邊界條件

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

2.3 水位邊界條件

通常在水位邊界上給定一個(gè)水位變化過(guò)程

z

=

z

(

t

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

h

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

2.4 單寬流量邊界條件

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

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

式中，

Q

為入口處斷面流量，

B

為入口寬度。

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

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

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

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

圖3 Toce河測(cè)點(diǎn)分布及深泓線示意圖Fig. 3 View of the gauging point and thalweg of Toce river

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

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

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

J

為原方法計(jì)算出的比降，

J

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

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

3.2 均勻流入流

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

R

=

h

，

J

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

式中，

A

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

C

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

z

時(shí)，則：

總斷面流量可表示為：

3.3 入口流量分配

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

m

個(gè)網(wǎng)格上：

式中，

l

為第

k

條入流邊的長(zhǎng)度，可得：

則第

i

條邊的流量為：

則第

i

條邊的單寬流量為：

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

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

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

圖7 Toce河入流流量過(guò)程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

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

B

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

B

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

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

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

對(duì)于不同的入流流量來(lái)說(shuō)，其過(guò)水?dāng)嗝娴乃鎸挾纫彩遣煌?，這就需要根據(jù)不同流量合理地確定水面寬度，所以本文采用了一種新的水面寬度自適應(yīng)確定方法。方法M3：在每個(gè)網(wǎng)格上以均勻流的方式入流，用二分法逼近求解對(duì)應(yīng)流量下的限定水位，用該水位計(jì)算出的水深去分配入流邊界上每個(gè)網(wǎng)格上的流量（高出限定水位的網(wǎng)格上不分流量），以確保流量守恒。再結(jié)合單寬流量邊界條件迭代求解邊界右側(cè)網(wǎng)格上的水力要素，進(jìn)而計(jì)算界面通量。方法M4：從方法M3已求得的入流邊界右側(cè)網(wǎng)格上的水位和單寬流量，不用迭代，就可直接通過(guò)式（11）計(jì)算流速，進(jìn)而計(jì)算通量。

在入流峰值時(shí)刻

t

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

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

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

5 模型應(yīng)用

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

n

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

t

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

圖11 Toce河洪水演進(jìn)過(guò)程水深分布Fig. 11 Water depth distribution of Toce river

圖12 Toce河洪水演進(jìn)過(guò)程流速分布Fig. 12 Velocity distribution of Toce river

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

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

圖14 Toce河道中各測(cè)點(diǎn)水位變化過(guò)程Fig. 14 Water level evolution process of gauging points in Toce river

6 結(jié) 論

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

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

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

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

通過(guò)本文對(duì)入流模塊的改進(jìn)與實(shí)現(xiàn)，該水動(dòng)力模型可以成功的應(yīng)用于流域洪水演進(jìn)過(guò)程的高效準(zhǔn)確模擬，為洪水災(zāi)害管理決策提供一個(gè)更直接有效的工具，更好地指導(dǎo)防洪減災(zāi)及預(yù)警預(yù)報(bào)工作。