廉靜靜，楊 曉，尹 勇

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

0 引 言

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

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

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

1.1 控制方程

求解各種計(jì)算模型流場(chǎng)及流體動(dòng)力所采用的基本方法是RANS方法[7]，其RANS方程的形式如下：

(1)

1.2 湍流模型

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

ρε-YM+Sκ

(2)

(3)

1.3 Q準(zhǔn)則

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

(4)

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

2 數(shù)值計(jì)算

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

圖1 KVLCC2 船模Fig. 1 KVLCC2 ship model

m

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

圖2 計(jì)算域Fig. 2 Computational domain

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

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

3 結(jié)果與分析

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

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

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

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

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

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

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

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

圖8 船體周圍的漩渦流場(chǎng)Fig. 8 Vorticity flow distribution around ship

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

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

4 結(jié) 語(yǔ)

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

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

參考文獻(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 船模斜航運(yùn)動(dòng)黏性流場(chǎng)及水動(dòng)力數(shù)值計(jì)算[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)水槽平面運(yùn)動(dòng)機(jī)構(gòu)試驗(yàn)及其斜航試驗(yàn)的數(shù)值模擬[J]. 大連海事大學(xué)學(xué)報(bào), 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] 王福軍.計(jì)算流體動(dòng)力學(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.