波浪對雙層圓弧型貫底式開孔介質(zhì)防波堤的繞射

林皋，劉俊

(1.大連理工大學(xué) 海岸和近海工程國家重點實驗室，遼寧 大連 116024;2.大連理工大學(xué) 建設(shè)工程學(xué)部 水利工程學(xué)院，遼寧 大連116024)

20世紀60年代由加拿大學(xué)者Jarlan率先提出了開孔板式防波堤［1］，該結(jié)構(gòu)可以有效地降低結(jié)構(gòu)波浪力和反射系數(shù)的值.隨后許多研究者提出并研究了許多改進型的開孔防波堤結(jié)構(gòu)，以滿足實際工程的不同需要［2-3］.與此同時，為了在更廣泛的頻率范圍內(nèi)降低結(jié)構(gòu)波浪力和反射系數(shù)的值，兩層或多層開孔板式防波堤是較佳的選擇.如意大利波爾圖Torres海港中修建了三層開孔板沉箱結(jié)構(gòu)［4］，最近，中國大連市化工港正在修建五層開孔板沉箱結(jié)構(gòu)［3］.在理論研究方面，Sawaragi等［5］研究了波浪與雙層開孔板沉箱作用條件下的衰減，結(jié)果表明，在特定的條件下，雙層開孔板沉箱對波浪能量的吸收比單層開孔板要好.Kondo［6］發(fā)展了一種線性長波條件下的解析模型.李玉成等［7］應(yīng)用特征函數(shù)展開法分析了雙層局部開孔板沉箱對波浪反射.汪宏等［8］通過實驗給出了雙層開孔直立式消浪板不同開孔形式在規(guī)則波作用下的透射系數(shù).劉勇等［9］利用特征函數(shù)展開方法研究了無窮層局部開孔板沉箱對波浪反射.有關(guān)更多的波浪與雙層或者多層開孔結(jié)構(gòu)的相互作用問題研究可詳見文獻［3］.但是，以上文獻大多局限于二維平面波問題.程建生等［10］用特征函數(shù)展開解析方式研究過平面波作用下開孔介質(zhì)圓弧貫底式防波堤的防浪效果的三維問題.然而，海浪實際上非常復(fù)雜，往往需要采用短峰波來描述［11］.迄今為止未見短峰波對雙層圓弧型貫底式開孔介質(zhì)防波堤的繞射問題的報道.

本文將采用一種新的半解析方法－比例邊界有限元方法(SBFEM)建立短峰波對雙層圓弧型貫底式開孔介質(zhì)防波堤的繞射問題的分析模型.在此基礎(chǔ)上分析了短峰波波浪參數(shù)和結(jié)構(gòu)的幾何參數(shù)等對整個結(jié)構(gòu)繞射的影響.

1 雙層弧型開孔介質(zhì)防波堤物理模型

如圖1所示，設(shè)在均勻水深h的海域中設(shè)置有內(nèi)、外環(huán)半徑分別為a和b的雙層圓弧形開孔介質(zhì)防波堤.取坐標系O-xyz，其中，O-xy平面位于靜水面上，原點在圓筒柱的圓心處，O-z軸垂直向上.將2個圓弧延伸構(gòu)建成2個虛擬圓柱面，整個流場被劃分成1個無限子域Ω3和2有限子域Ω1、Ω2.根據(jù)線性理論，無限子域Ω3中的速度總勢Φ3可表示為入射勢和散射勢ΦSA3之和，即Φ3=+;兩個有限域的速度總勢分別用Φ1和Φ2表達.為方便下面公式推導(dǎo)，用速度勢Φ統(tǒng)一代表Φ1、Φ2、Φ3、以及.本問題為三維問題，應(yīng)用分離變量法，可分離出時間項t和空間變量z，即速度勢Φ(x，y，z，t)可表達為

式中:Z(z)代表z方向的函數(shù)，ω為角頻率，t為時間項.滿足邊界條件的解為(k為波數(shù)):

圖1 雙層圓弧型開孔介質(zhì)防波堤及比例邊界有限元模型Fig.1 Double-layered arc-shaped porous breakwater and model of SBFEM

與Φ類似，無限子域Ω3中，速度總勢φ3表示為入射勢和散射勢之和，即φ3=+有限子域可直接應(yīng)用速度總勢φ1、φ2.其中，、φ1和φ2均可由二維Helmholtz方程進行描述可統(tǒng)一地表示為φ(x，y)，即:

兩有限域界面處的邊界條件可由下式表達:

其中，對于多孔介質(zhì)圓弧段標量G1等于其自身的孔隙影響系數(shù)，對非多孔介質(zhì)圓弧段G值取為無窮大.

有限域Ω2和無限域Ω3界面處的邊界條件可由下式表達:

類似地，G2為有限域Ω和無限域Ω3界面處的孔隙影響系數(shù).

入射速度勢可由下式表達:

總場速度勢φ求出之后，其他變量如波速度，波面高度和波動壓力等可通過以下公式求出:

2 比例邊界有限元方程推導(dǎo)

通過以上闡述可以發(fā)現(xiàn)，問題歸結(jié)為求解有限域的速度總勢φ1、φ2和無限域的散射勢，使之滿足方程和相應(yīng)的邊界條件.仍用φ統(tǒng)一表示φ1、φ2和，并設(shè)φ在各子域邊界上的法向?qū)?shù)為

方程(2)和(8)的等效泛函變分可表達為

圖2 比例邊界有限元方法計算Fig.2 The coordinate definition of SBFEM

應(yīng)用SBFEM求解泛函變分時，首先需要建立笛卡爾坐標系(x，y)和比例邊界有限元坐標系(ξ，s)之間的轉(zhuǎn)換關(guān)系.在域內(nèi)選擇一點О作為SBFEM的相似中心(通過此中心，該域的全部邊界可見，也就是說只有半徑為b圓邊界需要離散)，分別建立一個圍繞比例中心的環(huán)向坐標s(s在數(shù)量上可表示為從E點開始逆時針旋轉(zhuǎn)半徑為b的圓弧長度)和從比例中心出發(fā)向計算區(qū)域邊界輻射的輻射坐標ξ，如圖2所示.環(huán)向坐標的定義從側(cè)邊面x=sb開始，逆時針旋轉(zhuǎn)到側(cè)邊面s=sa，當sa和sb重合時，計算域為封閉的圓或圓環(huán)，本文中sa和sb是重合的.輻射坐標在比例中心處為ξ=0，在計算區(qū)域的半徑為b邊界上ξ=1.結(jié)合圖1和圖2，對于本文中計算情況:有限域Ω1中0=ξ0≤ξ=ξ1=a/b;有限域Ω2中a/b=ξ0≤ξ=ξ1≤1;無限域Ω2中1=ξ0≤ξ=ξ1=∞.在笛卡爾坐標系中計算域內(nèi)的一點可表示為(^x，^y)，在邊界處(ξ=1)上則表示為(x，y).相似中心處表示坐標為(x0，y0).圖2中，陰影部分代表一SBFEM單元.兩坐標系統(tǒng)的變換關(guān)系如下(s可以通過邊界節(jié)點插值得出，本文采用三節(jié)點二次曲線插值函數(shù)N(s)，即s=N(s)s，s，為三節(jié)點弧長):

兩坐標系統(tǒng)的求導(dǎo)，通過雅克比矩陣加以聯(lián)系:

其中，

在比例邊界坐標下，梯度算子可以表示為

其中:

根據(jù)等參變換概念，速度勢函數(shù)可采用同樣的插值函數(shù)進行離散:

將式(13)、(14)代入式(9)，經(jīng)整理后可得:

其中:

引入以下系數(shù)矩陣E0、E1、E2、M0和Fs(ξ):

式中，F(xiàn)s(ξ)為側(cè)面邊界(s=sa、s=sb)第二類邊界條件不為零時在s=sa、s=sb處產(chǎn)生的通量.當= 0時或s=sb、s=sb存在第一類邊界條件時，F(xiàn)s(ξ)= 0，注意到本文中Fs(ξ)=0.

對式中含δφ(ξ)，ξ的項做分部積分，考慮到δφ (ξ)任意性，并令ζ=kbξ，最終可得如下簡化方程［13］:

方程為各個子域的 SBFEM的基本方程，式(17)、(18)為其內(nèi)、外邊界條件.

3 比例邊界有限元方程求解

3.1 對無限域Ω3的求解

方程為貝塞爾(Bessel)微分方程組，考慮到無窮遠的 Sommerfeld輻射邊界條件，其解可以Hri(ζ)Ti作為基函數(shù)來表示:

其中，rj為漢克爾函數(shù)的階數(shù)，Tj為n階向量，ci為參與系數(shù)，Hrj(ζ)為第一類漢克爾(Hankel)函數(shù).

將式(19)代入方程(16)，并且考慮到漢克爾函數(shù)的導(dǎo)數(shù)性質(zhì)，經(jīng)整理后可得［13］:

式(21)為一個特征值問題.設(shè) λj為矩陣的特征值，則，并且T為的特征向量.

由于漢克爾函數(shù)在無窮遠處自動滿足Sommerfeld輻射邊界條件，所以只需考慮比例邊界有限元在ξ=1處的邊界條件:

3.2 對有限域Ω1、Ω2的求解及其和無限域Ω2耦合

類似地，有限域Ω1中的總場φ1可用Bessel函數(shù)作為解的基函數(shù):

其中，

同時考慮SBFEM的邊界條件可得:

其中:

同理有限域Ω2中的總場φ1可表示為:

考慮內(nèi)邊界條件(17)、外邊界條件(18)分別可得Ω2中的SBFEM內(nèi)外邊界條件:

其中:

考慮到有限域Ω1和Ω2之間的邊界條件可得:

其中:

其中，G1為多孔介質(zhì)的孔隙影響系數(shù)對角矩陣，對角元素為SBFEM單元對應(yīng)的孔隙影響系數(shù)，對于多孔介質(zhì)圓弧段選取自身的孔隙影響系數(shù)，而對非多孔介質(zhì)圓弧段選為無窮大.

并注意到有限域Ω2和無限域Ω3界面存在如下關(guān)系:

再考慮方程限域Ω2和無限域Ω3界面的協(xié)調(diào)邊界條件(4)可解出C3:

其中，對角矩陣G2定義類似T1，ˉvI3n為Ω3中邊界的入射波節(jié)點速度，而且:

將式(32)代入式(28)可得C2，進一步將C2C3代入方程(31)、(29)可分別求出C和C1.將系數(shù)向量C、C1、C2和C3代入各自的域的基本解，可得到各自域的速度總勢.

4 計算結(jié)果及分析

為了檢驗SBFEM的準確性和高效性.選取外接開孔結(jié)構(gòu)物為閉合的圓形作為算例.該問題對入射波為平面波(ky=0)時有解析解.選取計算條件為: h=15 m，A=0.5 m，b=10 m，a/b=0.2，G=1.0，kb≈0.8.對應(yīng)于本模型的參數(shù)為:雙層圓弧變成為全圓、G1=0(相當不透水壁)和G2=1.圖3給出了SBFEM和解析解［13］沿內(nèi)外柱壁歸一化波浪爬升值(|η|/A)的對比，結(jié)果表明SBFEM只需在外壁離散4個3節(jié)點單元就能很好的和解析解吻合.本文采用Intel Core Q8200(2.33GHz)計算機對算例進行了計算時間測試.即使在離散16個SBFEM單元的情況下，其消耗的時間不到1s，由此可以看出計算效率非常高.

圖4進一步給出了以上5種不同角度入射波(取值分別為β=0、π/8、π/4、3π/8、π/2.β定義如下:β=arctan(ky/kx).其中，β=0對應(yīng)為平面波; β=π/2則為駐波)作用下，掩護區(qū)中心點處(r=0)的歸一化的波浪爬升(|η/a|)的變化規(guī)律.模擬計算中的參數(shù)為:兩圓弧段的中心點在D，圓弧段的角度θ=120°(見圖1)，h=15 m，A=0.5 m，b=2 m，a/b=0.5，G1=0.1，G2=1.0.從圖4可以看出，波浪爬升呈現(xiàn)震蕩趨勢.隨著β不斷增大，掩護區(qū)域內(nèi)中心處的波浪爬升增大.

圖5進一步給出短峰波(kx=ky=sprt(2)/2)作用下，3種圓弧型開孔型防波堤結(jié)構(gòu)(θ=120°):1)存在雙層圓弧開孔結(jié)構(gòu);2)只存在單獨內(nèi)層圓弧開孔結(jié)構(gòu);3)只存在單獨內(nèi)層圓弧開孔結(jié)構(gòu)，中心點處(r=0)的歸一化的波浪爬升(|η/A|)對比結(jié)果.從圖中可以明顯發(fā)現(xiàn)，雙層開孔結(jié)構(gòu)比單層開孔結(jié)構(gòu)的在更廣泛的頻率范圍內(nèi)具有更好的防護效果，同時可以發(fā)現(xiàn)，在中間頻率段(大致ka=1.6～3.7)，單獨內(nèi)層和單獨外層情況的中心點處的波浪爬高非常接近.

圖3 平面波作用下的內(nèi)、外柱壁波浪爬升Fig.3 Run-up of an incident plane wave around the interior and exterior cylinders.

圖6分別表明了內(nèi)、外圓弧段孔隙影響系數(shù)G1(G2=1.0)、G2(G1=1.0)變化對中心點處(R=0)的歸一化的波浪爬升(|η/A|)的影響.除G1和G2變化外，其他計算參數(shù)和圖5一致.從圖6(a)可以看出，對隨著孔隙影響系數(shù)的逐漸增大，中心點處的波浪爬升先是快速下降到最小值后，再逐步緩慢增大并最終趨于穩(wěn)定值.而且隨著β(除駐波外)的增大，波浪爬升最小值點往更小G2的靠近.從圖6(b)可以看出，中心點處的波浪爬升的變化規(guī)律可以分成2種情況:1)對于較小β值(β=0，π/8)的波，隨著G2的增大，中心點處的波浪爬升也是先減小再增大并趨于一漸進值;2)對于較β大值(β≥π/4)的波，隨著G2的增大，中心點處的波浪爬升快速增大并趨于一漸近值.

圖4 隨相對波數(shù)變化下中心點處(r=0)波浪爬升Fig.4 Variation of wave run-up at r=0 versus

圖5 存在不同結(jié)構(gòu)層下的中心點處(r=0)波浪爬升對比Fig.5 Comparison of wave run-up at r=0 between different structures

圖6 隨孔隙影響系數(shù)G1、G2變化下中心點處波浪爬升Fig.6 Variation of wave run-up versus G1，G2

在以上兩種分析結(jié)果基礎(chǔ)上，圖7給出了不同波對域內(nèi)繞射波輪廓的影響.其中，分別選取3種波，即:平面波(kxa=1=1.0，kya=0)，短峰波(kxa=0，kya=1.0)，并且選擇了3種不同的孔隙影響系數(shù)G1=G2分別為0.1、1.0和10.0.其他計算參數(shù)和圖4一致.從圖中可以看出，所有的繞射波輪廓圖都關(guān)于x軸對稱;隨著G1和G2的減小，外圓弧段和內(nèi)筒間的波浪高度呈現(xiàn)明顯減小趨勢.垂直入射平面波能體現(xiàn)出明顯的反射特性，即掩護效果更好，但當孔隙影響系數(shù)G1和G2很大情況下，也表現(xiàn)得不是很明顯;而駐波和短峰波的繞射波相對平面波的更為復(fù)雜.

圖8表明了內(nèi)外圓弧半徑比a/b變化對中心點處(r=0)的歸一化的波浪爬升(|η/A|)的影響.其他計算參數(shù)和圖4一致.從圖中可以看出，隨著歸一化相對波數(shù)kb/2的增大，中心點處的波浪爬升先是快速下降，降到到一定值后又呈現(xiàn)震蕩變化.在低頻區(qū)，半徑比a/b對中心點處的波浪爬升影響非常小，而在中高頻率段，半徑比a/b只影響震蕩的波峰及波谷出現(xiàn)的位置.

圖9表明了3種組合圓弧張開角度θ變化對中心點處(r=0)歸一化的波浪爬升(|η|/A)的影響.除圓弧張開角度θ外，其他計算參數(shù)和圖4一致(兩圓弧段的中心點在D).3種組合分別為:1)內(nèi)圓弧段的角度變化而外圓弧段角度固定為120°;2)外圓弧段的角度變化而內(nèi)圓弧段角度固定為120°;3)內(nèi)、外圓弧段的角度同時變化.從圖中可以看出，總體來說，不管那種組合情況，隨著張開角度的增大，中心點處的波浪爬升震蕩的振幅增大;組合情況2)的震蕩的頻率更大;而且可以進一步證實，組合情況3)在更大的頻率范圍內(nèi)掩護效果更好.

圖10表明了3種組合圓弧位置變化γ對對中心點處(r=0)的歸一化的波浪爬升(|η|/A)的影響.圓弧中心點從A點逆時針轉(zhuǎn)到D點(總共180°).其他計算參數(shù)和圖6一致.3種組合分別為:1)內(nèi)圓弧段的位置變化而外圓弧段位置不變;2)外圓弧段的位置變化而內(nèi)圓弧段位置不變;3)內(nèi)、外圓弧段的位置同時變化.從圖中可以看出，不管哪種組合情況，隨著r的增大，中心點處的波浪爬升震蕩的振幅減小，即r值越大掩護效果更佳;組合情況2)的震蕩的頻率更大，同時組合情況3)的振幅更大.

圖7 不同短峰波在不同孔隙影響系數(shù)繞射波的波幅分布Fig.7 Diffracted wave amplitude distributions with different incident wave and different porous-effect parameters

圖8 隨內(nèi)、外半徑比變化下中心點處波浪爬升Fig.8 Variation of wave run-up versus a/b

圖9 隨圓弧段張開角度變化下中心點處波浪爬升Fig.9 Variation of wave run-up at r=0 versus opening angle

圖10 隨圓弧段位置變化下中心點處波浪爬升Fig.10 Variation of wave run-up at r=0 versus location of arcs

5 結(jié)論

開展了應(yīng)用比例邊界有限元法探討短峰波對雙層圓弧型貫底式開孔介質(zhì)防波堤的繞射問題.將開孔圓弧段延伸形成虛擬圓弧段，以不同的孔隙影響系數(shù)反映實際結(jié)構(gòu)段與虛擬圓弧段的不同特性，從而構(gòu)建出新的匹配邊界條件，能使該問題落入SBFEM求解范圍之內(nèi).基于比例邊界坐標變換和變分原理推導(dǎo)出SBFEM表示的本問題的基本方程.在求解過程中，針對有限域和無限域分別引入了系列貝塞爾函數(shù)和漢克爾函數(shù)作為其解的基函數(shù)，并且結(jié)合兩個有限域界面處以及有限域和無限域界面處的邊界條件求解出對應(yīng)的各系數(shù)向量.數(shù)值算例與解析解結(jié)果對比表明該方法具有很高的精度和計算效率.最后，詳細探討出諸如相對波數(shù)、半徑比、圓弧段的位置、張角以及孔隙影響系數(shù)G對整個結(jié)構(gòu)繞射的影響.結(jié)果表明:1)平面波的防浪效果最好，而駐波的效果最差;2)在更廣泛的頻率范圍內(nèi)，雙層開孔結(jié)構(gòu)對波浪的防浪效果比單層開孔結(jié)構(gòu)要好;3)在內(nèi)層的孔隙影響系數(shù)較小的情況下，內(nèi)層的孔隙影響系數(shù)增大將能有效地降低中心點處的波浪爬升;4)而外層的孔隙影響系數(shù)增大將會使中心點處的波浪爬升提高;5)在較低的頻率范圍內(nèi)，內(nèi)、外柱半徑比對波浪的防浪效果影響較小;6)內(nèi)、外柱的張角變化以及它們的位置變化對中心點處波浪爬升影響非常大.本研究可為工程上設(shè)計雙層圓弧型貫底式開孔介質(zhì)防波堤提供理論參考.

［1］JARLAN G E.A perforated vertical wall breakwater［J］.The Dock and Harbor Authority，1961，41(486):394-398.

［2］CHWANG A T，CHAN A T.Interaction between porous media and wave motion［J］.Annual Review of Fluid Mechanics，1998，30:53-84.

［3］HUANG Z H，LI Y C，LIU Y.Hydraulic performance and wave loadings of perforated/slotted coastal structures:a review［J］.Ocean Engineering，2011，38(10):1031-1053.

［4］FRANCE L.Vertical breakwaters:the Italian experience［J］.Coast Engineering，1994，22:31-55

［5］SAWARAGI T，IWATA K.Wave attenuation of a vertical breakwater with two air chambers［J］.Coastal Engineering in Japan，1978，21:63-74.

［6］KONDO H.Analysis of breakwaters having two porous walls［C］//Proceedings of Coastal Structures’79.Virginia，USA，1979.

［7］李玉成，劉洪杰，滕斌.雙層局部開孔板沉箱對波浪反射的理論研究［J］.海洋工程，2005，23(1):18-27.

LI Yucheng，LIU Hongjie，TENG Bin.Theoretical analysis of the wave reflection by caissons with double-layered perforated wall［J］.The Ocean Engineering，2005，23(1):18-27.

［8］汪宏，沈麗玉，王勇.雙層開孔直立式板結(jié)構(gòu)的消波性能試驗［J］.水運工程，2011，450(2):21-25.

WANG Hong，SHEN Liyu，WANG Yong.Test on wavedissipating performance of double-layer-hole vertical plate structure［J］.Port＆Waterway Engineering，2011，450 (2):21-25.

［9］LIU Y，LI Y C，TENG B.Interaction between obliquely incident waves and an infinite array of multi-chamber perforated caissons［J］.Journal of Engineering Mathematics，DOI 10.1007/s10665-011-9484-2.

［10］程建生，繆國平，王景全，等.圓弧型貫底式多孔介質(zhì)防波堤防浪效果的解析研究［J］.哈爾濱工程大學(xué)學(xué)報，2008，29(1):5-10.

CHENG Jiansheng，MIAO Guoping，WANG Jingquan，et al Wave-prevention efficiency of an arc-shaped bottom-mounted porous breakwater［J］.Journal of Harbin Engineering University，2008，29(1):5-10.

［11］TSAI C P，JENG D S，HSU J R C.Computations of the almost highest short-crested waves in deep water［J］.Applied Ocean Research，1994，16(6):317-326.

［12］WOLF J P，SONG C M.Consistent infinitesimal finite element cell method:three-dimensional vector wave equation［J］.International Journal for Numerical Methods in Engineering，1996，39(13):2189-2208.

［13］TAO L B，SONG H，CHAKRABARTI S.Scaled boundary FEM model for interaction of short-crested waves with a concentric porous cylindrical structure［J］.Journal of Waterway，Port，Coastal and Ocean Engineering，2009，135 (5):200-212.