CC BY
20
Вычислительные технологии Том 19, № 5, 2014
Прямое статистическое моделирование высокоскоростного течения около затупленного тела в смеси CO2/N2 с учётом поверхностной рекомбинации*
A. Н. Молчанова1'2, Е.А. Бондарь1'2, А. В. Кдтттковский1'2, М.С. Иванов1
Проведено численное исследование обтекания затупленного тела высокоэнталь-пийным химически реагирующим потоком для условий, соответствующих входу в атмосферу Марса, с использованием метода прямого статистического моделирования Монте-Карло. Предложена молекулярная модель поверхностной рекомбинации, позволяющая учесть влияние этого процесса на высотную аэротермодинамику. Исследованы влияние эффектов разреженности и область применимости континуального подхода сравнением результатов прямого статистического моделирования с решением уравнений Навье — Стокса для разреженного и околоконтинуального режимов течения. Показано существенное влияние поверхностной рекомбинации на аэротермодинамику тела в рассмотренном диапазоне параметров. Установлено, что для околоконтинуального течения (Кп ~ 0.002) уравнения Навье — Стокса хорошо предсказывают течение в ударном слое и распределение теплового потока по поверхности тела, а в разреженном режиме (Кп ~ 0.03) континуальный подход неверно описывает структуру течения, завышает скорость диссоциации и предсказывает заниженную величину теплового потока на поверхность.
Ключевые слова: аэротермодинамика космических аппаратов в атмосфере Марса, гетерогенная рекомбинация, метод прямого статистического моделирования Монте-Карло.
Введение
На начальных стадиях разработки перспективных космических аппаратов (КА) необходимы исследования их аэротермодинамических характеристик вдоль траектории входа в атмосферу. Характерная особенность высокоскоростных течений около КА состоит в существенной термохимической неравновесности потока. С технической точки зрения экспериментальное моделирование таких условий затруднительно, поэтому основным методом исследования высотной аэротермодинамики КА является численное моделирование.
Существуют два основных подхода к численному моделированию течений: сплош-носредный (континуальный) и кинетический. Континуальный подход, основанный чаще всего на уравнениях Навье — Стокса, не применим для больших высот полета, когда длина свободного пробега молекул сопоставима с размером обтекаемого тела. Для
1 Институт теоретической и прикладной механики СО РАН, Новосибирск, Россия
2Новосибирский национальный исследовательский государственный университет, Россия Контактный e-mail: shumakova@itam.nsc.ru
* Работа выполнена при поддержке СО РАН (Междисциплинарные интеграционные проекты № 39 и № 47), Правительства РФ (Договор № 14.Z50.31.0019) и РФФИ (гранты № 14-08-01252 и № 14-0700065).
моделирования таких течений следует использовать кинетический подход, наиболее развитой и эффективной реализацией которого в настоящее время является метод прямого статистического моделирования (ПСМ) [1]. В Институте теоретической и прикладной механики СО РАН создана программная система SMILE+—+ [2] для численного исследования высотной аэротермодинамики космических аппаратов методом ПСМ, позволющая проводить полный цикл численного исследования на многопроцессорных ЭВМ начиная от создания трёхмерной геометрической модели КА и заканчивая визуализацией результатов математического моделирования.
Одна из актуальных задач численного моделирования высотной аэротермодинамики КА состоит в учёте физико-химических процессов на поверхности аппарата. В высокотемпературном вязком ударном слое у поверхности КА протекают реакции диссоциации молекул набегающего потока на атомы и более простые молекулы, в результате которых происходит поглощение тепла. Однако поверхность КА является катализатором для обратных реакций рекомбинации, происходящих с выделением тепла, что может в разы увеличить тепловой поток на поверхность аппарата [3]. Отметим, что поверхностные химические процессы изучены гораздо хуже, чем газофазные, и степень неопределённости, связанная с моделированием этих процессов, существенно выше. В то время как для методов континуальной газовой динамики был предложен ряд моделей поверхностных химических процессов, различающихся сложностью реализации, детальностью и точностью, моделирование каталитических реакций в методе ПСМ практически не проводилось. В литературе представлено небольшое число публикаций, посвящённых данному вопросу [4-6], при этом влияние гетерогенных физико-химических процессов на высотную аэротермодинамику космических аппаратов на основе метода ПСМ пока не исследовалось.
Настоящая работа посвящена исследованию аэротермодинамики затупленного тела (цилиндр со скругленными кромками) в высокоэнтальпийном разреженном потоке смеси CO2/N2 на основе метода прямого статистического моделирования. В частности, для учёта эффектов каталитичности поверхности проводятся разработка, реализация в системе SMILE+—+ и верификация молекулярной модели поверхностной рекомбинации для использования в методе ПСМ. Рассматриваются условия, подробно исследованные в рамках проекта SACOMAR (Key Technologies for Safe and Controlled Martian Entry) седьмой рамочной программы Европейского Союза (FP7) на основе континуального подхода [7]. Изучается влияние каталитических свойств поверхности на структуру вы-сокоэнтальпийного разреженного CO2/N2-течения и тепловой поток к поверхности КА. Область применимости традиционных методов вычислительной газовой динамики, основанных на уравнениях Навье — Стокса, оценивается путём сравнения с результатами, полученными методом ПСМ с использованием системы SMILE+—+. Перекрёстная валидация кинетического и континуального подходов выполняется для условий, соответствующих умеренной высоте, на которой эффекты разреженности незначительны. В условиях, соответствующих большой высоте, исследуется влияние эффектов разреженности и термической неравновесности сравнением результатов, полученных с использованием кинетического и континуального подходов.
1. Моделирование гетерогенной рекомбинации
При обтекании космического аппарата на его поверхности протекают гетерогенные физико-химические процессы: атомы газа, попадая на поверхность тела, адсорбиру-
ются и затем могут рекомбинировать с образованием молекулы и последующей десорбцией. При этом рекомбинация проходит либо за счёт соударения падающего атома с адсорбированным (механизм Или — Ридила), либо путём взаимодействия двух адсорбированных на поверхности атомов (механизм Ленгмюра — Хиншельвуда).
Процессы гетерогенной рекомбинации в континуальной динамике учитываются в граничных условиях следующим образом [3,8,9]. В единицу времени на стенку поступает поток частиц ¿-го сорта
Р у 2птг'
где сг,тг — массовая концентрация и масса частиц, р,Т — плотность и температура потока, к — постоянная Больцмана. Пусть доля от числа падающих частиц вступает в гетерогенные реакции, что в установившемся течении компенсируется диффузионным потоком, а Зг — диффузионный поток частиц данного сорта. Согласно первому закону Фика
дс'
3 = -ри дП' •
Тогда на стенке выполняется условие
п дС' I кТ
- 3 = рА— = 7грс*
дп V 2птг
В этом случае закон сохранения массы данного компонента на стенке примет вид
п дСг / кТ
' дп V 2пт.
г
Уравнение (1) является граничным условием для концентраций компонентов ¿-го сорта на стенке.
Величину 7г называют вероятностью рекомбинации. Вероятность рекомбинации зависит от сорта частиц, от материала поверхности, на которую они попадают, и от температуры этой поверхности. Предполагается, что вероятность рекомбинации не зависит от давления, так как данный процесс на поверхности считается реакцией первого порядка.
При учёте поверхностных процессов в методе ПСМ параметры модели должны определять взаимодействие каждой отдельно взятой молекулы газа с поверхностью. Ниже приведена реализованная в системе БМИЕН—+ молекулярная модель поверхностных реакций, параметром которой является вероятность рекомбинации 7г. При этом значение вероятности рекомбинации — константа, одинаковая для всех частиц одного сорта. Для учёта взаимодействия частиц с поверхностью без химических превращений в работе используется диффузный закон отражения с полной аккомодацией энергии.
Рассмотрим процесс реакции рекомбинации на поверхности:
А + В —► АВ + АЕ, (2)
где А, В — некоторые модельные частицы, которые могут быть как атомами, так и молекулами, АВ — молекула, АЕ — энергия, выделяющаяся в результате реакции.
Поверхность в системе БМИЕН—+ представлена набором ячеек, которые в трёхмерном случае представляют собой треугольники, а в двухмерном — отрезки (рис. 1).
Рис. 1. Схема алгоритма для частицы А
Для каждой ячейки задан набор вероятностей рекомбинации 7А и 7В (для компонентов А и В соответственно) и коэффициентов аккомодации энергии рекомбинации вАВ (0 < 7а , 7В, вАВ < 1) для всех возможных реакций рекомбинации на поверхности вида (2). Вероятности рекомбинации 7А и 7В и коэффициенты аккомодации вАВ зависят лишь от материала поверхности и её температуры. Алгоритм взаимодействия частицы с поверхностью, описанный ниже, применяется к каждой ячейке поверхности на каждом временном шаге.
Пусть модельная частица А сталкивается с некоторой ячейкой поверхности. Тогда она с вероятностью 7А выбирается для реакции А + В —> АВ на поверхности или с вероятностью 1 — 7а отражается от поверхности согласно обычному закону отражения без рекомбинации. Частица, выбранная для реакции на поверхности, может либо вступить в реакцию рекомбинации с модельной частицей В, либо адсорбироваться поверхностью ячейки. Первый вариант возможен, если ячейка уже содержит частицы В, которые были адсорбированы на предыдущих временных шагах. Тогда частица А ре-комбинирует с одной из этих частиц, при этом доля вАВ выделяющейся энергии АЕ поглощается поверхностью, а оставшаяся часть (1 — вАВ)АЕ переходит к продукту реакции — молекуле АВ, которая тут же покидает поверхность. Если ячейка не содержит адсорбированных частиц В, реализуется второй вариант. Заметим, что задано также целое число С (ёмкость ячейки), ограничивающее максимальное число частиц любого
компонента, которые могут быть адсорбированы поверхностью одной ячейки (влияние данного параметра будет рассмотрено ниже). Это означает, что если для частицы A, выбранной для реакции на поверхности, реализуется второй вариант, но количество частиц A, уже адсорбированных на предыдущих временных шагах, равно C, то данная частица отражается от поверхности без рекомбинации согласно обычному закону отражения. Алгоритм для частицы B аналогичен. Схема алгоритма для частицы A приведена на рис. 1.
Таким образом, в рассматриваемой модели поверхностная рекомбинация происходит по механизму, аналогичному рекомбинации Или — Ридила, когда в реакцию вступают падающая газофазная частица и частица, адсорбированная поверхностью. При этом адсорбция происходит по упрощённому механизму с вероятностью, равной вероятности рекомбинации.
Здесь можно выделить два предельных случая:
— некаталитическая поверхность, для которой ya = Yb = 0, т. е. на ней реакция рекомбинации (2) не происходит вообще;
— полностью каталитическая поверхность, для которой ya = yb = 1, т. е. все частицы A и B, приходящие на поверхность, вступают в реакцию рекомбинации (2).
Отметим, что в общем случае коэффициенты ya и Yb зависят друг от друга и не всегда идеальная ситуация, когда рекомбинируют все приходящие на поверхность частицы, может быть достижима. Например, если рассматривается реакция на поверхности, в которой реагенты A и B различны, а поток частиц одного сорта (например A) на стенку больше, чем другого (B), то вступить в реакцию на поверхности смогут не все частицы A, так как для них просто не хватит партнеров по рекомбинации частиц B. Далее под полностью каталитической понимается поверхность, на которой рекомби-нируют все частицы, для которых находится партнер по рекомбинации. Отметим, что если на поверхности рассматриваются реакции с совпадающими реагентами (A = B), то такой проблемы не возникает, и в случае полностью каталитической стенки на ней рекомбинируют все приходящие частицы.
2. Модели для описания столкновительных процессов
Моделирование столкновений в системе SMILE+—+ основано на схеме мажорантной частоты метода ПСМ [10]. Рассмотрена восьмикомпонентная смесь газов: C, O, N, CO, NO, O2, N2, CO2. Для описания столкновительных процессов в газовой фазе использованы следующие модели. Для упругих столкновений была применена модель переменных мягких сфер (Variable Soft Sphere — VSS), предложенная Коурой и Мацумо-то и подробно описанная в работе [1]. Для описания поступательно-вращательного и поступательно-колебательного обмена использована модель Ларсена — Боргнакке (LB) [11] для дискретной внутренней энергии с зависящими от температуры вращательными [12] и колебательными [13, 14] релаксационными числами. Для всех двухатомных молекул применялись модели мягкого ротатора для представления вращений и ангармонического осциллятора с параметром ангармоничности 0.82 для представления колебаний. Для описания вращательной энергии CO2 использовалось приближение твёрдого ротатора. Все три колебательные моды CO2 моделировались как гармонические осцилляторы со значениями характеристических колебательных температур 959.66, 1918.6 и 3382 К соответственно для изгибной, симметричной и антисимметричной мод. Вырож-
Таблица 1. Реакции в газовой фазе
Реакции диссоциации СО2 + М -^ СО + О + М СО + М —► С + О + М МО + М —> М + О + М О2 + М О + О + М М2 + М —> М + М + М
Реакции обмена СО + О ^ О2 + С N2 + О <—► МО + N МО + О <—> О2 СО2 + О <—► О2 + СО
денная изгибная мода моделировалась как два идентичных ангармонических осциллятора, обладающих одинаковой колебательной энергией.
Для моделирования реакций в газовой фазе использовался механизм [15], содержащий 40 реакций диссоциации и восемь реакций обмена (табл. 1, здесь М — любой из восьми компонентов газа). Для описания реакций диссоциации и обмена применялась модель полной столкновительной энергии (ТСЕ), модифицированная для дискретной внутренней энергии. Вращательно-вращательный и колебательно-колебательный энергообмен в данном исследовании не рассматривался. При моделировании столкновений применялась методология, исключающая столкновения с одновременной релаксацией более двух энергетических мод [16,17].
3. Обсуждение результатов 3.1. Влияние параметров модели
Представленная в настоящей работе модель гетерогенной рекомбинации имеет несколько параметров: вероятность рекомбинации, коэффициент аккомодации энергии рекомбинации и ёмкость ячеек. Вероятность рекомбинации и коэффициент аккомодации являются физическими параметрами, зависящими от входных условий задачи (материала поверхности и её температуры), поэтому необходимо исследовать, как влияют эти параметры на само течение. Ёмкость ячейки — технический параметр, ограничивающий число частиц какого-либо сорта, которые могут быть адсорбированы поверхностью одной ячейки. Если данный параметр не ввести, то возможны ситуации, когда число адсорбированных частиц в процессе расчёта будет неограниченно расти. Если же этот параметр выбрать слишком малым, то доля рекомбинирующих частиц может быть меньше заданной. С точки зрения программной реализации ёмкость ячейки представляет собой размер массива, в который записываются адсорбированные поверхностью ячейки частицы какого-либо сорта. Ниже оценивается влияние этого параметра и определяется его оптимальная величина, при которой можно избежать указанных выше ситуаций и которая используется далее при проведении расчётов.
Влияние ёмкости ячейки С. Прежде всего определим, в каких ситуациях ёмкость ячейки будет оказывать влияние на рассматриваемый процесс. В случае такой реакции на поверхности, когда реагенты совпадают, т. е. А + А —у 2А + АЕ, в процессе моделирования на поверхности каждой ячейки будет адсорбироваться не более одной частицы сорта А, так как следующая пришедшая на эту ячейку частица сорта А либо отразится от неё (если вероятность рекомбинации меньше 1), либо вступит в реакцию с уже адсорбированной частицей. Таким образом, ёмкость ячейки влияет только в том случае, если реагенты в реакции различны. В случае, когда на поверхности реагируют частицы
разного сорта и при этом одного реагента больше, чем другого, т. е., например, реакция на поверхности имеет вид А + В —у АВ + ДЕ и поток частиц А на поверхность больше, чем частиц В, адсорбироваться поверхностью в ожидании более редких частиц В будут более частые частицы А. Тогда ёмкость ячейки при её величине больше нуля на рассматриваемый процесс оказывать влияния не будет. Если же в реакции участвуют частицы разных сортов и потоки частиц этих сортов на поверхность примерно равны, то возможны ситуации, когда некоторые пары частиц, которые должны были прореагировать (согласно заданной вероятности рекомбинации), отразятся от поверхности, не вступив в реакцию, так как массив частиц, адсорбированных поверхностью данной ячейки, оказался полностью заполнен. Таким образом, емкость ячейки влияет на работу алгоритма только в случае реакций с различными частицами (А и В не совпадают) и при условии, когда потоки этих частиц на поверхность примерно равны.
Таким образом, необходимо определить, какая величина ёмкости ячейки гарантирует учёт всех поверхностных реакций, которые должны произойти согласно выбранным значениям вероятности рекомбинации. Для исследования этого вопроса были проведены расчёты двухмерного обтекания затупленной пластины. Состав течения — 50 % СО, 50% О, поверхность полагалась полностью каталитической, на поверхности рассматривалась реакция СО + О —> СО2. Чтобы не нарушить состав приходящего на поверхность течения, реакций в газовой фазе не происходило. На рис. 2, а представлено влияние коэффициента С на тепловой поток к поверхности (здесь QNC тепловой поток в некаталитическом случае, Qc — тепловой поток в полностью каталитическом случае с ёмкостью ячейки в алгоритме, равной С). Видно, что с увеличением С от 1 до 8 тепловой поток возрастает, но при С > 8 величина теплового потока не меняется. Это указывает на то, что при С = 8 все возможные реакции на поверхности уже учтены. Во всех дальнейших расчётах и была использована величина ёмкости С = 8.
Влияние вероятности рекомбинации 7 и коэффициента аккомодации энергии рекомбинации в. Для исследования влияния вероятности рекомбинации 7 была проведена серия расчётов осесимметричного течения около цилиндра со скруглёнными кромками (модель цилиндра и поле поступательной температуры представлены на рис. 3). Параметры течения представлены в табл. 2 (условия II). На поверхности рассматривались следующие реакции:
СО + О СО2, К + К —► Ы2. (3)
Таблица 2. Параметры течения
Параметр Условия I Условия II
Давление, Па 22 0.96
Температура, К 222.92 155.2
Плотность, кг/м3 5.26 ■ 10-4 3.207 ■ 10-5
Скорость, м/с 4072 5941.2
Число Маха 17.32 30.17
Энтальпия, МДж/кг 8.5 18
Число Кнудсена
(по диаметру цилиндра) 1.89 ■ 10-3 2.68 ■ 10-2
Массовая доля С02 0.97 0.97
Массовая доля N2 0.03 0.03
Рис. 2. Влияние параметров модели: а — зависимость теплового потока в точке торможения от ёмкости ячейки С в полностью каталитическом случае (7 = 1), QC = (Я с — Яыс)/Яыс; б — массовая доля С02 вдоль линии торможения при разных значениях вероятности рекомбинации 7; в — зависимость поступательной температуры потока в точке торможения Т(в) от коэффициента аккомодации энергии рекомбинации в в полностью каталитическом случае (7 = 1), Т* = (Т(в) — Т(1))/Т(1); г — зависимость теплового потока в точке торможения Я от вероятности рекомбинации 7 (при в = 1) и коэффициента аккомодации энергии рекомбинации в (при 7 =1) для случая некаталитической поверхности, = (Я7 — Яыс)/Яыс,
а
в
г
Чтобы оценить влияние каталитических реакций на аэротермодинамику КА, в настоящей работе в основном рассматриваются два предельных случая: некаталитическая поверхность, на которой не происходят реакции рекомбинации (700 = То = TN = 0), и полностью каталитическая поверхность (7с0 = Т0 = TN = 1). Фактически при данных условиях последний случай совпадает со случаем так называемой суперкаталитической
Рис. 3. Поле температуры, условия II (см. табл. 2)
поверхности, на которой концентрации компонентов стремятся вернуться к своим исходным значениям в набегающем потоке. Реакция О + О —> О2 авторами не исследовалась.
Расчёты были выполнены с различными промежуточными значениями коэффициента 7 = 7со = То = TN: 7 = 1, 7 = 0.5, 7 = 0.25, 7 = 0.1, 7 = 0.01, 7 = 0.001 и 7 = 0. Коэффициент аккомодации химической энергии полагался равным единице (в = 1), т.е. энергия, выделяющаяся при каталитической реакции, полностью поглощалась поверхностью. На рис. 2, г сплошной линией приведена зависимость теплового потока в точке торможения ) относительно некаталитического случая ^N<3) от вероятности рекомбинации 7 = 7СО = 70 = 7N. Видно, что данная зависимость нелинейная: даже при малых значениях коэффициента 7 = 0.1 каталитические реакции увеличивают тепловой поток в точке торможения на 7%, тогда как при 7 = 1 этот вклад составляет около 13 %. На рис. 2, б представлены массовые доли компонента СО2 вдоль линии торможения для различных 7. Здесь также можно наблюдать сильную нелинейную зависимость массовой концентрации продукта каталитической реакции на стенке от вероятности рекомбинации 7: при изменении 7 от 1 до 0.1 падение массовой концентрации СО2 составляет примерно 6%, что почти совпадает с падением массовой концентрации этого компонента на стенке при изменении 7 от 0.1 до 0.
Для исследования влияния коэффициента аккомодации химической энергии в была проведена серия расчётов с теми же параметрами течения, что и в предыдущем случае, но поверхность полагалась полностью каталитической и варьировалось значение в: в = 0, в = 0.25, в = 0.5, в = 0.75, в = 1. На рис. 2, г пунктирной линией представлена зависимость теплового потока в точке торможения QFC в полностью каталитическом случае относительно некаталитического ^N(3) от коэффициента аккомодации энергии рекомбинации в. Зависимость поступательной температуры Т(в) в точке торможения относительно случая с полной аккомодацией энергии рекомбинации (Т(1), т. е. в = 1) при полностью каталитической стенке представлена на рис. 2, в. Видно, что
зависимость теплового потока и температуры от коэффициента аккомодации энергии рекомбинации в является линейной. Отметим также, что при любом ненулевом фиксированном значении вероятности рекомбинации 7 поступательная температура у стенки аналогично зависит от коэффициента в.
Похожим образом изменение коэффициента 7 влияет на температуру газа у стенки. Расчёты показали, что при в = 1 температура газа у стенки при варьировании коэффициента 7 практически не меняется, так как вся энергия рекомбинации поглощается поверхностью. Однако при в < 1 с варьированием 7 температура газа у стенки меняется, поскольку в данном случае меняется также величина поглощаемой газом энергии, выделяющейся в реакциях рекомбинации на поверхности.
Таким образом, показано влияние параметров модели гетерогенной рекомбинации на предсказываемую величину теплового потока. Отметим, что различным параметрам 7 и в может соответствовать практически одинаковый по тепловому потоку результат, что является хорошо известным фактом [3]. Интересно, что в описанном режиме течения тепловой поток качественно различным образом зависит от вероятности рекомбинации (нелинейно) и коэффициента аккомодации энергии рекомбинации (линейно).
Ниже при исследовании аэротермодинамики затупленного тела рассмотрены два предельных случая, соответствующих минимальному и максимальному значению теплового потока: некаталитическая стенка (поверхностные реакции не происходят) и полностью каталитическая стенка (происходит максимально возможное число поверхностных реакций с полной аккомодацией энергии).
3.2. Влияние гетерогенных физико-химических процессов на аэротермодинамику затупленных тел
Было проведено исследование влияния гетерогенных физико-химических процессов на аэротермодинамику цилиндра диаметром 50 мм со скруглённой кромкой (радиус закругления 11.5 мм) в высокоэнтальпийных С02/^-течениях (см. рис. 3). Такая форма может обеспечить большой отход головной ударной волны, что удобно при проведении экспериментальных измерений, и является достаточно простой и удобной для численных исследований. Расчёты проводились в осесимметричной постановке. Параметры течений представлены в табл. 2 (см. выше). Условия I и II соответствуют двум точкам траектории входа космического аппарата в атмосферу Марса — 1 и 2 соответственно [7]. Точка 1 соответствует максимуму полного давления, посчитанного по параметрам за фронтом прямого скачка уплотнения, вдоль траектории (высота 40 км) и представляет околоконтинуальные условия обтекания. Точка 2 соответствует большей высоте полета (55 км) и представляет высокоскоростные разреженные условия обтекания. Температура поверхности во всех случаях полагалась равной 300 К. Расчёты выполнялись на вычислительном кластере ИТПМ СО РАН, было задействовано от 35 до 96 процессоров. В вычислениях использовано до 40 млн модельных частиц.
Для условий I рассматривалась расчётная область размером 30 х 40 мм (—15 мм < х < 15мм, 0 < у < 40мм), прямоугольная однородная сетка для вычисления газодинамических макропараметров содержала 120 х 160 ячеек. Для условий II расчётная область имела размеры 400 х 200мм (—120мм < х < 280мм, 0 < у < 200мм), а сетка содержала 400 х 200 ячеек. В обоих случаях начало координат совпадало с точкой торможения, а для моделирования столкновений использовалась двухуровневая прямоугольная равномерная сетка. Первый уровень совпадал с сеткой для вычисления мак-
ропараметров. Прямоугольная сетка второго уровня автоматически адаптировалась в процессе расчёта таким образом, чтобы размер ячейки был в три раза меньше локальной длины свободного пробега. Также в обоих случаях передняя поверхность цилиндра содержала 27 поверхностных ячеек, скругленная кромка — 15, боковая поверхность — 80 и задняя поверхность — 20.
На поверхности моделировались две реакции рекомбинации (3). Для иллюстрации влияния гетерогенной рекомбинации было рассмотрено два типа поверхности: некаталитическая (тс0 = 70 = YN = 0) и полностью каталитическая (7С0 = 70 = 7N = 1). Коэффициенты аккомодации энергии рекомбинации полагались равными вС02 = в^ = 1.
Условия I соответствуют околоконтинуальному режиму течения и могут быть использованы для сравнения и перекрестной валидации кинетического и континуального подходов. Результаты для данных условий представлены на рис. 4, I. На рис. 4, I, а показаны поступательная, вращательная и колебательная температуры вдоль линии торможения, направление потока — слева направо, начало координат соответствует точке торможения. В набегающем потоке (левая часть графика, х < —8 мм) температуры различных мод совпадают. Далее по мере приближения к точке торможения, во фронте ударной волны, становится заметным рост поступательной, а потом вращательной и колебательной температур. В зоне фронта ударной волны наблюдается
Некаталитиче ская
--- С02
-х- СО
Полностью 0 каталитическая ---------------------- С02
-0.01
С0 0
-0.005
2000 г 2, кВт/м2
1000
Полностью " каталитическая
Некаталитическая
X, м
0.005 0.01 0.015 0.02 0.025
г, м
Ь
II
Некаталитиче ская С02 СО О
Полностью каталитическая ------------ С02
- СО
............ О
2000 г 2, кВт/м2
0.05 -0.04 -0.03 -0.02 -0.01 0
х, м
1000
Полностью каталитиче ская
- Некаталитиче ская
0.005 0.01 0.015 0.02 0.025
г, м
I
0
0
0
а
в
0
X
X
0
а
в
Рис. 4. Профили поступательной Т, вращательной Тг и колебательной Ту температуры вдоль линии торможения (а), профили массовых концентраций компонентов вдоль линии торможения для некаталитической и полностью каталитической поверхности (б) и распределение теплового потока на поверхности (в). I — условия I, II — условия II
сильная термическая неравновесность: максимальное значение поступательной температуры достигает 12000 К, вращательной — 7000 К, колебательной — 5000 К. Вращательная мода релаксирует быстрее колебательной, поэтому вращательная температура превышает колебательную. Толщина этой зоны составляет около 4 мм. Далее температуры выравниваются и на расстоянии от 4 до 1 мм до точки торможения наблюдается вязкий ударный слой, в котором все три температуры совпадают и происходит их падение примерно с 5000 до 4000 К. В окрестности точки торможения в пограничном слое ^ > — 1 мм) происходит резкое падение температур из-за взаимодействия газа с холодной стенкой.
В силу роста температуры в ударной волне в газовой фазе происходят реакции диссоциации. На рис. 4, I, б приведены массовые доли компонентов СО2, СО и О вдоль линии торможения. Видно, что диссоциация приводит к падению массовой концентрации СО2 до 75 % и росту массовых концентраций СО иОдо15и8% соответственно. Протекание реакций диссоциации объясняет падение температуры в вязком ударном слое. Из сравнения расчётов в некаталитическом и полностью каталитическом случаях следует, что в окрестности точки торможения реакции рекомбинации на поверхности при полностью каталитическом варианте значительно увеличивают концентрацию СО2 и уменьшают концентрации СО и О: значения концентраций компонентов в точке торможения практически совпадают с их значениями в набегающем потоке. Каталитические реакции вносят существенный вклад в тепловой поток на поверхность, что видно из рис. 4,1, в, на котором представлено распределение теплового потока на поверхности цилиндра по радиусу для условий I (см. табл. 2). Профиль теплового потока немонотонный и имеет максимум в точке г = 16 мм. График теплового потока, полученный для полностью каталитической стенки, лежит выше графика для некаталитического случая (отличие в точке торможения составляет около 40%).
Условия II соответствуют разреженному режиму течения. На рис. 4, II, а показаны профили поступательной, вращательной и колебательной температур вдоль линии торможения. Структура течения в этом случае качественно отличается от структуры течения для условий I. Ударная волна фактически неотделима от пограничного слоя, а зона вязкого ударного слоя, наблюдаемая в расчётах для условий I, отсутствует. Во фронте ударной волны наблюдается сильная термическая неравновесность течения, которая имеет значительно большую толщину по сравнению со случаем I: поступательная температура начинает расти на расстоянии 5 см до точки торможения. Поступательная, вращательная и колебательная температуры выравниваются только в окрестности точки торможения, так как частицы отражаются от поверхности с одинаковыми во всех модах энергиями, соответствующими температуре поверхности Tw.
Газ за головной ударной волной нагревается выше 10000 К — при таких условиях происходит диссоциация молекул СО2. На рис. 4, II, б приведены профили массовых долей СО2, СО и О вдоль линии торможения. Отметим, что при данных условиях скорость диссоциации несколько ниже, чем при условиях I: реакции диссоциации в газовой фазе приводят к уменьшению концентрации СО2 до 80 % и увеличению концентраций СО и О до 10 и 5% соответственно. В случае полностью каталитической стенки реакции рекомбинации на поверхности приводят к значительному увеличению концентрации СО2 и уменьшению концентраций СО и О около точки торможения. Значения массовых долей компонентов СО2, СО и О в точке торможения в случае полностью каталитической стенки близки к их значениям в набегающем потоке. На рис. 4, II, в представлено распределение теплового потока по поверхности цилиндра как функция
от его радиуса. Как и в предыдущем случае, профиль теплового потока немонотонный и имеет максимум в точке r =16 мм, а график, полученный с учётом поверхностной рекомбинации, лежит выше графика для некаталитической стенки (отличие составляет более 13 % в точке торможения).
На рис. 5 представлены результаты вычислений на основе континуального подхода с использованием различных моделей и расчётных программ [7] в сравнении с результатами, полученными методом ПСМ, для условий I и II. В работе [7] были использованы следующие расчётные программы.
• Программа TAU Немецкого аэрокосмического центра DLR. В программе TAU решаются нестационарные уравнения Навье — Стокса для многокомпонентной реагирующей газовой смеси на неструктурированных сетках. Применяемые в программе модели и алгоритмы описаны в [18,19], а также на сайте http://tau.dlr.de.
• Программа, разработанная в ЦНИИмаш. Данная программа основана на решении нестационарных уравнений Навье — Стокса для реагирующей многокомпонентной смеси на структурированных многоблочных сетках (информация об используемых моделях и алгоритмах представлена в [20]).
• Коммерческий пакет ANSYS Fluent [21], применяемый в расчётах сотрудниками Итальянского центра аэрокосмических исследований CIRA. В данных расчётах также решались нестационарные уравнения Навье — Стокса для многокомпонентной реагирующей смеси.
12000гT K
8000-
4000-
I
10"
2000 г Q, кВт/м:
1000-
DLR
TsNIImash
CIRA
DSMC
-0.01
-0.005
-0.01
-0.005
0.01
0.02
-T, K
20000
10000
DLR - TsNIImash CIRA DSMC
-0.05 -0.04 -0.03 -0.02 -0.01 0
x, м
100
§
я o o
lio-2
10
б II
DSMC
-■--CO2
-*--CO
-с--O
* г CIRA
-0.03 -0.02
б
2000 г Q, кВт/м2
-0.01 0
x, м
1000.=-
DLR TsNIImash CIRA DSMC
0.01
в
0.02
Рис. 5. Сравнение с расчётами на основе уравнений Навье — Стокса: а — температура, б — массовые доли, некаталитическая стенка, в — тепловой поток, полностью каталитическая стенка. I — условия I, II — условия II
0
0
0
0
X. м
X. м
r. м
а
в
0
0
r. м
а
Во всех трёх указанных программах течение моделируется на основе уравнений На-вье — Стокса, которые в настоящее время являются основными при исследовании аэротермодинамики космических аппаратов (см. [22], с. 365). Отметим, что в континуальной гиперзвуковой аэротермодинамике используется также ряд других моделей, таких как приближение вязкого ударного слоя или параболизованные уравнения Навье — Стокса (см., например, [9]). По мнению авторов, сравнение результатов расчётов на основе кинетического подхода с этими моделями также представляет интерес и планируется в будущем.
Отметим, что физико-химические модели, использованные в расчётах методом ПСМ, хорошо согласованы с моделями, применяемыми в континуальных расчётах перечисленными выше программами. В частности, 1 — в представленных расчётах методом ПСМ использован тот же набор газофазных реакций диссоциации и обмена [15], что и в континуальных расчётах [7], 2 — модель полной столкновительной энергии обеспечивает в условиях термического равновесия скорости реакций, совпадающие с аррениусовыми скоростями, используемыми при решении уравнений Навье — Стокса, 3 — модель Ларсе-на — Боргнакке, применяемая для описания поступательно-колебательного энергообмена, обеспечивает макроскопическую скорость колебательной релаксации, соответствующую данным Милликена — Уайта [13] с высокотемпературной поправкой Парка [14], использованым в континуальных расчётах [7].
Такое сравнение позволяет провести перекрестную валидацию континуального и кинетического подходов и оценить влияние эффектов разреженности на аэротермодинамику затупленного тела в атмосфере Марса. Отметим, что во всех континуальных расчётах использована модель однотемпературной кинетики (без учёта колебательной релаксации и колебательно-диссоциационного взаимодействия) и не учтены релаксационные процессы введением коэффициента объёмной вязкости, в то время как учёт данных процессов может заметно повлиять на результаты расчётов.
Для околоконтинуальных условий течения (условия I, соответствующие высоте 40 км) результаты, полученные с помощью метода ПСМ и на основе уравнений Навье — Стокса, в целом хорошо согласуются. На рис. 5, I, а представлены графики поступательной температуры, полученные методом ПСМ и решением уравнений Навье — Стокса. Видно их различие при x < —4 мм: график ПСМ имеет ярко выраженный пик и находится выше всех профилей, полученных на основе уравнений Навье — Стокса. Последнее объясняется тем, что метод ПСМ в отличие от приведенных расчётов на основе уравнений Навье— Стокса учитывает термическую неравновесность течения, которая присутствует в зоне ударной волны, где температуры различных мод существенно отличаются друг от друга (см. рис. 4, I, а): медленная релаксация вращательной и колебательной энергий приводит к очень большим значениям поступательной температуры. В зоне вязкого ударного слоя (при x > —4мм), где наступает термическое равновесие, результаты расчётов методом ПСМ хорошо совпадают с континуальными: температура для всех расчётов сначала становится примерно одинаковой и затем падает также с одинаковой скоростью.
Из графиков массовых долей (в логарифмической шкале) также видно небольшое различие в результатах, полученных континуальным и кинетическим подходами, в зоне ударной волны: в методе ПСМ реакции диссоциации начинаются раньше, чем при решении уравнений Навье — Стокса (рис. 5Д, б). Однако в ударном слое и в окрестности точки торможения массовые доли, предсказываемые разными подходами, согласуются достаточно хорошо для всех компонентов. По всей видимости, учёт объёмной вязко-
сти в уравнениях Навье — Стокса позволил бы добиться ещё лучшего согласия с методом ПСМ. Также лучшего согласия (особенно по профилям температуры) следовало бы ожидать и при использовании в уравнениях Навье — Стокса двух- или многотемпературной кинетики. На рис. 5, I, в представлены графики распределения теплового потока по поверхности в полностью каталитическом случае. Весь график, полученный методом ПСМ, находится между графиками, рассчитанными на основе уравнений На-вье — Стокса (отклонение составляет не более 5 %), и почти полностью совпадает с кривой, полученной с использованием программы ЦНИИмаш. Таким образом, результаты расчётов методом ПСМ находятся в пределах неопределённости, связанной с моделированием на основе уравнений Навье — Стокса.
Для разреженных условий течения (условия II, соответствующие высоте 55 км) зона различия результатов, полученных на основе метода ПСМ и уравнений Навье — Стокса, гораздо шире и фактически продолжается вплоть до самого цилиндра. Из рис. 5, II, а видно, что профиль поступательной температуры, рассчитанный методом ПСМ, существенно превышает аналогичные профили, полученные на основе уравнений Навье — Стокса, что связано с большей степенью термической неравновесности, которую предсказывает кинетический подход (на рис. 4, II, а температуры различных молекулярных мод для метода ПСМ различаются в разы вдоль большей части течения начиная от переднего фронта ударной волны, х = —40 мм, и заканчивая окрестностью точки торможения, х = —5мм). Такие существенные различия в структуре течения приводят к заметному различию в скоростях протекания газофазных реакций. Действительно, графики массовых долей в некаталитическом случае (рис. 5, II, б) для метода ПСМ и уравнений Навье — Стокса сильно различаются: уравнения Навье — Стокса предсказывают большую степень диссоциации молекул С02, при этом диссоциация в методе ПСМ начинает протекать гораздо раньше в набегающем потоке. Как и для околоконтинуальных условий, использование более сложной двух- или многотемпературной кинетики и учёт объёмной вязкости в уравнениях Навье — Стокса, возможно, позволили бы несколько улучшить согласие результатов между рассматриваемыми подходами.
В случае полностью каталитической поверхности происходит выравнивание графиков в окрестности точки торможения: значения массовых долей компонентов на стенке для метода ПСМ и уравнений Навье — Стокса совпадают друг с другом (данный результат на графике не представлен). Это происходит потому, что при таких условиях на стенке (суперкаталитическая поверхность) все частицы, образовавшиеся вследствие диссоциации молекул в газовой фазе, на поверхности вступают в реакцию рекомбинации с образованием молекулярных компонентов. Различие двух подходов можно наблюдать также на графиках теплового потока (рис. 5, II, в): график, полученный с помощью метода ПСМ, лежит выше всех графиков, рассчитанных на основе уравнений Навье — Стокса. В частности, в точке торможения расчёты методом ПСМ предсказывают величину теплового потока, превышающую вычисленную с использованием уравнений Навье — Стокса более чем на 25%. Отметим, что для некаталитического случая между континуальным и кинетическим подходами в величинах теплового потока были получены почти двукратные отличия.
Такие результаты для условий II можно объяснить влиянием эффектов разреженности, которые вызывают не только рост коэффициента теплопередачи, но и, как было показано, снижение скорости диссоциации вследствие термической неравновесности, что в свою очередь приводит к увеличению температуры и дополнительному росту теплового потока. В околоконтинуальном случае условий I влияние эффектов разре-
женности существенно меньше и результаты моделирования на основе континуального и кинетического подходов хорошо совпадают.
Заключение
Численно исследовано обтекание затупленных тел высокоэнтальпийным разреженным потоком для условий, соответствующих входу космического аппарата в атмосферу Марса с учётом поверхностной рекомбинации. Для этого разработана и реализована в системе БМ^ЕН—+, предназначенной для прямого статистического моделирования газовых течений, молекулярная модель гетерогенной рекомбинации, позволяющая учитывать влияние каталитических свойств поверхности на аэротермодинамику космического аппарата. Модель проанализирована в приложении к задаче обтекания тела в атмосфере Марса. В частности, исследована зависимость аэротермодинамических характеристик от физических параметров модели — вероятности рекомбинации и коэффициента аккомодации энергии рекомбинации. Установлено, что тепловой поток и значения массовых концентраций компонентов в точке торможения в рассмотренном режиме течения зависят от коэффициента аккомодации энергии рекомбинации линейно, а от вероятности рекомбинации, напротив, — нелинейно: даже малое значение вероятности рекомбинации оказывает существенное влияние на эти величины. Вклад каталитических реакций в тепловой поток на поверхность в полностью каталитическом случае может достигать 40 %. Несмотря на то что данная модель описывает поверхностные процессы упрощённо, она подходит для реальных аэрокосмических приложений, когда необходимо оценить влияние поверхностной рекомбинации на высотную аэротермодинамику космического аппарата. Вместе с тем необходима дальнейшая разработка более детальных моделей гетерогенной рекомбинации, описывающих различные поверхностные процессы: рекомбинацию Или — Ридила и Ленгмюра — Хиншельвуда, а также адсорбцию и десорбцию.
Метод ПСМ был использован для численного исследования аэротермодинамики цилиндра со скругленной кромкой в двух точках траектории спуска космического аппарата в атмосфере Марса: в точке с околоконтинуальными условиями течения (Кп ^ 0.002) и в точке с переходными условиями течения с высокой степенью разреженности (Кп ^ 0.03). Результаты вычислений методом ПСМ для околоконтинуальных условий близки к результатам, полученным на основе уравнений Навье — Стокса. Получено хорошее согласие в расстоянии отхода головной ударной волны, в массовых концентрациях компонентов и в распределении теплового потока. Все различия находятся в пределах неопределённости результатов численного моделирования на основе уравнений На-вье — Стокса, выполняемого различными расчётными программами. Для разреженных условий метод ПСМ предсказывает структуру течения, качественно отличающуюся от результатов, полученных на основе уравнений Навье — Стокса. Эффекты разреженности приводят к тому, что головная ударная волна не может быть полностью отделена от тела и толщина зоны термической неравновесности (т. е. зоны различия поступательной, вращательной и колебательной температур) сравнима с радиусом цилиндра. Скорость диссоциации, предсказываемая методом ПСМ в таких условиях, значительно ниже, чем в уравнениях Навье — Стокса. Это приводит к тому, что величина теплового потока, полученного методом ПСМ, примерно в полтора раза превышает аналогичную величину, вычисленную с использованием уравнений Навье — Стокса.
Авторами планируются как разработка, верификация и валидация более детальных моделей поверхностной рекомбинации для метода ПСМ, так и подробные исследования
области применимости континуального подхода сравнением результатов ПСМ с более точными континуальными моделями, учитывающими процессы релаксации внутренней энергии и влияние термической неравновесности на скорости физико-химических процессов за фронтом головной ударной волны.
Список литературы
[1] Bird G.A. Molecular Gas Dynamics and the Direct Simulation of Gas Flows. Oxford: Clarendon Press, 1994. 458 p.
[2] Ivanov M.S., Kashkovsky A.V., Vashchenkov P.V., Bondar Ye.A. Parallel object-oriented software system for dsmc modeling of high-altitude aerothermodynamic problems // Proc. of 27th Intern. Symp. on Rarefied Gas Dynamics. Pacific Grove, California, 2010. Vol. 1. P. 211-218.
[3] Ковалёв В.Л. Гетерогенные каталитические процессы в аэротермодинамике. М.: Физ-матлит, 2002. 224 с.
[4] Bergemann F. A detailed surface chemistry model for the DSMC method // Proc. of 19th Intern. Symp. on Rarefied Gas Dynamics. Oxford, 1994. Vol. 2. P. 947-953.
[5] Choquet I. A new approach to model and simulate numerically surface chemistry in rarefied flows // Phys. Fluids. 1999. Vol. 11, No. 6. P. 1650-1661.
[6] Simmons R.S., Lord R.G. DSMC simulation of hypersonic metal catalysis in a rarefied hypersonic nitrogen/oxygen flow // Proc. of 20th Intern. Symp. on Rarefied Gas Dynamics. Beijing, 1996. P. 416-421.
[7] GUlhan A. Technologies for Safe and Controlled Martian Entry (SACOMAR). Final Report, 2012, is. 1.1.
[8] Лунев В.В. Течение реальных газов с большими скоростями. М.: Физматлит, 2007. 760 с.
[9] Anderson J.D. Hypersonic and High-Temperature Gas Dynamics. 2nd Edit. Virginia: AIAA, 2006.
[10] Иванов М.С., РогАЗИНСКИй С.В. Сравнительный анализ алгоритмов метода прямого статистического моделирования в динамике разреженного газа // Журнал вычисл. математики и матем. физики. 1988. Т. 28, № 7. С. 1058-1070.
[11] Borgnakke С., Larsen S.P. Statistical collision model for Monte Carlo simulation of polyatomic gas mixture //J. Comp. Phys. 1975. Vol. 18. P. 405-420.
[12] Parker J.G. Rotational and vibrational relaxation in diatomic gases // Phys. Fluids. 1959. Vol. 2, No. 4. P. 449-462.
[13] Millikan R.C., White D.R. Systematics of vibrational relaxation // J. Chem. Phys. 1966. Vol. 39, No. 12. P. 3209-3213.
[14] Park C. Problems of rate chemistry in the flight regimes of aeroassisted orbital transfer vehicles // Progress in Astronautics and Aeronautics. 1985. Vol. 96. P. 511-537.
[15] Park C., Howe J.T., Jaffe R.L., Candler G.V. Review of chemical-kinetic problems of future NASA missions. II: Mars entries // J. of Thermophys. and Heat Transfer. 1994. Vol. 8, No. 1. P. 9-23.
[16] Haas B.L., Hash D., Bird G.A. et al. Rates of thermal relaxation in direct simulation Monte Carlo method // Phys. Fluids. 1994. Vol. 6, No. 6. P. 2191-2201.
[17] Gimelshein S.F., Gimelshein N.E., Levin D.A. et al. On the use of chemical reaction rates with discrete internal energies in the direct simulation Monte Carlo method // Ibid. 2004. Vol. 16, No. 7. P. 2442-2451.
[18] Schwamborn D., Gerhold T., Heinrich R. The DLR TAU-code: Recent applications in research and industry // ECCOMAS CDF 2006. The Netherlands, 2006.
[19] Mack A., Hannemann V. Validation of the unstructured DLR TAU-code for hypersonic flows // 32nd AIAA Fluid Dynamics Conf. and Exhibit. St. Louis, Missouri, 2002. AIAA. 2002-3111.
[20] Kovalev R.V., Gorshkov A.B., Rudin N.F. et al. Experimental and Numerical Simulation of Martian Entry Conditions. // 9th Intern. Planetary Probe Workshop. Toulouse, France, 2012.
[21] ANSYS FLUENT Theory Guide, Version 14.0. ANSYS Inc., 2011.
[22] Hirschel E.H. Basics of Aerothermodynamics. Progress in Astronautics and Aeronautics. AIAA, Reston, VA, 2005. Vol. 204. 412 p.
Поступила в 'редакцию 23 октября 2013 г., с доработки —18 августа 2014 г.
DSMC modeling of a high-enthalpy CO2/N2 flow around a blunt body with surface recombination _
Molchanova A.N.1'2, Bondar Ye.A.1'2, Kashkovsky A.V.1'2, lIvanov M.ST]
Direct Monte Carlo simulation of a high-enthalpy chemically reacting flow around a blunt body (flat-faced cylinder with rounded edges) was performed for conditions corresponding to entry to Martian atmosphere. At these conditions, surface chemical processes, such as surface molecular recombination, can significantly affect high-altitude aerothermodynam-ics. A molecular surface chemistry model was proposed to take into account the effects of heterogeneous chemical processes on high-altitude aerothermodynamics. The surface chemistry effects, as well as the rarefaction effects along with the domain for applicability of continuum approach, were investigated by comparing the results of Direct Simulation Monte Carlo and solutions of Navier —Stokes equations for rarefied and near-continuum flow regimes. We demonstrate that the effect of surface recombination to body aerother-modynamics is significant for the considered flow regimes. Surface reactions on a fully catalytic surface can increase the stagnation point heat flux up to 40%. The limitation for applicability of continuum approach due to rarefaction effects was demonstrated. For the case of the near-continuum conditions (Kn = 0.002) the Navier —Stokes equations correctly predict the flow in the shock layer and the heat flux distribution. For the rarefied conditions (Kn = 0.03), the continuum approach incorrectly describes the flow structure. It overestimates the dissociation rate, and predicts a lower heat flux to the surface. Keywords: spacecraft aerothermodynamics in Martian atmosphere, heterogeneous recombination, direct simulation Monte Carlo method.
Received 23 October 2013 Received in revised form 18 August 2014
1 Khristianovich Institute of Theoretical and Applied Mechanics SB RAS, 630090, Novosibirsk, Russia
2 Novosibirsk State University, 630090, Novosibirsk, Russia
Corresponding author: Molchanova Alexadra N., e-mail: shumakova@itam.nsc.ru
Aknowlegements: Research was supported by SB RAS Multidisciplinary Integration Grants No. 39, No. 47, Grant of Russian Government No. 14.Z50.31.0019 and RFBR, grant No. 14-08-01252 and No. 14-07-00065.
