[发明专利]一种计算高超声速风洞轴对称喷管内型面的方法

权利要求书

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

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

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

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

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

所述步骤一中给定的初始条件参数包括：型面马赫数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₃由边界条件确定，边界条件如下：

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

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

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

由边界条件得：

x_A＝3(y_A-y_*)(2tanθ_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型面的计算；

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

步骤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为止；

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