|
[Sponsors] |
Boundary condition for 2d lid driven cavity using ghost cells |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
April 7, 2009, 04:12 |
Boundary condition for 2d lid driven cavity using ghost cells
|
#1 |
Senior Member
TWB
Join Date: Mar 2009
Posts: 414
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, 04: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, 02:18 |
|
#3 |
Senior Member
TWB
Join Date: Mar 2009
Posts: 414
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, 15: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, 11:51 |
|
#5 |
Senior Member
TWB
Join Date: Mar 2009
Posts: 414
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, 18: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, 19:04 |
|
#7 |
Senior Member
TWB
Join Date: Mar 2009
Posts: 414
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, 16: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<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; } } |
|
August 26, 2011, 15: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, 06: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 |
|
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
inlet velocity boundary condition | murali | CFX | 5 | August 3, 2012 09:56 |
cubic/ lid driven cavity | Adrian | Main CFD Forum | 0 | October 24, 2008 19:41 |
unsteady results lid driven cavity problem | student | Main CFD Forum | 1 | July 20, 2008 14:50 |
lid driven cavity | sudha | Main CFD Forum | 5 | July 11, 2007 04:05 |
Ignore cells on partition boundary | Karl | FLUENT | 7 | May 11, 2002 23:12 |