简介一套专用于三维自然对流模拟的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 倍。下表是我在代码里常用的两组权重和方向组模型方向数权重使用对象D3Q1919静止 1/3面心 1/18棱心 1/36速度分布函数 fD3Q77静止 1/3轴向 1/9温度分布函数 g离散速度的具体编号需要和权重严格配合不能只抄权重不抄方向。D3Q19 的方向可以分成三组一个静止方向六个面中心方向正负 x/y/z 轴十二个棱中心方向两个坐标轴相组合。D3Q7 则只有一个静止方向加六个轴向方向。实现时最好写成独立模块不要手写 19 次 if用数组存方向后续碰撞和迁移都循环数组下标既清晰又方便换模型。2.3 用两个松弛时间控制 Prandtl 数在双分布函数中运动学粘性系数 nu 由 tau_f 决定热扩散系数 alpha 由 tau_g 决定Prandtl 数 Pr nu / alpha。所以在给定物理 Pr 后代码里只需要选取 tau_f再通过下面的公式计算 tau_gc_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 轴x0 为高温壁面xLx-1 为低温壁面 # Ra 目标值比如 1e4 Ra 1e4 g_beta Ra * alpha * nu / (Lx**3) # 分配分布函数数组 Qf, Qg 19, 7 f np.zeros((Qf, Lx, Ly, Lz), dtypenp.float64) g np.zeros((Qg, Lx, Ly, Lz), dtypenp.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 geqg_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 的符号与重力方向、坐标系有关。如果高温壁面在 x0浮力应指向 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 方向为绝热或周期边界。热壁面用等温边界直接把壁面上的温度分布函数覆盖为平衡态# x0 为高温壁面 T[0, :, :] T_hot g[:, 0, :, :] equilibrium_g(T_hot, u[0,:,:], v[0,:,:], w[0,:,:]) # xLx-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(fMa 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 数会持续下降。下表是我在三维自然对流调试中常用的经验参数范围参数经验范围调试倾向tau_f0.55 ~ 1.0小于 0.5 立即发散最大格子速度小于 0.1超过则加密网格或增大 tau_f初始温度扰动1e-4 ~ 1e-6越小越慢越大越易发散Ra 数1e3 ~ 1e6同网格下高 Ra 需要更小速度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 加密到 48Nu 变化超过 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 的七个方向单独存为连续的小数组可以显著提高缓存命中率。双分布函数的内存开销约为 (197) 个单精度三维数组128^3 时大约只有 100 MB单卡就能放下瓶颈反而在迁移的访存模式。如果让我重新搭一套我会先用 48^3 网格把边界和 Nu 验证跑通再放大网格比一上来就 128^3 省很多调试时间。本文还有配套的精品资源点击获取