【小白的CFD之旅】26 何为收敛
小白最近对流体计算的收敛产生了困惑。以前在学习高等数学的时候，小白接触过了级数的收敛，由于当时贪玩，并未将其放在心上，因此大学结束了小白也只是记住有这么一个名词罢了。现如今在利用CFD的过程中，小白又一次碰到了“收敛”这一名词。小白找了很多的资料，然而资料中关于收敛的介绍，无一不是一大堆的数学公式，小白觉得头很疼。
“出来混，总是要还的。”小白的心情很复杂。
“流体计算为什么要收敛？收敛是什么意思？不收敛又有什么后果？如何判断是会否收敛？如果不收敛该采取何种措施使其收敛？”小白带着满腔的疑惑，找到了小牛师兄。
“师兄，最近在看资料的时候，老是碰到收敛这个名词，我想问一下，什么叫收敛呢？我最近看了很多的资料，都是用一堆数学公式来表达，看了很久概念也很模糊。”小白问道。
“也好，既然你注意到了收敛，今天我们就来聊一聊这个话题。”小牛师兄说。
计算机解决工程问题
“在谈收敛之前，我们还是先来谈一谈数值方法。关于数值计算，我想你也了解了不少。你说计算机解决这些工程问题到底是如何实现的呢？”小牛师兄问道。
“我看了一些资料。利用计算机来解决工程问题，其核心是利用计算机求解物理模型，或者说是利用计算机求解数学方程。”小白说。
“基本上是这样。我们知道描述现实世界物理现象的数学模型大部分都是微分方程，利用计算机并不能直接求解微分方程，计算机通常需要利用数学方法将这些微分方程转化为代数方程，通过求解代数方程获取原微分方程的解。那么问题来了，计算机是如何将微分方程转化为代数方程的呢？”小牛师兄说。
“转化的方法其实有很多种，比如我们常说的有限差分法、有限元法、有限体积法等都是用来干这种工作的。为了实现这种软件，引入了网格。采用网格的目的是将连续的空间和时间转化为离散的空间与时间，从而在每一个网格单元或网格节点上获得代数方程。所有网格上的代数方程集合在一起，就构成了整个计算域上的代数方程组，求解方程组就可以得到每一个网格单元或网格节点上的物理量。”
“当然，求解方程组得到的物理量还是离散的，为了得到节点之间或单元之间的物理量分布，采用了一种称之为插值的数学方法。利用插值可以得到近似连续的物理量分布。当计算网格或计算时间间隔足够小的时候，通过插值可以得到相当精确的计算结果。”小牛师兄说。
“因此我们可以说，利用计算机求解工程问题，实际上归根结底是利用计算机求解工程问题所涉及的物理模型转化而成的代数方程。那么问题来了，利用计算机求解代数方程是怎么实现的呢？”小牛师兄说。
“我记得线
性代数里面有讲到，可以采用矩阵求逆、高斯消去法、三角分解法、克拉默法则等。”小白说。
“嗯，这些方法的确可以用于计算机求解线性方程组，这些方法可以统称为直接求解法。然而在工程应用上是存在很大问题的。我们知道，为了得到比较精确的计算结果，网格数量常常非常多，直接后果就是最终的代数方程组包含的方程数量非常非常多。打个比方，如果网格数量是100万个，那么最终求解的代数方程组就有100万个方程，在利用直接法求解的过程中，需要采用100万维的矩阵进行存储，如此庞大的矩阵，消耗的内存是非常巨大的，而且其求逆计算或三角分解计算都是非常困难的，计算效率是极其低下的。因此，工程上求解这类数量庞大的方程组通常采用迭代法。”小牛师兄说。
“迭代法其实也有很多种，比较常见的如雅克比迭代、高斯-赛德尔迭代、超松弛迭代等等。关于迭代法求解线性方程组的具体方法，可以找一本数值计算的数看看。收敛就是迭代法求解线性方程组中的一个概念。”小牛师兄最后说。
关于迭代残差
为了更好的说明收敛的问题，这里先举个例子来说明什么叫残差。
比如说要求解下面的线性方程组：
可将(1)式改写为迭代计算的形式：
给定一个初值，如x(0)=(0,0,0)T，代入到式(2)中，可以得到x(1)=(2.5,3,3)T
反复迭代，
可以得到一系列x1,x2及x3的值。
| 迭代次数 | x1 | x2 | x3 |
|---|---|---|---|
| 0 | 0 | 0 | 0 |
| 1 | 2.5 | 3 | 3 |
| 2 | 2.875 | 2.363636364 | 1 |
| 3 | 3.136363636 | 2.045454545 | 0.971590909 |
| 4 | 3.024147727 | 1.947830579 | 0.920454545 |
| 5 | 3.000322831 | 1.983987603 | 1.000968492 |
| 6 | 2.993753228 | 1.999970652 | 1.003841684 |
| 7 | 2.999028573 | 2.002620797 | 1.003130723 |
| 8 | 3.000200118 | 2.000637857 | 0.999830514 |
| 9 | 3.000281568 | 1.999911822 | 0.999740477 |
| 10 | 3.000031814 | 1.999874019 | 0.999881261 |
| 11 | 2.999982442 | 1.999977637 | 1.000015588 |
| 12 | 2.999987717 | 2.000007802 | 1.00001437 |
| 13 | 2.999999333 | 2.000005773 | 1.000004191 |
| 14 | 3.000001117 | 2.000000623 | 0.99999889 |
| 15 | 3.000000511 | 1.999999493 | 0.999999286 |
| 16 | 2.999999988 | 1.999999749 | 0.999999871 |
| 17 | 2.999999938 | 1.999999992 | 1.000000068 |
| 18 | 2.99999998 | 2.000000029 | 1.000000033 |
| 19 | 3.000000003 | 2.00000001 | 1.000000003 |
| 20 | 3.000000003 | 1.999999999 | 0.999999996 |
可以看出，随着迭代计算的进行，计算结果越来越接近解析解[3,2,1]。
我们可以将某一物理量两次迭代计算值的差值称之为该物理量的迭代残差。比如上面的迭代计算x1的第一次迭代计算残差为2.5，第二次迭代计算后残差变为0.375。
当然上面提到的残差为绝对残差，在实际计算过程中，有时也常常采用相对残差，即物理量随迭代进行的变化量。可以看到，当残差非常小时，可以认为不用再往下计算了，或者说再往下算已经没有太大意义了。
CFD计算中的残差
在CFD计算中，每一个网格上都会存储众多物理量，因此每一个网格上的任一个物理量在计算迭代过程中都会存在一个残差，这意味着在一次迭代过程中，同一物理量在不同的计算网格上有不同的计算残差，而实际上我们在进行CFD计算时，每一个迭代步只对应着一个残差值。
CFD中残差分为几种：
最大残差：在一次迭代中，取所有网格中的残差值的最大值作为本次计算的残差。
平均残差：在一次迭代中，计算所有网格中的算术平均值作为本次迭代的残差
均方根残差：在一次计算，计算所有网格中残差值的均方根作为本次迭代计算的残差
在CFD计算中，常常采用均方根残差（RMS）作为残差值。
CFD中的收敛
所谓迭代收敛，简单来说，就是在迭代计算过程中，物理量趋于某一值的情况。CFD计算中判断收敛通常有三种方法。
收敛判断规则之一：残差达到某一设定标准时可以认为迭代计算达到收敛。
这条规则最简单，在实际工程应用中也最常用。通常设定某一标准，当迭代计算过程中残差值低于此标准时则认为计算收敛。这也是几乎所有CFD软件用于判断收敛的基本方法。
然而此规则在实际工程应用中常常无法使用，有些复杂的问题，无论你怎么计算，其残差也不会下降，甚至有时候残差会出现周期性震荡。
残差稳定在某一位置不下降的原因有很多，常见的原因包括：
- 计算区域中存在低质量的网格。低质量的网格会造成计算残差增大及残差震荡
- 边界条件设置有误。错误的边界条件或边界类型搭配都会导致计算残差震荡。
- 利用稳态求解器计算瞬态问题也会造成残差的震荡。
在残差无法达到设定标准的时候该如何判断收敛呢？需要注意的是，此时的收敛并非数学意义上的收敛了，而只是意味着我们可以停止计算。
收敛判断规则之二：进出口物理量通量达到平衡
最常用的是判断进出口质量是否相等，这实际上是判断连续性条件是否达到满足的。实际上还有很多，如计算域中包含化学反应时，判断进出口组分是否守恒；如计算域中包含多相流时，判断进出口各相质量是否守恒等。此规则是一种非常弱的规则，实际上只是收敛的一个必要条件而已。但是在第一条判断规则无法达到时，也常常采用此规则来判断。
收敛判断规则三：计算域中的物理量随迭代进行不再发生变化
这条规则在实际应用中也很常用，甚至比第二条规则更常用。在实际工程中，经常监测某些敏感位置的物理量，当随着迭代进行，监测的量不再发生变化时，基本可以认为计算达到收敛。
注意后两种方法只是在残差无法达到标准时才采用的判断方法，并不意味着计算就收敛了。
稳态及瞬态计算的收敛
在CFD迭代计算过程中，稳态残差曲线与瞬态残差曲线的形态有很大的不同。
如下图所示为稳态计算残差，每一个物理量都有一条残差曲线，当残差曲线低于设定的标准时认为计算收敛，当所有的残差曲线低于设定的标准时，计算结束。
下图所示为典型的瞬态计算残差。可以看到其形状特征与稳态残差曲线有很大的不同。瞬态计算要求在每一个计算时间步内达到收敛。若不能达到收敛，则迭代次数达到所设定的内迭代次数后进入下一个时间步，重新开启迭代计算。因此书台计算残差呈现出下图所示的锯齿状。
注意：瞬态计算要求每一个时间步内均达到收敛。
