请教：Roe格式求解一维欧拉方程

pittsun

主要步骤在于求解通量。
我的算法完全参照教材，可是在求算Sod问题时，在激波处总是计算发散，不知道为什么。我怀疑是不是程序中的算法有什么问题，在这儿我把求解通量的函数贴在这儿，请高手们帮我看看。
int i,j,k;
// new items in Roe method
double *lambda = new double [3];
double vec[3][3];
double *alpha = new double [3];
double FL[3];
double FR[3];
double sound;
double rho;
double u;
double e[2];
double enthalpy;
double temp_l,temp_r;
double temp1,temp2;
//  calculation process of Flux[]
for(i=0;i<cells+1;i++)
{
//set values
rho = sqrt(U[0]*U[i+1][0]);
temp_l = sqrt(U[0]);
temp_r = sqrt(U[i+1][0]);
u = (temp_l*U[1]+temp_r*U[i+1][1])/(temp_l+temp_r);
e[0] = 0.5*G4*U[2]/U[0];
e[1] = 0.5*G4*U[i+1][2]/U[i+1][0];
temp1 = temp_r*(e[1]+0.5*U[i+1][1]*U[i+1][1]+U[i+1][2]/U[i+1][0]);
temp2 = temp_l*(e[0]+0.5*U[1]*U[1]+U[2]/U[0]);
enthalpy = (temp1+temp2)/(temp_l+temp_r);
sound = sqrt(G8*(enthalpy-0.5*u*u));
lambda[0] = u-sound;
lambda[1] = u;
lambda[2] = u+sound;
vec[0][0] = 1.0;              vec[0][1] = 1.0;      vec[0][2] = 1.0;
vec[1][0] = lambda[0];        vec[1][1] = u;        vec[1][2] = lambda[2];
vec[2][0] = enthalpy-u*sound; vec[2][1] = 0.5*u*u;  vec[2][2] = enthalpy+u*sound;
temp1 = (U[i+1][2]-U[2])-rho*sound*(U[i+1][1]-U[1]);
temp2 = (U[i+1][2]-U[2])+rho*sound*(U[i+1][1]-U[1]);
alpha[0] = 0.5*temp1/sound/sound;
alpha[1] = (U[i+1][0]-U[0])-(U[i+1][2]-U[2])/sound/sound;
alpha[2] = 0.5*temp2/sound/sound;
FL[0] = U[0]*U[1];
FL[1] = FL[0]*U[1]+U[2];
FL[2] = U[1]*(0.5*G3*U[2]+0.5*U[0]*U[1]*U[1]);
FR[0] = U[i+1][0]*U[i+1][1];
FR[1] = FL[0]*U[i+1][1]+U[i+1][2];
FR[2] = U[i+1][1]*(0.5*G3*U[i+1][2]+0.5*U[i+1][0]*U[i+1][1]*U[i+1][1]);

for(j=0;j<3;j++)
{
double temp = 0;
for(k=0;k<3;k++)
{
temp += fabs(lambda[k])*alpha[k]*vec[k][j];
}
Flux[j] = 0.5*(FL[j]+FR[j])-0.5*temp;
}
}
程序中U数组为基本变量数组，依次为密度，速度，压力。
多谢先。[br][br][以下内容由 pittsun 在 2007年12月30日 00:31pm 时添加] [br]
// G8 = gamma-1.0
// G4 = 2.0/(gamma-1.0); [br][br][以下内容由 pittsun 在 2007年12月30日 00:32pm 时添加] [br]
// G3 = 2.0*gamma/(gamma-1.0);

zongxing85

你看的那本书啊？这个程序完整吗？

zongxing85

你把熵条件加进去，再试试看看可不可以

pittsun

下面引用由zongxing85在 2008/01/15 04:47pm 发表的内容：
你把熵条件加进去，再试试看看可不可以

不是熵条件的问题，是我程序里一个i,j的顺序弄错了，现在全改好了，多谢哈~

shirazbj

好