flow over a cylinder with LBM

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