基于VOF法的風(fēng)暴潮洪水演進(jìn)三維數(shù)值模擬

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

(1.天津大學(xué)水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室，天津 300072;2.天津市水利勘測設(shè)計(jì)院，天津 300204)

1 前言

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

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

2 風(fēng)暴潮洪水演進(jìn)三維數(shù)學(xué)模型

2.1 控制方程

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

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

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

連續(xù)性方程:

對于第q相流體，體積分?jǐn)?shù)連續(xù)方程為:

動(dòng)量方程:

k方程:

ε方程:

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

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

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

當(dāng)氣相是紊流時(shí):

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

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

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

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

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

在淹沒出流條件下的計(jì)算公式為:

式(11)中，H2為淹沒區(qū)水位高出堰頂?shù)母叨?n為指數(shù);h下為淹沒區(qū)水位;p為淹沒區(qū)內(nèi)堰高[12]。

2.2 邊界條件

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

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

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

4)壓力邊界條件:在計(jì)算域的邊界上，壓力應(yīng)滿足Neumann條件。為保證計(jì)算的穩(wěn)定性，在規(guī)定的某一內(nèi)點(diǎn)上，壓力為給定值，定義為參考壓力pref。

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

3 模型驗(yàn)證

采用意大利ENEL機(jī)構(gòu)Toce河物理模型檢驗(yàn)數(shù)學(xué)模型的精確性和穩(wěn)定性[14]。選用平行式城區(qū)潰壩洪水試驗(yàn)資料對數(shù)學(xué)模型進(jìn)行驗(yàn)證，計(jì)算模型共劃分為8 904個(gè)網(wǎng)格，其中網(wǎng)格長X×Y×Z=1 m×1 m ×0.01 m和X×Y×Z=1 m ×1 m ×0.05 m，如圖1所示。物理模型區(qū)域內(nèi)共布置了10個(gè)觀測點(diǎn)，測點(diǎn)分布位置如圖2所示。將具有代表性的4個(gè)測點(diǎn)處實(shí)測水深與模型模擬水深進(jìn)行對比，如圖3所示。

圖1 模型網(wǎng)格劃分Fig.1 Mesh generation of the model

圖2 測點(diǎn)分布圖Fig.2 The distribution of measuring points

圖3 各個(gè)測點(diǎn)實(shí)測水深與模擬水深對比Fig.3 The comparison of measured depth and simulated depth at every test point

由圖3可知，總體上實(shí)測水深值與模擬值符合較好。計(jì)算結(jié)果表明，該數(shù)學(xué)模型可以用于模擬風(fēng)暴潮洪水演進(jìn)，具有較好的精確性和穩(wěn)定性。

4 風(fēng)暴潮洪水演進(jìn)數(shù)值模擬結(jié)果分析

4.1 網(wǎng)格劃分

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

4.2 風(fēng)暴潮洪水演進(jìn)數(shù)值模擬

建立了三維風(fēng)暴潮洪水演進(jìn)數(shù)學(xué)模型，對天津市濱海新區(qū)100年一遇風(fēng)暴潮歷時(shí)4 h的過程進(jìn)行了數(shù)值模擬，得到的研究區(qū)域中不同時(shí)刻VOF的分布如圖5所示。

圖4 風(fēng)暴潮洪水演進(jìn)數(shù)值模擬研究范圍Fig.4 Research area of storm flood routing surge simulation

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

由圖5可知，2 300 s時(shí)，東疆港區(qū)基本被洪水淹沒;3 600 s時(shí)，南疆港區(qū)洪水演進(jìn)基本穩(wěn)定;20 560 s時(shí)，演進(jìn)范圍達(dá)到最大;32 000 s時(shí)，除各壅水處外，洪水基本消退。

4.3 不同頻率風(fēng)暴潮洪水嚴(yán)重性分析

4.3.1 不同頻率風(fēng)暴潮最大淹沒范圍分析

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

4.3.2 不同頻率風(fēng)暴潮特征水深分析

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

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

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

4.3.3 不同頻率風(fēng)暴潮下典型地點(diǎn)水深分析

選取淹沒范圍內(nèi)的天津港客運(yùn)站、保稅區(qū)和海事法院3個(gè)典型地點(diǎn)對不同頻率風(fēng)暴潮下洪水水深進(jìn)行分析，3個(gè)地點(diǎn)的位置如圖4所示。

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

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

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

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

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

圖8 3種頻率風(fēng)暴潮下天津港客運(yùn)站在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種頻率風(fēng)暴潮下保稅區(qū)在最大水深時(shí)的水深分布，3種頻率風(fēng)暴潮下的水深具有相同的分布趨勢:保稅區(qū)的西南部分洪水最深，而西北部分水深最淺，洪水由南向北逐漸變淺。50年一遇、100年一遇和200年一遇風(fēng)暴潮下保稅區(qū)的最大平均水深分別為0.92 m、0.99 m 和1.05 m，水深隨著頻率的增加逐漸減少。

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

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

5 結(jié)語

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

[1]李 鑫.渤海風(fēng)暴潮增減水及流場時(shí)空特征初步研究[D].南京:南京水利科學(xué)研究院，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，F(xiàn)rance，2006.

[6]張大偉.堤壩潰決水流數(shù)學(xué)模型及其應(yīng)用研究[D].北京:清華大學(xué)，2008.

[7]姚志堅(jiān)，陳文龍.潰壩洪水演進(jìn)數(shù)學(xué)模型研究及其在大鏡山水庫在的應(yīng)用[J].人民珠江，2008(3):7-9.

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

[9]黃金池，何曉燕.潰壩洪水的統(tǒng)一二維數(shù)學(xué)模型[J].水利學(xué)報(bào)，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]王曉玲，段琦琦，佟大威，等.長距離無壓引水隧洞水氣兩相流數(shù)值模擬[J].水利學(xué)報(bào)，2009，40(5):596-602.

[12]張行南，張文婷，劉永志，等.風(fēng)暴潮洪水淹沒計(jì)算模型研究[J].系統(tǒng)仿真學(xué)報(bào)，2006，18(增刊2):20-23.

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

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