# Boundary condition for 2d lid driven cavity using ghost cells

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

 April 7, 2009, 03:12 Boundary condition for 2d lid driven cavity using ghost cells #1 Senior Member   TWB Join Date: Mar 2009 Posts: 402 Rep Power: 19 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!

 April 9, 2009, 03:03 Pressure condition #2 New Member   Join Date: Mar 2009 Posts: 5 Rep Power: 17 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.

 April 12, 2009, 01:18 #3 Senior Member   TWB Join Date: Mar 2009 Posts: 402 Rep Power: 19 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?

 April 12, 2009, 14:17 #4 Member   private Join Date: Mar 2009 Posts: 74 Rep Power: 17 Don't you mean [u(y+) + u(y-)]/2 = u(boundary) (1.0 for your problem)?

 April 13, 2009, 10:51 #5 Senior Member   TWB Join Date: Mar 2009 Posts: 402 Rep Power: 19 Ya tt's right otd! Anyway, I managed to get it working now, after correcting a small careless mistake. Tks!

 August 22, 2011, 17:16 #6 New Member   Vignesh Suresh Join Date: Jul 2011 Posts: 12 Rep Power: 15 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

 August 22, 2011, 18:04 #7 Senior Member   TWB Join Date: Mar 2009 Posts: 402 Rep Power: 19 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!

 August 23, 2011, 15:46 #8 New Member   Vignesh Suresh Join Date: Jul 2011 Posts: 12 Rep Power: 15 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 #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"<v_residual) velocity_residual=u_residual; else velocity_residual=v_residual; if(velocity_residual<=factor3) flag=1; } cout<<"."; }while(flag!=1); cout<=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=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=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=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=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=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=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=1;i--) { for(j=1;j=1;i--) { for(j=2;j

 August 26, 2011, 14:35 #9 New Member   Vignesh Suresh Join Date: Jul 2011 Posts: 12 Rep Power: 15 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

 January 20, 2013, 05:54 Reynolds Effects #10 New Member   Aryan Join Date: Jan 2013 Posts: 8 Rep Power: 13 Hi Vignesh It might be occured because of the high Re regime. I suggest check your program for the Re = 1 Good Luck