一维欧拉方程组的warming-beam差分求解
小弟的一个求解一维rieman问题数字解的程序，收敛不了，求教各位帮我看看这个程序。。。。。
//#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define GAMA 1.4//气体常数
#define PI 3.141592654
#define L 2.0//计算区域
#define TT 0.4//总时间
#define Sf 1//时间步长因子
#define J 1000//网格数
//全局变量
double U[J+2][3],Uf[J+2][3],F1[J+2][3],F2[J+2][3], Ff1[J+2][3], Ff2[J+2][3];
/*-------------------------------------------------------
计算时间步长
入口: U，当前物理量，dx，网格宽度；
返回: 时间步长。
---------------------------------------------------------*/
double CFL(double U[J+2][3],double dx)
{
int i;
double maxvel,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;
}
return Sf*dx/maxvel;
}
/*-------------------------------------------------------
初始化
入口: 无；
出口: U， 已经给定的初始值，
dx, 网格宽度。
---------------------------------------------------------*/
void Init(double U[J+2][3],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;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， 已经给定的边界。
---------------------------------------------------------*/
void bound(double U[J+2][3],double dx)
{
int k;
//左边界
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方程组。
---------------------------------------------------------*/
void U2E(double U[3],double E[3])
{
double u,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;
}
/*-------------------------------------------------------
根据Uf计算F
入口: Uf， 当前U矢量；
出口: F， 计算得到的F矢量，
Uf、F的定义见warming-beam查分格式，及主函数中的定义,计算中间变量。
---------------------------------------------------------*/
void U3F(double U[3],double f1[3],double f2[3])
{
double a,p,u;
u=U[1]/U[0];
p=(GAMA-1)*(U[2]-0.5*U[1]*U[1]/U[0]);
a=GAMA*p/U[0];
if (fabs(u)<a)
{
f1[0]=U[0]/2/GAMA*((2*GAMA-1)*u+a);
f1[1]=U[0]/2/GAMA*(2*(GAMA-1)*u*u+a);
f1[2]=U[0]/2/GAMA*((GAMA-1)*(u*u*u)+0.5*(u+a)*(u+a)*(u+a)+(3-GAMA)*(u+a)*a*a/2/(GAMA-1));
f2[0]=U[0]/2/GAMA*(u-a);
f2[1]=U[0]/2/GAMA*(u-a)*(u-a);
f2[2]=U[0]/2/GAMA*(0.5*(u-a)*(u-a)*(u-a)+(3-GAMA)*(u-a)*a*a/2/(GAMA-1));
}
else
{
U2E(U,f1);
f2[0]=1e-100;f2[1]=1e-100;f2[2]=1e-100;
}
}
/*-------------------------------------------------------
一维 差分格式求解器
入口: U， 上一时刻的U矢量，Uf、F,Ff，临时变量，
dx，网格宽度，dt, 时间步长；
出口: U， 计算得到的当前时刻U矢量。
---------------------------------------------------------*/
void warmingbeam_1DSolver(double U[J+2][3],double Uf[J+2][3],double Ff1[J+2][3],
double Ff2[J+2][3],double F1[J+2][3],double F2[J+2][3],double dx,double dt)
{
int i,k;
double r,p,a,u;
r=dt/dx;
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);
a=GAMA*p/U[i][0];
if (fabs(u)<a)
{
F1[i][0]=U[i][0]/2/GAMA*((2*GAMA-1)*u+a);
F1[i][1]=U[i][0]/2/GAMA*(2*(GAMA-1)*u*u+a);
F1[i][2]=U[i][0]/2/GAMA*((GAMA-1)*u*u*u+0.5*(u+a)*(u+a)*(u+a)+(3-GAMA)*(u+a)*a*a/2/(GAMA-1));
F2[i][0]=U[i][0]/2/GAMA*(u-a);
F2[i][1]=U[i][0]/2/GAMA*(u-a)*(u-a)*(u-a);
F2[i][2]=U[i][0]/2/GAMA*(0.5*(u-a)*(u-a)*(u-a)+(3-GAMA)*(u-a)*a*a/2/(GAMA-1));
}
else
{
U2E(U[i],F1[i]);
F2[i][0]=1e-100;F2[i][1]=1e-100;F2[i][2]=1e-100;
}
}
for(i=0;i<=J;i++)
{ for(k=0;k<3;k++)
{Uf[i][k]=U[i][k]-r*(F1[i][k]-F1[i-1][k]+F2[i+1][k]-F2[i][k]);}
U3F(Uf[i],Ff1[i],Ff2[i]);
}
for(i=2;i<=J;i++)//这里需要注意。。
for(k=0;k<3;k++)
U[i][k]=0.5*(U[i][k]+Uf[i][k])-0.5*r*((F1[i][k]-2*F1[i-1][k]+F1[i-2][k])+(Ff1[i][k]-Ff1[i-1][k])+
(F2[i][k]-2*F2[i-1][k]+F2[i-2][k])+(Ff2[i][k]-Ff2[i-1][k]));
}
/*-------------------------------------------------------
输出结果, 用 数据格式画图
入口: U， 当前时刻U矢量，dx， 网格宽度；
出口: 无。
---------------------------------------------------------*/
void Output(double U[J+2][3],double dx)
{
int i;
FILE *fp;
double rou,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);
}
/*-------------------------------------------------------
主函数
入口: 无；
出口: 无。
---------------------------------------------------------*/
void main()
{
double T,dx,dt;
Init(U,dx);
T=0;
while(T<TT)
{
dt=CFL(U,dx);
T+=dt;
printf("T=%10g dt=%10g\n",T,dt);
warmingbeam_1DSolver(U,Uf,Ff1,Ff2,F1,F2,dx,dt);
bound(U,dx);
}
Output(U,dx);
}
