多周期穩(wěn)定數(shù)值造波的參數(shù)設(shè)計(jì)方法

呂磊，陳作鋼

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

2 高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240

3 上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海 200240

0 引 言

在船舶與海洋工程領(lǐng)域，利用計(jì)算流體力學(xué)（CFD）方法研究波浪中結(jié)構(gòu)物的各類響應(yīng)一直是重要的研究方向。然而，在進(jìn)行CFD計(jì)算時(shí)，數(shù)值模擬收斂需要一定的物理時(shí)間，但數(shù)值波在計(jì)算過(guò)程中會(huì)發(fā)生波幅衰減和相位偏差的問(wèn)題，影響結(jié)構(gòu)物響應(yīng)的計(jì)算精度。因此，實(shí)現(xiàn)多周期穩(wěn)定的數(shù)值造波顯得尤為重要。

目前，較為常用的數(shù)值造波方法為有限體積法（finite volume method，F(xiàn)VM），通過(guò)求解雷諾平均（Reynolds averaged Navier-Stokes，RANS）數(shù)值方程，并結(jié)合流體體積（volume of fluid，VOF）多相流模型來(lái)模擬波浪傳播過(guò)程。對(duì)于此類數(shù)值造波問(wèn)題，國(guó)內(nèi)學(xué)者進(jìn)行了廣泛研究。例如：楊云濤等[1]基于不可壓縮黏性流體的Navier-Stokes方程，采用質(zhì)量源造波法（mass source-based wave generation）建立了三維無(wú)反射數(shù)值波浪水池，克服了邊界造波法無(wú)法消除二次反射波的困難，并成功模擬了波浪與浮體相互作用的問(wèn)題。王東旭等[2]選擇兩相不可壓縮黏性流求解器 InterFoam作為平臺(tái)進(jìn)行了數(shù)值波浪水槽的開發(fā)，使用動(dòng)邊界造波法和質(zhì)量源造波法進(jìn)行了規(guī)則波和不規(guī)則波的模擬計(jì)算，結(jié)果表明，所開發(fā)的數(shù)值波浪水槽具有較高的精度與良好的造波能力。余陽(yáng)等[3]基于推板造波理論和搖板造波理論，在Open-FOAM平臺(tái)上采用重疊網(wǎng)格技術(shù)建立黏性數(shù)值波浪水槽，對(duì)比研究了此數(shù)值模型分別嵌入層流和湍流兩種模型后的計(jì)算精度及效率，并對(duì)三維數(shù)值造波問(wèn)題進(jìn)行了進(jìn)一步研究，驗(yàn)證了在三維兩相流體域中求解運(yùn)動(dòng)物體與流場(chǎng)交互的可靠性和正確性。上述研究從各方面探討了數(shù)值造波問(wèn)題，并發(fā)現(xiàn)數(shù)值波在傳播過(guò)程中普遍存在波幅衰減和相位偏移的現(xiàn)象。總之，包括網(wǎng)格設(shè)置在內(nèi)的多個(gè)參數(shù)共同影響了數(shù)值造波的精度，而對(duì)于數(shù)值造波誤差的評(píng)估及造成誤差的關(guān)鍵參數(shù)仍缺乏定量的系統(tǒng)性研究。

因此，本文將基于STAR-CCM+商業(yè)求解器進(jìn)行數(shù)值模擬，以五階Stokes波作為研究對(duì)象，首先將CFD計(jì)算結(jié)果與理論解進(jìn)行對(duì)比，通過(guò)對(duì)一維波動(dòng)方程的離散形式的分析，得到影響造波精度的幾個(gè)關(guān)鍵參數(shù)，進(jìn)一步研究這些參數(shù)對(duì)造波效果的影響。然后，通過(guò)算例計(jì)算，驗(yàn)證內(nèi)迭代次數(shù)的選擇、網(wǎng)格分辨率的設(shè)置及庫(kù)朗數(shù)等對(duì)波幅衰減和相位偏移的影響，并合理選擇參數(shù)來(lái)顯著提高多周期造波精度，最后，總結(jié)得出一套合理的參數(shù)設(shè)計(jì)方法，基于此方法，以KCS(KRISO container ship)船型為例，模擬計(jì)算其在迎浪中阻力和運(yùn)動(dòng)響應(yīng)，再通過(guò)與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比，以驗(yàn)證本文參數(shù)設(shè)計(jì)方法的可行性。

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

1.1 控制方程

計(jì)算的控制方程為Navier-Stokes方程，包括連續(xù)性方程和動(dòng)量守恒方程，如式（1）和式（2）所示。式中：V為控制體；S為控制體的邊界；v為流體的速度矢量；ui速度矢量的笛卡爾分量；n為垂直于邊界S并指向邊界外的單位矢量；t為時(shí)間；p為壓力；ρ為流體的密度；τij為黏性應(yīng)力張量的分量；ii和ij分別為xi方向和xj方向的單位矢量；g為體積力，在此即為重力；qi為一個(gè)可選的動(dòng)量源項(xiàng)，將在之后的消波設(shè)置中具體說(shuō)明。

1.2 湍流模型

在對(duì)多周期的波浪進(jìn)行數(shù)值模擬時(shí)，隨著計(jì)算物理時(shí)間的增加，自由液面處的湍流黏度急劇升高，此時(shí)會(huì)對(duì)波幅產(chǎn)生削減作用，且不同湍流模型對(duì)波浪傳播產(chǎn)生的負(fù)面影響也不盡相同。因此，本文使用標(biāo)準(zhǔn)的低雷諾數(shù)K-Epsilon模型，將其中的湍動(dòng)能和湍流耗散率設(shè)為非常低的初始值，這會(huì)有助于改善湍流黏度對(duì)波浪傳播的不利影響[4]。本文所有算例也均采用該K-Epsilon模型，根據(jù)文獻(xiàn)[4] 湍動(dòng)能和湍流耗散率的初始值分別設(shè)為：1×10?5J/kg和1×10?4m2/s3，與默認(rèn)值相比分別小了2和3個(gè)數(shù)量級(jí)。

1.3 自由液面捕捉

使用VOF法追蹤流體自由液面的位置。該方法使用流體體積比ai表 示第i相流體在網(wǎng)格單元中所占的體積分?jǐn)?shù)。每個(gè)網(wǎng)格單元中的體積分?jǐn)?shù)ai之和必須等于1。體積分?jǐn)?shù)的輸運(yùn)方程為

式中，qai為體積分?jǐn)?shù)ai的一個(gè)體積源。

1.4 造 波

造波方法采用邊界輸入法。定義縱向坐標(biāo)x對(duì)應(yīng)于波浪傳播方向，垂向坐標(biāo)z為鉛垂向上的方向。根據(jù)五階Stokes波的解析解，在計(jì)算域的邊界分別設(shè)定水質(zhì)點(diǎn)x方向（波傳播方向）速度ux，z方向速度uz和波面瞬時(shí)高度η，如下所示。

式中：ω為圓頻率；c=coshkd，其中，d為水深，k為波數(shù)；n為Stokes波階數(shù)。其他符號(hào)含義見文獻(xiàn)[5]。

1.5 消 波

在對(duì)波浪傳播進(jìn)行數(shù)值模擬時(shí)，因計(jì)算域范圍有限，在各邊界處常會(huì)發(fā)生波浪反射的情況，影響造波精度，因此，在研究時(shí)使用了阻尼消波法消除邊界處的反射波，以使模擬計(jì)算環(huán)境更接近于真實(shí)情況。該方法在動(dòng)量守恒方程 (2) 中將動(dòng)量源項(xiàng)qi定義如下：

式中：fl為線性阻尼系數(shù)；f2為二次阻尼系數(shù);uz為垂向坐標(biāo)z方向的速度分量；eκ?1/e1?1為彎曲項(xiàng)，其中κ為彎曲指數(shù)；xsd和xed分別為阻尼層開始和結(jié)束時(shí)的x坐標(biāo)（波傳播方向），一般xed位于計(jì)算域末端，故阻尼消波區(qū)的長(zhǎng)度為xd=|xed?xsd|；nd控制阻尼消波區(qū)的波形衰減形狀，一般取nd=2。

為實(shí)現(xiàn)理想的消波效果，xd,fl和f2需根據(jù)波浪參數(shù)進(jìn)行精確調(diào)整[6]。推薦只使用線性消波項(xiàng)（f2=0），阻尼消波區(qū)的長(zhǎng)度設(shè)為2倍波長(zhǎng)（xd=2λ），線性阻尼常數(shù)f1=Ψ1ω，其中系數(shù)Ψ1=π。

1.6 離散格式及其他設(shè)置

本文使用有限體積法（finite volume method,FVM）將控制方程離散化。所有數(shù)值差分、微分和積分都是基于二階近似。每個(gè)代數(shù)方程包含來(lái)自網(wǎng)格中心及所有相鄰網(wǎng)格中心的未知值，獲得的耦合方程組被線性化，并由隱式非定常求解器求解。計(jì)算時(shí)，使用SIMPLE方法獲取壓力值并校正速度。壓力的松弛因子為0.4，其他所有變量的松弛因子均為0.9。使用高分辨率接觸面捕捉（high-resolution interface capturing，HRIC）方法提高自由液面的捕捉精度。同時(shí)，激活梯度平滑，減少在大網(wǎng)格長(zhǎng)高比下的波面振蕩。

2 誤差分析

數(shù)值模擬的誤差來(lái)源主要有：模型誤差，即數(shù)學(xué)模型理論解與物理情況實(shí)際值間的誤差；離散誤差，即離散方程的解與數(shù)學(xué)模型理論解間的誤差；迭代誤差，即迭代算法的解與離散方程解間的誤差。

本研究只考慮實(shí)際的數(shù)值造波結(jié)果與理論解間的差異，故無(wú)需考慮模型誤差。而通過(guò)對(duì)線性方程組的求解和利用SIMPLE法進(jìn)行足夠次數(shù)的迭代，迭代誤差可遠(yuǎn)遠(yuǎn)小于離散誤差[7]，因而，離散誤差是主要誤差來(lái)源。為便于討論，本文參考一維波動(dòng)方程來(lái)分析影響離散誤差的具體參數(shù)。波面控制方程為：

若以顯式一階迎風(fēng)格式對(duì)波面方程(9)進(jìn)行有限差分，離散誤差則為式(10)等號(hào)右端項(xiàng)之和。

式中：φ為速度勢(shì)；c為波速；σ為庫(kù)朗數(shù)（Courant–Friedrichs–Lewy number），σ=cΔt/Δx，其中，Δx為網(wǎng)格尺寸，Δt為時(shí)間步長(zhǎng)。由式(10)可知，波浪傳播的離散誤差主要由2個(gè)無(wú)量綱數(shù)影響：

1） 庫(kù)朗數(shù)σ，即單個(gè)時(shí)間步內(nèi)波浪傳播的網(wǎng)格數(shù)。

2）與波浪尺度相比的幾何分辨率Δx/λ（其中λ為波長(zhǎng)）。當(dāng)Δx越小，幾何分辨率越高，單個(gè)網(wǎng)格內(nèi)的梯度就會(huì)越小，離散誤差也越小[8]。

離散誤差的具體形式會(huì)因具體方程及所用的離散格式而不同，但無(wú)論哪種形式，誤差的大小都與庫(kù)朗數(shù)和幾何分辨率有關(guān)。對(duì)于二維造波問(wèn)題來(lái)說(shuō)，庫(kù)朗數(shù)定義為：

式中，uz,max為z方向的速度最大值。幾何分辨率則與長(zhǎng)度 方 向x和高度 方 向z的 尺度比Δx/λ和Δz/H有關(guān)（其中H為波高，λ 為波長(zhǎng)）。為控制網(wǎng)格幾何分辨率的大小，若定義波面附近加密區(qū)z方向的網(wǎng)格數(shù)為Nz，網(wǎng)格長(zhǎng)高比為r，則有Δz=H/Nz，Δx=rΔz=rH/Nz。

當(dāng)網(wǎng)格尺寸Δz和Δx確定后，可通過(guò)改變時(shí)間步長(zhǎng)Δt來(lái)確定庫(kù)朗數(shù)大小，具體如式（12）所示。

綜上，后續(xù)研究將圍繞最小加密區(qū)z方向的網(wǎng)格數(shù)Nz、網(wǎng)格長(zhǎng)高比r以及庫(kù)朗數(shù)σ這些參數(shù)展開。造波誤差則由波幅誤差em和相位偏差ep來(lái)衡量[9]。定義波幅的理論解為E=|E|ejθ，數(shù)值解為E′=|E′|ejθ（其中θ 為極坐標(biāo)系旋轉(zhuǎn)角，θ?為其估計(jì)值），則坐標(biāo)點(diǎn)(xN，zN)(N表示數(shù)據(jù)點(diǎn))的波幅百分比誤差em和相位偏差ep(單位：rad)為：

在后處理分析誤差時(shí)，分別對(duì)各位置波幅誤差em在波周期 (0～10)T，(10～20)T，(20～30)T和(30～40)T之間的時(shí)刻進(jìn)行時(shí)間平均，相位偏差ep則主要關(guān)注每個(gè)位置在t=10，20，30，40T時(shí)刻的相位偏差ep。

3 算例設(shè)置

3.1 計(jì)算域設(shè)置

計(jì)算模擬二維深水波情況，計(jì)算域尺寸根據(jù)所對(duì)應(yīng)的波浪波長(zhǎng)λ進(jìn)行調(diào)整，計(jì)算域如圖1所示。為研究波浪在空間上傳播可能存在的衰減，計(jì)算域長(zhǎng)度Lx=10λ。為模擬深水條件，消除上、下邊界對(duì)流場(chǎng)的影響，計(jì)算域高度Lz=6λ，其中水深d=4λ，其他為空氣相。坐標(biāo)原點(diǎn)設(shè)在計(jì)算域入口邊界靜水液面高度處，縱向坐標(biāo)x對(duì)應(yīng)于波浪傳播方向，垂向坐標(biāo)z對(duì)應(yīng)于鉛垂向上的方向。

圖1 計(jì)算域設(shè)置示意圖Fig.1 Schematic view of computational domain setting

入口邊界（x=0）設(shè)為速度入口，數(shù)值波采用五階Stokes波形式；出口邊界（x=Lx）設(shè)為壓力出口，并向前設(shè)置xd=2λ長(zhǎng)度的消波區(qū)；底部邊界（z=?4λ）設(shè)為滑移壁面，對(duì)應(yīng)底部不可穿透條件；頂部邊界（z=2λ）設(shè)為速度入口，對(duì)應(yīng)于開放的頂部條件。

3.2 網(wǎng)格設(shè)置

為了能夠精確地捕捉數(shù)值波的交界面，需對(duì)波面附近網(wǎng)格進(jìn)行局部加密。這里，選擇波高方向網(wǎng)格數(shù)Nz和網(wǎng)格長(zhǎng)高比r=Δx/Δz作為網(wǎng)格設(shè)計(jì)策略的參數(shù)。如圖2所示，在自由液面附近沿波浪傳播方向共設(shè)置5個(gè)網(wǎng)格加密區(qū)。

圖2 網(wǎng)格加密區(qū)示意圖Fig.2 Schematic view of the zones by grid refinement

1) 加密區(qū)①，高度h1=1.25H，使加密區(qū)覆蓋整個(gè)波高范圍。加密區(qū)內(nèi)使用各向異性網(wǎng)格尺寸，高度方向的網(wǎng)格尺寸Δz=1/Nz；長(zhǎng)度方向的網(wǎng)格尺寸Δx取Δz的r倍。Nz和r的取值將在下文說(shuō)明。

2) 加密區(qū)②和③的高度h2,h3分別為2.5H和5H，網(wǎng)格尺寸分別為加密區(qū)①的2和4倍，并保持長(zhǎng)高比r不變，以實(shí)現(xiàn)網(wǎng)格尺寸的平滑過(guò)渡。

根據(jù)線性波的理論解，波面以下z=?0.5λ處質(zhì)點(diǎn)運(yùn)動(dòng)速度小于波面處的5%，z=?λ處質(zhì)點(diǎn)運(yùn)動(dòng)速度小于波面處的0.5%。因此，加密區(qū)④和⑤的高度h4，h5分別為0.5λ和λ，同樣保持網(wǎng)格長(zhǎng)高比r不變，網(wǎng)格尺寸分別為加密區(qū)①的8倍和16倍。

3.3 波高監(jiān)測(cè)設(shè)置

所有算例模擬計(jì)算了各波形從物理時(shí)間t=0到t=40T的傳播過(guò)程。為監(jiān)測(cè)波幅及周期在時(shí)間和空間上的變化，在x=1,2,···,8λ的位置設(shè)置了8個(gè)探針。對(duì)于x>8λ的區(qū)域，因?qū)儆谙▍^(qū)，故不再進(jìn)行監(jiān)測(cè)。

4 各參數(shù)對(duì)造波精度的影響

根據(jù)上述分析，針對(duì)可能影響造波精度的各重要參數(shù)進(jìn)行了研究，參數(shù)包括：影響迭代誤差的內(nèi)迭代次數(shù)n以及影響離散誤差的幾何分辨率和庫(kù)朗數(shù)σ。其中，幾何分辨率由加密區(qū)①波高方向的網(wǎng)格數(shù)Nz和網(wǎng)格長(zhǎng)高比r來(lái)控制，庫(kù)朗數(shù)σ則通過(guò)改變時(shí)間步長(zhǎng)來(lái)控制。誤差主要關(guān)注x=8λ位置的探針在t=40T時(shí)刻的波幅誤差em和相位偏差ep，即每組算例中em和ep的最大值。作為對(duì)照算例A0，其參數(shù)設(shè)置如表1所示，其他算例則在A0算例的基礎(chǔ)上修改相關(guān)參數(shù)得到。

表1 對(duì)照算例A0的參數(shù)設(shè)置Table 1 Parameter setting of Case A0 for comparison

4.1 內(nèi)迭代次數(shù)

為研究?jī)?nèi)迭代次數(shù)n對(duì)迭代誤差的影響，在對(duì)照算例A0的基礎(chǔ)上修改了4個(gè)不同的內(nèi)迭代次數(shù)。分別設(shè)為5,8,10和20次，對(duì)應(yīng)算例編號(hào)為B1～B4，如表2所示。表中，加黑數(shù)字表示本組算例中的變量參數(shù)（以下同）。誤差計(jì)算結(jié)果如圖3所示。

表2 B1～B4算例參數(shù)設(shè)置Table 2 Parameter setting for Case B1-B4

圖3 不同內(nèi)迭代次數(shù)對(duì)應(yīng)的誤差變化Fig.3 Error variation corresponding to different inner iteration times

由圖3計(jì)算結(jié)果可知，波幅誤差em和相位偏差ep隨內(nèi)迭代次數(shù)n的增加而降低；當(dāng)內(nèi)迭代次數(shù)達(dá)到15次以上時(shí)，波幅誤差em降至3.18%，相位偏差ep降至0.17 rad，這兩個(gè)誤差已達(dá)到較低水平，且隨著內(nèi)迭代步數(shù)的繼續(xù)增大，并無(wú)明顯改善，因此，選擇15次作為以下算例的內(nèi)迭代次數(shù)。

4.2 庫(kù)朗數(shù)

在對(duì)照算例A0的基礎(chǔ)上修改庫(kù)朗數(shù)取值，并研究其對(duì)造波精度的影響。設(shè)置從0.2～0.5的4個(gè)不同庫(kù)朗數(shù)取值，算例編號(hào)對(duì)應(yīng)C1～C4，如表3所示。誤差計(jì)算結(jié)果如圖4所示。

表3 C1～C4算例參數(shù)設(shè)置Table 3 Parameter setting for Case C1-C4

圖4 不同庫(kù)朗數(shù)對(duì)應(yīng)的誤差變化Fig.4 Error variation corresponding to different Courant numbers

由圖4計(jì)算結(jié)果可知，波幅誤差em和相位偏差ep都隨庫(kù)朗數(shù)的增加而急劇增大。在庫(kù)朗數(shù)為0.5時(shí)，數(shù)值波波幅與理論值相比衰減了33%，相位偏差可達(dá)2.31 rad，對(duì)離散誤差有極大影響。當(dāng)波面處的庫(kù)朗數(shù)為0.1時(shí)，波幅誤差em可控制在5%以內(nèi)，相位偏差ep約為0.2 rad。繼續(xù)調(diào)低庫(kù)朗數(shù)雖可以進(jìn)一步降低誤差，但卻會(huì)使時(shí)間步過(guò)小，計(jì)算成本將成倍增加。鑒此，綜合考慮計(jì)算精度和成本后，本文使用0.1作為后續(xù)算例的庫(kù)朗數(shù)較為合理。

4.3 波高方向網(wǎng)格數(shù)

波高方向網(wǎng)格數(shù)Nz影響沿波高方向網(wǎng)格的幾何分辨率，進(jìn)而對(duì)離散誤差造成影響。設(shè)置4個(gè)不同的Nz值算例進(jìn)行誤差影響的驗(yàn)證。具體參數(shù)設(shè)置和誤差結(jié)果見表4和圖5。

表4 波高方向網(wǎng)格數(shù)影響算例參數(shù)設(shè)置Table 4 Parameter setting of cases with different mesh number in the wave height

圖5 波高方向網(wǎng)格數(shù)對(duì)應(yīng)的誤差變化Fig.5 Error variation corresponding to the mesh number in the wave height direction

由表4和圖5所示計(jì)算結(jié)果可知，隨著波高方向網(wǎng)格數(shù)Nz的增加，波幅誤差em和相位偏差ep都呈明顯的下降趨勢(shì)。當(dāng)Nz為20時(shí)，波幅誤差em降至3.5%以下，相位偏差ep降至0.2 rad以下，造波精度在可以接受的范圍內(nèi)。算例D3和D4相較于A0，雖然兩種誤差都有一定的降低，但過(guò)大的網(wǎng)格數(shù)會(huì)造成計(jì)算成本的急劇增大。鑒此，綜合考慮計(jì)算精度及計(jì)算成本，后續(xù)計(jì)算選擇20作為波高方向的網(wǎng)格數(shù)值。

4.4 網(wǎng)格長(zhǎng)高比

波高方向網(wǎng)格數(shù)Nz確定后，網(wǎng)格長(zhǎng)高比r決定了波傳播方向的網(wǎng)格分辨率，這是影響離散誤差的另一個(gè)重要因素。在試算中發(fā)現(xiàn)，網(wǎng)格長(zhǎng)高比r對(duì)誤差的影響在不同波陡下有很大差異。為研究其具體影響，在保持波長(zhǎng)不變的前提下，進(jìn)行了0.01～0.07共7組不同波陡下的數(shù)值計(jì)算，每組波陡下完成3～4個(gè)不同長(zhǎng)高比網(wǎng)格的算例，共計(jì)26個(gè)。各算例的具體參數(shù)設(shè)置和誤差結(jié)果如表5和圖6所示。

表5 不同網(wǎng)格長(zhǎng)高比算例參數(shù)設(shè)置Table 5 Parameter setting for cases with different mesh aspect ratios

波幅誤差em和相位偏差ep的對(duì)比分別在t=20T和t=40T這2個(gè)時(shí)刻及x=8λ位置處展開。由誤差結(jié)果分析可知，隨著網(wǎng)格長(zhǎng)高比r的增加，網(wǎng)格分辨率的降低使得em和ep都逐漸增大，且變化趨勢(shì)基本相同；不同波陡的波浪對(duì)于網(wǎng)格長(zhǎng)高比r的敏感性不同，波陡較小的算例（例如，E1和E2兩組）在使用較大的r來(lái)劃分網(wǎng)格時(shí)仍能保持較低的誤差，而對(duì)于大波陡的算例（例如，E6和E7兩組），則需要非常小的r才能滿足誤差要求。鑒此，綜合考慮計(jì)算精度及計(jì)算成本，對(duì)于波陡在0.01～0.06之間的波浪，本文通過(guò)對(duì)圖6及圖7所示的結(jié)果進(jìn)行擬合可得到網(wǎng)格長(zhǎng)高比，具體由下式計(jì)算：

圖6 em與ep隨網(wǎng)格長(zhǎng)高比的變化（t=20T，x=8λ）Fig.6 Variation of em and ep with grid aspect ratio at t=20T and x=8λ

圖7 em與ep隨網(wǎng)格長(zhǎng)高比的變化（t= 40T，x=8λ）Fig.7 Variation of em and ep with grid aspect ratio at t= 40T and x=8λ

根據(jù)式（15）計(jì)算的r可以滿足：t=20T時(shí)刻的波幅誤差em約5%，相位偏差ep約0.2 rad；t=40T時(shí)刻波幅誤差em約10%, 相位偏差ep約0.4 rad。

而對(duì)于波陡H/λ≥0.06的波浪而言，即使使用網(wǎng)格長(zhǎng)高比為1的網(wǎng)格進(jìn)行計(jì)算，仍然不能滿足上述誤差要求。對(duì)于此類大波陡波浪，可能需要改進(jìn)自由液面捕捉方法或者使用其他創(chuàng)新方法等才能夠得到較為理想的結(jié)果。

4.5 參數(shù)設(shè)置的總結(jié)

經(jīng)過(guò)上述數(shù)值計(jì)算及結(jié)果對(duì)比可知，為實(shí)現(xiàn)多周期穩(wěn)定數(shù)值造波，推薦采用以下設(shè)置：

1） 選用標(biāo)準(zhǔn)低雷諾數(shù)K-Epsilon湍流模型，湍動(dòng)能和湍流耗散率的初始值分別為1×10?5J/kg和1×10?4m2/s3，以降低湍流黏度對(duì)波面衰減的影響；

2） 在計(jì)算域下游邊界設(shè)置阻尼消波區(qū)，阻尼消波區(qū)的長(zhǎng)度設(shè)為兩倍波長(zhǎng)（xd=2λ）；僅使用線性消波項(xiàng)（f2=0），線性阻尼常數(shù)f1=Ψ1ω；

3） 使用HRIC方法，并激活梯度平滑選項(xiàng)；

4） 內(nèi)迭代次數(shù)設(shè)為15次或以上，自由液面處最大庫(kù)朗數(shù)維持在0.1以下，并據(jù)此調(diào)整時(shí)間步長(zhǎng)；

5） 對(duì)自由液面附近進(jìn)行網(wǎng)格加密，加密區(qū)①波高方向的網(wǎng)格數(shù)Nz設(shè)為20，網(wǎng)格長(zhǎng)高比r使用式(15)計(jì)算，其他加密區(qū)網(wǎng)格尺寸參考4.3節(jié)。

計(jì)算結(jié)果表明，使用上述參數(shù)設(shè)置方法，可以實(shí)現(xiàn)多周期的穩(wěn)定數(shù)值造波，波幅誤差em和相位偏差ep滿足工程應(yīng)用的精度要求。圖8所示為算例E4.4在x=8λ位置波高隨時(shí)間的變化曲線，圖9所示為t=40T時(shí)刻波高隨空間位置的變化曲線，其中波高、時(shí)刻和位置都進(jìn)行了無(wú)量綱化處理。由圖8和圖9可以看出，實(shí)際波幅線（黑色）與理論波幅線（紅色）基本重疊，實(shí)際波幅與理論波幅在時(shí)間上吻合較好，空間上波形分布均勻，但隨著傳播距離的增加存在一定的衰減，可見，消波區(qū)的波形消波效果良好。

圖8 算例E4.4在x=8λ位置波高隨時(shí)間變化Fig.8 Variation of wave height with time at x=8λ in Case E4.4

圖9 算例E4.4在t=40T時(shí)刻波高隨空間位置的變化Fig.9 Variation of wave height with spatial position at t=40T in Case E4.4

5 數(shù)值造波方法應(yīng)用

以上算例都是針對(duì)模型尺度下的二維數(shù)值造波研究，為驗(yàn)證本文所提參數(shù)設(shè)計(jì)方法在三維數(shù)值造波問(wèn)題上的適用性，計(jì)算了KCS船型在迎浪中的阻力和運(yùn)動(dòng)響應(yīng)，并與工程流體動(dòng)力學(xué)（EFD）試驗(yàn)值[10]進(jìn)行對(duì)比。

為與實(shí)驗(yàn)條件保持一致，波長(zhǎng)設(shè)為λ=0.5LBP=115 m，波陡H/λ=1/60，弗勞德數(shù)Fr=0.26。由于船體左右對(duì)稱，為減少計(jì)算量，僅對(duì)半船的流場(chǎng)進(jìn)行模擬，并使用重疊網(wǎng)格方法及DFBI (displaying dynamic fluid body interaction)模型實(shí)現(xiàn)船體迎浪時(shí)的升沉與縱傾運(yùn)動(dòng)。

背景區(qū)域及重疊域網(wǎng)格采用切割體非結(jié)構(gòu)網(wǎng)格形式，其中，背景區(qū)域網(wǎng)格基準(zhǔn)尺寸取為4.0H，重疊區(qū)域網(wǎng)格基準(zhǔn)尺寸取為0.2H，并在背景區(qū)域設(shè)置多層嵌套的加密區(qū)以保證體網(wǎng)格尺寸的過(guò)渡平滑；船體邊界層采用棱柱層網(wǎng)格來(lái)捕捉，邊界層共劃分有16層，增長(zhǎng)率為1.2。自由液面處則根據(jù)4.3節(jié)所述參數(shù)進(jìn)行加密，y軸方向網(wǎng)格大小取為4.0H。最終，得到的背景區(qū)域體網(wǎng)格總數(shù)為430萬(wàn)，重疊區(qū)域體網(wǎng)格總數(shù)為1 100萬(wàn)，圖10所示為中縱剖面船體附近的網(wǎng)格示意圖。

圖10 船體附近網(wǎng)格Fig.10 Grids near the hull

總共計(jì)算20個(gè)波周期，在中縱剖面船體前方一倍LBP處設(shè)置探針，監(jiān)測(cè)得到的波高變化曲線如圖11所示，由圖11可知，在整個(gè)計(jì)算時(shí)間內(nèi)來(lái)波的波幅都無(wú)明顯衰減。計(jì)算結(jié)果主要關(guān)注波浪增阻系數(shù)σaw、升沉和縱傾的傳遞函數(shù)（transfer function，TF），其計(jì)算公式分別如下：

圖11 船前一倍LBP處波高的變化Fig.11 Variation of wave height at a distance of onefold-LBP ahead the ship

式中：Fx,wave和Fx,calm為船體在波浪和靜水中的縱向受力的平均值；BWL為水線面寬度；LBP為垂線間長(zhǎng)；x31，x51和ζI1分別為升沉、縱搖和波高瞬時(shí)幅值的一階傅里葉系數(shù)；波數(shù)k=2π/λ。其中，升沉和縱搖以船體重心為參考點(diǎn)。

圖12所示為自由液面在初始和計(jì)算一段時(shí)間后的自由液面波形對(duì)比圖，由自由液面波高圖可清晰地捕捉到開爾文船行波。增阻系數(shù)、升沉和縱傾的傳遞函數(shù)的計(jì)算結(jié)果見表6。

圖12 自由液面波形對(duì)比Fig.12 Comparison of wave pattern for free surface

表6 CFD與EFD計(jì)算結(jié)果對(duì)比Table 6 Comparison of calculation results by CFD and EFD

由表6的的計(jì)算結(jié)果可知，CFD計(jì)算結(jié)果與EFD結(jié)果相比均偏低，誤差約為4%～10%，這進(jìn)一步驗(yàn)證了本文所提造波參數(shù)設(shè)計(jì)方法可以推廣到三維數(shù)值造波的相關(guān)問(wèn)題之中。

6 結(jié) 論

本文分析了一維波動(dòng)方程的離散形式，得到了影響數(shù)值造波精度的幾個(gè)關(guān)鍵參數(shù)及其對(duì)造波效果的具體影響。通過(guò)設(shè)置不同參數(shù)對(duì)深水二維水池進(jìn)行了波浪模擬，對(duì)每組算例的計(jì)算誤差進(jìn)行評(píng)估，綜合考慮計(jì)算精度和成本，確定了各參數(shù)的最佳值。而且，提出了包括湍流模型、消波方法、內(nèi)迭代次數(shù)、庫(kù)朗數(shù)及網(wǎng)格設(shè)置的一整套參數(shù)設(shè)計(jì)方法。通過(guò)分析計(jì)算，得到如下結(jié)論：

1） 數(shù)值造波誤差主要包括離散誤差和迭代誤差。離散誤差是影響數(shù)值造波精度的主要誤差因素，與之相比的迭代誤差可通過(guò)增加內(nèi)迭代次數(shù)來(lái)降低到遠(yuǎn)小于離散誤差的范圍；

2） 影響離散誤差大小的主要因素為庫(kù)朗數(shù)和網(wǎng)格分辨率。其中，庫(kù)朗數(shù)對(duì)誤差影響較大，需確保自由液面庫(kù)朗數(shù)在0.1以下，網(wǎng)格分辨率則需通過(guò)調(diào)整波高方向網(wǎng)格數(shù)及網(wǎng)格長(zhǎng)高比來(lái)實(shí)現(xiàn)；

3） 波陡H/λ對(duì)數(shù)值造波的結(jié)果影響很大。在H/λ≤0.06的波陡范圍內(nèi)，使用4.5節(jié)參數(shù)設(shè)計(jì)方法可使得20個(gè)波周期和40個(gè)波周期內(nèi)波幅誤差分別約為5%和10%，相位偏差分別約為0.2和0.4 rad，并可推廣到三維造波問(wèn)題中；

4） 對(duì)于五階Stokes波以外的其他規(guī)則的非線性波或線性波，計(jì)算結(jié)果表明，上述參數(shù)設(shè)計(jì)方法均可一定程度上減少造波誤差，使造波結(jié)果滿足精度要求；

5） 對(duì)于波陡H/λ≥0.06的大波陡波浪，即使使用網(wǎng)格長(zhǎng)高比r=1的網(wǎng)格計(jì)算仍有較大誤差，若為提高計(jì)算精度而進(jìn)一步增加網(wǎng)格分辨率，則計(jì)算成本將變得過(guò)高，不具有實(shí)際應(yīng)用的可行性。對(duì)于此類大波陡數(shù)值造波問(wèn)題，仍需通過(guò)進(jìn)一步研究和創(chuàng)新方法來(lái)解決。