Flow around circular cylinder

Lattice Boltzmann Method LBM: Algorithms
#1

Hi all

I’m trying to simulate 2D flow around circular cylinder in channel.
I am using Fortran for coding.
The algorithm summarized as follows:

1-Initialize density distribution function with equlibrum distribution
D=Diameter
RHO0LB=1.0
VELO=0.03
Nu=(VELOD)/Re
TAU=3.0Nu+0.5
NI=50D
NJ=8D

2- HYDROVAR: In this routine macroscopice properties (namly rho,u, and v) are computed.

3-COLLISION: In this routine local collision are performed in the flow domain.(COMPUTE EQ. DISTRIBUTION FUNCTION FOR NODES)

4-PROPAGATION:In this routine interior distribution functions from non-occupied nodes alongs the lattice connections lines to their next neighbours.(Bouzidi et.al interpolation-based models for curved boundary treatment)

5-Boundary Condition:

[code=“fortran”]
!%%%%%%%%%%%%%%%%%%%
!..LOWER WALL
!%%%%%%%%%%%%%%%%%%%
!
DO I=1,NI
!
FIN(2,I,2) = FOUT(4,I,2)

  FIN(5,I,2) = FOUT(7,I,2)
  
  FIN(6,I,2) = FOUT(8,I,2)

  ENDDO

!
!%%%%%%%%%%%%%%%%%%%
!..UPPER WALL
!%%%%%%%%%%%%%%%%%%%
!
DO I=1,NI
!
FIN(4,I,NJM) = FOUT(2,I,NJM)

  FIN(7,I,NJM) = FOUT(5,I,NJM)

  FIN(8,I,NJM) = FOUT(6,I,NJM)

  
  ENDDO

!
!
!%%%%%%%%%%%%%%%%%%%
!..Inlet
!%%%%%%%%%%%%%%%%%%%
!
DO J=2,NJM

!
U(1,J)=VELO !(VELO =0.03)
V(1,J)=0.0
!
RHO(1,J)=(FIN(0,1,J)+FIN(2,1,J)+FIN(4,1,J)+2.*(FIN(3,1,J)+FIN(6,1,J)+FIN(7,1,J)))/(1-U(1,J))
!
FIN(1,1,J) = FIN(3,1,J) + (2.0D0 / 3.0D0 * (RHO(1,J)*U(1,J)))
!
FIN(5,1,J) = FIN(7,1,J) + (RHO(1,J)*U(1,J) / 6.0D0)
!
FIN(8,1,J) = FIN(6,1,J) + (RHO(1,J)*U(1,J) / 6.0D0)

  ENDDO

!
!%%%%%%%%%%%%%%%%%%%
! …Outlet
!%%%%%%%%%%%%%%%%%%%
!
DO J=2,NJM

 V(NI,J)=0.0
 
 FIN(1,NI,J) = 1.5 * FIN(1,NIM,J) - 0.5 * FIN(1,NIM-1,J)
 FIN(5,NI,J) = 1.5 * FIN(5,NIM,J) - 0.5 * FIN(5,NIM-1,J)
 FIN(8,NI,J) = 1.5 * FIN(8,NIM,J) - 0.5 * FIN(8,NIM-1,J) 

  ENDDO   	      




6- Calculate Drag Coeff

[code="fortran"]
 DO J=2,NJM
 DO I=2,NIM
 IF (.NOT.OBST1(I+1,J)) THEN
 ELSE IF (.NOT.OBST1(I,J)) THEN
 DRAG_FORCE=((CX(3)*(FIN(3,I,J))-CX(1)*(FOUT(1,I,J))))+DRAG_FORCE
 END IF

 IF (.NOT.OBST1(I-1,J)) THEN
 ELSE IF (.NOT.OBST1(I,J)) THEN
 DRAG_FORCE=((CX(1)*(FIN(1,I,J))-CX(3)*(FOUT(3,I,J))))+DRAG_FORCE
 ENDIF
 
IF (.NOT.OBST1(I+1,J+1)) THEN
 ELSE IF (.NOT.OBST1(I,J)) THEN
 DRAG_FORCE=((CX(7)*(FIN(7,I,J))-CX(5)*(FOUT(5,I,J))))+DRAG_FORCE
 ENDIF

 IF (.NOT.OBST1(I-1,J+1)) THEN
 ELSE IF (.NOT.OBST1(I,J)) THEN
 DRAG_FORCE=((CX(8)*(FIN(8,I,J))-CX(6)*(FOUT(6,I,J))))+DRAG_FORCE
 ENDIF
 
IF (.NOT.OBST1(I-1,J-1)) THEN
 ELSE IF (.NOT.OBST1(I,J)) THEN
 DRAG_FORCE=((CX(5)*(FIN(5,I,J))-CX(7)*(FOUT(7,I,J))))+DRAG_FORCE
 ENDIF

 IF (.NOT.OBST1(I+1,J-1)) THEN
 ELSE IF (.NOT.OBST1(I,J)) THEN
 DRAG_FORCE=((CX(6)*(FIN(6,I,J))-CX(8)*(FOUT(8,I,J))))+DRAG_FORCE
 ENDIF
  
 END DO
 END DO
 
 CD=((DRAG_FORCE))/(0.5*RHO0LB*(VELO**2)*D)

In calculation of Cd from D=10 , value of Cd is equal to 1.85 for Re=100 but according to references Cd is equal1.3 …

which part is wrong?
Is my definition of Boundary Conditions wrong?

Thanks and regards