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); /*输出结果*/ } }