60 CALL BANrAC(B,NYPP,l)

61 С

62 X = 0.

63 XPR = 0.

64 NCT = О

65 SUMT = 0.

66 7 нет = NCT + I

67 с

68 с GENERATE R.H.S.

69 с

70 DO 8 К = 2,NY?

71 КМ = К - 1

72 КР = К + 1

73 8 R(KM) = ССАМТ(КМ)-2.*Т{К)+Т(КР)) - CA*(V(KP)*Т(КР)-У{КН)*Т(ИС))

74 С

75 CALL BANS0L(R,DT,B,NYPP,1)

76 С

77 DO 9 К = 2,NYP

78 9 Т(К) = Т(К) + DT{K-1)

79 X = X + DX

80 С

81 С EXACT С/И SOLUTION

82 С

83 CALL TEXCL(X,TEX,PR,ALF,DYFL)

84 с

85 DMP = T(NYH) - ТЕХ

86 IF(NCT .GT. 2)SUMT = SUMT + DMP*DMP

87 IF(X .LT. XPR)GOTO 11

88 WRITE(6,10)X,(T(K),K=1,NYH),TEX

89 10 FORMATC X=\F4.2, T=,6F6.3,* ТЕХ-,Гб,3)

90 XPR = XPR 4- DXP - 0.0001

91 llIFiX :gE. XMAX)G0T0 12

92 IF(NCT .GE. NXMAX)G0T0 12

93 GOTO 7

94 12 ANCT = NCT - 2

95 RMS p SQRT(SUMT/ANCT)

96 WRITE(6,13)NCT,RMS

97 13 FORMATC NCT=M5/ RMS=\E10.3>

эа 14 STOP

99 END

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

Полученное решение вдоль центральной линии (у = 0) сравнивается с полуаналитическим решением Брауна [Brown, 1960]. Браун получил разделение переменных в уравнениях (16.43) и (16.48), основанное на экспоненциально затухающем в направлении X решении и разложении по собственным числам/собственным (функциям по у. Первые десять членов решения Брауна на центральной линии вычисляются в подпрограмме TEXCL (рис. 16.5). В результате работы этой подпрограммы определяется точное решение ТЕХ.

Типичное решение при Дх = 0.05 и Ду = 0.2, полученное по программе THRED, приведено на рис. 16.6. Рассчитанное распределение температуры симметрично относительно у = О, поэтому значения температуры приведены лишь в области -1

у 0. Крайний правый столбец температур Т соответствует значениям на центральной линии (у = 0), и значения в нем можно сравнить с полуаналитическими значениями ТЕХ. На сравнительно грубой сетке в решении, приведенном на рис. 16.6, заметны осцилляции вблизи точки а: О, у -1.0. Эти осцилляции связаны с быстрым изменением Т в граничных условиях

вблизи точек (О, ±1). Амплитуда осцилляции уменьшается с увеличением х. На более мелких по х или у сетках осцилляции не возникают.

Можно заметить, что в решениях, например AF-FEM, полученных по программе THERM (рис. 9.13), на сетке 11X11 (Ах = 0.20, Ау = 0.20) нет существенных осцилляции вблизи точек (О, dzl). Однако программа THERM основана на решении уравнения (9.90), в которое входит член дТ/дх, Этот член обладает сглаживающим свойством и противодействует появлению осцилляции вблизи точек (О, ±1).

Решение на центральной линии (RED-FEM), полученное по программе THRED, сравнивается в табл. 16.3 с решением, полученным по программе THERM, и с полуаналитическим

1 2 3 4

5 б 7 8 9 10 11 12 13 14 15 16 17 18 19

SUBROUTINE TEXCL(XДЕХ, PR, ALF,DYFL)

FOR GIVEN X AND PR COMPUTE EXACT CENTRE-LINE TEMPERATURE DISTRIBUTION

DIMENSION ALF(10),DYFL(10)

ZD -3.2*X/PR/3.0

ТВ 0.

DO 1 I 1,10

DUM ZD*ALF(I)*ALr(I)

ir(DUM .LT. -20.)GOTO 1

DUM EXP(DUM)

CF = -2./ALF(I)/DYFL(I)

TB = ТВ + CF*DUM

CONTINUE

TEX = 1. - TB

RETURN

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

REDUCED THERMAL ENTRY PROBLEM BY C.N.-FEM

НУ 11 NXMAX= 50 DX= .500E-01 XMAX= 2.000 PR .700

X .05 X .20 X .40

X .60 X .80 X=1.00 X 1.20 X 1.40 X 1.60 X=1.80 X 2.00 NCT

t 1.000 t 1.000 t= 1.000 t= 1.000 t 1.000 t 1.000 t= 1.000 t 1.000 t= 1.000 t= 1.000 t 1.000 41 rms=

1.076 .701 .868 .937 .969 .984

.477 .767 .917 .972 .992 .999

.992 1.001 .995 1.001 .997 1.001 .999 1.001 .Э99 1.000 .350E-02

.188

.614

.827

.924

.966

.985

.993

.997

.998

.999 1

.000 1

.074 .046

.520 .495

.801 .793

.918 .913

.967 .963

.987 .984

.995 .993

.998 .997

.999 .999

.000 .999

.000 1.000

tex tex= tex= tex tex= tex tex tex tex tex tex

.058 .493 .786 .909 .962 .984 .993 .997 .999 .999 1.000

Рис. 16.6. Типичная выдача программы THRED.

Таблица 16.3. Решение на центральной линии для задачи ввода тепла.

At/ =0.20

0.000

0.200

0.400

0.600

0.800

1.000

1.200

1.400

1.600

1.800

2.000

Среднеквадратичное отклонение

<Полу)-точное AF-FEM, Ал: = 0.20 RED-FEM, Ах = 0.05 RED-FEM, Лл: = 0.010

0.000 0.000 0.000 0.000

0.493 0.462 0.495 0.497

.794 0.

0.786 0. 0,793 0.789

0.910 .910 0.913 0.912

0.962 0.963 0.963

,963 О

0.984 0.984 0.984 1.984

0.993 0.994 0.993 0.993

,997 О

,997 О

0.997 0.997

1.999 1.999 0.999 0.999

1.000 0.999 1.000 1.000

1.000 1.000 1.000 1.000

0.003

0.0035

0.0019