一维黎曼问题数值解与计算程序Word格式文档下载.docx
《一维黎曼问题数值解与计算程序Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《一维黎曼问题数值解与计算程序Word格式文档下载.docx（15页珍藏版）》请在冰豆网上搜索。
E」
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<
stdio.h>
stdlib.h>
math.h>
#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）
doublerou1=1.0,u1=0.0,p1=1.0;
//初始条件
doublerou2=0.125,u2=0.0,p2=0.1;
dx=L/J;
for（i=0;
=J/2;
U[i][0]=rou1;
U[i][1]=rou1*u1;
U[i][2]=p1/（GAMA-1）+rou1*u1*u1/2;
for（i=J/2+1;
=J+1;
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++）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;
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）;
//开关函数
k++）
Ef[i][k]=u[i][k]+0.5*nu*q*（u[i+1][k]-2*u[i][k]+u[i-1][k]）;
//人工黏性项}
k++）for（i=1;
i++）u[i][k]=Ef[i][k];
i++）u2E（u[i],Ef[i]）;
k++）uf[i][k]=u[i][k]-r*（Ef[i+1][k]-Ef[i][k]）;
//u（n+1/2）（i+1/2）
i++）u2E（uf[i],Ef[i]）;
//E（n+1/2）（i+1/2）
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）
FILE*fp;
doublerou,u,p;
fp=fopen（"
result.txt"
"
w"
）;
rou=U[i][0];
u=U[i][1]/rou;
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<
TT）
{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:
2）dimensionEf（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:
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：
初始条件
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
U，已经给定边界。
subroutinebound（U,dx）implicitdoubleprecision（a-h,o-z）common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU（0：
左边界
do30k=0,2
U（0,k）=U（1,k）
30continue
右边界
do31k=0,2
U（J+1,k）=U（J,k）
31continue
根据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：
2）,E（0：
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
一维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:
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
入口:
U，当前时刻U矢量，
dx，网格宽度；
出口:
subroutineOutput（U,dx）
implicitdoubleprecision（a-h,o-z）common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU（0:
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
升级会员 