[quarkz]

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]

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]

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]

Don't you mean

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

[quarkz]

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

[illuminati5288]

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]

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]

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]

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]

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

Good Luck