Lid Driven cavity Fortran code

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

 June 26, 2016, 00:42 Lid Driven cavity Fortran code #1 New Member   Vyas Join Date: Jun 2016 Posts: 7 Rep Power: 10 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 , collocated grid, upwind scheme and Rhie-Chow interpolation corrected by Majumdar in fortran and have been working on it since last two weeks but still am not getting the results. I have posted my code below, could someone please go through it and help me find where I am going wrong. 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

 Tags lid driven
 Thread Tools Search this Thread Search this Thread: Advanced Search 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 Off Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post josephlm Main CFD Forum 4 August 17, 2023 20:36 akin Main CFD Forum 1 June 16, 2012 07:53 gholamghar Main CFD Forum 0 August 1, 2010 01:55 [OpenFOAM] Paraview - Lid Driven Cavity kieranhood ParaView 0 February 13, 2010 16:28 KK FLUENT 1 December 16, 2009 18:10

All times are GMT -4. The time now is 17:36.

 Contact Us - CFD Online - Privacy Statement - Top