《可压缩准一维管道无粘流动_欧拉方程_MacCormack概要.docx》由会员分享，可在线阅读，更多相关《可压缩准一维管道无粘流动_欧拉方程_MacCormack概要.docx（11页珍藏版）》请在三一文库上搜索。
1、可压缩准一维无粘管道流动摘 要 本题利用一维欧拉方程求解可压缩一维无粘管道流动，并针对出口不同的条件，出口亚音速和出口超音速两种不同条件，将流道进行网格划分，利用 MacCormack 进行差分求解，得到管道内的总压、马赫数、总始、内能的分布，并给出计算过程中残 差收敛的过程。关键词 准一维；欧拉方程；MacCormack1问题提出如Figure 1所示，流动为流经变截面流道内、理想、不可压、定常、平面的流动。设进出口截面速度均匀分布。来流条件为M/ = 1. 5,1.2218kg/m3 ,P8 = 47892.4Pa。计算出口超音速和亚音速（4讥=119m/s）时，这个准一维流动的压强和马赫
2、数分布，并给出残差收敛过程Figure 1准一维管道示意图本题的分析思路：首先，建立数学模型，之后运用边界条件定解。接下来将模型问 题转化为数值求解问题，通过网格划分、差分离散，并运用 matlab编程进行求解。 2模型建立选定控制体，根据积分型方程导出该一维流动的欧拉方程。控制体如Figure 1中红框所示，得到连续性、动量和能量方程如（1）、（2）和（3）二（：A）dx =（：uA）dx=0二 t二 xCC一（PuA）dx （:?u2 A）dx - - pAdx二 t：xC个Id一 （： EA）dx （ ： EuA）dx =一 一（puA）dx二t二xdx以上方程对于任意控制体均成立，因此
3、可得到如下微分方程，即一维可压缩无粘流 动的欧拉方程紧凑形式(4)二至二S.:t:x式中,一 p1Q =A Pu小2+ p，S=dXL PHu(5)11其中A = 1.398 + 0.347tanh(0.8x-4.0)。对于此处的欧拉方程，有三个方程组，但是未知量分别为p, P,u, H,E共5个，因此为了方程组的封闭需要补充两个方程。补充方程如下:1 9p= 一1 ： Eu2,=1.4(6)H=E -至此，已实现方程组的封闭，可以进行求解。3边界条件及初始条件进口边界条件：进口为题中给定条件为乂小=1.5 , Pa = 1.2218kg/m3 ,Pg = 47892.4Pa。当地音速按公式a
4、 = 疗,则进口速度为uguaMg,总始和总内能 按照式(4)和式(5)进行求解。出口边界条件：出口边界条件分别为超音速和亚音速。超音速则不给出口条件，亚音速条件按题中给定u0ut=119m/s,在。初始条件：其余节点初始条件假设与进口边界条件相同。4网格划分由于本题所需求解的是一维管道流动问题，在同一截面位置，各个参数均匀分布， 因此只需要沿着管道周向进行网格的划分即可，划分500个节点，空间步长为Ax= 0.02 o5模型的差分离散MacCormack格式广泛应用于求解流动方程， 是一种显示预估/校正二步格式，对于 求解一维欧拉方程其离散格式为:n1n - t n nnQ i 1 - 2Q
5、i Q i jQiQi-(Fi1_ Fi)ttSi:(1 ct2 (8)xx). ,n： ；1n 1n -1n 1 1nn-1：tn-1n-1n_1Qi i -2Qi - QiQ i =Qi Qi -(Fi - Fi ) l tSi2(9)2xx两个方程中均加入了人工粘性进行修正，取人工粘性系数3=0.10.2 ,本题计算时取 a=0.15 0在求解过程中，方程求出的Qin*是Q在第n+1时间层的临时估算值，而第n+1时间层的Q最终值由校正方程(8)提供。需要注意到对空间进行差分离散时，预估步中采 用的是向前差分，而在校正步中采用的是向后差分。每步迭代过程中，先将连续性方程、动量方程和能量方程
6、进行预估迭代,之后将能量方程进行校正迭代,得到(PAE)的新值,根据(3)式求出各点的压力新值，带入动量方程，并进行校正迭代，得到新的(PAu),再将其带入连续性方程，并进行校正迭代，得到新的(PA),最终得到一组新的p, P,u,H,E 值。为了使算法稳定，时间步长取法为：x x ：t=0.5 (10)定义收敛残差为相邻两次迭代过程所得的压力列向量之差的1范数：err =max| pin 1 - pin |(11)计算时取err 10为收敛条件。6计算结果及分析(1)出口为超音速条件时，计算结果如 Figure 2出口超音速条件下计算结果所示(b)马赫数沿轴向分布(a)压力沿轴向分布(c)
7、pv 沿轴向分布(d)残差收敛曲线Figure 2出口超音速条件下计算结果从上述图中可以观察得知，当入口超音、出口超音时，流道内的压强p以及马赫数Ma变化较为平滑，符合超音流在扩张管中降压增速的物理事实，且内部不发生激波现象，为等嫡过程(c图中pv，的波动在1以内，可以认为不变)；但同时也注意到，x = 10处压强和马赫数的突变是因为出口给定条件和计算得到的出口不一致， 从物理上来说就 是出口出现了激波。收敛残差在迭代次数2000次之后开始呈减小趋势，经过4300步迭 代，达到判定收敛标准。根据等嫡流动变截面管道流量公式(12),结合进口条件，可求得此题理论解1 工m = A, p0：0Ma(
8、1 Ma2) 2(12)P0、%为总参数，等嫡流动中为常数AMa(1 VMa2)O=const(13)据(13)求出马赫数理论分布，如Figure 3所示，除了 x=10处的突变，其余完全重合马赫数分布曲线2.8rJ:1理论数值4/1/一2.62.42.221.81.61.4 012345678910流向位置XFigure 3 Ma理论和数值解对比(2)出口为亚音速条件时，4讥=119m/s,计算结果如所示:BIJLwd出募(a)压力沿轴向分布(b)马赫数沿轴向分布pF分布曲线3E流向位置区(c) pv 沿轴向分布(d)残差收敛曲线图1出口为亚音速条件下计算结果从上述结果可以看出，当入口超音、
9、出口亚音时，流场内的压力、马赫数均在流场的 中间某一位置发生了剧烈的突变现象，即流动过程中产生了激波，位置在 x = 4.94左右。 在激波之前为等端流动，经过激波之后，燧突增，之后仍为等端过程。经过 2000次迭代后 残差开始呈下降趋势，到 87000次迭代后收敛。若无激波，超音流在扩张管中持续加速，不可能实现亚音出口，所以激波是超音进口、 亚音出口条件下的必然产物；通过激波之后的亚音流减速增压，负荷亚音流在扩压管中的 流动特征。如果继续调小出口速度，由于需要有足够的空间留给亚音流减速，那么激波位 置必然提前；同理，若增大出口速度，激波位置会后移。下面给出激波位置的理论计算方法。式(14)为
10、激波前等嫡流动的连续方程，Ain、Main、Masi和人分别为进口面积、进口马赫数、激波前马赫数和激波处面积，同理，有式 (15), Aout、Maout、Mas2 分别为出口面积、出口马赫数和激波后马赫数。d _.L._.L 122(二122(AsMasi(1 Maj) 2(=%Main(1 Ma：) 2( )(14)-12 w 占12 才占A.Mas2(1 -Mas22) 2( )= AoutMaout(1 万Maut2) 2(15)驻正激波前后有(16)错误！未找到引用源。，P02为激波之后的总压1 1 Mas12Mas：二2;(16)Mas12 - 一12根据进出口能量守恒(总始守恒)
11、，可得出口温度和马赫数为流道面积方程为Toutp 二uin -UoutMaout =uut/RToutAs =1.398 0.347tanh(0.8xs -4)(17)(18)(19)根据(14)(15)(16)(17)(18)(19)共 6个方程，可以求出 As, Mas- Mas2、Maut、Tout和 xs共6个未知数。据此，可以做出xs随uout变化的关系，出口速度为119时，理论激波 位置4.954和之前的数值解4.94十分接近，激波位置随出口速度减小而前移，符合预期。7.576.565.554.543.5激波位置与出口速度关系X置位波激OOOOO0X: 119Y: 4.954Oft90100110120130140150出口速度UoutFigure 4 xs随uout变化关系7结论本文基于Euler方程，利用MacCormack格式，数值求解了进口超音，出口超音或 亚音时，准一维扩张管内部流动问题。理论和数值解均显示，超音出口下，管内为等嫡 流动；亚音出口下，管内会出现激波，激波位置与出口速度有关，随其减小而前移。8文件说明main_sub是亚音计算程序，main_sup是超音计算程序。