尖拱细长旋成体低速大攻角绕流数值模拟

张浩南，李国辉，徐 宇

(空军航空大学 航空作战勤务学院，长春 130000)

细长体绕流在大攻角下会出现一种现象，即旋涡对称流动发展为非对称流动。当旋涡表现为非对称流动时，细长体会产生一个侧向力，进而产生偏航力矩。Zilliac等[1]进行的实验研究显示了四种不同的流动模式：α<30°时的稳定对称流动、30°<α<50°时的不稳定对称流动、50°<α<65°时的双稳态模式以及α>65°时的类卡门涡街(或随机尾迹)流态。这种流动状态转换的另一种表现为侧向力的转换[1-3]，了解这种现象的原理将对于航空航天应用具有重要价值,例如对于导弹和高机动性飞机而言，可以通过流动控制进而达到减阻的目的。大攻角下非对称流发展的研究是一项极具挑战性的工作，由于缺乏对流动机理的认识，在实验和数值模拟时常常会遇到很多困难，包括模型缺陷、来流条件的选择、和测量技术的局限性等等。为研究这种现象，众多学者对影响流动非对称性的因素进行了研究分析，例如马赫数[4]、雷诺数[5]、钝度[6]和细长比[7]等等。由于在传统风洞实验研究中，实验数据重现性很差,随着计算机技术的发展，数值模拟方法开始逐渐显示出它的优势。大多数数值模拟使用RANS方法进行研究，其中引入了人工扰动[8]或几何缺陷[9]是产生流动非对称所必需的；同样，部分学者使用LES方法研究[10]也必须添加人工扰动，而也有学者则认为，在没有任何扰动以及缺陷的情况下也可以得到非对称的旋涡[11]。

本文的目的是使用RANS方法研究细长体在大攻角下尾涡的发展，而既不引入几何缺陷，也不引入人为扰动。在五个攻角α=30°，α=40°，α=50°，α=55°和α=60°上进行研究，使用压力分布分析流场沿细长体的截面法向力和侧面力，此外，提出一种验证细长旋成体截面绕流拓扑结构的方法，结果证明，数值仿真结果与理论分析保持一致[12]。

本文的目的是使用RANS方法研究细长体在大攻角下尾涡的发展，既不引入几何缺陷，也不引入人工扰动。在5个攻角α=30°，α=40°，α=50°，α=55°和α=60°下进行研究，使用压力分布分析流场沿细长体的截面法向力和侧面力，此外，提出一种验证细长旋成体截面绕流拓扑结构的方法。

1 数值模拟方法

控制方程为三维粘性N-S方程，其守恒通量形式为

(1)

式中：U为非定常项；E、F、G是为无粘项；Eμ、Fμ、Gμ为粘性项，本文对N-S方程进行时间平均进而模拟湍流流动。通过Reynolds平均法对流动控制方程式(1)进行时均处理得到的雷诺平均N-S方程(RANS)为封闭RANS方程，选用SSTk-ω湍流模型，该模型对网格适应性较好[13]，其湍流动能k和比耗散率ω的输运方程为

(2)

Gω-Yω+Dω+Sω

(3)

式中：Gk表示湍流的动能；Gω为ω的生成项；Γk和Γω分别代表k与ω的有效扩散项；Yk和Yω分别代表k与ω的发散项；Dω代表正交发散项，有效扩散方程为

(4)

(5)

式中：σk和σω是方程的湍流能量普朗特数；μt为湍流黏度，计算公式如下：

(6)

其中S为旋率，系数α*使得湍流粘度产生低雷诺数修正，a1=0.31。该模型在近壁区等价于k-ω模型，在远离壁面的区域自动转换为标准的k-ε模型。对控制方程采用有限体积法进行求解，对于控制方程的空间离散方式，扩散项采用二阶中心差分格式，对流项采用三阶MUSCL格式，速度与压力的耦合方式选用PISO算法。

2 网格及边界条件

研究中使用的几何模型由细长比为3的尖拱状前身和细长比为10的圆柱状后身组成。模型的总长度为26 cm，直径D=2 cm。分别在尖拱状前身的11个截面(OG段)和圆柱状后身的9个截面(CY段)收集数据以便于检查分析[如图1(a)所示]，截面位置与文献[14]所提供的位置相匹配，尖拱段部分截面用以监视初始旋涡的发展。在选择计算域时，细长体前端延长至6D处，后端与径向均延长至15D处[图1(b)]。整体网格选用结构网格，第一层网格高度为 0.000 1D，对模型头部网格进行了细化，总网格数为150万(图2)。计算采用Fluent软件，边界条件为速度入口与压力出口，壁面选用无滑移绝热壁面，此外，自由来流速度固定为0.2马赫，湍流度固定为0.2%。

图1 截面选取以及计算域示意图

3 模拟结果与分析

在本节中，通过绘制x涡量等值线图，定性地评价了计算结果；通过截面侧向力和法向力系数定量地分析了计算结果；并且根据截面上奇点的数量与位置，对细长体的截面拓扑结构进行了验证。

图3显示了所有计算攻角下的x涡量等值线。可以看出：非对称性沿轴向发展，并随着攻角的增加而增加，在30°≤α<50°时旋涡基本保持对称涡流型，非对称性的发展在α=50°时开始，并且起初较为微弱，在尖拱段时相对较低，到圆柱段时加强。在α=60°时旋涡已发展为高度非对称的，非对称性在尖拱段末端已经开始显现。

图3 平均x涡量等值线

侧向力和法向力系数(分别为Cy和Cz)相对于攻角的变化如图4所示。可以看出，Cz随着攻角的增加而增加，这是由于随着攻角的增加，背风侧的吸力会越来越大，因此Cz会有所增加。而且，随着攻角的增加，Cz最大值的位置沿轴向向下游移动。在α=30°时，Cz最大值出现在OG(6)截面；在α=40°到50°时，Cz最大值的位置移到OG(7)截面；在α=60°时，Cz最大值出现在OG(8)截面，并且这也是所有计算攻角中的最大值。在所有计算攻角下，当达Cz到最大幅值后均沿轴向减小。另一方面，在30°≤α≤40°时，整个模型Cy值保持在较低水平，这表明尾流旋涡仍然是对称的。但是，Cy在α=50°时沿轴向显着增加，在OG(11)截面达到最大值，然后，侧向力沿轴向在负向和正向之间波动。当攻角增加到α=60°时，Cy达到所有攻角下的最大幅值。与Cz相反，Cy最大值的位置随着攻角增加沿轴向向前移动，在α=60°时Cy最大值与Cz最大值发生在相同位置，即OG(8)截面。α=55°和α=60°的曲线趋势类似于α=50°，即Cy方向不断在正向与负向之间改变，这表明尾流旋涡具有强烈的非对称性，非对称性随着攻角的增加而增强，并沿轴向从圆柱段移动到尖拱段。

图4 侧向力和法向力系数曲线

图5为攻角α=30°时OG(3)截面、OG(6)截面以及CY(1)截面的速度矢量分布示意图。从图中可以看出，在OG(3)截面,流动还未产生分离,截面流型为附着结构(图5(a))；沿轴向往下游发展，这种附着结构开始发生变化，细长旋成体背涡出现了流动分离现象，到OG(6)截面，在细长旋成体左右两侧产生一对对称旋涡，由于此时对称旋涡刚开始出现，流动的分离还较为微弱，卷入旋涡中的涡量较少，使得产生的旋涡为体积小且强度弱的一对对称旋涡，故而在卷成旋涡之后很快就在背风侧的对称面附近发生再附(图5(b))，当旋涡处于对称流态时，两个头涡之间的诱导作用是相互平衡的，并且他们之间相距较远，相互挤压的作用微乎其微，这就保证了这种对称涡流态是一种稳定的平衡状态；随着流动沿轴向继续往下游发展，分离现象进一步增强，分离剪切层不断向两个头涡输入涡量，使得旋涡强度增强体积增大，并且相互靠近进而伴随着挤压现象的发生(图5(c))，此时旋涡仍保持对称涡流态，两个头涡之间的诱导作用依旧是相互平衡的，但是由于他们不断靠近，已经产生了较为强烈的挤压作用，使得此时的对称涡流态成为一种不稳定的平衡状态。

图5 速度矢量分布示意图

为探究细长旋成体对称涡流态的演变，基于Lowson和Ponton[15]给出的细长旋成体绕流截面流型拓扑结构(图6)，为获得绕流截面流型拓扑结构中鞍点的精确位置，图7给出了α=30°攻角下OG(3)截面、OG(6)截面以及CY(1)截面上壁面切向速度随方位角θ的变化曲线，进而确定近壁面奇点位置，除此之外，通过寻求速度矢量与截面法向量夹角的0点来确定空间鞍点的精确位置，图8(a)与图8(b)分别为该攻角下空间奇点沿轴向不同截面纵向与横向位置变化曲线，由于细长旋成体头部为尖拱型，不同截面对应的模型半径不同，为准确表示空间奇点与壁面的距离，在图8(a)中选用空间奇点纵向坐标与该截面模型半径的差值作为y轴坐标值，从图8可以看出，随着截面沿轴向往下游发展，空间奇点逐渐远离细长体壁面，线性程度较高，且一直处于背风侧对称面附近很小的范围之内。通过上述分析可得，对于OG(3)截面，切向速度存在两个0点，分别存在于θ=0°和θ=180°处，而由于壁面是无法穿过的，所以细长旋成体壁面的径向速度始终为0，进而可以得出，在θ=0°和θ=180°处存在两个奇点，并且通过壁面切向速度的流动方向可以判断θ=0°处为再附型半鞍点Rh1，θ=180°处为分离型半鞍点Sh1，空间中无奇点；同理，对于OG(6)截面，除去θ=0°和θ=180°处的两个奇点，在背风侧对称面两侧附近存在两个位置对称的再附型半鞍点Rh2和Rh3(θ=140°与θ=220°)，而在距离背风侧对称面两侧更远一些的位置存在两个位置对称的分离型半鞍Sh2和Sh3(θ=90°与θ=270°)，空间中无奇点，在整个细长旋成体壁面上，对称面两侧附近的两个再附型半鞍点在旋成体壁面上代表两条再附线；对于CY(1)截面，从图8可以看出，相比于OG(6)截面，θ=180°处的奇点由分离型半鞍点转变为再附型半鞍点Rh2，并且细长体背风侧对称面左右两侧的一对再附型半鞍点也不复存在，这代表细长旋成体壁面上将不再有再附线，此时空间中存在一个奇点S1。不难得出，OG(3)截面、OG(6)截面以及CY(1)截面的对称涡流型分别对应图6所示的3种流动拓扑机构。

图6 细长体绕流拓扑结构示意图

图7 切向速度曲线

图8 空间鞍点的纵向位置与横向位置曲线

在上述分析中，验证了两涡对称流型的两种拓扑结构，它们的主要区别在于两个再附型半鞍点的消失以及θ=180°处奇点类型的转变，为此对这3个奇点的演变进行进一步研究，由于对称涡流型的转变是一种渐变的过程，对OG(6)截面与CY(1)截面中间数个截面的壁面切向速度进行了分析，图9为其中部分截面上3个奇点的位置示意图，从图中可以看出细长体背风侧对称面左右两侧的一对再附型半鞍点Rh2和Rh3会不断的靠近θ=180°处的分离型半鞍点Sh1，直至3个奇点重合，转变为高阶奇点O，此时在θ=180°附近已再无任何奇点存在，称这个状态为临界流动状态，从上述分析可以得出，临界流动状态的高阶奇点O由图6(b)中Rh2、Rh3以及Sh1重合而成，此时壁面存在两个分离型半鞍Sh1和Sh2，一个再附型半鞍点Rh1和一个高阶奇点O，且空间中无奇点存在，拓扑结构如图10所示，与文献[12]提出的临界状态保持一致。

图9 鞍点位置示意图

图10 临界状态拓扑结构示意图

4 结论

1) 在无任何扰动与不规则性的前提下，可以得到非对称的旋涡，且非对称性随着攻角的增加而增强，并沿轴向逐渐增强。

2) 确定截面绕流中奇点的精确数量及位置对截面绕流拓扑结构的验证，结果与理论分析一致。

3) 在数值模拟中通过截面绕流中特殊奇点的位置变化可观察到对称涡流型的临界流动状态。