可压缩二维无粘流动_二维_欧拉方程_有限差分_MacCormack_Bump

可压缩二维无粘流动 摘要本题利用欧拉方程求解可压缩二维无粘流动，并将其与Numeca Fine/Turbo的计算结果对比。流道由上平板固壁和带有凸起的下固壁组成，进口给定总温、总压和速度方向，出口给定压力。自编代码求解时，基于有限差分方法，利用MacCormack格式对控制方程进行离散，根据黎曼不变量和边界条件由内层网格数据外推获得边界数据。文中给出了计算收敛残差历史、密度、速度、压力、马赫数和熵分布，并将其和Numeca计算结果对比，分析自编代码计算结果的合理性和误差来源。 关键词二维；欧拉方程；有限差分；MacCormack；Bump 1  问题提出 该问题是经典的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中，数据为已收敛的结果。 参考资料 [1] Rizzi AandViviand H. Numerical Methods for theComputation ofinviscid Transonic Flowswith Shock Waves. Springer Fachmedien Wiesbaden GmbH; 1981, p.52.