波浪水槽中Stokes五階波的數值生成

李邦華　鄭向遠　李　煒　榮維棟

(大連理工大學船舶工程學院1)　大連　11602)　(清華大學深圳研究生院2)　深圳　518055)

(中國電建集團華東勘測設計研究院有限公司3)　杭州　310014)

?

波浪水槽中Stokes五階波的數值生成

李邦華1)鄭向遠2)李煒3)榮維棟1)

(大連理工大學船舶工程學院1)大連11602)(清華大學深圳研究生院2)深圳518055)

(中國電建集團華東勘測設計研究院有限公司3)杭州310014)

摘要：在Fluent軟件平臺下運用邊界造波法，通過給定造波邊界處流體的五階速度和波面瞬時升高，實現了在數值波浪水槽中使用流體體積VOF法對Stokes五階波浪的精確模擬，得到的波浪與傳統理論值在時域和頻域內均吻合很好.此外，為了更好地確定邊界水質點輸入參數，對現有的若干種五階波理論解與流函數解在相同的波浪參數下進行對比,發現波浪參數在Ursell判據所定義的五階波使用范圍內時，理論解在波長、波速、波面瞬時升高等參數上區別甚微，但在接近破碎和橢圓余弦波的波浪條件下，波長和波速有較大的差別.

關鍵詞：斯托克斯五階波；數值波浪水槽；邊界造波法；VOF方法；流函數

0引言

在計算流體力學飛速發展的背景下，通過在計算機上建立成本低、易于改造和拓廣的二維或三維的多相數值波浪水槽已經得到了大量的使用.

劉霞等[1]運用FLUENT軟件中的UDF功能,通過定義造波邊界處波浪的速度和波形，曾成功模擬了二階Stokes波的產生和傳播.楊錦凌等[2]利用動量源方法開發了能夠消除二次反射波的消波方法.廉靜靜等[3]建立了三維規則波數值波浪水槽，通過搖板造波方法得到了質量較好的波形.路寬等[4]在三維數值波浪水槽中建立閉環控制系統，解決了數值水槽造波過程中輸入端與采集端數據不相符的問題.陸萍[5]使用源函數造波法并用VOF法對自由表面進行追蹤，成功地模擬出不同工況下的線性規則波.

對于海洋結構物波浪荷載的計算有時采用線性波理論或者低階的非線性波理論不能達到工程精度的要求.竺艷蓉[6]根據海洋結構物的受力特性，通過理論計算與水槽試驗結果進行對比分析，認為使用Stokes五階波理論可以獲得比低階波理論更為精確的結果，因此在數值波浪水槽中生成五階波是十分必要的.Skjelbreia等[7]發表的論文中，斯托克斯五階波公式有微小的錯誤，然而很多國內外出版物中均直接采用此公式，并沒有對錯誤項進行修正.當進行數值波浪水槽模擬時若直接采用此公式，將引起波形和水質點速度等誤差，從而影響水中結構物所受的波浪力.因此本文在二維和三維數值水槽中進行了Stokes五階波的生成，并饒有興趣地將修正前后的公式(Nishimura等[8]，邱強[9])、Fenton[10]傅里葉級數法,以及流函數法進行了對比.

1控制方程

對于求解三維不可壓縮，密度為常數的波浪運動，可將控制方程簡化如下.

根據質量守恒得到的連續性方程：

(1)

根據牛頓第二定律得到的動量方程：

(2)

式中：u，v，w分別為速度矢量x，y，z方向上的分量；Sx，Sy，Sz分別為x，y，z方向上的動量源項；p為流體內部壓強；μ為粘度；ρ為流體的密度;t為時間.

2VOF方法

流體體積法[11]其核心思想是通過網格單元的流體體積比φq來準確地追蹤流體的自由液面.φq是標量，表示第q相流體在網格單元中所占的體積分數.當φq=1時，則表示該網格單元內全部為第q相流體；當φq=0時，則表示該網格單元內沒有第q相流體；若0<φq<1時，則表示該網格單元中有一部分為q相流體.φq需滿足如下方程：

(3)

(4)

式中：u，v，w分別為速度矢量在x，y，z方向上的分量；t為時間.對于追蹤數值波浪水槽里的自由液面，只需設置氣液兩相即可.

3造波與消波

3.1邊界造波法

本文在FLUENT軟件平臺上進行了二次開發，通過給定造波邊界處的流體在x和z方向上的速度以及波面瞬時升高實現了波浪的生成.對于Stokes五階波浪，在造波邊界處的速度和波面瞬時升高滿足以下條件.

x方向速度：

(5)

z方向速度：

(6)

波面瞬時升高：

(7)

其中各項系數如下：

式中：ω，d，k分別為圓頻率、水深和波數.定義c=coshkd；s=sinhkd，其余各項系數參見文獻[7]，但應注意其中的C2系數在1977年被Nishimura等指出有錯誤.因此本文在計算時采用文獻[8]修正后的C2系數.

3.2消波方法

本文采用阻尼消波法，該方法是通過在特定的計算區域內的動量方程中添加阻尼項來削弱或消除該區域的波動.添加阻尼項后的動量方程為

(8)

式中：μ(x)為阻尼系數，沿波浪傳播方向為線性分布；υ為運動粘度；α為經驗系數；x1和x2為消波區起始位置和終止位置在x軸上的坐標值.

4波浪的數值模擬

4.1數值水槽模型的建立與網格劃分

本文中首先在二維數值波浪水槽實現了五階波的模擬，但由于二維數值水槽不能對具有復雜三維幾何構造的結構物進行實驗，因此又建立了三維數值波浪水槽，見圖1.為了減少網格數量、提高計算效率，對于本文所要生成的五階波(波浪參數為水深1.2 m，周期1 s，波高0.123 4 m、波長1.647 2 m)，造波區的長度約一個波長取1.7 m，試驗區的長度根據模型的大小進行設置，本文取0.8 m，過渡區的長度為0.7 m，消波區的長度約兩個波長取3.2 m.建立了二維和三維數值波浪水槽的幾何模型后，需要對其進行網格劃分，為了在試驗區得到較好模擬結果，本文在造波區和試驗區以及靜水面上下一個波高范圍內設置了較密的網格(見圖2)，所得到的二維和三維數值水槽的網格數量分別為21 952和439 040.

圖1　數值波浪水槽的縱向圖

圖2　數值波浪水槽網格劃分圖(側視圖)

4.2離散以及求解方法

本文中的數值波浪水槽使用了層流模型，水的密度取為1 024 kg/m3，粘度取0.001 003 Pa·s.對于計算域內偏微分方程組的離散，Fluent軟件使用了計算效率較高的有限體積法，對于控制體積界面上的物理量及其導數的插值本文均采用二階迎風格式.運用SIMPLE算法的分離式解法對離散后的方程組進行求解，求解時的時間步長取0.002 5 s.

4.3模擬結果與分析

為了驗證數值波浪水槽的造波和消波效果，在數值水槽的造波區x=1.4 m、試驗區x=2 m、消波區x=6.4 m處均設置了虛擬浪高儀來監測波面變化.由于二維和三維數值水槽得到的波形基本相同(圖3)，因此下文中的比較均使用了三維數值波浪水槽得到的結果.首先，將試驗區距離造波端2 m處虛擬浪高儀得到的波面時域曲線與文獻[8]得到的理論解進行了對比，見圖4.圖4中比較的是若干秒后的穩定波形；其次，本文進一步分析其頻譜特性(見圖5)，檢驗其頻率含量與相應的能量是否與理論值吻合；最后，圖6中比較的是造波區和消波區的波形，以檢驗消波效果如何.

圖3　二維與三維數值波浪水槽在試驗區的波形比較

圖4　造波區的波面瞬時升高與理論值的在時域上的比較

圖5　造波區波面瞬時升高與理論值在頻域譜上的比較

圖6　造波區與消波區波浪比較

為了定量對比數值水槽生成的五階波浪與Stokes五階波修正后解的差距，在此引入以下誤差評價指標對數值水槽生成的五階波浪的精確度進行評價.

1) 峰值誤差Erp是指數值水槽生成的波浪在穩定后的波面瞬時升高Eta′(t)的峰值和谷值相對于理論值Eta(t)的峰值和谷值的2個相對誤差的求和平均值，數學定義如下.

(9A)

2) 絕對誤差Ers是指數值水槽生成的波浪的波面瞬時升高Eta′(t)與理論值Eta(t)之差的絕對值之和與理論值Eta(t)的絕對值之和的比值，數學定義如下.

(9B)

3) 波峰的平均誤差Eap指數值水槽生成的波浪在穩定后的波峰值Wp′(i)與理論解的波峰值Wp之差的求和平均值與理論解波峰值Wp的比值.

(9C)

4) 波谷的平均誤差Eav指數值水槽生成的波浪在穩定后的波谷值Wv′(i)與理論解的波谷值Wv之差的求和平均值與理論解波峰值Wv的比值.

(9D)

使用4種誤差評價指標對數值水槽生成的五階波的波面瞬時升高與理論值的誤差進行了定量的分析，見表1.

表1　數值水槽模擬與理論值的誤差分析

由圖4、圖5，以及表1可知，數值波浪水槽生成的五階波與理論值在周期、波高等波浪參數上基本相同并且誤差較小，在頻域上的頻譜特性也十分一致，尤其是從圖5b)中可以看到5個顯著的峰值，對應5個不同頻率的諧波.表1中的誤差主要來自于粘性作用，而傳統的五階波解大多基于無粘的攝動解.由圖6造波區與消波區的波面瞬時升高對比知，所使用的阻尼消波能夠在消波區將波浪進行有效的衰減，消波效果良好.

5五階波理論解的比較

對于邊界造波法，輸入的水質點速度與波面瞬時升高等參數直接影響著水槽造波質量的好壞.文獻[7]認為文獻[8]中的C2系數有誤，即“+2 592”實際上是“-2 952”； 文獻[9]也按照Skjelbreia提出的方法進行了推導并修正了C2系數；文獻[10]將流函數展開為傅里葉級數的形式，并將波陡kH/2作為攝動參數進行展開得到了五階波的解.文獻[12]提出了流函數波浪理論，與Stokes五階波浪理論相比，能夠更好的滿足邊界條件，因此計算結果最為精確.本文對以上幾種五階波公式與流函數波浪理論在8種具有代表性的波浪參數下得到的理論解進行了對比，見表2.其中，RW7為的波浪參數為文獻[7]中波浪接近破碎和橢圓余弦波時的參數，RW8和RW3的波陡(在流函數解下)為最大和次大.此外，圖7～11更為顯著的展示了各浪理論解在波面上的差異.

圖7　RW3波浪條件下波面瞬時升高的對比(T=4 s，d=16 m，H=3.11 m)

圖8　RW5波浪條件下波面瞬時升高的對比(T=7 s，d=3.9 m，H=1 m)

圖9　RW6波浪條件下波面瞬時升高的對比(T=7 s，d=3 m，H=0.8 m)

對于各種波浪理論的適用范圍，很多學者進行了相應的研究.Le Méhauté[14]通過實驗畫出了目前工業界廣泛使用的各種波浪理論的適用范圍圖，見圖12. Fenton[15-16]使用Ursell判據得到了Stokes與橢圓余弦波浪理論的分界線，從而代替了1990年所確定的分界線.本文在Fention所確定的波浪理論適用范圍圖(見圖13)中畫出了以上2種界線,其中Fention(1990)和Ursell=40分別表示他在1990年和1998年給出的分界線.

表2　8種波浪參數下五階波與流函數的比較

注：T為周期；d為水深；H為波高；L為波長；H/L為波陡；ΔL為與用流函數波浪理論得到的波長值之差；ΔPeak為與用流函數波浪理論得到的波峰值之差；Δux為與用流函數波浪理論得到的水質點在x方向的最大速度值之差；Δuz為與用流函數波浪理論得到的水質點在z方向的最大速度值之差.

圖10　RW7波浪條件下波面瞬時升高的對比(T=7.72 s，d=9.144 m，H=5.689 6 m)

圖11　RW8波浪條件下波面瞬時升高的對比(T=6.1 s，d=13 m，H=6.5 m)

文獻[14]給出的波浪破碎界限公式為

淺水區：

(10)

有限水深采用Michell公式：

(11)

深水區：

(12)

文獻[14]給出的Stokes與橢圓余弦波浪理論的分界線公式為

(13)

文獻[15-16]擬合出的波浪破碎界限公式為

(14)

文獻[15]給出的Stokes與橢圓余弦波浪理論的分界線公式為

(15)

文獻[16]給出的Stokes與橢圓余弦波浪理論的分界線公式為

(16)

式中：H，L，k，d分別為為波高、波長、波數和水深.

本文將RW1～RW8所代表的波浪參數表示在以上2種常用的波浪理論適用范圍圖中(見圖12、圖13)，每種工況都有相應的幾何形狀.在圖13中，波長值均使用流函數波浪理論計算得到的波長值.

圖12　波浪理論適用范圍圖(Le Méhauté 1969)

圖13　波浪理論適用范圍圖(Fenton)

1) 由表2和圖7，11，12可見，波浪參數在文獻[14]得到的斯托克斯波范圍內時，Skjelbreia原文中和修正后的五階波公式得到的理論解區別甚微，同時與用Fenton推出的五階波公式得到的解也沒有太大區別，但用Fenton推出的五階波公式得到的解與用流函數波浪理論得到的解更為接近，對于RW1，2，3，4，8五種波浪條件均是如此.圖11表明，由于RW8的波陡大且接近破碎界限，各種五階波公式解的精度有所降低，但仍在工程誤差范圍內.

2) 圖11再度表明RW3在斯托克斯波的適用范圍內，雖然其波陡為次大，但在圖7中由各種斯托克斯五階波理論得到的解和流函數得到的解有著很高的吻合度.不同于圖12，RW5在圖13中落入斯托克斯波的范圍內，但從圖8中良好的吻合度可以看出， 用Ursell=40作為斯托克斯波與橢圓余弦波的分界線對RW5來說更為合理，而圖12的波浪理論適用范圍圖偏保守.

3) 對于RW6，圖13的Ursell=40判據表明其落入橢圓余弦波的范圍，而Fenton(1990)的判據卻將其納入斯托克斯波的范圍內，從圖9中和表2中的差異可以看出，Fenton得到的五階波公式對RW6并不適用，而其余3種五階波依然適用.

4) 對于RW7，從圖9中和表2中的差異可以看出，Fenton的五階波解與流函數解相差就更大了，其余3種五階波解也有不小的誤差，但修正后的斯托克斯五階波公式得到的解在波長和波速上誤差更小.

綜上所述，使用圖13中Ursell=40作為斯托克斯波與橢圓余弦波的分界線是最為合理的，其右域為橢圓余弦波，其左域為斯托克斯波.其次，Fenton的五階波理論僅適用于落入此判據左側的波浪參數，甚至不適合其1990年的直線判據.

6結論

1) 本文以FLUENT軟件為計算平臺，采用VOF法進行自由液面的追蹤，使用斯托克斯五階波浪理論、邊界造波法和阻尼消波法，建立了有效的數值波浪水槽并生成了精度良好的斯托克斯五階波浪.

2) 波浪參數在文獻[14]通過實驗確定的五階波范圍內時，Skjelbreia原文中和修正后的五階波公式得到的波理論解區別甚微，Fenton導出的五階波公式的解與流函數解最為接近.在接近破碎和橢圓余弦波的波浪參數下，幾種五階波公式得到的解雖然在波面上相差不大，但波長和波速以及波陡有較大的差別.Fenton導出的五階波公式甚至不適用于其1990年定義的斯托克斯波范圍.因此，在運用數值水槽進行模擬時，比如波浪對淺水中海洋風力發電機基礎的作用，針對不同的波浪參數應選擇正確的五階波公式或直接使用流函數計算出的水質點速度與波形作為造波端的邊界輸入，此時應使用Ursell=40作為斯托克斯波與橢圓余弦波的分界線.

參 考 文 獻

[1]劉霞，譚國煥，王大國.基于邊界造波法的二階Stokes波的數值生成[J].遼寧工程技術大學學報(自然科學版),2010(1)：107-111.

[2]楊錦凌，孫大鵬.基于FLUENT二次開發的數值波浪水槽[J].中國水運，2012(5):59-61.

[3]廉靜靜，尹勇，楊曉，等.基于粘性流船舶數值波浪水池造波和消波方法研究[J].船舶力學,2013(Z1)：56-62.

[4]路寬，齊連明，楊寧，等.閉環控制在波浪數值水槽中的應用[J].海洋技術學報,2014(1):80-83.

[5]陸萍.基于源函數造波法三維數值水槽波浪模擬研究[J].廣東造船,2014(6)：45-47.

[6]竺艷蓉.幾種波浪理論適用范圍的分析[J].海岸工程,1983(2)：11-27.

[7]SKJELBREIA L, HENDRICKSON J. Fifth order gracity wave theory[C]. Proceedings 7th Conference of Coastal Engineer,1961：184-196.

[8]NISHIMURA H, ISOBE M， HORIKAWA K. Higher order solutions of the stokes and cnoidalwaves[J]. Journal of the Faculty of Engineering. Univ. Tokyo, Series B,1977(2)：267-293.

[9]邱強.關于五階重力波系數C2的一點注記:更正文獻[1]中系數C2的表達式[J].海洋學報,1983(4)：539-541.

[10]FENTON J D. A fifth-order Stokes theory for steady waves[J]. Waterway Port Coastaland Ocean Engineering ASCE,1985,111:216-234.

[11]HIRT C W, HICHOLS B D. Volume of fluid (VOF) method for dynamics of free boundaries[J]. Journal of Computational Physics,1981,39:201-221.

[12]WOLF J P. Soil-structure interaction analysis in time domain[M]. New Jersey: Prentice-Hall,1987.

[13]DEAN R G. Stream function representation of nonlinear waveforce[J]. Journal of Geophysical Research,1965,70(18):4651-4572.

[14]LE MéHAUTéB. An introduction to hydrodynamics and water waves[R]. Water Wave Theories, U.S. Department of Commerce, ESSA, Washington, D.C.,1990,Vol. II, TR ERL 118-POL-3-2.

[15]FENTON J D. Non-linear wave theories[M]. Wiley：New York,1990.

[16]FENTON J D. The cnoidal theory of water waves [M].Gulf: Houston,1998.

Simulation of Fifth-order Stokes Waves in the Numerical Wave Flume

LI Banghua1)ZHENG Xiangyuan2)LI Wei3)RONG Weidong1)

(FacultyofVehicleEngineeringandMechanics,

DalianUniversityofTechnology,Dalian116024,China)1)(ShenzhenGraduateSchool,TsinghuaUniversity,Shenzhen518055,China)2)

(POWERCHINAHuadongEngineeringCorporation,OffshoreWindPowerR&DCenterof

POWERCHINAHuadong,Hangzhou310014,China)3)

Abstract：Using FLUENT, simulation of the fifth-order Stokes waves is achieved in the numerical wave flume by defining the inlet boundary conditions and using VOF solver. Comparison with theoretical wave profiles in both time- and frequency-domain shows that the simulated waves are accurate. In order to better determine the inlet water particle kinematics, the existing several models of fifth-order Stokes waves are compared with the stream function method being the most accurate method. These nonlinear solutions all agree well with each other when the wave conditions fall into the Stokes range specified by the Ursell criteria. However, they differ significantly from each other when waves approach the breaking limit and elliptical cosine wave conditions.

Key words：fifth-order stokes wave; numerical wave flume; wave-generation by defining inlet boundary conditions; VOF method; stream function

doi:10.3963/j.issn.2095-3844.2016.02.008

中圖法分類號：U66

收稿日期：2016-01-17

李邦華(1992- )：男，碩士生,主要研究領域為海洋結構物結構分析

*國家自然科學基金項目(51379035，41206075)、中國電建集團華東勘測設計研究院有限公司“波浪在單立柱三樁式風機基礎結構上的爬升及對工作平臺的砰擊作用研究”項目資助