计算传热学程序设计

1、中国石油大学 （华东 ）储运与建筑工程学院热能与动力工程系计算传热学程序设计设计报告学生姓名 : 学 号： 专业班级：指导教师2012 年 7 月 7 日1、设计题目有一房屋的砖墙厚 S = m ,入=W/（mC）, pc=xiO6J/（ m 3 - K），室内温度Tfi 保持20C不变，表面传热系数hi=6W/（mC ）。开始时墙的温度处于稳定状态，内墙表 面温度Tw为15C寒潮入侵后，室外温度Tf2下降为-10 C,外墙的表面传热系数为 35W（m2C ）。试分析寒潮入侵后多少时间内墙壁面方可感受到外界气温的变化。已知参数壁厚，墙壁导热系数，密度与比热容的乘积，室内和寒潮入侵后室外空气温度

2、，室 内空气和外墙的表面传热系数，开始时稳定状态下的内墙表面温度。求解寒潮入侵多少时间后内墙壁面可感受到外界气温的变化2物理与数学模型物理模型该墙面为常物性，可以假设：（1）其为无限大平面，（2）只有在厚度方向传热, 没有纵向传热，则该问题转化为一维常物性无限大平面非稳态导热问题。数学模型以墙外表面为坐标原点，沿厚度方向为坐标正方向，建立坐标系。基于上述模型， 取其在x方向上的微元作为研究对象，则该问题的数学模型可描述如下：(1a)T)初始条件:T | % - q = T.i - hj (Tf i - T刚)：6 - X)/ X(1b)在两侧相应的边界条件是第三类边界条件，分别由傅立叶定律可描

3、述如下:左边界：-1x0h2(Tx 0 Tf2)(1c)X右边界：|xhi(TxTfi)(1d)X3数值处理与程序设计数值处理采用外点法用均匀网格对求解区域进行离散化，得到的网格系统如图2所示。一共使用了 0N-1共N个节点。节点间距Sx为：%二占图2墙壁内的网格划分此例中墙壁导热系数为常值，无源项。则可采用有限体积法对控制方程离散化，得到离散方程为:apTpaETEawTwb(2a)式中:aP aEaW0aP(2b)(2c)0C X-，aPX(2d)0 0bapTp其中的上标“ 0”表示此为上一时刻的值,分别为节点所在控制容积左右边界上的导热系数，由于墙壁导热系数不变，故都等于入，t为时间步

4、长。由元体能量平衡法可以得知左右边界节点的离散方程分别为: 左边界节点:(3)(I + +右边界节点:(2 + h丄 + 臥)Tn - 1 二阪 - 2 + hf丄 4 製n - 1(4)离散方程的详细推导过程见附录。程序设计由物理模型可以知道本问题为一维导热问题，一维导热问题的离散方程在取遍所有节点后形成的是三对角的代数方程组，采用追赶法进行求解。程序构成和方法：程序由主程序和一个子程序构成。主程序进行变量定义和各已知 参数的输入，以及左右边界节点和内部节点控制方程的输入；子程序 tdma实现追赶法 用来计算每个节点新的温度。Thomas算法求解过程分为两步：消元和回代。消元是从系 数矩阵的

5、第二行起，逐一将每一行的非零元素消去一个，使原来的三元方程化为二元方 程。消元进行到最后一行时，二元方程就化为一元方程，直接得到最后一个未知数的值。 然后逐一往前回代，由各二元方程求出其它未知解。程序特点：该程序有很强的适应性，一维常物性非稳态平壁导热问题都可以使用此 程序，只要适当更改边值条件即可。还可以进行修改解决非常物性问题。程序中对输出节点，最大输出量都进行了控制，对计算结果的分析有很大帮助。而且 Thoms算法的优点需要内存小，工作量小，程序设计简单。程序流程图：首先对变量赋值，然后由初始条件建立初始温度场，接着从左边界， 内部节点，到右边界进行迭代，直到满足精度要求为止，最后输出结

6、果，程序结束。程 序流程如下图3。4、模型与程序验证模型本题简化为厚度为2 =的一维非稳态模型如图4所示，初始温度为15C，在其中 间建立坐标系，左两边为对流换热，且换热系数相同都为h=25 W/(m2C )，且流体温度Tf=-10 C对于x 0,列出其导热微分方程式及定解条件：2a -2 (0 x ,0)(5)x图3程序流程图a cT(x,O) T(0(6)T(x,)(7)引入过余温度:h T(,)T(x,)(8)1i、1i左侧：h TfIiiii iiIi iiIi-0右侧：2h=25 W/(m C)Tf=-10C厂x图4 一维导热简化模型To=15C直接根据公式得到解析解如下:(,)0.

7、Cnexp(n 12nFO)COS( n )(10)式中，FO ay,-，系数Cn应该使上述无穷级数在0是满足初始条件，由傅里叶级数理论可得：2sin n(11)nn COSnSin nn是超越方程的根，称为特征根。tan n Bi,nn1,2,（其中Bi程序验证（1）由模型可以得到相关信息然后进行编程，同等时间下计算出中心处温度的解 析解和数值解进行比较，数据记录在表 1。然后计算出相对误差，作图5,观察数值解 与分析解的比较曲线。由图表中可以发现，平壁中心不同时刻温度值的分析解和数值解相差不是很大，二 者吻合的比较好，可以说明所编制的数值解法的程序是正确的。 相对误差先增大后减小， 增大的

8、原因是此时温度接近零度，相对误差的基数比较小，所以造成相对误差较大，但 是此时的绝对误差并不大，在合理范围内，所以除去个别点外，都满足误差小于百分之 1。可以验证所编数值解法的程序是正确的。（2）空间步长对墙内壁的温度影响如图 6及表2。在程序编写过程中用网格节点根据数对空间步长进行控制，为了观察空间步长对墙内壁温度的影响，表中选择了三个不同 的空间步长，分别为选取51，101，201个网格节点，则相应的空间步长为， 。 不同步长时温度的变化曲线可以看出，空间步长对内墙壁的影响不大，当空间步长控制 在合理范围时可以忽略空间步长的影响。表1分析解与数值解比较时间（h）分析解（C）数值解（C）相对

9、误差（0151505. 00图6空间步长对墙内壁温度影响表2空间步长对温度影响数据时间（h）N=51N=101N=201图5平壁中心不同时刻的数值解和分析解（3）时间步长对温度的影响如图7和表3,根据图中曲线可以看出时间步长选择 50s, 100s, 200s时基本重合，对墙内壁温度影响不大。图7时间步长对温度的影响表3时间步长对温度的影响时间（h）T=50 （s） 时间（h）T=100（ s） 时间（h）T=200（ s）015001515152468图8墙内壁温度随时间的变化曲线5计算结果与分析墙内壁温度分析根据题目中要求，计算寒潮入侵多长时间后内墙壁可以感受到外界气温的变化，通 过建模，

10、方程离散化，最终通过程序求解方程，得到图 8和表4。由图可以看出，开始 阶段，内墙壁温不变，随着时间的进一步深入，内壁温度开始降低，当很长时间后，温度变化基本趋于平缓，直到再次平衡。根据图 8就可以得到墙内壁温度开始发生变化的 时间。表4墙内壁温度随时间变化数据表时间（h）温度C）时间（h）温度C）时间（h） 温度C）1515导热系数入对墙内壁温度的影响墙的导热系数对内表面的影响，在图 9和表5中发现，导热系数对内壁温度影响比 较大，入=时，温度下降的趋势会更快，要比 入=时下降快的多，下降速度更快，更短时 间内达到稳态，因此导热系数越大温度扩散越快，导热系数越小则温度变化越慢，需要 更长时间

11、到达稳态，但是这时对于要求恒温的空间有好处，受波动影响更小。表5导热系数对墙内壁温度的影响时刻（h）入二入二入二0151515151515图9导热系数对墙内壁温度的影响墙外换热系数h的影响墙外表面传热系数对温度分布的影响，如图 10和表6影响不大。图10墙外换热系数对温度影响表6墙外换热系数对温度影响h=502时刻 h=20h=35W/(mC(h) w/(mC ) w/(mC )015151510墙内壁传热热系数的影响由图11可以看出，墙内表面传热系数对内表面温度影响较大，当传热系数比较小 时受墙外流体温度影响明显，传热系数越大受墙外流体温度影响越小，当到了极限大的 情况下，内表面温度则等于墙

12、内流体温度，不再受墙外流体温度影响。图11墙内表面传热系数对温度影响表7墙内表面传热系数对温度影响时刻h=2h=6h=102 2 2(h) W/(m ) W/(m ) W/(m )墙壁厚度S对内壁温度的影响改变墙壁的厚度3,内壁温度将发生变化。在表8和图12中可以发现，不同厚度 的墙壁对外界温度的感应快慢是不一样的，随着墙壁厚度的增加，感应到外界温度变化 的时间越来越长，且温度变化越来越慢，墙壁越厚要越长的时间才能达到新的平衡图12 墙壁厚度对内壁温度的影响表8墙壁厚度对内壁温度的影响时间（h）0151515156结论(1) 整个墙壁经历了以下变化过程：首先外壁直接与室外冷空气接触，温度变化很

13、 快，随着时间的推移，墙内各点的温度也开始变化，并影响到右边内墙壁的温度，慢慢 降低。(2) 对于墙内壁的温度随时间的变化，变化趋势总是由快到慢，最后重新达到稳态。 当改变墙壁厚度的时候，随着墙壁厚度的增加，对于外界温度的感应是越来越慢。这对 于一些对温度有要求的地反很重要，根据温度的变化作出适当的调整措施。参考文献1 黄善波 刘中良编著 计算传热学基础2 杨世铭陶文铨编著传热学(第四版)北京高等教育出版社3 王元明编数学物理方程与特殊函数(第三版)高等教育出版社附录1数学模型的离散过程推导w1 Xe1 X1按照Taylor级数展开法，温度对时间的偏导有向刖差分格式，中心差分格式和向 后差分格

14、式，使用向后差分格式。向后差分格式：3)W PW XePEXaxxpp 1aPP由x x得x 2pxp E11p Pp wx ex ex wx wEwc x 0PxxawEw bap边界条件处理:如图1:其中为常数aPaEoaEawapaEaWx右边界根据元体能量平衡法处理:0aPb a；TP1Iy2/B-1，x x图1右边界的能量守恒BN 2EsBhr (TfrTn1),N 2N 1N2xN 2N 1. z-rhr (TfrTn1)c x0N 1N 1xx2将式（12-b）（12- c）入式（12-a）以得到:c x hhrN2xhrTfr同理可以得到左边界点离散方程：hTfl-T1x-T0

15、0(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)因此，离散方程为NeTeb式中,(14)(15)aE -，aW，aPxxb a；TP式中，上标“ 0”表示上一时刻值，为热扩散系数，为时间步长0ap Ne aw apo c x(17)由元体能量平衡法可以分别求出左右边界节点的离散方程, 左边界为：0ap0aP 0T02(Th2 )T。xT1h2Tf2x右边界为：00h1)T N 1-Tn 2h1Tf1aPtN2xx2(18)(19)2程序清单序号程序名称程序功能计算结果(1)设计程序用数值方法求解 温度 表1、4,图5、8 改变空间步长：表2,图6 改变时间步长：表3,图7 改变导热系数：表5,图9 改变墙外换热系数：表6,图10 改变墙内换热系数：表7,图11 改变墙壁厚度：表8,图12(2)验证程序用分析解求解温 度表1,图5设计程序主要程序段(1)初始时刻的输入dx=dlt/(N-1);h2+lw/dx;A0=;f0=h2*Tf2+ap0*T00/;CN-1=;fN-1=ap0*T0N-1/2.+h1*Tf1(2)追赶法程序tdma(N,A,B,C,f,Tw);fa=a0*ta n(a0)-Bi;fc=c*ta n(c)-Bi;if(