变加速运动圆球附加质量和阻力仿真分析

姚保太+康宁+郑伟奇

作者简介： 姚保太（1973—），男，河南获嘉人，高级工程师，研究方向为机械制造及仿真计算，(Email)yaobaotai@163.com

通信作者： 康宁（1963—），女，辽宁金州人，副教授，博士，研究领域方向为流体仿真计算，(Email)kangning@buaa.edu.cn0引言

附加质量和阻力计算是水下物体运动分析的关键问题.［1］以往对附加质量的计算主要是根据势流理论［14］，如在无界无黏流体中球体等简单形状物体的附加质量存在解析解,复杂形状物体的附加质量通过估算法或三维面元法求得.近年来，随着计算流体力学的发展，开始用数值模拟的方法计算物体的附加质量，发展基于动网格技术计算附加质量的新方法.［5］数值方法有更广泛的应用空间，不仅可以计算理想流体中的附加质量［56］，还可以计算黏性流体中的附加质量［710］.

人们已经熟知无界无黏流体中物体加速运动时的附加质量只与物体形状及流体密度有关，与运动加速度无关，目前还缺少系统地研究在有界区域中流体黏性和加速度对做变加速运动物体附加质量的影响.由于工程中实际问题的复杂性，比如水下导弹发射问题等，因此有必要研究在有界区域黏性流体中物体做变加速运动的附加质量问题.

本文采用数值模拟方法和动网格技术，计算管体域中黏性不可压流体中圆球做变加速运动时的附加质量、阻力和阻力因数，深入分析和研究在管体域中圆球变加速运动对附加质量、阻力和阻力因数的影响.

1计算模型和计算方法

1.1控制方程组

水中运动物体满足非定常不可压NS方程，

连续方程： vixi=0(1)动量方程：vit+vjvixj=Fi－1ρpxi+μρxjvixj(2)式中：vi和vj为速度；F为质量力；p为压力；ρ为流体密度；μ为流体动力黏性系数.

1.2计算域、边界条件和网格划分

1.2.1管体域

管体域由管体、喇叭口和水域组成，圆球直径为0.254 m，圆球在管体域中运动前的计算域见图1.通海水口边界条件为压强进口，出口为压强出口，压强均为57 m水深处压强574 826.8 Pa；根据圆球平移运动特点，动网格模型选择动态网格层变方法，圆球表面为无滑移运动壁面条件；管体和喇叭口为无滑移静止壁面条件，水域边界为滑移静止壁面条件.圆球在管体域中运动前、后网格分布见图2和3.采用六面体结构网格，网格数约为58万个.

图 1管体域中圆球运动前的计算域

Fig.1Computation domain in tube domain before ball

movement

图 2管体域中圆球运动前的网格分布

Fig.2Mesh distribution of tube domain before ball movement

图 3管体域中圆球运动后的网格分布

Fig.3Mesh distribution of tube domain after ball movement

1.2.2小计算域

作为对比，对圆球在小计算域中的运动情况进行计算.小计算域为圆柱形，取较大尺度，直径为3.6 m，长度为14 m.纵向对称面网格分布见图4.六面体结构网格，网格数约为139万个.计算域左侧为压强入口边界条件，右侧为压强出口条件，相对压强为0；其他边界和动网格模型同管体域.

图 4小计算域纵向对称面的网格分布

Fig.4Mesh distribution in longitudinal symmetry plane of

small computation domain

1.2.3大计算域

为验证附加质量仿真方法的可靠性，对圆球在大计算域中的运动情况进行计算.大计算域仍为圆柱形，直径加大到14.4 m，长度加大到20.86 m.采用六面体结构网格，网格数约为264万个.边界条件和动网格模型同管体域.1.3圆球变速运动参数

1.3.1先加速后减速

为模拟发射体在管体内直至出管的运动情况，对圆球先加速后减速运动情况进行计算.加速阶段模拟发射时在主动力作用下的运动情况，减速阶段模拟在惯性力作用下直至出管的运动情况，具体参数见表1.

表 1圆球先加速后减速的运动参数

Tab.1Movement parameters when ball is accelerated motion first and then decelerated初速度/

(m/s)加速度/

(m/s2)时间/s末速度/

(m/s)位移/

m加速阶段052.0000.32917.1002.814减速阶段17.108-5.6690.39014.8976.241总和0.71931.8979.055

1.3.2一直加速

作为对比，对圆球在管体域内一直加速运动工况进行计算，加速度和初速度与先加速后减速的加速阶段相同，位移与先加速后减速的总位移相近，见表2.

表 2圆球一直加速的运动参数

Tab.2Movement parameters when ball is accelerated motion all the time初速度/

(m/s)加速度/

(m/s2)时间/s末速度/

(m/s)位移/m052.0000.58930.6209.020

1.4求解方法和计算参数

速度和压强耦合方法采用SIMPLE方法，压力修正方程离散格式采用Standard格式，对流项离散格式采用2阶迎风格式，湍流模型为RNG kε湍流模型，近壁区域处理采用非平衡壁面函数.

海水密度为1 028 kg/m3，动力黏性系数为1.394×10-3 kg/(m•s)，参考压强为1标准大气压101 325 Pa.计算初始速度为0，初始压强为57 m水深处压强574 826.8 Pa.

2仿真方法验证

2.1阻力

分2种情况对圆球阻力因数进行验证，分别为水绕静止圆球运动和圆球在静水中运动.

圆球实验阻力因数的范围［3］为0～106，根据定义可得圆球速度范围为0～5.339 m/s，因此对最大阻力因数106和较小阻力因数104，其分别对应最大速度5.339 m/s和较小速度0.053 4 m/s的情况进行计算.采用小计算域进行计算.

2.1.1静止圆球

计算域左侧为速度入口边界条件，速度大小分别为5.339 0和0.053 4 m/s；右侧为压强出口边界条件，相对压强为0.

圆球表面为无滑移静止壁面边界条件；除进口、出口和圆球表面以外的边界条件为滑移静止壁面边界条件.计算结果见表3.

表 3圆球静止时的阻力因数

Tab.3Drag factor of stationary ball速度/(m/s)实验值计算值误差/%0.053 40.400.4164.005.339 00.170.1652.94

由表3可知,在2种速度情况下计算的圆球阻力因数与实验结果的误差都在4%以内,说明圆球静止时阻力计算结果可靠.

2.1.2运动圆球

计算域左侧为压强入口条件，右侧为压强出口条件，相对压强为0.动网格模型及圆球表面边界条件均同管体域中的运动圆球，速度大小分别为5.339 0和0.053 4 m/s；其余边界为滑移静止壁面条件；初始速度和压强均为0.计算结果见表4.

表 4圆球运动时的阻力因数

Tab.4Drag facctor of moving ball速度/(m/s)实验值计算值误差/%0.053 40.400.3717.255.339 00.170.18710.00

由表4可知,在2种速度情况下计算的圆球阻力因数与实验结果的误差都在10%以内,说明圆球运动时阻力计算结果可靠.

2.2附加质量

2.2.1附加质量计算公式

根据文献［2］，附加质量λ=ρv2iViv2o(3)式中：vi为第i控制体的流体速度；Vi为第i控制体的体积；vo为物体运动速度.

2.2.2圆球附加质量的理论解

根据文献［1］，在无界无黏流体情况下，圆球直线加速运动附加质量的理论解λ=23πρR3 (4)式中：R为圆球半径.

圆球直径、水密度和动力黏性因数同前，这时圆球附加质量的理论解为4.41 kg.

2.2.3圆球附加质量的计算结果

为验证附加质量仿真方法的可靠性，对小计算域中圆球在无黏流体中做不同加速度的运动进行计算，圆球的初速度均为0，加速时间为1.1 s，加速度分别为0.01，0.10，1.00和5.00 m/s2，并与圆球附加质量的理论解进行对比.所采用的网格、边界条件、初始条件、计算设置和动网格模型均与运动圆球阻力验证时一样.计算的附加质量见图5.

图 5小计算域中加速度对圆球附加质量的影响

Fig.5Effect of acceleration on added mass of ball in

small computation domain

由图5可知，时间和加速度越大附加质量越大.在加速度较小时（0.01，0.10和1.00 m/s2）计算的附加质量与理论解较为接近；在加速度较大时（5.00 m/s2），运动初期计算的附加质量与理论解较接近，随着时间增大其值开始大于理论解.这是由于随着加速度和时间的增大，运动速度增大，导致计算域边界效应的影响也随之增大，附加质量增加，因而与理论解的差别增大.

为进一步研究计算域边界效应的影响，在大计算域中计算加速度为5.00 m/s2的工况，初速度和加速时间均同前文，无黏流体中的附加质量的计算结果见图6，并与之前小计算域结果进行对比.由图6可知，除运动初期，无黏流体中小计算域附加质量均大于大计算域附加质量，大计算域附加质量均大于理论解，并且随着时间增加速度加大，差别也越大.大计算域附加质量与理论解较接近，小计算域附加质量与理论解差别较大.由图5和6的结果可知，加速度越小或计算域越大，边界效应越小，则计算的附加质量越接近理论解.图 6在a=5 m/s2时无黏流体中计算域尺度对圆球

附加质量的影响

Fig.6Effect of computation domain scale on added mass of

ball at a=5 m/s2 in inviscid fluid

2.2.4椭球体附加质量的实验验证

为验证附加质量的理论计算公式（3），将黏性流体中加速运动椭球体附加质量的仿真结果与文献［6］的实验结果进行对比.椭球长轴为0.4 m，短轴为0.2 m，与文献［6］中的完全一样.计算域为圆柱形，直径为6 m，长度为14 m.计算时水温度为5 ℃，密度为1 000 kg/m3，动力黏性系数为1.518×10-3 kg/(m•s).椭球体为六面体网格，网格数为106万个，网格分布见图7.

图 7椭球体纵向对称面的网格分布图

Fig.7Mesh distribution in longitudinal symmetry

plane of ellipsoid

边界条件、初始条件和其他设置均同阻力验证时的运动圆球，椭球动网格模型同运动圆球.椭球体运动初速度为0，加速度为1 m/s2，加速时间分别为0，0.5和0.9 s.根据式（3）计算的附加质量和与文献［6］的对比见表5.由表5可知，计算的附加质量与实验结果比较接近，误差在9%以内，说明附加质量的计算结果基本正确.

表 5椭球体的附加质量

Tab.5Added mass of ellipsoid加速时间/s实验值/kg计算值/kg误差/%0.11.701.742.350.51.701.774.120.91.701.858.82

综合前文分析结果可知，附加质量的仿真方法较为可靠.3计算结果及分析

对圆球在管体域黏性流体中先加速后减速和一直加速2种变速运动的附加质量、阻力和阻力因数进行仿真计算，并对小计算域黏性流体中圆球先加速后减速运动进行仿真计算，同时给出圆球在无界无黏流体中附加质量的理论解.

管体域和小计算域中圆球在黏性流体中附加质量、阻力和阻力因数随位移的变化情况见图8～10，在管体域中4 m处圆球先加速后减速和一直加速的速度矢量图和压强云图见图11和12.

图 8管体域和小计算域中圆球的附加质量

Fig.8Added masses of ball in tube domain and small

computation domain

图 9管体域和小计算域中圆球的阻力

Fig.9Drags of ball in tube domain and small computation

domain

图 10管体域和小计算域中圆球的阻力因数

Fig.10Drag factors of ball in tube domain and small

computation domain(a)速度矢量图,m/s

(b)压强云图,Pa

图 11在管体域中4 m处先加速后减速圆球的计算结果

Fig.11Calculation results at 4 m in tube domain while ball is accelerated first and then decelerated

(a)速度矢量图,m/s

(b)压强云图,Pa

图 12在管体域中4 m处一直加速圆球的计算结果

Fig.12Calculation results at 4 m in tube domain while ball is accelerated all the time

由图8可知，在运动开始时，小计算域中附加质量仿真值与理论解4.41 kg吻合较好；随着运动进行仿真值开始大于理论解，最后达到20.8 kg.管体域中附加质量均大于小计算域中附加质量，并随着运动进行差别逐渐增大，在出管体（7 m）后达到最大值，然后开始变小.先加速后减速运动的最大值为58.6 kg，最终降为48.1 kg；一直加速运动的最大值为44.7 kg，最终降为33.9 kg.由于管体域尺寸小于小计算域，边界效应较大，因而附加质量较大，出喇叭口后过流断面变大，边界效应减小，导致附加质量变小.管体域中圆球先加速后减速和一直加速的附加质量在2.8 m之前完全一样，因为该阶段2种工况运动状态完全一样；在2.8 m之后，一直加速运动附加质量仍保持连续增加态势，而先加速后减速运动附加质量在2.8 m处出现不连续现象，开始大于一直加速附加质量.先加速后减速运动由加速状态变为减速状态时，运动状态的改变远大于一直加速运动状态.由图11和12可看出，在4 m处，先加速后减速圆球周围速度虽然小于一直加速的速度，但尾迹流影响范围明显较大，在圆球尾部较远处还卷起较大尾涡，因而导致附加质量非但没有变小，反而增大.小计算域附加质量在加速阶段结束（2.8 m）时也出现不连续现象，只是由于域较大，现象没有管体域的明显.

由图9可知，先加速后减速运动圆球在管体域中的阻力大于小计算域中的阻力.随着运动进行，阻力不断增大，在加速阶段结束（2.8 m）后达到最大，管体域中达5 359.5 N，小计算域中达1 976.3 N，然后出现突然减小而后又增大的震荡现象.随着速度减小阻力也变小，最终管体域中阻力减小到1 354.9 N，小计算域中减小到1 268.8 N.

由文献［3］可知，无界黏性流体中加速运动的阻力F非取决于附加质量、加速度和匀速运动时的阻力F定，即F非=λa+F定 (5)阻力中一部分为匀速运动阻力，而该部分阻力与速度平方成正比，因此速度越小，匀速运动阻力越小.在加速阶段，速度和附加质量都增加，因此阻力增加；在减速阶段，速度减小，附加质量增加，但增加得不太多，最终使得阻力变小；加速阶段转换为减速阶段时出现震荡的原因是由于圆球运动加速度发生突变，导致流体运动也发生剧变，因而阻力发生突变.管体域中圆球在出管体处（7 m）阻力减小程度加剧，这是由于圆球进入喇叭口后，过流断面开始变大，因此除速度变小的因素外，又增加一个降低阻力的因素.

在2.8 m之前,管体域中圆球一直加速运动与先加速后减速运动的阻力相同；在2.8 m之后,一直加速运动的速度继续增大，因而阻力继续增加,直到圆球在出管体处（7 m）达最大值11 238.7 N，之后阻力开始变小，最终减小到6 832.4 N.

由图11和12可知，在4 m处先加速后减速圆球前后压差小于一直加速的压差，因此阻力较小.

由图10可知，管体域中圆球阻力因数均随着运动的进行不断减小，都大于小计算域中的阻力因数.在运动初期速度较小，阻力因数与速度平方成反比，因此阻力因数较大.管体域中圆球先加速后减速阻力因数最终降为0.24，一直加速时最终降为0.28，小计算域中圆球先加速后减速时其最终降为0.22.管体域中圆球先加速后减速与一直加速运动的阻力因数在2.8 m之后也比较接近，原因是虽然一直加速的阻力大于先加速后减速,但速度也较大.

4结论

管体域中圆球的附加质量均大于无限域中的附加质量，并随着运动的进行差别增大，在出管体后达到最大值，然后开始变小.管体域中圆球先加速后减速与一直加速的附加质量在加速阶段完全相同，减速阶段先加速后减速的附加质量大于一直加速的附加质量.

管体域中先加速后减速的圆球阻力大于无限域中的阻力.随着运动的进行阻力均不断增大，并在加速阶段结束后达到最大.随着速度的减小阻力也变小.管体域中圆球在出管体处阻力减小的程度加剧.在加速阶段,管体域中圆球一直加速运动与先加速后减速运动的阻力相同；在减速阶段,一直加速运动的阻力继续增加到圆球出管体处达最大，之后阻力开始变小.

管体域中圆球的阻力因数均随着运动的进行不断减小，在运动初期阻力因数较大.管体域中圆球在出管体处阻力因数减小的程度加剧.管体域中的阻力因数比较接近，都大于无限域中的阻力因数.

以上结论说明运动区域、加速度和流体黏性对附加质量、阻力和阻力因数均有较大的影响.参考文献：

［1］王献孚, 熊鳌魁. 高等流体力学［M］. 武汉: 华中科技大学出版社, 2003: 200.

［2］夏国泽. 船舶流体力学［M］. 武汉: 华中科技大学出版社, 2003: 117118.

［3］潘文全. 流体力学基础［M］. 北京: 机械工业出版社, 1982: 288293.

［4］李宝元, 姜昱汐. 船舶在进出船厢运动中的附加质量的计算［J］. 计算力学学报, 2003, 20(3): 350354.

LI Baoyuan, JIANG Yuxi. Calculation of added mass for ship entering and leaving a shipbox［J］. Chin J Comput Mech, 2003, 20(3): 350354.

［5］马烨, 单雪雄. 数值计算复杂外形物体附加质量的新方法［J］.