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

1、第四章 導(dǎo)熱問題的數(shù)值求解隨著計算機的普及應(yīng)用和性能的不斷改善，以及相關(guān)的數(shù)值計算方法的發(fā)展和應(yīng)用程序的開發(fā)，傳熱學(xué)數(shù)值計算方法作為數(shù)值求解傳熱問題的有效工具也得到了相應(yīng)的發(fā)展，利用計算機求解傳熱學(xué)問題愈來愈受到人們的普遍重視，而且在計算復(fù)雜傳熱問題中顯示出它的優(yōu)越性，因而成為傳熱學(xué)的一個重要的分支。數(shù)值傳熱的相關(guān)內(nèi)容也很自然地成為工程類學(xué)生學(xué)習傳熱學(xué)課程的不可缺少的部分。為了使學(xué)生能簡要地掌握傳熱學(xué)數(shù)值計算的基本方法，在這里我們以導(dǎo)熱問題為例對傳熱學(xué)數(shù)值計算方法做一個簡單的介紹。4-1導(dǎo)熱問題數(shù)值解概述在第二章和第三章中我們對較為簡單的導(dǎo)熱問題，如一維、二維簡單幾何形狀和邊界條件的穩(wěn)態(tài)導(dǎo)熱和

2、非穩(wěn)態(tài)導(dǎo)熱、以及一些特殊導(dǎo)熱問題，象通過肋片的導(dǎo)熱和忽略內(nèi)阻的集總導(dǎo)熱系統(tǒng)，進行了分析求解。然而對于一些更為復(fù)雜的導(dǎo)熱問題，如復(fù)雜的幾何形狀和邊界條件以及物性變化較大的情況，分析求解往往很復(fù)雜或者根本不可能。此時求解問題的唯一途徑是利用數(shù)值分析的辦法獲得數(shù)值解。x10T3x2x3T2T1T0x0 b d c axT圖4-1溫度場的有限差分表示數(shù)值求解通常是對微分方程直接進行數(shù)值積分或者把微分方程轉(zhuǎn)化為一組代數(shù)方程組再求解。這里要介紹的是后一種方法。如何實現(xiàn)從微分方程到代數(shù)方程的轉(zhuǎn)化又可以采用不同的數(shù)學(xué)方法，如有限差分法、有限元法和邊界元法等。作為一本入門的教材，這里僅向讀者簡要地介紹用有限差分

3、析方法從微分方程確立代數(shù)方程的處理過程。有限差分法的基本思想是把原來在時間和空間坐標中連續(xù)變化的物理量（如溫度、壓力、速度和熱流等），用有限個離散點上的數(shù)值集合來近似表示。有限差分的數(shù)學(xué)基礎(chǔ)是用差商代替微商（導(dǎo)數(shù)），而幾何意義是用函數(shù)在某區(qū)域內(nèi)的平均變化率代替函數(shù)的真實變化率。在圖4-1中可以看出有限差分表示的溫度場與真實溫度場的區(qū)別。圖中用T0、T1、T2表示連續(xù)的溫度場T；x為步長，它將區(qū)域的x方向劃分為有限個數(shù)的區(qū)域，x0、x1、x2，它們可以相等，也可以不相等。當x相等時，T1處的真實變化率a可以用平均變化率b、c或d來表示，其中b、c和d分別表示三種不同差分格式下的溫度隨時間的變化率

4、，即：b為向后差分格式；c為向前差分格式；d為中心差分格式。這種差分格式也可以推廣到高階微商的情形。對于二階微商的差分格式可以在一階差分格式的基礎(chǔ)上得出：。采用這樣的處理之后，反映溫度場隨時間、空間連續(xù)變化的微分方程就可以用反映離散點間溫度線性變化規(guī)律的代數(shù)方程來表示。當利用相應(yīng)的數(shù)學(xué)辦法求解這些代數(shù)方程組之后，我們就能獲得離散點上的溫度值。這些溫度值就可以近似表示溫度場的連續(xù)的溫度分布。從上面的分析不難看出，當我們要對導(dǎo)熱問題進行數(shù)值求解時一定要采取三個大的步驟，即研究區(qū)域的離散化；離散點（節(jié)點）差分方程的建立；節(jié)點方程（代數(shù)方程）的求解。下面我們將導(dǎo)熱問題的數(shù)值求解進行較為詳細的討論。4-

5、2研究區(qū)域溫度場的離散化SWENPSWENPSWENPxy圖4-2矩形長柱體截面區(qū)域離散化K-1時刻K時刻K+1時刻xy導(dǎo)熱問題的溫度場是假設(shè)為時間和空間的連續(xù)函數(shù)，當進行數(shù)值求解時首先要做的事情是在所研究的時間和空間區(qū)域內(nèi)把時間和空間分割成為有限大小的小區(qū)域，尤如地球被人為地劃分為不同的地域且冠以不同的名稱，時間被年、月、日和時、分、秒分割。如果在所分割的每一個時間間隔和空間區(qū)域內(nèi)均用同一個溫度值來表示，那么原來連續(xù)變化的溫度場就被一個離散的階躍變化的溫度分布所代替。這就是連續(xù)變化的溫度場離散化處理的基本思路。這里我們以一個矩形長柱體的非穩(wěn)態(tài)導(dǎo)熱過程為例來討論區(qū)域離散化問題。如果不考慮矩形長

6、柱體長度方向上的溫度變化，那么它是一個二維非穩(wěn)態(tài)導(dǎo)熱問題，圖4-2表示了長柱體矩形截面上區(qū)域離散化的情況。圖中可見，對于給定的空間區(qū)域，在x方向上的步長為x，在y方向上的步長為y，用它們作為空間尺度可以將矩形區(qū)域劃分成縱橫交錯的網(wǎng)格，交點稱為節(jié)點。然后以節(jié)點為中心，在兩個節(jié)點的中心處劃分界限，定出節(jié)點的控制面積，對于三維情況則為控制體積或控制容積，因而常在一般意義上稱之為節(jié)點的控制體?？刂企w的形狀是隨著坐標系的不同而改變的，這里的控制體是一個個的矩形面積。網(wǎng)格的步長在每一個方向上可以均勻劃分，也可以不均勻的劃分。因此，選用不同的步長和不同的劃分方法，可以將同一區(qū)域劃分出不同大小、不同數(shù)目的控制

7、區(qū)域，以及不同數(shù)目的節(jié)點數(shù)。顯然，隨著步長的不斷減小，節(jié)點數(shù)目的不斷增加，由節(jié)點溫度表示的離散的溫度場就會更加接近連續(xù)的溫度場，但計算工作量也會隨之增加。在時間方向上離散化的步長常用來表示，的選取也是可大可小的，也可以隨時間的進程而變化。顯然，無限小的時間步長亦會使得離散溫度變化接近連續(xù)的溫度改變，但隨之而來的是相應(yīng)的計算工作量的增加。4-3溫度場節(jié)點方程的建立為了得出所研究區(qū)域的節(jié)點溫度，必須建立相應(yīng)的節(jié)點方程。建立節(jié)點方程可以采用不同的方法，為了更好地理解節(jié)點方程的物理意義和掌握節(jié)點方程的建立方法，我們采用控制體熱平衡法來建立節(jié)點方程。下面我們以實例的形式介紹不同節(jié)點的節(jié)點方程的建立過程。

8、1 控制體的內(nèi)節(jié)點方程控制體熱平衡法建立節(jié)點方程的過程是將能量守恒方程應(yīng)用于控制體，建立該節(jié)點與周圍節(jié)點之間的能量平衡關(guān)系式，再利用傅立葉的導(dǎo)熱定律，最后獲得控制體節(jié)點溫度與周圍節(jié)點溫度之間的關(guān)系式。考察圖4-2中的節(jié)點P及其控制體，由能量平衡關(guān)系應(yīng)有，4-1式中，QW、QE、QS和QN分別為鄰近節(jié)點W、E、S和N通過傳導(dǎo)方式傳給節(jié)點P的熱流量；QV為單位時間控制體內(nèi)熱源的發(fā)熱量；為控制體單位時間內(nèi)熱能的增加量。由導(dǎo)熱傅立葉定律，在線性溫度分布的假設(shè)下，時刻K周圍節(jié)點傳給節(jié)點P的熱流量分別為：以及控制體的發(fā)熱流量,（qV為內(nèi)熱源強度，即單位時間單位體積的內(nèi)熱源發(fā)熱量。）控制體單位時間的內(nèi)能增加

9、量為，前者為時間上的向前差分，而后者為時間上向后差分。以上關(guān)系式中溫度T的上標為所在時刻，下標為所在空間位置。將以上關(guān)系式一并代入方程4-1中，且假設(shè)x=y，經(jīng)整理可以得出二維非穩(wěn)態(tài)導(dǎo)熱問題的內(nèi)節(jié)點的兩種差分格式的差分方程，即顯示差分格式4-2和隱示差分格式。4-3比較上面兩種差分格式可以看出，顯示差分格式最突出的優(yōu)點是節(jié)點溫度表達式的右邊只涉及K時刻的節(jié)點溫度值，那么只要知道這一時刻周圍節(jié)點的溫度值就可以求出該節(jié)點的下一時刻的溫度值；而隱示差分格式卻相反，溫度表達式的兩端都是同一K時刻的節(jié)點溫度值，這就意味著必須同時計算同一時刻所有節(jié)點的溫度值，即必須聯(lián)立求解K時刻所有節(jié)點的差分方程組，增大

10、計算工作量是顯而易見的。雖然顯示差分格式計算比較方便，但它卻存在著一個缺點，即計算式中a/x2值必須滿足一定的條件才不至于引起數(shù)值計算出現(xiàn)不收斂的問題，這在數(shù)值計算中稱為差分格式的不穩(wěn)定性。這里差分方程穩(wěn)定性的條件是方程4-2中的變量T前面的系數(shù)必須大于或等于零，分析一下差分方程中的各項系數(shù)，在TPK前的系數(shù)應(yīng)為，改寫成為。4-4此式稱為顯示差分格式的穩(wěn)定性判據(jù)，從中看出時間步長和空間步長是相互制約的。為了獲得較為精確的節(jié)點溫度值，空間步長x的選擇不能太小，按照穩(wěn)定性判據(jù)的要求勢必會使時間步長也要相應(yīng)地不能太大，因而必須在增加節(jié)點數(shù)目的同時增多時間間隔，從而使計算工作量加大。與顯示差分格式相反

11、，由于隱示差分格式的節(jié)點方程中沒有會使方程系數(shù)成為負值的系數(shù)項，因而不存在方程求解的不穩(wěn)定性的問題。也就是說，對于隱示差分格式，無論 a/x2中的x和取什么樣的數(shù)值，均不會出現(xiàn)數(shù)值計算結(jié)果的不收斂問題，因而是無條件穩(wěn)定的。這樣就使得我們能在滿足一定精度的情況下盡可能地加大時間步長或空間步長，亦可以在計算過程中隨意改變步長，而不必擔心會造成計算結(jié)果的不收斂。2 控制體的邊界節(jié)點方程在數(shù)值計算中所研究的區(qū)域的邊界條件是通過邊界節(jié)點的節(jié)點方程來反映的，因而邊界節(jié)點的差分方程的建立十分重要。這里同樣采用邊界節(jié)點的控制體熱平衡來確立邊界節(jié)點的差分方程。下面將具體討論一些典型的邊界節(jié)點的節(jié)點方程。對流換熱

12、邊界條件下的節(jié)點方程。 如圖4-3所示的邊界節(jié)點P，其控制體的熱平衡關(guān)系式為，4-5式中各項可以寫成：QWSNWPQNQSQCTyx圖4-3對流換熱邊界節(jié)點示意圖 將上述關(guān)系式一并代入4-5式中，經(jīng)整理可得出在二維非穩(wěn)態(tài)導(dǎo)熱過程中處于平直邊界上的邊界節(jié)點在對流換熱邊界條件下的節(jié)點差分方程（當x=y），即：對于顯示差分格式有；4-6對于隱示差分格式有。4-7分析公式4-6，可以得到其穩(wěn)定性判據(jù)為，不難看出條件要比內(nèi)節(jié)點的情形更加苛刻一些。絕熱邊界條件下的邊界節(jié)點方程。如圖4-4所示的絕熱邊界節(jié)點P，其控制體的熱平衡關(guān)系式為，4-8式中各項可寫成：圖4-4熱絕緣邊界節(jié)點示意圖QESNEPQNQSy

13、x絕熱邊界將上述關(guān)系式代入4-8式中，可以整理得到絕熱邊界條件下的節(jié)點方程（當x=y），即：對于顯示差分格式有；4-9對于隱示差分格式有。4-10對于其它邊界節(jié)點也可以采用上述的控制體熱平衡的辦法建立對應(yīng)的節(jié)點差分方程式。從以上的節(jié)點方程中可以發(fā)現(xiàn)兩個特別的系數(shù)項，即,它們分別是控制體（元體）的傅立葉數(shù)和畢歐數(shù)，其物理意義是前者表征控制體的導(dǎo)熱性能與熱儲蓄性能的對比關(guān)系，反映控制體溫度隨時間變化的動態(tài)特性；而后者則體現(xiàn)控制體和環(huán)境間的換熱性能與其導(dǎo)熱性能的對比關(guān)系。另外，還要一個常數(shù)項,它則反映了內(nèi)熱源產(chǎn)生的熱流引起的節(jié)點溫度升高值。下面我們歸納性地給出二維非穩(wěn)態(tài)導(dǎo)熱問題的兩種差分格式的節(jié)點方

14、程式，見表4-1。表4-1二維非穩(wěn)態(tài)導(dǎo)熱節(jié)點差分方程（式中Fo=a/x2,Bi=x/）節(jié)點內(nèi)型差分方程（x=y）穩(wěn)定性條件i,ji+1,ji-1,ji,j+1i,j-1內(nèi)節(jié)點 顯示差分格式Ti,jK+1=Fo(Ti+1K,j+Ti-1K,j+TiK,j-1+Ti,Kj+1)+(1-4Fo)Ti,jKFo1/4隱示差分格式Ti,jK=Fo(Ti+1K,j+Ti-1K,j+TiK,j-1+Ti,Kj+1)+ Ti,jK-1/ (1+4Fo)無條件i-1,ji,ji+1,ji-1,j平直絕熱邊界節(jié)點顯示差分格式Ti,jK+1=Fo(2Ti-1K,j+TiK,j-1+Ti,Kj+1)+(1-4Fo)T

15、i,jKFo1/4隱示差分格式Ti,jK=Fo(2Ti-1K,j+TiK,j-1+Ti,Kj+1)+ Ti,jK-1/ (1+4Fo)無條件i-1,ji,ji+1,ji-1,jT平直對流邊界節(jié)點顯示差分格式Ti,jK+1=Fo(2Ti-1K,j+TiK,j-1+Ti,Kj+1)+2FoBiTK +(1-2FoBi-4Fo)Ti,jKFo1/2(2+Bi)隱示差分格式Ti,jK=Fo(2Ti-1K,j+TiK,j-1+Ti,Kj+1)+ 2FoBiTK + Ti,jK-1/ (1+2FoBi+4Fo)無條件i-1,ji,ji+1,ji-1,jST平直輻射邊界節(jié)點顯示差分格式Ti,jK+1=Fo(

16、2Ti-1K,j+TiK,j-1+Ti,Kj+1)+(1-4Fo)Ti,jK-20S(Ti,jK)4-(TK)4/(cx)隱示差分格式Ti,jK (1-4Fo)+ 20S(Ti,jK)/(cx)=Fo(2Ti-1K,j+TiK,j-1+Ti,Kj+1)+ Ti,jK-1+ 20S (TK)4/(cx)無條件4-4節(jié)點方程的求解由上面的討論可以看出，對應(yīng)于離散溫度場的每一個節(jié)點均可以列出相應(yīng)的差分方程，這樣就可以得出與節(jié)點數(shù)目相同的一組代數(shù)方程組。當聯(lián)立求解這個代數(shù)方程組時，最后就可以得出每一個節(jié)點的溫度值。一般情況下，差分方程組是線性代數(shù)方程組，而線性代數(shù)方程組是可以用直接法和迭代法求解的。常

17、用的直接法有高斯消元法、列主元素消去法和矩陣求逆法，而迭代法常用的有高斯賽德爾迭代和超（欠）松弛迭代。下面我們僅介紹迭代法求解代數(shù)方程組的過程，并在例題中加以應(yīng)用。今有一線性代數(shù)方程組：迭代求解該方程組的思路為，尋找一個由（T1，T2，Tn）組成的列向量，使其收斂于某一個極限向量（T1*, T2*,，Tn*）,且該極限向量就是該方程的精確解。當這個線性代數(shù)方程組的系數(shù)項aii0（i=1,2,n）時，可將其改寫成迭代形式，有：以上各式可以用一個通用的形式來表示： i=1,2,n利用上式就可以進行迭代求解了，其步驟是，合理選擇（假設(shè)）各節(jié)點的初始溫度，將其作為第零次迭代的近似溫度值，記為Ti(0)

18、 (i=1,2,n);將Ti(0)代入上式的右端，得到第一次迭代的近似值Ti(1)；之后將Ti(1)再代入上式的右端，則得出第二次的近似值Ti(2)；如此反復(fù)進行下去，直至進行到K次，使相鄰的兩次近似解Ti(K+!)和Ti(K) (i=1,2,n)之間的偏差小于預(yù)先設(shè)定的小量時，即滿足Ti(K+1)Ti(K) (i=1,2,n)或(Ti(K+1)Ti(K)/Ti(K) (i=1,2,n)。此時各節(jié)點的溫度值T1(K), T2(K), Tn(K)已經(jīng)有足夠的精度用來表示代數(shù)方程組的解，從而可以結(jié)束方程求解的迭代過程。從上述的迭代過程不難發(fā)現(xiàn)，當我們用第零次迭代值去進行第一次迭代時，Ti(1)的值

19、已經(jīng)不斷地產(chǎn)生出來，當計算Tr(1 )時，到r-1的Ti(1)已經(jīng)求出。如果此時在計算Tr(1 )時涉及到的Ti(0) (i=1,2,r-1)全部用已求出的Ti(1) (i=1,2,r-1)代替，這勢必會加快迭代收斂的速度。這種改進后的迭代方法被稱之為高斯賽德爾迭代法。高斯賽德爾迭代法的方程組迭代形式為： i=1,2,n。歸納起來，高斯賽德爾迭代法的求解步驟可表述為，將代數(shù)方程組寫成迭代形式；設(shè)初始值經(jīng)迭代得出節(jié)點新值；有新值則去掉舊值，不斷以新?lián)Q舊，且在迭代過程中應(yīng)用；在迭代獲得滿足給定精度的節(jié)點溫度值后結(jié)束方程組的迭代。4-5導(dǎo)熱問題數(shù)值求解舉例圖4-5大平板非穩(wěn)態(tài)導(dǎo)熱數(shù)值解區(qū)域離散化示意

20、圖K+1KK-1i+1i-1ixTT0為了使讀者對導(dǎo)熱問題的數(shù)值求解過程有一個完整的認識，下面將以無限大平板對稱加熱的非穩(wěn)態(tài)導(dǎo)熱過程為例，闡明顯示差分格式和隱示差分格式的計算求解過程。有一厚度為0.12m的無限大平板，初始溫度為20,兩側(cè)表面同時受到溫度為150的流體加熱，流體與平板表面之間的對流換熱系數(shù)為24W/(m2),平板材料的導(dǎo)熱系數(shù)為0.24W/(m),而熱擴散系數(shù)為0.147×10-6m2/S。試計算平板溫度分布隨時間的變化。由于大平板是對稱受熱，利用溫度場的對稱性可以只對平板的一半進行數(shù)值計算。于是就可以給出該導(dǎo)熱問題的邊界條件為一側(cè)是平板與環(huán)境進行對流換熱的第三類邊界

21、條件，另一側(cè)（即對稱中心面）是處于絕熱狀態(tài)的第二類邊界條件。下面我們分別用顯示差分格式和隱示差分格式對其進行數(shù)值分析求解。1 顯示差分格式求解過程介紹針對該問題可從表4-1中簡化出中心節(jié)點和邊界節(jié)點的差分方程，即相應(yīng)的穩(wěn)定性條件在目前情況下有兩個，一個是內(nèi)節(jié)點和絕熱邊界節(jié)點應(yīng)滿足的F0=a/x21/2；一個是對流邊界節(jié)點應(yīng)滿足的F01/(2+Bi)。在這兩個條件中后一個要求要苛刻一些，故只需滿足后一個穩(wěn)定性條件即可。由穩(wěn)定性條件可見，空間步長x和時間步長是相互關(guān)聯(lián)的，當選定x之后也隨之被確定下來。在此我們選取x=0.006m，從穩(wěn)定性條件計算出76.3S,計算中選取=50S。圖4-6顯示差分格

22、式計算機程序框圖開始是結(jié)束CONTINUE打印AA,T(IJ,J),(IJ=1,N+1)AA=FLOAT(J-1)*DT計算T(1,J) 中心節(jié)點計算T(IJ,J)內(nèi)節(jié)點(2IJN計算T(N1,J)邊界節(jié)點DO 20 J=2,M+1輸入,H K, A, DT, TF, L, T0, N, MDX=L/N,BIO=HDX/K, FO=ADT/DX2FO1/(2BIO+2)T(I,J)=T0否計算程序的框圖如圖4-6所示，而程序中的變量標識符則說明于表4-2中。表4-2計算機程序中的標識符（含顯示和隱示差分格式用到的變量和常數(shù)）A平板材料的導(dǎo)溫系數(shù)K平板材料的導(dǎo)熱系數(shù)AB(I,J)節(jié)點溫度中間變量

23、I表示空間的下標變量AA計算的累計時間IJ表示空間的下標變量BIO網(wǎng)格畢歐數(shù)J表示時間的下標變量B1中間變量L平板厚度的一半B2中間變量M計算的時間間隔數(shù)目DT時間步長N計算的空間間隔數(shù)目DX空間步長T(I,J)節(jié)點溫度變量FO網(wǎng)格傅立葉數(shù)TF環(huán)境流體溫度H對流換熱系數(shù)T0平板初始溫度計算機程序采用FORTRAN語言編寫，詳細清單如下：REAL K，LDIMENSION T(50,50)READ(*,*) H,K,A,DT,TF,L,T0,N,MDX=L/FLOAT(N)BIO=(H*DX)/KFO=(A*DT)/(DX*2)IF(FO.GT.(1.0/(2.0+2.0*BIO) GOTO 6

24、0DO 10 I=1,N+110 T(I,1)=T0DO 20 J=2,M+1T(1,J)=2.0*FO*T(2, J-1)+(1.0-2.0*FO)*T(1,J-1)DO 30 IJ=2,NT(IJ,J)=FO*(T(IJ-1,J-1)+T(IJ+1,J-1)30 T(IJ,J)=T(IJ,J)+(1.0-2.0*FO)*T(IJ,J-1)T(N+1,J)=2.0*FO*(T(N,J-1)+TF*BIO)T(N+1,J)=T(N+1,J)+(1.0-2.0*FO-2.9*FO*BIO*T(N+1,J-1)AB=FLOAT(J-1)*DTWRITE(*,40) AB40 FORMAT(F9.2

25、)WRITE(*,50 (T(II,J),II=1,N+1)50 FORMAT(11F7.2)20 CONTINUE60 STOPEND利用以上程序經(jīng)過49次運算，即計算到非穩(wěn)態(tài)導(dǎo)熱過程進行2400S后的結(jié)果，將其顯示于表4-3中表中x為平板的坐標位置，x=0為對稱中心平面，x=0.06m為平板的外表面。從表中可見經(jīng)過2400S后，溫度為150的流體仍然在向平板加熱。進一步的計算表明，平板被均勻加熱到150的時間大約到96000S（26.7h），共需計算1920次。表4-3算例進行了49次迭代后所得的結(jié)果xm0.00.0060.0120.0180.0240.0300.0360.0420.048

26、0.0560.06T23.0223.5825.2528.3433.1740.1649.6962.0577.2995.21115.312 隱示差分格式計算機程序的編制從表4-1中同樣可以簡化出該算例的隱示差分格式的內(nèi)節(jié)點及邊界節(jié)點的差分方程式如下：TiK=Fo(Ti+1K+Ti-1K)+Ti-1K-1/(1+2Fo)（內(nèi)節(jié)點）；T1K=2FoT2K+T1K-1)/(1+2Fo)（對稱中心邊界節(jié)點）；TN+1K=(2FoTNK+2FoBiTK+TN+1K-1)/(1+2FoBi+2Fo)（對流換熱邊界節(jié)點）。由于隱示差分格式是無條件穩(wěn)定的，故在編制程序時不必設(shè)立穩(wěn)定性判別語句，但確要聯(lián)立求解所有的

27、節(jié)點方程組。這里我們采用高斯賽德爾迭代法，程序用FORTRAN語言編寫。程序中的EPS是精度控制的允許誤差，計算中取為0.01。若時間步長也取50S，計算到2000S時可以得到與表4-3所示的相近似的結(jié)果；若取2000S，只需計算48次即可接近平板被均勻加熱到150的穩(wěn)定工況。當然隨著的增大，計算結(jié)果的誤差也會相應(yīng)增大。計算程序框圖如圖4-7所示。開始是結(jié)束CONTINUE打印AA,T(IL,J),(IL=1,N+1)AA=FLOAT(J-1)*DT計算T(I,J) (1IN+1)DO 20 J=2,M+1輸入H，K, A, DT, TF, L, T0, EPS，N, M計算DX,BIO,FO

28、,B1,B2AB(I,J)-T(I,J)>EPST(I,1)=T0,1IN+1否圖4-7隱示差分格式計算機程序框圖T(I,1)=AB(I,J)=T0AB(I,J)=T(I,J)按上面的框圖編寫的計算機程序如下：REAL K, LDIMENSION T(50,50), AB(50,50)READ(*,*) H,K,A,DT,DF,L,T0,EPS,N,MDX=L/FLOAT(N)BIO=(H*DX)/KFO=(A*DT)/DX*2B1=1.0+2.0*FOB2=1.0+2.0*FO+2.0*FO*BIODO 10 I=1,N+110 T(I,1)=T0DO 20 J=2,M+1DO 30

29、IJ=1,N+1T(IJ,J)=T030 AB(IJ,J)=T035 T(1,J)=(T(1,J-1+2.0*FO*T(2,J)/B1DO 40 IK=2,N40 T(IK,J)=(T(N+1,J-1)+FO*(T(IK+1,J)+T(IK-1,J)/B1T(N+1,J)=(T(N+1,J-1)+2.0*FO*BIO*TF+2.0*FO*T(N,J)/B2DO 50 II=1,N+1IF(ABS(T(II,J)-AB(II,J)GTEPS) GOTO 6550 CONTINUEGOTO 7565 DO 70 IM=1.N=170 AB(IM,J)=T(IM,J)GOTO 3575 AA=FLO

30、AT(J-1)*DTWRITE(*,80)80FORMAT(F9.2)WRITE(*,90) (T(IL,J),IL=1,N+1)90 FORMAT(11F7.2)20 CONTINUE60 STOPEND思 考 題1.試簡要說明對導(dǎo)熱問題進行有限差分數(shù)值計算的基本思想與步驟。2.試說明用節(jié)點控制體熱平街法建立節(jié)點溫度差分方程的基本思想。3.推導(dǎo)導(dǎo)熱微分方程的步驟和過程與用熱平衡法建立節(jié)點溫度離散的差分方程的過程十分相似，什么前者得到的是溫度場的精確描寫，而由后者解出的卻是溫度場的近似分布。4.第三類邊界條件的邊界節(jié)點也可以采用將第三類邊界條件表達式中的一階導(dǎo)數(shù)用差分公式表示來建立其節(jié)點差分方

31、程。試比較這樣建立起來的差分方程與用熱平衡法建立起來的差分方程的異同與優(yōu)劣。5.什么是顯式格式，而什么是隱示差分格式？為什么顯式格式在計算中存在穩(wěn)定性問題而隱示差分格式卻不存在?6.用高斯-賽德爾迭代法求解代數(shù)方程組時是否一定可以得到收斂的解?在不能得出收斂的解時是否是因為初場的假設(shè)不合適造成的?習 題一般性數(shù)值計算 4-1試用數(shù)值計算證實，對方程組用高斯-賽德爾迭代法求解，其結(jié)果是發(fā)散的，并分析其原因。4-2試對附圖所示的常物性、無內(nèi)熱源的三維穩(wěn)態(tài)導(dǎo)熱問題用高斯-賽德爾迭代法計算t1，t2，t3，t4之值。4-24-3 4-3試對附圖所示的等截面直肋的穩(wěn)態(tài)導(dǎo)熱問題，用數(shù)值方法求解節(jié)點2、3的

32、溫度。圖中t0=85、tf=25、h=30W/(m2·K)。肋高H=4cm，縱剖面面積AL=4cm2，導(dǎo)熱系數(shù)=20W/(m·K)。 離散方程建立4-4試將直角坐標中的常物性、無內(nèi)熱源的三維非穩(wěn)態(tài)導(dǎo)熱微分方程化為顯式差分格式，并指出其穩(wěn)定性條件(xy)。4-5極坐標中常物性、無內(nèi)熱源的非穩(wěn)態(tài)導(dǎo)熱方程為 試利用本題附圖中的符號，列出節(jié)點(i，j)的差分方程式。4-6一金屬短圓柱在爐內(nèi)受熱后被豎直地移置到空氣中冷卻，底面可認為是絕熱的。為用數(shù)值法確定冷卻過程中柱體溫度的變化，取中心角為1rad的區(qū)域來研究(如本題附圖所示)。已知柱體表面發(fā)射率、自然對流表面?zhèn)鳠嵯禂?shù)從環(huán)境溫度t、

33、金屬的熱擴散率a,試列出圖中節(jié)點（1，1）、(M,l)、(M，n)及(M,N)的離散方程式。在r及z方向上網(wǎng)格是各自均分的。4-64-54-7一個三維物體的豎直表面受液體自然對流冷卻，為考慮局部表面?zhèn)鳠嵯禂?shù)的影響，表面?zhèn)鳠嵯禂?shù)采用h=c(t-tf)1.25來表示。試列出附圖所示的穩(wěn)態(tài)、無內(nèi)熱源物體邊界節(jié)點(M,n)的溫度方程，并對如何求解這一方程組提出你的看法。設(shè)網(wǎng)格均分。4-9在附圖所示的有內(nèi)熱源的三維導(dǎo)熱區(qū)域中，一個界面絕熱，一個界面等溫(包括節(jié)點4)，其余兩個界面與溫度為tf的流體對流換熱,h均勻，內(nèi)熱源強度為qv。試列出節(jié)點1、2、5、6、9、10的離散方程式。4-74-8一維穩(wěn)態(tài)導(dǎo)熱

34、計算4-9 一等截面直肋，高H，厚,肋根溫度為t0，流體溫度為tf,表面?zhèn)鳠嵯禂?shù)為h,肋片導(dǎo)熱系數(shù)為 。將它均分成4個節(jié)點(見附圖)，并對肋端為絕熱及為對流邊界條件(h同側(cè)面)的兩種情況列出節(jié)點2、3、4的離散方程式。設(shè)H=45mm，=l0mm，h=50W/(m2·K)=50W/(m·K),t0=100,tf=20，計算節(jié)點2、3、4的溫度(對于肋端的兩種邊界條件)。習題4-9附圖習題4-10附圖4-10復(fù)合材料在航空航天及化工等工業(yè)中日益得到廣泛的應(yīng)用。附圖所示為雙層圓筒壁，假設(shè)層間接觸緊密，無接觸熱阻存在。已知r1=l2.5mm，r2=l6mm，r3=18mm, 11=

35、40W/(m·K),2=120W/(m·K)，tf1=150,h1=l000W/(m2·K)，tf2=60,h2=380W/(m2·K)。試用數(shù)值方法確定穩(wěn)態(tài)時雙層圓筒壁截面上的溫度分布。一維非穩(wěn)態(tài)導(dǎo)熱計算4-11 一直徑為1cm、長4cm.的鋼制圓柱形肋片，初始溫度為25。其后，肋基溫度突然升高到200，同時溫度為25的氣流橫向掠過該肋片，肋片的端部及側(cè)面的表面?zhèn)鳠嵯禂?shù)均為l00W/(m2·K)。試將該肋片等分成兩段(見附圖)，并用有限差分法的顯式差分格式計算從開始加熱時刻起相鄰4個時刻上的溫度分布(以穩(wěn)定性條件所允許的時間間隔為計算依據(jù))。

36、已知=43W/(m·K)，a=l.333×10-5m2/S。(提示:節(jié)點4的離散方程可按端面的對流散熱與從節(jié)點3到節(jié)點4的導(dǎo)熱相平衡這一條件列出。)習題4-11附圖4-12 一厚為2.54cm的鋼板，初始溫度為650，后置于水中淬火，其表面溫度突然下降為93.5并保持不變。試用數(shù)值方法計算中心溫度下降到450所需的時間。已知。a=1.16x10-5m2/S。建議將平板8等分，取9個節(jié)點，并把數(shù)值計算的結(jié)果與按海斯勒圖計算的結(jié)果作比較。4-13一火箭燃燒器，殼體內(nèi)徑為400mm，厚l0mm，殼體內(nèi)壁上涂了一層厚2mm的包覆層?；鸺l(fā)動時，推進劑燃燒生成的溫度為3000的煙氣，經(jīng)燃燒器端部的噴管噴往大氣。大氣溫度為30。設(shè)包覆層內(nèi)壁與燃氣間的表面?zhèn)鳠嵯禂?shù)為2500W/(m2·K)，外殼表面與大氣間的表面?zhèn)鳠嵯禂?shù)為350W/(m2·K)，外殼材料的最高允許溫度為1500。試用數(shù)值法確定:為使外殼免受損壞，燃燒過程應(yīng)在多長時間內(nèi)完成。包覆材料的=0.3W/(m·K),a=2×l0-7m2/S;外殼的A=l0W/(