P1
P8
P7
F
6P2 B 0C E
5
P7 P4
P8
P3
P1
65
P2 P0
P3 C P4
Ji kiikijkik kim Ti Jj kkk k T jijjjkjm j
Tk kkkkkJkikjkkkm
Jm kmi kmi kmk kmm Tm
Ti （2） t
niinijnik nim T pi
j nnn n t pj jijjjkjm
nkinkjnkk nkm Tk pk
t pm ， nmi nmi nmk nmm
Tm t
式中，Ji、Jj、Jk、Jm分别表示该矩形单元对其四个
节点i、j、k、m所在的控制体积能量方程的贡献；Ti、Tj、Tk、Tm分别表示四个节点i、j、k、m的温 度.系数矩阵和列向量由表1计算.
表1 矩阵系数对比
Tab.1 comparison of matrix coefficients
矩阵及向量参数
有限单元法
新有限 体积法×3/4
常用有限 体积法 3 22
8A b c
（a） FVM1 （b） FVM2
图1 控制体积示意 Fig.1 control volume schematic
图1中P表示节点，单元的划分用实线表示，围
绕节点P0的控制体积用虚线表示.图1（a）、图1（b）中的控制体积均为多边形ABCDEFGH，其中A、C、E、G为其所在矩形单元的中心，不同之处在于图1（a）中的B、D、F、H分别为所在线段的中点；图1（b）中B、D、F、H至P0距离为所在边长的三分之二.
kii kjj kkk kmm
kij kji kkm kmk
kik kjm kki kmj kim kjk kkj kmi
3 A 2b2 c2 6 A b2 c2 6 A 2c2 b2 6
A b2 c2
3
A b2 c2
A 2b2 c2 6
A 3b2 c2 8
A b2 c2 6
nii njj nkk nmm
nik njm nki nmj
3 有限体积法的单元分析
在有限单元法中，采用温度场离散，将每个方
程分解为各个单元对它的贡献.这样就可采用单元搜索方法，计算每个单元对各个方程的贡献，最后总体合成为要求解的线性方程组.对于有限体积法，
11
CA CA
1636
nij njk nkm nmi 1 453
CA CA CA
16nji nkj nmk nim86418
111
pi pj pk pm qvAqvAqvA
444注：b、c分别为沿Y方向、X方向的单元的边长，A为四边形单元的面积，A bc.
1
CA9
A 2c2 b2 6
111
CA864 15
CA864
A b2 c2 8
A 3c2 b2 8
9
CA16
第5期 秦跃平，等：非稳态导热问题有限体积法
579
4 边界分析
在导热问题中，边界条件分为三类：第一类边界上的温度为已知值，第二类边界上边界内法线方向的热流密度分量q2为已知，第三类边界条件上相接触的流体温度Tf和换热系数α为已知.
由于第一类边界上的温度为已知值，不必求解，因而不必建立围绕该节点的控制体积的能量方程.
图2为边界节点及其控制体积圈划的示意图.可以看出，对于边界节点，除了单元对能量方程有贡献外，边界上的热流也有贡献.在单元的离散过程中，如果不区分边界节点和内部节点，则图2中各个单元对边界节点的贡献已经计算.以下只需补充第二类和第三类边界上热流对能量方程的贡献.
由表1~表3可以看出，对于稳态导热问题，本文提出的有限体积法与有限单元法形式完全相同，只是对应的所有系数均相差3/4倍.而反映非稳态的系数矩阵（[n]）的各系数没有必然联系，但有限单元法在处理非稳态问题时时间变量采用差分，而不是变分原理，理论上并不完备.常用的有限体积法与有限单元法没有必然联系，其原因在于不同的控制体积圈划方法造成了内热源、单元中的热流、边界上的热流对能量方程贡献的比例有差异，这种差异会造成计算结果明显不同.为了对比各种计算方法的优劣，说明控制体积圈划不当造成计算结果错误的严重性，以下以有内热源无限长方柱体非稳态温度场为例，进行理论计算和数值模拟计算.
B
B 有一截面为i
G
j
i
G
j
5 计算实例的理论解
（a）FVM1 （b） FVM2
图2 控制体积边界节点的示意 Fig.2 control volume schematic of boundary nodes
下面分析边界ij对两节点i和j周围控制体积
能量方程的贡献，
vi Ji uii uij Ti
. （3）
J uu jj j vj ji Tj
通过分析可得：对于第二类边界条件，式（3）中
各参数取值见表2.表2中，q2为边界热流密度，sb为边界ij长度.对于第三类边界条件，式（3）中各参数取值见表3.表3中，h为传热系数，Tf为介质温度.
表2 第二类边界矩阵系数
Tab.2 matrix coefficients of second boundary
矩阵及向量参数
有限单元法
新有限体积法
×3/4
1 q2sb2
0.2 m×0.3 m的无限长方柱体，初
始温度为100 ℃，放入水温为10 ℃的液体中冷却，其中该方柱体的体积密度ρ=1 500 kg/m3，比热容cp=1 000 J/(kg·K)，导热系数λ=40 W/(m·K)，方柱体表面和水之间的表面传热系数h=80 W/(m2·K)，分析方柱体的温度分布随时间的变化情况.
本例题可以近似的看作是一个第三类边界条件下的二维非稳态导热问题.根据热力学理论[11]，当τ＞168.8（即Fo＞0.2）时，计算得方柱体的温度分布为
t(x,y, ) 96.977753 cos(4.328x)
（4）
cos(3.478667y)exp(-0.000822 )+10.
6 编制数值计算软件
6.1 网格划分
由于该问题的对称性，取四分之一截面作为数值计算的区域，将其进行矩形网格剖分，见图3.
uii ujj uij=ujivi vj
常用有限 体积法
0 1
q2sb2
1
q2sb2
表3 第三类边界矩阵系数
Tab.3 matrix coefficients of third boundary
矩阵及向量 参数
有限单元法 1 hsb3
1hsb61 hsbTf2
uii ujj
新有限体积法
×3/4 1hsb3
1hsb6 1
hsbTf2
uij=uji vi vj
常用有限体
积法 1 hsb4
1hsb41 hsbTf2