如何在流体模拟MIC(0)预处理中整合Neumann边界条件防压力爆炸

问题背景 #

基于Robert Bridson的《Fluid Simulation Notes》实现欧拉流体模拟器，替换原Gauss-Seidel压力求解器为MICCG0（Modified Incomplete Cholesky Conjugate Gradient, Level Zero）后，液体接触固体边界时求解器崩溃：压力值变为-NaN(ind)，投影后速度也出现NaN，导致半拉格朗日平流步骤x_prev = x - v * dt无法执行。核心需求是将Neumann边界条件正确整合到Poisson方程的MIC(0)预处理流程中，避免压力值异常激增。

核心问题根源 #

流体模拟中，固体壁面的Neumann边界条件对应法向速度为0，等价于压力的法向导数为0。若预处理或CG迭代阶段未正确处理该条件，会导致矩阵奇异性、数值溢出，最终出现NaN。

具体实现方案 #

1. MIC(0)预处理阶段的边界处理 #

构建下三角矩阵L时，需跳过固体边界的邻居方向系数，并调整对角元确保非零：

• 跳过固体网格点，不对其进行预处理计算
• 对于流体网格点，若邻居为固体，跳过该方向的系数（符合Neumann条件的梯度要求）
• 对角元添加极小epsilon（如1e-6），避免因边界条件导致对角元为0，引发除零错误

2. CG迭代全流程的边界约束 #

在残差计算、搜索方向更新、预处理的前后向替换步骤中，均需跳过固体网格点，或固定其压力值以符合边界条件。

核心代码片段 #

// 构建MIC(0)下三角预处理矩阵L
void buildMIC0(const Grid& grid, const std::vector<int>& solid_mask, std::vector<double>& L) {
    const int N = grid.totalCells();
    for (int i = 0; i < N; ++i) {
        if (solid_mask[i]) continue; // 跳过固体网格
        
        double diagSum = 0.0;
        const auto& neighbors = grid.getFluidNeighbors(i);
        L[i * N + i] = 1.0; // 初始化对角元
        
        for (int j : neighbors) {
            if (solid_mask[j]) continue; // 跳过固体邻居，应用Neumann条件
            
            const double coeff = grid.getPoissonCoeff(i, j);
            diagSum += std::abs(coeff);
            if (j < i) { // 仅处理下三角元素
                L[i * N + j] = coeff / L[j * N + j];
            }
        }
        
        // 确保对角元非零，添加epsilon防止除零
        L[i * N + i] = std::sqrt(diagSum + 1e-6);
    }
}

// MIC(0)前向替换（求解Ly = r）
void forwardSubMIC0(const std::vector<double>& L, const std::vector<int>& solid_mask, 
                    std::vector<double>& y, const std::vector<double>& r, const Grid& grid) {
    const int N = y.size();
    for (int i = 0; i < N; ++i) {
        if (solid_mask[i]) {
            y[i] = 0.0; // 固体边界压力固定为0（符合Neumann条件）
            continue;
        }
        
        double sum = r[i];
        const auto& neighbors = grid.getFluidNeighbors(i);
        for (int j : neighbors) {
            if (j >= i || solid_mask[j]) continue;
            sum -= L[i * N + j] * y[j];
        }
        y[i] = sum / L[i * N + i];
    }
}

// MIC(0)后向替换（求解L^T z = y）
void backwardSubMIC0(const std::vector<double>& L, const std::vector<int>& solid_mask, 
                     std::vector<double>& z, const std::vector<double>& y, const Grid& grid) {
    const int N = z.size();
    for (int i = N - 1; i >= 0; --i) {
        if (solid_mask[i]) {
            z[i] = 0.0;
            continue;
        }
        
        double sum = y[i];
        const auto& neighbors = grid.getFluidNeighbors(i);
        for (int j : neighbors) {
            if (j <= i || solid_mask[j]) continue;
            sum -= L[j * N + i] * z[j]; // 转置矩阵取L[j][i]
        }
        z[i] = sum / L[i * N + i];
    }
}

// 计算Poisson方程残差r = b - A*p
void computeResidual(const Grid& grid, const std::vector<int>& solid_mask, 
                     std::vector<double>& r, const std::vector<double>& p, const std::vector<double>& b) {
    const int N = r.size();
    for (int i = 0; i < N; ++i) {
        if (solid_mask[i]) {
            r[i] = 0.0;
            continue;
        }
        
        double Ap = 0.0;
        const auto& neighbors = grid.getFluidNeighbors(i);
        for (int j : neighbors) {
            if (solid_mask[j]) continue;
            Ap += grid.getPoissonCoeff(i, j) * p[j];
        }
        r[i] = b[i] - Ap;
    }
}

关键注意事项 #

• 所有矩阵操作必须严格跳过固体网格点，避免无效计算引入数值错误
• 预处理矩阵的对角元必须添加epsilon，防止边界条件下对角元为0引发除零
• Poisson方程右端项（速度散度）的计算需排除固体边界方向的速度分量，确保与Neumann条件一致

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