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
```