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

1、第四章 導(dǎo)熱問題的數(shù)值解法1、重點內(nèi)容： 掌握導(dǎo)熱問題數(shù)值解法的基本思路； 利用熱平衡法和泰勒級數(shù)展開法建立節(jié)點的離散方程。2、掌握內(nèi)容：數(shù)值解法的實質(zhì)。3、了解內(nèi)容：了解非穩(wěn)態(tài)導(dǎo)熱問題的兩種差分格式及其穩(wěn)定性。由前述3可知，求解導(dǎo)熱問題實際上就是對導(dǎo)熱微分方程在定解條件下的積分求解，從而獲得分析解。但是，對于工程中幾何形狀及定解條件比較復(fù)雜的導(dǎo)熱問題，從數(shù)學(xué)上目前無法得出其分析解。隨著計算機技術(shù)的迅速發(fā)展，對物理問題進行離散求解的數(shù)值方法發(fā)展得十分迅速，并得到廣泛應(yīng)用，并形成為傳熱學(xué)的一個分支計算傳熱學(xué)（數(shù)值傳熱學(xué)），這些數(shù)值解法主要有以下幾種：（1） 有限差分法 （2）有限元方法 （3）邊

2、界元方法數(shù)值解法能解決的問題原則上是一切導(dǎo)熱問題，特別是分析解方法無法解決的問題。如：幾何形狀、邊界條件復(fù)雜、物性不均、多維導(dǎo)熱問題。分析解法與數(shù)值解法的異同點：1、 相同點：根本目的是相同的，即確定 t=f(x，y，z)； 。2、 不同點：數(shù)值解法求解的是區(qū)域或時間空間坐標系中離散點的溫度分布代替連續(xù)的溫度場；分析解法求解的是連續(xù)的溫度場的分布特征，而不是分散點的數(shù)值。 §41 數(shù)值求解的基本思路及穩(wěn)態(tài)導(dǎo)熱內(nèi)節(jié)點離散方程的建立一、 解法的基本概念1、 實質(zhì)對物理問題進行數(shù)值解法的基本思路可以概括為：把原來在時間、空間坐標系中連續(xù)的物理量的場，如導(dǎo)熱物體的溫度場等，用有限個離散點上的

3、值的集合來代替，通過求解按一定方法建立起來的關(guān)于這些值的代數(shù)方程，來獲得離散點上被求物理量的值。該方法稱為數(shù)值解法。這些離散點上被求物理量值的集合稱為該物理量的數(shù)值解。2、基本思路：數(shù)值解法的求解過程可用框圖4-1表示。由此可見：1）物理模型簡化成數(shù)學(xué)模型是基礎(chǔ)；2）建立節(jié)點離散方程是關(guān)鍵；3）一般情況微分方程中，某一變量在某一坐標方向所需邊界條件的個數(shù)等于該變量在該坐標方向最高階導(dǎo)數(shù)的階數(shù)。二、 數(shù)值求解的步驟如圖4-2（a），二維矩形域內(nèi)無內(nèi)熱源、穩(wěn)態(tài)、常物性的導(dǎo)熱問題采用數(shù)值解法的步驟如下：1、 建立控制方程及定解條件控制方程：是指描寫物理問題的微分方程針對圖示的導(dǎo)熱問題，它的控制方程（

4、即導(dǎo)熱微分方程）為： （a）邊界條件：x=0時，x=H時，當y=0時，當y=W時，2、 區(qū)域離散化（確立節(jié)點）用一系列與坐標軸平行的網(wǎng)格線把求解區(qū)域劃分成若干個子區(qū)域，用網(wǎng)格線的交點作為需要確定溫度值的空間位置，稱為節(jié)點(結(jié)點)，節(jié)點的位置用該節(jié)點在兩個方向上的標號m ，n表示。相鄰兩節(jié)點間的距離稱步長。x, y每個節(jié)點都可以看成是以它為中心的一個小區(qū)域的代表把節(jié)點代表的小區(qū)域稱為元體（又叫控制容積），如圖4-2(b)。3、 建立節(jié)點物理量的代數(shù)方程（離散方程）節(jié)點上物理量的代數(shù)方程稱離散方程。其過程如下：1） 首先劃分各節(jié)點的類型；2） 其次，建立節(jié)點離散方程；3） 最后，代數(shù)方程組的形成。

5、對節(jié)點(m,n)的代數(shù)方程，當x=y時，有： （b）4、 設(shè)立迭代初場代數(shù)方程組的求解方法有直接解法與迭代解法，傳熱問題的有限差分法中主要采用迭代法。采用迭代法求解時，需對被求的溫度場預(yù)先設(shè)定一個解，這個解稱為初場，并在求解過程中不斷改進。5、 求解代數(shù)方程組求解時遇到的問題：線性；非線性；收斂性等。如圖4-2（b），除m=1的左邊界上各節(jié)點的溫度已知外，其余(M-1)N個節(jié)點均需建立離散方程，共有(M-1)N個方程，則構(gòu)成一個封閉的代數(shù)方程組。1）線性代數(shù)方程組：代數(shù)方程一經(jīng)建立，其中各項系數(shù)在整個求解過程中不再變化；2）非線性代數(shù)方程組：代數(shù)方程一經(jīng)建立，其中各項系數(shù)在整個求解過程中不斷更

6、新。3）是否收斂判斷：是指用迭代法求解代數(shù)方程是否收斂，即本次迭代計算所得之解與上一次迭代計算所得之解的偏差是否小于允許值。關(guān)于變物性(物性為溫度的函數(shù))導(dǎo)熱問題，建立的離散方程，四個鄰點溫度的系數(shù)不是常數(shù)，而是溫度的函數(shù)。在迭代計算時，這些系數(shù)應(yīng)不斷更新，這是非線性問題。6、 解的分析通過求解代數(shù)方程，獲得物體中的溫度分布，根據(jù)溫度場應(yīng)進一步計算通過的熱流量，熱應(yīng)力及熱變形等。因此，對于數(shù)值分析計算所得的溫度場及其它物理量應(yīng)作詳細分析，以獲得定性或定量上的結(jié)論。三、穩(wěn)態(tài)導(dǎo)熱中位于計算區(qū)域內(nèi)部的節(jié)點離散方程的建立方法1、基本概念1） 內(nèi)節(jié)點：位于計算區(qū)域內(nèi)部的節(jié)點，稱內(nèi)節(jié)點；2） 差分格式：差

7、商中的差分可以用向前、向后、中心差分表示的格式稱差分格式。2、基本方法方法：泰勒級數(shù)展開法；熱平衡法，以下分述之。1）泰勒級數(shù)展開法如圖4-3所示，以節(jié)點(m,n)處的二階偏導(dǎo)數(shù)為例，對節(jié)點(m+1,n)及(m-1,n)分別寫出函數(shù)t對(m,n)點的泰勒級數(shù)展開式：對(m+1,n)： （a）對（m-1,n）： （b）（a）+（b）得： 變形為的表示式得：上式是用三個離散點上的值計算二階導(dǎo)數(shù) 的嚴格表達式，其中： 稱截斷誤差，誤差量級為，即表示未明確寫出的級數(shù)余項中 的最低階數(shù)為2。在數(shù)值計算時，用三個相鄰節(jié)點上的值近似表示二階導(dǎo)數(shù)的表達式即可，則相應(yīng)的略去。于是得： （4-1a）同理： （4-

8、1b）根據(jù)導(dǎo)熱問題的控制方程(導(dǎo)熱微分方程) 得： (4-2)若x=y 則有：2) 平衡法：其本質(zhì)是傅里葉導(dǎo)熱定律和能量守恒定律的體現(xiàn)。對每個元體，可用傅里葉導(dǎo)熱定律寫出其能量守恒的表達式。如圖4-3所示，元體在垂直紙面方向取單位長度，通過元體界面(w,e,n,s)所傳導(dǎo)的熱流量可以對有關(guān)的兩個節(jié)點根據(jù)傅里葉定律寫出：從節(jié)點(m-1,n)通過界面w傳導(dǎo)到節(jié)點(m,n)的熱流量： （a）同理：通過界面e,n,s傳導(dǎo)給節(jié)點（m,n）的熱流量： （b） （c） （d）對元體(m,n).根據(jù)能量守恒定律可知： （4-3）其中，規(guī)定：導(dǎo)入元體（m,n）的熱流量為正； 導(dǎo)出元體（m,n）的熱流量為負。說明

9、： 上述分析與推導(dǎo)是在笛卡兒坐標系中進行的； 熱平衡法概念清晰，過程簡捷； 熱平衡法與§22建立微分方程的思路與過程一致，但不同的是前者是有限大小的元體，后者是微元體。§42 穩(wěn)態(tài)導(dǎo)熱邊界節(jié)點離散方程的建立及代數(shù)方程的求解（1） 對于第一類邊界條件的導(dǎo)熱問題，所有內(nèi)節(jié)點的離散方程組成一個封閉的代數(shù)方程組，即可求解；（2） 第二類或第三類邊界條件的導(dǎo)熱問題，所有內(nèi)節(jié)點的離散方程組成的代數(shù)方程組是不封閉的，因未知邊界溫度，因而應(yīng)對位于該邊界上的節(jié)點補充相應(yīng)的代數(shù)方程，才能使方程組封閉，以便求解。一、用熱平衡法導(dǎo)出典型邊界點上的離散方程假設(shè)物體具有內(nèi)熱源 (不必均勻分布)，而且邊

10、界上有向該元體傳遞的熱流密度：1、位于平直邊界上的節(jié)點如圖所示4-4邊界節(jié)點(m,n)只能代表半個元體，若邊界上有向該元體傳遞的熱流密度為，據(jù)能量守恒定律對該元體有： （4-4a）若時，則： （4-4b）2、外部角點如圖4-5所示，二維墻角計算區(qū)域中，該節(jié)點外角點僅代表1/4個以為邊長的元體。假設(shè)邊界上有向該元體傳遞的熱流密度為，則據(jù)能量守恒定律得其熱平衡式為： （4-5a）若時，則： （4-5b）3、內(nèi)部角點：如圖4-5所示內(nèi)部角點代表了3/4個以為邊界長的元體。同理得： (4-6a)時，則： （4-6b）4、討論有關(guān)的三種情況：1）若是絕熱邊界則，即令上式即可。2）若時流入元體，取正，流出

11、元體，取負使用上述公式。3）若屬對流邊界則：，將代入上式即可。時，則：對于平直邊界： （4-7）對外角點： （4-8）對內(nèi)角點：（4-9）其中，無量綱數(shù)是以網(wǎng)格步長為特征長度的畢渥數(shù)，即為。二、代數(shù)方程的求解方法1、直接解法：通過有限次運算獲得精確解的方法，如：矩陣求解，高斯消元法。2、迭代法：先對要計算的場作出假設(shè)（設(shè)定初場），在迭代計算中不斷予以改進，直到計算前的假定值與計算結(jié)果相差小于允許值為止的方法，稱迭代計算收斂。目前應(yīng)用較多的是：1）高斯賽德爾迭代法：每次迭代計算，均是使用節(jié)點溫度的最新值。2）用雅可比迭代法：每次迭代計算，均用上一次迭代計算出的值。設(shè)有一三元方程組：其中 （i=1

12、,2,3；j=1,2,3）及均不為零。采用高斯賽德爾迭代法的步驟：（1）將三元方程變形為迭式方程：（2）假設(shè)一組解（迭代初場），記為：，并代入迭代方程求得第一次解，同理求得改進值，（注：再次計算應(yīng)該用新值）如：（3）以新的初場重復(fù)計算，直到相鄰兩次迭代值之差小于允許值，則稱迭代收斂，計算終止。三、判斷迭代收斂的準則1、2、3、其中上角標k,k+1 表示迭代次數(shù)，為第k次迭代計算所的計算區(qū)域中的最大值。若計算區(qū)域中有t0時，應(yīng)采用3判斷之。說明： 1）對于一個代數(shù)方程組，若選用的迭代方式不合適，有可能導(dǎo)致發(fā)散，即稱迭代過程發(fā)散；2）對于常物性導(dǎo)熱問題，組成的差分方程組，迭代公式的選擇應(yīng)使一個迭代

13、變量的系數(shù)總是大于或等于該式中其他變量系數(shù)絕對值的代數(shù)和，此時，結(jié)果一定收斂。這一條件數(shù)學(xué)上稱主對角線占優(yōu)（對角占優(yōu)）； ， ， 3）采用熱平衡法導(dǎo)出差分方程時，若每一個方程都選用導(dǎo)出該方程中心節(jié)點的溫度作為迭代變量，則上述條件必滿足，迭代一定收斂。§43 非穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解法由前可知：非穩(wěn)態(tài)導(dǎo)熱和穩(wěn)態(tài)導(dǎo)熱二者微分方程的區(qū)別在于控制方程中多了一個非穩(wěn)態(tài)項，其中擴散項的離散方法與穩(wěn)態(tài)導(dǎo)熱一樣。本節(jié)重點討論：（1）非穩(wěn)態(tài)項離散的方法； （2）擴散項離散時所取時間層的不同對計算帶來的影響。一、一維非穩(wěn)態(tài)導(dǎo)熱時間空間區(qū)域的離散化1、基本概念如圖4-8所示，x為空間坐標，為時間坐標。1）時

14、間步長：指從一個時間層到下一個時間層的間隔。2）節(jié)點（n, i）表示空間網(wǎng)格線與時間網(wǎng)格線的交點，即表示了時間空間區(qū)域中一個節(jié)點的位置，相應(yīng)的記為：。2、非穩(wěn)態(tài)項的離散非穩(wěn)態(tài)項的離散有三種不同的格式：1）向前差分2）向后差分3）中心差分1）向前差分將函數(shù)t在節(jié)點（n,i+1）對點（n,i）作泰勒展開，則有： 其中0（）截斷誤差表示余項中的最低階為一次。由上式得：函數(shù)t在節(jié)點（n,i+1）對點（n,i）處一階導(dǎo)數(shù)的向前差分公式： （4-10）2）向后差分 將函數(shù)t在節(jié)點（n,i-1）對點(n,i)作泰勒展開，可得的向后差分公式： （4-11）3）中心差分 的向前差分與向后差分之和，即得 的中心差

15、分表達式： （4-12）二、一維非穩(wěn)態(tài)導(dǎo)熱微分方程的離散方法1、泰勒級數(shù)展開法1）一維非穩(wěn)態(tài)導(dǎo)熱微分方程中的擴散項離散與穩(wěn)態(tài)導(dǎo)熱微分方程中的方法相同，則對一維非穩(wěn)態(tài)導(dǎo)熱微分方程中的擴散項中心差分； 非穩(wěn)態(tài)項向前差分。（1）非穩(wěn)態(tài)項： 采用向前差分為： （4-13）（2）穩(wěn)態(tài)項：采用中心差分則為： （4-14）由此可得：變形得： （4-15）由此可見，只要i時層上各節(jié)點的溫度已知，那么i+1時層上各節(jié)點的溫度即可算出，且不需設(shè)立方程組求解。此關(guān)系式即為顯式差分格式。2）顯示差分與隱式差分格式求解非穩(wěn)態(tài)導(dǎo)熱微分方程，是從已知的初始溫度分布出發(fā)，根據(jù)邊界條件依次求得以后各個時間層上的溫度值。 顯示差

16、分格式定義：就是指若已知i時層上各節(jié)點的溫度值，根據(jù)該差分格式即可算出（i+1）時層上各內(nèi)點的溫度，而不必求解聯(lián)立方程。即是前一時刻（i）n節(jié)點及相鄰兩節(jié)點溫度的顯函數(shù)。優(yōu)點：計算工作量?。蝗秉c：受時間及空間步長的限制。 隱式差分格式對一維非穩(wěn)態(tài)導(dǎo)熱微分方程中的擴散項在(i+1)時層上采用中心差分，非穩(wěn)態(tài)項將t在節(jié)點（n,i+1）處對節(jié)點（n,i）采用向前差分，得： (4-16)式中，已知的是i時層上的值，而未知量有3個，無法求解。定義：就是指已知i 時層上各節(jié)點的溫度值，根據(jù)差分格式不能直接算出(i+1)時層上各節(jié)點的溫度，而必須求解（i+1）時層上的一個聯(lián)立方程組，才能算出(i+1)時層各

17、節(jié)點的溫度，此種差分格式稱隱式差分格式。優(yōu)點：不受時間及空間的步長影響；缺點：計算工作量大。綜上可知： 非穩(wěn)態(tài)導(dǎo)熱微分方程中，擴散項采用中心差分，非穩(wěn)態(tài)項采用向前差分得到顯式差分格式；非穩(wěn)態(tài)導(dǎo)熱微分方程中，擴散項采用中心差分，非穩(wěn)態(tài)項采用向后差分得到隱式差分格式。2、熱平衡法1）優(yōu)點：（1）不受網(wǎng)格是否均勻限制； （2）不受物體是否為常數(shù)限制。2）求解方法如圖4-9所示，一無限大平板，右側(cè)面受周圍流體的冷卻，表面?zhèn)鳠釋?dǎo)數(shù)為h ，對于邊界節(jié)點N代表了寬為的元體。對于該元體，根據(jù)傅立葉定律和能量守恒定律得： =變形為 （4-17）其中是以為特征長度的傅里葉數(shù)，稱網(wǎng)絡(luò)傅里葉數(shù)，記為：。是以為特征長度

18、的畢渥數(shù)，稱網(wǎng)絡(luò)畢渥數(shù)，記為：。一項變形如下：所以 =（22）+2+2 （4-18）補充：4-17的推導(dǎo)過程對于的元體： 根據(jù)傅立葉定律，在i 時層上，從節(jié)點N1傳導(dǎo)給節(jié)點N的熱流量，即從N1傳給元體 單位面積的熱流量為： （a） 根據(jù)牛頓冷卻公式，平板右側(cè)被冷卻時，在i時層上其單位面積損失的熱流量為： （b） 在i時層上元體 熱力學(xué)能的增量： 其中 根據(jù)能量守恒定律可知：在I時層通過導(dǎo)熱和對流進入元體 的能量應(yīng)等于元體熱力學(xué)能的變化量，即=變形為：=（22）+2+2 （4-18）說明：對多維非穩(wěn)態(tài)導(dǎo)熱問題應(yīng)用熱平衡法來建立離散方程的原則與過程與之類似。三、討論一維導(dǎo)熱問題顯式差分格式穩(wěn)定性限制的物理意義從離散方程的結(jié)構(gòu)分析：1、 對于一維導(dǎo)熱顯式格式的內(nèi)節(jié)點方程， 其中由方程式得知，點n上i+1時刻的溫度是在該點i時刻溫度的基礎(chǔ)上計及了左右兩鄰點溫度的影響后得出的。若兩鄰點的影響保持不變，則合理的情況是：越高，則越高；越低，則越低。在上式中，滿足這種合理性是有條件的，即上式中 前的系數(shù)必大于等于零，即 (