風(fēng)暴潮三維數(shù)值計算模式的研究及在渤海灣的應(yīng)用

李大鳴, 范 玉,2, 徐亞男, 付慶軍, 李曉瑜

(1. 天津大學(xué) 建筑工程學(xué)院 港口與海洋工程教育部重點(diǎn)實驗室, 天津 300072; 2. 華北水利水電學(xué)院, 河南鄭州450011; 3. 交通運(yùn)輸部 天津水運(yùn)工程科學(xué)研究所, 天津 300456; 4. 國家海洋局 第一海洋研究所, 山東青島266061)

風(fēng)暴潮三維數(shù)值計算模式的研究及在渤海灣的應(yīng)用

李大鳴1, 范 玉1,2, 徐亞男3, 付慶軍4, 李曉瑜1

(1. 天津大學(xué) 建筑工程學(xué)院 港口與海洋工程教育部重點(diǎn)實驗室, 天津 300072; 2. 華北水利水電學(xué)院, 河南鄭州450011; 3. 交通運(yùn)輸部 天津水運(yùn)工程科學(xué)研究所, 天津 300456; 4. 國家海洋局 第一海洋研究所, 山東青島266061)

以實驗室二維溫帶風(fēng)暴潮數(shù)值模型為基礎(chǔ), 綜合考慮海洋潮波動力與風(fēng)應(yīng)力聯(lián)合作用, 建立溫帶風(fēng)暴潮三維數(shù)值計算模型。模型從推導(dǎo)三維風(fēng)暴潮基本控制方程出發(fā), 并應(yīng)用交替方向隱格式(ADI)方法對方程進(jìn)行離散求解。對于淺水動邊界, 模型采取局部深槽、縮小水域的活動邊界處理方法。利用擬三維數(shù)值計算方法, 并提出了非平面水深等分模式和平面等水深分布模式, 應(yīng)用這兩種計算模式分別對渤海灣2009年5月8~10日發(fā)生的風(fēng)暴潮過程進(jìn)行了數(shù)值模擬。將風(fēng)暴潮位計算結(jié)果和增水位計算結(jié)果與塘沽驗潮站的實際觀測數(shù)值進(jìn)行對比驗證, 結(jié)果顯示受風(fēng)應(yīng)力與潮波聯(lián)合作用的風(fēng)暴潮位和增水位與實測數(shù)據(jù)吻合良好; 通過比較得到了平面等水深分布模式的計算成果要比非平面水深等分模式的計算成果更接近觀測資料的結(jié)論, 為風(fēng)暴潮預(yù)報提供了理論依據(jù)。

風(fēng)暴潮; 數(shù)值模擬; 擬三維模型; 交替方向隱格式(ADI)方法; 渤海灣

風(fēng)暴潮的研究始于20世紀(jì)20年代, 最開始是對風(fēng)暴潮現(xiàn)象的成因進(jìn)行初步討論。從 50年代開始,風(fēng)暴潮數(shù)值預(yù)報的研究逐步發(fā)展起來, Hallsen[1]于1956年首次用電子計算機(jī)對歐洲北海的風(fēng)暴潮做了數(shù)值模擬, 并取得了初步成功, 但由于差分格式的不穩(wěn)定, 計算中出現(xiàn)了虛假的誤差短波。70年代, 日本的宮崎正衛(wèi)、宇野木早苗等[2]采用臺風(fēng)預(yù)報的多種模式計算結(jié)果相互對應(yīng)對照來進(jìn)行風(fēng)暴潮預(yù)報。美國提出了SPLASH(Special Program to List Amplitude of Surges from Hurricane)模式[3], 成為美國風(fēng)暴潮的業(yè)務(wù)預(yù)報模式, 至今仍為許多國家和地區(qū)使用。80年代后, 美國發(fā)布了新的風(fēng)暴潮數(shù)值預(yù)報模式SLOSH(Sea, Lake and Overland Surges from Hurricanes)模式[4], SLOSH模式是一個二維模式, 總是會低估淺水區(qū)域的風(fēng)暴潮程度。進(jìn)入90年代后, 風(fēng)暴潮與天文潮耦合的二維水流模型成為風(fēng)暴潮業(yè)務(wù)預(yù)報的主要手段, 同時風(fēng)暴潮三維水流數(shù)值模型也發(fā)展起來[5]。Jones用一個高分辨率的三維模型計算了1997年11月東愛爾蘭海的風(fēng)暴潮[6], 模式中還考慮了波浪和海流的相互作用。近年來, 風(fēng)暴潮運(yùn)動的三維數(shù)值計算模式得到很大發(fā)展, 如美國 Princeton大學(xué)開發(fā)的ECOM(Estuary Coast and Ocean Model)和POM(Princeton Ocean Model)模型、德國Hamburg大學(xué)開發(fā)的 HAMSOM(The Hydrodynamic Hamburg Shelf Ocean Model)陸架海洋動力模型、荷蘭Delft 科技大學(xué)的 Delft-3D模型、美國麻州大學(xué)海洋科學(xué)技術(shù)學(xué)院的 FVCOM(An Unstructured Grid Finite-Volume Coastal Ocean Model)模式等[7-11]。

我國在風(fēng)暴潮研究方面的工作開始得較晚, 進(jìn)入 70年代以后才逐漸開展了風(fēng)暴潮機(jī)制和預(yù)報的研究工作。1982 年科學(xué)出版社出版了馮士 筰[12]院士編著的國際上第一部風(fēng)暴潮專著——《風(fēng)暴潮導(dǎo)論》,孫文心[13]發(fā)表了國內(nèi)第一篇風(fēng)暴潮數(shù)值模擬的論文, 開創(chuàng)了國內(nèi)數(shù)值風(fēng)暴潮預(yù)報的先河。80年代以來, 我國風(fēng)暴潮的數(shù)值研究和應(yīng)用有了很大發(fā)展,對渤海、黃海、東海和南海陸架區(qū)以內(nèi)及臺灣海峽海域的風(fēng)暴潮數(shù)值模擬研究取得了較好的成果[14]。近年來, 風(fēng)暴潮運(yùn)動的三維數(shù)值模式得到很大的發(fā)展, 趙士清[15]采用與Leendertse類似的固定分層方法, 對長江口南槽和外海域的三維潮流進(jìn)行了數(shù)值模擬。竇振興[16]等采用σ坐標(biāo)系的三維模型和模態(tài)分離法對渤海灣潮流進(jìn)行了數(shù)值模擬。李艷蕓、李紹武[17]運(yùn)用 COHERHENS三維多功能大陸架水動力學(xué)模型建立了渤海灣天文潮數(shù)學(xué)模型。隨著計算技術(shù)的發(fā)展和三維水動力數(shù)學(xué)模型理論及計算方法的完善, 三維風(fēng)暴潮模型必將得到越來越廣泛的應(yīng)用。

本研究在實驗室二維風(fēng)暴潮數(shù)值預(yù)報模型的基礎(chǔ)上[18], 建立了三維風(fēng)暴潮數(shù)值計算模型。在模型中提出了兩種擬三維風(fēng)暴潮數(shù)值計算模式, 即非平面水深等分模式和平面等水深分布模式。利用該兩種模式對發(fā)生于渤海灣的一次風(fēng)暴潮過程進(jìn)行了數(shù)值模擬。計算得到天文潮與風(fēng)暴潮過程的流場流速、潮位與增水位值。將兩種模式的計算成果相互比較,表明平面等水深分布模式比非平面水深等分模式更接近實測數(shù)據(jù)。

1 擬三維計算模式的控制方程及模型的邊界條件

1.1 風(fēng)暴潮三維計算模式的控制方程

三維數(shù)值計算的控制方程是以靜水壓強(qiáng)假設(shè)下不可壓縮流體的三維流動的基本方程推導(dǎo)出的, 靜水壓強(qiáng)假設(shè)下不可壓縮流體的三維流動的基本方程如下:

連續(xù)方程:

運(yùn)動方程:

式中, u, v, w分別為x, y, z方向的速度分量; ρ為流體密度, g 為重力加速度, p 為壓強(qiáng)。τxx, τxy, τxz, τyx, τyy, τyz分別為流體受到的切應(yīng)力。f為柯氏(Coriolis)系數(shù),f= 2 ω sin φ, ω = 7 .29× 1 0-5rad/s為地球自轉(zhuǎn)的角速度, φ為緯度。對于大型水域如大洋, 往往需要考慮地球自轉(zhuǎn)引起的柯氏力。沿水深方向?qū)⒂嬎銋^(qū)域分為底層、中間層和頂層三層, 在每層中將不可壓縮流體的三維流動方程沿水深方向進(jìn)行積分, 積分后經(jīng)簡化得擬三維計算模式的控制方程如下:

底層連續(xù)方程:

底層運(yùn)動方程:

中間層連續(xù)方程:

中間層運(yùn)動方程:

頂層連續(xù)方程:

頂層運(yùn)動方程:

若H1,H2,H3為任意點(diǎn)處的底層、中間層、頂層靜水深。η,ξ,δ分別為底層、中間層、上層的增水位。D1,D2,D3任意點(diǎn)處的底層、中間層、頂層的厚度。有D1=H1+ ,D2=H2+ζ-η,D3=H3+δ-ζ。u1,v1,u2,v2,u3,v3分別為底層、中間層、底層的x,y方向的速度分量。γ1,γ2,γ3分別為底部、底層與中間層、中間層與頂層的摩擦系數(shù)。

1.2 模型的邊界條件及數(shù)值計算方法

岸邊界:vn= 0 (n為邊界法線方向)

初始條件: 當(dāng)t= 0時ξ=ξ0, 取ξ0=0;u1=v1=u2=v2=u3=v3= 0

應(yīng)用多分潮調(diào)和分析方法確定黃渤海海域水邊界條件, 然后利用所建立的三維數(shù)值模型計算黃渤海海域得出渤海灣海域的水邊界條件。在計算過程中采用了嵌套加密計算模式, 有限差分方法中的交替方向隱格式法(ADI)。ADI差分計算模式要求整個計算域應(yīng)保持在水深以下, 淺水岸邊界的露灘、淹沒變化應(yīng)是連續(xù)、穩(wěn)定過程。本文采用局部深槽、縮小水域的活動邊界處理方法, 當(dāng)全水深接近0.1 m時, 在淺水網(wǎng)格區(qū)全水深保持為H10=0.1 m, 保持流量、流速不變, 流量Q為:

Q=uH10Δs=u( |min(h,ξ)|+H10)Bs, 當(dāng)h+ξ≤H10時,變化后的水域?qū)挾菳s為

2 渤海灣風(fēng)暴潮數(shù)值模型的建立

本文關(guān)于風(fēng)暴潮的三維數(shù)值計算采用擬三維數(shù)值計算模式, 考慮兩種不同的計算模式, 即非平面水深等分模式和平面等水深分布模式。非平面水深等分模式是對整個水域的水深等距離地分層, 平面等水深分布模式是沿深度方向分層。非平面水深等分模式將整個計算區(qū)域平均分為三層, 在平面等水深分布模式中, 分層采用自由表面為第一層, 6 m深處為第二層, 12 m處為第三層。非平面水深等分模式及平面等水深分布模式的模型分層示意圖如圖 1所示。

圖1 兩種模式的分層示意圖Fig. 1 Schematic diagram of the two model delamination

3 三維風(fēng)暴潮數(shù)學(xué)模型在渤海灣的應(yīng)用

3.1 非平面水深等分模式的應(yīng)用

應(yīng)用非平面水深等分三維數(shù)值計算模式, 對渤海灣2009年5月8~10日發(fā)生的風(fēng)暴潮過程進(jìn)行數(shù)值模擬驗證, 計算得到天文潮、受潮波動力及風(fēng)應(yīng)力影響的風(fēng)暴潮流場, 潮位及增水過程。計算所得的流場圖如圖2、圖3所示, 其中圖2為天文潮作用下的各層流場示意圖, 圖 3為天文潮與風(fēng)應(yīng)力共同作用下的各層流場示意圖。

采用塘沽驗潮站的實測潮位資料對渤海灣計算潮位進(jìn)行驗證, 天文潮潮位過程、風(fēng)暴潮潮位及增水驗證曲線如圖4所示, 分析圖 4曲線趨勢可以看出,曲線擬合度較好, 計算風(fēng)暴潮峰值潮位為 5.33 m,而實測峰值潮位4.95 m, 相對誤差為7.6%。增水驗證效果良好, 在增水極值處, 計算值與實測值較為接近, 相對誤差4.3%。

3.2 平面等水深分布模式

利用平面等水深分布三維數(shù)值計算模式對渤海灣風(fēng)暴潮進(jìn)行模擬計算, 計算所得的潮流流場示意圖如圖5、圖6所示。其中圖5為天文潮作用下各層流場示意圖, 圖 6為天文潮與風(fēng)應(yīng)力共同作用下各層流場示意圖。

采用塘沽驗潮站的實測資料進(jìn)行潮位驗證, 得到平面等水深分布模式下的天文潮、風(fēng)暴潮及增水驗證曲線如圖 7所示, 計算曲線與實測曲線擬合良好。在風(fēng)暴潮曲線的潮位峰值處, 計算值為 4.72 m,實測數(shù)據(jù)是4.95 m, 相對誤差僅為4.5%。計算最大增水位1.3 m, 與實測最大增水的相對誤差是2.6%。表 1給出了兩種模式的風(fēng)暴潮計算值與測量值的數(shù)據(jù)比較。

圖2 非平面水深等分模式下天文潮流場圖Fig. 2 Flow field of astronomic tide under non-planar depth division mode

圖3 非平面水深等分模式下受潮波與風(fēng)應(yīng)力共同作用的風(fēng)暴潮流場圖Fig. 3 Flow field of storm surge by astronomic tide and wind stress under non-planar depth division mode

圖 4 非平面水深等分模式下天文潮、風(fēng)暴潮潮位及增水驗證Fig. 4 Tidal level verification of astronomic tide,storm surge and storm surge setup verification under non-planar depth division mode

圖5 平面等水深分布模式下天文潮作用下流場圖Fig. 5 Flow field of astronomic tide under equal depth distribution mode

圖6 平面等水深分布模式下天文潮與風(fēng)應(yīng)力共同作用下風(fēng)暴潮流場圖Fig. 6 Flow field of storm surge by astronomic tide and wind stress under equal depth distribution mode

圖 7 平面等水深分布模式下天文潮、風(fēng)暴潮潮位及增水驗證Fig. 7 Tidal level verification of astronomic tide,storm surge and storm surge setup verification under equal depth distribution mode

表1 兩種計算模式數(shù)據(jù)對比Tab. 1 Data comparison of the two calculation mode

4 結(jié)論

本文基于實驗室的二維溫帶風(fēng)暴潮數(shù)學(xué)模型,建立了三維風(fēng)暴潮數(shù)值預(yù)報模型。模型中考慮海洋潮波動力與風(fēng)應(yīng)力的聯(lián)合作用, 從 N-S方程出發(fā)推導(dǎo)了三維風(fēng)暴潮控制方程, 并應(yīng)用 ADI離散方程求解, 保持了計算的穩(wěn)定性。模型對淺水動邊界采取局部深槽、縮小水域的活動處理方法, 使得邊界能夠進(jìn)行連續(xù)穩(wěn)定的計算。利用兩種不同的三維數(shù)值計算模式, 分別對渤海灣2009年5月8~10日發(fā)生的風(fēng)暴潮過程進(jìn)行數(shù)值模擬。將計算得到的風(fēng)暴潮潮位計算結(jié)果和增水位計算結(jié)果與塘沽驗潮站的實際觀測資料進(jìn)行了對比分析。首先, 結(jié)合曲線的相關(guān)性進(jìn)行分析, 兩種計算模式均顯示良好的潮位曲線擬合度。其次, 分析潮位與增水峰值的相對誤差, 結(jié)果顯示平面等水深分布模式的計算成果要優(yōu)于非平面水深等分模式的計算成果。本文系統(tǒng)分析了形成這一結(jié)論的原因, 兩種模式均受數(shù)值計算模式本身、近岸水深和岸邊界動態(tài)變化的影響, 使得平面等水深分布模式的計算成果要比非平面水深等分模式的計算成果更接近實際觀測資料。

[1] Hansen W.Theorie zur Errechnuang des Wasserstancles und der Stromungen in randmeeren Nebst Anwendungen [J]. Tellus, 1956, 8(3): 287-300.

[2] 宮崎正衛(wèi), Jelesnianski C P.臺風(fēng)風(fēng)暴潮預(yù)報技術(shù)手冊[M]. 北京: 海洋出版社, 1982.

[3] Jelesnianski C P. A numerical calculation of storm tides induced by a tropical storm impinging on a continental shelf [J]. Mon Wea Rev, 1965, 93: 343-358.

[4] Jelesnianski C P, Chen J, Wilson A. SLOSH: Sea, Lake and ovedand Surges from Hurricanes[R]. NOAA: NWS 48, 1992.

[5] Yamashita T, Yoshito T, Yoshioka H. Quasi-threedimensional model for storm surges and its verification[J]. Coastal Engineering, 1994, 901-915.

[6] Jones J E, Davies A M. Influence of wave-current interaetion, and high frequeney forcing upon storm induced currents and elevations [J]. Estuarine, Coastal and ShelfScience, 2001, 53: 397-413.

[7] 湯立群, 申錦瑜, 劉大濱, 等. 中國河口海岸風(fēng)暴潮及海洋動力三維數(shù)值預(yù)報模型(CHINACOAST)研究I: 模型結(jié)構(gòu)與設(shè)置[J]. 水利水電技術(shù), 2009, 12(40): 1-7.

[8] 馬鋼峰, 劉曙光, 戚定滿. 長江口鹽水入侵?jǐn)?shù)值模型研究[J]. 水動力學(xué)研究與進(jìn)展, 2006, 1(21): 53-61.

[9] 姚遠(yuǎn), 蔡樹群, 王盛安. 臺灣海峽一次海嘯的初步數(shù)值模擬[J]. 熱帶海洋學(xué)報, 2009, 2(28): 1-6.

[10] 姬厚德. 廈門灣臺風(fēng)浪場數(shù)值模擬[J]. 海洋地質(zhì)動態(tài), 2009, 6(25): 15-19.

[11] 李希彬, 鮑獻(xiàn)文, 馬超, 等. 東海大堤對湛江灣水動力環(huán)境影響的研究[J]. 中國海洋大學(xué)學(xué)報, 2009, 39:287-296.

[12] 筰馮士. 風(fēng)暴潮導(dǎo)論[M]. 北京: 科學(xué)出版社, 1982:1-9.

[13] 孫文心, 馮士筰, 秦曾灝. 超淺海風(fēng)暴潮的數(shù)值模擬(I)[J]. 海洋學(xué)報, 1979, 1(2): 193-211.

[14] 孫英蘭, 張越美. 膠州灣三維變動邊界的潮流數(shù)值模擬[J]. 海洋與湖沼, 2001, 32(4): 355-362.

[15] 趙士清. 長江口三維潮流的數(shù)值模擬[J]. 水利水運(yùn)研究, 1985, 1: 8-20.

[16] 竇振興, 楊連武, Ozer J. 渤海三維潮流數(shù)值模擬[J].海洋學(xué)報, 1993, 15(5): 1-15.

[17] 李艷蕓, 李紹武. 風(fēng)暴潮預(yù)報模式理論及應(yīng)用研究[D]. 天津: 天津大學(xué),2005.

[18] 李大鳴, 徐亞男, 宋雙霞, 等. 波浪輻射應(yīng)力在渤海灣海域?qū)︼L(fēng)暴潮影響的研究[J]. 水動力學(xué)研究與進(jìn)展A輯, 2010, 25(3): 374-382.

Research on 3-D numerical calculation model of storm surge and application in the Bohai Bay

LI Da-ming1, Fan Yu1,2, XU Ya-nan3, FU Qing-jun4, LI Xiao-yu1
(1. Key Laboratory of Harbor and Ocean Engineering Ministry of Education, School of Civil Engineering, Tianjin University, Tianjin 300072, China; 2.North China University of Water Resources and Electric Power, Zhengzhou 450011, China 3.Tianjin Research Institute for Water Transport Engineering, Ministry of Transport of the Peple’s Republic of China, Tianjin 300456, China; 4. The first Institute of Oceangraphy, State Oceanic Administration,Qingdao 266061, China)

Feb.,17,2011

storm surge; numerical imitation; quasi-three-dimensional model; ADI (Alternating Direction Implicit) method; the Bohai Bay

This study is based on the 2-D numerical model of extratropical storm surge in our laboratory. The astronomical tidal force and wind stress were considerd as the main factors when we built the 3-D numerical model of storm surge. The model derives three-dimensional storm surge mathematical model equations and applies ADI(Alternating Direction Implicit) method to discretize and solve governing equations. To deal with the lateral boundary, the method of local-deepening and water range reducing were considered in the model. In the numerical calculation of extratropical storm surge, the numerical model of quasi-three-dimensional was applied, and then non-planar depth division model and planar equal depth distribution model were introduced and used in Bohai Bay for the computing of storm surge which occurred May 8 to 10, 2009. The comparison between the observed value of Tanggu tidal station and the calculated results of storm surge tidal level and setup level showed that the computed values which affected by astronomical tidal force and the wind stress are in good agreement with the measured data.According to the analysis, the results of the planar equal depth distribution model are closer to the practical situation than those of the non-planar depth division model, and it supplies theoretical basis for storm surge forecasting.

P753

A

1000-3096(2012)07-0007-07

2011-02-17;

2012-04-20

國家自然科學(xué)基金資助項目(51079095); 國家自然科學(xué)基金創(chuàng)新研究群體科學(xué)基金(51021004)

李大鳴(1957-), 男, 天津人, 博士, 教授, 博士生導(dǎo)師, 主要從事水力學(xué)及河流動力學(xué)、海岸工程, 電話: 022-87401579,E-mail: lidaming@tju.edu.cn

(本文編輯：劉珊珊)