流体力学与谱方法的精度挑战：探索计算边界

工程上进行CFD计算的时候，常用的离散方法是有限差分法和有限体积法。例如FLUENT采用的就是有限体积法。但是，在需要精度很高的一些基础研究中，常常采用另一种离散方法——谱方法（spectral method），它具有有限差分法和有限体积法无法比拟的计算精度。今天我们就来聊聊谱方法。

我们以一个具体问题为例。这个问题是常微分方程边值问题。未知函数u=f(x)满足

边界条件是

这个问题是有解析解的，其解析解为

下面我们分别用有限差分法和谱方法来求解这个方程的数值解。读者可能对有限差分法已经有一定的了解，所以这样讲更有助于读者对谱方法的理解。

什么是（1）、（2）的数值解呢？（1）、（2）的数值解就是满足边界条件（2）的一个函数u*=g(x)，它使得残差

尽可能地小。不同的数值解法的差别就在于构造数值解u*=g(x)的方法不同，以及评价残差“小”的方法不同。

首先我们用有限差分法来试试。我们在[-1,1]内划分n个等距网格节点

其中h=2/(n-1)是网格尺寸。x₁=-1和x_n=1是边界节点，其余的是内部节点。将各个节点处的数值解u*的值记为u₁,u₂,…,u_n。那相邻节点之间的区间里面u*的值呢？对于有限差分法来说，其做法是在每个节点附近构造局部的低阶插值多项式。例如，在节点x_i附近，让u*等于一个二阶拉格朗日插值多项式p(x)

式子有点复杂，但其实本质很简单，就是构造一个二次函数p(x)，让它的函数图像穿过(x_i-1,u_i-1)、(x_i,u_i)、(x_i+1,u_i+1)三个点。图1以节点数n=5为例画出了有限差分法中数值解u*的构造方法示意图。可以看出插值多项式的局部性，例如x₃附近的u*是用基于x₂、x₃、x₄为插值节点的拉格朗日插值多项式表示的。可以看出u*是分段光滑的函数。

图1  有限差分法中数值解u*的构造方法示意图。蓝色的圆点代表各个节点处的数值解u*的值u₁,u₂,…,u_n，红色、绿色和青色的粗线代表局部的低阶（这里是二阶）插值多项式，而黑色的细线代表数值解u*。可以看出u*是分段光滑的函数。

数值解u*构造出来了，那怎样评价残差的“小”呢？有限差分法的评价方法比较简单粗暴，即在每个内部网格节点处令残差等于零：

下面我们就以采用局部二阶插值多项式的有限差分法来具体求解一下（1）、（2）的数值解。将（6）中的p(x)求二阶导数，并在x_i处求值，得到数值解u*在节点x_i处的二阶导数

因此，残差在各内部网格节点等于零的条件（公式（7））转化为下面的差分方程（也叫做差分格式）

此外根据边界条件（2）有

联立（9）、（10）可以得到关于u₁,u₂,…,u_n的线性方程组，求解该方程组便可以得到所研究的问题的数值解。由于局部插值多项式是二阶的，插值的时候用到三个插值点，所以我们称之为三点有限差分格式。图2是节点数n=33时的数值解与精确解的比较。表1是计算域内的最大误差随着节点数变化的情况。可以看出，节点数大于33之后，每次将网格尺寸h减半都使得误差减小到原来的1/4左右。理论上，通过泰勒级数展开分析，可以知道差分格式（9）的误差是O(h²)的，即h→0时，误差正比于h的平方，现在数值计算结果说明确实是这样。

图2  有限差分法节点数n=33时的解答与精确解的比较。采用三点有限差分格式。

表1  误差随节点数增加的变化；采用三点有限差分格式。

如果我们想加快收敛速度，那么很自然的思路是提高局部插值多项式的阶数。例如，如果我们采用局部四阶插值多项式，则数值解u*在节点x_i处的二阶导数变成

用这个公式重新构造差分格式并算出数值解，计算域内的最大误差随着节点数变化的情况如表2所示。由于局部插值多项式是四阶的，插值的时候用到五个插值点，所以我们称之为五点有限差分格式。可以看出，网格节点数大于33之后，每次将网格尺寸h减半都使得误差减小到原来的1/16左右。理论分析也可知，这种格式的误差是O(h⁴)的。

表2  误差随节点数增加的变化；采用五点有限差分格式。

从上面的分析以及数值计算结果都可以看出，局部插值多项式的阶数越高，则计算误差随网格节点数增加的收敛速度越快。我们自然想到一个极限——能否以计算域的全部网格节点为插值节点构造插值多项式呢？显然，如果以计算域的全部网格节点为插值节点构造插值多项式，那么数值解u*就不再是分段光滑的函数了，而是全局光滑的函数。这就是谱方法的基本思想——用全局光滑的函数来构造数值解u*。

下面我们就来按照这个思想试一试。我们令数值解u*等于一个函数图像穿过全部n个点(x₁,u₁),(x₂,u₂),…,(x_n,u_n)的n-1次多项式p(x)，即p(x)满足

这可以通过拉格朗日插值法构造出来。图3以节点数节点数n=5为例画出了这时数值解u*的构造方法示意图。可以看出插值多项式是全局的，数值解u*是全局光滑的函数。

图3  谱方法中，数值解u*的构造方法示意图。蓝色的圆点代表各个节点处的数值解u*的值u₁,u₂,…,u_n，绿色粗线代表全局插值多项式，黑色的细线代表数值解u*。

这样，差分格式（9）变成

因为p(x)是全局的插值多项式，所以式中的p’’(x₂)、p’’(x₃)，…，p’’(x_n-1)每个都是全部节点上的数值解u₁,u₂,…,u_n的线性组合，由于表达式比较复杂，这里不一一列出。联立（13）、（10）可以得到关于u₁,u₂,…,u_n的线性方程组，求解该方程组便可以得到数值解。其计算结果如图4所示。可以看出，结果是十分失败的——随着节点数n的增加，数值解并没有越来越接近于精确解，反而是发散了。

（a）n=9

（b）n=17

（c）n=33

图4  采用全局插值多项式的计算结果。

这是为什么呢？？其实，数值解u*的构造方法并没有问题，问题就在于我们对于残差“小”的定义不再合适了——残差在一些离散的节点上等于零，并不一定意味着残差在整个计算域[-1,1]上就一定很小。在有限差分法中，我们使用分段光滑的函数来逼近未知函数，所以残差也是分段光滑的函数，实践证明“让残差在一些离散的节点上等于零”是可行的。但是现在我们使用全局光滑的函数来逼近未知函数，因此残差也是全局光滑的函数，实践证明以前衡量残差“小”的方法失效了。

那应该怎样重新定义这个“小”呢？答案就是将残差R(x)展开为某种正交多项式系的级数形式，然后令展开式中尽可能多的低阶项等于零。例如，最常见的就是将R(x)展开为切比雪夫（Chebyshev）级数：

其各项系数的计算公式为

其中T_i(x)是i阶切比雪夫多项式。读者可以类比周期函数傅里叶级数的各项系数的计算公式。

我们知道，对于周期函数来说，光滑函数的傅里叶级数的系数随着阶数的增加迅速衰减。相似地，在[-1,1]上的光滑函数的切比雪夫级数的系数也是随着阶数的增加迅速衰减的（实际上随着阶数的增加是指数衰减的），所以，我们只要强迫展开式（14）中尽可能多的低阶项等于零，那么R(x)就可以变得很小很小。注意这里我们充分利用了“残差是全局光滑的函数”这一性质。如果残差不是全局光滑的函数（例如在前面的有限差分法中），这一切无从谈起。

由于未知量共有n个（u₁,u₂,…,u_n），而边界条件将其中两个（u₁,u_n）的值固定了，因此需要n-2个条件才能确定剩余的未知量。因此我们自然地想到，令展开式（14）中的前n-2项的系数等于零：

总结一下，在这种新的对于“小”的定义下，我们的计算方法归结为：令数值解u*是函数图像通过全部点(x₁,u₁),(x₂,u₂),…,(x_n,u_n)的n-1次拉格朗日插值多项式，其中u₁和u_n满足边界条件（10），剩余的n-2个未知量u₂~u_n-1则通过这样的方法确定：将u*代入（4）得到残差的表达式，然后再将残差的表达式代入方程（16），将积分积出来之后，转化为关于u₂~u_n-1的线性方程组。这就是谱方法中的Chebyshev-Tau谱方法。当然，Chebyshev-Tau谱方法在具体实施的时候，常常在形式上稍稍改变，即把数值解u*用切比雪夫多项式系展开，将展开式中的系数作为待求解的未知量，而不是用各节点上的数值解u₁,u₂,…,u_n作为待求解的未知量。但这只是形式上的变化，本质上没有变，所以我们不再赘述。

图5是用Chebyshev-Tau谱方法，n=17时算出的数值解。可以看出，数值解与精确解很接近，前面遇到的发散问题得以解决。

图5  采用Chebyshev-Tau方法的计算结果。n=17。

从函数图像上只能看个大概，下面我们来定量看一下Chebyshev-Tau谱方法的精度。表3给出了计算域内的最大误差随着节点数变化的情况。

表3  误差随节点数增加的变化；Chebyshev-Tau方法。

可以看出，对于谱方法来说，计算误差随网格节点数增加的收敛速度是极快的。当n从33增大到65时，误差竟然减小到原来的将近1/600。理论上来说，n→∞时，谱方法的计算误差为O(K^-n)，其中K是一个大于1的常数。这意味着谱方法的误差是以指数规律减小的，这与有限差分法中误差以幂函数规律减小的特性是完全不同的，也是任何阶的有限差分法都无法比拟的。谱方法的这种收敛特性叫做谱精度（spectral accuracy）。从前面的分析中容易看出，谱方法能达到谱精度的根本原因在于采用了全局光滑的函数来逼近未知函数。采用全局光滑的函数逼近未知函数之后，残差也是全局光滑的函数，所以残差的Chebyshev级数的系数随着阶数的增加是指数衰减的。因此，我们强迫级数的前n-2项等于零，其结果就是让残差随着n的增加以指数规律减小。

上面结合具体的问题介绍了Chebyshev-Tau谱方法，这实际上只是谱方法的冰山一角。谱方法的种类其实很多，可以从离散过程以及展开式所用的基函数这两个方面来分类。从离散过程分类，有Tau方法、Galerkin方法和配置点法三种类型。至于基函数，当计算域是有限区域的时候，经常使用的就是上面介绍的Chebyshev多项式，所形成的谱方法称为Chebyshev谱方法，但是如果计算域是无限的且具有周期性，则往往选用三角函数来作为基函数，所形成的谱方法称为傅里叶（Fourier）谱方法。此外还有无限且不具有周期性的区域、半无限区域等等，也有相应的基函数，这里就不讨论了。无论哪种谱方法，由于都是用全局光滑的函数来逼近未知函数，所以其计算精度都是达到谱精度的。

谱方法的早期研究有1938年Lanczos的工作以及20世纪60年代Clenshaw, Elliott, Fox等人的工作^[2]。但是，由于谱方法计算量比较大，一直没有引起人们的兴趣。直到1965年Cooley 和Tukey发现了快速傅里叶变换算法，才改变了谱方法的命运。人们发现，对于Fourier谱方法和Chebyshev谱方法，实际计算的时候可以利用快速傅里叶变换算法来大大提高计算效率。于是，20世纪70年代谱方法得以迅速崛起。那时候，美国数学家Steven Alan Orszag是其中最有力的推动者。他在1971年发表于Journal of fluids mechanics的文章“Accurate solution of the Orr-Sommerfeld stability equation”中（文献[8]），利用Chebyshev谱方法对平面泊肃叶流动的临界雷诺数进行了精确的计算，通过求解Orr-Sommerfeld特征值问题，得出临界雷诺数为5772.22的结果。虽然整个研究的结果仅此一个数值而已，但是其影响却十分深远。这项工作不仅是谱方法应用于流动稳定性分析的开山之作，同时实际上也确立了谱方法在科学计算中的重要地位。谱方法作为一种求解微分方程的方法，已经广泛用于流体力学、量子力学、振动、线性和非线性波等等很多领域。下面我们主要谈谈流体力学方面的应用。

谱方法的特点是精度高，随着网格加密的收敛速度快；但也有一些明显的缺点，例如不能适应复杂形状的计算域等，由于上面讨论的例子是一维的，所以看不出这个缺点来。由于这个缺点，谱方法无法应用于实际工程中具有复杂形状的流体区域。但是在很多流动机理的研究中，计算域简单，且需要高精度，正适合谱方法施展拳脚。在流体力学中，谱方法最令人瞩目的应用是湍流直接数值模拟与流动稳定性分析。

湍流直接数值模拟（Direct numerical simulation，DNS）^[4,5]就是不使用任何湍流模型，而是通过直接求解Navier-Stokes方程组来求解湍流的运动。由于湍流中的旋涡运动的尺度范围极其宽广，所以DNS需要很高的网格分辨率，因此DNS自然成为了流体力学中最需要高精度算法的领域之一。湍流直接数值模拟的网格分辨率总是随着计算机的发展而逐步增加，每个时期从事DNS的研究者都总是希望用上当时最好的计算机。以均匀各向同性湍流为例，最初（1972年）Orszag和Patterson对于均匀各向同性湍流的模拟^[7]是在32³的网格上进行的计算，而21世纪初，对于均匀各向同性湍流的最大规模的计算已经用到了4096³的网格。尽管如此，离模拟工程上的高雷诺数湍流还差很远。在对物理现象的分辨率不变的前提下，提高算法的精度，就可以采用相对疏一些的网格，从而减小计算时间和内存消耗。

图6是She等^[6]1990年利用谱方法对均匀各向同性湍流进行直接数值模拟得到的涡结构图，在该模拟中，计算域为一个正方体，x、y、z每个方向上均采用周期性边界条件，网格分辨率为96³。由于边界条件为周期性，所以采用了Fourier谱方法。

图6  用谱方法进行DNS得到的各向同性湍流中的涡结构（复制自[6]）。该工作发表于1990年Nature。

谱方法在流体力学中还常常用于流动稳定性分析。流动稳定性研究的是在Navier-Stokes方程组和某些边界条件下，数学上适定问题的层流解答，在施加扰动的情形下，扰动衰减还是放大的问题。公众号之前的文章“层流为何会转变为湍流：托尔明-施利希廷波的故事”就是这类问题的例子，即平板层流边界层流动在施加扰动的情况下，扰动最终究竟是会衰减还是会放大。这个问题的重要性在于它揭示了湍流的产生机制，即为什么随着流动雷诺数的增大，稳定的层流流动会失稳转化为湍流流动。

流动稳定性分析往往最终归结为微分方程的特征值问题，例如平行剪切流线性稳定性问题最终归结为Orr-Sommerfeld特征值问题。对于特征值问题来说，计算精度十分关键，这是谱方法在这类问题中得到重要应用的原因。历史上，Orr-Sommerfeld特征值问题的求解曾经遇到过很大的困难，还吸引过像Heisenberg（海森堡）、Tollmien这样的大科学家为之努力。其解析求解方法在数学上十分复杂。但是，随着高速计算机的发展和数值方法的进步，求解Orr-Sommerfeld特征值问题已经不再困难^[9]。今天，谱方法已经称为求解Orr-Sommerfeld特征值问题的标准方法。

墨尔本大学的研究生刘丽媛和北航宇航学院的研究生胡涛阅读了本文的初稿并提出了很好的修改意见，在此表示感谢。

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删