Применение метода конечного элемента для расчета нестационарных температур в сечении тонкостенных конструкций Текст научной статьи по специальности «Математика»

УЧЕНЫЕ ЗАПИСКИ Ц А Г И

Том XIII 1 98 2 № 1

УДК 536.2

629.735.33.015.4-977

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

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

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

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

ЛЛААДАА

Рис. 1

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

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

1. Рассмотрим вначале задачу определения нестационарных температур в стержне переменной толщины о(х), длины I. Будем считать, что нагревание стержня происходит по боковым поверхностям посредством излучения и теплоотдачи от сред с температурами Те 1 и Те2 и коэффициентами теплоотдачи я2. Пусть эти величины, так же как и лучистый поток С2Л, являются функциями координаты х, отсчитываемой по длине стержня, и времени т. Тепловые потоки на концах стержня обозначим (г) и дг('г), температурные перепады по его толщине отсутствуют. Температурное поле в стержне определяется решением уравнения теплопроводности

^^=^й17+а<(7'.,-2г)+«,(т,г-г) + <ь (1) с начальными

Т {х, 0) = Т°(х)

и граничными условиямм

д Т I / > л-1 дТ ] / \

__&х— =?1(<г), оХ— = я-г (у;

дх |х=о дх \х~1.

здесь р — плотность материала, с, X — удельная теплоемкость и коэффициент теплопроводности, в общем случае зависящие от температуры.

Разделим стержень на А7— 1 отдельных элементов, для аппроксимации температур внутри них будем использовать линейные функции:

Г--=Ы[Т, + №тш,

Т{, Г,-+1 — температуры узлов, М, N\ — функции формы г-го элемента. Поток <3Л, коэффициенты теплоотдачи аь а2 и температуры восстановления для простоты будем считать постоянными в пределах каждого элемента. Для получения пространственных аппроксимаций с помощью конечных элементов используем метод Бубнова — Галеркина. Для этого рассмотрим А/'-мерное пространство в котором в качестве базисных функций выберем и кусочно-линейные функции, задаваемые соотношениями

О при А?-1 при <.*<**, при Хк^Х^Хк+1.

В любой момент времени приближенное решение задачи будем представлять в виде

N

Т(х,

г

Согласно методу Бубнова — Галеркина функция Тк(\) определяется-нз условия (граничные условия естественные):

■ $А(Т)^(х)<1х = 0, к= 1, ...,М

(2)

где

А (Т) = Вер

дт

д гх-^ + (а, + а2)Г-<Э,

Обозначив

дт дх дх

(2 = с?л -Ь <4 1 + а2 Те 2

р,=/,(«,+«л

(/,, 8; — длина и средняя толщина /-го элемента, <3{-— опреде-

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

’.(2^ + ^)-|Ч(7'2-Г1) + М2Г1 + Г2) + -^--<?,=(>,

дТ

У/-1 ~к^ + 2(у/-1 + *<) —|Г + (^ — 7"/_1) — гг)+

+ (3,_, Г,., + 2 (Р,_, + И Т, + р, Гц., + + =0)

дТ,

дГ

1+1

(*-2, ЛГ-1);

>ЛГ-1

(—^ + 2 ■ '] + Н-л/—>(Т/г — Тп-1) +

+ Р,у-1(/-д'-1 + 27») + ----92в0;

которую удобно переписать в матричном виде

(3)

в-^ + (Д + о)г=/, 0<т<т0, Г(0)=Г,

где

Н-1 — Iх! 0 0

Л = 5А1 У-1 + Н У* 2 ... 0

0 0 0 ••• Н-Л/-1

2у1 *1 0 . . . 0 Т\.

VI 2.(*1 + *2) 0 т2

0 0 о • 2>дг_1 Ты

2Р. ?1 0 * 0

о = Рх 2 (?! + ?,) Рг • 0

0 0 0 . . 2(3^-

/

?1

0»/1

?лг

*Л/-1

(4)

4—«Ученые записки» № 1

49

Таким образом, мы пришли к задаче Коши для системы линейных дифференциальных уравнений, причем количество уравнений и неизвестных зависит от разбиения стержня на элементы. Пространственная часть аппроксимации (5) при а1 = а2=0 полностью совпадает с аппроксимацией, получаемой методом балансов [1, 2]. Члены, связанные с конвективным теплообменом и неста-ционарностью, отличаются, в частности, недиагональностью матрицы В.

На отрезке 0<!х<;т0 введем равномерную сетку с шагом Дх:

ю, = |<Су=уДт, у =0, 1, . . . , /о, Дт = ~^-|

и будем рассматривать функции у (х), <р(т) дискретного аргумента х=уДт с значениями в соответствующие Т (х), /(т), так что для всех 'с=уДх^ш..

Заменив дифференциальные операторы конечно-разностными, рассмотрим зависящее от параметра а семейство двухслойных схем с весами

Ч- (-4 -ь г>) (ау + (1 — а)3/) = ср, 3>(0)=У°, (5)

где у — вектор, соответствующий вектору температур Т в момент л л

лДх, у — его значение при х = (я+1)Дх, у^(у — у)/Дх, После введения обозначений '

е1 = ех (а) =— 2v1 — аДх + 2М,

=~ ем(а) = — 2ул-1 — аДх ({4дг-1 -(- 2Рлг-1), е1 = е1 (в) = — 2 (у 1—1 +•**) — о Дх (р.._1 + + 2^г_1 +2^), г = 2, . . . , IV— 1

а1 = ах (о) = 0, ^ = аДо) = у*_1 — аДх^.^ — р._2), г = 2, . ... Л/; • ^1 = е,(а — 1)У1 — а2(а - 1)Уг + Дх(— Я\ + ,

Ь = — -1 )у!-г + е{(о~ 1)^ — а1+1 (а - 1)у1+1 4- Дх <**~1

I — 2, ЛГ-1;

Удг = — (а — 1)^лт-1 + ен (а — 1) Ул? + Дх (— ^ ®N-1^-1 I

система (5) может быть переписана в виде, удобном для ее решения методом прогонки

Л Л

азУ2 — ?1»

А А А

аі Уі—\— е1уі + аі+і уі+і = — оі}

А . А

аыуы-1 — Єн Уії =—-<рл'.

(6)

Для исследования устойчивости семейства схем (8), использовав равенство уі=“ + *Уі, перепишем (5) в каноническом виде [3]:

в\Уі + {А + £>)У = <Р, (7)

где Я^В + аД^Л +£>). Операторы Л, В, £>, Вх — самосопряженные; В, £) — положительно определены, А — неотрицателен. Дальнейшие рассуждения проведем в предположении, что исходная

задача является линейной, т. е. свойства материала и поток не зависят от искомых температур. Достаточным условием устойчивости семейства схем (7) является выполнение неравенства [3]

В, — 0,5 Ах(А + £>)>0.

Обозначив А? минимальное собственное число оператора В, можно записать

5, — 0,5Дх {А + О) > X?Е + Дх (а - 0,5) (А + О) >

г \в Г хв 1

^ ^-----Г. Дт; (а — 0,5) (А + П)> ------1-—-+Дт(а— 0,5) (Л + Я).

I. М + £>|] ~ К ’'Г \ 9 1И1 + |£>|| \ ’

Отсюда получаем условие устойчивости

М || +(!£>(!

+ Дх(а-0,5)>0. (8)

В работе [3] приведены оценки границ спектра операторов типа В и D при vf = const, — const; воспользовавшись ими, найдем

*?>», цоц<бр,

где v==minvt., {3 = тах(3г. Для оператора А аналогичные оценки

І г

даны в книге [4]: [| A(^ = тах^).

Условие (8) может быть переписано в виде

_1_+Дт(в-0,5)>0. (9)

6р + 4{*

Отсюда в случае явной схемы (о==0) при постоянных X, ср, 8, /, а получаем

Fo= — < Fo, = ,

с?Р * с?Р

где

с 1 *

Fo

* 12 1+ + ’

. 48Х

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

и 1 1

Ро * — -к-----

1 +

т

Асимптотически при /, Дх 0 приведенные условия устойчивости можно, с учетом [4], ослабить, а именно, при о = 0 получается

Ро *= для МКЭ, Ро # ж -у- для метода балансов.

В случае симметричной двухслойной схемы ^а==*|-------------схема

Кранка — Николсона [4]) и чисто неявной схемы (о= 1) неравенство (9) выполняется при любом шаге Дх, т. е. эти схемы безусловно

1 Л

устойчивы. Схема со значением о = — имеет при <р = 0,5(/ + /) или =/^л -+ Д^ порядок точности о (х2), порядок точности чисто

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

2. Рассмотрим задачу расчета нестационарных температурных полей в системе стержней, соединенных между собой в отдельных точках ~ узлах связи. Такая задача возникает при определении температур в сечениях тонкостенных авиационных конструкций [3]. Уравнение теплопроводности для каждого из стержней может быть записано в виде (1) с соответствующими граничными условиями на не соединенных с другими стержнями концах и дополнительными условиями сопряжения в узлах связи

т

2 5* = 0, = ... = Тт1,

к

где Пр т!—номера стержней, входящих в у'-й узел, Ткг -----------температура и тепловые потоки на концах х~0

и х-=^Ьк 6-го стержня.

Эта задача может быть решена МКЭ и методом подконструк-ций (суперэлементов). Применив, как описано выше, метод Бубнова — Галеркина для каждого стержня, получим системы обыкновенных дифференциальных уравнений вида (3) (число систем равно числу стержней), которые подчинены связям (10).

Система уравнений метода конечного элемента для всех стержней может быть записана в виде

В-^4-(А+0)Т =?, Г(0)=Г°,

где векторы Т и у составлены из векторов отдельных стержней Тк, глобальные матрицы В, А, О составлены из матриц отдельных стержней Вк, Ак, £)к. После замены дифференциальных операторов разностными рассмотрим, как и выше, семейство схем

+ +5)(зу + (1— 6)3;) = ?, У(0) = у°.

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

Исследование устойчивости проведем аналогично тому, как это было сделано для отдельного стержня. Достаточным условием является выполнение неравенства'(8), в котором следует заменить А, В и £) на А, В, £). Получим нужные оценки для границ спектра операторов, основываясь на соответствующих оценках для операторов отдельного стержня. Пусть Ц (Вк) — наименьшее собственное значение, оператора Вк и Х = шшХ|(5^), тогда

*

(у, •йу> = 20'*> у*)>Му. У)>».

А к

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

а® = гХ>0. Аналогичным образом оцениваются границы спектра операторов А и О. Норму оператора А можно оценить так:

(у, Ау) = £(у*, Ду*)<ЛЛ£(у*, .У*)<#ЛЛ(У, У),

£ &

где ЛА = тах|Л* |, R — наибольшее число стержней, встречающихся

к

в узле. Получаем ||Л»-</?ЛА, точно так же ||£)|| Достаточное

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

+Дх(а — 0,5) >0.

Л(Лл + Л°)

Отсюда следует, что двухслойные схемы — чисто неявная (а =1)

/ 1 \

и симметричная — — I — и в этом случае безусловно устойчивы.

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

3. Рассмотрим также задачу расчета нестационарных температурных полей в многослойной пластине толщины £, подвергающейся нестационарному аэродинамическому нагреву. Пусть подвод тепла к ией производится в течение времени х из внешней среды с параметрами Те(х), а(т), охлаждение осуществляется путем излучения с внешней поверхности и отвода тепла конвекцией на внутренней поверхности в хладагент с температурой Гхл и коэффициентом теплоотдачи ахл. Теплофизические свойства отдельных слоев являются известными функциями температуры. Температурное поле по толщине пластины определяется уравнением вида (1), в котором о=1, сс1 = а2 == 0, <2л = 0 с начальным Т(х, 0)~Т°{х) и граничными условиями вида (1), где

Я\ = ахл (Гхд — Т [д-—о)>

Я2 ==: ' ® {Т\х^Ь Те) ^о£(^4|х=^ 7оо),

(10)

с0 — постоянная Стефана — Больцмана, е —степень черноты,^ — температура среды на бесконечности.

Сформулированная таким образом задача имеет нелинейное граничное условие (10), кроме того, в уравнении (1) коэффициенты с и X разрывны и зависят от температуры. Используем для ее решения метод конечного элемента. Разделив каждый слой на отдельные элементы и применив, как описано выше, метод Бубнова —Галеркина для определения температур в узлах в момент (я+1)Дт, получим систему алгебраических уравнений

л л — е1у1 + а2у2 = -/„

л л- л

л, у/_!—е/у/ + а/+1.У*+1 = —/4, / = 2, ..., Ы— 1, (11)

aNУN_l-eN Ун = —/.V + - 71о),

ег = ех (о) = — 2v1 — сДх ({Л, + ахл);

е.=е£а) = — 2(у,_1 -|- V,) — оДх^.! + р.,), /==2, . . А' — 1;

(о) = — 2^_1 — аДх (р^_1 + а);

Я* = аь (а) = v^^1 — {Х._1 ах> / = 2, . . . , Л7;

Л = — ™хл Т’хл + ех (о — 1) ух —а2 (а — 1 ).у2;

/,= -а,(о- 1)у/_1+^(<з— 1)л — в(+1(»г-1)^|+1. г* = 2, ..., А/— 1;

\

f1) Ул?-1 “Ь (а ^).Улг

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

В качестве конкретного примера проведен расчет температурных полей в трехслойной пластине, внутренний и внешний слои которой имеют толщины 1 мм, а средний — 16 мм. Первый и третий слои выполнены из одного материала; теплофизические характеристики его и материала второго слоя при различных температурах приведены в табл. 1. В расчетах они определялись с помощью линейной интерполяции. Одна из поверхностей пластины предполагалась идеально теплоизолированной, на второй осуществляется конвективный теплообмен с потоком воздуха и отвод тепла излучением. Степень черноты считалась постоянной и равной 0,5. Коэффициент теплоотдачи изменялся линейно от значения 200 Вт/м2-град в начальный момент времени (х = 0 с) до значения 100 Вт/м3-град при х = 60 с и далее также линейно до 50 Вт/м2-град при х = 999 с.

Температура восстановления менялась аналогичным образом, причем ее соответствующие значения составляли 290 К, 860 К, 860 К. Начальная температура пластины считалась равной 290 К.

Таблица 1

т, к ^1. Хз, Вт/м2-град ^2» Вт/м2-град СРр ср2, Дж/м3-град ср2. Дж/м3-град

300 17,6 0,1580 3 569 600 36 638,37

400 19,7 0,2064 3 608 400 37 036,61

500 21,4 0,2739 4 035 200 41 417.29

600 | 23,1 0,3689 4 694 800 48 187,42

700 25,1 0,4988 5 393 200 55 355,80

800 26,6 0,6623 6 440 800 66 108,37

900 27,6 0,8640 6 790 000 69692,56

На рис, 2 показано полученное изменение температур внешней и внутренней поверхностей с течением времени, там же показано распределение температур по толщине при х = 100 с. Расчеты проводились как методом конечного элемента, так и методом тепловых балансов. При одинаковом количестве элементов и одинаковом шаге по времени быстродействие обеих программ на языке Фортран для ЭВМ БЭСМ-6 примерно одинаковое. За 2 мин 20 с машинного времени при разбиении пластины на 16 элементов

и шаге Дт = 0,1 с по обеим программам просчитываются температурные поля на временном интервале продолжительностью 810 с. Для сравнения результатов рассмотрим момент времени т=100с, когда различие в них наибольшее. В табл. 2 показаны значения температур на внутренней и

внешней границах в этот момент Таблица 3

в зависимости от шага Дч: при

решении задачи по неявной схеме N 7"внутр Т'внешн ■

г а о л и ц а г

(1, 1, 1) 349,50 674,52

X ^внутр Т’внешн (3, Ю, 3) 356,72 670,96

- МБ (3, 20, 3) 356,77 671,03

0,02 356,72 670,96 (6,'40, 6) 356,78 671,02

МБ 0,05 * 356,74 671,03 (9, 80, 9) 356,78 671,02

0,1 356,77 670,99

2 358,01 669,41 (1. 1. 1) 334,81 682,31

• мкч (3, Ю, 3) 354,49 672,11

- 0,02 354,49 11\ С/ 672,11 с *—• 1 (3, 20, 3) 355,76 671,57

мкэ 0,05 354,51 672,09 (6, 40, 6) 356,40 671,30

0 — 1 0,1 354,54 672,05 (9, 80, 9) 356,72 671,16

2 355,78 670,54

(1. 1. О 334,80 682,33

. 0,02 354,50 672’13 МКЭ (3, 10, 3) 354,50 672,13

МКЭ 0,05 354,52 1»А 672'14 пк о =0,5 (3, 20, 3) 355,77 671,59

ю сГ II о 0,1 354,56 672,15 (6, 40, 6) 356,40 671,32

2 356,09 672,47 (9, 80, 9) 356,72 671,18

метода баланса и методом конечного элемента по чисто неявной схеме (о = 1) и по симметричной схеме (а = 0,5). Во всех трех случаях принималось следующее разбиение слоев пластины и а элементы: А/~(3, 10, 3). Видно, что температуры, найденные методом конечного элемента, отличаются от подсчитанных методом баланса примерно на 2,2°С на внутренней поверхности и на 1,1°С — на внешней. При значениях Дт = 0,02с и 0,1 с результаты практически не отличаются.

Были проведены расчеты с различным разбиением пластины по толщине при Дт = 0,02с. Результаты даны в табл. 3. Видно, что сгущение сетки при использовании метода балансов не вносит существенных изменений. Сгущение сетки в расчетах методом конечного элемента приводит к уточнению, причем уточненные значения приближаются к результатам, полученным методом баланса. Таким образом, в данной задаче для достижения методом конечного элемента точности, соответствующей методу балансов, требуется использовать сетку, значительно более мелкую.

ЛИТЕРАТУРА

1. Иванов С. Н. Экономичный численный метод расчета температурных полей в тонкостенных цилиндрических конструкциях. „Инженерно-физический журнал*, т. XXVIII, № 1, 1975.

2. 3 а м у л а Г. Н., Иванов С. Н. Экономичный метод расчета температурных полей в тонкостенных авиационных конструкциях. „Ученые записки ЦАГИ", т. VII, Кг 3, 1976.

3. Самарский А. А. Теория разностных схем. М., „Наука*,

1977.

4. С т р е н г Г., Д ж. Ф и к с. Теория метода конечного элемента. М., „Мир*, 1977.

Рукопись поступила 171VI 1980 г.