水下航行体流体阻尼力系数的cFd计算研究

尤天庆，王占莹，程少华，穆 洲，鲍文春

（北京宇航系统工程研究所，北京，100076）

水下航行体流体阻尼力系数的cFd计算研究

尤天庆，王占莹，程少华，穆 洲，鲍文春

（北京宇航系统工程研究所，北京，100076）

水下航行体流体阻尼力关系到弹道的稳定性和操纵性，当航行体以较大角速度机动时阻尼力更是不可忽略，因此准确预示航行体流体阻尼力具有重要意义。以计算流体动力学（computational Fluid dynamic，cFd）方法为基础，结合动网格和动坐标系技术，进行基于雷诺平均N-S方程的阻尼力计算方法研究，发展针对水下航行体附加力系数和俯仰阻尼力矩系数的数值计算方法，通过与试验数据对比表明，该方法具有较好精度，在此基础上分析了空化现象对流体阻尼力系数的影响。

计算流体动力学；水下航行体；流体阻尼力；空化

0 引 言

潜射导弹、鱼雷、水雷和潜艇等水下航行体在水下运动过程中所受的流体力为弹道或控制系统的主要输入参数，关系到航行体弹道稳定性、操控性或航程等总体参数，是各类型航行体设计的关键参数。航行体所受流体动力主要取决于流场性质、航行体流体动力外形以及相对流体的运动，迄今为止，还难以用理论计算或模型试验的方法直接得到包含上述所有因素的总流体动力，只能在各种假设和简化下，把流体动力分解成若干部分，分别由计算或试验的手段获得，然后综合成总流体动力[1]。通常在航行体设计过程中，将流体动力分为位置力、阻尼力和惯性力3个部分。航行体以一定速度作平移运动的同时，还作旋转运动，由旋转运动而产生的流体力增量为阻尼力[2]。流体动力中阻尼力（如俯仰阻尼力矩、附加力等）关系到航行体弹道稳定性及操纵性[3，4]，当旋转运动角速度较大时阻尼力更是不可忽略，因此阻尼力计算在水动力设计中有非常重要的作用。

目前，无论是气动力还是水动力，针对流体阻尼力的计算主要有工程估算、试验以及计算流体力学（computational Fluid dynamic，cFd）等方法。其中，工程估算方法计算速度快但精度较差，只能作为初步估算[5]；试验所得的数据较精确，但周期长、费用高；cFd计算方法比势流理论更具从理论到工程实践的可扩展性，比试验方法更容易实现，成本低。

近年来，随着cFd及计算机技术发展，基于N-S方程求解流体阻尼力取得一定的发展。在气动力计算方面，分别采用稳态和非稳态求解方法计算飞行器在空气中飞行的流体阻尼力[6]。其中，稳态求解的方法即通过改造N-S方程，在方程中加入惯性力源项以求解匀速圆周运动飞行器受力情况，得到飞行器的定常拉升运动中阻尼力[7，8]；非稳态求解即计算航行体强迫振荡运动中的俯仰力矩迟滞曲线，通过辨识得到计算模型的俯仰稳定性参数[9]。不同于气动力，当水下航行体高速运动或通气时，航行体表面会附着有空泡，空泡发展过程带有很强的非定常和非线性特性，这同样给流体阻尼力的确定造成很大困难。由于水的密度较大，航行体水中运动的流体惯性力（附连水质量）不可忽略，计算或试验的阻尼力系数会有附连水质量的影响。

本文基于通用cFd计算程序，运用动坐标系稳态计算求解流体阻尼力，并通过与试验结果对比验证了计算方法的有效性。在此基础上对数值计算结果进行讨论，分析了空化现象对流体阻尼力系数的影响。

1 数值计算方法

1.1 阻尼力计算方法

航行体所受的流体阻尼力系数与坐标系选取紧密相关，因此首先介绍本体坐标系的定义。本文航行体本体坐标系与航行体固连，原点选在航行体质心，Ox1轴与航行体纵轴重和指向前方，Oy1轴位于航行体纵对称面内与Ox1轴垂直指向上方，Oz1轴与Ox1y1平面垂直并与Ox1轴和Oy1轴构成右手系，如图1所示。

图1 本体坐标系

无论采用试验还是数值计算手段获得航行体流体阻尼力，其基本思路是构造一种运动方式，使试验测得或数值计算得到航行体所受的流体力包含所关注的流体力，并能通过数据处理的手段得到所关注的流体力。这里构造的运动方式为匀速圆周运动，数值模拟实现方案为采用旋转坐标系的方法。

令航行体以攻角α绕指定圆心做匀速圆周运动，航行体质心与圆周运动圆心距离为R，转动角速度为ωz1，则匀速圆周运动线速度为V=Rωz1。在本体坐标系下，航行体所受到的流体力有如下表达式：

式中 q为动压；S，l分别为参考面积和参考长度；N为法向力；λ为附连水质量；M为俯仰力矩。

通过航行体定速圆周运动数值计算，可求得指定攻角α和无量纲旋转角速度1下的法向力N和俯仰力矩M。式（1）和式（2）等号右侧变量均为已知，由此可求得1和。

由于λ11V ωz1cosα项和λ26V ωz1cosα 项内均含有ωz1，因此在弹道计算中通常将含有λ11的λ11Vωz1cosα项记入内，将含有λ26的λ26Vωz1cosα 项记入1内。1和的计算公式为

数值计算控制方程为基于雷诺平均N-S方程的质量守恒方程、动量守恒方程和kε-湍流模型方程。该计算流场建立在旋转坐标系上，在此坐标系下，网格在计算中保持静止。其中：

质量守恒方程为

动量守恒方程为式中rU为相对旋转坐标系的速度矢量；ω为旋转坐标系旋转角速度矢量；r为流场质点在旋转坐标系中的位置矢量；p为压力；t为时间；effμ为有效粘度。

本文采用标准kε-湍流模型来建立流场中时均统计量和脉动量之间的关系。kε-模型是目前比较通用的一种湍流模型，在计算上易于处理，鲁棒性也较强。

本文数值离散方法采用有限体积法，湍流输运方程采用一阶差分格式，时间离散为一阶隐式。压力速度耦合方式采用SIMpLe算法，压力及动量方程离散采用二阶迎风格式。

1.2 空化模型

航行体在水下高速航行时，周围流体与航行体相对运动，当航行体表面压力下降至饱和蒸汽压力以下并维持一定时间后，空泡将产生。流场中是否存在空泡，依赖于流动条件和航行体的几何形状。在空化现象中，通常用无量纲参数σ来表征自然空化状态。

应力张量τ与应变率关系为

式中 p∞为流场中未受扰动处的压力；vp为水的饱和蒸汽压；V为航行体运动速度。

数值计算采用基于输运方程的空化模型。空化考虑为水与水蒸汽之间的质量转换，多相流中的输运方程为

式中 Γib为b相进入i相的单位体积质量源（本文中b相为水，i相为水蒸汽），Γib=m˙ibAib，其中，b为b到i相单位相间交界面质量流率；Aib为相间界面面积。

根据空泡动力学可知，气泡壁面的运动方程为

式中BR为气泡半径；τ为表面张力；p为流体参考压力；lρ为液体密度。忽略式（10）中表面张力项以及高次项，化简后得：

气体体积含量为γb=VBNB=4πNB/3，其中NB，VB分别为单位体积气泡数和气泡体积。

水和水蒸汽间单位体积总质量传递率为

式中 F为实验获得的经验系数；gρ为水蒸汽密度。

2 数值计算结果讨论

2.1 数值计算结果的有效性

为验证阻尼力计算方法的准确程度，计算一航行体在不同无量纲旋转角速度ωz1条件下的法向力N及俯仰力矩M，图2为航行体周围流场速度云图。通过对计算数据处理得到附加力系数Cωz1和俯仰阻尼力矩

。数值计算结果与实验数据无量纲化对比结果如表1所示。

图2 α=2˚1=0.46航行体周围流场速度云图

表1 和m数值计算结果与试验对比

表1 和m数值计算结果与试验对比

系数 攻角/(°) 数值计算 试验结果 误差0 0.660 0.628 5% Cωz1y12 0.655 0.618 6% 6 0.527 0.500 5% 0 -0.376 -0.385 2% mωz1z12 -0.394 -0.385 2% 6 -0.380 -0.394 4%

表1中的数值计算结果表明，所采用的流体阻尼力数值计算方法具有较高的精度。

以往研究结果表明圆柱航行体长径比，对航行体阻尼力有较大的影响，并且其影响近似为线性[1]。为验证数值计算结果规律的正确性，在保持质心相对位置cg/xL和航行体头部外形不变的情况下，计算了不同长径比/L d航行体的附加力系数和俯仰阻尼力矩系数（见图3、图4）。数值计算结果表明随着长径比的增大，附加力系数和俯仰阻尼力矩系数均逐渐增大。其中附加力系数增加幅度较小，俯仰阻尼力矩系数增幅较大，长径比/L d在7.1～14.2范围内逐渐增加时，附加力系数在0.905～0.951内缓慢递增，增幅约为5%，俯仰阻尼力矩系数绝对值在0.453～0.573内近似程线性逐渐递增的趋势，增幅约为26%。

图3 不同长径比的附加力系数

图4 不同长径比的俯仰阻尼力矩系数

2.2 空化现象的影响

航行体表面空泡发生在初生空化临界状态之后，水流从物面脱流，在壁面附近形成相对稳定的气水交接面。但空泡发展具有一定的非定常特性，即在来流稳定条件下，空泡形态也会存在一定的波动，同时航行体的受力也存在一定的波动。在此忽略这种空泡脉动的影响，采用均质多相流模型结合空化模型，分析了不同空化数条件下的附加力系数和俯仰阻尼力矩系数。

在数值计算中采用减小环境压力的方法，使航行体自然空化数逐渐降低，自然空泡覆盖的区域逐渐增大。在空泡覆盖航行体表面的情况下，空泡覆盖区域内航行体表面压力为恒定，空泡末端与航行体壁面接触时，水流损失径向速度形成压力峰值，如图5所示。

图5 0.24σ=头部空泡形态及周围流场压力云图

图6 不同空化条件附加力系数分布

图7 不同空化条件的俯仰阻尼力矩系数分布

针对空泡覆盖区域，由于空泡内部压力恒定，空泡覆盖区域内不存在分布法向力，且随着旋转角速度的增加法向力不会发生变化，因此在空泡覆盖区域，分布附加力系数和俯仰阻尼力矩系数均为零，在图6和图7中表现为分布阻尼力系数曲线的平直段。

针对回射压力影响区域，航行体圆周运动过程中，头部空泡覆盖区域存在局部攻角，使空泡末端闭合处迎、背水面空泡回射压力存在明显不同，迎水面压力较背水面压力峰值在位置上略靠前，但空间上分布范围较窄。因此在空泡末端，沿航行体轴线方向分布法向力在空间上先后出现负正两个峰值。随着旋转角速度变化，空泡区的局部攻角不断发生变化，空泡迎背水面回射压力的位置和强度随之而发生变化，即空泡末端的负正法向力峰随着旋转角速度变化而大幅度变化，因此在此区域的分布附加力系数和俯仰阻尼力矩系数也先后存在负正两个峰值。这种阻尼力系数分布情况，在图6和图7中体现为曲线先下降而后缓慢上升，最终与全湿流曲线变化规律一致。

随着空泡数降低，空泡覆盖的长度逐渐增大，表现为图6和图7中，阻尼力系数曲线平直段逐渐增大。不同空泡数条件下，空泡闭合区位置发生移动，相应的回射区局部攻角发生变化，从而导致分布流体阻尼力系数正负峰值存在差异。以上两方面因素共同作用，使得航行体整体的附加力系数和俯仰阻尼力矩系数随着空泡数降低而增大，如图8和图9所示。

图8 不同空化条件的附加力系数

图9 不同空化条件的俯仰阻尼力矩系数

3 结 论

本文基于cFd计算，利用动坐标系技术，发展水下航行体流体阻尼力的数值计算方法，在此基础上，研究了空化现象的影响，得到如下结论：

a）基于旋转坐标下雷诺平均N-S方程，求解水下航行体流体阻尼力系数，通过与试验数据以及以往规律性认识对比验证计算方法的准确性；本文数值计算方法能考虑流动非线性、涡以及粘性等因素，同时适用于任何几何外形的航行体，其对于工程问题更具有适用性和可扩展性。

b）针对空化现象的影响，在空泡内部，压力恒定并且随着旋转角速度增加分布法向力保持不变，分布流体阻尼力系数为零；在回射压力影响区域，分布阻尼力系数沿航行体轴线先下降而后缓慢上升；随着空泡数减小，空泡覆盖长度增加，航行体整体的附加力系数增大，俯仰阻尼力矩系数绝对值减小。

[1] 荣建徳. 水下运载器性能的分析与设计[M]. 北京: 国防工业出版社, 2008.

[2] 张宇文. 鱼雷总体设计原理与方法[M]. 西安: 西北工业大学出版社, 1998.

[3] 蒋胜矩, 刘玉琴, 党明利. 基于定常N-S方程的飞行器滚转阻尼力矩系数导数计算方法[J]. 弹箭与制导学报, 2008, 28(1): 180-182.

[4] Sen d. a study on sensitivity of maneuverability performance on the hydrodynamic coefficients for submerged bodies[J]. Journal of Ship Research, 2000, 44(3), 186-196.

[5] danberg J e, Weinacht p. approximate computation of pitch-damping coefficients[c]. Monterey: aIaa atmospheric Flight Mechanics conference and exhibit, 2002.

[6] park S H, Kwon J H. comparisons of steady and unsteady methods for pitch-damping predictions[c]. Orlando: 21stapplied aerodynamics conference, 2003.

[7] Weinacht p, Sturek W B, Schiff L B. Navier-stokes predictions of pitch-damping for axisymmetric shell using steady coning motion[R]. ad-a285866, 1994

[8] Weinacht p. Navier-stokes predictions of the individual components of the pitch-damping coefficient sum[R]. aRL-TR-3169, 2004.

[9] 陶洋, 等. 基于cFd的带控制舵导弹的动导数数计算[J]. 航空动力学报, 2010, 25(1): 102-105.

[10] Knapp R T, daily J W, Hammitt F G. cavitation[M]. New York: Mcgraw-Hill Book company, 1970.

[11] Kunz R F, Boger d a, chyczewski T S, et al. Multi-phase cFd analysis of natural and ventilated cavitation about submerged bodies[c]. San Francisco: Third aSMe/JSMe Joint Fluids engineering conference, 1999.

[12] Knapp R T, daily J W, Hammitt F G. cavitation[M]. New York: Mcgraw-Hill Book company, 1970.

CFd Research on Underwater Vehicle Hydrodynamic damping Force Coefficient

You Tian-qing, Wang Zhan-ying, cheng Shao-hua, Mu Zhou, Bao Wen-chun
(Beijing Institute of aerospace System engineering, Beijing, 100076)

Underwater vehicle hydrodynamic damping force is crucial to trajectory stability and maneuverability, especially for the vehicle with high angle velocity, thus the accurate prediction of hydrodynamic damping force is important. Based on the cFd method, Reynolds averaged N-S equations based hydrodynamic damping force calculation research is conducted. combining with dynamic mesh and rotating reference frames, pitch damping Force coefficient calculation methods are introduced. compared with experiment result, the validity and accuracy is proofed. Furthermore, the effect of cavitation is discussed.

computational Fluid dynamics; Underwater vehicle; Hydrodynamic damping force; cavitation

O351.2

a

1004-7182(2016)01-0066-05

10.7654/j.issn.1004-7182.20160115

2015-04-23；

2015-06-02

水动力学教育部重点试验室资助项目

尤天庆（1984-），男，工程师，主要研究方向为计算流体力学