# 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 Thread Tools Search this Thread Show Printable Version Email this Page Search this Thread: Advanced Search Display Modes Linear Mode Switch to Hybrid Mode Switch to Threaded 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 Off Pingbacks are On Refbacks are On Forum Rules Similar Threads Thread Thread Starter Forum Replies Last Post Dhairya CFD Freelancers 0 June 26, 2016 00:48 Dhairya CFD Freelancers 0 June 26, 2016 00:42 Mohan CFX 20 March 30, 2011 18:56 gholamghar Main CFD Forum 0 August 1, 2010 01:55 KK FLUENT 1 December 16, 2009 18:10

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