D3Q19模型LBM代码的ZouHe与反弹边界条件问题排查

我基于D3Q19模型开发了复杂几何流体动力学的Lattice Boltzmann Method（LBM）代码，入口边界采用Zou-He方法，壁面用完全反弹边界条件，出口维持恒定压力。目前遇到两个问题：

• Zou-He方法未按预期设置入口速度
• 壁面速度分量不为零，密度稳定在约1.03而非目标值1.0，推测是分布函数处理出错

核心代码片段 #

Collision函数 #

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));
                        }

                    }
                }
            }
        }
    }
}

Streaming函数 #

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);
            }
        }
    }

}

Zou-He入口边界函数 #

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]));
}

Stream与Equilibrium函数 #

void Node::Stream(int dir)
{
    m_neighbours[dir]->m_distributions[dir] = m_distributions[dir];
}

double Node::Equilibrium(std::shared_ptr<VelocitySet> velSet, int dir)
{
    double du = Velocity(velSet)[0] * velSet->GetDirection(dir)[0] 
        + Velocity(velSet)[1] * velSet->GetDirection(dir)[1]
        + Velocity(velSet)[2] * velSet->GetDirection(dir)[2];
    double u_sqr = Velocity(velSet)[0] * Velocity(velSet)[0]
        + Velocity(velSet)[1] * Velocity(velSet)[1] 
        + Velocity(velSet)[2] * Velocity(velSet)[2];

    return velSet->GetWeight(dir) * Density() * (1 + 3.0 * du + 9.0 / 2.0 * du * du - 3.0 / 2.0 * u_sqr);
}

问题排查与解决建议 #

1. Zou-He边界实现问题 #

• 密度计算顺序错误：当前先假设rho=1.0计算Ny、Nz，再重新计算rho，但新rho未回代到之前的计算步骤，导致边界用错误密度值。正确流程是先通过已知分布函数算rho，再用该rho计算Ny、Nz和未知分布函数。
• 方向索引不匹配：D3Q19的方向索引需严格对应（如dir=0为静止方向，dir=1_{6为轴向，dir=7}18为对角向），检查velSet的方向定义是否与Zou-He代码中使用的索引完全一致。
• 公式符号/系数错误：对照Zou-He原始3D入口边界公式，检查Ny、Nz的符号和分布函数更新公式，避免符号颠倒或系数错误。

2. 壁面反弹边界问题 #

• 反弹方向错误：确认OppositeDirection(dir)返回正确反向索引（如dir对应(vx,vy,vz)，反向应为(-vx,-vy,-vz)），索引映射错误会导致分布函数反弹异常。
• 边界节点遍历逻辑：检查WallBoundaryNode_是否包含所有壁面节点，且m_neighbours[dir] == nullptr的判断是否准确（即dir是指向域外的方向）。
• 边界时机错误：LBM中反弹边界需在Streaming之后、Collision之前应用，确认代码执行顺序是否正确（正确顺序：Streaming → 边界条件 → Collision）。

3. 核心LBM流程问题 #

• Streaming逻辑致命错误：当前Node::Stream直接将m_distributions[dir]赋值给邻居的m_distributions[dir]，正确逻辑应为：将当前节点的f_i传递给邻居节点的反向方向分布函数，即m_neighbours[dir]->m_distributions[opp_dir] = m_distributions[dir]（opp_dir为dir的反向），方向传递错误会直接导致整个流场计算失效。
• Collision函数跳过边界方向：Collision中跳过m_neighbours[dir] == nullptr的方向，需确认边界节点的碰撞处理逻辑是否正确，避免遗漏必要的碰撞步骤。
• 平衡态函数验证：检查Equilibrium的系数是否正确，同时确认velSet->GetWeight(dir)返回的权重符合D3Q19标准（静止方向w0=1/3，轴向w1_{w6=1/18，对角向w7}w18=1/36）。

4. 调试建议 #

• 输出入口边界节点的分布函数、rho、Ny、Nz值，对比理论值定位偏差步骤。
• 输出壁面节点的分布函数，检查反弹后是否满足f_i = f_{opp_i}（完全反弹条件）。
• 初始化所有分布函数为平衡态（rho=1.0，u=0），运行少量步数，观察密度和速度变化，确认是否为初始化或流程错误导致偏移。

内容的提问来源于stack exchange，提问作者Resa