朱志夏，熊偉

(1.中設(shè)設(shè)計(jì)集團(tuán)股份有限公司，江蘇 南京 210014; 2.浙江貴仁信息科技股份有限公司，浙江 杭州 310051)

臺(tái)風(fēng)引起動(dòng)力強(qiáng)大的臺(tái)風(fēng)浪和風(fēng)暴潮流很容易造成近岸淺灘大量泥沙的起動(dòng)，在短時(shí)間內(nèi)造成航道或港池的強(qiáng)淤，使堤防與防波堤等港口航道與海岸及近海工程建筑物遭到嚴(yán)重?fù)p壞，同時(shí)，風(fēng)暴增水也會(huì)給沿海地區(qū)帶來安全隱患，因此，對工程海域臺(tái)風(fēng)浪、風(fēng)暴潮以及泥沙輸移、港口航道強(qiáng)淤等的數(shù)值模擬和預(yù)測提出了更高的要求。為了更加準(zhǔn)確地模擬臺(tái)風(fēng)浪風(fēng)暴潮作用下的泥沙輸移、海床演變、港池航道強(qiáng)淤，有必要進(jìn)行臺(tái)風(fēng)浪風(fēng)暴潮作用下三維泥沙數(shù)值模擬技術(shù)的研究。

有關(guān)波浪潮流共同作用下三維泥沙數(shù)值模擬，已取得了不少成果。朱志夏[1]基于Navier-Stokes方程和質(zhì)量傳輸方程，應(yīng)用多重尺度法、小參數(shù)方法及變分原理，導(dǎo)出了波浪作用下的淺水環(huán)流方程、不恒定非均勻潮流場中隨機(jī)波的折繞射聯(lián)合的數(shù)值模式及波流聯(lián)合作用下的懸沙擴(kuò)散方程。開發(fā)并建立了波流共同作用下二、三維嵌套泥沙數(shù)學(xué)模型，成功應(yīng)用于天津新港附近海域疏浚棄土處理的工程。梁丙臣[2]基于國際著名的水動(dòng)力-生態(tài)-懸浮泥沙耦合模式COHERENS-SED和第3代波浪模式SWAN，對COHERENS-SED進(jìn)行了完善與發(fā)展，最終建立了波流聯(lián)合作用下的三維水動(dòng)力-懸浮泥沙耦合數(shù)學(xué)模型，并應(yīng)用結(jié)構(gòu)化網(wǎng)格二重嵌套的方法，成功應(yīng)用于黃河三角洲濱海區(qū)的潮流和懸沙輸運(yùn)規(guī)律的研究。王紅[3]采用大、中、小三重結(jié)構(gòu)化網(wǎng)格二、三維嵌套模式，結(jié)合第3代波浪模式SWAN和三維水動(dòng)力-多組分泥沙模型EFDC，將模型應(yīng)用于濰坊港海域?qū)嶋H工程問題的研究。由于采用結(jié)構(gòu)化網(wǎng)格，工程中關(guān)注的防波堤附近水動(dòng)力的模擬結(jié)果有些失真。堵盤軍[4]基于第3代波浪模式SWAN和ECOMSED模式的改進(jìn)和優(yōu)化，進(jìn)行了長江口、杭州灣海域三維懸沙數(shù)值模擬的初步研究。研究表明波浪在近岸潮灘區(qū)域?qū)δ嗌硢?dòng)的作用不可忽視。王平[5]根據(jù)考慮流場效應(yīng)的橢圓型緩坡方程，結(jié)合水動(dòng)力泥沙模型FVCOM(an unstructured grid, finite-volume coastal ocean model)，構(gòu)建了近岸大范圍三維波流耦合模型和波流聯(lián)合作用下的近岸物質(zhì)輸運(yùn)模型，分別應(yīng)用于大連灣海域和旅順琥珀灣海域，模擬灣內(nèi)的潮流和物質(zhì)輸運(yùn)。

綜上所述，臺(tái)風(fēng)作用下波浪潮流耦合三維泥沙數(shù)值模擬技術(shù)仍有待進(jìn)一步研發(fā)，例如高分辨率臺(tái)風(fēng)風(fēng)場模擬、波-流耦合計(jì)算、細(xì)顆粒泥沙的沉降與輸移的模擬等等。為此，本文應(yīng)用中尺度大氣模式WRF、第3代波浪模型MIKE21-SW、二維潮流模型MIKE21-HD和三維水動(dòng)力泥沙模型FVCOM，針對臺(tái)風(fēng)作用下泥沙輸移、海床演變、港池航道沖淤，通過二次開發(fā)，建立了臺(tái)風(fēng)作用下波浪潮流耦合三維泥沙數(shù)值模擬系統(tǒng)。著重討論臺(tái)風(fēng)作用下三維泥沙輸移的模擬，主要解決了波流耦合問題和三維模型中增加了波生近岸流、波浪對泥沙與運(yùn)動(dòng)的影響、泥沙的絮凝沉降和高含沙量時(shí)的制約沉降機(jī)制等影響因素，建立了臺(tái)風(fēng)作用下波流耦合作用下三維泥沙數(shù)學(xué)模型，并將模型應(yīng)用于連云港港30萬t級(jí)淺灘深開挖航道工程的研究。

1 三維泥沙數(shù)學(xué)模型FVCOM改進(jìn)

1.1 FVCOM泥沙擴(kuò)散方程

FVCOM采用三維波流耦合作用下的泥沙模型(ROMS_SED)的泥沙計(jì)算模式，將泥沙按粒徑大小分組，考慮黏性沙、非黏性沙，又稱為多組分泥沙模型。其懸沙濃度擴(kuò)散方程為：

(1)

式中：Ci為第i類泥沙的懸沙濃度；AH為水平渦黏性系數(shù)；KH為垂向渦黏性系數(shù)；wi為第i類泥沙的沉降速率。

邊界條件：在自由表面處不考慮大氣沉降的顆粒；在海底邊界處，泥沙的沉降和再懸浮過程分別作為水中泥沙的源和匯處理，即：

式中：Di為懸浮泥沙的沉降通量，Ei為海底泥沙侵蝕通量，Ei的表達(dá)式為：

式中:Qi為泥沙侵蝕系數(shù)；Pb為海底孔隙率；τb為海底剪切率；τci為泥沙起動(dòng)的臨界剪切率。Qi依賴于底部泥沙的物理化學(xué)特性，可在實(shí)驗(yàn)室中測量得到，根據(jù)英國HR實(shí)測資料，對于黏性土，Qi=2×10-4～4×10-3kg/m2s，α=1.16。

1.2 泥沙模型的改進(jìn)

針對連云港港淤泥質(zhì)海岸開敞海域，模型中需要考慮黏性泥沙的絮凝沉降過程和高含沙量時(shí)的制約沉降機(jī)制。當(dāng)泥沙濃度小于制約沉降臨界濃度時(shí)，絮凝沉降速度計(jì)算公式采用彭潤澤等基于長江口泥沙絮凝實(shí)驗(yàn)的表達(dá)式，即：

(2)

當(dāng)泥沙濃度大于或等于制約沉降臨界濃度時(shí)，泥沙沉速計(jì)算公式為：

(3)

式中：Cgel為膠凝濃度，根據(jù)細(xì)顆粒泥沙的相關(guān)研究，計(jì)算中Cgel取為250 kg/m3。

2 三維泥沙模型的建立與驗(yàn)證

2.1 模型建立

數(shù)值計(jì)算采用二維大、中模型和三維小模型三重嵌套的方法，模擬計(jì)算連云港港及附近海域“韋帕”臺(tái)風(fēng)作用下波流耦合三維泥沙輸移過程[6]。三維模型計(jì)算范圍(D03)北起嵐山頭，南至中山河附近，計(jì)算域?yàn)楸本?4°22.1′ ～35°6.6′, 東經(jīng)116°0.0′～120°8.8′，模型采用非結(jié)構(gòu)化三角形網(wǎng)格，網(wǎng)格最小尺度為20 m，位于連云港主港區(qū)和主航道處，網(wǎng)格最大尺度為2 000 m，位于計(jì)算域的深水開邊界處。三維模型的邊界條件由中范圍臺(tái)風(fēng)作用下二維波浪潮流耦合數(shù)學(xué)模型提供。

2.2 “韋帕”臺(tái)風(fēng)作用下波流耦合三維泥沙模型的驗(yàn)證

為了準(zhǔn)確地模擬連云港港及其附近海域“韋帕”臺(tái)風(fēng)作用下波流耦合的三維泥沙輸移過程，根據(jù)天津市氣象科學(xué)研究所應(yīng)用WRF中尺度預(yù)報(bào)模式和美國環(huán)境預(yù)報(bào)中心NCEP歷史再分析數(shù)據(jù)，在天河一號(hào)大型計(jì)算機(jī)上模擬計(jì)算了整個(gè)“韋帕”臺(tái)風(fēng)過程。計(jì)算中采用超高分辨率的三重網(wǎng)格嵌套技術(shù)，嵌套區(qū)域D01、D02、D03的分辨率分別為9、3、1 km，為“韋帕”臺(tái)風(fēng)作用下波浪、潮流、泥沙輸移的模擬計(jì)算提供了準(zhǔn)確的驅(qū)動(dòng)風(fēng)場[6]，計(jì)算時(shí)間為2007年9月17日0時(shí)～9月21日20時(shí)，共117 h。

1)三維含沙量驗(yàn)證。

在三維泥沙數(shù)值模擬計(jì)算中，潮位、潮流、波浪、風(fēng)場等均由二維嵌套模型提供[6]，三維模型外模計(jì)算時(shí)間步長取1.0 s，內(nèi)模計(jì)算時(shí)間步長選取外模的5倍，垂向平均分為10層。底部摩阻系數(shù)范圍取0.002 5～0.003 5。根據(jù)“韋帕”臺(tái)風(fēng)期間實(shí)測的波浪、水文資料[7]，進(jìn)行了“韋帕”臺(tái)風(fēng)作用下波流耦合的三維泥沙數(shù)學(xué)模型的率定與驗(yàn)證。其中1#點(diǎn)位于-3.0 m等深線，2#位于-5.0 m等深線(如圖1所示)，每個(gè)測站分為表下0.5 m、底上0.5 m共3層進(jìn)行觀測。其中，2#測站的分層含沙量過程的驗(yàn)證結(jié)果如圖2所示，從驗(yàn)證結(jié)果可以看出，計(jì)算值與實(shí)測值吻合較好。

圖1 水文觀測點(diǎn)位置Fig.1 Hydrological observation point location

圖2 2#觀測點(diǎn)含沙量驗(yàn)證Fig.2 Verification of sediment concentration at observation point 2#

受“韋帕”臺(tái)風(fēng)影響，2007年9月19日0時(shí)“韋帕”臺(tái)風(fēng)遠(yuǎn)在浙江東海海面上，距離連云港海域較遠(yuǎn)，水體含沙量較小。水體含沙量逐漸增大，到9月20日12點(diǎn)，含沙量增至最大值(如圖2、圖3所示)。然后逐漸減小，主要特征表現(xiàn)為：最大含沙量時(shí)刻落后于最大波高時(shí)刻。根據(jù)三維泥沙數(shù)值模擬結(jié)果得到，海州灣內(nèi)南部的高含沙量水體將隨著風(fēng)暴潮流橫向越過主航道后進(jìn)入徐圩海域，可成為連云港港主航道泥沙回淤的重要來源之一。因此，在“韋帕”臺(tái)風(fēng)期間，海州灣內(nèi)南部的高含沙水體對主航道回淤有比較重要的影響。綜上，本文建立的三維泥沙數(shù)學(xué)模型較好反映了“韋帕”臺(tái)風(fēng)期間連云港海域泥沙分布特征及輸移過程，可為連云港港30萬t級(jí)淺灘深開挖航道工程的設(shè)計(jì)提供技術(shù)支撐。

圖3 2007年9月20日12時(shí)含沙量分布Fig.3 Sediment concentration distribution at 12 o′clock on September 20, 2007

2)航道回淤驗(yàn)證。

根據(jù)2007年9月13日(“韋帕”臺(tái)風(fēng)前)和2007年9月21日(“韋帕”臺(tái)風(fēng)后)連云港港主航道斷面水深測量[8]，2次實(shí)測的航道水深沿程分布得到(如圖4所示)，風(fēng)前航道水深沿程分布在0～1.5 km和1.5～3.0 km呈階梯狀變化，3.0 km以外水深逐漸增大；由于臺(tái)風(fēng)期間航道產(chǎn)生淤積，風(fēng)后航道水深沿程分布在0～3.0 km內(nèi)基本一致，向外水深逐漸加大。

圖4 “韋帕”臺(tái)風(fēng)前后主航道沿程水深Fig.4 Water depth along main channel before and after typhoon Wipha

對比2次水深測量結(jié)果可知，在航道里程0～1.3 km內(nèi)，航道內(nèi)的淤積強(qiáng)度基本一致，在0.3～0.5 m；航道里程1.3～3.0 km的淤積強(qiáng)度較大，約為前者的2倍。模擬計(jì)算得到的“韋帕”臺(tái)風(fēng)作用下主航道回淤與實(shí)測較為一致，模擬計(jì)算的最大實(shí)質(zhì)回淤強(qiáng)度為0.40 m。

3 “韋帕”臺(tái)風(fēng)浪作用下三維泥沙數(shù)學(xué)模型應(yīng)用

3.1 航道工程方案

連云港港30萬t級(jí)航道工程由連云港區(qū)航道、連云港區(qū)外航道、徐圩進(jìn)港航道組成，其中，連云港區(qū)外航道外段寬350 m，底標(biāo)高-23.0 m，外航道內(nèi)段寬310 m，底標(biāo)高-22.5 m；徐圩進(jìn)港航道寬370 m，底標(biāo)高-22.0 m，徐圩港區(qū)內(nèi)航道寬210 m，底標(biāo)高-13.3 m；徐圩二港池航道寬170 m，底標(biāo)高-11.0 m。徐圩港區(qū)防波堤口門分為“平口”、“八字口”，結(jié)合不同港內(nèi)圍墾工程方案布置，分為4個(gè)方案，本文以方案2、方案4為例進(jìn)行分析，具體布置如圖5所示。

圖5 方案2、方案4布置圖Fig.5 Layout of plan 2 and 4

3.2 含沙量計(jì)算分析

圖6～圖11給出了工程方案2、方案4連云港海域在“韋帕”臺(tái)風(fēng)期間代表時(shí)刻的表層和底層含沙量分布情況。2007年9月19日0時(shí)，“韋帕”臺(tái)風(fēng)遠(yuǎn)在浙江東海海面上，距離連云港海域甚遠(yuǎn)，除灌河口含沙量較大外，連云港海域內(nèi)的含沙量均不大(如圖6、圖9所示)。隨著臺(tái)風(fēng)中心逐漸往北移動(dòng)，臺(tái)風(fēng)中心距離連云港越近，臺(tái)風(fēng)浪也隨之加大，導(dǎo)致近岸水體含沙量逐漸增加，波浪掀沙、波流共同輸沙作用明顯，9月20日0時(shí)表層、底層含沙量均有明顯增大(如圖7、圖10所示)。至9月20日12時(shí)，連云港海域含沙量達(dá)到最大(如圖8、圖11所示)，滯后最大臺(tái)風(fēng)浪出現(xiàn)時(shí)刻約5 h[6]。值得注意的是，由于該時(shí)刻的風(fēng)暴潮流為由北向南的沿岸流，海州灣內(nèi)的高含沙水體會(huì)越過連云港區(qū)主航道進(jìn)入徐圩海域，與實(shí)測資料較為相符，即-5.0 m測點(diǎn)的含沙量從最高濃度逐漸減小后，仍然有大幅度增大的趨勢。

圖11 2007年9月20日12時(shí)方案4含沙量分布Fig.11 Plan 4 sediment concentration distribution at 12 o′clock on September 20, 2007

圖10 2007年9月20日0時(shí)方案4含沙量分布Fig.10 Plan 4 sediment concentration distribution at 0 o′clock on September 20, 2007

圖8 2007年9月20日12時(shí)方案2含沙量分布Fig.8 Plan 2 sediment concentration distribution at 12 o′clock on September 20, 2007

圖7 2007年9月20日0時(shí)方案2含沙量分布Fig.7 Plan 2 sediment concentration distribution at 0 o′clock on September 20, 2007

由于工程方案2、方案4港區(qū)圍填面積的不同，港池內(nèi)部的水體含沙量也有差異，但與口門外相比，含沙量都較小。一方面，由于各方案中防波堤均已建成，防波堤的防浪效果較好，港池內(nèi)的波高均不大，口門處波高最大值約為2.2 m，波浪自口門處逐漸遞減，港池內(nèi)最大時(shí)刻波高介于1～2 m，和港池外相比小得多，難以引起港池內(nèi)部泥沙大面積啟動(dòng)和底部高濃度含沙層的產(chǎn)生。另一方面，圍填區(qū)域面積的不同，造成港池內(nèi)的環(huán)流強(qiáng)度均有不同，其特征表現(xiàn)為港池圍填面積越小，環(huán)流強(qiáng)度越大。如工程方案2環(huán)流強(qiáng)度較大的時(shí)刻，港池內(nèi)的環(huán)流區(qū)含沙量濃度稍高(圖6～8)；工程方案4港區(qū)圍填全部形成，環(huán)流強(qiáng)度最弱，港內(nèi)含沙量較小(圖9～11)。因此，徐圩大環(huán)抱防波堤建成后，港池內(nèi)的水體泥沙主要來源于港外輸入和港池內(nèi)大強(qiáng)度環(huán)流造成的淺灘沖刷。

圖9 2007年9月19日0時(shí)方案4含沙量分布Fig.9 Plan 4 sediment concentration distribution at 0 o′clock on September 19, 2007

圖6 2007年9月19日0時(shí)方案2含沙量分布Fig.6 Plan 2 sediment concentration distribution at 0 o′clock on September 19, 2007

3.3 回淤計(jì)算分析

連云港主航道內(nèi)外段的沿程回淤強(qiáng)度如圖12所示。由圖可知，航道最大回淤發(fā)生在航道里程3～5 km，相當(dāng)于旗臺(tái)北防波堤外1～3 km內(nèi)(旗臺(tái)北防波堤距離主航道拐點(diǎn)約2.2 km)，回淤強(qiáng)度約為0.40 m?？傮w來說，主航道外段回淤小于內(nèi)段，航道里程5～40 km測點(diǎn)隨著距離口門越遠(yuǎn)，回淤強(qiáng)度迅速減小，然后趨于平緩。由于受防波堤的阻擋作用，旗臺(tái)港區(qū)內(nèi)部回淤較小，如航道里程0～3 km隨著離口門越遠(yuǎn)回淤越小。

圖12 主航道內(nèi)外段回淤強(qiáng)度Fig.12 Siltation intensity of the main channel

徐圩進(jìn)港航道方案2、方案4沿程回淤強(qiáng)度如圖13所示?？梢姡燠走M(jìn)港航道回淤最大范圍出現(xiàn)在航道里程3～7 km，該段航道2個(gè)方案回淤強(qiáng)度差別不大；其中，最大回淤位于距離徐圩大環(huán)抱口門約5 km處，其風(fēng)后最大回淤強(qiáng)度約為0.30 m。位于口門內(nèi)的0、1 km處回淤比口門小，平均回淤約0.15 m。總體來說，徐圩進(jìn)港航道最大回淤強(qiáng)度比主航道略小，但沿程回淤強(qiáng)度變化較為平緩。

圖13 方案2、4徐圩進(jìn)港航道回淤強(qiáng)度Fig.13 Plan 2、4 siltation intensity of Xuwei entry channel

徐圩一港池航道方案2、方案4沿程回淤強(qiáng)度如圖14所示。其中，方案4與方案2相比，港區(qū)陸域圍填全部完成，港池內(nèi)水域面積最小，一港池航道回淤強(qiáng)度明顯減小，并隨著與口門之間距離的增大而減小，尤其是一港池內(nèi)航道里程0～4 km回淤強(qiáng)度明顯小于方案2。

圖14 方案2、4徐圩一港池航道回淤強(qiáng)度Fig.14 Plan 2、4 siltation intensity of Xuwei port #1

總之，與徐圩進(jìn)港航道回淤相比，徐圩港區(qū)港池內(nèi)回淤強(qiáng)度較小，主要是因?yàn)樵谂_(tái)風(fēng)期間，防波堤掩護(hù)作用較好，港內(nèi)水體含沙量較?。涣硗?，在臺(tái)風(fēng)風(fēng)速較大期間，港池口門處的潮流流向幾乎平行于港池口門，使得直接進(jìn)入港池內(nèi)的泥沙較少。雖然徐圩航道靠近南部的灌河口海域，由于“韋帕”臺(tái)風(fēng)期間，潮流的主流是由西北向東南沿岸流動(dòng)，所以，灌河口海域的高含沙水體不易將泥沙輸移至徐圩航道工程海域，進(jìn)一步說明波浪掀沙、波流共同輸沙的重要性。

4 結(jié)論

1)基于三重網(wǎng)格嵌套的超高分辨率(分辨率分別為9、3、1 km)“韋帕”臺(tái)風(fēng)風(fēng)場，應(yīng)用第3代波浪模型MIKE21-SW、二維潮流模型MIKE21-HD和三維水動(dòng)力泥沙模型FVCOM，以二維、三維大、中、小模型三重嵌套的方式，主要解決了波流耦合問題和三維模型中增加了波生近岸流、波浪對泥沙運(yùn)動(dòng)、泥沙的絮凝沉降和高含沙量時(shí)的制約沉降機(jī)制等影響因素，建立了臺(tái)風(fēng)作用下波浪潮流耦合三維泥沙數(shù)學(xué)模型。

2)應(yīng)用臺(tái)風(fēng)作用下波浪潮流耦合三維泥沙數(shù)學(xué)模型，模擬了2007年“韋帕”臺(tái)風(fēng)作用下連云港海域三維泥沙輸移過程。2#測站的分層含沙量過程的計(jì)算值與實(shí)測值吻合較好?！绊f帕”臺(tái)風(fēng)作用下連云港區(qū)主航道的回淤強(qiáng)度與實(shí)測較為一致，最大實(shí)質(zhì)回淤強(qiáng)度為0.40 m。較好地反映在“韋帕”臺(tái)風(fēng)作用下連云港海域的泥沙運(yùn)動(dòng)特征。

3)計(jì)算分析了方案2、方案4連云港港海域“韋帕”臺(tái)風(fēng)作用下波流耦合三維泥沙運(yùn)動(dòng)特征，得到：連云港主航道航道最大回淤發(fā)生在航道里程3～5 km，相當(dāng)于旗臺(tái)北防波堤外1 km～3km內(nèi)回淤強(qiáng)度約為0.40 m。受防波堤的阻擋作用，旗臺(tái)港區(qū)內(nèi)部回淤較小。徐圩進(jìn)港航道最大回淤值比主航道略小，最大回淤強(qiáng)度約為0.30 m，但其沿程回淤強(qiáng)度變化較為平緩。徐圩港區(qū)港池內(nèi)回淤強(qiáng)度較小，主要是因?yàn)榕_(tái)風(fēng)期間，防波堤掩護(hù)作用較好，港內(nèi)水體含沙量較小。另外，由于“韋帕”臺(tái)風(fēng)期間，風(fēng)暴潮流的主流是由西北向東南沿岸流動(dòng)，所以，灌河口海域的高含沙水體不易將泥沙輸運(yùn)至徐圩航道工程海域，進(jìn)一步說明波浪掀沙、波浪潮流共同輸沙的重要性。