【MFiX有反应算例设置】1 usr_rates.f和usr_rates_des.f（以tutorial的silane_pyrolysis_3d为例）

目录：一共四步

概述
1 导入模块部分（46-69行）
2 声明当地变量（73-120行）
3 给当地变量赋值（130行-170行）
4 给速率赋值（171行-结束）
后处理里看到反应

概述

以19.3.1版本为例
在MFiX里设置化学反应，有这么三步

在GUI（.mfix文件，或者旧版本里的mfix.dat）里设置反应方程
在usr_rates.f和usr_rates_des.f里设置速率
在GUI里编译并运行算例

第一步和第三步都不用说了，很简单。经常出问题的是第二步，因为usr_rates.f和usr_rates_des.f是要自己用Fortran语言写代码的。

usr_rates_des.f是用来编写异相反应的速率的
usr_rates.f是用来编写同相反应的速率的

两者差异不太大，先从难的开始讲，也就是usr_rates_des.f

很遗憾目前DEM的tutorial里面没有带反应的流动，所以我们以PIC的tut为例（另外注意，与DEM不同的是，TFM中所有的反应都在usr_rates.f文件里编写）

参考的算例在

<mfix的源文件所在地址>\mfix-19.3.1\tutorials\pic\silane_pyrolysis_3d

参考的例子是
silane_pyrolysis_3d

所用的看代码的工具是Visual studio code
该usr_rates_des.f一共有225行

除去前面的注释部分，从
在这里插入图片描述
开始，到

结束

一共可以分为四部分看

所导入的模块声明
所定义的局部变量声明
与速率相关参数的赋值
速率的赋值

1 导入模块部分（46-69行）

在这里插入图片描述首先看子程序的函数头。

这个子程序的名字叫做USR_RATES_DES
一共从外部引入了四个参数：

NP 代表颗粒的编号
pM 代表颗粒所在固体相的编号（固体相可以有很多个，从1开始，在solids选项卡中定义了几种颗粒就有几种。流体相只有一个，编号为0）
IJK 代表流体网格的编号
DES_RATES 代表异相反应的速率。这个子程序的最终目标就是给这个变量赋值。

继续往下看

往下有很多use，每use一个，代表引入了一个module。具体module是什么去学习一下Fortran语言就知道了，和Python的 module是差不多的概念。

由于module很多，这里就挑选几个讲解一下含义。以后有时间会再写篇文章专门分析每个module的作用。

compar，用于比较数字。我们都知道计算机里面浮点数都是不精确的，不能用==来比较，compar就是干这个用的。
constant 定义一些常数。比如通用气体常数，pi，代表极其小的数字如1e-10，代表极其大的数字如10^10(这里随便说说，不一定是这个数)之类的。
des_thermo 与颗粒相关的热物性变量，比如des_T_s, 表示颗粒的温度，des_C_ps表示颗粒的比热容。
discretelement 还是颗粒相关的参数，只不过更基础一些。比如pmass，代表颗粒的质量，就在这个module里
fldvar 全名我猜是field variables，与场相关的变量。这里定义的应该都是宏观的量。P_g即气体压力，rop_g即气体密度，ep_g即气体体积分数。

注：再次注意这个模块里定义的是宏观的变量，比如x_s虽说是代表固体的某组分分数，但是不能直接用在某个颗粒身上，如果你在运行中输出这个变量观察它，就会发现它永远都是0。我的理解是这个x_s是用来后处理用的。真正要用某个颗粒的某组分分数，要用des_x_s这个变量。

physprop 这个模块包含了很多东西，都是与物性相关的。全名应该是physical properties
其他的模块以后再分析，但是看名字应该就能大致知道是干什么的。比如parallel是用来并行的，rxn应该与反应有关，des_rxn应该与颗粒反应有关。

2 声明当地变量（73-120行）

在这里插入图片描述首先implicit none表示不使用Fortran默认定义变量的惯例
INTEGER INTENT(IN)有三行，是给引入的三个整型参数下定义。IN表示只读，不能更改。
DES_RATES则是用于输出（OUT）的变量，能读能写。并且它是一个双精度浮点数数组，数组大小是
NO_OF_DES_RXNS
很明显这个变量代表颗粒反应的数目。

要特别注意的是，颗粒反应的单位是kmol/s !!

2020/4/17 勘误
之前说的

“要特别注意的是，颗粒反应的单位是mol/s !! 这与气相反应时不同的，气相是kmol/s”

是错误的!
现已改正

在MFiX19.2以后，usr_rates.f和usr_rates_des.f中的单位均改为了kmol/s
旧的tutorial里面的算例使用的CGS单位制，所以用的是mol/s
而现在的版本里全部都是SI单位制，全都改为了kmol/s

接着往下看
Dp0代表颗粒初始直径
xTs xTg分别代表颗粒的和气体的受限后的温度
Pg_KPa代表气体压力转换成kPa 没啥乱用
c_SiH4, c_SiH2, c_Si2H6, c_H2是四种气体的摩尔浓度
p_SiH4, p_SiH2, p_H2是四种气体的分压 kPa
Sa是颗粒表面积，但是注释说错了，它说是单位体积的表面积

以下几个都和具体的反应有关，不具有通用性
DIFF_SiH2 SiH2的扩散速率（什么是扩散速率，看《燃烧学》）
R_SiH2 SiH2的气体常数
K_f 膜扩散阻力
k_so 一阶速率常数
K_H氢附着常数 Ks硅烷附着常数

FWD, RVS, NET分别是正反应，逆反应，净反应速率，可逆反应会用到

继续往下看

 DOUBLE PRECISION, parameter :: c_Limiter = 1.0d-7

这句是给摩尔浓度加一个限制，假如反应物浓度太低，就认为这种物质是不存在的，也就不会去计算反应，减少了计算量。parameter代表它是常变量，一旦给定就不能变了。

INCLUDE 'species.inc'

这句是引入species.inc这个头文件，相当于把这个头文件的内容直接复制粘贴到这里。这个头文件里都是啥呢？其实很容易看见，打开你算例文件夹下的build文件夹（编译之后出现），会发现其实是用来给组分名，反应名编号用的。
比如我随便打开一个我之前的算例（不是当前解读的算例）下的build/species.inc
在这里插入图片描述可以看见其实就是给各个你定义的组分，反应一个编号。这也是为什么user guide上说，你一旦变更了组分和反应，就一定要重新编译的原因。因为这个编号变了。

3 给当地变量赋值（130行-170行）

在这里插入图片描述
前两个赋值不用多说，分别是通过半径计算直径和表面积
xTs和xTg说明一下：

  xTg = max(min(T_g(IJK),   TMAX), 300.0d0)
  xTs = max(min(DES_T_s(NP), TMAX), 300.0d0)

这两句很显然是限制气体和颗粒温度不能超过TMAX，不能低于300K的
TMAX就是你定义组分的Cp的时候，所给的温度上限，默认是4000K

问题来了，为什么要限制温度呢？为什么不直接用温度呢？

之前我也不明白，但是后来我发现这与帮助收敛有关系。有反应的流动，尤其是在颗粒上，由于会放出大量的热，会导致局部的温度很高（我之前的一个例子里，其他地方都是1000K左右，颗粒却超过2700K）
假如温度太高了，超出Cp所定义的范围了，那我估计这个时候Cp会随机给一个不知道是什么的数，发散是必然的。

继续往下
Pg_KPa不用说了

下面这几句说一下

  c_SiH4  = ZERO;   p_SiH4 = ZERO;
  c_SiH2  = ZERO;   p_SiH2 = ZERO;
  c_H2    = ZERO;   p_H2   = ZERO;
  c_Si2H6 = ZERO;

这其实是与下面的

 IF(X_g(IJK,SiH4) > c_Limiter) then

等语句相对应的

之前说过 c_Limiter是用来减小计算量的。假如组分的质量分数低于这个c_Limiter，就认为它没有。

为什么要赋初值0呢？
因为一旦不满足(X_g(IJK,SiH4) > c_Limiter这个条件，计算机有可能会随便瞎给一个数！这就会导致错误了。

（官方这里又出一个小错误，x_g应该是质量分数，但是和浓度去比较了，不过无伤大雅，c_Limiter就是一个很小的数而已）

! Calculate the partial pressure and molar concentration of SiH4
      IF(X_g(IJK,SiH4) > c_Limiter) then
         p_SiH4 = Pg_KPa * X_g(IJK,SiH4) * (MW_MIX_g(IJK) / MW_g(SiH4))
         c_SiH4 = RO_g(IJK) * X_g(IJK,SiH4) / Mw_g(SiH4)
      ENDIF

先算分压，用总压力乘以该组分的摩尔分数。其中摩尔分数是通过质量分数乘以混合物的分子量与该物质分子量的比得到的。（MW代表Molecular Weight 分子量）
摩尔浓度，用气体混合物密度乘以该组分质量分数，再除以分子量

后面几行代码都类似

4 给速率赋值（171行-结束）

在这里插入图片描述注释写的很清楚，总共就俩异相反应：

! SiH4 --> Si + 2H2  (kg-mole/m^3.s)
! SiH2 --> Si + H2  (kg-mole/m^3.s)

具体反应速率怎么给的不说了，因为不具有代表性。
说两点

IF ELSE语句是干什么用的
des_rates这个数组

首先 IF ELSE很显然和之前一样，具有减少计算量的作用（否则即使很少量的反应物也要进行，因为舍入误差的原因可能永远无法终止反应）

else 赋0的作用也是，就是要让反应及时停止，否则会无止境的计算下去。

其次，des_rate这个数组。
它的大小是NO_OF_DES_RXNS个双精度数
每个双精度数存储的是一个异相反应的速率
下标则由species.inc里面的编号给定
比如DES_RATES(RX3)代表的是RX3这个反应的速率
而RX3其实是被species.inc里面的编号代替的。
RX3就是1
RX4就是2
(异相反应和同相反应的编号是互不干扰的)

后处理里看到反应

最后额外说一点官方算例里面没有的：

要想在后处理里看到反应，光有DES_RATES(RX3)是不够的！

即你此时打开mfix GUI，观察rrates数组，还是不能看到任何反应！

因为des_rates数组是默认不会输出出去的！
ReactionRates这个数组才能用于后处理！

所以还要在return之前加上这么几句

   LB= NO_OF_RXNS+1     
   UB= NO_OF_RXNS+NO_OF_DES_RXNS
   ReactionRates(IJK,LB) =  ReactionRates(IJK,LB)     + DES_RATES(RX3)
   ReactionRates(IJK,LB+1) =  ReactionRates(IJK,LB+1) + DES_RATES(RX4)

当然，在前面当地变量声明的位置，你还要声明一下LB和UB

integer LB,UB

ReactionRates数组包括同相和异相反应，所以要重新编号一下。LB就代表气相反应之后的第一个下标，也就是异相反应存储的第一个下标

为什么要叠加上去而不是直接赋值DES_RATES(RX3)呢？

因为这个USR_RATES_DES子程序，是针对每个单独的颗粒的！而一个网格内可能有很多颗粒都发生反应，所以要叠加！

结束！收工！

posted @ 2020-04-16 18:16 北斗123 阅读(767) 评论(0) 编辑收藏举报

刷新页面返回顶部

北斗123