simple.c

这是一个用SIMPLE法模拟的矩形流场的程序

#include "global.h"


void simple()
{
	n=1;
	//u方程系数设定
	
	//ae和aw设定
	for(j=2;j<=m-1;j++)
	{//左边界aw[3][j]
		d=dif[n]*dy/dx;
		flow=den*u(2,j)*dy;
		acf();
		aw[3][j][1]=d*adf+max(flow,0);
		//右边界ae[l-1][j]
		flow=den*u(l,j)*dy;
		acf();
		ae[l-1][j][1]=d*adf+max(-1*flow,0);
		for(i=3;i<=l-2;i++)
		{
		//u[i][j]的ae和u[i+1][j]的aw
			d=dif[n]*dy/dx;
			flow=den*dy*(u(i,j)+u(i+1,j))/2;
			acf();
			aw[i+1][j][1]=d*adf+max(flow,0);
			ae[i][j][1]=aw[i+1][j][1]-flow;
		}
	}

	//an和as设定
	for(i=3;i<=l-1;i++)
	{		
		//下边界as和上边界an
		if(i==3)
		{	//as[3][2]
			d=dif[n]*dx*3/dy;//d=dif[n]*(dx*3/2)/(dy/2)
			flow=den*dx*(v(i-1,2)+v(i,2)/2);
			acf();
			as[i][2][1]=d*adf+max(flow,0);
			//an[3][m-1]
			flow=den*dx*(v(i-1,m)+v(i,m)/2);
			acf();
			an[i][m-1][1]=d*adf+max(-1*flow,0);	
		}
			else
				if(i==l-1)
				{	//as[l-1][2]
					d=dif[n]*dx*3/dy;
					flow=den*dx*(v(i-1,2)/2+v(i,2));
					acf();
					as[i][2][1]=d*adf+max(flow,0);
					//an[l-1][m-1]
					flow=den*dx*(v(i-1,m)/2+v(i,m));
					acf();
					an[i][m-1][1]=d*adf+max(-1*flow,0);
				}
			
				else
				{
					d=dif[n]*dx*2/dy;
					//as[i][2]
					flow=den*dx*(v(i-1,2)+v(i,2))/2;
					acf();
					as[i][2][1]=d*adf+max(flow,0);
					//an[i][m-1]
					flow=den*dx*(v(i-1,m)+v(i,m))/2;
					acf();
					an[i][m-1][1]=d*adf+max(-1*flow,0);
					
				}

		//内部边界an[i][j]和as[i][j+1]
		for(j=2;j<=m-2;j++)
		{
			if(i==3)
			{
				d=dif[n]*dx*3/2/dy;
				flow=den*dx*(v(i-1,j+1)+v(i,j+1)/2);
			}
			else
				if(i==l-1)
				{
					d=dif[n]*dx*3/2/dy;
					flow=den*dx*(v(i-1,j+1)/2+v(i,j+1));
				}
				else
				{
					d=dif[n]*dx/dy;
					flow=den*(v(i-1,j+1)+v(i,j+1))/2*dx;
				}
			acf();
			as[i][j+1][1]=d*adf+max(flow,0);
			an[i][j][1]=as[i][j+1][1]-flow;	
		}//for j
	}//for i
		//边界条件
	for(j=2+inlen1;j<=m-1;j++)
		aw[3][j][1]=0;//左壁aw＝0
	for(j=2+inlen2;j<=m-1;j++)
		ae[l-1][j][1]=0;
	for(i=3;i<=l-1;i++)
		as[i][2][1]=an[i][m-1][1]=0;

	//ap和b设定
	for(i=3;i<=l-1;i++)
		for(j=2;j<=m-1;j++)
		{
			if(i==3||i==l-1) ap0[i][j][1]=den*dx*3/2*dy/dt;
			else ap0[i][j][1]=den*dx*dy/dt;
			ap[i][j][1]=ae[i][j][1]+aw[i][j][1]+an[i][j][1]+as[i][j][1]+ap0[i][j][1];
			ap[i][j][1]=ap[i][j][1]/relax;
			b[i][j][1]=ap0[i][j][1]*u0(i,j)+(p[i-1][j]-p[i][j])*dy+(1-relax)*ap[i][j][1]*u0(i,j);
		
		}

/*			i=5;
	printf("\nv方程系数\n");
		for(j=3;j<=m-1;j++)
		{printf("%d %f,%f,%f,%f,%f,%f,%f,\n",j,aw[i][j][1],ae[i][j][1],as[i][j][1],an[i][j][1],ap0[i][j][1],ap[i][j][1],b[i][j][1]);
	fprintf(filehandle,"%d %f,%f,%f,%f,%f,%f,%f,\n",j,aw[i][j][1],ae[i][j][1],as[i][j][1],an[i][j][1],ap0[i][j][1],ap[i][j][1],b[i][j][1]);
		}

*/


	//v方程系数设定

	//as和an设定
	for(i=2;i<=l-1;i++)
	{//下边界as[i][3]
		d=dif[n]*dx/dy;
		flow=den*v(i,2)*dx;
		acf();
		as[i][3][2]=d*adf+max(flow,0);
		//上边界an[i][m-1]
		flow=den*v(i,m)*dx;
		acf();
		an[i][m-1][2]=d*adf+max(-1*flow,0);
		for(j=3;j<=m-2;j++)
		{
		//v[i][j]的an和v[i][j+1]的as
			flow=den*dx*(v(i,j)+v(i,j+1))/2;
			acf();
			as[i][j+1][2]=d*adf+max(flow,0);
			an[i][j][2]=as[i][j+1][2]-flow;
		}
	}

	//ae和aw设定
	for(j=3;j<=m-1;j++)
	{		
		//左边界aw和右边界ae
		if(j==3)
		{	//aw[2][3]
			d=dif[n]*dy*3/dx;
			flow=den*dy*(u(2,j-1)+u(2,j)/2);
			acf();
			aw[2][j][2]=d*adf+max(flow,0);
			//ae[l-1][3]
			flow=den*dx*(u(l,j-1)+u(l,j)/2);
			acf();
			ae[l-1][j][2]=d*adf+max(-1*flow,0);	
		}
			else  
				if(j==m-1)
				{	//aw[2][m-1]
					d=dif[n]*dy*3/dx;
					flow=den*dy*(u(2,j-1)/2+u(2,j));
					acf();
					aw[2][j][2]=d*adf+max(flow,0);
					//ae[l-1][m-1]
					flow=den*dy*(u(l,j-1)/2+u(l,j));
					acf();
					ae[l-1][j][2]=d*adf+max(-1*flow,0);
				}
			
				else
				{
					d=dif[n]*dy*2/dx;
					//aw[2][j]
					flow=den*dy*(u(2,j-1)+u(2,j))/2;
					acf();
					aw[2][j][2]=d*adf+max(flow,0);
					//ae[l-1][j]
					flow=den*dy*(u(l,j-1)+u(l,j))/2;
					acf();
					ae[l-1][j][2]=d*adf+max(-1*flow,0);
					}
		//内部边界ae[i][j]和aw[i+1][j]
		for(i=2;i<=l-2;i++)
		{
			if(j==3)
			{
				d=dif[n]*dy*3/2/dx;
				flow=den*dy*(u(i+1,j-1)+u(i+1,j)/2);
			}
			else
				if(j==m-1)
				{
					d=dif[n]*dy*3/2/dx;
					flow=den*dy*(u(i+1,j-1)/2+u(i+1,j));
				}
				else
				{
					d=dif[n]*dy/dx;
					flow=den*dy*(u(i+1,j-1)+u(i+1,j))/2;
				}
			acf();
			aw[i+1][j][2]=d*adf+max(flow,0);
			ae[i][j][2]=aw[i+1][j][2]-flow;				
		}//for i
	}//for j
	//ap0和b设定
	for(j=3;j<=m-1;j++)
		for(i=2;i<=l-1;i++)
		{
			if(j==3||j==m-1) ap0[i][j][2]=den*dy*1.5*dx/dt;
			else ap0[i][j][2]=den*dy*dx/dt;
			ap[i][j][2]=an[i][j][2]+as[i][j][2]+ae[i][j][2]+aw[i][j][2]+ap0[i][j][2];
			ap[i][j][2]=ap[i][j][2]/relax;
			b[i][j][2]=ap0[i][j][2]*v0(i,j)+(p[i][j-1]-p[i][j])*dx+(1-relax)*ap[i][j][2]*v0(i,j);
		}
  //边界条件
	for(j=2+inlen1;j<=m-1;j++)
		aw[2][j][2]=0;
	for(j=2+inlen2;j<=m-1;j++)
		ae[l-1][j][2]=0;
	for(i=2;i<=l-1;i++)
		as[i][3][2]=an[i][m-1][2]=0;
	
/*		i=2;
	printf("\nv方程系数\n");
		for(j=3;j<=m-1;j++)
		{printf("%d %f,%f,%f,%f,%f,%f,%f,\n",j,aw[i][j][2],ae[i][j][2],as[i][j][2],an[i][j][2],ap0[i][j][2],ap[i][j][2],b[i][j][2]);
	fprintf(filehandle,"%d %f,%f,%f,%f,%f,%f,%f,\n",j,aw[i][j][2],ae[i][j][2],as[i][j][2],an[i][j][2],ap0[i][j][2],ap[i][j][2],b[i][j][2]);
		}
*/	


	//求u方程
	n=1;
	ist=3;
	jst=2;
	solve();
	n=0;

	//求v方程
	n=2;
	ist=2;
	jst=3;
	solve();
	n=0;

	//压力校正方程设定
	for(i=2;i<=l-1;i++)
	{	
		for(j=2;j<=m-1;j++)
		{
			ae[i][j][3]=den*dy*dy/ap[i+1][j][1];
			aw[i][j][3]=den*dy*dy/ap[i][j][1];
			an[i][j][3]=den*dx*dx/ap[i][j+1][2];
			as[i][j][3]=den*dx*dx/ap[i][j][2];
		}

	}
		//边界修正值
	for(i=2;i<=l-1;i++)
		as[i][2][3]=an[i][m-1][3]=0;
	for(j=2;j<=m-1;j++)
		aw[2][j][3]=ae[l-1][j][3]=0;

	remain=b[2][2][3];
	for(i=2;i<=l-1;i++)
		for(j=2;j<=m-1;j++)
		{	
			ap[i][j][3]=ae[i][j][3]+aw[i][j][3]+an[i][j][3]+as[i][j][3];
			b[i][j][3]=den*((u(i,j)-u(i+1,j))*dy+((v(i,j)-v(i,j+1))*dx));
			if(remain<fabs(b[i][j][3]))remain=b[i][j][3];
		}	
	
	
	//求压力校正方程
	n=3;
	ist=2;
	jst=2;
	solve();
	n=0;


	//压力和速度校正
	for(i=2;i<=l-1;i++)
		for(j=2;j<=m-1;j++)
		{
			p[i][j]=p[i][j]+dp(i,j)*relaxp;
			if(i!=2)u(i,j)=u(i,j)+(dp(i-1,j)-dp(i,j))*dy/ap[i][j][1];
			if(j!=2)v(i,j)=v(i,j)+(dp(i,j-1)-dp(i,j))*dx/ap[i][j][2];
		}
	//压力参考点
		minp=p[2][2];
	for(i=2;i<=l-1;i++)
		for(j=2;j<=m-1;j++)
			{
				p[i][j]=p[i][j]-minp;
			}
/*		printf("\np方程校正值\n");
		for(j=m;j>=1;j--)
		{printf("\n");
			for(i=1;i<=l;i++)
				printf("%f,",f[i][j][3]);
		}
*/
}