2 3 4 5 б 7 6 9 le 11 12 13 14 15 16 17 19 19 29 21 22 23 24 25 26 27 28 2S Зв 31

33 34 35

37 38 39 40 41 42 43 44 45 46 47 49 49 59 51 52 53

С 1 С 1

прогрши РАСЧЕТА С!Щ1тшт таиш)го шш систаш -ш и пошов шшосшш

SLOUTIKE SYSrr »(HI.NL.HK.P.F,TC,HIJ.HIL,HIK.IJ,a,IK. »SrJ,srL,SIK.NG,ML,GC,H.A,T)

ПАРАМЕГГРЫ: HI - ЧИСЛО ТЕЛ

NL - ЧИСЛО ОБаЮВ С ТЕПЛОНОСИТЕЛЯМИ

НК - ЧИСЛО ОВаЮВ с ЗАДАННШИ ТЕШЕРАТУРАИИ (СРЕД)

рснг) - нсядности источников тЕплота в шах

FCNL) - КОЭФМШПи НЕРАШОМЕРНОСТИ ДЛЯ ТЕПЛОНООГГЕЛЕЙ

ТССШ) - TOfilEPATmj СРЕД

ИЫ - число СВЯЗЕЙ МЕКДУ ТЕЛАМИ

NIL - число СЖЯЗЕЯ МЕКДУ ТЕЛАМИ И ТЕПЛ0Н0СИШЯМ4

НШ - число СЖЯЗЕЙ МЕКДУ ТЕЛАМИ И СРЕДАМИ

IJ(2»NIJ) - НОМЕРА СВЯЗАННЫХ ТЕЛ I К J

ILC2»HIL) - НОМЕРА СВЯЗАННЫХ ТЕЛА I Н ОЕЕНА L

1КС2»Н1К) - HfflffiPA СВЯЗАННЫХ ТЕЛА I И СРВДЫ К

siJCNij) - проводимоста "тало - тело"

SILCHIL) - ПРОВОДИМОСТИ "ТЕЛО - ОЕЕМ С ТЕПЛОНОСИТЕЛЕМ"

SIKCHIK) - ПРОВОДИМОСТИ "тало - среда"

NG - число ПОТОКОВ МЕКДУ ОБЕМАМИ J ТЕПЛОНОСИТЕЛЕМ MLC2»NG) - НОМЕРА Н И L СВЯЗАННЫХ ОБЕМОВ С ТЕПЛОНОСИТЕЛЯМИ

СИЗ U ПОТОК ВЫТЕКАЕТ, В L ВТЕКАЕТ) GCCHG) - ПРОИЗЫдаИЯ МАССОШХ РАСХОДОВ НА УДЕЛЬНЫЕ ТЕПЛОИИОСТИ N " HI+2»NL - ЧИСЛО ИСКОМЫХ ТИйЕРАТУР

АСН.Н) - РАБОЧИЙ МАССИВ ДЛЯ ХРАНЕНИЯ МАТРИШ СИСТЕМЫ УРАВНЕНИЙ

Т(Н) - ШХОДНОЙ ПАРАМЕТР: МАССИВ, СОДЕРВАЩИЯ N1 ТИйЕРАТУР TU И 2*NL ТИЙЕРАТУР ПОТОКОВ НА ШХОДЕ А НА ВХОДЕ В ОВЕМИ С ТШОНОаГГЕЛЯМИ. внутри ПОДПРОГРАММЫ ЭТОТ МАССИВ ИСПОЛЬЗУЕТСЯ для ХРАНЕНИЯ ПРАИК ЧАСТЕЙ СИСТЕМЫ УРАШНИЙ.

DIMENSION PCl).FCl),TCCl).IJCl).IL(l),rK(l), •SIJ(1),SILC1),SIKC1),MLC1),GCC1).ACN,N).TC1) . ФОРМИРОВАНИЕ MATPffiffi ACH,N) И ЫЖТОРА ПРАВЫХ ЧАСТЕЙ T{N) Л. ОБНУЛШЕ МАТРИЦЫ И ЫЖТОРА

DO 1 1=1,N

Т(1)=в.

DO 1 J-1,N

1 ACI,J)=e.

.2. ЦИКЛ ПО ТЕЛАМ. ЗАПИСЬ МСЩНОСТЕЯ В ПРАВЫЕ ЧАС1М. DO 2 Ы.Ш

2 ТСГ)-РСГ)

3. ЦИКЛ ПО связям "ТЕЛО - таж"

IF{SIJ.EQ.e)GCrrO 4 DO 3 M=1.NIJ

ОПРЕДЕЛЕНИЕ НОМЕРОВ СВЯЗАННЫХ ТЕЛ I.J И ПРОВОфМОСП! 9 MJC2*M-1) J-IJC2«H)

S=SIJCH)

Рис. 1.4

54 55 56 57 56 59 60 61 62 63 64

70 7i 72 73 74 75 76 77 78 79 8в 61 82 83 84 85 86 67

91 82 69 94 95 69 97 98 98

185 106

ЗАПИСЬ КОЭФФИЦИЕНТОВ УРАВНШЙ ДЛЯ TU A(I,r)-A(I.I)+S A(I,J)-S A(J.JM(J.J)+S 3 A{J.I)=-S С 1 4 ЦИКЛ ПО СЖЯЗЯМ "тало - СРЕДА" 4lFCNlK.Eq.e)G0T0 6 D0 5H-1,NIK

С ОПРЕДЕЛЕНИЕ НОМЕРОВ СВЯЗАННЫХ WLK I И СРВДЫ К. ПРОВОданОЛИ 3 1=1К{2»М-1) К=1К{2«Н) S=SIKCH)

С зАПиа КОЭФФИЦИЕНТОВ УРАВНЕНИЙ для тал

ACI,I)=ACI.I)+S

5 TCI)=TCI)+S»TCCK-NL)

С 1.5. ШКЛ ПО СВЯЗЯМ "таЛЭ - ТЕПЛОНОСИТЕЛЬ"

6 IF(NIL.EQ.e)GOTO 6 DO 7 M=1,NIL

С ОПРЕДЕЛЕНИЕ НОМЕРОВ СВЯЗАННЫХ ТЕЛА I И ОВЕНА С С ТШОЮСИШЕН L. ПРОВОДИМОСТИ S

I=IL(2«M-1)

UIL(2«H>

S=SIL(H)

С LL - НОМЕР УРАВНЕНИЯ ДЛЯ ТИйЕРАТУта НА ШХОДЕ С LUl - НОМЕР УРАВНЕНИЯ ДЛЯ ТЕШРРАТУИ) НА ВХОДЕ LL-NI+2»L-1

С зАПиа коэФФищтагов уравнений для тел

ACI.I)=ACI.I)+S ACI,LL)F(L)«S ACI,LL+1)C1.-F(L))> С ЗАПИСЬ КОЭФФИЦИЕНТОВ УРАВНЕНИЙ ДЛЯ ТЕПЛОНОСШЛЕЯ. ACLL.I)=-S

ACLL.LL)=ACLL.LL)+FCL)«S

7 ACLL,LUl)=A(LL,LL+1)+(1 .-F(L) )»S

С 1.6, цикл ПО ПОТОКАМ МЕЖДУ ОБЕМАМИ С ТЕПЛОНОСИТЕЛЯМИ

8 IF{NG.EQ.e)GOTO И DO 10 K=1.NG

С Я1РЩЛЕНИЕ НОМЕРОВ СВЯЗАННЫХ ОБЕМОВ СИЗ OBIMA М ЙПШЕТ.. С В ОВЕМ L ВТЕКАЕТ) И ПРОИЗВЕДЕНИЯ C*G

H-MLC2»K-1)

L-MLC2*i()

CG=GCCK)

С ЗАПИСЬ КОЭФФИЦИЕНТОВ УРАШНИЙ ДЛЯ таялоноситЕЛЕй

LL=NI+2»L-1

A(LL,LL)=ACLL.LL)+CG

ACLL,LUl )=A(LL.LL+1)-CG

ACLUl, LUl )=A(LU1 ,LU1 )+CG С ЕСЛИ М ЕОДЫЕ NL, ТО М - ЮМЕР СРЕДЫ С ЗАДАЖОЗ ТЕМПЕРАТУРОЙ.

IFCH.GT.NDGOTO 9

M№=NI+2»M-1

ACLUl,MM)=-CG

GOTO 10

С ЗАПИСЬ ПРАВЫХ ЧАСТЕЙ УРАВНЕНИЙ ДЛЯ ТЕПЛОНОСИтаЛЕЙ Рис. 1.4. Продолжение

т? 9 TCLUl)-T(LUl)+CG»TCCiHfl,)

108 10 COirriNUE 1в9 С 2. РЕЗВЕШЕ СИСТЕМЫ УРАВНЕНИЙ

110 И еда GELGCT.A,H,1.1.E-7.IER)

111 IFCrER.Iffi,0)PRIHT Ig.IER

112 12 KHaiATC5X.IER=,I2)

113 С 3. ПЕЧАТЬ РЕЗУЛЬТАТШ

114 PRINT 13

115 13 FOMOiTС/5Х.ТЕМПЕРАТУРЫ Т£Л) 118 ГЮ 14 1=1,HI

117 14 PRINT 15.1,TCI)

118 15 F0MIATC2X,I3.G12.3) 118 PRINT 18

1Йв 16 Р0М1АТС/5Х,ТИЙ1ЕРАТУИ) ТЕПЛОНОСИИЖЙ/

121 «IX,НОМЕР и.ВХОДА и.ВЫХОДА)

122 ГЮ 17 L-1,NL

123 LI>2»UHI

124 17 fHINT 16.L,TCLL),TCLH)

125 18 F0R1IATCI5.2G12.3)

126 шт

127 m

Рис. 1.4 Продолжение

ТОКОВ теплоносителей из одного объема в другой. Массив ML содержит номера связанных объемов: из объема номер ML (2 * К - 1) теплоноситель втекает в объем с номером ML (2 * К), К 1. NG. В массиве GC записаны произведения удельных теплоемкостей на массовые расходы Cjjjmi- Поток может вытекать либо нз среды с за-данЕюй температурой O/iGj в уравнении (L28)l, либо нз объема с неизвестной температурой {Gz)- Чтобы в программе различить эти две разные ситуации, использован следующий прием: номера сред с заданной температурой в массивах ML и 1К начинаются нес единицы, а с NL --f 1, где NL - число объемов с теплоносителями.

Для формирования матрицы А и столбца свободных членов Т организуются пять циклов: по телам; по связям между телами {оУ/)\ по связям между телами и теплоносителями (o]f); по связям .между телами и средами (аЦ) и по потокам теплоносителя между объемами (Gmb Oki)- В каждом цикле последовательно выбирается связь, определяются номера взаимодействующих тел н объемов, находятся номера строк и столбцов, определяющие коэффициенты матрицы и столбца свободных членов, в которые вносит «вклад» данная тепловая связь, и. наконец, к ранее вычисленному значению коэффициента матрицы прибавляется «вклад» от рассматриваемого взаимодействия.

Например, в цнкле по связям между телами и объемами с теплоносителями на основе массивов IL и SIL проводится учет вкладов

матрицу А проводимостей ajf, которые входят в ее первые гтоок [уравнения (1.26)] и в строки с номерами N.. + I, /V + 3, ... *дг }-2 * Лш [уравнения (1.27)1. Номер тела 1 совпадает с но-епом температуры Ti в столбце неизвестных, а по номеру объема L определяются номера температур теплоносителя: (N1 + 2 * L - 1) - для температуры УГ"" и (N1 + 2 * L) - для температуры UT-Далее в коэффициенты матрицы А с указанными номерами добавляются слагаемые ajf, о]Тfi, o1T{l-fi) согласно уравнениям (1.26), (L27).

После формирования матрицы А и столбца Т выполняется обращение к стандартной подпрограмме GELG, выходным параметром которой является массив Т. В этом массиве первые Л/ элементов соответствуют температурам тел Ti, а затем парами расположены температуры теплоносителей У""", У". Эти температуры выводятся на печать. Если в подпрограмме GELG параметр ошибки 1ER принимает значение, отличное от нуля, то печатается соответствующее сообщение об ошибке.

§ 1.5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ КОШИ ДЛЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

Расчет нестационарного теплового режима по моделям с сосредоточенными параметрами сводится к решению систем уравнений теплового баланса вида (1.2), (1.3) с начальными условиями (1-6), т. е. к решению задачи Коши для систем обыкновенных дифференциальных уравнений первого порядка. В случае лннейиых уравнений решение удается представить в аналитическом виде при числе уравнений Л/ < 4. Для нелинейных задач и в случае Л/ :> 4 точное решение в аналитическом виде получить ие удается, за исключением некоторых частных случаев. Поэтому при расчетах нестационарных тепловых режимов систем тел широко применяют численные методы, которые мы сначала рассмотрим применительно к одному уравнению вида

dT dT

Начальным условием

= f(T, Г), 0<т <тг

(1.29)

(L30)

Cl in - методы, изложенные на примере задачи (1.29), jO), будут обобщены для системы Л/ уравнений.

исновиые понятия теории численных методов решения дифферен-Ных уравнений будут достаточно подробно рассмотрены в гла-ча примере дифференциального уравнения теплопроводности. Ичас лишь кратко сформулируем ряд понятий, которые поиадо-

бятся для изложения практических вопросов численного решения обыкновенных дифференциальных уравнений.

Прн численном решении задачи непрерывная область изменения независимой переменной [О, тах! заменяется множеством значений которые будем называть узлами сетки. В случае равномер. ной сетки ij /Дт. / = I. /; ДтТщах-шаг по времени. Вместо задачи отыскания непрерывной функции Т (т) ставится задача определения дискретного множества значений функции в узлах сетки: TJ = T{tj). Величина называется сеточной функцией точного решения. Как мы увидим дальше, точные значения Ti найти не удается, а вместо них в результате численного решения задачи получаются приближенные значения искомой функции в узлах сетки, которые будем обозначать и и называть сеточной функцией разностного решения или просто разностным (численным) решением. Погрешность численного решения определим как разность сеточных функций точного и разностного решений: е/ = Т - и.

Для нахождения разностного решения вместо задачи (1.29), (1.30) рассматривают задачу решения системы алгебраических уравнений относительно искомых значений Такую систему алгебраических уравнений называют разностной схемой. При измельчении сетки (при Дт 0) погрешность численного решения должна в случае удачной схемы стремиться к нулю, т. е.

lim е/=0,

Д;-D

Это требование к разностной схеме называют условием сходимости. Для сходимости разностной схемы необходимо и достаточно выполнения двух других условий - аппроксимации и устойчивости, которые будут пояснены ниже на примере схем Эйлера.

Если при достаточно малых Дт начинает выполняться неравенство

е/]=.П -ul </1Дт, /l-=const,

(I.3I)

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

Теперь перейдем к рассмотрению наиболее распространенных методов численного решения задачи (1.29), (1.30). К ним относятся: методы, основанные на разложении функции Т (т) в ряд Тейлора (наиболее распространенная схема этого вида - схема Эйлера); методы Рунге-Кутта; линейные многошаговые методы.

Схемы Эйлера. Используя разложение решения Т (т) в ряд Тейлора в точке Tj, можно записать следуюш,ее выражение для значения решения в точке Tj- - ту + Дт:

Г/+ iTi + T (Tf) Дт -f Г (Tj) ДтV 2 + ....

(1.32)

Пренебрегая в (1.32) членами порядка О (Дт) и выше и учитывая, согласно исходному уравнению (1.29) производная Т (т;) равна / ("j нУ следуюш,ую схему:

1 u/ + At/(i;j, «о, (1.33)

называемую схемой Эйлера, Обратим внимание, что вместо мы записали «Л так как отбросив члены порядка О (Дт) в (1.32), мы уже ие получим точное решение Л(т).

Из (1-33) вытекает, что значение и" вычисляется на основе значения и Д. предыдуш,его момента времени по явной формуле. Подобные схемы называются явными. Зная из начального условия (1.30) значение и" в момент времени т - О, легко по формуле (1.33)

последовательно вычислить значения «\ и-.....и- во все последуюш,ие

моменты времени.

Кроме явных суш,ествуют неявные схемы, в которых значение искомой функции на новом временном слое находится в результате решения уравнения. включаюш,его это значение и значения для предыдуш,их моментов времени. Неявную схему Эйлера можно получить, если использовать разложение в ряд Тейлора в точке т; TJ1 - Т (T,-i.i) Дт. Тогда придем к схеме

(1.34)

при реализации которой на каждом шаге необходимо решать в об-ш,ем случае нелинейное уравнение относительно и-.

Проведем сопоставление явной и неявной схем Эйлера. С точки зрения объема вычислений для одного шага явная схема имеет пре-имуш,ество. Только в случае, когда функция / (т, Т) линейна относительно Т, т. е. f {%) а (т)Т -\- b (т), вычисления по неявной схеме не сложней, чем по явной, поскольку тогда уравнение (1.34) разрешается относительно

! [«/ + 6(T,.+i) Дт]/[1 -Q(T,-+i) Дт]. (1.35)

Если же / (т, Т) - нелинейная функция Т, то приходится решать на каждом шаге уравнение (1.34) каким-либо итерационным методом.

Рассмотрим теперь вопрос о погрешностях численных решений, получаемых по явной и неявной схемам Эйлера. Для этого введем понятия аппроксимации и устойчивости.

Вернемся к разностному уравнению (1.33), переписав его в виде

= /(ту, и).

[1.36)

Поскольку это уравнение получается нз соотношения (1.32) для путем отбрасывания малых членов, точное сеточное решение Р удовлетворяет уравнению (1.36). Величину я!)/, характеризующую п неудовлетворения» функции Т разностному аналогу (1.36) Дифференциального уравнения (1.29), очевидно, можно определить