Для решения этого уравнения необходимы начальные условия (3.3), роль которых состоит в задании начальных значений за-

Тп(х)

Задано Тд или 9Тд/Эх

ч\Х\<\\\\ Изоляционный слой \N\SSS

Задано Tg или дТ/дх

1=0 х=1

Рис. 3.12. Нестационарная теплопроводность в стержне.

висимого переменного (см. рис. 3.2). Один из простейших конечно-разностных аналогов уравнения диффузии достигается путем замены производной по времени аппроксимацией с разностью вперед, а производной по пространству - центрированной разностной аппроксимацией (см. п. 3.1.1). Это приводит через посредство формулы (3.5) к схеме с разностями, сдвинутыми вперед по времени и с центрированными по пространству (ВВЦП)

rsTU + (1 - 25)+ sT%u (3.41)

где S = aAx/At. Формула (3.41) применяется ко всем внутренним узлам, где / = 2, 7-1. В случае типичной задачи о теплопроводности граничные значения ГГ и Г/ задаются лри помощи граничных условий (3.2). Процесс решения повторяется по мере продвижения во времени (д=1, 2, ...), пока не достигается необходимое время окончания.

3.5.2. DIFF: нестационарная задача теории теплопроводности (диффузии)

В п. 3.5.1 было дано качественное описание реализации метода конечных разностей. Здесь будет описана соответствующая программа расчетов на ЭВМ под названием DIFF.

= 0°С. при = 0 два горячих резервуара (Г=100°С) приводятся в контакт с двумя концами стержня А и В. Задача состоит в том, чтобы определить (численно) температуру T(xJ) любой точки стержня в любой последующий момент.

Исходным уравнением для этой задачи является уравнение (3.1), т. е.

дТ дТ

2 С SOLVES ID TRANSIENT HEAT CONDUCTION EQUATION USING

3 с FTCS SCHEME,

5 DIMENSION TN(41),DUM(41),TD(41),X(41),TE(41)

7 С INPUT DATA;

9 С JMAX = THE NUMBER OF POINTS ALONG THE ROD

10 С MAXEX = THE NUMBER OF TERMS IN THE EXACT SOLUTION

11 С NMAX = THE MAXIMUM NUMBER OF TIME STEPS

12 С ALPH = THE THERMAL DIFFUSIVITY

13 С S = ALPH*DELT/DELX/DELX

14 С TMAX = THE MAXIMUM TIME

15 С

16 0PEN(1,FILE=DIFF.DAT)

17 OPEN(6,FILE= DIFF.OUT)

18 READ(1,1)JMAX,MAXEX,NMAX,ALPH,S,TMAX

19 1 FORMAT(3I5,E10.3,F5.2.F5.0)

20 PI = 3.1415927

21 С

22 С TD = DIMENSIONAL TEMPERATURE

23 С TN = NONDIMENSIONAL TEMPERATURE

24 С

25 JMAP = JMAX - 1

26 AJM = JMAP

27 DELX = l./AJM

28 DELT = DELX*DELX*S/ALPH

29 WRITE(6,2)JMAX,MAXEX, NMAX,TMAX

30 2 FORMATC JKAX=M5, MAXEX=M5, NMAX=M5, TMAX=,F8.2)

31 WRITE(6, 3)S, ALPH,DELT,DELX

32 3 FORMATC S=,F5.3, ALPH=4E10.3, DELT=,E10.3,

33 1 DELX ,E10.3, )

34 WRITE(6,4)S

35 4 FORMATC FTCS(EXPLICIT) SCHEME,5X, S =,F5.3, )

36 с

37 С SET INITIAL CONDITIONS

38 С

39 DO 5 J = 1,JMAP

40 5 TN(J) = 0.

41 N = 0

42 T = 0.

43 С

44 С SET BOUNDARY CONDITIONS

45 С

46 С EACH TIME STEP STARTS AT STATEMENT 2

47 С

48 6 TN(1) = 1.

49 TN(JMAX) = 1.

50 IF(T .LT. O.Ol)TN(l) = 0.5

51 IF(T .LT. O.Ol)TN(JMAX) = 0.5

52 TD(1) = 100.*TN(1)

53 TD(JMAX) = 100.*TN(JMAX)

54 С

55 С COMPUTE r.D. SOLUTION

56 С

57 DO 7 J = 2,JMAP

58 DUM(J) = (l.-2.*S)*TN(J) + S*(TN(J-1) + TN(J+1))

59 7 CONTINUE

60 DO 8 J = 2,JMAP

61 8 TN(J) DUM(J)

62 С

63 DO 9 J = 2,JMAP

64 9 TD(J) = 100.*TN(J)

65 T = T + DELT

66 WRITE(6,10)T,(TD(J),J=1,JMAX)

67 10 FORMAT С T= ,F5.0, TD Mir6.2)

68 С

69 С IF MAXIMUM TIME OR MAXIMUM NUMBER OF TIME-STEPS EXCEEDED

70 С EXIT FROM LOOP

71 С

72 IF(W .GE. NMAX)GOTO 11

73 IF(T .LT. TMAX)GOTO 6

74 С

75 С OBTAIN EXACT SOLUTION AND COMPARE

76 С

77 11 SUM = 0.

78 DO 13 J 1,JMAX

79 AJ = J - 1

80 X(J) DELX*AJ

81 TE(J) = 100.

82 DO 12 M lMAXEX

83 AM M

84 DAM (2.*AM - 1.)

85 DXM DAM*PI*X(J)

86 DTM = -ALPH*DAM*DAM*PI*PI*T

87 С

88 С LIMIT THE ARGUMENT SIZE OF EXP(DTM)

89 С

90 IFCDTM .LT. - 87.)DTM = -87.0

91 12 TE(J) = TE(J) - 400./DAM/PI*SIN(DXM)*EXP(DTM)

92 SUM = SUM + (TE(J) - TD(J))**2

93 13 CONTINUE

94 WRITE(6,14)T,(TE(J),J=1,JMAX)

95 14 FORMAT( T=,F5.0/ TE=M1F6.2, )

96 с

97 С RMS IS THE RMS ERROR

98 С

99 AVS = SUM/(1. + AJM)

100 RMS SQRT(AVS)

101 WRITE(6,15)RMS

102 15 FORMATC RMS DIF = ЧЕ11.4, )

103 STOP

104 END

Рис. 3.13. Распечатка программы DIFF. 7 К Флетчер, т. I