一维激波管的maccormack格式模拟

gghhjj

1. // CFD_1D_1st.cpp
   
2. // 一维激波管问题，采用二阶MacCormark格式，加入人工粘性以保证震荡稳定
   
3. //
   
4. // 最后修正2006.5.8
   
5. 
   
6. # include <math.h>
   
7. # include <iostream.h>
   
8. # include <iomanip.h>//引用setw()函数
   
9. //Max函数&Min函数，用于处理误差，保证健壮性
   
10. # define max_abs(a,b) ((fabs(a) > fabs(b)) ? (a) : (b))
    
11. # define min_abs(a,b) ((fabs(a) < fabs(b)) ? (a) : (b))
    
12. //常量定义
    
13. # define gama_ideal_gas 1.4//气体动力学中的常数，与气体类型有关，本值是相对理想气体而言
    
14. # define tolerance 1e-6//定义计算中的误差大小，小于之将被视为零，以保证程序的健壮性
    
15. # define L 0.2//经验的人工粘性系数，一般取0.6-3.0之间
    
16. //符号说明dens=density,velc=velocity,ener=energy,pres=pressure
    
17. 
    
18. void U_star_star (double *,double *,double *,double *,double *,double);
    
19. double ener_comp(double,double,double);
    
20. double pres_comp(double,double,double);
    
21. long U_norm(double,double,double,double,double *);
    
22. long F_fun_U(double *,double *);
    
23. 
    
24. int main ()
    
25. {
    
26.         // 网格参数 x_range[3],t_range[2],Stmp_X,CFL or Stmp_T
    
27.         const double Step_X=0.002,CFL=0.5,Step_T=Step_X*CFL;
    
28.         const double x_range[3]={-1,0,1};
    
29.         const double t_range[3]={0,0,0.3};
    
30.         
    
31.         const long MESH_X_P1=(long) floor( (x_range[1]-x_range[0])/Step_X )+1;
    
32.         const long MESH_X_P2=(long) floor( (x_range[2]-x_range[1])/Step_X )+1;
    
33.         const long MESH_X       =MESH_X_P1+MESH_X_P2;
    
34.         const long MESH_T       =(long) floor( (t_range[2]-t_range[0])/Step_T )+2;
    
35.         
    
36.         // 求解变量 dens[MESH_T][MESH_X],velc[][],ener[][],pres[][]
    
37.         double *dens=new double[ MESH_T * MESH_X ];
    
38.         double *velc=new double[ MESH_T * MESH_X ];
    
39.         double *ener=new double[ MESH_T * MESH_X ];
    
40.         double *pres=new double[ MESH_T * MESH_X ];
    
41. 
    
42.         //初始条件
    
43.         for(long i=0;i<MESH_X_P1;i++)
    
44.         {
    
45.                 //第一层计算区域网格
    
46.                 dens [i+MESH_X]=1;
    
47.                 velc [i+MESH_X]=max_abs(0,tolerance);
    
48.                 pres [i+MESH_X]=1;
    
49.                 ener [i+MESH_X]=max_abs(ener_comp(pres[i+MESH_X],velc[i+MESH_X],dens [i+MESH_X]),tolerance);
    
50.                 //第一层冗余网格
    
51.                 dens [i]=1;
    
52.                 velc [i]=max_abs(0,tolerance);
    
53.                 pres [i]=1;
    
54.                 ener [i]=max_abs(ener_comp(pres[i],velc[i],dens [i]),tolerance);
    
55.         }
    
56.         for(i=0;i<MESH_X_P2;i++)
    
57.         {
    
58.                 //第一层计算区域网格
    
59.                 dens [i+MESH_X_P1+MESH_X]=0.125;
    
60.                 velc [i+MESH_X_P1+MESH_X]=max_abs(0,tolerance);
    
61.                 pres [i+MESH_X_P1+MESH_X]=0.1;
    
62.                 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);
    
63.                 //第一层冗余网格
    
64.                 dens [i+MESH_X_P1]=0.125;
    
65.                 velc [i+MESH_X_P1]=max_abs(0,tolerance);
    
66.                 pres [i+MESH_X_P1]=0.1;
    
67.                 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);
    
68.         }
    
69. 
    
70.         //计算
    
71.         
    
72.         double *U_2_star=new double [3*MESH_X];//光滑前的U(n+1)(j)
    
73.         double Q;//开关函数，源于人工粘性滤波算法
    
74.         for(long n=1;n<MESH_T-1;n++)
    
75.         {      
    
76.                 if(n>1)
    
77.                 {       //开头与结尾的冗余网格，取前时刻本位置左（右）平均
    
78.                         dens [n*MESH_X]=max_abs(0.5*(dens [(n-1)*MESH_X]+dens [(n-1)*MESH_X+1]),tolerance);
    
79.                         velc [n*MESH_X]=max_abs(0.5*(velc [(n-1)*MESH_X]+velc [(n-1)*MESH_X+1]),tolerance);
    
80.                         pres [n*MESH_X]=max_abs(0.5*(pres [(n-1)*MESH_X]+pres [(n-1)*MESH_X+1]),tolerance);
    
81.                         ener [n*MESH_X]=max_abs(0.5*(ener [(n-1)*MESH_X]+ener [(n-1)*MESH_X+1]),tolerance);
    
82. 
    
83.                         dens [(n+1)*MESH_X-1]=max_abs(0.5*(dens [n*MESH_X-1]+dens [n*MESH_X-2]),tolerance);
    
84.                         velc [(n+1)*MESH_X-1]=max_abs(0.5*(velc [n*MESH_X-1]+velc [n*MESH_X-2]),tolerance);
    
85.                         pres [(n+1)*MESH_X-1]=max_abs(0.5*(pres [n*MESH_X-1]+pres [n*MESH_X-2]),tolerance);
    
86.                         ener [(n+1)*MESH_X-1]=max_abs(0.5*(ener [n*MESH_X-1]+ener [n*MESH_X-2]),tolerance);
    
87.                 }
    
88.                
    
89.                 //计算U_star_star向量
    
90.                 for(long j=1;j<MESH_X-1;j++)
    
91.                 {
    
92.                         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);
    
93.                 }
    
94.                
    
95.                 //U(n+1)(j)计算，由U_star_star的函数表达，抹平U_star_star的震荡
    
96.                 //其中U(n+1)(1)和U(n+1)(MESH_X-2)为两边边界，不可能有激波，不用抹平
    
97.                 //U(n+1)(1)
    
98.                         dens [(n+1)*MESH_X+1]=max_abs(U_2_star[3],tolerance);
    
99.                         velc [(n+1)*MESH_X+1]=max_abs(U_2_star[4]/U_2_star[3],tolerance);
    
100.                         ener [(n+1)*MESH_X+1]=max_abs(U_2_star[5],tolerance);
     
101.                         pres [(n+1)*MESH_X+1]=max_abs(pres_comp(U_2_star[5],U_2_star[4],U_2_star[3]),tolerance);
     
102.                 //U(n+1)(MESH_X-1)
     
103.                         dens [(n+2)*MESH_X-2]=max_abs(U_2_star[3*(MESH_X-2)],tolerance);
     
104.                         velc [(n+2)*MESH_X-2]=max_abs(U_2_star[3*(MESH_X-2)+1]/U_2_star[3*(MESH_X-2)],tolerance);
     
105.                         ener [(n+2)*MESH_X-2]=max_abs(U_2_star[3*(MESH_X-2)+2],tolerance);
     
106.                         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);
     
107.                 //U(n+1)(j)
     
108.                 for(j=2;j<MESH_X-2;j++)
     
109.                 {
     
110.                         //开关函数Q
     
111.                         //MacCormack定义的开关函数，比Harten好用，但仍会震荡，实际中用l限制其范围，有效果，但不明显
     
112.                         //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 ) );
     
113.                         //Q=min_abs(fabs(Q),L);Q=((fabs(Q) < fabs(tolerance)) ? (0) : (Q));//将Q限制到0-L之间
     
114.                         Q=L;
     
115.                         //Harten定义的开关函数，不好用，会出现无穷
     
116.                         //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]) );
     
117.                         //cout<<Q<<"            "<<endl;用于Q调错
     
118.                         double final_res[3];//U(n+1)(j)数据的临时保存
     
119.                         for(short t=0;t<3;t++)//计算光滑过的U(n+1)(j)
     
120.                         {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]);}
     
121.                         //返回所求n+1时刻的dens,velc,ener,pres
     
122.                         dens [(n+1)*MESH_X+j]=max_abs(final_res[0],tolerance);
     
123.                         velc [(n+1)*MESH_X+j]=max_abs(final_res[1]/final_res[0],tolerance);
     
124.                         ener [(n+1)*MESH_X+j]=max_abs(final_res[2],tolerance);
     
125.                         pres [(n+1)*MESH_X+j]=max_abs(pres_comp(final_res[2],final_res[1],final_res[0]),tolerance);                     
     
126.                 }
     
127.         }
     
128.         
     
129.         //结果输出
     
130.         for(n=1;n<MESH_T-1;n++)
     
131.         {
     
132.                 cout<<n<<"   时刻"<<endl;
     
133.                 for(long j=1;j<MESH_X-1;j++)
     
134.                 {
     
135.                         if(fabs(dens[n*MESH_X+j])<=tolerance){dens[n*MESH_X+j]=0;}
     
136.                         if(fabs(velc[n*MESH_X+j])<=tolerance){velc[n*MESH_X+j]=0;}
     
137.                         if(fabs(pres[n*MESH_X+j])<=tolerance){pres[n*MESH_X+j]=0;}
     
138.                         if(fabs(ener[n*MESH_X+j])<=tolerance){ener[n*MESH_X+j]=0;}
     
139. 
     
140.                   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;
     
141.                 }
     
142.                 cout<<""<<endl;
     
143.         }
     
144.         return 0;
     
145. }
     
146. 
     
147. void U_star_star (double *dens1,double *velc1,double *ener1,double *pres1,double *res,double CFL1)
     
148. {
     
149.         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
     
150.         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)
     
151.         double U_data_tmp[3];//临时存放从U_norm()拿到的数据
     
152.         double F_data_tmp[3];//临时存放从F_fun_U()拿到的数据
     
153.         
     
154.         //计算U(n)(j),数据临时存在U_data_tmp里
     
155.         long err_flag=U_norm(*dens1,*velc1,*ener1,*pres1,U_data_tmp);
     
156.         U_data[3]=U_data_tmp[0];U_data[4]=U_data_tmp[1];U_data[5]=U_data_tmp[2];
     
157. 
     
158.         //计算U(n)(j+1),数据临时存在U_data_tmp里
     
159.         err_flag=U_norm(*(dens1+1),*(velc1+1),*(ener1+1),*(pres1+1),U_data_tmp);
     
160.         U_data[9]=U_data_tmp[0];U_data[10]=U_data_tmp[1];U_data[11]=U_data_tmp[2];
     
161. 
     
162.         //计算U(n)(j-1),数据临时存在U_data_tmp里
     
163.         err_flag=U_norm(*(dens1-1),*(velc1-1),*(ener1-1),*(pres1-1),U_data_tmp);
     
164.         U_data[12]=U_data_tmp[0];U_data[13]=U_data_tmp[1];U_data[14]=U_data_tmp[2];
     
165. 
     
166.         //计算F(n)(j)
     
167.         err_flag=F_fun_U(U_data+3,F_data_tmp);
     
168.         F_data[6]=F_data_tmp[0];F_data[7]=F_data_tmp[1];F_data[8]=F_data_tmp[2];
     
169. 
     
170.         //计算F(n)(j+1)
     
171.         err_flag=F_fun_U(U_data+9,F_data_tmp);
     
172.         F_data[9]=F_data_tmp[0];F_data[10]=F_data_tmp[1];F_data[11]=F_data_tmp[2];
     
173.         
     
174.         //计算F(n)(j-1)
     
175.         err_flag=F_fun_U(U_data+12,F_data_tmp);
     
176.         F_data[12]=F_data_tmp[0];F_data[13]=F_data_tmp[1];F_data[14]=F_data_tmp[2];
     
177.                
     
178.         //计算U(n+1)(j)_star
     
179.         U_data[0]=U_data[3]-CFL1*(F_data[9] -F_data[6]);
     
180.         U_data[1]=U_data[4]-CFL1*(F_data[10]-F_data[7]);
     
181.         U_data[2]=U_data[5]-CFL1*(F_data[11]-F_data[8]);
     
182. 
     
183.         //计算U(n+1)(j-1)_star
     
184.         U_data[15]=U_data[12]-CFL1*(F_data[6]-F_data[12]);
     
185.         U_data[16]=U_data[13]-CFL1*(F_data[7]-F_data[13]);
     
186.         U_data[17]=U_data[14]-CFL1*(F_data[8]-F_data[14]);
     
187.         
     
188.         //计算F(n+1)(j)_star
     
189.         err_flag=F_fun_U(U_data+0,F_data_tmp);
     
190.         F_data[0]=F_data_tmp[0];F_data[1]=F_data_tmp[1];F_data[2]=F_data_tmp[2];
     
191.         
     
192.         //计算F(n+1)(j-1)_star
     
193.         err_flag=F_fun_U(U_data+15,F_data_tmp);
     
194.         F_data[3]=F_data_tmp[0];F_data[4]=F_data_tmp[1];F_data[5]=F_data_tmp[2];
     
195.         
     
196.         //计算U(n+1)(j)
     
197.         U_data[6]=0.5*(U_data[3]+U_data[0])-0.5*CFL1*(F_data[0]-F_data[3]);
     
198.         U_data[7]=0.5*(U_data[4]+U_data[1])-0.5*CFL1*(F_data[1]-F_data[4]);
     
199.         U_data[8]=0.5*(U_data[5]+U_data[2])-0.5*CFL1*(F_data[2]-F_data[5]);
     
200. 
     
201.         //返回U_star_star
     
202.         res[0]=U_data[6];//dens
     
203.         res[1]=U_data[7];//velc*dens
     
204.         res[2]=U_data[8];//ener
     
205. }
     
206. 
     
207. //用rou,u,p,gama表达e
     
208. double ener_comp(double pres_0,double velo_0,double dens_0)
     
209. {
     
210.         return (pres_0/(gama_ideal_gas-1))+(0.5*dens_0*velo_0*velo_0);
     
211. }
     
212. 
     
213. //用rou,rou*u,e,gama表达p
     
214. double pres_comp(double ener_1,double velo_multi_dens_1,double dens_1)
     
215. {
     
216.         return (gama_ideal_gas-1)*(ener_1-0.5*velo_multi_dens_1*velo_multi_dens_1/dens_1);
     
217. }
     
218. 
     
219. //用rou,u,p,gama,e表达U
     
220. long U_norm(double dens2,double velc2,double ener2,double pres2,double *U_data2)
     
221. {
     
222.         U_data2[0]=dens2;
     
223.         U_data2[1]=dens2*velc2;
     
224.         U_data2[2]=ener_comp(pres2,velc2,dens2);
     
225.         return 0;
     
226. }
     
227. 
     
228. //F_fun_U是以U表达F的函数
     
229. long F_fun_U(double *U_data3,double *F_data3)
     
230. {
     
231.         double P_tmp=(gama_ideal_gas-1)*(U_data3[2]-0.5*U_data3[1]*U_data3[1]/U_data3[0]);
     
232.         F_data3[0]=U_data3[1];
     
233.         F_data3[1]=U_data3[1]*U_data3[1]/U_data3[0]+P_tmp;
     
234.         F_data3[2]=U_data3[1]/U_data3[0]*(U_data3[2]+P_tmp);
     
235.         return 0;
     
236. }

复制代码