计算传热学程序
计算传热学课程报告
一、问题概述：
有限单元法是上个世纪五、六十年代首先在力学中发展起来的数值计算方法，由于它是基于变分原理，理论基础统一，对于复杂边界的适应程度比较好，所以很快的在其它领域得到运用，其中就包括了在传热学中的运用。本次计算传热学的课程就是对有限单元法在传热学中运用的一个学习与练习。
有限单元法处理问题的步骤，首先是建立有限元模型也即是将问题离散化，它的主要步骤之一就是将要计算的物体进行有限元的划分；第二步，进行单元分析也就是将变分原理运用到问题的方程与单元中，形成单元刚度矩阵；第三步，进行整体刚度矩阵的组集；最后就是引入边界条件进行求解的过程。
在计算传热学的课程中，主要完成了两个任务：第一，是将一个比较复杂的活塞进行了网格划分，并编译成一个通用性比较好的程序。第二，在前一个程序的基础上，加入计算过程，运用焓法，对一个比较简单的平面相变问题进行了计算。
二、划分单元网格：
划分单元网格是将问题进行有限元法分析的基础，但是如果在图纸上进行手工的单元划分，不但繁琐、容易出错，而且也不利于进一步计算程序的利用。因此有必要编辑一个程序，以自动完成划分网格的目的。网格的自动划分必须遵循以下的几条规则：(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
),H10(E3),H11(E3),H12(E3),H13(E3),HH1 (L0),HH2(L0),HH3(L0),HH32(L0),HH4(L 0),TT(L0))
allocate(J(V0),P(V0),Q(V0),Y(V0),Z(V0),I (C0),C(C0),S(C0),X(L0),R(L0),k2(L0,L0), n2(L0,L0),EH(L0,L0))
110 DO 125 J0=1,B0
120 READ(1,*) W(J0),F(J0)
125 CONTINUE
READ(1,*)(B(A0),A0=1,F1) 150 READ(1,*) (M(A0),A0=1,V7) 160 do 190 A0=1,V0
165 READ(1,*) J(A0),P(A0),Q(A0),Y(A0),Z(A0)
190 CONTINUE
195 DO 215 A0=1,C0
200 READ(1,*) I(A0),C(A0),S(A0) 215 CONTINUE
READ(1,*)KC,CP,MD,TS,PANBIE,ARF,T F,TS1,TM
READ(1,*)(HH1(I0),I0=1,L0)
B33=W(1)
DO I0=1,L0
TT(I0)=0
ENDDO
DO I0=1,L0
DO J0=1,L0
K2(I0,J0)=0
N2(I0,J0)=0
ENDDO
ENDDO
DO I0=1,L0
HH3(I0)=0
HH4(I0)=0
ENDDO
TS2=0
M2=1
250 DO I0=1,V0
252 IF (J(I0)==5) GOTO 266
254 IF(J(I0)==6) GOTO 262
256 IF(J(I0)==7) GOTO 262
258 D9=M(I0+1)-1
GO TO 270
262 D9=M(I0+1)-5
GO TO 270
266 D9=M(I0+1)-4
270 DO J0=M(I0),D9
272
X(J0)=Y(I0)+(J0-M(I0))*(Z(I0)-Y(I0))/(D 9-M(I0))
274
R(J0)=P(I0)+(J0-M(I0))*(Q(I0)-P(I0))/(D9 -M(I0))
276 ENDDO
ENDDO
290 Z0=0
292 DO I0=1,V0-1
294 D9=M(I0+1)-M(I0)
296 IF(J(I0)/=8) GOTO 302
298 CALL JK1(D9,M,I0,V7,Z0,H9,E3,M2,B0,F,F1,B, W,H10,H11,H12,H13)
GOTO 316
302 IF(J(I0)==2) GOTO 314
304 IF(J(I0)==4) GOTO 314
306 IF(J(I0)==10) GOTO 314
308 IF(J(I0)==12) GOTO 314
310 CALL JK7(D9,M,I0,V7,B0,Z0,H9,E3,M2,F,F1,B, W,J,V0,H10,H11,H12,H13)
GOTO 316
314CALL JK4(D9,M,I0,V7,Z0,H9,E3,M2,B0,F,F1,B, W,J,V0,H10,H11,H12,H13)
316 ENDDO
350 DO A0=1,C0
352 X(I(A0))=C(A0)
R(I(A0))=S(A0)
356 ENDDO
close(1)
open (2,file='out.dat')
380 write(2,"('coordinate of knots')")
385 DO A0=1,L0
390
write(2,500)A0,X(A0),A0,R(A0)
400 ENDDO
405 write(2,"('INFORMATION OF ELEMENTS')")
410 DO I0=1,E3
415
write(2,600)I0,H10(I0),H11(I0),H12(I0),H 13(I0)
420 ENDDO
close(2)
open (3,file='out2.dat')
do a0=1,L0
WRITE(3,*)X(A0)
WRITE(3,*)R(A0)
ENDDO
DO I0=1,E3
WRITE
(3,700)H10(I0),H11(I0),H12(I0),H13(I0)
ENDDO
CLOSE(3)
500
format(1x,'X',I4,'=',F16.4,'R',I4,'=',F16.4) 600 FORMA T(1X,'H',I4,'=',I4,I4,I4,I4)
700 FORMA T(1X,I4,I4,I4,I4)
open (4,file='out4.dat')
DO J0=1,E3
call
DANGANG(K1,X,R,H10,H11,H12,J0,L0, KC,E3,CP,N1,PANBIE,ARF,HH33,TF,W, B0,MD,B33)
CALL
ZHENGGANG(K1,K2,H10,H11,H12,L0,E 3,J0,N1,N2,HH33,HH4)
ENDDO
CALL ZUIHOU(K2,N2,EH,L0)
DO WHILE(TS2<=TS1)
CALL
YOUBIAN(N2,L0,HH1,HH3,PANBIE,HH 4)
IF(PANBIE==1)THEN
CALL
ZHENGLI(L0,EH,W,B0,H10,H11,H12,E3, EH2,HH32,HH3,HH1,B33)
CALL
GUASS(EH2,HH32,HH2,L0)
ELSE
CALL
GUASS(EH,HH3,HH2,L0)
ENDIF
CALL
WENDU(HH2,L0,CP,TM,TT)
IF
(TT(15)>=500.AND.(TT(7)<=1823)) THEN
DO J0=1,L0
WRITE(4,800)J0,TT(J0)
ENDDO
ELSE
CONTINUE
ENDIF
DO J0=1,L0
HH1(J0)=HH2(J0)
ENDDO
TS2=TS2+TS
ENDDO
CLOSE(4)
800 FORMA T(1X,'I=',I4,'H=',F16.4) 425 END
SUBROUTINE
WENDU(HH2,L0,CP,TM,TT)
INTEGER L0
REAL
HH2(L0),TT(L0),CP,LS,C,TM
LS=271100
C=CP*TM
DO I=1,L0
IF(HH2(I)<=C) THEN
TT(I)=HH2(I)/CP
ELSE IF(HH2(I)>=C+LS) THEN
TT(I)=(HH2(I)-LS)/CP
ELSE
TT(I)=TM
ENDIF
ENDDO
END
SUBROUTINE
GUASS(EH2,HH32,HH2,L0)
INTEGER L0,L
REAL
HH32(L0),HH2(L0),SU,C(L0)
double precision EH2(L0,L0),D(L0,L0),O
L=0
SU=0.0
DO K=1,L0-1
M=K
S=EH2(K,K)
L=K
DO WHILE(M<=L0-1)
M=M+1
S1=ABS(S)
S2=EH2(M,K)
S3=ABS(S2)
IF(S1<=S3)THEN
L=M
S=EH2(M,K)
ELSE
CONTINUE
ENDIF
ENDDO
M=0
DO J=K,L0
C(J)=EH2(L,J)
EH2(L,J)=EH2(K,J)
EH2(K,J)=C(J)
ENDDO
O=HH32(L)
HH32(L)=HH32(K)
HH32(K)=O
DO I=K+1,L0
D(I,K)=EH2(I,K)/EH2(K,K)
ENDDO
DO I=K+1,L0
DO J=K+1,L0 EH2(I,J)=EH2(I,J)-D(I,K)*EH2(K,J)
ENDDO
ENDDO
DO I=K+1,L0 HH32(I)=HH32(I)-D(I,K)*HH32(K)
ENDDO
ENDDO
S4=EH2(L0,L0)
HH2(L0)=HH32(L0)/S4
DO
K=L0-1,1,-1
DO
J=K+1,L0
SU=SU+EH2(K,J)*HH2(J)
ENDDO HH2(K)=(HH32(K)-SU)/EH2(K,K)
SU=0
ENDDO
END SUBROUTINE
ZHENGLI(L0,EH,W,B0,H10,H11,H12,E3, EH2,HH32,HH3,HH1,B33)
INTEGER
L0,B0,W(B0),B33,E3,H10(E3),H11(E3),H 12(E3),J33,M33
REAL HH32(L0),HH1(L0),HH3(L0)
double precision EH(L0,L0),EH2(L0,L0)
DO I0=1,L0
DO J0=1,L0
EH2(I0,J0)=EH(I0,J0)
ENDDO
HH32(I0)=HH3(I0)
ENDDO
DO I0=1,E3
IF(I0<=B33) THEN
CONTINUE
ELSE
J33=H11(I0)
M33=H12(I0)
EH2(J33,J33)=EH(J33,J33)*(10**8) EH2(M33,M33)=EH(M33,M33)*(10**8)
HH32(J33)=
EH2(J33,J33)*HH1(J33)
HH32(M33)=EH2(M33,M33)*HH1(M33)
ENDIF
ENDDO
END
SUBROUTINE
YOUBIAN(N2,L0,HH1,HH3,PANBIE,HH 4)
INTEGER L0,PANBIE
double precision X,HH(L0)
double precision N2(L0,L0)
REAL
HH1(L0),HH3(L0),HH4(L0)
I=0
X=10
DO I=1,L0
HH3(I)=0
HH(I)=0
ENDDO
DO I0=1,L0
DO J0=1,L0
HH(I0)=HH(I0)+N2(I0,J0)*HH1(J0)
ENDDO
ENDDO
DO I0=1,L0
HH3(I0)=HH(I0)/X
ENDDO
DO I0=1,L0
HH3(I0)=HH3(I0)+HH4(I0)
ENDDO
END
SUBROUTINE
ZUIHOU(K2,N2,EH,L0)
INTEGER L0
double precision
K2(L0,L0),N2(L0,L0),EH(L0,L0),EHX(L0 ,L0)
REAL X0
X0=10
DO I0=1,L0
DO J0=1,L0
EHX(I0,J0)=0
ENDDO
ENDDO
DO I0=1,L0
DO J0=1,L0
EHX(I0,J0)=N2(I0,J0)/X0
ENDDO
ENDDO
DO I0=1,L0
DO J0=1,L0
EH(I0,J0)=EHX(I0,J0)+K2(I0,J0)
ENDDO
ENDDO
END
SUBROUTINE
ZHENGGANG(K1,K2,H10,H11,H12,L0,E 3,J0,N1,N2,HH33,HH4)
INTEGER J0,I33,J33,M33,L0,E3
INTEGER
H10(E3),H11(E3),H12(E3)
double precision K1(3,3),K2(L0,L0),N1(3,3),N2(L0,L0)
REAL HH4(L0),HH33(3)
I33=H10(J0)
J33=H11(J0)
M33=H12(J0)
K2(I33,I33)=K2(I33,I33)+K1(1,1)
K2(I33,J33)=K2(I33,J33)+K1(1,2) K2(I33,M33)=K2(I33,M33)+K1(1,3)
K2(J33,I33)=K2(J33,I33)+K1(2,1)
K2(J33,J33)=K2(J33,J33)+K1(2,2)
K2(J33,M33)=K2(J33,M33)+K1(2,3)
K2(M33,I33)=K2(M33,I33)+K1(3,1)
K2(M33,J33)=K2(M33,J33)+K1(3,2)
K2(M33,M33)=K2(M33,M33)+K1(3,3)
N2(I33,I33)=N2(I33,I33)+N1(1,1)
N2(I33,J33)=N2(I33,J33)+N1(1,2)
N2(I33,M33)=N2(I33,M33)+N1(1,3)
N2(J33,I33)=N2(J33,I33)+N1(2,1) N2(J33,J33)=N2(J33,J33)+N1(2,2)
N2(J33,M33)=N2(J33,M33)+N1(2,3)
N2(M33,I33)=N2(M33,I33)+N1(3,1)
N2(M33,J33)=N2(M33,J33)+N1(3,2)
N2(M33,M33)=N2(M33,M33)+N1(3,3)
HH4(J33)=HH33(2)+HH4(J33)
HH4(M33)=HH33(3)+HH4(M33)
END
SUBROUTINE
DANGANG(K1,X,R,H10,H11,H12,J0,L0, KC,E3,CP,N1,PANBIE,ARF,HH33,TF,W, B0,MD,B33)
INTEGER
J0,I33,J33,M33,L0,E3,PANBIE,B0,B33
REAL
B1,B2,B3,C1,C2,C3,SM,KC,KCF,CP,SI,A RF,MD
INTEGER
H10(E3),H11(E3),H12(E3),W(B0)
REAL X(L0),R(L0),HH33(3)
double precision K1(3,3),N1(3,3)
I33=H10(J0)
J33=H11(J0)
M33=H12(J0)
B1=R(J33)-R(M33)
B2=R(M33)-R(I33)
B3=R(I33)-R(J33)
C1=X(M33)-X(J33)
C2=X(I33)-X(M33)
C3=X(J33)-X(I33)
SM=0.5*(B1*C2-B2*C1)
KCF=KC/(4*SM)
SI=SQRT(B1**2+C1**2)
IF((J0<=B33).and.(J0/=1).AND.(J0/=3).A ND.(J0/=5))THEN
K1(1,1)=KCF*(B1**2+C1**2)
K1(2,2)=KCF*(B2**2+C2**2)
K1(3,3)=KCF*(B3**2+C3**2)
K1(1,2)=KCF*(B1*B2+C1*C2)
K1(2,1)=K1(1,2)
K1(1,3)=KCF*(B1*B3+C1*C3)
K1(3,1)=K1(1,3)
K1(2,3)=KCF*(B2*B3+C2*C3)
K1(3,2)=K1(2,3)
N1(1,1)=(MD*CP*SM)/6.0
N1(2,2)=N1(1,1)
N1(3,3)=N1(1,1)
N1(1,2)=(MD*CP*SM)/12.0
N1(1,3)=N1(1,2)
N1(2,1)=N1(1,2)
N1(2,3)=N1(1,2)
N1(3,1)=N1(1,2)
N1(3,2)=N1(1,2)
HH33(1)=0
HH33(2)=0
HH33(3)=0
ELSE
K1(1,1)=KCF*(B1**2+C1**2)
K1(2,2)=KCF*(B2**2+C2**2)+(ARF*SI) /3.0
K1(3,3)=KCF*(B3**2+C3**2)+(ARF*SI) /3.0
K1(1,2)=KCF*(B1*B2+C1*C2)
K1(2,1)=K1(1,2)
K1(2,3)=KCF*(B2*B3+C2*C3)+(ARF*SI )/6.0
K1(3,2)=K1(2,3)
K1(1,3)=KCF*(B1*B3+C1*C3)
K1(3,1)=K1(1,3)
N1(1,1)=(MD*CP*SM)/6.0
N1(2,2)=N1(1,1)
N1(3,3)=N1(1,1)
N1(1,2)=(MD*CP*SM)/12.0
N1(1,3)=N1(1,2)
N1(2,1)=N1(1,2)
N1(2,3)=N1(1,2)
N1(3,1)=N1(1,2)
N1(3,2)=N1(1,2)
HH33(2)=(ARF*SI*TF)/2
HH33(3)=HH33(2)
HH33(1)=0
ENDIF
END
SUBROUTINE
JK1(D9,M,I0,V7,Z0,H9,E3,M2,B0,F,F1,B, W,H10,H11,H12,H13)
INTEGER
D9,I0,Z0,V7,J0,Z9,I1,M1,J1,E3,M2,B0,F1 DOUBLE PRECISION H9(E3)
INTEGER
M(V7),F(B0),B(F1),W(B0),H10(E3),H11( E3),H12(E3),H13(E3)
3802 D9=(D9+1)/2
3804 DO J0=2,D9
3806 I1=M(I0)+J0-1
J1=M(I0+1)+J0-2
3810 M1=M(I0)+J0-2
3812 IF(J0==2)GOTO 3822
3814 Z9=Z0
CALL
CM(H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12, H13)
Z0=Z0+1
3820 GOTO 3824
3822 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, H10,H11,H12,H13)
3824 I1=M(I0)+J0-1
J1=M(I0+1)+J0-1
3828 M1=M(I0+1)+J0-2
3830 CALL BCE(B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B ,W,H10,H11,H12,H13)
3832 ENDDO
3834 DO J0=D9,2*D9-2
3836 I1=M(I0)+J0-1
J1=M(I0)+J0
3840 M1=M(I0+1)+J0
3842 IF(J0==2*D9-2)GOTO 3852 3844 Z9=Z0
CALL
CM(H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12, H13)
Z0=Z0+1
3850 GOTO 3854
3852 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, H10,H11,H12,H13)
3854 I1=M(I0)+J0-1
J1=M(I0+1)+J0
3858 M1=M(I0+1)+J0-1
CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, H10,H11,H12,H13)
3862 ENDDO
END
subroutine CM (H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12,H1 3)
INTEGER E3,Z9,I1,J1,M1,M2
DOUBLE PRECISION H9(E3)
INTEGER
H10(E3),H11(E3),H12(E3),H13(E3)
3805
H9(Z9+1)=I1*1E-3+J1*1E-6+M1*1E-9
H10(Z9+1)=I1
H11(Z9+1)=J1
H12(Z9+1)=M1
H13(Z9+1)=M2
3810 END
SUBROUTINE BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, H10,H11,H12,H13)
INTEGER
B0,Z9,Z0,A0,Y8,S2,I1,J1,M1,M2,E3,F1 DOUBLE PRECISION H9(E3)
INTEGER F(B0),B(F1),W(B0)
INTEGER
H10(E3),H11(E3),H12(E3),H13(E3)
3702 Y8=0
S2=0
3704 DO A0=1,B0
3708 CALL MJ(S2,A0,F,B0,Y8,B,F1,H9,E3,I1,J1,M1, M2,W,H10,H11,H12,H13)
3710 IF(S2==2) GOTO 3720
3712 ENDDO
3714 Z9=Z0
CALL CM (H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12,H1 3)
Z0=Z0+1
3720 END
SUBROUTINE
MJ(S2,A0,F,B0,Y8,B,F1,H9,E3,I1,J1,M1, M2,W,H10,H11,H12,H13)
INTEGER
S2,E9,B0,A0,Y8,F1,Z9,E3,I1,J1,M1,M2
INTEGER F(B0),B(F1),W(B0)
DOUBLE PRECISION H9(E3)
INTEGER
H10(E3),H11(E3),H12(E3),H13(E3)
3602 S2=0
3604 DO E9=1,F(A0)
3608 Y8=Y8+1
3610 IF (B(Y8)==J1) THEN
S2=S2+1
ENDIF
IF(B(Y8)==M1) THEN
S2=S2+1
ENDIF
ENDDO
3620 IF(S2/=2)GOTO 3628
3622 Z9=W(A0)
CALL CM (H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12,H1 3)
3626 W(A0)=W(A0)+1
3628 END
subroutine JK7 (D9,M,I0,V7,B0,Z0,H9,E3,M2,F,F1,B,W,J, V0,H10,H11,H12,H13)
INTEGER
D9,I0,Z0,V7,J0,Z9,I1,M1,J1,E3,M2,B0,F1, V0
DOUBLE PRECISION H9(E3)
INTEGER
M(V7),F(B0),B(F1),W(B0),J(V0)
INTEGER
H10(E3),H11(E3),H12(E3),H13(E3)
4502 I1=M(I0+1)+1
J1=M(I0+1)
4506 M1=M(I0)
CALL
BCE(B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B ,W,H10,H11,H12,H13)
4510 IF(J(I0)==5) D9=D9-5
4512 IF(J(I0)==6) D9=D9-5
4518 IF(J(I0)/=7) GOTO 4522 4520 D9=D9-4
4522 DO J0=2,D9
4524 IF(J(I0)==1) GOTO 4532
4526 IF(J(I0)==7) GOTO 4532
4528 IF(J(I0)==9) GOTO 4532 4530 GOTO 4534
4532 IF(J0==D9) GOTO 4578
4534 I1=M(I0+1)+J0-1
J1=M(I0)+J0-2
4538 M1=M(I0)+J0-1
4540 IF(J(I0)==9) GOTO 4552
4542 IF(J(I0)==11) GOTO 4552 4544 Z9=Z0
CALL
CM(H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12, H13)
Z0=Z0+1
4550 GOTO 4554
4552 CALL BCE(B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B ,W,H10,H11,H12,H13)
4554 IF(J0 4556 GOTO 4578 4558 I1=M(I0)+J0-1 J1=M(I0+1)+J0 4562 M1=M(I0+1)+J0-1 4564 IF(J(I0)==9) GOTO 4576 4566 IF(J(I0)==11) GOTO 4576 4568 Z9=Z0 CALL CM (H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12,H1 3) Z0=Z0+1 4574 GOTO 4578 4576 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, H10,H11,H12,H13) 4578 ENDDO 4580 IF(J(I0)==1) GOTO 4586 4582 IF(J(I0)==9) GOTO 4586 4584 GOTO 4594 4586 I1=M(I0)+D9-2 J1=M(I0)+D9-1 4590 M1=M(I0+1)+D9-1 CALL BCE(B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B ,W,H10,H11,H12,H13) 4594 IF(J(I0)==3) GOTO 4600 4596 IF(J(I0)==11) GOTO 4600 4598 GOTO 4608 4600 I1=M(I0+1)+D9-1 J1=M(I0)+D9-1 4604 M1=M(I0+1)+D9 CALL BCE(B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B ,W,H10,H11,H12,H13) 4608 IF(J(I0)/=5) GOTO 4614 4610 CALL EAB(B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1, B,W,M,I0,J,V0,V7,H10,H11,H12,H13) 4612 CALL EDEC(B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1, B,W,M,I0,V7,H10,H11,H12,H13) 4614 IF(J(I0)==6) GOTO 4620 4616 IF(J(I0)==7) GOTO 4620 4618 GOTO 4624 4620 CALL EAB (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, M,I0,J,V0,V7,H10,H11,H12,H13) 4622 CALL EFGHI (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, M,I0,V7,H10,H11,H12,H13) 4624 END SUBROUTINE EAB (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, M,I0,J,V0,V7,H10,H11,H12,H13) INTEGER B0,Z9,Z0,I1,J1,M1,M2,E3,F1,I0,V0,V7 DOUBLE PRECISION H9(E3) INTEGER F(B0),B(F1),W(B0),J(V0),M(V7) INTEGER H10(E3),H11(E3),H12(E3),H13(E3) 4002 I1=M(I0+1)-6 4004 IF(J(I0)==5) J1=M(I0+1)-3 4008 IF(J(I0)/=5) J1=M(I0+1)-4 4012 IF(J(I0+1)<5) GOTO 4020 4014 IF(J(I0+1)>7) GOTO 4020 4016 M1=M(I0+2)-5 GOTO 4022 4020 M1=M(I0+2)-2 4022 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, H10,H11,H12,H13) 4024 I1=M(I0+1)-6 J1=M(I0+1)-5 4028 IF(J(I0)==5) M1=M(I0+1)-3 4032 IF(J(I0)/=5) M1=M(I0+1)-4 4036 Z9=Z0 CALL CM(H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12, H13) Z0=Z0+1 4042 END SUBROUTINE EFGHI(B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F 1,B,W,M,I0,V7,H10,H11,H12,H13) INTEGER B0,Z9,Z0,I1,J1,M1,M2,E3,F1,I0,V7 DOUBLE PRECISION H9(E3) INTEGER F(B0),B(F1),W(B0),M(V7) INTEGER H10(E3),H11(E3),H12(E3),H13(E3) 4102 I1=M(I0+1)-3 J1=M(I0+1)-1 4106 M1=M(I0+1)-4 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, H10,H11,H12,H13) 4110 I1=M(I0+1)-3 J1=M(I0+1)-2 4114 M1=M(I0+1)-1 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, H10,H11,H12,H13) 4118 I1=M(I0+1)-5 J1=M(I0+1)-3 4122 M1=M(I0+1)-4 4124 Z9=Z0 CALL CM (H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12,H1 3) Z0=Z0+1 4130 I1=M(I0+1)-3 J1=M(I0+1)-5 4134 M1=M(I0+1)-2 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, H10,H11,H12,H13) 4138 END SUBROUTINE EDEC (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, M,I0,V7,H10,H11,H12,H13) INTEGER B0,Z9,Z0,I1,J1,M1,M2,E3,F1,I0,V7 DOUBLE PRECISION H9(E3) INTEGER F(B0),B(F1),W(B0),M(V7) INTEGER H10(E3),H11(E3),H12(E3),H13(E3) 4210 I1=M(I0+1)-5 J1=M(I0+1)-1 4214 M1=M(I0+1)-3 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, H10,H11,H12,H13) 4218 I1=M(I0+1)-5 J1=M(I0+1)-2 4222 M1=M(I0+1)-1 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, H10,H11,H12,H13) 4226 I1=M(I0+1)-5 J1=M(I0+1)-4 4230 M1=M(I0+1)-2 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, H10,H11,H12,H13) 4234 END SUBROUTINE JK4(D9,M,I0,V7,Z0,H9,E3,M2,B0,F,F1,B, W,J,V0,H10,H11,H12,H13) INTEGER D9,I0,Z0,V7,J0,Z9,I1,M1,J1,E3,M2,B0,F1, V0 DOUBLE PRECISION H9(E3) INTEGER M(V7),F(B0),B(F1),W(B0),J(V0) INTEGER H10(E3),H11(E3),H12(E3),H13(E3) 4302 I1=M(I0)+1 J1=M(I0+1) 4306 M1=M(I0) CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, H10,H11,H12,H13) 4310 IF(J(I0)==12) D9=D9-1 4312 IF(J(I0)==4) D9=D9-1 4318 DO J0=2,D9 4320 IF(J0==D9) GOTO 4342 4322 I1=M(I0)+J0-1 J1=M(I0+1)+J0-1 4326 M1=M(I0+1)+J0-2 4328 IF(J(I0)==10) GOTO 4340 4330 IF(J(I0)==12) GOTO 4340 4332 Z9=Z0 CALL CM (H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12,H1 3) Z0=Z0+1 4338 GOTO 4342 4340 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, H10,H11,H12,H13) 4342 IF(J(I0)==4) GOTO 4350 4344 IF(J(I0)==12) GOTO 4350 4346 IF(J0/=D9) GOTO 4350 4348 GOTO 4370 4350 I1=M(I0+1)+J0-1 J1=M(I0)+J0-1 4354 M1=M(I0)+J0 4356 IF(J(I0)==10) GOTO 4368 4358 IF(J(I0)==12) GOTO 4368 4360 Z9=Z0 CALL CM (H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12,H1 3) Z0=Z0+1 4366 GOTO 4370 4368 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W, H10,H11,H12,H13) 4370 CONTINUE ENDDO 4372 IF(J(I0)==4) GOTO 4382 4374 IF(J(I0)==12) GOTO 4382 4376 I1=M(I0+1)+D9-2 J1=M(I0)+D9-1 4380 GOTO 4386 4382 I1=M(I0)+D9-1 J1=M(I0)+D9 4386 M1=M(I0+1)+D9-1 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,H10,H11,H12,H13) 4390 END 代码说明： （１）本程序大量的使用了子程序结构，以便今后可能的使用，其中从WENDU到DANGANG是计算部分的程序，主要是进行单元刚度矩阵的形成，整体刚度矩阵的组集、用列主元素高斯消去法求解焓以及利用（４）式求出温度。DANGANG以后的子程序都是为了划分网格而写的处理程序。 （２）本程序中的SCORE.DAT，是数据输入文件，本程序没有在程序中直接赋值。 （３）本程序中的OUT1.DAT与OUT2.DAT是网格划分结果输出文件，其中OUT2.DAT文件时为了方便VB的画图程序，用特定格式输出的。OUT4.DAT是计算结果输出文件，直接输出各个节点的温度。 （４）本程序中的时间步长选择为10秒。 结果分析： 程序运行结果较多，为了说明问题，选择两个时刻的结果复制如下： I= 1H= 1160.3767 I= 2H= 1510.4332 I= 3H= 1531.5004 I= 4H= 1354.9724 I= 5H= 1447.7742 I= 6H= 1790.0000 I= 7H= 1790.0000 I= 8H= 1790.0000 I= 9H= 1441.3875 I= 10H= 1488.2260 I= 11H= 1790.0000 I= 12H= 1817.0209 I= 13H= 1790.0000 I= 14H= 1452.3765 I= 15H= 1146.8425 I= 16H= 1513.0586 I= 17H= 1478.9375 I= 18H= 1515.7268 I= 19H= 1147.2844 这一组是经过两次迭代之后，也就是20秒之后的结果，下面一组是第四次迭代也就是40秒后的结果： I= 1H= 734.8505 I= 2H= 1094.4928 I= 3H= 1120.9365 I= 4H= 978.0062 I= 5H= 1050.5867 I= 6H= 1790.0000 I= 7H= 1711.4783 I= 8H= 1477.0049 I= 9H= 945.7737 I= 10H= 923.2153 I= 11H= 1610.6188 I= 12H= 1756.1940 I= 13H= 1609.3173 I= 14H= 1053.1027 I= 15H= 564.1854 I= 16H= 1154.9806 I= 17H= 1124.1097 I= 18H= 1125.2209 I= 19H= 638.7786 比较这两组结果并分析，我们可以发现以下规律： （１）40秒之后各个节点的温度显然比20秒的温度低了很多，这是因为有更长的热交换时间； （２）靠近边界的节点温度显著的比内部节点的温度降低的快，尤其是在转角处的节点温度更是低很多，这一点在１、１５、１９节点上表现的很明显，这也是符合客观规律的，因为在转角处的热量将散发的更快。（３）内部节点如６、７节点的温度降低的很慢，这显然是因为它们处于物体内部，远离边界的交换热量。 综合以上几点，可以基本认为，程序的结果是体现了客观的传热学规律的要求。 1、已知：一块厚度为0.1mm 的无限大平板，具有均匀内热源，q ＝50×103W/m 3，，导热系数K ＝10W/m.℃，一侧边界给定温度为75℃，另一侧对流换热，T f ＝25℃，，h=50W/m 2.℃,求解稳态分布。（边界条件用差分代替微分和能量平衡法），画图。(内,外节点) 2、试以下述一维非稳态导热问题为模型，编写求解一维非稳态扩散型问题的通用程序： 00 00000()()()() L L f x x x x L fL L x x x x T T k s c x x T k h T T W x T k h T T W x T T x τρτ =====???+=????=-+??-=-+?= 其中，x 是空间坐标变量，τ是时间坐标变量，T 是温度（分布），k 是材料的导热系数，s 是内热源强度，ρ是材料的密度，c 是材料的比热，h 0和h L 分别是x 0和x L 处流体与固体壁面间的换热系数，而T f0和T fL 分别是固体壁两侧流体的温度，W 0和W L 是x 0和x L 处(非对流换热)热流密度，T 0(x )是固体壁内初始温度分布。注意k 、ρ、c 、s 、h 0 、h L 、W 0和W L 均可以是温度T 和/或空间坐标x 的函数。 具体要求： 1) 将数学模型无量纲化； 2) 考虑各种可能的边界条件和初始条件组合 3) 提供完整的程序设计说明，包括数学推导过程和程序使用说明 3、对于有源项的一维稳态方程， s dx d T dx d u dx d +=)()(φφρ 已知 x=0,φ=0,x=1, φ=1.源项S=0.5-X 利用迎风格式、混合格式、乘方格式求解φ的分布. 数值计算大作业 一、用数值方法求解尺度为100mm×100mm 的二维矩形物体的稳态导热问题。物体的导热系数λ为1.0w/m·K。边界条件分别为： 1、上壁恒热流q=1000w/m2; 2、下壁温度t1=100℃； 3、右侧壁温度t2=0℃; 4、左侧壁与流体对流换热，流体温度tf=0℃，表面传热系数 h 分别为1w/m2·K、10 w/m2·K、100w/m2·K 和1000 w/m2·K; 要求： 1、写出问题的数学描述； 2、写出内部节点和边界节点的差分方程； 3、给出求解方法； 4、编写计算程序(自选程序语言)； 5、画出4个工况下的温度分布图及左、右、下三个边界的热流密度分布图； 6、就一个工况下（自选）对不同网格数下的计算结果进行讨论； 7、就一个工况下（自选）分别采用高斯迭代、高斯——赛德尔迭代及松弛法（亚松弛和超松弛）求解的收敛性（cpu 时间，迭代次数）进行讨论； 8、对4个不同表面传热系数的计算结果进行分析和讨论。 9、自选一种商业软件（fluent 、ansys 等）对问题进行分析，并与自己编程计算结果进行比较验证（一个工况）。（自选项） 1、写出问题的数学描述 设H=0.1m 微分方程 22220t t x y ??+=?? x=0，0 y=H ，0 传热学 二维稳态导热问题的数值解法 杨达文2011151419 赵树明2011151427 杨文晓2011151421 吴鸿毅2011151416 第一题： a=linspace(0,0.6,121); t1=[60+20*sin(pi*a/0.6)]; t2=repmat(60,[80 121]); s=[t1;t2]; %构造矩阵 for k=1:10000000 %理论最大迭代次数，想多大就设置多大S=s; for j=2:120 for i=2:80 S(i,j)=0.25*(S(i-1,j)+S(i+1,j)+S(i,j-1)+S(i,j+1)); end end if norm(S-s)<0.0001 break; %如果符合精度要求，提前结束迭代else s=S; end end S %输出数值解 数值解数据量太大，这里就不打印出来，只画出温度分布。 画出温度分布： figure(1) xx=linspace(0,0.6,121); yy=linspace(0.4,0,81); [x,y]=meshgrid(xx,yy); surf(x,y,S) axis([0 0.6 0 0.4 60 80]) grid on xlabel('L1') ylabel('L2') zlabel('t（温度）') .60.66666777778L 1 L 2t （温度） A0=[S(:,61)]; for k=1:81 B1(k)=A0(81-k+1); end B1 %x=L1/2时y方向的温度 A1=[S(41,:)] %y=L2/2时x方向的温度 x=0:0.005:0.6; y=0:0.005:0.4; A2=60+20*sin(pi*x/0.6)*((exp(pi*0.2/0.6)-exp(-pi*0.2/0.6))/2)/((exp(pi*0.4/0.6)-exp(-pi*0.4/0.6) )/2) %计算y=L2/2时x方向的解析温度 B2=60+20*sin(pi*0.3/0.6)*((exp(pi*y/0.6)-exp(-pi*y/0.6))/2)/((exp(pi*0.4/0.6)-exp(-pi*0.4/0.6))/ 2) %计算x=L1/2时y方向的解析温度 figure(2) subplot(2,2,1); plot(x,A1,'g-.',x,A2,'k:x'); %画出x=L1/2时y方向的温度场、画出x=L1/2时y方向的解析温度场曲线 xlabel('L1');ylabel('t温度'); title('y=L2/2'); legend('数值解','解析解'); subplot(2,2,2); plot(x,A1-A2); %画出具体温度场与解析温度场的差值曲线 xlabel('L1');ylabel('差值'); title('y=L2/2时，比较=数值解-解析解'); subplot(2,2,3); plot(y,B1,'g-.',y,B2,'k:x'); %画出y=L2/2时x方向的温度场、画出y=L2/2时x方向的解析温度场曲线 xlabel('L2');ylabel('t温度'); title('x=L1/2'); legend('数值解','解析解'); subplot(2,2,4); plot(y,B1-B2); %画出具体温度场与解析温度场的差值曲线 xlabel('L2');ylabel('差值'); title('x=L1/2时，比较=数值解-解析解'); y=L2/2时x方向的温度： 60 60.1635347276130 60.3269574318083 60.4901561107239 60.6530189159961 60.8154342294146 60.9772907394204 61.1384775173935 61.2988840936779 61.4584005332920 61.6169175112734 61.7743263876045 61.9305192816696 62.0853891461909 62.2388298405943 62.3907362037523 62.5410041260577 62.6895306207746 62.8362138946214 62.9809534175351 63.1236499915702 63.2642058188844 63.4025245687647 63.5385114436490 63.6720732440951 63.8031184326565 63.9315571966177 64.0573015095482 64.1802651916318 64.3003639687311 64.4175155301449 64.5316395850212 64.6426579173846 64.7504944397430 64.8550752452343 64.9563286582797 65.0541852837075 第四章循环流化床锅炉炉内传热计算 循环流化床锅炉炉膛中的传热是一个复杂的过程，传热系数的计算精度直接影响了受热面设计时的布置数量，从而影响锅炉的实际出力、蒸汽参数和燃烧温度。正确计算燃烧室受热面传热系数是循环流化床锅炉设计的关键之一，也是区别于煤粉炉的重要方面。 随着循环流化床燃烧技术的日益成熟，有关循环流化床锅炉的炉膛传热计算思想和方法的研究也在迅速发展。许多著名的循环流化床制造公司和研究部门在此方面也做了大量的工作，有的已经形成商业化产品使用的设计导则。 但由于技术保密的原因，目前国内外还没有公开的可以用于工程使用的循环流化床锅炉炉膛传热计算方法，因此对它的研究具有重要的学术价值和实践意义。 清华大学对CFB锅炉炉膛传热作了深入的研究，长江动力公司、华中理工大学、浙江大学等单位也对CFB锅炉炉膛中的传热过程进行了有益的探索。根据已公开发表的文献报导，考虑工程上的方便和可行，本章根椐清华大学提出的方法，进一步分析整理，作为我们研究的基础。为了了解CFB锅炉传热计算发展过程，也参看了巴苏的传热理论和计算方法，浙江大学和华中理工大学的传热计算与巴苏的相近似。 4.1 清华的传热理论及计算方法 4.1.1 循环流化床传热分析 CFB锅炉与煤粉锅炉的显著不同是CFB锅炉中的物料(包括煤灰、脱硫添加剂等)浓度C p 大大高于煤粉炉，而且炉内各处的浓度也不一样，它对炉内传热起着重要作用。为此首先需要计算出炉膛出口处的物料浓度C p，此处浓度可由外循环倍率求出。而炉膛不同高度的物料浓度则由内循环流率决定，它沿炉膛高度是逐渐变化的，底部高、上部低。近壁区贴壁下降流的温度比中心区温度低的趋势，使边壁下降流减少了辐射换热系数；水平截面方向上的横向搅混形成良好的近壁区物料与中心区物料的质交换，同时近壁区与中心区的对流和辐射的热交换使截面方向的温度趋于一致，综合作用的结果近壁区物料向壁面的辐射加强，总辐射换热系数明显提高。在计算水冷壁、双面水冷壁、屏式过热器和屏式再热器时需采用不同的计算式。物料浓度C p对辐射传热和对流传热都有显著影响。燃烧室的平均温度是床对受热面换热系数的另一个重要影响因素。床温的升高增加了烟气辐射换热并提高烟气的导热系数。虽然粒径的减小会提高颗粒对受热面的对流换热系数，在循环流化床锅炉条件下，燃烧室内部的物料颗粒粒径变化较小，在较小范围内的粒径变化时换热系数的变化不大，在进行满负荷传热计算时可以忽略，但在低负荷传热计算时，应该考虑小的颗粒有提高传热系数的能力。 炉内受热面的结构尺寸，如鳍片的净宽度、厚度等，对平均换热系数的影响也是非常明显的。鳍片宽度对物料颗粒的团聚产生影响；另一方面，宽度与扩展受热面的利用系数有关。根 传热学大作业报告二维稳态计算 院系：能源与环境学院 专业：核工程与核技术 姓名：杨予琪 学号：03311507 一、原始题目及要求 计算要求： 1. 写出各未知温度节点的代数方程 2. 分别给出G-S 迭代和Jacobi 迭代程序 3. 程序中给出两种自动判定收敛的方法 4. 考察三种不同初值时的收敛快慢 5. 上下边界的热流量（λ=1W/(m ℃）） 6. 绘出最终结果的等值线 报告要求： 1. 原始题目及要求 2. 各节点的离散化的代数方程 3. 源程序 4. 不同初值时的收敛快慢 5. 上下边界的热流量（λ=1W/(m ℃）） 6. 计算结果的等温线图 7. 计算小结 二、各节点的离散化的代数方程 左上角节点 )(21 1,22,11,1t t t += 右上角节点 )(2 15,24,15,1t t t += 左下角节点 C t ?=1001,5 右下角节点 )2(211,24,55,5λ λ x h t t x h t ?++?+= 左边界节点 C t i ?=1001,，42≤≤i 上边界节点 C t j ?=200,1，42≤≤j 右边界节点 )2(415,15,14,5,+-++= i i i i t t t t ，42≤≤i 下边界节点 )42()2(211,51,5,4,5∞+-?+++?+=t x h t t t x h t j j j j λλ ，42≤≤j 内部节点 )(2 1,1,11,1,,j i j i j i j i j i t t t t t +-+-+++= ，4,2≤≤j i 三、源程序 1、G-S 迭代法 t=zeros(5,5); t0=zeros(5,5); dteps=0.0001; for i=2:5 %左边界节点 t(i,1)=100; end for j=2:4 %上边界节点 t(1,j)=200; end t(1,1)=(t(1,2)+t(2,1))/2; t for k=1:100 for i=2:4 %内部节点 for j=2:4 t(i,j)=(t(i-1,j)+t(i+1,j)+t(i,j-1)+t(i,j+1))/4; end end t(1,5)=(t(1,4)+t(2,5))/2;%右上角节点 for i=2:4;%右边界节点 t(i,5)=(2*t(i,4)+t(i-1,5)+t(i+1,5))/4; end for j=2:4; %下边界节点 硕士研究生《高等工程热力学与传热学》作业 查阅相关资料,回答以下问题： 1、一滴水滴到120度和400度的板上，哪个先干？试从传热学的角度分析？ 答：在大气压下发生沸腾换热时,上述两滴水的过热度分别是△ t=tw–ts=20℃和△t=300℃,由大容器饱和沸腾曲线,前者表面发生的是泡态沸腾,后者发生膜态沸腾。虽然前者传热温差小,但其表面传热系数大,从而表面热流反而大于后者。所以水滴滴在120℃的铁板上先被烧干。 2、锅铲、汤勺、漏勺、铝锅等炊具的柄用木料制成，为什么？ 答：是因为木料是热的不良导体,以便在烹任过程中不烫手。 3、滚烫的砂锅放在湿地上易破裂。为什么？ 答：这是因为砂锅是热的不良导体, 如果把烧得滚热的砂锅,突然放到潮湿或冷的地方,砂锅外壁的热就很快地被传掉，而壁的热又一下子传不出来，外壁冷却很快的收缩，壁却还很热，没什么收缩，加以瓷特别脆，所以往往裂开。 或者：烫砂锅放在湿地上时,砂锅外壁迅速放热收缩而壁温度降低慢,砂锅外收缩不均匀,故易破裂。 4、往保温瓶灌开水时，不灌满能更好地保温。为什么？ 答：因为未灌满时,瓶口有一层空气,是热的不良导体,能更好地防止热量散失。 5、煮熟后滚烫的鸡蛋放入冷水中浸一会儿，容易剥壳。为什么？ 答：因为滚烫的鸡蛋壳与蛋白遇冷会收缩,但它们收缩的程度不一样,从而使两者脱离。 6、用焊锡的铁壶烧水，壶烧不坏，若不装水，把它放在火上一会儿就烧坏了。为什么？ 答：这是因为水的沸点在1标准大气压下是100℃,锡的熔点是232℃,装水烧时,只要水不干,壶的温度不会明显超过100℃,达不到锡的熔点,更达不到铁的熔点,故壶烧不坏.若不装水在火上烧,不一会儿壶的温度就会达到锡的熔点,焊锡熔化,壶就烧坏了。 7、冬壶里的水烧开后，在离壶嘴一定距离才能看见“白气”，而紧靠壶嘴的地方看不见“白气”。这是因为紧靠壶嘴的地方温度高，壶嘴出来的水蒸气不能液化，而距壶嘴一定距离的地方温度低；壶嘴出来的水蒸气放热液化成小水滴，即“白气”。 答：这是因为紧靠壶嘴的地方温度高,壶嘴出来的水蒸气不能液化,而距壶嘴一定距离的地方温度低；壶嘴出来的水蒸气放热液化成小水滴,即“白气”。 8、某些表演者赤脚踩过炽热的木炭，从传热学角度解释为何不会烫伤？不会烫伤的基本条件是什么？ 答：因为热量的传递和温度的升高需要一个过程，而表演者赤脚接触炽热木炭的时间极短，因此在这个极短的时间传递的温度有限，不足以达到令人烫伤的温度，所以不会烫伤。 基本条件：表演者接触炽热木炭的时间必须极短，以至于在这段时间所传递的热量不至于达到灼伤人的温度 取步长δx=0.02。已知x=0,Φ=0;x=1,Φ=1.令k=ρｕ/Γ计算结果图表： 程序及数据结果： 追赶法： #include a[i]=2+0.02*k; b[i]=4; c[i]=2-0.02*k; f[i]=0; } tdma(a,b,c,f,x); for(i=0;i 一维稳态导热的数值计算 1.1物理问题 一个等截面直肋，处于温度t ∞=80 的流体中。肋表面与流体之间的对流换热系数为 h =45W/(m 2?℃),肋基处温度t w =300℃，肋端绝热。肋片由铝合金制成，其导热系数为λ=110W/(m ?℃)，肋片厚度为δ=0.01m ，高度为H=0.1m 。试计算肋内的温度分布及肋的总换热量。 1.2数学描述及其解析解 引入无量纲过余温度θ = t?t ∞t w ?t ∞ ，则无量纲温度描述的肋片导热微分方程及其边界条件： 22 20d m dx θθ-= x=0,θ=θw =1 x=H, 0x θ?=? 其中m = 上述数学模型的解析解为：[()] ()() w ch m x H t t t t ch mH ∞∞--=-? ()()w hp t t th mH m ∞?= - 1.3数值离散 1.3.1区域离散 计算区域总节点数取N 。 1.3.2微分方程的离散 对任一借点i 有：22 2 0i d m dx θ θ??-= ??? 用θ在节点i 的二阶差分代替θ在节点i 的二阶导数，得：211 2 20i i i i m x θθθθ+--+-= 整理成迭代形式：()1122 1 2i i i m x θθθ+-=++ (i=2,3……,N-1) 1.3.3边界条件离散 补充方程为：11w θθ== 右边界为第二类边界条件，边界节点N 的向后差分得：1 0N N x θθ--= ，将此式整理为 迭代形式，得：N 1N θθ-= 1.3.4最终离散格式 11w θθ== ()1122 1 2i i i m x θθθ+-= ++ (i=2,3……,N-1) N 1N θθ-= 1.3.5代数方程组的求解及其程序 假定一个温度场的初始发布，给出各节点的温度初值：01θ，02θ，….，0 N θ。将这些初值代 入离散格式方程组进行迭代计算，直至收敛。假设第K 步迭代完成，则K+1次迭代计算式为： K 11w θθ+= () 11 11 2212i i K K K i m x θθθ+-++= ++ (i=2,3……,N-1) 1 11N K K N θθ-++= #include 二维导热物体温度场的数值模拟 一、物理问题 有一个用砖砌成的长方形截面的冷空气通道，其截面尺寸如下图1-1所示，假设在垂直于纸面方向上用冷空气及砖墙的温度变化很小，可以近似地予以忽略。在下列两种情况下试计算： 砖墙横截面上的温度分布；垂直于纸面方向的每米长度上通过砖墙的导热量。 第一种情况：内外壁分别均匀维持在0℃及30℃； 第二种情况：内外壁均为第三类边界条件，且已知： K m K m W h C t K m W h C t ?=?=?=?=?=∞∞/35.0/93.3,10/35.10,302 22211λ砖墙导热系数 二、数学描写 由对称的界面必是绝热面，可取左上方的四分之一墙角为研究对象，该问题为二维、稳态、无内热源的导热问题。 控制方程： 02 222=??+??y t x t 边界条件： 第一种情况： 由对称性知边界1绝热： 0=w q ； 边界2为等温边界，满足第一类边界条件： C t w ?=0； 边界3为等温边界，满足第一类边界条件： C t w ?=30。 第一种情况： 由对称性知边界1绝热： 0=w q ； 边界2为对流边界，满足第三类边界条件： )()( 2f w w w t t h n t q -=??-=λ； 边界3为对流边界，满足第三类边界条件： )()(2f w w w t t h n t q -=??-=λ。 1 -1图2 -1图 三、方程离散 用一系列与坐标轴平行的间隔0.1m 的二维网格线将温度区域划分为若干子区域，如图1-3所示。 采用热平衡法，利用傅里叶导热定律和能量守恒定律，按照以导入元体（m,n ）方向的热流量为正，列写每个节点代表的元体的代数方程， 第一种情况： 边界点： 边界1（绝热边界）： 5~2)2(4 1 1,11,12,1,m =++= +-m t t t t m m m ， 11~8)2(4 1 1,161,16,15,16=++=+-n t t t t n n n n ， 边界2（等温内边界）： 7,16~7;7~1,6,0,=====n m n m t n m 边界3（等温外边界）： 12,16~2;12~1,1,30,=====n m n m t n m 内节点： 11 ~8,15~6;11~2,5~2)(41 1,1,,1,1,====+++= -+-+n m n m t t t t t n m n m n m n m n m 第二种情况 边界点： 边界1（绝热边界）： 5~2)2(4 1 1,11,12,1 ,m =++=+-m t t t t m m m ， 11~8)2(4 1 1,161,16,15,16=++=+-n t t t t n n n n ， 边界2（内对流边界）： 6~1) 2(2221 11,61,6,5,6=++++= ??-+n Bi t Bi t t t t n n n n ， 3 -1图 中国石油大学(华东) 储建学院热能与动力工程系 《计算传热学程序设计》 设计报告 1引言 有关墙体传热量计算的方法是随着人们对房间负荷计算精度要求的不断提高而不断发展的.考虑辐射强度和周围空气温度综合作用，当外界温度发生周期性的变化时，屋顶内部的温度和热流密度也会发生周期性的变化。 计算题目 有一个用砖墙砌成的长方形截面的冷空气通道，其截面尺寸如图1所示。假设在垂直于纸面方向上冷空气及砖墙的温度变化相对较小，可近似地予以忽略。试计算稳态时砖墙截面的温度分布及垂直于纸面方向1米长度的冷量损失。设砖墙的导热系数为(m·℃)。内、外壁面均为第三类边界条件，外壁面：t f1=30℃,h1=10W(m2·℃);内壁面：t f2=10℃, h2=4W(m2·℃)。 图1 砖墙截面 已知参数 砖墙的基本尺寸，砖墙的导热系数，外壁面的表面传热系数，对应的流体温度，内壁面的表面传热系数，对应的流体温度。 2 物理与数学模型 物理模型 由题知垂直于纸面方向上冷空气及砖墙的温度变化相对较小，可近似予以忽略，墙面为常物性，可以假设： 1）砖墙在垂直于纸面方向上没有导热。 2）由于系统是几何形状与边界条件是对称的，它的中心对称面就是一个绝热边界，这时只需求解1/4个对称区域就可以得到整个区域的解。 数学模型 考虑到对称性，取右下的1/4为研究对象，建立如图2的坐标系。 a 图2 砖墙的稳态导热计算区域 由上述的物理模型与上面的坐标系，该问题的数学模型可直接由导热微分方程简化而来，即 22220T T x y ??+=?? (1) 相应的边界条件是： 1.1 0y T y =?=? 1.5 0x T x =?=? (2) 110 ()f x x T h T T x λ ==?-=-? (3) 111.1 1.1 ()f y y T h T T y λ ==?-=-? (4) 22(0.5,00.6)(0.5,00.6) ()f x y x y T h T T x λ =<<=< 《传热学》上机大作业 二维导热物体温度场的数值模拟 学校：西安交通大学 姓名：张晓璐 学号:10031133 班级：能动A06 一．问题（4-23） 有一个用砖砌成的长方形截面的冷空气通道，形状和截面尺寸如下图所示，假设在垂直纸面方向冷空气和砖墙的温度变化很小，差别可以近似的予以忽略。在下列两种情况下计算：砖墙横截面上的温度分布；垂直于纸面方向上的每米长度上通过墙砖上的导热量。 第一种情况：内外壁分别维持在10C ?和30C ? 第二种情况：内外壁与流体发生对流传热，且有C t f ?=101， )/(2021k m W h ?=，C t f ?=302，)/(422k m W h ?=，K m W ?=/53.0λ 二．问题分析 1.控制方程 02222=??+??y t x t 2.边界条件 所研究物体关于横轴和纵轴对称，所以只研究四分之一即可，如下图： 对上图所示各边界： 边界1：由对称性可知：此边界绝热，0=w q 。 边界2：情况一：第一类边界条件 C t w ?=10 情况二：第三类边界条件 )()( 11f w w w t t h n t q -=??-=λ 边界3：情况一：第一类边界条件 C t w ?=30 情况二：第三类边界条件 )()( 22f w w w t t h n t q -=??-=λ 三：区域离散化及公式推导 如下图所示，用一系列和坐标抽平行的相互间隔cm 10的网格线将所示区域离散化，每个交点可以看做节点，该节点的温度近似看做节点所在区域的平均温度。利用热平衡法列出各个节点温度的代数方程。 第一种情况： 内部角点： 、室内一根水平放置的无限长的蒸汽管道，其保温层外径d=583 mm，外表面 实测平均温度及空气温度分别为，此时空气与管道外 表面间的自然对流换热的表面传热系数h=3.42 W /(m2 K),墙壁的温度近似取为 室内空气的温度，保温层外表面的发射率 问：(1)此管道外壁的换热必须考虑哪些热量传递方式； (2)计算每米长度管道外壁的总散热量。(12分) 解： (1)此管道外壁的换热有辐射换热和自然对流换热两种方式。 (2)把管道每米长度上的散热量记为qi 当仅考虑自然对流时，单位长度上的自然对流散热 q i,c =二d h t =二dh (j - t f ) = 3.14 0.583 3.42 (48 - 23 ) 二156 .5(W / m) 近似地取墙壁的温度为室内空气温度，于是每米长度管道外表面与室内物体及墙壁 之间的辐射为： q i厂d (T； -T；) = 3.14 0.583 5.67 10》0.9 [(48 273)4-(23 273)4] = 274.7(W /m) 总的散热量为q i = q i,c +q i,r = 156.5 +274.7 = 431.2(W/m) 2、如图所示的墙壁，其导热系数为50W/(m- K),厚度为50mm在稳态情况下的 墙壁内的一维温度分布为：t=200-2000x 2，式中t的单位为°C, x单位为m 试 求： t (1) 墙壁两侧表面的热流密度； (2) 墙壁内单位体积的内热源生成的热量 2 t =200 —2000x 解：(1)由傅立叶定律: ① dt W q ' (―4000x) = 4000二x A dx 所以墙壁两侧的热流密度： q x _. =4000 50 0.05 =10000 (1)由导热微分方程 茫?生=0得: dx 扎 3、一根直径为1mm 勺铜导线，每米的电阻为2.22 10 。导线外包有厚度为 0.5mm 导热系数为0.15W/(m ? K)的绝缘层。限定绝缘层的最高温度为 65°C,绝 缘层的外表面温度受环境影响，假设为40°C 。试确定该导线的最大允许电流为多 少？ 解：(1)以长度为L 的导线为例，导线通电后生成的热量为I 2RL ,其中的一部分 热量用于导线的升温，其热量为心务中：一部分热量通过绝热层的 导热传到大气中，其热量为：门二 1 , d In 2 L d 1 根据能量守恒定律知：l 2RL -门 述二厶E = I 2RL -门 即 E = — L dT m = I 2RL - t w1 _tw2 4 di 1 , d 2 In 2 L d 1 q v 、d 2t ——' 2 dx =-(7000)= 4000 50 二 200000 W/m 3 t w1 - t w2 。 2 q x 卫=4000.： 0 = 0 传热学经典计算题 热传导 1. 用热电偶测量气罐中气体的温度。热电偶的初始温度为20℃，与气体的表面传热系数为()210/W m K ?。热电偶近似为球形，直径为0.2mm 。试计算插入10s 后，热电偶的过余温度为初始过余温度的百分之几？要使温度计过余温度不大于初始过余温度的1％，至少需要多长时间？已知热电偶焊锡丝的()67/W m K λ=?，7310ρ= 3/kg m ，()228/c J kg K =?。 解： 先判断本题能否利用集总参数法。 3 5100.110 1.491067hR Bi λ--??===?＜0.1 可用集总参数法。 时间常数 3 73102280.110 5.563103c cV c R hA h ρρτ-??===?= s 则10 s 的相对过余温度 0θθ=exp c ττ??-= ???exp 1016.65.56??-= ???％ 热电偶过余温度不大于初始过余温度1％所需的时间，由题意 0θθ=exp c ττ??- ??? ≤0.01 exp 5.56τ?? - ???≤0.01 解得 τ≥25.6 s 1、空气以10m/s 速度外掠0.8m 长的平板，C t f 080=,C t w 030=，计算 该平板在临界雷诺数c e R 下的c h 、全板平均表面传热系数以及换热量。 （层流时平板表面局部努塞尔数 3/12/1332.0r e x P R Nu =，紊流时平板表面局部努塞尔数3/15/40296.0r e x P R Nu =，板宽为1m ，已知5105?=c e R ，定性 温度C t m 055=时的物性参数为： )/(1087.22K m W ??=-λ，s m /1046.1826-?=ν，697.0=r P ） 解：(1)根据临界雷诺数求解由层流转变到紊流时的临界长度 C t t t w f m 055)(21=+=,此时空气得物性参数为： )/(1087.22K m W ??=-λ，s m /1046.1826-?=ν，697.0=r P )(92.0101046.1810565m u R X ul R c c e c e =???==?=-ν ν 由于板长是0.8m ，所以，整个平板表面的边界层的流态皆为层流 ? ==3/12/1332.0r e x P R hl Nu λ)/(41.7697.0)105(8.01087.2332.0332.023/12/1523/12 /1C m W P R l h r e c c ?=????==-λ (2)板长为0.8m 时，整个平板表面的边界层的雷诺数为： 561033.41046.188.010?=??==-νul R e 全板平均表面传热系数: )/(9.13697.0)1033.4(8.01087.2664.0664.023/12/1523/12 /1C m W P R l h r e c ?=????==-λ 全板平均表面换热量W t t hA w f 9.557)3080(18.09.13)(=-???=-=Φ 储运与建筑工程学院能源与动力工程系 计算传热学课程大作业报告 作业题目：代数方程组的求解 学生姓名：田 学号： 专业班级：能动1 2017年9月23日 目录 一、计算题目 (3) 二、离散方程 (3) 三、程序设计 (4) 3.1 高斯赛德尔迭代法 (4) 3.2 TDMA法 (5) 四、程序及计算结果验证 (6) 五、网格独立性考核.................... 错误！未定义书签。 3.1 高斯赛德尔迭代法 (7) 3.2 TDMA法 (8) 六、结果分析与结论 (8) 3.1 高斯赛德尔迭代法 (9) 3.2 TDMA法 (10) 一、计算题目 分别用高斯赛德尔迭代和TDMA 方法求解方程 2 2dx d dx d u φφρΓ= （1） 在Γ u ρ=-5，-1,0,1,5情况下的解，并表示在图中。 其中，x =0，φ=0；x =1，φ=1. 二、离散方程 采用控制容积法： 即??Γ=e 22w e w dx d dx d u φφ ρ（2） ) )()(()2 2 ( w W P e P E p w p e x x u δφφδφφφφφφρ---Γ=+- +（3） 假设均分网格，则有x x x w e ?==)()(δδ 上式则变为： )2(2)(W P E W E u x φφφφφρ+-Γ=-?（4） 即11)2()2(4-+?+Γ+?-Γ=Γi i i u x u x φρφρφ（5） 11)421()421(-+Γ ?-+Γ?-=i i i u x u x φρφρφ（6） 三、程序设计 3.1 高斯赛德尔迭代法 由已知公式 11)421()421(-+Γ ?-+Γ?-=i i i u x u x φρφρφ可设计高斯赛德尔迭代C 语言程序如下： #include 课程编号：13SD02010340 课程名称：传热学 上课时间：2014年春季 电子元器件散热方法研究 姓名： 学号： 班级： 所在学院： 任课教师： 摘要：随着电子器件的高频、高速以及集成电路技术的迅速发展和技术的进步，电子元器件的总功率密度大幅度增长而物理尺寸却越来越小，热流密度也随之增加，所以高温的 温度环境势必会影响电子元器件的性能，这就要求对其进行更加高效的热控制。因此，有 效解决电子元器件的散热问题已成为当前电子元器件和电子设备制造的关键技术。本文针 对电子元器件的散热与冷却问题，综述了当前应用研究中不同的散热和冷却方法，并进行 了适当的分析。 关键词热管理; 冷却; 电子器件 近些年来,电子技术的快速发展。电子器件的高频、高速以及集成电路的密集和小型化,使得单位容积电子器件的总功率密度和发热量大幅度地增长,从而使电子器件的冷却问题 变得越来越突出。如: 大型计算机的芯片热流量已达到了60 W/ cm2,到2000 年已经超过了,目前最高已达到200 W/ cm2。特别是由于MEMS技术突飞猛进,使得电子元器件的尺寸越来越小,已经从微米量级进入到了亚微米量级。尽管随着器件或系统尺寸的减小, 消耗功率也会有所减小, 但为了完成一定的任务,可减小的余地非常有限,这使得为系统内的热流密度非 常大, 据报道可达, 远远高出航天飞行器回归地球与大气摩擦时产生的惊人的高热流密度。在微系统中可能出现的高热流密度对于电子器件是致命的, 然而使用传统的冷却技术要使 如此高的热流密度在短时间内散去几乎是不现实的; 另一方面, 电子器件工作的可靠性对 温度十分敏感, 器件温度在70～80 水平上每增加1, 可靠性就会下降5%。因而电子产品的 开发、研制中必须要充分考虑到良好的散热手段, 才能保证产品的可靠性和表观。由于电 子元器件的小型化、微型化和集成化,所采用的散热和冷却手段必须要求具有紧凑性、可靠性、灵活性、高散热效率等特点。 1 电子元器件的散热或冷却方法 电子元器件的高效散热问题与传热学、流体力学等原理的应用密切相关。电子器件散 热的目的是对电子设备的运行温度进行控制,以保证其工作的稳定性和可靠性。这其中涉及了与传热有关的散热或冷却方式、材料等多方面内容。从应用的角度看,常用的方法主要有: 自然散热或冷却、强制散热或冷却、液体冷却、制冷方式、疏导方式、热隔离方式和PCM 温度控制方法等。 1.1 自然散热或冷却方法 自然散热或冷却方法是指不使用任何外部辅助能量的情况下,实现局部发热器件向周 围环境散热达到温度控制的目的,这其中通常都包含了导热、对流和辐射三种主要传热方式, 其中对流以自然对流方式为主。自然散热或冷却往往适用对温度控制要求不高、器件发热 的热流密度不大的低功耗器件和部件,以及密封或密集组装的器件不宜采用其它冷却技术 的情况下。有时,在对散热能力要求不高时也常常利用电子器件自身特点增强与邻近热沉的导热或辐射、通过结构设计强化自然对流,在一定程度上提高系统向环境散热能力。计算传热学
传热学数值计算大作业2014011673
哈工程传热学数值计算大作业
传热系数计算方法
传热学大作业报告 二维稳态导热
生活中的传热学(问答题整理答案)
计算传热学中国石油大学(华东)第四章大作业
传热学上机C程序源答案之一维稳态导热的数值计算
西安交通大学传热学大作业---二维温度场热电比拟实验
计算传热学程序设计
西安交通大学传热学大作业
传热学计算例题
传热学经典计算题
计算传热作业1
传热学大作业