CFD中的预处理方法

10月前浏览1137

新型冠状病毒肺炎疫情仍在持续，全国上下正在有序开展防控工作，疫情当前，减少外出活动、做好个人防护、远离人群是最好的预防手段。宅在家的这段时间也是充电的好时机，大家可以系统梳理、学习CFD，以下主要分享CFD中的预处理方法，欢迎留言交流~

为什么需要预处理？

现代CFD中，会经常听到「刚性」、「病态」以及「预处理」等术语？它们是什么意思呢？

「刚性」对应英文单词stiffness，「预处理」对应单词preconditioning，「病态」对应ill-conditioning。它们都与另外一个单词有关：条件数，即condition number，ϗ。

如果一个方程组的条件数很大，那么它就被认为是刚性的。如果一个代数方程组的系数矩阵的条件数很大，那么就可以说它是病态的（ill-conditioned）。条件数可以简单的认为是输入误差的放大倍数。求解一个方程组，如果对应的矩阵条件数很大，那么输入误差就会像放大镜一样放大。譬如其条件数是ϗ=10¹²，对于64位双精度计算过程，可以期望的输入误差一般不会小于10^-15，那么无论采用何种高精度的数值格式，多么密集的网格，结果的误差不会小于0.0001。对于大多数计算来说，这是不可忍受的。

一般的，如果一个矩阵的条件数非常大，那么它是病态的，或者刚性的。对于刚性或者病态系统，除了改变系统本身的条件数，没有其它特别奏效的数值方法。尽管本质相同，我们还是可以粗略的把刚性/病态问题分为两类。一类是微/积分方程组的刚性问题，另外一类则是针对一个特定的代数方程组。低速时的可压缩控制方程属于第一类，通常用的共轭梯度方法往往需要配合预处理才能有效的求解病态的代数方程组，这是第二类。

这里简单介绍第一类问题的产生根源和通用的处理方法。

一个例子

航空航天应用中，工程师面对的一般都是可压缩流动。然而其中有些情形，流动速度比较低，比如飞机起降时的绕机翼流场。还有些情形流动更加复杂，如直升机的旋翼附近的流场，流动是高、低速混合的。这种情况下，无论是不可压缩求解器，亦或是纯粹的可压缩流动求解器都不适用。

粗略的说，流动中存在两个典型的特征速度，一个是流动速度v，一个是声速c，后者压力脉动的传播速度。

上面所述流场，流动的速度和声速存在非常大的差异，二者的比值:

这就是可压缩流动控制方程刚性产生的物理根源。

控制可压缩流动的Euler方程组（以下简称Euler方程组）可以简写为如下的形式

这里假设流动是一维定常的，τ 是伪时间（pseudo-time）。其中Q = (ρ, ρu, ρE)是守恒变量， F = (ρ,ρuu+p, ρuE)是对应的通量，E是总能。将该方程改写为拟线性的形式：

其中是Jacobi矩阵，三个特征值分别是v+c, v-c, v。其条件数：

可见流动的速度越小，条件数越大，那么Euler方程组越发病态。

我们已经知道问题所在，但重要的是如何解决问题？

对于病态问题，没有完全有效的数值方法，除非降低问题的条件数。不过条件数是问题本身的特性，如果改变条件数，那么就改变了问题。这是一对矛盾。

观察一维定常Euler方程(1)，人们只关心最后的定常解，对于求解中间过程是否对应真实物理，并不关心。当流动趋近于真实的定常解时，。

因此，我们可以在伪时间导数项上乘以任意的矩阵，只要收敛时，该项等于零，就不影响最后的真实解。于是：

其中Γ即预处理矩阵，容易看出当伪时间导数项趋近于零时，方程(2)的解就是方程(1)的解，即真实的定常解。方程(2)整理后：

其中，。

不同的Γ对应不同的预处理方法[1]，如果预处理矩阵对应最理想的情况，能够有

，那么我们要求解的方程(3)不再「病态」，是一个「好方程」了。必须说明，预处理方法破坏了原方程组的特征值，因此中间过程没有物理意义，不能直接用于求解非定常问题。不过，幸运的是，借助A. Jameson提出的双时间步方法[2]，上述预处理方法稍加改进，即可应用于非定常问题的求解。

这里仅以一维无粘Euler方程为例，实际的三维粘性可压缩NS方程在细节的处理上会繁琐很多，但预处理方法的原理是通用的。从数值求解的角度，一个可压缩求解器添加预处理模块并不困难，将方程(1)中的矩阵A替换为方程(3)中的即可，即只需修改数值通量部分。

相关应用

图1中展示了一个典型的可压缩求解器给出的数值流场结果，Ma=0.01，均使用AUSM格式。构型是直管道带10%高度的鼓包构型。对比(a)图和(b)图，预处理后流场得到了明显的改善。该应用使用的预处理方法见文献[5]。

图1 二维Bump构型的数值模拟结果：等压线图，Ma=0.01，均使用ASUM格式。（a）无预处理；（b）预处理。

图2给展示了，数值模拟DLR-F4翼身半模的低速绕流的残值收敛历史[6]。非常明显，采用预处理后，残值收敛历史得到了非常大地改善。多重网格迭代至100步时，相比未预处理情形，预处理后的残值减小了3个数量级。

图2 预处理后的残值收敛历史[6]

图3中展示一个简化的转子附近的低速流场的涡量云图[7]，(a)图表明没有施加预处理技术时，求解器的低速分辨率很差，流场中充满了「污点」。(b)图表明，低速预处理技术应用后，流场光滑很多，「污点」几乎全部消失。

图3 涡量云图[7] （a）无预处理；（b）预处理。

延伸

无论是微/积分方程方程组还是代数方程组中矩阵的病态，本质上都是系统特征尺度的分离。这是数值模拟多尺度/跨尺度问题的本质困难之一。面对这些无法规避的困难，预处理方法是必须的，甚至是唯一的处理办法。跨尺度问题是普遍性的，因此预处理方法也具有普遍的意义。不难列举出一些其它需要预处理方法的情形：

多物理场耦合问题[4]。多物理场，也必然是多尺度的，如MHD问题[3]。
燃烧问题。燃烧过程诸多化学反应的特征尺度，通常会跨越多个数量级。
一般来说，代数方程组的规模越大，条件数也随之急剧增加，比如Laplace/Poisson方程。Krylov子空间方法在现代CFD中非常受欢迎，但在求解大规模代数方程组时，也必须结合预处理方法。
...

预处理方法原理虽然简单，但面对特定问题，如何给出最优的预处理矩阵？这需要学者对该问题具有非常深刻的物理洞见。针对交叉学科具体问题的预处理方法依然是现代CFD研究的一个核心课题。

参考文献

[1] Turkel E. Preconditioning techniques in computational fluid dynamics[J]. Annual Review of Fluid Mechanics, 1999, 31(1): 385-416.

[2] Jameson A, Schmidt W, Turkel E. Numerical solution of the Euler equations by finite volume methods using Runge Kutta time stepping schemes[C], 14th fluid and plasma dynamics conference. 1981: 1259.

[3] Ma Y, Hu K, Hu X, et al. Robust preconditioners for incompressible MHD models[J]. Journal of Computational Physics, 2016, 316: 721-746.

[4] Ferronato M, Franceschini A, Janna C, et al. A general preconditioning framework for coupled multiphysics problems with application to contact-and poro-mechanics[J]. Journal of Computational Physics, 2019, 398: 108887.

[5] Weiss J M, Smith W A. Preconditioning applied to variable and constant density flows[J]. AIAA journal, 1995, 33(11): 2050-2057.

[6] Han Z, He F, Song W, et al. A preconditioned multigrid method for efficient simulation of three-dimensional compressible and incompressible flows[J]. Chinese Journal of Aeronautics, 2007, 20(4): 289-296.

[7] Campobasso M S, Yan M, Bonfiglioli A, et al. Low-speed preconditioning for strongly coupled integration of Reynolds-averaged Navier–Stokes equations and two-equation turbulence models[J]. Aerospace Science and Technology, 2018, 77: 286-298.

来源：多相流在线