2.3　偏微分方程

在CFD和海洋数值模拟中，与其他许多包含流体流场计算的研究一样，我们主要使用偏微分方程求解，这些方程通常是非线性的（也就是会出现包含因变量和（或）因变量的导数项），它们不服从解析解，因此需要使用数值方法进行求解。对于由Navier-Stokes方程控制的流体流问题而言，非线性平流项是产生困难的主要根源，然而，这些项包含一阶导数，不是方程中的最高阶项，最高阶项通常是包含二阶导数的扩散项，因此是这些扩散项决定了方程的性质，也决定了对其进行求解的合适方法。这可以用一个简单的单一因变量φ（x，t）（其中t可以是时间也可以是空间）的线性二阶偏微分方程进行说明

该方程的性质只与包含最高阶导数的项有关，即前三项，其余项无关紧要（即使较低阶的项是非线性的也一样）。在数学上，x-t空间中存在两个基本方向，用特征方程的解对其进行定义为

这两个方向称为特征方向，由下面的式子得到

如果两个特征是不同实数［（B2-4AC）＞0］，则方程是双曲线；如果是两个相同实数（B2-4AC=0），则方程是抛物线；如果是复数［（B2-4AC）＜0］，则方程是椭圆。

在物理学中，特征定义了一条线，信息在物理区域中（x-t）沿着该线进行传播。波传播问题（超音速流是一个典型的实例）是由具有实特征的双曲线方程问题。稳定状态问题（Stommel问题和Munk问题是经典的例子）通常由具有虚特征的椭圆方程描述。特征的实数值或复数值性质对问题的性质非常重要，同时对用来进行求解的方法也很重要。

对于含有一个耦合偏微分方程组和所包含因变量不止一个（假设为n）的真实的流问题来说，其控制方程可以转化为一个一阶偏微分方程组如下

其中Ф是一个包含因变量的矢量，矩阵A和B是因变量x，y和时间t的函数，P为另一个矢量。该方程的性质与矩阵A和B的特征值有关，如果A的特征值是不同的实数，则问题是关于x和t的双曲线，如果是复数，则是关于x和t的椭圆，而B则控制着y-t空间中解的性质。

式（2.3.4）中只含自变量x和y的子集定义了二维空间中的稳定状态问题，它的特征方向由下式（Hoffman和Chiang，1993）的解来决定

该方程特征曲面法向量具有n个解的n阶代数方程，特征曲面就是特征方向my/mx。如果所有特征都是实数，则系统为双曲线；如果都是复数，则为椭圆。

2.3.1　一致性、收敛性以及稳定性

该问题的适定性（Hadamard定义）也是很重要的，对于与CFD相关的数学问题来说，适定性就是不仅要求解存在并是唯一的，还必须持续依赖于辅助数据。换句话说，控制方程并不能单独决定唯一解，而是和辅助条件（初始条件和边界条件）一起决定唯一解。一般数值问题要求满足一致性、收敛性和稳定性几个条件，也就是对控制微分方程的离散逼近必须在一种意义上是一致的，即令离散间隔任意趋于0时可以对原始微分方程进行还原，这一点保证所用的形式逼近所考虑的微分方程，而不是其他方程。可以通过离散过程求逆来对一致性进行检查，具体方法是用邻近格网点替换泰勒级数展开式，并要求时间步长Δt和格网大小Δx，Δy以任意方式趋于0，结果应该是原始偏微分方程，通常使用泰勒级数展开式求控制微分方程的有限差分近似。因此，通常情况下一致性不构成问题，只有在极少的形式中会遇到一致性问题，一个罕见个例是抛物线方程的Dufort-Frankel方案。

计算方案也必须在一种意义上满足收敛性，即离散间隔减小时，数值解应该收敛于真实解。通常情况下，对于复杂系统来说收敛性很难证明，虽然通过使用不同格网大小（只受限于增长的累计舍入误差）进行重复计算经验性地检查收敛性的原理比较简单，但由于资源需求大，一般很少从这方面进行验证。因此，很难保证收敛性。然而，Lax提出的定理对此有很大帮助，Lax等价定理表明，如果一个方案满足一致性和稳定性，那么它也满足收敛性。对一个满足一致性的方案来说，只需要对该方案的稳定性进行检查，推理如下：如果一个方案是一致的，那么可以使用泰勒级数展开式对离散过程求逆从而得到控制偏微分方程，其次，如果该方案也是稳定的，那么格网的逐步求精应该能使有限差分解至少在理论上收敛到精确解，因为截断误差单调递减到0。实际上，舍入误差会产生干预，但是，收敛性的检查依然依赖于该理论，而很少使用经验性手段。大多数时候，模型分辨率的选择与可用的计算机资源以及稳定性需求有关。

计算算法也必须满足稳定并不易受到舍入误差的积累和增大带来的数值问题的影响。在求解过程中，解的任何自发性扰动必须在任一点处无增大趋势，否则，将会导致计算误差增大，无法得到真实解，从而使最终结果失去意义。这通常对允许的时间步长强加了一个限制并降低了方案的效率。一个离散方案的稳定性非常重要，最常用的方案是Von Neumann方案，然而它的适用范围仅限于线性系统和域中的内点，非线性问题必然要调用局部离散化方法。

Von Neumann的方法依赖于一个事实：由于方程具有线性性质，解的误差即数值解和真实解的差异和数值解本身一样，也受到相同方程的控制。现在如果认为误差由具有不同波长nΔx的分量构成，并且如果在所有分量在时间上有界的条件下进行研究的话，将得到该方案的一个稳定状态。L.F.Richardson在第一次世界大战期间最早尝试了数值天气预测（Richardson，1922），他使用了非稳定方案（以及控制方程的不合适形式），因此虽然他的预测完全是开拓性的，却没有多大用处。1928年Courant，Friedrichs和Lewy推导出了显式方案稳定性的时间步长约束，20世纪30年代Carl-Gustaf Rossby的工作让大气预测涡度守恒方程的应用成为可能（取代原始方程惯性——重力波问题引起的困扰），该方程将棘手的重力波问题去除掉。然而，成功的数值天气预测直到40年代晚期电子计算机问世后才成为可能（Richardson估算过，别说预测，仅仅跟上全球范围内最新的天气情况就需要64,000个真人“计算机”），Charney、Fjortoft和von Neumann在普林斯顿使用涡度守恒方程和第一台电子计算机（ENIAC）首次非常成功地对天气进行了预测（虽然从某种程度上说只是对于重力势的预测）。

在控制变量u（x）的时间演变的一维问题中，Von Neumann技术包括=ξnexpi（jkΔx），k=1，2，…的替换，其中ξ（=ξn+1ξ-n）是一个复数，表示第k个空间模式的振幅，称为放大因子，如果ξ保持有界，即对于波数k的所有可能取值|ξ|≤1，计算误差或者减小，或者不再增大，此时解是稳定的。这一标准的检查非常简单，尽管其代数方法比较麻烦。需要注意的是，ξ通常是复数。对于因变量大于1的方程系来说（比如CFD中的方程），Von Neuman方法得到一个放大矩阵ξ，而不是一个标量放大因子，然后需要求取该矩阵的特征值λm，其中m为变量个数，稳定条件为：对于所有m和波数k的所有可能取值，|λm|≤1。

可以用一维线性平流方程对Von Neuman方法进行说明，该方程的精确解是u=f（x-ct），u=f（x）表示初始条件，随着时间的增长，初始形式在x的正方向上以平流速度c进行不变的传播。利用时间步长Δt和格网大小Δx，如果cΔt=Δx或者Courant数C=cΔt/Δx，很容易得到穿过一个格网的真实解等于1。很显然如果对数值解进行求取，那么将解从一个格网点到另一格网点进行传播的速度不可能超过真实解本身的传播速度，也就是说，在解变为非物理的之前，C不可能大于1，通俗地说，所考虑的点的计算依赖域必须包含依赖的物理域（图2.3.1），实际上，对于随机舍入误差来说变得不稳定的解会导致计算结果的失败。

要在数学上对其进行证明，则可以将解看成不同波数的傅里叶级数的重叠并确保所有级数都不随时间而增大。由于方程的线性性质，该误差也满足相同的方程，这也意味着数值解中的误差保持有界，可以用平流方程对此进行说明。考虑波数k的单一级数为

将其带入式（2.5.3）得到振荡方程为

它的解为U（t）=Re［U（0）exp（-ikct）］。将式（2.3.6）和式（2.3.7）结合很容易看出，解是传播保持不变的初始形式。在此，考虑式（2.3.6）的有限差分等价方程为

使用式（2.3.7）的中央（跳步）时间差分为

如果用U（n+1）=ξU（n）定义一个放大因子ξ，则数值解的稳定性要求|ξ|n=|U（n）|/| U（0）|有界，或者说当Δt/t较小时，|ξ|≤o（Δt/t）或|ξ|≤1+o（Δt/t），这是稳定性的Von Neumann必要条件。然而，对于随机的舍入误差和其他计算误差而言，充要条件是|ξ|≤1，因为ξ通常是复数，这意味着在复数平面（Reξ-Imξ）中，对于方案的稳定性，ξ的所有可能取值都必须位于以原点为中心的半径1的圆内。式（2.3.9）得到跳步中心的空间差分方案为

如果对于kΔx的所有容许值有| Csin（kΔx）|≤1，则满足稳定性的充要条件|ξ|≤1，因为kΔx可假设值为0～2π，则sin（kΔx）的最大可能取值为1，因此方案的稳定性要求C≤1，方案是条件稳定的。

虽然Von Neumann方法在概念上比较简单，但实际上，它只有在内部域并且不把边界条件的影响考虑在内时有效。矩阵方法确实依赖于将时间步长为n+1时的解用步长为n时的解进行表示的形式为

考虑了空间域中所有格网点处的解。此时，U是一个长度为N的向量，N为格网点个数，A是N×N矩阵，A的形式可以很容易地用边界条件下的差分方程得到，很显然U（n）=A（n）U（0），因此A是放大矩阵，该方案的稳定性取决于该矩阵的性质，解有界的充要条件是该矩阵的值大小不超过1，即|λn≤1|，n=1，…，N，N个特征值由特征方程|A-λI|的根得到，一般而言，就像上面对Neumann的分析那样，它们会通过较小的O（Δt/t）超过单位一。而且，如果放大的本征模在初始条件下不存在，则解是不稳定的，但事实上，求解过程中的舍入误差必不可少地要求特征值都不大于1。这个过程可以用扩散方程［见下文式（2.5.21）］进行说明，为了简单，使用两级FTCS（时间上为前向，空间上为中央）的方案，即

（N-1）×（N-1）阶放大矩阵A变为

其中D=vΔt/Δx2，使用的边界条件=0。该矩阵是三对角矩阵，它的元素为a（=D），b（=1-2D）和c（=D），它的特征值用下面式子得到（Haltiner和Williams，1980）

如果D≤1/2，则所有特征值有界，因此该方案是条件稳定的。

Nuemann方法和矩阵方法都是适合于线性方程的最好方法。第三种方法适合于非线性方程，它也是能量方法，对因变量（相当于能量）的平方和的有界性进行检查。如果能量有界，那么在每个点处的解也有界，并且具有稳定性。这可以用式（2.5.3）的FTUS（时间上向前，上游差分）进行说明，因为对它来说，与其他方案相比这种方法更为简单（Mesiniger和Arakawa, 1976；Haltiner和Williams, 1980）。

由此得到

使用循环边界条件得到利用它和式（2.3.16）中要求（1-C）C＞0的Schwarz不等式∑ab＜（∑a2）1/2（∑b2）1/2，得到

因此，只要满足该不等式，那么方案就是稳定的，当然，这要求CFL条件C＜1。

2.3.2　椭圆系统、双曲线系统及抛物线系统

椭圆方程问题包含物理域边界条件，在该邻域中对方程进行求解。这些边界条件可以是狄利克雷边界条件，其中规定因变量的值位于边界上，也可以是Neumann条件，对其导数（通常与边界垂直）进行描述，在Nuemann条件情形中，可以确定解是积分的一个未知常数。边界条件还可以是包含因变量与其导数线性结合的混合（或Robin）条件。所需边界条件的数量取决于方程的阶数。调和方程只包含二阶导数，它只需要因变量的一个条件。例如，在椭圆拉普拉斯方程描述下的稳定状态流体流问题中，拉普拉斯方程为

或者等式右侧被f（x，y）替换的泊松方程，它足以对φ或者其边界上的导数进行描述。此外，还必须为以下的双调和方程描述两个条件

需要在边界上对φ及其导数进行描述。

双曲线方程多少包含较少耗散或无耗散的传播。物理域中某一给定点处的扰动与两个特征一起，从该点以一个有限的速度向该域的“下游”进行传播，因此这种扰动的影响只有在被这两个特征划定界限的楔形域中才会产生。同样地，物理域中的一个点只会受到以两个特征为边界的楔形域中向“上游”传播的扰动的影响（图2.3.1，考虑x正方向上的超音速流，对此的理解便很清晰）。由于这种方向性，可以在受稳定性标准控制的一致向下游的适当增大条件下，用已知的上游条件和控制方程的逐步积分进行求解，此处的稳定性标准包含媒介中扰动（c）的传播速度。涉及时间和空间的双曲线问题是初始值问题，需要在已知的初始条件前向时间下用积分逐步对它进行求解，寻找所要求的满足物理域中边界条件的解，求解过程中的时间步长受到微扰动的传播速度（在超音速流中为声速）的影响，从物理学上来看，就是对于选定的格网大小来说，数值解的前向时间传播速度不能小于扰动本身在该格网中的传播速度：时间步长Δt以Δx/c为界，Δx为物理空间中的格网大小，相当于柯朗数C=cΔt/Δx必须小于1。更确切地说，计算域必须包含进行时间步长点的物理域，这就是著名的Courant-Friedrichs-Lewy（CFL）条件，它用包含偏微分方程的时间步进式问题对解的效率进行约束，在显式时间步进方案中，每到下一时间步长t+Δ t时使用一个物理点（一个x）对解进行确定。

可以使用隐式时间步进方法来避开这一限制，但是这一方法需要在下一个时间步长时求取所有物理点处（所有x）的解。因此，尽管时间步长不受稳定性的约束，仅受到过渡解的精度要求，但是解通常要求对一个大而稀疏的矩阵求逆，通常很少使用蛮力求逆手段比如传统的高斯消元法，而是使用迭代方法从初始估计开始，最终收敛到正确解。使用哪一种方法主要取决于简化程度和计算效率的要求。通常在涉及对方程的二阶导数项进行二阶差分的问题中（例如，涉及垂直方向上混合的问题），所得矩阵为三对角矩阵形式或者其他同样简单的矩阵形式，因为这种形式可以使用快速算法。否则，能否采用快速矩阵解算器或者能否在所用的迭代方案下使解快速地收敛就成为一个问题，可用的计算机体系结构也是一个问题。因为不同过程之间具有相互联系的需求，大型并行计算机的体系结构本质上不太适合对涉及全局依赖的问题进行求解，只包含局部依赖的显式方法使用这种体系结构则会更加有效。

抛物线问题是一种典型问题，它涉及空间域中的性质随时间增长而产生扩散的问题，这种扩散的经典实例有热量扩散、动量的黏性扩散和涡度的黏性扩散等。解的趋势是最终消除由物理过程的扩散性质引起的任何初始不连续性（而在双曲线问题中，低阶的非线性平流项能产生不连续性，即使初始状态下没有任何不连续性存在）。时间依赖问题的稳定状态限制涉及的空间方向不止一个，并受一个双曲线偏微分方程控制，将其降低为一个椭圆方程。

因此，抛物线问题和双曲线问题都是“时间依赖”问题，而椭圆问题是一个“稳定状态”问题。数值上，抛物线问题的求解方法与双曲线问题的求解方法没有多大差异，它们都是时间依赖初始值问题（初始值问题s），这是最重要的。而椭圆问题采用的求解方法则有很大不同，因为它是稳定状态边界值问题（边界值问题s），它要求整个物理域上的解在边界处服从一些条件，通常从一个初始估计开始进行迭代求解，最终收敛到真实解。在椭圆问题中，物理域中任何点处的一个扰动也会同时在该域中的其他所有点产生，这与时间依赖的抛物线问题和双曲线问题不同，这两个问题中可以每次对一个物理点进行求解。在数值上，这意味着将椭圆偏微分方程转化为一个代数方程组，因此也“转化”一个大的矩阵。实际上，迭代方法被用来进行求解，在此，解的收敛速度再次成为一个关键问题。

当进行数值求解时，初始值问题与边界值问题之间的差异通常变得模糊起来，因为它们都可以使用迭代矩阵解算器进行求解。实际上，通过对边界值问题添加含有因变量的时间导数的项以及从一个估计的初始状态开始进行积分直到一个稳定的渐进状态，可以将一个边界值问题转化为一个初始值问题，这称为伪瞬态方法。事实上，边界值问题的迭代松弛法可以看成是在时间上进行迭代的方法，换句话说，就是一个初始值问题。