# Lid driven cavity in Fortran

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

 June 26, 2016, 01:39 Lid driven cavity in Fortran #1 New Member   Vyas Join Date: Jun 2016 Posts: 7 Rep Power: 8 Hello, everyone I am a student who is trying to learn cfd on his own. I am trying to simulate lid riven cavity problem using fvm in fortran and have been working on it since last two weeks but still am not getting the results. I am using collocated grids along with Rhie-Chow interpolation modified by Majumdar. First order upwind scheme has been used for convective terms and for pressure velocity decoupling SIMPLE algorithm is used. I have posted my code below, could someone please go through it and help me find where I am going wrong. I am getting lot of NaN in output and in few iterations the pressure and velocity rise to extremely large values. program liddrivencavity implicit none integer i,j,k,nx,ny,iter real, dimension(100,100) :: u,v,p,ae,aw,an,as,ap,Fe,Fw,Fs,Fn,Atotu,Atotv,de,dn real, dimension(100,100) :: ueterm1,ueterm2,ueterm3,ue,vnterm1,vnterm2,vnterm3 ,vn real, dimension(100,100) :: p1,app,aep,awp,asp,anp,bp real, dimension(100):: x,y,tridu,tridl,tridd,tridb real dx,dy,alphau,alphav,alphap,rho,uref,mew rho=400. mew=1. uref=1. alphau=0.3 alphap=0.1 nx=5 !Num of pts in x dir ny=5 dx=1./(real(nx)-1.) dy=1./(real(ny)-1.) !GENERATING GRID do i=1,nx x(i)=(i-1.)*dx end do do j=1,ny y(j)=(j-1.)*dy end do !BOUNDARY CONDITIONS do i=1,nx !horizontal walls u(i,1)=0. v(i,1)=0. u(i,ny)=1. v(i,ny)=0. end do do j=1,ny !verticle walls u(1,j)=0. v(1,j)=0. u(nx,j)=0. v(nx,j)=0. end do !INITIALIZING do i=2,nx-1 do j=2,ny-1 u(i,j)=0.1 v(i,j)=0.1 p(i,j)=0. end do end do do i=1,nx-1 do j=1,ny-1 ue(i,j)=0. vn(i,j)=0. p1(i,j)=0. end do end do master: do iter=1,1000 !Fe, Fw, Fn and Fs do i=2,nx-1 do j=2,ny-1 Fe(i,j)=rho*dy*0.5*(u(i,j)+u(i+1,j)) Fw(i,j)=rho*dy*0.5*(u(i,j)+u(i-1,j)) Fn(i,j)=rho*dx*0.5*(v(i,j)+v(i,j+1)) Fs(i,j)=rho*dx*0.5*(v(i,j)+v(i,j-1)) end do end do !GENERATING anb do i=2,nx-1 do j=2,ny-1 ae(i,j)=(mew*dy/dx)+max(-Fe(i,j),0.) aw(i,j)=(mew*dy/dx)+max(Fw(i,j),0.) an(i,j)=(mew*dx/dy)+max(-Fn(i,j),0.) as(i,j)=(mew*dx/dy)+max(Fs(i,j),0.) ap(i,j)=ae(i,j)+aw(i,j)+as(i,j)+an(i,j)+Fe(i,j)-Fw(i,j)+Fn(i,j)-Fs(i,j) Atotu(i,j)=(ae(i,j)*u(i+1,j)+aw(i,j)*u(i-1,j)+as(i,j)*u(i,j-1)+an(i,j)*u(i,j+1))/ap(i,j) Atotv(i,j)=(ae(i,j)*v(i+1,j)+aw(i,j)*v(i-1,j)+as(i,j)*v(i,j-1)+an(i,j)*v(i,j+1))/ap(i,j) end do end do do i=1,nx-1 do j=1,ny-1 de(i,j)=0. dn(i,j)=0. end do end do !FINDING ue FOR i=2 TO Nx-2 AND j=2 TO Ny-1 USING RHIE-CHOW INTERPOLATION !ue=ueterm1+ueterm2+ueterm3 do i=2,nx-2 do j=2,ny-1 !0.5 is 1/2 de(i,j)=dy*alphau*((1/ap(i,j))+(1/ap(i+1,j)))*0.5 ueterm1(i,j)=0.5*(Atotu(i,j)+Atotu(i+1,j))*alphau ueterm2(i,j)=(p(i,j)-p(i+1,j))*de(i,j) do k=1,100 !MAJUMDAR'S IDEA TO MAKE IT INDEPENDENT OF VALUE OF ALPHA ueterm3(i,j)=(1.-alphau)*(ue(i,j)-0.5*(u(i+1,j)+u(i,j))) ue(i,j)=ueterm1(i,j)+ueterm2(i,j)+ueterm3(i,j) end do end do end do !FINDING vn FOR i=2 TO Nx-1 AND j=2 TO Ny-2 USING RHIE-CHOW INTERPOLATION !vn=vnterm1+vnterm2+vnterm3 do i=2,nx-1 do j=2,ny-2 !0.5 is 1/2 dn(i,j)=dx*alphav*((1/ap(i,j))+(1/ap(i,j+1)))*0.5 vnterm1(i,j)=0.5*(Atotv(i,j)+Atotv(i,j+1))*alphav vnterm2(i,j)=(p(i,j)-p(i,j+1))*dn(i,j) do k=1,100 !MAJUMDAR'S IDEA TO MAKE IT INDEPENDENT OF VALUE OF ALPHA vnterm3(i,j)=(1.-alphav)*(vn(i,j)-0.5*(v(i,j+1)+v(i,j))) vn(i,j)=vnterm1(i,j)+vnterm2(i,j)+vnterm3(i,j) end do end do end do !FINDINF p1 do i=2,nx-1 do j=2,ny-1 awp(i,j)=de(i-1,j)*dy aep(i,j)=de(i,j)*dy asp(i,j)=dn(i,j-1)*dx anp(i,j)=dn(i,j)*dx app(i,j)=awp(i,j)+aep(i,j)+asp(i,j)+anp(i,j) bp(i,j)=(ue(i-1,j)-ue(i,j))*dy+(vn(i,j-1)-vn(i,j))*dx end do end do !do k=1,1000 !do i=2,nx-1 ! do j=2,ny-1 ! p1(i,j)=(awp(i,j)*p1(i-1,j)+aep(i,j)*p1(i+1,j)+asp(i,j)*p1(i,j-1)+anp(i,j)*p1(i,j+1)+bp(i,j))/app(i,j) ! end do !end do !end do !TDMA do k=1,100 do j=2,ny-1 !generating trid matrix do i=2,nx-1 tridu(i)=-aep(i,j) tridd(i)=app(i,j) tridl(i)=-awp(i,j) tridb(i)=anp(i,j)*p1(i,j+1)+asp(i,j)*p1(i,j-1)+bp(i,j) end do tridu(nx-1)=0. tridl(2)=0. !TDMA do i=2,nx-2 tridd(i+1)=tridd(i+1)-(tridl(i+1)*tridu(i)/tridd(i)) tridb(i+1)=tridb(i+1)-(tridl(i+1)*tridb(i)/tridd(i)) end do p1(nx-1,j)=tridb(nx-1)/tridd(nx-1) i=nx-2 do p1(i,j)=(tridb(i)-tridu(i)*p1(i+1,j))/tridd(i) write(*,*) p1(i,j) i=i-1 if(i==1) then exit end if end do end do end do !correcting pressures do i=2,nx-1 do j=2,ny-1 p(i,j)=p(i,j)+alphap*p1(i,j) end do end do !correcting velocities !TAKING CARE OF POINS NEXT TO BOUNDARY BOUNDARIES do j=2,ny-1 p1(1,j)=p1(2,j) p1(nx,j)=p1(nx-1,j) end do do i=2,nx-1 p1(i,1)=p1(i,2) p1(i,ny)=p1(i,ny-1) end do do i=2,nx-1 do j=2,ny-1 u(i,j)=u(i,j)+alphau*dy*0.5*(p1(i-1,j)-p1(i+1,j))/ap(i,j) v(i,j)=v(i,j)+alphav*dx*0.5*(p1(i,j-1)-p1(i,j+1))/ap(i,j) end do end do end do master !OUTPUT open(unit=1,file='xvel.txt',status='unknown') open(unit=2,file='yvel.txt',status='unknown') open(unit=3,file='p.txt',status='unknown') open(unit=4,file='grid.txt',status='unknown') do i=1,nx do j=1,ny write(1,*) u(i,j) write(2,*) v(i,j) write(3,*) p(i,j) write(4,*) x(i),y(j) end do end do close(1) close(2) close(3) close(4) write(*,*) 'EXECUTION SUCESSFULL' end program

 June 28, 2016, 09:11 #2 Member   Mianzhi Wang Join Date: Jan 2015 Location: Columbus, IN Posts: 34 Rep Power: 10 first problem: you forgot to initialize alphav Dhairya likes this.

 June 28, 2016, 18:24 #3 Senior Member   Join Date: Oct 2011 Posts: 213 Rep Power: 14 +1 for alphav When developping a code, you should always enable check flags at compilation. This can tell you a lot about common mistakes. Here come some of them for the gnu compiler: gfortran -fbounds-check -fbacktrace -fcheck=all -ffpe-trap=invalid,zero,overflow,underflow -O0 -g -Wuninitialized In particular, the -Wuninitialized tells you alphav is not initialized before even running your executable Dhairya likes this.

 June 28, 2016, 22:05 #4 New Member   Vyas Join Date: Jun 2016 Posts: 7 Rep Power: 8 Thank you very much for helping me out, after reading the replies i immediately took alphav=0.3 but then also my solution is diverging. I even tried the same problem using stream-vorticity approach but there also I am facing the issue of divergence even by using first order upwind scheme.

 June 29, 2016, 01:25 #5 New Member   Join Date: Jun 2016 Posts: 2 Rep Power: 0 Did you check if your coefficient matrix is diagonally dominant?

 Tags collocated grid, fortran, lid driven cavity, rhie-chow simple