呂磊，陳作鋼

1 上海交通大學 海洋工程國家重點實驗室，上海 200240

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

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

0 引 言

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

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

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

1 數(shù)學模型

1.1 控制方程

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

1.2 湍流模型

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

1.3 自由液面捕捉

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

式中，qai為體積分數(shù)ai的一個體積源。

1.4 造 波

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

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

1.5 消 波

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

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

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

1.6 離散格式及其他設置

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

2 誤差分析

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

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

若以顯式一階迎風格式對波面方程(9)進行有限差分，離散誤差則為式(10)等號右端項之和。

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

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

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

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

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

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

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

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

3 算例設置

3.1 計算域設置

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

圖1 計算域設置示意圖Fig.1 Schematic view of computational domain setting

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

3.2 網(wǎng)格設置

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.2 庫朗數(shù)

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

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

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

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

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

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

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

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

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

4.4 網(wǎng)格長高比

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

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

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

圖6 em與ep隨網(wǎ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)格長高比的變化（t= 40T，x=8λ）Fig.7 Variation of em and ep with grid aspect ratio at t= 40T and x=8λ

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

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

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

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

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

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

3） 使用HRIC方法，并激活梯度平滑選項；

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

6 結(jié) 論

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

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

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

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

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

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