计算传热学课程设计

1、中国石油大学(华东)储建学院热能与动力工程系计算传热学程序设计设计报告学生姓名： 学 号：专业班级：热能与动力工程指导教师： 2011年 7 月 5 日目录1 设计题目22 物理与数学模型22.1 物理模型22.2 数学模型23 数值处理与程序设计33.1 数值处理33.1.1 数值处理与方程离散化33.1.2 离散方程的导出43.2 程序设计64 模型与程序验证75 计算结果与分析96 结论107 参考文献118 附录11附录一 数据表格11附录二 主要程序121 设计题目有一冷柜初始时处于均匀的温度20。开启压缩机后，冷冻室（包括柜门的内表面）以均匀地速度降温，温降速度为18/h。柜门尺寸

2、为1.2x1.2，保温材料厚0.08m，导热系数为0.04W/(m·)。假设外表面包层很薄，其热阻可以忽略。柜门外受空气自然对流和环境的辐射换热。自然对流可按下式计算， W/(·)其中L为门高。表面发射率为0.8。通过柜门的导热可按一维问题处理。试计算压缩机启动两个小时内通过柜门的冷量损失。2 物理与数学模型2.1 物理模型由题可知，柜门左边界温度均匀下降，右边界置于恒温中，并有对流和辐射散热，由于柜门沿高度方向没有散热，所以沿高度方向没有温度的变化，即只有沿厚度方向的温度变化，所以可以视为一维问题处理，同时可以得出左右边界都是第三类边界情况。2.2 数学模型以柜门左侧为原

3、点，沿柜门厚度方向为坐标轴正方向建立坐标轴。则有=0.08m，基于上述物理模型，只对厚度x方向分析，可以得到该问题的数学模型为 (1)初始条件为： (2)边界条件为： (3)X0 图1 柜门模型3 数值处理与程序设计3.1 数值处理3.1.1 数值处理与方程离散化采用外节点法用均匀网格对求解区域进行离散化，得到的网格系统如图2所示。图中为了和C语言中数组定义一致，节点标号从0到N-1一共有N个节点。节点间距x为 图2 宽度x方向的网格划分简化图采用有限差分法或者有限体积法对控制方程进行离散化，得到的离散方程为 (4a)式中： (4b)， , (4c) (4d)其中，上标“0”表示上一时刻值，为

4、节点所在控制容积界面上的导热系数，为时间步长。由于柜门的材料的为常数，在离散方程中取值时，w=e=0.04W/(m)。由于左侧界面柜门温度均匀下降，故有: (5)由题意知右边界为第三类边界条件，故有 (6)由元体能量平衡法与源项线性化可以得到右边界节点的离散方程，为： (7)其中 由于柜门的材料的为常数，在离散方程中取值时，w=e=。3.1.2 离散方程的导出-具体的推导过程如下为了得到离散方程，除了要对空间区域进行离散化之外，还需要对时间坐标进行离散化。对时间进行离散化：从初始时刻开始，取定一个时间步长构成如下的时间序列即可， （p=0，1，2，3，） (8)其中p=0对应于初始时刻，后面为

5、方便起见将p时刻简称为p时刻。内部节点P及其与邻点的关系仍如图2所示。求解区域内的任一点P，在p时刻，它仍然满足控制方程，即 (9)上述方程中，温度对空间坐标的导数项可以表示为 (10)按照Taylor级数展开法，温度对时间的偏导数采用向后差分格式： (11)将(12)带入到差分方程(11)中去，并定义 ， (12)整理后得到隐式差分格式的离散方程为： (13)上式即是该题目中一维非稳态导热方程的离散方程。下面来分别计算左右两边界的边界条件：对于左边界节点，温度随时间均匀变化，故有 (14)wB图3 右边界节点的能量守恒对于右边界节点，采用元体能量平衡法进行离散。如图3，根据能量守恒，有 (1

6、5)w为由左界面进入控制容积的热量。假设节点间的温度按线性分布，根据Fourier定律有: (16)B为由边界流入控制容积的热量。由对流和辐射换热方程有: (17)式中，s为控制容积在对应时间段内的热力学能变化率，它可以表示为 (18)将1618式代入15式中，整理得到 (19)令 (20)故，运用源项线性化可得 (21)将式（21）代入（19），可得 (24)3.2 程序设计对离散方程的求解，在此处是采用追赶法进行求解。运用追赶法编写程序，会使程序编写简单，同时，使程序中的每条语句的作用明确、容易读懂。图4为程序流程图，根据程序可以求出柜门的温度分布。4 模型与程序验证由于该题目难以求取解析

7、解，所以这里不直接求取该题目的解析解，而是通过其他方式来验证程序和输出结果的正确性。我们将题目中的柜门左侧设为突然变为-16，右侧柜面处恒定为20，由于原柜门初始温度都是20，所以当处于这种假设情况时，柜门最终会处于热平衡状态，这种情况时可以通过程序和解析方程式求出柜门的热量变图4一维非稳态导热问题求解程序流程图化的，所以用两种方法分别求出柜门的热量变化，比较这两个解，如果在误差允许的范围内是相等的，就能证明程序和输出结果在当前假设情况的正确性，也就从侧面表明在原题设的条件下是正确的。验证的过程为：在程序中令：,左边界改为f0=-16,右边界节点方程改为AN-1=0.0;BN-1=1.0;CN

8、-1=0.0;fN-1=20;此时可以得到程序的输出的温度曲线（即图5）中，在T=7200s时刻，柜门内的温度已经到达稳态，从初始时刻到T=7200s时刻冷量损失Q1=233.893KJ解析解过程：达到稳态后：所以，柜门内的温度分布为： 此时，柜门的热量变化为： 解得Q2=233.280KJ两者的误差：由计算得出的误差Eps=0.02627%<5%，可以得出在误差范围内Q2=Q1，所以验证了在假设情况下程序的计算结果是正确的，也就证明了，在题设条件下程序的计算结果是正确的。图5 柜门在验证程序情况下，不同时刻的温度分布5 计算结果与分析利用程序对问题进行问题模型进行计算，其中=100s,

9、 x=0.008m，可以得到柜门的温度场，利用2小时内通过柜门损失到空气中得冷量Q1=和2小时内存储在柜门中得冷量Q2=可以计算出2小时内的冷量损失Q=Q1+Q2=205.752171KJ，并将柜门部分温度分布结果列于附录的表中由表可见，受左表面均匀降温和受右边界对流及辐射换热的影响，柜门内各点的温度均随时间的变化而不断变化，主要特点是（1）由图6可得同一位置温度随时间变化的规律，由于左边界温度均匀下降，使得左侧柜门面温度均匀下降，所以柜门左侧的温度变化曲线应为一条斜直线，而其他位置温度变化因为有导热的影响而不为直线。图6 同一位置温度随时间的变化（2）由图7不同时刻的柜门温度场可得，柜门左侧

10、的温度最低，随着位置的变化温度不断的升高，而在0时刻，因为出于初始状态，所以柜门内的温度不变且都为初始状态的温度。（3）由图8不同导热系数下的冷量损失得，因为柜门为保温材料，所以总的导热过程热阻主要由柜门的导热系数决定，不同的导热系数下损失的冷量将会有很大的差异，取不同的导热系数，分别计算相应导热系数下的冷量损失。比较发现，在导热系数小于10w/(mk)时，冷量损失随导热系数变化而剧烈变化，之后冷量损失基本不再变化。这表明在导热系数小于10 w/(mk)时，对冷量损失的影响因素中，材料的导热系数其主导作用；而当导热系数大于10 w/(mk)时，在冷量的损失因素中，导热系数将不再其主要的作用，由

11、此可见冷柜中对冷量损失起主要影响因素的为柜门的保温材料种类。图7 图7 不同时刻的柜门温度场图8 导热系数不同时的冷量损失6 结论通过上述的分析,可以得到以下结论:1两个小时的冷量损失为Q=205.752171KJ2用有限差分法处理一维非稳态的导热问题，可得到与分析解吻合的数值解。可见数值方法具有一定的可靠性。3在冷柜保温材料中，材料本身的导热系数对冷量的损失有较大的影响。7 参考文献1 黄善波，刘中良编.计算传热学基础.东营：中国石油大学(华东)胶印，2008.2 杨世铭，陶文铨编著.传热学(第四版).北京：高等教育出版社，2007.3 谭浩强著.C程序设计(第二版).北京：清华大学出版社，

12、2006.8 附录附录一 程序列表1追赶法计算柜门温度分布及散热量 对应表1，图6，图7。改变程序算柜门温中的h值，得表2，图8。3验证程序 对应表3，图5。附录二 数据表格表1 =100s, x=0.008m,柜门内温度分布 x00.0080.0160.0240.0320.040.0480.0560.0640.0720.080202020202020202020202018001113.7915.8317.2918.2918.9519.3819.6519.8119.919.96360026.119.4312.0814.1615.771717.9418.6519.219.655400-7-1.

13、852.496.139.1511.6413.715.416.8218.0419.137200-16-9.93-4.67-0.143.767.119.9912.514.716.6818.52表2 =100s, x=0.008m，不同h下的冷量损失h/ w(mk)0.020.040.080.101.0010.0050.00100.00 150.00 200.00 Q/KJ157.59205.75259.89279.83580.09742.98763.76766.45 767.36 767.81 表3 =100s, x=0.008m验证程序的输出结果 x0.0000.0080.0160.0240.0

14、320.0400.0480.0560.0640.0720.080020.0020.0020.0020.0020.0020.0020.0020.0020.0020.0020.001800-16.00-9.50-3.332.206.9310.7613.7415.9917.6518.9220.003600-16.00-11.36-6.82-2.481.605.368.7911.9114.7617.4320.005400-16.00-12.00-8.05-4.16-0.383.286.8210.2413.5516.8020.007200-16.00-12.25-8.51-4.80-1.132.496.

15、079.6013.0916.5520.00附录三 主要程序部分主程序loop:/左边界节点A0=0.0;B0=1.0;C0=0.0;f0=293-18*Tao/3600; /左侧壁面温度从20度按每小时18度均匀下降/内部节点for(i=1;i<N;i+)Ai=lmda/dx;Ci=lmda/dx;Bi=Ai+Ci+ap0;fi=ap0*T0i;/右边界节点h=dlcrxs(TwN-1);AN-1=lmda/dx;BN-1=ap0/2+h+lmda/dx-fsl*htcs*T0N-1*T0N-1*T0N-1;CN-1=0;fN-1=h*Tf+ap0*T0N-1/2+fsl*htcs*Tf

16、*Tf*Tf*Tf-2*fsl*htcs*T0N-1*T0N-1*T0N-1*T0N-1;/调用TDMA计算新的温度tdma(N,A,B,C,f,Tw);/计算迭代间的最大误差Er=0.0;for(i=0;i<N;i+)Err=fabs(Twi-Tw0i);if(Err>Er)Er=Err;Iter=Iter+1;if(Er>Eps)for(i=0;i<N;i+)Tw0i=Twi;goto loop;/满足精度输出本时刻参考点的温度printf("Ta0=%6.3f Iter=%3d T2=%7.3f T8=%7.3fn",Tao/3600,Iter,Tw1,Tw7);/保存给定时刻的温度if(fmod(It,Nout)=0.0)Iout=Iout+1;TimeIout=Tao;for (i=0;i<N;i+)ToutIouti=Twi;/设置下一时刻计算的温度初始值for(i=0;i<N;i+)T0i=Twi;while(Tao<Taomax);追赶法void tdma(int N,double a,double b,double c,double f,double x)int i;double PNT,QNT;P0=c0/