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

(大連理工大學 海岸和近海工程國家重點實驗室,遼寧 大連 116024)

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

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

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

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

1 數(shù)學模型

1.1 控制方程

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

(1)

(2)

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

1.2 自由液面捕捉

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

(3)

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

(4)

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

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

1.3 孤立波方程

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

(5)

(6)

(7)

(8)

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

1.4 離散格式及邊界條件

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

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

2 模型驗證

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

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

表1 網(wǎng)格尺度

圖1 網(wǎng)格示意

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

2.2 模型驗證

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

表2 結構物尺寸

圖3 本次實驗示意

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

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

3 結果與討論

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

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

q*=q/q0

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

(9)

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

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

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

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

圖7 越浪流傳播情況

3.3 越浪流傳播速度

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

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

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

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

圖8 越浪流場

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

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

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

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

4 結論

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

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

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

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