本发明涉及超/高超声速流场的反设计技术领域，尤其涉及基于特征线法的无激波压力最速上升曲线及其对应的无激波压力最速上升流场的求解方法。
背景技术：
特征线方法作为求解非线性双曲型偏微分方程的精确方法，在超声速/高超声速流动中得到了广泛的应用，尤其是吸气式超声速/高超声速飞行器进气道及喷管的设计。
目前基于特征线法的流场求解主要有包括两种类型：给定壁面型线求解流场和给定流场中的某流动参数求解流场及流场对应的壁面。前者称为基于特征线法的流场正设计，后者称为基于特征线法的流场反设计。
基于特征线法的流场反设计通过给定流场参数可求解出相应的流场。这一功能是现行通用流场计算软件fluent所不具备的，因此在流场设计中表现出极大的优势。其中的基于压力分布的流场反设计通过给定壁面的压力分布实现对特定压力分布下的流场及壁面的求解。在具体求解中通过左右行马赫线相交求解下游点的单元过程和右(左)行马赫线与壁面相交求解下游壁面点的单元过程在空间中步进求解，从而得到完整的流场。
基于压力分布的流场反设计方法，需要预先给定压力分布。而在特征线设计方法中，若给定压力分布曲线斜率过大，会出现流场求解中特征线相交的情况，从而出现非物理解。因此在基于压力分布的流场反设计中确定出恰好使特征线不相交所对应的压力最快上升曲线显得很有必要。而目前并没有一种直接可以确定出压力最快上升曲线的方法。目前解决该问题常用的做法是先根据经验给定一条使特征线不相交的压力曲线，然后逐渐增大压力分布曲线，最终获得满足压缩要求的压力分布曲线。这种方法存在两个问题，其一设计效率低下，需要进行多次人工优化；其二，这种方法下得到的压力分布曲线往往并不是最优的。
技术实现要素：
针对现有技术存在的缺陷，本发明提供了一种基于压力的特征线法中无激波压力最速上升曲线的求解方法。
为了实现上述技术目的，本发明的技术方案是：
本发明提供一种基于压力的特征线法中无激波压力最速上升曲线的求解方法，首先通过二分法给定待求解壁面点的压力值，然后根据特征线法中的单元求解过程得到该壁面点的位置坐标和流动参数，并求解由该壁面点发出的马赫波线上的所有点的位置坐标以及流动参数。根据求解得到的马赫线与马赫波初始线的末端点的位置关系，判断初次给定的壁面点的压力值是否满足要求，如不满足要求，再次使用二分法重新确定该壁面点压力值，重复上述过程，直至满足上述要求。最后，采用上述相同的方法求解出全流场内所有壁面点和内部点的位置坐标和流动参数，从而求解得到整个无激波压力最速上升流场，其中无激波压力最速上升流场壁面的压力曲线为无激波壁面压力最速上升曲线。
具体地。一种基于压力的特征线法中无激波压力最速上升曲线的求解方法，包括以下步骤：
s1.给定马赫波初始线a-fp，其中，点a为马赫波初始线的起始点，同时也是无激波压力最速上升流场壁面的起始壁面点，为给定值。点fp为马赫波初始线的末端点，为给定值。
对线a-fp进行均匀离散得到一系列的离散点，离散精度为流场实际尺寸中离散点间距控制在10-2～10-1m量级，得到的离散点为a、点a1、点a2……点fp。离散得到的线a-fp上的各离散点的位置坐标和流动参数均为已知条件，来流条件为已知条件，其中，离散点的坐标(x,y)，流动参数包括离散点的静压p、速度v、密度ρ、马赫数ma、静温t、流动方向角θ。来流条件包括静压p，密度ρ，静温t以及马赫数ma。
s2.基于二分法初次确定第一个待求解的壁面点b0的压力值，其中：待求解的壁面点b0为起始壁面点a的下游壁面点。
s3.根据特征线法中基于压力的特征线壁面点的预估-校正方法，利用起始壁面点a和与起始壁面点a相邻的内部点a1求解壁面点b0的位置和流动参数。
s4.基于有旋特征线的内部点求解方法，利用壁面点b0以及点a2作为上游相邻特征点，求解由点b0发出的马赫波线上的点b1的位置坐标以及流动参数。接着利用点b1以及点a3作为上游相邻特征点，求解由点b0发出的马赫波线上的点b2的位置坐标以及流动参数，利用点b2以及点a4作为上游相邻特征点，求解由点b0发出的马赫波线上的点b3的位置坐标以及流动参数，依次类推，直至求得以点fp为上游点求解得到由点b0发出的马赫波线上的最后一个内部点bn的位置坐标以及流动参数，至此得到从壁面点b0发出的马赫波线上点b1、点b2、点b3……点bn的位置坐标和流动参数。依次光滑连接壁面点b0、点b1、点b2、点b3……点bn所得到的线即是从壁面点b0发出的马赫波线。
s5.通过s4中得到的由壁面点b0发出的马赫波线与马赫波初始线末端点fp的位置关系，判断s2中初次确定的壁面点b0的压力值是否满足要求，如果不满足要求，则返回s2，基于二分法重新给定壁面点b0压力值，并按照s3、s4中的方法重新计算由壁面点b0发出的马赫波线上的点，直至给定的壁面点b0的压力值满足要求为止，输出最终满足要求的壁面点b0的位置坐标和流动参数，以及由壁面点b0发出的马赫波线上所有点的位置坐标和流动参数。
s6.采用s2至s5中相同的方法求解下一个待求解的壁面点的位置坐标和流动参数以及由该壁面点发出的马赫波线上所有点的位置坐标和流动参数。依次类推，直至求解出流场内的所有点，该流场即为无激波压力最速上升流场，对应的流场壁面的压力曲线为无激波压力最速上升曲线。
求解完壁面点b0以及由壁面点b0发出的马赫波线后，对于下一个待求解的壁面点c0，同样的先给定壁面点c0其壁面点压力值的初始取值区间，按照s1中相同的方法初次确定壁面点c0的压力值，接着按照s2中相同的方法，利用壁面点b0和内部点b1求解壁面点c0，接着按照s3中相同的方法，求解从壁面点c0发出的马赫线上的点c1、c2…的位置坐标和流动参数。依次光滑连接壁面点c0、点c1、点c2……所得到的线即是从壁面点c0发出的马赫波线。按照s4中相同的方法对初次确定的壁面点c0的压力值是否满足要求，如果不满足要求，再次基于二分法重新给定壁面点c0压力值，并按照s2、s3中的方法重新计算由壁面点c0发出的马赫波线上的点，直至给定的壁面点c0的压力值满足要求为止。
依次类推，继续采用s2至s5中相同的方法求解下一个待求解的壁面点……，直至解出流场内的所有点。
在本发明s2中，首先给定待求解的壁面点b0其壁面点压力值的初始取值区间
在本发明中由于求解的无激波压力最速上升流场对气流起压缩作用，壁面压力不断增大，因此可选择
pa为给定值，一般取值范围在1×10～1×102量级。
由此，通过式(2)初次确定待求解的壁面点b0的压力值
本发明s3中，基于压力的特征线壁面点的预估-校正方法是本领域公知技术。具体可参见《气体动力学》，童秉纲，孔祥言，邓国华，高等教育出版社，2012年，p227-261；《气体动力学》(下册)，m.j.左克罗，j.d.霍夫曼，国防工业出版社，1984年，p126-136。
壁面点b0的上游壁面点即起始壁面点a，利用点a、点a1求解第一个待求解的壁面点b0的位置坐标以及流动参数(即利用上游壁面点a和与上游壁面点a相邻的内部点a1求解上游壁面点a的下游壁面点b0的位置坐标以及流动参数)，下游壁面点的求解采用特征线法中壁面点的顺处理方法，具体求解采用预估-校正方法。
本发明s4中，点b1、b2、b3……bn的求解方法相同。基于有旋特征线的内部点求解过程为本领域公知技术，具体可参见《气体动力学》，童秉纲，孔祥言，邓国华，高等教育出版社，2012年，p227-261；《气体动力学》(下册)，m.j.左克罗，j.d.霍夫曼，国防工业出版社，1984年，p126-136。
本发明s5中，判断s2中初次确定的壁面点b0的压力值是否满足要求的方法是：
通过式(19)、(20)求壁面点b0和马赫波初始线末端点fp的连线的斜率以及壁面点b0发出马赫波的起始点b0和终止点bn之间连线的斜率，其中壁面点b0的位置坐标(xb0，yb0)，终止点bn的位置坐标(xbn，ybn)，马赫波初始线末端点fp的位置坐标(xfp，yfp)。
km＝(yfp-yb0)/(xfp-xb0)(19)
kt＝(ybn-yb0)/(xbn-xb0)(20)
δ＝km-kt(21)
并通过式(21)对压力值是否合适进行判断，ε为给定值，取值范围设置为10-4～10-3；
若δ＞ε，则说明壁面点b0处的压力值太小，通过式(22)对壁面点b0的压力值进行调整，其中式中参数右上角标代表迭代次数，即将第k次迭代的值按照下述规则赋给第k+1次迭代作为迭代初值；
其中
若δ＜ε，则说明壁面点b0处的压力值太大，通过式(23)对壁面点b0的压力值进行调整，其中式中参数右上角标代表迭代次数，即将第k次迭代的值按照下述规则赋给第k+1次迭代作为迭代初值；
其中
与现有技术相比，本发明具有以下优点：
将特征线法中基于压力分布的流场反设计和二分法结合，通过二分法调整壁面压力值，使该点发出的马赫波末端与马赫波初始线末端点重合，从而求解得到无激波压力最速上升流场及无激波壁面压力最速上升曲线。
采用本发明方法可避免在特征线法基于给定压力分布的流场反设计中因压力分布给定不当而产生的特征线相交的情况；利用本发明可求解得到使特征线恰不相交的压力最速上升曲线。
基于压力分布的流场设计，一般在设计中给定的压力分布曲线多为解析式函数，限制了流场壁面压力分布的多样性且无法预知给定的压力分布是否合理。通过本发明生成的最速压力曲线对现行基于压力分布的特征线流场反设计中压力分布的给定可以起指导作用。
附图说明
图1为本发明无激波压力最速上升流场求解示意图
图中：
1：超/高超声速来流条件，为给定值
点a：马赫波初始线的起始点，同时也是无激波压力最速上升流场壁面的起始壁面点，为给定值
点fp：马赫波初始线的末端点，为给定值；
点a1：马赫波初始线上的离散点，为给定值；
点a2：马赫波初始线上的离散点，为给定值；
点a3：马赫波初始线上的离散点，为给定值；
点a4：马赫波初始线上的离散点，为给定值；
点b0：无激波压力最速上升流场壁面点，待求；
点c0：无激波压力最速上升流场壁面点，待求；
点d0：无激波压力最速上升流场壁面点，待求；
点e0：无激波压力最速上升流场壁面点，待求；
点b1：无激波压力最速上升流场内部点，待求。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。
本发明提供一种基于压力的特征线法中无激波压力最速上升曲线的求解方法，包括以下步骤：
s1.参照图1，给定马赫波初始线a-fp，其中，点a为马赫波初始线的起始点，同时也是无激波压力最速上升流场壁面的起始壁面点，为给定值。点fp为马赫波初始线的末端点，为给定值。
对线a-fp进行均匀离散得到一系列的离散点，离散精度为流场实际尺寸中离散点间距控制在10-2～10-1m量级，在马赫波初始线的起始点与马赫波初始线的末端点之间得到一系列的离散点为a、点a1、点a2……点fp。离散得到的各离散点的位置坐标和流动参数均为已知条件，来流条件为已知条件。其中，离散点的坐标(x,y)，流动参数包括离散点的静压p、速度v、密度ρ、马赫数ma、静温t、流动方向角θ。来流条件包括静压p，密度ρ，静温t以及马赫数ma。
s2.基于二分法初次确定第一个待求解的壁面点b0的压力值，待求解的壁面点b0为起始壁面点a的下游壁面点。
首先给定待求解的壁面点b0其壁面点压力值的初始取值区间
在本发明中由于求解的无激波压力最速上升流场对气流起压缩作用，壁面压力不断增大，因此可选择
pa为给定值，一般取值范围在1×10～1×103量级。
由此，通过式(2)初次确定待求解的壁面点b0的压力值
s3.根据特征线法中基于压力的特征线壁面点的预估-校正方法，求解壁面点b0的位置坐标和流动参数；
基于压力的特征线壁面点的预估-校正方法是本领域公知技术。具体可参见《气体动力学》，童秉纲，孔祥言，邓国华，高等教育出版社，2012年，p227-261；《气体动力学》(下册)，m.j.左克罗，j.d.霍夫曼，国防工业出版社，1984年，p126-136。
壁面点b0的上游壁面点即起始壁面点a，利用点a、点a1求解第一个待求解的壁面点b0的位置坐标以及流动参数(即利用上游壁面点a和与上游壁面点a相邻的内部点a1求解上游壁面点a的下游壁面点b0的位置坐标以及流动参数)，下游壁面点的求解采用特征线法中壁面点的顺处理方法。具体求解采用预估-校正方法，过程如下：
s3.1先对壁面点的位置坐标和流动参数进行预估；
由点a1发出的右行特征线a1-b0和壁面点a发出流线a-b0相交，交点为壁面点b0。壁面点b0的预估位置坐标由式(3)求解得到。
式(3)中，下标n代表该参数为待求解的壁面点，在式中代表壁面点b0的参数，下标m1代表该参数为上游相邻内部点(与上游壁面点相邻的内部点)，在式中代表点a1的参数，下标m0代表该参数为待求解的壁面点的上游壁面点，在式中代表a的参数。式中θ为该点(即θ其下标代表的点，下标m1代表该参数为上游相邻内部点，下标m0代表该参数为上待求解的壁面点的上游壁面点)的流动方向角，μ为该点的马赫角，式(3)中流动方向角μm1可通过式(4)中给定该点(即μm1其下标代表的点，下标m1代表该参数为上游相邻内部点，在式中代表点a1的参数)的速度v和静温t求得。
将s1中初次确定的壁面点b0的压力值代入方程(5)、(6)，即
pn-pm0＝(am0)2(ρn-ρm0)(5)
ρm0vm0(vn-vm0)+pn-pm0＝0(6)
式中下标m0的参数代表该参数为待求解的壁面点的上游壁面点的参数，在式中代表点a的参数，下标n的参数代表待求解的壁面点的参数，在式中代表b0的参数。pn代表为待求解的壁面点b0的压力值，pm0代表待求解的壁面点的上游壁面点a的压力值，am0为壁面点a的当地声速，ρn代表待求解的壁面点b0的密度，ρm0代表待求解的壁面点的上游壁面点a的密度，vn为待求解的壁面点b0的速度，vm0为待求解的壁面点的上游壁面点a的速度，vm1为与上游壁面点a相邻的内部点a1的速度。
再通过求解式(7)，得到壁面点b0的流动方向角θb0＝θn；
至此完成对壁面点b0位置坐标和流动参数的预估。
s3.2对s2.1中得到的壁面点b0的预估位置坐标和预估流动参数进行初次校正；
以
根据上式得到初次校正后的壁面点b0的位置坐标。
分别将式(7)中的差分因子以外的量用平均值代替，得到式(9)并计算，得到初次校正后的壁面点b0的速度、密度和流动方向角
s3.3将得到的初次校正后的壁面点b0的位置坐标、速度、密度和流动方向角按照s3.2中的校正方法迭代校正多次，当迭代至满足式(10)时认为计算收敛，停止计算，输出最后校正得到的壁面点b0位置坐标和流动参数，至此，完成对壁面点b0的求解；
式(10)中εb为给定值，取值范围一般设置为10-4～10-3。式(10)中各参数的参数上标代表校正的次数。
s4基于有旋特征线的内部点求解方法，利用壁面点b0以及点a2作为上游相邻特征点，求解由点b0发出的马赫波线上的点b1的位置坐标以及流动参数；接着利用点b1以及点a3作为上游相邻特征点，求解由点b0发出的马赫波线上的点b2的位置坐标以及流动参数，利用点b2以及点a4作为上游相邻特征点，求解由点b0发出的马赫波线上的点b3的位置坐标以及流动参数，依次类推，直至求得以点fp为上游点求解得到由点b0发出的马赫波线上的最后一个内部点bn的位置坐标以及流动参数，至此得到从壁面点b0发出的马赫波线上点b1、点b2、点b3……点bn的位置坐标和流动参数，依次光滑连接壁面点b0、点b1、点b2、点b3……点bn所得到的线即是从壁面点b0发出的马赫波线；
基于有旋特征线的内部点求解过程为本领域公知技术，具体可参见《气体动力学》，童秉纲，孔祥言，邓国华，高等教育出版社，2012年，p227-261；《气体动力学》(下册)，m.j.左克罗，j.d.霍夫曼，国防工业出版社，1984年，p126-136。
点b1、b2、b3……bn的求解方法相同。下面利用壁面点b0以及点a2作为上游相邻特征点，求解由点b0发出的马赫波线上的点b1的位置坐标以及流动参数为例来进行说明。
s4.1对点b1的位置坐标和流动参数进行预估：
先由式(11)求解点b1的预估位置坐标
yn-ym＝tan(θm±μm)(xn-xm)(11)
其中，m代表上游相邻特征点b0和点a2，n为下游所求特征线上的节点b1，对于左行特征线b0-b1取式(11)中“+”，对于右行特征线a2-b1取式(11)中“-”，分别将点b0、点a2的流动参数和位置参数带入式(12)联列求解，得到点b1的速度
然后求解点b1所在流线与线a2-b0交点，记为
其中kmm为线a2-b0的斜率。
再由公式(14)、(15)得到点b1的密度ρn和压力值pn。
由此，完成对点b1位置参数和流动参数的预估。
s4.2对s4.1中得到的点b1的预估位置坐标和预估流动参数进行校正步；
用
s4.3将得到的初次校正后的点b1的位置坐标、速度、密度和流动方向角按照s4.2中的校正方法迭代校正多次，当迭代至满足式(18)时认为计算收敛，停止计算，输出最后校正得到的点b1位置坐标和流动参数，至此，完成对点b1的求解；
式(18)中εb为给定值，取值范围一般设置为10-4～10-3。式(18)中各参数的参数上标代表校正的次数。
由此得到了点b1的位置坐标和流动参数，
s5.通过s3和s4中得到的由壁面点b0发出的马赫波线与马赫波初始线末端点fp的位置关系，判断s2中初次确定的壁面点b0的压力值是否满足要求，如果不满足要求，则返回s2，基于二分法重新给定壁面点b0压力值，并按照s3、s4中的方法重新计算由壁面点b0发出的马赫波线上的点，直至给定的壁面点b0的压力值满足要求为止，输出最终满足要求的壁面点b0的位置坐标和流动参数，以及由壁面点b0发出的马赫波线上所有点的位置坐标和流动参数。
下面开始判断通过s2中使用二分法确定的壁面点压力值是否合适。此处以判断壁面点b0的压力值为例进行说明。其他壁面点压力值的判断与该过程相同。
通过式(19)、(20)求壁面点b0和马赫波初始线末端点fp的连线的斜率以及壁面点b0发出马赫波的起始点b0和终止点bn之间连线的斜率，其中壁面点b0的位置坐标(xb0，yb0)，终止点bn的位置坐标(xbn，ybn)，马赫波初始线末端点fp的位置坐标(xfp，yfp)。
km＝(yfp-yb0)/(xfp-xb0)(19)
kt＝(ybn-yb0)/(xbn-xb0)(20)
δ＝km-kt(21)
并通过式(21)对压力值是否合适进行判断。ε为给定值，取值范围设置为10-4～10-3；
若δ＞ε，则说明壁面点b0处的压力值太小，通过式(22)对壁面点b0的压力值进行调整，其中式中参数右上角标代表迭代次数，即将第k次迭代的值按照下述规则赋给第k+1次迭代作为迭代初值；
其中
若δ＜ε，则说明壁面点b0处的压力值太大，通过式(23)对壁面点b0的压力值进行调整，其中式中参数右上角标代表迭代次数，即将第k次迭代的值按照下述规则赋给第k+1次迭代作为迭代初值；
其中
s5.采用s2至s5中相同的方法求解下一个待求解的壁面点的位置坐标和流动参数，以及由该壁面点发出的马赫波线上所有点的位置坐标和流动参数；依次类推，直至求解出流场内的所有点，该流场即为无激波压力最速上升流场，对应的流场壁面的压力曲线为无激波压力最速上升曲线。
求解完壁面点b0以及由壁面点b0发出的马赫波线后，对于下一个待求解的壁面点c0，同样的先给定壁面点c0其壁面点压力值的初始取值区间，按照s1中相同的方法初次确定壁面点c0的压力值，接着按照s2中相同的方法，利用壁面点b0和内部点b1求解壁面点c0，接着按照s3中相同的方法，求解从壁面点c0发出的马赫线上的点c1、c2…的位置坐标和流动参数。依次光滑连接壁面点c0、点c1、点c2……所得到的线即是从壁面点c0发出的马赫波线。按照s4中相同的方法对初次确定的壁面点c0的压力值是否满足要求，如果不满足要求，再次基于二分法重新给定壁面点c0压力值，并按照s2、s3中的方法重新计算由壁面点c0发出的马赫波线上的点，直至给定的壁面点c0的压力值满足要求为止。
求解完壁面点c0以及由壁面点c0发出的马赫波线后，对于下一个待求解的壁面点d0，同样的先给定壁面点d0其壁面点压力值的初始取值区间，按照s1中相同的方法初次确定壁面点d0的压力值，接着按照s2中相同的方法，利用壁面点c0和内部点c1求解壁面点d0，接着按照s3中相同的方法，求解从壁面点d0发出的马赫线，按照s4中相同的方法对初次确定的壁面点d0的压力值是否满足要求，如果不满足要求，再次基于二分法重新给定壁面点d0压力值，并按照s2、s3中的方法重新计算由壁面点d0发出的马赫波线上的点，直至给定的壁面点d0的压力值满足要求为止。
依次类推，继续求解下一个待求解的壁面点……，直至解出流场内的所有点。
以上所述仅是本发明的优选实施方式，本发明的保护范围并不局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。