非靜壓方程與波作用譜模型的波浪傳播嵌套模擬

鄒國良，張慶河

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

選取合理的數(shù)學(xué)模型來模擬不同范圍和水深的波浪運(yùn)動，對準(zhǔn)確描述波浪傳播變形、獲得建筑物附近的設(shè)計波浪要素有重要的工程意義。工程中在描述波浪從大范圍到小范圍傳播變形時通常將基于相位平均的波作用譜模型與基于求解相位的沿水深積分模型進(jìn)行嵌套計算，即在大范圍海域(以下簡稱大模型)采用計算效率較高的波作用譜模型來模擬波浪傳播變形，在工程關(guān)心的小范圍區(qū)域如港池附近(以下簡稱小模型)則采用求解相位的模型計算。小模型的邊界條件由大模型提供。工程中常用的做法是由大模型推算出嵌套邊界處代表波浪要素以作為小模型的嵌套邊界輸入條件。這與嵌套邊界處的真實波況往往是不一致的，由此得到的小模型計算結(jié)果與實際情況必然會有出入。Tozer等［1-5］提出了更為合理的嵌套方法，即將波作用譜模型與淺水方程模型進(jìn)行空間上譜嵌套，研究結(jié)果表明，嵌套邊界條件采用沿邊界變化的實際波浪譜時，復(fù)雜地形中的波浪傳播變形模擬才能具有較高的精度。

近年來，越來越多的學(xué)者開始通過求解含非靜壓項的非線性淺水方程來模擬波浪的傳播變形［6-11］。非靜壓模型在垂向往往只需要2～3層［6-8］即可較好地描述波浪的強(qiáng)非線性和強(qiáng)色散性，使得近岸波浪變形模擬的精度和計算效率提高。已有的針對求解相位模型與求解相位平均模型的嵌套工作主要是將波作用譜模型與平面二維的水深積分模型［1-4］及靜壓假定的淺水方程模型［5］進(jìn)行嵌套，尚未見到波作用譜模型與三維波浪模型及非靜壓模型嵌套的文獻(xiàn)。因此，為了更高效、準(zhǔn)確地描述波浪從大尺度范圍傳播至中小尺度范圍的傳播變形，并將嵌套模型拓展到三維，本文將波作用譜模型和非靜壓淺水方程波浪模型進(jìn)行嵌套，嵌套邊界條件采用由波作用譜模型提供的空間變化譜作為非靜壓方程的波浪輸入條件。

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

本文的波作用譜大模型模擬采用目前應(yīng)用非常廣泛的第3代風(fēng)浪模型SWAN模型，求解相位的小模型模擬采用考慮非靜壓項的非線性淺水方程模型SWASH。下面針對兩種模型做簡要的介紹。

1.1 波作用譜模型

SWAN模型采用波作用譜方程描述風(fēng)浪生成及其在近岸區(qū)的演化過程。關(guān)于SWAN模型的控制方程詳見文獻(xiàn)［13］以及 SWAN 技術(shù)手冊［14］，這里不再贅述。

SWAN模型計算域的邊界類型分為陸地和水邊界，其中陸地邊界為波能吸收邊界，水邊界可根據(jù)現(xiàn)場實測和計算結(jié)果給定。針對多向不規(guī)則波模擬，這里在迎浪水邊界采用JONSWAP頻譜［15］和光易型方向分布函數(shù)的乘積作為輸入的方向譜，即

式中:Hs為有效波波高;Tp和fp分別為譜峰周期和譜峰頻率;f為組成波的頻率;γ為譜峰因子;ms為方向集中度;θ0為波浪傳播的主波向。

1.2 非靜壓淺水方程模型

忽略黏性影響，同時假定水體密度恒定為ρ0，SWASH模型的控制方程可表示為

式中:u、v和w分別為沿x、y和z方向的流速;η為自由面;d為靜水深;h=η+d為總水深;g為重力加速度;壓力p分解成靜水壓力g(η-z)和非靜水壓力q。自由面和底部的運(yùn)動學(xué)邊界條件分別為

自由面可通過由連續(xù)性方程(4)沿水深積分，并結(jié)合式(8)獲得的自由面方程求解:

式中:Qx、Qy分別為x和y方向的流量，可表示為

忽略風(fēng)應(yīng)力和表面張力，自由面處的邊界條件為

底部邊界條件為

固邊界采用自由滑移的邊界條件，以平行于y軸的邊界為例:

入流邊界通過指定速度分布并結(jié)合弱反射邊界條件［16］吸收二次反射波。在模擬不規(guī)則波時，可將方向譜按照頻率分為M份，方向分為N份，共MN個諧波疊加，其速度邊界為

出流邊界采用改進(jìn)的海綿層吸收邊界［17-19］，海綿層內(nèi)每個時間步內(nèi)水平向速度采用顯式阻尼消波，阻尼系數(shù)可表示為

2 數(shù)值方法

2.1 波作用譜模型的數(shù)值方法

SWAN模型方程求解采用有限差分格式進(jìn)行離散。在空間上采用一階迎風(fēng)格式離散;在譜方向上采用迎風(fēng)格式與中心差分相結(jié)合的混合格式離散，并取等間距方向步長Δθ;在頻率上采用符合對數(shù)頻率分布的等間距相對頻率步長Δf/f。時間離散采用歐拉全隱格式，具體的離散形式可參照文獻(xiàn)［14］。

2.2 非靜壓淺水方程模型的數(shù)值方法

非靜壓淺水方程模型SWASH基于交錯網(wǎng)格系統(tǒng)對沿垂向各層積分的控制方程采用具有二階精度的顯式有限差分格式進(jìn)行空間離散。動量方程中的對流項采用具有守恒特性的對流項離散格式進(jìn)行離散［19］。為了減少垂向的分層數(shù)同時保證模型能較好的描述短周期波的色散特性，動量方程中的非靜壓梯度項采用緊致Keller-box格式進(jìn)行離散，并將壓力項定義在單元的邊中心(如圖1)，從而使表層自由面的壓力邊界自動滿足，而不需要進(jìn)行假定［7］。對于時間的積分，連續(xù)性方程和動量方程均采用二階蛙跳格式進(jìn)行時間積分。此外，模型還采用了干-濕處理技術(shù)來描述動邊界。對于SWASH模型控制方程的詳細(xì)數(shù)值離散過程可參照文獻(xiàn)［6，18］

圖1 SWASH交錯網(wǎng)格變量布置Fig.1 Arrangement of the unknown with a staggered grid system in SWASH

模型在求解壓力速度耦合方程時采用分步法。在靜壓步僅考慮非靜壓梯度項中的顯式部分，求解新時刻的過渡流速;在非靜壓步求解局部質(zhì)量守恒推導(dǎo)出的壓力泊松方程獲得非靜壓項修正項，并對過渡流速進(jìn)行修正，最終通過自由面方程更新自由面。

2.3 嵌套方法

本文對2種波浪模型進(jìn)行的嵌套計算并不考慮由SWASH模型所計算的反射波對SWAN模型的影響。嵌套過程需要解決2個問題:1)采用非靜壓淺水方程進(jìn)行多向不規(guī)則波模擬時，為了使模擬的方向譜在空間上較為均勻，同時為了滿足平穩(wěn)性要求，方向譜按照頻率分割時的數(shù)量一般要大于波作用譜模型所要求的數(shù)量;2)SWASH模型所需要的網(wǎng)格尺寸要遠(yuǎn)小于SWAN模型。為此，在嵌套時先采用3次樣條插值函數(shù)重新插值SWAN模型計算出嵌套邊界上的每一個方向上的頻率譜［2?，再將新插值的網(wǎng)格尺寸較粗節(jié)點處的方向譜線性插值到網(wǎng)格尺寸相對較小的SWASH模型的邊界處。

SWASH模型嵌套邊界處的方向譜采用傅里葉逆變換將其轉(zhuǎn)換成多組諧波，由此可獲得各組成諧波的波浪參數(shù)，如振幅、相位、圓頻率等，進(jìn)而根據(jù)式(15)可獲得SWASH模型的輸入條件。

3 數(shù)值計算結(jié)果與分析

為了說明所提出的嵌套方法的合理性，設(shè)計了兩組多向不規(guī)則波在理想地形上傳播的算例進(jìn)行數(shù)值模擬試驗。兩組算例中的波浪條件以及水深地形均不考慮破碎情況的發(fā)生。

3.1 等水深波浪傳播算例

算例1中大模型的計算地形為1 km×1 km，水深為-5m，網(wǎng)格尺寸為Δx=Δy=20m，計算地形如圖2所示。多向不規(guī)則波的有效波高Hs為1 m，譜峰周期為8 s，譜峰因子為3.3，頻率分布范圍為0.08～1.0 Hz，按照頻率分為28等份，波浪入射邊界為南邊界，主波向為90°，方向分布范圍為75°～105°，按照方向分為21份，方向譜集中度取2。小模型計算域為400 m×350 m，南邊界(圖2中bnd)為嵌套邊界，北邊界設(shè)置100 m寬的海綿層進(jìn)行消波，東西側(cè)邊界為自由滑移邊界。模型計算網(wǎng)格尺寸為Δx=Δy=2 m，時間步長Δt=0.1 s，總計算時間為45 min。

圖2 等水深計算地形Fig.2 Com putation bathymetry of hydroisobath

圖3比較了由SWAN模型以及嵌套后SWASH模型計算出點P1～P3位置處的能量譜。表1為計算出的有效波高和波譜平均周期(Tm01)對比結(jié)果。其中，SWASH模型計算出時序列的水位通過快速傅里葉變化進(jìn)行了能量譜估計，采樣時間間隔為0.1 s，采樣時長15 min。根據(jù)計算結(jié)果可看出，由采用譜作為嵌套邊界輸入條件的淺水方程模型計算出的能量譜、波高和波周期與用SWAN計算出的結(jié)果基本一致。

圖3 P1～P3處SWAN模型和SWSAH模型計算譜比較Fig.3 Comparison of SWAN and SWASH com puted energy density spectra at P1～P3

表1 P1～P3處SWAN和SWASH計算波浪參數(shù)比較Table 1 Comparison of SWAN and SWASH computed wave parameters at P1～P3

3.2 斜坡上的波浪傳播算例

為了進(jìn)一步體現(xiàn)嵌套模型在模擬波浪發(fā)生折射變形時的準(zhǔn)確性，本文設(shè)置了水深以1/400坡度從-10 m變化至-5 m的地形，如圖4所示。為了消除SWAN計算時側(cè)邊界的影響從而得到合理的嵌套邊界波浪要素，大模型計算范圍為10 km×2 km，網(wǎng)格尺寸為Δx=Δy=20 m。多向不規(guī)則波的有效波高、譜峰周期以及方向譜等相關(guān)參數(shù)同3.1節(jié)，波浪入射邊界為南邊界，入射主波向為45°，方向分布為0°～90°，方向集中度為15。

圖4 斜坡水深地形Fig.4 Sketch of linear beach bathymetry

圖5 兩種嵌套模式計算域Fig.5 Plan view of the com putation domain for two types nestingmode

模型嵌套采用2種模式:1)嵌套邊界取為南邊界bnd1和西邊界bnd2，北、東邊界為150 m寬的海綿層消波，如圖5(a)所示;2)采用工程中常見做法即波浪正向入射嵌套邊界，即取嵌套邊界bnd3垂直于入射波向，其他邊界均設(shè)置為150 m寬的海綿層，如圖5(b)所示。2種模式的計算范圍均為1.7 km×1 km，網(wǎng)格尺寸為Δx=1 m，Δy=1 m，垂向等分為2層，時間步長Δt=0.025 s，總計算時間為50 min。

圖6為大模型提供的沿嵌套邊界的波高分布圖以及平均波向圖。由圖可看出:采用第1種嵌套模式時，波高和波向沿邊界bnd1的變化不大，而沿邊界bnd2變化較大;采用第2種嵌套模式時，波高沿嵌套邊界bnd3的變化較明顯，而平均波向沿邊界bnd3基本在48°左右。

圖6 沿嵌套邊界的有效波高和平均波向分布Fig.6 Significant wave height and mean direction along nesting boundaries

圖7、8比較了2種嵌套模式條件下大模型以及小模型嵌套計算出的沿同一位置S1～S2以及S3～S4的有效波高分布以及能量譜的對比結(jié)果(圖中MODE1表示嵌套模式1，MODE2為嵌套模式2)。

圖7 2種嵌套模式下SWAN和SWASH模型計算出有效波高分布對比Fig.7 Comparison of SWAN and SWASH computed wave significant wave height for two types of nestingmode

圖8 S1～S4處SWAN模型和SWSAH模型計算譜比較Fig.8 Comparison of SWAN-and SWASH-computed energy density spectra at S1～S4

根據(jù)上述波浪在等水深地形傳播和在斜坡上進(jìn)行淺水折射變形的算例結(jié)果，嵌套后的非靜壓模型計算波高、周期與波作用譜模型接近，且兩者計算出的譜比較吻合。這說明所建立的嵌套模型可以合理描述波浪從大范圍區(qū)域向小范圍區(qū)域的傳播，所建立的嵌套模型可為工程區(qū)小范圍波浪模擬(如港內(nèi)波高分布等)提供一種新的計算模式。

表2則給出相應(yīng)的計算波參數(shù)對比結(jié)果.根據(jù)圖表結(jié)果可知，由兩種嵌套模式分別計算出的有效波高沿程分布與采用波作用譜模型計算出的有效波高均基本一致，能量譜也基本接近，部分位置處的波高變動主要是因為波浪在斜坡上傳播會產(chǎn)生部分反射所造成。

表2 S1～S4處SWAN和嵌套模式計算出的波浪參數(shù)比較Table 2 Comparison of SWAN and nesting mode computed wave parameters at S1～S4

4 結(jié)束語

通過波作用譜模型提供給非靜壓模型嵌套邊界處的波浪譜參數(shù)，并對其在空間和頻率方向上進(jìn)行插值使其滿足求解相位模型的邊界條件，較好地實現(xiàn)了非靜壓淺水方程模型與波作用譜模型的嵌套。波浪在等水深地形傳播和斜向入射波在斜坡上傳播的2個算例結(jié)果表明，嵌套模型可以獲得較為準(zhǔn)確的波高、周期以及波浪譜，可以合理描述波浪從大范圍區(qū)域向小范圍區(qū)域的傳播，在解決類似于港內(nèi)波高分布等工程波浪計算問題時，可以獲得更高精度的結(jié)果。非靜壓波浪模型與譜作用模型的嵌套有著十分廣闊的應(yīng)用前景。

［1］TOZER N P，DURAND N.Improved coastal wave modelling using coupled spectral wave transformation and wave disturbance models［C］//Proceedings of 28th International Conference on Coastal Engineering.Wales，UK，2002:448-458.

［2］OKIK，SAKAIT.Coupling of phase-averaging model and phase-resolvingmodel by using double connection boundaries［C］//Proceedings of 31st International Conference on Coastal Engineering.Hamburg，Germany，2008:302-313.

［3］CHENW，DEMIRBILEK Z，LIN L.Coupling phase-resolving nearshore wave model with phase-averaged spectral wave models in coastal applications［C］//Proceedings of 11st International Conference on Estuarine and CoastalModeling.Seattle，Washington，USA，2009:582-600.

［4］MCCUBEM.Coupled wave action and shallow-watermodelling for random wave runup on a slop［J］.Journal of Hydraulic Research，2011，49(4):515-522.

［5］RUSU E，SOARES C G.Modeling waves in open coastal areas and harbourswith phase-resolving and phase-averaged models［J］.Journal of Coastal Research，2013，29(1):1309-1325.

［6］ZIJLEMA M，STELLINGG S，SMITP.SWASH:An operational public domain code for simulating wave fields and rapidly varied flows in coastalwaters［J］.Coastal Engineering，2011，58(10):992-1012.

［7］ZIJLEMA M，STELLING G S.Further experiences with computing non-hydrostatic free-surface flows involvingwater waves［J］.International Journal for Numerical Methods in Fluids，2005，48(2):169-197.

［8］ZIJLEMA M，STELLING G S.Efficient computation of surf zone waves using the nonlinear shallow water equations with non-hydrostatic pressure［J］.Coastal Engineering，2008，55:780-790.

［9］AIC，JIN S，LV B.A new fully non-hydrostatic 3D free surface flowmodel forwaterwavemotions［J］.International Journal for Numerical Methods in Fluids，2011，66(11):1354-1370.

［10］STELLING G S，ZIJLEMA M.Numerical modeling of wave propagation，breaking and run-up on a beach［J］.Computational Science and Engineering，2010，71:373-401.

［11］SMIT P，ZIJLEMA M，STELLING G S.Depth-induced wave breaking in a non-hydrostatic，near-shorewavemodel［J］.Coastal Engineering，2013，76:1-16.

［12］BOOIJN，RISR C，HOLTHUIJSEN LH.A third generation wavemodel for coastal regions，part1:model description and validation［J］.Journal of Geophysical Research，1999，104(C4):7649-7666.

［13］HOLTHUIJSEN L H，A HERMAN，BOOIJN.Phase-decoupled refraction-diffraction for spectralwavemodels［J］.Coastal Engineering，2003，49:291-305.

［14］SWAN Team.SWAN Usermanual(version 40.85)［M］.Delft:Delft University of Technology，2011.

［15］GODA Y.A comparative review on the functional forms of directionalwave spectrum［J］.Coastal Engineering，1999，41(1):1-20.

［16］BLAYO E，DEBREU L.Revisiting open boundary conditions from the point of view of characteristic variables［J］.Ocean Modelling，2005，9:231-252.

［17］ZOU Guoliang，ZHANG Qinghe.Improvement of absorbing boundary conditions for non-hydrostatic wave-flow model SWSAH［J］.Applied Mechanics and Materials，2013(353/354/355/356):2676-2682.

［18］鄒國良，張慶河.非靜壓波浪模型無反射造波［J］.海洋工程，2012，30(4):55-61.ZOU Guoliang，ZHANG Qinghe.Wave generation without re-reflection for non-hydrostatic wavemodel［J］.The Ocean Engineering，2012，30(4):55-61.

［19］STELLING G S，DUINMMEIJER SP A.A staggered conservative scheme for every Froude number in rapidly varied shallow water flows［J］.International Journal for Numerical Methods in Fluids，2003，43:1-23.