130 DO 16 JM = 3,NXP

131 JMM = JM - 1

132 DUM = B(1,JM)/B(2,JMM)

133 B(2,JM) B(2,JM) - B(3,JMM)*DUM

134 B(3,JM) * B(3,JM) - B(4,JMM)*DU

135 B(1,JM) = 0.

136 RRT(JK) = RRT(JM) - RRT(JMM)*DUM

137 16 CONTINUE 133 B(l,2) 0.

139 С

140 17 CALL BANFAC(B,NXPД)

141 CALL BANSOL(RRT,DDT,B,NXP,1)

142 с

143 DO 18 J = 2,NX

144 JM = J - 1

145 18 R(K,J) = DDT(JM)

146 19 CONTINUE

147 с

148 С TRIDIAGONAL SYSTEMS IN THE Y-DIRECTIO№

149 С

150 DO 22 J = 2,NX

151 DO 20 К 2,NYP

152 KM = К - 1

153 B(2,KM) EMY(l) - CCYA - CYA*V(K ,J)

154 B(3,KM) EMY(2) + 2.*CCYA

155 B(4,KM) = EMY(3) - CCYA + CYA*V(K+1,J>

156 20 RRT(KM) = R(K,J)

157 B(2,l) = 0.

158 B(4,KM) 0.

159 С

160 CALL BANFAC(B,NYPP,1)

161 CALL BANS0L(RRT,DDT,B,NYPP,1)

162 С

163 С INCREMENT T

164 С

165 DO 21 К = 2,NYP

166 KM = К - 1

167 DT(K,J) DDT(KM)

168 21 T(K,J) = T(K,J) + DDT(KM)

169 22 CONTINUE

170 DO 23 К = 2,NYP

171 T(K,NXS) T(K,NXP)

172 23 DT(K,NXS) = DT(K,NXP)

173 ITER = ITER + 1

174 GAM = GAMH

175 BET = BETH

176 IFdTER .LE. ITMAX)GOTO 13

177 С

178 С COMPARE WITH SEMI-EXACT SOLUTION

179 С

180 24 SUM = 0.

181 DO 25 J = 2,NX

182 AJ = J - 1

183 X = AJ*DX

184 С

185 CALL TEXCL(X,TDX,PR,ALF,DYFL

186 С

187 TEX(J) = TDX

188 IF(J .EQ. 2)GOT0 25

189 DIF = T(NYH,J) - TDX

Рис. 9.13 (продолжение).

190 SUK SUM + DIF*DIF

191 25 CONTINUE

192 RMS = SQRT(SUM/(ANX-1))

193 WRITE(6,26)ITER,RMST,RMS

194 26 FORMATC AFTER M3, ITERATIONS, RMS-RHS =,E10.3,

195 1 RMS-ERR =,E10.3)

196 WRITE(6,27)(T(NYH,J),J=2,NX)

197 ,WRITE(6,28)(TEX(J),J=2,NX)

198 27 FORMATC C/L TEMP =M1F6.3)

199 28 FORMATC EX TEMP =M1F6.3)

200 STOP

201 END

Рис. 9.13 (окончание).

Применение группового метода конечных элементов (§ 10.3) [Fletcher, 1983] с линейной интерполяцией на прямоугольных элементах приводит к системе обыкновенных дифференциальных уравнений

Мх ® My []. = {М, ® [-L (иТ) + a.LxxT] +

+ Мх [-Ly (vT) + ауЬууТ]} ь (9.94)

которая отличается от (9.85) только тем, что в конвективные операторы Lx и Ly входят и, v. Следовательно, двухступенчатый алгоритм (9.88), (9.89) реализуется с помощью программы THERM (рис. 9.13) с операторами Lx и Ly, действующими на иТ и vT соответственно.

Программа THERM работает при любых заданных распределениях скорости. Однако для проведения сравнения расчетного профиля температуры с полуаналитическим решением [Brown, 1960] предполагалось, что профиль скорости имеет вид, характерный для полностью развитого течения:

1.5(1-у2), у = 0. (9.95)

Чтобы избежать разрыва решения в точке А, граничное условие в программе THERM задавалось в виде

Т{0, у) = у\ (9.96)

Это дает профиль температуры на входе, который равен нулю всюду, кроме угловых точек: л: = О, у = ± 1.

Вычисление величины RHS в (9.88) выполняет подпрограмма RETHE (рис. 9.14). Основные параметры, использованные в программе THERM и подпрограмме RETHE, описаны в табл. 9.11. Типичная выдача результатов программы THERM приведена на рис. 9.15 для конечно-разностного метода с приближенной факторизацией.

Решение на центральной линии сравнивается с полуаналитическим решением [Brown, 1960], которое вычисляется

27 К. Флетчер, т. I

418 Гл. 9. Линейные задачи с преобладающим влиянием конвекции 1

2 SUBROUTINE RETHE(RHST,GAM,DTIM)

4 С EVALUATES RIGHT-HAND SIDE OF THE

5 С TWO-DIMENSIONAL TRANSPORT EQUATION

7 REAL*8 RTD,SUMT,AN,DSQRT

8 DIMENSION DT{41,42),T(41,42),R(41,42),U(41,42),V(41,42),

9 1DMX{3,3),DMY{3,3),EMX{3),EMY(3),CX(3),CY(3),CXQ(4)

10 COMMON CX,CXQ,CY,CCX,CCY,EMX,EMY,NX,NY,R,T,DT,U,V

11 NXP = NX - 1

12 NYP = NY - 1

13 DO 1 I = 1,3 CCX*EMY(I) -2.*DMX(I,1) DMX(I,1) CCY*EMX(I) -2.*DMY(1,I) DMY(1,I)

14 DMX(I,1

15 DMX{I,2

16 DMX(I,3

17 DMY(1,I

18 DMY(2,I

19 1 DMY(3,I

20 С

21 SUNT = 0.

22 DO 6 J =2,NX

23 DO 5 К = 2,NYP

24 RTD = 0.

25 DO 4 N = 1,3

26 NK = К - 2 + N

27 DO 2 M = 1,3

28 MJ = J - 2 + M

29 DUM = DMX{N,M) + DMY(N,M) - CX(M)*EMY(N)*U(NK,MJ)

30 1 - CY(N)*EMX{M)*V(NK,MJ)

31 2 RTD = RTD + DUM*T(NK,MJ)*DTIM

32 1 + GAM*EMX(M)*EMY(N)*DT(NK,MJ)

33 MST = 1

34 IF(J .EQ. 2)MST = 2

35 DO 3 M = MST,4

36 MJ = J - 3 + M

37 3 RTD = RTD - CXQ(M)*EMY(N)*U(NK,MJ)*T{NK,MJ)*DTIM

38 IF(J .EQ. 2)RTD = RTD - CXQ(1)*EMY(N)*U(NK,1)*T(NK,1)*DTIH

39 4 CONTINUE

40 SUMT = SUNT + RTD*RTD

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

42 6 CONTINUE

43 AN = NXP*(NYP-1)

44 RMST = DSQRT(SUMT/AN)/DTIM

45 RETURN

46 END

Рис. 9.14. Описание программы RETHE.

С ПОМОЩЬЮ подпрограммы TEXCL (рис. 16.5, т. 2). Полуаналитическое решение строится для приведенной формы уравнения (9.90), в которой предполагается, что членом осевой диффузии €LxdT/dx можно пренебречь. Справедливость такого предположения обсуждается в п. 16.1.4. Следует заметить, что