MFiX源码解读：离散单元法颗粒温度计算子程序calc_thermo_des

1. 这个文件在 MFiX 里是什么角色初看 calc_thermo_des.f 的第一印象先说个结论如果你准备通读 MFiX 源码不巧走的是 DEM 这条线那 calc_thermo_des.f 早晚会出现在你的阅读清单里而且大概率是在你已经把颗粒碰撞、运动积分这些硬骨头啃完之后才会轮到它。我这个系列前面几篇主要在看整体框架和曳力计算这次算是顺着能量这条支线把离散颗粒的热量计算翻出来单独过一遍。calc_thermo_des 这个名字很直白calculate thermo of DES。DES 是 Discrete Element Simulation也就是离散单元法。整个 MFiX 在 DEM 模式下会把颗粒当作一个个独立的软球按牛顿第二定律积分位置和速度而这套运动之外还有一套能量状态——每个颗粒有自己的温度颗粒与颗粒之间会导热颗粒与周围的流体之间会对流换热高温工况下还要算辐射。负责把这些热量算清楚、然后更新颗粒温度的子程序就是 calc_thermo_des.f。这个文件在 DEM 模块里的位置很特殊。它不像 des_thermo_convection.f 或 des_thermo_radiation.f 那样只负责某一种换热机制而是把所有换热机制的贡献汇总在一起统一完成颗粒温度的状态更新再把颗粒给流体的能量源项挂到连续相方程上。换句话说真正干活的换热模型散落在各个 des_thermo 系列文件里但把它们组织起来、决定整个能量计算的主次关系和先后时序的是这个 calc_thermo_des.f。读它之前最好先搞明白一件事MFiX 的 DES 时间推进是子步制的一个流体时间步内要跑若干次 DEM 子步calc_thermo_des 在每个 DEM 子步都会被执行颗粒温度就是在这个频率下被刷新。搞清楚这个时间尺度后面读代码时对哪些变量要保存快照、哪些要清零重算才会敏感。写这篇番外之前我在自己电脑上又把这段源码翻了两遍。说实话源码里几乎没有什么注释大部分逻辑是从变量名和调用关系上反推出来的。我尽量按先看整体骨架、再看换热模型落点、最后聊坑的顺序来讲尽量让没有碰过 MFiX 的读者也能大致知道这段代码到底在算什么、算得有多准。2. 代码入手的正确姿势先理清模块归属与文件边界MFiX 的 DEM 相关源码集中在 des 相关的目录组织里命名规律比较统一。calc_thermo_des.f 从名字上属于 calc 系列和 calc_drag、calc_pic 这种并列。但想要真正读进去得先清楚它周围有哪些邻居文件因为在源代码里它的边界感和职责划分非常清晰。以我自己读源码的习惯拿到一个子程序文件第一步不是逐行读而是用文本编辑器打开后先搜索 SUBROUTINE 关键字确认这个文件里总共定义了几个子程序。calc_thermo_des.f 的情况很典型基本上一个文件里主要就一个顶层子程序内部按逻辑段切分。这个顶层子程序头几行会接收calc_energy、nphase等相关参数calc_energy 这个逻辑变量很关键它决定当前模拟算不算能量方程。如果 calc_energy 是 .FALSE.那整个子程序基本只做一件事情给颗粒温度导数数组赋零值然后直接返回。这种短路设计在 MFiX 里很常见目的是让不开启能量方程的用户跑纯等温流化床时没有任何多余开销。搞清楚文件边界之后还要大致知道它和哪些数组打交道。颗粒温度 T_S、颗粒比热容 DES_C_PS、颗粒导热系数 DES_K_TS、颗粒半径 DES_RADIUS这些都是颗粒固有属性。而流体侧需要的是能量源项数组以及网格体积、流体温度 T_G、流体密度等。calc_thermo_des 本质上就是把这两类量连接起来的中枢。读这种大型开源代码我强烈建议先用 grep 或者 ctags 把函数调用关系列出来。calc_thermo_des 直接调用的子程序不多主要包括 des_thermo_radiation、des_thermo_conduction、des_thermo_convection还有可能调 des_thermo_fluid_energy。把这几个函数在编辑器里并排打开再回来看主流程你会发现主流程其实非常短真正的工程量都在子函数里。我实际读的时候是先画了一份很粗糙的调用关系笔记——不会画太细的流程图但会把谁在什么地方被调用、更新了哪些数组记下来这样后面排查问题的时候能快速定位。3. 主循环内部在做什么逐段读懂颗粒热源项与温度推进calc_thermo_des 的主体结构是一个大循环遍历当前进程管辖的所有颗粒。每个颗粒在这里要经历四件事判断是否是活性颗粒、累计本颗粒的所有热源项、计算颗粒与流体的换热对连续相能量方程的贡献、最后用热源项更新颗粒温度。这四件事看似简单每件里面都有值得展开的说头。先说活性颗粒判断。MFiX 在并行计算时会做区域分解颗粒被分配到不同的进程管辖。这意味着循环变量不能简单地从 1 写到最大颗粒数而是要用颗粒索引数组来映射。代码里常见的手段是先从颗粒到网格的映射数组拿到当前颗粒所属的网格然后判断这个网格是否属于本进程。如果颗粒落在别的进程里就直接跳过。很多第一次读代码的人会被这块绕晕误以为代码有 bug 把颗粒漏算了其实是并行归属判断。我读的时候习惯把本进程负责的颗粒范围先打印出来确认一遍省得后面瞎猜。接下来就是热源项累计。一个颗粒的温升由三部分贡献组成颗粒与周围流体的对流传热方向取决于流体温度 T_G 与颗粒温度 T_S 的差值颗粒与颗粒之间以及颗粒与壁面之间的接触导热这个在稠密气固两相流里相当重要辐射传热通常在高温燃烧场景下才开启。每个颗粒的累计热源 Q_Source 会被清零然后三个换热模块往这个累计变量里面加各自的贡献。有的版本还会考虑化学反应产热比如煤粉热解或者气化模拟里需要加反应热但 MFiX 标准的 calc_thermo_des 并不会自己去算反应源项它只是预留了一个接口让其他模块在调用它之前先把反应热加进来。这又是一个典型的源码设计风格能量计算只管传热化学反应交给另外的模块处理。看完热源项重点就是温度推进。伪代码大致长这样IF 颗粒是活性颗粒 THEN Q_Source 0 IF calc_energy THEN 调用 des_thermo_radiation ! 累计辐射热源 调用 des_thermo_conduction ! 累计接触导热源项 调用 des_thermo_convection ! 累计对流换热源项 ENDIF ! 更新颗粒温度 T_S_new T_S_old DT_DES * Q_Source / (MASS * DES_C_PS) ! 累加流体能量方程源项 ENERGY_SOURCE(网格) ENERGY_SOURCE(网格) Q_Source ENDIF注意T_S_new T_S_old DT_DES * Q_Source / (MASS * DES_C_PS)这一行它就是颗粒温度更新的核心。分母里乘的是颗粒质量和比热容分子是净热源乘子步长时间。假如颗粒质量特别小而热源特别大一个子步下温度可能往上蹿几千度这在物理上显然是不合理的。代码里通常会有保护措施把温度限制在某个合理范围内或者强制截断到流体温度附近。MFiX 的做法在细节上可能因版本而异但显式推进 截断保护这个思路是共通的。最后说流体能量方程源项。这里有个很容易忽略的细节颗粒给流体的热量不是平均分摊到整个计算域的而是只贡献给颗粒所在的那个网格单元。能量源项数组按网格离散每个颗粒把它自己的 Q_Source 加到所在网格的累计值里。有些版本里这个源项还要求除以网格体积和流体密度因为连续相能量方程以温度为单位进行推进源项需要换算成温度变化率。算到这里calc_thermo_des 的主流程就完整了。主流程的逻辑并不复杂但代码在数组处理和单位换算上非常讲究稍微弄错一个量纲颗粒温度就会明显异常。4. 换热模型在代码里的落点对流、接触传导和辐射分别藏在哪里calc_thermo_des 虽说是solver但它自己不实现具体换热模型模型都在 des_thermo 系列子程序里。可如果你只读那些子函数而不知道主流程怎么调度会很难理解每个模型什么时候生效、什么时候被跳过。所以要把换热模型和主流程整个串起来看才能知道这几个热源项到底从哪来、精度大概在什么水平。4.1 对流换热决定传热系数的那几个关联式对流换热部分主要做一件事计算颗粒与流体之间的对流传热系数然后乘上换热面积和温差得到热流。传热系数一般不用自己瞎猜代码里会根据雷诺数 Re 和普朗特数 Pr 调用著名的努塞尔数关联式最常见的两个是 Ranz-Marshall 关联式和 Gunn 关联式都适用于球形颗粒在流体中的对流传热。Ranz-Marshall 在颗粒尺度 CFD-DEM 里非常经典Nu 2.0 0.6 * Re^0.5 * Pr^(1/3)。当颗粒静止且周围流体也静止时Nu 退化成 2.0对应纯导热极限。这个特征我记得非常清楚因为早年调试时我就犯过嘀咕为什么算出来的传热系数是个有限值而不是零后来想明白了静止流体里颗粒也会有导热这个2是物理下限。对流换热面积按球体实际表面积计算也就是 4 * pi * R^2。在这个环节代码还要处理一个重要问题颗粒表面温度用什么值去和流体温度做差。这个选择对热流计算影响很大有的实现会用上一时刻的颗粒温度有的会用更新前的旧值。MFiX 的标准做法是用当前时刻已知的颗粒温度然后显式推进不涉及隐式迭代。4.2 接触导热稠密堆积里最容易低估的一项接触导热是颗粒之间传递热量的主要机制在固定床、移动床这类颗粒紧密堆积的场景里甚至超过对流换热。des_thermo_conduction.f 的核心公式基于 Hertz 接触理论和颗粒材料导热系数。两个温度不同的颗粒接触时热量通过接触点附近的一个有限圆盘传导圆盘半径取决于接触法向重叠量、颗粒半径和材料弹性模量。代码里会循环当前颗粒的邻居列表对每个发生接触的邻居对计算一份接触导热然后累加到双方的 Q_Source。这里有个特别有意思的设计细节接触导热是对称的A 颗粒得到的热量刚好等于 B 颗粒失去的热量符号相反。所以代码里要么一个接触对同时更新两个颗粒的热源项要么把同一对接触编号做了正负处理。我读的时候专门核对过这个对称性因为如果不对称整个体系能量会有净产生整体温度漂移会暴露问题。壁面接触导热和颗粒间接触导热原理一样只是把其中一个颗粒换成壁面壁面的状态由固壁温度决定壁面温度是不变量。在流化床模拟里壁面换热对床层温度分布影响显著这部分代码不能随意跳过。4.3 辐射换热只在你需要高温时才会出现的分支辐射模型并不是 calc_thermo_des 的默认必选项只有当用户开启辐射开关时des_thermo_radiation 才会被调用。颗粒尺度辐射的模型相对简化通常基于斯特藩-玻尔兹曼定律把颗粒看作灰体用颗粒黑度和视在系数计算与周围环境的辐射热流。因为辐射是温度的四次方函数高温燃烧场景下辐射热源会迅速膨胀数值稳定性压力比前两者都大。如果发现开启辐射后颗粒温度出现振荡先查时间步长再查黑度参数最后查环境温度定义这基本是排查辐射问题的老三样。把三个模型的位置搞清楚之后你再看 calc_thermo_des 这个文件就会明白为什么主流程这么简单了。它要做的事情是在正确的时间、用正确的顺序调用正确的模型然后把结果汇总。换热模型本身的改进空间很大但如果你对整体流程没概念改进根本无从下手。5. 能量守恒与数值稳定性这段代码最容易藏 bug 的地方我读过不少开源代码calc_thermo_des 这种看起来简单但实际细节多的子程序是最容易让人翻车的。表面上一个循环加几个除法好像没什么可错的实际运行起来能量不守恒、温度跳变、结果发散这些坑一个比一个隐蔽。以下几个问题是我反复对比代码和运行结果后总结出来的几乎每条都对应过一次真实的调试经历。第一坑是累计数组没有清零。calc_thermo_des 是每个 DEM 子步都会被调用的如果某个热源累计数组在循环之前没有显式清零上一子步的残留会给下一个子步叠加出虚构热源颗粒温度会以一个离谱的斜率漂移。这种 bug 一旦出现温度场不是跳变而是缓慢爬升极具迷惑性。我的排查经验是在子程序入口处打印所有热源累计数组的初值如果入口初值非零那就说明清零漏了或者清零代码写在某个条件分支里没有执行到。第二坑是单位不一致。MFiX 内部大部分能量相关量采用国际单位制但处理颗粒物性时有的参数用摄氏度、有的用开尔文温差没问题但辐射这种涉及四次方的模型对温标极为敏感。实际代码里颗粒温度保存的基准温标在文件开头有定义读代码时不能想当然地认为整个文件从头到尾都一致。我吃过一次闷亏自己写的后处理脚本读取温度数组时默认是开尔文结果和源码里用的摄氏度混在一起画出的温度曲线看起来像随机噪声。第三坑是显式更新对时间步长的限制。颗粒温度更新用的是显式欧拉法稳定条件要求子步长必须小于某个由比热容、传热系数和颗粒尺寸决定的特征时间。如果用户为了加快计算把 DEM 子步长调得过大温度场会先出现震荡然后迅速发散。有些版本会在 calc_thermo_des 里面做温度截断限制颗粒温度不能低于某个值或高于某个值这能在一定程度掩盖稳定性问题但也可能掩盖物理上的错误。遇到温度场出现局部尖峰我建议先缩小步长验证再怀疑模型本身。第四坑是并行时源项归属处理。颗粒循环里每个颗粒自己累加的 Q_Source 不涉及跨进程通信但流体能量源项是按网格分布的一个网格在区域分解后可能被多个进程的颗粒贡献热量需要归约求和。这部分逻辑如果处理不对会出现周期性温度脉动而且脉动频率和颗粒跨进程迁移频率高度相关。排查时把进程数改成 1 跑同一算例如果温度场立刻恢复正常基本就能锁定是并行归约问题。除了这些坑calc_thermo_des 还要求颗粒温度初始值合理。如果初始温度给得过高第一个子步对流项就可能产生巨大热流如果初始温度给得过低接触导热项也会出现类似的问题。初始化和清零一样都是最不起眼但最容易出问题的地方。在数值上我还特别注意过能量守恒的全局检查。把 calc_thermo_des 算完后的颗粒总内能变化、流体能量源项总积分和壁面换热量加在一起应当等于零。如果误差超过一定阈值那说明这段代码里有变量引用顺序错了很可能是用了更新后的温度去算本该用更新前温度来计算的换热项产生了一个额外能量。这种做法并不总是错的有些格式在热源项里故意采用平均值以提高精度但你必须清楚它到底用的哪个温度否则检查守恒时就会一头雾水。6. 读懂 calc_thermo_des 之后的下一步改模型之前必须想清楚的三件事如果你已经能顺畅读懂 calc_thermo_des 的主流程和调用关系那你基本具备了修改或扩展这个文件的能力。但结合实际项目经验我想再说三件同样重要的事免得你花大力气改完代码之后发现方向本身有问题。第一件事是分清粒度和尺度。calc_thermo_des 是颗粒尺度的能量计算它假设每个颗粒内部温度均匀也就是所谓的 lumped capacitance 模型。这个假设对导热系数很高、颗粒尺寸很小的场景成立但对大颗粒或者低导热材料会带来显著误差。如果你要模拟木屑颗粒单靠 calc_thermo_des 的均匀温度模型可能精度不够得考虑引入颗粒内部导热或缩核模型那就不只是改这一个文件的事了。第二件事是注意能量方程的时间推进方式与流体耦合迭代的匹配。calc_thermo_des 在每个 DEM 子步里更新颗粒温度但流体能量方程一般是在流体步内更新的两个时间尺度相差很大。颗粒把热量传给流体流体温度变化又会反馈到下一子步的对流换热量上这种顺序耦合在温差大、热流大的工况下可能产生数值振荡。如果你发现温度场整体在呼吸颗粒和流体交替升温降温那也顺藤摸瓜查一查 calc_thermo_des 和流体能量方程之间的调用顺序是不是和你预期的一致。第三件事是验证比实现更重要。改完 calc_thermo_des 后先在简单算例上做解析解验证比如静止颗粒在均匀流体中的冷却过程温度应当按指数衰减。拿这种最少能手工算出解析解的算例去验证能快速暴露量纲、符号、时间推进错误。通过了这层验证再往复杂算例推进也不迟。我见过不少同行直接把新换热模型塞进去跑一个全尺寸流化床结果温度场完全不合物理最后发现只是某个单位转换小错误。先小后大这个顺序能省下大量排查时间。MFiX 的代码风格在开源 CFD 里算得上偏工程化变量命名直白但有时太长函数入口参数多到一眼看不完。读 calc_thermo_des 不需要什么高深技巧按模块归属—主流程—子函数—数值验证这个顺序来你会发现它并没有想象中那么难啃。整个过程最花时间的其实是把那些没有注释的隐式约定一点点试出来。希望这篇番外能帮你少走几步弯路少踩几个我踩过的坑。