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

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

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

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

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

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

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

風(fēng)暴潮的研究始于20世紀(jì)20年代, 最開始是對風(fēng)暴潮現(xiàn)象的成因進(jìn)行初步討論。從 50年代開始,風(fēng)暴潮數(shù)值預(yù)報(bào)的研究逐步發(fā)展起來, Hallsen[1]于1956年首次用電子計(jì)算機(jī)對歐洲北海的風(fēng)暴潮做了數(shù)值模擬, 并取得了初步成功, 但由于差分格式的不穩(wěn)定, 計(jì)算中出現(xiàn)了虛假的誤差短波。70年代, 日本的宮崎正衛(wèi)、宇野木早苗等[2]采用臺風(fēng)預(yù)報(bào)的多種模式計(jì)算結(jié)果相互對應(yīng)對照來進(jìn)行風(fēng)暴潮預(yù)報(bào)。美國提出了SPLASH(Special Program to List Amplitude of Surges from Hurricane)模式[3], 成為美國風(fēng)暴潮的業(yè)務(wù)預(yù)報(bào)模式, 至今仍為許多國家和地區(qū)使用。80年代后, 美國發(fā)布了新的風(fēng)暴潮數(shù)值預(yù)報(bào)模式SLOSH(Sea, Lake and Overland Surges from Hurricanes)模式[4], SLOSH模式是一個(gè)二維模式, 總是會低估淺水區(qū)域的風(fēng)暴潮程度。進(jìn)入90年代后, 風(fēng)暴潮與天文潮耦合的二維水流模型成為風(fēng)暴潮業(yè)務(wù)預(yù)報(bào)的主要手段, 同時(shí)風(fēng)暴潮三維水流數(shù)值模型也發(fā)展起來[5]。Jones用一個(gè)高分辨率的三維模型計(jì)算了1997年11月東愛爾蘭海的風(fēng)暴潮[6], 模式中還考慮了波浪和海流的相互作用。近年來, 風(fēng)暴潮運(yùn)動的三維數(shù)值計(jì)算模式得到很大發(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ù)報(bào)的研究工作。1982 年科學(xué)出版社出版了馮士 筰[12]院士編著的國際上第一部風(fēng)暴潮專著——《風(fēng)暴潮導(dǎo)論》,孫文心[13]發(fā)表了國內(nèi)第一篇風(fēng)暴潮數(shù)值模擬的論文, 開創(chuàng)了國內(nèi)數(shù)值風(fēng)暴潮預(yù)報(bào)的先河。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é)模型。隨著計(jì)算技術(shù)的發(fā)展和三維水動力數(shù)學(xué)模型理論及計(jì)算方法的完善, 三維風(fēng)暴潮模型必將得到越來越廣泛的應(yīng)用。

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

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

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

三維數(shù)值計(jì)算的控制方程是以靜水壓強(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ū)⒂?jì)算區(qū)域分為底層、中間層和頂層三層, 在每層中將不可壓縮流體的三維流動方程沿水深方向進(jìn)行積分, 積分后經(jīng)簡化得擬三維計(jì)算模式的控制方程如下:

底層連續(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ù)值計(jì)算方法

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

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

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

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

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

本文關(guān)于風(fēng)暴潮的三維數(shù)值計(jì)算采用擬三維數(shù)值計(jì)算模式, 考慮兩種不同的計(jì)算模式, 即非平面水深等分模式和平面等水深分布模式。非平面水深等分模式是對整個(gè)水域的水深等距離地分層, 平面等水深分布模式是沿深度方向分層。非平面水深等分模式將整個(gè)計(jì)算區(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ù)值計(jì)算模式, 對渤海灣2009年5月8~10日發(fā)生的風(fēng)暴潮過程進(jìn)行數(shù)值模擬驗(yàn)證, 計(jì)算得到天文潮、受潮波動力及風(fēng)應(yīng)力影響的風(fēng)暴潮流場, 潮位及增水過程。計(jì)算所得的流場圖如圖2、圖3所示, 其中圖2為天文潮作用下的各層流場示意圖, 圖 3為天文潮與風(fēng)應(yīng)力共同作用下的各層流場示意圖。

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

3.2 平面等水深分布模式

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

采用塘沽驗(yàn)潮站的實(shí)測資料進(jìn)行潮位驗(yàn)證, 得到平面等水深分布模式下的天文潮、風(fēng)暴潮及增水驗(yàn)證曲線如圖 7所示, 計(jì)算曲線與實(shí)測曲線擬合良好。在風(fēng)暴潮曲線的潮位峰值處, 計(jì)算值為 4.72 m,實(shí)測數(shù)據(jù)是4.95 m, 相對誤差僅為4.5%。計(jì)算最大增水位1.3 m, 與實(shí)測最大增水的相對誤差是2.6%。表 1給出了兩種模式的風(fēng)暴潮計(jì)算值與測量值的數(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)暴潮潮位及增水驗(yàn)證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)暴潮潮位及增水驗(yàn)證Fig. 7 Tidal level verification of astronomic tide,storm surge and storm surge setup verification under equal depth distribution mode

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

4 結(jié)論

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

[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ù)報(bào)技術(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ù)報(bào)模型(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é)報(bào), 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é)報(bào), 2009, 39:287-296.

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

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

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

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

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

[17] 李艷蕓, 李紹武. 風(fēng)暴潮預(yù)報(bào)模式理論及應(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é)基金資助項(xiàng)目(51079095); 國家自然科學(xué)基金創(chuàng)新研究群體科學(xué)基金(51021004)

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

(本文編輯：劉珊珊)