大范圍波浪傳播數(shù)學(xué)模型

楊耀中，馮衛(wèi)兵

（河海大學(xué)海洋學(xué)院，南京210098）

南黃海輻射沙脊群位于淺海區(qū)，海域開闊，有折射、繞射、破碎等波浪現(xiàn)象。為了研究波浪機(jī)理，科研人員開發(fā)出大量的波浪淺水變形數(shù)學(xué)模型［1］。一般說(shuō)來(lái)，可按照計(jì)算水域的范圍，把數(shù)學(xué)模型劃分為小、中、大三類。第一類是小范圍模型，以三維N-S方程為代表。此類模型可以求解完整的三維N-S方程，具有精度高的特點(diǎn)［2-3］。但由于計(jì)算機(jī)容量和速度的限制，目前以小范圍應(yīng)用為主。第二類是中等范圍模型，以Boussinesq方程和Berkhoff緩坡方程為代表。Boussinesq方程描述的是波面隨時(shí)間的變化過(guò)程，在數(shù)值求解時(shí)一般時(shí)間步長(zhǎng)取周期的1/24～1/30，空間步長(zhǎng)取波長(zhǎng)的1/8～1/12。因此此類方程只適合于中等范圍計(jì)算，如港池內(nèi)波浪場(chǎng)的計(jì)算［4］。同樣Berkhoff緩坡方程由于待求量為波勢(shì)函數(shù)，需要在一個(gè)波長(zhǎng)范圍內(nèi)布置8～10個(gè)計(jì)算點(diǎn)，因此在有限的內(nèi)存條件下，計(jì)算范圍也受到限制［5］。第三類是以動(dòng)譜平衡方程為代表的大范圍模型。此類模型雖然可以計(jì)算大范圍波浪場(chǎng)，但前提條件是加大空間步長(zhǎng)，降低網(wǎng)格分辨率［6］。而輻射沙脊群海域面積十分廣闊，達(dá)到140 km×200 km，同時(shí)地形變化劇烈。用加大空間步長(zhǎng)的方法計(jì)算很可能造成模擬結(jié)果失真。鑒于此，本文建立了大范圍、細(xì)網(wǎng)格波浪折射、繞射數(shù)學(xué)模型。

1 推廣的緩坡方程

本文在洪廣文推導(dǎo)的緩坡方程基礎(chǔ)上建立了模型［7］。洪廣文在傳統(tǒng)的緩坡方程中引入波數(shù)矢無(wú)旋性假定，把定常緩坡方程轉(zhuǎn)化為波作用方程和光程函數(shù)方程。這樣就把求解隨空間快速變化的波動(dòng)勢(shì)函數(shù)轉(zhuǎn)化為變化較慢的波作用和波數(shù)，使得求解緩坡方程時(shí)，空間步長(zhǎng)只取決于海底地形變化的劇烈程度，而不受波長(zhǎng)條件限制，增大了數(shù)學(xué)模型的計(jì)算范圍。

彌散方程

式中：ω為絕對(duì)頻率；g為重力加速度；k為波數(shù)；h為水深；W*為底摩阻能量損耗項(xiàng)。

波數(shù)矢無(wú)旋性方程

式中：K為波數(shù)矢；α為波浪傳播方向。

波作用守恒方程

光程函數(shù)方程

2 數(shù)值求解

波數(shù)矢無(wú)旋性方程和波作用守恒方程均屬于非線性雙曲型偏微分方程，需要選取一種有限差分格式進(jìn)行數(shù)值求解。本文采用后差分半隱格式，運(yùn)用這種差分格式編寫源代碼時(shí)，可以把方程對(duì)內(nèi)存的需求轉(zhuǎn)換為對(duì)硬盤的需求，從而大幅度提高計(jì)算范圍和網(wǎng)格分辨率。本文稱此方法為硬盤同步記錄法，適合計(jì)算像輻射沙脊群這種范圍大、網(wǎng)格分辨率要求高的海域。

2.1 硬盤同步記錄法

篇幅所限，這里僅以式（2）為例進(jìn)行有限差分，說(shuō)明硬盤同步記錄法的實(shí)現(xiàn)過(guò)程，其他方程解法類似。運(yùn)用后差分半隱式法，對(duì)式（2）進(jìn)行離散可得

式中：K、αi-1，j、Δx、Δy 為已知量；αi，j-1、αi，j、αi，j+1為未知量。式（5）屬于非線性雙曲型方程，用迭代法對(duì)其求解。由離散格式可以看出，只要知道 αi，j就可以求出 αi+1，j，αi-1，j的值不參與計(jì)算。這樣可以把 αi-1，j記錄到硬盤，同時(shí)釋放內(nèi)存。按以上步驟編程就可以用硬盤代替內(nèi)存記錄數(shù)據(jù)，從而實(shí)現(xiàn)硬盤同步記錄法，達(dá)到用高分辨率網(wǎng)格計(jì)算大范圍波浪場(chǎng)的目的。

2.2 初始條件和邊界條件

直接給定初始條件為外海原始波高H0、波周期T和波向α0。波高邊界條件

波向邊界條件

波浪破碎指標(biāo)為0.78H。

3 數(shù)學(xué)模型的驗(yàn)證

為了檢驗(yàn)本數(shù)學(xué)模型模擬波浪傳播變形的精確程度，采用Whalin［8］的物理模型實(shí)驗(yàn)值與數(shù)學(xué)模型計(jì)算值進(jìn)行比較。實(shí)驗(yàn)用的水槽長(zhǎng)25.5 m，寬6.096 m，水深分布如圖1所示，且可以用下式表示

式中：G（y）=［y（6.096-y）］1/2，0≤y≤6.096。

試驗(yàn)條件為波周期1 s，波高0.019 m。波浪入射方向?yàn)閤向，其空間步長(zhǎng)為0.05 m，y方向空間步長(zhǎng)為0.122 m，側(cè)邊界為全反射。在x軸5.8 m、17.0 m、18.76 m、20.15 m處，設(shè)置了4個(gè)測(cè)量斷面，用以獲得波高信息。圖2表示4個(gè)斷面的能量分布，由圖2可見數(shù)學(xué)模型與物理模型實(shí)驗(yàn)值能較好地吻合，模型的精確程度令人滿意。

4 數(shù)學(xué)模型在輻射沙脊群海域的應(yīng)用

4.1 數(shù)學(xué)模型在輻射沙脊群的適用性

輻射沙脊群海域范圍廣闊、地形復(fù)雜。它是呈輻射狀分布的出露于海面以上的沙洲與隱伏于海面以下的沙脊及其潮流水道的總稱。其地理位置處于江蘇岸外，黃海南陸架海域。北起射陽(yáng)河口，南至長(zhǎng)江口北部的蒿枝港，南北長(zhǎng)200 km，東西寬140 km。其水下地形十分復(fù)雜，大致以 港為頂點(diǎn)呈褶扇狀向海輻射（圖3）。由70多條沙脊和潮流水道組成，脊槽相間分布，水深介于0～25 m。分割沙脊的潮流水道主要有西洋、小北槽、陳家塢槽、草米樹洋、苦水洋、黃沙洋、爛沙洋、小廟洪等。這些大型的潮流水道，水深均大于10 m，深度向海遞增。

圖3用25個(gè)標(biāo)號(hào)點(diǎn)表示出8條主要潮流水道。同時(shí)為了方便討論各個(gè)潮流水道內(nèi)波浪場(chǎng)分布情況，把輻射沙脊群海域分成3塊區(qū)域。

計(jì)算輻射沙脊群波浪場(chǎng)所需網(wǎng)格數(shù)巨大，所以有必要使用硬盤同步記錄法。首先輻射沙脊群范圍廣闊，使得計(jì)算域面積大。其南北范圍介于 32°00＇N，33°48＇N，長(zhǎng)達(dá)199.6 km，東西范圍介于 120°40＇E～122°10＇E，寬度為140 km。其次水下地形變化劇烈，對(duì)網(wǎng)格分辨率要求高。在地形變化劇烈的海域，實(shí)測(cè)點(diǎn)距離小于30 m才能準(zhǔn)確反映出地形的起伏。因此如果網(wǎng)格步長(zhǎng)大于30 m則不能精確反映地形的變化。

綜上所述，由于輻射沙脊群范圍大、地形變化劇烈，為了準(zhǔn)確計(jì)算波浪場(chǎng)，需要按照x×y=140 km×200 km，dx×dy=30 m×30 m確定模型。這使得計(jì)算所需要的網(wǎng)格數(shù)巨大，超過(guò)了32位Fortran編譯器能夠編譯的范圍。如果采用硬盤同步記錄法則可以解決此問(wèn)題。

圖3 輻射沙脊群地形圖Fig.3 Relief map of radiative sandy ridge field

4.2 輻射沙脊群海域波浪場(chǎng)特征分析

應(yīng)用上述數(shù)學(xué)模型計(jì)算了輻射沙脊群海域波浪場(chǎng)。根據(jù)對(duì)此海域波況的觀測(cè)統(tǒng)計(jì)，本海區(qū)常年多出現(xiàn)NE向浪。限于篇幅，這里只給出NE方向平均潮位下，平均情況下波浪場(chǎng)的情況。

平均情況下的波浪場(chǎng)由于其出現(xiàn)頻率高，所以對(duì)輻射沙脊群海底地形的演變和沖淤平衡具有重要分析價(jià)值。本文取平均年最大風(fēng)速對(duì)應(yīng)的有效波高2.64 m［9］為平均情況下的波浪場(chǎng)計(jì)算原始條件。這里的平均年最大風(fēng)速相當(dāng)于2 a一遇的風(fēng)速。由于平均潮位歷時(shí)長(zhǎng)，波浪在此條件下對(duì)輻射沙脊群地形的作用時(shí)間也長(zhǎng)，因此對(duì)地形的演變和沖淤平衡具有重要分析價(jià)值。下面將以T=8 s，Ho=2.64 m作為分析探討沙脊群海域波浪場(chǎng)特征的原始條件，對(duì)平均潮位時(shí)波浪場(chǎng)各個(gè)特征進(jìn)行分析。

4.2.1 等波高線呈弧形分布

圖4實(shí)線表示1 m等波高線，虛線表示2 m等波高線。由圖4可以看出1 m等波高線環(huán)繞輻射沙脊群頂點(diǎn)弶港，呈弧形分布。在岸線與1 m等波高線之間形成范圍較大的、波高較小的掩蔽區(qū)。2 m等波高線環(huán)繞輻射沙脊群頂點(diǎn)弶港，呈鋸齒形分布。鋸齒形表示此片海域波高多為2 m。

圖4 輻射沙脊群海域等波高線分布Fig.4 Isoline of wave height in radiative sandy ridge field

為了清晰地表示出潮汐通道內(nèi)波高分布情況，在沙脊群間8條主要潮汐通道深槽中布置計(jì)算點(diǎn)25處，并將計(jì)算點(diǎn)連成等波高線（圖5）。其中細(xì)實(shí)線表示1 m等 波高線，粗實(shí)線表示2 m等波高線。由圖5可以看出在岸線與1 m等波高線之間形成范圍較大、波高較小的掩蔽區(qū)。這種獨(dú)特的波高分布狀態(tài)的形成是由于等水深線同樣大致呈弧形分布。弶港東側(cè)沙脊群頂點(diǎn)附近海域水深較小，波浪由外海傳來(lái)，由于多次破碎，波能消耗，波高減小。而弶港南北兩側(cè)水深較大，波浪傳至近岸海域，才會(huì)發(fā)生破碎，波高減小。因而形成了等波高線弧形分布的特征。

4.2.2 波浪多次破碎

波浪處于破碎狀態(tài)時(shí)，水體的紊動(dòng)與漩渦非常強(qiáng)烈。破碎帶是波能的主要消耗區(qū)域，也是海灘上泥沙運(yùn)動(dòng)最劇烈的地方。因此研究破波帶對(duì)輻射沙脊群的演變和沖淤平衡具有重要的分析價(jià)值。

圖6中近岸黑色區(qū)域是波浪的破碎點(diǎn)，灰色輪廓線是破波帶的外緣。由圖6可見，弶港東側(cè)海域破波帶外緣比其他海域離岸線要遠(yuǎn)。這主要是由于水下存在巨大的沙脊，使得該處水深較小。同時(shí)此海域?qū)儆谳椌蹍^(qū)，能量聚集使得波高增大。當(dāng)波高增大到一定程度，就產(chǎn)生了破碎現(xiàn)象，因此該海域破碎帶距離岸線比較遠(yuǎn)。

圖6 輻射沙脊群海域破波帶分布Fig.6 Surf zone of radiative sandy ridge field

隨著波浪傳播到近岸，又多次破碎，能量迅速減小，波高逐漸降低，在近岸海域形成范圍較大的、波高較小的掩蔽區(qū)。所以波浪多次破碎是輻射沙脊群響應(yīng)波浪場(chǎng)的重要特征之一。

4.2.3 以黃沙洋為界分為2個(gè)特征相異的波浪場(chǎng)

由圖 7可見，西洋（1～4號(hào)點(diǎn)）、陳家塢槽（8～11號(hào)點(diǎn)）、草米樹洋（12、13號(hào)點(diǎn)）、苦水洋（14～16號(hào)點(diǎn)）、黃沙洋（17～19號(hào)點(diǎn)）的波向隨著波浪向岸傳播逐漸轉(zhuǎn)向弶港。這主要是因?yàn)樵撈Ｓ虻壬罹€總體上呈以弶港為圓心的弧形分布（圖8）。其中黑色曲線表示概化的等水深線。由于地形對(duì)波浪的折射作用，波向線將集中，從而產(chǎn)生輻聚現(xiàn)象。波向線沿垂直于等深線方向傳播，于是形成了波向線以弶港為頂點(diǎn)呈輻射狀分布的現(xiàn)象。而爛沙洋（20～22號(hào)點(diǎn)）、小廟洪（23、24號(hào)點(diǎn)）海域波向并沒有明顯轉(zhuǎn)向弶港，這同樣是因?yàn)樗畹恼凵渥饔?。波浪受地形的折射影響，沿著與等深線垂直的方向傳播，而該海域等深線較為平順，總體上與岸線走向一致，所以波向線表現(xiàn)出垂直于岸線傳播的特征。

由以上分析可以看出地形對(duì)波向影響明顯。在黃沙洋以北地區(qū)，地形呈明顯的輻射狀，波浪受其影響轉(zhuǎn)向明顯。而在黃沙洋以南地區(qū)，由于深水區(qū)的入射波浪從NE方向入射，波峰線與等深線大體平行，同時(shí)這些位置水深較大，所以波向偏轉(zhuǎn)較小。

5 結(jié)論

黃海輻射沙脊群海域范圍廣闊，水下地形復(fù)雜，本文提出的波浪數(shù)學(xué)模型同時(shí)兼顧了大范圍和細(xì)網(wǎng)格的要求。模型通過(guò)硬盤同步記錄法，把數(shù)值求解時(shí)所需要的內(nèi)存用量轉(zhuǎn)化為硬盤用量，同時(shí)大幅度提高了計(jì)算范圍與網(wǎng)格分辨率。與N-S方程、Boussinesq方程、緩坡方程相比，本模型能夠計(jì)算大范圍波浪場(chǎng)；與動(dòng)譜平衡方程相比，本模型能夠用細(xì)網(wǎng)格計(jì)算大范圍波浪場(chǎng)。同時(shí)為了反映劇烈變化的地形對(duì)波浪場(chǎng)的影響，還引入了高階地形變化影響因子，使模型計(jì)算精度進(jìn)一步提高。綜上所述，本模型適合計(jì)算超大范圍下復(fù)雜地形波浪場(chǎng)。

模型在輻射沙脊群計(jì)算結(jié)果表明：

（1）在平均潮位下，外圍深水波浪傳播進(jìn)入沙脊群海域，由于水下地形阻礙作用，等波高線環(huán)繞輻射沙脊群頂點(diǎn)弶港呈弧形分布。

（2）在黃沙洋以北的潮汐通道中，波向變化顯著，大體指向弶港方向，而黃沙洋以南波峰線與岸線基本平行。

（3）隨著波浪進(jìn)入沙脊群頂點(diǎn)附近的淺水區(qū)，多次發(fā)生破碎，波能不斷消耗，波高逐漸減小，在近岸形成了范圍較大的掩蔽區(qū)。

［1］劉海成，劉海源，楊會(huì)利.近岸波浪變形數(shù)值模型的比較研究［J］.水道港口，2009，30（3）：153-158.LIU H C，LIU H Y，YANG H L.Comparison of numerical models of nearshore wave transformation［J］.Journal of waterway and harbor，2009，30（3）：153-158.

［2］ZHAO L H，LI T C.Fractional-step finite element method for calculation of 3-D free surface problem using level set method［J］.Journal of Hydrodynamics：Series B，2006，18（6）：742-747.

［3］WU C H，YUAN H L.Efficient non-hydrostatic modeling of surface waves interacting with structures［J］.Applied Mathematical Modelling，2007，31（4）：687-699.

［4］FUHRMAN D R，MADSEN P A.Simulation of nonlinear wave run-up with a high-order Boussinesq model［J］.Coastal Engineering，2008，55：139-154.

［5］JIN H，ZOU Z L.Hyperbolic mild slope equations with inclusion of amplitude dispersion effect：regular waves［J］.China Ocean Engineering，2008，22（3）：431-444.

［6］XU F M，Will P，ZHANG J L，et al.Simulation of typhoon-driven waves in the Yangtze estuary with multiple-nested wave models［J］.China Ocean Engineering，2005，19（4）：613-624.

［7］HONG G W.Mathematical models for combined refraction-diffraction of waves on non-uniform current and depth［J］.China Ocean Engineering，1996，10（4）：433-454.

［8］Whalin R W.Wave refraction theory in a convergence zone［C］//Coastal Engineering Research Council.Proc.13th Int.Conf.Coastal Engineering.Vancouver：ASCE，1972：451-470.

［9］張東生，張君倫，張長(zhǎng)寬，等.潮流塑造-風(fēng)暴破壞-潮流恢復(fù)-試釋黃海海底輻射沙脊群形成演變的動(dòng)力機(jī)制［J］.中國(guó)科學(xué)：D輯，1998，28（5）：394-402.ZHANG D S，ZHANG J L，ZHANG C K，et al.Tidal current work-surge destruction-tidal current restoration：A preliminary study on the dynamic mechanism of the development of the radical submarine sand ridges of the South Yellow Sea［J］.Science in China：Series D，1998，28（5）：394-402.