计算流体力学数值方法

1、3-1 计算流体力学计算流体力学第三章第三章 数值方法数值方法3-2 计算流体力学计算流体力学主要内容主要内容u空间离散技术空间离散技术 标量输运方程标量输运方程 动量方程动量方程u时间离散技术时间离散技术3-3 计算流体力学计算流体力学1 标量输运方程标量输运方程u有限体积法有限体积法u一维对流扩散方程一维对流扩散方程u扩散项的离散扩散项的离散u源项的离散源项的离散u代数方程的组装代数方程的组装u二维和三维问题二维和三维问题u对流项离散基础对流项离散基础u离散特性离散特性3-4 计算流体力学计算流体力学u对流项高级离散方法对流项高级离散方法 u高阶对流方法的实现高阶对流方法的实现u曲线网格曲

2、线网格u边界条件边界条件u代数方程的求解代数方程的求解u小结小结3-5 计算流体力学计算流体力学对于任意控制体内一个守恒物理量对于任意控制体内一个守恒物理量: 控制体内物理量的控制体内物理量的变化率变化率+经过控制经过控制体边界面的体边界面的净通量净通量控制体积内的控制体积内的源源(汇汇) 经过控制体边界面的总的通量由经过控制体边界面的总的通量由对流对流(随流动的迁移随流动的迁移)和和扩散扩散(由随机分子运动或由随机分子运动或湍流运动导致的净输运湍流运动导致的净输运)两部分组成。如果两部分组成。如果用用 表示单位质量流体的守恒量，则通用标表示单位质量流体的守恒量，则通用标量输运方程量输运方程(

3、或称对流扩散方程或称对流扩散方程)可表示可表示为：为：3-6 计算流体力学计算流体力学1，有限体积法直接对上式进行离散，有限体积法直接对上式进行离散2，本章只考虑稳态问题，即上式左边第一项为，本章只考虑稳态问题，即上式左边第一项为零零SVnACVdtdfaces)(）（AuCn3-7 计算流体力学计算流体力学u 有限体积法有限体积法(FVM) (1) 定义流场求解域几何形定义流场求解域几何形状状 (2) 将求解域划分为计算网将求解域划分为计算网格，即一组互不重叠的有限格，即一组互不重叠的有限体或单元。体或单元。 (3) 基于上述划分的单元对基于上述划分的单元对积分方程进行离散，即用节积分方程进

4、行离散，即用节点值来近似。点值来近似。 (4) 对得到的离散方程进行对得到的离散方程进行数值求解。数值求解。3-8 计算流体力学计算流体力学 计算网格可以是计算网格可以是结构网格结构网格或或非结构网格非结构网格，笛笛卡儿网格卡儿网格或或非笛卡儿网格非笛卡儿网格。 最常用的网格形式包括最常用的网格形式包括基于单元中心基于单元中心的存储的存储方式和方式和基于单元顶点基于单元顶点的存储方式两种。在后面的的存储方式两种。在后面的讨论中将会看到，不是所有的变量都必须存储在讨论中将会看到，不是所有的变量都必须存储在同样的位置。同样的位置。 3-9 计算流体力学计算流体力学 为了简单起见，本课程将只考虑为了

5、简单起见，本课程将只考虑结构化笛卡结构化笛卡儿网格儿网格，和采用，和采用基于单元中心基于单元中心的存储方式，即单的存储方式，即单元变量元变量 存储在单元中心节点位置上。存储在单元中心节点位置上。 右图为一典型三维控制右图为一典型三维控制体积，其中体积，其中P点为单元中心点为单元中心存储节点。相对于中心节点，存储节点。相对于中心节点，沿坐标方向通常表示为沿坐标方向通常表示为west, east, south, north, bottom, top. . 其中，小写字母其中，小写字母w, e, s, n, b, t 表示单元面，大写字母表示单元面，大写字母W,E, S, N, B, T表示中心节点

6、的表示中心节点的相邻节点。相邻节点。3-10 计算流体力学计算流体力学 于是单元各面的面积表示为于是单元各面的面积表示为Aw, Ae, As, An Ab, At 。体积为。体积为V。对于二维问题，可以视为单。对于二维问题，可以视为单位厚度为位厚度为1的一层单元的一层单元( ) 对于结构网格，可以交替使用对于结构网格，可以交替使用ijk下标表示单下标表示单元节点，譬如，元节点，譬如， 3-11 计算流体力学计算流体力学u 一维对流扩散方程一维对流扩散方程 我们首先讨论一维稳态对流扩散方程主要基我们首先讨论一维稳态对流扩散方程主要基于下面的考虑：于下面的考虑： (1) 它使问题分析大大简化它使问

7、题分析大大简化 (2) 离散方程可以进行手算。离散方程可以进行手算。 (3) 尽管只是一维的，但要扩展到二维或三维尽管只是一维的，但要扩展到二维或三维是非常直接的是非常直接的 (4) 实际上，通量实际上，通量( (对流和扩散对流和扩散) )的离散一般是的离散一般是沿坐标方向进行的，即分别沿沿坐标方向进行的，即分别沿i,j,k线进行。线进行。 (5) 有许多重要的理论问题是一维的。有许多重要的理论问题是一维的。3-12 计算流体力学计算流体力学 对于如图所示的一维控制体，物理量的守恒对于如图所示的一维控制体，物理量的守恒可表述为如下关系式：可表述为如下关系式： 通量通量 e (fluxe)- 通

8、量通量 w ( (fluxw) ) 源源( (source) ) 这里的通量是指穿过单这里的通量是指穿过单元表面的输运率。元表面的输运率。 如果如果 表示单位质量表示单位质量的输运量，则总的通量为的输运量，则总的通量为对流通量和扩散通量之和，对流通量和扩散通量之和，其中：其中： 对流通量对流通量 扩散通量扩散通量 3-13 计算流体力学计算流体力学 于是，于是， 的对流扩散方程为：的对流扩散方程为： 其中，其中，S表示单位长度的源。将等式两边除以表示单位长度的源。将等式两边除以 ，并取极限，并取极限 ，可得到相应的微分方程：，可得到相应的微分方程：注意：注意：1 1，这里面积，这里面积A A是

9、表示可变截面的准一维问题是表示可变截面的准一维问题，对于真正的一维问题，只需令，对于真正的一维问题，只需令A=1,A=1,即即 2 2，这里假设，这里假设 和和S均为常数，但在一般的均为常数，但在一般的CFDCFD问题中，问题中，u本身也是问题的解变量。本身也是问题的解变量。3-14 计算流体力学计算流体力学u 例子例子 1 1，纯扩散问题，纯扩散问题 如图所示的隔热棒，长度如图所示的隔热棒，长度1m，截面为，截面为1cmx1cm的方的方形截面，棒的两端为固定形截面，棒的两端为固定温度，分别为温度，分别为100度和度和500度。穿过任意截面度。穿过任意截面A的热通的热通量由下式给定：量由下式给

10、定： 其中，热传导系数其中，热传导系数u 将棒划分为将棒划分为5个控制体，并用有限体积分析沿棒个控制体，并用有限体积分析沿棒的温度分布的温度分布u 写出沿棒温度分布的微分方程写出沿棒温度分布的微分方程u 求出微分方程的解析解，并与求出微分方程的解析解，并与(a)(a)的解进行比较的解进行比较3-15 计算流体力学计算流体力学 这个问题是求解一维对流扩散方程：这个问题是求解一维对流扩散方程： 上面的方程也可写成如下的积分方程：上面的方程也可写成如下的积分方程： 由于本例只考虑扩散，即没有对流和源项：由于本例只考虑扩散，即没有对流和源项： 3-16 计算流体力学计算流体力学 对于本例的温度对于本例

11、的温度T的纯扩撒问题，最终有如下的纯扩撒问题，最终有如下的微分方程和积分方程：的微分方程和积分方程： 00)(ewdxdTkAdxdTkAdxd3-17 计算流体力学计算流体力学 非边界面上通量的计算：非边界面上通量的计算： 3-18 计算流体力学计算流体力学 边界面上通量的计算：边界面上通量的计算： 3-19 计算流体力学计算流体力学02/00002/52454534342323121211xTTxTTxTTxTTxTTxTTxTTxTTxTTxTTBB3-20 计算流体力学计算流体力学1000302020220035454343232121TTTTTTTTTTTTT460380300220

12、14014311TTTTT3-21 计算流体力学计算流体力学 解析解：解析解： 100400500) 1(100)0(00)(0122xTxTxTcxcTdxTddxdTkAdxd3-22 计算流体力学计算流体力学u 控制方程扩散项的离散控制方程扩散项的离散 梯度扩散项的离散几乎梯度扩散项的离散几乎总是采用中心差分格式：总是采用中心差分格式： 提示：提示： (1) 有限体积法中，通量是在面上计算的，而有限体积法中，通量是在面上计算的，而不是在节点处。不是在节点处。 (2) (2) 上述对梯度扩散项的近似在空间上具有二上述对梯度扩散项的近似在空间上具有二阶精度，后面将会给出证明。阶精度，后面将会

13、给出证明。3-23 计算流体力学计算流体力学 (3) 如果扩散系数也是如果扩散系数也是 变量的话，它在单元面变量的话，它在单元面 上的值必须通过插值得到。上的值必须通过插值得到。 (4) 采用中心差分格式意采用中心差分格式意 味着在单元面两边的节点上权重相同。这味着在单元面两边的节点上权重相同。这 和扩撒的物理意义是一致的，因为扩散在和扩撒的物理意义是一致的，因为扩散在 所有方向作用相同。这后面将要讨论的对流所有方向作用相同。这后面将要讨论的对流 是不一样的，对流是具有明确的方向性的。是不一样的，对流是具有明确的方向性的。 3-24 计算流体力学计算流体力学 (5) (5) 单元西面也应有相似

14、的单元西面也应有相似的表达式，这是守恒定律所要表达式，这是守恒定律所要求的。即从一个单元流出的求的。即从一个单元流出的通量一定等于流入相邻单元通量一定等于流入相邻单元的通量。这是有限体积法优越有限差分和有限的通量。这是有限体积法优越有限差分和有限元法的地方。元法的地方。 在有限体积法中，常常定义一个扩散输运在有限体积法中，常常定义一个扩散输运系数系数D D，譬如，对于单元东面有：，譬如，对于单元东面有：3-25 计算流体力学计算流体力学u 控制方程源项的离散控制方程源项的离散 当源与流体的量成比例的情况下，每个单当源与流体的量成比例的情况下，每个单元的总的源强度可简单的表示为：元的总的源强度可

15、简单的表示为： 单位体积的源单位体积的源x单元体积单元体积SxV 其中，其中，S表示平均源密度，一般取单元中心表示平均源密度，一般取单元中心节点处的值。节点处的值。 一般来说，一般来说，S常常是与变量常常是与变量 相关的，譬相关的，譬如，如果变量为温度的话，牛顿冷却定律认为如，如果变量为温度的话，牛顿冷却定律认为：单位时间热量的损失与相对与环境的温差成：单位时间热量的损失与相对与环境的温差成比例：比例： 一般将源项分解为与解相关和无关两部分一般将源项分解为与解相关和无关两部分：3-26 计算流体力学计算流体力学u 代数方程的组装代数方程的组装 当速度当速度u=0时，定常扩散时，定常扩散问题在每

16、个单元上的离散为：问题在每个单元上的离散为： 经过整理后：经过整理后： 或表示成如下标量输运方程规范形式：或表示成如下标量输运方程规范形式：3-27 计算流体力学计算流体力学 上述方程为离散后的标量输运方程规范形式上述方程为离散后的标量输运方程规范形式，每个求解变量的输运方程都有相同的形式，只，每个求解变量的输运方程都有相同的形式，只是各自的矩阵系数不同。对于纯扩散问题，矩阵是各自的矩阵系数不同。对于纯扩散问题，矩阵系数为：系数为： 3-28 计算流体力学计算流体力学 上述方程为只是一个控制体的离散方程，如上述方程为只是一个控制体的离散方程，如果将所有控制体单元中心节点处的变量值果将所有控制体

17、单元中心节点处的变量值 组装组装到一个向量中，就会得到如下形式的代数方程：到一个向量中，就会得到如下形式的代数方程： 求解这种代数方程可采用三对角矩阵算法求解这种代数方程可采用三对角矩阵算法3-29 计算流体力学计算流体力学u 二维和三维的扩展二维和三维的扩展 前面讨论的虽然只是一维前面讨论的虽然只是一维问题，但要推广到二维和三维问题，但要推广到二维和三维是非常直接的。是非常直接的。 对于多维问题，经过单元对于多维问题，经过单元面的净通量只需将所有面上的面的净通量只需将所有面上的通量求代数和即可！譬如，对于三维控制体，通量求代数和即可！譬如，对于三维控制体，其守恒方程可表示成：其守恒方程可表示

18、成： 3-30 计算流体力学计算流体力学 对于多维问题，离散后的代数方程仍然具有对于多维问题，离散后的代数方程仍然具有如下形式：如下形式：上述方程只是对于单个上述方程只是对于单个控制体单元的，而对所控制体单元的，而对所有单元进行组装便得到有单元进行组装便得到一组节点值向量一组节点值向量 的的矩阵方程。对于二维问矩阵方程。对于二维问题，最终得到的矩阵具有如图所示的带状形式题，最终得到的矩阵具有如图所示的带状形式。三维问题将增加非对角元素！。三维问题将增加非对角元素！3-31 计算流体力学计算流体力学u 对流项的离散对流项的离散(I) 纯扩散问题纯扩散问题( (u=0) )只是一只是一种理想情况，

19、对于典型的工种理想情况，对于典型的工程问题，对流项远远超过扩程问题，对流项远远超过扩散项，因为雷诺数非常大，即惯性力和粘性力散项，因为雷诺数非常大，即惯性力和粘性力的比值很大。的比值很大。 定常的对流扩散方程方程：定常的对流扩散方程方程： 3-32 计算流体力学计算流体力学 为了简化书写，将质量通量用为了简化书写，将质量通量用C表示，譬如表示，譬如，单元，单元e面的质量通量表示为：面的质量通量表示为： 扩散项和源项的离散采用与前面相同的方扩散项和源项的离散采用与前面相同的方法，于是有：法，于是有：现在的问题是如何离散对流项？现在的问题是如何离散对流项？也就是怎样利也就是怎样利用相邻节点值来确定

20、单元面上的变量值。这类用相邻节点值来确定单元面上的变量值。这类方法就称作对流方案。方法就称作对流方案。 3-33 计算流体力学计算流体力学u 中心差分法中心差分法 离散对流项的最明然离散对流项的最明然的方法就是线性插值，即的方法就是线性插值，即中心差分法：中心差分法： 于是将上述公式代入对流扩散方程后，得于是将上述公式代入对流扩散方程后，得到关于单个控制体单元如下形式的离散方程：到关于单个控制体单元如下形式的离散方程： 3-34 计算流体力学计算流体力学 上式中的求和是对中心点的所有相邻节点进行上式中的求和是对中心点的所有相邻节点进行的。各个系数表达式如下：的。各个系数表达式如下： 上式中，根

21、据质量守恒有上式中，根据质量守恒有 3-35 计算流体力学计算流体力学 为了说明中心差分格式是否适合离散对流项，为了说明中心差分格式是否适合离散对流项，下面给出了定常对流扩散问题下面给出了定常对流扩散问题( (无源项无源项) )在下面在下面两种情况下的解的图示：两种情况下的解的图示： Pe=1/2Pe=1/2 pe=4 pe=4其中其中PePe是是PecletPeclet的缩写，称作的缩写，称作PePe数。数。 3-36 计算流体力学计算流体力学 上图中红圈表示计算值，实线为解析解。第一上图中红圈表示计算值，实线为解析解。第一重情况吻合很好，但第二种情况却有显著的摆动重情况吻合很好，但第二种情

22、况却有显著的摆动，这是为什么呢？这是为什么呢？ 3-37 计算流体力学计算流体力学 从从数学数学的角度来说，当的角度来说，当Pe数大于数大于2 2，系数矩阵的一，系数矩阵的一个元素将是负值，即相应于该个元素将是负值，即相应于该元素的节点值的增加会使中心元素的节点值的增加会使中心节点值减小！节点值减小！ 从从物理物理的角度来看，是因的角度来看，是因为对流过程是有方向性的。而为对流过程是有方向性的。而输运属性仅仅在流动方向。但输运属性仅仅在流动方向。但中心差分格式没有反映处方向性，而是在上风中心差分格式没有反映处方向性，而是在上风(upwind)节点和下风节点和下风(downwind)节点上赋予了

23、相节点上赋予了相同的权重。同的权重。 3-38 计算流体力学计算流体力学u 上风差分法上风差分法 在上风差分格式种，在上风差分格式种，单元面上变量的值用其单元面上变量的值用其上风节点的值近似。譬上风节点的值近似。譬如，对于单元东如，对于单元东( (e) )面面: : 于是将上述公式可归结为：于是将上述公式可归结为： 3-39 计算流体力学计算流体力学 采用上风差分格式后，定常对流扩散问题采用上风差分格式后，定常对流扩散问题的规范离散格式为：的规范离散格式为： 其中：其中： 3-40 计算流体力学计算流体力学 当将上风差分法应用于前述与中心差分同样的当将上风差分法应用于前述与中心差分同样的问题后

24、可观察到：问题后可观察到： (1)(1) 当当 所得结果精度不如中心差分法，这是意料之所得结果精度不如中心差分法，这是意料之中的。中的。 (2) (2) 当当Pe=4Pe=4时，得到的解不是很精确，但不再时，得到的解不是很精确，但不再出现明显的摆动！出现明显的摆动！ 于是，便出现了精度和有界性于是，便出现了精度和有界性( (无摆动无摆动) )之间之间如何取舍的问题？下面将介绍一些离散方程的属如何取舍的问题？下面将介绍一些离散方程的属性。性。3-41 计算流体力学计算流体力学u 离散属性离散属性 1 1，一致性，一致性 一致性是指当网格间距趋向于一致性是指当网格间距趋向于0时，离散时，离散方程等

25、价于连续方程。譬如，方程等价于连续方程。譬如， 就是就是 的一致近似。的一致近似。 2，守恒性，守恒性 即从一个单元流出的即从一个单元流出的通量一定等于流入相邻单元通量一定等于流入相邻单元的通量。这是有限体积法优的通量。这是有限体积法优越有限差分和有限元法的地方。越有限差分和有限元法的地方。3-42 计算流体力学计算流体力学 3 3，输运性，输运性 方向性的影响反应在对流项的离散方法中方向性的影响反应在对流项的离散方法中。实际上，就是在上风节点赋予较高的权重。实际上，就是在上风节点赋予较高的权重。 4，有界性，有界性 在无源的对流扩散问题中的解一定介在无源的对流扩散问题中的解一定介于周围节点流

26、动变量最大值和最小值之间。于周围节点流动变量最大值和最小值之间。 5, 稳定性稳定性 稳定性决定了求解过程能否得到最终的稳定性决定了求解过程能否得到最终的解。它不涉及解的精度。稳定即意味着误差不解。它不涉及解的精度。稳定即意味着误差不随求解过程增加。随求解过程增加。3-43 计算流体力学计算流体力学 6 6，阶，阶( (解的精度解的精度) ) 阶阶是指随着空间网格尺寸的减小，误差衰是指随着空间网格尺寸的减小，误差衰减的速度。譬如，对于均匀空间网格减的速度。譬如，对于均匀空间网格 ，在某，在某种数值方法中，如果当空间网格尺寸趋向于种数值方法中，如果当空间网格尺寸趋向于0，而误差与空间网格尺寸的，

27、而误差与空间网格尺寸的n次方成次方成 正比正比，则称这种数值方法是，则称这种数值方法是n阶阶的。的。 某种数值方法的阶可以通过泰勒级数展开某种数值方法的阶可以通过泰勒级数展开式中的首阶误差项确定式中的首阶误差项确定(后面将讨论后面将讨论) 。 注意：此处的误差是指展开式中的理论截注意：此处的误差是指展开式中的理论截断误差，而不考虑计算机的舍入误差断误差，而不考虑计算机的舍入误差(即计算即计算机只能存储有限位数的数机只能存储有限位数的数)3-44 计算流体力学计算流体力学 高阶精度可通过采用更多的节点值近似来获高阶精度可通过采用更多的节点值近似来获得。一个节点允许的最高精度为得。一个节点允许的最

28、高精度为1 1阶，两个节点允阶，两个节点允许的最高精度为许的最高精度为2 2阶，依此类推。阶，依此类推。 理论上讲，某种数值方法的精度越高，随着理论上讲，某种数值方法的精度越高，随着网格的加密，误差减小的就越大网格的加密，误差减小的就越大 。也就是说，采。也就是说，采用高精度的数值方法，只需较少的网格数即可获用高精度的数值方法，只需较少的网格数即可获得要求的精度。得要求的精度。 但是，高阶精度的方法常常需要更多的计算但是，高阶精度的方法常常需要更多的计算时间，而且常常会导致解的有界性问题。时间，而且常常会导致解的有界性问题。3-45 计算流体力学计算流体力学 前面我们讲述过的对流和前面我们讲述

29、过的对流和扩散项的差分离散格式的阶可扩散项的差分离散格式的阶可以通过对单元面两边的节点值以通过对单元面两边的节点值在单元面上按泰勒级数展开来在单元面上按泰勒级数展开来获得：获得： 3-46 计算流体力学计算流体力学 将上述两式相减即得：将上述两式相减即得： 所以，当取所以，当取 为误差项时，用为误差项时，用 近似近似 的精度是二阶的精度。的精度是二阶的精度。3-47 计算流体力学计算流体力学 另外，将前面两式相加即得：另外，将前面两式相加即得： 所以，当取所以，当取 为误差项时，用中心差分格式为误差项时，用中心差分格式 近似近似 的精度是二阶的精度。而如果的精度是二阶的精度。而如果采用上风差分

30、近似采用上风差分近似 或或 (根据流动方向来根据流动方向来定定) 则是一阶精度。则是一阶精度。3-48 计算流体力学计算流体力学u 对流项高级离散方法对流项高级离散方法 在了解了前面对流项离散的基础方法的前在了解了前面对流项离散的基础方法的前提下，下面将进一步介绍一些常用的高级离散提下，下面将进一步介绍一些常用的高级离散方法发。方法发。 1 1，混合法，混合法(Splading,1972)(Splading,1972) 这种方法的思路是：这种方法的思路是：当当|Pe|Pe|小于或等于小于或等于2 2时，采用中心差分格式离时，采用中心差分格式离散对流项，当散对流项，当|Pe|Pe|大于大于2 2

31、时，采用上风差分格时，采用上风差分格式离散对流项式离散对流项 譬如，对于譬如，对于u0,0,3-49 计算流体力学计算流体力学其中，其中，对每个面对每个面(e(e或或w)w)：3-50 计算流体力学计算流体力学 对混合法的评价：对混合法的评价： a) a) 这种方法是满足对流离散特性中的守恒性、这种方法是满足对流离散特性中的守恒性、输运性和有界性的。输运性和有界性的。 b) b) 由于这种方法的非常稳定和健壮，使其在过去由于这种方法的非常稳定和健壮，使其在过去很长一段时间内成为商用很长一段时间内成为商用CFDCFD程序中最流行的程序中最流行的方法。方法。c) c) 但由于实际工程中的流动问题常

32、常是高对流低但由于实际工程中的流动问题常常是高对流低扩散的，因而，混合法实际上主要还是一阶的扩散的，因而，混合法实际上主要还是一阶的上风差分格式。所以，现代上风差分格式。所以，现代CFDCFD试图寻求更高试图寻求更高精度的对流离散方法。精度的对流离散方法。3-51 计算流体力学计算流体力学 2 2，QUICKQUICK法法(Leonard,1979)(Leonard,1979) QUICKQUICK是是QUadratic Interpolation for QUadratic Interpolation for Convective KinematicsConvective Kinematic

33、s的缩写。它通过三个节的缩写。它通过三个节点拟合二次多项式而得到的三阶对流离散方法点拟合二次多项式而得到的三阶对流离散方法。对于每一个单元面，。对于每一个单元面， QUICKQUICK法不仅采用了法不仅采用了单元两边的节点，并单元两边的节点，并增加了上风节点的上增加了上风节点的上风节点风节点, ,右面给出了东右面给出了东(e)面的情况：面的情况：3-52 计算流体力学计算流体力学 为了强调与单元面相关的通量的守恒特性，在为了强调与单元面相关的通量的守恒特性，在后面三点插值公式推导中将采用如下的记号：后面三点插值公式推导中将采用如下的记号： 和和 分别表示下风分别表示下风(Downwind)(D

34、ownwind)节点节点, ,上风节上风节点点(Upwind)(Upwind)和上上风节点和上上风节点(Upwind-Upwind)(Upwind-Upwind) ：3-53 计算流体力学计算流体力学 通过这些节点拟合二次通过这些节点拟合二次多项式，可得多项式，可得QUICKQUICK方法方法的公式：的公式： 譬如，当譬如，当u00时，时， 3-54 计算流体力学计算流体力学 于是有：于是有： 其中，其中， 3-55 计算流体力学计算流体力学 对对QUICKQUICK法的评价：法的评价： a) a) 这种方法具有这种方法具有3 3阶精度。阶精度。 b) b) 满足守恒性满足守恒性 c) c)满

35、足输运性，因为第三个节点是根据上风节点满足输运性，因为第三个节点是根据上风节点来选取的来选取的 d) d) 不满足有界性不满足有界性( (譬如，譬如，u0u0时，系数时，系数aE为负值为负值) ) 总之，除了不满足有界性外，总之，除了不满足有界性外， QUICKQUICK法是使法是使用最广泛的高阶离散方法。但需要注意的是，用最广泛的高阶离散方法。但需要注意的是，对于湍流问题，有界性是很重要的，譬如，对于湍流问题，有界性是很重要的，譬如， 和和 要求是正定的。要求是正定的。3-56 计算流体力学计算流体力学 3 3，通量限制法，通量限制法(Flux Limited)(Flux Limited)

36、到目前为止，我们所看到的方法都是针对到目前为止，我们所看到的方法都是针对矩阵系数矩阵系数aF为常数的情况为常数的情况( (即与解变量无关即与解变量无关) )，可以证明这些方法中，逆风差分格式是唯一一可以证明这些方法中，逆风差分格式是唯一一个无条件满足有界性的方法，但该方法只有个无条件满足有界性的方法，但该方法只有1 1阶精度。阶精度。QUICKQUICK法通过三点来拟合三阶多项式法通过三点来拟合三阶多项式来计算单元面上的通量，这样计算出来的单元来计算单元面上的通量，这样计算出来的单元面上的值会超出插值点值的范围，即不在面上的值会超出插值点值的范围，即不在 的最的最 大值和最小值之间。为了避免上

37、述情况，大值和最小值之间。为了避免上述情况，现代现代CFDCFD方法采用方法采用解相关的限制器解相关的限制器( (solution-dependent limiters) )来强制满足有界性，同时保来强制满足有界性，同时保持较高精度。持较高精度。3-57 计算流体力学计算流体力学 对于三点插值方法，如果对于三点插值方法，如果 或或 则称变量则称变量 是单调增加或单调减小。如果局部不是单调增加或单调减小。如果局部不是单调是单调( (增加或减少增加或减少) )的，则满足有界性的必要的，则满足有界性的必要条件是必须默认逆风方法条件是必须默认逆风方法( (即即 ) ) 。下面给出。下面给出两个自编程序

38、常用的两种通量限制器法。两个自编程序常用的两种通量限制器法。3-58 计算流体力学计算流体力学 在两种方法中，单调变化的校验是通过检查相在两种方法中，单调变化的校验是通过检查相邻两点的变化是否具有相同的符号，即邻两点的变化是否具有相同的符号，即 单调单调 对于单调变化，总变化的分数表示为：对于单调变化，总变化的分数表示为：3-59 计算流体力学计算流体力学1 1，UMIST UMIST 法法(Lien and Leschziner, 1993)Leschziner, 1993) UMIST UMIST是是Upstream Monotonic Interpolation Upstream Mon

39、otonic Interpolation for Scalar Transport for Scalar Transport 的缩写。该方法在单调情况的缩写。该方法在单调情况下就是二阶精度下就是二阶精度QUICKQUICK的限制变化：的限制变化： 3-60 计算流体力学计算流体力学2 2，Harmonic Harmonic 法法 (Van Leer, 1974) (Van Leer, 1974)该方法在单调情况下是二阶精度：该方法在单调情况下是二阶精度：说明：说明：a)a)通量限制器法有很多种，这里只给出两种。通量限制器法有很多种，这里只给出两种。b)b)通量限制器法满足有界性，但是非线性的，

40、即通量限制器法满足有界性，但是非线性的，即矩阵是随解变量变化，因而，数值求解必须采用矩阵是随解变量变化，因而，数值求解必须采用迭代方法。迭代方法。 3-61 计算流体力学计算流体力学u高阶对流项离散方法的实现高阶对流项离散方法的实现 一般的稳态对流扩散输运方程：一般的稳态对流扩散输运方程： 根据已经学习过的扩散项和源项的标准离散方根据已经学习过的扩散项和源项的标准离散方法可得：法可得： 其中其中F F表示相邻节点，表示相邻节点，P P表示单元中心节点。表示单元中心节点。3-62 计算流体力学计算流体力学由于质量守恒由于质量守恒( ),( ),为了方便处理，从上述方为了方便处理，从上述方程等式两

41、边减去程等式两边减去 ： 于是，需要选定一种对流项离散方法于是，需要选定一种对流项离散方法( (中心差中心差分、逆风差分，分、逆风差分，QUICKQUICK，通量限制法等，通量限制法等) )来计来计算单元面上的值算单元面上的值 。由于多数矩阵求解方法。由于多数矩阵求解方法均要求系数是正定的，一种比较通用的处理方均要求系数是正定的，一种比较通用的处理方法就是将法就是将 分解成逆风值分解成逆风值 ( (因为它能保证因为它能保证系数正定系数正定) )以及二者的差以及二者的差 ，即：，即：3-63 计算流体力学计算流体力学 方程右端第一项保证系数方程右端第一项保证系数aF正定，第二项正定，第二项移到源

42、项，称作延迟校正。移到源项，称作延迟校正。 3-64 计算流体力学计算流体力学于是于是 其中其中 延迟校正项延迟校正项(deferred correction)(deferred correction)必须根据当必须根据当前值计算前值计算( (即，在一个迭代中保持为常量即，在一个迭代中保持为常量) )3-65 计算流体力学计算流体力学u曲线网格曲线网格 非笛卡儿坐标系统称作非笛卡儿坐标系统称作曲线坐标，其坐标表示为曲线坐标，其坐标表示为 或或 ，坐标方向随，坐标方向随空间位置而变。空间位置而变。 坐标线相互垂直的坐标系统称作正交坐标系坐标线相互垂直的坐标系统称作正交坐标系，包括笛卡儿坐标，柱坐

43、标系统以及球坐标系，包括笛卡儿坐标，柱坐标系统以及球坐标系统。通过笛卡儿网格单元东面的通量可表示为统。通过笛卡儿网格单元东面的通量可表示为： 3-66 计算流体力学计算流体力学 正交网格可采用相同的方法，即单元面法线方正交网格可采用相同的方法，即单元面法线方向与坐标线一致，所以，单元面上通量的离散可向与坐标线一致，所以，单元面上通量的离散可按坐标方向进行。按坐标方向进行。 而对于大多数曲线网格，坐标线是不正交的而对于大多数曲线网格，坐标线是不正交的。单元面的法线不必与坐标线一致，于是，与。单元面的法线不必与坐标线一致，于是，与 相关的扩散通量不能仅仅由单元面两边的节点相关的扩散通量不能仅仅由单

44、元面两边的节点来近似。来近似。 譬如，二维问题中，如果控制体东面譬如，二维问题中，如果控制体东面( ( ) )与与 一致，则法向导数为：一致，则法向导数为：3-67 计算流体力学计算流体力学 于是，斜导数于是，斜导数 可离散为可离散为 ，但与单元面平行的导数，但与单元面平行的导数 还与其他节点有关。还与其他节点有关。 一般来说：一般来说：1 1，扩散通量的非对角，扩散通量的非对角分量分量( (譬如上例中包含譬如上例中包含 的项的项) ) 需要需要移到方程右端源项作显式处理；移到方程右端源项作显式处理； 2 2，需要存储全部度量分量需要存储全部度量分量 ( (三三维中有维中有9 9个个) ) 。

45、3 3，显然非正交网格需要更多的计，显然非正交网格需要更多的计算时间，但由于能适应比较复杂的几何外形而广算时间，但由于能适应比较复杂的几何外形而广泛用于一般的求解方法中。泛用于一般的求解方法中。 3-68 计算流体力学计算流体力学u边界条件边界条件 最常用的边界最常用的边界条件有两类：条件有两类：1 1，Dirichlet边界条件，即边界上变量值是已知边界条件，即边界上变量值是已知的。譬如，固壁面上的。譬如，固壁面上u=0，以及某些面上温度，以及某些面上温度为固定值。为固定值。2, 2, Neumann边界条件，即沿边界法线方向变量边界条件，即沿边界法线方向变量的导数值是已知的。譬如，对称平面

46、以及出流的导数值是已知的。譬如，对称平面以及出流边界等。边界等。 3-69 计算流体力学计算流体力学 对于边界上的单元，对于边界上的单元， 于是，对于边界单元离散方程需作两方面的修于是，对于边界单元离散方程需作两方面的修正：正：1 1，与边界相关节点的系数置为，与边界相关节点的系数置为0 02 2，将边界通量移到源项，将边界通量移到源项3-70 计算流体力学计算流体力学 如果边界上通量是给定的，如果边界上通量是给定的，则处理非常简单。而如果是给则处理非常简单。而如果是给选定变量值，则需作一些处理。选定变量值，则需作一些处理。譬如，边界上扩散通量可离散为：譬如，边界上扩散通量可离散为： 在将边界

47、通量移到源项之前，需要对系数作一在将边界通量移到源项之前，需要对系数作一些简单的处理：些简单的处理：3-71 计算流体力学计算流体力学u代数方程的求解方法代数方程的求解方法 采用有限体积法离散后，对于每个控制体采用有限体积法离散后，对于每个控制体单元有如下的离散方程：单元有如下的离散方程： 而对整个求解域的全部单元离散方程进行而对整个求解域的全部单元离散方程进行组装后得到是如下形式的矩阵方程：组装后得到是如下形式的矩阵方程： 其中，其中， 节点值向量节点值向量 ，而系数矩阵，而系数矩阵 A A 一般为稀疏矩阵。下面介绍一些常用的算法。一般为稀疏矩阵。下面介绍一些常用的算法。3-72 计算流体力

48、学计算流体力学1 1，GaussianGaussian消去法消去法a)a)是一种直接求解方法是一种直接求解方法b)b)只适合于小问题只适合于小问题 2 2，Gauss-SeidelGauss-Seidel迭代迭代算法算法 通过连续控制体积的迭代通过连续控制体积的迭代更新公式：更新公式：这种方法易于编程而且常用于这种方法易于编程而且常用于非结构网格。但收敛较慢，常非结构网格。但收敛较慢，常常需要引入亚松弛技术。常需要引入亚松弛技术。3-73 计算流体力学计算流体力学3 3，线性，线性Gauss-SeidelGauss-Seidel迭代迭代算法算法 沿任何一坐标方向，系统沿任何一坐标方向，系统是三

49、对角的，即沿是三对角的，即沿i方向有：方向有：整个坐标线可通过三对角算法一次迭代来更新，整个坐标线可通过三对角算法一次迭代来更新，即一个迭代中的信息可在域中传播，而不象即一个迭代中的信息可在域中传播，而不象Gauss-SeidelGauss-Seidel迭代迭代算法是对节点。一个典型的迭代算法是对节点。一个典型的迭代过程包括：首先对连续的过程包括：首先对连续的i i线更新，然后线更新，然后j j线，然后线，然后k k线。这种方法在块结构化网格中使用最广。也是线。这种方法在块结构化网格中使用最广。也是许多自编程序使用的算法基础。许多自编程序使用的算法基础。3-74 计算流体力学计算流体力学u收敛

50、准则收敛准则 在所有的迭代方法中，当迭代误差达到预在所有的迭代方法中，当迭代误差达到预先设定的精度时即停止迭代。残值误差是对所先设定的精度时即停止迭代。残值误差是对所有控制体积的残差求和，譬如：有控制体积的残差求和，譬如： 绝对误差：绝对误差： 均方根误差均方根误差(rms)(rms)：每个控制体积的残差为：每个控制体积的残差为：3-75 计算流体力学计算流体力学u亚松弛迭代亚松弛迭代 当前面介绍的迭代方法中应用于非线性、当前面介绍的迭代方法中应用于非线性、耦合的流体动力学变量耦合的流体动力学变量 时，每个时，每个迭代中变量的变化可能会很大导致方法失稳。迭代中变量的变化可能会很大导致方法失稳。

51、为了避免这种情况，采用亚松弛法减小变量的为了避免这种情况，采用亚松弛法减小变量的变化。变化。 对代数方程对代数方程迭代更新可写为：迭代更新可写为：3-76 计算流体力学计算流体力学当只采用更新值的分数时当只采用更新值的分数时重新整理后：重新整理后：所以，亚松弛算法，只需简单的改变系数即可所以，亚松弛算法，只需简单的改变系数即可实现：实现： 可以证明，亚松弛算法可使系数矩阵对角占可以证明，亚松弛算法可使系数矩阵对角占优，因为优，因为 更小。更小。3-77 计算流体力学计算流体力学u小结小结 1 1，标量，标量 的积分守恒形式为：的积分守恒形式为： 变化率净通量源变化率净通量源2 2，通量是指通过

52、表面的输运率，它由两部分，通量是指通过表面的输运率，它由两部分组成：组成：a) a) 对流随流动的迁移对流随流动的迁移b) b) 扩撒分子随机运动或湍流波动引起的输运扩撒分子随机运动或湍流波动引起的输运3 3， 对于每个控制体积，标量输运方程离散后对于每个控制体积，标量输运方程离散后的代数方程为：的代数方程为：其中，求和是对邻近节点进行的。其中，求和是对邻近节点进行的。3-78 计算流体力学计算流体力学4 4，将每个控制体积的离散方程进行组装后得，将每个控制体积的离散方程进行组装后得到的矩阵方程是有限带宽的，典型的求解这类到的矩阵方程是有限带宽的，典型的求解这类矩阵方程的方法包括矩阵方程的方法

53、包括Gauss-SeidelGauss-Seidel迭代迭代算法以算法以及线性及线性Gauss-SeidelGauss-Seidel迭代迭代算法。算法。5 5，源项一般线性化为如下形式：，源项一般线性化为如下形式：6 6， 扩撒项一般采用中心差分近似：扩撒项一般采用中心差分近似：3-79 计算流体力学计算流体力学7 7，对流方法是指为了计算对流通量而采用的，对流方法是指为了计算对流通量而采用的近似单元面上近似单元面上 值的方法，包括中心差分、逆值的方法，包括中心差分、逆风差分，混合法，风差分，混合法，QUICKQUICK，通量限制法等。，通量限制法等。8 8，离散方法需要满足的一般要求：，离散

54、方法需要满足的一般要求： 一致性一致性 守恒性守恒性 有界性有界性 稳定性稳定性 输运性输运性 精确性精确性3-80 计算流体力学计算流体力学9 9，由于上述离散特性要求，矩阵系数应满足，由于上述离散特性要求，矩阵系数应满足如下条件：如下条件： 10 10，边界条件的实现通过将边界通量移到源，边界条件的实现通过将边界通量移到源项。项。 11 11，在耦合非线性问题的求解一般需要采用，在耦合非线性问题的求解一般需要采用亚松弛技术。亚松弛技术。3-81 计算流体力学计算流体力学2 动量方程动量方程u动量的特征动量的特征u压力速度耦合压力速度耦合u压力校正法压力校正法u小结小结3-82 计算流体力学

55、计算流体力学 每一个流体流动问题都满足质量和动量守恒每一个流体流动问题都满足质量和动量守恒方程。方程。 而每个动量分量满足标量输运方程，对于单而每个动量分量满足标量输运方程，对于单个控制体积，标量输运方程为：个控制体积，标量输运方程为： 当输运量为动量时：当输运量为动量时： 3-83 计算流体力学计算流体力学u 动量方程的特殊属性动量方程的特殊属性 标量输运方程系数的最重要的特征就是其标量输运方程系数的最重要的特征就是其不是固定值，而是与解相关的，因为质量通量不是固定值，而是与解相关的，因为质量通量C C与速度相关，因此，运动方程必须采用：与速度相关，因此，运动方程必须采用： (1) 耦合求解

56、耦合求解 (2) 迭代方法迭代方法 特别是当特别是当 u,v,w时：时： 1，对流通量，对流通量(质量通量质量通量x速度，譬如，速度，譬如， )包含非线性项包含非线性项 ， 以及矩阵系数与以及矩阵系数与解相关。解相关。3-84 计算流体力学计算流体力学2，质量方程和动量方程是耦合的，因为质量，质量方程和动量方程是耦合的，因为质量守恒和动量守恒方程中都涉及速度变量。守恒和动量守恒方程中都涉及速度变量。3，压力在每个动量方程中都出现，于是，压力在每个动量方程中都出现，于是a)进进一步使方程耦合，一步使方程耦合，b)需要确定压力的方法。需要确定压力的方法。4，总的有，总的有4个方程个方程(一个质量守

57、恒方程和一个质量守恒方程和3个动个动量守恒方程量守恒方程)，要使方程有解，只能有四个未，要使方程有解，只能有四个未知量知量:u,v,w,+? 在可压缩流动中，质量守恒提供密度在可压缩流动中，质量守恒提供密度 的的输运方程，能量方程或熵方程提供温度输运方程，能量方程或熵方程提供温度T 的输的输运方程，而压力运方程，而压力 p 由热力学关系，即状态方程由热力学关系，即状态方程确定，譬如，理想气体的状态方程为确定，譬如，理想气体的状态方程为3-85 计算流体力学计算流体力学 在不可压缩流动中，密度的变化与压力没有在不可压缩流动中，密度的变化与压力没有联系。而对动量方程的解必须满足质量方程要求联系。而

58、对动量方程的解必须满足质量方程要求最终产生了压力方程。最终产生了压力方程。 在许多程序中，质量方程和动量方程采用迭在许多程序中，质量方程和动量方程采用迭代方法依次求解，这种求解方法称作分离解法，代方法依次求解，这种求解方法称作分离解法，下面是分离解法的伪代码：下面是分离解法的伪代码：3-86 计算流体力学计算流体力学u 压力和速度的耦合压力和速度的耦合 需要回答三个问题：需要回答三个问题： (1) 压力和速度是如何联系的？压力和速度是如何联系的？ (2) 压力方程是如何产生的？压力方程是如何产生的？ (3) (3) 压力和速度应该存储在同一位置吗？压力和速度应该存储在同一位置吗？1，压力与速度

59、的联系，压力与速度的联系 在动量方程中，压力力在动量方程中，压力力作为源项出现在动量方程作为源项出现在动量方程中，譬如，中，譬如，x动量方程，动量方程， 压力力压力力3-87 计算流体力学计算流体力学于是，离散后的动量方程具有如下形式：于是，离散后的动量方程具有如下形式： 其他力其他力 压力力压力力 净通量净通量因此，因此， 其中，其中，可以看出：可以看出：a) 动量方程中的力项提供了速度与压力的联动量方程中的力项提供了速度与压力的联系系3-88 计算流体力学计算流体力学b) 速度依赖于压力梯度，或者说，离散后速速度依赖于压力梯度，或者说，离散后速度依赖于离单元中心点距离度依赖于离单元中心点距

60、离 的两边的压力的两边的压力差。即差。即 动量方程动量方程 用符号表示就是：用符号表示就是： 其中，其中， 表示中心差分。表示中心差分。 2/x3-89 计算流体力学计算流体力学2，压力方程的产生，压力方程的产生 将上述速度与压力差的关系将上述速度与压力差的关系代入连续方程后有：代入连续方程后有：这和离散的标量输运方程具有相同的形式。这和离散的标量输运方程具有相同的形式。 所以，动量方程提供了速度和压力的关系所以，动量方程提供了速度和压力的关系，进而将这个关系代入质量方程就得到压力方，进而将这个关系代入质量方程就得到压力方程。或者说，质量守恒对动量方程的约束而产程。或者说，质量守恒对动量方程的