CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Lid driven cavity in Fortran

Register Blogs Community New Posts Updated Threads Search

Like Tree2Likes
  • 1 Post By wangmianzhi
  • 1 Post By naffrancois

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 26, 2016, 02:39
Default Lid driven cavity in Fortran
  #1
New Member
 
Vyas
Join Date: Jun 2016
Posts: 7
Rep Power: 10
Dhairya is on a distinguished road
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
Dhairya is offline   Reply With Quote

Old   June 28, 2016, 10:11
Default
  #2
Member
 
Mianzhi Wang
Join Date: Jan 2015
Location: Columbus, IN
Posts: 34
Rep Power: 11
wangmianzhi is on a distinguished road
first problem: you forgot to initialize alphav
Dhairya likes this.
wangmianzhi is offline   Reply With Quote

Old   June 28, 2016, 19:24
Default
  #3
Senior Member
 
Join Date: Oct 2011
Posts: 242
Rep Power: 17
naffrancois is on a distinguished road
+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.
naffrancois is offline   Reply With Quote

Old   June 28, 2016, 23:05
Default
  #4
New Member
 
Vyas
Join Date: Jun 2016
Posts: 7
Rep Power: 10
Dhairya is on a distinguished road
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.
Dhairya is offline   Reply With Quote

Old   June 29, 2016, 02:25
Default
  #5
New Member
 
Join Date: Jun 2016
Posts: 2
Rep Power: 0
ravachol is on a distinguished road
Did you check if your coefficient matrix is diagonally dominant?
ravachol is offline   Reply With Quote

Reply

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


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 Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


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


All times are GMT -4. The time now is 05:32.