【读MFiX源代码】3 run_fluid计算流体相 (主要步骤DO_ITERATION) (源文件iterate.f)

核心内容请直接跳到DO_TERATION节

文件层次

位置

文件位置: model/iterate.f

模块：module iterate

被谁调用/入口：run_mfix -> run_fluid (mfix.f)-> do_iteration

本文绝大部分内容（DO_ITERATION）均是下图186行这一行所运行的内容

注意观察，每次都是按一次F5（或者说F10）

运行后结果：产生一行残差输出（如下图）

内部结构

其中最主要的部分就是定义了三个子程序： ITERATE_INIT, DO_ITERATION和POST_ITERATION

我们的重点放在DO_ITERATION，简要介绍其他俩子程序

文件解读

ITERATE_INIT（118行-200行）

主要进行该时间步流体相计算的初始化（不是整个程序的初始化）

其中

! initializations
      DT_prev = DT
      NIT = 0
      MUSTIT = 0
      CONVERGED = .FALSE.
      DIVERGED  = .FALSE.
      RESG = ZERO
      RESS = ZERO

是初始化改时间步流体相所需要的全局变量。DT_prev应该是预测的下一个时间步，DT应该就是已经计算出来的上一个时间步。
NIT是目前流体相计算循环的步数
MUSTIT是判断收敛与否的一个flag，0代表收敛，1是未收敛，2是已经发散

! flag indicating convergence status with MUSTIT = 0,1,2 implying
! complete convergence, non-covergence and divergence respectively
      INTEGER :: MUSTIT

RESG与RESS分别是气相压力残差和固相压力残差

! gas & solids pressure residual
      DOUBLE PRECISION :: RESg, RESs

      IF(NORM_G == UNDEFINED) THEN
      以下忽略直到186行
......

133-186基本上可以忽略
133-147行几个IF可以忽略，基本上是用来残差归一化的。
往下调用了一个函数CALL INIT_RESID ()，按照注释是初始化残差的
再往下 IF(CYCLIC) CALL GoalSeekMassFlux也忽略，和循环边界条件有关
再往下一个大的条件嵌套IF (FULL_LOG) THEN，也忽略，和输出日志有关
再往下CALL CALC_RESID_MB(0, errorpercent)应该也和残差有关，忽略

一直到186行

! Calculate the face values of densities and mass fluxes for the first
! solve_vel_star call.
      CALL CONV_ROP()
      CALL CALC_MFLUX ()
      CALL SET_BC1
      CALL SET_EP_FACTORS

这几步用来计算V*，V*应该和simple算法有关系。是重点
注释说是用来为第一次计算vel_star而计算面上的密度和质量流量的

! JFD: modification for cartesian grid implementation
      IF(CARTESIAN_GRID) CALL CG_SET_OUTFLOW

不懂，但大致是处理笛卡尔网格的，应该比较重要

! Default/Generic Error message
      lMsg = 'Run diverged/stalled'

      RETURN
      END SUBROUTINE ITERATE_INIT

这句证明lMsg这个字符串数组变量是用来报告是否发散的，默认是给发散的。
结束

总结：主要就做了两件事：1.给一些全局变量赋0， 2.为第一次计算vel_star而计算面上的流量和密度

DO_ITERATION

从头截止到calc_coeff（211-242）

函数头
先看函数头，只引进了一个变量，就是MFIX_DAT，而且是只读的，这个变量是从.mfx文件读取而来的字符串数组

这个局部变量整数M目前不知道是干什么的，往后看

phip是Johnson & Jackson壁面边界条件所用到的，部分滑移的壁面边界，我们这里统一用no slip wall 因此不涉及。跳过

往下的IF (.NOT.SETG) THEN等依然忽略

USR_2
这个比较重要，我们知道了USR2这个子程序是在计算流体相的每个时间步的最开始调用的。USR2又是我们可以自定义的子程序。我们可以在此处定义一些与流体相有关的变量。比如当前时刻的Re，当前时刻的CFL数等。 当然这些实际上都是用上一时刻得到的值计算出来的，因为当前时刻的计算刚开始。

CALC_COEF
这个也比较重要，注释说的比较简单，是计算除了密度和反应之外的一些参数的。
如果计算出现问题，报错管理器就会使程序强行终止（return)

是什么参数还要具体去看子程序CALC_COEFF(MFIX_DAT, IER, 1)内容
位置:calc_coeff.f
模块：MODULE CALC_COEFF_MOD
子程序 SUBROUTINE CALC_COEFF(MFIX_DAT, IER, pLevel)
函数头注释
! Subroutine: CALC_COEFF !
! Purpose: This routine directs the calculation of all physical and !
! transport properties, and exchange rates.
用以直接计算所有物理和输运特性，并计算动量和能量交换
内部主要调用的函数为三个
PHYSICAL_PROP
TRANSPORT_PROP
EXCHANGE
如下图

具体信息按下不表

从SCLAR_PROP到CALC_RRATE：SIMPLE算法的前期准备

从这里开始往后，应该进入了SIMPLE算法的步骤，是核心的部分

往下每一步都是调用一个子程序

本节是SIMPLE算法的前期准备步骤

SCLAR_PROP计算用户自定义标量的扩散系数和源项，跳过
K_Epsilon_PROP计算K_Epsilon方程的扩散系数和源项，跳过
CALC_TRD_AND_TAU为MMS算例（？）更新应力张量（？），跳过

SOLVE_VEL_STAR计算v，重点*

DADIAL_VEL_CORRECTION与柱坐标有关，跳过

PHYSICAL_PROP计算密度，重点
并且如果报错会利用错误管理器返回错误值并直接返回

CALC_RRATE计算化学反应，重点

从 IF(.NOT.(DISCRETE_…到该ENDIF(269-309)：无用，直接跳过

为了解非DEM算例（如TFM等）的密填充而解体积分数的方程

该部分直接跳过

从CALC_P_STAR的IF到 k_ESPSILON的IF（312-341）：SIMPLE核心步骤，重中之重

第一个IF(关于calC_P_STAR的) 也是TFM等使用的，跳过

CONV_ROP 计算面上的气体密度值，重点，SIMPLE算法相关

下一个IF(RO_G0/=0)，极其重点，SIMPLE相关
其中
RO_G0是设定的气体恒定密度，一般不设常密度，默认用理想气体公式

下一个有关SOLVE_PP_G的IF，SIMPLE算法核心步骤之一
用于计算压力修正方程，

CORRECT_0()，也是SIMPLE核心步骤之一
用于修正压力，速度和密度

PHYSICAL_PROP，SIMPLE核心步骤之一
重新计算密度，并且如果出错用错误管理器直接return

从SET_WALL_BC到SOLVE_K_EPSILON_EQ(336-372)：SIMPLE的后续步骤

该节也较为重要

SET_WALL_BC，设定壁面边界，重要
注释中的NSW指的是no slip wall，即不设定默认为no slip wall

CALC_MFLUX
SET_BC1
SET_EP_FACTORS
这三步与质量流量边界有关，重要

下面的IF(CARTESIAN_GRID)跳过，与笛卡尔坐标系有关

SOLVE_ENERGY_EQ显然是解能量方程的，并附加错误管理器
值得一提的是可见ENERGY_EQ这个变量是用来控制能量方程开关的。应该能在deck file中指定

GRANULAR_ENERGY与TFM有关，跳过

SOLVE_SPECIES_EQ解组分方程，重要
附加了错误管理器IER_MANAGER()

SOLVE_SCALAR_EQ用于计算其他标量的输运，跳过

SOLVE_K_EPSILON_EQ计算湍流相关，目前先用不到

从IF (.not. cyclic)到contains(基本算作结束）：非重要后处理步骤，终端输出残差信息

IF (.not. cyclic) 用户自定义线性方程，目前跳过

下面的几步是检查是否收敛，残差

再下面是展示残差，DISPLAY_RESID()
终端输出的残差信息就是在这一步输出的
经过这一步之后，会输出如下

从此可以看出，每个DO_ITERATION只是输出了残差输出的一行，NIT的一次计算

最后调用了END_ITERATION

END_ITERATION和IER_MANAGER(定义在DO_ITERATION内部的子程序)

简单说下
子程序定义在DO_ITERATION内部，
意味着只有DO_ITERATION能调用这两个小工具

END_ITERATION（判断收敛，仅最后调用）

只在DO_ITERATION最后调用了
主要用于判断是收敛、发散还是还未收敛

IER_MANAGER（解方程出错时返回错误代码）

经常跟在解方程的后面，用于返回错误代码并且直接return掉DO_ITERATION
这是一个逻辑函数，就返回TRUE和FALSE两个值，true代表计算方程时有错误了

本质上就是一个很长的IF ELSE语句

POST_ITERATE（写出收敛/发散日志，错误信息）

很短的一个子程序，其作用就是写出日志，错误信息等

IER=1代表没发散也没收敛

收敛与否的信息输出是这个子程序处理的