# Trouble with open boundary

Hello,

I have some trouble with my open boundary while calculating poiseuille-flow in a 2-dimensional gap. The problem is that the resulting velocity profile on the outlet is kind of deformed and I have a loss of velocity in the x-direction. Unfortunately I can’t post pictures for a better discription.
To me it seems like the velocity of one of the diagonal directions ist wrong. But I guess that my mistake is somewhere in the application of the Zou/He-Boundary-Conditions. Maybe I made a mistake with the open boundary. I checked the propagation step a few time and think that it is implemented the right way.
I spent the whole day checking the code, but I can’t find my mistake. That’s why I posted the main parts of the code with the hope some of you have a few moments or maybe any idea what could be wrong.

Thanks a lot for any kind of hints.
Sven

! calculation of main values
nu = 0,024
nx = 251
nz =
Ma = 0.05
U_oo = 0.002
xi = U_oo / Ma ! Haenel Formel (9.50)
cs = xi / sqrt(3.d0)
cs_square = cs * cs
delta_t = dx / xi
gross_omega = (cs_square * delta_t) /(nu + cs_square*(delta_t/2.))

! Calculate velocity-profile at the inlet

``````    dp_dx = 0.0000001 * 100000 / lx
do z=1,nz_lb
if (z > node_wand_unten .and. z < node_wand_oben) then
u_inlet_2D(z) = h_inlet**2.d0 / (2.d0 * eta_phys) * dp_dx * ((z-j)*dz)/h_inlet *(1 - ((z-j)*dz)/h_inlet)
else
u_inlet_2D(z) = 0.d0
end if
end do
``````

! Calculate Inlet-Velocity

``````x = 1       ! left nodes at inlet
do z=1,nz_lb
if (spalt_2D(x,z) == 2) then
! Density Zou/He
rho_inlet = (fneq_2D(0,x,z) + fneq_2D(2,x,z) + fneq_2D(4,x,z) + 2 * (fneq_2D(3,x,z) + fneq_2D(6,x,z) + fneq_2D(7,x,z))) /(1-u_inlet_2D(z))

fneq_2D(1,x,z) = fneq_2D(3,x,z) + 4.0 * rho_inlet * u_inlet_2D(z) /6.d0
fneq_2D(5,x,z) = fneq_2D(7,x,z) - 0.5 *(fneq_2D(2,x,z) - fneq_2D(4,x,z)) + rho_inlet * u_inlet_2D(z) /(6.d0)
fneq_2D(8,x,z) = fneq_2D(6,x,z) + 0.5 *(fneq_2D(2,x,z) - fneq_2D(4,x,z)) + rho_inlet * u_inlet_2D(z) /(6.d0)
end if
end do
``````

! Bounce-Back

``````do x=1,nx_lb
do z=1,nz_lb
! spalt_2D(x,z) = 2 -> wall
if (spalt_2D(x,z) = 2.d0) then
fneq_2D(0,x,z) = fneq_2D(0,x,z)
fneq_2D(1,x,z) = fneq_2D(3,x,z)
fneq_2D(2,x,z) = fneq_2D(4,x,z)
fneq_2D(3,x,z) = fneq_2D(1,x,z)
fneq_2D(4,x,z) = fneq_2D(2,x,z)
fneq_2D(5,x,z) = fneq_2D(7,x,z)
fneq_2D(6,x,z) = fneq_2D(8,x,z)
fneq_2D(7,x,z) = fneq_2D(5,x,z)
fneq_2D(8,x,z) = fneq_2D(6,x,z)
end if
end do
end do
``````

! Calculate Density

``````do x=1,nx_lb
do z=1,nz_lb
rho_2D(x,z) = fneq_2D(0,x,z) + fneq_2D(1,x,z) + fneq_2D(2,x,z) + fneq_2D(3,x,z)  &
+ fneq_2D(4,x,z) + fneq_2D(5,x,z) + fneq_2D(6,x,z) + fneq_2D(7,x,z) + fneq_2D(8,x,z)
u_2D(x,z) = xi * (fneq_2D(1,x,z) + fneq_2D(5,x,z) + fneq_2D(8,x,z) - fneq_2D(3,x,z) - fneq_2D(6,x,z) - fneq_2D(7,x,z)) / rho_2D(x,z)
w_2D(x,z) = xi * (fneq_2D(2,x,z) + fneq_2D(5,x,z) + fneq_2D(6,x,z) - fneq_2D(4,x,z) - fneq_2D(7,x,z) - fneq_2D(8,x,z)) / rho_2D(x,z)
end do
end do
``````

! Calculate feq_2D

``````! x-direction
ex_2D(0) = 0
ex_2D(1) = 1
ex_2D(2) = 0
ex_2D(3) = -1
ex_2D(4) = 0
ex_2D(5) = 1
ex_2D(6) = -1
ex_2D(7) = -1
ex_2D(8) = 1
! z-direction
ez_2D(0) = 0
ez_2D(1) = 0
ez_2D(2) = 1
ez_2D(3) = 0
ez_2D(4) = -1
ez_2D(5) = 1
ez_2D(6) = 1
ez_2D(7) = -1
ez_2D(8) = -1

do x=1,nx
do z=1,nz
do i=0,8
! Calculate Velocity
vel_2D(x,z) = xi * (ex_2D(i) * u_2D(x,z) + ez_2D(i) * w_2D(x,z))
vel_2D_sqr(x,z) = u_2D(x,z) * u_2D(x,z) + w_2D(x,z) * w_2D(x,z)

feq_2D(i,x,z) = rho_2D(x,z) * tp_2D(i) * ((rho_2D(x,z) * cs**2.d0)/(rho_2D(x,z) * cs**2.d0) &
+ (1.d0/cs**2.d0) * vel_2D(x,z) + 1.d0/(2.d0 * cs**4.d0) * vel_2D(x,z)**2.d0 &
- 1.d0/(2.d0 * cs**2.d0) * vel_2D_sqr(x,z))
end do
end do
end do
``````

! Collision-Step

``````do x=1,nx_lb
do z=1,nz_lb
do i=0,8
fneq_2D(i,x,z) = fneq_2D(i,x,z) + gross_omega * (feq_2D(i,x,z) - fneq_2D(i,x,z))
end do
end do
end do
``````

! Propagation-Step

``````do z=1,nz_lb
do x=1,nx_lb-1
fneq_2D(1,x+1,z) = feq_2D(1,x,z)       ! east
end do
do x=nx_lb,2,-1
fneq_2D(3,x-1,z) = feq_2D(3,x,z)       ! west
end do
end do

do x=1,nx_lb
do z=1,nz_lb-1
fneq_2D(2,x,z+1) = feq_2D(2,x,z)       ! north
end do
do z=nz_lb,2,-1
fneq_2D(4,x,z-1) = feq_2D(4,x,z)       ! south
end do
end do

do z=1,nz_lb-1
do x=1,nx_lb-1
fneq_2D(5,x+1,z+1) = feq_2D(5,x,z)       ! north-east
end do
do x=nx_lb,2,-1
fneq_2D(6,x-1,z+1) = feq_2D(6,x,z)       ! north-west
end do
end do

do z=nz_lb,2,-1
do x=1,nx_lb-1
fneq_2D(8,x+1,z-1) = feq_2D(8,x,z)       ! south-east
end do
do x=nx_lb,2,-1
fneq_2D(7,x-1,z-1) = feq_2D(7,x,z)       ! south-west
end do
end do``````

Hello,

[Note: there is a button in the editor which you can use to format codes more beautifully]

After skimming through your code, the following two points come two my mind:

1. There’s a variable called fneq, which could be thought to correspond to non-equilibrium parts of the populations, given its name. But that’s not true: it stands simply for the plain populations, right?

2. The propagation step looks suspicious enough to me. For example:

``````
fneq_2D(1,x+1,z) = feq_2D(1,x,z)

``````

This means that you are just propagating the value of the equilibrium, and override the previously computed values “fneq” of the post-collision populations. I’d rather expect something like this:

``````
tmp(1,x+1,z) = fneq_2D(1,x,z)

``````

and after the loop, copy tmp to fneq (but again, all of this is only valid if fneq, in spite of its name, corresponds to f).

A final comment: you’ll find simple example codes in various programming languages on lbmethod.org, and some of them include Zou/He boundary conditions. It might be a good idea to take one of these as an example to get a first working LB program.

Hello jlatt,

thanks a lot for your hints. It helped me to improve the code. You’re absolutely right, fneq = f. I already changed it. And honestly I tried to copy one sample code and were sure I made it right, but I didn’t copy f to a temporary matrix, because I thought it doesn’t matter.
After changing that the results became much better and the velocity profile is symmetric.

Unfortunately I have another problem now: the velocity profile becomes higher at the outlet than it is at the inlet. I already changed from open boundary conditions to periodic boundary conditions, but then I become a fall-off of velocity close to the outlet. I don’t htink that this is due to the propagation step, since it is a stationary problem which doesn’t move through the gap.

Again, I would be very grateful, if somebody have a hint.
Sven