一维黎曼问题数值解与计算程序.docx
《一维黎曼问题数值解与计算程序.docx》由会员分享，可在线阅读，更多相关《一维黎曼问题数值解与计算程序.docx（19页珍藏版）》请在冰豆网上搜索。
一维黎曼问题数值解与计算程序
一维黎曼问题数值解与计算程序
一维
问题数值解与计算程序
一维
问题，即激波管问题，是一个典型的一维可压缩无黏气体动力学问题，并有解析解。
对它采用二阶精度
两步差分格式进行数值求解。
同时，为了初学者入门和练习方便,这里给出了用
语言和
编写的计算一维
问题的计算程序，供大家学习参考。
A-1利用
两步差分格式求解一维
问题
1.一维
问题
一维
问题实际上就是激波管问题。
激波管是一根两端封闭、内部充满气体的直管，如图A.1所示。
在直管中由一薄膜将激波管隔开，在薄膜两侧充有均匀理想气体（可以是同一种气体，也可以是不同种气体），薄膜两侧气体的压力、密度不同。
当
时，气体处于静止状态。
当
时，薄膜瞬时突然破裂，气体从高压端冲向低压端，同时在管内形成激波、稀疏波和接触间断等复杂波系。
2.基本方程组、初始条件和边界条件
设气体是理想气体。
一维
问题在数学上可以用一维可压缩无黏气体
方程组来描述。
在直角坐标系下量纲为一的一维
方程组为：
（A.1）
其中
（A.2）
这里
、
、
、
分别是流体的密度、速度、压力和单位体积总能。
理想气体状态方程：
（A.3）
，计算总时间为
。
计算得到在
时刻的密度、速度和压力分布如图A.2（
语言计算结果）和图A.3（
语言计算结果）所示。
采用两种不同语言编写程序所得到的计算结果完全吻合。
从图A.2和图A.3中可以发现，
两步差分格式能很好地捕捉激波，计算得到的激波面很陡、很窄，计算激波精度是很高的。
采用带开关函数的前置人工滤波法能消除激波附近的非物理振荡，计算效果很好。
从图A.2和图A.3中可以看出通过激波后气体的密度、压力和速度都是增加的；在压力分布中存在第二个台阶，表明在这里存在一个接触间断，在接触间断两侧压力是有间断的，而密度和速度是相等的。
这个计算结果正确地反映了一维
问题的物理特性，并被激波管实验所验证。
A-2一维
问题数值计算源程序
1.
语言源程序
//MacCormack1D.cpp:
定义控制台应用程序的入口点。
//
/*-----------------------------------------------------------------------------------------------------
*利用
差分格式求解一维激波管问题（
语言版本）
*-------------------------------------------------------------------------------------------------------*/
//#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;
}
/*-------------------------------------------------------
一维
差分格式求解器
入口:
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;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）
}
/*-------------------------------------------------------
输出结果,用
数据格式画图
入口:
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.
语言源程序
!
MacCormack1D.for
------------------------------------------------------------------------------------------------------------
!
利用
差分格式求解一维激波管问题（
语言版本）
--------------------------------------------------------------------------------------------------------------*/
programMacCormack1D
implicitdoubleprecision（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.4
PI=3.1415926
!
网格数
J=M
!
计算区域
dL=2.0
!
总时间
TT=0.4
!
时间步长因子
Sf=0.8
callInit（U,dx）
T=0
1dt=CFL（U,dx）
T=T+dt
write（*,*）'T=',T,'dt=',dt
callMacCormack_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,Sf
dimensionU（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,Sf
dimensionU（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）=rou1
U（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*u2
21continue
end
!
-------------------------------------------------------
!
边界条件
!
入口:
dx，网格宽度；
!
出口:
U，已经给定边界。
!
-------------------------------------------------------
subroutinebound（U,dx）
implicitdoubleprecision（a-h,o-z）
common/G_def/GAMA,PI,J,JJ,dL,TT,Sf
dimensionU（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,Sf
dimensionU（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
!
-------------------------------------------------------
!
一维
差分格式求解器
!
入口:
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,Sf
dimensionU（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
!
-------------------------------------------------------
!
输出结果,用
数据格式画图
!
入口:
U，当前时刻U矢量，
!
dx，网格宽度；
!
出口:
无。
!
-------------------------------------------------------
subroutineOutput（U,dx）
implicitdoubleprecision（a-h,o-z）
common/G_def/GAMA,PI,J,JJ,dL,TT,Sf
dimensionU（0:
J+1,0:
2）
open（1,file='result.txt',status='unknown'）
do80i=0,J+1
rou=U（i,0）
uu=U（i,1）/rou
p=（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
------------------------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------------------------
升级会员 