计算传热学程序设计.doc

中国石油大学(华东)储运与建筑工程学院热能与动力工程系计算传热学程序设计设计报告学生姓名:学 号：专业班级：指导教师 2012年 7 月 7 日1、设计题目有一房屋的砖墙厚=0.3 m ，=0.85 W/(m)，c=1.05106 J/( m3K)，室内温度Tf1保持20不变，表面传热系数h1=6W/(m2)。开始时墙的温度处于稳定状态，内墙表面温度Tw1为15寒潮入侵后，室外温度Tf2下降为-10，外墙的表面传热系数为35W/(m2)。试分析寒潮入侵后多少时间内墙壁面方可感受到外界气温的变化。 室外寒流入侵室内0 x 图1 墙壁简化图1.1已知参数壁厚，墙壁导热系数，密度与比热容的乘积，室内和寒潮入侵后室外空气温度，室内空气和外墙的表面传热系数，开始时稳定状态下的内墙表面温度。1.2 求解寒潮入侵多少时间后内墙壁面可感受到外界气温的变化？2 物理与数学模型2.1 物理模型该墙面为常物性，可以假设：（1）其为无限大平面，（2）只有在厚度方向传热，没有纵向传热，则该问题转化为一维常物性无限大平面非稳态导热问题。2.2 数学模型以墙外表面为坐标原点，沿厚度方向为坐标正方向，建立坐标系。基于上述模型，取其在x方向上的微元作为研究对象，则该问题的数学模型可描述如下： （1a）初始条件： T|=0=Tw1-h1(Tf1-Tw1)(-x)/ （1b）在两侧相应的边界条件是第三类边界条件，分别由傅立叶定律可描述如下：左边界： （1c） 右边界： （1d） 3 数值处理与程序设计3.1 数值处理采用外点法用均匀网格对求解区域进行离散化，得到的网格系统如图2所示。一共使用了0N-1共N个节点。节点间距x为: x=N-1 图2 墙壁内的网格划分此例中墙壁导热系数为常值，无源项。则可采用有限体积法对控制方程离散化，得到离散方程为： （2a）式中： （2b） ， （2c） （2d）其中的上标“0”表示此为上一时刻的值,分别为节点所在控制容积左右边界上的导热系数，由于墙壁导热系数不变，故都等于，为时间步长。由元体能量平衡法可以得知左右边界节点的离散方程分别为：左边界节点：ap02+h2+xT0=xT1+h2Tf2+ap02T00 （3）右边界节点： ap02+h1+xTN-1=xTN-2+h1Tf1+ap02TN-10 （4）离散方程的详细推导过程见附录。3.2程序设计由物理模型可以知道本问题为一维导热问题，一维导热问题的离散方程在取遍所有节点后形成的是三对角的代数方程组，采用追赶法进行求解。程序构成和方法:程序由主程序和一个子程序构成。主程序进行变量定义和各已知参数的输入，以及左右边界节点和内部节点控制方程的输入；子程序tdma实现追赶法用来计算每个节点新的温度。Thomas算法求解过程分为两步:消元和回代。消元是从系数矩阵的第二行起，逐一将每一行的非零元素消去一个，使原来的三元方程化为二元方程。消元进行到最后一行时，二元方程就化为一元方程，直接得到最后一个未知数的值。然后逐一往前回代，由各二元方程求出其它未知解。程序特点：该程序有很强的适应性，一维常物性非稳态平壁导热问题都可以使用此程序，只要适当更改边值条件即可。还可以进行修改解决非常物性问题。程序中对输出节点，最大输出量都进行了控制，对计算结果的分析有很大帮助。而且Thoms算法的优点需要内存小，工作量小，程序设计简单。程序流程图：首先对变量赋值，然后由初始条件建立初始温度场，接着从左边界，内部节点，到右边界进行迭代，直到满足精度要求为止，最后输出结果，程序结束。程序流程如下图3。4、模型与程序验证4.1模型本题简化为厚度为2=0.3m的一维非稳态模型如图4所示，初始温度为15，在其中间建立坐标系，左两边为对流换热，且换热系数相同都为h=25 W/(m2)，且流体温度Tf=-10对于x 0，列出其导热微分方程式及定解条件： （5） 开始进行参数赋值用迭代法求解温度场计数并判断是否完成温度场计算输入至屏幕输入至文件程序结束NY图3 程序流程图 （6） （7） （8）引入过余温度： （9）图4 一维导热简化模型0左侧：h Tf右侧：h=25 W/(m2)Tf=-10 x-T0=15直接根据公式得到解析解如下： （10）式中，系数应该使上述无穷级数在是满足初始条件，由傅里叶级数理论可得： （11）是超越方程的根，称为特征根。 （12）其中 。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.0669.0810.1654532.7786.9746.9990.3584743.3335.0825.1120.5903193.8893.3933.4250.9431184.4441.891.9231.7460325000.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.42144图5 平壁中心不同时刻的数值解和分析解 图6 空间步长对墙内壁温度影响表2 空间步长对温度影响数据时间(h) N=51 N=101 N=2010.000 15.000 15.000 15.000 0.500 15.000 15.000 15.000 1.000 14.998 14.998 14.998 1.500 14.980 14.981 14.981 2.000 14.924 14.924 14.924 2.500 14.820 14.820 14.820 3.000 14.675 14.675 14.676 3.500 14.501 14.502 14.502 4.000 14.310 14.311 14.311 4.500 14.111 14.112 14.112 5.000 13.911 13.911 13.912 5.500 13.715 13.715 13.715 6.000 13.525 13.525 13.525 6.500 13.343 13.344 13.344 7.000 13.171 13.172 13.172 7.500 13.009 13.010 13.010 8.000 12.858 12.858 12.858 8.500 12.716 12.716 12.716 9.000 12.583 12.583 12.583 9.500 12.460 12.460 12.460 9.917 12.363 12.363 12.364 （3） 时间步长对温度的影响如图7和表3，根据图中曲线可以看出时间步长选择50s，100s，200s时基本重合，对墙内壁温度影响不大。图7 时间步长对温度的影响表3 时间步长对温度的影响时间(h) T=50（s） 时间(h)T=100（s)时间（h） T=200（s）015015.000 0150.694150.69415.000 0.667151.38914.9881.38914.987 1.33314.9872.08314.9122.08314.910 214.9182.77814.7462.77814.744 2.66714.773.47214.5133.47214.511 3.33314.5584.16714.2454.16714.244 414.3084.86113.9664.86113.966 4.66714.0445.55613.6925.55613.693 5.33313.7816.2513.4316.2513.433 613.5276.94413.1886.94413.190 6.66713.2887.63912.9647.63912.966 7.33313.0668.33312.768.33312.762 812.8629.02812.5749.02812.576 8.66712.675图8 墙内壁温度随时间的变化曲线5 计算结果与分析5.1墙内壁温度分析根据题目中要求，计算寒潮入侵多长时间后内墙壁可以感受到外界气温的变化，通过建模，方程离散化，最终通过程序求解方程，得到图8和表4。由图可以看出，开始阶段，内墙壁温不变，随着时间的进一步深入，内壁温度开始降低，当很长时间后，温度变化基本趋于平缓，直到再次平衡。根据图8就可以得到墙内壁温度开始发生变化的时间。 表4 墙内壁温度随时间变化数据表 时间（h) 温度（） 时间（h) 温度（） 时间（h) 温度（）0.000 156.667 13.28513.333 11.7630.556 157.222 13.09813.889 11.6911.111 14.9977.778 12.92414.444 11.6261.667 14.9678.333 12.76215.000 11.5652.222 14.8838.889 12.61215.556 11.512.778 14.7449.444 12.47316.111 11.4583.333 14.56210.000 12.34516.667 11.4113.889 14.35410.556 12.22717.222 11.3684.444 14.13311.111 12.11817.778 11.3295.000 13.91111.667 12.01818.333 11.2925.556 13.69312.222 11.92618.889 11.2596.111 13.48412.778 11.84119.444 11.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.625.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.514图9 导热系数对墙内壁温度的影响5.3墙外换热系数h的影响墙外表面传热系数对温度分布的影响，如图10和表6影响不大。 图10 墙外换热系数对温度影响表6 墙外换热系数对温度影响时刻（h） h=20 W/(m2) h=35W/(m2) h=50W/(m2)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.03111.62611.44815.55611.91511.5111.33416.66711.81511.41111.23717.77811.7311.32911.15618.88911.65711.25911.0885.4 墙内壁传热热系数的影响由图11可以看出，墙内表面传热系数对内表面温度影响较大，当传热系数比较小时受墙外流体温度影响明显，传热系数越大受墙外流体温度影响越小，当到了极限大的情况下，内表面温度则等于墙内流体温度，不再受墙外流体温度影响。图11 墙内表面传热系数对温度影响表7 墙内表面传热系数对温度影响时刻（h）h=2 W/(m2)h=6 W/(m2)h=10 W/(m2)0.000 15.000 15.000 15.000 1.111 14.994 14.997 14.999 2.222 14.789 14.883 14.955 3.333 14.169 14.562 14.837 4.444 13.292 14.133 14.688 5.556 12.332 13.693 14.540 6.667 11.390 13.285 14.408 7.778 10.507 12.924 14.295 8.889 9.700 12.612 14.201 10.000 8.972 12.345 14.123 11.111 8.318 12.118 14.059 12.222 7.735 11.926 14.007 13.333 7.214 11.763 13.964 14.444 6.749 11.626 13.929 15.556 6.336 11.510 13.900 16.667 5.968 11.411 13.876 17.778 5.640 11.329 13.857 18.889 5.348 11.259 13.841 5.5 墙壁厚度对内壁温度的影响改变墙壁的厚度，内壁温度将发生变化。在表8和图12中可以发现，不同厚度的墙壁对外界温度的感应快慢是不一样的，随着墙壁厚度的增加，感应到外界温度变化的时间越来越长，且温度变化越来越慢，墙壁越厚要越长的时间才能达到新的平衡图12 墙壁厚度对内壁温度的影响表8 墙壁厚度对内壁温度的影响时间（h）0.15m0.3m0.45m01515151.11113.86414.997152.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.8116 结论 (1)整个墙壁经历了以下变化过程：首先外壁直接与室外冷空气接触，温度变化很快，随着时间的推移，墙内各点的温度也开始变化，并影响到右边内墙壁的温度，慢慢降低。(2)对于墙内壁的温度随时间的变化，变化趋势总是由快到慢，最后重新达到稳态。当改变墙壁厚度的时候，随着墙壁厚度的增加，对于外界温度的感应是越来越慢。这对于一些对温度有要求的地反很重要，根据温度的变化作出适当的调整措施。参考文献1 黄善波 刘中良编著 计算传热学基础 2 杨世铭 陶文铨编著 传热学（第四版）北京高等教育出版社3 王元明编 数学物理方程与特殊函数（第三版） 高等教育出版社附录1数学模型的离散过程推导 （1）按照Taylor级数展开法，温度对时间的偏导有向前差分格式，中心差分格式和向后差分格式，使用向后差分格式。向后差分格式: （2）所以 （3） （4）由 得 （5） （6） ， （7）边界条件处理：右边界根据元体能量平衡法处理： xN-2 N-1 xw B图1右边界的能量守恒如图1： （8）其中 ， （9） （10）为常数 , （11）将式（12-b）(12-c)入式（12-a）以得到: （12）同理可以得到左边界点离散方程: （13）因此，离散方程为: （14）式中， （15）， （16） （17）式中，上标“0”表示上一时刻值，为热扩散系数，为时间步长。 由元体能量平衡法可以分别求出左右边界节点的离散方程，左边界为: （18）右边界为: （19）2 程序清单序号程序名称程序功能计算结果（1）设计程序用数值方法求解温度 表1、4，图5、8 改变空间步长：表2，图6 改变时间步长：表3，图7 改变导热系数：表5，图9 改变墙外换热系数：表6，图10 改变墙内换热系数：表7，图11 改变墙壁厚度：表8，图12（2）验证程序用分析解求解温度 表1，图52.1 设计程序主要程序段（1） 初始时刻的输入 dx=dlt/(N-1);/空间步长ap0=dc*dx/dtao;/系数ap0for(i=0;iN;i+)T0i=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;/假设温度场（2