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

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

www.mai.ru/science/trudy/

УДК 629.7.063 (075.8)

Расчёт тепло - и массообмена на поверхности спускаемого

космического аппарата

Никитин П.В. 1*, Павлюк Е.А. 2**

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

111401, Россия *е-тай: petrunecha @ gmail. ги * *е-таИ: рау1уикеа @poweгmcb. ги

Аннотация.

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

По указанному алгоритму выполнен расчёт траекторных параметров и

теплообмена на поверхности спускаемого космического аппарата на гиперзвуковом

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

Ключевые слова: спускаемый космический аппарат, тепло-массообмен,

гиперзвуковая скорость, тепловая защита, траекторные параметры, энтальпия

торможения, массовый секундный расход газа-охладителя

Используемые сокращения:

КЛА - космический летательный аппарат;

ПКТ - передняя критическая точка

Введение.

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

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

Достижение указанной цели потребовало решения нескольких задач:

- составить алгоритм и рассчитать траекторные параметры в соответствии с заданным законом изменения высоты и дальности полёта КЛА от времени.

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

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

1. Расчёт траекторных параметров КЛА при управляемом (штатном) и баллистическом (аварийном) спуске в атмосфере Земли.

Траектория полета КЛА характеризуется законом изменения высоты и дальности полёта во времени.

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

Для расчёта теплообмена была выбрана достаточно сложная, произвольная траектория спуска КЛА в виде изменения высоты от времени полёта

представленная на рисунке 1. На начальном участке а-Ь траектории аппарат совершает свободный баллистический полёт. Далее в точке Ь с использованием двигательной установки аппарат совершает манёвр, достигая максимальной высоты в точке с. В дальнейшем, из точки с аппарат совершает планирующий полёт до завершающей стадии.

Н,км 100

80 60 40 20

а

с

0 0 200 400 600 800 1000 1200 т,сек

Рис. 1. Траектория спуска КЛА.

КЛА представляет собой аппарат конической формы с притупленной сферой, радиус которой =1м. Миделево сечение аппарата Бмид = 10,9м . Масса спускаемого аппарата т = 20т, коэффициент аэродинамического сопротивления аппарата Сх = 0.3

где ах -баллистический коэффициент КЛА;

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

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

V (т ) = Vвв • ехр

( _ °х § 0р0 ) е - вн ( т )

2 - 2вв п

(1)

Затем, в точке Ь КЛА переводится в режим управляемого полёта. При этом его скорость рассчитывается уже по другому алгоритму:

Расчет скорости полета КЛА был произведен графически. Для этого траектория была разбита на элементарные участки (Рисунок 2).

Рис. 2. Сетка расчёта траекторных параметров.

И по каждому участку был найдена производная по времени. Итогом процедуры являлась зависимость вертикальной составляющей скорости по высоте. Но ее было недостаточно, так как спускаемый летательный аппарат летит не вертикально вниз, следовательно, для достоверного знания распределения скорости по времени нужно знать и горизонтальную составляющую. Была выбрана произвольная траектория горизонтального перемещения КЛА при спуске, указанная на рисунке 3.

Рис. 3. Траектория горизонтального перемещения КЛА при спуске

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

У =

V У

+

V У

(2)

2

2

Рис. 4. Изменение скорости спуска КЛА от времени.

Как известно [1] , при гиперзвуковых скоростях полета, когда перед КЛА образуется отошедшая ударная волна, давление торможения в передней критической точке рассчитывается с помощью уравнения Гюгонио в виде:

Р о = Р н •

к +1

к+1 к-1

2

к -1

к-1

2 к

Мк-1

м н

2к к -1

• мн -1

к-1

(3)

1

1

Рис. 5. Изменение давления торможения на траектории

спуска КЛА.

При известной скорости полёта расчёт энтальпии торможения производился

по известной термодинамической зависимости:

1 о = 1н +

У

2

(4)

2

оо

где 1н - термодинамическая энтальпия атмосферы на высоте Н полёта, У^ -скорость полета КЛА.

Рис. 6. Изменение энтальпии торможения на траектории

спуска КЛА.

Энтальпия является функцией температуры и давления, то есть

1о=/(Т,Р)

(5)

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

температура торможения.

2. Алгоритм расчёта тепло - массообмена на поверхности КЛА.

Наиболее полный анализ теплообмена на поверхности притуплённого тела для случая "равновесного" и "замороженного" пограничных слоев, как в передней критической точке, так и ее окрестности, провели Фей и Ридделл [4]. Алгоритм

расчета конвективного теплообмена на поверхности притуплённого тела включает 3 этапа:

- расчет теплообмена в передней критической точке;

- расчет теплообмена в области передней критической точки;

- расчет теплообмена вне области передней критической точки.

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

2.1. Алгоритм расчёта тепло - массообмена в передней критической точке КЛА.

Передней критической точкой называется точка, которая появляется на поверхности тела при нулевом угле атаки потока, омывающего тело. Газ в этой точке полностью тормозится, а критическими параметрами являются температура и давление. При этом давление в этой точке равно давлению торможения, равенство же критической температуры температуре торможения возможно лишь при условии, что стенка теплоизолирована и критерий Прандтля Рг=1 [1].

Тепловой поток рассчитывается с использованием модифицированного уравнения конвективного теплообмена Ньютона в виде:

2 о = {а / Ср )о •(I о - К) Параметр теплообмена рассчитывается с

Ридделла в виде [1]:

{а/Ср)о = о,763 • Рг;0'6• (РоАоГ • (Р^)оД

(6)

использованием уравнения Фея-

У с1х

\ о,5

1 + (Ьвп -1 )-

Б

(7)

о

где ¡0 и цк -коэффициенты динамической вязкости, которые определяются из таблиц для многокомпонентной высокотемпературной смеси при давлении, Р0 и температуре Т0 и Тщ соответственно;

Ро

р0 и pw - значения плотности смеси, р0

R ■ T

0

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

du с \ъ(р'0-рн

— = — I—--- градиент скорости,

dx R0 V Р'о

где коэффициент С = ^ 1 — ; Я0 - радиус притупления сферы;

Ро , Р0 плотность и давление - за ударной волной, соответственно; Рн - давление на высоте полета;

1 • СРу

Prw - критерий Прандтля, РГ

r~w г w

w

я

w

Id - энтальпия диссоциации молекул, входящих в состав

многокомпонентной высокотемпературной смеси.

ID=Z?=iIj-(Ci0-Ciw), где

I-* - энтальпия образования атомов, Cio и С iw - массовые концентрации атомов при температуре торможения (индекс "0") и температуре поверхности (индекс "w");

Le - критерий Льюиса-Семенова для многокомпонентной высокотемпературной газовой смеси

Ьв =

Р • Б 12 • СР Л

Показатель п критерия Ье выбирается в зависимости от состояния пограничного слоя [1].

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

п

•С ■

1=1

Р1 - плотность смеси

Л- коэффициент теплопроводности смеси, который рассчитывается из

соотношения:

Л = о.5 •

I х, •Л +

/ л-1 Vх!

г=1

V г:

,=1 Лг

где мольные доля ,-го компонента X, = С{

А А,

А и а - мольные массы 1-ои компоненты и смеси, соответственно;

коэффициент бинарной диффузии [5]:

Б12 = 2.6 • 1о

- 7

3/2

А + М2 А1 • М2

р о2 а11

ро о 1,2 а 1,:

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

до сверх, а затем дозвукой, в условиях отсутствия химических реакций, для

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

q о = а ■ (T о - T w ) Где коэффициент теплоотдачи рассчитывается из соотношения:

а=0,77 ■

1+0,074

~w

T

V T0

-1

J

Í \1/3 PgMG

V Pw№w J

^wPwPo ■ CP ■ P rw

(8)

На рисунке 7 представлена зависимость рассчитанного теплового потока от времени, на которой видно, что в ПКТ тепловой поток имеет точки экстремума.

Рисунок 7. Изменение плотности теплового потока в ПКТ при спуске

КЛА.

Появление первой точки экстремума является результатом баллистического полета и доказывает, что при входе в плотные слои атмосферы многоразового КЛА, полет следует проектировать управляемым. Далее до времени 700 сек, плотность

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

2.2. Алгоритм расчёта тепло-массообмена в окрестности передней

критической точки.

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

Расчет теплового потока проводился по всей траектории спуска КЛА. Скоростной режим постепенно снижался от гиперзвукового до дозвукого. В связи с этим, алгоритм расчета делился на расчет при наличии химических реакций (при гиперзвуковых скоростях) и при их отсутствии при скоростях полета, меньших, чем гиперзвуковые.

Алгоритм расчета плотности теплового потока в окрестности ПКТ при отсутствии химических реакций.

Плотность теплового потока рассчитывалась из соотношения:

2 о = а •(Те - Т„ )

(9)

Где коэффициент теплоотдачи рассчитывается:

N4 • Л

а =- ,

х

эф

где для ламинарного режима:

(1о)

Ми = о.332 • Кеэ0ф5 • Рг1/3 К К,

(11)

К =

1 + о.16

Т 1+^

Л ' 2т ^1/3

V -о

V т +1 у

о.5

К =

Г л1/3 РА

V у

для турбулентного режима:

Ми = о.о296 • Яеэ0ф8 • Рго43 • КТ

К- =

Г- \

1 У

Т

V То у

о.4

1+• у.и2

2

о.11

К =

г л°-6

А

V у

(12)

Расчет плотности теплового потока в окрестности ПКТ при наличии

химических реакций проводился по следующему алгоритму.

Плотность теплового потока рассчитывается с использованием модифицированного уравнения конвективного теплообмена Ньютона в виде:

qa = (а/ Cp )0 • (1е -)

(13)

Параметр теплообмена при ламинарном режиме[1]:

(а / Cp )о = 0,763 • Pi;4

-0,6

Д

w

X,

эф

f Л0,4

РеДе

V Р;Д; J

1 + (Len -1 )-

D

(14)

как и ранее n=0,63 для замороженного и n=0,52 для равновесного

пограничных слоев.

Re эфw =

u (x) • X,

эф

В уравнении (10) характерный размер хэф рассчитывается с помощью

соотношения:

Хэф =

0

J р; u (x) • R 2(x) • dx

Pw • u (x) • R 2(x)

(15)

Поставив под интеграл соответствующие функции и(х) = ß• х и R(x) ~x в непосредственной близости от передней критической точки сферы, получим:

_ 1 :^эф _ 4 •x

где х - длина дуги окружности от передней критической точки до расчетной точки на поверхности сферы. Вдали от критической точки, где R(x) * x, определение хэф

осуществляется графическим решением выражения (15).

0

v

w

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

При реальных полетах космических летательных аппаратов с гиперзвуковой скоростью на поверхностях его элементов конструкции возникают течения, характеризуемые громадными значениями чисел Рейнольдса. Такие течения обладают особыми свойствами, известными как "турбулентность". Закономерности таких течений в большой степени отличаются от закономерностей ламинарных течений. Следует отметить, что методы расчета теплообмена в турбулентном, химически реагирующем пограничном слое в значительной степени еще не отработаны. Поэтому расчет теплообмена при турбулентном режиме течения, как правило, проводится с использованием различных эмпирических методов, допускающих серьезные упрощающие предположения. Например, в настоящее время главная концепция методологии расчета теплоотдачи в турбулентном химически реагирующем пограничном слое базируется на предположении о справедливости аналогичных методов, разработанных для химически нейтрального турбулентного пограничного слоя [1].

Расчет параметра теплообмена можно проводить с использованием соотношения:

(а / Ср )о = о, о296 • Яе оф8 • Рг;0 • 6 • А •

^ Х. л.

Гт Ло'4

эф

Т

V Те у

1+• у-ш2 2

о.11

где коэффициент восстановления г =

Эффективная длина хэф рассчитывается из следующего интегрального

соотношения:

x , _

x 5

4 ( x I •р

w

0 (17)

J R 4( x) P; ■ u (x) • dx

R4(x) P;•u(x)

Вдали от передней критической точки, где R(x) * x, значения хэф для каждой точки поверхности х, определяются графическим решением интеграла (17), где для сферической части КЛА R(x) рассчитывается из соотношения:

R( x) _ Ro ' VX •( 2 — X ) , а для конической части из:

_ r (Х -1)• smm+1 _ _ х

R( x ) Ro ^^ , где все Xi _

_______, где все » г»

СОБЮ К0

Расчет расстояний по горизонтальной оси от ПКТ, Х{, указан на рисунке 8.

Аз

Х<\

._Хв

Рис. 8. Процедура расчета эффективной длины, хэф, на поверхности притуплённого сферой конуса.

Расчетная зона от ПКТ до заданной расчетной точки разбивается на 5 участков. На сфере притупления проводятся три сечения через точки х!,х2,х3. Причем первое сечение проводится через точку, которая является ближайшей к ПКТ (х0). Второе сечение проходит через точку, лежащую на звуковой линии, третье сечение через точку, которая является точкой сопряжения сферы и конуса, далее на поверхности конуса проведены произвольно еще два сечения: одно через точку, х5, а второе между точками х3 и х5.

На рисунке 9 представлены результаты расчёта изменения траекторных параметров по траектории спуска КЛА.

Рис. 9. Изменение траекторных параметров на траектории

спуска КЛА.

На рисунках 10-14 видно, что экстремумы, выявленные в ПКТ, сглаживаются в области передней критической точки и на боковой поверхности КЛА.

с/, кВт/м 2 1000

800

600

400

200 0

'20° 0 200 400 600 800 1000 ¿200 Ъсек

Рис. 10. Зависимость теплового потока от времени в точке

с координатой х=0.00195.

с], кВт

8000\-

7000 —

6000-

5000-

4000-

3000-

2000 -

1 ООО —

о -

-1000 0 ш 400 б00 ш 1000 1Ш т сек

Рис. 11. Зависимость теплового потока от времени в точке

с координатой х=0.268.

q, кВт/м 1ООО

son 600 400

200 О

~:000 200 400 600 800 1000 1200 т^ек

Рис. 12. Зависимость теплового потока от времени в точке с координатой

х=0.658

21

ц. кИт/м2

16000 14000 • 12000 10000 НООО 6000 4000 2000 о

0

Рис. 13. Зависимость теплового потока от времени в точке с координатой

х=0.83.

Ч> кВт/м2 ¡4000 12000 10000 НООО 6000 4000 2000

О

О 200 400 600 800 1000 1200 %сек

Рис. 14. Зависимость теплового потока от времени в точке с координатой

х=1.057.

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

3. Выбор системы тепловой защиты. Расчет массового секундного расхода

Анализ полученных уровней расчетных параметров показал, что для данного класса КЛА в ПКТ можно использовать проницаемую стенку с пористым охлаждением. Наличие вдува газа в пограничный слой заметно снижает плотность теплового потока. Плотность теплового потока при вдуве газа описывается, как и на непроницаемой стенке, законом Ньютона, в виде соотношений:

При сопоставлении соотношений между плотностями тепловых потоков в проницаемую и непроницаемую стенки можно сделать вывод, что между ними должна существовать функциональная зависимость типа qw =/(до) [6,7].

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

qо - qw- qвд=0, где (20)

q0- плотность теплового потока при отсутствии вдува; qw - плотность теплового потока при наличии вдува; qвд - плотность теплового потока, блокированного вдувом.

охладителя.

для передней критической точки

(18)

(19)

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

двд = у • т „ • (Iе -1.„ ) ,

(21)

где У - коэффициент пропорциональности, называемый коэффициентом вдува газа-охладителя через проницаемую стенку. стенку. Согласно [1]:

у = 0,6

Г - Л0'24

Я

V Дохл J

г

х-0,03

ч>,охл

V I J

(22)

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

Преобразовав уравнение баланса (20) с учетом уравнений (6), (18) и (21) получится выражение:

СР ) „ = 1 - т ..

у--—__(22)

(а/Ср )0 (а/Ср )0

Отношение (а/ ^р) представляет собой безразмерную величину, которая называется безразмерным массовым расходом газа-охладителя, обозначается:

о= 171

" (а / Ср)

0

Тогда (1.21) принимает вид:

(а/ Ср)

w

(а/Ср )

= 1 -у-О

w

(23)

В свою очередь, используя (6), (18), выражение (23) можно преобразовать к

виду:

я

W _

я

= 1 -у-О

(24)

0

При уровне максимального теплового потока в стенку для достижения максимальной температуры поверхности, равной 1200 К через пористую стенку необходимо вдувать воздух с максимальным расходом 0,22 кг/м2 -с (Рисунок 15).

Рис. 15. Зависимость массового расхода газа-охладителя от времени.

На рисунке 16 представлена зависимость теплового потока с учетом вдува от

времени.

(], кВт/м2

400

350 300 250 200

150 100

50 О

-50

Рис. 16. Зависимость от времени плотности теплового потока в передней критической точке при наличии вдува.

\

0 4< Ю 6^0 ю ООО 12 00

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

Выводы:

1. Составлены алгоритмы и проведены расчеты траекторных параметров полета КЛА формы притупленного сферой конуса. Выявлены диапазоны изменения траекторных параметров, таких как температура торможения Т0, энтальпия торможения 10, давления торможения Р0', и скорости полета Уп0лета .

2. Составлен алгоритм и проведены расчёты тепло - и массообмена в непроницаемую стенку в передней критической точке, и отдельных точках в

окрестности точки торможения на высотах полета КЛА от 100км до 0км. Для указанных расчетных точек получены плотности тепловых потоков q0.

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

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

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты РФФИ 12-08-00145-а и 13-08-01328). Библиографический список.

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

2. Основы теории полета космических аппаратов / под ред. Г.С. Нариманова, М.К. Тихонравова. М.: Машиностроение, 1972. 608 с.

3. Варгафтик Н.Б., Справочник по теплофизическим свойствам газов и жидкостей. -М.: Наука, 1972. 720 с.

4. Авдуевский В.С., Галицейский Б.М., Глебов Г.А.. и др. Основы теплопередачи в авиационной и ракетно-космической технике / Под ред. В.С. Авдуевского, В.К. Кошкина - 2-е изд., перераб. и доп. - М.:Машиностроение,1992. 528 с.

5. Гиршфельдер Дж., Кертис Ч., Берд Р. «Молекулярная теория газов и жидкостей». М.: ИЛ, 1961, 929 с.

6. Анфимов H.A. «Тепло - и массообмен в окрестности критической точки при вдуве и отсосе различных разов через поверхность тела». - «Механика жидкости и газа». 1966, №1.

7. Анфимов H.A., Альтов В.В. «Теплообмен, трение и массообмен в ламинарном многокомпонентном пограничном слое при вдуве инородных газов». -«Теплофизика высоких температур», 1965,№3, С. 409-420.