Рекомендации по априорному выбору параметров дискретизации при решении уравнения теплопроводности методом конечных элементов Текст научной статьи по специальности «Физика»

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том XXIV 19 93

№ 1

УДК 629.7.015.4—977 519.63 : 512.643

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

А. В. Микневичюс

Разработаны в виде номограмм рекомендации по априорному выбору характерного размера конечных элементов (КЭ) А и шага счета А^, обеспечивающих при решении уравнения теплопроводности методом конечных элементов (МКЭ) уровень погрешности численного решения не выше заданного. Результаты, полученные при исследовании модельных задач, обобщаются на более сложные конструкции и условия теплообмена, приводится пример их использования в расчете теплового состояния реального объекта.

При использовании МКЭ в комплексах вычислительных программ (КВП) одним из ключевых для пользователей является вопрос выбора параметров дискретизации /г и М, позволяющих достичь нужной точности численного решения при допустимых затратах ресурсов ЭВМ. В основе разработки рекомендаций по априорному выбору Л и лежит определение связи между точностью численного решения и параметрами дискретизации. В настоящей работе эта связь устанавливается на основе обработки результатов вычислительного эксперимента с помощью аппарата анализа размерностей (см., например, [1]). Вычислительный эксперимент проводился с применением КВП «Транзит-2» [2], созданного на основе работы [3]. Рекомендации разработаны в виде номограмм для набора модельных задач, имеющих аналитические решения, приведенные в работе [4]. Показано, что погрешности численного решения при тепловых расчетах реальных конструкций будут заведомо ниже погрешностей, определенных из номограмм, если Л и А/ выбраны в соответствии с правилами, предложенными в настоящей работе.

При подборе модельных задач и исследовании погрешности численного решения с использованием программ [2, 4] вводились следующие ограничения:

— рассматривались одно-, двумерные задачи, имеющие аналитические решения. Геометрия расчетных областей и границы теплообмена показаны на рис. 1, а, б, где Ь — характерный размер расчетной области;

— исследовались только КЭ типа «прямоугольник». Разбиение расчетных областей на КЭ равномерное. На рис. 1 ,а, б показаны примеры разбиений на КЭ двумерной и одномерной областей расчета с числом разбиений N=4 по осям х, у. Цифрами без скобок обозначены номера узлов, цифрами в скобках — номера элементов;

— предполагалось, что внутренние источники тепла отсутствуют, начальное распределение температур равномерно и известно. Рассматривались стационарные и линейные по времени граничные условия 1, 2, 3-го родов;

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

12 — «Ученые записки» № I

177

12 3- , Л. ,

(V 12) <1

(15) (1В) л;

Г У

а)

б)

(1)

(2)

(3) •ч

Г?)

ф-У/’/’У/ 10

Рис. 1

1. Рассматривается смешанная схема дискретизации применительно к уравнению теплопроводности:

V 7' ~ а ДГ (1)

с начальными и граничными условиями (НУ) и (ГУ). В уравнении (1) Т —температура, а='к/рс — коэффициент температуропроводности, р — плотность, с—удельная теплоемкость. После интегрирования дифференциального уравнения в частных производных (1) по координатам МК.Э имеем обыкновенное дифференциальное уравнение (ОДУ) в матричной форме:

[С] + [АГ1 {Г} (2)

где [С], [/Г]—матрицы теплоемкости и теплопроводности; {Г}, {дТ1д(}, {/?) —векторы температур, градиентов, температур и сил. Применение к ОДУ (2) одношагового метода численного интегрирования с параметром неявности 0 приводит к системе линейных алгебраических уравнений:

г П + 1 ;_ I Т)П

[С] ——А(П ' 1 + 0[/С]"+1 {Г}"+1+ (1 — 0) [/С]" {Т1п = 6 1Р}п+1 + (1 — В) {/?}\ (3)

где индексы п и п+1 относятся к п-му и л+1-му шагам счета по времени соответственно. Численная схема (3), применяемая для решения уравнения (1), обладает, как следует из работ [3, 5—7], свойствами аппроксимации, сходимости и устойчивости. В настоящей работе исследуется ее точность.

Погрешностью численного решения е=\ Т—Г| будем называть разность между истинным решением уравнения (1) Т и его решением, полученным по численной схеме (3). Основной вклад в величину е вносят:

— погрешность аппроксимации ен, обусловленная природой функций формы, используемых внутри каждого КЭ;

— погрешность усечения б( от пренебрежения в численной схеме (3) членами высших порядков при разложении векторов {дТ/дЬ}, (Т), {Т7} в ряды Тейлора;

— погрешность округления, возникающая в связи с тем, что на любом этапе вычислений в ЭВМ можно запомнить только конечный объем информации. Сравнение результатов расчетов, полученных при решении уравнения (1) с помощью [2] с одинарной и двойной точностью, показывает, что в исследуемом в данной работе диапазоне к и М вклад погрешности округления в общую погрешность е пренебрежимо мал по сравнению с вкладами погрешностей аппроксимации и усечения; это — погрешности, связанные с неверной аппроксимацией границ теплообмена, неточным заданием или грубой интерполяцией ГУ и т. д.

В настоящей работе исследуется влияние на точность численной схемы (3) только погрешностей аппроксимации и усечения.

Оценка погрешности аппроксимации дана во многих работах, в частности в [5—7]. Авторы этих работ в случае равномерного разбиения расчетной области на КЭ приходят к неравенству для нормы ЦелII:

\\еН\ <Мк\ (4)

где Л=Л//.—относительный характерный размер КЭ.

Если дифференциальное уравнение представляет собой двумерное уравнение теплопроводности, норма Це/,11 определена как Н'—норма, т. е.

II II = К] |4 + (йгд/Л*)3 ■}■ (дещду)* йх йу

и используются билинейные КЭ, то V = 1, а коэффициент М может быть определен для задач, имеющих аналитическое решение.

Из работ [5—7] следует, что погрешность усечения при интегрировании уравнения (2) одношаговым методом с параметром неявности 0 имеет порядок

\е1\~0(М) при 0 >1/2,')

\е(\~ 0(Д<*) при 0=1/2.)

На основе соотношений (4), (5) может быть проведена оценка погрешности численного определения е, если она полностью определена лишь погрешностью аппроксимации или лишь погрешностью усечения. При рациональном выборе параметров дискретизации погрешности аппроксимации и усечения сопоставимы по величине. Алгоритма, позволяющего априорно оценить погрешность численного решения в этом случае, в данное время не существует. Попытка его разработки осуществлена в этой работе.

Применение аппарата анализа размерностей к уравнению (1) при различных ГУ с учетом оговоренных выше ограничений приводит к критериальным уравнениям, сведенным в таблицу, позволяющим исследовать поведение погрешности численного определения е.

ГУ ГУ-функцин времени Критериальное уравнение Рису нок

1-го рода Т s — const К = ?! (/?!. /"о. Ю 3, а

2-го рода q — const К - Ъ №, Р<и Ы) 3, б

3-го рода a = const ^ = Тз№, Ри, N. В!) 3, в

Tc — const

1-го рода T§ = bt &1Л — 1» М) 3, г

2-го рода q = b2t Лм — Фа о» М) 3, д

3-го рода a = const ^?зл = Фз (/Ъ Л7. В1) 3, е

О .11 w .

В таблице использованы следующие обозначения: Г0—начальная температура; Те — заданная на границе тела температура при ГУ 1-го рода; ^ — плотность теплового потока; а — коэффициент теплоотдачи; Тс—температура среды; t — текущее время; b^, Ьг, Ь3 — коэффициенты, характеризующие скорость изменения ГУ по времени; —сеточный критерий Фурье; р0 = а(*/Ь2— критерий Фурье;

Ri = шах /=1, N,

уз

Т, - Т0

R2— max|

\е |1

ЖIJi ’

R3 — max

относительные погрешности для постоянных ГУ I, 2, 3-го родов соответственно; А^уз — число узлов КЭМ; (* — время, начиная с которого уровни относительных погрешностей R1, /?2> Rз становятся ниже заданных;

R1

< = 1, лг,

уз

еа ) п i \е 1 i еа

bjJ |.> “зл = max j b2Lt j . ^зл = тах| ь7ц

относительные погрешности для линейных по времени ГУ 1, 2, 3-го родов соответственно; N = Ь/к — число разбиений расчетной области на КЭ; В1 = «¿Д — критерий Био; у2, уъ +1. 'Ь — функции, определяемые с помощью вычислительного эксперимента,

В таблице приведены зависимости ГУ от времени при />0. Предполагается, что при КО выполняются условия:

Г5 = Т0, д = О, Гс = Т0.

При стационарных ГУ зависимость погрешности численного решения от времени характеризуется всплеском ее величины в области малых времен с последующим асимптотическим спадом (рис. 2). При линейных ГУ величина погрешности е монотонно растет и с течением времени становится постоянной при ГУ 1, 3-го родов, либо линейно изменяющейся по времени при ГУ 2-го рода. Таким образом, зависимости е . от Ь при постоянных и линейных ГУ качественно отличны друг от друга. Критериальные уравнения из таблицы, полученные для постоянных ГУ, позволяют определять время ¿*, начиная с которого уровни относительных погрешностей Я2, Яз становятся ниже заданных, а критериальные уравнения, полученные для линейных ГУ, используются для определения стабилизировавшихся, т. е. не зависящих от времени, относительных погрешностей л, ^2 л, Яз л (заметим, что в знаменатель выражения для определения л входит ¿). Задача определения погрешностей численного расчета решалась в одно-, двумерной постановках. Критериальные уравнения для каждого типа ГУ в одно-, двумерной постановках имеют одинаковый вид, но в первом случае число КЭ и КЭМ равно Ы, а во втором — ./V2.

На рис. 3 приведена часть полученных результатов. В таблице дано соответствие между типами ГУ и номерами рисунков. Рис. 3, а, в соответствуют двумерной постановке (см. рис. 1 , а), а рис. 3,6, г, д, е — одномерной (см. рис. 1,6). Для ГУ 3-го рода результаты, приведенные на рис. 3, в, е, соответствуют критерию Био В1 =10 Границы варьирования остальных параметров приведены на рис. 3, а—3 е.

2. Необходимо дать некоторые рекомендации, позволяющие использовать номограммы, приведенные на рис. 3, в вопросах априорного выбора параметров дискретизации при расчетах тепловых режимов реальных конструкций. Эти рекомендации основаны на следующем предположении, требующем доказательства: кондуктивный отвод тепла с границ, показанных на рис. 1, а, б теплоизолированными, снижает уровень погрешности численного решения. Без ущерба для общности проведем доказательство для постоянных ГУ 2-го рода одномерной постановки задачи. Рассмотрим КЭМ с числом разбиений Д% и характерным размером ¿1 по оси у (рис. 4, а). Для определенности примем <71>92>0. Качественное распределение температур по координате в произвольный момент времени представлено на рис. 4, г. На расстоянии

дТ

от начала координат выполняется условие теплоизолированное™

ду

У=и

0.

КЭМ на рис. 4, б с числом разбиений N2 и характерным размером является частью КЭМ на рис. 4, а. Для нее применимы результаты, приведенные на номограммах

дТ I

рис. 3. КЭМ с условием теплоизолированности =0. характерным размером

Ьз(Ьз<Ьг) и числом разбиений N3, показанная на рис. 4, в, является частью КЭМ, представленной на рис. 4, б. Для доказательства выдвинутого выше предположения необходимо показать, что уровень погрешности численного решения для КЭМ на рис. 4, в выше уровня погрешности для части КЭМ на рис. 4, б, включающей в себя

Lг 0,6

ол

0,2

О

Ч)

* Т0-0,003

□ о, од

а,їв

п3=0,05 ; BL=10 -

-----□ —Q

12 N

Рис. 3

а)

tii

\Ь

(1) • m •

• m • «41 •Ч №

dt,у (Nt)

••• ^ Ч) q-0

(1)

4 « • •

(К)

S)

q>0

Рис. 4

КЭ с номерами от (1) до (N3). Перепишем критериальное уравнение из таблицы для постоянных ГУ 2-го рода в виде:

шах

4\h

) lat* аМ Л

= F' W ’ ~W' 17» ) ■ 1 = N>3 2

(6)

для КЭМ на рис. 4, б с числом узлов, равным iVy3 2; Хе3

qxh

_ qV*_ сШ_ Ji

ч ' Л2 ’ ¿з )' г ’ у:> 3

(7)

для КЭМ на рис. 4, в с числом узлов, равным iVy3 з

Сравнивая критериальные уравнения (6), (7), можно сделать вывод, что

тах{|е2|}; может отличаться от шах{|е3|}; только за,счет различия относительных размеров КЭ h/Li и hjL3. Если находиться в рабочем диапазоне параметров разбиения, где с уменьшением относительного размера КЭ уровень погрешности снижается (для нормы Н1 скорость снижения без учета влияния погрешности усечения иллюстрируется неравенством (4)), то, учитывая неравенство L3<L2, приходим к неравенству max{|e3|}i>max|{e2|}t-, справедливому как для всей КЭМ на рис. 4,6, так и для ее части. Таким образом, предположение доказано. Его косвенным доказательством может служить также тот факт, 'что кондуктивный отвод тепла с границы y — L3 (рис. 4, в) снижает градиент температуры по координате в области расчета. В случае (]z<0 характерный размер L2 не определен. Однако для части КЭМ на рис. 4, а, включающей в себя КЭ с номерами от (1) до (N3)', принципиально ничего не меняется: подвод тепла за счет ГУ с одной стороны и кондуктивный отвод тепла с другой стороны. Следовательно, и в этом случае погрешности в КЭМ на рис. 4, в будут выше, чем в части КЭМ на рис. 4, а, включающей в себя КЭ с номерами от (1) до если \qi\^>\qz\.

Доказанное предположение и некоторые другие рассуждения позволяют рекомендовать Следующий выбор правил по использованию номограмм на рис. 3 при выборе параметров дискретизации для численного расчета МКЭ теплового состояния реальных конструкций. Если тепловые потоки, подводимые к конструкции в направлениях осей координат х, у, сравнимы по величине, то используются номограммы, построенные для двумерной постановки задачи в противном случае — для одномерной. При этом всегда в качестве параметров ГУ (Ts, а, Тс, q), входящих в критериальные уравнения из таблицы, используются параметры, обеспечивающие максимальный тепловой поток на границе конструкции. Если конструкция неоднородна и содержит ряд включений из различных материалов, то каждое включение или часть его представляется в виде конструкций, показанных на рис. 1 , а или 1,6, теплоизолированных с одних сторон и подверженных заданным ГУ с других. Важным является вопрос правильного выбора характерного размера L, входящего в крите-

риальные уравнения. При |</i | > |д21 >0 (см. рис. 4) выполняется условие L2>/.3/2. Оценка погрешностей с помощью номограмм на рис. 3 для КЭМ на рис. 4, а может проводиться лишь при условии L<L2. Это условие обеспечивается равенством L = = Li/2. В общем случае характерный размер определяется из уравнения

А=£ш1п/2, (8)

где ¿min—минимальный размер включения в направлении распространения потока тепла. При двумерной постановке задачи Lmin—минимальный из размеров в направлениях осей х, у. В частном случае, когда одна из границ включения теплоизолирована при одномерном распространении теплового потока или две соседние границы теплоизолированы при двумерном распространении теплового потока, L —Lmm. Величина размеров КЭ, определенная из номограмм с учетом уравнения (8), не должна превышаться и в тех частях КЭМ, где размеры включения превышают ¿min> что, в общем случае, может потребовать различное число разбиений на различных участках КЭМ, моделирующей включение. Применение в какой-либо части КЭМ неравномерной сетки разбиения, сгущенной в области максимальных значений градиента температур, позволяет получить точность численного решения не ниже заданной, если число разбиений в этой части КЭМ определяется из номограмм на рис. 3. При определении с помощью номограмм требуемого шага счета для различных включений неоднородной конструкции получаются различные значения Ai, которые могут отличаться друг от друга на несколько порядков. Используя терминологию работы [4], будем называть термически тонкими или термически толстыми тела, перепад температур в которых относительно мал или велик. Поведение погрешности численного расчета в конструкции, состоящей из граничащих термически тонкой и толстой прослоек, в случае когда термически тонкая прослойка толщиной LT подвержена постоянный ГУ 3-го рода и моделируется одним КЭ, а термически толстая прослойка толщиной L теп ■ лоизолирована с тыльной стороны, описывается критериальным уравнением

/Ъ = Л(Яз. Fm N, Bi, S), (9)

(К)Г /-Т

где а = —; параметры с индексом «т» относятся к термически тонкому телу,

остальные параметры к термически толстому. КЭМ, изображенная на рис. 1,6, соответствует рассматриваемой конструкции при УУ=3, если полагать, что КЭ под номером (1) моделирует термически тонкое тело, а остальные — термически толстое. Аналитическое решение такой одномерной задачи при А,т->-оо, постоянных а, Тс без учета теплового излучения и зависимости теплофизических свойств от температуры или времени приведено в работе [8]. Исследование уравнения (9) с помощью КВП [2] и работы [8] при Р3 = 0,05; В! =10; р == 0,01; 0,1; 1; 5 дало результаты, близкие к тем, что приведены на рис. 3, в и .соответствуют уравнению (9) при ¡3 = 0. Кроме того, расчеты показали, что с ростом ¡3 в диапазоне 0<|3<5 погрешность численного решения имеет слабую тенденцию к снижению. Таким образом, наличие в неоднородной конструкции термически тонких тел оказывает слабое влияние на погрешность численного решения. Это свидетельствует о том, что те части КЭМ, которые могут рассматриваться как термически тонкие тела, не должны учитываться при выборе шага счета, а моделировать их следует минимально возможным числом КЭ. Реально термически тонким телам соответствуют те, для которых значения или /?1 л, ^2Л, /е3 л выходят за пределы верхней границы диапазона их изменений на номограммах рис. 3, кроме рис. 3, б. В качестве параметра А/ выбирается минимальный из числа полученных для термически толстых тел.

3. Приведем пример использования номограмм. Объектом теплового расчета является экспериментальная установка ЦАГИ вместе с испытываемым образцом. Установка, как показано на рис. 5, а, включает в себя держатель (позиции 1, 2) и гермоотсек (позиция 3). Гермоотсек представляет собой полый параллелепипед со средней толщиной стенок 5 мм, выполненный из АМГ-6. Роль верхней крышки гермоотсека играет испытываемый образец, состоящий из слоя ТЗИ толщиной 66 мм (позиция 4) и панели, выполненной из Д-16 (позиция 5). Панель толщиной 1,5 мм подкреплена пятью дюралевыми стрингерами. Внешняя поверхность слоя ТЗИ, левая и правая части держателя, а также днище гермоотсека поддерживаются при нестационарных температурах 7\±Дь Т-г + Д2, Т3±Д3, T^±A¡, соответственно, где Д1 = = 34,5 К, Дг = Дз = Д4 = 6 К — предельные погрешности замеров температур в различных температурных интервалах. Зависимости от времени температуры 7\ и давления р на внешней поверхности образца приведены на рис. 5, б. Зависимости от времени

Л)

670

S30

Маисшчамные расчетные j ¿7 температуры

ЪНКк

200

150

100

50

О

)'/гсne ■ JI

шмент^л

х^^-Минансиы

oacvemmt

S)

500

'Мининсибнш расчетные температуры L

то

t.c

Рис. 5

температур Т2, Тз, ТI, приведены на рис. 5, в. Остальные границы конструкции теплоизолированы. Как следует из сказанного выше и рис. 5, а, плотность теплового потока в направлении оси г существенно ниже плотностей тепловых потоков в направлениях осей х, у. Это обстоятельство позволяет использовать двумерную тепловую схему, приведенную на рис. 5, г. Схема соответствует сечению А—А на рис. 5, а. Держатель в нее не включен, так как его влияние моделируется заданием ГУ 1-го рода па боковых стенках гермоотсека. Сравнение расчетных и экспериментальных данных проведено в точке Э (рис. 5, г). Подвод тепла к точке Э осуществляется: кондуктивно от внешней поверхности образца, левой и правой сторон держателя; естественной конвекцией и излучением внутри объема гермоотсека. В расчете учитывается зависимость теплофизических свойств материалов от температуры. Кроме того, зависимость коэффициента теплопроводности материала ТЗИ от давления описывается эмпирическим соотношением

Хр — (А0 А,)

124 + Т 1 + 9- 10-< ——— р

где А.о, Яь КР — соответственно коэффициенты теплопроводности при нулевом, атмосферном давлениях и давлении р.

Панель с подкрепляющими стрингерами представляет собой в направлении оси у термически тонкое тело. В связи с этим в тепловую схему на рис. 5, г она включена как пластина постоянной эффективной толщины, равной 2,2 мм. В расчете эффективной толщины учитывались массы панели и стрингеров.

Предельная погрешность измерения температур в точке Э составляла 6 К. Выберем параметры дискретизации для составления КЭМ. Используем критериальное уравнение из таблицы для линейных ГУ 1-го рода и рис. 3, г. Для выбора параметров дискретизации в части тепловой схемы, моделирующей слой ТЗИ, последний представляется, в соответствии с предложенными рекомендациями, в виде прямоугольной области, подверженной заданным ГУ с одной стороны и теплоизолированной с другой (боковые стенки слоя ТЗИ по условию теплоизолированы). Характерный размер в направлении оси у определяется соотношением (8) и составляет ¿ = 0,033 м.. Параметр /;|, характеризующий скорость изменения ГУ по времени, определим для интервала времени от 700 до 1600 с, где он принимает максимальное значение (рис. 5,6) =

= 1,1 К/с. Зададимся предельно допустимой погрешностью численного расчета е = 6К, совпадающей с предельной погрешностью измерения температуры в точке Э. Осред-ненный в интервале температур от 323 до 1323 К и интервале давлений от 0 до 105 Па коэффициент температуропроводности ТЗИ имеет значение а=10-6 м2/с. Отсю да определяем допустимую относительную погрешность численного решения Я1 Л = = 5-10_3. Из рис. Зг следует, что такой уровень погрешности может быть обеспечен выбором N=4, Ро = 1 • Ю~3, откуда находим Л = 8 10~3 м, Л/=1 с. Допустимый размер КЭ /г и толщина слоя ТЗИ в направлении оси у (она равна удвоенному характерному размеру ¿) задают требуемое число разбиений N = 8.

Вычисления параметров Ri л в направлениях осей х, у для стенок гермоотсека и панели испытываемого образца показывают, что их значения выходят за пределы верхней границы диапазона изменений п на рис. 3, г. Как отмечалось выше, это означает, что рассматриваемые части конструкции являются термически тонкими телами, могут быть моделированы одним КЭ каждая и не учитываться при выборе шага счета. Разбиение по оси х областей КЭМ на рис. 5, г, соответствующих панели испытываемого образца (на 4 части) и слою ТЗИ (на 2 части), использовано не с целью снижения погрешности численного расчета, а для согласования частей КЭМ и удобства сопоставления результатов расчета с экспериментальными данными в точке Э. Расхождение экспериментальных и расчетных данных в рассматриваемом примере обусловлено, в основном, погрешностями: численного решения, замеров темпера-

тур в точке Э и воспроизведения ГУ 1-го. рода в эксперименте. На рис. 5, д приведены максимальный и минимальный уровни температур в точке Э, рассчитанные при задании в качестве ГУ температур 7\+Д/ и Г,-—А;, /=1,4 соответственно. Вертикальные штрихи, обозначающие экспериментальные температуры с предельными погрешностями замеров в точке Э, полностью или частично попадают в расчетную область, ограниченную максимальными и минимальными расчетными температурами. Следовательно, расхождение результатов численного решения с экспериментальными данными при выбранных параметрах дискретизации лежит в пределах погрешностей, обусловленных погрешностями замеров температур в точке Э и воспроизведениями ГУ 1-го рода в эксперименте. Сопоставление экспериментальных и расчетных данных в приведенном примере дает основания полагать, что предложенная в настоящей работе методика выбора параметров дискретизации Применима при определении тепловых режимов реальных конструкций, в том числе и при решении существенно нелинейных задач.

ЛИТЕРАТУРА

1. Гул мал А. А. Введение в теорию подобия.—М.: Высшая школа. 1973.

2. Замула Г. Н., Тес лен ко С. Ф., Микневичюс А. В. Комплекс вычислительных программ «Транзит-2» расчета нестационарных температурных полей и напряжений элементов конструкций летательных аппаратов методом конечных элементов. Специализированный раздел фонда алгоритмов и программ//Информационный сборник. ГосНИИАС.— 1990, № 5.

3. Замула Г. Н., Т е с л е н к о С. Ф. Численная реализация метода конечных элементов в задачах теплопроводности авиаконструкций//Труды ЦАГИ. — 1984. Вып. 2229.

4. П е х о в и ч А. И., Жидких В. М. Расчеты теплового режима твердых тел. — Л.: Энергия, 1968.

5. Ши Д. Численные методы в задачах теплообмена. — М.: Мир, 1988.

6. Флетчер К. Численные методы на основе метода Галерки-иа. — М.: Мир, 1988.

7. Зенкевич О., Морган К. Конечные элементы и аппроксимация.— М.: Мир, 1986.

8. Зарубин В. С. Температурные поля в конструкции летательных аппаратов. — М.: Машиностроение, 1966.

Рукопись поступила 25/1 1991