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

Димитриенко Ю.И.; Коряков М.Н.; Захаров А.А.; Строганов А.С.

УДК 539.3

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

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

Ключевые слова: сопряженные процессы, аэрогазодинамика, термомеханика, гиперзвук, тепломассоперенос, термодеструкция, композиты, композиционные материалы, тепловая деформация, поровое давление, термонапряжения, расслоение.

Введение. Вопросам исследования гиперзвуковых аэродинамических режимов полета летательных аппаратов посвящено значительное число работ (например, [1-13] и многие другие), менее исследованы вопросы теплообмена при гиперзвуковых скоростях. Более сложную и значительно менее изученную проблему представляет собой высокотемпературное термомеханическое поведение композиционных материалов на основе термостойких матриц и наполнителей; значительный вклад в исследование этих проблем внесли работы [14-21]. Комплексные сопряженные задачи аэротермодинамики, теплообмена, теплофизики и термопрочности конструкций гиперзвуковых летательных аппаратов (ГЛА) также относительно мало исследованы. Вместе с тем в реальных условиях эксплуатации ГЛА задачи аэротермодинамики, теплообмена и теплофизики конструкций являются связанными через граничные условия на поверхности конструкции, поэтому параметры теплового потока, воздействующего на материалы, зависят от свойств этих материалов. В свою очередь, теплофизические свойства материалов при высоких температурах могут зависеть от напряженно-деформированного состояния конструкций: так, значительный уровень

термонапряжений в композиционных материалах приводит к микрорастрескиванию их матрицы еще задолго до полного макроразрушения конструкции, вследствие чего меняются газопроницаемость и теплопроводность материалов, а следовательно, и температурное поле в конструкции. Таким образом, для исследования реальных процессов, происходящих в конструкциях ГЛА, возникает необходимость разработки методов решения сопряженной задачи аэротермодинамики, теплообмена, теплофизики и термомеханики конструкций. Один из подходов к решению этой сопряженной задачи предложен в настоящей работе. Этот подход является дальнейшим развитием методов, разработанных в [21].

Общая система уравнений сопряженной задачи аэротермодинамики и термомеханики. Рассмотрим процесс обтекания конструкции ГЛА, имеющей конусообразную неосесимметричную форму, высокоскоростным набегающим потоком. Общая постановка сопряженной задачи аэротермодинамики и термомеханики состоит из трех систем уравнений:

• уравнений Навье - Стокса внешнего газового потока, обтекающего конструкцию;

• уравнений внутреннего тепломассопереноса в конструкции;

• уравнений термоупругости оболочечной конструкции.

Химические реакции во внешнем газовом потоке учитывать не

будем, так как эти реакции становятся существенными при скоростях потока с числом М > 9...10. В данной работе рассмотрены скорости ГЛА, не превышающие этих значений.

Система уравнений аэротермодинамики. Рассмотрим систему уравнений вязкого теплопроводного газа (уравнения Навье - Стокса), состоящую из уравнения неразрывности, уравнений движения и уравнения энергии. Записанная в бескоординатной форме, эта система в рассматриваемой области V движения потока газа имеет следующий вид [11]:

К этим уравнениям присоединяются определяющие соотношения вязкого совершенного газа:

(1)

I |2

р = Яр9; е = ^9; Е = е + у /2;

(2)

Т = д1(У-у)Е + д2(У®у + У®уг ); (3)

Я = -А,У9. (4)

Здесь р — плотность газа; I — время; Е — полная энергия газа,

I |2 V

Е = ^9 + ^^; ^ — теплоёмкость при постоянном объеме; 0 — тем-I |2 I

пература газа; у = vvi — квадрат модуля скорости; р — давление; Я — газовая постоянная (Я= М /д, ц — молекулярная масса газа, М — универсальная газовая постоянная); Е — метрический тензор; IV — тензор вязких напряжений в газе; я — вектор потока теплоты; ^1, ^2 — коэффициенты вязкости газа; X — коэффициент теплопроводности газа; V — набла-оператор Гамильтона [22]. Коэффициенты вязкости и теплопроводности газа являются функциями температуры, зависимости д1(9), д2(9) и Х(9) для воздуха выбирались согласно классической модели, которая описана, например в [23].

Рассмотрим четыре случая граничных условий для системы уравнений (1)-(4). На твердой непроницаемой поверхности обтекаемого тела к системе (1) присоединяются граничное условие прилипания и условие теплового баланса:

у = 0; -^9-п = ^е +вза94 -в,, (5)

где qe — тепловой поток, который идет на нагрев конструкции; 9^ — температура твердой стенки (совпадает с температурой газа на этой стенке); 9е — температура внешней поверхности пограничного слоя;

в ё, в, — интегральные коэффициенты излучения нагретого газа и

твердой поверхности; а — коэффициент Стефана - Больцмана. Унос материала конструкции не учитывался.

На границе входа потока Е2 (сверхзвуковой и дозвуковой) задаются условия:

Р = Ре, у = у е; 9 = 9е, (6)

где ре, уе, 9е — заданные значения параметров.

На дозвуковой границе выхода потока Е3 задается одно условие,

и еще четыре условия формулируются при численной аппроксимации решения [11]:

5у 59

р = ре; — = 0; — = 0, (7)

оп оп

дv

где — = п • V ® V — нормальная производная вектора скорости. На

дп

сверхзвуковой границе выхода потока Ед граничные условия не задаются, но при численной реализации формулируются следующиее условия:

дV = 0; » = 0. (8)

дп дп

На плоскости симметрии Е5 (которая может и отсутствовать) задаются условия симметрии:

дР= 0; V • п = 0; ^ = 0, I = 1,2; д®= 0. (9)

дп дп дп

Начальные условия для системы (1)-(5) имеют вид:

г = 0: р = р°; V = 0;Е = е¥%, (10)

где р0,90. — заданные значения.

Система уравнений внутреннего тепломассопереноса в конструкции ГЛА. Будем рассматривать конструкцию ГЛА, обтекаемую внешним потоком, изготовленную из композиционного материала на полимерной матрице с термостойкими керамическими волокнами. В матрице такого композита при нагреве до высоких температур, характерных для аэродинамического нагрева, происходят физико-химические процессы термодеструкции, сопровождаемые образованием газообразных продуктов терморазложения, которые накапливаются в порах материала и отфильтровываются во внешний газовый поток, а также формированием новой твердой фазы — пиролитической фазы матрицы, которая обладает существенно более низкими упруго-прочностными характеристиками, чем исходная полимерная фаза. В работах [14-17] была предложена четырехфазная модель для описания внутреннего тепломассопереноса и деформирования такого композита. Эта модель состоит из:

уравнения изменения массы полимерной фазы матрицы

ръ % = -/; х ЕК; (11)

дг

уравнения фильтрации газообразных продуктов термодеструкции в порах композиционного материала

др „ ф „

-т^^фяVя = /Г; х е¥; (12)

уравнения теплопереноса в термодеструктирующем композите

09

рс — = -V- я - - 3Ае0; х е V; (13)

Здесь и далее приняты обозначения: ф/, фь, ф и ф — объемные

концентрации трех твердых фаз (армирующие волокна, фаза исходной полимерной матрицы, пиролитическая фаза матрицы) и газовой фазы; р/, рь, р — плотности твердых фаз (полагаются постоянными);

р г — среднее по поре значение плотности газовой фазы (переменная

величина); с% — удельная теплоемкость газовой фазы при постоянном

объеме; р, с — плотность и удельная теплоемкость композита в целом; Я — вектор плотности теплового потока; 9 — температура композита, общая для всех фаз; уё — вектор скорости движения газовой фазы в порах; Ае0 — удельная теплота термодеструкции матрицы; 3 — массовая скорость термодеструкции матрицы; Г — коэффициент газификации матрицы.

К уравнениям (11)—(13) присоединяются определяющие соотношения, связывающие вектор-функции я , у г с функциями V9, Vp

с помощью законов Фурье и Дарси, а также соотношение Аррениуса для /-массовой скорости термодеструкции:

Я = -к -V 9; (14)

р г фг У г =-К ^ Р'; (15)

' Е ^

3 = 30 ехр

л

Я9у

(16)

где р — поровое давление §-фазы, для которого выполняется уравнение Менделеева — Клапейрона р = Ярг9 ; 30 — предэкспоненциаль-

ный множитель; Ел — энергия активации процесса термодеструкции; / — универсальная газовая постоянная; к — тензор теплопроводности; К — тензор газопроницаемости композита; зависят от концентраций фаз [14].

Кроме того, к системе (11—13) присоединяются соотношения для плотности и удельной теплоемкости композита:

р=р/ф/+рь фь + рр фр+рвф%;

ср = с / р/ф/ +сьрь фь +сррр фр +сг рвф%, (17)

где с /, сь, с — удельные теплоемкости твердых фаз при постоянной

деформации, полагаемые постоянными, не зависящими от температуры.

Введем также обозначение для плотности и удельной теплоемкости каркаса композита (совокупности всех твердых фаз):

Ps =Р/Ф/ +рьфь +РрФp; ^ =—(сfРfФ/ + сьрьфь + срРрфp). (18)

Ps

Объемная концентрация Фp пиролитической фазы матрицы может быть выражена аналитически через фь : Фp = ^Ф° - фь ) (1 -Г)—.

Рр

На нагреваемой части поверхности Ец конструкции граничные условия для уравнений (11-13) выглядят следующим образом:

X еЕц : p = pe; 0 = 0,, (19)

где ре — местное давление внешнего газового потока на поверхности Ец; 0, — температура поверхности композита.

На остальной части Е, поверхности Е композита задаем граничные условия герметичности и теплоизоляции:

хеЕ, : п-Ур = 0; -п• к-У0 = 0. (20)

Имеют место следующие соотношения между различными частями

поверхности тела: Е = Е, и Е = Еа и Еи; Еа =Е^ и Еq; = и Еи.

Для системы уравнений (11-13) начальные условия записывают в виде

г = 0: Фf =Ф0; Фь = Ф°°; ря =р°8; 0 = 0о. (21)

Система уравнений термомеханики оболочечной конструкции ГЛА. В криволинейной системе координат Оцц2Цз, связанной со срединной поверхностью оболочечной конструкции, система уравнений термомеханики ГЛА имеет следующий вид [14, 16, 18-20]: уравнения равновесия оболочки:

д д 8Лр дЛ дРд -(Л Т ) +-(Л Т )--^ Тпп^-^ ТаР + ЛаЛЛО,-= 0;

дц Р аа дq аР дц РР дц ар а^а ^ дц

8 ^ д _ , , ч дЛ

(ЛМ )+ — (Л М )--^ МРР +

дq Р аа дц а ар дq РР

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

а р а

дЛ 8Ря

+ ^ МаР - ЛаЛР0а - ЛР = 0;

д аР а Рдц

1Р 1 а

(22)

- Л Л2{кхтп + к2т22)+ддцО-- РеЛ1 Л2 - (/1 + к2) Л Л2^ Р8 = 0,

дЦ1 дц2

а, Р = 1,2, а^Р;

кинематические соотношения:

1 Зи„ 1 5Л тт , ттг „ 1 ОЖ ,

еаа =--а +--а ив+ каЖ; 2еа3 =--+ Уа- каиа;

Л 0q Л Л 0q р а а3 Л 0q а а'

Ла Л1Л2 Ла 0qа

1 ои, 1 0и2 1 ,04 ^ 0Л2 тт , 2е12 =--1 +--2--(-^ и1 + —2 и2); (23)

Л2 0q2 Л1 Л1Л2 0q2

1 5уа 1 0Ла Каа =--^ + ; 2Ка3 = -каУа;

аа

Ла ^а Л1Л2 ^р

1 Оу. 1 0у2 1 ,зд 0Л2 . 2к12 =--^ +--^--(—^ у1 + —2 у 2);

Л2 0q2 Л1 Л1Л2 0q2

определяющие соотношения оболочки:

2 ° Таа = Р^1(^арерр + ^рКрр) - Рга- Та;

2 ° Маа = Е (^арерр + АхрКрр ) -М%а -Ма; (24)

Т12 = 2(С66е12 + ^66К12); М12 = 2(^66е12 + А>6К12); ба = Са+3, а+3eа3,

а = 1, 2.

Здесь Таа, Тар, Маа, Мар — усилия и моменты в оболочке; ба — перерезывающие усилия; еаа, еа3, е12 — деформации срединной поверхности оболочки; каа, ка3, к12 — искривления срединной поверхности; иа, уа ,Ж — перемещения, углы искривления и прогиб срединной поверхности; Ла, ка — параметры первой квадратичной формы и главные кривизны срединной поверхности оболочки [22]; Рг, Мг — усилие и момент порового давления в оболочке,

и/ и/

/2 /2

Рг = I ф&р^Ъ; Мг = I фgРqзdqз. (25)

- И/ -И/

/2 /2

Введены также обозначения для усилий и моментов тепловых

о о

напряжений Та,Ма, зависящих от температурных деформаций Ва оболочки:

о 3 0 о 3 о

Та = Е Сар вр ; Ма = Е СаР вр ; (26)

Р=1 р=1

о (Л) % о о(л) И2 о

вр = I «91 врqз,dq3; вэ = | взq3dq3, л=0, 1, р= 1, 2; (27)

_И/ -И/

/2 /2

о t О

S7= (a f ф fB1+ab фй Qy )(0 - 0О) + а р QyJ (0(t) - 0(т)) ф „dx - р „ ф „ Qy,

О

У = 1,2,3. (28)

Здесь а f, аь, а р — коэффициенты теплового расширения волокна, полимера и пиролитической фазы матрицы; р — коэффициент усадки; By, Qy — коэффициенты, зависящие от расположения волокон в композите [14].

Усилия и моменты межфазного взаимодействия Pga,Mga в оболочке определены следующим образом:

V V

/2 /2

Pga = J Р fad43; Mga = J Р faЧз^ (29)

-А/

/2 /2

где fa — коэффициент межфазного взаимодействия [14].

В соотношениях (24) обозначены мембранные, смешанные и из-гибные жесткости оболочки Сар, Nap, Dap :

С = C0 a(0)- N = C0 Л1)- D = C0 a(2)- (30)

Сар _ CaPa01 ; CaPa01 ; Da^ CaPa01; (30)

с = C 0 a(0)- N = C 0 aa). D = C 0 a (2)-C66 - C66a01 ; N66 " C66a01 ; D66 " C66a01 ;

h/

_ 0 (0) ( ') '

Са+3,a+3 = Ca+3,a+3a02 , a = 1, 2; a0k = J a0kq3dq3, к=1, 2;'=0, 1 2;

- h/ /2

Вследствие размягчения полимерной матрицы и её термодеструкции жесткости оболочки изменяются при нагреве, учет этого изменения для ортотропных композитных оболочек осуществляется с помощью двух функций: a^, [14-17].

Система уравнений равновесия (22), в которую подставлены определяющие соотношения (24) и кинематические соотношения (23), образует замкнутую систему из пяти уравнений относительно пяти неизвестных функций Ua, y a, W.

В качестве граничных условий, например на линии qa = const, к этой системе следует присоединить заданные значения пяти величин (по одному из каждой пары): (T11 ^gPg, ua ), (T22, Up ), (Qa ,W),

(M11 -фgMg, Уа ), (M22, Ур ).

Деформации вар и напряжения аар в оболочке вычисляются по следующим формулам:

8аР - еаР + #3КаР, а, Р-^ 2; 833 - 0; 8аз - ва3;

-ар

_ 3

°аа - ~/аР + а61 Е СаР р-1

о Л

8РР+^3КРР 8р

а -1, 2;

(31)

(32)

012 - а61С66(812 + ^3К12)-

Для поперечного нормального напряжения 033 и напряжений межслойного сдвига оа3 имеем следующие формулы [14, 15]:

°33 - 6Л

+ ]^С31(ае0)е11 + <кп -8(0)) + 2 п п

1 0 (0) (1) о(0) 1 0 0(0)"

+ — С32 (ае1 е22 + ае1 К22 - 82 ) -~Т С33 83

п п

+(Р2 - Р1 ) + Р1 2 Р2 + Ф

013 -

12Л(^3)С0 е а(0). . _ - 12Л(^3)С0 е а(0).

п

'55с23"62 >

(33)

^3) - 1 - ^ лЫ - "4 - (Т)2; а-1,2,

2 п 4 п

где Р1, Р2 — давления на внешних поверхностях оболочки. Максимальные значения касательных напряжений достигаются на срединной поверхности, где ^(0) - 4 :

(34)

0 -3С е а(0)- а -3С е а(°) 13тах _ , С44е13а62 ; а23шах " , С55е23а62 .

пп

Метод численного решения сопряженной задачи. Для численного решения сопряженной задачи использовался пошаговый метод, с двумя шагами по времени Д^, Д^. Максимальный крупный шаг

Д^1 по времени использовался для изменения граничных условий задачи — входных данных набегающего потока и решения задачи внутреннего тепломассопереноса. На каждом шаге Д^ решение осуществлялось в четыре этапа.

Этап 1. Решение системы уравнений внутреннего тепломассопереноса (11)-(20). Это решение осуществлялось численным конечно-разностным методом с использованием метода линеаризации и неявной разностной схемы, детали которого описаны в [11, 12]. Температура 6^ поверхности конструкции, взаимодействующей с набегаю-

щим газовым потоком, на этом этапе полагалась известной и бралась с предыдущего временного шага по Л^.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Этап 2. Осуществлялся цикл решения системы уравнений Навье -Стокса (1)-(4) с граничными условиями (5)-(9). После каждого приращения на Л^ осуществлялся численный расчет системы уравнений Навье - Стокса с граничными (5) и начальными (1)-(8) данными пошаговым методом с мелким шагом Л^, до установления потока. Численный расчет осуществлялся методом трехмерного пограничного слоя [10] с расщеплением по физическим переменным, при котором сначала производился расчет уравнений идеального газа методом ленточных адаптивных сеток с применением схемы ТУО, а затем — решение уравнений вязкого газа с применением факторизованных неявных разностных и методов матричной прогонки.

В граничном условии теплообмена (5) при численном расчете тепловой поток на твердой стенке = X5 • п , вообще говоря, не-

7 /Т7 \

известный, вычислялся по формуле = а(0^ - 0О), где а = —--

Н

коэффициент теплообмена; Н — толщина стенки материала; Х3 — коэффициент теплопроводности стенки оболочки в поперечном направлении; g (Бо) — безразмерный коэффициент теплообмена, за-

^ Т-. Xзto

висящий от параметра Фурье конструкции материала Бо = —; ^ —

РЛН

характерное «медленное» время задачи. Функция g(Fo) рассчитывалась заранее, с помощью предварительного численного решения одномерной задачи теплопроводности для пластины при различных значениях параметра Бо. С учетом такой аппроксимации теплового потока граничное условие (5) было выбрано в следующем виде:

Е„ : ХУ0-п = а(0,-0о) + 8^-Sgс04. (35)

Далее это условие присоединялось к системе (1)-(4), в результате совместно с вычислением параметров V, р, 0, р в области внешнего газового потока, обтекающего конструкцию, вычислялась температура стенки 0К.

Этап 3. Осуществлялось решение системы уравнений термоупругости оболочечной конструкции (22)-(24) с помощью метода конечного элемента, подробности этого метода решения задачи изложены в [18—21].

Этап 4. После решения задачи термоупругости осуществлялся расчет термонапряжений в оболочке с помощью формул (32), (33). Далее осуществлялся переход на следующий шаг по времени Л^.

Результаты численного моделирования. На рис. 1-4 представлены результаты численного моделирования обтекания фрагмента корпуса модельного летательного аппарата гиперзвуковым потоком газа на высоте 15 км. Начальные условия и скорость набегающего потока (граничное условие на входной границе) были заданы равными и имели вид: р-0,195 кг/м ; ух - уу - 0 м/с; -1800 м/с; Р -12346 Па. Значения размерных параметров задачи были приняты следующими: 60 - 293К; 8! - 0,8 и 82 - 0,3; Xя - 0,3Вт/(м • К); р^ -1800 кг/м3; - 0,8 кДж/(кг• К).

3

Рис. 1. Распределение плотности газового потока р кг/м , набегающего на конструкцию ГЛА в момент времени t = 0

Рис. 2. Распределение давления газового потока P, Па, набегающего на конструкцию ГЛА в момент времени t = 0

Рис. 3. Распределение температуры газового потока 0, К, набегающего на конструкцию ГЛА в момент времени t = О

Рис. 4. Распределение продольной скорости газового потока V, м/с, набегающего на конструкцию ГЛА в момент времени t = О

Численное решение с хорошим разрешением воспроизводит картину обтекания фрагмента корпуса гиперзвуковым потоком: в критической точке ЛА формируется скачок уплотнения, а в окрестности поверхности конструкции образуется ударная волна. Максимум плотности, давления и температуры приходится на критическую точку носка аппарата, в которой температура потока достигает 2000 К.

На рис. 5-7 представлены результаты численных расчетов полей внутреннего тепломасспереноса в оболочке конструкции ГЛА на нижней его поверхности на удалении 30 Я от критической (контрольной) точки, где Я - радиус затупления ГЛА в критической точке. На

рис. 5 показано распределение температуры по толщине оболочки для семи различных моментов времени 1Ъ..., Ц полета ГЛА, где

х = 0,5 + qзlh — безразмерная координата по толщине оболочки.

0,6 0,5 0,4 0,3 0,2 ОД

- 4

0,2 0,4 0,6 0,8 1,0

0 0,2 0,4 0,6 0,8 1,0

Рис. 5. Распределение перепада температуры 9-90, °С, по безразмерной координате х в контрольной точке оболочки в различные моменты времени полета ГЛА

Рис. 6. Распределение концентрации полимерной фазы фь

композита по безразмерной координате х оболочки для различных моментов времени полета ГЛА

На рис. 6 показано распределение концентрации полимерной фазы по толщине оболочки для тех же семи моментов времени. На рис. 7 представлено распределение концентрации пиролитической фазы по толщине оболочки для различных моментов времени.

0,35 0,30 0,25 0,20 0,15 0,10 0,05

Ч

и

Н

, 8.347е+002 7 754е+002 И 7.162е+002 > И 6 570а+002 М Щ 5.977е+002 В И 5.385е+002 ^^

0 0,2 0,4 0,6 0,8 1,0

Рис. 7. Распределение концентрации пиролитической фазы

Фь композита по безразмерной

координате х в контрольной точке оболочки в различные моменты времени полета ГЛА

Рис. 8. Распределение температуры, К, поверхности оболочки ГЛА в момент времени / = 17

На рис. 8 показано распределение температуры по поверхности оболочки в момент времени ^ полета.

Терморазложение полимерной фазы композита приводит к образованию большого количества газообразных продуктов в порах композиционного полимерного материала оболочки корпуса. Ввиду низкой газопроницаемости композита образующиеся газы не успевают отфильтровываться во внешний газовый поток и создают внутреннее поровое давление. Распределение внутрипорового давления в оболочке показано на рис. 9 и 10 в различные моменты времени.

Рис. 9. Распределение максимального значения порового давления рт, ГПа, в момент времени t = t7

Рис. 10. Распределение порового давления р, атм, газообразных продуктов терморазложения полимерной фазы ф композита по

безразмерной координате х в контрольной точке оболочки в различные моменты времени

На рис. 11 показано распределение коэффициента изменения упругих свойств композита «01 по толщине оболочки для разных моментов времени, на рис. 12 — распределение значений тепловой про-

о

дольной деформации 81 по толщине оболочки для четырех контрольных точек в различные моменты времени. Характерной особенностью этих функций является их немонотонность: в зоне повышенной температуры (по толщине оболочки) происходит рост зна-

чений тепловой деформации в!, однако при переходе через условную точку начала интенсивного терморазложения происходит уменьшение тепловой деформации, обусловленное образованием пиролитиче-ского остатка. Это уменьшение сказывается на значениях тепловых напряжений оп, о22 в оболочке.

Рис. 11. Распределение коэффициента изменения упругих свойств

композита а по безразмерной

координате х в контрольной точке оболочки в различные моменты времени полета ГЛА

Рис.

0,2 0,4 0,6 0,8 1,0 12. Распределение

тепловой деформации в1

по безразмерной координате х в контрольной точке В оболочки в различные моменты времени полета ГЛА

На рис. 13 показано распределение значений тепловой поперечной деформа-

о

ции 83 по толщине оболочки для четырех

контрольных точек в различные моменты

времени. Значения этой деформации пре-

о

восходят значения 81 , поскольку матрица

обладает значительно большим коэффициентом теплового расширения, чем во- ' о 0,2 0,4 0,6 0,8 1,0

о

локна. Так же как и функция в1, деформа- Рис. 13. Распределение тепло-

О

о

ция в3 немонотонно распределена по вой деформации 83 по без-толщине оболочки: в зоне повышенной размерной к^рдтате х т°л-

температуры (по толщине оболочки) про- щины в контрольной точке

о оболочки в различные момен-

исходит рост значений 83 , при переходе ты времени полета ГЛА

через условную точку начала интенсивного терморазложения проис-

о

ходит уменьшение значений 83 , обусловленное образованием пиро-

литического остатка, причем, в отличие от в!, тепловая деформация

о

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

В процессе нагрева сжимающее окружное напряжение 022 (И / 2) во всех контрольных точках постепенно увеличивает свои значения, вместе с этим увеличиваются и значения максимальных растягивающих значений окружного напряжения на периферийной части оболочки ближе к кромкам. В момент времени Ш возникает пик положительных растягивающих напряжений 022 (И / 2), обусловленный термодеструкцией композита, вследствие которой возрастает поровое давление газообразных продуктов терморазложения матрицы (рис. 14).

Рис. 14. Распределение максимального окружного напря-женийя ст22 (И/2), ГПа, на внутренней поверхности оболочки в момент времени t = /7

Рис. 15. Распределение максимального поперечного напряжения О33, ГПа, в оболочке в момент времени t = ^

На рис. 15, 16 показано распределение поперечного напряжения 033 по

толщине оболочки в различные моменты времени полета ГЛА. Графики показывают, что в подповерхностном слое после начала терморазложения полимера возникает значительный максимум напряжения 033, обусловленный ростом поро-

вого давления газообразных продуктов термодеструкции.

Из этих графиков видно, что в зоне наибольшего нагрева происходит значительное возрастание положительных (растягивающих) значений поперечных напряжений 033. Причиной такого роста

поперечных напряжений является, главным образом, рост порового давления газообразных продуктов термодеструкции в подповерхностном слое оболочки, а также возникновение усадочных деформаций. К моменту времени ^ значения поперечных напряжений на нижней части оболочки достигают значений в 0,13 ГПа, что значительно превышает предел прочности композитной оболочки в поперечном направлении. В результате, в этой части оболочки может возникнуть разрушение по типу расслоения, при котором верхние слои ткани композита отслоятся от остальной части материала. Аналогичное явление может произойти и на верхней части оболочки к моменту т поскольку поперечные напряжения на верхней части оболочки начинают расти и достигают значений в 0,06 ГПа, что хотя и меньше значений в нижней части оболочки, но также сравнимо с пределом прочности материала в поперечном направлении.

Выводы. Предложенный метод численного моделирования позволяет осуществлять расчеты сопряженных аэрогазодинамических и термомеханических процессов в композитных конструкциях высокоскоростных летательных аппаратов с учетом процессов термодеструкции в полимерных композитах, процессов внутреннего тепло-массопереноса, изменения упругих характеристик при нагреве. Проведенные численные расчеты сопряженных процессов для модельной конструкции ГЛА показали, что наиболее вероятным механизмом нарушения работоспособности полимерной композитной оболочки ГЛА является расслоение его конструкции. Причиной вероятного расслоения является терморазложение полимерной матрицы композиционного материала, происходящее при температуре 400.. .600 оС, а также низкая пористость материала и, как следствие,

0 0,2 0,4 0,6 0,8 1,0

Рис. 16. Распределение поперечного напряжения О33 по

безразмерной координате х в контрольной точке оболочки в различные моменты времени

низкая проницаемость по отношению к газообразным продуктам терморазложения, фильтрующимся по порам.

Исследование выполнено за счет средств Задания № 1.445.2014/К Ми-нобрнауки РФ, гранта Российского научного фонда (проект №14-19-00847) и гранта РФФИ № 12-08-000998-а.

ЛИТЕРАТУРА

[1] Anderson J. D. Hypersonic and high-temperature gas dynamics. 2th ed. American Institute of Aeronautics and Astronautics, Reston, Virginia, 2006, 232 p.

[2] Лунёв В.В. Гиперзвуковая аэродинамика. Москва, Машиностроение, 1975, 330 с.

[3] Тирский Г.А. Гиперзвуковая аэродинамика и тепломассообмен спускаемых космических аппаратов и планетных зондов. Москва, Физматлит, 2011, 548 с.

[4] Лесин А.Б., Лунёв В.В. Аномальный теплообмен на треугольной пластине с затупленным носком в гиперзвуковом потоке. Механика жидкости и газа, 1994, № 2.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

[5] McNamara J., Friedmann P. Aeroelastic and aerothermoelastic analysis of hypersonic vehicles: Current status and future trends. 48th AI-AA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 23-26 April 2007, Honolulu, Hawaii. URL: http://www.mecheng.osu.edu/lab/cael/sites/default/files/AIAA-2007-2013

[6] Crowell A.R., McNamara J.J., Miller B.A. Hypersonic aerothermoelastic response prediction of skin panels using computational fluid dynamic surrogates. ASDJournal, 2011, vol. 2, no. 2, pp. 3-30.

[7] В.П. Котенев, В.А. Сысенко. Аналитические формулы повышенной точности для расчета распределения давления на поверхности выпуклых затупленных тел вращения произвольного очертания. Математическое моделирование и численные методы, 2014, № 1, с. 68-81.

[8] Братчев А.В., Забарко Д.А., Ватолина Е.Г., Коробков А.А., Сахаров В.И. Вопросы теплотехнического проектирования перспективных гиперзвуковых летательных аппаратов аэробаллистического типа. Изв. ин-та инженерной физики, 2009, т. 2, № 12, с .42-49.

[9] Полежаев Ю.В., Юревич Ф.Б. Тепловая защита. Москва, Энергия, 1976, 368 с.

[10] Димитриенко Ю.И., Захаров А.А., Коряков М.Н. Модель трехмерного пограничного слоя и ее численный анализ. ВестникМГТУ им. Н.Э.Баумана. Сер. Естественные науки, 2011, отец. выпуск, с. 136-150.

[11] Димитриенко Ю.И., Котенев В.П., Захаров А.А. Метод ленточных адаптивных сеток для численного моделирования в газовой динамике. Москва, Физматлит, 2011, 286 с.

[12] Димитриенко Ю.И., Коряков М.Н., Захаров А.А., Сыздыков Е.К. Развитие метода ленточно-адаптивных сеток на основе схем TVD для решения задач газовой динамики. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2011, № 2, с. 87-97.

[13] Гильманов А.Н. Методы адаптивных сеток в задачах газовой динамики. Москва, Физматлит, 2000, 248 с.

[14] Димитриенко Ю.И. Механика композиционных материалов при высоких температурах. Москва, Машиностроение, 1997, 366 с.

[15] Dimitrienko Yu.I. Thermal stresses and heat mass-transfer in ablating composite materials. International Journal of Heat Mass Transfer, 1995, vol. 38, no. 1, pp. 139-146.

[16] Dimitrienko Yu.I. Thermal stresses in ablative composite thin-walled structures under intensive heat flows. International Journal of Engineering Science, 1997, vol. 35, no. 1, pp. 15-31.

[17] Dimitrienko Yu.I. A structural thermomechanical model of textile composite materials at high temperatures. Composite science and technologies, 1999, vol. 59, pp. 1041-1053.

[18] Димитриенко Ю.И., Минин В.В., Сыздыков Е.К. Моделирование термомеханических процессов в композитных оболочках при локальном нагреве излучением. Механика композиционных материалов и конструкций,

2011, т. 17, №1, с. 71-91.

[19] Димитриенко Ю.И., Минин В.В., Сыздыков Е.К. Моделирование внутреннего тепломассопереноса и термонапряжений в композитных оболочках при локальном нагреве. Математическое моделирование, 2011, т. 23, № 9, с.14-32.

[20] Димитриенко Ю.И., Минин В.В., Сыздыков Е.К. Численное моделирование процессов тепломассопереноса и кинетики напряжений в термоде-структирующих композитных оболочках. Вычислительные технологии,

2012, т. 17, № 2, с. 44-60.

[21] Димитриенко Ю.И., Захаров А.А., Коряков М.Н., Сыздыков Е.К. Численное решение сопряженной задачи аэрогазодинамики и внутреннего теп-лопереноса в когатрукциях гиперзвуковых летательных аппаратов. Инженерный журнал: наука и инновации, 2012, вып. 11. URL: http://engjournal.ru/articles/426/426.pdf

[22] Димитриенко Ю.И. Тензорное исчисление. Москва, Высшая школа, 2001, 575 с.

[23] Краснов Н.Ф. Аэродинамика. В 2 т. Москва, Высшая школа, 1980, т. 1, 495 с., т. 2, 416 с.

Статья поступила 05.11.2014

Ссылку на эту статью просим оформлять следующим образом:

Димитриенко Ю.И., Коряков М.Н., Захаров А.А., Строганов А.С. Численное моделирование сопряженных аэрогазодинамических и термомеханических процессов в композитных конструкциях высокоскоростных летательных аппаратов. Математическое моделирование и численные методы, 2014, № 3, с. 3-24.

Димитриенко Юрий Иванович родился в 1962 г., окончил МГУ им. М.В. Ломоносова. Д-р физ.-мат наук, профессор, зав. кафедрой «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана, директор Научно-образовательного центра «Суперкомпьютерное инженерное моделирование и разработка программных комплексов» МГТУ им. Н.Э. Баумана. e-mail: dimit.bmstu@gmail.com

Коряков Михаил Николаевич родился в 1987 г., младший научный сотрудник Научно-образовательного центра «Суперкомпьютерное инженерное моделирование и разработка программных комплексов» МГТУ им. Н.Э. Баумана, ассистент кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. e-mail: mkoryakov@bmstu.ru

Захаров Андрей Алексеевич родился в 1982 г., окончил МГТУ им. Н.Э. Баумана. Канд. физ.-мат. наук, доцент кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана, старший научный сотрудник Научно-образовательного центра «Суперкомпьютерное инженерное моделирование и разработка программных комплексов» МГТУ им.Н.Э. Баумана. e-mail: azaharov@bmstu.ru

Строганов Александр Сергеевич родился в 1990 г., окончил МГТУ им. Н.Э. Баумана, младший научный сотрудник Научно-образовательного центра «Суперкомпьютерное инженерное моделирование и разработка программных комплексов» МГТУ им. Н.Э. Баумана. e-mail: alt90@rambler.ru

Computational modeling of conjugated gasdynamic and thermomechanical processes in composite structures of high speed aircraft

In the article we propose an algorithm for the numerical simulation of conjugate gasdynamic and thermomechanical processes in composite structures of high-speed aircraft. The algorithm allows calculating all parameters of the three-dimensional gasdynamic flow near the surface of the aircraft, heat exchange on the surface, heat and mass transfer processes in the internal structure of thermodestructive polymer composite, as well as processes of composite construction thermodeformation, including the effects of changes in the elastic characteristics of the composite, variable thermal deformation, shrinkage caused by thermal degradation, building up interstitial gas pressure in the composite. An example of numerical simulation of conjugated processes in a model composite construction of high-speed aircraft illustrates the possibilities of the proposed algorithm.

Keywords: conjugated processes, aerogasdynamics, thermomechanics, hypersonic flows, heat and mass transfer, thermodestruction, composites, composite materials, thermal deformation, pore pressure, thermal strains, bundle.

REFERENCES

[1] Anderson J.D. Hypersonic and high-temperature gas dynamics. 2nd ed. American Institute of Aeronautics and Astronautics, Reston, Virginia, 2006, 232 p.

[2] Lunev V.V. Giperzvukovaya aerodinamika [Hypersonic gasdynamics]. Moscow, Mashinostroenie Publ., 1975, 330 p.

[3] Tirskiy G.A., et. al. Giperzvukovaya aerodinamika i teplomassoobmen spus-kaemykh kosmicheskikh apparatov i planetnykh zondov [Hypersonic gasdynam-ics and heat and mass transfer in reentry spacecrafts and planetary probes]. Moscow, Fizmatlit Publ., 2011, 548 p.

[4] Lesin A.B., Lunev V.V. Mekhanika zhidkosti i gaza - Fluid Mechanics, 1994, no. 2.

[5] McNamara J., Friedmann P. Aeroelastic and aerothermoelastic analysis of hypersonic vehicles: Current status and future trends. 48th AI-AA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 23-26 April 2007, Honolulu, Hawaii. Available at: http://www.mecheng.osu.edu/lab/cael/sites/default/files/AIAA-2007-2013

[6] Crowell A.R., McNamara J.J., Miller B.A. Hypersonic aerothermoelastic response prediction of skin panels using computational fluid dynamic surrogates. ASDJournal, 2011, vol. 2, no. 2, pp. 3-30.

[7] Kotenev V.P., Sysenko V.A. Matematicheskoe modelirovanie i chislennye meno-dy - Mathematical Modeling and Numerical Methods, 2014, no. 1, pp. 68-81.

[8] Bratchev A.V., Zabarko D.A., Vatolina E.G., Korobkov A.A., Sakharov V.I. Izvestiya Instituta inzhenernoy fiziki - Proceedings of the Engineering Physics Institute, 2009, vol. 2, no. 12, pp .42-49.

[9] Polezhaev Yu.V., Yurevich F.B. Teplovaya zashchita [Thermal protection]. Moscow, Energiya Publ., 1976, 368 p.

[10] Dimitrienko Yu.I., Zakharov A.A., Koryakov, M.N. Vestnic MGTU im. N.E. Baumana. Seria Estestvennye nauki - Herald of the Bauman Moscow State Technical University. Series: Natural Sciences, 2011, special issue, pp. 136-150.

[11] Dimitrienko Yu.I., Kotenev V.P., Zakharov A.A. Metod lentochnykh adap-tivnykh setok dlya chislennogo modelirovaniya v gazovoy dinamike [The adaptive banded grid method for numerical simulation in gas dynamics]. Moscow, Fizmatlit Publ., 2011, 280 p.

[12] Dimitrienko Yu.I., Koryakov, M.N., Zakharov A.A., Syzdykov E.K. Vestnic MGTU im. N.E. Baumana. Seria Estestvennye nauki - Herald of the Bauman Moscow State Technical University. Series: Natural Sciences, 2011, no. 2, p. 87-97.

[13] Gilmanov A.N. Metody adaptivnykh setok v zadachakh gazovoy dinamiki [The adaptive grid methods in gas dynamics problems]. Moscow, Fizmatlit Publ., 2000, 248 p.

[14] Dimitrienko Yu.I., Mekhanika kompozitsionnykh materialov pri vysokikh tem-peraturakh [Composite material mechanics under high temperature]. Moscow, Mashinostroenie Publ., 1997, 366 p.

[15] Dimitrienko Yu.I. Thermal stresses and heat mass-transfer in ablating composite materials. International Journal of Heat Mass Transfer, 1995, vol. 38, no. 1, pp. 139-146.

[16] Dimitrienko Yu.I. Thermal stresses in ablative composite thin-walled structures under intensive heat flows. International Journal of Engineering Science, 1997, vol. 35, no. 1, pp. 15-31.

[17] Dimitrienko Yu.I. Composite science and technologies, 1999, vol. 59, pp. 1041-1053.

[18] Dimitrienko Yu.I., Minin V.V., Syzdykov E.K. Mekhanika kompozitsionnykh materialov i konstruktsiy - Mechanics of Composite Materials and Structures,

2011, vol. 17, no. 1, pp. 71-91.

[19] Dimitrienko Yu.I., Minin V.V., Syzdykov E.K. Matematicheskoe modelirovanie -MathematicalModeling. 2011, vol. 23, no. 9, pp.14-32

[20] Dimitrienko Yu.I., Minin V.V., Syzdykov E.K. Vychislitelnye tekhnologii -Computational Technologies, 2012, vol. 17, no. 2, pp. 44-60.

[21] Dimitrienko Yu.I., Zakharov A.A., Koriakov M.N., Syzdykov E.K. Inzhenernyi zhurnal: nauka i innovatsii — Engineering Journal: Science and Innovations,

2012, issue 11. Available at: http://engjournal.ru/articles/426/426.pdf

[22] Dimitrienko Yu.I. Tenzornoe ischislenie [Tensor calculus]. Moscow, Vysshaya Shkola Publ., 2001, 575 p.

[23] Krasnov N.F., Aerodinamika [Aerodynamics]. In 2 volumes. Moscow, Vysshaya Shkola Publ., 1980, vol.1, 495 p., vol. 2, 416 p.

Dimitrienko Yu.I. (b. 1962) graduated from Lomonosov Moscow State University in 1984. Dr. Sci. (Phys & Math.), professor, head of the Computational Mathematics and Mathematical Physics Department, director of Scientific-educational Center of Supercomputer Engineering Modeling and Program Software Development at Bauman Moscow State Technical University. Member of the Russian Academy of Engineering Science. Author of over 300 publications in the field of computational mechanics, gasdynamics, thermomechanics of composite materials, mathematical simulations in material science. e-mail: dimit.bmstu@gmail.com

Koryakov M.N. (b. 1987) junior member of teaching staff of the Computational Mathematics and Mathematical Physics Department, associate scientist of Scientific-educational Center of Supercomputer Engineering Modeling and Program Software Development at Bauman Moscow State Technical University. e-mail: mkoryakov@bmstu.ru

Zakharov A.A. (b. 1982) graduated from Bauman Moscow State Technical University. Ph.D., assoc. professor of the Computational Mathematics and Mathematical Physics Department, senior scientist of Scientific-educational Center of Supercomputer Engineering Modeling and Program Software Development at Bauman Moscow State Technical University. e-mail: azaharov@bmstu.ru

Stroganov A.S. (b. 1990) graduated from Bauman Moscow State Technical University, associate scientist of Scientific-educational Center of Supercomputer Engineering Modeling and Program Software Development at Bauman Moscow State Technical University. e-mail: alt90@rambler.ru

Аннотация научной статьи по механике и машиностроению, автор научной работы — Димитриенко Ю.И., Коряков М.Н., Захаров А.А., Строганов А.С.

Похожие темы научных работ по механике и машиностроению , автор научной работы — Димитриенко Ю.И., Коряков М.Н., Захаров А.А., Строганов А.С.