Серия «Математика»
Том 2 (2009), №1, С. 63-71
Онлайн-доступ к журналу: http://isu.ru/izvestia
УДК 517.958 / 517.977 Оптимизация форм тел, аблирующих в атмосфере планеты *
М. А. Аргучинцева
Иркутский государственный университет
Аннотация. В статье исследована задача о нахождении оптимальной начальной формы аблирующего осесимметричного тела, движущегося по баллистической траектории и обладающего минимальным суммарным радиационным нагревом поверхности.
Ключевые слова: оптимальное аэродинамическое проектирование, интегро-диф-ференциальная система.
Оптимальное аэродинамическое проектирование космических аппаратов является одной из интересных и перспективных областей приложения методов оптимального управления и вариационного исчисления. Известно, что при движении с гиперзвуковыми скоростями в атмосферах планет поверхности тел подвергаются интенсивному радиационному и конвективному нагреву. В результате этого происходит разрушение теплозащитного покрытия, сопровождающееся уносом массы тела. Изучение уноса массы космического аппарата представляет большой интерес в связи с созданием зондов, способных входить в атмосферы внешних планет Солнечной системы.
Одним из эффективных путей снижения тепловых потоков к телу и, следовательно, уменьшения уноса массы покрытия является соответствующий выбор формы спускаемого аппарата, которая должна удовлетворять целому ряду аэродинамических и конструктивных требований. В данной работе исследуется вопрос нахождения оптимальной начальной формы тела, для которой суммарный радиационный нагрев поверхности вдоль баллистической траектории движения будет минимальным.
* Работа выполнена при финансовой поддержке РФФИ, грант 08-01-00709.
1. Введение
Нужно отметить, что унос массы теплозащитного покрытия приводит к изменению первоначальной формы тела, его аэродинамических коэффициентов и, в конечном счете, влияет на траекторию движения тела. Поэтому важно совместно решать задачи обтекания и нагрева тела с учетом абляции поверхности и движения тела по траектории. В наиболее общей постановке проблема поиска оптимальной начальной формы тела, обгорающего в атмосфере планеты, представляет собой сложную задачу оптимального управления системой нестационарных уравнений Навье-Стокса с управляемыми границами. Исследование этой задачи сопряжено со значительными трудностями.
В настоящее время в научной литературе рассматриваются различные упрощения исходной постановки задачи. Так, в работах [1,2] исследовалось уравнение абляции, когда унос массы происходил под воздействием конвективного нагрева. Соответствующие решения для стационарных форм аблирующих тел получены в работах [3-5]. Исследованию уноса массы осесимметричных и пространственных тел за счет радиационного теплового потока в атмосферах Земли и Юпитера при неизменных условиях в набегающем потоке посвящены работы [6-9]. В [8] рассмотрена постановка сопряженной задачи, которая включает уравнения движения тела переменной формы и уравнения радиационной газовой динамики для химически равновесных течений невязкого газа в ударном слое. В работе [9] была решена задача об оптимальной стационарной форме аблирующего тела с точки зрения минимума уноса массы при неизменных условиях в набегающем потоке. При этом предполагалось, что радиационный нагрев зависит от локального угла наклона поверхности.
Данная работа посвящена более общей постановке задачи о нахождении оптимальной начальной формы аблирующего осесимметричного тела, движущегося по баллистической траектории и обладающего минимальным суммарным радиационным нагревом поверхности. Общая нестационарная задача обтекания разбивается на две взаимосвязанные части:
1) нестационарную задачу о движении аблирующего тела по баллистической траектории в атмосфере планеты. Данная задача описывается системой обыкновенных дифференциальных уравнений движения, зависящей от аэродинамических коэффициентов в каждой точке траектории, и интегро-дифференциальным уравнением абляции тела за счет теплового нагрева;
2) стационарную задачу обтекания тела в заданной точке траектории с учетом теплообмена. Асимптотическое решение данной задачи позволило получить функционалы для радиационного теплового потока [10] и волнового сопротивления [11], явно зависящие от формы тела и параметров обтекания.
Интегрирование уравнений движения с учетом найденных коэффициентов сопротивления и теплообмена позволяет поставить и эффективно решить соответствующие вариационные задачи нахождения оптимальных начальных форм аблирующих тел.
2. Постановка задачи
Рассмотрим гиперзвуковое движение осесимметричного тела в атмосфере планеты вдоль баллистической траектории [8]
M(t)dv(t) = -2p(t)v2(t)CB(t)S(t)+ M(t)g sin Y(t), Ml =cos Y(t)(-^-V<tl)
dt v(t) Re ^
dH() = -v(t)sin y (t). (1)
Здесь t - время (t € [0,T]); v(t), CB(t), M(t), R(t), S(t) (S(t) =
nR2(t))- соответственно скорость, коэффициент волнового сопротивления, масса, радиус и площадь донного сечения тела; p(t) - плотность газа на высоте H(t); Y(t) - угол наклона траектории движения; g - ускорение свободного падения; Re - средний радиус планеты. Для изотермической атмосферы имеем
p(t) = ро exp(-AH (t)),
где ро - плотность атмосферы на уровне поверхности планеты; A-1 -шкала высот для плотности.
При больших сверхзвуковых скоростях входа в плотные слои атмосферы (порядка второй космической скорости и выше) преобладающую роль в нагреве тела играет радиационный тепловой поток от сжатого ударного слоя. Тогда уравнение абляции, описывающее обгар осесимметричного тела за счет радиационного нагрева, имеет вид [6]:
= qn(x,t) (2)
He,pT
где х - координата вдоль оси тела; у(х, Ь) - функция формы тела; Чк(х, Ь) - величина локального радиационного теплового потока в заданную точку поверхности тела; Hef - эффективная энтальпия разрушения теплозащитного покрытия; рт - плотность тела.
Уравнения (1), (2) решаются при начальных условиях
Н (0) = Но, у(0) = Уо, 7 (0) = 7о; (3)
у(х, 0) = Уо(х), уо(0) = 0, уо(Ь) = Ко, (4)
где Но, Уо, 7о - заданные высота, скорость и угол входа в атмосферу;
Уо(х), Ь, Ко - начальная форма тела, ее длина и радиус донного сечения.
Выражения для коэффициента волнового сопротивления и локального радиационного теплового потока осесимметричного тела получены в работах [10,11] из асимптотического исследования задач стационарного гиперзвукового обтекания тел в заданной точке траектории с учетом теплообмена
4 ЬГ У(х,і) () Ах
св«> = Км / 1 + (дм^у • (5)
хо(і) + \ дх )
( л г т р(() ”3(і) 1 X ;у(М) (^)2К+П
т(х’1) = Г“^^^ -Ш) } Л , (д,(„) ^+9)/2 Х
хо(і)
1 +
х I—7дГ(—т2 (да^)2М+8
(ду^щ ^ I д^ )
(і + (4й )2)'
1+ И' !1 + {-^Г) ё‘1 Л /„„,„) )2Л(2«+7)/2
N+5 ' N+4
(6)
Здесь N - параметр излучения газа [10], хо(і) - координата критической точки тела
і
хо(і) = I Но(Т) Ат.
] Не/Рт
Яе/рт
0
Для вычисления коэффициента радиационного теплообмена в критической точке затупленного тела С^(Ь) можно использовать известные эмпирические зависимости [10] для различных составов атмосфер планет. Например, для воздуха при 9 км/с < -и(£) < 20 км/с; 10-7г/см3 < р(£) < 10-4г/см3; 0.3 м < Е(Ь) < 40 м имеем
Г
СЯ0(^ =
2 + 8.9Г5/6 ’
Г = 0.00344(0.1г>(і))6(р(і) 106)о-35Яа5(і), Ь = 2(п + 4)Г.
Нужно отметить, что в классе тонких осесимметричных тел, когда (Мх^) << 1, выражение для радиационного теплового потока (6)
можно упростить, тогда уравнение абляции примет вид
ду(х,і) ^(x,^) (7)
ді у(х,і)’
х
х
2М+11
р(х,і)= К(і) ! у(8,і){^
хо(і)
Масса тела равна
ь
М (і) = рт п
/
у2(х, і) Ах.
(8)
хо(і)
Поставим следующую оптимизационную задачу: в классе гладких функций, удовлетворяющих граничным условиям (4), найти начальную форму аблирующего тела уо(х) = у(х, 0), которая подчиняется системе уравнений (1)—(3), дополнительному ограничению на массу аппарата
(М* - заданное минимально возможное значение массы аппарата) и минимизирует функционал суммарного радиационного нагрева
Таким образом, поставленная задача является задачей оптимального управления начальными условиями системы, состоящей из нелинейного интегро-дифференциального уравнения абляции (2) и системы обыкновенных дифференциальных уравнений движения тела по баллистической траектории (1). Для ее решения использовался численный метод локальных вариаций [12]. Опишем общую схему метода.
Для дальнейших рассуждений функцию формы тела у(х, і) удобно переписать в безразмерном виде п(С,і) (С Є [0,1], і Є [0,Т]):
Тогда искомая начальная форма тела По (С) = уо(х)/К будет удовлетворять граничным условиям по(0) = 0, По(1) = 1. На плоскости (С,і) ведем равномерную сетку узлов
М (і) > М„
(9)
о хо(і)
3. Метод решения
п = у/к, С = х/Ь, т = К/Ь.
іі = і Аі, Аі = Т/п, і = 0,1,...,п; Сі = І АС, АС = 1/т, і = 0,1,...,т.
Первоначальное число узлов по временной переменной (п = 80) затем уточнялось в процедуре автоматического выбора шага сетки АЬ метода Рунге-Кутта интегрирования системы уравнений движения (1). Соотношение между количеством узлов т и п выбиралось из условия выполнения соответствующего критерия Куранта для конечно-разностной схемы уравнения абляции (2).
В качестве начального приближения искомого решения По(С) задавалась одна из следующих зависимостей:
п0 = {1/4, По = «1/2, п0 = С3/4, По = СКак показали расчеты, все эти начальные приближения в процессе варьирования сходились к одной и той же оптимизирующей кривой, однако, в случае конического начального приближения п0 = С скорость сходимости метода была значительно быстрее. Для выбранного начального приближения по квадратурным формулам Симпсона производился подсчет интегралов С в (5), М (8), стоящих в правых частях уравнений движения (1). Далее, на каждом временном слое проводилось совместное интегрирование указанной системы интегро-дифференциальных уравнений (1)-(2). Система обыкновенных дифференциальных уравнений движения (1) исследовалась методом Рунге-Кутта 4-го порядка.
Для решения уравнения абляции в случае тонких тел вращения (7) использовалась неявная разностная схема, которая на каждом временном слое приводила к решению нелинейного уравнения относительно
п(6- ,и)
'П(6' ,Ьг) - П(6' ,Ьг-1) '
П(& ,Ьг) + ^(Лз-1,и)+
+К(Ьг)п(Сз ,Ьг)
АЬ
п(С3 Л) - п(С- ,и)+11
т2И+10АС = 0.
АС
Для более общего случая нетонких осесимметричных тел использование неявного подхода было сопряжено с рядом трудностей, поэтому для уравнения абляции (2),(6) применялась явная разностная схема
П(& ,Ьг) - П(& ,Ьг-1) =
АЬ
= _ ди(С] ,Ьг-1) /. . 2 ГУ(С3 ,и-1) - П(С^-1,Ьг-1) I2
РТ Е(Ьг-1Н I- АС . •
На рис. 1 представлены примеры расчетов начальной п° = С3/4 и конечных форм аблирующих тел при Уо = {11; 15; 19} км/с, Но = 80 км, 7о = п/3, ро = 1.225 г/м3, Л-1 = 0.149 км, КЕ = 6371 км, т = 0.5,
рТ = 2200 кг/м3, Иef = 2 * 107 Дж/кг.
В конце процедуры вычислялось значение суммарного радиационного нагрева поверхности вдоль траектории полета Qя (10). Затем задавался шаг варьирования Н по ординате ц из условия [12]: Д£ ^ 0, Н/(Д£)2 ^ 0, необходимого для сходимости метода локальных вариаций. Процесс варьирования ординат По(С^') заключался в следующем: значения По(С^') последовательно заменялись на (п°(0) ± Н). Для новых точек находились приближенные значения функционалов сопротивления и массы (5), (8), проверялось условие (9) на массу аппара-та,интегрировалась система уравнений (1)—(2), и находилось значение целевого функционала Qя (10). Если значение функционала Qя уменьшалось, то в таблицу решений записывалось соответствующее значение ординаты (п°(&) ± Н). Таким образом, проводился полный перебор точек (^, По(^)±Н), ] =0,1,... ,т, и находилось первое приближение экстремали (По), целевой функционал которого не больше вычисленного на предыдущем приближении.
В процессе счета итерации с фиксированными Н и Д^ проводились до тех пор, пока в результате очередной итерации не изменялось ни одно из значений (пк(^)}. Затем шаг Н уменьшался вдвое, и процесс варьирования продолжался. Уменьшение Н при фиксированном Д£ проводилось до тех пор пока не удовлетворялось неравенство Н < 10-1°. Затем рассматривалось решение задачи на более мелкой сетке, т.е. чис-
ло разбиений т удваивалось. Функция ц = п(£) во вновь образованных точках находилось интерполяцией по соседним узлам. Процесс варьирования продолжался до тех пор, пока не выполнялось неравенство Д£ < 10-5.
Рис. 2.
4. Анализ результатов расчетов
Была проведена серия расчетов оптимальных начальных форм осесим-
метричных тел в широком спектре условий входа в атмосферу Земли и заданных геометрических характеристиках тела. Например, на рис.2 представлены оптимальные начальные формы аблирующих тел при следующих параметрах входа в атмосферу Земли: Уо = 15 км/с, Но = 80 км, 7о = п/3, р0 = 1.225 г/м3, Л-1 = 0.149 км, Ее = 6371 км. Соответствующие заданные геометрические и теплофизические характеристики тел: т = {1. — 0.5; 2. — 0.8}, рт = 2200 кг/м3, Не/ = 2 * 107 Дж/кг. В качестве минимально возможной массы аппарата М* в ограничении (9) бралась масса конуса с теми же геометрическими характеристиками.
Сравнение полученных начальных форм тел с традиционными телами (степенными телами и острым конусом) показало, что использование
оптимальных форм тел позволяет снизить суммарный радиационный нагрев поверхности до 50 % и, соответственно, уменьшить унос массы теплозащитного покрытия тела. Полученные результаты подтверждают необходимость совместного решения задач оптимизации, обтекания и теплообмена при изучении движения тел, аблирующих вдоль траектории полета в атмосфере.
Список литературы
1. Лунев В. В. Некоторые свойства уравнения абляции / В. В. Лунев // Изв. АН СССР. МЖГ. - 1977. - № 3. - С. 96-102.
2. Знаменский В. В. Численное решение уравнения уноса / В. В. Знаменский // Изв. АН СССР. МЖГ. - 1978. - № 2. - С. 147-154.
3. Воронкин В. Г. О стационарной форме тел при их разрушении за счет аэродинамического нагрева / В. Г. Воронкин, В. В. Лунев, А. Н. Никулин // Изв. АН СССР. МЖГ. - 1978. - № 2. - С. 138-146.
4. Мурзинов И. Н. Аналитическое решение о стационарных формах тел в условиях абляции / И. Н. Мурзинов // Газовая и волновая динамика. — М.: Изд-во Моск. ун-та, 1979. — С. 64-74.
5. Коняев В. Г. Аналитическое исследование изменения формы аблирующих тел при их движении в атмосфере со сверхзвуковыми скоростями /
B. Г. Коняев // Уч. записки ЦАГИ. - 1974. - № 6. - С. 125-133.
6. Левин В. А. Форма тонких тел при уносе вещества под воздействием лучистого теплового потока из ударного слоя / В. А. Левин, В. П. Марков, Н. Н. Пилюгин // Современные газодинамические и физико-химические модели гиперзвуковой аэродинамики и теплообмена. Часть 2. — М.: Изд-во Моск. ун-та, 1995. - С. 100-106.
7. Апштейн Э. З. Стационарная форма тел при их разрушении под действием тепловых потоков, зависящих от локального угла наклона поверхности / Э. З. Апштейн, Н. Н. Пилюгин // Изв. АН СССР. МЖГ. — 1981. — № 6. —
C. 137-143.
8. Апштейн Э. З. Унос массы и изменение формы трехмерного тела при движении по траектории в атмосфере Земли / Э. З. Апштейн, Н. Н. Пилюгин, Г. А. Тирский // Космические исследования. — 1979. — № 2. — С. 246-255.
9. Аргучинцева М. А. Экстремальные задачи радиационной газовой динамики / М. А. Аргучинцева, Н.Н. Пилюгин. — М.: Изд-во Моск. ун-та, 1997. — 197 с.
10. Пилюгин Н. Н. Динамика ионизированного излучающего газа / Н. Н. Пилюгин, Г. А. Тирский. — М.: Изд-во Моск. ун-та, 1989. — 312 с.
11. Теория оптимальных аэродинамических форм. — М.: Мир, 1969. — 507 с.
12. Черноусько Ф. Л. Вариационные задачи механики и управления / Ф. Л. Черноусько, Н. В. Баничук. — М.: Наука, 1973. — 238 с.
M. A. Arguchintseva
Shape optimization of bodies subliming in a planet atmosphere
Abstract. In this paper it is considered the optimal shape problem for subliming bodies moving along a ballistic trajectory with minimal radiation surface heat.
Keywords: optimal aerodynamic design, integro-differential system.
Аргучинцева Маргарита Александровна, кандидат физико-математических наук, доцент, Институт математики,экономики и информатики, Иркутский государственный университет, 664003, Иркутск, ул. К. Маркса, 1, тел.: (3952) 24-22-12, (marguch@math.isu.ru).
Arguchintseva Margarita, Irkutsk State University, 1, K. Marks St., Irkutsk, 664003, Phone: (3952) 24-22-12, (marguch@math.isu.ru).