FVM in CFD 学习笔记_第15章_流动计算：不可压缩流动_4_SIMPLE家族算法

学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - An Advanced Introduction with OpenFOAM and Matlab
Chapter 15 Fluid Flow Computation: Incompressible Flows

前面几章讲解的关于变量$\varphi$的一般输运方程的离散和求解流程，均是建立在事先已知速度场的前提下。但是一般情况下，速度场是未知的，且需要通过求解Navier-Stokes方程组来获取。对于不可压缩流动而言，该项工作尤为复杂，因为压力和速度是强烈耦合在一起的，而且压力并不以主变量的形式出现在动量或是连续方程中。本章的重点在于展示解决上述两个问题的方法，以及不可压缩流动问题的流场计算方法。先讲解一维交错网格，然后是一维同位网格，最后是同位三维非结构网格。除了阐明SIMPLE、SIMPLEC、PRIME和PISO算法的来龙去脉，还清晰阐明了Rhie-Chow插值方法，以及将其扩展到瞬态、松弛和体积力项的方法。最后，展现了一些常见边界条件的添加细节。

由于本章内容繁杂，篇幅较长，故分成了四部分来讲解，各部分主要内容分别为：交错网格、同位网格、边界条件、SIMPLE家族算法。

这里是第四部分，主要讲解在SIMPLE算法基础上衍生出的其他SIMPLE家族算法、最佳欠松弛因子、Rhie-Chow插值对不同项的处理、代码讲解。

7 SIMPLE家族算法

在SIMPLE算法中，速度和压力是用分离式（顺序式）方式来处理的，压力场是通过压力修正方程算出的，而压力修正方程的推导则是利用离散的动量方程，将连续方程中的速度场替换成压力项来得到的（即，从离散动量方程推导出速度和压力的关系，然后把该关系代入到连续方程中，得到用压力项而非速度项表示的连续方程，即压力修正方程）。在推导过程中，速度修正项，$\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]$，是被忽略掉了的，因为如果保留它的话会使得方程难以处理。

抛掉该项并不会影响到最终的解，因为该项在收敛的时候值为零。只是，这么做影响了收敛的过程，因为在最初的迭代过程中该项的值是相当大的。这个很大的值可能会导致发散或是减缓收敛速度，因为高估了压力修正场。为了平衡该作用，在SIMPLE算法中，压力要做欠松弛处理，即$p=p^{\ast }+{\lambda }^p{p}^{\prime }$，其中${\lambda }^p$是压力欠松弛因子。为了获得最佳收敛效果，${\lambda }^p$通常取成($1-{\lambda }^{\mathbf{v}}$)，其中${\lambda }^{\mathbf{v}}$为动量欠松弛因子，这个后面会详细探讨。

尽管使用了欠松弛处理，SIMPLE算法的收敛速率仍然依据问题的不同而不同，因此学者们寻求各种替代方法来做更好的改进。他们的努力最终发展了一系列SIMPLE类型的家族算法，例如SIMPLEC、SIMPLER、PISO、SIMPLEX、PRIME、SIMPLEM和SIMPLEST。Moukalled和Darwish统一了针对不可压缩和可压缩流动的这些算法的公式，而Darwish等人和Jang等人则评估了它们的性能。这里并不打算给出这些算法的完整描述，而是，聚焦在两种最常用的算法上，它们是Van Doormal和Raithby的SIMPLEC（SIMPLE Consistent）算法，以及Issa的PISO（Pressure-Implicit Split Operator）算法。这两种算法展示了对于$\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]$项的两种不同的处理方法。在SIMPLEC算法中，主网格节点的速度修正近似为 邻近位置处速度修正的加权平均，将$\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]$项转变为修正项，$\widetilde{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]$，其幅值更小（比$\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]$更小），可以忽略。在PISO算法中，$\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]$项将作为分裂算子方法中的一部分来考虑。在其它算法中，SIMPLE算法中直接忽略掉的$\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]$项将进行修改并引入到动量方程或算子${\mathbf{D}}^{\mathbf{v}}$中。因为PISO算法等效于一步SIMPLE算法和一步或多步PRIME算法的结合，所以PRIME算法也将进行介绍。

在PRIME算法中，动量方程是显式求解的。对动量方程的该显式处理是合适的，因为其对整体流场的收敛贡献很小。另一方面，找寻压力场的正确解则是影响整场收敛的重要因素。

在SIMPLEST算法中，动量方程的系数分成对流部分和扩散部分，对流部分显式处理，扩散部分隐式处理，这样影响着${\mathbf{D}}^{\mathbf{v}}$和$\mathbf{H}$。对流部分显式处理的正当性是基于，在纯对流情况下扰动的传播是以幅值不变的有限速度进行的，与，显式迭代方法的单次迭代中从某个特定点到邻近网格点的误差的传播，是相似的。而扩散部分的隐式处理的证明则源自于，在纯扩散情况下扰动的传播是同时向着各个方向进行的且幅值是快速衰减的，与，隐式求解方法中的单次迭代步内误差是在整个计算域上衰减的，是相似的。

在SIMPLEM（SIMPLE-Modified）算法中，压力修正方程是先于动量方程求解的，表明压力场是使用旧的速度场求解的。这样会比速度修正产生更好的压力修正，相较于SIMPLE算法，两者的优点和缺点是互换的。

在SIMPLER（SIMPLE-Revised）算法中，提出了一个额外的方程，从中可直接解出压力，而类似SIMPLE的压力修正方程则用于更新速度场。采用分离的压力方程的原因是，一旦速度场使用预测压力修正场更新后，其不再满足动量方程。因此，压力需要从另一个满足速度的方程来求解，这样可同时满足动量方程。

SIMPLEX算法的发展目的是确保收敛速度不会随着网格的加密而衰减，其与SIMPLE算法的不同是${\mathbf{D}}^{\mathbf{v}}$场是要计算的。这是使用额外的方程组来求解的，其发展和求解所基于的假设是，压力差值的空间分布的影响随着网格的加密变化很小。因此，如果压力差值的影响限制在一个单元上，那么做出下面的假设是合适的，通过外插，主网格节点的压力差值足够代表单元面的压力差值。

尽管所有上述算法最初是用于交错网格的，然而他们在同位网格框架下依然适用。

7.1 SIMPLEC算法

SIMPLEC（SIMPLE-Consistent）算法是SIMPLE算法的修改变型，其改动在于，简单地假设在点$C$处的速度修正是邻近网格节点的修正值的加权平均。数学上表示为
$${\mathbf{v}}_C^{\prime }\approx \frac{\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_F^{\prime }}{\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}}\Rightarrow \sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_F^{\prime }\approx {\mathbf{v}}_C^{\prime }\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}$$
使用$\mathbf{H}$算子，上式可写为
$$\sum _{F\sim NB(C)}\frac{a_F^{\mathbf{v}}}{a_C^{\mathbf{v}}}{\mathbf{v}}_F^{\prime }\approx {\mathbf{v}}_C^{\prime }\sum _{F\sim NB(C)}\frac{a_F^{\mathbf{v}}}{a_C^{\mathbf{v}}}\Rightarrow {\mathbf{H}}_C[{\mathbf{v}}^{\prime }]\approx {\mathbf{v}}_C^{\prime }{\mathbf{H}}_C[1]$$
因此，并非像在SIMPLE中那样忽略掉$\overline{{\mathbf{H}}_C[{\mathbf{v}}^{\prime }]}$项，把它用上式的近似值来代替。这样，速度修正方程${\mathbf{v}}_C^{\prime }+{\mathbf{H}}_C[{\mathbf{v}}^{\prime }]=-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p^{\prime }{)}_C$变为了
$$(1+{\mathbf{H}}_C[1]){\mathbf{v}}_C^{\prime }=-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p^{\prime }{)}_C\Rightarrow {\mathbf{v}}_C^{\prime }=-\widetilde{{\mathbf{D}}_C^{\mathbf{v}}}(\nabla p^{\prime }{)}_C$$
接下来，可用上式推导出压力修正方程。

通过在动量方程中添加和减去$\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_C$项，也可获得同样的结果，即
$$\left(a_C^{\mathbf{v}}+\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}\right){\mathbf{v}}_C+\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}({\mathbf{v}}_F+{\mathbf{v}}_C)=-{\mathbf{v}}_C(\nabla p{)}_C+{\hat{\mathbf{b}}}_C^{\mathbf{v}}$$
可以写成
$${\mathbf{v}}_C+\widetilde{{\mathbf{H}}_C}[\mathbf{v}-{\mathbf{v}}_C]=-\widetilde{{\mathbf{D}}_C^{\mathbf{v}}}(\nabla p{)}_C+{\hat{\mathbf{B}}}_C^{\mathbf{v}}$$
根据上式，速度修正方程变为
$${\mathbf{v}}_C^{\prime }=-\widetilde{{\mathbf{H}}_C}[{\mathbf{v}}^{\prime }-{\mathbf{v}}_C^{\prime }]-\widetilde{{\mathbf{D}}_C^{\mathbf{v}}}(\nabla p^{\prime }{)}_C$$
忽略掉$\widetilde{{\mathbf{H}}_C}[{\mathbf{v}}^{\prime }-{\mathbf{v}}_C^{\prime }]$项，就和前面推导的近似形式一样了，同样用这个修正速度可推导出压力修正方程。

在SIMPLEC算法中，得益于更好的近似估算（即，忽略的项更小），与SIMPLE算法相比，压力的松弛处理不再需要了，所得速度也更加符合动量方程。于是乎，可得到更快的收敛速度。如此看来，忽略掉$\widetilde{{\mathbf{H}}_C}[{\mathbf{v}}^{\prime }-{\mathbf{v}}_C^{\prime }]$项，而非${\mathbf{H}}_C[{\mathbf{v}}^{\prime }]$项，并且把${\mathbf{D}}_C^{\mathbf{v}}$替换成$\widetilde{{\mathbf{D}}_C^{\mathbf{v}}}$，SIMPLEC算法中的步骤和SIMPLE中的步骤是相似的。

7.2 PRIME算法

在PRIME（Pressure Implicit Momentum Explicit）算法中，动量方程是显式求解的，该显式处理是恰当的，因为在迭代求解过程中动量方程对于整体流场的收敛贡献较小，而另一方面，压力场的正确解则对整体收敛的贡献至关重要。基于该认知，将PRIME算法汇总如下：

动量方程显式求解，获得新的速度场${\mathbf{v}}^{\ast }$
$${\mathbf{v}}_C^{\ast }=-{\mathbf{H}}_C[{\mathbf{v}}^{(n)}]-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p^{(n)}{)}_C+{\mathbf{B}}_C^{\mathbf{v}}$$
该速度场将用于推导压力修正方程。定义修正场如下
$${\mathbf{v}}_C^{\ast \ast }={\mathbf{v}}_C^{\ast }+{\mathbf{v}}_C^{\prime }{p}_C^{\ast }=p_C^{(n)}+p_C^{\prime }$$
修正场应满足
$$\begin{array}{rl}{\mathbf{v}}_C^{\ast \ast }& =-{\mathbf{H}}_C[{\mathbf{v}}^{\ast \ast }]-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p^{\ast }{)}_C+{\mathbf{B}}_C^{\mathbf{v}}\\ & =-{\mathbf{H}}_C[{\mathbf{v}}^{\ast }+{\mathbf{v}}^{\prime }]-{\mathbf{D}}_C^{\mathbf{v}}[\nabla (p^{(n)}+p^{\prime }){]}_C+{\mathbf{B}}_C^{\mathbf{v}}\end{array}$$
上式减去上上上式，可得到如下的速度和压力关系
$${\mathbf{v}}_C^{\prime }=\underline{-\left({\mathbf{H}}_C[{\mathbf{v}}^{\ast }-{\mathbf{v}}^{(n)}]+{\mathbf{H}}_C[{\mathbf{v}}^{\prime }]\right)}-{\mathbf{D}}_C^{\mathbf{v}}\nabla p_C^{\prime }$$
上式再代入到连续方程，得
$$-\sum _{f\sim nb(C)}{\rho }_f\overline{{\mathbf{D}}_f}\nabla p_f^{\prime }\cdot {\mathbf{S}}_f=-\sum _{f\sim nb(C)}{\dot{m}}_f^{\ast }+\underline{\sum _{f\sim nb(C)}\left[{\rho }_f\left(\overline{{\mathbf{H}}_f}\left[{\mathbf{v}}^{\ast }-{\mathbf{v}}^{(n)}\right]+\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]\right)\cdot {\mathbf{S}}_f\right]}$$
其中的下划线项可以忽略。

如果${\mathbf{H}}_C[{\mathbf{v}}^{\prime }]$与${\mathbf{H}}_C[{\mathbf{v}}^{\ast }-{\mathbf{v}}^{(n)}]$两项符号相反的话，那么在PRIME算法中忽略的项${\mathbf{H}}_C[{\mathbf{v}}^{\ast }-{\mathbf{v}}^{(n)}]+{\mathbf{H}}_C[{\mathbf{v}}^{\prime }]$是比SIMPLE算法中忽略的项${\mathbf{H}}_C[{\mathbf{v}}^{\prime }]$要来的小的。注意到${\mathbf{H}}_C[{\mathbf{v}}^{\prime }]={\mathbf{H}}_C[{\mathbf{v}}_C^{\ast \ast }-{\mathbf{v}}_C^{\ast }]$是为了满足连续方程而带来的修正，而${\mathbf{H}}_C[{\mathbf{v}}^{\ast }-{\mathbf{v}}^{(n)}]$是为了满足动量方程而带来的修正，通常加到动量方程中的修正与加到连续方程中的修正是符号相反的，因此，PRIME算法中忽略掉的项${\mathbf{H}}_C[{\mathbf{v}}^{\ast }-{\mathbf{v}}^{(n)}]+{\mathbf{H}}_C[{\mathbf{v}}^{\prime }]$是较小的。此外，由于动量方程是显示求解的，所以不需要欠松弛处理，这样带来的好处是增加了算法的稳定性。

7.3 PISO算法

在PISO算法中，修正过程由两步或者更多的步骤构成，而${\mathbf{H}}_C[{\mathbf{v}}^{\prime }]$项则在修正过程的部分步骤中考虑。第一步与SIMPLE算法类似，其中${\mathbf{v}}^{\prime }$是用忽略掉${\mathbf{H}}_C[{\mathbf{v}}^{\prime }]$的方法算出的。满足连续方程的速度${\mathbf{v}}^{\ast \ast }$和压力$p^{\ast }$场将用于重新计算动量方程的系数，并显式求解该方程。新的速度场${\mathbf{v}}^{\ast \ast \ast }$将用于计算单元面处的质量流量场${\dot{m}}^{\ast \ast \ast }$，使用Rhie-Chow插值。然后，${\mathbf{H}}_C[{\mathbf{v}}^{\prime }]$在第二个修正步中将做部分恢复，其中速度修正写为
$$\begin{array}{rl}{\mathbf{v}}_C^{\ast \ast \ast \ast }& ={\mathbf{v}}_C^{\ast \ast \ast }+{\mathbf{v}}_C^{\prime \prime }\\ & =-{\mathbf{H}}_C^{\ast \ast }[{\mathbf{v}}^{\ast \ast }]-({\mathbf{D}}_C^{\mathbf{v}}{)}^{\ast \ast }(\nabla p^{\ast }{)}_C+{\mathbf{v}}_C^{\prime \prime }\\ & =-{\mathbf{H}}_C^{\ast \ast }[{\mathbf{v}}^{\ast }+{\mathbf{v}}^{\prime }]-({\mathbf{D}}_C^{\mathbf{v}}{)}^{\ast \ast }(\nabla p^{\ast }{)}_C+{\mathbf{v}}_C^{\prime \prime }\\ & =-{\mathbf{H}}_C^{\ast \ast }[{\mathbf{v}}^{\ast }]-{\mathbf{H}}_C^{\ast \ast }[{\mathbf{v}}^{\prime }]-({\mathbf{D}}_C^{\mathbf{v}}{)}^{\ast \ast }(\nabla p^{\ast }{)}_C+{\mathbf{v}}_C^{\prime \prime }\\ & =\underset{\approx {\mathbf{v}}_C^{\ast \ast }}{\underbrace{-{\mathbf{H}}_C^{\ast \ast }[{\mathbf{v}}^{\ast }]-({\mathbf{D}}_C^{\mathbf{v}}{)}^{\ast \ast }(\nabla p^{\ast }{)}_C}}-{\mathbf{H}}_C^{\ast \ast }[-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p^{\prime }{)}_C]+{\mathbf{v}}_C^{\prime \prime }\\ & \approx {\mathbf{v}}_C^{\ast \ast }+\underline{{\mathbf{v}}_C^{\prime \prime }-{\mathbf{H}}_C^{\ast \ast }[-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p^{\prime }{)}_C]}\end{array}$$
上式中的下划线部分代表着${\mathbf{H}}_C[{\mathbf{v}}^{\prime }]$的部分，其在第二个修正步中有所恢复。第二个速度修正满足
$${\mathbf{v}}_C^{\prime \prime }=-\underline{{\mathbf{H}}_C^{\ast \ast }[{\mathbf{v}}^{\prime \prime }]}-({\mathbf{D}}_C^{\mathbf{v}}{)}^{\ast \ast }(\nabla p^{\prime \prime }{)}_C$$
使用上上式，以及$C$和$F$节点之间的Rhie-Chow插值，可得到新的压力修正方程
$$-\sum _{f\sim nb(C)}{\rho }_f\overline{{\mathbf{D}}_f}\nabla p_f^{\prime \prime }\cdot {\mathbf{S}}_f=-\sum _{f\sim nb(C)}{\dot{m}}_f^{\ast }+\underline{\sum _{f\sim nb(C)}\left({\rho }_f\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime \prime }]\cdot {\mathbf{S}}_f\right)}$$
其中下划线项可以忽略，该修正步可以想重复多少次就重复多少次，每次恢复${\mathbf{H}}_C[{\mathbf{v}}^{\prime }]$新的部分。

通过对该流程的分析，不难发现，PISO可以视为一个SIMPLE步和一个或多个PRIME步的结合，因此结合了SIMPLE算法的隐式特性和PRIME算法的稳定性。同位网格上的PISO算法的流程可汇总如下：

1. 为计算时刻$t+\Delta t$的解，用时刻$t$的压力、速度和质量流量场$p_C^t$、${\mathbf{v}}_C^t$、${\dot{m}}_f^t$作为初始猜测值$p_C^{(n)}$、${\mathbf{v}}_C^{(n)}$、${\dot{m}}_f^{(n)}$；

SIMPLE 一步

2. 隐式求解动量方程获得新速度场${\mathbf{v}}_C^{\ast }$；
   $${\mathbf{v}}_C^{\ast }+{\mathbf{H}}_C[{\mathbf{v}}^{\ast }]=-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p^{(n)}{)}_C+{\mathbf{B}}_C^{\mathbf{v}}$$
3. 更新单元面处的质量流量，使用Rhie-Chow差值技术获得满足动量方程的质量流量场${\dot{m}}_f^{\ast }$；
   $${\dot{m}}_f^{\ast }={\rho }_f{\mathbf{v}}_f^{\ast }\cdot {\mathbf{S}}_f={\rho }_f\overline{{\mathbf{v}}_f^{\ast }}\cdot {\mathbf{S}}_f-{\rho }_f\overline{{\mathbf{D}}_f^{\mathbf{v}}}\left(\nabla p_f^{(n)}-\overline{\nabla p_f^{(n)}}\right)\cdot {\mathbf{S}}_f$$
4. 使用新的质量流量${\dot{m}}_f^{\ast }$组装压力修正方程，并求解该方程以获得压力修正场$p_C^{\prime }$；
   $$a_C^{p^{\prime }}{p}_C^{\prime }+\sum _{F\sim NB(C)}{a}_F^{p^{\prime }}{p}_F^{\prime }=b_C^{p^{\prime }}$$
   其中
   $$\begin{array}{rl}{a}_F^{p^{\prime }}& =Flux{F}_f=-{\rho }_f{\mathcal{D}}_f\\ a_C^{p^{\prime }}& =-\sum _{f\sim nb(c)}Flux{F}_f=-\sum _{F\sim NB(C)}{a}_F^{p^{\prime }}\\ b_C^{p^{\prime }}& =-\sum _{f\sim nb(C)}Flux{V}_f+\underline{\sum _{f\sim nb(C)}\left({\rho }_f\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]\cdot {\mathbf{S}}_f\right)}\\ & =-\sum _{f\sim nb(C)}{\dot{m}}_f^{\ast }+\underline{\sum _{f\sim nb(C)}\left({\rho }_f\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime }]\cdot {\mathbf{S}}_f\right)}\end{array}$$
5. 基于压力修正场，更新单元形心处的压力和速度场，以及单元面处的质量流量场，以获得满足连续方程的场，这些场分别用$p_C^{\ast }$、${\mathbf{v}}_C^{\ast \ast }$、${\dot{m}}_f^{\ast \ast }$表示；
   $$\begin{array}{rlrl}{\mathbf{v}}_C^{\ast \ast }& ={\mathbf{v}}_C^{\ast }+{\mathbf{v}}_C^{\prime }& {\mathbf{v}}_C^{\prime }& =-{\mathbf{D}}_C^{\mathbf{v}}(\nabla p^{\prime }{)}_C\\ {\dot{m}}_f^{\ast \ast }& ={\dot{m}}_f^{\ast }+{\dot{m}}_f^{\prime }& {\dot{m}}_f^{\prime }& =-{\rho }_f\overline{{\mathbf{D}}_f^{\mathbf{v}}}\nabla p_f^{\prime }\cdot {\mathbf{S}}_f\\ p_C^{\ast }& =p_C^{(n)}+{\lambda }^p{p}_C^{\prime }& & \end{array}$$

PRIME一步或多步

6. 使用最新可用的速度和压力场，计算动量方程的系数并显式求解该方程；
7. 使用Rhie-Chow差值技术，更新单元面处的质量流量；
8. 使用新的质量流量，组装压力修正方程，并求解以获得压力修正场；
   $$-\sum _{f\sim nb(C)}{\rho }_f\overline{{\mathbf{D}}_f}\nabla p_f^{\prime \prime }\cdot {\mathbf{S}}_f=-\sum _{f\sim nb(C)}{\dot{m}}_f^{\ast }+\underline{\sum _{f\sim nb(C)}\left({\rho }_f\overline{{\mathbf{H}}_f}[{\mathbf{v}}^{\prime \prime }]\cdot {\mathbf{S}}_f\right)}$$
9. 更新压力、速度、质量流量场；
10. 返回第6步，依据想要的修正步数，多次重复该流程；
11. 将最新的速度、压力、质量流量值作为初值猜测值；
12. 返回第2步，并重复该流程直到收敛；
13. 将收敛值设为$t+\Delta t$时刻的解；
14. 推进到下一时刻，$t\leftarrow t+\Delta t$；
15. 返回第1步并重复该流程直至到达最后时间步。

PISO算法的流程图如下图所示

8 $\mathbf{v}$和$p^{\prime }$的最佳欠松弛因子

为推动SIMPLE算法的收敛，动量和连续方程要做欠松弛处理，分别使用欠松弛因子${\lambda }^{\mathbf{v}}$和${\lambda }^p$，那么找寻能获得最佳收敛速度的欠松弛因子值就变得很有必要了，回想一下，不含欠松弛处理的速度修正是这样子的
$${\mathbf{v}}_C^{\prime }=-{\mathbf{D}}_C(\nabla p^{\prime }{)}_C$$
此外，在计算压力场时，压力修正要做欠松弛处理，以便让速度修正场满族精确的速度修正方程
$${\mathbf{v}}_C^{\prime }=-{\mathbf{H}}_C[{\mathbf{v}}^{\prime }]-{\lambda }^p{\mathbf{D}}_C(\nabla p^{\prime }{)}_C$$
结合上式和上上式，可得${\lambda }^p$的表达式
$$\begin{array}{rl}& -{\mathbf{D}}_C(\nabla p^{\prime }{)}_C=-{\mathbf{H}}_C[{\mathbf{v}}^{\prime }]-{\lambda }^p{\mathbf{D}}_C(\nabla p^{\prime }{)}_C\\ \Rightarrow & {\lambda }^p=1+\frac{{\mathbf{H}}_C[{\mathbf{v}}^{\prime }]}{{\mathbf{v}}_C^{\prime }}=1+\frac{\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_F^{\prime }}{a_C^{\mathbf{v}}{\mathbf{v}}_C^{\prime }}\end{array}$$
SIMPLEC算法无需对压力修正做欠松弛处理，其有着最佳的加速速度。因此，使用SIMPLEC算法中引入的近似，点$C$处的速度修正可以写成邻近网格节点速度修正的加权平均，即
$${\mathbf{v}}_C^{\prime }\approx \frac{\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_F^{\prime }}{\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}}$$
在前面的推导中，系数$a_C^{\mathbf{v}}$可展开为
$$a_C^{\mathbf{v}}=\frac{1}{{\lambda }^{\mathbf{v}}}\left(a_C^{\mathbf{v}}-\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}+\sum _{f\sim nb(C)}{\dot{m}}_f\right)$$
其在稳态解的限制下，简化为
$$a_C^{\mathbf{v}}=-\frac{1}{{\lambda }^{\mathbf{v}}}\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}$$
上式代入到上上上式，速度修正可近似为
$${\mathbf{v}}_C^{\prime }\approx -\frac{\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_F^{\prime }}{{\lambda }^{\mathbf{v}}{a}_C^{\mathbf{v}}}\Rightarrow a_C^{\mathbf{v}}{\mathbf{v}}_C^{\prime }\approx -\frac{\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_F^{\prime }}{{\lambda }^{\mathbf{v}}}$$
结合前面的式子${\lambda }^p=1+\frac{\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_F^{\prime }}{a_C^{\mathbf{v}}{\mathbf{v}}_C^{\prime }}$，得
$${\lambda }^p\approx 1-{\lambda }^{\mathbf{v}}$$
经验表明使用满足上式的欠松弛因子的SIMPLE算法的性能，与SIMPLEC算法，是相似的。

9 Rhie-Chow插值对各种项的处理

9.1 欠松弛项的处理

基于Rhie-Chow插值，使用同位网格所得到的解是依赖于动量方程中的欠松弛因子的。为了消除该依赖性，需要对Rhie-Chow插值做修改，欠松弛动量方程为
$$\frac{1}{{\lambda }^{\mathbf{v}}}{a}_C^{\mathbf{v}}{\mathbf{v}}_C=-\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_F+{\mathbf{b}}_C^{\mathbf{v}}-V_C\nabla p_C+\left(\frac{1-{\lambda }^{\mathbf{v}}}{{\lambda }^{\mathbf{v}}}\right)a_C^{\mathbf{v}}{\mathbf{v}}_C^{(n)}$$
其中${\mathbf{b}}_C^{\mathbf{v}}$是动量方程中的源项，从中可提取出压力和欠松弛源项，${\mathbf{v}}_C^{(n)}$是单元形心$C$处速度的前次迭代值。对应的欠松弛动量方程，使用交错网格形式，为
$$\frac{1}{{\lambda }^{\mathbf{v}}}{a}_f^{\mathbf{v}}{\mathbf{v}}_f=-\sum _{nb\sim NB(C)}{a}_{nb}^{\mathbf{v}}{\mathbf{v}}_{nb}+{\mathbf{b}}_f^{\mathbf{v}}-V_f\nabla p_f+\left(\frac{1-{\lambda }^{\mathbf{v}}}{{\lambda }^{\mathbf{v}}}\right)a_f^{\mathbf{v}}{\mathbf{v}}_f^{(n)}$$
通过构建单元面上的伪动量方程，Rhie-Chow插值方法模仿了交错网格方程。正是由于该特性，Rhie-Chow才能成功应用。因此作为指导准则，对Rhie-Chow插值的任何改动的衡量标准应该是，改动方程是否和交错网格方程形式类似。因此，使用Rhie-Chow插值的欠松弛方程的形式应该是
$$\frac{1}{{\lambda }^{\mathbf{v}}}\overline{a_f^{\mathbf{v}}}{\mathbf{v}}_f=-\overline{\sum _{nb\sim NB(C)}{a}_{nb}^{\mathbf{v}}{\mathbf{v}}_{nb}+{\mathbf{b}}_f^{\mathbf{v}}}-\overline{V_f}\nabla p_f+\left(\frac{1-{\lambda }^{\mathbf{v}}}{{\lambda }^{\mathbf{v}}}\right)\overline{a_f^{\mathbf{v}}}{\mathbf{v}}_f^{(n)}$$
右端项的第一项平均的获取方法为
$$\begin{array}{rl}-\overline{\sum _{nb\sim NB(C)}{a}_{nb}^{\mathbf{v}}{\mathbf{v}}_{nb}+{\mathbf{b}}_f^{\mathbf{v}}}=& -g_C\left(\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_F+{\mathbf{b}}_C^{\mathbf{v}}\right)-g_F\left(\sum _{N\sim NB(C)}{a}_N^{\mathbf{v}}{\mathbf{v}}_N+{\mathbf{b}}_F^{\mathbf{v}}\right)\\ =& g_C\left[\frac{1}{{\lambda }^{\mathbf{v}}}{a}_C^{\mathbf{v}}{\mathbf{v}}_C+V_C\nabla p_C-\left(\frac{1-{\lambda }^{\mathbf{v}}}{{\lambda }^{\mathbf{v}}}\right)a_C^{\mathbf{v}}{\mathbf{v}}_C^{(n)}\right]+\\ & g_F\left[\frac{1}{{\lambda }^{\mathbf{v}}}{a}_F^{\mathbf{v}}{\mathbf{v}}_F+V_F\nabla p_F-\left(\frac{1-{\lambda }^{\mathbf{v}}}{{\lambda }^{\mathbf{v}}}\right)a_F^{\mathbf{v}}{\mathbf{v}}_F^{(n)}\right]\\ =& \frac{1}{{\lambda }^{\mathbf{v}}}\overline{a_f^{\mathbf{v}}{\mathbf{v}}_f}+\overline{V_f\nabla p_f}-\left(\frac{1-{\lambda }^{\mathbf{v}}}{{\lambda }^{\mathbf{v}}}\right)\overline{a_f^{\mathbf{v}}{\mathbf{v}}_f^{(n)}}\end{array}$$
上式代入到上上式中，可得单元面速度${\mathbf{v}}_f$的扩展Rhie-Chow插值
$${\mathbf{v}}_f=\overline{{\mathbf{v}}_f}-\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\nabla p_f-\overline{\nabla p_f})+(1-{\lambda }^{\mathbf{v}})({\mathbf{v}}_f^{(n)}-\overline{{\mathbf{v}}_f^{(n)}})$$
没有计及欠松弛效应对面速度的影响，所导致的解对欠松弛因子的依赖性。

9.2 瞬态项的处理

当使用后向Euler瞬态格式求解瞬态问题时，离散动量方程为
$$a_C^{\mathbf{v}}{\mathbf{v}}_C=-\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_F+{\mathbf{b}}_C^{\mathbf{v}}-V_C\nabla p_C+a_C^{\circ }{\mathbf{v}}_C^{\circ }$$
其中${\mathbf{b}}_C^{\mathbf{v}}$是动量方程中的源项，从中可提取出压力和瞬态源项。交错网格变量构架下的等效方程形式如下
$$a_f^{\mathbf{v}}{\mathbf{v}}_f=-\sum _{nb\sim NB(f)}{a}_{nb}^{\mathbf{v}}{\mathbf{v}}_{nb}+{\mathbf{b}}_f^{\mathbf{v}}-V_f\nabla p_f+a_f^{\circ }{\mathbf{v}}_f^{\circ }$$
使用Rhie-Chow插值方法，伪单元面方程构建为
$$\overline{a_f^{\mathbf{v}}}{\mathbf{v}}_f=-\overline{\sum _{nb\sim NB(C)}{a}_{nb}^{\mathbf{v}}{\mathbf{v}}_{nb}+{\mathbf{b}}_f^{\mathbf{v}}}-\overline{V_f}\nabla p_f+\overline{a_f^{\circ }}{\mathbf{v}}_f^{\circ }$$
右端项的第一项平均的获取方法为
$$\begin{array}{rl}-\overline{\sum _{nb\sim NB(C)}{a}_{nb}^{\mathbf{v}}{\mathbf{v}}_{nb}+{\mathbf{b}}_f^{\mathbf{v}}}=& -g_C\left(\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_F+{\mathbf{b}}_C^{\mathbf{v}}\right)-g_F\left(\sum _{N\sim NB(C)}{a}_N^{\mathbf{v}}{\mathbf{v}}_N+{\mathbf{b}}_F^{\mathbf{v}}\right)\\ =& g_C\left[a_C^{\mathbf{v}}{\mathbf{v}}_C+V_C\nabla p_C-a_C^{\circ }{\mathbf{v}}_C^{\circ }\right]+\\ & g_F\left[a_F^{\mathbf{v}}{\mathbf{v}}_F+V_F\nabla p_F-a_F^{\circ }{\mathbf{v}}_F^{\circ }\right]\\ =& \overline{a_f^{\mathbf{v}}{\mathbf{v}}_f}+\overline{V_f\nabla p_f}-\overline{a_f^{\circ }{\mathbf{v}}_f^{\circ }}\end{array}$$
可得单元面速度${\mathbf{v}}_f$的扩展Rhie-Chow插值
$${\mathbf{v}}_f=\overline{{\mathbf{v}}_f}-\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\nabla p_f-\overline{\nabla p_f})+\frac{\overline{a_f^{\circ }{\mathbf{D}}_f^{\circ }}}{V_f}({\mathbf{v}}_f^{\circ }-\overline{{\mathbf{v}}_f^{\circ }})$$
没有计及非定常项对面速度的影响，所导致的解对时间步的依赖性，以及小时间步的振荡特性。该修正仅对一阶Euler离散格式有效，对于更加精确的时间离散格式，可遵循同样的原则推导类似的修正。

9.3 体积力项的处理

当在交错网格上处理体积力的时候，体积力项的架构和压力梯度项的架构是一样的。在同位网格情况下，体积力、速度、动量变量是在相同位置上计算的。这样，为了让体积力的离散保持和压力同样的架构，需要把体积力做重新分配。离散动量方程为
$$a_C^{\mathbf{v}}{\mathbf{v}}_C=-\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_F+{\mathbf{b}}_C^{\mathbf{v}}-V_C\nabla p_C+V_C\overline{\overline{{\mathbf{B}}_C^{\mathbf{v}}}}$$
其中双上划线表明两个平均步骤，第一个平均步骤是计算在单元面处的$\overline{{\mathbf{B}}_C^{\mathbf{v}}}$
$$\overline{{\mathbf{B}}_f^{\mathbf{v}}}=g_C{\mathbf{B}}_C^{\mathbf{v}}+(1-g_C){\mathbf{B}}_F^{\mathbf{v}}$$
如下图所示

第二步是对这些面值进行平均，算得单元形心值，如下图所示

单元形心的平均值的推导，最好是通过考虑如下图所示的一维问题来实现。

对于静止流体，压力梯度应该和体积力平衡，推出
$$0=-\nabla p_f+{\mathbf{B}}_f^{\mathbf{v}}$$
展开上述方程，可得在点$C$和$F$压力间的关系
$$p_C=p_F+B_f\delta y$$
或者更加通有的形式
$$p_C=p_F+{\mathbf{B}}_f\cdot {\mathbf{d}}_{CF}$$
其中$B_f$是${\mathbf{B}}_f$的幅值，比如重力加速度情况下是
$$B_f={\rho }_{f}g$$
对于不可压缩流动，体积力中密度随温度的变化是用Boussinesq假设来近似的。

再次对于单元$C$，压力梯度应该等于体积力，即
$$0=-\nabla p_C+{\mathbf{B}}_C^{\mathbf{v}}\Rightarrow \nabla p_C={\mathbf{B}}_C^{\mathbf{v}}$$
然而单元$C$的压力梯度是这样算得的
$$\begin{array}{rl}\nabla p_C& =\frac{\sum _f{p}_f{\mathbf{S}}_f}{V_C}\\ & =\frac{\sum _f(g_C{p}_C+(1-g_C)p_F){\mathbf{S}}_f}{V_C}\end{array}$$
把$p_C=p_F+{\mathbf{B}}_f\cdot {\mathbf{d}}_{CF}$代入，有
$$\begin{array}{rl}\nabla p_C& =\frac{\sum _f(g_C(p_F+\overline{{\mathbf{B}}_f^{\mathbf{v}}}\cdot {\mathbf{d}}_{CF})+(1-g_C)p_F){\mathbf{S}}_f}{V_C}\\ & =\frac{\sum _f{p}_F{\mathbf{S}}_f}{V_C}+\frac{\sum _f{g}_C(\overline{{\mathbf{B}}_f^{\mathbf{v}}}\cdot {\mathbf{d}}_{CF}){\mathbf{S}}_f}{V_C}\\ & =p_F\underset{=0}{\underbrace{\frac{\sum _f{\mathbf{S}}_f}{V_C}}}+\frac{\sum _f{g}_C(\overline{{\mathbf{B}}_f^{\mathbf{v}}}\cdot {\mathbf{d}}_f){\mathbf{S}}_f}{V_C}\\ & =\frac{\sum _f{g}_C(\overline{{\mathbf{B}}_f^{\mathbf{v}}}\cdot {\mathbf{d}}_f){\mathbf{S}}_f}{V_C}\\ & =\overline{\overline{{\mathbf{B}}_C^{\mathbf{v}}}}\end{array}$$
这表明
$$\overline{\overline{{\mathbf{B}}_C^{\mathbf{v}}}}=\frac{\sum _f{g}_C(\overline{{\mathbf{B}}_f^{\mathbf{v}}}\cdot {\mathbf{d}}_f){\mathbf{S}}_f}{V_C}$$
第二个要求是单元面速度与交错网格情况下的方程类似
$$\overline{a_f^{\mathbf{v}}}{\mathbf{v}}_f=\overline{-\sum _{nb\sim NB(C)}{a}_{nb}^{\mathbf{v}}{\mathbf{v}}_{nb}+{\mathbf{b}}_f^{\mathbf{v}}}-\overline{V_f}\nabla p_f+\overline{V_f{\mathbf{B}}_f^{\mathbf{v}}}$$
其中${\mathbf{b}}_C^{\mathbf{v}}$是动量方程中的源项，从中可提取出压力和体积力项。右端项中的第一项为
$$\begin{array}{rl}-\overline{\sum _{nb\sim NB(C)}{a}_{nb}^{\mathbf{v}}{\mathbf{v}}_{nb}+{\mathbf{b}}_f^{\mathbf{v}}}=& -g_C\left(\sum _{F\sim NB(C)}{a}_F^{\mathbf{v}}{\mathbf{v}}_F+{\mathbf{b}}_C^{\mathbf{v}}\right)-g_F\left(\sum _{N\sim NB(C)}{a}_N^{\mathbf{v}}{\mathbf{v}}_N+{\mathbf{b}}_F^{\mathbf{v}}\right)\\ =& g_C\left[a_C^{\mathbf{v}}{\mathbf{v}}_C+V_C\nabla p_C-V_C\overline{\overline{{\mathbf{B}}_C^{\mathbf{v}}}}\right]+\\ & g_F\left[a_F^{\mathbf{v}}{\mathbf{v}}_F+V_F\nabla p_F-V_F\overline{\overline{{\mathbf{B}}_F^{\mathbf{v}}}}\right]\\ =& \overline{a_f^{\mathbf{v}}{\mathbf{v}}_f}+\overline{V_f\nabla p_f}-\overline{V_f\overline{\overline{{\mathbf{B}}_f^{\mathbf{v}}}}}\end{array}$$
可得单元面速度${\mathbf{v}}_f$的扩展Rhie-Chow插值
$${\mathbf{v}}_f=\overline{{\mathbf{v}}_f}-\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\nabla p_f-\overline{\nabla p_f})+\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\overline{{\mathbf{B}}_f^{\mathbf{v}}}-\overline{\overline{\overline{{\mathbf{B}}_f^{\mathbf{v}}}}})$$
其中$\overline{\overline{\overline{{\mathbf{B}}_f^{\mathbf{v}}}}}$为
$$\overline{\overline{\overline{{\mathbf{B}}_f^{\mathbf{v}}}}}=g_C\overline{\overline{{\mathbf{B}}_C^{\mathbf{v}}}}+(1-g_C)\overline{\overline{{\mathbf{B}}_F^{\mathbf{v}}}}$$
上述对单元面速度的额外处理增加了求解过程整体的健壮性，尤其是对于那些体积力的变化影响较大的情况（比如，自由表面流动问题）。

9.4 欠松弛、瞬态和体积力项的联合处理

一般来说，前面三项应该一并处理，这样就要把Rhie-Chow插值改成同时考虑三种效应的形式。幸运的是，使用叠加原理可以很容易地推出如下的界面速度形式
$$\begin{array}{rl}{\mathbf{v}}_f=& \overline{{\mathbf{v}}_f}-\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\nabla p_f-\overline{\nabla p_f})+\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\overline{{\mathbf{B}}_f^{\mathbf{v}}}-\overline{\overline{\overline{{\mathbf{B}}_f^{\mathbf{v}}}}})\\ & +\frac{\overline{a_f^{\circ }{\mathbf{D}}_f^{\circ }}}{V_f}({\mathbf{v}}_f^{\circ }-\overline{{\mathbf{v}}_f^{\circ }})+(1-{\lambda }^{\mathbf{v}})({\mathbf{v}}_f^{(n)}-\overline{{\mathbf{v}}_f^{(n)}})\end{array}$$
其中在计算$\overline{{\mathbf{D}}_f^{\mathbf{v}}}$时，使用的是系数$a_C^{\mathbf{v}}$的欠松弛值。

10 代码讲解

10.1 uFVM

在uFVM中，压力修正方程的实现是在cfdAssembleMdotTerm中进行的，下面给出了算法的核心，其中压力修正方程的系数是通过线性化每个内部面通量来组装的，此外，mdot场（面上的质量流量）是基于下式计算的
$$\begin{array}{rl}{\mathbf{v}}_f=& \overline{{\mathbf{v}}_f}-\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\nabla p_f-\overline{\nabla p_f})+\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\overline{{\mathbf{B}}_f^{\mathbf{v}}}-\overline{\overline{\overline{{\mathbf{B}}_f^{\mathbf{v}}}}})\\ & +\frac{\overline{a_f^{\circ }{\mathbf{D}}_f^{\circ }}}{V_f}({\mathbf{v}}_f^{\circ }-\overline{{\mathbf{v}}_f^{\circ }})+(1-{\lambda }^{\mathbf{v}})({\mathbf{v}}_f^{(n)}-\overline{{\mathbf{v}}_f^{(n)}})\end{array}$$
共分成了9项（即，$\text{I}$到$\text{IX}$），并一步步组装来计算面速度，它们分别是

第$\text{I}$项：插值速度场$\overline{{\mathbf{v}}_f}$；

第$\text{II}$项和第$\text{III}$项：面和平均压力场$-\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\nabla p_f-\overline{\nabla p_f})$；

第$\text{IV}$项和第$\text{V}$项：平均和重新分配体积力$\overline{{\mathbf{D}}_f^{\mathbf{v}}}(\overline{{\mathbf{B}}_f^{\mathbf{v}}}-\overline{\overline{\overline{{\mathbf{B}}_f^{\mathbf{v}}}}})$；

第$\text{VI}$项和第$\text{VII}$项：瞬态通量$\frac{\overline{a_f^{\circ }{\mathbf{D}}_f^{\circ }}}{V_f}({\mathbf{v}}_f^{\circ }-\overline{{\mathbf{v}}_f^{\circ }})$；

第$\text{VIII}$项和第$\text{IX}$项：松弛修正项$(1-{\lambda }^{\mathbf{v}})({\mathbf{v}}_f^{(n)}-\overline{{\mathbf{v}}_f^{(n)}})$。

%
% assemble term I
%     density_f [v]_f.Sf
%
U_bar_f = (dot(vel_bar_f(:,:)',Sf(:,:)'))';
local_FLUXVf = local_FLUXVf + density_f.*U_bar_f;
%
% Assemble term II and linearize it
%     - density_f ([DPVOL]_f.P_grad_f).Sf
%
DUSf = [DU1_f.*Sf(:,1),DU2_f.*Sf(:,2),DU3_f.*Sf(:,3)];
geoDiff =( dot(Sf(:,:)',DUSf') ./ dot(CN(:,:)',Sf(:,:)'))';
local_FLUXCf1 = local_FLUXCf1 + density_f.*geoDiff;
local_FLUXCf2 = local_FLUXCf2 - density_f.*geoDiff;
local_FLUXVf = local_FLUXVf - density_f.*dot(p_grad_f',T')';
%
% assemble term III
%     density_f ([P_grad]_f.([DPVOL]_f.Sf))
%
local_FLUXVf = local_FLUXVf +
    density_f.*dot(p_grad_bar_f(iFaces,:)',DUSf(iFaces,:)')';
%
% assemble terms IV and V
%     density_f [DBVOL]_f.([B]_f -[[B]]_f).S_f
%
local_FLUXVf = local_FLUXVf +
    density_f[iFace]*FVVectorDotProduct(FVTensorVectorDotProduct(DB_f,FVVe
    ctorSubstract(FVMakeVector(bf1_bar_f[iFace],bf2_bar_f[iFace]),FVMakeVe
    ctor(bf1_redistributed_f[iFace],bf2_redistributed_f[iFace]))),S_f) ;
%
% assemble terms VI and VII
%     [Dt]_f (U_Old_f -[v_old]_f.S_f)
%
U_bar_old_f = [velx_old_bar_f[iFace] vely_old_bar_f[iFace])] * S_f';
local_FLUXVf = local_FLUXVf + DT_f*(mdot_old_f[iFace] -
    density_old_f[iFace]*U_bar_old_f);
%
% assemble terms VIII and IX
%     (1-URF)(U_f -[v]_f.S_f)
%
local_FLUXVf = local_FLUXVf + (1.0-Mdot_URF)*(mdot_previous_f -
    density_f.*U_bar_f);
%
% assemble the flow term dot for the face
%
local_mdot_f = local_FLUXCf_1*(pressure[iElement1]+ Pref) +
    local_FLUXCf_2*(pressure[iElement2]+ Pref) + local_FLUXVf;
%
%
% Assemble in Global Fluxes
%
theFluxes.FLUXC1f(iFaces,1) = local_FLUXCf1;
theFluxes.FLUXC2f(iFaces,1) = local_FLUXCf2;
theFluxes.FLUXVf(iFaces,1) = local_FLUXVf;
%
theFluxes.FLUXTf(iFaces,1) = theFluxes.FLUXC1f.*pressureC(iFaces)
    + theFluxes.FLUXC2f.*pressureN(iFaces) + theFluxes.FLUXVf(iFaces);
%
mdot_f = theFluxes.FLUXTf(iFaces);

10.2 OpenFOAM

To be continued…