1、1,數(shù)值計算在熱工中應(yīng)用 B,天津理工大學(xué)自動化學(xué)院熱能系,第六章 代數(shù)方程組的求解方法,2,6.1 代數(shù)方程組求解方法概述,6.2 線性代數(shù)方程組迭代方法的構(gòu)造,6.3 線性代數(shù)方程組迭代收斂條件與加速方法,第6章 代數(shù)方程組的求解方法,6.4 促進守恒條件滿足的塊修正技術(shù),6.5 促進各種諧波分量同步衰減的多重網(wǎng)格技術(shù),3,6.1.1 求解代數(shù)方程組的直接法與迭代法,6.1 代數(shù)方程組求解方法概述,6.1.2 迭代解法的基本思想及關(guān)鍵問題,6.1.3 終止代數(shù)方程組迭代的判據(jù),4,6.1.1 求解代數(shù)方程組的直接法與迭代法 1. 直接解法 (direct method),通過有限次運算可以

2、獲得代數(shù)方程組的精確解的,求解方法；如TDMA，PDMA,沒有舍入誤差就得精確解。,通過一組假定的初場，由代數(shù)方程組本身不斷,加以改進以獲得近似解的求解方法。,2. 迭代法(iterative method),6.1 代數(shù)方程組求解方法概述,5,工程流動與傳熱計算代數(shù)方程求解大多采用迭代法： 1)問題多為非線性的，在獲得收斂解之前，各層次代 數(shù)方程的系數(shù)均是臨時的，不必求出其真解；,2)直接解法計算次數(shù)正比于變量數(shù)目的2.53次方， 變量多時，計算次數(shù)十分可觀；迭代法可適時終止迭 代。,6,6.1.1求解代數(shù)方程的TDMA及ADI 方法 一、 求解一維導(dǎo)熱問題代數(shù)方程的三對角陣算法,1. 一維

3、導(dǎo)熱問題代數(shù)方程通用形式,穩(wěn)態(tài)及非穩(wěn)態(tài)隱式 （f0)都要聯(lián)立求解一 組代數(shù)方程： aPTP = aETE + aW TW + b 其系數(shù)矩陣是一三對角陣 (Tri-diagonal matrix )。,每行三個未知數(shù),7,i,系數(shù)除到等號后的項,2. Thomas算法的一般形式 將上式改寫為： ATi = BiTi +1 + CiTi 1 + Di , i = 1, 2,.M 1 (a) 端點條件：i=1, Ci=0; iM1, Bi0 (1) 消元過程把每行的未知數(shù)由三個減少為二個。 設(shè)消元后方程形式為： Ti 1 = Pi1Ti + Qi 1 (b),8,所謂消元就是要找出系數(shù)Pi,Qi與

4、Ai,Bi,Ci,Di間的關(guān)系。 將(b)乘以Ci,并與(a)式相加：,(a) (b),對照,9,i,上式特點：數(shù)學(xué)上是遞歸的（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 時兩點上 未知量的關(guān)系式，將它與(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算法的實施 第一類邊界條件下，求解區(qū)域為i=2,.M1-1=M2。 將消元公式用于i=1, 注意T1是給定的：,T1 = P1T2 + Q1,P1 = 0;,Q1 = T1,因TM1已知，消元從M2開始： TM 2 = PM 2TM 1 + QM 2 注意：采用附加源項法來處理第二類，第三類邊 界條件時，均將第二類，第三類邊界條件問題視為第 一類邊界條件問題，數(shù)學(xué)上的處理與此相同。,12,二、 求解多維導(dǎo)熱問

6、題代數(shù)方程的ADI方法 1. 求解二維非穩(wěn)態(tài)導(dǎo)熱全隱格式代數(shù)方程的方法,變量一維存儲順序與矩陣系數(shù)的關(guān)系,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方向?qū)嵤╇[式;,2D交替方向隱式,14,設(shè) ui,j,k, vi,j,k為兩個中間子時層上的值；,表示n時層x方向二階導(dǎo)數(shù)的中心差分；,第一個 子時層：

7、 第二個 子時層： 第三個 子時層：,15,用von Neumann分析方法可以證明穩(wěn)定性條件為：,iteration)方法極為相似。,表面上看，相對于一維問題允許時間步長放大了3倍； 實際上并不！ 對二維問題PR方法絕對穩(wěn)定。 3. 這種求解非穩(wěn)態(tài)全隱格式的交替方向隱式(ADI- implicit)與求解多維穩(wěn)態(tài)問題的交替方向迭代(ADI-,16,三、 多維問題離散方程系數(shù)矩陣特點 采用二階截差的格式時，二維與三維流動與傳熱 問題離散方程為：,二維,三維 采用一維存儲方式來表示，對2D 具有L1*M1個未 知數(shù)的問題第k個變量的代數(shù)方程為：,17,采用二階截差的格式時二維問題代數(shù)方程等號前

8、只有五項系數(shù)不為零，系數(shù)矩陣為準五對角陣; 由于在 L1*M1項中只有5項不為零，因而稱為稀疏矩陣；因 為未知數(shù)成千上萬，稱為大型稀疏矩陣（large scale sparse matrix)。,當(dāng)采用右圖所示 方式構(gòu)成未知數(shù) 的一維數(shù)組時,18,19,多維流動傳熱問題控制方程離散形成的代數(shù)方程特點： 1) 常物性導(dǎo)熱問題均分網(wǎng)格：系數(shù)矩陣對稱、正定： 2) 其它情形生成的系數(shù)矩陣一般既非正定又非對稱。,構(gòu)成大型系數(shù)矩陣的代數(shù)方程組一般采用迭代求,解方法。,20,則其真解為：,1. 迭代法基本思想,設(shè)所求解的方程組為：,迭代法時要在一個多維空間R (其空間維數(shù)為變量的個,，要求：,數(shù))中構(gòu)造一

9、個序列,當(dāng) k ,一般地，在第k次迭代中有：,2. 迭代法的關(guān)鍵問題 1) 怎樣構(gòu)造迭代序列？(迭代方式) 2) 所構(gòu)造的迭代方式是否收斂？ 3) 怎樣加快迭代收斂速度？,6.1.2 迭代解法的基本思想及關(guān)鍵問題,21,6.1.3 終止代數(shù)方程組迭代的判據(jù) (1)簡單地規(guī)定迭代的輪數(shù)； (2)規(guī)定p方程余量范數(shù)小 于某值； (3)規(guī)定p方程余量范數(shù)相 對值小于某值； (4)規(guī)定變量所謂相對偏差 小于允許值：,22,6.2.1 點(顯式)迭代法explicit iteration,6.2.2 塊(隱式)迭代法implicit iteration,6.2 線性代數(shù)方程組迭代方法的構(gòu)造,6.2.3

10、交替方向塊迭代法ADI,23,6.2 線性代數(shù)方程組迭代方法的構(gòu)造,6.2.1 點迭代法,每一步計算只能改進求解區(qū)域中一點之值；將求解 區(qū)域中各點之值都更新一次的計算過程稱為一輪迭代。 每一點上被更新之值均與其它各點之值顯式相關(guān)。,未知值的更新均用上一輪計算的鄰點之值收斂,1. Jakob迭代,速度與迭代方向無關(guān)。 2. GaussSeidel迭代,未知值的更新均用鄰點的最新值來進行。,24,3. SOR/SUR迭代 1 超松弛 注意：這是代數(shù)方程求解的松弛，不是非線性迭代的 松弛。 6.2.2 塊(隱式)迭代法implicit iteration 1. 基本思想 將求解區(qū)域分為若干塊，每一塊

11、內(nèi)各點采用直接 解法，塊與塊之間的推進采用迭代方式，又稱隱式迭代,法。,25,Jakob:,G-S:,2. 塊迭代法的最基本形式線迭代（Line iteration) 塊的最小單位為線，同一條線上各點用TDMA方 法直接求解，線間推進用迭代方式。 東西掃描南北求解時：,26,6.2.3 交替方向塊迭代法ADI 1. 基本思想,先按行（或列）直接求解，再按列（或行）直接,求解；兩次全場更新，組成一輪迭代。,Alternative direction iteration(ADI)與交替方向 隱式(ADI)之間的聯(lián)系：,27,（1） 交替方向Jakob式迭代可以表示為：,（2）,這與非穩(wěn)態(tài)導(dǎo)熱問題的

12、Peaceman-Rachford的交替 方向隱式十分相似：,28,2-D Peaceman-Rachford方法 將 t 兩等分： 第一個 t / 2 X方向為隱式， y方向為顯式； 第二個 t / 2 在y 方向?qū)嵤╇[式、X方向為顯式 2D交替方向隱式,表示k時層x方向二階導(dǎo)數(shù)的中心差分；,設(shè),為中間子時層上的值；,29,(1),第一個 子時層：,第二個 子時層：,將第一式展開：,非穩(wěn)態(tài)步進一個時層相當(dāng)于穩(wěn)態(tài)完成一輪迭代。,(2),30,凡結(jié)構(gòu)化網(wǎng)格上的代數(shù)方程求解均可采用：,2. 交替方向線迭代廣泛應(yīng)用于傳熱與流動問題的求解,這種求解求解多維穩(wěn)態(tài)問題的交替方向迭代,(ADI-iterat

13、ion)方法與非穩(wěn)態(tài)全隱格式的交替方向 隱式(ADI-implicit)與極為相似。,31,6.3.1 Jakob,G-S 迭代收斂充分條件及其分析,6.3.2 影響迭代收斂速度的因素分析,6.3 線性代數(shù)方程組迭代收斂條件與加速方法,6.2.3 加速邊界條件影響傳入計算區(qū)域的方法,32,predominant :,6.3 線性代數(shù)方程組迭代收斂條件與加速方法 6.3.1 Jakob,G-S 迭代收斂充分條件及其分析 1. 充分條件Scarborough準則 系數(shù)矩陣不可約且按行或按列對角占優(yōu)-diagonal,對于各行成立，且至少對一行不等號成立。 2. 對流擴散問題離散方程系數(shù)特性分析,3

14、3,1)系數(shù)矩陣不可約所謂可約是指由矩陣下標元素組 成的集合W一定可以分為兩個非空子集R與S，它們 滿足WRS，而且從R和S中各取任一元素k及 l，,一定有：ak ,l, 0 ；如果W不滿足這樣的條件則稱不,可約non-reducible. 分析：離散方程的系數(shù)代表了節(jié)點間的相互影響；對于 空間上為橢圓型的問題，區(qū)域中的任一點，必對其周 圍的點有影響，而矩陣可約就相當(dāng)于一定可以將計算 區(qū)域分割為兩個互不影響的兩塊，這是不可能的。,34,矩陣不可約是流動 與熱量傳遞過程總是互 相影響的這一物理事實 的數(shù)學(xué)反映。 2) 矩陣對角占優(yōu)特性分析本課程推薦的方法一定 滿足對角占優(yōu)，而且至少一行不等號成立

15、！,(1) 非穩(wěn)態(tài)全隱格式：,(2) 穩(wěn)態(tài)有源項：,35,(3) 穩(wěn)態(tài)無源項：對內(nèi)點,從邊界點必可找出至少一點滿足不等號： 設(shè)Tw已知，則,TW,求解時變?yōu)椋?故此點有： 對第三類邊界條件的邊界點，則有附件源項,36,如果邊界點均為第二類邊界條件，不可能！至少有一 點為第一類或者第三類；否則無確定解！,本課程推薦的方法一定滿足對角占優(yōu)，而且至少一,行不等號成立。,6.3.2 影響迭代收斂速度的因素分析 1. 從邊界條件影響傳入內(nèi)部的角度：,常物性無內(nèi)熱源物體穩(wěn)態(tài)導(dǎo)熱，控制方程為Laplace 方程 ；從均勻初場開始迭代；均勻場滿足,Laplace方程，但不是問題的解，原因在于不滿足邊,37,界

16、條件；因此邊界條件影響傳入的快慢必影響迭代收斂 速度。,2. 從滿足守恒條件的角度：,對于第一類邊界條件的問題，可以將全部邊界條,件組織到迭代初場中，但是這樣的場不滿足守恒定律， 因此能加快守恒條件滿足的方法必可加快收斂；,3. 從誤差矢量衰減的角度：,迭代過程是使誤差矢量逐漸衰減的過程，誤差矢量 由各諧波分量組成，能使各種諧波分量協(xié)調(diào)地衰減的方 法能促使迭代過程地收斂；,38,4. 從增加直接求解分量的角度：,直接求解是能使各點同時滿足守恒條件與邊界條件 的最強方式，適當(dāng)增加直接求解的分量有助于提高迭代 收斂速度。,6.2.3 加速邊界條件影響傳入計算區(qū)域的方法,Jakob迭代：每一輪迭代

17、邊界條件的直接影響只能 傳入一個網(wǎng)格，收斂速度 最慢。,39,GS迭代：每一輪迭代 起始邊的影響傳到整個計 算區(qū)域內(nèi)，收斂速度加快。,線迭代：每一輪迭代起始 邊、上下邊界條件的影響 全部傳入計算區(qū)域，收斂 進一步加快。,40,ADI線迭代：每一輪迭全 部邊界條件的影響均傳入 計算區(qū)域，收斂進一步最 快。,ADI線迭代線迭代G-S迭代Jakob迭代,Jakob迭代收斂速度最慢，相鄰兩次迭代間的變化 最小，有時被用來克服外迭代不收斂的困難；SIMPLEST 算法中對流項采用Jakob迭代，在規(guī)定迭代次數(shù)的情況 有利于非線性強烈的問題外迭代收斂。,41,6.4.1 塊修正技術(shù)的需要,6.4.2 塊修

18、正技術(shù)的基本思想,6.4 促進守恒條件滿足的塊修正技術(shù),6.4.3 單塊修正值的代數(shù)方程及其邊界條件,6.4.4 應(yīng)用塊修正技術(shù)注意事項,42,6.4 促進守恒條件滿足的塊修正技術(shù),6.4.1 塊修正技術(shù)的需要,對下圖示二維穩(wěn)態(tài)無內(nèi)熱源的導(dǎo)熱問題，采用常 規(guī)的ADI方法求解時，迭代收斂速度大大減慢：東西 方向影響最強，但離散方程系數(shù)很??；南北絕熱邊界 沒有提供什么信息，但系數(shù)很大因此需要一種能加 速收斂的方法。,43,而是一個塊加一個共同的修正值， 記為 或者,要求 或者 能滿足該塊的守恒要求。,6.4.2 塊修正技術(shù)的基本思想,從物理上看，所謂迭代收斂就是計算結(jié)果逐漸滿 足守恒條件的過程；所

19、謂迭代一輪，就是在原來不滿,是代數(shù)方程求解過程。,足守恒特性的解 上加上一個修正值 使得 更能滿足守恒要求。獲得 的過程就,對二維問題，修正值也是二維的；為了只求解一 個一維的未知數(shù)，不在每個節(jié)點上加單獨的修正值，,44,1.修正值方程： 要求,滿足該塊的守恒條件：,6.4.3 單塊修正值的代數(shù)方程及其邊界條件,IST為X方向開始迭代的下標；L2為右端第一個內(nèi)接點,45,46,我們采用附加源項法來處理第二類，第 三類邊界條件，相當(dāng)于所有的邊界條件 均是第一類的，邊界上的修正值永遠為 零；所以對于修正值求和時JJST以及 J=M2均為零，導(dǎo)致AJM需不計JST， 而AJP需不計M2。 但該兩個位

20、置上 值均存在，因此在常數(shù)項BLC中均 應(yīng)計及。,47,2.修正值邊界條件,48,6.4.4 應(yīng)用塊修正技術(shù)注意事項,1.塊修正不是一個獨立的求解方法，須與其它方法，如 ADI聯(lián)合使用；,2.為進一步加速收斂可以交替方向使用塊修正； 3.對于物理意義上恒大于零的量， 如湍流脈動動能、組分等有時不宜 采用塊修正技術(shù)。因為塊修正相當(dāng) 于對某一個坐標系均勻地加上或者 減去一個常數(shù)，可能導(dǎo)致不符合物 理意義的解。,49,6.5.1 代數(shù)方程迭代求解過程是誤差矢量的衰減過程,6.5.2 多重網(wǎng)格的基本思想及關(guān)鍵問題,6.5 促進各種諧波分量同步衰減的多重網(wǎng)格技術(shù),6.5.3 不同網(wǎng)格上解的傳遞,6.5.

21、4 不同網(wǎng)格上的循環(huán)方式,50,6.5 促進各種諧波分量同步衰減的多重網(wǎng)格技術(shù) 6.5.1 代數(shù)方程迭代求解過程是誤差矢量的衰減過程 1.迭代過程中誤差矢量的衰減規(guī)律 采用von Neumann分析法研究迭代過程中誤差 衰減情況。設(shè)有一維穩(wěn)態(tài)有源項導(dǎo)熱問題：,在均分網(wǎng)格上離散后有： 采用自左至右的G-S迭代，則有：,51,設(shè)在第k輪迭代中引入的誤差矢量為 ，分量為 則 ，根據(jù)第三章的推導(dǎo)知道誤差矢 量傳遞規(guī)律為：,2. 諧波分量衰減過程分析,將,表示成為：,代入上式，得,增長因子,52,據(jù)不同輻角來分析增長因子的數(shù)值：,53,按 von Neumann分析：,在一定網(wǎng)格步長下，短波分量輻角大，

22、很快就衰減， 長波分量輻角小，難以衰減，使迭代收斂速度減慢。 但短波或者高頻是相對于一定步長而言的；如果迭 代幾次增加網(wǎng)格步長，則長波就會稱為短波，很快 衰減，如此幾次增加網(wǎng)格步長，可望將不同頻率的 諧波分量均勻地衰減掉，促進迭代收斂。 一般將 /2 范圍內(nèi)的波長稱為短波分量。,54,6.5.2 多重網(wǎng)格的基本思想及關(guān)鍵問題,1. 基本思想在幾套不同疏密的網(wǎng)格上對密網(wǎng)格上 的解進行求解，使不同頻率的誤差分量得以衰減。 1. 關(guān)鍵問題,(1) 相鄰兩套網(wǎng)格間密網(wǎng)格的解怎樣傳遞？ (2) 不同疏密網(wǎng)格間密網(wǎng)格的解如何循環(huán)？,6.5.3 不同網(wǎng)格上密網(wǎng)格的解解的傳遞,1. 基本觀點不同網(wǎng)格上傳遞的是最密網(wǎng)格上的解,55,所謂多重網(wǎng)格方法就是為了克服固定網(wǎng)格的缺點而發(fā)展起來的迭代解法。采用此法時可先在較細的網(wǎng)格上進行迭