计算流体力学数值方法

3-31
计算流体力学
为了简化书写，将质量通量用C表示，譬如 ，单元e面的质量通量表示为：
扩散项和源项的离散采用与前面相同的方 法，于是有：
现在的问题是如何离散对流项？也就是怎样利 用相邻节点值来确定单元面上的变量值。这类 方法就称作对流方案。
3-32
计算流体力学
中心差分法 离散对流项的最明然 的方法就是线性插值，即 中心差分法：
计算流体力学
第三章
数值方法
3-1
计算流体力学
主要内容
空间离散技术 标量输运方程 动量方程 时间离散技术

3-2
计算流体力学
1 标量输运方程

有限体积法 一维对流扩散方程 扩散项的离散 源项的离散 代数方程的组装 二维和三维问题 对流项离散基础 离散特性
3-3
计算流体力学
3-25
计算流体力学
代数方程的组装 当速度u=0时，定常扩散 问题在每个单元上的离散为：
经过整理后：
或表示成如下标量输运方程规范形式：
3-26
计算流体力学
上述方程为离散后的标量输运方程规范形式 ，每个求解变量的输运方程都有相同的形式，只 是各自的矩阵系数不同。对于纯扩散问题，矩阵 系数为：
3-27
3-37
计算流体力学
上风差分法 在上风差分格式种， 单元面上变量的值用其 上风节点的值近似。譬 如，对于单元东(e)面:
于是将上述公式可归结为：
3-38
计算流体力学
采用上风差分格式后，定常对流－扩散问题 的规范离散格式为：
其中：
3-39
计算流体力学
当将上风差分法应用于前述与中心差分同样的 问题后可观察到： (1) 当
计算流体力学
解析解：
d dT (kA ) 0 dx dx d 2T 0 2 dx T c1 x c0 T ( x 0) 100 T ( x 1) 500 T 400x 100
3-21
计算流体力学
控制方程扩散项的离散 梯度扩散项的离散几乎 总是采用中心差分格式：
3-44
计算流体力学
提示： (1) 有限体积法中，通量是在面上计算的，而 不是在节点处。 (2) 上述对梯度扩散项的近似在空间上具有二 阶精度，后面将会给出证明。
3-22
计算流体力学
(3) 如果扩散系数也是 变量的话，它在单元面 上的值必须通过插值得到。
(4) 采用中心差分格式意 味着在单元面两边的节点上权重相同。这 和扩撒的物理意义是一致的，因为扩散在 所有方向作用相同。这后面将要讨论的对流 是不一样的，对流是具有明确的方向性的。
所得结果精度不如中心差分法，这是意料之 中的。 (2) 当Pe=4时，得到的解不是很精确，但不再 出现明显的摆动！ 于是，便出现了精度和有界性(无摆动)之间 如何取舍的问题？下面将介绍一些离散方程的属 性。
3-40
计算流体力学
离散属性 1，一致性 一致性是指当网格间距趋向于0时，离散 方程等价于连续方程。譬如， 就是 的一致近似。 2，守恒性 即从一个单元流出的 通量一定等于流入相邻单元 的通量。这是有限体积法优 越有限差分和有限元法的地方。
e
3-16
计算流体力学
非边界面上通量的计算：
3-17
计算流体力学
边界面上通量的计算：
3-18
计算流体力学
T1 TB1 T2 T1 0 x / 2 x T2 T1 T3 T2 0 x x T3 T2 T4 T3 0 x x T4 T3 T5 T4 0 x x T5 T4 TB 2 T5 0 x x / 2
3-9
计算流体力学
于是单元各面的面积表示为Aw, Ae, As, An Ab, At 。体积为V。对于二维问题，可以视为单 位厚度为1的一层单元( ) 对于结构网格，可以交替使用ijk下标表示单 元节点，譬如，
3-10
计算流体力学
一维对流－扩散方程 我们首先讨论一维稳态对流扩散方程主要基 于下面的考虑： (1) 它使问题分析大大简化 (2) 离散方程可以进行手算。 (3) 尽管只是一维的，但要扩展到二维或三维 是非常直接的 (4) 实际上，通量(对流和扩散)的离散一般是 沿坐标方向进行的，即分别沿i,j,k线进行。 (5) 有许多重要的理论问题是一维的。
3-12
计算流体力学
于是， 的对流扩散方程为：
其中，S表示单位长度的源。将等式两边除以 ，并取极限 ，可得到相应的微分方程：
注意：1，这里面积A是表示可变截面的准一维问题， 对于真正的一维问题，只需令A=1,即 2，这里假设 和S均为常数，但在一般的 CFD问题中，u本身也是问题的解变量。
3-13
计算流体力学
例子
1，纯扩散问题 如图所示的隔热棒，长度 1m，截面为1cmx1cm的方 形截面，棒的两端为固定 温度，分别为100度和500 度。穿过任意截面A的热通 量由下式给定： 其中，热传导系数 (a) 将棒划分为5个控制体，并用有限体积分析沿棒 的温度分布 (b) 写出沿棒温度分布的微分方程 (c) 求出微分方程的解析解，并与(a)的解进行比较
3-7
计算流体力学
计算网格可以是结构网格或非结构网格，笛 卡儿网格或非笛卡儿网格。 最常用的网格形式包括基于单元中心的存储 方式和基于单元顶点的存储方式两种。在后面的 讨论中将会看到，不是所有的变量都必须存储在 同样的位置。
3-8
计算流体力学
为了简单起见，本课程将只考虑结构化笛卡 儿网格，和采用基于单元中心的存储方式，即单 元变量 存储在单元中心节点位置上。 右图为一典型三维控制 体积，其中P点为单元中心 存储节点。相对于中心节点， 沿坐标方向通常表示为west, east, south, north, bottom, top. 其中，小写字母w, e, s, n, b, t 表示单元面，大写字母W, E, S, N, B, T表示中心节点的 相邻节点。
3-29
计算流体力学
对于多维问题，离散后的代数方程仍然具有 如下形式：
上述方程只是对于单个 控制体单元的，而对所 有单元进行组装便得到 一组节点值向量 的 矩阵方程。对于二维问 题，最终得到的矩阵具有如图所示的带状形式 。三维问题将增加非对角元素！
3-30
计算流体力学
对流项的离散(I) 纯扩散问题(u=0)只是一 种理想情况，对于典型的工 程问题，对流项远远超过扩 散项，因为雷诺数非常大，即惯性力和粘性力 的比值很大。 定常的对流－扩散方程方程：
3-23
计达式，这是守恒定律所要 求的。即从一个单元流出的 通量一定等于流入相邻单元 的通量。这是有限体积法优越有限差分和有限 元法的地方。 在有限体积法中，常常定义一个扩散输运 系数D，譬如，对于单元东面有：
3-24
计算流体力学
控制方程源项的离散 当源与流体的量成比例的情况下，每个单 元的总的源强度可简单的表示为： 单位体积的源x单元体积＝SxV 其中，S表示平均源密度，一般取单元中心 节点处的值。 一般来说，S常常是与变量 相关的，譬 如，如果变量为温度的话，牛顿冷却定律认为 ：单位时间热量的损失与相对与环境的温差成 比例： 一般将源项分解为与解相关和无关两部分 ：
3-42
计算流体力学
6，阶(解的精度) 阶是指随着空间网格尺寸的减小，误差衰 减的速度。譬如，对于均匀空间网格 ，在某 种数值方法中，如果当空间网格尺寸趋向于0 ，而误差与空间网格尺寸的n次方成 正比 ，则称这种数值方法是n阶的。 某种数值方法的阶可以通过泰勒级数展开 式中的首阶误差项确定(后面将讨论) 。 注意：此处的误差是指展开式中的理论截 断误差，而不考虑计算机的舍入误差(即计算 机只能存储有限位数的数)
3-14
计算流体力学
这个问题是求解一维对流－扩散方程：
上面的方程也可写成如下的积分方程：
由于本例只考虑扩散，即没有对流和源项：
3-15
计算流体力学
对于本例的温度T的纯扩撒问题，最终有如下 的微分方程和积分方程：
d dT (kA ) 0 dx dx dT kA 0 dx w
计算流体力学
上述方程为只是一个控制体的离散方程，如 果将所有控制体单元中心节点处的变量值 组装 到一个向量中，就会得到如下形式的代数方程：
求解这种代数方程可采用三对角矩阵算法
3-28
计算流体力学
二维和三维的扩展 前面讨论的虽然只是一维 问题，但要推广到二维和三维 是非常直接的。 对于多维问题，经过单元 面的净通量只需将所有面上的 通量求代数和即可！譬如，对于三维控制体， 其守恒方程可表示成：
3-19
计算流体力学
3T1 T2 200 T1 2T2 T3 0 T2 2T3 T4 0 T3 2T4 T5 0 T4 3T5 1000
T1 140 T1 220 T3 300 T4 380 T1 460
3-20
3-41
计算流体力学
3，输运性 方向性的影响反应在对流项的离散方法中 。实际上，就是在上风节点赋予较高的权重。
4，有界性 在无源的对流－扩散问题中的解一定介 于周围节点流动变量最大值和最小值之间。 5, 稳定性 稳定性决定了求解过程能否得到最终的 解。它不涉及解的精度。稳定即意味着误差不 随求解过程增加。
3-43
计算流体力学
高阶精度可通过采用更多的节点值近似来获 得。一个节点允许的最高精度为1阶，两个节点允 许的最高精度为2阶，依此类推。 理论上讲，某种数值方法的精度越高，随着 网格的加密，误差减小的就越大 。也就是说，采 用高精度的数值方法，只需较少的网格数即可获 得要求的精度。 但是，高阶精度的方法常常需要更多的计算 时间，而且常常会导致解的有界性问题。