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

Boundary condition for 2d lid driven cavity using ghost cells

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

Reply
 
LinkBack Thread Tools Display Modes
Old   April 7, 2009, 03:12
Default Boundary condition for 2d lid driven cavity using ghost cells
  #1
Senior Member
 
TWB
Join Date: Mar 2009
Posts: 105
Rep Power: 8
quarkz is on a distinguished road
Hi,

I'm trying to use ghost cells to implement boundary condition for the 2d lid driven cavity problem. I'm using staggered grid and I wonder if my implementation is correct because the ans diverged after 20+ time steps.

My BC are:

u(1,:)=-u(-1,:) where u/v(-1,:) is the ghost cell for left wall
v(1,:)=-v(-1,:)
p(1,:)=-p(-1,:)

u(size_x,:)=-u(size_x+1,:) where u/v(size_x+1,:) is the ghost cell for right wall
v(size_x,:)=-v(size_x+1,:)
p(size_x,:)=-p(size_x+1,:)

u(:,1)=-u(:,-1) where u/v(:,-1) is the ghost cell for bottom wall
v(:,1)=-v(:,-1)
p(:,1)=-p(:,-1)

u(:,size_y)=-u(:,size_y+1)+2. where u/v(:,size_y+1) is the ghost cell for top wall
v(:,size_y)=-v(:,size_y+1)
p(:,size_y)=-p(:,size_y+1)

There is a "+2" becos at the top wall, u=1. At other walls, dp/dn=0,

For left/right wall

du/dx=0,dv/dx=0

For bottom wall

du/dy=0,dv/dy=0

Is my implementation correcto? tks!
quarkz is offline   Reply With Quote

Old   April 9, 2009, 03:03
Default Pressure condition
  #2
New Member
 
Join Date: Mar 2009
Posts: 5
Rep Power: 8
CfdLunak is on a distinguished road
Hallo,

velocities look correct, but try no negative value of pressure anywhere e.g. p(1,=p(-1, instead of p(1,=-p(-1,.

For upper wall(the moving one) I used simple velocity prescription e.g. u=1, v=0.

Good luck.

L.
CfdLunak is offline   Reply With Quote

Old   April 12, 2009, 01:18
Default
  #3
Senior Member
 
TWB
Join Date: Mar 2009
Posts: 105
Rep Power: 8
quarkz is on a distinguished road
Oh ya, the pressure was a typo error. it's just dp/dn=0. For the upper wall, due to the staggered grid implementation,the wall is actually positioned in between the u(:,size_y) and u(:,size_y+1), hence the implementation is u(y)+u(y+1)/2 = 1.0, where 1.0 is the wall velocity.

Also, I'm using the fractional method. Hence, there's the momentum eqn and poisson eqn. I also need to update the velocity at some point in time. Does it mean that I must update the ghost pts constantly?
quarkz is offline   Reply With Quote

Old   April 12, 2009, 14:17
Default
  #4
otd
Member
 
private
Join Date: Mar 2009
Posts: 74
Rep Power: 8
otd is on a distinguished road
Don't you mean

[u(y+) + u(y-)]/2 = u(boundary) (1.0 for your problem)?
otd is offline   Reply With Quote

Old   April 13, 2009, 10:51
Default
  #5
Senior Member
 
TWB
Join Date: Mar 2009
Posts: 105
Rep Power: 8
quarkz is on a distinguished road
Ya tt's right otd! Anyway, I managed to get it working now, after correcting a small careless mistake. Tks!
quarkz is offline   Reply With Quote

Old   August 22, 2011, 17:16
Default
  #6
New Member
 
Vignesh Suresh
Join Date: Jul 2011
Posts: 10
Rep Power: 5
illuminati5288 is on a distinguished road
Hey quarkz....can you please tell me what the mistake was since I am having the same problem. My residuals keep increasing after the 15th iteration.

Thanks
illuminati5288 is offline   Reply With Quote

Old   August 22, 2011, 18:04
Default
  #7
Senior Member
 
TWB
Join Date: Mar 2009
Posts: 105
Rep Power: 8
quarkz is on a distinguished road
Hi illuminati5288, I am sorry I can't remember since it's 2 yrs ago. I am sure you will find your mistake if you looked carefully. Good luck!
quarkz is offline   Reply With Quote

Old   August 23, 2011, 15:46
Default
  #8
New Member
 
Vignesh Suresh
Join Date: Jul 2011
Posts: 10
Rep Power: 5
illuminati5288 is on a distinguished road
Hey quarkz can you remember if it was due to some mistake in the boundary conditions or did have something to do with the signs? I am just unable to figure out why my code diverges. I have followed the steps as illustrated in some of the papers and in the books.

If you dont mind can you please check my code and tell me where I am going wrong.

Thanks.


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

using namespace std;

double Re=100.0;
double nu;
int grid_size=32;
int counter,i,j;
double dx,dy,dt,old,residual_L2,residual,g,pressure_resid ual,u_residual,v_residual,velocity_residual,x,y;
double ux,vy,Fxx,Gyy,Hx_xy,Hy_xy;
double Lx=1.0;
double Ly=1.0;
double U=1.0;
double w=1.0;
double factor=pow(10.0,-9.0);
double factor2=pow(10.0,-4.0);
double factor3=pow(10.0,-4.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],old_v[1000][1000];
double p_ref,q,qh,phi;
double nodal_u[1000][1000],nodal_v[1000][1000],nodal_p[1000][1000];
double stream[1000][1000];
int flag=0;

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 nodal_velocities();
void stream_fucntion();
void output();

ofstream fout ("OUTPUT.dat");

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

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

void calculation()
{
system("CLS");
cout<<"SOLVING PRESSURE POISSON EQUATION"<<endl<<endl;
cout<<"PLEASE WAIT ";
do
{
pressure_residual=0.0;
u_residual=0.0;
v_residual=0.0;
velocity_residual=0.0;

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

for(i=1;i<grid_size+1;i++)
{
for(j=2;j<grid_size+1;j++)
{
old_u[i][j]=u[i][j];
}
}

for(i=2;i<grid_size+1;i++)
{
for(j=1;j<grid_size+1;j++)
{
old_v[i][j]=v[i][j];
}
}

interior_pressure_poisson();
pressure_correction();
update_uv();
pressure_bc();
time+=dt;

for(i=0;i<grid_size+2;i++)
pressure_residual+=fabs(p[i][0]-old_p[i][0])+fabs(p[i][grid_size+1]-old_p[i][grid_size+1]);

for(j=1;j<grid_size+1;j++)
pressure_residual+=fabs(p[0][j]-old_p[0][j])+fabs(p[grid_size+1][j]-old_p[grid_size+1][j]);

pressure_residual=pressure_residual*(1/((grid_size*4.0)+4.0));

cout<<endl<<pressure_residual;

if(pressure_residual<=factor2)
{
for(i=1;i<grid_size+1;i++)
{
for(j=2;j<grid_size+1;j++)
{
u_residual+=fabs(u[i][j]-old_u[i][j]);
}
}

u_residual=u_residual*(1/((grid_size)*(grid_size-1)));

for(i=2;i<grid_size+1;i++)
{
for(j=1;j<grid_size+1;j++)
{
v_residual+=fabs(v[i][j]-old_v[i][j]);
}
}

v_residual=v_residual*(1/((grid_size)*(grid_size-1)));

if(u_residual>v_residual)
velocity_residual=u_residual;
else
velocity_residual=v_residual;

if(velocity_residual<=factor3)
flag=1;
}
cout<<".";
}while(flag!=1);
cout<<endl<<endl<<"CONVERGENCE IS ATTAINED AT TIME t="<<time<<" seconds"<<endl;
}

void input()
{
nu=(U*Lx)/Re;
dx=Lx/grid_size;
dy=Ly/grid_size;

for(i=0;i<grid_size+2;i++)
for(j=0;j<grid_size+2;j++)
{
F[i][j]=0.0;
G[i][j]=0.0;
}

for(i=0;i<grid_size+3;i++)
for(j=0;j<grid_size+3;j++)
{
stream[i][j]=0.0;
nodal_p[i][j]=0.0;
nodal_u[i][j]=0.0;
nodal_v[i][j]=0.0;
Hx[i][j]=0.0;
Hy[i][j]=0.0;
}
}

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 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=1;j<grid_size+1;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 F_QUICK()
{
for(j=0;j<grid_size+2;j++)
{
for(i=grid_size+1;i>=0;i--)
{
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;
}
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;
}
else
{
ff=0.0;
}

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

void G_QUICK()
{
for(j=0;j<grid_size+2;j++)
{
for(i=grid_size+1;i>=0;i--)
{
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;
}
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;
}
else
{
gf=0.0;
}

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

void Hx_QUICK()
{
for(j=1;j<grid_size+2;j++)
{
for(i=grid_size+1;i>=1;i--)
{
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;
}
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;
}
else
{
hx=0.0;
}

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

void Hy_QUICK()
{
for(j=1;j<grid_size+2;j++)
{
for(i=grid_size+1;i>=1;i--)
{
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;
}
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;
}
else
{
hy=0.0;
}

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

void interior_pressure_poisson()
{
for(counter=1;counter<=1000000;counter++)
{
residual_L2=0.0;
for(j=1;j<grid_size+1;j++)
{
for(i=grid_size;i>=1;i--)
{
old=p[i][j];
ux=(u[i][j+1]-u[i][j]);
vy=(v[i][j]-v[i+1][j]);
Fxx=(F[i][j-1]-(2.0*F[i][j])+F[i][j+1]);
Hx_xy=(((Hx[i][j+1]-Hx[i][j]))-((Hx[i+1][j+1]-Hx[i+1][j])));
Hy_xy=(((Hy[i][j+1]-Hy[i][j]))-((Hy[i+1][j+1]-Hy[i+1][j])));
Gyy=(G[i-1][j]-(2.0*G[i][j])+G[i+1][j]);
g=((1.0/dt)*((ux*dx)+(vy*dy)))-(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 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 update_uv()
{
for(j=2;j<grid_size+1;j++)
{
for(i=grid_size;i>=1;i--)
{
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(j=1;j<grid_size+1;j++)
{
for(i=grid_size;i>=2;i--)
{
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 nodal_velocities()
{
for(i=1;i<grid_size+2;i++)
{
for(j=1;j<grid_size+2;j++)
{
nodal_u[i][j]=(u[i-1][j]+u[i][j])/2.0;
nodal_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*nodal_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*nodal_v[i][j]);
}
}
}

void output()
{
nodal_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"<<nodal_u[i][j]<<"\t"<<nodal_v[i][j]<<"\t"<<nodal_p[i][j]<<"\t"<<stream[i][j]<<endl;
x=x+dx;
}
y=y-dy;
}
}
illuminati5288 is offline   Reply With Quote

Old   August 26, 2011, 14:35
Default
  #9
New Member
 
Vignesh Suresh
Join Date: Jul 2011
Posts: 10
Rep Power: 5
illuminati5288 is on a distinguished road
hey quarkz,

I have kind of narrowed down the cause for divergence. When a moving lid is horizontal then the divergence occurs in the x-velocity but no divergence in the y-velocity. The vice-versa occurs when the moving lid is vertical. Any ideas why this is happening.

Thanks
illuminati5288 is offline   Reply With Quote

Old   January 20, 2013, 06:54
Default Reynolds Effects
  #10
New Member
 
Aryan
Join Date: Jan 2013
Posts: 8
Rep Power: 4
aryan is on a distinguished road
Hi Vignesh
It might be occured because of the high Re regime. I suggest check your program for the Re = 1

Good Luck
aryan 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


Similar Threads
Thread Thread Starter Forum Replies Last Post
inlet velocity boundary condition murali CFX 5 August 3, 2012 08:56
lid driven cavity sudha Main CFD Forum 6 November 21, 2010 06:11
cubic/ lid driven cavity Adrian Main CFD Forum 0 October 24, 2008 18:41
unsteady results lid driven cavity problem student Main CFD Forum 1 July 20, 2008 13:50
Ignore cells on partition boundary Karl FLUENT 7 May 11, 2002 22:12


All times are GMT -4. The time now is 23:15.