CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   Main CFD Forum (http://www.cfd-online.com/Forums/main/)
-   -   Boundary condition for 2d lid driven cavity using ghost cells (http://www.cfd-online.com/Forums/main/63384-boundary-condition-2d-lid-driven-cavity-using-ghost-cells.html)

 quarkz April 7, 2009 03:12

Boundary condition for 2d lid driven cavity using ghost cells

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!

 CfdLunak April 9, 2009 03:03

Pressure condition

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.

 quarkz April 12, 2009 01:18

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?

 otd April 12, 2009 14:17

Don't you mean

[u(y+) + u(y-)]/2 = u(boundary) (1.0 for your problem)?

 quarkz April 13, 2009 10:51

Ya tt's right otd! Anyway, I managed to get it working now, after correcting a small careless mistake. Tks!

 illuminati5288 August 22, 2011 17:16

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

 quarkz August 22, 2011 18:04

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!

 illuminati5288 August 23, 2011 15:46

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 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;
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);
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 August 26, 2011 14:35

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

 aryan January 20, 2013 06:54

Reynolds Effects

Hi Vignesh
It might be occured because of the high Re regime. I suggest check your program for the Re = 1

Good Luck

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