非結(jié)構(gòu)化網(wǎng)格下大范圍波生流模擬和應用

王 平，張寧川

(大連理工大學海岸和近海工程國家重點實驗室，遼寧大連 116024)

波浪在向近岸海域傳播的過程中，由于地形的影響會發(fā)生變形、破碎等現(xiàn)象，伴隨著質(zhì)量、動量和能量的輸運及傳遞，進而形成強度和影響范圍都較大的近岸流。在近岸海域，波生沿岸流的現(xiàn)象非常普遍，尤其在淺水和河口區(qū)域更加顯著，對近岸地貌變遷、泥沙輸運以及近海污染物輸運等都有重要的影響。故在近岸潮流、泥沙和污染物輸運的模擬時考慮波浪的影響必不可少。

目前波生沿岸流數(shù)值模擬方法主要有兩種，一種是基于Boussinesq方程的波生流模型［1］，該模型將波浪和流場的計算融合到一個控制方程中;第二種是建立在Longuet-Higgins和Stewart［2］提出的波浪輻射應力理論的基礎上，該方法是將波浪場計算得到的輻射應力等波浪影響參量作為驅(qū)動因素添加到潮流場的控制方程中，很多學者對該方法進行了大量的數(shù)值模擬和研究，得到了一定的結(jié)果。已有的基于輻射應力理論模擬近岸波生流的方法，較多采用有限差分法對控制方程在矩形網(wǎng)格下進行數(shù)值離散，如包西林［3］、Sun［4］、Yoon［5］和 Tang［6］等，但實際工程中邊界多為不規(guī)則形狀，矩形網(wǎng)格對其擬合的精度較差。Wu［7］、李孟國［8］和唐軍［9］等建立了基于三角形網(wǎng)格的波生流模型，但其波浪場都是基于緩坡方程，計算時網(wǎng)格步長受制于波長的限制，造成計算大范圍波浪場時的效率較低，且不能和流場同步耦合。Choi［10］、朱首賢［11］、王彪［12］等將SWAN模型模擬得到的輻射應力添加到潮流模型中得到一種考慮波浪作用的流場模型，該方法解決了波生流計算時波浪場和潮流場不同尺度的問題，使得波浪和流場可以同步耦合，但受模式的限制在模擬中均未考慮波浪紊動的影響，同時在非結(jié)構(gòu)化網(wǎng)格中SWAN模型并未考慮繞射的影響。

從含流的緩坡方程［13］出發(fā)，推出光程函數(shù)方程，聯(lián)立光程函數(shù)方程和波作用守恒方程建立了基于非結(jié)構(gòu)化網(wǎng)格下考慮繞射效應的波浪模型，將波浪場模擬得到的輻射應力和紊動系數(shù)等同步添加到流場模型(FVCOM［14］)中得到基于非結(jié)構(gòu)化網(wǎng)格下的大范圍波生流數(shù)值模式。其中波浪場模型的計算不受波長的限制，適用于近岸大范圍波生流計算，同時使得波浪場和流場可以同步耦合計算。對所建立的模型進行了多個算例驗證，數(shù)值模式的結(jié)果與實測吻合較好，驗證了模型精度、可靠性以及對復雜邊界的適應性。同時將模型應用于大連近岸區(qū)域，對常見浪在大連近岸區(qū)域產(chǎn)生的波生流進行了研究，驗證了模式對大范圍區(qū)域波生流計算的適用性。

1 數(shù)值模式

1.1 近岸波浪場控制方程

Kirby［12］推導出的含流非定常波浪方程，其表達式為:

其中方程左端第三項和第四項為波能密度在θ方向和σ方向的傳播，右端為源匯項，可以包括波浪破碎，底摩阻等對能量的影響。利用波數(shù)矢K替代波數(shù)k的求解，可以計算繞射對波高的影響。改進后波速、波群速以及波能在θ和σ方向的傳播速度分別為:

聯(lián)立式(2)、式(4)建立包含繞射的波浪計算模式，該模式在空間網(wǎng)格上的劃分不再受波長的限制，可以用于大范圍波浪傳播計算;同時利用有限體積法對控制方程進行離散，離散基于非結(jié)構(gòu)化網(wǎng)格可以和目前比較通用的流場模型同步耦合，且對復雜的岸線能很好的擬合;而且光程函數(shù)方程可以有效改善模型在計算近岸復雜海域的波浪繞射效應［16］，提高了模型計算近岸波浪傳播的精度。

1.2 近岸波生流控制方程

流場模型采用修改后的三維水動力模型FVCOM，其水動力方程采用在ζ坐標系下的三維Navier-Stokes方程組，同時引入三維波生時均剩余動量和波浪的紊動摻混效應對水體的影響，方程基本形式見式(5)、式(7):

式中:η為自由水面;x和y為水平坐標;ζ=(z-η)/D為轉(zhuǎn)換后的垂向坐標;w為ζ坐標系下的垂向流速分量;u和v為水平流速分量;Kmc為波流共存的垂向紊動粘性系數(shù);Amc為波流共存的水平紊動粘性系數(shù);R為波浪輻射應力項。

對波流共同作用下的紊動系數(shù)采用文獻［17］的方法，即分別單獨求解水流與波浪引起的紊動系數(shù)A和K，并將其線性疊加，可表述為:

波浪引起的紊動系數(shù)Aw和Kw計算均采用Larson-Kraus公式計算，Aw=Kw=λumaxH=2λH2cosh(1+σ)kD)/Tsinh(kD)，式中umax為波浪底部質(zhì)點最大流速;λ為無因次系數(shù)。

水流引起的水平紊動系數(shù)采用Smagorinsky公式(10)。式中C為常數(shù);Ω為控制單元體面積。水流引起的垂向紊動系數(shù)Km采用Mellor-Yamada紊流閉合模型求解。

FVCOM模型采用外模和內(nèi)模分別計算的方法，其中外模為沿深度平均的二維流場計算，內(nèi)模為沿深度分層的三維流場計算。計算波浪對流場的影響時，在外模和內(nèi)模均考慮輻射應力項。波浪輻射應力分項采用波能E等參數(shù)計算，二維輻射應力形式具體計算公式如下:

三維輻射應力采用鄭金海［18］推導的公式:

上述公式將輻射應力的計算沿深度分為三個計算區(qū)域，R01、R02和R03的具體計算形式見文獻［18］。

流場和波浪的耦合過程分為，流場和波浪單獨計算以及參量相互傳遞兩個過程。其中流場為波浪提供流速和水位參量，而波浪場則為流場提供輻射應力項和紊動系數(shù)。流場和波浪場在完成每一步計算后以及下一步計算開始之前進行參量傳遞，最終實現(xiàn)波流的相互耦合計算。

1.3 控制方程的離散

基于三角形網(wǎng)格，波浪場采用的有限體積法對波作用守恒方程進行離散;對光程函數(shù)方程采用網(wǎng)格中心格式的有限體積法離散。時間離散采用歐拉向前格式，空間采用格林公式，具體形式如下:

對式(2)采用網(wǎng)格頂點格式的有限體積法在控制體內(nèi)積分，空間求導采用將面積分轉(zhuǎn)化為線積分，并將線積分寫成各控制邊求和的形式，整理得到離散后的方程為:

式中:MT(i)為第i個控制體的邊數(shù)，對三角形控制體取為3;Ai為控制體面積;標注m為控制體第m邊上的參量;δ表示控制體邊上參量在x或y方向上的梯度，采用上述格林公式計算;ζ表示參量在時間上的梯度，采用向后差分的一階格式離散。

對式(3)分四步離散求解，空間離散采用單元中心頂點格式的有限體積法在控制體內(nèi)積分;對波譜在頻率方向上的離散采用FCT［19］離散方法;對波譜在傳播方程上的離散采用二階隱身Crank-Nicolson差分方法;對源項對波譜的影響采用二階隱式中心差分離散求解。離散后見式(15a)～(15d)。

式中:NT(i)為第i個控制體的邊數(shù);Ai為控制體面積;標注m為控制體第m邊上的參量;FCT中的φ1和A*計算公式見文獻［19］;Cgx，i，m和 Cgy，i，m為控制體邊 m 上 U+(Cg/k)K 在 x，y 方向上的分量。

流場控制方程的離散見文獻［14］，對于增加的輻射應力項同樣采用單元中心頂點格式的有限體積法在控制體內(nèi)積分，同時用格林公式將面積分轉(zhuǎn)換為各控制邊的求和形式，具體形式見式(16a)～(16b)。

2 數(shù)值模型驗證及應用

2.1 沿岸流計算及實驗驗證

波浪斜向入射近岸時，破碎引起的輻射應力會引起水體的順岸移動而發(fā)生沿岸流。王淑平［20］進行了沿岸流的模擬實驗，實驗分別研究了1∶40和1∶100兩種坡度下波浪斜向入射產(chǎn)生的沿岸流，實驗工況見表1。

表1 沿岸流實驗計算工況參數(shù)Tab.1 The calculation parameters of longshore current experiment

其中D為波浪入射處水深，H0為入射波高，T為入射波周期，γ為波浪破碎指標，Bfric為波流作用下的底摩擦系數(shù)。波浪沿正向入射，圖1給出了兩種工況計算穩(wěn)定后的波生流平均流場圖，圖2和圖3給出了2種計算工況下的波高、增減水和沿岸流平均流速的數(shù)值解和實測值的比較。從計算結(jié)果可以看出，所建立的模式可以很好地模擬波浪斜向入射近岸時引起的沿岸流現(xiàn)象。

圖1 兩種工況下波生沿岸流流場Fig.1 The wave-induced longshore current velocity field of two cases

圖2 工況1數(shù)值解與實測值對比Fig.2 Comparison between the numerical and experimental results in case 1

圖3 工況2數(shù)值解與實測值對比Fig.3 Comparison between the numerical and experimental results in case 2

在破碎帶內(nèi)外選取4個點得到其流速的垂線分布，見圖4，Dx為點距岸邊的垂線距離。從圖中可以看出波生流的垂線分布區(qū)別于潮流的近對數(shù)分布形態(tài)，而垂向分布較為均勻。根據(jù)Visser［21］的研究，這種特征與波浪引起的附加垂向摻混有關，其中較大的附加粘性可使流速垂向剖面更加均勻。

圖4 兩種工況下四個點的流速垂向分布Fig.4 The vertical velocity distribution of four points in two cases

2.2 人工島后波生環(huán)流的計算及驗證

波浪傳播過程中遇到人工島等障礙物時，會在人工島后發(fā)生繞射;由于波高的變化從而產(chǎn)生輻射應力，在輻射應力的驅(qū)動下會形成島后環(huán)流現(xiàn)象。日本東京電力研究所于1994年利用大型水池實驗研究了柏崎刈羽核電站人工島內(nèi)側(cè)取水口周圍波浪場以及波浪破碎而引起的波生流［22］。

試驗中人工島布置在坡度為1∶50的斜坡上，整個實驗區(qū)域長為23.5 m，寬為19.2 m;人工島大小為6.6 m×3.3 m。模型計算域選取19.2 m×15.0 m，入射波高為0.03 m，周期為1.0 s，波浪破碎指標選為0.76，流場底摩擦系數(shù)選為0.025。

圖5給出實驗得到的波生流實測流場分布，表2給出了A～J10個點的平均流場實測值和計算值對比，圖6給出了波生流模型的計算網(wǎng)格和平均流場分布計算結(jié)果。根據(jù)對比可以看出本模型流場流速的數(shù)值解和實測值比較接近;模擬得到在人工島內(nèi)側(cè)的對稱軸附近產(chǎn)生了回流現(xiàn)象，同時在島內(nèi)側(cè)區(qū)域形成一對方向相反的環(huán)流現(xiàn)象，該形態(tài)和實測得到人工島后的波生流結(jié)果一致。

圖5 人工島后波生流實測流場Fig.5 The measured result of wave-induced current velocity field behind the artificial island

表2 人工島后波生流場沿深度平均流速的計算值和實測值對比Tab.2 Comparison between the numerical and experimental results of the average velocity along the depth

圖6 人工島后波生流計算網(wǎng)格和流場計算結(jié)果Fig.6 Computational grids and wave-induced current velocity field

2.3 連續(xù)凹槽地形上的波生流計算及驗證

Dingemans［23］和李孟國［8］等在驗證近岸波生流模型時，均采用波浪在正向和斜向入射連續(xù)凹槽地形來驗證模式的適定性。連續(xù)凹槽海區(qū)的水深可以用以下表達式表示:

模型的計算水深圖和計算網(wǎng)格見圖7，分別計算波浪正向入射(T=4.0 s，H0=0.92 m，θ0=90°)和斜向入射(T=4.0 s，H0=1.0 m，θ0=63°)兩種工況。李孟國模擬得到的波生流場見圖8。

圖7 模型計算網(wǎng)格和水深Fig.7 Computational grids and the bottom topography

本模式模擬得到的平均波生流場見圖8，得到波浪在正向入射連續(xù)凹槽地形時，會在凹槽的兩邊產(chǎn)生相向的沿岸流，沿岸流匯聚到凹槽頂點時會形成水體的雍高，同時雍高的水體則會通過破碎帶回流海洋從而產(chǎn)生裂流;波浪在斜向入射連續(xù)凹槽時則會產(chǎn)生順凹槽形狀的沿岸流。對比圖8和圖9可知本模式模擬的結(jié)果和已有研究結(jié)果較吻合。

2.4 大連灣區(qū)域的波流耦合計算

由于本波浪模式的空間步長不再受制于波長的限制，波浪模型計算與流場計算可以共用一套網(wǎng)格，波浪和流場模型同步計算，每步計算完成后將流場實時得到的水深、流速提供給波浪場;同時將波浪場得到的輻射應力、紊動系數(shù)等反饋給流場，即實現(xiàn)了波浪場和流場之間參量的雙向傳遞和同步耦合。

模型計算選擇大連灣及附近海域作為研究對象，其中波浪采用不規(guī)則波，不規(guī)則波的能量采用JONSWAP譜計算。圖10(a)和10(b)分別為計算區(qū)域地形和計算網(wǎng)格圖。波浪采用S和SE向入射(為大連灣區(qū)域的常見浪向)，波高為2.0和4.0 m，周期為8 s;計算中波浪和流場的步長均選為5 s，最小空間步長為100 m。

模型分兩步計算完成。首先計算波浪場，等波浪場計算穩(wěn)定后，流場開始計算，并將波浪場計算所得的輻射應力、紊動系數(shù)等參量添加進流場;此時波浪和流場同步計算直至波生流場達到穩(wěn)定形態(tài)。由于波生流多發(fā)生在近岸區(qū)域，模型選取圖10(a)中方框中的波生流場分析。

圖8 李孟國模擬得到的波生流場Fig.8 The wave-induced current velocity field simulated by LI Meng-guo

圖9 本模式模擬得到的波生流場Fig.9 The wave-induced current velocity field computed by the present model

圖10 大連灣海域分布及模型計算網(wǎng)格Fig.10 The area distribution and computational mesh of the Dalian Bay

圖11(a)和11(b)為波浪S向入射下波高分別為2.0 m和4.0 m產(chǎn)生的波生流平均流場;圖12(a)和12(b)為波浪SE向入射下波高分別為2.0 m和4.0 m產(chǎn)生的波生流場。從圖可知波浪斜向入射破碎時會在波浪破碎處至岸邊區(qū)域發(fā)生沿岸流，沿岸流的方向同波浪入射成小角度;當遇到凹槽地形時，若產(chǎn)生兩股相對的沿岸流時(圖11)，則會在相匯處形成逆流;波高越大，波浪破碎處離岸越遠，沿岸流的影響范圍就會越大，強度也會越大;同時沿岸流的強度和近岸潮流強度為同一數(shù)量級，在考慮近岸污染物、泥沙輸運時，波生流為主要影響因素之一。

圖11 S向波浪入射在近岸產(chǎn)生的沿岸流Fig.11 The wave-induced current velocity field by the south incident wave

圖12 SE向波浪入射在近岸產(chǎn)生的沿岸流Fig.12 The wave-induced current velocity field by the south-east incident wave

圖13為4.0 m有效波高的波浪S向和SE向入射后，形成的近岸增減水圖。從圖中可知波浪在近岸破碎時不僅會形成沿岸的波生流，同時在波浪破碎區(qū)形成增水，而在波浪未破區(qū)域形成減水。

從上述結(jié)果可以得到，波浪對流場的影響主要分為兩個方面，一是會在近岸破碎帶內(nèi)產(chǎn)生沿岸流，二是會在近岸形成增減水;而流對波浪的影響主要表現(xiàn)為，波流同向時波高減小、周期變大，而波流逆向時波高增大、周期變小，當波浪與水流斜向相交時，波浪則會發(fā)生折射現(xiàn)象。

圖13 S向和SE向波浪入射后在近岸形成的增減水Fig.13 The distribution of wave set-up and set-down by the south and south-east incident waves

3 結(jié)語

在非結(jié)構(gòu)化網(wǎng)格下基于光程函數(shù)方程和波作用守恒方程，建立了考慮繞射效應的大范圍波浪計算模型;利用波浪模型計算得到的輻射應力、波浪紊動系數(shù)等參量以及三維水動力模型(FVCOM)建立了非結(jié)構(gòu)化網(wǎng)格下的大范圍波生流數(shù)值模型。模型分別對波浪斜向入射產(chǎn)生的沿岸流、人工島后的波生環(huán)流以及連續(xù)凸起地形下的逆流進行了模擬，數(shù)值結(jié)果和實測結(jié)果對比表明，本波生流數(shù)值模式可以高效精確地實現(xiàn)近岸波浪破碎產(chǎn)生的近岸流的模擬計算。

利用本模型對實際海域的計算結(jié)果，表明該模型可以用于大范圍的波流同步耦合計算;同時非結(jié)構(gòu)化網(wǎng)格可以很好地擬合復雜岸線的變化;且說明波生流為近岸水動力不可缺少的影響因素之一。由于缺少實際海域波生流的實測資料，本模式對實際海域波生流的模擬計算需要進一步的驗證。

［1］ 盧 吉，余錫平.基于Boussinesq方程的近岸波流統(tǒng)一模型［J］.水動力學研究與進展:A輯，2008，23(3):314-320.

［2］ Longuet-Higgins M S，Stewart R W.Radiation stress and mass transport in gravity waves with application to surf beats［J］.Journal of Marine Research，1962，13:481-504.

［3］ 包西林，西村仁嗣.Hardy-Cross法波生流數(shù)值解析［M］.北京:海洋出版社，2006.

［4］ Sun T，Tao J H.Experimental and numerical study of wave-induced long-shore currents on a mild slope beach［J］.China Ocean Engineering，2005，19(3):469-484.

［5］ Yoon S B，Cho Y，Lee C.Effects of breaking-induced currents on refraction-diffraction of irregular waves over submerged shoal［J］.Ocean Engineering，2004，31(5-6):633-652.

［6］ Tang J，Shen Y M，Qiu D H.Numerical study of pollutant movement in waves and wave-induced long-shore currents in surf zone［J］.Acta Oceanologica Sinica，2008，27(1):122-131.

［7］ Wu C S，Liu P L F.Finite element modeling of nonlinear coastal currents［J］.Journal of Waterway，Port，Coastal，and Ocean Engineering，1985，111(2):417-432.

［8］ 李孟國.波生近岸流的數(shù)學模型研究［J］.水道港口，2003，24(4):161-166.

［9］ 唐軍，魏美芳.非結(jié)構(gòu)化網(wǎng)格下近岸波生流數(shù)值模擬［J］.海洋學報，2010，32(6):41-46.

［10］ Choi J，Lim C，Lee J，et al.Evolution of waves and currents over a submerged laboratory shoal［J］.Coastal Engineering，2009，56(3):297-312.

［11］朱首賢.流、浪模式和物質(zhì)長期輸運分離研究［D］.上海:華東師范大學，2005.

［12］王 彪，沈永明，王 亮.長興島海區(qū)波流相互作用數(shù)值模擬研究［J］.海洋工程，2012，30(3):87-96.

［13］ Kirby J T.A note on linear surface wave-current interaction over slowly varying topography［J］.J.Geophys Res.，1984，89:745-747.

［14］Chen C，Liu H，Beardsley R C.An unstructured grid，finite-volume，three-dimensional，primitive equations ocean model:Application to coastal ocean and estuaries［J］.J.Atm.＆ Oceanic Tech.，2003，20:159-186.

［15］ Mei C C.The Applied Dynamics of Ocean Surface Waves［M］.New York:Word Scientific Co.Pte.Ltd.，1983:740.

［16］ Hong Guang-wen.Mathematical models for combined refraction-diffraction of waves on non-uniform current and depth［J］.China Ocean Eng.，1996，10(4):433-454.

［17］ Xia H Y，Xia Z W，Zhu L S.Vertical variation in radiation stress and wave-induced current［J］.Coast.Eng.，2004，51(4):309-321.

［18］鄭金海，嚴以新，彭世根.波浪剩余動量流垂向分布研究［J］.河海大學學報，2000，28(1):8-13.

［19］ Boris J P，Book D L.Flux corrected transportⅠ，SHASTA，a fluid transport algorithm that works［J］.J.Comp.Phys. ，1973，11:38-39.

［20］王淑平.沿岸流的研究［D］.大連:大連理工大學，2001.

［21］ Visser P J.Laboratory measurements of uniform longshore currents［J］.Coastal Engineering，1991，15:563-593.

［22］李紹武，柴山知也.近岸區(qū)人工島周圍波生流系統(tǒng)數(shù)學模型的驗證［J］.天津大學學報，1998，31(5):584-589.

［23］ Dingemans M W，Radder A C，De Vriend H J.Computation of the driving forces of wave-induced currents［J］.Coastal Engineering，1987，11:539-563.