操纵运动船舶的水动力计算研究.doc

-专业文档，值得下载!-专业文档，值得珍藏！-操纵运动船舶的水动力计算研究李冬荔，杨亮，聂武摘要：以船舶操纵水动力预报为研究背景，通过对商用计算流体力学软件FLUENT的二次开发，采用其动网格技术以及后处理系统，对大型船舶操纵水动力导数进行了数值计算.船体按照斜航、不同舵角、纯横荡和纯首摇等状态做运动，得出随船坐标系下作用于其上的水动力及力矩.通过进行基于最小二乘法的曲线拟合，最终求得船舶操纵水动力导数.计算结果与势流理论计算结果一致，表明了所提出的计算方法适用于复杂船舶运动的水动力导数计算.关键词：船舶操纵；水动力导数；FLUENT；动网格中图分类号：U661.33文献标识码：A文章编号：1000-6982(2009)02-0008-05ComputationalinvestigationofhydrodynamicforcesaroundmaneuveringshipLIDong-li,YANGLiang,NIEWu(StateKeyLaboratoryofAutonomousUnderwaterVehicle,HarbinEngineeringUniversity,Harbin150001,China)Abstract:Inordertopredicttheshipmaneuverabilityinviscousflows,hydrodynamicderivativesrelatedtoshipmaneuveringwereanalyzedusingadynamicmeshmethodandpost-processingsystem,whichwasbasedonthecommercialsoftwareFLUENT.Shipmaneuveringmotionsofthedifferentheadingangles,rudderangles,swayingandyawingareconsidered.Hydrodynamicforcesandmomentsactingonamaneuveringshipareobtainedinthebody-fixedcoordinatesystem.Thushydrodynamicderivativesoftheshipmaneuveringmotionsarecalculatedbyleastsquarecurvefitting.Thecomputationalresultsareconfirmedbythedataofpotentialtheorymethod,whichshowthatthismethodcanbeusedtocalculatethehydrodynamicderivativesofshipcomplicatedmotions.Keywords:shipmaneuvering；hydrodynamicderivatives；FLUENT；dynamicmesh0引言操纵运动船体水动力的数值计算方法可分为两大类，一是势流方法，二是粘性流方法.其中，势流方法对水动力预报精度能基本满足工程要求，成为日趋成熟且实用的方法.而三维船舶粘性流的数值研究是从20世纪60年代后期才开始的.在1996年第21届ITTC会议上1，求解船舶粘性流的RANS方法被认为已经成熟.随着现代船舶科技水平的不断提高，对船舶水动力方面的研究也提出了更高要求.许多新技术，如网格生成技术、并行计算技术、多重网格加速等，被广泛应用于复杂的流场计算中，形成了国际上船舶操纵水动力预报的热门和前沿课题.本文通过FLEUNT的非定常动网格的方法，数值求解了船舶操纵水动力导数.此方法可应用于一系列和船舶操纵问题相关的粘性流动与水动力计算，得到了令人满意的结果.1船舶操纵运动水动力导数的数值计算模型在普通的长条形流场区域内，船体不前进，水以一定速度流向船体，船体就可以作各种规定的运动以及频率很低的振荡.测量船体所受的水动力，可求得计算船舶操纵运动所需的各种加速度导数和速度导数.船舶操纵运动是低频运动，要求接近于零频率的水动收稿日期：2008-12-16；修回日期：2009-01-14基金项目：国家自然科学基金资助项目（50879014）9力导数.船舶操纵水动力导数计算是根据船舶所做各种运动计算水动力，可以分为以下五种工况（各种运动情况如图1所示，各运动参数都在随船坐标系下获得，计算公式见参考文献2）.数值计算船体所受水动力需达到一个稳定周期，此时方可将水动力（随船坐标系下船体所受侧向力Y与绕船中的力矩N）与运动状态的函数关系进行非线性最小二乘法拟合，即可得到各个水动力导数.此处稳定的含义为：当船体周期运动到同一相位时，相邻两次受力大小相差甚微.图1船体所做的各种运动2FLUENT的求解原理及应用应用FLUENT软件研究船体操纵运动的非定常水动力性能是本文的主要工作，这项工作的核心方法是基于有限体积法对雷诺平均纳维斯托克斯方程（RANS方程）及其补充、修正方程的离散的分区网格法.将该方法直接应用到基于VOF模型的PISO算法是数值计算的最主要内容3.2.1三维非结构网格中动网格技术的实现数值模拟船舶操纵运动中，所采用的动网格技术为弹簧近似法4.该方法是将非结构网格的整个区域看作一张具有弹性的网，每条边都看作是一根独立的具有倔强系数的弹簧，于是节点i，j间的弹簧张力为：ijijjiFKxx（1）式中，Kij为连结节点的位置矢量.认为移动后网格点的位置满足各点在初始状态的受力情况，即：ijijFS（2）式中，Si为初始状态时点i所受的合力，对所有与节点i相连的节点进行求和.图2为船体表面网格以及对整个流场取剖面后的网格结构示意图，船体表面由型线图画出，其表面网格按外形布置疏密可有效保证运动物体边界附近区域的网格质量.此外，为了有效减少三维动网格计算量，将整个流场沿船体中纵剖面（ZX平面）与标准吃水处的自由表面（YX平面）划分为四个子区域，并且对接近船体的中纵剖面处和自由表面处进行了单独的网格划分与加密，从而有效控制整个区间网格的疏密程度，并能达到计算要求.图2网格结构示意图2.2边界条件的设定计算区域内自由液面以下的入口、出口处采用给定速度分布；自由液面以上均为标准大气压，计算区域的侧面以及底面均为不可穿透的固壁面，船体的运动规律由UDF程序给出.该处所采用的计算模式与控制参数主要是针对基于PISO算法的流场数值计算方法.考虑到非定常运动，采用一阶隐式格式方法离散时间项，此外，考虑到自由液面的VOF处理方法，压力项的离散格式采用将体积力加权的处理方法，体积分数的离散格式采用几何重构的处理方法，湍流模式采用RNGk-模式，并对壁面处使用标准壁面函数.在RNGk-模型中，通过在大尺度运动和修正后的粘度项体现小尺度的影响，使这些小尺度运动有系统地从控制方程中去除.3数值计算结果及分析3.1数值方法验证为验证文中计算方法的有效性，应用前面介绍的求解RANS方程的数值方法对6:1长椭球体在垂直平面内做定常斜航运动时的流场和水动力进行计算，计算网格结构图如图3所示.并将计算结果与Wetzel，和图36:1长椭球体网格结构图不同舵角斜航运动纯横荡运动纯艏摇运动组合运动舵角=0漂角x0y0yyyx0x0x0y0y0y0yyyyyyyyyyyyyxxyxxxxxxxyxxx漂角固定xyzxyz10Simpson6的试验数据及Rhee和Hino7,8的计算结果进行比较，图4为长椭球体升力系数Cl与力矩系数Cm计算结果.结果表明，本文的数值计算方法能够反映以湍流分离流为主的复杂粘性流动的主要物理特征.图4长椭球体升力系数与力矩系数计算结果比较3.2水动力导数的计算结果数值模拟采用实船大小，设计水线间长286m，型宽36.7m，型深27.2m，设计吃水9.4m，航速18kn.1）变舵角运动计算控制导数，每个舵角下采用的时间步长均为0.01s，数值计算6480时间步后船体所受侧向力与力矩结果稳定.图5为不同舵角下（弧度）船体所受无量纲侧向力Y与力矩N曲线，得到无量纲化的控制导数.拟合采用如下matlab命令流：；(,2)；(1)*.3(2)*2(3)*(4),)；(,0,0,)1(1)*3(2)*2(3)*(4)；(,1)typpolyfittyfuninlinextxtxtxxtxlsqcurvefitfuntyyxtxtxtxplottyty（3）图5不同舵角下侧向力与力矩曲线2）斜航运动计算位置导数的方法同上，图6为不同漂角下船体所受无量纲侧向力Y与力矩N曲线（不同的漂角得到不同的侧向速度），线性导数结果有表一给出，非线性导数vvvY=-0.30948，vvvN=0.00131.图6不同侧向速度下侧向力与力矩曲线3）纯横荡运动时，振幅a=3.3m，横荡频率取0.02、0.04、0.05、0.08、0.1，来流速度V=9.26m/s（18kn），将一个数值周期取1600个时间步进行计算，计算4个周期，发现在第3个周期船体所受平均阻力、侧向力与力矩结果均收敛.图7为f=0.04Hz计算收敛后一个周期内的无量纲化的受力曲线.图7纯横荡运动时船体受侧向力与力矩曲线纯横荡运动时船体在随体坐标系下的水动力表达式的无量纲化形式为：222321.111222VVYYYvvVVVLLVLL（4）将运动表达式sinyat、cosvyat、2sinvyat代入式（4），可得22()sin()cossincosVVaLaYYtYtVVAtBt（5）其中，22()VaLAYV，()VaBYV，对各横荡频率下的受力曲线进行非线性拟合得到各系数项A、B.由于步骤2）斜航试验给出的位置导数较准确，所以只用该方法确定线加速度导数Yv，即对公式22()VaLAYV进行第二次拟合，直线斜率即为VY,同理求得VN值，图8给出vY拟合结果.0.0120.0100.0080.0060.0040.0020.000-0.002-0.004-0.006-0.008-0.0100102030攻角/(°)0.060.050.040.030.020.010.00-0.01-0.02-0.03-0.04-0.05力矩系数Cm升力系数ClRhee等力矩系数计算值,CmWetzel等力矩系数计算值,Cm本文力矩系数计算值,CmRhee等升力系数计算值,ClWetzel等升力系数计算值,Cl本文升力系数计算值,Cl0.0150.0100.0050.000-0.005-0.010-0.015-0.6-0.40.00.4舵角/(°)侧向力Y力矩N船体侧向力数值解,Y船体力矩数值解,Ny=-0.0000095+0.0023x+0.0000047x2-0.0019x3,Y0.0150.0100.0050.000-0.005-0.010-0.015-0.20.20.6y=0.0000036-0.00105x+0.00002x2+0.0075x3,N0.030.020.010.00-0.01-0.02-0.5-0.40.00.4侧向速度V侧向力Y力矩N船体侧向力数值解,Y船体力矩数值解,Ny=-0.00009-0.03752x+0.00007x2-0.05158x3,Y-0.20.20.5y=0.000002-0.0052x+0.00025x2+0.00022x3,N0.0030.0020.0010.000-0.001-0.002-0.3-0.10.10.30.0100.0050.000-0.0050时间t/s侧向力Y力矩N船体侧向力数值解,Y船体力矩数值解,Ny=0.03192sint-0.002475cost,Y20y=0.00186sint-0.000525cost,N100.0100.0050.000-0.005YN11图8纯横荡运动水动力导数vY拟合结果4）纯首摇运动时，数值求解过程同步骤3）.图9为频率为0.04Hz时计算稳定后一个周期内受力曲线，图10为加速度导数拟合结果，数值表一给出，非线性导数rrrY=-4.52E-4，rrrN=-9.06E-5.图9纯首摇运动时船体受力曲线图10纯首摇水动力导数rY、rN拟合结果为验证本文计算结果的准确性，另对上述中的线性水动力导数采用势流理论研究，其中，计算惯性力采用无升力势流理论及两因次切片法，计算粘性力采用细长体理论、短翼理论和分离绕流理论，并且将裸船体和附体所受到的水动力分开来计算，再进行叠加，具体计算方法参见文献9，计算结果如表1所示，二者具有明显的一致性.表1本文数值方法和势流方法计算值比较水动力导数(×104)本文数值计算值势流理论计算值Y23266N-10.5-12.4vY-375-380vN-52-55vY-26-22vN0.971.2rY-14.3-14.6rN-15.4-12.7rY-3.2-2.8rN-2.44-2.295）漂角加首摇组合试验确定交叉耦合导数，由于常数和2220sinrrt在一个周期内区间变号积分之值为零.以0sinrrt，20arV代入水动力表达式无量纲化，并在一个周期内进行变号积分得：20236322620236322611(,)()(,)()44111(sin)()()29611(,)()(,)()44111(sin)()()296rvvrrrrrvvrrrrYvrdtYvrdtaLaLYYYVVNvrdtNvrdtaLaLNNNVV（6）取频率为0.02Hz，分别计算了漂角为0.2°、0.4°、0.6°、0.8°、1.2°、1.4°、1.6°、1.8°几种情况的组合运动下的水动力，采用离散化函数积分形式计算等式左边的积分值，即一个周期1600步，步长0.03125s，得到的一个周期内每个步长对应的水动力值计算区间变号积分，带入已求得的rY,rrrY，可以得到vvrY的值.同理求出力矩的交叉耦合导数vvrN.扣除斜航运动求得的Yv、Yvvv，即可求得Yvrr.表2交叉耦合导数求解结果水动力导数vvrYvvrNvrYvrN本文数值计算值0.02850.010080.02770.030844结论基于CFD软件FLUENT二次开发，采用RNGk-模型，数值计算了实船尺度下船体作操纵运动的水动力导数，计算结果表明：1）非结构网格的变形移动，该方法可以保证在物面附近的网格质量，同时适用船体斜航、不同舵角、纯横移和纯摇首等状态大幅度运动.0.0150.0050.000-5-aL2/V2侧向力系数项A(Y)数值解-2y=0.00221-0.0026x-30.010-4-100.0030.0020.000-0.00230时间t/s侧向力Y力矩N船体侧向力数值解,Y船体力矩数值解,Ny=0.0018sint+0.00144cost-0.000064sin3t,Y50y=0.00068sint+0.0014cost-0.0006sin3t,N400.0060.0020.000-0.004YN0.001-0.001-0.0020.0040.0150.005-30-a3L2/V3侧向力系数项B(Y)力矩系数B项(N)侧向力系数,B(Y)力矩系数,B(N)y=0.00032x,B(Y)-50y=-0.000244x,B(N)-400.0100.0000.0150.0050.0100.000-20-100