寧波舟山近海三維潮汐潮流數(shù)值模擬

熊 偉 ，劉必勁 ，，孫昭晨 ，梁書秀 ，張亦飛

（1.大連理工大學(xué)海岸和近海工程國家重點實驗室，大連 116023；2.國家海洋局杭州海洋工程勘測設(shè)計研究中心，杭州 310012）

寧波舟山近海三維潮汐潮流數(shù)值模擬

熊 偉1，劉必勁1，2，孫昭晨1，梁書秀1，張亦飛2

（1.大連理工大學(xué)海岸和近海工程國家重點實驗室，大連 116023；2.國家海洋局杭州海洋工程勘測設(shè)計研究中心，杭州 310012）

利用FVCOM模型對寧波舟山近海的潮汐潮流進行了三維數(shù)值模擬，并對其水動力特性作了相應(yīng)分析。FVCOM模型采用非結(jié)構(gòu)化三角形網(wǎng)格，很好地解決了精確擬合寧波舟山群島的復(fù)雜岸線的問題。基于計算區(qū)域內(nèi)的海床性質(zhì)，采用Koutitas公式對FVCOM模型的中海底摩阻系數(shù)的計算進行了改進。通過計算域內(nèi)多個潮位站和海流站的實測資料驗證表明，改進的摩阻系數(shù)計算公式是合理的，流場的計算結(jié)果與實測符合良好，可以用于三維污染物擴散和泥沙輸移計算。

摩阻系數(shù)；潮汐；潮流；數(shù)值模擬；FVCOM模型

Biography：XIONG Wei（1988-），male，master student.

寧波舟山及其近海水域島嶼密布，漕灘相間，深水航道眾多，地形極其復(fù)雜。由北往南的海灣主要有杭州灣、象山港和三門灣，這些都是強潮海灣且灣內(nèi)有大面積的淺灘。寧波舟山群島內(nèi)的主要深水航道有金塘水道、冊子水道、螺頭水道、佛渡水道等。

國內(nèi)學(xué)者對于杭州灣及浙江近海的潮流特性和數(shù)值模擬做了一定的研究。曹德明等［1-2］用有限差分法對杭州灣的潮汐、潮流進行了二維數(shù)值模擬；李身鐸［3］采用垂向σ坐標(biāo)下的三維數(shù)值模式模擬了杭州灣三維潮波運動；陳倩等［4］采用改進后的HAMSOM模式對浙江近海潮汐潮流進行了三維數(shù)值模擬；壽瑋瑋［5］采用ECOM模型對舟山群島附近海域水動力場進行了三維數(shù)值模擬。但是大多數(shù)采用矩形網(wǎng)格和調(diào)和常數(shù)來模擬該區(qū)域內(nèi)的流場，由于矩形網(wǎng)格不能很好地擬合這一區(qū)域的復(fù)雜岸線，加之采用調(diào)和常數(shù)進行潮流數(shù)值模擬，潮位和潮流的計算結(jié)果缺乏一定的精度，而且重要的是沒有詳細(xì)的討論這一區(qū)域的海底摩阻系數(shù)的計算及其對水動力場的影響。采用基于有限體積法的三維FVCOM模型，利用寧波市海洋環(huán)境檢測中心2009年對該區(qū)域進行的枯水期和豐水期水文調(diào)查實測資料，建立了寧波舟山及其近海的三維水動力模型，發(fā)現(xiàn)海底摩阻系數(shù)對這一區(qū)域的水動力場影響較大，并采用Koutitas公式計算這一區(qū)域的海底摩阻系數(shù)，經(jīng)驗證計算值與實測值符合良好。

1 FVCOM水動力模型

1.1 水動力控制方程組

FVCOM模型的全稱是非結(jié)構(gòu)網(wǎng)格的有限體積法海洋模式。該模型采用非結(jié)構(gòu)化的三角形網(wǎng)格，能夠精確的擬合復(fù)雜曲率的岸界；采用有限體積法，能夠在整個計算區(qū)域內(nèi)更好的保證了各物理量的守恒。該模型采用垂向靜壓假定，垂直方向采用σ坐標(biāo)變換，使垂向沿水深的分層更趨合理。σ坐標(biāo)下的連續(xù)性方程、動量方程和溫鹽方程如式（1）～式（5）

式（1）～式（3）分別是連續(xù)性方程、x方向的動量方程和y方向的動量方程，式（4）、式（5）分別是溫度和鹽度方程。其中u，v，ω分別是σ坐標(biāo)系下的水平速度和垂向速度。σ取值從海底處-1到海表面處0。此模型采用新的MY-2.5湍流閉合模型求解垂向湍流黏滯系數(shù)和擴散系數(shù)。具體的湍動能擴散方程和詳細(xì)描述見文獻(xiàn)［6］，具體的數(shù)值離散格式詳見文獻(xiàn)［7］，這里不再一一贅述。

1.2 模型定解條件

1.2.1 自由表面條件

運動學(xué)邊界條件

動力學(xué)邊界條件

式中：τsx，τsy分別為風(fēng)應(yīng)力矢量在x和y方向上的分量。

溫鹽邊界條件

式中：Qn（x，y，t）為表面靜熱通量；SW（x，y，0，t）為在海表面處短波輻射通量；cp為海水比熱系數(shù)。

1.2.2 海底條件

運動學(xué)邊界條件

動力學(xué)邊界條件

式中：τbx，τby分別為底摩擦應(yīng)力在x和y方向上的分量。

溫鹽邊界條件

1.2.3 側(cè)邊界條件

側(cè)邊界可分為岸邊界與開邊界兩種。

（1）岸邊界：海岸線或沿岸建筑物即為不透水的岸邊界，水質(zhì)點沿岸線切向自由滑移，其邊界條件為

式中：n為岸邊界的法向量。

（2）開邊界：FVCOM模型中有2種開邊界條件，即采用調(diào)和常數(shù)或?qū)崟r水位，本模型中采用實測潮位資料作為開邊界。

2 模型的計算

2.1 計算區(qū)域與網(wǎng)格

由于杭州灣漲落潮流受到寧波舟山群島的影響，并考慮已有實測潮位資料的站點位置，確定計算范圍的開邊界為北起蘆潮港，往東至嵊山，然后轉(zhuǎn)向西南，延伸至漁山列島，最后與三門灣健跳相連；另外，將余姚和澉浦相連作為杭州灣內(nèi)另一開邊界。模型經(jīng)緯度范圍為120.902°E～122.800°E，28.885°N～30.855°N，計算范圍如圖1所示。網(wǎng)格數(shù)量為165 450，節(jié)點總數(shù)為86 659。網(wǎng)格步長最大3 500 m，最小為150 m，計算網(wǎng)格如圖2所示。

圖1 計算區(qū)域及驗證站位圖（Δ為潮位站，#為潮流站）Fig.1 Calculated range and distribution of stations

圖2 計算網(wǎng)格圖Fig.2 Gird of calculated range

2.2 模型計算參數(shù)的設(shè)置

本模型計算中垂向均勻分5層，時間步長取為1.0 s?？紤]了研究區(qū)域內(nèi)2009年枯水期（2月25日～3月11日）和豐水期（6月21日～7月5日）各15 d的潮位變化，開邊界上的潮位根據(jù)主要潮位測站的實測潮位資料進行線性插值得到，所有測站的潮位高程基準(zhǔn)面和水深基準(zhǔn)面統(tǒng)一采用國家85高程。由于缺乏溫度鹽度的實測資料，本文中溫度鹽度不予計算。

2.2.1 海底摩阻系數(shù)的設(shè)置與討論

在研究中發(fā)現(xiàn)這一區(qū)域的水動力場對海底摩阻系數(shù)甚為敏感，特在此做詳細(xì)討論。FVCOM模型中采用的海底摩阻系數(shù)計算表達(dá)式為

式中：k為卡門常數(shù)，一般取為0.4；zab為最底部水層的厚度；zo為海床糙率高度，其值的大小主要與床面泥沙的粒徑、級配和床面幾何形狀有關(guān)。因此，Cd應(yīng)該是與水深和底質(zhì)類型有關(guān)的一個參數(shù)。由式（14）計算出的Cd≥0.002 5，而對于淤泥質(zhì)海床，其zo一般很小，摩阻系數(shù)遠(yuǎn)小于0.002 5，在這種情況下，式（14）中取max值計算Cd并不合適。

在模擬過程中，針對杭州灣內(nèi)1#測流站局部區(qū)域做了0.002 5和0.000 8兩種摩阻系數(shù)下水動力的計算，將計算結(jié)果與實測流速大小進行了比較分析（圖3）。由圖3可以看出，摩阻系數(shù)取為0.002 5時計算流速與實測值比較明顯偏小，摩阻系數(shù)取為0.000 8更為合適。同理，對象山縣外海至象山港內(nèi)做了0.002 5和0.001 0兩種摩阻系數(shù)下水動力計算，將象山港內(nèi)烏沙山站的計算潮位值與實測值比較（圖4）。由圖4可以看出，摩阻系數(shù)的選取對象山港這一區(qū)域的潮位過程曲線的影響也較大，當(dāng)Cd=0.002 5時象山港內(nèi)烏沙山站的計算潮位較實測值滯后現(xiàn)象較為明顯，而Cd=0.001 0時，滯后現(xiàn)象并不存在。

圖3 兩種流速過程驗證曲線的對比Fig.3 Comparison between two tidal current processes

圖4 兩種潮位過程驗證曲線的對比Fig.4 Comparison between two tidal level processes

因此，對于淤泥質(zhì)海岸，其海底摩阻系數(shù)一般很小且對水動力場的準(zhǔn)確模擬有著很大影響，選取計算公式時應(yīng)慎重。對多種摩阻計算公式的模擬結(jié)果進行比較，選取對數(shù)形式的摩阻系數(shù)計算式（Koutitas，1988）

式中：k為卡門常數(shù)，取為0.4。D為當(dāng)?shù)厮钪?；zo為海床糙率高度，其值的大小主要與床面泥沙的粒徑和級配有關(guān)。在模擬過程中，根據(jù)流場計算值與實測值的比較，通過修改參數(shù)zo來改變摩阻系數(shù)的大小。本文的數(shù)學(xué)模型中，根據(jù)計算區(qū)域內(nèi)表層底質(zhì)采樣的粒度分析數(shù)據(jù)顯示，其中值粒徑絕大多數(shù)都在0.01 mm左右，屬于粘土質(zhì)粉砂類，海床基本為淤泥質(zhì)一類。經(jīng)研究后，文中zo取為0.002 mm，計算的Cd值隨水深地形在0.000 6～0.001 4范圍內(nèi)變化。

另外，海底的摩阻系數(shù)除反映海床粗糙度外，還包括了其他阻力因素對水流的綜合影響，所以已不是原有意義的海床糙率，應(yīng)當(dāng)將其看成是一個綜合阻力的影響因子。在水動力數(shù)值模擬計算中，由于計算范圍的限制，對海岸工程建筑物會有不同程度的近似，此時可以根據(jù)需要在局部區(qū)域加大或減小摩阻系數(shù)來考慮流場阻力情況。

3 模型的驗證

3.1 潮位計算驗證

首先對研究區(qū)域的枯水期和豐水期各進行了15 d的水動力場計算。計算的初始條件為：當(dāng)t=0時，u=v=w=ζ=0，其中 u，v分別為水平流速的x方向的分量和y方向的分量；w為垂向流速；ζ為未擾動海面上的潮位。整個計算域由靜止開始計算，大約6 h后全場流態(tài)達(dá)到穩(wěn)定。計算時按小時輸出驗證潮位站的潮位變化過程，用于對大小潮進行潮位驗證。限于篇幅，文中只給出部分潮位測站豐水期大小潮的潮位驗證過程曲線（圖5）。由圖5可知，豐水期時，寧波舟山及其近海海域的潮汐日不等現(xiàn)象在大潮期間表現(xiàn)的更為顯著；枯水期（文中未給圖）間，潮汐日不等現(xiàn)象在小潮期間表現(xiàn)的更為顯著。其中高高潮與低高潮差值又以杭州灣和三門灣為最大，如洋山和烏沙山大潮期間這種差值達(dá)到近70～90 cm，桃花島與北侖一帶這種差值也在50 cm左右。從計算域內(nèi)的潮差上來看，以杭州灣灣內(nèi)的潮差最大，豐水期能達(dá)到8.0 m左右；其次是三門灣灣內(nèi)，豐水期潮差達(dá)到6.0 m；以鎮(zhèn)海北侖的潮差為最小。另外，南韭山站的潮位較蝦峙門站的潮位提前20 min左右，這說明寧波舟山近海的潮汐由東南向西北方向傳播。

3.2 潮流計算驗證

圖5 潮位過程驗證曲線Fig.5 Validation of tidal level process

計算時按小時輸出驗證潮流站的流速流向變化過程，用于對大小潮進行流速流向驗證。豐水期大潮期間的流速流向驗證過程如圖6所示，驗證結(jié)果較為理想，只有2#站流速落潮流速稍大，2#測站位于渦旋附近，取與其相鄰2個網(wǎng)格的計算速度，則沒有這種落潮流速偏大的情況，因此可能是測流中測點位置移動的緣故。從圖6可以看出，1#和2#潮流站的流速較大，峰值在1.2～2.0 m/s，其他的站點流速較小。從1#潮流站的流速流向可以看出，杭州灣漲潮流經(jīng)過洋山以后以正西方向向灣內(nèi)挺進；寧波舟山群島內(nèi)部海域水體往復(fù)流明顯，其流向基本與岸線平行，由于島嶼的阻擋和地形約束，漲落潮流在途經(jīng)這些水道時流速會陡然增大，最大流速在漲潮時偏向北岸，落潮時偏于南岸，這主要是由水道的開口方向和科氏力作用而引起的。

3.3 三維水動力分析

3.3.1 流場漲落過程及其特性分析

從計算的流場結(jié)果分析，外海潮波沿東南至西北方向向?qū)幉ㄖ凵浇逗Ｓ騻鞑?。?dāng)潮波傳過南韭山后，一部分直接進入象山港海域，另外一部分又分為兩股，一股以較強流速穿越舟山群島途經(jīng)大衡山和大戢山向西進入杭州灣，一股穿過寧波舟山群島內(nèi)部海域（鎮(zhèn)海北侖一帶）的諸條水道繼續(xù)北上進入杭州灣。計算全域內(nèi)漲潮流場和落潮流場、穿山半島附近漲潮流場和落潮流場如圖7所示。

計算結(jié)果表明，杭州灣內(nèi)的流向基本與其灣口走向一致，漲潮時平均流向為230°～300°，落潮時平均流向為60°～110°，漲落潮平均流向約為180°，近乎往復(fù)流性質(zhì)。杭州灣北岸水域受漲潮流控制，水量為凈輸入；南岸水域受落潮流控制，水量為凈輸出。這一特點與杭州灣內(nèi)高懸浮泥沙的成因即長江口向杭州灣內(nèi)凈輸沙一致。在寧波舟山群島內(nèi)部的金塘水道和螺頭水道內(nèi)，漲落潮流形成明顯的往復(fù)流，漲落流場結(jié)構(gòu)基本一致，只是流動方向相反。從圖7-c、圖7-d中可以看出，漲潮時流速分布以北岸較大，落潮時流速分布以南岸較大。另外，在寧波舟山群島內(nèi)部眾多島嶼附近，漲落潮時流場出現(xiàn)渦旋，與實測海流資料相符，這主要是由于群島地形約束，過水?dāng)嗝嬲魉俅蠖鸬摹Ａ硗?，潮波自外海傳入象山港口門后，由于不斷受地形影響及邊界的反射作用，逐漸由前進潮波向駐波轉(zhuǎn)變，其潮差由口門向灣內(nèi)逐漸增大。高潮位沿程抬高，低潮位沿程降低。

3.3.2 流場垂向結(jié)構(gòu)分析

圖6 流速流向驗證曲線Fig.6 Validation of tidal current process

圖7 三維水動力場計算結(jié)果Fig.7 Simulated results of 3D flow field

選取有代表性的三點即1#、2#和4#站的實測垂向流速結(jié)構(gòu)進行分析（圖8）。處在杭州灣的1#站分層流速隨著水深的增大而逐漸減小，并且流速減小的幅度較大，表層比底層大近80 cm/s，可能由于杭州灣水深較淺（大部水深在10～15 m），流速受底部阻力因素影響所致；漲落時刻分層流向基本一致。2#站位于北侖附近的深水航道內(nèi)，流速大小介于120～240 cm/s，其漲潮時刻流速大于落潮時刻，且漲潮時刻分層流速隨水深的增大呈現(xiàn)先增大后減小的趨勢，即最大流速并不是出現(xiàn)在流場表層；落潮時刻分層流速隨水深的增大而減小，表層與底層流速差別較??；漲落時刻分層流向基本一致，與岸線平行。4#站位于象山港外的牛鼻山水道附近，此處是外海水進出象山港的主要通道，4#漲落時刻的分層流向基本一致；漲落潮時刻分層流速都隨水深的增大而減小，但是漲潮時這種減小程度并不明顯，底層比表層小20 cm/s左右；而落潮時的減小程度較大，底層比表層小80 cm/s左右；這可能受地形和過水面積的約束所致，漲潮時外海水進入水道，過水面積急劇減小；落潮時象山港內(nèi)的水退入外海，過水面積急劇增大。由上可知，水深變化劇烈的深水航道內(nèi)的垂向流速分布有別于一般開闊水域或地形平坦水域。

圖8 測站實測垂向流速流向Fig.8 Measured tidal current at stations

4 小結(jié)

在準(zhǔn)確擬合寧波舟山及其近海復(fù)雜岸線和改進底部摩阻系數(shù)計算公式的基礎(chǔ)上，利用FVCOM模型對寧波舟山及其近海海域進行了三維水動力數(shù)值模擬。通過枯水期和豐水期實測水文數(shù)據(jù)和計算值的比較，驗證了該模型的精確計算效果。這是首次將實測潮位資料和FVCOM模型應(yīng)用于寧波舟山大范圍海域內(nèi)的三維潮汐潮流數(shù)值模擬，從計算效果看，該模型在研究近岸水動力學(xué)上有著獨特的優(yōu)越性。除模型的成功應(yīng)用之外，本文得出以下結(jié)論：

（1）寧波舟山及其近海海域內(nèi)海床基本為淤泥質(zhì)，底部摩阻系數(shù)很小且其大小對水動力場影響較大，表現(xiàn)為流速明顯偏小和局部潮位明顯滯后，F(xiàn)VCOM中原摩阻計算公式并不適用于這一區(qū)域。文章經(jīng)研究采用Koutitas公式計算摩阻系數(shù)，經(jīng)模型驗證該公式是合理的。

（2）杭州灣內(nèi)的流向基本與其灣口走向一致，漲落潮平均流向夾角約為180°，近乎往復(fù)流性質(zhì)。寧波舟山海域內(nèi)潮汐日不等現(xiàn)象顯著，在枯水期間以小潮表現(xiàn)更為顯著，豐水期間以大潮表現(xiàn)更為顯著，日不等現(xiàn)象由外海至近岸逐漸增強。

（3）寧波舟山群島內(nèi)部海域（鎮(zhèn)海北侖一帶）漲落潮流形成明顯的往復(fù)流，漲落潮流場基本結(jié)構(gòu)一致。

（4）受地形的影響及邊界的反射作用，象山港內(nèi)的潮差由口門向灣內(nèi)逐漸增大，高潮位沿程抬高，低潮位沿程降低。

（5）對于三維流場垂向分層流速結(jié)構(gòu)，水深變化大的深水航道內(nèi)，其流速分布隨水深先增大后減小，這有別于一般開闊水域或水深變化小的水域。

本次數(shù)值模擬的精度達(dá)到了要求，取得了好的模擬效果，但是對寧波舟山及其近海這一大范圍區(qū)域來說，地形復(fù)雜，深水航道處較多而且水深變化劇烈，其海底底部糙率系數(shù)不僅僅只是和水深、海底糙率有關(guān)，因此這一區(qū)域內(nèi)海底摩阻系數(shù)的計算和選取仍然是一個有待于進一步研究的問題。

致謝：本文采用的實測潮位潮流資料均來源于寧波市海洋環(huán)境監(jiān)測中心2009年的海洋水文調(diào)查，特致感謝。

［1］曹德明，方國洪.杭州灣潮汐潮流的數(shù)值計算［J］.海洋和湖沼，1986，17（2）：93-101.

CAO D M，F(xiàn)ANG G H.A Numerical computation of the tides and tidal currents in Hangzhou bay［J］.Oceanalogia et limnologia sinica，1986，17（2）：93-101.

［2］曹德明，方國洪.杭州灣和錢塘江潮波的聯(lián)合數(shù)值模型［J］.海洋學(xué)報，1988，10（5）：521-530.

［3］李身鐸，顧思美.杭州灣潮波三維數(shù)值模擬［J］.海洋與湖沼，1993，24（1）：7-15.

LI S D，GU S M.3D Numerical simulation of tidal waves in Hangzhou bay［J］.Oceanalogia et limnologia sinica，1993，24（1）：7-15.

［4］陳倩，黃大吉，章本照.浙江近海潮汐潮流的數(shù)值模擬［J］.海洋學(xué)報，2003，25（5）：9-19.

CHEN Q，HUANG D J，ZHANG B Z.Numerical simulation of tide and tidal currents in the seas adjacent to Zhejiang［J］.Acta Oceanologica Sinica，2003，25（5）：9-19.

［5］壽瑋瑋，吳建政，胡日軍，等.舟山群島附近海域三維水動力數(shù)值模擬［J］.海洋地質(zhì)動態(tài)，2009，25（11）：1-9.

［6］吳倫宇.基于FVCOM的浪、流、泥沙模型耦合及應(yīng)用［D］.青島：中國海洋大學(xué)，2009.

［7］朱建榮.海洋數(shù)值計算方法和數(shù)值模擬［M］.北京：海洋出版社，2003.

［8］宋志堯.平面潮流數(shù)值模擬中底床摩阻系數(shù)的修正［J］.水動力學(xué)研究與進展，2001，A16（1）：56-61.

SONG Z Y.A correction of formula for bed resistance coefficient in numerical simulation of plane tidal current［J］.Journal of hydrodynamics，2001，A16（1）：56-61.

［9］Koutitas C G.Mathematical models in coastal engineering［M］.London：Pentech Press，1988.

3D numerical simulation of tide and tidal currents in sea adjacent to Ningbo and Zhoushan

XIONG Wei1，LIU Bi-jin1，2，SUN Zhao-chen1，LIANG Shu-xiu1，ZHANG Yi-fei2
（1.State Key Laboratory of Coastal and Offshore Engineering，Dalian University of Technology，Dalian116023，China；2.Hangzhou Investigation Design Institute of Ocean Engineering，State Oceanic Administration，Hangzhou310012，China）

By means of FVCOM model，the tide and tidal currents in the seas adjacent to Ningbo and Zhoushan were simulated.Furthermore，the dynamic characteristics of tide were analyzed.The grid system of FVCOM model was unstructured triangle mesh，which solved the accurate fitting of complex of Ningbo and Zhoushan shoreline.Based on the property of seabed in the calculation zone，Koutitas formula was used to instead of the original calculation method of drag coefficient in FVCOM.Comparing the computed values with those tidal observatories，the improved formula proved to be reasonable and the two values are in good agreement.The computed model can be used for 3D pollutant concentration and sediment transport calculation.

drag coefficient；tide;tidal currents；numerical simulation；FVCOM model

P 731.23；O 242.1

A

1005-8443（2011）06-0399-09

2011-06-20；

2011-07-13

寧波市海洋泥沙水動力環(huán)境調(diào)查與建設(shè)用海區(qū)域評價項目（1102-082028）；水沙科學(xué)及水災(zāi)害防治湖南省重點實驗室開放課題基金資助項目（2011SS01）

熊偉（1988-），男，湖北省荊州市人，碩士研究生，主要從事海岸水動力學(xué)及泥沙運動研究。