寧波舟山近海三維潮汐潮流數值模擬

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

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

寧波舟山近海三維潮汐潮流數值模擬

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

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

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

摩阻系數；潮汐；潮流；數值模擬；FVCOM模型

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

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

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

1 FVCOM水動力模型

1.1 水動力控制方程組

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

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

1.2 模型定解條件

1.2.1 自由表面條件

運動學邊界條件

動力學邊界條件

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

溫鹽邊界條件

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

1.2.2 海底條件

運動學邊界條件

動力學邊界條件

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

溫鹽邊界條件

1.2.3 側邊界條件

側邊界可分為岸邊界與開邊界兩種。

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

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

（2）開邊界：FVCOM模型中有2種開邊界條件，即采用調和常數或實時水位，本模型中采用實測潮位資料作為開邊界。

2 模型的計算

2.1 計算區域與網格

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

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

圖2 計算網格圖Fig.2 Gird of calculated range

2.2 模型計算參數的設置

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

2.2.1 海底摩阻系數的設置與討論

在研究中發現這一區域的水動力場對海底摩阻系數甚為敏感，特在此做詳細討論。FVCOM模型中采用的海底摩阻系數計算表達式為

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

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

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

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

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

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

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

3 模型的驗證

3.1 潮位計算驗證

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

3.2 潮流計算驗證

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

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

3.3 三維水動力分析

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

從計算的流場結果分析，外海潮波沿東南至西北方向向寧波舟山近岸海域傳播。當潮波傳過南韭山后，一部分直接進入象山港海域，另外一部分又分為兩股，一股以較強流速穿越舟山群島途經大衡山和大戢山向西進入杭州灣，一股穿過寧波舟山群島內部海域（鎮海北侖一帶）的諸條水道繼續北上進入杭州灣。計算全域內漲潮流場和落潮流場、穿山半島附近漲潮流場和落潮流場如圖7所示。

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

3.3.2 流場垂向結構分析

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

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

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

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

4 小結

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

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

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

（3）寧波舟山群島內部海域（鎮海北侖一帶）漲落潮流形成明顯的往復流，漲落潮流場基本結構一致。

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

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

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

致謝：本文采用的實測潮位潮流資料均來源于寧波市海洋環境監測中心2009年的海洋水文調查，特致感謝。

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

CAO D M，FANG 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］曹德明，方國洪.杭州灣和錢塘江潮波的聯合數值模型［J］.海洋學報，1988，10（5）：521-530.

［3］李身鐸，顧思美.杭州灣潮波三維數值模擬［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］陳倩，黃大吉，章本照.浙江近海潮汐潮流的數值模擬［J］.海洋學報，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］壽瑋瑋，吳建政，胡日軍，等.舟山群島附近海域三維水動力數值模擬［J］.海洋地質動態，2009，25（11）：1-9.

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

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

［8］宋志堯.平面潮流數值模擬中底床摩阻系數的修正［J］.水動力學研究與進展，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

寧波市海洋泥沙水動力環境調查與建設用海區域評價項目（1102-082028）；水沙科學及水災害防治湖南省重點實驗室開放課題基金資助項目（2011SS01）

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