# Lid-Driven cavity flow with SIMPLE method

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 August 9, 2005, 13:55 Lid-Driven cavity flow with SIMPLE method #1 luckyxu Guest   Posts: n/a Sponsored Links Hi, I am a newbie of CFD. I am writing FORTRAN code using SIMPLE algorithm to solve a lid-driven cavity flow problem. I followed the procedure in Anderson's book(Computational Fluid Dynamics:The Basics with Applications),pp.435-442,pp.254-264. I used staggered grids, Finite Difference Method, 50x50 uniform grids. But I just could not get the correct results. The velocity vector I plotted has no vortex at all. I have been working on this for several weeks and could not figure out what's wrong. Could you please help me check my code? I was almost driven mad by this. The notations are as followed: vb: v_bar vbb:v_bar_bar ux:u* vx:v* px* pp' rux:rho_u* rvx_rho_v* Code: PROGRAM driven_cavity IMPLICIT NONE REAL,DIMENSION(51,51) :: stream,p,px,x(51),y(51),d0,& pp,ppm,ppb,u(52,51),v(53,52),ux(52,51),vx(53,52),r ux(52,51),rvx(53,52) INTEGER :: is_beg,is_end,js_beg,js_end,nstep,cp,i,j,nstep1,ns tep2 REAL :: delx,dely,delt,rho,mu,& a0,b0,c0,fct1,fct2,fct3,delx2,dely2,vb,vbb,& ruu53,ruu33,ruu43,ruu42,uu1,uu2,ax,bx,ub,ubb,rvv45 ,rvv43,rvv54,rvv34,vv1,vv2,& ruv44,ruv42,rvu54,rvu34,height,width,Error0,Error1 ,Error2,Error00,Error3,& Error_stream,temp DATA is_beg,is_end,js_beg,js_end,rho,mu,delt,fct1,fct2, fct3,height,width& /1,51,1,51,998.2,1.005e-3,0.0001,0.8,0.8,0.5,0.1,0.1/ delx=width/REAL(is_end-1) dely=height/REAL(js_end-1) delx2=delx**2 dely2=dely**2 !assign initial values ux=0.0 u=0.0 u(is_beg:is_end+1,js_end)=2.0e-2 rux=rho*ux vx=0.0 v=0.0 rvx=rho*vx px=0.0 p=0.0 pp=0.0 nstep=0 ppm=0.0 stream=0.0 Error1=1. DO WHILE (Error1>9.5e-7) nstep=nstep+1 ! ! Pressure correction (SIMPLE) algorithm ! step 1. solve ru* and rv* ! a. for ru and u, pp.439(9.58) DO i=is_beg,is_end-1 DO j=js_beg+1,js_end-1 vb=0.5*(v(i+1,j+1)+v(i+2,j+1)) vbb=0.5*(v(i+1,j)+v(i+2,j)) ruu53=rho*u(i+2,j)*u(i+2,j) ruu33=rho*u(i,j)*u(i,j) ruv44=rho*u(i+1,j+1)*vb ruv42=rho*u(i+1,j-1)*vbb uu1=(u(i+2,j)-2.*u(i+1,j)+u(i,j))/delx2 uu2=(u(i+1,j+1)-2.*u(i+1,j)+u(i+1,j-1))/dely2 ax=-((ruu53-ruu33)/(2.*delx)+(ruv44-ruv42)/(2.*dely))+mu*(uu1+uu2) rux(i+1,j)=rux(i+1,j)+ax*delt-delt/delx*(px(i+1,j)-px(i,j)) ux(i+1,j)=rux(i+1,j)/rho END DO END DO !Boudary values ux(is_beg,js_beg+1:js_end-1)=0.0 ux(is_end+1,js_beg+1:js_end-1)=0.0 ux(is_beg:is_end+1,js_beg)=0.0 ux(is_beg:is_end+1,js_end)=2.0e-2 ! !b. for rv and v, pp.440(9.59) DO i=is_beg,is_end DO j=js_beg,js_end-1 ub=0.5*(u(i+1,j)+u(i+1,j+1)) ubb=0.5*(u(i,j)+u(i,j+1)) rvv45=rho*v(i+1,j+2)*v(i+1,j+2) rvv43=rho*v(i+1,j)*v(i+1,j) rvu54=rho*v(i+2,j+1)*ub rvu34=rho*v(i,j+1)*ubb vv1=(v(i+2,j+1)-2.*v(i+1,j+1)+v(i,j+1))/delx2 vv2=(v(i+1,j+2)-2.*v(i+1,j+1)+v(i+1,j))/dely2 bx=-((rvv45-rvv43)/(2.*dely)+(rvu54-rvu34)/(2.*delx))+mu*(vv1+vv2) rvx(i+1,j+1)=rvx(i+1,j+1)+bx*delt-delt/dely*(px(i,j+1)-px(i,j)) vx(i+1,j+1)=rvx(i+1,j+1)/rho END DO END DO !Boundary values vx(is_end+2,js_beg+1:js_end)=0.0 vx(is_beg,js_beg+1:js_end)=0.0 vx(is_beg:is_end+2,js_beg)=0.0 vx(is_beg:is_end+2,js_end+1)=0.0 ! ! step 2. solve p' from pressure correction formula ! by relaxation method (inner iteration) ! !calculate pressure corrector, pp.441(9.60) a0=2.0*(delt/delx2+delt/dely2) b0=-delt/delx2 c0=-delt/dely2 DO i=is_beg+1,is_end-1 DO j=js_beg+1,js_end-1 pp(i,j)=0.0 d0(i,j)=(rux(i+1,j)-rux(i,j))/delx+(rvx(i+1,j+1)-rvx(i+1,j))/dely END DO END DO ! calculate error of d0 Error1=0.0 DO i=is_beg+1,is_end-1 DO j=js_beg+1,js_end-1 Error0=ABS(ux(i,j)-u(i,j))!/ABS(u(i,j)) Error1=MAX(Error1,Error0) END DO END DO !calculate p' using relaxation method Error2=1.0 nstep1=0 DO WHILE (Error2>1.0e-7) nstep1=nstep1+1 Error2=0.0 DO i=is_beg+1,is_end-1 DO j=js_beg+1,js_end-1 ppm(i,j)=-1./a0*(b0*pp(i+1,j)+b0*ppm(i-1,j)+c0*pp(i,j+1)+& c0*ppm(i,j-1)+d0(i,j)) !ppm(i,j)=pp(i,j)+fct1*(ppb(i,j)-pp(i,j)) Error00=ABS(ppm(i,j)-pp(i,j))/ABS(pp(i,j)) Error2=MAX(Error2,Error00) pp(i,j)=ppm(i,j) END DO END DO END DO ! ! ! step 3. calculate p at n+1 step ! DO i=is_beg+1,is_end-1 DO j=js_beg+1,js_end-1 p(i,j)=px(i,j)+fct2*ppm(i,j) END DO END DO !update boundary values p(is_beg:is_end,js_beg)=p(is_beg:is_end,js_beg+1) p(is_beg:is_end,js_end)=p(is_beg:is_end,js_end-1) p(is_beg,js_beg+1:js_end-1)=p(is_beg+1,js_beg+1:js_end-1) p(is_end,js_beg+1:js_end-1)=p(is_end-1,js_beg+1:js_end-1) ! step 4. calculate u and v at n+1 step,from pp.239(6.100,6,101) ! DO i=is_beg,is_end-1 DO j=js_beg,js_end u(i+1,j)=ux(i+1,j)-delt/delx/rho*(pp(i+1,j)-pp(i,j)) END DO END DO DO i=is_beg,is_end DO j=js_beg,js_end-1 v(i+1,j+1)=vx(i+1,j+1)-delt/dely/rho*(pp(i,j+1)-pp(i,j)) END DO END DO ! ! update boundary conditions ! u(is_beg:is_end+1,js_end)=2.0e-2 u(is_beg:is_end+1,js_beg)=0. v(is_beg:is_end+2,js_end+1)=0. v(is_beg:is_end+2,js_beg)=0. ! u(is_beg,js_beg:js_end)=0.0 u(is_end+1,js_beg:js_end)=0.0 v(is_beg,js_beg:js_end+1)=0.0 v(is_end+2,js_beg:js_end+1)=0.0 ! step 5. calculate stream function Error3=1.0 nstep2=0 DO WHILE (Error3>1.0e-2) nstep2=nstep2+1 Error3=0.0 DO i=is_beg+1,is_end DO j=js_beg+1,js_end temp=stream(i,j) stream(i,j)=0.5*(stream(i-1,j)+stream(i,j-1))+0.125*dely*& (u(i,j)+u(i+1,j)+u(i,j-1)+u(i+1,j-1))-0.125*delx*&(v(i,j)+v(i+1,j)+v(i,j+1)+v(i+1,j+1)) Error_stream=ABS(stream(i,j)-temp)/ABS(temp) Error3=MAX(Error_stream,Error3) END DO END DO stream(is_beg,js_beg:js_end)=stream(is_beg+1,js_be g:js_end) stream(is_beg:is_end,js_beg)=stream(is_beg:is_end, js_beg+1) END DO ! ! step 6. update pressure field px=p WRITE (*,*) nstep,Error1,nstep1,Error2,nstep2,Error3 END DO OPEN(25,file="C:\Research\LSM\data.dat") WRITE (25,*) 'VARIABLES="X","Y","P"' WRITE (25,*) 'ZONE I=51, J=51 , F=POINT' DO i=is_beg,is_end DO j=js_beg,js_end x(i)=REAL(i)*delx y(j)=REAL(j)*dely WRITE(25,"(2f10.3,f36.30)") x(i),y(j),p(i,j) END DO END DO OPEN(26,file="C:\Research\LSM\data_u1.dat") DO i=is_beg,is_end DO j=js_beg,js_end x(i)=REAL(i)*delx y(j)=REAL(j)*dely WRITE(26,"(2f10.3,2f36.30)") x(i),y(j),u(i,j),v(i,j) END DO END DO OPEN(27,file="C:\Research\LSM\stream.dat") WRITE (27,*) 'VARIABLES="X","Y","STREAM"' WRITE (27,*) 'ZONE I=51, J=51 , F=POINT' DO i=is_beg,is_end DO j=js_beg,js_end x(i)=REAL(i)*delx y(j)=REAL(j)*dely WRITE(27,"(2f10.3,f36.30)") x(i),y(j),stream(i,j) END DO END DO END PROGRAM driven_cavity Thank you in advance.

 June 22, 2011, 03:16 fnm simple #2 New Member   iman Join Date: Dec 2010 Location: Iran Posts: 5 Rep Power: 8 hi I need a fortran for simple method @ cavity. if you find any code that had been correct result please send me. thanks

 June 22, 2011, 23:42 doubt #3 Member   kabilan.B Join Date: Jul 2010 Location: chennai Posts: 79 Rep Power: 9 in the above program from bottom 3rd line "f13.36" what is indicate that

 June 23, 2011, 08:22 #4 New Member   anil Join Date: Jun 2011 Location: chennai Posts: 20 Rep Power: 8 i hav done d research work on square lid cavity using fluent n gambit..but instead of using SIMPLE scheme i used SIMPLEC scheme...i even compared it wid ghia n ghia....

June 24, 2011, 12:29
#5
New Member

iman
Join Date: Dec 2010
Location: Iran
Posts: 5
Rep Power: 8
Quote:
 Originally Posted by kabilan in the above program from bottom 3rd line "f13.36" what is indicate that
hi.
the above program dont have a correct result.
I need another code.

 June 25, 2011, 07:38 #6 New Member   anil Join Date: Jun 2011 Location: chennai Posts: 20 Rep Power: 8 i hav word on lid driven cavity but with fluent.....n got d same results as gvn by ghia n ghia..

 November 9, 2011, 08:17 #7 New Member   Join Date: Sep 2011 Posts: 3 Rep Power: 7 I hope someone can upload a softcopy of the following work: Ghia, U., Ghia, K.N. & Shin, C.T. 1982, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, Journal of Computational Physics vol. 48, pp. 387-411 I don't have free access to this Sciencedirect article from my university's subscription.

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post klw Main CFD Forum 8 May 21, 2011 04:57 abcdef123 Main CFD Forum 0 February 26, 2010 09:24 kieranhood OpenFOAM Paraview & paraFoam 0 February 13, 2010 17:28 Mathias Krause Main CFD Forum 5 October 17, 2006 06:34 yfyap Main CFD Forum 3 September 19, 2001 18:16