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

1、5 熱傳導(dǎo)問題的數(shù)值方法5.1一維穩(wěn)態(tài)導(dǎo)熱一維穩(wěn)態(tài)導(dǎo)熱在直角坐標(biāo)系下的控制方程可表示為： (5-1)式中k為導(dǎo)熱系數(shù)，T是溫度，s是單位容積的熱產(chǎn)生率。 首先選定控制體和網(wǎng)格，如圖5.1所示，并對方程（5-1）在所選定的控制體進(jìn)行積分，即得： (5-2)圖5.1 控制體和網(wǎng)格然后進(jìn)行離散化。如果用分線段性分布來計算方程（5-2）中的微商，那么最終的方程為: (5-3)假設(shè)源項s在任一控制體中之值可以表示為溫度的線性函數(shù)，即，則導(dǎo)出的離散化方程為： (5-4)式中 (5-5)式(5-4)就是一維穩(wěn)態(tài)導(dǎo)熱方程的離散形式，系數(shù)aE和aW分別代表了節(jié)點(diǎn)P與E間及W與P間導(dǎo)熱阻力的倒數(shù)，它們的大小反映了

2、節(jié)點(diǎn)W和E處的溫度對P點(diǎn)的影響程度。式中的ke和kw是控制容積中的e和w界面上的當(dāng)量導(dǎo)熱系數(shù)。進(jìn)行計算時，物理參數(shù)值存儲在節(jié)點(diǎn)的位置上。為了確定ke和kw，還需規(guī)定由節(jié)點(diǎn)上的物理量來計算相應(yīng)界面上的量的方法。常用的方法由兩種，即算術(shù)平均法與調(diào)和平均法。1、算術(shù)平均法假定k與x呈線性關(guān)系，由P與E點(diǎn)的導(dǎo)數(shù)系數(shù)確定的公式為： （5-6）2、調(diào)和平均法利用傳熱學(xué)的基本公式可以導(dǎo)出確定界面上當(dāng)量導(dǎo)熱系數(shù)的調(diào)和平均公式?？刂迫莘e中P和E的導(dǎo)熱系數(shù)不相等，但界面上熱流密度應(yīng)該連續(xù)，則由Fourier定律可得： （5-7）而 則 （5-8）這就是確定界面上當(dāng)量導(dǎo)熱系數(shù)的調(diào)和平均公式，它反映了串聯(lián)過程熱阻的迭

3、加原則。3、兩種方法比較設(shè)，網(wǎng)格均分時， （5-9）即P和E兩點(diǎn)間的導(dǎo)熱阻力為，表明此時P和E間的熱阻主要是由導(dǎo)熱系數(shù)大的物體所決定，這顯然不符合傳熱學(xué)的基本原理。實際上，此時控制體E構(gòu)成了熱阻的主要部分。P和E間的熱阻可表示為： （5-10）從中可以看出與調(diào)和平均一致。令 ， 則 （5-11）因此，總體上看，調(diào)和平均要比算術(shù)平均更好一些。5.2 邊界條件與源項的處理式（5-4）導(dǎo)出的離散方程只適用于內(nèi)部節(jié)點(diǎn)。為了對某個特定的導(dǎo)熱問題進(jìn)行計算，還應(yīng)加上邊界條件。傳熱的邊界條件有三類，即（1）給定邊界溫度；（2）給出邊界熱流量；（3）通過換熱系數(shù)以及周圍流體溫度給定邊界熱流量。 如果是第一類邊界

4、條件，則就比較簡單。但如果是第二或第三類邊界，則需要對邊界條件進(jìn)行處理。在有限差分范圍內(nèi)，有兩種處理方法，即對邊界點(diǎn)補(bǔ)充代數(shù)方程和附加源項法。首先介紹邊界點(diǎn)補(bǔ)充代數(shù)方程的方法。先討論區(qū)域離散外節(jié)點(diǎn)法的情形。在邊界給定熱流量，如圖5-2所示，即給定第二類邊界條件，可表示為： （5-12）上式可化為： （5-13）即 (5-14)圖5-2 邊界節(jié)點(diǎn)控制體上式的截斷誤差為一階，而內(nèi)節(jié)點(diǎn)上如采用中心差分，則截斷誤差為二階。一般希望內(nèi)節(jié)點(diǎn)與邊界點(diǎn)離散方程截斷誤差等級保持一致，否則會影響計算結(jié)果的準(zhǔn)確度。為得出具有二階截斷誤差的公式，可采用虛擬點(diǎn)法。即在右邊界外虛設(shè)M+1點(diǎn)，這樣節(jié)點(diǎn)M就可視為內(nèi)節(jié)點(diǎn)，其一

5、階導(dǎo)數(shù)就可采用中心差分，即 （5-15）為消去，由一維穩(wěn)態(tài)含內(nèi)熱源的控制方程（5-1）可得到在 M點(diǎn)的離散形式，即 （5-16）從以上兩式消去可得 (5-17)其中，是節(jié)點(diǎn)M所代表的控制體的厚度。對于第三類邊界條件，以代入式(5-14)和(5-17)可得相應(yīng)于一階與二階截斷誤差的節(jié)點(diǎn)離散方程： （5-18） （5-19）用控制容積平衡法來推導(dǎo)。對圖5-2所示邊界節(jié)點(diǎn)的控制體作能量平衡，即可得 （5-20）解出即得(5-17)。為了使第二和第三類邊界條件的離散方程具有統(tǒng)一的形式，把邊界上的熱流密度表示成以下形式，其中為邊界溫度。a等于給定的熱流密度(進(jìn)入為+)，b=0相當(dāng)于第二類邊界條件；a =

6、，b=相當(dāng)于第三類邊界條件。，0 （5-21），0 （5-22）當(dāng)區(qū)域離散化采用第二種方法（即內(nèi)節(jié)點(diǎn)法）時，邊界節(jié)點(diǎn)可以看成是第一種區(qū)域離散法中當(dāng)邊界條件所代表的控制容積厚度趨近于0時的極限，式（5-22）變?yōu)椋? （5-23）其中為邊界節(jié)點(diǎn)與第一個內(nèi)節(jié)點(diǎn)之間的距離。上式雖然在形式上與區(qū)域離散外節(jié)點(diǎn)方法具有一階截斷誤差的公式一樣，但它卻是區(qū)域離散內(nèi)節(jié)點(diǎn)方法中具有二階截斷誤差的公式。離散化過程中源項的處理問題：如果源項是常數(shù)，則在離散方程的建立過程中不帶來任何困難。當(dāng)源項是新求解未知量的函數(shù)時，源項處理顯得十分重要，有時甚至是數(shù)值求解成敗的關(guān)鍵所在。目前應(yīng)用較廣泛的一種處理方法是把源項局部線性化

7、，即把源項寫成，其中為常數(shù)，為S隨T而變化的曲線在P點(diǎn)的斜率。源項的線性化處理需要注意以下幾點(diǎn)：1、當(dāng)源項為未知量函數(shù)時，線性化的處理比假定源項為常數(shù)更為合理。因為如果S=f(T)，那么把各個控制體中的S作為常數(shù)處理就是將上一次迭代計算得到的T*計算S，也就是說源項的計算相對于當(dāng)前T的計算有一個滯后，而如果按線性化處理，就不存在這個問題。2、線性化處理是建立線性代數(shù)方程所必須的。因為如果采用更高階的多項式，則所得到的離散化方程就不可能成為線性代數(shù)方程。3、為了保證代表方程迭代求解的收斂，要求0。因為，（下標(biāo)nb表示鄰點(diǎn)，為控制體的體積），線性代數(shù)迭代求解收斂的一個充分條件是對角占優(yōu)，這就要求0

8、。4、由代數(shù)方程迭代求解的公式為 (5-24)的大小影響到迭代過程中溫度的變化速度，的絕對值越大，所對應(yīng)系統(tǒng)的慣性越大，相鄰兩次迭代之間的變化越小，因而收斂速度下降，但有利于克服迭代過程的發(fā)散。例題：設(shè)有一導(dǎo)熱型方程，邊界條件為x=0，T=0；x=1，。試將該區(qū)域三等分，分別用外節(jié)點(diǎn)法和內(nèi)節(jié)點(diǎn)法求解該問題。解：（1）采用區(qū)域離散方法外節(jié)點(diǎn)法時，網(wǎng)格劃分如圖5-3(a)所示。內(nèi)節(jié)點(diǎn)采用中心差分，即。右端點(diǎn)采用一階截差時，離散方程為： 圖5-3 兩種離散方法的網(wǎng)格對節(jié)點(diǎn)： => 對： => 邊界條件： (x=0) (x=1) => 右端點(diǎn)采用二階截差時,，即 => 此問題精

9、確解為： 邊界條件采用不同離散表達(dá)式，對計算結(jié)果會產(chǎn)生影響，由表5-1可知，右端點(diǎn)采用二階截差時的結(jié)果遠(yuǎn)比采用一階截差時的結(jié)果準(zhǔn)確。 表5-1 采用外節(jié)點(diǎn)法時的數(shù)值解格式精確解0.22000.46480.7616一階截差0.24770.52290.8563二階截差0.21680.48670.7497（2）采用離散方法內(nèi)節(jié)點(diǎn)法求解此問題。網(wǎng)格選取如圖5-3(b)，則對于節(jié)點(diǎn)2，3，4和5的離散化方程分別為：數(shù)值解與精確解的比較見表5-2。表5-2 采用內(nèi)節(jié)點(diǎn)法時的數(shù)值解格式精確解0.10850.33770.60480.7616數(shù)值解0.10340.33720.60350.77025.3 一維非穩(wěn)

10、態(tài)導(dǎo)熱一維非穩(wěn)態(tài)導(dǎo)熱控制方程可寫為: (5-25)為了建立其離散化方程，在t, t+t時間間隔內(nèi)對控制體P作積分，為簡便起見，設(shè)與時間無關(guān)，則可得到:(5-26)為了完成上式的積分，需對上式右端項中T如何隨時間而變的方式作出選擇。常用的有顯式、隱式和Crank-Nicolson三種，可用下式來表示： （5-27）式中f是在0與1之間的加權(quán)因子。上述積分式最后可化為: （5-28）整理上式可得: （5-29）其中，式(5-29)就是一維非穩(wěn)態(tài)導(dǎo)熱的離散化方程，取f=0，1，以及1/2可依次得到顯式、隱式以及C-N格式。在直角坐標(biāo)系中，當(dāng)網(wǎng)格均勻時，無內(nèi)熱源、常物性導(dǎo)熱問題有三種格式，它們分別為：

11、1、顯式： (5-30)2、隱式： （5-31）3、C-N格式： （5-32）采用von Neumann分析方法可以證明，對于源項不隨時間而變的問題，對于式(5-29)，當(dāng)0.5f1時，絕對穩(wěn)定；而當(dāng)0f<1時，穩(wěn)定的條件則為，其中。當(dāng)采用顯式時，由式（5-30）可知，TP并不與TE和TW等其它未知數(shù)有關(guān)。任何f0的格式都將是隱式，TP必定和未知數(shù)TE和TW有關(guān)，并且必須求解一組聯(lián)立方程。相比之下，顯式就比較方便，但這一點(diǎn)將被它局限性所抵消。因為式(5-30)中的系數(shù)可能為負(fù)。事實上，為了使這個系數(shù)為正，必須使時間步長小到足以使大于。對于均勻?qū)?，均勻網(wǎng)格，可以表達(dá)為：。從中可以看出，當(dāng)

12、為了提高計算精度而減小時，只能采用更小的。至此，可以總結(jié)出保持離散化方程計算收斂和不失真的四個基本規(guī)則:1、控制容積界面上的連續(xù)性原則。在兩個控制容積的離散化方程中，通過界面的流量必須用相同的表達(dá)式來表示。如對于擴(kuò)散項，由于取了分段線性分布假設(shè)，因而保證了在控制容積界面上熱量和熱通量的連續(xù)性。如果如圖5-4所示的二次曲線來計算界面上的熱通量（-kdT/dx），則在公共界面上會造成梯度的不連續(xù)性。圖5-4 由二次曲線分布引發(fā)的熱通量不連續(xù)示意圖1-右邊斜率；2-左邊斜率2、系數(shù)為正原則。某個網(wǎng)格節(jié)點(diǎn)處的因變量值只是通過對流和擴(kuò)散的過程才受到相鄰節(jié)點(diǎn)上值的影響。如對于一維導(dǎo)熱問題，如果TE的增加必

13、然導(dǎo)致TP的增加，則要求TE的系數(shù)aE和TP的系數(shù)aP具有相同的符號。換言之，中心系數(shù)與相鄰系數(shù)的符號必須相同。3、源項的負(fù)斜率線性化原則。由式（5-5）可知，當(dāng)SP為正值時，中心節(jié)點(diǎn)的系數(shù)仍有可能為負(fù)值。只有當(dāng)SP不為正值時，才能完全避免這種情況的出現(xiàn)。因此，當(dāng)源項線性化處理時，。從計算方法上講，保持SP負(fù)值，使之不致產(chǎn)生計算過程中的不穩(wěn)定以及導(dǎo)致產(chǎn)生物理意義上的不真實解，這是至關(guān)重要的。4、相鄰節(jié)點(diǎn)系數(shù)加和原則?；疚⒎址匠掏话幸蜃兞康膶?dǎo)數(shù)項，于是，如果T代表因變量，則函數(shù)T與T+C（C為任一常數(shù)）兩者均滿足微分方程，微分方程所具有的這個特性必定要反映到與之相對應(yīng)的離散化方程中。因

14、此，當(dāng)TP以及所有的Tnb都增加同一常數(shù)時，離散化方程應(yīng)該仍然適用，這就要求=。但當(dāng)源項與T有關(guān)時，特別進(jìn)行線性化處理時，相鄰節(jié)點(diǎn)系數(shù)不服從此原則，這只是一個特例，如SP為零，仍滿足此原則。需要指出的是在數(shù)學(xué)上認(rèn)為是穩(wěn)定的初值問題的差分格式，未必能保證在所有的時間步長下均獲得具有物理意義的解。C-N格式就屬于這種情況。當(dāng)時，方程中的系數(shù)變成。對于均勻?qū)嵯禂?shù)和均勻網(wǎng)格問題，可以得出該系數(shù)為。從中可以看出只要時間步長不是足夠小，該系數(shù)就有可能變?yōu)樨?fù)，從而有可能出現(xiàn)物理上不真實的解。如果要求方程中的系數(shù)永遠(yuǎn)不為負(fù)的話，只有f為1才能確保這個要求。因此，全隱身格式(f=1)能滿足格式簡單且在物理上不

15、失真的要求。正是因為此，我們通常采用全隱式格式。但必須承認(rèn)，對小的時間步長，全隱式格式不如C-N格式精確。全隱式離散化方程可表示為: （5-33）式中, , , , 全隱式格式的主要原理是新值在整個時間步長占控制地位。5.4二維非穩(wěn)態(tài)導(dǎo)熱問題的全隱式格式在直角坐標(biāo)系下二維非穩(wěn)態(tài)導(dǎo)熱方程為: （5-34）在時間間隔t, t+t內(nèi)，對圖5-4所示的控制容積P做積分，假設(shè)在控制容積界面上熱流密度均勻，采用全隱式格式，則有：1、時間項： （5-35）2、擴(kuò)散項: （5-36） 圖5-5 二維直角坐標(biāo)系下的網(wǎng)格體系3、源項: （5-37）將（5-35）、（5-36）和（5-37）代入（5-34）并整理可

16、得: （5-38）其中: , , , , , 為控制體的容積，界面上當(dāng)量導(dǎo)熱系數(shù)按調(diào)解和平均計算。軸對稱的圓柱坐標(biāo)中，非穩(wěn)態(tài)導(dǎo)熱問題的控制方程為: （5-39）取一弧度為中心角所包含的范圍作為研究對象，網(wǎng)格體系見圖5-6所示，采用類似的推導(dǎo)方法可得形式上與上式完全相同的離散方程式，即 （5-40）, , , , , , ,圖5-6 二維柱坐標(biāo)系下的網(wǎng)格體系在極坐標(biāo)中, 非穩(wěn)態(tài)導(dǎo)熱問題的控制方程為： （5-41） （5-42）, , , , , ,圖5-7 極坐標(biāo)系下的網(wǎng)格體系5.5三維問題的離散化方程在z方向上再加兩個相鄰近節(jié)點(diǎn)T和B(頂和底),即可完成三維情況的離散化，離散化方程的形式與二維

17、的類似，即 （5-43）式中, , , , , , , 5.6多維導(dǎo)熱問題邊界條件的處理對于第一類邊界條件的問題，由于邊界溫度為已知值，這樣就不必補(bǔ)充邊界節(jié)點(diǎn)方程，求解的區(qū)域僅限于內(nèi)部節(jié)點(diǎn)。對于二維、三維問題的第二、三類邊界條件問題，如果也能做到這樣，則不僅能節(jié)省計算時間，而且有助于編程。目前比較成功的處理方法就是附加源項法。所謂附加源項法就是把由第二類或第三類邊界條件所規(guī)定的進(jìn)入或?qū)С鲇嬎銋^(qū)域的熱量作為與邊界相鄰的控制容積的當(dāng)量源項。以直角坐標(biāo)中的情形為例，如圖5-8所示。與邊界相鄰的控制容積中的節(jié)點(diǎn)為P。此控制容積的離散化方程可表達(dá)為： （5-44）圖5-8 二維邊界控制體為了在的代數(shù)方程中不出現(xiàn)未知的邊界溫度，就需要利用已知的邊界條件把消去，為此將上式作如下變換：由于其中為邊界節(jié)點(diǎn)的導(dǎo)熱系數(shù)，為進(jìn)入該控制容積的熱流密度，以進(jìn)入為正。這樣關(guān)于P點(diǎn)的方程可化為： （5-45）對于第二類邊界條件，為已知，故可把它與b組成一個新項：同時， 這就是說，對第二類邊界條件，如果把作為與邊界相鄰的控制容積的附加常數(shù)源項，記為，同時令，則所得的離散方程既