船舶在長峰不規(guī)則波中頂浪縱向運動的數(shù)值模擬

劉 可，吳 明，楊 波

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

0 引言

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

1 波浪的數(shù)值模擬

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

1.1 數(shù)值波浪水池的關鍵技術

1.1.1 數(shù)值造波

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

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

1.1.2 數(shù)值消波

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

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

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

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

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

1.2 長峰不規(guī)則波數(shù)值波浪水池的構建

1.2.1 網(wǎng)格劃分

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

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

1.2.2 算例描述

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

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

1.2.3 數(shù)值模擬結果

長峰不規(guī)則波數(shù)值模擬瞬時波面場和局部速度矢量如圖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時，長峰不規(guī)則波數(shù)值波浪水池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

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

1.2.4 數(shù)值模擬精度誤差計算

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

圖6 數(shù)值模擬海浪譜Fig.6 Numerical simulation spectrum

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

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

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

2 船模在長峰不規(guī)則波中頂浪縱向運動的數(shù)值模擬

2.1 船模的選擇

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

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

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

2.2 數(shù)值模擬方案

2.2.1 試驗計劃

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

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

2.2.2 計算域劃分

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

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

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

圖9 耐波性波浪數(shù)值水池網(wǎng)格劃分Fig.9 Build grid for sea-keeping numerical simulation wave tank

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

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

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

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

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

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

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

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

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

8)更新船體位置和浮態(tài);

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

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

通過以上迭代、循環(huán)，可實現(xiàn)流體與船體運動的耦合。

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

2.4 數(shù)值模擬結果

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

1)改進船模貼體網(wǎng)格質量、數(shù)量及分區(qū)網(wǎng)格的匹配;

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

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

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

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

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

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

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

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

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

3 數(shù)值模擬結果驗證

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

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

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

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

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

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

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

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

4 結語

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

1)船舶搖蕩、阻力、操縱、推進等課題的數(shù)值研究，技術細節(jié)上的一個主要不同，體現(xiàn)為網(wǎng)格的布設上。

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

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

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

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

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

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

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

［1］吳乘勝，朱德祥，顧民.數(shù)值波浪水池及頂浪中船舶水動力計算［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方法的數(shù)值波浪水池及造波、消波方法研究［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］吳乘勝，朱德祥，顧民.數(shù)值波浪水池中船舶頂浪運動數(shù)值模擬研究［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］郭海強，朱仁傳，繆國平，余建偉.數(shù)值波浪水池中船舶水動力系數(shù)測試與分析技術［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.