摘要：探究蓄滯洪區(qū)洪水演進(jìn)對(duì)蓄滯洪區(qū)的防洪工作有著重要意義。為在洪水來(lái)臨時(shí)，短時(shí)間內(nèi)對(duì)洪水運(yùn)動(dòng)進(jìn)行預(yù)測(cè)。運(yùn)用 Godunov 型有限體積法建立了1個(gè)二維淺水水動(dòng)力數(shù)值模型，該模型的的優(yōu)勢(shì)在于采用非結(jié)構(gòu)性三角形網(wǎng)格，能夠較好地?cái)M合不規(guī)則邊界，并采用 GPU 異構(gòu)并行技術(shù)，有效解決了非結(jié)構(gòu)化網(wǎng)計(jì)算耗時(shí)長(zhǎng)的問(wèn)題，運(yùn)用限制器減少二階精度數(shù)值的擴(kuò)散和振蕩，引入負(fù)水深和干濕邊界處理技術(shù)，使模型模擬的精確度得到了有效提高。模型應(yīng)用多個(gè)潰壩水流實(shí)驗(yàn)算例對(duì)模型進(jìn)行率定，在驗(yàn)證其具有良好的干濕邊界計(jì)算能力和靜水和諧性的基礎(chǔ)上，采用嫩江流域1998年洪水?dāng)?shù)據(jù)對(duì)嫩江流域的胖頭泡蓄滯洪區(qū)的洪水演進(jìn)進(jìn)行模擬計(jì)算。結(jié)果顯示：模型對(duì)胖頭泡蓄滯洪區(qū)模擬的蓄水量為34.9億 m3，與實(shí)測(cè)值的相對(duì)誤差為1.7%，洪水淹沒面積為1267.27 km2，與實(shí)測(cè)值的相對(duì)誤差為9.24%。對(duì)研究區(qū)域不同的土地類型進(jìn)行糙率分區(qū)的情況下，模型模擬的蓄水量結(jié)果更趨近于實(shí)測(cè)結(jié)果。并行化后的程序能使模型計(jì)算速度提升33%，因此模型能快速、精準(zhǔn)計(jì)算復(fù)雜地形潰決水流的洪水演進(jìn)過(guò)程。研究結(jié)果可以為蓄滯洪區(qū)的防洪規(guī)劃、災(zāi)情預(yù)警提供相關(guān)信息。

關(guān)鍵詞：淺水方程；有限體積法；蓄滯洪區(qū)；洪水演進(jìn)

中圖分類號(hào)：TV133 文獻(xiàn)標(biāo)識(shí)碼：A文章編號(hào)：1001-9235（2024）06-0082-10

Numerical Simulation of Flood Evolution in Pangtoupao Flood Storage Area Based on Two-Dimensional Shallow Water Equation

SUN Zhenyu1， ZHOU Yunhao1， ZHANG Qinxu1， QIN Mengen1， ZHANG Mingliang1，2

（1. School of Ocean Science and Environment， Dalian Ocean University， Dalian 116023， China;2. Liaoning Coastal EcologicalEnvironment and Disaster Protection Engineering Technology Innovation Center， Dalian 116023， China）

Abstract： Exploring the flood evolution in flood storage areas is of great significance for flood control work in flood storage areas. To predict the dynamics of incoming floods within a short period of time， The Godunov-type finite volume method was used to establish a two-dimensional hydrodynamic numerical model of shallow water. The model had the advantages of adopting unstructured triangular mesh， which could better fit irregular boundaries， and it adopted GPU heterogeneous parallel technology， which effectively solved the problem of slow computation speed of an unstructured network. The model used a limiter to reduce the diffusion and oscillation of second-order precision values. The accuracy of the model simulation was improved effectively by using negative water depth and dry and wet boundary processing techniques. The model was calibrated using multiple dam-break flow experiments and validated with an excellent capability in calculating dry and wet boundaries and static water balance. Flood data from 1998 in the Nenjiang River Basinwas used to simulate and calculate flood evolution in the Pangtoupo flood storage area of the Nenjiang River Basin. The results show that the model simulates a storage capacity of 34.9 billion m3 for the Pangtoupo flood storage area， with a relative error of 1.7% compared to the measured value， and a flooded area of 1267.27 km2， with a relative error of 9.24% compared to the measured value. When different land types in the study area are classified based on roughness coefficients， the simulated storage capacity results of the model approach to the measured values. The parallelized program improves the model′s calculation speed by 33%. The model thus can rapidly and accurately calculate the flood evolution process of dam-break flows in complex terrain. The results provide relevant information for flood control planning and disaster warning in flood storage areas.

Keywords： shallow water equation; finite volume method; flood storage area; flood evolution

洪水一直以來(lái)都是全球最主要的自然災(zāi)害之一，目前很多專家和學(xué)者依據(jù)水動(dòng)力學(xué)原理建立了一維和二維水動(dòng)力模型，用于在洪水來(lái)臨時(shí)提供有效的預(yù)報(bào)和預(yù)警。在研究一維水動(dòng)力模型方面，大部分學(xué)者［1-4］是對(duì)河道、河網(wǎng)和水庫(kù)等區(qū)域進(jìn)行洪水演進(jìn)模擬和驗(yàn)證，但在研究大尺度水體區(qū)域過(guò)程中，如湖泊、蓄滯洪區(qū)、河口等，其水深、流動(dòng)速度等要素在水平方向的改變會(huì)遠(yuǎn)大于垂直方向上的改變，所以僅用一維模型很難精準(zhǔn)的模擬洪水演進(jìn)過(guò)程；若采用三維模型模擬，雖然可以更真實(shí)的反映水流運(yùn)動(dòng)特性，但對(duì)于垂向沿程變化不明顯的水體，例如陸面洪水等，其垂直方向上的水深、流速等各要素值接近平均分布［5］，所以使用三維模型計(jì)算會(huì)較為復(fù)雜且計(jì)算速度效率將大幅降低，因此目前多采用二維模型對(duì)研究區(qū)域進(jìn)行模擬計(jì)算。現(xiàn)如今對(duì)二維淺水方程組的二維模型研究和應(yīng)用越來(lái)越廣泛，因其在洪水模擬中的計(jì)算效率以及計(jì)算精度遠(yuǎn)高于簡(jiǎn)化模型，可以更精準(zhǔn)地描述洪水的運(yùn)動(dòng)情況［6-7］，尤其是基于 Godunov 型有限體積法建立的淺水方程模型，其具有良好的捕捉動(dòng)邊界的能力，在模擬洪水運(yùn)動(dòng)中的瞬時(shí)陸地流量等方面具有很大的潛力［8-9］。目前眾多國(guó)內(nèi)外學(xué)者［10-16］都基于有限體積法對(duì)二維淺水方程模型進(jìn)行改進(jìn)，用于洪水演進(jìn)數(shù)值模擬當(dāng)中。平面二維淺水方程的應(yīng)用可以較為恰當(dāng)?shù)孛枋鰷\水水流運(yùn)動(dòng)，在滿足模擬精度的同時(shí)，計(jì)算效率也能得到相應(yīng)的提高。由于洪水水流其流速快以及強(qiáng)間斷性的特點(diǎn)，模型更適合于采用顯格式計(jì)算，而顯格式計(jì)算受到時(shí)間步長(zhǎng)及網(wǎng)格精度的限制；在研究地形不規(guī)則的研究區(qū)域中，相對(duì)于結(jié)構(gòu)化網(wǎng)格而言，非結(jié)構(gòu)化網(wǎng)格具有更好的擬合邊界的能力，網(wǎng)格質(zhì)量相對(duì)更好，但對(duì)計(jì)算機(jī)的性能要求以及計(jì)算成本較高，需要耗費(fèi)更長(zhǎng)的時(shí)間。

基于非結(jié)構(gòu)性三角形網(wǎng)格，并運(yùn)用 GPU異構(gòu)并行技術(shù)解決運(yùn)用非結(jié)構(gòu)化網(wǎng)格模型計(jì)算緩慢問(wèn)題。使用 Roe 格式計(jì)算界面通量，運(yùn)用 Godunov 有限體積法建立出1個(gè)二維淺水動(dòng)力模型，并且引入了負(fù)水深和干濕邊界處理技術(shù)，應(yīng)用模型模擬多個(gè)潰壩水流實(shí)驗(yàn)室算例，對(duì)模型進(jìn)行率定，驗(yàn)證模型穩(wěn)定性以及準(zhǔn)確性。利用衛(wèi)星遙感技術(shù)識(shí)別胖頭泡蓄滯洪區(qū)的地貌景觀特征，對(duì)不同的區(qū)域設(shè)置其相應(yīng)的糙率值，從而保證模型的精準(zhǔn)度，模擬在相同或不同的糙率條件下，洪水在蓄滯洪區(qū)的運(yùn)動(dòng)狀況及洪水的衰減作用。

1水動(dòng)力模型

1.1二維淺水控制方程

遵循淺水Boussinesq假設(shè)與靜壓假定，對(duì) Navier-Stokes方程沿垂直方向流速進(jìn)行深度平均近似，在忽略風(fēng)應(yīng)力、科氏力以及和二階擴(kuò)散項(xiàng)的情況下，二維淺水控制方程最終表達(dá)式見式（1）、（2）。

+ + = S（1）

式中：U 為守恒性變量的向量；F 為 x 的對(duì)流通量；G 為 y 的對(duì)流通量；S 為源項(xiàng)。

?0?

U =（ v（u）），F(xiàn) =（ v（2）），G =（ 2（v）），S = -- g（g）h（h） --τ（τ）by（bx）（2）

式中：h 為水深；u 為 x 方向的流速，v 為 y 方向的流速；τbx 為 x 方向的底部摩阻項(xiàng)，τby 為 y 方向的底部摩阻項(xiàng)；g 為重力加速度；η為水位；n 為曼寧系數(shù)。其表達(dá)式見式（3）。

τbx = gn2 uu2 + v2， τby = gn2 vu2 + v2（3）

1.2有限體積離散

用非結(jié)構(gòu)性三角網(wǎng)格對(duì)研究區(qū)域進(jìn)行劃分，用 Godunov 型有限體積法對(duì)控制方程進(jìn)行數(shù)值離散和求解，使淺水方程在非結(jié)構(gòu)三角網(wǎng)格上進(jìn)行積分可得式（4）。

dv +?·Eadv dv =?·Edif dv + Sdv（4）

式中：E adv 為對(duì)流通量；Eaif為擴(kuò)散通量。運(yùn)用 Green 公式將方程進(jìn)行積分整理后可得式（5）。

ΔUi = （E ij（*）?nij）?lij +Sdv（5）

式中：Ai 為第i個(gè)單元的面積；m 為單元邊的個(gè)數(shù)； E*ij·nij為第i個(gè)單元第j 條邊的通量；lij為單元每條邊的長(zhǎng)度；Vi 為控制單元；Ui 為控制單元平均值，其表達(dá)式見式（6）。

Ui =Udv（6）

1.3動(dòng)量通量計(jì)算

由于非結(jié)構(gòu)網(wǎng)格在每個(gè)控制單元邊界存在有間斷現(xiàn)象，從而構(gòu)成 Riemann 問(wèn)題。本文采用 Roe 格式求解 Riemann 界面法向數(shù)值通量問(wèn)題［17］，數(shù)值通量表達(dá)式見式（7）。

E* ? n = （F，G ）R? n +（F，G ）L? n -| J（?）|（UL -

UR ）]（7）

式中：J（?）為 Roe 平均的 Jacobian 矩陣；UL、UR 為單元邊界兩側(cè)的 Riemann 守恒性變量。

1.4動(dòng)邊界處理及負(fù)水深處理

干濕界面是指在計(jì)算區(qū)域內(nèi)有水和無(wú)水交替變化或部分淹沒的區(qū)域，即動(dòng)邊界。為了真實(shí)的模擬在天然河流，河口以及海灣處發(fā)生的洪水演進(jìn)、潰壩水流演進(jìn)等現(xiàn)象，從而引入干濕界面處理技術(shù)。針對(duì)復(fù)雜地形引入了干濕界面處理技術(shù)［18］，以h0為最小限制水深，若水深低于 h0則為干邊界；高于h0則為濕邊界。圖1為4種干濕邊界情況。

根據(jù)干濕邊界的4種情況判斷，可以把計(jì)算網(wǎng)格分為3種網(wǎng)格：網(wǎng)格的3個(gè)頂點(diǎn)全為濕點(diǎn)，且3個(gè)邊界全為濕邊界或半濕邊界組成，則此網(wǎng)格為濕網(wǎng)格；網(wǎng)格的3個(gè)邊界全為干邊界或者半濕邊界組成，則此網(wǎng)格為干網(wǎng)格；除濕網(wǎng)格與干網(wǎng)格外其余均為半濕網(wǎng)格。

基于淺水方程的模擬過(guò)程當(dāng)中，邊界通量的計(jì)算是由時(shí)間步長(zhǎng)、單寬通量及邊界長(zhǎng)度決定。如果時(shí)間步長(zhǎng)設(shè)置過(guò)小，會(huì)使模型計(jì)算速度大幅降低。在較小水深的情況下，時(shí)間步長(zhǎng)如果設(shè)置過(guò)大，會(huì)導(dǎo)致計(jì)算不穩(wěn)定，并且會(huì)造成網(wǎng)格邊界的通量高于上游網(wǎng)格的水量，從而導(dǎo)致上游網(wǎng)格形成負(fù)水深，最后導(dǎo)致質(zhì)量不守恒。針對(duì)負(fù)水深問(wèn)題，在計(jì)算過(guò)程中記錄下負(fù)水深單元網(wǎng)格 Ai，并將這些單元網(wǎng)格的水深調(diào)整為 h0（最小限制水深），流速調(diào)置為0，之后將單元網(wǎng)格Ai 周圍的共邊網(wǎng)格中高于最小限制水深的水量補(bǔ)入單元網(wǎng)格 Ai 當(dāng)中，直至負(fù)水深成為最小限制水深。

1.5并行化處理

由于二維水動(dòng)力模型數(shù)據(jù)復(fù)雜，計(jì)算速度緩慢，因此本模型采用 GPU異構(gòu)并行技術(shù)實(shí)現(xiàn)高速運(yùn)算。在計(jì)算過(guò)程中，CPU 作為主要領(lǐng)導(dǎo)，負(fù)責(zé)產(chǎn)生和交付多線程任務(wù)，而 GPU 作為計(jì)算任務(wù)的執(zhí)行者，將任務(wù)分配至 GPU 硬件的計(jì)算單元中進(jìn)行計(jì)算，兩者進(jìn)行必要的信息傳遞，完成整個(gè)計(jì)算流程。OpenACC是一種高級(jí)編程模型，可通過(guò)指令來(lái)加速C/C++和 Fortran語(yǔ)言中的并行部分。其原理是先在主機(jī)端的加速器上分配空間，然后將數(shù)據(jù)從 CPU 的內(nèi)存復(fù)制到 GPU 的內(nèi)存，并發(fā)送相應(yīng)的代碼。這些數(shù)據(jù)和代碼在加速器上排隊(duì)等待空閑資源來(lái)進(jìn)行加速計(jì)算，其并行過(guò)程主要包括3個(gè)部分：首先是對(duì)變量參數(shù)、網(wǎng)格信息和流量等數(shù)據(jù)進(jìn)行初始化；其次是通量計(jì)算和項(xiàng)源求解；最后輸出網(wǎng)格單元的水深、流速等信息，模型計(jì)算框架見圖2。在本模型中，使用 kernel、loop、private 和 routine 指令來(lái)完成并行計(jì)算，模型中每個(gè)附入代碼的區(qū)域都被注釋（！ acc）來(lái)控制在 GPU 上的并行計(jì)算。在內(nèi)部循環(huán)計(jì)算過(guò)程中，使用 kernels loop 組合指令將計(jì)算區(qū)域加載到 GPU 上，成為一個(gè)能執(zhí)行的加速函數(shù)，使用！ acc routine seq 指令用來(lái)分配加速設(shè)備上的內(nèi)存。由于本模型使用的是顯格式求解控制方程，所以在網(wǎng)格計(jì)算過(guò)程中，每個(gè)單元網(wǎng)格的計(jì)算只需前一步時(shí)間步長(zhǎng)的結(jié)果，因此使用OpenACC模式的并行算法可以實(shí)現(xiàn)對(duì)本模型的并行計(jì)算。

2數(shù)值模擬驗(yàn)證與應(yīng)用

2.1在 90°彎道下的潰壩面模擬

在潰壩模型應(yīng)用中，特別是大尺度區(qū)域進(jìn)行潰壩水流數(shù)值模擬時(shí)，可能會(huì)遇到在山區(qū)河道中出現(xiàn)急彎的情況，為了研究本模型是否能處理地形存在突變的潰壩水流，模型模擬了上游矩形水庫(kù)潰壩實(shí)驗(yàn)，模型示意見圖3。上游水庫(kù)幾何尺寸為244 cm×239 cm，下游河道為寬49.5 cm 的 L 型矩形斷面，河道上端長(zhǎng)約400 cm，下端長(zhǎng)約300 cm，下游河道高程比上游水庫(kù)高程高33 cm。河道與水庫(kù)內(nèi)無(wú)坡度，河道出口端為自由出口，試驗(yàn)前河道設(shè)置為濕河床，水庫(kù)初始水位高于河道20 cm，初始底部粗糙值設(shè)置為0.012。試驗(yàn)時(shí)將河道與水庫(kù)的連接閘門快速打開從而產(chǎn)生潰壩水流。圖4為在不同測(cè)點(diǎn)（G1—G6）水位的實(shí)測(cè)值和模擬值隨時(shí)間的變化圖像。

本算例總計(jì)算時(shí)長(zhǎng)為40 s，圖中 G1點(diǎn)為上游水庫(kù)當(dāng)中的測(cè)點(diǎn)，由圖4a可以看到其測(cè)點(diǎn)的水位會(huì)隨著閘門打開呈現(xiàn)持續(xù)降低的趨勢(shì)，水位下降的幅度隨著時(shí)間的增加而逐步減少。G2、G3和 G4測(cè)點(diǎn)都位于河道上端，都呈現(xiàn)出了一種水位突然增大的現(xiàn)象，其原因是下游彎道對(duì)水波的反射作用所造成的，其水位均是到達(dá)極大值后經(jīng)過(guò)一段時(shí)間再突然增大到最大值，之后再緩慢降低；由圖4b—4d 可看出測(cè)點(diǎn) G4的水位最大值出現(xiàn)得較早，其原因是 G4點(diǎn)相對(duì)于 G2和 G3測(cè)點(diǎn)離 L 形轉(zhuǎn)彎處的距離更近。 G5、G6均在彎道后方的河道下游，由于水流受到彎道的調(diào)整作用，其水位的變化過(guò)程相對(duì)于其他測(cè)點(diǎn)來(lái)說(shuō)更平緩了些，模型計(jì)算結(jié)果與實(shí)測(cè)結(jié)果大致相同。從整體結(jié)果來(lái)看，本模型能夠模擬類似于這種地形存在突變的潰壩水流實(shí)驗(yàn)。

2.2過(guò)駝峰的潰壩波傳播模擬

本算例是為了驗(yàn)證模型的穩(wěn)定性，以及動(dòng)態(tài)邊界處理的合理性與有效性。計(jì)算區(qū)域?yàn)?5 m×30 m 的矩形封閉水槽，其中大壩位于 x=16 m 處，大壩的厚度忽略不計(jì)。水槽底部的曼寧系數(shù)為0.018，四周采用固壁邊界。水槽有3個(gè)駝峰地形，其中2個(gè)高度為1 m 的低駝峰分別位于（30 m，6 m）和（30 m，24 m），一個(gè)高度為3 m 的高駝峰位于（47.5 m，15 m）。駝峰高程的表達(dá)式見式（8）。

試驗(yàn)時(shí)，大壩上游初始水位為1.875 m。試驗(yàn)?zāi)M計(jì)算了300 s 內(nèi)水流運(yùn)動(dòng)情況。t=2、6、12、24、30、300 s 的二維水深結(jié)果見圖5。由結(jié)果可看出，水流特征具有良好的對(duì)稱性，符合水流運(yùn)動(dòng)規(guī)律。當(dāng)潰壩波傳播時(shí)高駝峰未曾被水淹沒，但高度為1 m 的低駝峰曾被水完全淹沒，低駝峰的漲水和退水現(xiàn)象比較明顯。t=2 s 時(shí)，水流通過(guò)小駝峰，并處于漲水狀態(tài)；t=6 s 時(shí)，小駝峰已被完全淹沒，潰壩波傳播至高峰處；t=12 s 時(shí)，洪水繞高峰兩側(cè)流過(guò)，且繞流現(xiàn)象較為明顯；t=24、30 s 時(shí)，洪水已完全淹沒平底區(qū)域，水槽末端固壁邊界對(duì)水流有明顯的反射現(xiàn)象，t=300 s 時(shí)，由于潰壩波與邊界相互作用和能量的消散，研究區(qū)域內(nèi)的水流幾乎成靜止?fàn)顟B(tài)。

2.3胖頭泡蓄滯洪區(qū)洪水演進(jìn)

胖頭泡蓄滯洪區(qū)作為哈爾濱市重要的防洪體系之一，其作用是承擔(dān)松花江哈爾濱市發(fā)生巨大洪水時(shí)的防洪工作，從而減輕松花江的防洪壓力。利用1998年胖頭泡蓄滯洪區(qū)的歷史資料，使用本模型對(duì)1998年大洪水中胖頭泡蓄滯洪區(qū)的洪水演進(jìn)過(guò)程進(jìn)行了數(shù)值模擬，其中還包括有網(wǎng)格分布、地形插值、分區(qū)域糙率設(shè)置、邊界的條件以及初始條件的設(shè)置等。從淹沒面積、淹沒最大水深，蓄水量等方面，將模擬結(jié)果與1998年洪水?dāng)?shù)據(jù)進(jìn)行對(duì)比，驗(yàn)證模型的合理性。

2.3.1網(wǎng)格敏感度分析

相對(duì)于結(jié)構(gòu)化網(wǎng)格，非結(jié)構(gòu)化網(wǎng)格能更好的模擬不規(guī)則邊界，且網(wǎng)格的質(zhì)量更好，但非結(jié)構(gòu)網(wǎng)格在計(jì)算速度上會(huì)低于結(jié)構(gòu)化網(wǎng)格，使其在使用過(guò)程中對(duì)計(jì)算機(jī)的性能要求較高，且計(jì)算的成本也較大。因此本文在保證模型計(jì)算精準(zhǔn)度的前提下，為了提高模型計(jì)算速度；節(jié)約計(jì)算成本，本文對(duì)網(wǎng)格的收斂性進(jìn)行了分析。最終選出網(wǎng)格數(shù)量分別為13240、16877、20303個(gè)的3套網(wǎng)格，分別對(duì)其地理位置相同的4個(gè)測(cè)點(diǎn)的洪水水深和洪水淹沒面積進(jìn)行分析比對(duì)，4個(gè)測(cè)點(diǎn)地理位置見圖6，淹沒水深變化圖和洪水淹沒面積見圖7、8。從結(jié)果來(lái)看，隨著網(wǎng)格數(shù)量的增加，模型計(jì)算時(shí)間也隨之增加，模擬時(shí)長(zhǎng)為200 h 時(shí)，13240、16877、20303個(gè)的網(wǎng)格所用計(jì)算時(shí)間分別需為54、74、90 min，但其水深及淹沒面積的計(jì)算結(jié)果并未發(fā)生較大變化。因此在不同網(wǎng)格精度對(duì)模擬結(jié)果影響結(jié)果不大的前提下，選擇計(jì)算速度較快且計(jì)算成本較低的13240個(gè)網(wǎng)格完成模擬計(jì)算。

2.3.2變化曼寧系數(shù)分析

糙率系數(shù)是指研究區(qū)域中各種地物對(duì)水流阻力產(chǎn)生作用的1個(gè)綜合系數(shù)，因此糙率系數(shù)的取值要根據(jù)研究區(qū)域的地形地貌、各地區(qū)土地實(shí)際利用情況、植被覆蓋情況及洪水泛濫的季節(jié)等進(jìn)行初步分析，從而劃定不同土地類型的糙率取值。

由于胖頭泡蓄滯洪區(qū)的土地情況較為復(fù)雜，故本文將地理信息系統(tǒng)（GIS）和遙感影像處理等技術(shù)應(yīng)用于胖頭泡蓄滯洪區(qū)洪水淹沒的研究當(dāng)中。根據(jù)衛(wèi)星遙感影像，可將研究區(qū)域大致分為四類：水體、濕地、植被和人工建筑地區(qū)。經(jīng)過(guò)查證，胖頭泡蓄滯洪區(qū)植被主要為玉米等農(nóng)作物，而位于東北的玉米一般在每年8月份收獲，根據(jù)《1998年松花江大洪水》的記載，洪水發(fā)生的時(shí)間為八月至九月初。因此，當(dāng)洪水發(fā)生時(shí)玉米正處于成熟階段，故在本文中將植被的糙率值設(shè)置為0.065。其他區(qū)域的糙率值參考了實(shí)際情況及文獻(xiàn)［19-22］，設(shè)置如下：水體：0.03；濕地：0.04；人工建筑區(qū)：0.09。由于裸地面積相對(duì)于大尺度的研究區(qū)面積極小，對(duì)研究區(qū)域洪水演進(jìn)影響也非常小，同時(shí)考慮到小區(qū)域土地類型特征的設(shè)置無(wú)法在大尺度網(wǎng)格上得到滿足，故在此次模擬中忽略了裸地糙率值的設(shè)定。衛(wèi)星影像及土地利用情況見圖9。為了探究在胖頭泡蓄滯洪區(qū)設(shè)置全域糙率和分區(qū)糙率對(duì)模型模擬的影響作用，根據(jù)其地貌特點(diǎn)結(jié)合《松花江流域防汛資料匯編》對(duì)全域糙率設(shè)置為 n=0.055和 n=0.065。將模擬結(jié)果與糙率分區(qū)模擬結(jié)果結(jié)合《1998年松花江大洪水》中提供的分洪和退洪數(shù)據(jù)推導(dǎo)出的各個(gè)時(shí)刻蓄水量的結(jié)果進(jìn)行對(duì)比，其對(duì)比結(jié)果見圖10。結(jié)果顯示，在糙率分區(qū)的情況下，模擬洪水在340 h 的蓄水量為34.9億 m3，與《1998年松花江大洪水》記錄的35.5億 m3較為接近，相對(duì)誤差為1.7%。

2.3.3胖頭泡洪水演進(jìn)模擬分析

本模型以1998年松花江實(shí)測(cè)洪水資料為依據(jù)，在糙率分區(qū)的情況下對(duì)胖頭泡蓄滯洪區(qū)進(jìn)行數(shù)值模擬計(jì)算。由資料可知，本次洪水決口位于大安市上游、蓄滯洪區(qū)國(guó)營(yíng)農(nóng)場(chǎng)附近的嫩江左岸堤防處附近，其流量過(guò)程及決口位置見圖11。本次模擬總時(shí)長(zhǎng)為340 h，計(jì)算結(jié)果選取100、150、200和340 h 4個(gè)時(shí)刻的洪水淹沒范圍以及流場(chǎng)分布圖來(lái)演示洪水的演進(jìn)過(guò)程，胖頭泡蓄滯洪區(qū)地形高程（Z）變化及流場(chǎng)分布見圖12。從圖12中可以看出，洪水從胖頭泡決口進(jìn)入蓄滯洪區(qū)后，分洪初期洪水向著東和東北方向流入低洼地區(qū)，后期洪水主要向著南和東南方向演進(jìn)，洪水演進(jìn)過(guò)程符合胖頭泡蓄滯洪區(qū)西北高、東南低的地勢(shì)地貌特征。在模擬過(guò)程中，串行計(jì)算結(jié)果與并行計(jì)算結(jié)果基本一致，在使用并行計(jì)算時(shí)，計(jì)算效率提升了33%。不同時(shí)刻淹沒面積計(jì)算結(jié)果與其他文獻(xiàn)［22］中的 MIKE、EFDC模型的計(jì)算結(jié)果對(duì)比圖13，以《1998年松花江大洪水》中洪水調(diào)查量算的淹沒面積為依據(jù)，3種模型模擬淹沒面積的結(jié)果與實(shí)測(cè)數(shù)據(jù)的相對(duì)差值見表1。從結(jié)果數(shù)據(jù)來(lái)看，本模型與實(shí)測(cè)值1160 km2的相對(duì)誤差為9.24%，在誤差的允許范圍之內(nèi)，故本模型計(jì)算結(jié)果合理。模型計(jì)算所得的最大水深分布面積結(jié)果，與其他文獻(xiàn)［22］中 MIKE、EFDC模型的計(jì)算結(jié)果對(duì)比結(jié)果見圖14。圖13、14可知，二維淺水控制方程模擬結(jié)果與王靜［22］的模型對(duì)1998年大洪水進(jìn)行建模計(jì)算的結(jié)果比較一致?？傮w來(lái)說(shuō)，模型能夠模擬復(fù)雜地形上的洪水運(yùn)動(dòng)過(guò)程，可以為復(fù)雜地形的洪水演進(jìn)進(jìn)行預(yù)測(cè)。

3結(jié)論

基于有限體積法以 Roe格式求解 Riemann界面法向數(shù)值通量問(wèn)題，模型采用了干濕邊界處理以及負(fù)水深處理技術(shù)，進(jìn)一步提高了模型地穩(wěn)定性以及靜水和諧性；在非結(jié)構(gòu)三角網(wǎng)格控制體上用有限體積法對(duì)方程進(jìn)行離散，使模型能夠較好的擬合類似蓄滯洪區(qū)等復(fù)雜地形及不規(guī)則邊界，通過(guò)引入OpenACC模式的并行算法，成功實(shí)現(xiàn)了對(duì)該模型的加速計(jì)算，使模型計(jì)算效率提升了33%。將模型應(yīng)用于胖頭泡蓄滯洪區(qū)洪水演進(jìn)數(shù)值模擬后發(fā)現(xiàn)：洪水由潰口流入胖頭泡蓄滯洪區(qū)后，洪水先向著東和東北方向流入低洼地區(qū)，有部分洪水順勢(shì)向南和東南方向順勢(shì)而下，后期洪水主體主要向東南和南演進(jìn)，到達(dá)排水口后有少部分洪水向西擴(kuò)散，分洪212 h后洪水由排水口回歸嫩江，這與胖頭泡蓄滯洪區(qū)西北高東南低的地勢(shì)相符合，其洪水流量與淹沒面積呈正相關(guān)。在分洪340 h 的過(guò)程當(dāng)中，模型模擬洪水淹沒面積和蓄水量的結(jié)果與測(cè)量數(shù)據(jù)基本符合，且洪水演進(jìn)過(guò)程中的流場(chǎng)分布合理。此模型可以對(duì)大區(qū)域復(fù)雜地形的洪水運(yùn)動(dòng)進(jìn)行快速預(yù)測(cè)，其結(jié)果具有可靠性、合理性，本模型可以為蓄滯洪區(qū)的防洪工作、預(yù)報(bào)預(yù)警提供科學(xué)的依據(jù)。

參考文獻(xiàn)：

［1］季益柱，丁全林，王玲玲，等．三峽水庫(kù)一維水動(dòng)力數(shù)值模擬及可視化研究［J］．水利水電技術(shù)，2012，43（11）：21-24.

［2］劉恒恒，齊鵬云，萬(wàn)能勝，等．巢湖閘下河網(wǎng)洪水演進(jìn)及閘門調(diào)度數(shù)值模擬［J］．人民長(zhǎng)江，2022，53（5）：41-46.

［3］胡嘉鏜，李適宇．珠江三角洲一維鹽度與三維斜壓耦合模型［J］．水利學(xué)報(bào)，2008（11）：1174-1182.

［4］王盼，何洋，杜志水．基于水動(dòng)力數(shù)值計(jì)算的城市設(shè)計(jì)洪水模擬研究［J］．西安理工大學(xué)學(xué)報(bào)，2020，36（3）：362-366.

［5］張明亮．近海及河流環(huán)境水動(dòng)力數(shù)值模擬方法與應(yīng)用［M］.北京：科學(xué)出版社，2015.

［6］ FAN Y Y，AO T Q，YU H J，et al． A coupled 1D-2D hydrodynamic model for urban flood inundation［J］． Advances in Meteorology，2017． DOI：10.1155/2017/2819308.

［7］ GUAN M，WRIGHT N G，SLEIGH P A． A robust 2D shallow water model for solving flow over complex topography using homogenousfluxmethod ［J］．InternationalJournal for Numerical Methods in Fluids，2013，73（3）：225-249.

［8］ OALBARRIEAT M，MEDINA R，GONZALEZ M，et al． A finite volume-finite difference hybrid model for tsunami propagation and run up［J］．Computers and GeoSciences，2010，37（8）：1003-1014.

［9］張明亮．濱海鹽沼濕地退化機(jī)制及生態(tài)修復(fù)技術(shù)研究進(jìn)展［J］．大連：大連海洋大學(xué)學(xué)報(bào)，2022，37（4）：539-549.

［10］LI D M，ZHEN Z，ZHANG H Q，et al． Hydrodynamic model of Daya Bay based on finite element method［J］. Environmental Earth Sciences，2020，79（1）. DOI：10.1007/S12665-020-09019-X.

［11］GALLEGOS H A，SCHUBERT J E，SANDERS B F． Two- dimensional， high-resolutionmodelingofurbandam-break flooding：A case study of Baldwin Hills，California［J］． Advances in Water Resources，2009，32（8）：1323-1335.

［12］陳海鑫．二維淺水流動(dòng)數(shù)值模擬研究及其應(yīng)用［D］．大連：大連理工大學(xué)，2013.

［13］LIANG Q H． Flood simulation using a well-balanced shallow flow model［J］． Journal of Hydraulic Engineering，2010，136（9）：669-675.

［14］宋利祥，周建中，王光謙，等．潰壩水流數(shù)值計(jì)算的非結(jié)構(gòu)有限體積模型［J］．水科學(xué)進(jìn)展，2011，22（3）：373-381.

［15］張大偉，程曉陶，黃金池，等．基于 Godunov格式的潰壩水流數(shù)學(xué)模型［J］．水科學(xué)進(jìn)展，2010，21（2）：167-172.

［16］侯精明，張兆安，馬利平，等．基于 GPU 加速技術(shù)的非結(jié)構(gòu)流域雨洪數(shù)值模型［J］．水科學(xué)進(jìn)展，2021，32（4）：567-576.

［17］冀永鵬，張洪興，王運(yùn)濤，等．基于二維淺水方程的城市地面洪水演進(jìn)數(shù)值模擬研究［J］．水資源與水工程學(xué)報(bào)，2020，31（2）：42-49，56.

［18］HUANG Y X，ZHANG N C，PEI Y G． Well-balanced finite volume scheme for shallow water flooding and drying over arbitrary topography［J］． Engineering Applications of Computational Fluid Mechanics，2013（7）：40-54.

［19］王秀杰，胡冰，苑希民，等．洪水與風(fēng)暴潮共同作用下的潰堤洪水一維、二維耦合模型及應(yīng)用［J］．南水北調(diào)與水利科技，2017，15（5）：43-49.

［20］李大鳴，羅珊，范麗虹，等．西三洼洪水演進(jìn)數(shù)值模擬及洪水風(fēng)險(xiǎn)分析［J］．水利水電技術(shù)，2018，49（8）：78-86.

［21］肖玉紅，李大鳴，白玲，等．蓄滯洪區(qū)洪水演進(jìn)過(guò)程數(shù)值模擬及洪災(zāi)損失分析［J］．中國(guó)農(nóng)村水利水電，2012（12）：143-147，149.

［22］王靜．基于 EFDC 的胖頭泡蓄滯洪區(qū)洪水演進(jìn)模擬研究［D］.大連：大連理工大學(xué)，2017.

（責(zé)任編輯：李燕珊）