声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2638|回复: 1

[近似分析] 请教:怎么用MATLAB编程解MacCormack格式的声波连续方程

[复制链接]
发表于 2006-10-19 21:57 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
哪位高人知道怎么用MATLAB编程解MacCormack格式的声波连续方程。

[ 本帖最后由 xinyuxf 于 2006-12-26 11:21 编辑 ]
回复
分享到:

使用道具 举报

发表于 2007-5-2 05:15 | 显示全部楼层
有一个C语言的例子可以参考
  1. // CFD_1D_1st.cpp
  2. // 一维激波管问题,采用二阶MacCormark格式,加入人工粘性以保证震荡稳定
  3. //
  4. // 最后修正2006.5.8

  5. # include <math.h>
  6. # include <iostream.h>
  7. # include <iomanip.h>//引用setw()函数
  8. //Max函数&Min函数,用于处理误差,保证健壮性
  9. # define max_abs(a,b) ((fabs(a) > fabs(b)) ? (a) : (b))
  10. # define min_abs(a,b) ((fabs(a) < fabs(b)) ? (a) : (b))
  11. //常量定义
  12. # define gama_ideal_gas 1.4//气体动力学中的常数,与气体类型有关,本值是相对理想气体而言
  13. # define tolerance 1e-6//定义计算中的误差大小,小于之将被视为零,以保证程序的健壮性
  14. # define L 0.2//经验的人工粘性系数,一般取0.6-3.0之间
  15. //符号说明dens=density,velc=velocity,ener=energy,pres=pressure

  16. void U_star_star (double *,double *,double *,double *,double *,double);
  17. double ener_comp(double,double,double);
  18. double pres_comp(double,double,double);
  19. long U_norm(double,double,double,double,double *);
  20. long F_fun_U(double *,double *);

  21. int main ()
  22. {
  23.         // 网格参数 x_range[3],t_range[2],Stmp_X,CFL or Stmp_T
  24.         const double Step_X=0.002,CFL=0.5,Step_T=Step_X*CFL;
  25.         const double x_range[3]={-1,0,1};
  26.         const double t_range[3]={0,0,0.3};
  27.         
  28.         const long MESH_X_P1=(long) floor( (x_range[1]-x_range[0])/Step_X )+1;
  29.         const long MESH_X_P2=(long) floor( (x_range[2]-x_range[1])/Step_X )+1;
  30.         const long MESH_X       =MESH_X_P1+MESH_X_P2;
  31.         const long MESH_T       =(long) floor( (t_range[2]-t_range[0])/Step_T )+2;
  32.         
  33.         // 求解变量 dens[MESH_T][MESH_X],velc[][],ener[][],pres[][]
  34.         double *dens=new double[ MESH_T * MESH_X ];
  35.         double *velc=new double[ MESH_T * MESH_X ];
  36.         double *ener=new double[ MESH_T * MESH_X ];
  37.         double *pres=new double[ MESH_T * MESH_X ];

  38.         //初始条件
  39.         for(long i=0;i<MESH_X_P1;i++)
  40.         {
  41.                 //第一层计算区域网格
  42.                 dens [i+MESH_X]=1;
  43.                 velc [i+MESH_X]=max_abs(0,tolerance);
  44.                 pres [i+MESH_X]=1;
  45.                 ener [i+MESH_X]=max_abs(ener_comp(pres[i+MESH_X],velc[i+MESH_X],dens [i+MESH_X]),tolerance);
  46.                 //第一层冗余网格
  47.                 dens [i]=1;
  48.                 velc [i]=max_abs(0,tolerance);
  49.                 pres [i]=1;
  50.                 ener [i]=max_abs(ener_comp(pres[i],velc[i],dens [i]),tolerance);
  51.         }
  52.         for(i=0;i<MESH_X_P2;i++)
  53.         {
  54.                 //第一层计算区域网格
  55.                 dens [i+MESH_X_P1+MESH_X]=0.125;
  56.                 velc [i+MESH_X_P1+MESH_X]=max_abs(0,tolerance);
  57.                 pres [i+MESH_X_P1+MESH_X]=0.1;
  58.                 ener [i+MESH_X_P1+MESH_X]=max_abs(ener_comp(pres[i+MESH_X_P1+MESH_X],velc[i+MESH_X_P1+MESH_X],dens [i+MESH_X_P1+MESH_X]),tolerance);
  59.                 //第一层冗余网格
  60.                 dens [i+MESH_X_P1]=0.125;
  61.                 velc [i+MESH_X_P1]=max_abs(0,tolerance);
  62.                 pres [i+MESH_X_P1]=0.1;
  63.                 ener [i+MESH_X_P1]=max_abs(ener_comp(pres[i+MESH_X_P1],velc[i+MESH_X_P1],dens [i+MESH_X_P1]),tolerance);
  64.         }

  65.         //计算
  66.         
  67.         double *U_2_star=new double [3*MESH_X];//光滑前的U(n+1)(j)
  68.         double Q;//开关函数,源于人工粘性滤波算法
  69.         for(long n=1;n<MESH_T-1;n++)
  70.         {      
  71.                 if(n>1)
  72.                 {       //开头与结尾的冗余网格,取前时刻本位置左(右)平均
  73.                         dens [n*MESH_X]=max_abs(0.5*(dens [(n-1)*MESH_X]+dens [(n-1)*MESH_X+1]),tolerance);
  74.                         velc [n*MESH_X]=max_abs(0.5*(velc [(n-1)*MESH_X]+velc [(n-1)*MESH_X+1]),tolerance);
  75.                         pres [n*MESH_X]=max_abs(0.5*(pres [(n-1)*MESH_X]+pres [(n-1)*MESH_X+1]),tolerance);
  76.                         ener [n*MESH_X]=max_abs(0.5*(ener [(n-1)*MESH_X]+ener [(n-1)*MESH_X+1]),tolerance);

  77.                         dens [(n+1)*MESH_X-1]=max_abs(0.5*(dens [n*MESH_X-1]+dens [n*MESH_X-2]),tolerance);
  78.                         velc [(n+1)*MESH_X-1]=max_abs(0.5*(velc [n*MESH_X-1]+velc [n*MESH_X-2]),tolerance);
  79.                         pres [(n+1)*MESH_X-1]=max_abs(0.5*(pres [n*MESH_X-1]+pres [n*MESH_X-2]),tolerance);
  80.                         ener [(n+1)*MESH_X-1]=max_abs(0.5*(ener [n*MESH_X-1]+ener [n*MESH_X-2]),tolerance);
  81.                 }
  82.                
  83.                 //计算U_star_star向量
  84.                 for(long j=1;j<MESH_X-1;j++)
  85.                 {
  86.                         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);
  87.                 }
  88.                
  89.                 //U(n+1)(j)计算,由U_star_star的函数表达,抹平U_star_star的震荡
  90.                 //其中U(n+1)(1)和U(n+1)(MESH_X-2)为两边边界,不可能有激波,不用抹平
  91.                 //U(n+1)(1)
  92.                         dens [(n+1)*MESH_X+1]=max_abs(U_2_star[3],tolerance);
  93.                         velc [(n+1)*MESH_X+1]=max_abs(U_2_star[4]/U_2_star[3],tolerance);
  94.                         ener [(n+1)*MESH_X+1]=max_abs(U_2_star[5],tolerance);
  95.                         pres [(n+1)*MESH_X+1]=max_abs(pres_comp(U_2_star[5],U_2_star[4],U_2_star[3]),tolerance);
  96.                 //U(n+1)(MESH_X-1)
  97.                         dens [(n+2)*MESH_X-2]=max_abs(U_2_star[3*(MESH_X-2)],tolerance);
  98.                         velc [(n+2)*MESH_X-2]=max_abs(U_2_star[3*(MESH_X-2)+1]/U_2_star[3*(MESH_X-2)],tolerance);
  99.                         ener [(n+2)*MESH_X-2]=max_abs(U_2_star[3*(MESH_X-2)+2],tolerance);
  100.                         pres [(n+2)*MESH_X-2]=max_abs(pres_comp(U_2_star[3*(MESH_X-2)+2],U_2_star[3*(MESH_X-2)+1],U_2_star[3*(MESH_X-2)]),tolerance);
  101.                 //U(n+1)(j)
  102.                 for(j=2;j<MESH_X-2;j++)
  103.                 {
  104.                         //开关函数Q
  105.                         //MacCormack定义的开关函数,比Harten好用,但仍会震荡,实际中用l限制其范围,有效果,但不明显
  106.                         //Q=fabs( ( pres[n*MESH_X+j+1] + pres[n*MESH_X+j-1] - 2*pres[n*MESH_X+j] ) / max_abs( (pres[n*MESH_X+j+1] + pres[n*MESH_X+j-1] + 2*pres[n*MESH_X+j]),1e-307 ) );
  107.                         //Q=min_abs(fabs(Q),L);Q=((fabs(Q) < fabs(tolerance)) ? (0) : (Q));//将Q限制到0-L之间
  108.                         Q=L;
  109.                         //Harten定义的开关函数,不好用,会出现无穷
  110.                         //Q=( fabs(pres[n*MESH_X+j+1]-pres[n*MESH_X+j])-fabs(pres[n*MESH_X+j]-pres[n*MESH_X+j-1]) ) / ( fabs(pres[n*MESH_X+j+1]-pres[n*MESH_X+j])+fabs(pres[n*MESH_X+j]-pres[n*MESH_X+j-1]) );
  111.                         //cout<<Q<<"            "<<endl;用于Q调错
  112.                         double final_res[3];//U(n+1)(j)数据的临时保存
  113.                         for(short t=0;t<3;t++)//计算光滑过的U(n+1)(j)
  114.                         {final_res[t]=U_2_star[3*j+t]+0.5*Q*(U_2_star[3*(j+1)+t]+U_2_star[3*(j-1)+t]-2*U_2_star[3*j+t]);}
  115.                         //返回所求n+1时刻的dens,velc,ener,pres
  116.                         dens [(n+1)*MESH_X+j]=max_abs(final_res[0],tolerance);
  117.                         velc [(n+1)*MESH_X+j]=max_abs(final_res[1]/final_res[0],tolerance);
  118.                         ener [(n+1)*MESH_X+j]=max_abs(final_res[2],tolerance);
  119.                         pres [(n+1)*MESH_X+j]=max_abs(pres_comp(final_res[2],final_res[1],final_res[0]),tolerance);                     
  120.                 }
  121.         }
  122.         
  123.         //结果输出
  124.         for(n=1;n<MESH_T-1;n++)
  125.         {
  126.                 cout<<n<<"   时刻"<<endl;
  127.                 for(long j=1;j<MESH_X-1;j++)
  128.                 {
  129.                         if(fabs(dens[n*MESH_X+j])<=tolerance){dens[n*MESH_X+j]=0;}
  130.                         if(fabs(velc[n*MESH_X+j])<=tolerance){velc[n*MESH_X+j]=0;}
  131.                         if(fabs(pres[n*MESH_X+j])<=tolerance){pres[n*MESH_X+j]=0;}
  132.                         if(fabs(ener[n*MESH_X+j])<=tolerance){ener[n*MESH_X+j]=0;}

  133.                   cout<<setw(15)<<dens[n*MESH_X+j]<<setw(15)<<velc[n*MESH_X+j]<<setw(15)<<pres[n*MESH_X+j]<<setw(15)<<ener[n*MESH_X+j]<<endl;
  134.                 }
  135.                 cout<<""<<endl;
  136.         }
  137.         return 0;
  138. }

  139. void U_star_star (double *dens1,double *velc1,double *ener1,double *pres1,double *res,double CFL1)
  140. {
  141.         double U_data[18];//分别存放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
  142.         double F_data[15];//分别存放F(n+1)(j)_star,F(n+1)(j-1)_star,F(n)(j),F(n)(j+1),F(n)(j-1)
  143.         double U_data_tmp[3];//临时存放从U_norm()拿到的数据
  144.         double F_data_tmp[3];//临时存放从F_fun_U()拿到的数据
  145.         
  146.         //计算U(n)(j),数据临时存在U_data_tmp里
  147.         long err_flag=U_norm(*dens1,*velc1,*ener1,*pres1,U_data_tmp);
  148.         U_data[3]=U_data_tmp[0];U_data[4]=U_data_tmp[1];U_data[5]=U_data_tmp[2];

  149.         //计算U(n)(j+1),数据临时存在U_data_tmp里
  150.         err_flag=U_norm(*(dens1+1),*(velc1+1),*(ener1+1),*(pres1+1),U_data_tmp);
  151.         U_data[9]=U_data_tmp[0];U_data[10]=U_data_tmp[1];U_data[11]=U_data_tmp[2];

  152.         //计算U(n)(j-1),数据临时存在U_data_tmp里
  153.         err_flag=U_norm(*(dens1-1),*(velc1-1),*(ener1-1),*(pres1-1),U_data_tmp);
  154.         U_data[12]=U_data_tmp[0];U_data[13]=U_data_tmp[1];U_data[14]=U_data_tmp[2];

  155.         //计算F(n)(j)
  156.         err_flag=F_fun_U(U_data+3,F_data_tmp);
  157.         F_data[6]=F_data_tmp[0];F_data[7]=F_data_tmp[1];F_data[8]=F_data_tmp[2];

  158.         //计算F(n)(j+1)
  159.         err_flag=F_fun_U(U_data+9,F_data_tmp);
  160.         F_data[9]=F_data_tmp[0];F_data[10]=F_data_tmp[1];F_data[11]=F_data_tmp[2];
  161.         
  162.         //计算F(n)(j-1)
  163.         err_flag=F_fun_U(U_data+12,F_data_tmp);
  164.         F_data[12]=F_data_tmp[0];F_data[13]=F_data_tmp[1];F_data[14]=F_data_tmp[2];
  165.                
  166.         //计算U(n+1)(j)_star
  167.         U_data[0]=U_data[3]-CFL1*(F_data[9] -F_data[6]);
  168.         U_data[1]=U_data[4]-CFL1*(F_data[10]-F_data[7]);
  169.         U_data[2]=U_data[5]-CFL1*(F_data[11]-F_data[8]);

  170.         //计算U(n+1)(j-1)_star
  171.         U_data[15]=U_data[12]-CFL1*(F_data[6]-F_data[12]);
  172.         U_data[16]=U_data[13]-CFL1*(F_data[7]-F_data[13]);
  173.         U_data[17]=U_data[14]-CFL1*(F_data[8]-F_data[14]);
  174.         
  175.         //计算F(n+1)(j)_star
  176.         err_flag=F_fun_U(U_data+0,F_data_tmp);
  177.         F_data[0]=F_data_tmp[0];F_data[1]=F_data_tmp[1];F_data[2]=F_data_tmp[2];
  178.         
  179.         //计算F(n+1)(j-1)_star
  180.         err_flag=F_fun_U(U_data+15,F_data_tmp);
  181.         F_data[3]=F_data_tmp[0];F_data[4]=F_data_tmp[1];F_data[5]=F_data_tmp[2];
  182.         
  183.         //计算U(n+1)(j)
  184.         U_data[6]=0.5*(U_data[3]+U_data[0])-0.5*CFL1*(F_data[0]-F_data[3]);
  185.         U_data[7]=0.5*(U_data[4]+U_data[1])-0.5*CFL1*(F_data[1]-F_data[4]);
  186.         U_data[8]=0.5*(U_data[5]+U_data[2])-0.5*CFL1*(F_data[2]-F_data[5]);

  187.         //返回U_star_star
  188.         res[0]=U_data[6];//dens
  189.         res[1]=U_data[7];//velc*dens
  190.         res[2]=U_data[8];//ener
  191. }

  192. //用rou,u,p,gama表达e
  193. double ener_comp(double pres_0,double velo_0,double dens_0)
  194. {
  195.         return (pres_0/(gama_ideal_gas-1))+(0.5*dens_0*velo_0*velo_0);
  196. }

  197. //用rou,rou*u,e,gama表达p
  198. double pres_comp(double ener_1,double velo_multi_dens_1,double dens_1)
  199. {
  200.         return (gama_ideal_gas-1)*(ener_1-0.5*velo_multi_dens_1*velo_multi_dens_1/dens_1);
  201. }

  202. //用rou,u,p,gama,e表达U
  203. long U_norm(double dens2,double velc2,double ener2,double pres2,double *U_data2)
  204. {
  205.         U_data2[0]=dens2;
  206.         U_data2[1]=dens2*velc2;
  207.         U_data2[2]=ener_comp(pres2,velc2,dens2);
  208.         return 0;
  209. }

  210. //F_fun_U是以U表达F的函数
  211. long F_fun_U(double *U_data3,double *F_data3)
  212. {
  213.         double P_tmp=(gama_ideal_gas-1)*(U_data3[2]-0.5*U_data3[1]*U_data3[1]/U_data3[0]);
  214.         F_data3[0]=U_data3[1];
  215.         F_data3[1]=U_data3[1]*U_data3[1]/U_data3[0]+P_tmp;
  216.         F_data3[2]=U_data3[1]/U_data3[0]*(U_data3[2]+P_tmp);
  217.         return 0;
  218. }// CFD_1D_1st.cpp
  219. // 一维激波管问题,采用二阶MacCormark格式,加入人工粘性以保证震荡稳定
  220. //
  221. // 最后修正2006.5.8

  222. # include <math.h>
  223. # include <iostream.h>
  224. # include <iomanip.h>//引用setw()函数
  225. //Max函数&Min函数,用于处理误差,保证健壮性
  226. # define max_abs(a,b) ((fabs(a) > fabs(b)) ? (a) : (b))
  227. # define min_abs(a,b) ((fabs(a) < fabs(b)) ? (a) : (b))
  228. //常量定义
  229. # define gama_ideal_gas 1.4//气体动力学中的常数,与气体类型有关,本值是相对理想气体而言
  230. # define tolerance 1e-6//定义计算中的误差大小,小于之将被视为零,以保证程序的健壮性
  231. # define L 0.2//经验的人工粘性系数,一般取0.6-3.0之间
  232. //符号说明dens=density,velc=velocity,ener=energy,pres=pressure

  233. void U_star_star (double *,double *,double *,double *,double *,double);
  234. double ener_comp(double,double,double);
  235. double pres_comp(double,double,double);
  236. long U_norm(double,double,double,double,double *);
  237. long F_fun_U(double *,double *);

  238. int main ()
  239. {
  240.         // 网格参数 x_range[3],t_range[2],Stmp_X,CFL or Stmp_T
  241.         const double Step_X=0.002,CFL=0.5,Step_T=Step_X*CFL;
  242.         const double x_range[3]={-1,0,1};
  243.         const double t_range[3]={0,0,0.3};
  244.         
  245.         const long MESH_X_P1=(long) floor( (x_range[1]-x_range[0])/Step_X )+1;
  246.         const long MESH_X_P2=(long) floor( (x_range[2]-x_range[1])/Step_X )+1;
  247.         const long MESH_X       =MESH_X_P1+MESH_X_P2;
  248.         const long MESH_T       =(long) floor( (t_range[2]-t_range[0])/Step_T )+2;
  249.         
  250.         // 求解变量 dens[MESH_T][MESH_X],velc[][],ener[][],pres[][]
  251.         double *dens=new double[ MESH_T * MESH_X ];
  252.         double *velc=new double[ MESH_T * MESH_X ];
  253.         double *ener=new double[ MESH_T * MESH_X ];
  254.         double *pres=new double[ MESH_T * MESH_X ];

  255.         //初始条件
  256.         for(long i=0;i<MESH_X_P1;i++)
  257.         {
  258.                 //第一层计算区域网格
  259.                 dens [i+MESH_X]=1;
  260.                 velc [i+MESH_X]=max_abs(0,tolerance);
  261.                 pres [i+MESH_X]=1;
  262.                 ener [i+MESH_X]=max_abs(ener_comp(pres[i+MESH_X],velc[i+MESH_X],dens [i+MESH_X]),tolerance);
  263.                 //第一层冗余网格
  264.                 dens [i]=1;
  265.                 velc [i]=max_abs(0,tolerance);
  266.                 pres [i]=1;
  267.                 ener [i]=max_abs(ener_comp(pres[i],velc[i],dens [i]),tolerance);
  268.         }
  269.         for(i=0;i<MESH_X_P2;i++)
  270.         {
  271.                 //第一层计算区域网格
  272.                 dens [i+MESH_X_P1+MESH_X]=0.125;
  273.                 velc [i+MESH_X_P1+MESH_X]=max_abs(0,tolerance);
  274.                 pres [i+MESH_X_P1+MESH_X]=0.1;
  275.                 ener [i+MESH_X_P1+MESH_X]=max_abs(ener_comp(pres[i+MESH_X_P1+MESH_X],velc[i+MESH_X_P1+MESH_X],dens [i+MESH_X_P1+MESH_X]),tolerance);
  276.                 //第一层冗余网格
  277.                 dens [i+MESH_X_P1]=0.125;
  278.                 velc [i+MESH_X_P1]=max_abs(0,tolerance);
  279.                 pres [i+MESH_X_P1]=0.1;
  280.                 ener [i+MESH_X_P1]=max_abs(ener_comp(pres[i+MESH_X_P1],velc[i+MESH_X_P1],dens [i+MESH_X_P1]),tolerance);
  281.         }

  282.         //计算
  283.         
  284.         double *U_2_star=new double [3*MESH_X];//光滑前的U(n+1)(j)
  285.         double Q;//开关函数,源于人工粘性滤波算法
  286.         for(long n=1;n<MESH_T-1;n++)
  287.         {      
  288.                 if(n>1)
  289.                 {       //开头与结尾的冗余网格,取前时刻本位置左(右)平均
  290.                         dens [n*MESH_X]=max_abs(0.5*(dens [(n-1)*MESH_X]+dens [(n-1)*MESH_X+1]),tolerance);
  291.                         velc [n*MESH_X]=max_abs(0.5*(velc [(n-1)*MESH_X]+velc [(n-1)*MESH_X+1]),tolerance);
  292.                         pres [n*MESH_X]=max_abs(0.5*(pres [(n-1)*MESH_X]+pres [(n-1)*MESH_X+1]),tolerance);
  293.                         ener [n*MESH_X]=max_abs(0.5*(ener [(n-1)*MESH_X]+ener [(n-1)*MESH_X+1]),tolerance);

  294.                         dens [(n+1)*MESH_X-1]=max_abs(0.5*(dens [n*MESH_X-1]+dens [n*MESH_X-2]),tolerance);
  295.                         velc [(n+1)*MESH_X-1]=max_abs(0.5*(velc [n*MESH_X-1]+velc [n*MESH_X-2]),tolerance);
  296.                         pres [(n+1)*MESH_X-1]=max_abs(0.5*(pres [n*MESH_X-1]+pres [n*MESH_X-2]),tolerance);
  297.                         ener [(n+1)*MESH_X-1]=max_abs(0.5*(ener [n*MESH_X-1]+ener [n*MESH_X-2]),tolerance);
  298.                 }
  299.                
  300.                 //计算U_star_star向量
  301.                 for(long j=1;j<MESH_X-1;j++)
  302.                 {
  303.                         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);
  304.                 }
  305.                
  306.                 //U(n+1)(j)计算,由U_star_star的函数表达,抹平U_star_star的震荡
  307.                 //其中U(n+1)(1)和U(n+1)(MESH_X-2)为两边边界,不可能有激波,不用抹平
  308.                 //U(n+1)(1)
  309.                         dens [(n+1)*MESH_X+1]=max_abs(U_2_star[3],tolerance);
  310.                         velc [(n+1)*MESH_X+1]=max_abs(U_2_star[4]/U_2_star[3],tolerance);
  311.                         ener [(n+1)*MESH_X+1]=max_abs(U_2_star[5],tolerance);
  312.                         pres [(n+1)*MESH_X+1]=max_abs(pres_comp(U_2_star[5],U_2_star[4],U_2_star[3]),tolerance);
  313.                 //U(n+1)(MESH_X-1)
  314.                         dens [(n+2)*MESH_X-2]=max_abs(U_2_star[3*(MESH_X-2)],tolerance);
  315.                         velc [(n+2)*MESH_X-2]=max_abs(U_2_star[3*(MESH_X-2)+1]/U_2_star[3*(MESH_X-2)],tolerance);
  316.                         ener [(n+2)*MESH_X-2]=max_abs(U_2_star[3*(MESH_X-2)+2],tolerance);
  317.                         pres [(n+2)*MESH_X-2]=max_abs(pres_comp(U_2_star[3*(MESH_X-2)+2],U_2_star[3*(MESH_X-2)+1],U_2_star[3*(MESH_X-2)]),tolerance);
  318.                 //U(n+1)(j)
  319.                 for(j=2;j<MESH_X-2;j++)
  320.                 {
  321.                         //开关函数Q
  322.                         //MacCormack定义的开关函数,比Harten好用,但仍会震荡,实际中用l限制其范围,有效果,但不明显
  323.                         //Q=fabs( ( pres[n*MESH_X+j+1] + pres[n*MESH_X+j-1] - 2*pres[n*MESH_X+j] ) / max_abs( (pres[n*MESH_X+j+1] + pres[n*MESH_X+j-1] + 2*pres[n*MESH_X+j]),1e-307 ) );
  324.                         //Q=min_abs(fabs(Q),L);Q=((fabs(Q) < fabs(tolerance)) ? (0) : (Q));//将Q限制到0-L之间
  325.                         Q=L;
  326.                         //Harten定义的开关函数,不好用,会出现无穷
  327.                         //Q=( fabs(pres[n*MESH_X+j+1]-pres[n*MESH_X+j])-fabs(pres[n*MESH_X+j]-pres[n*MESH_X+j-1]) ) / ( fabs(pres[n*MESH_X+j+1]-pres[n*MESH_X+j])+fabs(pres[n*MESH_X+j]-pres[n*MESH_X+j-1]) );
  328.                         //cout<<Q<<"            "<<endl;用于Q调错
  329.                         double final_res[3];//U(n+1)(j)数据的临时保存
  330.                         for(short t=0;t<3;t++)//计算光滑过的U(n+1)(j)
  331.                         {final_res[t]=U_2_star[3*j+t]+0.5*Q*(U_2_star[3*(j+1)+t]+U_2_star[3*(j-1)+t]-2*U_2_star[3*j+t]);}
  332.                         //返回所求n+1时刻的dens,velc,ener,pres
  333.                         dens [(n+1)*MESH_X+j]=max_abs(final_res[0],tolerance);
  334.                         velc [(n+1)*MESH_X+j]=max_abs(final_res[1]/final_res[0],tolerance);
  335.                         ener [(n+1)*MESH_X+j]=max_abs(final_res[2],tolerance);
  336.                         pres [(n+1)*MESH_X+j]=max_abs(pres_comp(final_res[2],final_res[1],final_res[0]),tolerance);                     
  337.                 }
  338.         }
  339.         
  340.         //结果输出
  341.         for(n=1;n<MESH_T-1;n++)
  342.         {
  343.                 cout<<n<<"   时刻"<<endl;
  344.                 for(long j=1;j<MESH_X-1;j++)
  345.                 {
  346.                         if(fabs(dens[n*MESH_X+j])<=tolerance){dens[n*MESH_X+j]=0;}
  347.                         if(fabs(velc[n*MESH_X+j])<=tolerance){velc[n*MESH_X+j]=0;}
  348.                         if(fabs(pres[n*MESH_X+j])<=tolerance){pres[n*MESH_X+j]=0;}
  349.                         if(fabs(ener[n*MESH_X+j])<=tolerance){ener[n*MESH_X+j]=0;}

  350.                   cout<<setw(15)<<dens[n*MESH_X+j]<<setw(15)<<velc[n*MESH_X+j]<<setw(15)<<pres[n*MESH_X+j]<<setw(15)<<ener[n*MESH_X+j]<<endl;
  351.                 }
  352.                 cout<<""<<endl;
  353.         }
  354.         return 0;
  355. }

  356. void U_star_star (double *dens1,double *velc1,double *ener1,double *pres1,double *res,double CFL1)
  357. {
  358.         double U_data[18];//分别存放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
  359.         double F_data[15];//分别存放F(n+1)(j)_star,F(n+1)(j-1)_star,F(n)(j),F(n)(j+1),F(n)(j-1)
  360.         double U_data_tmp[3];//临时存放从U_norm()拿到的数据
  361.         double F_data_tmp[3];//临时存放从F_fun_U()拿到的数据
  362.         
  363.         //计算U(n)(j),数据临时存在U_data_tmp里
  364.         long err_flag=U_norm(*dens1,*velc1,*ener1,*pres1,U_data_tmp);
  365.         U_data[3]=U_data_tmp[0];U_data[4]=U_data_tmp[1];U_data[5]=U_data_tmp[2];

  366.         //计算U(n)(j+1),数据临时存在U_data_tmp里
  367.         err_flag=U_norm(*(dens1+1),*(velc1+1),*(ener1+1),*(pres1+1),U_data_tmp);
  368.         U_data[9]=U_data_tmp[0];U_data[10]=U_data_tmp[1];U_data[11]=U_data_tmp[2];

  369.         //计算U(n)(j-1),数据临时存在U_data_tmp里
  370.         err_flag=U_norm(*(dens1-1),*(velc1-1),*(ener1-1),*(pres1-1),U_data_tmp);
  371.         U_data[12]=U_data_tmp[0];U_data[13]=U_data_tmp[1];U_data[14]=U_data_tmp[2];

  372.         //计算F(n)(j)
  373.         err_flag=F_fun_U(U_data+3,F_data_tmp);
  374.         F_data[6]=F_data_tmp[0];F_data[7]=F_data_tmp[1];F_data[8]=F_data_tmp[2];

  375.         //计算F(n)(j+1)
  376.         err_flag=F_fun_U(U_data+9,F_data_tmp);
  377.         F_data[9]=F_data_tmp[0];F_data[10]=F_data_tmp[1];F_data[11]=F_data_tmp[2];
  378.         
  379.         //计算F(n)(j-1)
  380.         err_flag=F_fun_U(U_data+12,F_data_tmp);
  381.         F_data[12]=F_data_tmp[0];F_data[13]=F_data_tmp[1];F_data[14]=F_data_tmp[2];
  382.                
  383.         //计算U(n+1)(j)_star
  384.         U_data[0]=U_data[3]-CFL1*(F_data[9] -F_data[6]);
  385.         U_data[1]=U_data[4]-CFL1*(F_data[10]-F_data[7]);
  386.         U_data[2]=U_data[5]-CFL1*(F_data[11]-F_data[8]);

  387.         //计算U(n+1)(j-1)_star
  388.         U_data[15]=U_data[12]-CFL1*(F_data[6]-F_data[12]);
  389.         U_data[16]=U_data[13]-CFL1*(F_data[7]-F_data[13]);
  390.         U_data[17]=U_data[14]-CFL1*(F_data[8]-F_data[14]);
  391.         
  392.         //计算F(n+1)(j)_star
  393.         err_flag=F_fun_U(U_data+0,F_data_tmp);
  394.         F_data[0]=F_data_tmp[0];F_data[1]=F_data_tmp[1];F_data[2]=F_data_tmp[2];
  395.         
  396.         //计算F(n+1)(j-1)_star
  397.         err_flag=F_fun_U(U_data+15,F_data_tmp);
  398.         F_data[3]=F_data_tmp[0];F_data[4]=F_data_tmp[1];F_data[5]=F_data_tmp[2];
  399.         
  400.         //计算U(n+1)(j)
  401.         U_data[6]=0.5*(U_data[3]+U_data[0])-0.5*CFL1*(F_data[0]-F_data[3]);
  402.         U_data[7]=0.5*(U_data[4]+U_data[1])-0.5*CFL1*(F_data[1]-F_data[4]);
  403.         U_data[8]=0.5*(U_data[5]+U_data[2])-0.5*CFL1*(F_data[2]-F_data[5]);

  404.         //返回U_star_star
  405.         res[0]=U_data[6];//dens
  406.         res[1]=U_data[7];//velc*dens
  407.         res[2]=U_data[8];//ener
  408. }

  409. //用rou,u,p,gama表达e
  410. double ener_comp(double pres_0,double velo_0,double dens_0)
  411. {
  412.         return (pres_0/(gama_ideal_gas-1))+(0.5*dens_0*velo_0*velo_0);
  413. }

  414. //用rou,rou*u,e,gama表达p
  415. double pres_comp(double ener_1,double velo_multi_dens_1,double dens_1)
  416. {
  417.         return (gama_ideal_gas-1)*(ener_1-0.5*velo_multi_dens_1*velo_multi_dens_1/dens_1);
  418. }

  419. //用rou,u,p,gama,e表达U
  420. long U_norm(double dens2,double velc2,double ener2,double pres2,double *U_data2)
  421. {
  422.         U_data2[0]=dens2;
  423.         U_data2[1]=dens2*velc2;
  424.         U_data2[2]=ener_comp(pres2,velc2,dens2);
  425.         return 0;
  426. }

  427. //F_fun_U是以U表达F的函数
  428. long F_fun_U(double *U_data3,double *F_data3)
  429. {
  430.         double P_tmp=(gama_ideal_gas-1)*(U_data3[2]-0.5*U_data3[1]*U_data3[1]/U_data3[0]);
  431.         F_data3[0]=U_data3[1];
  432.         F_data3[1]=U_data3[1]*U_data3[1]/U_data3[0]+P_tmp;
  433.         F_data3[2]=U_data3[1]/U_data3[0]*(U_data3[2]+P_tmp);
  434.         return 0;
  435. }
复制代码


来自:http://swell2005.yculblog.com/post.1267582.html
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-12-4 18:44 , Processed in 0.060558 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表