CC BY
38
УДК533.6.011.6
РАЗРАБОТКА ТЕХНОЛОГИЙ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ СОПРЯЖЕННЫХ ЗАДАЧ ГАЗОВОЙ ДИНАМИКИ И ТЕРМОДИНАМИКИ КОМПОЗИТНЫХ КОНСТРУКЦИЙ ПЕРСПЕКТИВНЫХ ВЫСОКОСКОРОСТНЫХ ЛЕТАТЕЛЬНЫХ АППАРАТОВ
© 2016 Ю.И. Димитриенко, М.Н. Коряков, А.А. Захаров
Московский государственный технический университет имени Н.Э. Баумана
Статья поступила в редакцию 03.11.2016
Предложен алгоритм численного моделирования сопряженных аэрогазодинамических и термомеханических процессов в композитных конструкциях высокоскоростных летательных аппаратов, который позволяет рассчитывать все параметры трехмерного аэрогазодинамического потока в окрестности поверхности аппарата, теплообмен на поверхности, процессы внутреннего тепломассопереноса в конструкции из термодеструктирующего полимерного композитного материала, а также динамические процессы термодеформирования композитной конструкции, включающие в себя эффекты изменения упругих характеристик композита, переменную тепловую деформацию, усадку, вызванную термодеструкцией, образование внутрипорового давления газов в композите. Приведен пример численного моделирования сопряженных процессов в модельной композитной конструкции высокоскоростного летательного аппарата. Показано, что вследствие высоких температур аэродинамического нагрева конструкций, в них может происходить терморазрушение из-за термодеструкции полимерной фазы и интенсивного внутреннего газообразования в материалах конструкции.
Ключевые слова: сопряжённый процесс, аэрогазодинамика, термомеханика, гиперзвук, композиционный материал, термонапряжение
В настоящее время в России и за рубежом проявляется большой интерес к созданию гиперзвуковых летательных аппаратов, способных двигаться в атмосфере Земли продолжительное время со скоростью, которая превышает скорость звука более чем в 6 раз (число Маха M > 6). Среди первых разработок в данной области можно отметить аппарат X-51, созданный в США, успешный запуск которого был произведен в августе 2012 г. (во время экспериментального полета X-51 развил скорость в 5 Махов). Движение летательных аппаратов в атмосфере Земли с гиперзвуковой скоростью неизбежно приводит к появлению высоких температур (2000°К и выше) газа в окрестности затупленных элементов конструкции летательного аппарата, что в свою очередь обуславливает наличие высоких тепловых потоков (от 0,5 МВт/м2) к поверхности тела. Для теплозащиты элементов конструкции летательного аппарата часто используют полимерные композиционные материалы (угле-, стекло-,
Димитриенко Юрий Иванович, доктор физико-
математических наук, профессор. E-mail: dimit@bmstu.ru
Коряков Михаил Николаевич, ассистент. E-mail: mkoryakov@bmstu. ru
Захаров Андрей Алексеевич, кандидат физико-
математических наук, доцент. E-mail: azaharov@bmstu. ru
органопластики), которые в процессе эксплуатации термодеструктируют и подвергаются механической эрозии. Выбор оптимального теплозащитного материала является важной задачей на этапе проектирования летательного аппарата.
Освоение гиперзвуковых скоростей требует проведения комплексных исследований, в которых можно выделить такие составные части, как исследование гиперзвуковой аэродинамики полета и теплообмена на поверхности конструкций гиперзвуковых летательных аппаратов, исследование теплофизики и термопрочности материалов конструкций, разработка новых материалов конструкций, а также проблемы гиперзвуковой аэроупругости, управления и др. Экспериментальные методы исследований этих задач во многих случаях затруднены, поэтому основным инструментом проектировщика является численное моделирование. Задача анализа поведения теплозащитного материала должна решаться в совокупности с задачей газовой динамики для того, чтобы учесть взаимное влияние набегающего газового потока и продуктов терморазложения композиционного материала. Существующие коммерческие пакеты не позволяют решить задачу в такой сложной сопряженной постановке, поэтому необходимо разрабатывать специализированное математическое и программное обеспечение.
Математическая постановка задачи.
Рассмотрим процесс обтекания конструкции летательного аппарата, имеющего конусообразную неосесимметричную форму, высокоскоростным набегающим потоком. Общая постановка сопряженной задачи аэродинамики и термомеханики состоит из 3-х систем уравнений:
1) уравнений Навье-Стокса внешнего газового потока, обтекающего конструкцию;
2) уравнений внутреннего тепломассо-переноса в конструкции;
3) уравнений термоупругости оболочечной конструкции.
Влияние газообразных продуктов терморазложения и изменения геометрии конструкции вследствие термодеформации на течение внешнего газового потока не учитывается.
A. Система уравнений аэротермодинамики
Система уравнений вязкого теплопроводного сжимаемого газа состоит из уравнения неразрывности, уравнений движения и уравнения энергии [1]. Граничное условие на твёрдой стенке, которая является границей раздела газовой и твёрдой сред имеет вид:
V = о, е = е„,
-хуе • п = -х„уе„ • п+ае! - в„ае:,
где V - вектор скорости, 0 - температура газа, -температура твердой стенки (совпадает с температурой газа на этой стенке), X и ^ - коэффициенты теплопроводности газа и твёрдой поверхности соответственно, 0тах - максимальная
температура газа в пограничном слое, уе: -градиент температуры со стороны твердой поверхности, £g и £w - интегральные коэффициенты излучения нагретого газа и твердой поверхности и с - коэффициент Стефана-Больцмана.
B. Система уравнений внутреннего тепло-массопереноса
Будем рассматривать конструкцию, обтекаемую внешним потоком, изготовленную из композиционного материала на полимерной матрице с термостойкими керамическими волокнами. В матрице такого композита при нагреве до высоких температур, характерных для аэродинамического нагрева, происходят физико-химические процессы термодеструкции, сопровождающиеся образованием газообразных продуктов терморазложения, которые накапливаются в порах материала и отфильтровываются во внешний газовый поток, а также образованием новой твердой фазы - пиролитической фазы матрицы, которая обладает существенно более низкими упруго-прочностными характеристиками, чем исходная полимерная фаза. Четырех-фазная модель для описания внутреннего
тепломассопереноса и деформирования такого композита имеет вид [2]:
1) уравнения изменения массы полимерной фазы матрицы:
дфь
Ы
= - J;
(1)
2) уравнения фильтрации газообразных продуктов термодеструкции в порах композиционного материала:
Фв Фв
ы
+ УРв Фе ^ = J Г;
(2)
3) уравнения теплопереноса в термодеструк-тирующем композите:
де
рс ¥=•4 - свуе-рв Фвв- -1 Ле°;
(3)
где фь, фg - объемные концентрации исходной полимерной матрицы и газовой фазы; рь -плотность исходной полимерной матрицы (полагается постоянной); рg - среднее по поре значение плотности газовой фазы (переменная величина); Cg - удельная теплоемкость газовой фазы при постоянном объеме, рис- плотность и удельная теплоемкость композита в целом, q -вектор плотности теплового потока, 0 - температура композита, общая для всех фаз; Vg -вектор скорости движения газовой фазы в порах; Де° - удельная теплота термодеструкции матрицы; ] - массовая скорость термодеструкции матрицы и Г - коэффициент газификации матрицы.
К уравнениям (1)-(3) присоединяются определяющие соотношения, связывающие вектор-функции q и Vg с функциями Уе и Ур с помощью законов Фурье и Дарси, а также соотношение Аррениуса для ] и уравнения Менделеева-Клайперона для порового давления газовой фазы р:
я = -л • уе, рй фв Vг = -к •Ур,
Ел
-1=^М-^> р=
Рв ^е
(4)
где /° - предэкспоненциальный множитель, Еа -энергия активации процесса термодеструкции, а Л - тензор теплопроводности и К - тензор газопроницаемости композита зависят от концентраций фаз [2].
На нагреваемой части поверхности конструкции граничные условия для уравнений(1)-(3) выглядят следующим образом: р = ре, 0 = 0w, где ре - местное давление внешнего газового потока на поверхности, а 0w - температура поверхности композита. На остальной части поверхности
композита задаются граничные условия герметичности и теплоизоляции: п -Ур = 0, -п - Л-У0 = 0.
C. Система уравнений термомеханики оболочечной конструкции
В криволинейной системе координат Оффф, связанной со срединной поверхностью оболочечной конструкции летательного аппарата, система уравнений термомеханики включает в себя [2]:
1) уравнения равновесия оболочки:
дг( 4Тш )+^( АТ )-дт т +
91а
дА дР
+^ т+ал ш- л =о,
( АвМ аа ( АМ "в " ^ М вв +
дА дР
+—мав - ЛАва - А —~ = 0, - АА (КТ + к2т22 )+дАа +дАа - РеАА -
дд1 дд2
-(k + k) аАф8 р = 0;
2) кинематические соотношения:
1 dU а 1 dAa
ваа =--+---'
Ш А ^а A1A2 ^ß 1
2eа3 ~ ^ + Ya kaUa
Aa дч*
1 3U1 1 dU2 1 f 9Д SA
(5)
U + ^ ^,
2e12 = — ^+—^--—I U1 + U.
A dq2 A dq1 АA \dq2 dq1
_ 1 dYa , 1 dAa
A dq A A2 dqß
2k a3 = ka Ya >
2Kn = ± К+1 dYi f di Y +dA Yi
A2 3q2 A dqL AA2 ^dq2 3qL
(6)
3) определяющие соотношения оболочки:
2
Taa = ^ (Caßßßß + NaßKßß ) — Pga - Ta, ß=1
^2 = 2(C66ei2 + N66 Ki2),
2
Maa =Z (Naß^ßß + DaßKßß ) - Mga - Ma,
ß=1
M12 = 2( N66ei2 + D66 K12)' öa = C+3,a+3e a 3;
Ла, ка - параметры первой квадратичной формы и главные кривизны срединной поверхности оболочки, а,в = 1,2.
Введены также обозначения для усилий и
моментов тепловых напряжений та , Ма, зависящих от температурных деформаций Ба оболочки:
Т =ус б(0) ММ =Ус Б(1) в=1 в=1
А/
SO) = fcß "
J aeisßqjdq.
з' ^з^) = I ae2S3i3™43>
S у = (af 9f +аьфьПу )(e-e0) +
t
+apQy J (e(t) - е(т))фpd t - ßpфpQy
(8)
где а, аь, ар - коэффициенты теплового расширения волокна, полимера и пиролитической фазы матрицы, вр - коэффициент усадки, Бу, Оу -коэффициенты, зависящие от расположения волокон в композите, у = 1,2,3 [2].
Вследствие размягчения полимерной матрицы и её термодеструкции, жесткости оболочки изменяются при нагреве. Учет этого изменения для ортотропных композитных оболочек осуществляется с помощью 2-х функций ал, а02 [2]. Деформации еар и напряжения сар в оболочке вычисляются по следующим формулам:
Бав = еав + 43Кав , Б33 = 0, Ба3 = ва3, °аа = "/аР + Й01 + ^3КРР - БР
СТ12 =а01С66 (Б12 +^3К12 ), Р=1,2.
Для поперечного нормального напряжения сзз и напряжений межслойного сдвига саз имеем следующие формулы [2]:
(Pl + p2 ) + Lr \J0)„ _(1)K
О (0)
°зз = 6h| V"L + hC3i I аме11 + ae1) K11 -S1 |-
P1 1 f О (0) ^ 1 О (0) ^
—11 +—C0 | a„(iV, + a(1)^ -S2 |--C0 S3 +
(7)
где Гаа, raß, M„a, Maß - усилия и моменты в оболочке; Oa - перерезывающие усилия; eaa, ea3, ei2 - деформации срединной поверхности оболочки; Kaa, ксз, К12 - искривления срединной поверхности; Ua, ya , W - перемещения, углы искривления и прогиб срединной поверхности;
+(P2 -Pi)^(q3) + Pl 2 P2 +ФеP>
a = 12h(q3) Cо e a(0) o = 12h(q3) Cо e (0) h h
«q.) = i-q^. h(q3) = 4-i^
(9)
где р1 и р2 - давления на внешних поверхностях оболочки.
Разработка метода решения сопряженной задачи. Для численного решения
сопряжённой задачи используется пошаговый метод с двумя шагами по времени Ati и At2. Максимальный крупный шаг Ati по времени используется для изменения граничных условий задачи — входных данных набегающего потока и решения задачи внутреннего тепломассопере-носа. На каждом шаге Ati решение осуществляется в четыре этапа.
Этап 1. Решение системы уравнений внутреннего тепломассопереноса (1)-(3). Это решение осуществляется численным конечно-разностным методом с использованием метода линеаризации и неявной разностной схемы, детали которого описаны в [2]. Температура поверхности конструкции, взаимодействующей с набегающим газовым потоком, на этом этапе полагается известной и берется с предыдущего временного шага Ati.
Этап 2. Осуществляется цикл решения системы уравнений Навье-Стокса с мелким шагом At2 до установления потока. Для решения систем дифференциальных уравнений Навье-Стокса используется метод конечных объемов на основе элементов, центрированных относительно узлов сетки (vertex-centered volume) [3].
Этап 3. Осуществляется решение системы уравнений термоупругости оболочечной конструкции аппарата (5)-(7) с помощью метода конечного элемента. Входными данными для этой задачи являются поля давлений на внешней pi и внутренней p2 поверхности оболочки, которые определяются после решения уравнений газовой динамики, а также распределение температуры 0, объёмных концентраций фаз и порового давления p газообразных продуктов терморазложения композитной оболочки, которые рассчитываются при решении уравнений внутреннего тепломассопереноса (i)-(3).
Этап 4. После решения задачи термоупругости осуществляется расчёт термонапряжений в оболочке с помощью формул (9).
Результаты численного решения. В работе представлены результаты численного моделирования обтекания фрагмента корпуса модельного летательного аппарата эллипсоидальной формы гиперзвуковым потоком газа (М=6) на высоте 15 км [4, 5]. На рис. 1-4 представлены распределения тепловых деформаций ei и £з, порового давления p газообразных продуктов терморазложения полимерной фазы и поперечного напряжения сзз вдоль безразмерной координаты по толщине оболочки на нижней поверхности в точке, удаленной от критической на расстояние 30Я, где R - радиус затупления носовой части.
Характерной особенностью функций ei и ез является их немонотонность: в зоне повышенной температуры (по толщине оболочки) происходит
рост значений тепловых деформаций, однако при переходе через условную точку начала интенсивного терморазложения происходит уменьшение тепловой деформации, обусловленное образованием пиролитического остатка. График распределения сзз показывает, что в подповерхностном слое после начала терморазложения полимера возникает значительный максимум напряжения, обусловленный ростом порового давления газообразных продуктов термодеструкции, а также возникновением усадочных деформаций. Значения поперечных напряжений на нижней части оболочки достигают значений в 0,13 ГПа, что значительно превышает предел прочности композитной оболочки в поперечном направлении. В результате, в этой части оболочки может возникнуть разрушение по типу расслоения, при котором верхние слои ткани композита отслоятся от остальной части материала. Аналогичное явление может произойти и на верхней части оболочки, поскольку поперечные напряжения начинают расти и достигают значений в 0,06 ГПа, что хотя и меньше значений в нижней части оболочки, но также сравнимо с пределом прочности материала в поперечном направлении.
Исследование выполнено при поддержке гранта Президента РФМК-3007.2015.8.
012т...............]...............;...............т...............г...............г
0.1 0.08 0-06 0.04
0.02 4................|................;................•................)................
0 0.2 0.4 0.6 0.8 1 Рис. 1. Тепловая деформация £1
0 0 2 0 4 0 С 0LS V
Рис. 2. Тепловая деформация ез
СПИСОК ЛИТЕРАТУРЫ:
Димитриенко, Ю.И. Метод ленточных адаптивных сеток для численного моделирования в газовой динамике / Ю.И. Димитриенко, В.П. Котенев, А.А. Захаров. - М.: Физматлит, 2011. 280 с. Димитриенко, Ю.И. Механика композиционных материалов при высоких температурах. - М.: Машиностроение, 1997. 366 с. Волков, К.Н. Течения и теплообмен в каналах и вращающихся полостях / К.Н. Волков, В.Н. Емельянов. - М.: Физматлит, 2010. 488 с. Димитриенко, Ю.И. Моделирование сопряженных процессов аэрогазодинамики и теплообмена на поверхности теплозащиты перспективных гиперзвуковых летательных аппаратов / Ю.И. Димитриенко, А.А. Захаров, М.Н. Коряков, Е.К. Сьз-дыков // Известия высших учебных заведений. Машиностроение. 2014. № 3(648), С. 23-34. Димитриенко, Ю.И. Численное моделирование сопряженных аэрогазодинамических и термомеханических процессов в композитных конструкциях высокоскоростных летательных аппаратов / Ю.И. Димитриенко, М.Н. Коряков, А.А. Захаров, А.С. Строганов // Математическое моделирование и численные методы. 2014. № 3. C. 3-24.
DEVELOPMENT OF COMPUTATIONAL METHODS FOR NUMERICAL SIMULATION OF CONJUGATED PROBLEMS OF GAS DYNAMICS AND THERMOMECHANICS IN COMPOSITE STRUCTURES OF PROMISING HIGH-SPEED AIRCRAFTS
© 2016 Yu.I. Dimitrienko, M.N. Koryakov, A.A. Zakharov
1 Bauman Moscow State Technical University
An algorithm for the numerical simulation of coupled aerodynamic and thermomechanical processes in composite structures of high-speed aircraft is proposed, which allows to calculate all parameters of the three-dimensional aerogasdynamic flow near the surface of the aircraft, the heat exchange on the surface, the internal heat-and-mass transfer processes in the construction of thermodecomposition polymer composite material, and dynamic processes of heat displacement in the composite construction, including the effects of changes in the elastic characteristics of the composite, variable thermal deformation, shrinkage caused by thermodecomposition, the formation of the internal pore pressure in the composite. An example of the numerical solution of the coupled problem of a model composite construction of high-speed aircraft, showing the possibility of the proposed algorithm is given. It is shown that due to the high temperatures of the aerodynamic heating of the structures there can appear a polymer phase thermodecomposition because of the polymer phase thermodecomposition and intensive internal gas generation into the structure of the material.
Key words: conjugated problem, computational fluid dynamics, aerothermodynamics, hypersonic flows, thermomechanics, polymer composites, thermodecomposition structure, heat shield
Yuriy Dimitrienko, Doctor of Physics and Mathematics, Professor. E-mail: dimit@bmstu.ru Mikhail Koryakov, Assistant. E-mail: mkoryakov@bmstu. ru
Andrey Zakharov, Candidate of Physics and Mathematics Associate Professor. E-mail: azaharov@bmstu.ru
Рис. 3. Поровое давление p, атм
0.1 0.08 0.06 0.04 Ü.02 0
-0.02 -0.04
(
....................
0 0 2 ok 0 6 0 8
1
1.
3.
4.
Рис. 4. Поперечное напряжение сзз, ГПа
