水流界面追踪VOF模型的CPP源码

#include “Dendritic.h“

#define min（ab） （（（a） < （b）） ? （a） : （b））
#define max（ab） （（（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 f1f2f3f4;
double u1u2u3u4;

double rnxrny;

int itype;
//c   distance of point（x1y1） and point（x2y2）
double dist（double x1double y1double x2double y2） 
{              
      return sqrt（（x1-x2）*（x1-x2）+（y1-y2）*（y1-y2））;
}


//c    calculate the transportation from neighbour cell       
void	transport（double cint itypedouble tanalfadouble tanbeta \
				  double u1double u3double u4double u2 double dxdouble dy\
				  double dtdouble f1double f3double f4double f2）

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

	double s1s2s3s4temp1;
	
	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*cot