CN111859520A - 一种计算高超声速风洞轴对称喷管内型面的方法 - Google Patents

Abstract

本发明公开了一种计算高超声速风洞轴对称喷管内型面的方法，包括以下步骤：步骤一、给定初始条件参数；步骤二、计算轴对称喷管的亚声速型面坐标；步骤三、计算轴对称喷管的超声速型面坐标，包括：计算喉道区型面参数计算；沿径向流波前气流参数计算；转换区中心线上的马赫数分布计算；消波区参数计算；利用轴对称特征线网格计算型面；计算轴对称喷管壁面坐标；步骤四、对轴对称喷管的超声速段位流型面附面层进行修正计算。通过本方法设计的喷管型面，流场所际需要的马赫数较一致，喷管出口马赫数分布均匀，并且能比较方便的采用MATLAB软件编程实现，快速的得出型面坐标，在计算轴对称喷管型面时具有普适性。

Description

一种计算高超声速风洞轴对称喷管内型面的方法

技术领域

本发明属于流体力学技术领域，更具体地说，本发明涉及一种计算高超声速风洞轴对称喷管内型面的方法。

背景技术

轴对称喷管的气动设计主要包括两部分：位流型面(亚音速段和超音速段)和超音速位流型面附面层修正。位流又分为三个区域：第一区域是亚音速区；第二个区域是由喉道到拐点的喉道区；第三个区域是由拐点到喷管出口的转换区。

在喉道，外形由单调递增的三次曲线给出。在转换区，外形主要由Cresci方法计算，如图1所示，假定流动从0′源点出发，沿径向流一直膨胀到波前AB，在拐点A处具有最大膨胀角—半锥角θ_A的锥形流。CE下游是均匀平行流区，在ABCE区域中流动由径向流过渡到平行流，利用波前AB和BC上的已知条件，用轴对称特征线网格计算。最后得到TE位流型面。

在给定θ_A的值，如果B、C重合，这样得到最短喷管。不过在这样的喷管中，沿轴向马赫数梯度在B点不连续，其曲率半径在拐点处不连续，理论上，曲率不连续会导致在做附面层修正时喷管的附面层为负；如曲率是连续的，附面层的增长也就是光滑的，况且，即使附面层增长被准确计算出来，其喷管型面的粘性修正也只是针对设计运转条件的，那么修正量也只有在计运转条件下才能够应用到非粘性外形上去，所以离开设计运转条件，附面层增长的偏差将产生不连续，这些不连续将顺着喷管因机械加工不光滑的地方向下移动，从而导致试验段的流动中存在扰动。若外形有连续的曲率，这些设计运转条件就可以减弱一些。

发明内容

本发明的一个目的是解决至少上述问题和/或缺陷，并提供至少后面将说明的优点。

为了实现根据本发明的这些目的和其它优点，提供了一种计算高超声速风洞轴对称喷管内型面的方法，包括以下步骤：

步骤一、给定内型面的初始条件参数，利用初始条件参数进行轴对称喷管内型面计算；轴对称喷管内型面计算主要包括位流型面计算和超声速位流型面附面层修正；其中位流型面又包括：亚音速段型面和超声速段型面；流经喷管内型面的位流包括：亚音速区、从喉道到拐点的喉道区和从拐点到喷管出口的转换区；

步骤二、计算轴对称喷管的亚声速型面坐标；

步骤三、计算轴对称喷管的超声速型面坐标，包括：计算喉道区型面参数计算；沿径向流波前气流参数计算；转换区中心线上的马赫数分布计算；消波区参数计算；利用轴对称特征线网格计算型面；计算轴对称喷管壁面坐标；

步骤四、对轴对称喷管的超声速段位流型面附面层进行修正计算。

优选的是，其中，所述步骤一中给定的初始条件参数包括：型面马赫数M，喷管入口高度Y_in，喷管出口高度Y_in，喷管长度L，型面马赫角α，喷管壁面总压强P₀，喷管壁面总温度T₀，亚声速段移轴后的喷管进出口半径高度比n，喷管壁温T_w；同时在喷管上引入转换区，标记为ABCD；喉道区型面标记为TA；消波区，标记为BCE。

优选的是，其中，所述步骤二中计算轴对称喷管亚声速型面坐标的方法包括以下步骤：

采用维托辛斯基公式计算亚声速端型面，计算按以下公式进行：

上式中，

Y₁′＝Y₁+b；Y_*′＝Y_*+b；Y＝Y′-b；

其中，Y₁、Y_*和Y分别为收缩段的入口、喉道的半高度和X轴上的截面半径，Y′₁、Y′_*和Y′分别为移轴后收缩段的入口、喉道的半高度和X′轴上的截面半径，a为特征长度，b为对X的移轴量、n表示亚音速段移轴后的进出口半径高度比，即n＝Y′₁/Y′_*。

优选的是，其中，所述步骤三中计算轴对称喷管的超声速型面坐标的方法包括以下步骤：

步骤S31、计算喉道区型面TA参数；假设喉道区为一元等熵流，经验证明喉道区型面用一元三次方程能满意的获得要求的气流流动状态，该一元三次方程如下：

y＝a₀+a₁x+a₂x²+a₃x³ (2)

系数a₀、a₁、a₂、a₃由边界条件确定，边界条件如下：

x＝0y＝y_*；

x＝x_Ay＝y_A；

将(3)式代入(2)式经整理得：

a₀＝y_*；a₁＝0；

将(4)式代入(2)式整理后得：

(6)式为所求的喉道区型面曲线，其中θ_A为给定的经验值，x_A由下述关系式确定：

由边界条件得：

x_A＝3(y_A-y_*)(2 tanθ_A)^-1 (6)

由

可知：

于是将(7)、(8)式代入(6)式，求出x_A为：

从连续方程：

将(8)式代入(11)式得：

(12)式中：

r为从0′到任一点的距离；

引入函数F(M_E)对喷管喉道面积做不完全气体影响的修正，引入函数：

则：

到此，只要知道M_A和F(M_E)，喉道区型面曲线TA可由(5)式计算出；

步骤S311、计算F(M_E)：

不完全气体的比热比有：

a_tp ²＝RTγ_tp (16)

a_*tp ²＝RT_*γ_*tp (20)

将(16)～(21)式代入下式：

而

由M_tp＝1和(M_E)_tp为喷管出口设计马赫数代入(19)式分别求出T_*和T_E，代入(22)获得(A/A_*)_tp和(23)式一起代入(14)求得：

对给定的T₀由M_E就可以计算得F(M_E)；

步骤S312、计算M_A：

在轴对称等熵流中，P-M膨胀角有以下关系式：

ψ＝ψ_B-θ (26)

为求ψ_B必须选择M_B的值，由于M_B选择有任意性，有两种方法：

步骤(a)由Cresci的建议选取的M_B值比M_C低0.2即：

M_B＝M_C-0.2 (27)

步骤(b)根据BC线上的速度系数成三次方分布，并利用边界条件最后推得B点与C点的速度系数关系

式中：M_C＝M_E，或W_C＝W_E；

确定了M_B之后，按以下步骤求解M_A：

步骤(c)将M_B(W_B)代入(25)式求出ψ_B；

步骤(d)以ψ_B代入(26)式并将ψ及θ换成A点的ψ_A和θ_A的值，求出ψ_A；

步骤(e)将ψ_A代入(25)式求出M_A；

步骤S32、计算沿径向流波前AB气流参数：

在波前AB上有几何分割关系：

(29)式中：M_A≤M_P≤M_B；

N为做特征线网格时沿径向流波前AB的分割总数，从B到A，P＝0，1，2…N；

由M_P代入(25)、(26)式可求得θ_P值；

由几何关系：

(30)式中：

步骤S33、计算ABCD转换区中心线BC上的M数分布：

BC上的M数分布是连续变化的，假定速度系数成三次方分布：

(31)式中：

系数a₀、a₁、a₂、a₃由边界条件决定，边界条件为：

在B点：

在C点：

由(32)、(33)式解得：

a₀＝W_B；a₁＝3(W_C-W_B)；a₂＝-3(W_C-W_B)；a₃＝W_C-W_B (34)

将(34)式代入(31)式得：

且M和W的关系如下：

此处

当γ＝1.4时，α＝6，则

利用(32)和(35)式求得：

上式对W微分，且x＝x_B得：

最后有：

(40)式中：

步骤S34、计算消波区BCE参数：

利用步骤二中的初始条件，做特征线网格求解型面AD,当特征线网格做到DC，则以DC上的参数作初始条件，只要求得CE上的初始条件，用同样的方法作特征线网格，就可以解出型面DE，于是整个型面AE就确定了；

要使试验段获得均匀的平行气流，DCE为减波区，CE必须为一条直线，在CE上的M数均为M_E且CE上的参数为：

(41)式中，N为做特征线网格时沿径向流波前AB的分割总数，P＝0，1，2…N；

步骤S35、计算轴对称特征线网格：

轴对称喷管中，气流特性是对称于中心线的，可以只研究通过中心轴的xy平面的流态就能决定整个喷管的流态；

这里用轴对称特征线网格计算型面AD和DE的方法，喷管的对称轴为x轴，气流沿x轴是对称的，在物理平面xy上设已知点P₁(x₁，y₁)和P₂(x₂，y₂)，P₁(x₁，y₁)的左特征线和P₂(x₂，y₂)的右特征线交于P₃(x₃，y₃)；

由几何关系：

(42)式中：

(42)式经简化整理得：

y₃＝y₁+(x₃-x₁)tan(θ₁+α₁) (44)

或

y₃＝y₂+(x₃-x₂)tan(θ₂-α₂) (45)

由轴对称特征线方程：

(47)式中：上边符号对应于左特征线，下边符号对应于右特征线，ds_1,2为左右特征线微元长度，把上式改写为点1、2、3的差分方程为：

(47)式中：

由(47)式解得：

(48)式的下标13，23是指取P₁和P₃，P₂和P₃的平均值；

当P1和P2位于X_BC上，则有y₁＝y₂＝θ₁＝θ₂＝0，(48)式右边第二项为不定式，出现奇异点，若在X_BC附近近似为径向流有：

故(47)、(48)式变为

为了提高计算精度，采用迭代法计算，迭代式如下：第一先用P₁、P₂点代入以上各式计算出P₃的值，以后各次可用下列各式子进行迭代计算：

将(50)式代入(48)或(49)，直到达到5×10^-6精度为止；

步骤S36、计算轴对称喷管壁面坐标：

当特征线做到边界上的点A后，为了要确定流线上的点R，必须对边界条件作一些处理，求出壁面上的坐标；利用边界上流线和特征线的关系，设流线为AR，13和23为1、2点发出的左右特征线；

线性插值：

得：

由流线和特征线的关系：

由(53)解出：

将x_R代入

得：

线性关系：

x_R由(54)式给出；其中，下标1为前一个壁面点，初始壁面点为A点，下一个壁面点是以R为初始条件进行的，如此类推，最后完成AE型面的计算。

优选的是，其中，所述步骤四轴对称喷管超声速段位流型面附面层进行修正的计算方法包括以下步骤：

步骤S41、建立动量方程：

在高超音速喷管中，附面层属于湍流附面层，湍流附面层增长用冯·卡门轴对称动量积分方程来描述：

方程(59)左边最后一项只有在轴对称的情况下才出现，这里C_f、H可由湍流附面层的关系给出；将喷管的曲线坐标s换成轴对称坐标x，得：

方程(59)变为：

利用stwartson变换：

由动量厚度θ定义，并以(61)式代入得：

(62)式中：

由位移厚度δ^*定义，并利用温度型关系式：

其中，

当Pr＝1时，在绝热壁面条件下，壁面总温T_s＝T₀，则δ^* _tr由H定义及(62)、(66)式得：

最后把(63)、(67)代入(60)简化后得：

上列各式中：下标“0”表示自由驻点条件；下标“e”表示附面层底部条件；下标“tr”表示经Stewartson转换后的参数，初始条件：当x＝0时，θ_tr＝0；

为了求方程(68)，必须先求出气流与喷管壁面的摩擦系数C_f和附面层形状因子H_tr：

步骤S411、计算附面层形状因子系数H_tr：

采用修正的Crocco二次定律，可给出考虑了绝热壁温和Pr≠1影响的可压缩流附面层形状因子系数H_tr，附面层内的温度分布为：

其中(69)式中：T_w为壁温，T_aw为绝热壁温；由(61)、(69)式代入(65)式并简化得：

其中，

把(62)式代入(70)式得：

其中，

(67)、(71)联合解得：

只要给出H_i，就可以求出H_tr，然而H_i与M无关，只和不可压摩擦系数有关，由半经验公式给出：

这里系数7是实验值；

步骤S412、计算气流与喷管壁面的摩擦系数C_f：

利用Eckert给出的参考温度法，即参考温度T′为：

其中C′_f＝F(Re_x′)使用的是不可压情况下的C_fi＝F(Re_x)关系式；

在此

其中，

由无压力梯度的Karmán—Schoenherr的平均摩擦系数公式：

由此式从Re_xi解C_Fi是不容易的，引入它的近似关系式，及求局部摩擦系数：

把不可压的Re_xi换成Re_x′并代入(75)式得：

为了把不可压的C_fi应用到可压的情况，必须对方程(68)作如下变换：

把(78)、(79)代入上式并作近似计算得：

其中，

将(68)转换成不可压的情况，由Re→Re_x′→Re_X(＝Re_xi)求C_fi，Re_x′→C_f；

带“'”的参数是以T′为特征温度的值，因为高马赫喷管要进行强迫水冷，所以设计时，对于M_E＝10时，取T_w＝573K；M_E＝11时，取T_w＝623K，绝热壁温按下式计算：

式中：σ为复温系数，对于M_E＝5、6、7取σ＝0.88，对于M_E＝8、9、10、11、12取σ＝0.896，均为湍流状态；

步骤S42、实际计算要考虑的因素；

(1)Re数的范围

得到的马赫数与实Karmán—Schoeherr的公式，只适合没有压力梯度的情形，因此，应注意到(79)式的结果也是针对没有压力梯度的，而且(78)式为近似公式，当logRe_xi＝1.5或logRe_xi＝2.3686时，出现奇异点，以上公式适用在Re＝10⁵～10⁹范围，计算是满足要求的；

(2)假想原点

附面层的计算通常是以喉道为初始条件的，当x＝0时，求得C_f是发散的，为了消除奇异点，以喉部上游为假想原点，则喉部处x≠0，假想原点的x_*由下式求得：

这样原来的x坐标从喉道往上游移动x_*的距离，移轴后的坐标应为x_s＝x_*+x；

(3)附面层修正因子

为了得到预先给定的出口直径，在进行附面层修正时，要引入附面层修正因子f，修正后喷管型面坐标如下：

计算时，第一步利用无粘位流型面参数，计算出喷管出口的

然后通过(84)式计算出f₁，利用f₁进行计算，得到修正后的喷管型面参数Y₁、X₁，在利用Y₁、X₁求出

f₂直到满足y_E-Y_n≤10^-6mm为止；

通过以上步骤，可以得出高超声速喷管的气动型面，为喷管的设计加工提供内型面坐标。

本发明至少包括以下有益效果：通过本方法设计的喷管型面，流场所际需要的马赫数较一致，喷管出口马赫数分布均匀，并且该方法能比较方便的采用MATLAB软件编程实现，快速的得出型面坐标，在计算轴对称喷管型面时具有普适性。计算出可以投入生产加工的内型面曲线，该方法设计的喷管具有马赫数精准度高，流场品质好的优点。

本发明的其它优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。

附图说明：

图1为本发明喷管气动设计计算图；

图2为本发明亚声速段型面坐标；

图3为本发明物理平面上特征线与参数的关系；

图4为本发明便捷上流线与特征线的关系；

图5为本发明坐标转换关系示意图；

图6为本发明马赫数为5的喷管型面示意图

具体实施方式：

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

应当理解，本文所使用的诸如“具有”、“包含”以及“包括”术语并不配出一个或多个其它元件或其组合的存在或添加。

如图1-6所示：本发明的一种计算高超声速风洞轴对称喷管内型面的方法，包括以下步骤：

步骤一、给定内型面的初始条件参数，利用初始条件参数进行轴对称喷管内型面计算；轴对称喷管内型面计算主要包括位流型面计算和超声速位流型面附面层修正；其中位流型面又包括：亚音速段型面和超声速段型面；流经喷管内型面的位流包括：亚音速区、从喉道到拐点的喉道区和从拐点到喷管出口的转换区；

步骤二、计算轴对称喷管的亚声速型面坐标；

步骤三、计算轴对称喷管的超声速型面坐标，包括：计算喉道区型面参数计算；沿径向流波前气流参数计算；转换区中心线上的马赫数分布计算；消波区参数计算；利用轴对称特征线网格计算型面；计算轴对称喷管壁面坐标；

步骤四、对轴对称喷管的超声速段位流型面附面层进行修正计算。

在上述技术方案中，所述步骤一中给定的初始条件参数包括：型面马赫数M，喷管入口高度Y_in，喷管出口高度Y_in，喷管长度L，型面马赫角α，喷管壁面总压强P₀，喷管壁面总温度T₀，亚声速段移轴后的喷管进出口半径高度比n，喷管壁温T_w；同时在喷管上引入转换区，标记为ABCD；喉道区型面标记为TA；消波区，标记为BCE。引入转换区ABCD，喷管增长了，曲率连续可以满足，流场性能可以改善。在B点马赫数M_B是可以任意给定的(但要比C点低)，BC之间的M数分布由三次方程给出，方程中系数由连续条件确定。

喉道区外形用三次曲线计算，它是单调递增到具有θ_A的切点上，在A点上二阶导数为零，在喉道处满足连续条件。

亚音速收缩段是将稳定段来的气流均匀加速至声速。根据高超声速喷管的设计要求，到达喉部的声速流必须是均匀的。经验证明，如果稳定段来流是均匀的，只要有一条光滑连续而又渐变的收缩曲线就能基本满足要求。风洞试验段M数不同，喉部面积也不同，因而收缩比是随M数而变化的，M数越高，收缩比越大。亚音速段外形采用移轴的维托辛斯基公式，并由喉道处的连续性条件来确定移轴量。表1为几组高超声速风洞喷管内型面计算条件参数。

表1高超声速风洞喷管内型面计算条件参数

在上述技术方案中，所述步骤二中计算轴对称喷管亚声速型面坐标的方法包括以下步骤：

采用维托辛斯基公式计算亚声速端型面，高马赫数时，型面曲线很陡，必须进行适当的移轴，才能很好的满足要求；如图2所示，为了保证喉道处气流均匀和连续，在应满足连续条件：即收缩段喉道处的曲率半径等于附面层修正后的超声速段喉道处的曲率半径，由此确定移轴量，计算按以下公式进行：

上式中，

Y₁′＝Y₁+b；Y_*′＝Y_*+b；Y＝Y′-b；

其中，Y₁、Y_*和Y分别为收缩段的入口、喉道的半高度和X轴上的截面半径，Y′₁、Y′_*和Y′分别为移轴后收缩段的入口、喉道的半高度和X′轴上的截面半径，a为特征长度，b为对X的移轴量、n表示亚音速段移轴后的进出口半径高度比，即n＝Y′₁/Y′_*。

在上述技术方案中，所述步骤三中计算轴对称喷管的超声速型面坐标的方法包括以下步骤：

步骤S31、计算喉道区型面TA参数；假设喉道区为一元等熵流，经验证明喉道区型面用一元三次方程能满意的获得要求的气流流动状态，该一元三次方程如下：

y＝a₀+a₁x+a₂x²+a₃x³ (2)

系数a₀、a₁、a₂、a₃由边界条件确定，边界条件如下：

x＝0y＝y_*；

x＝x_Ay＝y_A；

将(3)式代入(2)式经整理得：

a₀＝y_*；a₁＝0；

将(4)式代入(2)式整理后得：

(7)式为所求的喉道区型面曲线，其中θ_A为给定的经验值，x_A由下述关系式确定：

由边界条件得：

x_A＝3(y_A-y_*)(2 tanθ_A)^-1 (6)

由

可知：

于是将(7)、(8)式代入(6)式，求出x_A为：

从连续方程：

将(8)式代入(11)式得：

(12)式中：

r为从0′到任一点的距离；

常规高超声速风洞，气体温度在1000K左右，气体仍然可以作为量热完全气体处理，在1000K以上，则要考虑真实气体效应；

量热不完全的真实气体喷管设计方法和理想气体设计方法的步骤相同，但原先量热完全的气体的速度比、密度和马赫数的关系不再能适用；对对常规高超声速风洞高马赫数要利用NACA TR1135的修正量，即同样马赫数下不同总温下真实气体物理量和理想气体物理量的比值，进行修正；

为了对喷管喉道面积做量热不完全气体影响的修正，引入函数F(M_E)对喷管喉道面积做不完全气体影响的修正，引入函数：

则：

到此，只要知道M_A和F(M_E)，喉道区型面曲线TA可由(5)式计算出；

步骤S311、计算F(M_E)：

不完全气体的比热比有：

a_tp ²＝RTγ_tp (16)

a_*tp ²＝RT_*γ_*tp (20)

将(16)～(21)式代入下式：

而

由M_tp＝1和(M_E)_tp为喷管出口设计马赫数代入(19)式分别求出T_*和T_E，代入(22)获得(A/A_*)_tp和(23)式一起代入(14)求得：

对给定的T₀由M_E就可以计算得F(M_E)；

步骤S312、计算M_A：

在轴对称等熵流中，P-M膨胀角有以下关系式：

ψ＝ψ_B-θ (26)

为求ψ_B必须选择M_B的值，由于M_B选择有任意性，有两种方法：

步骤(a)由Cresci的建议选取的M_B值比M_C低0.2即：

M_B＝M_C-0.2 (27)

步骤(b)根据BC线上的速度系数成三次方分布，并利用边界条件最后推得B点与C点的速度系数关系

式中：M_C＝M_E，或W_C＝W_E；

确定了M_B之后，按以下步骤求解M_A：

步骤(c)将M_B(W_B)代入(25)式求出ψ_B；

步骤(d)以ψ_B代入(26)式并将ψ及θ换成A点的ψ_A和θ_A的值，求出ψ_A；

步骤(e)将ψ_A代入(25)式求出M_A；

步骤S32、计算沿径向流波前AB气流参数：

在波前AB上有几何分割关系：

(29)式中：M_A≤M_P≤M_B；

N为做特征线网格时沿径向流波前AB的分割总数，从B到A，P＝0，1，2…N；

由M_P代入(25)、(26)式可求得θ_P值；

由几何关系：

(30)式中：

步骤S33、计算ABCD转换区中心线BC上的M数分布：

BC上的M数分布是连续变化的，假定速度系数成三次方分布：

(31)式中：

系数a₀、a₁、a₂、a₃由边界条件决定，边界条件为：

在B点：

在C点：

由(32)、(33)式解得：

a₀＝W_B ；a₁＝3(W_C-W_B) ；a₂＝-3(W_C-W_B) ；a₃＝W_C-W_B (34)

将(34)式代入(31)式得：

且M和W的关系如下：

此处

当γ＝1.4时，α＝6，则

利用(32)和(35)式求得：

上式对W微分，且x＝x_B得：

最后有：

(40)式中：

步骤S34、计算消波区BCE参数：

利用步骤二中的初始条件，做特征线网格求解型面AD,当特征线网格做到DC，则以DC上的参数作初始条件，只要求得CE上的初始条件，用同样的方法作特征线网格，就可以解出型面DE，于是整个型面AE就确定了；

要使试验段获得均匀的平行气流，DCE为减波区，CE必须为一条直线，在CE上的M数均为M_E且CE上的参数为：

(41)式中，N为做特征线网格时沿径向流波前AB的分割总数，P＝0，1，2…N；

步骤S35、计算轴对称特征线网格：

轴对称喷管中，气流特性是对称于中心线的，可以只研究通过中心轴的xy平面的流态就能决定整个喷管的流态。

这里用轴对称特征线网格计算型面AD和DE的方法，喷管的对称轴为x轴，气流沿x轴是对称的，在物理平面xy上设已知点P₁(x₁，y₁)和P₂(x₂，y₂)，P₁(x₁，y₁)的左特征线和P₂(x₂，y₂)的右特征线交于P₃(x₃，y₃)，如图3所示；

由几何关系：

(42)式中：

(42)式经简化整理得：

y₃＝y₁+(x₃-x₁)tan(θ₁+α₁) (44)

或

y₃＝y₂+(x₃-x₂)tan(θ₂-α₂) (45)

由轴对称特征线方程：

(48)式中：上边符号对应于左特征线，下边符号对应于右特征线，ds_1,2为左右特征线微元长度，把上式改写为点1、2、3的差分方程为：

(47)式中：

由(47)式解得：

(48)式的下标13，23是指取P₁和P₃，P₂和P₃的平均值；

当P1和P2位于X_BC上，则有y₁＝y₂＝θ₁＝θ₂＝0，(48)式右边第二项为不定式，出现奇异点，若在X_BC附近近似为径向流有：

故(47)、(48)式变为

为了对上述各式提高其计算精度，采用迭代法计算。迭代式如下：第一先用P₁、P₂点代入以上各式计算出P₃的值，以后各次可用下列各式子进行迭代计算：

将(50)式代入(48)或(49)，直到达到5×10^-6精度为止；

步骤S36、计算轴对称喷管壁面坐标：

当特征线做到边界上的点A后，为了要确定流线上的点R，必须对边界条件作一些处理，求出壁面上的坐标；利用边界上流线和特征线的关系，如图4所示，图中AR为流线，13和23为1、2点发出的左右特征线；

线性插值：

得：

由流线和特征线的关系：

由(53)解出：

将x_R代入

得：

线性关系：

x_R由(54)式给出；其中，下标1为前一个壁面点，初始壁面点为A点，下一个壁面点是以R为初始条件进行的，如此类推，最后完成AE型面的计算。

在上述技术方案中，所述步骤四轴对称喷管超声速段位流型面附面层进行修正的计算方法包括以下步骤：

用无粘流特征线法计算得到的喷管壁面坐标，理论上应该就是喷管的型面坐标，但实际上沿喷管壁面有边界层增长，愈到下游边界层愈厚。边界层的增长，使得喷管气流有效通道截面积减少，因而喷管出口不能得到设计的马赫数，并且也不能得到满意的均匀流场。为防止这一现象，首先用无粘特征线法求出喷管无粘型面及沿壁面马赫数分布，然后把它作为求边界层位移厚度方程的初始值，在无粘流型面的坐标上，加上边界层位移厚度，最终得到喷管实际型面坐标。采用这种方法是为了防止因边界层的存在而影响试验段马赫数的分布的均匀性和准确性。高超声速风洞喷管壁面边界层一般(低密度风洞除外)为湍流边界层。而且，在风洞中为了防止空气组分液化，需要把空气加热到数百度至数千度，因此，有的风洞喷管壁面还需要冷却(常规高超声速风洞的高马赫数喷管)。由于壁面冷却，在壁面上存在热交换，这样，喷管边界层计算就变得复杂。

对于湍流附面层的近似计算方法，在低速时，有种种方法，而用于总温很高的高超音速的近似计算方法一般有四种，即Sivells-Payne方法，Reshtko-Tucker方法，Persh-Lee方法和Bartz方法。附面层修正计算，即在喷管的每一个横截面上给以一个质量补偿量，它对超音速喷管也同样适用，不同之处在于高超音流动需要考虑温度效应。

对有热交换和压力梯度的高超声速边界层，特别是高焓气体边界层计算方法现在还不太成熟，大多数是把低速湍流边界层计算方法移植到高超声速边界层，用得比较多的是Sivells-Payne方法。本文采用Sivells-Payne方法进行边界层修正。

步骤S41、建立动量方程：

在高超音速喷管中，附面层属于湍流附面层，湍流附面层增长用冯·卡门轴对称动量积分方程来描述：

方程(59)左边最后一项只有在轴对称的情况下才出现，这里Cf、H可由湍流附面层的关系给出；将喷管的曲线坐标s换成轴对称坐标x，如图5所示，得：

方程(59)变为：

利用stwartson变换：

由动量厚度θ定义，并以(61)式代入得：

(62)式中：

由位移厚度δ^*定义，并利用温度型关系式：

其中，

当Pr＝1时，在绝热壁面条件下，壁面总温T_s＝T₀，则δ^* _tr由H定义及(62)、(66)式得：

最后把(63)、(67)代入(60)简化后得：

上列各式中：下标“0”表示自由驻点条件；下标“e”表示附面层底部条件；下标“tr”表示经Stewartson转换后的参数。初始条件：当x＝0时，θ_tr＝0；

为了求方程(68)，必须先求出气流与喷管壁面的摩擦系数C_f和附面层形状因子H_tr：

步骤S411、计算附面层形状因子系数H_tr：

采用修正的Crocco二次定律，可给出考虑了绝热壁温和Pr≠1影响的可压缩流附面层形状因子系数H_tr，附面层内的温度分布为：

其中(69)式中：T_w为壁温，T_aw为绝热壁温；由(61)、(69)式代入(65)式并简化得：

其中，

把(62)式代入(70)式得：

其中，

(67)、(71)联合解得：

只要给出H_i，就可以求出H_tr，然而H_i与M无关，只和不可压摩擦系数有关，由半经验公式给出：

这里系数7是实验值；

步骤S412、计算气流与喷管壁面的摩擦系数C_f：

压缩流比不可压缩流的知识不完善得多，因此，使用所谓参考温度的概念，从而可以利用不可压流的C_fi求出可压流的C_f，所谓参考温度法，就是即使有压缩性存在的场合，在壁面附近也形成不可压缩流，所以把适当的温度作为基准来评价ρ，μ，那么，不可压缩流情况下的关系就可以用于可压缩流，取壁温作为参考温度时，可以知道评价压缩性大小程度；利用Eckert给出的参考温度法，即参考温度T′为：

其中C′_f＝F(Re_x′)使用的是不可压情况下的C_fi＝F(Re_x)关系式；

在此

其中，

由无压力梯度的Karmán—Schoenherr的平均摩擦系数公式：

由此式从Re_xi解C_Fi是不容易的，引入它的近似关系式，及求局部摩擦系数：

把不可压的Re_xi换成Re_x′并代入(75)式得：

为了把不可压的C_fi应用到可压的情况，必须对方程(68)作如下变换：

把(78)、(79)代入上式并作近似计算得：

其中，

将(68)转换成不可压的情况，由Re→Re_x′→Re_X(＝Re_xi)求C_fi，Re_x′→C_f。

以上各式带“'”的参数是以T′为特征温度的值，因为高马赫喷管要进行强迫水冷，所以设计时，对于M_E＝10时，取T_w＝573K；M_E＝11时，取T_w＝623K。绝热壁温按下式计算：

式中：σ为复温系数，对于M_E＝5、6、7取σ＝0.88，对于M_E＝8、9、10、11、12取σ＝0.896，均为湍流状态；

步骤S42、实际计算要考虑的因素；

(1)Re数的范围

得到的马赫数与实Karmán—Schoeherr的公式，只适合没有压力梯度的情形，因此，应注意到(79)式的结果也是针对没有压力梯度的。而且(78)式为近似公式，当logRe_xi＝1.5或logRe_xi＝2.3686时，出现奇异点，以上公式适用在Re＝10⁵～10⁹范围，计算是满足要求的。

(2)假想原点

附面层的计算通常是以喉道为初始条件的，当x＝0时，求得Cf是发散的，为了消除奇异点，以喉部上游为假想原点，则喉部处x≠0，假想原点的x_*由下式求得。

这样原来的x坐标从喉道往上游移动x_*的距离，移轴后的坐标应为x_s＝x_*+x。

(3)附面层修正因子

为了得到预先给定的出口直径，在进行附面层修正时，要引入附面层修正因子f，修正后喷管型面坐标如下：

计算时，第一步利用无粘位流型面参数，计算出喷管出口的

然后通过(84)式计算出f₁，利用f₁进行计算，得到修正后的喷管型面参数Y₁、X₁，在利用Y₁、X₁求出

f₂直到满足y_E-Y_n≤10^-6mm为止。

通过以上步骤，可以得出高超声速喷管的气动型面，为喷管的设计加工提供内型面坐标。

这里说明的设备数量和处理规模是用来简化本发明的说明的。对本发明的应用、修改和变化对本领域的技术人员来说是显而易见的。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。

Claims (5)

1.一种计算高超声速风洞轴对称喷管内型面的方法，其特征在于，包括以下步骤：

步骤一、给定内型面的初始条件参数，利用初始条件参数进行轴对称喷管内型面计算；轴对称喷管内型面计算主要包括位流型面计算和超声速位流型面附面层修正；其中位流型面又包括：亚音速段型面和超声速段型面；流经喷管内型面的位流包括：亚音速区、从喉道到拐点的喉道区和从拐点到喷管出口的转换区；

步骤二、计算轴对称喷管的亚声速型面坐标；

步骤三、计算轴对称喷管的超声速型面坐标，包括：计算喉道区型面参数计算；沿径向流波前气流参数计算；转换区中心线上的马赫数分布计算；消波区参数计算；利用轴对称特征线网格计算型面；计算轴对称喷管壁面坐标；

步骤四、对轴对称喷管的超声速段位流型面附面层进行修正计算。

2.如权利要求1所述的计算高超声速风洞轴对称喷管内型面的方法，其特征在于，所述步骤一中给定的初始条件参数包括：型面马赫数M，喷管入口高度Y_in，喷管出口高度Y_in，喷管长度L，型面马赫角α，喷管壁面总压强P₀，喷管壁面总温度T₀，亚声速段移轴后的喷管进出口半径高度比n，喷管壁温T_w；同时在喷管上引入转换区，标记为ABCD；喉道区型面标记为TA；消波区，标记为BCE。

3.如权利要求2所述的计算高超声速风洞轴对称喷管内型面的方法，其特征在于，所述步骤二中计算轴对称喷管亚声速型面坐标的方法包括以下步骤：

采用维托辛斯基公式计算亚声速端型面，计算按以下公式进行：

上式中，

Y₁′＝Y₁+b；Y_*′＝Y_*+b；Y＝Y′-b；

其中，Y₁、Y_*和Y分别为收缩段的入口、喉道的半高度和X轴上的截面半径，Y₁′、Y_*′和Y′分别为移轴后收缩段的入口、喉道的半高度和X′轴上的截面半径，a为特征长度，b为对X的移轴量、n表示亚音速段移轴后的进出口半径高度比，即n＝Y₁′/Y_*′。

4.如权利要求2所述的计算高超声速风洞轴对称喷管内型面的方法，其特征在于，所述步骤三中计算轴对称喷管的超声速型面坐标的方法包括以下步骤：

步骤S31、计算喉道区型面TA参数；假设喉道区为一元等熵流，经验证明喉道区型面用一元三次方程能满意的获得要求的气流流动状态，该一元三次方程如下：

y＝a₀+a₁x+a₂x²+a₃x³ (2)

系数a₀、a₁、a₂、a₃由边界条件确定，边界条件如下：

x＝0y＝y_*；

将(3)式代入(2)式经整理得：

将(4)式代入(2)式整理后得：

(5)式为所求的喉道区型面曲线，其中θ_A为给定的经验值，x_A由下述关系式确定：

由边界条件得：

x_A＝3(y_A-y_*)(2 tanθ_A)^-1 (6)

由

可知：

于是将(7)、(8)式代入(6)式，求出x_A为：

从连续方程：

将(8)式代入(11)式得：

(12)式中：

r为从0′到任一点的距离；

引入函数F(M_E)对喷管喉道面积做不完全气体影响的修正，引入函数：

则：

到此，只要知道M_A和F(M_E)，喉道区型面曲线TA可由(5)式计算出；

步骤S311、计算F(M_E)：

不完全气体的比热比有：

a_tp ²＝RTγ_tp (16)

a_*tp ²＝RT_*γ_*tp (20)

将(16)～(21)式代入下式：

而

由M_tp＝1和(M_E)_tp为喷管出口设计马赫数代入(19)式分别求出T_*和T_E，代入(22)获得(A/A_*)_tp和(23)式一起代入(14)求得：

对给定的T₀由M_E就可以计算得F(M_E)；

步骤S312、计算A点马赫数M_A：

在轴对称等熵流中，P-M膨胀角有以下关系式：

ψ＝ψ_B-θ (26)

为求ψ_B必须选择M_B的值，由于M_B选择有任意性，有两种方法：

步骤(a)由Cresci的建议选取的M_B值比M_C低0.2即：

M_B＝M_C-0.2 (27)

步骤(b)根据BC线上的速度系数成三次方分布，并利用边界条件最后推得B点与C点的速度系数关系

式中：M_C＝M_E，或W_C＝W_E；

确定了M_B之后，按以下步骤求解M_A：

步骤(c)将M_B(W_B)代入(25)式求出ψ_B；

步骤(d)以ψ_B代入(26)式并将ψ及θ换成A点的ψ_A和θ_A的值，求出ψ_A；

步骤(e)将ψ_A代入(25)式求出M_A；

步骤S32、计算沿径向流波前AB气流参数：

在波前AB上有几何分割关系：

(29)式中：M_A≤M_P≤M_B；

N为做特征线网格时沿径向流波前AB的分割总数，从B到A，P＝0，1，2…N；

由M_P代入(25)、(26)式可求得θ_P值；

由几何关系：

(30)式中：

步骤S33、计算ABCD转换区中心线BC上的M数分布：

BC上的M数分布是连续变化的，假定速度系数成三次方分布：

(31)式中：

系数a₀、a₁、a₂、a₃由边界条件决定，边界条件为：

在B点：

在C点：

由(32)、(33)式解得：

a₀＝W_B；a₁＝3(W_C-W_B)；a₂＝-3(W_C-W_B)；a₃＝W_C-W_B (34)

将(34)式代入(31)式得：

且M和W的关系如下：

此处

当γ＝1.4时，α＝6，则

利用(32)和(35)式求得：

上式对W微分，且x＝x_B得：

最后有：

(40)式中：

步骤S34、计算消波区BCE参数：

利用步骤二中的初始条件，做特征线网格求解型面AD,当特征线网格做到DC，则以DC上的参数作初始条件，只要求得CE上的初始条件，用同样的方法作特征线网格，就可以解出型面DE，于是整个型面AE就确定了；

要使试验段获得均匀的平行气流，DCE为减波区，CE必须为一条直线，在CE上的M数均为M_E且CE上的参数为：

(41)式中，N为做特征线网格时沿径向流波前AB的分割总数，P＝0，1，2…N；

步骤S35、计算轴对称特征线网格：

轴对称喷管中，气流特性是对称于中心线的，可以只研究通过中心轴的xy平面的流态就能决定整个喷管的流态；

这里用轴对称特征线网格计算型面AD和DE的方法，喷管的对称轴为x轴，气流沿x轴是对称的，在物理平面xy上设已知点P₁(x₁，y₁)和P₂(x₂，y₂)，P₁(x₁，y₁)的左特征线和P₂(x₂，y₂)的右特征线交于P₃(x₃，y₃)；

由几何关系：

(42)式中：

(42)式经简化整理得：

y₃＝y₁+(x₃-x₁)tan(θ₁+α₁) (44)

或

y₃＝y₂+(x₃-x₂)tan(θ₂-α₂) (45)

由轴对称特征线方程：

(46)式中：上边符号对应于左特征线，下边符号对应于右特征线，ds_1,2为左右特征线微元长度，把上式改写为点1、2、3的差分方程为：

(47)式中：

由(47)式解得：

(48)式的下标13，23是指取P₁和P₃，P₂和P₃的平均值；

当P1和P2位于X_BC上，则有y₁＝y₂＝θ₁＝θ₂＝0，(48)式右边第二项为不定式，出现奇异点，若在X_BC附近近似为径向流有：

故(47)、(48)式变为

为了提高计算精度，采用迭代法计算，迭代式如下：第一先用P₁、P₂点代入以上各式计算出P₃的值，以后各次可用下列各式子进行迭代计算：

将(50)式代入(48)或(49)，直到达到5×10^-6精度为止；

步骤S36、计算轴对称喷管壁面坐标：

当特征线做到边界上的点A后，为了要确定流线上的点R，必须对边界条件作一些处理，求出壁面上的坐标；利用边界上流线和特征线的关系，设流线为AR，13和23为1、2点发出的左右特征线；

线性插值：

得：

由流线和特征线的关系：

由(53)解出：

将x_R代入

得：

线性关系：

x_R由(54)式给出；其中，下标1为前一个壁面点，初始壁面点为A点，下一个壁面点是以R为初始条件进行的，如此类推，最后完成AE型面的计算。

5.如权利要求2所述的计算高超声速风洞轴对称喷管内型面的方法，其特征在于，所述步骤四轴对称喷管超声速段位流型面附面层进行修正的计算方法包括以下步骤：

步骤S41、建立动量方程：

在高超音速喷管中，附面层属于湍流附面层，湍流附面层增长用冯·卡门轴对称动量积分方程来描述：

方程(59)左边最后一项只有在轴对称的情况下才出现，这里C_f、H可由湍流附面层的关系给出；将喷管的曲线坐标s换成轴对称坐标x，得：

方程(59)变为：

利用stwartson变换：

由动量厚度θ定义，并以(61)式代入得：

(62)式中：

由位移厚度δ^*定义，并利用温度型关系式：

其中，

当Pr＝1时，在绝热壁面条件下，壁面总温T_s＝T₀，则δ^* _tr由H定义及(62)、(66)式得：

最后把(63)、(67)代入(60)简化后得：

上列各式中：下标“0”表示自由驻点条件；下标“e”表示附面层底部条件；下标“tr”表示经Stewartson转换后的参数，初始条件：当x＝0时，θ_tr＝0；

为了求方程(68)，必须先求出气流与喷管壁面的摩擦系数C_f和附面层形状因子H_tr：

步骤S411、计算附面层形状因子系数H_tr：

采用修正的Crocco二次定律，可给出考虑了绝热壁温和Pr≠1影响的可压缩流附面层形状因子系数H_tr，附面层内的温度分布为：

其中(69)式中：T_w为壁温，T_aw为绝热壁温；由(61)、(69)式代入(65)式并简化得：

其中，

把(62)式代入(70)式得：

其中，

(67)、(71)联合解得：

只要给出H_i，就可以求出H_tr，然而H_i与M无关，只和不可压摩擦系数有关，由半经验公式给出：

这里系数7是实验值；

步骤S412、计算气流与喷管壁面的摩擦系数C_f：

利用Eckert给出的参考温度法，即参考温度T′为：

其中C′_f＝F(Re_x′)使用的是不可压情况下的C_fi＝F(Re_x)关系式；

在此

其中：

由无压力梯度的Karmán—Schoenherr的平均摩擦系数公式：

由此式从Re_xi解C_Fi是不容易的，引入它的近似关系式，及求局部摩擦系数：

把不可压的Re_xi换成Re_x′并代入(75)式得：

为了把不可压的C_fi应用到可压的情况，必须对方程(68)作如下变换：

把(78)、(79)代入上式并作近似计算得：

其中，

将(68)转换成不可压的情况，由Re→Re_x′→Re_X(＝Re_xi)求C_fi，Re_x′→C_f；

带“'”的参数是以T′为特征温度的值，因为高马赫喷管要进行强迫水冷，所以设计时，对于M_E＝10时，取T_w＝573K；M_E＝11时，取T_w＝623K，绝热壁温按下式计算：

式中：σ为复温系数，对于M_E＝5、6、7取σ＝0.88，对于M_E＝8、9、10、11、12取σ＝0.896，均为湍流状态；

步骤S42、实际计算要考虑的因素；

(1)Re数的范围

得到的马赫数与实Karmán—Schoeherr的公式，只适合没有压力梯度的情形，因此，应注意到(79)式的结果也是针对没有压力梯度的，而且(78)式为近似公式，当logRe_xi＝1.5或logRe_xi＝2.3686时，出现奇异点，以上公式适用在Re＝10⁵～10⁹范围，计算是满足要求的；

(2)假想原点

附面层的计算通常是以喉道为初始条件的，当x＝0时，求得C_f是发散的，为了消除奇异点，以喉部上游为假想原点，则喉部处x≠0，假想原点的x_*由下式求得：

这样原来的x坐标从喉道往上游移动x_*的距离，移轴后的坐标应为x_s＝x_*+x；

(3)附面层修正因子

为了得到预先给定的出口直径，在进行附面层修正时，要引入附面层修正因子f，修正后喷管型面坐标如下：

计算时，第一步利用无粘位流型面参数，计算出喷管出口的

然后通过(84)式计算出f₁，利用f₁进行计算，得到修正后的喷管型面参数Y₁、X₁，在利用Y₁、X₁求出

f₂直到满足|y_E-Y_n|≤10^-6mm为止；

通过以上步骤，可以得出高超声速喷管的气动型面，为喷管的设计加工提供内型面坐标。

Images