计算传热学程序设计

中国石油大学(华东)

储运与建筑工程学院热能与动力工程系

《计算传热学程序设计》

设计报告

学生姓名:

学号：

专业班级：

指导教师

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 c

x x ρλτ∂∂∂=∂∂∂ （1a ）

初始条件：

（1b ）

在两侧相应的边界条件是第三类边界条件，分别由傅立叶定律可描述如下： 左边界：

020

2()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 T

a 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）

直接根据公式得到解析解如下：

2

1

0.(,)exp()cos()n n n n C Fo θητμμηθ∞==-∑ （10） 式中，2

,a x

Fo τ

ηδ

δ

=

=，系数n C 应该使上述无穷级数在0τ=是满足初始条件，由傅里叶

级数理论可得：

2sin cos sin n

n n n n

C μμμμ=

+ （11）

n μ是超越方程的根，称为特征根。

tan ,1,2,n n

Bi

n μμ=

= (12)

其中 h Bi δ

λ

=。

程序验证

（1） 由模型可以得到相关信息然后进行编程，同等时间下计算出中心处温度的解析解和数值解进行比较，数据记录在表1。然后计算出相对误差，作图5，观察数值解与分析解的比较曲线。

由图表中可以发现，平壁中心不同时刻温度值的分析解和数值解相差不是很大，二者吻合的比较好，可以说明所编制的数值解法的程序是正确的。相对误差先增大后减小，增大的原因是此时温度接近零度，相对误差的基数比较小，所以造成相对误差较大，但

2·℃)

10℃

x

T 0图4 一维导热简化模型