Hi all.
I run flow over a cylinder with lattice boltzmann method. I want to implement pressure boundry condition . I need a paper to valid my result with it. can any one help me to find it?
I tried the 2d cylinder flow code in fortran provided from http://wiki.palabos.org/numerics:codes
and plot it with tecplot.
Unfortunately the plot is wrong since the streamlines pass through the cylinder, instead of circumventing it and creating vortex behind.
even after I fill the cylinder with “0” value, the streamlines are still wrong,eg, the velocity on the cylinder is not “0”.
I compared the fortran code with the matlab one but cant figure out the reason.The boundary treatment seems reasonable and the way to handle obstacle is quite general.
Anyone can simply copy the code below and run it.Then it will generate the tecplot format output.
Hope someone can point out the bug.
[code="fortran"]
! unsteady.f90:
! This example examines an unsteady flow past a cylinder placed in a channel.
! The cylinder is offset somewhat from the center of the flow to make the
! steady-state symmetrical flow unstable. At the inlet and outlet, a Poiseuille
! profile is imposed on the velocity. At Reynolds numbers around 100,
! an unstable periodic pattern is created, the Karman vortex street.
! Note that with the implemented Zou/He boundary condition, you must
! increase the resolution to keep the simulation stable if you increase
! the Reynolds number.
! ========================================================
! Constants that identify different cell-types according
! to the dynamics they implement
! ========================================================
MODULE cellConst
integer, parameter:: fluid = 0, wall = 100, inlet = 50, outlet = 80
END MODULE cellConst
! ========================================================
! Lattice constants for the D2Q9 lattice
! ========================================================
MODULE D2Q9Const
! D2Q9 Weights
double precision,parameter:: t(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:: v(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:: xDim = 520
integer, parameter:: yDim = 180
integer, parameter:: obstX = xDim/4
integer, parameter:: obstY = yDim/2
integer, parameter:: obstR = yDim/9
integer, parameter:: tMax = 200000
double precision, parameter:: uMax = 0.40d0
double precision, parameter:: Re = 20.0d0
END MODULE simParam
! ========================================================
! The main program, implementing a flow past a cylinder
! ========================================================
PROGRAM unsteady
USE simParam, ONLY: xDim, yDim, tMax
use cellConst
implicit none
double precision:: omega, time1, time2, timeTot
double precision, dimension(:,:,:), allocatable:: f,f_previous, fEq, u,u0
double precision, dimension(:,:), allocatable:: rho, uSqr, lift, drag
integer, dimension(:,:), allocatable:: image
integer:: tStep
integer:: i, j,k ,x, y
double precision :: Err
allocate(f(yDim,xDim,0:8))
allocate(f_previous(yDim,xDim,0:8))
allocate(fEq(yDim,xDim,0:8))
allocate(u(yDim,xDim,0:1))
allocate(u0(yDim,xDim,0:1))
allocate(uSqr(yDim,xDim))
allocate(rho(yDim,xDim))
allocate(lift(yDim,xDim))
allocate(drag(yDim,xDim))
allocate(image(yDim,xDim))
CALL constructImage(image)
CALL computeOmega(omega)
CALL writeInput(omega)
CALL initMacro(rho,u,u0,uSqr,image)
CALL computeFeq(fEq,rho,u,uSqr)
f = fEq
timeTot = 0.0d0
do tStep = 1, tMax
CALL CPU_TIME(time1)
CALL collide(f,fEq,omega,image)
CALL stream(f)
CALL inletOutlet(f,rho,u,image)
CALL boundaries(f,image,u)
CALL computeMacros(f,rho,u,uSqr)
CALL computeFeq(fEq,rho,u,uSqr)
!f_previous(1:yDim,1:xDim,0:8)=f(1:yDim,1:xDim,0:8)
!call computeForce(f,rho,u,uSqr,f_previous,lift,drag)
CALL CPU_TIME(time2)
timeTot = timeTot + (time2-time1)
if( mod(tStep,1000) == 0) then
call output(tStep,image,u,rho,lift,drag)
write(*,*) "tStep ",tStep," done "," err ",Err(u,u0,image)
end if
if( Err(u,u0,image) < 1.0d-6 .or. isnan(Err(u,u0,image))) then
exit
end if
end do
CALL writeImage(image)
!CALL writeOutput(u,0)
write(*,*) dble(tMax) * (dble(yDim * xDim)) / timeTot ,'cells per second'
deallocate(f)
deallocate(fEq)
deallocate(u)
deallocate(uSqr)
deallocate(rho)
deallocate(image)
END PROGRAM unsteady
double precision function Err(u,u0,image)
use cellConst
USE simParam
double precision, INTENT(INOUT):: u(yDim,xDim,0:1)
double precision, INTENT(INOUT):: u0(yDim,xDim,0:1)
integer, INTENT(INOUT)::image(yDim, xDim)
integer:: i, j,k ,x, y
double precision :: e1,e2
do i=1,yDim
do j=1,xDim
if(image(i,j)/=wall) then
e1=e1+ sqrt( ( u(i,j,0)-u0(i,j,0) )**2+ ( u(i,j,1)-u0(i,j,1) )**2 )
!e2+=sqrt( ux[j][i]*ux[j][i]+uy[j][i]*uy[j][i]);
e2=e2+sqrt( u(i,j,0) **2+ u(i,j,1)**2)
!u0[j][i]=ux[j][i];v0[j][i]=uy[j][i];
u0(i,j,0)=u(i,j,0)
u0(i,j,1)=u(i,j,1)
end if
end do
end do
Err =e1/e2
end function Err
subroutine output(tStep,image,u,rho,lift,drag)
use cellConst
USE simParam, ONLY: xDIm, yDim
double precision, INTENT(INOUT):: u(yDim,xDim,0:1), rho(yDim, xDim),lift(yDim, xDim),drag(yDim, xDim)
integer, INTENT(INOUT)::image(yDim, xDim)
integer:: i,j
integer:: tStep
Character( Len = 1200 ) :: filename
write(filename,"(I7,A4)") tStep,".plt"
open(1,file=filename)
write(1,*)"VARIABLES =X, Y, U, V, RHO, IMAGE, Vel"
write(1,*)"ZONE ","I=",xDim,"J=",yDim,",","F=POINT"
do i=1,yDim
do j=1,xDim
if(image(i,j)==wall) then
u(i,j,0)=0
u(i,j,1)=0
rho(i,j)=wall
!! lift(i,j)=0
end if
write(1,*) j, i,u(i,j,0), u(i,j,1), rho(i,j),image(i,j),sqrt(u(i,j,0)**2+u(i,j,1)**2)
end do
end do
end subroutine output
! ========================================================
! Compute the relaxation parameter from the Reynolds number
! ========================================================
SUBROUTINE computeOmega(omega)
USE simParam, ONLY: Re,uMax,obstR
implicit none
double precision, INTENT(INOUT):: omega
double precision:: nu
nu = uMax * 2.0d0 * dble(obstR) / Re
omega = 1.0d0 / (3.0d0*nu+0.5d0)
END SUBROUTINE computeOmega
! ========================================================
! Construct an array the defines the flow geometry
! ========================================================
SUBROUTINE constructImage(image)
USE cellConst
USE simParam, ONLY: xDim, yDim, obstX, obstY, obstR
USE D2Q9Const, ONLY: v
implicit none
integer, INTENT(INOUT):: image(yDim,xDim)
integer:: x,y
v(0:8,0) = (/0,1,0,-1,0,1,-1,-1,1/)
v(0:8,1) = (/0,0,1,0,-1,1,1,-1,-1/)
image = fluid
image(:,1) = inlet
image(:,xDim) = outlet
image(1,:) = wall
image(yDim,:) = wall
do x = 1, xDim
do y = 1, yDim
if (((x-obstX)**2 + (y-obstY)**2) <= (obstR**2) ) image(y,x) = wall
end do
end do
END SUBROUTINE constructImage
! ========================================================
! Initialize the simulation to Poiseuille profile at
! an equilibrium distribution
! ========================================================
SUBROUTINE initMacro(rho,u,u0,uSqr,image)
USE simParam!, ONLY: xDim, yDim
use cellConst
implicit none
integer, INTENT(INOUT):: image(yDim,xDim)
double precision, INTENT(INOUT):: rho(yDim,xDim), u(yDim,xDim,0:1), u0(yDim,xDim,0:1),uSqr(yDim,xDim)
double precision:: uProf
integer:: x,y
do y = 1, yDim
!u(y,:,0) = uProf(y)
u(y,:,0) = 0.0d0
u(y,:,1) = 0.0d0
u0(y,:,0) = 0.0d0
u0(y,:,1) = 0.0d0
end do
do x=1,xDim
do y=1,yDim
if(image(y,x)/=wall) u(y,:,0) = 0.4d0!uProf(y)
end do
end do
rho = 1.0d0
uSqr = u(:,:,0) * u(:,:,0) + u(:,:,1) * u(:,:,1)
END SUBROUTINE initMacro
! ========================================================
! Compute equilibrium distribution
! ========================================================
SUBROUTINE computeFeq(fEq,rho,u,uSqr)
USE D2Q9COnst, ONLY: t, v
USE simParam, ONLY: xDim, yDim
implicit none
double precision, INTENT(IN):: rho(yDim,xDim), uSqr(yDim,xDim), u(yDim,xDim,0:1)
double precision, INTENT(INOUT):: fEq(yDim,xDim,0:8)
integer:: i, x, y
double precision:: uxy
do i = 0, 8
do x = 1, xDim
do y = 1, yDim
uxy = u(y,x,0) * v(i,0) + u(y,x,1) * v(i,1)
fEq(y,x,i) = t(i) * rho(y,x) * (1.0d0 + 3.0d0 * uxy + 4.5d0 * uxy * uxy - 1.5d0 * uSqr(y,x))
end do
end do
end do
END SUBROUTINE computeFeq
! ========================================================
! Compute drag and lift from distribution functions
! ========================================================
SUBROUTINE computeForce(f,rho,u,uSqr,f_previous,lift,drag)
USE simParam, ONLY: xDIm, yDim
implicit none
double precision, INTENT(IN):: f(yDim,xDim,0:8),f_previous(yDim,xDim,0:8)
double precision, INTENT(INOUT):: u(yDim,xDim,0:1), rho(yDim, xDim), uSqr(yDim, xDim), lift(yDim, xDim),drag(yDim, xDim)
integer:: x,y
do x = 1, xDim-1
do y = 1, yDim-1
!rho(y,x) = f(y,x,0) + f(y,x,1) + f(y,x,2) + f(y,x,3) + f(y,x,4) + f(y,x,5) + f(y,x,6) + f(y,x,7) + f(y,x,8)
!u(y,x,0) = (f(y,x,1) - f(y,x,3) + f(y,x,5) - f(y,x,6) - f(y,x,7) + f(y,x,8)) / rho(y,x)
!u(y,x,1) = (f(y,x,2) - f(y,x,4) + f(y,x,5) + f(y,x,6) - f(y,x,7) - f(y,x,8)) / rho(y,x)
!uSqr(y,x) = u(y,x,0) * u(y,x,0) + u(y,x,1) * u(y,x,1)
drag(y,x) = (f(y,x+1,1) - f(y,x+1,3) + f(y,x+1,5) - f(y,x+1,6) - f(y,x+1,7) + f(y,x+1,8)) + (f_previous(y,x,1) - f_previous(y,x,3) + f_previous(y,x,5) - f_previous(y,x,6) - f_previous(y,x,7) + f_previous(y,x,8))
lift(y,x) = (f(y+1,x,2) - f(y+1,x,4) + f(y+1,x,5) + f(y+1,x,6) - f(y+1,x,7) - f(y+1,x,8)) + (f_previous(y,x,2) - f_previous(y,x,4) + f_previous(y,x,5) + f_previous(y,x,6) - f_previous(y,x,7) - f_previous(y,x,8))
end do
end do
END SUBROUTINE computeForce
! ========================================================
! Compute density and velocity from distribution functions
! ========================================================
SUBROUTINE computeMacros(f,rho,u,uSqr)
USE simParam, ONLY: xDIm, yDim
implicit none
double precision, INTENT(IN):: f(yDim,xDim,0:8)
double precision, INTENT(INOUT):: u(yDim,xDim,0:1), rho(yDim, xDim), uSqr(yDim, xDim)
integer:: x,y
do x = 1, xDim
do y = 1, yDim
!if (image(y,x) /= wall) then
rho(y,x) = f(y,x,0) + f(y,x,1) + f(y,x,2) + f(y,x,3) + f(y,x,4) + f(y,x,5) + f(y,x,6) + f(y,x,7) + f(y,x,8)
u(y,x,0) = (f(y,x,1) - f(y,x,3) + f(y,x,5) - f(y,x,6) - f(y,x,7) + f(y,x,8)) / rho(y,x)
u(y,x,1) = (f(y,x,2) - f(y,x,4) + f(y,x,5) + f(y,x,6) - f(y,x,7) - f(y,x,8)) / rho(y,x)
uSqr(y,x) = u(y,x,0) * u(y,x,0) + u(y,x,1) * u(y,x,1)
!end if
end do
end do
END SUBROUTINE computeMacros
! ========================================================
! Implement Bounce-back on upper/lower boundaries
! ========================================================
SUBROUTINE boundaries(f,image,u)
USE D2Q9Const, ONLY: opposite
USE cellConst, ONLY: wall
USE simParam, ONLY: xDim, yDim
implicit none
integer, INTENT(IN):: image(yDim,xDim)
double precision, INTENT(INOUT):: f(yDim,xDim,0:8), u(yDim,xDim,0:1)
!double precision, INTENT(INOUT):: rho(yDim,xDim)
double precision:: fTmp(0:8)
integer:: i, x, y
do x = 1, xDim
do y = 1, yDim
if (image(y,x) == wall) then
u(y,x,0)=0.0
u(y,x,1)=0.0
do i = 0, 8
fTmp(i) = f(y,x,opposite(i))
end do
do i = 0, 8
f(y,x,i) = fTmp(i)
end do
!rho(y,x)=wall
end if
end do
end do
END SUBROUTINE boundaries
! ========================================================
! Use Zou/He boundary condition to implement Dirichlet
! boundaries on inlet/outlet
! ========================================================
SUBROUTINE inletOutlet(f,rho,u,image)
USE cellConst, ONLY: inlet, outlet
USE simParam
implicit none
double precision, INTENT(INOUT):: f(yDim,xDim,0:8), u(yDim,xDim,0:1), rho(yDim,xDim)
integer, INTENT(IN):: image(yDim,xDim)
double precision:: uProf
integer:: x, y
do x = 1, xDim
do y = 1, yDim
if (image(y,x) == inlet) then
u(y,x,0) = uProf(y)
u(y,x,1) = 0.0d0
CALL inletZou(f(y,x,:),u(y,x,:),rho(y,x))
else if (image(y,x) == outlet) then
u(y,x,0) = uProf(y)
u(y,x,1) = 0.0d0
CALL outletZou(f(y,x,:),u(y,x,:),rho(y,x))
end if
end do
end do
CONTAINS
! ========================================================
! Zou/He boundary on inlet
! ========================================================
SUBROUTINE inletZou(f,u,rho)
implicit none
double precision, INTENT(INOUT):: f(0:8),rho
double precision, INTENT(IN):: u(0:1)
double precision:: fInt, fInt2
fInt = f(0) + f(2) + f(4)
fInt2 = f(3) + f(6) + f(7)
rho = (fInt + 2.0d0 * fInt2) / (1.0d0 - u(0))
CALL zouWestWall(f,rho,u)
END SUBROUTINE inletZou
SUBROUTINE zouWestWall(f,rho,u)
implicit none
double precision, INTENT(INOUT):: f(0:8)
double precision, INTENT(IN):: rho, u(0:1)
double precision:: fDiff, rhoUx, rhoUy
fDiff = 0.5d0 * (f(2) - f(4))
rhoUx = rho * u(0) / 6.0d0
rhoUy = 0.5d0 * rho * u(1)
f(1) = f(3) + 4.0d0 * rhoUx
f(5) = f(7) - fDiff + rhoUx + rhoUy
f(8) = f(6) + fDiff + rhoUx - rhoUy
END SUBROUTINE zouWestWall
! ========================================================
! Zou/He boundary on outlet
! ========================================================
SUBROUTINE outletZou(f,u,rho)
implicit none
double precision, INTENT(INOUT):: f(0:8),rho,u(0:1)
double precision:: fInt, fInt2
fInt = f(0) + f(2) + f(4)
fInt2 = f(1) + f(8) + f(5)
rho = (fInt + 2.0d0 * fInt2) / (1.0d0 + u(0))
CALL zouEastWall(f,rho,u)
END SUBROUTINE outletZou
SUBROUTINE zouEastWall(f,rho,u)
implicit none
double precision, INTENT(INOUT):: f(0:8)
double precision, INTENT(IN):: rho, u(0:1)
double precision:: fDiff, rhoUx, rhoUy
fDiff = 0.5d0 * (f(2) - f(4))
rhoUx = rho * u(0) / 6.0d0
rhoUy = 0.5d0 * rho * u(1)
f(3) = f(1) - 4.0d0 * rhoUx
f(7) = f(5) + fDiff - rhoUx - rhoUy
f(6) = f(8) - fDiff - rhoUx + rhoUy
END SUBROUTINE zouEastWall
END SUBROUTINE inletOutlet
! ========================================================
! Computation of Poiseuille profile for the inlet/outlet
! ========================================================
FUNCTION uProf(y)
USE simParam, ONLY: yDIm, uMax
implicit none
integer, INTENT(IN):: y
double precision:: radius, uProf
radius = dble(yDim-1) * 0.5d0
uProf = -uMax * ((abs(1 - dble(y-1) / radius))**2 - 1.0d0)
END FUNCTION uProf
! ========================================================
! Streaming step: the population functions are shifted
! one site along their corresponding lattice direction
! (no temporary memory is needed)
! ========================================================
SUBROUTINE stream(f)
USE simParam
implicit none
double precision, INTENT(INOUT):: f(yDim,xDim,0:8)
double precision:: periodicHor(yDim), periodicVert(xDim)
! -------------------------------------
! right direction
periodicHor = f(:,xDim,1)
f(:,2:xDim,1) = f(:,1:xDim-1,1)
f(:,1,1) = periodicHor
! -------------------------------------
! up direction
periodicVert = f(1,:,2)
f(1:yDim-1,:,2) = f(2:yDim,:,2)
f(yDim,:,2) = periodicVert
! -------------------------------------
! left direction
periodicHor = f(:,1,3)
f(:,1:xDim-1,3) = f(:,2:xDim,3)
f(:,xDim,3) = periodicHor
! -------------------------------------
! down direction
periodicVert = f(yDim,:,4)
f(2:yDim,:,4) = f(1:yDim-1,:,4)
f(1,:,4) = periodicVert
! -------------------------------------
! up-right direction
periodicVert = f(1,:,5)
periodicHor = f(:,xDim,5)
f(1:yDim-1,2:xDim,5) = f(2:yDim,1:xDim-1,5)
f(yDim,2:xDim,5) = periodicVert(1:xDim-1)
f(yDim,1,5) = periodicVert(xDim)
f(1:yDim-1,1,5) = periodicHor(2:yDim)
! -------------------------------------
! up-left direction
periodicVert = f(1,:,6)
periodicHor = f(:,1,6)
f(1:yDim-1,1:xDim-1,6) = f(2:yDim,2:xDim,6)
f(yDim,1:xDim-1,6) = periodicVert(2:xDim)
f(yDim,xDim,6) = periodicVert(1)
f(1:yDim-1,xDim,6) = periodicHor(2:yDim)
! -------------------------------------
! down-left direction
periodicVert = f(yDim,:,7)
periodicHor = f(:,1,7)
f(2:yDim,1:xDim-1,7) = f(1:yDim-1,2:xDim,7)
f(1,1:xDim-1,7) = periodicVert(2:xDim)
f(1,xDim,7) = periodicVert(1)
f(2:yDim,xDim,7) = periodicHor(1:yDim-1)
! -------------------------------------
! down-right direction
periodicVert = f(yDim,:,8)
periodicHor = f(:,xDim,8)
f(2:yDim,2:xDim,8) = f(1:yDim-1,1:xDim-1,8)
f(1,2:xDim,8) = periodicVert(1:xDim-1)
f(1,1,8) = periodicVert(xDim)
f(2:yDim,1,8) = periodicHor(1:yDim-1)
END SUBROUTINE stream
! ========================================================
! LBGK collision step
! ========================================================
SUBROUTINE collide(f,fEq,omega,image)
USE simParam, ONLY: xDim, yDim
USE cellConst, ONLY: wall
implicit none
integer, INTENT(IN):: image(yDim,xDim)
double precision, INTENT(IN):: fEq(yDim,xDim,0:8), omega
double precision, INTENT(INOUT):: f(yDim,xDim,0:8)
integer:: x,y,i
do i = 0, 8
do x = 1, xDim
do y = 1, yDim
!if (image(y,x) /= wall)
f(y,x,i) = (1.0d0 - omega) * f(y,x,i) + omega * feq(y,x,i)
end do
end do
end do
END SUBROUTINE collide
! ========================================================
! Write the flow geometry to a file
! ========================================================
SUBROUTINE writeImage(image)
USE simParam, ONLY: xDim, yDim
implicit none
integer, INTENT(IN):: image(yDim,xDim)
integer:: x,y
open(13,file='outputImage.dat')
do x=1, xDim
do y=1, yDim
write(13,102) image(y,x)
end do
end do
102 format(3i10)
close(15)
END SUBROUTINE writeImage
! ========================================================
! Print out simulation parameters to screen
! ========================================================
SUBROUTINE writeInput(omega)
USE simParam
implicit none
double precision, INTENT(IN):: omega
write(*,*) 'xDim = ', xDim
write(*,*) 'yDim = ', yDim
write(*,*) 'Obstacle X = ', obstX
write(*,*) 'Obstacle Y = ', obstY
write(*,*) 'Obstacle Radius = ', obstR
write(*,*) 'tMax = ', tMax
write(*,*) 'uMax = ', uMax
write(*,*) 'Re = ', Re
write(*,*) 'omega = ', omega
END SUBROUTINE writeInput