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

計(jì)算傳熱學(xué)第4講擴(kuò)散方程的數(shù)值解NumericalSolutionofDiffusionEquations主要內(nèi)容一維穩(wěn)態(tài)問題的數(shù)值解一維非穩(wěn)態(tài)問題多維非穩(wěn)態(tài)問題的離散化差分方程的求解主要目的掌握用數(shù)值方法求解傳熱問題的整體步驟數(shù)值方法的計(jì)算機(jī)實(shí)現(xiàn)邊界條件的處理閱讀與作業(yè)閱讀要求：陶文銓《數(shù)值傳熱學(xué)》第4章作業(yè)：P124題4－1；P125題4－7完成課外作業(yè)第一題和第二題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)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ū)域的離散化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）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ū)域的離散化4.1一維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE圖1

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

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

(x)1也不等于

(x)2所以要對第一個(gè)內(nèi)部節(jié)點(diǎn)給予特別注意。邊界條件的處理例如，對于直角坐標(biāo)系，對于節(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！邊界條件的處理或?qū)懗桑吔鐥l件的處理附加源項(xiàng)法（Additionalsourcetermmethod)以內(nèi)節(jié)點(diǎn)法為例由方程（19）解出邊界節(jié)點(diǎn)上的待求變量T1，代入與第1個(gè)近邊界節(jié)點(diǎn)的差分方程（21），邊界條件的處理代入與第1個(gè)近邊界節(jié)點(diǎn)的差分方程（21），整理后得到，邊界條件的處理整理后得到，或者寫成，其中，邊界條件的處理其中Additionalsourceterm!邊界條件的處理對于第3類邊界條件，也可以做類似的處理，但是這時(shí)，請大家證明，邊界條件的處理附加源項(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)對應(yīng)的系數(shù)（aW）等于0特別提示邊界條件的處理是傳熱問題數(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)的相同4.1.7差分方程的求解將前面得到的差分方程（8）改寫為，或者簡單的寫成矩陣的形式，其中，4.1.7差分方程的求解4.1.7差分方程的求解與方程（33）對比，知，由方程（33）系數(shù)陣[A]的特殊性，通常稱之為三對角方程（Tri-diagonalequation)三對角方程可以采用非常高效的追趕法（TDMA法）求解基于矩陣分解屬于必須掌握的內(nèi)容4.1.7差分方程的求解TDMA法Fortran源程序4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例求解下面的一維穩(wěn)態(tài)導(dǎo)熱問題：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)的差分方程4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例其中，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í)4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例所以，當(dāng)i＝2時(shí)，4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例當(dāng)i＝3，4，。。。，N-2時(shí)，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)e4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例所以，當(dāng)i＝N-1時(shí)，4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例最后得到由（N－2）個(gè)方程構(gòu)成的方程組為求解上面的方程，即可得到（N－2）個(gè)未知數(shù)，即，T2,T3,T4,…….,TN-1。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é)束4.1.8計(jì)算機(jī)實(shí)現(xiàn)：算例希望大家用計(jì)算機(jī)完成上面的計(jì)算，并與下面的分析解結(jié)果比較：特別提示計(jì)算機(jī)實(shí)現(xiàn)的基礎(chǔ)地位關(guān)鍵：掌握循環(huán)變量的使用基礎(chǔ)：對算法清晰透徹的把握保障：細(xì)心細(xì)心再細(xì)心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ū)別控制方程: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)4.2.1非穩(wěn)態(tài)項(xiàng)的處理非穩(wěn)態(tài)項(xiàng)的的離散化方案2、向前差分（Forwarddifference)方案3、向后差分（Backwarddifference)4.2.2控制方程的離散化參照圖示的節(jié)點(diǎn)組，對于任意一個(gè)CV，x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE圖1

一維問題空間區(qū)域的離散化4.2.2控制方程的離散化將方程（49）代入方程（48）4.2.2控制方程的離散化假定節(jié)點(diǎn)間按線性分布，則x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE圖1

一維問題空間區(qū)域的離散化而按式（47c），4.2.2控制方程的離散化整理后得到，將式（51）代入方程（50），4.2.2控制方程的離散化其中，4.2.2控制方程的離散化將之與穩(wěn)態(tài)問題對比發(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í)間步長都要求解一個(gè)方程組4.2.3邊界條件的處理處理方法與穩(wěn)態(tài)問題類似元體能量平衡法：考慮CV內(nèi)能隨時(shí)間的變化考慮圖示的左邊界，按能量守恒，顯然有123eeqBx=0(x)1(x)21?(x)1Sq1其中Et是CV單位時(shí)間內(nèi)能的變化，顯然，邊界條件的處理而，將q1和Et代入方程（54），有整理后得到，邊界條件的處理它與穩(wěn)態(tài)問題相比，多出了內(nèi)能變化項(xiàng)4.3.4求解過程開始，t＝0t=t+t計(jì)算有關(guān)系數(shù)形成方程組，進(jìn)行求解輸出結(jié)果t<tcal是STOP否特別提示非穩(wěn)態(tài)問題的差分格式有顯式格式（Explicitformulations）時(shí)間步進(jìn)法：Time-marching不需要求解方程組程序簡單，對計(jì)算機(jī)的內(nèi)存要求低穩(wěn)定性差隱式格式（Implicitformulations）每推進(jìn)一個(gè)時(shí)間步長，都需要求解一個(gè)方程組程序復(fù)雜，要求計(jì)算機(jī)有較大的內(nèi)存穩(wěn)定性好4.3多維非穩(wěn)態(tài)導(dǎo)熱問題的離散化4.3.1控制方程以直角坐標(biāo)系中的二維導(dǎo)熱為例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)+w4.3多維非穩(wěn)態(tài)導(dǎo)熱問題的離散化4.3.3控制方程的離散化基本思路與一維問題完全相同對方程（59）在控制容積CV上積分，4.3.3控制方程的離散化其中SNWEPsnwexy(y)n(y)s(y)-n(y)+s(x)w(x)e(x)-e(x)+w4.3.3控制方程的離散化節(jié)點(diǎn)間按線性分布，則代入（63）得到，4.3.3控制方程的離散