计算流体力学基础简介

计算流体力学基础简介本硕111班王鹏主要内容：

流体力学的控制方程组适合CFD使用的控制方程组偏微分方程的离散化——有限差分计算流体力学的基本方法流体力学的控制方程组一.粘性流动的N-S方程

1.连续性方程

非守恒形式守恒形式

2.动量方程

非守恒形式

流体力学的控制方程组守恒形式3.能量方程非守恒形式

流体力学的控制方程组

守恒形式

流体力学的控制方程组一.无粘流Euler方程

1.连续性方程

非守恒形式守恒形式

2.动量方程

非守恒形式

流体力学的控制方程组守恒形式3.能量方程非守恒形式守恒形式

流体力学的控制方程组对于牛顿流体，斯托克斯在1845年得到，在上面的方程组中：

根据傅里叶热传导定律，热传导产生的热流与当地的温度梯度成正比:其中k为热导率。适合CFD使用的控制方程组

守恒形式的连续性方程、动量方程和能量方程可以使用同一个通用方程来表达，这有助于计算程序的简化和结构组织，为算法设计和编程提供了方便。

守恒形式的通用方程

（1）

如果将上式中的U、F、G、H和J看成列向量，则上式就可代表整个守恒形式的控制方程。这些列向量为：

适合CFD使用的控制方程组

（1）式中的F、G和H称为通量项(通量向量），J代表源项，U称为解向量。对于一般的时间推进方法，我们可将（1）式变换一下以方便求解

适合CFD使用的控制方程组

针对Euler方程也可进行类似分析。偏微分方程的离散化——有限差分

偏微分方程的解析解是封闭形式的表达式，它给出了未知函数在区域内的连续变化。偏微分方程的离散化是用另外一个类似的表达式来近似这些偏微分方程，但是这个近似表达式只在去域内有限多个离散点或控制体上规定了取值。

利用泰勒级数展开可以推导出各阶导数的有限差分的一般形式。（2）（3）利用（2）、（3）式可以得到以下一些差分公式：偏微分方程的离散化——有限差分一阶导数的一阶向前差分一阶导数的二阶中心差分一阶导数的一阶向后差分二阶导数的二阶中心差分一阶混合导数的二阶中心差分……..偏微分方程的离散化——有限差分利用多项式也可以实现有限差分。例如，如右图所示的边界上，u可以表示成多项式：

（4）（5）（6）（7）偏微分方程的离散化——有限差分解出b,再对（4）式求导，边界上y=0,代入上式有，最终得到上式就是单侧差分表达式，之所以被称为单侧，是因为它只用到了边界一侧的网格点上的信息。计算流体力学的基本方法

计算流体力学包括很多方法，其中较基础的有拉克斯-温德罗夫方法、麦考马克方法、空间推进方法、松弛法、交替方向隐式方法、压力修正法等。1.拉克斯-温德罗夫方法

以非定常二维无粘流为例，控制方程如下

（8）（9）（10）（11）计算流体力学的基本方法

（12）（13）（14）（15）

计算流体力学的基本方法

（17）

再将方程（8）对时间求导：

上式中出现二阶混合导数，可由方程式（8）到（11）对相应的空间自变量求导得到。例如（16）

将右边所有的项还是由t时刻的二阶中心差分表示，即计算流体力学的基本方法

计算流体力学的基本方法

二.麦考马克（Macormack）方法

麦考马克方法是从拉克斯-温德罗夫方法中变化过来的一种方法，但更简单，在1969年一经提出，就成为解决流动问题最流行的显示有限差分方法，一直流行了15年，现在它已被更先进的方法取代了。

（18）计算流体力学的基本方法

预估步：在连续方程（8）中，用向前差分代替方程右边的空间导数：

（19）同理可求出u、v和e的预估值，即计算流体力学的基本方法（20）

式（18）中密度的时间导数的平均值

计算流体力学的基本方法（21）

松弛法是一种适合于求解椭圆型偏微分方程的迭代法。例如，亚声速无粘流动由椭圆型偏微分方程控制松弛法经常被用来求解亚声速的低速流动。松弛法可以是显式的也可以是隐式的，这里介绍显式松弛法。

我们将在右图所示的网格上数值求解方程（21），利用二阶导数的二阶中心差分式代替上式的偏导数，得

上式中，上标n和n+1表示迭代次数。（22）（23）

三.松弛法计算流体力学的基本方法（24）

在对所有内部网格点都使用了方程（23）后，就完成了第一次迭代，n=1。接着进行第二次迭代，n=2。这个迭代过程反复不断，直至收敛到解。（为什么一定是朝解的方向收敛？）

再具体些，将方程（23）应用到网格点21，进行第n+1次迭代，有

（25）计算流体力学的基本方法

计算流体力学的基本方法（26）

四.交替方向隐式（ADI）方法讲这个方法之前，有必要先介绍一下显示方法和隐式方法。所有的有限差分方法都可以归结为这两种不同的通用方法——显示方法或隐式方法。显示方法。以一位热传导方程为例等式左边用向前差分，等是右边用二阶中心差分，即整理此式后可写成方程（26）是抛物线型偏微分方程，可推进求解。此时，推进变量是时间t，为更简明起见，考虑右图所示的有限差分网格。

假设在第n个时间层上的每个网格点上温度T是已知的，时间推进就意味着第n+1个时间层每个网格点上的T值都要用第n个时间层上的已知量计算出来。（28）（27）计算流体力学的基本方法（29）

例如，在网格点2处，式（28）将写成

隐示方法

我们可以再大胆些，把式（27）右边的空间差分写成第n个时间层的量与第n+1个时间层上的量的平均值，即

计算流体力学的基本方法（30）

考虑右图所示的有限差分网格，以七点空间网格点作为具体例子。整理式（29），将未知量放在等式左边，已知量放到等式右边，得

于是，式（30）可写成

（31）计算流体力学的基本方法

把方程（31）~（35）写成矩阵形式：（32）（33）（34）（35）

这是一个三对角矩阵，可用托马斯算法（国内称为追赶法）求解。解出之后即得到第n+1个时间层上的温度T值。重复以上过程，即可得到最终的数值解。计算流体力学的基本方法（37）

现在，介绍交替方向隐式（AlternatingDirectionImplicit,ADI）方法。以二维热传导方程为例，

对上式运用克兰克-尼克尔森方法，有

（36）（38）计算流体力学的基本方法（39）

方程（38）可化为三对角形式

式中

计算流体力学的基本方法（40）

方程（40）可化为三对角形式

式中

（41）计算流体力学的基本方法（46）五.压力修正法

压力修正法本质上也是一种迭代法，思路如下：

1）迭代开始时，先给定流场压力的初始近似值p*。

2）用p*的值从动量方程中求解出u,v,w.因为这些速度都与p*有关，所以用u*,v*,w*表示它们。

3）

u*,v*,w*不一定满足连续方程。因此，要用连续方程构造压力的修正量p’，加到p*上，使速度满足连续方程。设修正后的压力为

相应的速度修正量u’,v’,w’可以从p’得到，使得

压力修正法的基本原理（42）（43）（44）（45）计算流体力学的基本方法（47）

此时我们要做到就是求出压力修正量p’，为简单起见，考虑二维流动，并忽略体积力。对不可压流动，x方向和y方向的动量方程分别是式（43）和（44），这些方程都是非守恒型的。若写成守恒形式就是4）用方程（46）左边的p作为新的p*，回到步骤2）。重复这个过程，直到速度场满足连续方程为止。压力修正公式（48）计算流体力学的基本方法（49）

考虑如右图所示的交错网格，压力在实心圆点上计算，速度在空心圆点上计算。在右图中的点（i+1/2,j)点对x方向动量方程（47）进行中心差分（右图中阴影部分相当于有限体积法中的控制体）。中心差分要用到阴影区上边a点和下边b点处的v的平均值，我们用相邻两点的线性插值来定义这些平均值，即

以（i+1/2,j)点为中心，方程（47）的差分方程为

即式中计算流体力学的基本方法（50）

y方向动量方程（48）的中心差分方程也可以用相同的方式得到。中心差分要用到阴影区左边的c点和右边的d点处u的平均值，定义为

以（i,j+1/2)点为中心，对方程（48）进行差分，有式中计算流体力学的基本方法（51）根据前面给出的步骤，在迭代开始时，给定p=p*。此时，式（49）和（50）分别为（52）从式（49）中减去式（51），得式中从式（50）中减去式（52），得（53）（54）计算流体力学的基本方法（55）（56）式中

式（53）和式（54），就是用压力修正量p’和速度修正量u’,v’表示的动量方程。

则，式（55）可表达为计算流体力学的基本方法（58）（57）

则，式（56）可表达为

回到二维连续方程

在网格点（i,j)处写出相应的差分方程，即（59）

去掉式（57）（58）的上标，然后代入式（59），得到

整理上式，得到计算流体力学的基本方法（60）式中方程（60）就是压力修正公式，它具有椭圆性质。因此用前面讲过的松弛法