劉文龍，毛晨浩，陳 維，蔣茗韜，宋春梅

（1.浙江海洋大學(xué)海洋工程裝備學(xué)院，浙江舟山 316022；2.上海大學(xué)力學(xué)與工程科學(xué)學(xué)院，上海 200444）

安海灣是位于福建省東南部泉州市轄區(qū)圍頭灣內(nèi)的一個(gè)狹長(zhǎng)半封閉型小海灣，東側(cè)為晉江市、西側(cè)為南安市。水陸交通方便，東南方向與金門(mén)島隔海相望。然而安海灣內(nèi)灘涂面積較多，尤其在安海灣北半部，低平潮時(shí)幾乎全是潮灘露出，水深條件不足。受水深條件限制，僅在石井作業(yè)區(qū)有萬(wàn)噸級(jí)以上深水泊位和建設(shè)條件，其余均為小泊位，東石、安海等現(xiàn)有碼頭均需乘高潮進(jìn)出港。灣內(nèi)有陸源排污口3 個(gè)，分別位于水頭、安海及東石。海灣周?chē)€有數(shù)條排洪溝，積有大量的廢棄物、垃圾、泥沙等，洪水季節(jié)會(huì)隨暴雨沖刷進(jìn)入安海灣。持續(xù)惡化的生態(tài)環(huán)境以及較弱的水體交換能力，已經(jīng)嚴(yán)重影響到安海灣的港口運(yùn)輸業(yè)發(fā)展及海灣納潮量和環(huán)境容量。此次疏浚工程旨在通過(guò)對(duì)灣底泥沙的深挖疏浚，清除安海灣已受污染的灣底沉積物，徹底阻止其二次污染及擴(kuò)散，擴(kuò)大安海灣的納潮量與資源環(huán)境容量。

許多學(xué)者對(duì)地形變化情況下水動(dòng)力條件變化進(jìn)行了研究，孫芹芹等[1]通過(guò)對(duì)不同疏浚方案的對(duì)比，從資源環(huán)境效益方面為安海灣疏浚方案的選擇提供了參考依據(jù)。崔崢等[2]通過(guò)研究認(rèn)為安海灣清淤后能夠增加納潮量、改善石井作業(yè)區(qū)動(dòng)力及安海灣水質(zhì)。張雪松[3]對(duì)半封閉海灣的水動(dòng)力環(huán)境進(jìn)行了研究，龔旭東等[4]從潮流、納潮量和水交換這幾個(gè)方面探討灣內(nèi)圍填海對(duì)東山灣水動(dòng)力環(huán)境的影響，得到隨著填海面積的增加，影響范圍及程度不斷增大的結(jié)論。錢(qián)明霞等[5]對(duì)南港河段地形與水動(dòng)力的相互響應(yīng)特征進(jìn)行了研究。何昆等[6]通過(guò)波浪潮流耦合作用下的數(shù)學(xué)模型，得到清淤工程的實(shí)施對(duì)唐島灣內(nèi)水動(dòng)力及泥沙輸移有一定影響結(jié)論。林國(guó)堯等[7]發(fā)現(xiàn)在研究區(qū)內(nèi)歐拉余流表現(xiàn)為多個(gè)順時(shí)針與逆時(shí)針的渦流，渦流分布與地形具有較好的對(duì)應(yīng)關(guān)系，摩擦力在渦流的發(fā)育中起著重要作用。

由于疏浚工程完成后將要開(kāi)展安海灣內(nèi)南安一側(cè)的五里橋畔休閑慢道景觀深化項(xiàng)目，該項(xiàng)目北接五里橋文化公園，南接石井鎮(zhèn)，全長(zhǎng)約為7.2 km，疏浚工程會(huì)導(dǎo)致安海灣內(nèi)的地形條件被改變，考慮到疏浚工程導(dǎo)致安海灣海域水動(dòng)力環(huán)境改變對(duì)后續(xù)安海灣工程項(xiàng)目實(shí)施對(duì)該項(xiàng)目以及后續(xù)可能開(kāi)展的其他項(xiàng)目可能產(chǎn)生的影響，在已有的水文觀測(cè)資料的基礎(chǔ)上，本文采用MIKE21 Flow Model 建立安海灣潮流模型方法對(duì)灣內(nèi)工程前后潮汐、潮流特性展開(kāi)研究。

1 工程概況

安海灣口門(mén)寬度0.8 km，海灣東西寬1.88 km，南北長(zhǎng)9 km，屬狹長(zhǎng)半封閉型小海灣，面積13.13 km2。低平潮時(shí)灘涂面積9.79 km2，占海灣總面積的75%，海灣北側(cè)有大約100 hm2的水草灌木灘涂分布，水深較大處位于口門(mén)西側(cè)以及進(jìn)入灣內(nèi)的狹長(zhǎng)航道處。疏浚區(qū)總面積約7.9 km2，清理沙場(chǎng)區(qū)面積約0.24 km2。疏浚區(qū)北至安海灣灣頂，南至安海灣口門(mén)白沙頭附近，共分為4 個(gè)區(qū)，其中疏浚區(qū)范圍如圖2 所示，疏浚區(qū)1 的面積為7.45 km2，疏浚底高程為-11.956 m（基準(zhǔn)面為85 國(guó)家高程，下同）；疏浚區(qū)2 的面積為0.32 km2，疏浚底高程為-4.956 m；疏浚區(qū)3 的面積為0.12 km2，疏浚底高程為-2.956 m；疏浚區(qū)4 的面積為0.03 km2，疏浚底高程為-0.956 m。晉江一側(cè)疏浚區(qū)邊界線（坡頂線）距現(xiàn)有海堤前沿線的距離為50 m，南安一側(cè)疏浚區(qū)邊界線（坡頂線）距現(xiàn)有海堤前沿線的距離為140～270 m，疏浚邊坡坡度取1:10，總疏浚量7.7×107m3。

圖2 疏浚區(qū)域Fig.2 Dredging area

2 數(shù)學(xué)模型的建立及驗(yàn)證

本文采用丹麥水力學(xué)研究所（DHI）研發(fā)的MIKE21 Flow Model 模型[8]中的水動(dòng)力模塊（HD）建立安海灣潮流數(shù)學(xué)模型。模型通過(guò)有限體積法進(jìn)行離散，顯式歐拉法求解，具備計(jì)算速度快、易收斂、計(jì)算精度高等優(yōu)點(diǎn)[9]。

2.1 潮流模型控制方程

潮流模型建立在Navier-Stokes 方程[10]的基礎(chǔ)上。在笛卡爾坐標(biāo)系下，通過(guò)對(duì)三維水平動(dòng)量方程和連續(xù)方程沿水深積分得到二維淺水方程[11]

式中：t 為時(shí)間，s；x，y 為笛卡爾坐標(biāo)；h 為總水深，m，η 為水位高程，m，d 為靜止水深，h=η+d，m；f 為科氏力系數(shù)，s-1；ρ 為水的密度，kg·m-3；u，v 分別為x，y 方向的水深平均流速，m·s-1；τsx，τsy分別為x，y 方向的表面風(fēng)應(yīng)力N·m-2；τbx=ρgv，其中M 為曼寧數(shù)，m1/3·s-1，τbx，τby分別為x，y 方向的底床切應(yīng)力，N·m-2；側(cè)應(yīng)力Tij包括黏性摩擦和紊動(dòng)摩擦等，通過(guò)水深平均流速沿梯度方向的變化計(jì)算如下：

亞網(wǎng)格尺度渦流黏度由Samagorinsky 公式確定：

式中：cs為常數(shù)，在本文中取值為0.28；l 為特征長(zhǎng)度，m；變形率

2.2 計(jì)算網(wǎng)格、邊界條件及參數(shù)設(shè)定

模型計(jì)算范圍如圖3 所示，范圍約4 361 km2。模型對(duì)計(jì)算區(qū)域采用非結(jié)構(gòu)網(wǎng)格進(jìn)行漸變剖分，網(wǎng)格的節(jié)點(diǎn)總數(shù)為10 810 個(gè)，網(wǎng)格單元總數(shù)為19 810 個(gè)，在安海灣區(qū)域進(jìn)行了局部加密，最小的網(wǎng)格步長(zhǎng)約33 m，至外海網(wǎng)格逐漸擴(kuò)至9 470 m。

模型開(kāi)邊界包括東、南、西外海開(kāi)邊界以及河口邊界（見(jiàn)圖3），潮流模型的外海開(kāi)邊界采用潮位過(guò)程控制，潮位過(guò)程由MIKE21 軟件包中的全球潮汐預(yù)報(bào)模型提供,流量根據(jù)實(shí)測(cè)資料河口1 取10 m3·s-1，河口2 取15 m3·s-1，河口3 取15 m3·s-1。研究區(qū)域隨著潮位漲落，部分區(qū)域在落潮時(shí)岸灘出露，因此邊界區(qū)域采用動(dòng)邊界處理潮間帶干濕交換過(guò)程，干水深hdry、淹沒(méi)水深hflood和濕水深hwet分別取0.005 m、0.05 m 和0.1 m，曼寧數(shù)根據(jù)底部泥沙粒徑和水深分布取值26～120。在時(shí)間上采用一階顯性歐拉法進(jìn)行離散[12]，一階顯式歐拉法有較高的計(jì)算效率，但時(shí)間步長(zhǎng)，需受克朗數(shù)的限制。本文在保證計(jì)算精度的前提下，選擇計(jì)算快速的一階顯式歐拉法。模型的計(jì)算時(shí)間步長(zhǎng)在0.1～30 s 的范圍內(nèi)自行調(diào)節(jié)，從而保證模型計(jì)算時(shí)始終保證克朗數(shù)CFL＜1。

2.3 模型驗(yàn)證

安海灣潮流模型潮位驗(yàn)證資料采用文獻(xiàn)[13]中后石以及大嶝驗(yàn)證點(diǎn)，驗(yàn)證時(shí)段為2013 年9 月11 日00:00 至2013 年9 月15 日00:00，潮流驗(yàn)證資料采用文獻(xiàn)[14]中1#和2#驗(yàn)證點(diǎn)，驗(yàn)證時(shí)段為2014 年10月24 日8:00 至2014 年10 月25 日10:00，后石、大嶝、1#和2#驗(yàn)證點(diǎn)位置如圖3 所示。

為了定量的評(píng)價(jià)模型模擬結(jié)果，本文采用WILLMOTT[15]提出的統(tǒng)計(jì)學(xué)方法對(duì)模型的模擬結(jié)果進(jìn)行評(píng)價(jià)：

式中：Mi為模擬值；Di為實(shí)測(cè)值；為實(shí)測(cè)平均值。skill 值代表了實(shí)測(cè)值與實(shí)測(cè)平均值的偏差，1～0.65 為極好，0.65～0.5 為非常好，0.5～0.2 為好，0.2～0 為差，計(jì)算結(jié)果如表1 所示。

圖5 流速流向驗(yàn)證Fig.5 Verification of flow velocity and direction

由表1 可以看出，潮流模型評(píng)價(jià)均為極好，說(shuō)明建立的安海灣潮流模型是合理的，能夠很好的模擬該區(qū)域的潮流特征。

表1 模型效率系數(shù)Tab.1 The index of agreement

3 工程影響分析

為研究安海灣疏浚工程對(duì)工程前后潮流的影響，數(shù)學(xué)模型修改了工程區(qū)域的高程以反映設(shè)計(jì)工況以及因水深變化而改變的曼寧數(shù)，其余計(jì)算參數(shù)與驗(yàn)證模型相同。

3.1 潮流響應(yīng)分析

圖6 為疏浚工程前后典型時(shí)刻流速變化（工程后減工程前），潮流為半日潮。潮流的運(yùn)動(dòng)形式通?？梢苑譃樾D(zhuǎn)流和往復(fù)流兩種，旋轉(zhuǎn)流是潮流運(yùn)動(dòng)的基本形式，一般在外海和開(kāi)闊的海域?yàn)樾D(zhuǎn)流，在近岸區(qū)、海峽、水道、河口等以往復(fù)流為主[16]。潮流運(yùn)動(dòng)形式可依主要分潮M2 的橢圓率K 予以判定，對(duì)D1、D3 以及D4 分析點(diǎn)（圖1）的潮流橢圓率K 如表2 所示，其中D1、D3 橢圓率較小，說(shuō)明這2 個(gè)點(diǎn)位處潮流為往復(fù)流，D4 橢圓率較大，是由于該點(diǎn)位距岸線較近，受岸線影響較大，研究區(qū)域總體表現(xiàn)為順岸往復(fù)流，漲潮流基本指向?yàn)檎狈较?，落潮流基本指向?yàn)檎戏较?。由于疏浚工程使得工程區(qū)域整體得到浚深，使得流速的變化主要發(fā)生在原航道處以及口門(mén)處，航道處流速有了明顯的下降，流速減幅最高達(dá)0.44 m·s-1，灣內(nèi)整體流態(tài)變得更加平穩(wěn)。口門(mén)處由于口門(mén)寬度不變而灣內(nèi)地形浚深，使得流速增幅最高達(dá)0.41 m·s-1。

圖1 工程區(qū)域Fig.1 Project area

表2 調(diào)和分析結(jié)果Tab.2 Results of tidal current harmonic analysis

圖6 典型時(shí)刻流速變化Fig.6 Distribution contours of velocity magnitude change at representative time

工程實(shí)施前后，大潮期灣口斷面H-H（圖1）潮量的變化情況見(jiàn)表3。從表中可以看出，疏浚后，通過(guò)灣口斷面的漲、落潮量均有所增加。疏浚后，灣口斷面漲潮時(shí)潮量增加8.58×106m3，落潮時(shí)潮量增加8.61×106m3。疏浚后，安海灣內(nèi)漲、落潮量增加，這主要是由于疏浚后，水深增加較多而使得納潮量增加。

D1、D2、D3 以及D4 分析點(diǎn)（圖1）潮位過(guò)程線對(duì)比見(jiàn)圖7 和圖8，潮波從口門(mén)向?yàn)硟?nèi)傳播，由于受到底摩擦和地形的反射作用，口門(mén)處與灣內(nèi)相位差逐漸變大（圖7a）。在疏浚工程完成后，由于底摩擦變小，潮波由口門(mén)向安海灣內(nèi)傳播過(guò)程中所形成的相位差減少（圖7b）。由圖8a、圖8b、圖8c 可知，潮流從外海向安海灣傳播的過(guò)程中，漲落潮歷時(shí)不對(duì)稱(chēng)性在疏浚之后得到明顯增強(qiáng)，疏浚過(guò)后的落潮歷時(shí)從圖中可以明顯看出較疏浚前有所增加，相應(yīng)的漲潮歷時(shí)有所減小，振幅較疏浚前有所減小。圖8d 所示D4 點(diǎn)位由原來(lái)的低潮時(shí)潮灘外露，經(jīng)過(guò)疏浚，此處水深增大，改善了這一區(qū)域的通航條件。

圖7 疏浚工程前后D1 與D4 潮位對(duì)比Fig.7 Comparison of D1 and D4 tidal levels before and after dredging project

圖8 疏浚工程前后潮位對(duì)比Fig.8 Tide level comparison before and after dredging project

3.2 余流特征分析

將1 個(gè)固定測(cè)站的25 h 等間距流矢量相加，再除以25 h，剩余的那部分就叫歐拉余流[17]。從圖9 安海灣歐拉余流場(chǎng)對(duì)比圖可以看出，水的凈輸移大部分是朝向?yàn)惩獾摹Ｔ跒硟?nèi)，余流整體強(qiáng)度較弱，疏浚前主要在航道處有朝向外的余流，其余位置預(yù)留較小，在疏浚后航道處余流減小，而在口門(mén)處，由于口門(mén)較窄，且疏浚工程使得灣內(nèi)地形改變較大，口門(mén)處漲落潮時(shí)斷面通量較疏浚前有了很大的提升，余流在口門(mén)處變化較大。從圖中可以明顯看出，疏浚后余流在口門(mén)右側(cè)形成了1 個(gè)渦旋，口門(mén)處的地形以及灣內(nèi)疏浚使得口門(mén)斷面通量增大是渦旋較疏浚前增強(qiáng)明顯的主要驅(qū)動(dòng)力。

圖9 工程實(shí)施前后歐拉余流場(chǎng)Fig.9 Eulerian residual current field without and with works

4 結(jié)論

本文采用MIKE21 Flow Model 建立了安海灣潮流數(shù)學(xué)模型，計(jì)算和分析了疏浚工程前后地形變化對(duì)水動(dòng)力環(huán)境的影響，得出以下主要結(jié)論：

（1）由于疏浚工程使得工程區(qū)域潮灘減少，潮波從口門(mén)向?yàn)硟?nèi)傳播過(guò)程中受到的底摩擦作用減小，使得口門(mén)測(cè)量點(diǎn)至灣內(nèi)測(cè)量點(diǎn)之間的相位差較疏浚前減少明顯。

（2）由于疏浚工程對(duì)工程區(qū)域整體的挖深，安海灣內(nèi)整體的流速較之疏浚前有所下降，尤其在航道處流速下降明顯，口門(mén)處則由于狹窄口門(mén)地形的影響，加之安海灣整體納潮量的增加，使得口門(mén)處流速較疏浚工程前增加明顯。

（3）灣內(nèi)余流場(chǎng)整體變化不大，由于疏浚工程對(duì)工程區(qū)域地形的改變，疏浚區(qū)域底摩擦減少，口門(mén)處渦旋有所增強(qiáng)。