计算传热学程序

计算传热学课程报告

一、问题概述：

有限单元法是上个世纪五、六十年代首先在力学中发展起来的数值计算方法，由于它是基于变分原理，理论基础统一，对于复杂边界的适应程度比较好，所以很快的在其它领域得到运用，其中就包括了在传热学中的运用。本次计算传热学的课程就是对有限单元法在传热学中运用的一个学习与练习。

有限单元法处理问题的步骤，首先是建立有限元模型也即是将问题离散化，它的主要步骤之一就是将要计算的物体进行有限元的划分；第二步，进行单元分析也就是将变分原理运用到问题的方程与单元中，形成单元刚度矩阵；第三步，进行整体刚度矩阵的组集；最后就是引入边界条件进行求解的过程。

在计算传热学的课程中，主要完成了两个任务：第一，是将一个比较复杂的活塞进行了网格划分，并编译成一个通用性比较好的程序。第二，在前一个程序的基础上，加入计算过程，运用焓法，对一个比较简单的平面相变问题进行了计算。

二、划分单元网格：

划分单元网格是将问题进行有限元法分析的基础，但是如果在图纸上进行手工的单元划分，不但繁琐、容易出错，而且也不利于进一步计算程序的利用。因此有必要编辑一个程序，以自动完成划分网格的目的。网格的自动划分必须遵循以下的几条规则：(1).要严格区分边界单元与内部单元，并且严格区分边界单元不同的组；(2).单元标号必须先标志内部单元，然后依次标志第一类边界条件，第二类边界条件，第三类边界条件，如果同一类边界条件中有不同的组，那么也必须严格先划分第一组，然后第二组，第三组；(3). 对于边界单元，每一个边界单元必须只有一条边在边界上，而且为了程序的简单，一般是j,m边作为边界；(4).内部单元节点标号必须遵循逆时针方向的规则；(5). 一个单元中只能有一种材料组成。

遵循以上的规则，用FORTRAN 90编制了一个对形状比较复杂的活塞的网格划分，由于在编制过程中考虑了多种情况，所以这个程序有比较好的通用性，只需要输入不同的数据，程序也可以对许多其它情况进行划分。

需要指出的是，由于FORTRAN 90程序对于制图功能比较弱，所以下面的图是用VB 6.0的程序做出的，由于该网格划分程序集成了后续对第一类边界条件和第三类边界条件的焓法计算程序，故该程序源代码将在最后统一给出。网格划分的结果如图(1)。

需要输入的初始数据主要有：边界单元分组总数、边界单元分组中前一组的最后一个单元号、各组边界单元节点数、各边界单元边界节点号、

每一条层线的左右端点Ｘ、Ｙ坐标，每一层单元划分所属的类型。

图（１）

三、焓法有限单元法原理：

带有相变的传热问题，又被称为斯蒂芬问题，在冶金、铸造、建筑、冷冻、航天和医疗等领域有着广泛的运用。由于是这类问题存在的相变过程，使得求解区域中存在着一个随着时间移动的固液或者固气界面，这一界面使得这种问题的求解非常困难。

一般说来，处理这种问题有两种不同的思路，一种思路是，首先着眼于相变界面的求解，确定相变界面以后再分别处理固相或者液相的温度分布；另一种思路是将该问题看作是“单相”的非线性导热问题，首先确定整个求解区域上的温度分布或者焓的分布，然后把达到相变温度的位置定为相变界面。在实际运用中，后一种思路比较简单和实用而得到了广泛的使用；目前后一种处理方式中主要有焓法和显热容法，这里主要讨论焓法。

焓法有限元法的主要思路是，不把温度作为求解的变量，而是把焓作为求解变量，因此可以在固相和液相的整个区域建立统一的方程进行求解。然后根据焓与温度的关系确定整个区域的温度与相变的界面。焓法的优点就是不需要跟踪相变的界面而可以对整个区域进行统一的求解。焓法的主要推导如下：

在一般过程中，传热微分方程由(1)式表示：

T k t

T

C

2∇=∂∂ρ （１）

其中k 为导热系数，ρ为物体的密度，C 为物体的定压热容，T 为物体的温度，而t 表示时间。由于焓h 与温度T 存在如下的关系：

C t

h

=∂∂ 所以（１）式可以化为（２）式：

h k t

h

C

2∇=∂∂ρ （２） 该式子即是在整个区域上都适用的焓法的基本方程，对于这个方程，在空间上用有限单元法离散，在时间上用向后差分格式离散，就可以得到如下的(3)式：

t t t t P h t N h t N K }{}{/][}{)/][]([+⋅∆=⋅∆+∆- （３）

其中：][K ——整体刚度矩阵 ][N ——变焓矩阵

t ∆——时间步长

t h }{——t 时刻热焓 t t h ∆-}{——t t ∆-时刻热焓

t P }{——t 时刻右端列向量

值得注意的是，首先使用变分法的使用要直接对焓而不是温度变分；其次整体刚度矩阵和变焓矩阵是与时间无关的的矩阵，而右端向量则与边界条件有关。最后焓与温度的关系由（４）式给出：

s m s m m m s m l CT h l CT h CT cT h C l h T C

h T +≥+≤≤≤⎪⎩

⎪

⎨⎧-=/)(/ （４）

四、焓法实例：

为了简单起见，运用上述的焓法对如图(2)的简单区域的相变问题进行求解，其条件如下：边长为80cm ⨯60cm ，周围用温度为500K ，铁水温度为1833K ，比热为711.62J/kg.k ，潜热为271100J/kg ，密度为7800kg/m^3，

导热系数为33.5W/m. 相温度为1790K，使用第三类边界条件，换热系数=

α1950W/K·m2。

图（2）

图中1-4，1-15等边都是第三类边界，其余单元则作为边界单元处理。程序源代码如下：

PROGRAM MAIN

INTEGER

L0,V0,E3,C0,B0,F1,V7,M2,A0,D9,Z0,PA NBIE,B33

80

integer,dimension(:),allocatable::B(:),F(:), W(:),M(:),J(:),I(:),H10(:),H11(:),H12(:),H1 3(:)

90

real,dimension(:),allocatable::P(:),Q(:),X(:) ,Z(:),Y(:),HH1(:),HH2(:),HH3(:),HH32(:), HH4(:),TT(:)

DOUBLE

PRECISION,DIMENSION(:),ALLOCATA BLE::H9(:) double precision K1(3,3),N1(3,3) double

precision,dimension(:,:),allocatable::K2(:,:) ,N2(:,:),EH(:,:),EH2(:,:)

100

REAL,dimension(:),allocatable::C(:),S(:),R (:)

REAL

KC,CP,MD,TS1 ,TS2,ARF,TF,TM

REAL HH33(3)

open (1,file='score.dat')

70 read(1,*) L0,V0,E3,C0,B0,F1 75 V7=V0+1

allocate(B(F1),F(B0),W(B0),H9(E3),M(V7