可压缩二维无粘流动二维欧拉方程有限差分maccormackbump.docxVIP

可压缩二维无粘流动摘要本题利用欧拉方程求解可压缩二维无粘流动，并将其与Numeca Fine/Turbo的计算结果对比。流道由上平板固壁和带有凸起的下固壁组成，进口给定总温、总压和速度方向，出口给定压力。自编代码求解时，基于有限差分方法，利用MacCormack格式对控制方程进行离散，根据黎曼不变量和边界条件由内层网格数据外推获得边界数据。文中给出了计算收敛残差历史、密度、速度、压力、马赫数和熵分布，并将其和Numeca计算结果对比，分析自编代码计算结果的合理性和误差来源。关键词二维；欧拉方程；有限差分；MacCormack；Bump1 问题提出该问题是经典的Bump计算问题[1]，如图1所示，上壁为平板，下壁带有凸起，均为滑移边界。进口为轴向进气，且给定总参数为和，出口为。图1准一维管道示意图本题的分析思路：首先，建立计算域中的主控方程，然后根据MacCormack格式对方程进行离散，最后通过边界条件和黎曼不变量确定边界数据。收敛条件为相邻时间步的压力差的最大值小于Pa。2模型建立物理域中的主控方程为二维欧拉方程，如式所示。将物理域x-y变换到计算域ξ-η，控制方程变为式，为坐标变换的雅克比行列式，、、和均为物理域坐标对计算域的偏导数。未知物理量为共6个，因此为了方程组的封闭西需要补充两个方程。补充方程如下：至此，已实现方程组的封闭，可以进行求解。3网格划分计算网格已经给定，分析网格数据可知，第一个数字表示网格块的数目，本题中是1；第二到四个数字表示xyz三个方向的节点数，以L1为例，分别是65、17和1，说明网格为长方向有65个节点、宽方向有17个节点的单层网格；剩余的数字均为网格坐标，前65×17个数是所有点的x坐标，中间65×17个数是所有点的y坐标，最后的是z坐标（本题中全是0）。计算时将物理域转换成矩形的计算域ξ-η，取Δξ=Δη=1。5差分离散和边界条件MacCormack格式广泛应用于求解流动方程，本题中其预估和校正格式分别如式和式，式中i和j分别对应η和ξ节点顺序。以上离散格式只能用于计算非边界网格点，边界网格点数据需由相邻非边界点外推得到。判断本问题为二维亚音入流，因此在进口处需要保持相邻两层的黎曼不变量相同，再结合总温、总压和速度方向，计算出；出口给法类似于进口，不过出口只需一个已知量（如压力）即可；在固壁处，让其压力和密度和相邻内层的相等，将内层速度的法向分量去掉，赋于边界，作为边界处速度。6计算结果以L1为例展示计算结果，计算时人工粘性为0，如图2所示。经过约5万步计算达到收敛，在第4万步时计算曾被暂停，优化代码后继续进行。图2L1计算结果由于是亚音无粘流动，因此流动过程等熵。考虑其物理过程，因为几何完全对称，进口给定总问总压和速度方向，出口给定静压，故以中心（x= 0.5）为对称轴的标量对称流动必然是其解。即温度、压力、速度大小和马赫数均以中心对称分布。另外，由伯努利定律可知，在顶部流速大，压强小，计算结果验证了这一点。6仿真结果本题的仿真软件使用的是Numeca Fine/Turbo。图3为计算1000步，残差到10-5的仿真结果，各参量的分布和数值范围与图2基本一致，说明自编代码数值求解的正确性，验证了本题流动的对称性和最高点高速低压的性质。图3仿真结果7结论本文使用MacCormack格式进行有限差分求解，计算结果与商业软件仿真结果一致，表明本题中二维无粘可压缩管流是一种标量物理量对称分布的流动。编写代码时，尤其需要注意的是边界网格的数据需要遵循黎曼不变量进行外推，否则计算无法收敛。8附录文件说明mainEx2是计算程序，matlab工作空间的计算数据全部保存在prj2DataL1e6中，数据为已收敛的结果。参考资料Rizzi AandViviand H. Numerical Methods for theComputation ofinviscid Transonic Flowswith Shock Waves. Springer Fachmedien Wiesbaden GmbH; 1981, p.52.