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

22 篇文章 7 订阅
14 篇文章 22 订阅

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 代码结构

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

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

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

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

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

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

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

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

  1. des_time_march.f(整体控制)
  2. calc_thermo_des.f(整体控制)
  3. des_thermo_newvalues.f(总热源项置0)
  4. des_thermo_conv.f(对流相关)
  5. des_thermo_rad.f(辐射相关)
  6. rxns_gs_des1.f(反应热相关)
  7. des_thermo_cond_mod.f(导热相关)
  8. calc_force_dem.f(导热相关,颗粒之间导热)
  9. 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_Sourcedes_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_Qsdes_usr_var(4,:)
5反应热RXNS_Qsdes_usr_var(5,:)
6辐射传热Qraddes_usr_var(6,:)

注:斜线后面标注的是本地变量,前面是共享变量(封装在module里)。加粗的是我实际使用的变量。

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

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

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

3.1 总传热量

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

  1. 在des_thermo_newvalues.f implicit none语句前添加
  !######CL_200914#######
  USE discretelement, only: des_usr_var
  !######################

如图

在这里插入图片描述

  1. 同一个源文件,在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_Qs
     des_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 算例要求

  1. 打开传热传质和化学反应
  2. 打开des_usr_var并给出数组上限(在这里是6)
  3. 如果考察墙壁导热,应该设置墙壁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_Sourcedes_usr_var(1,:)13
2颗粒颗粒导热QQ_TMP(本地)des_usr_var(2,:)14
3颗粒墙壁导热QSWALL(本地)des_usr_var(3,:)15
4对流传热CONV_Qsdes_usr_var(4,:)16
5反应热RXNS_Qsdes_usr_var(5,:)17
6辐射传热Qraddes_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

  • 12
    点赞
  • 11
    收藏
    觉得还不错? 一键收藏
  • 2
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值