計算傳熱學第4講擴散方程的數值解NumericalSolutionofDiffusionEquations主要內容一維穩態問題的數值解一維非穩態問題多維非穩態問題的離散化差分方程的求解主要目的掌握用數值方法求解傳熱問題的整體步驟數值方法的計算機實現邊界條件的處理閱讀與作業閱讀要求：陶文銓《數值傳熱學》第4章作業：P124題4－1；P125題4－7完成課外作業第一題和第二題4.1一維穩態導熱問題的數值解控制方程：其中，A(x)是面積函數。定義如下：直角坐標系：A(x)＝1（無限大平板導熱問題）柱坐標系：A(x)＝x（極坐標系中的一維問題，無限長圓筒壁導熱問題）球坐標系：A(x)＝x2

（通過球壁的導熱）變截面肋片：A(x)4.1一維穩態導熱問題的數值解4.1.1求解區域的離散化x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE圖1

一維問題空間區域的離散化4.1一維穩態導熱問題的數值解4.1.2源項的線性化在通常情況下，S＝S（T）線性化：

S＝Sc＋SpT

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

Sp＝Sp(T*)0（3）4.1一維穩態導熱問題的數值解4.1.3控制方程的離散化將方程（1）兩邊通乘A（x），并對x從w到e積分：x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE圖1

一維問題空間區域的離散化4.1一維穩態導熱問題的數值解x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE圖1

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

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

(x)1也不等于

(x)2所以要對第一個內部節點給予特別注意。邊界條件的處理例如，對于直角坐標系，對于節點2，qBx=0(x)2(x)3123(x)1(x)2注意：CV2的左控制面w與節點W（節點1）重合，即與左邊界重合！控制面e！邊界條件的處理或寫成，邊界條件的處理附加源項法（Additionalsourcetermmethod)以內節點法為例由方程（19）解出邊界節點上的待求變量T1，代入與第1個近邊界節點的差分方程（21），邊界條件的處理代入與第1個近邊界節點的差分方程（21），整理后得到，邊界條件的處理整理后得到，或者寫成，其中，邊界條件的處理其中Additionalsourceterm!邊界條件的處理對于第3類邊界條件，也可以做類似的處理，但是這時，請大家證明，邊界條件的處理附加源項法的實質邊界節點消去法不僅能用于內節點網格，也能用于外節點網格實施方法：計算附加源項：Sc,ad，Sp,ad把附加源項計入該控制容積中的源項中令與邊界節點對應的系數（aW）等于0特別提示邊界條件的處理是傳熱問題數值計算最重要的環節之一元體能量平衡法的基礎地位盡可能采用外節點法劃分網格邊界節點消去法廣泛應用提高收斂速度盡可能采用均勻網格邊界節點的離散化方程在形式上與內部節點的相同4.1.7差分方程的求解將前面得到的差分方程（8）改寫為，或者簡單的寫成矩陣的形式，其中，4.1.7差分方程的求解4.1.7差分方程的求解與方程（33）對比，知，由方程（33）系數陣[A]的特殊性，通常稱之為三對角方程（Tri-diagonalequation)三對角方程可以采用非常高效的追趕法（TDMA法）求解基于矩陣分解屬于必須掌握的內容4.1.7差分方程的求解TDMA法Fortran源程序4.1.8計算機實現：算例求解下面的一維穩態導熱問題：4.1.8計算機實現：算例求解區域的離散化：內節點法：先劃分控制容積，在確定節點均勻網格：x=x將整個求解區域劃分為（N-2）個控制容積，N個節點（包括2個邊界節點）內部節點的差分方程4.1.8計算機實現：算例其中，4.1.8計算機實現：算例注意：采用內節點法劃分網格時，近邊界節點與其它內部節點不盡相同，所以必須單獨考慮。123NN-1(x)w(x)e(x)e當i＝2時(x)w=?xw=W=1所以，當i＝2時4.1.8計算機實現：算例所以，當i＝2時，4.1.8計算機實現：算例當i＝3，4，。。。，N-2時，4.1.8計算機實現：算例同樣，當i＝N-1時(x)e=?xe=E=N所以，當i＝N-1時123NN-1(x)w(x)e(x)e4.1.8計算機實現：算例所以，當i＝N-1時，4.1.8計算機實現：算例最后得到由（N－2）個方程構成的方程組為求解上面的方程，即可得到（N－2）個未知數，即，T2,T3,T4,…….,TN-1。4.1.8計算機實現：算例注意：上面的方程組是非線性的，必須用迭代法求解求解方法：假定一個溫度分布：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’，重復上面的計算滿足精度要求：計算結束4.1.8計算機實現：算例希望大家用計算機完成上面的計算，并與下面的分析解結果比較：特別提示計算機實現的基礎地位關鍵：掌握循環變量的使用基礎：對算法清晰透徹的把握保障：細心細心再細心4.2一維非穩態導熱問題與穩態問題的區別：增加了時間項數學上：常微分方程變為偏微分方程數值方法上：與穩態問題沒有本質的區別控制方程:4.2.1非穩態項的處理時間坐標的離散化tk=(k-1)t,k=1,2,3,……(46)非穩態項的的離散化至少有三種方案：方案1、二階中心差分（Centraldifference)4.2.1非穩態項的處理非穩態項的的離散化方案2、向前差分（Forwarddifference)方案3、向后差分（Backwarddifference)4.2.2控制方程的離散化參照圖示的節點組，對于任意一個CV，x(x)+e(x)-w(x)-e(x)+w(x)w(x)eweWPE圖1

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

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