Конечно-элементное моделирование термонапряжений в композитных термодеструктирующих конструкциях при аэродинамическом нагреве Текст научной статьи по специальности «Физика»

УДК 539.3

Конечно-элементное моделирование термонапряжений в композитных термодеструктирующих конструкциях при аэродинамическом нагреве

© Ю.И. Димитриенко, М.Н. Коряков, Ю.В. Юрин, А.А. Захаров МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

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

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

Введение. В ракетно-космической технике для теплозащиты внешних теплонагруженных частей летательных аппаратов часто применяется пассивный метод защиты, который заключается в использовании аблирующих композитных материалов [1-4]. К преимуществам пассивного метода теплозащиты относятся: относительная низкая стоимость, простота использования, надёжность, отсутствие ограничений по максимальному тепловому потоку и большой выбор материалов в связи с развитием органической химии. Согласно введенной в [5] классификации, существуют 3 основных типа аблирующих материалов: 1) материалы с линейной абляцией, у которых при аэродинамическом нагреве происходит изменение геометрических размеров при неизменной плотности, 2) материалы с объемной абляцией, у которых при высоких температурах происходит изменение плотности, при неизменном объеме, 3) материалы, у которых линейная и объемная абляция реализуются одновременно. В данной работе рассматриваются материалы 2-го класса, к которому относятся, например, композиты на основе термодеструктирующих (термореактивных) полимерных

матриц. При умеренно высоких температурах (примерно до 2000 К) и умеренных скоростных напорах аэродинамического потока эти материалы терморазрушаются по модели 2-го типа.

В работах [5-7] была разработана математическая модель процессов внутреннего тепломассопереноса в полимерных термодеструктирующих композиционных материалах при нестационарном нагреве, учитывающая процессы термодеструкции, газообразования и фильтрации газов в порах материала. В работах [810] были разработаны методы решения задач расчета напряженно -деформированного состояния конструкций из термодеструктирующих композитов, основанные на обобщенной теории пластин и оболочек и конечно-элементном анализе этих уравнений. Поскольку в различные моменты движения летательного аппарата в атмосфере на его конструкцию воздействуют меняющиеся силовые и тепловые нагрузки, то для задания граничных условий для задачи термомеханики необходимо предварительно решить задачу аэродинамики, в результате возникает сопряженная задача аэротермомеханики. В работе [11] разработан метод решения неполной сопряженной задачи — аэрогазодинамики и внутреннего тепломассопереноса.

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

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

- уравнений Навье-Стокса для внешнего высокоскоростного потока;

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

- уравнений термоупругости композитной термодеструкти-рующей конструкции;

с соответствующими граничными и начальными условиями.

Рассмотрим систему уравнений высокоскоростного потока вязкого теплопроводного газа (уравнения Навье-Стокса) [12, 13, 14], обтекающего рассматриваемую композитную конструкцию:

^•ру = 0,

ы

дрУ + V•(ру ® У + pE - Ту) = 0, (1)

ы

дре

У-((ре + р) у - т • у + q) = 0,

+ V • ((ре + р) у - Т • у + q ) =

где обозначены следующие параметры высокоскоростного потока: р — плотность газа, у — вектор скорости, р — давление, Е — единичный тензор, е — плотность полной энергии газа, Ту — тензор вязких напряжений, q — вектор теплового потока, I — время, V — на-бла-оператор, ® — знак тензорного произведения [15].

Присоединим к этой системе (1) определяющие соотношения совершенного линейно-вязкого газа:

Я , ,2

Р = Р — °, е = св, е = е + у /2, М

Т =^1(V• у )Е + ^2(V® у + ^ ут ), (2)

q = ^в,

где Я — универсальная газовая постоянная, М — молекулярная масса, в — температура, е — внутренняя энергия, ср — теплоёмкость при постоянном объёме, / и /2 — коэффициенты вязкости, X — теплопроводность газа. Коэффициенты вязкости и теплопроводности являются функциями температуры и вычисляются по формулам [12, 13, 14]. Система уравнений рассматривается в области Р1, занятой газом.

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

у = 0, в = в„, (3)

где в^ — температура поверхности конструкции.

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

условия невозмущенного потока.

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

конструкции. Таким образом, исходный композиционный материал рассматривается как 4-х фазная структура [5-7].

Математическая модель внутреннего тепломассообмена в термодеструктирующем композите состоит из системы уравнений изменения массы полимерной фазы, уравнения фильтрации газообразных продуктов в порах материала и уравнения теплопроводности [5]

Рь -( = - 3,

ь дt

дР( *7 гт

ее + Ур( = 3Г, (4)

дt

дО

здесь введены обозначения: (ъ, ( — объёмные концентрации полимерной и газовой фаз композита; ръ — плотность исходной полимерной матрицы; — плотность газовой фазы в порах; ся — теплоёмкость газа в порах при постоянном объёме, р, и с, — плотность и теплоёмкость композита, ч, — вектор теплового потока в композите, О,, — температура композита; уя — скорость фильтрации газовой

фазы; Ле0 — удельная теплота термодеструкции матрицы; 3 — массовая скорость термодеструкции матрицы; Г — коэффициент газификации матрицы. Концентрацию пиролитической фазы композита вычисляем при помощи соотношения :

( = (((-(ъ)(1 - Г)Р.

Рр

Уравнения тепломассопереноса (4) замыкаются определяющими соотношениями, которые связывают вектор-функции ч, и уя с градиентами температуры У О и давления Ур при помощи законов

Фурье и Дарси, а также соотношением Аррениуса для массовой

скорости термодеструкции и уравнением Менделеева-Клапейрона:

К О

Ч = -Л • УО, (5)

Р( у Я =-К •УРЯ,

3 = 3о(ъ еХР (-Еа/К08 ) ,

где J0 — предэкспоненциальный множитель, ЕА — энергия активации процессов термодеструкции, — молекулярная масса газа, р

— поровое давление газов, Л — тензор теплопроводности, К — тензор газопроницаемости композита.

Для композита, как для многофазной системы, выполняются соотношения смеси и условие нормировки (индекс / относится к армирующему наполнителю, который полагаем термоустойчивым в рассматриваемом диапазоне температур, например, это могут быть углеродные волокна):

Р.*с. = Ргсг^г + РъсъРъ +Ррсрфр, Р. =РгЬ + РъРъ + Ррфр + РР, (6)

<Рг + Ръ + Рр + Р. = 1,

где рг — пористость.

Граничные условия для системы уравнений (4)-(5) на нагреваемой поверхности Е1 следующие:

р. = Р. • ,7)

п-Л-Vвs = ,

где ре — давление внешнего аэродинамического газового потока на поверхности тела, qe — тепловой поток на внешней поверхности Е1

qe = Ж в •п + е%а5Ввв - е^^ - ус8у. • ^. (8)

Здесь введены обозначения: ве — максимальная температура газового потока в пограничном слое, Vве — градиент температуры внешнего потока на поверхности конструкции, еи е. — коэффициенты излучения газа во внешнего потоке и твердого тела, а5Б — константа Сте-фана-Больцмана, у— коэффициента вдува внутреннего газа из пор во внешний поток.

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

п Vpg = 0,

(9)

пЛ^ Vв = 0.

Для расчета термонапряжений в теплозащитной конструкции используется уравнение равновесия упругих композитных сред с наличием термодеструктирующей матрицы [5]:

У-о-У((р) = 0, (10)

где о — тензор упругих напряжений.

К уравнениям (10) добавляются определяющие соотношения и соотношения Коши:

0

0 = 4 С-• (8 - 8),

0

8 = а(0-00)-Р, (11)

1 Т

8 = -(У® и + У® иТ ), 2

здесь обозначены: 8 — тензор малых деформаций, и — вектор пере-

0

" 4/-1

мещений, С — тензор модулей упругости, 8 — тензор тепловых деформаций, а — тензор теплового расширения, р — тензор химической усадки — все эти характеристики, а также тензоры Л и К зависят от температуры, концентраций фаз и матрицы преобразования Q/ от локального базиса анизотропии ri к базису единых осей координат [15]:

4С = (((( АМ'аГ&О!*. ®е;, ®е*. ®е,, (]2)

а = С((((,0,Щ е, ®ег, а = (Р,а,Л,К},

где С■ и с■ — компоненты тензоров в локальном базисе .

Конкретные выражения для функций (12) представлены в [5]. Композиционный материал конструкции полагается криволинейно-ортотропным.

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

п-о = -(Ре )п. (13)

На остальных частях поверхности задана либо свободная поверхность, либо жесткая заделка

п-о-п(( р) = 0 или и=0. (14)

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

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

о = Се. ®е= ссг'7'г., ®г,.,.

i ] ' ]

Численные методы решения сопряженной задачи. Для решения задачи аэродинамики (1)-(3) используется конечно-объёмный метод второго порядка аппроксимации с TVD реконструкцией [16, 17], подробно реализация этой численной методики решения описана в [18, 19].

Для решения задачи внутреннего тепломассопереноса (4)-(6) запишем вариационную постановку этой задачи. Домножим первое уравнение системы (4) на вариацию концентрации полимерной фазы 5<ь, второе — на вариацию порового давления 5р& и третье — на вариацию температуры 50,, и проинтегрируем полученные уравнения по объёму конечного элемента:

\рь5«ь <dV = -\5<ьШ,

{5р dV -\5р8 У • (К • Ур )dV = \5pjJYdV,

—^■V V

-0

{pscs50s—^dV

V -

(15)

{50У • (Л • )dV +1св50у0, • К • Ур^ -{50sJ^e0dV.

V

Далее, используя формулу Остроградского-Гаусса и граничные условия (7), преобразуем интегралы, содержащие производные тензоров газопроницаемости и теплопроводности:

\5ргV• (К• Ур)dV = {V • (5р • К• Ур)dV

V

-\У5рг • К• Ур^ИУ = {п5р • К • Ур^ -

V Е

-\У5р8 • К• Ур^ = \У5р8 • К • Ур^,

V V

{50У • (Л• )dV = {У • (50s • Л• )dV -

V V

-{У50х • Л • = { п50, • Л • -

V Е

~{У50Х • Л • = { 50^8 - {У50, • Л • .

(16)

Подставляя (16) в (15) и используя уравнение Менделеева-Клапейрона, получим вариационную постановку:

\рь5(рь др-сУ = -\&рьШ,

dV =

\5p-g- -\5Рз?М ^

V 8 Rв дг V 8 Rв дг

¡Чр8 • К • Ч8р^У + ¡5р8ШУ,

¡р^ = |• л• че^у +

(17)

¡ зерр8е8 чех • к •чр^у - ¡ е Ле^у

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

Рь = [Ф}т(и), Р8 = [Ф]т{у}, е,=[Ф]т

(18)

где и, у, 3 — значения соответствующих функций в узлах конечного элемента, а (Ф)т = (Ц,, L2, Ц3, L4) — строка функций форм. Аналогично будут выглядеть вариации неизвестных функций:

ръ = (Ф}т 3{и}, Р8 = (Ф}т 3{.у}, е = (Ф}т 3{3}.

Градиенты представим так:

ЧР8 = {Ц{Ф}т {у = [ В]{ у},

3 х 4

че, = {¿}{Ф}т {3} = [ В]{3},

(19)

(20)

где {Ц} =

(_д _д_ дх' ду' д2

— столбец дифференциальных операторов,

[ В] — матрица из производных функций форм.

3 х4

Тогда, подставляя соотношения (18)-(20) в вариационную постановку (17), вынося вариации неизвестных функций, получаем локальную систему уравнений для каждого конечного элемента (КЭ) в виде:

[М ] дЫ} + [ Я]{и] = 0,

дt

,д{ У} [6 ] д{Г} Ы

[Мв] ^Ыг+[кв]{Щ+[ШМ - [ Я»] = 0, дt

здесь введены обозначения для следующих матриц:

[Мр - [6 + [ Кр ]{у} - [ Яр ]{и} = 0, дt

тМе

Яв

dV,

[Мр ] = | {Ф}{Ф}

V

[Мв] = | {Ф}{Ф}Т рсЛУ,

V

[М ] = ръ | {Ф}{Ф}тdV,

[8 ] = | {Ф}{Ф}

Т РМ,

яв2

dV,

(21)

(22)

[Кр ] = | [ Б]т-[К ]•[ В] dV,

V 3x4 3x3 3 х4

[ К»] = | [ В]т •[Л] •[ В] dV,

V 3x4 3x3 3x4

[Я] = J01 ехр |-]{Ф}{Ф}^,

[ Яр ] = П01 ехр |»]{Ф} {Ф}т dV,

[0в ] = Ае0 J01 ехр Я» {Ф} {Ф}т dV, [ Я»] = | Че {Ф^8 + | {Ф}^dV,

^ = сяV » • К • ^ = ся {Г}т • [В]т • [К] • [В] • {у}.

3 x4 3 x3 3 x 4

Полученная система уравнений (21) является нелинейной системой обыкновенных дифференциальных уравнений, поскольку все матрицы М, Мр, Мв, Б, Кр, Кв, О», Я, Яр и вектор правых частей Яв

зависят от неизвестных функций сръ, р, в .

Для решения системы (21) воспользуемся методом линеаризации и неявной схемой Эйлера:

6'

([М ]п + М[К]" ){и}п+1 = [М ]"{и]п,

~[Яр ]п{и}п+1 + ([Мр ]п +М[Кр ]п К^Г - [5 ]п {¿Г = \ (23)

= [Мр]п{у}п -[5]п{3}п,

тп Мп+1 + ([М,]п + А[ Кв]п ){3}п+1 = [Мв]п {3}п + А[ ЯЛ

п — индекс шага интегрирования по времени. Таким образом, в (23) получаем локальную систему линейных уравнений. Собирая её для всех конечных элементов, получим систему линейных алгебраических уравнений (СЛАУ) с разреженной, несимметричной матрицей. Для решения СЛАУ используется численный метод бисопряжённого градиента.

Конечно-элементный метод решения задачи термомеханики.

При решении задачи термомеханики (10)-(14), как и для системы уравнений тепломассопереноса, будем использовать метод конечного элемента (МКЭ). Выпишем вариационную постановку уравнения уравнения (11), для чего умножим его на вариацию перемещения ¿и и проинтегрируем по объёму конечного элемента:

(о -ад Е) = 0. (24)

V

Далее рассмотрим цепочку равенств, используя формулу интегрирования по частям, формулу Остроградского-Гаусса и граничные условия (13), (14):

^и -V- (о -(р Е)^ =

V

= IV - [(о - адЕ) - 5и]<ЗУ -1'V ® ¿и - (о - адЕ)с1У =

"г ? (25)

= | п - (о - Е) - ¿иЕ -1'V ® ¿и - (о - Е)^" =

I V

= -1 ¿и - решЖ -1V ® ¿и - (о - (р8р8Е)^ = 0.

11 V

Затем в последнем интеграле формулы (22) перейдём от вариации перемещений к вариации тензора деформаций Коши 8г, что возможно в силу симметричности тензора напряжений о и единичного тензора Е

-1 ¿и - репй$ - - (о - адЕ)^ = 0. (26)

11 V

Далее подставляем определяющее соотношение (11) для тензора напряжений и получаем искомую вариационную постановку:

рпс® -^е--4 С к1У + • 4 С-е ёУ

Е! У У

ШУ = 0.

+

(27)

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

и = [Ф](у) =

3x12 12

0 0 Ь2 0 0 Ьъ 0 0 Ь4 0 0 0 ^00 Ь2 0 0 Ьъ 0 0 Ь4 0 0 0 Ц 0 0 Ь2 0 0 Ьъ 0 0 ь4

{V}, (28)

где {V} представляет собой столбец перемещений в узлах 4-х узлового

12

тетраэдального КЭ. Введём дифференциальный оператор [Б]:

[Б] =

д дх 0 0

0 д ду 0

0 0 д дг

д д 0

ду дх

д 0 д

дг дх

0 д д

дг ду

с помощью которого компоненты тензора деформаций можно представить в виде:

е = [Б]и = [ В][Ф]^} = [ B]{v}.

6х3 3x12 12 6x12 12

(29)

Независимые компоненты тензора модулей упругости можно представить как квадратную матрицу размера 6х6 [15]. Подставляя (28), (29) в вариационную постановку (27) и вынося за скобки вариации тензора деформаций, получаем следующую систему уравнений:

I [ В]Т [С ][ Б][у] dV = Г р [Ф]» dS -Г [ В]т [С ]{*} dV

,, 6х12 6x6 6x12 12 V 12х3 3 ,, 6x12 6x6 6

V V (30)

ЧЧ^[В]т {I}dV.

V 6x12 6

где {I} = (1,1,1,0,0,0)т — единичная матрица в векторном представлении.

Линеаризуя систему уравнений (30), получаем следующую локальную СЛАУ:

[ к ]{у} = {/}, (31)

12 x12 12 12

где обозначены матрицы

[к] = Г[В]т [С][В]dV,

12x12 V 6x12 6x6 6x12

0 (32)

{Г} = I Р [Ф] • {«} dS - I[В]т [С]{*} dV - Г(ргр [В]т {I} dV.

12 £ 12x3 3 V 6x12 6x6 6 V 6x12 6

Собирая все локальные СЛАУ (31), получаем глобальную СЛАУ с симметричной матрицей, которая решается методом сопряжённого градиента.

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

локальным векторам базиса в системе координат, связанной с поверхностью конструкции. Введем адаптивную систему координат

X', i = 1,3, которая согласована с базовыми поверхностями, ограничивающими область V рассматриваемой конструкции, и построим криволинейные блоки, охватывающие геометрию конструкции (рис. 1.), подобно тому, как строится расчетная область для генерации блочно-структурированных геометрически-адаптивных сеток [12]:

x = Г {X1), (33)

где x1 — исходные декартовы координаты, /' — известное преобразование перехода, у = 1,3.

С помощью преобразования координат (33) найдем вычислим локальные векторы базиса вдоль криволинейных направлений анизотропии в любой заданной точке (рис. 1.):

дх1

Г = , & =дХ-. (34)

дХ

Рис. 1. Визуализация поля локальных векторов базиса для учета свойств анизотропии материала

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

Результаты численного моделирования. Численное решение задачи было проведено для модельного летательного аппарата (ЛА), движущегося со скоростью М = 7 . Конструкция ЛА представляет собой конус, затупленный по сфере (рис. 2). В качестве композиционных материалов конструкции были выбраны: фенольный стеклопластик (сферическая область затупления) и эпоксидный стеклопластик (коническая оболочка). Константы материалов были взяты из [5].

На рис. 2 показано распределение температуры в\ в элементах

конструкции модельного ЛА. Наиболее теплонагруженной областью является критическое затупление, где температура достигала значений 1926 К. Под действием высоких температур после определенного времени началось интенсивное терморазложение фенольной и эпоксидной матриц (рис. 3), в результате чего образовывалось большое количество газов в порах обоих типов композиционных материалов. На рис. 4 показано распределение порового давления внутри элементов

конструкции в некоторый момент времени ^ . Максимальные значения

порового давления (свыше 35 МПа) достигались в конической оболочке из композита на эпоксидной матрице, несмотря на то, что для данной области уровень температур ве ниже, чем в области сферического затупления. Этот эффект объясняется наличием более высокой первичной пористостью в композите на фенольной матрице, по сравнению с эпоксидной и более высокой термостойкостью данного типа композита.

I 1,926**03 1-В315С+03 1,7Э67е*03 иыюе+ОЗ 1,54716 »03 1Л52Э14Ю 1.357Зе »03 й 1.26274МОЗ 1.1679е*03 || 1-07Э1*»03 № 9.783 №»02

7-8872е+02 6.9392**02 5.991Л+02 51МЗЭе*02 4Л93Э|«02 З.Н7Эе*02 7.1094**02

Рис. 2. Распределение температуры ^ (К) в модельной композитной конструкции

Рис. 3. Распределение концентрации полимерной матрицы фъ в композитных элементах модельной конструкции

Рис. 4. Распределение порового давления р (Па) в композитных элементах модельной конструкции

На рис. 5 показаны перемещения по оси г . Видно, что в области затупления перемещения отрицательные, что связано с тепловыми расширениями композита под действием высоких температур. На рис. 6 представлены распределения поперечных напряжений <х33 в конструкции модельного ЛА в некоторый момент времени ^ . Их максимальное значение достигается в области сферического затупления на некоторой глубине от нагреваемой поверхности.

Рис. 5. Распределение перемещения по оси г (м) в модельной конструкции

12.179Л-С1 UMJtfll ™ ] hSHrOI

lAXSMtl l.lSXtOl 5ЛТМ0)>

шшв

1

■I.KiBrOi

ri ueoif o;

I ci

I -3JN69C

P -BJJ747(-CZ 1ЛШ>Ч1|

У i ar,.v : i I UHlrai

I-ijoxrci -iJBMHU ■ZJOOSrOl

Рис. 6. Распределение поперечного напряжения сг3 3 (ГПа) в модельной композитной конструкции

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

ЛИТЕРАТУРА

[1] Riccio A., Raimondo F., Sellitto A., Carandente V., Scigliano R., Tescione D. Optimum design of ablative thermal protection systems for atmospheric entry vehicle. Applied Thermal Engineering, 2017, no. 119, pp. 541-552.

[2] Lancelle D., Bozic O. Simulation of an ablative thermal protection system for the hypersonic ascend of an electromagnetically launched payload carrier. Proceedings of 5th European Conference for Aeronautics and Space Sciences (EUCASS), 2013, 12 p.

[3] Eekelen T., Bouilly J.-M., Hudrisier S., Dupillier J.-M., Aspa Y. Design and numerical modelling of charring material ablators for re-entry applications. Proceedings of the Sixth European Workshop on Thermal Protection Systems and Hot Structures, University Stuttgart, Germany, 2009, European Space Agency -WPP319.

[4] Liu Z., Hao A., Zhang S., Dessureault Y.-S., Liang R. Lightweight carbon nanotube surface thermal shielding for carbon fiber/bismaleimide composites. Carbon, 2019, vol. 153, pp.320- 329. DOI: 10.1016/j.carbon.2019.07.018

[5] Димитриенко Ю.И. Механика композитных конструкций при высоких температурах. Москва, Физматлит, 2018, 448 с.

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

[7] Dimitrienko Yu.I. Thermomechanical behaviour of composite materials and structures under high temperatures. Part 2. Structures. Composites. Part A: Applied Science and Manufacturing, 1997, vol. 28A, рр. 463-471.

[8] Димитриенко Ю.И., Захаров А.А., Коряков М.Н., Сыздыков Е.К., Минин В.В. Численное решение сопряженной задачи гиперзвуковой аэродинамики и термомеханики термодеструктирующих конструкций. Инженерный журнал: наука и инновации. 2013, вып. 9 (21).

URL: http://engjournal.ru/catalog/mathmodel/aero/1114.html

[9] Димитриенко Ю.И., Захаров А.А., Коряков М.Н., Сыздыков Е.К. Моделирование сопряженных процессов аэрогазодинамики и теплообмена на поверхности теплозащиты перспективных гиперзвуковых летательных аппаратов. Известия высших учебных заведений. Машиностроение, 2014, № 3, с. 2334.

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

[11] Димитриенко Ю.И., Коряков М.Н., Захаров А.А. Разработка технологий численного моделирования сопряженных задач газовой динамики и термодинамики композитных конструкций перспективных высокоскоростных летательных аппаратов. Известия Самарского научного центра Российской академии наук, 2016, т. 18, № 2 (3), с. 891-895.

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

[13] Краснов Н.Ф. Аэродинамика. Т.1. Основы теории. Аэродинамика профиля и крыла. Москва, Высшая школа, 1980, 496 с.

[14] Краснов Н.Ф. Аэродинамика. Т.2. Методы аэродинамического расчета. Москва, Высшая школа, 1980, 416 с.

[15] Димитриенко Ю.И. Механика сплошной среды. В 4 т. Т. 1. Тензорный анализ. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2011, 367 с.

[16] Anderson J.D. Hypersonic and High-Temperature Gas Dynamics. 2nd edition. American Institute of Aeronautics and Astronautics, Reston, Virginia, 2006, 232 p.

[17] Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. Москва, Физматлит, 2012, 656 с.

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

[19] Димитриенко Ю.И., Коряков М.Н., Захаров А.А. Применение метода RKDG для численного решения трехмерных уравнений газовой динамики на неструктурированных сетках. Математическое моделирование и численные методы, 2015, № 4, с. 75-91.

Статья поступила в редакцию 30.09.2019

Димитриенко Ю.И., Коряков М.Н., Юрин Ю.В., Захаров А.А. Конечно-элементное моделирование термонапряжений в композитных термодеструктирующих конструкциях при аэродинамическом нагреве. Математическое моделирование и численные методы, 2019, № 2, с. 15-34.

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

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

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

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

Finite-element modeling of thermal stresses in composite structures with thermal decomposition under aerodynamic

heating

© Yu.I. Dimitrienko, M.N. Koryakov, Yu.V. Yurin, A.A. Zakharov Bauman Moscow State Technical University, Moscow, 105005, Russia

The coupled task of aero-thermo-mechanics of heat-loaded structures from thermally destructive polymer composite materials under the influence of an intense aerodynamic flow is considered. The mathematical formulation of the conjugate problem is formulated and algorithms for the numerical solution of this problem are proposed. The algorithms are based on an iterative solution of three types ofproblems: aerodynamics, internal heat and mass transfer, and thermomechanics of the modeling aircraft structure. An example of a numerical solution to the problem for an aircraft structural element in the form of a blunt cone is presented. It is shown that the effect of high temperatures of aerodynamic heating of the structure leads to thermal degradation of the polymer composite material, which results in the generation of a large amount of gases in the pores and thermo-chemical shrinkage, which under certain conditions can lead to premature destruction of the heat-loaded composite structure.

Keywords: conjugate problems, aero-thermomechanics, heat and mass transfer, thermal destruction, thermomechanics, composite materials, thermal deformation, pore pressure, thermal stresses, finite element method

REFERENCES

[1] Riccio A., Raimondo F., Sellitto A., Carandente V., Scigliano R., Tescione D. Optimum design of ablative thermal protection systems for atmospheric entry vehicle. Applied Thermal Engineering, 2017, no. 119, pp. 541-552.

[2] Lancelle D., Bozic O. Simulation of an ablative thermal protection system for the hypersonic ascend of an electromagnetically launched payload carrier. Proceedings of 5th European Conference for Aeronautics and Space Sciences (EUCASS), 2013, 12 p.

[3] Eekelen T., Bouilly J.-M., Hudrisier S., Dupillier J.-M., Aspa Y. Design and numerical modelling of charring material ablators for re-entry applications. Proceedings of the Sixth European Workshop on Thermal Protection Systems and Hot Structures, University Stuttgart, Germany, 2009, European Space Agency -WPP319.

[4] Liu Z., Hao A., Zhang S., Dessureault Y.-S., Liang R. Lightweight carbon nanotube surface thermal shielding for carbon fiber/bismaleimide composites. Carbon, 2019, vol. 153, pp.320- 329. DOI: 10.1016/j.carbon.2019.07.018

[5] Dimitrienko Yu.I. Mekhanika kompozitnyh konstrukcij pri vysokih temperaturah [Mechanics of composite structures at high temperatures]. Moscow, Fizmatlit Publ., 2009, 624 p.

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

[7] Dimitrienko Yu.I. Thermomechanical behaviour of composite materials and structures under high temperatures. Part 2. Structures. Composites. Part A: Applied Science and Manufacturing, 1997, vol. 28A, рр. 463-471.

[8] Dimitrienko Yu.I., Zakharov A.A., Koryakov M.N., Syzdykov E.K., Minin V.V. Inzhenernyy zhurnal: nauka i innovatsii — Engineering Journal: Science and Innovation, 2013, no. 9. Available at: http://engjournal.ru/catalog/mathmodel/aero/! 114.html

[9] Dimitrienko Yu.I., Zakharov A.A., Koryakov M.N., Syzdykov E.K. Izvestiya vys-shikh uchebnykh zavedeniy. Mashinostroenie — Proceedings of Higher Educational Institutions. Machine Building, 2014, no. 3, pp. 23-34.

[10] Dimitrienko Yu.I., Koryakov M.N., Zakharov A.A., Stroganov A.S. Matematich-eskoe modelirovanie i chislennye metody — Mathematical modeling and Computational Methods, 2014, no. 3 (3), pp. 3-24.

[11] Dimitrienko Yu.I., Koryakov M.N., Zakharov A.A. Izvestiya Samarskogo nauch-nogo centra Rossijskoj akademii — nauk Izvestia of Samara Scientific Center of the Russian Academy of Sciences, 2016, vol. 18, no. 2 (3), pp. 891-895.

[12] Dimitrienko Yu.I., Kotenev V.P., Zakharov A.A. Metodlentochnykh adaptivnykh setok dlya chislennogo modelirovaniya v gazovoy dinamike [Adaptive tape grid method for computational simulation in gas dynamics]. Moscow, FIZMATLIT PUBL., 2011, 280 p.

[13] Krasnov N.F. Aerodinamika. T.1. Osnovy teorii. Aerodinamika profilya i kryla [Aerodynamics. Vol. 1. Fundamentals of theory. Airfoil and wing aerodynamics]. Moscow, Vysshaya shkola Publ., 1980, 496 p.

[14] Krasnov N.F. Aerodinamika. T.2. Metody aerodinamicheskogo rascheta [Aerodynamics. Vol. Methods of aerodynamic calculation]. Moscow, Vysshaya shkola Publ., 1980, 416 p.

W.M. ffuMumpueuKO, M.H. KopnKoe, W.B. WpuH, A.A. 3axapoe

[15] Dimitrienko Yu.I. Tensor analysis and nonlinear tensor functions. Springer, 2002, 662 p.

[16] Anderson J.D. Hypersonic and High-Temperature Gas Dynamics. 2nd edition. American Institute of Aeronautics and Astronautics, Reston, Virginia, 2006, 232 p.

[17] Kulikovskiy A.G., Pogorelov N.V., Semenov A.Yu. Matematicheskie voprosy chislennogo resheniya giperbolicheskikh sistem uravneniy [Mathematical problems on numerical solution of hyperbolic equation systems]. Moscow, Fizmatlit Publ., 2012, 656 p.

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

[19] Dimitrienko Yu.I., Koryakov M.N., Zakharov A.A. Matematicheskoe modeliro-vanie i chislennye metody — Mathematical modeling and Computational Methods, 2015, no. 4, pp. 75-91.

Dimitrienko Yu. I., Professor, Head of the Department of Computational Mathematics and Mathematical Physics, Director of Scientific-Educational Center of Supercomputer Engineering Modeling and Program Software Development, Bauman Moscow State Technical University. Member of the Russian Academy of Engineering Science. Author of over 400 research publications in the field of computational mechanics, gasdynamics, thermo-mechanics of composite materials, mathematical simulations in material science. e-mail: dimit.bmtstu@gmail.com

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

Yurin Y.V., Cand. Sc. (Phys.-Math.), Assoc. Professor, Department of Computational Mathematics and Mathematical Physics, Bauman Moscow State Technical University; Senior Research Fellow, Scientific and Educational Centre for Supercomputer-based Engineering Simulation and Software Development, Bauman Moscow State Technical University. e-mail: yvyurin@yandex.ru

Zakharov A.A., Cand. Sci. (Phys. & Math.), Associated 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