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

計算傳熱學(xué)第講擴散方程的數(shù)值解第1頁，共83頁，2023年，2月20日，星期四主要內(nèi)容一維穩(wěn)態(tài)問題的數(shù)值解一維非穩(wěn)態(tài)問題多維非穩(wěn)態(tài)問題的離散化差分方程的求解第2頁，共83頁，2023年，2月20日，星期四主要目的掌握用數(shù)值方法求解傳熱問題的整體步驟數(shù)值方法的計算機實現(xiàn)邊界條件的處理第3頁，共83頁，2023年，2月20日，星期四閱讀與作業(yè)閱讀要求：陶文銓《數(shù)值傳熱學(xué)》第4章作業(yè)：P124題4－1；P125題4－7完成課外作業(yè)第一題和第二題第4頁，共83頁，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)系中的一維問題，無限長圓筒壁導(dǎo)熱問題）球坐標(biāo)系：A(x)＝x2（通過球壁的導(dǎo)熱）變截面肋片：A(x)第5頁，共83頁，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頁，共83頁，2023年，2月20日，星期四4.1一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解4.1.2源項的線性化在通常情況下，S＝S（T）線性化：

S＝Sc＋SpT

（2）其中，按負斜率源項原則，

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

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

(x)1也不等于

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