一种基于声速解和特征线逆推的喷管设计方法与流程

1.本发明涉及风洞设计领域，具体地，涉及一种基于声速解和特征线逆推的喷管设计方法。


背景技术：

2.现有半柔壁喷管型面设计中均采用sivells设计方法，在喷管流场均匀性和跨超声速试验段共用性上还存在提高的空间。
3.第一，现有sivells设计方法在喉道处采用泉流假设，认为气流沿径向呈等加速运动，忽略了喷管壁面对流动的影响，造成喷管消波区的初始条件与实际情况存在较大偏差，进而导致喷管膨胀段消波不彻底，影响了马赫数均匀性；第二，现有喷管设计方法在半消波区和消波区之间过渡未考虑喷管马赫数轴线分布，使得流场在试验段内会形成明显的菱形过渡区。该过渡区不仅影响流场均匀性，也导致试验段流场难以适应跨、超声速状态共用需求。


技术实现要素：

4.为解决现有sivells设计方法存在的缺陷，本发明提供了一种基于声速解和特征线逆推的喷管设计方法。
5.本方法具体包括：一种基于声速解和特征线逆推的喷管设计方法，所述方法包括：步骤1：生成多个马赫数下的喷管初始型面；步骤2：基于喷管初始型面，获得喷管的关键控制点信息，基于关键控制点信息获得喷管的特征线关系图；步骤3：基于所述特征线关系图和喷管喉道声速近似解，计算获得喷管喉部声速线参数信息和喷管固壁边界；步骤4：基于喷管喉部声速线参数信息和喷管固壁边界进行特征线推进，获得喷管初始膨胀区的下游边界信息；步骤5：基于喷管初始膨胀区的下游边界信息进行特征线推进，获得喷管转折推进区的下游边界信息；步骤6：基于贝塞尔曲线获得喷管转折推进区轴线终点至喷管特征线逆推区轴线终点之间的马赫数分布信息，基于所述马赫数分布信息和所述喷管转折推进区的下游边界信息，获得喷管最大转折点至喷管特征线逆推区型面终点之间的第一喷管型面坐标信息，以及获得喷管特征线逆推区型面终点至喷管特征线逆推区轴线终点之间出口边界的位置信息和参数分布信息；步骤7：基于喷管特征线逆推区壁面终点至喷管特征线逆推区轴线终点之间出口边界以及喷管特征线逆推区轴线终点至喷管理论出口点之间的出口边界进行特征线推进，获得喷管特征线逆推区壁面终点至喷管理论出口点之间的第二型面坐标信息；步骤8：分别利用第一喷管型面坐标信息和第二型面坐标信息对喷管初始型面中
对应的型面进行替换获得喷管最终型面。
6.其中，本方法能够进一步提升半柔壁喷管流场均匀性和跨超声速流场区域共用性，满足大型超声速风洞设备对于超声速半柔壁喷管的建设需求。
7.本发明提供的半柔壁喷管设计方法，采用喉道声速近似解，能够充分考虑喷管喉道处真实流动分布，取消了常规设计方法所引入的泉流假设，可为喷管消波区的特征线迭代以及波系结构控制提供真实流动条件；基于beizer（贝塞尔）函数设定了轴向加速区的马赫数分布，实现了马赫数沿轴向单调上升，避免了压缩波集中产生于喷管轴线，并通过特征线线逆推方法，削弱了喷管菱形区过渡区的范围和强度，实现喷管流场均匀性的提升。
8.优选的，所述步骤1具体为采用sivells方法生成各马赫数下的喷管初始型面。利用sivells方法可以快速生成初始型面，但是sivells方法存在相应的缺陷，因此利用后续步骤进行调整和优化。
9.优选的，所述步骤1包括：从喷管喉道至喷管最大转折点，使用多项式确定喷管型面，通过型面计算可确定喷管最大转折点和喷管理论出口点位置，通过马赫角关系计算得到特征线逆推区轴线终点的位置。
10.优选的，在所述特征线关系图中，o点为喷管喉道坐标点，i点为起始点，h点为初始膨胀区壁面起点，a点为喷管最大转折点，v点为初始膨胀区轴线终点，s点为转折推进区轴线终点，w点为特征线逆推区壁面终点，t点为特征线逆推区轴线终点，u点为喷管理论出口点；h点、i点、v点和a点作为区域hiva的4个顶点形成喷管初始膨胀区hiva，a点、v点和s点作为区域avs的3个顶点形成喷管转折推进区avs，a点、s点、w点和t点作为区域aswt的4个顶点形成喷管特征线逆推区aswt，w点、t点和u点作为区域wtu的4个顶点形成喷管反射波消波区wtu，其中，沿喷管轴线喷管初始膨胀区hiva、喷管转折推进区avs、喷管特征线逆推区aswt和喷管反射波消波区wtu依次分布。
11.优选的，喷管喉部声速线采用以下方程表示：其中，、、和均为常数；其中，，为坐标原点，为起始线原点；其中，为喉道半高；起始线采用以下方程表示：其中，采用以下方程表示：
其中，为喉部曲率半径。
12.优选的，所述步骤4包括：根据喉道声速近似解，计算出i点的位置，计算通过起始线oi的流量；确定一个等分n_hi，作为hi的离散点数，hi为i点的左伸特征线；预估h点的高度，计算得到等分高度dy，从i点开始，由左伸特征线方程计算下一点的位置和参数，结合喉道声速近似解推进，并计算每一推进步的流量积分，根据流量守恒，最终计算得到与固壁边界oa的交点h；据计算获得的h点位置，计算更新等分高度dy，计算入口边界hi上的离散点获得喉道初始边界；在区域hiva中，利用入口边界hi和喷管壁面条件oa，利用特征线方法进行推进计算获得av边界参数，av边界为区域hiva的出口边界。
13.优选的，所述步骤5包括：计算得到avs区域的所有离散点参数和坐标，as为a点的右行特征线，与喷管中心线相交于s点，计算得到在s点上的马赫数一阶和二阶导数。
14.优选的，所述步骤6基于5次贝塞尔曲线获得喷管转折推进区轴线终点至喷管特征线逆推区轴线终点之间的马赫数分布信息。在st马赫数分布具体计算过程中，考虑到对s点二阶导数的敏感性，s点只考虑一阶导数连续，t点一、二阶导数连续，则可采用5次beizer曲线拟合，有一个自由控制顶点作为优化参数使用。
15.优选的，所述步骤7包括：通过边界wt参数和出口边界ut参数，通过特征线推进，利用质量守恒关系计算出型面wu坐标，边界ut为直线，边界ut线上气流参数均相同，计算每一层的流量积分方式为利用总流量减去从t点沿tu方向到所在层位置的部分流量。
16.优选的，在完成每一层的推进计算后，均对本层的离散点进行重构。在每推进一层的计算过程中，离散点的个数增加一个，因此，为控制每层的离散点数量，在完成每一层的推进计算后，均需要对本层的离散点采用线性差值或二阶差值进行重构。
17.本发明提供的一个或多个技术方案，至少具有如下技术效果或优点：本方法在喉道处采用跨声速近似解来模拟声速条件下近壁面真实流动。首先，通过设置近似微分方程，将流场扰动势函数展开，形成以轴向距离x的函数为系数，以型面高度为因变量的幂级数，并通过边界条件得到喉道声速条件下的真实流动分布。然后，以该流动分布为入口条件，通过设定喷管轴向及出口马赫数分布，采用特征线逆推方法，实现消波区和全消波区的有效组合。
18.通过本方法能够消除半柔壁喷管喉道处泉流假设，使得喷管喉道处流动分布更加符合真实情况；同时，能够改进喷管轴向马赫数分布形式，降低喷管不同马赫数下菱形过渡区的范围和强度，从而进一步提升喷管流场均匀性，提升跨声速和超声速流场共用性。
附图说明
19.此处所说明的附图用来提供对本发明实施例的进一步理解，构成本发明的一部分，并不构成对本发明实施例的限定；图1是本方法的流程示意图；图2是sivells方法特征线示意图；图3是一般特征线示意图；图4是喷管喉道区几何形状示意图；图5是hiva区域特征点推进示意图；图6是avs区域特征点推进示意图；图7是astw、wtu区域特征点推进示意图；图8是本方法所得到的曲线簇示意图。
具体实施方式
20.为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明进行进一步的详细描述。需要说明的是，在相互不冲突的情况下，本发明的实施例及实施例中的特征可以相互组合。
21.在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述范围内的其他方式来实施，因此，本发明的保护范围并不受下面公开的具体实施例的限制。
22.实施例一请参考图1，图1为本方法的流程示意图，本方法具体包括如下步骤：a、使用sivells方法生成各马赫数下的喷管初始型面：从喷管喉道至转折点a，使用多项式确定喷管型面，通过型面计算可确定a、u点位置，通过马赫角关系计算得到点t的位置（见图2，图2为sivells方法特征线示意图）。
23.其中，图2中，a点为喷管最大转折点；b点为部分消波区控制点；c点为通过pt、qr特征线交点s的流线与特征线ar的交点，与b点相互确定，为转折消波关系点/控制点；r点为泉流区终点，为最大转折点a的特征线与中轴线的交点；o点为泉流区起始点，即泉流区与中轴线交点；p点为简单波区壁面起点，也可叫为特征点；s点为特征线pt和rq的交点，为反射膨胀波分隔点/波分割点；t点为p点发出的特征线与中轴线的交点，为试验区起始点；u点为喷管理论出口点；q点为简单波区壁面终点，为特征线rq与喷管壁面的交点；y点为过s点的流线与特征线ut的交点，为延伸点；喉道至oa：泉流区/径向流发展区；ap：部分消波区；pq：简单波区；qu：反射消波区；y*：为喉道高度；h：为喷管出口高度。
24.b、得到初始喷管型面后，定义喷管各关键控制点，确定一般特征线关系图，一般特征线关系图如图3所示。曲线oi为基于声速近似解方法确定的计算起始线，i点为起始线与喷管轴线的交点，hi为i点的左伸特征线，h点为i点左伸特征线与喷管壁面的交点；a点为喷管最大转折点，由sivells喷管型面计算方法得到，av为区域hiva的出口边界，v点的左伸特征线与喷管壁面相交于a点；oh和ha段，由sivells方法得到，为一段多项式曲线；以oi为入口条件，oa为固壁边界条件，使用特征线法，确定边界iv和av的参数值；s点为a点右伸特征线与喷管轴线的交点。图3中，o点为喷管喉道坐标点；i点为起始点；h点为初始膨胀区壁面
起点；a点为喷管最大转折点；v点为初始膨胀区轴线终点；s点为转折推进区轴线终点；w点为特征线逆推区壁面终点；t点为特征线逆推区轴线终点；u点为喷管理论出口点；hiva为喷管初始膨胀区；avs为喷管转折推进区；aswt为喷管轴向马赫数设定区/喷管膨胀过渡区/喷管特征线逆推区；wtu为喷管反射波消波区。
25.c、确定喷管喉部声速线oi，计算声速线上的参数，作为特征线往下推进的起始边界。
26.d、oa段作为喷管固壁边界，与喉部声速线一起，通过特征线往下推进，可求出区域hiva的下游边界av（av为v点的左伸特征线，与喷管壁面交于a点），确定av边界的坐标点位置以及参数值。
27.e、通过特征线推进，求出区域avs的下游边界as。as为a点的右行特征线，与喷管中心线相交于s点。通过计算，可得到中心线vs上的马赫数分布，可确定s点的马赫数以及一阶、二阶导数，以用于ｓｔ段马赫数分布设定的函数边界条件。
28.f、区域aswt，出口边界wt的参数计算。基于5次beizer曲线设定st中心线上的马赫数分布，由边界as和st，通过特征线推进，利用通过每条特征线的质量守恒关系式，计算出喷管壁面aw的位置坐标，以及出口边界wt的位置和参数分布。
29.g、区域wtu，计算wu段喷管壁面坐标。基于wt边界和ut出口边界，通过特征线推进，利用质量守恒关系计算出wu型面坐标。
30.h、以计算得到的aw、wu型面代替sivells方法得到的ap、pq、qu型面，oa保持原有型面，则得到了新的经过修正后的喷管位流型面。
31.实施例二在实施例一的基础上，本实施例结合具体的数据和公式，对本发明中的方法进行具体登记介绍和说明。
32.本实施例中的方法包括：a．使用sivells方法生成设定马赫数下的喷管初始型面：从喷管喉道至转折点a，使用多项式确定喷管型面，通过型面计算可确定a、u点位置，通过马赫角关系计算得到点t的位置（见图2）。
33.（1）初始曲线
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（1）边界条件：其中，a、b、c、d、e为公式系数；x、y为型面坐标；xa为喷管转折点轴向坐标；ya为喷管
最大转折点纵向坐标；为喷管最大转折点气流方向；为初始曲线的一阶导数；为初始曲线的二阶导数；为初始曲线在初始段终点时的三阶导数；为初始曲线在终止段处的三阶导数；为部分消波区控制点处的气流方向；通过前几个初始条件，可确定各系数。ya和的选定，要根据特征线法确定。
34.引入可做出喉道曲率与转折点处曲率斜率的关系图。对应每一个膨胀角，存在一个，使得常数a为零，从而三次多项式能满足所有马赫数和膨胀角的条件，而气动型面曲线可用弹性板曲线精确地吻合。
35.（2）计算半消波区ap之间曲线，引入： (2)；其中，为部分消波区初始边界处特征线上的任意一点d的气流方向；为d点的沿x的轴向坐标；为最大转折点a点的轴向坐标；为试验段半高；为t点处与喉道半高面积比；为d点处与喉道半高面积比；为a点处与喉道半高面积比；为g点处与喉道半高面积比；为半消波区内任意一点g点的马赫数；为g点处气流方向；为积分符号；为泉流区d点气流方向；为d点处壁面y向坐标。
36.其中，满足：；公式（2）中的未定义系数由公式（3）～（7）求解： (3)； (4)； (5)； (6)； (7)；为常数，设计输入，为某点处马赫数，为t点处普朗特迈耶函数，为a点处普朗特迈耶函数，为某点处与喉道面积比值，为比热比，常数，1.4，为普
朗特迈耶函数。
37.（3）计算完全消波区pq之间曲线：(8)；其中，为d点处x向坐标，为a点处与喉道面积比，为k点处气流方向，为w点处与喉道面积比，为w点处马赫数，为w点处气流方向，为d点处y向坐标；其中，mg见公式（4）；公式（8）中的未定义系数由公式（11）～（13）求解： (9)； (10)； (11)；（4）计算之间参数曲线：(12)；其中，为d点处x向坐标，为r点处与喉道面积比，为k点处气流方向，为w点处与喉道面积比，为w点处马赫数，为半消波区内任一点v点处与喉道面积比，为g点处气流方向，为d点处气流方向，为d点处y向坐标，为v点处气流方向，为v点处马赫数，为g点处马赫数。
38.其中，满足：，公式（12）未定义系数见公式（13）～（20）； (13)； (14)； (15)； (16)； (17)； (18)； (19)； (20)；（5）计算之间型面曲线：
ꢀꢀꢀꢀ
(21)；其中，系数满足边界条件：
ꢀꢀꢀꢀ
(22)；其中，为ut段y坐标，为ut段以x为自变量的y函数坐标，至分别为多项式系数，为ut段函数的导数，为y坐标的微分形式，为x坐标的微分形式，为与差值处的x坐标，为与差值处的y坐
标，为ut段以t点处的x坐标为自变量的y坐标函数，为ut段以t点处的x坐标为自变量的y坐标函数的导数，为ut段以u点处的x坐标为自变量的y坐标函数的二阶导数。
39.（6）分别对oa、ap、pq、qu、ut等曲线进行线积分，求和后得到ot长度s；（7）调整试验段马赫数m
t
，确保所有设计m
t
中所得到ot长度s与最高马赫数条件下一致。如果不一致，则调整当前m
t
数下值并重新计算，直至所有设计马赫数下ot长度一致。
40.b．得到初始喷管型面后，定义喷管各关键控制点。曲线oi为基于声速近似解方法确定的计算起始线，i点为起始线与喷管轴线的交点，hi为i点的左伸特征线，h点为i点左伸特征线与喷管壁面的交点；a点为喷管转折点，由sivells喷管型面计算方法得到，av为区域hiva的出口边界，v点的左伸特征线与喷管壁面相交于a点。oh和ha段，由sivells方法得到，为一段多项式曲线。以oi为入口条件，oa为固壁边界条件，使用特征线法，确定边界iv和av的参数值。特征线关系如图3所示。
41.c．确定喷管喉部声速线oi，计算声速线上的参数，作为特征线往下推进的起始边界。
42.图4是喷管喉部示意图，喉道高度为，喉道型面曲率半径为。声速线与喷管中心线交于o点，o点为分析流场所用的坐标系原点，流场的全速势方程为： (23)其中，为全速度势函数对x的偏微分，为全速度势函数对y的偏微分，为全速度势函数对x和y的偏微分，为常数，为全速度势函数对x的二阶偏微分，为常数，为全速度势函数对y的二阶偏微分，式中：，对应平面流动，对应轴对称流动。
43.根据喷管喉部的流动特点，即流动速度与音速相近，可将两个速度分量表示为： (24)其中，为速度在x方向的分量，为速度势函数的偏微分，为对x的偏微分，为常数，声速，为速度势函数对x的偏微分，为速度在y方向的分量，为对y的偏微分，为速度势函数对y的偏微分。
44.式中，，。此外，，对x，y的偏导数也是小量。引进扰动速度势，它与扰动速度分量的关系为： (25)将全速势的偏导数用扰动势的偏导数表示，利用小扰动线化方法，忽略二阶以上的小量，得到喷管喉部近声速区的扰动速度势方程：
ꢀ
(26)将上式用幂级数展开，可以得到扰动速度分量： (27) (28)其中，为扰动速度在x方向的分量的无量纲化，为扰动速度在x方向的分量的无量纲化，为常数，半径和喉道高度之和，为势函数的无量纲化，为势函数，为势函数微分的无量纲化，为坐标原点，为起始线原点，为曲率半径，为喉道半高。
45.式中，为无量纲参数，其表达式为：为无量纲参数，其表达式为：在声速线上，ma=1，可以确定声速线的位置：在声速线上，，于是，由式（27）得以抛物线表示的声速线方程： (29)令根据喉部边界条件，由（28）式得到：
ꢀꢀ
(30)在具体计算跨声速流场参数时，需要在喉部区选择一个ma》1的起始位置作为计算起始线。一种简单的办法就是选择所对应的曲线作为起始线。则根据式（28），起始线方程为：(31)比较式（30）和（31），可以看出起始线为一条通过t点和o点的抛物线，位于声速线下游。起始线上的速度值由式（27）和（28）完全确定。
46.设喉部的一小段型面可用圆弧表示，则给定型面的曲率半径，可求出待定系数。图7为推导与待定系数关系的示意图，根据几何关系，喉部型面的曲率可表示为：(32)其中，ds为喉道定点处的一段弧线，t为喉道顶点。
47.一般情况下，喉部型面的曲率半径比喉道高度大得多，则有：
根据式（28）可简化得到喉部型面曲率与常数α的关系式： (33)一般情况下，是未知的，而喉部曲率半径是已知的，则： (34)由此，喷管喉部区域的速度分布式（27）和（28），以及声速线方程和起始线方程均可由喉部参数和表示，这样就可得到起始线为界的特征线入口边界条件。
48.d．oa段作为喷管固壁边界，与喉部声速线一起，通过特征线往下推进，可求出区域hiva的下游边界av（见图4，图4为 喷管喉道区几何形状示意图，av为v点的左伸特征线，与喷管壁面交于a点），确定av边界的坐标点位置以及参数值。
49.(1)根据喉道声速近似解，计算出i点的位置，计算通过起始线oi的流量；(2)确定一个等分n_hi，作为hi的离散点数；预估一个h点的高度，得到等分高度dy，从i点开始，由左伸特征线方程（35）计算下一点的位置和参数，结合喉道声速近似解，逐步推进，并计算每一推进步的流量积分，根据流量守恒，最终计算得到与喷管壁面oa的交点h；在特征线上，满足特征线方程和相容性条件方程：左伸特征线： (35)右伸特征线： (36)为马赫数角；(3)根据确定的h点位置，重新计算等分高度dy，ih边界上的离散点重新计算一遍，以作为喉道初始边界。
50.在区域hiva中，利用入口边界hi和喷管壁面条件oa，可使用特征线方法进行层层推进计算（如图5所示），最终求得av边界参数。
51.下面以hi往下游推进一层为例：（1）边界iv上离散点参数的求解。边界hi上是等距离分布的离散点，从标号为1的点，向喷管轴线上发出右伸特征线与轴线交于点2，通过下列方程组的迭代计算可求出点2的参数： (37)式中：，，。两个方程，两个未知数、，一
般迭代10次以内，结果可满足计算精度要求，可确定点2的所有参数。
52.其中，为1点处x坐标，为2点处x坐标，为1点处y坐标，为2点处y坐标，为12点连线的气流方向，为12点马赫角差值，为2点处气流方向，为1点处气流方向，为2点处弦长；为1点处弦长；为1点弦长和2点弦长之和。
53.（2）内点4参数的求解。在点2参数求出以后，由于hi上的离散点参数已知，通过迭代，可由点3和点2求出未知内点点4参数：右伸特征线： (38) 左伸特征线： (39)式中下标24为点2、4对应参数的平均值。
54.（3）边界点6参数的求解。点5参数由上述方法求出以后，通过左伸特征线关系式： (40)可得两个独立方程。oa型面的每一局部区域，可由多次曲线表示，由sivells得到的oa段离散点，可作为局部拟合的控制点，通过两步法（第一步：第一次试算，通过点5的位置即点5-6连线角度，初步确定点6在oa上的x坐标的大概值；其中，为3点处x坐标，为4点处x坐标，为3点处y坐标，为4点处y坐标，为3、4点方向角之和，为3、4点马赫数角之和，为2、3点方向角之和，为2、4点马赫角之和，为2、3点马赫数角之和，为4点方向角，为3点方向角，为5点y坐标，为6点y坐标，为5、6点方向角之和，为5、6点马赫数角之和，为6点方向角，为5点方向角，为6点弦长，为5点弦长，为56点弦长之和。
55.第二步：第二次校正计算，通过x坐标确定点6 附近的几个离散点，拟合多次曲线，作为进一步迭代求解的方程），求得oa段局部拟合方程。由于点6在oa上的切向角即为气流的方向角，因此，四个方程，四个未知量，通过迭代可求出点6的参数。
56.e、通过特征线推进，求出区域avs的下游边界as。as为a点的右行特征线，与喷管中心线相交于s点。通过计算，可得到中心线vs上的马赫数分布，可确定s点的马赫数以及一阶、二阶导数，以用于st段马赫数分布设定的函数边界条件。
57.avs区域为一三角形区域，从超声速区影响域来说，该区域内参数由av参数分布确定。av边界的离散数确定了，则边vs和as的离散点数也确定了。
58.下面简述一下计算流程，以av边界下一推进层的计算为例（见图6）：(1)边界点2，由av边界上的点1右伸特征线与喷管轴线相交得到；内点3，由上一层的点4右伸特征线与点2的左伸特征线相交得到，该层的其它内点以此类推；该层的最后一点点6，由边界点a的右伸特征线与内点5的左伸特征线相交得到。需要注意，该层的离散点，
比上一层少一个点。
59.(2)as上的边界点9，由7的左伸特征线和点8的右伸特征线相交得到；端点s由点9右伸特征线与轴线相交得到。
60.通过上述计算流程，可得到avs区域的所有离散点参数和坐标。边界av上的离散点数和离散点分布，决定了边界as和vs上的离散点数和离散点的分布，该三角区域边界上的离散点数相同。在vs边界上，s点及其附近的点位置及参数均已知，可得到在s点上的马赫数一阶、二阶导数，用于计算st上的马赫数分布。
61.f． 区域astw，出口边界wt的参数计算。设定st中心线上的马赫数分布，由边界as和st，通过特征线推进，利用通过每条特征线的质量守恒关系式，计算出喷管壁面aw的位置坐标，以及出口边界wt的位置和参数分布。
62.astw区域的入口边界as上的参数已知，结合喷管中心线上的st边界，通过特征线法和流量守恒关系，可计算出喷管型面坐标。如图7所示。
63.中心轴线上的st边界，马赫数分布需要额外确定。为了保证轴线上马赫数分布单调上升，并保证出口马赫数过渡到设计值，一般要求在st之间二阶导数须小于0。在端点s和t处，要求保持马赫数的一阶、二阶连续，确保在喷管轴线上没有压缩波的集中，并保证马赫数分布的光滑连续性。拟合st上的马赫数分布采用n次beizer函数，构造喷管轴线上的马赫数分布。beizer函数能保证其二阶导数在整个轴线上连续，这将使得喷管内部参数不会发生突变，并且由于beizer曲线本身的特性，不仅可以通过调节控制顶点，保证曲线的单调性，而且可以灵活地控制曲线的形状。beizer曲线采用特征多边形顶点的位置矢量与伯恩斯坦基函数的线性组合表达。
64.定义：给定n 1个空间向量，称n次参数曲线段 (41)为一条n次的beizer曲线，称为控制顶点。特别当伯恩斯坦系数为实数时，可以取为控制顶点。
[0065] (42)
ꢀꢀ
(43)beizer函数在端点有如下性质：在t=0时，，,；在t=1时，，，
。
[0066]
其中，为t点空间向量，为t点伯恩斯坦基函数，为阶乘系数，为比例参数，为控制顶点起点；为0次伯恩斯坦系数；为控制顶点导数；为1次控制顶点伯恩斯坦系数；为端点控制顶点二阶导数；为2次伯恩斯坦系数； 为控制顶点终点，为n次伯恩斯坦系数，为控制顶点终点导数，为n-1次伯恩斯坦系数，为控制顶点终点二阶导数，为n-2次伯恩斯坦系数。
[0067]
由于在端点须保证马赫数导数的相同，可通过上式求得相邻控制顶点参数。关于beizer函数的次数，有报告建议采用6次beizer函数，有7个控制顶点，若考虑端点一、二阶导数连续，则占去6个控制点，剩下中间一个控制点可自由调节，以优化参数。在st马赫数分布具体计算过程中，考虑到对s点二阶导数的敏感性，s点只考虑一阶导数连续，t点一、二阶导数连续，则可采用5次beizer曲线拟合，有一个自由控制顶点作为优化参数使用。
[0068]
下面就astw区域的参数计算进行描述（以as的下一层线为例）：（1）内点2，由s点的左伸特征线，与1、2点连线相交得到，点1为点2右伸特征线与轴线交点，点s和点1坐标参数为已知。根据两点具体坐标，可通过拟合公式得到马赫数值。通过下述方程，迭代计算可得到点2参数：右伸特征线： (44) 左伸特征线： (45)其中，为s点x坐标，为s点y坐标，为2点和s点方向角之和，为2点和s点弦长之和，为s点弦长，为2点和s点马赫角之和，为s点方向角。
[0069]
内点4可经由点2、点3通过上述同样方法计算得到。
[0070]
（2）边界点8，为喷管型面上的边界点。由a点的左伸特征线，与6、7点连线相交，点6为内点，由上一步计算得，由点a和点6得到点7。由于通过入口边界as的流量与推进的每一层的流量积分相同，通过流量插值关系，按流量比例迭代求出点6、7之间的壁面点8。
[0071]
需要注意的是，在每推进一层的计算过程中，离散点的个数增加一个，因此，为控制每层的离散点数量，在完成每一层的推进计算后，均需要对本层的离散点采用线性差值或二阶差值进行重构。
[0072]
g．区域wtu，计算wu段喷管壁面坐标。通过wt边界和ut出口边界，通过特征线推进，利用质量守恒关系计算出wu型面坐标。如图7所示。
[0073]
边界ut为一直线，线上气流参数均相同，t点和u点已经确定，计算每一层的流量积分时，应该是总流量减去从t点沿tu方向到所在层位置这部分流量。
[0074]
（1）内点11，由入口边界9点的左伸特征线，与11、10点连线相交得到，点10为点11右伸特征线与边界ut交点，点9和点10坐标参数为已知。通过下述方程，迭代计算可得到点
11参数：右伸特征线： (46) 左伸特征线： (47)其中，为10点处x坐标，为10点处y坐标，为11点处x坐标，为11点处y坐标，为10点和11点马赫角之和，为10点和11点方向角之和，为10点处方向角， 为11点处方向角，为9点处x坐标，为9点处y坐标，为11点处x坐标，为11点处y坐标，为9点和11点方向角之和，为9点和11点马赫数角之和， 为9点和11点弦长之和，为10点和11点弦长之和，为9点方向角。
[0075]
（2）边界点15，为喷管型面上的边界点。由w点的左伸特征线，与14、13点连线相交，点13为内点，由上一步计算得，由点w和点13得到点14。由于通过入口边界wt的流量与推进的每一层的流量积分相同，通过流量插值关系，按流量比例迭代求出点13、14之间的壁面点15相关参数。
[0076]
需要注意的是，在每推进一层的计算过程中，离散点的个数减少一个，因此，为控制每层的离散点数量和间距，在完成每一层的推进计算后，均需要对本层的离散点进行重构。重构的方法很多，采用线性插值或二阶插值等方法都行。
[0077]
h．以计算得到的aw、wu型面代替sivells方法得到的ap、pq、qu型面，oa保持原有型面，则得到了新的经过修正后的喷管位流型面。图8给出了改进设计方法所得到的曲线簇。
[0078]
尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。
[0079]
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。