三维自然对流模拟的双分布函数：D3Q19动量与D3Q7温度耦合实现

简介一份面向流体力学与数值传热方向的研究者、工程技术人员及高年级本科生的三维自然对流模拟C源码适用于瑞利数小于一千万的RB自然对流问题分析与教学参考。程序围绕leave7pj与strugglemnm两个关键框架搭建通过双分布函数将流场分解为平均与波动部分以捕捉对流与扩散相互影响并采用迭代法求解连续性、动量及能量方程。压缩包内共有一个文件即三维双分布函数.cpp包体大小仅2KB代码精简、结构清晰便于逐行阅读。读者可从中掌握自然对流数值模拟的基本框架、双分布函数建模思路及低瑞利数条件下的求解策略也可作为教学示例或科研对比的轻量级参考实现。现有155人学习浏览适合对三维自然对流计算感兴趣的中高级学习者。1. 三维自然对流模拟的双分布函数为什么一个分布函数带不动温度你照着二维热对流的教程把代码改成三维换上 D3Q19 后加热底面温度场却纹丝不动Nusselt 数一直是 1。别怀疑是浮力没加上——问题出在“只有一套分布函数”的等温 LBM 根本无法表达温度输运。标题里的“三维双分布函数”正是为 natural convection 这类浮力驱动流动准备的一套分布函数演化速度与压力另一套演化温度浮力项把两者耦合起来。这篇文章会带你拆开这套耦合逻辑、落实到可运行的 Python 代码再把尺度参数和边界坑一个个讲透。适合正在写三维 LBM 还差点火候的工程师也适合刚入热流仿真的研究新手。2. 双分布函数模型怎么搭D3Q19 动量格子与 D3Q7 温度格子的选择2.1 两个分布函数的分工与 Boussinesq 耦合三维自然对流不像等温流动那样只有一个未知场。流体运动由速度与压力决定温度场上还叠了一个浮力源项温差改变局部密度密度差产生体积力体积力驱动流动流动再反过来改变温度分布。用格子 Boltzmann 方法模拟这个闭环最常用的方案就是双分布函数——动量场用 f 分布温度场用 g 分布每组分布各自完成碰撞和流更新只在浮力项上交换信息。这个耦合之所以成立靠的是 Boussinesq 近似除密度随温度线性变化外其余物性当作常数。浮力项写成 F ρ₀ g β (T − T_ref) e_ze_z 是重力方向的单位向量。这里的 β 是热膨胀系数T_ref 一般取冷热壁平均温度。在 LBM 里这个力不是直接加到宏观速度上而是通过离散力项进到分布函数的碰撞中后面第 3 章会给出具体形式。工程上要注意Boussinesq 只适用于小温差场景。电子散热里元件表面与空气温差动辄几十度严格说已经越界需要用变物性模型或低马赫数可压求解但这不影响方法本身的验证价值绝大多数论文在做算法对比时仍用 Boussinesq。三维双分布函数的真正优势在于温度场的演化独立于动量场格子模型能按各自需求裁剪不必强行共用一套 D3Q27。2.2 离散速度和权重一边 19 个方向一边 7 个方向动量场要恢复完整的 Navier-Stokes 方程离散速度必须满足中心对称和足够的各向同性。三维里 D3Q19 是平衡点19 个方向包括中心点、6 个轴向和 12 个面对角质量上能覆盖三阶各向同性需求内存又比 D3Q27 省 30% 左右。D3Q15 虽然更省但剪切应力的各向同性误差会在转捩流动里被放大做自然对流不太划算。温度场就简单得多。温度方程是对流扩散方程宏观守恒量只有零阶矩温度 T和一阶矩对流热流 uT不涉及应力张量所以格子的离散方向越多越浪费。D3Q7 只保留中心点和 ±x、±y、±z 六个方向刚好把对流速度的三个分量耦合进去。温度场的格子声速平方 c_sT² 1/4这个值会在下面参数换算时持续用到。为了让你后面写代码时不至于对着纸面数方向我习惯直接让数组替我们算import numpy as np # D3Q19动量场19 个方向 c19 np.array([ [ 0, 0, 0], [ 1, 0, 0], [-1, 0, 0], [ 0, 1, 0], [ 0,-1, 0], [ 0, 0, 1], [ 0, 0,-1], [ 1, 1, 0], [-1,-1, 0], [ 1,-1, 0], [-1, 1, 0], [ 1, 0, 1], [-1, 0,-1], [ 1, 0,-1], [-1, 0, 1], [ 0, 1, 1], [ 0,-1,-1], [ 0, 1,-1], [ 0,-1, 1] ], dtypenp.float64) w19 np.array([1/3] [1/18]*6 [1/36]*12) # D3Q7温度场7 个方向 c7 np.array([ [0, 0, 0], [1, 0, 0], [-1, 0, 0], [0, 1, 0], [0,-1, 0], [0, 0, 1], [0, 0,-1] ], dtypenp.float64) w7 np.array([1/4] [1/8]*6)c19 的布局顺序不是随便排的第 1、2 个方向互为相反方向第 3、4 个也是以此类推到最后 12 个面对角方向两两配对。这个顺序直接决定了反弹边界里“反方向索引”的写法后面避坑章还会提到。权重之和都为 1这一点检查代码时务必验证一下任何权重数组写错都不会立刻报错但会给宏观量带入无法解释的偏差。2.3 松弛时间换算热扩散系数必须用对格子声速双分布函数最重要的是把松弛时间和物理输运系数对齐。动量场松弛时间 τ_f 与运动粘度 ν 的关系是 ν c_s² (τ_f − 1/2)动量格子 c_s² 1/3。温度场同理α c_sT² (τ_g − 1/2)但 c_sT² 是 1/4 而不是 1/3这是最容易抄错的地方。无量纲数 Pr 和 Ra 在格子单位下直接构造Pr ν / αRa g β ΔT L³ / (ν α)LBM 实践里的常见做法是设 β 1、ΔT 1、L nz一个方向上的格点数把重力加速度参数反解出来tau_f 0.6 nu (tau_f - 0.5) / 3.0 # 动量格子 c_s^2 1/3 alpha nu / Pr # 由 Pr 反推热扩散系数 tau_g 4.0 * alpha 0.5 # 温度格子 c_sT^2 1/4 g_acc Ra * nu * alpha / (nz ** 3)tau_f 选 0.6 属于中庸值既避开 0.5 附近的压缩性误差也不会让扩散过强拖慢收敛。tau_g 由 Pr 决定不要手工指定。我见过有人把 tau_g 也写 0.6结果 Pr 变成了 1Nu 和基准解差了 40%查半天才发现是松弛时间没按格子声速换算。3. 用 Python 跑通三维自然对流初始化、主循环与热边界3.1 初始化平衡态分布与线性温度场三维自然对流的经典设置是底部热壁、顶部冷壁四周绝热。模拟域取 nz 个格子沿 z 方向热壁在 z0冷壁在 znz−1。初始温度场给一个线性剖面比起全常值场更接近稳态解能显著减少前期瞬态振荡。初始速度全部为零这样平衡态分布函数可以简化但代码里我习惯把一阶速度项也写上便于以后从续算字段恢复计算。nx ny nz 32 rho0 1.0 u np.zeros((nx, ny, nz, 3)) T np.zeros((nx, ny, nz)) for k in range(nz): T[:, :, k] 1.0 - k / (nz - 1) # z0 热znz-1 冷 f np.zeros((19, nx, ny, nz)) g np.zeros((7, nx, ny, nz)) for i in range(19): cu c19[i, 0] * u[..., 0] c19[i, 1] * u[..., 1] c19[i, 2] * u[..., 2] f[i] w19[i] * rho0 * (1.0 3.0 * cu) for i in range(7): cu c7[i, 0] * u[..., 0] c7[i, 1] * u[..., 1] c7[i, 2] * u[..., 2] g[i] w7[i] * T * (1.0 4.0 * cu)f 数组的形状是 (19, nx, ny, nz) 而不是 (nx, ny, nz, 19)这个细节在 Python 里影响不大但在 C/C 里按轴遍历时前者对f[i]整面切片更友好。后续的宏观量求和np.sum(f, axis0)也能直接得到密度场。温度场的平衡态带 4.0 倍速度项对应 c_sT² 1/4这里的系数一旦写成 3.0温度场就输运不动了。3.2 主循环碰撞、流更新和浮力的 Guo 格式主循环的顺序建议是流更新 → 算宏观量 → 计算浮力并碰撞 → 再算带外力修正的宏观量 → 热碰撞 → 边界处理。先流后撞是标准的“流-撞”顺序在瞬态计算中时间语义明确。我不建议把碰撞写在流更新前面两种写法最终都能收敛但后者的时间层混在一起续算和边界处理都容易出偏差。total_steps 50000 T_ref 0.5 for step in range(total_steps): # 1. 流更新用 np.roll 做示意真实代码要把边界方向改成反弹 for i in range(19): f[i] np.roll(np.roll(np.roll(f[i], c19[i, 0], axis0), c19[i, 1], axis1), c19[i, 2], axis2) for i in range(7): g[i] np.roll(np.roll(np.roll(g[i], c7[i, 0], axis0), c7[i, 1], axis1), c7[i, 2], axis2) # 2. 宏观量 rho f.sum(axis0) u np.einsum(i...,i...-..., c19, f) / rho T g.sum(axis0) # 3. 浮力与动量碰撞Guo 力项 Fz rho * g_acc * (T - T_ref) for i in range(19): cu (c19[i, 0] * u[..., 0] c19[i, 1] * u[..., 1] c19[i, 2] * u[..., 2]) feq w19[i] * rho * (1.0 3.0 * cu 4.5 * cu * cu - 1.5 * (u[..., 0]**2 u[..., 1]**2 u[..., 2]**2)) force (1.0 - 0.5 / tau_f) * w19[i] * ( (c19[i, 2] - u[..., 2]) 3.0 * cu * c19[i, 2]) * Fz f[i] - (f[i] - feq) / tau_f f[i] force # 4. 外力修正宏观速度 rho f.sum(axis0) u np.einsum(i...,i...-..., c19, f) / rho u[..., 2] Fz / (2.0 * rho) # 5. 温度碰撞 for i in range(7): cu (c7[i, 0] * u[..., 0] c7[i, 1] * u[..., 1] c7[i, 2] * u[..., 2]) geq w7[i] * T * (1.0 4.0 * cu) g[i] - (g[i] - geq) / tau_g # 6. 边界处理见 3.3Guo 力项中的force写法是标准离散格式第三行(c_i - u) 3 cu * c_i里的 3 就是动量格子的 1/c_s²。前置因子(1 − 0.5/τ_f)是确保多尺度展开后外力项系数正确所必需如果漏掉这个系数浮力会被人为放大低速区出现假振荡。宏观速度的修正项Fz / (2ρ)也来自同一套格式不加的话稳态度会偏移Nusselt 数对不上基准解。np.roll 是环形移动会把本该留在边界外的分布滚回对面这里只是给你一个能跑起来的骨架。真实生产代码应当在流更新前先把边界方向单独取出用反弹格式重赋值第 5 章第 4 条避坑会细说。3.3 热边界反反弹格式与绝热外推速度边界用标准反弹格式温度边界要看类型。固定温度壁面用反反弹anti-bounce-back格式绝热壁则先把相邻内点的温度外推到边界格点上再套反反弹公式。反反弹的递推关系是g_i(x_f, tδt) −g_i(x_f, t) 2 w_i T_wall其中 i 指向壁面内侧i 是它的相反方向。写成 Python 轮廓# z0 热壁T_wall 1.0 for i in range(7): if c7[i, 2] 0: # 这个方向指向域外 i_opp OPP7[i] g[i, :, :, 0] -g[i_opp, :, :, 0] 2.0 * w7[i] * 1.0 # 四周绝热先外推再按固定温度写 # T[:, :, 0] T[:, :, 1] 之类实际实现时按 y0 / yny-1 / x0 / xnx-1 四个面分别处理OPP7 是 7 个方向的反向索引表建议写成列表常量而不是每次算。绝热外推的要点是外推后的边界温度是动态的每个时间步都要重新算。有人图省事在初始化时固定了绝热壁温度结果边界从第一步开始就在漏热。4. 参数怎么设Ra、Pr、tau 与续算收敛控制4.1 无量纲数到格子单位的翻译自然对流 LBM 最常见的翻车点不在代码逻辑而在单位换算。物理世界里你要输入重力加速度 9.8 m/s²、空气的 β 和 α转换到格子单位时差一个数量级是常事。三维双分布函数的默认做法是反过来直接以格子单位构造无量纲数。设置 β1、ΔT1、Lnz从 Ra 的定义反解 g_acc Ra·ν·α / nz³。这个 g_acc 只是动量方程里的一个体积力参数不代表物理重力加速度。以 Ra1e4、Pr0.71、nz32、τ_f0.6 为例的完整参数组参数公式示例值ν(τ_f − 0.5)/30.0333αν/Pr0.0469τ_g4α 0.50.6876g_accRa·ν·α / nz³4.76e-4初始最大速度0—这套参数能让 Ra1e4 的方腔对流在 3 万步内收敛到稳定解。如果你改 Ra只动 g_acc其他按表重算一遍。写代码时把这些参数做成字典或 dataclass避免在十几个函数里散落硬编码。4.2 松弛时间 tau从稳定性到对流强度τ_f 的可选范围理论上只要大于 0.5 就满足 LBM 的数值稳定性下限但工程上 0.5 到 0.62 才是舒适区。τ_f 越接近 0.5数值声速效应越强高 Ra 时容易出现高频振荡τ_f 偏大则粘性耗散增强Nu 会偏低。我用过最稳的组合是 τ_f0.6τ_g 按 Pr 算出来落在 0.65~0.8 区间。一个常见的误解是把 τ_f 和 τ_g 都加到很大的值来压振荡。这能延迟 NaN 出现但也把对流强度削弱了算出来的流型可能从分叉态退化成对称稳态。如果高 Ra 下不稳定正确的处理顺序是提高网格分辨率 → 让 τ_f 往 0.58 靠近并配合续算 → 最后才考虑调大 τ_f。4.3 收敛判据与续算策略三维自然对流是否收敛不能只看某一点的速度要用全局量。平均 Nusselt 数是公认的指标def avg_nusselt(u, T, alpha, nz): dTdz np.gradient(T, axis2) q u[..., 2] * T - alpha * dTdz return q.mean() * nz / alpha每 100 步算一次当前 10 次平均值的相对变化小于 1e-4 就认为进入稳态。注意这里的热流 q 包含对流项 u_z·T 和扩散项 α·∂T/∂z两项都要算。只看扩散项会把 Nu 低估不少。针对 Ra 从 1e3 往 1e6 推的情况我一般用续算先在 1e3 跑通把 f、g、u、T 存成 npy 文件再加 10 倍 Ra 重新启动。这算是在为高 Ra 买保险能省下反复追发散初场的几天时间。如果你一开始就往 Ra1e6 怼大概率会撞上第 5 章要说的 NaN 坑。5. 三维自然对流模拟避坑五个让代码翻车的细节5.1 浮力项加错位置温度场白算现象程序能跑、速度场有动静但温度分布始终接近纯导热Nu 在 1.2 附近上不去。原因把 Fz 直接加到了宏观速度 u_z 上没有通过分布函数力项进碰撞。外力只改宏观量等于在每一时间步做了一个瞬时的速度增量多尺度展开恢复不了连续方程里的动量源项。解决按 3.2 的 Guo 格式把力写进分布函数并在宏观速度里补 Fz/(2ρ)。这是双分布函数耦合浮力唯一推荐的做法不要自己发明“更简单”的加力方式。5.2 D3Q7 的格子声速用错Pr 和 Nu 一起崩现象程序完全稳定温度场形状也对但 Nu 和基准解稳定偏差 15%~25%。原因把 D3Q19 的 c_s²1/3 惯性思维带到温度格子上τ_g 用了 3α0.5导致热扩散率偏高、Pr 偏低。解决连算三遍 c_sT² Σ w_i c_i² / 3 1/4把 τ_g 4α 0.5 写进公共头部并在参数表里打印出来。任何格子改动后先验证 α 的取值用纯导热问题跑 100 步对解析解。5.3 Ra 调太高程序直接 NaN现象初始化正常跑到几百步温度或密度出现 NaN且位置不固定。原因初始温度场给予全域线性梯度浮力初值在边界层处产生局部高速格子单位下的马赫数超过 0.1BGK 模型的密度扰动剧烈放大。解决分阶段续算Ra 每次只加 10 倍且用上一阶段末的完整字段初始化。如果仍发散把 τ_f 从 0.6 调到 0.65再不行就加网格数。不要靠降 Pr 掩盖问题。5.4 反弹边界写成滑移边界现象壁面附近的速度剖面在壁面处不为零甚至沿壁面出现抛物线状分布。原因np.roll 的周期性回卷自带绕环边界如果你只在 z0 和 znz−1 两个面上做反弹、而把另外四个面留给 np.roll就等于给流动开了四条环形通道也可能是反弹赋值时边界位置比流体节点偏移了半格。解决显式处理六个面的出域方向。每个方向 i 找出 i_opp让边界格点的出域分布换向补回域内。角点和棱线不要用同一个循环覆盖两次否则最后一次赋值会覆盖掉先前的反弹结果。我建议写一个通用遍历扫描所有格点对该格点的 19/7 个方向逐一看是否出域出域则反弹到目标位置省去对面循环的重复赋值。5.5 Nusselt 数对不上基准解现象流型和文献差不多局部 Nu 曲线在热壁两端有异常的尖峰或凹陷。原因宏观速度漏了外力半修正或温度梯度用了前向差分或温度边界几何位置不对齐。前向差分的截断误差在边界层里显著放大端部尖峰尤其明显。解决用中心差分计算 ∂T/∂z边界层内至少预留 3 个格点。Nu 计算同时保留对流项和扩散项然后对比 32³ 网格下的方腔基准值。如果差 2% 以内边界处理基本正确如果差 5% 以上优先怀疑第 5.2 条而不是网格。6. 用方腔基准解验收你的双分布函数代码代码跑通只是第一步。老办法是拿 de Vahl Davis 的二维方腔自然对流结果做标尺Ra1e3 时 Nu≈1.118Ra1e4 时 2.243Ra1e5 时 4.519。三维代码想验证物理部分最简单的技巧是让 z 方向只放 2 层网格两侧设为绝热这样模拟结果近似二维可以直接对上表。我把验证脚本这么组织for res in [24, 32, 48]: nu_mean run_sim(nxres, nyres, nz2, Ra1e4, Pr0.71) print(res, nu_mean)24³ 到 32³ 的 Nu 差距应在 1% 以内如果 32³ 和 24³ 差 3% 以上说明边界层分辨率不足不是代码 bug。除了平均 Nu画热壁局部 Nu 沿 x 方向的曲线更能看出问题曲线应当平滑两端有边界层尖峰如果出现锯齿或周期性波动基本就是反弹边界在角点重复赋值。网格无关性检查时不要只调网格尺寸还要按 4.1 的公式同步更新 g_acc因为 Lnz 进到 Ra 的定义里。很多人换网格后只改了数组大小g_acc 不变Nu 自然乱飘。正确做法是每次改分辨率都重新算一遍参数表并把 g_acc 打印到日志里留痕。这算我的一条血泪经验早期跑 Ra1e6直接上 128³没先做 32³ 的网格无关性检查结果局部 Nu 曲线在热壁中心出现锯齿我以为是物理振荡调了一周参数最后发现是 z 方向绝热外推在边界处少写了两个面。现在我的习惯是每次换 Ra 或换网格先跑 24³ 的短算例画局部 Nu 曲线确认平滑度再上规模。希望帮到你。本文还有配套的精品资源点击获取