计算流体力学(computational fluid dynamics, CFD)数值计算相对于经验公式计算和风洞实验具有特有的优势。CFD数值模拟可以很好的预测复杂外形模型的流场结构^[1-4]，能使更加深刻地理解问题产生的机理，为实验提供指导，模拟风洞实验无法模拟的条件，节省实验所需的人力、物力和时间，并对实验结果的整理和规律的得出起到很好的指导作用。Baum等^[5]分别对F-16战机副油箱下落过程进行风洞实验、飞行实验和CFD数值计算，得到副油箱下落过程随时间变化的曲线。从文献[5]中可以看出，CFD数值模拟结果相比风洞实验更接近飞行实验结果，说明CFD数值模拟具有很高模拟精度和工程应用价值。文献[6]分别采用风洞实验、飞行实验和计算软件3种方法预测弹箭气动参数，并作了误差统计。得到的结论是风洞实验无法准确预测俯仰阻尼力矩、马格努斯力矩、滚转阻尼力矩等系数。

CFD模拟结果的好坏很大程度上取决于湍流模型的准确度。Spalart等^[7-8]在1997年提出的DES (detached eddy simulation)方法是一种(reynolds average navier-stokes / large eddy simulation)RANS / LES混合方法，DES整场还是求解湍流模型的输运方程，只是将实时变化的长度尺度作为转换器，在边界层内流动仍由RANS模化，并未显著增大计算量；在远离壁面区域达到类似LES的滤波效果。DES模型的提出一开始是为了改进RANS对非定常大尺度分离湍流的模拟，但DES可能出现雷诺应力提供不足且边界层内模化应力不足缺陷, 导致了严重的速度型偏离对数律现象或者网格诱导分离现象^[9]。混合模型具有计算量比LES小且计算精度高于RANS的显著优势。Menter^[10-11]借鉴DES中长度尺度的思路，基于剪切应力传输SST湍流模型，并在SST模型中引入可根据当地流动拓扑自动调整的von Karman尺度作为湍流模型的第二长度尺度, 而不是使用LES的滤波尺度, 由此提出了尺度自适应模拟(scale adaptive simulation, SAS)的概念。它是一种新的RANS/LES混合模型，SAS模型长度尺度不依赖于计算单元的大小，可以依靠流动状态动态自适应调整湍流黏性，并克服了RANS/LES交界面问题，且SAS模型远场处耗散小，从而可以自适应的分辨湍流结构。从2003年至今，这种极具吸引力和潜力的尺度自适应模拟方法，越来越得到研究者的重视。

1 动力学模型

本文选用M910旋转稳定弹丸为物理模型，几何参数、飞行速度、旋转速度等参数详见文献[14]。

1.1 控制方程和湍流模型

直接求解高度非线性的三维瞬态纳维斯托克斯(N-S)方程，需要采用对计算机内存和速度要求很高的直接模拟方法，在实际工程中目前还不能采用此方法。工程上对瞬态的N-S方程做时间平均处理，得到雷诺平均的N-S方程(RANS):

$ \frac{{\partial \rho }}{{\partial t}} + \frac{\partial }{{\partial {x_i}}}(\rho {u_i}) = 0 $

(1)

$ \begin{array}{l} \frac{{\partial (\rho {u_i})}}{{\partial t}} + \frac{\partial }{{\partial {x_j}}}(\rho {u_i}{u_j}) = - \frac{{\partial p}}{{\partial {x_i}}} + \\ \frac{\partial }{{\partial {x_j}}}\left[{\mu \frac{{\partial {\mu _i}}}{{\partial {x_j}}}-\rho \overline {{{u'}_i}{{u'}_j}} } \right]+ {S_i} \end{array} $

(2)

式中：t为时间；下标i、j的取值范围是(1、2、3)，代表坐标系方向；ρ、p和μ分别是流体密度、压力和动力黏度；S_i是动量守恒方程的广义源项。定义-ρ $ \overline {{{u'}_i}{{u'}_j}} $为Reynolds应力，表示湍流的影响，其中符号“—”代表时间平均值，“′”代表脉动值。

RANS方程不封闭，因此需要引入湍流模型来封闭方程组。标准的SST湍流模型^[13]，具有在边界层内能很好模拟低雷诺数流动，在边界层外能很好模拟完全湍流流动的优点。k和ω的输运方程为

$ \frac{{{\rm{d}}\left( {\rho k} \right)}}{{{\rm{d}}t}} = {\tau _{ij}}\frac{{\partial {u_i}}}{{\partial {x_j}}} - {\beta ^*}\rho \omega k + \frac{\partial }{{\partial {x_j}}}\left[{(\mu + {\sigma _k}{\mu _t})\frac{{\partial k}}{{\partial {x_j}}}} \right] $

(3)

$ \begin{array}{l} \frac{{{\rm{d}}\left( {\rho \omega } \right)}}{{{\rm{d}}t}} = \frac{{\gamma \rho }}{{{\mu _t}}}{\tau _{ij}}\frac{{\partial {u_i}}}{{\partial {x_j}}} - \beta \rho {\omega ^2} + \\ \frac{\partial }{{\partial {x_j}}}\left[{(\mu + {\sigma _\omega }{\mu _t})\frac{{\partial \omega }}{{\partial {x_j}}}} \right] + 2\rho (1 - {F_1}){\sigma _{{\omega ^2}}}\frac{1}{\omega }\frac{{\partial k}}{{\partial {x_j}}}\frac{{\partial \omega }}{{\partial {x_j}}} \end{array} $

(4)

具体湍流模型动力学方程、模型中u_r为湍流黏度，它是空间坐标的函数，取决于流动状态，而不是物性参数，下标“t”表示湍流，常数β^*、β、σ_κ、γ及函数F₁的具体形式详见文献[13]。

传统的湍流模型产生项除了包含湍流的某一个尺度外并没有其他尺度，而对于湍流长度尺度的计算最少需要2个湍流尺度，因此不具备根据湍流长度尺度来进行调节生成和毁灭的能力。在Menter的改进模型中，产生项包含了von Karman长度尺度，该尺度在流动未失稳时远远大于流动失稳，具备根据当地流场结构调节当地耗散的能力。同时von Karman长度尺度仅仅只与当地的流场速度变量相关，也就使得湍流尺度的自适应性极大提高，同时von Karman长度尺度引入了速度导数，而速度变量属于输运量，因此该尺度是可以反映一定的湍流输运效应，从而改善湍流模型的预测精度。Menter将这种尺度自适应特性引入了SST模型，并形成了SST-SAS湍流模型。SST-SAS湍流模型实际上只在标准SST湍流模型的基础上修改其ω输运方程，在ω方程中增添一个自适应源项Q_SAS^[12]

$ \begin{array}{l} {Q_{{\rm{SAS}}}} = {\rm{max}}\left[{\rho {\zeta _2}\kappa {S^2}\left( {\frac{L}{{{L_{vk}}}}} \right){^2}-} \right.\\ \left. {C\frac{{2\rho k}}{{{\sigma _\phi }}}{\rm{max}}\left( {\frac{{|\nabla \omega {|^2}}}{{{\omega ^2}}}, {\rm{ }}\frac{{{\rm{|}}\nabla k{|^2}}}{{{k^2}}}} \right), 0} \right] \end{array} $

(5)

$ {L_{vk}} = \frac{{\kappa \sqrt {2{S_{ij}}{S_{ij}}} }}{{\sqrt {{{({\nabla ^2}u)}^2} + {{({\nabla ^2}v)}^2} + {{({\nabla ^2}w)}^2}} }} $

(6)

$ L = \sqrt k /(c_\mu ^{1/4}\omega ) $

(7)

式中：${L_{vk}} = \frac{{\kappa \sqrt {2{S_{ij}}{S_{ij}}} }}{{\sqrt {{{({\nabla ^2}u)}^2} + {{({\nabla ^2}v)}^2} + {{({\nabla ^2}w)}^2}} }} $和$L = \sqrt k /(c_\mu ^{1/4}\omega ) $分别为von Karman尺度和模化湍流力尺度，具体动力学模型方程、模型常数ζ₂、κ、C、$ {\sigma _\phi }$、C_μ详见文献[10-12]。

1.2 计算网格以及边界条件

网格的合理设计和高质量生成是CFD计算的前提条件^[16]。对于亚音速到超音速粘性流动，物体近壁面处的流场会产生急剧的变化，因此在近壁面处和激波处需要加密网格。同时，为了准确预测滚转阻尼，以及马格努斯力矩等气动参数，在边界层内对网格进行加密，保证第一层网格高度的无量纲量y⁺≤1。采用多Block方法生成流场拓扑结构，然后生成了合理的且高质量的流场结构化计算网格，网格数量为116万。图 1(a)、(b)、(c)分别为流场拓扑结构、流场计算网格及流场网格剖面图。远场边界采用速度进口和压力出口，在物面上采用无滑移条件来处理，对计算模型采用moving-wall设置壁面旋转。

1.3 数值格式

采用基于有限元的有限体积法，在控制体内引入插值点，实现双重网格化，对六面体网格单元采用24点插值，而单纯的有限体积法仅仅采用6点插值，在保证了有限体积法的守恒特性的基础上, 吸收了有限元法的数值精确性。采用全隐式完全耦合求解技术，克服了传统算法需要“假设压力项-求解-修正压力项”的反复迭代过程, 而同时求解动量方程和连续性方程, 再加上其采用多网格技术，计算速度和稳定性相比传统方法提高1~2个数量级^[19]。

2 气动计算公式及批文件处理方法 2.1 计算公式

CFD求解气动特性系数的原理是，先对流场进行求解，得到流场每一个网格上的速度、密度、压力等参数，然后根据计算公式求出气动系数:

$ \left\{ \begin{array}{l} {C_D} = [{F_X}{\rm{cos}}\alpha + {F_Y}{\rm{sin}}\alpha]/(0.5\rho {v^2}S)\\ {C_L} = [{F_Y}{\rm{cos}}\alpha-{F_X}{\rm{sin}}\alpha]/(0.5\rho {v^2}S)\\ {C_N} = {F_Y}/(0.5\rho {v^2}S)\\ {C_m} = {M_Z}/(0.5\rho {v^2}Sl)\\ {C_l} = {M_X}/(0.5\rho {v^2}Sl)\\ {C_n} = {M_Y}/(0.5\rho {v^2}Sl)\\ {x_{cp}} = {x_{cg}} - ({C_m}/{C_N}) \end{array} \right. $

(8)

式中：ρ为密度，v为来流速度，S为截面积，l为参考长度，此处取弹径为参考长度。由CFD求解得到弹体轴向力F_X和法向力F_Y。M_X, M_Y, M_Z分别是弹体滚转力矩、偏航力矩、俯仰力矩。x_cg为重心位置系数。通过式(8)可以求出阻力系数C_D、法向力系数C_N、俯仰力矩系数C_m、滚转力矩系数C_l、偏航力矩系数C_n、压心位置系数x_cp。

对于外形简单，小攻角飞行的旋转弹丸，可采用CFD与经验公式(9)结合的方法计算其气动导数^[17-18]

$ \left\{ \begin{array}{l} {C_{N\alpha }} = {C_N}/{\rm{sin}}\alpha \\ {C_{m\alpha }} = {C_m}/{\rm{sin}}\alpha \\ {C_{lp}} = - {C_l}/{{\bar \omega }_x}\\ {C_{np\alpha }} = - {C_n}/[{{\bar \omega }_x}{\rm{sin}}\alpha] \end{array} \right. $

(9)

式中: $ {{\bar \omega }_x} = ({\omega _x}l)/\left( {2v} \right)$，ω_x为弹丸绕弹轴的定角速度。

2.2 批文件计算方法

为了节省人为操作的大量时间，无需每次读入网格、设置边界条件、保存计算结果等人为操作。本文采用CFX软件的二次开发功能编制了CCL文件和批文件计算的脚本，计算过程中不需要暂停，并且可以基于上一次计算结果作为初始流场，进行下一个case的计算，提高了批文件计算的效率。批文件计算的过程如下：

1) 首先采用CEL语言编写包含来流参数及气动系数的一个CCL文件：

LIBRARY：

CEL：

  EXPRESSIONS：

   Uinf=136.1[m/s]

   ……

  cnpa=cnp/sin(AOA)

   cYpa=cY/(wx*sin(AOA))

  END

 END

END

  COMMAND FILE:

    Version=16.0.0

END

2) 将改CCL文件、网格读入CFX软件，定义边界条件，选择计算方法，然后生成包含所有计算信息的case1.def文件。修改该CCL文件中的来流速度、攻角、转速等信息然后另存为case2.def、case3.def、case4.def.....

3) 编写脚本语言并将该文档后缀名改为.bat。脚本语言为

PATH C:\\Program Files\\ANSYS Inc\\v160\\CFX\\bin start/wait cfx5solve -case1.def -ini-file initial1.resdouble -par-local -partition 16 -fullname case1result start/wait cfx5solve -case2.def -ini-file initial1.resdouble -par-local -partition 16 -fullname case2result ……

4) 将所有.def文件、一个.dat文件放到一个文件夹下，双击.bat文件便可以计算所有case且保存了所有case的计算结果，从而提高了仿真工作效率。

3 旋转弹丸气动特性仿真结果及分析

本文首先采用基于SST湍流模型的定常CFD对不同马赫数Ma下的M910旋转稳定弹丸进行数值计算，并与实验数据以及文献[14]数值计算结果对比。从图 2可以看出基于SST湍流模型的定常CFD数值模拟能力要比PRODAS软件模拟结果更好, 模拟精度与文献[14]采用的Realizable k-ε湍流模型定常CFD模拟精度相当。与实验数据相比，阻力系数、法向力系数导数、俯仰力矩系数导数以及压心系数误差均在10%以内，模拟结果较好。滚转阻尼系数在超音速部分计算误差较小，在亚音速和跨音速部分计算误差较大。马格努斯力矩系数导数在超声速部分计算误差较小，在亚音速和跨音速部分模拟误差很大。说明基于两方程涡粘模型的定常CFD模拟不能满足滚转阻尼系数和马格努斯力矩系数导数在亚音速和跨音速部分的模拟要求。

在SST湍流模型的模拟结果基础上，将定常模拟结果作为初始流场，采用SST-SAS湍流模型对对M910旋转稳定弹丸进行尺度自适应模拟，这样可以大大节省计算时间。从图 2可以看出，采用基于SST-SAS湍流的非定常CFD模拟可以很好的计算从亚音速到超音速全段的阻力系数、法向力系数导数、俯仰力矩系数导数、压心系数、滚转阻尼系数以及马格努斯力矩系数导数。相比文献[14]采用的基于DES湍流模型的非定常CFD模拟，采用SST-SAS计算旋转弹丸在亚音速和跨音速段的滚转阻尼系数和马格努斯力矩系数导数有更好的模拟精度。

本文采用惠普16核32G内存工作站进行计算，数值模拟基于ANSYS CFX商业软件及CEL语言二次开发，编写了CCL文件以及用于批文件计算的脚本。采用基于SST湍流模型的定常CFD模拟，亚音速段单个文件只需要5~10 min，气动系数便可以收敛。如图 3(a)所示，计算Ma=0.5, α=3°, ω=1 789 rad/s时M910旋转弹丸的气动参数，只需要累计70个迭代步就可以收敛，计算时间为8 min，相同计算平台情况下，相同的计算网格计算该例子，FLUENT软件(压力基求解器，即分离求解)需要几千步的迭代，一个多小时的计算才能收敛，采用FLUENT软件(密度基求解器，即耦合求解)也需要上千步的迭代，半小时的计算。这说明本文采用的多网格完全耦合计算方法的计算速度快。在跨音速和超音速段计算单个文件需要15~30 min，气动系数便可以收敛。如图 3(b)所示，计算Ma=3.5, α=3°, ω=12 522 rad/s时M910旋转弹丸的气动参数，累计300个迭代步气动系数收敛，计算时间为25 min。采用本文数值方法并同时采用批文件计算方法，计算完13组工况，只需要55 min。在SST湍流模型计算结果的基础上，进行M910旋转弹丸尺度自适应模拟，时间步长为1E-6 s，计算一组参数需要大约2~12 h，其中超音速段计算需要时间少，亚音速和跨音速段需要计算时间较长。图 4为部分马格努斯力矩系数导数随仿真时间的变化曲线。

为了精确模拟滚转阻尼系数以及马格努斯效应，必须确保y⁺≤1。如图 5所示，给出了从亚音速到超音速沿着弹体长度分布的y⁺值分布。可以看出，在亚音速到超音速全速度段范围内，y⁺≤0.6。

图 6(a)、(b)、(c)分别为基于非定常SST-SAS湍流模型，Ma=2.5, α=3°, ω=8 944 rad/s时数值计算得到的马赫数云图、压力云图、纵向截面速度矢量云图，从图 6(a)可以看出弹头部和尾部形成了清晰的斜激波，从图 6(b)可以看到在弹肩和弹底部区域形成了低压区。由于存在攻角，M910弹丸周围的马赫数云图和压力云图在弹轴两侧不对称。弹体下方压力值和马赫数略大于弹体上方。仿真结果符合空气动力学规律。

图 6(c)中也可以看到，贴近弹体壁面处的速度矢量是沿着弹体旋转的方向，由于空气粘性的作用，然后渐变为轴向方向，形成沿弹丸表面的外法线方向的速度梯度。

图 7为Ma=2.5, α=3°, ω=8 944 rad/s时基于非定常SST-SAS模型，模拟得到重心横截面处的速度矢量图和局部处放大图。从图 7(a)中可以看出，由于弹体旋转造成重心横截面处的速度矢量图分布不对称。图 7(c)和图 7(e)分别为弹丸右侧背风面Ⅱ处和迎风面Ⅳ处局部速度矢量图。弹丸右侧表面切向速度方向与横流方向相同。图 7(b)和图 7(d)分别为弹丸左侧背风面Ⅰ处和迎风面Ⅲ处的局部速度矢量放大图。从图中可以看到在边界层内弹体表面切向速度方向与横流方向相反，发生了畸变。由于速度矢量不对称、边界层畸变导致由速度梯度产生的周向切应力分布不对称，而周向切应力不对称就会导致弹体产生一个侧向力，这也是马格努斯力的成因。

图 8~10分别给出了亚音速、跨音速、超音速两种湍流模型模拟结果的流场结构。图 8、9为亚音速模拟，从图中可以看出采用非定常SST-SAS湍流模型模拟得到的流场包含丰富的流动结构, 而采用定常SST模拟则基本将小尺度的结构抹平了, 从滚转阻尼系数和马格努斯力矩系数导数计算结果也可以看出，非定常尾流是否能准确模拟对计算亚音速和跨音速段的滚转阻尼和马格努效应影响很大。图 10为超音速段的数值模拟结果，从图中可以看出，采用定常SST模型和非定常SST-SAS模型的模拟结果几乎一致，说明在超音速段两种模型对可以较好的计算滚转阻尼和马格努斯力矩系数导数。

4 结论

1) 采用全隐式多网格完全耦合求解技术，同时求解动量方程和连续性方程，在计算速度上相比传统的“假设压力项-求解-修正压力项”的反复迭代算法速度要快许多。

2) 基于SST湍流模型的定常CFD数值模拟方法可以较好的计算阻力、法向力、俯仰力矩和压心系数，不能准确预测亚音速和跨音速阶段旋转弹丸的滚转阻尼系数以及马格努斯力矩系数导数。

3) 基于SST-SAS湍流模型的非定常CFD数值模拟方法可以较好的计算旋转弹丸亚音速到超音速全段的滚转阻尼系数和马格努斯力矩系数导数。同时，计算结果表明，能否准确捕捉非定常尾涡效应对滚转阻尼及马格努斯力矩计算影响很大。