《一维黎曼问题数值解与计算程序》由会员分享，可在线阅读，更多相关《一维黎曼问题数值解与计算程序（12页珍藏版）》请在金锄头文库上搜索。
1、图A.1激波管问题示意图其中du dfdtdx=0, 1 x 1p、pu 、u =pu,f =pu 2 + pE丿、(E + p)u 丿(A.1)(A.2)一维Riemann问题数值解与计算程序一维Riemann问题，即激波管问题，是一个典型的一维可压缩无黏气体动力 学问题，并有 解析解。对它采用二阶精度MacCormac两步差分格式进行数值求 解。同时，为了初学者入门和练习方便，这里给出了用C语言和Fortran77编写的 计算一维Riemann问题的计算程序，供大家学习参考。A-1利用 MacCormack两步差分格式求解一维Riemam问题1. 一维 Riemann问题一维Riemann
2、问题实际上就是激波管问题。激波管是一根两端封闭、内部充满气体的直管，如图A.1所示。在直管中由一薄膜将激波管隔开，在薄膜两侧充 有均匀理想气体（可以是同一种气体, 也可以是不同种气体），薄膜两侧气体 的压力、密度不同。当t 0时，气体 处于静止状态。当t二0时，薄膜瞬时 突然破裂，气体从高压端冲向低压端, 同时在管内形成激波、稀疏波和接触间 断等复杂波系。2基本方程组、初始条件和边界条件设气体是理想气体。一维Riemann问题在数学上可以用一维可压缩无黏气体Euler方程组来描述。在直角坐标系下量纲为一的一维Euler方程组为:这里p、u、p、e分别是流体的密度、速度、压力和单位体积总能。理想
3、气体状态方程:p =(Y -l)pe = (丫 -1)(U 2 + V 2)(A.3)初始条件：p = 1, u = 0, p = 1 ； p = 0.125, u = 0, p = 0.1。1 1 1 2 2 2边界条件：x = -1和x = 1处为自由输出条件,u = u , u = u01 NN-13.二阶精度MacCormac差分格式MacCormac两步差分格式:Un+ 2 = Un 一 r (f n f n )jjjj-11Un+1 =j 2(丄)1(1丄）Un + Un+ 2一一rf n+2-f n+2#include vstdlib.h#include #define GAMA
4、 1.4/气 体常数#define PI 3.141592654 #define L 2.0计算区域#define TT 0.4总时间#define Sf 0.8/时间步长因子#define J 1000/ 网格数全局变量double UJ+23,UfJ+23,EfJ+23;/*计算时间步长入口： U,当前物理量，dx,网格宽度； 返回：时间步长。*/ double CFL(double UJ+23,double dx) int i;double maxvel,p,u,vel;maxvel=1e-100;for(i=1;iv=J;i+)u=Ui1/Ui0; p=(GAMA-1)*(Ui2-0.
5、5*Ui0*u*u); vel=sqrt(GAMA*p/UiO)+fabs(u);if(velmaxvel)maxvel=vel;return Sf*dx/maxvel;/*初始化入口：无；出口： U，已经给定的初始值， dx,网格宽度。*/ void Init(double UJ+23,double & dx)int i;double rou1=1.0 ,u1=0.0,p1=1.0; /初始条件 double rou2=0.125,u2=0.0,p2=0.1;dx=L/J;for(i=0;iv=J/2;i+)Ui0=roul;Uil=roul*ul; Ui2=p1/(GAMA-1)+rou1
6、*u1*u1/2;for(i=J/2+1;iv=J+1;i+)Ui0=rou2;Ui1=rou2*u2; Ui2=p2/(GAMA-1)+rou2*u2*u2/2;/*边界条件入口： dx,网格宽度； 出口： U,已经给定的边界。 */ void bound(double UJ+23,double dx)int k;左边界 for(k=0;kv3;k+)U0k=U1k;右边界 for(k=0;k3;k+)UJ+1k=UJk;/*根据U计算E入口： U，当前U矢量；出口： E,计算得到的E矢量， U、E的定义见Euler方程组。*/ void U2E(double U3,double E3)do
7、uble u,p;u=U1/U0; p=(GAMA-1)*(U2-0.5*U1*U1/U0);E0=U1;E1=U0*u*u+p; E2=(U2+p)*u;/*一维MacCormack差分格式求解器 入口： U,上一时刻的U矢量，Uf、Ef,临时变量，dx,网格宽度，dt,时间步长；出口： U，计算得到的当前时刻U矢量。*/void MacCormack_lDSolver(double UJ+23,double UfJ+23,double EfJ+23,double dx,double dt)int i,k;double r,nu,q;r=dt/dx;nu=0.25;for(i=1;iv=J;
8、i+)q=fabs(fabs(Ui+10-Ui0)-fabs(Ui0-Ui-10) /(fabs(Ui+10-Ui0)+fabs(Ui0-Ui-10)+1e-100); /开关函数for(k=0;k3;k+)Efik=Uik+0.5*nu*q*(Ui+1k-2*Uik+Ui-1k);/人工黏性项for(k=0;k3;k+)for(i=1;i=J;i+)Uik=Efik;for(i=0;i=J+1;i+)U2E(Ui,Efi);for(i=0;i=J;i+)for(k=0;k3;k+)Ufik=Uik-r*(Efi+1k-Efik); /U(n+1/2)(i+1/2)for(i=0;i=J;i+
9、)U2E(Ufi,Efi); /E(n+1/2)(i+1/2)for(i=1;i=J;i+)for(k=0;k3;k+)Uik=0.5*(Uik+Ufik)-0.5*r*(Efik-Efi-1k); /U(n+1)(i)/*输出结果，用Origin数据格式画图入口： U，当前时刻U矢量，dx，网格宽度；出口：无。*/void Output(double UJ+23,double dx)int i;FILE *fp;double rou,u,p;fp=fopen(result.txt,w);for(i=O;iv=J+l;i+)rou=Ui0;u=Ui1/rou;p=(GAMA-1)*(Ui2-0
10、.5*Ui0*u*u);fprintf(fp,%20f%20.10e%20.10e%20.10e%20.10en,i*dx,rou,u,p,Ui2);fclose(fp);/*主函数入口：无；出口：无。*/void main()double T,dx,dt;Init(U,dx);T=0;while(TvTT)dt=CFL(U,dx);T+=dt;printf(T=%10g dt=%10gn,T,dt);MacCormack_1DSolver(U,Uf,Ef,dx,dt);bound(U,dx);Output(U,dx);2. Fortran77吾言源程序! MacCormacklD.for利用MacCormac差分格式求解一维激波
