1、南昌航空大学飞行器工程学院课程设计课 程 设 计题 目：二维超声速普朗特-迈耶系数波流场的数值解学 院： 飞行器工程学院 专业名称： 飞行器设计与工程 班级学号： 07034211 学生姓名： 李 桂 平 指导教师： 刘 勇 二 O 一 O 年 十一 月南昌航空大学飞行器工程学院课程设计第一部分1物理问题简介：普朗特迈耶稀疏波的解析解图-1 中，超声速流围绕着一个尖的扩张角膨胀，无数个无限弱的马赫波组成了稀疏波，在尖角处展开成扇形。扇形稀疏波的波头与來流方向的夹角 1，而 2 是其波尾与下游方向的夹角。 1和 2称为 马赫角，定义为： 和1=11 2=12Ma1和 Ma2分别为上下游的马赫数。
2、通过稀疏波的流动是等熵流动。当流体通过稀疏波后，马赫数增加，压力、温度和密度降低；图-1 中标明了这些变化趋势。在中心稀疏波前的流动是均匀的，马赫数为 Ma1，而且流动平行与波前的壁面。稀疏波后的流动是均匀的，马赫数为 Ma2，并且平行于下游的壁面。在稀疏波内，流动参数光滑变化，流线弯曲，如图-1 所示。稀疏波内的流动是二维的，唯一的例外是折角的定点，它是一个奇点，壁面流线的方向在此处有一个突然的变化，而且此处的流动参数也是不连续的。这个奇点对流动的数值解会产生影响。给定超声速来流条件和拐角处的偏转角 ，下游参数是唯一确定的。对于完全气体，在稀疏波后的流动有精确的解析解，下面给出这个解。流过中
3、心稀疏波的流动，其解析解取决于简单的关系式12=1+式中，f 是普朗特迈耶函数； 是流动偏转角。对于完全气体，普朗特迈耶函数是 Ma 和 的函数，定义为2=+11+11（ 2） 21解析解中如下依次得到。对给定的 Ma1，从式(2)计算函数 f1。然后，对给定的偏转角 ，从式(1)得到 f2。用这样得到 f2 的值，通过求解式(2)求出 Ma2。式(2) 是关于 Ma2 的隐式关系式，需要用试凑法求解。波后的压力温度和密度都可以由等熵流动关系式：南昌航空大学飞行器工程学院课程设计 32=1 1+(1)/2211+(1)/222 1 41=21+(1)/2211+(1)/222和状态方程： 52
4、=22得到。借助是式(1)式(5)，中心稀疏波后的流动就完全确定了。2.问题的提法考虑图-2 所示的物理平面。来流马赫数为 2，来流的压力、密度、温度分别为：1.01x10 5N/m2、1.23kg/m 3、286.1K。超声速流动的扩张角 =5.325。 ，计算区域为：x=0 到 x=60m，壁面到 y=40m，如图-4 所示。计算区域内的压力、密度、温度、马赫数等。扩张角定点位置是 x=10m。此时，h=h(x)为=40 （ 010） 40+（ 10） tan （ 1065） 初值线。初值线在x=0 处，在位于这条铅垂线的网格点上，初值由来流给定。计算从这条线开始并以 x 为步长向下游推进
5、。为更好的解决这个问题，下面就对这个问 图-2题的解决办法提出一些理论上的内容，做好准备，以便更明确这个问题的求解。第二部分 普朗特迈耶稀疏波流场的数值解南昌航空大学飞行器工程学院课程设计1控制方程定常二维流强守恒形式的控制方程组可以表示如下的通用形式：6=F 和 G 为列向量，其中：7 =2+(+22)+8=2+(+22)+考虑没有体积力的等熵流动。(6)式中的源项 J 等于零。把(7)式中的列向量每一个分量记作：1= 72=2+ 73= 74=(+22)+ 7对于完全气体=1= 11因此可消去式(7d)的 e，最后得到7e4=1+2+22同样的，可以得到：1= 82= 83=2+ 84=(
6、+22)+=1+2+22 8沿流向推进方法的基本思路。南昌航空大学飞行器工程学院课程设计在方程(6)中，将 x 写在了导数的左边，y 的导数写在了右边。看下图3，如果沿着位于 x0出的处之线给定流场变量是 y 的函数，那么沿着这条线可以求出方程(6)中的 G 的 y 方向导数，进而得到 F 的 x 方向导数。再由这些方向导数，就可以得到位于 x0+x 处的下一条铅垂线上的流场变量。按这种方式，可以从沿着初值线给定的流场开始，通过沿 x 方向以x 为步长的推进得到全流场的解，如图-2。关于强守恒的试验表明，数值求解这种形式的方程，会出现一些额外的问题，即： 需要将通量 F1、F 2、F 3、F
7、4分解，才能求的原变量； 向量 G 的元 1 2素 G1、 G2、 G3、 G4只能用 F1、F 2、F 3、F 4来表示，而不是像式（8a8d）那样，用原始变量来表示。下面讨论下这两个问题。对于第一个问题，从通量变量中求出原始变量，结果如下：9=+242其中：; ; =23214 =112 =+12(1)3110=111=3112=21 以及状态方程 13=如果用强守恒形式的控制方程进行数值求解，也就是求解方程(6)，直接得到的是通量 F1、F 2、F 3、F 4的值，不是原始变量的值。、u、v、p 和 T 的值必南昌航空大学飞行器工程学院课程设计须由式(9)式(13)得出。对于第二个问题，
8、即如何计算方程(6)中的 G。在给定的网格点上，F1、F 2、F 3、F 4的值可以直接从方程(6)的数值解得到，所以对于下一个网格点处的计算，将这些值用在 G1、G 2、G 3、G 4的计算中是有道理的。也就是说，应该用 F1、F 2、F 3、F 4的值直接计算 G1、G 2、G 3、G 4。而不是先从式(9)式(13)中求的原变量，然后再用这些变量根据式（8ad）计算 G1、G 2、G 3、G 4。因为，G 显然是 F 的函数。下面给出这种函数关系。首先，由式(8a)和式(11)，得到141=31由方程(7c)和式(8b)，可以得到 G2，即：G2=F3 15由式(8c)和式(11)，有1
9、63=2+=(31)2+利用式(7b)和式(10)，可以从式(16)中消去 p，因为：17=22=221将式(17)代入式(16)，得183=(31)2+221最后，G4 的表达式如下：194=1+2+22 =1(221)31+231(1)2+(31)22.网格生成与坐标变换南昌航空大学飞行器工程学院课程设计为了给流过扩张角的流动建立有限差分解法，必须使用贴体网格系统，如图-4 所示，用 xy 笛卡尔坐标系统表示的物理平面如图-4a 所示。带有扩张角的壁面构成了物理平面的下边界。入流边界位于 x=0，出流边界位于 x=L。上边界为水平线 y=H。很明显，由于扩张角下游的壁面向下倾斜，物理平面不
10、能构成矩形网格。因此，必须把物理平面转化成计算平面，后者的有限差分网格是正交的，如图-4b 所示。计算平面自变量用 、 表示。分析图-4a，用 h 表示物理平面内从下表面到上表面边界的距离，有 h=h(x)。用 ys 表示壁面的纵坐标，这里 ys=ys(x)。由此，可以定义如下变换：=x 2021=()()通过这个变换，计算平面内 在 0 到 L 之间变化， 在 0 到 1 之间变化；=0 对应物理平面内的壁面，=1 对应上边界。 和 为常数的线组成了计算平面内的有规则的正交网格。在 平面的正交网格上实施有限差分计算，控制流动的偏微分方程组是在转换后的平面内进行数值求解的。因此，为了能在计算平
11、面内使用他们，必须做适当的变换。即，方程（6）必须变换成 和 表示的形式。导数的变换：a=()+()南昌航空大学飞行器工程学院课程设计b=()+()上两式中的度量可以通过变换关系式（20）和式（21）得到，即：22=123=024=1 25=1式（24）中的 可以可以写成更简单的形式，由图-4a，用 x=E 表示扩张角的位置，则对 x E： y s=0 h=常数对 x E： y s=-（x-E）tan h=H+（x-E）tan对这些表达式进行微分，得到：对 x E： =0 =0对 x E： =tan =tan因此，度量 可以表示为：=0 （ 当 ） 26（ 1） tan （当 ） 26将式（2
12、2） 、式（23） 、式（24） 、式（25）和式（26）代入式（a） 、式（b） ，就得到了导数的变换：27=+()南昌航空大学飞行器工程学院课程设计28=1式（27）中， 由式（26a）或式（26b）给出。再来看看方程（6）给出的物理平面内的守恒形式的流动控制方程。由于J=0，方程写为：29=利用式（27） 、式（28）对方程（29）进行变换，得到：或 30+()=1 =()+1其中度量 由式（ 26a）或式（26b）给出。用向量 F 和 G 的分量来写，方程（30）就是下面一组方程：连续性方程： 311=()1+11X 动量方程: 322=()2+12Y 动量方程: 333=()3+13
13、能量方程： 344=()4+14方程（31）到方程（34）就是要在图-3b 所示的计算平面内求解的流动控制方程。3.推进步长的计算定常无粘超声速流的控制方程是双曲线方程，所以沿流向推进求解才是合适的。对于时间推进，根据 CFL 条件可以得到可允许的最大步长。指出：从物理概念上讲，显式时间推进可允许的的最大时间步长应该小于或等于，声波从一个网格点运动到相邻的网格点所需的时间。这种声波传播解释，使我们能够南昌航空大学飞行器工程学院课程设计直观的确定定长流动的 CFL 条件。如图-5 中显示了 x 站位上垂直排列的网格点。点 1 处的小扰动沿着该点处的两条特征线向外传播。特征线就是流动的马赫线，它和
14、流线的夹角就是马赫角 。假设点 1 处流线与 x 轴的夹角为 ，那么左行和右行马赫波与 x 轴的夹角分别是 + 和 -。图-5 只给出了点 1 处的左行马赫波。设有一条通过点 2 的水平线，点 1 的左行特征线与水平线相交与点 a。于是，点 a 和点 2 之间的水平距离 (x)1为35()1= tan(+)1根据点 2 处的 CFL 条件，为了稳定性，推进步长x 的值不应该大雨x1；因而点 2 和点 a 之间的距离应小于或至多等于，声波从点 1 传到与点2 同样高的位置所需的距离。对于点3 处的右行马赫波，也有类似的结果，设它与过点 2 的水平线交于点 b。点b 和点 2 之间的水平距离 为(
15、)336()3= tan()3为保证点 2 处沿流向推进计算稳定性，步长 的值不能大于 和 ()1两者中较小的一个。将这种分析应用于 x0处垂直排列的所有网格点，就()3能给出 x0处沿流向推进的步长x，表达式为：37= |tan()|上式中 是 x=x0 处垂直排列的所有网格点上 的绝|tan()| tan()对中最大的。由于式(20)和(21)定义的坐标变换给出的 =x，那么图-4b 所示的计算区域中，沿流向推进的步长为=x 38南昌航空大学飞行器工程学院课程设计式中x 由式（ 37）确定。联立式（ 37）和式（38 ） ，并引入柯朗数 C，得到 应满足的稳定条件 39= |tan()|C
16、FL 条件要求上式中的 C 1。第三部分 最终结果用空间推进计算得到的结果如下：图-6 列出的是 x=8.9375m，11.375m，23.5625m，31.6875m ，47.7375m和 64.1875m 处的 x 方向速度分量 u 对纵坐标 y 的函数图。图 -6 速度 u 对纵坐标 y 的函数图-7 列出的是 x=8.9375m，11.375m，23.5625m，31.6875m ，47.7375m 和64.1875m 处的压力 p 对纵坐标 y 的函数图。南昌航空大学飞行器工程学院课程设计图-7 压力 p 对纵坐标 y 的函数图-8 列出的是 x=8.9375m，11.375m，23
17、.5625m，31.6875m ，47.7375m 和64.1875m 处的温度 T 对纵坐标 y 的函数图。图-8 温度 T 对纵坐标 y 的函数图-9 列出的是 x=8.9375m，11.375m，23.5625m，31.6875m ，47.7375m 和64.1875m 处的马赫数 Ma 对纵坐标 y 的函数图。南昌航空大学飞行器工程学院课程设计图 -9 马赫数 Ma 对纵坐标 y 的函数第四部分 小结通过这次课程设计，对空间推进的原理有了较深的理解。空间推进用的是定常流的守恒方程组。根据求解区域的形状，使用贴体坐标系。练习了网格生成的某些方法，并使用变换后的控制方程组。还用到了扑捉波的技术。捕捉波应该使用守恒式的控制方程组，还要添加适当的人工粘性，使得解更加光滑。人工粘性主要用在扩张角附近，因为扩张角的定点是数学额上的奇点。最后，对于无粘流的壁面，使用了阿比特的数值边界条件。这种边界条件用到了局部的普朗特迈耶稀疏波，以便将速度转到与壁面平行的方向。计算这个局域内的数值时，我们采用的是 MATLAB 辅助软件，使得计算简便快速。在这次课程设计中了解了 MATLAB 的基本应用。对以前掌握不好的知识点有了更全面的了解和更深层次的掌握，同时对一些以前不熟悉或陌生的知识也有了一定的了解，使得自己的知识面得到了进一步的扩展。
