基于OpenFOAM的擋板透空式防波堤水動(dòng)力特性數(shù)值分析

王國(guó)玉，顧新紅，秦世杰，于 明，張恩銘

(1.大連理工大學(xué)海岸和近海工程國(guó)家重點(diǎn)試驗(yàn)室，遼寧 大連 116024;2.山東港通工程管理咨詢有限公司， 山東 煙臺(tái) 264000)

近年來，隨著人們對(duì)海洋生態(tài)環(huán)境的重視，防波堤除了要抵御外海波浪、形成一個(gè)供船舶?？孔鳂I(yè)的掩蔽水域外，還需具有實(shí)現(xiàn)港內(nèi)外海水交換、發(fā)展生態(tài)漁業(yè)和休閑觀光等功能。因此，樁基擋板式透空防波堤越來越受到研究者和工程界的關(guān)注。

Ursell[1]于1947年首先理論分析了無限水深情況下直立薄板防波堤的透射系數(shù)與反射系數(shù)。Weigel[2]于1960利用微幅波理論推導(dǎo)了有限水深情況下直立薄板防波堤透反射系數(shù)的解析解。國(guó)內(nèi)關(guān)于透空式防波堤的研究始于1986年，邱大洪等[3]在假定波浪遇結(jié)構(gòu)物發(fā)生繞射的前提下，推導(dǎo)了任意水深下透射系數(shù)與反射系數(shù)的解析解，計(jì)算結(jié)果與理論值和試驗(yàn)值吻合較好，為透空式防波堤的研究奠定了基礎(chǔ)。王國(guó)玉等[4]利用理論分析方法對(duì)T型透空式防波堤進(jìn)行了研究，發(fā)現(xiàn)水平板的潛深和豎直板的高度對(duì)其消波性能有較大影響。嚴(yán)以新等[5]利用斷面模型試驗(yàn)，對(duì)多層擋板樁基透空堤消浪特性進(jìn)行了試驗(yàn)研究，得出了影響波浪透射率的主要因素，并證明了該種結(jié)構(gòu)可以建于水深浪大潮差大的海區(qū)。鐘瑚穗等[6]對(duì)小間距直樁式、柵欄式和直立擋板式透空堤開展物理模型試驗(yàn)研究，對(duì)比了3種模型的透射系數(shù)，并提出了消浪效果較好的新型組合式透空堤。王國(guó)玉等[7-8]提出了一種多層水平板式防波堤結(jié)構(gòu)，并通過物理模型試驗(yàn)、理論分析等方法重點(diǎn)研究了影響其消浪性能的因素。Laju等[9]設(shè)計(jì)了一種帶有雙側(cè)擋板的樁基式透空堤，分別應(yīng)用特征函數(shù)展開法和物理模型試驗(yàn)得到了結(jié)構(gòu)的透射系數(shù)，并證明了該結(jié)構(gòu)物適用于深海。Shin等[10]通過將理論值與試驗(yàn)值對(duì)比，分析了板間距、吃水深度對(duì)雙垂板透反射系數(shù)的影響，結(jié)果表明，當(dāng)兩板間距一定時(shí)，隨著板的相對(duì)入水深度增大，結(jié)構(gòu)物的反射系數(shù)增大，透射系數(shù)減小，同時(shí)也研究了板間的液面振蕩問題。程永舟等[11]研究了上下層平板相對(duì)間距以及格柵板間隙比等參數(shù)對(duì)防波堤透射系數(shù)和反射系數(shù)的影響。Alsaydalani 等[12]研究了三塊豎直板防波堤與波浪作用后的反射、透射、能量損耗、波浪力等水動(dòng)力特性。

現(xiàn)有關(guān)于擋板式透空防波堤的文獻(xiàn)研究結(jié)果表明，學(xué)者們主要集中于對(duì)防波堤透射系數(shù)和反射系數(shù)的研究，而對(duì)于擋板周圍流場(chǎng)的研究工作相對(duì)缺乏。特別是在追求較佳的掩護(hù)效果時(shí)，擋板往往入水深度較大，此時(shí)擋板底部海床表面的流速是工程人員關(guān)注的主要問題之一。因此，開展透空式防波堤底部流速的分析工作，將對(duì)擋板透空式防波堤的優(yōu)化設(shè)計(jì)以及減緩防波堤底部海床的沖刷具有重要的參考價(jià)值。

1 數(shù)學(xué)模型

二維數(shù)值波浪水槽的建立是基于開源CFD軟件OpenFOAM。流體控制方程的離散采用有限體積法，該方法可靈活地運(yùn)用結(jié)構(gòu)化網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格，便于模擬具有復(fù)雜邊界形狀結(jié)構(gòu)物周圍的流體運(yùn)動(dòng)。

1.1 控制方程

不可壓縮黏性流體的運(yùn)動(dòng)應(yīng)遵循如下連續(xù)性方程和動(dòng)量方程：

(1)

(2)

式中：ui為流體質(zhì)點(diǎn)在3個(gè)方向上的速度分量；ρ為密度；t為時(shí)間；p為壓強(qiáng)；gi為重力加速度；μe為有效動(dòng)力黏性系數(shù)。

1.2 邊界條件

1.2.1自由水面條件

VOF方法是追蹤自由液面的一個(gè)有效途徑，其基本思想是，在整個(gè)流場(chǎng)中定義一種流體體積與網(wǎng)格體積比值α的函數(shù)，其中α=1處的網(wǎng)格充滿該流體，α=0處的網(wǎng)格則不含有該流體，在0<α<1的網(wǎng)格區(qū)域含有自由面。體積分?jǐn)?shù)α需滿足輸運(yùn)方程[13]：

(3)

為避免自由液面附近出現(xiàn)兩相流界面模糊不清的問題，OpenFOAM對(duì)上述輸運(yùn)方程進(jìn)行了修正，通過添加人工壓縮項(xiàng)來提高體積分?jǐn)?shù)α的計(jì)算精度，修正后的方程為

(4)

式中：uri為垂直于兩相流界面的速度分量。

1.2.2造波、消波邊界

造波邊界波浪的產(chǎn)生方法為，依據(jù)Stokes 5階波浪理論，在造波邊界直接給出隨時(shí)間變化的自由表面高程和流體質(zhì)點(diǎn)速度分布，生成目標(biāo)規(guī)則波浪。水槽左側(cè)造波端設(shè)置一定厚度的海綿層，其作用為穩(wěn)定造波邊界產(chǎn)生的波浪并消除波浪二次反射的影響；類似的，水槽右側(cè)尾端設(shè)置一定厚度的海綿層，其作用為吸收向右傳播的波浪從而消除水槽尾端壁面產(chǎn)生的反射波浪。數(shù)值海綿層[14]的基本思想是，在每一個(gè)計(jì)算時(shí)間步長(zhǎng)內(nèi)將自由水面高程η、流體質(zhì)點(diǎn)速度u從計(jì)算值修正至目標(biāo)值，修正的公式如下：

φ=αRφcomputed+(1-αR)φtarget

(5)

式中：φ為需要被修正的物理量；φcomputed為由控制方程求出未經(jīng)修正的計(jì)算值；φtarget為目標(biāo)值；αR為修正權(quán)重，其在海綿層區(qū)域內(nèi)的分布用如下公式計(jì)算：

(6)

式中：χR為網(wǎng)格單元在海綿層內(nèi)的相對(duì)水平位置。

αR與χR的關(guān)系如圖1所示?？梢钥闯?，在海綿層區(qū)域內(nèi)，αR平滑的由0過渡到1，可以保證計(jì)算值的修正平滑進(jìn)行。

圖1 消波區(qū)αR(χR)變化示意圖

1.3 控制方程的離散及求解

數(shù)值計(jì)算采用有限體積法進(jìn)行，將計(jì)算區(qū)域劃分成一系列控制體積，再通過守恒型控制方程對(duì)控制體積做積分來導(dǎo)出一組離散方程，各項(xiàng)數(shù)值積分通過插值來計(jì)算。OpenFOAM提供了多種插值格式供選擇。本文所采用的插值格式如表1所示。壓力、速度耦合方程采用壓力的隱式算子分割算法(pressure implicit with splitting of operators， PISO)進(jìn)行求解。

表1 插值格式

2 數(shù)值波浪水槽

2.1 數(shù)值水槽建立

為了驗(yàn)證上述數(shù)學(xué)模型，在大連理工大學(xué)海岸和近海工程實(shí)驗(yàn)室的海洋環(huán)境波浪水槽(50.0 m×3.0 m×1.0 m)中，筆者設(shè)計(jì)并完成了規(guī)則波與樁基擋板透空式防波堤作用的物理模型試驗(yàn)(圖2)。試驗(yàn)水深為0.25 m，入射波高為0.04 m和0.08 m，對(duì)應(yīng)的波浪周期分別為0.91 s、1.10 s、1.28 s、1.46 s和1.64 s。結(jié)構(gòu)物布置在距離造波板27.0 m處，前側(cè)擋板入水深度為0.2 m，流速儀布置在擋板前沿處，流速測(cè)點(diǎn)位置在擋板下方距離水槽底部0.02 m位置處。試驗(yàn)中布置了5個(gè)浪高儀G1～G5用來觀測(cè)結(jié)構(gòu)后方的波面歷時(shí)，第一個(gè)浪高儀布置在距離結(jié)構(gòu)物后方1.35 m處，相鄰兩個(gè)浪高儀間距為0.50 m。在擋板的迎浪側(cè)設(shè)置了5個(gè)波浪壓力測(cè)點(diǎn)，測(cè)點(diǎn)位置如圖3所示。

圖2 擋板透空式防波堤模型試驗(yàn)布置(單位：m)

圖3 物理模型試驗(yàn)壓力測(cè)點(diǎn)布置示意圖(單位：mm)

圖4給出了基于OpenFOAM建立的二維數(shù)值波浪水槽。水槽長(zhǎng)度取為15.0 m，結(jié)構(gòu)物放在距造波板6.0 m處，距左端約3～5倍波長(zhǎng)。波面、流速和波壓力測(cè)點(diǎn)位置與物理模型試驗(yàn)中保持一致。水槽左端造波邊界位置處設(shè)置了一個(gè)寬度為2.0 m的海綿層，用來吸收由結(jié)構(gòu)物產(chǎn)生的反射波浪，從而消除數(shù)值波浪水槽中的波浪二次反射問題。類似的，水槽右側(cè)吸收邊界位置處設(shè)置了一個(gè)寬度為3.0 m的海綿層，用來吸收向右傳播的波浪。

圖4 數(shù)值水槽計(jì)算域布置示意圖(單位：m)

兼顧計(jì)算精度與時(shí)間，經(jīng)過一系列網(wǎng)格收斂性驗(yàn)證，背景網(wǎng)格的剖分依據(jù)如下[15]：波長(zhǎng)范圍內(nèi)100個(gè)網(wǎng)格，波高范圍內(nèi)10個(gè)網(wǎng)格，并對(duì)波面附近及結(jié)構(gòu)物周圍進(jìn)行了網(wǎng)格局部加密，以提高波面信息的捕捉精度。時(shí)間步長(zhǎng)的控制因素主要為Courant數(shù)，表示一個(gè)時(shí)間步長(zhǎng)內(nèi)允許流體質(zhì)點(diǎn)穿過的網(wǎng)格數(shù)量，Courant數(shù)越大，步長(zhǎng)越大，計(jì)算收斂的速度越快。本文中最大Courant數(shù)設(shè)置為0.25，并設(shè)置允許時(shí)間步長(zhǎng)根據(jù)Courant數(shù)進(jìn)行調(diào)節(jié)。

圖5為水槽中未布置結(jié)構(gòu)物(空水槽)時(shí)，左側(cè)造波邊界處海綿層內(nèi)G0處、右側(cè)吸收邊界海綿層內(nèi)部G6處和放置模型的位置Gm處3點(diǎn)在入射波浪的波高H=0.08 m和周期T=1.46 s時(shí)的波面η歷時(shí)過程。G0處的波高幾乎不發(fā)生變化，說明造波邊界處的海綿層可以很好地消除波浪的二次反射作用。同樣的，G6處的波高幾乎為0，這表明右側(cè)吸收邊界處的海綿層能夠有效地吸收向右傳播的波浪。而結(jié)構(gòu)物處Gm點(diǎn)的波面過程線分布均勻，波高隨時(shí)間衰減的相對(duì)誤差在2%以內(nèi)。

圖5 空水槽各測(cè)點(diǎn)波面歷時(shí)曲線(H=0.08 m，T=1.46 s)

2.2 數(shù)值水槽驗(yàn)證

2.2.1波面驗(yàn)證

圖6 不同周期數(shù)值模擬和物理模型試驗(yàn)堤后1.35 m處(G1)波面歷時(shí)曲線

基于上述數(shù)值波浪水槽，圖6分別為入射波高H=0.04 m、0.08 m時(shí)，不同周期的入射波浪向右傳播經(jīng)過擋板結(jié)構(gòu)物后在G1測(cè)點(diǎn)位置處的數(shù)值計(jì)算和物理模型試驗(yàn)觀測(cè)的波面η的歷時(shí)過程。各個(gè)波浪周期條件下，G1測(cè)點(diǎn)波面歷時(shí)過程的數(shù)值計(jì)算結(jié)果和試驗(yàn)觀測(cè)結(jié)果吻合良好。數(shù)值計(jì)算和試驗(yàn)觀測(cè)結(jié)果均表明，隨著入射波浪周期的增大，擋板后方G1測(cè)點(diǎn)的透射波高逐步增大。這反映出在擋板透空的條件下，入射波浪周期變長(zhǎng)，擋板結(jié)構(gòu)后方的掩護(hù)效果降低。

2.2.2波壓力驗(yàn)證

圖7 不同壓力測(cè)點(diǎn)數(shù)值模擬和物理模型試驗(yàn)的壓力歷時(shí)曲線(T=1.10 s)

圖7給出了入射波高H=0.04 m、0.08 m且波浪周期T=1.10 s時(shí)，擋板迎浪側(cè)的1～5號(hào)測(cè)點(diǎn)位置處波浪壓力歷時(shí)過程的數(shù)值計(jì)算結(jié)果和試驗(yàn)觀測(cè)結(jié)果的比較。數(shù)據(jù)結(jié)果顯示各測(cè)點(diǎn)壓力的數(shù)值計(jì)算結(jié)果與試驗(yàn)觀測(cè)結(jié)果吻合較好。水面以上的測(cè)點(diǎn)壓力值均存在數(shù)值為0的時(shí)間段，這是由于水面以上的壓力測(cè)點(diǎn)會(huì)出現(xiàn)出水而數(shù)據(jù)回零的情況。

2.2.3底流速驗(yàn)證

圖8給出了入射波高H=0.04 m、0.08 m時(shí)，不同周期的入射波浪向右傳播時(shí)擋板底部流速v的數(shù)值計(jì)算與物理模型試驗(yàn)觀測(cè)值的比較。各波浪周期條件下，擋板底部流速的數(shù)值計(jì)算結(jié)果與試驗(yàn)觀測(cè)結(jié)果吻合良好。數(shù)值計(jì)算結(jié)果和試驗(yàn)觀測(cè)數(shù)據(jù)均表明，入射波高或周期的增加均會(huì)增大擋板底部的流速。

上述分析表明，本文建立的數(shù)值波浪水槽所模擬得到的波浪與擋板式透空結(jié)構(gòu)相互作用的波面、波壓力和底流速歷時(shí)過程與實(shí)測(cè)結(jié)果均吻合較好。該數(shù)學(xué)模型可應(yīng)用于進(jìn)一步分析波浪作用下?lián)醢迨酵缚战Y(jié)構(gòu)的水動(dòng)力學(xué)特性。

3 擋板透空結(jié)構(gòu)水動(dòng)力特性分析

擋板透空結(jié)構(gòu)對(duì)波浪的掩護(hù)效果主要取決于擋板的入水深度。然而，隨著擋板入水深度的增加，結(jié)構(gòu)對(duì)波浪的反射作用增強(qiáng)，透射波高降低，但由于在水深方向上，水體在水平方向上的往復(fù)運(yùn)動(dòng)受到擋板的阻礙作用，擋板兩側(cè)水體不可避免的只能通過擋板下方與海床之間的孔隙進(jìn)行交換，這勢(shì)必會(huì)在很大程度上增加擋板下方的底部流速。為進(jìn)一步分析擋板透空結(jié)構(gòu)的掩護(hù)效果和底部流場(chǎng)特性，基于上述數(shù)值波浪水槽，通過改變結(jié)構(gòu)物的相對(duì)入水深度(擋板入水深度d/水深h)，對(duì)規(guī)則波與透空結(jié)構(gòu)相互作用問題展開了數(shù)值計(jì)算，重點(diǎn)分析擋板的相對(duì)入水深度對(duì)結(jié)構(gòu)掩護(hù)效果和結(jié)構(gòu)下方最大底流速的影響。數(shù)值計(jì)算的條件和波浪要素與上述物理模型試驗(yàn)中一致，并選取相對(duì)入水深度分別為d/h=0.2、0.4和0.8進(jìn)行水動(dòng)力特性分析。

圖8 不同周期數(shù)值模擬和物理模型試驗(yàn)的流速歷時(shí)曲線

3.1 透射系數(shù)

圖9 不同相對(duì)入水深度時(shí)G1處波面歷時(shí)曲線比較(T=1.28 s)

防波堤的透反射系數(shù)是衡量防波堤性能的重要指標(biāo)之一，而影響其透反射系數(shù)的因素很多，如相對(duì)水深、相對(duì)板寬、相對(duì)波長(zhǎng)和擋板相對(duì)入水深度等，這里主要討論擋板相對(duì)入水深度對(duì)透射系數(shù)的影響。

圖9為在入射波高H=0.04 m、0.08 m和波浪周期T=1.28 s時(shí)，不同相對(duì)入水深度情況下G1處的波面歷時(shí)曲線比較。計(jì)算結(jié)果表明，隨著擋板相對(duì)入水深度的增加，擋板后方的透射波高逐漸降低，說明增加入水深度可以有效地增強(qiáng)結(jié)構(gòu)的掩護(hù)效果。

圖10為H=0.04 m、0.08 m時(shí)，不同相對(duì)入水深度條件下透射系數(shù)Ct隨波數(shù)kh的變化規(guī)律，其中k=2π/L?？梢钥闯?，當(dāng)波數(shù)一定時(shí)，隨著擋板相對(duì)入水深度的增加，透射系數(shù)逐漸減小，說明擋板入水深度增加，結(jié)構(gòu)的掩護(hù)效果增強(qiáng)。另一方面，隨著波數(shù)增加，入射波浪的周期減小，透射系數(shù)隨之明顯降低，這說明對(duì)于波長(zhǎng)較小或者周期較小的短波情形，擋板式透空防波堤的效果更佳。數(shù)值模擬結(jié)果顯示，當(dāng)擋板相對(duì)入水深度增大到0.8時(shí)，對(duì)于絕大部分入射波浪情形，透射系數(shù)均可降低至0.5以下。

圖10 不同相對(duì)入水深度的擋板式透空防波堤的透射系數(shù)隨波數(shù)的變化

3.2 底流速

對(duì)于工程實(shí)際而言，當(dāng)?shù)琢魉佥^大時(shí)海床會(huì)發(fā)生沖刷侵蝕，這將影響海洋生態(tài)環(huán)境同時(shí)也降低結(jié)構(gòu)物的穩(wěn)定性，因此本節(jié)主要討論擋板透空式防波堤的底部流速問題，包括正向流速和反向流速。

圖11為在入射波高H=0.04 m、0.08 m和波浪周期T=1.28 s時(shí)，不同相對(duì)入水深度情況下，擋板的底流速歷時(shí)曲線。計(jì)算結(jié)果表明，不同擋板相對(duì)入水深度條件下，擋板底流速變化趨勢(shì)一致，隨著擋板相對(duì)入水深度的增加，底流速明顯增大。

圖11 不同相對(duì)入水深度時(shí)擋板底流速歷時(shí)曲線比較(T=1.28 s)

圖12為入射波高H=0.04 m、0.08 m時(shí)，擋板不同相對(duì)入水深度情形下，最大正向和負(fù)向無量綱底流速V隨著波數(shù)的變化關(guān)系。底流速的無量綱化是選取正向或負(fù)向流速的最大值與波浪傳播速度的比值。其中，正向表示流速方向與波浪傳播方向一致，負(fù)向則表示相反方向。

圖12 不同入水深度情況下底部無量綱流速隨波數(shù)的變化關(guān)系

圖12中的最大無量綱底流速數(shù)據(jù)結(jié)果顯示最大流速的分布在正向和負(fù)向上具有一定的不對(duì)稱性。最大流速的幅值隨著波數(shù)的增加而減小，這表明對(duì)于高頻或短周期波浪而言，其傳播經(jīng)過擋板透空結(jié)構(gòu)時(shí)，在擋板下方形成的往復(fù)流的底流速相對(duì)較小。另外，隨著擋板相對(duì)入水深度的增加，最大底流速的幅值明顯增大。特別是在擋板相對(duì)入水深度較大時(shí)，如d/h=0.8，最大底流速的幅值成倍增加。而且，隨著入射波高的增加，底流速也大幅增加。這一點(diǎn)在工程設(shè)計(jì)上需要特別注意。一般而言，為了謀求更佳的掩護(hù)效果，工程人員通常需要增大擋板的入水深度，以達(dá)到有效降低透射波高的目的。但同時(shí)會(huì)使得底部流速大幅增加，這將會(huì)對(duì)海床形成局部沖刷，進(jìn)而對(duì)結(jié)構(gòu)的安全形成威脅。實(shí)際應(yīng)用時(shí)，工程設(shè)計(jì)人員需結(jié)合底流速的具體情況，增設(shè)人工護(hù)底，以保證海床的穩(wěn)定性。

4 結(jié) 語

本文應(yīng)用開源計(jì)算流體力學(xué)軟件OpenFOAM，基于有限體積法建立了波浪與擋板透空結(jié)構(gòu)相互作用的二維數(shù)值波浪水槽，并通過規(guī)則波與擋板透空式防波堤相互作用的物理模型試驗(yàn)進(jìn)行了驗(yàn)證，結(jié)果顯示波面、壓力和流速的數(shù)據(jù)結(jié)果吻合良好。數(shù)值計(jì)算結(jié)果表明，增加擋板的相對(duì)入水深度可有效降低波浪的透射系數(shù)。當(dāng)相對(duì)入水深度為0.8時(shí)，擋板透空結(jié)構(gòu)對(duì)波浪的透射系數(shù)可降低至0.5以下。另外，隨著擋板相對(duì)入水深度的增加，底部水體往復(fù)流動(dòng)的速度也明顯增大，這需要引起工程技術(shù)人員的額外關(guān)注。