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

-A.1-一维Riemann问题数值解与计算程序一维Riemann问题，即激波管问题，是一个典型的一维可压缩无黏气体动力学问题，并有解析解。对它采用二阶精度MacCormack两步差分格式进行数值求解。同时，为了初学者入门和练习方便,这里给出了用C语言和Fortran77编写的计算一维Riemann问题的计算程序，供大家学习参考。A-1利用MacCormack两步差分格式求解一维Riemann问题1.一维Riemann问题一维Riemann问题实际上就是激波管问题。激波管是一根两端封闭、内部充满气体的直管，如图A.1所示。在直管中由一薄膜将激波管隔开，在薄膜两侧充有均匀理想气体（可以是同一种气体，也可以是不同种气体），薄膜两侧气体的压力、密度不同。当0t时，气体处于静止状态。当0t时，薄膜瞬时突然破裂，气体从高压端冲向低压端，同时在管内形成激波、稀疏波和接触间断等复杂波系。2.基本方程组、初始条件和边界条件设气体是理想气体。一维Riemann问题在数学上可以用一维可压缩无黏气体Euler方程组来描述。在直角坐标系下量纲为一的一维Euler方程组为：,11xtx0uf(A.1)其中2,()uuupEEpuuf(A.2)这里、u、p、E分别是流体的密度、速度、压力和单位体积总能。理想气体状态方程：221112peEuv（A.3）初始条件：1111,0,1up；2220.125,0,0.1up。图A.1激波管问题示意图-A.2-边界条件：1x和1x处为自由输出条件，01uu,1NNuu。3.二阶精度MacCormack差分格式MacCormack两步差分格式：121111122211122nnnnjjjjnnnnnjjjjjrruuffuuuff(A.4)其中trx。计算实践表明，MacCormack两步差分格式不能抑制激波附近非物理振荡。因此在计算激波时，必须采用人工黏性滤波方法：,,1,,1,122nnnnnijijijijijuuuuu(A.5)为了在激波附近人工黏性起作用，而在光滑区域人工黏性为零，需要引入一个与密度（或者压力）相关的开关函数：1111iiiiiiii(A.6)由式(A.6)可以看出，在光滑区域，密度变化很缓，因此值也很小；而在激波附近密度变化很陡，值就很大。带有开关函数的前置人工黏性滤波方法为：,,1,,1,122nnnnnijijijijijuuuuu(A.7)其中参数往往需要通过实际试算来确定，也可采用线性近似方法得到：||||1tataxx(A.8)由于声速不会超过3，所以取||3a，在本计算中取0.25。4.计算结果分析计算分别采用标准的C语言和Fortran77语言编写程序。计算中网格数取1000，计算总时间为0.4T。计算得到在0.4T时刻的密度、速度和压力分布如图A.2（C语言计算结果）和图A.3（Fortran77语言计算结果）所示。采用两种不同语言编写程序所得到的计算结果完全吻合。从图A.2和图A.3中可以发现，MacCormack两步差分格式能很好地捕捉激波，计算得到的激波面很陡、很窄，计算激波精度是很高的。采用带开关函数的前置-A.3-人工滤波法能消除激波附近的非物理振荡，计算效果很好。从图A.2和图A.3中可以看出通过激波后气体的密度、压力和速度都是增加的；在压力分布中存在第二个台阶，表明在这里存在一个接触间断，在接触间断两侧压力是有间断的，而密度和速度是相等的。这个计算结果正确地反映了一维Riemann问题的物理特性，并被激波管实验所验证。A-2一维Riemann问题数值计算源程序1.C语言源程序//MacCormack1D.cpp:定义控制台应用程序的入口点。///*-----------------------------------------------------------------------------------------------------*利用MacCormack差分格式求解一维激波管问题（C语言版本）*-------------------------------------------------------------------------------------------------------*///#includestdafx.h#includestdio.h#includestdlib.h#includemath.h#defineGAMA1.4//气体常数#definePI3.141592654#defineL2.0//计算区域图A.2采用C语言程序得到的一维Riemann问题密度、速度和压力分布图A.3采用Fortran77语言程序得到的一维Riemann问题密度、速度和压力分布-A.4-#defineTT0.4//总时间#defineSf0.8//时间步长因子#defineJ1000//网格数//全局变量doubleU[J+2][3],Uf[J+2][3],Ef[J+2][3];/*-------------------------------------------------------计算时间步长入口:U，当前物理量，dx，网格宽度；返回:时间步长。---------------------------------------------------------*/doubleCFL(doubleU[J+2][3],doubledx){inti;doublemaxvel,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);if(velmaxvel)maxvel=vel;}returnSf*dx/maxvel;}/*-------------------------------------------------------初始化入口:无；出口:U，已经给定的初始值，dx,网格宽度。---------------------------------------------------------*/voidInit(doubleU[J+2][3],double&dx){inti;doublerou1=1.0,u1=0.0,p1=1.0;//初始条件doublerou2=0.125,u2=0.0,p2=0.1;dx=L/J;for(i=0;i=J/2;i++)-A.5-{U[i][0]=rou1;U[i][1]=rou1*u1;U[i][2]=p1/(GAMA-1)+rou1*u1*u1/2;}for(i=J/2+1;i=J+1;i++){U[i][0]=rou2;U[i][1]=rou2*u2;U[i][2]=p2/(GAMA-1)+rou2*u2*u2/2;}}/*-------------------------------------------------------边界条件入口:dx，网格宽度；出口:U，已经给定的边界。---------------------------------------------------------*/voidbound(doubleU[J+2][3],doubledx){intk;//左边界for(k=0;k3;k++)U[0][k]=U[1][k];//右边界for(k=0;k3;k++)U[J+1][k]=U[J][k];}/*-------------------------------------------------------根据U计算E入口:U，当前U矢量；出口:E，计算得到的E矢量，U、E的定义见Euler方程组。---------------------------------------------------------*/voidU2E(doubleU[3],doubleE[3]){doubleu,p;u=U[1]/U[0];p=(GAMA-1)*(U[2]-0.5*U[1]*U[1]/U[0]);E[0]=U[1];E[1]=U[0]*u*u+p;E[2]=(U[2]+p)*u;}-A.6-/*-------------------------------------------------------一维MacCormack差分格式求解器入口:U，上一时刻的U矢量，Uf、Ef，临时变量，dx，网格宽度，dt,时间步长；出口:U，计算得到的当前时刻U矢量。---------------------------------------------------------*/voidMacCormack_1DSolver(doubleU[J+2][3],doubleUf[J+2][3],doubleEf[J+2][3],doubledx,doubledt){inti,k;doubler,nu,q;r=dt/dx;nu=0.25;for(i=1;i=J;i++){q=fabs(fabs(U[i+1][0]-U[i][0])-fabs(U[i][0]-U[i-1][0]))/(fabs(U[i+1][0]-U[i][0])+fabs(U[i][0]-U[i-1][0])+1e-100);//开关函数for(k=0;k3;k++)Ef[i][k]=U[i][k]+0.5*nu*q*(U[i+1][k]-2*U[i][k]+U[i-1][k]);//人工黏性项}for(k=0;k3;k++)for(i=1;i=J;i++)U[i][k]=Ef[i][k];for(i=0;i=J+1;i++)U2E(U[i],Ef[i]);for(i=0;i=J;i++)for(k=0;k3;k++)Uf[i][k]=U[i][k]-r*(Ef[i+1][k]-Ef[i][k]);//U(n+1/2)(i+1/2)for(i=0;i=J;i++)U2E(Uf[i],Ef[i]);//E(n+1/2)(i+1/2)for(i=1;i=J;i++)for(k=0;k3;k++)U[i][k]=0.5*(U[i][k]+Uf[i][k])-0.5*r*(Ef[i][k]-Ef[i-1][k]);//U(n+1)(i)}/*-------------------------------------------------------输出结果,用Origin数据格式画图入口:U，当前时刻U矢量，dx，网格宽度；出口:无。-A.7----------------------------------------------------------*/voidOutput(doubleU[J+2][3],doubledx){inti;FILE*fp;doublerou,u,p;fp=fopen(result.txt,w);for(i=0;i=J+1;i++){rou=U[i][0];u=U[i][1]/rou;p=(GAMA-1)*(U[i][2]-0.5*U[i][0]*u*u);fprintf(fp,%20f%20.10e%20.10e%20.10e%20.10e\n,i*dx,rou,u,p,U[i][2]);}fclose(fp);}/*-------------------------------------------------------主函数入口:无；出口:无。----------------