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

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том XHF 1982

М 3

УДК 629.735.33.015.4-977

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

Г. Н. Замула, С. И. Иванов, С. Ф. Тесленко

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

1. Обычно при дискретизации задач теплопроводности методом конечных элементов (МКЭ) используют процедуру метода Бубнова— Галеркина на локальных координатных функциях для соответствующей континуальной краевой задачи. Иногда применяют также процедуру метода Ритца для эквивалентной вариационной задачи [1, 2]. В этом случае рассматривают либо стационарный случай, либо осуществляют только пространственную дискретизацию. Дадим несколько более общий подход к этому вопросу, основанный на аналогии задачи теплопроводности и теории упругости и на вариационных принципах теории теплообмена, разработанных М. Био [3].

Определим вначале температурно-механическую аналогию, ограничившись рассмотрением стационарных задач. Уравнение стационарной теплопроводности для трехмерного изотропного тела v можно представить в виде трех групп соотношений, а именно: уравнений баланса тепла (закона сохранения энергии)

div {<?[ + ^ = 0, (1)

зависимостей, определяющих „деформации“ температурного поля через температуру

{а} = grad Т , (2)

и закона Фурье (уравнений состояния) для материала

\я) — к Is}« (3)

Здесь (д)—вектор, противоположный вектору плотности теплового потока, F — объемная плотность тепловыделений, {s}—введенный вектор „деформаций“ температурного поля, >^>0—коэффи-

циент теплопроводности материала.

Граничные условия на поверхности 2 тела запишем в виде:

= Мт {?} к =/» (4)

где й1} “ части поверхности 2, на которых заданы граничные

условия соответственно 1-го и 2-го рода, /—поверхностная плотность подводимого теплового потока, {«} — вектор единичной внешней нормали к поверхности 2.

Как видим, соотношения (1) — (4) полностью аналогичны соответственно уравнениям равновесия, зависимостям Коши, закону Гука, кинематическим и статическим граничным условиям теории упругости. При понижении на единицу тензорного ранга величин температура Т является аналогом вектора перемещений {«} теории упругости, вектор {</} соответствует тензору напряжений [з/у.], р и У соответствуют векторам распределенных объемных и поверхностных нагрузок. Пользуясь указанной аналогией, можно сразу сформулировать вариационный принцип для задачи стационарной теплопроводности, идентичный принципу возможных перемещений в механике. Из всех возможных распределений температур в теле только для действительного распределения, удовлетворяющего закону сохранения, выполняется соотношение

Здесь возможные распределения температур непрерывны и принимают заданные значения g(х) на

ода= о 'Ц)с1у — ^{ог}т{(у}. (6)

Действительно, с использованием формулы Гаусса—Остроград-ского и (2)—(6) соотношение (5) преобразуется к виду:

— Г(сНу {д\ + Р)йч + Т({п}т{д} ~ /)(12=0,

с е»

откуда в силу произвольности о Т приходим к уравнениям (1) и граничным условиям (4) на 22.

Важно отметить, что в этом доказательстве не используются формулы (3), т. е. принцип и соотношение (5) справедливы при

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

{<?} = [£] И, (7)

где [/^] — I, / = 1, 2, 3 —симметричный тензор коэффициентов

теплопроводности материала, X, Х^ — могут зависеть от температуры и т. п.

Вариационное уравнение (5) справедливо также при любом виде функций / тепловых „нагрузок“, включая их зависимость, линей-

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

Ьйе = РЪТ, Ы, = /ЪТ, ¿Г=^-{е)т )<?!==^-{е)т [О] (,|>0, (8)

и из (5) приходим к формулировке принципа, аналогичного вариационному принципу Лагранжа теории упругости

о (т — и) = 0, (9)

где функционал от Т

J = т — — ир— | <тФо — Ирйъ —2— (10)

аналог полной потенциальной энергии системы, ч® — аналог удельной энергии деформации тела. При этом в стационарной точке функционала, удовлетворяющей формуле (9), достигается минимум.

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

Р = /=а(Те~ Т) +5?пад — ьс0(Т‘— аТ— гСцТ4, .

где Г,/—заданные функции координат х = (хи х2>

£>0—коэффициент теплоотдачи и степень черноты поверхности тела, не зависящие от температуры, с0 — постоянная Стефана—Больцмана. В этом случае получаем

ар=ТГ, и/=тТ-^-^ (11)

^ о

и после подстановки (8), (11) в (10) приходим к следующему выражению для функционала У:

J = -1|"(е)т [ ТГйъ-(■ Т(7- ~ Т-'-2 Г‘)<Й , (12)

V V -а

в котором {е}, {<?} определяются через Т по формулам (2), (7).

^ у

Учтем теперь нестационарный член с т — уравнения теплопро-

дt

водности, где с ?>0— удельная объемная теплоемкость материала. По аналогии с правилом д'Аламбера можно считать соотношения (1)—(6) справедливыми и для этого случая при учете дополнитель-

(3 т*

ных объемных тепловых „нагрузок“ — В частности,

д(

вариационный принцип (5) для задачи нестационарной теплопроводности, в том числе нелинейной, запишем в виде:

Мт{?})*>- ¡¡ьт/ли=о, (¡3)

V V 8^

совпадающем с полученным М. Био [3] (в несколько другой форме и трактовке) дополнительным вариационным принципом. Из (13), так же как в [3], после введения обобщенных координат у1 (/), ¿=1, 2,. . Т = Т (уи у2, . . у*, х, ¿)> Для случая не зависящих

от температуры коэффициентов теплопроводности приходим к уравнениям для неизвестных уг

(М)

ду1 ду1

аналогичным уравнениям Лагранжа в механике (у* — — ).

(И

Здесь •ш=\<тс1‘0-—введенный выше аналог потенциальной энер-

V

гии деформации тела:

<15>

внешние обобщенные тепловые „силы“;

■ “*>0~ (16)

V

аналог диссипативной функции.

Для случая потенциальных тепловых „нагрузок“ и

уравнения (14) преобразуются к виду

^ + ^- = 0. (17)

¿Ус

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

Остановимся еще на двух вопросах, важных для расчета сложных составных авиаконструкций. Первый из них касается необходимости для некоторых элементов использования приближенных расчетных моделей (термически тонкой оболочки, теплового стержня) с понижением пространственной мерности задачи. Вариационные принципы записаны нами в виде, удобном для такого упрощения. Рассмотрим, например, тонкую оболочку переменной толщины 8(£), £ = (£), 12)> показанную на рис. 1, в ортогональной криволинейной системе координат ?2> $з« Функционал /(12) при предположении постоянства температур по толщине оболочки (в направлении 13) преобразуется к виду

- [ Т (7- -1 Г)еИ, (18)

где {г}—двумерный вектор (2) „деформации“ температурного поля в оболочке, {<7} —двумерный вектор, выражающийся формулой [7|

при [£)] = (8Х^], г, у — 1, 2, р = о F — Ьc ^ ~, / =8/, а = оа. Теплофизические свойства и тепловые „нагрузки“ приняты постоянными по толщине оболочки, 3 = 0 на торцах пластины вдоль контура ¿2 срединной поверхности 9 оболочки. Можно убедиться, что с помощью вариационного принципа (9) из (18) получаются уравнения

нестационарной теплопроводности с учетом внешнего излучения для оболочки

Здесь все функции и дифференциальные операторы определены в двух криволинейных координатах 53 срединной поверхности оболочки, краевая задача для Т (I, Ь) должна быть дополнена соответствующими начальными условиями, /о, ¡хо представляют собой суммы соответствующих величин по двум сторонам (внешней и внутренней) поверхности 2.

Второй вопрос касается необходимости учета в ряде случаев контактных термических сопротивлений в соединениях конструктивных элементов. Получим расширение рассмотренных вариационных принципов на этот случай. Пусть тело V состоит из двух взаимосвязанных частей ги1 и VI, разделенных поверхностью контакта 23, на которой имеется термическое сопротивление

и, следовательно, возможен разрыв температурного поля ДГЛ = = 7+ — Т-, где Т- и Т+—предельные значения температур на £3 слева и справа от Й3 (рис. 2). Введем в рассмотрение плотность теплового потока цп через поверхность контакта, связанную, как известно, с перепадом А Тп формулой

X.

Рис.

be, div[D]grad Т -f- оF -\-fs— aQT —

и граничные условия 3-го рода на

{п)ЦО] grad Г = /~аГ.

qn — А АТп = А(Т + Т_).

(19)

Выделяя из областей Vх, V2 малые объемы Д^,, Дг>2, прилегающие к 23, как показано на рис. 2, и выполняя в первом члене соотношения (13) интегрирование и предельный переход при Дт^, А — 0, получаем вместо (13) следующую формулировку вариационного принципа — аналога принципа возможных перемещений в механике

[ («г + [ш^®-С£77ай+[8(Г+-Г_)?„<*2=0, (20)

где | | | | дп выражается формулой (19).

Рис. 2

Легко убедиться, что с использованием формулы Гаусса—Остроградского для каждой из областей vi} v2 соотношение (20) преобразуется к виду

v J ЪТ,(С^~div {9,}_^ + 2 |зг({/г)т{?,) -f)dQ +

I =1 Vi i- 1

+ J ВТ’, ({«)r{?.)-?n)rf2-j гГг({*Н?2)-<?„)<* 2 = 0,

О о

откуда в силу произвольности оТг = оТ\х^ b7'2~oT\X£Vt получаются обычные уравнения теплопроводности и граничные условия на 22/ Для каждой из областей vif i ~ 1, 2 и условия сопряжения на Q3

WPki} = {п} т{?2} =Яю

где

Ш = №А &rad 71;, /= 1, 2, qn-=\ (Т2 — Тг).

Вариационный принцип (9), (10) и уравнения (14)—(17) Лагранжа также сохраняют силу и получаются из (20) при оговоренных выше условиях и дополнительном требовании независимости термического сопротивления от температуры. Учет контактного термического сопротивления сводится лишь к добавлению в выражении для функционала w члена

¡0ДОn = J®„nd2 = -i-J (r+-r_)9„dS = -l-J.V(7+-r_)y2>0. (21)

£‘з О 8

. Понятно, что все сказанное переносится на случай любого числа А/ взаимосвязанных частей vit i = 1, 2, . . TV, тела v, разделенных произвольным количеством частей 28;., У= 1, 2, . . п поверхности контакта.

2. Перейдем теперь непосредственно к дискретизации задач расчета температурных полей в сложных пространственных конструкциях методом конечных элементов, ограничившись вначале случаем не зависящих от температуры коэффициентов теплопроводности, теплоотдачи и степени черноты. Скобками { }, [ ] ниже

будем обозначать, в отличие от п. 1, вектор-столбцы и матрицы произвольного размера соответственно. Разобьем конструкцию (тело V), как обычно, на конечных элементов так, чтобы контактные поверхности £3 служили границами раздела элементов. В

Рис. 3

качестве обобщенных координат у1 (¿) возьмем значения температуры в узлах элементов, для каждого из которых можно записать

Т' = [о{х)У\уУ, (22)

где [V (*)]' = [®/ (.*), ъгт{х)у...]— матрица заданных непрерывных

(У/1

функций формы г-го элемента; {у}'— \ут } — вектор-столбец температур в узлах элемента; у, т, . . . — номера узлов элемента (рис. 3).

При этом во всех узлах, лежащих на контактных поверхностях, имеются не одно, а два значения температур ут_, ут+, соответствующих смежным элементам, лежащим слева и справа от Q3.

Теперь воспользуемся уравнениями Лагранжа (17), получив предварительно выражения для входящих в них функционалов /г {формула (16)], J [формулы (12), (21)] через уу1 с помощью представления (22) и замены интегралов суммами по конечным элементам

к = У

Г=1

= 2 л

Г=1

где в кг, ]г интегрирование осуществляется по соответствующим объемам 1)гу граничным и контактным поверхностям &г2, конечных элементов. Для трехмерных элементов получаем:

б'=4-{у)'тМ'М'.

[с]г= | [х/]гТ С ? [и]г (¡V —

(23)

хматрица теплоемкости элемента;

У = - иг = -1- [у}'* щг {у}' - {у Г {/К + Ег >

где

[*|- = + [Л,]'- | [Я]'т [О] [Я]'Л> -I-1 №т * ММ9 ,

®2

|/Г= {/[)' + !Л}'= [ М'т?*-+[ М'т/й 2, Е'= р-^{Т'Г л о.

(24)

—соответственно матрица „жесткости“, вектор тепловых „нагрузок“ и член функционала, связанный с излучением, причем

¿>>с1

[В]'=И!И*)]'( И] = {£[;

_д' дх3)

«4>»=-у((У)?- |У)-Т) 1*.1' «УГ+ - !У(~), [А,]'— f К+Т Л М+<< 2=| М-Т А И-а 2 -

где

(25)

дополнительный член, связанный с проводимостью Л через поверхности контакта. При этом важно отметить, что пара индексов г, г' здесь пробегает значения номеров всех конечных элементов, соприкасающихся по 23, а функции формы должны обязательно обеспечивать непрерывность температуры на границах элементов при ее непрерывности в узлах, так что [*»]+ = [<о]1 на Йз = 2з ; индексы „ + “ и „—“ у матриц и векторов означают выделение подматриц и подвекторов, относящихся к узлам элементов, лежащим на 23.

Осуществляя суммирование по всем конечным элементам и перегруппировывая члены, приходим к выражениям

А=уЫтМ(А ^-¿-Мтт [у]-{уг{/) + £, (26)

где [с], [к\ — суммарные матрицы теплоемкости и „жесткости“ конструкции; {у}, {/}—векторы искомых температур в узлах и тепловых

Л'э

„нагрузок“ конструкции; £= — член, связанный с излучением.

Г=1

Подстановка (26) в (17) дает систему N обыкновенных диффе-

ренциальных уравнений для определения {_у (¿)) в модели

И М + [*] (у) = (/). (27)

где в правой части учтено внешнее теплоизлучение, так что

{/)'={/}'+ {/з)г» причем

Уравнения (27) нелинейны, так как [с]г и {/}г согласно предыдущему и формулам (23), (28), (22) могут зависеть от {у)г\ начальные условия для них:

{У (0)} = Ы- (29)

Полученные соотношения справедливы, строго говоря, при условии Й = Й2, обычно выполняющемся для задач расчета температурных полей в авиаконструкциях. В противном случае их также можно использовать, дополнительно преобразовав строки и столбцы матриц и компоненты векторов уравнений (27), (29), соответствующие узлам на поверхности1 £1} как это обычно делается в методе конечных элементов [1].

Остановимся кратко на особенностях, связанных с понижением мерности отдельных элементов и с учетом контактных термических сопротивлений. Первая особенность, а значит, и пространственная дискретизация задач расчета температурных полей на разномерных конечных элементах, легко разрешима на основе описанного подхода. Действительно, переход к одномерным конечным элементам теплового стержня, двумерным—термически тонкой оболочки и т. п., можно выполнить в соответствии с и. 1, лишь несколько преобразовав выражения для функционалов &, J на соответствующих частях конструкции. Например, конструкция может состоять как из трехмерных элементов, так и из связанных с ними конечных элементов тонкой оболочки. Для последних узлы выбираются на срединной поверхности 2, а аппроксимация (22) осуществляется по ее координатам 1=(^, ?2), т. е. Тт = [г> (;)]г{.у}г, причем так, чтобы выполнялось условие непрерывности поля температур на стыках разномерных элементов. Описанная выше процедура пространственной дискретизации, но с учетом того обстоятельства, что на элементах оболочки для J справедлива формула (18), а выражение (16) для /г преобразуется к виду:

л = ~ а\ха\г, (30)

приводит к тем же уравнениям (27), (29) для всей конструкции с несколько другими формулами для конечных элементов оболочки

[С]'_| 1*]г=| ([в]'т \о\ \ву+

2Г 2Г

-ь [г>]'т Яо \ю}г) ¿2 + | №та \v\dL ;

_

{/}'=[[фГ(8/7 + /«)(/9+ СМгТ7^. [Я]“[*М* *■ /=Ь 2, I

Т 2

• 1 д

д Ъ

1 д •

. ¿2 ді2

Здесь ¿1? ¿2 — коэффициенты первой квадратичной формы срединной поверхности оболочки. Аналогичный результат получается при добавлении одномерных конечных элементов и в общем случае конструкций из элементов разной мерности.

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

I где г"—

М-г

номер элемента сопротивления на общей границе смежных элементов г и г'. При этом все узлы, лежащие на 23, будем считать двойными и имеющими разные номера т_, т+ и т. д. слева и справа от 23 (см. рис. 3). Полученные ранее выражения для тгЛОп тогда преобразуются к виду

=~-ь-гт ту г,

при векторе обобщенных координат в узлах {у}г" —

в котором \к]г"

—матрица „жесткости“ элемента тер-

Ы -[*.]' -Ыг Ыг

мического сопротивления — суммируется с матрицами „жесткости“ остальных конечных элементов в процессе формирования [&] по обычным правилам. Матрица теплоемкости и вектор тепловых „нагрузок“ этого элемента — нулевые. В каждом узле конечно-элементной модели конструкции при таком подходе, как и в случае отсутствия контактных термических сопротивлений, имеется лишь одна обобщенная координата — значение температуры в узле и, следовательно, №= /Уузл, где ЛГу3д — число узлов модели.

Для решения задачи (27), (29) используются различные шагово-итерационные схемы. Интегрирование в (23) — (31) для конечных элементов может осуществляться численно, в частности, это необходимо при учете излучения. Указанные вычислительные задачи будут рассмотрены нами в отдельной статье.

При учете зависимости коэффициента теплоотдачи я и степени черноты е от температуры для формирования соотношений МКЭ воспользуемся уравнениями Лагранжа в форме (14), (15) вместо (17).

7—.Ученые записки ЦАГИ“ № 3

97

Легко видеть с учетом (22), что конечные уравнения (27) для модели и выражения (23), (24), (28), (30), (31) для элементов при этом не претерпят изменений, в них нужно подразумевать лишь зависимость соответствующих членов от температур {у}г в узлах. К аналогичному результату приходим и в наиболее общем случае зависимости от температуры коэффициентов теплопроводности и контактного термического сопротивления для этого нужно воспользоваться вместо уравнений Лагранжа общим вариационным принципом (20).

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

(Г)4=Г + 4Г3 (ТГ—Т) = 4ТЧ'Г — ЗТ* ,

где Т — температура линеаризации, зависящая или не зависящая от х, г, причем }Тг—7| <С Г.

Подстановка этого соотношения и (22) в (28) после преобразования дает

{/„)'=-[*.!' М'-М/зГ,

где

(Й,1'=]>Г4 гсаТ°М<Ш, 17„1'=)'1яГЗгс9Г<*9. (32)

Таким образом, уравнения (27) сохраняют свой вид при замене (/3}г на {/г3}г и добавлении в матрицах жесткости конечных элементов члена [64]г согласно (32), аналогичного члену [¿2]г от конвективного теплообмена. В простейшем случае температуры линеаризации для элементов задаются априорно и задачи (27), (29) решаются, как обычно. При использовании шагово-итерационных схем линеаризация может осуществляться относительно значений температуры на предыдущем шаге расчета либо на предыдущей итерации.

ЛИТЕРАТУРА

1. Зенкевич О. Метод конечных элементов в технике, М-, „Мнр“, 1975.

2. Постнов В. А. Численные методы расчета судовых конструкций. Л., „Судостроение", 1977.

3. Б и о М. Вариационные принципы и теория теплообмена. М., „Энергия“, 1975.

Рукопись поступила 121X11 1980 г. Переработанный вариант поступил 2! VI 1980 г.