潜艇旋臂回转试验数值模拟

肖昌润，刘瑞杰，许可，刘洋，徐亚运

(1.海军工程大学舰船工程系，湖北武汉430033)(2.中国人民解放军91039部队，北京102401)

(3.中国人民解放军92339部队，广东湛江524000)

潜艇是海军重要的作战武器，对于维护国家海权、协同水面舰船作战、进行战略威慑具有重要意义.获取潜艇水动力系数是评价其操纵性的重要途径.以节省经费和缩短研究周期为目的，计算流体力学(computational fluid dynamics，CFD)预报潜艇操纵性越来越受到重视，且已经成为重要学科.随着应用和研究的深入，用计算流体力学方法获得的水动力系数越来越准确，计算精度越来越受到认同.

文献［1-5］中对船舶定常计算进行了研究，分析了网格和湍流模型对计算精度的影响，为模拟回转运动起到了指导作用.国内通过数值计算方法模拟回转运动的文献较少，其中文献［6］中通过MRF方法对回转运动进行了数值模拟，分析了回转运动的原理，并进行了初步研究.国外对于回转运动的研究已经比较成熟，文献［7］中通过变形方法预报流动分离情况，为研究回转运动提供了新的思路.文献［8-10］中通过自编程序对回转运动进行了研究，分析了湍流模式、网格结构、边界层等对计算精度的影响.文献［11］中通过稳流NS方程和旋转坐标系的方法计算了最大漂角达到18°的回转运动，计算结果误差在11%以内.

基于6自由度的潜艇操纵性方程涉及到固定于地球的固定坐标系和固联于潜艇的运动坐标系.通过CFD对潜艇非定常运动进行计算时需采用动网格.动网格参考坐标系可以是固定坐标系也可以是固联于潜艇的运动坐标系.文中的Mesh Motion方法采用地球坐标系，而基于源项的方法采用位于潜艇的运动坐标系.文中以潜艇SUBOFF模型为研究对象，通过流体计算软件FLUENT，对潜艇回转运动进行了研究，分析了两种方法对潜艇水动力计算精度的影响，结果表明计算精度满足工程预报的要求.

1 基本理论与控制方程

1.1 控制方程

基于雷诺时均方法，用张量形式表述的控制方程［12］为:

式中:ui为速度;t为时间;ρ ui′uj′为雷诺应力项;ρ为水的密度.

1.2 Mesh Motion方法模拟旋臂回转运动

当采用地球坐标系时，网格整体随潜艇做回转运动.Fluent软件的Mesh Motion功能可以帮助实现这个过程.回转运动是非定常的，通过这种方式进行迭代时，网格不会发生任何变形，避免了计算时网格发生拉伸、重构，使网格质量变差，甚至产生负网格，影响计算精度.

基于6自由度潜艇操纵性方程，关于无因次化的横向力Y′和偏航力矩N′，有方程组:

式中:r为角速度;Y′r，N′r，N′r|r|为旋转水动力导数.

1.3 基于源项法的潜艇旋臂回转试验

由于潜艇做回转运动，在运动坐标系下的速度可以分解为线速度V和角速度Ω.此时，对于任意流体单元，其在地球坐标系下的速度可以分解为U=Ur+Ue，Ur为相对速度，Ue为牵连速度.Ue可以表示为Ue=V+Ω×r′，r′为相对于运动坐标系的坐标向量.则任意流体单元在固定坐标系下的绝对加速度为a=ar+ae+ac.

又Ue是关于时间的函数，故2U≡0，则上式可以简化为:

将 -ρae-ρac定义为源项，用MS表示.

通过引入MS源项，即可将运动坐标系下非定常的旋臂回转运动等效为地球坐标系下的定常运动.

令Ur= ［u，v，w］T，V= ［u0，v0，w0］T，Ω = ［p，q，r］T，r′=［x，y，z］T.其中u，v，w分别是纵向速度、横向速度、垂向速度;p，q，r分别是横倾角速度、纵倾角速度、偏航角速度;u0，v0，w0表示初始速度.潜艇绕Z轴做回转运动，则u0=v0=w0=0，u=v=0=0，p=q=0=0.消去这些项可得源项的分量形式:

通过相似的步骤也可以得到绕X，Y轴做回转运动的源项公式.将式(8)编译，借助UDF的DEFINE-SOURCE，函数代入Fluent即可对回转运动进行模拟.

2 计算区域网格划分

文中两种方法使用的网格及边界条件的设置是相同的.计算区域是一个旋转中心距离内壁面1.5L(L=4.356 m是模型的长度)的半环形.为了降低近壁面网格Y+值，提高网格质量，同时减少网格数量，采用分块网格划分技术.使用RNG k-ε湍流模型计算时，网格总数为450万，Y+值在30左右，并选用标准壁面函数.Mesh Motion方法中，入口及四周壁面均设为速度入口，整个计算区域设定为旋转域，潜艇壁面相对旋转域静止.旋转域的旋转角速度 r分别为 0.08，0.1，0.15，0.2 rad/s.源项法中由于入口处速度是相对速度，且随半径变化，Uinlet=-Ω ×r′，速度，故需要用到DEFINE-PROFILE函数，同时将DEFINESOURCE函数通过UDF接口添加到流场中，实现回转运动模拟.计算区域网格如图1.

图1 计算网格示意Fig.1 Schematic diagram of grids

3 计算结果校验与分析

3.1 Mesh Motion方法计算结果分析

由于网格绕旋转中心旋转，故力矩中心坐标也是旋转的，通过已知点坐标，换算可得到力矩中心坐标.仿真得到力Y和力矩N后，通过软件Matlab拟合三次样条曲线，得到相关水动力系数.力和力矩及其旋转导数的无因次化表达式分别为:.式中:ρ为流体的密度;U为潜艇旋转的线速度;L为艇长;Y为回转试验中沿Z方向的力;N为绕Y轴的力矩.

图2 拟合曲线Fig.2 Fitted curves

由图2的拟合曲线可得旋转导数计算结果，并与试验结果进行对比，如表1所示，误差满足工程计算需要.图3显示不同时刻流场的速度场分布示意图.

表1力和力矩计算结果与试验结果对比Table 1 Force and torque calculation results compare with the test results

图3 基于Mesh Motion方法不同时刻速度场分布及艇体周围流场细节Fig.3 Velocity distribution at different time and flow field detail around the submarine based on Mesh Motion

3.2基于源项法的潜艇旋臂回转试验计算

Mesh Motion功能计算是一个非定常的过程，其对计算机的性能要求较高，耗费大量时间.基于源项法的方式可以有效拟补Mesh Motion方法的不足.基于源项法的潜艇旋臂回转模拟是通过Fluent软件的源项功能，加入UDF自编程序实现的.计算是定常的，极大减少了工作量.

对比RNG k-ε湍流模型和SST k-ω湍流模型的计算精度和计算量.使用SST k-ω湍流模型时，由于SST k-ω要求Y+＜1，故网格数大增，网格数目为930万.计算结果如表2所示.

表2 计算结果与试验结果比较Tab.2 Comparison between calculation results and test results

从计算时间上看，基于源项法的回转运动计算效率明显高于Mesh Motion方法.由于网格数目大增，SST k-ω湍流模型大大增加了计算量，但是并没有明显提高精度.

4 结论

文中通过Mesh Motion方法和基于源项法旋臂试验对全附体的SUBOFF潜艇回转运动进行了模拟.对比了RNG k-ε湍流模型和SST k-ω湍流模型的计算精度，并与试验结果进行对比.结果表明:

1)Mesh Motion方法和基于源项法旋臂试验的计算精度满足工程需求.但基于源项法悬臂试验计算量更小，计算精度更高.

2)SST k-ω湍流模型对于网格Y+值要求较高，造成网格数增加，与RNG k-ε湍流模型对比，计算耗时更长，但并没有明显提高计算精度.

References)

［1］ 刘志华，熊鹰，韩宝玉.潜艇流场数值计算网格与湍流模型选取［J］.华中科技大学:自然科学版，2009，39(7):98-101.Liu Zhihua，Xiong Ying，Han Baoyu.Computational grid and turbulent model for calculating submarine viscous flow field［J］.Journal of Huazhong University of Science and Technology:Natural Science Edition，2009，39(7):98-101.(in Chinese)

［2］ 肖昌润，刘巨斌，朱建华.DARPA2潜艇模型定常绕流水动力数值计算［J］.华中科技大学:自然科学版，2007，35(8):115-118.Xiao Changrun，Liu Jubin，Zhu Jianhua.Numerical computation of hydrodynamic force of DARPA2 submarine model［J］.Journal of Huazhong University of Science and Technology:Natural Science Edition，2007，35(8):115-118.(in Chinese)

［3］ Xing T，Bhushan S，Stern F.Vortical and turbulent structures for KVLCC2 at drift angle 0，12 and 30 degrees［J］.Ocean Engineering，2012，55:23-43.

［4］ 郑小龙，黄胜，尚秀敏.基于CFD的船舶阻力预报方法研究［J］.江苏科技大学学报:自然科学版，2014，28(2):109-113.Zheng Xiaolong，Huang Sheng，Shang Xiumin.Study of ship resistance prediction method based on CFD［J］.Journal of Jiangsu University of Science and Technology:Natural Science Edition，2014，28(2):109-113.(in Chinese)

［5］ 陈淑玲，杨松林，刘智.基于Fluent的五体船静水中水动力特性数值模拟［J］.江苏科技大学学报:自然科学版，2012，26(6):541-545.Chen Shuling，Yang Songlin，Liu Zhi.Numerical simulation of hydrodynamic performance of pentamaran in calm water based on Fluent［J］.Journal of Jiangsu University of Science and Technology:Natural Science Edition，2012，26(6):541-545.(in Chinese)

［6］ 刘帅，葛彤，赵敏.基于源项法的潜艇旋臂试验模拟［J］.大连海事大学学报，2011，37(2):1-4.Liu Shuai，Ge tong，Zhao Min.Simulation for submarine rotating-arm test based on added momentum source method［J］.Journal of Dalian Maritime University:Natural Science Edition，2011，37(2):1-4.(in Chinese)

［7］ Zhang J T，Maxell J A，Gerber A G，et al.Simulation of the flow over axisymmetric submarine hulls in steady turning［J］.Ocean Engineering，2013，57:180-196.

［8］ Toxopeus S，Atsavapranee P，Wolf E，et al.Collaborative CFD exercise for a submarine in a steady turn［C］∥International Conference on Ocean，Offshore and Arctic Engineering.Rio de Janeiro:ASME，2012:761-772.

［9］ Gregory P A，Joubert P N，Chong M S.Flow over a body of revolution in a steady turn［R］.Rockingham:Defence Science and Technology Organisation Victoria Platform Sience Lab，2004.

［10］ Phillips A B，Turnock S R，Furlong M.The use of computational fluid dynamic to aid cost-effective hydrodynamic design of autonomous underwater vehicles［J］.Journal of Engineering for the Maritime Environment，224(4):239-254.

［11］ Marshallsay P G，Eriksson A M.Use of computational fluid dynamics as a tool to assess the hydrodynamic performance of a submarine［C］∥Australasian Fluid Mechanics Conference.Launceston:［s.n.］2012.

［12］ SAS IP.Ansys fluent user′s guide，version 14.5.2013.