聊聊CFD计算中的残差

本文简单聊聊CFD计算中的残差问题。

群里经常有道友讨论计算收敛的问题。问得最多的问题大致可以归纳成下面的几类：

• 残差曲线平了，但没有降低到残差标准，这样能算收敛么？
• 残差曲线在周期性振荡，这样是收敛了么？
• 残差曲线未降低到指定标准，计算结果是否可以使用？

要回答这些问题，还得回到问题的起点。

1 残差是什么？

首先要弄明白残差是什么。在CFD离散的过程中，需要将描述物理世界的偏微分方程转化为可供计算机求解的代数方程组。

求解代数方程组的方法主要有两种：直接法与迭代法。

1.1 直接法

直接法指的是计算机通过算法一次性将方程组求解出来，在求解的过程中不会有任何中间解。因此直接法不会有残差的概念。典型的直接法如高斯消去法、矩阵求逆法、LU分解法等，这些在线性代数中应该都有学过。举个简单的例子。如下面的方程组：

写成矩阵形式：

采用矩阵求逆方法可以解出：

看，一下子就能把结果算出来了。采用直接法求解代数方程组，只有能求解和不能求解，没有求解精确不精确一说，只要能求解的，得到都是精确解。

然而CFD求解的问题规模可能会很大，大到远远超出计算机内存所能存储的极限，这时候直接法就比较难使用了，此时常采用迭代法。

1.2 迭代法

还是上面的方程组，可以使用最简单的迭代法（如Jacobi迭代），将方程组改造成下面的迭代式：

迭代法需要初始值，否则迭代没有办法进行下去。如上面的方程组中，为了演示，给定一个不怎么靠谱的初始值：

这个例子比较简单，下面是迭代20步计算得到的各变量的值。

可以看出，随着迭代的进行，计算结果越来越趋近于式（3）中的真实解。

1.2.1 残差评估方法1

将每一步迭代计算得到的x的值与正确值[64.227，36.911，26.504，22.602，21.301]相减，可以得到每次迭代后各变量的残差，如下表所示。

为方便以对数坐标轴绘制残差曲线，这里取绝对值。

 

可以看出，随着迭代的进行，残差值越来越小，迭代20次后，每个偏离与正确值的偏差接近0.001。

1.2.2 残差评估方法2

然而在实际的计算过程中，我们并不知道变量的正确值，这时候就没有办法直接与正确值进行比较来获取每次迭代计算的残差。此时可以使用相邻两次迭代计算的差值作为当前迭代步的残差，即：

若采用这种方式，可以得到残差如下表所示。

随着计算的进行，残差也是逐渐降低且趋于零。

 

可以看到，采用相邻迭代步计算得到的残差要比与正确值比较得到的残差值更小。

需要注意的是，采用此中方式反映的是本次迭代值相对于相邻迭代步上的变量值的变化量，并不能表征本次迭代值与正确值的接近程度。

那如果要评价本次迭代与正确值的差异的话，该怎么办呢？看方法3。

1.2.3 残差评估方法3

当正确值未知时，需要转换一下。

令：

其中    为当前迭代步上的    的值，    为    与正确值    的偏差。则方程(2)可表示为：

转换一下形式可以得到：

可将上式的右侧视为残差，即残差    可表示为：

式(9)可以表征当前迭代值    的接近程度，若    为正确值，则残差    为零。    越接近正确值，残差越小。采用此方法计算得到的残差值如下表所示。

图形显示为：

 

可以看到采用此方法计算的残差要比方法2保守得多，而且残差存在振荡，不过总体趋势是残差在减小。需要注意，残差是否振荡或持续减小，取决于原方程组的特性以及所采用的迭代方法。

那么残差到底是什么？可以这么认为：迭代残差是用于评价迭代解与方程组正确解之间差异的数值，残差越小，意味着迭代解越接近方程组的正确解。

2 前面的问题

认识到残差的含义之后，再来回答前面的问题。

• 残差曲线平了，但没有降低到残差标准，这样能算收敛么？
• 残差曲线在周期性振荡，这样是收敛了么？

残差不减小意味着当前迭代值与方程组的正确解之间存在偏差。残差曲线为一条水平线，表示每次迭代值都存在相同的偏差，如果偏差值比较大，显然结果不能算作收敛。残差曲线周期性振荡，可能的原因很多，这可能是选用了不合适的迭代算法，也有可能是采用稳态方法求解了瞬态问题导致。只要残差没有降低到我们认可的标准，就不可以说计算已经收敛。

瞬态计算与稳态计算的要求一样，需要每个时间步内迭代残差下降到一定程度。

再来看另外一个问题。

• 残差曲线未降低到指定标准，计算结果是否可以使用？

需要注意的是：残差仅用于评价迭代法求解代数方程组时中间解与正确解的逼近程度，就算残差为零，也只能说我们得到了方程组的正确解而已。然而，得到了方程组的正确解，也未必就能说得到了正确的物理值，因为从物理现象到代数方程组，中间还经历了大量简化和近似处理。残差未降到零，表示迭代计算结果与方程组的正确值存在偏差。但这个偏差到底是哪个网格位置的迭代值残差呢？

想象一下加入有100万网格，那么对于某一个变量的离散方程组求解后就会有100万个迭代值与残差值，这时需要在这100万个迭代残差值中挑选一个最具有代表性的值（如取最大值、均方根等），加入以整个区域内所有网格上的最大残差值作为班次迭代计算的变量残差值，这个残差值有可能是在感兴趣区域的网格上，当然也有可能是在非重要区域的网格上。如果问题出在感兴趣区域，自然大残差不可接受，但如果不感兴趣区域由于网格质量或网格尺寸造成的大残差，此时是否可以接受？有时候是可以接受的。然而计算过程中显然不太可能输出每个网格上的残差值（当然可以做，但是太浪费计算资源，没人会这么做），所以这时候经常用监测物理量的方式对计算结果进行辅助判断。

通过在感兴趣区域设置若干监测点，查看迭代计算过程中监测点位置感兴趣物理量的变化趋势。若随着迭代进行，物理量几乎没有变化，此时可以认为计算可用。

另外，计算结果是否可用，还需要结合其他数据（如理论分析及实验测试）综合判断，不能仅仅只看残差或监测数据。