尹勝強(qiáng),劉 勇,王心玉
(中國海洋大學(xué) 海岸與海洋工程研究所,青島 266100)
圖1 幕墻式防波堤結(jié)構(gòu)示意圖Fig.1 Sketch of a curtain breakwater
幕墻式防波堤是一種新型防波堤,由直立板防波堤逐漸發(fā)展而來。如圖1所示,幕墻式防波堤由前下垂板和后直墻兩部分組成,下垂板與直墻之間形成消浪室。消浪室促使下垂板周圍的渦旋明顯增大,入射波能量耗散增加。與傳統(tǒng)防波堤相比,幕墻式防波堤可以有效降低波浪反射,保證結(jié)構(gòu)前水域平穩(wěn)。
許多學(xué)者對(duì)直立板結(jié)構(gòu)的水動(dòng)力特性進(jìn)行了研究。Ursell[1]最早研究了無限水深中波浪與直立板結(jié)構(gòu)的相互作用問題,給出了反射系數(shù)和透射系數(shù)的變化規(guī)律。Evans[2]研究了波浪與淹沒式直立板結(jié)構(gòu)的相互作用,得到了一階波浪力與二階波浪力的計(jì)算公式。Morris[3]建立了一種新的變分方法,并應(yīng)用于分析波浪與非對(duì)稱直立板結(jié)構(gòu)的相互作用問題。Losada等[4]和Abul-Azm[5]采用匹配特征函數(shù)展開法分別研究了正向波和斜向波與直立板結(jié)構(gòu)的相互作用問題,給出了反射系數(shù)和透射系數(shù)的變化規(guī)律。Kriebel和Bollmann[6]利用能量透射理論和匹配特征函數(shù)展開法研究了波浪與出水實(shí)體薄板結(jié)構(gòu)的相互作用,發(fā)現(xiàn)兩種方法計(jì)算得到的反射系數(shù)基本一致。陳雪峰等[7]研究了波浪與開孔板結(jié)構(gòu)的相互作用,給出了開孔率等因素對(duì)開孔板點(diǎn)壓力分布和反射系數(shù)的影響規(guī)律。
在實(shí)際工程中,為了提高掩護(hù)效果,經(jīng)常采用雙排或多排直立板結(jié)構(gòu)。與單排直立板結(jié)構(gòu)相比,雙排或多排直立板結(jié)構(gòu)可以更有效地反射入射波能量。Porter[8]采用多項(xiàng)伽遼金方法研究了波浪與雙排實(shí)體板的相互作用,給出了反射系數(shù)和透射系數(shù)的精確計(jì)算結(jié)果。Liu等[9]結(jié)合理論分析和物理模型試驗(yàn),研究了波浪對(duì)Jarlan型透空潛堤(前開孔板和后實(shí)體板組成)的作用,理論計(jì)算結(jié)果與試驗(yàn)結(jié)果符合良好,當(dāng)相對(duì)板間距為 0.3~0.4時(shí),潛堤掩護(hù)效果最佳。Clauss等[10]結(jié)合物理模型試驗(yàn)和理論分析,研究了多排垂直開孔板結(jié)構(gòu)的水動(dòng)力特性。
Nakamura等[11-12]在研究雙排直立板結(jié)構(gòu)水動(dòng)力特性的基礎(chǔ)上,采用直墻代替后側(cè)垂直板,提出了幕墻式防波堤的概念,并通過物理模型試驗(yàn)研究了消浪室寬度、下垂板吃水深度對(duì)幕墻式防波堤反射系數(shù)的影響。此后,Nakamura團(tuán)隊(duì)針對(duì)幕墻式防波堤及其改進(jìn)結(jié)構(gòu)開展了大量物理模型試驗(yàn),深入研究了幕墻式防波堤的水動(dòng)力特性。Nishikawa等[13]對(duì)幕墻式防波堤的受力特性進(jìn)行了理論和試驗(yàn)研究,分析了入射波高和消浪室相對(duì)寬度對(duì)下垂板與直墻受力的影響。Ono等[14]采用數(shù)值模擬結(jié)合物理模型試驗(yàn)的方法,分析了波浪周期對(duì)幕墻式防波堤周圍流場(chǎng)特性的影響。耿寶磊等[15]研究了波浪與直墻前多層透空薄板結(jié)構(gòu)的相互作用,得到了反射系數(shù)、透射系數(shù)與波能吸收率的變化規(guī)律。
綜上所述,現(xiàn)有關(guān)于幕墻式防波堤水動(dòng)力特性的研究以理論分析和物理模型試驗(yàn)為主,關(guān)于幕墻式防波堤水動(dòng)力特性的數(shù)值模擬研究較少。與理論分析和物理模型試驗(yàn)相比,數(shù)值模擬可以有效分析幕墻式防波堤周圍的流場(chǎng)特性,從而深入理解幕墻式防波堤的消浪機(jī)理。為此,本文基于CFD開源程序庫OpenFOAM建立二維數(shù)值波浪水槽,模擬波浪對(duì)幕墻式防波堤的作用,研究波長(zhǎng)、消浪室寬度、下垂板吃水深度和下垂板底角對(duì)幕墻式防波堤反射系數(shù)、波浪力和流場(chǎng)特性的影響,分析幕墻式防波堤的消浪機(jī)理,研究結(jié)果可為工程設(shè)計(jì)提供科學(xué)指導(dǎo)。
將流體假設(shè)為不可壓縮粘性流體,在計(jì)算過程中滿足連續(xù)性方程和動(dòng)量方程。
連續(xù)性方程
(1)
動(dòng)量方程
(2)
式中:ρ為流體密度;xi為空間點(diǎn)的坐標(biāo);ui為在t時(shí)刻坐標(biāo)xi點(diǎn)的速度分量;p為壓力;μ為流體動(dòng)力粘性系數(shù);fi為源項(xiàng);針對(duì)本文的二維問題而言,i、j的取值范圍為1、2,表示直角坐標(biāo)系下,各個(gè)矢量在x、y方向上的分量。
采用VOF方法(流體體積法)捕捉自由面。定義體積分?jǐn)?shù)函數(shù)α表示計(jì)算網(wǎng)格單元中液相所占體積分?jǐn)?shù)。
體積分?jǐn)?shù)函數(shù)α滿足連續(xù)性方程
(3)
式中:ur是為了實(shí)現(xiàn)對(duì)自由液面的壓縮而引入的只在自由面起作用的相對(duì)速度。由下式確定
(4)
式中:φ為通過自由面單元表面的體積流量;Sf為對(duì)應(yīng)的單元表面面積;nf為自由面單位法向向量;Cα為自由面壓縮系數(shù)。
通過公式(3)計(jì)算得到體積分?jǐn)?shù)α后,氣液兩相流場(chǎng)內(nèi)各點(diǎn)的密度和粘性系數(shù)可以通過體積分?jǐn)?shù)函數(shù)α表示為
ρ=αwaterρwater+(1+αwater)ρa(bǔ)ir
(5)
μ=αwaterμwater+(1+αwater)μair
(6)
采用速度入口邊界造波方法造波,通過在入口邊界上設(shè)置波面和流體速度實(shí)現(xiàn)造波。本文二維數(shù)值波浪水槽采用二階斯托克斯波,二階斯托克斯波的波面和流體速度如下式所示
(7)
(8)
(9)
式中:η為波面變化;H為波高;k為波數(shù);ω為波浪圓頻率;u為水平質(zhì)點(diǎn)速度;v為垂直質(zhì)點(diǎn)速度;h為水深。
采用松弛區(qū)域法消波,在二維數(shù)值波浪水槽的入口區(qū)域和出口區(qū)域設(shè)置松弛區(qū)域,在松弛區(qū)域內(nèi)每一時(shí)刻對(duì)速度進(jìn)行如下修正
u=αRumodel+(1-αR)utarget
(10)
式中:umodel為通過求解N-S方程得到的速度值;utarget為期望得到的目標(biāo)速度值;α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
采用有限體積法對(duì)控制方程進(jìn)行離散求解,計(jì)算過程中將計(jì)算域劃分為若干網(wǎng)格單元,在每個(gè)網(wǎng)格單元內(nèi)對(duì)控制方程進(jìn)行積分,再通過高斯變化得到離散后的控制方程。
采用GAMG(Geometric-Algebraic Multi-Grid)方法求解離散后的線性代數(shù)方程組,先將計(jì)算域網(wǎng)格合并為少數(shù)粗網(wǎng)格,在粗網(wǎng)格上快速求解得到方程組的解,然后將粗網(wǎng)格上得到的方程組的解作為細(xì)網(wǎng)格上的預(yù)測(cè)值,在細(xì)網(wǎng)格上對(duì)線性代數(shù)方程組迭代求解,最終得到精確的數(shù)值解。
求解過程中采用同位網(wǎng)格存儲(chǔ),速度、壓力等均存儲(chǔ)在計(jì)算網(wǎng)格中心處,計(jì)算得到網(wǎng)格中心處的流場(chǎng)變量值后,再通過插值方法獲得其他位置的流場(chǎng)信息。
圖3 二維數(shù)值波浪水槽Fig.3 Sketch of two-dimensional numerical tank
如圖3所示,二維數(shù)值波浪水槽長(zhǎng)12 m,高0.8 m。幕墻式防波堤由前下垂板和后直墻兩部分組成,其中下垂板厚0.042 m,直墻厚0.15 m。試驗(yàn)布置參考Nakamura[11]的物理模型試驗(yàn)(比尺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è)置。所有的數(shù)值試驗(yàn)條件均列于表1。
表1 數(shù)值試驗(yàn)條件Tab.1 Numerical test conditions
本文建立的二維數(shù)值水槽采用結(jié)構(gòu)化網(wǎng)格,以矩形網(wǎng)格為主,存在少量的不規(guī)則四邊形網(wǎng)格。經(jīng)網(wǎng)格收斂性檢驗(yàn),最終選用圖4所示的網(wǎng)格劃分方式,計(jì)算網(wǎng)格縱橫比為1/4,沿水平方向計(jì)算網(wǎng)格尺寸不超過0.02 m, 沿豎直方向計(jì)算網(wǎng)格尺寸不超過0.005 m。為了更精確地模擬自由面與結(jié)構(gòu)物周圍的流場(chǎng)變化,對(duì)自由面與結(jié)構(gòu)物附近的網(wǎng)格做局部加密。網(wǎng)格數(shù)量保持在15萬左右。
圖4 計(jì)算網(wǎng)格劃分Fig.4 Sketch of grid generation
對(duì)于數(shù)值模擬而言,波浪在傳播過程中,隨著傳播距離的增大,由于流體粘性、數(shù)值截?cái)嗾`差以及邊界作用的影響波高會(huì)逐漸衰減。因此,在數(shù)值模擬試驗(yàn)開始之前,需要在二維數(shù)值波浪水槽中進(jìn)行率波。
選取典型波浪參數(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處模擬值與理論值的對(duì)比結(jié)果進(jìn)行說明。如圖5~圖7所示,點(diǎn)為基于二階斯托克斯波的數(shù)值模擬結(jié)果,線為二階斯托克斯波的理論結(jié)果??梢钥闯?,數(shù)值模擬結(jié)果波形規(guī)則、穩(wěn)定,說明二維數(shù)值波浪水槽可以有效消除造波邊界和水槽末端波浪反射的影響。數(shù)值模擬結(jié)果與理論結(jié)果存在很小的差異,這是由于在數(shù)值模擬過程中考慮了流體的粘性,存在耗散現(xiàn)象,而理論結(jié)果沒有考慮流體粘性??傮w而言,數(shù)值模擬結(jié)果與理論結(jié)果符合良好。
圖5 波面歷時(shí)曲線(T=0.9 s,H=6 cm)Fig.5 Duration curve of wave surface(T=0.9 s and H= 6 cm) 圖6 波面歷時(shí)曲線(T= 1.4 s,H= 6 cm)Fig.6 Duration curve of wave surface(T= 1.4 s and H= 6 cm)圖7 波面歷時(shí)曲線(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)相對(duì)波高Hc/H的數(shù)值模擬結(jié)果與Nakamura等[10]的物理模型試驗(yàn)結(jié)果進(jìn)行對(duì)比驗(yàn)證。在計(jì)算反射系數(shù)與消浪室內(nèi)相對(duì)波高時(shí),以至少十個(gè)穩(wěn)定波為準(zhǔn)。圖8~圖10給出了三種不同工況下反射系數(shù)的數(shù)值模擬結(jié)果與Nakamura等[11]的物理模型試驗(yàn)結(jié)果對(duì)比,圖中實(shí)心點(diǎn)表示數(shù)值模擬結(jié)果,空心點(diǎn)表示物理模型試驗(yàn)結(jié)果??梢钥闯?,數(shù)值模擬結(jié)果與物理模型試驗(yàn)結(jié)果符合良好,此外,隨著波長(zhǎng)的變化,反射系數(shù)存在最小值(約為0.1),說明通過合理選取結(jié)構(gòu)參數(shù),幕墻式防波堤可以非常有效地耗散入射波能量。
圖8 反射系數(shù)的數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果[11]對(duì)比(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é)果與試驗(yàn)結(jié)果[11]對(duì)比(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é)果與試驗(yàn)結(jié)果[11]對(duì)比(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)相對(duì)波高的數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果[11]對(duì)比(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)相對(duì)波高的數(shù)值模擬結(jié)果與Nakamura等[11]的物理模型試驗(yàn)結(jié)果對(duì)比,圖中實(shí)心點(diǎn)表示數(shù)值模擬結(jié)果,空心點(diǎn)表示物理模型試驗(yàn)結(jié)果??梢钥闯觯瑪?shù)值模擬結(jié)果與物理模型試驗(yàn)結(jié)果符合良好。
3.1.1 消浪室寬度與下垂板吃水深度的影響
圖12給出了不同消浪室寬度下,反射系數(shù)Cr的數(shù)值計(jì)算結(jié)果隨下垂板相對(duì)吃水深度d/L的變化曲線。從圖12可以看出:對(duì)于長(zhǎng)周期波而言,隨著消浪室寬度增大,波浪能量耗散增加,反射系數(shù)降低;但是對(duì)于短周期波而言,變化規(guī)律正好相反。圖13給出了不同下垂板吃水深度下,反射系數(shù)Cr的數(shù)值計(jì)算結(jié)果隨消浪室相對(duì)寬度B/L的變化曲線。從圖13可以看出:對(duì)于長(zhǎng)周期波而言,隨著下垂板吃水深度增大,防波堤的消浪效果增強(qiáng);對(duì)于短周期波而言,隨著下垂板吃水深度增大,防波堤的消浪效果減弱。此外,當(dāng)消浪室相對(duì)寬度為0.08~0.10,下垂板相對(duì)吃水深度為0.09~0.12時(shí),幕墻式防波堤消浪效果最好。從圖13還可以看出,隨著下垂板吃水深度的增大,反射系數(shù)最小值對(duì)應(yīng)的消浪室相對(duì)寬度逐漸減小。
圖12 反射系數(shù)隨下垂板相對(duì)吃水深度的變化(d=21 cm)Fig.12 Variation of the reflection coefficient versus the relative immersed depth of the drooping plate(d=21 cm) 圖13 反射系數(shù)隨消浪室相對(duì)寬度的變化(B=29 cm)Fig.13 Variation of the reflection coefficient versus the relative wave chamber width(B=29 cm)
3.1.2 下垂板底角的影響
圖14 反射系數(shù)隨消浪室相對(duì)寬度的變化(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ù)值計(jì)算結(jié)果隨消浪室相對(duì)寬度B/L的變化曲線。從圖14可以看出:下垂板底角越尖銳,耗散的波浪能量越多,反射系數(shù)越?。慌c長(zhǎng)周期波相比,短周期波條件下反射系數(shù)對(duì)于下垂板底角的變化更加敏感。此外,下垂板底角的方向(θ=45°或θ=-45°)對(duì)于反射系數(shù)的變化幾乎沒有影響。
3.2.1 消浪室寬度的影響
圖15給出了不同消浪室寬度下,下垂板波浪力的數(shù)值計(jì)算結(jié)果隨下垂板相對(duì)吃水深度d/L的變化曲線。從圖15可以看出:對(duì)于長(zhǎng)周期波而言,隨著消浪室寬度的增加,下垂板受力增加;對(duì)于短周期波而言,隨著消浪室寬度增加,下垂板受力減小。
圖16給出了不同消浪室寬度下,直墻波浪力的數(shù)值計(jì)算結(jié)果隨下垂板相對(duì)吃水深度d/L的變化曲線。從圖16可以看出:隨著消浪室寬度增加,直墻受力減少;隨著下垂板相對(duì)吃水深度增大(波長(zhǎng)減小),直墻受力減小。
圖15 下垂板受力隨下垂板相對(duì)吃水深度的變化(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 直墻受力隨下垂板相對(duì)吃水深度的變化 (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ù)值計(jì)算結(jié)果隨消浪室相對(duì)寬度B/L的變化曲線。從圖17可以看出:隨著下垂板吃水深度增加,下垂板受力明顯增大,這是因?yàn)殡S著下垂板吃水深度增加,下垂板的受力面積增大;長(zhǎng)周期波條件下,下垂板受力對(duì)下垂板吃水深度的變化更加敏感。另外,當(dāng)消浪室相對(duì)寬度B/L=0.15時(shí),下垂板受力最大,約為0.6。
圖18給出了不同下垂板吃水深度下,直墻波浪力的數(shù)值計(jì)算結(jié)果隨消浪室相對(duì)寬度B/L的變化曲線。從圖18可以看出:隨著下垂板吃水深度增加,直墻受力減??;隨著消浪室相對(duì)寬度增大(波長(zhǎng)減小),直墻受力減小。
圖17 下垂板受力隨消浪室相對(duì)寬度的變化(B=29 cm)Fig.17 Variation of wave force on drooping plate versus the relative wave chamber width (B=29 cm)圖18 直墻受力隨消浪室相對(duì)寬度的變化(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ù)值計(jì)算結(jié)果隨消浪室相對(duì)寬度B/L的變化曲線。從圖19可以看出:隨著下垂板底角的增大,下垂板受力增大;θ=45°(右下端)時(shí)下垂板受力明顯小于θ=-45°(左下端)。
圖20給出了不同下垂板底角下,直墻波浪力的數(shù)值計(jì)算結(jié)果隨消浪室相對(duì)寬度B/L的變化曲線。從圖20可以看出:隨著下垂板底角增大,直墻受力的變化減小,最大變化值約為0.1(變化率僅為3%)。因此在實(shí)際工程應(yīng)用中,下垂板底角對(duì)直墻受力的影響基本可以忽略。
圖19 下垂板受力隨下垂板相對(duì)吃水深度的變化(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 直墻受力隨消浪室相對(duì)寬度的變化(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)
本節(jié)以T=1.4 s,B=29 cm,d=21 cm,θ= 45°的工況為例,分析下垂板周圍的流場(chǎng)特性。圖21-a~21-e給出了一個(gè)波浪周期內(nèi)典型時(shí)刻幕墻式防波堤周圍的流場(chǎng)變化圖,圖中箭頭方向?yàn)樗俣确较?,箭線長(zhǎng)度代表速度大小。從圖中可以看出,在下垂板周圍存在明顯的渦旋,渦旋使得入射波能量迅速耗散,反射減小。如圖21-a所示,t=0時(shí)刻前一個(gè)波浪傳播到直墻,由于直墻的阻礙作用,水體以一定的速度回流,此時(shí)下垂板前水體流速較小,對(duì)回流水體阻礙作用較弱,回流水體通過下垂板底部流出消浪室。如圖21-b所示,t=0.25T時(shí)刻隨著波浪的傳播,下垂板前水體流速增大,與回流水體相互作用分別在下垂板前與消浪室內(nèi)形成一個(gè)順時(shí)針方向的渦旋與一個(gè)逆時(shí)針方向的渦旋。如圖21-c所示,t=0.5T時(shí)刻消浪室內(nèi)的渦旋耗散了大量的波浪能量,消浪室內(nèi)水體流速降低,渦旋消失,阻礙作用減弱,水體通過下垂板底部流入消浪室,消浪室內(nèi)水位開始雍高。如圖22-d所示,t=0.75T時(shí)刻消浪室內(nèi)水位進(jìn)一步雍高,通過下垂板底部流入消浪室的水體受到直墻的阻礙作用,以一定的速度回流與流入消浪室的水體相互作用在消浪室內(nèi)形成一個(gè)逆時(shí)針方向的渦旋。如圖21-e所示,t=T時(shí)刻流場(chǎng)進(jìn)入下一個(gè)周期的循環(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 幕墻式防波堤周圍流場(chǎng)變化圖Fig.21 Flow field around the curtain breakwater
本文基于OpenFOAM建立了二維數(shù)值波浪水槽,研究了規(guī)則波與幕墻式防波堤的相互作用問題,分析了波長(zhǎng)、消浪室寬度、下垂板吃水深度和下垂板底角對(duì)幕墻式防波堤水動(dòng)力特性的影響規(guī)律。研究發(fā)現(xiàn):波長(zhǎng),消浪室寬度,下垂板吃水深度和下垂板底角對(duì)幕墻式防波堤的反射系數(shù)與下垂板受力影響顯著,但是對(duì)直墻受力影響較??;綜合考慮反射系數(shù)與波浪力,推薦消浪室相對(duì)寬度B/L為0.08~0.10,下垂板相對(duì)吃水深度d/L為0.09~0.12,下垂板底角位于右下端。