+ - - + Ч--- +12 )-а

Б (9.73) uAx/Rceu = аАх. Видно, что для потоков с /?ceii 1 член с дисперсией дает наибольший вклад в ошибку аппрокси-

Рис. 9.7. Начальные условия в задаче о распространении температурного

фронта.

мации. Этот член можно исключить выбором = 0.5 -f 0.25С2. При таком выборе q и больших значениях Rceu диссипативный член имеет наименьший порядок и вносит положительную диссипацию.

Приведенные выше схемы (9.69) и (9.71), а также схемы, рассмотренные в п. 9.4.1 и 9.4.2, были использованы для решения задачи о распространении фронта температуры, иллюстрируемой на рис. 9.7. При / = О задан скачок температуры в точке X = 0. В последующие моменты времени фронт движется со скоростью и вправо и его профиль размывается под дей- ствием термодиффузии с коэффициентом а. Следовательно, в

Подстановка (9.72) в (9.71) приводит к четырехдиагональной системе уравнений, которая должна решаться на каждом шаге по времени. Для этого требуется использовать обобщенный алгоритм Томаса (п. 6.2.4). В программе TRAN обобщение алгоритма Томаса достигается за счет дополнительной прямой прогонки (строки 126-133), сводящей систему к трехдиагональной форме. В данном случае обобщенный алгоритм Томаса требует выполнения почти на 80 % операций больше, чем обычный алгоритм Томаса. Схема (9.71), (9.72) устойчива при q -3/jRceii. Практически это условие не является ограничением, поскольку обычно интерес представляют только положительные значения q.

Эквивалентное уравнению (9.71) модифицированное уравнение (п. 9.2.2) имеет вид

Рис. 9.8. Описание программы TRAN (начало).

2 С TRAN SOLVES THE LINEAR TRANSPORT EQUATION USING

3 С VARIOUS EXPLICIT AND IMPLICIT SCHEMES

5 DIMENSION R(65),T(65),TD(65),TEX(65),A(5,65)

6 1,D(4),X(65)

7 OPEN(1, FILE= TRAN.DAT)

8 OPEN(6, FILE= TRAN.OUT)

9 READ(1,1)KE,JMAX,NTIM,NEX,JPR,LN

10 READ{1,2)C,U,S,Q,EM

11 1 FORMAT(615)

12 2 FORMAT(2F5.2,3E10.3)

13 С

14 IF(JPR .EQ. 1)WRITE(6,3)

15 IF(JPR .EQ. 2)VRITE(6,4)

16 3 FORMATC PROPAGATING SINE-WAVE)

17 4 FORMATC PROPAGATING TEMP-FRONT*)

18 JMAF = JMAX - 2

19 JMAP = JMAX - 1

20 AMP = JMAP

21 DX = 1.0/AMP

22 IF(JPR .EQ. 2)DX = 4./АИР

23 DT = C*DX/U

24 EL = LN

25 ALPH = S*DX*DX/DT

26 ir(ALPH .LT..1.0E-10)ALPH l.OE-10

27 RCEL = U*DX/ALPH

28 QQ = Q*C/3.

29 IF(ME .LT. 3}QQ = 0.

30 MQ = 0

31 IF(ABS{QQ) .GT. O.OOODMQ 1

32 ATIM = NTIM

33 TIM = 0.

34 TIMAX = DT ATIM

35 PI = 3.141592654

36 С

37 IF(ME .EQ. 1)WRITE(6,5)ME

38 IF(ME .EQ. 2}VRITE(6,6)ME

39 IFCME .EQ. 3)VRITE(6,7)ME

40 IF(ME .EQ. 4}WRITE{6,8)ME

41 5 FORMATC ME =M2, FTCS DIFFERENCING*)

42 6 FORMATC ME =M2, LAX-VENDROFF*)

43 7 FORMATC ME =M2, EXPLICIT 4PT UPWIND)

44 8 FORMATC ME =M2, GENERAL CRANK-NICOLSON*)

45 WRITE(6,9)JMAX,NTIM,C,U,DX,DT

46 9 FORMATC JMAX=M3, NTIM=M3/ C=\F5.2/ U=\F5.2

47 1* DX=\F5.3, DT=,F5.3)

48 WRITE(6,10)S,ALPH,RCEL,Q,QQ

49 WRITE(6,11)NEX,LN,EM

50 10 FORMATC S=,F5.2, ALPH=,E10.3, RCEL=. r?6.3,

51 1* Q=,F5.2, QQ=\F6.3)

52 11 FORMATC- NEX=M5, EL=M5,* EM=\E10.3)

IF(KE .GT. 3)G0T0 12 SS = S

IF(KE .EQ. 2)SS = S + 0.5*C*C AA = (0.5*C + SS) + 3.*0Q BB = 1. - 2.*SS - 3.*QQ CC = -0.5*C + SS + QQ GOTO 13

12 AA = EM - 0.25*C - 0.5*S - 1.5*QQ BB 1.0 - 2.0*EM + S + 1.5*QQ

CC = EK + 0.25*C - 0.5*S - 0.5*QQ AE ЕИ +0.25*C + 1.5*QQ + 0.5*S BE 1. - 2.*ЕИ - 1.5*QQ - S CE = EM - 0.25*C + 0.5*00 + 0.5*S

13 VRITE(6,14)AA,BB,CC,AE,BE,CE

14 FORKATC AA=,F8.5, BB=,F8.5, CC=4F8.5, AE=\F8.5, BE 1F8.5, CE=,F8.5,/)

INITIALISE T AND EVALUATE TEX

CALL EXSOL{JPR, JMAX, X, T, TEX,NEX, DX,U, ALPH,TIMAX, EL)

DO 16 J = 1,JMAX DO 15 К = 1,5

15 A{K,J) f 0.

16 CONTINUE VRITE(6,17)TIM

17 FORMATC INITIAL SOLUTION, TIM =\F5.3) VRITE(6,18)(X(J),J=1,JMAX) VRITE(6,19)(T(J),J=1,JMAX)

18 FORMATC X=M2F6.3)

19 FORMATC T=M2F6.3)

MARCH SOLUTION IN TIME

DO 26 N = 1,NTIM

IF(ME .GT. 3)G0T0 21

EXPLICIT SCHEMES

D(l) = T(l)

D(2) = T(l)

D(3) = T(2)

DO 20 J = 2,аМАР

IF(ME .EQ. 3)D(4) = D(l)

D(l) = D(2)

D(2) = D{3)

D(3) = T(J+1)

T(J) = AA*D(1) + BB*D(2) + CC*D(3) IF(ME .EQ. 3)T(J) = T(J) - QQ*D(4) 20 CONTINUE GOTO 26