【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里设置化学反应，有这么三步

1. 在GUI（.mfix文件，或者旧版本里的mfix.dat）里设置反应方程
2. 在usr_rates.f和usr_rates_des.f里设置速率
3. 在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	utorialspicsilane_pyrolysis_3d

参考的例子是
silane_pyrolysis_3d

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

除去前面的注释部分，从

开始， 到

结束

一共可以分为四部分看

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

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)

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

1. IF ELSE语句是干什么用的
2. 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子程序，是针对每个单独的颗粒的！而一个网格内可能有很多颗粒都发生反应，所以要叠加！

结束！ 收工！