[Dhairya]

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