基于間斷有限元方法的三維近岸波生流數值模型建立

趙張益，張慶河

(天津大學 水利工程仿真與安全國家重點實驗室，天津 300072)

波浪及其破碎所產生的波生流對近岸區域污染物和泥沙輸運等具有顯著影響。掌握近岸海域波生流場的分布規律，對海岸工程規劃設計以及海岸帶污染治理具有重要意義［1-2］。由于波生流垂向結構的非均勻性，發展三維波生流或近岸波流相互作用數值模型是近年來的主要趨勢。

目前，學者們通常采用比較常用的數值格式如有限差分(FD)法，有限元(FE)法，有限體積(FV)法等［3-5］來建立和發展三維波生流模型，其中，Warner等［3］在ROMS模型基礎上發展了三維波生流模型，但并未對其進行充分驗證，同時也沒有考慮波浪附加紊動效應;Wang等［4］在FVCOM模型基礎上建立了三維波生流模型，并對模型進行了驗證，不足在于模型沒有考慮波浪附加紊動效應;Xie［5］在ECOM模型基礎上建立了三維波生流模型，其優點在于考慮了波浪附加紊動效應，并對模型進行了大量的驗證。近年來，間斷有限元(DG)法因為其在守恒性、穩定性、處理復雜地形、網格自適應以及并行計算等方面所具有的優勢［6］，在三維水動力模型的數值求解中越來越得到重視［7-10］。然而，采用DG法的近岸波生流模型目前尚未見諸報道。為了給近岸波生流模擬提供高效準確的模擬工具，將采用DG法構建一個非結構化網格的三維波生流數值模型，并將其應用于三維近岸波生流模擬中以說明模型的合理性。

1 數學模型

1.1 三維波生流控制方程

三維波生流控制方程包含以下動量方程和連續方程:

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

三維輻射應力采用如下形式的公式［11］:

根據Warner等［3］的工作，波浪水滾剪切應力公式:

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

采用如下形式［14］計算波浪破碎引起的附加水平紊動系數

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

式中:α為水流水平紊動渦粘系數，在模型中取0.1;Λ為波浪水平摻混程度的經驗系數，將在后面的算例中給出其取值;Umax為近底波浪水質點軌跡速度水平分量的幅值;H為波高。

1.2 邊界條件

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

或Neumann邊界條件:

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

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

自由表面邊界條件由表面的Reynolds應力和風應力平衡可得:

式中:風應力可以表示為:

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

底邊界條件由底部Reynolds應力和底部切應力平衡可得:

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

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

1.3 波浪模型

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

2 數值離散

2.1 空間離散

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

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

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

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

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

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

推進最嚴格的水資源管理制度落實。抓好實行最嚴格水資源管理制度有關規定的落實和考核工作，開展中水回用項目建設規劃、水資源保護規劃和水中長期供求規劃工作，完成地下水的調查評估，劃定禁采區和限采區，并向社會公告。

式中:s為單元內部邊界的集合。

由此可以求得:

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

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

2.2 時間離散

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

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

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

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

通過上述時間推進，可以得到近岸波生流的數值解。

3 模型驗證及結果分析

在FLUIDITY純水動力模型［7］基礎上，添加波浪影響項，建立了三維近岸波生流數值模型，關于原始純水動力模型的驗證可以參閱文獻［25］，這里主要對所建立的近岸波生流模型進行驗證。

3.1 近底回流

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

采用SWAN模型模擬波浪傳播，由于在破波區域波高模擬結果偏小，所以在該區域采用擬合實測數據的方式對波高進行確定［27］，擬合曲線如圖2所示。

圖1 試驗布置示意Fig.1 Sketch of experiment

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

波浪增減水驗證如圖3所示。可以看出考慮波浪水滾剪切應力的計算值與實測值吻合較好，反映了波浪增減水過程，而不考慮波浪水滾剪切應力得到的最大減水點位置向破波帶外延伸，且增水值偏大。流速垂向分布驗證如圖4所示，其中測點坐標如表1所示。從流速垂向分布的計算結果可以看出，其數量級、垂向分布形態與實測值較為一致，但中間水體的計算流速梯度分布則相對比較均勻，其原因可能是由于數學模型中采用的三維輻射應力表達式是基于線性波理論，而實際情況中，由于較大波浪及破碎引起的強非線性效應使得二者之間得到的流速垂向結構有一定區別。模型整體上反映了近底回流的實際現象。

表1 流速測點坐標Tab.1 Coordinates of velocity stations

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

圖3 波浪增減水驗證Fig.3 Comparison between the simulated and observed wave setup

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

圖4 流速驗證Fig.4 Comparisons between the simulated and observed velocities

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

3.2 沿岸流

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

數值模型中，平面計算域34 m×16.6 m，布置三角形單元8 010個，節點4 131個;垂向每層高0.035 m。Λ取0.1，cf取0.018。造波端設置為開邊界，坐標原點位于岸線處。

圖6 試驗布置示意Fig.6 Sketch of experiment

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

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

沿岸流沿水深平均流速驗證結果如圖8所示，其中實測流速的測點坐標如表2所示，考慮波浪水滾剪切應力的沿岸流速計算值與實測值吻合程度較好，而不考慮波浪水滾剪切應力的沿岸流速峰值偏小，且峰值位置向離岸方向推移。沿岸流速垂向分布驗證如圖9所示，其中x軸坐標如表3所示，從圖中可以看出，個別測點處表層附近流速有一定差異，原因可能是由于模型沒有考慮波浪運動引起的附加垂向紊動效應。但總的來說，模擬得到的沿岸流速垂向分布與實測值較為相近，模型能夠合理地模擬波浪斜向入射而形成的沿岸流現象。

表2 水深平均流速測點坐標Tab.2 Coordinates of depth-averaged velocity stations

表3 流速垂向分布測點坐標Tab.3 Coordinates of velocity stations

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

4 結語

基于DG方法構建了非結構化網格三維波生流數值模型，并與近底回流和沿岸流物模實驗結果進行了對比，主要結論如下:

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

2)近底回流算例的模擬結果顯示，所建立的數值模型能夠合理反映波生流的垂向結構，模擬得到的破波帶附近流速的垂向變化規律與已有文獻對近底回流現象的描述一致，數值解與實測值基本吻合。

3)沿岸流算例的模擬結果表明，模型能夠合理模擬由于波浪斜向入射而形成的沿岸流現象，沿岸流速計算值與實測值吻合程度良好。

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

5)模擬結果表明，算法的構建是比較合理的，但在模型中需要進一步引入波浪運動引起的附加垂向紊動效應，使其能夠更好地描述近岸波生流場。另外，模型中有關參數如波浪紊動摻混強度的系數、底部摩擦系數等如何取值，還有待進一步研究。

［1］ 辛文杰.潮流、波浪綜合作用下河口二維懸沙數學模型［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.波流共存場中多向隨機波浪傳播變形數學模型［J］.水科學進展，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.