FVM in CFD 学习笔记_第14章_源项离散，松弛，其它细节

最新推荐文章于 2022-08-13 20:42:38 发布

梅冠华

最新推荐文章于 2022-08-13 20:42:38 发布

阅读量2.7k

点赞数 5

分类专栏： FVM in CFD 学习笔记文章标签：源项欠松弛处理收敛判据残差方程均方根残差

本文链接：https://blog.csdn.net/meiguanhua/article/details/105480945

版权

FVM in CFD 学习笔记专栏收录该内容

14 篇文章 129 订阅

订阅专栏

学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - An Advanced Introduction with OpenFOAM and Matlab
Chapter 14 Discretization of the Source Term, Relaxation, and Other Details

本章咱们讲讲那些可以对解的特性产生“较大”影响的“细小”数值细节。首先，考察源项的处理，在通常状况下它是依赖于解的，即 $Q^\phi=Q^\phi(\phi)$ ，源项将线性化处理成因变量的形式，并分成两部分，一部分显式处理，另一部分隐式处理。紧接着，讨论对代数方程组做欠松弛处理的显式和隐式技术，并展现一些隐式欠松弛方法，有著名的Patankar隐式欠松弛方法、van Doormaal和Raithby的E-因子方法、Mallinson和de Vahl Davis的伪瞬态方法。然后，介绍离散代数方程组的残差形式。最后是评估解收敛与否的收敛指标。

1 源项离散

源项（sink and source，源和汇）出现在众多流动和输运问题的控制方程中。例如，湍流模型、化学反应、辐射和传热、质量传递、多相流等，不逐一列举。这些源项不仅影响着问题的物理特性，同时也影响着计算的数值稳定性。如果处理得当，源项会提高健壮性（robustness）。通常建议把负值（sinks，汇）做隐式处理，把正值（sources，源）做显式处理。

为厘清源项的处理方法，考虑在体积为 $V_C$ 的单元 $C$ 上含源项（显式表达）的关于变量 $\phi$ 的通有守恒方程的离散形式
在这里插入图片描述

$a_C\phi_C+\sum_{F\sim NB(C)}a_F\phi_F=Q_C^\phi V_C$
其中 $Q_C^\phi V_C$ 代表在单元 $C$ 上的源项积分。

通常来说，源项是因变量 $\phi$ 的函数，其函数关系为
$Q_C^\phi = Q(\phi_C)$
在这种形式中，源项可基于可用的 $\phi$ 值来显式计算，即，在迭代过程中的上一次迭代值。然而，如果 $Q_C^\phi$ 的值为常数或相对较小时该方法尚可接受，倘若若 $Q_C^\phi$ 的变化跟方程中的其他项相比较大时，收敛速率将会受到极为不利的负面影响。此时，可采用类似Taylor级数展开的方法对 $Q_C^\phi$ 线性化处理，以提高其收敛速率。把前次迭代的值用上标 $^*$ 标识，当前迭代的源项 $Q_C^\phi$ 值可展开为
$\begin{aligned} Q(\phi_C) &=Q(\phi_C^*)+\left( \frac{\partial Q}{\partial \phi_C}\right)^*(\phi_C-\phi_C^*) \\ &= \underbrace{\left( \frac{\partial Q}{\partial \phi_C}\right)^*\phi_C}_{\text{Implicit~part}}+ \underbrace{Q(\phi_C^*)-\left( \frac{\partial Q}{\partial \phi_C}\right)^*\phi_C^*}_{\text{Explicit part}} \end{aligned}$
其中的显式部分是基于前次迭代值算出的，隐式部分则并入对角系数中，于代数方程组中求解。

此时，代数方程 $a_C\phi_C+\sum_{F\sim NB(C)}a_F\phi_F=Q_C^\phi V_C$ 的右端项变为
$\begin{aligned} Q_C^\phi V_C &= \iint_{V_C}Q^\phi dV \\ &= \iint_{V_C}\left( \frac{\partial Q_C^*}{\partial \phi_C}\phi_C \right)dV +\iint_{V_C} \left( Q_C^*-\frac{\partial Q_C^*}{\partial \phi_C}\phi_C^*\right) dV \\ &=\left( \frac{\partial Q_C^*}{\partial \phi_C}V_C \right)\phi_C + \left( Q_C^*-\frac{\partial Q_C^*}{\partial \phi_C}\phi_C^*\right)V_C \\ &= FluxC_C\phi_C+FluxV_C \end{aligned}$
则离散代数方程变为
$(a_C-FluxC_C)\phi_C+\sum_{F\sim NB(C)}a_F\phi_F=FluxV_C$
在该方程中，源项的隐式部分， $FluxC_C$ 需要是负数，以保证对角占优或Scarborough准则，若不满足会造成计算发散。此外，当变量 $\phi$ 为正定的时候，显式部分， $FluxV_C$ 必须是正的，以保证预测出正的 $\phi$ 值。

例1
在涉及辐射传热的问题中，源项在能量方程中的形式为
$Q^T=A(T_\infin^4-T^4)$
其中 $A$ 为常数， $T$ 是任意网格点处的温度， $T_\infin$ 代表不变的环境温度。对该源项在体积为 $V_C$ 的单元 $C$ 上做体积分，然后使用不同的方法来线性化处理，并解释他们在收敛性能上的差别。

解
$\int_{V_C}Q^TdV=Q_C^TV_C=FluxV_C+FluxC_CT_C$
一些选项可被用于线性化 $Q^T$

选择1
$\begin{aligned} & FluxC_C=0 \\ & FluxV_C=A(T_\infin^4-T_C^4)V_C \end{aligned}$
采用这种方法的话不是很好，因为，当 $T_C>T_\infin$ 的时候， $FluxV_C$ 是负值，而前面讲过，如果 $\phi$ 是正值（这里的温度 $T$ 必须是正值）那么 $FluxV_C$ 也得是正值才对的，所以 $FluxV_C$ 是负值的话，会在迭代过程中产生不符合物理意义的负的温度值，还会造成解的发散。

选择2
以前次迭代的温度值 $T_C^*$ 为参考做线性展开，将源项写成
$\begin{aligned} Q_C^T &=(Q_C^T)^*+\left( \frac{dQ_C^T}{dT_C} \right)^*(T_C-T_C^*) \\ &= A(T_\infin^{*4}-T_C^{*4})-4AT_C^{*3}(T_C-T_C^*) \end{aligned}$
于是， $FluxC_C$ 和 $FluxV_C$ 为
$\begin{aligned} & FluxC_C=-4AT_C^{*3}V_C \\ & FluxV_C=A(T_\infin^4+3T_C^{*4})V_C \end{aligned}$
这才是理想的方法啊，其 $FluxV_C$ 是正值，这样迭代时就不会出现负的温度值，而 $FluxC_C$ 是负值，可保证最好的收敛速率，惬意。

2 代数方程组的欠松弛处理

前面章节中一再讲到，离散过程最终生成的是一个线性代数方程组，如下图，其形式为 $a_C\phi_C+\displaystyle\sum_{F\sim NB(C)}a_F\phi_F=b_C$ ，其中 $a_F$ 为邻居系数，代表邻近变量 $\phi_F$ 对单元变量 $\phi_C$ 的影响， $b_C$ 是方程右端项，常代表源项和其它变量的影响，而 $a_C$ 为代数方程的主系数，其包含了众多影响，如空间离散效应、瞬态效应等。该方程通常是对角占优的。

在这里插入图片描述

在该代数方程组的迭代求解过程中，经常希望减慢因变量在迭代过程中的变化。这在提高非线性问题的收敛性上是必须的，同时也是为了避免从一个猜测的初始场（有可能跟解的差距非常远）开始迭代所造成的发散。由于网格系统的非正交性、源项的存在、模型方程的非线性属性，等，会产生各种非线性特性。经常用来促进收敛性的一种方法是松弛方法，其通过“减慢”（“松弛”）待求解变量的（有时候过于剧烈的）变化来实现。在许多CFD代码中使用的标准松弛方法是Patankar的隐式欠松弛方法，其在第8章中已经简要介绍过了。其他欠松弛方法有E-Factor松弛方法，伪瞬态方法，等。Van Doormaal和Raithby表明，这些不同的松弛方法是有某些相关性的，而且这些方法中任意一个的欠松弛与其他方法中的欠松弛也是相关的，因为它们本质上都是减缓了邻近单元和源项对欠松弛单元值的影响。换句话说，欠松弛均等地影响着相关单元的源项和空间系数，一些松弛方法展示如下。

2.1 欠松弛方法

松弛方法的实现，可通过在获得每次迭代解后，显式处理，也可以在解获得前，隐式地将其影响添加到方程中去。两种方法都将介绍。

2.2 显式欠松弛

在显示欠松弛方法中，每次迭代的最后，当获得新解后，访问计算域中的所有单元，并把任意单元 $C$ 上的预测值 $\phi_C^{new,predicted}$ 按下式进行修正
$\phi_C^{new,used}=\phi_C^{old}+\lambda^\phi(\phi_C^{new,predicted}-\phi_C^{old})$
其中 $\lambda^\phi$ 为松弛因子，对于显式或者隐式松弛皆有该因子存在，其符号不同代表的意义也不同。

$\lambda^\phi<1$ 代表欠松弛，这会减缓收敛的速度，但是会增加求解的稳定性，即，其减弱了解中可能出现的发散或震荡。
$\lambda^\phi=0$ 代表无松弛，迭代中的预测值不加修正，直接用于下次迭代。
$\lambda^\phi>1$ 代表超松弛，某些情况下其可用于加速收敛，但是通常会降低计算的稳定性（欲速则不达，过犹不及）。

显式欠松弛将用于SIMPLE算法中的压力欠松弛处理，这将在下一章详细讲解。还有，在流动特性依赖于解且需要迭代更新的情况下，显式欠松弛也许是促进收敛的必要操作，这些例子包含但不限于，湍流流动中的湍流粘性、可压缩流动中的粘性、使用高阶格式计算的界面值，等。此外，其也可用于对守恒方程中的单独项（比如源项）进行欠松弛处理，以及某些情况下对求解变量的梯度进行欠松弛处理。

2.3 隐式欠松弛

Patankar的标准欠松弛方法在第8章讲过了，但是为了完整起见，这里给出其表达式。此外，还将讨论E-factor方法和伪瞬态技术。

2.3.1 Patankar欠松弛

如前所述，可通过引入松弛因子 $\lambda^\phi$ ，并通过表达式 $\phi_C^{new,used}=\phi_C^{old}+\lambda^\phi(\phi_C^{new,predicted}-\phi_C^{old})$ 来实现方程组系统迭代求解过程中的欠松弛处理。将该式修改为
$\phi_C^{}=\phi_C^{*}+\lambda^\phi(\phi_C^{new,iteration}-\phi_C^{*})$
其中 $\phi_C^{*}$ 为前次迭代所得值 $\phi_C$ 。在Patankar松弛方法中， $\phi_C^{new,iteration}$ 用其从 $a_C\phi_C+\displaystyle\sum_{F\sim NB(C)}a_F\phi_F=b_C$ 推出的等效表达式代入，得
$\phi_C^{}=\phi_C^{*}+\lambda^\phi\left(\frac{-\displaystyle\sum_{F\sim NB(C)}a_F\phi_F+b_C}{a_C}-\phi_C^{*}\right)$
重新整理后，有
$\frac{a_C}{\lambda^\phi}\phi_C+\sum_{F\sim NB(C)}a_F\phi_F=b_C+\frac{1-\lambda^\phi}{\lambda^\phi}a_C\phi_C^{*}$
上式中，松弛因子 $\lambda^\phi$ 修正了对角系数和右端项，而并没有修改方程的数学形式。因为 $\lambda^\phi<1$ ，欠松弛增加了代数方程组系统的对角占优特性，从而提高了迭代线性求解器的稳定性。这与显式方法相比是一个重大优势。

然而值得注意的是，隐式松弛施加的关系是与对角线系数成正比的，这样，对角线系数较大时松弛也较大，这就使得较小体积的单元松弛较大，影响较大，这将在下一节进行讲解。

2.3.2 E-Factor松弛

E-Factor方法是Patankar方法的重组，其通过把 $a_C\phi_C+\displaystyle\sum_{F\sim NB(C)}a_F\phi_F=b_C$ 重写为如下形式
$a_C\phi_C=b_C-\displaystyle\sum_{F\sim NB(C)}a_F\phi_F$
然后进行欠松弛处理，将上式的右端项转化为
$a_C\phi_C=\lambda^\phi\left(b_C-\displaystyle\sum_{F\sim NB(C)}a_F\phi_F\right)+ (1-\lambda^\phi)a_C\phi_C^*$
再将欠松弛因子替换为 $\dfrac{E^\phi}{1+E^\phi}$ ，得
$a_C\phi_C=\dfrac{E^\phi}{1+E^\phi}\left(b_C-\displaystyle\sum_{F\sim NB(C)}a_F\phi_F\right)+ \left(1-\dfrac{E^\phi}{1+E^\phi}\right)a_C\phi_C^*$
整理为
$a_C\left(1+\frac{1}{E^\phi}\right)\phi_C+\sum_{F\sim NB(C)}a_F\phi_F= b_C+\frac{1}{E^\phi}a_C\phi_C^*$
可以把上式中的欠松弛效应理解为某些人工瞬态时间尺度项，用于推进 $\phi_C$ 在每个迭代步的求解。时间步 $\Delta t$ 与特征时间间隔 $\Delta t^*$ 成正比
$\Delta t=E^\phi \Delta t^*$
其中
$\Delta t^*=\frac{\rho_C V_C}{a_C}$
上式中， $\rho_C$ 为体积为 $V_C$ 的单元 $C$ 的流体密度。特征时间间隔，是与 $\phi_C$ 的对流和扩散穿过一个单元的所需时间相关的。所以E-factor等效于单元CFL数。

从上式中清晰可见，E-Factor松弛的推进时间步是和单元体积相关的，小单元上推进得比粗糙单元上要慢。这对于稳态解的收敛速度是不利的，因为在边界附近经常会用到高拉伸率的单元且体积很小（那么这么小的时间步算到猴年马元才能收敛啊？！），这样就需要在计算域中开辟一块临界区域，和计算域的其余部分相比，该部分区域要用非常小的时间步长。这也是Patankar松弛方法的一个特性。

$E^\phi$ 和 $\lambda^\phi$ 的关系为
$E^\phi=\frac{1}{1-\lambda^\phi}$
一般来说， $E^\phi$ 选值范围是4-10，对应的 $\lambda^\phi$ 为0.75-0.9。

例2
在下图中示意性地展示了边界网格，单元 $A$ 和单元 $D$ 代表壁面边界附近的拉伸网格，其体积比为 $V_{C_A}/V_{C_D}\approx0.1$ ，这些网格的对角系数值通常是由扩散系数主导的，在该例中 $a_{C_A}/a_{C_D}\approx2$ ，因为 $A$ 有个边界面。

计算单元 $A$ 和单元 $D$ 的相对伪瞬态时间步长，欠松弛因子选为0.8。
在这里插入图片描述

解
求解等效 $E$ 因子
$E=\frac{1}{1-\lambda}=\frac{1}{1-0.8}=5$
这样，对每个单元的伪时间步为
$\Delta t=E^\phi \Delta t^*=5\frac{\rho_C V_C}{a_C}$
因此，单元 $A$ 和单元 $D$ 的相对伪瞬态时间步长为
$\frac{\Delta t_A}{\Delta t_B}= \left(\frac{V_{C_A}}{a_{C_A}}\right)/\left(\frac{V_{C_D}}{a_{C_D}}\right)= \left(\frac{V_{C_A}}{V_{C_D}}\right)/\left(\frac{a_{C_D}}{a_{C_A}}\right)\approx 0.1\left(\frac{1}{2}\right)=0.05$
表明单元 $D$ 的伪时间步长是单元 $A$ 的将近20倍。

2.3.3 伪瞬态松弛

伪瞬态松弛方法从Euler一阶瞬态隐式方法修改而来，用前次迭代值替换了旧时刻值。和Euler方法中一样，通过将伪瞬态项， $a_C^\circ\phi_C$ ，添加到对角系数中，以及将伪旧时间步项， $a_C^\circ\phi_C^*$ ，添加到右端项中，增加了代数方程的对角占优特性。修改后的方程变为
$(a_C+a_C^\circ)\phi_C+\sum_{F\sim NB(C)}a_F\phi_F=b_C+a_C^\circ\phi_C^*$
其中系数 $a_C^\circ$ 由下式计算
$a_C^\circ=\frac{\rho_C V_C}{\Delta t}$
其本质上等同于瞬态项的一阶隐式Euler离散所得到的瞬态系数， $\rho_C$ 是密度， $V_C$ 是单元体积，而 $\Delta t$ 是用户定义的伪时间步长。若 $\Delta t$ 较大，则添加项可忽略不计，欠松弛效应可忽略不计。若 $\Delta t$ 值非常小，则 $a_C^\circ$ 的值较大，且比其他项占优，从而使解有强烈的欠松弛特性，导致 $\phi_C$ 值的变化非常小（即， $\phi_C\approx\phi_C^*$ ）。

除了允许解的推进在整个计算域上是连续进行的之外，伪瞬态方法还能确保给对角系数添加的是非零贡献，即便在最极端的对角线系数是0的情况下亦如是。

并没有给定最佳欠松弛因子的通有准则，因为适用于一种情形的值未必适用于另一种情形。此外，不同的方程可能需要赋给不同的欠松弛因子。再者，在整个计算域上也未必要使用同一个松弛因子。更有甚者，欠松弛值也可能在迭代过程中是变化的。例如第15和16章将要讲到的SIMPLE算法，Raithby和Schneider推出了速度场和压力场的欠松弛因子的最佳关系，下一章将详细阐述。

3 方程的残差形式

把离散代数方程写成其“直接”或“标准”形式
$a_C\phi_C+\displaystyle\sum_{F\sim NB(C)}a_F\phi_F=b_C$
这也是OpenFOAM中使用的形式。

上式也可写成“修正”或“残差”形式，将其整理成求解修正值的形式。这样一来，若 $\phi_C^*$ 和 $\phi_C'$ 是 $\phi_C$ 的前次迭代值和所需满足上式的修正值，则解为
$\phi_C=\phi_C^*+\phi_C'$
从而把 $a_C\phi_C+\displaystyle\sum_{F\sim NB(C)}a_F\phi_F=b_C$ 改写为
$a_C(\phi_C^*+\phi_C')+\displaystyle\sum_{F\sim NB(C)}a_F(\phi_F^*+\phi_F')=b_C$
或
$a_C\phi_C'+\sum_{F\sim NB(C)}a_F\phi_F'=b_C-\left( a_C\phi_C^*+\sum_{F\sim NB(C)}a_F\phi_F^* \right)$
注意，上式右端项代表着 $\phi_C^*$ 场的方程的残差，把 $C$ 单元的残差用 $Res_C^\phi$ 表示，上式变为
$a_C\phi_C'+\sum_{F\sim NB(C)}a_F\phi_F'=Res_C^\phi$
注意，如果是精确场的话， $Res_C^\phi$ 应为零。

尽管在数学上来说，上式和最初的离散代数方程是等效的。然而上式的数值优势在于，在该形式下，方程在求解过程中的数值误差相对较小，而标准形式下的数值误差相对是较大的，即， $\phi$ 值较大时，其微小的变化可能会导致较大的数值误差。

3.2 Patankar欠松弛的残差形式

隐式Patankar松弛方程的残差形式可通过将方程 $a_C\phi_C+\displaystyle\sum_{F\sim NB(C)}a_F\phi_F=b_C$ 重写为如下形式来获得
$a_C(\phi_C^*+\phi_C')=\lambda^\phi\left(b_C-\sum_{F\sim NB(C)}a_F(\phi_F^*+\phi_F')\right)+ (1-\lambda^\phi)a_C\phi_C^*$
可简化为
$a_C\phi_C'+\lambda^\phi\sum_{F\sim NB(C)}a_F\phi_F'=\lambda^\phi\left[ b_C-\left( a_C\phi_C^* +\sum_{F\sim NB(C)}a_F\phi_F^* \right) \right]$
注意上式的右端项代表最初方程的残差，将上式重写为
$\frac{a_C}{\lambda^\phi}\phi_C'+\sum_{F\sim NB(C)}a_F\phi_F'=Res_C^\phi$
表明欠松弛处理对方程残差形式的改变只需要更改对角系数即可（不需要跟标准形式方程那样去改变右端项了）。

4 残差和解的收敛

在任何迭代求解过程中，一个非常重要的事情便是，啥时候解可被认为是足够好了，或者何时误差可被认为是在某个可容忍的范围以下，或者何时守恒方程的精度被认为是满足了。回答上面问题的工具是任何CFD代码中都非常重要的部分。这其实就是说，如何评估求解过程中流场的收敛程度，而事先并不知道最终解。基于此，数年来人们提出了许多指标，从最简单的监测单一位置点在迭代过程中的表现，到监测积分值比如阻力系数、整体质量流量、壁面剪切应力，等，或更普遍意义上的监测某些类型方程的残差。这些方法所面临的挑战是它们必须适用于众多不同的流动参数、几何构型和边界条件。

4.1 残差

在求解离散代数方程组 $a_C\phi_C+\displaystyle\sum_{F\sim NB(C)}a_F\phi_F=b_C$ 的过程中，可用下述形式的单元残差来衡量平衡方程的误差
$Res_C^\phi=b_C-\left(a_C\phi_C+\sum_{F\sim NB(C)}a_F\phi_F\right)$
显然，当求得了满足方程的解后， $Res_C^\phi$ 将会是零。使用 $Res_C^\phi$ ，可推导出适用于整个计算域的一系列残差指标，如下。

4.2 绝对残差

上式中定义的残差，其符号可能是正的，也可能是负的。因为符号实际上无关紧要，所以把 $Res_C^\phi$ 的绝对值，定义为 $R_C^\phi$ ，用于判断解是否收敛。如果 $R_C^\phi$ 随着迭代是减小的，那么解是收敛的，否则解是发散的。在点 $C$ 处的 $R_C^\phi$ 定义为
$R_C^\phi=\left| b_C-\left(a_C\phi_C+\sum_{F\sim NB(C)}a_F\phi_F\right) \right|$

4.3 最大残差

当在整个区域中的绝对残差的最大值低于某个非常小的量时，可以认为解是收敛的。绝对残差的最大值定义为
$R_{C,max}^\phi=\max_{all~cells}\left| b_C-\left(a_C\phi_C+\sum_{F\sim NB(C)}a_F\phi_F\right) \right|=\max_{all~cells}R_C^\phi$
非常小的量用 $\epsilon$ 表示，则收敛标准为
$R_{C,max}^\phi \le \epsilon \Rightarrow \text{solution has converged}$

4.4 均方根残差

另一个用于收敛判断的指标，是把所有节点的绝对残差 $R_C^\phi$ 平方后加和，再取平均后开方，即可得到所谓的均方根残差 $R_{C,rms}^\phi$ ，其数学形式为
$\begin{aligned} R_{C,rms}^\phi &= \sqrt{\frac{\displaystyle\sum_{C\sim all~elements}\left[ b_C-\left(a_C\phi_C+\sum_{F\sim NB(C)}a_F\phi_F\right) \right]^2} {number~of~elements}} \\ &= \sqrt{\frac{\displaystyle\sum_{C\sim all~elements}\left(R_C^\phi\right)^2} {number~of~elements}} \end{aligned}$
将收敛指标写成
$R_{C,rms}^\phi \le \epsilon \Rightarrow \text{solution has converged}$

4.5 残差的归一化

绝对残差是变量 $\phi$ 的函数，因此不同的 $\phi$ 值会产生不同的 $R_C^\phi$ ，这使其非常难判断解到底收敛与否（举个例子，如果采用有量纲的变量表示残差，如果精确解是1e10，那么残差1e4可能就OK了，如果精确解是1e-6，那么残差得是1e-12才能OK，即，绝对残差是跟变量值的大小相关的）。此时，较好的处理方法是，将不同的残差除以它们当中的最大通量，来做缩放处理（相当于向量的归一化（单位化，无量纲化）处理）。因为 $a_C$ 代表单元上通量的加和，那么残差的缩放，可以把它除以整个区域中的最大 $a_C\phi_C$ ，以获得残差相对于局部值 $\phi$ 缩放处理后的相对值，数学形式为
$R_{C,scaled}^\phi=\frac{\left| b_C-\left(a_C\phi_C+\displaystyle\sum_{F\sim NB(C)}a_F\phi_F\right) \right|} {\displaystyle\max_{all~cells}|a_C\phi_C|}$
解达到收敛的标准为，缩放绝对残差的最大值小于一个很小的量，即
$\max_{all~cells}\left(R_{C,scaled}^\phi\right)\le\epsilon \Rightarrow \text{solution has converged}$
一般来说，缩放残差的收敛阈值 $\epsilon$ 量级应为 $10^{-3}\sim 10^{-5}$ 或者更小，以保证收敛。（我一般喜欢用 $10^{-6}$ ）。

除了使用绝对或缩放残差外，在解收敛过程中，实时监测如前所述的积分量（阻力系数、整体质量流量、壁面剪切应力等）也是蛮有借鉴意义的。在认定解是收敛前一定要确认是否为正确的合适的收敛，因为没有收敛的解经常是毫无物理意义的，还容易误导人们的理解（即，把数值误差产生的乱七八糟的混账玩意儿误认为是发现了新大陆而欣喜若狂……）。

5 代码讲解

5.1 uFVM

5.1.1 源项线性化

在uFVM中，对源项的线性化和组装是在函数cfdAssembleSourceTerm中进行的。函数需要用户指定线性化集合，即，提供源的常数部分，Sb，和线性化部分，Sc。这些项将添加进FLUXV和FLUXC中。

注意在uFVM中，方程是用残差形式求解的，这也解释了为何将整体源添加到FLUXTE，而非只有其常数项。

theEquationField = cfdGetMeshField(theEquationName);
phi = theEquationField.phi(iElements);
%
Sb = cfdComputeFormulaAtLocale(theTerm.Sb,'Interior Elements')';
Sc = cfdComputeFormulaAtLocale(theTerm.Sc,'Interior Elements')';
%
volume = [theMesh.elements.volume];
%
% Assemble Source Term
%
pos = zeros(1,size(phi));
pos(Sc<0) = 1;
theFluxes.FLUXCE = -pos .* Sc .* volume;
theFluxes.FLUXTE = -(Sb + Sc .*phi) .* volume;

5.1.2 欠松弛

uFVM中使用的是Patankar的隐式欠松弛方法。因为方程是用残差形式来求解的，所以其欠松弛只需要修改对角系数，即可。

function cfdApplyURF(theEquationName)
%===================================================
% written by the CFD Group @ AUB, Fall 2006
%===================================================
theEquation = cfdGetModel(theEquationName);
urf = theEquation.urf;
theCoefficients = cfdGetCoefficients;
theCoefficients.ac = theCoefficients.ac/urf;
cfdSetCoefficients(theCoefficients);

5.2 OpenFOAM

to be continued…

梅冠华

关注

5
点赞
踩
16

收藏

觉得还不错? 一键收藏
0
评论
FVM in CFD 学习笔记_第14章_源项离散，松弛，其它细节

本章咱们讲讲那些可以对解的特性产生“较大”影响的“细小”数值细节。首先，考察源项的处理，在通常状况下它是依赖于解的，源项将线性化处理成因变量的形式，并分成两部分，一部分显式处理，另一部分隐式处理。紧接着，讨论对代数方程组做欠松弛处理的显式和隐式技术，并展现一些隐式欠松弛方法，有著名的Patankar隐式欠松弛方法、van Doormaal和Raithby的E-因子方法、Mallinson和de Vahl Davis的伪瞬态方法。然后，介绍离散代数方程组的残差形式。最后是评估解收敛与否的收敛指标。
复制链接

扫一扫