Lid Driven Cavity using Ghost Cell Method and in C++
 

Hi,
I am programming the lid driven cavity using the ghost cell method and using the 9 point laplacian to solve the pressure equation. I have programmed just the laplacian before coding the driven cavity and the laplacian part is fine. However the entire program seems to be diverging for every iteration instead of converging. And also only the y-velocity seems to be updating for each iteration. I am using the QUICK scheme to calculate the fluxes. Could someone please have a look at the program and tell me where I am going wrong.

http://www.menet.umn.edu/~kaustubh/R...UMNME5351.html

I am using the above site as a reference for my results. L=0.4m, H=0.3m, Re=600

#include<iostream>
#include"fstream"
#include"conio.h"
#include"stdlib.h"
#include"math.h"
#include"iomanip"

using namespace std;

double Re=600.0;
double nu;
int grid_size=16;
int counter,i,j;
double dx,dy,dt,old,residual_L2,residual,g,outer_residual ,u_residual,x,y;
double ux,vy,Fxx,Gyy,Hx_xy,Hy_xy;
double Lx=0.4;
double Ly=0.3;
double U=5.0;
double w=1.0;
double factor=pow(10.0,-10.0);
double factor2=pow(10.0,-3.0);
double factor3=pow(10.0,-2.0);
double diffusion_time,advection_time;
double time=0.0;
double u[1000][1000],v[1000][1000],p[1000][1000],F[1000][1000],G[1000][1000],Hx[1000][1000],Hy[1000][1000];
double ff,gf,hx,hy;
double old_p[1000][1000],old_u[1000][1000];
double p_ref,q,qh,phi;
double cell_centered_u[1000][1000],cell_centered_v[1000][1000],nodal_p[1000][1000];
double stream[1000][1000];

void input();
void ic_u();
void ic_v();
void ghost_u();
void ghost_v();
void F_QUICK();
void G_QUICK();
void Hx_QUICK();
void Hy_QUICK();
void time_step();
void update_uv();
void pressure_ic();
void pressure_bc();
void interior_pressure_poisson();
void pressure_correction();
void initial();
void calculation();
void cell_centered_velocities();
void stream_fucntion();
void output();

ofstream fout ("OUTPUT.dat");
ofstream fout2 ("flux.dat");

int main()
{
initial();
calculation();
output();
}

void input()
{
//cout<<"ENTER REYNOLDS NUMBER: ";
//cin>>Re;
nu=(U*Lx)/Re;
//cout<<"ENTER GRID SIZE: ";
//cin>>grid_size;
dx=Lx/grid_size;
dy=Ly/grid_size;
for(i=0;i<grid_size+2;i++)
for(j=0;j<grid_size+2;j++)
{
cell_centered_u[i][j]=0.0;
cell_centered_v[i][j]=0.0;
stream[i][j]=0.0;
nodal_p[i][j]=0.0;
}
}

void ic_u()
{
for (i=0;i<grid_size+2;i++)
{
for (j=0;j<grid_size+3;j++)
{
u[i][j]=0.0;
}
}
}

void ghost_u()
{
for(j=0;j<grid_size+3;j++)
{
u[0][j]=(2.0*U)-u[1][j];
u[grid_size+1][j]=-u[grid_size][j];
}

for(i=0;i<grid_size+2;i++)
{
u[i][0]=u[i][2];
u[i][grid_size+2]=u[i][grid_size];
}
}

void ic_v()
{
for (i=0;i<grid_size+3;i++)
{
for (j=0;j<grid_size+2;j++)
{
v[i][j]=0.0;
}
}
}

void ghost_v()
{
for(i=0;i<grid_size+3;i++)
{
v[i][0]=-v[i][1];
v[i][grid_size+1]=-v[i][grid_size];
}

for(j=0;j<grid_size+2;j++)
{
v[0][j]=v[2][j];
v[grid_size+2][j]=v[grid_size][j];
}

}

void pressure_ic()
{
for(i=0;i<grid_size+2;i++)
for(j=0;j<grid_size+2;j++)
{
p[i][j]=0.0;
}
}

void pressure_bc()
{
for(j=0;j<grid_size+2;j++)
{
p[0][j]=p[1][j]-(2.0*nu*v[2][j]/dy);
p[grid_size+1][j]=p[grid_size][j]+(2.0*nu*v[grid_size][j]/dy);
}

for(i=0;i<grid_size+2;i++)
{
p[i][0]=p[i][1]-(2.0*nu*u[i][2]/dx);
p[i][grid_size+1]=p[i][grid_size]+(2.0*nu*u[i][grid_size]/dx);
}
}

void pressure_correction()
{
p_ref=p[grid_size][1];
for(i=1;i<grid_size+1;i++)
{
for(j=1;j<grid_size+1;j++)
p[i][j]=p[i][j]-p_ref;
}
}

void interior_pressure_poisson()
{
for(counter=1;counter<=1000000;counter++)
{
residual_L2=0.0;
for(i=grid_size;i>=1;i--)
{
for(j=1;j<grid_size+1;j++)
{
old=p[i][j];
ux=(u[i][j+1]-u[i][j])/(dx);
vy=(v[i][j]-v[i+1][j])/(dy);
Fxx=(F[i][j-1]-(2.0*F[i][j])+F[i][j+1])/pow(dx,2.0);
Hx_xy=(((Hx[i][j+1]-Hx[i][j]))-((Hx[i+1][j+1]-Hx[i+1][j])))/(dy*dx);
Hy_xy=(((Hy[i][j+1]-Hy[i][j]))-((Hy[i+1][j+1]-Hy[i+1][j])))/(dy*dx);
Gyy=(G[i-1][j]-(2.0*G[i][j])+G[i+1][j])/pow(dy,2.0);
g=(-(dx/dt)*(ux+vy))-(Fxx+Hx_xy+Hy_xy+Gyy);
residual=w*(((1.0/5.0)*(p[i][j-1]+p[i][j+1]+p[i-1][j]+p[i+1][j]))+((1.0/20.0)*(p[i-1][j-1]+p[i+1][j-1]+p[i-1][j+1]+p[i+1][j+1]))-(old)+((6.0/20.0)*g*dx*dy));
p[i][j]=(old+residual);
residual_L2+=pow(residual,2.0);
}
}
residual_L2=pow((1.0/(grid_size*grid_size))*residual_L2,0.5);
if(residual_L2<=factor)
{
break;
}
}
}

void F_QUICK()
{
for(i=grid_size;i>=1;i--)
{
for(j=1;j<grid_size+1;j++)
{
q=(u[i][j]+u[i][j+1])/2.0;
if(q<0)
{
phi=(1.0/8.0)*((3.0*u[i][j])+(6.0*u[i][j+1])-(u[i][j+2]));
ff=q*phi;
//ff=pow(u[i][j+1],2.0);
}
else if(q>0)
{
phi=(1.0/8.0)*((3.0*u[i][j+1])+(6.0*u[i][j])-(u[i][j-1]));
ff=q*phi;
//ff=pow(u[i][j],2.0);
}
else
{
ff=0.0;
}

qh=(u[i][j+1]-u[i][j])/dx;
F[i][j]=(ff)-(nu*qh);
}
}
}

void G_QUICK()
{
for(i=grid_size;i>=1;i--)
{
for(j=1;j<grid_size+1;j++)
{
q=(v[i][j]+v[i+1][j])/2.0;
if(q<0)
{
phi=(1.0/8.0)*((3.0*v[i+1][j])+(6.0*v[i][j])-(v[i-1][j]));
gf=q*phi;
//gf=pow(v[i][j],2.0);
}
else if(q>0)
{
phi=(1.0/8.0)*((3.0*v[i][j])+(6.0*v[i+1][j])-(v[i+2][j]));
gf=q*phi;
//gf=pow(v[i+1][j],2.0);
}
else
{
gf=0.0;
}

qh=(v[i][j]-v[i+1][j])/dy;
G[i][j]=(gf)-(nu*qh);
}
}
}

void Hx_QUICK()
{
for(i=grid_size;i>=2;i--)
{
for(j=2;j<grid_size+1;j++)
{
q=(v[i][j]+v[i][j-1])/2.0;
if(q<0)
{
phi=(1.0/8.0)*((3.0*u[i][j])+(6.0*u[i-1][j])-(u[i-2][j]));
hx=q*phi;
//hx=pow(u[i-1][j],2.0);
}
else if(q>0)
{
phi=(1.0/8.0)*((3.0*u[i-1][j])+(6.0*u[i][j])-(u[i+1][j]));
hx=q*phi;
//hx=pow(u[i][j],2.0);
}
else
{
hx=0.0;
}

qh=(v[i][j]-v[i][j-1])/dx;
Hx[i][j]=(hx)-(nu*qh);
}
}
}

void Hy_QUICK()
{
for(i=grid_size;i>=2;i--)
{
for(j=2;j<grid_size+1;j++)
{
q=(u[i][j]+u[i-1][j])/2.0;
if(q<0)
{
phi=(1.0/8.0)*((3.0*v[i][j-1])+(6.0*v[i][j])-(v[i][j+1]));
hy=q*phi;
//hy=pow(v[i][j],2.0);
}
else if(q>0)
{
phi=(1.0/8.0)*((3.0*v[i][j])+(6.0*v[i][j-1])-(v[i][j-2]));
hy=q*phi;
//hy=pow(v[i][j-1],2.0);
}
else
{
hy=0.0;
}

qh=(u[i-1][j]-u[i][j])/dy;
Hy[i][j]=(hy)-(nu*qh);
}
}
}

void time_step()
{
diffusion_time=dx*dy/(2.0*nu);
advection_time=dx/(abs(U));
if(advection_time<diffusion_time)
dt=advection_time;
else
dt=diffusion_time;
}

void update_uv()
{

for(i=grid_size;i>=1;i--)
{
for(j=2;j<grid_size+1;j++)
{
u[i][j]=(u[i][j]-((dt/dx)*((F[i][j]-F[i][j-1])+(p[i][j]-p[i][j-1])+(Hx[i][j]-Hx[i+1][j]))));
}
}

for(i=grid_size;i>=2;i--)
{
for(j=1;j<grid_size+1;j++)
{
v[i][j]=(v[i][j]-((dt/dy)*((G[i-1][j]-G[i][j])+(p[i-1][j]-p[i][j])+(Hy[i][j+1]-Hy[i][j]))));
}
}

ghost_u();
ghost_v();
F_QUICK();
G_QUICK();
Hx_QUICK();
Hy_QUICK();
}

void cell_centered_velocities()
{
for(i=1;i<grid_size+2;i++)
{
for(j=1;j<grid_size+2;j++)
{
cell_centered_u[i][j]=(u[i-1][j]+u[i][j])/2.0;
cell_centered_v[i][j]=(v[i][j-1]+v[i][j])/2.0;
}
}
}

void nodal_pressure()
{
for(i=1;i<grid_size+2;i++)
{
for(j=1;j<grid_size+2;j++)
{
nodal_p[i][j]=(p[i-1][j-1]+p[i-1][j]+p[i][j-1]+p[i][j])/4.0;
}
}
}

void stream_function()
{

for(i=grid_size;i>=1;i--)
{
for(j=1;j<grid_size+2;j++)
{
stream[i][j]=stream[i+1][j]+(dy*cell_centered_u[i][j]);
}
}

for(i=grid_size+1;i>=1;i--)
{
for(j=2;j<grid_size+2;j++)
{
stream[i][j]=stream[i][j-1]-(dx*cell_centered_v[i][j]);
}
}
}

void initial()
{
input();
time_step();
pressure_ic();
ic_u();
ghost_u();
ic_v();
ghost_v();
F_QUICK();
G_QUICK();
Hx_QUICK();
Hy_QUICK();
}

void calculation()
{
system("CLS");
cout<<"SOLVING PRESSURE POISSON EQUATION"<<endl<<endl;
cout<<"PLEASE WAIT ";
outer_residual=1.0;
while(time<(70.0*dt))
{
pressure_bc();
interior_pressure_poisson();
pressure_correction();
update_uv();
time+=dt;
cout<<".";
}
cout<<endl;
}

void output()
{
cell_centered_velocities();
nodal_pressure();
stream_function();

fout<<"\"Variables= x\""<<"\t"<<"\"y\""<<"\t"<<"\"u (m/s)\""<<"\t"<<"\"v (m/s)\""<<"\t"<<"\"Pressure (atm)\""<<"\t"<<"\"Stream Function\""<<endl;
fout<<"zone i="<<grid_size+1<<" j="<<grid_size+1<<" f=point"<<endl;

y=Ly;
for(i=1;i<grid_size+2;i++)
{
x=0;
for(j=1;j<grid_size+2;j++)
{
fout<<x<<"\t"<<y<<"\t"<<cell_centered_u[i][j]<<"\t"<<cell_centered_v[i][j]<<"\t"<<nodal_p[i][j]<<"\t"<<stream[i][j]<<endl;
fout2<<i<<"\t"<<j<<"\t"<<F[i][j]<<"\t"<<G[i][j]<<endl;
x=x+dx;
}
y=y-dy;
}
}


Thanks