CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   Main CFD Forum (https://www.cfd-online.com/Forums/main/)
-   -   Lid-Driven cavity flow with SIMPLE method (https://www.cfd-online.com/Forums/main/9669-lid-driven-cavity-flow-simple-method.html)

luckyxu August 9, 2005 13:55
Lid-Driven cavity flow with SIMPLE method
 
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:p* pp:p' 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.

iman_d June 22, 2011 03:16
fnm simple
 
hi
I need a fortran for simple method @ cavity.
if you find any code that had been correct result please send me.
thanks

kabilan June 22, 2011 23:42
doubt
 
in the above program from bottom 3rd line "f13.36" what is indicate that

akj June 23, 2011 08:22
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....

iman_d June 24, 2011 12:29
Quote:
Originally Posted by kabilan (Post 313186)
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.
:)

akj June 25, 2011 07:38
i hav word on lid driven cavity but with fluent.....n got d same results as gvn by ghia n ghia..

student88 November 9, 2011 07:17
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.


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