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

李玉婷，宋志堯，李瑞杰，彭國強(qiáng)，陳 鵬

(1.南京師范大學(xué) 虛擬地理環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，南京 210023； 2.江蘇省地理信息資源開發(fā)與利用協(xié)同創(chuàng)新中心，南京 210023；3.大規(guī)模復(fù)雜系統(tǒng)數(shù)值模擬江蘇省重點(diǎn)實(shí)驗(yàn)室，南京 210023；4.河海大學(xué) 海岸災(zāi)害及防護(hù)教育部重點(diǎn)實(shí)驗(yàn)室，南京 210098； 5.河海大學(xué) 環(huán)境海洋實(shí)驗(yàn)室，南京 210098)

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

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

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

1 工程海域概況

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

1.1 氣象條件

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

1.2 潮汐潮流條件

工程海域潮汐以M2分潮為主。工程海域的潮汐性質(zhì)為不規(guī)則半日潮，潮汐不對稱現(xiàn)象明顯，三個臨時潮位站均有淺水分潮存在。工程海域潮流運(yùn)動形式為往復(fù)流。實(shí)測最大漲潮流流速為2.28 m/s，流向?yàn)?16°；最大落潮流流速為2.52 m/s，流向?yàn)?7°；垂向平均的最大漲潮流流速為2.17 m/s，流向?yàn)?37°；垂向平均的最大落潮流流速為2.34 m/s，流向?yàn)?8°。

1.3 泥沙條件

工程海域底質(zhì)以細(xì)顆粒的粉砂和粘土為主，其中粉砂約占69%、粘土約占25%，此外還有少量的砂約占6%。底質(zhì)中值粒徑范圍在3.72～6.51 μm之間，平均為5.17 μm。底質(zhì)中值粒徑的空間分布有東面粗、西面細(xì)的特征。實(shí)測最大含沙量為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 水沙動力學(xué)方程

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

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

動力學(xué)方程包括：

連續(xù)性方程

(1)

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

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

(2)

(3)

式中：g為重力加速度；f=2ωesinφ為科氏參數(shù)，ωe為地球自轉(zhuǎn)角速度，φ為當(dāng)?shù)鼐暥?；εxx、εxy、εyx和εyy為不同方向上渦粘系數(shù)；CZ為謝才系數(shù)。

二維懸沙輸移擴(kuò)散方程為

(4)

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

(5)

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

地形變化方程為

(6)

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

不考慮海面風(fēng)作用,陸邊界

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

(7)

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

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

(8)

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

3 模型建立及水沙模擬

3.1 模型的構(gòu)建

圖3 大、小范圍數(shù)學(xué)模型計算區(qū)域Fig.3 The large and small range of calculated area with mathematical model

模型計算區(qū)域如圖3所示，計算區(qū)域北面以長江口以南、南匯嘴以東的東海海域?yàn)楸遍_邊界，南面以六橫島以南約10 km的東海海域?yàn)槟祥_邊界，東面以嵊山島以東約40 km的東海海域?yàn)闁|開邊界，西面以黃灣鎮(zhèn)南杭州灣口為西開邊界，其他邊界均為固邊界。

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

3.2 模型的驗(yàn)證

圖4 大模型計算區(qū)域網(wǎng)格剖分圖 圖5 小模型計算區(qū)域網(wǎng)格剖分圖Fig.4 Mesh map of large model areaFig.5 Mesh map of small model area

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

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

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

3.3 水動力模擬

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

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

4 工程應(yīng)用計算

4.1 工程方案情況

根據(jù)舟山市交通規(guī)劃設(shè)計院提供舟山北向客運(yùn)樞紐工程工程可行性研究報告(2017年)[10]中的燕窩山陸島交通碼頭工程總體布置圖，擬建的燕窩山碼頭位于舟山市岱山縣岱山島北部，地理位置約北緯30°20′48″，東經(jīng)121°10′23″?；静贾梅桨笧閮煞N如下：

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

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

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

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

具體建設(shè)方案如表1所示。

表1 各方案建設(shè)內(nèi)容Tab.1 Construction content of each plan

4.2 方案實(shí)施后流場變化

4.2.1 工程實(shí)施后潮流運(yùn)動特征

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

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

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

4.2.2 工程建設(shè)前后流速變化比較分析

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

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

4.3 方案實(shí)施后淤積變化

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

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

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

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

4.4 方案實(shí)施后航道水流條件分析

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

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

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

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

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

5 結(jié)論

綜合泥沙數(shù)學(xué)模型計算結(jié)果的主要結(jié)論如下：

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

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

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

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