CC BY
33
АВИАЦИОННАЯ И РАКЕТНО-КОСМИЧЕСКАЯ ТЕХНИКА
УДК 629.7.042.2.001.24:622.998
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ТЕПЛОВОГО СОСТОЯНИЯ
ГЕРМЕТИЗИРОВАННОГО ОТСЕКА ГИПЕРЗВУКОВОГО САМОЛЁТА
© 2018 С. А. Гусев12, В. Н. Николаев3
1 Институт вычислительной математики и математической геофизики Сибирского отделения Российской академии наук, г. Новосибирск
2 Новосибирский государственный технический университет 3 ФГУП «Сибирский научно-исследовательский институт авиации имени C.A. Чаплыгина», г. Новосибирск
Статья поступила в редакцию 26.06.2018
Разработан метод определения теплового состояния герметизированного отсека гиперзвукового самолёта, основанный на использовании математической модели теплового состояния отсека с металлическими сотовыми конструкциями. Математическая модель системы герметизированного теплоизолированного отсека с системой охлаждения обшивки представлена системой дифференциальных и алгебраических уравнений, описывающих конвективный и радиационный теплообмен наружной поверхности обшивки, конвективный теплообмен в канале системы охлаждения и теплообмен теплопроводностью сотовой конструкции. Проведена разработка стохастического метода решения прямой задачи теплообмена металлической сотовой конструкции с использованием стохастических дифференциальных уравнений методом Эйлера и методом случайного блуждания по сферам. Моделирование траекторий случайного процесса осуществлялось с использованием параллельного датчика гауссовских случайных величин из библиотеки Intel MKL. При решении жёсткой системы обыкновенных дифференциальных уравнений используется неявный метод Ро-зенброка второго порядка. Обратная задача для металлической сотовой конструкции решена по стохастическому квазиградиентному алгоритму с переменной метрикой. Кроме этого параметрическую идентификацию модели теплового состояния отсеков предложено прводить композицией метода наискорейшего спуска, метода Ньютона и квазиньютоновского метода Бройдена-Флетчера-Гольдфарба-Шэнно. Вектор произведения погрешностей оценок коэффициентов нелинейной математической модели определён по выражению, являющемся функцией квантиль %2 - распределения и матрицы Грамма. При этом используется метод проецирования совместной доверительной области оценок на координатные оси пространства коэффициентов. Исследования гиперзвукового самолёта проводились с использованием Норм лётной годности. Получены потребные значения толщины металлической сотовой конструкции герметизированного отсека, температуры и расхода хладагента системы охлаждения обшивки гиперзвукового самолёта. Оптимизация характеристик металлической сотовой конструкции герметизированного отсека гиперзвукового самолёт и системы охлаждения обшивки и проводилась для обеспечения надёжности работы самолёта, его бортового оборудования, комфортных условий в кабине экипажа и салоне пассажиров. Ключевые слова: Гиперзвуковой самолёт, сотовая конструкция, система охлаждения обшивки, математическая модель.
Работа С.А. Гусева выполнена в рамках государственного задания ИВМиМГ СО РАН
(проект 0315-2016-0002).
ВВЕДЕНИЕ
Оптимизация характеристик металлической сотовой конструкции герметизированного отсека гиперзвукового самолёт и системы охлаждения обшивки и проводилась для обеспече-
Гусев Сергей Анатольевич, доктор физико-математических наук, старший научный сотрудник Института вычислительной математики и математической геофизики СО РАН, профессор Новосибирского государственного технического университета E-mail: sag@osmf.sscc.ru Николаев Владимир Николаевич, доктор технических наук, начальник отдела. E-mail: nikvla50@mail.ru
ния надежности работы самолёта, его бортового оборудования, комфортных условий в кабине экипажа и салоне пассажиров.
Параметры режима полёта и воздушной среды за бортом модели применения самолёта приведены на рисунке 1.
Панели фюзеляжа самолёта (рисунок 2) состоят из титановой обшивки и дюралюминиевых сотовых конструкций. Титановая обшивка самолёта и дюралюминиевые сотовые конструкции разделены друг от друга каналами. В качестве активной теплозащиты используется жидкий хладагент.
290 -270 1 250 230 -210 -
0 3600 7200 10800 14400 18000 с Рис. 1. Параметры режима полёта и воздушной среды за бортом модели применения гиперзвукового самолета: Р ¡гШ ~ давление воздуха за бортом за пределами теплового пограничного слоя; М - число Маха на внешней границе пограничного слоя;
Т.
■ температура воздуха за бортом за пределами теплового пограничного слоя
ру, V,
а)г,оиЬ ' ай-/>1Л
Внешняя обшивка
Труба
Хладаген
Сотовая теплоизоляция
Внутренняя обшивка
Герметизированный отсек
Рис. 2. Схема расположения элементов сотовой теплоизоляции и активнойтеплозащиты гиперзвукового самолёта: - воздушная скорость полёта; - плотность воздуха за бортом; 7\ - температура воздуха за бортом за пределами теплового пограничного слоя; рс1 - плотность хладагента; ]¥с1 - скорость хладагента; Г, - температура хладагента; ро!> - плотность воздуха в отсеке; — скорость воздуха в отсеке; Т( - температура воздуха в отсеке
Критерии оптимизации теплофизических характеристик герметизированного отсека определены в соответствии с Нормами летной годности АП 25 [1].
В качестве математической модели теплового состояния обшивки, дюралюминиевых сотовых конструкций и системы охлаждения обшивки для герметизированного отсека примем систему уравнений, описывающих конвективный и радиационный теплообмен наружной поверхности обшивки, конвективный теплообмен в канале системы охлаждения и теплообмен теплопроводностью сотовой конструкции.
1. ВЫНУЖДЕННАЯ КОНВЕКЦИЯ СНАРУЖИ ГЕРМЕТИЗИРОВАННОГО ОТСЕКА
Температура воздуха за бортом у поверхности обшивки Г определяется по формуле
Т* =Т + 0 72 (Т -Т )
1 спг.оШ ^ ' ^У1 е 1тг,оШ>>
(1)
где Та1гоМ — температура воздуха за бортом за пределами теплового пограничного слоя.
Температура восстановления Т в уравнении (1) определяется по выражению
т„ = /:
(1+М г(к-1) / 2),
(2)
где г - коэффициент восстановления температуры (для ламинарного пограничного слоя
ОМ
г = д/Рг0((Г , для турбулентного - г _М - чисто Маха на внешней границе пограничного слоя; к=с[/су - отношение удельной теплоёмкости газа при постоянном давлении и объёме.
Коэффициент теплоотдачи асуоцГ наружной поверхности обшивки будем вычислять по формуле [2]
(Хщтй A26.RC.,,) ' - (PC)-0'666 ¿4 й- JГair out (3) для ламинарного пограничного слоя при критерии Рейнольдса ReOH¡ < 1.106 и по формуле [2]
=0,181(lg Relj-?'584 (РС,Г2/3 pout ср Vairout (4)
для турбулентного пограничного слоя при критерии Рейнольдса ЫО6 <Re0MÍ< ЫО9.
В уравнениях (3), (4) использованы следую-
¡j. ¡ц
щие обозначения: ReOIrf, Prouí — критерии соответственно Рейнольдса и Прандтля при темпе-
* *
ратуре воздуха Т ; Pout — плотность воздуха
*
за бортом при температуре Т; ср — удельная теплоёмкость воздуха при температуре Т; V — воздушная скорость полёта.
Критерий Reolrt вычисляется по формуле
air.out
R
air, out Р out cv,out' № air > (5)
где 1суоМ — характерный размер для наружной
*
поверхности обшивки; //ш>. — динамическая вязкость воздуха.
Характерный размер £суон, принимается равным расстоянию от носовой точки фюзеляжа до середины рассчитываемого отсека или его части.
*
Динамическая вязкость // определяется по выражению
//*,,. =0,1222229-10 6+ 0,682674-10-8 Г* -
0,313155-10 11 {Т*)\ (6)
Критерий Рг';, вычисляется по формуле
Р' «« = Ср Айг I ^СПГ' О)
где X' — теплопроводность воздуха.
1*
Теплопроводность Ла!Г определяется по выражению
£т = 0,141483-Ю 2 + 0,896161-104 Г* -
-0,204759-10-7 (Г*)2. (8)
Плотность воздуха р ш за бортом вычисляется по выражению
Р*аа = 3,4826 • 10~3 /Т*, (9)
где рфоМ — давление воздуха за бортом.
Тепловой поток ОаоМ , поступающий к наружной поверхности обшивки от прямого и отраженного излучения Солнца и прямого излучения Земли определяется по уравнению
Оа%Ш ~ С0 &с\',оШ ^су,оШ ^ср,01й (*")» (Ю)
где с0 - постоянная Стефана-Больцмана; -коэффициент черноты излучения наружной поверхности обшивки; площадь обшивки при наружном конвективном теплообмене; -
температура наружного слоя обшивки.
2. ТЕПЛООБМЕН В КАНАЛЕ СИСТЕМЫ ОХЛАЖДЕНИЯ
Теплообмен в канале системы охлаждения будет определяться теплообменом наружной поверхности обшивки, теплообменом сотовой конструкции и обыкновенным дифференциальным уравнением, описывающим конвективный теплообмен внутренней поверхности каната и перенос теплоты из одной части каната в другую
^Ы,к,1 = Щщт О ^сг.оикМ I ^с1.к (Т к (4 ' О - 1 с],к ) +
+ аЬс,с1,к(Г)Ркс,с1,к ¡СС! к {Т,К с1 д-(4,/) ~Тс1к) (11) + Ср,с1^с1,с ^с ¡Сс1,к 1
где Тс1к1, Тс1к- температуры хладагента соответственно в (к— 1 )-ой и к-ой частях каната; ТНсс1к -температура поверхности сотовой конструкции на границе с к-ойчастью канала; ак- коэффициент теплоотдачи внутренней поверхности /<-ой части канала; й - площадь теплоотдающей поверхности к-ой части каната; ]а с - массовая скорость хладагента в канале; с с1 - удельная теплоёмкость хладагента; /•' - площадь поверхности сечения канала; Сс1к - теплоёмкость хладагента в к-ой части каната.
Тс1к в выражении (11) с индексом t означает её дифференцирование по времени г:.
Теплоёмкость хладагента Сс1к определяется по выражению
Л ~ ср Рал
где рсгк — плотность хладагента в к-ой части каната; Шс1с — скорость хладагента в канале; At — интервал дискретизации времени при решении системы дифференциальных уравнений; У (к — объём хладагента в к-ой части канала.
Коэффициент теплоотдачи ак внутренней поверхности канала определяется уравнениями, приведёнными в [3].
Ламинарное течение при критерии Рейнольдса Иег( < 2100 описывается уравнением
Nit
cl,d,av
= 0,15 Re
°'33 Рг °'33х
cl.d 1 'el л
x (Grdd Prcl )ол (Prcí / Prk .
(13)
Критерий Нуссельта ,\:и , А т представляет собой безразмерный средний коэффициент теплоотдачи и определяется из выражения
М'сЫт = ас1М^1 Кп (14)
где о-с1Лт - средний коэффициент конвективной теплоотдачи; с1 - эквивалентный диаметр канала, равный учетверённой площади поперечного сечения канала, делённой на его полный «смоченный» периметр [4]; Хс{ - теплопроводность хладагента.
Критерий Рейнольдса а характеризует режим движения хладагента и определяется из выражения
где \\' ау - средняя скорость хладагента; рсг -плотность хладагента; - динамическая вязкость хладагента.
Критерий Прандтля Рг.( определяет физические свойства хладагента при определяющей средней температуре хладагента и вычисляется из выражения
Р гс/ =
Ср,с1 (16)
Критерий Прандтля Рг(. определяет физические свойства хладагента при определяющей средней температуре поверхности канала.
Критерий Грасгофа Сг^ характеризует подъёмную силу, возникшую вследствие разности плотности хладагента, и определяется из выражения
Gr
cl M
■ flclgd3ATd/v:
(17)
где Рс( — коэффициент объёмного расширения хладагента; & — ускорение в поле тяготения; АТс[ — разность между температурой поверхности каната и хладагента за пределами пограничного слоя; V — кинематическая вязкость хладагента.
Коэффициент е; учитывает изменение среднего коэффициента теплоотдачи по длине канала и зависит от отношения длины канала к его эквивалентному диаметру Щ. Значения этого коэффициента представлены в таблице 1.
Переходной решим течения при критерии Рей-нольдса 2100 < < МО' описывается уравнением
=К0 РГс/МЗ(Ргсг/Рг,)0-25^, (18) Таблица 1. Значения е.
где значения коэффициента К0 зависят от критерия Рейнольдса и представлены в таблице 2.
Турбулентное течение при критерии Рейнольдса Р.ег( > МО* описывается уравнением
= 0,021 Ш4/'щ Рг/43 (19)
Расход хладагента Сс] предлагается представить линейной регрессией вида
^(0 = %0+%1А-(0М2(0, (20)
где За 0, 9(; | - коэффициенты модели, которые будут оцениваться.
3. ТЕПЛООБМЕН ТЕПЛОПРОВОДНОСТЬЮ СОТОВОЙ КОНСТРУКЦИИ
Систему уравнений теплообмена сотовой конструкции представим в виде одномерных уравнений теплопроводности, описывающих процесс передачи тепла в многослойной сотовой конструкции [2, 5]
€ы{х)Тш = (4(1)Гк1), 0 < х </, 0 < Г <&; (21)
/Сш(Ш4)~Тм),Х=0; (22) (23)
Т1к(0,х) = Т0(х), 0 <х<1, (24)
х е dur ;
Сф, хе air ,
Âdur, х е dur ; Яаг, хе air, при ламинарном течении
l/d 1 2 5 10 15 20 30 40 50
El 1,9 1,7 1,44 1,28 1,18 1,13 1,05 1,02 1,0
Таблица 2. Значения коэффициента К0 от критерия Рейнольдса Re
Rec/xl03 2,1 2,2 2,3 2,4 2,5 3 4 5 6 8 10
к0 1,9 2,2 3,3 3,8 4,4 6,0 10,3 15,5 19,5 27,0 33,3
Таблица 3. Значения е; при турбулентном течении
l/d Rec/ = 2-103 Rec/ = 2-104 Rec/ = 2-105
1 1,9 1,51 1,28
2 1,70 1,40 1,22
5 1,44 1,27 1,15
10 1,28 1,18 1,10
15 1,18 1,13 1,08
20 1,13 1,11 1,06
30 1,05 1,05 1,03
40 1,02 1,02 1,02
50 1,00 1,00 1,00
то есть коэффициенты С)к, Х)к зависят от того, в каком слое рассматривается перенос тепла.
В уравнениях (21) - (24) использованы следующие обозначения:
%А х) - температура сотовой конструкции; Т. г - первая производная Т1к по V, ТНсх - первая производная Г, по х ; Щ1СХХ - вторая производная ТНс по х ; С/1с(х) - объёмная теплоёмкость сотовой конструкции обшивки, определяемая теплоёмкостью дюралюминия С и теплоёмко-
1 сотро
стью воздуха С ,„; /., (/) - теплопроводность сотовой конструкции, определяемая теплопроводностью дюралюминия %Л1г и теплопроводностью воздуха\.г; а1кк - коэффициент теплоотдачи на-ружной поверхности сотовой конструкции; ас. - коэффициент теплоотдачи внутренней поверхности сотовой конструкции; Р,х к - площадь сотовой конструкции при наружном и внутреннем теплообмене; 7\. - температура воздушной среды в отсеке или в части отсека; I - толщина сотовой конструкции.
4. ВЫНУЖДЕННАЯ КОНВЕКЦИЯ В ГЕРМЕТИЗИРОВАННОМ ОТСЕКЕ
Коэффициент теплоотдачи ® внутренней поверхности сотовой конструкции можно вычислить по формуле [5]
а
he. in
= 0,66 /L, Re°i5Pr,°;43/L
he,in
для ламинарного пограничного слоя при критерии Рейнольдса Ке.я < 4-Ю4 и по формуле [6]:
(26)
аШп= 0,037 4rRe«Pr^3/L
® ~ [Lhc' Tci, j,
%1 (28)
и включает в себя необходимые характеристики (значения толщины сотовой конструкции 1(, вм, температуры хладагента Тс1 в К и коэффициенты
0, 3( ; 1 , описывающие расход хладагента Сс] в кг/с).
Для решения прямой задачи теплового состояния сотовой конструкции фюзеляжа рассмотрим математическую модель теплообмена в гетерогенной структуре, как параболическую краевую задачу с разрывными коэффициентами. Описание и физические свойства гетерогенных структур можно найти в [7].
При этом делалось сглаживание разрывных коэффициентов уравнения теплопроводности на основе интегрального усреднения, а численное решение стохастических дифференциальных уравнений осуществлялось методом Эйлера и методом случайного блуждания по сферам [8].
Предполагается, что рассматриваемый процесс теплопередачи происходит на отрезке времени [0, г^] и описывается следующей краевой задачей для уравнения теплопроводности
дТг
he
dt
i
^Эх
Ъ{х)
дТ,
he
i V
Эх
i J
(25) где b(x) =
bA,,xe Gx Pair ' X G G2 ' rhc(xfi) = (p(x).
hc,in
при критерии Рейнольдса Re. > 4-10я.
Теплопроводность % и динамическая вязкость ц в формулах (25), (26) в критериях Re, Рг определяются соответственно по выражениям (8), (6) для температуры воздуха в отсеке Т\ ..
Произведение Jw плотности р и скорости Wair воздуха в отсеке или массовая скорость Jw в критерии
in = Р air ^air ^ in ^ № air (27)
принимается из результатов эксперимента.
Характерный размер /. , принимается по аналогии с i.( ш равным расстоянию от начала отсека до середины рассчитываемой части отсека.
5. ОПТИМИЗАЦИЯ толщины сотовой КОНСТРУКЦИИ, ТЕМПЕРАТУРЫ И РАСХОДА ХЛАДАГЕНТА СИСТЕМЫ ОХЛАЖДЕНИЯ ОБШИВКИ
При оптимизации тошцины сотовой консгрук-ции,температуры и расхода хладагента системы охлаждения обшивки вектор коэффициентов предложенной модели (1)-(27) имеет следующий вид:
ЗГг
he
Эхх
Л
Эх2 Э Т,
Х2 —¡2
= 0,
= 0,
ЭГг
he
Эх,
Щш
= 0,
■Щ-1!
Эх,
= 0,
(29)
(30)
(31)
(32)
he
'м
Эх,
х3=0
3 dThc
— А
ш
Эх,
(t)(Thc-Td{t)\ (зз)
А'З =/
В (29) - (34) использованы следующие обозначения: ф - начальное распределение температуры в сотовой конструкции; Хм - коэффициент теплопроводности дюралюминия; Ъм, Ъф - коэффициенты температуропроводности дюралюминия и воздуха, соответственно; осй ,, а)к а.г - коэффициенты теплоотдачи поверхности сотовой конструкции в канале, и между поверхности сотовой конструкции в отсеке, соответ-
ственно; Тс1, Тф - температура хладагента и воздушной среды в отсеке, соответственно.
Поскольку теплозащитная сотовая конструкция состоит из двух материалов с различными теплофизическими свойствами, то в задаче (29) - (34) мы имеем уравнение теплопроводности с разрывным коэффициентом температуропроводности [9].
Решение задачи (29) - (34) с разрывным коэффициентом температуропроводности можно получить, если решать задачу, в которой в уравнении (29) этот коэффициент Ъ(х) заменить на его гладкую аппроксимацию />'"". В качестве такой замены можно, например, рассматривать его сглаживание в окрестности поверхностей разрыва с помощью интегрального усреднения с бесконечно дифференцируемым финитным ядром [10]
Ьш'( х) = \%т\х-у\Ь(у)11у =
\х~у\<Рш
/I
4 I
со.
■г'-у|<рт
"Я
(35)
Ыу)с1(у),
В результате замены в (29) - (34) коэффициента температуропроводности функцией Ь^, приходим к следующей краевой задаче
дПт)
ы
^ эьш дт
-I
Ьс
7=1
ЭХ/ ^ ЭХ; Т1>")(х,0) = <р(х),
Эх,
э?Г
ЭХ!
ЭУГ
= 0,
эгГ
Эх,
я д1"с
= 0,
Эх!
э т^
= о.
Эх,
= 0,
\41
Эх,
лгр(т)
Эх,
=аъс,ытк-тс1т
- Щр.тг
(36)
(37)
(38)
(39)
(40)
(41)
Л'; =/
Для решений краевых задач для параболических уравнений известны вероятностные представления [11], [12].
Вероятностным представлением решения задачи (36) - (41) является математическое ожидание
Ты(х,0 = Ех%т_( {(р(Хт)¥т + 1Т), (42)
где символом Е щ обозначено математическое ожидание относительно вероятностной меры Р , соответствующей случайному процессу,
стартующему в момент времени 74 из точки х, т.е.Хг=х.
Приближённые значения можно
получить методом Монте-Карло. Для этой цели в работах [7, 13, 14] предлагалось использовать метод Эйлера. Моделирование траекторий методом Эйлера для стохастических дифференциальных уравнений (36) - (41) осуществляется по следующей схеме
(43)
х;+1 = Х;А+1 + (Л/+1 *)//,. г з - (ХгА+1), (44)
Ум = У1 Г1"+ (аг),Ъо г)К +
+ ^ (Х- Л Х;>0 ( X, К).
?+1
*г. + У:(ЫХX (Щ )%Ш (X, ; л,+1 к +
+ у,(а2), (ц/2)х-бо(х,)А1+1К),
(45)
(46)
(47)
к1+1=к, +А МК,
где г, г+1 - номера узлов сетки; |. - вектор трех независимых N (0, 1) случайных величин;
Аг+1 К = с! ( Л ) - расстояние от Ж^ до дС
в случае выхода из области.
Известно, что точка первого выхода из шара винеровского процесса, начинающего движение из центра этого шара, имеет на границе этого шара равномерное распределение. При этом эта точка не зависит от времени первого выхода на границу. В связи со сказанным, предлагается [15] для ускорения моделирования траекторий стохастических дифференциальных уравнений при их прохождении внутри ячейки использовать метод случайного блуждания по сферам, а движение по каркасу и некоторой его окрестности моделировать с помощью метода Эйлера с малым шагом.
Решение обратной задачи, то есть оценивание коэффициентов 0 модели сводится к минимизации взвешенной суммы квадратов невязок между заданными по принятому критерию значениями У и соответствующими значениями
©), полученными в ходе расчётов по уравнениям модели:
N
Ф(0) = -^А-©) гО7,-%(Ч,в)Х (48)
к=1
где tk — моменты времени при к = 1,...,М.
Как было отмечено в работе [16, 17], для минимизации функции (48) целесообразно использовать композицию метода наискорейшего спуска, квазиньютоновского метода Бройде-на-Флетчера-Гольдфарба-Шэнно [18] и метода Ньютона [19], которые реализуются в соответствии с формулой:
0;+1 -да,), (49)
где и — коэффициент, характеризующий длину шага на /-Й итерации; Б — параметр, указывающий направление поиска вектора 0(! действительных значений коэффициентов ©.
Для определения вектора коэффициентов 0 модели теплового состояния сотовой конструкции будем определять минимум функции Ф(0) взвешенной суммы квадратов невязок (уравнение (48)) с помощью итерационного алгоритма минимизации, использующего производные функции Ф (С)) . Для этой цели предлагается использовать вариант стохастического квазиградиентного алгоритма с переменной метрикой [18], в котором приближения к точке минимума строятся по правилу
0А+1 = 0А — рк НкУк Ф, к = ОД,..., (50) где № — случайная квадратная матрица размера /х/; УФ — градиент целевой функции в точке 0 к; рк — параметр шага.
Решение системы обыкновенных дифференциальных уравнений (1)...(4), выполнялось по численной схеме типа Розенброка второго порядка аппроксимации для неавтономных систем, изложенной в [20].
6. РЕЗУЛЬТАТЫ РАСЧЁТОВ
Панели фюзеляжа (рисунок 2) состоят из дюралюминиевых сотовых конструкций. Между внешней обшивкой и сотовыми конструкциями расположены титановые трубы с квадратным сечением и толщиной стенки 0,5 мм. Трубы имеют размеры 15 на 15 мм. В качестве хладагента активной теплозащиты в охлаждаемых канатах используется раствор алкоголя с водой в соотношении 30 к 70.
Внешняя обшивка выполнена из титанового сплава и имеет толщину 2 мм. Внутренняя обшивка выполнена из дюралюминия и имеет толщину 0,5 мм.
Толщина сотового заполнителя (рисунок 3) составляет 0,05 мм.
Параметры режима полёта и воздушной среды за бортом модели применения самолёта приведены на рисунке 1.
Температура внутренней поверхности сотовой конструкции в отсеке при граничных условиях первого рода были приняты 293 К.
Начальные значения толщины дюралюминиевой сотовой конструкции были приняты 50 мм, температура хладагента в канале - 233 К, скорость перемещения хладагента в канале -0,08 м/с.
Границы при параметрической идентификации следующие:
- толщина сотовой конструкции 40 -100 мм;
Рис. 3. Металлическая сотовая конструкция герметизированного отсека
- температура хладагента в канате 233 - 313 К;
- скорость перемещения хладагента в канале 0,05-0,20 м/с.
Эффективная теплопроводность сотовой конструкции была получена с помощью тепло-физического эксперимента при стационарном процессе теплопередачи. Было получено следующее значение 0,12 Вт/(м*К). Температуропроводность многослойной сотовой конструкции была определена также экспериментально и равнялась 3 10 3 м2/с.
Для расчёта теплового состояния панели была разработана параллельная программа на языке Fortran 90. Распараллеливание в программе осуществляется по схеме ведущий-ведомые (Master-Slave). В этой схеме одно вычислительное ядро считается главным, и оно распределяет весь объём работы по моделированию случайных траекторий по всем ядрам, участвующим в работе. По окончании моделирования всех траекторий все ядра передают ведущему ядру полученные результаты расчётов для вычисления математического ожидания функционала, дающего оценку температуры.
При написании параллельной программы использовалось программное обеспечение Intel MPI, Version 4.1. Моделирование траекторий случайного процесса осуществлялось с использованием параллельного датчика гауссовских случайных величии из библиотеки Intel MKL. Вычисления проводились в Сибирском Суперкомпьютерном центре на гибридном кластере HKC-30T+GPU с использованием 36 4-х ядерных процессоров Е5540 2,53 GHz.
Оценки коэффициентов модели (28) для толщины сотовой конструкции, температуры и расхода хладагента соответственно равны
0 =[0,052; 233,6; 0,002; 0,0006 f.
Вектор произведения PD погрешностей
оценок коэффициентов 0 определялся по алгоритму, приведённому в [19]. Он составляет при доверительной вероятности р = 0,99 соответственно
PD = [ 0,012; 3,2; 0,0005; 0,00008 f.
ЗАКЛЮЧЕНИЕ
Предложен теоретический метод определения потребных характеристик толщины металлической сотовой конструкции герметизированного отсека, температуры и расхода хладагента системы охлаждения обшивки гиперзвукового самолёта. Разработанная математическая модель системы герметизированного теплоизолированного отсека с системой охлаждения обшивки представлена системой дифференциальных и алгебраических уравнений.
В качестве методов параметрической идентификации модели теплового состояния отсеков предложено использовать композицию метода наискорейшего спуска, метода Ньютона и квазиньютоновского метода Бройдена-Флетчера-Голь-дфарба-Шэнно и вариант стохастического квазиградиентного алгоритма с переменной метрикой. Для решения прямой задачи сотовых конструкций применялся метод, основанный на численном решении стохастических дифференциальных уравнений, с использованием метода Эйлера и метода случайного блуждания по сферам.
Вектор произведения погрешностей оценок коэффициентов нелинейной математической модели вычислялся по атгоритму, применяющему функцию квантиль %2 - распределения и матрицы Грамма.
Определены потребные толщина металлической сотовой конструкции герметизированного отсека, температура и расход хладагента системы охлаждения обшивки гиперзвукового самолёта.
СПИСОК ЛИТЕРАТУРЫ
1. Авиационные правила. Часть 25. Нормы лётной годности самолётов транспортной категории. М.: Межгосударственный авиационный комитет, 2004. 236 с.
2. Воронин Г.И. Системы кондиционирования на летательных аппаратах. М.: Машиностроение, 1973, 443 с.
3. Михеев М.А., Михеева И.М. Основы теплопередачи. М.: Энергия, 1973, 320 с.
4. Ярышев H.A. Теоретические основы измерения нестационарных температур. Л.: Энергоатомиз-дат, 1990, 255 с.
5. Дульнев Г.Н., Тарновский H.H. Тепловые режимы электронной аппаратуры. Л.: Энергия, 1971, 248 с.
6. Малозёмов В.В., Кудрявцева Н.С. Системы терморегулирования космических аппаратов. М.: Машиностроение, 1995,107 с.
7. Calculation of heat transfer in heterogeneous
structures such as honeycomb by using numerical solution of stochastic differential equations / Gasev, S.A., Nikolaev, V.N. // Advanced Materials Research. 2014. No 1016. P. 758-763.
8. Some continuous Monte Carlo methods for the Dirichlet problem / Muller M.E. // The Annals of Mathematical Statistics. 1956. Vol. 27, No. 3. P. 569-589.
9. Ладыженская O.A., СолонниковВ.А., Уральцева H.H. Линейные и квазилинейные уравнения параболического типа. М.: Наука, 1967, 736 с.
10. Соболев С.Л. Некоторые применения функционального анализа в математической физике. М.: Наука, Изд. 3-е, 1988, 336 с.
11. Гнхман И.И., Скороход A.B. Стохастические дифференциальные уравнения и их приложения. Киев: Наукова думка, 1982.
12. Кушнер Г.Дж. Вероятностные методы аппроксимации в стохастических задачах управления и теории эллиптических уравнений. М.: Наука, 1985.
13. Применение СДУ к оценке решения уравнения теплопроводности с разрывными коэффициентами / Гусев С.А // Сиб. журн. вычисл. Математики: РАН. Сиб. отд-ние. 2015. Т. 18, № 2. С. 147-161.
14. Gusev S., Nikolayev V. Heat condition compartments of aircraft with a honeycomb structure. Saarbrucken: LAMBERT, 2017,133 p.
15. Giisev S.A., Nikolayev V.N. Estimation of the Thermal Process in the Honeycomb Panel by a Monte Carlo Method. ATCES 2017 ЮР Publishing. ЮР Conf. Series: Materials Science and Engineering 302. 2017. P. 1-6.012045 doi: 10.1088/1757-899X/302/1/012045.
16. Математическая модель конвективно-лучистого теплообмена продуваемого теплоизолированного негерметичного отсека летательного аппарата / Николаев В.Н., Гусев СЛ., Махоткин O.A. Прочность летательных аппаратов. Расчёт на прочность элементов авиационных конструкций. Науч. - техн. сб. Новосибирск: СибНИА. 1996. Вып. 1.С. 98-108.
17. Gusev S.A. andNikolayevV.N. Optimization parameters of air-conditioning and heat insulation systems of a pressurized cabins of long-distance airplanes. ATCES 2017 ЮР Publishing. ЮР Conf. Series: Materials Science and Engineering 302. 2017. P. 1-6. 012042 https://10.1088/1757-899X/302/l/012042
18. Ouasi-Newton methods for unconstrained optimization/ Gill, P., Murray, E. //J. of the institute of mathematics and its applications. 1971. Vol. 9, No 1. P. 91-108.
19. Himmelblau D.M. Applied Nonlinear Programming. Texas: McGraw-Hill Book Company, 1972.
20. Артемьев, С.С., Демидов, Г.В., Новиков, Е.А. Минимизация овражных функций численным методом для решения жёстких систем уравнений. Новосибирск: ВЦ СО АН СССР. Препринт № 74.1980.
MATHEMATICAL SIMULATION OF THE HYPERSONIC AIRCRAFT PRESSURIZED COMPARTMENT THERMAL STATE
© 2018 S.A.Gusev '», V.N. Nikolaev'
institute of Computational Mathematics and Mathematical Geophysics SB RAS, Novosibirsk
2 Novosibirsk State Technical University 3 FSUE «S.A. Chaplygin Siberian Aeronautical Research Novosibirsk Institute», Novosibirsk
A method for determining the hypersonic airplane pressurized compartment thermal state is developed based on the use of a mathematical model of the compartment with metal honeycomb constructions thermal state. The mathematical model of the system of the pressurized thermal insulated compartment with a surface cooling system is represented by a system of differential and algebraic equations describing the convective and radiation heat transfer of the skin outer surface, convective heat transfer in the cooling system channel, and heat exchange by direct conduction of the honeycomb structure. The development of a stochastic method for solvingthe direct heat exchange problem of a metal honeycomb construction using stochastic differential equations by the Euler method and random walk by spheres method is earned out. The random process trajectories simulation was carried out using the parallel sensor of Gaussian random variables from the Intel MKL library. To solve the stiff set of ordinary differential equations the implicit Rosenbrock method of second-order is used. The inverse problem for the metal honeycomb construction is solved by a stochastic quasigradient algorithm with a variable metric. In addition, parametric identification of the compartments thermal state model is proposed to be carried out by the composition of the steepest descent method, Newton's method and the quasi-Newton method of Broyden-Fletcher-Goldfarb-Shanno. The product vector of estimated coefficients errors of a nonlinear mathematical model is determined by the expression that is a quantile function %2 - distributions and the Gram matrix. At the same time, the method of joint confidence evaluation region projection on the coordinate axes of the coefficients space is used. Investigations of the hypersonic aircraft were conducted in accordance with the Airworthiness standards. The required thickness of the metal honeycomb construction of the pressurized compartment, the temperature and the refrigerant consumption of the hypersonic airplane surface cooling system are obtained. Optimization of the metal honeycomb structure of the hypersonic aircraft pressurized compartment and the surface cooling system was carried out to ensure the reliability of t he aircraft and its on-board equipment operation, comfortable conditions in the cockpit and passenger cabin. Keywords', Hypersonic airplane, honeycomb structure, surface cooling system, mathematical model.
Sergey Gasev, Doctor of Pltys. Math. Sciences, Senior Researcher of Institute of Computational Mathematics and Mathematical Geophysics SB RAS, Professor of Novosibirsk State Technical University. E-mail: sag@osmf.sscc.ra Vladimir Nikolaev, Doctor of Engineering Sciences, Head of Department. E-mail: nikvla50@mail.ni
