63 DO 10 к е 1,NY

64 AK = К - 1

65 Y(K) = AK*DY

66 SIYCK) SIN(PIH*Y{K))

67 10 T(K,J) 80.*(Y(K) - SIX(J)*SIY(K)) + 20.

68 11 CONTINUE

69 IFdBC .EQ. DGOTO 14

70 DO 12 J s 1,NX

71 12 T{NYS,J) = T(NYP,J) + 160.*DY

72 DO 13 К 1,NYS

73 13 T(K,NXS) T(K,NXP)

74 14 CONTINUE

75 GAMH > GAM

76 GAM 0.

77 BETH = BET

78 BET 0.5

79 TIM = 0.

80 TMCH = -0.00001 + DTMCH

81 DIM = - 0.5*ALE*PI*PI

82 DO 17 J 1,NX

83 DO 15 К = 1,NY

84 DT(K,J) > 0.

85 15 R(K,J) 0.

86 DO 16 К = 1,5

87 16 B(K,J) * 0.

88 17 CONTINUE

89 18 CCXA = BET*DTIM*ALF/DX/DX/(1.+GAM)

90 CCYA BET*DTIM*ALF/DY/DY/(1.+GAM)

91 TFAC = EXP(DIM*(TIM+DTIM)) - EXP(DIM*TIM)

92 С

93 CALL REDIF(ALF,GAM,DTIM,NXN,NYN)

94 С

95 С TRIDIAGONAL SYSTEMS IN THE X-DIRECTIOM

96 С

97 DO 21 К = 2,NYM

98 DO 19 J 2,NXM

99 JM = J - 1

100 B(2,JM) EMXd) - CCXA

101 B<3,JM) s EMX(2) + 2.*CCXA

102 B(4,JM) EMX(3) - CCXA

103 19 RRT(JM) * R{K,J)

104 B(2,l) = 0.

105 IFdBC .EQ. 1)DT(K,NX) = - 80.*SIY(K)*TFAC

106 IFdBC .EQ. 1)DUM (EMY(1)-CCYA)*SIY(K-1) + (EMY(2)+2*CCYA)*SIY(K>

107 1 + (EMY(3)-CCYA)*SIY(K+1)

108 IFdBC .EQ. l)RRT(JM) = RRT(JM) + В(4, JM) *80.*TFAC*DUM

109 IFdBC .EQ. 2)DT(K,NXS) = DT(K,NXP)

110 IFdBC .EQ. 2)B(2,JM) = B(2,JM) + B(4,JM)

111 B(4,JM) = 0.

112 С

113 CALL BANFAC(B,NXPP,1)

114 CALL BANS0L(RRT,DDT,B,NXPP,1)

115 С

116 DO 20 J = 2,NXN

117 JM = J - 1

118 20 R(KcJ) = DDT(JM)

119 21 CONTINUE

120 С

121 С TRIDIAGONAL SYSTEMS IN THE Y-DIRECTION

122 С

123 DO 24 J 2,NXN

124 DO 22 К = 2,NYN

125 KM = К - 1

126 В(2,КН) EHY(l) - CCYA

127 ВОДИ) = E;MY(2) + 2.*CCYA

128 B(4,KH) EMY(3) - CCYA

129 22 RRT(KH) = R(K,J)

130 В(2Д) = 0.

131 IFdBC .EQ. 1)DT(NY,J) = - 80.*SIX(J) *TFAC

1>2 IFdBC .EQ. l)RRT(KM) = RRT(KM) - 3(4 ДМ) *DT(MY, J)

133 IFdBC .EQ. 2)DT(NYS,J) DT(NYP,J)

134 IFdBC .EQ. 2)В(2,КИ) B(2,1CH) -i- B(4,KH)

135 B(4,KH) 0.

136 с

137 CALL BANFAC(B,NYPPД)

138 CALL BANSOL(RRT,DDT,B,NYPPД)

139 с

140 с INCREMENT T

141 с

142 DO 23 К = 2,NYN

143 KM = К - 1

144 DT(K,J) = DDT(KM)

145 23 T(K,J) = T(K,J) + DDT(KH)

146 24 CONTINUE

147 IFdBC .EQ. 1)DT(NY,NX) - 80.*TFAC

148 IFdBC .EQ. 2)DT(NYS,NXS) = DT(NYP,NXP)

149 TIM TIM + DTIM

150 TFAC = EXP(DIM*TIM)

151 IFdBC .EQ. 2)GOTO 26

152 DO 25 К 1,NY

153 T(K,NX) = 20. + 80.*(Y(K) - TFAC*SIY(K))

154 25 T(NY,K) = 20. + 80.*(1. - TFAC*SIY(K)) 155. GOTO 29

156 26 DO 27 J * 1,NX

157 27 T(NYS,J) = T(NYP,J) + 160.*DY

158 DO 28 К = 1,NYS JL59 28 T(K,NXS) = T(K,NXP)

160 29 GAM GAMH

161 BET BETH

162 IF(TIM .GE. TMCH)GOTO 30

163 IF(TIM .LE. TMAX)GOTO 18

164 GOTO 37

165 С

166 С COMPARE VITH EXACT SOLUTION

167 С

168 30 TMCH = ТИСН + DTHCH

169 SUNT = 0.

170 DO 33 KA = 1,NY

171 К 1 + NY - KA

172 DO 31 J = 1,NX

173 ТЕ 20. + 80.MY(K) - TFAC*SIX(J) *SIY(K))

174 ERR(K,J) = T(K,J) - ТЕ

175 31 SUNT = SUNT + (T(K,J) - TE)**2

176 WRITE(6,32)(T(K,J),J=1,NX)

177 32 FORMATC T=M1F7.3)

178 33 CONTINUE

179 DO 34 KA = 1,NY

180 К = 1 + NY - KA

181 34 WRITE(6,35)(ERR(K,J),J=1,NX)

182 35 rORMATC E M1F7.3)

183 RMST DSQRT(SUNT/AN)

184 WRITE(6,36)TIM,RrtST

185 36 FORMATC TIME =4E11.4, RMST =4D11.4)

186 IF(TIM .LE. TMAX)GOTO 18

187 37 CONTINUE

188 STOP

189 END

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

SUBROUTINE REDIF(ALF, GAM,DTIM,NXN, NYN)

EVALUATES RIGHT-HAND SIDE OF THE UNSTEADY HEAT CONDUCTION EQUATION

DIMENSION DT(41,41),T(41,41),R{41,41), 1DMX(3,3),DMY(3,3),EMX(3),EMY(3) COMMON DX,DY,EMX,EMY,NX,NY,R,T,DT CCX = ALF/DX/DX CCY = ALF/DY/DY

DO 1 I = DMX(I,1) DMX(I,2) DMX{I,3) DMY(1,I) DMY(2,I) DMY{3,I)

DO 5 J = DO 4 К = RTD = 0. DO 3 M MJ = J -DO 2 N NK = К -

= CCX*EMY(I) = -2.*DMX{I,1) = DMXd.l) = CCY*EMX{I) = -2.*DMY(1,I) = DMY(1,I)

2, NXN 2, NYN

1,3 2 + M 1,3 2 + N

2 RTD = RTD + (DMX(N,M)+DMY(N,M))*T(NK,MJ)*DTIK 1 + GAM*EMX(M)*EMY(N)*DT(NK,MJ)

3 CONTINUE

4 R(K,J) = RTD/(1. + GAM)

5 CONTINUE RETURN END

Рис. 8.5. Распечатка подпрограммы redif.

UNSTEADY HEAT CONDUCTION WITH NX,NY = 6 б

ME = 1 GAM = .00 BETA = .50 SX,SY = 1.000 1.000

APPROX. FACT., 3PT FDM, EMX= .OOOE+00 .lOOE+01 .OOOE+00

DIRICHLET AND NEUMANN B.C.

DTIM = .400E+04 PRINT INT. = .240E+05 TMAX .250E+O5

Рис. 8.6. Типовая выдача результатов, полученных с помощью программы

TWDIF.