sola.cpp

sola方法的源代码

#include <stdio.h>
#include <iostream.h>
#include <math.h>

                                  //注意，取计算区域的下半部分，根据对称性可以算出全流场 
#define Re  100       //雷诺数 
#define I   80    //x方向空间步数         
#define J   20     //y方向空间步数         
#define dt  0.015     //时间步长              
#define alpha 0.4    //alpha系数             
#define D0  1.0       //D参考值
#define epsl 1e-4     //判断收敛的小数epsilon
#define om   1.8      //超松弛因子omega       

static double   dx=4.0/I;
static double   dy=1.0/J;
static int  i,j,k,K,flag,M=11,P=20,Q=40;    //M为BC上y上的纵坐标 ；P为C点对应的横坐标、Q为E点对应的横坐标,flag为n+1时间层全流场未收敛点数标志
static double u[I][J],u1[I][J],v[I][J],v1[I][J],p[I][J],D[I][J],dp[I][J],max=0,x,y;
static double  FUX,FUY,FVX,FVY,VISX,VISY,thet=om/2/dt/(1/dx/dx+1/dy/dy);  //thet:thet因子

void initialize();
void margin();
void comput();      //求u1(i,j),v1(i,j)
void modify();      //修正各值   
void save_data();





void initialize()
{
  for(i=0;i<=I;i++)                         //初始化(流场加边界)
	  for(j=0;j<J;j++)
	  {
		  u[i][j]=0;
		  v[i][j]=0;
		  if(j>=M||j<M&&j>0&&i>=P&&i<=Q)
			  p[i][j]=1;
	  }
}


void margin()
{
     for(j=M;j<J;j++)
  {
	if(j<=int(0.6/dy))                   //输入边界条件AB边 
	   u1[1][j]=u1[0][j]=10*(j*dy-dy/2-0.5);       //y=j*dy-dy/2
   else
	   u1[1][j]=u1[0][j]=1;
       v1[1][j]=v1[0][j]=0;

    u1[I-1][j]=u1[I-2][j];                    //输出边界，u1j=u2j
	v1[I-1][j]=v1[I-2][j];
  }
  for(j=1;j<M;j++)                            //壁面
   {
    u1[P-1][j]=-u1[P][j];                        //CD
	v1[P][j]=0;
	p[P][j]=p[P+1][j]+dp[P+1][j];         //;
	u1[Q][j]=-u1[Q-1][j];                     //EF
	v1[Q][j]=0;
	p[Q][j]=p[Q-1][j]+dp[Q-1][j];          //;
}
for(i=1;i<I;i++)
  {
  v1[i][J-2]=0;                                 //AH
  u1[i][J-1]=u1[i][J-2];
  p[i][J-1]=p[i][J-2];
     if(i<P||i>=Q)                               //BC、FG
   {
	 v1[i][M-1]=0;
     u1[i][M-1]=-u1[i][M];
   }
   if(i>=P&&i<Q)                                 //DE
   {
	v1[i][0]=0;
    u1[i][0]=-u1[i][1];
   }
  }
	                               //拐角点
//	v1[P][0]=0;                    //D
	v1[Q][0]=0;                    //E
}

void comput()
{
  for(i=1;i<I-1;i++)           //对流项和粘性项
      for(j=1;j<J-1;j++)
	  {
		if(j>=M||j<M&&i>P&&i<Q)                //流场里求u1(i,j)    
		{
	 FUX=((u[i][j]+u[i+1][j])*(u[i][j]+u[i+1][j])+alpha*fabs(u[i][j]+u[i+1][j])*(u[i][j]-u[i+1][j])-(u[i-1][j]+u[i][j])*(u[i-1][j]+u[i][j])-alpha*fabs(u[i-1][j]+u[i][j])*(u[i-1][j]-u[i][j]))/(4.0*dx);
     FUY=((v[i][j]+v[i+1][j])*(u[i][j]+u[i][j+1])+alpha*fabs(v[i][j]+v[i+1][j])*(u[i][j]-u[i][j+1])-(v[i][j-1]+v[i+1][j-1])*(u[i][j-1]+u[i][j])-alpha*fabs(v[i][j-1]+v[i+1][j-1])*(u[i][j-1]+u[i][j]))/(4.0*dy);
     VISX=((u[i+1][j]-2*u[i][j]+u[i-1][j])/dx/dx+(u[i][j+1]-2*u[i][j]+u[i][j-1])/dy/dy)/Re;
	   u1[i][j]=u[i][j]+dt*((p[i][j]-p[i+1][j])/dx-FUX-FUY+VISX);
		}
         if(j>=M||j<M&&i>P&&i<Q)                //流场里求v1(i,j)
		 {
     FVX=((u[i][j]+u[i][j+1])*(v[i][j]+v[i+1][j])+alpha*fabs(u[i][j]+u[i][j+1])*(v[i][j]-v[i+1][j])-(u[i-1][j]+u[i-1][j+1])*(v[i-1][j]+v[i][j])-alpha*fabs(u[i-1][j]+u[i-1][j+1])*(v[i-1][j]+v[i][j]))/(4.0*dx);
	 FVY=((v[i][j]+v[i][j+1])*(v[i][j]+v[i][j+1])+alpha*fabs(v[i][j]+v[i][j+1])*(v[i][j]-v[i][j+1])-(v[i][j-1]+v[i][j])*(v[i][j-1]+v[i][j])-alpha*fabs(v[i][j-1]+v[i][j])*(v[i][j-1]-v[i][j]))/(4.0*dy);
     VISY=((v[i+1][j]-2*v[i][j]+v[i-1][j])/dx/dx+(v[i][j+1]-2*v[i][j]+v[i][j-1])/dy/dy)/Re;
       v1[i][j]=v[i][j]+dt*((p[i][j]-p[i][j+1])/dy-FVX-FVY+VISY);
		 }
	  }
}

void modify()
{
      for(i=1;i<I;i++)
	    for(j=J-2;j>0;j--)
		   if(j>=M||j<M&&i>P&&i<Q)         //流场内散度D的点
		   {
	   D[i][j]=(u1[i][j]-u1[i-1][j])/dx+(v1[i][j]-v1[i][j-1])/dy;
	     if(fabs(D[i][j]/D0)>epsl)              //不满足散度条件则修正
		 {
		   flag++;
	    #if 0
		   if(k==40)
	   cout <<"D["<<i<<"]["<<j<<"]="<<D[i][j]<<"  ";
       #endif
	   dp[i][j]=-D[i][j]*thet;
		 }
		 else 
	    dp[i][j]=0;
       p[i][j]+=dp[i][j];
       u1[i-1][j]-=dt*dp[i][j]/dx;    
	   v1[i][j-1]-=dt*dp[i][j]/dy;
	   u1[i][j]+=dt*dp[i][j]/dx;
	   v1[i][j]+=dt*dp[i][j]/dy;
		   }
}

void main()
{
  initialize();                      //对t=0时赋值
         k=0;
    do     
	 {
      if(k)
	  
			if(j>=M||j<M&&i>=P&&i<Q)
		{
	      u[i][j]=u1[i][j];
	      v[i][j]=v1[i][j];
		}

       modify();              //求u(n+1),v(n+1),即：u1,v1
	   ////////////////
	//     循环迭代     //
	   ///////////////
	   do
	   {
		 flag=0;
   	  margin();                 //边界处理
	  modify();	   //求散度D(ij),修正压力p(ij),速度u1(ij),v1(ij)
	  cout <<" now,time-step为"<<k<<", 不满足的点的个数为"<<flag<<'\n';
	   }
        while(flag);
		 for(i=1;i<I;i++)
			 for(j=1;j<J-1;j++)
			   if(j>=M||j<M&&i>=P&&i<Q)                //流场里切向速度u的点
			 {
			 if(fabs(u1[i][j]-u[i][j])>=epsl)
				 if(fabs(u1[i][j]-u[i][j])>max)          //求稳定场  
				 max=fabs(u1[i][j]-u[i][j]);
			 }
		cout <<"误差是"<<max<<'\n';
		k++;
	 	if(k==40) break;
  	}
	 while(max);
	 save_data();

}





		//////////////////////
	//  结果输出   
	   ///////////////////////
void save_data()
{
      for(i=1;i<I;i++)
	   for(j=1;j<J-1;j++) 
		   if(j>=M||j<M&&i>P&&i<Q)         //流场内的速度点，也是散度压力点
		   {
  u[i][j]=(u1[i-1][j]+u1[i][j])/2;      
  v[i][j]=(v1[i][j-1]+v1[i][j])/2;
		   }
		   else
		   {
			   u[i][j]=0;
			   v[i][j]=0;
			   p[i][j]=0;
		   }
		  
FILE *fp;
  fp=fopen("./sola_result.txt","w");        
   fprintf(fp,"VARIABLES = \"X\", \"Y\", \"u\", \"v\",\"Pressure\"\nZONE i=%d, j=%d,F=POINT\n",I,J-2);
for(j=1;j<=J;j++)
  for(i=0;i<I;i++)
	  {
	  x=i*dx;
	  y=j*dy-dy/2;
   fprintf(fp,"%f\t%f\t%f\t%f\t%f\n",x,y,u[i][j],v[i][j],p[i][j]);
	  }
  fclose(fp);
}