二维超声速普朗特

右g卷女亭课程设计题目：二维超声速普朗特-迈耶系数波流场的数值解学院：飞行器工程学院专业名称：飞行器设计与工程班级学号：学生姓名：指导教师：

第一部分1.物理问题简介：普朗特一一迈耶稀疏波的解析解图-1中，超声速流围绕着一个尖的扩张角膨胀，无数个无限弱的马赫波组成了稀疏波，在尖角处展开成扇形。扇形稀疏波的波头与来流方向的夹角..],而..2是其波尾与下游方向的夹角。“1和..2称为马赫角，定义为：=罚■:sin亍和和二=—占■fillMa1和Ma2分别为上下游的马赫数。通过稀疏波的流动是等熵流动。当流体通过稀疏波后，马赫数增加，压力、温度和密度降低；图-1中标明了这些变化趋势。在中心稀疏波前的流动是均匀的，马赫数为Ma1，而且流动平行与波前的壁面。稀疏波后的流动是均匀的，马赫数为Ma2,并且平行于下游的壁面。在稀疏波内，流动参数光滑变化，流线弯曲，如图-1所示。稀疏波内的流动是二维的，唯一的例外是折角的定点，它是一个奇点，壁面流线的方向在此处有一个突然的变化，而且此处的流动参数也是不连续的。这个奇点对流动的数值解会产生影响。给定超声速来流条件和拐角处的偏转角三，下游参数是唯一确定的。对于完全气体，在稀疏波后的流动有精确的解析解，下面给出这个解。流过中心稀疏波的流动，其解析解取决于简单的关系式式中，f是普朗特一一迈耶函数；E是流动偏转角。对于完全气体，普朗特一一迈耶函数是Ma和Y的函数，定义为f=arctan(Ma2)-arrtanVMa2-12解析解中如下依次得到。对给定的Ma1，从式(2)计算函数匕。然后，对给定的偏转角。，从式(1)得到f2。用这样得到f2的值，通过求解式(2)求出Ma2。式(2)是关于Ma2的隐式关系式，需要用试凑法求解。波后的压力温度和密度都可以由

等熵流动关系式：P2-P1［不和状态方程：得到。借助是式(1)〜式(5)，中心稀疏波后的流动就完全确定了。2.问题的提法考虑图-2所示的物理平面。来流马赫数为2,来流的压力、密度、温度分别为：1.01x105N/m2、1.23kg/m3、286.1K。超声速流动的扩张角。=5.325。，计算区域为：x=0到x=60m，壁面到y=40m，如图-4所示。计算区域内的压力、密度、温度、马赫数等。扩张角定点位置是x=10m温度、马赫数等。扩张角定点位置是x=10m。，h=h(x)为40m40+(x-101tan9初值线。初值线在65nn初值线。初值线在65nnx=0处，在位于这条铅七=->垂线的网格点上，初值，"-二\「-pi=1.23kg/nn3由来流给定。计算从这-条线开始并以X为步以便更明确这个问题的求解。长向下游推进。以便更明确这个问题的求解。为更好的解决这个问题，下面就对这个问题的解决办法提出一些理论上的内容，做好准备，第二部分普朗特一一迈耶稀疏波流场的数值解

1.控制方程定常二维流强守恒形式的控制方程组可以表示如下的通用形式:F和G为列向量，其中:考虑没有体积力的等熵流动。(6)式中的源项J等于零。把(7)式中的列向量每一个分量记作：=puF?=pu'十p%=puv(V2匕=pnIe-F—对于完全气体_RT1pefT-因此可消去式(7d)的e，最后得到pu+pupu+pu——7e同样的，可以得到:沿流向推进方法的基本思路。在方程(6)中，将x写在了导数的左边，y的导数写在了右边。看下图一3,如果沿着位于x。出的处之线给定流场变量是y的函数，那么沿着这条线可以求出方程(6)中的G的y方向导数，进而得到F的x方向导数。再由这些方向导数，就可以得到位于xo+Ax处的下一条铅垂线上的流场变量。按这种方式，可以从沿着初值线给定的流场开始，通过沿x方向以Ax为步长的推进得到全流场的解，如图-2。关于强守恒的试验表明，数值求解这种形式的方程，会出现一些额外的问题，即：①需要将通量F、F、F、F分解，才能求的原变量；②向量G的元素G、12341G、G、G只能用F、F、F、F来表示，而不是像式(8a〜8d)那样，用原始变2341234量来表示。下面讨论下这两个问题。对于第一个问题，从通量变量中求出原始变量，结果如下：-B+VBa-4AC2A其中:101112以及状态方程13如果用强守恒形式的控制方程进行数值求解，也就是求解方程(6)，直接得到的是通量F、F、F、F的值，不是原始变量的值。P、u、v、p和T的值必1234须由式⑼〜式(13)得出。对于第二个问题，即如何计算方程(6)中的G。在给定的网格点上，F「F2、%、F4的值可以直接从方程(6)的数值解得到，所以对于下一个网格点处的计算，将这些值用在G、G、G、G的计算中是有道理的。也就是说，应该用F、F、123412

TOC\o"1-5"\h\zF、F的值直接计算G、G、G、G。而不是先从式⑼〜式(13)中求的原变量，341234然后再用这些变量根据式(8a〜d)计算G、G、G、G。因为，G显然是F的函1234数。下面给出这种函数关系。首先，由式(8a)和式(11)，得到二=二「=瑚14由方程(7c)和式(8b)，可以得到G2，即：G2=F315由式(8c)和式(11)，有部=：■■■-p=矿予,—洪16利用式(7b)和式(10)，可以从式(16)中消去p，因为:§="-二二-="-曾17将式(17)代入式(16)，得==，号；_二……18最后，G4的表达式如下：乌"』隹丫川乌丫］192.网格生成与坐标变换为了给流过扩张角的流动建立有限差分解法，必须使用贴体网格系统，如图-4所示，用xy笛卡尔坐标系统表示的物理平面如图-4a所示。带有扩张角的壁面构成了物理平面的下边界。入流边界位于x=0，出流边界位于x=L。上边界为水平线y=H。很明显，由于扩张角下游的壁面乌"』隹丫川乌丫］192.网格生成与坐标变换为了给流过扩张角的流动建立有限差分解法，必须使用贴体网格系统，如图-4所示，用xy笛卡尔坐标系统表示的物理平面如图-4a所示。带有扩张角的壁面构成了物理平面的下边界。入流边界位于x=0，出流边界位于x=L。上边界为水平线y=H。很明显，由于扩张角下游的壁面□心计算平而图一4中心稀疏波数值解的计算网格分析图-4a，用h表示物理平面内从下表面到上表面边界的距离，有h=h（x）。用ys表示壁面的纵坐标，这里ys=ys（x）。由此，可以定义如下变换：g=x……20通过这个变换，计算平面内E在0到L之间变化，n在0到1之间变化；n二0对应物理平面内的壁面，n=i对应上边界。n和E为常数的线组成了计算平面内的有规则的正交网格。在E—n平面的正交网格上实施有限差分计算，控制流动的偏微分方程组是在转换后的平面内进行数值求解的。因此，为了能在计算平面内使用他们，必须做适当的变换。即，方程（6）必须变换成E和n表示的形式。导数的变换：葛=冒（蒙）+亲偿）上两式中的度量可以通过变换关系式（20）和式（21）得到，即:222423——X吕曳hdxhd瞥2524式（24）中的三可以可以写成更简单的形式，由图-4a，用x=E表示扩张角的位置，则

对x<E：y=0h二常数对x王E：y=-（x-E）tan。h=H+（x-E）tan。对这些表达式进行微分，得到：TOC\o"1-5"\h\z对x三E：M二。二二：□xdx对x王E：—=-三—=:dx.日is因此，度量h可以表示为：t?xP（当ME）….26a\o"CurrentDocument"独一j（l—ri）兰四（当kEE）……26hVh将式（22）、式（23）、式（24）、式（25）和式（26）代入式（a）、式（b）,就得到了导数的变换：\o"CurrentDocument"7-=--：-7.：-……2728式（27）中，二由式（26a）或式（26b）给出。dx.再来看看方程（6）给出的物理平面内的守恒形式的流动控制方程。由于J=0,方程写为：3FdG2930—=3dy2930利用式（27）、式（28）对方程（29）进行变换，得到：--—或—=-；—英\阪_/加hdv\小显[\djtJdv\h其中度量三由式（26a）或式项）给出。用向量F和G的分量来写，方程

（30）就是下面一组方程:……31……32……33……31……32……33……34X动量方程：告=-建专-：吝Y动量方程:号=_「三籍_：吝

能量方程：言=「三三1斗方程（31）到方程（34）就是要在图-3b所示的计算平面内求解的流动控制方程。3.推进步长的计算定常无粘超声速流的控制方程是双曲线方程，所以沿流向推进求解才是合适的。对于时间推进，根据CFL条件可以得到可允许的最大步长。指出：从物理概念上讲，显式时间推进可允许的的最大时间步长应该小于或等于，声波从一个网格点运动到相邻的网格点所需的时间。这种声波传播解释，使我们能够直观的确定定长流动的CFL条件。如图-5中显示了x站位上垂直排列的网格点。点1处的小扰动沿着该点处的两条特征线向外传播。特征线就是流动的马赫线，它和流线的夹角就是马赫角四。假设点1处流线与x轴的夹角为。，那么左行和右行马赫波与x轴的夹角分别是。+四和。-四。图-5只给出了点1处的左行马赫波。设有一条通过点2的水平线，点1的左行特征线与水平线相交与点a。于是，点a和点2之间的水平距离（Ax'为35TOC\o"1-5"\h\z根据点2处的CFL条件，为了稳定/A性，推进步长Ax的值不应该大雨Ax1；“因而点2和点a之间的距离应小于或至?多等于，声波从点1传到与点2同样高'的位置所需的距离。对于点3处的右行a/'"■顷囹方仙宜椎繇年的示意囹【物理甲面）马赫波，也有类似的结果，设它与过点2的水平线交于点b。点b和点2之间的水平距离二二"为二二：：二36、ta.n(S-p>B为保证点2处沿流向推进计算稳定性，步长-二的值不能大于二二"和危二"两者中较小的一个。将这种分析应用于x。处垂直排列的所有网格点，就能给出X0处沿流向推进的步长Ax，表达式为：二二=37|0舶±31球1上式中［浏二三二基是x=x0处垂直排列的所有网格点上:"技二的绝对中最大的。由于式(20)和(21)定义的坐标变换给出的E=x，那么图-4b所示的计算区域中，沿流向推进的步长为A(=Ax38式中Ax由式(37)确定。联立式(37)和式(38)，并引入柯朗数C，得到A£应满足的稳定条件39CFL条件要求上式中的C31。第三部分最终结果用空间推进计算得到的结果如下：图-6列出的是x=8.9375m，11.375m，23.5625m，31.6875m，47.7375m和64.1875m处的x方向速度分量u对纵坐标y的函数图。

7207LQ70069Qd=.-x=3].6B75m5——x=47.M75m7207LQ70069Qd=.-x=3].6B75m5——x=47.M75m6-mgz即mx=II.J75m——■>x=23.5tf2.j-m速度u/(m/S)730图-7列出的是x=8.9375m，11.375m，23.5625m，31.6875m，47.7375m和图-7压力p对纵坐标y的函数图-8列出的是x=8.9375m，11.375m，23.5625m，31.6875m，47.7375m和64.1875m处的温度T对纵坐标y的函数图。

290T/k28026D25DM3■3-一^23.5625m—--s=31.6S75ehi.…»=47.9375id1&75芋290T/k28026D25DM3■3-一^23.5625m—--s=31.6S75ehi.…»=47.9375id1&75芋131?Z1£5Mns37日】y/(y-L)in230图-9列出的是x=8.9375m，11.375m，23.5625m，31.6875m，47.7375m和64.1875m处的马赫数Ma对纵坐标y的函数图。1…….x=89375mL3-I1图-9马赫数Ma1…….x=89375mL3-I1第四部分小结通过这次课程设计，对空间推进的原理有了较深的理解。空间推进用的是定常流的守恒方程组。根据求解区域的形状，使用贴体坐标系。练习了网格生成的某些方法，并使用变换后的控制方程组。还用到了扑捉波的技