基于數(shù)值分析的船舶斜航運動水動力研究

廉靜靜，楊 曉，尹 勇

(1.大連海事大學(xué) 航海學(xué)院，遼寧 大連 116026； 2. 大連海事大學(xué) 航海動態(tài)仿真和控制實驗室，遼寧 大連 116026)

0 引 言

船舶操縱性能與船舶航行的安全性緊密相關(guān)，為了提高船舶的航行安全，負責(zé)海事安全的國際海事組織IMO(International Maritime Organization)于1993年頒布實施了InterimStandardsforShipManeuverability[1]以及2002年通過了StandardsforShipManeuverability[2]，這使得船舶的操縱性能越來越受到關(guān)注和重視。近年來，由于采用數(shù)值計算方法能夠提供船體和螺旋槳葉面周圍獨特的定性和定量流場信息，它已成為求取船舶操縱性水動力導(dǎo)數(shù)主流方法之一，可將取得船舶操縱性水動力導(dǎo)數(shù)代入船舶操縱運動方程，對船舶的操縱性實現(xiàn)預(yù)報，同時，采用數(shù)值計算方法也逐步成為船模實驗方法的主要輔助手段。

利用計算流體動力學(xué)(computational fluid dynamics，CFD)方法研究船舶水動力性能，國內(nèi)外學(xué)者均開展了相關(guān)研究工作 ， B. J. GUO等[3]研究了在頂浪中KVLCC2船舶運動的附加阻力；Y. T. JIN等[4]考慮了船舶尺度大小對KVLCC2船舶水動力系數(shù)的影響；田喜民等[5]，劉晗等[6]對KVLCC2不同漂角的斜航運動進行了模擬，但這些研究對船體周圍漩渦流場分布情況并未關(guān)注，且船體周圍漩渦對船體的阻力、噪聲密切相關(guān)，從而對船舶操縱性產(chǎn)生影響。筆者選用國際操縱性標(biāo)準(zhǔn)船型KVLCC2為研究對象，采用ANSYS FLUENT 14.5軟件對該船模的斜航運動的黏性流場進行了數(shù)值模擬，對黏性流場中數(shù)值模擬中的計算網(wǎng)格、湍流模型等關(guān)鍵技術(shù)問題進行了深入研究，獲得作用在船舶上的縱向水動力、橫向水動力、艏搖力矩，并將計算的數(shù)值結(jié)果與日本海上技術(shù)安全研究所(National Maritime Research Institute，NMRI)實驗結(jié)果進行比較分析，同時采用最小二乘法進行曲線擬合，求得船舶操縱水動力導(dǎo)數(shù)，采用Q準(zhǔn)則獲得了船舶周圍的漩渦流場。

1 數(shù)學(xué)模型

1.1 控制方程

求解各種計算模型流場及流體動力所采用的基本方法是RANS方法[7]，其RANS方程的形式如下：

(1)

1.2 湍流模型

分析比較了幾種湍流模型的適用性，筆者選用的湍流模型為重整化群RNG(Renormalization Group)κ-ε模型，與求解高雷諾數(shù)的標(biāo)準(zhǔn)Standard κ-ε湍流模型相比，重整化群RNG κ-ε湍流模型考慮了低雷諾數(shù)的影響，通過在大尺度運動和修正后的湍流黏度項體現(xiàn)了小尺度的影響，使得在平均流動中的旋轉(zhuǎn)及漩流流動情況得以處理，同時，在重整化群RNG κ-ε模型中引入時均應(yīng)變率Eij，從而使得高應(yīng)變率和流線彎曲程度較大的流動該模型能夠更好地處理。該模型由V. YAKHOT等[8]于1986年提出的。該湍流模型是通過重整化數(shù)學(xué)方法推導(dǎo)出瞬時的NS方程，重整化群RNG κ-ε模型方程[9]為

ρε-YM+Sκ

(2)

(3)

1.3 Q準(zhǔn)則

Q準(zhǔn)則是由J. C. R. HUNT等[10]提出的，通過使用正值Q來作為流場漩渦的識別方法，當(dāng)Q值為正時，說明該流場區(qū)域內(nèi)的流體微元以旋轉(zhuǎn)運動為主導(dǎo)，當(dāng)Q值為負時，則說明該流場區(qū)域內(nèi)的流體微元的旋轉(zhuǎn)運動較弱。定義流場中速度梯度張量的第二矩陣不變量Q可表達為

(4)

式中：Wij為速度梯度張量的旋轉(zhuǎn)率張量，代表純粹旋轉(zhuǎn)運動；Sij為速度梯度張量應(yīng)變率張量。

2 數(shù)值計算

筆者研究對象為國際標(biāo)準(zhǔn)船型KVLCC2油輪，船舶模型如圖1，其船模的主要參數(shù)如表1。

圖1 KVLCC2 船模Fig. 1 KVLCC2 ship model

m

首先構(gòu)建計算域的幾何模型，即求解計算時的離散積分空間，依據(jù)船舶斜航運動的物理問題，確定計算域的大小。計算域如圖2(以下船長記為L)，距離船首1.5L設(shè)置為上游邊界，距離船尾3L設(shè)置為下游邊界，向下距離1.5L處設(shè)置為水深方向，左右兩側(cè)方向為1.5L。只有根據(jù)物理條件，給出合理的邊界條件，才能獲得流場的正確解，上游邊界是速度進口邊界，船舶航速是0.994 m/s，對應(yīng)的Froude數(shù)為0.142 4，為低速船，自由面興波影響較小，故采用“疊?！蹦Ｐ?，將自由面視為對稱邊界，船舶的兩側(cè)和計算域的底面視為無滑移壁面邊界條件，船體看作無滑移壁面邊界條件，下游邊界看作出流outflow邊界。

圖2 計算域Fig. 2 Computational domain

采用有限體積法對微分方程進行離散時，需要將控制方程在空間區(qū)域上進行離散，隨后求解得到離散方程組，需要生成網(wǎng)格將控制方程在空間區(qū)域上進行離散，且網(wǎng)格質(zhì)量和數(shù)量直接影響計算精度與計算速度，可通過前處理軟件ICEM CFD對整個計算域進行網(wǎng)格劃分，如圖3船舶表面網(wǎng)格。在整個網(wǎng)格劃分過程中，盡量使用結(jié)構(gòu)網(wǎng)格以減少網(wǎng)格數(shù)量，提高計算速度，在船舶球鼻首和船舶尾部的地方采用結(jié)構(gòu)網(wǎng)格比較困難，可以利用非結(jié)構(gòu)網(wǎng)格進行劃分，網(wǎng)格生成方法采用自上而下的方法，船舶首尾部曲率大的線條部分線網(wǎng)格加密尺寸小，使得該部分網(wǎng)格較為密集，相比之下，船舶首尾部和舭部面網(wǎng)格加密尺寸較大，以便減少整個計算域內(nèi)網(wǎng)格數(shù)量。當(dāng)處理近壁面流動時，邊界層內(nèi)的位置通常用無量綱形式y(tǒng)+來表示，y+是以垂直于壁面的距離為長度尺度，即為第一個網(wǎng)格點到壁面距離，此處將邊界層棱柱網(wǎng)格y+可設(shè)置為30，整個計算域網(wǎng)格總數(shù)為70萬左右，圖4為整個計算域的網(wǎng)格劃分遠觀圖。為提高計算效率，計算域內(nèi)網(wǎng)格數(shù)量不能過多，但考慮精度，某些區(qū)域的網(wǎng)格數(shù)量還需加密，主要依據(jù)計算時耗時和計算的精度，從兩者中折中確定網(wǎng)格數(shù)量。采用壓力耦合方程組的半隱式SIMPLE算法求解連續(xù)方程中建立的壓力修正方程，得到壓力場的修正值，利用壓力修正值更新速度場和壓力場，不收斂時重復(fù)迭代，最終得到壓力場和速度場的收斂解，為了保證計算過程中空間離散不發(fā)散，在初始階段，采用一階迎風(fēng)格式離散動量、湍流動能和湍流耗散率等參數(shù)，待計算穩(wěn)定后采用二階迎風(fēng)格式即可，這樣即可以保證穩(wěn)定性又提高準(zhǔn)確性，同時，為保證計算的收斂性，改變亞松弛因子。

圖3 船舶體表網(wǎng)格Fig. 3 Ship hull surface mesh

3 結(jié)果與分析

筆者計算了KVLCC2斜航船模漂角分別為0°、3°、6°、9°和12°時的縱向力系數(shù)、橫向力系數(shù)和艏搖力矩系數(shù)，如圖5，并將計算得到的結(jié)果與日本NMRI獲得的試驗結(jié)果[11]進行對比，從圖5中可以看出，實驗結(jié)果吻合較好，漂角較小時，數(shù)值計算的船舶的橫向力系數(shù)與試驗結(jié)果比較誤差約為6%左右，在整個計算過程中，船舶縱向力系數(shù)與試驗結(jié)果誤差較大，約為12%左右，誤差與離散時間步長、離散方法和網(wǎng)格疏密等因素有關(guān)，湍流動能值和湍流耗散值會對橫向阻力系數(shù)和縱向阻力系數(shù)產(chǎn)生影響，選擇合適的湍流動能值和湍流耗散值是非常重要的。

圖5 縱向力系數(shù)、橫向力系數(shù)及艏搖力矩系數(shù)Fig. 5 Longitudinal force coefficientlateral force coefficient and yaw moment coefficient

根據(jù)最小二乘法原理進行曲線擬合，求得船舶水動力導(dǎo)數(shù)，如表2，圖6為不同側(cè)向速度下水動力導(dǎo)數(shù)的曲線擬合結(jié)果。

表2 船舶水動力導(dǎo)數(shù)Table 2 Ship hydrodynamic derivatives

圖曲線和曲線擬合Fig. curve and curve fitting

將漂角為12°時船體艏部和尾部表面壓力系數(shù)分布與日本NMRI[12]實驗數(shù)據(jù)進行比較如圖7。從圖7中可以看出采用數(shù)值計算的結(jié)果和實驗數(shù)據(jù)結(jié)果相似，在船首部迎風(fēng)面所受的壓力較大，背風(fēng)面所受壓力較小，而船尾則相反。兩者不同之處為，實驗數(shù)據(jù)在船舶尾部背風(fēng)面和迎風(fēng)面壓力系數(shù)曲線中某條線出現(xiàn)了閉合，數(shù)值模擬并非發(fā)現(xiàn)該現(xiàn)象。

圖7 β=12°時船體表面壓力系數(shù)分布Fig. 7 Hull pressure coefficient distribution(β=12°)

通常Q準(zhǔn)則可觀測到船舶周圍的渦結(jié)構(gòu)，圖8是漂角為0°、3°、6°、9°和12°時船舶周圍的漩渦流場。

圖8 船體周圍的漩渦流場Fig. 8 Vorticity flow distribution around ship

從圖8中可以看出，隨著漂角逐漸增加，船舶周圍的漩渦流場變得更加復(fù)雜，尤其是船舶尾部錐形部分和船中部周圍流場高度復(fù)雜，這是因為隨著漂角逐漸增加從船舶中部側(cè)漩渦越發(fā)明顯，逐漸出現(xiàn)分離，這使得漩渦結(jié)構(gòu)得到了充分發(fā)展。當(dāng)漂角為0°時，船舶KVLCC2周圍的漩渦流場受到船首和船尾錐形部分邊界層發(fā)展支配。當(dāng)漂角為12°時，船首舭部漩渦、船首側(cè)渦、船尾漩渦、側(cè)漩渦逐漸形成，同時，船舶尾部漩渦在迎風(fēng)面與船舶尾部側(cè)渦合并在一起，如圖9，在背風(fēng)面與船舶尾部發(fā)卡渦合并在一起，在船首舭部漩渦看起來有些卡門渦脫落，而船舶首部側(cè)漩渦有些螺旋不穩(wěn)定，產(chǎn)生了一些螺旋波動和螺旋漩渦。

圖9 船尾局部漩渦流場(β=12°)Fig. 9 Vorticity flow distribution around ship stern(β=12°)

4 結(jié) 語

通過求解不可壓縮黏性流的RANS方程，結(jié)合重整化群RNGκ-ε湍流模型，對KVLCC2船舶斜航運動的黏性流場進行數(shù)值模擬，計算不同漂角下船舶所受的縱向水動力、橫向水動力、艏搖力矩，并將其計算獲得的結(jié)果與實驗數(shù)據(jù)相比較，吻合較好，計算結(jié)果產(chǎn)生的誤差在允許的范圍內(nèi)，湍流動能值、湍流耗散值、離散時間步長、離散方法和網(wǎng)格疏密對計算結(jié)果有影響，同時采用最小二乘法進行曲線擬合，求得船舶操縱水動力導(dǎo)數(shù)，采用Q準(zhǔn)則獲得了船舶周圍的漩渦流場，當(dāng)船舶漂角逐漸增加時，船體周圍

的漩渦越發(fā)復(fù)雜。采用數(shù)值分析的方法求取水動力導(dǎo)數(shù)，將計算的結(jié)果帶入船舶運動方程中，為船舶CFD實現(xiàn)實用化提供一種有效的思路。

參考文獻(References)：

[1] International Maritime Organization (IMO).ResolutionA. 751(18):InterimStandardsforShipManeuverability[S]. London: IMO, 1993.

[2] International Maritime Organization (IMO).ResolutionMSC.137 (76):StandardsforShipManeuverability[S]. London: IMO, 2002.

[3] GUO B J, STEEN S, DENG G B. Seakeeping prediction of KVLCC2 in head waves with RANS[J].AppliedOceanResearch, 2012, 35:56-67.

[4] JIN Y T, DUFFY J, CHAI S H,et al. URANS study of scale effects on hydrodynamic manoeuving coefficients of KVLCC2[J].OceanEngineering, 2016, 118:93-106.

[5] 田喜民, 鄒早建, 王化明. KVLCC2 船模斜航運動黏性流場及水動力數(shù)值計算[J]. 船舶力學(xué), 2010,14(8):834-840.

TIAN Ximin, ZOU Zaojian, WANG Huaming. Computation of the viscous flow and hydrodynamic forces on a KVLCC2 model in oblique motion[J].JournalofShipMechanics,2010,14(8): 834-840.

[6] 劉晗,馬寧,鄧德衡,等. 循環(huán)水槽平面運動機構(gòu)試驗及其斜航試驗的數(shù)值模擬[J]. 大連海事大學(xué)學(xué)報, 2016,42(2):8-13.

LIU Han, MA Ning DENG Deheng ,et al. Planar motion mechanism tests in circulating water channel and numerical simulation of oblique towing[J].JournalofDalianMaritimeUniversity, 2016, 42(2):8-13.

[7] 王福軍.計算流體動力學(xué)分析——CFD軟件原理與應(yīng)用[M].北京:清華大學(xué)出版社,2004:7-10.

WANG FuJun.ComputationalFluidDynamicsanalyzer:softwareprincipleandapplication[M].Beijing: Tsinghua University Press, 2004:7-10.

[8] YAKHOT V, ORSZAG S A, PANDA R. Computational test of the renormalization group theory of turbulence[J].JournalofScientificComputing, 1988, 3(2):139-147.

[9] ORSZAG S A, YAKHOT V, FLANNERY W S, et al. Renormalization group modeling and turbulence simulations[C]//InternationalConferenceonNear-WallTurbulentFlows,Tempe,Arizona. Tempe, Arizona:[s.n.], 1993.

[10] HUNT J C R, WRAY A A, MOIN P. Eddies, stream and convergence zones in turbulent flows[C]//StudyingTurbulenceUsingNumericalSimulationDatabases,Proceedingsofthe1988SummerProgram. United States:[s.n.], 1988: 193-208.

[11] JACQUIN E, GUILLERM P E, DROUET A, et al. Simulation of unsteady ship maneuvering using free-surface RANS solver[C]//Proceedingof26thSymposiumonNavalHydrodynamics. Rome,Italy:[s.n.],2006.

[12] KUME K, HASEGAWA J, TSUKADA Y, et al. Measurements of hydrodynamic forces, surface pressure, and wake for obliquely towed tanker model and uncertainty analysis for CFD validation[J].JournalofMarineScienceandTechnology, 2006, 11(2): 65-75.