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

劉 可，吳 明，楊 波

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

0 引言

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

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

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

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

1.1.1 數(shù)值造波

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

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

1.1.2 數(shù)值消波

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

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

在阻尼消波段內(nèi)，動量方程寫為:

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

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

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

1.2.1 網(wǎng)格劃分

三維數(shù)值水槽的網(wǎng)格劃分如圖1所示，本文構(gòu)建的長峰不規(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)格尺度為基準(zhǔn)向右邊界逐漸擴大。

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

1.2.2 算例描述

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

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

1.2.3 數(shù)值模擬結(jié)果

長峰不規(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軸正方向向下游傳播過程中，波能有一定程度的衰減，這是由于波浪在傳播過程中，高頻子波衰減導(dǎo)致的。

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

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

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

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

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

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

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

2.1 船模的選擇

選取具有球鼻首和方位的DTMB 5512船模作為研究對象，該船模是ITTC(國際船模試驗水池會議)推薦的瘦削型標(biāo)準(zhǔn)船模DTMB 5415的全相似幾何模型，以美國海軍DDG－51型驅(qū)逐艦為模板，船型數(shù)據(jù)詳實，而且與我海軍艦船船型相似 (見圖7)，適合軍艦參考。表4為該船模及其所對應(yīng)的實船尺度的主要數(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ī)程》，將計算域設(shè)置成長方體形狀，如圖8所示。船模與計算域各邊界的位置關(guān)系如下:入口距船首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)，當(dāng)波浪作用于船體時，其運動的速度、角速度以及位置、姿態(tài)等可以通過控制方程求解、積分得到。對于流浮耦合運動，波浪作用在船體上的力和力矩使船體產(chǎn)生運動，同時船體的運動又對其周圍流場產(chǎn)生影響。因此本文在數(shù)值模擬中，分段計算流體與船體運動的耦合，步驟如下:

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

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

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

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

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

6)求解體積分?jǐn)?shù)方程重構(gòu)自由面;

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

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

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

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

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

2.4 數(shù)值模擬結(jié)果

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

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

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

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

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

通過上述方法，解決了DTMB 5512船模復(fù)雜的幾何船形與波浪作用的數(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ù)值模擬結(jié)果驗證

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

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

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

圖15 縱搖運動響應(yīng)譜Fig.15 Responsive spectrum for pitch movement

圖16 垂蕩運動響應(yīng)譜Fig.16 Responsive spectrum for heave movement

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

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

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

4 結(jié)語

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

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

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

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

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

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

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

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

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

［1］吳乘勝，朱德祥，顧民.數(shù)值波浪水池及頂浪中船舶水動力計算［J］.船舶力學(xué)，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］.水動力學(xué)研究與進展(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］.船舶力學(xué)，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ù)測試與分析技術(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］.哈爾濱:哈爾濱工程大學(xué)出版社，1992:1－10.