请教：怎么用MATLAB编程解MacCormack格式的声波连续方程
哪位高人知道怎么用MATLAB编程解MacCormack格式的声波连续方程。[ 本帖最后由 xinyuxf 于 2006-12-26 11:21 编辑 ] 有一个C语言的例子可以参考// CFD_1D_1st.cpp
// 一维激波管问题，采用二阶MacCormark格式，加入人工粘性以保证震荡稳定
//
// 最后修正2006.5.8
# include <math.h>
# include <iostream.h>
# include <iomanip.h>//引用setw()函数
//Max函数&Min函数，用于处理误差，保证健壮性
# define max_abs(a,b) ((fabs(a) > fabs(b)) ? (a) : (b))
# define min_abs(a,b) ((fabs(a) < fabs(b)) ? (a) : (b))
//常量定义
# define gama_ideal_gas 1.4//气体动力学中的常数，与气体类型有关，本值是相对理想气体而言
# define tolerance 1e-6//定义计算中的误差大小，小于之将被视为零，以保证程序的健壮性
# define L 0.2//经验的人工粘性系数，一般取0.6-3.0之间
//符号说明dens=density,velc=velocity,ener=energy,pres=pressure
void U_star_star (double *,double *,double *,double *,double *,double);
double ener_comp(double,double,double);
double pres_comp(double,double,double);
long U_norm(double,double,double,double,double *);
long F_fun_U(double *,double *);
int main ()
{
// 网格参数 x_range,t_range,Stmp_X,CFL or Stmp_T
const double Step_X=0.002,CFL=0.5,Step_T=Step_X*CFL;
const double x_range={-1,0,1};
const double t_range={0,0,0.3};
const long MESH_X_P1=(long) floor( (x_range-x_range)/Step_X )+1;
const long MESH_X_P2=(long) floor( (x_range-x_range)/Step_X )+1;
const long MESH_X =MESH_X_P1+MESH_X_P2;
const long MESH_T =(long) floor( (t_range-t_range)/Step_T )+2;
// 求解变量 dens,velc[][],ener[][],pres[][]
double *dens=new double[ MESH_T * MESH_X ];
double *velc=new double[ MESH_T * MESH_X ];
double *ener=new double[ MESH_T * MESH_X ];
double *pres=new double[ MESH_T * MESH_X ];
//初始条件
for(long i=0;i<MESH_X_P1;i++)
{
//第一层计算区域网格
dens =1;
velc =max_abs(0,tolerance);
pres =1;
ener =max_abs(ener_comp(pres,velc,dens ),tolerance);
//第一层冗余网格
dens =1;
velc =max_abs(0,tolerance);
pres =1;
ener =max_abs(ener_comp(pres,velc,dens ),tolerance);
}
for(i=0;i<MESH_X_P2;i++)
{
//第一层计算区域网格
dens =0.125;
velc =max_abs(0,tolerance);
pres =0.1;
ener =max_abs(ener_comp(pres,velc,dens ),tolerance);
//第一层冗余网格
dens =0.125;
velc =max_abs(0,tolerance);
pres =0.1;
ener =max_abs(ener_comp(pres,velc,dens ),tolerance);
}
//计算
double *U_2_star=new double ;//光滑前的U(n+1)(j)
double Q;//开关函数，源于人工粘性滤波算法
for(long n=1;n<MESH_T-1;n++)
{
if(n>1)
{ //开头与结尾的冗余网格，取前时刻本位置左（右）平均
dens =max_abs(0.5*(dens [(n-1)*MESH_X]+dens [(n-1)*MESH_X+1]),tolerance);
velc =max_abs(0.5*(velc [(n-1)*MESH_X]+velc [(n-1)*MESH_X+1]),tolerance);
pres =max_abs(0.5*(pres [(n-1)*MESH_X]+pres [(n-1)*MESH_X+1]),tolerance);
ener =max_abs(0.5*(ener [(n-1)*MESH_X]+ener [(n-1)*MESH_X+1]),tolerance);
dens [(n+1)*MESH_X-1]=max_abs(0.5*(dens +dens ),tolerance);
velc [(n+1)*MESH_X-1]=max_abs(0.5*(velc +velc ),tolerance);
pres [(n+1)*MESH_X-1]=max_abs(0.5*(pres +pres ),tolerance);
ener [(n+1)*MESH_X-1]=max_abs(0.5*(ener +ener ),tolerance);
}
//计算U_star_star向量
for(long j=1;j<MESH_X-1;j++)
{
U_star_star(dens+j+n*MESH_X,velc+j+n*MESH_X,ener+j+n*MESH_X,pres+j+n*MESH_X,U_2_star+3*j,CFL);
}
//U(n+1)(j)计算，由U_star_star的函数表达，抹平U_star_star的震荡
//其中U(n+1)(1)和U(n+1)(MESH_X-2)为两边边界，不可能有激波，不用抹平
//U(n+1)(1)
dens [(n+1)*MESH_X+1]=max_abs(U_2_star,tolerance);
velc [(n+1)*MESH_X+1]=max_abs(U_2_star/U_2_star,tolerance);
ener [(n+1)*MESH_X+1]=max_abs(U_2_star,tolerance);
pres [(n+1)*MESH_X+1]=max_abs(pres_comp(U_2_star,U_2_star,U_2_star),tolerance);
//U(n+1)(MESH_X-1)
dens [(n+2)*MESH_X-2]=max_abs(U_2_star,tolerance);
velc [(n+2)*MESH_X-2]=max_abs(U_2_star/U_2_star,tolerance);
ener [(n+2)*MESH_X-2]=max_abs(U_2_star,tolerance);
pres [(n+2)*MESH_X-2]=max_abs(pres_comp(U_2_star,U_2_star,U_2_star),tolerance);
//U(n+1)(j)
for(j=2;j<MESH_X-2;j++)
{
//开关函数Q
//MacCormack定义的开关函数，比Harten好用，但仍会震荡，实际中用l限制其范围，有效果，但不明显
//Q=fabs( ( pres + pres - 2*pres ) / max_abs( (pres + pres + 2*pres),1e-307 ) );
//Q=min_abs(fabs(Q),L);Q=((fabs(Q) < fabs(tolerance)) ? (0) : (Q));//将Q限制到0-L之间
Q=L;
//Harten定义的开关函数，不好用，会出现无穷
//Q=( fabs(pres-pres)-fabs(pres-pres) ) / ( fabs(pres-pres)+fabs(pres-pres) );
//cout<<Q<<" "<<endl;用于Q调错
double final_res;//U(n+1)(j)数据的临时保存
for(short t=0;t<3;t++)//计算光滑过的U(n+1)(j)
{final_res=U_2_star+0.5*Q*(U_2_star+U_2_star-2*U_2_star);}
//返回所求n+1时刻的dens,velc,ener,pres
dens [(n+1)*MESH_X+j]=max_abs(final_res,tolerance);
velc [(n+1)*MESH_X+j]=max_abs(final_res/final_res,tolerance);
ener [(n+1)*MESH_X+j]=max_abs(final_res,tolerance);
pres [(n+1)*MESH_X+j]=max_abs(pres_comp(final_res,final_res,final_res),tolerance);
}
}
//结果输出
for(n=1;n<MESH_T-1;n++)
{
cout<<n<<" 时刻"<<endl;
for(long j=1;j<MESH_X-1;j++)
{
if(fabs(dens)<=tolerance){dens=0;}
if(fabs(velc)<=tolerance){velc=0;}
if(fabs(pres)<=tolerance){pres=0;}
if(fabs(ener)<=tolerance){ener=0;}
cout<<setw(15)<<dens<<setw(15)<<velc<<setw(15)<<pres<<setw(15)<<ener<<endl;
}
cout<<""<<endl;
}
return 0;
}
void U_star_star (double *dens1,double *velc1,double *ener1,double *pres1,double *res,double CFL1)
{
double U_data;//分别存放U(n+1)(j)_star,U(n)(j),U(n+1)(j),U(n)(j+1),U(n)(j-1),U(n+1)(j-1)_star
double F_data;//分别存放F(n+1)(j)_star,F(n+1)(j-1)_star,F(n)(j),F(n)(j+1),F(n)(j-1)
double U_data_tmp;//临时存放从U_norm()拿到的数据
double F_data_tmp;//临时存放从F_fun_U()拿到的数据
//计算U(n)(j),数据临时存在U_data_tmp里
long err_flag=U_norm(*dens1,*velc1,*ener1,*pres1,U_data_tmp);
U_data=U_data_tmp;U_data=U_data_tmp;U_data=U_data_tmp;
//计算U(n)(j+1),数据临时存在U_data_tmp里
err_flag=U_norm(*(dens1+1),*(velc1+1),*(ener1+1),*(pres1+1),U_data_tmp);
U_data=U_data_tmp;U_data=U_data_tmp;U_data=U_data_tmp;
//计算U(n)(j-1),数据临时存在U_data_tmp里
err_flag=U_norm(*(dens1-1),*(velc1-1),*(ener1-1),*(pres1-1),U_data_tmp);
U_data=U_data_tmp;U_data=U_data_tmp;U_data=U_data_tmp;
//计算F(n)(j)
err_flag=F_fun_U(U_data+3,F_data_tmp);
F_data=F_data_tmp;F_data=F_data_tmp;F_data=F_data_tmp;
//计算F(n)(j+1)
err_flag=F_fun_U(U_data+9,F_data_tmp);
F_data=F_data_tmp;F_data=F_data_tmp;F_data=F_data_tmp;
//计算F(n)(j-1)
err_flag=F_fun_U(U_data+12,F_data_tmp);
F_data=F_data_tmp;F_data=F_data_tmp;F_data=F_data_tmp;
//计算U(n+1)(j)_star
U_data=U_data-CFL1*(F_data -F_data);
U_data=U_data-CFL1*(F_data-F_data);
U_data=U_data-CFL1*(F_data-F_data);
//计算U(n+1)(j-1)_star
U_data=U_data-CFL1*(F_data-F_data);
U_data=U_data-CFL1*(F_data-F_data);
U_data=U_data-CFL1*(F_data-F_data);
//计算F(n+1)(j)_star
err_flag=F_fun_U(U_data+0,F_data_tmp);
F_data=F_data_tmp;F_data=F_data_tmp;F_data=F_data_tmp;
//计算F(n+1)(j-1)_star
err_flag=F_fun_U(U_data+15,F_data_tmp);
F_data=F_data_tmp;F_data=F_data_tmp;F_data=F_data_tmp;
//计算U(n+1)(j)
U_data=0.5*(U_data+U_data)-0.5*CFL1*(F_data-F_data);
U_data=0.5*(U_data+U_data)-0.5*CFL1*(F_data-F_data);
U_data=0.5*(U_data+U_data)-0.5*CFL1*(F_data-F_data);
//返回U_star_star
res=U_data;//dens
res=U_data;//velc*dens
res=U_data;//ener
}
//用rou,u,p,gama表达e
double ener_comp(double pres_0,double velo_0,double dens_0)
{
return (pres_0/(gama_ideal_gas-1))+(0.5*dens_0*velo_0*velo_0);
}
//用rou,rou*u,e,gama表达p
double pres_comp(double ener_1,double velo_multi_dens_1,double dens_1)
{
return (gama_ideal_gas-1)*(ener_1-0.5*velo_multi_dens_1*velo_multi_dens_1/dens_1);
}
//用rou,u,p,gama,e表达U
long U_norm(double dens2,double velc2,double ener2,double pres2,double *U_data2)
{
U_data2=dens2;
U_data2=dens2*velc2;
U_data2=ener_comp(pres2,velc2,dens2);
return 0;
}
//F_fun_U是以U表达F的函数
long F_fun_U(double *U_data3,double *F_data3)
{
double P_tmp=(gama_ideal_gas-1)*(U_data3-0.5*U_data3*U_data3/U_data3);
F_data3=U_data3;
F_data3=U_data3*U_data3/U_data3+P_tmp;
F_data3=U_data3/U_data3*(U_data3+P_tmp);
return 0;
}// CFD_1D_1st.cpp
// 一维激波管问题，采用二阶MacCormark格式，加入人工粘性以保证震荡稳定
//
// 最后修正2006.5.8
# include <math.h>
# include <iostream.h>
# include <iomanip.h>//引用setw()函数
//Max函数&Min函数，用于处理误差，保证健壮性
# define max_abs(a,b) ((fabs(a) > fabs(b)) ? (a) : (b))
# define min_abs(a,b) ((fabs(a) < fabs(b)) ? (a) : (b))
//常量定义
# define gama_ideal_gas 1.4//气体动力学中的常数，与气体类型有关，本值是相对理想气体而言
# define tolerance 1e-6//定义计算中的误差大小，小于之将被视为零，以保证程序的健壮性
# define L 0.2//经验的人工粘性系数，一般取0.6-3.0之间
//符号说明dens=density,velc=velocity,ener=energy,pres=pressure
void U_star_star (double *,double *,double *,double *,double *,double);
double ener_comp(double,double,double);
double pres_comp(double,double,double);
long U_norm(double,double,double,double,double *);
long F_fun_U(double *,double *);
int main ()
{
// 网格参数 x_range,t_range,Stmp_X,CFL or Stmp_T
const double Step_X=0.002,CFL=0.5,Step_T=Step_X*CFL;
const double x_range={-1,0,1};
const double t_range={0,0,0.3};
const long MESH_X_P1=(long) floor( (x_range-x_range)/Step_X )+1;
const long MESH_X_P2=(long) floor( (x_range-x_range)/Step_X )+1;
const long MESH_X =MESH_X_P1+MESH_X_P2;
const long MESH_T =(long) floor( (t_range-t_range)/Step_T )+2;
// 求解变量 dens,velc[][],ener[][],pres[][]
double *dens=new double[ MESH_T * MESH_X ];
double *velc=new double[ MESH_T * MESH_X ];
double *ener=new double[ MESH_T * MESH_X ];
double *pres=new double[ MESH_T * MESH_X ];
//初始条件
for(long i=0;i<MESH_X_P1;i++)
{
//第一层计算区域网格
dens =1;
velc =max_abs(0,tolerance);
pres =1;
ener =max_abs(ener_comp(pres,velc,dens ),tolerance);
//第一层冗余网格
dens =1;
velc =max_abs(0,tolerance);
pres =1;
ener =max_abs(ener_comp(pres,velc,dens ),tolerance);
}
for(i=0;i<MESH_X_P2;i++)
{
//第一层计算区域网格
dens =0.125;
velc =max_abs(0,tolerance);
pres =0.1;
ener =max_abs(ener_comp(pres,velc,dens ),tolerance);
//第一层冗余网格
dens =0.125;
velc =max_abs(0,tolerance);
pres =0.1;
ener =max_abs(ener_comp(pres,velc,dens ),tolerance);
}
//计算
double *U_2_star=new double ;//光滑前的U(n+1)(j)
double Q;//开关函数，源于人工粘性滤波算法
for(long n=1;n<MESH_T-1;n++)
{
if(n>1)
{ //开头与结尾的冗余网格，取前时刻本位置左（右）平均
dens =max_abs(0.5*(dens [(n-1)*MESH_X]+dens [(n-1)*MESH_X+1]),tolerance);
velc =max_abs(0.5*(velc [(n-1)*MESH_X]+velc [(n-1)*MESH_X+1]),tolerance);
pres =max_abs(0.5*(pres [(n-1)*MESH_X]+pres [(n-1)*MESH_X+1]),tolerance);
ener =max_abs(0.5*(ener [(n-1)*MESH_X]+ener [(n-1)*MESH_X+1]),tolerance);
dens [(n+1)*MESH_X-1]=max_abs(0.5*(dens +dens ),tolerance);
velc [(n+1)*MESH_X-1]=max_abs(0.5*(velc +velc ),tolerance);
pres [(n+1)*MESH_X-1]=max_abs(0.5*(pres +pres ),tolerance);
ener [(n+1)*MESH_X-1]=max_abs(0.5*(ener +ener ),tolerance);
}
//计算U_star_star向量
for(long j=1;j<MESH_X-1;j++)
{
U_star_star(dens+j+n*MESH_X,velc+j+n*MESH_X,ener+j+n*MESH_X,pres+j+n*MESH_X,U_2_star+3*j,CFL);
}
//U(n+1)(j)计算，由U_star_star的函数表达，抹平U_star_star的震荡
//其中U(n+1)(1)和U(n+1)(MESH_X-2)为两边边界，不可能有激波，不用抹平
//U(n+1)(1)
dens [(n+1)*MESH_X+1]=max_abs(U_2_star,tolerance);
velc [(n+1)*MESH_X+1]=max_abs(U_2_star/U_2_star,tolerance);
ener [(n+1)*MESH_X+1]=max_abs(U_2_star,tolerance);
pres [(n+1)*MESH_X+1]=max_abs(pres_comp(U_2_star,U_2_star,U_2_star),tolerance);
//U(n+1)(MESH_X-1)
dens [(n+2)*MESH_X-2]=max_abs(U_2_star,tolerance);
velc [(n+2)*MESH_X-2]=max_abs(U_2_star/U_2_star,tolerance);
ener [(n+2)*MESH_X-2]=max_abs(U_2_star,tolerance);
pres [(n+2)*MESH_X-2]=max_abs(pres_comp(U_2_star,U_2_star,U_2_star),tolerance);
//U(n+1)(j)
for(j=2;j<MESH_X-2;j++)
{
//开关函数Q
//MacCormack定义的开关函数，比Harten好用，但仍会震荡，实际中用l限制其范围，有效果，但不明显
//Q=fabs( ( pres + pres - 2*pres ) / max_abs( (pres + pres + 2*pres),1e-307 ) );
//Q=min_abs(fabs(Q),L);Q=((fabs(Q) < fabs(tolerance)) ? (0) : (Q));//将Q限制到0-L之间
Q=L;
//Harten定义的开关函数，不好用，会出现无穷
//Q=( fabs(pres-pres)-fabs(pres-pres) ) / ( fabs(pres-pres)+fabs(pres-pres) );
//cout<<Q<<" "<<endl;用于Q调错
double final_res;//U(n+1)(j)数据的临时保存
for(short t=0;t<3;t++)//计算光滑过的U(n+1)(j)
{final_res=U_2_star+0.5*Q*(U_2_star+U_2_star-2*U_2_star);}
//返回所求n+1时刻的dens,velc,ener,pres
dens [(n+1)*MESH_X+j]=max_abs(final_res,tolerance);
velc [(n+1)*MESH_X+j]=max_abs(final_res/final_res,tolerance);
ener [(n+1)*MESH_X+j]=max_abs(final_res,tolerance);
pres [(n+1)*MESH_X+j]=max_abs(pres_comp(final_res,final_res,final_res),tolerance);
}
}
//结果输出
for(n=1;n<MESH_T-1;n++)
{
cout<<n<<" 时刻"<<endl;
for(long j=1;j<MESH_X-1;j++)
{
if(fabs(dens)<=tolerance){dens=0;}
if(fabs(velc)<=tolerance){velc=0;}
if(fabs(pres)<=tolerance){pres=0;}
if(fabs(ener)<=tolerance){ener=0;}
cout<<setw(15)<<dens<<setw(15)<<velc<<setw(15)<<pres<<setw(15)<<ener<<endl;
}
cout<<""<<endl;
}
return 0;
}
void U_star_star (double *dens1,double *velc1,double *ener1,double *pres1,double *res,double CFL1)
{
double U_data;//分别存放U(n+1)(j)_star,U(n)(j),U(n+1)(j),U(n)(j+1),U(n)(j-1),U(n+1)(j-1)_star
double F_data;//分别存放F(n+1)(j)_star,F(n+1)(j-1)_star,F(n)(j),F(n)(j+1),F(n)(j-1)
double U_data_tmp;//临时存放从U_norm()拿到的数据
double F_data_tmp;//临时存放从F_fun_U()拿到的数据
//计算U(n)(j),数据临时存在U_data_tmp里
long err_flag=U_norm(*dens1,*velc1,*ener1,*pres1,U_data_tmp);
U_data=U_data_tmp;U_data=U_data_tmp;U_data=U_data_tmp;
//计算U(n)(j+1),数据临时存在U_data_tmp里
err_flag=U_norm(*(dens1+1),*(velc1+1),*(ener1+1),*(pres1+1),U_data_tmp);
U_data=U_data_tmp;U_data=U_data_tmp;U_data=U_data_tmp;
//计算U(n)(j-1),数据临时存在U_data_tmp里
err_flag=U_norm(*(dens1-1),*(velc1-1),*(ener1-1),*(pres1-1),U_data_tmp);
U_data=U_data_tmp;U_data=U_data_tmp;U_data=U_data_tmp;
//计算F(n)(j)
err_flag=F_fun_U(U_data+3,F_data_tmp);
F_data=F_data_tmp;F_data=F_data_tmp;F_data=F_data_tmp;
//计算F(n)(j+1)
err_flag=F_fun_U(U_data+9,F_data_tmp);
F_data=F_data_tmp;F_data=F_data_tmp;F_data=F_data_tmp;
//计算F(n)(j-1)
err_flag=F_fun_U(U_data+12,F_data_tmp);
F_data=F_data_tmp;F_data=F_data_tmp;F_data=F_data_tmp;
//计算U(n+1)(j)_star
U_data=U_data-CFL1*(F_data -F_data);
U_data=U_data-CFL1*(F_data-F_data);
U_data=U_data-CFL1*(F_data-F_data);
//计算U(n+1)(j-1)_star
U_data=U_data-CFL1*(F_data-F_data);
U_data=U_data-CFL1*(F_data-F_data);
U_data=U_data-CFL1*(F_data-F_data);
//计算F(n+1)(j)_star
err_flag=F_fun_U(U_data+0,F_data_tmp);
F_data=F_data_tmp;F_data=F_data_tmp;F_data=F_data_tmp;
//计算F(n+1)(j-1)_star
err_flag=F_fun_U(U_data+15,F_data_tmp);
F_data=F_data_tmp;F_data=F_data_tmp;F_data=F_data_tmp;
//计算U(n+1)(j)
U_data=0.5*(U_data+U_data)-0.5*CFL1*(F_data-F_data);
U_data=0.5*(U_data+U_data)-0.5*CFL1*(F_data-F_data);
U_data=0.5*(U_data+U_data)-0.5*CFL1*(F_data-F_data);
//返回U_star_star
res=U_data;//dens
res=U_data;//velc*dens
res=U_data;//ener
}
//用rou,u,p,gama表达e
double ener_comp(double pres_0,double velo_0,double dens_0)
{
return (pres_0/(gama_ideal_gas-1))+(0.5*dens_0*velo_0*velo_0);
}
//用rou,rou*u,e,gama表达p
double pres_comp(double ener_1,double velo_multi_dens_1,double dens_1)
{
return (gama_ideal_gas-1)*(ener_1-0.5*velo_multi_dens_1*velo_multi_dens_1/dens_1);
}
//用rou,u,p,gama,e表达U
long U_norm(double dens2,double velc2,double ener2,double pres2,double *U_data2)
{
U_data2=dens2;
U_data2=dens2*velc2;
U_data2=ener_comp(pres2,velc2,dens2);
return 0;
}
//F_fun_U是以U表达F的函数
long F_fun_U(double *U_data3,double *F_data3)
{
double P_tmp=(gama_ideal_gas-1)*(U_data3-0.5*U_data3*U_data3/U_data3);
F_data3=U_data3;
F_data3=U_data3*U_data3/U_data3+P_tmp;
F_data3=U_data3/U_data3*(U_data3+P_tmp);
return 0;
}
来自：http://swell2005.yculblog.com/post.1267582.html
页:
[1]