基于cfd方法的moeri集裝箱船自航性能數(shù)值模擬

基于cfd方法的moeri集裝箱船自航性能數(shù)值模擬

目前，研究船槳干擾的數(shù)值方法可分為兩類：一般建模和體積力法。在使用船槳通用建模方法研究船槳干擾問題時，應(yīng)全面報告信息，以反映真流場中螺旋作用區(qū)域的復(fù)雜性，因此使用該方法是廣泛的，但計(jì)算往往需要很長時間。體積力法較整體建模方法而言,建模簡單,計(jì)算效率高,能預(yù)報流場的基本信息,因此得到一定的應(yīng)用.體積力源項(xiàng)通常可由基于勢流理論的螺旋槳模型計(jì)算得到,依據(jù)它與雷諾平均納維-斯托克斯(RANS)方程的耦合方式,體積力法可以分為迭代型和描述型兩種.迭代型體積力法有較多應(yīng)用,但該方法計(jì)算過程較為繁鎖,需要在螺旋槳性能計(jì)算和RANS方程計(jì)算之間進(jìn)行相互迭代.描述型體積力法則相對簡便,對實(shí)際螺旋槳的幾何參數(shù)要求十分簡單,且不需要螺旋槳性能計(jì)算,但該方法無法直接得到螺旋槳性能參數(shù)而須采用其他方法.本文通過計(jì)算流體力學(xué)(CFD)軟件Fluent的RANS求解器,采用文獻(xiàn)中的描述型體積力模型,分別在無舵實(shí)船自航點(diǎn)和有舵模型自航點(diǎn)2種條件下對MOERI(KoreaResearchInstituteforShipsandOceanEngineering)集裝箱船(KCS)的自航性能進(jìn)行了數(shù)值預(yù)報,并通過等推力法原理、敞水計(jì)算和相對應(yīng)的螺旋槳敞水曲線獲得船體自航性能相關(guān)參數(shù),通過與試驗(yàn)結(jié)果的對比驗(yàn)證了該方法的有效性.1數(shù)學(xué)模型1.1均速度分量a流場的控制方程為連續(xù)方程和動量方程:?ui?xi=0(1)?ui?t+uj?ui?xj=-??p?xi+1Re?2ui?xj?xj-??xjˉu′iu′j+f*i(2)?ui?xi=0(1)?ui?t+uj?ui?xj=??p??xi+1Re?2ui?xj?xj???xju′iu′jˉˉˉˉˉˉˉˉ+f?i(2)式中:ui(i=1,2,3)為雷諾時均速度矢量的分量;?pp?為無量綱系數(shù);ˉu′iu′ju′iu′jˉˉˉˉˉˉˉˉ為雷諾應(yīng)力;f*i為無量綱體積力,代表螺旋槳力源項(xiàng).1.2槳葉的推力及轉(zhuǎn)速本文采用文獻(xiàn)中的方法施加體積力源項(xiàng),該方法通過對螺旋槳所在區(qū)域施加推力和轉(zhuǎn)矩代替槳葉表面載荷,公式如下:f*b?x=Axr*√1-r*(3)f*b?θ=Aθr*√1-r*(1-RΗ)r*+YΗ(4)f?b?x=Axr?1?r??????√(3)f?b?θ=Aθr?1?r?√(1?RH)r?+YH(4)式中:r*=r-RΗRΡ-RΗ(5)Ax=CtΔx10516(4+3YΗ)(1-YΗ)(6)Aθ=ΚΜJ2Δx105π(4+3YΗ)(1-YΗ)(7)r?=r?RHRP?RH(5)Ax=CtΔx10516(4+3YH)(1?YH)(6)Aθ=KMJ2Δx105π(4+3YH)(1?YH)(7)f*b,x和f*b,θ分別為推力和轉(zhuǎn)矩在螺旋槳區(qū)域內(nèi)的分布;r為螺旋槳區(qū)域內(nèi)任意一點(diǎn)至螺旋槳軸線的距離;YH=RH/RP為轂徑比,RP、RH分別為螺旋槳半徑和槳轂半徑;Δx為螺旋槳在軸向的投影弦長,一般取槳轂厚度δ;Ct=4T/(ρu2πDP)為一種無量綱推力系數(shù),T為螺旋槳推力,ρ為流體密度,u為船速,DP為螺旋槳直徑;KM=M/(ρn2D5P)為轉(zhuǎn)矩系數(shù),M為轉(zhuǎn)矩,n為轉(zhuǎn)速;J=va/(nDP)為螺旋槳進(jìn)速系數(shù),va為螺旋槳進(jìn)速.在實(shí)際計(jì)算中,施加體積力的具體過程是將推力和轉(zhuǎn)矩施加給螺旋槳作用區(qū)域的網(wǎng)格單元進(jìn)而對控制方程進(jìn)行求解計(jì)算.由于船艉形狀復(fù)雜,網(wǎng)格單元與螺旋槳區(qū)域邊界不易重合導(dǎo)致實(shí)際體積力大小與指定值相比會有一定誤差.為此,本文使用UDF(UserDefinedFunction)施加體積力時,在UDF源代碼中添加計(jì)算推力的DEFINE_ADJUST宏和修正項(xiàng),用以查看實(shí)際推力大小并進(jìn)行修正以達(dá)到指定值.需要說明的是,修正后推力T和轉(zhuǎn)矩M的變化率相同,這表明通過修正項(xiàng)對推力進(jìn)行修正時,也能對轉(zhuǎn)矩進(jìn)行精確的修正,進(jìn)而達(dá)到準(zhǔn)確地施加體積力的目的.2計(jì)算條件2.1設(shè)計(jì)水線和船模的點(diǎn)坐標(biāo)本文的研究對象為KCS集裝箱船模型.該船型關(guān)于實(shí)船自航點(diǎn)和模型自航點(diǎn)的試驗(yàn)數(shù)據(jù)取自文獻(xiàn).計(jì)算中參考坐標(biāo)系的坐標(biāo)原點(diǎn)位于船艏中心線與設(shè)計(jì)水線的交點(diǎn),x軸為中縱剖面和設(shè)計(jì)水線面的交線,正方向指向船艉;y軸正方向指向船體右弦;z軸正方向向上.本文共研究了2個算例.算例1為靜水中實(shí)船自航點(diǎn)的自航試驗(yàn),船體附體只包括槳;算例2則是靜水中船模自航點(diǎn)的自航試驗(yàn),其中船?？勺杂缮梁涂v傾,船體附體包括槳和舵.以上2個算例的具體模擬計(jì)算條件如表1所示.表中:LPP為船體垂線間長;Fo為傅里葉數(shù);Re為雷諾數(shù);θ為溫度;Sw為濕表面積.2.2網(wǎng)格劃分參數(shù)本文的主要內(nèi)容是預(yù)報船體自航性能.為提高計(jì)算精度和收斂性,采用分塊結(jié)構(gòu)化網(wǎng)格,使用CFD網(wǎng)格生成軟件Gambit2.4生成網(wǎng)格,其中2個算例的計(jì)算域均為長方形.具體的網(wǎng)格劃分參數(shù)如表2所示.表中:y+為壁單位;Δyp為最小網(wǎng)格間距.網(wǎng)格劃分的難點(diǎn)在于船艉區(qū)域,2個算例的船艉區(qū)域網(wǎng)格劃分如圖1所示.3近壁流動動力學(xué)方程本節(jié)采用商用RANS求解器FLUENT6.3中的VOF方法在船模尺度下對自由液面情況的三維黏性流場進(jìn)行數(shù)值求解,即用有限體積法離散RANS方程,并采用壓力基求解器進(jìn)行求解計(jì)算.計(jì)算中PISO(PressureImplicitwithSplittingofOperators)算法和SIMPLE(Semi-ImplicitforPressure-LinkedEquations)算法先后被用于速度和壓力間的耦合計(jì)算.湍流方程采用帶旋流修正k-ε湍流模式封閉方程,并使用非平衡壁面函數(shù)簡化處理近壁流動.數(shù)值離散方面,動量方程中的對流項(xiàng)采用三階QUICK(QuadraticUpwindInterpolationofConvectiveKinematics)格式,k-ε方程和體積分?jǐn)?shù)方程均采用二階上風(fēng)格式.計(jì)算域的邊界條件設(shè)置為:進(jìn)流面,表面壓力為0;上表面,滑移邊界條件;船體表面,無滑移邊界條件;側(cè)表面,滑移邊界條件;下表面,滑移邊界條件;出流面,表面壓力為0.在非定常計(jì)算過程中,取時間步長Δt=0.5ms,最大內(nèi)循環(huán)步數(shù)為10步.3.1試驗(yàn)數(shù)據(jù)分析在應(yīng)用體積力法預(yù)報船體自航性能之前需先計(jì)算出裸船體的總阻力系數(shù)CT、摩擦阻力系數(shù)CF、興波阻力系數(shù)CP及標(biāo)稱伴流系數(shù)φ.如表3所示,阻力系數(shù)和標(biāo)稱伴流系數(shù)的計(jì)算值與試驗(yàn)數(shù)據(jù)吻合良好,相對偏差均在5%以內(nèi).3.2預(yù)報忽略轉(zhuǎn)速的船體自航性能為獲得船體自航性能參數(shù),本文的計(jì)算流程分為2步:(1)預(yù)報忽略轉(zhuǎn)矩的船體自航性能,體積力中僅包含推力T;(2)預(yù)報考慮轉(zhuǎn)矩影響的船體自航性能,體積力中同時包含推力T和轉(zhuǎn)矩M.3.2.1旋轉(zhuǎn)進(jìn)速系數(shù)(1)依據(jù)等推力法原理求取平衡點(diǎn)推力TB.由于模擬計(jì)算的是實(shí)船自航點(diǎn)的自航試驗(yàn),考慮到摩擦阻力修正,平衡點(diǎn)推力的求取依據(jù)下式:Τe=Τ-RΤ,SΡ+RF(8)Te=T?RT,SP+RF(8)式中:T為所加入的螺旋槳推力;RT,SP為加入推力T后計(jì)算得到船體自航總阻力;RF為實(shí)船和船模間的摩擦阻力修正值,本文取RF=30.25N,與該自航試驗(yàn)中的取值相同;Te=0所對應(yīng)的T即為平衡點(diǎn)推力TB.實(shí)際求取Te時是根據(jù)不同的推力T′、T″得出不同的T′e、T″e,這兩者不一定等于零,但要保證異號并且其絕對值應(yīng)盡可能小.然后通過T′、T″線性插值得出平衡點(diǎn)推力TB.本文中的計(jì)算,先后嘗試了T=59,57,56,55N才最終確定當(dāng)T′=55N,T′e≈-0.646N以及T″=56N,T″e≈0.438N時,兩者異號,再通過線性插值得出平衡點(diǎn)推力TB=55.6N,如圖2所示.(2)求解螺旋槳進(jìn)速va.其原理是先將僅包含平衡點(diǎn)推力TB=55.6N的體積力加入到3.1節(jié)得到的計(jì)算結(jié)果中再次計(jì)算,待計(jì)算收斂導(dǎo)出船后槳盤面處平均速度的平方,即ˉv2diskvˉ2disk.本文得出的ˉv2disk=3.6636m2/s2vˉ2disk=3.6636m2/s2,相應(yīng)的RT,SP=86.376N,故推力減額系數(shù)Δk=(RT,SP-RT,Tow)/TB=0.870.其中,RT,Tow為3.1節(jié)中裸船體總阻力.為了得到va,運(yùn)用動量守恒原理進(jìn)行敞水計(jì)算,即無自由液面的情況下,賦予螺旋槳體積力為平衡點(diǎn)推力TB,改變來流速度v,使得槳盤面處ˉv2diskvˉ2disk與上文中ˉv2diskvˉ2disk相等,則此時的來流速度v即為螺旋槳進(jìn)速va.通過敞水計(jì)算得出ˉv2disk=3.6634m2/s2vˉ2disk=3.6634m2/s2對應(yīng)的進(jìn)速va=1.7617m/s.(3)求解螺旋槳進(jìn)速系數(shù).通過下式得到K′T與J的二次曲線,Κ′Τ=ΤBJ2ρv2aD2Ρ(9)K′T=TBJ2ρv2aD2P(9)式中:ρ=999.01kg/m3;va=1.7617m/s;DP=0.25m;TB=55.6N.然后將K′T的曲線與螺旋槳敞水曲線KT繪制在同一圖中并取兩者交點(diǎn),即可求得進(jìn)速系數(shù)JB.具體做法如圖3所示,其中進(jìn)速系數(shù)JB=0.749.(4)通過下式求解忽略轉(zhuǎn)矩的船體自航性能:推力系數(shù)ΚΤ=Τ/(ρn2D4Ρ)(10)KT=T/(ρn2D4P)(10)轉(zhuǎn)矩系數(shù)ΚΜ=Μ/(ρn2D5Ρ)(11)KM=M/(ρn2D5P)(11)螺旋槳轉(zhuǎn)數(shù)n=va/(JDΡ)(12)n=va/(JDP)(12)敞水效率ηΟ=ΚΤJ/(2πΚΜ)(13)船身效率ηΗ=Δk/(1-φ)(14)推進(jìn)效率η=ηΟηRηΗ(15)式中,ηR為相對旋轉(zhuǎn)效率,取ηR=1.通過上述公式計(jì)算得到的忽略轉(zhuǎn)矩影響的船體自航性能相關(guān)參數(shù)如表4所示.表中,e為誤差.3.2.2基于扭矩的自航性能計(jì)算若要考慮轉(zhuǎn)矩的影響,則需要將同時包含表5中的轉(zhuǎn)矩M和平衡點(diǎn)推力TB=55.6N的體積力再一次加入到3.1節(jié)所得到的裸船體計(jì)算結(jié)果中進(jìn)行計(jì)算,然后按照上節(jié)的步驟依次求取螺旋槳進(jìn)速va和進(jìn)速系數(shù)J,再通過式(10)～(15)進(jìn)行計(jì)算便可得到考慮轉(zhuǎn)矩影響的船體自航性能相關(guān)參數(shù).在考慮轉(zhuǎn)矩情況下,經(jīng)計(jì)算可分別得到進(jìn)速va=1.768m/s和進(jìn)速系數(shù)J=0.748(見圖4),其他自航性能參數(shù)如表5所示.圖5、6為槳盤面后方0.25DP處軸向速度分布的計(jì)算值和試驗(yàn)值.由圖5,6可見,考慮轉(zhuǎn)矩影響時,槳盤面軸向速度分布體現(xiàn)了螺旋槳的旋轉(zhuǎn)效應(yīng),但與試驗(yàn)結(jié)果仍存在一定差異.對比表4、5可知,考慮轉(zhuǎn)矩較忽略轉(zhuǎn)矩而言計(jì)算精度只有微小提高,且都與試驗(yàn)結(jié)果較為接近,表明在預(yù)報船體自航性能時,僅考慮推力而忽略轉(zhuǎn)矩影響也能得到具有相當(dāng)精度的計(jì)算結(jié)果,而此時整個計(jì)算域流場沿船體中縱剖面對稱.因此,在建模計(jì)算時可以僅建立半船模型,使得網(wǎng)格數(shù)量減少50%,以進(jìn)一步減少計(jì)算時間,從而實(shí)現(xiàn)對船體的自航性能進(jìn)行更為快速和準(zhǔn)確的預(yù)報.4計(jì)算結(jié)果與分析本節(jié)的主要內(nèi)容是應(yīng)用體積力法預(yù)報KCS集裝箱船模型在帶舵船模自航點(diǎn)條件下的船體自航性能.若假設(shè)船體為固定姿勢,則計(jì)算流程可完全參照第3節(jié)的計(jì)算步驟.但是,此次模擬計(jì)算參考的試驗(yàn)所采用的螺旋槳型號為SVP1193,與第3節(jié)中所用的螺旋槳型號不同.試驗(yàn)在船體自由升沉和縱傾條件下進(jìn)行,而使用CFD預(yù)測船體的姿態(tài)相當(dāng)耗時,為此,本節(jié)采用的計(jì)算策略是:采用試驗(yàn)結(jié)果的升沉值σ和縱傾值π固定船體姿勢,求取船模平衡點(diǎn)推力TB.同時在忽略轉(zhuǎn)矩影響的條件下,所得到船體自航性能相關(guān)參數(shù)仍具有相當(dāng)高的精度.因此,本節(jié)僅對比船體總阻力、摩擦阻力(由ITTC公式估計(jì))以及螺旋槳推力,來考察體積力法在求解船槳舵干擾問題時的計(jì)算精度.建立帶舵的整船模型時考慮到船體的升沉和縱傾,首先對船體航行姿態(tài)進(jìn)行調(diào)整,整體下沉8.33mm和艏傾0.143°,然后分區(qū)劃分網(wǎng)格,網(wǎng)格單元總數(shù)約為1.15×106.本文的模擬計(jì)算所采用的計(jì)算策略和邊界條件與算例1基本相同,但考慮到舵的影響,為提高計(jì)算精度,作了如下改動:動量方程中的對流項(xiàng)采用三階MUSCL(MonotoneUpstream-CenteredSchemesforConservationLaws)格式,體積分?jǐn)?shù)方程采用修正HRIC(HighResolutionInterfaceCapturing)格式,湍流方程采用k-ωSST(Shear-StressTransport)湍流模式封閉方程.計(jì)算策略仍是首先在非定常條件下采用PISO進(jìn)行計(jì)算,當(dāng)自由面的形狀相對穩(wěn)定后,將非定常條件改為定常條件,采用SIMPLE算法并適當(dāng)減小松弛因子進(jìn)行計(jì)算.計(jì)算時首先加入推力T′B=25N,得出總阻力R′T,SP=24.153N,進(jìn)而可以得出T′e=T′B-R′T,SP=0.847N;然后再加入推力T″B=23N,得出R″T,SP=23.868N,進(jìn)而再得出T″e=T″B-R″T,SP=-0.868N,并最終插值得出平衡點(diǎn)推力TB=24.02N,自航阻力RT,SP=24.134N.計(jì)算結(jié)果與試驗(yàn)值的對比如表6所示.表中,RF,SP為摩擦阻力.由表6可見,總阻力誤差較大,最大達(dá)到-7.01%.一般而言,對整船帶舵的自由面黏性流計(jì)算,網(wǎng)格總數(shù)1.15×106不夠充分;另外,由于SV