基于Boussinesq 水波模型的聚焦波模擬

劉必勁，張振偉，劉忠波，傅丹娟，陳小云

(1.廈門理工學(xué)院 土木工程與建筑學(xué)院，福建 廈門 361024；2.大連海事大學(xué) 交通運(yùn)輸工程學(xué)院，遼寧 大連 116026)

1 引言

畸形波是一種非線性很強(qiáng)的波浪，其波高很大，通常是有效波高的2.2 倍以上，這種波浪對海上建筑物或過往的船只帶來嚴(yán)重威脅。為了深入掌握畸形波的形成機(jī)制，國內(nèi)外學(xué)者通過理論分析、數(shù)值模擬或?qū)嶒?yàn)?zāi)M等方式展開了數(shù)十年的研究。2003 年，Kharif 和Pelinovsky[1]較為全面地分析了其形成機(jī)制，認(rèn)為在相位聚焦、波?流相互作用、地形變化、表面風(fēng)、波浪不穩(wěn)定性等影響下都可能產(chǎn)生畸形波。

為了研究畸形波，實(shí)驗(yàn)室中一般都采用相位聚焦的方式產(chǎn)生較大的波浪，如Baldock 等[2]在二維水槽中采用線性相位聚焦的方式產(chǎn)生了聚焦波；柳淑學(xué)和洪起庸[3]在三維水池中開展了三維極限波浪的實(shí)驗(yàn)研究。在理論分析方面，裴玉國等[4]采用線性模型模擬了定點(diǎn)生成畸形波的方法；趙西增等[5]則借助于聚焦波和不同形式波浪的組合模擬畸形波。在數(shù)值模擬方面，國內(nèi)外學(xué)者則通過利用N-S 方程模型、非靜壓模型、Boussinesq 方程數(shù)值模型、勢流方程模型、高階譜模型等不同的數(shù)學(xué)模型模擬聚焦波。如趙西增等[6]利用高階譜方法建立了模型，并應(yīng)用于畸形波的模擬；Li 和Liu[7]建立了非周期邊界的高階譜聚焦波模型，并研究設(shè)定聚焦波幅、中心頻率、頻率寬度、波幅分布以及方向角等對聚焦波的影響，研究表明，寬的方向譜減弱了非線性聚焦波幅；Ai 等[8]建立了不同網(wǎng)格下的非靜壓模型，并模擬了水槽和水池中的聚焦波；Ning 等[9]基于勢流理論建立了高階邊界元數(shù)值模型，模擬了深水聚焦波；Li 等[10]在CIP 模型（利用CIP 求解對流項(xiàng)）中引入內(nèi)部造波法，模擬了潛堤上聚焦波的傳播變形問題。

關(guān)于Boussineq 模型在聚焦波模擬的研究文獻(xiàn)相對較少，能成功模擬聚焦波速度分布特征的僅有少數(shù)Boussinesq 水波方程，這是因?yàn)槎鄶?shù)Boussinesq 方程本身在色散性、非線性，特別是速度分布特征方面存在嚴(yán)重不足。分析相關(guān)數(shù)值模擬聚焦波的研究可知，準(zhǔn)確模擬波面?zhèn)鞑プ冃涡枰獢?shù)值模型具有良好的色散性與非線性特征，準(zhǔn)確模擬速度分布必須要求模型具備良好的速度分布特征。當(dāng)前能具備這些性能的Boussineq 數(shù)值模型并不多，在分析大量Boussineq方程文獻(xiàn)基礎(chǔ)上，本文選用了一組單層Boussinesq 水波方程，即Liu 和Fang[11]給出的雙層Boussineq 模型的簡化版。事實(shí)上，該簡化版本的模型在平底情況下與Madsen 等[12]給出的Boussineq 水波方程完全一致。本文將針對單層Boussinesq 水波方程建立數(shù)值模型，并考察其是否能適用于模擬強(qiáng)非線性波浪傳播和聚焦波的演化問題。

2 Boussinesq 數(shù)值模型

2.1 控制方程

假設(shè)水體無旋、無黏，在保留垂向速度下，Liu 和Fang[11]經(jīng)嚴(yán)格數(shù)學(xué)推導(dǎo)給出了適合緩變地形的雙層Boussinesq 方程，將其簡化成為單層方程。立面二維模型由7 個(gè)方程組成，含有7 個(gè)待求變量：波面位移η，自由表面處水平和垂向速度uη和wη，靜止水位處的水平和垂向速度u0和w0以及計(jì)算水平和垂向速度u*和w*，方程的具體表達(dá)形式如下：

（1）自由面 η處的運(yùn)動(dòng)學(xué)邊界條件和動(dòng)力學(xué)邊界條件為

（2）自由面速度（uη，wη）與靜止水位處速度（u0，w0）的關(guān)系為

（3）靜止水位處速度（u0，w0）與計(jì)算速度（u*，w*）的關(guān)系為

（4）水底處滿足的運(yùn)動(dòng)學(xué)邊界條件為

式中，g是重力加速度；下標(biāo)x表示對x求導(dǎo)；h表示水深；其中(α,β12,β13)=(0.5,0.95,0.55)。

整個(gè)流場的表達(dá)式分為兩部分，從靜止水位到波面處的速度表達(dá)式為

式中，z表示原點(diǎn)位于靜水面時(shí)向上為正的垂向位置。從水底到靜止水位的速度表達(dá)式為

方程（1）至方程（7）構(gòu)成單層Boussinesq 水波方程，圖1 給出了變淺梯度與解析解的對比以及無因次相速度。由圖可見，在2%誤差范圍內(nèi)，相速度可適用于最大水深k0h=10（k0為深水波數(shù)，h為水深）[12]；速度分布特征在1%誤差范圍內(nèi)適用的最大水深為k0h=3.5；變淺性能可適用0＜k0h＜10 范圍①林鵬程，劉忠波，劉勇.基于Boussinesq 數(shù)值模型的波浪速度垂向分布模擬研究.海洋湖沼通報(bào)，已采用.。

2.2 數(shù)值模型的求解

以上垂向二維Boussinesq 方程中的最高空間導(dǎo)數(shù)是3 階，模型本身與開源代碼FUNWAVE 的最高空間導(dǎo)數(shù)是一致的[13]，為此時(shí)空差分格式均采用與該開源代碼一樣的格式。在時(shí)間差分格式上采用混合4 階Adams-Bashforth-Moulton 格式，預(yù)報(bào)時(shí)利用3 階Adams-Bashforth 格式求解方程（1）和（2），可以得到波面位移以及相應(yīng)的波面處水平速度的預(yù)報(bào)值；分別求解方程（3）和（5）可得到2 個(gè)水平速度的預(yù)報(bào)值；求解方程（7）、（6）和（4），可依次得到3 個(gè)垂向速度的預(yù)報(bào)值。校正時(shí)利用4 階Adams-Moulton 時(shí)間步進(jìn)格式，求解得到所有變量的校正值。當(dāng)7 個(gè)變量的校正值與預(yù)報(bào)值之間的誤差在0.000 1 內(nèi)，當(dāng)前迭代結(jié)束，否則用校正值和預(yù)報(bào)值的加權(quán)重新賦值給變量的預(yù)報(bào)值，本文算例中校正加權(quán)值R為0.1～0.2，而預(yù)報(bào)加權(quán)為1-R，并多次重復(fù)校正過程，直到滿足條件為止。以上過程可以參見文獻(xiàn)[13-14]。

圖1 方程的變淺梯度（a）和無因次相速度（b）Fig.1 The shoaling gradient (a) and non-dimensional (b)phase celerity of the present model

3 數(shù)值模型的驗(yàn)證

為了驗(yàn)證數(shù)值模型，我們采用流函數(shù)波浪的數(shù)值解作為近似的解析解來考察模型的適用情況。數(shù)值計(jì)算長度設(shè)置為10L（L為流函數(shù)波浪的波長），在入射邊界2 個(gè)波長范圍內(nèi)用流函數(shù)波浪的解析解作為模擬的驅(qū)動(dòng)條件，并在該范圍內(nèi)采用了Fuhrman[15]提出的邊界松弛造波法，這相當(dāng)于將近似解析解逐漸轉(zhuǎn)化為方程數(shù)值解的一種方式。在出口邊界，設(shè)置了2L的海綿邊界層進(jìn)行消波。在固定波浪周期為6 s情況下，針對水深為20 m、30 m 和40 m 情況下，對波高為6 m 的流函數(shù)波浪進(jìn)行了數(shù)值模擬研究。3 種情況對應(yīng)的流函數(shù)波長分別為60.469 m、61.524 m 和61.679 m，無因次水深kh分別為2.078、3.064 和4.074；對應(yīng)的波陡H/L=0.099、0.098 和0.097；這幾組工況的波高約是極限波高的0.7、0.65 和0.64 倍左右，因此屬于強(qiáng)非線性波浪算例。

圖2 水深為20 m 的計(jì)算結(jié)果與解析結(jié)果的比較Fig.2 Comparisons of analytical solution and numerical simulation for water depth 20 m

圖3 水深為20 m 波峰面下的水平速度分布與解析解的比較（x=4L）Fig.3 Comparisons of velocity profile between numerical simulation and analytical solution for water depth 20 m (x=4L)

每一個(gè)工況計(jì)算時(shí)間為300 s，空間計(jì)算步長取為L/32，時(shí)間計(jì)算步長取為0.05 s。將計(jì)算得到的波面位移、波面處水平速度與垂向速度以及沿垂向分布的水平速度與近似解析解（流函數(shù)波浪的數(shù)值解析解）進(jìn)行了比較，具體結(jié)果見圖2 至圖7。圖2、圖4和圖6 給出了計(jì)算波面位移、波面處水平速度與垂向速度與近似解析解的比較；圖3、圖5 和圖7 給出了x=4L波峰面下水平速度的垂向分布與解析解的比較。由圖2、圖4 和圖6 可見，波面以及波面處速度與解析解誤差較小，為了明確誤差的具體值，計(jì)算了波長與解析解的誤差（從2L～8L范圍內(nèi)標(biāo)定），3 種情況下僅第一組誤差是0.52%，其他組幾乎不存在誤差；以2L～8L范圍內(nèi)變量最大值的平均值為統(tǒng)計(jì)，最大波峰面與解析解的誤差、最大水平速度、最大垂向速度與相應(yīng)的解析解誤差均在1%范圍內(nèi)。進(jìn)一步，從x=4L波峰面下水平速度（圖3，圖5，圖7）對比來看，伴隨水深的加大，波峰面水平速度分布的精確度逐漸變差，這與文獻(xiàn)中給出的線性Stokes 波的速度分布與解析解在kh=3.5 的積分誤差是1%的規(guī)律是類似的，這里不再進(jìn)行速度分布的積分誤差的詳細(xì)計(jì)算。

圖4 水深為30 m 的計(jì)算結(jié)果與解析結(jié)果的比較Fig.4 Comparisons of analytical solution and numerical simulation for water depth 30 m

4 聚焦波演化的模擬

圖5 水深為30 m 波峰面下的水平速度分布與解析解的比較（x=4L）Fig.5 Comparisons of velocity profile between calculated and analytical solution for water depth 30 m (x=4L)

受水槽長度的限制，實(shí)驗(yàn)室中一般采用線性相位聚焦方式產(chǎn)生較大的聚焦波進(jìn)行物理模型實(shí)驗(yàn)。在深水情況下，由于受波浪非線性的影響導(dǎo)致實(shí)際聚焦位置與設(shè)置的線性聚焦位置存在向后偏移，對應(yīng)的聚焦時(shí)間也會向后偏移。這種偏移是因?yàn)榫劢共ǖ牟ǚ龃蠛螅€性波的波長不能詮釋聚焦波波幅增加后引起波長增加的現(xiàn)象，其本質(zhì)就是線性色散關(guān)系表達(dá)式已經(jīng)不適用于闡釋強(qiáng)非線性的波長，此時(shí)應(yīng)考慮波幅引起色散關(guān)系的變化，即波幅離散引起頻散。關(guān)于這一點(diǎn)在上面的算例中已有相關(guān)印證，以第3 個(gè)工況為例，線性波長為56.192 m，而流函數(shù)波浪的波長為61.679 m，二者誤差為9.7%。聚焦波的非線性增強(qiáng)，則這種偏移量越大[2]。

圖6 水深為40 m 的計(jì)算結(jié)果與解析結(jié)果的比較Fig.6 Comparisons of analytical solution and numerical simulation for water depth 40 m

圖7 水深為40 m 波峰面下的水平速度分布與解析解的比較（x=4L）Fig.7 Comparisons of velocity profile between numerical simulation and analytical solution for water depth 40 m (x=4L)

Baldock 等[2]開展的深水實(shí)驗(yàn)含有不同波浪頻率段范圍，在不同周期內(nèi)等分周期成29 份，每個(gè)頻率的波幅取值相等，其中B 組（T=0.6～1.4 s，寬譜）和D 組（T=0.8～1.2 s，窄譜）作為主要參考用來考察本文數(shù)值模型，B 組和D 組的線性無因次水深kh分別為1.568～7.825 和2.026～4.403。為了模擬聚焦波的生成過程，本文在入射邊界處引入線性解疊加方式，具體做法：給定每一個(gè)組成成分的波面位移的時(shí)間歷程，利用線性化的Boussinesq 方程確定出其他所有變量的時(shí)間歷程值，詳細(xì)過程可參見文獻(xiàn)[14,16]。設(shè)置計(jì)算域?yàn)?4 m，空間計(jì)算步長采用0.04 m，時(shí)間計(jì)算步長為0.01 s，末端設(shè)置4 m 的海綿邊界層消波，計(jì)算時(shí)間為40 s。對B 組和D 組的3 種不同設(shè)定的聚焦波峰面為22 mm、38 mm 和55 mm 工況進(jìn)行數(shù)值模擬，計(jì)算的聚焦位置處的波面位移和最大波峰面下的速度分布與實(shí)驗(yàn)結(jié)果進(jìn)行了比較，具體見圖8 和圖9。整體而論，計(jì)算波面位移與實(shí)驗(yàn)結(jié)果吻合程度較好，水平速度的分布與實(shí)驗(yàn)結(jié)果基本吻合。寬譜情況下（B 組），計(jì)算的水平速度出現(xiàn)“S”形，圖9 與圖7 的結(jié)果有相類似的規(guī)律，這是因?yàn)閿?shù)值模型在線性無因次水深kh＞3.5 后誤差將逐漸增大，在kh＜3.5 只有18 個(gè)組成波在這個(gè)范圍內(nèi)；窄譜情況下（D 組），線性無因次水深kh＜3.5 范圍內(nèi)22 個(gè)組成波，因此D 組工況比B 組工況模擬效果要好。關(guān)于D22 算例數(shù)值結(jié)果與實(shí)驗(yàn)結(jié)果吻合程度不好的原因，本文也嘗試了采用兩個(gè)不同的雙層Boussinesq 水波方程模型[14,16]；Liu 等[14]數(shù)值模擬結(jié)果與本文計(jì)算結(jié)果相差不大，也存在與實(shí)驗(yàn)偏差較大的情況，鑒于這一組的波浪非線性相對較弱，因此可以初步認(rèn)為有可能是實(shí)驗(yàn)結(jié)果存在一定的偏差。

對于設(shè)定聚焦波面為55 mm 的情況，固定波浪中心頻率不變，即保持1 s 周期不變，設(shè)定線性聚焦位置為8 m，30 s。在周期范圍△T分別為0.8 s、0.7 s、0.6 s、0.5 s 和0.4 s 情況下聚焦波時(shí)的波峰面、波峰面速度、聚焦時(shí)間和空間差等在表1 中給出，時(shí)空偏移是數(shù)值模型（或?qū)嶒?yàn)?zāi)Ｐ停┑玫降膶?shí)際聚焦位置與線性波群設(shè)定的聚焦時(shí)間和位置之間的偏移量。

5 結(jié)論

本文建立了基于單層Boussinesq 水波方程的垂向二維數(shù)值模型，首先與流函數(shù)波浪解析解進(jìn)行了對比，并將這個(gè)數(shù)值模型應(yīng)用于聚焦波的模擬中，研究得出主要結(jié)論如下：

圖8 聚焦的計(jì)算波面位移與實(shí)驗(yàn)波面位移的比較Fig.8 Comparisons of surface elevations between modeled and experimental data

圖9 聚焦波峰面下的水平速度與實(shí)驗(yàn)結(jié)果的比較Fig.9 Comparisons of velocity profiles between modeled and experimental data

表1 計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的比較Table 1 Comparisons of results between modeled and experimental data

（1）數(shù)值模型計(jì)算的波面位移、波面位移處的水平速度與垂向速度均與流函數(shù)的解析解有較好的吻合，同時(shí)垂向速度在適用范圍內(nèi)與解析解的比較，也展示出本文建立的數(shù)值模型能夠勝任模擬強(qiáng)非線性波浪。

（2）聚焦波演化中，波浪非線性性起重要作用，本文Boussinesq 型方程具有良好的色散性與非線性性能，保證了方程具有模擬聚焦波的能力；數(shù)值模擬波面與實(shí)驗(yàn)結(jié)果的對比，展現(xiàn)了本文Boussinesq 數(shù)值模型能夠模擬聚焦波波面演化。

（3）與精確模擬的波面相比而言，在模擬速度分布的精度方面，在理論方面，方程的色散適用范圍明顯大于水深速度分布的適用范圍。數(shù)值模型受限于線性速度分布的精度，其模擬非線性速度分布的能力也明顯受限于kh＜3.5。盡管數(shù)值模型可以較為精確模擬波面以及波峰面處的水平速度與垂向速度，但對于kh＞3.5 的組成成分較多的聚焦波，基于本文Boussinesq 水波方程的數(shù)值模型無法精準(zhǔn)模擬水平速度分布情況。

本文明確了水平速度分布的非線性適用范圍，在此基礎(chǔ)上可期望將模型拓展到研究聚焦波與建筑物相互作用問題。此時(shí)只需將速度表達(dá)式（公式（9）至公式（12））代入垂向歐拉方程中，進(jìn)行兩次垂向積分，最終可以得到波浪力，關(guān)于這一方面將在下一步進(jìn)行研究。