基于cfd方法的moeri集装箱船自航性能数值模拟

基于cfd方法的moeri集装箱船自航性能数值模拟

目前，研究船桨干扰的数值方法可分为两类：一般建模和体积力法。在使用船桨通用建模方法研究船桨干扰问题时，应全面报告信息，以反映真流场中螺旋作用区域的复杂性，因此使用该方法是广泛的，但计算往往需要很长时间。体积力法较整体建模方法而言,建模简单,计算效率高,能预报流场的基本信息,因此得到一定的应用.体积力源项通常可由基于势流理论的螺旋桨模型计算得到,依据它与雷诺平均纳维-斯托克斯(RANS)方程的耦合方式,体积力法可以分为迭代型和描述型两种.迭代型体积力法有较多应用,但该方法计算过程较为繁锁,需要在螺旋桨性能计算和RANS方程计算之间进行相互迭代.描述型体积力法则相对简便,对实际螺旋桨的几何参数要求十分简单,且不需要螺旋桨性能计算,但该方法无法直接得到螺旋桨性能参数而须采用其他方法.本文通过计算流体力学(CFD)软件Fluent的RANS求解器,采用文献中的描述型体积力模型,分别在无舵实船自航点和有舵模型自航点2种条件下对MOERI(KoreaResearchInstituteforShipsandOceanEngineering)集装箱船(KCS)的自航性能进行了数值预报,并通过等推力法原理、敞水计算和相对应的螺旋桨敞水曲线获得船体自航性能相关参数,通过与试验结果的对比验证了该方法的有效性.1数学模型1.1均速度分量a流场的控制方程为连续方程和动量方程:∂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ˆ为无量纲系数;¯u′iu′ju′iu′j¯¯¯¯¯¯¯¯为雷诺应力;f*i为无量纲体积力,代表螺旋桨力源项.1.2桨叶的推力及转速本文采用文献中的方法施加体积力源项,该方法通过对螺旋桨所在区域施加推力和转矩代替桨叶表面载荷,公式如下: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,θ分别为推力和转矩在螺旋桨区域内的分布;r为螺旋桨区域内任意一点至螺旋桨轴线的距离;YH=RH/RP为毂径比,RP、RH分别为螺旋桨半径和桨毂半径;Δx为螺旋桨在轴向的投影弦长,一般取桨毂厚度δ;Ct=4T/(ρu2πDP)为一种无量纲推力系数,T为螺旋桨推力,ρ为流体密度,u为船速,DP为螺旋桨直径;KM=M/(ρn2D5P)为转矩系数,M为转矩,n为转速;J=va/(nDP)为螺旋桨进速系数,va为螺旋桨进速.在实际计算中,施加体积力的具体过程是将推力和转矩施加给螺旋桨作用区域的网格单元进而对控制方程进行求解计算.由于船艉形状复杂,网格单元与螺旋桨区域边界不易重合导致实际体积力大小与指定值相比会有一定误差.为此,本文使用UDF(UserDefinedFunction)施加体积力时,在UDF源代码中添加计算推力的DEFINE_ADJUST宏和修正项,用以查看实际推力大小并进行修正以达到指定值.需要说明的是,修正后推力T和转矩M的变化率相同,这表明通过修正项对推力进行修正时,也能对转矩进行精确的修正,进而达到准确地施加体积力的目的.2计算条件2.1设计水线和船模的点坐标本文的研究对象为KCS集装箱船模型.该船型关于实船自航点和模型自航点的试验数据取自文献.计算中参考坐标系的坐标原点位于船艏中心线与设计水线的交点,x轴为中纵剖面和设计水线面的交线,正方向指向船艉;y轴正方向指向船体右弦;z轴正方向向上.本文共研究了2个算例.算例1为静水中实船自航点的自航试验,船体附体只包括桨;算例2则是静水中船模自航点的自航试验,其中船模可自由升沉和纵倾,船体附体包括桨和舵.以上2个算例的具体模拟计算条件如表1所示.表中:LPP为船体垂线间长;Fo为傅里叶数;Re为雷诺数;θ为温度;Sw为湿表面积.2.2网格划分参数本文的主要内容是预报船体自航性能.为提高计算精度和收敛性,采用分块结构化网格,使用CFD网格生成软件Gambit2.4生成网格,其中2个算例的计算域均为长方形.具体的网格划分参数如表2所示.表中:y+为壁单位;Δyp为最小网格间距.网格划分的难点在于船艉区域,2个算例的船艉区域网格划分如图1所示.3近壁流动动力学方程本节采用商用RANS求解器FLUENT6.3中的VOF方法在船模尺度下对自由液面情况的三维黏性流场进行数值求解,即用有限体积法离散RANS方程,并采用压力基求解器进行求解计算.计算中PISO(PressureImplicitwithSplittingofOperators)算法和SIMPLE(Semi-ImplicitforPressure-LinkedEquations)算法先后被用于速度和压力间的耦合计算.湍流方程采用带旋流修正k-ε湍流模式封闭方程,并使用非平衡壁面函数简化处理近壁流动.数值离散方面,动量方程中的对流项采用三阶QUICK(QuadraticUpwindInterpolationofConvectiveKinematics)格式,k-ε方程和体积分数方程均采用二阶上风格式.计算域的边界条件设置为:进流面,表面压力为0;上表面,滑移边界条件;船体表面,无滑移边界条件;侧表面,滑移边界条件;下表面,滑移边界条件;出流面,表面压力为0.在非定常计算过程中,取时间步长Δt=0.5ms,最大内循环步数为10步.3.1试验数据分析在应用体积力法预报船体自航性能之前需先计算出裸船体的总阻力系数CT、摩擦阻力系数CF、兴波阻力系数CP及标称伴流系数φ.如表3所示,阻力系数和标称伴流系数的计算值与试验数据吻合良好,相对偏差均在5%以内.3.2预报忽略转速的船体自航性能为获得船体自航性能参数,本文的计算流程分为2步:(1)预报忽略转矩的船体自航性能,体积力中仅包含推力T;(2)预报考虑转矩影响的船体自航性能,体积力中同时包含推力T和转矩M.3.2.1旋转进速系数(1)依据等推力法原理求取平衡点推力TB.由于模拟计算的是实船自航点的自航试验,考虑到摩擦阻力修正,平衡点推力的求取依据下式:Τe=Τ-RΤ,SΡ+RF(8)Te=T−RT,SP+RF(8)式中:T为所加入的螺旋桨推力;RT,SP为加入推力T后计算得到船体自航总阻力;RF为实船和船模间的摩擦阻力修正值,本文取RF=30.25N,与该自航试验中的取值相同;Te=0所对应的T即为平衡点推力TB.实际求取Te时是根据不同的推力T′、T″得出不同的T′e、T″e,这两者不一定等于零,但要保证异号并且其绝对值应尽可能小.然后通过T′、T″线性插值得出平衡点推力TB.本文中的计算,先后尝试了T=59,57,56,55N才最终确定当T′=55N,T′e≈-0.646N以及T″=56N,T″e≈0.438N时,两者异号,再通过线性插值得出平衡点推力TB=55.6N,如图2所示.(2)求解螺旋桨进速va.其原理是先将仅包含平衡点推力TB=55.6N的体积力加入到3.1节得到的计算结果中再次计算,待计算收敛导出船后桨盘面处平均速度的平方,即ˉv2diskv¯2disk.本文得出的ˉv2disk=3.6636m2/s2v¯2disk=3.6636m2/s2,相应的RT,SP=86.376N,故推力减额系数Δk=(RT,SP-RT,Tow)/TB=0.870.其中,RT,Tow为3.1节中裸船体总阻力.为了得到va,运用动量守恒原理进行敞水计算,即无自由液面的情况下,赋予螺旋桨体积力为平衡点推力TB,改变来流速度v,使得桨盘面处ˉv2diskv¯2disk与上文中ˉv2diskv¯2disk相等,则此时的来流速度v即为螺旋桨进速va.通过敞水计算得出ˉv2disk=3.6634m2/s2v¯2disk=3.6634m2/s2对应的进速va=1.7617m/s.(3)求解螺旋桨进速系数.通过下式得到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绘制在同一图中并取两者交点,即可求得进速系数JB.具体做法如图3所示,其中进速系数JB=0.749.(4)通过下式求解忽略转矩的船体自航性能:推力系数ΚΤ=Τ/(ρn2D4Ρ)(10)KT=T/(ρn2D4P)(10)转矩系数ΚΜ=Μ/(ρn2D5Ρ)(11)KM=M/(ρn2D5P)(11)螺旋桨转数n=va/(JDΡ)(12)n=va/(JDP)(12)敞水效率ηΟ=ΚΤJ/(2πΚΜ)(13)船身效率ηΗ=Δk/(1-φ)(14)推进效率η=ηΟηRηΗ(15)式中,ηR为相对旋转效率,取ηR=1.通过上述公式计算得到的忽略转矩影响的船体自航性能相关参数如表4所示.表中,e为误差.3.2.2基于扭矩的自航性能计算若要考虑转矩的影响,则需要将同时包含表5中的转矩M和平衡点推力TB=55.6N的体积力再一次加入到3.1节所得到的裸船体计算结果中进行计算,然后按照上节的步骤依次求取螺旋桨进速va和进速系数J,再通过式(10)～(15)进行计算便可得到考虑转矩影响的船体自航性能相关参数.在考虑转矩情况下,经计算可分别得到进速va=1.768m/s和进速系数J=0.748(见图4),其他自航性能参数如表5所示.图5、6为桨盘面后方0.25DP处轴向速度分布的计算值和试验值.由图5,6可见,考虑转矩影响时,桨盘面轴向速度分布体现了螺旋桨的旋转效应,但与试验结果仍存在一定差异.对比表4、5可知,考虑转矩较忽略转矩而言计算精度只有微小提高,且都与试验结果较为接近,表明在预报船体自航性能时,仅考虑推力而忽略转矩影响也能得到具有相当精度的计算结果,而此时整个计算域流场沿船体中纵剖面对称.因此,在建模计算时可以仅建立半船模型,使得网格数量减少50%,以进一步减少计算时间,从而实现对船体的自航性能进行更为快速和准确的预报.4计算结果与分析本节的主要内容是应用体积力法预报KCS集装箱船模型在带舵船模自航点条件下的船体自航性能.若假设船体为固定姿势,则计算流程可完全参照第3节的计算步骤.但是,此次模拟计算参考的试验所采用的螺旋桨型号为SVP1193,与第3节中所用的螺旋桨型号不同.试验在船体自由升沉和纵倾条件下进行,而使用CFD预测船体的姿态相当耗时,为此,本节采用的计算策略是:采用试验结果的升沉值σ和纵倾值π固定船体姿势,求取船模平衡点推力TB.同时在忽略转矩影响的条件下,所得到船体自航性能相关参数仍具有相当高的精度.因此,本节仅对比船体总阻力、摩擦阻力(由ITTC公式估计)以及螺旋桨推力,来考察体积力法在求解船桨舵干扰问题时的计算精度.建立带舵的整船模型时考虑到船体的升沉和纵倾,首先对船体航行姿态进行调整,整体下沉8.33mm和艏倾0.143°,然后分区划分网格,网格单元总数约为1.15×106.本文的模拟计算所采用的计算策略和边界条件与算例1基本相同,但考虑到舵的影响,为提高计算精度,作了如下改动:动量方程中的对流项采用三阶MUSCL(MonotoneUpstream-CenteredSchemesforConservationLaws)格式,体积分数方程采用修正HRIC(HighResolutionInterfaceCapturing)格式,湍流方程采用k-ωSST(Shear-StressTransport)湍流模式封闭方程.计算策略仍是首先在非定常条件下采用PISO进行计算,当自由面的形状相对稳定后,将非定常条件改为定常条件,采用SIMPLE算法并适当减小松弛因子进行计算.计算时首先加入推力T′B=25N,得出总阻力R′T,SP=24.153N,进而可以得出T′e=T′B-R′T,SP=0.847N;然后再加入推力T″B=23N,得出R″T,SP=23.868N,进而再得出T″e=T″B-R″T,SP=-0.868N,并最终插值得出平衡点推力TB=24.02N,自航阻力RT,SP=24.134N.计算结果与试验值的对比如表6所示.表中,RF,SP为摩擦阻力.由表6可见,总阻力误差较大,最大达到-7.01%.一般而言,对整船带舵的自由面黏性流计算,网格总数1.15×106不够充分;另外,由于SV