圆锥前体大攻角绕流的数值计算与分析

王中一，郑博睿，高 超，方 洪，熊俊涛，刘 锋，罗时钧

(1西北工业大学翼型/叶栅空气动力学国防科技重点实验室，西安 710072;2北京临近空间飞行器系统工程研究所，北京 100076;3美国加州大学尓湾分校，美国加州尔湾 92697-3975)

0 引言

现代飞行器前体和导弹头部，通常采用细长尖头旋成体的设计外形。在大攻角下，细长旋成体背风面会产生一对分离涡，当攻角增大到一定程度时，原本对称的分离涡突然变得非对称，同时伴随有方向和大小均无法预估的侧向力和力矩，这对飞行器的操纵性和稳定性有很大影响[1－2]。Keener等发现侧力大小和方向对顶点处的几何外形很敏感，受雷诺数或马赫数的影响并不大[3]。研究发现，细长前体分离涡对尖头处小的扰动非常敏感[4]，而尖鼻点在形状上近似于圆锥，当地流场相当于一个正切圆锥处的流场，所以可以通过研究圆锥体绕流来了解细长体上非对称流场的特征。杨宇等[5]对圆锥前体流动进行了数值模拟，小于25°攻角得到的计算结果与实验对比较好，但大攻角时差异较大。文中在文献[5]的基础上进一步优化网格，并引入了Gamma theta转捩模型进行数值模拟。

1 数值方法

根据研究内容，在CFX中采用3个计算模型进行分析，分别是层流模型，Gamma theta转捩模型和SST湍流模型。其中层流模型以非定常N-S方程作为控制方程，只适合于层流和低雷诺数的情况;Gamma theta转捩模型多用于转捩预测，它集合了转捩经验关系式和低雷诺数湍流模型的优势，通过求解两个变量(间歇因子和动量厚度雷诺数)的标准输运方程来实现该目的;SST湍流模型的控制方程为RANS方程，该模型集合了k-ε模型可较好的模拟充分发展湍流流动的优点和k-ω模型可广泛应用于各种压力梯度下的边界层问题的优点，在求解边界层流动时有很高精度，但其对转捩的预测受到质疑，因为该模型中阻尼函数的标定依据是再现粘性底层的行为，而未考虑转捩的物理机理。

计算初始流场为均匀来流，来流速度为10～70m/s，攻角为15°～35°，基于圆锥段前体底面直径的雷诺数为 0.1 ×106～0.7 ×106，壁面光滑且满足无滑移条件，计算分析类型采用定常态。计算结果与Meng XS[6]的实验数据进行了对比。计算模型与文献[6]尺寸一致，采用圆锥圆柱组合体模型，依次由圆锥段、过渡段和圆柱段三部分组成。圆锥半顶角为10°，为避免数值计算时尾迹流动对圆锥前体流场造成影响，圆柱段延伸至远场边界。圆锥尖端第一个轴向网格长度是圆锥段长度的0.1%，远离物面的第一个径向网格高度是0.001mm，为模型圆柱段半径的10－5倍，径向远场边界距轴线的距离是圆柱段半径的40倍。实验模型上的测压孔分布在沿轴线方向的9个截面上。由于圆锥绕流在相邻截面的压力分布差异很小，而第一和第二截面处的测压孔较少，不能精细地捕捉到压力分布趋势，所以文中选用第三截面的实验数据进行对比。

2 计算结果及讨论

2.1 风洞实验压力分布

细长圆锥前体在大攻角时测得的40个滚转角的压力系数存在一些内在的规律，以自由来流速度30m/s、攻角35°、第三截面测得的数据为例进行分析。图1(a)为压力系数随周向角的变化，图1(b)为Cp，min，p(左舷最小压力系数)、Cp，min，s(右舷最小压力系数)和Cp，min，ave(左右舷最小压力系数平均值)随滚转角的变化，其关系如下:

可以看出:随着滚转角从0°变到360°，Cp，min，p和 Cp，min，s沿相反的趋势进行变化，基本成镜像状态，而它们在每一个滚转角下的平均值接近于一个常数，近似为 －0.85。通过对其它风速和攻角下的实验数据进行分析，也发现了同样的规律:40个滚转角下的Cp，min，ave接近于一个常数或在一个很小的范围内波动，但不同工况下得到的平均值的量值有所差异。

图1 30m/s，α =35°，第三截面

2.2 网格无关性与相关性

前期数值计算发现，轴向和径向网格数的密度不同对计算结果的影响很小，但圆锥前体的压力分布复杂，周向网格须分布的足够密。文中对四种不同的周向网格进行了计算，轴向和径向网格数均为81和121，周向网格数分别为120、180、240和360，称之为网格1到网格4。当计算步数为500步时，计算残差可降低2至3个量级。采用 Gamma theta转捩模型对网格无关性进行了分析。使用四套网格模拟了30m/s，α=25°和α=35°的流场。四套网格在攻角25°时的收敛解几乎重合，也就是说25°时的数值计算存在网格无关性，如图2(a)，不同的网格计算结果与实验数据吻合得都很好。而四套网格在35°时的收敛解相差较大，也就是说35°时不同的计算网格会得到不同的解，其原因可能是使用不同网格计算时产生了不同的截断误差，如图2(b)，计算结果与实验数据有一定偏差，但4个计算解的左右舷吸力峰处的压力系数平均值与实验数据一致，在 －0.85左右(参考2.1)。因此，采用四套网格得到的计算结果均视为合理的。

综合考虑计算残差和网格无关性等因素，文中采用网格3(81×121×240)进行后续计算。

图2 30m/s时，对不同周向网格第三截面的压力预测

2.3 不同模型计算结果比较

2.3.1 压力分布比较

对每一工况下40组不同滚转角的实验数据，从中分别挑选出与相同条件下第三截面计算结果最为吻合的一组来对比。此处只给出来流速度为30m/s的压力分布对比图作为典型进行分析，如图3。分析得到，层流模型对自由来流速度大于50m/s和攻角大于25°的流场的计算结果与实验差异较大;SST湍流模型对自由来流速度小于50m/s和攻角小于25°的流场的计算结果与实验差异较大;但对文中所计算的各风速和攻角下的流场，Gamma theta转捩模型的计算结果与实验数据均吻合较好。这说明，文中各工况下的流动主要处于转捩流动状态，适合采用Gamma theta转捩模型来模拟。

图3 30m/s，计算与实验压力分布对比，第三截面

2.3.2 云图分析

文中分析了不同风速和攻角下，圆锥前体第三截面上的总压系数云图、速度矢量和轴向涡量云图。这里只给出了来流速度为30m/s，攻角25°和35°时的云图(图4、图5)，分别用Cp0和ωxd/U∞来代表总压系数和涡量，计算模型为Gamma theta转捩模型。为了更好的显示主涡处的涡量，对云图中涡量的范围进行了人工选择，实际涡量范围要大一些。分析得出:每一个主涡在涡核处均存在一个总压系数最小值，但这并不是整个横截面上的最小总压系数，横截面最小总压系数处于紧贴物面的一层很薄的边界层内。同时，每一个主涡在涡核处还存在一个轴向涡量的最大绝对值，该最大值同样不是整个横截面上的最大值，横截面上的最大轴向涡量绝对值也存在于紧贴物面的一层很薄的边界层内。

进一步分析发现，左右舷主涡涡核处的最小总压系数随着攻角的增大而减小，随着自由来流速度的增大而增大，其纵坐标则随着攻角的增大而增大;涡核处的轴向涡量绝对值的最大值一般会随着攻角和自由来流速度的增大而减小，其位置与最小总压系数的位置非常接近。

图4 30 m/s，α =25°，Gamma theta 转捩模型，第三截面

图5 30 m/s，α =35°，Gamma theta 转捩模型，第三截面

3 结论

通过对圆锥前体上的绕流进行数值模拟，并与实验数据比较，可得出以下结论:

1)文中涉及的各种流场主要处于转捩区域，适合于用Gamma theta转捩模型进行模拟。

2)攻角小于25°的流动主要是对称分离流，计算时存在网格无关性。攻角大于25°的流动是非对称的，采用不同的网格进行计算会得到不同的数值解，其原因可能是采用不同网格计算时产生的截断误差不同。在数值计算中采用不同网格得到的这种效果和实验中转动滚转角得到的效果类似。

3)对于非对称流动，不同滚转角下的实验得到的左右舷最小压力系数的平均值接近于一个常数，而采用不同网格进行计算时得到的左右舷最小压力系数的平均值也接近于一个常数，这两个值大小相近。因此，可以说文中在大攻角下的数值模拟结果与实验数据是相符的。

4)攻角大于10°时，在左右舷处出现主涡，每个主涡在涡核处都有一个总压系数最小值和一个无量纲的轴向涡量绝对值的最大值。该最小总压系数随着攻角的增大而减小，随着自由来流速度的增大而增大，其纵坐标随着攻角的增大而增大。该轴向涡量最大绝对值一般会随着攻角和自由来流速度的增大而减小，其位置与最小总压系数的位置相近。

[1]Ericsson L E. Sources of high alpha vortex asymmetry at zero sideslip[J]. Journal of Aircraft，1992，29(6):1086－1090.

[2]Lowson M V，Ponton A J C.Symmetry breaking in vortex flows on conical bodies[J]. AIAA Journal，1992，30(6):1576－1583.

[3]Keener E，Chapman G，Cohen L，et al. Side forces on forebodies at high angles of attack and mach numbers from 0.1 to 0.7:Two tangent ogives，paraboloid and cone，NASA TM，X －3438[R].1977.

[4]Zilliac G G，Degani D，Tobak M.Asymmetric vortices on a slender body of revolution[J]. AIAA Journal，1991，29(5):667－675.

[5]Yang Y，Zheng B R，Yang W，et al.Numerical computation of pressure distributions over conical forebody at high angles of attack，AIAA 2012－1237[R].2012.

[6]Meng X S，Jia C，Qiao Z D，et al. Development of flow asymmetry over a 20-Degree circular cone at low speed，AIAA，2007－4118[R].2007.