0 引言

高速旋转会产生很强的壁面剪切层^[1-2], 湍流剪切应力对流线弯曲和模型旋转很敏感, 在凹角区域及旋转管道的抽吸侧, 湍流都得到了增强, 在凸表面及旋转管道的受压侧, 湍流都得到了削弱, 传统的基于Boussinesq线性涡黏性假设的湍流模型(如一方程Spalart-Allmaras(SA)模型、两方程k-ω Shear-Srtess-Transport(SST)模型、k-epsilon(k-ε)模型)不能准确的捕捉到这种特性^[3-4]，因为当模型旋转时，科氏力(Coriolis forces正比于二阶应力分量)以旋转生成项的形式在雷诺应力张量中体现，但是输运方程中应力张量的对角项是相互抵消的，旋转效应体现不明显。雷诺应力方程中可以显含弯曲和旋转效应的生成项，但是其昂贵的计算代价，使得在实际应用中还鲜有见到。对涡黏性模型进行弯曲和旋转效应修正，使其实现或近于实现：1)单纯的旋转运动不引入额外的湍流；2)近壁面的流线弯曲在湍流特性中有体现；3)不增加太多的计算计时；4)能应用到工程问题中，是近半个世纪高速旋转问题研究中众多学者不懈奋斗的方向。

早期的弯曲和旋转效应修正中存在着几种思想：一种是修正湍流长度尺度^[5-6]，是最直接和直观的方法，因为涡黏性μ_T正比湍流长度尺度，修改后将直接影响湍流黏性的计算，如修正湍流长度尺度为$l = \frac{{{l_0}}}{{l + {C_{rc}}{R_i}}}$或$l = {l_0}\left( {1 - {C_{rc}}{R_i}} \right)$；第二种思想是修正输运方程ω方程或ε方程中的衰减项^[7-8]，如修正ω方程中衰减项为${D_\omega } = - {F_4}\beta \rho {\omega ^2}, {F_4}{\rm{ = }}\frac{1}{{l + {C_{rc}}{R_i}}}$，通过调节Ri(Richardson Number)可实现不同的效果；第三中思想是修正输运方程k方程中的生成项，比如Wilcox^[9]在k方程中加入源项9μ_T(∂u/∂y)实现修正。上述修正都存在一个问题，就是只限定于某类流动或某类坐标系统，并不能适用于任意的三维流动。

1997年Spalart和Shur^[10-11]提出了一种半经验的通用三维修正模型(SARC，RC代表弯曲和旋转修正)，在生成项前引入fr1函数，函数中引入应变率张量拉格朗日导数(Lagrangian Derivative)对湍流中的弯曲和旋转量进行简单模化，为保证伽利略不变性(Galiean-invariant)，求导沿着主应变轴进行。随后对雷诺数Re=11500~36000、滚转数Ro=0.0~0.21管道旋转流动进行了计算，与实验和直接数值模拟(DNS)结果比较后发现壁面速度型、分离区、凹凸区域的摩擦系数都比SA和SST模型有了很大改善，但存在在小滚转数下摩擦速度偏大，大雷诺数下弯曲特性预估不足、再附区偏长等问题。1998年Hellsten^[12]通过对三维流动中弯曲和旋转效应敏感的物理量进行分析，在k-ω SST模型中引入弯曲函数F4，同时修正壁面边界的ω值，提出了一种通用修正模型(SSTRC)，并对马赫数Ma=0.2，Re=35000旋转管道流进行了计算，发现Ro小于0.21时与实验符合很好，大滚转数时失效，随后对ONERA-A翼形计算发现与原始的k-ω SST结果相差不大。2004年Mani^[13-14]使用Wind代码，对Ma=0.2，Re=3.28×10⁶亚声速U型管算例对比了SARC和SSTRC修正模型，发现在分离区前都与实验结果符合很好，但是分离区出现压力过估、再附区出现摩阻预估不足的问题。2012年Sunil^[15]使用FOAM代码，通过模化湍流输运方程中的湮灭项实现对旋转和流线弯曲的控制，并对旋转的槽道、凹曲面的边界层进行了数值模拟，发现比SSTRC模型预测结果略优。2013年Nash^[16]使用FUN3D代码，采用SA和SARC模型，对Ma=0.15、Re=4.6×10⁶翼尖涡脱落及向后传播的历程进行了数值模拟，发现SARC模拟的涡核轴向速度和压力与实验结果符合更好。2015年Rumsey^[17]使用CFL3D和FUN3D代码，对Ma=0.2，Re=4.2×10⁶、迎角α=6°~37°的NASA梯形翼进行了计算，发现SARC、SSTRC模型预测的升力更大、对翼尖涡及向下游传播的历程刻画更精细，但没有改善翼尖附近的压力分布。

国内对湍流模型中弯曲和旋转效应修正的研究还比较少见，只调研到1997年李卫东^[18]采用不同的紊流黏性系数构造了一种修正的k-ε模型，并计算了液-液旋流分离器内锥变截面流场。弯曲和旋转修正的湍流模型在弹箭旋转气动特性方面的研究，国内、外还未见到。为此，本文采用完全时间相关的非定常Navior-Stokes(N-S)方程，以一方程SA湍流模型和两方程SST湍流模型为基础，对典型带翼弹箭开展研究，从气动特性和流场结构分析了SARC及SSTRC修正模型对计算精度的影响。

1 数值方法 1.1 控制方程

积分形式的时间相关三维可压缩N-S方程为：

$ \frac{\partial }{{\partial t}}\int_\mathit{\Omega } \mathit{\boldsymbol{Q}} {\rm{d}}V + \int_{\partial \mathit{\Omega }} {\left( {\mathit{\boldsymbol{\hat F}} - \mathit{\boldsymbol{\hat G}}} \right)} \cdot \mathit{\boldsymbol{\hat n}}{\rm{d}}S = 0 $

(1)

Ω是任意形状的控制体，dS是控制体上的微元面积，${\mathit{\boldsymbol{\hat n}}}$是微元面的外法向单位矢量。Q为守恒变量，${\mathit{\boldsymbol{\hat F}}}$、${\mathit{\boldsymbol{\hat G}}}$分别为无黏通量、黏性通量。

本文计算采用的湍流模型包括SA一方程湍流模型：

$ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {\rho \tilde \upsilon } \right) + \frac{\partial }{{\partial {x_i}}}\left( {\rho \tilde \upsilon {u_i}} \right) = {f_{r1}}{C_{b1}}\rho \bar S\tilde \upsilon {\rm{ + }}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{{{\sigma _v}}}\left[ {\frac{\partial }{{\partial {x_i}}}\left( {\left( {\mu + \rho \tilde \upsilon } \right)\frac{{\partial \tilde \upsilon }}{{\partial {x_i}}}} \right)} \right.{\rm{ + }}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{C_{b2}}\rho \frac{{\partial \tilde \upsilon }}{{\partial {x_i}}}\frac{{\partial \tilde \upsilon }}{{\partial {x_i}}}} \right] - {C_{w1}}\rho {f_w}{\left( {\frac{{\tilde \upsilon }}{d}} \right)^2} \end{array} $

(2)

以及k-ω SST两方程湍流模型：

$ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {\rho k} \right) + \frac{\partial }{{\partial {x_i}}}\left( {\rho {u_i}k} \right) = \\ \;\;\;\;\;\frac{\partial }{{\partial {x_i}}}\left( {\left( {\mu + {\sigma _k}{\mu _t}} \right)\frac{{\partial k}}{{\partial {x_i}}}} \right) + {P_k} - {\beta _k}\rho \omega k \end{array} $

(3)

$ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {\rho \omega } \right) + \frac{\partial }{{\partial {x_i}}}\left( {\rho {u_i}\omega } \right) = \frac{\partial }{{\partial {x_i}}}\left( {\left( {\mu + {\sigma _\omega }{\mu _t}} \right)\frac{{\partial \omega }}{{\partial {x_i}}}} \right) + \\ \;\;\;\;\;\;\;{P_\omega } - {F_4}\rho {\beta _\omega }{\omega ^2} + \frac{{2\rho \left( {1 - {F_1}} \right){\sigma _{\omega 2}}}}{\omega }\frac{{\partial k}}{{\partial {x_i}}}\frac{{\partial \omega }}{{\partial {x_i}}} \end{array} $

(4)

第一种弯曲和旋转修正采用文献[10]中的Spalart发展的SARC方法：

$ \begin{array}{l} {f_{r1}}\left( {{r^*}, \bar r} \right) = \left( {1 + {c_{r1}}} \right)\frac{{2{r^*}}}{{1 + {r^*}}} \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left[ {1 - {c_{r3}}{{\tan }^{ - 1}}\left( {{c_{r2}}\bar r} \right)} \right] - {c_{r1}} \end{array} $

(5)

$ {r^*} = S/\omega , \bar r = 2{\omega _{ik}}{S_{jk}}\left( {\frac{{D{S_{ij}}}}{{{D_t}}}} \right)/{D^4} $

(6)

其中，f_r1为修正函数，通过调节生成项${C_{b1}}\rho \bar S\tilde \upsilon $来体现弯曲和旋转效应的影响；${S_{ij}} = 0.5\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right)$为对称张量，${\omega _{ij}} = 0.5\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right)$为反对称张量；D²=ω_ij²+S_ij²，系数c_r1=12，c_r2=1.0，c_r3=1.0。

第二种弯曲和旋转修正采用文献[12]中Hellsten发展的SSTRC方法：

$ {F_4} = 1/\left( {1 + {C_{RC}}{R_i}} \right), {R_i} = \frac{{\left| {{\omega _{ij}}} \right|}}{{\left| {{S_{ij}}} \right|}}\left( {\frac{{\left| {{\omega _{ij}}} \right|}}{{\left| {{S_{ij}}} \right|}} - 1} \right) $

(7)

其中，F₄为弯曲函数，作用于ω方程中的衰减项ρβ_ωω²来实现对弯曲和旋转效应的修正，R_i下界为-1/4，系数C_RC参考文献[12]中取为3.6。

1.2 数值求解方法

求解采用格心格式的结构有限体积法，Roe格式计算无黏通量，反距离权重最小二乘法计算单元梯度，Venkatakrishnan限制器抑制间断附近的过冲和振荡。基于Riemann条件的远场边界条件，Riemann不变量构造时需要计及网格速度；壁面边界条件要求壁面流体速度和网格运动速度相同。非定常计算采用双时间步法，物理时间层采用二阶向后差分离散，伪时间层采用Low Upper Symmetric Gauss-Seidel(LU-SGS)隐式时间推进。

1.3 计算模型及计算网格

计算模型是美国空军标模Air Force Modified Basic Finner Missile(AFF)^[19-20]，由弧形段、圆柱段及尾翼(舵)组成，典型的细长体正常式布局导弹，尾翼采用前缘后掠角为53.13°的梯形翼，绕弹身呈“十”字形分布，全弹长细比为10。计算状态如表 1，共4种组合，每组7个迎角，转速Ω^*=771.6 rad/s，无量纲旋转速度Ω=0.03，Ω定义为$\mathit{\Omega = }\frac{{{\mathit{\Omega }^*}D}}{{2{V_\infty }}}$, V_∞为来流速度，D为直径。一个计算周期采用1440个物理时间步，对应Δt=5.655×10^-6s。气动力和力矩参考点位于质心(5D，0，0)，参考长度为直径，参考面积为最大横截面积。

本文参考文献[20]的网格分布，生成的计算网格(图 1)流向×法向×周向为：250×200×250，网格单元数共1255万，法向第一层网格间距均为5×10^-6m，保证壁面y⁺≤1。

1.4 旋转周期和角度定义

图 2给出了旋转周期及角度定义，从弹体底部向前看，逆时针旋转为正，周向角θ=0°~90°、270°~360°为背风区，周向角θ=90°~270°为迎风区。沿x、y、z轴分别为轴向力、法向力、侧向力，绕x、y、z轴分别为滚转力矩、侧向力矩、俯仰力矩。可知侧向力系数C_Z即为马格努斯力系数，偏航力矩系数C_MY即为马格努斯力矩系数。

2 计算结果与分析 2.1 马格努斯力和力矩系数数值模拟结果

图 3给出了迎角α=-5°~40°动态气动特性随迎角变化规律，并与AEDC实验^[19]及美国陆军研究实验室(ARL)计算结果^[20]进行了对比。从图中可知侧向力矩旋转导数C_MYΩ随迎角增加呈减小趋势，侧向力旋转导数C_ZΩ随迎角增加先减小后增大，在a=14°附近出现拐点，拐点附近的计算值比实验值略小。总体而言本文计算结果与实验及文献结果吻合很好，验证了数值方法的可靠性。

对比SA、SST湍流模型及修正后的SARC、SSTRC模型可知，弯曲和旋转修正后的计算结果基本与原始模型重合，在实验和ARL结果的分散区间内。侧向力和力矩旋转导数差异主要在α=5°~20°，侧向力矩旋转导数SA偏大，与SARC最大偏差5.6%；侧向力旋转导数SA偏小，与SARC最大相差4.5%，由于5°~20°迎角范围侧向力和力矩旋转导数在零值附近，绝对差量很小，在图中体现不明显。

2.2 工程估算与数值模拟结果对比

工程估算方法是通过风洞实验数据或数值模拟数据，整理归纳出具有理论解析表达式的形式，应用于工程设计中。早期的学者对超声速旋成体马格努斯特性总结了很多工程估算方法，得到了很多规律性结论，如马格努斯力矩导数随马赫数增加(Ma=1-2.5)呈线性变化^[21]，随船尾角增加呈非线性变化等^[22]，对于本文这种带翼外形，选取四组典型估算公式进行分析，以考量工程估算方法的适用性。

1) Martin公式：

$ C{z_\mathit{\Omega }} = \frac{{52.6}}{{Re_l^{1/2}}}{\left( {\frac{L}{D}} \right)^2}\alpha $

(8)

2) Kelly和Thacker公式：

$ C{z_\mathit{\Omega }} = \frac{{63.68}}{{Re_l^{1/2}}}{\left( {\frac{L}{D}} \right)^2}\left[ {1 - 6.36{\mathit{\Omega }^2}{{\left( {\frac{L}{D}} \right)}^2}} \right]\alpha $

(9)

3) Vaughn和Reis公式：

$ \begin{array}{l} C{z_\mathit{\Omega }} = \frac{{20.8f\left( {T, Ma, \mathit{\Omega }} \right)\left( {1 - 1.52\left| \alpha \right|} \right)}}{{Re_l^{1/2}}}{\left( {\frac{L}{D}} \right)^2}\\ \;\;\;\;\;\;\;\;\;\;\left\{ {1 - \left[ {1 - \frac{{0.231}}{{1 - 1.52\left| \alpha \right|}}} \right]} \right\}{\left( {\frac{{{L_h}}}{L}} \right)^{3/2}}\alpha \end{array} $

(10)

4) 吴承清提出的公式：

$ C{z_\mathit{\Omega }} = \frac{A}{{Re_l^{1/2}}}{\left( {\frac{L}{D}} \right)^2}{\eta _M} \cdot {\eta _t} \cdot {\eta _{\theta \alpha }} $

(11)

式(8)-(11)中的各参数定义见文献[20-22]，图 4给出了工程估算方法与本文计算的侧向力旋转导数C_ZΩ随迎角变化曲线，可以发现：1)工程算法中C_ZΩ随迎角α多接近线性分布，数值计算的C_ZΩ随迎角α非线性变化未能体现；2)数值模拟与工程估算结果在整个迎角范围α=5~40°差异都很大，不满足工程估算公式可用性标准(误差 < 15%)，且随着迎角的增加差异越来越大；3)分析差异原因可能是工程估算多是针对旋成体(尖拱+圆柱、尖拱+圆柱+船尾)外形，对于本文翼-身组合外形，弹翼的存在引入很大误差，弹翼的影响不可忽略，导致工程估算方法失效。

2.3 沿轴向弹身表面压力及摩阻分布

图 5给出了t=T、α=5°、40°背风区典型周向角下弹身表面沿流向的摩阻分布，从图中可知α=5°时摩阻系数变化平缓，系数范围[0.75，4]，从头部到尾翼前缘逐渐减小，只在肩部有一次转折，翼-身干扰使尾翼所在截面摩阻略有增加；α=40°时摩阻抖动很大，系数范围[-0.2，1]，大迎角和旋转使弹身表面流场变化剧烈，摩阻系数无明显规律。

比较两湍流模型可知，SA计算的摩阻系数比SST偏大；比较修正模型与原始模型可知，左右两侧摩阻系数差异都不大，只在弹体尾部区域有微小差异，而且摩阻本来就是小量，对全弹的影响很小。

图 6给出了t=T、α=5°、40°背风区典型周向角下弹身表面沿流向的压力分布，从图中可知α=5°时从头部到尾部为膨胀-压缩-膨胀过程，p/p_∞范围[0.75，1]；而α=40°时从头部到尾部先膨胀-压缩，随后恒压段，最后再膨胀-压缩，由于弹体对来流的遮挡效应更强，压力系数明显减小，p/p_∞范围[0.18，0.35]。

比较修正模型与原始模型可知，α=5°时SSTRC模型比SST模型计算的压力略小；α=40°时SARC模型比SA模型计算的压力稍大，差异主要在弹身后部及尾部，由于背风区压力量值很小、左右两侧压力抵消作用，全弹侧向力系数原始模型与改进模型差异不大。

2.4 沿周向截面和对称面空间流场分布

图 7给出了t=T、α=40°、x/D=8截面SST和SSTRC模型空间压力云图及流线分布。从压力分布可知头部激波向后发展，在弹身尾部迎风区存在很大的弧型激波，背风区为膨胀区域，压力云图左右两侧差异不明显。而流线分布不再对称，正向旋转使弹体左侧速度向下，而正迎角时来流速度分量向上，两者相互阻碍使流线更脱体，易于分离；右侧旋转速度与来流同向，流线更贴体，不易分离，分析认为这是分离涡左侧比右侧大的原因。对比两模型结果可知，背风区SSTRC捕捉到的分离涡比SST大，可见弯曲和旋转修正使模型对分离流动的抑制能力减弱。

图 8给出了t=T、α=40°、对称面SA和SARC模型计算的马赫数云图分布。可以看到两模型模拟到的头部激波、背风区膨胀波及舵上二次压缩激波基本相同，主要差异在弹体底部，SARC计算的底部低速区更大，即修正模型预测的底部分离区更大，恢复区更长。

3 结论

本文采用完全时间相关的非定常N-S方程，对超声速带翼旋转弹箭开展计算，研究了弯曲和旋转修正的湍流模型SARC和SSTRC对弹箭旋转气动特性和流场结构产生的影响，通过数据分析得出以下结论：

1) 对全弹侧向动态特性而言，弯曲和旋转修正的湍流模型与原始模型精度相当，侧向力和力矩旋转导数最大差异 < 6%，工程估算方法不适用于本文翼身组合外形。

2) 对流场结构而言，修正模型使物面压力左右两侧同时偏大或偏小，与原始模型相比并没有加剧或削弱不对称效应，因此全弹马格努斯特性变化不大，弯曲和旋转修正模型预测的分离区更大，对分离流动的抑制能力减弱。

3) 修正模型与原始模型结果差异不大的分析有两点：一是原始模型预测结果已经很好，改进模型提升效果不明显；二是当前计算工况下逆压梯度引起的流线弯曲占主导，大于旋转和弯曲效应影响。修正模型对旋转弹箭的适用性还需做更多的研究。