吳明 王驍 應(yīng)榮镕 侯建軍 石愛國(guó)

（海軍大連艦艇學(xué)院航海系 大連 116018）

船舶淺水航行下沉量和縱傾的數(shù)值計(jì)算

吳明 王驍 應(yīng)榮镕 侯建軍 石愛國(guó)

（海軍大連艦艇學(xué)院航海系 大連 116018）

基于CFD方法對(duì)船舶淺水航行航態(tài)變化問題進(jìn)行研究。通過同時(shí)求解RANS方程和剛體運(yùn)動(dòng)學(xué)方程對(duì)船舶淺水航行流場(chǎng)進(jìn)行數(shù)值模擬，并根據(jù)淺水流場(chǎng)空間受限的特點(diǎn)，采用三種動(dòng)網(wǎng)格相結(jié)合的方法，解決船舶在淺水中航行時(shí)航態(tài)變化的網(wǎng)格更新問題，在數(shù)值模擬中還考慮了螺旋槳旋轉(zhuǎn)的影響。文中以S60船模（Cb=0.6）為對(duì)象，計(jì)算了其淺水航行的下沉量和縱傾值，計(jì)算結(jié)果與船模水池試驗(yàn)數(shù)據(jù)對(duì)比，吻合良好。

計(jì)算流體力學(xué)；螺旋槳；淺水；下沉；縱傾

0 引言

船舶在淺水航行時(shí)，由于船體周圍流場(chǎng)的壓力和速度分布不均勻，會(huì)產(chǎn)生比較明顯的下沉現(xiàn)象，尤其是大型船舶在淺水高速航行時(shí)，船體下沉現(xiàn)象更加明顯，這已經(jīng)成為影響船舶淺水航行安全的重要因素之一。為了確保船舶航行安全，許多學(xué)者對(duì)船舶淺水航行的下沉現(xiàn)象進(jìn)行了理論研究和試驗(yàn)，最具代表性的是1967年，Tuck利用細(xì)長(zhǎng)體理論提出的計(jì)算公式［1］。后期在Tuck公式的基礎(chǔ)上，學(xué)者們又發(fā)展了許多解析式及半經(jīng)驗(yàn)的計(jì)算公式。第23屆ITTC操縱性委員會(huì)專門對(duì)淺水及受限水域船舶的操縱性問題進(jìn)行研究總結(jié)，并推薦一些實(shí)際廣泛應(yīng)用的計(jì)算公式［2-5］，對(duì)保證船舶淺水航行安全起到一定的指導(dǎo)作用。

雖然針對(duì)船舶淺水域航行下沉量的研究已取得不少成果，但這些研究多基于勢(shì)流理論方法，忽略了流體的粘性影響，而船舶淺水流場(chǎng)與船體及水底的壁面效應(yīng)密切相關(guān)；因此，采用勢(shì)流方法研究船舶淺水流場(chǎng)有一定的局限性。近年來，計(jì)算機(jī)技術(shù)的飛速發(fā)展為采用粘性流方法研究船舶水動(dòng)力學(xué)問題提供了有利的條件，粘性流方法在船舶阻力、耐波性、操縱性等方面取得了不少成果。然而，目前關(guān)于考慮浮態(tài)變化船舶淺水航行數(shù)值模擬研究比較少，主要原因是在淺水流場(chǎng)中，船底與水底之間的空間狹小，采用簡(jiǎn)單的動(dòng)網(wǎng)格方式，難以實(shí)現(xiàn)船體浮態(tài)變化的模擬。此外，現(xiàn)有的船舶淺水粘性流場(chǎng)數(shù)值模擬尚未考慮螺旋槳的影響，而受螺旋槳的排出流影響，船艉部流場(chǎng)會(huì)發(fā)生變化，這將會(huì)影響船舶在淺水中的航態(tài)，因此也是需要深入研究的問題。

針對(duì)上述問題，本文借助CFD商用軟件Fluent，采用同時(shí)求解描述船舶流場(chǎng)的RANS方程和描述船舶自由運(yùn)動(dòng)的剛體運(yùn)動(dòng)學(xué)方程的方法，結(jié)合軟件提供的三種動(dòng)態(tài)網(wǎng)格方法，實(shí)現(xiàn)帶螺旋槳船模淺水中航行流場(chǎng)的數(shù)值模擬，給出船模淺水航行時(shí)的下沉量和縱傾值，并與船模水池試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比驗(yàn)證。

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

船舶淺水流場(chǎng)可視為不可壓縮粘性流場(chǎng)，采用雷諾平均方法對(duì)船舶淺水域流場(chǎng)進(jìn)行數(shù)值模擬，流場(chǎng)連續(xù)性方程和動(dòng)量方程為：

式中：ui為流體時(shí)均速度分量，m/s；

p為流體壓強(qiáng)，Pa；

fi為流體體積力分量，N；

ρ為流體密度，kg/m3；

ν為流體的運(yùn)動(dòng)粘性系數(shù)，m2/s；

ui′為相對(duì)于時(shí)均流速的湍流脈動(dòng)速度分量，m/s；為雷諾應(yīng)力。

采用船舶計(jì)算流體力學(xué)中應(yīng)用較為廣泛的RNG k-ε二方程［6］湍流模型封閉上述方程組。該模型是由重正化群理論推導(dǎo)得出，其湍流動(dòng)能和湍流耗散率方程為：

式（3）和式（4）中各符號(hào)的物理意義詳見文獻(xiàn)［7］。

自由面的模擬采用VOF方法［8］來追蹤。設(shè)第q項(xiàng)流體的體積分?jǐn)?shù)為Cq，其輸運(yùn)方程為：

每一控制體的密度ρ為體積分?jǐn)?shù)的平均值：

式中：ρq為第q項(xiàng)流體的密度，單位kg/m3。

2 船舶運(yùn)動(dòng)數(shù)學(xué)模型

2.1 船舶運(yùn)動(dòng)控制方程

為了計(jì)算船舶在波浪中的升沉和縱傾運(yùn)動(dòng)，需定義三個(gè)坐標(biāo)系：第一個(gè)是固定坐標(biāo)系，其空間位置不變，流體的速度、壓力等在該坐標(biāo)系求解；第二個(gè)是平移坐標(biāo)系，其隨船體作平移運(yùn)動(dòng)；第三個(gè)是隨船坐標(biāo)系，其隨船體作旋轉(zhuǎn)運(yùn)動(dòng)，船體的搖蕩運(yùn)動(dòng)在該坐標(biāo)系求解。通常將船體坐標(biāo)系的原點(diǎn)取在船的重心上，把坐標(biāo)軸取作與船的中心慣性主軸相重合，這樣描述船舶搖蕩運(yùn)動(dòng)簡(jiǎn)單方便。

設(shè)Λ（φθψ）為船舶旋轉(zhuǎn)的歐拉角，對(duì)應(yīng)的分別為船舶的橫搖角、縱搖角和首搖角。如果只考慮角位置關(guān)系，向量在固定坐標(biāo)系與隨船坐標(biāo)系之間的轉(zhuǎn)換矩陣為：

設(shè)Ω=（p q r）為角速度在運(yùn)動(dòng)坐標(biāo)系的投影，則角速度在運(yùn)動(dòng)坐標(biāo)系的投影與歐拉角的導(dǎo)數(shù)之間的轉(zhuǎn)換關(guān)系為：

船體可視為剛體，船舶搖蕩運(yùn)動(dòng)遵循剛體運(yùn)動(dòng)的動(dòng)量定理和動(dòng)量矩定理。當(dāng)重心位于原點(diǎn)時(shí)，在運(yùn)動(dòng)坐標(biāo)系，船舶空間運(yùn)動(dòng)一般方程為：

式中：F=（X Y Z）為船體所受外力，N；L=（K M N）為船體所受外力合力對(duì)原點(diǎn)的力矩，N·m。

2.2 螺旋槳運(yùn)動(dòng)控制方程

在考慮螺旋槳的船舶搖蕩運(yùn)動(dòng)數(shù)值模擬中，螺旋槳的運(yùn)動(dòng)比較復(fù)雜，可分解為兩種運(yùn)動(dòng)：一是螺旋槳隨船體作搖蕩運(yùn)動(dòng)，二是螺旋槳繞轉(zhuǎn)軸做旋轉(zhuǎn)運(yùn)動(dòng)。由于螺旋槳與船體的相對(duì)位置固定，螺旋槳與船體的搖蕩運(yùn)動(dòng)是同步的，所以前者可以通過求解船體運(yùn)動(dòng)控制方程獲得。而螺旋槳繞轉(zhuǎn)軸的旋轉(zhuǎn)時(shí)，螺旋槳的軸向和旋轉(zhuǎn)中心是變量，因此需要給出瞬態(tài)的控制方程。

在全局坐標(biāo)系，設(shè)初始時(shí)刻螺旋槳旋轉(zhuǎn)中心的坐標(biāo)為A，船體繞轉(zhuǎn)心旋轉(zhuǎn)后，新的螺旋槳旋轉(zhuǎn)中心坐標(biāo)為A1=AS。

另取螺旋槳軸線上一點(diǎn)坐標(biāo)B，船體繞轉(zhuǎn)心旋轉(zhuǎn)后，該點(diǎn)坐標(biāo)為B1=BS，則新的螺旋槳軸向向量為：

2.3 數(shù)值求解步驟

船舶升沉和縱傾運(yùn)動(dòng)模擬的關(guān)鍵是如何實(shí)現(xiàn)流場(chǎng)與船舶運(yùn)動(dòng)的耦合求解，本文數(shù)值模擬的步驟如下：

（1）將船模按初始浮態(tài)固定，生成計(jì)算域和網(wǎng)格，對(duì)計(jì)算域進(jìn)行初始化；

（2）開始迭代計(jì)算，待流場(chǎng)基本穩(wěn)定時(shí)，加載船體和螺旋槳運(yùn)動(dòng)的方程；

（3）在當(dāng)前時(shí)間步長(zhǎng)Δt內(nèi)，根據(jù)船體所受垂向力Z和縱搖力矩M，通過求解船體運(yùn)動(dòng)控制方程，得到船體垂蕩位移的變化量ΔxCG和縱搖角的變化量ΔθCG；

（4）根據(jù)螺旋槳自轉(zhuǎn)角速度更新螺旋槳繞軸的旋轉(zhuǎn)角度，同時(shí)根據(jù)ΔxCG和ΔθCG更新船體和螺旋槳的位置和姿態(tài)；

（5）返回第（3）步，進(jìn)行下一時(shí)間段的步長(zhǎng)計(jì)算，求解航態(tài)改變后的流場(chǎng)，并根據(jù)計(jì)算結(jié)果再次更新螺旋槳繞軸的旋轉(zhuǎn)角度，改變船體和螺旋槳的位置和姿態(tài)，直到設(shè)定的時(shí)間步數(shù)計(jì)算完畢。

通過上述步驟可實(shí)現(xiàn)流場(chǎng)與船體運(yùn)動(dòng)的耦合求解，實(shí)現(xiàn)帶螺旋槳船模淺水航行升沉和縱傾運(yùn)動(dòng)的數(shù)值模擬。

3 數(shù)值模擬及結(jié)果分析

3.1 算例描述

選用S60船模進(jìn)行數(shù)值模擬，文獻(xiàn)［9］對(duì)該船模進(jìn)行淺水試驗(yàn)，給出阻力、縱傾和升沉測(cè)試結(jié)果。表1為S60船模的主尺度數(shù)據(jù)，圖1為本文采用GAMBIT軟件生成的三維船體模型。

表1 S60船模船型參數(shù)

圖1 S60船模

為了比較螺旋槳對(duì)船舶在淺水中航態(tài)的影響，本文分別對(duì)不帶螺旋槳和帶螺旋槳S60船模淺水域航行進(jìn)行數(shù)值模擬，不帶螺旋槳的數(shù)值試驗(yàn)7組，船模前進(jìn)速度U=1.107 4~1.701 8m/s（對(duì)應(yīng)的淺水傅汝德數(shù)Fh=0.5~0.8）；帶螺旋槳的數(shù)值試驗(yàn)，選用前蘇聯(lián)克雷洛夫船舶研究院的CK3-110型螺旋槳，船模前進(jìn)速度為U=1.328 8（對(duì)應(yīng)的淺水傅汝德數(shù)Fh=0.6），試驗(yàn)中參考物理水池船模自航試驗(yàn)方法，根據(jù)不同航速下的船模自航點(diǎn)，確定該航速下螺旋槳的轉(zhuǎn)速。

3.2 計(jì)算域及網(wǎng)格劃分

根據(jù)水池試驗(yàn)的條件設(shè)定計(jì)算域的尺寸，深吃水比H/T=2，水池寬度為2.09L，計(jì)算域的長(zhǎng)、寬、深分別為：4L×2.09L×2T，其中入口距船艏1L，出口距船艉2L，自由面距上、下邊界分別約0.25L和2T。

船舶搖蕩運(yùn)動(dòng)的模擬需要依靠動(dòng)網(wǎng)格技術(shù)來實(shí)現(xiàn)。Fluent軟件提供了三種基本的網(wǎng)格更新算法：彈簧近似光滑法（Spring Based Smoothing）、動(dòng)態(tài)分層法（Dynamic Layering）和局部重構(gòu)法（Local Remeshing），這三種算法通過拉伸、壓縮網(wǎng)格或者增加、減少網(wǎng)格以及局部重新生成網(wǎng)格來解決邊界運(yùn)動(dòng)的問題。

受限水域船舶下沉量數(shù)值計(jì)算的難點(diǎn)在于：船舶的升沉和縱傾運(yùn)動(dòng)，需要通過動(dòng)網(wǎng)格技術(shù)來實(shí)現(xiàn)；而限制水域船體底部與航道底面的空間區(qū)域狹小，使網(wǎng)格彈性形變和重構(gòu)的余地較小，在數(shù)值計(jì)算中極容易引起網(wǎng)格質(zhì)量的下降，導(dǎo)致計(jì)算發(fā)散。因此如何減小或者避免數(shù)值計(jì)算中的網(wǎng)格質(zhì)量下降，是受限水域船舶下沉量數(shù)值計(jì)算的關(guān)鍵。

本文采用三種動(dòng)網(wǎng)格方案相結(jié)合的方法解決上述問題。如圖2所示，跟據(jù)淺水流場(chǎng)及船?？v傾和升沉運(yùn)動(dòng)的特點(diǎn)，采用分區(qū)劃分網(wǎng)格的方法，生成了結(jié)構(gòu)/非結(jié)構(gòu)混合網(wǎng)格。在船體附近Z1區(qū)域生成非結(jié)構(gòu)網(wǎng)格，在該區(qū)域采用彈性形變和局部重構(gòu)法實(shí)現(xiàn)船體的縱搖運(yùn)動(dòng)；在Z2區(qū)域生成結(jié)構(gòu)網(wǎng)格，在該區(qū)域采用動(dòng)態(tài)分層法實(shí)現(xiàn)船體的升沉運(yùn)動(dòng)，同時(shí)Z1區(qū)域整體也隨Z2區(qū)域作升沉運(yùn)動(dòng)。采用這樣的動(dòng)網(wǎng)格方案可以減小網(wǎng)格彈性形變和重構(gòu)，從而避免數(shù)值計(jì)算過程中因網(wǎng)格質(zhì)量下降而引起計(jì)算發(fā)散。

圖2 網(wǎng)格劃分

3.3 離散格式及求解算法

采用有限體積法對(duì)微分方程進(jìn)行離散，壓力插值采用PRESTO格式，動(dòng)量、湍流動(dòng)能和湍流耗散率均采用二階迎風(fēng)格式離散，體積分?jǐn)?shù)采用Modified HRIC格式離散。采用SIMPLE算法對(duì)流場(chǎng)進(jìn)行求解。

3.4 數(shù)值模擬結(jié)果及分析

圖3為S60裸船模淺水中不同淺水傅汝德數(shù)下的升沉量計(jì)算值與船模試驗(yàn)值的比較圖，從圖中可以看出，隨著航速不斷增加，船體淺水航行的下沉量不斷增大。本文對(duì)淺水下沉量的計(jì)算與船模水池試驗(yàn)值相比吻合較好。

圖4為S60裸船模淺水中不同淺水傅汝德數(shù)下的縱傾計(jì)算值與船模試驗(yàn)值的比較圖。定義正的縱傾值表示船體產(chǎn)生向艉縱傾。從圖中可以看出，在Fh＜0.75時(shí)，船體有微小的向艏縱傾；Fh＞0.75時(shí)，船體開始產(chǎn)生向艉縱傾，通過比較可以看出，本文方法對(duì)縱傾的計(jì)算值與水池試驗(yàn)值吻合得較好。

圖3 淺水中的升沉

圖4 淺水中的縱傾

圖5 升沉?xí)r歷（Fh=0.6）

圖6 縱傾時(shí)歷（Fh=0.6）

圖5和圖6中虛線為S60裸船模的淺水航行下沉量和縱傾的水池試驗(yàn)值，從圖中可以看出，與S60裸船模淺水水池試驗(yàn)數(shù)據(jù)相比，螺旋槳的旋轉(zhuǎn)使船舶淺水航行中的船體下沉量和艉傾增大。

圖7和圖8為淺水傅汝德數(shù)為0.6，船體航態(tài)固定時(shí)，S60裸船模與帶螺旋槳船模艉部船底的壓力分布對(duì)比圖。從圖中可以看出，盡管帶槳船體在淺水中運(yùn)動(dòng)的流場(chǎng)十分復(fù)雜，但螺旋槳的排出流使緊貼船艉濕表面的作用力降低。由于船艉的壓力是表面積分的結(jié)果，因此會(huì)導(dǎo)致船體向艉縱傾力矩增大，進(jìn)而使船體下沉量和艉傾增大。淺水流場(chǎng)中船體的壓力分布變化是引起船體的下沉量和縱傾變化的重要原因之一，本文的數(shù)值模擬可以捕捉到淺水流場(chǎng)變化的細(xì)節(jié)，為研究船舶淺水流場(chǎng)提供了一種有效的方法。

圖7 裸船模艉底部表面壓力等值線

圖8 帶螺旋槳船模艉底部表面壓力等值線

4 結(jié)論

本文基于粘性流理論，實(shí)現(xiàn)帶螺旋槳船模淺水域航行升沉和縱傾運(yùn)動(dòng)的數(shù)值模擬。數(shù)值計(jì)算結(jié)果表明：

（1）借助CFD軟件，采用三種動(dòng)網(wǎng)格相結(jié)合的方法，可以實(shí)現(xiàn)考慮船體浮態(tài)變化的船舶淺水航行的數(shù)值模擬；計(jì)算得到的S60裸船模淺水航行下沉量和縱傾與船模水池結(jié)果吻合良好，可為淺水域航行下沉量計(jì)算提供參考。

（2）螺旋槳的旋轉(zhuǎn)會(huì)使船艉附近流場(chǎng)發(fā)生變化，引起船艉底部壓力下降，這是船體艉縱傾及下沉量增大的重要原因之一。

本文對(duì)帶螺旋槳船舶淺水域航行下沉量問題進(jìn)行了初步研究。模擬的船舶航行工況較少，在今后的工作中將進(jìn)一步研究螺旋槳轉(zhuǎn)速同船舶航行下沉量及縱傾之間的關(guān)系。

［1］TUCK E O.Sinkage and trim in shallow water of finite width［J］.Schiffstechnik，1967，14：92-94.

［2］HOOFT JP.The be haviour of a ship in head waves at restricted water depth［J］.International Shipbuilding Progress，1974，21：367.

［3］BARRASSC B.A unified approach to squat calculations for ships［J］.PIANC Bulletin，1979，32：3-10.

［4］DAND IW，F(xiàn)erguson A.The squat of full ships in shallow water［J］.The Naval Architect，1973（4）：237-255.

［5］ANKUDINOV V，DAGGETT L，HUVALL C，et al.Squat predictions for maneuvering applications［C］//Proceedings ofMARSIM’96，Copenhagen，［s.n.］，1996：467-495.

［6］YAKHOT V，ORSZAG SA，THANGAM S，et al.Development of turbulence models for shear flows by a double expansion technique［J］.Physics of Fluids A，1992（4）：1510-1520.

［7］王福軍.計(jì)算流體動(dòng)力學(xué)分析［M］.北京：清華大學(xué)出版社，2004：124.

［8］HIRTCW，NICHOLSBD.Volume of fluid（VOF）method for the dynamics of free boundaries［J］.Journal of Computational Physics，1981，39：201-225.

［9］JIANG T.Investigation of waves generated by ships in shallow water［C］//Proceeding of the 22nd Symposium on Naval Hydrodynamics，Washington DC：［s.n.］，1998：601-612.

Numerical calculation of ship sink age and trim in shallow water

WU Ming WANG Xiao YING Rong-rong HOU Jian-jun SHIAi-guo
（Dalian Naval Academy，Department of Navigation，Dalian 116018，China）

A study on the navigation behavior for ships in shallow water has been carried out based on CFD.This paper numerically simulates flow field around ship in shallow water by simultaneously solving equations of Reynolds Averaged Navier-Stokes（RANS）and kinematic equations of a rigid body.To settle the mesh update issue for the ships in shallow water when the navigation behavior is changed，this paper adopts the method of coupling three dynamic mesh according to the characteristics of limited space of shallow water.It also considers the influence of the propeller rotation.The numerical simulation is carried out on a S60model（Cb=0.6）to calculate the ship sinkage and trim travelling in shallow water.The numerical results are in good agreement with the experimental data.

CFD；ship；propeller；shallow water；sinkage；trim

U661.3

A

1001-9855（2013）03-0007-05

遼寧省博士科研啟動(dòng)基金項(xiàng)目（20111037），遼寧省優(yōu)秀青年人才基金項(xiàng)目（2012J21DW027）。

2013-01-10；

2013-03-25

吳明（1984-），男，博士，研究方向：船舶耐波性。

王驍（1980-），男，博士，講師，研究方向：艦船操縱性。

應(yīng)榮镕（1984-），男，碩士，講師，研究方向：艦船操縱性。

侯建軍（1973-），男，博士，高級(jí)工程師，研究方向：艦船耐波性。

石愛國(guó)（1956-），男，教授，研究方向：艦船操縱性。