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

1、3 3时间导数项的差分时间导数项的差分在时间坐标上，一个显著的特点是只有一阶导数项，没有高阶导数项，考虑定常二维边界层方程其中在 方向也只有一阶导数，可称 为类时间坐标此处的讨论同样适用于有类时间坐标的问题。yyyuxxx时间导数项 在控制容积 上积分，得（45）这里假定 、在整个控制容积中是均匀的，网格不随时间而变在时间导数项的差分式子中，包含了两个时间层（在类时间坐标上是两个截面）的参数值与此相对应，其它坐标方向的对流项、扩散项以及源项中的参数，在时间层上可以有不同的取法，因而可以得到显式、隐式和克兰克-尼科尔森（Crank-Nicholson）格式t 11nPnnPnttdt1 1）显式

2、差分格式）显式差分格式（46）2 2）隐式差分格式）隐式差分格式（47）3 3）克兰克）克兰克-尼科尔森格式尼科尔森格式 若对流项、扩散项和源项中的参数都取两个时间层上的算术平均值，可得差分方程（48）克兰克-尼科尔森格式的差分方程也需要与其它结点的差分方程联立求解。当时间步长较小时，其精度比隐式格式高。nbnnbnbnpPnnPnBBtt1111111nPnnbnnbnbnPnPtBtBnbnnbnbnPpnnbnnbnbnPpnBBtBBt1111224 4、源项的差分、源项的差分 设源项在控制容积内均匀，并等于中心点的值，则有（49）微分方程离散化后得到的是一个拟线性代数方程，但源项往往

3、是参数的非线性函数为能够用解线性代数方程的方法求解差分方程，需要将源项线性化源项线性化有两种办法:一是用前一时刻的参数代入源项，算得的值作为常数参加下一时刻（或下一循环）的计算；二是假定源项与参数 的函数关系可以近似用式子（50）表示 SSd pcSSS 代入差分方程得（51）其中Sc和Sp可能也是 的函数，因此式（50）也是拟线性关系。第二种办法优于第一种办法为了使拟线性代数方程组收敛，得到有意义的解，Sc和Sp必须满足 ,（52）111npnbncnnbnbnpppntSBSBt0PS0cS 在满足条件（52）的基础上，和 的确定有一定的任意性。一个方便的办法是用切线使其线性化，即令（53

4、）这时 ，按照切线法得到的 和 往往很复杂，而且不一定满足条件（52）比较多的是用经验办法来确定cSPS 000SSS0000SSS000SSSc0SSpcSPS五、代表方程组的求解五、代表方程组的求解1 1、概述、概述采用有限容积法在已生成的网格上将所求解变量的控制方程离散后，就形成了各求解变量的代数方程组求解离散化所得的代数方程组是流动与燃烧过程数值计算的最后一个环节求解的方法有两大类：直接解法与迭代法迭代法的研究主要是研究收敛性以及如何加快收敛速度求解代数方程的解法包括Gauss消元法、三对角阵（TDMA）算法、五对角阵（PDMA）算法等现有文献中求解由QUICK等高阶格式形成的代数方程

5、组时多采用TDMA方法求解，TDMA算法是Gauss消元法的一种特例（每一行上仅三个非零元素）2.Gauss消元法3.TDMA算法 对图1中S-N方向的网格线进行计算时，差分方程需重新整理为（56）方程右端各项是邻近网格线上的结点值（或源项），取前一循环得到的值，故为已知值每个节点的代数方程中最多只包含三个节点的未知值，可以认为其它节点上未知值的系数均为零如果把上述有限差分离散方程组写成矩阵的形式，其系数阵是一个三对角阵 仅对角元素及其上下邻位上的元素不为零，而其它元素均为零把Gauss消元法应用于这种情形，便构成了称为三对角阵算法的有效求解方法，简记为TDMA（Tridiagonal Mat

6、rix Algorithm）bBBBBSNnbnbnbSSPPNN,（57）MMMMMMMMMDDDDDBAOOCBAOOCBAOOOCB1321132111132211对于这种系数矩阵为三对角线矩阵的方程组，对于这种系数矩阵为三对角线矩阵的方程组，TDMA求解的具体步骤是：求解的具体步骤是：（1）第一个方程的各项除以B1，得 ，。第二个方程减去化简后的第一个方程乘以A2，并将所得新方程的各项除以（），最后得 对第三个以后的每个方程都作同样的处理，得到方程组系数的通式为：（58）11B111BCC 111BDD 212ACB2122222,1,0ACBCCBA mmmmmmmmmmmmmmAC

7、BDADDACBCCBA11110（2）自下而上的解方程 解的通式为（59）若有两个或三个坐标是椭圆型的，则需对各椭圆型方向都进行逐线迭代扫描对两个或三个方向各作了一次逐线计算，称为进行了一次双重或三重扫描，也叫一个迭代循环重复多重扫描，直至两次相继迭代得出的值差别不大为止。把最后结果作为n时间层的值对每条网格线进行计算时，边界条件很快传入流场内部，因此收敛速率比逐点迭代快得多。隐式差分方程的逐线迭代解法对时间步长也无限制，因此在实际中得到广泛的应用MMMMMMCDD1111mmmmCD六、压力和速度之间的耦合六、压力和速度之间的耦合1、引言、引言若知道流场中压力的分布，则动量方程的求解与其它

8、方程完全一样，不会产生困难但是压力值一般都是预先不知道的在可压流中，压力与密度间的关系由状态方程确定，需通过连续方程和状态方程确定在不可压流中，流场中的压力分布对速度场有很强的影响，压力梯度以源项形式出现在动量方程中，但压力却没有独立的控制方程在计算流体力学的发展过程中，提出了多种不同的方法，来解决在以速度、压力为求解变量的原始变量法中的这个问题大致可分为：涡量-流函数方法、压力修正方法、压力与速度之间的迭代算法等几种。以压力修正方法中的SIMPLE系列方法应用最为广泛2压力修正方法的基本思想压力修正方法的基本思想在求解不可压缩流体的流场问题时，如果我们把从动量方程与连续性方程离散得到的代数方

9、程组联立起来直接求解，就可以得到各速度分量及相应的压力值这样的直接解法要占用大量的计算机内存，对于目前大多数工程应用场合还不适用如果采用分离式的迭代求解方法（segregated method），），即先求解u速度场，再求解v速度场，则对压力场因其无独立的方程而无法对其求解（或无法改进其原先的假定值）另一方面，上述分离式求解过程中只利用了u,v动量方程的离散形式而未用到连续性方程要解决的问题:如何使得用质量守恒方程假定的压力场能不断地随迭代过程的进行而得到改进，这就是所谓的压力修正算法压力修正算法源于1972年由Patankar与Spalding提出的SIMPLE算法16，在20世纪80年代初

10、期，又相继提出了SIMPLER与SIMPLEC等方法，并由此形成了SIMPLE系列算法3SIMPLE算法算法l以二维直角坐标中的对流换热求解为例，其控制方程为（60）（61）（62）l在交错网格（图3）上动量离散方程为：（63）（64）l质量守恒方程为（65）2,10ixutiiijijijijiBxuxxpxuutujpjjjxTckxxTutTEPenbnbeePPAbuauau:NPnnbnbnnPPAbvavav:00nsneweppAvvAuutyx采用采用SIMPLE算法来求解时计算步骤如下：算法来求解时计算步骤如下：l（1）假定一个速度场，记为u(0)，v(0)，由此计算ae，a

11、n，anb及 b；l（2）假定一个压力场，记为 ；l（3）求解动量离散方程（63）、（64），得 ；l（4）计算压力的修正值 ，要求与 相对应的 、能满足连续性方程pvu,pppuuvv（5）计算速度的修正值 。要求 仍满足线性化了的动量方程，即 与式（63）相减，得（69）此式表明，要据 值确定 ，需要解一个代数方程组。为能利用 值显式地求解 ，此处略去式（69）右端第一项，于是得（70）（6）将 、作为本迭代层次之解，开始下一层次的计算，重复15步直到流场收敛，即所解得之流速能同时满足连续性方程及动量方程为止 EEPPenbnbnbeeeppppAbuuauuaEPenbnbeeppAua

12、uaEPeEPeeeppdppaAuvu,nneevvuu及pupnevu,nbnbaauu papvvP及SIMPLE算法中引入了以下假定或简化处理：算法中引入了以下假定或简化处理：（1）速度场的假定（u（0），v（0）与压力场的假定（）是相互独立地进行的，u(0)，v(0)与 间没有任何联系；（2）在导出速度修正值计算式（70）时，未计及邻点速度修正值的影响；（3）动量离散方程中的b(控制容积P的剩余质量流量)在速度修正前后保持不变；（4）由式（63）、（64）解出的 满足动量守恒但未必满足质量守恒，而由式（69）解出的 决定 时，保证 、满足质量守恒，但动量守恒则未必满足。注意，虽然在导

13、出式（69）过程中曾要求 、满足线性化了的动量方程，但由于以后略去了式（69）中的 ，因而所得的 ，就未必可使 、满足上述动量守恒方程。l以上四条假设或近似处理，是SIMPLESIMPLE算法提出之后所出现的一些改进方案的着眼之处。ppvu,pvu,vuuvuuuvnbnbaau vuuvv4、SIMPLER算法算法 lSIMPLER（1980）算法主要用以改进SIMPLE方法中的第一项近似处理方法。一旦速度场给定，压力场就可以从动量离散方程中予以求解，而不再任意假定，其主要计算步骤如下：（1）假定一个速度场，记为u(0)，v(0)，由此计算动量离散方程系数ae，an，anb，b及 ，：（71

14、）注意，引入 ，后，动量离散方程便可写为（72）u v nnbnbnenbnbeabvavabuau)0()0(eu nv NPnnnEPeeeppdvuppduu（2）据 计算相应的压力场 。为此，将式（72）代入质量守恒方程的离散形式（65），得（73）其中，anb及b的计算式形式上与SIMPLE算法中的方程的一样，只要将 代替 即可。（3）求解动量离散方程，得替 。（4）求解压力修正方程，得 。（5）用 修正速度得 （但不修正压力）。（6）用 ，开始下一层次的迭代，重复15步直到收敛。为有利于非线性迭代的收敛，SIMPLER算法中流速应予以亚松弛，但 及p则不作亚松弛。vu,pbpapa

15、nbnbPPvu,vu,vu,ppvu,uuvvp5、SIMPLEC算法算法6SIMPLEX算法算法7、SIMPLET算法算法8SIMPLE系列算法小结系列算法小结lSIMPLE，SIMPLEC，SIMPLER，SIMPLET及SIMPLEX的主要计算步骤可以归纳为：（1）由假定的或上一次计算得到的速度u(0)，v(0)确定动量离散方程的系数及常数项；（2）确定压力场：对于SIMPLE，SIMPLEC，SIMPLET及SIMPLEX采用假定的或用上一次的计算值 ；对于SIMPLER则由已知的速度场通过求解压力Poisson方程而获得；（3）求解动量离散方程以获得 ；（4）求解压力修正值 方程，

16、要求与 相应的 能 使 、满足连续性方程。其中除SIMPLET外，其它算法只考虑压力修正值对 的影响，而SIMPLET中还考虑了温度修正对 的影响，该影响最后体现在方程的源项之中；pvu,ppvu,uuvvvu,vu,（5）确定速度修正值：lSIMPLE，SIMPLER，SIMPLET：，邻点速度影响略去不计；lSIMPLEC：，邻点速度修正值的影响已基本考虑在内；lSIMPLEX：，de等是通过求解一个Poisson方程得出的，相当于考虑了邻点的影响。（6）采用下列速度、压力值作为初值开始下一层次的迭代计算（设速度的亚松弛已组合在u，v的求解过程中）：lSIMPLE，SIMPLET：lSIM

17、PLEC，SIMPLEX：lSIMLER：以上计算过程可以大致分为两步：预估步（1）（3）及校正步（4）（6）。因而可以说上述五种算法都是两步算法（一步预估，一步校正）。eeeEPeeaAdppdu,nbeeeEPeeaaAdppdu,EPeeppdu papvvuuP,ppvvuu,vvuu,9加速迭代收敛的一些方法加速迭代收敛的一些方法1）选择合适的松驰因子之值）选择合适的松驰因子之值l对于正交的网格以及虽非正交但网格倾斜不是严重的情况，建议au与ap的取值应满足以下关系：（76）式中常数c取为123或1.124，同时au之值应尽可能地取得大，一般可取0.70.8。上述取值原则对交错网格、

18、同位网格、正交网格与非正交网格原则上都适用，但对网格严重倾斜的情况，ap之值应小于按式（76）得出之值。caapu2）其它加速迭代收敛的方法）其它加速迭代收敛的方法迭代循环方式对收敛速率的影响描述燃烧过程的是一组互相耦合的方程组，在用迭代法求解时，迭代循环方式对收敛速率有很大的影响可将流体力学方程（动量和修正压力方程）和其它标量方程（能量、湍流量方程等）一起迭代另一种选择是将流体力学方程首先迭代若干次，然后把其它标量方程（能量、湍流量方程等）一起迭代；流体力学方程首先迭代若干次，然后把其它标量方程再加进去一起循环两者的收敛速率是不一样的。但目前还给不出最佳循环方式的一般原则网格均匀度对收敛速率

19、的影响网格均匀度对收敛速率有很大的影响计算时间与结点数成正比但增加结点数，若能使差分网格更均匀，由于收敛速率加快，有时仍有可能减少总的计算时间七、边界条件的数值处理方法七、边界条件的数值处理方法1、引言、引言 初始条件和边界条件初始条件和边界条件:对抛物型问题，只需给出一端的边界条件，如时间坐标上给出初始的参数分布，边界层流动中给定进口截面上的参数分布等等对椭圆型问题，则需给出求解域两端的边界条件，如果三个空间坐标都是椭圆型坐标，则求解域的四周各点都需给出边界条件2、进口边界、进口边界 模型计算湍流时，一般无法获得由实验测得的进口处k及 k:（1）（2）:(1)(2)进口截面上k与 取值的这种

20、不确定性对计算结果的影响（亦即计算结果对进口条件的敏感性）应当在数值计算中进行考察。k172)21%)(5.15.0(muk 322)(5.1imTuk 18234307.0lLLkc来确定从其中中计算而得从1000100,tmtlukc3、出口边界、出口边界l描写有回流的流动与燃烧的控制方程是关于空间坐标的椭圆型方程，从数学上要求在每个坐标方向上都应当有两个边界条件。出口边界应该置在什么位置及怎样确定这些人为设定的计算边界条件，就成了计算流体力学及燃烧过程模拟中的一个重要研究课题。l出口边界应设置在什么位置上出口边界应设置在什么位置上:Patankar曾认为，出口边界不能设置在有回流的截面上

21、，否则得到的解是没有意义的数值实践证明，对出口边界位置选择的这种苛刻的要求是不必要的文献27提出对这种出口边界条件定性的要求是：（1）应当对出口边界附近，尤其是离开出口边界较远处的流场没有影响；（2）出口边界应当是“透明的”,又称为“无反射”（3）无论出口边界设置在何处，应对其它区域中的数值解没有影 响图11求解流场时出口边界位置的确定l目前，文献中已提出了多种关于设置出口边界条件的方法，归纳起来大约有以下几种情形：1）局部单向化假设）局部单向化假设只有当出口截面无回流时才能近似地成立2）充分发展的假设）充分发展的假设假设n为出口截面的法线方向，这一条件的数学表达式为（84）3）法向速度局部质

22、量守恒、切向速度齐次）法向速度局部质量守恒、切向速度齐次Neumann条件条件以图11（b）中出口边界相邻的控制容积为例，取出后放大而示于图12中。在本例中速度u是与出口截面平行的流速，而v为与出口截面垂直的流速。对所示控制容积作质量守恒运算，得（85a）而按齐次Neumann条件有（85b）0n1,11,1,MiMiMiMiuuxyvv1,MiMiuu图12 局部质量守恒的图示 图13 说明切向速度按局部质量守恒的计算图示4）出口截面上的法向速度与切向速度都按局部质量守恒确定）出口截面上的法向速度与切向速度都按局部质量守恒确定l在这一方法中，对图12所示情形，计算式同式（85a），但ui,M

23、的更新则通过对图13所示半个控制容积的质量守恒来获得。以0.5(vi,M+vi,M-1)作为该半个控制容积南侧流速，则可以写出（86a）将边界上的切向流速取本迭代层次之值而其余速度取上一迭代层次的值，得 （86b）l由图11（b）可见，在计算区域右边界，之值为零，由此可推得 ，这样依次计算可以得出出口边界切向速度的更新值。Miv,04211,1,11,1,MiMiMiMiMiMiMiuuuuyvvvx02144,1,1,11,1,MiMiMiMiMiMiMivvvxuuyuuyMiu,1Miu,5）对流边界条件）对流边界条件l用以下公式来更新出口截面法向速度及其它标量：（87a）其中，c为对流

24、速度或相速度，一般取为通道平均流速。上式离散后，可得（87b）这里，vm为出口截面平均法向流速。相当于假定在出口边界上无扩散作用,为一个纯对流问题。0 xct)(1,)(,)(,)1(,kMikMimkMikMivxt6）内点插值方法）内点插值方法 l在文献36中，提出采用由内点的速度来决定出口流速的方法：（88a）其中，为与出口边界相邻的控制容积的宽度，x为出口边界上游第四个截面位置起算的坐标，当网格均分时（如图14所示），有（88b）2xxCxxBAu161227,3,2,1,jLjLjLjLuuuu图14 插值公式（88a）的图示7）各种高阶导数为零的方法）各种高阶导数为零的方法（89）

25、其中，x为出口方向，上角标j及q表示求导的阶数，要求 ，并建议取j=3，q=2。这一类方法实际上相当于提出了由内点来插值的一个方式。导数的阶数越高，涉及的内点数越多，该阶导数取零对解的影响也越微弱。0qqjjxvxu1,1qj 出口边界上边界条件的设定有多种看来相差很大的方法。在从事工程流动与传热问题的数值计算中实用的处理出口边界条件的方法：出口边界上无回流时，可采用局部单向化或法向一阶导数为零的方法；当出口边界有回流时，可以采用法向流速按质量守恒、切向流速按齐次Neumann条件的方法；无论采用哪一种方法，出口截面上的法向速度分布必须满足总体质量守恒的条件。4固体边界固体边界 1）对于温度一

26、类的变量 用附加源项法来处理 2）对于速度固体壁面是受限流动中最常见的边界，原则上应采用无滑移条件，即流体的速度等于壁面的速度（第一类边界条件），层流时无论是对于动量方程中的流速还是连续性方程中的流速都是这样（后者表现在SIMPLE算法压力修正方程中的剩余质量流量的计算上）当采用标准 模型时，平行于壁面的流速及垂直于壁面的流速还存在不同处理方法。对于平行于壁面的流速（例如u），采用无滑移条件，同时壁面上的动力粘性系数由壁面函数法的公式来计算；但对垂直于该壁面的流速（v）则应区分动量方程还是连续方程。在连续性方程中仍然采用无滑移条件（静止壁面上取v=0），而对v动量方程则采用法向一阶导数为零的条

27、件：k0nv 3）对于）对于k采用低Reynolds数 模型时，可直接令壁面上k=0采用标准 模型时一般取第一个内节点处 ，可通过令该点与壁面之间的k方程的扩散系数为零来实现。4）对于）对于O在低Reynolds数 模型中，可以取壁面上为零，也有取法向一阶导数为零的在采用标准模型及壁面函数法时，一般对第一个内节点赋予规定值，同时令第一个内节点与边界之间的扩散系数为零5）组分的浓度如果固体壁面对化学反应无催化作用，固体壁面上各组分的浓度决定于通过壁面的质量渗透率，以及渗透介质的化学成分。如壁面对组分是不可渗透的，则该组分在固体壁面法线方向的导数等于零。kk0nkk5对称边界对称边界对于一切标量及

28、与对称线平行的速度分量，对称边界上采取法向一阶导数为零的条件对于垂直于对称线的速度分量，则取对称线上该分量为零。6壁面函数壁面函数用上述给出的边界条件，在用数值求解燃烧问题时，对壁面附近的流场要特别小心，才能保证计算有足够的精度。这是因为在壁面附近，输运系数变化很陡在层流流动中，是由于通过壁面经常有很强的热通量，因此温度变化很陡造成的；在湍流流动中，是由于湍流输运系数在接近壁面时，迅速趋于零所致如果象建立内部结点的差分方程那样，用线性插值办法计算非结点位置的输运系数，在这种变化很陡的区域，则需相当密的网格，并仔细安排，才能合理计算输运系数值，但所需的计算机容量和计算时间都大大增加在湍流流动中，由于壁面附近湍流雷诺数比较小，在湍流模型中需考虑湍流雷诺数的影响，求解也比流场内部复杂若我们主要关心的是流场内部的情况，对壁面附近的流场可以采用近似的分析方法，建立起内部结点的参数值和通过壁面的通量之间的定量关系 壁面函数，并通过壁面函数给出内部流场的边界条件当近壁流动为相对简单或可简化的流动时，可采用数学分析比较容易地推导其壁面函数。但是，如果近壁为复杂流动，则需要通过实验或数值积分方法来确定壁面函数。