申请/专利权人：中国空气动力研究与发展中心高速空气动力研究所

申请日：2020-08-04

公开（公告）日：2023-05-26

公开（公告）号：CN111859520B

主分类号：G06F30/13

分类号：G06F30/13;G06F30/17;G06F30/23;G06F30/28;G06F113/08;G06F119/14

优先权：

专利状态码：有效-授权

法律状态：2023.05.26#授权;2020.11.17#实质审查的生效;2020.10.30#公开

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

主权项：1.一种计算高超声速风洞轴对称喷管内型面的方法，其特征在于，包括以下步骤：步骤一、给定内型面的初始条件参数，利用初始条件参数进行轴对称喷管内型面计算；轴对称喷管内型面计算包括位流型面计算和超声速位流型面附面层修正；其中位流型面又包括：亚音速段型面和超声速段型面；流经喷管内型面的位流包括：亚音速区、从喉道到拐点的喉道区和从拐点到喷管出口的转换区；步骤二、计算轴对称喷管的亚声速型面坐标；步骤三、计算轴对称喷管的超声速型面坐标，包括：计算喉道区型面参数计算；沿径向流波前气流参数计算；转换区中心线上的马赫数分布计算；消波区参数计算；利用轴对称特征线网格计算型面；计算轴对称喷管壁面坐标；步骤四、对轴对称喷管的超声速段位流型面附面层进行修正计算；所述步骤一中给定的初始条件参数包括：型面马赫数M，喷管入口高度Y_in，喷管出口高度Y_in，喷管长度L，型面马赫角α，喷管壁面总压强P0，喷管壁面总温度T0，亚声速段移轴后的喷管进出口半径高度比n，喷管壁温Tw；同时在喷管上引入转换区，标记为ABCD；喉道区型面标记为TA；消波区，标记为BCE；所述步骤二中计算轴对称喷管亚声速型面坐标的方法包括以下步骤：采用维托辛斯基公式计算亚声速端型面，计算按以下公式进行： 上式中，Y1′＝Y1+b；Y*′＝Y*+b；Y＝Y′-b； 其中，Y1、Y*和Y分别为收缩段的入口、喉道的半高度和X轴上的截面半径，Y1′、Y*′和Y′分别为移轴后收缩段的入口、喉道的半高度和X′轴上的截面半径，a为特征长度，b为对X的移轴量、n表示亚音速段移轴后的进出口半径高度比，即n＝Y1′Y*′；所述步骤三中计算轴对称喷管的超声速型面坐标的方法包括以下步骤：步骤S31、计算喉道区型面TA参数；假设喉道区为一元等熵流，经验证明喉道区型面用一元三次方程能满意的获得要求的气流流动状态，该一元三次方程如下：y＝a0+a1x+a2x2+a3x32系数a0、a1、a2、a3由边界条件确定，边界条件如下： 将3式代入2式经整理得： 将4式代入2式整理后得： 5式为所求的喉道区型面曲线，其中θA为给定的经验值，xA由下述关系式确定：由边界条件得：xA＝3yA-y*2tanθA-16由可知： 于是将7、8式代入6式，求出xA为： 从连续方程： 将8式代入11式得： 12式中：r为从0′到任一点的距离；引入函数FME对喷管喉道面积做不完全气体影响的修正，引入函数： 则： 到此，只要知道MA和FME，喉道区型面曲线TA可由5式计算出；步骤S311、计算FME：不完全气体的比热比有：atp2＝RTγtp16 a*tp2＝RT*γ*tp20 将16～21式代入下式： 而 由Mtp＝1和MEtp为喷管出口设计马赫数代入19式分别求出T*和TE，代入22获得AA*tp和23式一起代入14求得： 对给定的T0由ME就可以计算得FME；步骤S312、计算A点马赫数MA：在轴对称等熵流中，P-M膨胀角有以下关系式： ψ＝ψB-θ26为求ψB必须选择MB的值，由于MB选择有任意性，有两种方法：步骤a由Cresci的建议选取的MB值比MC低0.2即：MB＝MC-0.227步骤b根据BC线上的速度系数成三次方分布，并利用边界条件最后推得B点与C点的速度系数关系 式中：MC＝ME，或WC＝WE；确定了MB之后，按以下步骤求解MA：步骤c将MBWB代入25式求出ψB；步骤d以ψB代入26式并将ψ及θ换成A点的ψA和θA的值，求出ψA；步骤e将ψA代入25式求出MA；步骤S32、计算沿径向流波前AB气流参数：在波前AB上有几何分割关系： 29式中：MA≤MP≤MB；N为做特征线网格时沿径向流波前AB的分割总数，从B到A，P＝0，1，2…N；由MP代入25、26式可求得θP值；由几何关系： 30式中：步骤S33、计算ABCD转换区中心线BC上的M数分布：BC上的M数分布是连续变化的，假定速度系数成三次方分布： 31式中：系数a0、a1、a2、a3由边界条件决定，边界条件为：在B点： 在C点： 由32、33式解得：a0＝WB；a1＝3WC-WB；a2＝-3WC-WB；a3＝WC-WB34将34式代入31式得： 且M和W的关系如下： 此处当γ＝1.4时，α＝6，则 利用32和35式求得： 上式对W微分，且x＝xB得： 最后有： 40式中：步骤S34、计算消波区BCE参数：利用步骤二中的初始条件，做特征线网格求解型面AD,当特征线网格做到DC，则以DC上的参数作初始条件，只要求得CE上的初始条件，用同样的方法作特征线网格，就可以解出型面DE，于是整个型面AE就确定了；要使试验段获得均匀的平行气流，DCE为减波区，CE必须为一条直线，在CE上的M数均为ME且CE上的参数为： 41式中，N为做特征线网格时沿径向流波前AB的分割总数，P＝0，1，2…N；步骤S35、计算轴对称特征线网格：轴对称喷管中，气流特性是对称于中心线的，可以只研究通过中心轴的xy平面的流态就能决定整个喷管的流态；这里用轴对称特征线网格计算型面AD和DE的方法，喷管的对称轴为x轴，气流沿x轴是对称的，在物理平面xy上设已知点P1x1，y1和P2x2，y2，P1x1，y1的左特征线和P2x2，y2的右特征线交于P3x3，y3；由几何关系： 42式中： 42式经简化整理得： y3＝y1+x3-x1tanθ1+α144或y3＝y2+x3-x2tanθ2-α245由轴对称特征线方程： 46式中：上边符号对应于左特征线，下边符号对应于右特征线，ds1,2为左右特征线微元长度，把上式改写为点1、2、3的差分方程为： 47式中：由47式解得： 48式的下标13，23是指取P1和P3，P2和P3的平均值；当P1和P2位于XBC上，则有y1＝y2＝θ1＝θ2＝0，48式右边第二项为不定式，出现奇异点，若在XBC附近近似为径向流有： 故47、48式变为 为了提高计算精度，采用迭代法计算，迭代式如下：第一先用P1、P2点代入以上各式计算出P3的值，以后各次可用下列各式子进行迭代计算： 将50式代入48或49，直到达到5×10-6精度为止；步骤S36、计算轴对称喷管壁面坐标：当特征线做到边界上的点A后，为了要确定流线上的点R，必须对边界条件作一些处理，求出壁面上的坐标；利用边界上流线和特征线的关系，设流线为AR，13和23为1、2点发出的左右特征线；线性插值： 得： 由流线和特征线的关系： 由53解出： 将xR代入得： 线性关系： xR由54式给出；其中，下标1为前一个壁面点，初始壁面点为A点，下一个壁面点是以R为初始条件进行的，如此类推，最后完成AE型面的计算；所述步骤四轴对称喷管超声速段位流型面附面层进行修正的计算方法包括以下步骤：步骤S41、建立动量方程：在高超音速喷管中，附面层属于湍流附面层，湍流附面层增长用冯·卡门轴对称动量积分方程来描述： 方程59左边最后一项只有在轴对称的情况下才出现，这里Cf、H可由湍流附面层的关系给出；将喷管的曲线坐标s换成轴对称坐标x，得：方程59变为： 利用stwartson变换： 由动量厚度θ定义，并以61式代入得： 62式中： 由位移厚度δ*定义，并利用温度型关系式： 其中，当Pr＝1时，在绝热壁面条件下，壁面总温Ts＝T0，则δ*tr由H定义及62、66式得： 最后把63、67代入60简化后得： 上列各式中：下标“0”表示自由驻点条件；下标“e”表示附面层底部条件；下标“tr”表示经Stewartson转换后的参数，初始条件：当x＝0时，θtr＝0；为了求方程68，必须先求出气流与喷管壁面的摩擦系数Cf和附面层形状因子Htr：步骤S411、计算附面层形状因子系数Htr：采用修正的Crocco二次定律，可给出考虑了绝热壁温和Pr≠1影响的可压缩流附面层形状因子系数Htr，附面层内的温度分布为： 其中69式中：Tw为壁温，Taw为绝热壁温；由61、69式代入65式并简化得： 其中，把62式代入70式得： 其中，67、71联合解得： 只要给出Hi，就可以求出Htr，然而Hi与M无关，只和不可压摩擦系数有关，由半经验公式给出： 这里系数7是实验值；步骤S412、计算气流与喷管壁面的摩擦系数Cf：利用Eckert给出的参考温度法，即参考温度T′为： 其中C′f＝FRex′使用的是不可压情况下的Cfi＝FRex关系式；在此 其中：由无压力梯度的Karmán—Schoenherr的平均摩擦系数公式： 由此式从Rexi解CFi是不容易的，引入它的近似关系式，及求局部摩擦系数： 把不可压的Rexi换成Rex′并代入75式得： 为了把不可压的Cfi应用到可压的情况，必须对方程68作如下变换： 把78、79代入上式并作近似计算得： 其中，将68转换成不可压的情况，由Re→Rex′→ReX＝Rexi求Cfi，Rex′→Cf；带“'”的参数是以T′为特征温度的值，因为高马赫喷管要进行强迫水冷，所以设计时，对于ME＝10时，取Tw＝573K；ME＝11时，取Tw＝623K，绝热壁温按下式计算： 式中：σ为复温系数，对于ME＝5、6、7取σ＝0.88，对于ME＝8、9、10、11、12取σ＝0.896，均为湍流状态；步骤S42、实际计算要考虑的因素；1Re数的范围得到的马赫数与实Karmán—Schoeherr的公式，只适合没有压力梯度的情形，因此，应注意到79式的结果也是针对没有压力梯度的，而且78式为近似公式，当logRexi＝1.5或logRexi＝2.3686时，出现奇异点，以上公式适用在Re＝105～109范围，计算是满足要求的；2假想原点附面层的计算通常是以喉道为初始条件的，当x＝0时，求得Cf是发散的，为了消除奇异点，以喉部上游为假想原点，则喉部处x≠0，假想原点的x*由下式求得： 这样原来的x坐标从喉道往上游移动x*的距离，移轴后的坐标应为xs＝x*+x；3附面层修正因子为了得到预先给定的出口直径，在进行附面层修正时，要引入附面层修正因子f，修正后喷管型面坐标如下： 计算时，第一步利用无粘位流型面参数，计算出喷管出口的然后通过84式计算出f1，利用f1进行计算，得到修正后的喷管型面参数Y1、X1，在利用Y1、X1求出f2直到满足|yE-Yn|≤10-6mm为止；通过以上步骤，可以得出高超声速喷管的气动型面，为喷管的设计加工提供内型面坐标。

全文数据：

权利要求：

百度查询： 中国空气动力研究与发展中心高速空气动力研究所 一种计算高超声速风洞轴对称喷管内型面的方法

免责声明
1、本报告根据公开、合法渠道获得相关数据和信息，力求客观、公正，但并不保证数据的最终完整性和准确性。
2、报告中的分析和结论仅反映本公司于发布本报告当日的职业理解，仅供参考使用，不能作为本公司承担任何法律责任的依据或者凭证。

阅读全文 双屏查看 官方信息 专利公告 收藏专利 下载PDF 下载WORD