CN115618168A - 高马赫数高焓风洞中皮托压数据处理方法 - Google Patents

Abstract

本发明公开了一种高马赫数高焓风洞中皮托压数据处理方法，包括步骤：S1、输入风洞来流条件和皮托测点数据；S2、调用临界状态处理流程，判断皮托测点是否处于超声速流动，若是超声速流动则进入步骤S3，否则进入步骤S4；S3、调用超声速环境下的皮托压数据处理流程，得到测点当地的状态参数；S4、调用非超声速环境下的总压数据处理流程，得到测点当地的状态参数；S5、根据步骤S3得到的测点当地的状态参数或步骤S4得到的测点当地的状态参数为已知条件，得到流场其他信息。本发明解决了真实气体效应环境下的流场试验数据处理问题，对于单点的流场参数测量、多测点情况下的边界层参数测量及截面参数测量具有重要意义。

Description

高马赫数高焓风洞中皮托压数据处理方法

技术领域

本发明涉及实验空气动力学技术领域，具体涉及马赫数达到M10量级、焓值达到

量级的高马赫数高焓风洞中皮托压数据处理方法。

背景技术

随着先进航天飞行器的快速发展,与航天飞行器紧密相关的高超声速气动力和气动加热问题越发受到高度重视，当来流总温在2500K以上时（如大气中马赫数10飞行时的总温接近4000K，中国空气动力研究与发展中心的FD-14A激波风洞总温可达3500K以上，中科院力学所的JF系列激波风洞来流总温高达8000K），气流尤其是与壁面接触处的气流存在高温真实气体效应。

高马赫数进气道内流试验研究中常常需要测量皮托压（如测量截面的流量或边界层参数分布等），由于高温真实气体效应的存在，传统的完全气体假设不再成立，在进行数据处理时必须考虑比热随温度等的变化效应。一个感兴趣的基本问题是：如何在高马赫数高焓风洞试验数据处理中，通过壁面静压p ₁、皮托总压P _t ^*测量值求解得到对应测点处的流状态参数，如当地马赫数、密度、速度、总压等。

关于高焓真实气体效应方面的研究，目前主要集中在CFD气体模型、外流气动力/热研究试验、单变量隐式方程求解等方面：

(1)在CFD数值模拟领域，针对气体模型开展的研究

北京交通大学的胡雨濛基于化学非平衡模型,引入了对于能量方程的变换方法和等效比热比的概念,开发了一套高超声速条件下化学平衡流和化学非平衡流气体的数值模拟程序,为近空间高超声速气动热分析提供了稳健的数值模拟工具；

南京理工大学的代春良等人为研究高马赫数内转进气道在真实气体效应下的工作特性，设计了额定工作状态Ma12的高超声速内转进气道，结合不同气体模型对其进行了数值模拟。

(2)在外流气动力与气动热试验领域

中国科学院力学研究所的刘云峰等人采用六分量应变天平在JF12长实验时间激波风洞上开展的10°尖锥标模的天平测力实验研究，研究了真实气体效应对该模型气动力特性的影响；

西北工业大学的王小朋等人以钝化锥导乘波体为研究对象,开展了高焓激波风洞测热试验以及高温化学非平衡气动加热数值验证,对乘波布局滑翔飞行器前缘线和下壁面热流分布特征进行了研究。

(3)工业控制领域隐式方程求解方面的研究工作

中国移动集团福建有限公司福州分公司的江磊等人开展了工业控制领域隐式方程在线求解器的设计研究，以单一变量隐式方程作为研究对象，利用新型组态软件PineCAD作为仿真编程平台，实现对该方程的在线求解。从数学方程本身特性以及反馈控制原理两方面入手，设计出两种各具特色的求解模型，解决了工业现场控制以及仿真研究中普遍存在的隐式方程在线求解的难题。

通过调研发现，到目前为止尚未见真实气体效应环境下风洞试验皮托压数据处理技术的报道。

现有技术中，无法根据皮托压数据处理得到皮托测点的状态参数。因为：

高马赫数高焓风洞试验中，来流高温超过2500K将会产生高温真实气体效应，传统的完全气体假设不再成立。在高焓2500K及以上环境下，定压比热C _p 已不再是常数，而是温度的函数（如多项式函数），于是比热比γ不再是常数(如对空气γ=1.4，对燃气γ=1.33等)，而是与温度相关的一个表达式。相应地，焓值是Cp与温度的积分表达式，熵值则是静压、C _p 以及温度的积分表达式。这些表达式的复杂性使得类似f(x₁、x₂)=y₀的隐式方程组不像完全气体条件下γ是常数时的单一方程那样易于求解。

发明内容

为了解决上述问题，本发明旨在提供一种高马赫数高焓风洞中皮托压数据处理方法，以来流条件和皮托测点数据为入口，针对高马赫数高焓风洞试验中皮托测点数据满足五变量五方程且其中两个方程是隐式方程的特征，调用变比热等熵过程和变比热正激波过程，以求解五变量方程组为研究对象，通过合理的预估-迭代校正数值方法使得隐式方程显式化，按照方程组解耦求解与耦合迭代的思路，建立皮托压数据处理方法，据此可处理得到测点处气流应考虑真实气体效应而更为准确更符合实际的波前信息。

为达到上述目的，本发明具有如下技术方案：

一种高马赫数高焓风洞中皮托压数据处理方法，其特征在于，包括如下步骤：

S1、输入风洞来流条件和皮托测点数据，作为步骤S2-S4的入口参数；

其中，S1中的风洞来流条件和皮托测点数据包括风洞来流总温T ₀*、测点所在截面的壁面静压p ₁、测点皮托压P _t ^*；

S2、调用临界状态处理流程，判断皮托测点是否处于超声速流动，若是超声速流动则进入步骤S3，否则进入步骤S4；

其中，S2所述的临界状态处理流程，包括步骤S21-S24：

S21、设定Ma ₁=1，然后根据公式(1)预估T ₁：

(1)

其中，γ为常数，T ₁ 、Ma ₁分别为当地静温、当地马赫数，T ₁*为当地总温，T ₁*＝T ₀*，τ（Ma ₁）是关于Ma ₁的τ函数，γ为比热比，对空气γ=1.4；

S22、以步骤S21中所述的Ma ₁、T ₁以及当地静压P ₁为入口，调用变比热等熵过程求得T ₁*即T _cr ^*，其中P ₁=p ₁，T _cr ^*为测点处气流在声速状态时对应的总温即临界总温；采用二分法迭代校正T ₁，以Ma ₁与校正的T ₁以及当地静压P ₁为入口重新求解T _cr ^*，迭代校正T ₁直到T _cr ^*= T ₀*时完毕，得到T ₁校正值和T _cr ^*；

S23、以步骤S22中所述的T ₁校正值、T _cr ^*以及当地静压P ₁为入口，调用变比热等熵过程求得P ₁*即P _cr ^*，其中P ₁*为当地总压，P _cr ^*为测点处气流在声速状态时对应的总压即临界总压；

S24、判断P _t ^*与P _cr ^*的大小，若P _t ^* >P _cr ^*则进入步骤S3，否则进入步骤S4；

S3、调用超声速环境下的皮托压数据处理流程，得到测点当地的状态参数；

其中，S3所述的超声速环境下的皮托压数据处理流程，包括步骤S31-S34：

S31、采用定γ法，根据公式(2)预估Ma ₁，根据步骤S21中所述的公式(1)预估T ₁：

(2)

其中，f(Ma ₁ )是关于Ma ₁的函数，P ₂*为波后总压，P ₂*= P _t ^*，π(Ma ₁ )、σ(Ma ₁ )分别是等熵气动函数和正激波总压恢复系数；

S32、以步骤S31中所述的T ₁、Ma ₁及当地静压P ₁为入口，调用变比热等熵过程求解对应的T ₁*、P ₁*，固定T ₁采用二分法迭代校正Ma ₁，迭代校正Ma ₁直到T ₁*= T ₀*为止，得到Ma ₁校正值；

S33、以步骤S32中所述的T ₁、Ma ₁校正值及当地静压P ₁为入口，调用变比热正激波过程求得波后静压P ₂、波后静温T ₂、波后马赫数Ma ₂；

S34、以步骤S33中所述的P ₂、T ₂、Ma ₂为入口，调用变比热等熵过程求解P ₂*、T ₂*，当P ₂ *≠P _t ^* 时采用二分法迭代校正T ₁，返回步骤S32-S34进行迭代，迭代校正T ₁直到满足P ₂ *=P _t ^* ，得到测点当地的状态参数P ₁、T ₁、Ma ₁，其中T ₂*为波后总温；

S4、调用非超声速环境下的总压数据处理流程，得到测点当地的状态参数；

其中，S4所述的非超声速环境下的总压数据处理流程，包括步骤S41-S43：

S41、采用定γ法，根据公式(3)预估Ma ₁，根据步骤S21中所述的公式(1)预估T ₁：

(3)；

S42、以步骤S41中所述的T ₁、Ma ₁及当地静压P ₁为入口，调用变比热等熵过程求解对应的P ₁*、T ₁*，包括：

S4201、固定T ₁采用二分法迭代校正Ma ₁，迭代校正Ma ₁直到T ₁ *= T ₀*为止，得到T ₁对应的Ma ₁与P ₁*；

S4202、采用二分法迭代校正T ₁，根据步骤S41中所述的公式(3)预估Ma ₁，返回步骤S4201，迭代校正T ₁直到P ₁ *=P _t ^* 为止；

S43、得到测点当地的状态参数P ₁、T ₁、Ma ₁；

S5、根据步骤S3得到的测点当地的状态参数或步骤S4得到的测点当地的状态参数为已知条件，得到流场的其他信息。

作为优选方式，变比热等熵过程包括步骤：

变比热等熵过程包括步骤：

(11)、根据已知点参数P ₁ ^(E)、T ₁ ^(E)、Ma ₁ ^(E)求得V ₁ ^(E)，其中P ₁ ^(E)、T ₁ ^(E)、Ma ₁ ^(E)、V ₁ ^(E)分别是等熵过程已知点的当地静压、当地静温、当地马赫数、当地气流速度；

(12)、定γ法预估等熵过程待求点的当地静温T ₂ ^(E)；

(13)、根据能量方程(4)求得V ₂ ^(E)，并求得Ma ₂ ^(E)：

(4)

其中，Ma ₂ ^(E)、V ₂ ^(E)分别为等熵过程待求点的当地马赫数、当地气流速度，C _P (T)为与温度T相关的定压比热，dT为温度T的微分；

(14)、根据熵方程(5)且等熵时ds=0，对熵函数两侧积分得到公式(6)，根据公式(6)求得P ₂ ^(E)：

且ds=0 (5)

(6)

其中，P ₂ ^(E)为待求点的当地静压，ds为熵的微分，R为气体常数，dP为压力P的微分；

(15)、判断P ₂ ^(E)或Ma ₂ ^(E) 是否满足迭代结束条件，如满足迭代结束条件则进入步骤(16)，否则采用二分法迭代校正T ₂ ^(E)，返回步骤(13)-步骤(15)进行迭代，迭代校正T ₂ ^(E)直到P ₂ ^(E) 或 Ma ₂ ^(E)满足迭代结束条件为止；

(16)、输出等熵过程待求点状态参数P ₂ ^(E)、T ₂ ^(E)、Ma ₂ ^(E)；

(17)、将步骤(15)中所述的迭代结束条件换成Ma ₂ ^(E)=0，重复步骤(13)-步骤(16)求得当地总压P ₁* 、当地总温T ₁*；

(18)、将步骤(15)中所述的迭代结束条件换成Ma ₂ ^(E)=1，重复步骤(13)-步骤(16)求得临界状态压力P _cr 、临界状态温度T _cr 。

作为优选方式，变比热正激波过程包括步骤：

(21)、根据已知波前参数P ₁ ^(S)、T ₁ ^(S)、Ma ₁ ^(S)求得V ₁ ^(S)，其中P ₁ ^(S)、T ₁ ^(S)、Ma ₁ ^(S)、V ₁ ^(S)分别是波前静压、波前静温、波前马赫数、波前气流速度；

(22)、定γ法预估波后静温T ₂ ^{(S )}；

(23)、根据能量方程(7)求得 V ₂ ^{(S )}，并求得Ma ₂ ^(S)：

(7)

其中，Ma ₂ ^(S)、V ₂ ^{(S )}分别为波后马赫数、波后气流速度；

(24)、根据基本方程(8)结合状态方程(9)可得到连续方程(10)，通过连续方程(10)求得 P ₂ ^(S)并记为P ₂ ^{(S )(1)}：

(8)

(9)

(10)

其中，ρ ₁ ^(S)、ρ ₂ ^(S)分别为波前气流密度和波后气流密度，ρ为气流密度，P ₂ ^{(S )}为波后静压， P ₂ ^{(S )(1)}为采用连续方程(10)得到的波后静压；

(25)、根据基本方程(11)结合状态方程(9)可得到动量方程(12)，通过动量方程(12)求得P ₂ ^(S)并记为P ₂ ^{(S )(2)}：

(11)

(12)

其中，P ₂ ^{(S )(2)}为采用动量方程(12)得到的波后静压；

(26)、根据公式(13)得到T ₂ ^{(S )}对应的静压差∆P：

(13)

(27)、判断∆P=0是否成立，若成立则进入步骤(28)，否则采用二分法迭代校正T ₂ ^{(S )}，返回步骤(23)-步骤(27)进行迭代，迭代校正T ₂ ^{(S )}直到满足∆P=0为止；

(28)、输出波后参数P ₂ ^(S)、T ₂ ^(S)、Ma ₂ ^(S)。

作为优选方式，步骤S5中，得到流场的其他信息，包括步骤：

S51、根据步骤S34中得到的测点当地的状态参数或步骤S43中得到的测点当地的状态参数P ₁、T ₁、Ma ₁，由T ₁求得T ₁ 时的定压比热Cp ₁和比热比γ₁，由 T ₁、γ₁求得 T ₁时的声速a₁，由a₁、Ma ₁求得当地的气流速度V₁，由P ₁、T ₁并结合状态方程求得当地的气流密度ρ₁；

S52、流场的其他信息还包括步骤S3中得到的当地总压 P ₁*、当地总温T ₁*或步骤S4中得到的当地总压P ₁*、当地总温T ₁* 。

相比现有技术，本发明的有益效果：

(1)本发明不是气体模型的计算与对比，研究对象也不是外流气动力/热，方程的求解既不是利用工业控制单元，也不是求解单一隐式变量方程；

(2)本发明针对高马赫数高焓风洞试验中皮托测点数据满足五变量五方程且其中两个方程是隐式方程的特征，调用变比热等熵过程和变比热正激波过程，以求解五变量方程组为研究对象，通过合理的预估-迭代校正数值方法使得隐式方程显式化；

(3)根据绝热过程气流总焓不变的原理，可知皮托测点当地的总温、波后总温、临界总温等均与风洞总温相等。使得风洞总温T ₀*必须作为已知量参与到数据处理中，测点气流各个状态的总温必须与风洞总温一致，实际上给出的是一个等熵过程的能量方程，而等熵过程能量方程的引入实际上又给出了测点当地气流某任意状态即波前状态或波后状态的温度与马赫数之间存在着相互联系的定量关系；

(4)本发明技术中，未知待求的状态参数虽然只有2个，分别为当地静温T ₁与当地马赫数Ma ₁，但实际求解时，由于出现了状态参数以外的已知量即空间测点的皮托压亦即波后总压P _t ^*，因此需要引入了正激波过程，而正激波过程又产生了额外的未知参数即波后的3个状态参数P ₂、T ₂、Ma ₂，这就使得未知参数变为5个即T ₁、Ma ₁，P ₂、T ₂、Ma ₂，由连续方程、动量方程、变比热能量方程3个正激波关系式，外加波前等熵关系式、波后等熵关系式组成的方程组共有5个方程，因此方程是封闭的，从而解也是唯一的；

(5)变比热等熵过程主要调用了能量方程与等熵方程，但都考虑了变比热的影响。对于已知某一点状态参数的等熵过程，在其他任意一点的状态参数中，压力、温度、马赫数只要知道其一，就可得到其余两个：如果已知温度，变比热等熵过程最为简单；如果已知压力或马赫数，则需要采用二分法校正温度，然后返回已知温度处迭代，直到压力或马赫数满足要求为止，这正是将隐式方程显式化的关键举措；

(6)变比热正激波过程求解连续方程与动量方程时没有耦合求解而是解耦求解并将两个压力结果进行耦合作差，根据耦合结果进行迭代求解的思路符合计算机数值求解的实际，同时也达到了能够快速求解复杂方程组的一个较好效果；

(7)本发明技术解决了真实气体效应环境下的流场试验数据处理问题，对于单点的流场参数测量、多测点情况下的边界层参数测量及截面参数测量具有重要意义。

附图说明

图1是本发明实施例提供的高马赫数高焓激波风洞中皮托压数据处理流程主框图。

图2是本发明实施例提供的临界状态处理流程框图。

图3是本发明实施例提供的超声速环境下的皮托压数据处理流程框图。

图4是本发明实施例提供的非超声速环境下的总压数据处理流程框图。

图5是本发明实施例提供的变比热等熵过程框图。

图6是本发明实施例提供的变比热正激波过程框图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细的说明：

如图1所示的，一种高马赫数高焓风洞中皮托压数据处理方法，包括如下步骤：

S1、输入风洞来流条件和皮托测点数据，作为步骤S2-S4的入口参数；

S2、调用临界状态处理流程，判断皮托测点是否处于超声速流动，若是超声速流动则进入步骤S3，否则进入步骤S4；

S3、调用超声速环境下的皮托压数据处理流程，得到测点当地的状态参数；

S4、调用非超声速环境下的总压数据处理流程，得到测点当地的状态参数；

S5、根据步骤S3得到的测点当地的状态参数或步骤S4得到的测点当地的状态参数为已知条件，得到流场的其他信息。

读入风洞来流参数：主要是风洞总温T ₀*、空间皮托测点对应的壁面静压p ₁、空间皮托测点测得的皮托压P _t ^*，然后调用临界状态处理流程获得P _cr ^*，并用测得的皮托压与之比对来确定皮托测点是否处于超声速流动，并据此调用超声速环境下的皮托压数据处理流程或非超声速环境下的总压处理流程，求得测点当地的静压T ₁与马赫数Ma ₁。

其中涉及的三个主要流程分别为临界状态处理流程、超声速环境下的皮托压数据处理流程、非超声速环境下的总压数据处理流程。这三个主要流程还调用了两个基本过程，即变比热等熵过程、变比热正激波过程。

判断测点处气流是否为超声速流动的判据是P _t ^*＞P _cr ^*还是P _t ^*≤P _cr ^*，前者对应超声速流动，后者对应非超声速流动。之所以区分皮托测点所处的状态，就是因为超声速时必须既调用变比热等熵过程又调用变比热正激波过程，而非超声速状态不能调用变比热正激波过程，只需调用变比热等熵过程即可。

求得测点当地的状态参数P ₁、T ₁、Ma ₁后，可以由此求得皮托测点的其他辅助气流信息。

S1的风洞来流条件和皮托测点数据包括：风洞来流总温T ₀*、测点所在截面的壁面静压p ₁、测点皮托压P _t ^*。

具体的，风洞来流条件和皮托测点数据为：

来流条件：需要用到的是来流总温T ₀*；

皮托测点数据包括：空间测点对应的壁面静压p ₁、空间测点皮托压P _t ^*。

如图2所示的，S2的临界状态处理流程，包括：

S21、设定Ma ₁=1，然后根据公式(1)预估T ₁：

(1)

其中，γ为常数 T ₁ 、Ma ₁分别为当地静温、当地马赫数，T ₁*为当地总温，T ₁*＝T ₀*，τ (Ma ₁ )是关于Ma ₁的τ函数，γ为比热比，对空气γ=1.4；

S22、以步骤S21中的Ma ₁、T ₁以及当地静压P ₁为入口，调用变比热等熵过程求得T ₁*即T _cr ^*，其中P ₁=p ₁，T _cr ^*为测点处气流在声速状态时对应的总温即临界总温；采用二分法迭代校正T ₁，以Ma ₁与校正的T ₁以及当地静压P ₁为入口重新求解T _cr ^*，迭代校正T ₁直到T _cr ^*= T ₀*时完毕，得到T ₁校正值和T _cr ^*；

S23、以步骤S22中的T ₁校正值、T _cr ^*以及当地静压P ₁为入口，调用变比热等熵过程求得P ₁*即P _cr ^*，其中P ₁*为当地总压，P _cr ^*为测点处气流在声速状态时对应的总压即临界总压；

S24、判断P _t ^*与P _cr ^*的大小，若P _t ^* >P _cr ^*则进入步骤S3，否则进入步骤S4。

具体的，首先该流程的入口参数就是流程主框图的入口参数：T ₀*、p ₁、P _t ^*，出口参数是用于流程主框图中控制流程走向的P _cr ^*。

①设定Ma ₁=1，然后根据公式(1)预估T ₁：

(1)

其中：γ为常数， T ₁*＝T ₀*，τ（Ma ₁）是关于Ma ₁的τ函数。显然，将Ma ₁=1和定值γ以及T ₁*代入公式(1)求解出T ₁。以Ma ₁=1与预估得到的T ₁作为入口进入下一步。

②以Ma ₁=1、T ₁以及当地静压P ₁为入口，其中P ₁=p ₁，调用变比热等熵过程求得T ₁*，记为T _cr ^*。由于入口的T ₁是通过定γ法预估得到的，此时求得的T ₁*必然与风洞总温T ₀*不相等，因此，必须迭代校正T ₁，校正T ₁的数学方法很多，这里选择二分法。迭代是指，以Ma ₁=1、校正的T ₁以及当地静压P ₁为入口重新求解T _cr ^*；迭代校正T ₁直到T _cr ^*= T ₀*时完毕。转入下一步。

③以p ₁、T ₁、T _cr ^*为入口，利用变比热等熵过程求得P _cr ^*。

至此，临界状态处理流程结束，得到的P _cr ^*用来与P _t ^*比较以确定调用超声速还是非超声速处理流程。

临界状态处理流程中待求量为P _cr ^*，需要使用等熵过程的熵方程求解，为此增加了一个变量T ₁，从而还需要使用等熵过程的能量方程。这样，由熵方程和能量方程这两个方程解得两个未知数，也就得到了P _cr ^*。

如图3所示的，S3的超声速环境下的皮托压数据处理流程，包括：

S31、采用定γ法，根据公式(2)预估Ma ₁，根据步骤S21中的公式(1)预估T ₁：

(2)

其中，f(Ma ₁ )是关于Ma ₁的函数，P ₂*为波后总压，P ₂*= P _t ^*，π(Ma ₁ )、σ(Ma ₁ )分别是等熵气动函数和正激波总压恢复系数；

S32、以步骤S31中的T ₁、Ma ₁及当地静压P ₁为入口，调用变比热等熵过程求解对应的T ₁*、P ₁*，固定T ₁采用二分法迭代校正Ma ₁，迭代校正Ma ₁直到T ₁*= T ₀*为止，得到Ma ₁校正值；

S33、以步骤S32中的T ₁、Ma ₁校正值及当地静压P ₁为入口，调用变比热正激波过程求得波后静压P ₂、波后静温T ₂、波后马赫数Ma ₂；

S34、以步骤S33中的P ₂、T ₂、Ma ₂为入口，调用变比热等熵过程求解P ₂*、T ₂*，当P ₂ *≠P _t ^* 时采用二分法迭代校正T ₁，返回步骤S32-S34进行迭代，迭代校正T ₁直到满足P ₂ *=P _t ^* ，得到测点当地的状态参数P ₁、T ₁、Ma ₁，其中T ₂*为波后总温。

具体的，确定皮托测点处于超声速流动后，首先预估皮托测点当地Ma ₁与T ₁；然后固定T ₁，调用变比热等熵过程并迭代校正Ma ₁，得到T ₁对应的Ma ₁；再调用变比热正激波过程，获得波后参数P ₂、T ₂、Ma ₂；以P ₂、T ₂、Ma ₂为入口，再次调用变比热等熵过程，求得波后总压，并与皮托测点的皮托压相比，如果不相等就对T ₁进行二分法迭代校正，然后返回到固定T ₁的地方，循环执行；对T ₁的迭代校正一直进行到波后总压与皮托总压相等才停止；这就求得了皮托测点的波前状态参数。据此还可根据波前状态参数求得测点位置的总压P ₁*、速度V₁、密度ρ₁等其他信息。

超声速环境下的皮托压数据处理流程中待求量为T ₁与Ma ₁，由于涉及到正激波波后总压（即测点皮托压），因此未知量还包括波后P ₂、T ₂、Ma ₂，因此有5个未知量；为此，需要使用正激波过程的连续方程、动量方程、能量方程，以及波前等熵过程的能量方程、波后等熵过程的熵方程，这样方程的数目也是5个，于是方程封闭，解唯一。

详细过程解释如下：

(1)采用定γ法预估波前静温T ₁和波前马赫数Ma ₁

首先预估波前马赫数Ma ₁。γ为定值时，根据公式(2)求得Ma ₁：

(2)

其中：P ₁=p ₁，P ₂*= P _t ^*，π(Ma ₁ )、σ(Ma ₁ )分别是等熵气动函数和正激波总压恢复系数，当比热比γ为定值(如对空气γ=1.4)时，通过公式(2)，利用单自变量方程求解方法容易求解得到Ma ₁。

其次，预估波前静温T ₁。将Ma ₁、定值γ、T ₁*= T ₀*代入步骤S21中的公式(1)可以求解出T ₁。

(2)调用变比热等熵过程求解对应的T ₁*，迭代校正Ma ₁，直到T ₁*= T ₀*为止

由使用定γ预估的入口参数T ₁、Ma ₁通过变比热等熵过程可求得当地总压P ₁*、当地总温T ₁*，但这时必定无法满足T ₁*= T ₀*，因此迭代校正T ₁、Ma ₁尤为必要：

固定T ₁，调用变比热等熵过程求解总压P ₁*及总温T ₁*，采用二分法迭代校正Ma ₁，迭代校正Ma ₁直至T ₁*= T ₀*为止，这样就得到了T ₁对应的Ma ₁。

(3)以P ₁=p ₁、T ₁、Ma ₁的校正值为入口，调用变比热正激波过程，求得波后参数P ₂、T ₂、Ma ₂。

(4)以P ₂、T ₂、Ma ₂为入口，调用变比热等熵过程，求解P ₂*、T ₂*；当P ₂ *≠P _t ^* 时采用二分法迭代校正T ₁，返回步骤(2)迭代，直到满足P ₂ *=P _t ^* 。

(5)得到测点当地的状态参数P ₁、T ₁、Ma ₁，其中P ₁=p ₁。

如图4所示的，S4的非超声速环境下的总压数据处理流程，包括：

S41、采用定γ法，根据公式(3)预估Ma ₁，根据步骤S21中的公式(1)预估T ₁：

(3)；

S42、以步骤S41中的T ₁、Ma ₁及当地静压P ₁为入口，调用变比热等熵过程求解对应的P ₁*、T ₁*，包括：

S4201、固定T ₁采用二分法迭代校正Ma ₁，迭代校正Ma ₁直到T ₁ *= T ₀*为止，得到T ₁对应的Ma ₁与P ₁*；

S4202、采用二分法迭代校正T ₁，根据步骤S41中的公式(3)预估Ma ₁，返回步骤S4201，迭代校正T ₁直到P ₁ *=P _t ^* 为止；

S43、得到测点当地的状态参数P ₁、T ₁、Ma ₁。

具体的，确定皮托测点处于非超声速流动后，①预估皮托测点当地Ma ₁，这可以通过定γ时皮托压公式预估，然后再通过定γ预估得到T ₁；②固定T ₁，迭代校正Ma ₁，调用变比热等熵过程求得测点的满足T ₁ *＝T ₀ ^* 的T ₁、Ma ₁；③迭代校正T ₁，定γ法预估Ma ₁，返回步骤②，直到P ₁ *=P _t ^* 为止，得到了满足T ₁ *＝T ₀ ^* 、P ₁ *=P _t ^* 的T ₁、Ma ₁。

非超声速环境下的总压数据处理流程中的待求量为 T ₁、Ma ₁，需要使用的方程也有两个，分别为等熵过程的能量方程与等熵过程的熵方程。同样方程组是封闭的，其解也唯一且可求。

具体说明如下：

(1) 采用定γ法根据公式(3)预估Ma ₁：

(3)

然后根据步骤S21中的公式(1)预估波前静温T ₁；

(2) 调用变比热等熵过程求解P ₁ *、T ₁ *，迭代校正T ₁、Ma ₁直到满足基本状态参数T ₁ *＝T ₀ ^* 、P ₁ *=P _t ^* 为止

①固定T ₁，调用变比热等熵过程求解P ₁ *、T ₁ *，采用二分法迭代校正Ma ₁，直到T ₁ *=T ₀ ^* 为止，得到T ₁对应的Ma ₁与P ₁ *；

②当P ₁ *≠P _t ^* 时采用二分法迭代校正T ₁，并用定γ预估Ma ₁，返回步骤①，直到P ₁ *=P _t ^* 为止，这就得到了满足T ₁ *＝T ₀ ^* 、P ₁ *=P _t ^* 的T ₁、Ma ₁；

(3)得到测点当地的状态参数P ₁、T ₁、Ma ₁。

如图5所示的，一种高马赫数高焓风洞中皮托压数据处理方法的变比热等熵过程，包括：

(11)、根据已知点参数P ₁ ^(E)、T ₁ ^(E)、Ma ₁ ^(E)求得V ₁ ^(E)，其中P ₁ ^(E)、T ₁ ^(E)、Ma ₁ ^(E)、V ₁ ^(E)分别是等熵过程已知点的当地静压、当地静温、当地马赫数、当地气流速度；

具体的，根据已知点参数P ₁ ^(E)、T ₁ ^(E)、Ma ₁ ^(E)，由T ₁ ^(E)求得T ₁ ^(E)时的比热比γ₁ ^(E)，由T ₁ ^(E)、γ₁ ^(E)求得T ₁ ^(E)时的声速a ₁ ^(E)，由a ₁ ^(E)、Ma ₁ ^(E)求得当地的气流速度V ₁ ^(E)；

(12)、定γ法预估等熵过程待求点的当地静温T ₂ ^(E)；

(13)、根据能量方程(4)求得V ₂ ^(E)，并求得Ma ₂ ^(E)：

(4)

其中，Ma ₂ ^(E)、V ₂ ^(E)分别为等熵过程待求点的当地马赫数、当地气流速度，C _P (T)为与温度T相关的定压比热，dT为温度T的微分；

具体的，根据能量方程(4)求出V ₂ ^(E)后，由预估值T ₂ ^(E)求得T ₂ ^(E)时的比热比γ₂ ^(E)，由T ₂ ^(E)、γ₂ ^(E)求得T ₂ ^(E)时的声速a ₂ ^(E)，由a ₂ ^(E)、V ₂ ^(E)求得当地的气流速度Ma ₂ ^(E)；

(14)、根据熵方程(5)且等熵时ds=0，对熵函数两侧积分得到公式(6)，根据公式(6)求得P ₂ ^(E)：

且ds=0 (5)

(6)

其中，P ₂ ^(E)为待求点的当地静压，ds为熵的微分，R为气体常数，dP为压力P的微分；

(15)、判断P ₂ ^(E)或Ma ₂ ^(E) 是否满足迭代结束条件，如满足迭代结束条件则进入步骤(16)，否则采用二分法迭代校正T ₂ ^(E)，返回步骤(13)-步骤(15)进行迭代，迭代校正T ₂ ^(E)直到P ₂ ^(E) 或 Ma ₂ ^(E)满足迭代结束条件为止；

(16)、输出等熵过程待求点状态参数P ₂ ^(E)、T ₂ ^(E)、Ma ₂ ^(E)；

(17)、将步骤(15)中所述的迭代结束条件换成Ma ₂ ^(E)=0，重复步骤(13)-步骤(16)求得当地总压P ₁* 、当地总温T ₁*；

(18)、将步骤(15)中所述的迭代结束条件换成Ma ₂ ^(E)=1，重复步骤(13)-步骤(16)求得临界状态压力P _cr 、临界状态温度T _cr 。

具体的，变比热等熵过程主要调用了能量方程与等熵方程，但都考虑了变比热的影响。对于已知某一点基本状态参数的等熵过程，在其他任意一点的状态参数中，压力、温度、马赫数只要知道其一，就可得到其余两个：如果已知温度，变比热等熵过程最为简单；如果已知压力或马赫数，则需要校正温度（采用二分法），然后返回已知温度处迭代，直到压力或马赫数满足要求为止。这正是将隐式方程显式化的关键举措。

变比热等熵过程还可求得滞止状态参数P ₁*、T ₁*或临界状态参数P _cr 、T _cr ，只需将迭代结束的条件换成Ma ₂ ^(E)=0或Ma ₂ ^(E)=1即可。

如图6所示的，一种高马赫数高焓风洞中皮托压数据处理方法的变比热正激波过程，包括：

(21)、根据已知波前参数P ₁ ^(S)、T ₁ ^(S)、Ma ₁ ^(S)求得V ₁ ^(S)，其中P ₁ ^(S)、T ₁ ^(S)、Ma ₁ ^(S)、V ₁ ^(S)分别是波前静压、波前静温、波前马赫数、波前气流速度；

具体的，根据已知点参数P ₁ ^(S)、T ₁ ^(S)、Ma ₁ ^(S)，由T ₁ ^(S)求得T ₁ ^(S)时的比热比γ₁ ^(S)，由T ₁ ^(S)、γ₁ ^(S)求得T ₁ ^(S)时的声速a ₁ ^(S)，由a ₁ ^(S)、Ma ₁ ^(S)求得当地的气流速度V ₁ ^(S)；

(22)、定γ法预估波后静温T ₂ ^{(S )}；

(23)、根据能量方程(7)求得 V ₂ ^{(S )}，并求得Ma ₂ ^(S)：

(7)

其中，Ma ₂ ^(S)、V ₂ ^{(S )}分别为波后马赫数、波后气流速度；

具体的，根据能量方程(7)求出V ₂ ^(S)后，由预估值T ₂ ^(S)求得T ₂ ^(S)时的比热比γ₂ ^(S)，由T ₂ ^(S)、γ₂ ^(S)求得T ₂ ^(S)时的声速a ₂ ^(S)，由a ₂ ^(S)、V ₂ ^(S)求得当地的气流速度Ma ₂ ^(S)；

(24)、根据基本方程(8)结合状态方程(9)可得到连续方程(10)，通过连续方程(10)求得 P ₂ ^(S)并记为P ₂ ^{(S )(1)}：

(8)

(9)

(10)

其中，ρ ₁ ^(S)、ρ ₂ ^(S)分别为波前气流密度和波后气流密度，ρ为气流密度，P ₂ ^{(S )}为波后静压， P ₂ ^{(S )(1)}为采用连续方程(10)得到的波后静压；

(25)、根据基本方程(11)结合状态方程(9)可得到动量方程(12)，通过动量方程(12)求得P ₂ ^(S)并记为P ₂ ^{(S )(2)}：

(11)

(12)

其中，P ₂ ^{(S )(2)}为采用动量方程(12)得到的波后静压；

(26)、根据公式(13)得到T ₂ ^{(S )}对应的静压差∆P：

(13)

(27)、判断∆P=0是否成立，若成立则进入步骤(28)，否则采用二分法迭代校正T ₂ ^{(S )}，返回步骤(23)-步骤(27)进行迭代，迭代校正T ₂ ^{(S )}直到满足∆P=0为止；

(28)、输出波后参数P ₂ ^(S)、T ₂ ^(S)、Ma ₂ ^(S)。

具体的，变比热正激波过程首先预估了波后静温，然后调用了变比热能量方程来求波后速度或波后马赫数，再调用变比热连续方程与变比热动量方程分别求得波后静压并作差，得到波后静压差，当静压差不为0时，通过二分法校正波后静温，返回到调用变比热能量方程的地方进行迭代，直到得到的静压差为0为止。变比热过程中除预估时的初值使用了定γ外，其余的求解过程均考虑了变比热的影响。

上述在求解连续方程与动量方程时没有耦合求解而是解耦求解并将两个压力结果进行耦合作差，根据耦合结果进行迭代求解的思路符合计算机数值求解的实际，同时也达到了能够快速求解复杂方程组的一个较好效果。

步骤S5中，得到流场的其他信息，包括：

S51、根据步骤S34中得到的测点当地的状态参数或步骤S43中得到的测点当地的状态参数P ₁、T ₁、Ma ₁，由T ₁求得T ₁ 时的定压比热Cp ₁和比热比γ₁，由 T ₁、γ₁求得 T ₁时的声速a₁，由a₁、Ma ₁求得当地的气流速度V₁，由P ₁、T ₁并结合状态方程求得当地的气流密度ρ₁；

S52、流场的其他信息还包括步骤S3中得到的当地总压 P ₁*、当地总温T ₁*或步骤S4中得到的当地总压P ₁*、当地总温T ₁* 。

具体的，根据测点当地的波前参数得到流场的其他信息：

①通过超声速环境下的皮托压数据处理流程，已得到测点处的波前信息P ₁、T ₁、Ma ₁，其中P ₁=p ₁；

②由T ₁求得Cp ₁、γ₁，由T ₁、γ ₁求得a₁，由a₁、Ma ₁求得V₁，由P ₁、T ₁并结合状态方程求得ρ₁；

③超声速环境下的皮托压数据处理流程也同时求得了波前的总压 P ₁*、总温T ₁* ，很显然T ₁*＝T ₀*。

根据测点的基本状态参数得到流场的其他信息：

①通过非超声速环境下的总压数据处理流程，已得到测点的基本状态参数P ₁、T ₁、Ma ₁，其中P ₁直接采用测量值p ₁；

②由T ₁求得Cp ₁、γ₁，由T ₁、γ ₁求得a₁，由a₁、Ma ₁求得V₁，由P ₁、T ₁并结合状态方程求得ρ₁；

③非超声速环境下的总压数据处理流程也同时求得了测点处的气流总压P ₁*、总温T ₁*。

以上所述的仅是本发明的实施例，方案中公知的具体结构和/或特性等常识在此未作过多描述。应当指出，对于本领域的技术人员来说，在不脱离本发明技术方案的前提下，还可以作出若干变形和改进，这些也应该视为本发明的保护范围，这些都不会影响本发明实施的效果和专利的实用性。本申请要求的保护范围应当以其权利要求的内容为准，说明书的具体实施方式等记载可以用于解释权利要求的内容。

Claims (4)

1.一种高马赫数高焓风洞中皮托压数据处理方法，其特征在于，包括如下步骤：

S1、输入风洞来流条件和皮托测点数据，作为步骤S2-S4的入口参数；

其中，S1中的风洞来流条件和皮托测点数据包括风洞来流总温T ₀*、测点所在截面的壁面静压p ₁、测点皮托压P _t ^*；

S2、调用临界状态处理流程，判断皮托测点是否处于超声速流动，若是超声速流动则进入步骤S3，否则进入步骤S4；

其中，S2所述的临界状态处理流程，包括步骤S21-S24：

S21、设定Ma ₁=1，然后根据公式(1)预估T ₁：

(1)

其中，γ为常数，T ₁ 、Ma ₁分别为当地静温、当地马赫数，T ₁*为当地总温，T ₁*＝T ₀*，τ（Ma ₁）是关于Ma ₁的τ函数，γ为比热比，对空气γ=1.4；

S22、以步骤S21中所述的Ma ₁、T ₁以及当地静压P ₁为入口，调用变比热等熵过程求得T ₁*即T _cr ^*，其中P ₁=p ₁，T _cr ^*为测点处气流在声速状态时对应的总温即临界总温；采用二分法迭代校正T ₁，以Ma ₁与校正的T ₁以及当地静压P ₁为入口重新求解T _cr ^*，迭代校正T ₁直到T _cr ^*= T ₀*时完毕，得到T ₁校正值和T _cr ^*；

S23、以步骤S22中所述的T ₁校正值、T _cr ^*以及当地静压P ₁为入口，调用变比热等熵过程求得P ₁*即P _cr ^*，其中P ₁*为当地总压，P _cr ^*为测点处气流在声速状态时对应的总压即临界总压；

S24、判断P _t ^*与P _cr ^*的大小，若P _t ^* >P _cr ^*则进入步骤S3，否则进入步骤S4；

S3、调用超声速环境下的皮托压数据处理流程，得到测点当地的状态参数；

其中，S3所述的超声速环境下的皮托压数据处理流程，包括步骤S31-S34：

S31、采用定γ法，根据公式(2)预估Ma ₁，根据步骤S21中所述的公式(1)预估T ₁：

(2)

其中，f(Ma ₁ )是关于Ma ₁的函数，P ₂*为波后总压，P ₂*= P _t ^*，π(Ma ₁ )、σ(Ma ₁ )分别是等熵气动函数和正激波总压恢复系数；

S32、以步骤S31中所述的T ₁、Ma ₁及当地静压P ₁为入口，调用变比热等熵过程求解对应的T ₁*、P ₁*，固定T ₁采用二分法迭代校正Ma ₁，迭代校正Ma ₁直到T ₁*= T ₀*为止，得到Ma ₁校正值；

S33、以步骤S32中所述的T ₁、Ma ₁校正值及当地静压P ₁为入口，调用变比热正激波过程求得波后静压P ₂、波后静温T ₂、波后马赫数Ma ₂；

S34、以步骤S33中所述的P ₂、T ₂、Ma ₂为入口，调用变比热等熵过程求解P ₂*、T ₂*，当P ₂ *≠P _t ^* 时采用二分法迭代校正T ₁，返回步骤S32-S34进行迭代，迭代校正T ₁直到满足P ₂ *=P _t ^* ，得到测点当地的状态参数P ₁、T ₁、Ma ₁，其中T ₂*为波后总温；

S4、调用非超声速环境下的总压数据处理流程，得到测点当地的状态参数；

其中，S4所述的非超声速环境下的总压数据处理流程，包括步骤S41-S43：

S41、采用定γ法，根据公式(3)预估Ma ₁，根据步骤S21中所述的公式(1)预估T ₁：

(3)；

S42、以步骤S41中所述的T ₁、Ma ₁及当地静压P ₁为入口，调用变比热等熵过程求解对应的P ₁*、T ₁*，包括：

S4201、固定T ₁采用二分法迭代校正Ma ₁，迭代校正Ma ₁直到T ₁ *= T ₀*为止，得到T ₁对应的Ma ₁与P ₁*；

S4202、采用二分法迭代校正T ₁，根据步骤S41中所述的公式(3)预估Ma ₁，返回步骤S4201，迭代校正T ₁直到P ₁ *=P _t ^* 为止；

S43、得到测点当地的状态参数P ₁、T ₁、Ma ₁；

S5、根据步骤S3得到的测点当地的状态参数或步骤S4得到的测点当地的状态参数为已知条件，得到流场的其他信息。

2.根据权利要求1所述的一种高马赫数高焓风洞中皮托压数据处理方法，其特征在于：变比热等熵过程包括步骤：

(11)、根据已知点参数P ₁ ^(E)、T ₁ ^(E)、Ma ₁ ^(E)求得V ₁ ^(E)，其中P ₁ ^(E)、T ₁ ^(E)、Ma ₁ ^(E)、V ₁ ^(E)分别是等熵过程已知点的当地静压、当地静温、当地马赫数、当地气流速度；

(12)、定γ法预估等熵过程待求点的当地静温T ₂ ^(E)；

(13)、根据能量方程(4)求得V ₂ ^(E)，并求得Ma ₂ ^(E)：

(4)

其中，Ma ₂ ^(E)、V ₂ ^(E)分别为等熵过程待求点的当地马赫数、当地气流速度，C _P (T)为与温度T相关的定压比热，dT为温度T的微分；

(14)、根据熵方程(5)且等熵时ds=0，对熵函数两侧积分得到公式(6)，根据公式(6)求得P ₂ ^(E)：

且ds=0 (5)

(6)

其中，P ₂ ^(E)为待求点的当地静压，ds为熵的微分，R为气体常数，dP为压力P的微分；

(15)、判断P ₂ ^(E)或Ma ₂ ^(E) 是否满足迭代结束条件，如满足迭代结束条件则进入步骤(16)，否则采用二分法迭代校正T ₂ ^(E)，返回步骤(13)-步骤(15)进行迭代，迭代校正T ₂ ^(E)直到P ₂ ^(E) 或 Ma ₂ ^(E)满足迭代结束条件为止；

(16)、输出等熵过程待求点状态参数P ₂ ^(E)、T ₂ ^(E)、Ma ₂ ^(E)；

(17)、将步骤(15)中所述的迭代结束条件换成Ma ₂ ^(E)=0，重复步骤(13)-步骤(16)求得当地总压P ₁* 、当地总温T ₁*；

(18)、将步骤(15)中所述的迭代结束条件换成Ma ₂ ^(E)=1，重复步骤(13)-步骤(16)求得临界状态压力P _cr 、临界状态温度T _cr 。

3.根据权利要求1所述的一种高马赫数高焓风洞中皮托压数据处理方法，其特征在于：变比热正激波过程包括步骤：

(21)、根据已知波前参数P ₁ ^(S)、T ₁ ^(S)、Ma ₁ ^(S)求得V ₁ ^(S)，其中P ₁ ^(S)、T ₁ ^(S)、Ma ₁ ^(S)、V ₁ ^(S)分别是波前静压、波前静温、波前马赫数、波前气流速度；

(22)、定γ法预估波后静温T ₂ ^{(S )}；

(23)、根据能量方程(7)求得 V ₂ ^{(S )}，并求得Ma ₂ ^(S)：

(7)

其中，Ma ₂ ^(S)、V ₂ ^{(S )}分别为波后马赫数、波后气流速度；

(24)、根据基本方程(8)结合状态方程(9)可得到连续方程(10)，通过连续方程(10)求得 P ₂ ^(S)并记为P ₂ ^{(S )(1)}：

(8)

(9)

(10)

其中，ρ ₁ ^(S)、ρ ₂ ^(S)分别为波前气流密度和波后气流密度，ρ为气流密度，P ₂ ^{(S )}为波后静压，P ₂ ^{(S )(1)}为采用连续方程(10)得到的波后静压；

(25)、根据基本方程(11)结合状态方程(9)可得到动量方程(12)，通过动量方程(12)求得P ₂ ^(S)并记为P ₂ ^{(S )(2)}：

(11)

(12)

其中，P ₂ ^{(S )(2)}为采用动量方程(12)得到的波后静压；

(26)、根据公式(13)得到T ₂ ^{(S )}对应的静压差∆P：

(13)

(27)、判断∆P=0是否成立，若成立则进入步骤(28)，否则采用二分法迭代校正T ₂ ^{(S )}，返回步骤(23)-步骤(27)进行迭代，迭代校正T ₂ ^{(S )}直到满足∆P=0为止；

(28)、输出波后参数P ₂ ^(S)、T ₂ ^(S)、Ma ₂ ^(S)。

4.根据权利要求1所述的一种高马赫数高焓风洞中皮托压数据处理方法，其特征在于：步骤S5中，得到流场的其他信息，包括步骤S51-S52：

S51、根据步骤S34中得到的测点当地的状态参数或步骤S43中得到的测点当地的状态参数P ₁、T ₁、Ma ₁，由T ₁求得T ₁ 时的定压比热Cp ₁和比热比γ₁，由 T ₁、γ₁求得 T ₁时的声速a₁，由a₁、Ma ₁求得当地的气流速度V₁，由P ₁、T ₁并结合状态方程求得当地的气流密度ρ₁；

S52、流场的其他信息还包括步骤S3中得到的当地总压 P ₁*、当地总温T ₁*或步骤S4中得到的当地总压P ₁*、当地总温T ₁* 。

Images