一维欧拉方程matlab代码,一维欧拉方程组的warming-beam差分求解

已结贴√

问题点数：20 回复次数：5

一维欧拉方程组的warming-beam差分求解

小弟的一个求解一维rieman问题数字解的程序，收敛不了，求教各位帮我看看这个程序。。。

。。

//#include "stdafx.h"

#include

#include

#include

#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)

{

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)