1в 11

12 13 14 15 16 17 16 18 2в 21

С ГОЛОВНАЯ ПРОГРАММА ДЛЯ РИЗШЯ ОЛЧ УРА8НШЙ С ТЕПЛОВОГО БАЛАНСА ПО СХЕМЕ РУНГЕ-ЙТТА EXTERNAL FCT.Oin?

DIMEJJSION T(2a),DERY(2e),PRW(5),iJX(i6e) COMMON /CFCr/NI.HK,P(2a),C(2e),-Ki5),[JIJ NIK. »IJ(1W),IH(50),SIJ(M),SIK(25) COMMON /СОиТР/ L.TV(20)

CALL WODCNl.NK.P.C.T.TCfJIJ.NIHJJ IK •SIJ.SlK.TAU.mX.TV) L=l

PRJfi(l)-0. PRJfr(2)-TMAX PRW(3)-TAU PR«T(4M.ei Ю 1 W,(JI

1 DERy(I)-l./Wl

CALL RKGS(PRl(T,T.DHllf.(JI,lHLF,FCT,OUlP,AUX) IFdHLF.CE.ielPRHfr 2,IHLF

2 FORMATC IHLF-, 13, ИЗШЗОПЪШАГШ П0ГРИН0С1Ъ> STOP

Рис. 1.9

Однако если или Ci зависят от TCMnepaiyp, то на каждом шаге / приходится вычислять их новые значения на основе температур и, следовательно, снова формировать матрицу А н решать систему уравнений. Поэтому выше рассмотрена структура программы, которая пригодна и для таких нелинейных задач, решаемых путем пересчета матрицы на основе температур предыдущего шага по времени.

Решение по схеме Рунге-Кутта. Перейдем к программе решения задачи (1.63), (1.64) по схеме Рунге-Кутта четвертого порядка. Она строится на основе описанной в §1.5 стандартной подпрограммы R KGS, в которую уже «заложен» цикл по времени. Поэтому в головной программе (рис. 1.9) реализуется лишь задание размерности массивов, ввод исходных данных и обращение к RKGS. Форма представления исходных данных совпадает с использованной в предыдущей программе для схемы Эйлера, а для ввода применяется та же подпрограмма WOD. Входными данными для подпрограммы RKGS являются; начальные значения температур (массив Т). весовые коэффициенты <pj, полагаемые равными для все неизвестнь[Х (массив DERY), число неизвестных (Nl). а также массив PRMT, содержащи четыре значения: начальное и конечное значения времени (О и ТМАХ). начальный шаг (TAU), допустимую локальную погрешность (0.01)-Для испатьзовання подпрограммы R KGS требуются две подпрограммы: FCT - расчета правых частей в системе уравнений вида (1.61) н ОиТР - вывода на печать результатов расчета (рис. 1.10). Для передачи в эти подпрограммы из головной прогрзм-

.я исходных данных, которые не являются параметрами под-Р** у R KGS, нспользукугся две общие области COMMON с име-

"PuCFCT н COUPT.

R подпрограмме FCT вычисляются значения производных dTi/dx ментов массива DERY) согласно формулам (1.63) на основе мас-в мощностей, температур тел и сред, а также массивов U, SIJ. kSI к. указывающих номера взаимодействующих элементов и теп-ггпвые проводимости.

Б подпрограмме OUTP проводится сравнение текущего момента пемени TIME с временами вывода результатов, которые заданы массивом TV. В случае TIME > TV (L) реализуется вывод на печать, да-

2 3 4 5 6 7 6 9 1в И 12 13 14 15 16 17 16 18 29 21 22 23 24 25 26 27 26 28 30 31 23 33 34 35 36 37 36

ПОДПРОГРАММА ВЫЧИСЛЕНИЯ ПРОИЗВОДНЫХ ДЛЯ РЕШЕНИЯ ПО СХЕМЕ РУНГБ-КША

SUBROUTINE FCT(TIHE,T,Dm)

DIHEKSION Т(1),СШ(1)

COMMON /CFCT/NI,NK,P(2e},C(2«).TC(5),NIJ,NIK. •IJ{1W),IK(50),SIJ(50).SIK(25) ЗАПИСЬ УОЦЮОЕЛ DO 1 Ы.Я1

1 DERY(n-P(I)

ЗАПИСЬ ШЛОВЫХ ПОТОКОВ Шт ТЕЛАШ К) 2 M-l.NIJ I-IJ(2»M-1) J-IJ(2»H)

Dm(I)-DERy(I)-SIJ(H)»(T(I)-T(J))

2 №(J)-DERY(J)-SIJ(M)#(T(J)-T(I)) ЗАПИСЬ ТЕПЛОВЫХ ПОТЖОВ В СРЕДЫ

DO 3 M-l.NIK

I-IK{2»M-1)

K-IK(2»M)

3 DIRY(I)-DERY(I)-SIK(M)*(T(I)-TC(K)) ДЕЛЕНИЕ НА ТШОИКОСТЙ

DO 4 I-l.NI

4 DERY(:)-=DERY(I)/C(I) REmiRH

ПОДПРОГРАММА ПЕЧАТИ РЕЗУЛЬТАТОВ ПРИ РЕШШ ПО СХЕМЕ РУНГЕ-КУТГА SUBROUTINE OUTP(TIME,T,DERY,IHLF,NI.PWfr) DIMHJSION T(1),DERY(1),PR«T(1) COMMON /СОиТР/ L,TV(2e) IF(TIME.LT.TV(L))GOTO 2 UL+1

Ffl-PRUr(3)/2»»IHLF

PRINT l.TIME,re,(T(I).I"l.NI)

1 FORMATC/ ШМЯ-,С1в.З. ШАГ-.Gift.3/ . TOmXPAIVPU ТЕЛДбСИ.З))

2 REmiRH

Рис. 1.10

лее сравнение проводится со следующим элементом массива TV (L + 1) и т. д. На печать выводится значение текущего моуец та времени, текущего шага интегрирования и jV. температур -j. Значение шага (PR) определяется на основе начального значение (PRjMT (3) TAU) и целого параметра IHLF, который указыва сколько раз начальное значение шага делилось пополам.

Если начальный шаг в подпрограмме RKGS уменьшался вдвое 11 раз (IHLF II), то управление передается в головную програц,. му, расчеты прекращаются и печатается соответствующее coo6ai,eHHf В этом случае следует в исходных данных изменить начальное значение m;ira (TAU) или погрешность (PRMT (4)) и повторить расчет

ГЛАВА

РЕАЛИЗАЦИЯ НА ЭВМ ТОЧНЫХ АНАЛИТИЧЕСКИХ РЕШЕНИЙ

Для многих задач расчета пространственных температурных полей в телах канонической формы могут быть получены точные аналитические решения. Однако для нестационарных одномерных и любых дву- и трехмерных задач эти решения записываются в виде рядов, интегралов, часто содержат специальные функции. Во многих случаях в аналитические выражения входят параметры, являющиеся корнями трансцендентных уравнений и систем таких уравнений, которые могут быть решены лишь численно. Поэтому расчеты пространственных температурных полей на основе точных аналитических решений также требуют применения ЭВМ.

Рассматриваемые в главах 3-5 численные методы расчета позволяют решать значительно более широкие классы задач по сравнению с аналитическими методами. Однако тем не менее использование ночных аналитических решений при расчетах на ЭВМ температурных полей в ряде случаев весьма полезно, Это вызвано следующими обстоятельствами. Во-первых, эти решения используют в качестве тестовых при анализе различных численных схем. Во-вторых, применение аналитических решений часто позволяет существенно сократить затраты машинного времени и памяти, так как число пространственно-временных точек, в которых находятся значения искомой функции, опреде-тяется только объемом требуемой информации об исследуемом процессе. При испатьзовании же численных методов число узлов пространственно-временной сетки, необходимое для получения разностного решения с удовлетворительной точностью, как правило, оказывается существенно б6льи]им. Кроме того, реализация многих раз-

hix схем требует больших дополнип Д1Я хранения рабочих массивов при

ительных затрат машнннон па-решении систем разностных

мятн Д упавненни.

В данном разделе на примере одной задачи теплопроводности рас-пн-М типичные задачи численного анализа, возникающие при реа-Jiym точных аналитических решений, и методы их решения, гвопросамн построения точных решений задач теории теплопроводности и конвективного теплообмена можно познакомиться по учебным пособиям 13, 13].

§ 2.1. ПОСТАНОВКА ЗАДАЧИ

Рассмотрим задачу расчета нестационарного одномерного темпе-патурного поля в неограниченной пластине толщиной /. В пластине распределен источник теплоты, имеющий объемную плотность мощности (/ц (.V). Поверхность пластины л О теплоизолирована, а на поверхности л- - I происходит теплообмен со средой по закону Ньютона. Начальное распределение температуры равномерное, и эта температура отлична от температуры среды. При такой постановке задачи уравнение теплопроводности и краевые условия имеют вид [311

= 0,

= 0,

(2.1)

\2.2) 2.3)

здесь < t - t,

перегрев над температурой среды. Решение (.v, т) задачи (2.1)-(2.3). исходя из принципа суперпозиции, представим в виде суммы решения , (х, т) однородного уравнения при q,, о с начальным условием (2.3) и решения 2 {х, т) неоднородного уравиения (2.1) с нулевым начальным условием.

Решение Э-, (х, т), описывающее охлаждение неограниченной пластины без источников теплоты, получено в учебнике [9] и имеет вид

(-V.

0 I sin \in fin II r-

(2.4)

Где [ij собственные числа краевой задачи, являющиеся положительными корнями трансцендентного уравнения

=~- а А - критерий Био; ЦДЛ И cos (.1,.л );

\\fnf~-\rnix)dx О

- нормы

(2.5)

«.обственныX функций

(2.6)

Второе слагаемое в решении - 0 (х, т), описывающее нагрев тастииы внутренним источником (х), получим методом конечных интегральных преобразований 131. Применим к уравнению (2.1) интегральное преобразование

(2.7) \

Тогда для изображения e„ (т) получим обыкновенное дифферен-циальиое уравнение

е„(т)

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

е„ (0)=0;

здесь W„ - изображение функции q„ (х):

Решение задачи (2.8), (2.9) записывается в виде

1 - ехр

(2.8) (2.9)

(2.10) (2.11)

Переходя к оригиналу (х, т) по формуле обращения, получим 2{х, -г) = 2 п(г)!пШ\[п\\ =

п= I

COS р„

. (2.12)

Окончательное выражение для решения задачи (2.1)-(2.3) имеет вид

ехр(ряРо) +

.JEn-exp( PFo),jl

здесь Ро = ахП - число Фурье.

Таким образом, в рассмотренном примере для проведения р"* четов по точному решению (2.13) необходимо, ю-первых, нитй значений собственных чисел р из трансцендентного уравнения {N - число учитываемых членов ряда); ю-вторых, вычислить

(2 10) для изображения W, которые для некоторых (?„ ix) геГрзУрожаться помощью элементарных функций, и, наконец, ланных пространственных точек Xt, i = 1, / в заданные мо-дя g gpeHH т,-, / -= 1, J провести суммирование членов ряда, "далогичные с позиций вычислительной математики задачи воз-т для .многих точных решений задач теории теплопроводности и "ектиБНОГО теплообмена. Поэтому далее рассмотрим методы реше-"нелинейных уравнений, методы численного интегрирования, а "*же приведем некоторые реко.мендации по программной реализации точных аналитических решений.

2 2 РЕШЕНИЕ НЕЛИНЕЙНЫХ УРАВНЕНИИ

Нелинейные уравнения, содержащие тригонометрические, экспоненциальные, показательные и другие функции, называются трансцендентными. Собственные числа краевых задач определяются пут№1 решения трансцендентных уравнений

/(р)0, (2,14)

где Д1я) содержит тригонометрические илн некоторые специальные функции, например функции Бесселя, полиномы Лежандра и т. д. Для рассматриваемого примера

f(p) = p/Bi-ctgp-0, (2.15>

и уравнение (2.15) имеет бесчисленное множество корней pi, р., которые графически определяются как точки пересечения прямой f/p/Bi с графиками периодической функции ctgp (рис. 2.1). Видно,, что эти корни расположены в интервалах (О, л/2), (л, Зл/2) и т. д., т. е. корень р„ лежит в интервале {{п - \)п, (2п - ~1)п12)].

В ойцем случае нелинейной функции Ми) уравнение (2.14) может иметь любое число корней, поэтому во всех численных методах предпачагается, что ожио найтн интервал [а, Ь], содержа-JJH только одни искомый корень. Для оиска этого интервала в случае, когда телТ" /(М) непрерывна, последова- Ряс. 2Л

ато* числяются значения функции

Деля Ро-ониых через равные интервалы на осн р. Это ных "°Р будут найдены два последователь-

кн ря ДрО и f (р*"), имеющие противоположные зна-ные 3?" точках р(> и р(*+> непрерывная функция [ (р) имеет раз-/(ц)Ао интервале [р(*\ р(+>1 лежит корень уравнения tqijijj " нашем случае такой способ не проходит, так как ctgp имеет разрыва, ио из графических соображений мы и так смогли