# Doubts in Zou-He boundary condtion

Hello,
I am very new to the world of LBM. I am trying to simulate the benchmark problems with different boundary conditions.
I am confused with the algorithm of the Zou-He boundary condition. I want to use inlet velocity and zerogradient at outlet BC. Could you have a look at the following code and correct me where am I making mistake. Thank you

// 8 3 5 ^y
// |/ | x
// 2-0-1 —>
// /|
// 6 4 7

void LBM(int time) {

``````/// Swap populations
// The present code used old and new populations which are swapped at the beginning of each time step.
// This is sometimes called 'double-buffered' or 'ping-pong' algorithm.
// This way, the old populations are not overwritten during propagation.
// The resulting code is easier to write and to debug.
// The memory requirement for the populations is twice as large.

double ***swap_temp = pop_old;
pop_old = pop;
pop = swap_temp;

/// Lattice Boltzmann equation
// The lattice Boltzmann equation is solved in the following.
// The algorithm includes
// - computation of the lattice force
// - combined collision and propagation (faster than first collision and then propagation)

for (int X = 1; X < Nx-2; ++X) {
for (int Y = 1; Y < Ny - 1; ++Y) {

/// Compute equilibrium
// The equilibrium populations are computed.

equilibrium(density[X][Y], (velocity_x[X][Y] + (force_x[X][Y] + gravity) * tau / density[X][Y]), (velocity_y[X][Y] + (force_y[X][Y]) * tau / density[X][Y]));

/// Compute new populations
// This is the lattice Boltzmann equation (combined collision and propagation) including external forcing.

pop[X][Y] = pop_old[X][Y] * (1 - omega) + pop_eq * omega;
pop[(X + 1)][Y] = pop_old[X][Y] * (1 - omega) + pop_eq * omega;
pop[X - 1][Y] = pop_old[X][Y] * (1 - omega) + pop_eq * omega;
pop[X][Y + 1] = pop_old[X][Y] * (1 - omega) + pop_eq * omega;
pop[X][Y - 1] = pop_old[X][Y] * (1 - omega) + pop_eq * omega;
pop[X + 1][Y + 1] = pop_old[X][Y] * (1 - omega) + pop_eq * omega;
pop[X - 1][Y - 1] = pop_old[X][Y] * (1 - omega) + pop_eq * omega;
pop[X + 1][Y - 1] = pop_old[X][Y] * (1 - omega) + pop_eq * omega;
pop[X - 1][Y + 1] = pop_old[X][Y] * (1 - omega) + pop_eq * omega;

}
``````

//Bounce back condition for top and bottom wall
for (int X = 1; X < Nx-1; X++)
{

``````// Bottom wall (y = 0)

pop[X] = pop[X];
pop[X] = pop[X];
pop[X] = pop[X];

//Top wall (y = Ny - 1)

pop[X][Ny - 2] = pop[X][Ny - 1];
pop[X][Ny - 2] = pop[X][Ny - 1];
pop[X][Ny - 2] = pop[X][Ny - 1];
``````

}
//Zou-He BC for inlet and outlet

for (int Y = 1; Y < Ny-1; Y++) {

``````//Left wall (X = 0)

pop[Y] = pop[Y] + (2 / 3)* inletdensity * inletvel;
pop[Y] = pop[Y] - 0.5*(pop[Y] - pop[Y]) + (1 / 6)*inletdensity* inletvel;
pop[Y] = pop[Y] + 0.5*(pop[Y] - pop[Y]) + (1 / 6)*inletdensity * inletvel;

//Right wall (X = Nx - 1)

pop[Nx - 1][Y] = pop[Nx - 2][Y];
pop[Nx - 1][Y] = pop[Nx - 2][Y];
pop[Nx - 1][Y] = pop[Nx - 2][Y];
``````

}