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

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

3、三類(lèi)邊界條件，分別由傅立葉定律可描述如下:左邊界：X 0 h2仃X 0 Tf2)( 1c)X右邊界：xhi(T<Tfi)(1d)X3數(shù)值處理與程序設(shè)計(jì)3.1數(shù)值處理2所示。一共采用外點(diǎn)法用均勻網(wǎng)格對(duì)求解區(qū)域進(jìn)行離散化，得到的網(wǎng)格系統(tǒng)如圖 使用了 0N-1共N個(gè)節(jié)點(diǎn)。節(jié)點(diǎn)間距Sx為:?Sx =污L & J< >JE>圖2墻壁的網(wǎng)格劃分此例中墻壁導(dǎo)熱系數(shù)為常值，無(wú)源項(xiàng)。則可采用有限體積法對(duì)控制方程離散化，得到離散方程為:apTpaETEawTwb(2a)式中:aP aEaW0aP(2b)aEawX，ap(2c)(2d)其中的上標(biāo)“ 0”表示此為上一時(shí)刻的值,分別為節(jié)

4、點(diǎn)所在控制容積左右邊界上的導(dǎo)熱系數(shù)，由于墻壁導(dǎo)熱系數(shù)不變，故都等于入，At為時(shí)間步長(zhǎng)。由元體能量平衡法0 0bapTp可以得知左右邊界節(jié)點(diǎn)的離散方程分別為:左邊界節(jié)點(diǎn):(?2?+ ?2 + 爲(wèi))?0 =爲(wèi)？1 + ?2?2 + ?2?0(3) 右邊界節(jié)點(diǎn)：(-2?+ ?1 + ?) ?- 1 = ?- 2 + ?1?1 + ¥?- 1(4) 離散方程的詳細(xì)推導(dǎo)過(guò)程見(jiàn)附錄。3.2程序設(shè)計(jì)由物理模型可以知道本問(wèn)題為一維導(dǎo)熱問(wèn)題，一維導(dǎo)熱問(wèn)題的離散方程在取遍所有節(jié)點(diǎn)后形成的是三對(duì)角的代數(shù)方程組，采用追趕法進(jìn)行求解。程序構(gòu)成和方法：程序由主程序和一個(gè)子程序構(gòu)成。主程序進(jìn)行變量定義和各已知 參

5、數(shù)的輸入，以及左右邊界節(jié)點(diǎn)和部節(jié)點(diǎn)控制方程的輸入；子程序tdma實(shí)現(xiàn)追趕法用來(lái)計(jì)算每個(gè)節(jié)點(diǎn)新的溫度。Thomas算法求解過(guò)程分為兩步：消元和回代。消元是從系數(shù) 矩陣的第二行起，逐一將每一行的非零元素消去一個(gè)，使原來(lái)的三元方程化為二元方程。 消元進(jìn)行到最后一行時(shí)，二元方程就化為一元方程，直接得到最后一個(gè)未知數(shù)的值。然 后逐一往前回代，由各二元方程求出其它未知解。程序特點(diǎn)：該程序有很強(qiáng)的適應(yīng)性，一維常物性非穩(wěn)態(tài)平壁導(dǎo)熱問(wèn)題都可以使用此 程序，只要適當(dāng)更改邊值條件即可。還可以進(jìn)行修改解決非常物性問(wèn)題。程序中對(duì)輸出節(jié)點(diǎn)，最大輸出量都進(jìn)行了控制，對(duì)計(jì)算結(jié)果的分析有很大幫助。而且 Thoms算法的優(yōu)點(diǎn)需要

6、存小，工作量小，程序設(shè)計(jì)簡(jiǎn)單。程序流程圖：首先對(duì)變量賦值，然后由初始條件建立初始溫度場(chǎng)，接著從左邊界， 部節(jié)點(diǎn)，到右邊界進(jìn)行迭代，直到滿(mǎn)足精度要求為止，最后輸出結(jié)果，程序結(jié)束。程序 流程如下圖3。4、模型與程序驗(yàn)證4.1模型本題簡(jiǎn)化為厚度為2 =0.3m的一維非穩(wěn)態(tài)模型如圖4所示，初始溫度為15C，在 其中間建立坐標(biāo)系，左兩邊為對(duì)流換熱，且換熱系數(shù)相同都為h=25 W/(m2C )，且流體溫度Tf=-10 C對(duì)于x 0,列出其導(dǎo)熱微分方程式及定解條件：Q(0x0)(5)圖3程序流程圖aT(x,O) To(0x )cT(x,)x 0x引入過(guò)余溫度:h T( , ) TT(x,)(6)(7)(8)

7、T(x, ) T(9)圖4 一維導(dǎo)熱簡(jiǎn)化模型2 -C)直接根據(jù)公式得到解析解如下:(,)0.Cnexp(n 12nFO)COS( n )(10)式中，F(xiàn)o ay ,，系數(shù)Cn應(yīng)該使上述無(wú)窮級(jí)數(shù)在0是滿(mǎn)足初始條件，由傅里葉級(jí)數(shù)理論可得：2sin n(11)Cn(II丿n COS nSi n nn是超越方程的根，稱(chēng)為特征根。丄Bitan n,n 1,2,n（其中Bi4.2程序驗(yàn)證（1）由模型可以得到相關(guān)信息然后進(jìn)行編程，同等時(shí)間下計(jì)算出中心處溫度的解 析解和數(shù)值解進(jìn)行比較，數(shù)據(jù)記錄在表 1。然后計(jì)算出相對(duì)誤差，作圖5,觀察數(shù)值解 與分析解的比較曲線。由圖表中可以發(fā)現(xiàn)，平壁中心不同時(shí)刻溫度值的分析解

8、和數(shù)值解相差不是很大，二 者吻合的比較好，可以說(shuō)明所編制的數(shù)值解法的程序是正確的。 相對(duì)誤差先增大后減小， 增大的原因是此時(shí)溫度接近零度，相對(duì)誤差的基數(shù)比較小，所以造成相對(duì)誤差較大，但 是此時(shí)的絕對(duì)誤差并不大，在合理圍，所以除去個(gè)別點(diǎn)外，都滿(mǎn)足誤差小于百分之1。可以驗(yàn)證所編數(shù)值解法的程序是正確的（2）空間步長(zhǎng)對(duì)墻壁的溫度影響如圖 6及表2。在程序編寫(xiě)過(guò)程中用網(wǎng)格節(jié)點(diǎn)數(shù) 對(duì)空間步長(zhǎng)進(jìn)行控制，為了觀察空間步長(zhǎng)對(duì)墻壁溫度的影響，表中選擇了三個(gè)不同的空 間步長(zhǎng)，分別為選取51，101, 201個(gè)網(wǎng)格節(jié)點(diǎn)，貝U相應(yīng)的空間步長(zhǎng)為 0.006，0.003， 0.0015。根據(jù)不同步長(zhǎng)時(shí)溫度的變化曲線可以看出

9、，空間步長(zhǎng)對(duì)墻壁的影響不大，當(dāng) 空間步長(zhǎng)控制在合理圍時(shí)可以忽略空間步長(zhǎng)的影響。表1分析解與數(shù)值解比較時(shí)間（h）分析解（C）數(shù)值解（C）相對(duì)誤差（%0151500.55614.85414.811-0.289481.11113.46313.425-0.282261.66711.30611.3-0.053072.2654532.7786.9746.9990.3584743.3335.0825.1120.5903193.8893.3933.4250.9431184.4441.891.9231.7460325. 000.5540.5886.1371845.556-0.631-0.598

10、-5.229796.111-1.684-1.651-1.959626.667-2.618-2.586-1.222317.222-3.447-3.417-0.870327.778-4.184-4.154-0.717028.333-4.837-4.809-0.578878.889-5.417-5.391-0.479979.444-5.932-5.907-0.42144圖5平壁中心不同時(shí)刻的數(shù)值解和分析解表2空間步長(zhǎng)對(duì)溫度影響數(shù)據(jù)時(shí)間（h）N=51N=101N=2010.00015.00015.00015.0000.50015.00015.00015.0001.00014.99814.99814.9

11、981.50014.98014.98114.9812.00014.92414.92414.9242.50014.82014.82014.8203.00014.67514.67514.6763.50014.50114.50214.5024.00014.31014.31114.3114.50014.11114.11214.1125.00013.91113.91113.9125.50013.71513.71513.7156.00013.52513.52513.5256.50013.34313.34413.3447.00013.17113.17213.1727.50013.00900

12、12.85812.85812.8588.50012.71612.71612.7169.00012.58312.58312.5839.50012.46012.46012.4609.91712.36312.36312.36415.014.5)°C14.0(度13.5溫13.012.512.0時(shí)間（h）圖6空間步長(zhǎng)對(duì)墻壁溫度影響（3）時(shí)間步長(zhǎng)對(duì)溫度的影響如圖7和表3,根據(jù)圖中曲線可以看出時(shí)間步長(zhǎng)選擇 50s, 100s, 200s時(shí)基本重合，對(duì)墻壁溫度影響不大55 4 5 3 54 13 12T T TT=50 （s）T=100sT=200s246時(shí)間（h）8 10圖7時(shí)間步長(zhǎng)對(duì)溫度的影響

13、時(shí)間（h）T=50 （s）時(shí)間（h）T=100（ s）時(shí)間（h）T=200 （s）015015.0000150.694150.69415.0000.667151.38914.9881.38914.9871.33314.9872.14.9122.14.910214.9182.77814.7462.77814.7442.66714.773.47214.5133.47214.5113.33314.5584.16714.2454.16714.244414.3084.86113.9664.86113.9664.66714.0445.55613.6925.55613.6935.33313.7816.251

14、3.4316.2513.433613.5276.94413.1886.94413.1906.66713.2887.63912.9647.63912.9667.33313.8.33312.768.33312.762812.8629.12.5749.12.5768.66712.675表3時(shí)間步長(zhǎng)對(duì)溫度的影響16100246810121416182022時(shí)間（h）圖8墻壁溫度隨時(shí)間的變化曲線5計(jì)算結(jié)果與分析5.1墻壁溫度分析根據(jù)題目中要求，計(jì)算寒潮入侵多長(zhǎng)時(shí)間后墻壁可以感受到外界氣溫的變化，通過(guò) 建模，方程離散化，最終通過(guò)程序求解方程，得到圖8和表4。由圖可以看出，開(kāi)始階段，墻壁溫不變，隨著時(shí)間的進(jìn)

15、一步深入，壁溫度開(kāi)始降低，當(dāng)很長(zhǎng)時(shí)間后，溫度變化 基本趨于平緩，直到再次平衡。根據(jù)圖 8就可以得到墻壁溫度開(kāi)始發(fā)生變化的時(shí)間。表4墻壁溫度隨時(shí)間變化數(shù)據(jù)表時(shí)間（h）溫度C）時(shí)間（h）溫度（C）時(shí)間（h）溫度（C）0.000156.66713.28513.33311.7630.556157.22213.09813.88911.6911.11114.9977.77812.92414.44411.6261.66714.9678.33312.76215.00011.5652.22214.8838.88912.61215.55611.512.77814.7449.44412.47316.11111.45

16、83.33314.56210.00012.34516.66711.4113.88914.35410.55612.22717.22211.3684.44414.13311.11112.11817.77811.3295.00013.91111.66712.01818.33311.2925.55613.69312.22211.92618.88911.2596.11113.48412.77811.84119.44411.2285.2導(dǎo)熱系數(shù)入對(duì)墻壁溫度的影響墻的導(dǎo)熱系數(shù)對(duì)表面的影響，在圖9和表5中發(fā)現(xiàn)，導(dǎo)熱系數(shù)對(duì)壁溫度影響比較大，入=1.2時(shí)，溫度下降的趨勢(shì)會(huì)更快，要比入=0.85時(shí)下降快的多，下降速

17、度更快，更短時(shí)間達(dá)到穩(wěn)態(tài)，因此導(dǎo)熱系數(shù)越大溫度擴(kuò)散越快，導(dǎo)熱系數(shù)越小則溫度變化越慢，需要 更長(zhǎng)時(shí)間到達(dá)穩(wěn)態(tài)，但是這時(shí)對(duì)于要求恒溫的空間有好處，受波動(dòng)影響更小。表5導(dǎo)熱系數(shù)對(duì)墻壁溫度的影響時(shí)刻（h）入=0.85入=1.0入=1.201515150.5561515151.11114.99714.9914.9711.66714.96714.92414.8332.22214.88314.7714.5642.77814.74414.54214.2073.33314.56214.26613.8083.88914.35413.96813.3984.44414.13313.66612.999513.91113

18、.37112.62圖9導(dǎo)熱系數(shù)對(duì)墻壁溫度的影響5.55613.69313.08812.2676.11113.48412.82311.9426.66713.28512.57511.6447.22213.09812.34611.3727.77812.92412.13611.1258.33312.76211.94210.9018.88912.61211.76410.6989.44412.47311.60210.5145.3墻外換熱系數(shù)h的影響墻外表面?zhèn)鳠嵯禂?shù)對(duì)溫度分布的影響，如圖10和表6影響不大。16151413121110024681012141618時(shí)間（h）圖10墻外換熱系數(shù)對(duì)溫度影響表6墻

19、外換熱系數(shù)對(duì)溫度影響時(shí)刻 h=20h=35h=50（h） W/（mfC ） W/（mC） W/（mC)01515151.11114.99814.99714.9962.22214.91414.88314.8653.33314.66214.56214.5074.44414.31114.13314.0435.55613.93713.69313.5736.66713.5813.28513.1447.77813.25712.92412.7688.88912.97212.61212.4461012.72412.34512.17311.11112.5112.11811.94212.22212.32511.9

20、2611.74813.33312.16711.76311.58514.44412.11.62611.44815.55611.91511.5111.33416.66711.81511.41111.23717.77811.7311.32911.15618.88911.65711.25911.0885.4墻壁傳熱熱系數(shù)的影響由圖11可以看出，墻表面?zhèn)鳠嵯禂?shù)對(duì)表面溫度影響較大，當(dāng)傳熱系數(shù)比較小時(shí)受 墻外流體溫度影響明顯，傳熱系數(shù)越大受墻外流體溫度影響越小，當(dāng)?shù)搅藰O限大的情況 下，表面溫度則等于墻流體溫度，不再受墻外流體溫度影響。402468101214161820時(shí)刻（h）圖11墻表面?zhèn)鳠嵯禂?shù)對(duì)溫度影

21、響表7墻表面?zhèn)鳠嵯禂?shù)對(duì)溫度影響時(shí)刻h=2h=6h=10（h） W/（mC ） W/（mfC ） W/（mfC ）0.00015.00015.00015.0001.11114.99414.99714.9992.22214.78914.88314.9553.33314.16914.56214.8374.44413.29214.13314.6885.55612.33213.69314.5406.66711.39013.28514.4087.77810.50712.92414.2958.8899.70012.61214.20110.0008.97212.34514.12311.1118.31812.1

22、1814.12.2227.73511.92614.00713.3337.21411.76313.96414.4446.74911.62613.92915.5566.33611.51013.90016.6675.96811.41113.87617.7785.64011.32913.85718.8895.34811.25913.8415.5墻壁厚度S對(duì)壁溫度的影響改變墻壁的厚度3,壁溫度將發(fā)生變化。在表8和圖12中可以發(fā)現(xiàn)，不同厚度的 墻壁對(duì)外界溫度的感應(yīng)快慢是不一樣的，隨著墻壁厚度的增加，感應(yīng)到外界溫度變化的 時(shí)間越來(lái)越長(zhǎng)，且溫度變化越來(lái)越慢，墻壁越厚要越長(zhǎng)的時(shí)間才能達(dá)到新的平衡時(shí)間（h）圖12

23、 墻壁厚度對(duì)壁溫度的影響表8墻壁厚度對(duì)壁溫度的影響時(shí)間（h）0.15m0.3m0.45m01515151.1112.22211.35414.88314.9993.3339.56814.56214.9894.4448.43614.13314.9555.5567.72713.69314.8946.6677.28513.28514.817.7787.00912.92414.7138.8896.83612.61214.608106.72812.34514.50111.1116.66112.11814.39712.2226.61911.92614.29613.3336.59211.76314.214.4

24、446.57611.62614.1115.5566.56611.5114.02616.6676.55911.41113.94917.7786.55511.32913.87718.8896.55311.25913.8116結(jié)論(1) 整個(gè)墻壁經(jīng)歷了以下變化過(guò)程：首先外壁直接與室外冷空氣接觸，溫度變化很 快，隨著時(shí)間的推移，墻各點(diǎn)的溫度也開(kāi)始變化，并影響到右邊墻壁的溫度，慢慢降低(2) 對(duì)于墻壁的溫度隨時(shí)間的變化，變化趨勢(shì)總是由快到慢，最后重新達(dá)到穩(wěn)態(tài)。 當(dāng)改變墻壁厚度的時(shí)候，隨著墻壁厚度的增加，對(duì)于外界溫度的感應(yīng)是越來(lái)越慢。這對(duì) 于一些對(duì)溫度有要求的地反很重要，根據(jù)溫度的變化作出適當(dāng)?shù)恼{(diào)整措施。

25、參考文獻(xiàn)1 黃善波中良編著計(jì)算傳熱學(xué)基礎(chǔ)2 世銘 陶文銓編著 傳熱學(xué)(第四版)北京高等教育3 王元明編數(shù)學(xué)物理方程與特殊函數(shù)(第三版)高等教育附錄2 P 彳PP 彳p_ 丄 _ 丄e21數(shù)學(xué)模型的離散過(guò)程推導(dǎo)11P1 P WPXXX按照Taylor級(jí)數(shù)展開(kāi)法，溫度對(duì)時(shí)間的偏導(dǎo)有向刖差分格式，中心差分格式和向后差分格式，使用向后差分格式向后差分格式：pp Pp 1Pp Pp 1PaPp E11pp WXXeXeXwPXwXpp 1p E11pp W3)得XX由aXP2pXeXecXXw0PXPXEXWPaPaEEa3WWb3eaw03eawap0 c X aP5)ap邊界條件處理:x右邊界根據(jù)元

26、體能量平衡法處理:b aPTP0(7)如圖1:其中LI2 /Ob <1 1N-2N-1:><-Ax圖1右邊界的能量守恒EsB h(Tf T n1),hr(TfrTN 1)(8)(9)為常數(shù)將式（12-b） （12- c）入式（12-a）以得到:c x hhrN 12xN 2 hrTfrc xToIn 12同理可以得到左邊界點(diǎn)離散方程：(hihTfixTo0因此，離散方程為：apTpaETEaWTW式中,0ap aE aw apaE，awxo，aPx0 0apTp(12)(13)(14)(15)(16)(17)式中，上標(biāo)“ 0”表示上一時(shí)刻值，為熱擴(kuò)散系數(shù)，為時(shí)間步長(zhǎng)由元體能量

27、平衡法可以分別求出左右邊界節(jié)點(diǎn)的離散方程，左邊界為：00h2 )T。T1h2Tf2ap 丁 0c10(18)2xx2右邊界為：00ap(寸h1)T N 1-Tn 2h1Tf1ap 丁0 1 N 1(19)2xx22程序清單序號(hào)程序名稱(chēng)程序功能計(jì)算結(jié)果(1)設(shè)計(jì)程序用數(shù)值方法求解 溫度 表1、4,圖5、8 改變空間步長(zhǎng)：表2,圖6 改變時(shí)間步長(zhǎng)：表3,圖7 改變導(dǎo)熱系數(shù)：表5,圖9 改變墻外換熱系數(shù)：表6,圖10 改變墻換熱系數(shù)：表7,圖11 改變墻壁厚度：表8,圖12（2）驗(yàn)證程序用分析解求解溫 度表1,圖52.1設(shè)計(jì)程序主要程序段(1) 初始時(shí)刻的輸入dx=dlt/(N-1); 空間步長(zhǎng)apO=dc*dx/dtao; 系數(shù) apOfor(i=0;i<N;i+)TOi=Tr-h1*(Tf1-Tr)*(dlt-i*dx)/lw;初始溫度Twi=Tr-h1*(Tf1-Tr)*(dlt-i*dx)/lw;Tw0i=Tr-h1*(Tf1-Tr)*(dlt-i*dx)/lw;/假設(shè)溫度場(chǎng)(2) 邊界條件及中間節(jié)點(diǎn)程序輸入/左邊界節(jié)點(diǎn)(墻外)C0=lw/dx;B0