|
[Sponsors] |
June 26, 2016, 02:39 |
Lid driven cavity in Fortran
|
#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 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, 10:11 |
|
#2 |
Member
Mianzhi Wang
Join Date: Jan 2015
Location: Columbus, IN
Posts: 34
Rep Power: 11 |
first problem: you forgot to initialize alphav
|
|
June 28, 2016, 19:24 |
|
#3 |
Senior Member
Join Date: Oct 2011
Posts: 242
Rep Power: 17 |
+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 |
|
June 28, 2016, 23:05 |
|
#4 |
New Member
Vyas
Join Date: Jun 2016
Posts: 7
Rep Power: 10 |
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, 02: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 |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
Lid Driven cavity Fortran code | Dhairya | CFD Freelancers | 0 | June 26, 2016 01:48 |
Lid Driven cavity Fortran code | Dhairya | CFD Freelancers | 0 | June 26, 2016 01:42 |
CFX11 + Fortran compiler ? | Mohan | CFX | 20 | March 30, 2011 19:56 |
is there any parallel code for the famous Lid Driven Cavity flow? | gholamghar | Main CFD Forum | 0 | August 1, 2010 02:55 |
Lid driven cavity flow | KK | FLUENT | 1 | December 16, 2009 19:10 |