航天器发射技术：射流计算方法

作为流体力学的一个分支，与流体力学一样，射流的研究方法主要有3种，即试验研究、理论分析和数值模拟，有时人们也将理论分析和数值模拟统称为理论分析或理论研究，这里理论分析指单纯使用流体力学基本理论研究射流的方法。

试验是自然科学的基础，理论如果没有实验的证明是没有意义的。力学是以试验为基础的科学，流体力学中绝大多数重要的概念和原理都源于试验。但是，试验往往受到模型尺寸、流动扰动、人身安全和测量精度的限制，有时可能很难通过试验方法得到结果。此外，试验还会遇到经费、人力和物力的巨大耗费及周期长等许多困难。

理论分析方法的优点在于所得结果具有普遍性，各种影响因素清晰可见，是指导试验研究和验证新的数值计算方法的理论基础。但是，它往往要求对计算对象逐项抽象和简化，才有可能得出理论解。对于自然界中普遍存在的非线性情况，只有极少数流动能给出解析结果。

作为数值模拟工具的计算流体力学（CFD）方法有效地克服了前面两种方法的弱点。CFD的长处是适应性强、应用面广，其不受物理模型和试验模型的限制，计算周期短，经费投入少，灵活性高，很容易模拟那些试验中只能接近而无法实现的理想条件，并得到满足工程需要的数值解。当然，CFD也存在一定的局限性，主要体现在对人员经验与技巧的依赖性和对计算机计算能力的高度依赖性上。但随着计算机计算能力的飞速发展，CFD的计算精度和计算效率已经有了显著的提高。

半个多世纪以来，特别是第二次世界大战末期以来，在火箭与导弹、航空与航天以及火炮武器等技术领域，经常遇到燃气射流问题，因此它在各种具体条件下的流场计算问题早为有关专业人员所关心。但由于它在很多情况下的流动结构极其复杂而且多变，给理论分析与计算造成了很大的困难与麻烦。因此，在燃气射流动力学发展早期（暂称它为第一阶段），除一些极其简单或被高度简化的射流可找到其理论计算方法外，一般以试验为主预估射流流场的参数及其对流经物体与周围环境的作用。

随着航空航天科技的发展，燃气射流动力学问题已到了急需更好地解决的时候，迫使人们不得不在某些限定条件下经过大量简化，从而提出一些较为系统的计算方法，但这往往会使计算精度大为降低。为弥补此缺陷，工程上往往尚需辅以较多的试验以对计算结果进行适当的修正。这种办法一般称为燃气射流的工程计算方法。这一时期称为燃气射流动力学发展的第二阶段。这个阶段的特点是计算与试验相结合，试验的目的是修正工程计算的算法和计算结果。

计算机和计算流体力学获得高度发展以来，燃气射流的某些问题可以用数值方法来计算求解，这时求解问题的范围和精度都有很大提高。这一时期称为燃气射流动力学发展的第三阶段。在这个阶段，虽然也需做一些试验，但其目的主要是验证计算结果。目前，该阶段仍在发展之中。

本小节简单介绍射流的工程计算方法和数值模拟方法。

1.工程计算方法

尽管工程计算方法存在精度不高的缺点，而且一般总需辅以试验，但在某些尚不具备数值计算条件和尚不能通过试验解决问题的条件下，要想取得相关工程上的某些概略数据，以适应工程设计之需，燃气射流的工程计算方法仍不失为一种有效的、可资利用的方法和手段。

1）轴对称亚声速等温自由射流的工程计算方法

气体向无限空间喷射称为自由射流；形象地看，它是被淹没在无限空间的介质之中的，所以又称它为淹没射流。严格地说，射流所射向的空间中的介质应与射流介质相同，才能称为自由射流。但一般工程计算中，只要喷向无限空间，即便介质不同，也视为自由射流。这种射流一般都是呈现为紊流（湍流）流动。

（1）亚声速自由射流的流动结构。

如图3.34所示，气体以均匀速度u0自喷口射出，由于它具有黏性和紊流横向脉动，故它在流动过程中不断地将周围的介质粘连和裹挟走一部分，亦即射流介质与周围介质二者不断发生质量和动量交换，结果是射流的质量流量和横截面积沿x方向不断增加，形成图3.34所示的锥形流场。一般认为喷口速度呈均匀分布，记为u0。由于周周介质的不断混入，均匀速度场呈现为图示的AO′D三角形，称为核心区。这时O′点所在截面BO′E便称为过渡截面或转捩截面。而自喷口截面至该截面的一段称为初始段，自该截面往后的一段称为主体段。图中ABC和DEF线均称为射流外边界或外边界线。而AO′和DO′线则称为内边界或内边界线。CA、FD二线向后延长在轴线上的交点O称为射流极点。∠AOD的一半称为极角，记为α。整个流场除去核心区的部分均称为边界层。显然，核心区内轴心线上以及全区内的速度均为u0，而主体段轴心线上的速度则沿x方向（原点为O）不断下降，且主体段完全为边界层所占据。

图3.34　亚声速自由射流的流动结构

（2）流动特性参数计算。

①紊流系数a。

紊流系数a是射流计算的一个关键参数，是具体表示射流流动结构的特征系数。它一般由试验确定，其值的大小与出口截面上的紊流强度有关，强度越大，则a值也越大。这说明射流与周围介质的混合能力强。其结果是使射流扩散角即极角α增大，因此被带动的周围介质增多，射流速度沿x轴的下降加快。此外，紊流系数a还与出口截面上速度分布的均匀性有关，对均匀分布射流，a＝0.066，若不太均匀，当＝1.25（为平均速度）时，则a＝0.076。

在工程计算中，a的近似取值如下：

对于轴对称收缩喷管，a＝0.066～0.071，紊流强度和分布不均匀度小者取小值，大者取大值；对于圆柱形喷管，a＝0.076～0.08，取值原则同上。

a对确定射流的几何参数起关键作用。参照图3.34，此处所谓的几何参数指射流的外边界走向即外边界方程，射流极点深度h0，射流初始段长度s0及射流内边界收缩角θ。下面逐步找到它们的求解公式。

②射流外边界方程。

根据试验和理论得知，射流的内、外边界都可以认为是直线。从图3.34可看出，BO′为过渡面上边界层的厚度（对轴对称射流而言，一般指圆截面半径R），它与从O点计起的x成正比，亦即R＝Kx，此处K为外边界的斜率，在射流理论中，它是一个试验系数，对轴对称射流而言，K＝3.4a，而a是紊流系数，已如上述。

由图3.34可得

式（3.2.59）或式（3.2.60）即射流的外边界方程。在工程计算中，通过极点O画极角为α的斜直线即得射流外边界。另外，外边界方程也可表示如下：

若用r0进行无因次化，又有

式（3.2.62）表明射流无因次半径与无因次距离成正比。

③射流极点深度h0。

由tanα＝r0／h0得

④速度场的求解。

自模性函数的半经验表达式（以速度场的自模性函数为例）为

令＝η，则式（3.2.64）简化为

式（3.2.65）对初始段和主体段都适用，但在不同段上式中符号的意义有所不同。对初始段，参见图3.34，式中：

r：边界层截面上任意点至内边界的距离；

R：同截面上的边界层厚度；

u：截面上边界层中r点的速度；

um：核心区速度，即u0。

对主体段，参见图3.34，式中：

r：截面上任意点至射流轴心线的距离；

R：同截面上的射流半径（半宽度）；

u：r点的速度；

um：同截面轴心线上的速度。

从式（3.2.64）或式（3.2.65）可看出，r／R或η的值从内边界或轴心线到射流外边界的变化范围为0～1，而u／um相应的变化范围为1～0。利用式（3.2.65）可直接求出初始段上核心区以外流场各点的无因次速度。在此，um＝u0＝常数，已知u0即可求得速度的真值。同理，可求主体段点的速度，但此段上的um是随x的增加而减小的，需要单独求解。

实际上，无因次速度的等速线就是一些由A或D点（对初始段而言）和由极点O（对主体段而言）引出的辐射线。因为无因次速度取决于η值，而每一条射线上的η值都相等。

⑤轴心速度um。

前已述及，当已知无因次速度值欲求流场的各点速度时，需知道um值。现根据动量守恒原理求解um的表达式。首先，亚声速射流流场各点的压强相等且等于周围介质的压强。这时，轴对称射流的动量守恒方程可这样表述，即射流出口的动量等于射流任意横截面（为简便起见，使用主体段上的横截面）上的动量，其表达式（参见图3.34）为

以除上式两端，得

利用式（3.2.65）可将上式右端的积分写为

由于η是无因次量，积分上式可得

式（3.2.70）说明，轴心速度um与无因次距离成反比例变化。

⑥初始段核心区长度s0及核心区边界（内边界）收缩角θ。

如图3.34所示，核心区长度s0可由式（3.2.69）或式（3.2.70）解出。因核心区端点O′处的速度仍为u0，而主体段轴心速度um在此点也恰好等于u0，所以a＝0.96，但＝0.294／a，所以有

而初始段内边界即核心区边界收缩角θ为

⑦主体段截面流量Q。

以射流出口流量为基准，取无因次流量进行计算，这时有

代入式（3.2.74），则无因次积分项的值为0.098 5，同时将式（3.2.69）所表达的和式（3.2.62）所表达的一并代入式（3.2.74），则得

(www.xing528.com)

式（3.2.76）说明，射流的截面流量是与无因次距离成正比例变化的。这是因为主射流不断地从周围环境中吸入一部分介质。

2）射流对物体作用力的近似确定

火箭、导弹在采用自力发射时，不可避免地要出现燃气射流对发射装置和发射载体产生冲击载荷的情况，而在多数情况下，为了减小这种冲击载荷的危害而迎着射流加装导流器。这就引出了求解导流器所受作用力的问题。

（1）物体在正面迎面冲击下气动载荷确定的一般原理。

不管物体的形状多么复杂，只要物体相对气流方向的正面迎风面积并结合气流动压头在整个面积上的分布规律，那么，该物体垂直于气流方向的总气动载荷就可按下述办法计算出来。

轴对称燃气射流的动压头可以表示成x和r的函数，即动压头是沿射流轴线距离和环绕轴线的半径的函数。当x为某一确定值时，动压头仅是环绕轴线的半径的函数。如图3.35所示，设F是置于射流中的距喷口为x的某平板的正面迎风面积；O点是射流轴线的投影；r是从轴线引出的半径。由于动压头随r的变化规律已知，气流流动对平板的作用力随r的变化规律是已知的。设它们的关系式为ρu2＝q（r），则平板所受的总气动载荷为

图3.35　射流冲击平板时气动载荷的确定

若在射流中放置的是图3.36所示的锥形导流器，由于气流动压头是（x，r）的函数，当导流器表面上的动压头自x0（这时x0可视为常数）往下仍能化为r的函数时，则锥体所受的总气动载荷仍可计算出来。此时dF仍是正面迎风微面积。若动压头化不成r的完整的连续函数，但可分成几个连续段，则可分段积分。除此之外，不管导流器表面的母线形状多么复杂，都可把它视作由很多层阶梯折线组成，这时每一层阶梯都对应一组相关的x，r和环形面积f。因此，在每一环面上可求得一个平均动压头，以它乘以相应的环形面积则得该环面上的气动载荷。把各环面上的气动载荷加起来，就得到锥体导流器垂直于气流方向的总气动载荷。这时有下列计算式：

式中，p0为锥体导流器的总气动载荷；n为微环面的数目（包括中心圆在内）；qi为第i个环面上的平均动压；fi为第i个环面的面积。

图3.36　锥形导流器气动载荷的确定

如果导流器不是圆锥体而是多面锥体，只要仿照上述原理在每个锥面上求出各阶层的平均动压，然后乘以各相应阶层的面积，最后求和，仍可近似地求得总气动载荷。

如果要求板面的压力分布，可近似地将板面上各点的ρu2乘以各点顺气流方向的切线与气流轴线夹角的正弦的平方。设夹角为αi，则各点上垂直于板面的压力为（ρu2）i·sin2αi。还要指出的一点是，射流中物体的气动载荷随着距喷口距离的变化而变化。在某些情况下，气动载荷的最大值出现在某一定的距离或时间上，这就要求计算出不同距离上的气动载荷，从中找出最大值。

（2）典型导流器受射流冲击时气动载荷的近似计算法。

①单面楔形导流器表面作用力的确定。

图3.37所示是一种单面楔形导流器，其特点是气流最后以垂直于原来气流的方向导出。设发动机的推力p＝mu，此处m为燃气流每秒的质量流量，u是燃气流的有效排气速度。对一定的火箭发动机，推力p是已知值，该力以气流冲击力的形式作用于导流器表面。表面直线段和圆弧段的水平分力和垂直分力分别按下列各式确定：

图3.37　单面楔形导流器表面作用力的确定

有了这些力，再根据射流边界的尺寸，以及射流横截面上的动压分布规律，即可把它们分布在导流器表面的相应面积上。据此就可进行导流器的结构强度分析。

上面给出的是单面楔形导流器，因此，水平分力Q1和Q2将使导流器产生横移。为此，在使用中要对导流器进行约束固定。如果是双面楔形式，而且燃气流对称地作用在两面上，则横移力可互相抵消。

根据图3.37所给出的载荷计算原理，对多发动机（如4个）以及单发动机其燃气流从双面楔形的顶部对称地喷下等各种情况，都可以类似地确定导流器所受的载荷。

②双面楔形平直导流器表面作用力的确定。

图3.38所示是一种双面楔形平直导流器及其受力情况。在某些发射装置上，为减小燃气流对发射装置及其后方设备的冲击和烧蚀作用，有时采用这种导流器。按实际作用，也可称它为分流器。其特点是气流被一分为二，而且各自折转的角度正好等于楔顶角的一半。其所受载荷的确定可按下述原则进行：

发动机推力分成两部分，分别作用在导流器的两个平直板面上。按公式即可算出Q和R。注意这时要用p／2代替原来的p。同时，因为导流板面没有弯曲段，燃气流的部分，即（1／2）mu cosω沿板面一直冲向导流器的后部空间，可以把它看作对导流器不产生作用力，因此没有类似图3.37所示的Q2和R2诸分力。

图3.38　双面楔形平直导流器表面作用力的确定

2.数值模拟方法

目前燃气射流的研究方法主要采用数值模拟方法，所以采用合适的分析模型和数值方法研究更符合燃气射流实际情况的流场，这仍然是重要的研究课题。篇幅所限，本小节仅对数值模拟方法作简单介绍。

燃气射流的计算方法大致可分为以下几种：特征线法、空间推进法、时间相关法和自适应网格法。用特征线法求解时，略去了黏性混合效应、化学动力学效应及两相流效应，控制方程简化为无黏双曲型方程，进而用特征线法求解。由于它只适用于双曲型方程，故应用受到了限制。20世纪70年代末80年代初，燃气射流的计算大多采用空间推进法，用该法求解时，首先把Navier－Stokes方程（简称为N－S方程）对应于流向坐标抛物化，得到抛物化的N－S方程。空间推进法对射流远场应用得很成功，但它无法模拟射流近场的复杂波系结构。时间相关法可直接对整个射流流场内的N－S方程进行数值求解，可以很好地模拟射流的近场流动，并可以对非定常流动进行计算，但它对计算机要求较高。自适应网格法适用于网格总数不变的情况，物理梯度大的地方网格能自动加密，物理量变化缓慢的地方网格相应变疏。针对燃气射流流场的数值计算，国内外工作者已做过大量的研究工作。

1）控制方程离散

流体力学基本方程可以通过不同的方式离散。最常用的是有限差分法、有限元法以及有限体积法，目前主流的计算流体力学软件包均采用有限体积法。

有限差分法是将求解域划分为差分网格，用有限个网格节点代替连续的求解域，然后将控制方程的导数用差商代替，推导出含有有限个未知数的差分方程组，求差分方程组的解，也就是微分方程定解问题的数值近似解，这是一种直接将微分问题变为代数问题的近似数值解法。这种方法发展较早，比较成熟，较多用于求解双曲型和抛物型问题，用它求解边界条件复杂，尤其是椭圆形问题不如有限元法或有限体积法方便。

有限元法是将一个连续的求解域任意分成适当形状的许多微小单元，并于各小微单元分片构造插值函数，然后根据极值原理（变分或加权余量法），将问题的控制方程转化为所有单元上的有限元方程，把总体的极值作为各单元极值之和，即将局部单元总体合成，形成嵌入指定边界条件的代数方程组，求解该方程组就得到各节点上待求的函数值。有限元法具有广泛的适应性，特别适用于几何及物理条件比较复杂的问题，而且便于程序的标准化，对椭圆形方程问题具有更好的适用性。但是其求解速度较有限差分法和有限体积法慢。

将求解域划分为有限小体积单元的思想也可以引入有限差分解法中，使之更适合求解复杂几何形状边界的流场。但在每个单元内不妨只做简单的积分，而不必像有限元法那样进行较烦琐的加权积分。这种离散法称为有限体积法。

有限体积法和有限单元法一样，可用于任意形状的体积单元，它无须有限差分法所要求的结构化计算网格，这使它具有应用的广泛性和灵活性；而有限体积法的基本方程表述了单个体积内的力学守恒关系式，它比有限单元法所采用的加权余量法更贴近力学连续介质的运动规律。

2）流场数值计算方法

（1）算法。

流场计算的基本过程是在空间上用有限体积法或其他类似方法将计算域离散成许多小的体积单元，在每个体积单元上对离散后的控制方程组进行求解。流场计算方法的本质就是对离散后的控制方程组的求解。对离散后的控制方程组的求解可分为耦合式解法（coupledmethod）和分离式解法（segregatedmethod），归纳后如图3.39所示。

（2）线性化方法。

无论是用分离式算法还是用耦合式算法求解，离散的非线性控制方程首先要被线性化以获得每个计算网格内的独立变量的线性方程组。通过求解线性方程组，获得最新的变量值。线性化的方法有两种，即隐式（implicit）方法和显式（explicit）方法。

隐式方法指的是对于一个给定的变量，使用相邻网格单元内的已知和未知变量求解每个网格单元内的未知量。这样一来，每个未知变量将会出现在不止一个方程中。因此，这些方程必须同时求解以获得未知量的值。

图3.39　流场数值计算方法分类

显式方法指的是对于一个给定的变量，每个网格单元内的未知变量的计算只使用相邻网格单元内的已知变量。因此，每个未知变量将只出现在一个方程中。因此，一次可求解一个方程以获得未知变量的值。

3）初始条件和边界条件

初始条件的给定和边界条件的选取是求解流体力学问题的必要条件。边界可以分为物理边界和人工边界两种。物理边界由问题的性质决定，从而是固定的；人工边界是由于计算所采用的范围永远小于实际区域而人为引入的。下面给出部分初始条件和边界条件。

（1）初始条件。

在对射流进行数值计算时，初始条件可以按以下方法给出：

①对于包含具体形状喷管的流动，可以直接取燃烧室的总温、总压值作为喷管入口的初始条件进行计算；

②对于喷管出口开始的流动情况，首先根据具体的喷管形状以及燃烧室的总温T0、总压值p0，通过一维定常等熵流动关系式，得到喷管各个截面的物理量分布，将这些值作为喷管出口的初始条件；

③在计算的其他区域，一般采用伴随流条件作为初始条件，即初始条件取上游边界值。

（2）边界条件。

在对射流进行数值计算时，涉及的边界条件有：壁面边界条件、压力入口边界条件、压力出口边界条件、轴边界条件、对称面边界条件、压力远场边界条件等。

①壁面边界条件。

火箭、导弹表面，发动机表面，导流器，发射装置等处采用壁面边界条件。在壁面边界条件中，壁面边界采用无滑移壁面边界条件，近壁面紊流计算采用标准壁面函数法处理。

标准壁面函数是基于Launder和Spalding的提议建立的一种半经验公式，考虑了壁面和完全紊流区域之间的黏性影响，避免了修改紊流模型，并满足一定的精度要求，因此得到广泛的应用。

②压力入口边界条件。

对于不可压流，压力入口总压和静压的关系由伯努利方程确定：

对于可压流，压力入口总压和静压的关系通过如下等熵关系式确定：

③压力出口边界条件。

喷管射流计算过程中所取的外场区域的边界采用压力出口边界条件。对于亚声速情况来讲，只需指定一个出口静压；对于超声速情况来讲，压力远场的压力通过流场内部点的数值采用二阶外推的方法获得。

④轴边界条件。

对于轴对称几何体，其轴线指定为轴边界条件。轴上任意一点处的变量的值取其相邻网格单元的值，无须为它定义边界处的数值。

4）常用CFD软件

目前常用的CFD商业软件有FLUENT、CFX、STAR－CD、PHOENICS等，它们的特点是功能比较全面、适用性强，几乎可以求解工程中的各种复杂问题；具有比较易用的前后处理系统和与其他CAD及CFD软件的接口能力，便于用户快速完成造型、网格划分等工作；同时，还可让用户扩展自己的开发模块；具有比较完备的容错机制和操作界面，稳定性高；可在多种计算机、多种操作系统，包括并行环境下运行。

此外还有包括OPENFOAM、GERRIS、BASILISK等开源代码可供选择，这些代码已经可以提供接近商业软件的使用体验，同时用户可从底层开始开发自己的模块，灵活性更大。