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

尹勝強，劉 勇，王心玉

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

圖1 幕墻式防波堤結(jié)構(gòu)示意圖Fig.1 Sketch of a curtain breakwater

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

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

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

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

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

1 基于OpenFOAM的數(shù)值模型

1.1 基本控制方程

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

連續(xù)性方程

(1)

動量方程

(2)

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

1.2 自由面處理方法

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

體積分數(shù)函數(shù)α滿足連續(xù)性方程

(3)

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

(4)

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

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

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

(5)

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

(6)

1.3 數(shù)值造波方法

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

(7)

(8)

(9)

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

1.4 數(shù)值消波方法

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

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

(10)

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

(11)

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

圖2 松弛區(qū)域內(nèi)αR和XR的變化圖Fig.2 Variation of αR versus XR in relaxation zone

1.5 數(shù)值求解方法

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

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

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

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

2 數(shù)值模型驗證

2.1 數(shù)值水槽設(shè)置

如圖3所示，二維數(shù)值波浪水槽長12 m，高0.8 m。幕墻式防波堤由前下垂板和后直墻兩部分組成，其中下垂板厚0.042 m，直墻厚0.15 m。試驗布置參考Nakamura[11]的物理模型試驗(比尺1:12)，水深42 cm，波高6 cm，波浪周期0.9～1.8 s；幕墻式防波堤距離二維數(shù)值水槽末端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。出口區(qū)域的消波區(qū)為率定入射波時使用，放入防波堤后不再設(shè)置。所有的數(shù)值試驗條件均列于表1。

表1 數(shù)值試驗條件Tab.1 Numerical test conditions

2.2 計算網(wǎng)格

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

圖4 計算網(wǎng)格劃分Fig.4 Sketch of grid generation

2.3 數(shù)值結(jié)果驗證

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

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

圖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)

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

圖8 反射系數(shù)的數(shù)值模擬結(jié)果與試驗結(jié)果[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 反射系數(shù)的數(shù)值模擬結(jié)果與試驗結(jié)果[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 反射系數(shù)的數(shù)值模擬結(jié)果與試驗結(jié)果[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 消浪室內(nèi)相對波高的數(shù)值模擬結(jié)果與試驗結(jié)果[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給出了消浪室內(nèi)相對波高的數(shù)值模擬結(jié)果與Nakamura等[11]的物理模型試驗結(jié)果對比，圖中實心點表示數(shù)值模擬結(jié)果，空心點表示物理模型試驗結(jié)果。可以看出，數(shù)值模擬結(jié)果與物理模型試驗結(jié)果符合良好。

3 數(shù)值算例與討論

3.1 反射系數(shù)

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

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

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

3.1.2 下垂板底角的影響

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

3.2 波浪力

3.2.1 消浪室寬度的影響

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

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

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

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

圖20給出了不同下垂板底角下，直墻波浪力的數(shù)值計算結(jié)果隨消浪室相對寬度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 流場特性分析

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

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 結(jié)論

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