l(E+p)u』
这里「、u、p、E分别是流体的密度、速度、压力和单位体积总能。
理想气体
状态方程:
p=-1◎-1E-Ju2v2(A.3)
初始条件:
J=1,5=0,山=1;=0.125,u2=0,p2=0.1。
边界条件:
X--1和X=1处为自由输出条件,UohU—UnhUn/。
3.二阶精度MacCormac差分格式
n-2
MacCormac两步差分格式:
2nn-n
Uj-Uj-rfj-fjj
其中r二卫。
计算实践表明,MacCormac两步差分格式不能抑制激波附近非物
Ax
理振荡。
因此在计算激波时,必须采用人工黏性滤波方法:
(A.5)
需要引入一个与
-nn1nQnn
Ui,j=Ui,j-U「1,j-2嘔-UiJ,j
为了在激波附近人工黏性起作用,而在光滑区域人工黏性为零,
密度(或者压力)相关的开关函数:
由式(A.6)可以看出,在光滑区域,密度变化很缓,因此二值也很小;而在激波附近密度变化很陡,V值就很大。
带有开关函数的前置人工黏性滤波方法为:
(A.7)
_1
Uinj=unj+罗日(ut-2uinj+ulj
其中参数往往需要通过实际试算来确定,也可采用线性近似方法得到:
(A.8)
由于声速不会超过3,所以取|a|=3,在本计算中取=0.25。
4.计算结果分析
计算分别采用标准的C语言和Fortran77语言编写程序。
计算中网格数取1000,计算总时间为T=0.4。
计算得到在T=0.4时刻的密度、速度和压力分布
如图A.2(C语言计算结果)和图A.3(Fortran77语言计算结果)所示。
采用两
种不同语言编写程序所得到的计算结果完全吻合。
从图A.2和图A.3中可以发现,MacCormac两步差分格式能很好地捕捉激波,计算得到的激波面很陡、很窄,计算激波精度是很高的。
采用带开关函数的前置
人工滤波法能消除激波附近的非物理振荡,计算效果很好
从图A.2和图A.3中可以看出通过激波后气体的密度、压力和速度都是增加的;在压力分布中存在第二个台阶,表明在这里存在一个接触间断,在接触间断两侧压力是有间断的,而密度和速度是相等的。
这个计算结果正确地反映了一维Riemann问题的物理特性,并被激波管实验所验证。
A-2一维Riemanr问题数值计算源程序
1.C语言源程序
//MacCormackID.cpp:
定义控制台应用程序的入口点。
//
/*
*利用MacCormac差分格式求解一维激波管问题(C语言版本)
*
*/
〃#include"stdafx.h"
#include
#include
#include#defineGAMA1.4〃气体常数
#definePI3.141592654
#defineL2.0〃计算区域
#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(vel>maxvel)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++)
{
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;k<3;k++)U[0][k]=U[1][k];
//右边界
for(k=0;k<3;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;
}
/*
一维MacCormac差分格式求解器
入口:
u,上一时刻的u矢量,Uf、Ef,临时变量,
dx,网格宽度,dt,时间步长;
出口:
u,计算得到的当前时刻u矢量。
*/
Ef[J+2][3],double
voidMacCormack_1DSolver(doubleu[J+2][3],doubleuf[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;k<3;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;k<3;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;k<3;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;k<3;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,网格宽度;
出口:
无。
*/
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);
}
/*
主函数入口:
无;出口:
无。
*/
voidmain()
{
doubleT,dx,dt;
Init(U,dx);
T=0;while(T
{dt=CFL(U,dx);
T+=dt;
printf("T=%10gdt=%10g\n",T,dt);
MacCormack_1DSolver(U,Uf,Ef,dx,dt);
bound(U,dx);
}
Output(U,dx);
}
2.Fortran7语言源程序
!
MacCormack1D.for
!
利用MacCormac差分格式求解一维激波管问题(Fortran77语言版本)*/programMacCormack1Dimplicitdoubleprecision(a-h,o-z)parameter(M=1000)common/G_def/GAMA,PI,J,JJ,dL,TT,Sf
dimensionU(0:
M+1,0:
2),Uf(0:
M+1,0:
2)dimensionEf(0:
M+1,0:
2)
!
气体常数
GAMA=1.4PI=3.1415926
!
网格数
J=M
!
计算区域
dL=2.0
!
总时间
TT=0.4
!
时间步长因子
Sf=0.8
callInit(U,dx)
T=0
1dt=CFL(U,dx)
T=T+dtwrite(*,*)'T=',T,'dt=',dtcallMacCormack_1D_Solver(U,Uf,Ef,dx,dt)callbound(U,dx)
if(T.lt.TT)goto1
callOutput(U,dx)end
!
!
计算时间步长
!
入口:
U,当前物理量,dx,网格宽度;
!
返回:
时间步长。
!
doubleprecisionfunctionCFL(U,dx)implicitdoubleprecision(a-h,o-z)common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU(0:
J+1,0:
2)
dmaxvel=1e-10
do10i=1,J
uu=U(i,1)/U(i,0)p=(GAMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu)vel=dsqrt(GAMA*p/U(i,0))+dabs(uu)if(vel.gt.dmaxvel)dmaxvel=vel
10continue
CFL=Sf*dx/dmaxvel
end
!
!
初始化
!
入口:
无;
!
出口:
U,已经给定的初始值,dx,网格宽度。
!
subroutineInit(U,dx)implicitdoubleprecision(a-h,o-z)common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU(0:
J+1,0:
2)
!
初始条件
rou1=1.0
u1=0
v1=0
p1=1.0
rou2=0.125
u2=0
v2=0
p2=0.1
dx=dL/J
do20i=0,J/2
U(i,0)=rou1U(i,1)=rou1*u1
U(i,2)=p1/(GAMA-1)+0.5*rou1*u1*u1
20continue
do21i=J/2+1,J+1
U(i,0)=rou2
U(i,1)=rou2*u2
U(i,2)=p2/(GAMA-1)+0.5*rou2*u2*u221continue
end
!
!
边界条件
!
入口:
dx,网格宽度;
!
出口:
U,已经给定边界。
!
subroutinebound(U,dx)implicitdoubleprecision(a-h,o-z)common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU(0:
J+1,0:
2)
!
左边界
do30k=0,2
U(0,k)=U(1,k)
30continue
!
右边界
do31k=0,2
U(J+1,k)=U(J,k)
31continue
end
!
!
根据U计算E
!
入口:
U,当前U矢量;
!
出口:
E,计算得到的E矢量,
!
U、E定义见Euler方程组。
!
subroutineU2E(U,E,is,in)implicitdoubleprecision(a-h,o-z)common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU(0:
J+1,0:
2),E(0:
J+1,0:
2)
do40i=is,in
uu=U(i,1)/U(i,0)p=(GAMA-1)*(U(i,2)
$-0.5*U(i,1)*U(i,1)/U(i,0))
E(i,0)=U(i,1)
E(i,1)=U(i,0)*uu*uu+p
E(i,2)=(U(i,2)+p)*uu
40continue
end
!
!
一维MacCormac差分格式求解器
!
入口:
U,上一时刻U矢量,
!
Uf、Ef,临时变量,
!
dx,网格宽度,dt,,时间步长;
!
出口:
U,计算得到得当前时刻U矢量。
!
subroutineMacCormack_1D_Solver(U,Uf,Ef,dx,dt)implicitdoubleprecision(a-h,o-z)common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU(0:
J+1,0:
2),Uf(0:
J+1,0:
2)dimensionEf(0:
J+1,0:
2)
r=dt/dx
dnu=0.25
do60i=1,J
do60k=0,2
!
开关函数
q=dabs(dabs(U(i+1,0)-U(i,0))-dabs(U(i,0)-U(i-1,0)))
$/(dabs(U(i+1,0)-U(i,0))+dabs(U(i,0)-U(i-1,0))+1e-10)!
人工黏性项
Ef(i,k)=U(i,k)+0.5*dnu*q*(U(i+1,k)-2*U(i,k)+U(i-1,k))
60continue
do61k=0,2
do61i=1,J
U(i,k)=Ef(i,k)
61continue
callU2E(U,Ef,0,J+1)
do63i=0,J
do63k=0,2
!
U(n+1/2)(i+1/2)
Uf(i,k)=U(i,k)-r*(Ef(i+1,k)-Ef(i,k))
63continue
!
E(n+1/2)(i+1/2)
callU2E(Uf,Ef,0,J)
do64i=1,J
do64k=0,2
!
U(n+1)(i)
U(i,k)=0.5*(U(i,k)+Uf(i,k))-0.5*r*(Ef(i,k)-Ef(i-1,k))
64continue
end
!
!
输出结果,用Origin数据格式画图
!
入口:
U,当前时刻U矢量,
!
dx,网格宽度;
!
出口:
无。
!
subroutineOutput(U,dx)
implicitdoubleprecision(a-h,o-z)common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU(0:
J+1,0:
2)
open(1,file='result.txt',status='unknown')
do80i=0,J+1
rou=U(i,0)
uu=U(i,1)/roup=(GAMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu)write(1,81)i*dx,rou,uu,p,U(i,2)
80continue
close
(1)
81format(D20.10,D20.10,D20.10,D20.10,D20.10)end