舟山岱山島燕窩山碼頭工程潮流泥沙數值模擬

李玉婷，宋志堯，李瑞杰，彭國強，陳 鵬

(1.南京師范大學 虛擬地理環境教育部重點實驗室，南京 210023； 2.江蘇省地理信息資源開發與利用協同創新中心，南京 210023；3.大規模復雜系統數值模擬江蘇省重點實驗室，南京 210023；4.河海大學 海岸災害及防護教育部重點實驗室，南京 210098； 5.河海大學 環境海洋實驗室，南京 210098)

圖1 燕窩山碼頭工程所在位置示意圖Fig.1 Topographic map of simulated domain in Yanwoshan

近年來，為了開發利用深水岸線資源，修建碼頭和防護工程較為普遍，尤其在受風浪作用明顯的海岸帶，防護工程需要突出于海岸線一定距離才能滿足防浪擋沙的要求，保證船舶靠泊條件[1]。在海區建設碼頭工程與防波堤，需要首先進行泥沙沖淤問題的研究。我國在這方面的研究成果很多。例如，李瑞杰等[2]通過建立潮流泥沙運動二維數學模型模擬河口的懸沙運動沖淤變化，為河口整治工程提供科學依據；宋志堯等[3]通過建立海岸河口二維水沙數學模型研究了長江口懸沙的分布特征和各主要入海通道的分沙比變化，為海岸河口相關工程的建設服務；柯杰等[4]在前人基礎上通過建立多重嵌套的二維潮流泥沙數學模型，對枸杞島擬建防波堤工程所引起的潮流場、含沙量場以及海床沖淤變化進行了分析和預測；曹成林等[5]以煙臺八角海域為例，分別采用數值模擬計算與實測值對比兩種方式分析評價了防波堤工程； 解鳴曉等[6]建立波流共同作用下的泥沙數學模型，研究了連云港口門防波堤建設對年平均含沙量場及進港航道回淤的影響。

岱山島位于舟山群島中部，在北緯30°13′～30°21′，東經121°3′～122°13′之間，北接長江口，西臨杭州灣，東瀕東海，南鄰舟山本島，是舟山市第二大島，面積119.3 km2。岱山島位于長江口外國際航線范圍內，地理位置優越, 具有通江達海的區位優勢。為滿足舟山北向的水路客運需求，擬在岱山島北側建設燕窩山碼頭工程(圖1)，包括碼頭和防波堤。根據燕窩山海港工程的規劃，擬建碼頭為高樁式碼頭，棧橋為透空過水式，碼頭與陸域通過透空棧橋相連，具體方案詳見本文第四部分。本文采用數值模擬計算分析碼頭工程實施后的流場以及泥沙淤積情況, 對現有防波堤平面布置方案對水流和泥沙輸運的影響進行分析，選出合理解決方案。

1 工程海域概況

舟山市岱山縣岱山島屬于北亞熱帶南緣季風海洋型氣候。島上存在沙灘、海涂、低山等多種地貌類型，受冷暖空氣交替影響，災害性天氣頻繁。

1.1 氣象條件

此區域多年平均降水量為873.6 mm，6月份月降雨量最多，達121.9 mm，最大年降水量為128.6 mm(1996年)，全年大于25 mm的降水天數為6.6 d。多年平均風速為6.6 m/s，最大風速為38.6m/s。臺風是影響工程區域的主要災害性天氣系統。

1.2 潮汐潮流條件

工程海域潮汐以M2分潮為主。工程海域的潮汐性質為不規則半日潮，潮汐不對稱現象明顯，三個臨時潮位站均有淺水分潮存在。工程海域潮流運動形式為往復流。實測最大漲潮流流速為2.28 m/s，流向為216°；最大落潮流流速為2.52 m/s，流向為77°；垂向平均的最大漲潮流流速為2.17 m/s，流向為237°；垂向平均的最大落潮流流速為2.34 m/s，流向為78°。

1.3 泥沙條件

工程海域底質以細顆粒的粉砂和粘土為主，其中粉砂約占69%、粘土約占25%，此外還有少量的砂約占6%。底質中值粒徑范圍在3.72～6.51 μm之間，平均為5.17 μm。底質中值粒徑的空間分布有東面粗、西面細的特征。實測最大含沙量為1.874 kg/m3，最小含沙量為0.171 kg/m3。垂向平均含沙量最大值為1.495 kg/m3，最小值為0.342 kg/m3(圖2)。

圖2 岱山島北部海洋工程水文泥沙工程站位Fig.2 Engineering station of northern Daishan island

2 水沙動力學方程

采用建立大、小兩套模型進行計算，大模型為岱山島海域，小模型為碼頭工程局部海域。大模型為小模型提供水動力模型所需要的邊界條件，以保證碼頭工程局部流場計算符合岱山島海域潮流場的整體物理特征。

模型采用垂向平均二維淺水方程，采用 ADI法(Alternating Direction Implicit)和DS法(Double Sweep)格式離散，采用非結構三角形網格和有限體積方法進行數值離散和求解。高樁碼頭支撐樁對工程海域流場阻流作用的模擬采用等效阻力的方法。

動力學方程包括：

連續性方程

(1)

式中：d=h+η為總水深，η為水位，h為靜水深；t為時間；u、v為垂線平均流速分別在x、y方向上的分量。

水深平均x、y方向的動量方程分別為

(2)

(3)

式中：g為重力加速度；f=2ωesinφ為科氏參數，ωe為地球自轉角速度，φ為當地緯度；εxx、εxy、εyx和εyy為不同方向上渦粘系數；CZ為謝才系數。

二維懸沙輸移擴散方程為

(4)

式中：C為垂向平均含沙量；εx、εy為擴散系數；Fc為水沙界面通量，由下式給出

(5)

式中：α為沉積系數；ω為泥沙沉速；M為沖刷系數；τb為底部切應力；τe為臨界沖刷切應力；τd為臨界淤積切應力。

地形變化方程為

(6)

式中：γ0為泥沙干密度，Zb為床面高程。

不考慮海面風作用,陸邊界

Vn(x,y,z,t)=0(n為陸邊法向)

(7)

開邊界處給定水位,由東中國海潮波模型[7]計算得到。

懸沙閉邊界條件采用法向零通量邊界條件，懸沙開邊界條件如下式

(8)

式中：Γ為水域開邊界；C*(x,y,t)為已知懸沙濃度。

3 模型建立及水沙模擬

3.1 模型的構建

圖3 大、小范圍數學模型計算區域Fig.3 The large and small range of calculated area with mathematical model

模型計算區域如圖3所示，計算區域北面以長江口以南、南匯嘴以東的東海海域為北開邊界，南面以六橫島以南約10 km的東海海域為南開邊界，東面以嵊山島以東約40 km的東海海域為東開邊界，西面以黃灣鎮南杭州灣口為西開邊界，其他邊界均為固邊界。

大、小模型均采用非結構三角形網格,其中大模型最小網格尺寸為50 m，節點數為15 468，網格數為28 987。小模型最小網格尺寸為5 m , 網格單元34 808個，網格節點17 289個。結果輸出時間步長為1 800 s，計算時段為2017年3月21日00:00時～2017年4月7日00:00；經率定，舟山海域糙率取為0.02；本文利用考慮亞尺度網格效應的Smagororinsky(1963)公式計算水平渦粘系數，取為0.28；由東中國海潮波數學模型確定模型開邊界的外海開邊界條件。大、小模型的計算域及網格剖分見圖4、圖5所示。

3.2 模型的驗證

圖4 大模型計算區域網格剖分圖 圖5 小模型計算區域網格剖分圖Fig.4 Mesh map of large model areaFig.5 Mesh map of small model area

模型計算的水位、流速、振幅以及過程線，均與原型驗證良好，符合《海岸與河口潮流泥沙模擬技術規程》[8]的要求，限于篇幅，不再一一列出。關于含沙量，實測值依據舟山北向客運樞紐(燕窩山)工程水文泥沙測驗分析報告(2017年)[9]，圖6是6個測站的含沙量驗證圖。含沙量的驗證較好地反映工程海域泥沙的輸運特征。

6-a SW1測站6-b SW2測站6-c SW3測站

6-d SW4測站6-e SW5測站6-f SW6測站圖6 各測站大潮含沙量驗證圖Fig.6 Sediment concentration verification

3.3 水動力模擬

工程區潮流數學模型計算流場如圖7所示，以大潮為例，大潮漲潮時，潮流由東南方向進入計算區域，流經衢山島和長涂山時發生分流。水流繞過燕窩島后進入工程海域，近岸海域流速較小。大潮落潮時，來自杭州灣的水流繞過大魚山向西北方向流去，部分水流進入工程海域，進一步繞過燕窩島和鲞蓬山，穿過岱山、長涂山與衢山島中間的海域與南側落潮流匯合進入東海。

7-a 漲急7-b 落急 圖7 小模型工程實施前大潮漲落急流場圖Fig.7 Rapid flow field of tidal fluctuation before the implementation of small model engineering

4 工程應用計算

4.1 工程方案情況

根據舟山市交通規劃設計院提供舟山北向客運樞紐工程工程可行性研究報告(2017年)[10]中的燕窩山陸島交通碼頭工程總體布置圖，擬建的燕窩山碼頭位于舟山市岱山縣岱山島北部，地理位置約北緯30°20′48″，東經121°10′23″。基本布置方案為兩種如下：

方案一：防波堤為折線型布置，堤長1 200 m，防波堤西側折線段長度100 m(方位角為N46°～N226°)，東側折線段長度1 100 m(方位角為N76°～N256°)，擬建碼頭與擬建防波堤的東側折線段以N76°～N256°的方位角平行布置。碼頭陸域位于防波堤南側陸地，碼頭與陸域通過透空棧橋相連，陸域場地位于擬建碼頭東南側(圖8)。

方案二：防波堤軸線布置同方案一。碼頭陸域位于燕窩山西側，碼頭與陸域通過透空棧橋相連，陸域場地位于擬建碼頭東側(圖9)。

圖8 方案一平面布置示意圖Fig.8 Layout of plan I

圖9 方案二平面布置示意圖Fig.9 Layout of plan II

具體建設方案如表1所示。

表1 各方案建設內容Tab.1 Construction content of each plan

4.2 方案實施后流場變化

4.2.1 工程實施后潮流運動特征

圖10為方案一建設后工程海域潮流的大潮漲落急流場圖，圖11為方案二建設后工程海域潮流的大潮漲落急流場圖。由圖可知，大潮漲潮時，水流由燕窩山東側流入工程海域，流經燕窩山后水流發生分離，在燕窩山西側形成回流區。由于防波堤對水流的阻擋，水流由防波堤東側口門處進入港池及近岸海域，水流流量與防波堤建設前相比大幅減少，同時在防波堤后側港池水域形成回流區，流速均小于0.5 m/s。水流流至東墾山時再次發生分離，主流部分繼續向西運動，小部分水流向南側近岸海域運動形成大面積回流。大潮落潮時水流由工程海域西北側流入防波堤附近海域，流經東墾山時主流繼續向燕窩山運動，小部分水流流入近岸淺水海域及港池水域。由于陸地邊界處水深較淺，落急時刻部分邊界區域露出水面。防波堤的建設導致水流流經防波堤時發生分離，并在緊靠防波堤北側邊界區域形成回流區。由于工程海域地形較為復雜，導致防波堤南北兩側流速差異較大。防波堤北側海域流速較大，西側堤頭以北2 km處漲急流速達1.75 m/s，落急流速達2.01 m/s。防波堤以南淺水海域流速基本小于0.5 m/s。中潮與小潮漲落潮流情況與大潮相似，但潮流動力較弱，并呈現落潮流流速略大于漲潮流流速的趨勢。由于岱山北側附近地形復雜，水深梯度大、島嶼礁石眾多，無論漲潮流還是落潮流均會形成許多水流漩渦，工程海域附近潮流仍以往復流流態為主。

10-a 漲急10-b 落急 圖10 方案一工程實施后大潮漲落急局部流場圖Fig.10 The tide flood maximum map after the project I implementation

11-a 漲急11-b 落急 圖11 方案二工程實施后大潮漲落急局部流場圖Fig.11 The tide flood maximum map after the project II implementation

4.2.2 工程建設前后流速變化比較分析

工程前大潮漲潮時，防波堤南、北側特征站位平均流速和最大流速均呈現由東向西遞增的趨勢。大潮落潮時，防波堤南、北側平均流速和最大流速均呈現由西向東遞減的趨勢。圖12為方案一建設后工程海域潮流的大潮漲落急流速等值線圖，圖13為方案二建設后工程海域潮流的大潮漲落急流速等值線圖。根據圖12～圖13和方案一、二實施前后小范圍潮流數學模型的數值模擬結果，以大潮為例，除了部分站位流速有所增大以外，其余流速均有所減小，最大流速變化趨勢和平均流速變化趨勢基本一致。

12-a 漲急12-b 落急13-a 漲急13-b 落急圖12 方案一工程實施后大潮漲落急局部流速等值線圖Fig.12 The velocity contour map after the plan I implementation圖13 方案二工程實施后大潮漲落急局部流速等值線圖Fig.13 The velocity contour map after the plan II implementation

4.3 方案實施后淤積變化

各主要方案年淤積強度如圖14所示，其中“正”表示淤積增加或沖刷減弱，“負”表示淤積減弱或沖刷加強。不同方案工程實施后航道碼頭區域年最大淤積強度和年回淤量見表2。由圖表可知，方案一工程實施后，防波堤附近海域泥沙沖淤強度有所變化，由于防波堤縱向軸線與等深線存在一定夾角，水流流至防波堤附近時流速減小，泥沙形成落淤。落潮時由于防波堤的建設，回流區流速減小，泥沙淤積強度較大；回流區兩側底床有輕微沖刷現象。防波堤東側口門區域由于防波堤的束水作用，導致防波堤堤頭區域流速增大，泥沙沖刷強度增大。防波堤工程實施后航道和港池水域由于開挖疏浚將導致泥沙形成落淤，淤積強度最大處位于航道彎曲段。方案二防波堤的布置方案與方案一一致，僅陸域場地設置在東側燕窩山附近，對工程區域泥沙沖淤強度影響較小，其沖淤強度變化趨勢與方案一類似。

由上述分析可以看出，方案一實施后航道和港池水域年淤積總量為8.65萬m3，比選方案中，方案二防波堤布置形式引起的航道和港池水域回淤量為9.53萬m3。

圖14 不同方案工程實施前后淤積強度變化分布圖Fig.14 Distribution of siltation intensity before and after implementation of different plans

表2 各方案工程實施后航道和港池水域最大淤積強度及年淤積總量Tab.2 The maximum siltation intensity and the total annual sedimentation volume of waterways and harbor basins after each plan implementation

4.4 方案實施后航道水流條件分析

對工程前和各方案工程實施后航道口門段沿航道各特征站位的漲潮和落潮期間橫流流速進行統計分析，結果見表3和表4。防波堤工程實施前，落潮期間的橫流起控制作用。方案一航道口門段橫流流速沿航道向港池方向逐漸減小。在漲潮期間，方案二航道口門段橫流流速沿航道向港池方向逐漸減小。工程實施后，各方案仍是落潮期間的最大橫流大于漲潮期間的最大橫流，落潮期間的橫流起控制作用。

方案一實施后，在漲潮期間，航道口門段橫流流速沿航道向港池方向先增大后減小，在C2處取最大值1.11 m/s，C5處取最小值0.24 m/s；在落潮期間，航道口門段橫流流速沿航道向港池方向逐漸減小，大小為0.56～1.33 m/s；連續兩個潮周期內，橫流大于0.50 m/s的時長沿航道向港池方向逐漸減小，最大為17.50 h，最小為2.50 h。與方案實施前相比，方案實施后漲潮期間航道口門前半段(C1、C2、C3)橫流流速增大，最大增量為0.21 m/s(C2)，后半段減小，最大減小量為0.57 m/s(C5)；漲潮期間各點均減小，最大減小量為0.48 m/s(C5)。

方案二實施后，航道口門段在漲落潮期間橫流流速沿航道向港池方向變化趨勢與方案一相似，漲潮期間在C3處取最大值1.11 m/s，C5處取最小值0.34 m/s；落潮期間在C1處取最大值1.40 m/s，C5處取最小值0.65 m/s；連續兩個潮周期內，橫流大于0.50 m/s的時長沿航道向港池方向逐漸減小，最大為17.50 h，最小為4.00 h。與方案一相比，方案二航道口門段橫流流速除漲潮期間點C1和C2外均大于方案一，最大差值為0.34 m/s，且橫流流速大于0.50 m/s的時長也均大于方案一。

表3 工程前航道特征站位橫流流速Tab.3 Cross current velocity of channel characteristic stations before project implementation

表4 各方案實施后航道特征站位橫流流速Tab.4 Cross current velocity of channel characteristic stations after project implementation

5 結論

綜合泥沙數學模型計算結果的主要結論如下：

(1)工程海域附近潮流以往復流為主，燕窩山和東墾山的北海域流速較大，超過2.0 m/s，近岸海域水深變化梯度較大，5 m等深線以內水流流速較小，在0.6 m/s以內。

(2)各方案防波堤工程實施后，對防波堤附近海域潮位的影響很小。方案一實施后，流速變化最大的特征站位與工程前相比增大32.88%，各特征站位最大流速均在1.0 m/s以內。碼頭前沿特征站位最大橫流流速為0.15 m/s，港池水域回流強度為0.37 m/s，航道特征站位最大橫流流速為1.33 m/s；方案二實施后，碼頭前沿特征站位最大橫流流速為0.12 m/s，港池水域回流強度為0.49 m/s，航道特征站位最大橫流流速為1.40 m/s。

(3)計算結果表明，方案一實施后，航道和港池水域最大年淤積強度為0.71 m/a，年淤積總量為8.65萬m3，方案二實施后, 航道和港池水域最大年淤積強度為0.69 m/a，年淤積總量為9.53萬m3。

綜上所述，從水流、泥沙計算結果看，燕窩山碼頭工程實施后并未對岱山北部海域的整體流態產生較大影響。但根據航道橫流流速變化特征調整航道設計方案，應增大航道設計寬度，減小航道軸線與航道口門附近水流流向的夾角，保證船舶航行安全；從水深條件角度考慮，方案二碼頭位置水深較方案一好,挖泥浚深量較小。從淤積角度考慮，方案一年淤積總量略小于方案二。從文中模擬出的相關流速數據角度考慮，方案二僅碼頭前沿特征站位最大橫流略小于方案一，而方案二港池水域回流強度和航道特征站位最大橫流流速均大于方案一。從建設角度考慮，方案一利用燕窩山遮擋了東向風浪，故減少了防波堤長度，從而減少工程造價；陸域東側為消浪區，碼頭泊穩條件好。而方案二，港區全由防波堤掩護，防波堤較長，工程造價較高；且碼頭前沿受反射波浪影響，不利于泊穩。故綜合評價為方案一優于方案二。