船舶在長峰不規則波中頂浪縱向運動的數值模擬

劉 可，吳 明，楊 波

(1.中國人民解放軍91550部隊，遼寧 大連 116026;2.海軍大連艦艇學院，遼寧 大連 116018)

0 引言

船舶航行的縱向運動主要包括船舶的縱搖及垂蕩。嚴重的縱搖和垂蕩將引起甲板上浪、砰擊和失速等一系列后果，對船舶在海上的航行安全以及艦載武器的使用影響極大。因此，本文選取船舶的縱向運動作為研究對象，基于CFD方法對船模在長峰不規則波中頂浪運動進行數值模擬，具有迫切的現實意義。

1 波浪的數值模擬

要實現對船模在長峰不規則波中頂浪縱向運動的數值模擬，首要前提就是構建精度滿足耐波性計算要求的長峰不規則波數值波浪水池。所謂的數值波浪水池，是指對非線性波浪水動力以及浮體運動數值模擬設施的統稱。它能夠通過實驗觀測為各種波浪理論的研究奠定堅實的物理基礎，還能為海洋、船舶等工程設計提供可靠而高效的試驗數據。本文基于粘性流理論構建長峰波數值波浪水池，采用有限體積法 (FVM)對RANS方程和連續控制方程進行離散求解，利用Fluent軟件的二次開發功能UDF，編寫邊界條件完成數值造波和阻尼消波功能。選用ITTC海浪譜作為目標海浪譜，構建長峰波數值波浪水池，并對其計算精度進行誤差計算。

1.1 數值波浪水池的關鍵技術

1.1.1 數值造波

數值造波是指用數值方法模擬波浪的生成過程，為數值模型實驗提供各種形式的波浪環境條件。本文采用的數值造波方法是邊界條件造波法，即在數值波浪水池［1－2］入口邊界設置3個方向的速度。X方向沿水池向下游為正，Z方向向上為正，Y軸與X軸和Z軸符合右手法則。3個方向的入射速度分別為:

式中:η為波動水面相對于靜止水面的瞬時高度;Ai，ki，ωi和εi分別為第i個組成波的波幅、波數、圓頻率和初始相位，εi是在(0，2π)范圍內的隨機相位;X軸為波浪傳播方向;U，V，W分別為波浪X軸、Y軸和Z軸的速度分量。

1.1.2 數值消波

當波浪到達水池末端開邊界處，會引發水波的二次反射，反射波與入射波疊加，會導致波場的失真。因此，消除反射波的影響也是數值波浪水池的一項重要技術。數值波浪水池常用的消波技術主要有輻射邊界條件消波、主動消波器消波、阻尼消波3種。本文選取的是阻尼消波。

阻尼消波法是指在流場中添加人工粘性，因其對來波的頻率和波長不敏感，可以有效地消除各種頻率和波長的來波，因而被廣泛采用。在計算域出口邊界前設置1～2倍波長的阻尼消波段，利用Fluent中的 UDF宏 DEFINE_SOURCE(mom_source，cell，thread，dS，eqn)編程實現消波。

在阻尼消波段內，動量方程寫為:

其中μ(x)為在阻尼段起點為0的單調遞增函數，可以取為線性遞增、指數遞增等形式。

式中:Xmin_D和Xmax_D分別為消波區的最小、最大X坐標。

1.2 長峰不規則波數值波浪水池的構建

1.2.1 網格劃分

三維數值水槽的網格劃分如圖1所示，本文構建的長峰不規則波數值水池［3－5］長18 m，其中12～18 m為消波區，寬3 m，深2.5 m，自由面以上1.5 m，整個水槽的網格數為413820。垂直自由面方向的網格尺寸取為有義波高1/5。造波區沿X軸正方向網格尺寸與垂直于自由面Z方向的最小網格尺寸相同，從自由面到水槽頂部網格按1∶1.1的比例等比分布;從自由面到水池底部網格按1∶1.05的比列等比分布。網格基本上是離自由面越遠尺寸越大。對于消波區的網格劃分，垂直方向劃分與造波區的一致，水平方向網格以造波區網格尺度為基準向右邊界逐漸擴大。

圖1 數值波浪水池示意圖Fig.1 Sketch map of the numerical wave tank

1.2.2 算例描述

選取ITTC海浪譜作為目標靶譜，對3種海況下的長峰不規則波進行數值模擬。波浪的目標參數見表1。

表1 ITTC海浪譜數值模擬參數Tab.1 Parameters for the simulation of ITTC wave spectrum

1.2.3 數值模擬結果

長峰不規則波數值模擬瞬時波面場和局部速度矢量如圖2和圖3所示。

圖2 瞬時波面 (t=87.36 s)Fig.2 Instantaneous wave line

圖3 局部速度矢量圖Fig.3 Vectorgraph of partial speed

圖4和圖5是有義波高分別取H1/3=0.08 m，H1/3=0.125 m時，長峰不規則波數值波浪水池X=1 m，X=6 m，X=11 m處的波面時歷曲線對比圖。

圖4 時歷曲線 (H1/3=0.08 m，Ts=1.092 s)Fig.4 Time history of wave line

圖5 時歷曲線 (H1/3=0.125 m，Ts=1.365 s)Fig.5 Time history of wave line

從上述數值造波水池不同位置的波面時歷曲線對比可看出，長峰不規則波在沿X軸正方向向下游傳播過程中，波能有一定程度的衰減，這是由于波浪在傳播過程中，高頻子波衰減導致的。

1.2.4 數值模擬精度誤差計算

對上面監測得到長峰不規則波波面時歷進行譜分析與目標譜對比如圖6所示。

圖6 數值模擬海浪譜Fig.6 Numerical simulation spectrum

分別從譜面積 m0、譜峰頻率 ωp和有義波高H1/33個方面對上述數值模擬海浪譜進行誤差分析，誤差計算結果如表2和表3所示。

表2 長峰不規則波相關誤差計算 (ITTC譜X=1 m處波面時歷)Tab.2 Relative errors for long-crested irregular waves

表3 長峰不規則波相關誤差計算(ITTC譜X=6 m處波面時歷)Tab.3 Relative errors for long-crested irregular waves

2 船模在長峰不規則波中頂浪縱向運動的數值模擬

2.1 船模的選擇

選取具有球鼻首和方位的DTMB 5512船模作為研究對象，該船模是ITTC(國際船模試驗水池會議)推薦的瘦削型標準船模DTMB 5415的全相似幾何模型，以美國海軍DDG－51型驅逐艦為模板，船型數據詳實，而且與我海軍艦船船型相似 (見圖7)，適合軍艦參考。表4為該船模及其所對應的實船尺度的主要數據。

圖7 DTMB 5512船模Fig.7 Ship model of DTMB5512

表4 船模和全尺度船的主要參數Tab.4 The primary parameter of DTMB5512 ship model

2.2 數值模擬方案

2.2.1 試驗計劃

本文對DTMB 5512船模在長峰不規則波中頂浪縱向運動數值模擬［6－8］試驗中，只考慮垂蕩和縱搖2個自由度的運動。船模CFD耐波性數值模擬試驗對計算資源的要求較高，考慮到本文研究所使用的計算機配置的實際情況，以及試驗水池網格劃分所帶來的計算效率等問題，試驗計劃如表5所示。

表5 DTMB 5512型船模頂浪縱向運動數值模擬試驗Tab.5 Numerical simulation test of DTMB 5512 ship model

2.2.2 計算域劃分

參照《水面船模耐波性實驗規程》，將計算域設置成長方體形狀，如圖8所示。船模與計算域各邊界的位置關系如下:入口距船首1倍船長，出口距船尾2倍船長，頂部邊界距水線0.5倍船長，底部邊界距水線1倍船長，左、右邊界距船中縱剖面0.5倍船長。

圖8 水池計算域劃分示意圖Fig.8 Outline of computational domain for wave tank

耐波性數值波浪水池分成5個區域進行網格劃分，即近船體區域、近流場區域、自由面區域、上下遠流場區域及消波區域、各區之間互不重疊，且連接處選擇connected方式，如圖9所示。近船體區網格劃分如圖10所示。

圖9 耐波性波浪數值水池網格劃分Fig.9 Build grid for sea-keeping numerical simulation wave tank

圖10 DTMB 5512型船模近船體網格劃分Fig.10 Build grid for the ship model of DTMB 5512

2.3 船模頂浪縱向運動數值模擬的求解過程

根據船舶六自由度運動的控制方程(5)，當波浪作用于船體時，其運動的速度、角速度以及位置、姿態等可以通過控制方程求解、積分得到。對于流浮耦合運動，波浪作用在船體上的力和力矩使船體產生運動，同時船體的運動又對其周圍流場產生影響。因此本文在數值模擬中，分段計算流體與船體運動的耦合，步驟如下:

1)將船模按初始浮態固定，原點與重心重合;

2)對流場進行初始化，設定初始航速后造波;

3)以時間步長Δt=0.001 s步進;

4)通過當前流場變量迭代求解流場的速度矢量;

5)通過壓力—速度耦合算法獲得壓力場;

6)求解體積分數方程重構自由面;

8)更新船體位置和浮態;

9)返回第3步，求解改變浮態后各量，并根據計算再次改變浮態，按此迭代求解，直到方程組的殘差小于設定值或迭代次數達到設定值;

10)返回第2步并重復以下步驟，直到設定的時間步數計算完畢。

通過以上迭代、循環，可實現流體與船體運動的耦合。

步驟2按照1.2節構建的數值波浪水池進行造波和消波。試驗過程中可以將事先保存好的穩定流場導入耐波性數值模擬水池，進行數值模擬，這樣可以提高計算效率，縮短試驗時間。步驟4～6，按前文設置的數值方法進行計算。步驟7則通過UDF編程實現。步驟9殘差標準可取軟件默認設定值，迭代次數上限為40步。

2.4 數值模擬結果

對于DTMB 5512船模由于帶球鼻首、具有方尾，線性復雜，所以實現其耐波性的數值模擬，要比一般商用船模困難得多。可采用以下方法對實驗進行改進:

1)改進船模貼體網格質量、數量及分區網格的匹配;

2)采用遞增方法實現船模的縱向運動，即先對船模進行單自由度縱搖，穩定之后再增加垂蕩的數值模擬;

3)逐漸增加船模質量 (將船模質量降低到原船模的一半進行數值模擬，當殘差穩定后再逐漸增大質量直至原值);

4)減小欠松弛因子，控制單元網格內速度的變化量。

通過上述方法，解決了DTMB 5512船模復雜的幾何船形與波浪作用的數值問題。

DTMB 5512型船模在長峰不規則波中縱向運動數值模擬的壓力場和速度場，如圖11和圖12所示。

圖12 DTMB 5512船模數值模擬速度場Fig.12 Numerical simulation speed field for the ship model of DTMB 5512

DTMB 5512型船模在長峰不規則波中縱搖及垂蕩的時歷曲線如圖13和圖14所示。

圖13 縱搖運動時歷曲線Fig.13 Time history of pitch movement

圖14 垂蕩運動時歷曲線Fig.14 Time history of heave movement

3 數值模擬結果驗證

由于本文對船模在長峰不規則波中頂浪縱向運動的數值模擬研究目前在國內外還處于起步階段，未找到具體的水池實驗數據。考慮到基于勢流理論艦船六自由度計算，在理論上是成熟的，并在工程應用上取得了很多成果，得到了水動力學界認可。一般來說，SCFD方法由于考慮到流體粘性，計算精度應略高于勢流理論計算結果，但沒有本質上的差異。至于非線性搖蕩則需要另作考慮，本課題局限于線性搖蕩，所以用勢流理論計算結果進行驗證［9］是可行的。

由于勢流理論切片法計算過程相對繁瑣，可參考文獻［1］，本文直接給出DTMB 5512型船模在算例波浪環境下縱搖和垂蕩的響應方差。

將數值模擬的DTMB 5512型船模搖蕩運動時歷曲線運用線性譜分析方法進行分析，獲得縱搖及垂蕩的搖蕩譜，如圖15和圖16所示。

圖15 縱搖運動響應譜Fig.15 Responsive spectrum for pitch movement

圖16 垂蕩運動響應譜Fig.16 Responsive spectrum for heave movement

對上述搖蕩譜進行積分即可得到本文數值模擬DTMB 5512船模頂浪運動縱搖和垂蕩的響應方差。

對DTMB 5512型船模在長峰不規則波中頂浪縱向運動的數值結果與勢流理論計算結果進行對比得出相對誤差如表6所示。

表6 DTMB 5512型船模相對誤差計算Tab.6 Relative errors for the ship model of DTMB 5512

4 結語

本文運用SCFD方法對艦船在長峰不規則波中頂浪運動進行數值模擬，得到如下可供參考的經驗:

1)船舶搖蕩、阻力、操縱、推進等課題的數值研究，技術細節上的一個主要不同，體現為網格的布設上。

2)船舶耐波性研究的網格，必須將造波和船體網格、動網格三者進行很好的協調、匹配，做到三者的有機結合，否則會出現計算發散及非物理現象等不合理現象的發生。

3)船舶多自由度搖蕩試驗目前的困難主要集中在計算資源上，由于其計算量巨大，PC機以及低端的工作站、服務器已經不能滿足其正常情況下的計算需要。如果計算資源等硬件設施有限，則需要在離散方法、格式，控制方程、調節參數、UDF開發等“軟件”上下功夫。

4)在三維空間內完成船模非規則波中的搖蕩試驗，其流場及自由面要比船模在規則波中的復雜得多，因此在船模貼體網格的布設上，要求網格質量非常高，同時在interface交界面處，左右網格要尺寸一致，且在船模外表面盡可能多地使用結構性網格，特別是在阻力計算上，還要盡可能多地布設邊界層，以提高計算的精度。

5)耐波性流場的高度復雜性還表現在輸入、輸出及船模的響應及其變化率上。具體表現為1個網格上，至少要輸入3個方向的線速度 u，v，ω，壓力P，參數k，ε，流體體積分數Cq等7個參數，及其輸出 u'，v'，ω'，P'，k'，ε'，C'q，并且要求它們之間的時空變化率在一定的范圍內。

6)船模外形的復雜程度也對船模多自由度搖蕩試驗產生一定的影響，當計算資源無法滿足大計算量的計算時，應盡量考慮用外形簡單、對稱性好的船模，如Wigley系列船模，以防止計算發散及非物理等不合理現象的發生。

7)耐波性波浪數值水池的網格布設同船模搖蕩的維數密切相關，而且同等條件下，在傅汝德數適中的情況下，計算效果較好。

8)在波浪的選擇及其他相關參數的選擇上，盡量參照《耐波性試驗章程》的相關要求，以防止非物理等不合理現象的發生。

［1］吳乘勝，朱德祥，顧民.數值波浪水池及頂浪中船舶水動力計算［J］.船舶力學，2008，12(2):168 －179.WU Cheng-sheng，ZHU De-xiang，GU Min.Computation of hydrodynamic forces for a ship in regular heading waves by a viscous wave tank［J］.Journal of Ship Mechanics，2008，12(2):168－179.

［2］董志，詹杰民.基于VOF方法的數值波浪水池及造波、消波方法研究［J］.水動力學研究與進展(A輯)，2009，24(1):15－21.DONG Zhi，ZHAN Jie-min.Comparison of existing methods for wave generating and absorbing in VOF-based numerical tank［J］.Journal of Hydrodynamics，2009，24(1):15 －21.

［3］GRILLI S T，VOGELMANN S，WATTS P.Development of a 3D numerical wave tank for modeling tsunami generation by underwater landslides［J］.Eng Analwith Boundary Elements，2002，26:301 －313.

［4］ARAI M，PAUL U K，CHENG L Y，et al.A technique for boundary treatment in numerical wave tanks［J］.Journal of the Society of Naval Architects of Japan，1993.

［5］吳乘勝，朱德祥，顧民.數值波浪水池中船舶頂浪運動數值模擬研究［J］.船舶力學，2008，12(5):692 －696.WU Cheng-sheng，ZHU De-xiang，GU Min.N-S CFD simulation of wave-induced ship motions in regular head waves［J］.Journal of Ship Mechanics，2008，12(5):692 －696.

［6］郭海強，朱仁傳，繆國平，余建偉.數值波浪水池中船舶水動力系數測試與分析技術［J］.中國造船，2008，49(S):58－65.

［7］RICCARDO B，ADREA D M.Unsteady RANS calculation of the flow around a moving ship hull［C］.In:Proc.of 8th International Conference in Numerical Ship Hydrodynamics，Busan:2003，A:153 －165.

［8］LUQUET R，JACQUIN E，GUILLERM P E，GENTZ L，et al.RANSE with free surface computations around fixed and free DTMB 5415 model in still water and in waves.Proc.the CFD Workshop Tokyo，2005.

［9］李積德.船舶耐波性［M］.哈爾濱:哈爾濱工程大學出版社，1992:1－10.