Trouble with open boundary

gamlastan78 · June 24, 2009, 6:04pm

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

jlatt · June 24, 2009, 11:11pm

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:

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?
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.

gamlastan78 · June 30, 2009, 5:22pm

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