文章目录

1 目标
2 代码结构
- 2.1 des_time_step外层总控制
- 2.2 颗粒导热
- - 2.2.1 颗粒-颗粒导热：calc_force_dem.f和des_thermo_cond_mod.f
  - 2.2.2 颗粒-墙壁导热：calc_collision_wall_mod.f和des_thermo_cond_mod.f
- 2.3 除导热外其他三种传热的控制：calc_thermo_des
- - 2.3.1 CONV_GS_DES1 对流传热量
  - 2.3.2 DES_RADIATION计算辐射量
  - 2.3.3 RXNS_GS_DES1计算反应传热
- 2.4 总热源项置零：des_thermo_newvalues
- 2.5 显式耦合计算传热方式（待完善）
3 修改代码以输出传热量
- 3.1 总传热量
- 3.2 颗粒颗粒导热量
- 3.3 颗粒墙壁导热量
- 3.4 对流传热量和反应传热量（显式）
- 3.5 辐射传热
4 算例验证
- 4.1 算例要求
- 4.2 结果展示（待更正）
- - 4.2.1 总源项 des_usr_var 1
  - 4.2.2 颗粒颗粒导热 des_usr_var 2
  - 4.2.3 颗粒墙壁导热 des_usr_var 3
  - 4.2.4 对流 des_usr_var 4
  - 4.2.5 反应热 des_usr_var 5
  - 4.2.6 辐射des_usr_var 6
- 4.3 验证结果分析
附录：子程序调用关系图（显式待完善）
- 2020-10-28：出现问题，总传热量不等于五项和
- 2020-10-29：问题有可能出在输出时刻不同上
补充3：输出传热时刻S_TIME
- 操作方法
- 结果
源代码下载
- 第1版
- 第2版

1 目标

分别输出mfix-dem中的颗粒的1.总传热量、2.导热传热量、3.对流传热量、4.辐射传热量和5.化学反应传热量

一共五项，存储到des_usr_var的第1到第6个数组中（后面会看到，导热分为两种：颗粒颗粒以及颗粒墙壁）。

关于des_usr_var的用法见前面的博文。

前面博文的地址：

https://blog.csdn.net/weixin_43940314/article/details/108138825

2 代码结构

首先搞清楚计算传热的大致代码结构。

计算传热方面有两条路径：隐式耦合和显式耦合。两者耦合方式我现在还没有完全弄明白，但是大概知道几点：

这里的耦合指的是气体和颗粒之间的耦合。
隐式耦合在每个固体时间步都进行传热量的耦合。但是显式耦合只在流体时间步和颗粒时间步的开始进行传热耦合。
尽管user guide上面说有传热传质和反应的模拟无法使用显式耦合，但是这个说法显然过时了。在目前的版本里是可以用的。

显然隐式的计算量要大很多，但是计算会更加精准，而且按理来说应该更容易收敛（因为传热传质每次加入的源项小了）。

调用关系图我放在附录了，可以先去瞟一眼有个大致印象。（不一定准确，是我自己debug一步步画出来的）

隐式和显式耦合计算四种传热方式的代码都是一致的。但是调用关系和路径不一样。

计算四种传热量隐式耦合的路径更加清晰简单，我们就以隐式为例子。（虽然我实际计算用的显式耦合，但是我目前显式耦合的路径还没研究透彻，以后待完善）

以下要用到这几个源文件：

des_time_march.f(整体控制)
calc_thermo_des.f（整体控制）
des_thermo_newvalues.f（总热源项置0）
des_thermo_conv.f（对流相关）
des_thermo_rad.f（辐射相关）
rxns_gs_des1.f（反应热相关）
des_thermo_cond_mod.f（导热相关）
calc_force_dem.f（导热相关，颗粒之间导热）
calc_collision_wall_mod.f（导热相关，颗粒-墙壁导热）

为了结构清楚，我们按照调用路径来讲解。

（下面我们主要指的都是函数或者子程序）

最外层:run_dem中调用des_time_step。这里就不多说了。

2.1 des_time_step外层总控制

次外层：des_time_step

位置des_time_march.f

关注des_time_march中的174行（颗粒颗粒导热）、178行（颗粒墙壁导热）、183行calc_thermo_des（除了导热外其他传热方式的总控制）和193行des_thermo_newvalues（把总热源项置零）

2.2 颗粒导热

导热是比较特殊的，因为它只有当颗粒相接触或者很近的时候才发生。因此把它写在了碰撞里。

分为两种：颗粒颗粒导热和颗粒墙壁导热

2.2.1 颗粒-颗粒导热：calc_force_dem.f和des_thermo_cond_mod.f

calc_force_dem.f第165-175行

QQ_TMP是存储导热量的临时变量

des_conduction是定义在 module des_thermo_cond的一个函数

进入des_thermo_cond函数，位置在des_thermo_cond_mod.f第41行

Fortran中的函数，传出的变量默认是与函数同名的变量，因此直接看DES_CONDUCTION这个变量的赋值语句部分，在199行

就是颗粒颗粒直接传热加上颗粒-流体-颗粒传热。具体传热机制请看Musser的博士论文。

根据函数头给定的信息，把形参和实参对照，可以发现实参变量LL就是颗粒I的编号，I就是颗粒J的编号。sqrt(DIST_MAG)是颗粒中心距离，PIJK(LL,5)是颗粒的相，PIJK(LL，4)是颗粒所在流体网格编号。

因此LL应该就是当前的颗粒，而I是LL的相邻的一个颗粒。假定方向为：热量从I传递给LL(可以为负号)。

2020-09-17补充
偶然在des_granular_temperature.f中发现了PIJK数组的用法

第一个下标指的是颗粒编号。第二个下标从1到5分别代表：网格I编号J编号K编号，IJK编号，最后是相编号。

注：

165-175这一段是在两个嵌套的循环内的

do_index就是颗粒的编号NP，只不过是当地变量。start_index和end_index分别是1和最后一个颗粒编号
CC就是邻域颗粒的编号。内层循环就是在邻域颗粒内循环。

2021-5-10补充
目前碰撞导热的时间问题还没有弄清楚。
其他几种传热形式如热对流等均为连续的传热，因此没有时间的问题。只要统一用W做单位，在能量方程里取单位时间的传热量就行了
然而碰撞导热是离散的，只在碰撞期间传热。因此就出现时间的问题。
目前可以确定的是，从des_thermo_cond_mod.f中计算的DES_CONDUCTION函数的值，得到的是单位时间的传热量
如图所示

可见注释中明确标注了是Rate。

2.2.2 颗粒-墙壁导热：calc_collision_wall_mod.f和des_thermo_cond_mod.f

calc_collision_wall_mod.f中

934-936行调用DES_conduction_wall函数计算导热量。这个函数同样在des_thermo_cond_mod.f里

939行把墙壁颗粒导热加入到Q_Source

2.3 除导热外其他三种传热的控制：calc_thermo_des

在calc_thermo_des.f文件的第50-64行

这段翻译一下就是

首先看是不是显式耦合的，如果是显式耦合，直接把传热量加进热源项Q_source之中

如果不是显式耦合，分别调用三个子程序计算对流传热、辐射传热和反应传热

于是我们分别看这三个调用的子程序。

2.3.1 CONV_GS_DES1 对流传热量

位置：des_thermo_conv.f

模块名：RXNS_GS_DES1_MOD

看69-88行

和之前的逻辑一样，也是个IF结构

如果是显式耦合，直接把对流传热量给到CONV_QS这个变量中

如果是隐式耦合，给到Qcv这个变量（这是个当地变量，出了子程序就会消失），然后再加和到热源项中

2.3.2 DES_RADIATION计算辐射量

位置：des_thermo_rad.f

模块名：des_radiation_mod

看67-94行

有点儿长，前面是具体的怎么计算的，着重看最后一个IF ELSE语句，也就是87-91行

MPPIC的话就多乘以一个颗粒统计重量（就是一个颗粒团中有几个实际颗粒）

不是PIC的话就直接加上辐射传热量Q_rad

2.3.3 RXNS_GS_DES1计算反应传热

在RXNS_GS_DES1.f的第188行，调用子程序calc_rrates_des

跳转到calc_rrates_des.f以后，找到第222行-227行，这一行是存储反应传热量的（这几行往上是计算反应传热的，这里为了简略没给出），变量为lQSource。

如果是显式耦合，把反应热赋值给RXNS_QS这个数组

如果是隐式耦合，把反应热直接加进Q_source

这里也就解释了为什么之前calc_thermo_des只有显示耦合把RXNS_QS加进去了。因为隐式耦合在这里就已经把反应热加入源项了。

2.4 总热源项置零：des_thermo_newvalues

这个文件在之前的关于des_usr_var的博客https://blog.csdn.net/weixin_43940314/article/details/108138825 中已经说过了，就是更新颗粒的温度、热源项等等热力学相关量的。和cfnewvalues类似。

这里再复述一下关键行

看第67行，这里是把热源项置为0，意思是开始新一轮计算前先情况变量。所以可以认为到此为止所有的传热量都被加进到了热源项了。这一行之前的Q_source就是总的热源项。

2.5 显式耦合计算传热方式（待完善）

显式耦合的路径只有每个流体时间步才计算一次热源项。

目前了解的有这么几种：

导热：和隐式完全一样，在des_time_step中（也就是每个固体时间步）调用碰撞的时候计算导热。

对流：在DES_TIME_INIT，也就是固体时间步循环的开始，流体时间步刚结束之后。在des_time_march.f的115行

反应热：在流体时间步（run_fluid），计算反应的同时调用了计算反应热的子程序。具体关系看调用关系图。中间隔了好几层有些复杂。

3 修改代码以输出传热量

注：修改源代码前先把源代码复制到自己的文件夹下。

因为导热有两种，所以一个六个变量

编号	传热方式	原本存储的变量	我们存储的变量
1	总传热量	Q_Source	des_usr_var(1,:)
2	颗粒颗粒导热	DES_CONDUCTION/QQ_TMP(本地)	des_usr_var(2,:)
3	颗粒墙壁导热	DES_Qw_cond（共享）/DES_CONDUCTION_WALL(共享)/QSWALL(本地)	des_usr_var(3,:)
4	对流传热	CONV_Qs	des_usr_var(4,:)
5	反应热	RXNS_Qs	des_usr_var(5,:)
6	辐射传热	Qrad	des_usr_var(6,:)

注：斜线后面标注的是本地变量，前面是共享变量（封装在module里）。加粗的是我实际使用的变量。

我们既可以把本地变量存储起来，也可以存储共享变量。不过共享变量有可能在后面某些地方会被改变，而我们不知道在哪，所以无法确定是否就真的是这个值。本地变量还有一个好处，就是不用关心它是显式还是隐式的路径，只要在实际计算的地方添加就行了。

注：添加的位置请看图中行号和上下文。

注：有些代码前面直接整个USE discretelement了，所以不用再单独添加 USE discretelement, only: des_usr_var。如果没有的话要加上。

3.1 总传热量

之前的那篇博客已经说过，这里再重复一下

在des_thermo_newvalues.f implicit none语句前添加

  !######CL_200914#######USE discretelement, only: des_usr_var!######################

如图

同一个源文件，在Q_source(:)=0前添加

         !#########CL_200912##############des_usr_var(1,:)=Q_source(:)!################################

如图
68行（大概）后

3.2 颗粒颗粒导热量

实际上还可以细分成颗粒颗粒直接导热和颗粒-流体-颗粒的间接导热，这里就合在一起了。

在calc_force_dem.f添加

              !#######CL_200912#########des_usr_var(2,LL)=QQ_TMP;!#######################

位置如图

170行后

3.3 颗粒墙壁导热量

在calc_collision_wall_mod.f添加

               !########CL_200914############des_usr_var(3,LL)=QSWALL!#############################

如图
937行后

3.4 对流传热量和反应传热量（显式）

这里有些许的不同，显式和隐式必须分开。因为即使是本地变量的存储方式也不一样。隐式的直接就在计算的位置加入到热源项里面了（所以他才只需要call子程序就行了）。

我们的算例是显式的。先添加显式。隐式有机会再测试。

在calc_thermo_des.f

  !###########CL_200912###########des_usr_var(4,:)=CONV_Qsdes_usr_var(5,:)=RXNS_Qs!###############################

65行后

3.5 辐射传热

辐射略有不同，显式的我没搞清楚他在哪里调用的des_conduction（或许显式压根就没算辐射？以后再研究），隐式的在前面那张图里调用了。但是不管显式还是隐式，直接存本地变量应该都没问题。

在des_thermo_rad.f

           !#########CL_200912#########des_usr_var(6,NP)=Qrad!###########################

91行后

2020-9-15补充
没有输出辐射传热，可能是mfix bug显式没有调用计算传热的代码？

由于没有发现显式耦合中调用了计算辐射的子程序，所以和朋友讨论以后决定自己调用一下

在calc_thermo_des.f中加入

 !#########CL_200915##############IF(any(CALC_RADT_DES)) CALL DES_RADIATION!#######################

如图

代码是模仿下面那一行隐式调用的。至于为什么不模仿显式调用呢？因为des_thermo_rad里面不管显式还是隐式都把辐射加入到源项了（如下图），所以只要调用，就自动加入源项。

之前des_usr_var(6,NP)那一行还是不变。这样把本地变量Qrad存储进来。

4 算例验证

4.1 算例要求

打开传热传质和化学反应
打开des_usr_var并给出数组上限（在这里是6）
如果考察墙壁导热，应该设置墙壁BC（这里给确定温度）。

一般来说我们不单独设置墙壁BC的话，默认会设置为绝热的无滑移边界，这样就不会产生任何热量交换也就看不见墙壁导热了。

4.2 结果展示（待更正）

4.2.1 总源项 des_usr_var 1

4.2.2 颗粒颗粒导热 des_usr_var 2

4.2.3 颗粒墙壁导热 des_usr_var 3

由于一开始设置算例没有设置墙壁BC，导致默认是绝热壁，因此没有值。
2020-9-15已更新，可以看到贴壁的颗粒是有传热量的。

4.2.4 对流 des_usr_var 4

4.2.5 反应热 des_usr_var 5

4.2.6 辐射des_usr_var 6

2020-9-17
这个就是自己手动在calc_thermo_des.f里面的显式耦合路径下调用辐射子程序的结果

4.3 验证结果分析

在某一瞬时时刻，对某一生物质颗粒，加和各个输出传热。

做平均传热量（所有生物质颗粒平均）随不同时刻变化线图
先mark一下，以后有空再做。

附录：子程序调用关系图（显式待完善）

2020-10-28：出现问题，总传热量不等于五项和

测试了一下，发现了一些问题。主要就是总传热量并不等于后面五项之和

下面这个表格表示的是所有颗粒的各项传热，对全体颗粒平均后，随时间的输出。

sum2to6表示后面五项之和，最后一列表示五项之和减去总传热量
标黄的表示差值的绝对值大于1E-6的，标红的表示大于1E-4（总传热的量级在1E-4左右）

标黄的占了1951个，标红的占了15个

可见总传热量绝对不等于后面几种传热量之和！

这里因为没把墙壁设为定温BC所以没有墙壁颗粒导热。
问题出在哪呢？目前还没找到解决方案

2020-10-29：问题有可能出在输出时刻不同上

看这篇的最后结果

https://blog.csdn.net/weixin_43940314/article/details/109353312

时刻还是存储到des_usr_var当中
在这里插入图片描述从左到右依次为：
颗粒编号，
总传热计算时间，
导热计算实际，
与墙壁导热计算时间（因为没开墙壁导热所以没计算），
对流换热计算时间，
反应热计算时间，
辐射热计算时间
总传热计算时间减去导热计算时间

发现其他几项都还好，导热相差时间很大！
计算总传热量和计算导热的时间差了几十秒

补充3：输出传热时刻S_TIME

操作方法

编号	传热方式	原本存储的变量	传热存储的变量	传热时间S_TIME存储序号
1	总传热量	Q_Source	des_usr_var(1,:)	13
2	颗粒颗粒导热	QQ_TMP(本地)	des_usr_var(2,:)	14
3	颗粒墙壁导热	QSWALL(本地)	des_usr_var(3,:)	15
4	对流传热	CONV_Qs	des_usr_var(4,:)	16
5	反应热	RXNS_Qs	des_usr_var(5,:)	17
6	辐射传热	Qrad	des_usr_var(6,:)	18

增加传热时间输出量S_TIME
仍然采用des_usr_var的方式

在每个源代码前面增加导入模块
增加代码如图
1 Q_source处

2 QQ_TMP处

3 QSWALL处

4和5 CONV_QS和 RXNS_QS处

6 Qrad处

结果

这是在t=5.00s的时候的生物质颗粒的情况

可以看见，7001所有时刻均一致
但是7002开始，均有变化。而且16\17\18总是和13保持一致，但是14和15常常和13差距很大

也就是说，两种导热–颗粒的碰撞导热（QQ_TMP）和颗粒与墙壁的接触导热(QSWALL)，与其他时间不能保持一致。

这也是应该的，因为碰撞发生的时间是随机的，间断的。而对流、辐射、反应是无时无刻不在发生的。

所以这也是为什么碰撞导热的时刻总是小于总传热的时间（13）。因为加和总传热值的时候，只能考虑到前面发生了的碰撞。

有些颗粒甚至15（QSWALL）为0，证明它从未与墙壁接触过。

源代码下载

第1版

9个源文件
链接：https://pan.baidu.com/s/1O0NNRTBh62H40IsNiOa72Q
提取码：abcd

第2版

输出传热所需时间后的源代码（2020-10-29）

链接：https://pan.baidu.com/s/1VODCZ3xwgG3s3Vwkp3hCrg
提取码：nqcs

【读MFiX源代码】MFiX中四种传热方式全面详解（对流、导热、辐射、反应热）并且输出以供后处理（2020-12-15更新）相关推荐

计算机系统的四种启动方式是,详解电脑为何有四种不同的启动设备教程
通常我们开机的时候都是按F2进BIOS的boot菜单设置启动优先级,或者说开机按F12在启动菜单选择启动设备,下面小编就为大家讲解有关启动设备的教程. 首先我们知道电脑是由很多部件组成的,什么CPU, ...
华人科学家发现第四种传热方式！真空声子传热将改写物理教科书
早在上中学时,我们就知道传热一共有三种方式,热传导.热对流.热辐射.如今,这一教科书上的知识要被改写了. 近期,一篇发表在<自然>杂志上的论文让不少传热学.物理学.量子力学等领域的科学家们 ...
计算机曝光模式有哪些,摄影：单反相机中P、A、S、M四种曝光模式的用法详解 -电脑资料...
这篇教程是向脚本之家的朋友介绍单反相机中P.A.S.M四种曝光模式的用法,对于摄影爱好者非常值得学习,推荐到脚本之家,喜欢的朋友一起来看看吧很多朋友在初接触单反相机时对相机的P.A.S.M四种曝光模 ...
python可变参数_Python 的四种共享传参详解
点击上方"Python数据之道",选择"星标公众号" 精品文章,第一时间送达作者 | 杨仁聪编辑 | Lemon 出品 | Python数据之道本文来自公 ...
Plotly中4种文本类型设置详解
作者:Lemon 来源:Python数据之道 Plotly中4种文本类型设置详解大家好,我是 Lemon . 在过去的一段时间里,我写了几篇用 Python 的交互式可视化工具 Plotly 来演示 ...
php读音量大小,Android_Android中实时获取音量分贝值详解，基础知识度量声音强度，大 - phpStudy...
Android中实时获取音量分贝值详解基础知识度量声音强度,大家最熟悉的单位就是分贝(decibel,缩写为dB).这是一个无纲量的相对单位,计算公式如下: 分子是测量值的声压,分母是参考值的声压 ...
C++中四种类型转换方式
C风格的强制类型转换(Type Cast)很简单,不管什么类型的转换统统是:TYPE b = (TYPE)a,但是c 风格的类型转换有不少的缺点,有的时候用c风格的转换是不合适的,因为它可以在任意类型 ...
Burp Suite爆破模块中的四种模式的区别详解和演示(暴力破解)
BrupSuite爆破的四种模式详解最近看了好多关于暴力破解的博客,其中用的最多的工具就是bp了,但是好多都是一上来给了执行步骤,却没有对爆破的这几个模式选择进行解释,所以今天萌新写个纪录,来阐明这 ...
spi协议时序图和四种模式实际应用详解
大家好,我是无际. 上个章节我们讲解了spi接口定义,今天我们更加深入讲解下spi协议时序图和spi四种模式的用法. 刚开始接触单片机开发时,最怕就是看时序图,对于我来说就是奇怪的知识. 特别是SPI ...

【读MFiX源代码】MFiX中四种传热方式全面详解（对流、导热、辐射、反应热）并且输出以供后处理（2020-12-15更新）