SOS,计算程序需要修改,通过一个时间的控制,求一维轴上201个点的速度压力密度,能运行,但是计算结果溢出,求修改
程序代码:#include<stdio.h>
#include<math.h>
float compare(float *b,int m) /*比较函数声明*/
{
float max;
int j;
max=*b;
for(j=0;j<m;j++)
{
if(max<b[j])
max=b[j];
else
continue;
}
return max;
}
void main()
{
int i=0,R;
float U1new[201],U2new[201],U3new[201],
U1[201],U2[201],U3[201],
F1[201],F2[201],F3[201],
Pre[201],Den[201],Vel[201],Tem[201],Ene[201],Ent[201],H[201],a[201],M[201],T,
dx=0.05,dt=0.00,CFL=0.98,PreL=10e5,PreR=1e5,DenL,DenR,Tem0=300.00,Ene0,H0,u0=0.00,a0,Vmax,amax,k=1.40;
R=287; /* 数值初始化*/
DenL=PreL/(R*Tem0);
DenR=PreR/(R*Tem0);
Ene0=R*Tem0/(k-1);
H0=Ene0+R*Tem0;
a0=(float)sqrt(k*R*Tem0);
for(i=1;i<=101;i++) /*通过循环将值赋到左边和右边的点 */
{
U1[i]=DenL;
U2[i]=DenL*u0;
U3[i]=DenL*Ene0;
F1[i]=DenL*u0;
F2[i]=DenL*u0*u0+PreL;
F3[i]=DenL*u0*H0;
}
for(i=101;i<=201;i++)
{
U1[i]=DenR;
U2[i]=DenR*u0;
U3[i]=DenR*Ene0;
F1[i]=DenR*u0;
F2[i]=DenR*u0*u0+PreR;
F3[i]=DenR*u0*H0;
}
T=0;
dt=CFL*dx/a0;
while(T<=0.005) /* 主循环,通过CFL数,对时间dt进行控制 */
{
for(i=1;i<=200;i++)
{
U1new[i]=(U1[i+1]+U1[i-1])/2-(dt/dx)*(F1[i+1]-F1[i-1])/2;
U2new[i]=(U2[i+1]+U2[i-1])/2-(dt/dx)*(F2[i+1]-F2[i-1])/2;/* 通过U1 U2 U3 F1 F2 F3 计算出新的U1new U2new U3new*/
U3new[i]=(U3[i+1]+U3[i-1])/2-(dt/dx)*(F3[i+1]-F3[i-1])/2;
}
for(i=1;i<=200;i++) /*通过循环将计算的U1new U2new U3new赋到原先的U1 U2 U3 F1 F2 F3*/
{
U1[i]=U1new[i];
U2[i]=U2new[i];
U3[i]=U3new[i];
Den[i]=U1new[i];
Vel[i]=U2new[i]/Den[i];
Pre[i]=U3new[i]*(k-1)/(Den[i]*R);
Ene[i]=U3new[i]/Den[i];
Tem[i]=(float)(Ene[i]-0.5*Vel[i]*Vel[i])*(k-1)/R;
Ent[i]=Ene[i]+R*Tem[i];
a[i]= (float)sqrt(k*R*Tem[i]);
M[i]=Vel[i]/a[i];
F1[i]=Den[i]*Vel[i];
F2[i]=Den[i]*Vel[i]*Vel[i]+Pre[i];
F3[i]=Den[i]*Vel[i]*H[i];
}
Vmax=compare(Vel,i);
amax=compare(a,i);
dt=CFL*dx/(Vmax+amax);
T=T+dt; /*根据计算的Vmax和amax计算新的dt*/
}
for(i=0;i<=201;i++)
{
printf("%f %f %f %f\n",Den[i],Pre[i],Vel[i],T); /*输出结果*/
}
}









