趙張益，張慶河

(天津大學(xué) 水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072)

波浪及其破碎所產(chǎn)生的波生流對(duì)近岸區(qū)域污染物和泥沙輸運(yùn)等具有顯著影響。掌握近岸海域波生流場(chǎng)的分布規(guī)律，對(duì)海岸工程規(guī)劃設(shè)計(jì)以及海岸帶污染治理具有重要意義［1-2］。由于波生流垂向結(jié)構(gòu)的非均勻性，發(fā)展三維波生流或近岸波流相互作用數(shù)值模型是近年來(lái)的主要趨勢(shì)。

目前，學(xué)者們通常采用比較常用的數(shù)值格式如有限差分(FD)法，有限元(FE)法，有限體積(FV)法等［3-5］來(lái)建立和發(fā)展三維波生流模型，其中，Warner等［3］在ROMS模型基礎(chǔ)上發(fā)展了三維波生流模型，但并未對(duì)其進(jìn)行充分驗(yàn)證，同時(shí)也沒(méi)有考慮波浪附加紊動(dòng)效應(yīng);Wang等［4］在FVCOM模型基礎(chǔ)上建立了三維波生流模型，并對(duì)模型進(jìn)行了驗(yàn)證，不足在于模型沒(méi)有考慮波浪附加紊動(dòng)效應(yīng);Xie［5］在ECOM模型基礎(chǔ)上建立了三維波生流模型，其優(yōu)點(diǎn)在于考慮了波浪附加紊動(dòng)效應(yīng)，并對(duì)模型進(jìn)行了大量的驗(yàn)證。近年來(lái)，間斷有限元(DG)法因?yàn)槠湓谑睾阈浴⒎€(wěn)定性、處理復(fù)雜地形、網(wǎng)格自適應(yīng)以及并行計(jì)算等方面所具有的優(yōu)勢(shì)［6］，在三維水動(dòng)力模型的數(shù)值求解中越來(lái)越得到重視［7-10］。然而，采用DG法的近岸波生流模型目前尚未見(jiàn)諸報(bào)道。為了給近岸波生流模擬提供高效準(zhǔn)確的模擬工具，將采用DG法構(gòu)建一個(gè)非結(jié)構(gòu)化網(wǎng)格的三維波生流數(shù)值模型，并將其應(yīng)用于三維近岸波生流模擬中以說(shuō)明模型的合理性。

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

1.1 三維波生流控制方程

三維波生流控制方程包含以下動(dòng)量方程和連續(xù)方程:

式中:Am為水平渦粘系數(shù)，Kmv為垂向渦粘系數(shù)，Mij為三維輻射應(yīng)力，Rij為波浪水滾剪切應(yīng)力。

三維輻射應(yīng)力采用如下形式的公式［11］:

根據(jù)Warner等［3］的工作，波浪水滾剪切應(yīng)力公式:

垂向渦粘系數(shù)采用k-ε模型［12］計(jì)算，水平渦粘系數(shù)同時(shí)考慮水流和波浪破碎影響。采用Smagorinsky渦粘模型［13］計(jì)算水流渦粘系數(shù)

采用如下形式［14］計(jì)算波浪破碎引起的附加水平紊動(dòng)系數(shù)

考慮波流共同影響的水平渦粘系數(shù):

式中:α為水流水平紊動(dòng)渦粘系數(shù)，在模型中取0.1;Λ為波浪水平摻混程度的經(jīng)驗(yàn)系數(shù)，將在后面的算例中給出其取值;Umax為近底波浪水質(zhì)點(diǎn)軌跡速度水平分量的幅值;H為波高。

1.2 邊界條件

根據(jù)具體情況，開邊界采用Dirichlet邊界條件:

或Neumann邊界條件:

式中:U0為給定流速矢量;n為邊界法向量。

岸線或建筑物等固邊界采用流速法向分量為0，即Un=0。

自由表面邊界條件由表面的Reynolds應(yīng)力和風(fēng)應(yīng)力平衡可得:

式中:風(fēng)應(yīng)力可以表示為:

式中:ρa(bǔ)為空氣密度;W為水面以上垂直高度10 m處風(fēng)速;Wx和Wy分別為x和y方向上的風(fēng)速分量;CDs為風(fēng)應(yīng)力系數(shù)，采用如下形式［15］:

底邊界條件由底部Reynolds應(yīng)力和底部切應(yīng)力平衡可得:

式中:τbx和 τby為底部切應(yīng)力，采用如下形式［16］:

式中:cf為波流共同作用下的底摩阻系數(shù)，將在后面的算例中給出其取值。

1.3 波浪模型

式中:N為動(dòng)譜密度;σ為波浪相對(duì)頻率;θ為波向;Cx和Cy分別為x和y方向的波浪傳播速度;Cσ和Cθ分別為σ和θ空間的波浪傳播速度。

2 數(shù)值離散

2.1 空間離散

為了離散方便起見(jiàn)，將式(1)寫為如下形式:

式中:fi(U)為非粘性通量;fv(2U)為粘性通量;S包含壓力、柯氏力、輻射應(yīng)力和波浪水滾剪切應(yīng)力等。

引入輔佐變量m，式(19)可以寫為:將計(jì)算域Ω劃分為不重合的N個(gè)單元的集合，e為集合中任一單元，?e為單元邊界。取任意光滑試驗(yàn)函數(shù)φ與式(20)相乘，在單元e上積分，并在計(jì)算域Ω上求和，得到

對(duì)式(21a)做分部積分，得到:

對(duì)式(22)再次做分部積分，則有:

由于DG格式允許單元交界處存在間斷，導(dǎo)致單元內(nèi)部流速梯度與流速近似解的梯度不等價(jià)，為此，引入局部提升算子 r?e和總體提升算子 r［18］。

推進(jìn)最嚴(yán)格的水資源管理制度落實(shí)。抓好實(shí)行最嚴(yán)格水資源管理制度有關(guān)規(guī)定的落實(shí)和考核工作，開展中水回用項(xiàng)目建設(shè)規(guī)劃、水資源保護(hù)規(guī)劃和水中長(zhǎng)期供求規(guī)劃工作，完成地下水的調(diào)查評(píng)估，劃定禁采區(qū)和限采區(qū)，并向社會(huì)公告。

式中:s為單元內(nèi)部邊界的集合。

由此可以求得:

對(duì)式(21b)做分部積分，并代入m，可得:

連續(xù)方程的離散可以類似處理，此處不再贅述。

2.2 時(shí)間離散

動(dòng)量方程通過(guò)上述空間離散，可以得到半離散形式的矩陣方程:

式中:M為質(zhì)量矩陣;I為非粘性項(xiàng)矩陣;V為粘性項(xiàng)矩陣;R包含邊界，壓力，柯氏力，輻射應(yīng)力和波浪水滾剪切應(yīng)力項(xiàng)等。

采用θ格式［24］對(duì)式(28)進(jìn)行時(shí)間離散，可以得到:

式中:θ在［0，1］之間取值。

通過(guò)上述時(shí)間推進(jìn)，可以得到近岸波生流的數(shù)值解。

3 模型驗(yàn)證及結(jié)果分析

在FLUIDITY純水動(dòng)力模型［7］基礎(chǔ)上，添加波浪影響項(xiàng)，建立了三維近岸波生流數(shù)值模型，關(guān)于原始純水動(dòng)力模型的驗(yàn)證可以參閱文獻(xiàn)［25］，這里主要對(duì)所建立的近岸波生流模型進(jìn)行驗(yàn)證。

3.1 近底回流

Ting和Kirby［26］在如圖1所示水槽中進(jìn)行了近底回流實(shí)驗(yàn)。試驗(yàn)水槽尺寸為40 m×0.6 m×1.0 m(長(zhǎng)×寬×高)，海岸坡度為1∶35，平底區(qū)域靜水深0.4 m。規(guī)則波向岸正向入射，波高0.125 m，周期2.0 s。

采用SWAN模型模擬波浪傳播，由于在破波區(qū)域波高模擬結(jié)果偏小，所以在該區(qū)域采用擬合實(shí)測(cè)數(shù)據(jù)的方式對(duì)波高進(jìn)行確定［27］，擬合曲線如圖2所示。

圖1 試驗(yàn)布置示意Fig.1 Sketch of experiment

數(shù)值模型中，平面計(jì)算域40 m×0.6 m，布置三角形單元5 678個(gè)，節(jié)點(diǎn)3 244個(gè);垂向每層高0.01 m。由于近底回流算例對(duì)底部切應(yīng)力比較敏感，故在模型中采取分段設(shè)置的方式，破波點(diǎn)以外cf取0.005，破波點(diǎn)以內(nèi)cf在0.005～0.001之間取值，且隨水深減小線性遞減。Λ取0.01。造波端設(shè)置為開邊界，坐標(biāo)原點(diǎn)位于造波機(jī)處。

波浪增減水驗(yàn)證如圖3所示?？梢钥闯隹紤]波浪水滾剪切應(yīng)力的計(jì)算值與實(shí)測(cè)值吻合較好，反映了波浪增減水過(guò)程，而不考慮波浪水滾剪切應(yīng)力得到的最大減水點(diǎn)位置向破波帶外延伸，且增水值偏大。流速垂向分布驗(yàn)證如圖4所示，其中測(cè)點(diǎn)坐標(biāo)如表1所示。從流速垂向分布的計(jì)算結(jié)果可以看出，其數(shù)量級(jí)、垂向分布形態(tài)與實(shí)測(cè)值較為一致，但中間水體的計(jì)算流速梯度分布則相對(duì)比較均勻，其原因可能是由于數(shù)學(xué)模型中采用的三維輻射應(yīng)力表達(dá)式是基于線性波理論，而實(shí)際情況中，由于較大波浪及破碎引起的強(qiáng)非線性效應(yīng)使得二者之間得到的流速垂向結(jié)構(gòu)有一定區(qū)別。模型整體上反映了近底回流的實(shí)際現(xiàn)象。

表1 流速測(cè)點(diǎn)坐標(biāo)Tab.1 Coordinates of velocity stations

圖2 波高擬合結(jié)果Fig.2 Comparison between the fitted and observed wave height

圖3 波浪增減水驗(yàn)證Fig.3 Comparison between the simulated and observed wave setup

岸灘剖面不同位置處的流場(chǎng)垂向分布如圖5所示。由圖可知:在破波帶內(nèi)表層附近水體沿向岸方向流動(dòng)，底層水體沿離岸方向流動(dòng);在破波點(diǎn)外表層水體沿離岸方向流動(dòng)，而底層水體沿向岸方向流動(dòng)。計(jì)算得到的流場(chǎng)特征與已有文獻(xiàn)［28］對(duì)近底回流現(xiàn)象的描述一致。

圖4 流速驗(yàn)證Fig.4 Comparisons between the simulated and observed velocities

圖5 流速垂向分布Fig.5 Flow field structure of the undertow

3.2 沿岸流

Visser［29］在Delft理工大學(xué)進(jìn)行了一系列沿岸流試驗(yàn)，模型平面布置如圖6所示。試驗(yàn)水槽尺寸為34 m×16.6 m×0.68 m(長(zhǎng)×寬×高)。這里選取其中一組試驗(yàn)進(jìn)行模型的驗(yàn)證，海岸坡度為1∶20，平底區(qū)域靜水深0.35 m，波浪入射方向與平直海岸的法向夾角為15.4°，入射波高0.078 m，周期1.02 s。同樣對(duì)破波區(qū)波高進(jìn)行擬合，如圖7所示。

數(shù)值模型中，平面計(jì)算域34 m×16.6 m，布置三角形單元8 010個(gè)，節(jié)點(diǎn)4 131個(gè);垂向每層高0.035 m。Λ取0.1，cf取0.018。造波端設(shè)置為開邊界，坐標(biāo)原點(diǎn)位于岸線處。

圖6 試驗(yàn)布置示意Fig.6 Sketch of experiment

圖7 波高擬合結(jié)果Fig.7 Comparison between the fitted and observed wave height

圖8 水深平均沿岸流速驗(yàn)證Fig.8 Comparison between the simulated and observed longshore current speeds

沿岸流沿水深平均流速驗(yàn)證結(jié)果如圖8所示，其中實(shí)測(cè)流速的測(cè)點(diǎn)坐標(biāo)如表2所示，考慮波浪水滾剪切應(yīng)力的沿岸流速計(jì)算值與實(shí)測(cè)值吻合程度較好，而不考慮波浪水滾剪切應(yīng)力的沿岸流速峰值偏小，且峰值位置向離岸方向推移。沿岸流速垂向分布驗(yàn)證如圖9所示，其中x軸坐標(biāo)如表3所示，從圖中可以看出，個(gè)別測(cè)點(diǎn)處表層附近流速有一定差異，原因可能是由于模型沒(méi)有考慮波浪運(yùn)動(dòng)引起的附加垂向紊動(dòng)效應(yīng)。但總的來(lái)說(shuō)，模擬得到的沿岸流速垂向分布與實(shí)測(cè)值較為相近，模型能夠合理地模擬波浪斜向入射而形成的沿岸流現(xiàn)象。

表2 水深平均流速測(cè)點(diǎn)坐標(biāo)Tab.2 Coordinates of depth-averaged velocity stations

表3 流速垂向分布測(cè)點(diǎn)坐標(biāo)Tab.3 Coordinates of velocity stations

圖9 沿岸流垂向分布驗(yàn)證Fig.9 Comparisons between the simulated and observed vertical profiles of the longshore current

4 結(jié)語(yǔ)

基于DG方法構(gòu)建了非結(jié)構(gòu)化網(wǎng)格三維波生流數(shù)值模型，并與近底回流和沿岸流物模實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比，主要結(jié)論如下:

1)在三維水動(dòng)力模型中引入三維輻射應(yīng)力、波浪水滾剪切應(yīng)力和波浪附加紊動(dòng)系數(shù)的概念，并首次基于DG方法建立了三維波生流數(shù)值模型。

2)近底回流算例的模擬結(jié)果顯示，所建立的數(shù)值模型能夠合理反映波生流的垂向結(jié)構(gòu)，模擬得到的破波帶附近流速的垂向變化規(guī)律與已有文獻(xiàn)對(duì)近底回流現(xiàn)象的描述一致，數(shù)值解與實(shí)測(cè)值基本吻合。

3)沿岸流算例的模擬結(jié)果表明，模型能夠合理模擬由于波浪斜向入射而形成的沿岸流現(xiàn)象，沿岸流速計(jì)算值與實(shí)測(cè)值吻合程度良好。

4)在波生流數(shù)值模擬中，不能忽略由于波浪破碎所引起的表面水滾的影響。針對(duì)近底回流算例，波浪水滾剪切應(yīng)力使得最大減水點(diǎn)位置向破波帶內(nèi)延伸;針對(duì)沿岸流算例，波浪水滾剪切應(yīng)力使得沿岸流峰值位置向岸推移，且峰值流速增大。

5)模擬結(jié)果表明，算法的構(gòu)建是比較合理的，但在模型中需要進(jìn)一步引入波浪運(yùn)動(dòng)引起的附加垂向紊動(dòng)效應(yīng)，使其能夠更好地描述近岸波生流場(chǎng)。另外，模型中有關(guān)參數(shù)如波浪紊動(dòng)摻混強(qiáng)度的系數(shù)、底部摩擦系數(shù)等如何取值，還有待進(jìn)一步研究。

［1］ 辛文杰.潮流、波浪綜合作用下河口二維懸沙數(shù)學(xué)模型［J］.海洋工程，1997，151(1):30-47.(XIN Wen-jie.Numerical model of 2D estuarial suspended sediment motion under the interaction of tidal flow and waves［J］.The Ocean Engineering，1997，151(1):30-47.(in Chinese))

［2］ 鄭金海，Hajime Mase.波流共存場(chǎng)中多向隨機(jī)波浪傳播變形數(shù)學(xué)模型［J］.水科學(xué)進(jìn)展，2008，19(1):78-83.(ZHENG Jin-hai，Hajime Mase.Numerical modelling of multi-directional random wave transformation in wave-current coexisting water areas［J］.Advances in Water Science，2008，19(1):78-83.(in Chinese))

［3］ Warner J C，Sherwood C R，Signell R P，et al.Development of a three-dimensional，regional，coupled wave，current，and sediment-transportmodel［J］.Computers ＆ Geosciences，2008，34(10):1284-1306.

［4］ Wang J H，Shen Y M.Development and validation of a three-dimensional，wave-current coupled model on unstructured meshes［J］.Science China Physics，Mechanics and Astronomy，2011，54(1):42-58.

［5］ Xie M X.Establishment，validation and discussions of a three dimensional wave-induced current model［J］.Ocean Modelling，2011，38(3):230-243.

［6］ Cockburm B，Shu C W.Runge-Kutta discontinuous Galerkin methods for convection-dominated problems［J］.Journal of Scientific Computing，2001，16(3):173-261.

［7］ Ford R，Pain C C，Piggott M D，et al.A nonhydrostatic finite-element model for three-dimensional stratified oceanic flows.Part I:model formulation［J］.Monthly Weather Review，2004，132(12):2816-2831.

［8］ Dawson C，Aizinger V.A discontinuous Galerkin method for three-dimensional shallow water equations［J］.Journal of Scientific Computing，2005，23(1):245-267.

［9］ Aizinger V，Dawson C.The local discontinuous Galerkin method for three dimensional shallow water Flow［J］.Computer Methods in Applied Mechanics and Engineering，2007，196(1):734-746.

［10］ Karna T，Legat V，Deleersnijder E.A baroclinic discontinuous Galerkin finite element model for coastal flows［J］.Ocean Modelling，2013，61:1-20.

［11］ Lin P Z，Zhang D.The depth-dependent radiation stresses and their effect on coastal currents［C］//Proceedings of the 6th International Conference of Hydrodynamics:Hydrodynamics VI Theory and Applications.2004，247-253.

［12］ Rodi W.Turbulence models and their application in hydraulics［M］.IAHR，Monograph Series，the Netherlands，1984.

［13］ Smagorinsky J S.General circulation experiments with the primitive equations I，the basic experiment［J］.Monthly Weather Review，1963，91:99-164.

［14］ Larson M，Kraus N C.Numerical model of longshore current for bar and trough beaches［J］.Journal of Waterway，Port，Coastal and Ocean Engineering，1991，117(4):326-347.

［15］ Oey L Y，Ezer T，Wang D P，et al.Loop current warming by hurricane Wilma［J］.Geophysical Research Letters，2006，33:L08613.

［16］ Longuet-Higgins M S.Longshore currents generated by obliquely incident sea waves，1［J］.Journal of Geophysical Research，1970，75(33):6778-6789.

［17］Booij N，Ris R C，Holthuijsen L H.A third-generation wave model for coastal regions:1.Model description and validation［J］.Journal of Geophysical Research，1999，104(C4):7649-7666.

［18］ Bassi F，Crivellini A，Rebay S，et al.Discontinuous Galerkin solution of the Reynolds-averaged Navier-Stokes and k-ω turbulence model equations［J］.Computers ＆ Fluids，2005，34:507-540.

［19］ Qiu J X，Khoo B C，Shu C W.A numerical study for the performance of the Runge-Kutta discontinuous Galerkin method based on different numerical fluxes［J］.Journal of Computational Physics，2006，212:540-565.

［20］ Bassi F，Rebay S.A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations［J］.Journal of Computational Physics，1997，131(1):267-279.

［21］ Bassi F，Rebay S.Numerical evaluation of two discontinuous Galerkin methods for the compressible Navier-Stokes equations［J］.International Journal of Numerical Methods in Fluids，2002，40(1-2):197-207.

［22］ Cockburn B，Shu C W.The local discontinuous Galerkin method for time-dependent convection-diffusion systems［J］.SIAM Journal on Numerical Analysis，1998，35(6):2440-2463.

［23］ Peraire J，Persson P O.The compact discontinuous Galerkin(CDG)method for elliptic problems［J］.SIAM Journal of Scientific Computing.2008，30(4):1806-1824.

［24］ Casulli V，Walters R A.An unstructured grid，three-dimensional model based on the shallow water equations［J］.International Journal for Numerical Methods in Fluids，2000，32(3):331-348.

［25］Ford R，Pain C C，Piggott M D，et al.A nonhydrostatic finite-element model for three-dimensional stratified oceanic flows.Part II:model validation［J］.Monthly Weather Review，2004，132(12):2832-2844.

［26］ Ting F C K，Kirby J T.Observation of undertow and turbulence in a laboratory surf zone［J］.Coastal Engineering，1994，24:51-80.

［27］ Dally W R，Brown C A.A modeling investigation of the breaking wave roller with application to cross-shore currents［J］.Journal of Geophysical Research，1995，100(C12):24873-24883.

［28］ Svendsen I A.Mass flux and undertow in a surf zone［J］.Coastal Engineering，1984，8(4):347-365.

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