高俊亮， 張一兆， 何志偉， 王崗

(1.江蘇科技大學(xué) 船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江 212003； 2.河海大學(xué) 港口海岸與近海工程學(xué)院，江蘇 南京 210098)

平板式結(jié)構(gòu)物廣泛應(yīng)用于沿海和海洋工程中，但是由于海嘯波的沖擊和設(shè)計不當(dāng)?shù)仍颍?jīng)常會導(dǎo)致結(jié)構(gòu)物受損造成巨大的災(zāi)難。眾所周知，孤立波是研究海嘯波沖擊海岸結(jié)構(gòu)物時模擬海嘯型波的理想選擇，因此，孤立波作用下固定式水平板的水動力特性研究對于海洋工程結(jié)構(gòu)物設(shè)計和人們的生產(chǎn)活動都有重要意義和參考價值。

關(guān)于板式結(jié)構(gòu)物的研究已獲得諸多成果，許多學(xué)者采用典型結(jié)構(gòu)水平板為研究對象開展了一系列關(guān)于波浪作用于水下板[1-4]和水上板[5-7]的實(shí)驗(yàn)研究。近年來，隨著計算流體力學(xué)(computational fluid dynamics,CFD)方法的發(fā)展，一些學(xué)者開始采用CFD的方法研究孤立波與水平板之間的相互作用問題[8-10]。You等[11]采用CIP的方法研究了孤立波作用于水下板的水動力特性。Hayatdavoodi等[12]利用Green-Naghdi方程研究了孤立波和橢圓余弦波作用下，波高、淹沒深度和平板寬度等因素對水平板波散、透射系數(shù)和反射系數(shù)的影響。Seiffert等[13]基于OpenFOAM?中的interFoam求解器研究了波幅、淹沒深度和水上高度等因素對平板波浪力的影響，但其研究只涉及到了3個凈空高度且高度較小，和Xu等[14]一樣未能全面地分析凈空高度對波浪力的影響。綜合來看，關(guān)于孤立波作用下水平板的水動力特性研究還不夠全面，水上板的凈空高度和板寬等因素對水平板的水動力特性的影響分析還不夠全面，下文將根據(jù)水上板的凈空高度和板寬2個因素對水平板水動力特性進(jìn)行研究。

本文基于OpenFOAM?建立二維數(shù)值模型，以Navier-Stokes方程為控制方程，采用流體體積方法進(jìn)行對自由面進(jìn)行捕捉描述，運(yùn)用waves2Foam進(jìn)行造波和消波。應(yīng)用該數(shù)值模型，先是對淹沒固定平板上的波浪力、表面高程進(jìn)行了數(shù)值驗(yàn)證，證明其計算波浪力結(jié)果和模擬波浪運(yùn)動狀態(tài)的準(zhǔn)確性。進(jìn)而對不同的板寬和相對凈空比參數(shù)下平板上的波浪力和波面變化情況進(jìn)行計算和分析，討論了板寬和相對凈空比對波浪力、透射系數(shù)和反射系數(shù)的影響規(guī)律。

1 數(shù)值模型描述

本文求解器采用粘性流體求解器。基于OpenFOAM?多相求解器interFoam,利用Jacobsen等[15]提出的基于松弛方法的波浪生成工具箱“waves2Foam”進(jìn)行造波和消波。

1.1 控制方程

數(shù)值模型采用連續(xù)性和Navier-Stokes方程求解水和空氣兩相流動的控制方程：

(1)

(2)

數(shù)值模型采用流體體積(VOF)方法對自由面進(jìn)行捕捉。α為數(shù)值模型中每個計算單元的體積分?jǐn)?shù)，其定義為：

(3)

α的分布滿足對流輸運(yùn)方程：

(4)

式中ur=uwater-uair表示為水和空氣之間的相對速度。通過使用α，任何流體性質(zhì)的空間變化φ(例如流體密度ρ和動態(tài)粘度μ)都可以通過加權(quán)平均法來表示:

φ=αφwater+(1-α)φair

(5)

1.2 邊界條件與數(shù)值實(shí)現(xiàn)

使用Jacobsen等[15]提出的“waves2Foam”工具箱在入口邊界處進(jìn)行造波，速度設(shè)置為入射波波速，波速為零。在入口和出口邊界設(shè)置了2個松弛區(qū)，分別用來吸收反射波和透射波。在水平板的上部，邊界條件設(shè)置為“大氣”；水平板的底部和水槽箱的邊界條件設(shè)置為“無滑移”；本模型為二維模型、三維墻上的邊界條件設(shè)置為“空”。

基于有限體積法求解控制方程 (1) 和 (2) 以及對流輸運(yùn)方程 (4) ，方程在每個時間步求解時，平板上的波浪力表示為：

(6)

式中：F表示波浪力的矢量；uτ表示切線速度分量；n表示單位法向矢量。

2 數(shù)值水槽

圖1展示了本研究中使用的2-D數(shù)值波浪水槽的示意圖，坐標(biāo)系的定義、浪高儀的位置不按比例。波浪水槽長9 m，高0.292 m，寬(w) 0.149 m，坐標(biāo)系原點(diǎn)位于左側(cè)松弛消波區(qū)邊界的端點(diǎn)位置。x軸定義為波浪傳播的方向，z軸定義為垂直向上的方向，y方向的水槽厚度對應(yīng)于一個單元。水平板放在波浪水槽的中間,板長L=0.149 m，寬B=0.305 m，高d=0.012 7 m。水深定義為h=0.086 m，波高定義為A=0.025 m，ζ為自由面高程。

圖1 數(shù)值波浪水槽Fig.1 Sketch of numerical wave tank

實(shí)驗(yàn)進(jìn)行了5組模擬，水平板板寬分別設(shè)置為B=0.305，0.400，0.500，0.600，0.800 m，每組的數(shù)值模擬中都考慮了6種相對凈空比z*/h= 0，0.05，0.1，0.139，0.213，0.287(z*表示凈空高度，即水平板的下表面到靜水面的距離，z*/h=0時板下表面與靜水面相平齊，z*/h= 0.139時板的上表面與波峰相平齊，z*/h= 0.287時板的下表面與波峰相平齊)，共30個算例。設(shè)置了3個浪高儀，用G1～G3用來記錄自由面高程，分別位于x=2.010，3.535，4.75 m處。在波浪水槽的入出口邊界分別設(shè)置了2、2.245 m長的松弛區(qū)吸收反射波和入射波。

使用OpenFOAM?提供的內(nèi)置網(wǎng)格生成程序“blockMesh”生成網(wǎng)格，B=0.305 m，z*/h=0工況下一個典型的計算網(wǎng)格如圖2所示，平板四周的網(wǎng)格區(qū)域分別對應(yīng)水、空氣兩相環(huán)境。本數(shù)值實(shí)驗(yàn)采用非均勻網(wǎng)格，為了節(jié)省計算時間和準(zhǔn)確地捕捉水和空氣之間的界面，在平板周圍使用分辨率更高的細(xì)網(wǎng)格。

圖2 計算域中典型網(wǎng)格的側(cè)視圖Fig.2 Side view of typical meshes in the computational domain

為了檢驗(yàn)數(shù)值結(jié)果對網(wǎng)格分辨率的準(zhǔn)確性，使用了3種不同分辨率的網(wǎng)格(即粗網(wǎng)格、中網(wǎng)格和細(xì)網(wǎng)格)進(jìn)行數(shù)值模擬結(jié)果對比，網(wǎng)格單元數(shù)分別為164 510、370 215和658 040(見表1)。計算的參數(shù)：h=0.086 m，A=0.025 m，B=0.305 m，z*/h=0.1。圖3顯示了入射波作用于平板引起的G1處自由面變化和波浪力的時間歷程情況。

表1 網(wǎng)格參數(shù)Table 1 Mesh parameters

由圖3可以看出，3種不同網(wǎng)格分辨率下得到的自由表面高程、水平波浪力和垂直波浪力的時間歷程曲線幾乎吻合?？紤]到與粗網(wǎng)格、細(xì)網(wǎng)格相比，中網(wǎng)格能更加精確地模擬入射波所引的波場和節(jié)省數(shù)據(jù)計算的時間，所以接下來的算例中均采用中網(wǎng)格的網(wǎng)格配置進(jìn)行計算。

圖3 h=0.086 m，A=0.025 m，B=0.305 m，z*/h=0.1,不同網(wǎng)格分辨率下時間歷程Fig.3 h=0.086 m，A=0.025 m，B=0.305 m，z*/h= 0.1 time history at different grid resolutions

3 數(shù)值驗(yàn)證

為了驗(yàn)證本文數(shù)值模型的可靠性以及得到的數(shù)值結(jié)果的準(zhǔn)確性，將OpenFOAM?得到的計算結(jié)果與Seiffert等[13]的實(shí)驗(yàn)結(jié)果進(jìn)行了比較。其實(shí)驗(yàn)?zāi)Ｐ团c本文的數(shù)值模型大致相同(見圖1)，板長L=0.149 m，寬B=0.305 m，高d=0.012 7 m，水深h=0.086 m，波高A=0.025 m，改變了水平板的位置，Seiffert等[13]實(shí)驗(yàn)?zāi)Ｐ偷乃桨逦挥谒聑/h=0.2處(z/h表示水平板的淹沒比，z表示水平板的上表面距靜水面的垂直距離)。

通過圖4可以看到，通過OpenFOAM?模擬得到的G1～G3處的自由面變化情況、水平波浪力Fx和垂直波浪力Fz的時間歷程情況與Seiffert等[13]得到的實(shí)驗(yàn)結(jié)果大致吻合。驗(yàn)證了本文數(shù)值模型對模擬孤立波作用于固定水平板得到的波面變化情況和提取波浪力的準(zhǔn)確性。

圖4 OpenFOAM?計算和實(shí)驗(yàn)數(shù)據(jù)得到的自由面高程、波浪力的比較Fig.4 Comparison of free surface elevation and wave force from OpenFOAM? calculation and experimental data

4 算例分析與討論

4.1 波浪力的時間歷程

圖5展示了0.305≤B≤0.8 m, 0≤z*/h≤0.287工況下，孤立波作用于水平板上的水平波浪力和垂向波浪力的時間歷程。對于水平波浪力Fx的時間序列(圖5(a))，圖中的正值和負(fù)值分別指沿x軸正向和負(fù)向，隨著板寬、相對凈空比的增加,水平波浪力在曲線形狀上沒有較為明顯的差異,隨著孤立波的作用,波浪力先增大后減小，后期由于波浪完整越過水平板，波板之間相互作用復(fù)雜，殘余波的振蕩引起板的力振蕩。對于垂直波浪力Fz，隨著相對凈空比的增加，垂向波浪力的時間歷程越來越復(fù)雜，垂向正向力明顯大于反向力。從圖5(b)可以清楚地看到，z*/h=0.05工況下的垂向波浪力Fz較為明顯地大于其他相對凈空比狀態(tài)下的力，結(jié)合圖6可知是由于其凈空高度的原因，z*/h=0.05工況下的平板高于靜水面，其下表面亦能全部受到波浪的作用，而當(dāng)平板繼續(xù)升高，相對凈空比z*/h=0.1時，由于平板過高，波浪在上升的過程中發(fā)生破碎，能量逐漸減小，導(dǎo)致板的下表面不能全部受到波浪作用，當(dāng)平板凈空繼續(xù)增大時，垂向波浪力變得更?。黄浯瓜虿ɡ肆σ嗍馨鍖挼挠绊?，隨著板寬的增加而增大。針對垂向波浪力Fz復(fù)雜的時間變化序列，下面將根據(jù)B=0.305 m，0≤z*/h≤0.213工況下的波浪力時間歷程圖(見圖7)進(jìn)行具體描述。

圖5 B=0.305～0.800 m，z*/h=0～0.287工況下的波浪力時間歷程圖Fig.5 Time series of wave forces with B=0.305～0.800 m, z*/h=0～0.287

圖6 B=0.8 m，z*/h=0，0.05和0.1垂向波浪力到達(dá)最大值時板周圍的自由面情況Fig.6 Free surface around the plate when B=0.8 m, z*/h=0, 0.05 and 0.1 vertical wave force reach the maximum value

由圖7(a) 可以看出，在整個階段垂向波浪力只出現(xiàn)了一次正向極值，并隨著孤立波的運(yùn)動先增大后減小，后期由于殘余波振蕩的緣故，導(dǎo)致出現(xiàn)了時間較長的波浪力震蕩現(xiàn)象。通過圖7可知，隨著平板的相對凈空比增大，垂向波浪力的時間歷程階段開始發(fā)生了較為明顯的改變。圖7 (b)顯示了平板開始離開靜水面，在孤立波作用于平板的過程中發(fā)生了2次砰擊，首先是當(dāng)孤立波的波峰接近平板左側(cè)時波浪對平板下端的沖擊作用得到了一個極值，之后波浪主體發(fā)生破碎，一部分波浪在平板下方繼續(xù)沖擊向前，并向平板的下表面作用垂直向上的力，另一部分波浪受到平板左側(cè)部分阻擋向上激起，落下后作用于平板的上表面，結(jié)合平板下表面所受的垂直向上的波浪力，再次得到了一個極值，從而得到了垂向波浪力的時間歷程上具有2個正向極值的特點(diǎn)。這種垂向波浪力的變化情況也被Ren等[16]、Cuomo[17]提到過。隨著平板的相對凈空比繼續(xù)增加(圖7(c))，當(dāng)發(fā)生砰擊取得第1次極值之后，波浪發(fā)生破碎，能量減少，不足以支撐板下的波浪繼續(xù)作用于板的下表面，波浪逐漸削減，垂向波浪力逐漸減小，不能得到較為明顯的第2次垂向極值波浪力。當(dāng)平板的相對凈空比繼續(xù)增加，平板所受的垂向波浪力就變得更小(圖7(d))。

圖7 B=0.305 m，z*/h=0～0.13垂向波浪力的時間歷程圖Fig.7 Time series of vertical wave force with B=0.305 m, z*/h=0～0.213

4.2 極值波浪力

圖8 不同板寬條件下，正、負(fù)水平極值波浪力相對于相對凈空比的變化Fig.8 Variation of positive and negative horizontal extreme wave forces with respect to relative clearance ratio under different plate width conditions

圖9顯示了多組板寬條件下，正、負(fù)垂向極值波浪力相對于相對凈空比的變化情況。由圖9 (a)可得知，隨著相對凈空比的增加，正垂向極值波浪力呈先增大后減小的趨勢，0.5≤B≤0.8，z*/h=0.05工況下的平板，其正垂向極值波浪力達(dá)到最大，為整體變化趨勢的轉(zhuǎn)折點(diǎn)，并且隨著板寬的增加這種現(xiàn)象更加明顯，結(jié)合圖5和圖6自由面高程圖了解到是由于凈空高度和板寬的原因?qū)е?。?dāng)相對凈空比一定時，正垂向極值波浪力隨著板寬的增加而逐漸增加；當(dāng)相對凈空比大于0.139時(即波峰處于平板下表面到上表面之間)，板寬變化對正垂向極值波浪力的影響逐漸減弱最后趨于集中，此情況與對正水平極值波浪力的影響一致(見圖8(a))。隨著相對凈空比的增加，負(fù)垂向極值波浪力亦呈先增大后減小的趨勢(見圖9(b))。比較圖9 (a) 和 (b) 可以發(fā)現(xiàn)正垂向極值波浪力均大于負(fù)極值波浪力，主要是因?yàn)椴ɡ嗽谧饔糜谄桨鍟r，一部分波浪向上涌起發(fā)生破碎，導(dǎo)致負(fù)向的極值波浪力小于正向的極值波浪力。

圖9 不同板寬條件下，正、負(fù)垂向極值波浪力相對于相對凈空比的變化Fig.9 Variation of wave forces with respect to the relative clearance ratio of positive and negative vertical extremum under different plate width conditions

4.3 透射、反射系數(shù)

基于Ning等[18]中關(guān)于單孤立波的運(yùn)動描述，反射波的形狀由一個波峰和一個波谷組成，通過G1處的自由表面高程可獲得入射波、反射波的波高(見圖1，圖10(a))，反射波的波高定義為波峰到波谷的垂直距離，反射系數(shù)Cr定義計算為反射波高與入射波高A的比值。透射波高可由G2處的自由面高程求得(見圖10(b))，透射系數(shù)Ct定義計算為透射波高與入射波高A的比值。相對凈空比和板寬對反射系數(shù)Cr、透射系數(shù)Ct的影響如圖11所示。

圖10 G1，G2的自由面變化歷程Fig.10 Free surface change course of G1 and G2

通過比較圖11 (a) 和 (b) 可以發(fā)現(xiàn)，對于本文考慮的所有板寬參數(shù)，反射系數(shù)、透射系數(shù)隨著相對凈空比的增加具有不同的變化趨勢。對于反射系數(shù)(圖11(a))，隨著相對凈空比的增加，不同板寬下的Cr呈逐漸減小的趨勢；當(dāng)相對凈空比一定時，Cr隨著板寬的增加而增加，且當(dāng)相對凈空比逐漸增大時，其隨板寬增加對反射系數(shù)的影響逐漸減弱。相對凈空比和板寬參數(shù)對反射系數(shù)變化特性的影響與對水平波浪極值力的影響非常相似(圖8(a))。對于透射系數(shù)(圖11(b))，Cr與Ct的變化趨勢相反，在不同的板寬條件下，透射系數(shù)隨著相對凈空比的增加皆呈單調(diào)遞增趨勢；當(dāng)相對凈空比一定時，透射系數(shù)隨著板寬的增加而逐漸減小。

圖11 不同板寬條件下，反射系數(shù)、透射系數(shù)相對于凈空比的變化Fig.11 Variation of reflection coefficient and transmission coefficient with respect to clearance ratio under different plate width conditions

5 結(jié)論

1) 水平板板寬和相對凈空比的變化對水平波浪力的變化趨勢影響不大，對垂向波浪力的影響較為明顯，且隨著相對凈空比的增加，其影響變得更加顯著。

2) 隨著相對凈空比的增加，平板的水平極值波浪力逐漸減小，垂向極值波浪力先增大后減??；板寬變化對極值波浪力的影響隨著相對凈空比的增加逐漸減弱。

3) 透射、反射系數(shù)受相對凈空比的影響，隨著相對凈空比的增加反射系數(shù)逐漸減小，透射系數(shù)逐漸增大，并且當(dāng)平板的相對凈空比處于0.139～0.287時(即波峰高度處于板的上下表面高度之間)，板寬變化對反射系數(shù)的影響逐漸減弱，反射系數(shù)曲線逐漸趨于集中；當(dāng)相對凈空比一定時，隨著板寬的增加，反射系數(shù)逐漸增大,透射系數(shù)逐漸減小。