蘇紹娟，王有志，王天霖，李新飛

(1. 大連海事大學 船舶與海洋工程學院，遼寧 大連 116026；2. 哈爾濱工程大學 船舶與海洋工程學院，黑龍江 哈爾濱 150001)

0 引 言

隨著計算機能力以及對圖形、數(shù)字處理能力的提高，數(shù)值水槽模型成為國內外重點研究對象之一，與實驗相比具有成本低，效率高的優(yōu)勢，且能避免縮比效應帶來的影響。目前數(shù)值仿真軟件的拓展模塊越來越豐富，功能越來越全，進一步提高了計算能力。當波浪與結構物的相互作用時，波面會出現(xiàn)破碎、攀爬以及非線性響等現(xiàn)象，同時周圍流場會發(fā)生變化，與大型船舶和海洋平臺相比，其結果對小型結構物將產(chǎn)生更大的影響，所以流體的粘性一般不能忽視。由于浮標類結構物尺寸較小，主要受影響的波浪其波高波長較小，周期較短，并不需要建造較大的數(shù)值水槽進行數(shù)值模擬，所以如何保證小型數(shù)值水槽的波浪傳遞的穩(wěn)定成為計算結果準確性的關鍵。

Brorsen M和Larsen J[1]根據(jù)源造波理論，通過在計算域內設置的造波源可分別向源區(qū)域兩邊同時生成2列方向相反的波；劉加海[2]通過設置邊界造波版的運動來模擬物理造波板的運動過程，實現(xiàn)了二維規(guī)則波的數(shù)值模擬，并進一步分析了模擬波高隨造波板運動周期、振幅和坐標的關系；Lal和Elangovan[3]通過設置不同的波浪周期等參數(shù)分析波高波長的影響因素，并研究了基于阻尼消波的沙灘斜坡方法的影響，結果表明當沙灘斜坡比例為1∶3時消波效果最佳；黃華和鄧冰[4]通過建立二維數(shù)值水槽對2階stokes波進行仿真模擬，運動搖板造波和增設消波區(qū)域保證了短時間內水槽內波浪的穩(wěn)定性；Kim[5]通過調整消波區(qū)域的網(wǎng)格尺寸使得數(shù)值計算衰減來進行波浪反射的消除，通過改變網(wǎng)格尺寸的大小找到最佳的求解設置；Tian[6]等人為研究圓柱體與波浪的相互作用建造數(shù)值水槽，使用入口速度法進行數(shù)值造波，并通過尾部設置阻尼進行消波，得到了較為精確的規(guī)則波，并能保證一段時間內的穩(wěn)定性；Fabio M和 Marques Machado[7]則對不同邊界速度法和推板造波法進行了波浪模擬對比，結果表明邊界速度法受波浪阻尼的影響更大，而推板造波法能造出更加準確地波浪。

本文使用STARCCM+軟件，分別基于傳統(tǒng)純粘流以及勢粘流結合的方法建立數(shù)值波浪水槽，使用邊界流體函數(shù)直接輸入的方式進行造波，通過對比波浪模擬結果，分析2種數(shù)值波浪水槽的特點，為數(shù)值波浪水槽的建立提供科學的參考依據(jù)。

1 數(shù)值方法

1.1 流體體積法

在數(shù)值造波時需要跟蹤氣體和液體的界面，流體體積法(VOF)是由美國學者Hirt和Nichols等[8]提出的。在流場中的每個網(wǎng)格，這個函數(shù)定義為目標流體的體積與網(wǎng)格體積的比值，只要知道這個函數(shù)在每個網(wǎng)格上的值，就可以實現(xiàn)對運動界面的追蹤。因此在這里選用VOF方法處理氣體和液體的界面跟蹤問題，該方法跟蹤界面的連續(xù)性方程為：

在每個控制體積內，所有相的體積分數(shù)總和為1，所有變量及其屬性在控制容積內各項共享，表示第q項流體在單元中的體積分數(shù)，在數(shù)值水槽模擬中，因為只有氣液兩相介質，q取2；造波時將VOF模型相數(shù)設置為2，第一相為空氣，第二相為水。

1.2 基本控制方程

目前雷諾平均方程（RANS）方法是計算流體力學中最常用的方法之一，從N-S方程出發(fā)對湍流進行直接的數(shù)值模擬（DNS）難以解決工程中遇到的復雜湍流問題，依靠實驗取得經(jīng)驗數(shù)據(jù)，不僅耗資巨大，周期很長，而且對于某些實際工程問題，完全相似的實驗室模擬不可能實現(xiàn)。在這種情況下，求解雷諾平均的N-S方程（RANS）方法成為解決工程問題比較有效和切實可行的手段，其中空間離散采用有限體積法，采用二階時間差分方法保證計算的準確性。其連續(xù)性方程和動量方程為：

式中：ui表示略去平均符號的雷諾平均速度分量，ρ為密 度 ，P為 壓 強 ，為脈 動 速 度，為應力張量分 量 。

1.3 波浪水槽參數(shù)

目前常見的數(shù)值造波方法主要可以分為兩大類：源造波法和邊界造波法。源造波法又分為質量源造波法和動量源造波法[9]；邊界造波法主要包括物理造波法和邊界流體速度函數(shù)直接輸入法；常用的消波方法包括常用的方法有輻射邊界條件，主動消波器，阻尼消波區(qū)[10]；波浪水槽邊界條件的設置則要保證至少一個的速度入口和壓力出口，前后邊界選用對稱平面，頂部距水面要足夠高，以防止空氣流動對波高的影響，底部設為壁面，水深設置為6 m。其他邊界條件及具體設置如下：

1）模擬1

選用邊界流體函數(shù)直接輸入進行波浪的生成，在出口處增加阻尼區(qū)域進行消波，全流域為粘性流體；因為出口處設有阻尼區(qū)域，如果將出口邊界定義為速度入口將定義出口區(qū)域的波高，這與消波意愿相違背，所以將出口邊界設置為唯一壓力出口，頂部設置為速度入口，實施開放邊界允許空氣雙向流過邊界。具體如圖1（a）所示。

2）模擬2

選用邊界流體函數(shù)直接輸入進行波浪的生成，在出入口兩端設置一定長度的勢流區(qū)域防止波的反射，整個流域中間為粘流區(qū)域，兩邊為勢流區(qū)域；因為不需要考慮出口處反射波的影響，出口邊界設置為速度入口定義出口區(qū)域波高，進而更加有利于保持中間粘流區(qū)域的波浪精度，頂部設置為唯一壓力出口，保證壓力的釋放，具體如圖1（b）所示。

文中模擬的波浪參數(shù)如表1所示。

1.4 網(wǎng)格劃分

圖 1 數(shù)值波浪水槽示意圖Fig. 1 Numerical wave tank schematic diagram

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

通過切割體網(wǎng)格生成器生成六面體網(wǎng)格，為了更好地捕捉水線面，需要對水線面附近的網(wǎng)格進行多層加密以減少壓力傳遞過程中的耗散，全域網(wǎng)格劃分如圖2（a）所示，水線面網(wǎng)格加密情況如圖2（b）所示。

圖 2 網(wǎng)格劃分情況Fig. 2 Meshing situation

表2為Havn[11]提出的網(wǎng)格劃分的具體標準，為了具體研究網(wǎng)格對波浪傳遞的影響，劃分了3套質量不同的網(wǎng)格，其中第2和第3套符合Havn提出的網(wǎng)格劃分要求，而第1套低于標準，具體網(wǎng)格質量見表3。

在其他一些參數(shù)設置中，選取k-ε湍流模型進一步還原真實環(huán)境。對于時間步長的選擇標準并不統(tǒng)一，大的時間步長影響計算的精度，而小的時間步長則影響計算的總體速度；Finnegan和Goggins給出最佳的時間步長為波周期的1/50（T/50），Havn[11]指出時間步長值小于或等于T/100，通過對不同時間步長的收斂性驗證，根據(jù)庫朗數(shù)值及殘差最終選取的時間步長間隔為0.004 s，約為T/500。

表 2 網(wǎng)格劃分標準Tab. 2 Meshing standards

表 3 網(wǎng)格參數(shù)Tab. 3 Mesh parameters

2 結果及結論

在T=0 s時刻初始化定義全域波面，對不同時刻的全域波面以及各個監(jiān)測位置波高進行監(jiān)測。

2.1 網(wǎng)格和時間步長的選取

對于網(wǎng)格的選取，通過基于純粘流數(shù)值水槽對不同質量的網(wǎng)格進行評價，3套網(wǎng)格的同一時刻波面高度情況如圖3所示。

圖 3 不同網(wǎng)格波面高度Fig. 3 Wave front heights of different mesh

通過對結果的分析表明，增加網(wǎng)格尺寸會加劇導致波浪的衰減，第1套網(wǎng)格相較于第2和第3套網(wǎng)格數(shù)值衰減更明顯，且第2和第3套網(wǎng)格模擬下的波浪結果相差不大，綜合考慮計算準確度與計算速度，認定網(wǎng)格2為最優(yōu)網(wǎng)格劃分方案，其基本劃分尺寸符合Havn提出的網(wǎng)格劃分標準。

對于時間步長的選取，使用第2套網(wǎng)格并基于勢粘流結合理論的數(shù)值水槽進行驗證，選取L3監(jiān)測點對不同時間步長下的波面高度結果進行對比，結果如圖4所示。

圖 4 不同時間步長L3位置對應波面高度Fig. 4 Wave front height at L3 position at different time steps

通過對結果的分析表明，時間步長為0.005 s（T/400）時計算結果與0.004 s（T/500）和0.003 s（T/650）有著很明顯的誤差，且0.004 s和0.003 s兩種時間步長下的結果十分接近。圖5所示計算過程中0.004 s時間步長下的庫朗數(shù)穩(wěn)定在0.19左右，而0.005 s時間步長下的庫朗數(shù)穩(wěn)定在0.28左右，而這2個時間步長都遠遠小于T/100，所以時間步長不能僅根據(jù)波浪周期的大小來選取，還要參考庫朗數(shù)和網(wǎng)格質量來綜合確定。

圖 5 全域庫朗數(shù)標量場景Fig. 5 Scalar scene of global courant number

2.2 波浪模擬結果

為了更直觀地觀察和對比2種數(shù)值水槽造波的結果，分別將第20 s，40 s，60 s的波面情況與理論值進行對比，如圖6所示。

結果表明，基于純粘流的水質波浪水槽，沿著波浪傳遞的方向波高逐漸減小，隨著模擬時間的增加現(xiàn)象也越來越明顯。這是因為波浪傳遞過程中由于流體存在粘性，部分動能轉化為內能耗散，如圖7所示。通過增加阻尼區(qū)域防止波浪在出口處的反射，但增設的阻尼對上游波浪的傳遞也會造成影響，且從第一波到達阻尼區(qū)域開始由近至遠對波浪造成影響，離阻尼區(qū)域越近反射到來的越快，波面受其影響也會越大，且這種影響伴隨著整個計算過程將越來越明顯。增加計算域的長度可以一定程度減小阻尼區(qū)域對計算區(qū)域的影響，但由此導致網(wǎng)格的增加會大大增加模擬時間。也可以通過減少模擬的時間來確保一段時間波浪的穩(wěn)定性，但是這只適合較短時間的數(shù)值模擬，對于較長時間的模擬并不可行。

圖 6 波面數(shù)值解與理論值比較Fig. 6 Comparison of wave front numerical solutions and theoretical values

圖 7 基于純粘流數(shù)值水槽不同位置波面高度隨時間變化情況Fig. 7 Viscous flow numerical water tank wave surface height changes with time

基于勢粘流結合建造的數(shù)值水槽的波浪與理論值相比更為接近，且波高不隨著時間而衰減，因為波浪在非粘流區(qū)域傳遞并不會造成動能轉化為內能的損耗，出口處的勢流區(qū)域保證了波浪不受出口反射波的影響，波浪的精度只會因數(shù)值計算的耗散以及中間粘流區(qū)域將少量動能轉化為內能而影響，相比較而言具有更高的精度。

圖 8 兩種數(shù)值水槽L3位置波面隨時間變化情況Fig. 8 Wave front changes with time of two numerical wave tanks

3 結 語

為了探究規(guī)則波在傳遞過程中的衰減問題，運用STARCCM+軟件基于不同方式建立了2種數(shù)值波浪水槽，通過收斂性驗證確定了最優(yōu)的網(wǎng)格大小以及時間步長，結果顯示造波網(wǎng)格可按照Havn提出的標準來劃分，波高方向網(wǎng)格最小尺寸不小于H/25，波長方向網(wǎng)格最小尺寸不小于L/100；時間步長的選取應兼顧庫朗數(shù)以及波浪周期，庫朗數(shù)值應小于0.2，時間步長參考值不小于T/500。

基于純粘流的數(shù)值水槽因為受流體粘性的影響，會造成波浪傳遞過程中部分動能的耗散，由于初始化定義了T=0 s的波面，模擬從第一時間就會受到消波阻尼反射的影響，并從下游逐漸影響至上游，且影響越來越大，導致不能進行長時間的計算模擬。為了保證計算區(qū)域波浪的精度，可以通過加長數(shù)值水槽使阻尼反射影響更晚到達計算區(qū)域，但這必然會增加網(wǎng)格數(shù)目進而導致計算速度的減慢。

對于基于勢粘流結合的數(shù)值水槽，其中間計算區(qū)域流體為粘流，兩側流體為勢流，使得上游波浪并不會受到出口處反射波的影響，且波浪在勢流中的傳遞也不會造成動能的損失，保證了中間粘流區(qū)域波浪精度的同時進而提高了計算結果的準確性。

綜合上述分析，傳統(tǒng)的基于純粘流的數(shù)值水槽適用于時間較短或者水槽長度較長的模擬，且波浪傳遞過程中衰減明顯?；趧菡沉鹘Y合的新型數(shù)值水槽傳遞過程中并無明顯衰減，滿足長時間計算模擬的要求。