牟迪, 王榮泉, 寧德志, 陳麗芬
(大連理工大學(xué) 海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,遼寧 大連 116024)
海嘯常指地球地殼運(yùn)動(dòng)引起的海浪。通常引起海嘯的原因包括但不限于:地震、風(fēng)暴潮、水下爆炸等。當(dāng)以上事件發(fā)生時(shí),水面的突然抬升或者降低將產(chǎn)生沿各個(gè)方向傳播的波浪,并產(chǎn)生海嘯。結(jié)構(gòu)物在海嘯波作用下遭受猛烈的極端載荷并引起流體諸多強(qiáng)非線性現(xiàn)象。全球在過去的100年中,每年大約發(fā)生一次毀滅性海嘯,在世界范圍內(nèi)造成嚴(yán)重的人員損失與財(cái)產(chǎn)破壞[1]。海嘯在產(chǎn)生時(shí),波幅較小,通常僅有1~2 m,但與當(dāng)?shù)厮钕啾?波長(zhǎng)很長(zhǎng),通常在100 km左右。孤立波只有一個(gè)較大的波峰或波谷,其波形和水質(zhì)點(diǎn)的運(yùn)動(dòng)與海嘯極為相似,因此常被用來(lái)模擬海嘯等淺水大波[2]。
從大量海洋災(zāi)害所造成破壞的案例看來(lái),海嘯引發(fā)的越浪現(xiàn)象,是造成海洋結(jié)構(gòu)物破壞及財(cái)產(chǎn)損失的重要原因之一。自20世紀(jì)起,國(guó)內(nèi)外眾多學(xué)者已經(jīng)對(duì)于越浪問題進(jìn)行了大量的研究并取得了大量研究成果[3-5],但以往的研究多集中于海嘯波與不同結(jié)構(gòu)形式的海堤或者近岸跨海橋梁等海岸結(jié)構(gòu)物作用時(shí)產(chǎn)生的越浪現(xiàn)象,而關(guān)于孤立波與振蕩水柱式(oscillation water column, OWC)波浪能轉(zhuǎn)換裝置作用下引起的越浪現(xiàn)象研究較少。
OWC波能裝置的主體帶有一個(gè)下端開口的氣室結(jié)構(gòu),氣室上部有通氣管道連接空氣透平。在波浪作用下,水柱上下振蕩迫使氣室內(nèi)部空氣在通氣管道處形成往復(fù)氣流,氣流通過推動(dòng)空氣透平帶動(dòng)發(fā)電機(jī)工作產(chǎn)生電能。其與普通海洋結(jié)構(gòu)物不同之處主要在于結(jié)構(gòu)內(nèi)部氣室的存在,氣室內(nèi)部上方氣流由于氣孔的阻礙作用使氣壓增大,進(jìn)而抑制氣室內(nèi)部自由液面的上升。被抑制的水柱通過壓強(qiáng)的傳遞進(jìn)而影響波浪在前墻上的爬坡運(yùn)動(dòng)最終對(duì)越浪產(chǎn)生影響。以往對(duì)該裝置的研究多從水動(dòng)力轉(zhuǎn)換效率角度出發(fā),優(yōu)化該裝置結(jié)構(gòu),提升其波能轉(zhuǎn)換能力[6-8]。而對(duì)于該裝置在極端波況下的水動(dòng)力特性,國(guó)內(nèi)外大多學(xué)者關(guān)注點(diǎn)集中在裝置的受力情況,如馬子然[9]研究了在聚焦波作用下,離岸式OWC裝置的受力情況。盡管如此,目前未見有關(guān)于極端波浪作用下,OWC波能裝置越浪特性的研究。
所以本文基于以上背景,利用開源軟件OpenFOAM 中waveIsoFoam求解器求解三維Navier-Stokes方程,對(duì)孤立波與OWC裝置作用時(shí)發(fā)生的越浪現(xiàn)象進(jìn)行數(shù)值模擬研究,分析與越浪相關(guān)的物理量如越浪量、越浪厚度、越浪流速度等隨相對(duì)波高改變的變化規(guī)律,并展示越浪發(fā)生時(shí)的流場(chǎng)特性。為OWC波能裝置的優(yōu)化及安全設(shè)計(jì)提供理論指導(dǎo)。
本文假定流體不可壓,無(wú)熱傳導(dǎo)。基于上述假定,N-S方程可分為連續(xù)性方程及動(dòng)量方程:

(1)
(2)
式中ρ、μ表示混合流體密度及流體粘度。u=(u,v,w)和x=(x,y,z)為笛卡爾坐標(biāo)系下流體速度矢量與位置坐標(biāo),p*為動(dòng)壓,g為重力加速度。
本文通過VOF方法獲得自由液面運(yùn)動(dòng)。該方法求解運(yùn)輸方程(3)以獲得體積分?jǐn)?shù):
(3)
式中:α為每個(gè)網(wǎng)格體心的體積分?jǐn)?shù),如果α的值界于0~1,則這個(gè)網(wǎng)格則被標(biāo)記為包含自由液面的計(jì)算網(wǎng)格。當(dāng)α=0時(shí),對(duì)應(yīng)為氣相;當(dāng)α=1時(shí),對(duì)應(yīng)水相。單元內(nèi)流體物理屬性為:
(4)
式中:角標(biāo)a與w分別代表空氣與水,盡管基于不可壓縮假定,但是交界面處密度隨計(jì)算時(shí)間步推進(jìn)而時(shí)刻變化。
在本文中,用于重構(gòu)自由液面的VOF方法為Iso-advector 方法,該方法是一種基于幾何特征的重構(gòu)方法。相比于基于代數(shù)重構(gòu)的MULES方法,該方法在計(jì)算效率、收斂性、穩(wěn)定性有著較強(qiáng)的優(yōu)勢(shì),關(guān)于該方法詳細(xì)描述見文獻(xiàn)[10-11]。
本文采用waves2Foam中提供的一階孤立波波面方程進(jìn)行數(shù)值造波,并在計(jì)算域兩側(cè)分別布置被動(dòng)消波區(qū)域。其中孤立波波面方程為:
(5)
(6)
(7)
(8)
式中:η為孤立波波面;H為孤立波波高;d為水深;c為孤立波波速。
本文采用有限體積法將計(jì)算域離散成眾多結(jié)構(gòu)化網(wǎng)格。對(duì)于方程(2)中時(shí)間項(xiàng)采用一階Euler隱式離散,對(duì)流項(xiàng)采用linearUpwind 二階TVD格式離散,在保證計(jì)算精度的同時(shí)保證其求解的穩(wěn)定性,速度梯度通過Gauss linear interpolation 方法進(jìn)行插值計(jì)算。速度壓力耦合通過pimple算法求解,pimple算法為simple算法與piso算法的結(jié)合,允許用戶根據(jù)所設(shè)定的最大庫(kù)郎數(shù)進(jìn)行時(shí)間步的自適應(yīng)調(diào)節(jié)[12],本文計(jì)算中庫(kù)郎數(shù)為0.3。
本文中所采用的邊界條件如下:對(duì)于所有結(jié)構(gòu)物物面及底部地面,采用noslip無(wú)滑移邊界條件并配合zeroGradient 零梯度壓力邊界條件;對(duì)于計(jì)算域頂部采用totalPressure及pressureInletOutlet邊界,即在頂部指定總壓,并僅允許空氣流出,不允許回流。
為排除空間離散對(duì)數(shù)值結(jié)果產(chǎn)生的影響,本節(jié)將進(jìn)行網(wǎng)格收斂性驗(yàn)證。表1列出了不同網(wǎng)格參數(shù),結(jié)構(gòu)物外部水平方向網(wǎng)格與垂向網(wǎng)格尺寸相等。由于氣孔處通量變化較大,對(duì)其附近區(qū)域進(jìn)行了局部加密(見圖1)。數(shù)值模型中計(jì)算域尺寸與物模實(shí)驗(yàn)中保持一致(見下節(jié))。圖2展示了不同網(wǎng)格下由OpenFOAM模擬得到的x=9.3 m處波面歷程,本文模擬中水深保持d=0.4 m。圖2中,不同網(wǎng)格得到的結(jié)果差距較小,為平衡數(shù)值結(jié)果與計(jì)算量,最終選擇中等程度的網(wǎng)格進(jìn)行計(jì)算。

表1 網(wǎng)格尺度

圖1 網(wǎng)格示意

圖2 H*=H/d=0.425時(shí)不同網(wǎng)格下波面歷程(η*=η/d,t*=t(g/d)1/2)
試驗(yàn)在大連理工大學(xué)海岸及近海工程實(shí)驗(yàn)室PIV水槽進(jìn)行,水槽長(zhǎng)22 m、寬0.45 m、高0.6 m。如圖3所示,模型由厚度為0.012 m的有機(jī)玻璃制成,裝置前墻固定于距造波板9.3 m處,以保證孤立波在傳播過程進(jìn)行充分演化。裝置的具體幾何尺寸見表2。由于水槽寬度較小,并且在孤立波作用時(shí),自由液面伴隨著強(qiáng)非線性破碎現(xiàn)象,難以準(zhǔn)確測(cè)量距前墻較近處自由液面變化,故無(wú)接觸式超聲波浪高儀安置在距裝置前墻2.6 m處。為避免測(cè)力計(jì)應(yīng)變片微小運(yùn)動(dòng)而引起受力的不準(zhǔn)確測(cè)量,2個(gè)測(cè)力計(jì)分別安裝于裝置頂部與后方底部并通過一鋼架固定于水槽側(cè)壁。實(shí)驗(yàn)結(jié)果表明,通過雙測(cè)力計(jì)測(cè)量得到的總力,其結(jié)果誤差小于3%。

表2 結(jié)構(gòu)物尺寸

圖3 本次實(shí)驗(yàn)示意
圖4分別展示數(shù)值模擬結(jié)果與模型實(shí)驗(yàn)結(jié)果的對(duì)比。
圖4(a)第1個(gè)峰值為入射波,第2個(gè)峰值為反射波。OpenFOAM計(jì)算得到的反射波峰值略大于實(shí)驗(yàn)結(jié)果,其原因?yàn)閷?shí)驗(yàn)水槽側(cè)壁高度受限,波浪與前墻作用時(shí)部分水體由水槽兩側(cè)溢出而導(dǎo)致能量部分損失,最終使反射波波高有所降低。圖4(b)展示氣室內(nèi)氣壓對(duì)比,其中正壓對(duì)應(yīng)于OWC運(yùn)動(dòng)的上升階段,當(dāng)氣室內(nèi)液面高度達(dá)到最大時(shí),氣室內(nèi)氣壓降至零。此后,OWC開始做下降運(yùn)動(dòng),氣室內(nèi)部產(chǎn)生負(fù)壓,外部空氣被吸入至氣室內(nèi)部。圖4(c)展示結(jié)構(gòu)總水平力對(duì)比,受力曲線呈現(xiàn)單峰狀,受力峰值與實(shí)驗(yàn)結(jié)果誤差小于2%。總體來(lái)說(shuō),數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比良好,驗(yàn)證了OpenFOAM模型的有效性。
采用Baldock等[13]提出的方法,選取孤立波在靜水面以上的單寬水體體積作為特征量定義無(wú)量綱越浪量即:
q*=q/q0
式中:q為數(shù)值模擬得到的單寬越浪量;q0為孤立波在靜水面以上的單寬水體體積,對(duì)方程(5)積分可得:
(9)
圖5給出無(wú)量綱越浪量與相對(duì)波高的關(guān)系,所監(jiān)測(cè)斷面為前墻上方,其中定義相對(duì)波高H*=H/d。可以看出越浪量隨相對(duì)波高增大而明顯增大,曲線呈現(xiàn)出線性變化趨勢(shì),該結(jié)果與Baldock等[13]、張金牛等[14]實(shí)驗(yàn)所測(cè)孤立波與防波堤作用所得越浪量的變化規(guī)律類似。對(duì)比數(shù)值結(jié)果與不同經(jīng)驗(yàn)公式,可以看出,3種結(jié)果雖然趨勢(shì)保持一致,但越浪量大小有著明顯區(qū)別,Baldock等[13]所提出的經(jīng)驗(yàn)公式是基于孤立波與緩坡作用模型,而張金牛等[14]則通過模型實(shí)驗(yàn)提出了適用于孤立波與陡坡模型作用的經(jīng)驗(yàn)公式,因此Baldock等[13]公式所得結(jié)果最大。相比于以上結(jié)果,本文的OWC模型更類似于以上適用條件的極限情況,因此所得結(jié)果最小。

圖5 相對(duì)波高對(duì)越浪量影響(H*=H/d,q*=q/q0)
圖6(a) 給出無(wú)量綱越浪厚度沿程衰減規(guī)律,其中,C0為OWC前墻上方越浪流厚度;Ci為OWC上頂板任意位置處的越浪流厚度;xi為上頂板任意位置距OWC前墻距離;以順流方向?yàn)檎膱D7中可以看出,除x*≈1.1處外,越浪流厚度均沿程變小,衰減系數(shù)隨相對(duì)波高增大而增大,且在上頂板前半部分衰減速度大于后半部分。孤立波與OWC作用時(shí),先沿前墻進(jìn)行爬坡,靠近前墻附近的水體僅有垂直方向流速,但隨著爬坡水體厚度增加,部分水體由于左側(cè)水體作用仍具有部分水平速度(見4.4節(jié)),因而,越浪水體沿一角度入射至OWC上頂板,在上頂板底部形成一密閉空腔。由于空腔無(wú)法從兩側(cè)排出,隨越浪流逐漸順流運(yùn)動(dòng),在x*≈1.1處破碎釋放,所以引起水體厚度局部增大。盡管存在OWC氣室的非線性作用,但本文所得結(jié)果趨勢(shì)與Oumeraci等[15]提出的規(guī)則波與防波堤作用下的越浪厚度衰減經(jīng)驗(yàn)公式趨勢(shì)保持一致(圖6(b)),但對(duì)比張金牛[14]針對(duì)于孤立波越浪提出的改進(jìn)經(jīng)驗(yàn)公式,在堤頂后方與本文模擬結(jié)果具有較大差別,可能是由于實(shí)驗(yàn)中防波堤后方水流較淺,難以準(zhǔn)確測(cè)量導(dǎo)致。

圖7 越浪流傳播情況
越浪現(xiàn)象發(fā)生時(shí),越浪流在結(jié)構(gòu)物頂部傳播的速度稱為越浪流速度(水舌速度),越浪流速度大小將間接影響岸邊結(jié)構(gòu)物所受到的越浪沖擊力。圖7展示越浪水舌位置的時(shí)間歷程,及越浪流速度隨相對(duì)波高改變的變化規(guī)律。在越浪流傳播后期,其主要受到底部切應(yīng)力影響,最終自由液面會(huì)產(chǎn)生不連續(xù)現(xiàn)象,故圖7(a)中水舌位置僅選取連續(xù)的水面過程。由圖7中可以看出,在整個(gè)時(shí)間歷程中,越浪流主要以勻速傳播。圖中在初始階段有一明顯的速度轉(zhuǎn)折,原因與上節(jié)類似(越浪流呈一角度入射至上頂板)。越浪流傳播速度通過勻速運(yùn)動(dòng)過程中位移曲線斜率獲得。圖7(b)給出該速度隨相對(duì)波高變化關(guān)系,發(fā)現(xiàn)隨相對(duì)波高增大越浪流傳播速度增大,且在相對(duì)波高為0.375時(shí)存在一微小轉(zhuǎn)折,此后速度增大程度稍有變緩的趨勢(shì)。
借鑒Oumeraci等[15]研究越浪時(shí)把越浪劃分成若干個(gè)階段的方法,本文將整個(gè)過程劃分為4個(gè)特征階段。云圖中,箭頭僅代表流動(dòng)方向,流速大小由云圖顏色表示。
1)波浪到達(dá)前墻并與前墻作用產(chǎn)生變形,水體沿前墻繼續(xù)爬升,由于前墻不透水,緊貼前墻一側(cè)水體僅有垂直方向流速,而外側(cè)水體仍帶有部分水平速度分量(圖8(a))。

圖8 越浪流場(chǎng)
2)波浪越過前墻,失去前墻作用的水體由于慣性繼續(xù)向斜上方運(yùn)動(dòng),最終受到重力作用產(chǎn)生回落,回落的水舌以一角度入射至上頂板,在與上頂板作用瞬間產(chǎn)生了較為明顯的砰擊現(xiàn)象,引起局部受力增大(圖8(b))。
3)在砰擊作用發(fā)生的瞬間,由于存在入射角度,小部分空氣被卷入并由水體包裹隨越浪流的演化順流運(yùn)動(dòng)(圖8(c) 方框線),此時(shí)隨著越浪厚度逐漸減小,該空腔逐步被釋放。當(dāng)空腔經(jīng)過采集位置時(shí),造成瞬時(shí)越浪量減小與越浪厚度的突變。
4)隨著越浪的繼續(xù)進(jìn)行,左側(cè)波浪開始后退,波浪不再向上頂板輸送水體,越浪流左側(cè)動(dòng)壓降低,由于上頂板仍有水體存在,阻礙了水體繼續(xù)向前流動(dòng),故在重力作用下,部分水體從左側(cè)回流,此時(shí)出現(xiàn)了頂板后方水體厚度大于前方的情況(圖8(d))。
圖9展示上頂板在不同相對(duì)波高作用下,各測(cè)點(diǎn)波壓力的極值沿水平位置的變化情況。由圖9中可以看出,波壓力極值隨相對(duì)波高增大而明顯減小,不同相對(duì)波高下,頂板波壓力極值存在一個(gè)突變點(diǎn),突變點(diǎn)處波壓力明顯增大并且隨相對(duì)波高增大,突變量變大,作用位置逐漸向后方移動(dòng),突變點(diǎn)位置對(duì)應(yīng)于砰擊現(xiàn)象所發(fā)生的位置,即砰擊點(diǎn)。根據(jù)Li[16]所述,越浪流在與前墻作用后,與上頂板作用前,垂直方向作自由落體運(yùn)動(dòng),水平方向做勻速運(yùn)動(dòng)。由于初始入射速度隨相對(duì)波高增加而增加,進(jìn)而導(dǎo)致回落時(shí)間增加,所以當(dāng)發(fā)生越浪現(xiàn)象時(shí),水舌在空中水平運(yùn)動(dòng)距離變大,最終導(dǎo)致砰擊點(diǎn)逐漸向后方移動(dòng)。
1)以往基于規(guī)則波或不規(guī)則波與防波堤作用模型所提出的越浪量經(jīng)驗(yàn)公式與如今的數(shù)值模擬結(jié)果趨勢(shì)相同,即越浪量均隨相對(duì)波高的增大而線性增大,但其結(jié)果明顯偏大,因此不適用于OWC模型。
2)越浪厚度與經(jīng)驗(yàn)公式對(duì)比良好,隨相對(duì)波高增大而增大,因此OWC氣室的存在對(duì)越浪厚度影響較小。越浪流在頂部傳播時(shí)以勻速階段為主,該速度隨相對(duì)波高增大而增加。
3)通過流場(chǎng)分析發(fā)現(xiàn),在越浪現(xiàn)象發(fā)生的同時(shí),由于越浪流沿一角度與上頂板作用,導(dǎo)致了頂板砰擊現(xiàn)象的發(fā)生。砰擊壓強(qiáng)隨相對(duì)波高增大而降低,砰擊點(diǎn)也隨之后移。因此需要在設(shè)計(jì)時(shí)考慮采取相應(yīng)措施減小該砰擊現(xiàn)象所帶來(lái)的危害。
本文以期為振蕩水柱波能轉(zhuǎn)換裝置的優(yōu)化及安全設(shè)計(jì)提供理論指導(dǎo)。