波流混合作用的完全非線性數(shù)值水槽模型

寧德志,陳麗芬,田宏光

(1.大連理工大學(xué)海岸和近海工程國家重點實驗室,遼寧大連 116024;2.南京水利科學(xué)研究院水資源與水利工程科學(xué)國家重點實驗室,江蘇南京 210029)

隨著計算機技術(shù)和現(xiàn)代數(shù)值方法的發(fā)展,完全非線性數(shù)值水槽技術(shù)作為一種非常有效的工具已經(jīng)解決了許多復(fù)雜的水波問題[1-5].從以往的研究成果可以發(fā)現(xiàn),完全非線性數(shù)值數(shù)槽模型大多被應(yīng)用于純波浪問題.然而,在實際的工程問題中,波浪和水流往往同時存在,由于水流的存在,物面上波動勢的邊界條件發(fā)生變化,波浪對結(jié)構(gòu)物的繞射也相應(yīng)發(fā)生變化,從而波浪對結(jié)構(gòu)物的荷載以及在結(jié)構(gòu)物上的爬高也受到水流的影響.這樣再用純波浪理論和模型來計算就會導(dǎo)致很大的誤差,無法指導(dǎo)工程設(shè)計和作業(yè).在許多海洋工程項目中,準(zhǔn)確預(yù)測結(jié)構(gòu)物在波浪和水流共存環(huán)境下的荷載和運動響應(yīng)是非常重要的.國內(nèi)外很多專家學(xué)者在此方面也開展了大量的工作,目前比較成熟的研究大多還是基于線性和攝動展開理論,如采用頻域方法研究波流混合與結(jié)構(gòu)物的作用問題[6-8],采用時域方法研究波流混合與結(jié)構(gòu)物作用的線性問題[9].

高階邊界元法作為目前國際上先進的算法之一,已經(jīng)在很多領(lǐng)域得到廣泛應(yīng)用,特別是在非線性水波問題的時域模擬方面.然而在完全非線性波流混合作用方面的應(yīng)用還不多見.本文將采用時域高階邊界元方法建立波流混合作用的完全非線性數(shù)值水槽模型,研究波流混合作用與純波浪相比的波面變化及非線性變化情況.

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

1.1 控制方程和邊界條件

在假定流體無粘和流動無旋的前提下考慮波流混合傳播問題,整個流域可用速度勢來描述.建立一個笛卡爾坐標(biāo)系Oxyz如圖1所示,使得z=0位于靜水面上,且z軸向上為正,波浪傳播方向與 x軸正向保持一致.這里,自由水面、固體邊界分別用Sf、SN(包括水底Sd、水槽兩個側(cè)面Sc和入射邊界面 S1)表示.考慮到水流的作用,流域內(nèi)的整體速度勢可以表示為Φ=φ(x,y,z,t)+U0x,其中U0代表水流速度,φ(x,y,z,t)表示由波浪引起的速度勢,速度勢Φ和φ都滿足Laplace方程和下列邊界條件.

圖1 波流水槽示意圖Fig.1 Definition sketch of the wave-current tank

在瞬時自由水面Sf上,滿足完全非線性運動學(xué)和動力學(xué)邊界條件.本文采用混合歐拉-拉格朗日法更新自由水面,故自由水面邊界條件可以寫成以下形式[10]:

式中:XF(t)=(xF(t),yF(t),zF(t))為瞬時自由表面上任意流體質(zhì)點的位置矢量,zF=η(x,y,t)代表自由水面高程,g是重力加速度,物質(zhì)導(dǎo)數(shù)流體質(zhì)點速度v=▽Φ,也就是說波面隨著流體質(zhì)點的速度進行更新.

在水底和水槽側(cè)壁上滿足固壁不可滲透邊界條件

在水槽入射邊界上給定波浪傳播速度的二階斯托克斯理論解[11]:

式中:k是波數(shù),ω是波浪角頻率,h是水深,A是入射波幅.由于水流的存在,波數(shù)k和角頻率ω將滿足如下色散方程關(guān)系:

與其對應(yīng)的波流混合傳播的線性波面理論表達式為

從式中可以看出,當(dāng)水流速度為0時,式(4)和(5)與純波浪的表達式一致.

由于本問題是在時域內(nèi)求解,還要給定初始邊界條件,下列靜水面邊界條件被滿足:

在水槽的末端,采用在下游水面區(qū)加一人工阻尼層來吸收向右傳播的波浪,防止波浪發(fā)生反射,通過在自由面運動學(xué)和動力學(xué)邊界條件式(1)中加入阻尼項來實現(xiàn)[12],即

其中,

式中:X0=(x0,y0,0)是指水質(zhì)點靜止時的位置,x0為阻尼層起始坐標(biāo),α為阻尼系數(shù),βλ為阻尼層長度,λ為波長,β為岸灘寬度系數(shù),本文取α=1,β=1.

1.2 數(shù)值求解

在整個流域內(nèi)應(yīng)用格林第二定理,則上述邊值問題可轉(zhuǎn)化為如下的邊界積分方程:

式中:p=(x,y,z)為源點,q=(x0,y0,z0)為場點, α(p)為固角系數(shù).邊界S包括自由水面Sf和固體表面 SN,G是簡單格林函數(shù),考慮到計算域的對稱性,可以寫出如下形式:

其中,

這樣積分區(qū)域僅僅在一半的計算域上進行,且不包括水底.

本文用高階邊界元離散計算域成一些曲面單元,對每個單元通過數(shù)學(xué)變換,將其轉(zhuǎn)換成參數(shù)坐標(biāo)(ξ,ζ)下的等參單元,采用二次形狀函數(shù)插值方法保證單元內(nèi)物理量分布的連續(xù)性.單元內(nèi)任一點的幾何坐標(biāo)和速度勢等物理量都可以由各節(jié)點上的值和形狀函數(shù)表示成如下形式[13]:

通過配點法,將源點p分別取在各個節(jié)點上,可建立如下線性方程組:

由于積分邊界是不斷地隨著時間變化的,在每一計算時刻都要重新建立系數(shù)矩陣 A和B,并且在每一計算時刻都要對方程進行求解.

計算中認(rèn)為當(dāng)前時刻物面 SN上的速度勢法向?qū)?shù)和自由水面Sf上的速度勢是已知的,根據(jù)積分方程(9)計算當(dāng)前時刻物面SN上的速度勢和自由水面Sf上的法向速度,然后應(yīng)用四階Runge-Kutta法,根據(jù)自由水面條件式(7)計算下一時刻的水質(zhì)點位置和自由水面Sf上的速度勢,再對自由水面重新劃分網(wǎng)格,重新應(yīng)用積分方程計算下一時刻物面上的速度勢和自由水面上的法向速度.這樣計算周而復(fù)始,直到計算結(jié)束.

1.3 網(wǎng)格重新劃分

由于自由水面的變化,計算域的大小每一時間步都在改變,網(wǎng)格相應(yīng)地跟著改變,為了防止由于水質(zhì)點運動導(dǎo)致網(wǎng)格大小不均勻引起數(shù)值不穩(wěn)定問題,需要定期進行網(wǎng)格重新劃分來調(diào)整這種網(wǎng)格不均勻問題.網(wǎng)格重新劃分在直立面(如水槽入射面和側(cè)壁)非常容易,對于自由水面上的網(wǎng)格重新劃分,首先通過下式確定新節(jié)點計屬于重新劃分前哪一個單元:

式中:S0是未重新劃分單元的面積,Si是新節(jié)點與未重新劃分單元中 2個相鄰節(jié)點所組成的三角形的面積,M是環(huán)繞新節(jié)點所形成的三角形單元總數(shù),當(dāng)參數(shù) ε趨于 0時,認(rèn)為該節(jié)點隸屬于這個單元.然后求出新節(jié)點在所在舊單元中的形狀函數(shù),最后通過式(11)得到新節(jié)點處的空間位置和速度勢,詳細過程可參見文獻[14].

2 數(shù)值計算及討論

作為算例,考慮靜水深h=1.0m水槽中完全非線性波流混合傳播情況,研究波流同向、反向和純波浪3種工況.取波浪角頻率w=1.506 76,入射波幅A=0.05m,水流速度U0=0.313m/s.經(jīng)過開展數(shù)值收斂性試驗選定時間步長Δt=T/80(T為波浪周期),空間步長Δx=Δy=λ/20(λ為波長),水槽尺度定義為 5λ×0.1λ,垂向布置 4個單元,水槽末端布置長度為1.5λ的人工阻尼層.

從圖2中可以看出,波(流)傳到8T時已達到穩(wěn)定,3種工況下 8T和 10T的波面已經(jīng)完全重合,說明本模型已達到數(shù)值收斂,對波浪和波流混合問題都可以得到收斂的數(shù)值結(jié)果.

圖2 t=6T、8T和10T時水槽中線波面分布圖Fig.2 Distribution of wave elevation along the symmetry of wave flume at t=6T,8T and 10T

圖3～5是分別距離水槽入射邊界 5.0 m和10.0m兩點處波流同向、純波浪和波流反向 3種工況下的波面時間歷程及本文數(shù)值結(jié)果與已發(fā)表Boussinesq模型結(jié)果[10,15]的對比.從圖中可以看出,各點在各工況下的波面都保持穩(wěn)定傳播,且本文數(shù)值結(jié)果與Boussinesq模型結(jié)果吻合的很好,進而證明本模型的準(zhǔn)確性.

下面研究水流對波浪非線性的影響.圖 6是x=5.0 m和 10.0 m處,水流同向、無水流和波流反向時的波面時間歷程圖.從圖中可以看出,各工況下兩點處都經(jīng)歷一定的緩沖時間達到穩(wěn)定狀態(tài),波面幅值保持不變,但是由于本算例屬于淺水問題,各點之間存在一定的能量交換,以至于各點之間的波幅不同;波流同向與純波浪相比波浪的傳播有超前性,而波流反向與純波浪相比有滯后性,也就是說波浪同向時的傳播速度最小,即波數(shù)最大,這與色散關(guān)系式(5)得到的結(jié)論是一致的;3種工況下波流同向時波峰最小,波流反向時波峰最大,進而說明波流反向時使得波浪非線性增強,波流同向時非線性減弱.

圖3 波流同向時兩點(x=5.0m和10.0m)處的波面時間歷程及與Boussinesq模型結(jié)果對比Fig.3 Time series of non linear free surface wave with cop lanar current at x=5.0m and 10.0m and comparisons of present and Boussinesq model results.

圖4 無水流時兩點(x=5.0m和10.0m)處的波面時間歷程及與Boussinesq模型結(jié)果對比Fig.4 Time series of nonlinear free surface wave without current at x=5.0m and 10.0m and comparisons of present and Boussinesqmodel results.

圖5 波流反向時兩點(x=5.0m和10.0m)處的波面時間歷程及與Boussinesq模型結(jié)果對比Fig.5 Time series of nonlinear free surface wave with adverse currentat x=5.0m and 10.0m and comparisons of present and Boussinesqmodel results.

圖6 x=5.0m和10.0m兩點處的波面歷程圖Fig.6 Time histories of wave elevation at x=5.0m and 10.0m

圖 7是對圖 6中波面歷程進行傅里葉變換得到頻譜變化關(guān)系,圖中橫坐標(biāo)通過除以波浪頻率f0=1/T進行無量綱化處理.從圖中可以清楚看出波流同向、純波浪和波流反向 3種工況下非線性變化關(guān)系,以及 1階、2階甚至 3階波浪的貢獻.通過對比發(fā)現(xiàn)與圖 6相同的結(jié)論,波流反向時非線性最強,純波浪時次之,波流同向時非線性最弱.

圖7 x=5.0m和10.0m 3點處的頻譜變化關(guān)系Fig.7 Variantion ofwave speactra with frequencey at x=5.0m、and 10.0m

接下來是保持水流速度不變,研究不同入射波幅A情況下波流同向、純波浪和波流反向時的波面變化情況.圖8是x=5.0m和10.0m兩點在波流混合作用下的波峰 A*隨波陡的變化情況,圖中縱坐標(biāo)通過除以入射波幅 A做無量綱化處理,為了進行比較,各工況下線性理論結(jié)果也在圖中給出.從圖中可以看出,非線性結(jié)果都隨著波陡的增加而增大,但其增加的幅度波流反向情況最為明顯,純波浪次之,波流同向時最小,這與圖7中非線性分析的結(jié)論是一致的.

圖9是x=5.0 m和10.0m兩點處的平均波高一半 A*隨波陡的變化情況圖,以及線性結(jié)果和非線性結(jié)果的比較.從圖中可以看出與圖 8相同的變化規(guī)律,但是變化幅度明顯小于波峰情況,這說明波谷的起到減弱的作用.線性和非線性結(jié)果的區(qū)別也進一步說明,對于非線性較強情況,無論是波流反向還是波流同向,線性理論值都會導(dǎo)致很大的不準(zhǔn)確性.

圖 8 兩點處不同波陡情況下波峰的變化關(guān)系Fig.8 Variation ofwave crestwith wave slope at two positions

圖 9 兩點處不同波陡情況下平均波高一半的變化關(guān)系Fig.9 Variation ofaverge wave height with wave slope at two positions

3 結(jié)束語

本文基于時域高階邊界元方法建立波流混合作用的完全非線性數(shù)值水槽模型,并且通過不同水流條件下與Boussinesq方程模型結(jié)果對比進行驗證,結(jié)果表明本文所建立數(shù)學(xué)模型可以準(zhǔn)確模擬強非線性條件下波流混合作用問題,這是線性理論解很難做到的,特別是對于強非線性情況.研究發(fā)現(xiàn),相對于純波浪問題,波流反向作用時波峰變得更高和陡峭,波谷變得更低而平坦,也即波浪的非線性特點更加突出,對于波流同向作用則得到相反的結(jié)論,對于如本文中所研究的淺水問題,在水槽下游各點處會體現(xiàn)出更強的非線性特性,特別是波流反向時.本文模型的建立將有助于進一步開展波流混合與海洋工程結(jié)構(gòu)物作用問題的模擬研究,為海洋工程安全設(shè)計提供參考.

[1]SCORPIO S M.Fully nonlinear ship-wave computations using amultipole accelerated,desingularizedmethod[D]. Michigan:The University ofMichigan,1997:1-106.

[2]WU G X,HU Z Z.Simulation of non linear interactions between waves and floating bodies through a finite-elementbased numerical tank[J].Proc R Soc Lond A,2004, 460:797-817.

[3]寧德志,滕斌,周斌珍,等.源造波法在三維完全非線性波浪水槽中的應(yīng)用[J].大連海事大學(xué)學(xué)報,2008,34:1-5.

[4]NING D Z,TENGB,EATOCK TAYLOR R,etal.Numerical simulation ofnonlinear regular and focused waves in an infinite water-depth[J].Ocean Engineering,2008,35:887-899.

[5]BAIW,EATOCK T R.Numerical simulation of fully nonlinear regu lar and focused wave diffraction around a vertical cylinder using domain decomposition[J].Applied Ocean Research,2007,29(1-2):55-71.

[6]NOSSEN J,GRUE J,PALM E.Wave forces on 3-D floating bodies with samll forward speed[J].J Fluid Mech, 1991,242:135-160.

[7]TENG B,EATOCK TAYLORR.App lication of a higher order BEM in the calculation of wave run-up in a weak current [C]//Proceedings of the Conference ISOPE.Osaka,Japan, 1994.

[8]滕斌.深水中波浪與弱流對結(jié)構(gòu)物的作用[J].海洋學(xué)報,1996,18:117-125.

TENG Bin.Wave and weak current actions on object in deep water[J].Acta Oceanologica Sinica,1996,18:117-125.

[9]ISAACSONM,NG Y T.Time-domain Second-order wave interaction with 3-D floating bodies[J].International Journal of Offshore and Polar Engineering,1995,5(3):171-179.

[10]RYU S,KIM M H,LYNETT P J.Fully nonlinear wavecurrent interactions and kinematics by a BEM based numerical wave tank[J].Computational Mechanics,2003, 32:336-346.

[11]KOOW,KIM M.Current effects on nonlinear wave body interactions by a 2-D fully non linear numerical wave tank [J].Journalof Waterway,Port,Coastaland Ocean Engineering,2007,133(2):136-146.

[12]李玉成,滕斌.波浪對海上建筑物的作用 [M].2版,北京:海洋出版社,2002:1-295.

[13]寧德志.快速多極子邊界元方法在完全非線性水波問題中的應(yīng)用[D].大連:大連理工大學(xué),2005:1-166.

NING Dezhi.Application of the fastmulitpole boundary elementmethod to fu lly nonlinear water wave problems[D]. Dalian:Dalian University of Technology,2005:1-166.

[14]周斌珍,寧德志,滕斌.造波板運動造波實時模擬[J].水動力學(xué)研究與進展:A輯,2009,24(4):1-12.

ZHOU Binzhen,NING Dezhi,TENG Bin.Realtime simulaiton of waves generated by a wavemaker[J].Chinese Journalof Hydrodynamics:Ser.A,2009,24(4):1-12.

[15]KOO W,KIM M H.Current effects on non linear wavebody interactions by a 2-D fully nonlinear numerical wave tank[J].Journal of Waterway,Port,Coastal and Ocean Engineering,2007,133(2):136-146.