梁修鋒，楊建民，李 欣，李 俊

（上海交通大學(xué)海洋工程國家重點(diǎn)實(shí)驗(yàn)室，上海200030）

不規(guī)則波浪的數(shù)值模擬

梁修鋒，楊建民，李 欣，李 俊

（上海交通大學(xué)海洋工程國家重點(diǎn)實(shí)驗(yàn)室，上海200030）

優(yōu)良的波浪環(huán)境條件是研究非線性波浪和浮式結(jié)構(gòu)物相互作用問題的基礎(chǔ)。文章在數(shù)值波浪水池中對(duì)不規(guī)則波浪進(jìn)行了數(shù)值模擬，造波通過搖板的運(yùn)動(dòng)實(shí)現(xiàn)，消波通過在水池后段設(shè)置消波區(qū)實(shí)現(xiàn)，自由液面由VOF方法捕捉。對(duì)計(jì)算所得的波浪時(shí)歷進(jìn)行譜分析并將所得譜與理論譜進(jìn)行了比較，吻合良好。

數(shù)值波浪水池；動(dòng)網(wǎng)格；VOF方法；JONSWAP譜

1 引 言

隨著計(jì)算機(jī)技術(shù)的發(fā)展，建立數(shù)值水池應(yīng)用于數(shù)值模擬的條件越來越成熟。與物模實(shí)驗(yàn)相比數(shù)值水池具有經(jīng)濟(jì)性高、無觸點(diǎn)測(cè)量流場(chǎng)、比尺效應(yīng)小、消除物模中由傳感器尺寸及模型變形等因素對(duì)流場(chǎng)的影響、可獲得較詳細(xì)的流場(chǎng)信息等優(yōu)點(diǎn)。但是目前數(shù)值水池的端部反射問題，尤其是造波板的二次反射限制了數(shù)值水池的應(yīng)用。通常的處理方法是在二次反射到達(dá)結(jié)構(gòu)物之前停止實(shí)驗(yàn)，或模型寬度設(shè)計(jì)得比水池的寬度小很多，對(duì)數(shù)值模擬意味著需要大量增加計(jì)算區(qū)域，這受到計(jì)算機(jī)內(nèi)存及計(jì)算機(jī)速度的限制。王永學(xué)[1]基于線性造波機(jī)理論在水池的一端設(shè)置了可吸收式造波邊界條件，造波板的運(yùn)動(dòng)除了產(chǎn)生行進(jìn)波外，同時(shí)還產(chǎn)生一個(gè)抵消反射波的局部波動(dòng)。Brorsen和Larsen[2]提出了適合于邊界積分方程方法（BIEM）的源造波方法（source generation），即在計(jì)算域內(nèi)設(shè)置一造波源，在源項(xiàng)兩邊同時(shí)產(chǎn)生方向相反的兩列波，源項(xiàng)處可透過波浪遇建筑物形成的反射波。

為了消除邊界處的反射，Larsen和Dancy[3]提出了用海綿層方法吸收傳遞到邊界處的波能，通過海綿層消除波能的大部分，同時(shí)在出流邊界處利用Sommerfled條件，使未能衰減的部分波浪透過邊界傳到域外。孫大鵬和李玉成[4]利用上述思想采用時(shí)域內(nèi)對(duì)波面運(yùn)動(dòng)位置追蹤的邊界元方法，建立了一種非線性波浪變形計(jì)算的三維數(shù)值模式，開發(fā)了一個(gè)三維非線性波的數(shù)值造波水池，進(jìn)而對(duì)水池內(nèi)的Stokes波進(jìn)行了波浪變形計(jì)算。高學(xué)平等[5]利用MAC（marker and cell）法直接數(shù)值求解連續(xù)方程和NS方程，為模擬不規(guī)則波和長(zhǎng)時(shí)段內(nèi)連續(xù)造波及消除波浪遇結(jié)構(gòu)物形成的反射，采用了源造波方法，在開敞邊界處采用了海綿層阻尼消波和Sommerfled條件相結(jié)合的處理方式。

本文首先推導(dǎo)出用于不規(guī)則波浪模擬的搖板信號(hào)，結(jié)合FLUENT軟件中的用戶自定義函數(shù)（UDF）功能和動(dòng)網(wǎng)格（Dynamic Mesh）功能進(jìn)行造波，同時(shí)利用FLUENT軟件中可以加入源項(xiàng)（Source Terms）的特征，在水池后段加入了消波源項(xiàng)，使得水質(zhì)點(diǎn)速度在水池的后段逐步衰減為零，從而消除了反射波對(duì)水池內(nèi)波形的影響，為長(zhǎng)時(shí)間的波浪模擬提供了保證，自由液面由兩相流中的VOF方法捕捉。基于以上對(duì)不規(guī)則波浪的數(shù)值模擬進(jìn)行了研究，得到了高品質(zhì)的數(shù)值波浪。借助于FLUENT軟件可以提供復(fù)雜網(wǎng)格的特性，該水池可以為復(fù)雜浮式結(jié)構(gòu)物水動(dòng)力特性的模擬提供數(shù)值波浪環(huán)境，具有廣闊的應(yīng)用前景。

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

控制方程包括連續(xù)性方程：

上式中，u、v分別為 x、y方向的速度，ρ為流體的密度，μ 為動(dòng)力學(xué)粘性系數(shù)，f1=0，f2=-g，F(xiàn)1和 F2為附加源項(xiàng)。

本文采用VOF方法捕捉自由液面，因此還需要引入另外一個(gè)控制方程，即流體體積分?jǐn)?shù)的輸運(yùn)方程

本文所采用的Geometric Reconstruction方法是由Youngs[6]的方法發(fā)展而來的。該方法可以用于非結(jié)構(gòu)網(wǎng)格，計(jì)算精度較高。幾何重構(gòu)的第一步是計(jì)算交界面在網(wǎng)格內(nèi)相對(duì)于網(wǎng)格中心的位置，第二步由速度矢量和交界面的位置確定通過網(wǎng)格面的流體流量，第三步由前兩步的結(jié)果計(jì)算每一個(gè)網(wǎng)格的體積分?jǐn)?shù)。它在每一時(shí)間步內(nèi)根據(jù)網(wǎng)格內(nèi)的流體體積分?jǐn)?shù)將網(wǎng)格標(biāo)記為流體網(wǎng)格（aq=1），空網(wǎng)格（aq=0）或自由液面網(wǎng)格（0＜aq＜1）。并根據(jù)自由液面網(wǎng)格內(nèi)的體積分?jǐn)?shù)aq重構(gòu)自由液面的幾何形狀。然后利用離散后的傳輸方程確定體積分?jǐn)?shù)得到下一個(gè)時(shí)間步的自由液面形狀，基于交錯(cuò)網(wǎng)格或同位網(wǎng)格求

文中q=1表示空氣相，q=2表示水相。aq是一個(gè)標(biāo)量，表示q相流體在網(wǎng)格內(nèi)占的體積分?jǐn)?shù)解不可壓流體的N-S方程，得到下一個(gè)時(shí)間步流域內(nèi)各點(diǎn)的速度值和壓力值，如此循環(huán)而求出流體在各個(gè)時(shí)刻的運(yùn)動(dòng)情形。本文默認(rèn)體積分?jǐn)?shù)aq=0.5處為自由液面。

3 幾何模型

采用GAMBIT建立計(jì)算所需的幾何模型，見圖1。坐標(biāo)原點(diǎn)O位于水池左端的初始自由液面處，水池長(zhǎng)40m，水深2m，消波區(qū)長(zhǎng)6m，初始自由液面以上的高度為0.4m，造波板繞R點(diǎn)做定軸轉(zhuǎn)動(dòng)。在x=15m處設(shè)置了虛擬浪高儀器以監(jiān)測(cè)波浪升高時(shí)歷。流場(chǎng)網(wǎng)格劃分如圖2所示，在流場(chǎng)的大部分區(qū)域都采用了結(jié)構(gòu)網(wǎng)格，這樣既可以減少網(wǎng)格數(shù)量，同時(shí)又可以得到光順的自由液面，在自由液面附近加密了網(wǎng)格以提高計(jì)算精度，見圖2，自由液面附近網(wǎng)格大小為0.015m，縱向?yàn)?.03m，網(wǎng)格總數(shù)約為75 000。

4 搖板的運(yùn)動(dòng)規(guī)律

本文的目標(biāo)波波譜為JONSWAP譜[7]，水深為128m，有義波高為14.3m，譜峰周期為15.3s，譜型參數(shù)γ=2.4，縮尺比為64:1，計(jì)算時(shí)采用了模型值，最后將計(jì)算結(jié)果換算成了實(shí)型值。

由波高和波浪周期表示的譜公式為：

式中a=0.062 4/[0.023 0+0.336γ-0.185（ 1.9+γ ）-1]；ω≤ωP時(shí) σ=0.07，ω＞ωP時(shí) σ=0.09。

零階譜距（譜密度曲線下的面積）為：

為了得到不規(guī)則波需要對(duì)規(guī)則波進(jìn)行波浪疊加[7]，劃分頻率區(qū)間的方法有等分頻率法和等分能量法兩種，本文采用等分能量法。定義累積譜為

如按等分能量法分成M份，則分界頻率ωi可按式

為了給定造波板的運(yùn)動(dòng)規(guī)律，需要知道水力傳遞函數(shù)，即不同周期時(shí)板前的波浪振幅a0與造波板振幅e的關(guān)系。

為了適用于一般情況，假定造波板的運(yùn)動(dòng)為推搖混合，距離池底距離d1，d2處板的振幅分別為e1，e2，水面處振幅為 e，見圖 3，由勢(shì)流理論可得[8]。

對(duì)于搖板，e1=e，e2=0，d1=d，d2=0，則

本文采用搖板方式進(jìn)行造波，因此采用（11）式進(jìn)行搖板振幅的計(jì)算，并利用圖3的關(guān)系求解搖板的幅角 θi=arctan（ e/d），各個(gè)頻率對(duì)應(yīng)的搖板初相位 εi在（0～2π ）內(nèi)隨機(jī)取值，見圖5。 最后將各個(gè)頻率下的搖板運(yùn)動(dòng)進(jìn)行疊加即可得到搖板在任意時(shí)刻的角速度信號(hào)

5 邊界條件、初始條件及求解方案

搖板為動(dòng)邊界，模擬JONSWAP譜時(shí)，其運(yùn)動(dòng)規(guī)律見（12）式，頂部為壓力入口邊界條件，其余各邊為無滑移的固壁邊界條件（見圖1），初始時(shí)刻水池內(nèi)的水體呈靜水壓力分布

流場(chǎng)采用兩相流中的VOF模塊和層流模式求解，壓力—速度項(xiàng)采用PISO算法進(jìn)行迭代計(jì)算，動(dòng)量方程中的瞬態(tài)項(xiàng)采用一階隱格式，為了同時(shí)保證計(jì)算方案的穩(wěn)定性和精確性，本文對(duì)流項(xiàng)和擴(kuò)散項(xiàng)的離散在最初的五個(gè)時(shí)間步內(nèi)采用一階迎風(fēng)格式，此后的計(jì)算采用二階迎風(fēng)格式，時(shí)間步進(jìn)方案為NITA算法。由于造波板附近主要采用結(jié)構(gòu)網(wǎng)格（見圖2），并且造波板做周期性的搖動(dòng)，所以采用動(dòng)網(wǎng)格模塊中的鋪層（Layering）功能進(jìn)行網(wǎng)格更新來實(shí)現(xiàn)造波板的運(yùn)動(dòng)。最后將動(dòng)邊界的運(yùn)動(dòng)方程及消波區(qū)中的源項(xiàng)[9-10]寫成用戶自定義函數(shù)UDF在FLUENT軟件里的相關(guān)接口激活進(jìn)行造波和消波。

6 計(jì)算結(jié)果與分析

按照前述的原理和計(jì)算方法，對(duì)JONSWAP譜進(jìn)行了數(shù)值模擬，然后對(duì)所得的波浪時(shí)歷進(jìn)行譜分析，最初水池中生成的不規(guī)則波波譜不一定是所尋求的目標(biāo)譜，若兩者差異較大，就需要按下式進(jìn)行迭代修正[11]：

7 結(jié) 論

本文基于FLUENT中的用戶自定義函數(shù)（UDF）功能，采用搖板方式進(jìn)行造波，在水池后段的動(dòng)量方程中添加了源項(xiàng)進(jìn)行消波，從而建立了數(shù)值波浪水池。對(duì)JONSWAP譜進(jìn)行數(shù)值模擬的結(jié)果表明，消波區(qū)可以很好地減小反射波對(duì)水池內(nèi)波形的影響，所得波浪的運(yùn)動(dòng)規(guī)律與理論值比較接近。該水池可以非常方便地?cái)U(kuò)展到三維情形，由于計(jì)算時(shí)間的限制本文沒有進(jìn)行三維隨機(jī)波浪的模擬。借助于FLUENT軟件強(qiáng)大的建模和后處理功能，該水池可以為結(jié)構(gòu)復(fù)雜的海洋浮式結(jié)構(gòu)物的水動(dòng)力特性的數(shù)值模擬提供數(shù)值波浪環(huán)境，具有廣泛的應(yīng)用前景。

[1]王永學(xué).無反射造波數(shù)值波浪水池[J].水動(dòng)力學(xué)研究與進(jìn)展,1994,19(2):205-214.

[2]Brorsen M,Larsen J.Source generation of nonlinear gravity waves with the boundary integral equation method[J].Coastal Engineering,1987,11:93-113.

[3]Larsen J,Dancy H.Open boundaries in short wave simulation-a new approach[J].Coast.Eng.,1983,7:285-297.

[4]孫大鵬,李玉成.數(shù)值水池內(nèi)的阻尼消波和波浪變形計(jì)算[J].海洋工程,2000,18(2):46-50.

[5]高學(xué)平,曾廣東,張 亞.不規(guī)則波浪數(shù)值水池的造波和阻尼消波[J].海洋學(xué)報(bào),2002,24(2):127-132.

[6]Youngs D L.Time-dependent multi-material flow with large fluid distortion[C]//Proceedings of Numerical Methods for Fluid Dynamics.New York,1982.

[7]盛振邦,劉應(yīng)中.船舶原理[M].上海:上海交通大學(xué)出版社,2004:360-362.

[8]孫昭晨.推搖混合式造波機(jī)理論曲線[J].港口工程,1988(4):30-32.

[9]梁修鋒,楊建民,李 欣,于 洋.FPSO甲板上浪數(shù)值模擬[J].水動(dòng)力學(xué)研究與進(jìn)展,A輯,2007,22(2):229-236.

[10]周勤俊,王本龍,蘭雅梅,劉 樺.海浪越浪的數(shù)值模擬[J].力學(xué)季刊,2005,26(4):629-633.

[11]俞湘三,陳澤梁,樓連根,遲云鵬.船舶性能試驗(yàn)技術(shù)[M].上海:上海交通大學(xué)出版社,1991:188.

Numerical simulation of irregular wave

LIANG Xiu-feng,YANG Jian-min,LI xin,LI jun

(State Key Laboratory of Ocean Engineering,Shanghai Jiao Tong University,Shanghai 200030,China)

Quality of wave environment is the basis for studying the nonlinear interactions between water waves and floating structures.In this paper,a numerical wave tank was used to simulate the irregular wave.Flap type wave-maker was adopted for generating waves and wave reflection was handled properly through setting up wave absorption zone.VOF method was employed for tracking the free surface.Spectrum analysis was performed on time trace of waves and the calculated spectrum agreed well with the theoretical spectrum.

numerical wave tank;dynamic mesh;VOF method;JONSWAP spectrum

U661.7

A

1007-7294（2010）05-0481-06

2009-12-24

深水油氣開發(fā)工程試驗(yàn)技術(shù)（2006AA09A107）；自然科學(xué)基金（50709079）

梁修鋒（1981-），男，博士研究生，E-mail:liangxiufeng@sjtu.edu.cn；

楊建民（聯(lián)系人），男，教授，博士生導(dǎo)師，E-mail:jmyang@sjtu.edu.cn。