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

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

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

1 推廣的緩坡方程

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

彌散方程

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

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

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

波作用守恒方程

光程函數(shù)方程

2 數(shù)值求解

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

2.1 硬盤同步記錄法

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

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

2.2 初始條件和邊界條件

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

波向邊界條件

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

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

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

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

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

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

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

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

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

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

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

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

4.2 輻射沙脊群海域波浪場特征分析

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

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

4.2.1 等波高線呈弧形分布

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

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

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

4.2.2 波浪多次破碎

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

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

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

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

4.2.3 以黃沙洋為界分為2個特征相異的波浪場

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

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

5 結(jié)論

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

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

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

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

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

［1］劉海成，劉海源，楊會利.近岸波浪變形數(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］張東生，張君倫，張長寬，等.潮流塑造-風(fēng)暴破壞-潮流恢復(fù)-試釋黃海海底輻射沙脊群形成演變的動力機制［J］.中國科學(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.