梭型車客渡船水動力特性數值預報研究

呂全亮，劉雯玉，徐 犇，劉鵬鑫

(1.江蘇省鎮江船廠(集團)有限公司，江蘇 鎮江 212000；2.江蘇科技大學 船舶與海洋工程學院，江蘇 鎮江 212003；3.江蘇省船舶設計研究所有限公司，江蘇 鎮江 212003)

0 引言

隨著計算機和信息技術的發展、湍流模型的突破、非線性波浪數字模擬技術與動網格技術的運用，基于計算流體動力學(CFD)的船舶阻力性能研究技術愈加成熟。

許煬塏等[1]利用CFD軟件，針對3種不同湍流模型，對KCS船的阻力性能進行了預報，通過與船模試驗數據進行比較，得到了最適合速度范圍內的湍流模型。倪崇本[2]基于粘流理論，將三因次船體阻力預報方法與CFD數值模擬手段進行了有效的結合，提出了基于CFD數值模擬計算的一種快速預報實船阻力的方法。王金寶等[3]對低速肥大型船舶進行了阻力性能研究，并對尾流場進行了數值模擬，將其方法應用于其他實船，結果較為精確。蘇紹娟[4]等采用數值模擬進行船舶阻力性能預報，發現在恰當的模型與高質量的網格情況下，CFD數值模擬能很好且快速地模擬船舶航行狀態。王健等[5]利用Fluent軟件分別對高速單體船和高速三體船進行了數值模擬，預報了各阻力并分析了各流場的變化，為多體船的研究做出了貢獻。魏可可等[6]以5415船模為研究對象，用STAR-CCM+對其船模阻力、升沉及橫傾角進行了計算，并與試驗值進行了對比，得到了總阻力系數、剩余阻力系數隨弗勞德數Fr增大而增大的結論。查若思等[7]基于RANS方程的黏性數值模擬方法，利用OpenFOAM對多種船型在不同航速下靜水阻力和興波情況進行了模擬，分析了各船型快速性。JOSIP[8]利用流體體積表面捕獲技術的RANS方程，預測了受損船舶在靜水中前行時周圍的流場，并將結果與Brodarski研究所的模型實驗進行對比，預測了受損船舶的總阻力性能。

以往對于船舶阻力性能多限于實船，且對于船模數據與實船數據的轉換多依靠于二因次法。本文首先基于粘性流體理論，采用隨體動網格技術模擬船舶航行狀態；其次通過三因次法進行船模數據換算，與試驗數據對比驗證了數值算法及換算方法的可靠性，完成了不同載況下客渡船靜水中阻力性能的數值預報研究，采用理論計算給出了不同載況下車客渡船的有效功率，對車客渡船的推進器選型具有指導意義。

1 數值計算方法

1.1 基本控制方程

流體流動要受物理守恒定律的支配，最基本的守恒定律包括：質量守恒定律、動量守恒定律及能量守恒定律，在流體力學中具體體現為連續性方程、動量方程和能量方程。水通常被當做不可壓縮流體，且本文也不考慮能量的交換，因此，數值計算僅需要滿足連續方程和動量方程。

(1)連續性方程

(1)

(2)

(2)動量守恒方程

(3)

式中：p為流體微元體上的壓力；τxx、τxy、τxz、τyx、τyy、τyz、τzx、τzy、τzz為因分子粘性作用在微元體表面上的粘性應力τ的分量；Fx、Fy、Fz為微元體上的體力，若只受重力，且z軸豎直向上，則Fx=0,Fy=0,Fz=-ρg。

對于不可壓縮流體，該方程可進一步簡化為：

(4)

式中：μ為k的粘性系數。

1.2 湍流模型

本文采用混合了計算近壁區域粘性流動的k-ω模型和用來計算遠場自由流動的k-ε模型的SSTk-ω模型。在這個模型中，湍流動能k和特定耗散率ω的輸運方程為：

(5)

Yω+Dω+Sω

(6)

與標準k-ω模型相同，其中Γk、Γω分別為k和ω的擴散系數，且表達式相同。其中，μt的定義為：

這里，普朗特數σk與σω不再為常數，其定義為：

上列兩式中F1與F2的定義為：

(7)

Gk和Gω為k和ω方程的產生項，其定義為：

其中，

α∞=F1α∞,1+(1-F1)α∞,2；

Yk、Yω為k和ω方程的耗散項，其定義為：

Yk=ρβ*kω；Yω=ρβω2；β1=F1βi,1+(1-F1)βi,2

Dω為橫向擴散項，是將k-ε模型轉換成了k-ω模型形式產生的，定義為：

1.3 離散方法

計算流體動力學的核心內容就是離散方法，常用的離散方法為有限差分法、有限元法和有限體積法。本文采用有限體積法進行計算，其控制方程使得任何一組控制體積都能滿足特征變量的積分守恒，即使在網格較粗的條件下，也能準確表現出積分守恒。

1.4 數值方法驗證

在進行正式計算之前，需要對以上數值方法進行驗證。選取KCS標準船型進行縮尺計算，縮尺比為31.6。其幾何模型見圖1，主要參數見表1。

將縮尺后的KCS船體模型導入HEXPRESS軟件進行網格劃分。建立計算域長寬高分別為4Lpp、3Lpp、2.5Lpp，船模距前方入口1Lpp，距出口3Lpp，距上下邊界分別0.5Lpp和2Lpp，左右邊界各1.5Lpp。計算域內船體采用標準壁面函數，上下邊界設置為指定壓強(液體靜壓)，入口邊界、出口邊界及兩側邊界設定為遠場邊界。為更精確地得到流場數據，需要對船體球鼻艏、水線面、尾流場進行局部加密。計算域劃分與船體網格分別見圖2、圖3。

圖1 KCS船體模型

表1 KCS船主要參數

圖2 計算域網格劃分

圖3 船體表面網格

得到縮尺船模數據之后需要將船模阻力換算成實船阻力，計算方法分別為二因次法與三因次法。

二因次法具體步驟如下：

(1)由CFD模擬求得船模的總阻力系數Ctm

(8)

式中：Rtm為船模阻力；ρ為淡水密度；Sm為船模濕表面積；Vm為船模速度。

(2)利用桑海公式計算船模和實船摩擦阻力系數Cf

(9)

式中：Re為雷諾數。

(3)求出船模的剩余阻力系數Crm

Crm=Ctm-Cfm

(10)

式中：Ctm、Cfm分別為船模阻力系數、船模摩擦阻力系數。

(4)計算實船總阻力系數Cts

Cts=Crs+Cfs+Ca+ΔCf

(11)

式中：Crs為根據弗勞德數假設Crs=Crm求得的實船剩余阻力系數；Cfs為實船摩擦阻力系數；Ca為空氣阻力系數，壓載狀態下取0.142×10-3；ΔCf為粗糙度補貼系數，壓載狀態下取值為0.046×10-3。

(5)計算實船阻力公式

(12)

式中：Rts為實船阻力；Ss為實船濕表面面積；Vs為實船速度。

由于三因次法認為粘壓阻力系數Cpv與摩擦阻力系數Cf之比為常數ξ，即

(13)

實船總阻力系數為：

Cts=(1+k)Cfs+Cws+ΔCf

(14)

式中：Cws為實船興波阻力系數。

船模總阻力系數為：

Ctm=(1+k)Cfm+Cwm

(15)

式中：Cwm為船模興波阻力系數；k為形狀系數，(1+k)為形狀因子，僅與船體形狀有關。

又因：

Cwm=Cws

故實船總阻力系數為：

Cts=(1+k)(Cfs-Cfm)+Ctm+ΔCf

(16)

式中(1+k)形狀因子可通過普魯哈斯卡的假設求得。普魯哈斯卡認為弗勞德數在0.1～0.2之間時，可近似認為興波阻力Cwm與弗勞德數Fr的4次方成正比關系。故有：

(17)

通過線性擬合可求得(1+k)值。

摩擦阻力系數可用1957ITTC公式計算：

(18)

(19)

式中：ks為船體表面粗糙度，通常取1.50×10-4m。

因此，實船的總阻力為：

(20)

由以上公式將船模換算成實船數據后再與試驗值進行對比，結果見表2。表中結果證明，該數值方法與船模數據換算方法能較為準確地預報船舶阻力性能。

表2 KCS船模阻力計算結果

2 計算模型

本文目標船型為60 m鎮揚車客渡船模型，選用Solidworks軟件完成三維建模，縮尺后導入Hexpress軟件里進行網格劃分，縮尺比為1∶16。在船體附近及自由液面進行特殊加密，其中：上部分為空氣域，下部分為水域。計算船舶阻力時，建立計算域的長、寬、高分別為15、11.25、9.375 m。船模距前方入口1LOA(LOA為總長)，距出口3LOA，距上下邊界分別0.5LOA和2LOA，左右邊界各1.5LOA。計算域內船體采用標準壁面函數，上下邊界設置為指定壓強(液體靜壓)，入口邊界、出口邊界及兩側邊界設定為遠場邊界。圖4為船體模型三視圖，圖5為表面網格劃分示意圖，實船與船模主要參數見表3。

圖4 船體三視圖

圖5 計算域及網格剖面

表3 車客渡船主要參數

本次計算5種吃水狀態、實船航速Vs=2～12 kn共30種工況，見表4。湍流模型選用k-omega(SST-Menter)模型，k與ω的值與雷諾數有關，采用中心差分(AVLSMART)格式離散；自由面采用混合自由面捕捉與重構相結合(BRICS)格式離散；采用壓力速度耦合算法求解。

3 結果及分析

3.1 車客渡船阻力及有效功率的計算與分析

本文在預報了船舶分別在A航區滿載到港、滿載出港、空載到港，以及B航區滿載出港、滿載到港的情況下，以實船航速2、4、6、8、12 kn和設計航速10 kn共30種工況下的運動狀態，船模模擬速度Vm分別為0.257、0.514、0.772、1.029、1.286、1.543 m/s。當流場穩定，阻力變化穩定后，取穩定時間內數據的平均值。A航區、B航區所得船模阻力數值計算結果分別見表5、表6。

表4 計算工況

表5 A航區不同狀態船模阻力系數表

本文所研究的車客渡船設計航速為10 kn，弗勞德數為0.210，屬于低速船。而對于低速船來說，粘性阻力所占的比例往往比較大(可達70%以上)；興波阻力所占的比例則比較小。此船的數據也正符合這一規律，并隨著速度的增加，興波阻力所占比例逐漸變大，總阻力值增長速度也增大。

表6 B航區不同狀態船模阻力系數表

對該船采用二因次法換算其實船數據，并與三因次法換算所得數據進行比較，繪出阻力隨航速變化關系曲線，見圖6。

圖6 三因次法和二因次法的計算結果對比

從圖6中可以看出，隨著速度的增加，阻力在逐漸增大。實船阻力換算方法中二因次法小于三因次法。首先，二因次法把船體阻力分成了互不相關的兩個獨立部分：一部分僅與重力或弗勞德數有關；另一部分僅與粘性或雷諾數有關，而忽略了兩者的相互影響，事實上這兩種影響是存在的。其次，弗勞德數將興波阻力和粘壓阻力這兩種不同性質的阻力合并為剩余阻力，在理論上是不恰當的。最后，船體形狀是相對復雜的三因次物體，其周圍流動情況與平板相比顯然有一定差別，因而用相當平板的摩擦阻力來代替船體摩擦阻力，必然是有誤差的。因此，三因次法計算阻力更加符合實際情況。本文采用三因次法給出船舶不同工況下的阻力，見圖7。船舶附體阻力取船舶阻力的10%，得到總阻力后乘以對應航速得到船舶的有效功率，見圖8。

圖7 車客渡船阻力特性曲線

從圖7、圖8中可以看出：

(1)船舶阻力隨著航速和排水量的增加而增加，呈二次型增加趨勢。

(2)空載到港阻力最小，設計航速下的阻力約為67.6 kN，有效功率約為347 kW。

(3)A航區航行中滿載出港阻力最大，設計航速下的阻力約為85 kN，有效功率約為450 kW。

圖8 車客渡船有效功率特性曲線

3.2 車客渡船流場分析

圖9給出了不同航速下車客渡船自由面波形圖。隨著航速增加船體興波愈加顯著，與常規單體船開爾文波形不同。由于車客渡船下潛體呈梭型，具有分流作用，在船體尾流場重新匯聚，尾波形較為分散且不斷向外擴散。

車客渡船船體表面壓力分布云圖見圖10。從圖中看出，船底梭型潛體艏部為壓力駐點形成高壓區，對流體形成分流作用，降低了高壓區的面積，流體繞過梭型最寬區域后，梭型體后半部分由于流體分離形成渦流區域，出現一個低壓區，所以整體梭型結構相對于方形結構具有明顯的導流減阻的效果。

圖9 車客渡船自由面波形圖

圖10 車客渡船船體表面壓力分布云圖

4 結論

本文基于粘性流體理論，完成了梭型車客渡船不同工況下的阻力數值預報，采用三因次法將船模阻力換算至實船阻力，并對車客渡船尾流場進行分析，得到以下結論：

(1)本文所建立數值預報方法可較為準確地預報船舶阻力及流場分布規律。

(2)三因次法船模阻力換算至實船阻力考慮因素更加貼合實際情況，準確性較高。

(3)車客渡下潛體采用梭型結構相對于常規方形底，具有分流減阻效果。