CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

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

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   August 12, 2011, 23:05
Default Lid Driven Cavity using Ghost Cell Method and in C++
  #1
New Member
 
Vignesh Suresh
Join Date: Jul 2011
Posts: 10
Rep Power: 5
illuminati5288 is on a distinguished road
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
illuminati5288 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On



All times are GMT -4. The time now is 09:20.