湍黏系數(shù)對淺灘海域三維風暴潮的影響

熊 偉,朱志夏,董 佳,齊慶輝

(1.江蘇省交通規(guī)劃設計院股份有限公司,江蘇 南京 210014;2.江蘇省水運工程技術研究中心,江蘇 南京 210014)

湍黏系數(shù)對淺灘海域三維風暴潮的影響

熊 偉1,2,朱志夏1,2,董 佳1,2,齊慶輝1,2

(1.江蘇省交通規(guī)劃設計院股份有限公司,江蘇 南京 210014;2.江蘇省水運工程技術研究中心,江蘇 南京 210014)

利用WRF大氣模式和基于有限體積法的三維海洋數(shù)值模型FVCOM,建立了連云港及其附近淺灘海域的三維風暴潮流數(shù)學模型。運用k-ε和MY-2.5兩種不同的湍流模型,對“韋帕”臺風作用下的三維風暴潮流進行了數(shù)值模擬,并對計算得到的三維風暴潮流流場結(jié)果和垂向湍流黏性系數(shù)結(jié)果進行了對比分析。計算分析結(jié)果表明,兩種不同湍流模型計算所得的風暴潮增水差別很小,可以忽略;但在風生流較強的時間段內(nèi),MY-2.5湍流模型計算所得的水平流速流向沿垂向分布并不一致,而k-ε湍流模型的計算結(jié)果較好,與實測分層潮流資料更為符合。研究結(jié)果表明,垂向湍流黏性系數(shù)沿水深的垂向分布對淺灘海域風暴潮流水平流速的垂向結(jié)構至關重要,建議選用k-ε湍流模型計算垂向湍流黏性系數(shù)和湍流擴散系數(shù)。

湍黏系數(shù);淺灘海域;三維數(shù)值模型;風暴潮;垂向分布

風暴潮是指強烈的大氣擾動(強風和氣壓驟變)引起的海面異常升高或降低的現(xiàn)象[1]。風暴潮位若遇到天文大潮,就會發(fā)生大的潮災。而伴隨強風產(chǎn)生的波浪和水流流速很大,使得淺灘海岸泥沙的啟動變得異常容易,其運動也變得較為復雜,極易形成港口航道工程的淤積,如具有淺灘淤泥質(zhì)海岸性質(zhì)的連云港海域。在大風浪作用下,水流和泥沙具有明顯的分層現(xiàn)象,因此加強風暴潮特別是三維風暴潮的數(shù)值模擬研究也顯得尤為必要。

20世紀末,風暴潮數(shù)值模擬研究大多采用二維數(shù)值模型計算增減水[2-4]。隨著計算機的發(fā)展,近年來三維風暴潮數(shù)值模式得到了較大發(fā)展,出現(xiàn)了許多海洋數(shù)值模型如POM,ECOM和FVCOM等[5]。從物理角度看,三維模型可以用來描述風暴潮的水流垂向結(jié)構,并且三維模型在近底處的底部摩阻假設更為真實[6]。強風天氣下的水流運動、泥沙濃度分層現(xiàn)象較為明顯。因此,對風暴潮作用下的三維水流進行準確的模擬,是研究港口航道工程驟淤問題的前提條件。李大鳴等[7]對渤海灣的風暴潮進行了三維數(shù)值模擬,曹叢華等[8]利用FVCOM模型在膠州灣進行了三維風暴潮漫灘數(shù)值模擬研究,但對三維風暴潮研究都只進行了水位的驗證,缺乏三維流速流向的驗證分析;張娜等[9]利用SWAN和MOHID聯(lián)合模型對連云港深水航道進行了風暴潮驟淤研究,取得了較好的成果,但未對三維水流結(jié)構進行詳細分析研究。陳永平等[10]討論了垂向紊動黏性系數(shù)對三維潮流數(shù)值計算結(jié)果的影響,結(jié)果表明垂向紊動黏性系數(shù)越大,流速分層越不明顯,流速的絕對值也越小,垂向紊動黏性系數(shù)的梯度變化則對水平流速的垂向分布起著決定性的作用。本文利用基于有限體積法的FVCOM模型和已有的三維潮流資料,對連云港海域在“韋帕”臺風作用下的三維風暴潮進行數(shù)值模擬,著重研究風暴潮作用下的水流垂向結(jié)構和不同湍流模型的應用對比情況。

1 數(shù)學模型的建立

“韋帕”(WIPHA)臺風于2007年9月16日在西北太平洋洋面上生成,在浙江省蒼南縣霞關鎮(zhèn)登陸,登陸時臺風中心附近最大風速達45 m/s,臺風中心向北移動,過境江蘇省連云港市,對該海域造成一定的影響。有關研究單位在連云港徐圩海域?qū)Υ舜闻_風進行了連續(xù)現(xiàn)場觀測,獲得了較為詳細的波浪潮流泥沙等觀測資料。本文建立臺風模型和三維潮流模型,對該時段的三維風暴潮流進行數(shù)值模擬研究。

1.1 大氣模式

WRF(Weather Research Forecast)大氣模式為新一代中尺度預報模式與同化系統(tǒng),具有可移植、易維護、可擴充、高效率等特性,可用于數(shù)值天氣預報的研究、大氣-海洋模式的耦合等。目前,WRF模式包括2個動力框架ARW(the Advanced Research WRF)和NMM(the Nonhydrostatic Mesoscale Model),分別由美國國家大氣研究中心(NCAR)和美國國家環(huán)境預報中心(NCEP)主要開發(fā)并維護更新。本文使用的是WRF-ARW動力框架,WRF-ARW系統(tǒng)主要由4部分組成:WRF預處理系統(tǒng)(WPS)、WRF數(shù)據(jù)同化系統(tǒng)(WRFDA)、ARW求解、程序后處理以及可視化工具[11]。

WRF模式是一個完全可壓的非靜力模式,控制方程都寫為通量的形式。垂直方向采用地形跟隨的質(zhì)量坐標,水平方向采用Arakawa C網(wǎng)格[12]。

1.2 風暴潮數(shù)值模型

FVCOM采用非結(jié)構化的三角形網(wǎng)格單元,能更精確擬合復雜不規(guī)則的岸界,如處理連云港東西連島海域岸線;采用有限體積法,能夠在整個計算區(qū)域內(nèi)更好地保證各物理量的守恒;采用垂向σ坐標,可以體現(xiàn)不規(guī)則變化的海底地形,并在淺水區(qū)域具有更高的垂向分辨率。另外,該數(shù)值模型程序代碼的并行化計算,大大縮短了數(shù)值計算求解的時間。其具體控制方程和數(shù)值離散求解模式詳見文獻[5]。

適當?shù)孛枋龊陀嬎阃牧?對于三維水流的模擬、預測泥沙輸運和海床變化具有重要的影響, L.O.Amoudry等人詳細研究了湍流模型對泥沙分層和輸運的影響[13]。本文在研究過程中發(fā)現(xiàn),選用不同的湍流模型,得到的三維水流計算結(jié)果差異較大。在此,對湍流模型作簡要介紹。FVCOM海洋數(shù)值模型提供了多樣性的湍流封閉模型,用于計算垂向渦黏性系數(shù)和垂向擴散系數(shù)。FVCOM自帶了改進過的MY-2.5湍流封閉模型,此外,它通過耦合GOTM湍流模式,也包含了k-ε湍流封閉模型供研究人員選擇。

1.2.1 MY-2.5湍流模型 目前比較常用的MY-2.5湍流模型,它包括湍動能方程和混合長度方程,可表示如下:

式中:~W=1+E2l2/(κL)2為固壁迫近函數(shù),κ=0.4為卡曼常數(shù);E1,E3和B1為模型閉合常數(shù);Km,Kh和Kq分別為垂向渦黏系數(shù)、垂向擴散系數(shù)和垂向湍流擴散系數(shù);q2/2為湍動能;l為湍流混合長度。MY-2.5模型中,表面及底面的邊界條件分別為:

式中:η和H分別為潮位和水深。

1.2.2 k-ε湍流模型 GOTM湍流模式實現(xiàn)了眾多的從簡單的Richardson參數(shù)化到復雜的雷諾應力等湍流模型,其中k-ε湍流模型可表示為:

式中:Km為垂向渦黏系數(shù);k為湍流動能;ε為湍流耗散;經(jīng)驗常數(shù)cμ,c1,c2,σk,σε的取值分別為0.09,1.44, 1.92,1.00,1.30。表面及底面邊界條件分別為:

2 數(shù)學模型的驗證

2.1 韋帕臺風驗證

為了給潮流計算提供更高精度的風場數(shù)據(jù),減小邊界對計算結(jié)果的影響,并模擬出臺風條件下的臺風路徑。本次研究中,“韋帕”臺風的風場模型計算區(qū)域采用大模型、中模型與小模型三重網(wǎng)格雙向嵌套技術,即大模型的計算結(jié)果直接影響小模型,且小模型的計算結(jié)果也會對大模型進行反饋。模型的計算范圍與網(wǎng)格精度如表1所示,模型垂直方向分為34層。模擬的臺風路徑與實際路徑對比見圖1,風速、風向與連云港大西山海洋站實測數(shù)據(jù)對比見圖2。由圖2可見,建立的WRF大氣模式能較好地反映連云港海域風場狀況。

圖1 模擬臺風路徑與實際路徑對比Fig.1 Comparison between simulated typhoon path and actual path

表1 WRF模型計算范圍與網(wǎng)格精度Tab.1 Range and accuracy of grid in WRF model

圖2 臺風風速風向?qū)Ρ闰炞CFig.2 Validation of typhoon speed and direction

2.2 三維風暴潮驗證

本次風暴潮數(shù)值模擬研究采用二維大模型、三維小模型雙層嵌套的方法計算。二維大模型主要為三維小模型的計算提供合理的邊界條件。二維、三維模型的計算范圍如圖3所示。模型均采用非結(jié)構化的三角形網(wǎng)格,網(wǎng)格最小尺度為30 m,位于連云港主港區(qū)和主航道;網(wǎng)格最大尺度為2 000 m,位于計算域的離岸開邊界處。測站相對位置如圖4所示,1#高程為-3.0 m,2#高程為-5.0 m(高程基準為當?shù)乩碚摶?,但1#點沒有流速驗證資料,因此采用2#的實測潮位流速資料進行驗證。

圖3 雙重嵌套模型計算范圍Fig.3 Calculated range of nested models

圖4 計算驗證點位置Fig.4 Position of measured points

計算時間為2007-09-17T0:00/2007-09-21T20:00,整個過程共計117 h。三維模型的計算時間步長內(nèi)模取1.0 s,外模為內(nèi)模的5倍;垂向均勻分為10層。由于連云港海域基本為淤泥質(zhì)海岸性質(zhì),摩擦阻力系數(shù)相對較小,底部粗糙厚度取0.5 mm。風場數(shù)據(jù)源于WRF模型計算得到的結(jié)果。需要說明的是,連云港海域水深較淺,水流垂向結(jié)構較為復雜。因此,本次數(shù)值模擬研究中湍流模型選MY-2.5和k-ε兩種模型分別計算,計算驗證結(jié)果如圖5和6所示。

圖5 2#觀測點潮位驗證Fig.5 Verification of tidal level for 2#measured point

圖6 2#觀測點流速流向驗證Fig.6 Verification of velocity and direction of current for 2#measured point

從潮位計算結(jié)果可知,兩種湍流模型都能夠很好地模擬風暴潮位的過程,潮位過程差別不大,這說明湍流黏性系數(shù)數(shù)值變化對風暴潮位的影響微弱;從分層水流計算結(jié)果可知,相對于MY-2.5模型來說,k-ε模型的計算結(jié)果和實測資料更相符。很顯然,在其他基本參數(shù)一致的情況下,三維風暴潮流結(jié)果差異顯著的原因是采用的湍流模型不同造成的。下面通過MY-2.5模型和k-ε模型的計算,對流速過程和垂向渦黏性系數(shù)等進行對比分析。

3 三維水流結(jié)構與湍黏系數(shù)

選取2個特征性的代表時刻進行對比分析。風前時刻為2007-09-18T06:00,模擬風速為0.94 m/s;大風時刻為09-20T08:00,模擬風速為16.31 m/s。2個時刻的表層與底層流場分別如圖7和8所示?？梢钥闯?風前時刻流場保持著正常天氣下的潮流漲落特性,但是在風中時刻,受大風影響,k-ε模型此刻近岸水體表底層流向基本一致,全變?yōu)檠匕读鞯男再|(zhì),這與實測資料基本吻合。這說明在淺灘區(qū)域,淺灘海域的水體在垂向上作為一個整體對臺風風暴作出響應。而大風時刻MY-2.5模型的計算結(jié)果中,表底層流向差異較大,表現(xiàn)為表層是沿岸流,底層近乎為離岸流,這與實測情況不符。

圖7 大風時刻MY-2.5模型計算的表、底層流場Fig.7 Surface and bottom flow field calculated by MY-2.5 model at strong wind moment

圖8 大風時刻k-ε模型計算的表、底層流場Fig.8 Surface and bottom flow field calculated by k-ε model during strong wind

計算輸出2個測點在風前、風中兩個不同時刻的垂向渦黏性系數(shù)Km和垂向擴散系數(shù)Kh(見圖9)。由于Km與Kh具有相關性,僅對Km作詳細比較,Kh類同。

從不同水深點的Km值比較來看,風前時刻與大風時刻,-5.0與-3.0 m的Km值垂向分布趨勢基本一致,不同的是-5.0 m的湍流黏性系數(shù)值比-3.0 m稍大,這與兩者水流流速大小有關。另外,k-ε模型計算得到表底面處的Km為0,這與模型的表底面邊界條件式(8)和(9)有關。湍流模型中沒有考慮波浪的影響,因此表底面的Km計算值與實際情況存在差別。

從不同時刻來看,如圖9(a)和(b)所示,風前時刻的兩種模型計算的湍流黏性系數(shù)分布具有些許相似,都呈現(xiàn)中間層大、表底層小的特征,但是Km的最大值分布位置不同,MY-2.5模型偏底層,而k-ε模型偏表層。圖9(c)和(d)則說明兩種模型計算結(jié)果差異很大,MY-2.5模型的Km表現(xiàn)出了間斷性不連續(xù)現(xiàn)象,且Km最大值分布在表層,數(shù)值達1.0左右;而k-ε模型的結(jié)果依舊與風前時刻相似,但Km的最大值在垂向上的位置分布向下移,更趨于水深中間層,數(shù)值也大了許多,-3.0和-5.0 m處的Km的最大值分別為0.014和0.019 m2/s。需要強調(diào)的是,MY-2.5模型的Km表現(xiàn)出了間斷性不連續(xù)現(xiàn)象,造成了三維風暴潮流在垂向流速分布上不一致,且間斷層與水流垂向流速流向的突變層均位于0.6H層。為此,本文嘗試在MY-2.5模型計算中將垂向分層層數(shù)加大為20,但Km值依舊在0.6H層出現(xiàn)間斷性分布。這說明,在淺灘海域的三維風暴潮流數(shù)值模擬計算中,垂向湍流黏性系數(shù)這一參數(shù)的計算至關重要,它關系到風暴潮流水平流速的垂向結(jié)構,與文獻[10]中的結(jié)論一致。另外,由于Km與Kh具有相關性,而垂向擴散系數(shù)Kh值也關系到風暴潮流作用下泥沙、污染物等物質(zhì)的輸運擴散;因此應選用適當?shù)耐牧髂Ｐ陀嬎愦瓜驕u黏性系數(shù)和擴散系數(shù)。

圖9 風前時刻和大風時刻的測站垂向湍流黏性系數(shù)和擴散系數(shù)分布Fig.9 Vertical distribution of turbulent viscosity coefficients and diffusion coefficients before and during strong wind

4 結(jié) 語

本文對“韋帕”臺風作用下連云港淺灘海域的三維風暴潮流進行數(shù)值模擬和研究,采用k-ε模型得到了比較好的計算結(jié)果,與實測風暴潮流資料符合較好,并得出以下結(jié)論:

(1)在淺灘海域的三維風暴潮流數(shù)值模擬計算中,垂向湍流黏性系數(shù)這一參數(shù)的計算和取值至關重要,它關系到風暴潮流的垂向結(jié)構;應選用適當?shù)耐牧髂Ｐ?建議采用k-ε湍流模型進行計算。

(2)在大風作用達到一定時間之后,淺灘海域的水體在垂向上作為一個整體對臺風風暴作出響應。

[1]馮士筰.風暴潮導論[M].北京:科學出版社,1982:1-2.(FENG Shi-zuo.Storm surge introduction[M].Beijing:Science Press,1982:1-2.(in Chinese))

[2]李樹華,陳文廣,梁善任,等.廣西沿海8007號臺風暴潮數(shù)值模擬及臺風暴潮某些特性的分析[J].廣西科學,1994,1 (2):14-20.(LI Shu-hua,CHEN Weng-guang,LIANG Shan-ren.Numerical simulation and certain characteristics analysis on the storm surge No.8007 along the coast of Guangxi[J].Guangxi Sciences,1994,1(2):14-20.(in Chinese))

[3]余曉軍,呂開華.9618號臺風風暴潮特征及其數(shù)值模擬[J].海洋通報,1998,17(6):7-11.(YU Xiao-jun,LV Kai-hua. Characteristics and numerical simulation of the storm surge caused by typhone No.9618[J].Marine Science Bulletin,1998,17 (6):7-11.(in Chinese))

[4]柴扉,汪景庸.東海風暴潮與天文潮的非線性相互作用[J].青島海洋大學學報,1990,20(3):56-62.(CAI Fei,WANG Jing-yong.The nonlinear interaction of surge and tide in the East China Sea[J].Journal of Ocean University of Qingdao,1990, 20(3):56-62.(in Chinese))

[5]熊偉.非結(jié)構網(wǎng)格下的三維潮流泥沙數(shù)值模擬[D].大連:大連理工大學,2012:3-4.(XIONG Wei.3D numericalsimulation of tidal currents and sediment based on unstructured mesh model[D].Dalian:Dalian University of Technology,2012: 3-4.(in Chinese))

[6]黃華,朱建榮,吳輝.長江口與杭州灣風暴潮三維數(shù)值模擬[J].華東師范大學學報:自然科學版,2007(4):9-19. (HUANG Hua,ZHU Jian-rong,WU Hui.3-D numerical simulation of storm surge in the Changjiang estuary and the Hangzhou bay[J].Journal of East China Normal University(Natural Science),2007(4):9-19.(in Chinese))

[7]李大鳴,范玉海,徐亞男,等.風暴潮三維數(shù)值計算模式的研究及在渤海灣的應用[J].海洋科學,2012,36(7):7-13. (LI Da-ming,FAN Yu-hai,XU Ya-nan.Research on 3-D numerical calculation model of storm surge and application in the Bohai Bay[J].Marine Sciences,2012,36(7):7-13.(in Chinese))

[8]曹叢華,白濤,高松,等.膠州灣高分辨率三維風暴潮漫灘數(shù)值模擬[J].海洋科學,2013,37(2):118-125.(CAO Cong-hua,BAI Tao,GAO Song,et al.High resolution 3D storm surge and inundation numerical model used in the Jiaozhou bay [J].Marine Sciences,2013,37(2):118-125.(in Chinese))

[9]張娜,楊華,嚴冰,等.SWAN和MOHID聯(lián)合模型在連云港港30萬噸級航道建設中的應用[J].水運工程,2012(2): 128-133.(ZHANG Na,YANG Hua,YAN Bing,et al.Application of SWN/MOHID-based three-dimensional numerical model in 300 000 DWT deep-water channel of Lianyungang port[J].Port&Waterway Engineering,2012(2):128-133.(in Chinese))

[10]陳永平,劉家駒,喻國華.潮流數(shù)值模擬中紊動粘性系數(shù)的研究[J].河海大學學報:自然科學版,2002,30(1):39-43.(CHEN Yong-ping,LIU Jia-ju,YU Guo-hua.A study of eddy viscosity coefficient in numerical tidal simulation[J]. Journal of Hohai University(Natural Sciences),2002,30(1):39-43.(in Chinese))

[11]WANG W,BRUYèRE C,DUDA M,et al.ARW Version 3 modeling system user′s guide:mesoscale and microscale meteorology division[M].Boulder:National Center for Atmosphere Research,2011.

[12]王勇.基于WRF模式的酒泉風電基地近地層風場特征模擬研究[D].蘭州:蘭州大學,2012:25-26.(WANG Yong. Simulation and analysis of wind characteristics near surface layer in Jiuquan with WRF[D].Lanzhou:Lanzhou University, 2012:25-26.(in Chinese))

[13]AMOUDRY L O,SOUZA A J.Impact of sediment-induced stratification and turbulence closures on sediment transport and morphological modeling[J].Continental Shelf Research,2011,31(9):912-928.

Effects of turbulent viscosity coefficient on 3-D storm surge within shallow seas

XIONG Wei1,2,ZHU Zhi-xia1,2,DONG Jia1,2,QI Qing-hui1,2
(1.Jiangsu Province Communications Planning and Design Institute Limited Company,Nanjing 210014,China; 2.Jiangsu Province Waterway Engineering Technology Research Center,Nanjing 210014,China)

Using a WRF atmospheric model and a 3-D ocean numerical model FVCOM based on finite volume method,a three-dimensional numerical model of storm surge within Lianyungang shallow waters is developed.The 3-D storm surge during typhoon WIPHA is simulated by use of different turbulence models such as k-ε and MY-2.5.The simulation results of the 3-D storm surge flow and vertical turbulence viscosity coefficient are analyzed in the study.The analysis results show that,compared with the poor results from MY-2.5 model,the vertical velocity and direction distribution of the storm surge calculated by k-ε model has a good agreement with the measured data,and the values of the vertical turbulence viscosity coefficient are crucial for the vertical structure of the storm surge,therefore a right turbulence model is preferred.The k-ε turbulence model is suggested for simulation of the 3-D storm surge in the future.

viscosity coefficients;shallow seas;3D numerical simulation;storm surge;vertical distribution

P731.23

A

1009-640X(2014)06-0045-07

2014-04-26

國家高新技術研究發(fā)展計劃(863計劃)資助項目(2012AA112509)

熊 偉(1988-),男,湖北荊州人,碩士,主要從事海岸近海水動力與泥沙研究。E-mail:xw_7816@126.com通信作者:朱志夏(E-mail:zhixiazhu@sina.com)