Shan-Chen LBM两相流C++实现：从伪势力到VTK可视化

简介：本资源是一份面向计算流体力学初学者与C++编程学习者的两相流数值模拟实践代码，聚焦Lattice Boltzmann Method（LBM）与Shan-Chen多相模型的工程实现。它解决了二维两相流中界面演化、表面张力建模等关键问题，适用于高校课程设计、科研入门及CFD算法验证场景。压缩包为3KB的RAR文件，仅含1个核心源码文件shanchen.cpp，完整实现了D2Q9格子模型、势能计算、碰撞-传播迭代、密度场更新及基础边界处理，代码结构清晰、注释充分，便于理解LBM离散动力学框架与Shan-Chen相互作用力的嵌入逻辑。已有1846人学习下载，读者可直接编译运行，观察气液相分离、液滴合并等典型现象，掌握从理论公式到C++可执行代码的关键转化路径，并为拓展三维模拟或耦合传热模块提供坚实基础。

1. Shan-Chen 模型不是“加个力就能分相”的黑箱，它是用格子玻尔兹曼方法（LBM）在 C++ 中显式编码分子间作用的两相流模拟核心

很多刚接触计算流体力学的工程师看到“Shan-Chen 模型”第一反应是：不就是 LBM 里加个伪势函数吗？改几行 force 计算就完事？结果一跑代码，密度场震荡发散、液滴不聚并、气液界面模糊成一片灰——问题不在编译器，而在对模型物理本质的误读。Shan-Chen 模型的本质，是将连续介质中复杂的分子间作用（如范德华力）离散化为格点邻域内的密度加权相互作用，它不求解纳维-斯托克斯方程，而是通过分布函数演化+非局部力耦合，让宏观两相行为从微观碰撞规则中自然涌现。这套机制对 C++ 实现提出刚性要求：必须严格控制内存布局（避免 cache miss 拖慢每步碰撞）、精确管理浮点精度（密度比超 100:1 时单精度易溢出）、显式分离流场更新与力计算时序（否则出现非物理振荡）。本文面向已掌握基础 LBM 概念、正用 VSCode 或 Visual Studio 编写 C++ 数值模拟代码的从业者，不讲推导，只拆解一个可本地编译、带边界验证、能输出 VTK 可视化数据的真实项目骨架——从shanchen_两相流Shan-Chen模型_C++这个标题出发，把LBMShan-Chen_LBM_shanchen模型落到每一行#include和每个for (int i = 0; i < nx; ++i)里。

2. 用 C++ 构建 Shan-Chen LBM 的最小可运行骨架：从格点定义、分布函数初始化到平衡态计算

Shan-Chen 模型的 C++ 实现绝非在标准 D2Q9 LBM 上简单叠加 force term。它的结构刚性体现在三个不可简化的层级：格点拓扑定义 → 分布函数内存布局 → 平衡态与伪势力的耦合时序。跳过任一层，都会导致后续所有优化失效。

2.1 定义 D2Q9 格子拓扑与内存对齐的格点数组

Shan-Chen 模型必须使用固定速度集（如 D2Q9），其方向向量和权重是硬编码常量。C++ 中若用std::vector<std::array<double, 9>>存储分布函数，会因动态分配引入 cache 不友好访问。生产级实现应采用一维连续数组 + 手动索引映射：

// constants.h constexpr int Q = 9; constexpr int cx[Q] = {0, 1, 0, -1, 0, 1, -1, -1, 1}; // x-direction of velocity vectors constexpr int cy[Q] = {0, 0, 1, 0, -1, 1, 1, -1, -1}; // y-direction constexpr double w[Q] = {4.0/9.0, 1.0/9.0, 1.0/9.0, 1.0/9.0, 1.0/9.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0};

提示：cx,cy,w必须声明为constexpr，确保编译期常量折叠，避免运行时查表开销。Q=9是 D2Q9 的硬约束，不可改为#define——C++20 要求模板参数必须是字面类型。

格点密度与分布函数需严格分离存储，且密度数组必须支持快速邻域求和（伪势计算核心）：

// lattice.h class Lattice2D { public: const int nx, ny; std::vector<double> rho; // size = nx * ny, density at each node std::vector<double> feq; // size = nx * ny * Q, equilibrium f_i std::vector<double> f; // size = nx * ny * Q, current f_i std::vector<double> f_new; // size = nx * ny * Q, post-collision f_i Lattice2D(int _nx, int _ny) : nx(_nx), ny(_ny), rho(nx * ny, 1.0), // initial uniform density feq(nx * ny * Q, 0.0), f(nx * ny * Q, 0.0), f_new(nx * ny * Q, 0.0) {} // 一维索引转二维坐标，避免除法（性能关键） inline int idx(int i, int j) const { return j * nx + i; } inline int f_idx(int i, int j, int q) const { return (j * nx + i) * Q + q; } };

2.2 实现 Shan-Chen 平衡态与伪势力：密度加权与力项分离

Shan-Chen 的核心创新在于：平衡态feq不仅依赖局部密度rho[i][j]和速度u，还隐含了非局部力的影响；而力本身由邻域密度加权和生成。二者必须解耦计算，否则产生自引用循环。标准做法是：先用当前rho计算feq，再用feq推出宏观速度u，最后用rho邻域和计算力F，并将F注入碰撞项。

// shanchen_force.h #include "lattice.h" #include "constants.h" class ShanChenForce { private: const double G; // interaction strength, typically -1.0 to -5.0 const double psi0; // reference density for pseudo-potential, e.g., 0.5 // Pseudo-potential function: psi(rho) = rho0 * (1 - exp(-rho/rho0)) inline double psi(double rho) const { return psi0 * (1.0 - std::exp(-rho / psi0)); } public: ShanChenForce(double _G = -1.0, double _psi0 = 0.5) : G(_G), psi0(_psi0) {} // Compute force F_x, F_y at node (i,j) using 8-neighbour sum void computeForce(const Lattice2D& lat, std::vector<double>& Fx, std::vector<double>& Fy) { const int total_nodes = lat.nx * lat.ny; for (int j = 0; j < lat.ny; ++j) { for (int i = 0; i < lat.nx; ++i) { double fx = 0.0, fy = 0.0; const int center = lat.idx(i, j); const double psi_center = psi(lat.rho[center]); // Sum over 8 neighbours (exclude center itself) for (int dq = 1; dq < Q; ++dq) { // skip q=0 (rest particle) int ni = i + cx[dq]; int nj = j + cy[dq]; if (ni >= 0 && ni < lat.nx && nj >= 0 && nj < lat.ny) { const int nidx = lat.idx(ni, nj); const double psi_neigh = psi(lat.rho[nidx]); fx += cx[dq] * psi_center * psi_neigh; fy += cy[dq] * psi_center * psi_neigh; } } Fx[center] = G * fx; Fy[center] = G * fy; } } } };

注意：computeForce中psi函数必须用std::exp而非近似多项式——在低密度区（rho << psi0）多项式会严重失真，导致气相力计算错误。G为负值才产生吸引力，这是两相分离的物理前提；若设为正，系统将坍缩成单点。

2.3 碰撞与传播：将力项注入 BGK 碰撞算子

Shan-Chen 的力不直接修改分布函数，而是作为额外项加入 BGK 碰撞项。标准 BGK 为f_i^{new} = f_i - 1/tau * (f_i - f_i^{eq})，Shan-Chen 扩展为：

f_i^{new} = f_i - 1/tau * (f_i - f_i^{eq}) + (1 - 1/(2*tau)) * F_i

其中F_i = w_i * (c_i - u) · F / (c_s^2 * rho)是力在第i个速度方向的投影。c_s^2 = 1/3（D2Q9），w_i为权重：

// lbm_solver.h void collideAndStream(Lattice2D& lat, const std::vector<double>& Fx, const std::vector<double>& Fy, double tau) { const int total_nodes = lat.nx * lat.ny; const double cs2 = 1.0 / 3.0; // Step 1: Compute macroscopic velocity u_x, u_y at each node std::vector<double> ux(total_nodes, 0.0), uy(total_nodes, 0.0); for (int j = 0; j < lat.ny; ++j) { for (int i = 0; i < lat.nx; ++i) { const int idx2d = lat.idx(i, j); double sum_f_cx = 0.0, sum_f_cy = 0.0; for (int q = 0; q < Q; ++q) { const int fidx = lat.f_idx(i, j, q); sum_f_cx += lat.f[fidx] * cx[q]; sum_f_cy += lat.f[fidx] * cy[q]; } ux[idx2d] = sum_f_cx / lat.rho[idx2d]; uy[idx2d] = sum_f_cy / lat.rho[idx2d]; } } // Step 2: Compute equilibrium f_eq and force term F_i for (int j = 0; j < lat.ny; ++j) { for (int i = 0; i < lat.nx; ++i) { const int idx2d = lat.idx(i, j); const double rho_local = lat.rho[idx2d]; const double ux_local = ux[idx2d]; const double uy_local = uy[idx2d]; // Equilibrium: f_i^eq = w_i * rho * [1 + 3*(c_i·u)/c_s^2 + 4.5*(c_i·u)^2/c_s^4 - 1.5*u^2/c_s^2] for (int q = 0; q < Q; ++q) { const double cu = cx[q] * ux_local + cy[q] * uy_local; const double u2 = ux_local * ux_local + uy_local * uy_local; const double feq_val = w[q] * rho_local * (1.0 + 3.0 * cu / cs2 + 4.5 * cu * cu / (cs2 * cs2) - 1.5 * u2 / cs2); lat.feq[lat.f_idx(i, j, q)] = feq_val; } // Force term: F_i = w_i * (c_i - u) · F / (c_s^2 * rho) const double fx_local = Fx[idx2d]; const double fy_local = Fy[idx2d]; for (int q = 0; q < Q; ++q) { const double cx_q = static_cast<double>(cx[q]); const double cy_q = static_cast<double>(cy[q]); const double c_minus_u_x = cx_q - ux_local; const double c_minus_u_y = cy_q - uy_local; const double force_proj = w[q] * (c_minus_u_x * fx_local + c_minus_u_y * fy_local) / (cs2 * rho_local); // Apply BGK with force: f_new = f - (f - f_eq)/tau + (1 - 0.5/tau) * force_proj const int fidx = lat.f_idx(i, j, q); lat.f_new[fidx] = lat.f[fidx] - (lat.f[fidx] - lat.feq[fidx]) / tau + (1.0 - 0.5 / tau) * force_proj; } } } // Step 3: Streaming (bounce-back boundaries handled separately) for (int j = 0; j < lat.ny; ++j) { for (int i = 0; i < lat.nx; ++i) { for (int q = 0; q < Q; ++q) { int ni = i - cx[q]; // reverse direction for streaming int nj = j - cy[q]; if (ni >= 0 && ni < lat.nx && nj >= 0 && nj < lat.ny) { lat.f[lat.f_idx(ni, nj, q)] = lat.f_new[lat.f_idx(i, j, q)]; } // Bounce-back for boundaries: set f_i(new) = f_{opp(i)}(old) at wall else { const int opp_q = getOpposite(q); // defined as [0,3,4,1,2,7,8,5,6] for D2Q9 lat.f[lat.f_idx(i, j, opp_q)] = lat.f_new[lat.f_idx(i, j, q)]; } } } } }

提示：getOpposite(q)是 D2Q9 的固定映射：q=0→0（静止），q=1→3（东↔西），q=2→4（北↔南），q=5→7（东北↔西南），q=6→8（东南↔西北）。必须硬编码，不可用公式推导——避免分支预测失败。

3. 在 VSCode 中配置 C++ 编译环境并验证 Shan-Chen 模型的两相分离行为

VSCode 本身不编译代码，它依赖外部构建系统。对数值模拟项目，必须放弃tasks.json的简单命令拼接，改用 CMakeLists.txt 驱动 Ninja 构建——因为 LBM 涉及 OpenMP 并行、VTK 输出、高精度数学库链接，手动写g++命令极易遗漏-march=native -O3 -ffast-math等关键 flag。

3.1 编写 CMakeLists.txt：启用 OpenMP 与 VTK 支持

# CMakeLists.txt cmake_minimum_required(VERSION 3.10) project(ShanChenLBM CXX) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) # Find required packages find_package(OpenMP REQUIRED) find_package(VTK REQUIRED COMPONENTS vtkIOXML vtkCommonCore) # Compiler flags for HPC if(MSVC) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /O2 /arch:AVX2 /fp:fast") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /D_VCRT_SECURE_NO_WARNINGS") else() set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -march=native -ffast-math -funroll-loops") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp") endif() # Add executable add_executable(shanchen_main main.cpp lattice.h constants.h shanchen_force.h lbm_solver.h) # Link libraries target_link_libraries(shanchen_main ${OpenMP_CXX_LIBRARIES} ${VTK_LIBRARIES}) target_include_directories(shanchen_main PRIVATE ${VTK_INCLUDE_DIRS}) # For Windows: ensure VTK DLLs are in output dir if(WIN32) add_custom_command(TARGET shanchen_main POST_BUILD COMMAND ${CMAKE_COMMAND} -E copy_if_different $<TARGET_FILE_DIR:${VTK_LIBRARIES}>/vtkCommonCore-9.2.dll $<TARGET_FILE_DIR:shanchen_main>/vtkCommonCore-9.2.dll) endif()

3.2 主程序main.cpp：初始化、迭代、输出 VTK 文件

一个可验证的最小主程序必须包含：初始密度扰动（如中心高密度斑块）、足够迭代步数（>10000）、VTK 格式输出（供 Paraview 查看相分离）。不能只打印“simulation done”。

// main.cpp #include <iostream> #include <vector> #include <cmath> #include <fstream> #include "lattice.h" #include "shanchen_force.h" #include "lbm_solver.h" // Write VTK image data file for Paraview void writeVTK(const Lattice2D& lat, int step) { std::string filename = "shanchen_" + std::to_string(step) + ".vti"; std::ofstream file(filename); file << "<?xml version=\"1.0\"?>\n"; file << "<VTKFile type=\"ImageData\" version=\"0.1\" byte_order=\"LittleEndian\">\n"; file << " <ImageData WholeExtent=\"0 " << (lat.nx-1) << " 0 " << (lat.ny-1) << " 0 0\" " << "Origin=\"0 0 0\" Spacing=\"1 1 1\">\n"; file << " <Piece Extent=\"0 " << (lat.nx-1) << " 0 " << (lat.ny-1) << " 0 0\">\n"; file << " <PointData Scalars=\"density\">\n"; file << " <DataArray type=\"Float64\" Name=\"density\" Format=\"ascii\">\n"; for (int j = 0; j < lat.ny; ++j) { for (int i = 0; i < lat.nx; ++i) { file << lat.rho[lat.idx(i, j)] << " "; } file << "\n"; } file << " </DataArray>\n"; file << " </PointData>\n"; file << " </Piece>\n"; file << " </ImageData>\n"; file << "</VTKFile>\n"; file.close(); } int main() { const int nx = 128, ny = 128; const int max_iter = 20000; const double tau = 0.6; // relaxation time, must be > 0.5 for stability Lattice2D lat(nx, ny); ShanChenForce force(-2.0, 0.5); // G=-2.0, psi0=0.5 // Initialize: central high-density region (liquid droplet) const int r0 = 15; for (int j = ny/2 - r0; j <= ny/2 + r0; ++j) { for (int i = nx/2 - r0; i <= nx/2 + r0; ++i) { if ((i - nx/2)*(i - nx/2) + (j - ny/2)*(j - ny/2) <= r0*r0) { lat.rho[lat.idx(i, j)] = 2.0; // liquid phase } else { lat.rho[lat.idx(i, j)] = 0.1; // vapor phase } } } // Pre-allocate force arrays std::vector<double> Fx(lat.nx * lat.ny, 0.0), Fy(lat.nx * lat.ny, 0.0); std::cout << "Starting Shan-Chen LBM simulation...\n"; for (int iter = 0; iter < max_iter; ++iter) { // Update density from distribution function for (int j = 0; j < lat.ny; ++j) { for (int i = 0; i < lat.nx; ++i) { double sum_f = 0.0; for (int q = 0; q < Q; ++q) { sum_f += lat.f[lat.f_idx(i, j, q)]; } lat.rho[lat.idx(i, j)] = sum_f; } } // Compute force force.computeForce(lat, Fx, Fy); // Collide & stream collideAndStream(lat, Fx, Fy, tau); // Output every 1000 steps if (iter % 1000 == 0) { std::cout << "Step " << iter << ", max rho = " << *std::max_element(lat.rho.begin(), lat.rho.end()) << "\n"; writeVTK(lat, iter); } } std::cout << "Simulation finished.\n"; return 0; }

3.3 VSCode 配置：c_cpp_properties.json与构建流程

在 VSCode 中按Ctrl+Shift+P→ “C/C++: Edit Configurations (UI)”，设置以下关键项：

构建流程：

1. 在终端执行mkdir build && cd build
2. cmake -G "Ninja" -DCMAKE_BUILD_TYPE=Release ..
3. ninja（生成shanchen_main.exe）
4. 运行./shanchen_main，生成shanchen_0000.vti等文件
5. 用 Paraview 打开.vti文件，添加Warp By Scalar滤镜观察液滴形变

注意：若遇到undefined reference to 'vtkCommonCore::Initialize()'，说明 VTK 链接顺序错误。在CMakeLists.txt中将target_link_libraries改为target_link_libraries(shanchen_main ${VTK_LIBRARIES} ${OpenMP_CXX_LIBRARIES})，VTK 库必须在前。

4. 调优 Shan-Chen 模型的 3 个必调参数：G、psi0、tau 与两相密度比的定量关系

Shan-Chen 模型的物理真实性完全由三个参数控制：相互作用强度G、伪势参考密度psi0、松弛时间tau。它们不独立，而是共同决定两相共存密度比rho_liq / rho_vap。盲目调参只会得到无物理意义的“漂亮图”，而非可复现的相图。

4.1 G 与 psi0 共同决定两相密度比，tau 控制动力学粘度

理论表明，在 D2Q9 Shan-Chen 模型中，两相共存密度满足隐式方程：

rho_vap = psi0 * W0( exp(-1) * exp( -G * psi0 * (1 - rho_vap/psi0) ) )

rho_liq = psi0 * W_{-1}( exp(-1) * exp( -G * psi0 * (1 - rho_liq/psi0) ) )

其中W0,W_{-1}是朗伯 W 函数的两个实数分支。实际应用中，我们不求解该方程，而是通过预计算表格建立G-psi0到rho_ratio的映射：

提示：rho_ratio > 30时，单精度float会因rho_liq - rho_vap有效位不足导致界面弥散。必须用double，且tau需同步增大以维持稳定性。

4.2 tau 的物理意义与稳定边界：从粘度公式反推安全范围

tau直接控制流体运动粘度nu = cs^2 * (tau - 0.5)。但更重要的是，tau决定了数值稳定性上限。当G很大（强相互作用）时，力项放大误差，要求tau > 0.5 + |G| * psi0 * 0.1。经验公式：

tau_min ≈ 0.5 + 0.05 * |G| * psi0

例如G = -3.0, psi0 = 0.5→tau_min ≈ 0.575。若仍用tau = 0.6，则迭代 5000 步后密度场开始高频震荡；此时应设为tau = 0.7，虽降低nu，但保证收敛。

4.3 验证两相分离的 3 个量化指标：必须写进日志

每次运行后，不能只看图片，必须输出以下三项到log.txt：

1. 界面厚度（Interface Width）：沿液滴直径取线，计算rho从0.1*rho_liq到0.9*rho_liq的像素数。理想值为 4~6 格点（D2Q9 精度极限）。
2. 液滴圆度（Circularity）：4π * Area / Perimeter²，>0.95 为良好数值各向同性。
3. 质量守恒误差（Mass Error）：|sum(rho_final) - sum(rho_initial)| / sum(rho_initial)，应 < 1e-10（双精度下）。

// 在 main.cpp 结尾添加 double total_mass_init = 0.0, total_mass_final = 0.0; for (double r : lat.rho) total_mass_init += r; // ... after simulation ... for (double r : lat.rho) total_mass_final += r; std::cout << "Mass error: " << std::abs(total_mass_final - total_mass_init) / total_mass_init << "\n"; // Interface width estimation (simplified) int width_count = 0; for (int i = 0; i < lat.nx; ++i) { double r = lat.rho[lat.idx(i, lat.ny/2)]; if (r > 0.1 * 1.6 && r < 0.9 * 1.6) width_count++; } std::cout << "Interface width (px): " << width_count << "\n";

5. 解决 Shan-Chen 模型在 C++ 实现中最常见的 3 类崩溃与发散：内存越界、力项符号错误、tau 设置失当

生产环境中，90% 的 Shan-Chen 模拟失败并非算法错误，而是 C++ 层面的工程疏漏。以下三类问题在 VSCode 调试器中表现为Segmentation fault或NaN密度，必须逐条排查。

5.1 内存越界：f_idx计算错误导致f数组写爆

最隐蔽的 bug 是f_idx(i, j, q)公式错误。常见错误写法：

// 错误！会导致 j*nx+i 超出 nx*ny*Q 范围 int f_idx(int i, int j, int q) { return j * ny + i + q; } // 混淆 nx/ny // 更危险的错误：未检查 i,j 边界就计算 int f_idx(int i, int j, int q) { return (j * nx + i) * Q + q; } // i,j 越界时仍计算

正确做法是：在debug模式下启用断言，并在每次f_idx调用前校验：

inline int f_idx(int i, int j, int q) const { assert(i >= 0 && i < nx && j >= 0 && j < ny && q >= 0 && q < Q); return (j * nx + i) * Q + q; }

编译时加-DDEBUG -g，运行时报错位置直指越界坐标。

5.2 力项符号错误：G 为正或psi函数返回负值

G必须为负，否则力为排斥力，液滴炸裂。但更隐蔽的是psi(rho)在rho=0时返回0，导致psi_center * psi_neigh = 0，力全为零。psi0=0.5时，rho=0.01的psi=0.01，但rho=0的psi=0，造成气相无力。解决方案：给rho设下限：

// 在 computeForce 中 const double rho_safe = std::max(lat.rho[center], 1e-10); const double psi_center = psi0 * (1.0 - std::exp(-rho_safe / psi0));

5.3 tau 设置失当引发的 NaN 瀑布

当tau过小（如tau=0.51）且G很大时，1/tau项主导，f_new迅速溢出为inf，下一步inf/inf得NaN，随后全数组污染。检测方法：在每次collideAndStream后插入：

// 在 collideAndStream 结尾添加 for (double v : lat.f_new) { if (std::isnan(v) || std::isinf(v)) { std::cerr << "NaN/Inf detected in f_new at iter " << iter << "\n"; exit(1); } }

修复策略：若报错，立即增大tau至0.5 + 0.1*|G|*psi0，并检查G是否误设为正。

提示：Windows 下若遇0xC0000094 Integer division by zero，不是除零，而是rho为0导致ux = sum_f / 0。必须在速度计算前加if (rho_local < 1e-10) { ux_local = 0; uy_local = 0; continue; }。

本文还有配套的精品资源，点击获取