CN102999697A - 一种高超声速喷管多阶连续气动型面的确定方法 - Google Patents

Abstract

本发明提供一种高超声速喷管多阶连续气动型面的确定方法，给定喷管出口参数，通过计算出边界条件，建立特征线网格，计算通过特征线的流量，利用质量守恒原理得到喷管型面曲线，再计算出各点处附面层位移厚度，无粘型面加上附面层位移厚度，即得到喷管的物理型面，根据本发明方法设计的喷管，不仅可以保证气动型面的多阶连续，且没有经验公式的影响，喷管出口气流均匀分布。

Description

一种高超声速喷管多阶连续气动型面的确定方法

技术领域

本发明涉及一种高超声速喷管多阶连续气动型面的确定方法，属于流体动力技术领域。

背景技术

喷管型面设计方法最先由普朗特(Prandtl)和贝斯曼(Busemann)提出，他们依据特征线方法完成了图解法，为喷管型面设计技术打下基础。稍后，派克特(Puckett)将图解法修改为半图解法。该方法设计喷管需要大量精确的绘图工作，同时还含有未定量的误差，尤其当图线接近对称轴的时候更加明显，设计过程枯燥且复杂。为解决图解法设计喷管的弊端，喷管设计方法经过了一系列的改进，如解析法，经验公式法、半理论半经验公式法，分析设计法等。现常采用的喷管设计技术方案主要有两种，Foelsch设计法和Cresci设计法。

Foelsch继承前人的思想，提出了喷管的解析设计法，该方法是一种近似简化的方法，喷管外形坐标可以由两个简单方程迅速求出。其主要思想是通过假定在转折点区域的流动是源流而简化型面的计算。源流下游的边界是转折点发出到对称轴上的特征线。沿着这条特征线流动参数很容易计算出来；并且因为在源流区域的下游的所有左行特征线都是直的，下游的全部型线可以通过解析的方法决定出。在转折点的上游，通过采用一个简单的经验曲线(圆弧加直线或三次曲线)，使得源流可以产生。其示意图如图1所示。

Cresci设计法也采用源流假设，通过在源流的下游区域引进一个过渡区获得连续的曲率型线。在转折点的上游，采用一条三次曲线使得源流产生，同时转折点下游型线舍弃Foelsch法采用的近似的解析解，而是通过已知边界条件，建立特征线网格，根据质量守恒原理确定型面坐标。其示意图如图2所示。

现有的技术中，Foelsch设计方法型面曲率的不连续会影响流场品质且不适合挠性喷管的设计。Cresci设计方法虽然可以获得连续曲率的型面，但转折点前的曲线是经验曲线。选择的经验曲线不一定能保证源流的形成或者源流区域很小，从而影响特征线参数的求取及消波区中消除膨胀波效果，最终的流场品质有所下降。

发明内容

本发明的目的：克服现有技术的不足，提供了一种高超声速喷管多阶连续气动型面的确定方法，克服现有技术在喷管设计中型面曲率不连续和需要选择经验曲线来完成喷管设计的问题，完全通过计算获取多阶连续曲率的喷管型面。

本发明的技术解决方案：

一种高超声速喷管多阶连续气动型面的确定方法，包括如下步骤：

(1)建立坐标系，以喷管入口和出口中心点的连线为x轴，喷管出口的方向为x轴的正方向，以过喉道顶点T垂直于x轴的直线为y轴，喉道顶点T位于y轴的正半轴上，x轴和y轴的交点为原点O，

G点和A点为所述高超声速喷管多阶连续气动型面上的点，A点位于G点和D点之间，D为喷管出口点，T、G、A、D四点的横坐标值依次增大，其中D点为喷管型面上横坐标值最大的点，E点、B点和C点均为x轴上的点，且E点、B点和C点的坐标依次增大；

(2)确定边界右行特征线TI的参数，I为x轴上的点，具体为：

(2.1)通过公式 ${\left(\frac{y_D}{y^{*}}\right)}^{1+\mathrm{\σ}}=\frac{1}{M_D}{(\frac{\mathrm{\γ}-1}{\mathrm{\γ}+1}{M_D}^2+\frac{2}{\mathrm{\γ}+1})}^{\frac{\mathrm{\γ}+1}{2(\mathrm{\γ}-1)}}$ 计算出喉道半高y^*，其中，y_D是预设的出口半径或半高，M_D为预设的喷管出口马赫数，γ为比热比，当喷管为平面喷管时，σ＝0；当喷管为轴对称喷管时，σ＝1；

(2.2)把喉道半高y^*分成N份，依次得到每点的纵坐标值y_j，j＝1，2，3，4…N，N为正整数，y₁即为喉道顶点T的纵坐标；

(2.3)通过跨声速流动解计算出T点的气流角θ₁、马赫角μ₁和马赫数M₁：

(2.4)通过公式 $x_2=x_1+\frac{2(y_2-y_1)}{\tan ({\mathrm{\θ}}_1-{\mathrm{\μ}}_1)+\tan ({\mathrm{\θ}}_2-{\mathrm{\μ}}_2)}$ 计算第2个点的横坐标x₂，其中，x₁为T点的横坐标，θ₂、μ₂分别为第2个点的假设气流角和假设马赫角；

(2.5)根据步骤(2.4)中确定的第2点的坐标，通过跨声速理论计算出第2点的马赫数M₂和气流角θ₂，继而根据公式 $v_2=\sqrt{\frac{\mathrm{\γ}+1}{\mathrm{\γ}-1}}\arctan \sqrt{\frac{\mathrm{\γ}-1}{\mathrm{\γ}+1}({M_2}^2-1)}-\arctan \sqrt{{M_2}^2-1}$ 计算得到第2点的普朗特-梅耶角v₂；

(2.6)通过跨声速流动解结合特征线方程的方式再次求解普朗特-梅耶角v₂，即通过公式

$v_2=v_1+{\mathrm{\θ}}_1-{\mathrm{\θ}}_2+\frac{\mathrm{\σ}}{2}[\frac{v_1}{w_1{y}_1{M}_1}+\frac{v_2}{w_2{y}_2{M}_2}]{[{(y_2-y_1)}^2+{(x_2-x_1)}^2]}^{\frac{1}{2}}$ 再次计算普朗特-梅耶角v₂，其中，v₁为T点的普朗特-梅耶角，且

$v_1=\sqrt{\frac{\mathrm{\γ}+1}{\mathrm{\γ}-1}}\arctan \sqrt{\frac{\mathrm{\γ}-1}{\mathrm{\γ}+1}({M_1}^2-1)}-\arctan \sqrt{{M_1}^2-1},$ θ₂为步骤(2.5)计算得到的气流角θ₂，w₁为T点的速度比，w₂为第2点的速度比，可分别通过M₁和M₂计算得到；

(2.7)比较步骤(2.5)中得到的普朗特-梅耶角和步骤(2.6)中得到的普朗特-梅耶角是否相等，若相等，则假设的气流角θ₂和马赫角μ₂正确，且假设的气流角θ₂与步骤(2.5)中计算得到的气流角θ₂相等，则确定得到了第2点的坐标、马赫数和气流角；

若不相等，则重新假设气流角θ₂和马赫角μ₂的值并且返回步骤(2.4)；

(2.8)以步骤(2.4)-(2.7)相同的方式依次确定第3，4，5…N点坐标、马赫数和气流角，即确定了由喉道顶点T发出的右行特征线TI上的参数，其中，第N点即为x轴上的I点。

(3)通过公式 $\left\{\begin{array}{cc}{r}_1=\frac{y^{*}}{{\mathrm{\θ}}_G}& \mathrm{\σ}=0\\ r_1=\frac{y^{*}}{2\sin \frac{{\mathrm{\θ}}_G}{2}}& \mathrm{\σ}=1\end{array}\right.$ 计算源流区马赫数为1的圆的半径r₁；其中，θ_G为预设气流角，即最大膨胀角；

(4)确定边界IE的参数，通过如下步骤进行：

(4.1)令边界IE的速度分布满足w＝C₀+C₁ζ+C₂ζ²+C₃ζ³+C₄ξ⁴，其中，w为速度比，ζ＝(x-x_I)/(x_E-x_I)，x_I为I点的横坐标，x_E为E点的横坐标，x为x_I和x_E之间的任一点的横坐标，四次多项式的系数C₀～C₄依次为：C₀＝W_I，C₁＝X_EW_I′，

C₃＝4(W_E-C₀)-W_E′X_E-3C₁-2C₂，C₄＝-3(W_E-C₀)+W_E′X_E+2C₁+C₂；

其中I点的速度比W_I、速度比的一阶导数W_I′和二阶导数W″使用跨声速理论得到；通过公式 $({W_I}^{\′\′}-{W_E}^{\′\′})X_E^2+6(W_I^{\′}+W_E^{\′})X_E-12(W_E-W_I)=0$ 计算得到X_E，W_E′和W_E″分别为E点速度比的一阶导数和二阶导数，由面积比和速度比之间的关系式求出；

(4.2)通过公式X_E＝(x_E-x_I)/r₁计算E点的横坐标x_E，

(4.3)根据x的坐标，通过公式w＝C₀+C₁ζ+C₂ζ²+C₃ζ³+C₄ζ⁴可以计算得到边界IE上的速度比，坐标，马赫数，气流角，且其中气流角为0；

(4.4)坐标原点O和源流源点O₁之间的距离x₀由下式求出：

x₀＝r_E-X_E×r₁-x_I

r_E在选取E点马赫数之后通过公式 ${\left(\frac{r}{r_1}\right)}^{1+\mathrm{\σ}}=\frac{1}{M}{(\frac{\mathrm{\γ}-1}{\mathrm{\γ}+1}{M}^2+\frac{2}{\mathrm{\γ}+1})}^{\frac{\mathrm{\γ}+1}{2(\mathrm{\γ}-1)}}$ 确定；

(5)确定特征线EG的参数；具体为步骤如下：

(5.1)令G点的预设气流角θ_G为最大膨胀角，将最大膨胀角分成K份，则分成K份之后，特征线EG上每点的气流角为已知，K为正整数；

(5.2)根据公式v-v_E＝(σ+1)θ计算得到气流角θ所对应的点的普朗特-梅耶角v，其中，v_E为E点的普朗特-梅耶角，可通过预设的E点的马赫数计算得到；

(5.3)根据步骤(5.2)中得到的所述气流角θ所对应的点的普朗特-梅耶角v，计算得到该点的马赫数M；

(5.4)通过公式 ${\left(\frac{r}{r_1}\right)}^{1+\mathrm{\σ}}=\frac{1}{M}{(\frac{\mathrm{\γ}-1}{\mathrm{\γ}+1}{M}^2+\frac{2}{\mathrm{\γ}+1})}^{\frac{\mathrm{\γ}+1}{2(\mathrm{\γ}-1)}}$ 计算所述气流角θ所对应的点到源流源点O₁的距离r；其中，源流源点O₁为A、G两点连线的延长线与x轴的交点，r₁为源流区马赫数为1的圆的半径；

(5.5)通过公式x＝r×cosθ-x₀和y＝r×sinθ将r，θ转换成该点的直角坐标系上的坐标(x，y)，x₀为坐标原点和源流源点的距离，则特征线EG上的参数坐标、马赫数和气流角均已确定；

(6)确定特征线AB的参数；具体为步骤如下：

(6.1)令A点的预设气流角θ_A为最大膨胀角，且θ_A＝θ_G将最大膨胀角分成K份，则分成K份之后，特征线AB上每点的气流角为已知，K为正整数；

(6.2)根据公式v_B-v＝(σ+1)θ计算得到气流角θ所对应的点的普朗特-梅耶角v，其中，v_B为B点的普朗特-梅耶角，可通过预设的B点的马赫数计算得到；

(6.3)根据步骤(6.2)中得到的所述气流角θ所对应的点的普朗特-梅耶角v，计算得到该点的马赫数M；

(6.4)通过公式 ${\left(\frac{r}{r_1}\right)}^{1+\mathrm{\σ}}=\frac{1}{M}{(\frac{\mathrm{\γ}-1}{\mathrm{\γ}+1}{M}^2+\frac{2}{\mathrm{\γ}+1})}^{\frac{\mathrm{\γ}+1}{2(\mathrm{\γ}-1)}}$ 计算所述气流角θ所对应的点到源流源点O₁的距离r；其中，源流源点O₁为A、G两点连线的延长线与x轴的交点，r₁为源流区马赫数为1的圆的半径；

(6.5)通过公式x＝r×cosθ-x₀和y＝r×sinθ将r，θ转换成该点的直角坐标系上的坐标(x，y)，x₀为坐标原点和源流源点的距离，则特征线AB上的参数坐标、马赫数和气流角均已确定；

(7)确定边界BC的参数，通过如下步骤进行：

(7.1)令边界BC的速度分布满足M＝D₀+D₁ξ+D₂ξ²+D₃ξ³+D₄ξ⁴，其中，M为马赫数，ξ＝(x-x_B)/(x_C-x_B)，x_B为B点的横坐标，x_C为C点的横坐标，x为x_B和x_C之间的任一点的横坐标，四次多项式的系数D₀～D₄依次为：D₀＝M_B，D₁＝X_CM_B′，

D₃＝4(M_C-D₀)-M_C′X_C-3D₁-2D₂，D₄＝-3(M_C-D₀)+M_C′X_C+2D₁+D₂；

其中B点的马赫数M_B、马赫数的一阶导数M_B′和二阶导数M″由面积比和马赫数之间的关系式求出；M_C等于设计马赫数M_D；

通过公式 $({W_B}^{\′\′}-{W_C}^{\′\′})X_C^2+6({W_B}^{\′}+{W_C}^{\′})X_C-12(W_C-W_B)=0$ 计算得到X_C，M_C′和M_C″分别为C点马赫数的一阶导数和二阶导数，均为0；

(7.2)通过公式X_C＝(x_C-x_B)/r₁计算C点的横坐标x_C，

(7.3)根据x的坐标，通过公式M＝D₀+D₁ξ+D₂ξ²+D₃ξ³+D₄ξ⁴可以计算得到边界BC上的坐标，马赫数，气流角，且其中气流角为0；

(8)确定边界CD的参数，具体为：CD为直特征线，特征线上的马赫数为设计马赫数M_D，气流角为0，CD的倾斜角为出口马赫数的马赫角，即得到CD的斜率，则得到CD直线方程，即得到了CD上的点的坐标；

(9)在确定边界条件TI、IE、EG、AB、BC、CD之后，建立特征线网格，特征线网格TIEG由E点出发，向上游推进建立，特征线网格ABCD由B点出发，向下游推进建立，然后由质量守恒原理确定出无粘型面坐标；

(10)通过公式

计算喷管出口的边界层位移厚度

令喉道边界层位移厚度为0，喉道到出口间的边界层位移厚度增长为线性关系，即得到各点处的边界层位移厚度，其中p_D为D点运动粘性，u_D为D点速度x轴分量，x_D为D点横坐标，无粘型面坐标加上位移厚度，即得到所述高超声速喷管多阶连续气动型面。

本发明的有益效果：

本发明方法通过假设先存在一个流场性能优越的喷管，该喷管存在一个合理的轴向马赫数或者速度分布，通过反推的办法确定出喷管的气动型线，不仅可以保证气动型面的多阶连续，且没有经验公式的影响，喷管出口气流均匀分布。

附图说明

图1为Foelsch法示意图。

图2为Cresci法示意图。

图3为本发明的喷管设计法示意图。

图4为在不同最大膨胀角及曲率半径比下，E点马赫数可选择的范围，其中图4a为在最大膨胀率为8deg时，曲率半径比和E点马赫数的关系；图4b为在最大膨胀率为12deg时，曲率半径比和E点马赫数的关系。

图5常见E点马赫数选取。

图6为本发明流程图。

具体实施方式

对于一个喷管，必定有对应的轴向马赫数或者速度分布，若给出一个合理的轴向马赫数或者速度分布，则可以确定出一个喷管的气动型线。本方案就是通过这个原理，通过设置轴向马赫数或者速度分布而计算出具有多阶连续曲率的型面。其示意图如图3所示。

用跨声速流动解结合特征线方程来确定特征线TI上的参数，假定区域GEBA为源流流动区，利用源流理论结合特征线理论计算EG，AB上的参数，并按照设计的马赫数计算特征线CD上的流动参数，沿IE的4次速度多项式分布和沿BC的4次马赫数多项式分布按沿喷管中心线的距离函数给出，并保证速度分布或者马赫数分布多项式在端点I、E、B、C处函数、一阶导数、二阶导数连续。E点速度、速度导数由源流解求出，B点马赫数、马赫数导数由源流解求出，I点处的速度、速度导数由跨声速流动解求得，C点处的马赫数为设计马赫数，马赫数的一阶导数、二阶导数为0，那么可以确定IE、BC上的参数。这样根据这些设定的边界条件，建立特征线网格TIEG和ABCD，利用质量守恒原理求解这两个区域的流场，从而得到喷管型面曲线TG和AD，GA在源流区，为一条直线，连接G点和A点即得到GA。结合TG、GA、AD即获得喷管的无粘型面TD。再计算出各点处附面层位移厚度，无粘型面加上附面层位移厚度，即得到喷管的物理型面。

一种高超声速喷管多阶连续气动型面的确定方法，如图6所示，包括如下步骤：

(1)建立坐标系，以喷管入口和出口中心点的连线为x轴，喷管出口的方向为x轴的正方向，以过喉道顶点T垂直于x轴的直线为y轴，喉道顶点T位于y轴的正半轴上，x轴和y轴的交点为原点0，

G点和A点为所述高超声速喷管多阶连续气动型面上的点，A点位于G点和D点之间，D为喷管出口点，T、G、A、D四点的横坐标值依次增大，其中D点为喷管型面上横坐标值最大的点，E点、B点和C点均为x轴上的点，且E点、B点和C点的坐标依次增大；

(2)确定边界右行特征线TI的参数，I为x轴上的点。特征线TI上的参数由跨声速流动解结合特征线方程来确定。在跨声速流动区域，利用跨声速流动解，区域中任一点的速度分量仅是喉道曲率半径与喉道半高之比R和该点坐标的函数，因此在选择好喉道曲率半径与喉道半高之比后，若已知喉道附近任意一点的坐标位置，则可得到该点的速度分量。

具体为：

(2.1)通过公式 ${\left(\frac{y_D}{y^{*}}\right)}^{1+\mathrm{\σ}}=\frac{1}{M_D}{(\frac{\mathrm{\γ}-1}{\mathrm{\γ}+1}{M_D}^2+\frac{2}{\mathrm{\γ}+1})}^{\frac{\mathrm{\γ}+1}{2(\mathrm{\γ}-1)}}$ 计算出喉道半高y^*，其中，y_D是预设的出口半径或半高，M_D为预设的喷管出口马赫数，γ为比热比，当喷管为平面喷管时，σ＝0；当喷管为轴对称喷管时，σ＝1；

(2.2)把喉道半高y^*分成N份，依次得到每点的纵坐标值y_j，j＝1，2，3，4…N，N为正整数，y₁即为喉道顶点T的纵坐标；

(2.3)通过跨声速流动解计算出T点的气流角θ₁、马赫角μ₁和马赫数M₁；

(2.4)通过公式 $x_2=x_1+\frac{2(y_2-y_1)}{\tan ({\mathrm{\θ}}_1-{\mathrm{\μ}}_1)+\tan ({\mathrm{\θ}}_2-{\mathrm{\μ}}_2)}$ 计算第2个点的横坐标x₂，其中，x₁为T点的横坐标，θ₂、μ₂分别为第2个点的假设气流角和假设马赫角；

(2.5)根据步骤(2.4)中确定的第2点的坐标，通过跨声速理论计算出第2点的马赫数M₂和气流角θ₂，继而根据公式 $v_2=\sqrt{\frac{\mathrm{\γ}+1}{\mathrm{\γ}-1}}\arctan \sqrt{\frac{\mathrm{\γ}-1}{\mathrm{\γ}+1}({M_2}^2-1)}-\arctan \sqrt{{M_2}^2-1}$ 计算得到第2点的普朗特-梅耶角v₂；

(2.6)通过跨声速流动解结合特征线方程的方式再次求解普朗特-梅耶角v₂，即通过公式

$v_2=v_1+{\mathrm{\θ}}_1-{\mathrm{\θ}}_2+\frac{\mathrm{\σ}}{2}[\frac{v_1}{w_1{y}_1{M}_1}+\frac{v_2}{w_2{y}_2{M}_2}]{[{(y_2-y_1)}^2+{(x_2-x_1)}^2]}^{\frac{1}{2}}$ 再次计算普朗特-梅耶角v₂，其中，v₁为T点的普朗特-梅耶角，且

$v_1=\sqrt{\frac{\mathrm{\γ}+1}{\mathrm{\γ}-1}}\arctan \sqrt{\frac{\mathrm{\γ}-1}{\mathrm{\γ}+1}({M_1}^2-1)}-\arctan \sqrt{{M_1}^2-1},$ θ₂为步骤(2.5)计算得到的气流角θ₂，w₁为T点的速度比，w₂为第2点的速度比，可分别通过M₁和M₂计算得到；

(2.7)比较步骤(2.5)中得到的普朗特-梅耶角和步骤(2.6)中得到的普朗特-梅耶角是否相等，若相等，则假设的气流角θ₂和马赫角μ₂正确，且假设的气流角θ₂与步骤(2.5)中计算得到的气流角θ₂相等，则确定得到了第2点的坐标、马赫数和气流角；

若不相等，则重新假设气流角θ₂和马赫角μ₂的值并且返回步骤(2.4)；

(2.8)以步骤(2.4)-(2.7)相同的方式依次确定第3，4，5…N点坐标、马赫数和气流角，即确定了由喉道顶点T发出的右行特征线TI上的参数，其中，第N点即为x轴上的I点。

(3)通过公式 $\left\{\begin{array}{cc}{r}_1=\frac{y^{*}}{{\mathrm{\θ}}_G}& \mathrm{\σ}=0\\ r_1=\frac{y^{*}}{2\sin \frac{{\mathrm{\θ}}_G}{2}}& \mathrm{\σ}=1\end{array}\right.$ 计算源流区马赫数为1的圆的半径r₁；其中，θ_G为预设气流角，即最大膨胀角；

(4)确定边界IE的参数。

通过如下步骤进行：

(4.1)令边界IE的速度分布满足w＝C₀+C₁ζ+C₂ζ²+C₃ζ³+C₄ζ⁴，其中，w为速度比，ζ＝(x-x_I)/(x_E-x_I)，x_I为I点的横坐标，x_E为E点的横坐标，x为x_I和x_E之间的任-点的横坐标，四次多项式的系数C₀～C₄依次为：C₀＝W_I，C₁＝X_EW_I′，

C₃＝4(W_E-C₀)-W_E′X_E-3C₁-2C₂，C₄＝-3(W_E-C₀)+W_E′X_E+2C₁+C₂；

其中I点的速度比W_I、速度比的一阶导数W_I′和二阶导数W″使用跨声速理论得到；通过公式 $({W_I}^{\′\′}-{W_E}^{\′\′})X_E^2+6(W_I^{\′}+W_E^{\′})X_E-12(W_E-W_I)=0$ 计算得到X_E，W″_E和W_E″分别为E点速度比的一阶导数和二阶导数，由面积比和速度比之间的关系式求出；

(4.2)通过公式X_E＝(x_E-x_I)/r₁计算E点的横坐标x_E，

(4.3)根据x的坐标，通过公式w＝C₀+C₁ζ+C₂ζ₂+C₃ζ₃+C₄ζ₄可以计算得到边界IE上的速度比，坐标，马赫数，气流角，且其中气流角为0；

(4.4)坐标原点O和源流源点O₁之间的距离x₀由下式求出：

x₀＝r_E-X_E×r₁-x_I

r_E在选取E点马赫数之后通过公式 ${\left(\frac{r}{r_1}\right)}^{1+\mathrm{\σ}}=\frac{1}{M}{(\frac{\mathrm{\γ}-1}{\mathrm{\γ}+1}{M}^2+\frac{2}{\mathrm{\γ}+1})}^{\frac{\mathrm{\γ}+1}{2(\mathrm{\γ}-1)}}$ 确定；

(5)确定特征线EG的参数；具体为步骤如下：

(5.1)令G点的预设气流角θ_G为最大膨胀角，将最大膨胀角分成K份，则分成K份之后，特征线EG上每点的气流角为已知，K为正整数；

(5.2)根据公式v-v_E＝(σ+1)θ计算得到气流角θ所对应的点的普朗特-梅耶角v，其中，v_E为E点的普朗特-梅耶角，可通过预设的E点的马赫数计算得到；

(5.3)根据步骤(5.2)中得到的所述气流角θ所对应的点的普朗特-梅耶角v，计算得到该点的马赫数M；

(5.4)通过公式 ${\left(\frac{r}{r_1}\right)}^{1+\mathrm{\σ}}=\frac{1}{M}{(\frac{\mathrm{\γ}-1}{\mathrm{\γ}+1}{M}^2+\frac{2}{\mathrm{\γ}+1})}^{\frac{\mathrm{\γ}+1}{2(\mathrm{\γ}-1)}}$ 计算所述气流角θ所对应的点到源流源点O₁的距离r；其中，源流源点O₁为A、G两点连线的延长线与x轴的交点，r₁为源流区马赫数为1的圆的半径，即 $\left\{\begin{array}{cc}{r}_1=\frac{y^{*}}{{\mathrm{\θ}}_G}& \mathrm{\σ}=0\\ r_1=\frac{y^{*}}{2\sin \frac{{\mathrm{\θ}}_G}{2}}& \mathrm{\σ}=1\end{array}\right.,$ y^*为喉道半径；

(5.5)通过公式x＝r×cosθ-x₀和y＝r×sinθ将r，θ转换成该点的直角坐标系上的坐标(x，y)，x₀为坐标原点和源流源点的距离，则特征线EG上的参数坐标、马赫数和气流角均已确定；

(6)确定特征线AB的参数；具体为步骤如下：

(6.1)令A点的预设气流角θ_A为最大膨胀角，且θ_A＝θ_G将最大膨胀角分成K份，则分成K份之后，特征线AB上每点的气流角为已知，K为正整数；

(6.2)根据公式v_B-v＝(σ+1)θ计算得到气流角θ所对应的点的普朗特-梅耶角v，其中，v_B为B点的普朗特-梅耶角，可通过预设的B点的马赫数计算得到；

(6.3)根据步骤(6.2)中得到的所述气流角θ所对应的点的普朗特-梅耶角v，计算得到该点的马赫数M；

(6.4)通过公式 ${\left(\frac{r}{r_1}\right)}^{1+\mathrm{\σ}}=\frac{1}{M}{(\frac{\mathrm{\γ}-1}{\mathrm{\γ}+1}{M}^2+\frac{2}{\mathrm{\γ}+1})}^{\frac{\mathrm{\γ}+1}{2(\mathrm{\γ}-1)}}$ 计算所述气流角θ所对应的点到源流源点O1的距离r；其中，源流源点O₁为A、G两点连线的延长线与x轴的交点，r₁为源流区马赫数为1的圆的半径，即 $\left\{\begin{array}{cc}{r}_1=\frac{y^{*}}{{\mathrm{\θ}}_A}& \mathrm{\σ}=0\\ r_1=\frac{y^{*}}{2\sin \frac{{\mathrm{\θ}}_A}{2}}& \mathrm{\σ}=1\end{array}\right.,$ y^*为喉道半径；

(6.5)通过公式x＝r×cosθ-x₀和y＝r×sinθ将r，θ转换成该点的直角坐标系上的坐标(x，y)，x₀为坐标原点和源流源点的距离，则特征线AB上的参数坐标、马赫数和气流角均已确定；

(7)确定边界BC的参数，通过如下步骤进行：

(7.1)令边界BC的速度分布满足M＝D₀+D₁ξ+D₂ξ²+D₃ξ³+D₄ξ⁴，其中，M为马赫数，ξ＝(x-x_B)/(x_C-x_B)，x_B为B点的横坐标，x_C为C点的横坐标，x为x_B和x_C之间的任一点的横坐标，四次多项式的系数D₀～D₄依次为：D₀＝M_B，D₁＝X_CM_B′，D₃＝4(M_C-D₀)-M_C′X_C-3D₁-2D₂，D₄＝-3(M_C-D₀)+M_C′X_C+2D₁+D₂；

其中B点的马赫数M_B、马赫数的一阶导数M_B′和二阶导数M″由面积比和马赫数之间的关系式求出；M_C等于设计马赫数M_D；

通过公式 $({W_B}^{\′\′}-{W_C}^{\′\′})X_C^2+6({W_B}^{\′}+{W_C}^{\′})X_C-12(W_C-W_B)=0$ 计算得到X_C，M_C′和M_C″分别为C点马赫数的一阶导数和二阶导数，均为0；

(7.2)通过公式X_C＝(x_C-x_B)/r₁计算C点的横坐标x_C，

(7.3)根据x的坐标，通过公式M＝D₀+D₁ξ+D₂ξ²+D₃ξ³+D₄ξ⁴可以计算得到边界BC上的坐标，马赫数，气流角，且其中气流角为0；

(8)确定边界CD的参数，具体为：CD为直特征线，特征线上的马赫数为设计马赫数M_D，气流角为0，CD的倾斜角为出口马赫数的马赫角，即得到CD的斜率，则得到CD直线方程，即得到了CD上的点的坐标；

(9)在确定边界条件TI、IE、EG、AB、BC、CD之后，建立特征线网格，特征线网格TIEG由E点出发，向上游推进建立，特征线网格ABCD由B点出发，向下游推进建立，然后由质量守恒原理确定出无粘型面坐标；

(10)通过公式

计算喷管出口的边界层位移厚度

令喉道边界层位移厚度为0，喉道到出口间的边界层位移厚度增长为线性关系，即得到各点处的边界层位移厚度，其中p_D为D点运动粘性，u_D为D点速度x轴分量，x_D为D点横坐标，无粘型面坐标加上位移厚度，即得到所述高超声速喷管多阶连续气动型面。

B、E两点马赫数的选择要在合理的范围之内。E点马赫数选择和喉道曲率半径与喉道半高之比R、最大膨胀角η(θ_A)相关，E点常见的马赫数选择如图4，5所示。

对于B点马赫数的选取，在轴对称喷管中，通常取M_B＝(0.75～0.85)M_C，平面喷管中，M_B＝(0.85～0.95)M_C。选取完毕，还需检查马赫数的二阶导数在B和C之间必须小于零。

Claims (8)

1.一种高超声速喷管多阶连续气动型面的确定方法，其特征在于，包括如下步骤：

(1)建立坐标系，以喷管入口和出口中心点的连线为x轴，喷管出口的方向为x轴的正方向，以过喉道顶点T垂直于x轴的直线为y轴，喉道顶点T位于y轴的正半轴上，x轴和y轴的交点为原点O，

G点和A点为所述高超声速喷管多阶连续气动型面上的点，A点位于G点和D点之间，D为喷管出口点，T、G、A、D四点的横坐标值依次增大，其中D点为喷管型面上横坐标值最大的点，E点、B点和C点均为x轴上的点，且E点、B点和C点的坐标依次增大；

(2)确定边界右行特征线TI的参数，I为x轴上的点；

(3)计算源流区马赫数为1的圆的半径r₁；

(4)确定边界IE的参数；

(5)确定特征线EG的参数；

(6)确定特征线AB的参数；

(7)确定边界BC的参数；

(8)确定边界CD的参数；

(9)在确定边界条件TI、IE、EG、AB、BC、CD之后，建立特征线网格，特征线网格TIEG由E点出发，向上游推进建立，特征线网格ABCD由B点出发，向下游推进建立，然后由质量守恒原理确定出无粘型面坐标；

(10)通过公式 

计算喷管出口的边界层位移厚度 

令喉道边界层位移厚度为0，喉道到出口间的边界层位移厚度增长为线性关系，即得到各点处的边界层位移厚度，其中p_D为D点运动粘性，u_D为D点速度x轴分量，x_D为D点横坐标，无粘型面坐标加上位移厚度，即得到所述高超声速喷管多阶连续气动型面。

2.根据权利要求1所述的一种高超声速喷管多阶连续气动型面的确定方法，其特征在于：所述步骤(2)确定边界右行特征线TI的参数，具体为步骤如下：

(2.1)通过公式

计算出喉道半高y^*，其中，y_D 是预设的出口半径或半高，M_D为预设的喷管出口马赫数，γ为比热比，当喷管为平面喷管时，σ＝0；当喷管为轴对称喷管时，σ＝1；

(2.2)把喉道半高y^*分成N份，依次得到每点的纵坐标值y_j，j＝1，2，3，4…N，N为正整数，y₁即为喉道顶点T的纵坐标；

(2.3)通过跨声速流动解计算出T点的气流角θ₁、马赫角μ₁和马赫数M₁；

(2.4)通过公式

计算第2个点的横坐标x₂，其中，x₁为T点的横坐标，θ₂、μ₂分别为第2个点的假设气流角和假设马赫角；

(2.5)根据步骤(2.4)中确定的第2点的坐标，通过跨声速理论计算出第2点的马赫数M₂和气流角θ₂，继而根据公式 

计算得到第2点的普朗特-梅耶角v₂；

(2.6)通过跨声速流动解结合特征线方程的方式再次求解普朗特-梅耶角v₂，即通过公式

再次计算普朗特-梅耶角v₂，其中，v₁为T点的普朗特-梅耶角，且

θ₂为步骤(2.5)计算得到的气流角θ₂，w₁为T点的速度比，w₂为第2点的速度比，可分别通过M₁和M₂计算得到；

(2.7)比较步骤(2.5)中得到的普朗特-梅耶角和步骤(2.6)中得到的普朗特-梅耶角是否相等，若相等，则假设的气流角θ₂和马赫角μ₂正确，且假设的气流角θ₂与步骤(2.5)中计算得到的气流角θ₂相等，则确定得到了第2点的坐标、马赫数和气流角；

若不相等，则重新假设气流角θ₂和马赫角μ₂的值并且返回步骤(2.4)；

(2.8)以步骤(2.4)-(2.7)相同的方式依次确定第3，4，5…N点坐标、马赫数和气流角，即确定了由喉道顶点T发出的右行特征线TI上的参数，其中，第N点即为x轴上的I点。

3.根据权利要求1所述的一种高超声速喷管多阶连续气动型面的确定方法，其特征在于：所述步骤(3)中计算源流区马赫数为1的圆的半径r₁具体为：通过公式 计算得到；其中，θ_G为预设气流角，即最大膨胀角，y^*为喉道半径。

4.根据权利要求1所述的一种高超声速喷管多阶连续气动型面的确定方法，其特征在于：所述步骤(4)中确定边界工E的参数，通过如下步骤进行：

(4.1)令边界IE的速度分布满足w＝C₀+C₁ζ+C₂ζ²+C₃ζ³+C₄ζ⁴，其中，w为速度比，ζ＝(x-x_I)/(x_E-x_I)，x_I为I点的横坐标，x_E为E点的横坐标，x为x_I和x_E之间的任一点的横坐标，四次多项式的系数C₀～C₄依次为：C₀＝W_I，C₁＝X_EW_I′， 

C₃＝4(W_E-C₀)-W_E′X_E-3C₁-2C₂，C₄＝-3(W_E-C₀)+W_E′X_E+2C₁+C₂；

其中I点的速度比W_I、速度比的一阶导数W_I′和二阶导数W″使用跨声速理论得到；通过公式

计算得到X_E，W_E′和W_E″分别为E点速度比的一阶导数和二阶导数，由面积比和速度比之间的关系式求出；

(4.2)通过公式X_E＝(x_E-x_I)/r₁计算E点的横坐标x_E，

(4.3)根据x的坐标，通过公式w＝C₀+C₁ζ+C₂ζ²+C₃ζ³+C₄ζ⁴可以计算得到边界IE上的速度比，坐标，马赫数，气流角，且其中气流角为0；

(4.4)坐标原点O和源流源点O₁之间的距离x₀由下式求出：x₀＝r_E-X_E×r₁-x_Ir_E在选取E点马赫数之后通过公式

确定。

5.根据权利要求1所述的一种高超声速喷管多阶连续气动型面的确定方法，其特征在于：所述步骤(5)中确定特征线EG的参数，具体为步骤如下：

(5.1)令G点的预设气流角θ_G为最大膨胀角，将最大膨胀角分成K份，则分成K份之后，特征线EG上每点的气流角为已知，K为正整数；

(5.2)根据公式v-v_E＝(σ+1)θ计算得到气流角θ所对应的点的普朗特-梅耶角v，其中，v_E为E点的普朗特-梅耶角，可通过预设的E点的马赫数计算得到； 

(5.3)根据步骤(5.2)中得到的所述气流角θ所对应的点的普朗特-梅耶角v，计算得到该点的马赫数M；

(5.4)通过公式

计算所述气流角θ所对应的点到源流源点O₁的距离r；其中，源流源点O₁为A、G两点连线的延长线与x轴的交点，r₁为源流区马赫数为1的圆的半径。

(5.5)通过公式x＝r×cosθ-x₀和y＝r×sinθ将r，θ转换成该点的直角坐标系上的坐标(x，y)，x₀为坐标原点和源流源点的距离，则特征线EG上的参数坐标、马赫数和气流角均已确定。

6.根据权利要求1所述的一种高超声速喷管多阶连续气动型面的确定方法，其特征在于：所述步骤(6)中确定特征线AB的参数，具体为步骤如下：

(6.1)令A点的预设气流角θ_A为最大膨胀角，且θ_A＝θ_G将最大膨胀角分成K份，则分成K份之后，特征线AB上每点的气流角为已知，K为正整数；

(6.2)根据公式v_B-v＝(σ+1)θ计算得到气流角θ所对应的点的普朗特-梅耶角v，其中，v_B为B点的普朗特-梅耶角，可通过预设的B点的马赫数计算得到；

(6.3)根据步骤(6.2)中得到的所述气流角θ所对应的点的普朗特-梅耶角v，计算得到该点的马赫数M；

(6.4)通过公式

计算所述气流角θ所对应的点到源流源点O₁的距离r；其中，源流源点O₁为A、G两点连线的延长线与x轴的交点，r₁为源流区马赫数为1的圆的半径；

(6.5)通过公式x＝r×cosθ-x₀和y＝r×sinθ将r，θ转换成该点的直角坐标系上的坐标(x，y)，x₀为坐标原点和源流源点的距离，则特征线AB上的参数坐标、马赫数和气流角均已确定。

7.根据权利要求1所述的一种高超声速喷管多阶连续气动型面的确定方法，其特征在于：所述步骤(7)中确定边界BC的参数，通过如下步骤进行：

(7.1)令边界BC的速度分布满足M＝D₀+D₁ξ+D₂ξ²+D₃ξ³+D₄ξ⁴，其中，M为马赫数，ξ＝(x-x_B)/(x_C-x_B)，x_B为B点的横坐标，x_C为C点的横坐标，x为x_B和x_C之间 的任一点的横坐标，四次多项式的系数D₀～D₄依次为：D₀＝M_B，D₁＝X_CM_B′，

D₃＝4(M_C-D₀)-M_C′X_C-3D₁-2D₂，D₄＝-3(M_C-D₀)+M_C′X_C+2D₁+D₂；

其中B点的马赫数M_B、马赫数的一阶导数M_B′和二阶导数M″由面积比和马赫数之间的关系式求出；M_C等于设计马赫数M_D；

通过公式

计算得到X_C，M_C′和M_C″分别为C点马赫数的一阶导数和二阶导数，均为0；

(7.2)通过公式X_C＝(x_C-x_B)/r₁计算C点的横坐标x_C，

(7.3)根据x的坐标，通过公式M＝D₀+D₁ξ+D₂ξ²+D₃ξ³+D₄ξ⁴可以计算得到边界BC上的坐标，马赫数，气流角，且其中气流角为0。

8.根据权利要求1所述的一种高超声速喷管多阶连续气动型面的确定方法，其特征在于：所述步骤(8)确定边界CD的参数，具体为：CD为直特征线，特征线上的马赫数为设计马赫数M_D，气流角为0，CD的倾斜角为出口马赫数的马赫角，即得到CD的斜率，则得到CD直线方程，即得到了CD上的点的坐标。 

Images