劉文帥 ，鄭建國* ，張翀岳，金智濤，許國輝

(1.中國海洋大學 海洋環(huán)境與生態(tài)教育部重點實驗室，山東 青島266100；2.中國海洋大學 山東省海洋環(huán)境地質工程重點實驗室，山東 青島266100；3.中國海洋大學 環(huán)境科學與工程學院，山東 青島266100)

引 言

波浪由外海傳至近岸，隨水深變淺波高逐漸增大，無法維持較大的波陡時波浪就會發(fā)生破碎。波浪破碎引起強烈的水體紊動，能量損耗大。波浪破碎后會以新的波高繼續(xù)向前傳播，甚至可能發(fā)生多次破碎，直至能量消耗殆盡。波浪破碎涉及到劇烈的能量和動量交換，波浪能量的衰減及近岸水體運輸都與這種交換作用有關。波浪的破碎易造成破碎區(qū)結構物的破壞，航道和港區(qū)的淤積等問題[1]，因此針對波浪破碎相關問題進行研究，能夠更好地指導實際工程，也可以為岸線侵蝕的評估等實際問題提供一定參照。

目前許多國內外學者針對波浪破碎相關的能量問題作了大量理論模擬和試驗探索，得到一些寶貴的實測資料和經驗公式,為求解破后波、近岸區(qū)增減水等一系列問題提供了依據。鄒志利[2]總結前人提出的破碎波能量損失波能方程，指出波能損耗率是計算關鍵，通過將計算波高與實測波高吻合，試圖找到符合實際的波能損耗率。許多研究都在此基礎上尋找更好的波能損耗率表達式，主要包括涌模型和水滾模型兩種。

涌模型可描述波浪破碎過程及其能量的變化情況，破碎過程可分波浪變形劇烈的外區(qū)和滾動旋渦消耗能量為主的內區(qū)[3-4]，利用涌模型進行能量損失計算對波浪崩破及破碎內區(qū)有一定的模擬效果[5]，但涌模型無法解釋破碎機理；而水滾模型是將波峰附近表層形成的混亂水體稱為水滾，水滾的水平速度和波浪的相速度相同[6]，水滾和下面波浪之間的剪力摩擦是破碎能量損失的主要因素，單位面積的能量損失率與水滾的剖面面積有關[7]，不同的剖面面積可有不同的算法[4,8]，但眾多參數(shù)之間的相互關系仍有爭議。除上述兩模型外，也有學者通過其他模型研究波浪破碎過程中的能量問題，主要包括Boussinesq方程[9-11]和緩坡方程[12-14]。

波浪破碎后，水體狀態(tài)與紊流相似，因此可利用紊流模型來研究波浪破碎造成的能量損失問題。HORIKAWA等[15]根據各項同性紊流及孤立波理論，計算了破波能量損失問題，SAWARAGI等[16]指出波浪破碎形成的紊亂狀態(tài)不是各向同性的，且Horikawa的計算中未能考慮到增水的影響，Mizuchi[17]進一步優(yōu)化了破碎帶內的能量方程。

盡管許多學者對于波浪能量等問題提出了很多見解，但大都未考慮過波浪淺化、破碎過程中灘面形態(tài)的影響。波浪破碎造成水體紊動和能量損失，并伴隨著灘形變化，平直灘形最終發(fā)展形成沙紋灘形，不同灘形在同樣波況下的紊動能量分布和破波能損可能是不同的。本文基于FLOW-3D軟件，建立并驗證了數(shù)值波浪水槽內的沙灘沖刷模型，通過固定歐拉網格求解Navier-Stokes方程，得到了平直灘形和沙紋灘形在多種工況卷破波作用下的紊動能量分布，并對比了多種工況下的波浪破碎能耗。

1 室內水槽試驗

試驗所用波浪水槽長寬高尺寸分別為14、0.5、0.6 m，水槽右端配有推板造波系統(tǒng)，靠近左端鋪設試驗用平直沙灘模型(試驗用砂為取自青島周邊海灘篩分得到的中砂，D50=0.375 mm，各組分質量分數(shù)：0～0.25 mm占2.8%，0.25～0.5 mm占81.9%，0.5～1 mm占15.3%)，水槽兩端均配備消浪設施以保證波浪穩(wěn)定性。在沙灘坡腳設置1個波高儀，記錄波高變化，具體試驗設置見圖1。試驗主要內容為，觀察特定波浪向左傳播并在砂質岸灘上破碎過程中，破波點附近的波浪形態(tài)以及破碎對于岸灘造成的影響。在試驗開始前和結束后，通過地形儀分別測得平直灘形和沙紋灘形(圖2)。試驗工況見表1，所有工況波浪破碎類型均為卷破波。

圖1 試驗設置示意圖

表1 工況設置表

圖2 各工況灘形

2 模型建立與驗證

2.1 理論模型

FLOW-3D計算性能高，可以自定義多種物理場，廣泛應用于不同領域，能夠精準預測自由液面流動。FLOW-3D通過其獨特的FAVORTM技術，將其各理論模型表達為體積和面積、孔隙度的函數(shù)，在模擬復雜幾何區(qū)域時具有一定優(yōu)勢。本文所用計算模型主要包括連續(xù)方程、動量方程、流體自由面方程(VOF方法)及泥沙沖刷模型，紊流模型選取RNGk-ε紊流模型(求解速度快,精度高[18])。

2.1.1 連續(xù)性方程

(1)

2.1.2 動量方程

(2)

(3)

(4)

式中，p是壓力；VF是流體所占體積分數(shù)，由FLOW-3D建模時確定；ρ為流體密度；u、v、w為在t時刻對應坐標上流速分量；Ax、Ay、Az為FLOW-3D在t時刻3個方向FAVORTM面積函數(shù)；Gx、Gy、Gz為加速度分量；fx、fy、fz為3個方向的粘性力。

2.1.3 流體自由面方程

流體自由面捕獲采用流體體積法(VOF)，其方程如下：

(5)

(6)

式中,VF為擴散系數(shù)。

2.1.4 泥沙沖刷模型

FLOW-3D中的泥沙模型可用于計算泥沙的輸運、沖刷和沉降，采用VOF方法模擬流體和泥沙的分界面[19]，泥沙以底沙(堆積在河床上的泥沙)和懸沙(流體中的泥沙)兩種狀態(tài)存在，懸沙濃度通過對流擴散方程控制。泥沙模型通過追蹤底沙和懸沙的含沙濃度(Cp，Cs)，計算泥沙體積分數(shù)?？偰嗌丑w積分數(shù)αs是網格中懸沙和底沙的總體積分數(shù)，表達式為：

(7)

ρ=αfρf+(1-αf)ρs

(8)

式中，ρf為流體密度；ρs為泥沙密度；αf為網格中流體的體積分數(shù)。用當?shù)孛芏忍鎿QNavier-Stokes方程中的流體密度，根據總含沙量在每個網格中計算出一對一的當?shù)孛芏?，粘性增量或者阻力增量引起的泥沙體積分數(shù)取決于αc0。

2.2 模型建立

采用FLOW-3D軟件建立三維沖刷模型(寬度方向為單一網格)，F(xiàn)LOW-3D可以通過在固定歐拉網格中采用有限體積差分法求解三維瞬態(tài)Navier-Stokes方程。使用嵌套網格對灘形附近和自由水面波動范圍內進行局部加密。數(shù)值波浪水槽一端(Xmin)為造波邊界，為了消除回波反射的影響，在數(shù)值水槽的另一端(Xmax)設置多孔介質消波裝置配合出流邊界消除反射，具體設置及參數(shù)主要參考張婷[20]的相關設置，其他邊界條件為Ymin,Ymax，Zmax：對稱邊界；Zmin：墻邊界。計算時長為30 s，步長為0.05 s，網格總數(shù)為50萬個，Z方向和X方向網格最小尺寸分別為0.002 75、0.004 m。

圖3 網格劃分圖(局部)

2.3 模型驗證

為驗證數(shù)值波浪水槽的可靠性，在波形穩(wěn)定階段(15～25 s)從平直灘形坡腳自由水面時序曲線(圖4)、波浪破碎形態(tài)(圖5，平直灘形與沙紋灘形各選一例用于簡要說明)兩個方面進行驗證。

(a) H=6.7 cm (b) H=8.6 cm (c) H=13.2 cm

圖5 波浪破碎形態(tài)驗證圖(部分)

驗證結果表明，穩(wěn)定后自由水面時序曲線基本吻合(波峰、波谷最大誤差分別不超過3.8%和4.2%)，波浪破碎形態(tài)模擬與實驗對應關系良好,軟件對于波浪在岸灘上的變形與運動具有良好的描述能力，能夠較為準確地捕捉岸灘附近的水體，可用于進一步研究岸灘上的波浪破碎相關問題。

3 模擬結果

使用FLOW-3D提取了模擬結果中波形穩(wěn)定時段(15～25 s)平直灘形與沙紋灘形波浪破碎前后水體自由水面時序曲線(圖6)。

圖6 破碎前后自由水面時序曲線圖

計算結果表明，波浪破碎前，由于水深不斷變淺，波浪表現(xiàn)出一定的非線性特征，隨水深進一步變淺，波浪發(fā)生破碎，破碎后波浪非線性明顯增強，且波高變化明顯，波形與波高的變化表明了破碎造成波浪能量的變化，關于能量變化的具體分析詳見下節(jié)。

4 分析與討論

4.1 紊動能量

卷破波造成水體紊亂的基本模式為：波浪卷破前形成水舌，在水舌形成后至水舌濺落至水體前過程中，除底層水體因與表層沙灘相互作用造成較弱的水體紊動能量外，水舌表層也存在較弱的紊動能量(圖7a,7b)。水舌形成后，由于重力、慣性的作用，向前方濺落，水舌頭部濺落后迅速提升水體紊動程度(圖7c,7d)，在波浪破碎的全過程中此時水體的紊動能量達到最大值，以濺落點為中心，紊動能量向附近擴散，周圍水體也逐漸受到紊動影響。當波高較大時，波浪往往會發(fā)生多次破碎，其模式基本與前述大致相同，都在水舌頭部濺落時達到紊動能量高值，并以濺落點為中心向周圍水體擴散，直到水體只能形成微弱起伏，無法形成水舌，也無法再次破碎，因此后方水體基本不產生較大紊動能量。

圖7 水舌形成、濺落圖

本文統(tǒng)計了所有工況中的最大紊動能量值(圖8)，發(fā)現(xiàn)最大紊動能量值與波高呈明顯的正向相關性，這是因為波高增大時，破碎形成的水舌也隨之變大，水舌濺落后在水體表面造成的紊動更強。而相同波況作用下，沙紋灘形水面最大紊動強度要小于平直灘形，在波高較大時，這一現(xiàn)象更加明顯，這表明沙紋灘形會影響卷破波的破碎過程，尤其是水舌的形成和濺落過程，從而影響卷破波的紊動能量的大小與分布。所以，在分析討論實際海灘對于卷破波紊動能量的大小和分布的影響時，不應當將實際灘形簡化為表面平直具有一定坡度的斜坡，而是要充分考慮到凹凸不平的灘面形態(tài)對于水體流場的影響和作用。

圖8 最大紊動能量統(tǒng)計圖

4.2 破波能損

將破波前后水體的自由水面時序曲線通過快速傅里葉變換將數(shù)據從時域轉換為頻域，分離出各個子波，通過白志剛等[21]提出的波能計算方法分別求解各子波的波能后在一定時間內(取1個周期，即1.68 s)疊加，得到波浪破碎前后的波能，并計算得到波浪的能量破損率(圖9)。

圖9 各工況波浪能量破損率圖

平直灘形，1:8坡度隨著波高的增加,破波能損率逐漸降低；而1:10坡度的破波能損率呈相反趨勢，隨波高增加而逐漸提高。且在波高較小(H=6.7,8.6 cm)時，1:8坡度破波能損率要高于1:10坡度；而波高較大(H=13.2 cm)時，1:8坡度破波能損率低于1:10坡度。

沙紋灘形破波能損率變化較為復雜，伴隨著波高的提高，破波能損率會出現(xiàn)陡增。其原因應與沙壩在破波區(qū)域中是否存在具有一定的關聯(lián)性。1:8沙紋灘形在波高較小(H=6.7 cm)時，破波能損率較低，而波高增大至H=8.6,13.2 cm后，破波能損率分別提升12.29%和15.99%，而對于相同坡度波高較大的工況(H=8.6,13.2 cm)波高增大后破波能損率變化不大(3.7%);1:10灘形情況相似，在H=6.7,8.6 cm工況，破波能損率較低，分別為55.17%和48.38%，當波高增至H=13.2 cm時，破波能損率均有大幅提升，分別提升15.82%和22.61%。在所有破波能損率明顯增大的工況中(1:8 H=8.6,13.2 cm；1:10 H=13.2 cm)，均在破波區(qū)域內出現(xiàn)明顯沙壩，而破波能損率較小的工況(1:8 H=6.7 cm；1:10 H=6.7,8.6 cm) 則未出現(xiàn)明顯沙壩。沙壩的存在，主要是由于波浪在沙灘上破碎后造成大量泥沙起動，在回落水流等的作用下沉積而最終形成。波浪傳播至破波線附近時，由于水深變淺，波浪開始發(fā)生破碎，水舌初步形成，水舌向前傳播過程中體積不斷增大，當遭遇水下沙壩，水舌會越過沙壩向前運動，在水舌跨越沙壩過程中，由于水深急劇變淺，波浪難以維持水舌并進一步向前傳播，在沙壩后方破碎，沙壩對于波浪破碎的作用與水下潛堤的效果相似。

此外，1:10坡度H=13.2 cm波況破波能損率除了受到沙壩的較大影響外，造成這一波況下破波能損率遠高于同一坡度沙紋灘形其他兩種波況的另外一個重要原因在于，水舌墜落后，由于沙灘地形的原因二次彈起的水舌分為兩個部分，一部分水舌與其他工況相同，向前繼續(xù)傳播，而另一部分則是向反方向傳播，造成了向前傳播的水舌體積小于正常情況，較小體積的水舌在濺落時無法引起前方水體的較大幅度波動，因此前方水體波浪能量較小(圖10)。

圖10 波浪破碎模擬圖 (1:10 H=13.2 cm)

5 結論

本文基于FLOW-3D軟件，采用VOF方法描述水體的自由表面，建立并驗證了三維岸灘波浪沖刷模型，探討了卷破波在平直灘形和沙紋灘形上破碎的紊動能量分布與波浪破碎能耗，得到以下結論：

(1)波浪時序曲線、波浪破碎形態(tài)與實驗結果基本吻合，表明本模型的合理性與可靠性。

(2)波浪卷破后，水體紊動能量在水面附近達到最大值，在深度方向上逐漸減小，沙紋灘形的紊動能量最大值小于同波況下的平直灘形，沙灘坡度與灘形會影響紊動能量在水體中的分布。在研究波浪破碎造成的紊動能量分布相關問題時，不能將灘形簡化為表面平直具有一定坡度的斜面。

(3)平直灘形卷破波的破波能損與波高的關系受坡度影響，1:8平直灘形破波能損隨波高增大而逐漸減小，1:10平直灘形破波能損隨波高增大而逐漸增大；沙紋灘形卷破波的破波能損與灘形在破波區(qū)域內是否發(fā)育出現(xiàn)沙壩有較大關系，沙壩的出現(xiàn)提高了卷破波的能量耗損率。灘形的起伏通過阻礙水舌彈起后的向前傳播或改變部分水舌的傳播方向，一定程度上改變了破波能損率。