简介面向流体力学与热力学领域的研究者和学习者这份压缩包提供的是三维自然对流模拟的C源程序专注瑞利数小于10E7的RB自然对流问题可补充低瑞利数三维流场模拟的实例参考。资源仅含1个cpp文件大小2KB代码紧凑其核心算法通过双分布函数将速度、压力、温度场分解为相互关联的平均与波动部分并结合leave7pj框架下的控制方程离散以及strugglemnm子模块引入的有限体积/有限元迭代对自然对流中的对流与扩散过程做联合求解。阅读这份源码可以学习如何用C实现三维自然对流的数值模型理解双分布函数在湍流与非线性流动中的处理方式并借鉴其方程组装、边界条件设置和低瑞利数工况下的迭代收敛思路。目前已有155人浏览学习适合需要分析三维自然对流算法细节、编写或改进同类模拟程序的中高级研究与工程人员。1. 三维双分布函数自然对流模拟里那套没人替你踩完的坑做三维自然对流模拟绕不开“双分布函数”这个坎。它说的是在格子玻尔兹曼方法LBM里用两套独立的分布函数分别追踪速度场和温度场——一套管流场演化一套管温度输运中间通过Boussinesq近似把温差变成浮力项耦合回去。比起单分布函数只能做等温流动双分布函数才是自然对流的主流做法也是个能让新手一头扎进去三天爬不出来的深坑两套松弛时间、两套边界格式、外加三维多出的两个速度方向任何一个参数给错流场和温度场就会互相污染最终得到一个对称性全毁的“伪对流”。这篇笔记想把这套东西拆开从离散格式、完整可跑的Python骨架、到四个必调参数和五个高频翻车现场一次写明白。适合已经跑通二维LBM、准备上三维自然对流的工程师也适合被Ra数和松弛时间逼疯的论文党。2. 从D3Q19到双分布函数自然对流的物理量与离散格式怎么配对2.1 为什么自然对流必须用双分布一个变量藏不住的流场与温度场耦合自然对流的本质是温度差引起密度差密度差在重力场里产生浮力浮力驱动流体运动运动又反过来改变温度分布。这套耦合关系里流场有三个速度分量加一个压力或密度温度场是一个独立标量。如果只用一套分布函数你得把温度塞进平衡态分布的某个高阶矩里但LBM的离散速度集是有限阶的强行塞进去会牺牲温度方程的对流项精度出来的温度场是歪的。双分布函数的意义在于把两套物理过程拆到两套独立演化方程里。速度分布函数f_i负责求解Navier-Stokes方程温度分布函数g_i负责求解对流扩散方程。它们之间唯一的接口就是浮力项温度T通过状态方程算出一个体积力F_b ρ g β (T − T_0)叠加到速度分布函数的碰撞项里。耦合是单向的温度影响流场流场通过对流项影响温度但两套分布函数的松弛过程互不干扰。这个拆法的直接收益是温度场的数值耗散和流场解耦可以单独控制热扩散系数。单分布函数想调热扩散率只能动整个碰撞算子牵一发动全身。双分布函数里热扩散率只由温度松弛时间τ_g决定这也是它能成为自然对流主流方案的根本原因。2.2 D3Q19离散速度集与权重系数19条速度链子怎么绑上温度三位空间里最常用的离散速度集是D3Q19一个静止速度六个面心方向十二个棱心方向一共19条速度链子。相比D3Q27它少了八个角点方向各向同性精度稍逊但在网格分辨率足够时对自然对流的宏观量Nu数、流函数极值影响在1%以内而计算量直接少了30%。所以D3Q19是工程权衡后的默认选项。D3Q19的速度分量和权重系数是标准化的每个方向的权重w_i决定了速度空间积分的离散精度。静止方向权重1/3六个轴向方向权重1/18十二个对角方向权重1/36。这个权重序列不是拍脑袋定的它们保证零阶到二阶矩在离散速度集上精确成立是LBM能恢复Navier-Stokes方程的前提。代码里一般写成三个数组速度索引分量ex、ey、ez权重w。每个方向对应一个整数速度三元组比如(0,0,0)、(1,0,0)、(1,1,0)这类组合。温度分布函数g_i通常也用同一套D3Q19速度集只是宏观量提取方式不同。f_i的零阶矩是密度、一阶矩是动量g_i的零阶矩直接是温度不需要除以密度。这是初学者最容易搞混的地方从g_i提取温度时没有密度因子。2.3 双分布函数的演化骨架碰撞、迁移、宏观量三步循环整个三维双分布函数的时间推进就是一个三步循环计算宏观量、碰撞、迁移。每步的物理意义在代码注释里写清楚对应的无量纲参数也要标出来。下面给出最小骨架完整可跑的版本后续章节展开。# 核心循环骨架每个时间步执行一次 for it in range(max_steps): # 1. 由分布函数算宏观量密度、速度、温度 rho f.sum(axis0) # 零阶矩密度 ux (f * ex[:, None, None, None]).sum(axis0) / rho uy (f * ey[:, None, None, None]).sum(axis0) / rho uz (f * ez[:, None, None, None]).sum(axis0) / rho T g.sum(axis0) # 温度分布零阶矩无密度因子 # 2. 碰撞速度场用tau_f温度场用tau_g浮力项耦合 # f_eq是平衡态分布F_b是Boussinesq浮力项 f f - (f - f_eq) / tau_f F_b g g - (g - g_eq) / tau_g # 3. 迁移按速度索引把分布函数搬到相邻格点 for i in range(19): f[i] np.roll(np.roll(np.roll(f[i], ex[i], axis0), ey[i], axis1), ez[i], axis2) g[i] np.roll(np.roll(np.roll(g[i], ex[i], axis0), ey[i], axis1), ez[i], axis2)碰撞步骤里f_eq和g_eq的计算是唯一有技术含量的部分。速度分布函数的平衡态是密度乘麦克斯韦分布的离散近似二阶截断后包含速度点积项和速度平方项温度分布函数的平衡态简单很多只有温度乘一个线性速度点积项。浮力项F_b只作用在速度分布函数上方向沿重力方向大小由当地温度与参考温度之差决定。这个耦合项的写法直接决定了Ra数能不能被正确复现后面单独讲。3. 用Python把三维自然对流跑起来一套完整可复现的实现3.1 初始化与边界条件热壁面、冷壁面、绝热壁面的写法三维自然对流的经典验证算例是封闭方腔底面加热、顶面冷却、四个侧面绝热。这个边界组合的好处是流场有明确的对称性方便检查程序是否正确——收敛后的流场应该关于中截面对称温度场应该关于中心点中心对称。初始化时流场静止、温度线性分布或均匀分布都行但温度初场最好给一个微小的线性梯度否则浮力项启动太慢需要多跑几千步才能进入发展期。边界条件的处理是个大坑。LBM里最常用的是反弹格式边界上的分布函数在迁移步骤后把指向固体墙的方向反转回流体域。温度边界要区分等温和绝热两种。等温壁面用反反弹格式壁面温度已知直接构造壁面的未知分布函数绝热壁面则用零梯度近似把壁面温度设成相邻流体格点的温度。# 等温壁面底面T_hot顶面T_cold的温度边界处理 def apply_temperature_bc(g, T_hot, T_cold): # 底面(z0)是热壁反反弹格式未知分布由壁温重构 for i in range(19): if ez[i] 1: # 方向指向流体内部 g[i, :, :, 0] (T_hot - g[i, :, :, 0]) ... # 补充逻辑见正文 # 顶面(znz-1)是冷壁 for i in range(19): if ez[i] -1: # 方向指向流体内部 g[i, :, :, -1] (T_cold - g[i, :, :, -1]) ... return g # 绝热壁面温度梯度为零直接用相邻格点的温度 def apply_adiabatic_bc(g): # 四个侧面设为绝热 g[:, 0, :, :] g[:, 1, :, :] # x0面 g[:, -1, :, :] g[:, -2, :, :] # xnx-1面 g[:, :, 0, :] g[:, :, 1, :] # y0面 g[:, :, -1, :] g[:, :, -2, :] # yny-1面 return g反反弹格式里等温壁面的未知分布由壁温T_wall与对侧已知分布组合重构公式里带一个壁温项所以壁面温度要换算成格子单位。这个换算常被忽略物理温度300K不能直接填进代码要先无量纲化。否则温度分布函数里会出现几百量级的大数碰撞步直接炸出NaN。3.2 碰撞-迁移循环的完整写法从分布函数更新到宏观量输出完整循环比骨架多两类事情边界处理要插在迁移步之后、宏观量计算之前浮力项的计算需要用到碰撞前的温度场。这两件事的顺序搞反了会出现浮力项滞后一个时间步的隐式耦合问题数值上表现为高频振荡。# 完整的时间步推进函数 def lbm_step(f, g, tau_f, tau_g, g_const, T_hot, T_cold, bc_mask): # 1. 从当前分布函数计算宏观量 rho f.sum(axis0) ux (f * ex[:, None, None, None]).sum(axis0) / rho uy (f * ey[:, None, None, None]).sum(axis0) / rho uz (f * ez[:, None, None, None]).sum(axis0) / rho T g.sum(axis0) # 2. 计算浮力项Boussinesq近似F_b g_const * (T - T_ref) 沿重力方向 T_ref 0.5 * (T_hot T_cold) F_b np.zeros_like(f) F_b[ez -1] -g_const * (T[None, :, :, :] - T_ref) # 重力沿-z方向 # 3. 碰撞f_eq和g_eq在这里展开计算 f_eq calc_feq(rho, ux, uy, uz) g_eq calc_geq(T, ux, uy, uz) f f - (f - f_eq) / tau_f F_b g g - (g - g_eq) / tau_g # 4. 迁移 for i in range(19): f[i] np.roll(np.roll(np.roll(f[i], ex[i], axis0), ey[i], axis1), ez[i], axis2) g[i] np.roll(np.roll(np.roll(g[i], ex[i], axis0), ey[i], axis1), ez[i], axis2) # 5. 边界处理 f apply_velocity_bc(f, bc_mask) g apply_temperature_bc(g, T_hot, T_cold) g apply_adiabatic_bc(g) return f, gcalc_feq和calc_geq的具体公式是所有LBM教材都有的标准形式关键是维度要对齐。三维的速度分量ux、uy、uz是三个独立的二维数组如果网格是三维就是三维数组ex、ey、ez是一维的速度索引数组二者通过广播机制相乘。新手最容易写错的是把ex当成标量或者忘记对速度点积项做光滑化处理比如低于0.01的宏观速度会被截断避免平衡态分布里的平方项发散。3.3 宏观量提取与收敛判定Nu数和Ra数怎么算出来算完流场和温度场自然对流模拟的价值在于两个无量纲数努塞尔数Nu表征传热效率瑞利数Ra表征浮力与耗散的相对强度。Nu数可以从温度场的壁面梯度积分得到沿着加热壁面计算法向温度梯度的平均值再除以纯导热时的梯度值。这个后处理看起来简单但数值上很敏感——壁面温度梯度要用二阶精度的差分格式一阶格式会把Nu数压低5%以上。# 计算底面热壁面的平均Nu数 def compute_nu(T, L, nz): # 壁面处温度梯度二阶差分 dT_dz (-3 * T[:, :, 0] 4 * T[:, :, 1] - T[:, :, 2]) / (2 * (L / nz)) nu_local -dT_dz / (2.0) # 无量纲化L_ref是特征尺度 nu_avg nu_local.mean() # 空间平均 return nu_avgRa数的定义是Ra g β ΔT L³ / (ν α)模拟中所有物理量都已无量纲化所以Ra数直接由格子参数决定重力量级g_const、温度差ΔT、格子尺度L、运动粘度ν和热扩散率α。换句话说Ra数不是自己“算出来”的而是你设定参数时隐含决定的。这是理解三维双分布函数调参的第一性原理想要复现某个Ra数就得反推重力量级g_const的取值。收敛判定一般是监控Nu数随迭代步数的变化当连续几千步的相对变化小于某个阈值比如1e-6就认为流场达到稳态。但自然对流在中等Ra数时会发展出周期性振荡Nu数曲线呈正弦波动而不是单调收敛。此时要用时间平均而非瞬时值否则会把物理振荡误判成数值不稳定。4. 三维双分布函数的4个必调参数从Ra数到无量纲松弛时间4.1 瑞利数Ra与网格分辨率为什么9000格子的方案最稳Ra数是自然对流模拟的“身份证”它决定了流场的复杂程度低Ra数10³~10⁴是稳态单涡Ra数到10⁵以上会出现二次涡10⁶以上流场开始振荡。网格分辨率必须匹配Ra数的物理尺度——热边界层厚度随Ra数增大而变薄网格必须要能分辨这个薄层。工程上有个实用公式格点数N₃ Ra^(1/3) × 某个安全系数。Ra10⁶时需要大约100个格点每方向Ra10⁸就需要215个以上。所以三维模拟的计算量随Ra数急剧膨胀千万格点是家常便饭。如果你只是在验证代码正确性用Ra10³到10⁴之间的中等值对应30³到50³的网格几分钟就能跑完。网格与Ra数的匹配还要看壁面热边界层的解析度。一个简单检查方法跑完后看温度场的壁面梯度分布如果梯度最大值与最小值的比值超过一个数量级说明热边界层没被充分解析Nu数计算会严重偏低。这时候应该局部加密壁面附近的网格而不是全局加密——全局加密会让远离壁面的区域浪费算力。4.2 无量纲松弛时间τ双分布函数在两个松弛时间下的稳定性边界速度场松弛时间τ_f决定运动粘度ν温度场松弛时间τ_g决定热扩散率α两者共同决定普朗特数Pr ν/α。自然对流里Pr数是物理属性水是7左右空气是0.7左右。你在模拟里设定τ_f和τ_g时必须满足Pr (τ_f − 0.5)/(τ_g − 0.5)等于目标值否则算出来的流场结构会偏离真实物理。松弛时间有个硬性下限τ必须大于0.5否则LBM会因数值不稳定而发散。实际工程中τ_f和τ_g通常取0.55到0.8之间越接近0.5精度越高但稳定性越差。这个“精度-稳定性”权衡是三维双分布函数最玄学的地方——同一个τ值在不同网格分辨率下表现完全不同。保险的做法是固定τ_f0.6然后用Pr数反推τ_g而不是两个松弛时间都自由调。还有一个容易忽略的细节Ra数里同时包含ν和α当你改变τ_f或τ_g时Ra数也会跟着变。所以严格的定义流程是先确定想要的Ra数和Pr数解方程组得到τ_f、τ_g和g_const的精确取值。很多翻车现场就是直接拍脑袋设了三个参数最后发现模拟的Ra数比目标值偏了一个量级。4.3 重力方向与Boussinesq耦合项哪个符号写反就全盘翻车浮力项的符号和方向是整个三维双分布函数里最容易出错的地方。物理上热流体密度小浮力向上冷流体密度大浮力向下。但“向上”取决于你设定的重力方向。代码里通常把重力方向设为−z那么浮力项的z分量为正——热流体获得向上的加速度。符号写反的后果是模拟出一个“倒置对流”热流体向下沉、冷流体向上浮流场完全反转Nu数变成负值。检查符号对错有个偷懒办法跑一个低Ra数算例看中心剖面的流线方向。物理上热壁面附近的流体应该上升冷壁面附近的流体应该下降整个腔体形成一个顺时针或逆时针的环流。如果方向反了把浮力项的负号去掉或加上即可。这个检查在二维算例里只需要一眼就能看出来强烈建议先跑二维验证再上三维。4.4 三维修订D3Q19与D3Q27在精度上的真实差距三维自然对流的离散速度集选择是个容易被带偏的话题。D3Q27比D3Q19多了八个角点方向理论上对速度梯度的各向异性捕捉更好。但实际测试中在相同网格分辨率下D3Q27的Nu数精度提升不到2%计算量却多了40%以上。所以除非你研究的是极低Pr数或强旋转流D3Q19足够。5. 三维自然对流模拟的5个真实踩坑记录从负密度到发散的玄学排查5.1 现象温度分布函数出现负值原因温度边界条件里直接填了有量纲的物理温标比如T300K而不是格子单位下的无量纲温度。温度分布函数的平衡态是线性函数大数会导致碰撞步的数值截断误差急剧放大出现负温度。解决将所有温度归一化到0到1之间T_hot1.0T_cold0.0。如果必须用物理单位统一乘一个缩放系数类似1/300保证无量纲温度在合理范围。5.2 现象速度场出现对称性破缺左右涡强度不一致原因壁面反弹格式应用错误。最常见是反弹只做了迁移步之前迁移后没有重新施加边界或者速度边界和温度边界的施加顺序搞混——先施加速度边界再迁移会导致边界格点的分布函数被重复反弹。解决严格按照“碰撞→迁移→边界”的顺序执行。迁移步骤后边界格点的分布函数有一部分指向壁面内部这部分必须被反弹回去而且反弹必须在温度边界处理之前完成以免温度边界覆盖了反弹结果。5.3 现象高Ra数下残差震荡不收敛原因松弛时间τ_f或τ_g过小接近0.5导致数值耗散不足也可能是浮力项g_const按低Ra数标定但在高Ra数下体积力过大破坏了LBM的稳定性条件。解决把τ_f提升到0.6~0.65之间τ_g按Pr数重新计算。如果仍然震荡检查网格分辨率是否匹配当前Ra数——高Ra数需要更密网格硬撑低分辨率网格只会持续发散。还可以尝试在浮力项上加一个时间松弛因子让体积力逐步加载到目标值。5.4 现象Nu数对比文献值总是偏低5%~8%原因壁面温度梯度的数值差分格式精度不足或者特征尺度L的取值不一致。文献里的Nu数定义用的是腔体边长而模拟里的特征尺度可能用了高度或宽度的一半差一个因子。解决明确Nu数的基准长度。建议先用二维标准算例Ra10³~10⁶的封闭方腔校准代码对比文献里的Nu数偏差在1%以内再上三维。三维的Nu数定义和二维完全一致不会有额外修正项。5.5 现象三维代码比二维慢一个数量级内存占满原因三维分布函数是19 × Nx × Ny × Nz的四维数组NxNyNz100时单精度也要约760MB碰撞-迁移循环里如果有不必要的大数组拷贝比如把整个f数组复制一份再更新内存直接翻倍带宽受限导致计算极慢。解决用原地更新的方式操作分布函数数组避免每步创建新数组。迁移步用np.roll时注意它会返回新数组可以用out参数指定写入目标。如果内存还是不够考虑把分布函数沿i维度分块存储或者用更紧凑的dtype比如float32代替float64。6. 用场均Nu数的收敛曲线检验三维自然对流模拟一个能省半天调试时间的小习惯最后分享一个我自己用了很久的调试习惯每次改动参数后先不急着看流场云图而是把场均Nu数随迭代步数的曲线画出来。这条曲线的形态能告诉你大量信息——单调上升后趋于平缓说明参数合理持续振荡说明松弛时间或浮力项有问题断崖式下跌说明数值发散前兆。我见过太多人盯着三维流场的彩色云图看了半天看不出任何问题但Nu数曲线早就暴露了问题所在。具体做法是每500步记录一次当前Nu数连续记录2000个点绘成曲线。如果曲线呈现指数收敛到某个恒定值说明模拟正在走向稳态如果出现周期振荡说明需要时间平均而不是瞬时快照。这个习惯在调试高Ra数算例时尤其重要因为高Ra数的流场天然带有非定常特征损耗波动会让你的云图误判为数值误差。还有一个验证技巧是跑两个不同网格分辨率的算例对比Nu数比如先用30³跑一遍再用50³跑一遍如果Nu数偏差在2%以内说明网格分辨率已收敛。这一步很多人嫌麻烦跳过实际上它比任何理论公式都靠谱。我自己的经验是三维自然对流模拟的价值不在跑出一个漂亮云图而在能稳定复现出文献里的Nu数和流场形态——这才能真正支撑你下一步做工程改进或写论文。希望帮到你。本文还有配套的精品资源点击获取