甌江口波流耦合作用下二維流場數(shù)值模擬研究

，

(大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室， 遼寧 大連 116023)

河口近岸海區(qū)[1]，尤其在淺水區(qū)域，波流相互作用對(duì)水動(dòng)力有著重要影響，當(dāng)波浪、水流聯(lián)合作用時(shí)，其流場結(jié)構(gòu)以及剪應(yīng)力特性與純波或純流相比發(fā)生了很大變化[2]，所以研究波流共同作用比研究單純的波浪或潮流更能真實(shí)地描述河口近岸海區(qū)的水動(dòng)力特性。

甌江口海區(qū)島嶼林立、沙灘密布、水道縱橫、地形復(fù)雜，波浪潮流相互作用顯著。近年來有很多關(guān)于該海區(qū)水動(dòng)力環(huán)境的數(shù)值研究[3-4]，其中以李孟國等[5-6]的研究最具代表性。但以上研究一般都是將波浪和潮流分開來處理，考慮波流耦合作用的少；即使是耦合計(jì)算，也是單向耦合的多，雙向耦合的少。因此建立雙向波流耦合作用數(shù)值模型對(duì)該海區(qū)的水動(dòng)力進(jìn)行研究就顯得非常必要。

1 數(shù)學(xué)模型

1.1 波浪模型

MIKE21 SW 是基于波浪作用平衡方程的譜波浪模型，模型可以處理近海和海岸地區(qū)大范圍的波浪的成長、衰減及變形；可以求解準(zhǔn)定常和非定常波浪作用平衡方程的全譜公式和方向解耦參數(shù)公式，并可與2D/3D水流模型完全耦合，是目前有代表性的波浪模型。在笛卡爾坐標(biāo)系下，模型的控制方程為

(1)

式中：σ——波浪相對(duì)頻率；

v——波群速度，v=(cx,cy,cσ,cθ),

其中：cx、cy——波浪在x、y方向的傳播速度；

cσ、cθ——波浪在σ、θ方向的傳播(變形)速度，均采用線性波理論計(jì)算。

(2)

(3)

(4)

(5)

式中：k——波數(shù),k=(kx,ky)；

u——流速,u=(ux,uy)；

s——沿θ方向空間坐標(biāo)；

m——垂直于s的坐標(biāo)；

h——總水深。

方程(1)左邊第一項(xiàng)代表波能隨時(shí)間的變化；第二項(xiàng)代表波能隨地理空間、頻率空間和波向空間的變化；右邊項(xiàng)S=Sin+Snl+Sds+Sbot+Ssurf代表使方程能量守恒的源項(xiàng)，Sin指風(fēng)輸入的能量;Snl指波與波之間的非線性作用引起的能量損耗;Sds指由白帽引起的能量損耗;Sbot指底摩阻引起的能量損耗;Ssurf指由于水深變化引起的波浪破碎產(chǎn)生的能量損耗。

1.2 潮流模型

1.2.1 水流模型

對(duì)于平面大范圍的自由表面流動(dòng)、水深尺度遠(yuǎn)小于平面尺度、無明顯垂直環(huán)流和垂向流速小的淺水流動(dòng)，可用靜水壓力代替動(dòng)水壓力，并沿水深方向進(jìn)行積分來簡化方程，整合水平動(dòng)量方程和連續(xù)方程，得到下面的二維淺水方程。

1)連續(xù)方程。

(6)

2)動(dòng)量方程。

(7)

(8)

式中：t——時(shí)間；

x、y、z——右手Cartesian坐標(biāo)系；

d——靜止水深；

h——總水深，h=d+η；

η——水位；

u、v、w——流速在x、y、z方向上的分量；

ρ——水的密度；

ρo——海水密度；

pa——當(dāng)?shù)氐拇髿鈮海?/p>

f——Coriolis參數(shù)(′Ω是地球自轉(zhuǎn)角速率，Φ為地理緯度)，f=2′sinΦ；

fu,fv——地球自轉(zhuǎn)引起的加速度；

Ax、Ay——應(yīng)力項(xiàng)，由應(yīng)力模型提供；

S——源匯項(xiàng)，(us,vs)源匯項(xiàng)水流流速；

u、v——平均深度x、y方向的速度分量，被定義為

(9)

1.2.2 應(yīng)力模型

應(yīng)力模型用于計(jì)算水流模型中的應(yīng)力項(xiàng)Ax、Ay，是波浪模型和水流模型實(shí)現(xiàn)耦合的關(guān)鍵部分，它包括水平粘滯應(yīng)力、表面風(fēng)應(yīng)力矢量、波流共同作用下的底部切應(yīng)力、波浪輻射應(yīng)力。其方程如下。

(10)

(11)

1)水平粘滯應(yīng)力。Txx、Txy、Tyy為水平粘滯應(yīng)力，計(jì)算采用一種基于平均深度流速梯度的渦粘性公式

式中：A——水平粘滯系數(shù)。

2)表面風(fēng)應(yīng)力矢量。τs為表面風(fēng)應(yīng)力矢量，τs=(τsx,τsy)=ρa(bǔ)cd|uw|uw，ρa(bǔ)為大氣密度;uw=(uw,vw)為海面上10 m處的風(fēng)速矢量；cd為風(fēng)的拖曳系數(shù)

cd=

(12)

3)底部切應(yīng)力。τb為波浪、潮流引起的底部切應(yīng)力[7]，由波浪切應(yīng)力矢量和潮流切應(yīng)力矢量疊加得到

τb=(τbx,τby)=ρocf|ub|ub+лρofw|Uw|Uw/8

(13)

式中：ρo——水的密度；

ub——平均深度的速度,ub=(u,v)；

cf=——水流的底摩阻系數(shù),

cf=h(ksh)1/3/645.16；

ks——河床阻力長度；

T——波周期；

v——水的運(yùn)動(dòng)粘性系數(shù),v=10-6m2/s；

Uw=(Uw,Vw)——近底波浪軌跡速度矢量，

其中：H——波高；

T——波周期；

k——波數(shù)。

4)輻射應(yīng)力。Sxx、Sxy、Syy為輻射應(yīng)力分量[8]，波浪輻射應(yīng)力分量表達(dá)式為

Sxx=E[(2n-1/2)-nsin2θ]

Sxy=Syx=Ensinθcosθ

(14)

Syy=E[(n-1/2)+nsin2θ]

式中：E——波能，E=ρgH2/8；

θ——波向。

n=cg/c

其中：cg——波群速度；

c——波速。

1.3 動(dòng)邊界處理

河口淺灘部位，漲潮時(shí)灘面被逐漸“淹沒”，落潮時(shí)逐漸“干出”。數(shù)值計(jì)算中，選定一標(biāo)準(zhǔn)水深Ho，當(dāng)在某一時(shí)刻某一網(wǎng)格點(diǎn)的實(shí)際水深H≤Ho時(shí)，認(rèn)為該節(jié)點(diǎn)“干出”，令該點(diǎn)的流速值為零。在以后的每個(gè)時(shí)間步長的計(jì)算中，“干出”點(diǎn)潮位值由周圍非“干出”點(diǎn)的潮位值的線性插值得到；當(dāng)在某一時(shí)刻某一網(wǎng)格點(diǎn)的實(shí)際水深H≥Ho時(shí)，則認(rèn)為該節(jié)點(diǎn)被“淹沒”，恢復(fù)程序計(jì)算。

1.4 耦合模型實(shí)現(xiàn)過程

首先在風(fēng)場來臨前一段時(shí)間進(jìn)行潮流的模擬，在得到穩(wěn)定的潮流場后，加上風(fēng)應(yīng)力項(xiàng)，進(jìn)行數(shù)值模擬，得到水位與流速。再將水位和流速輸入到波浪模型，進(jìn)行下一時(shí)間步的波浪場數(shù)值模擬，利用計(jì)算得到的波參數(shù)輸入到應(yīng)力模型，可以計(jì)算出輻射應(yīng)力和底應(yīng)力，然后將計(jì)算得到的輻射應(yīng)力和底應(yīng)力代入水流模型，得到新的水位、流速，然后再把水位、流速代入波浪模型進(jìn)行模擬，見圖1。兩個(gè)模型以此方式循環(huán)交換數(shù)據(jù)實(shí)現(xiàn)雙向耦合，得到考慮波浪和潮流相互作用機(jī)制下的波高、水位和流速等參數(shù)。在這個(gè)耦合模型中，潮流模型的時(shí)間步長和波浪模型相同，使數(shù)據(jù)得到同步交換，保證了計(jì)算精度。

圖1 耦合模型變量傳遞示意

2 模型驗(yàn)證

2.1 模型設(shè)置

以2005年6～7月全潮實(shí)測數(shù)據(jù)為依據(jù)，分別在耦合和未耦合情況下對(duì)甌江口海區(qū)進(jìn)行大、中、小3個(gè)潮過程(從2005年6月14日零時(shí)～7月15日零時(shí))的驗(yàn)證。海洋開邊界以潮位作為驅(qū)動(dòng)力，潮位值取自2005年潮汐表；河口邊界來源于實(shí)測徑流資料，入射波浪條件來源于1994年南麓島波浪觀測站的波浪資料的統(tǒng)計(jì)分析數(shù)據(jù)；風(fēng)場資料來源于ncep/qscat再分析風(fēng)場數(shù)據(jù)，空間分辨率30′×30′，時(shí)間分辨率6 h。

模型的計(jì)算范圍27°46′N～28°24′N,120°45′～121°20′E，采用UTM-51投影將大地坐標(biāo)轉(zhuǎn)化為平面坐標(biāo)，計(jì)算范圍總面積約2 600 km2。網(wǎng)格剖分采用三角形網(wǎng)格，對(duì)復(fù)雜地形附近的計(jì)算網(wǎng)格作進(jìn)一步加密，保證了區(qū)域流場的模擬精度。計(jì)算區(qū)域共布有16 223個(gè)計(jì)算節(jié)點(diǎn)和30 489個(gè)三角形單元，域內(nèi)最大水深達(dá)100 m以上，位于玉環(huán)島正東部。測站模型計(jì)算網(wǎng)格、地形及測站布置見圖2。

圖2 模型計(jì)算網(wǎng)格、地形及測站布置示意

2.2 驗(yàn)證結(jié)果與分析

2.2.1 風(fēng)成浪驗(yàn)證

根據(jù)所得的風(fēng)場資料，分別研究了耦合潮流場及不耦合潮流場情況下的波浪場。表1為波高、周期的驗(yàn)證情況，圖3為某時(shí)刻部分計(jì)算區(qū)域的波高分布圖。

表1 風(fēng)成浪驗(yàn)證結(jié)果

圖3 某時(shí)刻部分計(jì)算區(qū)域波高分布示意(單位:m/s)

計(jì)算結(jié)果顯示，考慮波、流耦合效應(yīng)的波高與周期都比沒有考慮耦合效應(yīng)的情況增加明顯，更為接近觀測結(jié)果。

2.2.2 潮位、流速驗(yàn)證

從小潮型部分測站的潮位、流速情況看，計(jì)算所得的潮位過程和流速過程均與實(shí)測資料吻合良好，并且耦合模型的計(jì)算結(jié)果與實(shí)測資料吻合程度比未耦合的好。圖4為耦合情況下甌江口海區(qū)的漲潮、落潮時(shí)刻流場分布圖。

圖4 耦合模型計(jì)算區(qū)域流場

為了更確切地描述波浪對(duì)流速的影響，表2給出了各測站波流耦合情況下波浪對(duì)流速的影響幅度的平均值與最大值。圖5為落潮憩息時(shí)刻甌江口部分區(qū)域耦合前后的流速等值線分布圖。

表2 波流耦合情況下波浪對(duì)流速的影響幅度表 %

注：最大影響幅度均發(fā)生在落潮憩息時(shí)刻。

圖5 落潮憩息時(shí)刻部分計(jì)算區(qū)域流速等值線分布圖(單位:m)

2.2.3 結(jié)果分析

1)波流耦合作用下甌江口潮流場總體特征。如圖4所示，受甌江口外群島影響，漲、落潮時(shí)水流均被分為幾股水流進(jìn)入、流出甌江口和樂清灣，近岸和江口及樂清灣大部分區(qū)域往復(fù)流明顯，離岸越遠(yuǎn)旋轉(zhuǎn)流越明顯。樂清灣的漲潮水流來自玉環(huán)島—橫趾山水道和大門島—鹿棲島水道，落潮水流沿漲潮路徑流出。甌江口漲落潮水流與樂清灣漲落潮水流不直接交換，南口的漲潮水流主要來自霓嶼島西側(cè)的水道，北口的漲潮水流主要來自大門島兩側(cè)的各個(gè)水道，為水流出海和入海的主要通道。南北口的落潮水流均沿各自沿漲潮路徑流出；南口的漲落潮先于北口的漲落潮。溫州淺灘處于漲潮水流緩流區(qū)，其上下兩側(cè)有比較明顯的出入流，漲潮時(shí)水流從兩側(cè)進(jìn)入，淹沒淺灘；落潮時(shí)水流從兩側(cè)流出，露出淺灘。在小潮型的潮流過程中，由于甌江口外島嶼的阻攔，樂清灣東側(cè)的華歧淺灘淹沒面積較小；而甌江口南部的毆飛淺灘處在開敞區(qū)域，淹沒面積較大。這些現(xiàn)象也和文獻(xiàn)[9]對(duì)甌江口潮流特性的分析相一致。

2)甌江口潮流場對(duì)波浪的作用。由表1可見，具有耦合效應(yīng)的情況改善了波高和周期的計(jì)算結(jié)果，更為接近觀測結(jié)果。尤其是在潮流作用較強(qiáng)的區(qū)域，如中山水道、黃大岱水道等，考慮波流耦合效應(yīng)能夠產(chǎn)生更強(qiáng)的波浪條件，如圖2所示。所以在甌江口圍涂工程的波浪條件計(jì)算中，考慮波、流耦合的波浪計(jì)算尤為重要。

3)波浪對(duì)甌江口潮流場的作用。波浪對(duì)潮位的影響很小，但對(duì)流速的影響較大。波流同向時(shí)，使得結(jié)果增大；波流逆向時(shí)，則使得流速減小。在甌江口北部樂清灣地區(qū)，無論是漲潮時(shí)段還是落潮時(shí)段，波浪一般都使流速增大。在甌江口南口、黃大岙水道、重山水道、黃大岙水道與重山水道之間的區(qū)域及狀元岙深水區(qū),漲潮時(shí)段波浪使流速減小，落潮時(shí)段波浪使流速增大。

從表2可以看到隨著水深、地形的變化，波浪對(duì)流速的影響幅度也不同。由于甌江口外眾多島嶼阻攔了外海波浪的傳播，掩護(hù)了河口及其海域，所以波浪對(duì)測站1JHJ的流速影響幅度較??；測站2JHJ、3JHJ、4JHJ位于岸線較復(fù)雜的區(qū)域，波浪容易發(fā)生變形破碎，流速受波浪影響較大；測站6JHJ位于外部開敞深水區(qū)，波浪對(duì)其流速的影響較小。各個(gè)測站波浪的最大影響幅度均發(fā)生在落潮憩息時(shí)刻，因?yàn)樵诼涑睍r(shí)段末期，潮流流速較小，波浪對(duì)潮流的的作用較大，進(jìn)而使流速發(fā)生較大幅度的改變。由此可見，波浪從深水向淺水傳播的過程中，波浪對(duì)潮流的作用也在逐漸變大，尤其是在近岸區(qū)，形成波浪增水和波生近岸流，疊加到潮流中，從而形成復(fù)雜的近岸流動(dòng)系統(tǒng)，見圖5。因此當(dāng)對(duì)淺水區(qū)進(jìn)行流場模擬研究與計(jì)算時(shí)，應(yīng)考慮波浪對(duì)流的影響。

3 結(jié)論

1)在潮流作用較強(qiáng)的區(qū)域，如甌江北口水道、中水道等，潮流對(duì)波浪有較大的影響。考慮波流耦合效應(yīng)能夠產(chǎn)生更強(qiáng)的波浪條件，因此在這些區(qū)域的工程項(xiàng)目波浪計(jì)算中應(yīng)該考慮潮流對(duì)波浪的影響。

2)由于波周期與潮周期相差甚大，二者可差幾個(gè)數(shù)量階，所以波浪對(duì)潮位的影響甚小，幾乎可以忽略不計(jì)。所以在一般波浪條件下，潮位計(jì)算時(shí)可以不考慮波浪的影響。

3)在近岸淺水區(qū)域和地形變化復(fù)雜的區(qū)域，波浪在傳播過程中更易發(fā)生變形和破碎，形成波浪增水和波生近岸流，疊加到潮流中，進(jìn)而影響流速的大小。但在深水區(qū)，波浪對(duì)流速的影響并不明顯，并且影響幅度很小。因此當(dāng)對(duì)地形復(fù)雜的淺水區(qū)的水動(dòng)力進(jìn)行計(jì)算時(shí)，應(yīng)考慮波浪對(duì)流的影響。

波浪對(duì)潮流的作用主要是通過波流邊界層作用的，故在深水區(qū)，波浪對(duì)潮流的作用在水深方向并不是平均的；而在二維耦合模型中，由于二維平均流的采用降低了波流底應(yīng)力的效應(yīng)，所以波浪對(duì)流速并不明顯，其影響幅度比實(shí)際流速小一個(gè)數(shù)量級(jí)。因此發(fā)展三維波流耦合模型將有助于這一問題的深入研究。

[1] 劉應(yīng)中,時(shí) 鐘.河口波流相互作用的數(shù)學(xué)模型[J].海洋科學(xué),2001,25(5):42-46.

[2] 吳永勝,練繼建,張慶河,等.波浪和水流共同作用下水流時(shí)均流速分布[J].水利學(xué)報(bào),2001,(1):35-41.

[3] 姜恒志,沈永明,汪守東.甌江口甌江口三維潮流和鹽度數(shù)值模擬研究[J].水動(dòng)力學(xué)研究與進(jìn)展,A輯,2009,24(1):63-70.

[4] 陸永軍,李浩麟,董 壯,等.強(qiáng)潮河口圍海工程對(duì)水動(dòng)力環(huán)境的影響[J].海洋工程,2002,22(4):17-25.

[5] 李孟國,王正林.甌江口潮流數(shù)值模擬[J].長江科學(xué)院院報(bào),2002,19(2):19-22.

[6] 李孟國,劉百橋.甌江口海區(qū)波浪數(shù)學(xué)模型研究[J].水道港口,2002,22(1):1-8.

[7] VOULGARIS G,WALLBRIDGE S,TOMLINSON B N,et al.Laboratory investigations into wave period effects on sand bed erodibility,under the combined action of waves and currents[J].Coastal Engineering,1995,26(3):117-134.

[8] LU Y,ZUO L,SHAO X,et al.A 2D mathematical model for sediment transport by waves and tidal currents[J].China Ocean Engineering,2005,19(4):571-586.

[9] 鄭敬云,李孟國,麥 苗,等.甌江口水文泥沙特征分析[J].水道港口,2008,29(1):1-7.