[arotester]

Hello,

I am trying to write a code for flow-through a channel(2D) using a pressure correction approach(SIMPLE). But, I am not able to get the right solution, can you tell me what is the mistake with the code?(code attached)
The code gives the unknown spikes in the velocity profile which for me is difficult to resolve.

Regards,

[arotester]

I have done some debugging but now the problem is different. Velocity profile becomes flatter at the end of the pipe.

//Primitive variable formulation with colocated grid

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

int i,j,nx,ny,count;
double L,H,dx,dy,dt,beta,alp,rho,x[500],y[500],u[500][500],uc[500][500],unew[500][500];
double v[500][500],vc[500][500],vnew[500][500],p[500][500],pc[500][500],pnew[500][500];
double as,ae,an,aw,ap,sum,sum1,sum2,diff,diff1,sqr,sqr1,sqr2,RMS,RMS1,RMS2;
double uw,ue,vn,vs,x_vel,y_vel,pe,pn,ps,pw,px,py,u_adv,u_diff,v_adv,v_diff,up[500][500],vp[500][500];
double upe,upw,vpn,vps,s[500][500],u0,corr,Res,uold;
double upnew[500][500],vpnew[500][500],ups,upn,vpw,vpe,uup[500][500],vup[500][500];
main()
{
L=10.0;
H=1.0;
nx=200;
ny=30;
dx=L/nx;
dy=H/ny;
dt=0.01;
alp=0.00001;//kinematic viscosity
rho=1.2;
beta=1.0;//relaxation factor
aw=dy/dx;
ae=dy/dx;
an=dx/dy;
as=dx/dy;
ap=aw+as+ae+aw;

//_____________________________________________________Grid__________________________________________________________
x[1]=0.5*dx;
for(i=2;i<=nx;i++)
	{
		x[i]=x[i-1]+dx;
	//	printf("x=%f\n",x[i]);
	}
y[1]=0.5*dy;
for(j=2;j<=ny;j++)
	{
		y[j]=y[j-1]+dy;
	//	printf("y=%f\n",y[j]);
	}
//________________________________________________Initialization_____________________________________________________
for(i=1;i<=nx;i++)
	{
	for(j=1;j<=ny;j++)
		{
		u[i][j]=0.0;
		v[i][j]=0.0;
		up[i][j]=0.0;
		vp[i][j]=0.0;
		p[i][j]=0.0;
		pc[i][j]=0.0;
		corr=0.0;
		Res=0.0;
		}
	}
u0=0.01;		
//__________________________________________________________________________________________________________________	
//________________________________________________Solution Loops____________________________________________________|	

for(count=1;count<=75000;count++)
//do
{

//_______________________________________Solving Momentum Equations__________________________________________________

//_______________________________Boundary Conditions_________________________________________________________________
for(i=1;i<=nx;i++)
	{	
		u[i][0]=-u[i][1];
		v[i][0]=-v[i][1];
		
		u[i][ny+1]=-u[i][ny];
		v[i][ny+1]=-v[i][ny];
	
		up[i][0]=-up[i][1];
		vp[i][0]=-vp[i][1];
		
		up[i][ny+1]=-up[i][ny];
		vp[i][ny+1]=-vp[i][ny];
		
	//	p[i][0]=p[i][1];
	//	p[i][ny+1]=p[i][ny];
		
	}
for(j=1;j<=ny;j++)
	{
		u[0][j]=2.0*u0-u[1][j];
		v[0][j]=-v[1][j];
		
		u[nx+1][j]=u[nx][j];
		v[nx+1][j]=v[nx][j];
	
		up[0][j]=2.0*u0-up[1][j];
		vp[0][j]=-vp[1][j];
		
		up[nx+1][j]=up[nx][j];
		vp[nx+1][j]=vp[nx][j];
		
	//	p[nx+1][j]=-p[nx][j];
	//	p[0][j]=p[1][j];
	}

	for(i=1;i<=nx;i++)
		{
		for(j=1;j<=ny;j++)
			{		
			double x_velw,x_vele,x_vels,x_veln,y_velw,y_vele,y_veln,y_vels;
				
					uw=(u[i][j]+u[i-1][j])/2.0;
					ue=(u[i][j]+u[i+1][j])/2.0;
					vs=(v[i][j]+v[i][j-1])/2.0;
					vn=(v[i][j]+v[i][j+1])/2.0;
					
					pe=(p[i][j]+p[i+1][j])/2.0;
					pw=(p[i-1][j]+p[i][j])/2.0;
					pn=(p[i][j+1]+p[i][j])/2.0;
					ps=(p[i][j-1]+p[i][j])/2.0;	
			//---------------FOU Routine------------------------------//
					if(uw>0)
						{
							uw=u[i-1][j];
							x_velw=u[i-1][j];
							y_velw=v[i-1][j];
							upw=u[i-1][j];
							vpw=v[i-1][j];
						}else  
							{
								uw=u[i][j];
								x_velw=u[i][j];
								y_velw=v[i][j];
								upw=u[i][j];
								vpw=v[i][j];
							}
					if(ue>0)
						{
							ue=u[i][j];
							x_vele=u[i][j];
							y_vele=v[i][j];
							upe=u[i][j];
							vpe=v[i][j];
						}else 
							{
								ue=u[i+1][j];
								x_vele=u[i+1][j];	
								y_vele=v[i+1][j];
								upe=u[i+1][j];
								vpe=v[i+1][j];
							}
					if(vs>0)
						{
							vs=v[i][j-1];
							y_vels=v[i][j-1];
							x_vels=u[i][j-1];
							ups=u[i][j-1];
							vps=v[i][j-1];
						}else 
							{
							  	vs=v[i][j];
								x_vels=u[i][j];
								y_vels=v[i][j];
								ups=u[i][j];
								vps=v[i][j];
							}
					if(vn>0)
						{
							vn=v[i][j];
							y_veln=v[i][j];
							x_veln=u[i][j]; 
							upn=u[i][j];
							vpn=v[i][j];
						}else 
							{
								vn=v[i][j+1];
								x_veln=u[i][j+1];
								y_veln=v[i][j+1];
								upn=u[i][j+1];
								vpn=v[i][j+1];
							}
							
					 	
			//---------------------------x-momentum------------------------------------//
			//u_adv=-uw*x_velw*dy+ue*x_vele*dy+vn*x_veln*dx-vs*x_vels*dx; 		
			//u_diff=alp*(-ap*u[i][j]+ae*u[i+1][j]+aw*u[i-1][j]+an*u[i][j+1]+as*u[i][j-1]);
			u_adv=-uw*upw*dy+ue*upe*dy+vn*upn*dx-vs*ups*dx; 		
			u_diff=alp*(-ap*up[i][j]+ae*up[i+1][j]+aw*up[i-1][j]+an*up[i][j+1]+as*up[i][j-1]);
			px=(1.0/rho)*((pe-pw))*(dy);
			upnew[i][j]=up[i][j]+dt/(dx*dy)*(u_diff-u_adv-px);
			//------------------------------y-momentum---------------------------------//
			v_adv=-uw*vpw*dy+ue*vpe*dy+vn*vpn*dx-vs*vps*dx; 		
			v_diff=alp*(-ap*vp[i][j]+ae*vp[i+1][j]+aw*vp[i-1][j]+an*vp[i][j+1]+as*vp[i][j-1]);
			py=(1.0/rho)*((pn-ps))*(dx);
			vpnew[i][j]=vp[i][j]+dt/(dx*dy)*(v_diff-v_adv-py);
			}
		}

			
//___________________________________Mass source term(check)____________________________________________________________
for(i=1;i<=nx;i++)
	{
	for(j=1;j<=ny;j++)
		{	
			upe=(upnew[i+1][j]+upnew[i][j])/2.0;
			upw=(upnew[i-1][j]+upnew[i][j])/2.0;
			vpn=(vpnew[i][j+1]+vpnew[i][j])/2.0;
			vps=(vpnew[i][j-1]+vpnew[i][j])/2.0;
			s[i][j]=(upe*dy-upw*dy+vpn*dx-vps*dx);
			
		}
	}	

//__________________________________Pressure correction__________________________________________________________
//GSSOR

do
	{
sum2=0.0;		
//-------------------------Boundary conditions for pressure correction-------------------------------
for(i=1;i<=nx;i++)
	{
		pc[i][0]=pc[i][1];
		pc[i][ny+1]=pc[i][ny];

	}
for(j=1;j<=ny;j++)
	{
		pc[0][j]=pc[1][j];
		pc[nx+1][j]=-pc[nx][j];
	
	}
//_________________________________________________________________________________________________________________
		for(i=1;i<=nx;i++)
			{
			for(j=1;j<=ny;j++)
				{
					double p_corr,source;
					p_corr=(-ap*pc[i][j]+aw*pc[i-1][j]+as*pc[i][j-1]+an*pc[i][j+1]+ae*pc[i+1][j])*dt/rho;
					source=s[i][j];
					Res=(source-p_corr);
					corr=Res/(-beta*ap);
					pc[i][j]=pc[i][j]+corr;
					sqr2=Res*Res;
					sum2=sum2+sqr2;
					RMS2=sqrt(sum2/((nx)*(ny)));
					
				}
			}
		}while(RMS2>0.00001);
		printf("pc=%f\tu=%f\tc=%d\n",pnew[nx/2][ny/2],uup[nx/2][ny/2],count);
//----------------------------------Velocity Corrections----------------------------------------
for(i=1;i<=nx;i++)
	{
	for(j=1;j<=ny;j++)
		{		
			uc[i][j]=-dt/(2.0*rho*dx)*(pc[i+1][j]-pc[i][j]);
			vc[i][j]=-dt/(2.0*rho*dy)*(pc[i][j+1]-pc[i][j]);	
		}
	}
//--------------------------------Find new values from correction----------------------------------------------------
for(i=1;i<=nx;i++)
	{
	for(j=1;j<=ny;j++)
		{		
			
			pnew[i][j]=p[i][j]+pc[i][j];
			vup[i][j]=vpnew[i][j]+vc[i][j];
			uup[i][j]=upnew[i][j]+uc[i][j];
		}
	}	
/*	
sum=0.0;	
//RMS Calculations	
for(i=1;i<=nx;i++)
	{
	for(j=1;j<=ny;j++)
		{	
			diff=up[i][j]-upnew[i][j];
			sqr=diff*diff;
			sum=sum+sqr;
			RMS=sqrt(sum/(nx*ny));
		}
	}
*/
//---------------------------Updating the values of pressure and velocity----------------------------------------- 		
for(i=1;i<=nx;i++)
	{
	for(j=1;j<=ny;j++)
		{	
			p[i][j]=pnew[i][j];
			up[i][j]=uup[i][j];
			vp[i][j]=vup[i][j];
		
		}
	}	

//---------------------------Making correction zero at the end of step--------------------------------------		
for(i=1;i<=nx;i++)
	{
	for(j=1;j<=ny;j++)
		{		
			pc[i][j]=0.0;
			uc[i][j]=0.0;
			vc[i][j]=0.0;
		}
	}	

	//printf("p=%f\tu=%f\tRMS=%f\n",p[nx/2][ny/2],uup[nx/2][ny/2],RMS);				
//}while(RMS>0.000001);
}
//--------------------------File Writing----------------------------------------------------------------------------
FILE *f1;
f1=fopen("Temp_2.0_1.dat","w");
	
	for(j=1;j<=ny;j++)
		{	
			fprintf(f1,"%f\t%f\n",uup[50][j],y[j]);
		}
fclose(f1);		
FILE *f2;
f2=fopen("Temp_2.0_2.dat","w");
	
	for(j=1;j<=ny;j++)
		{	
			fprintf(f2,"%f\t%f\n",uup[100][j],y[j]);
		}
fclose(f2);		
FILE *f3;
f3=fopen("Temp_2.0_3.dat","w");
	
	for(j=1;j<=ny;j++)
		{	
			fprintf(f3,"%f\t%f\n",uup[150][j],y[j]);
		}
fclose(f3);		
FILE *f4;
f4=fopen("Temp_2.0_4.dat","w");
	
	for(j=1;j<=ny;j++)
		{	
			fprintf(f4,"%f\t%f\n",uup[195][j],y[j]);
		}
fclose(f4);		



}

[FMDenaro]

Originally Posted by arotester

I have done some debugging but now the problem is different. Velocity profile becomes flatter at the end of the pipe.

I see wiggles in the inlet profile and then the continuity is not fulfilled.
Check the correct BC.

[arotester]

Thank you for your reply, I appreciate it.
I debugged the code and extended the code to 3D and this was the result.....

[randolph]

It looks like there might be some interpolation issue. Your BC looks right.