计算传热学程序的设计说明

中国石油大学(华东)储运与建筑工程学院热能与动力工程系《计算传热学程序设计》设计报告学生:学 号：专业班级：指导教师2012 年7 月7 日. . .1、设计题目有一房屋的砖墙厚δ=0.3m，λ=0.85W/(m·℃)，ρc6J/(m3·K)，=1.05×10室温度Tf1保持℃不变，表面传热系数h12·℃。开始时墙的温度处于稳定状20=6W/(m)态，墙表面温度Tw为℃寒潮入侵后，室外温度T下降为-10℃，外墙的表面传热系115f2数为35W/(m2·℃)。试分析寒潮入侵后多少时间墙壁面方可感受到外界气温的变化。室外室内寒流入侵0 x图1墙壁简化图1.1已知参数壁厚，墙壁导热系数，密度与比热容的乘积，室和寒潮入侵后室外空气温度，室空气和外墙的表面传热系数，开始时稳定状态下的墙表面温度。1.2 求解寒潮入侵多少时间后墙壁面可感受到外界气温的变化？物理与数学模型2.1 物理模型该墙面为常物性，可以假设：（ 1）其为无限大平面，（ 2）只有在厚度方向传热，没有纵向传热，则该问题转化为一维常物性无限大平面非稳态导热问题。2.2 数学模型以墙外表面为坐标原点，沿厚度方向为坐标正方向，建立坐标系。基于上述模型，取其在x方向上的微元作为研究对象，则该问题的数学模型可描述如下：cT(T)（1a）xx初始条件：??|??=0=????1-??1(????1-????1)(δ-x)/λ（1b）. . .在两侧相应的边界条件是第三类边界条件，分别由傅立叶定律可描述如下：左边界：T（1c）x0h2(Tx0Tf2)X右边界：Txh1(TxTf1)（1d）X数值处理与程序设计3.1 数值处理采用外点法用均匀网格对求解区域进行离散化，得到的网格系统如图 2所示。一共使用了0~N-1共N个节点。节点间距δx为:δx=

????-1图2墙壁的网格划分此例中墙壁导热系数为常值，无源项。则可采用有限体积法对控制方程离散化，得到离散方程为：apTpaETEaWTWb（2a）式中：aPaEaWaP0（2b）aE，aWa0cx（）xxba0pTp0（2d）其中的上标“0”表示此为上一时刻的值,分别为节点所在控制容积左右边界上的导热系数，由于墙壁导热系数不变，故都等于λ，△τ为时间步长。由元体能量平衡法可以得知左右边界节点的离散方程分别为：. . .左边界节点：0λλ0??????????????0(+??2+)??0=+2????2????（3）右边界节点：(????0+??1+λλ????+????0??02)????-1=??????2+??2??1????（4）离散方程的详细推导过程见附录。3.2程序设计由物理模型可以知道本问题为一维导热问题，一维导热问题的离散方程在取遍所有节点后形成的是三对角的代数方程组，采用追赶法进行求解。程序构成和方法:程序由主程序和一个子程序构成。主程序进行变量定义和各已知参数的输入，以及左右边界节点和部节点控制方程的输入；子程序tdma实现追赶法用来计算每个节点新的温度。Thomas算法求解过程分为两步:消元和回代。消元是从系数矩阵的第二行起，逐一将每一行的非零元素消去一个，使原来的三元方程化为二元方程。消元进行到最后一行时，二元方程就化为一元方程，直接得到最后一个未知数的值。然后逐一往前回代，由各二元方程求出其它未知解。程序特点：该程序有很强的适应性，一维常物性非稳态平壁导热问题都可以使用此程序，只要适当更改边值条件即可。还可以进行修改解决非常物性问题。程序中对输出节点，最大输出量都进行了控制，对计算结果的分析有很大帮助。而且Thoms算法的优点需要存小，工作量小，程序设计简单。程序流程图：首先对变量赋值，然后由初始条件建立初始温度场，接着从左边界，部节点，到右边界进行迭代，直到满足精度要求为止，最后输出结果，程序结束。程序流程如下图3。4、模型与程序验证4.1模型本题简化为厚度为 2 =0.3m的一维非稳态模型如图 4所示，初始温度为 15℃，在其中间建立坐标系，左两边为对流换热，且换热系数相同都为 h=25W/(m2·℃)，且流体温度Tf=-10℃对于x 0，列出其导热微分方程式及定解条件：. . .Ta2T(0x,0)（5）x2开始进行参数赋值

输入至屏幕用迭代法 输入至文件求解温度N程序结束计数并判断是否完成温度场计算Y图3程序流程图aT(x,0)T0(0x)（6）cT(x,)x00（7）xT(x,)hT(,)（8）Txx引入过余温度：. . .T(x,)T（9）T0=15℃左侧：hTf

右侧：h=25W/(m2·℃)Tf=-10℃-δ 0 δx图4一维导热简化模型直接根据公式得到解析解如下：(, ) Cnexp( n2Fo)cos( n ) （10）n1式中，Fo a2, x，系数Cn应该使上述无穷级数在 0是满足初始条件，由傅里叶级数理论可得：Cn2sinn（11）cosnsinnn是超越方程的根，称为特征根。tannBi,n1,2,⋯（12）nh其中 Bi 。4.2程序验证（1）由模型可以得到相关信息然后进行编程，同等时间下计算出中心处温度的解析解和数值解进行比较，数据记录在表 1。然后计算出相对误差，作图 5，观察数值解与分析解的比较曲线。由图表中可以发现，平壁中心不同时刻温度值的分析解和数值解相差不是很大，二者吻合的比较好，可以说明所编制的数值解法的程序是正确的。 相对误差先增大后减小，增大的原因是此时温度接近零度，相对误差的基数比较小，所以造成相对误差较大，但是此时的绝对误差并不大，在合理围，所以除去个别点外，都满足误差小于百分之 1。. . .可以验证所编数值解法的程序是正确的。（2）空间步长对墙壁的温度影响如图 6及表2。在程序编写过程中用网格节点数对空间步长进行控制，为了观察空间步长对墙壁温度的影响，表中选择了三个不同的空间步长，分别为选取 51，101，201个网格节点，则相应的空间步长为 0.006，0.003，0.0015。根据不同步长时温度的变化曲线可以看出，空间步长对墙壁的影响不大，当空间步长控制在合理围时可以忽略空间步长的影响。表1分析解与数值解比较时间（h）分析解（℃）数值解（℃）相对误差（%）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（数值解（℃）度4温20-2 0 1 2 3 4 5 6 7 8 9 10-4时间（h）-6-8图5平壁中心不同时刻的数值解和分析解. . .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空间步长对温度影响数据时间(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，根据图中曲线可以看出时间步长选择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墙壁温度随时间的变化曲线. . .计算结果与分析5.1墙壁温度分析根据题目中要求，计算寒潮入侵多长时间后墙壁可以感受到外界气温的变化，通过建模，方程离散化，最终通过程序求解方程，得到图 8和表4。由图可以看出，开始阶段，墙壁温不变，随着时间的进一步深入，壁温度开始降低，当很长时间后，温度变化基本趋于平缓，直到再次平衡。根据图 8就可以得到墙壁温度开始发生变化的时间。表4墙壁温度随时间变化数据表时间（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 导热系数λ对墙壁温度的影响墙的导热系数对表面的影响，在图9和表5中发现，导热系数对壁温度影响比较大，=1.2时，温度下降的趋势会更快，要比λ=0.85时下降快的多，下降速度更快，更短时间达到稳态，因此导热系数越大温度扩散越快，导热系数越小则温度变化越慢，需要更长时间到达稳态，但是这时对于要求恒温的空间有好处，受波动影响更小。表5导热系数对墙壁温度的影响时刻（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导热系数对墙壁温度的影响5.3墙外换热系数 h的影响墙外表面传热系数对温度分布的影响，如图 10和表6影响不大。1615h=20）14h=35℃h=50（13度温12111002468101214161820时间（h）图10墙外换热系数对温度影响表6墙外换热系数对温度影响时刻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 墙壁传热热系数的影响由图11可以看出，墙表面传热系数对表面温度影响较大，当传热系数比较小时受墙外流体温度影响明显，传热系数越大受墙外流体温度影响越小，当到了极限大的情况下，表面温度则等于墙流体温度，不再受墙外流体温度影响。1614）12℃（10h=2度温h=68h=106402468101214161820时刻（h）图11墙表面传热系数对温度影响表7墙表面传热系数对温度影响时刻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 墙壁厚度δ对壁温度的影响改变墙壁的厚度δ，壁温度将发生变化。在表 8和图12中可以发现，不同厚度的墙壁对外界温度的感应快慢是不一样的，随着墙壁厚度的增加，感应到外界温度变化的时间越来越长，且温度变化越来越慢，墙壁越厚要越长的时间才能达到新的平衡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结论整个墙壁经历了以下变化过程：首先外壁直接与室外冷空气接触，温度变化很快，随着时间的推移，墙各点的温度也开始变化，并影响到右边墙壁的温度，慢慢降低。对于墙壁的温度随时间的变化，变化趋势总是由快到慢，最后重新达到稳态。当改变墙壁厚度的时候，随着墙壁厚度的增加，对于外界温度的感应是越来越慢。这对于一些对温度有要求的地反很重要，根据温度的变化作出适当的调整措施。参考文献黄善波中良编著计算传热学基础世铭陶文铨编著传热学（第四版）北京高等教育王元明编数学物理方程与特殊函数（第三版）高等教育附录数学模型的离散过程推导2

p

p1

p

1p11pEpW（1）x2 P

x x e

x w

xxexexwPxw. . .按照Taylor级数展开法，温度对时间的偏导有向前差分格式，中心差分格式和向后差分格式，使用向后差分格式。向后差分格式:ppp1PP（2）P所以pp1ap11ppPPEW（3）xxexexwPxwxpp1p11ppEW（4）aPPxexexPxww由xx得x2cx0（5）xPxExWPaPPaEEaWWb（6）aPaEaWaP0aE，aW，aP0cxbaP0TP0（7）xx边界条件处理：右边界根据元体能量平衡法处理：Φw ΦBN-2 N-1 xx图1右边界的能量守恒如图1：BN2Es（）8其中Bhr(TfrTN1)，N2N2N1（9）xcx0N2N1N1N1（10）hr(TfrTN1)x. . .为常数ac,xx2将式（12-b）(12-c)入式（12-a）以得到:cxxhrN1N2hrTfrcxTN012x2同理可以得到左边界点离散方程:(hlxcx)T0hTlflT1cxT002x2因此，离散方程为:aPTPaETEaWTWb式中，aPaEaWaP0aEx，aW，aP0cxxbaP0TP0式中，上标“0”表示上一时刻值，为热扩散系数，为时间步长。由元体能量平衡法可以分别求出左右边界节点的离散方程，左边界为:(a0ph2)T0T1h2Tf2ap0T002xx2右边界为:(a0ph1)TN1TN2h1Tf1ap0TN012xx2

（11）（12）（13）（14）（15）（16）（17）（18）（19）程序清单序号程序名称程序功能计算结果①表1、4，图5、8②改变空间步长：表2，图6用数值方法求解③改变时间步长：表3，图7（1）设计程序④改变导热系数：表5，图9温度⑤改变墙外换热系数：表6，图10⑥改变墙换热系数：表7，图11⑦改变墙壁厚度：表8，图12. . .（2）验证程序用分析解求解温①表1，图5度2.1 设计程序主要程序段（1）初始时刻的输入dx=dlt/(N-1);// 空间步长ap0=dc*dx/dtao;// 系数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;// 假设温