企业官网建设流程全解析

简介：这是围绕QSGS方法重构二维多孔介质并开展LBM流动模拟的C++源代码，面向从事孔隙尺度数值计算、多孔介质微观重构及流体仿真研究的工程师、科研人员或相关专业研究生。压缩包仅含1个cpp文件，大小约2KB，代码紧凑，省略了工程框架，便于聚焦核心算法。目前已有451人学习下载。程序以Lattice Boltzmann Method为核心，结合Quasi-Spectral Ghost-Fluid Scheme处理复杂边界，演示了从多孔结构生成、边界条件设置到流场求解的完整流程。通过研读代码，可以快速掌握LBM离散模型的搭建细节、QSGS边界处理的实现逻辑，以及孔隙率、渗透率等结构参数如何影响流动模拟结果。对于希望验证算法、搭建二维多孔介质模拟模型或进一步扩展三维与多物理场应用的读者，这是一份简洁实用的入门参考模板，能够缩短前期编码调试时间，帮助将精力集中在物理机制分析与结果解释上。

1. 二维多孔介质重构为什么值得用C++重新写一遍

打开一个名为“二维多孔介质.cpp.zip”的压缩包，里面大概率是一个读入参数、生成网格、输出掩膜文件的控制台工程。这类代码在学术圈里流传很广，但真正要拿来跑自己的模拟时，往往卡在第一步：孔隙率对不上、连通性太差、生长方向不自然。QSGS（四参数随机生长法）是目前重构多孔介质最常用的方法之一，它用“随机撒核 + 定向生长”两个步骤，从统计上逼近真实岩石或纤维材料的几何结构。用 C++ 实现 QSGS，不是因为别的语言不行，而是当你需要在 1000×1000 网格上做几百次蒙特卡洛重复时，C++ 的循环和内存局部性优势就会直接变成时间成本。本文不做理论综述，只讲我用 C++ 落地 QSGS 时的数据表示、参数标定、连通性校验和工程化做法，让重构结果不仅“看着像”，还能算得准。

2. QSGS 重构算法的成核、生长与边界处理方式

2.1 成核阶段：核心概率与实际核心数目的换算

QSGS 的第一步是在网格上随机选取“核心相”单元，这些单元是后续生长的种子。假设网格总单元数为 (N)，用户设置核心概率 (p_c)，实际撒下的核心数期望值是 (N \times p_c)。在 C++ 里实现非常简单，但有一个细节容易引发偏差：随机数引擎的分布和网格遍历顺序。

std::mt19937_64 rng(seed); std::uniform_real_distribution<double> dist(0.0, 1.0); std::vector<unsigned char> grid(N, 0); // 0=孔隙, 1=固相 for (size_t i = 0; i < N; ++i) { if (dist(rng) < pc) { grid[i] = 1; } }

这段代码用均匀分布判断每个单元是否成为核心。注意这里的核心概率是“单元级概率”，不是“网格中必须有多少核心”。由于随机性，实际核心数会上下波动，对于 100×100 网格误差可能到 ±5%，在小网格上做对照实验时建议记录实际核心数，而不是直接用期望值算。另一个常见错误是随机数引擎的种子固定后，每次运行结果完全相同，若需要批量统计，应把种子作为命令行参数传入。

2.2 生长阶段：邻域遍历、定向概率与迭代轮数

生长阶段逐轮扫描网格，当前单元若为孔隙，且其邻域存在固相核心或已生长单元，则按生长概率 (p_g) 转为固相。邻域一般选四邻域或八邻域，前者生成的骨架更接近细长通道，后者更容易产生大面积团簇。QSGS 的核心参数是方向生长概率，常见的做法是给六个方向分别设置概率 (p_{g,1}\sim p_{g,6})，但二维下通常只区分四个方向或两个主轴。

// 方向向量: 上、下、左、右 const int dx[4] = {0, 0, -1, 1}; const int dy[4] = {-1, 1, 0, 0}; for (int step = 0; step < max_steps; ++step) { for (int idx = 0; idx < N; ++idx) { if (grid[idx]) continue; int x = idx % cols, y = idx / cols; for (int d = 0; d < 4; ++d) { int nx = x + dx[d], ny = y + dy[d]; if (nx < 0 || nx >= cols || ny < 0 || ny >= rows) continue; if (grid[ny * cols + nx] && dist(rng) < pg[d]) { grid[idx] = 1; break; } } } }

这段代码的问题在于扫描顺序固定，导致左上角先扫描的孔隙更容易被生长，生成的固相分布带有方向偏差。我常用的做法是生成一个索引数组，每轮开始前std::shuffle打乱索引，再按随机顺序生长。另一个需要注意的问题是迭代轮数不能设成“直到填满”，固相分数在生长后期会快速逼近 1，孔隙结构会被完全吞没，正确做法是每轮统计当前固相分数，达到目标值后立即停止。

2.3 周期性边界：避免重构样本边缘效应

多孔介质重构通常用于后续的渗透率、扩散率模拟，边界条件如果处理不好，统计结果会失真。QSGS 的邻居查找可以支持周期性边界，也就是右侧邻居回绕到左侧。实现时不需要真的把网格复制三份，只需要在坐标越界时取模：

int nx = (x + dx[0] + cols) % cols; int ny = (y + dy[0] + rows) % rows;

采用周期性边界后，生成的结构在几何上等价于无限大介质的单周期重复，做两相相对渗透率模拟时不会因为边界上的死孔隙造成误差。我建议在生成阶段就同步使用周期性邻域，否则生长出的边界区域通常更致密，肉眼不易察觉，但会在傅里叶变换统计中形成一个方向性伪影。

3. 用 C++ 实现 QSGS 二维重构的最小可运行工程

3.1 数据布局与随机数引擎选型

二维网格在 C++ 里最常见的存储方式有四种：二维vector<vector<char>>、一维vector<char>、std::array固定大小数组、以及裸指针new char[]。对于重构这种密集循环任务，二维vector的内存连续性不如一维数组，每次grid[x][y]都隐含两次指针跳转，在 500×500 网格上迭代 50 轮后差距明显。我一般用一维vector<unsigned char>，索引计算为y * cols + x，这样不仅能减少 cache miss，后续写文件时也能直接memcpy整块输出。

随机数生成建议使用 C++11 提供的<random>库，而不是旧的rand()配合srand。mt19937_64生成速度快，周期长，配合正态分布和均匀分布都能稳定复现。在 Windows 上用 MinGW 编译时，random_device可能每次返回相同序列，最好允许用户通过--seed指定随机种子。

3.2 初始化核心相与首轮生长

完整生成逻辑可以封装成一个类，下面是一个最小可运行版本的核心循环：

class QSGSGenerator { public: QSGSGenerator(int rows, int cols, double pc, double pg, int max_steps, double target_porosity) : rows_(rows), cols_(cols), pc_(pc), pg_(pg), max_steps_(max_steps), target_(target_porosity) { N_ = rows * cols; grid_.assign(N_, 0); } void Run(uint64_t seed) { std::mt19937_64 rng(seed); std::uniform_real_distribution<double> dist(0.0, 1.0); // 成核 for (auto& cell : grid_) { if (dist(rng) < pc_) cell = 1; } // 生长 std::vector<size_t> index(N_); for (size_t i = 0; i < N_; ++i) index[i] = i; for (int step = 0; step < max_steps_; ++step) { std::shuffle(index.begin(), index.end(), rng); bool changed = false; for (size_t pos : index) { if (grid_[pos]) continue; int x = static_cast<int>(pos % cols_); int y = static_cast<int>(pos / cols_); for (int d = 0; d < 4; ++d) { int nx = x + dx_[d], ny = y + dy_[d]; if (nx < 0 || nx >= cols_) nx = (nx + cols_) % cols_; if (ny < 0 || ny >= rows_) ny = (ny + rows_) % rows_; if (grid_[ny * cols_ + nx]) { if (dist(rng) < pg_[d]) { grid_[pos] = 1; changed = true; } break; } } } double solid_frac = 1.0 - GetPorosity(); if (solid_frac >= 1.0 - target_) break; if (!changed) break; } } double GetPorosity() const { size_t sum = 0; for (unsigned char v : grid_) sum += (v == 0); return static_cast<double>(sum) / N_; } private: int rows_, cols_, N_; double pc_, target_; int max_steps_; double pg_[4] = {0.3, 0.3, 0.3, 0.3}; std::vector<unsigned char> grid_; const int dx_[4] = {0, 0, -1, 1}; const int dy_[4] = {-1, 1, 0, 0}; };

这段代码有两个关键点。第一，每轮生长前用std::shuffle打乱访问顺序，避免网格扫描方向引入偏差。第二，提前判断target_，固相分数到达目标即跳出，不会把孔隙全部吞掉。参数数组pg_放在类里便于按方向分别调试，比如想构造横向裂缝，可以把左右方向概率调成 0.6，上下调成 0.1。

3.3 输出 PGM 掩膜文件与 CSV 统计

重构结果需要让下游程序使用，常见输出格式是 PGM 或 CSV。PGM 是二进制 P5 格式，一个 8 位灰度值对应一个网格单元，既能直接用图片查看器打开，又方便 ImageJ 做孔隙分析。

void WritePGM(const std::string& path) { FILE* fp = fopen(path.c_str(), "wb"); fprintf(fp, "P5\n%d %d\n255\n", cols_, rows_); // grid_ 中 0=孔隙, 1=固相; 转成 0/255 方便显示 std::vector<unsigned char> out(N_); for (size_t i = 0; i < N_; ++i) out[i] = grid_[i] ? 255 : 0; fwrite(out.data(), 1, N_, fp); fclose(fp); }

PGM 写入时要注意fprintf与fwrite混用：先写文本头，再用二进制写像素数据，中间不需要换行符，如果加多余空格，很多查看器会报错。CSV 输出则不推荐，因为 500×500 网格转成逗号分隔文本会膨胀到十几 MB，除非下游是 Python 脚本，否则直接用二进制文件更划算。

4. 参数标定与结果验证：孔隙率、连通度与界面密度

4.1 核心概率、生长概率和迭代轮数的关系

QSGS 参数不是互相独立的。核心概率决定固相核心的密度，核心少而生长概率高时，生成的是类似雪花扩散的粗大簇；核心多而生长概率低时，结构接近孤立沙堆。为了让结果能用，通常需要做一组小规模网格的预实验。

以下是我在 200×200 网格上得到的典型参数对照关系，固相分数由代码直接统计：

不需要把上表当成精确标准，不同随机种子和边界条件会带来 5% 左右的浮动，但它说明一个规律：核心概率对最终结构的决定性远大于生长概率。当你发现生成的孔隙被孤立成一个个小岛时，通常先提高pc，而不是盲调pg。还有一种常见做法是固定目标孔隙率，用二分法搜索合适的生长概率，每轮只跑 30 步以内，速度很快。

4.2 连通性检测：用并查集区分有效孔隙和死孔隙

多孔介质的“重构质量”不能只看孔隙率。两个样本的孔隙率都是 0.3，一个孔隙完全连通，流体可以穿越；另一个孔隙彼此孤立，渗透率为零。QSGS 生成的结构能否用于渗流模拟，取决于孔隙相的连通度。最常见的判断工具是并查集（Union-Find），对孔隙相做连通分量标号。

class UnionFind { public: UnionFind(int n) : parent(n) { iota(parent.begin(), parent.end(), 0); } int Find(int x) { while (parent[x] != x) { parent[x] = parent[parent[x]]; x = parent[x]; } return x; } void Union(int a, int b) { int ra = Find(a), rb = Find(b); if (ra != rb) parent[rb] = ra; } private: vector<int> parent; }; double ComputePoreConnectivity(const vector<unsigned char>& grid, int rows, int cols) { int n = rows * cols; UnionFind uf(n); auto idx = [&](int x, int y) { return y * cols + x; }; for (int y = 0; y < rows; ++y) { for (int x = 0; x < cols; ++x) { if (grid[idx(x, y)] == 1) continue; if (x + 1 < cols && grid[idx(x + 1, y)] == 0) uf.Union(idx(x, y), idx(x + 1, y)); if (y + 1 < rows && grid[idx(x, y + 1)] == 0) uf.Union(idx(x, y), idx(x, y + 1)); } } // 统计最大孔隙团簇占比 vector<int> region(n, 0); for (int i = 0; i < n; ++i) if (grid[i] == 0) region[uf.Find(i)]++; int max_cluster = *max_element(region.begin(), region.end()); int sum_pores = count(grid.begin(), grid.end(), (unsigned char)0); return static_cast<double>(max_cluster) / sum_pores; }

连通度 (C_{pore}) 定义为最大孔隙团簇单元数除以总孔隙单元数。(C_{pore}) 越接近 1，说明孔隙越能形成穿越网络。如果 (C_{pore}) 小于 0.5，即使孔隙率达到 0.35，这个结构在实际模拟中也没有传导能力。注意并查集路径压缩写法必须配合按秩合并，否则在大网格上递归过深会导致栈溢出，上面代码用迭代式while (parent[x] != x)避免递归。

4.3 两点相关函数：用统计函数验证重构合理性

孔隙率只描述浓度，连通度描述拓扑，若要与其他重构算法比较，需要用两点相关函数 (S_2(r)) 判断几何分布是否合理。定义是随机取两个距离为 (r) 的点，两点都在同一相的概率。实际计算时可随机投放大量点对，或用傅里叶变换快速求自相关。C++ 里的朴素实现是双重循环，复杂度 (O(N^2))，但 200×200 网格只需采样几万点对就能得到平滑曲线。

一个快速验证技巧是计算沿 (x) 轴和沿 (y) 轴的两点相关函数，如果两条曲线差异超过 10%，说明重构结构存在明显各向异性。QSGS 的方向生长概率在默认全等时，理论上应生成各向同性结构；若你的代码逐轮固定扫描顺序，就会在两点相关函数上暴露出一条人工纹理。我在调参时经常同时打印孔隙率、连通度、沿 x/y 方向的相关长度，这三项比肉眼看图更能判断是否需要增大核心概率或降低迭代轮数。

5. 从 zip 工程到可复现实验：文件组织、编译与批量参数扫描

5.1 解压后先确认工程形态

拿到一个 .cpp.zip，我一般不会直接全选编译，而是先看压缩包根目录有没有 CMakeLists.txt 或 Makefile。只有单个 .cpp 文件的工程最容易处理，但通常需要手动补上编译选项。推荐命令：

unzip er_wei_duo_kong_ji_zhi.cpp.zip -d qsgs_src cd qsgs_src g++ -O2 -std=c++17 -o qsgs main.cpp -lpthread ./qsgs --rows 200 --cols 200 --pc 0.03 --pg 0.3 --steps 20 --target 0.7 --seed 42 --out sample.pgm

-O2是必要的，QSGS 生长循环对 CPU 寄存器分配非常敏感，O0 和 O2 在 1000×1000 网格上速度差距可达 5 倍。-std=c++17是为了使用std::shuffle和std::filesystem，如果是老式代码，可能只支持 C++11，那就去掉 filesystem 相关部分。

5.2 命令行参数解析与批量运行

生产级程序不应该在 main 里写死参数。用getopt或简单的手写解析都可以，我喜欢用一个轻量级结构体集中接收参数，然后通过一个批量脚本循环遍历参数组合。避免多个结果文件相互覆盖，输出文件名必须带上种子值和目标孔隙率：

for pc in 0.02 0.04 0.06; do for seed in 1 2 3; do ./qsgs --rows 300 --cols 300 --pc $pc --pg 0.3 \ --steps 15 --target 0.7 --seed $seed \ --out "result_pc${pc}_seed${seed}.pgm" done done

批量执行后，同一个pc参数在不同种子下会得到不同孔隙率。把它们全部算平均，才是这个参数组合的统计预期。若只跑一个种子就下结论，常常会被随机波动误导。

5.3 一个快速判定重构质量的技巧：同时检查压缩包内的统计文件

我每次跑重构实验，都会让程序额外输出一个同名.txt统计文件，里面只放三列：孔隙率、最大孔隙团簇占比、界面密度。界面密度指的是固相与孔隙相邻的边数除以网格面积，它反映了比表面积的大小。

double ComputeSurfaceDensity() { int interface_edges = 0; for (int y = 0; y < rows_; ++y) { for (int x = 0; x < cols_; ++x) { if (grid_[y * cols_ + x] == 0) continue; if (x + 1 < cols_ && grid_[y * cols_ + x + 1] == 0) interface_edges++; if (y + 1 < rows_ && grid_[(y + 1) * cols_ + x] == 0) interface_edges++; } } return interface_edges / static_cast<double>(N_); }

一句话经验：孔隙率在 0.25 到 0.45 之间时，界面密度通常在 0.15 到 0.45 之间。如果你的结果孔隙率 0.4，界面密度却不到 0.05，说明生成的是几个巨大固相块，QSGS 的随机撒核可能没有生效，大概率是随机数分布范围写错，把所有单元都设成了固相。反过来，界面密度高于 0.6 时，结构多半是棋盘格状的噪声，这时降低生长概率，让固相更团簇一些。把这个统计值作为重构质量的“指纹”，比任何肉眼观察都可靠。

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