CC BY
4Труды МАИ. Выпуск № 101
http://trudymai.ru/
УДК: 533, 536.2, 519.63
Численное моделирование взаимодействия аргоновой плазмы с углеродным образцом теплозащитного покрытия
Савицкий Д.В.*, Аксёнов А.А.**, Жлуктов С.В.***
Объединённый институт высоких температур РАН, Ижорская ул. 13, Москва, 125412, Россия *e-mail: dmvlsav@yandex. ru **e-mail: andrey@tesis.com.ru ***e-mail: sz@flowvision.ru
Аннотация
В работе представлена комплексная математическая модель, описывающая течение горячей смеси газов около летательного аппарата (ЛА), унос массы с поверхности ЛА, влияние продуктов абляции на процессы, протекающие в газовой фазе около ЛА, прогрев теплозащитного покрытия (ТЗП) и изменение формы поверхности ЛА. Моделируется обтекание образца углеродного ТЗП аргоновой плазмой. Условия соответствуют эксперименту, проведённому в ОИВТ РАН. Предполагается, что основным механизмом абляции (уноса массы) является сублимация. Расчёты проводятся в программном комплексе FlowVision. Приводятся результаты численного моделирования. Проводится сравнение численных результатов с экспериментальными данными.
Ключевые слова: гиперзвуковой летательный аппарат, теплозащитное покрытие, плазма, абляция, сопряжённый теплообмен, уравнения Навье-Стокса, численное интегрирование, метод расщепления по физическим процессам.
Введение
В общем случае реальная картина физико-химических процессов около ГЛА очень сложная. При торможении газа в головной ударной волне почти вся кинетическая энергия набегающего потока переходит в тепловую, приводящую к нагреву газа до высокой температуры (десятки тысяч градусов) непосредственно за ударной волной, которая далее, при приближении к обтекаемому телу, спадает за счёт эндотермических реакций диссоциации и ионизации [1]. При столь высоких температурах времена релаксации колебательных и электронных степеней свободы молекул сравниваются с временами реакций и становится необходимым учитывать колебательно-диссоциационное и электронно-ионизационное взаимодействия [2]. Таким образом, в ударном слое реализуется течение многокомпонентного термохимически неравновесного газа с разными температурами поступательных, колебательных и электронных степеней свободы. При большой скорости входа в атмосферу Земли (> 8 км/^ и большом радиусе затупления ГЛА (> 5 м) необходимо учитывать радиационный поток энергии, идущий от горячей газовой смеси к ТЗП [1].
Взаимодействие нагретого газа с ТЗП инициирует многочисленные дополнительные взаимосвязанные процессы [3, 4]: гетерогенные химические
реакции на поверхности ГЛА, гомогенные химические реакции с участием продуктов разрушения ТЗП, сопряжённый теплообмен между газом и ТЗП, нестационарный прогрев многослойного ТЗП. Необходимость описывать эти процессы приводит к существенному усложнению математической модели. Моделирование взаимодействия горячего газа с ТЗП определяется типом ТЗП. В книгах [3, 4] рассматриваются различные ТЗП и предлагаются простые модели для расчётов скорости уноса массы с поверхности ТЗП (модели абляции поверхности ГЛА). Одним из способов защиты ГЛА от тепловых нагрузок является покрытие его поверхности углеродом. В работе [5] приводятся условия полёта экспериментального ГЛА IRV-2. Передняя часть аппарата была покрыта графитом. Скорость входа в атмосферу Земли была меньше первой космической. При спуске ГЛА в атмосфере фиксировалась скорость уноса массы с поверхности ГЛА и форма поверхности ГЛА в различных точках траектории. В той же работе предлагается модель абляции углеродного ТЗП при его взаимодействии с частично диссоциированным воздухом. Углеродные ТЗП отличаются кристаллической структурой. В работе [6] предложена математическая модель процессов в вязком ударном слое и на поверхности спускаемого аппарата, покрытого пиролитическим графитом. Модель включает 19 газовых компонентов, 40 реакций в газовой фазе и 12 гетерогенных реакций на поверхности графита (с участием адсорбированных атомов O и N). Работа выполнялась в рамках проекта "MUSES C Asteroid Sample Return Mission" [7]. Возвращаемая капсула вошла в атмосферу Земли со второй космической скоростью. Пиролитический графит имеет наиболее
структурированную (многослойную) кристаллическую решётку и поэтому способен выдерживать большие тепловые нагрузки при небольшой скорости абляции.
В настоящей работе рассматривается обтекание образца углеродного теплозащитного покрытия аргоновой плазмой. Условия соответствуют эксперименту, проведённому в ОИВТ РАН [8]. Образец изготовлен из графита марки МПГ-6. Данный тип искусственного графита образуется прессованием непрокаленного нефтяного кокса и последующей графитацией при температуре около 3000 К. В результате получается мелкозернистый графит, с размером кристаллитов 30-150 мкм. Это обеспечивает изотропность физико-химических свойств образца. Исследуется механизм сублимации (испарения) такого образца[9]. Именно поэтому в эксперименте [8] в качестве газа использовался аргон. Расчёты проводятся в программном комплексе (ПК) FlowVision [10].
1. Математическая модель
Задача решается в сопряжённой постановке. В газовой фазе интегрируются полные уравнения Навье-Стокса. В образце ТЗП интегрируется уравнение энергии. На контактных поверхностях устанавливаются сопряжённые граничные условия, моделирующие балансы массы и энергии. Предполагается, что основным механизмом абляции образца является сублимация. Такая постановка позволяет упростить модель газовой фазы. Можно предположить, что газовая фаза состоит из двух компонентов: аргона и продуктов абляции. В настоящей работе продукты
абляции ТЗП рассматриваются, как одно вещество. Модель позволяет рассчитывать нестационарный прогрев образца.
1.1. Модель газовой фазы Математическая модель газовой фазы основана на полных уравнениях Навье-Стокса для сжимаемого газа. Она включает уравнения (конвективно-диффузионного типа) в частных производных и алгебраические соотношения. Уравнение неразрывности:
|+уО,к)=о (1)
дг
Здесь р - плотность, V - скорость, г - время. Уравнение импульсов:
+ ® У) = -Ур + V• Т (2)
Т = /2 Л - 2 (V V)11 (3)
(дК 8У, 1 (4)
2
—^ + —]
Здесь р - статическое давление, т - тензор вязких напряжений, § - тензор скоростей деформации, р - динамический коэффициент вязкости, 1 - метрический тензор. Уравнения массопереноса (в общем случае):
^•(л у)^ I — а <5)
дг
Здесь ^ - массовая доля компонента i. В настоящей работе предполагается, что уравнения для аргона и для продуктов абляции - однородные. Поэтому а — 0. В
предположении закона Фика выражение для диффузионного потока компонента I имеет вид:
= - рПУТ1 уг1 (6)
Ье
Здесь Б - коэффициент диффузии, Ье = /и/ pD - число Шмидта. Уравнение энергии решается относительно полной энтальпии.
(7)
+У(р ун )=|-V-1, + У-
01 СЯ
Л 2£- 2(У-У)1 1-У
Здесь Н = к + V2 / 2 - полная энтальпия, И - термодинамическая энтальпия. Выражение для теплового потока имеет вид:
N
19 =- хуТ + £ к 1г 1=1
Или
1 = - А (УН - V - (V - У)) + £ ^ (1 - Ье)к1УУ1 (9)
С р ¿=1 С р
Здесь 1 - коэффициент теплопроводности смеси, Ср - теплоёмкость смеси, Ьв = Рг/Ье - число Льюиса, Рг = цС /Я - число Прандтля. Система (1) - (9)
замыкается уравнением состояния идеального газа
рт_ (10)
*лТ
и общими соотношениями
(11)
(8)
N
Тъ = 1
1=1
1 N У - =
т *~=- т,
Здесь т - молярная масса компонента i, т - молярная масса смеси, N - число компонентов газовой фазы.
Излучение плазмы не учитывается. Ионизация аргона учитывается в его свойствах - предполагается, что свойства аргона определены условием локального термо-химического равновесия. Молярная масса, теплоёмкость, вязкость и теплопроводность задаются в интерфейсе FlowVision в виде таблиц (р, Т). Плазмотрон создаёт ламинарный поток плазмы [8]. Поэтому реализованные в ПК FlowVision модели турбулентности [10] в расчётах не используются. Конвективно-диффузионные уравнения (1-9) записаны в консервативной форме. Такой формы записи требует метод конечных объёмов - см. раздел 3.
1.2. Модель твёрдого тела
В твёрдом теле (в образце) решается уравнение энергии относительно термодинамической энтальпии:
(13)
дг 9
Jq =-^н. (14)
ср
1.3. Модель контактной поверхности
Математическая модель процессов на поверхности ТЗП включает алгебраические соотношения, описывающие баланс массы и баланс энергии. В этих балансах учитывается унос массы с поверхности образца. Предполагается, что
основным механизмом абляции является сублимация. Известно, что при сублимации углерода образуются молекулярные комплексы С, С2, С3, С4, С5, С6, С8. Спектрометрия [11] показывает, что наиболее массовым компонентом газовой смеси, образующейся на поверхности углеродного ТЗП, является С2. В настоящей работе предполагается, что единственным продуктом абляции является газ С2. Балансы массы и энергии определяют сопряжённое граничное условие для уравнений энергии в газовой и твёрдой фазах, а также граничные условия для уравнений массопереноса в газовой фазе.
1.3.1. Баланс массы Полный поток /-го компонента с поверхности тела определяется соотношением:
т = Р, У „ (15)
Здесь р,,„ и ум - плотность и нормальная составляющая скорости /-го компонента
на поверхности ТЗП (на стенке). Суммарный полный поток массы (скорость уноса массы с поверхности тела):
(16)
N N
т = 2т, =^Рг;Уг,„ = РмУм
г=1 г=1
Здесь р„ и Ум - плотность и нормальная составляющая скорости газовой смеси на поверхности ТЗП. Диффузионный поток /-го компонента с поверхности ТЗП по определению выражается следующим образом:
( т Л лг ^ (17)
— - — п V — ГР1 --— = /4? — --—
г, м
у = р (у _ у )=р. У _ — =р у _ т Рг,м = т _ т
° гИг,г,м ум) Иг,м у г,м Иг,М г,м "1г
V р„ У рм р„
Рм Рм
Поскольку в газе
= Х М± = у (18) р пМ 1 М 1'
получаем
3 . = тг - тУ1, . . (19)
-5
В выражении (18) п - концентрация частиц /-го компонента [м- ], Мг - масса /-го компонента (молекулы, атома) [кг], М - средняя масса смеси [кг], X - мольная доля /-го компонента (= объёмная в газе). Если предположить, что между стенкой и центром примыкающей к стенке ячейки выполняется закон Фика, то выражение для диффузионного потока /-го компонента приобретает вид
=-рБУУг = щ -ту (20)
Ув
Здесь ув - расстояние от центра пристенной ячейки до стенки по нормали. В
результате получили необходимое балансовое соотношение. Способ вычисления удельной скорости уноса массы /-го компонента определяется режимом (моделью) абляции. Суммарный поток массы т определяет локальную скорость изменения формы разрушающейся поверхности.
1.3.2. Баланс энергии Баланс энергии на разрушающейся поверхности можно записать так:
3д, В + 1 • ИВ Т ) + аЪя = 3д,з + 1 - К (Т. ) (21)
Здесь 3я - нормальная составляющая конвективного теплового потока, к(тм,) -
термодинамическая энтальпия, е - степень черноты стенки, индексами g и 8, соответственно, обозначены величины, относящиеся к газовой и твёрдой фазам.
Радиационный поток энергии, падающий на стенку из газовой фазы, в данной работе не учитывается. По определению
Т - Т N (22)
У в 1=1
Ивв Т)=ТКв ТХ (23)
г=1
у _ ^е, з - Т* (24)
д 3 Уз
N
1 • Км Т )=£ к1 , ,(Т* )т 1
(25)
к1 , з Т )т 1
1=1
Преобразуем балансовое соотношение (21) с использованием выражений (22) - (25):
Т - Т N N (26)
Яв +£КТ 1т -т ^+т ^)+оеТ:ь = +2(Т*т
Ув 1=1 1=1
В результате получаем другую форму записи баланса энергии на разрушающейся поверхности:
Т - Т Т - Т N (И\
+ сеПдЫ (Т* )т (27)
в Ув Уз 1=1
Здесь е - степень черноты поверхности образца (в представленных ниже расчётах
использовалось значение е = 1), АИг (Т* ) = к1в (Т*)- к1з (Т*) - теплота фазового перехода
вещества / при температуре поверхности ТЗП.
1.3.3. Модель уноса массы
В настоящей работе рассматривается один механизм абляции ТЗП -
сублимация. Предполагается, что в поток газа со стенки вдувается один компонент.
2 1
Локальная интенсивность вдува массы т [кг м- с- ] определяется давлением
насыщенного пара продукта абляции при локальной температуре поверхности ТЗП Т^. Интенсивность вдува и температура поверхности - искомые величины. Поскольку в поток газа вдувается один компонент, баланс массы на поверхности образца упрощается:
и ^АВ^ _ т ( _ 1 )
о _ ГПV1 IABL,w)
& у „
(28)
Здесь УАВЬ, - массовая доля продукта абляции на стенке, 1АВЬс - массовая доля
продукта абляции в центре примыкающей к стенке ячейки. Соотношение (28) используется для вычисления интенсивности вдува т. При этом используется ограничение [3]
т <
Р ABL Т)
(29)
\2ж
т
'ABL
По определению
т
У у mABL
IABLХ ABL
т
(30)
т _ Е тХ^
1_1
(31)
pABL (Т..) р
у _ -Г ABL,SaТ \ .
Х АВВЬ,,.
(32)
Зависимость парциального давления паров продукта абляции от температуры
рАВ1,!
Т (Т.) задаётся в интерфейсе FlowVision в табличном виде - это одно из
задаваемых свойств газообразного вещества.
Зная т, можно вычислить температуру поверхности, используя балансовое
соотношение (27), которое в случае одного продукта абляции также упрощается:
N
. Tw Tc,g rp 4 1 Tc, s Tw Л1 (33)
*g -- + abs = К —-- AhABL (Tw )m
yg ys
При этом значения температуры в центрах примыкающих к поверхности ячеек Tcg, TCs s полагаются известными с предыдущего шага по времени.
Как видно из соотношений (32), (29), (33), m зависит от T, а T, в свою очередь, зависит от mm. Поэтому для нахождения этих величин организуется итерационный процесс.
2. Моделирование изменения формы образца
Для моделирования изменения формы образца используется технология VOF. Метод VOF (Volume Of Fluid) традиционно используется для моделирования движения свободных границ и контактных поверхностей между двумя несмешивающимися жидкостями. В данной работе предлагается распространить этот метод на контактную границу между газом и твёрдым телом.
2.1. Стандартный метод VOF Метод VOF, реализованный в ПК FlowVision, является усовершенствованной версией метода [12]. Перенос границы сплошной фазы описывается уравнением для объёмной доли этой фазы в расчётной ячейке (переменной VOF). Переменная VOF принимает значения от 0 (газ) до 1 (жидкость). Ячейка, в которой 0 < VOF < 1, содержит межфазную границу. В программе межфазная граница представляется набором многоугольников. Перенос контактной поверхности описывается конвективным уравнением для переменной VOF (F):
^ + у.^ _ 0 (34)
ат
После решения уравнения (34) строится контактная поверхность. Затем решаются уравнения неразрывности, импульса, турбулентности, массопереноса, энергии. На контактной поверхности автоматически устанавливаются условия сопряжения всех искомых переменных. Таким образом, реализованный в ПК FlowVision алгоритм не требует "взвешивания" свойств контактирующих сплошных фаз (плотности, вязкости, теплопроводности и т. д.) в ячейках, через которые проходит контактная поверхность (0 < т < 1).
2.2. Модифицированный метод УОР
В настоящей работе предлагается версия метода VOF, позволяющая моделировать движение поверхности ТЗП вследствие уноса массы. Эта версия предполагает решение неоднородного уравнения для твёрдого тела, поверхность которого разрушается:
& _ 0 (35)
аТ
Здесь Оавь - источниковый член, описывающий изменение объёма твёрдой фазы за счёт абляции в ячейках, где 0 < т < 1. Выражение для него имеет следующий вид:
0 _ тS (36)
Оав,ь р т п
Здесь £ - площадь контактной поверхности в данной ячейке, п - объём ячейки.
3. Численные технологии FlowVision
Численные технологии, реализованные в ПК FlowVision [10], основаны на конечно-объемном подходе к аппроксимации уравнений в частных производных. Метод конечных объёмов предполагает интегрирование уравнений, описывающих течение жидкости и протекающих в ней процессов, по объёмам ячеек расчётной сетки. По теореме Гаусса для произвольной векторной или тензорной величины ¥
Здесь Псе11 - объём ячейки, ASt - площадь /-ой грани, щ - внешняя нормаль к /-ой грани. Как видим, при интегрировании решаемых уравнений в ячейке производится суммирование потоков массы, импульса, энергии и других характеристик течения, вычисленных на гранях ячеек. Поскольку каждая грань разделяет две соседние ячейки, соответствующий поток входит с противоположными знаками в дискретные уравнения для обеих ячеек. Этим обеспечивается консервативность массы, импульса, энергии и других искомых величин в расчётной области.
Уравнения, представленные в разделе 1.1, интегрируются последовательно. Последовательность операций, производимых солвером FlowVision на шаге по времени, определяется методом расщепления по физическим процессам [13]. Расчётная область импортируется из произвольной системы автоматического проектирования в форматах STL, WRL и достраивается с использованием инструментов препроцессора FlowVision. Расчётная сетка строится автоматически. Начальная сетка - всегда трёхмерная декартова. Она охватывает всю расчётную область. Пользователь задаёт критерии адаптации сетки (по поверхности, в объёме,
(37)
по решению). После запуска задачи на расчёт генератор сетки FlowVision автоматически разбивает и укрупняет ячейки в соответствующих зонах расчётной области. При задании адаптации по решению сетка перестраивается на каждом шаге по времени, "отслеживая" движение ударной волны, вихря, или какого-то другого элемента течения. Триангулированная граница расчётной области может иметь произвольную криволинейную форму. Прямоугольные ячейки (параллелепипеды) рассекаются границей естественным образом. В результате, ячейки, примыкающие к границе, представляют собой многогранники произвольной формы. Несмотря на геометрическую сложность расчетной сетки FlowVision, порядок аппроксимации решаемых уравнений по пространству - второй. Это достигается использованием высокоточной расчетной схемы с реконструкцией решения внутри ячейки [14]. Уравнения в частных производных решаются с использованием неявного алгоритма [13]. Алгебраические уравнения, возникающие при аппроксимации уравнений, решаются методами крыловского типа. ПК FlowVision работает на компьютерах, имеющих смешанную архитектуру: межузловое МР1-распараллеливание на них производится одновременно с распараллеливанием на узле как на компьютере с общей памятью. Использование смешанного распараллеливания позволяет добиться высококачественного масштабирования программного комплекса при работе на большом числе процессоров [15], [16].
Обзор технологий FlowVision представлен в работе [10].
4. Тестовая задача
Для проверки модели сублимации и алгоритма расчёта изменяющейся формы образца проведено численное моделирования взаимодействия струи аргоновой плазмы с графитовым образцом [8]. В работе [8] изменение формы образца отслеживается с помощью метода лазерной профилометрии [17]. В эксперименте [8] плазмотрон [18] создает поток аргоновой плазмы с температурой 11000 К и скоростью 600 м/с на расстоянии 20 мм от среза сопла. Температура измеряется спектральными методами, скорость - с помощью охлаждаемой трубки Пито. Углеродный образец устанавливается на расстоянии 20 мм от среза сопла плазмотрона. Образец имеет форму параллелепипеда с размерами 23мм, 23мм, 16мм.
Рис. 1 - Фото образца во время эксперимента (слева) и по окончании (справа). При воздействии аргоновой плазмы на образец происходит сублимация углерода. В результате его форма меняется. Используя данные лазерного профилометрирования,
авторы работы [8] получают зависимость значения удельной скорости уноса массы от времени.
Зависимость парциального давления паров графита взята из ресурса [19] - см. Таблицу 2.
Таблица 2. - Зависимость давления насыщенного пара от температуры.
Абсолютная температура, K Давление, Па
0 0
2839 10
3048 100
3289 1000
3572 10000
3908 100000
Свойства газов Аг и С2 были рассчитаны в программе ИТВАН ТЕРМО в предположении индивидуального термохимического равновесия.
5. Расчёты
Для численного моделирования создана точная геометрическая модель экспериментальной установки - см. рис.2.
Рис. 2 - Расчетная область Перед проведением систематических расчетов было проведено исследование сеточной сходимости результатов. Задавалась адаптация сетки по поверхности образца и в объёме струи. Расчеты показали, что оптимальной является расчётная сетка, состоящая примерно из 600 тыс. ячеек. Сечение этой сетки плоскостью, проходящей через ось симметрии струи, показано на рис.3.
Рис. 3 - Расчетная сетка
В результате численного моделирования получено распределение температуры и скорости в потоке плазмы - см. рис. 4 и 5.
X
Рис. 4 - Распределение температуры через 60с, К
Рис. 5 - Распределение скорости через 60с, м/с
Из приведенных рисунков видно, что температура и скорость плазмы около поверхности образца через 60 секунд после начала эксперимента составляют 11000 К и 600 м/с. Это соответствует экспериментальным данным [8].
Предполагается, что основным продуктом абляции является газ С2 [11]. На рис. 6 показано распределение удельной скорости уноса массы в начале эксперимента и через 60 секунд, когда уже образовался заметный кратер на верхней поверхности образца.
2 1
Рис. 6 - Распределение скорости уноса массы [кг м- с- ] в начальный момент времени (слева) и через 60 с. (справа). На рис. 7 показано распределение массовой доли С2 через 60 секунд после начала эксперимента.
Рис. 7 - Массовая доля продуктов сублимации (С2) через 60 с.
Интегрирование по площади кратера в разные моменты времени дает среднее
-2 -1
значение скорости уноса массы от 5 до 6 мг см- с- . Несколько значений показано на рис.8.
Время, с
-2 -1
Рис. 8 - Средняя удельная скорость уноса массы [кг м- с- ]: результаты расчётов - маркеры, экспериментальные данные - линия.
Авторы работы [8] утверждают, что пик на экспериментальном графике связан с тем, что интенсивный неравномерный прогрев образца в первые 40 секунд приводит к возникновению в нём значительных внутренних напряжения, что, в свою очередь, инициирует унос массы посредством выкрашивания. Через 40 секунд после начала эксперимента основным механизмом разрушения образца становится сублимация. Таким образом, результаты численного моделирования с использованием реализованной модели абляции согласуются с экспериментальными данными в пределах погрешности эксперимента.
Выводы
Необходимость в новых теплозащитных материалах мотивирует разработку и реализацию комплексных математических моделей, описывающих течение горячей смеси газов около ЛА, унос массы с поверхности ЛА, влияние продуктов абляции на процессы, протекающие в вязком ударном слое около ЛА, прогрев ТЗП и изменение формы поверхности ЛА. При этом такие модели не должны быть слишком сложными - в противном случае ресурсов даже очень мощного кластера не хватит для проведения систематических инженерных расчётов характеристик реальных летательных аппаратов. По этой причине модели отдельных процессов должны использовать эмпирические (экспериментальные) данные и, при необходимости, настраиваться на условия решаемой задачи. В данной работе приводится математическая модель, описывающая основные процессы, протекающие в газовой
фазе около ЛА, прогрев ТЗП, процессы на поверхности ТЗП, а также изменение геометрии ЛА вследствие уноса массы. Условия соответствуют эксперименту [8]. В этом эксперименте образец углеродного ТЗП помещается в струю аргоновой плазмы. Предполагается, что основным механизмом абляции поверхности образца является сублимация. Такая постановка позволяет упростить модель газовой фазы, а именно, предположить, что газовая фаза состоит из двух компонентов: аргона и продуктов абляции. Расчёты проводятся в ПК FlowVision. Полученные в расчётах значения скорости уноса массы с поверхности образца отличаются от экспериментальных не более чем на 20%. Это примерно соответствует погрешности эксперимента [8].
Работа была выполнена с использованием оборудования центра коллективного пользования «Комплекс моделирования и обработки данных исследовательских установок мега-класса» НИЦ «Курчатовский институт» (субсидия Минобрнауки, идентификатор работ RFMEFI62117X0016), http://ckp.nrcki.ru/.
Библиографический список
1. Тирский Г.А., Сахаров В.И., Ковалев В.Л., Власов В.И., Горшков А.Б. Гиперзвуковая аэродинамика и тепломассообмен спускаемых космических аппаратов и планетных зондов. - М.: Физматлит, 2011. - 545 с.
2. Никитин П.В., ШкуратенкоА.А. Влияние каталитически активной поверхности на интенсивность конвективного теплообмена // Труды МАИ. 2016. № 88. URL: http://trudymai.ru/published.php?ID=70574
3. Панкратов Б.М., Полежаев Ю.В., Рудько А.К. Взаимодействие материалов с газовыми потоками. - М.: Машиностроение, 1976. - 224 с.
4. Полежаев Ю.В., Юревич Ф.Б. Тепловая защита. - М.: Энергия, 1976, - 392 с.
5. Kuntz D., Hassan B., Potter D. Predictions of Ablating Hypersonic Vehicles Using an Iterative Coupled Fluid/Thermal Approach // Journal of Thermophysics and Heat Transfer, 2001, vol. 15, no. 2, pp. 129 - 139.
6. Zhluktov S.V., Abe T. Viscous Shock Layer Simulation of Air Flow Past Ablating Blunt Body with Carbon Surface // Journal of Thermophysics and Heat Transfer, 1999, vol. 13, no. 1, pp. 50 - 59.
7. Jones R., Weinstein S., Wilcox B., Yeomans D. NASA ISAS Collaboration on the ISAS MUSES C Asteroid Sample Return Mission. The Jet Propulsion Laboratory, Pasadena California, USA, URL: https://trs.jpl.nasa.gov/bitstream/handle/2014/19211/98-0488.pdf?sequence=1
8. Ageev A.G., Kavyrshin D.I., Sargsyan M.A., Gadzhiev M.Kh., Chinnov V.F. Determination of graphite sublimation rate in high enthalpy plasma flow using 'laser knife' method // International Journal of Heat and Mass Transfer, 2017, vol. 107, pp. 146 -153.
9. Chinnov V.F., Tyuftyaev A.S., Kavyrshin D.I., Ageev A.G., Sargsyan M.A., and Gadzhiev M.Kh. Comprehensive Study of the Effect of Plasma Stream on Heat-Resistant Materials // High Temperature, 2018, vol. 56, no.1, pp. 25 - 32.
10. Аксёнов А.А. FlowVision: Индустриальная вычислительная гидродинамика // Компьютерные исследования и моделирование. 2017. Т. 9. № 1. С. 5 - 20.
11. Ashraf A., Yaqub K., Javeed S., Zeeshan S. Sublimation of graphite in continuous and pulsed arcdischarges // Turkish Journal of Physics, 2010, vol. 34, pp. 33 - 42.
12. Hirt C.W., Nicholls B. D. Volume of fluid (VOF) method for the dynamics of free boundaries // Journal of Computational Physics, 1981, vol. 39, pp. 201.
13. Aksenov A.A., Zhluktov S.V., Savitskiy D.V., Bartenev G.Y., Pokhilko V.I. Simulation of 3D flows past hypersonic vehicles in FlowVision software // Journal of Physics: Conference Series, 2015, vol. 653, no. 012072, available at: http://iopscience.iop.org/1742-6596/Y2015.
14. Aksenov A.A., Gudzovsky A.V., Serebrov A.A. Electrohydrodynamic Instability of Fluid Jet in Microgravity // Proc. of 5th Int. Symposium on Computational Fluid Dynamics (ISCFD), Aug. 31 - Sept. 3, 1993, Sendai, Japan, 1993, vol.1, pp 19 - 24.
15. Сушко Г.Б., Харченко С.А. Экспериментальное исследование на СКИФ МГУ "Чебышев" комбинированной MPI+threads реализации алгоритма решения систем линейных уравнений, возникающих во FlowVision при моделировании задач вычислительной гидродинамики // Труды международной научной конференции Параллельные вычислительные технологии. (ПаВТ'-2009), (Нижний Новгород, 30 марта - 3 апреля 2009). - Челябинск: Изд-во ЮУрГУ, 2009. C. 316 - 324.
16. Харченко С.А. Влияние распараллеливания вычислений с поверхностными межпроцессорными границами на масштабируемость параллельного итерационного алгоритма решения систем линейных уравнений на примере уравнений вычислительной гидродинамики // Труды Международной научной конференции "Параллельные вычислительные технологии", (Санкт-Петербург, 28 января - 1 февраля 2008,). - Челябинск: Изд-во ЮУрГУ, 2008. - С. 494 - 499.
17. Гаджиев М.Х., Тюфтяев А.С., Саргсян М.А., Демиров Н.А., Щербаков В.В. Диагностический комплекс для исследования взаимодействия плазменной струи с термостойкими материалами // Вестник Дагестанского государственного университета. Естественные науки. 2016. Т. 31. № 1. С. 22 - 27.
18 Tyuftyaev A.S., Gadzhiev M.Kh., Sargsyan M.A., Chinnov V.F., Demirov N.A., Kavyrshin D.I., Ageev A.G. and Khromov M.A. Research methods of plasma stream interaction with heat-resistant materials // Journal of Physics: Conference Series, 2016, vol. 774, no. 012204, URL: http://iopscience.iop.org/article/10.1088/1742-6596/774/1/012204.
19. Справочные данные по давлению насыщенных паров, URL: http://www.physics.nyu.edu/kentlab/How to/ChemicalInfo/VaporPressure/morepressure
