young.cpp

VOF追踪界面cpp语言的实现。 VOF追踪界面cpp语言的实现。

#include "Dendritic.h"

#define min(a,b) (((a) < (b)) ? (a) : (b))
#define max(a,b) (((a) > (b)) ? (a) : (b))
const	int isize=100;
const double pi=3.14159265;



double supt=0.5;

double h=1.0/isize;  
double dt=0.1*h;
double dx=h;
double dy=h;
double emikro=1.0e-10;
double emk2=0;
double x_0=0.5;
double y_0=0.3;
double r1=0.2;
double r2=0.15;


double rnx1;
double rny1;
double tanbeta;
double cotbeta;
double tanalfa;
double cotalfa;


double c0[isize][isize];
double c1[isize][isize];
double u[isize][isize];
double v[isize][isize];
double x[isize];
double y[isize]; 

double ft=0;
double fb=0;
double fl=0;
double fr=0 ;
double c;
double ut;
double ub;
double ul;
double ur;

double f1,f2,f3,f4;
double u1,u2,u3,u4;

double rnx,rny;

int itype;
//c   distance of point(x1,y1) and point(x2,y2)
double dist(double x1,double y1,double x2,double y2) 
{              
      return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}


//c    calculate the transportation from neighbour cell       
void	transport(double c,int itype,double tanalfa,double tanbeta, \
				  double u1,double u3,double u4,double u2, double dx,double dy,\
				  double dt,double f1,double f3,double f4,double f2)

{
	double cotalfa=1.0/tanalfa;
	double cotbeta=1.0/tanbeta;  
	f1=0;
	f3=0;
	f4=0;
    f2=0; 

	double s1,s2,s3,s4,temp1;
	
	if(itype==1)
	{ 

        s1=0;
        s3=sqrt(2.0*c*cotalfa);
        s4=0;
        s2=sqrt(2.0*c*tanalfa);

        if(u1>0) 
		{
			if(u1*dt<=(1.0-s2)*dy)
			{
				f1=0;
			}
			else
			{
				temp1=u1*dt-(1.0-s2)*dy;
				f1=0.5*temp1*temp1*cotbeta;
			}
        }  

        if(u2>0)
		{
			if(u2*dt>=s3*dx)
			{
				f2=c*dx*dy;
			}
			else
			{
				f2=0.5*u2*dt*(2.0-u2*dt/(s3*dx))*s2*dy;
			} 
        }
		if(u3<0)
		{
			if(fabs(u3)*dt>=s2*dy) 
			{
				f3=c*dx*dy;
			}
			else
			{
				f3=0.5*fabs(u3)*dt*(2.0-fabs(u3)*dt/(s2*dy))*s3*dx;
			}
        }  
	
		if(u4<0)
		{
			if(fabs(u4)*dt<=(1.0-s3)*dx) 
			{
				f4=0;
			}
			else   
			{
				temp1=fabs(u4)*dt-(1.0-s3)*dx;
				f4=0.5*temp1*temp1*tanbeta;
			}
		}
	}
	else if(itype==2)
	{
        s1=0;
		s3=1.0;
        s4=c-0.5*tanalfa;
        s2=c+0.5*tanalfa;
        if(u1>0)
		{
			if(u1*dt<=(1.0-s2)*dy)
			{      
				f1=0;
			}
			else if(u1*dt<=(1.0-s4)*dy)
			{
				temp1=u1*dt-(1.0-s2)*dy;
				f1=0.5*temp1*temp1*cotbeta;
			}
			else
			{
				f1=u1*dt*dx-(1.0-c)*dx*dy;
			}
        }  
        if(u2>0) 
		{
			f2=u2*dt*(s2*dy-0.5*u2*dt*tanbeta);
        }  
        if(u3<0)
		{
			if(fabs(u3)*dt<=s4*dy) 
			{
				f3=fabs(u3)*dt*dx;
			}
			else if(fabs(u3)*dt<=s2*dy) 
			{ 
				temp1=fabs(u3)*dt-s4*dy;
				f3=fabs(u3)*dt*dx-0.5*temp1*temp1*cotbeta;
			}
			else
			{
				f3=c*dx*dy;
			}
        }
        if(u4<0)
		{
			f4=fabs(u4)*dt*(s4*dy+0.5*fabs(u4)*dt*tanbeta);
        }
	}
	
	else if(itype==3)
	{
        s1=c-0.5*cotalfa;
        s3=c+0.5*cotalfa;
        s4=0;
        s2=1.0;
        if(u1>0) 
		{
          f1=u1*dt*(s1*dx+0.5*u1*dt*cotbeta);
        }
        if(u2>0) 
		{
          if(u2*dt<=s1*dx)
		  {
            f2=u2*dt*dy;
		  }
          else if(u2*dt<=s3*dx)
		  {
            temp1=u2*dt-s1*dx;
            f2=u2*dt*dy-0.5*temp1*temp1*tanbeta;
		  }
          else
		  {
            f2=c*dx*dy;
          }
        }
        if(u3<0) 
		{
          f3=fabs(u3)*dt*(s3*dx-0.5*fabs(u3)*dt*cotbeta);
        }
        if(u4<0) 
		{
          if(fabs(u4)*dt<=(1.0-s3)*dx)
		  {
            f4=0;
		  }
          else if(fabs(u4)*dt<=(1.0-s1)*dx)
		  {
            temp1=fabs(u4)*dt-(1.0-s3)*dx;
            f4=0.5*temp1*temp1*tanbeta;
		  }
          else
		  {
            f4=fabs(u4)*dt*dy-(1.0-c)*dx*dy;
          }
        }
	}
	else if(itype==4) 
	{
        s1=1.0-sqrt(2.0*(1.0-c)*cotalfa);
        s3=1.0;
        s4=1.0-sqrt(2.0*(1.0-c)*tanalfa);
        s2=1.0;
        if(u1>0)
		{
			if(u1*dt>=(1.0-s4)*dy)
			{
				f1=u1*dt*dx-(1.0-c)*dx*dy;
			}
			else  
			{
				f1=u1*dt*(s1*dx+0.5*u1*dt*cotbeta);
			}
        }
        if(u2>0)
		{
			if(u2*dt<=s1*dx) 
			{
				f2=u2*dt*dy;
			}
			else      
			{
				temp1=u2*dt-s1*dx;
				f2=u2*dt*dy-0.5*tanbeta*temp1*temp1;
			}
        }
        if(u3<0) 
		{
			if(fabs(u3)*dt<=s4*dy)
			{
				f3=fabs(u3)*dt*dx;
			}
			else   
			{
				temp1=fabs(u3)*dt-s4*dy;
				f3=fabs(u3)*dt*dx-0.5*temp1*temp1*cotbeta;
			}
        }
            
        if(u4<0) 
		{
          if(fabs(u4)*dt>=(1.0-s1)*dx) 
		  {
            f4=fabs(u4)*dt*dy-(1.0-c)*dx*dy;
		  }

          else
		  {
            f4=fabs(u4)*dt*(s4*dy+0.5*fabs(u4)*dt*tanbeta);
          }
        }
	  }
      else
	  {
	  
        ;//write(*,*)'error in subroutine transport';
	  }
				
       
      
  
      
}
	 
//c     program main    
//c    solve fluid volume function equation using Youngs' VOF method       
//     program main


void WriteSuperbee(int iter)
{
	//c.....local variables
	int  i,j;
	
	
	FILE    *out;
	char    fn[256];
	sprintf(fn, "Res%05d.plt", iter);
    out = fopen(fn,"wt");
	
	fprintf(out,"TITLE     = \" Superbee \"\n");
	fprintf(out,"VARIABLES = \"X\", \"Y\", \"U\", \"V\", \"C0\",\n");
	fprintf(out,"ZONE I= %d , J= %d , F=POINT\n",isize,isize);
	

	
	for(i = 0;i<isize;i++)
	{
        for(j = 0;j<isize; j++)
		{
			
			fprintf(out,"%d\t%d\t%g\t%g\t%f\n",i,j,u[i][j],v[i][j],c0[i][j]);	
			
			
		}
	}
	fclose(out);
	
	
}
    
void main()
{

 
      


	  int i,j;
	  for(i=0;i<isize;i++)
	  {
		  x[i]=dx*i;
		  y[i]=dy*i;
	  }
   	  
      for( i=0;i<isize;i++)
	  {
		  for( j=0;j<isize;j++)
		  {

			  u[i][j]=-pi*cos(pi*(x[i]-0.5))*sin(pi*(y[i]-0.5));
			  v[i][j]= pi*sin(pi*(x[i]-0.5))*cos(pi*(y[i]-0.5)); 

			  //case 1:
			  u[i][j]=-pi*(y[j]-0.5);
			  v[i][j]= pi*(x[i]-0.5);
		  }
	  }

	  double radius=0.4;

	  for( i=0;i<isize;i++)
	  {
		  for( j=0;j<isize;j++)
		  {

			  if(dist(x[i],y[j],x_0,y_0)<=r1)
			  {
				  c0[i][j]=1.0;
			  }
			  else
			  {			 
				  c0[i][j]=0;
			  }	
			  

			  //case 1:
			  if(dist(x[i],y[j],0.5,0.5) < radius)
			  {
				  c0[i][j]=1.0;
			  }
			  else
			  {
				  c0[i][j]=0.0;
			  }
			  
			  if(y[j]<0.5 && x[i]>0.425 && x[i]<0.575)
			  {
				  c0[i][j]=0.0;
			  }


		  }
	  }
    


	  WriteSuperbee(0);
	  

 //   do ll=1,2 //???

      double t=0;
      double it=0;
	  while(t<supt)
	  {
		  
		  for( i=0;i<isize;i++)
		  {
			  for( j=0;j<isize;j++)
			  {
				  c1[i][j]=c0[i][j];
			  }
		  }
		  
		  for( i=1;i<isize-1;i++)
		  {
			  for( j=1;j<isize-1;j++)
			  { 
				  ft=0;
				  fb=0; 
				  fl=0;
				  fr=0 ;
				  c=c0[i][j];  
				  ut=v[i][j];
				  ub=c1[i][j+1];
				  ul=u[i-1][j];
				  ur=u[i][j]; 
				  
				  if(c>=1.0-emikro)
				  {					  
					  if(ut>0) 
					  {
						  ft=ut*dt*dx*c;
					  }  
					  if(ub<0)
					  {
						  fb=fabs(ub)*dt*dx*c;
					  }
					  if(ul<0)
					  {
						  fl=fabs(ul)*dt*dy*c;
					  }
					  if(ur>0) 
					  {
						  fr=ur*dt*dy*c;
					  }
				  } 
				  else if(c>0) 
				  {
					  rnx=(c0[i+1][j+1]+2.0*c0[i+1][j]+c0[i+1][j-1] \
						  -c0[i-1][j+1]-2.0*c0[i-1][j]-c0[i-1][j-1])/dx;
					 
					  rny=(c0[i+1][j+1]+2.0*c0[i][j+1]+c0[i-1][j+1] \
						  -c0[i+1][j-1]-2.0*c0[i][j-1]-c0[i-1][j-1])/dy;			  
					  
					  if(fabs(rny)<=emk2)
					  {					  
						  if(rnx>=0) 
						  {
							  u1=ut;
							  u2=ur;
							  u3=ub;
							  u4=ul;
						  }
						  else
						  {
							  u1=ut;
							  u2=-ul;   
							  u3=ub;
							  u4=-ur;
						  } 
						  
						  f1=0;
						  f2=0;
						  f3=0;
						  f4=0  ;
						  
						  if(u1>0)
						  {
							  f1=fabs(u1)*dt*c*dx;
						  }
						  if(u3<0) 
						  {
							  f3=fabs(u3)*dt*c*dx;
						  }
						  if(u2>0) 
						  {
							  if(fabs(u2)*dt<=c*dx) 
							  {
								  f2=fabs(u2)*dt*dy;
							  }
							  else
							  {
								  f2=c*dx*dy;
							  }
						  }
						  
						  if(u4<0) 
						  {
							  if(fabs(u4)*dt <= (1.0-c)*dx)
							  {
								  f4=0;
							  }
							  else
							  {
								  f4=(fabs(u4)*dt-(1.0-c)*dx)*dy;
							  }
						  }
						  
						  if(rnx>=0) 
						  {
							  ft=f1;   
							  fb=f3;
							  fr=f2;
							  fl=f4;
						  }
						  else
						  {
							  ft=f1;
							  fb=f3;
							  fr=f4;
							  fl=f2;
						  }
					  }
					  else if(fabs(rnx)<=emk2) 
					  {
						  if(rny>=0)  
						  {
							  u1=ut;
							  u2=ur;
							  u3=ub;
							  u4=ul;
						  }
						  else
						  {
							  u1=-ub;
							  u2=ur;
							  u3=-ut;
							  u4=ul;
						  }  
						  f1=0;
						  f2=0;
						  f3=0;
						  f4=0;  
						  if(u1>0)
						  {
							  if(u1*dt<=c*dy) 
							  {
								  f1=u1*dt*dx;
							  }
							  else
							  {
								  f1=c*dx*dy;								  
							  }
						  } 
						  if(u3<0) 
						  {
							  if(fabs(u3)*dt<=(1.0-c)*dy)
							  {
								  f3=0;
							  }
							  else
							  {
								  f3=(fabs(u3)*dt-(1.0-c)*dy)*dx;
							  }  
						  }
						  if(u2>0) 
						  {
							  f2=u2*dt*c*dy;
						  }
						  if(u4<0) 
						  {
							  f4=fabs(u4)*dt*c*dy;
						  }    
						  
						  if(rny>=0)
						  {
							  ft=f1; 
							  fb=f3;
							  fr=f2;
							  fl=f4;
						  }
						  else
						  {							  
							  ft=f3;
							  fb=f1;
							  fr=f2;
							  fl=f4;
						  }
					  }
					  else
					  {
						  if(rnx>0&&rny<0) 
						  {
							  u1=ut;
							  u2=ur;
							  u3=ub;
							  u4=ul;
						  }	
						  else if(rnx<0&&rny>0)
						  {
							  u1=-ub;
							  u2=-ul;
							  u3=-ut;
							  u4=-ur;
						  }
						  else if(rnx<0&&rny<0) 
						  {
							  u1=ut;
							  u2=-ul;
							  u3=ub;
							  u4=-ur; 
						  }
						  else
						  {
							  u1=-ub;
							  u2=ur;
							  u3=-ut;
							  u4=ul;
						  }
						  rnx1=fabs(rnx);
						  rny1=-fabs(rny) ;
						  tanbeta=-rnx1/rny1;
						  cotbeta=1.0/tanbeta;
						  tanalfa=dx/dy*tanbeta;
						  cotalfa=1.0/tanalfa;
						  
						  if(tanalfa<=1.0)
						  {							  
							  if(c<=0.5*tanalfa)
							  {
								  itype=1;
							  } 
							  else if(c<=1.0-0.5*tanalfa)
							  {
								  itype=2;
							  }
							  else
							  {
								  itype=4;
							  }
						  }
						  else
						  {
							  if(c<=0.5*cotalfa)
							  {
								  itype=1;
							  }
							  else if(c<=1.0-0.5*cotalfa)
							  {
								  itype=3;
							  }
							  else
							  {
								  itype=4;
							  }
						  }
					  }
					  transport(c,itype,tanalfa,tanbeta,u1,u3,u4,u2, dx,dy,dt,f1,f3,f4,f2);
					  
					  if(rnx>0&&rny<0) 
					  {
						  ft=f1;
						  fb=f3;
						  fl=f4;
						  fr=f2;
					  }
					  else if(rnx<0&&rny>0) 
					  {
						  ft=f3;
						  fb=f1;
						  fl=f2;
						  fr=f4;
					  }
					  else if(rnx<0&&rny<0)
					  {
						  ft=f1;
						  fb=f3;
						  fl=f2; 
						  fr=f4 ;
					  }
					  else
					  {
						  ft=f3;
						  fb=f1;
						  fl=f4;
						  fr=f2;
					  }
              }		  
			  c1[i][j+1]=c1[i][j+1]+ft/dx/dy;
			  c1[i][j-1]=c1[i][j-1]+fb/dx/dy;
			  c1[i+1][j]=c1[i+1][j]+fr/dx/dy;
			  c1[i-1][j]=c1[i-1][j]+fl/dx/dy;
			  c1[i][j]=c1[i][j]-(ft+fb+fl+fr)/dx/dy;
          } //...i end
        }  //...j end
                  
                 

        for( i=1;i<isize;i++)
		{
          for( j=1;j<isize;j++)
		  {
			  c0[i][j]=max(min(c1[i][j],1.0),0);
          }
        }

        t=t+dt;
        it=it+1;
      }  //..  t end
  

   /*
      for( i=0;i<isize;i++)
	  {	  
        for( j=0;j<isize;j++)
		{
          u[i][j]=pi*cos(pi*(x[i]-0.5))*sin(pi*(y[i]-0.5));
          v[i][j]=-pi*sin(pi*(x[i]-0.5))*cos(pi*(y[i]-0.5)) ;
        }
      }
*/
	  WriteSuperbee(1);
	  
}