水下航行体回转水动力数值计算研究

卢锦国 梁中刚 吴方良 周轶美

中国舰船研究设计中心，湖北武汉 430064

水下航行体回转水动力数值计算研究

卢锦国 梁中刚 吴方良 周轶美

中国舰船研究设计中心，湖北武汉 430064

准确预报水动力对水下航行体的安全操纵设计而言至关重要。基于商用软件FLUENT数值求解RANS方程，采用定常旋转坐标系，运用相对运动理论及运动叠加原理，并对控制方程中向心力源项空间离散误差进行修正，预报水下航行体做单平面回转运动的受力及力矩。在Re＝11.7×106条件下，力和力矩预报精度的偏差在12%以内，表明所采用的数值预报方法有效、可行，具有较好的工程实用价值。

回转运动；水动力；向心力修正；FLUENT；水下航行体

1 引言

操纵性能预报是水下航行体设计的重要环节。数值方法应用于操纵性能预报由来已久，针对不同的模型及不同的航行姿态，国内外众多研究人员开展了这方面的研究。公开发表的文献大多针对直航或偏航状态下水动力性能进行研究，对比试验数据来源于直线拖曳水池的试验结果。水下航行体旋臂水池试验通过在一定范围内变换半径、攻角及舵角来确定旋转导数、非线性及耦合水动力导数及力矩导数。试验研究表明，在较大半径旋臂水池试验条件下得到旋转导数，通常情况下比相应的PMM试验得到的旋转导数具有更高的精度［1］。然而，在对旋臂水池试验的数值模拟研究方面，由于流动复杂及试验数据来源少等原因，已经开展的工作较为有限，针对这方面的研究主要由美国泰勒水池研究人员发起。Sung等［2］融合了多块网格、多重网格、局部加密等技术求解三维不可压缩RANS方程，采用标准k－ω及修正的B－L湍流模型，数值模拟了回转主体在不同偏航角条件下的水下定常回转运动，作者对影响数值计算结果的条件进行了系统研究，求得的垂向力及纵倾力矩与试验值的误差在20%以内。2002年，Sung 等［3］人采用 standard k－ ω 及 realizable k － ε湍流模型，对ONR Body－1潜艇模型进行回转运动模拟，并与试验数据进行验证，对比较结果进行了合理分析。随着旋臂水池试验报告的逐步公开［4－6］，模拟水下航行体回转运动的研究将渐渐成为新的研究热点。

本文基于对旋臂水池试验的数值模拟，采用商用软件FLUENT，首先获得稳定的纯旋转运动流场，提出控制方程中向心力源项空间离散误差修正方法，进而预报某型水下航行体做单平面回转运动水动力性能，通过与已有试验数据进行比较，对影响数值计算结果的条件进行优选，最终建立起一套可行的回转水动力数值预报方法。

2 数值计算方法

2.1 物理描述

对于固定于大地的惯性坐标系而言，由于水下航行体的定常回转运动存在向心加速度，因此该运动在此坐标系下为非定常问题，将其转化为旋转运动坐标系下进行求解，可以使问题得到简化，如图1所示。

图1 水下航行体回转运动示意图Fig.1 Schematic of the submerge d vehicle in turningmotion

当运动方程在旋转系中求解时，流体的加速度（向心加速度及科氏加速度）将作为源项出现在动量方程中。

2.2 控制方程

在角速度一定的旋转运动坐标系下的连续性方程及不可压缩RANS方程为：

式中，xj与uj分别为笛卡尔坐标系下的坐标及速度分量；δij为克罗内克符号记法；ν为运动学粘性系数；τij为雷诺应力张量。

2.3 湍流模型

为了实现对回转运动水流横向分离复杂流动的准确捕捉，结合文献［7］给出的水下航行体在偏航状态下力和力矩的数值预报研究成果，可知SST k－ω湍流模型对带攻角来流条件下的垂向力和纵倾力矩预报具有较好的预报精度。本文采用SST k－ω湍流模型进行数值模拟研究。

2.4 模型及计算域

2.4.1 模型描述

旋臂试验模型为流线型回转体的58系列模型，编号4621，采用6阶多项式数学表达式描述。体长为4.572m、最大直径为0.623m、排水体积0.836 2 m3，形心位置位于回转轴线上，距首部2.037 3 m。进行模型试验时，采用侧装模型的方法，绕OY轴转动，取角速度Ωy＜0，重心线速度大小为2.572 m／s，以回转体长为特征长度，相应的雷诺数大小为Re＝11.7×106。

图2 计算模型示意图Fig.2 Schematic of the comput ation model

2.4.2 计算域选取

计算域的选取需要根据水下航行体的主尺度和回转半径来确定［2－3］。选取来流截面距航行体前端约1倍体长，外周边界距航行体中轴线约1倍体长，出流截面处距航行体尾端约3倍体长，如图3所示。

2.5 边界条件

在FLUENT软件中，复杂的数值模拟边界条件可以通过用户自定义函数 （UDF）功能编程实现。

图3 计算域示意图Fig.3 Schematic of the comput ation domain

水下航行体模型表面设置为无滑移壁面边界条件（u＝v＝w＝0），在旋转运动坐标系下，将此条件代入式（1），并忽略粘性影响，可以得到模型表面外法线方向应满足压力梯度为0条件，即

在流场中未放置模型的条件下，由描述流体运动的欧拉方程可知，外周边界的速度分布在某一固定时刻，同一半径曲线上任一点的速度方向与曲线在该点的切线方向重合，即外周边界为由一系列沿同一半径的流线组成的流管。因此，在得到稳定数值解的条件下，可以得到稳定的纯旋转运动流场。

2.6 网格划分及数值求解

2.6.1 网格划分

本文使用ICEM－CFD软件对流域进行多块结构化网格划分。对航行体边界层网格采用O型网格划分，主要绕流区域采用C型网格划分，外部流域采用H型网格划分，如图4所示。

图4 网格划分示意图Fig.4 Schematic of the grid partition

边界层网格考虑各湍流模型及壁面函数对粘性层的不同要求，参考FLUENT手册对y＋进行合理取值，y＋＜5，50%以上分布于［0，1］之间；层间增长率均设置为1.15。沿纵向及周向的网格布置考虑背流面、首部及尾段区域流动的复杂性，在这些区域划分较细密的网格节点；同时也在这些区域采用局部网格加密技术。

2.6.2 数值求解

流体运动的控制方程及湍流运动方程采用有限体积法进行离散，压力项采用标准差分格式，速度、动量及湍流参量采用2阶迎风差分格式。应用Simplec法处理压力与速度耦合问题，离散方程运用Gauss－Seidel迭代方法求解，并以MutilGrid技术加速计算收敛。

3 计算结果与对比分析

3.1 稳定的旋转流场

首先验证上述边界条件下求解三维不可压缩RANS方程可以获得稳定的旋转流场。本文采用环半径为10m、圆半径为5m的1／4圆环计算模型。比较了高度方向Z＝0 m平面内，沿半径方向R＝8、9、10、11、12m 的速度分布， 结果如表 1 所示。Z＝0截面各半径处速度分布的不均匀度（（最大值－最小值）／平均值）小于1%，数值耗散、截断误差等可能是造成这一差异的原因。图5所示为Z＝0截面速度随半径的分布曲线。

表1 Z＝0截面各半径处速度分布的不均匀度分析Tab.1 Nonuniform analysis of velocity distribution on the plane Z＝0 at different radii

图5 Z＝0截面速度随半径变化分布曲线Fig.5 C urves of velocity distribution on the plane Z＝0 at different radi i

由此可见，在本文前述边界条件下获得稳定旋转流场是完全可行的，符合流体运动规律。

3.2 向心力源项修正方法

对式（4）分析，当矢量点乘积为0时，两矢量或者相互垂直，或者至少其一为零矢量，由此可得：

对该式两边做体积积分，并对左端项应用奥高公式转化为面积积分，可以得到壁面压强分布与向心力源项的关系。但采用有限体积法建立离散方程时，很重要的步骤是将控制体积界面上的物理量及其导数通过节点物理量插值求出。虽然本文采用的离散格式在很大程度上减小了伪扩散现象，但针对本文的研究对象，仍至少存在以下2个问题：动量方程的求解基于直角坐标系下的各速度分量，其方向与网格线呈倾斜交叉；建立离散格式时没有考虑复杂形式源项的影响，把源项局部线性化的常规处理方法并不适用。为提高预报精度，需要对控制方程中的向心力源项做补偿修正。

在得到稳定旋转流场的基础上，本文选取与58－4621模型具有相同几何形状的流体介质体作为计算对象。理论上，在稳定的均质旋转流场中，这部分流体介质体所受的向心力Fc＝ρVrω2由介质所处的静压场提供。同时，对这部分流体介质体形心处的力矩应当为0。依据力的相互作用原理，由介质体表面压强积分得到的合力和合力矩应与航行体所受向心力及形心力矩一致。

应用上述理论方法计算得到的向心力及形心力矩未考虑离散格式所引入的误差、耗散。为此，首先对整个流域进行数值求解，在数值收敛的条件下，选取相同几何形状的包络边界作静压场积分，可以得到考虑误差、耗散的向心力及形心力矩。通过与理论值比较，可以对数值解算引起的误差做定量估算。运用这一原理，对R＝30.48 m条件下的多个湍流模型做计算，结果如表2所示。

表2 R＝30.48m条件下不同湍流模型向心力与形心力矩预报值Tab.2 Predictions of centripetal force and moment in different turbulencemodels at 30.48m turning radi us

可以看出，对于不同的湍流模型，计算结果与理论值的偏差各不相同，但相对偏差量很小，小于1%。该方法对边界层网格非常细密的情况，很难实现原有边界面两侧网格的均匀过渡，计算收敛困难，因此需要粗化原有的边界层网格，从而导致修正值在一定程度上失真；但是，该方法从另一个角度验证了本文选用的计算方法、边界条件及网格形式等对获得稳定旋转流场是准确、可行的。

另一种求取静压场的方法是假定水下航行体与流体以同一角速度做回转运动，流体与壁面间不存在相对运动，只受到流体由于回转运动而产生的静压场的作用。对 R ＝30.48 m、22.86 m、18.288 m、15.24 m 条件下 SST k － ω 模型进行计算，结果如表3所示。由表3可见，计算值与理论值的偏差均在5%以内。由于该方法不改变原有的网格分布，能较真实地反映出离散格式引入的误差，更适合作为本文向心力源项修正的方法。

表3 不同回转半径条件下向心力与形心力矩预报值Tab.3 Predictions of centripetal force and moment at different turning radi i

3.3 不同半径下的结果对比

在前述计算条件优选的基础上，数值模拟不同半径下航行体回转运动流场，对所受垂向力与纵倾力矩进行预报，并与试验值比较，如图6、图7所示。由图可见，数值预报的垂向力和纵倾力矩与试验值基本吻合，随着回转半径的减小、角速度值的增加，垂向力和纵倾力矩的预报值与试验值的偏差逐渐增大，且在数值上与试验值相比始终偏大，最大偏差不超过12%，但变化趋势与试验值相吻合。

图6 预报与试验垂向力系数随无因次回转角速度变化分布Fig.6 Effect of nondimensional turning radius，r＇＝（L ／R turn） on the computed and measured coefficients of normal force

图7 预报与试验纵倾力矩系数随无因次回转角速度变化分布Fig.7 Effect of nondimensional turning radius，r＇＝（L ／R turn） on the computed and measured coefficients of pitchingmoment

为了进一步分析预报结果，同时与分段模型试验结果进行比较，将航行体不均等地划分成10分段，划分方式与模型试验相同，见图2。模拟预报回转半径为R＝22.86m条件下，各个分段所受垂向力值、通过比较无因次垂向力沿纵向位置分布如图8所示。由图可见，首尾分段预报值与试验值偏差较大，导致总的垂向力与力矩存在偏差，尤其对纵倾力矩有较大影响。首部分段的偏差可能由于物理试验时模型首部装有激流装置，而数值计算则通过人工设定湍流参数（经验公式估算，尚无定论）来模拟，两者可能存在较大差异。进行几何建模时，考虑结构化网格划分的需要，对尾端采用小圆弧光顺，与模型的尖细尾端有较大差异，直接导致尾部分段值的偏差。同时，由于两方程模型的固有缺陷，SST k－ω模型对尾部3.5 m附近区段背水面复杂分离流动的捕捉能力存在一定的不足，可能导致预报值偏离试验值，湍流模型参数的调整可能对结果有所改善。

图8 R＝22.86m条件下预报与试验垂向力系数沿纵向位置分布曲线Fig.8 Distribution curves of the computed and measured coefficients of normal force along the longitudinal position of the vehicle at 22.86m turning radi us

此外，依据流体运动规律，随着回转半径变小，航行体首端与来流夹角逐渐增大，尾部去流段流动分离点前移，且由于角速度增大，航行体周围横向分离流逐渐变强，分离产生的涡旋运动愈发强烈，并且可能随流动下泄，出现严重的非定常分离流动，力和力矩则会出现不稳定波动。与文献［8］中对大攻角直航运动的研究相类似，此时需要考虑对流动进行非定常计算，这也是今后需要进一步研究的内容。

4 结 论

本文采用数值求解RANS方程及湍流模型模拟水下航行体旋臂水池试验，在选择适当计算条件的基础上，可以实现Re＝11.7×106条件下回转水动力的准确预报。预报结果与试验值相比，偏差在12%以内，具有一定的工程实用价值。本文提出的向心力源项空间离散误差修正方法，在一定简化条件下，通过物理对象的转换，较好地实现了水动力计算中对由旋转流动静压场产生的向心力修正，避免对控制方程中源项的复杂数学处理，具有一定的新颖性。

SST k－ω湍流模型对回转运动复杂流动具有良好的表征能力，但对首尾段强烈分离流动的捕捉能力仍不足。对较大无因次回转角速度条件下，流动非定常特性的描绘，还需要具有更优良性能的湍流模型来实现。

将本文数值预报方法进一步推广应用，可以对操纵运动方程中与定常回转运动相关的角速度导数、舵角导数及系列耦合水动力导数进行预报，最终建立水下航行体操纵性水动力的数值回转水池。

［1］ FELDMAN J.Straight－line and rotating arm captive－model experiments to investigate the stability and control characteristics of submarine and other submerged vehicles［R］.Taylor Naval Ship Research and Development Center，1989.

［2］ SUNG C，FU T，GRIFFIN M，et al.Validation of incompressible flow computation of forces and moments on axisymmetric bodies undergoing constant radius turning［C］//Twenty－First Symposium on Naval Hydrodynamics，1997.

［3］ SUNG C H，JIANG M Y，RHEE B，et al.Validation of the flow around a turning submarine ［C］//The Twenty－Fourth Symposium on Naval Hydrodynamics， Fukuoka， Japan：2002.

［4］ FU T C，ATSAVAPRANEE P，HESSD E.Measurements of the c ross－f low w ake of a t urning s ubmarine m odel（ONR Body－1）［C］//25th ONR Symposium on Naval Hydrodynamics，Fukuoka， Japan，2002.

［5］ ATSAVAPRANEE P，FORLINIT，FUREY D.Experimental m easurements for CFD v alidation of the f low about a s ubmarine m odel （ONR Body－1）［C］//25th ONR Symposium on Naval Hydrodynamics，St John’s，Canada，2004.

［6］ ETEBARI A， ATSAVAPRANEE P， et al.Experimental measurements on a SUBOFFmodel in a turningmaneuver［C］//27th Symposium on Naval Hydrodynamics，Seoul，Korea，2008.

［7］ 柏铁朝，梁中刚，周轶美，等.潜艇操纵性水动力数值计算中湍流模式的比较与运用［J］.中国舰船研究，2010，5（2）： 22－28.

BAIT C，LIANG Z G，ZHOU Y M，et al.Comparison and a pplication of t urbulence m odes in s ubmarine m aneuvering h ydrodynamic f orces c omputation ［J］.Chinese Journal of Ship Research，2010，5（2）： 22－28.

［8］ BRIDGES D H.A detailed study of the flow field of a submarine at large angle drifts ［R］.US：Office of Naval Research， 2001.

Numerical Calculation on Hydrodynam ic Performance of the Submerged Vehicle in Turning Motion

Lu Jin－guo Liang Zhong－gang Wu Fang－liang Zhou Yi－mei
China Ship Developmentand Design Center，Wuhan 430064，China

A num erical procedure based on solving the incompressible RANSequations in a steady rotating coordinate system was developed for the prediction of forces and moments on submerged vehicle in turningmotion with the commercial code FLUENT.New method considered numerical dissipation was developed for the correction of centripetal force，and several key elements influenced the solution were also taken into account.In the case of Re＝11.7×106， themaximum relative deviation between the calculation results and measured data is less than 12%.The good agreement proves that the method is practical for numerical prediction of the hydrodynamic performance on the submerged vehicle in turningmotion.

turnin gmotion； hydrodynamics； centripetal force correction； FLUENT； submerged vehicle

U661.3

A

1673－3185（2011）06－08－05

10.3969/j.issn.1673-3185.2011.06.002

2011-06-28

卢锦国（1986－），男，硕士研究生。研究方向：潜艇操纵性水动力数值计算。E-mail：ljgfatherland＠126.com

梁中刚（1966－），男，研究员，硕士生导师。研究方向：船舶总体设计。E-mail：lzg701papers＠126.com