高等计算流体力学-03

1、求解流体力学方程组的有限差分方法,第三讲,1,控制方程,2,贴体网格 H型 C型 O型,3,一、贴体坐标中的基本方程,正变换,逆变换,假定已知上述变换的离散形式 即已经生成好贴体网格,4,导数的变换,一阶导数,5,二阶导数,6,二阶导数,度量系数 metrics,7,度量系数及其计算方法,困难,容易,如果能够找到 和 之间的关系， 就可以得到 等的计算方法,8,Jacobi行列式 （jacobian）,9,10,11,12,任意曲线坐标系中流体力学方程组的守恒形式,13,14,15,曲线坐标系中NavierStokes方程的强守恒形式,16,二、守恒型Euler方程,在任意曲线坐标系中，二维守

2、恒型Euler方程为：,其中,17,二维守恒型Euler方程可以改写为拟线性的形式：,其中,令 或 ，则 可以写为统一形式：,18,其中， 为绝热指数，对于完全气体， ， 和 的定义为：,的特征值为：,其中 是音速。,19,三、守恒型的NavierStokes方程,平面上NavierStokes方程的强守恒形式为：,其中,注意到，Euler方程可以看作NavierStokes方程 的特例。一般而言，Euler方程中无粘通量 的离散 方法，同样可用于NavierStokes方程；而为了求 解NavierStokes方程，我们还必须提供粘性通量 的离散方法。,20,有限差分方法,21,有限差分方法

3、,对 式直接进行半离散 差分近似,注意到 ，采用中心差分，则有：,令，,则可以得到守恒格式：,22,在有限差分方法中，在网格点上的度量系数(如 )也需要通过差分 方法计算，因此要求计算网格是充分光滑的。,23,有限差分方法的主要特点是： 1）有限差分方法只需构造偏导数的离散方法，这使得它比较容易推广到高阶精度，对于多维问题也是如此。 2）有限差分方法对网格的光滑性有较高要求，不易在复杂形状求解域上实施。 3）在曲线坐标系中，有限差分方法要对几何量和物理量的确定组合进行差分离散，这样做的后果之一是有限差分方法可能产生所谓几何诱导误差（geometrically induced errors） 。

4、例如，考虑一无 界均匀流场，易知，该流动应是定常的。即 。在由一般曲 线坐标构成的网格上采用有限差分方法计算此流动，可能有 （三维 以及动网格中的有限差分方法更容易出现这种情况）。而采用有限体积方 法，则恒有 。有限差分方法在某些情况下不能复现均匀流 场，是几何诱导误差的表现之一。,24,如何计算数值通量？,数值耗散 数值色散,线性: 分辨率, 精度 非线性: 激波捕捉,分辨率, 相位畸变,简单中心差分一般来讲不是好的选择!,25,激波间断和广义解,26,可压缩流动中一个非常重要而又复杂的现象是激波现象。 在NavierStokes方程的范畴内，激波可以看作是一个连续但物理量的梯度非常大的有限

5、厚度的结构。 在Euler方程的范畴内，由于方程本身缺乏必要的耗散机制，激波厚度为零，激波两侧物理量存在间断。 在数值求解Euler方程时，必须解决包含间断的流场的计算问题。 在数值求解NavierStokes方程时，虽然理论上激波附近物理量是连续的，但是在通常情况下，激波的厚度非常薄（为分子平均自由程的量级），远远小于计算网格的尺度；所以在计算中，我们仍需把激波当作间断来处理。,在Euler方程的范畴内如何表示含激波间断的流场?,27,一、激波的形成,激波是一种非线性现象。 考虑线性对流方程的初值问题，如果其初始值是连续的，则解中不会自发产生间断。 但是，对于非线性问题，即使初值是连续的，在

6、后续时间中，也可能发展出包含间断的解。 Bergers方程可以看作Euler方程的非线性模型方程；和Euler方程一样，是双曲型方程。 下面，我们以Bergers方程为例，介绍间断解或者激波的形成过程。,28,考虑Bergers方程,的初值问题，初始条件为 。根据特征理论，Bergers方程的解析解可以写为下列形式：,初始条件,Bergers方程的解为：,29,图1 Bergers方程初值问题的解,图1显示了解随时间的发展过程。容易知道， 当 时，Bergers方程的解为,可见,对于非线性的Bergers方程,即使初始值 是连续的,其解仍然可能出现间断。,事实上，一般的非线性双曲型守恒方程（组

7、）或称非线性双曲型守恒律的解中都有可能存在间断。 通常，我们认为偏微分方程的解是连续可微的。显然，对于非线性双曲型守恒律，这一点并不成立。 因此，必须拓展双曲型守恒律解的概念。,30,二、广义解,对于一般的一维双曲型守恒律 的初值问题，我们引入广 义解（generalized solution）或称弱解(weak solution)的概念，来表 示包含间断的解。,定义：设 是分片连续可微的函数，在 的半平面，如果对于与 的间断线只有有限个交点的任意分段光滑的闭曲线 ，都有: 则称 为方程 在初值 下的广义解或弱解。,31,下面讨论一下广义解的意义。 如果已知 是光滑的，设 围成的 区域为 ，则

8、由 式利用Green公式知,由于闭曲线可以在光滑区内任取，可得：,在光滑区，弱解就是通常的连续可微解。,32,如果 是由一条间断线 分隔开的分片连续可微函数，取如图所示的由 （ 是正数）围成的闭曲线 ，在 上应用 ,有:,33,间断关系的推导,令 ，并考虑到在 两侧 有间断， 则上式可简化为：,分别表示从左、右趋近间断线时守恒变量的值。,令,在间断两侧，显然 。 考虑到 可以任意取值，可得：,其中,Rankine-Hogoniot（R-H）关系,34,以一维Euler方程为例，R-H关系为：,注意到，当 时，如果 ，则上面的关系成立，此时只有密 度场存在间断。这种类型的间断称为接触间断。不属于

9、接触间断的间断称 为激波。,综上所述，双曲型守恒律的广义解 是被有限个间断线分开的分片光滑 函数。在光滑区， 是满足微分方程式 的连续可微解；在间 断线的两侧， 满足R-H关系。,35,三、熵条件,广义解推广了偏微分方程初值问题连续可微解的概念，但其后果是 导致了弱解不唯一。为了说明这一问题，我们举一个例子：考虑 Burgers方程在初值为 时的解,此时，初值在 处有一个间断。,把Bergers方程写为守恒形式,可知 处的Rankine-Hogoniot条件为：,由上式知 。,在间断处满足Rankine-Hogoniot条件，在其他地方满足微分方程， 即 是Bergers方程的一个广义解。,3

10、6,另外，考虑,分片连续可微函数，其导数在 处存在间断(或者认为函数值存在幅度 为零的间断)。,根据广义解的定义，容易验证此函数也是Bergers方程的一个广义解。所以 广义解一般不唯一，但是对于有明确物理意义的守恒律，其中只有一个解 是有物理意义的，我们称之为物理解。为了得到我们关心的物理解，广义解 除了必须满足式外，还必须满足附加的条件，这个条件因为与热力学第二定 律所起的作用相同，被称为熵条件（entropy condition）。,37,如何在不唯一的广义解中挑选所谓物理解呢？最直观的方法是认为物理解是 一个光滑的粘性问题的解在粘性趋于零时的极限。也就是说：,方程 的解如果当 时，几乎

11、处处有界的收敛到 分片连续可微函数 ，则 是 的物理解。,这种方法的物理基础是：在真实的物理过程中，一般会包含某种耗散机制。 而双曲型守恒律往往是在建立数学模型时忽略这种耗散机制的结果。例如， 我们前面说过，对Navier-Stokes方程而言，激波是有一定厚度的连续结构， 但忽略粘性效应以后的Euler方程，解中则可能出现理想的间断。很明显， Euler方程的物理解应该是Navier-Stokes方程在粘性趋于零时得解。,38,这种选择物理解的方法，虽然非常直观，但并不实用，因为粘性问题的解 一般也是未知的。对于标量非线性双曲型守恒律，Oleinik 提出了下面的熵 条件：,熵条件：设 是定

12、义在t0上半平面上，存在有限条光滑间断线的分片连续 可微函数;且是标量守恒律： 的弱解,若 在间断线附近满足,其中 。 则弱解是唯一的，并且就是物理解。,39,40,41,含激波的流场的计算方法,42,包含激波的流场的数值计算方法可以大致分为两类，即激波装配（shock fitting）方法和激波捕捉（shock capturing）方法,43,激波装配方法中，激波把求解域分为两个或者多个子域，激波作为各子域 边界的一部分，激波的两侧通过R-H关系相联系。 在每个子域内，流场是光滑的，可以用任何计算可压缩流动的数值方法求解； 在激波处，用R-H关系和必要的补充条件计算出激波运动的速度，从而可以

13、在计算中不断更新激波的位置。 激波装配方法的优点是：可以比较准确地计算出激波位置和光滑区的流场。但缺点是计算过程比较复杂，不容易设计通用的计算软件。 显然，如果在各个子域分别生成计算网格，则网格一般是随激波运动的；如果各个子域共用一套静止的网格（称为浮动激波装配方法），则激波附近的计算格式要进行特殊处理。 另外，当流场中激波的结构比较复杂，或者存在激波的产生和消失等过程时，激波装配方法将更难应用。,44,激波捕捉方法的基本思想是：在计算包含激波的流场时，采用统一的计算 格式，不对激波进行任何特殊处理。 当计算格式满足一定要求时，可以自动计算出流场中的间断。 这种方法的优点是：编程计算比较简单，

14、且可以适用于具有任意激波结构（包括激波的产生、消失、运动等各种情况）的流场计算。 其缺点是：此时激波不再是理想的间断，而是厚度为一至数个网格宽度的连续变化的结构；另外，由于激波内部流场的梯度较大，常规的计算格式或者会把激波“抹平”，或者在激波附近产生虚假数值振荡。 因此，为了高分辨率、无振荡的计算激波，对数值计算方法提出了很高的要求。 事实上，能有效计算激波的高分辨率格式，是最近30年以来，计算流体力学研究的热点问题之一。 随着高分辨率格式研究的进展，目前，激波捕捉方法已经成为计算包含激波的可压缩流场的主流方法。,45,一、守恒格式和Lax-Wendroff定理,在求解守恒律,我们可以定义一类

15、特殊的格式，称为守恒格式。,的有限差分和有限体积格式中,定义：一维守恒律 式的差分或者有限体积格式,称为守恒型差分格式，如果,且,等称为数值通量。,46,守恒格式的特点是： 对于任意的 ，数值通量用相同的函数关系计算， 且函数涉及的自变量（离散点上的数值解）与 的相对位置关系不变。 这一特点使得相邻网格点，如 和 点上的计算格式在公共界面 处的数值通量是相同的。稍后我们将说明，正是这个特点使守恒格式具有自 动捕捉间断的能力。,47,守恒格式的相容性条件,48,49,因此，,50,整理上式，可得：,当数值通量满足条件,数值通量的相容条件,51,相容的守恒格式具有自动捕捉激波和接触间断的能力。,5

16、2,由于 的任意性， 可以代表 平面的任意闭曲线；这样，当 时，相容的守恒格式的数值解满足弱解的定义。即当 时，相容的守恒格式的解就是守恒律的广义解或弱解。 这一结论可以概括为下面的定理：,由于守恒格式的数值通量的特点，内部的数值通量可以相互抵消，从而有：,53,Lax-Wendroff 定理:如果 是守恒律 初值问题的相容守恒格 式的离散的数值解，且当 时， 在某种范数的意义下趋于 ，则 是守恒律的一个弱解。,提示 （1）Lax-Wendroff 定理并不保证弱解的存在性和相容的守恒格式的收敛性。同时，该定理也并不保证得到的弱解满足熵条件。但是，通过该定理我们可以确信：用相容的守恒格式可以自

17、动捕捉激波和接触间断。 （2）Lax-Wendroff 定理给出了用数值方法计算弱解的一个充分条件，也就是说，该定理并未排除用非守恒型格式计算出弱解的可能性。事实上，在某些情况下，非守恒格式也可以得到弱解的良好近似。但是，非守恒型格式是否可以计算弱解，对具体的格式要进行具体分析，而这种分析通常是复杂的。所以，我们一般用守恒格式计算包含间断的流场。,54,二、人工粘性和格式粘性,在计算包含激波的流场时，考察计算格式的“数值粘性”特征是非常重要的。 这主要出于两个原因：（1）物理解是所谓“粘性消失解”，所以，在计算包 含激波间断的流场时，计算格式包含足够的数值粘性是必要的；（2）某些 计算格式在间

18、断附近会产生非物理振荡，合理地调整计算方法的数值粘性， 对于消除计算结果在间断附近的数值振荡是有益的。,55,为了保证格式是相容的，应该有：,在能够保证计算稳定和抑制间断附近振荡的前提下, 人工粘性和格式粘性越小越好!,56,中心型格式需要添加人工粘性; 迎风型格式自带格式粘性!,先考虑中心格式,57,考虑Euler显式格式,其修正方程为：,这是一个一阶精度的格式，二阶耗散项对数值解的幅值有主要影响。 该格式二阶耗散项是负的，即格式粘性为负，因而是无条件不稳定的。 为了使计算稳定，我们可以在格式的右端添加一个二阶导数项，如,此格式就是前面介绍过的Lax格式，其修正方程为：,则当 时，格式是稳定

19、的。,58,注意到：,而且上述两个格式修正方程之间的,，即,差的主项正好为,显然在Euler显式格式的右端增加,相当于增加了Euler显式,格式的耗散，所以这一项通常称为,Euler显式格式的“人工粘性”。,由于添加了这种人工粘性，使得格式是稳定的，而且间断会被抹平。,推广到Euler方程，这种添加了人工粘性的Euler显式格式（实际上就是 Lax格式），可以写为,上述格式可以改写为：,其中,所以，Lax格式是守恒格 式，可以自动捕捉流场中 的间断。,59,线性波动方程的Lax-Wendroff格式为,其修正方程为：,可以看到，这是一个二阶格式，其格式粘性是四阶的。四阶格式粘性 通常不足以抑制

20、间断附近的振荡。为了解决这一问题，我们也可以采 取添加人工粘性的方法。,例如，可以把差分格式改写为：,其中 为二阶粘性系数。如果 ，虽然可以有效的消除间断附近的 振荡，但是，间断被严重抹平，而且格式变为一阶精度，在光滑区的精度 也大为降低。因此，这种方法显然不是一个好的选择。,60,比较好的取法是取,其中 为量级为 的正数， 是一个小的正数，如 ，以防止出现分 母为零的情况。对于这种取法，在光滑区，有 ，则人工粘性项为 的量级，从而不会影响光滑区格式的精度。在间断附近， ， 人工粘性为 的量级， 激波附近的振荡可以得到有效控制。注意到， 在人工粘性项中一个自由参数 ，一般可取1/2，也可以根据

21、试算调整 ， 以达到更好的效果。但是， 对于不同的问题，最佳的 值是不同的。,在添加上述人工粘性项以后，我们并不能从理论上保证间断附近解没有振 荡。而且由于可调参数 的存在，使得这种方法的实施在一定程度上依赖 于经验。,61,推广到Euler方程，有,其中，,是Jacobi矩阵， 是压力， 为Jacobi矩阵 中所有特征值中绝对值最 大者。注意到，为了保证格式是守恒格式，Euler方程人工粘性的形式与 线性波动方程略有不同。,62,中心型格式需要添加人工粘性; 迎风型格式自带格式粘性!,再考虑迎风格式格式,63,迎风格式如何应用于方程组?,64,四、一维线性波动方程组的迎风格式,考虑线性双曲型

22、方程组,由于Jacobi矩阵是常数矩阵，上式也可以写为守恒形式,由于双曲性假设，矩阵 可以对角化：,其中 是 的特征值， 是由特征值组成的对角阵， 和 分别是 的左特征 向量和右特征向量组成的矩阵，且二者互为逆矩阵：,65,令,其中W 称为特征变量，则有,写成分量形式,其中,为W的第i个分量。,在形式上与线性波动方程相同,因此，线性波动方程的迎风格式可以直接应用于,以一阶迎风格式为例，有：,66,写成向量形式，有,其中,是特征值的绝对值组成的对角阵。,左右同乘以,写成守恒变量的形式，有：,67,68,69,注意到，对于常系数的双曲型方程组，有,70,SCM：Splitting Coeffici

23、ent Matrix,71,五、 Euler方程的迎风型有限差分格式,考虑一维Euler方程,Euler方程是非线性的方程，如何把上面所述迎风格式推广到求解非线性双曲型方程组特别是Euler方程呢？我们下面介绍几种方法。,72,1. SCM方法,和线性情况不同，上式中的矩阵A不是常数，而是U的函数。但是，我们仍然可以借鉴线性双曲型方程组迎风格式的思路，把矩阵A进行分裂。其做法与线性的情况完全类似，唯一的区别是由于A在各个网格点是不同的，所以，在每个网格点均需进行矩阵分裂。这种方法称为分裂系数矩阵方法（Split Coefficient Matrix Scheme），简称SCM方法。,把Euler方程写成拟线性形式,应用SCM方法的一阶迎风格式为：,，,73,2. FVS方法 FVS是矢通量分裂（Flux Vector Splitting） 的简称，最早由Steger-Warming提出。我们注意到，虽然SCM方法可以看作线性双曲型方程组迎风格式的直接推广，但是却不是守恒格式。为了构造守恒的迎风格式，Steger 和Warming注意到Euler方程的通量F是守恒变量U的一次齐次函数，即满足,在这种条件下，有,这一关系可证明如下：,则,74,利用这一性质，我们可以把通量分解为：,其中,所以，Euler方程可以改写为：,同线性双曲型