1. 为什么会有这篇番外最近在啃MFiX的源代码主线的颗粒流体力学的部分还没啃完中途遇到一个不算起眼的小程序文件calc_thermo_des.f本来打算扫一眼就跳过去结果越看越觉得有意思。它不涉及复杂的算法推导也不涉及高深的数值格式但它是颗粒相能量计算里绕不开的一个环节。更重要的是它作为代码里的“附带题”把我之前零零散散看到的热传导、对流传热、辐射换热这些概念全都串了起来。先说点题外话。读源代码这件事很多刚入门的朋友容易陷入一个误区把主程序、核心模块啃得特别细遇到边边角角的子程序就直接跳。但实际模拟工程里真正让结果对不上的往往就是这些边边角角的细节。比如颗粒温度迭代发散、床层温度场分布不对、热流密度异常跳变排查到最后通常都在一些容易被忽略的辅助子程序里找到根因。我个人的习惯是主程序大框架看两三遍即可反而要把每个被调用的子程序彻底弄清“输入、输出、内部逻辑三步”因为代码的骨架是恒定不变的血肉才是各项目差异的根源。这篇文章就以calc_thermo_des.f为起点讲一讲它在MFiX DES离散元法模块里的定位、子程序的结构逻辑、关键参数怎么理解以及在调试和后续修改中常见的坑。为了便于阅读文章不会贴出完整源码而是基于该文件的整体逻辑按行理解的方式拆解。如果你正打算阅读MFiX的源代码或者正在使用DEM方法做颗粒流动复合传热方面的模拟这篇文章应该能帮你节约不少时间。2. 定位这个文件在MFiX代码体系中处于什么角色2.1 DES模块整体结构里的坐标MFiX里面和离散元相关的内容主要集中在des/目录下。calc_thermo_des.f在DES模块里的地位可以类比为“颗粒相能量计算的中转站”。为什么要强调“中转站”这个说法因为它是整个颗粒相热量计算中被调用的入口子程序但真正算东西的还在后面的具体函数里。我读代码时有个习惯拿到一个子程序先别急着看内部实现而是先查调用关系。calc_thermo_des.f的调用情况比较单一通常在 DES 的时间步进主循环中被显式调用。也就是说在每个DEM时间步长里程序都会调用一次这个子程序然后对每个颗粒进行热量计算。这个文件之所以命名为“calc_thermo_des”理解起来很简单calc 是计算的意思thermo 指热力学/热物理量des 标记它属于离散元模块。组合起来它就是“计算离散元颗粒相热物性相关量的子程序”。需要注意的一点这个文件在较新的MFiX版本中命名可能是calc_thermo_des_mod.f或类似带_mod后缀的模块化文件取决于版本。旧版本中以calc_thermo_des.f单独存在的形式比较常见后续版本为了统一代码结构把很多子程序都模块化了。读老版本代码和新版本代码之间最大的坑就是文件名的变化但内部逻辑基本是一脉相承的。2.2 它算的到底是什么既然叫 thermo那它计算的内容自然和温度、热量相关。具体来说这个子程序主要负责的是颗粒相的能量方程组装与求解准备。DEM方法中每个颗粒都有自己独立的温度这个温度受以下几种传热机制影响颗粒与流体之间的对流传热颗粒与颗粒直接接触时的导热颗粒与壁面之间的导热和对流辐射换热化学反应热如涉热反应颗粒MFiX将热量源项显式地分门别类在计算颗粒温度变化率DES_COND_T等变量之前先把相关的热量量纲数组都正确组装好。calc_thermo_des.f在这个组装过程中承担了“总体调度”的责任。它本身的逻辑并不复杂核心就三步对所有颗粒做初始化或清零操作判断哪种传热机制需要被激活调用对应的子程序将计算出的热流量分配到每个颗粒上理解了这个定位再看整个文件就会顺畅很多。它不是物理模型的实现者而是调度员。2.3 为什么还要单独一个“番外”很多主线的代码阅读文章会集中讲力学的部分接触力模型、阻尼、切向力、摩擦力以及颗粒在流场里的运动轨迹更新。热学部分在DEM里往往不是初学者的第一站但工程应用中几乎躲不开。说个我实际参与过的例子。某个处理颗粒干燥的模拟项目刚开始只看流动和混合不考虑颗粒温度。后来发现颗粒在设备里停留时间较长热风换热对颗粒含水率有决定性影响就必须加上热能计算。这个时候calc_thermo_des.f这类子程序就成了调试的主力对象因为颗粒温度算不对干燥速率就完全失真后面的组分输运全部白搭。所以这篇番外想换个视角不按照结构力学的路线走而是单独把热能计算这条线拉出来捋清楚。这样对于正在做干燥、燃烧、热解、气化这类主题的人会更有针对性。3. 参数逐个拆解进入子程序之前先认识这些变量3.1 子程序的入参与出参在MFiX的DES模块中calc_thermo_des.f的入参通常包含颗粒数、颗粒直径、位置坐标、速度、温度等基础信息以及流体相的温度、速度、体积分数等。出参则主要体现在温度变化率数组和热流数组上它们会在之后的时间积分中被用来更新颗粒温度。这里有个比较容易混淆的点MFiX中有大量的全局数组和局部数组。calc_thermo_des.f内部很多变量并不通过子程序参数列表传入传出直接在模块变量层面使用。这一点和纯函数式编程的思路差异很大但这也是老牌Fortran代码的常见风格。首次阅读MFiX代码的人最需要适应的地方就在这里一个子程序的“真实输入输出”要比参数列表多很多必须对照模块定义才能完整了解。为了方便理解我将该文件中涉及的核心数组做一个归类说明。需要注意的是不同版本的变量名可能略有差异但含义是稳定的。颗粒温度数组存储每个颗粒当前时刻的温度值。单位通常是K。热量源项数组存储每个颗粒在当前时间步内接收到的总热量单位通常是W。传热系数相关数组存储颗粒与流体之间的对流传热系数或者颗粒与颗粒接触时的等效导热系数单位W/(m²·K)或W/(m·K)。体积分数数组流体网格中的颗粒体积分数会影响传热面积和换热系数的计算。时间步长控制变量DES子循环中实际采用的时间步长在热流计算中用于后续显式更新。这些数组的关系可以这样理解颗粒温度是“状态量”传热系数是“传递能力”热量源项是“结果量”时间步长是“推进尺度”。四者构成完整的显式温度更新闭环。3.2 传热面积怎么来的颗粒与流体之间的对流传热工业上最常用的基础公式是牛顿冷却定律乘以传热面积。在DEM代码里每个颗粒的传热面积并不等于球体真实表面积因为不是所有表面都参与换热比如颗粒被其他颗粒紧挨着遮挡的部分就不参与。calc_thermo_des.f里有一种常见的传热面积计算方式如果颗粒自由悬浮在流体中默认传热面积等于球体外表面面积。如果颗粒与周围颗粒形成接触则会对传热面积做修正。实际代码中这个修正可能并不在calc_thermo_des.f内部完成可能由接触检测子程序算好之后传递给热学计算模块。但调度逻辑里必须存在可供判断的接触状态。从代码阅读角度我建议先把“传热面积从哪来”这个问题抓住。因为后续所有的热流量计算结果都直接正比于传热面积面积算错一处颗粒温度就整体偏移。在MFiX的DEM中单个颗粒默认传热面积通常取颗粒表面积也就是 π*d_p^2。有些版本中的颗粒直径是实际非球形等效直径有些则直接使用颗粒直径需要仔细看颗粒物性初始化部分的定义。3.3 努塞尔数和传热系数之间的桥你可能好奇对流传热系数本身是怎么确定的。这就涉及颗粒雷诺数和颗粒普朗特数再借助关联式计算努塞尔数。MFiX的DEM模块中的颗粒对流传热大多沿用经典的球形颗粒在流场中强迫对流传热关联式。例如常见的里兹-马歇尔型关联式Nu 2 0.6 * Re_p^0.5 * Pr^0.33这个式子的物理含义是当颗粒与流体的相对速度很小、颗粒雷诺数趋近于0时努塞尔数下限是2对应静止球体在静止流体中的分子导热极限第二项则代表强迫对流带来的额外增强。从代码角度看calc_thermo_des.f未必直接写这个关联式它可能调用一个独立函数来计算颗粒努塞尔数。我读该文件时比较关注的是以下几点颗粒雷诺数使用的是颗粒直径不是水力直径相对速度取的是流体速度与颗粒速度的矢量差绝对值气体在颗粒表面温度下的物性而不是来流平均温度下的物性在源文件中颗粒努塞尔数的计算一般发生在传热系数赋值之前。传热系数 h 则按照 Nu * λ / d_p 来还原λ 是流体热导率。这里有个容易被忽略的细节如果颗粒速度和流体速度完全相等相对速度为零那么 Nu 2传热系数退化为纯导热模式。这个退化极限很重要因为它在数值上提供了稳定性下限。4. 内部逻辑逐步拆解calc_thermo_des.f的代码流程4.1 第一步初始化或清零进入子程序后首先进行的操作通常是对所有颗粒的传热相关累加量做清零或初始化。这一步在很多代码里看起来枯燥但在数值稳定性上却至关重要。如果不做清零上一个时间步的残留热流量会被带到下一个时间步导致颗粒温度持续累积误差最终温度漂移、发散甚至产生非物理的负温度。在旧版本代码中这一部分可能是直接的循环赋值DO NP 1, MAX_PIP Q_SOURCE(NP) ZERO DES_COND_T(NP) ZERO END DO新版本代码中一般会引入数组整体初始化的写法但逻辑本质相同。我的建议是读源代码时先确认清零循环覆盖了哪些数组。如果后续新增了自定义传热数组却忘了添加清零操作那么颗粒温度会在几十万步时间推进后慢慢失控。这个坑我自己就踩过排查了几个小时最后发现只是初始化数组少了一行。4.2 第二步判断是否需要进行传热计算MFiX在设计上并不是每个网格、每个颗粒都必然执行全量传热计算它会通过逻辑变量和开关来控制哪些物理过程处于激活状态。这部分的判断常见有以下几类流体相计算是否开启传热颗粒相温度是否被纳入求解辐射模型是否被激活化学反应是否释放或吸收热量颗粒与壁面之间是否允许换热从调度逻辑角度看这些开关好比是控制闸门决定热量要不要流入某个通道。我之前遇到一个情况模拟结果中颗粒温度始终不变化排查到最后发现是颗粒能量方程开关没有开启连热量计算都没有进入。这类问题发生的原因通常不是因为源代码逻辑复杂而是因为输入文件中的DES_ENERGY_EQ或者ENERGY_EQ标志位设置不正确导致子程序被调用但内部提前返回。读源代码前先看输入开关往往事半功倍。4.3 第三步类型判断与分支选择MFiX中支持多种传热模型。calc_thermo_des.f内部最核心的分支结构是根据颗粒与颗粒的接触状态以及周围流场状态决定采用哪条传热计算路径。这个分支结构是我认为该子程序中最值得细读的部分。我先梳理几种典型路径。路径A自由悬浮颗粒只和流体换热不涉及颗粒间接触导热。这种情形下热量仅通过对流输运进入颗粒。路径B颗粒与颗粒相互接触存在接触导热。此时需要额外考虑两颗粒接触点附近的导热效应通常会涉及接触半径、等效热导率等参数。路径C颗粒与壁面接触涉及壁面导热。处理方式与颗粒间接触导热相似但壁面的热物性参数与颗粒不同。路径D同时存在上述多种传热机制需要叠加计算。在实际工程模拟中一个时间步内单个颗粒可能同时与多个颗粒接触并同时处于流体对流之中所以路径D才是最常见的。因此分支代码里往往并不是互斥的 if-elseif 结构而是并行的多个独立 if 判断每条判断都在给热量源项数组做累加。这里插一句经验把接触导热与对流传热耦合起来的时候最怕的是时间步长过大导致显式更新不稳定。热扩散的特征时间通常比颗粒碰撞的特征时间短如果时间步长按力学标准取温度更新很容易振荡。所以后续在时间积分环节常常需要对能量方程另做时间步长限制代码中一般会有单独的热时间步长控制。只不过这个控制不一定在calc_thermo_des.f中可能位于主DES时间步进循环里。4.4 第四步调用底层传热计算函数calc_thermo_des.f自身承担的“粗活”是调度和累加真正的“细活”在各底层函数或子程序中完成。常见的被调用对象包括颗粒对流传热系数计算函数颗粒间接触导热的计算子程序辐射换热的计算子程序壁面换热计算子程序这些底层的函数中有的可以直接独立测试。我的建议是阅读完calc_thermo_des.f后再单独去跟入其中一个底层函数比如颗粒间导热计算函数看看它在不同接触半径下的行为是否合理。曾有同行跟我交流过说他在读这段代码时觉得函数跳来跳去很麻烦线索容易断。我的习惯是手绘调用树理清楚calc_thermo_des.f一共调用了几层每层函数的输入是什么输出是什么。理完这个调用树整条颗粒相热量计算链路就装进脑子了。之后不管遇到多复杂的传热bug都能快速定位到具体某一层。5. 关键模型细节不能跳过的物理量含义5.1 颗粒间接触导热的等效处理接触导热在DEM代码里是一个比较微妙的问题。宏观直觉上两个颗粒相互接触热量会通过接触点传导。但真实颗粒表面粗糙接触并不是理想球面与球面的点接触而是包含若干微凸体的小面积接触。DEM为了工程可行通常使用等效接触半径和等效热导率的简化模型。MFiX在这部分的常见处理是根据两个颗粒的接触重叠量和杨氏模量等参数计算出接触半径再用颗粒材料的导热系数和接触半径构造接触热阻。从数值实现上看接触导热量往往正比于温度差乘以一个与接触几何相关的系数。这个系数形式类似于Q_cond H_c * (T_j - T_i)其中 H_c 为接触传热系数单位与热导乘以特征长度的形式有关。不同文献给出的模型表达式差异较大所以源代码中具体采用哪种形式只能逐个版本确认。在阅读这部分代码时我建议重点关注接触半径的计算方法。有些版本使用简化的几何关系有些版本则从法向接触力反推接触面积。两者在结果上的差异不容忽视尤其在颗粒堆积比较密实的床层中颗粒间接触点多接触导热路径复杂对温度分布的影响非常显著。5.2 壁面导热的特殊之处MFiX DEM中壁面通常被视为特殊的颗粒或半无限大固体。calc_thermo_des.f在该情形下的逻辑往往在颗粒循环中额外判断颗粒是否与某壁面单元重叠然后单独计算壁面导热量。壁面导热和颗粒间导热的差异主要体现在热物性参数赋值上。壁面的导热系数可能是常数也可能是温度的函数壁面温度可能是固定壁温也可能是给定热流边界条件。这些都会改变热流表达式的形式。我自己调试时有一个教训壁面温度如果在代码中直接使用初始温度而不是每步更新热量累积就会失真。这类问题不易发现因为初看温度场变化趋势是合理的但定量对照实验结果时差异明显。所以读源码时务必确认壁面温度的更新位置。5.3 辐射换热的简化颗粒相在高温场景下的辐射换热不能忽略比如燃烧、热解、气化。但DEM中颗粒数量巨大直接做逐颗粒的辐射角系数计算几乎不现实所以MFiX中往往采用各种简化辐射模型。常见的简化方式包括将颗粒视为灰体引入等效辐射温度场忽略颗粒间的辐射遮挡直接按经验系数修正把辐射热流作为体积源项在流体相中处理calc_thermo_des.f在辐射处理环节通常只是调用模型函数真正的实现放在对应的辐射子程序中。我从个人使用角度看如果模拟温度不超过700K左右辐射换热的贡献相对有限可以暂缓但如果涉及高温燃烧辐射权重非常大必须确认代码中辐射源项的量纲和符号否则可能把净热流方向搞反。量纲和符号问题是源代码阅读里最容易被忽视、却也最容易出大错的细节。热量传入颗粒为正传出为负这个约定看似简单但在多个子程序多层调用后很容易在某个环节乘一个负号导致颗粒温度不升反降。6. 实操过程中最常见的问题与排查建议6.1 颗粒温度整体偏高或偏低出现这种情况首先检查对流传热系数模型和颗粒面积的定义。一般来说整体温度偏高说明热量输入偏多要么传热面积偏大要么热流密度偏大。整体温度偏低则反过来。我实践中的排查顺序是这样的先检查颗粒初始化温度是否合理检查传热面积是否采用了表面积而不是投影面积检查努塞尔数关联式中的物性取值温度是否正确检查时间步长是否过大导致显式更新的数值耗散多数情况下问题出在物性参数的温度取值基准不统一。比如气体热导率如果取的是入口温度下的值而实际颗粒表面温度远高于入口温度那么两者差异会导致计算出的传热系数和真实偏差很大。6.2 颗粒温度出现不正常的局部尖峰局部尖峰通常和接触导热有关。当两个颗粒发生碰撞接触半径若被高估接触导热量会在一瞬间显著增大表现为局部温度尖峰。还有一个容易被忽略的原因是接触判断的逻辑利用了接触重叠量而重叠量在碰撞力模型中可能因为时间步长过大出现振荡进而导致接触导热热流也振荡。处理思路是降低DEM时间步长观察尖峰是否减弱如果减弱说明是接触力-导热耦合导致的数值问题。如果降低步长无效则检查接触半径计算中是否合理限幅。6.3 颗粒温度比流体温度还平稳如果颗粒温度基本不变而流体温度在变化检查颗粒与流体之间是否真的启用了对流传热。还有可能是因为颗粒热容设置得非常大比如误将颗粒体积热容写成了比热容导致热量输入对温度变化贡献极小。颗粒热容的量纲在MFiX中需要特别注意。有些模型输入的是体积热容J/m³·K有些输入的是质量热容J/kg·K弄混之后温度响应会差好几个数量级。我建议在读代码时顺手把这个量纲链条推导一遍热流量单位W时间步长单位s热量单位J温度变化量单位K温度变化量 热量 / (颗粒质量 × 颗粒比热容)。如果发现计算结果和这个关系对不上那一定是代码中定义的某个热容或温度相关变量理解有误。6.4 传热开关全开了但就是不算热这类问题多半是逻辑开关的命名理解有偏差。比如MFiX中某些开关是默认关闭需要在输入文件中显式开启而另一些开关默认是开启的。还有一部分开关之间存在从属关系比如总体能量方程开关不开即使颗粒能量方程开关开启也无效。排查这类问题我通常会在calc_thermo_des.f的入口处加一个临时输出打印当前步数、颗粒数、开关状态然后观察日志。读取Fortran代码时在关键位置临时插入WRITE(*,*)是非常常用的调试手法虽然原始项目代码中不推荐留这类输出但临时用一下确实高效。7. 读这个文件时值得顺手做的备忘7.1 版本差异备忘MFiX发展过程中热学计算模块经历了较大的结构重构。老版本里calc_thermo_des.f是独立子程序文件新版本里可能被整合进模块中作为subroutine保留。文件名后缀_mod不是简单改名而是从子程序转成了模块化封装调用方式也发生了变化。如果你正在读的版本中找不到calc_thermo_des.f不要意外搜索calc_thermo_des或者calc_thermo这样的关键词大概率能找到对应的模块。源码阅读工具的全局搜索在这里比手动翻目录高效得多。7.2 数组维度问题MFiX中的颗粒数组维度通常不是固定常数而是动态分配。这意味着循环边界条件往往和颗粒最大编号相关而不是和实际当前颗粒数相关。在calc_thermo_des.f中循环边界如果写的是MAX_PIP那就会遍历所有可能的颗粒槽位包括尚未激活的颗粒。未激活颗粒的物理量可能包含垃圾值因此代码内部几乎都通过IF (SOLIDS ...)或者颗粒存在性判断来保护防止对空颗粒做热量累加。阅读时如果发现循环体里面没有存在性判断就要警惕。这是少数会导致结果莫名包含NaN或无穷大的原因。7.3 并行与串行差异MFiX支持并行计算calc_thermo_des.f的逻辑在并行情况下执行时可能会涉及颗粒分布式存储即每个进程只负责自己区域内的颗粒。这样一来相邻进程的边界颗粒换热需要依赖通讯而不是简单的本地循环。如果你只是在串行环境下调试则不需要过度关注并行通讯逻辑。但如果要做大规模并行模拟就必须搞清楚颗粒热量更新时的边界交换策略否则颗粒温度在跨界时会出现跳变。关于并行版本的代码我现在还没有完全吃透后续如果继续深入读完MPI通讯相关模块再另外成文分享。8. 从代码阅读到实际修改的经验小结读源代码最大的收获不是记住了某个公式是怎么写的而是建立起一种“物理量在代码中如何流转”的直觉。算例出问题时能快速判断是物性参数问题、边界条件问题还是模型选择问题。calc_thermo_des.f这个文件虽然简单但它的调度逻辑几乎涵盖了DEM热能计算的全部入口。把它彻底弄懂之后再去看具体的传热模型细节就会发现不过是在这个框架里填入不同的计算表达式而已。下面是我在为一个流化床干燥算例修改代码时走过的一小段真实流程供参考。第一步确认模拟目标是颗粒含水率变化而不是单纯温度变化。因此需要在能量计算之外再引入组分输运和蒸发模型但蒸发模型的源项必须和calc_thermo_des.f产生的热量源项做叠加。第二步在calc_thermo_des.f中找到颗粒热量累加的循环在循环末尾单独加上一个蒸发潜热的热量扣除项。这个扣除项是负的因为蒸发会带走热量。第三步调试中发现颗粒温度在水分快蒸干时出现轻微振荡。原因是最小含湿率判断和蒸发源项切换之间的迟滞。解决办法是添加滞回区间让热源项平滑过渡。这段经历想说明的是读源代码不是目的达成自己的物理模型目标才是目的。但如果不读源代码连往哪里加蒸发源项都无从下手。9. 给新手的一些更实际的阅读方法9.1 先跑通一个算例再读代码硬啃源代码尤其是没有运行过案例的情况下很容易迷失在变量和循环里。我建议先找一个带传热过程的标准案例MFiX自带有许多测试案例跑通一次再带着“模拟结果是怎么算出来”的问题去读calc_thermo_des.f。这时候变量不再是纸上的符号而是和实际模拟曲线联系起来的量。比如跑一个颗粒静止床层的导热案例观察颗粒温度从初值趋向壁面温度的过程。再到代码里找哪个变量控制了这个导热行为会很直观。9.2 打日志比看逻辑快Fortran代码的阅读有一个“笨办法”在关键子程序入口和出口加上临时输出打印颗粒编号、颗粒温度、累计热流、接触状态等信息。虽然每次运行日志量很大但对于单个子程序的逻辑验证这种方式远比反复推敲代码逻辑快。我以前的调试习惯是加一个条件输出只打印特定颗粒编号的信息比如IF (NP 100) THEN WRITE(*,*) NP100 T, DES_TEMP_NEW(NP), Q, Q_SOURCE(NP) END IF这样就能锁定单个颗粒观察它在整个时间步进中的温度和热流变化。这个方法对排查颗粒温度异常特别有效。9.3 善用差异对比MFiX支持同一算例用不同模型开关跑出不同结果。直接对比这两种结果的温度场差异再回到代码里找对应条件分支比凭空猜测分支逻辑更直接。我之前对比过开启接触导热和关闭接触导热两种设置下的床层温度分布发现差异集中出现在堆积区说明接触导热模型在堆积区域起主导作用。据此再反推代码理解就快多了。10. 最后想说的源文件虽小串联的知识面却不小calc_thermo_des.f在MFiX整个项目里算是一个规模不大的子程序。但它串联起了颗粒相运动、流体与颗粒对流传热、颗粒间接触导热、壁面换热和辐射换热这些模型。如果只从代码规模上判断它不重要就很容易错过DEN模块热能计算体系的关键入口。我现在读到这个阶段越来越觉得读源代码像在拆一个复杂的机械设备看起来只有一台机器但每个零件都有自己的角色拆得越多对整个系统的理解越完整。calc_thermo_des.f就是其中一个不大但绝对不能说无关紧要的零件。如果你也在读MFiX源代码建议以后遇到这个文件时不要跳过花一两个小时把它的调用关系和变量流转弄清楚。之后再看任何和颗粒传热相关的代码都会有一种豁然开朗的感觉。