一维黎曼问题数值解与计算程序.docVIP

实用标准 PAGE 文档 一维问题数值解与计算程序 一维问题，即激波管问题，是一个典型的一维可压缩无黏气体动力学问题，并有 解析解。对它采用二阶精度两步差分格式进行数值求解。同时，为了初学者入门和练习方便,这里给出了用语言和编写的计算一维问题的计算程序，供大家学习参考。 A-1利用两步差分格式求解一维问题 1.一维问题 图A.1 激波管问题示意图一维问题实际上就是激波管问题。激波管是一根两端封闭、内部充满气体的直管，如图A.1所示。在直管中由一薄膜将激波管隔开，在薄膜两侧充有均匀理想气体（可以是同一种气体，也可以是不同种气体），薄膜两侧气体的压力、密度不同。当时，气体处于静止状态。当时，薄膜瞬时突然破裂，气体从高压端冲向低压端，同时在管内形成激波、稀疏波和接触间断等复杂波系。 图A.1 激波管问题示意图 2.基本方程组、初始条件和边界条件 设气体是理想气体。一维问题在数学上可以用一维可压缩无黏气体方程组来描述。在直角坐标系下量纲为一的一维方程组为： (A.1) 其中 (A.2) 这里、、、分别是流体的密度、速度、压力和单位体积总能。理想气体状态方程： （A.3） 初始条件：；。 边界条件：和处为自由输出条件，,。 3.二阶精度差分格式 两步差分格式： (A.4) 其中。计算实践表明，两步差分格式不能抑制激波附近非物 理振荡。因此在计算激波时，必须采用人工黏性滤波方法： (A.5) 为了在激波附近人工黏性起作用，而在光滑区域人工黏性为零，需要引入一个与密度（或者压力）相关的开关函数： (A.6) 由式(A.6) 可以看出，在光滑区域，密度变化很缓，因此值也很小；而在激波附近密度变化很陡，值就很大。带有开关函数的前置人工黏性滤波方法为： (A.7) 其中参数往往需要通过实际试算来确定，也可采用线性近似方法得到： (A.8) 由于声速不会超过3，所以取，在本计算中取。 4.计算结果分析 计算分别采用标准的语言和语言编写程序。计算中网格数取，计算总时间为。计算得到在时刻的密度、速度和压力分布如图A.2（语言计算结果）和图A.3（语言计算结果）所示。采用两种不同语言编写程序所得到的计算结果完全吻合。 从图A.2和图A.3中可以发现，两步差分格式能很好地捕捉激波，计算得到的激波面很陡、很窄，计算激波精度是很高的。采用带开关函数的前置人工滤波法能消除激波附近的非物理振荡，计算效果很好。 从图A.2和图A.3中可以看出通过激波后气体的密度、压力和速度都是增加的；在压力分布中存在第二个台阶，表明在这里存在一个接触间断，在接触间断两侧压力是有间断的，而密度和速度是相等的。这个计算结果正确地反映了一维问题的物理特性，并被激波管实验所验证。 图A.2采用语言程序得到的一 图A.2采用语言程序得到的一维问题密度、速度和压力分布 图A.3采用语言程序得到的一维问题密度、速度和压力分布 A-2 一维问题数值计算源程序 1. 语言源程序 // MacCormack1D.cpp : 定义控制台应用程序的入口点。 // /* *利用差分格式求解一维激波管问题（语言版本） * */ //#include stdafx.h #include stdio.h #include stdlib.h #include math.h #define GAMA 1.4//气体常数 #define PI 3.141592654 #define L 2.0//计算区域 #define TT 0.4//总时间 #define Sf 0.8//时间步长因子 #define J 1000//网格数 //全局变量 double U[J+2][3],Uf[J+2][3],Ef[J+2][3]; /* 计算时间步长 入口: U，当前物理量，dx，网格宽度； 返回: 时间步长。 */ double CFL(double U[J+2][3],double dx) { int i; double maxvel,p,u,vel; maxvel=1e-100; for(i=1;i=J;i++) { u=U[i][1]/U[i][0]; p=(GAMA-1)*(U[i][2]-0.5*U[i][0]*u*u); vel=sqrt(GAMA*p/U[i][0])+fabs(u);