105 С

106 С TRIDIAGONAL SYSTEM FOR IMPLICIT SCHEMES

107 С

108 21 IF(MQ .EQ. DDIM = 1.

109 DO 22 J = 2,JMAP

110 ,:л = J - 1

111 A(1,JM) = 0.5*QQ

112 A(2,JH) = AA

113 A(3,JM) = BB

114 A(4,JM) = CC

115 R{JM) = A£ T(JM) + BE*T{J) + CE*T(J+1)

116 IF(MQ .EQ. 1 .MD. J .GT. 2)DIM = T(JM-l)

117 IFdlQ .EQ. l)R(JM) = R{JM) - 0.5*QQ*DIM

118 22 CONTINUE

119 R(l) = R{1) - A(2,1)*T(1)

120 IF(MQ .EQ. 0)GOTO 24

121 R(l) = R(l) - A(1,1)*T(1)

122 R(2) = R(2) - A(1,2)*T(1)

123 С

124 С REDUCE A TO TRIDIAGONAL FORM

125 С

126 DO 23 JM = 3,JMAF

127 JMM = JM - 1

128 DUM = A(1,JM)/A{2,JMM)

129 A(2,JM) = A(2,JM) - A(3,JMM)*DUM

130 A(3,JM) = A(3,JM) - A(4,JMM)*DUM

131 Ad.JM) = 0.

132 R(JM) = R(JM) - R(JMM)*DUM

133 23 CONTINUE

134 А(1Д) = 0.

135 A(l,2) = 0.

136 24 А(2Д) = 0.

137 A(4,JMAF) = 0.

138 С

139 CALL BANFAC(A,JMAF,1)

140 CALL BANSOL(R,TD,A,JMAF,1)

141 С

142 DO 25 J = 2,JMAP

143 25 T(J) = TD(J-l)

144 26 CONTINUE

145 С

146 WRITE(6,27)TIMAX

147 27 FORMATC FINAL SOLUTION, TIM =,F5.3)

148 WRITE(6,18){X(J),J=l,JMAX)

149 WRITE(6,19){T(J),J=1,JKAX)

150 WRITE(6,28)(TEX(J),J=1,JMAX)

151 28 FORMATC TEX=M2F6.3)

152 SUM = 0.

153 DO 29 J = 2,JMAP

154 29 SUM = SUM + (T(J) - TEX(J))**2

155 RMS = SQRT(SUM/(AMP-1.))

156 VRITE(6,30)RMS

157 30 FORMATC RMS ERR=M12.5)

158 STOP

159 END

SUBROUTINE EXSOL(JPR, JMAX,X,Т, TEX,NEX, DX,U, ALPH, TIMAX, EL)

SETS THE INITIAL T SOLUTION AND FINAL EXACT (TEX) SOLUTIOK

DIMENSION T(65),TEX(65),X(65) JMAP = JMAX - 1 PI = 3.141592654 IF(JPR .EQ. DXST = 0. IF(JPR .EQ. 2)XST = - 2.0 DO 1 J = l.JMAX AJ = J - 1 X(J) = XST + AJ*DX T(J) = 0.

1 TEX(J) = 0. IF(JPR .EQ. 2)GOTO 3

EXACT SOLUTION FOR PROPAGATING SINE-VAVE

JM = 0.1001/DX + 1.0 INC = U*TIMAX/DX + 0.001 DO 2 J = 1,JM T(J) = SIN(10.*PI*X(J)) JP = J + INC

2 TEX(J?) = T(J) RETURN

EXACT SOLUTION FOR PROPAGATING TEMP-FRONT

3 T(l) = 1.0 TEX(l) = 1.0 DO 5 J = 2,JMAP IF(X(J) .LT. O.)T(J) = 1.0 IF(ABS(X(J)) .LT. 1.0E-04)T(J) = 0.5 DEM = 0. DO 4 К = 1,NEX AK = 2*K - 1 DUM = AK*PI/EL

SNE = SIN(DUMMX(J)-U*TIMAX)) DIM = - ALPH*DUM*DUM*TIMAX IF(DIM .LT. -20.)GOTO 5

4 DEM = DEM + (SKE/AK)*EXP(DIM)

5 TEX(J) = 0.5 - 2.*DEM/PI RETURN END

Рис. 9.9. Описание программы EXSOL.

каждый момент времени чем больше jRceii = uAx/a, тем круче профиль температурного фронта.

Исходным уравнением в этой задаче является уравнение (9.56). Для достаточно малых значений t можно задать следующие граничные условия:

Т (~2, О =1.0, Т (2, О = 0.0. (9.74)

Т (jc, О = 0.5 ~ sin [(2k ~ 1)

п(х- Ut)

ехр[-а(2А;- Оя/!]

2k-\

(9.75)

Различные схемы, которые заложены в программу TRAN, приведены в табл. 9.4. Распечатка программы TRAN дана на рис. 9.8, а основные параметры указаны в табл. 9.5. Точное решение дает подпрограмма EXSOL (рис. 9.9). Типичная выдача результатов программы TRAN при использовании схемы Лакса - Вендроффа показана на рис. 9.10.

PROPAGATING TEMP-FRONT HE 2 LAX-WENDROFF

JMAX 21 NTIM= 10 C .25 U= .50 DX .200 DT .100 S= .25 ALPH= .lOOE+00 RCEL= 1.000 Q= .00 QQ= .000 HEX: 100 EL= 20 ENs .OOOE-i-OO

ЛА= .40625 BB= .43750 CC= .15625 AE> .00000 BE .00000 CE= .00000

INITIAL SOLUTION, TIM .000

X -2.000-1.800-1.600-1.400-1.200-1.ООО -.800 -.600 -.400 -.200 .000 .200

X= .400 .600 .800 1.000 1.200 1.400 1.600 1.800 2.000

T 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 .500 .000

T .000 .000 .000 .000 .000 .000 .000 .000 .000 riNAL SOLUTION, TIM 1.000

X=-2.000-1.800-1.600-1.400-1.200-1.ООО -.800 -.600 -.400 -.200 .000 .200

X= .400 .600 .800 1.000 1.200 1.400 1.600 1.800 2.000

T= 1.000 1.000 1.000 1.000 1.000 .999 .997 .991 .973 .934 .862 .747

T .593 .422 .261 .137 .059 .020 .005 .001 .000

TEX= 1.000 1.000 1.000 1.000 1.000 1.000 .998 .993 .978 .941 .868 .749

ТЕХ .588 .412 .251 .132 .059 .022 .007 .002 .000 KHS ERR .44842E-02

Рис. 9.10. Типичная выдача по программе TRAN.

При /?се11=1.0 В табл. 9.6 дано сравнение решений, полученных как по явным, так и по неявным схемам. Все методы приводят к гладким решениям, тогда как схема Лакса - Вендроффа оказывается наиболее точной из явных схем. Наиболее точной из неявных схем является схема (9.71) с = 0.5, что соответствует 3-му порядку аппроксимации члена идТ/дх.

Расчеты по неявным схемам Кранка - Николсона для более высоких значений jRceii показаны в табл. 9.7 (Rceu = 3.33) и табл. 9.8 (/?ceii= 100). В качестве эталона принята конечно-разностная схема Кранка - Николсона. Решение по этой схеме является осциллирующим и выявляет добавочный скачок за фронтом, который особенно хорошо виден при Rceu = 100 (табл. 9.8 и рис. 9.11). Видно, что две низкодисперсные схемы

Точное решение этой задачи, найденное методом разделения переменных, имеет вид ы