数值计算在热工中应用B第六章-1.ppt

1、1,数值计算在热工中应用 B,天津理工大学自动化学院热能系,第六章 代数方程组的求解方法,2,6.1 代数方程组求解方法概述,6.2 线性代数方程组迭代方法的构造,6.3 线性代数方程组迭代收敛条件与加速方法,第6章 代数方程组的求解方法,6.4 促进守恒条件满足的块修正技术,6.5 促进各种谐波分量同步衰减的多重网格技术,3,6.1.1 求解代数方程组的直接法与迭代法,6.1 代数方程组求解方法概述,6.1.2 迭代解法的基本思想及关键问题,6.1.3 终止代数方程组迭代的判据,4,6.1.1 求解代数方程组的直接法与迭代法 1. 直接解法 (direct method),通过有限次运算可以

2、获得代数方程组的精确解的,求解方法；如TDMA，PDMA,没有舍入误差就得精确解。,通过一组假定的初场，由代数方程组本身不断,加以改进以获得近似解的求解方法。,2. 迭代法(iterative method),6.1 代数方程组求解方法概述,5,工程流动与传热计算代数方程求解大多采用迭代法： 1)问题多为非线性的，在获得收敛解之前，各层次代 数方程的系数均是临时的，不必求出其真解；,2)直接解法计算次数正比于变量数目的2.53次方， 变量多时，计算次数十分可观；迭代法可适时终止迭 代。,6,6.1.1求解代数方程的TDMA及ADI 方法 一、 求解一维导热问题代数方程的三对角阵算法,1. 一维

3、导热问题代数方程通用形式,稳态及非稳态隐式 （f0)都要联立求解一 组代数方程： aPTP = aETE + aW TW + b 其系数矩阵是一三对角阵 (Tri-diagonal matrix )。,每行三个未知数,7,i,系数除到等号后的项,2. Thomas算法的一般形式 将上式改写为： ATi = BiTi +1 + CiTi 1 + Di , i = 1, 2,.M 1 (a) 端点条件：i=1, Ci=0; iM1, Bi0 (1) 消元过程把每行的未知数由三个减少为二个。 设消元后方程形式为： Ti 1 = Pi1Ti + Qi 1 (b),8,所谓消元就是要找出系数Pi,Qi与

4、Ai,Bi,Ci,Di间的关系。 将(b)乘以Ci,并与(a)式相加：,(a) (b),对照,9,i,上式特点：数学上是递归的（recurrent)首先 必须知道P1，Q1。 为此，试重新审视式(a) ATi = BiTi +1 + CiTi 1 + Di , i = 1, 2,.M 1 (a) 端点条件：i=1, Ci=0; iM1, Bi0 如果将(a)式用于i=1,则立即可得出i=1 时两点上 未知量的关系式，将它与(b) 相比就能得出P1，Q1。,10,i = 1, C1 = 0,A1T1 = B1T2 + D1,(2) 回代过程从M1点开始，利用式(b) 逐一得出Ti。,TM 1 =

5、 PM 1TM 1+1 + QM 1 ,端点条件：iM1, Bi0,TM 1 = QM 1,逐一得出：TM11，.T2,T1 。,PM 1 = 0,11,3. 第一类边界条件下Thomas算法的实施 第一类边界条件下，求解区域为i=2,.M1-1=M2。 将消元公式用于i=1, 注意T1是给定的：,T1 = P1T2 + Q1,P1 = 0;,Q1 = T1,因TM1已知，消元从M2开始： TM 2 = PM 2TM 1 + QM 2 注意：采用附加源项法来处理第二类，第三类边 界条件时，均将第二类，第三类边界条件问题视为第 一类边界条件问题，数学上的处理与此相同。,12,二、 求解多维导热问

6、题代数方程的ADI方法 1. 求解二维非稳态导热全隐格式代数方程的方法,变量一维存储顺序与矩阵系数的关系,13,(1) 五对角阵算法(Penta-diagonal ,PDMA) (2) 交替方向隐式方法 ( Alternative-direction Implicit, ADI) 2. 3-D Peaceman-Rachford方法 将 t 三等分： 第一个 t / 3 X方向为隐式， y,z方向为显式； 第二，三个 t / 3 分别在,y,z方向实施隐式;,2D交替方向隐式,14,设 ui,j,k, vi,j,k为两个中间子时层上的值；,表示n时层x方向二阶导数的中心差分；,第一个 子时层：

7、 第二个 子时层： 第三个 子时层：,15,用von Neumann分析方法可以证明稳定性条件为：,iteration)方法极为相似。,表面上看，相对于一维问题允许时间步长放大了3倍； 实际上并不！ 对二维问题PR方法绝对稳定。 3. 这种求解非稳态全隐格式的交替方向隐式(ADI- implicit)与求解多维稳态问题的交替方向迭代(ADI-,16,三、 多维问题离散方程系数矩阵特点 采用二阶截差的格式时，二维与三维流动与传热 问题离散方程为：,二维,三维 采用一维存储方式来表示，对2D 具有L1*M1个未 知数的问题第k个变量的代数方程为：,17,采用二阶截差的格式时二维问题代数方程等号前

8、只有五项系数不为零，系数矩阵为准五对角阵; 由于在 L1*M1项中只有5项不为零，因而称为稀疏矩阵；因 为未知数成千上万，称为大型稀疏矩阵（large scale sparse matrix)。,当采用右图所示 方式构成未知数 的一维数组时,18,19,多维流动传热问题控制方程离散形成的代数方程特点： 1) 常物性导热问题均分网格：系数矩阵对称、正定： 2) 其它情形生成的系数矩阵一般既非正定又非对称。,构成大型系数矩阵的代数方程组一般采用迭代求,解方法。,20,则其真解为：,1. 迭代法基本思想,设所求解的方程组为：,迭代法时要在一个多维空间R (其空间维数为变量的个,，要求：,数)中构造一

9、个序列,当 k ,一般地，在第k次迭代中有：,2. 迭代法的关键问题 1) 怎样构造迭代序列？(迭代方式) 2) 所构造的迭代方式是否收敛？ 3) 怎样加快迭代收敛速度？,6.1.2 迭代解法的基本思想及关键问题,21,6.1.3 终止代数方程组迭代的判据 (1)简单地规定迭代的轮数； (2)规定p方程余量范数小 于某值； (3)规定p方程余量范数相 对值小于某值； (4)规定变量所谓相对偏差 小于允许值：,22,6.2.1 点(显式)迭代法explicit iteration,6.2.2 块(隐式)迭代法implicit iteration,6.2 线性代数方程组迭代方法的构造,6.2.3

10、交替方向块迭代法ADI,23,6.2 线性代数方程组迭代方法的构造,6.2.1 点迭代法,每一步计算只能改进求解区域中一点之值；将求解 区域中各点之值都更新一次的计算过程称为一轮迭代。 每一点上被更新之值均与其它各点之值显式相关。,未知值的更新均用上一轮计算的邻点之值收敛,1. Jakob迭代,速度与迭代方向无关。 2. GaussSeidel迭代,未知值的更新均用邻点的最新值来进行。,24,3. SOR/SUR迭代 1 超松弛 注意：这是代数方程求解的松弛，不是非线性迭代的 松弛。 6.2.2 块(隐式)迭代法implicit iteration 1. 基本思想 将求解区域分为若干块，每一块

11、内各点采用直接 解法，块与块之间的推进采用迭代方式，又称隐式迭代,法。,25,Jakob:,G-S:,2. 块迭代法的最基本形式线迭代（Line iteration) 块的最小单位为线，同一条线上各点用TDMA方 法直接求解，线间推进用迭代方式。 东西扫描南北求解时：,26,6.2.3 交替方向块迭代法ADI 1. 基本思想,先按行（或列）直接求解，再按列（或行）直接,求解；两次全场更新，组成一轮迭代。,Alternative direction iteration(ADI)与交替方向 隐式(ADI)之间的联系：,27,（1） 交替方向Jakob式迭代可以表示为：,（2）,这与非稳态导热问题的

12、Peaceman-Rachford的交替 方向隐式十分相似：,28,2-D Peaceman-Rachford方法 将 t 两等分： 第一个 t / 2 X方向为隐式， y方向为显式； 第二个 t / 2 在y 方向实施隐式、X方向为显式 2D交替方向隐式,表示k时层x方向二阶导数的中心差分；,设,为中间子时层上的值；,29,(1),第一个 子时层：,第二个 子时层：,将第一式展开：,非稳态步进一个时层相当于稳态完成一轮迭代。,(2),30,凡结构化网格上的代数方程求解均可采用：,2. 交替方向线迭代广泛应用于传热与流动问题的求解,这种求解求解多维稳态问题的交替方向迭代,(ADI-iterat

13、ion)方法与非稳态全隐格式的交替方向 隐式(ADI-implicit)与极为相似。,31,6.3.1 Jakob,G-S 迭代收敛充分条件及其分析,6.3.2 影响迭代收敛速度的因素分析,6.3 线性代数方程组迭代收敛条件与加速方法,6.2.3 加速边界条件影响传入计算区域的方法,32,predominant :,6.3 线性代数方程组迭代收敛条件与加速方法 6.3.1 Jakob,G-S 迭代收敛充分条件及其分析 1. 充分条件Scarborough准则 系数矩阵不可约且按行或按列对角占优-diagonal,对于各行成立，且至少对一行不等号成立。 2. 对流扩散问题离散方程系数特性分析,3

14、3,1)系数矩阵不可约所谓可约是指由矩阵下标元素组 成的集合W一定可以分为两个非空子集R与S，它们 满足WRS，而且从R和S中各取任一元素k及 l，,一定有：ak ,l, 0 ；如果W不满足这样的条件则称不,可约non-reducible. 分析：离散方程的系数代表了节点间的相互影响；对于 空间上为椭圆型的问题，区域中的任一点，必对其周 围的点有影响，而矩阵可约就相当于一定可以将计算 区域分割为两个互不影响的两块，这是不可能的。,34,矩阵不可约是流动 与热量传递过程总是互 相影响的这一物理事实 的数学反映。 2) 矩阵对角占优特性分析本课程推荐的方法一定 满足对角占优，而且至少一行不等号成立

15、！,(1) 非稳态全隐格式：,(2) 稳态有源项：,35,(3) 稳态无源项：对内点,从边界点必可找出至少一点满足不等号： 设Tw已知，则,TW,求解时变为： 故此点有： 对第三类边界条件的边界点，则有附件源项,36,如果边界点均为第二类边界条件，不可能！至少有一 点为第一类或者第三类；否则无确定解！,本课程推荐的方法一定满足对角占优，而且至少一,行不等号成立。,6.3.2 影响迭代收敛速度的因素分析 1. 从边界条件影响传入内部的角度：,常物性无内热源物体稳态导热，控制方程为Laplace 方程 ；从均匀初场开始迭代；均匀场满足,Laplace方程，但不是问题的解，原因在于不满足边,37,界

16、条件；因此边界条件影响传入的快慢必影响迭代收敛 速度。,2. 从满足守恒条件的角度：,对于第一类边界条件的问题，可以将全部边界条,件组织到迭代初场中，但是这样的场不满足守恒定律， 因此能加快守恒条件满足的方法必可加快收敛；,3. 从误差矢量衰减的角度：,迭代过程是使误差矢量逐渐衰减的过程，误差矢量 由各谐波分量组成，能使各种谐波分量协调地衰减的方 法能促使迭代过程地收敛；,38,4. 从增加直接求解分量的角度：,直接求解是能使各点同时满足守恒条件与边界条件 的最强方式，适当增加直接求解的分量有助于提高迭代 收敛速度。,6.2.3 加速边界条件影响传入计算区域的方法,Jakob迭代：每一轮迭代

17、边界条件的直接影响只能 传入一个网格，收敛速度 最慢。,39,GS迭代：每一轮迭代 起始边的影响传到整个计 算区域内，收敛速度加快。,线迭代：每一轮迭代起始 边、上下边界条件的影响 全部传入计算区域，收敛 进一步加快。,40,ADI线迭代：每一轮迭全 部边界条件的影响均传入 计算区域，收敛进一步最 快。,ADI线迭代线迭代G-S迭代Jakob迭代,Jakob迭代收敛速度最慢，相邻两次迭代间的变化 最小，有时被用来克服外迭代不收敛的困难；SIMPLEST 算法中对流项采用Jakob迭代，在规定迭代次数的情况 有利于非线性强烈的问题外迭代收敛。,41,6.4.1 块修正技术的需要,6.4.2 块修

18、正技术的基本思想,6.4 促进守恒条件满足的块修正技术,6.4.3 单块修正值的代数方程及其边界条件,6.4.4 应用块修正技术注意事项,42,6.4 促进守恒条件满足的块修正技术,6.4.1 块修正技术的需要,对下图示二维稳态无内热源的导热问题，采用常 规的ADI方法求解时，迭代收敛速度大大减慢：东西 方向影响最强，但离散方程系数很小；南北绝热边界 没有提供什么信息，但系数很大因此需要一种能加 速收敛的方法。,43,而是一个块加一个共同的修正值， 记为 或者,要求 或者 能满足该块的守恒要求。,6.4.2 块修正技术的基本思想,从物理上看，所谓迭代收敛就是计算结果逐渐满 足守恒条件的过程；所

19、谓迭代一轮，就是在原来不满,是代数方程求解过程。,足守恒特性的解 上加上一个修正值 使得 更能满足守恒要求。获得 的过程就,对二维问题，修正值也是二维的；为了只求解一 个一维的未知数，不在每个节点上加单独的修正值，,44,1.修正值方程： 要求,满足该块的守恒条件：,6.4.3 单块修正值的代数方程及其边界条件,IST为X方向开始迭代的下标；L2为右端第一个内接点,45,46,我们采用附加源项法来处理第二类，第 三类边界条件，相当于所有的边界条件 均是第一类的，边界上的修正值永远为 零；所以对于修正值求和时JJST以及 J=M2均为零，导致AJM需不计JST， 而AJP需不计M2。 但该两个位

20、置上 值均存在，因此在常数项BLC中均 应计及。,47,2.修正值边界条件,48,6.4.4 应用块修正技术注意事项,1.块修正不是一个独立的求解方法，须与其它方法，如 ADI联合使用；,2.为进一步加速收敛可以交替方向使用块修正； 3.对于物理意义上恒大于零的量， 如湍流脉动动能、组分等有时不宜 采用块修正技术。因为块修正相当 于对某一个坐标系均匀地加上或者 减去一个常数，可能导致不符合物 理意义的解。,49,6.5.1 代数方程迭代求解过程是误差矢量的衰减过程,6.5.2 多重网格的基本思想及关键问题,6.5 促进各种谐波分量同步衰减的多重网格技术,6.5.3 不同网格上解的传递,6.5.

21、4 不同网格上的循环方式,50,6.5 促进各种谐波分量同步衰减的多重网格技术 6.5.1 代数方程迭代求解过程是误差矢量的衰减过程 1.迭代过程中误差矢量的衰减规律 采用von Neumann分析法研究迭代过程中误差 衰减情况。设有一维稳态有源项导热问题：,在均分网格上离散后有： 采用自左至右的G-S迭代，则有：,51,设在第k轮迭代中引入的误差矢量为 ，分量为 则 ，根据第三章的推导知道误差矢 量传递规律为：,2. 谐波分量衰减过程分析,将,表示成为：,代入上式，得,增长因子,52,据不同辐角来分析增长因子的数值：,53,按 von Neumann分析：,在一定网格步长下，短波分量辐角大，

22、很快就衰减， 长波分量辐角小，难以衰减，使迭代收敛速度减慢。 但短波或者高频是相对于一定步长而言的；如果迭 代几次增加网格步长，则长波就会称为短波，很快 衰减，如此几次增加网格步长，可望将不同频率的 谐波分量均匀地衰减掉，促进迭代收敛。 一般将 /2 范围内的波长称为短波分量。,54,6.5.2 多重网格的基本思想及关键问题,1. 基本思想在几套不同疏密的网格上对密网格上 的解进行求解，使不同频率的误差分量得以衰减。 1. 关键问题,(1) 相邻两套网格间密网格的解怎样传递？ (2) 不同疏密网格间密网格的解如何循环？,6.5.3 不同网格上密网格的解解的传递,1. 基本观点不同网格上传递的是最密网格上的解,55,所谓多重网格方法就是为了克服固定网格的缺点而发展起来的迭代解法。采用此法时可先在较细的网格上进行迭