1 с BIFEX SOLVES THE DIFFUSION (ID TRANSIENT HEAT CONDUCTION)

2 С EQUATION USING VARIOUS EXPLICIT SCHEMES

4 DIMENSION TN(41) ,DUH(41) ,TD(41) Д (41) Д£(41) 0L(41) EL(3)

5 REAL*8 SUM,AVS,RMS,DSQRT,DMP

7 OPEN(1, FILE* DIFEX.DAT) $ OPEN(б,FILE * DIFEX.OUT)

$ 1tEAD(l,l)ME,IPR.ОМАХ,MAXEX,NMAX,ALPH,S,TMAX,TST,GAM

10 1 rORHAT(5I5,E10.3,F5.3,3F5.2)

11 С

12 PI 3.14Ш27

13 JNAP = JMAX - 1

14 AJM = JMAP

15 DELX ж l./AJM

16 DELT = DELX*DEU*S/ALPH

17 С

18 ir(MC .EQ. 1)VRITE(6,2)

19 ir<ME .EQ. 2}VRITE(6,3)

20 ir(KE .EQ. 3)VRITE(6,4).

21 2 FORHATC FTCS SCHEME*, )

22 3 FORMATC DUFORT-FRANKEL SCHEME, )

23 4 FORMATC 3-LEVBL, 4TH-0RDER SCHEME, )

24 VRITE(6,$)JMAX,MAXEX,NMAX,TMAX,TST,GAK

25 5 FORMATC JMAX-M5,* MAXEX=M5,* NMAX-M5,

26 1 TMAX *,F5.2, TST ,F5.2, GAM-4F5.2)

27 VRITE< 6,6)S,ALPH,DELT,DELX

2$ 6 FORHATC S=\F5.3,* ALPH = ЧЕЮ.З, DELT ЧВЮ.З,

29 ! DELX = ЧБЮ.З, )

30 С

31 IF(NE .EQ. 2)AS 2.*S/(l.-l>2.*S)

32 IFCME .EQ. 2)BS (1. - 2.*S)/(1. 4- 2.*S)

33 IF(ME .NE. 3)G0T0 7

34 AS (l.+2.*GAM)/(l.+GAM)

35 as = -GAM/Cl.+GAM)

36 CS = S/d.-bGAN)

37 DS = 1.5 + GAM - 1./12./S

38 ES = -0.5 - GAM + 1./12./S

39 EL(1) = 1.

40 EL(2) = -2.

41 EL(3) = 1.

42 7 CONTINUE

43 С

44 С OBTAIN INITIAL CONDITIONS FROM EXACT SOLUTION

45 С

46 T = TST - 2.*DELT

47 DO 9 I 1,2

48 T = T + DELT

49 С

50 CALL EXTRA(JMAX,MAXEX,DELX,PI,ALPH,T,TE)

51 С

52 DO 8 J s 2,JMAP

53 IF(I .EQ. l)TOL(J) TE(J)/100.

54 IF(I .EQ. 2)TN(J> = TE(J)/100.

55 8 CONTINUE

56 9 CONTINUE

57 С

58 С SET BOUNDARY CONDITIONS

59 С

Рис. 7.4. Распечатка программы DIFEX (начало.)

60 TOL(l) = 1.

61 TOL(JMAX) 1.

62 TN(1) = 1.

63 TN(JMAX) = 1.

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

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

66 N 0

67 С

68 С EACH TIME STEP STARTS AT STATEMEKT 10

69 С

70 10 N N + 1

71 IF (ME .EQ. DGOTO 15

72 IF(ME .EQ. 2)GOTO 13

73 С

74 С ME 3, 3-LEV., 4TH ORDER

75 С

76 DO 12 J =2,JMAP

77 DUM(J) AS*TN(J) + BS*TOL(J)

78 DO 11 К 1,3

79 KJ J - 2 + К

80 DUM(J) = DUM(J) CS*EL(K)*(DS*TN(KJ) ES*TOb(KJ)>

81 11 CONTINUE

82 12 CONTINUE

83 GOTO 17

84 С

85 С ME = 2, DUFORT-FRANKEL

86 С

87 13 DO 14 J 2,JMAP

8$ 14 DUM(J) = AS*(TN(J-1) + TN(J+1)) + BS*TOL(J)

89 GOTO 17

90 С

91 С ME = 1, FTCS

92 С

93 15 DO 16 J = 2,JMAP

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

95 16 CONTINUE

96 17 DO 18 J = 2,JMAP

97 IF(ME .GT. l)TOL(J) = TN(J)

98 18 TN(J) = DUM(J)

99 С

100 DO 19 J = 2,JMAP

101 19 TD(J) = 100.*TN(J)

102 T = T + DELT

103 IFdPR .EQ. l)WRITE(6,20)T,(TDr(J),J=l,JMAX)

104 20 FORMATC T= ,F5.2, TD=M1F6.2)

105 С

106 С IF MAXIMUM TIME OR MAXIMUM NUMBER OF TIME-STEPS

107 С EXCEEDED EXIT FROM LOOP

108 С

109 IF(N .GE. NMAX)GOTO 21

110 IF(T .LT. TMAX)GOTO 10

111 С

112 С OBTAIN EXACT SOLUTION AND COMPARE

113 С

114 21 SUM = 0.

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

115 С

116 CALL EXTRA (JMAX, МАХЕХ, DELX, PI, ALPH ДДЕ)

117 С

118 DO 22 J = 1,JMAX

119 DMP TE(J) - TD(J)

120 SUM SUM DMP*DMP

121 22 CONTINUE

122 IFdPR .NE. 1)WRITE(6,20)T,(TD(J),J=1,JMAX)

123 WRITE(6,23)T,(TE(J),J 1,JMAX)

124 23 FORMAT(/, T ,F5.2, TE M1F6.2, )

125 С

126 С RMS IS THE RMS ERROR

127 С

128 AVS SUM/(1. + AJM)

129 RMS DSQRT(AVS)

130 WRITE(6,24)RMS

131 24 FORMATC RMS DIF ,D11.4, )

132 STOP in END

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

в тексте программы DIFEX, распечатка которой приведена на рис. 7.4 и которая является обобщением программы DIFF (рис. 3.13).

Решения ищутся в вычислительной области 0x1.0 и 2.00 f 9.00 при начальных условиях, заданных при t = = 2.00, с помощью точного решения (3.42), деленного на 100. Граничные условия имеют вид

7(0, t) = T{l, 0=1.0. (7.18)

Представление о точности различных схем можно получить путем оценки среднеквадратичной разности между численным и точным решениями при Т = 9.00. Точное решение рассчитывается с помощью подпрограммы EXTRA (рис. 7.5). Схема Дюфорта - Франкела и трехслойная схема нуждаются в задании начальных данных на двух слоях при t = 2.00 и t = 2.00 - А.

Основные параметры, используемые в программе DIFEX (см. рис. 7.4), приводятся в табл. 7.2. Характерная форма представления решения, построенного с помощью трехслойной схемы на сравнительно грубой сетке, показана на рис. 7.6.

Точности различных вариантов решения, полученных при использовании программы DIFEX, сравниваются в табл. 7.3. Чтобы получить возможность сравнения трехслойной схемы четвертого порядка в виде (7.13) со схемами ВВЦП и Дюфорта- Франкела, использовалось два значения 5: 5 = 0.3 и 0.41. Кроме того, в таблицу включены два частных случая - ВВЦП с 5=1/6 и схема Дюфорта - Франкела с s =(1/12)2. Оба