计算传热学程序设计

中国石油大学(华东)储运与建筑工程学院热能与动力工程系《计算传热学程序设计》设计报告学生姓名:学号：专业班级：指导教师2012年 7 月 7 日1、设计题目有一房屋的砖墙厚δ= m ，λ= W/(m·℃)，ρc=×106 J/( m3·K)，室内温度T f1保持20℃不变，表面传热系数h1=6W/(m2·℃)。

开始时墙的温度处于稳定状态，内墙表面温度Tw1为15℃寒潮入侵后，室外温度T f2下降为-10℃，外墙的表面传热系数为35W/(m2·℃)。

试分析寒潮入侵后多少时间内墙壁面方可感受到外界气温的变化。

图1 墙壁简化图 已知参数壁厚，墙壁导热系数，密度与比热容的乘积，室内和寒潮入侵后室外空气温度，室内空气和外墙的表面传热系数，开始时稳定状态下的内墙表面温度。

求解寒潮入侵多少时间后内墙壁面可感受到外界气温的变化2 物理与数学模型物理模型该墙面为常物性，可以假设：（1）其为无限大平面，（2）只有在厚度方向传热，没有纵向传热，则该问题转化为一维常物性无限大平面非稳态导热问题。

数学模型以墙外表面为坐标原点，沿厚度方向为坐标正方向，建立坐标系。

基于上述模型，取其在x 方向上的微元作为研究对象，则该问题的数学模型可描述如下：T()T cx x ρλτ∂∂∂=∂∂∂ （1a ）初始条件：（1b ）在两侧相应的边界条件是第三类边界条件，分别由傅立叶定律可描述如下： 左边界：0202()x f x T h T T X==∂-λ=-∂ （1c ）室外寒流入侵室内0 x右边界：11()x f x T h T T X=δ=δ∂-λ=-∂ （1d ）3 数值处理与程序设计数值处理采用外点法用均匀网格对求解区域进行离散化，得到的网格系统如图2所示。

一共使用了0~N-1共N 个节点。

节点间距δx 为:图2 墙壁内的网格划分此例中墙壁导热系数为常值，无源项。

则可采用有限体积法对控制方程离散化，得到离散方程为：p p E E W W a T a T a T b =++ （2a ） 式中：P W E P a a a a ++= （2b ）x a E δλ=，x a W δλ=，τδρ∆=x c a P 0（2c ） 00p p b a T = （2d ）其中的上标“0”表示此为上一时刻的值,分别为节点所在控制容积左右边界上的导热系数，由于墙壁导热系数不变，故都等于λ，△τ为时间步长。

由元体能量平衡法可以得知左右边界节点的离散方程分别为： 左边界节点：（3）右边界节点：（4）离散方程的详细推导过程见附录。

程序设计由物理模型可以知道本问题为一维导热问题，一维导热问题的离散方程在取遍所有节点后形成的是三对角的代数方程组，采用追赶法进行求解。

程序构成和方法:程序由主程序和一个子程序构成。

主程序进行变量定义和各已知参数的输入，以及左右边界节点和内部节点控制方程的输入；子程序tdma实现追赶法用来计算每个节点新的温度。

Thomas算法求解过程分为两步:消元和回代。

消元是从系数矩阵的第二行起，逐一将每一行的非零元素消去一个，使原来的三元方程化为二元方程。

消元进行到最后一行时，二元方程就化为一元方程，直接得到最后一个未知数的值。

然后逐一往前回代，由各二元方程求出其它未知解。

程序特点：该程序有很强的适应性，一维常物性非稳态平壁导热问题都可以使用此程序，只要适当更改边值条件即可。

还可以进行修改解决非常物性问题。

程序中对输出节点，最大输出量都进行了控制，对计算结果的分析有很大帮助。

而且Thoms算法的优点需要内存小，工作量小，程序设计简单。

程序流程图：首先对变量赋值，然后由初始条件建立初始温度场，接着从左边界，内部节点，到右边界进行迭代，直到满足精度要求为止，最后输出结果，程序结束。

程序流程如下图3。

4、模型与程序验证模型本题简化为厚度为2δ=的一维非稳态模型如图4所示，初始温度为15℃，在其中间建立坐标系，左两边为对流换热，且换热系数相同都为h=25 W/(m2·℃)，且流体温度T f=-10℃对于x≥0，列出其导热微分方程式及定解条件：22(0,0)T Ta x xδττ∂∂=<<>∂∂ （5）a cλρ=0(,0)(0)T x T x δ=≤≤ （6）(,)0x T x x τ=∂=∂ （7）[](,)(,)x T x h T T xδτδτλ∞=∂-=∂ （8）引入过余温度：图3 程序流程图(,)T x T θτ∞=- （9）直接根据公式得到解析解如下：210.(,)exp()cos()n n n n C Fo θητμμηθ∞==-∑ （10） 式中，2,a xFo τηδδ==，系数n C 应该使上述无穷级数在0τ=是满足初始条件，由傅里叶级数理论可得：2sin cos sin nn n n nC μμμμ=+ （11）n μ是超越方程的根，称为特征根。

tan ,1,2,n nBin μμ== (12)其中 h Bi δλ=。

程序验证（1） 由模型可以得到相关信息然后进行编程，同等时间下计算出中心处温度的解析解和数值解进行比较，数据记录在表1。

然后计算出相对误差，作图5，观察数值解与分析解的比较曲线。

由图表中可以发现，平壁中心不同时刻温度值的分析解和数值解相差不是很大，二者吻合的比较好，可以说明所编制的数值解法的程序是正确的。

相对误差先增大后减小，增大的原因是此时温度接近零度，相对误差的基数比较小，所以造成相对误差较大，但2·℃)10℃xT 0图4 一维导热简化模型是此时的绝对误差并不大，在合理范围内，所以除去个别点外，都满足误差小于百分之1。

可以验证所编数值解法的程序是正确的。

（2）空间步长对墙内壁的温度影响如图6及表2。

在程序编写过程中用网格节点数对空间步长进行控制，为了观察空间步长对墙内壁温度的影响，表中选择了三个不同的空间步长，分别为选取51，101，201个网格节点，则相应的空间步长为，，。

根据不同步长时温度的变化曲线可以看出，空间步长对内墙壁的影响不大，当空间步长控制在合理范围时可以忽略空间步长的影响。

表1 分析解与数值解比较时间（h）分析解（℃）数值解（℃）相对误差（%）0151505．00图5 平壁中心不同时刻的数值解和分析解图6 空间步长对墙内壁温度影响表2 空间步长对温度影响数据时间(h) N=51 N=101 N=201（3）时间步长对温度的影响如图7和表3，根据图中曲线可以看出时间步长选择50s，100s，200s时基本重合，对墙内壁温度影响不大。

图7 时间步长对温度的影响表3 时间步长对温度的影响时间(h) T=50（s）时间(h)T=100（s)时间（h） T=200（s）015001515152468图8 墙内壁温度随时间的变化曲线5 计算结果与分析墙内壁温度分析根据题目中要求，计算寒潮入侵多长时间后内墙壁可以感受到外界气温的变化，通过建模，方程离散化，最终通过程序求解方程，得到图8和表4。

由图可以看出，开始阶段，内墙壁温不变，随着时间的进一步深入，内壁温度开始降低，当很长时间后，温度变化基本趋于平缓，直到再次平衡。

根据图8就可以得到墙内壁温度开始发生变化的时间。

表4 墙内壁温度随时间变化数据表时间（h)温度（℃）时间（h)温度（℃）时间（h)温度（℃）1515导热系数λ对墙内壁温度的影响墙的导热系数对内表面的影响，在图9和表5中发现，导热系数对内壁温度影响比较大，λ=时，温度下降的趋势会更快，要比λ=时下降快的多，下降速度更快，更短时间内达到稳态，因此导热系数越大温度扩散越快，导热系数越小则温度变化越慢，需要更长时间到达稳态，但是这时对于要求恒温的空间有好处，受波动影响更小。

表5 导热系数对墙内壁温度的影响时刻（h）λ=λ=λ=01515151515155图9 导热系数对墙内壁温度的影响墙外换热系数h的影响墙外表面传热系数对温度分布的影响，如图10和表6影响不大。

图10 墙外换热系数对温度影响表6 墙外换热系数对温度影响时刻（h） h=20W/(m2·℃)h=35W/(m2·℃)h=50W/(m2·℃)015151510墙内壁传热热系数的影响由图11可以看出，墙内表面传热系数对内表面温度影响较大，当传热系数比较小时受墙外流体温度影响明显，传热系数越大受墙外流体温度影响越小，当到了极限大的情况下，内表面温度则等于墙内流体温度，不再受墙外流体温度影响。

图11 墙内表面传热系数对温度影响表7 墙内表面传热系数对温度影响时刻（h）h=2W/(m2·℃)h=6W/(m2·℃)h=10W/(m2·℃)墙壁厚度δ对内壁温度的影响改变墙壁的厚度δ，内壁温度将发生变化。

在表8和图12中可以发现，不同厚度的墙壁对外界温度的感应快慢是不一样的，随着墙壁厚度的增加，感应到外界温度变化的时间越来越长，且温度变化越来越慢，墙壁越厚要越长的时间才能达到新的平衡图12 墙壁厚度对内壁温度的影响表8 墙壁厚度对内壁温度的影响时间（h）015151515106 结论(1)整个墙壁经历了以下变化过程：首先外壁直接与室外冷空气接触，温度变化很快，随着时间的推移，墙内各点的温度也开始变化，并影响到右边内墙壁的温度，慢慢降低。

(2)对于墙内壁的温度随时间的变化，变化趋势总是由快到慢，最后重新达到稳态。

当改变墙壁厚度的时候，随着墙壁厚度的增加，对于外界温度的感应是越来越慢。

这对于一些对温度有要求的地反很重要，根据温度的变化作出适当的调整措施。

参考文献[1] 黄善波 刘中良编著 计算传热学基础[2] 杨世铭 陶文铨编著 传热学（第四版）北京高等教育出版社 [3] 王元明编 数学物理方程与特殊函数（第三版） 高等教育出版社附录1数学模型的离散过程推导()()()()221111p p p p p pW E P e w P e e w w x x x x x x x x x δδδδ⎡⎤⎡⎤⎡⎤⎡⎤T T ∂T ∂T ∂T ⎛⎫⎛⎫=-=-+T +⎢⎥⎢⎥⎢⎥ ⎪ ⎪⎢⎥∂∆∂∂∆⎝⎭⎝⎭⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦（1） 按照Taylor 级数展开法，温度对时间的偏导有向前差分格式，中心差分格式和向 后差分格式，使用向后差分格式。

向后差分格式:1pp p P P Pττ-T -T ∂T ⎛⎫= ⎪∂∂⎝⎭ （2）所以()()()()111p p p pp WP P EP ee w w a x x x x x τδδδδ-⎡⎡⎤⎤T T -T T =-+T +⎢⎢⎥⎥∆∆⎢⎢⎥⎥⎣⎣⎦⎦ （3）()()()()()111p p p p p W E P P P e e ww x a x x x x τδδδδ-⎡⎡⎤⎤T T ∆T -T =-+T +⎢⎢⎥⎥∆⎢⎢⎥⎥⎣⎣⎦⎦（4）由x x δ∆= 得02PE W P x c x x x x δααρδτδδδτ∂⎛⎫+T =T +T +T ⎪∆∆⎝⎭（5） P P E E W W a a a b T =T +T + （6）P E W P a a a a =++ E a x αδ=，W a x αδ=，0P c x a ρδτ=∆ 00P P b a T = （7） 边界条件处理：右边界根据元体能量平衡法处理：图1右边界的能量守恒如图1： 2B N s E -Φ+Φ=∆ （8） 其中1()B r fr N h T T -Φ=-，221N N N xλδ----T ΦT = （9）()011211()N N N N r fr N c x h T T x ρλδτ-----∆T -T T -T +-=∆ （10） λ为常数Δxa c λρ= , 2xx δ∆=（11）将式（12-b ）(12-c )入式（12-a ）以得到:012122r N N r fr N c x c x h h T T x x ρδλλρδτδδτ---⎛⎫++T =T ++ ⎪∆∆⎝⎭（12）同理可以得到左边界点离散方程:010()22l l fl c x c x h T hT T T x x λρδλρδδτδτ++=++∆∆ （13） 因此，离散方程为:P P E E W W a T a T a T b =++ （14）式中，0P E W P a a a a =++ （15）E a x αδ=，W a x αδ=，0P c x a ρδτ=∆（16） 00P P b a T = （17）式中，上标“0”表示上一时刻值，α为热扩散系数，τ∆为时间步长。