基于SST湍流模型的模拟SRM内流场数值仿真①

李映坤，韩珺礼，2，陈 雄，沈振华

(1．南京理工大学机械工程学院，南京 210094;2．北京机电研究所，北京 100012)

0 引言

准确计算固体火箭发动机的内流场，对于固体火箭发动机的结构设计与优化有着非常重要的作用。固体火箭发动机工作过程中的燃气流动大多呈湍流，推进剂表面的燃烧反应会受到湍流流动的影响，燃烧室内的湍流流动不仅对推进剂的燃烧产生直接影响，而且对喷管内跨声速流动及喷管外羽流产生影响［1］。

为研究固体火箭发动机燃烧室内流场，1990年Dunlap R等对直径为101．6 mm、长度为1 453 mm的燃烧室进行了冷流实验［2］，利用氮气流入含有444个圆孔的圆柱表面模拟发动机推进剂表面燃烧时的侧向加质，对诸如径向速度、脉动速度等重要物理参量进行了测量，为燃烧室内流场的数值研究［3－6］提供了重要参考，但尚未将 Menter F R提出的 k-ω SST(shearstress-transport)湍流模型［7－8］应用在结构网格固体火箭发动机燃烧室内流场数值模拟中的报道。

Menter F R提出的k-ω SST(shear-stress-transport)湍流模型，能充分发挥k-ε模型对自由流和k-ω模型对壁面受限流动的处理优势，在空气动力学［8］、航空发动机［9－10］等领域有着广泛的应用。

本文基于有限体积法求解Navier-Stokes流动控制方程组，将该湍流模型应用于固体火箭发动机内流场的数值模拟，并将计算结果与实验值及Spalart-Allmaras湍流模型、Wilcox的k-ω两方程湍流模型进行了对比，为固体火箭发动机内流场仿真过程中湍流模型的选择提供了参考。

1 控制方程

1．1 流体控制方程

积分形式的可压缩非定常Navier-Stokes方程为

式中 ρ为气体密度;u、v为流体运动速度矢量的2个分量;T为气体的温度;E为单位体积气体的总能量;τ表示应力张量，其具体形式参考文献［11］。

1．2 湍流模型

为便于对比分析，本文共采用了4种湍流模型。Spalart-Allmaras湍流模型［12］基于经验和量纲分析，考虑了自由剪切流、高雷诺数时的近壁区流动、有限雷诺数的近壁区流动等，直接针对涡粘性建立的方程。

Wilcox的k-ω两方程湍流模型［13］是应用最广的两方程涡粘性模式之一，为积分到壁面的不可压缩或可压缩湍流两方程涡粘性模式，该模式不需要显示的壁面衰减函数。

Menter提出的 k-ω SST(shear-stress-transport)剪切应力输运模型(以下简称为SST k-ω湍流模型)，该模型通过混合函数F1将k-ε模型和k-ω模型结合起来，这样充分发挥了k-ε模型对自由流和k-ω模型对壁面受限流动的处理优势。具体描述如下:

其中

式中 k为湍动能;ω为比耗散率;μt为湍流粘性系数，其他参数的具体形式见参考文献［7］。

目前，经性传播已成为我国艾滋病传播的主要方式，而家庭内配偶间经性传播已成为艾滋病进一步蔓延的重要因素之一，我国2011年估计的78万艾滋病患者中经异性传播占46.5%，其中约1/4为配偶间性传播［1］。因此，了解配偶间人类免疫缺陷病毒（human immunodeficiency virus，HIV）传播状况及其相关影响因素，采取相应措施降低配偶间HIV传播尤为重要，现将相关研究进展综述如下。

2003年，为减小原湍流模型对近壁处网格尺寸的依赖性［8］，Menter F R对原始的湍流模型进行了改进(以下简称为SST k-ω-2003湍流模型)，其改进后形式如下:

2 数值求解方法

本文自编程序，采用基于格心的多块结构网格迎风型有限体积法，求解上述Navier-Stokes流动控制方程组，对流项的离散是计算的关键，对于图1所示的控制体ABCD，流过AB边通量的计算步骤为

(1)采用具有保单调性的三阶MUSCL(Monotone Upstream-centered schemes for Conservation Laws)数值格式，利用界面左右两侧4个点的信息，计算出控制体界面AB处的物理量(密度、速度、压力等)，并使用Van abada限制器，以避免间断处的非物理震荡，具体形式如下:

其中，S为Van Albada限制器，具体形式为

式中 ε=10－6为一小量，防止上式分母为0。

(2)利用控制体界面处的物理量，采用AUSM-PW通量技术［14］，计算出单位时间内通过控制体边界的流通量，以i方向上的通量F为例，I+1/2界面上的通量可写为

式中 c为单元界面声速;Φ为守恒通量;p为压力项。

(3)利用积分型的守恒方程，计算出下一个时间步控制体内物理量的平均值。

图1 控制体示意图系Fig．1 The diagram of a control volume

粘性项采用Jameson中心差分法离散，时间推进采用三阶三步TVD型Runge-Kutta显式方法，此方法计算过程简单，内存需求小，但时间步长受稳定性条件限制，而必须取得很小。因此，本文采用局部时间步长加速收敛技术来加快定常流场计算的收敛速度。该方法的基本思想是定常问题在求解过程中，时间和空间离散是非耦合的，流场最终的计算结果不受各点的发展历程的影响，最终的定常计算结果与所采用的时间步长无关。因此，可适当地改变计算的中间过程，以加快得到最终定常解的目的。

3 结果与分析

为说明SST湍流模型模拟固体火箭发动机内流场的能力，本文选择Dunlap R的圆柱冷流实验模型作为对比，该实验利用氮气流入含有444个圆孔的圆柱表面，模拟发动机推进剂表面典型的入射速度和雷诺数，Dunlap R对此模型进行了大量的研究，得到了许多有用的实验数据，本文将这些实验数据与数值计算结果进行了对比。

该实验模型包括一个等截面的圆柱和一个收敛扩张的喷管，如图2所示，圆柱段的长度L=1 453 mm，直径D=101．6 mm，且L/D=14．3。喷管的收敛扩张段之间有一段半角为20°的圆弧，且该圆弧半径等于喷管喉部半径，扩张段下游的形状由一段抛物线段确定，喷管出口直径与圆柱段直径相等。因此，该喷管的扩张比为 2．27。

图2 冷流实验模型及边界条件Fig．2 of the cylindrical-port cold flow model and boundary conditions

根据文献［2］的实验，按照文献［4－6］的处理方法，忽略圆柱壁面上圆孔之间的间隙，将整个圆柱表面作为入口边界，圆柱表面注入氮气的速度 uinj=0．001 8 Ma，温度Tinj=303 K，基于圆柱半径rw和注入速度的雷诺数为Reinj=180 00。文献［4］的计算结果表明，入口处湍动能k和湍动能比扩散率ω的取值对计算结果有很大的影响。因此，本文根据文献［15］的方法，采用式(24)计算入口处的k和ω:

式中 σv=2．5 ×10－6;lw=1．0。

固体壁面假设为绝热壁，采用无滑移边界条件;出口边界根据马赫数判定，当出口为超声速时，此时所有物理量外推，当出口为亚声速时，给定环境反压101 325 Pa，其他参数由内向外插值。

3．1 网格影响分析

计算区域网格划分的疏密对数值计算有非常大的影响。因此，有必要开展网格相关性研究，排除网格对计算结果的影响。本文采用多块对接网格，如图3所示。

图3 模拟固体发动机计算网格Fig．3 The computational grid of model SRM

网格分为粗网格(454×60)、中等网格(539×84)和细网格(654×120)3种，壁面第一层网格的间距为1．0×10－5m。图4是采用 SST k-ω 湍流模型计算的x/D=10．30处轴向无量纲速度沿着径向的分布(图中轴向速度u以中心处速度uc无量纲化，燃烧室半径rw=50．8 mm)，D为燃烧室直径，x为距离燃烧室头部的距离。从图4可看出，在中心线附近粗网格与实验值吻合得很好，但在靠近壁面附近，粗网格与实验值相差很大，而中等网格与细网格整体上都与实验值吻合得很好，尤其是在靠近中心线处，中等网格计算的结果几乎与细网格一致。因此，本文后面的计算均采用中等网格。

图4 x/D=10．30处轴向速度网格影响分析Fig．4 Grid independence study for axial mean velocity at x/D=10．30

3．2 轴线上马赫数分布

采用不同的湍流模型计算得到的马赫数沿着轴向的分布如图5所示。

图5 不同湍流模型计算的马赫数沿着轴线的分布Fig．5 Comparison of computed and experimental axial mach number with different turbulent models

由图5可见，马赫数沿着轴线呈线性增加趋势，在x/D小于7．5时，SST k-ω-2003湍流模型与实验值吻合得非常好;其次是 SST k-ω湍流模型，Wilcox k-ω和Spalart-Allmaras湍流模型与实验值相差较大;当x/D超过7．5时，4种湍流模型的计算结果均与实验值有一定差别。这是因为实验中气体的流动状态迅速转变为湍流，而数值模拟中这一过程却相对较慢，但SST kω-2003湍流模型与实验值的误差是最小的，Spalart-Allmaras湍流模型次之。

3．3 不同位置处轴向速度对比

图6是采用不同的湍流模型计算的轴向速度与实验值的对比，沿着轴线图中，共给出了10个位置处的轴向速度沿着径向分布(图中轴向速度u以中心处速度uc无量纲化)。整体上来看，本文计算的结果与实验值的趋势是一致的，但在x/D=0．62处，4种湍流模型计算的结果都小于实验值，0．5≤r/rw≤1．0范围内，Spalart-Allmaras湍流模型计算的结果最接近实验值，但仍有很大的差距，靠近中心线附近SST k-ω-2003湍流模型计算的结果几乎与实验值一致。造成这一差异的原因是x/D=0．62处燃烧室前端轴向速度较小。此时，燃烧室表面注入的径向气流正在向轴向方向转变，如图7所示，燃烧室前端x/D=0．62处，最大轴向速度仅为3．07 m/s，而燃烧室尾部x/D=12．72处的最大轴向速度是前端的17倍。同时，文献［2］的实验结果表明，此处形成了比较强的涡流，并沿着轴线，涡流向燃烧室中心附近汇聚，湍流的效果越来越明显。所以，这就解释了x/D=1．80处的计算结果相比于x/D=0．62处更加接近实验值的原因。另外，从图7还可看出，靠近壁面附近SST k-ω-2003湍流模型计算的结果明显优于其他几个湍流模型，Wilcox k-ω湍流模型的结果与实验值相差甚远;在中心轴线附近，靠近燃烧室尾部的位置(r/rw≥7．88)，SST k-ω湍流模型的结果与实验值吻合得很好，而SST k-ω-2003湍流模型计算的结果相比实验值略微偏小，在x/D=7．88位置处达，到最大误差，约为5．1%。产生这一差异的原因为燃烧表面注入的径向气流往下游移动，流动状态从层流转变为湍流，相比径向与周向运动速度分量，轴向速度分量占据主导地位;同时，文献［2］中的实验分析表明，靠近燃烧室尾部处，涡量被限制在中心线附近。因此，基于涡量的SST k-ω湍流模型的计算值与实验值吻合得很好，而基于剪切率的SST k-ω-2003湍流模型的计算值与实验值略微有点差距，但总体趋势还是一致的。

3．4 不同位置处湍流强度对比

基于Boussinesq涡粘性假设的两方程湍流模型很难准确预测燃烧室内的湍流强度分布［6］，而大涡模拟方法和直接数值模拟方法预测精度较高，但其计算量较大，以目前的计算机资源很难在工程上广泛应用。本文计算的不同位置处的湍流强度I沿着径向的分布如图8所示。由图8可见，几个湍流模型计算的结果与实验值的趋势是一致的，但都有一定的差异，特别是Wilcox k-ω湍流模型与实验值相差很大，而SST湍流模型计算的结果更接近于实验值，尤其是 SST k-ω-2003湍流模型在x/D=10．30位置处，几乎与实验值是重合的，但在加质壁面附近，SST k-ω-2003湍流模型计算的结果与实验值还略微有些差距。文献［3］认为，这与加质壁面湍流边界条件的取值有关，而在中心线附近，SST k-ω-2003湍流模型的计算结果与实验值吻合得很好。

图6 不同位置处轴向速度与实验值对比Fig．6 Comparison of computed and experimental axial velocity profiles

图7 燃烧室前、后端速度沿径向的分布Fig．7 Velocity distribution at the head end and the aft of combustion chamber

由图8进一步分析可知，靠近燃烧室前端湍流强度很小，随着气流向下游移动，湍流强度逐渐增大，本文SST k-ω-2003湍流模型也成功地模拟出了这一趋势，相比于 x/D=3．04位置，x/D=5．46处湍流强度明显增大，在x/D=6．64处增大得更加明显，这说明在此位置，流动正在由层流转变为湍流。

4 结论

(1)采用不同湍流模型计算得到的马赫数沿轴线均呈线性增加趋势，当x/D小于7．5时，SST k-ω-2003湍流模型与实验值吻合得非常好;当x/D超过7．5时，4种湍流模型的计算结果均与实验值有一定差别，但SST k-ω-2003湍流模型与实验值的误差最小。

(2)4种湍流模型都能准确模拟出固体火箭发动机燃烧室内的径向速度分布，且计算结果差别不大，SST k-ω-2003湍流模型的计算结果与实验值吻合得最好，最大误差约为5．1%。

(3)SST k-ω-2003湍流模型计算的固体火箭发动机燃烧室内湍流强度分布与实验的规律一致，而其余湍流模型计算的结果与实验值有很大差异。因此，与其他湍流相比，SST k-ω-2003湍流模型能较准确地模拟固体火箭发动机的内流场，具有应用于固体火箭发动机内流场仿真的可行性和相对优势。

图8 不同位置处湍流强度与实验值对比Fig．8 Turbulence intensity in vertical direction at various axial locations

［1］ 武晓松，陈军，王栋，等．固体火箭发动机工作过程数值仿真［M］．北京:高等教育出版社，2005．

［2］ Dunlap R，Blackner A M，Waugh R C，et al．Internal flow field studies in a simulated cylindrical port rocket chamber［J］．Journal of Propulsion and Power，1990，6(6):690-704．

［3］ Sabnis J S，Gibeling H J，McDonald H．Navier-Stokes analysis of solid propellant rocket motor internal flows［J］．Journal of Propulsion and Power，1989，5(6):657-664．

［4］ Arabshahi A，Webster R S，Briley W R，et al．Numerical analysis of solid propellant rocket motor internal flows［R］．AIAA 2006-5114．

［5］ Arabshahi A，Sreenivas K，Nichols D S，et al．Computational analysis of turbulent internal flow in ballistic solid rocket motors［R］．AIAA 2007-1449．

［6］ Yumusak Mine．Analysis and design optimization of solid rocket motors in viscous flows［J］．Computers ＆ Fluids，2013，75:22-34．

［7］ Menter F R．Two-equation eddy-viscosity turbulence models for engineering applications［J］．AIAA Journal，1994，32(8):1598-1605．

［8］ Menter F R，Kuntz M，Langtry R．Ten years of industrial experience with the SST turbulence model［J］．Turbulence，Heat and Mass Transfer，2003(4):625-632．

［9］ 姚玉，张靖周，郭文．气膜孔角度对导叶冷却效果影响的数值研究［J］．航空动力学报，2009，24(3):507-512．

［10］ 曾军，王彬，卿雄杰．某双级高压涡轮全三维计算［J］．航空动力学报，2012，27(11):2553-2561．

［11］ 张涵信，陈坚强，高树椿．H2/O2燃烧的超声速非平衡流动的数值模拟［J］．宇航学报，1994，15(02):14-23．

［12］ Spalart P R，Allmaras S R．A one-equation turbulence model for aerodynamic flows［R］．AIAA 92-1439．

［13］ Wilcox D C．Turbulence modeling for CFD［M］．La Canada:DCW Industries，1998．

［14］ Meng-sing Liou．Methods for the Accurate Computations of Hypersonic Flows［J］．Journal of Computational Physics，2001，174(1):38-80．

［15］ Sachdev J S．Parallel solution-adaptive method for predicting solid propellant rocket motor core flows［D］．Doctor of Philosophy Graduate Department of Aerospace Science and Engineering University of Toronto，2007．