基本方程组的数值求解3.ppt

3 时间导数项的差分 在时间坐标上 一个显著的特点是只有一阶导数项 没有高阶导数项 考虑定常二维边界层方程其中在方向也只有一阶导数 可称为类时间坐标此处的讨论同样适用于有类时间坐标的问题 时间导数项在控制容积上积分 得 45 这里假定 在整个控制容积中是均匀的 网格不随时间而变在时间导数项的差分式子中 包含了两个时间层 在类时间坐标上是两个截面 的参数值与此相对应 其它坐标方向的对流项 扩散项以及源项中的参数 在时间层上可以有不同的取法 因而可以得到显式 隐式和克兰克 尼科尔森 Crank Nicholson 格式 1 显式差分格式 46 2 隐式差分格式 47 3 克兰克 尼科尔森格式若对流项 扩散项和源项中的参数都取两个时间层上的算术平均值 可得差分方程 48 克兰克 尼科尔森格式的差分方程也需要与其它结点的差分方程联立求解 当时间步长较小时 其精度比隐式格式高 4 源项的差分 设源项在控制容积内均匀 并等于中心点的值 则有 49 微分方程离散化后得到的是一个拟线性代数方程 但源项往往是参数的非线性函数为能够用解线性代数方程的方法求解差分方程 需要将源项线性化源项线性化有两种办法 一是用前一时刻的参数代入源项 算得的值作为常数参加下一时刻 或下一循环 的计算 二是假定源项与参数的函数关系可以近似用式子 50 表示 代入差分方程得 51 其中Sc和Sp可能也是的函数 因此式 50 也是拟线性关系 第二种办法优于第一种办法为了使拟线性代数方程组收敛 得到有意义的解 Sc和Sp必须满足 52 在满足条件 52 的基础上 和的确定有一定的任意性 一个方便的办法是用切线使其线性化 即令 53 这时 按照切线法得到的和往往很复杂 而且不一定满足条件 52 比较多的是用经验办法来确定 五 代表方程组的求解 1 概述采用有限容积法在已生成的网格上将所求解变量的控制方程离散后 就形成了各求解变量的代数方程组求解离散化所得的代数方程组是流动与燃烧过程数值计算的最后一个环节求解的方法有两大类 直接解法与迭代法迭代法的研究主要是研究收敛性以及如何加快收敛速度求解代数方程的解法包括Gauss消元法 三对角阵 TDMA 算法 五对角阵 PDMA 算法等现有文献中求解由QUICK等高阶格式形成的代数方程组时多采用TDMA方法求解 TDMA算法是Gauss消元法的一种特例 每一行上仅三个非零元素 2 Gauss消元法3 TDMA算法对图1中S N方向的网格线进行计算时 差分方程需重新整理为 56 方程右端各项是邻近网格线上的结点值 或源项 取前一循环得到的值 故为已知值每个节点的代数方程中最多只包含三个节点的未知值 可以认为其它节点上未知值的系数均为零如果把上述有限差分离散方程组写成矩阵的形式 其系数阵是一个三对角阵 仅对角元素及其上下邻位上的元素不为零 而其它元素均为零把Gauss消元法应用于这种情形 便构成了称为三对角阵算法的有效求解方法 简记为TDMA TridiagonalMatrixAlgorithm 57 对于这种系数矩阵为三对角线矩阵的方程组 TDMA求解的具体步骤是 1 第一个方程的各项除以B1 得 第二个方程减去化简后的第一个方程乘以A2 并将所得新方程的各项除以 最后得对第三个以后的每个方程都作同样的处理 得到方程组系数的通式为 58 2 自下而上的解方程 解的通式为 59 若有两个或三个坐标是椭圆型的 则需对各椭圆型方向都进行逐线迭代扫描对两个或三个方向各作了一次逐线计算 称为进行了一次双重或三重扫描 也叫一个迭代循环重复多重扫描 直至两次相继迭代得出的值差别不大为止 把最后结果作为n时间层的值对每条网格线进行计算时 边界条件很快传入流场内部 因此收敛速率比逐点迭代快得多 隐式差分方程的逐线迭代解法对时间步长也无限制 因此在实际中得到广泛的应用 六 压力和速度之间的耦合 1 引言若知道流场中压力的分布 则动量方程的求解与其它方程完全一样 不会产生困难但是压力值一般都是预先不知道的在可压流中 压力与密度间的关系由状态方程确定 需通过连续方程和状态方程确定在不可压流中 流场中的压力分布对速度场有很强的影响 压力梯度以源项形式出现在动量方程中 但压力却没有独立的控制方程在计算流体力学的发展过程中 提出了多种不同的方法 来解决在以速度 压力为求解变量的原始变量法中的这个问题大致可分为 涡量 流函数方法 压力修正方法 压力与速度之间的迭代算法等几种 以压力修正方法中的SIMPLE系列方法应用最为广泛 2 压力修正方法的基本思想 在求解不可压缩流体的流场问题时 如果我们把从动量方程与连续性方程离散得到的代数方程组联立起来直接求解 就可以得到各速度分量及相应的压力值这样的直接解法要占用大量的计算机内存 对于目前大多数工程应用场合还不适用如果采用分离式的迭代求解方法 segregatedmethod 即先求解u速度场 再求解v速度场 则对压力场因其无独立的方程而无法对其求解 或无法改进其原先的假定值 另一方面 上述分离式求解过程中只利用了u v动量方程的离散形式而未用到连续性方程要解决的问题 如何使得用质量守恒方程假定的压力场能不断地随迭代过程的进行而得到改进 这就是所谓的压力修正算法压力修正算法源于1972年由Patankar与Spalding提出的SIMPLE算法 16 在20世纪80年代初期 又相继提出了SIMPLER与SIMPLEC等方法 并由此形成了SIMPLE系列算法 3 SIMPLE算法 以二维直角坐标中的对流换热求解为例 其控制方程为 60 61 62 在交错网格 图3 上动量离散方程为 63 64 质量守恒方程为 65 采用SIMPLE算法来求解时计算步骤如下 1 假定一个速度场 记为u 0 v 0 由此计算ae an anb及b 2 假定一个压力场 记为 3 求解动量离散方程 63 64 得 4 计算压力的修正值 要求与相对应的 能满足连续性方程 5 计算速度的修正值 要求仍满足线性化了的动量方程 即与式 63 相减 得 69 此式表明 要据值确定 需要解一个代数方程组 为能利用值显式地求解 此处略去式 69 右端第一项 于是得 70 6 将 作为本迭代层次之解 开始下一层次的计算 重复1 5步直到流场收敛 即所解得之流速能同时满足连续性方程及动量方程为止 SIMPLE算法中引入了以下假定或简化处理 1 速度场的假定 u 0 v 0 与压力场的假定 是相互独立地进行的 u 0 v 0 与间没有任何联系 2 在导出速度修正值计算式 70 时 未计及邻点速度修正值的影响 3 动量离散方程中的b 控制容积P的剩余质量流量 在速度修正前后保持不变 4 由式 63 64 解出的满足动量守恒但未必满足质量守恒 而由式 69 解出的决定时 保证 满足质量守恒 但动量守恒则未必满足 注意 虽然在导出式 69 过程中曾要求 满足线性化了的动量方程 但由于以后略去了式 69 中的 因而所得的 就未必可使 满足上述动量守恒方程 以上四条假设或近似处理 是SIMPLE算法提出之后所出现的一些改进方案的着眼之处 4 SIMPLER算法 SIMPLER 1980 算法主要用以改进SIMPLE方法中的第一项近似处理方法 一旦速度场给定 压力场就可以从动量离散方程中予以求解 而不再任意假定 其主要计算步骤如下 1 假定一个速度场 记为u 0 v 0 由此计算动量离散方程系数ae an anb b及 71 注意 引入 后 动量离散方程便可写为 72 2 据计算相应的压力场 为此 将式 72 代入质量守恒方程的离散形式 65 得 73 其中 anb及b的计算式形式上与SIMPLE算法中的方程的一样 只要将代替即可 3 求解动量离散方程 得替 4 求解压力修正方程 得 5 用修正速度得 但不修正压力 6 用 开始下一层次的迭代 重复1 5步直到收敛 为有利于非线性迭代的收敛 SIMPLER算法中流速应予以亚松弛 但及p则不作亚松弛 5 SIMPLEC算法6 SIMPLEX算法7 SIMPLET算法 8 SIMPLE系列算法小结 SIMPLE SIMPLEC SIMPLER SIMPLET及SIMPLEX的主要计算步骤可以归纳为 1 由假定的或上一次计算得到的速度u 0 v 0 确定动量离散方程的系数及常数项 2 确定压力场 对于SIMPLE SIMPLEC SIMPLET及SIMPLEX采用假定的或用上一次的计算值 对于SIMPLER则由已知的速度场通过求解压力Poisson方程而获得 3 求解动量离散方程以获得 4 求解压力修正值方程 要求与相应的能使 满足连续性方程 其中除SIMPLET外 其它算法只考虑压力修正值对的影响 而SIMPLET中还考虑了温度修正对的影响 该影响最后体现在方程的源项之中 5 确定速度修正值 SIMPLE SIMPLER SIMPLET 邻点速度影响略去不计 SIMPLEC 邻点速度修正值的影响已基本考虑在内 SIMPLEX de等是通过求解一个Poisson方程得出的 相当于考虑了邻点的影响 6 采用下列速度 压力值作为初值开始下一层次的迭代计算 设速度的亚松弛已组合在u v的求解过程中 SIMPLE SIMPLET SIMPLEC SIMPLEX SIMLER 以上计算过程可以大致分为两步 预估步 1 3 及校正步 4 6 因而可以说上述五种算法都是两步算法 一步预估 一步校正 9 加速迭代收敛的一些方法 1 选择合适的松驰因子之值对于正交的网格以及虽非正交但网格倾斜不是严重的情况 建议au与ap的取值应满足以下关系 76 式中常数c取为1 23 或1 1 24 同时au之值应尽可能地取得大 一般可取0 7 0 8 上述取值原则对交错网格 同位网格 正交网格与非正交网格原则上都适用 但对网格严重倾斜的情况 ap之值应小于按式 76 得出之值 2 其它加速迭代收敛的方法 迭代循环方式对收敛速率的影响描述燃烧过程的是一组互相耦合的方程组 在用迭代法求解时 迭代循环方式对收敛速率有很大的影响可将流体力学方程 动量和修正压力方程 和其它标量方程 能量 湍流量方程等 一起迭代另一种选择是将流体力学方程首先迭代若干次 然后把其它标量方程 能量 湍流量方程等 一起迭代 流体力学方程首先迭代若干次 然后把其它标量方程再加进去一起循环两者的收敛速率是不一样的 但目前还给不出最佳循环方式的一般原则 网格均匀度对收敛速率的影响网格均匀度对收敛速率有很大的影响计算时间与结点数成正比但增加结点数 若能使差分网格更均匀 由于收敛速率加快 有时仍有可能减少总的计算时间 七 边界条件的数值处理方法 1 引言初始条件和边界条件 对抛物型问题 只需给出一端的边界条件 如时间坐标上给出初始的参数分布 边界层流动中给定进口截面上的参数分布等等对椭圆型问题 则需给出求解域两端的边界条件 如果三个空间坐标都是椭圆型坐标 则求解域的四周各点都需给出边界条件 2 进口边界 模型计算湍流时 一般无法获得由实验测得的进口处k及k 1 2 1 2 进口截面上k与取值的这种不确定性对计算结果的影响 亦即计算结果对进口条件的敏感性 应当在数值计算中进行考察 3 出口边界 描写有回流的流动与燃烧的控制方程是关于空间坐标的椭圆型方程 从数学上要求在每个坐标方向上都应当有两个边界条件 出口边界应该置在什么位置及怎样确定这些人为设定的计算边界条件 就成了计算流体力学及燃烧过程模拟中的一个重要研究课题 出口边界应设置在什么位置上 Patankar曾认为 出口边界不能设置在有回流的截面上 否则得到的解是没有意义的数值实践证明 对出口边界位置选择的这种苛刻的要求是不必要的文献 27 提出对这种出口边界条件定性的要求是 1 应当对出口边界附近 尤其是离开出口边界较远处的流场没有影响 2 出口边界应当是 透明的 又称为 无反射 3 无论出口边界设置在何处 应对其它区域中的数值解没有影响 图11求解流场时出口边界位置的确定 目前 文献中已提出了多种关于设置出口边界条件的方法 归纳起来大约有以下几种情形 1 局部单向化假设只有当出口截面无回流时才能近似地成立2 充分发展的假设假设n为出口截面的法线方向 这一条件的数学表达式为 84 3 法向速度局部质量守恒 切向速度齐次Neumann条件以图11 b 中出口边界相邻的控制容积为例 取出后放大而示于图12中 在本例中速度u是与出口截面平行的流速 而v为与出口截面垂直的流速 对所示控制容积作质量守恒运算 得 85a 而按齐次Neumann条件有 85b 图12局部质量守恒的图示图13说明切向速度按局部质量守恒的计算图示 4 出口截面上的法向速度与切向速度都按局部质量守恒确定 在这一方法中 对图12所示情形 计算式同式 85a 但ui M的更新则通过对图13所示半个控制容积的质量守恒来获得 以0 5 vi M vi M 1 作为该半个控制容积南侧流速 则可以写出 86a 将边界上的切向流速取本迭代层次之值而其余速度取上一迭代层次的值 得 86b 由图11 b 可见 在计算区域右边界 之值为零 由此可推得 这样依次计算可以得出出口边界切向速度的更新值 5 对流边界条件 用以下公式来更新出口截面法向速度及其它标量 87a 其中 c为对流速度或相速度 一般取为通道平均流速 上式离散后 可得 87b 这里 vm为出口截面平均法向流速 相当于假定在出口边界上无扩散作用 为一个纯对流问题 6 内点插值方法 在文献 36 中 提出采用由内点的速度来决定出口流速的方法 88a 其中 为与出口边界相邻的控制容积的宽度 x为出口边界上游第四个截面位置起算的坐标 当网格均分时 如图14所示 有 88b 图14插值公式 88a 的图示 7 各种高阶导数为零的方法 89 其中 x为出口方向 上角标j及q表示求导的阶数 要求 并建议取j 3 q 2 这一类方法实际上相当于提出了由内点来插值的一个方式 导数的阶数越高 涉及的内点数越多 该阶导数取零对解的影响也越微弱 出口边界上边界条件的设定有多种看来相差很大的方法 在从事工程流动与传热问题的数值计算中实用的处理出口边界条件的方法 出口边界上无回流时 可采用局部单向化或法向一阶导数为零的方法 当出口边界有回流时 可以采用法向流速按质量守恒 切向流速按齐次Neumann条件的方法 无论采用哪一种方法 出口截面上的法向速度分布必须满足总体质量守恒的条件 4 固体边界 1 对于温度一类的变量用附加源项法来处理2 对于速度固体壁面是受限流动中最常见的边界 原则上应采用无滑移条件 即流体的速度等于壁面的速度 第一类边界条件 层流时无论是对于动量方程中的流速还是连续性方程中的流速都是这样 后者表现在SIMPLE算法压力修正方程中的剩余质量流量的计算上 当采用标准模型时 平行于壁面的流速及垂直于壁面的流速还存在不同处理方法 对于平行于壁面的流速 例如u 采用无滑移条件 同时壁面上的动力粘性系数由壁面函数法的公式来计算 但对垂直于该壁面的流速 v 则应区分动量方程还是连续方程 在连续性方程中仍然采用无滑移条件 静止壁面上取v 0 而对v动量方程则采用法向一阶导数为零的条件 3 对于k采用低Reynolds数模型时 可直接令壁面上k 0采用标准模型时一般取第一个内节点处 可通过令该点与壁面之间的k方程的扩散系数为