孫亮， 那鑫宇， 陳柯宇， 李芬， 陳明勝

(1.武漢理工大學(xué) 交通學(xué)院，湖北 武漢 430063； 2.大連理工大學(xué) 海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024； 3.河海大學(xué) 港口海岸與近海工程學(xué)院，江蘇 南京 210098)

波浪是設(shè)計(jì)近岸結(jié)構(gòu)物需要考慮的主要環(huán)境因素，波浪從深水傳播到近岸區(qū)域會(huì)由于水下結(jié)構(gòu)物或者地形的影響而發(fā)生形態(tài)的改變。潛堤是一種建造在海岸線外、完全淹沒(méi)在水面下的護(hù)岸建筑物[1]，其主要功能是阻擋和減弱外海波浪的作用，為港口、防波堤、閘門等提供掩護(hù)作用[2]。波浪在潛堤上傳播時(shí)，不同頻率的波浪會(huì)發(fā)生相互作用，產(chǎn)生高階諧波，也可能發(fā)生破碎等現(xiàn)象。分析波浪在潛堤和潛堤地形上的傳播變形的方法主要有理論分析、物理模型試驗(yàn)和數(shù)值模擬3種研究方法。理論分析方法對(duì)試驗(yàn)研究及數(shù)值模擬具有指導(dǎo)意義，但在解決實(shí)際工程問(wèn)題時(shí)具有較大的局限性。物理模型試驗(yàn)可通過(guò)試驗(yàn)波浪水槽進(jìn)行準(zhǔn)確的現(xiàn)場(chǎng)觀測(cè)，得到較為可靠的實(shí)驗(yàn)數(shù)據(jù)，但成本偏高，流程復(fù)雜。隨著數(shù)值分析相關(guān)理論的發(fā)展和計(jì)算機(jī)硬件性能的提高，數(shù)值模擬方法成本低，重復(fù)性強(qiáng)，相比物理模型試驗(yàn)可以節(jié)約大量人力、資源和時(shí)間。

分析近岸波浪傳播變形的數(shù)學(xué)模型主要包括：Boussinesq方程、緩坡方程、N-S方程等。Beji等[3]基于一維的Boussinesq方程模擬了非破碎波在潛提上的傳播并與試驗(yàn)結(jié)果相比較。Losada等[4]使用緩坡方程模擬了可滲透梯形潛堤上的波浪演變與附近流場(chǎng)情況。滕愛(ài)國(guó)等[5]以完整N-S方程為基礎(chǔ)，改進(jìn)湍流模型，較為精確地模擬出橢圓余弦波在潛堤上的傳播過(guò)程。周援衡等[6]同樣地基于N-S方程，采用VOF法追蹤自由液面和湍流模型，對(duì)橢圓余弦波在滲透潛堤上的傳播變形以及潛堤前后的波形變化進(jìn)行數(shù)值模擬。劉誠(chéng)等[7]建立了比前人更為精確的PLIC-VOF數(shù)值方法并使用有限差分法求解N-S方程，研究了波浪越過(guò)斜坡潛堤和垂直潛堤時(shí)的波浪場(chǎng)。張婷婷[8]利用Boussinesq方程模擬了規(guī)則波和自由長(zhǎng)波在潛堤地形上的傳播變形。吳迪[9]采用Boussinesq方程，模擬了波浪在非滲透潛堤和滲透潛堤上的傳播并進(jìn)行對(duì)比，顯示出滲透作用對(duì)波浪的巨大消波作用。馬小舟等[10]利用完全非線性Boussinesq方程數(shù)值模型研究了不規(guī)則波浪在緩坡潛堤上傳播時(shí)非線性參數(shù)的變化。廖斌等[11]利用五階WENO有限差分格式以及三階TVD Runge-Kutta格式離散求解RANS方程分析了潛堤上波浪的傳播問(wèn)題。

本文基于開(kāi)源計(jì)算流體力學(xué)工具OpenFOAM建立了二維數(shù)值波浪水槽[12]，模擬規(guī)則波和橢圓余弦波在潛堤及潛堤地形上的傳播變形并與已發(fā)表的結(jié)果對(duì)比驗(yàn)證，為進(jìn)一步分析波浪對(duì)海岸結(jié)構(gòu)物的作用提供了堅(jiān)實(shí)基礎(chǔ)。

1 數(shù)值波浪水槽的建立

OpenFOAM利用有限體積法求解N-S方程分析可壓縮以及不可壓縮流體問(wèn)題。對(duì)于不可壓縮流體，N-S方程由連續(xù)性方程和動(dòng)量守恒方程組成[13]：

·U=0

(1)

?ρU/?t+·(ρUU)-·(μU)-ρg=

(2)

式中：ρ為流體密度；p為流體壓力；μ為粘滯系數(shù)；U為速度矢量；fσ為表面張力。

在建立數(shù)值波浪水槽的左端邊界，根據(jù)線性波浪理論計(jì)算波面、水質(zhì)點(diǎn)速度和壓力分布從而產(chǎn)生入射波浪。在左端入射邊界和右端邊界都設(shè)置了相應(yīng)的消波海綿層[14]。

2 規(guī)則波在潛堤上的傳播

2.1 物理模型試驗(yàn)及數(shù)值波浪水槽的布置

為了驗(yàn)證數(shù)值模型的精度，本文首先使用Gobbi等[15]研究波浪在淹沒(méi)潛堤上的傳播問(wèn)題時(shí)所采用的代爾伏特物理模型試驗(yàn)數(shù)據(jù)，入射波浪的參數(shù)如表1所示，潛堤的形狀和本文對(duì)比中所選取的浪高儀位置如圖1所示。

表1 代爾伏特物理模型試驗(yàn)的入射波浪參數(shù)

圖1 代爾伏特物理模型試驗(yàn)中潛堤與浪高儀的位置(單位：m)

本文建立了如圖2所示的二維數(shù)值波浪水槽，長(zhǎng)30 m，高0.5 m，水深設(shè)為0.4 m。數(shù)值波浪水槽的左端為造波邊界，波浪由左向右傳播，右端為出流邊界。在造波邊界后設(shè)置1.0 m長(zhǎng)的海綿層，避免由地形產(chǎn)生的反射波到達(dá)左端造波邊界。在右端出流邊界前設(shè)置6.0 m長(zhǎng)的海綿層，避免越過(guò)潛堤的波浪到達(dá)右端邊界在數(shù)值水槽內(nèi)形成反射。對(duì)于如圖2所示的計(jì)算區(qū)域，沿水槽寬度方向只劃分了1個(gè)網(wǎng)格。在水槽不同垂直截面處劃分了相同的單元數(shù)，這樣在潛堤斜坡上和坡頂網(wǎng)格會(huì)自然加密，有利于精確模擬潛提附近區(qū)域復(fù)雜的物理過(guò)程。初始網(wǎng)格的控制尺寸取為0.01 m，每個(gè)垂直截面上有50個(gè)單元，水平方向有3 000個(gè)單元，網(wǎng)格1在計(jì)算區(qū)域的總單元數(shù)為15萬(wàn)。網(wǎng)格2在每個(gè)垂直截面上有100個(gè)單元，計(jì)算區(qū)域單元總數(shù)30萬(wàn)。網(wǎng)格3在每個(gè)垂直截面上有200個(gè)單元，計(jì)算區(qū)域單元總數(shù)為60萬(wàn)。

2.2 波面驗(yàn)證與精度分析

圖3給出了10個(gè)斷面處(如圖3所示)的波面隨時(shí)間變化的計(jì)算和試驗(yàn)測(cè)量值的比較。圖3(a)和(b)中x=2.0 m和x=4.0 m處的結(jié)果用于研究入射波浪遇到潛堤之前的變化情況，由圖3(a)和(b)中的對(duì)比可以看出基于3種網(wǎng)格的數(shù)值模擬結(jié)果和測(cè)量值吻合較好，數(shù)值波浪水槽可以較好地重復(fù)試驗(yàn)中的入射波。

圖2 數(shù)值波浪水槽的設(shè)置(單位：m)

圖3 代爾伏特試驗(yàn)中波面測(cè)量值與數(shù)值結(jié)果

圖3(c)～(g)中斷面x=10.5，12.5，13.5，14.5，15.7 m位于潛堤斜坡上和坡頂。在x=10.5 m(位于潛堤前坡)，x=12.5和13.5 m(位于潛堤頂)處基于3種網(wǎng)格的數(shù)值結(jié)果幾乎沒(méi)有區(qū)別并且與測(cè)量值吻合得較好，表明目前的3種單元?jiǎng)澐侄紳M足了精度要求。在x=14.5，15.7 m(位于潛堤后坡)處基于3種網(wǎng)格的數(shù)值結(jié)果出現(xiàn)了差別，網(wǎng)格2和3的結(jié)果更加接近測(cè)量值。正如Gobbi等[15]所分析所得結(jié)論：波浪經(jīng)過(guò)潛堤后高階諧波發(fā)生分離使得波浪場(chǎng)變得更加復(fù)雜。所以更加精細(xì)的網(wǎng)格才能滿足計(jì)算精度的要求。圖3(h)～(j)中斷面x=17.3，19.0，21.0 m完全位于潛堤后，基于網(wǎng)格3的數(shù)值結(jié)果與測(cè)量值吻合得更好。

3 橢圓余弦波在潛堤地形上的傳播

3.1 數(shù)值波浪水槽的布置

張婷婷[8]進(jìn)行了規(guī)則波在潛提地形上的試驗(yàn)研究,潛堤地形由橢圓型曲線和指數(shù)型曲線組成，其數(shù)學(xué)表達(dá)式為：

z+h=

(3)

式中：h是常水深h=0.5 m。hs是潛堤最高點(diǎn)的高度hs=0.4 m，a1=4.0 m，a2=2.5 m，θ1=arctan(a1z1/hs)，z1=0.125，θ2=arctan(a2z2/hs)，z2=0.2。如圖4所示曲線地形分為4個(gè)部分，其中x0=20.2 m，為潛堤最高點(diǎn)對(duì)應(yīng)的橫坐標(biāo)，此處z+h=hs=0.4 m。x1、x2為潛堤最大坡度處對(duì)應(yīng)的橫坐標(biāo)，x1=17.077 m，x2=22.152 m。張婷婷[8]將基于Boussinesq方程的數(shù)值結(jié)果與模型試驗(yàn)數(shù)據(jù)比較，驗(yàn)證了數(shù)值結(jié)果的可靠性。

張婷婷[8]進(jìn)一步利用Boussinesq方程模擬了橢圓余弦波在潛堤地形上的傳播變形，分析了如圖4所示x=10.0，16.0，18.0，19.2，20.2，22.5，24.8 m處的波面變化。本研究中所采用的橢圓余弦波的波浪參數(shù)如表2所示。

圖4 大連理工大學(xué)試驗(yàn)水槽中的潛堤地形

表2 橢圓余弦波的波浪參數(shù)

本文中建立了與圖2類似的總長(zhǎng)130 m、高0.6 m的二維數(shù)值波浪水槽，左端吸波海綿層長(zhǎng)40 m，右端吸波海綿層長(zhǎng)45 m。二維數(shù)值波浪水槽中潛堤地形頂部的水平坐標(biāo)由圖4中的20.2 m變?yōu)?5.2 m，數(shù)值浪高儀的位置變?yōu)椋簒=55.0，61.0，63.0，64.2，65.2，67.5，69.8 m。計(jì)算區(qū)域中每個(gè)垂直截面上有300個(gè)單元，水平方向有1 300個(gè)單元，共計(jì)39萬(wàn)個(gè)單元。

3.2 波面驗(yàn)證與精度分析

圖5給出了數(shù)值波浪水槽中7個(gè)位置處基于Boussinesq方程和N-S方程的波面計(jì)算結(jié)果。

圖5 橢圓余弦波在潛堤地形的傳播過(guò)程

圖5(a)中x=55.0 m(對(duì)應(yīng)圖4中x=10.0 m)位于潛堤地形前，用于研究入射波浪遇到潛堤之前的變化情況；圖5(b)中斷面x=61.0 m(對(duì)應(yīng)圖4中x=16.0 m)，圖5(c)中x=63.0 m(對(duì)應(yīng)圖4中x=18.0 m)，圖5(d)中x=64.2 m(對(duì)應(yīng)圖4中x=19.2 m)位于潛堤前坡處，由圖5(a)～(d)中的對(duì)比可以看出：基于Boussinesq的數(shù)值結(jié)果和基于N-S方程的數(shù)值結(jié)果吻合較好。基于Boussinesq的數(shù)值結(jié)果和基于N-S方程的數(shù)值結(jié)果吻合較好。圖5(e)中x=65.2 m(對(duì)應(yīng)圖4中x=20.2 m)位于潛堤坡頂，圖5(f)中x=67.5 m(對(duì)應(yīng)圖4中x=22.5 m)和圖5(g)中x=69.8 m(對(duì)應(yīng)圖4中x=24.8 m)位于潛堤后，OpenFOAM數(shù)值模擬結(jié)果和Boussinesq方程數(shù)值模擬結(jié)果出現(xiàn)微小差別。

值得注意的是在圖5(g)中斷面x=69.8 m(對(duì)應(yīng)圖4中x=24.8 m)處，在基于Boussinesq方程的數(shù)值模擬結(jié)果中發(fā)現(xiàn)了毛刺波[8]。張婷婷[8]認(rèn)為：自由長(zhǎng)波本身含有高階諧波，波浪越過(guò)潛堤后發(fā)生反淺水效應(yīng)，原來(lái)被約束的諧波發(fā)生分離成為自由波。然而基于OpenFOAM的數(shù)值模擬結(jié)果中并沒(méi)有體現(xiàn)這一現(xiàn)象，需要進(jìn)一步分析。

4 結(jié)論

1)本文建立了基于OpenFOAM的二維數(shù)值波浪水槽，分別模擬了規(guī)則波在潛堤上的傳播、橢圓余弦波在潛堤地形上的傳播，所得的波面結(jié)果與試驗(yàn)數(shù)據(jù)、基于Boussinesq方程的數(shù)值結(jié)果都吻合較為理想。

2)本文的研究驗(yàn)證了基于OpenFOAM的數(shù)值波浪水槽在預(yù)測(cè)波浪與近岸結(jié)構(gòu)物相互作用時(shí)的可靠性，為近岸結(jié)構(gòu)物的相關(guān)研究及設(shè)計(jì)提供可靠依據(jù)。