Hello!
I recently modify LBM code with FORTRAN by Cheng Ming, but I can not gain the correct results, and do not know wher there is a problem. Please master hand help me see and correct it, thank you.
! File : LBM_D2Q9.f90
! Author : Cheng Ming (chengm@ihpc.a-star.edu.sg)
! Date : 01-09-2006
! Revision : 1.0
! DESCRIPTION
! This is a Fortan 77 LBM_D2Q9 code for 2D viscous incompressible flow
! COPYRIGHT & LICENSE INFORMATION
! Copyright 2006 Cheng Ming and The Institute of High
! Performance Computing (A*STAR)
! This file is part of lbmSolver.
! lbmSolver is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 2 of the License, or
! (at your option) any later version.
! lbmSolver is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
! You should have received a copy of the GNU General Public License
! along with depSolver; if not, write to the Free Software
! Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
MODULE D2Q9Const
! D2Q9 Weights
Double precision,parameter :: w(0:8)=(/4.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0,&
&1.0d0/9.0d0,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0/)
! D2Q9 Directions
Integer,parameter :: dr(0:8,0:1) =(/(/0,1,0,-1,0,1,-1,-1,1/),(/0,0,1,0,-1,1,1,-1,-1/)/)
Integer,parameter :: opposite(0:8) = (/0,3,4,1,2,7,8,5,6/)
END MODULE D2Q9Const
! =======================================================
! Constants for simulation setup
! =======================================================
MODULE simParam
Integer, parameter :: Xmax =301 !Grid size in x dimension
Integer, parameter :: Ymax = 81 !Grid size in x dimension
Integer, parameter :: a = 21 !length of rectangular cylinder in x direction
Integer, parameter :: b = 21 !length of rectangular cylinder in y direction
Integer, parameter :: x1 = 80 !Distance from inlet boundary to cylinder
Integer, parameter :: x2 =x1+a !Distance from inlet boundary to cylinder
Integer, parameter :: y1 =(Ymax+1-b)/2 !Distance from inlet boundary to cylinder
Integer, parameter :: y2 =(Ymax+1+b)/2 !Distance from inlet boundary to cylinder
Integer, parameter :: MaxStep= 80000 ! Maximum number of iterations
Integer, parameter :: NUMAX =1000 ! Output frequency
Double precision, parameter :: Uw = 0.01 ! Incoming velocity (Uw < 0.1)
Double precision, parameter :: Re = 100 ! Reynolds number (for high Re a&b should be increased)
END MODULE simParam
PROGRAM LBM
! Varibles and parameters
USE simParam, ONLY : Xmax,Ymax,MaxStep, NUMAX
implicit none
Double precision :: omega !omega...Collision parameter
Double precision :: tau !Collision time
Double precision :: nu,s8,s9 !nu the kinetice viscosity of the fluid
Double precision, parameter :: density=0.1d0
Double precision, dimension(:,:,:), allocatable :: fIn,fEq,fOut ! density function, equib,array of temporare storage
Double precision, dimension(:,:), allocatable :: Dxy,u,v
Integer, dimension(:,:), allocatable :: wall ! Solid wall array
Integer :: istep !interation counter
Integer :: NUMB,x,y !computing control
!
! ==================================================================
! | Begin initialisation |
! ==================================================================
!
Allocate(fIn(0:8,Xmax,Ymax))
Allocate(fEq(0:8,Xmax,Ymax))
Allocate(fOut(0:8,Xmax,Ymax))
Allocate(Dxy(Xmax,Ymax))
Allocate(u(Xmax,Ymax))
Allocate(v(Xmax,Ymax))
Allocate(wall(Xmax,Ymax))
DO x = 1, Xmax
DO y = 1, Ymax
Dxy(x,y)=0.10d0
u(x,y) = 0.0d0
v(x,y) = 0.0d0
END DO
END DO
! define wall boundary
CALL Wall_coordinate(wall)
CALL computeOmega(omega)
! initialize function f
CALL Init_density(density,Dxy,u,v,fIn)
!
fEq = fIn
!
! ==================================================================
! | END initialisation |
! ==================================================================
! ==================================================================
! | Begin iterations |
! ==================================================================
!
NUMB=1
DO istep = 1, MaxStep
IF (NUMB<=NUMAX) THEN
NUMB=NUMB+1
END IF
!
! CALL computeMacros(fIn,Dxy,u,v)
! CALL computeFeq(Dxy,u,v,fEq)
! Inlet-outlet boundary condition
CALL Collision(omega,fIn,fEq,wall,fOut)
CALL Inlet_outlet_BC(fEq,u,v,Dxy,fIn)
! density Collision
!
! density propagation
CALL Stream(fOut,fIn)
! bounc back from solid wall
CALL Bounceback(wall,fIn,fOut)
CALL computeMacros(fIn,Dxy,u,v)
CALL computeFeq(Dxy,u,v,fEq)
! each MaxStep/10 iteration
IF (NUMB>NUMAX) THEN
! Output the data velocity, vorticity and pressure
CALL Dump(wall,Dxy,u,v,istep)
NUMB=1
END IF
END DO
!
Deallocate(fIn)
Deallocate(fOut)
Deallocate(fEq)
Deallocate(u)
Deallocate(v)
Deallocate(Dxy)
Deallocate(wall)
! ==================================================================
! | END iterations |
! ==================================================================
END PROGRAM LBM
!
! ==================================================================
! | Compute the relaxation parameter from the Reynolds nmber |
! ==================================================================
SUBROUTINE computeOmega(omega)
USE simParam, ONLY: Re, Uw, x1,x2
implicit none
Double precision:: omega !omega...Collision parameter
Double precision :: tau !Collision time
Double precision :: nu,s8,s9 !nu the kinetice viscosity of the fluid
Integer XYMAX
XYMAX=x2-x1
tau=Uw*XYMAX*3.0/Re + 0.5
omega=1.0/tau
! Re=UwXYMAX3.0/(1.0/omega-0.5)
nu=(2tau-1)/6
s8=2/(6nu+1)
s9=s8
write(*,*)'tau=',tau
PRINT*,'Re=',Re,'omega=',omega
PRINT*, 'nu=',nu,'s8=',s8,'s9=',s9
IF(tau.lt.0.506) THEN
PRINT*,'Instablility will take place at such tau < 0.506 !!!'
PRINT*,'tau=Uw*a*3.0/Re + 0.5'
PRINT*,'Pleae increase value of a or Uw.'
STOP
END IF
END SUBROUTINE computeOmega
! ****************************************************************
! * *
! * Input solid wall coordinates *
! * *
! ****************************************************************
!
SUBROUTINE Wall_coordinate(wall)
USE simParam, ONLY: Xmax,Ymax,x1,x2,y1,y2,a,b
IMPLICIT NONE
INTEGER :: wall(Xmax,Ymax)
INTEGER :: x,y
! no solid wall in wall array
DO y = 1, Ymax
DO x = 1, Xmax
wall(x,y) = 0
END DO
END DO
OPEN (13,file= 'outputImage.dat')
DO x = x1, x2
DO y = y1, y2
wall(x,y) = 1
wall(x,y) = 1
END DO
END DO
do X =1,Xmax
write(13,"(51I1,/)") wall(x,:)
END DO
close(13)
END SUBROUTINE Wall_coordinate
!
! ****************************************************************
! * *
! * Initialize density distribution function n with equilibrium *
! * for zero velocity *
! * *
! ****************************************************************
!
SUBROUTINE Init_density(density,Dxy,u,v,fIn)
USE D2Q9Const, ONLY: w,dr
USE simParam, ONLY: Xmax,Ymax
IMPLICIT NONE
Double precision:: density, Dxy(Xmax,Ymax),u(Xmax,Ymax),v(Xmax,Ymax),uxy,uSqr
Double precision:: fIn(0:8,Xmax,Ymax)
! local variables
INTEGER x,y,i
! loop over computational Domain
DO x = 1, Xmax
DO y = 1, Ymax
uSqr = u(x,y) * u(x,y) + v(x,y) * v(x,y)
DO i = 0, 8
uxy=u(x,y)dr(i,0)+v(x,y)dr(i,1)
fIn(i,x,y)= w(i)Dxy(x,y)(1.0d0+3.0d0uxy+4.5d0uxyuxy-1.5d0uSqr)
END DO
END DO
END DO
END SUBROUTINE Init_density
!
! ========================================================
! Compute density and velocity from distribution functions
! ========================================================
SUBROUTINE computeMacros(fIn,Dxy,u,v)
USE simParam, ONLY: Xmax,Ymax
implicit none
double precision :: fIn(0:8,Xmax,Ymax)
double precision :: u(Xmax,Ymax), v(Xmax,Ymax), Dxy(Xmax,Ymax)
integer x, y
do x = 1, Xmax
do y = 1, Ymax
Dxy(x,y) = fIn(0,x,y) + fIn(1,x,y) + fIn(2,x,y) + fIn(3,x,y) + fIn(4,x,y) + fIn(5,x,y) + fIn(6,x,y) + fIn(7,x,y) + fIn(8,x,y)
u(x,y) = (fIn(1,x,y) - fIn(3,x,y) + fIn(5,x,y) - fIn(6,x,y) - fIn(7,x,y) + fIn(8,x,y)) / Dxy(x,y)
v(x,y) = (fIn(2,x,y) - fIn(4,x,y) + fIn(5,x,y) + fIn(6,x,y) - fIn(7,x,y) - fIn(8,x,y)) / Dxy(x,y)
end do
end do
END SUBROUTINE computeMacros
!
! ****************************************************************
! * *
! * Compute equilibrium density function
! *
! * *
! ****************************************************************
SUBROUTINE computeFeq(Dxy,u,v,fEq)
USE D2Q9Const, ONLY: w,dr
USE simParam, ONLY: Xmax,Ymax
IMPLICIT NONE
Double precision :: Dxy(Xmax,Ymax),u(Xmax,Ymax),v(Xmax,Ymax)
Double precision :: fEq(0:8,Xmax,Ymax)
Integer i,x,y
Double precision uxy, uSqr
DO x = 1, Xmax
DO y = 1, Ymax
uSqr = u(x,y) * u(x,y) + v(x,y) * v(x,y)
DO i = 0, 8
uxy=u(x,y)*dr(i,0)+v(x,y)*dr(i,1)
fEq(i,x,y)= w(i)*Dxy(x,y)*(1.0d0+3.0d0*uxy+4.5d0*uxy*uxy-1.5d0*uSqr)
END DO
END DO
END DO
END SUBROUTINE computeFeq
!
!
! ****************************************************************
! * *
! * Fluid densities are rotated. By the next propagation step, *
! * this results in a bounce back from obstacle nodes. *
! * *
! ****************************************************************
!
SUBROUTINE Bounceback(wall,fIn,fOut)
USE simParam, ONLY: Xmax,Ymax
USE D2Q9Const, ONLY: opposite
IMPLICIT NONE
Integer, INTENT(IN) :: wall(Xmax,Ymax)
Double precision,INTENT(INOUT):: fIn(0:8,Xmax,Ymax),fOut(0:8,Xmax,Ymax)
INTEGER :: x,y,i
! loop over all f
DO x = 1, Xmax
DO y = 1, Ymax
IF (wall(x,y)==1) THEN
DO i = 0, 8
fOut(i,x,y) = fIn(opposite(i),x,y)
! f(1,x,y) = ftemp(3,x,y)
! f(2,x,y) = ftemp(4,x,y)
! f(3,x,y) = ftemp(1,x,y)
! f(4,x,y) = ftemp(2,x,y)
! f(5,x,y) = ftemp(7,x,y)
! f(6,x,y) = ftemp(8,x,y)
! f(7,x,y) = ftemp(5,x,y)
! f(8,x,y) = ftemp(6,x,y)
END DO
END IF
END DO
END DO
END SUBROUTINE Bounceback
!
!
! ****************************************************************
! * *
! * Propagate fluid densities to their next neighbour nodes *
! * *
! * *
! ****************************************************************
!
SUBROUTINE Stream(fOut,fIn)
USE simParam, ONLY: Xmax,Ymax
IMPLICIT NONE
Double precision :: fIn(0:8,Xmax,Ymax),fOut(0:8,Xmax,Ymax)
! local variables
INTEGER :: x,y,xr,xl,yt,yb
! loop over all f
DO x = 1, Xmax
DO y = 1, Ymax
yt = MOD(y,Ymax) + 1
xr = MOD(x,Xmax) + 1
yb = Ymax - MOD(Ymax + 1 - y, Ymax)
xl = Xmax - MOD(Xmax + 1 - x, Xmax)
fIn(0,x ,y ) = fOut(0,x,y)
fIn(1,xr,y )= fOut(1,x,y)
fIn(2,x,yt)= fOut(2,x ,y)
fIn(3,xl,y) = fOut(3,x,y)
fIn(4,x,yb) = fOut(4,x ,y)
fIn(5,xr,yt) = fOut(5,x,y)
fIn(6,xl,yt) = fOut(6,x,y)
fIn(7,xl,yb) = fOut(7,x,y)
fIn(8,xr,yb) = fOut(8,x,y)
END DO
END DO
END SUBROUTINE Stream
!
!
! ****************************************************************
! * *
! * One-step density Collision process *
! * *
! ****************************************************************
!
SUBROUTINE Collision(omega,fIn,fEq,wall,fOut)
USE simParam, ONLY: Xmax,Ymax
IMPLICIT NONE
Integer :: wall(Xmax,Ymax)
Double precision :: fEq(0:8,Xmax,Ymax),omega
Double precision :: fIn(0:8,Xmax,Ymax),fOut(0:8,Xmax,Ymax)
!
! local variables
INTEGER x,y,i
DO x = 1, Xmax
DO y = 1, Ymax
IF(wall(x,y)==0) THEN
DO i = 0, 8
fOut(i,x,y) = fIn(i,x,y)-omega*(fIn(i,x,y)-fEq(i,x,y))!(1.0-omega)f(i,x,y)+omegafEq(i,x,y)
END DO
END IF
END DO
END DO
! print *, fEq
END SUBROUTINE Collision
!
! ****************************************************************
! * *
! * Inlet & outlet velocity/pressure bounday condition
! * *
! ****************************************************************
! Non-equilibrium extrapolation method
! pressure condition
SUBROUTINE Inlet_outlet_BC(fEq,u,v,Dxy,fIn)
! SUBROUTINE Inlet_outlet_BC(f,u,v,Dxy,uSqr)
USE simParam, ONLY: Xmax,Ymax,Uw
USE D2Q9Const, ONLY: w,dr
IMPLICIT NONE
Double precision :: u(Xmax,Ymax),v(Xmax,Ymax),Dxy(Xmax,Ymax),fEq(0:8,Xmax,Ymax)
Double precision :: fIn(0:8,Xmax,Ymax)
Double precision :: fEqq(0:8,Xmax,Ymax)
Double precision :: density
Integer :: x, y, i
Double precision :: uxy, uSqr
! DO x = 1,Xmax,Xmax-1
! DO y = 1, Ymax
! Dxy(x,y) = 1.0d0
! u(x,y) = Uw
! v(x,y) = 0.0
! uSqr = u(x,y) * u(x,y) + v(x,y) * v(x,y)
!
! DO i = 0, 8
! uxy=u(x,y)dr(i,0)+v(x,y)dr(i,1)
! f(i,x,y)= w(i)Dxy(x,y)(1.0d0+3.0d0uxy+4.5d0uxyuxy-1.5d0uSqr)
! END DO
! END DO
! END DO
! PRINT*, uSqr
! only fluid nodes are considered here
!-----------------------------------------------------
! Inlet pressure conditon
!-----------------------------------------------------
DO y =1,Ymax
DO i = 0, 8
Dxy(1,y)= 1.1
u(1,y) = u(2,y)
v(1,y) = v(2,y)
uSqr = u(1,y) * u(1,y) + v(1,y) * v(1,y)
uxy=u(1,y)*dr(i,0)+v(1,y)*dr(i,1)
fEqq(i,1,y)= w(i)*Dxy(1,y)*(1.0d0+3.0d0*uxy+4.5d0*uxy*uxy-1.5d0*uSqr)
fIn(i,1,y) = fEqq(i,1,y)+(fIn(i,2,y)-fEq(i,2,y))
Dxy(Xmax,y)= 1.0000
u(Xmax,y) = u(Xmax-1,y)
v(Xmax,y) = v(Xmax-1,y)
uSqr = u(Xmax,y) * u(Xmax,y) + v(Xmax,y) * v(Xmax,y)
uxy=u(Xmax,y)*dr(i,0)+v(Xmax,y)*dr(i,1)
fEqq(i,Xmax,y)= w(i)*Dxy(Xmax,y)*(1.0d0+3.0d0*uxy+4.5d0*uxy*uxy-1.5d0*uSqr)
fIn(i,Xmax,y) = fEqq(i,Xmax,y)+(fIn(i,Xmax-1,y)-fEq(i,Xmax-1,y))
END DO
END DO
END SUBROUTINE Inlet_outlet_BC
! ****************************************************************
! * *
! * Output of rsults to file ‘T100XXX.dat’ *
! * *
! ****************************************************************
!
SUBROUTINE Dump(wall,Dxy,u,v,istep)
USE simParam, ONLY: Xmax,Ymax,x1,x2,y1,y2,a,b,NUMAX
IMPLICIT NONE
Integer :: wall(Xmax,Ymax)
Double precision :: u(Xmax,Ymax),v(Xmax,Ymax),Dxy(Xmax,Ymax)
Integer :: istep
! local variables
INTEGER x,y,i
REAL*8 t(Xmax,Ymax), p(Xmax,Ymax),cs2
CHARACTER *6 rr
! square speed of sound
cs2 = 1.d0 / 3.d0
WRITE (rr,200) 100000+istep/NUMAX
200 FORMAT (I6)
PRINT*,'istep= ', istep
! file T100XXX.dat for velocity, vorticity and pressure
open(11,file='T'//TRIM(rr)//'.dat')
!
! WRITE header for postprocessing with TECPLOT
!
WRITE(11,) 'VARIABLES = X, Y, U, V, W, P ’
WRITE(11,) ‘ZONE I=’, Xmax, ‘, J=’, Ymax, ‘, F=POINT’
!
DO y = 1, Ymax
DO x = 1, Xmax
! IF obstacle node, nothing is to DO
IF (wall(x,y)==1) THEN
! velocity components = 0
u(x,y) = 0.d0
v(x,y) = 0.d0
Dxy(x,y) = 0.d0
! pressure = average pressure
p(x,y) = 0
!
ELSE
! calculate local density Dxy
! Dxy=f(0,x,y)+f(1,x,y)+f(2,x,y)+f(3,x,y)+f(4,x,y)+f(5,x,y)+f(6,x,y)+f(7,x,y)+f(8,x,y)
! calculate velocity components
!
! U(x,y) = (f(1,x,y) + f(5,x,y) + f(8,x,y)-(f(3,x,y) + f(6,x,y) + f(7,x,y))) / Dxy
!
! V(x,y) = (f(2,x,y) + f(5,x,y) + f(6,x,y)-(f(4,x,y) + f(7,x,y) + f(8,x,y))) / Dxy
!
! calculate pressure
! P(x,y) = Dxy*cs
p(x,y)=Dxy(x,y)*cs2
!
END IF
END DO
END DO
!-------------------- Calculate vorticity function ----------------------
DO x=1,Xmax
DO y=1,Ymax
t(x,y)=0.0d0
ENDDO
ENDDO
DO x=2, Xmax-1
DO y=2, Ymax-1
IF (wall(x,y)==0) THEN
t(x,y)=(v(x+1,y)-v(x-1,y))/(Xmax)-(u(x,y+1)-u(x,y-1))/(Xmax)
ELSE
t(x,y)= 0.5*(v(x+1,y)-v(x-1,y))/(Xmax)-0.5*(u(x,y+1)-u(x,y-1))/(Xmax)
ENDIF
END DO
END DO
! ____________________________________________________________________
! WRITE results to file
DO y=1,Ymax
DO x=1,Xmax
WRITE(11,99) x,y,u(x,y),v(x,y),t(x,y),p(x,y)
99 format(2I8,4E15.6)
END DO
END DO
!
! ________________________________________________
WRITE(11,) ‘GEOMETRY M=GRID, FC=WHITE,F=POINT’
WRITE(11,) ‘2’
WRITE(11,) ‘4’
WRITE(11,) x1,y1
WRITE(11,) x1,y2
WRITE(11,) x2,y2
WRITE(11,) x2,y1
WRITE(11,) ‘2’
WRITE(11,) x1,y1
WRITE(11,) x2,y1
! _________________________________________________
CLOSE(11)
END SUBROUTINE Dump