計(jì)算傳熱學(xué)第講擴(kuò)散方程的數(shù)值解

計(jì)算傳熱學(xué)第講擴(kuò)散方程的數(shù)值解第1頁(yè)，共83頁(yè)，2023年，2月20日，星期四主要內(nèi)容一維穩(wěn)態(tài)問題的數(shù)值解一維非穩(wěn)態(tài)問題多維非穩(wěn)態(tài)問題的離散化差分方程的求解第2頁(yè)，共83頁(yè)，2023年，2月20日，星期四主要目的掌握用數(shù)值方法求解傳熱問題的整體步驟數(shù)值方法的計(jì)算機(jī)實(shí)現(xiàn)邊界條件的處理第3頁(yè)，共83頁(yè)，2023年，2月20日，星期四閱讀與作業(yè)閱讀要求：陶文銓《數(shù)值傳熱學(xué)》第4章作業(yè)：P124題4－1；P125題4－7完成課外作業(yè)第一題和第二題第4頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解控制方程：其中，A(x)是面積函數(shù)。定義如下：直角坐標(biāo)系：A(x)＝1（無限大平板導(dǎo)熱問題）柱坐標(biāo)系：A(x)＝x（極坐標(biāo)系中的一維問題，無限長(zhǎng)圓筒壁導(dǎo)熱問題）球坐標(biāo)系：A(x)＝x2（通過球壁的導(dǎo)熱）變截面肋片：A(x)第5頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解4.1.1求解區(qū)域的離散化x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE圖1一維問題空間區(qū)域的離散化第6頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解4.1.2源項(xiàng)的線性化在通常情況下，S＝S（T）線性化：

S＝Sc＋SpT

（2）其中，按負(fù)斜率源項(xiàng)原則，

Sp＝Sp(T*)0（3）第7頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解4.1.3控制方程的離散化將方程（1）兩邊通乘A（x），并對(duì)x從w到e積分：x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE圖1一維問題空間區(qū)域的離散化第8頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE圖1一維問題空間區(qū)域的離散化第9頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解在上面的積分過程中，我們假定：待求變量T在控制容積P上為常數(shù)整個(gè)控制容積的A(x)為常數(shù)，且等于P點(diǎn)的值。第10頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解將（5）和（6）代入方程（4），整理后得到，其中，第11頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解其中，下標(biāo)：大寫字母表示在節(jié)點(diǎn)處取值，小寫字母表示在相應(yīng)的控制面處取值第12頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解可能的改進(jìn)方案：對(duì)源項(xiàng)積分時(shí)采用線性分布第13頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解4.1.4交界面參數(shù)的計(jì)算線性插值法（算術(shù)平均）調(diào)和平均法待求變量插值Kirchhoff變換法第14頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解4.1.5躍變界面的處理把躍變界面作為控制面調(diào)和平均法把躍變界面作為節(jié)點(diǎn)算術(shù)平均法Kirchhoff變換法待求變量插值法把躍變界面放置在其它位置所有方法都適用把躍變界面作為邊界第15頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解把躍變界面作為邊界可以考慮接觸熱阻rc

(W·m2)/K滿足流的唯一性原則，第16頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解4.1.6邊界條件的處理直角坐標(biāo)左邊界，第二類邊界條件qBx=0注意：qB的正方向與x軸的正方向一致??！第17頁(yè)，共83頁(yè)，2023年，2月20日，星期四ee邊界條件的處理網(wǎng)格是用外節(jié)點(diǎn)法劃分的邊界上出現(xiàn)半個(gè)控制容積qBx=0123(x)1(x)2邊界節(jié)點(diǎn)的差分方程可以用下述方法推出：一階精度的Taylor級(jí)數(shù)展開法第18頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理整理后得到：特點(diǎn)：最簡(jiǎn)單的處理方法只有一階精度與控制方程的精度不匹配第19頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理元體能量平衡法：在研究邊界節(jié)點(diǎn)所代表的控制容積（元體）的能量平衡流入CV的能量＋內(nèi)熱源發(fā)出的熱量＝流出CV的能量123eeqBx=0(x)1(x)21?(x)1Sq1流入CV的能量通過邊界流入的熱量qB通過控制面流入的熱量q1內(nèi)熱源發(fā)出的熱量＝?(x)1S流出CV的能量＝0第20頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理代入能量守恒關(guān)系，整理后得到，特點(diǎn)靈活，便于處理各種復(fù)雜的邊界條件二階精度，與內(nèi)部節(jié)點(diǎn)精度等級(jí)匹配（請(qǐng)證明?。。┩茖?dǎo)過程較繁第21頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理控制方程法在直角坐標(biāo)的條件下，方程（1）變?yōu)椋?23eeqBx=0(x)1(x)21假定在節(jié)點(diǎn)1～2之間的導(dǎo)熱系數(shù)為常數(shù)，且恒等于e，則有，對(duì)于點(diǎn)e，第22頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理另一方面，參照附圖，123eeqBx=0(x)1(x)21第23頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理所以，將之代入式（8）第24頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理整理后得到，特點(diǎn)二階精度不具有一般性推導(dǎo)繁瑣第25頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理二階精度的Taylor級(jí)數(shù)展開法123eeqBx=0(x)1(x)21按二階精度的差商公式第26頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理代入式（16），整理后得到，代入方程（8），整理后得到，第27頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理二階精度具有一般性增加計(jì)算工作量一般很少采用第28頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理求解區(qū)域是用內(nèi)節(jié)點(diǎn)法離散化的qBx=0(x)2(x)3123(x)1(x)2邊界節(jié)點(diǎn)的控制容積或它所代表的求解區(qū)域？邊界節(jié)點(diǎn)的控制容積＝0于是，按元體能量平衡法或其它二階精度的方法，令與源項(xiàng)對(duì)應(yīng)的項(xiàng)等于0，得到，第29頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理說明：盡管它與一階Taylor級(jí)數(shù)展開法的結(jié)果形式上相同，但它卻是二階精度的！請(qǐng)大家證明這一結(jié)論。采用內(nèi)節(jié)點(diǎn)法劃分網(wǎng)格時(shí)，即使在均勻網(wǎng)絡(luò)的前提下，第1個(gè)近邊界節(jié)點(diǎn)也不是等步長(zhǎng)的。qBx=0(x)2(x)3123(x)1(x)2從圖中可以清楚地看出這一點(diǎn)即使(x)2=(x)3

(x)1也不等于

(x)2所以要對(duì)第一個(gè)內(nèi)部節(jié)點(diǎn)給予特別注意。第30頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理例如，對(duì)于直角坐標(biāo)系，對(duì)于節(jié)點(diǎn)2，qBx=0(x)2(x)3123(x)1(x)2注意：CV2的左控制面w與節(jié)點(diǎn)W（節(jié)點(diǎn)1）重合，即與左邊界重合！控制面e！第31頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理或?qū)懗?，?2頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理附加源項(xiàng)法（Additionalsourcetermmethod)以內(nèi)節(jié)點(diǎn)法為例由方程（19）解出邊界節(jié)點(diǎn)上的待求變量T1，代入與第1個(gè)近邊界節(jié)點(diǎn)的差分方程（21），第33頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理代入與第1個(gè)近邊界節(jié)點(diǎn)的差分方程（21），整理后得到，第34頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理整理后得到，或者寫成，其中，第35頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理其中Additionalsourceterm!第36頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理對(duì)于第3類邊界條件，也可以做類似的處理，但是這時(shí)，請(qǐng)大家證明，第37頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理附加源項(xiàng)法的實(shí)質(zhì)邊界節(jié)點(diǎn)消去法不僅能用于內(nèi)節(jié)點(diǎn)網(wǎng)格，也能用于外節(jié)點(diǎn)網(wǎng)格實(shí)施方法：計(jì)算附加源項(xiàng)：Sc,ad，Sp,ad把附加源項(xiàng)計(jì)入該控制容積中的源項(xiàng)中令與邊界節(jié)點(diǎn)對(duì)應(yīng)的系數(shù)（aW）等于0第38頁(yè)，共83頁(yè)，2023年，2月20日，星期四特別提示邊界條件的處理是傳熱問題數(shù)值計(jì)算最重要的環(huán)節(jié)之一元體能量平衡法的基礎(chǔ)地位盡可能采用外節(jié)點(diǎn)法劃分網(wǎng)格邊界節(jié)點(diǎn)消去法廣泛應(yīng)用提高收斂速度盡可能采用均勻網(wǎng)格邊界節(jié)點(diǎn)的離散化方程在形式上與內(nèi)部節(jié)點(diǎn)的相同第39頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1.7差分方程的求解將前面得到的差分方程（8）改寫為，或者簡(jiǎn)單的寫成矩陣的形式，其中，第40頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1.7差分方程的求解第41頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1.7差分方程的求解與方程（33）對(duì)比，知，由方程（33）系數(shù)陣[A]的特殊性，通常稱之為三對(duì)角方程（Tri-diagonalequation)三對(duì)角方程可以采用非常高效的追趕法（TDMA法）求解基于矩陣分解屬于必須掌握的內(nèi)容第42頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1.7差分方程的求解TDMA法Fortran源程序第43頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例求解下面的一維穩(wěn)態(tài)導(dǎo)熱問題：第44頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例求解區(qū)域的離散化：內(nèi)節(jié)點(diǎn)法：先劃分控制容積，在確定節(jié)點(diǎn)均勻網(wǎng)格：x=x將整個(gè)求解區(qū)域劃分為（N-2）個(gè)控制容積，N個(gè)節(jié)點(diǎn)（包括2個(gè)邊界節(jié)點(diǎn)）內(nèi)部節(jié)點(diǎn)的差分方程第45頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例其中，第46頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例注意：采用內(nèi)節(jié)點(diǎn)法劃分網(wǎng)格時(shí)，近邊界節(jié)點(diǎn)與其它內(nèi)部節(jié)點(diǎn)不盡相同，所以必須單獨(dú)考慮。123NN-1(x)w(x)e(x)e當(dāng)i＝2時(shí)(x)w=?xw=W=1所以，當(dāng)i＝2時(shí)第47頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例所以，當(dāng)i＝2時(shí)，第48頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例當(dāng)i＝3，4，。。。，N-2時(shí)，第49頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例同樣，當(dāng)i＝N-1時(shí)(x)e=?xe=E=N所以，當(dāng)i＝N-1時(shí)123NN-1(x)w(x)e(x)e第50頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例所以，當(dāng)i＝N-1時(shí)，第51頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例最后得到由（N－2）個(gè)方程構(gòu)成的方程組為求解上面的方程，即可得到（N－2）個(gè)未知數(shù)，即，T2,T3,T4,…….,TN-1。第52頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例注意：上面的方程組是非線性的，必須用迭代法求解求解方法：假定一個(gè)溫度分布：Ti，i＝1，2，3，。。。，N計(jì)算i，i＝1，2，3，。。。，N計(jì)算a,b,c,d用TDMA法求解方程組，得到新的溫度分布：Ti’計(jì)算：Max{abs(Ti-Ti’),i＝1,2,3,……,N}判斷：abs(Ti-Ti’)是否小于（精度要求）如果不能滿足精度要求，令Ti＝Ti’，重復(fù)上面的計(jì)算滿足精度要求：計(jì)算結(jié)束第53頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例希望大家用計(jì)算機(jī)完成上面的計(jì)算，并與下面的分析解結(jié)果比較：第54頁(yè)，共83頁(yè)，2023年，2月20日，星期四特別提示計(jì)算機(jī)實(shí)現(xiàn)的基礎(chǔ)地位關(guān)鍵：掌握循環(huán)變量的使用基礎(chǔ)：對(duì)算法清晰透徹的把握保障：細(xì)心細(xì)心再細(xì)心第55頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.2一維非穩(wěn)態(tài)導(dǎo)熱問題與穩(wěn)態(tài)問題的區(qū)別：增加了時(shí)間項(xiàng)數(shù)學(xué)上：常微分方程變?yōu)槠⒎址匠虜?shù)值方法上：與穩(wěn)態(tài)問題沒有本質(zhì)的區(qū)別控制方程:第56頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.2.1非穩(wěn)態(tài)項(xiàng)的處理時(shí)間坐標(biāo)的離散化tk=(k-1)t,k=1,2,3,……(46)非穩(wěn)態(tài)項(xiàng)的的離散化至少有三種方案：方案1、二階中心差分（Centraldifference)第57頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.2.1非穩(wěn)態(tài)項(xiàng)的處理非穩(wěn)態(tài)項(xiàng)的的離散化方案2、向前差分（Forwarddifference)方案3、向后差分（Backwarddifference)第58頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.2.2控制方程的離散化參照?qǐng)D示的節(jié)點(diǎn)組，對(duì)于任意一個(gè)CV，x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE圖1一維問題空間區(qū)域的離散化第59頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.2.2控制方程的離散化將方程（49）代入方程（48）第60頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.2.2控制方程的離散化假定節(jié)點(diǎn)間按線性分布，則x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE圖1一維問題空間區(qū)域的離散化而按式（47c），第61頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.2.2控制方程的離散化整理后得到，將式（51）代入方程（50），第62頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.2.2控制方程的離散化其中，第63頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.2.2控制方程的離散化將之與穩(wěn)態(tài)問題對(duì)比發(fā)現(xiàn)：非穩(wěn)態(tài)項(xiàng)的存在并沒有改變差分方程的基本形式發(fā)生變化的只是：aP中多出了a0P項(xiàng)bP中多出了a0PTk-1P項(xiàng)求解方法與穩(wěn)態(tài)問題一樣每前進(jìn)一個(gè)時(shí)間步長(zhǎng)都要求解一個(gè)方程組第64頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.2.3邊界條件的處理處理方法與穩(wěn)態(tài)問題類似元體能量平衡法：考慮CV內(nèi)能隨時(shí)間的變化考慮圖示的左邊界，按能量守恒，顯然有123eeqBx=0(x)1(x)21?(x)1Sq1其中Et是CV單位時(shí)間內(nèi)能的變化，顯然，第65頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理而，將q1和Et代入方程（54），有整理后得到，第66頁(yè)，共83頁(yè)，2023年，2月20日，星期四邊界條件的處理它與穩(wěn)態(tài)問題相比，多出了內(nèi)能變化項(xiàng)第67頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.3.4求解過程開始，t＝0t=t+t計(jì)算有關(guān)系數(shù)形成方程組，進(jìn)行求解輸出結(jié)果t<tcal是STOP否第68頁(yè)，共83頁(yè)，2023年，2月20日，星期四特別提示非穩(wěn)態(tài)問題的差分格式有顯式格式（Explicitformulations）時(shí)間步進(jìn)法：Time-marching不需要求解方程組程序簡(jiǎn)單，對(duì)計(jì)算機(jī)的內(nèi)存要求低穩(wěn)定性差隱式格式（Implicitformulations）每推進(jìn)一個(gè)時(shí)間步長(zhǎng)，都需要求解一個(gè)方程組程序復(fù)雜，要求計(jì)算機(jī)有較大的內(nèi)存穩(wěn)定性好第69頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.3多維非穩(wěn)態(tài)導(dǎo)熱問題的離散化4.3.1控制方程以直角坐標(biāo)系中的二維導(dǎo)熱為例第70頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.3多維非穩(wěn)態(tài)導(dǎo)熱問題的離散化4.3.2求解區(qū)域的離散化SNWEPsnwexy(y)n(y)s(y)-n(y)+s(x)w(x)e(x)-e(x)+w第71頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.3多維非穩(wěn)態(tài)導(dǎo)熱問題的離散化4.3.3控制方程的離散化基本思路與一維問題完全相同對(duì)方程（59）在控制容積CV上積分，第72頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.3.3控制方程的離散化其中SNWEPsnwexy(y)n(y)s(y)-n(y)+s(x)w(x)e(x)-e(x)+w第73頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.3.3控制方程的離散化節(jié)點(diǎn)間按線性分布，則代入（63）得到，第74頁(yè)，共83頁(yè)，2023年，2月20日，星期四4.3.3控制方程的離散化同樣SNWEPsnwexy(