FVCOM模型關(guān)鍵參數(shù)的處理及其在涌潮模擬中的應用

程文龍,潘存鴻*,吳修廣

(1. 浙江省水利河口研究院，浙江 杭州 310020；2. 浙江省河口海岸重點實驗室，浙江 杭州 310016)

FVCOM模型關(guān)鍵參數(shù)的處理及其在涌潮模擬中的應用

程文龍1,2,潘存鴻*1,2,吳修廣1,2

(1. 浙江省水利河口研究院，浙江 杭州 310020；2. 浙江省河口海岸重點實驗室，浙江 杭州 310016)

通過改進海床阻力系數(shù)和設置合適的垂向紊動背景系數(shù)，應用FVCOM模型成功再現(xiàn)了錢塘江河口強涌潮的演進過程。海床阻力系數(shù)采用Manning公式形式，取值隨水深、地形在0.000 2～0.002 9之間變化；垂向紊動背景系數(shù)取1×10-4m2/s。模擬結(jié)果較好地復演了涌潮到達時刻、涌潮高度及涌潮抬升過程、涌潮水平流速以及其沿垂向分布規(guī)律，表明阻力系數(shù)及垂向紊動背景系數(shù)等關(guān)鍵參數(shù)的改進和處理是合理的，可應用于涌潮三維潮流運動特征模擬。

涌潮；FVCOM模型；三維數(shù)值模擬；阻力系數(shù)；垂向紊動系數(shù)

0 引言

潮波傳播到大陸架以后，產(chǎn)生非線性畸變，進入河口后，變形更加劇烈。在一定條件下，會形成水位驟然抬升的漲潮潮波前鋒線，即為涌潮[1]。涌潮是特殊的淺水間斷流運動，流速大、破壞力強、水流特性復雜，對這一問題的數(shù)值模擬一直是計算水動力學的難點之一，具有很高的學術(shù)價值和實際應用價值。

近10 a來，在涌潮一維和二維大尺度數(shù)值模擬方面取得了很大的進展[2-5]，尤其是采用ZHOU et al[6]和HUI et al[7]分別提出的水面梯度法和水位底床法等“和諧”方法后，二維涌潮數(shù)值模擬已相當成熟[8]。但是目前涌潮大尺度三維模擬研究較少，SIMON et al[9]采用LES方法對水槽產(chǎn)生的弱涌潮傳播進行了三維數(shù)值模擬，其尚未應用到實際河口模擬涌潮；王燦星 等[10]采用商用Fluent軟件，通過VOF方法處理自由面，對錢塘江彎道水流特征進行了分析；謝東風 等[11]應用FVCOM模型進行了錢塘江河口實際涌潮的三維模擬，再現(xiàn)了涌潮到達時流速及潮位的突變過程，但水平流速的垂向差異不夠明顯，與實測數(shù)據(jù)存在一定的差異。本文主要通過改進FVCOM模型的阻力系數(shù)，同時選取合適的垂向紊動背景系數(shù)，明顯改善了錢塘江河口涌潮三維數(shù)值模擬結(jié)果。

1 FVCOM模型簡介及其涌潮模型的建立

1.1 FVCOM模型簡介

FVCOM模型原始方程為雷諾平均的三維海洋控制方程組，包含水體連續(xù)性方程及動量方程,溫度、鹽度方程和紊流方程等[12-13]，紊流方程采用修正Mellor-Yamada 2.5階紊流模型計算，其在σ坐標下形式如下：

(1)

(2)

FVCOM原始方程組可選用半隱方法或者模式分裂求解，半隱方法計算時需要調(diào)用PETSC庫。模式分裂法中，外模(正壓)方程是二維的，基于CFL條件和重力外波波速，時間步長較短；內(nèi)模(斜壓)方程是三維的，基于CFL條件和內(nèi)波波速，時間步長較長。另外，F(xiàn)VCOM水平方向采用非結(jié)構(gòu)三角形網(wǎng)格，在垂直方向上可采用σ坐標、s坐標和混合坐標系統(tǒng)，同時FVCOM可以較好地處理漲潮和落潮帶來的干濕邊界轉(zhuǎn)換問題。

1.2 涌潮觀測資料

2010年10月9日至17日在錢塘江鹽官河段開展了一個完整潮汛期的涌潮水文觀測。共布置了3個斷面5個水文測點，用ADCP觀測流速和流向，另外布置6個潮位站，分別位于3個水文斷面的兩岸。具體為：鹽官斷面3個水文測點(左、中、右各1個點)、丁橋(大缺口上游2～3 km)主槽1個水文點和胡斗(老鹽倉下游約3～4 km處)主槽1個水文點。水文點位布置如圖1所示。流速、流向在涌潮到達后15 min內(nèi)每隔1 min記錄1次，15 min至30 min內(nèi)每隔2 min記錄1次，30 min至60 min內(nèi)每隔5 min記錄1次，其它時間每隔30 min記錄1次。同時,涌潮到達后30 min內(nèi)每隔1 min記錄1次潮位，30 min至60 min內(nèi)每隔5 min記錄1次潮位，其余時間每隔15 min記錄1次。

1.3 涌潮模型的建立

計算區(qū)域上邊界取在閘口，下邊界為澉浦。該區(qū)域包含涌潮形成和涌潮最強的河段(圖1)，涌潮特征明顯，同時該區(qū)域已有不少測量和計算成果可供驗證參考。

模型水平方向采用無結(jié)構(gòu)三角形網(wǎng)格對計算域進行剖分，共布置了25 785個單元和13 580個節(jié)點，最大步長532 m，網(wǎng)格最小步長56 m，對鹽官北岸丁壩等水工建筑物密集區(qū)域進行了加密。垂向上采用sigma坐標，共劃分20層，按水深均勻分布。壓力計算采用靜壓模式。紊流模型采用FVCOM自帶的修正Mellor-Yamada 2.5階紊流模型。計算采用內(nèi)、外模分裂模型，因涌潮的強非線性效應，對模型穩(wěn)定性要求較高，外模時間步長很小，僅為0.1 s，內(nèi)模時間步長1.0 s。

水下地形資料采用涌潮觀測前后同期的實測水下地形，上、下游缺乏同期實測地形的區(qū)域采用2010年11月實測水下地形。

圖1 錢塘江河口形勢及站位分布圖Fig.1 Outline of Qiantang Estuary and distribution of stations

2 FVCOM模型關(guān)鍵參數(shù)的改進及其三維涌潮模擬

FVCOM模型本身是一個不斷發(fā)展的、非常優(yōu)秀的大洋模型，當應用于錢塘江河口地區(qū)時，因該區(qū)域水動力過程與大洋相比有很多不同，如：強對流、水淺、河床往往灘槽交互等特征使得其水動力過程復雜多變，特別是涌潮到達時，數(shù)秒內(nèi)，水位驟然上漲2 m左右，高者3 m以上；水流從落潮狀態(tài)急速轉(zhuǎn)化為漲潮狀態(tài)，且垂向?qū)α鹘粨Q劇烈，實測最大垂向流速為0.73 m/s。涌潮過后的數(shù)分鐘至數(shù)十分鐘間，流速達到極值，一般水平流速為6～7 m/s。因此，將FVCOM模型應用到錢塘江河口涌潮模擬時，需對能顯著影響水動力結(jié)果的重要參數(shù)進行敏感性分析以確定合適的取值，必要時還需要對這些參數(shù)的處理方式進行改進。

2.1 海床阻力系數(shù)

FVCOM模型中默認的海床阻力系數(shù)計算公式[12-13]為：

(3)

其中：k為卡門常數(shù)，zab為最底部水層的厚度。z0為海床粗糙高度，其值的大小主要與床面泥沙的粒徑和級配有關(guān)，對于床面由均勻沙組成的情況，z0即可取泥沙粒徑；對于非均勻沙組成的床面，常選一個較粗的代表粒徑(如D90等)作為當量粗糙高度；還可以取更大一些的值，如2.5D90來表征床面起伏后大顆粒突出的影響。一般來說，粗糙高度z0與二維淺水動力學模型常用的Chezy系數(shù)C、Manning系數(shù)n相比，其物理意義明顯，有利于直觀地選擇其量值[14]。但是該式主要在床面平整、無沙波形態(tài)消長、床面阻力以沙粒阻力為主的大洋區(qū)域應用較好，而錢塘江強涌潮區(qū)域則顯著不同，河床灘槽交互，動床阻力主要以沙波阻力為主，而且其阻力系數(shù)普遍小于0.002 5，因此采用默認的阻力系數(shù)計算公式很難在涌潮計算中得到滿意的結(jié)果。

本文采用Manning系數(shù)來處理河(海)床底部阻力[14]，即將阻力系數(shù)一般表達式中的Chezy系數(shù)改用Manning系數(shù)，采用下式計算：

(4)

式中：n為Manning系數(shù)，D為水深，g為重力加速度。

計算中分別采用恒定阻力系數(shù)0.002 5和改進后的阻力系數(shù)(式4)兩種參數(shù)進行對比，計算結(jié)果顯示，盡管模型上、下游邊界均為給定實測潮位過程，阻力系數(shù)為0.002 5時，與實測值相比，鹽官河段沿程低潮位計算值偏高、高潮位計算值偏低，潮差計算值偏小，潮流流速減小也更明顯。以鹽官1號站為例(圖2)，式(3)計算得到的表層流速比式(4)普遍小10%～35%，式(4)結(jié)果與實測值吻合更好。另外，汪亞平 等[15]認為在小型潮汐汊道、潮溝、淺水河口、海灣、潟湖中，水深較小，潮流、波浪等動力活躍，整個水層可視為邊界層，這也說明采用式(4)總水深計算阻力系數(shù)理論上是適宜的。

圖2 不同阻力系數(shù)下鹽官1號站表層水平流速過程對比圖Fig.2 The comparison of surface horizontal velocity with different roughness coefficient at Yanguan 1 station

圖3 不同垂向背景紊流系數(shù)下鹽官1號站漲、落急水平流速垂向分布對比圖Fig.3 The comparison of vertical distribution of the maximum flood and ebb velocity with different background vertical eddy viscosity at Yanguan 1 station

2.2 垂向紊動系數(shù)

垂向紊動系數(shù)為二階項中的系數(shù)，其對水平流速的大小和垂向分布影響較大[16]，二階項對流速的垂向分布起均勻化作用，垂向紊動系數(shù)越大，流速垂向分布越均勻，且流速越小。FVCOM湍流模型采用修正Mellor-Yamada 2.5階紊流模型，其垂向紊流粘性系數(shù)Km通過湍流模型計算后再加上垂向背景紊流系數(shù)UMOL得到，即Km?Km+UMOL，在本研究河段Km值一般為10-4～10-3m2/s。經(jīng)數(shù)值試驗，UMOL取較小值(10-5m2/s左右)容易導致計算不穩(wěn)定，取較大值(10-3m2/s左右)水平流速垂向分布會被勻化，且表層水平流速也減少5%～10%(圖3)。因此，涌潮計算中UMOL一般取10-4量級，本文取1×10-4m2/s。

2.3 三維涌潮計算分析

2.3.1 大范圍潮位及潮流模擬結(jié)果

FVCOM模型通過式(4)計算海床阻力系數(shù)，背景垂向紊動系數(shù)取1×10-4m2/s，對鹽官河段大潮期涌潮開展驗證。其中Manning系數(shù)取值隨河床可動性、灘槽的不同一般為0.006～0.020，與二維淺水動力學模型取值接近[5]，較為簡單方便。對應阻力系數(shù)Cd值在0.000 2～0.002 9之間變化，其上限比下限大一個數(shù)量級，接近式(3)上限，且可以通過調(diào)整每個單元的Manning系數(shù)值來綜合反映河床深槽和邊灘的影響。圖4為大潮期鹽官河段潮位驗證圖，可見高低潮位值、潮差和相位計算與實測均吻合較好。潮位驗證也反映了該河段大潮期漲潮時潮位變化劇烈，涌潮高度自下而上逐漸增大的趨勢。

將水平流速模擬結(jié)果與實測值按照平均水深、水面以下0.5 m、水面以下1.0 m等至河底以每層層寬0.5 m為間隔繪制成流速、流向過程(圖5)，可見流速過程驗證吻合較好，特別是結(jié)果中各站點垂向上漲急流速、落急流速均能捕捉到，流向過程驗證值與實測吻合良好。在空間上，自下而上，鹽官河段漲潮最大流速從丁橋至鹽官達到最大，到胡斗又開始減小。

2.3.2 涌潮高度及抬升過程模擬結(jié)果

將涌潮到達前后，鹽官北站潮位過程的模擬結(jié)果以每分鐘1次的頻率輸出，與同頻率的實測過程進行比較，如圖6所示。當涌潮到達時，水位在1 min內(nèi)抬升接近3 m，再經(jīng)過幾次小幅抬升后達到高潮位，繼而開始衰減，潮位逐漸降低。這一潮位過程在模型中被復演，模擬結(jié)果在涌潮起漲時間、涌潮高度、小幅抬升以及衰減過程上都較好地吻合了實測值，說明FVCOM模型對流項不作梯度限制處理也能很好地捕捉到涌潮間斷過程，僅在涌潮到達后初次抬升幅度(涌潮高度)模擬結(jié)果比實測值小約0.8 m，這可能與靜壓模式中潮頭動能轉(zhuǎn)化為勢能不充分有關(guān)。

圖4 潮位驗證圖Fig.4 Validation diagram of tidal level

圖5 鹽官1號站分層流速過程驗證圖Fig.5 Validation diagram of current in different depth at Yanguan 1 station

圖6 涌潮到達前后鹽官北潮位過程Fig.6 The tidal level at Yanguanbei before and after tidal bore arrival

圖7 涌潮到達前后鹽官1號站表層、中層及底層流速、流向過程Fig.7 The surface, middle and bottom velocity magnitude & direction at Yanguan 1 station before and after tidal bore arrival

2.3.3 涌潮過后不同深度流速過程模擬結(jié)果

涌潮到達前后鹽官1號水文站表層、中層以及底層附近水平流速過程與實測值對比如圖7，當涌潮到達時，表層水平流速從落潮1.41 m/s立刻轉(zhuǎn)變?yōu)闈q潮2.65 m/s，很快漲潮流速回落到2.0 m/s以內(nèi)，之后流速繼續(xù)增大，反復一、二次后漲潮流速達到最大，此時離初漲時刻約20 min，模擬潮流結(jié)果仍可以復演該過程，且從圖中可見，表、中、底層流速均與實測值趨勢一致。

總體而言，F(xiàn)VCOM可以較好地模擬涌潮起漲、發(fā)展以及消失等整個階段的潮位和水平方向潮流過程，本次模擬結(jié)果較好地反映了本區(qū)域涌潮傳播的過程，可用于涌潮三維潮流運動特征分析。

3 小結(jié)

采用國際、國內(nèi)廣泛應用的FVCOM水動力模型，通過改進海床阻力系數(shù)和設置合適的垂向紊動背景系數(shù)，在靜壓假定下，較好地復演了錢塘江河口強涌潮的演進過程，對涌潮到達時刻、涌潮高度、涌潮抬升過程、涌潮水平流速以及其沿垂向分布規(guī)律等的模擬結(jié)果與實測吻合較好。不過靜壓假定忽略了水深方向的流動影響，弱化了涌潮潮頭“水滾”模擬效果，因此下一步工作擬加強對非靜壓假定的涌潮的全過程三維數(shù)值模擬。

致謝 感謝陳長勝教授領(lǐng)導的SMAST/UMASSD及MEDM研究團隊開發(fā)、提供FVCOM源代碼。

[1] PAN Cun-hong. Numerical simulation for discontinuous shallow water flow and its application to the analysis of the tidal bore at the Qiantang Estuary [D]. Shanghai: Shanghai University,2007.

潘存鴻.淺水間斷流動數(shù)值模擬及其在錢塘江河口涌潮分析中的應用[D].上海:上海大學,2007.

[2] DU Shan-shan, XUE Lei-ping. Research on one-dimension numerical simulation in tide flow computation[J]. Shanghai Water,2006，22(2):44-47,29.

杜珊珊,薛雷平.長江口北支涌潮的一維數(shù)值模擬[J].上海水務，2006，22(2):44-47,29.

[3] YU Pu-bing, PAN Cun-hong. 1D numerical simulation of tidal bore in Qiantang River[J].Chinese Journal of Hydrodynamics,2010，25(5)：669-675.

于普兵，潘存鴻.錢塘江涌潮一維數(shù)值模擬[J].水動力學研究與進展:A輯，2010，25(5)：669-675.

[4] PAN C H, LIN B Y, MAO X Z. Case study: Numerical modeling of the tidal bore on the Qiantang River, China [J]. Journal of Hydraulic Engineering,2007,133(2):130-138.

[5] PAN Cun-hong, LU Hai-yan, ZENG Jian. Characteristic and numerical simulation of tidal bore in Qiantang River[J]. Hydro-Science and Engineering,2008(2):1-9.

潘存鴻,魯海燕,曾劍.錢塘江涌潮特性及其數(shù)值模擬[J].水利水運工程學報,2008(2):1-9.

[6] ZHOU J G, CAUSON D M, MINGHAM C G, et al. The surface gradient method for the treatment of source terms in the shallow water equations[J]. Journal of Computational Physics,2001,168(1):1-25.

[7] HUI W H, PAN C H. Water level-bottom topography formulation for the shallow-water flow with application to the tidal bores on the Qiantang river[J]. Computational Fluid Dynamics Journal,2003,12(3):549-554.

[8] PAN Cun-hong. Advances in numerical simulation of discontinuous shallow water flows[J]. Advances in Science and Technology of Water Resources,2010,30(05):77-84.

潘存鴻.淺水間斷流動數(shù)值模擬研究進展[J].水利水電科技進展,2010,30(05):77-84.

[9] SIMON B, LUBIN P, GLOCKNER S, et al. Three-dimensional numerical simulation of the hydrodynamics generated by a weak breaking tidal bore[C]//Proc. 34th IAHR World Congress, Brisbane, Australia,26 June-1 July,2011:1 133-1 140.

[10] WANG Can-xing, CHEN Ju-fang, JIN Han-hui, et al. Three-dimensional numerical simulation of tidal bore of Qiantang River[J]. Chinese Journal of Hydrodynamics,2012,27(4):367-375.

王燦星，陳菊芳，金晗輝,等.涌潮對錢塘江河道流場影響的三維數(shù)值模擬研究[J].水動力學研究與進展:A輯,2012,27(4):367-375.

[11] XIE Dong-feng, PAN Cun-hong, WU Xiu-guang. Three-dimensional numerical modeling of tidal bore in Qiantang based on FVCOM[J]. The Ocean Engineering,2011，29(1)：47-52.

謝東風，潘存鴻，吳修廣.基于FVCOM模式的錢塘江河口涌潮三維數(shù)值模擬研究[J].海洋工程，2011，29(1)：47-52.

[12] CHEN C S, LIU H D. An unstructured grid, finite-volume, three-dimensional, primitive equations ocean model: Application to coastal ocean and estuaries[J]. Journal of Atmospheric and Oceanic Technology,2003,20(01):159-186.

[13] CHEN C S, BEARDSLEY R C, COWLESET G, et al. An unstructured grid, finite-volume community ocean model FVCOM user manual [M].2013.

[14] SHAO Xue-jun, WANG Xing-kui. Introduction to river mechanics[M]. Beijing: Tsinghua University Press,2005:61-62.

邵學軍,王興奎.河流動力學概論[M].北京：清華大學出版社有限公司,2005:61-62.

[15] WANG Ya-ping, GAO Shu, JIA Jian-jun. Flow structure in the marine boundary layer and Bedload Transport: A Review[J]. Marine Geology & Quaternary Geology,2000,20(3):101-106.

汪亞平,高抒，賈建軍.海底邊界層水流結(jié)構(gòu)及底移質(zhì)搬運研究進展[J].海洋地質(zhì)與第四紀地質(zhì),2000,20(3):101-106.

[16] CHEN Yong-ping, LIU Jia-ju, YU Guo-hua. A Study on eddy viscosity coefficient in numerical tidal simulation[J]. Journal of Hehai University: Natural Sciences,2002,30(01):39-43.

陳永平,劉家駒,喻國華.潮流數(shù)值模擬中紊動粘性系數(shù)的研究[J].河海大學學報:自然科學版,2002,30(01):39-43.

Processing on key parameter in FVCOM and its application on tidal bore simulation

CHENG Wen-long1,2, PAN Cun-hong*1,2, WU Xiu-guang1,2

(1.ZhejiangInstituteofHydraulicsandEstuary,Hangzhou310020,China; 2.KeyLaboratoryofEstuaryandCoastofZhejiangProvince,Hangzhou310016,China)

The evolution of tidal bore in Qiantang Estuary were successfully reproduced by applying FVCOM model through improving the sea bed roughness coefficient and setting suitable background vertical eddy viscosity. Sea bed roughness coefficient，calculated using the Manning formulation, is between 0.000 2～0.002 9 with the varying depth. The background vertical eddy viscosity is 1×10-4m2/s. The simulation well reproduces the tidal bore arriving time, tidal level, setup height, horizontal velocity magnitude and vertical profile, which implies the choice of parameters such as the bed roughness coefficient and background vertical eddy viscosity are reasonable. The procedure could be extended to simulate the 3D tidal movement for tidal bores.

tidal bore; FVCOM; 3D mathematic simulation; roughness coefficient; vertical eddy viscosity

10.3969/j.issn.1001-909X.2017.01.004.

2016-02-25

2017-01-05

國家自然科學基金項目資助(51379190)；水利部公益性行業(yè)專項項目資助(201001072)；浙江省省屬科研院所專項計劃項目資助(2016F50018);浙江省科技廳創(chuàng)新團隊與人才培養(yǎng)項目資助(2012F20031)；浙江省公益技術(shù)應用研究計劃項目資助(2016C33095)

程文龍(1981-)，男，湖北公安縣人，高級工程師，主要從事河口海岸動力學研究。E-mail:chengwl@zjwater.gov.cn

*通訊作者：潘存鴻(1963-)，男，教授級高級工程師，主要從事河口海岸水動力學、泥沙及水環(huán)境研究。E-mail：panch@zjwater.gov.cn

TV131.2

A

1001-909X(2017)01-0033-08

10.3969/j.issn.1001-909X.2017.01.004

程文龍，潘存鴻，吳修廣.FVCOM模型關(guān)鍵參數(shù)的處理及其在涌潮模擬中的應用[J].海洋學研究,2017,35(1):33-40,

CHENG Wen-long, PAN Cun-hong, WU Xiu-guang. Processing on key parameter in FVCOM and its application on tidal bore simulation[J].Journal of Marine Sciences,2017,35(1):33-40, doi:10.3969/j.issn.1001-909X.2017.01.004.