I’ve developed a Lattice Boltzmann code utilizing the D3Q19 model to tackle fluid dynamics within complex geometries. For the inlet boundary condition, I’ve implemented the ZouHe method as outlined in the publication https://iopscience.iop.org/article/10.1088/1742-5468/2010/01/P01018. I’ve enforced a fully bounce-back condition for the walls and maintained a constant pressure at the outlet. However, I’ve encountered two issues with obtaining the desired results. Firstly, the ZouHe method isn’t setting the velocity as anticipated. Secondly, while I expect velocity components at the wall to be precisely zero, they instead exhibit small values. On wall also the density value is around 1.03 instead of being exactly 1 and this is suggesting there should be something wrong with the distribution function.
I’ll share a portion of my code and welcome any insights to resolve these issues. I can’t generate a reproduceable code but I can share any function if required.

void Simulation::Collision()
{
    for (int k = 0; k < SubDomain_.my_Nz_; k++) {
        for (int j = 0; j < SubDomain_.my_Ny_; j++) {
            for (int i = 0; i < SubDomain_.my_Nx_; i++) {

                if (SubDomain_.lattice_[i][j][k] != nullptr) {

                    for (int dir = 0; dir < _nLatNodes; dir++) {

                        if (SubDomain_.lattice_[i][j][k]->m_neighbours[dir] != nullptr) {
                            SubDomain_.lattice_[i][j][k]->m_distributions[dir] -=
                                1 / m_relaxation * (SubDomain_.lattice_[i][j][k]->m_distributions[dir] -
                                    SubDomain_.lattice_[i][j][k]->Equilibrium(SubDomain_.GetVelSet(), dir));
                        }

                    }
                }
            }
        }
    }
}

void Simulation::Streaming()
{
    for (int k = 0; k < SubDomain_.my_Nz_; k++) {
        for (int j = 0; j < SubDomain_.my_Ny_; j++) {
            for (int i = 0; i < SubDomain_.my_Nx_; i++) {
                if (SubDomain_.lattice_[i][j][k] != nullptr) {

                    for (int dir = 0; dir < _nLatNodes; dir++) {
                        if (SubDomain_.lattice_[i][j][k]->m_neighbours[dir] != nullptr) {
                            SubDomain_.lattice_[i][j][k]->Stream(dir);
                        }
                    }
                }
            }
        }
    }
}

void Simulation::ApplyWallBc() {

    for (int lat = 0; lat < SubDomain_.WallBoundaryNode_.size(); lat++) {

        int x_pos = SubDomain_.WallBoundaryNode_[lat]->x_position;
        int y_pos = SubDomain_.WallBoundaryNode_[lat]->y_position;
        int z_pos = SubDomain_.WallBoundaryNode_[lat]->z_position;

        for (int dir = 1; dir < _nLatNodes; dir++) {

            if (SubDomain_.lattice_[x_pos][y_pos][z_pos]->m_neighbours[dir] == nullptr) {
                int opp_dir = SubDomain_.GetVelSet()->OppositeDirection(dir);
                  Boundary_.Apply_BounceBack(*(SubDomain_.lattice_[x_pos][y_pos][z_pos]), dir, opp_dir);
            }
        }
    }

}


void Boundary::Apply_ZouHe_Bc(std::shared_ptr<VelocitySet> velSet, Node& lat, std::vector<double>& vel)
{
    double rho = 1.0;

    double Ny = (1.0 / 2) * (lat.m_distributions[2] + lat.m_distributions[14] + lat.m_distributions[15]
        - lat.m_distributions[3] - lat.m_distributions[16] - lat.m_distributions[17]) - (1.0 / 3) * (rho * vel[1]);

    double Nz = (1.0 / 2) * (lat.m_distributions[4] + lat.m_distributions[10] + lat.m_distributions[14]
        - lat.m_distributions[5] - lat.m_distributions[15] - lat.m_distributions[17]) - (1.0 / 3) * (rho * vel[2]);


    lat.m_distributions[0] = lat.m_distributions[1]   + (1 / 3.0 * rho) * vel[0];

    lat.m_distributions[6] = lat.m_distributions[11]  + (1 / 6.0 * rho) * (vel[0] + vel[1]) - Ny;

    lat.m_distributions[7] = lat.m_distributions[10]  + (1 / 6.0 * rho) * (vel[0] - vel[1]) + Ny;

    lat.m_distributions[8] = lat.m_distributions[13]  + (1 / 6.0 * rho) * (vel[0] + vel[2]) - Nz;

    lat.m_distributions[9] = lat.m_distributions[12]  + (1 / 6.0 * rho) * (vel[0] - vel[2]) + Nz;

    rho = (1.0 / (vel[0] + 1.0)) * (lat.m_distributions[2] + lat.m_distributions[3] + lat.m_distributions[4]
        + lat.m_distributions[5] + lat.m_distributions[14] + lat.m_distributions[15] + lat.m_distributions[16]
        + lat.m_distributions[17] + lat.m_distributions[18] + 2 * (lat.m_distributions[1] + lat.m_distributions[10]
            + lat.m_distributions[11] + lat.m_distributions[12] + lat.m_distributions[13]));
}