可压缩准一维管道无粘流动-欧拉方程-MacCormack

1可压缩准一维无粘管道流动摘要本题利用一维欧拉方程求解可压缩一维无粘管道流动，并针对出口不同的条件，出口亚音速和出口超音速两种不同条件，将流道进行网格划分，利用MacCormack进行差分求解，得到管道内的总压、马赫数、总焓、内能的分布，并给出计算过程中残差收敛的过程。关键词准一维；欧拉方程；MacCormack1问题提出如Figure1所示，流动为流经变截面流道内、理想、不可压、定常、平面的流动。设进出口截面速度均匀分布。来流条件为1.5M，31.2218/kgm，47892.4pPa。计算出口超音速和亚音速(119/outums)时，这个准一维流动的压强和马赫数分布，并给出残差收敛过程。Figure1准一维管道示意图本题的分析思路：首先，建立数学模型，之后运用边界条件定解。接下来将模型问题转化为数值求解问题，通过网格划分、差分离散，并运用matlab编程进行求解。2模型建立选定控制体，根据积分型方程导出该一维流动的欧拉方程。控制体如Figure1中红框所示，得到连续性、动量和能量方程如(1)、(2)和(3)()()0AdxuAdxtx(1)2()()uAdxuAdxpAdxtx(2)()()()dEAdxEuAdxpuAdxtxdx(3)控制体2以上方程对于任意控制体均成立，因此可得到如下微分方程，即一维可压缩无粘流动的欧拉方程紧凑形式(4)txQFS(4)式中，20,,0dAAuAuppdxEHuQFS(5)其中1.3980.347tanh(0.84.0)Ax。对于此处的欧拉方程，有三个方程组，但是未知量分别为,,,,puHE共5个，因此为了方程组的封闭需要补充两个方程。补充方程如下：211,1.42pEu(6)pHE(7)至此，已实现方程组的封闭，可以进行求解。3边界条件及初始条件进口边界条件：进口为题中给定条件为1.5M，31.2218/kgm，47892.4pPa。当地音速按公式pa，则进口速度为uaM，总焓和总内能按照式(4)和式(5)进行求解。出口边界条件：出口边界条件分别为超音速和亚音速。超音速则不给出口条件，亚音速条件按题中给定119/outums，在。初始条件：其余节点初始条件假设与进口边界条件相同。4网格划分由于本题所需求解的是一维管道流动问题，在同一截面位置，各个参数均匀分布，因此只需要沿着管道周向进行网格的划分即可，划分500个节点，空间步长为0.02x。5模型的差分离散MacCormack格式广泛应用于求解流动方程，是一种显示预估/校正二步格式，对于求解一维欧拉方程其离散格式为：3111122()nnnnnnnniiiiiiiitttxxQQQQQFFS(8)11111111111221[()]2nnnnnnnnniiiiiiiiitttxxQQQQQQFFS(9)两个方程中均加入了人工粘性进行修正，取人工粘性系数=0.1~0.2，本题计算时取=0.15。在求解过程中，方程(7)求出的1niQ是Q在第1n时间层的临时估算值，而第1n时间层的Q最终值由校正方程(8)提供。需要注意到对空间进行差分离散时，预估步中采用的是向前差分，而在校正步中采用的是向后差分。每步迭代过程中，先将连续性方程、动量方程和能量方程进行预估迭代，之后将能量方程进行校正迭代，得到()AE的新值，根据(3)式求出各点的压力新值，带入动量方程，并进行校正迭代，得到新的()Au，再将其带入连续性方程，并进行校正迭代，得到新的()A，最终得到一组新的,,,,puHE值。为了使算法稳定，时间步长取法为：0.5xxtuaua(10)定义收敛残差为相邻两次迭代过程所得的压力列向量之差的1范数：1max{||}nniierrpp(11)计算时取810err为收敛条件。6计算结果及分析(1)出口为超音速条件时，计算结果如Figure2出口超音速条件下计算结果所示4(a)压力沿轴向分布(b)马赫数沿轴向分布(c)pv沿轴向分布(d)残差收敛曲线Figure2出口超音速条件下计算结果从上述图中可以观察得知，当入口超音、出口超音时，流道内的压强p以及马赫数Ma变化较为平滑，符合超音流在扩张管中降压增速的物理事实，且内部不发生激波现象，为等熵过程（c图中pv的波动在1以内，可以认为不变）；但同时也注意到，10x处压强和马赫数的突变是因为出口给定条件和计算得到的出口不一致，从物理上来说就是出口出现了激波。收敛残差在迭代次数2000次之后开始呈减小趋势，经过4300步迭代，达到判定收敛标准。根据等熵流动变截面管道流量公式(12)，结合进口条件，可求得此题理论解122(1)001(1)2mApMaMa(12)50p、0为总参数，等熵流动中为常数122(1)1(1)2AMaMaconst(13)据(13)求出马赫数理论分布，如Figure3所示，除了x=10处的突变，其余完全重合。0123456789101.41.61.822.22.42.62.8流向位置xMa马赫数分布曲线理论数值Figure3Ma理论和数值解对比(2)出口为亚音速条件时，119/outums，计算结果如所示：(a)压力沿轴向分布(b)马赫数沿轴向分布6(c)pv沿轴向分布(d)残差收敛曲线图1出口为亚音速条件下计算结果从上述结果可以看出，当入口超音、出口亚音时，流场内的压力、马赫数均在流场的中间某一位置发生了剧烈的突变现象，即流动过程中产生了激波，位置在4.94x左右。在激波之前为等熵流动，经过激波之后，熵突增，之后仍为等熵过程。经过2000次迭代后残差开始呈下降趋势，到87000次迭代后收敛。若无激波，超音流在扩张管中持续加速，不可能实现亚音出口，所以激波是超音进口、亚音出口条件下的必然产物；通过激波之后的亚音流减速增压，负荷亚音流在扩压管中的流动特征。如果继续调小出口速度，由于需要有足够的空间留给亚音流减速，那么激波位置必然提前；同理，若增大出口速度，激波位置会后移。下面给出激波位置的理论计算方法。式(14)为激波前等熵流动的连续方程，inA、inMa、1sMa和sA分别为进口面积、进口马赫数、激波前马赫数和激波处面积，同理，有式(15)，outA、outMa、2sMa分别为出口面积、出口马赫数和激波后马赫数。11222(1)2(1)1111(1)(1)22sssinininAMaMaAMaMa(14)11222(1)2(1)2211(1)(1)22sssoutoutoutAMaMaAMaMa(15)驻正激波前后有(16)错误!未找到引用源。，02p为激波之后的总压21222111212sssMaMaMa(16)7根据进出口能量守恒（总焓守恒），可得出口温度和马赫数为222inoutoutpuupTRc(17)/outoutoutMauRT(18)流道面积方程为1.3980.347tanh(0.8x4)ssA(19)根据(14)(15)(16)(17)(18)(19)共6个方程，可以求出sA，1sMa，2sMa、outMa、outT和xs共6个未知数。据此，可以做出xs随uout变化的关系，出口速度为119时，理论激波位置4.954和之前的数值解4.94十分接近，激波位置随出口速度减小而前移，符合预期。901001101201301401503.544.555.566.577.5X:119Y:4.954激波位置与出口速度关系出口速度Uout激波位置xFigure4xs随uout变化关系7结论本文基于Euler方程，利用MacCormack格式，数值求解了进口超音，出口超音或亚音时，准一维扩张管内部流动问题。理论和数值解均显示，超音出口下，管内为等熵流动；亚音出口下，管内会出现激波，激波位置与出口速度有关，随其减小而前移。88文件说明main_sub是亚音计算程序，main_sup是超音计算程序。