孤立波作用下振蕩水柱波能轉(zhuǎn)換裝置的越浪研究

牟迪, 王榮泉, 寧德志, 陳麗芬

(大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室,遼寧 大連 116024)

海嘯常指地球地殼運(yùn)動引起的海浪。通常引起海嘯的原因包括但不限于:地震、風(fēng)暴潮、水下爆炸等。當(dāng)以上事件發(fā)生時,水面的突然抬升或者降低將產(chǎn)生沿各個方向傳播的波浪,并產(chǎn)生海嘯。結(jié)構(gòu)物在海嘯波作用下遭受猛烈的極端載荷并引起流體諸多強(qiáng)非線性現(xiàn)象。全球在過去的100年中,每年大約發(fā)生一次毀滅性海嘯,在世界范圍內(nèi)造成嚴(yán)重的人員損失與財產(chǎn)破壞[1]。海嘯在產(chǎn)生時,波幅較小,通常僅有1～2 m,但與當(dāng)?shù)厮钕啾?波長很長,通常在100 km左右。孤立波只有一個較大的波峰或波谷,其波形和水質(zhì)點(diǎn)的運(yùn)動與海嘯極為相似,因此常被用來模擬海嘯等淺水大波[2]。

從大量海洋災(zāi)害所造成破壞的案例看來,海嘯引發(fā)的越浪現(xiàn)象,是造成海洋結(jié)構(gòu)物破壞及財產(chǎn)損失的重要原因之一。自20世紀(jì)起,國內(nèi)外眾多學(xué)者已經(jīng)對于越浪問題進(jìn)行了大量的研究并取得了大量研究成果[3-5],但以往的研究多集中于海嘯波與不同結(jié)構(gòu)形式的海堤或者近岸跨海橋梁等海岸結(jié)構(gòu)物作用時產(chǎn)生的越浪現(xiàn)象,而關(guān)于孤立波與振蕩水柱式(oscillation water column, OWC)波浪能轉(zhuǎn)換裝置作用下引起的越浪現(xiàn)象研究較少。

OWC波能裝置的主體帶有一個下端開口的氣室結(jié)構(gòu),氣室上部有通氣管道連接空氣透平。在波浪作用下,水柱上下振蕩迫使氣室內(nèi)部空氣在通氣管道處形成往復(fù)氣流,氣流通過推動空氣透平帶動發(fā)電機(jī)工作產(chǎn)生電能。其與普通海洋結(jié)構(gòu)物不同之處主要在于結(jié)構(gòu)內(nèi)部氣室的存在,氣室內(nèi)部上方氣流由于氣孔的阻礙作用使氣壓增大,進(jìn)而抑制氣室內(nèi)部自由液面的上升。被抑制的水柱通過壓強(qiáng)的傳遞進(jìn)而影響波浪在前墻上的爬坡運(yùn)動最終對越浪產(chǎn)生影響。以往對該裝置的研究多從水動力轉(zhuǎn)換效率角度出發(fā),優(yōu)化該裝置結(jié)構(gòu),提升其波能轉(zhuǎn)換能力[6-8]。而對于該裝置在極端波況下的水動力特性,國內(nèi)外大多學(xué)者關(guān)注點(diǎn)集中在裝置的受力情況,如馬子然[9]研究了在聚焦波作用下,離岸式OWC裝置的受力情況。盡管如此,目前未見有關(guān)于極端波浪作用下,OWC波能裝置越浪特性的研究。

所以本文基于以上背景,利用開源軟件OpenFOAM 中waveIsoFoam求解器求解三維Navier-Stokes方程,對孤立波與OWC裝置作用時發(fā)生的越浪現(xiàn)象進(jìn)行數(shù)值模擬研究,分析與越浪相關(guān)的物理量如越浪量、越浪厚度、越浪流速度等隨相對波高改變的變化規(guī)律,并展示越浪發(fā)生時的流場特性。為OWC波能裝置的優(yōu)化及安全設(shè)計提供理論指導(dǎo)。

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

1.1 控制方程

本文假定流體不可壓,無熱傳導(dǎo)?；谏鲜黾俣?N-S方程可分為連續(xù)性方程及動量方程:

(1)

(2)

式中ρ、μ表示混合流體密度及流體粘度。u=(u,v,w)和x=(x,y,z)為笛卡爾坐標(biāo)系下流體速度矢量與位置坐標(biāo),p*為動壓,g為重力加速度。

1.2 自由液面捕捉

本文通過VOF方法獲得自由液面運(yùn)動。該方法求解運(yùn)輸方程(3)以獲得體積分?jǐn)?shù):

(3)

式中:α為每個網(wǎng)格體心的體積分?jǐn)?shù),如果α的值界于0～1,則這個網(wǎng)格則被標(biāo)記為包含自由液面的計算網(wǎng)格。當(dāng)α=0時,對應(yīng)為氣相;當(dāng)α=1時,對應(yīng)水相。單元內(nèi)流體物理屬性為:

(4)

式中:角標(biāo)a與w分別代表空氣與水,盡管基于不可壓縮假定,但是交界面處密度隨計算時間步推進(jìn)而時刻變化。

在本文中,用于重構(gòu)自由液面的VOF方法為Iso-advector 方法,該方法是一種基于幾何特征的重構(gòu)方法。相比于基于代數(shù)重構(gòu)的MULES方法,該方法在計算效率、收斂性、穩(wěn)定性有著較強(qiáng)的優(yōu)勢,關(guān)于該方法詳細(xì)描述見文獻(xiàn)[10-11]。

1.3 孤立波方程

本文采用waves2Foam中提供的一階孤立波波面方程進(jìn)行數(shù)值造波,并在計算域兩側(cè)分別布置被動消波區(qū)域。其中孤立波波面方程為:

(5)

(6)

(7)

(8)

式中:η為孤立波波面;H為孤立波波高;d為水深;c為孤立波波速。

1.4 離散格式及邊界條件

本文采用有限體積法將計算域離散成眾多結(jié)構(gòu)化網(wǎng)格。對于方程(2)中時間項(xiàng)采用一階Euler隱式離散,對流項(xiàng)采用linearUpwind 二階TVD格式離散,在保證計算精度的同時保證其求解的穩(wěn)定性,速度梯度通過Gauss linear interpolation 方法進(jìn)行插值計算。速度壓力耦合通過pimple算法求解,pimple算法為simple算法與piso算法的結(jié)合,允許用戶根據(jù)所設(shè)定的最大庫郎數(shù)進(jìn)行時間步的自適應(yīng)調(diào)節(jié)[12],本文計算中庫郎數(shù)為0.3。

本文中所采用的邊界條件如下:對于所有結(jié)構(gòu)物物面及底部地面,采用noslip無滑移邊界條件并配合zeroGradient 零梯度壓力邊界條件;對于計算域頂部采用totalPressure及pressureInletOutlet邊界,即在頂部指定總壓,并僅允許空氣流出,不允許回流。

2 模型驗(yàn)證

2.1 網(wǎng)格收斂性驗(yàn)證

為排除空間離散對數(shù)值結(jié)果產(chǎn)生的影響,本節(jié)將進(jìn)行網(wǎng)格收斂性驗(yàn)證。表1列出了不同網(wǎng)格參數(shù),結(jié)構(gòu)物外部水平方向網(wǎng)格與垂向網(wǎng)格尺寸相等。由于氣孔處通量變化較大,對其附近區(qū)域進(jìn)行了局部加密(見圖1)。數(shù)值模型中計算域尺寸與物模實(shí)驗(yàn)中保持一致(見下節(jié))。圖2展示了不同網(wǎng)格下由OpenFOAM模擬得到的x=9.3 m處波面歷程,本文模擬中水深保持d=0.4 m。圖2中,不同網(wǎng)格得到的結(jié)果差距較小,為平衡數(shù)值結(jié)果與計算量,最終選擇中等程度的網(wǎng)格進(jìn)行計算。

表1 網(wǎng)格尺度

圖1 網(wǎng)格示意

圖2 H*=H/d=0.425時不同網(wǎng)格下波面歷程(η*=η/d,t*=t(g/d)1/2)

2.2 模型驗(yàn)證

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

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

圖3 本次實(shí)驗(yàn)示意

圖4分別展示數(shù)值模擬結(jié)果與模型實(shí)驗(yàn)結(jié)果的對比。

圖4(a)第1個峰值為入射波,第2個峰值為反射波。OpenFOAM計算得到的反射波峰值略大于實(shí)驗(yàn)結(jié)果,其原因?yàn)閷?shí)驗(yàn)水槽側(cè)壁高度受限,波浪與前墻作用時部分水體由水槽兩側(cè)溢出而導(dǎo)致能量部分損失,最終使反射波波高有所降低。圖4(b)展示氣室內(nèi)氣壓對比,其中正壓對應(yīng)于OWC運(yùn)動的上升階段,當(dāng)氣室內(nèi)液面高度達(dá)到最大時,氣室內(nèi)氣壓降至零。此后,OWC開始做下降運(yùn)動,氣室內(nèi)部產(chǎn)生負(fù)壓,外部空氣被吸入至氣室內(nèi)部。圖4(c)展示結(jié)構(gòu)總水平力對比,受力曲線呈現(xiàn)單峰狀,受力峰值與實(shí)驗(yàn)結(jié)果誤差小于2%?？傮w來說,數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果對比良好,驗(yàn)證了OpenFOAM模型的有效性。

3 結(jié)果與討論

3.1 越浪量變化規(guī)律

采用Baldock等[13]提出的方法,選取孤立波在靜水面以上的單寬水體體積作為特征量定義無量綱越浪量即:

q*=q/q0

式中:q為數(shù)值模擬得到的單寬越浪量;q0為孤立波在靜水面以上的單寬水體體積,對方程(5)積分可得:

(9)

圖5給出無量綱越浪量與相對波高的關(guān)系,所監(jiān)測斷面為前墻上方,其中定義相對波高H*=H/d。可以看出越浪量隨相對波高增大而明顯增大,曲線呈現(xiàn)出線性變化趨勢,該結(jié)果與Baldock等[13]、張金牛等[14]實(shí)驗(yàn)所測孤立波與防波堤作用所得越浪量的變化規(guī)律類似。對比數(shù)值結(jié)果與不同經(jīng)驗(yàn)公式,可以看出,3種結(jié)果雖然趨勢保持一致,但越浪量大小有著明顯區(qū)別,Baldock等[13]所提出的經(jīng)驗(yàn)公式是基于孤立波與緩坡作用模型,而張金牛等[14]則通過模型實(shí)驗(yàn)提出了適用于孤立波與陡坡模型作用的經(jīng)驗(yàn)公式,因此Baldock等[13]公式所得結(jié)果最大。相比于以上結(jié)果,本文的OWC模型更類似于以上適用條件的極限情況,因此所得結(jié)果最小。

圖5 相對波高對越浪量影響(H*=H/d,q*=q/q0)

3.2 越浪厚度沿程變化規(guī)律

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

圖7 越浪流傳播情況

3.3 越浪流傳播速度

越浪現(xiàn)象發(fā)生時,越浪流在結(jié)構(gòu)物頂部傳播的速度稱為越浪流速度(水舌速度),越浪流速度大小將間接影響岸邊結(jié)構(gòu)物所受到的越浪沖擊力。圖7展示越浪水舌位置的時間歷程,及越浪流速度隨相對波高改變的變化規(guī)律。在越浪流傳播后期,其主要受到底部切應(yīng)力影響,最終自由液面會產(chǎn)生不連續(xù)現(xiàn)象,故圖7(a)中水舌位置僅選取連續(xù)的水面過程。由圖7中可以看出,在整個時間歷程中,越浪流主要以勻速傳播。圖中在初始階段有一明顯的速度轉(zhuǎn)折,原因與上節(jié)類似(越浪流呈一角度入射至上頂板)。越浪流傳播速度通過勻速運(yùn)動過程中位移曲線斜率獲得。圖7(b)給出該速度隨相對波高變化關(guān)系,發(fā)現(xiàn)隨相對波高增大越浪流傳播速度增大,且在相對波高為0.375時存在一微小轉(zhuǎn)折,此后速度增大程度稍有變緩的趨勢。

3.4 越浪水體形態(tài)及砰擊現(xiàn)象

借鑒Oumeraci等[15]研究越浪時把越浪劃分成若干個階段的方法,本文將整個過程劃分為4個特征階段。云圖中,箭頭僅代表流動方向,流速大小由云圖顏色表示。

1)波浪到達(dá)前墻并與前墻作用產(chǎn)生變形,水體沿前墻繼續(xù)爬升,由于前墻不透水,緊貼前墻一側(cè)水體僅有垂直方向流速,而外側(cè)水體仍帶有部分水平速度分量(圖8(a))。

圖8 越浪流場

2)波浪越過前墻,失去前墻作用的水體由于慣性繼續(xù)向斜上方運(yùn)動,最終受到重力作用產(chǎn)生回落,回落的水舌以一角度入射至上頂板,在與上頂板作用瞬間產(chǎn)生了較為明顯的砰擊現(xiàn)象,引起局部受力增大(圖8(b))。

3)在砰擊作用發(fā)生的瞬間,由于存在入射角度,小部分空氣被卷入并由水體包裹隨越浪流的演化順流運(yùn)動(圖8(c) 方框線),此時隨著越浪厚度逐漸減小,該空腔逐步被釋放。當(dāng)空腔經(jīng)過采集位置時,造成瞬時越浪量減小與越浪厚度的突變。

4)隨著越浪的繼續(xù)進(jìn)行,左側(cè)波浪開始后退,波浪不再向上頂板輸送水體,越浪流左側(cè)動壓降低,由于上頂板仍有水體存在,阻礙了水體繼續(xù)向前流動,故在重力作用下,部分水體從左側(cè)回流,此時出現(xiàn)了頂板后方水體厚度大于前方的情況(圖8(d))。

圖9展示上頂板在不同相對波高作用下,各測點(diǎn)波壓力的極值沿水平位置的變化情況。由圖9中可以看出,波壓力極值隨相對波高增大而明顯減小,不同相對波高下,頂板波壓力極值存在一個突變點(diǎn),突變點(diǎn)處波壓力明顯增大并且隨相對波高增大,突變量變大,作用位置逐漸向后方移動,突變點(diǎn)位置對應(yīng)于砰擊現(xiàn)象所發(fā)生的位置,即砰擊點(diǎn)。根據(jù)Li[16]所述,越浪流在與前墻作用后,與上頂板作用前,垂直方向作自由落體運(yùn)動,水平方向做勻速運(yùn)動。由于初始入射速度隨相對波高增加而增加,進(jìn)而導(dǎo)致回落時間增加,所以當(dāng)發(fā)生越浪現(xiàn)象時,水舌在空中水平運(yùn)動距離變大,最終導(dǎo)致砰擊點(diǎn)逐漸向后方移動。

4 結(jié)論

1)以往基于規(guī)則波或不規(guī)則波與防波堤作用模型所提出的越浪量經(jīng)驗(yàn)公式與如今的數(shù)值模擬結(jié)果趨勢相同,即越浪量均隨相對波高的增大而線性增大,但其結(jié)果明顯偏大,因此不適用于OWC模型。

2)越浪厚度與經(jīng)驗(yàn)公式對比良好,隨相對波高增大而增大,因此OWC氣室的存在對越浪厚度影響較小。越浪流在頂部傳播時以勻速階段為主,該速度隨相對波高增大而增加。

3)通過流場分析發(fā)現(xiàn),在越浪現(xiàn)象發(fā)生的同時,由于越浪流沿一角度與上頂板作用,導(dǎo)致了頂板砰擊現(xiàn)象的發(fā)生。砰擊壓強(qiáng)隨相對波高增大而降低,砰擊點(diǎn)也隨之后移。因此需要在設(shè)計時考慮采取相應(yīng)措施減小該砰擊現(xiàn)象所帶來的危害。

本文以期為振蕩水柱波能轉(zhuǎn)換裝置的優(yōu)化及安全設(shè)計提供理論指導(dǎo)。