58 С

59 С ITERATE USING SOR

60 С

61 DO 13 I i;iTHX

62 SUM 0.

63 DO 12 J 2,NXP

64 JM J - 1

65 JP J + 1

66 DO 11 К 2,NYP

67 KM = К - 1

68 KP К + 1

69 IF(ME .EQ. 2)G0T0 9

70 DUM = 1. + PAR1/6.*(VA(JM,KP) + VA(JP,KP) + WA(JM,KM) + VA(JP,KM)>

71 DUM = DUM + PAR2*(VA(J,KP) + WA(J,KM)) + PAR3* (WA( JM,K) +VA(JP,K))

72 VDD 0.75*DUM/PAR1

73 GOTO 10

74 9 DUM 1. + BAS*(VA(JM,K) + WA(JP,K)) + (VA(J,KM) + VA(J,KP))/DYS

75 VDD 0.5*DUM/PAR1

76 10 Dir VDD - VA(J,K)

77 SUM SUM + Dir*DIF

78 VA(J,K) VA(J,K) + OM*Dir

79 11 CONTINUE 30 12 CONTINUE

81 RMS = SQRT(SUM/(ANX-1.)/(ANX-1.))

82 IF(RMS .LT. EPS)G0T0 15 3 13 CONTINUE

84 VRITE(6,14)ITMX,RMS

85 14 FORMATC CONVERGENCE NOT ACHIEVED IN* ,15/ STEPS ,5X,

86 1* ITER RMS ,E12.5)

87 с

B8 С COMPARE SOLUTION VITH EXACT 89 с

SO 15 SUM 0.

91 QL 0.

92 DO 18 J 2,NXP

93 DO 16 К 2,NYP

94 DIF VA(J,K) - VEX(J,K)

95 SUM SUM + DIF*DIF

96 QL QL + VA(J,K)

97 16 CONTINUE

98 IFdPR .EQ. OGOTO 18

99 VRITE(6,17)(VA(J,K),K 2,NYP)

100 17 FORMATC VA=M1F7.4)

101 18 CONTINUE

102 QL = QL*DX*DY

103 RMS SQRT(SUM/(ANX-1.)/(ANX-1.))

104 VRITE(6,19)I,RMS

105 VRITE(6,21)VA(NH,NH),VEX(NH,NH)

106 VRITE(6,20)QL,QLX

107 19 FORMAT(/, CONVERGED AFTER M3, STEPS,5X,SOL RKS ,E12.5, )

108 20 FORMATC COMP. FLOV RATE ,£12.5, EXACT FLOV RATE ,E12.5)

109 21 FORMATC V(C/L) ,E12.5, VEX(C/L}*\E12.5, )

110 STOP

111 END

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

DUCT считывает (и выводит на печать) различные управляющие параметры (строки 10-21). Задаются точное решение и его начальная форма (строки 37-51). При решении данной задачи число итераций по схеме ПВР не представляет интереса, поэтому в качестве начальной формы решения по ПВР применялось точное решение. Подпрограмма VEL (строка 47 на

2 3 4 5 б 7 8 9 О 1 2 3 4 5 6 7 8 9 20 21 22 23 24 25 26 27 28 29 30

SUBROUTINE VEL(X,У,N,PI,BAR,W,QL)

FOR GIVEN X,Y COMPUTES AXIAL VELOCITY IN CHANNEL AND COMPUTES THE FLOW RATE

CON = (8./PI/PI)**2 XAN 0.5*PI*X YAN 0.5*PI*Y W = 0.

QL DO AI CX DO AJ CY IJ IFdJ

I * 1,N,2 I

COS(AI*XAN) J = 1,N,2 J

COS(AJ*YAN) (I+J)/2 - 1 .EQ. 0}SIG

IFdJ .GT. 0)SIG

(-l.)**IJ

DEN AI*AJM(AI*BAR)**2 ♦ AJ*AJ) DUM = SIG*CX*CY/DEN W = W + DUM QL QL + l./DEN/AJ/AI

1 CONTINUE

2 CONTINUE W CON*W

QL = (16./PI/PI)*CON*QL RETURN END

Рис. 5.23. Распечатка подпрограммы VEL.

рис. 5.22 и рис. 5.23) вычисляет точное решение, получаемое по методу разделения переменных (см. п. 2.5.2), и имеющее вид

-Ы) L L [ц(1Ь/аГ + Р со8(0.5/ял:)со8(0.5/я)),

t = I,3, 5. ... / = 1,3, 5, ...

(5.110)

NEX NEX

i = l,3, 5, ... / = 1. 3, 5, ...

где полный расход потока равен wdxdy.

Далее программа DUCT (строки 58-83) ведет счет по формулам (5.108) и (5.109) с помощью итераций, пока величина I йу* - 1 не станет меньше заданной малой величины

EPS. После этого (строки 87-106) решение сравнивается с точным решением и выводится на печать. Основные параметры, используемые в программе DUCT, приводятся в табл. 5.11.

Таблица 5.11. Параметры, используемые в программе DUCT

Таблица 5.12. Изменение точности решений по программе DUCT по мере измельчения сетки {Ь/а=\)