王志東，凌宏杰

(江蘇科技大學(xué)船舶與海洋工程學(xué)院，江蘇鎮(zhèn)江212003)

滑行艇的流體動力特性研究早在上世紀(jì)60年代就引起了眾多學(xué)者的關(guān)注，并且在滑行艇阻力性能預(yù)報(bào)方面取得了顯著的研究成果，提出了一系列滑行艇阻力計(jì)算的經(jīng)驗(yàn)或半經(jīng)驗(yàn)公式［1-4］.近年來以實(shí)時求解雷諾時均方程RANS為目標(biāo)的現(xiàn)代CFD技術(shù)被應(yīng)用到水面高速艇流體動力性能的精確預(yù)報(bào)，文獻(xiàn)［5］基于有限體積方法研究了高速艇的運(yùn)動特性，其中壓力與速度的耦合采用分裂算法，水氣交界面的模擬采用VOF方法，采用貼體網(wǎng)格系統(tǒng)數(shù)值計(jì)算了定常直航時流體作用在高速艇上的力和力矩以及6自由度運(yùn)動響應(yīng)，給出了高速艇的運(yùn)動速度、加速度及其位置，通過與試驗(yàn)對比表明，基于數(shù)值求解RANS方程的數(shù)值方法可用于指導(dǎo)高速艇的流體動力設(shè)計(jì)及運(yùn)動性能預(yù)報(bào).文獻(xiàn)［6-8］利用非定常RANS方法編制了水面艦艇6自由度運(yùn)動響應(yīng)預(yù)報(bào)程序CFDSHIP-IOWA，求解器采用高階上風(fēng)離散格式、PISO壓力速度耦合方法k-ε/k-ω兩方程湍流模型，其中自由表面追蹤技術(shù)分別采用VOF方法和Level Set方法，數(shù)值預(yù)報(bào)了水面艦艇在多種航行狀態(tài)下的自由面繞流場及其運(yùn)動響應(yīng)特性，結(jié)果優(yōu)于傳統(tǒng)的線性切片理論，同時指出水面波形破碎是引起船舶非線性流體動力載荷的主要原因.

邊界元方法［9］也被用于滑行艇模型的水動力性能數(shù)值預(yù)報(bào)，通過開展滑行平板、楔形體、底部斜升體等3種滑行艇模型的計(jì)算分析，探討了不同航速下的壓力分布、阻力、升力及波形，與相關(guān)試驗(yàn)結(jié)果的吻合度較好.

文獻(xiàn)［10］利用計(jì)算機(jī)流體力學(xué)軟件FLUENT數(shù)值模擬了滑行艇在靜水中的直航運(yùn)動，給出了滑行艇阻力隨航行速度的變化曲線，探討了滑行艇底部的壓力分布及尾流場;文獻(xiàn)［11］建立了高速滑行艇6自由度運(yùn)動方程，利用Matlab軟件開展了滑行艇的運(yùn)動與姿態(tài)控制的仿真預(yù)報(bào)研究.

文中基于C語言編寫了滑行艇升沉縱搖耦合運(yùn)動預(yù)報(bào)用戶自定義程序，采用預(yù)測/修正法求解運(yùn)動方程，將求解方程嵌入到FLUENT軟件中，對5 種不同體積傅汝德數(shù) Fr2=0.9，1.8，2.7，3.6，5.0條件下滑行艇的縱搖與升沉耦合運(yùn)動特性進(jìn)行研究.在計(jì)算過程中為保證當(dāng)艇體產(chǎn)生大位移運(yùn)動時的網(wǎng)格質(zhì)量，采用彈簧法和局部重構(gòu)聯(lián)合作用的運(yùn)動網(wǎng)格技術(shù)，控制艇體表面及外部網(wǎng)格的生長速度和尺度;同時為保證計(jì)算精度，對網(wǎng)格重構(gòu)的時間步長進(jìn)行了合理地控制.

1 計(jì)算模型

1.1 計(jì)算模型和網(wǎng)格劃分

針對某滑行艇艇型模型，采用Maxsurf軟件進(jìn)行三維建模，導(dǎo)入Gambit軟件進(jìn)行網(wǎng)格劃分.計(jì)算模型的網(wǎng)格劃分采用分區(qū)混合網(wǎng)格系統(tǒng)，并在艇體附近進(jìn)行網(wǎng)格加密，其中艇體附近的球形和正方體網(wǎng)格區(qū)域用于保證艇體縱搖、升沉運(yùn)動時的網(wǎng)格質(zhì)量.圖1為計(jì)算模型及網(wǎng)格劃分示意圖，其中上部分為空氣域，下部分為水域.水域?yàn)?0L×4L×3L，空氣域?yàn)?0L×4L×1.5L(L為艇長).計(jì)算區(qū)域底部為固壁無滑移邊界條件，左側(cè)、上邊界與兩側(cè)壁為速度入口，右邊界為壓力出口.滑行艇模型的主要參數(shù)見表1.

圖1 三維模型計(jì)算控制區(qū)域及網(wǎng)格劃分Fig.1 Calculation control zone and grid partition of the three-dimensional model

表1 滑行艇主尺度Table 1 Principal dimensions of planning boat

1.2 數(shù)值計(jì)算方法

采用三維非定常分離隱式求解器，利用VOF方法追蹤自由液面，自由面重構(gòu)格式采用Geo-Reconstruct格式，選用k-ε湍流模型，控制方程的擴(kuò)散項(xiàng)采用中心差分格式離散，對流項(xiàng)采用二階迎風(fēng)格式;壓力方程采用Body Force Weighted格式，壓力速度合求解采用PISO算法，動網(wǎng)格采用彈性光順法和局部重構(gòu)，彈性常數(shù)參數(shù)設(shè)為默認(rèn)值，局部重構(gòu)算法中的最大網(wǎng)格扭曲率(maxim-um cell skewness)設(shè)為 0.85，Size Remesh Interval設(shè)為 5;啟動Size Function函數(shù)，調(diào)整尺寸函數(shù)變化率β(size function rate)為1.5，控制邊界處網(wǎng)格和內(nèi)部網(wǎng)格的生長速度.

1.3 滑行艇運(yùn)動預(yù)報(bào)算法

不計(jì)滑行艇的橫蕩、艏搖及橫搖，在滑行艇縱向三自由度的基礎(chǔ)上添加人工阻尼項(xiàng)，得到滑行艇的升沉、縱搖運(yùn)動控制方程［7］如下:

式中:m為艇體的質(zhì)量(kg);ω為艇體垂向速度(m·s-1);u為滑行艇前進(jìn)速度(m·s-1);Z為艇體垂向合力(N);θ為縱搖角(rad);My為艇體縱搖力矩;Iy為艇體對y軸的轉(zhuǎn)動慣量;Dh為垂蕩人工阻尼系數(shù)，取0.1;Dp為縱搖人工阻尼系數(shù)，取0.1.

采用預(yù)測/修正的隱式方法求解運(yùn)動方程，并通過求解滑行艇的運(yùn)動控制方程預(yù)測線速度和角速度，然后采用三階迎風(fēng)格式進(jìn)行修正:

判斷計(jì)算結(jié)果是否收斂.若收斂則將艇體運(yùn)動賦予網(wǎng)格系統(tǒng)［7］，滑行艇運(yùn)動響應(yīng)預(yù)報(bào)算法流程如圖2.

圖2 滑行艇運(yùn)動響應(yīng)預(yù)報(bào)算法流程Fig.2 Flowchart for the respone prediction of planing

2 計(jì)算結(jié)果與分析

從水動力學(xué)的觀點(diǎn)出發(fā)，若船舶的體積傅汝德數(shù)Fr2≥1.0，則屬于高速船，其中包括高速排水型船舶和動力增升型船舶.對于1.0＜Fr2＜3.0航速區(qū)間的高速排水型船舶，在其支持力中起主要作用的是靜浮力;而對于Fr2≥3.0的流體動力增升型船舶來說，其支持力起主要作用的是流體動升力.高速的流體作用于艇體表面，在水氣交界面處產(chǎn)生劇烈的噴濺現(xiàn)象，噴濺對航行中高速船的影響不容忽視.

2.1 不同F(xiàn)r2數(shù)下滑行艇的水動力性能

為了分析不同航速下滑行艇的水動力性能，選取與模型試驗(yàn)相同的計(jì)算工況，航速V=2，4，6，8，11 m·s-1，對應(yīng)的體積傅汝德數(shù) Fr2分別為 0.9，1.8，2.7，3.6，5.0.當(dāng) V=11 m·s-1時，滑行艇出現(xiàn)“海豚運(yùn)動”，在此不作深入研究.

圖3，4給出了滑行艇在不同航速下滑行艇阻力計(jì)算值與試驗(yàn)阻力的對比曲線和誤差曲線，結(jié)果表明:在低航速時誤差較小(10%左右)，航速較高時誤差較大(20%左右).出現(xiàn)這種現(xiàn)象的原因主要有:低航速時滑行艇處于排水航行狀態(tài)，總阻力主要由摩擦阻力貢獻(xiàn)，壓差阻力貢獻(xiàn)較小;高航速時情況恰好相反.由于滑行艇模型的轉(zhuǎn)動慣量不能準(zhǔn)確給定，引起滑行艇穩(wěn)定時的縱傾角偏大，縱傾角的大小直接影響壓差阻力，因此出現(xiàn)計(jì)算值大于試驗(yàn)值的現(xiàn)象.同時，由于計(jì)算機(jī)硬件條件的限制，艇體周圍的網(wǎng)格劃分不夠密，單位網(wǎng)格覆蓋艇體的面積過大，影響計(jì)算精度與準(zhǔn)確性.

航速 V=2，4，6，8 m·s-1，4 種工況下滑行艇能夠達(dá)到“動態(tài)平衡”狀態(tài).

圖3 艇體阻力曲線Fig.3 Resisitance curves of vessel

圖4 滑行艇阻力誤差曲線Fig.4 Resisitance error curves of vessel

圖5給出了滑行艇在不同航速下的力矩變化曲線，低速航行時滑行艇處于排水航行狀態(tài)，滑行艇縱傾角較小，艇體所受力矩較小;隨著航速提升滑行艇由排水航行向滑行艇狀態(tài)過渡，縱傾角變大，艇體所受力矩變大;滑行艇進(jìn)入滑行狀態(tài)后，隨著航速進(jìn)一步提升，滑行艇的縱傾角下降，艇體所受力矩出現(xiàn)下降.

圖5 滑行艇力矩曲線Fig.5 Moment curves of vessel

2.2 滑行艇的升沉縱搖運(yùn)動特性分析

圖6給出了滑行艇在不同航速下的滑行姿態(tài)的數(shù)值模擬結(jié)果.結(jié)果表明，滑行艇在航行過程中，在給定的主機(jī)輸入功率下滑行艇將產(chǎn)生初始航速，并出現(xiàn)初始攻角，產(chǎn)生的升力將使艇體抬升到一定的高度，同時升力對艇體重心產(chǎn)生俯仰力矩，從而使滑行艇原來的平衡狀態(tài)被打破.為達(dá)到新的動態(tài)平衡，艇體在不斷變化的垂向力和力矩同時作用下產(chǎn)生升沉和縱搖運(yùn)動，進(jìn)而引起艇體的排水體積與浸濕面積同時發(fā)生變化，改變了滑行艇的航行阻力，從而使滑行艇的航速發(fā)生變化，而航速的變化又會引起滑行艇升力及力矩的變化.當(dāng)滑行艇所受到的垂向合力與合力矩為零時，滑行艇的航行狀態(tài)達(dá)到了動態(tài)平衡［6］，滑行艇將保持這一姿態(tài)穩(wěn)定航行，直到再次受到外界的干擾.文中的數(shù)值計(jì)算模型中，由于滑行艇的航行速度為定值，所以滑行艇最終趨于動態(tài)平衡.

圖6 滑行艇不同航速下的滑行姿態(tài)Fig.6 Planning attribute of vessel at different speeds

圖7，8 給出了航速 V=2，4，6，8，11 m·s-1時艇體重心位置的變化曲線和縱傾角變化曲線，反映了滑行艇升沉和縱搖運(yùn)動的變化規(guī)律.從圖7a)，8a)可以看出，當(dāng) V=11 m·s-1時，艇體垂蕩速度和縱搖角速度呈周期性變化，因此艇體升沉量和縱傾角也呈周期性變化，出現(xiàn)“海豚運(yùn)動”現(xiàn)象，無法達(dá)到動態(tài)平衡.其他工況下，艇體垂蕩速度和縱搖角速度趨于零，艇體升沉量和縱傾角趨于穩(wěn)定值，趨于動態(tài)平衡.當(dāng)滑行艇低航速時，出現(xiàn)下沉現(xiàn)象，當(dāng)航速提高時，艇體抬升進(jìn)入滑行狀態(tài).

圖7 不同航速下滑行艇升沉變化特性Fig.7 Characteristics of heave at different speeds

圖8 不同航速下滑行艇縱搖變化特性Fig.8 Characteristics of pitch at different speeds

圖9，10給出了不同航速下滑行艇航行達(dá)到穩(wěn)定時升沉量和縱傾角計(jì)算值與試驗(yàn)值的對比曲線.升沉量與縱傾角的總體走勢計(jì)算值與試驗(yàn)值一致;當(dāng)V=2，8 m/s時誤差較大，升沉量誤差20%左右;當(dāng) V=4，6 m/s時誤差較小，升沉量誤差8.5%左右.滑行艇模型的轉(zhuǎn)動慣量Iy直接影響艇體縱傾角，由于無法準(zhǔn)確給定滑行艇模型的轉(zhuǎn)動慣量，因此縱傾角誤差較大，上文中提到的滑行艇總阻力隨航速增加，計(jì)算誤差增大的原因也在于此.滑行艇高速航行時，縱傾角對壓差阻力的影響較大.

圖9 升沉量曲線Fig.9 Curves of heave

圖10 縱傾角曲線Fig.10 Curves of pitch

2.3 滑行艇艇底動壓分布規(guī)律

圖11 滑行艇艇底動壓Fig.11 Dynamic pressure at bottom of vessel

2.4 滑行艇尾流場波高變化規(guī)律

均勻來流繞過滑行艇，由于滑行艇的存在流場發(fā)生相應(yīng)的變化，圖12給出了不同航速下滑行艇航行達(dá)到穩(wěn)定時尾流場的波高云圖.從圖中可以看出，隨著航速不斷增加船體艏部抬出水面，水流從兩側(cè)流出，產(chǎn)生大量的飛濺和水花，尾流場波高不斷增加，艇體尾緣“空穴”變大，形成“雞尾流”，波峰后移.

3 結(jié)論

基于計(jì)算流體力學(xué)軟件FLUENT，完成了滑行艇升沉縱搖運(yùn)動的實(shí)時數(shù)值預(yù)報(bào)研究，主要得出如下結(jié)論:

1)轉(zhuǎn)動慣量Iy對滑行艇運(yùn)動響應(yīng)數(shù)值預(yù)報(bào)的準(zhǔn)確性具有至關(guān)重要的影響，滑行艇轉(zhuǎn)動慣量Iy計(jì)算的經(jīng)驗(yàn)公式有待修正;

2)基于計(jì)算流體力學(xué)軟件FLUENT及運(yùn)動網(wǎng)絡(luò)技術(shù)，耦合求解滑行艇縱向運(yùn)動預(yù)報(bào)程序，能夠預(yù)報(bào)滑行艇的“海豚運(yùn)動”現(xiàn)象.

艇底滑行面上的動壓是艇體產(chǎn)生動升力的主要原因，圖11給出了不同航速下滑行艇航行達(dá)到穩(wěn)定狀態(tài)時，艇底部動壓分布云圖.從圖中可以看出，隨著航速的增大艇底部動壓值不斷增大且動壓力中心位置后移.

圖12 滑行艇尾流場波高Fig.12 High of weak wave of vessel

References)

［1］Savitsky D，Brown P W.Procedures for hydrodynamic evaluation of planning hulls in smooth and rough water［J］.Marine Technology，1976，13(4):381 -400.

［2］John M A.Resistance prediction of planing hulls:state of the art［J］.Marine Technology，1993，30(4):296 -307.

［3］Nagai T，Yoshida Y.Estimation of resistance trim and draft of planing craft［J］.Ship Technology Research，1993，40(3):133-137.

［4］邵世明，王云才.滑行艇阻力試驗(yàn)［R］.上海:上海交通大學(xué)，1980.

［5］ Panahi R，Jahanbakhsh E.Toward simulation of 3D nonlinear high-speed vessels motion［J］.Ocean Engineering，2009，36:256 -265.

［6］Wilson R V，Carrica P M，Stern F.Unsteady RANS method for ship motions with application to roll for a surface combatant［J］.Computers ＆ Fluids，2006，35:501-524.

［7］Carrica P M，Wilson R V，Stern F.Unsteady RANS simulation of the ship forward speed diffraction problem［J］.Computers ＆ Fluids，2006，35:545 -570.

［8］Carrica P M，Wilson R V，Noack R W.Ship motions using single-phase level set with dynamic overset grids［J］.Computers ＆ Fluids，2007，36:1415 -1433.

［9］Ghassemi H，Kohansal A R.Higher order boundary element method applied to the hydrofoil beneath the free surface［C］∥Proceedings of the 28th International Conference on O-cean，Offshore and Engineering(OMAE2009).Honolulu，Hawaii USA:［s.n.］，2009.

［10］王兆立，牛江龍，秦再白，等.基于CFD理論的滑行艇阻力數(shù)值計(jì)算［C］∥第十四屆中國海洋(岸)工程學(xué)術(shù)討論會論文集.內(nèi)蒙古:［s.n.］，2009，8:309 -315.

［11］高雙，朱齊丹，李磊.基于神經(jīng)網(wǎng)絡(luò)的高速無人艇模糊PID 控制［J］.系統(tǒng)仿真學(xué)報(bào)，2007，19(4):776 -779.