[arotester]

Hello,

I have written a code for flow-through a channel(2D) using a pressure correction approach(SIMPLE) on collocated grid. 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(it's not parabolic) and also the center velocity is lower than the inlet.

Regards,

//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[200],y[200],u[200][200],uc[200][200],unew[200][200];
double v[200][200],vc[200][200],vnew[200][200],p[200][200],pc[200][200],pnew[200][200];
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[200][200],vp[200][200];
double upe,upw,vpn,vps,s[200][200],u0,corr,Res,uold;

main()
{
L=10.0;
H=1.0;
nx=100;
ny=10;
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]);
    }
//printf("code is here");
//________________________________________________Initialization_____________________________________________________
for(i=1;i<=nx;i++)
    {
    for(j=1;j<=ny;j++)
        {
        u[i][j]=0.0;
        v[i][j]=0.0;
      // 	p[i][j]=0.0;
       	pc[i][j]=0.0;
        corr=0.0;
        Res=0.0;
        }
    }

u0=0.01;  
//uold=0.0;      
//__________________________________________________________________________________________________________________    
//________________________________________________Solution Loops____________________________________________________|    

for(count=1;count<=50000;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];
    //    p[i][0]=p[i][1];
        u[i][ny+1]=-u[i][ny];
        v[i][ny+1]=-v[i][ny];
    //    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];
    //	p[0][j]=p[1][j];
        u[nx+1][j]=u[nx][j];
        v[nx+1][j]=v[nx][j];
	//	p[nx+1][j]=-p[nx][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];
                            
                        }else  
                            {
                                uw=u[i][j];
                                x_velw=u[i][j];
                                y_velw=v[i][j];
                            }
                    if(ue>0)
                        {
                            ue=u[i][j];
                            x_vele=u[i][j];
                            y_vele=v[i][j];
                        }else 
                            {
                                ue=u[i+1][j];
                                x_vele=u[i+1][j];    
                                y_vele=v[i+1][j];
                            }
                    if(vs>0)
                        {
                            vs=v[i][j-1];
                            y_vels=v[i][j-1];
                            x_vels=u[i][j-1];
                        }else 
                            {
                                  vs=v[i][j];
                                x_vels=u[i][j];
                                y_vels=v[i][j];
                            }
                    if(vn>0)
                        {
                            vn=v[i][j];
                            y_veln=v[i][j];
                            x_veln=u[i][j]; 
                        }else 
                            {
                                vn=v[i][j+1];
                                x_veln=u[i][j+1];
                                y_veln=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]);
            px=(1.0/rho)*((pe-pw))*(dy);
            up[i][j]=u[i][j]+dt/(dx*dy)*(u_diff-u_adv-px);
    	  
                        
		    //------------------------------y-momentum---------------------------------//
            v_adv=-uw*y_velw*dy+ue*y_vele*dy+vn*y_veln*dx-vs*y_vels*dx;         
            v_diff=alp*(-ap*v[i][j]+ae*v[i+1][j]+aw*v[i-1][j]+an*v[i][j+1]+as*v[i][j-1]);
            py=(1.0/rho)*((pn-ps))*(dx);
            vp[i][j]=v[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=(up[i+1][j]+up[i][j])/2.0;
            upw=(up[i-1][j]+up[i][j])/2.0;
            vpn=(vp[i][j+1]+vp[i][j])/2.0;
            vps=(vp[i][j-1]+vp[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.0001);
        printf("pc=%f\tu=%f\tc=%d\n",pc[nx/2][ny/2],u[nx/2][ny/2],count);
//----------------------------------Velocity Corrections----------------------------------------
for(i=1;i<=nx;i++)
    {
    for(j=1;j<=ny;j++)
        {    
		//	u[i][j]=up[i][j]-dt/(2.0*dx*rho)*(pc[i+1][j]-pc[i-1][j]);
		//    v[i][j]=vp[i][j]-dt/(2.0*dy*rho)*(pc[i][j+1]-pc[i][j-1]);
            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];
            vnew[i][j]=vp[i][j]+vc[i][j];
            unew[i][j]=up[i][j]+uc[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;
        }
    }    
//---------------------------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];
           u[i][j]=unew[i][j];
            v[i][j]=vnew[i][j];
        
        }
    }    

}
//--------------------------File Writing----------------------------------------------------------------------------
FILE *f1;
f1=fopen("Temp.dat","w");
    
    for(j=1;j<=ny;j++)
        {    
            fprintf(f1,"%f\t%f\n",u[nx/2][j],y[j]);
        }
fclose(f1);        
}

[arotester]

Hello,
Did anyone get the mistake because I am still trying to debug it and it's not working?