MATHEMATICAL SIMULATION OF THE HYPERSONIC AIRCRAFT PRESSURIZED COMPARTMENT THERMAL STATE

Full Text

ВВЕДЕНИЕ Оптимизация характеристик металлической сотовой конструкции герметизированного отсека гиперзвукового самолёт и системы охлаждения обшивки и проводилась для обеспечения надежности работы самолёта, его бортового оборудования, комфортных условий в кабине экипажа и салоне пассажиров. Параметры режима полёта и воздушной среды за бортом модели применения самолёта приведены на рисунке 1. Панели фюзеляжа самолёта (рисунок 2) состоят из титановой обшивки и дюралюминиевых сотовых конструкций. Титановая обшивка самолёта и дюралюминиевые сотовые конструкции разделены друг от друга каналами. В качестве активной теплозащиты используется жидкий хладагент. Критерии оптимизации теплофизических характеристик герметизированного отсека определены в соответствии с Нормами летной годности АП 25 [1]. В качестве математической модели теплового состояния обшивки, дюралюминиевых сотовых конструкций и системы охлаждения обшивки для герметизированного отсека примем систему уравнений, описывающих конвективный и радиационный теплообмен наружной поверхности обшивки, конвективный теплообмен в канале системы охлаждения и теплообмен теплопроводностью сотовой конструкции. 1. ВЫНУЖДЕННАЯ КОНВЕКЦИЯ СНАРУЖИ ГЕРМЕТИЗИРОВАННОГО ОТСЕКА Температура воздуха за бортом у поверхности обшивки T* определяется по формуле , (1) где Tair,out - температура воздуха за бортом за пределами теплового пограничного слоя. Температура восстановления Te в уравнении (1) определяется по выражению (2) где r - коэффициент восстановления температуры (для ламинарного пограничного слоя , для турбулентного - ; - число Маха на внешней границе пограничного слоя; k=cp/cV - отношение удельной теплоёмкости газа при постоянном давлении и объёме. Коэффициент теплоотдачи αcv,out наружной поверхности обшивки будем вычислять по формуле [2] (3) для ламинарного пограничного слоя при критерии Рейнольдса < 1.106 и по формуле [2] (4) для турбулентного пограничного слоя при критерии Рейнольдса 1.106 1.109. В уравнениях (3), (4) использованы следующие обозначения: , - критерии соответственно Рейнольдса и Прандтля при температуре воздуха ; - плотность воздуха за бортом при температуре T*; - удельная теплоёмкость воздуха при температуре T*; Vair,out - воздушная скорость полёта. Критерий вычисляется по формуле , (5) где Lcv,out - характерный размер для наружной поверхности обшивки; - динамическая вязкость воздуха. Характерный размер Lcv,out принимается равным расстоянию от носовой точки фюзеляжа до середины рассчитываемого отсека или его части. Динамическая вязкость определяется по выражению = 0,1222229.10-6 + 0,682674.10-8 - 0,313155.10-11 ()2. (6) Критерий вычисляется по формуле , (7) где - теплопроводность воздуха. Теплопроводность определяется по выражению = 0,141483.10-2 + 0,896161.10-4 - - 0,204759.10-7 ()2 . (8) Плотность воздуха за бортом вычисляется по выражению , (9) где pair,out - давление воздуха за бортом. Тепловой поток Qcv,out , поступающий к наружной поверхности обшивки от прямого и отраженного излучения Солнца и прямого излучения Земли определяется по уравнению (10) где c0 - постоянная Стефана-Больцмана; εcv,out - коэффициент черноты излучения наружной поверхности обшивки; Fcv- площадь обшивки при наружном конвективном теплообмене; Tcv,out - температура наружного слоя обшивки. 2. ТЕПЛООБМЕН В КАНАЛЕ СИСТЕМЫ ОХЛАЖДЕНИЯ Теплообмен в канале системы охлаждения будет определяться теплообменом наружной поверхности обшивки, теплообменом сотовой конструкции и обыкновенным дифференциальным уравнением, описывающим конвективный теплообмен внутренней поверхности канала и перенос теплоты из одной части канала в другую (11) где Tcl,k-1, Tcl,k - температуры хладагента соответственно в (k-1)-ой и k-ой частях канала; Thc,cl,k - температура поверхности сотовой конструкции на границе с k-ой частью канала; αk - коэффициент теплоотдачи внутренней поверхности k-ой части канала; Fk - площадь теплоотдающей поверхности k-ой части канала; Jcl,c - массовая скорость хладагента в канале; cp,cl - удельная теплоёмкость хладагента; Fc - площадь поверхности сечения канала; Ccl,k - теплоёмкость хладагента в k-ой части канала. Tcl,k в выражении (11) с индексом t означает её дифференцирование по времени t. Теплоёмкость хладагента Ccl,k определяется по выражению (12) где ρcl,k - плотность хладагента в k-ой части канала; Wcl,c - скорость хладагента в канале; Δt - интервал дискретизации времени при решении системы дифференциальных уравнений; Vcl,k - объём хладагента в k-ой части канала. Коэффициент теплоотдачи αk внутренней поверхности канала определяется уравнениями, приведёнными в [3]. Ламинарное течение при критерии Рейнольдса Recl < 2100 описывается уравнением (13) Критерий Нуссельта Nucl,d,av представляет собой безразмерный средний коэффициент теплоотдачи и определяется из выражения (14) где αcl,d,av - средний коэффициент конвективной теплоотдачи; d - эквивалентный диаметр канала, равный учетверённой площади поперечного сечения канала, делённой на его полный «смоченный» периметр [4]; λcl - теплопроводность хладагента. Критерий Рейнольдса Recl,d характеризует режим движения хладагента и определяется из выражения , (15) где Wcl,,av - средняя скорость хладагента; ρcl - плотность хладагента; μcl - динамическая вязкость хладагента. Критерий Прандтля Prcl определяет физические свойства хладагента при определяющей средней температуре хладагента и вычисляется из выражения (16) Критерий Прандтля Prk определяет физические свойства хладагента при определяющей средней температуре поверхности канала. Критерий Грасгофа Grcl,d характеризует подъёмную силу, возникшую вследствие разности плотности хладагента, и определяется из выражения (17) где βcl - коэффициент объёмного расширения хладагента; g - ускорение в поле тяготения; ΔTcl - разность между температурой поверхности канала и хладагента за пределами пограничного слоя; v - кинематическая вязкость хладагента. Коэффициент εl учитывает изменение среднего коэффициента теплоотдачи по длине канала и зависит от отношения длины канала к его эквивалентному диаметру l/d. Значения этого коэффициента представлены в таблице 1. Переходной режим течения при критерии Рейнольдса 2100 < Recl < 1.104 описывается уравнением (18) где значения коэффициента K0 зависят от критерия Рейнольдса Recl и представлены в таблице 2. Турбулентное течение при критерии Рейнольдса Recl > 1.104 описывается уравнением (19) Расход хладагента Gcl предлагается представить линейной регрессией вида , (20) где - коэффициенты модели, которые будут оцениваться. 3. ТЕПЛООБМЕН ТЕПЛОПРОВОДНОСТЬЮ СОТОВОЙ КОНСТРУКЦИИ Систему уравнений теплообмена сотовой конструкции представим в виде одномерных уравнений теплопроводности, описывающих процесс передачи тепла в многослойной сотовой конструкции [2, 5] где то есть коэффициенты Chc, λhc зависят от того, в каком слое рассматривается перенос тепла. В уравнениях (21) - (24) использованы следующие обозначения: Thc (t, x) - температура сотовой конструкции; Thc,t - первая производная Thc по t; Thc,x - первая производная Thc по ; Thc,x,x - вторая производная Thc по ; Chc(x) - объёмная теплоёмкость сотовой конструкции обшивки, определяемая теплоёмкостью дюралюминия Ccompo и теплоёмкостью воздуха Cair; λhc(l) - теплопроводность сотовой конструкции, определяемая теплопроводностью дюралюминия λdur и теплопроводностью воздуха λair; αhc,k - коэффициент теплоотдачи наружной поверхности сотовой конструкции; αcv,in - коэффициент теплоотдачи внутренней поверхности сотовой конструкции; Fhc,k - площадь сотовой конструкции при наружном и внутреннем теплообмене; Tair - температура воздушной среды в отсеке или в части отсека; l - толщина сотовой конструкции. 4. ВЫНУЖДЕННАЯ КОНВЕКЦИЯ В ГЕРМЕТИЗИРОВАННОМ ОТСЕКЕ Коэффициент теплоотдачи αhc,in внутренней поверхности сотовой конструкции можно вычислить по формуле [5] (25) для ламинарного пограничного слоя при критерии Рейнольдса Rein < 4.104 и по формуле [6]: (26) при критерии Рейнольдса Rein 4.108 . Теплопроводность λair и динамическая вязкость μin в формулах (25), (26) в критериях Re, Pr определяются соответственно по выражениям (8), (6) для температуры воздуха в отсеке Tair. Произведение JW плотности ρair и скорости Wair воздуха в отсеке или массовая скорость JW в критерии (27) принимается из результатов эксперимента. Характерный размер Lhc,in принимается по аналогии с Lhc,out равным расстоянию от начала отсека до середины рассчитываемой части отсека. 5. ОПТИМИЗАЦИЯ ТОЛЩИНЫ СОТОВОЙ КОНСТРУКЦИИ, ТЕМПЕРАТУРЫ И РАСХОДА ХЛАДАГЕНТА СИСТЕМЫ ОХЛАЖДЕНИЯ ОБШИВКИ При оптимизации толщины сотовой конструкции, температуры и расхода хладагента системы охлаждения обшивки вектор коэффициентов предложенной модели (1)-(27) имеет следующий вид: (28) и включает в себя необходимые характеристики (значения толщины сотовой конструкции Lhc в м, температуры хладагента Tcl в К и коэффициенты , описывающие расход хладагента Gcl в кг/с). Для решения прямой задачи теплового состояния сотовой конструкции фюзеляжа рассмотрим математическую модель теплообмена в гетерогенной структуре, как параболическую краевую задачу с разрывными коэффициентами. Описание и физические свойства гетерогенных структур можно найти в [7]. При этом делалось сглаживание разрывных коэффициентов уравнения теплопроводности на основе интегрального усреднения, а численное решение стохастических дифференциальных уравнений осуществлялось методом Эйлера и методом случайного блуждания по сферам [8]. Предполагается, что рассматриваемый процесс теплопередачи происходит на отрезке времени [0, tend] и описывается следующей краевой задачей для уравнения теплопроводности где (29) В (29) - (34) использованы следующие обозначения: φ - начальное распределение температуры в сотовой конструкции; λAl - коэффициент теплопроводности дюралюминия; bAl, bair - коэффициенты температуропроводности дюралюминия и воздуха, соответственно; αhc,cl, αhc,air - коэффициенты теплоотдачи поверхности сотовой конструкции в канале, и между поверхности сотовой конструкции в отсеке, соответственно; Tcl , Tair - температура хладагента и воздушной среды в отсеке, соответственно. Поскольку теплозащитная сотовая конструкция состоит из двух материалов с различными теплофизическими свойствами, то в задаче (29) - (34) мы имеем уравнение теплопроводности с разрывным коэффициентом температуропроводности [9]. Решение задачи (29) - (34) с разрывным коэффициентом температуропроводности можно получить, если решать задачу, в которой в уравнении (29) этот коэффициент b(x) заменить на его гладкую аппроксимацию . В качестве такой замены можно, например, рассматривать его сглаживание в окрестности поверхностей разрыва с помощью интегрального усреднения с бесконечно дифференцируемым финитным ядром [10] (35) В результате замены в (29) - (34) коэффициента температуропроводности функцией , приходим к следующей краевой задаче Для решений краевых задач для параболических уравнений известны вероятностные представления [11], [12]. Вероятностным представлением решения задачи (36) - (41) является математическое ожидание (42) где символом Ex,T-t обозначено математическое ожидание относительно вероятностной меры Px,T-t , соответствующей случайному процессу, стартующему в момент времени T-t из точки x, т.е. XT-t=x. Приближённые значения можно получить методом Монте-Карло. Для этой цели в работах [7, 13, 14] предлагалось использовать метод Эйлера. Моделирование траекторий методом Эйлера для стохастических дифференциальных уравнений (36) - (41) осуществляется по следующей схеме где i, i+1 - номера узлов сетки; ξi - вектор трех независимых N (0, 1) случайных величин; - расстояние от до ∂G в случае выхода из области. Известно, что точка первого выхода из шара винеровского процесса, начинающего движение из центра этого шара, имеет на границе этого шара равномерное распределение. При этом эта точка не зависит от времени первого выхода на границу. В связи со сказанным, предлагается [15] для ускорения моделирования траекторий стохастических дифференциальных уравнений при их прохождении внутри ячейки использовать метод случайного блуждания по сферам, а движение по каркасу и некоторой его окрестности моделировать с помощью метода Эйлера с малым шагом. Решение обратной задачи, то есть оценивание коэффициентов Θ модели сводится к минимизации взвешенной суммы квадратов невязок между заданными по принятому критерию значениями и соответствующими значениями , полученными в ходе расчётов по уравнениям модели: (48) где tk - моменты времени при k = 1,...,N. Как было отмечено в работе [16, 17], для минимизации функции (48) целесообразно использовать композицию метода наискорейшего спуска, квазиньютоновского метода Бройдена-Флетчера-Гольдфарба-Шэнно [18] и метода Ньютона [19], которые реализуются в соответствии с формулой: (49) где - коэффициент, характеризующий длину шага на j-й итерации; S - параметр, указывающий направление поиска вектора действительных значений коэффициентов . Для определения вектора коэффициентов модели теплового состояния сотовой конструкции будем определять минимум функции взвешенной суммы квадратов невязок (уравнение (48)) с помощью итерационного алгоритма минимизации, использующего производные функции . Для этой цели предлагается использовать вариант стохастического квазиградиентного алгоритма с переменной метрикой [18], в котором приближения к точке минимума строятся по правилу (50) где Hk - случайная квадратная матрица размера - градиент целевой функции в точке k; ρk - параметр шага. Решение системы обыкновенных дифференциальных уравнений (1)…(4), выполнялось по численной схеме типа Розенброка второго порядка аппроксимации для неавтономных систем, изложенной в [20]. 6. РЕЗУЛЬТАТЫ РАСЧЁТОВ Панели фюзеляжа (рисунок 2) состоят из дюралюминиевых сотовых конструкций. Между внешней обшивкой и сотовыми конструкциями расположены титановые трубы с квадратным сечением и толщиной стенки 0,5 мм. Трубы имеют размеры 15 на 15 мм. В качестве хладагента активной теплозащиты в охлаждаемых каналах используется раствор алкоголя с водой в соотношении 30 к 70. Внешняя обшивка выполнена из титанового сплава и имеет толщину 2 мм. Внутренняя обшивка выполнена из дюралюминия и имеет толщину 0,5 мм. Толщина сотового заполнителя (рисунок 3) составляет 0,05 мм. Параметры режима полёта и воздушной среды за бортом модели применения самолёта приведены на рисунке 1. Температура внутренней поверхности сотовой конструкции в отсеке при граничных условиях первого рода были приняты 293 К. Начальные значения толщины дюралюминиевой сотовой конструкции были приняты 50 мм, температура хладагента в канале - 233 K, скорость перемещения хладагента в канале - 0,08 м / c. Границы при параметрической идентификации следующие: - толщина сотовой конструкции 40 - 100 мм; - температура хладагента в канале 233 - 313 K; - скорость перемещения хладагента в канале 0,05- 0,20 м / c. Эффективная теплопроводность сотовой конструкции была получена с помощью теплофизического эксперимента при стационарном процессе теплопередачи. Было получено следующее значение 0,12 Вт/(м*К). Температуропроводность многослойной сотовой конструкции была определена также экспериментально и равнялась 3 10-3 м2/c. Для расчёта теплового состояния панели была разработана параллельная программа на языке Fortran 90. Распараллеливание в программе осуществляется по схеме ведущий-ведомые (Master-Slave). В этой схеме одно вычислительное ядро считается главным, и оно распределяет весь объём работы по моделированию случайных траекторий по всем ядрам, участвующим в работе. По окончании моделирования всех траекторий все ядра передают ведущему ядру полученные результаты расчётов для вычисления математического ожидания функционала, дающего оценку температуры. При написании параллельной программы использовалось программное обеспечение Intel MPI, Version 4.1. Моделирование траекторий случайного процесса осуществлялось с использованием параллельного датчика гауссовских случайных величин из библиотеки Intel MKL. Вычисления проводились в Сибирском Суперкомпьютерном центре на гибридном кластере HKC-30T+GPU с использованием 36 4-х ядерных процессоров E5540 2,53 GHz. Оценки коэффициентов модели (28) для толщины сотовой конструкции, температуры и расхода хладагента соответственно равны = . Вектор произведения погрешностей оценок коэффициентов определялся по алгоритму, приведённому в [19]. Он составляет при доверительной вероятности β = 0,99 соответственно = . ЗАКЛЮЧЕНИЕ Предложен теоретический метод определения потребных характеристик толщины металлической сотовой конструкции герметизированного отсека, температуры и расхода хладагента системы охлаждения обшивки гиперзвукового самолёта. Разработанная математическая модель системы герметизированного теплоизолированного отсека с системой охлаждения обшивки представлена системой дифференциальных и алгебраических уравнений. В качестве методов параметрической идентификации модели теплового состояния отсеков предложено использовать композицию метода наискорейшего спуска, метода Ньютона и квазиньютоновского метода Бройдена-Флетчера-Гольдфарба-Шэнно и вариант стохастического квазиградиентного алгоритма с переменной метрикой. Для решения прямой задачи сотовых конструкций применялся метод, основанный на численном решении стохастических дифференциальных уравнений, c использованием метода Эйлера и метода случайного блуждания по сферам. Вектор произведения погрешностей оценок коэффициентов нелинейной математической модели вычислялся по алгоритму, применяющему функцию квантиль χ2 - распределения и матрицы Грамма. Определены потребные толщина металлической сотовой конструкции герметизированного отсека, температура и расход хладагента системы охлаждения обшивки гиперзвукового самолёта. Рис. 1. Параметры режима полёта и воздушной среды за бортом модели применения гиперзвукового самолета: Ppair,out - давление воздуха за бортом за пределами теплового пограничного слоя; M - число Маха на внешней границе пограничного слоя; Tair,out - температура воздуха за бортом за пределами теплового пограничного слоя Рис. 2. Схема расположения элементов сотовой теплоизоляции и активнойтеплозащиты гиперзвукового самолёта: Vair,out - воздушная скорость полёта; ρV - плотность воздуха за бортом; Tair,out - температура воздуха за бортом за пределами теплового пограничного слоя; ρcl - плотность хладагента; Wcl - скорость хладагента; Tcl - температура хладагента; ρair - плотность воздуха в отсеке; Wair - скорость воздуха в отсеке; Tair - температура воздуха в отсеке Таблица 1. Значения εl при ламинарном течении Таблица 2. Значения коэффициента K0 от критерия Рейнольдса Recl Таблица 3. Значения εl при турбулентном течении (21) (22) (23) (24) (30) (31) (32) (33) (34) (36) (37) (38) (39) (40) (41) (43) (44) (45) (46) (47) Рис. 3. Металлическая сотовая конструкция герметизированного отсека

About the authors

Sergey Anatol'evich Gusev

Institute of Computational Mathematics and Mathematical Geophysics SB RAS; Novosibirsk State Technical University

Email: sag@osmf.sscc.ru
Doctor of Phys.-Math. Sciences, Senior Researcher of Institute of Computational Mathematics and Mathematical Geophysics SB RAS, Professor of Novosibirsk State Technical University Novosibirsk

Vladimir Nikolaevich Nikolaev

FSUE «S.A. Chaplygin Siberian Aeronautical Research Novosibirsk Institute»

Email: nikvla50@mail.ru
Doctor of Engineering Sciences, Head of Department Novosibirsk

References

Supplementary files