CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Lid-Driven cavity flow with SIMPLE method

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   August 9, 2005, 14:55
Default Lid-Driven cavity flow with SIMPLE method
  #1
luckyxu
Guest
 
Posts: n/a
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.

  Reply With Quote

Old   June 22, 2011, 04:16
Default fnm simple
  #2
New Member
 
iman
Join Date: Dec 2010
Location: Iran
Posts: 5
Rep Power: 15
iman_d is on a distinguished road
hi
I need a fortran for simple method @ cavity.
if you find any code that had been correct result please send me.
thanks
iman_d is offline   Reply With Quote

Old   June 23, 2011, 00:42
Default doubt
  #3
Member
 
kabilan.B
Join Date: Jul 2010
Location: chennai
Posts: 79
Rep Power: 15
kabilan is on a distinguished road
in the above program from bottom 3rd line "f13.36" what is indicate that
kabilan is offline   Reply With Quote

Old   June 23, 2011, 09:22
Default
  #4
akj
New Member
 
anil
Join Date: Jun 2011
Location: chennai
Posts: 20
Rep Power: 14
akj is on a distinguished road
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....
akj is offline   Reply With Quote

Old   June 24, 2011, 13:29
Default
  #5
New Member
 
iman
Join Date: Dec 2010
Location: Iran
Posts: 5
Rep Power: 15
iman_d is on a distinguished road
Quote:
Originally Posted by kabilan View Post
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.
iman_d is offline   Reply With Quote

Old   June 25, 2011, 08:38
Default
  #6
akj
New Member
 
anil
Join Date: Jun 2011
Location: chennai
Posts: 20
Rep Power: 14
akj is on a distinguished road
i hav word on lid driven cavity but with fluent.....n got d same results as gvn by ghia n ghia..
akj is offline   Reply With Quote

Old   November 9, 2011, 08:17
Default
  #7
New Member
 
Join Date: Sep 2011
Posts: 3
Rep Power: 14
student88 is on a distinguished road
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.

student88 is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

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 Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Lack of Recirculation for Lid Driven Cavity Flow klw Main CFD Forum 8 May 21, 2011 05:57
SIMPLE method for inviscid flow abcdef123 Main CFD Forum 0 February 26, 2010 09:24
[OpenFOAM] Paraview - Lid Driven Cavity kieranhood ParaView 0 February 13, 2010 17:28
Lid Driven Cavity 3D Mathias Krause Main CFD Forum 5 October 17, 2006 07:34
on cavity driven flow yfyap Main CFD Forum 3 September 19, 2001 19:16


All times are GMT -4. The time now is 18:09.