計(jì)算傳熱學(xué)程序設(shè)計(jì)

1、中國石油大學(xué) （華東 ）儲(chǔ)運(yùn)與建筑工程學(xué)院熱能與動(dòng)力工程系計(jì)算傳熱學(xué)程序設(shè)計(jì)設(shè)計(jì)報(bào)告學(xué)生姓名 : 學(xué) 號(hào)： 專業(yè)班級(jí)：指導(dǎo)教師2012 年 7 月 7 日1、設(shè)計(jì)題目有一房屋的磚墻厚 S = m ,入=W/（mC）, pc=xiO6J/（ m 3 - K），室內(nèi)溫度Tfi 保持20C不變，表面?zhèn)鳠嵯禂?shù)hi=6W/（mC ）。開始時(shí)墻的溫度處于穩(wěn)定狀態(tài)，內(nèi)墻表 面溫度Tw為15C寒潮入侵后，室外溫度Tf2下降為-10 C,外墻的表面?zhèn)鳠嵯禂?shù)為 35W（m2C ）。試分析寒潮入侵后多少時(shí)間內(nèi)墻壁面方可感受到外界氣溫的變化。已知參數(shù)壁厚，墻壁導(dǎo)熱系數(shù)，密度與比熱容的乘積，室內(nèi)和寒潮入侵后室外空氣溫度

2、，室 內(nèi)空氣和外墻的表面?zhèn)鳠嵯禂?shù)，開始時(shí)穩(wěn)定狀態(tài)下的內(nèi)墻表面溫度。求解寒潮入侵多少時(shí)間后內(nèi)墻壁面可感受到外界氣溫的變化2物理與數(shù)學(xué)模型物理模型該墻面為常物性，可以假設(shè)：（1）其為無限大平面，（2）只有在厚度方向傳熱, 沒有縱向傳熱，則該問題轉(zhuǎn)化為一維常物性無限大平面非穩(wěn)態(tài)導(dǎo)熱問題。數(shù)學(xué)模型以墻外表面為坐標(biāo)原點(diǎn)，沿厚度方向?yàn)樽鴺?biāo)正方向，建立坐標(biāo)系?；谏鲜瞿Ｐ停?取其在x方向上的微元作為研究對(duì)象，則該問題的數(shù)學(xué)模型可描述如下：(1a)T)初始條件:T | % - q = T.i - hj (Tf i - T剛)：6 - X)/ X(1b)在兩側(cè)相應(yīng)的邊界條件是第三類邊界條件，分別由傅立葉定律可描

3、述如下:左邊界：-1x0h2(Tx 0 Tf2)(1c)X右邊界：|xhi(TxTfi)(1d)X3數(shù)值處理與程序設(shè)計(jì)數(shù)值處理采用外點(diǎn)法用均勻網(wǎng)格對(duì)求解區(qū)域進(jìn)行離散化，得到的網(wǎng)格系統(tǒng)如圖2所示。一共使用了 0N-1共N個(gè)節(jié)點(diǎn)。節(jié)點(diǎn)間距Sx為：%二占圖2墻壁內(nèi)的網(wǎng)格劃分此例中墻壁導(dǎo)熱系數(shù)為常值，無源項(xiàng)。則可采用有限體積法對(duì)控制方程離散化，得到離散方程為:apTpaETEawTwb(2a)式中:aP aEaW0aP(2b)(2c)0C X-，aPX(2d)0 0bapTp其中的上標(biāo)“ 0”表示此為上一時(shí)刻的值,分別為節(jié)點(diǎn)所在控制容積左右邊界上的導(dǎo)熱系數(shù)，由于墻壁導(dǎo)熱系數(shù)不變，故都等于入，t為時(shí)間步

4、長。由元體能量平衡法可以得知左右邊界節(jié)點(diǎn)的離散方程分別為: 左邊界節(jié)點(diǎn):(3)(I + +右邊界節(jié)點(diǎn):(2 + h丄 + 臥)Tn - 1 二阪 - 2 + hf丄 4 製n - 1(4)離散方程的詳細(xì)推導(dǎo)過程見附錄。程序設(shè)計(jì)由物理模型可以知道本問題為一維導(dǎo)熱問題，一維導(dǎo)熱問題的離散方程在取遍所有節(jié)點(diǎn)后形成的是三對(duì)角的代數(shù)方程組，采用追趕法進(jìn)行求解。程序構(gòu)成和方法：程序由主程序和一個(gè)子程序構(gòu)成。主程序進(jìn)行變量定義和各已知 參數(shù)的輸入，以及左右邊界節(jié)點(diǎn)和內(nèi)部節(jié)點(diǎn)控制方程的輸入；子程序 tdma實(shí)現(xiàn)追趕法 用來計(jì)算每個(gè)節(jié)點(diǎn)新的溫度。Thomas算法求解過程分為兩步：消元和回代。消元是從系 數(shù)矩陣的

5、第二行起，逐一將每一行的非零元素消去一個(gè)，使原來的三元方程化為二元方 程。消元進(jìn)行到最后一行時(shí)，二元方程就化為一元方程，直接得到最后一個(gè)未知數(shù)的值。 然后逐一往前回代，由各二元方程求出其它未知解。程序特點(diǎn)：該程序有很強(qiáng)的適應(yīng)性，一維常物性非穩(wěn)態(tài)平壁導(dǎo)熱問題都可以使用此 程序，只要適當(dāng)更改邊值條件即可。還可以進(jìn)行修改解決非常物性問題。程序中對(duì)輸出節(jié)點(diǎn)，最大輸出量都進(jìn)行了控制，對(duì)計(jì)算結(jié)果的分析有很大幫助。而且 Thoms算法的優(yōu)點(diǎn)需要內(nèi)存小，工作量小，程序設(shè)計(jì)簡單。程序流程圖：首先對(duì)變量賦值，然后由初始條件建立初始溫度場，接著從左邊界， 內(nèi)部節(jié)點(diǎn)，到右邊界進(jìn)行迭代，直到滿足精度要求為止，最后輸出結(jié)

6、果，程序結(jié)束。程 序流程如下圖3。4、模型與程序驗(yàn)證模型本題簡化為厚度為2 =的一維非穩(wěn)態(tài)模型如圖4所示，初始溫度為15C，在其中 間建立坐標(biāo)系，左兩邊為對(duì)流換熱，且換熱系數(shù)相同都為h=25 W/(m2C )，且流體溫度Tf=-10 C對(duì)于x 0,列出其導(dǎo)熱微分方程式及定解條件：2a -2 (0 x ,0)(5)x圖3程序流程圖a cT(x,O) T(0(6)T(x,)(7)引入過余溫度:h T(,)T(x,)(8)1i、1i左側(cè)：h TfIiiii iiIi iiIi-0右側(cè)：2h=25 W/(m C)Tf=-10C廠x圖4 一維導(dǎo)熱簡化模型To=15C直接根據(jù)公式得到解析解如下:(,)0.

7、Cnexp(n 12nFO)COS( n )(10)式中，F(xiàn)O ay,-，系數(shù)Cn應(yīng)該使上述無窮級(jí)數(shù)在0是滿足初始條件，由傅里葉級(jí)數(shù)理論可得：2sin n(11)nn COSnSin nn是超越方程的根，稱為特征根。tan n Bi,nn1,2,（其中Bi程序驗(yàn)證（1）由模型可以得到相關(guān)信息然后進(jìn)行編程，同等時(shí)間下計(jì)算出中心處溫度的解 析解和數(shù)值解進(jìn)行比較，數(shù)據(jù)記錄在表 1。然后計(jì)算出相對(duì)誤差，作圖5,觀察數(shù)值解 與分析解的比較曲線。由圖表中可以發(fā)現(xiàn)，平壁中心不同時(shí)刻溫度值的分析解和數(shù)值解相差不是很大，二 者吻合的比較好，可以說明所編制的數(shù)值解法的程序是正確的。 相對(duì)誤差先增大后減小， 增大的

8、原因是此時(shí)溫度接近零度，相對(duì)誤差的基數(shù)比較小，所以造成相對(duì)誤差較大，但 是此時(shí)的絕對(duì)誤差并不大，在合理范圍內(nèi)，所以除去個(gè)別點(diǎn)外，都滿足誤差小于百分之 1?？梢则?yàn)證所編數(shù)值解法的程序是正確的。（2）空間步長對(duì)墻內(nèi)壁的溫度影響如圖 6及表2。在程序編寫過程中用網(wǎng)格節(jié)點(diǎn)根據(jù)數(shù)對(duì)空間步長進(jìn)行控制，為了觀察空間步長對(duì)墻內(nèi)壁溫度的影響，表中選擇了三個(gè)不同 的空間步長，分別為選取51，101，201個(gè)網(wǎng)格節(jié)點(diǎn)，則相應(yīng)的空間步長為， 。 不同步長時(shí)溫度的變化曲線可以看出，空間步長對(duì)內(nèi)墻壁的影響不大，當(dāng)空間步長控制 在合理范圍時(shí)可以忽略空間步長的影響。表1分析解與數(shù)值解比較時(shí)間（h）分析解（C）數(shù)值解（C）相對(duì)

9、誤差（0151505. 00圖6空間步長對(duì)墻內(nèi)壁溫度影響表2空間步長對(duì)溫度影響數(shù)據(jù)時(shí)間（h）N=51N=101N=201圖5平壁中心不同時(shí)刻的數(shù)值解和分析解（3）時(shí)間步長對(duì)溫度的影響如圖7和表3,根據(jù)圖中曲線可以看出時(shí)間步長選擇 50s, 100s, 200s時(shí)基本重合，對(duì)墻內(nèi)壁溫度影響不大。圖7時(shí)間步長對(duì)溫度的影響表3時(shí)間步長對(duì)溫度的影響時(shí)間（h）T=50 （s） 時(shí)間（h）T=100（ s） 時(shí)間（h）T=200（ s）015001515152468圖8墻內(nèi)壁溫度隨時(shí)間的變化曲線5計(jì)算結(jié)果與分析墻內(nèi)壁溫度分析根據(jù)題目中要求，計(jì)算寒潮入侵多長時(shí)間后內(nèi)墻壁可以感受到外界氣溫的變化，通 過建模，

10、方程離散化，最終通過程序求解方程，得到圖 8和表4。由圖可以看出，開始 階段，內(nèi)墻壁溫不變，隨著時(shí)間的進(jìn)一步深入，內(nèi)壁溫度開始降低，當(dāng)很長時(shí)間后，溫度變化基本趨于平緩，直到再次平衡。根據(jù)圖 8就可以得到墻內(nèi)壁溫度開始發(fā)生變化的 時(shí)間。表4墻內(nèi)壁溫度隨時(shí)間變化數(shù)據(jù)表時(shí)間（h）溫度C）時(shí)間（h）溫度C）時(shí)間（h） 溫度C）1515導(dǎo)熱系數(shù)入對(duì)墻內(nèi)壁溫度的影響墻的導(dǎo)熱系數(shù)對(duì)內(nèi)表面的影響，在圖 9和表5中發(fā)現(xiàn)，導(dǎo)熱系數(shù)對(duì)內(nèi)壁溫度影響比 較大，入=時(shí)，溫度下降的趨勢會(huì)更快，要比 入=時(shí)下降快的多，下降速度更快，更短時(shí) 間內(nèi)達(dá)到穩(wěn)態(tài)，因此導(dǎo)熱系數(shù)越大溫度擴(kuò)散越快，導(dǎo)熱系數(shù)越小則溫度變化越慢，需要 更長時(shí)間

11、到達(dá)穩(wěn)態(tài)，但是這時(shí)對(duì)于要求恒溫的空間有好處，受波動(dòng)影響更小。表5導(dǎo)熱系數(shù)對(duì)墻內(nèi)壁溫度的影響時(shí)刻（h）入二入二入二0151515151515圖9導(dǎo)熱系數(shù)對(duì)墻內(nèi)壁溫度的影響墻外換熱系數(shù)h的影響墻外表面?zhèn)鳠嵯禂?shù)對(duì)溫度分布的影響，如圖 10和表6影響不大。圖10墻外換熱系數(shù)對(duì)溫度影響表6墻外換熱系數(shù)對(duì)溫度影響h=502時(shí)刻 h=20h=35W/(mC(h) w/(mC ) w/(mC )015151510墻內(nèi)壁傳熱熱系數(shù)的影響由圖11可以看出，墻內(nèi)表面?zhèn)鳠嵯禂?shù)對(duì)內(nèi)表面溫度影響較大，當(dāng)傳熱系數(shù)比較小 時(shí)受墻外流體溫度影響明顯，傳熱系數(shù)越大受墻外流體溫度影響越小，當(dāng)?shù)搅藰O限大的 情況下，內(nèi)表面溫度則等于墻

12、內(nèi)流體溫度，不再受墻外流體溫度影響。圖11墻內(nèi)表面?zhèn)鳠嵯禂?shù)對(duì)溫度影響表7墻內(nèi)表面?zhèn)鳠嵯禂?shù)對(duì)溫度影響時(shí)刻h=2h=6h=102 2 2(h) W/(m ) W/(m ) W/(m )墻壁厚度S對(duì)內(nèi)壁溫度的影響改變墻壁的厚度3,內(nèi)壁溫度將發(fā)生變化。在表8和圖12中可以發(fā)現(xiàn)，不同厚度 的墻壁對(duì)外界溫度的感應(yīng)快慢是不一樣的，隨著墻壁厚度的增加，感應(yīng)到外界溫度變化 的時(shí)間越來越長，且溫度變化越來越慢，墻壁越厚要越長的時(shí)間才能達(dá)到新的平衡圖12 墻壁厚度對(duì)內(nèi)壁溫度的影響表8墻壁厚度對(duì)內(nèi)壁溫度的影響時(shí)間（h）0151515156結(jié)論(1) 整個(gè)墻壁經(jīng)歷了以下變化過程：首先外壁直接與室外冷空氣接觸，溫度變化很

13、 快，隨著時(shí)間的推移，墻內(nèi)各點(diǎn)的溫度也開始變化，并影響到右邊內(nèi)墻壁的溫度，慢慢 降低。(2) 對(duì)于墻內(nèi)壁的溫度隨時(shí)間的變化，變化趨勢總是由快到慢，最后重新達(dá)到穩(wěn)態(tài)。 當(dāng)改變墻壁厚度的時(shí)候，隨著墻壁厚度的增加，對(duì)于外界溫度的感應(yīng)是越來越慢。這對(duì) 于一些對(duì)溫度有要求的地反很重要，根據(jù)溫度的變化作出適當(dāng)?shù)恼{(diào)整措施。參考文獻(xiàn)1 黃善波 劉中良編著 計(jì)算傳熱學(xué)基礎(chǔ)2 楊世銘陶文銓編著傳熱學(xué)(第四版)北京高等教育出版社3 王元明編數(shù)學(xué)物理方程與特殊函數(shù)(第三版)高等教育出版社附錄1數(shù)學(xué)模型的離散過程推導(dǎo)w1 Xe1 X1按照Taylor級(jí)數(shù)展開法，溫度對(duì)時(shí)間的偏導(dǎo)有向刖差分格式，中心差分格式和向 后差分格

14、式，使用向后差分格式。向后差分格式：3)W PW XePEXaxxpp 1aPP由x x得x 2pxp E11p Pp wx ex ex wx wEwc x 0PxxawEw bap邊界條件處理:如圖1:其中為常數(shù)aPaEoaEawapaEaWx右邊界根據(jù)元體能量平衡法處理:0aPb a；TP1Iy2/B-1，x x圖1右邊界的能量守恒BN 2EsBhr (TfrTn1),N 2N 1N2xN 2N 1. z-rhr (TfrTn1)c x0N 1N 1xx2將式（12-b）（12- c）入式（12-a）以得到:c x hhrN2xhrTfr同理可以得到左邊界點(diǎn)離散方程：hTfl-T1x-T0

15、0(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)因此，離散方程為NeTeb式中,(14)(15)aE -，aW，aPxxb a；TP式中，上標(biāo)“ 0”表示上一時(shí)刻值，為熱擴(kuò)散系數(shù)，為時(shí)間步長0ap Ne aw apo c x(17)由元體能量平衡法可以分別求出左右邊界節(jié)點(diǎn)的離散方程, 左邊界為：0ap0aP 0T02(Th2 )T。xT1h2Tf2x右邊界為：00h1)T N 1-Tn 2h1Tf1aPtN2xx2(18)(19)2程序清單序號(hào)程序名稱程序功能計(jì)算結(jié)果(1)設(shè)計(jì)程序用數(shù)值方法求解 溫度 表1、4,圖5、8 改變空間步長：表2,圖6 改變時(shí)間步長：表3,圖7 改變導(dǎo)熱系數(shù)：表5,圖9 改變墻外換熱系數(shù)：表6,圖10 改變墻內(nèi)換熱系數(shù)：表7,圖11 改變墻壁厚度：表8,圖12(2)驗(yàn)證程序用分析解求解溫 度表1,圖5設(shè)計(jì)程序主要程序段(1)初始時(shí)刻的輸入dx=dlt/(N-1);h2+lw/dx;A0=;f0=h2*Tf2+ap0*T00/;CN-1=;fN-1=ap0*T0N-1/2.+h1*Tf1(2)追趕法程序tdma(N,A,B,C,f,Tw);fa=a0*ta n(a0)-Bi;fc=c*ta n(c)-Bi;if(