吳 明，應(yīng)榮熔，蔡 烽，楊 波，王 驍

(1. 海軍大連艦艇學(xué)院航海系，遼寧 大連 116018；2. 海軍大連艦艇學(xué)院教務(wù)處，遼寧 大連 116018)

0 引 言

波浪水池是船舶與海洋工程問題研究的基礎(chǔ)試驗設(shè)施，在許多復(fù)雜水動力學(xué)問題研究中發(fā)揮著重要的作用。近年來，隨著計算機技術(shù)的進(jìn)步，以計算流體動力學(xué)方法為基礎(chǔ)的數(shù)值波浪水池技術(shù)得到了蓬勃發(fā)展，該技術(shù)在計算機上即可仿真物理模型水池的功能，并且可獲得豐富的流場細(xì)節(jié)信息，這使其逐漸在波浪演化[1-6]、船舶性能評估[7-12]等問題的研究中發(fā)揮重要作用，已成為當(dāng)前水動力學(xué)技術(shù)發(fā)展的重要方向之一。

數(shù)值波浪水池的造波精度，是衡量其性能的標(biāo)準(zhǔn)。在構(gòu)建數(shù)值波浪水池的過程中，發(fā)現(xiàn)波浪在傳播過程中非常容易出現(xiàn)衰減，即由數(shù)值模擬得到的波浪時歷，經(jīng)譜分析得到的譜密度值明顯小于相應(yīng)的目標(biāo)

譜的譜密度值，而減小的幅度遠(yuǎn)遠(yuǎn)超過了由于流體粘性而導(dǎo)致的物理衰減。文獻(xiàn)[13]對規(guī)則波的數(shù)值衰減問題進(jìn)行了探討，給出了在數(shù)值波浪水池中生成高精度規(guī)則波的方法。但不規(guī)則波的構(gòu)成更加復(fù)雜，該方法不完全適用于不規(guī)則波的模擬。本文從對不規(guī)則波模擬精度影響較大的幾個因素：網(wǎng)格尺寸、時間步長、欠松弛因子和頻率等分?jǐn)?shù)量4個方面，對不規(guī)則波數(shù)值模擬的精度進(jìn)行研究，提出不規(guī)則波數(shù)值模擬的建議。

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

1.1 RANS方程

波浪的傳播屬于低速流動，因此波浪的演化可視為粘性不可壓縮流體流動，遵循質(zhì)量守恒方程和動量守恒方程[7]：

式中：ρ為流體密度；p為流體壓強；ui為流體時均速度分量；fi為流體體積力分量；v為相對于時均流速的湍流脈動速度分量；ν為流體的運動粘性系數(shù)；為雷諾應(yīng)力。

1.2 湍流方程

上述方程組并不封閉，因此還需要建立與雷諾應(yīng)力有關(guān)的表達(dá)式，即建立湍流方程來封閉方程組求解，本文采用的是RNGk-ε湍流模型[14]，其表達(dá)式為：

1.3 自由面方程

采用VOF（Volume of Fluid）方法[15]模擬波面的傳播。該方法用體積分?jǐn)?shù)Cq表示網(wǎng)格單元內(nèi)某項流體體積占總體積的比例，其輸運方程為：

式中：q1為空氣項；q2為水項。

2 數(shù)值波浪水池的構(gòu)建

2.1 數(shù)值造波

采用邊界條件造波法生成波浪，即在造波邊界處給出波速和波高的解析解實現(xiàn)數(shù)值造波[7]。根據(jù)線性疊加原理，長峰不規(guī)則波的波高和速度方程為：

式中：Ai，ωi，εi，ki，S（ωi）分別為第i個子波的波幅，波浪圓頻率，初始相位，波數(shù)和波浪頻譜密度；U為流速。

2.2 數(shù)值消波

采用阻尼消波方法消波，該方法在數(shù)值水池末端區(qū)域加載的動量方程中增加阻尼源項[9]，可以有效耗散不同波長和頻率波浪的能量。為了增強消波效果，在消波區(qū)內(nèi)沿布設(shè)的網(wǎng)格沿波浪傳播方向逐漸稀疏，這樣可以進(jìn)一步加快波浪的數(shù)值耗散。阻尼消波區(qū)的動量方程為：

式中：μ（x）為一單調(diào)遞增函數(shù)，消波區(qū)內(nèi)越靠近出口消波功能越強。

3 不規(guī)則波模擬精度影響因素分析

為了便于對影響不規(guī)則波數(shù)值模擬精度的各個因素進(jìn)行研究，建立一個長12 m，高4 m，深3 m的二維數(shù)值波浪水池，如圖1所示。選擇ITTC雙參數(shù)譜作為目標(biāo)譜，不規(guī)則波的三一波高為0.1 m，平均周期為1.22 s，平均波長為1.55 m，在數(shù)值波浪水池的工作區(qū)內(nèi)距離入口4 m的位置，設(shè)置一個波高時歷監(jiān)測點。

圖1 數(shù)值波浪水池計算域Fig. 1 Computation domain of numerical wave tank

3.1 網(wǎng)格尺寸影響

網(wǎng)格是影響數(shù)值模擬精度的關(guān)鍵因素之一，為了研究網(wǎng)格尺寸對不規(guī)則波數(shù)值模擬精度的影響，進(jìn)行兩組對比試驗，在第1組數(shù)值試驗中，進(jìn)行波長方向網(wǎng)格數(shù)目對不規(guī)則波數(shù)值模擬精度的影響研究。在此基礎(chǔ)上，在第2組數(shù)值試驗中保持波長方向的網(wǎng)格數(shù)目，進(jìn)行波高方向不同網(wǎng)格數(shù)目的影響研究，最終優(yōu)選出適合不規(guī)則波數(shù)值模擬的網(wǎng)格尺寸。

在第1組數(shù)值試驗中，保持三一波高內(nèi)的網(wǎng)格數(shù)目為12個，在平均波長內(nèi)分別布設(shè)20，30，40個網(wǎng)格。圖2(a)為采用3種不同波長方向網(wǎng)格尺寸，所獲得的波高時歷與理論值的對比，其中線1為模擬值，線2為理論值。采用周期圖法對波浪時歷進(jìn)行譜分析，圖2(b)為波浪譜的對比圖。由對比圖可以看出，平均波長內(nèi)網(wǎng)格數(shù)Nx=20時，模擬譜與目標(biāo)譜誤差很大，平均波長內(nèi)網(wǎng)格數(shù)Nx=30時，模擬譜與目標(biāo)譜吻合較好，平均波長內(nèi)網(wǎng)格數(shù)Nx=40時，模擬譜與目標(biāo)譜吻合最好。

圖2 波長方向網(wǎng)格數(shù)的影響Fig. 2 The influence of grid number along wave length direction

在第2組數(shù)值試驗中，保持三一波長內(nèi)的網(wǎng)格數(shù)目為40個，在三一波高內(nèi)分別布設(shè)6，12，20個網(wǎng)格。圖3(a)為采用3種不同波高方向網(wǎng)格尺寸所獲取的波高時歷，圖3(b)為波浪譜的對比圖。由對比圖可以看出，當(dāng)ω＜5時，3種情況下的仿真精度相當(dāng)，當(dāng)ω＞5時，三一波高內(nèi)布設(shè)6個網(wǎng)格時，模擬譜有一定衰減，布設(shè)12個網(wǎng)格時，模擬譜與目標(biāo)譜吻合較好，布設(shè)20個網(wǎng)格時，模擬譜與目標(biāo)譜吻合度最好。

3.2 時間步長的影響

由于數(shù)值造波問題需要采用全隱式時間積分方案，該方案在時間域上只具有一階截差精度，因此需要使用較小的時間步長，以保證獲得較高精度的解。為了研究時間步長對不規(guī)則波數(shù)值模擬精度的影響，進(jìn)行不同時間步長的對比試驗。試驗中保持平均波長內(nèi)的網(wǎng)格數(shù)目為40個，三一波高內(nèi)的網(wǎng)格數(shù)目為12個，時間步長分別為0.000 5 s，0.001 s，0.002 s。

圖4(a)為采用3種不同時間步長所獲得的波高時歷，圖4(b)為波浪譜的對比圖。由對比圖可以看出，當(dāng)時間步長為0.002 s時，模擬譜的中高頻部分有一定的衰減，時間步長為0.001 s和0.000 5 s時，模擬譜與目標(biāo)譜均吻合較好，其中時間步長為0.001 s時，目標(biāo)譜高頻部分略有衰減，時間步長為0.000 5 s時，模擬譜與目標(biāo)譜吻合最好。

3.3 欠松弛因子的影響

欠松弛技術(shù)是使數(shù)值方程穩(wěn)定地趨于收斂的一種有效方法。欠松弛因子α的取值范圍為0～1。總體來言，增大的α值可以加快收斂速度，但若α取值過大時，可能導(dǎo)致系統(tǒng)很難得到穩(wěn)定的解，甚至出現(xiàn)發(fā)散。減小α值會使計算的穩(wěn)定性增加，但是若α取得過小，數(shù)值計算收斂的速度可能過慢，甚至影響計算的精度。波浪數(shù)值模擬屬于氣液兩相流問題，因此，體積分?jǐn)?shù)欠松弛因子對數(shù)值造波的精度具有一定影響，為了研究壓力和體積分?jǐn)?shù)欠松弛因子對不規(guī)則波數(shù)值模擬精度的影響，進(jìn)行2組對比試驗，分別選擇不同的壓力欠松弛因子和體積分?jǐn)?shù)欠松弛因子進(jìn)行模擬。

在第1組數(shù)值試驗中，將壓力欠松弛因子設(shè)置為0.5，而體積分?jǐn)?shù)欠松弛因子分別取0.3，0.5，0.8。圖5(a)為采用3種不同體積分?jǐn)?shù)欠松弛因子所獲取的波高時歷，圖5(b)為波浪譜的對比圖。由對比圖可以看出，當(dāng)體積分?jǐn)?shù)欠松弛因子為0.3時，模擬譜中頻部分衰減比較明顯，當(dāng)體積分?jǐn)?shù)欠松弛因子為0.5和0.8時，模擬譜與目標(biāo)譜均吻合良好，且此時得到的2種模擬譜幾乎沒有差別，這說明體積分?jǐn)?shù)欠松弛因子大于0.5時，即可保證數(shù)值仿真的精度。

取體積分?jǐn)?shù)欠松弛因子為0.5，而壓力欠松弛因子分別取0.1，0.2，0.3。圖6(a)為采用3種不同壓力欠松弛因子所獲得的波高時歷。圖6(b)為波浪譜的對比圖。由對比圖可以看出，當(dāng)壓力欠松弛因子為0.1時，模擬譜中頻部分產(chǎn)生了一定的誤差；當(dāng)壓力欠松弛因子為0.2和0.3時，模擬譜與目標(biāo)譜均吻合良好，且此時得到的2種模擬譜幾乎沒有差別，這說明壓力欠松弛因子大于0.2時，即可保證數(shù)值造波的精度。

圖4 時間步長的影響Fig. 4 The influence of time step

圖5 積分?jǐn)?shù)欠松弛因子的影響Fig. 5 The influence of relaxation factor of volume fraction

3.4 頻率等分?jǐn)?shù)量的影響

海浪譜的頻率范圍為0～∞，海浪譜能量大部分集中在某一頻率范圍，數(shù)值模擬時，需要根據(jù)線性疊加原理，將該頻率范圍分割為一定數(shù)量的區(qū)間，通常采用頻率等分法進(jìn)行分割，頻率等分?jǐn)?shù)量過少可能會影響數(shù)值造波的精度，過多會增加計算時間，因此需要對頻率等分?jǐn)?shù)量開展比較研究。分別將頻率等分為：15,30和60份，數(shù)值模擬中其他的參數(shù)設(shè)置為：網(wǎng)格數(shù)目在平均波長布設(shè)40個，在三一波高內(nèi)布設(shè)20個，時間步長設(shè)置為0.001 s，體積分?jǐn)?shù)欠松弛因子和壓力欠松弛因子分別設(shè)置為0.5和0.3。圖7(a)為設(shè)置3種不同頻率等分?jǐn)?shù)量所獲得的波高時歷，圖7(b)為波浪譜的對比圖。從圖中可以看出，頻率等分為15個時，模擬譜中低頻部分有所增大，中頻部分有所減?。活l率等分為30和60時，模擬譜與目標(biāo)譜均吻合良好，這說明頻率等分?jǐn)?shù)量超過30時，即可保證不規(guī)則波數(shù)值模擬的精度。

圖6 壓力欠松弛因子的影響Fig. 6 The influence of relaxation factor of pressure

3.5 不規(guī)則波造波模式

通過上述對比研究可以看出，通過加密網(wǎng)格、縮小時間步長、增大欠松弛因子，增加頻率等分?jǐn)?shù)量，可以提高使數(shù)值水池生成不規(guī)則波的精度，但數(shù)值造波的精度并不會隨著上述因素的細(xì)化和提高而不斷提高，即精度的提高是有上限的。與此同時，過度細(xì)化網(wǎng)格、縮小時間步長、增加頻率等分?jǐn)?shù)量，會使數(shù)值計算時間明顯增加，且數(shù)值模擬中采用過大的欠松弛因子還可能導(dǎo)致數(shù)值計算的不穩(wěn)定。因此，應(yīng)根據(jù)需求平衡計算精度和計算效率關(guān)系，選用合適的數(shù)值計算參數(shù)。根據(jù)上述的對比試驗結(jié)果，提出以下幾點不規(guī)則波數(shù)值模擬的建議：

1）網(wǎng)格數(shù)量在平均波長方向為30～40個，三一波高方向為12～20個；

2）時間步長取值為波浪平均周期的1/1 000～1/2000；

3）壓力和體積分?jǐn)?shù)欠松弛因子分別大于0.2和0.5；

4）頻率等分區(qū)間數(shù)量為30～60個。

圖7 頻率等分?jǐn)?shù)量的影響Fig. 7 The influence of the number of discrete frequencies

4 驗 證

為了驗證本文數(shù)值模擬模式的有效性，選擇模擬ITTC雙參數(shù)譜的模型尺度長峰不規(guī)則波和實測譜的實尺度長峰不規(guī)則波進(jìn)行數(shù)值模擬，與目標(biāo)譜參數(shù)進(jìn)行比較，驗證試驗的波浪參數(shù)設(shè)置如表1所示。

表1 波浪譜參數(shù)Tab. 1 Wave parameters

圖8為模擬的長峰不規(guī)則波的瞬時波面波高分布圖。算例1將監(jiān)測點設(shè)置在工作區(qū)距入口4 m和6 m處，算例2將監(jiān)測點設(shè)置在工作區(qū)距離入口40 m和60 m處，圖9(a)和圖10(a)為所測得的波高時歷，圖9(b)和圖10(b)為兩算例的波浪譜的對比圖。通過對比圖可以看出，模擬值與目標(biāo)值吻合較好。

中華人民共和國交通部發(fā)布JTJ/T 234-2001《波浪模型試驗規(guī)程》[16]規(guī)定了單項不規(guī)則波模擬的允許偏差：

1）波譜總能量允許偏差為±10%；

2）譜峰頻率模擬值允許偏差為±5%；

圖8 長峰不規(guī)則波瞬時波高圖Fig. 8 Instantaneous picture of irregular wave height

圖9 算例1模擬結(jié)果Fig. 9 Sumulation result of case1

3）在譜密度大于或等于0.5倍譜密度峰值的范圍內(nèi)，譜密度分布允許偏差為±15%

4）有效波高、有效周期或譜峰周期允許偏差為±5%。

對監(jiān)測的波面時歷進(jìn)行譜分析，可以得到本文數(shù)值模擬長峰不規(guī)則波的各項統(tǒng)計特征值，與目標(biāo)值的對比如表2所示。表中所列數(shù)據(jù)為工作區(qū)2個監(jiān)測點數(shù)據(jù)平均值，其中算例2實測譜有2個譜峰值，因此對比了2處峰值的譜峰頻率。從表2可以看出，數(shù)值模擬的三一波高略偏小，波浪平均周期及譜峰頻率與目標(biāo)值吻合良好，上述3項特征值的相對誤差均在4%以下，譜面積的相對誤差低于6%。與規(guī)范中的標(biāo)準(zhǔn)相比，本文的模擬相關(guān)特征值均低于規(guī)范允許偏差。另外，對于實測譜的模擬，本文模擬再現(xiàn)了雙峰特征，并且具有較高的精度，說明本文提出的數(shù)值模擬方案同樣適用于實測海浪的模擬。

圖10 算例2模擬結(jié)果Fig. 10 Sumulation result of case2

表2 數(shù)值仿真誤差Tab. 2 Relative errors of simulation result

5 結(jié) 語

本文采用計算流體動力學(xué)方法，在數(shù)值水池中生成長峰不規(guī)則波。在此基礎(chǔ)上，對影響長峰不規(guī)則波數(shù)值模擬精度的幾項因素：網(wǎng)格尺寸、時間步長、欠松弛因子及頻率等分?jǐn)?shù)量開展對比試驗，研究不規(guī)則波數(shù)值仿真精度隨上述因素的變化規(guī)律，給出了不規(guī)則波數(shù)值模擬的建議。最后，采用本文提出的方法，生成了模型尺度和實尺度的不規(guī)則波，符合相關(guān)行業(yè)標(biāo)準(zhǔn)的要求，驗證了本文提出的數(shù)值模擬模式的有效性。