第六章+计算流体力学的基本方法

计算流体动力学

Computationalfluidmechanics机械与动力工程学院凌祥2第六章流体力学的控制方程拉克斯-温德罗夫（Lax-Wendroff）方法麦考马克（MacCormack）方法粘性流动、守恒形式和空间推进松弛法及其在低速无粘流动中的应用数值耗散、色散及人工粘性交替方向隐式方法压力修正法及其在不可压粘性流动中的应用用于CFD的计算机绘图技术3拉克斯-温德罗夫方法拉克斯-温德罗夫方法是一种显示有限差分方法，特别适用于推进求解。通过沿时间或空间逐步推进而得到数值解的接法是与双曲型和抛物型偏微分方程的求解相关联的。4拉克斯-温德罗夫方法举例：非定常二维无粘流这个问题的控制方程：连续性方程：（6-1）X方向的动量方程：（6-2）Y方向的动量方程：（6-3）能量方程：（6-4）5拉克斯-温德罗夫方法的基础是时间导数的泰勒展开式。任意选择一个流动参量，为明确起见，选择密度。时刻，同一网格点(i,j)处的密度可由泰勒级数给出：拉克斯-温德罗夫方法（6-5）6拉克斯-温德罗夫方法

应该从连续性方程（6-1）得到，将方程中的空间导数由二阶中心差分给出，即

也可用类似的方法得到，但需要更复杂的计算，将方程（6-1）对时间求导，得到（6-9）（6-10）7拉克斯-温德罗夫方法

等二阶混合导数，可由方程式（6-1）~式（6-4)对相应的空间自变量求导得到。例如，可由方程（6-2）对x求导得到（6-11）8拉克斯-温德罗夫方法式（6-11）中，右边所有的项还是由t时刻的二阶中心有限差分表示，即（6-12）9（6-10）中二阶混合导数都可以用类似的方法得到。一阶空间导数仍由二阶中心差分代替，例如：拉克斯-温德罗夫方法有了这些值，终于可以从式（6-10）中得到的值。的值已从式（6-9）得到，所以我们知道了式（6-5）右边三项在t时刻的值。这样就可以用（6-5）计算出t+△t时刻的密度值10图6-2给出了上述计算过程中用到的网格，图中显示了t时刻和t+△t时刻两个时间平面上的空间网格。拉克斯-温德罗夫方法11麦考马克方法还是讨论欧拉方程（6-1）~方程（6-4）的求解。还是先考虑密度，在麦考马克方法中，密度用下式得到（6-13）表示在t和t+△t时刻之间的平均值12时间导数的平均值可通过下面的预估-校正原理得到麦考马克方法预估步在连续性方程中，用向前差分代替方程右边的空间导数（6-17）13取泰勒级数的前两项来求密度的估计值麦考马克方法（6-18）校正步首先用向后差分代替连续性方程右边的空间导数，然后用、u和v的预估值进行计算，得到时间导数在t+△t时刻的预估值14麦考马克方法（6-21）密度时间导数的平均值是和的算术平均值，即（6-22）从式（6-13）中得到t+△t时刻密度最终的“校正”值，即（6-13）15麦考马克方法与拉克斯-温德罗夫方法对比：麦考马克方法在预估步中用向前差分在校正步中用向后差分，具有二阶精度，与拉克斯-温德罗夫方法具有同样的精度。但是麦考马克方法不像拉克斯-温德罗夫方法那样需要计算二阶时间导数，所以麦考马克方法更容易应用。麦考马克方法16粘性流动的控制方程是汇总在2.8.1小节中的纳维-斯托克斯方程。粘性流动对于定常流动，纳维-斯托克斯方程的数学性质更多的表现为椭圆形的，麦考马克方法和拉克斯-温德罗夫方法都不适合于椭圆形偏微分方程的求解。对于非定常纳维-斯托克斯方程具有抛物型和椭圆型混合的性质，麦考马克方法和拉克斯-温德罗夫方法这时都是适用的。17还用欧拉方程进行讨论。二维流动的适用于CFD计算的守恒型方程：守恒形式（6-23）显然，用麦考马克方法和拉克斯-温德罗夫方法，都可以计算U的分量、、、在各时间步的值。注意:方程（6－23）中未知函数是守恒变量，因此在每个时间步结束时，要用式（2－100）～式（2－104）各式计算出每一个原始变量。18举例：将麦考马克方法应用与图6－3所示的二维流动。空间推进假设流动是无粘的，则控制方程为欧拉方程。二维的情形简化为：（6-24）19对于亚声速流动，方程是椭圆型方程，麦考马克方法不适应。处于一个处处超声速的流动来说，方程（4－24）是双曲型方程。麦考马克方法适用。运用麦考马克方法，可由网格点（i，j＋1）、(i，j）和（i，j－1）处的已知值计算网格点（i＋1，j）处的流场变量。方程（6－24）中的向量F，它在网格点（i＋1，j）处的值可从下式求出：空间推进（6-24）20空间推进通过预估－校正法得到这个平均值预估步用向前差分替代对y的导数:（6-26）再由泰勒级数计算F在网格点（i＋1，j）处的预估值（6-27）21空间推进F代表了它的每一个分量的预估值（6-28）这些原始变量的预估值在校正步中将用来计算通量向量G22空间推进校正步将J和G的预估值带入方程（6－24），用向后差分计算在x＋△x处的预估值，即（6-29）平均值由算术平均给出（6-30）23

最后的校正值可从（6－25）得到（6-25）空间推进沿流向的推进方法与时间推进方法的两个值得注意的区别：采用守恒方程进行时间推进求解时，原始变量的计算很简单，守恒方程的空间推进解法中，原始变量的计算就复杂多了。至少对于显示求解而言，沿流向的推进一定要用守恒型方程。24松弛法及其在低速无粘流动中的应用松弛法是一种特别适用于求解椭圆型偏微分方程的有限差分法亚声速无粘流动由椭圆型偏微分方程控制，因此松弛法经常被用来求解亚声速的低速流动。举例：无粘不可压流体的二维无旋流动控制方程：（6-31）25将在图6-4所示的网格上数值求解方程（6-31）。松弛法及其在低速无粘流动中的应用26松弛法及其在低速无粘流动中的应用原则上，在每个内部网格点上都可以写出式（6-23），15个线性方程组成了包含15个未知数的方程组。使用克莱默法则是不切实际的，另一种更合理的直接方法是Gauss消去法。然而，最简单的方法是运用下面的松弛法。利用二阶导数的二阶中心差分式（4-12）和式（4-13）代替方程（6-31）中的偏导数，得：（6-32）27松弛法及其在低速无粘流动中的应用松弛法是一种迭代法，方程（6-32）中有四个量的值被看成已知的第n次迭代值，只有一个量被看成是未知的第n+1次迭代值。从方程（6-32）中解出，有：（6-33）再具体些，将方程（6-33）应用到图6-4中的网格点21，假设已经进行了n次迭代。为了进行n+1次，方程（6-33）变换为：28松弛法及其在低速无粘流动中的应用（6-34）未知量上次迭代中已经求出来的边界条件给定的有人建议将更新的值尽快地用到方程（6-33）右边。例如，从方程（6-34）中计算出后，转移到网格点22，由方程（6-33）得到：29松弛法及其在低速无粘流动中的应用（6-35）刚刚从方程（6-34）得到的在这种方式中，可沿着一条给定的水平线，从左到右或从上到下扫描，扫描的方向并不重要。逐个计算第n+1次迭代时的未知量（这种方法叫做Gauss-Seidel迭代法）。当所有网格点处的变得都小于一个预定的值时，就认为迭代收敛了。30逐次超松弛法可加快收敛的过程，这是一种外推算法。方程（6-33）的结果可以看成是的中间值，记做是，即松弛法及其在低速无粘流动中的应用（6-36）最后，用上一次迭代得到的和由方程（6-36）得到的，外推出的值：（6-37）松弛因子31松弛因子的值通常根据试算的经验确定。

>1，上面的过程叫做逐次超松弛法，对超松弛法通常在1<<2之间取值。

<1，叫做逐次低松弛法,当收敛过程表现为在某个值附近来回摆动时，通常会用到低松弛法。

松弛法及其在低速无粘流动中的应用326.6数值耗散、色散及人工粘性概念的导出数值耗散，人工粘性的影响33一维波动方程 （6-38）用时间的一阶向前差分和空间的一阶向后差分将它离散（6-39）1、数值耗散的概念34将方程（6-39）中的和展开成泰勒级数，即（6-40）

（6-41）35将（6-40）中（6-41）代入方程（6-39）中，整理得到

（6-43）

截断误差36为了用x的导数替换方程右边的时间导数，将方程进一步转化把式（6-43）对t求导得到（6-44）把式（6-43）对x求导并在两边同时乘以a得到（6-45）两式相减，只写出一阶项使得方程更紧密，得到

（6-47）37

方程（6-47）中给出了的表达式，但仍没有消除掉时间导数，转化还没有结束。将（6-47）对时间求导，得到方程（6-48）对（6-43）求x的二阶导数，两边乘以a平方，得到（6-49）两式相加，得到

（6-50）38

方程（6-50）中给出了的表达式，但（6-47）还有两个关于x和t的混合导数项需要处理，转化还没有结束。将（6-47）对x求导，整理后得到整理式（6-48），有将（6-50）代入到（6-52）中，得到

（6-53）

（6-52）

（6-51）39将（6-50），（6-51）和（6-53）代入方程（6-47）中，得到

（6-54）我们前面已经得到

（6-50）40将（6-54）和（6-50）代入方程（6-43）中，并且定义v=a△t/△x，得到修正方程的表达式

（6-56）41

推导修正方程（6-56）达到了两个目的。一，对差分方程精确解的意义作出了一种新的解释。 方程（6-56）是（6-39）的精确解，若把它作为（6-38）解就产生了截断误差。二，给出了与差分方程解的性质有关的信息，导出了数值耗散的概念。

42纳维-斯托克斯方程中有的项，它们表示了物理粘性对流动的耗散。方程（6-56）中的项是数值离散出来的，完全来源于数值过程，没有物理意义。数值解中出现的这一项及其类似项称为数值耗散，而这一项的系数，因其作用很像物理粘性，称为人工粘性。43图6-5数值耗散的影响2、数值耗散，人工粘性的影响a，t=0时刻的波形b，t>0时刻的波形

一维无粘流体中传播的波的数值耗散如下图44图6-5数值色散的影响a，t=0时刻的波形b，t>0时刻的波形数值色散是奇数阶导数（）的直接结果45

现在来研究一种特定形式的人工粘性，假设考虑的是非定常二维流，即

（6-57）

在时间推进法的每一步，都加一个小的人工粘性，如下

（6-58）46

假设正在使用麦考马克法，在预估步是根据t时刻的已知计算的;在校正步，式（6-58）右边的量是预估量，这样就得到了与对应的值，记作

（6-59）47

和的值要加到麦考马克法的计算过程中去。在方程（6-18）中加入预估步，方程变成

（6-60）

校正步用式（6-59）的到人工粘性项，而t+△t时刻的密度修正值可由方程（6-13）加上这个人工粘性得到，即

（6-58）

人工粘性的形式完全是经验性的。48

把人工粘性添加到麦考马克方法的计算过程中去。在方程（6-18）中加入这个人工粘性项。

（6-60）

校正步用式（6-59）的到人工粘性项，而t+△t时刻的密度修正值可由方程（6-13）加上这个人工粘性得到，即

（6-61）

人工粘性的形式完全是经验性的。49

添加人工粘性对解的精度有多大程度的影响，这个问题取决于流动问题本身的性质，这里给出一些结果。 研究的问题是绕后台阶的超声速粘性流，如图6-7所示。50

图6-8给出了用麦考马克方法计算得到的压力等值线，其中包括四个不同的等值线，每一个对应的Cx和Cy在0-0.3之间取值。可以发现随着Cx和Cy增加，也就是人工粘性的增加，流体的性质受到影响。51人工粘性对流场解的影响类似于物理粘性的影响，和加大μ一样，通过增加人工粘性，可以使激波抹平、变宽。人工粘性的作用：使数值解稳定、光滑。526.7交替方向隐式（ADI）方法问题的提出交替方法隐式（ADI）方法

53

考虑一个非定常二维热传导方程 1、问题的提出

（6-62）

将4.4节中的克兰克-尼克尔森方法平行推广到这里，方程有限差分形式为 54

（6-63）55方程（4-40）可以转化为三对角形式的方程组（4-42）方程（6-63）却含有五个未知量，不能使用托马斯算法如果能发展一种格式，使得方程（6-62）用三对角形式就可以求解，将具有明显的优势。这种格式就是：交替方向隐式（ADI）方法

在xy空间中，方程（6-63）是一维问题的方程（4-40）对应形式，但是56

分成两步求解，第一步的时间步长为t+△t/2

，只对x的导数采用隐身处理 1、交替方向隐式（ADI）方法

（6-64）57

方程（6-64）可化简为三对角形式

（6-65）

对每一个固定的j，对所有的i写出方程（6-65），联立形成一个方程组。用追赶法得到的解（j固定）。58

如图，对于固定的j，我们在x方向扫描，得到n+1/2时刻的值。 假如在x方向有N个网格点，那么就从i=1扫描到N。扫描一次用一次托马斯算法。

y方向有M个网格点，这个过程就要重复M次。

59

第二步是利用t+△t/2时刻的已知值求解t+△t时刻的解

，仍用中心差分替代方程（6-62）中的空间导数，但这次只对只对y的导数采用隐式处理，即

（6-66）60

方程（6-66）可化简为三对角形式

（6-67）

对每一个固定的i，对所有的j写出方程（6-67），联立形成一个方程组。用托马斯算法得到的解（i固定）。61

第二步是同样在y方向扫描。 在这两个过程结束后，未知函数T在时间方向推进了步长△t。 在这个格式中，第一步差分方程的x方向是隐式的，第二步y方向是隐式的，所以称为交替方向隐式方法。626.8压力修正法及其在不可压粘性流动中的应用关于不可压Stokes方程的注释交错网格的应用压力修正法的基本原理压力修正公式数值方法：SIMPLE算法压力修正法的边界条件63

对于不可压Stokes方程，假设整个流动过程中的μ等于常数，方程（2-50a-c）变为 1、关于不可压Stokes方程的注释

（6-69）

（6-70）

（6-71）64

下面对方程（6-69-71）进行简化。在不可压流动中，由ρ=常数，有

将上式对x求导

（6-72）

（6-73）

（6-74）65

等式两边同时加上，并乘以μ，得

用此式代替方程（6-69）右边的第二项，展开方程（6-69）其他项并消去相应的项，得到不可压粘性流动x方向动量方程的一种简明的形式

（6-75）

（6-74）66

方程式（6-70）和式（6-71）可用同样的方式化简。最终得到的方程组就是不可压纳维-斯托克斯方程，汇总如下。（6-77）（6-78）（6-79）（6-80）连续性方程X方向动量方程Y方向动量方程Z方向动量方程67四个方程四个未知数，表明求解不可压流动的速度场和压力场，只需用到连续性方程和动量方程。CFD中不可压纳维-斯托克斯方程的解法通常不同于可压纳维-斯托克斯方程的解法。压力修正方程克服了这个困难，并以成功应用于可压流动。68

方程（6-77）给出了不可压流体的连续性方程。对于二维情形，有 2、交错网格的应用

（6-82）

相应的中心差分格式为

（6-83）69

这个差分方程给出的速度场，可能会出现棋盘式分布。 若将这些数值代入方程中去，则在每个网格点上，方程中的两项都为0。 这样的离散速度满足连续性方程，但是对于实际流动毫无意义。 70

在动量方程（6-77-80）中也有同样的问题。考虑图中的棋盘式的二维离散压力分布，并考虑压力梯度的中心差分公式

（6-84a）

（6-84a）71

如图，斯托克斯方程感觉不到这样的压力变化，而是把这种压力分布误认为是x方向和y方向的均匀分布。 两种修补方法，一是改用迎风差分，另一种是在交错网格上使用中心差分。72

如图，斯托克斯方程感觉不到这样的压力变化，而是把这种压力分布误认为是x方向和y方向的均匀分布。 两种修补方法，一是改用迎风差分，另一种是在交错网格上使用中心差分。73

如图，是一个交错网格，在标号为（i-1，j）、（i，j）。（i+1，j）、（i，j-1）等网格点上计算压力，在标号为（i-1/2，j）、（i+1/2，j）、（i，j+1/2）、（i，j-1/2）的点上计算速度。 这里的关键是在不同的网格点上计算压力和速度。74

交错网格的优点在于：计算时，压力梯度的中心差 分是基于两个相邻的压力网格点，这就消除了出现图

6-14所示的棋盘式的压力分布的可能性。 在网格点（i，j）上，连续性方程的中心差分表达式变成

（6-85）

同样地，由于方程（6-85）是基于相邻的速度网格点，所以图6-13所示的棋盘式速度分布也就不可能出现了。75压力修正法本质上是一种迭代法，思路如下给定压力的初始近似。用求解出，，。将，，代入连续性方程，构造压力修正量，使速度场满足连续性方程，其中用方程（6-86）左边的p值作为新的，回到步骤2，重复这个过程直到速度场满足连续性方程为止。 3、压力修正法的基本原理

（6-86）

（6-87a）

（6-87b）

（6-87c）76

对于不可压流体，把x，y方向的动量方程（6-78）和（6-79）写成守恒形式，就是4、压力修正公式

（6-88）

（6-89）77

考虑图6-16所示的交错网格，压力在实心圆点上计算，速度在空心圆点上计算，图中的（i+1/2，j）点对方程（6-88）进行中心差分。78

中心差分要用到阴影上边a点和下边b点处v的平均值，我们用相邻两点的线性插值来定义这些平均值，即

（6-90a）

（6-90b）在a点在b点79

以（i+1/2，j）点为中心，方程（6-88）的差分方程为

（6-92）式中80

方程（6-92）是x方向动量方程的差分方程。 同样的，y方向的动量方程的差分方程也可以用同样的方式得到。81

同理，方程（6-89）的差分方程为

（6-92）式中82

开始迭代，方程（6-92）和（6-93）分别为

（6-94）

（6-95）83

从式（6-92）中减去式（6-94），得

（6-96）

式中84

同理，从式（6-93）中减去式（6-95），得

（6-97）

式（6-96）和式（6-97）就是用压力和速度修正量表示的x方向和y方向的动量方程。 我们设计的的计算公式只要满足两点 （1）的计算公式最终要得到合理的、收敛的解。 （2）收敛后，的计算公式必须还能还原成物理上真正的连线性方程。85

下面我们开始构成压力修正方程

（6-100）设式（6-96）和式（6-97）中的都为零根据的定义，，于是得到

（6-101）

同理86

回到连续性方程

（6-102）

在网格点（i，j）处写出相应的中心差分方程，为

把式（6-100）和式（6-101）去掉括号外的上标，代入方程（6-102）中，整理得到87

（6-104）

式（6-104）就是压力修正方程，它具有椭圆型的性质。因此用6.5节描述的松弛法数值求解方程（6-104），就可以得到p’。88现在总结一下压力修正法的计算步骤按照交错网格，给定压力网格点的，和相应速度网格点的和。从方程（6-94）和（6-95）中解出，。将其代入