УЧЕНЫЕ ЗАПИСКИ Ц А Г И Т о м VII 197 6
№3
УДК 536.2
ЭКОНОМИЧНЫЙ МЕТОД РАСЧЕТА НЕСТАЦИОНАРНЫХ ТЕМПЕРАТУРНЫХ ПОЛЕЙ В ТОНКОСТЕННЫХ АВИАЦИОННЫХ конструкциях
Г. Н. Замула, С. Н. Иванов
Изложен экономичный численный метод расчета нестационарных температурных полей в плоских сечениях тонкостенных авиационных конструкций, построенный па основе неявной схемы метода балансов и на методе прогонки.
При расчетах нестационарных температурных полей в конструкциях летательных аппаратов наибольшее распространение получил метод элементарных тепловых балансов [1] и его модификации [2, 3], идейно близкие интегро-интерполяционному методу (методу баланса) построения однородных разностных схем для уравнения теплопроводности с переменными коэффициентами [4]. При этом время делят на отрезки Дт, а конструкцию разбивают на элементы, для каждого из которых записывают уравнение теплового баланса при определенных предположениях относительно тепловых потоков между элементами. Последние либо вычисляют по распределению температур по элементам на предыдущем шаге расчета, вследствие чего приходят к явным формулам [1—3] для подсчета температур, либо выражают неявно через искомые температуры элементов (неявная схема расчета [4]). Использование явной схемы при расчете реальных конструкций часто бывает затруднено и малоэффективно в связи с требованием устойчивости счета, приводящим к введению излишне мелкого шага по времени Дт и к значительной затрате машинного времени. Этого недостатка лишены безусловно устойчивые неявные схемы, однако в этом случае на каждом шаге расчета приходится решать систему алгебраических уравнений определенной структуры, соответствующей типу конструкции. В одномерном случае, когда эта система имеет диагональный вид, устойчивым и экономичным (т. е. требующим для перехода от шага к шагу число арифметических действий и затраты машинного времени того же порядка, что и по явной схеме) методом ее решения является метод прогонки [4]. При определении температурных по-
лей в плоских сечениях тонкостенных конструкций приходим к совокупности одномерных задач, связанных между собой в отдельных точках. В настоящей статье для решения получаемых при этом алгебраических систем предлагается общий экономичный метод, основанный на методе прогонки. Рассмотрены два типа условий нагревания конструкции: соответствующие конвективному нагреванию в полете и лучистому — при воспроизведении его в лабораторных испытаниях. Кроме того, дан способ учета внутреннего лучистого теплообмена, зависимости теплофизических свойств от температуры и термических сопротивлений в соединениях.
Рассмотрим нагревание тонкостенного цилиндрического тела. Тепловыми потоками в направлении образующих и температурными перепадами по толщине элементов этого тела можно пренебречь. Задача расчета нестационарных температурных полей в этом случае сводится к задаче определения температур в системе тепловых стержней переменной толщины, соединенных между собой в отдельных точках (узлах, см. фигуру). Обычно будем пренебрегать тепло-
Узлы
выми потоками на концах стержней, не входящих в узлы. Пусть система содержит т стержней, соединенных между собой в п узлах.
Введем разбиение каждого из стержней на ^ ^ = ,/и| эле-
ментов так, чтобы все разрывы характеристик, внешних воздействий, угловые и узловые точки совпадали с границами элементов.
Уравнение теплового баланса для типового г-го элемента по чисто неявной схеме имеет вид
'Т,Р + 1 _ 'ГР
I 8 С Р ----!_ = 0$' + + 0??\ (1)
Дт
Здесь <Зкг — тепловой поток от омывающих две боковые стороны элемента сред с температурами Тс1, Тс1 и коэффициентами теплоотдачи а1(- а2;,
<2к/ = [«11(ТЫ- Т,) + а2г (Гс, - 7^)] /г;
Ял! — лучистый поток,
Qл г == (*7реэ. г ~|~ <7рез. г) 1ц (2)
где <7Рез.1, ^Рез.< — плотности результирующего (на двух сторонах элемента) излучения с учетом лучистого теплообмена между элементами;
(Зт г — тепловой поток вследствие теплопроводности от соседних элементов I— 1, /-ЬЬ
Qт I — а1 {Т1-1 — Тг) -|- а,-+1 (Т!+1 — Т,), (3>
где аг, а,+1 — тепловые проводимости между центрами соседних элементов, определяемые с учетом граничных условий;
7? — температура в момент времени иДт;
— длина;
8г — толщина;
— удельная теплоемкость;
РI — плотность г-го элемента.
Внешние тепловые воздействия на конструкцию могут соответствовать как случаю полета с заданным изменением по траектории а и Тс, так и случаю лабораторных испытаний с заданными в каждый момент времени температурами некоторых из ее точек Р, по которым ведется управление лучистым нагревом. Ограничимся для простоты изложения одной точкой Р, при этом на часть стержней системы падает поток излучения плотностью <7, подлежащей определению. Точку Р будем считать узлом.
Запишем уравнения теплового баланса в местах стыка стержней (в узлах):
тк
£д<*> = 0, к = \, . . ., п. (4)
5=1
В этом выражении 0^ — тепловой поток вследствие теплопроводности, поступающий по х-му стержню в к-й узел, в котором сходится тк стержней. В случае, если система состоит из одного замкнутого стержня, узлом назовем точку стыка первого и последнего элементов.
Вначале будем пренебрегать собственным излучением конструкции и зависимостью теплофизических свойств материалов от температуры. Для определения температур элементов и узлов на каждом (л+1)-ом шаге расчета по времени возникает необходимость решения системы М + я линейных алгебраических уравнений (1), (4), которую можно рассматривать как т трехдиагональных подсистем, связанных между собой через неизвестные температуры узлов 0(А) (&=1, - . и). Для каждого стержня подсистему пред-
ставим в виде
Щ Г,-, — е1 Т1 -)- а<+1 7\+1 + ег /г <7 = —/,-, г=1, . . ., п (5)
с краевыми условиями
Т0 — Т1 7^+1 = Тг Т’лгЧ- 02- (6)
Здесь Т0, Тц+\ — температуры концов стержня; гг—поглощательная способность г-го элемента;
а,
2Аг 5; 2ХД-—1 В/—1
. __Ь_________II-1
2X; 2X1 —! Вг—1
// /г-1
ЙЛЧ-1 =
2Хдг бд^
‘я
= (а1г + ?гЛ^/ + ^ + а/+1 +
Ь в| Р/ . Дт ’
А = ^_£т тп + ^ {ли Тс1 + ^ Тс ) _ /г (?рез . + ^рез /)}
где Х/ — теплопроводность г-го элемента.
Если тг, 0, (/=1, 2) известны, система (5), (6) может быть решена методом прогонки:
Т1 — В1+1 7г + 1 + /^+1 + '■'/4-1 + Сг+1 Я, I — А/, . . ., 1.
(7)
Здесь 5/+1, ^1+1, '*1+1, чг+1 — прогоночные коэффициенты, подсчитываемые по рекуррентным соотношениям:
1
в\ = чй Л = °; ^ = 1; ^ = 0;
В1+1 = (е, — а, в,)-1 <*1+15 /^+1 = (е,— аг В,)~] (/< + а1 У7,);
»/+! = (г, — Щ В,)-1 V, а,; С|+1 = (е?г — аг В;)-1 (е, /г + а1 С;), *=1, . . М
(8)
Температура ТУ^ определяется с помощью второго краевого условия (6):
7^+1 = (1 — Тг Влг+1)-1[тГ2 (^лг+1 + %+10! + Слг+1 С[) + 02]- (9)
Возможен также способ решения, использующий прогонку в противоположном направлении:
Т1 — Ах Т 1-\ Ф( + ц* ®2 “Ь ^+1^1 I — 1, - ■ ч
(10)
где
■<4лч-1 = 72’. ®лч-1 *= 0; ^+1 = 1; ^лч-1 — 0;
А1 = (е1 — а,+1А,+1)-1а1] Ф, = (е1 — а1+х Л,*,)-1 (А + л»+1 Ф/+1); \
Нч = (£;--+ + 1 Н-/ + Г> ^ — (&1+\ а‘ + 1 •^<+0-1(£гЛ + Я/ + 1&г + 1),
/ = Л/, . . 1;
То —(1 — ТИ1)-1 [71 (Ф, + ^ е2 + «! д) + 0,]-
(П>
(12)
Для стержней, на которые лучистый поток <7 не падает, и в случае полетных условий (<7 = 0) коэффициенты С,-, 9г в соотношениях (7) — (12) тождественно равны нулю. Запишем формулы для тепловых потоков на концах стержня:
Qo — а\(Т1 Т0) — ■— а1 ------1- 0! + а1 -----^ 02 +
1 — ьА 1 — тп А
+ Ь,д + а, Ф1?
1 - 71 А
1 — 71
1 -в
N+1
0,
(?ЛН-1 = ЙЛЧ-1 (ТN + 7’л/+1 ) = — <^N+1 ------------------------------ ^2 *г
1 - 72 -о лг+1
+ аы+1 -—-—^— ^N+1 ©1 -г Ллг+1
1 —72
1 — 72^ + 1
+ aN +1
1 —72
1 — 72 в,
1 ~ 72 Вдг+1 р Ы+1 •
— Сл?+ г д +
ЛГ+1
Выражения (13) линейны относительно температур концов 0ц
02 и плотности теплового потока д. Воспользуемся этим свойством для построения метода расчета температур элементов и узлов в системе тепловых стержней. Проведем прямую прогонку по (8) или (11) для тех стержней, у которых один конец свободен, в направлении от свободного конца (где 7=1, 0 = 0) к узлу (где 7 = 0). Далее проведем прямую прогонку по остальным стержням слева направо согласно (8) и справа налево согласно (11), считая 71 = 72 = 0. Выразим в виде (13) тепловые потоки на концах стержней, входящих в узлы, через температуры концов 0^ 02 и величину д и подставим их в уравнения баланса (4). В результате этого, с учетом равенства температур 0Ь 02 температурам соответствующих узлов, получаем систему п линейных уравнений относительно п температур узлов 0(6), если условия нагревания соответствуют полетным. В случае лабораторных испытаний температура одной из узловых точек известна, но неизвестна плотность теплового потока д и температуры п—1 остальных узлов. После решения полученной системы уравнений температуры всех элементов определяются обратными прогонками по формулам (7) или (10).
Для системы стержней, сходящихся в одной точке, и полетных условий нагрева после подстановки потоков (13) в условие баланса <4) имеем одно уравнение с неизвестной температурой узла 0(1), после определения которой найдем температуры элементов обратной прогонкой. В этом простейшем случае, если в узле сходятся лишь два стержня, рассмотренный алгоритм сводится к встречной прогонке [4].
Предложенный метод может быть распространен на случай, когда в системе содержатся двухслойные тепловые стержни, представляющие собой два обычных стержня, соединенных по боковым сторонам так, что на поверхности контакта имеется термическое сопротивление. Такой расчетной схемой пользуются при описании реальных соединений тонкостенных элементов, например, обшивки с подкрепляющими элементами и т. п. Будем считать, что концы двухслойных стержней принадлежат двойным узлам, имеющим две точки с разными температурами, и для них записываются два уравнения теплового баланса типа (4).
Проведем разбиение каждого двухслойного стержня на элементы так, чтобы границы соприкасающихся элементов совпадали, и выпишем для них расчетные уравнения в векторном виде, обобщающем (5) и (6):
аг7\_1—е1 7^ а-1 + 1 Т^х + <78г^г =—/;> * —1. ■ • •> №\ (14)
Т0=Ъ + 01. 7^+1== 7з 7^+02. (15)
Здесь индексы „1“ и „2“ относят соответствующую величину к первому или второму из двух стержней; /?г — величина термического сопротивления; ПЛОТНОСТИ тепловых ПОТОКОВ <7,, могут быть равны либо <7, либо нулю. Решение уравнений (14) с граничными условиями (15) может быть записано формулами (8) — (12), если считать, что в последних Т;, fi, /,., 0., С,., Fi, Ф,— векторы, а а1> еи Я> Ть Н7. #г> ~ квадратные матрицы второго порядка.
Л
Единица заменяется единичной матрицей Е, у которой диагональные элементы равны единице, а остальные —нулю. Тепловые потоки на концах такого стержня определяются выражениями, линейными
относительно 8„ 02» Я'-
-*■ Л ( Л Л А ДА ГА Л Л -*■ I
Qo — — а11(^ — — "Ь А)~х 1.71 (®1 1^1 02 + Я ®1) + 0)] —
-*■ Л -*■ А -► )
— (®1 + н-1 02 -Ь я ;
"*■ II Л Л \/Л Л Л \-1 ГЛ /-> Л -*• Л \ ]
$N+1 = алг+т |\Ду+1 —Е) \Я — 72 Вы+\) [тг'.^л+х + ^лг-н 01 + <7Слч-1) + 02] +
: ■' + {^N+1 + >Л7+1 01 + 9СЛ/+1)} ,
где (Зо
Ою
^20
QN +1 =
0.1, л+1
0?, N+1
Процедура определения температур в узлах и элементах системы, содержащей двухслойные стержни, отличается от описанной выше лишь тем, что для таких стержней проводится матричная прогонка.
Отметим еще два важных обобщения изложенного метода.
; 1. Приведенные выше формулы были записаны для чисто неяв-
ной схемы. Очевидно, метод применим для любой неявной схемы, включая симметричную двухслойную схему, при этом изменятся лишь выражения для коэффициентов и правых частей в соотношениях (5) и (14).
2. Формулы (13) связывают тепловые потоки с температурами на концах стержней, поэтому после подстановки (13) в соотношения баланса (4) для тепловых потоков приходим к разрешающим уравнениям относительно температур в узлах. Такой подход имеет определенное сходство с методом перемещений в строительной механике стержневых систем. Очевидно, могут быть расссмотрены и другие способы расчленения стержневой системы на подсистемы и подходы, аналогичные методу сил и смешанному методу, использующие в качестве основных неизвестных наряду с температурами и тепловые потоки в местах сечений.
Рассмотрим теперь вопросы, связанные с учетом лучистого теплообмена, и зависимости теплофизических свойств материалов от температуры при использовании изложенного метода расчета температурных полей в тонкостенных авиационных конструкциях. В этом случае в соотношениях (2)
*7рез / = 811 (со Яи)< Ярез1 — £2( (^о Ян)* (17)
где с0 — постоянная Стефана — Больцмана; е1 о 41 — степени черноты двух боковых сторон элемента;
ЯЯн — плотности падающего на них излучения, определяемые для каждой невогнутой полости конструкции системой линейных алгебраических уравнений вида [3]
N1 \ т_т / //.. \
= 2 [£^»^ + (1-^)^]“7г + ^с(1_2“гг]’/=1- • • -м'-(18)
Здесь г, /' — новые номера всех элементов числом Ыи участвующих в лучистом теплообмене в данной полости;
<7С — плотность внешнего излучения, падающего на замыкающий полость отрезок;
— взаимная поверхность излучения элементов с номерами г и /, подсчитываемая по методу натянутых нитей:
_ (^* / “Ь /—1 ■^'1, /—1 —1, /)» ^ -ф- /, И ^ I с — 1 *
= V (*< - *,)2 + (у, - уу,
где х, у — координаты конечных точек элементов в выбранном направлении обхода; хоу У о — координаты начальной точки первого элемента.
При заданных температурах элементов система (18) достаточно просто и быстро решается методом итераций, например вида
Я\к) = 2 ^0 Т' - ^ + (1 -8у) Я?-"] + <7с,
У= 1 '
1 = 1, 2, . . М, к, = 1, 2, . . .
до выполнения условия шах
ч\к)
< £0 при <7<°>=0 или
д(0)П+1 — дп> сходящимся ввиду выполнения достаточного условия [5]
м■ я-
2(1-.,) ^<1.
/=1 ‘
Сами температуры в (17) и (18), как и зависящие от них теплофизические характеристики с, X, е в (1), (3), (17), (18), могут задаваться
как по предыдущему, так и по последующему (с итерациями) шагу расчета, либо экстраполяцией по предыдущим шагам, т. е. различными рассмотренными для квазилинейных уравнений теплопроводности в [4] способами. Безусловная устойчивость схемы при этом сохраняется.
Расчеты температурных полей в авиационных конструкциях с использованием предложенного метода и созданных на его основе
программ для ЭВМ показали его высокую эффективность и экономичность. В частности, применение его сокращает примерно на порядок время счета на ЭВМ по сравнению с временем счета по программам, составленным на основе явной схемы, при достаточной для практики точности.
ЛИТЕРАТУРА
1. Ваничев А. П. Приближенный метод решения задач теплопроводности при переменных константах. „Изв. АН СССР. ОТН“, 1946,
№ 12.
2. Марченко В. М. Температурные поля и напряжения в конструкции летательных аппаратов. М., Машиностроение, 1965.
3. 3 а м у л а Г. Н., Марченко В. М. Современное состояние теории лучистого теплообмена и ее приложения к расчету температурных полей в некоторых инженерных конструкциях. Сб. „Тепловые напряжения в элементах конструкций”, вып. 5. Киев, „Наукова думка', 1965.
4. С а м а р с к и й А. А. Введение в теорию разностных схем. М., „Наука”, 1971.
5. Березин И. С., Жидков Н. П. Методы вычислений. Т. 2,
М., Физматгиз, 1962.
Рукопись поступила 10/IV 1975 г.