计算流体力学数值方法讲解

3-1计算流体力学第三章数值方法3-2计算流体力学主要内容空间离散技术标量输运方程动量方程时间离散技术3-3计算流体力学1标量输运方程有限体积法一维对流扩散方程扩散项的离散源项的离散代数方程的组装二维和三维问题对流项离散基础离散特性3-4计算流体力学对流项高级离散方法高阶对流方法的实现曲线网格边界条件代数方程的求解小结3-5计算流体力学对于任意控制体内一个守恒物理量:控制体内物理量的变化率+经过控制体边界面的净通量＝控制体积内的源(汇)经过控制体边界面的总的通量由对流(随流动的迁移)和扩散(由随机分子运动或湍流运动导致的净输运)两部分组成。如果用表示单位质量流体的守恒量，则通用标量输运方程(或称对流－扩散方程)可表示为：3-6计算流体力学1，有限体积法直接对上式进行离散2，本章只考虑稳态问题，即上式左边第一项为零SVnACVdtdfaces)(）（AuCn3-7计算流体力学有限体积法(FVM)(1)定义流场求解域几何形状(2)将求解域划分为计算网格，即一组互不重叠的有限体或单元。(3)基于上述划分的单元对积分方程进行离散，即用节点值来近似。(4)对得到的离散方程进行数值求解。3-8计算流体力学计算网格可以是结构网格或非结构网格，笛卡儿网格或非笛卡儿网格。最常用的网格形式包括基于单元中心的存储方式和基于单元顶点的存储方式两种。在后面的讨论中将会看到，不是所有的变量都必须存储在同样的位置。3-9计算流体力学为了简单起见，本课程将只考虑结构化笛卡儿网格，和采用基于单元中心的存储方式，即单元变量存储在单元中心节点位置上。右图为一典型三维控制体积，其中P点为单元中心存储节点。相对于中心节点，沿坐标方向通常表示为west,east,south,north,bottom,top.其中，小写字母w,e,s,n,b,t表示单元面，大写字母W,E,S,N,B,T表示中心节点的相邻节点。3-10计算流体力学于是单元各面的面积表示为Aw,Ae,As,AnAb,At。体积为V。对于二维问题，可以视为单位厚度为1的一层单元()对于结构网格，可以交替使用ijk下标表示单元节点，譬如，3-11计算流体力学一维对流－扩散方程我们首先讨论一维稳态对流扩散方程主要基于下面的考虑：(1)它使问题分析大大简化(2)离散方程可以进行手算。(3)尽管只是一维的，但要扩展到二维或三维是非常直接的(4)实际上，通量(对流和扩散)的离散一般是沿坐标方向进行的，即分别沿i,j,k线进行。(5)有许多重要的理论问题是一维的。3-12计算流体力学对于如图所示的一维控制体，物理量的守恒可表述为如下关系式：[通量]e(fluxe)-[通量]w(fluxw)＝源(source)这里的通量是指穿过单元表面的输运率。如果表示单位质量的输运量，则总的通量为对流通量和扩散通量之和，其中：对流通量＝扩散通量＝3-13计算流体力学于是，的对流扩散方程为：其中，S表示单位长度的源。将等式两边除以，并取极限，可得到相应的微分方程：注意：1，这里面积A是表示可变截面的准一维问题，对于真正的一维问题，只需令A=1,即2，这里假设和S均为常数，但在一般的CFD问题中，u本身也是问题的解变量。3-14计算流体力学例子1，纯扩散问题如图所示的隔热棒，长度1m，截面为1cmx1cm的方形截面，棒的两端为固定温度，分别为100度和500度。穿过任意截面A的热通量由下式给定：其中，热传导系数(a)将棒划分为5个控制体，并用有限体积分析沿棒的温度分布(b)写出沿棒温度分布的微分方程(c)求出微分方程的解析解，并与(a)的解进行比较3-15计算流体力学这个问题是求解一维对流－扩散方程：上面的方程也可写成如下的积分方程：由于本例只考虑扩散，即没有对流和源项：3-16计算流体力学对于本例的温度T的纯扩撒问题，最终有如下的微分方程和积分方程：00)(ewdxdTkAdxdTkAdxd3-17计算流体力学非边界面上通量的计算：3-18计算流体力学边界面上通量的计算：3-19计算流体力学02/00002/52454534342323121211xTTxTTxTTxTTxTTxTTxTTxTTxTTxTTBB3-20计算流体力学1000302020220035454343232121TTTTTTTTTTTTT46038030022014014311TTTTT3-21计算流体力学解析解：100400500)1(100)0(00)(0122xTxTxTcxcTdxTddxdTkAdxd3-22计算流体力学控制方程扩散项的离散梯度扩散项的离散几乎总是采用中心差分格式：提示：(1)有限体积法中，通量是在面上计算的，而不是在节点处。(2)上述对梯度扩散项的近似在空间上具有二阶精度，后面将会给出证明。3-23计算流体力学(3)如果扩散系数也是变量的话，它在单元面上的值必须通过插值得到。(4)采用中心差分格式意味着在单元面两边的节点上权重相同。这和扩撒的物理意义是一致的，因为扩散在所有方向作用相同。这后面将要讨论的对流是不一样的，对流是具有明确的方向性的。3-24计算流体力学(5)单元西面也应有相似的表达式，这是守恒定律所要求的。即从一个单元流出的通量一定等于流入相邻单元的通量。这是有限体积法优越有限差分和有限元法的地方。在有限体积法中，常常定义一个扩散输运系数D，譬如，对于单元东面有：3-25计算流体力学控制方程源项的离散当源与流体的量成比例的情况下，每个单元的总的源强度可简单的表示为：单位体积的源x单元体积＝SxV其中，S表示平均源密度，一般取单元中心节点处的值。一般来说，S常常是与变量相关的，譬如，如果变量为温度的话，牛顿冷却定律认为：单位时间热量的损失与相对与环境的温差成比例：一般将源项分解为与解相关和无关两部分：3-26计算流体力学代数方程的组装当速度u=0时，定常扩散问题在每个单元上的离散为：经过整理后：或表示成如下标量输运方程规范形式：3-27计算流体力学上述方程为离散后的标量输运方程规范形式，每个求解变量的输运方程都有相同的形式，只是各自的矩阵系数不同。对于纯扩散问题，矩阵系数为：3-28计算流体力学上述方程为只是一个控制体的离散方程，如果将所有控制体单元中心节点处的变量值组装到一个向量中，就会得到如下形式的代数方程：求解这种代数方程可采用三对角矩阵算法3-29计算流体力学二维和三维的扩展前面讨论的虽然只是一维问题，但要推广到二维和三维是非常直接的。对于多维问题，经过单元面的净通量只需将所有面上的通量求代数和即可！譬如，对于三维控制体，其守恒方程可表示成：3-30计算流体力学对于多维问题，离散后的代数方程仍然具有如下形式：上述方程只是对于单个控制体单元的，而对所有单元进行组装便得到一组节点值向量的矩阵方程。对于二维问题，最终得到的矩阵具有如图所示的带状形式。三维问题将增加非对角元素！3-31计算流体力学对流项的离散(I)纯扩散问题(u=0)只是一种理想情况，对于典型的工程问题，对流项远远超过扩散项，因为雷诺数非常大，即惯性力和粘性力的比值很大。定常的对流－扩散方程方程：3-32计算流体力学为了简化书写，将质量通量用C表示，譬如，单元e面的质量通量表示为：扩散项和源项的离散采用与前面相同的方法，于是有：现在的问题是如何离散对流项？也就是怎样利用相邻节点值来确定单元面上的变量值。这类方法就称作对流方案。3-33计算流体力学中心差分法离散对流项的最明然的方法就是线性插值，即中心差分法：于是将上述公式代入对流扩散方程后，得到关于单个控制体单元如下形式的离散方程：3-34计算流体力学上式中的求和是对中心点的所有相邻节点进行的。各个系数表达式如下：上式中，根据质量守恒有3-35计算流体力学为了说明中心差分格式是否适合离散对流项，下面给出了定常对流－扩散问题(无源项)在下面两种情况下的解的图示：Pe=1/2pe=4其中Pe是Peclet的缩写，称作Pe数。3-36计算流体力学上图中红圈表示计算值，实线为解析解。第一重情况吻合很好，但第二种情况却有显著的摆动，这是为什么呢？3-37计算流体力学从数学的角度来说，当Pe数大于2，系数矩阵的一个元素将是负值，即相应于该元素的节点值的增加会使中心节点值减小！从物理的角度来看，是因为对流过程是有方向性的。而输运属性仅仅在流动方向。但中心差分格式没有反映处方向性，而是在上风(upwind)节点和下风(downwind)节点上赋予了相同的权重。3-38计算流体力学上风差分法在上风差分格式种，单元面上变量的值用其上风节点的值近似。譬如，对于单元东(e)面:于是将上述公式可归结为：3-39计算流体力学采用上风差分格式后，定常对流－扩散问题的规范离散格式为：其中：3-40计算流体力学当将上风差分法应用于前述与中心差分同样的问题后可观察到：(1)当所得结果精度不如中心差分法，这是意料之中的。(2)当Pe=4时，得到的解不是很精确，但不再出现明显的摆动！于是，便出现了精度和有界性(无摆动)之间如何取舍的问题？下面将介绍一些离散方程的属性。3-41计算流体力学离散属性1，一致性一致性是指当网格间距趋向于0时，离散方程等价于连续方程。譬如，就是的一致近似。2，守恒性即从一个单元流出的通量一定等于流入相邻单元的通量。这是有限体积法优越有限差分和有限元法的地方。3-42计算流体力学3，输运性方向性的影响反应在对流项的离散方法中。实际上，就是在上风节点赋予较高的权重。4，有界性在无源的对流－扩散问题中的解一定介于周围节点流动变量最大值和最小值之间。5,稳定性稳定性决定了求解过程能否得到最终的解。它不涉及解的精度。稳定即意味着误差不随求解过程增加。3-43计算流体力学6，阶(解的精度)阶是指随着空间网格尺寸的减小，误差衰减的速度。譬如，对于均匀空间网格，在某种数值方法中，如果当空间网格尺寸趋向于0，而误差与空间网格尺寸的n次方成正比，则称这种数值方法是n阶的。某种数值方法的阶可以通过泰勒级数展开式中的首阶误差项确定(后面将讨论)。注意：此处的误差是指展开式中的理论截断误差，而不考虑计算机的舍入误差(即计算机只能存储有限位数的数)3-44计算流体力学高阶精度可通过采用更多的节点值近似来获得。一个节点允许的最高精度为1阶，两个节点允许的最高精度为2阶，依此类推。理论上讲，某种数值方法的精度越高，随着网格的加密，误差减小的就越大。也就是说，采用高精度的数值方法，只需较少的网格数即可获得要求的精度。但是，高阶精度的方法常常需要更多的计算时间，而且常常会导致解的有界性问题。3-45计算流体力学前面我们讲述过的对流和扩散项的差分离散格式的阶可以通过对单元面两边的节点值在单元面上按泰勒级数展开来获得：3-46计算流体力学将上述两式相减即得：所以，当取为误差项时，用近似的精度是二阶的精度。3-47计算流体力学另外，将前面两式相加即得：所以，当取为误差项时，用中心差分格式近似的精度是二阶的精度。而如果采用上风差分近似或(根据流动方向来定)则是一阶精度。3-48计算流体力学对流项高级离散方法在了解了前面对流项离散的基础方法的前提下，下面将进一步介绍一些常用的高级离散方法发。1，混合法(Splading,1972)这种方法的思路是：当|Pe|小于或等于2时，采用中心差分格式离散对流项，当|Pe|大于2时，采用上风差分格式离散对流项譬如，对于u0,3-49计算流体力学其中，对每个面(e或w)：3-50计算流体力学对混合法的评价：a)这种方法是满足对流离散特性中的守恒性、输运性和有界性的。b)由于这种方法的非常稳定和健壮，使其在过去很长一段时间内成为商用CFD程序中最流行的方法。c)但由于实际工程中的流动问题常常是高对流低扩散的，因而，混合法实际上主要还是一阶的上风差分格式。所以，现代CFD试图寻求更高精度的对流离