第九章流动过程数值分析

发表于： 2019年11月16日 2021年4月30日
分类：固体火箭冲压发动机

被阅次数： 1,235

9.1 补燃室掺混反应流场数值计算

9.1.1 补燃室数值计算的现状

在固体火箭冲压发动机补燃室中，贫氧燃气从头部特定位置以亚音速或（超）音速喷入，与侧边或头部进入的空气进行掺混和燃烧，它们与燃烧产物一起形成复杂的三维流场，在流动过程中伴随着热量的释放和成份的急剧变化。补燃室结构、燃气性质、空燃比等都强烈地影响着补燃室内的流场分布，仅凭经验、半经验设计方法不能满足现代先进燃烧室的设计要求，加上燃烧室的加工和试验费用十分昂贵，迫切需要发展一种新型的设计计算方法。这种计算方法的特点是以计算流体力学（CFD）、数值传热学(CHT)、计算燃烧学(CCD)为核心，结合经验、半经验关系式，进行数值模拟计算，而后把计算结果经过试验验证。

燃烧室数值计算在近30年取得了重要进展。1976年，Spalding等人编制了PHOENICS（Parabolic，Hyperbolic,or Elliptic Numerical Integration Code System）大型通用计算机程序，并于1981年投放市场。PHOENICS软件可以用来模拟传热、流动、化学反应及燃烧过程。该软件中所采用的一些基本算法，如SIMPLE方法、混合格式等，正是由该软件的创始人及其合作者所提出的，对以后开发的商用软件有较大的影响。由于该软件投放市场较早，因而在工业界得到较广泛的应用，在算例库中收入了600多个例子。近年来它的功能不断得到改善。在燃烧界应用比较广泛的还有fluent公司研制的FLUENT软件，它于1983年投放市场。该软件在网格生成，物理模型，数值方法，前后置处理方面吸收了当代计算流体力学许多先进技术，如非结构网格，自适应网格，大涡模拟，并行处理等。可计算的物理问题类型有：定常与非定常流动，可压与不可压流动（对高Ma下的流动，专门有RAMPANT软件），含有粒子/液滴的蒸发、燃烧的过程，多组分介质的化学反应过程等。其他可用于燃烧过程数值模拟的商业软件还有CFX,STAR-CD,CFDRC等。勿庸置疑，这些商业软件在一定程度上都可以用来对燃烧室中气动热力过程进行数值模拟。然而，由于燃烧过程的复杂性，目前这些商业软件距离先进燃烧室设计方法的要求尚有很大一段距离。由于任何商业CFD软件都不可能提供先进燃烧室设计所必需的技术支持，因此世界各大发动机公司都投入巨大的人力、物力、财力开发自己的燃烧室设计软件。如NASA-Lewis研究中心使用SIMPLE法编制而成的大型通用程序TEACH（Teach Elliptic Axisymmetric Characteristics Heuristically）用来计算各种燃烧流动问题。

对固体火箭冲压发动机内部燃烧过程进行分析研究的文献，最早见于上世纪70年代，在此之前，对固体火箭冲压发动机内部的燃烧过程所知甚少。在补燃室数值计算的初期，主要是把补燃室内的掺混燃烧简化为一种气体的燃烧，如文献【51,52,53】等；由于一次燃气中含有大量凝聚相粒子，有的学者对补燃室内的两相反应流场进行了数值模拟，如文献【54】等；为了追踪凝聚相粒子的燃烧过程，有的学者用颗粒轨道模型对粒子轨迹进行追踪，分析凝聚相粒子的燃烧效率，这方面的文献如【55,56,57】等。在我国，随着固体火箭冲压发动机的研制，对相应的反应流场也开展了一些数值研究，绝大部分研究者都采用了当今比较流行的CFD软件，有的研究者对补燃室反应流场的数值计算方法进行了研究，如【58】。

对补燃室进行数值计算的主要功能有以下几点：

（1）模拟燃烧过程，如点火、火焰稳定，燃烧室中的气流结构，浓度场，温度场，金属粒子的运动、燃烧等物理化学过程；

（2）预估不同条件下的燃烧室性能，如总压损失，燃烧效率，壁温等；

（3）用于燃烧室优化设计，在初步设计阶段用于方案选择，在技术设计阶段用于产品性能评估与定型；

（4）指导研制试验，减少试验次数，缩短燃烧室研制周期；

（5）对燃烧室中复杂的物理化学现象提供更深刻的认识，导致新的设计概念。

9.1.2 湍流反应流场的计算模型

在自然界及工程实际中，遇到的流动几乎都是湍流。例如，地球大气附面层的流动，各种燃烧室和锅炉内的流动，江河水流等都是湍流。作为现代流体力学和水力学重要组成部分的湍流，对许多科学技术领域，如航空、航天、造船、气象、化工、环保、冶金等都是有直接或间接的关系。正因如此，近百年来，许多学者对湍流问题进行了大量试验与理论研究，并取得了丰硕的成果，但由于湍流过程十分复杂，目前对湍流结构、机理等许多理论问题尚未完全解决。为了解决工程问题，必须研究湍流对平均流场的影响，计算在各种复杂流动条件下的湍流平均流场。

9.1.2.1 湍流流动的性质

湍流的基本特征在于其具有随机性质的涡旋结构以及这些涡旋在流体内部的随机运动，从而引起流速、压力、温度等各种参数的脉动。这些涡旋形成以后，经历着发展和衰减的过程，直到消失。大尺度的涡旋破坏后形成尺度较小的涡旋，较小的涡旋破坏后形成更小的涡旋。这些尺寸大小不等的涡旋往往是大涡旋包含着小涡旋，小涡旋中包含着更小的涡旋。各种不同尺度的涡旋充满着整个湍流，组成连续的涡旋谱，大涡旋的尺度与湍流条件密切相关，它们决定了湍流的主要力学特征。而小尺度的涡旋在黏性影响下将湍流能量转换为热而消失。由于这些旋涡在湍流内部作随机运动，不断平移和转动，使得湍流各点的速度随时间不断变化。一个涡流的转动，常常会引起相邻流体相反方向的转动；一个涡旋离开后，必有另一个涡旋来补充。因而尽管湍流是由各种不同尺度涡团组成，但瞬时湍流仍然符合连续条件，湍流内部各点瞬时流速仍然可认为是连续的。一般来说。湍流的基本性质有以下几个方面：

（1）随机性

在湍流流动中所有参数都是随时间和空间而变化的随机量，但是它们仍具有一定规律的统计学特征，因而它们在空间某一点上既有平均值，又有脉动值，脉动是旋涡气团不规则的乱动。因此可认为各气流参数的瞬时值等于平均值加脉动值，例如：

$u=\bar{u}+{u}’$ $p=\bar{p}+{p}’$ $T=\bar{T}+{T}’$

其中$\bar{u}$、$\bar{p}$、$\bar{T}$为变量$u$、$p$、$T$的时均值，${u}’$、${p}’$、${T}’$为脉动值。脉动值是随机变量，但所有脉动值的时均值应为零。

（2）湍流强度

湍流流动中各参数的脉动值的变化是不规则的，由于所有脉动值的时均值为零，故常用其均方根值来表示脉动量的大小，均方根值可定义为：

${{\left( \overline{{{{{u}’}}^{2}}} \right)}^{{1}/{2}\;}}={{\left( \int_{{{t}_{0}}}^{{{t}_{0}}+{{t}_{1}}}{{{{{u}’}}^{2}}dt} \right)}^{{1}/{2}\;}}$ （9.1-1）

其中${{t}_{1}}$为时间间隔，它应取较大值，至少应大于脉动量的最大周期。

在工程实际中，把湍流强度定义为：

$I\equiv \frac{{{\left( \overline{{{{{u}’}}^{2}}} \right)}^{{1}/{2}\;}}}{{{u}_{\infty }}}$ （9.1-2）

其中${{u}_{\infty }}$为气流平均速度或自由流速度，${{\left(\overline{{{{{u}’}}^{2}}} \right)}^{{1}/{2}\;}}$为湍流速度。对一般管内流动湍流强度约为5％，在火焰稳定器后尾流区内湍流强度是变化的，大约为10％。

（3）时间相关函数

湍流流动中在同一点的两个脉动参数可能互相关联，例如$\overline{{u}'{v}’}\ne 0$，则称${u}’$与${v}’$相关联；如果$\overline{{u}'{v}’}=0$，则称${u}’$与${v}’$不相关联。相关联的两个脉动参数的脉动形式应大致类似，否则不可能相关联。在统计理论中，把两个脉动分量乘积的平均值$\overline{{u}'{v}’}$

称为二阶相关矩，通常可用它来描述湍流内部结构。

两个脉动参数${u}’$和${v}’$的相关程度，可用相关函数R来表示：

${{R}_{E}}\left( \tau \right)=\frac{\overline{{u}’\left( t \right){u}’\left( t+\tau \right)}}{\overline{{{{{u}’}}^{2}}}}$ （9.1-4）

它表示同一脉动量在同一空间点，但在不同时间测出量之间的关联。由时间相关函数可导出湍流的时间尺度。

（4）空间相关函数

为了更好地了解湍流结构，需要了解空间一点与其周围脉动的关系，为此必须处理脉动量在不同空间点之间的关联。假设在同一时间，距离为的两点脉动参数分别为${u}’\left( x \right)$和${u}’\left( x+r \right)$，取两者的乘积的时均值得：

$\overline{{u}'(x){u}'(x+r)}=\frac{1}{{{t}_{1}}}\int_{{{t}_{0}}}^{{{t}_{0}}+{{t}_{1}}}{{u}'(x){u}'(x+r)}dt$ （9.1-5）

与时间相关函数一样，空间相关函数可定义为：

$R\left( r \right)=\frac{\overline{{u}'(x){u}'(x+r)}}{{{\left[ \overline{{u}'{{\left( x \right)}^{2}}} \right]}^{{1}/{2}\;}}{{\left[ \overline{{u}'{{\left( x+r \right)}^{2}}} \right]}^{{1}/{2}\;}}}$ （9.1-6）

由空间相关函数可导出湍流长度尺度。

（5）扩散性

湍流扩散比层流扩散快得多，因为层流扩散是靠分子的迁移，而湍流扩散是靠湍流气流的迁移。由于湍流气团的尺寸要比分子大得多，所以湍流扩散系数也比层流扩散系数大得多。此外，湍流扩散系数不是流体的一种性质，而是取决于流场的各种因素。湍流扩散性能强的优点是可使混合过程加快和完善。

（6）湍流运动都是三维有旋

湍流中有大小不同尺度的涡旋气团，既有大气团的宏观运动，又有小气团的微观乱动，由此产生不规则的速度脉动，使湍流流动具有三维性，并有一定强度涡旋和涡线，其中涡旋不断产生和消失。

以上这些湍流的基本性质已经成为现代湍流理论的基础。但随着对湍流的深入研究，近年来对湍流的结构有了一些新的认识。为了突出重点和节省篇幅，本书对此理论不作过多介绍。

9.1.2.2 湍流流动的模拟方法

关于湍流流动与换热的数值计算，是目前计算流体力学与计算传热学中困难最多也是研究最活跃的领域之一。目前采用的数值计算方法大致可分三类：

（1）直接模拟

直接数值模拟方法是直接求解瞬时湍流控制方程。要对复杂的湍流流动进行直接的数值计算，必须采用很小的时间和空间步长，这样才能分辨出湍流中详细的空间结构及变化剧烈的时间特性。直接模拟对内存空间及计算速度要求非常高，目前只有少数能使用超级计算机的研究者才能从事这类计算，无法用于工程计算。

（2）大涡模拟

由湍流的性质可知，湍流的脉动与混合主要是由大尺度的涡造成的。大尺度的涡从主流中获得能量，大尺度的涡旋破坏后形成尺度较小的涡，较小的涡旋破坏后形成更小的涡旋。小尺度涡的主要作用是耗散能量，它们几乎是各向同性的。关于涡旋的上述认识就导致了大尺度涡模拟的数值方法。这种方法旨在用非稳定的N－S方程来直接模拟大尺度涡，但不直接计算小尺度涡，小涡对大涡的影响通过近似的模型来考虑，这种影响称为亚格子雷诺应力。大多数亚格子雷诺应力模型都是在涡黏性的基础上，即把湍流脉动所造成的影响用一个湍流黏性系数来描述。大涡模拟方法对计算机内存即速度要求虽仍然较高，但远低于直接模拟方法对计算资源的要求，在工作站或PC机上都可以进行一定的研究工作^【59】。

（3）应用雷诺时均方程的模拟方法

在这类方法中，将非稳态控制方程对时间作平均，在所得出的关于时均物理量的控制方程中包含了脉动量乘积的时均值等未知量，于是所得方程的个数就小于未知量的个数。要使方程组封闭，必须作出假设，即建立模型。这种模型把未知的更高阶的时间平均值表示成较低阶的计算中可以确定的量的函数，这是目前工程湍流计算中所采用的基本方法。

在雷诺时均方程法中，又有雷诺应力方程法及湍流黏性系数法两大类，其中的湍流黏性系数法是目前工程流动与数值计算中应用最广的方法。所谓湍流黏性系数法，就是把湍动黏度${{\mu }_{t}}$与湍流时均参数联系起来的关系式。依据确定${{\mu }_{t}}$的微分方程数目的多少，湍流黏性系数法包括零方程模型、一方程模型、两方程模型。目前两方程模型在工程中使用最广泛，最基本的两方程模型是标准$k-\varepsilon $模型，此外还有各种改进的$k-\varepsilon $模型，比较著名的是RNG $k-\varepsilon $模型和Realizable$k-\varepsilon $模型.

标准$k-\varepsilon $模型的优点是收敛稳定、精度适当、所需计算资源较少，因此该模型使用广泛。其缺点是对含有较大压力梯度或旋涡的复杂流场，计算结果与实际会产生偏差。为了消除标准$k-\varepsilon $模型的这些缺点，Yakhot及Orzag提出了RNG$k-\varepsilon $模型，该模型对较为复杂的流场，如气流冲击、流体分离、旋涡等较为适用。对于补燃室内的流动，可以采用RNG$k-\varepsilon $模型。

RNG$k-\varepsilon $模型的湍流动能方程和湍能耗散率方程分别如下：

$\frac{\partial }{\partial t}\left( \rho k \right)+\frac{\partial }{\partial {{x}_{i}}}(\rho k{{u}_{i}})=\frac{\partial }{\partial {{x}_{j}}}\left[ {{\alpha }_{k}}{{\mu }_{eff}}\frac{\partial k}{\partial {{x}_{j}}} \right]+{{G}_{k}}-\rho \varepsilon $ (9.1-7)

其中：${{\mu }_{eff}}=\mu +{{\mu }_{t}}$，${{\mu }_{t}}=\rho {{C}_{\mu }}\frac{{{k}^{2}}}{\varepsilon }$，${{C}_{\mu }}=0.0845$，${{\alpha }_{\text{k}}}=1.39$。

${{G}_{k}}$表示由层流速度梯度而产生的湍流动能，$\rho \varepsilon $为耗散项。

${{G}_{k}}={{\mu }_{t}}\left[ 2{{\left( \frac{\partial u}{\partial x} \right)}^{2}}+2{{\left( \frac{\partial v}{\partial y} \right)}^{2}}+2{{\left( \frac{\partial w}{\partial z} \right)}^{2}}+{{\left( \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \right)}^{2}}+{{\left( \frac{\partial w}{\partial x}+\frac{\partial u}{\partial z} \right)}^{2}}+{{\left( \frac{\partial v}{\partial z}+\frac{\partial w}{\partial y} \right)}^{2}} \right]$

$\frac{\partial }{\partial t}\left( \rho \varepsilon \right)+\frac{\partial }{\partial {{x}_{i}}}(\rho \varepsilon {{u}_{i}})=\frac{\partial }{\partial {{x}_{j}}}\left[ {{\alpha }_{\varepsilon }}{{\mu }_{eff}}\frac{\partial \varepsilon }{\partial {{x}_{j}}} \right]+C_{1\varepsilon }^{*}\frac{\varepsilon }{k}{{G}_{\text{k}}}-{{C}_{2\varepsilon }}\rho \frac{{{\varepsilon }^{2}}}{k}$ （9.1-8）

其中：

$C_{1\varepsilon }^{*}={{C}_{1\varepsilon }}-\frac{\eta \left( 1-\eta /{{\eta }_{0}} \right)}{1+\beta {{\eta }^{3}}}$

${{C}_{1\varepsilon }}=1.42$，${{C}_{2\varepsilon }}=1.68$

$\eta ={{\left( 2{{E}_{ij}}\cdot {{E}_{ij}} \right)}^{1/2}}\frac{k}{\varepsilon }$

${{E}_{ij}}=\frac{1}{2}\left( \frac{\partial {{u}_{i}}}{\partial {{x}_{j}}}+\frac{\partial {{u}_{j}}}{\partial {{x}_{i}}} \right)$

${{\eta }_{0}}=4.377$ $\beta =0.012$

除过上述描述湍流流动的模型外，还有雷诺应力微分方程模型，雷诺应力代数方程模型，其中目前已经用于工程数值计算的是二阶矩模型以及在此基础上经简化而得出的代数应力模型。关于这些模型的介绍可参考文献【59】。

9.1.2.3 湍流燃烧模型

燃烧现象是物理与化学过程相互作用的结果。燃烧反应的速率不仅与化学动力学因素有关，而且还与反应体系内能量的输运有密切关系。这种输运取决于系统的宏观运动以及系统内微观运动所引起的热传导和扩散现象。在实际燃烧过程中，燃烧速率往往受控与某些物理因素。对于混合气扩散速度比化学反应速度低得多的燃烧过程，称之为扩散燃烧。相反，对于混合均匀的气体，如果化学反应速率比扩散速率小得多的燃烧过程，称之为动力燃烧。

在湍流燃烧体系中，混气情况很复杂，不仅有气团的不规则运动，而且气团有大有小，气团中组分也不相同，有的含燃料多些，有的含氧多些，有的含高温燃烧产物多些。各种不同的气团在气流中互相混合、互相渗透，同时大气团不断转化为小气团；而大气团又通过雷诺应力对主气流作功后不断产生。在层流燃烧中，某种组分的产生率是由分子间化学反应造成的。分子间化学反应速率常用化学动力学中Arrhenius定律来表示。如燃料与氧气的反应速率为

${{R}_{fu}}=\rho \frac{d{{m}_{fu}}}{dt}={{A}_{0}}{{\rho }^{2}}{{m}_{ox}}{{m}_{fu}}\exp \left( -E/RT \right)$ （9.1-9）

式中，${{A}_{0}}$为前置因子，$E/R$为活化温度，随混气成分和温度而变，T为静温，$\rho $为混气密度，这里${{m}_{ox}}$、${{m}_{fu}}$、T和$\rho $

都是瞬时值。

在湍流燃烧流场中，化学反应速率与湍流流动和化学动力学因素有关。为了使基本方程得到封闭，必须求出平均化学反应速率，即

${{R}_{fu}}=\overline{{{A}_{0}}{{\rho }^{2}}{{m}_{ox}}{{m}_{fu}}\exp \left( -E/RT \right)}$ （9.1-10）

如果对式（9.1-10）进行雷诺平均，产生脉动值二阶相关项（如$\overline{{{{{m}’}}_{fu}}{{{{m}’}}_{ox}}}$、$\overline{{{{{T}’}}^{2}}}$、$\overline{{{{{m}’}}_{ox}}{T}’}$等）和三阶以及高阶相关项，计算${{R}_{fu}}$可通过对这些相关项加以模拟，使方程封闭。但这种方法涉及模拟量太多，而且高阶相关项如何计算至今也未完全解决。所以不得不采用近似计算方法，即找出影响${{R}_{fu}}$的主要因素，提出简化表达式，这就是所谓的湍流燃烧模型。

(1) 湍流扩散火焰模型

一般来说，燃料和氧化剂未预先混合的燃烧过程称为扩散燃烧，该火焰称为扩散火焰。实验表明，扩散火焰与层流火焰很不相同。以湍流射流火焰为例，其主要区别是，湍流扩散火焰高度与喷嘴出口速度无关，湍流脉动使火焰面变厚，出现噪声，在火焰区同时存在燃料与氧化剂的时均值，这并不意味着在同一地点同一时刻燃料与氧化剂并存，而是由于湍流脉动造成在同一地点某时刻出现过多的燃料，而在另一时刻出现过多的氧，可见在分析湍流扩散火焰时必须考虑湍流脉动的影响。

对于简化反应系统的单步不可逆反应：

$1kgSkg\to \left( \text{1}+S \right)kg$ (9.1-11)

定义$f\equiv \frac{{{\left( {{m}_{fu}}-{{{m}_{ox}}}/{S}\; \right)}_{M}}+{{{m}_{ox,A}}}/{S}\;}{1+{{{m}_{ox,A}}}/{S}\;}$，下标$M$表示混合气。$f$称为混合分数，物理意义是指1kg混合气中含有$f$kg未燃烧的燃料。定义的实际意义在于它可以度量混合进行的程度。当$f=0$时，混气中只有氧；$f=1$时，混气中只有燃料。当燃料与氧混合时，值在0与1之间。

分别列出燃料与氧的质量守恒方程，并利用定义的$f$，可将两个方程合二为一，即得到关于$f$\的方程：

$\frac{\partial }{\partial {{x}_{j}}}\left( \rho {{u}_{j}}f \right)-\frac{\partial }{\partial {{x}_{j}}}\left( {{\Gamma }_{j}}\frac{\partial f}{\partial {{x}_{j}}} \right)=0$ (9.1-12)

由于该方程的右端为0，因此是无源方程，f是个守恒标量。

由上述混合分数f的瞬时值方程，可得到时均值$\bar{f}$的守恒方程，它的脉动量的均方值$g$，可通过微分方程求得。

$\frac{\partial }{\partial t}\left( \rho \bar{f} \right)\frac{\partial }{{{x}_{j}}}\left( \rho {{{\bar{u}}}_{j}}\bar{f} \right)=\frac{\partial }{\partial {{x}_{j}}}\left( \frac{{{\mu }_{e}}}{{{\sigma }_{f}}}\frac{\partial \bar{f}}{\partial {{x}_{j}}} \right)$ (9.1-13)

$\frac{\partial }{\partial t}\left( \rho g \right)\frac{\partial }{{{x}_{j}}}\left( \rho {{{\bar{u}}}_{j}}g \right)=\frac{\partial }{\partial {{x}_{j}}}\left( \frac{{{\mu }_{g}}}{{{\sigma }_{g}}}\frac{\partial g}{\partial {{x}_{j}}} \right)+{{C}_{g1}}{{\mu }_{e}}{{\left( \frac{\partial \bar{f}}{\partial {{x}_{j}}} \right)}^{2}}-{{C}_{g2}}\rho g\frac{\varepsilon }{k}$ (9.1-14)

在湍流中，f是脉动的，它在火焰中脉动性质是不知道的，但它是个随机变量，可运用概率密度函数$P$($f$)描述其脉动性质，如果$P$($f$)已知，假设任意标量$\varphi $的瞬时值是f的单值函数$\varphi =\varphi \left( f \right)$，则任一标量的时均值和脉动均方根值，可由下列定义求出：

$\bar{\phi }=\int_{0}^{1}{\phi \left( f \right)p\left( f \right)}df$ (9.1-15)

$\overline{{{{{\phi }’}}^{2}}}=\int_{0}^{1}{{{\left[ \phi \left( f \right)-\bar{\phi } \right]}^{2}}p\left( f \right)}df$ (9.1-16)

可见，通过$f$的概率密度函数$P$($f$)及其基本公式就可以得到流场的各参数的时均值。

一般概率密度函数是未知的，确定方法有两种：一是直接建立以$P$($m$)

为因变量的微分方程，通过求解该输运方程，找出$P$$d$$f$的分布形式，这就是概率密度分布函数输运方程模型；另一种方法是根据对湍流脉动的设想，事先给出$P$$d$$f$的分布形式，可称之为简单概率密度函数模型或简化的$P$$d$$f$模型。常用的概率密度函数表示形式有城墙式分布、截尾正态分布概率分布函数、$\beta $函数等。

这样，用$k-\varepsilon $模型来描述湍流特性，燃烧按简单化学反应系统和快速反应假设进行简化，建立$\bar{f}$和$g$方程，利用假设的$P$($f$)，各标量的时均值都可求出。采用这种燃烧模型时，无需求解组分方程，而只求解$f$的均值和均方根方程，这样便大大提高了求解速度，这也是这种模型的特点和优点。

(2) 旋涡破碎模型

关于湍流预混燃烧反应速率的计算模型，最早由Spalding提出，即众所周知的旋涡破碎模型（Eddy Break-Up Model）。该模型的基本假设是：对于管内湍流均匀湍流火焰，在气流雷诺数很大、湍流强度也较大时，化学反应速率很快，因此混合气体中燃料的消耗率主要取决于气动力学，而与化学动力学关系不大。因此可以认为，在湍流燃烧区中充满了已燃和未燃的气团，化学反应在这两种气团交接面上进行，当湍流速度梯度增大，使气团进一步分裂，交接面也进一步增加，因此气团尺度减小快慢是决定混气反应速率的主要因素。

在发展过程中，旋涡破碎模型有多种表达式，但都应用了上述的基本思想。1976年，由Spalding提出的反应速率表达式：

${{R}_{fu}}=C{}_{EBU}{{\bar{m}}_{fu}}\rho {\varepsilon }/{k}\;$ (9.1-17)

其中${{C}_{EBU}}\approx 5$是可调的模型系数。该式是目前国内外常用的EBU模型表达式。

（3）涡团耗散模型

涡团耗散模型（Eddy-Dissipation Modal,简称EDM模型）是1976年由Magnussen等人提出的。其基本思想是：当气流涡团因耗散而变小时，分子之间碰撞机会增多，反应才容易进行并迅速完成，故化学反应速率在很大程度上受湍流的影响，而且反应速率还取决于涡团中包含燃料、氧化剂和产物中浓度最小的一个。该模型的特点是物理意义比较明确，反应速率取决与湍流脉动衰减速率，并能自动选择成分来控制速率，因此该模型既能用于预混火焰，也能用于扩散火焰。

${{R}_{fu}}=-\left( \bar{\rho }{\varepsilon }/{k}\; \right)\min \left[ A{{{\bar{m}}}_{fu}},A{{{\bar{m}}}_{ox}}/S,B{{{\bar{m}}}_{pr}}/\left( 1+S \right) \right]$ (9.1-18)

式中$A\approx 4,B\approx 0.5,S$是化学当量比。

（4）有限反应速率模型

在湍流燃烧模型研究中，为了使基本方程得到封闭，需要模拟平均化学反应速率。除上述各种湍流模型外，还可采用湍流燃烧关联矩模型或概率密度函数反应速率模型。

在湍流燃烧中若两个反应物为燃料合氧气，其瞬时化学反应速率可用Arrhenius公式来表示，见式（9.1-9）。为了求时均反应速率，令

${{m}_{fu}}={{\bar{m}}_{fu}}+{{{m}’}_{fu}}$，${{m}_{ox}}={{\bar{m}}_{ox}}+{{{m}’}_{ox}}$，$T=\bar{T}+{T}’$

先对非线性指数项，利用泰勒级数近似展开：

$\exp \left( -{E}/{RT}\; \right)=\exp \left[ -\frac{E}{R\left( \bar{T}+{T}’ \right)} \right]=\exp \left[ -\frac{E}{R\bar{T}}{{\left( 1+\frac{{{T}’}}{{\bar{T}}} \right)}^{-1}} \right]$

当$\frac{{{T}’}}{{\bar{T}}}\langle \langle 1$时，$\exp \left( -{E}/{RT}\; \right)\approx \exp \left( -\frac{E}{R\bar{T}} \right)\exp \left( \frac{E}{R{{{\bar{T}}}^{2}}}{T}’ \right)$

如果$\frac{E}{R\bar{T}}$不太大，即$\frac{E}{R\bar{T}}\frac{{{T}’}}{R\bar{T}}$仍相当小时，可以认为：$\exp \left( -{E}/{RT}\; \right)\approx \exp \left( -\frac{E}{R\bar{T}} \right)\left[ 1+\frac{E}{R\bar{T}}\frac{{{T}’}}{{\bar{T}}}+\frac{1}{2}{{\left( \frac{E}{R\bar{T}}\frac{{{T}’}}{{\bar{T}}} \right)}^{2}} \right]$

忽略密度脉动，对式（9.1-9）进行雷诺平均，同时忽略三阶和高阶相关项，便可得到时均反应速率：

$\begin{align}& {{{\bar{R}}}_{fu}}=\overline{{{A}_{0}}{{\rho }^{2}}{{m}_{ox}}{{m}_{fu}}\exp \left( -E/RT \right)} \\& ={{A}_{0}}{{{\bar{\rho }}}^{2}}\left( {{{\bar{m}}}_{ox}}+{{{{m}’}}_{ox}} \right)\left( {{{\bar{m}}}_{fu}}+{{{{m}’}}_{fu}} \right)\exp \left( -\frac{E}{R\bar{T}} \right)\left[ 1+\frac{E}{R\bar{T}}\frac{{{T}’}}{{\bar{T}}}+\frac{1}{2}{{\left( \frac{E}{R\bar{T}}\frac{{{T}’}}{{\bar{T}}} \right)}^{2}} \right] \\& ={{A}_{0}}{{{\bar{\rho }}}^{2}}{{{\bar{m}}}_{ox}}{{{\bar{m}}}_{fu}}\exp \left( -\frac{E}{R\bar{T}} \right)\left( 1+F \right) \\\end{align}$ （9.1-19）

式中

$F=\frac{\overline{{{{{m}’}}_{fu}}{{{{m}’}}_{ox}}}}{{{{\bar{m}}}_{fu}}{{{\bar{m}}}_{ox}}}+\frac{E}{R\bar{T}}\left( \frac{\overline{{T}'{{{{m}’}}_{fu}}}}{\bar{T}{{{\bar{m}}}_{fu}}}+\frac{{T}'{{{{m}’}}_{ox}}}{\bar{T}{{{\bar{m}}}_{ox}}} \right)+\frac{1}{2}{{\left( \frac{E}{R\bar{T}} \right)}^{2}}{{\left( \frac{{{\bar{T}}’}}{{\bar{T}}} \right)}^{2}}$

为了封闭式（9.1-19）应当求解$\overline{{{{{m}’}}_{fu}}{{{{m}’}}_{ox}}}$，$\overline{{T}'{{{{m}’}}_{ox}}}$，$\overline{{T}'{{{{m}’}}_{fu}}}$，$\overline{{{{{T}’}}^{2}}}$。这可以通过求解关于它们的二阶标量关联矩的微分方程来获得，这时，把这种有限反应速率模型又称为湍流燃烧关联矩模型，或平均反应速率输运方程模型。虽然该模型物理概念较清楚，但因求解方程增多，使计算机的存贮量和计算时均大大增加。为此文献【60】提出一种代数关联矩模型，该模型的基本思想是把对湍流瞬时反应速率时间平均时出现的相关项用关联的代数式来表示：

$\overline{{{{{m}’}}_{fu}}{{{{m}’}}_{ox}}}={{C}_{y}}\frac{{{k}^{3}}}{{{\varepsilon }^{2}}}\frac{\partial \overline{{{m}_{fu}}}}{\partial {{x}_{j}}}\frac{\partial \overline{{{m}_{ox}}}}{\partial {{x}_{j}}}$

$\overline{{T}'{{{{m}’}}_{fu}}}={{C}_{y1}}\frac{{{k}^{3}}}{{{\varepsilon }^{2}}}\frac{\partial \overline{T}}{\partial {{x}_{j}}}\frac{\partial \overline{{{m}_{fu}}}}{\partial {{x}_{j}}}$

$\overline{{T}'{{{{m}’}}_{ox}}}={{C}_{y2}}\frac{{{k}^{3}}}{{{\varepsilon }^{2}}}\frac{\partial \overline{T}}{\partial {{x}_{j}}}\frac{\partial \overline{{{m}_{ox}}}}{\partial {{x}_{j}}}$

$\overline{{{{{T}’}}^{2}}}={{C}_{T}}\frac{{{k}^{3}}}{{{\varepsilon }^{2}}}{{\left( \frac{\partial \bar{T}}{\partial {{x}_{j}}} \right)}^{2}}$

代数关联矩模型物理意义明确，它表示各二阶相关项及湍流尺度与相应的时均量梯度的乘积成正比，并且计算工作量相应较少，适合于工程应用。

在补燃室掺混反应流场的模拟中，采用的湍流燃烧模型主要有湍流扩散火焰模型、旋涡破碎模型、涡团耗散模型，有限反应速率模型使用得较少。

9.1.3 湍流两相燃烧模型

湍流两相燃烧，如液雾和煤粉的燃烧、含凝聚相粒子的贫氧燃气的燃烧等，广泛存在于燃气轮机燃烧室、航空发动机加力燃烧室、冲压发动机燃烧室、火箭发动机燃烧室中。对于有颗粒反应（或液雾蒸发）和气相反应的湍流两相流动和燃烧的数值模拟，大多数研究者采用欧拉－拉氏模型，即对颗粒相采用轨道模型。轨道模型虽然可以较好地模拟颗粒反应经历，免去颗粒相的伪扩散，但是要给出三维空间详细的颗粒速度场和浓度场，需要很大的计算量，而且难以处理颗粒进入和离开回流区的问题。少数研究者采用颗粒相的欧拉模型，实质上是单流体或无滑移模型。双流体模型的优点是易于给出三维空间颗粒的速度分布、浓度和温度分布，可以用统一的数值方法求解气粒两相方程组，处理颗粒进入和离开回流区没有困难。然而用双流体模型必须解决颗粒湍流模型和颗粒反应经历问题。

为了解决用双流体模型模拟颗粒反应经历的问题，有的学者提出了湍流两相燃烧的双流体－轨道模型（continuum-trajectory model, CT model of particle phase）和全双流体模型(pure two-fluid model, PTF model)。两者都用同样的气相欧拉方程组、颗粒相欧拉连续方程和动量方程，包括同样的两相湍流模型。PTF模型只用欧拉坐标的偏微分方程来模拟颗粒的各种经历效应——包括由于水分蒸发、热解挥发和焦炭燃烧引起的颗粒质量变化和由于反应放热和两相间对流传热引起的颗粒温度变化。而CT模型则用拉氏微分方程组模拟上述各种经历效应^【61】。

颗粒轨道模型是在拉氏坐标系下处理颗粒相，在欧拉坐标系下处理气相，并把颗粒看成与流体有滑移、沿轨道运动的分散群，忽略颗粒群间的相互扩散。在补燃室的湍流两相燃烧反应流场的模拟中，采用较多的是颗粒轨道模型，这种模型的控制方程如下。

（1）气相控制方程

a.连续方程

$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial {{x}_{i}}}\left( \rho {{u}_{i}} \right)={{S}_{m}}$ （9.1-20）

其中${{S}_{m}}$是离散相对气相的源项。

b.动量方程

$\frac{\partial }{\partial t}\left( \rho {{u}_{i}} \right)+\frac{\partial }{\partial {{x}_{j}}}\left( \rho {{u}_{i}}{{u}_{j}} \right)=-\frac{\partial p}{\partial {{x}_{i}}}+\frac{\partial {{\tau }_{ij}}}{\partial {{x}_{j}}}+{{F}_{i}}$ （9.1-21）

其中p是静压，${{F}_{i}}$是离散相对气相的作用力。${{\tau }_{ij}}$是因分子黏性作用而产生的作用在微元体表面上的黏性应力，由下式给出：

${{\tau }_{ij}}=\left[ \mu \left( \frac{\partial {{u}_{i}}}{\partial {{x}_{j}}}+\frac{\partial {{u}_{j}}}{\partial {{x}_{i}}} \right) \right]-\frac{2}{3}\mu \frac{\partial {{u}_{j}}}{\partial {{x}_{j}}}{{\delta }_{ij}}$ （9.1-22）

其中${{\delta }_{ij}}$为单位张量，$i=j$时，${{\delta }_{ij}}=1$，$i\ne j$时；${{\delta }_{ij}}=0$。

c.能量方程

$\frac{\partial }{\partial t}\left( \rho E \right)+\nabla \cdot \left[ \overset{\scriptscriptstyle\rightharpoonup}{v}\left( \rho E+p \right) \right]=\nabla \cdot \left[ {{k}_{eff}}\nabla T-\sum\limits_{j}{{{h}_{j}}}{{{\overset{\scriptscriptstyle\rightharpoonup}{J}}}_{j}} \right]+{{S}_{h}}$ （9.1-23）

其中：$E=h-\frac{p}{\rho }+\frac{{{v}^{2}}}{2}$

$h=\sum\limits_{j}{{{Y}_{j}}{{h}_{j}}}$，${{h}_{j}}=\int_{{{T}_{ref}}}^{T}{{{c}_{p,j}}dT}$

${{k}_{eff}}=k+{{k}_{t}}$

式（9.1-23）中右端第一项表示因导热引起的能量输运，第二项表示因组分扩散而引起的能量输运。能量源项${{S}_{h}}$包含气相化学反应带来的热量及离散相产生的热量。

d.组分方程

$\frac{\partial }{\partial t}\left( \rho {{Y}_{i}} \right)+\nabla \cdot \left( \rho \overset{\scriptscriptstyle\rightharpoonup}{v}{{Y}_{i}} \right)=-\nabla \cdot {{\overset{\scriptscriptstyle\rightharpoonup}{J}}_{i}}+{{R}_{i}}+{{S}_{i}}$ （9.1-24）

其中${{\overset{\scriptscriptstyle\rightharpoonup}{J}}_{i}}=-\left( \rho {{D}_{i,m}}+\frac{{{\mu }_{t}}}{S{{c}_{t}}} \right)\nabla {{Y}_{i}}$，$S{{c}_{t}}$是湍流施密特数，$S{{c}_{t}}=\frac{{{\mu }_{t}}}{\rho {{D}_{t}}}$。${{R}_{i}}$是化学反应的净产生速率，${{S}_{i}}$为离散相的额外产生速率。

（2）颗粒相基本方程

a.运动方程

在笛卡儿坐标系下，x方向颗粒受力平衡方程（颗粒惯性＝作用在颗粒上的各种力）为：

$\frac{\text{d}{{\text{u}}_{\text{P}}}}{\text{dt}}={{\text{F}}_{\text{D}}}\text{(u}-{{\text{u}}_{\text{P}}}\text{)}$ (9.1-25)

${{\text{F}}_{\text{D}}}\text{(u}-{{\text{u}}_{\text{P}}}\text{)}$为颗粒的单位质量拖曳力。

其中：${{\text{F}}_{\text{D}}}=\frac{\text{18 }\!\!\mu\!\!\text{ }}{{{\text{ }\!\!\rho\!\!\text{}}_{\text{P}}}\text{d}_{\text{P}}^{\text{2}}}\frac{{{\text{C}}_{\text{D}}}\text{Re}}{\text{24}}$，$\text{Re}=\frac{\text{ }\!\!\rho\!\!\text{ }{{\text{d}}_{\text{P}}}\text{ }\!\!|\!\!\text{ }{{\text{u}}_{\text{P}}}-\text{u }\!\!|\!\!\text{ }}{\text{ }\!\!\mu\!\!\text{ }}$

${{\text{C}}_{\text{D}}}=\frac{\text{24}}{\text{Re}}\text{ }\!\![\!\!\text{ 1}\text{.0}+\text{R}{{\text{e}}^{\text{2/3}}}\text{/6 }\!\!]\!\!\text{ Re}\le \text{100}0$

${{\text{C}}_{\text{D}}}=\text{0}\text{.44 Re}>\text{1000}$

b.轨道方程

在三个方向，各有一个受力平衡方程，通过积分这些微分方程可求得颗粒的速度

u_p、v_p、w_p。再由$\frac{dx}{dt}={{u}_{P}}$ ，$\frac{dy}{dt}={{v}_{P}}$ ， $\frac{dz}{dt}={{w}_{P}}$ ，对速度积分得到颗粒的轨道。

除上述运动方程和轨道方程外，还需要有描述颗粒相物理变化和化学变化过程的方程。例如对于油珠的蒸发及燃烧，还应有蒸发方程和油珠加热方程。对于含硼推进剂的补燃室反应流场，应有描述硼粒子点火及燃烧的控制方程。

硼粒子的燃烧过程极其特别，主要原因是硼的燃烧产物具有很高的沸点，而使其包裹在硼粒子的表面，这样使氧气不能到达硼粒子表面而继续燃烧。20世纪70到80年代，美国、前苏联、德国等国对单个硼粒子的点火燃烧模型进行了深入地试验研究，得到了大量的试验数据并提出了硼粒子的点火燃烧模型^{【62，13，14】}，文献【65】对硼粒子的点火燃烧模型进行了较好的综述。1999年，普林斯顿大学的W.ZHOU提出了硼粒子点火及燃烧的多相模型（Multi-phase model）^【66】，这种模型，考虑了粒子表面的化学反应，可以模拟硼在含氟或不含氟的推进剂燃烧环境中的点火与燃烧。在这些有代表性的模型中，King模型的物理意义明确，数学表达式简单易懂，一些物理化学参数易于获得，因而被较多地用于补燃室反应流场的模拟。

c.硼粒子点火模型

${{d}_{P}}={{d}_{B}}+2\delta $

$\frac{d{{d}_{B}}}{dt}=-\frac{2{{R}_{B}}{{M}_{B}}}{\pi d_{B}^{2}{{\rho }_{B}}}$

$\frac{d\delta }{dt}=-\frac{({{R}_{B}}/2-{{R}_{E}}-{{R}_{H}}){{M}_{B2O3}}}{\pi d_{B}^{2}{{\rho }_{B2O3}}}$

$\frac{d{{T}_{P}}}{dt}=\frac{{\dot{Q}}}{(d_{B}^{3}\pi /6){{c}_{p,}}_{BL}{{\rho }_{BL}}+\pi d_{B}^{2}\delta {{c}_{p,}}_{B2O3}{{\rho }_{B2O3}}} \text{ (}{{T}_{\text{P}}}>\text{2450K)}$

$\frac{d{{T}_{P}}}{dt}=\frac{{\dot{Q}}}{(d_{B}^{3}\pi /6){{c}_{p,}}_{B}{{\rho }_{B}}+\pi d_{B}^{2}\delta {{c}_{p,}}_{B2O3}{{\rho }_{B2O3}}}\text{ (}{{T}_{\text{P}}}<\text{2450K)}$

$\frac{df}{dt}=\frac{{\dot{Q}}}{(d_{B}^{3}\pi /6){{\rho }_{B}}\Delta {{H}_{M}}} \text{ (}{{T}_{\text{P}}}=\text{2450K)}$

${{R}_{B}}=f({{d}_{P}},{{T}_{P}},{{P}_{O2}},\delta )$

${{R}_{E}}=f({{d}_{P}},{{T}_{P}},P)$

${{R}_{H}}=f({{d}_{P}},{{T}_{P}},{{P}_{H2O}})$

其中dB、δ、f、TP分别是硼粒子直径、氧化层厚度、熔化分数、粒子温度。R_B、R_E、R_H分别是硼消耗速率、B₂O₃蒸发速率、B₂O₃与水的反应速率,其具体表达式及其余符号的含义见文献【62】。上述模型中有关理化性能参数的取值见表9.1-1。

表9.1-1 有关理化性能参数

图9.1-1 算例1的计算结果

图9.1-2 算例2的计算结果

图9.1-1是算例1（见表9.1-2）的d_P、δ、f、T_P的变化曲线。从左图中可看到，初温为1800K的硼粒子置于温度为2100K的环境中，温度逐渐上升，在6ms时达到熔点；同时，氧化层厚度先增厚，然后又由于蒸发和反应消耗，氧化层厚度逐渐减小。从右图可看到，硼粒子直径开始下降得较慢，当氧化层厚度减小到一定程度，粒子直径迅速减小。在6ms时，硼粒子开始熔化，6.7ms时完全熔化，硼粒子点燃并进入燃烧阶段。硼粒子的点火时间（从粒子置于环境中到氧化层厚度为零且硼粒子全部熔化的时间）为6.7ms。数值计算中，δ的值不能为零，一般当δ值小于初始值的20%后，便认为氧化层已经完全去除^【66】。

图9.1-2是算例2的δ、T_P曲线。可以看到，其它条件不变而降低环境温度后，硼粒子温度开始逐渐上升，但当上升至约2060K左右后，趋于平缓。而液化层厚度的变化也是逐渐趋于平缓，以至于到20ms时还有约70％的氧化层没有去除。可见，当环境温度太低时，硼粒子不能点火。

d. 硼粒子燃烧模型

硼粒子燃烧的总反应为^【67】：

2B(s)＋O₂(g) → B₂O₂(g)+110kcal/mol (3－26)

燃烧速率为：$\dot{m}=kd_{P}^{2}\pi {{P}_{O2}}$，其中$k=6.662\times {{10}^{-5}}kg/s\cdot {{m}^{2}}\cdot Pa$，为反应速率常数。

（3）离散相与气相相互作用源项表达式

通过上述离散相各方程求得在计算区域内颗粒尺寸、速度和温度后，再计算每个控制容积中颗粒质量、动量和能量的增减，求得颗粒源项。将每个起始尺寸组的颗粒各自源项叠加起来，就得出每个控制容积中离散相和气相之间相互作用的总源项。

a.质量源项

${{S}_{m}}=\sum\limits_{L=1}^{N}{{{\left( {{m}_{di}}-{{m}_{do}} \right)}_{L}}}$ (9.1-27)

式中${{m}_{di}}$和${{m}_{do}}$分别为控制容积进口和出口处颗粒质量，代表颗粒尺寸组总数。

b.动量源项

${{F}_{x}}=\sum\limits_{L=1}^{N}{{{\left( m{{u}_{do}}-m{{u}_{di}} \right)}_{L}}}$ (9.1-28)

${{F}_{y}}=\sum\limits_{L=1}^{N}{{{\left( m{{v}_{do}}-m{{v}_{di}} \right)}_{L}}}$ (9.1-29)

${{F}_{z}}=\sum\limits_{L=1}^{N}{{{\left( m{{w}_{do}}-m{{w}_{di}} \right)}_{L}}}$ (9.1-30)

式中${{m}_{di}}$和${{m}_{do}}$分别代表网格出口和进口粒子的$x$向动量，$m{{v}_{do}}$和$m{{v}_{di}}$以及$m{{w}_{do}}$和$m{{w}_{di}}$和分别为网格出口和进口粒子的$y$、$z$向动量。

c.总焓变化率

${{S}_{ph}}=\sum\limits_{L=1}^{N}{\left( {{m}_{di}}-{{m}_{do}} \right)\left( {{H}_{fu}}-L \right)-{{Q}_{oi}}{{c}_{Pp}}\left( {{T}_{p}}-{{T}_{po}} \right)}$ (9.1-31)

9.1.4 多相湍流反应流动的数值解法

一般来说，描述多相湍流反应流动的方程都是非线性偏微分方程组，不能用解析方法求解，必须用数值方法求解。为此首先要把微分方程离散化，可以用有限元或有限差分或有限分析法进行离散。在计算流体力学方面，主要用离散为差分方程的办法。在流动，传热和燃烧的数值模拟方面，比较普遍和流行的是SIMPLE（半隐式压力相关法）算法，关于这种算法的原理可参考文献【9】。本节只简要介绍多相湍流反应流动的数值解法。

求解颗粒轨道模型的方法是“单元内颗粒源法”（Particle Source in Cell Method），简称PSIC法。该法的基本思想是：在拉格朗日坐标下考察各组颗粒群沿各自轨道运动、质量损失及能量变化的过程，在欧拉坐标系下处理气相场，由颗粒群沿轨道中阻力、蒸发、挥发或燃烧产生颗粒群速度、温度、尺寸变化作为质量、动量及能量源加入气相场中成为气相与颗粒相之间相互耦合。国际上一些著名的模拟流动、传热和燃烧的商业软件，如FLUENT，STAR-CD等也都采用了这种算法。在PSIC算法中，对气相流场可以采用通常的SIMPLE系列算法（SIMPLE,SIMPLER,SIMPLEC等），即用压力－速度修正，逐线或逐面迭代的对三角矩阵计算和低松弛法求解。求解离散相和气相耦合流场的流程图如图3－3 。

图9.1-3 颗粒轨道模型求解过程流程图

9.1.5 应用举例

算例一：

由于固体推进剂中添加金属燃料，如铝、镁、硼等，补燃室内的燃气中含有大量的粒子相，是典型的两相或多相反应流场。燃烧室的构型不同，在解法上也带来许多特殊问题。限于篇幅，此处仅以含硼贫氧推进剂在补燃室中的掺混反应流场为例来说明数值模拟的结果及应用。

一种典型的非壅塞固体火箭冲压发动机的补燃室结构如图2所示。进气道与发动机轴向的夹角为60°，进气道周向夹角为120°。贫氧燃气从头部的三个喷口喷入补燃室内。

图3－4 典型补燃室结构示意图

从燃气发生器产生的一次燃气的成分极为复杂，它们在补燃室中与氧气的反应机理就更为复杂。为了能够近似模拟补燃室中的化学反应，可通过热力计算，选取其中的主要成分进行模拟。对于气相反应，可用涡耗散模型来确定反应速率，用颗粒轨道模型来处理离散相的燃烧过程，其中硼粒子的点火过程可采用如上所述的King模型。按照这种思路，文献【68】用FLUENT软件对某固体火箭冲压发动机补燃室的掺混反应流场进行了数值计算。

图9.1-5是温度分布，补燃室头部温度约2000K，后部的温度稍高，最高温度约2280K左右。

图9.1-6是硼粒子的点火位置图，当轨道为蓝色时，表示粒子没有点火，当轨道为红色时表示粒子已经点火，两种颜色交界处为硼粒子的点火位置。从图中可见，所有硼粒子在进气道之前已经完成了点火。

图9.1-7是硼粒子质量变化图。从中可见，在点火过程中，粒子质量几乎没有发生变化。在燃烧阶段，大部分粒子的质量都逐渐减小，直至质量为零，有少量粒子在喷出燃烧室时还没有完全燃烧。

图9.1-5 温度分布图9.1-6 硼粒子点火位置

图9.1-7 硼粒子质量变化

通过分析燃料流率的变化可得出燃烧效率沿长度方向的变化。其他参数相同，空气入射角度分别为30°、45°、60°、90°时的燃烧效率曲线如图6。

（a）补燃室燃烧效率（b）硼粒子燃烧效率

图9.1-8 空气入射角度对燃烧效率的影响

从图9.1-8（a）可见，在x=1.0m位置处，空气入射角度为45°和60°时的燃烧效率相等，达到86％；入射角度为30°时的燃烧效率较低，为80％；90°时的燃烧效率最差，仅有65％。图（b）所示的硼粒子燃烧效率比较曲线也符合这种规律。可见，空气入射角度为45°或60°时，补燃效率较高；入射角度大于60°时，补燃效率将下降。

不同直径粒子的补燃室燃烧效率、硼粒子燃烧效率分别见图9.1-9(a)、（b）。从图9.1-9（a）可以看到：在补燃室后部，当粒子的直径为2$\mu m$时，补燃效率达到85%，而直径10$\mu m$的仅有65%。由图（b）可见，粒子流遇到前进气道（x＝0.132m）的空气后，直径较小的粒子燃烧效率急剧增加，而直径较大的粒子增加缓慢。在流场的下游，较小粒子的燃烧效率增加明显，最后趋于平缓，而直径为5$\mu m$和10$\mu m$的粒子效率增加得一直较缓慢。在补燃室后部，直径2$\mu m$的粒子燃烧效率达到了75%，而直径10$\mu m$的仅有15%。

图9.1-9 硼粒子直径对燃烧效率的影响

从以上论述可见，对补燃室进行数值模拟，不仅可了解流场分布，而且可分析补燃效率的变化，这对补燃室的优化设计是十分有益的。

算例二：

两次进气是为了提高硼在补燃室内的燃烧效率而提出的进气方式，但是人们还没有清楚地认识两次进气对燃烧效率的影响机理，以及如何设计二次进气间距以达到最佳的燃烧效率。本文针对典型结构的整体式固体冲压发动机补燃室（图9.1-10），运用数值模拟方法研究了纯气相条件下二次进气间距对燃烧效率的影响，为进一步开展带凝聚相的数值模拟工作奠定基础。本文对二次进气间距${{d}_{12}}$分别为100、150、250、350mm共4种工况进行了数值模拟，主要计算结果见图9.1-11~3。在单次进气条件下采用相同的燃气喷嘴构形和相等的进气道面积，经过多次分析计算获得最高燃烧效率为81%。

图9.1-10 采用两次进气的固体火箭冲压发动机构型及计算网格

图9.1-11 不同二次进气间距下燃烧效率

从图9.1-12可以看出燃烧效率随进气间距的变化趋势：燃烧效率先随进气间距的增大而增大，在达到某一最大值后，燃烧效率将随进气间距的增大而减小，这表明存在一个最优的二次进气间距使燃烧效率最高。对于本文算例来说，该最优进气间距约为250mm，使补燃室内燃烧效率达到93.5%，与单次进气条件下的最高燃烧效率（81%）相比提高了12.5%。

(a) 无二次进气 (b) 进气间距mm

图9.1-12 流线图

图9.1-12显示的是不同进间气距条件下燃气（黑色线）和来流空气（灰色线）的流线，图9.1-13显示的是发动机和进气道轴线平面上的燃料分布，从中可以分析燃烧效率随进气间距变化的主要原因。本文采用三个喷口的燃气发生器，单次进气时补燃室内存在燃气的交汇区域（图9.1-12a），此处为燃气浓度较高的区域，若第二次进气在该区域切入（图9.1-12c），可增强空气与燃气的掺混程度，提高燃烧效率（图9.1-13c）。当两次进气间距过短时，第二次进气不会对补燃室中的掺混和燃烧产生显著影响（图9.1-12b），经过补燃后补燃室内仍存在燃气浓度较高区域（图9.1-13b），因此过短的进气间距对提高燃烧效率没有明显的作用（图9.1-13b）。如果第二次进气距冲压喷管入口太近，则没有足够的长度完成掺混燃烧，导致燃烧效率下降（图9.1-12d、图9.1-13d）。

图9.1-13 进气道对称面上甲烷质量分数等值线

数值模拟结果表明，不同的进气间距将促使二次进气在补燃室内呈现不同的流动特性，由此将对燃烧效率产生显著影响。通过对补燃室内流态的分析发现，二次进气在燃气交汇区域切入，可增强空气与燃气的掺混，有利于提高燃烧效率，即存在一个最优的二次进气间距使补燃室内燃烧效率达到最大。研究结果可以为两次进气结构的固体火箭冲压发动机设计提供参考。

9.2 进气道数值分析

目前，Navier-Stokes方程（简称N-S）是描述连续介质流动的最完整形式。它是以考虑流体黏性和热传导等因素，建立于质量、动量和能量守恒基础上的数学模型。由于非线性偏微分方程组数学理论还不完善，人们无法对N-S方程进行更为深入的数学分析，所以数值模拟的方法仍然是求解N-S方程的主要手段。但直接对全N-S方程进行数值模拟还存在着一些理论上和计算条件上的限制，因而简化N-S方程颇受人们青睐。

9.2.1 控制方程

在黏性流体力学中，描述流体运动规律的是N-S方程，它是连续方程、动量方程和能量方程的联立方程，在直角坐标系下的三维可压N-S方程的守恒形式为：

$\frac{\partial U}{\partial t}+\frac{\partial {{F}_{1}}(U)}{\partial x}+\frac{\partial {{F}_{\text{2}}}(U)}{\partial y}+\frac{\partial {{F}_{3}}(U)}{\partial z}=\frac{\partial {{G}_{1}}}{\partial x}+\frac{\partial {{G}_{2}}}{\partial y}+\frac{\partial {{G}_{3}}}{\partial z}$ 　　 (9.2-1)

$E=\rho (e+\frac{{{u}^{2}}+{{v}^{2}}+{{w}^{2}}}{2})$

在上述公式中，x_i（i＝1,2,3）对应于坐标x，y和z，u_i对应于u，v，w，分别为三个方向上的速度分量。ρ为密度，T为绝对温度，p为压强，μ为黏性系数，k为传热系数，τ_ij为黏性应力张量的分量；雷诺数Re＝ρUL/μ；普郎特数Pr＝c_pμ/k。为了使方程封闭，再增加一个状态方程。对于完全气体，状态方程可以写为：

_{$p=\rho T/(\gamma Ma_{\infty }^{2})$ (9.2-2)}

且有e＝c_VT，这里e为内能，c_V为质量定容比热容，γ为气体比热容比，对于完全气体，γ＝1.4，Ma_∞为来流马赫数。黏性系数μ由萨瑟兰（Sutherland）公式给出

$\frac{\mu }{{{\mu }_{\infty }}}\frac{{{T}^{{3}/{{}}\;2}}(1+{{T}_{s}}/{{T}_{\infty }})}{T+{{T}_{s}}/{{T}_{\infty }}}$ (9.2-3)

式中：T_s为参考温度，对于空气一般取T_s＝110.4K。

从上面的方程形式看Navier-Stokes方程加上气体状态方程是可以封闭的，但用于求解湍流时由于湍流中某一点处的各物理量的随机性质，直接求解瞬时的流动状况是不可能的，对湍流的各物理量进行时间平均（即通常所说的雷诺平均）是解决湍流问题的重要途径。通过对控制方程组进行时间平均处理，得到雷诺平均控制方程组，该方程组在三维笛卡尔坐标系中的微分形式的守恒形式如下：

（1）连续方程

$\frac{\partial \overline{\rho }}{\partial t}+\frac{\partial (\overline{\rho }\overline{{{\omega }_{i}}})}{\partial {{x}_{i}}}=0$ (9.2-4)

（2）动量方程

$\frac{\partial (\overline{\rho }\overline{{{\omega }_{i}}})}{\partial t}+\frac{\partial (\overline{\rho }\overline{{{\omega }_{i}}}\overline{{{\omega }_{j}}})}{\partial {{x}_{i}}}+\frac{\partial \overline{p}}{\partial {{x}_{i}}}-\frac{\partial (\overline{{{\tau }_{ij}}}-\overline{\rho }\overline{\omega _{i}^{‘}}\overline{\omega _{j}^{‘}})}{\partial {{x}_{i}}}=0$ (9.2-5)

（3）能量方程

$\frac{\partial e}{\partial t}+\frac{\partial (e\overline{{{\omega }_{i}}})}{\partial {{x}_{i}}}+\frac{\partial (\overline{p}\overline{{{\omega }_{i}}})}{\partial {{x}_{i}}}+\frac{\partial (\overline{{{\tau }_{ij}}}\overline{{{\omega }_{j}}}-\overline{{{q}_{i}}})}{\partial {{x}_{i}}}=0$ (9.2-6)

其中$e=\overline{\rho }{{c}_{V}}\overline{T}+\frac{1}{2}(\sum\limits_{i=1}^{3}{\overline{\omega }_{i}^{2}}+\sum\limits_{i=1}^{3}{\overline{\omega _{i}^{‘}\omega _{j}^{‘}}})$为总能量，包括内能、时均动能和湍流动能；

${{\overline{\tau }}_{ij}}=\mu (\frac{\partial \overline{{{\omega }_{i}}}}{\partial {{x}_{j}}}+\frac{\partial \overline{{{\omega }_{j}}}}{\partial {{x}_{i}}})-\frac{2}{3}{{\delta }_{ij}}\mu \frac{\partial \overline{{{\omega }_{l}}}}{\partial {{x}_{l}}}-\overline{\rho }\overline{{{\omega }_{i}}{{\omega }_{j}}}$ （i=1，2，3）

$\overline{{{q}_{i}}}=K{{\frac{\partial \overline{T}}{\partial {{x}_{i}}}}^{{}}}^{{}}(K={{c}_{P}}(\frac{\mu }{\Pr }+\frac{{{\mu }_{t}}}{{{\Pr }_{t}}})$，$\mu$和${{\mu }_{t}}$分别时分子黏度和湍流黏度，$Pr$和${{\Pr }_{t}}$分别为层流和湍流Prandtl数）；i＝1，2，3分别表示笛卡儿坐标系的三个坐标分量x，y，z。

9.2.2 湍流模型

雷诺平均方程组中除了时均项外，还存在一些脉动参数的相关量，这增加了方程的未知数数量，使得方程不再封闭。根据湍流的流动规律以寻求脉动量与时均量之间的关系，从而使方程封闭可解，就是建立湍流模型的过程。

目前两方程湍流模型在工程中使用最为广泛，最基本的两方程模型是标准k-ε模型，它是通过求解两个独立的输运方程，来分别确定湍流的速度尺度和长度尺度。由Boussinesgue涡黏模型，雷诺应力与应变率存在线性关系，则有效黏度近似湍流应力：

$\overline{{{u}_{i}}{{u}_{j}}}=\frac{2}{3}{{\delta }_{ij}}k-\frac{{{\mu }_{t}}}{\rho }\left( \frac{\partial {{U}_{i}}}{\partial {{x}_{j}}}-\frac{\partial {{U}_{j}}}{\partial {{x}_{i}}} \right)$ (9.2-7)

式中，湍流黏度${{\mu }_{t}}=\rho {{C}_{\mu }}{{k}^{2}}/\varepsilon $，${{U}_{i}}$为时均速度，${{\delta }_{ij}}$为“Kronecker delta”符号（当$i=j$时，${{\delta }_{ij}}=1$，当$i\ne j$$时，${{\delta }_{ij}}=0$）；其中湍流动能k以及它的耗散率可由Jones和Launder提出的微分方程求解：

$\rho \frac{Dk}{Dt}=\frac{\partial }{\partial {{x}_{j}}}\left[ \left( \mu +\frac{{{\mu }_{t}}}{{{\sigma }_{k}}} \right)\frac{\partial k}{\partial {{x}_{j}}} \right]+{{G}_{k}}-\rho \varepsilon -{{Y}_{M}}$ (9.2-8)

$\rho \frac{D\varepsilon }{Dt}=\frac{\partial }{\partial {{x}_{j}}}\left[ \left( \mu +\frac{{{\mu }_{t}}}{{{\sigma }_{\varepsilon }}} \right)\frac{\partial \varepsilon }{\partial {{x}_{j}}} \right]+{{C}_{\varepsilon 1}}\frac{\varepsilon }{k}{{G}_{k}}-{{C}_{\varepsilon 2}}\rho \frac{{{\varepsilon }^{2}}}{k}$ (9.2-9)

式中，${{\sigma }_{\varphi }}$为普朗特－施密特数（下标$\phi $代表k和$\varepsilon $）；湍动能生成率${{G}_{k}}=-\rho \overline{{{u}_{i}}{{u}_{j}}}\partial {{U}_{i}}/\partial {{x}_{j}}$，它反映了通过大尺度的湍流涡和平均应变之间相互作用引起的从平均过程到脉动过程的能量传递。Y_M代表在可压湍流里波动扩张对整个耗散率的贡献，由Sarkar^【22】提出的建议：

${{Y}_{M}}=2\rho \varepsilon M_{t}^{2}$ (9.2-10)

式中，M_t是湍流马赫数，定义为：${{M}_{t}}=\sqrt{\frac{k}{{{a}^{2}}}}$，这里α(=$\sqrt{\gamma RT}$)是声速。方程中的常数C_μ＝0.09, ${{\sigma }_{k}}$＝1，${{C}_{\varepsilon 1}}$＝1.44，${{C}_{\varepsilon 2}}$＝1.92，${{\sigma }_{\varepsilon }}$＝1.3。

标准$k-\varepsilon $模型是针对湍流发展非常充分的湍流流动来建立的，也就是说，它是一种针对高Re数的湍流计算模型，而当Re比较低时，例如，在近壁区内的流动，湍流发展并不充分，湍流的脉动影响可能不如分子黏性的影响大，在更贴壁面的底层内，流动可能处于层流状态。因此，对Re数较低的流动使用标准的$k-\varepsilon $模型进行计算就会出现问题。常用的解决方法有两种，一种是采用壁面函数法，另一种是采用低Re数的$k-\varepsilon $模型，本文采用壁面函数法来解决这个问题。

9.2.3 湍流的壁面处理

壁面是造成漩涡状态和湍流的来源，在可压流中，亚音流域的温度分布与近壁面区域的温度分布有很大的不同，由于壁面附近速度梯度很大和湍流动能的产生，在近壁区的较外区域湍流速度增强，采用不同的壁面处理模式对数值模拟精度影响很大。在壁面边界层内，温度、压力等各变量的梯度很大，动量等的输运剧烈，因此，正确地处理近壁区域的流动非常重要。

由于考虑到气体的黏性，因此，要求解边界层流域，包括划分网格节点和在次节点上求解有限差分方程。对于湍流流动，靠近壁面层区的速度梯度很大，需要很密的节点网格，在其余的区域可用较大的网格间距。

在近壁面区，可大致分为三个区域：在最内层，流动近似为层流，分子黏性在动量、热量、质量传递中起主要作用，在最外层，成为完全湍流层，湍流起主要作用，在两层之间，是过渡层，分子黏性和湍流作用同等重要。

$k-\varepsilon $ 两方程模型为高Re数模型，适用于离开壁面一定距离的湍流区域，在高Re数区域，分子黏性系数相对湍流黏性系数可略而不计。在与壁面相邻的黏性层中，湍流Re数很低，必须考虑分子黏性的影响，对于壁面附近的区域需特别处理。

壁面函数法在黏性底层内不布置任何节点，将第一个与壁面相邻的节点布置在旺盛湍流区内，然后对该节点与壁面之间的区域的当量黏性系数μ_t及K_p，ε_p做出规定。不求解层流底层及混合区，通过半经验公式获得层流底层与湍流之间区域信息。

对于非平衡的壁面函数法，近壁层流体的速度场控制为：

$\frac{\text{\tilde{U}C}_{\mu }^{1/4}{{k}^{1/2}}}{{{\tau }_{w}}/\rho }=\frac{1}{\kappa }\ln \left( E\frac{\rho \text{C}_{\mu }^{1/4}{{k}^{1/2}}y}{\mu } \right)$ （9.2-11)

其中：$\tilde{U}=U-\frac{1}{2}\frac{dp}{dx}\left[ \frac{{{y}_{v}}}{\rho \kappa \sqrt{k}}\ln \left( \frac{y}{{{y}_{v}}} \right)+\frac{y-{{y}_{v}}}{\rho \kappa \sqrt{k}}+\frac{y_{v}^{2}}{\mu } \right]$

${{y}_{v}}=\frac{\mu y_{v}^{*}}{\rho C_{n}^{1/4}k_{p}^{1/2}}$

$y_{v}^{*}=11.225$

κ为卡门常数，κ＝0.42，E为经验常数，E＝9.81，U为流体的平均速度，为p点流体的湍流动能，y为与壁面间的距离，μ为流体的动力黏度。

能量和动量的转换过程是非常复杂的，通常是应用雷诺比拟准则的相似关系得到它们之间的转换关系，这个相似关系为：对于近壁层的平均速度、温度，在占优的区域内，认为温度是线性分布的，在湍流占优的区域内，认为温度是对数分布的。

近壁层流体的温度分布规律为：

其中：$p=9.24\left[ {{\left( \frac{\Pr }{{{\Pr }_{t}}} \right)}^{\frac{3}{4}}}-1 \right](1+0.28{{e}^{-0.007\Pr /{{\Pr }_{t}}}})$

k_p为流体的导热系数，ρ为流体密度，c_P为流体的等压比热容容，$\dot{q}$为壁面上的热流密度。T_p是与壁面相邻的控制体积的节点P处的温度，T_W为壁面的温度，Pr为分子普朗特数，Pr_t为湍动分子普朗特数,A为冯.D常数，A＝26，$y_{r}^{*}$为线性分布和对数分布交界处的$y_{{}}^{*}$值，U_c是在$y_{{}}^{*}=y_{r}^{*}$处的平均速度值。

在近壁面区域湍流动能k和湍流耗散率ε用下面的方程计算：

G_k为由于速度梯度的湍流动能生成项。

对于分离流场，运用非平衡的壁面函数法，能部分解决压力梯度大和分离流动问题，改善对壁面和热传递的计算结果。

9.2.4 离散方法

对于在求解域内所建立的偏微分方程，必须通过数值方法把计算域内有限数量位置（网格节点或网格中心点）上的因变量值当做基本未知量来处理，从而建立一组关于这些未知量的代数方程组，然后通过求解代数方程组来得到这些节点值，而计算域内其他位置上的值则根据节点位置上的值来确定。根据应变量在节点之间的分布假设以及推导离散化方程的方法，对控制方程的离散主要分为：有限差分法、有限元法和有限体积法，其中有限体积法是目前CFD领域使用最为广泛的离散化方法。

　本文选用有限体积法，该方法的基本思路是：在每一个控制体内积分控制方程，从而产生基于控制体的每一个变量都守恒的离散方程。将通用变量φ的定常状态下守恒方程写成对于控制体积V的积分形式：

$\oint{\rho \phi \overrightarrow{\nu }}\cdot d\overrightarrow{A}=\oint{\Gamma \nabla \phi }\cdot d\overrightarrow{A}+\int_{V}{S\cdot dV}$ （9.2-14）

　其中ρ为密度，$\overrightarrow{\nu }$为速度矢量，$\overrightarrow{A}$为曲面面积矢量，$\Gamma $为广义扩散系数，S为广义源项。式中各项依次为对流项、扩散项和源项。

　本文中控制方程的扩散项采用具有二阶精度的中心差分格式离散，对流项采用二阶迎风格式离散，黏性项采用中心差分格式离散，并在此基础上构建有限体积计算框架。当使用耦合求解方法时，流动方程选用二阶迎风格式离散，其他方程选用一阶迎风格式离散。

　一阶迎风格式假定描述单元内变量的单元中心变量就是整个单元内各个变量的值，而且单元表面的变量值等于单元内的变量值。因此，当选择一阶迎风格式时，表面值${{\phi }_{f}}$等于迎风单元中心值。

　二阶迎风格式使用多维线性重建方法来计算单元表面处的值。在这种方法中，通过单元中心解在单元中心处的泰勒展开来得到单元表面的二阶精度。因此，当使用二阶格式时，用下面的方程来计算表面值${{\phi }_{f}}$：

${{\phi }_{f}}=\phi +\nabla \phi \cdot \Delta S$ 　　　（9.2-15）

　其中$\phi $和$\nabla \phi $分别是单元中心值和迎风单元的梯度值，$\Delta S$是从迎风单元中心到表面中心的位移矢量。在这种情况下需要确定每个单元内的梯度$\nabla \phi $。我们使用散度定理来计算这个梯度，其离散格式如下：

$\nabla \phi =\frac{1}{V}\sum\limits_{f}^{{{N}_{faces}}}{{{\overline{\phi }}_{f}}A}$ 　（9.2-16）

在这里，表面处的值$\overline{{{\phi }_{f}}}$由临近表面的两个单元的$\phi $的平均值来计算。最后，应限制梯度$\nabla \phi $以保证不会引进新的最大值和最小值。

9.2.5 边界条件的处理

边界条件的选取不仅会影响计算精度，而且也关系到求解过程的稳定和收敛过程。计算中涉及的边界条件主要有：

a. 固体壁边界条件

在壁面上速度采用无滑移条件：

u＝0 　　　　　　　　　（9.2-17）

v＝0

定常情况下，假设壁面法线方向上的压力梯度为零，即

$\frac{\partial p}{\partial n}=0$ （9.2-18）

假设固体壁绝热，温度梯度为零，即

$\frac{\partial T}{\partial n}=0$ （9.2-19）

对于$k-\varepsilon $方程，固体壁边界条件为

k＝0 　　　　　　　　　（9.2-20)

ε＝0

b. 出口边界条件

进气道出口边界为补燃室入口截面，气流速度降为亚声速，选用压力出口边界条件，进气道上壁面外延形成一个超声速气流出口，也选用压力出口边界条件。

对于超声速出口，不需给定任何边界条件，全部气流参数由内部流场二阶外推算出。

对于亚声速出口，由于受外界影响，需要给定反压；该边界条件可以处理出口回流问题，如果出口有回流存在，则回流条件需要给定：回流总温（能量方程）、湍流参数（湍流模型）等，合理的给定出口回流条件，有利于解决有回流出口问题的收敛难问题。回流流动方向与出口边界垂直。

c. 对称面边界条件

图中的滑移固体壁采用对称面边界条件。此时，既无质量的交换，也无热量等其他物理量的交换，因此对称边界的法向速度为零，即

v＝0 　　　　　　　（9.2-21）

在对称面上，所有物理量在其垂直方向上的梯度为零，即对于任意变量Φ有

$\frac{\partial \Phi }{\partial y}=0$ 　　　（9.2-22）

计算中不需要给定任何参数，只需要确定合理的位置。

d. 压力远场边界条件

采用法线方向上的黎曼不变量来建立无反射条件，用于模拟无穷远处的自由可压流动，需要给定静压、温度、马赫数、流动方向和湍流强度；为了满足无穷远边界条件，需要把边界放在我们关心区域足够远的地方。

9.2.6 算例分析

(1) 计算说明

本书算里采用FLUENT计算软件，计算中采用耦合隐式求解器，耦合求解器是同时求解连续方程、动量方程、能量方程及组分输运方程的耦合方程组，适用于高速可压缩流动，藕合隐式求解能量和动量方程，可较快的得到收敛解，湍流模型选用标准湍流$k-\varepsilon $模型，湍流模型的经验系数见9.2.2节所示。计算的收敛准则为：连续方程、动量方程、能量方程及k、ε方程的残差下降3个数量级，且进气道出口流量稳定。

在计算收敛后为了进一步改进网格中不合理的网格分布，通过改进网格使其更适合于流动计算，可以在先前求解的基础上，以速度梯度为基点来改善网格质量，然后继续进行计算，以期得到更准确的模拟结果。

(2) 物理特性

对于可压缩流体，由理想气体定律得

$\rho {}^{\left( {{p}_{OP}}+p \right)}/{}_{\frac{R}{M\omega }T}$ （9.2-23）

定压比热容容取定值，即${{c}_{P}}=1006.43\text{J/}\left( \text{kg}\centerdot \text{K} \right)$，热传导率也取为定值，即为$\lambda =0.0242$W/(m∙K)，黏性大小由Sutherland定律得到。

Sutherland定律是由萨瑟兰在1893年依据分子运动论，理想化分子间力得来，关系式可以是两系数或者三系数：

两系数关系式：

$\mu \frac{{{C}_{1}}{{T}^{1.5}}}{T+{{C}_{2}}}$，式中μ是黏度，单位为kg/m-s，T为静温，单位为K，C₁和C₂是系数。对于气体，在中等温度和压力下C₁ = 1.458×10^-6 kg/(m∙s∙K^1/2)，C ₂ = 110.4 K。

三系数关系式：
$\mu ={{\mu }_{0}}{{\left( \frac{T}{{{T}_{0}}} \right)}^{1.5}}\frac{{{T}_{0}}+S}{T+S}$ 　　　（9.2-24）

其中是$\mu $黏性，单位为kg/m∙s，T是静温，单位为K，${{\mu }_{0}}$是参考黏度值，单位为kg/m∙s, T₀是参考温度，单位为K，S是有效的温度，单位是K，被称为Sutherland常数，它是气体的特征值。对于适当的温度和压力：${{\mu }_{0}}$=1.7894×10^-5 kg/m∙s，T₀ = 273.11 K，S = 110.56 K。本文中选择三系数关系式。

(3) 进气道二维数值模拟

计算区域的确定

合理的选择计算区域可以更真实的模拟实际情况，有效的降低计算误差。考虑到进气道唇口区域在激波不脱体时会出现斜激波，在激波脱体时会出现弓形激波，对进气道特性参数会产生影响，因此作为进气道性能研究的一部分包含在计算域中；进气道楔面之前设定为气流入口，为了避免入口参数对附面层和斜激波的干扰，可将入口边界前移，增加一段滑移固体壁，具体的计算区域如图9.2-1所示。

网格划分及计算结果

网格划分是数值模拟好坏的关键环节之一，不合理的流场网格会造成数值解不精确、不真实、不易收敛甚至会引起计算的不稳定。本进气道采用结构化网格，设置附面层，在唇口上方，进气道入口部分，由于可能产生激波，物理量变化剧烈，采用局部加密；同时，为了能在进气道出口部分自后向前推出正激波，也采用局部加密，计算网格见图9.2-l。在进气道数值模拟过程中，通过改变边界条件，可以模拟进气道不同飞行高度、不同飞行马赫数和不同补燃室不同工况下进气道的流动情况。采用以上的计算方法得出二维流场计算结果见图9.2-2～9.2-5所示。

(4) 进气道与弹体一体化数值模拟

计算区域的确定

导弹与动力装置的整体化布局，要求进气道设计时应综合考虑导弹总体对进气道的气动影响。采用一体化流场数值模拟可较为准确地模拟在绕弹体非均匀来流条件下，弹体流场的变化情况和超声速进气道内激波、附面层的变化情况。混压式超声速进气道与弹体一体化流场数值模拟具有很多独特的特点和复杂性，需同时计算不同工况下的内外流场，计算域拓朴结构复杂。进气道二维型面确定后，根据设计给定的入口宽高比，三维型面就已经确定，需要注意的是，进气道入口处的上下唇口和两个侧壁均采用尖角过渡以减小阻力，防止低马赫数飞行时在唇口产生弓形激波。由于双下腹部进气道布局时两个进气道沿中心平面对称，作攻角飞行时两个进气道之间的相互影响比较小，从减少计算节点数量、节省计算时间考虑，只取一个进气道和一半的弹体流场进行分析是可取的，进气道与弹体一体化模型如图9.2-6所示。

网格划分

目前生成的网格主要有两类，即结构网格和非结构网格。结构网格的节点以阵列形式排列，与计算机语言自然匹配，便于矩阵验算与操作；非结构网格对网格节点没有结构性的限制，其节点与单元的分布具有任意性，因此具有优越的几何灵活性，尤其适用于具有复杂边界的流场计算问题，但它相对于结构网格需要更多的计算机内存及CPU时间。

为获得高质量的计算网格，进气道内部采用结构化网格，设置附面层，在入口和出口部分采用局部加密，保证数值模拟过程中由后向前逐渐推出正激波，弹体流场及隔道部分采用非结构化网格。另外，在计算过程中，还采用了自适应网格技术，由于激波前后气流的压力、速度等物理量发生突变，而自适应网格法正是从物理量的梯度出发构成网格，在物理量梯度大的地区，对网格加密，从而提高捕捉激波的分辨率，达到提高计算精度的要求。具体网格如图9.2-6所示，计算网格总数37万左右，唇口及进气道壁面网格分布见图9.2-7所示。在对进气道与弹体一体化数值模拟过程中，通过改变边界条件，不但可以模拟二维数值模拟所能模拟的情况，还可以模拟不同攻角和侧滑情况下弹体及进气道流场的变化情况。一体化流场数值模拟结果见图9.2-8～9.2-10。