第五章__導(dǎo)熱問題的數(shù)值解

1、 高等傳熱學內(nèi)容4第一章 導(dǎo)熱理論和導(dǎo)熱微分方程4第二章 穩(wěn)態(tài)導(dǎo)熱 4第三章 非穩(wěn)態(tài)導(dǎo)熱 4第四章 凝固和熔化時的導(dǎo)熱 4第五章 導(dǎo)熱問題的數(shù)值解 4第六章 對流換熱基本方程 4第七章 層流邊界層的流動與換熱 4第八章 槽道內(nèi)層流流動與換熱 4第九章 湍流流動與換熱 4第十章 自然對流 4第十一章 熱輻射基礎(chǔ) 4第十二章 輻射換熱計算 4第十三章 復(fù)合換熱 第五章第五章 導(dǎo)熱問題的數(shù)值解導(dǎo)熱問題的數(shù)值解 4數(shù)值求解方法是以離散數(shù)學為基礎(chǔ)，以計算機為工具的一種求解方法。與各種分析求解的方法相比，它在應(yīng)用方面表現(xiàn)出很大的適用性，對于處理諸如非線性、復(fù)雜幾何形狀、復(fù)雜邊界條件的問題以及藕合的偏微分方

2、程組等都能較好地解決。尤其是隨著計算機技術(shù)的迅速發(fā)展，數(shù)值計算的精度和速度都大大提高。因此，稍微復(fù)雜一些的導(dǎo)熱問題現(xiàn)在主要依靠數(shù)值法求解。4目前用于求解偏微分方程的數(shù)值方法主要有有限差分法(FDM)、有限元法(FEM) 1,2和邊界元法(BEM)。有限元法起源于固體力學和結(jié)構(gòu)分析。與有限差分法相比，有限元法在整個區(qū)域內(nèi)單元的劃分比較隨意，更易于處理不規(guī)則幾何形狀的問題。有限差分法涉及的數(shù)學基礎(chǔ)及其表達式比較簡單，本章將主要介紹有限差分法及其在導(dǎo)熱中的應(yīng)用。 5-1 導(dǎo)數(shù)的有限差分近似表達式導(dǎo)數(shù)的有限差分近似表達式 4有限差分的數(shù)學基礎(chǔ)是用差商代替微商，即用有限差分代替導(dǎo)數(shù)。若f(x)是連續(xù)函數(shù)

3、，則它的導(dǎo)數(shù)定義為4 (5-1-1) 4在這里，df /dx稱為微商(導(dǎo)數(shù))，f/x 稱為有限差商。微商是有限差商當x趨于零時的極限。在x沒有達到零以前，f/x只是df/dx的近似，而兩者的差值就是用差商代替微商的偏差。 4利用泰勒級數(shù)可以得到一個函數(shù)的導(dǎo)數(shù)的有限差分近似表達式，如圖5-1所示 4朝前差分 (5-l-7 )4朝后差分 (5-1-8) 00()( )limlimxxdtf xxf xfdxxx ()()()dff xxf xxOxdxx( )()()dff xf xxOxdxx5-1 導(dǎo)數(shù)的有限差分近似表達式導(dǎo)數(shù)的有限差分近似表達式 圖5-1 函數(shù)的微商與差商 5-1 導(dǎo)數(shù)的有限

4、差分近似表達式導(dǎo)數(shù)的有限差分近似表達式 4式中的O(x)是用差商代替微商的誤差，即舍去泰勒級數(shù)的高階項引起的誤差，稱為截斷誤差。這里O 表示數(shù)量級，O(x)表示該截斷誤差是與x 同量級的小量。4一階導(dǎo)數(shù)的中心差分表達式4 (5-1-9) 4注意中心差分格式的截斷誤差是(x)2量級的。4二階導(dǎo)數(shù)的中心差分格式4 (5-l-10)4前面介紹的一階導(dǎo)數(shù)的差分格式都只用到兩個點的函數(shù)值，為了得到更高精度的差分格式，可以利用更多點的函數(shù)值。 二階導(dǎo)數(shù)的朝前差分：4 (5-1-14) 4將以上方法推廣，可以得到更高精度的各種有限差分表達式。 2()()()2dff xxf xxOxdxx2222()2 (

5、 )()()()d ff xxf xf xxOxdxx22222 ( )5 ()4 (2)(3)()()d ff xf xxf xxf xxOxdxx 4這里以二維穩(wěn)態(tài)導(dǎo)熱為例討論有限差分法的應(yīng)用。為了進行數(shù)值求解，首先要把求解的區(qū)域離散化。有限差分法要求來用正交的網(wǎng)格，在直角坐標系中就是矩形網(wǎng)格，如圖5-2 所示。網(wǎng)線的交點稱為節(jié)點，在區(qū)域內(nèi)的節(jié)點稱為內(nèi)節(jié)點，落在邊界上的節(jié)點為邊界節(jié)點。網(wǎng)線之間的距離稱為步長。x 方向的步長記為x，y 方向的步長記為y。為了便于計算，需要對節(jié)點編號。若節(jié)點P的坐標為(x, y)，xix，yiy，則用(i, j )表示節(jié)點的位置，節(jié)點的溫度t(x，y)則記為t

6、i,j。4數(shù)值求解的第二步是建立節(jié)點方程，即對每一個節(jié)點寫出一個代數(shù)方程。建立節(jié)點方程的方法之一是從微分形式出發(fā)，也就是在微分方程或邊界條件中用差商近似微商，可以得到以節(jié)點溫度為未知量的代數(shù)方程。 5-2 穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 5-2 穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 圖5-2 差分區(qū)域的離散化、節(jié)點和步長 4有內(nèi)熱源的二維穩(wěn)態(tài)問題，在常物性條件下可由泊松方程描述：4 (5-2-1)4區(qū)域內(nèi)的所有點，包括內(nèi)節(jié)點(i, j)都應(yīng)滿足以上的方程。把節(jié)點(i, j)處的二階偏導(dǎo)數(shù)用對應(yīng)的差商來近似，有4 (5-2-2)4 (5-2 -3) 4把以上兩式舍去截斷誤差并代入式(5-2

7、-1)，就得到內(nèi)節(jié)點的差分方程4 (5-2-4) 5-2 穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 22220Vqttxy21,1,222,2()()iji jiji jtttd tOxdxx2,1,1222,2()()i ji ji ji jtttd tOydyy1,1,1,1, ,22220()()iji jiji ji ji jV i jttttttqxy4如果采用正方形的網(wǎng)格，即xy ，且無內(nèi)熱源(qV0) ，則式(5-2-4 ) 簡化為 4 ( 5-2-5 ) 4注意到該差分方程的截斷誤差是O(x)2 + (y)2。利用邊界條件也可導(dǎo)得邊界節(jié)點的差分方程。若邊界節(jié)點(1, j)和(n,

8、j)分別滿足絕熱邊界條件和第三類邊界條件： 4 (5-2-6) 4 (5-2-7) 4把式(5-l-7 )的朝前差分格式代入式(5-2-6)，可得邊界節(jié)點(1，j)的差分方程4 (5-2-8) 5-2 穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 ,1,1,1,114i jijiji ji jttttt0,0txx,()0ftxLh ttx1,2,0jjtt4注意到，該節(jié)點方程的截斷誤差為x量級，與內(nèi)節(jié)點差分方程的截斷誤差相比，其精度低于一個量級。為了使各節(jié)點方程的精度能夠均衡，可以利用“虛節(jié)點”的概念對邊界節(jié)點進行適當?shù)奶幚怼?假設(shè)在邊界節(jié)點(l, j)的外面還有個節(jié)點(0, j )，且 。注意該處

9、的邊界還是維持絕熱，而節(jié)點(1, j)可以按內(nèi)節(jié)點處理，得到4 (5-2-9)4很顯然，該節(jié)點方程的截斷誤差是(x)2量級的。4對于第三類邊界條件的節(jié)點(n, j ) ，出于同樣的考慮，也可以假設(shè)有一個虛節(jié)點(n+1, j )。這樣，一方面邊界節(jié)點(n, j )可以按內(nèi)節(jié)點處理，得到4 (5-2-10) 5-2 穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 0,2,jjtt1,2,1,1,1124jjji jtttt,1,1,1,114n jnjnjn jn jttttt4另一方面該節(jié)點還應(yīng)滿足邊界條件式(5-2-14)。采用截斷誤差為O(x)2的中心差分格式的式(5-l-9)代入邊界條件，得4 (5

10、-2-11) 4聯(lián)立式(5-2-10)、(5-2-11)消去虛節(jié)點的溫度，整理得4 (5-2-12)4其他邊界節(jié)點，如兩個邊界相交處的節(jié)點，也可以按同樣的思路加以處理。 5-2 穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 1,1,02njnjn jftth ttx,1,1,12224n jnjn jn jfh xh xttttt4建立節(jié)點方程的另一個方法是“元體熱平衡法”。這一方法不是從偏微分方程的邊值問題出發(fā)考慮問題，而是以積分形式對一些有限的元體從能量守恒關(guān)系建立方程，因此方程的物理概念更加清晰。兩者處理問題的步驟基本相同，對大多數(shù)情況也得到形式相同的差分方程。但是，從微分方程出發(fā)的方法要求溫度

11、場在節(jié)點處的函數(shù)值和一階、二階導(dǎo)數(shù)值連續(xù)，而元體熱平衡法則沒有這樣的要求。因此，在處理位于邊界上、復(fù)合介質(zhì)的結(jié)合面上、有接觸熱阻處等溫度分布不光滑、不連續(xù)的節(jié)點時，元體熱平衡法就更為方便和靈活。4用元體熱平衡法建立節(jié)點方程，第一步仍是把區(qū)域離散化，需要把整個區(qū)域劃分為有限個小單元體。在每個元體內(nèi)近似地取一點的溫度代表整個元體的溫度，該點也稱為節(jié)點。從原則上說，節(jié)點位置的選取可以是任意的，但從建立節(jié)點方程的方便考慮，內(nèi)部單元的節(jié)點總是取在它的中心，邊界單元的節(jié)點則取在邊界上。 5-2 穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4把能量守恒關(guān)系應(yīng)用于每個元體，在穩(wěn)態(tài)導(dǎo)熱的情況下，從所有相鄰的元體導(dǎo)入的

12、熱量和該元體本身的發(fā)熱量之代數(shù)和應(yīng)為零。這樣，對于圖5-3 所示的內(nèi)部元體，有以下的關(guān)系：4 ( 5-2-13 )4在根據(jù)傅里葉定律計算相鄰兩個元體的導(dǎo)熱量時，假設(shè)兩個節(jié)點間的溫度是線性分布的。這實際上也是采用了以一階差商近似偏導(dǎo)數(shù)的概念。例如：4 ， (5-2-14) 4代入式(5-2-13)并整理，得內(nèi)部元體的節(jié)點方程為4與從微分方程出發(fā)得到的節(jié)點方程完全相同。 5-2 穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 0EWNSVQQQQqx y 1,iji jEttQyx1,iji jWttQyx1,1,1,1, ,22220()()iji jiji ji ji jV i jttttttqxy5-

13、2 穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4下面再以圖5-4中位于角上的單元(1, l)為例，說明怎樣用元體熱平衡法處理外邊界從屬于兩種不同的邊界條件時的混合邊界元體。參見圖5-4 ，邊界元體(l, l)的一個邊界處于對流換熱，另一個邊界受到給定熱流qw的作用。用元體熱平衡法可以寫出4 (5-2-15)4如果x y ，且qV0上式可改寫為 4 (5-2-16)4各種不同組合的混合邊界元體以及處于復(fù)合介質(zhì)交界處的元體都可按同樣的思路列出節(jié)點方程。4如果微分方程和邊界條件都是線性的，則得到的內(nèi)節(jié)點和邊界節(jié)點的差分方程都將是線性代數(shù)方程。由全部節(jié)點的差分方程構(gòu)成一個線性代數(shù)方程組，其中未知數(shù)的個數(shù)與

14、方程的個數(shù)相等。求解這一線性代數(shù)方程組就可以得到節(jié)點處的溫度。 5-2 穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 2,11,11,21,11,1()022224fwVttttyxyxx yh ttqqxy 1,12,11,22wffqyh yh yttttt5-2 穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4下面介紹求解線性代數(shù)方程組的方法。4從物理意義上說，物體具有穩(wěn)態(tài)溫度分布的條件是，單位時間里在全部邊界上流出的熱量應(yīng)等于物體內(nèi)部發(fā)出的熱量；或者，在沒有內(nèi)熱源的情況下，在全部邊界上流出的總熱量應(yīng)等于零。如果在全部邊界上給定熱流邊界條件，則穩(wěn)態(tài)導(dǎo)熱問題可能無解(不滿足上述條件時)，或溫度場的解不確定

15、(滿足上述條件時)。所以，全部給定第一類邊界條件的穩(wěn)態(tài)導(dǎo)熱的提法是不充分的。4線性代數(shù)方程組的數(shù)值解法可分為直接法和迭代法兩類。直接法是指在沒有舍入誤差的條件下經(jīng)過有限次的運算即可得到方程組的精確解的方法。高斯消元法和系數(shù)矩陣求逆的方法就是常用的直接法。對于階數(shù)不是很高的方程組，采用直接法是很有效的。所謂迭代法，是把求解方程組的問題化為構(gòu)造一個無限序列，來逐步逼近所求的精確解，因而在有限步的迭代中將只能得到方程組的近似解。注意到用有限差分法求解穩(wěn)態(tài)導(dǎo)熱問題時，不管總節(jié)點數(shù)的多少，每個節(jié)點方程中未知數(shù)的個數(shù)對二維問題不超過5個，對三維問題不超過7個。因此，產(chǎn)生的代數(shù)方程組的系數(shù)知陣中包含有大量的

16、零元素，這樣的系數(shù)矩陣稱為稀疏矩陣。對于穩(wěn)態(tài)導(dǎo)熱的有限差分法得到的線性方程組，用迭代法求解時收斂較快，計算程序也比較簡單，因此被廣泛采用。 5-2 穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 5-2 穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4用迭代法求解線性代數(shù)方程組的基本思路如下。有線性代數(shù)方程組 4 (5-2-17) 4可以把其中任一方程改寫為4 (5-2-18) 11 112 21121 122 2221 12 2n nn nnnnn nna ta ta tba ta ta tba ta ta tbLLL L L LL11(),1,2,niiij jjiij itba tinaL5-2 穩(wěn)態(tài)導(dǎo)熱

17、的數(shù)值分析穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4任意給定各節(jié)點的溫度值 作為解的初次近似值。把這些近似值代入式(5-2-18 )的右端，可計算得到解的第一次近似 值 。把第一次近似值再次代入式(5-2-18)的右端，可得到解的第二次近似值。一般地說，在得到解的第k 次近似值后，可由式(5-2- 18 )得到第k1次近似值：4 (5-2-19) 4在實際計算中，當然只能迭代有限的次數(shù)，而取迭代的結(jié)果 作為方程組的近似解。通常以相鄰的兩次迭代值之間的差值不大于預(yù)先設(shè)定的小量作為結(jié)束迭代的判據(jù)，即 (0)1,2,itinL(1)1,2,itinL(1)( )11(),1,2,nkhiiij jjiij itba t

18、inaL( )1,2,kitinL( )(1),(1,2, )kkiittinL5-2 穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4以上介紹的迭代過程稱為簡單迭代，或同步迭代。它在計算k1 次迭代值時全部采用各節(jié)點的第k 次近似值。實際上，在進行第kl 次迭代的過程中，有一部分先計算的節(jié)點的溫度已經(jīng)得到了第k十1 次近似值。顯然它們優(yōu)于第k次近似值，應(yīng)優(yōu)先于第k次近似值而被采用。因此，可對式(5-2-19)進行修改，得4 (5-2-21) 4這種迭代過程稱為異步迭代，或高斯一賽德爾迭代。與同步迭代相比，它的收斂速度加快。此外，由于異步迭代不需要用兩套工作單元存放節(jié)點溫度的老值和新值，而可以只用一套

19、工作單元，因此可以節(jié)省計算機的存儲單元。 1(1)(1)( )111(),1,2,inkkkiiij jij jjj iiitba ta tina L5-2 穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4當節(jié)點的個數(shù)很多時，通常收斂速度較慢，需要的迭代次數(shù)很大。為此提出了各種改進收斂速度的措施。其中之一是在高斯一賽德爾迭代基礎(chǔ)上的超松弛迭代法。在每一輪迭代中，首先按高斯一賽德爾迭代計算得到節(jié)點溫度的新值 ，然后再用一個適當選取的松弛因子來改善這一結(jié)果，即4 (5-2-22)4當1 時，就相當于普通的高斯賽德爾迭代；0l 稱為欠松弛，不利于收斂；l 稱為超松弛。對于特定的問題，可以找出一個收斂速度最快

20、的的值，稱為最佳松弛因子*。最佳松弛因子的范圍為l *2 。4在處理節(jié)點個數(shù)很多的問題時，另一個加快收斂的方法是先采用大網(wǎng)格進行粗略的計算。用計算結(jié)果對小網(wǎng)格的節(jié)點溫度進行插值，作為迭代的初值。由于選取的迭代初值比較接近精確解，迭代過程能大大加快。 (1)kT(1)( )(1)1,1,2,kkkiiittTin5-3 非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4首先，求解區(qū)域的離散化不僅涉及空間坐標的離散化，而且需要將時間離散化。例如，對于一維的溫度響應(yīng)t = t(x, )，取適當?shù)目臻g步長x 和時間步長，數(shù)值分析就是把一個求連續(xù)函數(shù)t = t(x, )的問題轉(zhuǎn)化為求特定的節(jié)點上在特定的時間間

21、隔上 的溫度值 離散問題。4以一維非穩(wěn)態(tài)導(dǎo)熱為例，在直角坐標系中，區(qū)域內(nèi)的溫度滿足以下偏微分方程：4 (5-3-1 ) 4如果采用朝前差分，p時刻對時間的偏導(dǎo)數(shù)近似為4 ( 5-3-2 ) ,xi xp (1,2, ;0,1,2,)pitin pLL22ttax1ppiittt5-3 非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4對坐標x的二階偏導(dǎo)數(shù)用p 時刻的中心差分近似為4 ( 5-3-3 )4在微分方程中用有限差分代替導(dǎo)數(shù)，經(jīng)整理可以得到內(nèi)節(jié)點的朝前差分形式的節(jié)點方程4 ( 5-3-4 ) 4其截斷誤差為 。因為在p+1 時刻的節(jié)點溫度 可以由p 時刻的已知溫度 、 和 根據(jù)式(5-3-4

22、)直接求得，所以這樣的差分格式稱為顯式格式。 211222pppiiittttxx11122()1 2ppppiiiiaattttxx2()OOx1pit1pitpit1pit5-3 非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4如果p + l 時刻對時間的偏導(dǎo)數(shù)采用朝后差分近似為4 ( 5-3-5)4對坐標x 的二階偏導(dǎo)數(shù)用p + 1時刻的中心差分近似為4 ( 5-3-6)4在微分方程中用有限差分代替導(dǎo)數(shù)，經(jīng)整理可以得到內(nèi)節(jié)點的朝后差分形式的節(jié)點方程4 (5-3-7) 1ppiittt111211222pppiiittttxx11111(12 )ppppiiiiftf tftt5-3 非穩(wěn)態(tài)導(dǎo)

23、熱的數(shù)值分析非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4其中， ，其截斷誤差也是 ?，F(xiàn)在，在一個節(jié)點方程中有3 個未知的節(jié)點溫度。因此，為了確定p + l 時刻各節(jié)點的溫度，需要求解一組聯(lián)立代數(shù)方程。這樣的差分格式稱為隱式格式。4如果對空間變量x的偏導(dǎo)數(shù)的有限差分取顯式格式式(5-3-3 )和隱式格式式(5-3-6 )的算術(shù)平均值來表示，則得到另一種差分格式，稱為克蘭克-尼科爾森(Crank -Nicolson )格式，或六點差分格式。內(nèi)節(jié)點的差分方程變?yōu)? (5-3-8) 4或利用以上定義的符號f，上式可整理為4 (5-3-9 ) 2afx2()OOx1111111122222ppppppppiiiiiiiit

24、tttttttaxx11111112(1)2(1)ppppppiiiiiiftf tftftf tft5-3 非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4在這一差分格式中，對時間的差分相當于在 時刻的中心差分，因此這一差分格式的截斷誤差是 ，精度高于朝前差分和朝后差分格式。4對于邊界節(jié)點，同樣地利用元體熱平衡法建立節(jié)點方程更加方便。注意，與穩(wěn)態(tài)導(dǎo)熱問題不同，在熱平衡關(guān)系中增加了一項因元體溫度改變而引起的內(nèi)能的變化。例如，對于圖5-5 中給出的一維問題的邊界節(jié)點n，滿足對流邊界條件。4 (5-3-10) 4如果溫度對時間的變化率采用朝前差分計算，則對于單位面積的邊界元體的熱平衡方程可寫作4 (5

25、-3-11) 1 2p2()OOx,()0ftxh ttx 11()2ppppppnnnnfnttttxch ttx5-3 非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 圖5-5 時空區(qū)域的離散化 5-3 非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4整理得4 (5-3-12) 4其中， ， ，這一節(jié)點方程是顯式格式。如果采用朝后差分格式或克蘭克-尼科爾森格式，則可導(dǎo)得同一個邊界元體的節(jié)點方程分別為4 (5-3-13)4 (5-3-14)4以上兩個節(jié)點方程都屬隱式格式。其他邊界條件下的邊界節(jié)點，二維問題的邊界節(jié)點、二維混合邊界節(jié)點都可以按同樣的思路導(dǎo)得不同差分格式的節(jié)點方程。 112 ()(1 22

26、)ppppnnfntf tbtffb t2afxh xb11112(1 22)2ppppnnfnftffb tfbtft1111(1)()(1)pppppnnffnftffb tfb ttffb t5-3 非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 45-3-1 差分格式的穩(wěn)定性差分格式的穩(wěn)定性4用有限差分法求解非穩(wěn)態(tài)導(dǎo)熱問題時存在差分格式穩(wěn)定性的問題。4用有限差分法求導(dǎo)熱問題的數(shù)值解時，存在兩類誤差。其一是用有限差分近似導(dǎo)數(shù)時帶來的截斷誤差。減小差分步長可以減小截斷誤差。此外，數(shù)值計算只能進行到有限位有效數(shù)字，在每一步計算中都會帶來舍入誤差。因此，計算得到的代數(shù)方程組的解通常不可能是它的精確解

27、，而差分方程組的精確解一般地也不等同于微分方程的精確解。4對于非穩(wěn)態(tài)導(dǎo)熱問題(較一般地說，對于求解拋物線型的“擴散方程”) ，由于問題本身的特性，在數(shù)值求解時存在解的穩(wěn)定性問題。數(shù)值求解非穩(wěn)態(tài)導(dǎo)熱問題的過程是隨時間逐步推進的。即：先由初始溫度分布算出時刻的溫度分布，再算出2時刻的溫度分布，然后漸次推進。這就有一個以前各時刻計算的溫度產(chǎn)生誤差對后面的計算有什么影響的問題。如果誤差趨向于不斷積累和放大，就會引起解的不穩(wěn)定性。因此，保證解的穩(wěn)定性在實際計算中是十分重要的。 5-3 非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4這里僅介紹一種稱為傅里葉級數(shù)法的比較簡單的方法。這種方法的不足之處是僅適用于

28、線性齊次的差分方程，因此不能給出邊界條件對穩(wěn)定性的影響。4以一維問題為例，如果在時刻計算得到的溫度分布為 4 (5-3-15) 4其中 (x, )為精確的溫度分布，以(x, )為溫度分布的誤差。當然，若 = 0，(x, 0)可以看作初始溫度分布的誤差，由前面章節(jié)的分離變量法的分析中知，假設(shè)函數(shù)在0，區(qū)間上逐段可積，就可以在該區(qū)間上用特征函數(shù)展開成廣義傅里葉級數(shù)，級數(shù)的每一項都滿足一組特定的邊界條件。作為一個例子，假定傅里葉級數(shù)有以下的形式：4 (5-3-16) ( , )( , )( , )t xxx1( , )cosmmmxAx5-3 非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4由于微分方程

29、是線性的，級數(shù)中每一項對其后各時刻解的影響可以疊加，因此在這里只需研究級數(shù)中的任一項。以下去掉求和符號和下標m，且常數(shù)A 也可不加考慮。實際上，這就是可以簡化地考慮一個簡諧波形式的分布誤差隨計算過程的傳遞。進行離散化處理，則節(jié)點i, i+1, i-1 處的誤差分別為4 (5-3-17) 4經(jīng)過一次步進運算，到時刻，誤差分布變?yōu)?。記4 (5-3-18) 4稱為誤差的放大倍數(shù)。如果某種差分格式在運算過程中能始終保證 ，則誤差分布在計算過程中不會被放大，這種差分格式就是穩(wěn)定的；反之，若 ，則差分格式就是不穩(wěn)定的。 11cos()cos(1)cos()cos()sin()sin()cos(1)cos

30、()cos()sin()sin(pipipii xixi xxi xxixi xxi xx1,1,2,piin1pipi115-3 非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4由于我們討論的差分方程是線性齊次的，而計算的溫度分布看成是精確的溫度分布和分布誤差的疊加，因而和兩者應(yīng)各自滿足差分方程。對于內(nèi)節(jié)點的朝前差分格式，即式(5-3-4) ，有 4 (5-3-19) 4把的節(jié)點值式(5-3-17 )代入，可得4 (5-3-20) 111(12 )ppppiiiifff1cos()cos()sin()sin()(1 2 )cos()cos()cos()sin()sin()2cos()cos()(

31、1 2 )cos()pifi xxi xxfi xfi xxi xxfi xxfi x122cos()(1 2 )1 21 cos()1 4sin2pipifxffxxf 5-3 非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4根據(jù)以上討論的穩(wěn)定性條件，為了使顯式差分方程(5-3-4 )穩(wěn)定，需有4 (5-3-22) 4這一穩(wěn)定性條件意味著在選定空間步長x以后，時間步長不能隨意選取，必須滿足式(5-3-22)規(guī)定的條件，否則差分格式將不穩(wěn)定。4同樣地，對于隱式的朝前差分格式式(5-3-7)可以導(dǎo)得其誤差放大倍數(shù) 4 (5-3-23) 4對于任何f 0，很明顯 的穩(wěn)定性判據(jù)總是成立的。因此隱式差分格

32、式式(5-3-7 )是無條件穩(wěn)定的。 214sin12xf212afx2114sin2xf15-3 非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4對于克蘭克-尼科爾森差分格式，可以導(dǎo)得4 (5-3-24)4可知，對于任何f 0 , 的穩(wěn)定性判據(jù)總是成立的。因此克蘭克-尼科爾森差分格式也是無條件穩(wěn)定的。4以上是用數(shù)學的方法，根據(jù)溫度分布誤差在步進運算的過程中是否被放大為判據(jù)來確定某一差分格式的穩(wěn)定性。另一方面，從非穩(wěn)態(tài)導(dǎo)熱的物理概念上也可以得到一些有用的結(jié)論。還是以一維問題為例，對內(nèi)節(jié)點的顯式格式有221 2sin212sin2xfxf1111()(12 )ppppiiiitf ttf t5-3

33、非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4該差分方程表明，i節(jié)點在p+l 時刻的溫度可以用i-1、i、i+1三個節(jié)點在p 時刻的溫度的加權(quán)平均來表示。要使這一計算符合非穩(wěn)態(tài)導(dǎo) 熱過程的物理意義，必須使三個節(jié)點溫度 , , 的系數(shù)(即權(quán)重)都不小于零。在這一具體情況下，即要求1-2f 0，亦即f 1/2。如果不滿足這一條件，即的系數(shù)小于零，那就意味著i 節(jié)點在p時刻的溫度越高，在p + l 時刻它的溫度就反而越低，進而在p + 2 時刻它的溫度又會變得更高。這樣就會引起節(jié)點溫度的不穩(wěn)定振蕩，顯然不符合導(dǎo)熱的規(guī)律。同樣地，對于邊界節(jié)點的顯式格式差分方程，如4 (5-3-25) 1pitpit1p

34、it2122afbx5-3 非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 45-3-2 一維差分方程組的求解一維差分方程組的求解 追趕法追趕法4在用隱式格式，包括朝后差分和六點差分格式求解非穩(wěn)態(tài)導(dǎo)熱問題時，所導(dǎo)得的代數(shù)方程組在每一次步進計算中都必須聯(lián)立求解。但是對于一維非穩(wěn)態(tài)導(dǎo)熱問題，可以采用一種稱為“追趕法”的更為簡捷的直接求解方法。從以上推導(dǎo)的隱式格式的節(jié)點方程可以看到，每一個內(nèi)節(jié)點的方程中只包含有3個未知溫度，邊界節(jié)點的方程中最多只有2個未知溫度。在方程組的系數(shù)矩陣中只有主對角線及相鄰的兩條對角線上有非零元素，方程組可以整理成式以下的形式：4 (5-3-26) 1 11 212 12 22 32111i ii ii iin nn nnbtctda tb tc tdatbtctda tb td5-3 非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析非穩(wěn)態(tài)導(dǎo)熱的數(shù)值分析 4這樣的系數(shù)矩陣稱為“三對角線矩陣”，可以用一種比較簡單的直接法，即追趕法求解。整個求解分為兩段，第一部分是“順追趕”求系數(shù)。對于方程組中的第一個方程，可以改寫為4 (5-3-27) 4其中4 ， (5-3-28) 4把式(5-3-27)代入方程組的第二個方程可以消去其中的t1，經(jīng)整理得 4 (5-3-29) 11