МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССОВ ОБТЕКАНИЯ ЗАТУПЛЕННОГО ТЕЛА ВЫСОКОСКОРОСТНЫМ ПОТОКОМ Текст научной статьи по специальности «Химические технологии»

Электронный журнал «Труды МАИ». Выпуск № 78

www.mai.ru/science/trudy/

УДК 004.94, 533.6.011

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

тела высокоскоростным потоком

Быков Л.В.,* Никитин П.В.,** Пашков О. А.***

Московский авиационный институт (национальный исследовательский университет), МАИ, Волоколамское шоссе, 4, Москва, A-80, ГСП-3, 125993, Россия

*e-mail: bykov@mai.ru **e mail: petrunecha@gmail.com ***e -mail: gfon2@yandex.ru

Аннотация

В работе представлена математическая модель, описывающая процессы тепломассообмена, протекающие на поверхности затупленного конуса лобовой части летательного аппарата при полёте в атмосфере с гиперзвуковой скоростью. Модель основана на решении дискретных аналогов уравнений Навье-Стокса на нерегулярной расчётной сетке совместно с уравнениями конвекции-диффузии для каждого компонента газовой смеси и уравнениями модели турбулентности Лангтри и Ментера.

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

Ключевые слова: математическое моделирование, уравнения Навье-Стокса, гиперзвуковой летательный аппарат, газовая динамика, теплообмен, многокомпонентная смесь газов

Введение

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

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

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

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

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

Таким образом, моделирование обтекания тел гиперзвуковым потоком сводится к решению сложной многопараметрической задачи.

Постановка задачи.

Для решения задачи обтекания затупленного тела гиперзвуковым потоком выбран затупленный конус со следующими параметрами: конусность - 6°, радиус затупления - 0,0381м, скорость набегающего потока М<= 25,0. Статические параметры состояния газа в потоке: температура Т< - 265,86 К, давление Р< - 53,85 Па. Массовая концентрация молекул: С(02) = 0,233, С (N2) = 0,767. Геометрия исследуемого тела представлена на рисунке 1.

У, т

0.05 0.10 0.15 0.20 X, т

Рис. 1. Геометрия исследуемого тела.

Расчёт процессов теплообмена на поверхности проводился с учётом каталитической активности поверхности. Были рассмотрены два предельных случая. В первом, поверхность конуса принималась абсолютно каталитической К ^ ю). Во втором, на поверхности конуса формировалось покрытие, обладающее нулевой каталитической активностью К ^ 0). В том и другом случае твердая поверхность принималась химически нейтральной к компонентам набегающего потока, т.е. считалась непроницаемой.

Известно, что при обтекании тела гиперзвуковым потоком в сжатом и пограничном слое, газ становится химически активным [1]. В таком газе реализуются как реакции диссоциации молекул, так и реакции рекомбинации [2,3]. Наличие этих реакций изменяет механизм переноса теплоты и массы в пограничном слое, т.е. интенсифицирует процессы тепло-массообмена между газовым потоком и

Термодинамические свойства

поверхностью тела. Уровень критерия Маха (М), при котором в сжатом и пограничном слоях реализуются эти процессы, определяет газодинамическую природу течения, переводя его из сверхзвукового в гиперзвуковое. Начало такого перехода соответствует течению со значением М > 6. В рассматриваемом в данной работе случае, в условиях гиперзвукового полёта, воздух необходимо рассматривать как многокомпонентный газ (смесь из пяти компонентов О2, N2, О, К, N0), в связи с чем возникает необходимость совместного решения системы уравнений Навье-Стокса с уравнениями Максвелла.

Парциальная плотность каждой 1-ой компоненты воздуха вычисляется с использованием уравнения состояния в виде:

А = ^ (1)

где - парциальное давление ¿-ой компоненты в составе смеси; ^ = R^/Mi - газовая постоянная, Ям - универсальная газовая постоянная, Mi - мольная масса ¿-го компонента. Плотность смеси:

_ Рсм

Рем _ R Т, (2)

см

П

где давление смеси Рсм _ ^ р ,

¿_1

Rем = R^/Mем - газовая постоянная смеси, Мем - мольная масса смеси,

Мсм = Х С • м, (3)

С = -массовая концентрация г-ой компоненты в смеси;

Рсм

Термодинамическая энтальпия смеси:

Г' = ЁС Ч, (4)

г =

г=1

2

где ii = |ерЛ ■ dT- термодинамическая энтальпия г-ой компоненты.

т1

Удельная теплоёмкость Ср,г каждой г-ой компоненты задана по кусочно-линейному закону в виде функции от температуры и давления. В результате средняя удельная теплоёмкость газовой смеси вычисляется с использованием соотношения:

п

с = Х С ■ с . (5)

р см / 1 г р,г 1 V /

г=1

где: Сг - массовая концентрация каждого компонента;

Ср,г - удельная изобарная теплоёмкость каждого компонента.

Теплопроводность Хг каждой г-ой компоненты определяется с использованием соотношения из кинетической теории газов:

; 15 *м

Л =--г

г 4 М г

4 Срг ' Мг +1

15 3

(6)

где:^г - динамическая вязкость г-го компонента; Теплопроводность газовой смеси определяется по формуле:

. ^ хд

Дм=^=1(7)

1=1

где: X - мольная концентрация ¿-го компонента; Х-теплопроводность ¿-го компонента; параметр фу рассчитывался с использованием соотношения:

ф _

1 +

Г Л1/2

К

уК У

/ N1/4

ГМ, л

w, ]

VMw,¿ У

( »^ > Мм,,

8 1 + w,¿

М

_ V w, з _

1/2

(8)

Вязкость каждого ¿-го компонента вычислялась по формуле Сатерленда в

виде:

К _Кг

(Т т/2 Т + 5

норм

Т

V10

Т + 5

(9)

где: ^ог - динамическая вязкость ¿-го компонента при нормальных условиях;

Т - статическая температура, К;

Тнорм - температура при нормальных условиях;

5 - эффективная температура (константа Сатерленда).

Динамическая вязкость газовой смеси вычислялась следующим образом:

ХК

Км

(10)

г_1 V™ х -ф.

где:Х - мольная концентрация компонента ¿.

Следует отметить, что при высоких температурах результаты расчёта по формуле Сатерленда (9) несколько отличаются от известных табличных данных [4]. Однако применение формулы Сатерленда в первом приближении можно считать вполне оправданным.

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

Для вычисления коэффициента бинарной диффузии использовалось модифицированное соотношение Чемпена-Энскога.

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

Химическая кинетика

При решении задачи использовались дискретные аналоги системы уравнений Навье-Стокса для модели вязкой сжимаемой теплопроводной среды, а также уравнение переноса излучения.

С учетом характерных времен течения и протекания химических реакций при моделировании применялась модель неравновесной химии.

Для воздуха при высоких температурах известно довольно большое количество моделей химической кинетики, например, схема Или-Райдила, схема Ленгмюра-Хиншельвуда и др. В рамках данной работы рассматривалась модель, состоящая из пяти основных неравновесных химических реакций, три из которых реализуются с участием третьих тел (М) (Таблица 1).

Таблица 1.Основные неравновесные химические реакции

1 02+М<20+М

2 К2+М<2К+М

3 Ш+М<К+0+М

4 Ш+0<02+К

5

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

д

—угСг ) +У- {ргУСг ) = -У-ёг +4 , (11)

где: gi - диффузионный поток г-го компонента;

ю, - скорость образования г-го компонента в химических реакциях.

Скорость образования г-го компонента в химических реакциях вычисляется с использованием соотношения вида:

®,= М„ ЕЙгг (12)

г =1

где :Мм>,г - мольная масса г-го компонента; Ык - количество химических реакций в расчёте;

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

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

Для неравновесной химической реакции молярная скорость образования (распада) ¿-го компонента в реакции г, записывается в виде:

где Хз,г - мольная концентрация компонента 3 в реакции г(Кмоль/м3); П ),г - показатель степени для реагента3 в реакции г;

V 'з,г - стехиометрический коэффициент для реагента3 в реакции г;

V ''з,г - показатель степени для продукта 3 в реакции г (всегда равен стехиометрическому коэффициенту продукта реакции);

Г - коэффициент, учитывающий влияние третьих тел на скорость реакции; к/,г - константа скорости прямой реакции; къ,г - константа скорости обратной реакции.

В уравнении (13) коэффициент Г, учитывающий влияние третьих тел на скорость химической реакции, вычислялся следующим образом:

(13)

N

Г _Ег,Х

(14)

где: у/,г - эффективность компонента3 в реакции г как третьего тела;

X)■ - мольная концентрация компонента 3.

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

Таблица 2. Эффективность химических компонентов в качестве третьего тела

№ реакции 02 N2 N0 N 0

1 5 2 2 2 25

2 2 5 2 3 5

3 2 2 2 5 5

Константа скорости каждой прямой или обратной реакции г вычислялась по выражению Аррениуса (/- прямая реакция, Ь - обратная реакция):

к/Лг = 4,Ь, гТр! 'Ш (15)

где:Л/ь,г- предэкспоненциальный фактор; в/,ь,г - температурный показатель; Е/,ь,г - энергия активации реакции; Я - универсальная газовая постоянная.

Для вычисления константы скорости каждой реакции применялись эмпирические коэффициенты, приведённые в таблице 3.

Таблица 3. Эмпирические коэффициенты для вычисления констант скорости

реакций

№ реакции * К/.Г Аь, г вь, г Кь, г

1 2.5005е+13 -0.5 4.9365е+08 8.90е+11 -0.44 0.0

2 2.0004е+18 -1.5 9.4177е+08 1.91е+17 -1.57 0.0

3 5.5042е+17 -1.5 6.2782е+08 1.67е+17 -1.52 0.0

4 3.1999е+06 1.0 1.6365е+08 2.67е+07 0.92 2.9486е+07

5 6.8027е+10 0.0 3.1395е+08 2.10е+10 -0.04 0.0

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

ОТ{рК) + ч\у{рК + р)) = V- (к*УТ - X+ Т ■ Р)! + Б„ (16)

ОТ ^ у V I)

В свою очередь, член - источник $>ъ может быть представлен в виде:

^ У

^ (17)

где: ^-энтальпия образования компонента у.

Для решения системы уравнений Навье-Стокса, а также уравнений конвекции и диффузии применён так называемый связанный решатель. То есть, основные уравнения решались в виде единого связанного набора. Уравнения дополнительных моделей решались отдельно от связанного набора последовательно, одно за другим.

Система уравнений Навье-Стокса замыкалась с помощью модели турбулентности Лангтри и Ментера, которая основана на дополнении модели турбулентности ББТ к-ю двумя дифференциальными уравнениями: одним для перемежаемости и одним для критерия начала ламинарно-турбулентного перехода [10].

Для моделирования лучистого теплообмена применялась модель дискретных ординат в её классической формулировке [6], которая обеспечивает получение достоверных данных по уровню лучистых тепловых потоков в задачах с оптически тонкими средами. Эта модель представляется уравнением лучистого переноса.

Поверхность тела считалась полностью непрозрачной с диффузным отражением и коэффициентом черноты =1.0. На поверхности задавалась фиксированная температура 1260 К, то есть граничное условие 1-го рода.

Построение расчетной сетки

Для получения достоверных расчётных данных при наличии в потоке скачков уплотнения принципиально важно использовать расчетную сетку с высоким разрешением в области больших градиентов температуры и давления.

Задача решалась в двухмерной осесимметричной постановке. То есть, рассматривалось только одно меридиональное сечение тела. При этом считалось, что течение во всех меридиональных сечениях идентично.

Была построена расчётная сетка из четырехугольных ячеек, размер которой составил 5733 ячейки. В ходе измельчения в области больших градиентов температуры ее размер увеличился до 32652 тысяч ячеек.

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

Результаты расчета

Достоверность полученных при численных исследованиях результатов была проверена путем их сопоставления с результатами опубликованных ранее исследований. Рассматривались результаты исследований процессов газовой динамики и теплообмена при обтекании затупленных тел различных конфигураций в широком диапазоне изменения чисел Маха. Например, в работе [7] опубликованы данные экспериментальных исследований на ударных трубах и сверхзвуковых высокоэнтальпийных газодинамических стендах. В работах [8, 9] опубликованы результаты расчётов численными методами, с помощью решения дискретных аналогов системы Навье-Стокса. Наличие таких данных позволило провести верификацию представленной математической модели.

На рисунке 2 представлено распределение теплового потока по поверхности затупленного конуса в сравнении с данными работы [8].

1E+008 -q

q [w/m2] 1E+007

1E+006

1E+005

0 1 2 3 4 5 6 _S^Rn_

Wood

- Абсолютно каталитическая поверхность

- Некаталитическая поверхность

Рис. 2. Результаты расчёта распределения теплового потока

по поверхности тела.

Сопоставление результатов показало, что полученные данные хорошо согласуются с данными работы [8] в плане распределения теплового потока по поверхности тела. Следует отметить, что в работе [8] расчёт проводился только для случая некаталитической твердой стенки.

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

экзотермическими, они приводят к увеличению теплового потока в стенку.

Следует отметить, что кривая распределения теплового потока имеет классическую форму кривой при ламинарном теплообмене. Это объясняется малым характерным линейным размером исследуемого тела. Число Рейнольдса, вычисленное по характерному линейному размеру (диаметру затупления), составляет около 4000.

На рисунке 3 представлены значения распределения безразмерного давления по поверхности затупленного конуса в сравнении с данными работы [8].

1Е+000 -п

р/(р»о :

1Е-001 -

1Е-002

0 1 2 3 4 5 6 _SARn_

Wood

- Абсолютно каталитическая поверхность

- Некаталитическая поверхность

Рис. 3. Результатов расчёта распределения безразмерного давления по

поверхности тела.

Сопоставление результатов показало, что полученные данные хорошо согласуются с данными работы [8] в плане распределения статического давления по поверхности тела.

На рисунке 4 представлено распределение массовой концентрации атомарного кислорода (O) по оси ординат, проходящей через точку с координатой X/Rn = 5 в сравнении с данными работы [8].

1 -

0,8 -

0,6 -Со -0,4 -

0,2 -

0 —|-1--1--1--г--г--1

0 0,05 0,1 0,15 0,2 0,25

с0

Wood

- Абсолютно каталитическая поверхность

- Некаталитическая поверхность

Рис. 4. Результаты расчёта распределения массовой концентрации атомарного кислорода O по оси ординат в точке X/Rn = 5. Координаты от 0 (поверхность тела) до 1,0 - пограничный и сжатый слой.

Видно, что результаты, полученные для некаталитической поверхности, также хорошо согласуются данными работы [8].

Выводы

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

Показано, что полученные результаты по распределению теплового потока и распределению статической температуры по поверхности тела хорошо согласуются с результатами работы [8].

Получена хорошая корреляция результатов распределения массовых концентраций компонента O по нормали к поверхности тела в точке X/Rn=5 с результатами работы [8].

Таким образом, предложенная математическая модель может быть использована для решения газодинамических и теплотехнических задач при проектировании теплонапряженных элементов конструкции гиперзвуковых ЛА.

Работа выполнена при поддержке РФФИ, грант №11-08-00828, № 13-08-01328 а, 14-08-00982.

Библиографический список

1. Никитин П.В. Тепловая защита: Учебник. - М.: Изд-во МАИ, 2006. - 512 с.

2. Molchanov A. M., Bykov L. V., Donskikh V. V. A Numerical Method for Solving the Equations of Supersonic Chemically Reacting Turbulent Flows // Тепловые процессы в технике. 2012. Т.4, № 11. С.496-500.

3. Быков Л.В., Молчанов А.М., Янышев Д.С. Численный метод расчета сверхзвуковых турбулентных течений с химическими реакциями // Вестник Московского авиационного института. 2010. №3, Т.17. С.108-119.

4. McBride B. J., Gordon S., Reno M. A. Coefficients for Calculating Thermodynamic and Transport Properties of Individual Species // National Aeronautics and Space Administration. Office of Management Scientific and Technical Information Program. 1993.

5. Modest M. F. The Weighted-Sum-of-Gray-Gases Model for Arbitrary Solution Methods in Radiative Transfer // J. Heat Transfer. 1991. №113. pp. 650-656.

6. Chui E. H., Raithby G. D. Computation of Radiant Heat Transfer on a NonOrthogonal Mesh Using the Finite-Volume Method // Numerical Heat Transfer. 1993. Part B. №23. pp. 269-288.

7. Hollis B.R., Perkins J.N. Hypervelocity heat-transfer measurements in an expansion tube // AIAA Paper. 1996. № 96-2240. pp. 12.

8. Wood W. A., Eberhardt S. Dual-Code Solution Strategy for Chemically-Reacting Hypersonic Flows // AIAA Paper. 1995. № 95-0158.

9. Widhopf G. F., Wang J. C. T. A TVD Finite-Volume Technique for Nonequilibrium Chemically Reacting Flows // AIAA Paper. 1988. № 88-2711.

10. Menter F. R., Langtry R. B., Likki S. R., Suzen Y. B., Huang P. G., and VolkerS."A Correlation-Based Transition Model Using Local Variables: Part I — Model Formulation". (ASME-GT2004-53452), 2004.