搞CFD的人十有八九都听过MFiX做气固两相流、流化床、反应器模拟的基本绕不开它。但很多初学者只把它当黑箱跑算例、调参数、看云图一旦遇到非预期结果或者想往里边塞自定义模型就完全抓瞎。我自己以前也是这样直到某次被一个颗粒温度发散的问题逼到墙角才下决心去啃源码。这个“番外”就是我读MFiX源码时关于calc_thermo_des.f这篇记录聊的是离散元框架里热物性是怎么算的、怎么被调用的以及读它时最值得注意的地方。如果你也在用离散颗粒方法做传热相关模拟或者正准备改MFiX源码这份笔记应该能帮你省不少时间。1. 先说清楚这个文件是干嘛的它管着颗粒的“温度账本”读源码最忌讳拿到文件就从上往下逐行硬啃先搞清楚它在工程里扮演什么角色效率会高很多。calc_thermo_des.f这个名字已经很直白了des对应MFiX里的Discrete Element Section也就是离散元模块calc_thermo就是计算热力学量。所以它在MFiX里的职责就是负责更新每一个宏观颗粒的温度、颗粒的定压比热容以及颗粒与周围流体、壁面之间的传热量。换句话说你用MFiX做流化床内颗粒加热、干燥、化学反应放热耦合这类模拟时最后落到颗粒层面上的那个温度变化基本都是在这里算的。1.1 在MFiX源码里怎么找到它以及它管到哪一层MFiX的源码目录结构按功能模块分开DES相关代码通常集中在model/des/这类目录下不同版本的路径可能有差异但文件命名一般比较稳定。calc_thermo_des.f就是其中专门负责DES热物性计算的单元。这里要先理解MFiX的算法框架。MFiX同时支持传统的欧拉双流体模型以及欧拉-拉格朗日理念下的离散元模型。在离散元模式下流体相仍然用网格、场变量来描述颗粒相则是一个个独立的拉格朗日粒子每个粒子有自己的位置、速度、角速度、温度、质量、直径等属性。颗粒之间的碰撞用接触力学解决颗粒与流体的作用力通过相间作用在流体网格上求解而颗粒的温度则通过能量平衡方程单独更新。calc_thermo_des.f做的工作正是能量平衡那一层。它在每个DEM时间步里被调用根据当前颗粒与周围环境之间的换热条件计算每个颗粒的温度增量然后更新到颗粒的状态数组里。所以它本身不负责求解流体温度场那是欧拉部分的energy_eq这类模块的事。这个“颗粒温度计算”和流体湍流里的“温度方程”不是一个东西理解错了后面会乱。1.2 从调用链看它的位置我自己阅读时会习惯先去翻调用关系搞清楚谁在调用当前文件。对于calc_thermo_des.f典型的调用路径是这样的MFiX主时间循环逐时间步推进在离散元时间步内部des_time_loop里依次处理颗粒运动、碰撞、作用力、热传递等颗粒运动处理完之后会调用calc_thermo_des更新颗粒温度等热物性更新后的温度再反馈给化学反应、颗粒物性、或者下一时间步的传热计算。也就是说calc_thermo_des.f在时间轴上的位置通常在颗粒位置速度更新之后而又在下一轮的流体颗粒相互作用之前。这个位置很关键温度更新早了换热来不及用上最新的颗粒位置更新晚了温度反馈又会晚一个时间步在强耦合工况下容易振荡。1.3 它和周围模块的关系读这个文件时还会看到它频繁调用或使用其他模块的数据。比如颗粒的基本物理属性数组包括颗粒温度DES_T_s、颗粒质量、直径流体相在颗粒所在位置插值出来的温度、组分浓度、速度传热系数计算要用到的无量纲数比如雷诺数、普朗特数、努塞尔数化学反应模块可能需要的颗粒温度值。所以读这个文件前至少要对MFiX里常用的几个全局数组名有印象不然满屏幕DES_T_s、Q_source、DES_C_ps看下来很容易迷失。我的经验是先花半小时把颗粒属性数组的命名规律捋一遍再回头读这个文件几乎事半功倍。2. 打开文件后先抓主线温度更新、传热系数与热物性真正打开calc_thermo_des.f之后代码不会太长算上注释一般几百行到上千行左右但不要把精力浪费在格式和注释的细枝末节上。我的阅读顺序是先抓三样东西一是颗粒温度到底怎么更新二是传热系数从哪来三是热物性、初始条件这些常量是怎么处理的。2.1 颗粒温度更新的基本方程从物理上说单个颗粒的能量守恒方程可以简化为[ m_p c_{p,p} \frac{dT_p}{dt} Q_{conv} Q_{rad} Q_{cond} Q_{react} ]其中m_p是颗粒质量c_{p,p}是颗粒比热容T_p是颗粒温度等号右侧分别是颗粒与流体的对流换热、颗粒与壁面或其他颗粒的辐射与导热、以及化学反应放热或吸热项。在calc_thermo_des.f里温度更新的代码逻辑就是在组装这个方程的各部分。常见代码体现为先对每个颗粒初始化累积换热量的数组然后依次加入对流换热项、辐射项、以及颗粒间接触导热的贡献最后用总换热量除以颗粒质量与比热容的乘积得到温度变化量再乘上时间步长更新到新的DES_T_s里。我读的版本中温度更新主体大概是这种思路DO NP 1, MAX_PIP IF(IS_DEM_ACTIVE(NP)) THEN Q_INTERNAL Q_SOURCE(NP) Q_TOTAL Q_CONV(NP) Q_RAD(NP) Q_INTERNAL IF(DES_C_PS(NP) ZERO) THEN DES_T_S_NEW(NP) DES_T_S(NP) Q_TOTAL * DTSOLID / (DES_M_PP(NP)*DES_C_PS(NP)) ENDIF ENDIF ENDDO这段是简化示意不同版本变量名会变但整个结构基本就这样。注意这里用的是DTSOLID也就是DEM子时间步长不是外层流体的时间步长。颗粒温度更新必须是“小步长、多次更新”的节奏因为颗粒传热的特征时间往往比流体流动小得多。如果你在改代码时直接用外层时间步长去更新颗粒温度大概率温度发散。2.2 传热系数拆开看别被公式吓到对流换热项(Q_{conv})的通用形式是[ Q_{conv} h_p A_p (T_g - T_p) ]其中(h_p)是颗粒表面对流传热系数(A_p)是颗粒表面积(T_g)是颗粒所在位置的流体温度。calc_thermo_des.f本身往往不会从头计算所有传热学公式通常会调用一个专门的传热关联式函数或者直接在文件内通过雷诺数和普朗特数计算努塞尔数。最常见的球形颗粒对流换热关联式是Ranz-Marshall形式[ Nu 2.0 0.6 Re^{1/2} Pr^{1/3} ]然后通过h Nu * k_g / d_p得到传热系数其中k_g是流体的导热系数d_p是颗粒直径。这个公式看着简单但坑不少。第一普朗特数在不同温度下变化很大必须用当前流体膜温度下的物性不能拿常温物性硬套第二雷诺数里的相对速度涉及颗粒速度、流体速度插值还要考虑颗粒所在网格位置第三如果颗粒不是球形还要额外修正表面积和体积。这些在源码里都会体现为各种判断条件和系数修正读的时候要慢下来逐个确认。2.3 热物性参数与初始条件处理calc_thermo_des.f里另一块大头是热物性与初始条件处理。每个颗粒的温度初值从哪里来DES_T_S数组在初始化时是否都被赋值比热容能不能随组分变化而变化这些问题在文件中都会有对应逻辑。我读到的版本里文件开头通常有一段对颗粒比热容和温度的初始化或重置逻辑。比如当颗粒从某个反应器中进入系统时需要设置入口颗粒温度当颗粒比热容由组分决定时需要在温度更新之前先根据组分更新比热容。特别注意MFiX里颗粒组分相关的数组是DES_X_s这类如果化学反应模块改了颗粒的组分而这里的比热容没有同步更新温度计算就会失真。另外calc_thermo_des.f里经常出现“是否考虑颗粒内部导热”的开关。对于大颗粒内部温差可能不可忽略需要引入一维热传导模型对于小颗粒通常默认颗粒温度均匀也就是集总参数法。这个开关通常由一个运行时输入参数控制读代码时一定要搞明白默认值是什么否则你明明想模拟大颗粒结果代码走的是集总参数分支算出来的温度自然对不上实验。3. 实操路径我是怎么把calc_thermo_des.f读顺的读这种古老又庞大的Fortran代码库忌讳硬啃一定要有实操路径。下面是我这次读calc_thermo_des.f时实际采用的几个步骤基本可以复用到MFiX其它源代码文件上。3.1 第一步先建调用关系表再碰代码我在本地把与calc_thermo_des有关的调用点、被调用的函数名、共享的全局数组列成了一张表手动整理成类似下面这种关系模块/函数与 calc_thermo_des 的关系主要共享数组des_time_loop调用 calc_thermo_desDES_T_S、DTSOLID流体相能量方程提供流体温度、换热系数T_G、K_G颗粒碰撞模块提供接触导热判断接触面积、重叠量化学反应模块读取更新后的颗粒温度DES_T_S、Q_SOURCE不用画多复杂的图就一张表格足够你把数据流理清楚。做完这一步你就知道这个文件在时间步里的位置以及它使用哪些外部数据。后面看代码时遇到一个数组名马上能想起它的来源不会读着读着脑袋打结。3.2 第二步逐行读时按“状态更新类、系数计算类、数据交换类”分类Fortran源码一长语句就杂。我读calc_thermo_des.f时习惯把所有语句按功能分成三类状态更新类直接改变颗粒温度、比热容等状态变量这类是核心优先看系数计算类计算传热系数、辐射系数、努塞尔数等这类要对照公式看单位和公式版本数据交换类从流场数组取值、向颗粒数组写值、或者处理与MPI通信相关的数据收发这类要特别注意数组边界和并行条件下的随机性。例如这个简化片段就属于典型的“系数计算类”RE_G RO_G * ABS(V_G - DES_V_S(NP)) * DP(NP) / MU_G PR_G MU_G * C_PG / K_G NU_G 2.0D0 0.6D0 * SQRT(RE_G) * PR_G**(1.0D0/3.0D0) H_CONV NU_G * K_G / DP(NP)这类语句里最容易出错的就是单位。MFiX中默认单位制是SI但有些老版本在使用过程中会因为历史遗留原因在输出时换算成摄氏度在内部计算时却用开尔文。如果你只盯着公式看不追踪单位很容易在后处理时发现温度整体偏移几百开尔文。3.3 第三步改完代码后用一个最小算例验证我这个人的习惯是不亲手改一遍代码等于没读懂。读calc_thermo_des.f的过程中我搭了一个最简单的单颗粒算例一个颗粒放在静止流场里给定初始温度打开对流换热让颗粒逐步与流体达到平衡。颗粒温度曲线应该符合经典牛顿冷却公式也就是指数衰减趋于流体温度。改代码后跑这个算例观察颗粒温度曲线是否平滑收敛速度是否合理。这能用来验证很多问题时间步长是否足够小、传热系数计算公式是否正确、更新逻辑是否被执行到。如果用任何复杂的流化床算例去验证变量太多出了奇怪现象根本分不清是温度更新的问题还是流体场的问题。4. 读这个文件时我踩过的几个坑这部分是重点。单纯的代码阅读很枯燥但当你真踩过坑再回头看代码印象就深了。下面几个坑都是我或者周围同事实际遇到过的问题。4.1 数组下标与颗粒活动状态的判断MFiX里的颗粒并不是所有都始终活跃有些颗粒可能已经离开计算域有些被删除有些在边界外等待进入。calc_thermo_des.f的循环里往往会加IF(IS_DEM_ACTIVE(NP))之类的判断。如果漏掉这种判断或者判断条件写错了容易出现“幽灵颗粒温度”问题——一个已经离开系统的颗粒还在更新温度或者一个尚未进入系统的颗粒参与换热计算导致温度数组出现NaN进而影响整个流场。有一次我调代码发现计算域内某个区域温度出现尖锐峰值找了半天没找到原因最后发现是删除颗粒的时候没有清理DES_T_S数组那个颗粒位置残留了一个极端温度值在下一次传热计算时又被“复活”了。这种事在Fortran静态数组里特别容易发生因为数组不会自动清理。4.2 单位制的连环坑MFiX内部基本统一使用SI单位制但你在看这个文件时要特别小心温度的单位。SI下温度是开尔文但很多后处理脚本、化学反应数据或者物性表里喜欢用摄氏度甚至有些算例文件的初始条件里直接填摄氏温度。如果代码里某处把摄氏温度赋值给了内部温度数组那整个能量守恒就全乱了。我判断单位是否一致时习惯看一个关键点在温度更新代码里是否有类似T_STORE T_REF T_OFFSET的处理。如果有多半是单位转换或参考点修正。没有的话默认全代码走开尔文你在输入文件里也要用开尔文否则算错不怪人家。4.3 编译选项与宏开关的影响MFiX源码里有大量预处理指令也就是#ifdef、#if这类。calc_thermo_des.f虽然不大但同样会受编译选项影响。比如是否编译化学反应模块、是否启用辐射传热、是否启用颗粒内部导热模型这些都会显著改变代码走向。我第一次读这个文件时直接在编辑器里搜索“辐射”关键词居然没搜到多少内容后来才发现辐射模块的代码被宏开关包裹着默认编译配置下根本没编译进去。开启对应的宏开关后量突然多出一大截。所以读代码前先看一下你当前要分析的这个算例用的是哪个编译配置确认哪些宏被定义否则你看到的“完整源码”其实是残片。4.4 MPI并行下颗粒温度的随机性如果你跑的是多进程并行算例还要注意calc_thermo_des.f里颗粒数据的分布方式。MFiX在并行下会按空间域分配颗粒每个进程只负责自己管辖区域内的颗粒计算。可是颗粒会移动可能从一个进程边界穿到另一个进程的域里这时如果热物性数组没有做同步就会出现同一颗粒在两个进程里的温度不一样模拟结果不稳定。排查这种问题很难因为每次运行的结果都略有差异。我的经验是先把并行进程数设为1跑一遍如果结果稳定再去跑并行对比。如果单进程稳、并行不稳多半就是并行数据交换或数组同步的问题重点检查calc_thermo_des.f对外部颗粒数据通讯部分的代码。5. 常见问题排查速查表读源码的过程里很多问题会反复出现。我把它们整理成了一个速查表方便你事后对照排查。问题现象可能原因排查方向颗粒温度直接变成NaNDES_C_PS为0或负值除以0检查初始比热容设置确认颗粒质量数组有效颗粒温度长期不更新颗粒不满足活动状态判断检查IS_DEM_ACTIVE(NP)逻辑和颗粒进出流场逻辑温度曲线振荡发散DEM时间步长过大缩小DTSOLID必要时做子循环稳定性判断单颗粒温度与解析解不符传热系数公式、单位、物性参数对比努塞尔数关联式的参数版本核对流体导热系数并行结果与单进程不一致颗粒数据跨进程同步问题检查涉及颗粒数组的MPI通信或临时用单进程复现改完代码后报错未定义变量编译宏未打开检查生成Makefile时的编译配置和宏开关颗粒温度异常升高但流场温度正常化学反应放热项叠加多次检查Q_SOURCE是否在多个模块里被重复累加这张表不是万能的但它覆盖了我在实际项目中遇到的一大半问题。你如果也卡在类似现象上可以先按表格里的排查方向走基本能少走弯路。6. 从calc_thermo_des.f延伸出去的几条路读这个文件后其实可以顺势延伸出几个方向对后续建模会很有帮助。6.1 与化学反应模块的耦合注意点如果模拟的是煤燃烧、生物质热解这类过程颗粒温度会直接影响反应速率是否符合阿伦尼乌斯关系。在calc_thermo_des.f更新温度之前化学反应模块通常会先计算当前温度下的反应速率然后生成热量源项Q_SOURCE再送到温度更新里。这个顺序很关键温度先更新还是反应热先计算结果差别很大。我一般会把化学反应产热项看作“温度显式源项”也就是用当前时间步的旧温度去近似这一个时间步内的反应放热温度更新后再进入下一个反应步。时间步足够小的情况下这种显式处理是稳定的。一旦你发现反应放热导致温度步进振荡第一反应不是改温度更新格式而是去检查反应步长是否超过特征时间的约束。6.2 辐射传热与大颗粒内部导热辐射传热在高温和颗粒稀疏区域尤其明显。calc_thermo_des.f里处理辐射项的方式常常会假设颗粒表面灰体、按斯蒂芬-玻尔兹曼定律计算。颗粒内部导热则需要把颗粒离散成多层同心壳层每一层有自己的温度。别看文件不长涉及这些模型时代码复杂度会成倍上升。我的建议是先默认关闭这些高级模型跑通基础对流换热后再逐步打开这样便于定位每个模型带来的变化。6.3 热物性与组分变化的联动颗粒比热容不是固定不变的尤其当颗粒发生反应后组分变了DES_C_PS也要跟着变。calc_thermo_des.f里如果只在初始化时设置了一次比热容后面都不再更新那反应过程中颗粒温度就会偏真值越来越远。读代码时留意这个数组是否有在循环体内被赋予新值如果没有那么在耦合化学反应的地方就得手动补充更新逻辑。这一点项目里真正踩过坑的人都会认同比热容不更新的模型只适合经典的“无反应传热”验证算例跑反应模拟很容易翻车。写在阅读之后的一点心得我最初读calc_thermo_des.f的时候满脑子想的都是“怎么快点看完”现在回头看真正有价值的是阅读过程中查公式、翻调用关系、做最小算例的那几个小时。尤其是一开始我总把注意力放在颗粒温度更新的那一两行代码上后来才发现真正决定模拟成败的往往是被忽略的物性更新逻辑和单位换算。给还在啃MFiX源码的朋友一个建议不要追求一次读懂所有分支先沿着一条最简单的路径走通——单个颗粒的对流换热、温度变化、后处理输出曲线把这条链路搞清楚再去碰化学反应和辐射耦合会顺很多。至少我自己现在拿到新版本的MFiX看到可疑的温度现象时第一反应依然是回到calc_thermo_des.f这个原点把颗粒的能量收支持平了再说。
