企业官网建设流程全解析

简介：一套专用于三维自然对流模拟的C++程序源码，面向流体力学、热工与CFD方向的科研人员和高年级学生，适合研究瑞利数低于10E7的RB自然对流问题。包内共1个cpp源文件，压缩包仅2KB，代码精简，方便直接阅读、修改与二次开发。目前已有155人学习下载。程序以“leave7pj”和“strugglemnm”两个核心概念为主线，将流场的速度、压力、温度等物理量拆分为双分布函数，分别捕捉对流与扩散过程，并耦合连续性方程、动量方程、能量方程及状态方程来描述流体运动与热传递。针对自然对流中常见的非线性与湍流特性，程序采用有限体积或有限元类数值方法，并借助Gauss-Seidel或Jacobi等迭代格式逼近稳定解，帮助读者理解从层流到湍流过渡阶段的关键计算细节。研读后能够掌握双分布函数在低瑞利数自然对流模拟中的具体落地方式，包括方程离散、迭代收敛条件以及边界处理思路，对于开展三维自然对流独立编程、算法验证或教学演示都具有参考价值。

1. 三维自然对流为什么需要双分布函数

做三维自然对流模拟的人大概都有这种经历：直接解 Navier-Stokes 方程时压力-速度耦合很麻烦，加温度场后还要再解一个能量方程，网格稍密就迭代不收敛。改用三维双分布函数之后，速度场和温度场分别用两个分布函数演化，压力不再是需要反复修正的全局约束，浮力项也可以直接写进碰撞算子。这个方案本质上是格子玻尔兹曼方法在传热问题上的标准扩展，在处理方腔对流、电子散热和建筑通风这类边界规则的三维自然对流（natural convection）问题时非常顺手。下面我不打算复述某个现成项目，而是把一套最常用的 D3Q19/D3Q7 双分布函数求解器从控制方程到验证方法完整讲一遍，新手能照着写代码，熟手可以直接跳到第 4 章看稳定性边界和第 5 章的 Nusselt 数验证。

2. 双分布函数的控制方程与三维格子模型

2.1 从 BGK 方程到密度与温度两个分布函数

双分布函数并不是一个抽象口号，而是把速度场和温度场拆成两个独立演化过程的格式。流场用密度分布函数 f 描述，温度场用温度分布函数 g 描述，两者都满足带碰撞项的格子 BGK 方程。在三维坐标下，演化格式可以简洁地表示为：

# f 的演化：迁移 + 碰撞 # f_i(x + e_i*dt, t + dt) - f_i(x, t) = - (f_i - f_i_eq) / tau_f + F_i # g 的演化：同样结构但碰撞项更简单 # g_i(x + e_i*dt, t + dt) - g_i(x, t) = - (g_i - g_i_eq) / tau_g

这里 e_i 是第 i 个离散速度方向，tau_f 和 tau_g 分别控制动量和热量的扩散率，F_i 是浮力项。为什么要把温度单独放到 g 里而不是直接在 f 里加一个温度浓度项？因为自然对流的特征时间尺度由热扩散决定，而流动的时间尺度由动量扩散决定，两者在 Prandtl 数不为 1 时并不一致。双分布函数允许两个松弛时间各自独立设置，这样 Pr 只是一个派生参数，而不需要强制压缩时间步。

从物理图像上看，f 的零阶矩是密度，一阶矩是动量；g 的零阶矩是温度。温度场通过宏观速度被 f 输运，反过来温度又通过浮力项影响 f 的碰撞。这个单向耦合是 Boussinesq 近似的核心：密度只在浮力项里随温度线性变化，其他位置仍按不可压处理。因此双分布函数天然适配低马赫数自然对流，不需要显式求解泊松方程，也不必在每个时间步做压力修正。

2.2 D3Q19 与 D3Q7 离散速度及权重参数

三维自然对流里流场最常用 D3Q19，温度场最常用 D3Q7。D3Q19 有 19 个离散速度，速度组能恢复 Navier-Stokes 方程；D3Q7 只有 7 个方向，对温度场的对流扩散方程已经能给出二阶各向同性结果。如果追求更高精度，流场换 D3Q27，温度场换 D3Q15，但内存占用会明显上升，三维模拟里每一层方向数组都是完整的三维张量，方向数差 8 个，内存就差 8 倍。

下表是我在代码里常用的两组权重和方向组：

离散速度的具体编号需要和权重严格配合，不能只抄权重不抄方向。D3Q19 的方向可以分成三组：一个静止方向，六个面中心方向正负 x/y/z 轴，十二个棱中心方向两个坐标轴相组合。D3Q7 则只有一个静止方向加六个轴向方向。实现时最好写成独立模块，不要手写 19 次 if，用数组存方向，后续碰撞和迁移都循环数组下标，既清晰又方便换模型。

2.3 用两个松弛时间控制 Prandtl 数

在双分布函数中，运动学粘性系数 nu 由 tau_f 决定，热扩散系数 alpha 由 tau_g 决定，Prandtl 数 Pr = nu / alpha。所以在给定物理 Pr 后，代码里只需要选取 tau_f，再通过下面的公式计算 tau_g：

c_s2 = 1.0 / 3.0 nu = (tau_f - 0.5) * c_s2 alpha = nu / Pr tau_g = alpha / c_s2 + 0.5

参数说明：c_s2 是格子声速平方，LBM 的不可压极限依赖它；tau_f 必须大于 0.5，否则 nu 为负，模拟必然发散；tau_g 同理。实际调试中，我会先固定 tau_f = 0.8 这样偏稳的值，再按 Pr 算出 tau_g。如果 Pr 很小，tau_g 会很接近 0.5，此时每一步碰撞都会放大温度场的高频扰动，需要同步减小每步速度或加密网格，而不是硬抗。

提示：很多三维自然对流的发散事故都出在 Pr 数失真上。检查代码前先打印 nu 和 alpha 的实际值，确认它们与设定的 Pr 一致，再排查边界条件。

3. 用 Python 搭建三维双分布函数自然对流求解器

3.1 初始化网格与分布函数数组

写代码前要把无量纲单位定清楚。三维自然对流通常用 Rayleigh 数 Ra 和 Prandtl 数 Pr 描述，这两个量在格子单位里通过对流项和扩散项的比例关系体现。网格尺寸 L 以格子数为单位，Ra 的表达式为 g_beta * deltaT * L^3 / (nu * alpha)，其中 g_beta 是重力与热膨胀系数的乘积，需要从目标 Ra 反推。

初始化代码我一般这样写：

import numpy as np Lx, Ly, Lz = 48, 48, 48 tau_f = 0.8 Pr = 0.71 c_s2 = 1.0 / 3.0 nu = (tau_f - 0.5) * c_s2 alpha = nu / Pr tau_g = alpha / c_s2 + 0.5 # 重力方向设为 x 轴，x=0 为高温壁面，x=Lx-1 为低温壁面 # Ra 目标值，比如 1e4 Ra = 1e4 g_beta = Ra * alpha * nu / (Lx**3) # 分配分布函数数组 Qf, Qg = 19, 7 f = np.zeros((Qf, Lx, Ly, Lz), dtype=np.float64) g = np.zeros((Qg, Lx, Ly, Lz), dtype=np.float64)

这里有一个容易踩的坑：g_beta 的表达式依赖 Lx 的三次方，网格一变就要重新算。如果直接从物理加速度换算成格子加速度，往往会忽略空间步长和时间步长的缩比，导致浮力远大于数值稳定极限。我总是先把 nu 和 alpha 定下来，再反推 g_beta，这样 Ra 不会因网格加密而漂移。

3.2 碰撞模块：平衡态分布函数与浮力

碰撞步骤要同时更新 f 和 g。平衡态分布函数是碰撞的锚点，必须保证矩匹配。f_eq 需要保留速度的二阶项，g_eq 只保留一阶项，这是双分布函数处理自然对流的常见做法。

def equilibrium_f(rho, u, v, w): feq = np.zeros_like(f) for i in range(Qf): cu = ex[i]*u + ey[i]*v + ez[i]*w u2 = u*u + v*v + w*w feq[i] = wf[i] * rho * (1.0 + cu / c_s2 + 0.5 * cu**2 / c_s2**2 - 0.5 * u2 / c_s2) return feq def equilibrium_g(T, u, v, w): geq = np.zeros_like(T) for i in range(Qg): cu = exg[i]*u + eyg[i]*v + ezg[i]*w geq[i] = wg[i] * T * (1.0 + cu / c_s2) return geq

g_eq 只保留到一阶项并不会导致明显的各向异性误差，因为温度方程是标量扩散，二阶项只影响大温差可压缩效应，在 Boussinesq 近似下可以忽略。浮力项需要加到 f 的碰撞结果里，我用 Guo 力模型实现，这样比简单地把加速度加到宏观速度上更稳：

# 浮力方向沿 x 轴正方向，由高温指向低温 F = g_beta * (T - T_cold) * rho # 每个方向的力贡献 for i in range(Qf): f[i] += (1.0 - 0.5 / tau_f) * wf[i] * F * ex[i] / c_s2

参数说明：F 的符号与重力方向、坐标系有关。如果高温壁面在 x=0，浮力应指向 x 轴正方向，那么 ex[i] 为正的方向获得正向动量。写错符号的最直接后果是产生反向对流，Nusselt 数会小于 1，表现为完全没有传热增强。

3.3 迁移与热边界：三维数组的滚动操作

迁移步骤在离散网格上就是把每个方向的分布函数沿速度方向搬移一位。用 numpy 的 roll 可以实现，但要注意 roll 是周期搬移，对于方腔边界需要事后覆盖。

def stream_f(f): for i in range(Qf): f[i] = np.roll(f[i], shift=(ex[i], ey[i], ez[i]), axis=(0, 1, 2)) return f def stream_g(g): for i in range(Qg): g[i] = np.roll(g[i], shift=(exg[i], eyg[i], ezg[i]), axis=(0, 1, 2)) return g

对于三维自然对流的常见边界设置：x 方向为热壁面和冷壁面，y、z 方向为绝热或周期边界。热壁面用等温边界，直接把壁面上的温度分布函数覆盖为平衡态：

# x=0 为高温壁面 T[0, :, :] = T_hot g[:, 0, :, :] = equilibrium_g(T_hot, u[0,:,:], v[0,:,:], w[0,:,:]) # x=Lx-1 为低温壁面 T[Lx-1, :, :] = T_cold g[:, Lx-1, :, :] = equilibrium_g(T_cold, u[Lx-1,:,:], v[Lx-1,:,:], w[Lx-1,:,:])

速度边界要更严格一些。固定在壁面上速度为零，用反弹格式处理 f：迁移后，把从壁面外弹回的未知方向直接替换成相反方向的分布函数。对三维方腔，这等价于在 y、z 方向使用周期边界，在 x 方向使用标准反弹。注意热壁面的 f 也需要反弹，而 g 使用平衡态覆盖，两者不能混用。

4. 自然对流模拟的参数标定与稳定性调试

4.1 从物理参数到格子参数：Ra 数和 Ma 数约束

三维双分布函数模拟的稳定性不只是 tau 的取值问题，更关键的是格子 Mach 数。LBM 依赖低马赫数假设，格子单位下最大宏观速度通常要控制在 0.1 以下。三维自然对流在高 Ra 数下会产生强烈的羽流，如果 Ra 数设定过高而网格分辨率不足，局部速度会轻松突破 0.1，随后数值不稳定性会从羽流顶端扩散开。

因此我每跑一个新 Ra 数都会先做一个短时间预跑，打印最大速度：

max_u = np.max(np.sqrt(u**2 + v**2 + w**2)) if max_u > 0.1: print(f"Ma too high: {max_u:.3f}")

如果超限，第一选择是增大 tau_f，把 nu 调大，但这会改变有效 Ra 数，所以要用 3.1 节的反推公式重新计算 g_beta。第二选择是加密网格，L 变大后同样 Ra 数对应的 g_beta 变小，最大速度也会下降。第三选择是调整初场，去掉不必要的随机扰动幅度。一般顺序是网格、tau、扰动，不要一开始就缩小每一步的物理时间，那会直接破坏无量纲对应关系。

4.2 初始温度扰动与收敛判据

自然对流问题虽然控制方程是确定的，但静止热传导态在 Ra 超过临界值后是不稳定的。如果初始温度场只做线性层结，数值噪声不足以触发对流，模拟会停留在一个伪稳态。常见的做法是在初始温度场上加小振幅随机扰动，幅度取 1e-4 到 1e-6 量级：

T = T_cold + (T_hot - T_cold) * x / Lx T += 1e-5 * np.random.rand(Lx, Ly, Lz)

扰动的具体分布不会影响最终稳态，但幅度过大会直接给温度场注入非物理能量，导致一开始就发散；幅度过小则要等几万步才能看到对流发展。我一般同时监测全局平均 Nusselt 数和温度场方差，当两者在连续 1000 步内变化小于 1e-6 时认为达到稳态。这里要注意，LBM 是非定常格式，瞬时 Nu 有波动，只看相邻两步的残差没有意义，要看窗口内的平均值。

4.3 模拟发散时该检查哪些量

发散定位要按照从宏观到微观的顺序。第一步在所有节点上检查密度和温度是否为 NaN 或 Inf，第二步看 NaN 出现的坐标，第三步根据坐标判断是边界问题还是体区域问题。

if np.isnan(T).any() or np.isinf(u).any(): bad_coords = np.argwhere(np.isnan(T)) print("NaN at", bad_coords[:5], "step", step)

如果 NaN 集中在热壁面，优先检查 g 在壁面上的覆盖逻辑；如果集中在角点，很可能是两个方向的反弹格式叠加冲突；如果集中在中心区域，那就是浮力项或松弛时间的问题。另一个常见现象是温度整体漂移，低温壁温度越来越高，这是因为 g 的平衡态里没有包含源项，壁面覆盖频率不足导致能量守恒被破坏。这种问题不会立即 NaN，但 Nusselt 数会持续下降。

下表是我在三维自然对流调试中常用的经验参数范围：

5. 用 Nusselt 数验证三维双分布函数模拟结果

5.1 Nusselt 数的格子单位计算

自然对流模拟是否正确，最终要用 Nusselt 数来验证。Nu 表示对流换热与导热之比，在三维热壁面上可以写成局部热通量与纯导热热通量的比值。我通常在每个时间步后用热壁面内部的两个格点计算温度梯度，避免直接用壁面一阶差分造成过大误差：

gradT = (-3.0 * T[0, :, :] + 4.0 * T[1, :, :] - T[2, :, :]) / (2.0 * dx) Nu_local = gradT * Lx / (T_hot - T_cold) Nu_avg = np.mean(Nu_local)

从局部热通量计算而不是从温度场后处理插值，可以直接复用边界上的分布函数数据。对于三维方腔，这个 Nu 值应该与基准解在 5% 以内吻合。如果网格从 32 加密到 48，Nu 变化超过 5%，就要回查 g 的壁面覆盖和速度边界的反弹实现。

5.2 常见的三个验证偏差原因

第一个原因是热壁面的方向与重力方向不一致，导致浮力方向错误，Nu 会略小于 1。第二个原因是 g 的松弛时间格式不对，Pr 数失真，Nu 偏差会随 Ra 增大而增大。第三个原因是统计窗口太短，Nu 还在振荡期就被记录下来。我习惯把模拟分成预热和统计两阶段，前一万步只演化不统计，之后每 100 步记录一次 Nu，取窗口平均值。

5.3 用 GPU 数组替代 numpy 滚动

当网格达到 128^3 时，Python 里逐方向循环滚动的耗时很可观。我一般把 f 和 g 放到 GPU 上，用 numba 的 CUDA 内核实现碰撞，迁移用数组切片而不是 roll，因为 roll 会触发整包内存拷贝。这里有一个具体的优化技巧：把 g 的七个方向单独存为连续的小数组，可以显著提高缓存命中率。双分布函数的内存开销约为 (19+7) 个单精度三维数组，128^3 时大约只有 100 MB，单卡就能放下，瓶颈反而在迁移的访存模式。如果让我重新搭一套，我会先用 48^3 网格把边界和 Nu 验证跑通，再放大网格，比一上来就 128^3 省很多调试时间。

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