基于VOF法的風暴潮洪水演進三維數值模擬

王曉玲，孫宜超，陳華鴻，楊麗美，孫蕊蕊

(1.天津大學水利工程仿真與安全國家重點實驗室，天津 300072;2.天津市水利勘測設計院，天津 300204)

1 前言

風暴潮是指由于大風以及氣壓急劇變化等因素造成的沿海或河口水位的異常升降現象。風暴潮往往導致沿海地區水位暴漲，造成嚴重的淹沒損失。風暴潮數值模擬的研究開始于20世紀50年代[1]，國外一些學者采用數值模擬的手段研究了密集城區的潰堤洪水問題，Haider等和Inoue等采用二維淺水方程分別模擬了法國Nimes城和日本Osaka城的潰堤大洪水，較好地模擬了阻水區域和蓄水區域，均取得了理想的計算結果[2，3];Mignot等和 Soares 等通過城區潰壩物理模型試驗對二維淺水方程數學模型的精確性進行了驗證，通過對試驗數據和模擬數據的比較，發現洪水淹沒水深、流速和演進時間均較吻合，然而由于地形條件的不確定性和數學模型模擬的局限性，仍然存在一些偏差[4，5]。對于密集城區內潰堤洪水計算方面，國內涉及較少。張大偉等在堤壩潰決水流數學模型中，采用真實地形法對建筑物進行處理，較真實地模擬了建筑物對潰決水流產生的影響[6];姚志堅等提出在潰壩洪水演進計算中用“等效糙率”模擬建筑群的方法及“等效糙率”取值的水槽試驗手段，并成功應用在某城市水庫潰壩洪水演進研究中[7];張大偉通過考慮社區和樓房內部的容水性，引入侵入水量的概念并給出其合理的估算方法，建立了能夠模擬建筑物密集城區潰堤水流運動的二維淺水方程數學模型，對哈爾濱市可能發生的潰堤洪水進行了計算研究[8]。

以上研究多是針對水庫潰壩后洪水的數值模擬，而且多是一、二維[9]，國內外基于風暴潮洪水演進三維數值模擬的相關研究成果較少，缺乏對含有復雜地貌的大范圍洪水演進的數值模擬。文章建立了耦合流體體積函數(volume of fluid，VOF)法的三維非穩態水氣兩相流k-ε模型，采用等效糙率的方法處理城市密集建筑群，重點對天津市濱海新區海河與永定新河之間區域100年一遇的風暴潮洪水淹沒情況進行了數值模擬與分析，并對不同頻率的風暴潮洪水的嚴重性做了比較分析。

2 風暴潮洪水演進三維數學模型

2.1 控制方程

VOF方法由 Hirt和 Nichols提出[10]，是一種處理自由表面的有效方法，模型對每一項引入體積分數，通過求解每一控制單元內的體積分數，確定相間界面。用aq表示單元內第q相流體的體積分數，如果aq=1，說明該單元內僅含第q相流體;如果aq=0，說明該單元內不含第q相流體;如果0＜aq＜1，說明該單元包含部分該流體，也即該單元內包含第q相流體與其他流體的交接面。在每一個控制體中，所有的相體積分數之和為1，即:

式(2)中，u、v、w分別為主場沿x、y、z方向的分速度;Saq為該相的質量源，無源情況下為零;ρq為該相流體的密度，kg/m3。

VOF模型采用三維k-ε紊流模型，基本控制方程如下:

連續性方程:

對于第q相流體，體積分數連續方程為:

動量方程:

k方程:

ε方程:

式(3)～(6)中，t為時間，s;ui和uj為速度分量，m/s;xi和xj為坐標分量，m;ρ為密度，kg/m3;μ為分子動力粘性系數，N·m/s;k為紊動動能，m2/s2;G為紊動能生成率;ε為紊動耗散率，m2/s2;p為修正壓力，Pa;μt為紊流粘性系數;σk、σε分別為k、ε的紊流普朗特數，無因次;C1ε、C2ε為經驗常數，無因次;控制方程中的常數值 σk=1.0，σε=1.3，C1ε=1.44，C2ε=1.92[11]。

氣相單位面積上壁面摩阻力 g按式(7)計算:

式(7)中，ρg為氣相密度，kg/m3;μg為氣相斷面平均流速，m/s;fg為氣相壁面摩阻系數，無量綱;Qg為氣相斷面流量，m3/s;Ag為氣相斷面面積，m2。

當氣相是紊流時:

式(8)中，χg為氣相斷面的濕周，m;BL為水相的斷面寬度，m;νg為氣相運動粘滯系數，m2/s。

氣水界面單位面積上的摩阻力 i按式(9)計算:

式(9)中，uL為水相斷面平均流速，m/s;fi為水氣界面摩阻系數，無量綱;QL為水相斷面流量，m3/s;AL為水相斷面面積，m2。

進口流量過程分為兩種情況進行計算，在自由出流條件計算的計算公式為:

式(10)中，Q為堰前水頭為H時的自由出流流量;μ為潰堤流量系數，對于無坎寬頂堰取值為0.2～0.3;B為漫灘岸線長度，m;H為堰上水頭，即外部潮位與灘地邊緣高程的高度差，m。

在淹沒出流條件下的計算公式為:

式(11)中，H2為淹沒區水位高出堰頂的高度;n為指數;h下為淹沒區水位;p為淹沒區內堰高[12]。

2.2 邊界條件

1)進口邊界條件:進口流量按照式(10)和式(11)確定，k和ε根據經驗公式給出[13]。

2)出口邊界條件:通常在計算域的出口，各速度分量(u，v，w)以及k和ε均取為第二類齊次邊界條件。

3)固壁邊界條件:在壁面上采用無滑移條件，即:u=v=w=0。為計算近壁區的紊流，文章采用壁面函數法[12]。

4)壓力邊界條件:在計算域的邊界上，壓力應滿足Neumann條件。為保證計算的穩定性，在規定的某一內點上，壓力為給定值，定義為參考壓力pref。

5)城市建筑物:在處理建筑群糙率的問題上，將多個建筑物作為一個整體，采用等效糙率的方法來模擬其作用。

3 模型驗證

采用意大利ENEL機構Toce河物理模型檢驗數學模型的精確性和穩定性[14]。選用平行式城區潰壩洪水試驗資料對數學模型進行驗證，計算模型共劃分為8 904個網格，其中網格長X×Y×Z=1 m×1 m ×0.01 m和X×Y×Z=1 m ×1 m ×0.05 m，如圖1所示。物理模型區域內共布置了10個觀測點，測點分布位置如圖2所示。將具有代表性的4個測點處實測水深與模型模擬水深進行對比，如圖3所示。

圖1 模型網格劃分Fig.1 Mesh generation of the model

圖2 測點分布圖Fig.2 The distribution of measuring points

圖3 各個測點實測水深與模擬水深對比Fig.3 The comparison of measured depth and simulated depth at every test point

由圖3可知，總體上實測水深值與模擬值符合較好。計算結果表明，該數學模型可以用于模擬風暴潮洪水演進，具有較好的精確性和穩定性。

4 風暴潮洪水演進數值模擬結果分析

4.1 網格劃分

以海河與永定新河之間范圍確定計算區域，根據濱海新區地形圖劃分出風暴潮洪水演進計算區域，其最大長度為42.28 km，最大寬度為26.78 km，區域面積約為525 km2，計算區域如圖4(陰影區域)所示。

4.2 風暴潮洪水演進數值模擬

建立了三維風暴潮洪水演進數學模型，對天津市濱海新區100年一遇風暴潮歷時4 h的過程進行了數值模擬，得到的研究區域中不同時刻VOF的分布如圖5所示。

圖4 風暴潮洪水演進數值模擬研究范圍Fig.4 Research area of storm flood routing surge simulation

圖5 研究區域不同時刻VOF分布Fig.5 The VOF distribution of different time at the research area

由圖5可知，2 300 s時，東疆港區基本被洪水淹沒;3 600 s時，南疆港區洪水演進基本穩定;20 560 s時，演進范圍達到最大;32 000 s時，除各壅水處外，洪水基本消退。

4.3 不同頻率風暴潮洪水嚴重性分析

4.3.1 不同頻率風暴潮最大淹沒范圍分析

圖6為不同頻率風暴潮最大淹沒面積分析。50年一遇風暴潮潮位為5.83 m，100年一遇風暴潮潮位為6.01 m，200年一遇風暴潮潮位為6.19 m。隨著風暴潮發生頻率的增加，風暴潮淹沒面積逐漸減小，且200年一遇和100年一遇風暴潮的最大淹沒面積相差較大，而100年一遇和50年一遇風暴的最大淹沒面積相差比較小，這說明雖然200年一遇、100年一遇和50年一遇的風暴潮潮位均差0.18 m，但是200年一遇相對于100年一遇和50年一遇風暴潮來說，對濱海新區造成的損失將會大大增加。

4.3.2 不同頻率風暴潮特征水深分析

表1為不同頻率風暴潮最大淹沒范圍下的特征水深。3種典型風暴潮的特征水深基本沒有差別，這主要是由于3種風暴潮的潮位本來就差別不大，而且200年一遇、100年一遇和50年一遇風暴潮的淹沒范圍是逐漸減小的，所以造成了各特征水深差別不大。

圖6 不同頻率風暴潮最大淹沒范圍Fig.6 The maximum flooded area at different frequencies of the storm surge

表1 不同頻率風暴潮最大淹沒范圍下的水深Table 1 The water depth of the maximum flooded area at different frequencies of the storm surge m

4.3.3 不同頻率風暴潮下典型地點水深分析

選取淹沒范圍內的天津港客運站、保稅區和海事法院3個典型地點對不同頻率風暴潮下洪水水深進行分析，3個地點的位置如圖4所示。

圖7為典型地點不同頻率風暴潮下洪水水深隨時間的變化圖。天津港客運站在0～1 h時，3種頻率風暴潮下洪水水深隨時間逐漸增加。在1～4 h時，洪水在天津港客運站達到最大且趨于穩定;在4 h之后，風暴潮結束，洪水回流海洋，水深急劇下降，最后降為零。保稅區和天津港客運站有相同的水深變化趨勢，只是天津港客運站比保稅區更靠近海岸，所以天津港客運站水深比保稅區水深增加得快:在0～1.5 h 時，水深逐漸增大;在1.5 ～4 h時洪水水深達到最大;4 h后洪水退卻，水深逐漸變淺，最后趨于零。海事法院在0～1.5 h時水深為零;在1.5 ～3.5 h 內，水深逐漸增加;在 3.5 ～5.5 h 時段內水深趨于穩定并達到最大。5.5 h后由于風暴潮結束，沒有后續洪水供應，水深有所下降。

圖7 典型地點不同頻率風暴潮下洪水水深分析Fig.7 The water depth analysis of typical regions at different frequencies of the storm surge

天津港客運站和保稅區在不同頻率風暴潮下的最大淹沒水深有明顯差別，但在水位的上升和下降階段差別不大;天津港客運站緊挨著海岸，保稅區距離海岸也比較近，當風暴潮結束后，洪水開始回流海洋，并最終趨于零。海事法院位于內陸，距海岸相對較遠，所以風暴潮發生一段時間后才被淹沒，水深逐漸增加，并最后趨于穩定值;當風暴潮結束后，洪水積存于原處，并沒有流回海洋。

由于風暴潮在水深最大的時候所造成的損失是最大的，所以有必要分析三地在最大平均水深時的水深分布。

圖8為3種頻率風暴潮下天津港客運站在最大水深時的水深分布，x方向為實際地圖中正東方向，y方向為正北方向。3種頻率風暴潮下的水深具有相同的分布趨勢。天津港客運站的西南部分洪水最深，而西北部分水深最淺，洪水由南向北逐漸變淺。50年一遇、100年一遇和200年一遇風暴潮下天津港客運站的最大平均水深分別為1.13 m、1.21 m和1.26 m，水深隨著頻率的增加是逐漸減少的。

圖8 3種頻率風暴潮下天津港客運站在3 h的水深分布Fig.8 The water depth distribution at 3 h of Tianjin Port passenger station at different frequencies of the storm surge

圖9為3種頻率風暴潮下保稅區在最大水深時的水深分布，3種頻率風暴潮下的水深具有相同的分布趨勢:保稅區的西南部分洪水最深，而西北部分水深最淺，洪水由南向北逐漸變淺。50年一遇、100年一遇和200年一遇風暴潮下保稅區的最大平均水深分別為0.92 m、0.99 m 和1.05 m，水深隨著頻率的增加逐漸減少。

3種頻率風暴潮下海事法院在最大水深時的水深分布同天津港客運站和保稅區一樣，具有相同的分布趨勢，水深隨著頻率的增加逐漸減少。

圖9 3種頻率風暴潮下保稅區在3 h的水深分布Fig.9 The water depth distribution at 3 h of Bonded Area at different frequencies of the storm surge

5 結語

文章建立了耦合VOF法的三維非穩態水氣兩相流k-ε模型，在處理建筑群糙率的問題上，將多個建筑物作為一個整體，采用等效糙率的方法處理城市密集建筑群，既考慮了其阻水作用，又考慮了其蓄水作用，針對天津市濱海新區海河與永定新河之間區域的風暴潮洪水演進數值模擬與分析，對100年一遇風暴潮洪水淹沒情況進行分析，并對不同頻率的風暴潮洪水的嚴重性進行了比較。結果表明，隨著風暴潮發生頻率的增加，風暴潮淹沒范圍逐漸減小;水深隨著頻率的增加逐漸減少。研究為海堤安全管理、風暴潮災害的快速科學評估提供了理論依據和技術支持。

[1]李 鑫.渤海風暴潮增減水及流場時空特征初步研究[D].南京:南京水利科學研究院，2007.

[2]Haider S，Paquier A，Morel R，et al.Urban flood modeling using computational fluid dynamics[J].Water and Maritime Engineering，2003，156:129-135.

[3]Inoue K，Kawaike K，Hayashi H.Numerical simulation models on inundation flow in urban area[J].Journal of Hydroscience and Hydraulic Engineering，2000，18(1):119-126.

[4]Mignot E，Paquier A，Ishigaki T.Comparison of numerical and experimental simulations of a flood in a dense urban area[J].Water Science and Technology，2006，54(6-7):65-73.

[5]Soares S，Spinewine B，Duthoit A，et al. Dam-break flow experiments in simplified city layouts[C]∥7thInternational Conference on Hydro-informatics，Nice，France，2006.

[6]張大偉.堤壩潰決水流數學模型及其應用研究[D].北京:清華大學，2008.

[7]姚志堅，陳文龍.潰壩洪水演進數學模型研究及其在大鏡山水庫在的應用[J].人民珠江，2008(3):7-9.

[8]張大偉，程曉陶，黃金池.建筑物密集城區潰堤水流二維數值模擬[J].水利學報，2010，41(3):272-277.

[9]黃金池，何曉燕.潰壩洪水的統一二維數學模型[J].水利學報，2007，37(2):222-226.

[10]Hirt C W，Nichols B D.Volume of fluid(VOF)method for the dynamics of free boundaries[J].Journal of Computational Physics，1981，39(1):201-225.

[11]王曉玲，段琦琦，佟大威，等.長距離無壓引水隧洞水氣兩相流數值模擬[J].水利學報，2009，40(5):596-602.

[12]張行南，張文婷，劉永志，等.風暴潮洪水淹沒計算模型研究[J].系統仿真學報，2006，18(增刊2):20-23.

[13]周曉泉.復雜邊界條件下的三維紊流數值模擬研究[D].成都:四川大學，2003.

[14]Guido Testa，David Zuccalà，Francisco Alcrudo，et al.Flash flood flow experiment in a simplified urban district[J].Journal of Hydraulic Research，2007，45:37-44.