计算流体力学数值方法讲解ppt课件.ppt

第三章 数值方法,主要内容,空间离散技术 标量输运方程 动量方程时间离散技术,1 标量输运方程,有限体积法一维对流扩散方程扩散项的离散源项的离散代数方程的组装二维和三维问题对流项离散基础离散特性,对流项高级离散方法 高阶对流方法的实现曲线网格边界条件代数方程的求解小结,对于任意控制体内一个守恒物理量: 控制体内物理量的变化率+经过控制体边界面的净通量控制体积内的源(汇) 经过控制体边界面的总的通量由对流(随流动的迁移)和扩散(由随机分子运动或湍流运动导致的净输运)两部分组成。如果用 表示单位质量流体的守恒量，则通用标量输运方程(或称对流扩散方程)可表示为：,1，有限体积法直接对上式进行离散2，本章只考虑稳态问题，即上式左边第一项为零,有限体积法(FVM) (1) 定义流场求解域几何形状 (2) 将求解域划分为计算网格，即一组互不重叠的有限体或单元。 (3) 基于上述划分的单元对积分方程进行离散，即用节点值来近似。 (4) 对得到的离散方程进行数值求解。,计算网格可以是结构网格或非结构网格，笛卡儿网格或非笛卡儿网格。 最常用的网格形式包括基于单元中心的存储方式和基于单元顶点的存储方式两种。在后面的讨论中将会看到，不是所有的变量都必须存储在同样的位置。,为了简单起见，本课程将只考虑结构化笛卡儿网格，和采用基于单元中心的存储方式，即单元变量 存储在单元中心节点位置上。 右图为一典型三维控制体积，其中P点为单元中心存储节点。相对于中心节点，沿坐标方向通常表示为west, east, south, north, bottom, top. 其中，小写字母w, e, s, n, b, t 表示单元面，大写字母W,E, S, N, B, T表示中心节点的相邻节点。,于是单元各面的面积表示为Aw, Ae, As, An Ab, At 。体积为V。对于二维问题，可以视为单位厚度为1的一层单元( ) 对于结构网格，可以交替使用ijk下标表示单元节点，譬如，,一维对流扩散方程 我们首先讨论一维稳态对流扩散方程主要基于下面的考虑： (1) 它使问题分析大大简化 (2) 离散方程可以进行手算。 (3) 尽管只是一维的，但要扩展到二维或三维是非常直接的 (4) 实际上，通量(对流和扩散)的离散一般是沿坐标方向进行的，即分别沿i,j,k线进行。 (5) 有许多重要的理论问题是一维的。,对于如图所示的一维控制体，物理量的守恒可表述为如下关系式： 通量e (fluxe)- 通量w (fluxw) 源(source) 这里的通量是指穿过单元表面的输运率。 如果 表示单位质量的输运量，则总的通量为对流通量和扩散通量之和，其中： 对流通量 扩散通量,于是， 的对流扩散方程为： 其中，S表示单位长度的源。将等式两边除以 ，并取极限 ，可得到相应的微分方程：注意：1，这里面积A是表示可变截面的准一维问题，对于真正的一维问题，只需令A=1,即 2，这里假设 和S均为常数，但在一般的CFD问题中，u本身也是问题的解变量。,例子 1，纯扩散问题 如图所示的隔热棒，长度1m，截面为1cmx1cm的方形截面，棒的两端为固定温度，分别为100度和500度。穿过任意截面A的热通量由下式给定： 其中，热传导系数将棒划分为5个控制体，并用有限体积分析沿棒的温度分布写出沿棒温度分布的微分方程 求出微分方程的解析解，并与(a)的解进行比较,这个问题是求解一维对流扩散方程： 上面的方程也可写成如下的积分方程： 由于本例只考虑扩散，即没有对流和源项：,对于本例的温度T的纯扩撒问题，最终有如下的微分方程和积分方程：,非边界面上通量的计算：,边界面上通量的计算：,解析解：,控制方程扩散项的离散 梯度扩散项的离散几乎总是采用中心差分格式： 提示： (1) 有限体积法中，通量是在面上计算的，而不是在节点处。 (2) 上述对梯度扩散项的近似在空间上具有二阶精度，后面将会给出证明。,(3) 如果扩散系数也是 变量的话，它在单元面 上的值必须通过插值得到。 (4) 采用中心差分格式意 味着在单元面两边的节点上权重相同。这 和扩撒的物理意义是一致的，因为扩散在 所有方向作用相同。这后面将要讨论的对流 是不一样的，对流是具有明确的方向性的。,(5) 单元西面也应有相似的表达式，这是守恒定律所要求的。即从一个单元流出的通量一定等于流入相邻单元的通量。这是有限体积法优越有限差分和有限元法的地方。 在有限体积法中，常常定义一个扩散输运系数D，譬如，对于单元东面有：,控制方程源项的离散 当源与流体的量成比例的情况下，每个单元的总的源强度可简单的表示为： 单位体积的源x单元体积SxV 其中，S表示平均源密度，一般取单元中心节点处的值。 一般来说，S常常是与变量 相关的，譬如，如果变量为温度的话，牛顿冷却定律认为：单位时间热量的损失与相对与环境的温差成比例： 一般将源项分解为与解相关和无关两部分：,代数方程的组装 当速度u=0时，定常扩散问题在每个单元上的离散为： 经过整理后： 或表示成如下标量输运方程规范形式：,上述方程为离散后的标量输运方程规范形式，每个求解变量的输运方程都有相同的形式，只是各自的矩阵系数不同。对于纯扩散问题，矩阵系数为：,上述方程为只是一个控制体的离散方程，如果将所有控制体单元中心节点处的变量值 组装到一个向量中，就会得到如下形式的代数方程： 求解这种代数方程可采用三对角矩阵算法,二维和三维的扩展 前面讨论的虽然只是一维问题，但要推广到二维和三维是非常直接的。 对于多维问题，经过单元面的净通量只需将所有面上的通量求代数和即可！譬如，对于三维控制体，其守恒方程可表示成：,对于多维问题，离散后的代数方程仍然具有如下形式：上述方程只是对于单个控制体单元的，而对所有单元进行组装便得到一组节点值向量 的矩阵方程。对于二维问题，最终得到的矩阵具有如图所示的带状形式。三维问题将增加非对角元素！,对流项的离散(I) 纯扩散问题(u=0)只是一种理想情况，对于典型的工程问题，对流项远远超过扩散项，因为雷诺数非常大，即惯性力和粘性力的比值很大。 定常的对流扩散方程方程：,为了简化书写，将质量通量用C表示，譬如，单元e面的质量通量表示为： 扩散项和源项的离散采用与前面相同的方法，于是有：现在的问题是如何离散对流项？也就是怎样利用相邻节点值来确定单元面上的变量值。这类方法就称作对流方案。,中心差分法 离散对流项的最明然的方法就是线性插值，即中心差分法： 于是将上述公式代入对流扩散方程后，得到关于单个控制体单元如下形式的离散方程：,上式中的求和是对中心点的所有相邻节点进行的。各个系数表达式如下： 上式中，根据质量守恒有,为了说明中心差分格式是否适合离散对流项，下面给出了定常对流扩散问题(无源项)在下面两种情况下的解的图示： Pe=1/2 pe=4其中Pe是Peclet的缩写，称作Pe数。,上图中红圈表示计算值，实线为解析解。第一重情况吻合很好，但第二种情况却有显著的摆动，这是为什么呢？,从数学的角度来说，当Pe数大于2，系数矩阵的一个元素将是负值，即相应于该元素的节点值的增加会使中心节点值减小！ 从物理的角度来看，是因为对流过程是有方向性的。而输运属性仅仅在流动方向。但中心差分格式没有反映处方向性，而是在上风(upwind)节点和下风(downwind)节点上赋予了相同的权重。,上风差分法 在上风差分格式种，单元面上变量的值用其上风节点的值近似。譬如，对于单元东(e)面: 于是将上述公式可归结为：,采用上风差分格式后，定常对流扩散问题的规范离散格式为： 其中：,当将上风差分法应用于前述与中心差分同样的问题后可观察到： (1) 当 所得结果精度不如中心差分法，这是意料之中的。 (2) 当Pe=4时，得到的解不是很精确，但不再出现明显的摆动！ 于是，便出现了精度和有界性(无摆动)之间如何取舍的问题？下面将介绍一些离散方程的属性。,离散属性 1，一致性 一致性是指当网格间距趋向于0时，离散方程等价于连续方程。譬如， 就是 的一致近似。 2，守恒性 即从一个单元流出的通量一定等于流入相邻单元的通量。这是有限体积法优越有限差分和有限元法的地方。,3，输运性 方向性的影响反应在对流项的离散方法中。实际上，就是在上风节点赋予较高的权重。 4，有界性 在无源的对流扩散问题中的解一定介于周围节点流动变量最大值和最小值之间。 5, 稳定性 稳定性决定了求解过程能否得到最终的解。它不涉及解的精度。稳定即意味着误差不随求解过程增加。,6，阶(解的精度) 阶是指随着空间网格尺寸的减小，误差衰减的速度。譬如，对于均匀空间网格 ，在某种数值方法中，如果当空间网格尺寸趋向于0，而误差与空间网格尺寸的n次方成 正比，则称这种数值方法是n阶的。 某种数值方法的阶可以通过泰勒级数展开式中的首阶误差项确定(后面将讨论) 。 注意：此处的误差是指展开式中的理论截断误差，而不考虑计算机的舍入误差(即计算机只能存储有限位数的数),高阶精度可通过采用更多的节点值近似来获得。一个节点允许的最高精度为1阶，两个节点允许的最高精度为2阶，依此类推。 理论上讲，某种数值方法的精度越高，随着网格的加密，误差减小的就越大 。也就是说，采用高精度的数值方法，只需较少的网格数即可获得要求的精度。 但是，高阶精度的方法常常需要更多的计算时间，而且常常会导致解的有界性问题。,前面我们讲述过的对流和扩散项的差分离散格式的阶可以通过对单元面两边的节点值在单元面上按泰勒级数展开来获得：,将上述两式相减即得： 所以，当取 为误差项时，用 近似 的精度是二阶的精度。,另外，将前面两式相加即得： 所以，当取 为误差项时，用中心差分格式 近似 的精度是二阶的精度。而如果采用上风差分近似 或 (根据流动方向来定) 则是一阶精度。,对流项高级离散方法 在了解了前面对流项离散的基础方法的前提下，下面将进一步介绍一些常用的高级离散方法发。 1，混合法(Splading,1972) 这种方法的思路是：当|Pe|小于或等于2时，采用中心差分格式离散对流项，当|Pe|大于2时，采用上风差分格式离散对流项 譬如，对于u0,其中，对每个面(e或w)：,对混合法的评价： a) 这种方法是满足对流离散特性中的守恒性、输运性和有界性的。 b) 由于这种方法的非常稳定和健壮，使其在过去很长一段时间内成为商用CFD程序中最流行的方法。c) 但由于实际工程中的流动问题常常是高对流低扩散的，因而，混合法实际上主要还是一阶的上风差分格式。所以，现代CFD试图寻求更高精度的对流离散方法。,2，QUICK法(Leonard,1979) QUICK是QUadratic Interpolation for Convective Kinematics的缩写。它通过三个节点拟合二次多项式而得到的三阶对流离散方法。对于每一个单元面， QUICK法不仅采用了单元两边的节点，并增加了上风节点的上风节点,右面给出了东(e)面的情况：,为了强调与单元面相关的通量的守恒特性，在后面三点插值公式推导中将采用如下的记号： 和 分别表示下风(Downwind)节点,上风节点(Upwind)和上上风节点(Upwind-Upwind) ：,通过这些节点拟合二次多项式，可得QUICK方法的公式： 譬如，当u0时，,于是有： 其中，,对QUICK法的评价： a) 这种方法具有3阶精度。 b) 满足守恒性 c)满足输运性，因为第三个节点是根据上风节点来选取的 d) 不满足有界性(譬如，u0时，系数aE为负值) 总之，除了不满足有界性外， QUICK法是使用最广泛的高阶离散方法。但需要注意的是，对于湍流问题，有界性是很重要的，譬如， 和 要求是正定的。,3，通量限制法(Flux Limited) 到目前为止，我们所看到的方法都是针对矩阵系数aF为常数的情况(即与解变量无关)，可以证明这些方法中，逆风差分格式是唯一一个无条件满足有界性的方法，但该方法只有1阶精度。QUICK法通过三点来拟合三阶多项式来计算单元面上的通量，这样计算出来的单元面上的值会超出插值点值的范围，即不在 的最 大值和最小值之间。为了避免上述情况，现代CFD方法采用解相关的限制器(solution-dependent limiters)来强制满足有界性，同时保持较高精度。,对于三点插值方法，如果 或 则称变量 是单调增加或单调减小。如果局部不是单调(增加或减少)的，则满足有界性的必要条件是必须默认逆风方法(即 ) 。下面给出两个自编程序常用的两种通量限制器法。,在两种方法中，单调变化的校验是通过检查相邻两点的变化是否具有相同的符号，即 单调 对于单调变化，总变化的分数表示为：,1，UMIST 法(Lien and Leschziner, 1993) UMIST是Upstream Monotonic Interpolation for Scalar Transport 的缩写。该方法在单调情况下就是二阶精度QUICK的限制变化：,2，Harmonic 法 (Van Leer, 1974)该方法在单调情况下是二阶精度：说明：a)通量限制器法有很多种，这里只给出两种。b)通量限制器法满足有界性，但是非线性的，即矩阵是随解变量变化，因而，数值求解必须采用迭代方法。,高阶对流项离散方法的实现 一般的稳态对流扩散输运方程： 根据已经学习过的扩散项和源项的标准离散方法可得： 其中F表示相邻节点，P表示单元中心节点。,由于质量守恒( ),为了方便处理，从上述方程等式两边减去 ： 于是，需要选定一种对流项离散方法(中心差分、逆风差分，QUICK，通量限制法等)来计算单元面上的值 。由于多数矩阵求解方法均要求系数是正定的，一种比较通用的处理方法就是将 分解成逆风值 (因为它能保证系数正定)以及二者的差 ，即：,方程右端第一项保证系数aF正定，第二项移到源项，称作延迟校正。,于是 其中 延迟校正项(deferred correction)必须根据当前值计算(即，在一个迭代中保持为常量),曲线网格 非笛卡儿坐标系统称作曲线坐标，其坐标表示为 或 ，坐标方向随空间位置而变。 坐标线相互垂直的坐标系统称作正交坐标系，包括笛卡儿坐标，柱坐标系统以及球坐标系统。通过笛卡儿网格单元东面的通量可表示为：,正交网格可采用相同的方法，即单元面法线方向与坐标线一致，所以，单元面上通量的离散可按坐标方向进行。 而对于大多数曲线网格，坐标线是不正交的。单元面的法线不必与坐标线一致，于是，与 相关的扩散通量不能仅仅由单元面两边的节点来近似。 譬如，二维问题中，如果控制体东面( )与 一致，则法向导数为：,于是，斜导数 可离散为 ，但与单元面平行的导数 还与其他节点有关。 一般来说：1，扩散通量的非对角分量(譬如上例中包含 的项) 需要移到方程右端源项作显式处理； 2，需要存储全部度量分量 (三维中有9个) 。3，显然非正交网格需要更多的计算时间，但由于能适应比较复杂的几何外形而广泛用于一般的求解方法中。,边界条件 最常用的边界条件有两类：1，Dirichlet边界条件，即边界上变量值是已知的。譬如，固壁面上u=0，以及某些面上温度为固定值。2, Neumann边界条件，即沿边界法线方向变量的导数值是已知的。譬如，对称平面以及出流边界等。,对于边界上的单元， 于是，对于边界单元离散方程需作两方面的修正：1，与边界相关节点的系数置为02，将边界通量移到源项,如果边界上通量是给定的，则处理非常简单。而如果是给选定变量值，则需作一些处理。譬如，边界上扩散通量可离散为： 在将边界通量移到源项之前，需要对系数作一些简单的处理：,代数方程的求解方法 采用有限体积法离散后，对于每个控制体单元有如下的离散方程： 而对整个求解域的全部单元离散方程进行组装后得到是如下形式的矩阵方程： 其中， 节点值向量 ，而系数矩阵 A 一般为稀疏矩阵。下面介绍一些常用的算法。,1，Gaussian消去法a)是一种直接求解方法b)只适合于小问题 2，Gauss-Seidel迭代算法 通过连续控制体积的迭代更新公式：这种方法易于编程而且常用于非结构网格。但收敛较慢，常常需要引入亚松弛技术。,3，线性Gauss-Seidel迭代算法 沿任何一坐标方向，系统是三对角的，即沿i方向有：整个坐标线可通过三对角算法一次迭代来更新，即一个迭代中的信息可在域中传播，而不象Gauss-Seidel迭代算法是对节点。一个典型的迭代过程包括：首先对连续的i线更新，然后j线，然后k线。这种方法在块结构化网格中使用最广。也是许多自编程序使用的算法基础。,收敛准则 在所有的迭代方法中，当迭代误差达到预先设定的精度时即停止迭代。残值误差是对所有控制体积的残差求和，譬如： 绝对误差： 均方根误差(rms)：每个控制体积的残差为：,亚松弛迭代 当前面介绍的迭代方法中应用于非线性、耦合的流体动力学变量 时，每个迭代中变量的变化可能会很大导致方法失稳。为了避免这种情况，采用亚松弛法减小变量的变化。 对代数方程迭代更新可写为：,当只采用更新值的分数时重新整理后：所以，亚松弛算法，只需简单的改变系数即可实现： 可以证明，亚松弛算法可使系数矩阵对角占优，因为 更小。,小结 1，标量 的积分守恒形式为： 变化率净通量源2，通量是指通过表面的输运率，它由两部分组成：a) 对流随流动的迁移b) 扩撒分子随机运动或湍流波动引起的输运3， 对于每个控制体积，标量输运方程离散后的代数方程为：其中，求和是对邻近节点进行的。,4，将每个控制体积的离散方程进行组装后得到的矩阵方程是有限带宽的，典型的求解这类矩阵方程的方法包括Gauss-Seidel迭代算法以及线性Gauss-Seidel迭代算法。5，源项一般线性化为如下形式：6， 扩撒项一般采用中心差分近似：,7，对流方法是指为了计算对流通量而采用的近似单元面上 值的方法，包括中心差分、逆风差分，混合法，QUICK，通量限制法等。8，离散方法需要满足的一般要求： 一致性 守恒性 有界性 稳定性 输运性 精确性,9，由于上述离散特性要求，矩阵系数应满足如下条件： 10，边界条件的实现通过将边界通量移到源项。 11，在耦合非线性问题的求解一般需要采用亚松弛技术。,2 动量方程,动量的特征压力速度耦合压力校正法小结,每一个流体流动问题都满足质量和动量守恒方程。 而每个动量分量满足标量输运方程，对于单个控制体积，标量输运方程为： 当输运量为动量时：,动量方程的特殊属性 标量输运方程系数的最重要的特征就是其不是固定值，而是与解相关的，因为质量通量C与速度相关，因此，运动方程必须采用： (1) 耦合求解 (2) 迭代方法 特别是当 u,v,w时： 1，对流通量(质量通量x速度，譬如， )包含非线性项 ， 以及矩阵系数与解相关。,2，质量方程和动量方程是耦合的，因为质量守恒和动量守恒方程中都涉及速度变量。3，压力在每个动量方程中都出现，于是a)进一步使方程耦合，b)需要确定压力的方法。4，总的有4个方程(一个质量守恒方程和3个动量守恒方程)，要使方程有解，只能有四个未知量:u,v,w,+? 在可压缩流动中，质量守恒提供密度 的输运方程，能量方程或熵方程提供温度T 的输运方程，而压力 p 由热力学关系，即状态方程确定，譬如，理想气体的状态方程为,在不可压缩流动中，密度的变化与压力没有联系。而对动量方程的解必须满足质量方程要求最终产生了压力方程。 在许多程序中，质量方程和动量方程采用迭代方法依次求解，这种求解方法称作分离解法，下面是分离解法的伪代码：,压力和速度的耦合 需要回答三个问题： (1) 压力和速度是如何联系的？ (2) 压力方程是如何产生的？ (3) 压力和速度应该存储在同一位置吗？1，压力与速度的联系 在动量方程中，压力力作为源项出现在动量方程中，譬如，x动量方程， 压力力,于是，离散后的动量方程具有如下形式： 其他力 压力力 净通量因此， 其中，可以看出：a) 动量方程中的力项提供了速度与压力的联系,b) 速度依赖于压力梯度，或者说，离散后速度依赖于离单元中心点距离 的两边的压力差。即 动量方程 用符号表示就是： 其中， 表示中心差分。,2，压力方程的产生 将上述速度与压力差的关系代入连续方程后有：这和离散的标量输运方程具有相同的形式。 所以，动量方程提供了速度和压力的关系，进而将这个关系代入质量方程就得到压力方程。或者说，质量守恒对动量方程的约束而产生了压力方程。,3，同位存储与线性插值 首先假设速度变量和压力储存在同一位置，而质量通量通过单元面的速度线性插值计算。 a)动量方程 单元的净压力力依赖于单元面上的压力pw和pe，而这些面上的值都需要通过插值得到：于是，i方向的净压力力为：,从上面的式子可以看出，i点压力被抵消。所以，动量方程离散后，i点的速度只与相距 的两点的压力差有关，而与i点的压力无关，事实上，pi可以取任何值而不影响相应单元的动量平衡。 b) 质量方程 控制体积单元面上的对流速度ue和uw通过插值得到：,于是质量方程可写成： 根据，速度与压力差的关系可得到下面的压力公式： 可以看出，每隔一个点的压力值之间才有联系。,换句话说，速度u与压力p存储在同一位置以及对流速度采用线性插值后导致奇数节点压力值p1,p3,p5,与偶数节点压力值p2,p4,p6解耦。这种奇偶解耦或称作跳棋盘(checkerboard)效应会使压力场产生不确定的振荡。,为了避免上述奇偶解耦效应造成的压力振荡，一般有两种常用的措施： (1) 采用交错网格技术(即速度和压力储存在不同的位置) (2) 采用同位存储，但对流速度用不同的插值方法(质量通量的计算采用单元面上的速度值) 上述两种技术都能使相邻节点的压力值相互联系，从而避免了奇偶解耦的问题。,(c) 交错网格(Harlow and Welch, 1965) 在交错网格技术中，速度分量储存在压力节点之间的中间位置上，这导致了压力和速度分属不同的控制体积单元。 而在节点变量编号上的一般惯例是速度节点和压力节点采用相同指标p(或ijk),交错网格的优点： 动量方程:现在压力正好储存在需要计算压力的节点上： 质量方程:现在速度分量正好储存在需要计算质量通量的点上：,可以看出，在交错网格技术中，pi的贡献不再被抵消，所以，也就不会出现奇偶解耦的问题。而且，得到的压力方程与前面的标量输运方程具有相同的形式 ，所以可以使用相同的子程序。 交错网格的缺点： 由于有多个不同的节点和控制体积单元组，增加了几何复杂性； 对于非笛卡儿网格，速度节点会不在压力节点之间，如右图所示。,(d) Rhie-Chow 速度插值 (Rhie 和 Chow,1983) 另外一种方法就是仍然采用速度和压力同位储存，但对用来计算质量通量的对流速度采用不同的插值方法。 动量方程使得单元中心的速度与单元面上的压力联系起来：即，,Rhie-Chow的思路就是对压力项和非压力项分别插值。 式中的上划线表示平均节点值得到上的值(严格来说是线性插值)。等价地(或者更方便地)加减中心压力差 基本上只是加减了一个与压力梯度成比例的量：,除了方程右端项第一项采用单元中心离散值的平均外，第二项采用经过单元面离散值的中心差分。采用上述插值技术后，质量方程为： 可以看出： 中心压力值不再被抵消，因而没有奇偶节点解耦问题；,这是一个压力方程，而且与标量输运方程的离散方程形式相同： 压力速度耦合是纳维司托克斯方程的重要特征，对于笛卡儿网格，交错网格技术是最有效的办法，但是，对于非笛卡儿网格(曲线网格)同位存储网格是最常用的，同时也是几乎所有通用CFD软件采用的方法。,压力校正法 什么是压力校正法： (1) 它是一种求解不可压缩流的数值迭代方法 (2) 用于导出同时满足质量和动量守恒的速度场和压力场 (3) 包括下列速度场和压力场的更新： 首先，固定压力，通过动量方程求解速度场。 然后，计算压力校正值， 并用于校正速度值使之满足质量守恒方程 (4) 有两种常用的校正方法:SIMPLE和PISO,压力的变化可用于强制满足质量守恒： 净流入质量通量 净流出质量通量 使单元压力增加 使单元压力减小,1, SIMPLE法(Patankar,Spalding,1972) SIMPLE Semi-Implicit Methodfor Pressure-Linked Equations 对于稳态情况： 第一步：用当前压力值求解动量方程：这样求得的速度一般是步满足质量守恒方程的。,第二步：导出压力校正公式：(i)使速度的变化与压力的变化联系起来：这里假设其他源项的解是精确的不需要校正。(ii) 进行SIMPLE近似，忽略,上述近似是合理的，因为我们只关心收敛解，最终所有校正值都将为零。 (iii)对压力控制体积应用质量守恒：即，净质量通量由当前速度 和校正速度 合成：净质量通量因此，,或者写成：作后的压力校正公式为：第三步：求解压力校正方程(与标量输运方程有相同的形式)：第四步：校正压力和速度,2, SIMPLE法的各种变型与改进SIMPLER (SIMPLE Revised Patankar, 1980),SIMPLEC (Van Doormaal and Raithby, 1984),SIMPLEX (Raithby and Schneider, 1988),3, PISO法(Issa, 1986) PISO Pressure Implicit with Splitting of Operators 这种方法最初是用于求解时间相关的非迭代压力校正法。每个时间步由下列三个阶段构成： 第一步：用上一步的压力值求解时间相关的动量方程。 第二步：用压力校正方程和前面的SIMPLE过程得到满足质量守恒的流场。 第三步：用第二校正步形成第二个与质量守恒一致的流场。,小结 1，每个动量分量都满足自己的标量输运方程：2，动量方程的特征： 非线性 耦合的 压力相关的所以，求解应是： 耦合的，迭代方法,3，压力梯度项会导致奇偶解耦。结决办法： 交错网格技术 Rhie-Chow插值4，广泛使用的两种压力校正方法： SIMPLE PISO,3 时间相关的方法,时间相关的标量输运方程单步法多步法CFD中时间步进法的使用 小结,时间相关的标量输运方程 对于任意控制体积，时间相关的标量输运方程可表示为： 单元体积内输运量总量的变化率通过单元表面的净通量源 其中， 单元内的总量质量x浓度 通量是指通过单元面的输运率 在前面的讨论中，已经讲述了通量和源项的离散，即 净通量源,我们首先考察一般的一阶常微分方程的数值解法，然后再将该方法推广到CFD中。 其中，F是时间t和 的任意函数。求解这类初值问题一般有两种主要的方法： 1,单步法：用前一个时步的值 2,多步法：用前几个时步的值,单步法 对于一阶微分方程：单步法问题就是：给定t(n-1)时刻的 值，计算t(n)时刻的 值。 公式中将用到下面的记号：t(n-1)用上标old表示已知的值t(n)用上标new表示待求的值。,将微分方程积分后得或这是一个精确算式，但由于平均导数Fav与解相关是未知的，需要估算！1,导数的简单估算法 导数的简单估算法是CFD中使用最为广泛的方法，包括向前差分(欧拉法)，向后差分(向后欧拉法)，以及中心差分法(Crank-Nicolson),向前差分法：取初始时步对应的导数作为平均导数Fav 优点：容易实现，为显式方法，即方程右边是已知的。 缺点：时间上只有一阶精度CFD中稳定性受时间步大小的限制。,向后差分法：取时步末对应的导数作为平均导数Fav 优点：在CFD中，没有时间步大小的限制。 缺点：时间上只有一阶精度为隐式方法。,中心差分法：取时步开始和结束对应的导数的平均值作为平均导数Fav 优点：时间上具有2阶精度。 缺点：隐式方法CFD中稳定性受时间步大小的限制。,2,改进的单步法 对于方程 的单步法，可通过连续估算平均梯度来改进算法。下面是两个重要的例子：改进欧拉法(2个估算)：,Runge-Kutta法(4个估算)对于标量输运量， Runge-Kutta法是工程中使用最为广泛的方法。但在CFD中，由于 和F一般都代表节点值的向量，因而计算导数F(即计算通量和源项)将非常昂贵，所以，主要的CFD还是采用前面介绍的三种简单的算法。,3, 单步法在CFD中的应用 一般的标量输运方程： 对于单步时间法，时间导数总是离散为： 而通量和源项可以在任意具体的时间水平上离散为下面的形式：,对上式采用不同的时间步进法：向前差分：重新整理后并略去上标new有评价： 显式方法，不需要求解方程组 受时间步长限制，为了稳定，要求 的系数为正定，即,向后差分：重新整理后并略去上标new有评价： 直接应用标量输运方程的离散格式，只需对系数作简单的修改： 隐式方法，没有时间步长限制。,Crank-Nicolson法：重新整理后并略去上标new有为方便，也可等式两边同时乘以2：,评价： 直接应用标量输运方程的离散格式，只需对系数作简单的修改： 受时间步长限制，为了稳定，要求 的系数为正定，即,上述几种时间步进法可写成如下的通用形式：其中， 为向前差分， 为向后差分， 为Crank-Nicolson法。 当 称作 法，当 (即除完全隐式向后差分法的任何其他算法)，稳定性要求时间步满足：,Courant数：Courant数定义为可理解为一个时间步内的行进距离 与网格间距的 比值。 对于完全显式的方法，要求c1,其意思是一个时间步内信息对流传播的距离不能超过一个网格间距。 对于Crank-Nicolson法， Courant数的限制可稍微放松一点。,多步法单步法只用到t(n-1)和t(n)时刻的值计算时间导数值 。 而多步法将用到更早时间步的 值，如t(n-2)用t(n3)，.等。 一个例子就是Gear法这是一个二阶精度的方法。,有一大类时间步进法是通过一个或多个校正来改进初值，这些方法称作预测校正法。其中，Adams-Bashforth-Moulton是最为流行的方法。 预测步 校正步,显然，多步法比单步法在时间上精度更高，但多步法存在许多缺点，使之在CFD中应用较少： 存储量大，每个计算变量在每个节点每个时步都要存储 起步困难，因为只有t=0时刻的值是已知的，其他启动步的值也需要单步法来计算。,CFD中时间步进法的使用 一般有两种情况需要使用时间相关的算法： (1) 一是本来就是时间相关的问题 (2) 另一种情况是稳态问题的时间步进法 其中，对于第一种情况，精度和稳定性要求对时间步有限制。由于所有节点值都以相同的时间长进行更新，因而， 是全局性的。 对于第二种情况，由于只是关心最终的结果，而中间时步的解并不重要，因此，对精度要求不高，只要选择稳定算法即可，一般选向后差分。,小结 时间相关的流动问题在时间上是一阶的，采用时间步进法求解 时间步进法可以是显式的(时间导数在时步开始式就已知)或隐式的(需要每个时步进行迭代) 常用的单步法包括向前差分(完全显式法)，向后差分(完全隐式法)，以及 Crank-Nicolson法(半隐式法),单步法一般通过改变矩阵系数就很容易实现，对向后差分法，只需这样修改： 唯一无条件稳定的方法是向后差分法，它是完全隐式法。其他方法均有时间步长的限制。一般来说由Courant 数决定：多步法可以改进精度和稳定性，但CFD中很少用到，因为需要巨大的存储空间。,