本文简单聊聊CFD计算中的残差问题。
群里经常有道友讨论计算收敛的问题。问得最多的问题大致可以归纳成下面的几类：
残差曲线平了，但没有降低到残差标准，这样能算收敛么？
残差曲线在周期性振荡，这样是收敛了么？
残差曲线未降低到指定标准，计算结果是否可以使用？
要回答这些问题，还得回到问题的起点。
1 残差是什么？
首先要弄明白残差是什么。在CFD离散的过程中，需要将描述物理世界的偏微分方程转化为可供计算机求解的代数方程组。
求解代数方程组的方法主要有两种：直接法与迭代法。
1.1 直接法
直接法指的是计算机通过算法一次性将方程组求解出来，在求解的过程中不会有任何中间解。因此直接法不会有残差的概念。典型的直接法如高斯消去法、矩阵求逆法、LU分解法等，这些在线性代数中应该都有学过。举个简单的例子。 如下面的方程组：
写成矩阵形式：
采用矩阵求逆方法可以解出：
看，一下子就能把结果算出来了。采用直接法求解代数方程组，只有能求解和不能求解，没有求解精确不精确一说，只要能求解的，得到都是精确解。
然而CFD求解的问题规模可能会很大，大到远远超出计算机内存所能存储的极限，这时候直接法就比较难使用了，此时常采用迭代法。
1.2 迭代法
还是上面的方程组，可以使用最简单的迭代法（如Jacobi迭代），将方程组改造成下面的迭代式：
迭代法需要初始值，否则迭代没有办法进行下去。如上面的方程组中，为了演示，给定一个不怎么靠谱的初始值：
这个例子比较简单，下面是迭代20步计算得到的各变量的值。
迭代
x1
x2
x3
x4
x5
0
0.0
0.0
0.0
0.0
0.0
1
55.00000
6.66667
6.66667
6.66667
10.00000
2
56.66667
27.22222
11.11111
12.22222
13.33333
3
61.80556
29.25926
19.81481
14.81481
16.11111
4
62.31481
33.87346
21.35802
18.64198
17.40741
5
63.46836
34.55761
24.17181
19.58848
19.32099
6
63.63940
35.88006
24.71536
21.16427
19.79424
7
63.97001
36.11826
25.68144
21.50320
20.58213
8
64.02956
36.55049
25.87382
22.08786
20.75160
9
64.13762
36.63446
26.21278
22.20847
21.04393
10
64.15862
36.78347
26.28098
22.41890
21.10424
11
64.19587
36.81320
26.40079
22.46174
21.20945
12
64.20330
36.86555
26.42498
22.53675
21.23087
13
64.21639
36.87609
26.46743
22.55195
21.26837
14
64.21902
36.89461
26.47601
22.57860
21.27597
15
64.22365
36.89835
26.49107
22.58400
21.28930
16
64.22459
36.90491
26.49411
22.59346
21.29200
17
64.22623
36.90623
26.49945
22.59537
21.29673
18
64.22656
36.90856
26.50053
22.59873
21.29769
19
64.22714
36.90903
26.50243
22.59941
21.29936
20
64.22726
36.90986
26.50281
22.60060
21.29970
可以看出，随着迭代的进行，计算结果越来越趋近于式（3）中的真实解。
1.2.1 残差评估方法1
将每一步迭代计算得到的x的值与正确值[64.227，36.911，26.504，22.602，21.301]相减，可以得到每次迭代后各变量的残差，如下表所示。
Res(x1)
Res(x2)
Res(x3)
Res(x4)
Res(x5)
0
-64.22700
-36.91100
-26.50400
-22.602
-21.301
1
-9.22700
-30.24433
-19.83733
-15.93533
-11.301
2
-7.56033
-9.68878
-15.39289
-10.37978
-7.96767
3
-2.42144
-7.65174
-6.68919
-7.78719
-5.18989
4
-1.91219
-3.03754
-5.14598
-3.96002
-3.89359
5
-0.75864
-2.35339
-2.33219
-3.01352
-1.98001
6
-0.58760
-1.03094
-1.78864
-1.43773
-1.50676
7
-0.25699
-0.79274
-0.82256
-1.09880
-0.71887
8
-0.19744
-0.36051
-0.63018
-0.51414
-0.54940
9
-0.08938
-0.27654
-0.29122
-0.39353
-0.25707
10
-0.06838
-0.12753
-0.22302
-0.18310
-0.19676
11
-0.03113
-0.09780
-0.10321
-0.14026
-0.09155
12
-0.02370
-0.04545
-0.07902
-0.06525
-0.07013
13
-0.01061
-0.03491
-0.03657
-0.05005
-0.03263
14
-0.00798
-0.01639
-0.02799
-0.02340
-0.02503
15
-0.00335
-0.01265
-0.01293
-0.01800
-0.01170
16
-0.00241
-0.00609
-0.00989
-0.00854
-0.00900
17
-0.00077
-0.00477
-0.00455
-0.00663
-0.00427
18
-0.00044
-0.00244
-0.00347
-0.00327
-0.00331
19
0.00014
-0.00197
-0.00157
-0.00259
-0.00164
20
0.00026
-0.00114
-0.00119
-0.00140
-0.00130
为方便以对数坐标轴绘制残差曲线，这里取绝对值。
可以看出，随着迭代的进行，残差值越来越小，迭代20次后，每个偏离与正确值的偏差接近0.001。
1.2.2 残差评估方法2
然而在实际的计算过程中，我们并不知道变量的正确值，这时候就没有办法直接与正确值进行比较来获取每次迭代计算的残差。此时可以使用相邻两次迭代计算的差值作为当前迭代步的残差，即：
若采用这种方式，可以得到残差如下表所示。
迭代
Res(x1)
Res(x2)
Res(x3)
Res(x4)
Res(x5)
0
0
0
0
0
0
1
55
6.666667
6.666667
6.666667
10
2
1.666667
20.55556
4.444444
5.555556
3.333333
3
5.138889
2.037037
8.703704
2.592593
2.777778
4
0.509259
4.614198
1.54321
3.82716
1.296296
5
1.153549
0.684156
2.813786
0.946502
1.91358
6
0.171039
1.322445
0.543553
1.575789
0.473251
7
0.330611
0.238197
0.966078
0.338935
0.787894
8
0.059549
0.43223
0.192377
0.584657
0.169467
9
0.108057
0.083976
0.338962
0.120615
0.292329
10
0.020994
0.149007
0.068197
0.21043
0.060307
11
0.037252
0.02973
0.119812
0.042835
0.105215
12
0.007433
0.052355
0.024188
0.075009
0.021417
13
0.013089
0.01054
0.042455
0.015202
0.037505
14
0.002635
0.018514
0.008581
0.026653
0.007601
15
0.004629
0.003739
0.015056
0.005394
0.013327
16
0.000935
0.006561
0.003044
0.009461
0.002697
17
0.00164
0.001326
0.005341
0.001914
0.00473
18
0.000332
0.002327
0.00108
0.003357
0.000957
19
0.000582
0.000471
0.001895
0.000679
0.001679
20
0.000118
0.000825
0.000383
0.001191
0.000339
随着计算的进行，残差也是逐渐降低且趋于零。
可以看到，采用相邻迭代步计算得到的残差要比与正确值比较得到的残差值更小。
需要注意的是，采用此中方式反映的是本次迭代值相对于相邻迭代步上的变量值的变化量，并不能表征本次迭代值与正确值的接近程度。
那如果要评价本次迭代与正确值的差异的话，该怎么办呢？看方法3。
1.2.3 残差评估方法3
当正确值未知时，需要转换一下。
令：
其中为当前迭代步上的的值，为与正确值的偏差。则方程(2)可表示为：
转换一下形式可以得到：
可将上式的右侧视为残差，即残差可表示为：
式(9)可以表征当前迭代值的接近程度，若为正确值，则残差为零。越接近正确值，残差越小。 采用此方法计算得到的残差值如下表所示。
迭代次数
Res(x1)
Res(x2)
Res(x3)
Res(x4)
Res(x5)
0
1100
100
100
100
100
1
33.33333
308.33333
66.66667
83.33333
33.33333
2
102.77778
30.55556
130.55556
38.88889
27.77778
3
10.18519
69.21296
23.14815
57.40741
12.96296
4
23.07099
10.26235
42.20679
14.19753
19.13580
5
3.42078
19.83668
8.15329
23.63683
4.73251
6
6.61223
3.57296
14.49117
5.08402
7.87894
7
1.19099
6.48345
2.88566
8.76986
1.69467
8
2.16115
1.25963
5.08444
1.80922
2.92329
9
0.41988
2.23510
1.02295
3.15646
0.60307
10
0.74503
0.44595
1.79718
0.64252
1.05215
11
0.14865
0.78532
0.36282
1.12514
0.21417
12
0.26177
0.15810
0.63682
0.22803
0.37505
13
0.05270
0.27772
0.12871
0.39980
0.07601
14
0.09257
0.05608
0.22584
0.08091
0.13327
15
0.01869
0.09842
0.04566
0.14191
0.02697
16
0.03281
0.01989
0.08011
0.02871
0.04730
17
0.00663
0.03491
0.01620
0.05036
0.00957
18
0.01164
0.00706
0.02842
0.01018
0.01679
19
0.00235
0.01238
0.00575
0.01787
0.00339
20
0.00413
0.00250
0.01008
0.00361
0.00596
图形显示为：
可以看到采用此方法计算的残差要比方法2保守得多，而且残差存在振荡，不过总体趋势是残差在减小。需要注意，残差是否振荡或持续减小，取决于原方程组的特性以及所采用的迭代方法。
那么残差到底是什么？可以这么认为：迭代残差是用于评价迭代解与方程组正确解之间差异的数值，残差越小，意味着迭代解越接近方程组的正确解。
2 前面的问题
认识到残差的含义之后，再来回答前面的问题。
残差曲线平了，但没有降低到残差标准，这样能算收敛么？
残差曲线在周期性振荡，这样是收敛了么？
残差不减小意味着当前迭代值与方程组的正确解之间存在偏差。残差曲线为一条水平线，表示每次迭代值都存在相同的偏差，如果偏差值比较大，显然结果不能算作收敛。残差曲线周期性振荡，可能的原因很多，这可能是选用了不合适的迭代算法，也有可能是采用稳态方法求解了瞬态问题导致。只要残差没有降低到我们认可的标准，就不可以说计算已经收敛。
瞬态计算与稳态计算的要求一样，需要每个时间步内迭代残差下降到一定程度。
再来看另外一个问题。
残差曲线未降低到指定标准，计算结果是否可以使用？
需要注意的是：残差仅用于评价迭代法求解代数方程组时中间解与正确解的逼近程度，就算残差为零，也只能说我们得到了方程组的正确解而已。然而，得到了方程组的正确解，也未必就能说得到了正确的物理值，因为从物理现象到代数方程组，中间还经历了大量简化和近似处理。残差未降到零，表示迭代计算结果与方程组的正确值存在偏差。但这个偏差到底是哪个网格位置的迭代值残差呢？
想象一下加入有100万网格，那么对于某一个变量的离散方程组求解后就会有100万个迭代值与残差值，这时需要在这100万个迭代残差值中挑选一个最具有代表性的值（如取最大值、均方根等），加入以整个区域内所有网格上的最大残差值作为班次迭代计算的变量残差值，这个残差值有可能是在感兴趣区域的网格上，当然也有可能是在非重要区域的网格上。如果问题出在感兴趣区域，自然大残差不可接受，但如果不感兴趣区域由于网格质量或网格尺寸造成的大残差，此时是否可以接受？有时候是可以接受的。然而计算过程中显然不太可能输出每个网格上的残差值（当然可以做，但是太浪费计算资源，没人会这么做），所以这时候经常用监测物理量的方式对计算结果进行辅助判断。
通过在感兴趣区域设置若干监测点，查看迭代计算过程中监测点位置感兴趣物理量的变化趋势。若随着迭代进行，物理量几乎没有变化，此时可以认为计算可用。
另外，计算结果是否可用，还需要结合其他数据（如理论分析及实验测试）综合判断，不能仅仅只看残差或监测数据。