計算傳熱學程序的設計說明

中國石油大學(華東)儲運與建筑工程學院熱能與動力工程系《計算傳熱學程序設計》設計報告學生:學 號：專業(yè)班級：指導教師2012 年7 月7 日. . .1、設計題目有一房屋的磚墻厚δ=0.3m，λ=0.85W/(m·℃)，ρc6J/(m3·K)，=1.05×10室溫度Tf1保持℃不變，表面?zhèn)鳠嵯禂?shù)h12·℃。開始時墻的溫度處于穩(wěn)定狀20=6W/(m)態(tài)，墻表面溫度Tw為℃寒潮入侵后，室外溫度T下降為-10℃，外墻的表面?zhèn)鳠嵯?15f2數(shù)為35W/(m2·℃)。試分析寒潮入侵后多少時間墻壁面方可感受到外界氣溫的變化。室外室內(nèi)寒流入侵0 x圖1墻壁簡化圖1.1已知參數(shù)壁厚，墻壁導熱系數(shù)，密度與比熱容的乘積，室和寒潮入侵后室外空氣溫度，室空氣和外墻的表面?zhèn)鳠嵯禂?shù)，開始時穩(wěn)定狀態(tài)下的墻表面溫度。1.2 求解寒潮入侵多少時間后墻壁面可感受到外界氣溫的變化？物理與數(shù)學模型2.1 物理模型該墻面為常物性，可以假設：（ 1）其為無限大平面，（ 2）只有在厚度方向傳熱，沒有縱向傳熱，則該問題轉(zhuǎn)化為一維常物性無限大平面非穩(wěn)態(tài)導熱問題。2.2 數(shù)學模型以墻外表面為坐標原點，沿厚度方向為坐標正方向，建立坐標系?；谏鲜瞿Ｐ停∑湓趚方向上的微元作為研究對象，則該問題的數(shù)學模型可描述如下：cT(T)（1a）xx初始條件：??|??=0=????1-??1(????1-????1)(δ-x)/λ（1b）. . .在兩側(cè)相應的邊界條件是第三類邊界條件，分別由傅立葉定律可描述如下：左邊界：T（1c）x0h2(Tx0Tf2)X右邊界：Txh1(TxTf1)（1d）X數(shù)值處理與程序設計3.1 數(shù)值處理采用外點法用均勻網(wǎng)格對求解區(qū)域進行離散化，得到的網(wǎng)格系統(tǒng)如圖 2所示。一共使用了0~N-1共N個節(jié)點。節(jié)點間距δx為:δx=

????-1圖2墻壁的網(wǎng)格劃分此例中墻壁導熱系數(shù)為常值，無源項。則可采用有限體積法對控制方程離散化，得到離散方程為：apTpaETEaWTWb（2a）式中：aPaEaWaP0（2b）aE，aWa0cx（）xxba0pTp0（2d）其中的上標“0”表示此為上一時刻的值,分別為節(jié)點所在控制容積左右邊界上的導熱系數(shù)，由于墻壁導熱系數(shù)不變，故都等于λ，△τ為時間步長。由元體能量平衡法可以得知左右邊界節(jié)點的離散方程分別為：. . .左邊界節(jié)點：0λλ0??????????????0(+??2+)??0=+2????2????（3）右邊界節(jié)點：(????0+??1+λλ????+????0??02)????-1=??????2+??2??1????（4）離散方程的詳細推導過程見附錄。3.2程序設計由物理模型可以知道本問題為一維導熱問題，一維導熱問題的離散方程在取遍所有節(jié)點后形成的是三對角的代數(shù)方程組，采用追趕法進行求解。程序構(gòu)成和方法:程序由主程序和一個子程序構(gòu)成。主程序進行變量定義和各已知參數(shù)的輸入，以及左右邊界節(jié)點和部節(jié)點控制方程的輸入；子程序tdma實現(xiàn)追趕法用來計算每個節(jié)點新的溫度。Thomas算法求解過程分為兩步:消元和回代。消元是從系數(shù)矩陣的第二行起，逐一將每一行的非零元素消去一個，使原來的三元方程化為二元方程。消元進行到最后一行時，二元方程就化為一元方程，直接得到最后一個未知數(shù)的值。然后逐一往前回代，由各二元方程求出其它未知解。程序特點：該程序有很強的適應性，一維常物性非穩(wěn)態(tài)平壁導熱問題都可以使用此程序，只要適當更改邊值條件即可。還可以進行修改解決非常物性問題。程序中對輸出節(jié)點，最大輸出量都進行了控制，對計算結(jié)果的分析有很大幫助。而且Thoms算法的優(yōu)點需要存小，工作量小，程序設計簡單。程序流程圖：首先對變量賦值，然后由初始條件建立初始溫度場，接著從左邊界，部節(jié)點，到右邊界進行迭代，直到滿足精度要求為止，最后輸出結(jié)果，程序結(jié)束。程序流程如下圖3。4、模型與程序驗證4.1模型本題簡化為厚度為 2 =0.3m的一維非穩(wěn)態(tài)模型如圖 4所示，初始溫度為 15℃，在其中間建立坐標系，左兩邊為對流換熱，且換熱系數(shù)相同都為 h=25W/(m2·℃)，且流體溫度Tf=-10℃對于x 0，列出其導熱微分方程式及定解條件：. . .Ta2T(0x,0)（5）x2開始進行參數(shù)賦值

輸入至屏幕用迭代法 輸入至文件求解溫度N程序結(jié)束計數(shù)并判斷是否完成溫度場計算Y圖3程序流程圖aT(x,0)T0(0x)（6）cT(x,)x00（7）xT(x,)hT(,)（8）Txx引入過余溫度：. . .T(x,)T（9）T0=15℃左側(cè)：hTf

右側(cè)：h=25W/(m2·℃)Tf=-10℃-δ 0 δx圖4一維導熱簡化模型直接根據(jù)公式得到解析解如下：(, ) Cnexp( n2Fo)cos( n ) （10）n1式中，F(xiàn)o a2, x，系數(shù)Cn應該使上述無窮級數(shù)在 0是滿足初始條件，由傅里葉級數(shù)理論可得：Cn2sinn（11）cosnsinnn是超越方程的根，稱為特征根。tannBi,n1,2,?（12）nh其中 Bi 。4.2程序驗證（1）由模型可以得到相關信息然后進行編程，同等時間下計算出中心處溫度的解析解和數(shù)值解進行比較，數(shù)據(jù)記錄在表 1。然后計算出相對誤差，作圖 5，觀察數(shù)值解與分析解的比較曲線。由圖表中可以發(fā)現(xiàn)，平壁中心不同時刻溫度值的分析解和數(shù)值解相差不是很大，二者吻合的比較好，可以說明所編制的數(shù)值解法的程序是正確的。 相對誤差先增大后減小，增大的原因是此時溫度接近零度，相對誤差的基數(shù)比較小，所以造成相對誤差較大，但是此時的絕對誤差并不大，在合理圍，所以除去個別點外，都滿足誤差小于百分之 1。. . .可以驗證所編數(shù)值解法的程序是正確的。（2）空間步長對墻壁的溫度影響如圖 6及表2。在程序編寫過程中用網(wǎng)格節(jié)點數(shù)對空間步長進行控制，為了觀察空間步長對墻壁溫度的影響，表中選擇了三個不同的空間步長，分別為選取 51，101，201個網(wǎng)格節(jié)點，則相應的空間步長為 0.006，0.003，0.0015。根據(jù)不同步長時溫度的變化曲線可以看出，空間步長對墻壁的影響不大，當空間步長控制在合理圍時可以忽略空間步長的影響。表1分析解與數(shù)值解比較時間（h）分析解（℃）數(shù)值解（℃）相對誤差（%）0151500.55614.85414.811-0.289481.11113.46313.425-0.282261.66711.30611.3-0.053072.2229.9.0.1654532.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-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.421441816141210）8分析解（℃）℃6（數(shù)值解（℃）度4溫20-2 0 1 2 3 4 5 6 7 8 9 10-4時間（h）-6-8圖5平壁中心不同時刻的數(shù)值解和分析解. . .15.515.0N=5114.5N=101）℃14.0 N=201（度13.5溫13.012.512.00 1 2 3 4 5 6 7 8 9 10 11時間（h）圖6空間步長對墻壁溫度影響表2空間步長對溫度影響數(shù)據(jù)時間(h)N=51N=101N=2010.00015.00015.00015.0000.50015.00015.00015.0001.00014.99814.99814.9981.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.00913.13.8.00012.85812.85812.8588.50012.71612.71612.7169.00012.58312.58312.5839.50012.46012.46012.4609.91712.36312.36312.364（3）時間步長對溫度的影響如圖 7和表3，根據(jù)圖中曲線可以看出時間步長選擇50s，100s，200s時基本重合，對墻壁溫度影響不大。. . .15.51514.5℃14T=50（s）度溫13.5T=100s13T=200s12.5120246810時間（h）圖7時間步長對溫度的影響表3時間步長對溫度的影響時間(h)T=50（s）時間(h)T=100（s)時間（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.2513.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.6751615）℃14（度13溫1211100 2 4 6 8 10 12 14 16 18 20 22時間（h）圖8墻壁溫度隨時間的變化曲線. . .計算結(jié)果與分析5.1墻壁溫度分析根據(jù)題目中要求，計算寒潮入侵多長時間后墻壁可以感受到外界氣溫的變化，通過建模，方程離散化，最終通過程序求解方程，得到圖 8和表4。由圖可以看出，開始階段，墻壁溫不變，隨著時間的進一步深入，壁溫度開始降低，當很長時間后，溫度變化基本趨于平緩，直到再次平衡。根據(jù)圖 8就可以得到墻壁溫度開始發(fā)生變化的時間。表4墻壁溫度隨時間變化數(shù)據(jù)表時間（h)溫度（℃）時間（h)溫度（℃）時間（h)溫度（℃）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.4583.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 導熱系數(shù)λ對墻壁溫度的影響墻的導熱系數(shù)對表面的影響，在圖9和表5中發(fā)現(xiàn)，導熱系數(shù)對壁溫度影響比較大，=1.2時，溫度下降的趨勢會更快，要比λ=0.85時下降快的多，下降速度更快，更短時間達到穩(wěn)態(tài)，因此導熱系數(shù)越大溫度擴散越快，導熱系數(shù)越小則溫度變化越慢，需要更長時間到達穩(wěn)態(tài)，但是這時對于要求恒溫的空間有好處，受波動影響更小。表5導熱系數(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.37112.62. . .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.5141615λ=0.85λ=1.0） 14 λ=1.2℃13度溫1211100 2 4 6 8 10時刻（h）圖9導熱系數(shù)對墻壁溫度的影響5.3墻外換熱系數(shù) h的影響墻外表面?zhèn)鳠嵯禂?shù)對溫度分布的影響，如圖 10和表6影響不大。1615h=20）14h=35℃h=50（13度溫12111002468101214161820時間（h）圖10墻外換熱系數(shù)對溫度影響表6墻外換熱系數(shù)對溫度影響時刻h=20h=35h=50（h）222·℃W/(m·℃)W/(m·℃)W/(m. . .)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.92611.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ù)對表面溫度影響較大，當傳熱系數(shù)比較小時受墻外流體溫度影響明顯，傳熱系數(shù)越大受墻外流體溫度影響越小，當?shù)搅藰O限大的情況下，表面溫度則等于墻流體溫度，不再受墻外流體溫度影響。1614）12℃（10h=2度溫h=68h=106402468101214161820時刻（h）圖11墻表面?zhèn)鳠嵯禂?shù)對溫度影響表7墻表面?zhèn)鳠嵯禂?shù)對溫度影響時刻h=2h=6h=10（h）222·℃)W/(m·℃)W/(m·℃)W/(m. . .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.11814.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 墻壁厚度δ對壁溫度的影響改變墻壁的厚度δ，壁溫度將發(fā)生變化。在表 8和圖12中可以發(fā)現(xiàn)，不同厚度的墻壁對外界溫度的感應快慢是不一樣的，隨著墻壁厚度的增加，感應到外界溫度變化的時間越來越長，且溫度變化越來越慢，墻壁越厚要越長的時間才能達到新的平衡16151413）12℃（11度100.15m溫0.3m980.45m7602468101214161820時間（h）圖12 墻壁厚度對壁溫度的影響表8墻壁厚度對壁溫度的影響時間（h）0.15m0.3m0.45m01515151.11113.86414.99715. . .2.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.4446.57611.62614.1115.5566.56611.5114.02616.6676.55911.41113.94917.7786.55511.32913.87718.8896.55311.25913.811結(jié)論整個墻壁經(jīng)歷了以下變化過程：首先外壁直接與室外冷空氣接觸，溫度變化很快，隨著時間的推移，墻各點的溫度也開始變化，并影響到右邊墻壁的溫度，慢慢降低。對于墻壁的溫度隨時間的變化，變化趨勢總是由快到慢，最后重新達到穩(wěn)態(tài)。當改變墻壁厚度的時候，隨著墻壁厚度的增加，對于外界溫度的感應是越來越慢。這對于一些對溫度有要求的地反很重要，根據(jù)溫度的變化作出適當?shù)恼{(diào)整措施。參考文獻黃善波中良編著計算傳熱學基礎世銘陶文銓編著傳熱學（第四版）北京高等教育王元明編數(shù)學物理方程與特殊函數(shù)（第三版）高等教育附錄數(shù)學模型的離散過程推導2

p

p1

p

1p11pEpW（1）x2 P

x x e

x w

xxexexwPxw. . .按照Taylor級數(shù)展開法，溫度對時間的偏導有向前差分格式，中心差分格式和向后差分格式，使用向后差分格式。向后差分格式:ppp1PP（2）P所以pp1ap11ppPPEW（3）xxexexwPxwxpp1p11ppEW（4）aPPxexexPxww由xx得x2cx0（5）xPxExWPaPPaEEaWWb（6）aPaEaWaP0aE，aW，aP0cxbaP0TP0（7）xx邊界條件處理：右邊界根據(jù)元體能量平衡法處理：Φw ΦBN-2 N-1 xx圖1右邊界的能量守恒如圖1：BN2Es（）8其中Bhr(TfrTN1)，N2N2N1（9）xcx0N2N1N1N1（10）hr(TfrTN1)x. . .為常數(shù)ac,xx2將式（12-b）(12-c)入式（12-a）以得到:cxxhrN1N2hrTfrcxTN012x2同理可以得到左邊界點離散方程:(hlxcx)T0hTlflT1cxT002x2因此，離散方程為:aPTPaETEaWTWb式中，aPaEaWaP0aEx，aW，aP0cxxbaP0TP0式中，上標“0”表示上一時刻值，為熱擴散系數(shù)，為時間步長。由元體能量平衡法可以分別求出左右邊界節(jié)點的離散方程，左邊界為:(a0ph2)T0T1h2Tf2ap0T002xx2右邊界為:(a0ph1)TN1TN2h1Tf1ap0TN012xx2

（11）（12）（13）（14）（15）（16）（17）（18）（19）程序清單序號程序名稱程序功能計算結(jié)果①表1、4，圖5、8②改變空間步長：表2，圖6用數(shù)值方法求解③改變時間步長：表3，圖7（1）設計程序④改變導熱系數(shù)：表5，圖9溫度⑤改變墻外換熱系數(shù)：表6，圖10⑥改變墻換熱系數(shù)：表7，圖11⑦改變墻壁厚度：表8，圖12. . .（2）驗證程序用分析解求解溫①表1，圖5度2.1 設計程序主要程序段（1）初始時刻的輸入dx=dlt/(N-1);// 空間步長ap0=dc*dx/dtao;// 系數(shù)ap0for(i=0;i<N;i++){T0[i]=Tr-h1*(Tf1-Tr)*(dlt-i*dx)/lw;//初始溫度Tw[i]=Tr-h1*(Tf1-Tr)*(dlt-i*dx)/lw;Tw0[i]=Tr-h1*(Tf1-Tr)*(dlt-i*dx)/lw;// 假設溫