高等气体动力学第7章－超声速流动

高等气体动力学2008.3.04第七章 超声速流动7 超声速流动（6学时）7.1 二维定常无黏、无旋、等熵超声速流动的特征线法7.2 二维定常无黏、有旋、等熵超声速流动的特征线法7.3 三维定常无黏、等熵超声速流动的特征线法7.4 二维非定常无黏、等熵超声速流动的特征线法7.5 特征线数值计算的预估校正算法u 本章以三维定常无黏、等熵超声速流动的特征线法为例，讨论超声速流动的特征线理论、求解方法以及在微推偏喷管中的应用。u 有关特征线的一般理论已在第4章引入，不再重复。u 对其它流动情况，如二维定常无黏无旋等熵超声速流动、二维定常无黏有旋等熵超声速流动以及二维非定常无黏等熵超声速流动的特征线法，请参阅霍夫曼气体动力学。73三维定常无黏等熵超声速流动的特征线法一、原始控制方程组p 在笛卡尔坐标系中，无黏可压流的三维定常等熵流动的控制方程组由动量方程（3个）、连续方程和声速方程（代替能量方程）所组成，即(7-3-1)式中(7-3-2)为流体中的声速，r、p、u、v和w分别是流体的密度、压强和三个速度分量。p 方程组(7-3-1)含有五个一阶齐次偏微分方程，当流体作超声速流动时，它是一个双曲型系统，因此可以应用特征线理论。二、特征方程和相容性方程p 因为五个未知量r、p、u、v和w是耦合在一起的，所以必须联立求解。p 仿照含两个方程的方程组情形，作一个线性组合的微分算子，即(7-3-3)式中s1，s2、，s5均为待定的任意函数。p 将(7-3-1)代入上式，合并同类项，整理得(7-3-4)p 以上式中偏导数的系数作分量，构造如下的矢量(7-3-5)于是，式(7-3-4)中的各项可以看成r、p、u、v和w的梯度在相应矢量上的投影，或称它的方向导数。即(7-3-4)式可改写成(7-3-6)这里表示u在矢量方向上的方向导数，等等。上式就是三维定常等熵流动的相容性方程。p 现在提出这样的问题，即三维空间中是否存在这样的面，使所有的(j=1, 2, , 5)都落在这个面上？如果这样的面存在，则称它为特征面，其法线称为特征法线，而方程(7-3-6)就变成了此面上的一个内算符，称为与这一特征面相对应的相容性方程。p 特征法线是一个三维矢量，定义为(7-3-7)因为所有的矢量，都在特征面上，故有(7-3-8)将5个矢量(7-3-5)式代入上式，展开后以任意参数为未知数，并写成矩阵形式为(7-3-9)式中，(7-3-10)是流体的速度矢量。p 为了使方程组(7-3-9)有非平凡解，其系数行列式必须为零，即(7-3-11)这就是确定特征面的方程。三、特征面1. 流面p 令方程(7-3-11)的第一个因子为零，有 (三重根)(7-3-12)因为特征法线n的实际大小并不重要，所以可将其当成单位矢量，即(7-3-13)p 以上两个方程，式(7-3-12)和(7-3-13)是确定特征法线所需的两个条件，而另一个条件则是任意的，因此特征法线在给定空间点上不是唯一的，而是有一个无限多的特征法线族，每一个特征法线对应一个特征面，称为流面，所以也就有无限多的流面。p 方程(7-3-13)是一个球面方程，其球心在流线上。方程(7-3-12)则是一个过球心的空间平面方程，平面的取向取决于三个速度分量u、v和w的大小。所以两个方程联立起来的解是此平面与球相交形成的以球心为圆心的圆，它的每一条半径代表一个特征法线，而每一个特征法线又对应一个特征面。p 平面方程(7-3-12)是所有特征法线所在的平面，称为法线平面，它垂直于流线；法线平面上的每一条过球心的线，即特征法线对应一个流面，该流面是垂直于特征法线、并包含流线的面。 图7-3-1 三维定常流中流面的法线平面 图7-3-2 与一条特征法线对应的流面p 流面的包络是流线，若设t表示流体微团沿流线运动的时间，则流线的方程可写成，(7-3-14)p 在dt时间内，流体微团沿流线运动的微元距离矢量为(7-3-15)2. 波面p 令方程(7-3-11)的第二个因子为零，解之得(7-3-16)式中的正负号是任选的。p 由此可见，确定n的条件也只有两个，即方程(7-3-16)和(7-3-13)，而另一个条件是任意的。所以，过空间一点有无限多的特征法线，它们形成了一个锥，称为法线锥。p 和流面一样，方程(7-3-13)仍是一个球面方程，方程(7-3-16)仍是一个空间平面方程，但该平面不通过球心，而且平面的取向不仅取决于三个速度分量u、v和w的大小，还取决于声速a的大小。所以两个方程联立起来的解是此平面与球相交形成的圆，其圆心与球心不重合。p 从球心发出、并指向此圆的每一条半径代表一个特征法线，所以全部法线组成了以球心为顶点的锥，即法线锥。这说明，法线锥的每一条母线都是特征法线，而每一条特征法线又对应一个特征面，这个特征面称为波面。p 显然，无数条特征法线就对应于无数个特征面，即无数个波面。这些波面的包络称为马赫劈锥，它的顶点也在流线上。与马赫劈锥在锥顶点相切的锥是马赫锥。马赫锥代表均匀流场中的马赫劈锥，所以它是一个直立圆锥，并且是法线锥的倒锥。p 方程(7-3-16)表明，垂直于特征面的流动速度的大小等于声速a。 图7-3-3 波面的法线锥 图7-3-4 与两条特定特征法线对应的波面p 以微分形式表示的马赫劈锥方程可写成(7-3-17)式中，。p 不难验证，只有在超声速流动时，即V2-a20时，由上式给出的马赫劈锥才是实的。p 马赫劈锥与波面相切的曲线称为双特征线。图7-3-5表示在类空间曲线上任一点P(x,y,z)上的两族波面、流面、马赫劈锥和双特征线之间的关系。类空间曲线是指马赫劈锥以外的任意曲线。如果此类空间曲线旋转360，则无限多个流面的包络就形成了流线，而无限多个波面的包络形成了马赫劈锥。p 马赫锥与波面相切的线称为双特征射线，它是直线。显然，双特征射线是双特征线的切线。p 如图7-3-6所示。若以表示流体微团沿流线在dt时间内运动的距离矢量（它是马赫锥的轴），表示声波在dt时间内沿传播的距离矢量，则在此时间内，沿双特征射线运动的距离矢量为。矢量是沿双特征射线的方向，定义为(7-3-18)所以，沿双特征射线运动的距离矢量可写成(7-3-19)这里dt是流体微团沿流线运动的时间。 图7-3-5 三维定常流中的流面、波面、 图7-3-6 运动或传播距离间的流线、马赫劈锥和双特征线 几何关系四、相容性方程1.在流面上成立的相容性方程p 将流面的法线方程(7-3-12)代入(7-3-9)式得(7-3-20)p 该方程组系数矩阵的秩为2，因此它有三个独立的解，其中s5可任意指定，s4=0，而s1，s2和s3则必须满足如下关系(7-3-21)p 满足方程(7-3-20)的三个可能的解是，(7-3-22)，s5任意(7-3-23)，(7-3-24)式中，是流面上与速度矢量不平行的任一矢量。p 在这三个解中，前两个，即式(7-3-22)和(7-3-23)对所有流面都成立，即在流线上是成立的。p 由第三个解(7-3-24)给出的si (i=1, 2, , 5)依赖于特定的，因为对每一个与特定流面对应的所得到的不同，所以相应的相容性方程只适用于这一特定的流面，因此通过改变可以得到无限多个相容性方程。p 将方程(7-3-22)(7-3-24)代入(7-3-4)式，得到下面三个在流面上成立的相容性方程(7-3-25)(7-3-26)(7-3-27)以上方程，式(7-3-25)(7-3-27)是流面内的微分算符。p 因为在三维空间中的一个面内，仅有两个独立的微分方向，所以用对三个笛卡尔坐标方向的偏导数写成的(7-3-25)(7-3-27)式可以用对流面内任意两个独立方向的偏导数来改写。事实上，方程(7-3-25)和(7-3-26)可用流线上的全导数来表示。考虑一般函数f(x, y, z)，其全微分为(7-3-28)由(7-3-15)式，在流线上有，所以(7-3-29)或(7-3-30)于是，对比上式，可将(7-3-25)和(7-3-26)改写成 (沿流线)(7-3-31) (沿流线)(7-3-32)式中，等等，t表示流体微团沿流线的运动时间。p 方程(7-3-27)不能写成上述沿流线的形式。因为根据定义，是流面内与速度矢量不平行的方向。2. 在波面上成立的相容性方程p 将方程(7-3-16)代入(7-3-9)式(选取正根)得(7-3-33)p 此方程组系数矩阵的秩是4，所以只有一个独立解。令s5 = 1可以得如下解，(7-3-34)p 将此解代入方程(7-3-4)可得波面上的相容性方程(7-3-35)该方程表示了在波面上的微分算符。p 与流面相容性方程一样，可以将它改写成对波面内任意两个独立方向的偏导数，这两个独立方向可以选取双特征线方向(即与马赫劈锥的射线相切的线)和波面内与双特征线垂直的方向。将方程(7-3-19)代入(7-3-28)得一般函数f对t的全导数(7-3-36)这里t仅表示流体微团沿流线运动的时间，而ft则表示f沿双特征射线的方向导数。p 在方程(7-3-35)中加减如下的项整理后与(7-3-36)比较可得(7-3-37)这里，等等。在上式中，方括号中的项是用对三个笛卡尔坐标的偏导数来表示的双特征线垂直方向上的方向导数，称为交叉导数。五、小结p 通过推导，得到了在每一个流面上成立的、独立的三个相容性方程(7-3-27)、(7-3-31)和(7-3-32)，在每一个波面上成立的、独立的一个相容性方程(7-3-37)。p 由于在空间中的一个点上有无数的流面和波面，因而也就有无数的相容性方程。p 但是，在用相容性方程组代替原始控制方程组时，独立的相容性方程的个数不能超过独立的原始控制方程的个数。p 因此，有必要确定相容性方程的可能组合哪些是独立的。Rusanov给出了三维定常等熵超声速流动的可能组合，几个可能的独立的相容性方程组之一是由沿流线的两个方程(7-3-31)和(7-3-32)及沿任意三个波面(三条双特征线)的相容性方程(7-3-37)组成的，这也就是在数值模拟喷管流动时所采用的相容性方程组。为了方便，将它们列在下面并重新编号。1. 特征方程流线，(7-3-38)马赫劈锥(7-3-39)2. 相容性方程沿流线(7-3-40)(7-3-41)沿双特征线(7-3-42)75特征线数值计算的预估校正算法p 在用特征线方法计算流场时，数值积分网格或者数值格式可以用直接步进法，也可以用逆步进法。p 在直接步进法中，解点的位置是从初值面上已知基点向前(下游)延伸数值积分网格的生成元来确定的；p 在逆步进法中，解点位置用某种方法预先给定，而将数值积分网格的生成元向后延伸与初值面相交，基点上的流动参数通过在初值面上进行内插来确定。p 一般情况下，直接步进法精度较高，计算省时，而逆步进法则编制程序容易，且解面位置可以按需要设置，这对分析流场是有利的。p 在计算中，数值积分网格可以用特征面也可以用双特征线来构成。特征面网格所得到的有限差分网格比较简单，而双特征线网格通常有较高的精度。p 本节讨论双特征线逆步进解法，求解方程组由流线上的两个方程(7-3-40)、(7-3-41)和三条相互间隔120的双特征线上的相容性方程(7-3-42)所构成。计算是在一系列平行平面上进行的，对喷管流场而言，这些平行平面(即解面)垂直于喷管对称轴。一、数值计算方法p 数值积分网格的构成如图7-5-1所示。在未知解面上，解点5的轴向位置即x坐标由稳定性条件来确定，然后数值积分网格生成元向后(上游)延伸与初值面(或先前的解面)相交，三条间隔120的双特征线与初值面的交点1、2和3称为基点。在此网格中，流线用4点和5点的连线来代替，而三条双特征线则用连线1-5、2-5和3-5来近似表示。图7-5-1 逆置双特征线网格p 为使计算达到二阶精度，采用欧拉预测校正法求解相容性方程组，并进行全场迭代。1. 预测步在预测步中，进行如下步骤的计算。(1) 确定解点位置 因为解面是一系列平行平面，且垂直于x轴，其间距Dx是由稳定性条件确定的，所以解点5的x坐标为(7-5-1) 由于连线4-5是流线，所以根据流线方程(7-3-38)，流体微团从4点运动到5点的时间为(7-5-2) 于是，由流线的另外两个方程可得解点的y和z两个坐标，或，(7-5-3)(2) 预定基点位置 基点1、2和3的位置可以通过求解马赫劈锥方程(7-3-39)来得到，即(7-5-4)这里，Qk是初值面上第k个基点与4点连续的斜率，是根据双特征线1-5、2-5和3-5相互间隔120来确定的。(3) 确定基点的流动参数 已知基点的位置后，用4点及其周围的点构造一个双变量(y和z)二次插值公式，通过插值计算出基点上的流场参数值，即uk、vk、wk、pk和rk，k=1,2,3。同时，对插值公式微分，得到基点上的交叉导数，即uy、uz、vy和vz等。(4) 修正基点位置 已知基点上的流动参数后，再用马赫劈锥方程对基点的位置做一次修正，即(7-5-5)(5) 确定解点流动参数 解点5上的流动参数是通过求解相容性方程组来得到的。将方程(7-3-40)(7-3-42)进行离散得沿流线(7-5-6)(7-5-7)沿双特征线(7-5-8)式中，。于是，联立求解方程(7-5-6)(7-5-8)可得到u5、v5、w5、p5和r5。2. 校正步p 校正步的计算步骤与预测步的相同，所不同的是在方程(7-5-2)(7-5-8)中，以基点参数表示的系数均应换成解点5和相应基点1、2、3、4之间的平均值。3. 迭代p 在完成一个解面计算后，根据计算需要可以进行全场迭代，即对每一点重复进行上述过程，以提高计算精度。一般情况下，不超过5次迭代即可满足精度要求。二、稳定性条件与计算的细节问题1. 稳定性条件p 在数值积分的步进过程中，两个相邻解面之间的距离Dx必须调整到使解面上的所有网格点都满足稳定性条件，即CFL稳定性准则。p 对解面上的每个网格点，可允许的步长大小是当地流动参数和网格点间隔的函数，可写成(7-5-9)式中，V是流动速度，u是步进方向即x方向的流动速度分量，Rmin是流线与初值面的交点即4点与其周围的点中离得最近的那个点之间的距离。参数C是马赫锥在与流线相垂直的那个平面中的扩张速度，即(7-5-10)p 因此，在求解时，每一个解面上的每个点都对应有一个Dx，为了既保证稳定又要使解面相互平行，应选取解面上所有Dx的最小值作为该解面与初值面的间距。通常，在确定解面位置即Dx时，解面上的流动参数是未知的，所以(7-5-9)和(7-5-10)也需要参予整个迭代过程。2. 计算交叉导数p 前面已提到，有关、的交叉导数是通过微分插值公式得到的，而对x的偏导数则需要从原始方程中求出。p 由方程(7-3-1)可直接导出(7-5-11)将u乘以x方向的动量方程再减去声速方程可得(7-5-12)上述方程对预测步和校正步都适用。3. 边界条件p 在计算过程中，当解点5位于固体壁面(如喷管壁)上时，基点1已不在流场之内，因此相容性方程中不再计及k=1的情况。同时由方程(7-5-3)所确定的解点坐标还必须满足边界的曲面方程。由于解点5位于边界上，其速度矢量必与边界相切，因此速度矢量V5垂直于边界的法矢量，即或(7-5-13)这里是固体壁面的法线矢量。p 于是，用上式代替k=1的相容性方程，从而使相容性方程组在边界上仍是封闭的。三、用三维特征线方法计算喷管的超声速非对称流场p 在固体火箭发动机中，由于推进剂燃烧的不均匀性或者喷管几何型面的不对称性，都会使喷管内的燃气流场呈现非对称流动。这种非对称流动将使火箭推力矢量偏离火箭的对称轴线，从而出现作用于火箭的侧向力和推力矢量绕火箭质心的力矩，这就是推力偏心现象。它是影响无控火箭射击精度的主要因素之一。p 用特征线方法数值模拟喷管中的三维定常等熵非对称超声速流场，并通过流场分析获得推力偏心特性的数据，在微推偏喷管的研制中已经得到了成功的应用。p 三维特征线方法只适用于超声速流场，因此只能计算喷管扩张段的流动。在扩张段中选取某一马赫数大于1、垂直于对称轴的截面作为初值面，则其下游的三维流动可以用特征线法进行求解。p 为了模拟非对称流动，假设各种非对称因素引起的非对称扰动等效于初值面上的流动速度方向偏离其均匀对称流动一个角度db，此速度偏离将引起y方向的推力分量，而在z方向流动仍是对称的。1. 锥型喷管p 计算时假设燃气为理想气体，气体参数为：比热比g =1.2216，速度倾角db =1，气体常数R=320.0J/(kg.K)。p 初值面上网格点数为181个，所以喷管扩张段的超声速流场由181条流线组成。根据流场计算结果，积分计算出喷管的推力矢量分量。x方向的主推力和y方向的侧向力沿喷管轴向的变化如图7-5-2和图7-5-3所示。图中，轴向坐标是，纵坐标分别为和，这里为喷管的质量流率。 图7-5-2 发动机主推力沿轴向的分布 图7-5-3 y方向侧向力沿轴向的分布p 从图中可见，侧向力沿轴向具有明显的摆动特性(对不同扩张半角喷管进行的计算表明，扩张半角越大其摆动周期越长)。随着喷管长度的增加，侧向力的数值迅速衰减。从这种发展趋势上看，当喷管扩张段为无限长时，侧向力将趋向于零，从而也就没有推力偏心，至少在理论上可以得出这样的结论。p 但是，在实用上喷管长度总是有限的，并且若设计不当其扩张段长度可能对应着较大的侧向力数值，从而火箭的射击密集度下降。p 因此如果在设计喷管时，将喷管扩张段长度选择在其出口截面恰好位于侧向力为零的位置，就可以极大地减小推力偏心，这正是微推偏喷管设计的基本原理。所以，从理论上并最终通过试验，确定侧向力的零点位置(特别是第二零点)具有重大的实际意义。p 图7-5-3中也同时给出了用六分力推力试验台测出的侧向力试验数据。从对比情况来看，侧向力计算值的大小与试验数据的差异较大。这是因为，在计算中所施加的非对称扰动很难与真实喷管流动的非对称性在数值上保证相当。根据已往的实验研究，流场非对称性大小只对侧向力的数值大小即侧向力曲线的振幅有决定性影响，而对侧向力曲线的零点位置的影响是极微弱的。因此，在对比计算结果和试验数据时，只关心零点位置的吻合情况。p 从图上看，两者在零点位置上的吻合程度是令人满意的。其微小偏差就数值计算方法而言，可以认为是由以下原因造成的：数值积分网格不够细；理想气体假设与真实燃烧产物的偏差；特别是初值面的布置与真实情况有很大不同。p 关于非对称扰动的大小对侧向力曲线零点位置仅有微弱影响的结论，从理论上也可以证明。图7-5-4给出了速度偏离角分别为db =1、1.5、2.0、2.5、3.0和3.5时的侧向力曲线。从图上可以看出，db的大小对零点位置几乎没有影响。而且，当db=2.0左右时，侧向力计算值与试验值相当接近。图7-5-4不同速度偏离角时侧向力沿轴向的分布这里2. 特型喷管p 给出了两种特型喷管的