Поясним назначение параметров подпрограммы. Сначала отметим что подпрограмма R KGS обращается к двум подпрограммам: Fcf для вычисления правых частей системы уравнений, записанной в виде (1.61), и ОиТР для вывода результатов расчета. Эти подпрограммы составляются пользователем для конкретной задачи. Их имена являются формальными параметрами RKGS, и они должны быть описаны в головной программе с помощью оператора EXTERNAL.

В массив PRMT длиной пять или более элементов записываются следующие входные параметры: PRMT (1) - начальное значение аргумента тц, PRMT (2) - конечное значение аргумента т., PRMT (3) - начальное значение шага интегрирования Дт, PR.MT (4) - допустимое значение локальной погрешности е / .

Оценка локальной погрешности к решения системы уравнений (и}, на шаге / находится иа основе оценок погрешностей е1 каждой из искомых функций следующим образом:

елок = 2 i =--1

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

Т - массив для хранения N текущих значений искомой сеточной функции, в который должны быть записаны начальные значения {i= 1, .... iV).

DERY - массив для хранения N текущих значений производных /, (Tj, u(, uQ. При обращении к R KGS в него должны быть записаны весовые коэффициенты ф,- для расчета локальной погрешности по N значениям сеточной функции.

IHLF - выходной параметр, который показывает, сколько раз первоначальное значение шага PRMT (3) делилось на два в процессе автоматического выбора шага с целью достижения заданной локальной погрешности. При IHLF-11 происходит передача управления в головную программу.

AUX -рабочий массив длиной 8*N.

Обращения к подпрограммам, которые должен составить пользователь, имеют следующий вид. Для подпрограммы вычисления производных: CALL ЕСТ (TIME, Т, DERY), где TIME - текущее значение аргумента т; Т - массив текущих значений функции uf {{ -= = 1, jV); DERY - выходной массив значений производных

Л {-h "У-

Для подпрограммы вывода результатов:

CALL OUTP(TlME. Т, DERY. IHLF, N, PRMT),

смысл всех параметров был пояснен выше. Ззметим, что в процессе автоматического уменьшения или увели-ения шага Дт первоначальное значение РРЛ\Т (3) ие может быть превышено.

6 ПРОГРАММНАЯ РЕАПИЗАЦИЯ РАСЧЕТА НЕСТАЦИОНАРНЫХ СРЕДНИХ ТЕМПЕРАТУР

Рассмотрим более простую по сравнению с описанной в § 1.1 систему, включающую только твердые тела и среды с постоянными температурами. В этом случае нестационарный тепловой режим описывается системой уравнений

/=1 k=]

1=1, N,

С начальными условиями

(1.63)

(1.64)

Решение по схеме Эйлера. Сначала остановимся на программе решения задачи (1.63). (1.64) по неявной схеме Эйлера, которая имеет вид

/4-1 / -

2 а-К+-е,). N„ / = 0, 1, (1.65)

k= 1

здесь - сеточная функция, соответствующая значению температуры Tt в момент времени т, = /Дт; Дт - величина шага по времени. Отличие неявных схем для одного линейного уравнения и для системы уравнений состоит в том. что разностная схема (1.35) разрешалась в явном виде относительно м+, а в данном случае мы имеем систему N- линейных алгебраических уравнений для определения значений сеточной фуикции 1)+ мЧ. м+}.

Запишем эту систему в матричной форме

AU+ =В+, (1.66)

где А - матрица размером Лт X Лт. коэффициенты которой вычисляют по формулам

а В/+1 - вектор-столбец с элементами.

(1.68)

Таким образом расчет по неявной схеме Эйлера сводится к решению на каждом шаге по времени системы линейных уравнении (1.66), которое может быть выполнено с помощью какой-либо стандартной подпрограммы. В рассматриваемой задаче матрица А является симметричной, так как согласно (1.67) a-i

12 13 14 !5 16 17 16 19 2в 21 22 23

ГСЖШНАЯ ПРОГРАША РШНИЯ СИСТЕМЫ НЕСТАЦИОНАРНЫХ УРАВНЕНИЙ

ТОНОВОГО БАЛАНСА ПО НЕЯВНОЯ СХИЕ ЭЙЛЕРА

DiarXSIOH Р(28).С(г8).Т(2в).ТС(5).Ы(гвв). «IK(5e),Sy(5e),SIK(25).TV(28).A(4ee).AUX(2«)

ела WOD(Nl,HK.P.C.T.TC.NIJ,NIK.IJ,IH, »SIJ,SIH.TAU.11£AX,TV)

ЦИКЛ ПО ЗРОЕНИ

1 TlUE-TItfE+TAU

ела MATOI(NI,HK,P.C,T.TC.NIJ,NIK,IJ,IK. »SU.SIH.TAU.A.T) CALL GELS{r.A.NI.l,i.!>-7.1B?.AU)() IF(IER.fffi.0)PRIW 2.Ш

2 FORMATC ОШИБКА Ш=,12) IF(TI11E.LT.TV(L))G0T0 3 L-UI

reiNT 4.T1UE.(T(I).I-1,NI) 4 FCMAT{/ ВРЕЛЯ-.0!в.З/

» гшшРАТзте ш/(вси.з))

3 IFdlHE.LT.TVAXJGOTO 1

ггор ш

Рис. 1.6

поэтому используется подпрограмма GELS (см. § 1.3).

Приводимый ниже пример программной реализации включает головную программу и подпрограммы VVOD, .MATRN, GELS (рис. 1,6-1.8). В головной Программе устанавливается максимальная размерность для всех используемых .массивов, выполняется обращение к подпрограмме ввода данных VVOD и организуется цикл по времени. В этом цикле по времени на каждом шаге проводится форми-

И 12 13 14 15 16 17 IS

2в 21 22 23 24 25 26 27 28 29

4в 4[

42 43 44 -45 46 47 48 48 5в Ъ1 52 53 54

С ПОДПРОГРАММА ВВОДА ДАННЫХ ДЛЯ РЕШЕНИЯ СИСТЕМЫ С НЕСТАЦИОНАРНЫХ УРАВНЕНИЙ ТЕПЛОВОГО БАЛАНСА PUBROUTlJffi WOD

i»{?JI,NK.P.C,T.TC.NIJ,NlH.IJ,IH, »SIJ,SlK.TAU.T!ilAX,TV)

DIMENSION Р(П,С<1),Т{1).ТС{1),Ы{(). «lK(!).SIJ(i).SIK(i),TV(I)

READ J.NI.NK

1 F0RMAT(iei5) PRINT S.NI.NK

2 FORMATC ЧИСЛО ТЕЛЛЗ. ЧИСЛО СРЕД.13) READ 3.(Р(1),1.1.МП

3 F0RMAT(5Fie.3) PRINT 4.(P(I),I=!,NI)

4 FORMATC M0H0CTH7(5Fie.3)) READ 3,(C(I).I=1,NI)

PRINT 5.(C(n,I=l.NI)

5 FORMATC ТЕПЛ0Ш0СП)7!&Г1г1.3)) READ 3.(1(1).I-l.NI)

PRINT 6,(T(I).I=1.NI)

6 FORKATC НАЧАЛЬНЫЕ ТШ1ЕРАТУРЫ7(5Е!0.3)) READ 3,(TC(I).I=!.NH)

PRINT 7.(TC{I),I=!.NK)

7 FORMATC ТЕМПЕРА17РУ CPEa7(5Fie.3)) READ l.NIJ.NIK

PRINT S.NIJ.NIK 6 FORMATC NIJ=M3. Nll{=,13) N=2*NIJ

READ J.IIJ(]),I=1.N) READ 3,(SIJ{I).I=1,H1J) PRINT 8

9 FORKATf ШДОВЫЕ СВЯЗИ ШШ ТШШ/

. # НОМЕР I НОМЕР J проводамостьi

DO 10 1=1.NIJ le PRINT (I.IJ(2#I-n.IJ(2«l),SIJ(I) ti F0RMAT(14.11I,F15.3)

N=a»NIH

READ l.dHdJ.M.NJ READ 3,(SIK(I).I=1.NIH) PRI.VT J2

!2 FORMATC ТЕПЛОВЫЕ СВЯЗИ МЕЯДУ ТЕЛАМИ И СРЕДАМИ7 » НОМЕР 1 HOMLP К ПРОВОДИМОСТЬ) DO (3 I=i,NIK

13 RUNT U.IH(2#I-1).IK(2»I).SIK(I> READ S.TAU.TVAX

PRINT J4.rAU.Tm

14 FORMATC ШАГ ПО ВРЕМЕНИ.Е!в.З/ * МАКС. HEJW.FIO.S)

READ Ijv

READ 3i(TV(I),M.IV) PRINT i5.(TV(I).I-I.IV)

15 FORMATC ВРЕМЕНА ВЫВОДА НА ПЕЧАТЪ/(5Р1в.З)) . RETURN

Рис. 1.7

1 с П0Д1РОГРАМНА ФОРМИРОВАНИЯ МАТРИЦЫ ДЛЯ РЕШЕНИЯ ОКЛВЫ

2 С НЕСТАШЮНЛРНЫХ УРАВНЕНИЙ ТЕПЛОВОГО БАЛАНСА

3 С ПО НЕЯВНМ) СХЕМЕ ЭЙЛЕРА

4 SUBROUTINE MATRN

5 •(NI.KK.P,C,T,TC,NU,NIH.IJ.IH.

6 «SIJ.SIH.TAU.А,В)

7 С ВХОДНЫЕ ПАРАМЕТРУ:

8 С N1 - ЧИСЛО ТЕЛ

9 С NK - ЧИСЛО СРЕД С ЗАДАННОЙ ТЕМПЕРАТУРОЙ

1в С P(NI) - мощности ИСТОЧНИКОВ ТЕПЛОТЫ В ТЕЛАХ

И С C(NI) - ПОЛНЫЕ ТЕПЛОЕМКОСТИ ТЕЛ

12 С T(NI) - та1ПЕРАТУРЫ ТЕЛ НА ПРЕДЫДЛПЕМ ШАГЕ ПО ВРЕМЕНИ

13 С ТС(т) - ТШ1ЕРА1УРЫ СРЕД

14 С NIJ - ЧИСЛО СВЯЗЕЙ ИЕВДУ ТЕЛАМИ

15 С Ш - ЧИСЛО СВЯЗЕЙ ИЕЙДУ ТЕЛАМИ И СРЕДАМИ

16 С IJ(2»HIJ) - НОМЕРА I И J СМЗАНШа ТЕЛ

17 С I1{{2»N1H) - НСЖРА СВЯЗАННЫХ ТЕЛА I И СРЕЦЫ Н

18 С sij(niJ) - тволтхгги межпу телами

19 с SIHCNIH) - ПРОВОДИМОСТИ МЕЙДУ ТЕЛАМИ И СРЕДАМИ

20 С TAU - ШАГ ПО ВРЕМЕНИ

21 С ШХОДНЫЕ ПАРАМЕТРЫ:

22 С А - МАССИВ ДЛЯ ХРАНЕНИЯ ВЕРХНЕЙ ТРЕУГОЛЬНА! ЧАСТИ СИММЕТРИЧНОЙ

23 С МАТРИЦЫ СИСТЕМЫ УРАВНЕНИЙ НЕЯВНОЙ СХЕМЫ ЭЙЛЕРА. ДЛИНА МАССИВА

24 С РАВНА NU(NI+I)/2

25 С B(NI) - ВЕКТОР-СТОЛБЕЦ СВОВОДНЫХ ЧШОВ

26 С

27 DIMENSION Р(1),С(1).Т(i).ТС( 1).Ud). 1И( 1).SIJ{1),SIH( 1),

28 »А(1),В(1)

29 С 1. ОБНУЛЕНИЕ МАТРИЦЫ Зв NA=I*(NI+l)/2

31 DO I r=I.NA

32 1 A(I)=e.

33 С 2. ЦИКЛ ПО ТЕЛАМ. ЗАПИСЬ Б ДИАГОНАЛЬНЫЕ ЭЛЕМЕНШ МАТИВЩ А

34 С И Б ЗЕКГОР В ПАРАМЕТРОВ ТЕЛ

35 DO 2 1=1.N1

36 K=I«(I+l)/2

37 A(K)=C(I)/TAU

38 2 B(I)=P(I)+C(I)/TAU.T(I)

38 С 3. ЦИКЛ ПО СВЯЗЯМ "ТЕЛО - ТЕЛО". ЗАПИСЬ В МАТРИЦУ А"

40 С ТЕПЛОВЫХ ПРОВОДИМОСТЕЙ

41 IF(NIJ.EQ.e)GOTO 4

42 ЮЗ H=I,fJIJ

43 l=IJ(2»M-l)

44 J»IJ(2»M)

45 S=SIJ(M)

48 K=I#n+l)/2

47 A(K)=A(K)fS

48 K«J»(J+l)/2 43 A(K)=.A(K)+S

50 IF(I.LE.J)K=J«(J-l)/e+I

51 IF(I.GT.J)K-b(I-i)A+J

52 3 А(К)-А(И)-3

53 С 4. mm ПО СВЯЗЯМ "ТЕЛО - СРЕДА". ЗАПИСЬ ПРСвОДвЮСГПЯ

54 55

57 58 58 60 81 62 63 64

В ДИАГОНАЛЫШЕ ЭШЕНШ МАТРИВД А И ТШЮВНХ ПОТОКС»

В ВШОР В

4 1Е(ШК.ВД.0)С»ТО 6 I>0 5M=1,NIK 1=Ш2»»*-1) К=1К(2»М) J=I»(Ul)/2 A(J)=A(J)+S]K{W

5 B(I)=B(I)+SIKtM)*TC(K)

6 ШШ

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

оованне матрицы А и столбца В путем обращения к подпрограмме MATRN, решение системы (1.66) с помощью подпрограммы GELS, а в заданные моменты времени выполняется вывод температур {"jlfjj на печать. Начальное распределение и" = Т;о задается в подпрограмме ввода путем заполнения массива Т значениями начальных температур тел. Заметим, что при обращении к подпрограмме MATRN массивы температур Т и свободных членов В совмеи1.ены.

Все исходные данные описаны в комментариях к тексту подпро-гpaммыMATRN. Способ описания взаимодействий между телами и межд> телами и средами {aV- и аТ) полностью идентичен рассмотренному в§1.4. Формирование массивов А и В, соотвегствуюн1Н\ матрице и столбцу свободных членов, проводится после их преднарительного обнуления также путем последовательного суммирования «вкладов» от отдельных тел и тепловых связей согласно формулам (1.67), (1.68). Отличие от рассмотренной ранее подпрограммы SYSTT состоит в том, что симметричной матрице А соответствует одномерный массив, содержащий лишь верхнюю треугольную ее часть, .(аиисанную по столб1.ам, как это требуется для подпрограммы GELS. Поэтому после выяснения индексов i. / коэффициента матрицы Оц выполняется расчет индекса п для одномерного массива по формуле п - / (/ -- I) 2 -Ь i. Замети.м, что если положить С; = О и выполнить один шаг по времени, то сразу патучим решение соответствуюп],ей стапио-нарной задачи.

Отметим одно обстоятельство, обусловленное ли[ейностью рас-"матриваемой задачи (L63), (1.64). Поскольку компоненты матрицы не изменяются во времени при const, Ci const, целесообразно один раз перед началом цикла по времени вычислить обратную Матрицу А-1 а затем дальнейшие расчеты в цикле свести к формированию столбца В+ и определению температур U" путем умноже-"я обратной матрицы на вектор-столбец; 1)+ -•= А-В+\ Обратная матрица может быть найдена, например, с помощью подпрограммы MINV (см. §-1.3). Таким путем можно достичь значительной ономии машинного времени по сравнению с формированием матри-ьt и решением системы на каждом шаге по времени.

Рис. 1.8