幕墻式防波堤水動力特性的數值模擬研究

尹勝強，劉 勇，王心玉

(中國海洋大學 海岸與海洋工程研究所，青島 266100)

圖1 幕墻式防波堤結構示意圖Fig.1 Sketch of a curtain breakwater

幕墻式防波堤是一種新型防波堤，由直立板防波堤逐漸發展而來。如圖1所示，幕墻式防波堤由前下垂板和后直墻兩部分組成，下垂板與直墻之間形成消浪室。消浪室促使下垂板周圍的渦旋明顯增大，入射波能量耗散增加。與傳統防波堤相比，幕墻式防波堤可以有效降低波浪反射，保證結構前水域平穩。

許多學者對直立板結構的水動力特性進行了研究。Ursell[1]最早研究了無限水深中波浪與直立板結構的相互作用問題，給出了反射系數和透射系數的變化規律。Evans[2]研究了波浪與淹沒式直立板結構的相互作用，得到了一階波浪力與二階波浪力的計算公式。Morris[3]建立了一種新的變分方法，并應用于分析波浪與非對稱直立板結構的相互作用問題。Losada等[4]和Abul-Azm[5]采用匹配特征函數展開法分別研究了正向波和斜向波與直立板結構的相互作用問題，給出了反射系數和透射系數的變化規律。Kriebel和Bollmann[6]利用能量透射理論和匹配特征函數展開法研究了波浪與出水實體薄板結構的相互作用，發現兩種方法計算得到的反射系數基本一致。陳雪峰等[7]研究了波浪與開孔板結構的相互作用，給出了開孔率等因素對開孔板點壓力分布和反射系數的影響規律。

在實際工程中，為了提高掩護效果，經常采用雙排或多排直立板結構。與單排直立板結構相比，雙排或多排直立板結構可以更有效地反射入射波能量。Porter[8]采用多項伽遼金方法研究了波浪與雙排實體板的相互作用，給出了反射系數和透射系數的精確計算結果。Liu等[9]結合理論分析和物理模型試驗，研究了波浪對Jarlan型透空潛堤(前開孔板和后實體板組成)的作用，理論計算結果與試驗結果符合良好，當相對板間距為 0.3～0.4時，潛堤掩護效果最佳。Clauss等[10]結合物理模型試驗和理論分析，研究了多排垂直開孔板結構的水動力特性。

Nakamura等[11-12]在研究雙排直立板結構水動力特性的基礎上，采用直墻代替后側垂直板，提出了幕墻式防波堤的概念，并通過物理模型試驗研究了消浪室寬度、下垂板吃水深度對幕墻式防波堤反射系數的影響。此后，Nakamura團隊針對幕墻式防波堤及其改進結構開展了大量物理模型試驗，深入研究了幕墻式防波堤的水動力特性。Nishikawa等[13]對幕墻式防波堤的受力特性進行了理論和試驗研究，分析了入射波高和消浪室相對寬度對下垂板與直墻受力的影響。Ono等[14]采用數值模擬結合物理模型試驗的方法，分析了波浪周期對幕墻式防波堤周圍流場特性的影響。耿寶磊等[15]研究了波浪與直墻前多層透空薄板結構的相互作用，得到了反射系數、透射系數與波能吸收率的變化規律。

綜上所述，現有關于幕墻式防波堤水動力特性的研究以理論分析和物理模型試驗為主，關于幕墻式防波堤水動力特性的數值模擬研究較少。與理論分析和物理模型試驗相比，數值模擬可以有效分析幕墻式防波堤周圍的流場特性，從而深入理解幕墻式防波堤的消浪機理。為此，本文基于CFD開源程序庫OpenFOAM建立二維數值波浪水槽，模擬波浪對幕墻式防波堤的作用，研究波長、消浪室寬度、下垂板吃水深度和下垂板底角對幕墻式防波堤反射系數、波浪力和流場特性的影響，分析幕墻式防波堤的消浪機理，研究結果可為工程設計提供科學指導。

1 基于OpenFOAM的數值模型

1.1 基本控制方程

將流體假設為不可壓縮粘性流體，在計算過程中滿足連續性方程和動量方程。

連續性方程

(1)

動量方程

(2)

式中：ρ為流體密度；xi為空間點的坐標；ui為在t時刻坐標xi點的速度分量；p為壓力；μ為流體動力粘性系數；fi為源項；針對本文的二維問題而言，i、j的取值范圍為1、2，表示直角坐標系下，各個矢量在x、y方向上的分量。

1.2 自由面處理方法

采用VOF方法(流體體積法)捕捉自由面。定義體積分數函數α表示計算網格單元中液相所占體積分數。

體積分數函數α滿足連續性方程

(3)

式中：ur是為了實現對自由液面的壓縮而引入的只在自由面起作用的相對速度。由下式確定

(4)

式中：φ為通過自由面單元表面的體積流量；Sf為對應的單元表面面積；nf為自由面單位法向向量；Cα為自由面壓縮系數。

通過公式(3)計算得到體積分數α后，氣液兩相流場內各點的密度和粘性系數可以通過體積分數函數α表示為

ρ=αwaterρwater+(1+αwater)ρair

(5)

μ=αwaterμwater+(1+αwater)μair

(6)

1.3 數值造波方法

采用速度入口邊界造波方法造波，通過在入口邊界上設置波面和流體速度實現造波。本文二維數值波浪水槽采用二階斯托克斯波，二階斯托克斯波的波面和流體速度如下式所示

(7)

(8)

(9)

式中：η為波面變化；H為波高；k為波數；ω為波浪圓頻率；u為水平質點速度；v為垂直質點速度；h為水深。

1.4 數值消波方法

采用松弛區域法消波，在二維數值波浪水槽的入口區域和出口區域設置松弛區域，在松弛區域內每一時刻對速度進行如下修正

u=αRumodel+(1-αR)utarget

(10)

式中：umodel為通過求解N-S方程得到的速度值；utarget為期望得到的目標速度值；αR為與空間位置有關的加權函數，αR滿足下式

(11)

式中：XR∈[0,1]，在松弛區域邊界XR=1，αR=0；而在松弛區域與非松弛區域的交界XR= 0，αR=1。圖2給出入口和出口區域松弛區域內αR和XR的變化。

圖2 松弛區域內αR和XR的變化圖Fig.2 Variation of αR versus XR in relaxation zone

1.5 數值求解方法

采用有限體積法對控制方程進行離散求解，計算過程中將計算域劃分為若干網格單元，在每個網格單元內對控制方程進行積分，再通過高斯變化得到離散后的控制方程。

采用GAMG(Geometric-Algebraic Multi-Grid)方法求解離散后的線性代數方程組，先將計算域網格合并為少數粗網格，在粗網格上快速求解得到方程組的解，然后將粗網格上得到的方程組的解作為細網格上的預測值，在細網格上對線性代數方程組迭代求解，最終得到精確的數值解。

求解過程中采用同位網格存儲，速度、壓力等均存儲在計算網格中心處，計算得到網格中心處的流場變量值后，再通過插值方法獲得其他位置的流場信息。

圖3 二維數值波浪水槽Fig.3 Sketch of two-dimensional numerical tank

2 數值模型驗證

2.1 數值水槽設置

如圖3所示，二維數值波浪水槽長12 m，高0.8 m。幕墻式防波堤由前下垂板和后直墻兩部分組成，其中下垂板厚0.042 m，直墻厚0.15 m。試驗布置參考Nakamura[11]的物理模型試驗(比尺1:12)，水深42 cm，波高6 cm，波浪周期0.9～1.8 s；幕墻式防波堤距離二維數值水槽末端2 m，消浪室寬度B分別為21 cm、23 cm、25 cm、 27 cm、29 cm；下垂板吃水深度d分別為21 cm、23 cm、25 cm、27 cm、29 cm；下垂板底角θ分別為30°、45°、60°、90°、-45°(左下端)。波高儀放置的位置為：x= 5.2 m、x=5.4 m、x=5.6 m、x= 9 m、x=10.2 m。出口區域的消波區為率定入射波時使用，放入防波堤后不再設置。所有的數值試驗條件均列于表1。

表1 數值試驗條件Tab.1 Numerical test conditions

2.2 計算網格

本文建立的二維數值水槽采用結構化網格，以矩形網格為主，存在少量的不規則四邊形網格。經網格收斂性檢驗，最終選用圖4所示的網格劃分方式，計算網格縱橫比為1/4，沿水平方向計算網格尺寸不超過0.02 m, 沿豎直方向計算網格尺寸不超過0.005 m。為了更精確地模擬自由面與結構物周圍的流場變化，對自由面與結構物附近的網格做局部加密。網格數量保持在15萬左右。

圖4 計算網格劃分Fig.4 Sketch of grid generation

2.3 數值結果驗證

對于數值模擬而言，波浪在傳播過程中，隨著傳播距離的增大，由于流體粘性、數值截斷誤差以及邊界作用的影響波高會逐漸衰減。因此，在數值模擬試驗開始之前，需要在二維數值波浪水槽中進行率波。

選取典型波浪參數(T=0.9 s，H=6 cm；T=1.4 s，H= 6 cm；T=1.8 s，H=6 cm)下x=10 m處模擬值與理論值的對比結果進行說明。如圖5～圖7所示，點為基于二階斯托克斯波的數值模擬結果，線為二階斯托克斯波的理論結果。可以看出，數值模擬結果波形規則、穩定，說明二維數值波浪水槽可以有效消除造波邊界和水槽末端波浪反射的影響。數值模擬結果與理論結果存在很小的差異，這是由于在數值模擬過程中考慮了流體的粘性，存在耗散現象，而理論結果沒有考慮流體粘性。總體而言，數值模擬結果與理論結果符合良好。

圖5 波面歷時曲線(T=0.9 s,H=6 cm)Fig.5 Duration curve of wave surface(T=0.9 s and H= 6 cm) 圖6 波面歷時曲線(T= 1.4 s,H= 6 cm)Fig.6 Duration curve of wave surface(T= 1.4 s and H= 6 cm)圖7 波面歷時曲線(T=1.8 s,H= 6 cm)Fig.7 Duration curve of wave surface(T=1.8 s and H=6 cm)

將反射系數Cr和消浪室內相對波高Hc/H的數值模擬結果與Nakamura等[10]的物理模型試驗結果進行對比驗證。在計算反射系數與消浪室內相對波高時，以至少十個穩定波為準。圖8～圖10給出了三種不同工況下反射系數的數值模擬結果與Nakamura等[11]的物理模型試驗結果對比，圖中實心點表示數值模擬結果，空心點表示物理模型試驗結果。可以看出，數值模擬結果與物理模型試驗結果符合良好，此外，隨著波長的變化，反射系數存在最小值(約為0.1)，說明通過合理選取結構參數，幕墻式防波堤可以非常有效地耗散入射波能量。

圖8 反射系數的數值模擬結果與試驗結果[11]對比(B=21 cm,d=21 cm)Fig.8 Comparison of reflection coefficient between numerical results and experimental data [11](B=21 cm and d=21 cm)圖9 反射系數的數值模擬結果與試驗結果[11]對比(B=29 cm,d=21 cm)Fig.9 Comparison of reflection coefficient between numerical results and experimental data [11](B=29 cm and d=21 cm)

圖10 反射系數的數值模擬結果與試驗結果[11]對比(B=29 cm,d=29 cm)Fig.10 Comparison of reflection coefficient between numerical results and experimental data[11](B=29 cm and d=29 cm)圖11 消浪室內相對波高的數值模擬結果與試驗結果[11]對比(B=29 cm,d=29 cm)Fig.11 Comparison of relative wave height in the wave chamber between numerical results and experimental data[11](B=29 cm and d=29 cm)

圖11給出了消浪室內相對波高的數值模擬結果與Nakamura等[11]的物理模型試驗結果對比，圖中實心點表示數值模擬結果，空心點表示物理模型試驗結果。可以看出，數值模擬結果與物理模型試驗結果符合良好。

3 數值算例與討論

3.1 反射系數

3.1.1 消浪室寬度與下垂板吃水深度的影響

圖12給出了不同消浪室寬度下，反射系數Cr的數值計算結果隨下垂板相對吃水深度d/L的變化曲線。從圖12可以看出：對于長周期波而言，隨著消浪室寬度增大，波浪能量耗散增加，反射系數降低；但是對于短周期波而言，變化規律正好相反。圖13給出了不同下垂板吃水深度下，反射系數Cr的數值計算結果隨消浪室相對寬度B/L的變化曲線。從圖13可以看出：對于長周期波而言，隨著下垂板吃水深度增大，防波堤的消浪效果增強；對于短周期波而言，隨著下垂板吃水深度增大，防波堤的消浪效果減弱。此外，當消浪室相對寬度為0.08～0.10，下垂板相對吃水深度為0.09～0.12時，幕墻式防波堤消浪效果最好。從圖13還可以看出，隨著下垂板吃水深度的增大，反射系數最小值對應的消浪室相對寬度逐漸減小。

圖12 反射系數隨下垂板相對吃水深度的變化(d=21 cm)Fig.12 Variation of the reflection coefficient versus the relative immersed depth of the drooping plate(d=21 cm) 圖13 反射系數隨消浪室相對寬度的變化(B=29 cm)Fig.13 Variation of the reflection coefficient versus the relative wave chamber width(B=29 cm)

3.1.2 下垂板底角的影響

圖14 反射系數隨消浪室相對寬度的變化(B=29 cm, d=21 cm)Fig.14 Variation of the reflection coefficient versus the relative wave chamber width (B=29 cm and d=21 cm)

圖14給出了不同下垂板底角下，反射系數Cr的數值計算結果隨消浪室相對寬度B/L的變化曲線。從圖14可以看出：下垂板底角越尖銳，耗散的波浪能量越多，反射系數越小；與長周期波相比，短周期波條件下反射系數對于下垂板底角的變化更加敏感。此外，下垂板底角的方向(θ=45°或θ=-45°)對于反射系數的變化幾乎沒有影響。

3.2 波浪力

3.2.1 消浪室寬度的影響

圖15給出了不同消浪室寬度下，下垂板波浪力的數值計算結果隨下垂板相對吃水深度d/L的變化曲線。從圖15可以看出：對于長周期波而言，隨著消浪室寬度的增加，下垂板受力增加；對于短周期波而言，隨著消浪室寬度增加，下垂板受力減小。

圖16給出了不同消浪室寬度下，直墻波浪力的數值計算結果隨下垂板相對吃水深度d/L的變化曲線。從圖16可以看出：隨著消浪室寬度增加，直墻受力減少；隨著下垂板相對吃水深度增大(波長減小)，直墻受力減小。

圖15 下垂板受力隨下垂板相對吃水深度的變化(d=21 cm)Fig.15 Variation of wave force on drooping plate versus the relative immersed depth of the drooping plate (d=21 cm) 圖16 直墻受力隨下垂板相對吃水深度的變化 (d=21 cm)Fig.16 Variation of wave force on vertical wall versus the relative immersed depth of the drooping plate (d=21 cm)

3.2.2 下垂板吃水深度的影響

圖17給出了不同下垂板吃水深度下，下垂板波浪力的數值計算結果隨消浪室相對寬度B/L的變化曲線。從圖17可以看出：隨著下垂板吃水深度增加，下垂板受力明顯增大，這是因為隨著下垂板吃水深度增加，下垂板的受力面積增大；長周期波條件下，下垂板受力對下垂板吃水深度的變化更加敏感。另外，當消浪室相對寬度B/L=0.15時，下垂板受力最大，約為0.6。

圖18給出了不同下垂板吃水深度下，直墻波浪力的數值計算結果隨消浪室相對寬度B/L的變化曲線。從圖18可以看出：隨著下垂板吃水深度增加，直墻受力減小；隨著消浪室相對寬度增大(波長減小)，直墻受力減小。

圖17 下垂板受力隨消浪室相對寬度的變化(B=29 cm)Fig.17 Variation of wave force on drooping plate versus the relative wave chamber width (B=29 cm)圖18 直墻受力隨消浪室相對寬度的變化(B=29 cm)Fig.18 Variation of wave force on vertical wall versus the relative wave chamber width (B=29 cm)

3.2.3 下垂板底角的影響

圖19給出了不同下垂板底角下，下垂板波浪力的數值計算結果隨消浪室相對寬度B/L的變化曲線。從圖19可以看出：隨著下垂板底角的增大，下垂板受力增大；θ=45°(右下端)時下垂板受力明顯小于θ=-45°(左下端)。

圖20給出了不同下垂板底角下，直墻波浪力的數值計算結果隨消浪室相對寬度B/L的變化曲線。從圖20可以看出：隨著下垂板底角增大，直墻受力的變化減小，最大變化值約為0.1(變化率僅為3%)。因此在實際工程應用中，下垂板底角對直墻受力的影響基本可以忽略。

圖19 下垂板受力隨下垂板相對吃水深度的變化(B=29 cm,d=21 cm)Fig.19 Variation of wave force on drooping plate versus the relative wave chamber width(B=29 cm and d=21 cm)圖20 直墻受力隨消浪室相對寬度的變化(B=29 cm,d=21 cm)Fig.20 Variation of wave force on vertical wall versus the relative wave chamber width(B=29 cm and d=21 cm)

3.3 流場特性分析

本節以T=1.4 s，B=29 cm，d=21 cm，θ= 45°的工況為例，分析下垂板周圍的流場特性。圖21-a～21-e給出了一個波浪周期內典型時刻幕墻式防波堤周圍的流場變化圖，圖中箭頭方向為速度方向，箭線長度代表速度大小。從圖中可以看出，在下垂板周圍存在明顯的渦旋，渦旋使得入射波能量迅速耗散，反射減小。如圖21-a所示，t=0時刻前一個波浪傳播到直墻，由于直墻的阻礙作用，水體以一定的速度回流，此時下垂板前水體流速較小，對回流水體阻礙作用較弱，回流水體通過下垂板底部流出消浪室。如圖21-b所示，t=0.25T時刻隨著波浪的傳播，下垂板前水體流速增大，與回流水體相互作用分別在下垂板前與消浪室內形成一個順時針方向的渦旋與一個逆時針方向的渦旋。如圖21-c所示，t=0.5T時刻消浪室內的渦旋耗散了大量的波浪能量，消浪室內水體流速降低，渦旋消失，阻礙作用減弱，水體通過下垂板底部流入消浪室，消浪室內水位開始雍高。如圖22-d所示，t=0.75T時刻消浪室內水位進一步雍高，通過下垂板底部流入消浪室的水體受到直墻的阻礙作用，以一定的速度回流與流入消浪室的水體相互作用在消浪室內形成一個逆時針方向的渦旋。如圖21-e所示，t=T時刻流場進入下一個周期的循環。

21-a t=021-b t=0.25 T21-c t=0.5 T

21-d t=0.75 T21-e t=T圖21 幕墻式防波堤周圍流場變化圖Fig.21 Flow field around the curtain breakwater

4 結論

本文基于OpenFOAM建立了二維數值波浪水槽，研究了規則波與幕墻式防波堤的相互作用問題，分析了波長、消浪室寬度、下垂板吃水深度和下垂板底角對幕墻式防波堤水動力特性的影響規律。研究發現：波長，消浪室寬度，下垂板吃水深度和下垂板底角對幕墻式防波堤的反射系數與下垂板受力影響顯著，但是對直墻受力影響較小；綜合考慮反射系數與波浪力，推薦消浪室相對寬度B/L為0.08～0.10，下垂板相對吃水深度d/L為0.09～0.12，下垂板底角位于右下端。