CC BY
164
МОДЕЛИРОВАНИЕ ПРОЦЕССОВ
УДК 533.6
А. Л. Железняков а, С. Т. Суржиков
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПОЛЯ ТЕЧЕНИЯ ПРИ ВХОДЕ В АТМОСФЕРУ ЗЕМЛИ СПУСКАЕМОГО АППАРАТА С АЭРОДИНАМИЧЕСКИМ КАЧЕСТВОМ
Приведены результаты численного моделирования обтекания сверхзвуковым потоком идеального газа спускаемого космического аппарата "Fire II" при входе в атмосферу Земли. Моделирование проводилось на основе численного решения системы уравнений Навье-Стокса на неструктурированных сетках в трехмерной постановке с учетом сжимаемости газа.
Методами вычислительной аэродинамики исследованы особенности структуры полей течения и теплофизические процессы во всей возмущенной области вокруг спускаемого аппарата от головной ударной волны до дальнего следа, поведение местных аэродинамических и тепловых характеристик, а также структура сложных трехмерных вихревых течений у поверхности аппарата в условиях невозмущенного набегающего потока, соответствующих различным точкам возможных траекторий входа спускаемого аппарата в атмосферу Земли. Вычисления проведены в диапазоне изменения углов атаки а от 0 до 20°.
В рамках ряда моделей турбулентности исследовано влияние угла атаки а на поле течения, вихревые приповерхностные течения и местные аэродинамические и тепловые характеристики вдоль передней (носовой) и задней поверхностей. Результаты расчетов сопоставлены с данными летных испытаний, лабораторных экспериментов и приведенными в литературе данными расчетно-теоретического анализа.
Ключевые слова: численное моделирование, аэротермодинамика, газовая динамика, теплообмен, сверхзвуковые течения, спускаемые космические аппараты, уравнения Навье-Стокса, турбулентность, неструктурированные сетки.
При движении в атмосфере спускаемые космические аппараты (КА) подвергаются большим тепловым и динамическим нагрузкам, и при проектировании перспективных спукаемых аппаратов большое значение имеют расчетно-теоретические методы анализа полей течения и тепловых нагрузок, адекватность которых реальным условиям может быть установлена при всестороннем анализе возможных физических моделей и расчетных схем и сопоставлении их с данными натурных испытаний. Наиболее полная информация по характеристикам поля течения и тепловым нагрузкам получена в рамках программы полетов КА "Fire II" при входе в атмосферу Земли со скоростью,
превышающей вторую космическую. В ходе эксперимента проведены также измерения радиационного и суммарного тепловых потоков к передней и задней поверхностям.
В расчетном анализе условия невозмущенного набегающего потока соответствовали последнему участку траектории входа КА "Fire II" в атмосферу [1]. Вычисление тепловых потоков сильно усложнялось тем, что на заключительном этапе спуска течение около КА может переходить в турбулентное, наиболее интенсивное в области следа. Причиной значительного повышения теплового потока к поверхности аппарата может стать изменение режима течения. Отметим, что вычисление тепловых потоков для турбулентного режима течения более сложная задача, чем для ламинарного. Основная проблема заключается в выборе турбулентной модели, наиболее подходящей для исследуемых процессов обтекания. Обычно, при проектировании тепловой защиты неточность, связанная с неопределенностью в выборе турбулентной модели, компенсируется большим коэффициентом запаса. Одна из целей настоящей работы — это анализ различных турбулентных моделей для потока совершенного газа при небольших числах Рейнольдса и больших числах Маха, порядка M = 16, что соответствует условиям спуска. Для замыкания исходной системы уравнений применялись двухпараметрические к—ш [2, 3], BSL (baseline) к—ш [4, 5] и SST (shear stress transport) k—ш [5, 6] модели турбулентности, а также SSG (Sarkar, Speziale and Gatski) и LRR (Launder, Reece and Rodi) модели рейнольдсовых напряжений [7, 8] и проводились расчеты с использованием ламинарной модели. Данные летных испытаний [9-12], с которыми сравнивались результаты моделирования, включают в себя как термохимически неравновесные, так и равновесные режимы течения.
Наибольшим тепловым нагрузкам подвержена лобовая часть спускаемого аппарата; проведенные расчеты показывают, что для КА "Fire II" конвективный тепловой поток, приходящийся на заднюю поверхность, составляет всего 2... 3 % максимального теплового потока в окрестности передней критической точки. Несмотря на это, численное моделирование теплового состояния задней поверхности КА также необходимо при проектировании тепловой защиты.
Моделирование теплового состояния космического аппарата "Fire II" проводилось в ряде предыдущих исследований методами вычислительной аэродинамики [13, 14]. Авторы работы [15], используя расчетный код, основанный на уравнениях Навье-Стокса и учитывающий излучение, получили тепловые характеристики КА "Fire II" на участке траектории между 1634-й и 1651-й секундами спуска, что соответствует высотам от 77 до 37 км. В работе [16] проведено сравнение результатов расчетов по теплообмену для нескольких моделей турбулентности на задней и передней поверхностях летательного аппарата,
схожего по конфигурации с КА "Mars Pathfinder". Исследование КА такой же конфигурации проведено в работе [17], но с учетом перехода к турбулентному режиму течения и последующим использованием турбулентной модели.
В ряде работ изучаются переходное и турбулентное течения около КА "Fire II". В частности, в работе [18] был использован метод усреднения уравнений Навье-Стокса по Рейнольдсу (RANS, Reynolds-averaged Navier-Stokes) в совокупности с k—ш и Spalart-Allmaras турбулентными моделями для моделирования поля течения вокруг КА "FireII". В последующих работах [19-21] исследовались возможности и ограничения метода моделирования отсоединенных вихрей (DES, detached eddy simulation) при расчетах теплового состояния спускаемых КА. В работе [22] приведены результаты RANS моделирования по тепловым нагрузкам с учетом неравновесной химической кинетики, в результате которых получено хорошее соответствие с данными летных испытаний.
Угол входа КА "Fire II" в атмосферу составлял ~3°, скорость входа — 11327 м/с и превышала вторую космическую. В процессе движения в атмосфере угол атаки а возрастал, принимая значения 7,7°, 13° и достигая максимума 19,5° на заключительном этапе эксперимента. Увеличение угла атаки является причиной изменения физических характеристик КА "Fire II" в процессе спуска. Поэтому необходимы моделирование полей течения около КА и расчет местных характеристик при угле атаки, отличном от нуля. Результаты расчетов полей течения при а = 0 показывают, что течение в следе за аппаратом крайне асимметрично.
В большинстве предшествующих работ по вычислительной аэродинамике спускаемых КА основное внимание уделялось исследованию тепловых и динамических нагрузок на заднюю поверхность аппарата и аэротермодинамики следа за ним; кроме того, вычисления проводились в основном для а = 0. В настоящей статье приведены описания геометрической модели спускаемого аппарата и условий летного эксперимента, применяемых вычислительных моделей и методик построения расчетных сеток, представлены результаты трехмерного RANS моделирования поля течения для КА "Fire II" при различных углах атаки с использованием разных моделей турбулентности и проанализировано их влияние на структуру сложных трехмерных приповерхностных течений, местные тепловые и аэродинамические характеристики как задней, так и передней поверхностей КА "Fire II".
Геометрическая модель КА и условия набегающего потока. Цель программы полетов КА "Fire II" [10] состояла в исследовании нагрева летательного аппарата, имеющего характерную форму для спускаемых КА "Аполлон", входившего в земную атмосферу на скорости,
Рис. 1. Спускаемый КА "Fire II" с трехслойной тепловой защитой передней поверхности и держателем:
1... 3 — бериллиевые тепловые щиты; 4 — передняя поверхность аппарата; 5 — радиометрические окна; 6 — фенольно-асбестовые прослойки; 7, 15 — теплоизоляция; 8 — задняя поверхность (конус с углом 66°); 9 — радиометрическое окно пропускания волн оптического диапазона; 10, 14 — телеметрические антенны; 11 — антенна С -диапазона; 12 — механизм отделения спускаемого модуля; 13 — держатель; 16 — фенольно-асбестовая защита
превышающей вторую космическую скорость. Изучение поля течения около задней поверхности этого КА проводилось в условиях реального спуска в атмосфере Земли, что предпочтительнее проведения экспериментов в аэродинамической трубе (например, из-за отсутствия влияния хвостовой державки).
Космический аппарат "Fire II" представляет собой полностью автоматизированный летательный аппарат, состоящий из двигательной установки и спускаемого модуля. Компоновочная схема спускаемого модуля [10] показана на рис. 1. Модуль состоит из затупленной передней 4 и конической задней 8 поверхностей. Передняя часть модуля имеет многослойную структуру и включает три бериллиевых тепловых экрана 1, 2, 3, чередующихся с тремя фенольно-асбестовыми теплозащитными прослойками 6. Первые два слоя тепловой защиты отстреливаются в процессе спуска. Третий тепловой щит не отстреливается и определяет форму летательного аппарата на последнем участке траектории спуска, который исследуется в настоящей работе. Задняя поверхность представляет собой усеченный конус с углом 66°.
Задняя поверхность аппарата, выполненная в виде конуса (рис. 2, а), оснащена двенадцатью калориметрами, расположенными по окружности вдоль трех радиальных направлений, соответствующих углам
а
VR 0.609
Рис. 2. Приборное оснащение спускаемого КА "Fire II" (а) и внешний контур геометрической модели модуля, используемой для вычислений (б) (все размеры даны в сантиметрах)
р = 0; 120 и 240°. Значения температуры, зафиксированные датчиками на высоте 35 км, обрабатываются для получения теплового потока на поверхности и используются для сравнения с результатами, полученными методами вычислительной аэродинамики. Изменение угла атаки КА "Fire II" во время спуска определяли с помощью гироскопа и акселерометра, находящихся на борту спускаемого аппарата.
В работах [23-25] представлены результаты расчета изменения угла атаки во времени и основные показатели перемещения модуля, а также проанализированы колебания, возникающие при движении. Для определения углов атаки использовались данные, полученные с помощью гироскопа и акселерометра, а также результаты статических и аэродинамических измерений, выполненных в аэродинамической трубе. Углы атаки для некоторых периодов определялись методами компьютерного моделирования. Аэродинамические силы и моменты вы-
числялись на основе данных, полученных с помощью акселерометра и гироскопа.
Линия, очерчивающая внешний контур КА "Fire II", моделируемый в настоящей работе, приведена на рис. 2, б. Геометрическая модель модуля (см. рис. 2, б) имеет некоторые упрощения, аналогичные упрощениям, принятым в работе [1]; в частности, антенна 11 (см. рис. 1) в расчетах заменена плоской поверхностью.
В настоящей работе моделировались условия, задающие невозмущенный набегающий поток в точке конечного участка реальной траектории входа КА в атмосферу, которая соответствует одной из наименьших высот (35 км) и для которой были доступны измерения.
Условия набегающего потока, используемые при моделировании течения
около КА "Fire II"
Высота, км..........................................................................................................35
Число Маха........................................................................................................16
Температура набегающего потока, К........................................................237
Плотность набегающего потока, кг/м3 ......................................................0,0082
Скорость потока, км/с......................................................................................5
Температура поверхности модуля, К........................................................553,3
Угол атаки, град..................................................................................................0; 5; 10; 15; 20
Число Рейнольдса (х106)..............................................................................1,76
Вычислительные модели. Исходная система уравнений. Поле течения вокруг спускаемого КА моделировали, используя уравнения Навье-Стокса в трехмерной постановке, уравнение неразрывности и уравнение сохранения энергии.
Уравнения Навье-Стокса и уравнение неразрывности:
2 д i ди\ д дх \ дх) ду
dv ди дх ду
д_
дг
ды ди дх дг
; (1)
£ + div (pvV ) = -1 - 31 0" divV) +
^ д ( д^ д ду \ ду) дх
v и ^ дх ду J
д_
дг
ды дv ду дг
(2)
d.v (pwv )=- др - з I (, divv)+
2 д / ды\ д
дг у дг I дх
i ды ди
х г
д_ ду
ды дгиN ду дг/
(3)
дР + а1у (ру) = о. (4)
Уравнение сохранения энергии
дТ др
Рср-д£ + рс'рУ §гай Т = ^ (Л §гай Т) + + У gradр + /Ф, (5)
где диссипативная функция определяется как
Ф 2 I ди\2 2 (ду\2 2 (ди,\2 ( дь ди\2 дх) \ду ) \дг ) \дх ду)
дь дь\2 ( ди дь\2 2 (ди дь дь
+ \ду + дг) + \дг + дх) 3 \дх + ду + д^ ^
Для замыкания уравнений Навье-Стокса применялись двухпара-метрические модели турбулентности и модели рейнольдсовых напряжений.
В двухпараметрических моделях решались два раздельных уравнения переноса, определяющих турбулентную скорость и пространственный масштаб турбулентного движения (длину пути смешения). В таких моделях используется градиентная гипотеза Буссинеска, связывающая рейнольдсовые напряжения с усредненными градиентами скоростей и турбулентной вязкостью. В рамках настоящей гипотезы вводится предположение об изотропности турбулентной вязкости.
Турбулентная вязкость, в свою очередь, принимается пропорциональной величине пути перемешивания и масштабу турбулентной скорости. В рассматриваемых двухпараметрических моделях масштаб турбулентной скорости вычисляется через турбулентную кинетическую энергию к, которая находится из решения соответствующего уравнения переноса. Пространственный масштаб турбулентности определяется с использованием двух параметров потока, обычно это турбулентная кинетическая энергия к и скорость турбулентной диссипации е. Скорость диссипации турбулентной кинетической энергии определяется из уравнения переноса для е.
к—ш модель [2, 3]. Полуэмпирическая модель турбулентности выполняется на основе уравнений переноса турбулентной кинетической энергии к и специфической скорости диссипации (частоты турбулентных пульсаций) ш.
Одно из преимуществ к—ш модели — это применимость данной формулировки к моделированию пристеночных течений при больших продольных перепадах давления для низкорейнольдсовых вычислений. Модель не включает в себя сложную нелинейную демпфирующую функцию, учитывающую влияние вязкости на характеристики турбулентности в окрестности стенки, необходимую для к—е модели,
поэтому к—ш модель имеет менее высокую точность решения, но при этом метод не такой трудоемкий.
Модель предполагает, что турбулентная вязкость связана с турбулентной кинетической энергией и частотой турбулентных пульсаций (псевдозавихренностью) соотношением
к (7) / = р-. (7)
ш
Модель предполагает решение двух уравнений переноса: одно — для турбулентной кинетической энергии, другое — для частоты турбулентных пульсаций. Тензор турбулентных напряжений вычисляется из концепции изотропной турбулентной вязкости.
Уравнение для кинетической энергии турбулентности:
д (рк)
+ у • (pUk) = V
dt
Уравнение для частоты пульсаций:
р + р J Vk
Ok
+ Pk - ß'ркш. (8)
^ + V • (pUw) = V •
р + р ) Vw о.
+ <x-Pk - ßpw2. (9)
В этих уравнениях Pk представляет собой производство турбулентной кинетической энергии:
Pk = ptVU • (VU + VUT) - 2V • U(3p*V • U + pk); (10) константы модели имеют следующие значения:
в = 0,09; а = 5/9; в = 0,075; ак = 2; = 2. (11)
BSL (Baseline) k—ш модель [4, 5]. Главный недостаток k—ш модели — это сильная чувствительность к условиям набегающего потока. Вследствие зависимости от параметра частоты турбулентных пульсаций ш, принимаемого на входе в расчетную область, может быть получен значительный разброс в результатах моделирования. Этот нежелательный эффект был устранен благодаря применению гибридного метода. Метод, заключающийся в использовании k-ш модели в пристеночных областях и k—e модели [26] в областях, находящихся на достаточном удалении от стенки, был применен Ментером [4, 5]. Он заключается в некотором преобразовании k—e модели в k—ш формулировку и в последующем добавлении соответствующих уравнений. Таким образом, k—e модель дополняется стыковочной функцией 1 — F1. Данная функция принимает значение Fi = 1 вблизи поверхности и обращается в ноль за пределами пограничного слоя, т.е. на линии границы пограничного слоя и за его пределами k—e модель возвращается к первоначальной стандартной формулировке.
k-ш модель:
д (pk)
öt д (рш) т
+ V • (pUk) = V
+ V • (риш) = V
р + ^ Vk Okl J
+ Pk - ß'ркш;
(12)
р + ) Vw
ш
+ ai kPk - ßipw2. (13)
Преобразованная к-е модель:
д (рк)
т
+ V • (pUk) = V
р + ^ Vk
Ok,2)
+ Pk - ß'pku;
(14)
д (рш) m
+ V • (риш) = V •
р + ) Vw Ош2
+
+ 2р—1— УкУш + а2тРк - в2рш2. (15) Ош2 ш к
Теперь, если уравнения видоизмененной к—е модели дополнить функцией 1 — и добавить соответствующие уравнения к—ш модели, то получим ББЬ модель турбулентности:
д (рк)
т
д (рш)
т
+ V • (pUk) = V
+ V • (риш) = V •
р + Vk
Ok3 J
+ Pk - ß'pkw;
(16)
р + — )Уш +
+ (1 - ^\)2р—1— УкУш + аз ШРк - взрш2. (17) Ош2 ш к
Коэффициенты новой модели — это линейная комбинация соответствующих коэффициентов моделей, лежащих в основе метода:
Фз = ЗД + (1 - *\)Ф2. (18)
Следовательно, константы модели определяются следующим образом:
в' = 0,09; а1 = 5/9; в1 = 0,075; ак1 = 2; аш1 = 2;
a2 = 0,44; ß2 = 0,0828; Ok2 = 1; = 1/0,856.
(19)
Полученная модель сочетает в себе преимущества k—e (нечувствительность к параметрам свободного потока) и k—ш моделей (дает хорошие результаты в пограничном слое в непосредственной близости к стенке, особенно при наличии большого неблагоприятного градиента давления).
Shear Stress Transport (SST) k—ш модель [5, 6]. При применении гибридных моделей часть проблем удается избежать, однако, характеристики таких моделей существенно ухудшаются при наличии отрывных
течений. При использовании гибридной BSL модели турбулентности плохо предсказывается отрыв потока под действием неблагоприятного градиента давления. Модель дает завышенные значения турбулентной вязкости: k—£ и k—ш модели не учитывают перенос касательных напряжений.
Для устранения перечисленных недостатков была разработана модификация модели Ментера, так называемая SST (shear stress transport — перенос касательных напряжений) модель. SST модель переноса касательных напряжений является столь же экономичной, как и k-£ модель, но обеспечивает более высокое качество воспроизведения процессов для широкого ряда потоков и граничных условий. Она дает очень точное предсказание отрыва потока при неблагоприятных градиентах давления.
В этой модели точное поведение параметров переноса получено с помощью ограничения в формулировке турбулентной вязкости
Vt =--,-WWT, (20)
max(a1w, sF2)
где
Vt = pt/p; (21)
F2 — стыковочная функция, аналогичная функции F1 предыдущей модели.
Основу модели составляют те же два уравнения переноса, что и в случае BSL модели: уравнение переноса для частоты турбулентных пульсаций видоизмененной k-£ модели, дополненное стыковочной функцией 1—F1 (17), и уравнение переноса турбулентной кинетической энергии модели Уилкокса (16). Коэффициенты SST модели так же, как и в предыдущем случае, получаются как линейная комбинация коэффициентов k—ш и k—£ моделей (18), (19).
Существенное различие заключается в определении стыковочных функций (причина успеха метода), которые в данной модели представляют собой сложные функции безразмерной координаты стенки у и параметров турбулентного потока:
F1 = tanh(argl); (22)
( ( Vk 500v\ 4pk \
argi = min max ——, , —-- ; (23)
y ye шу у2ш J CDkwаШ2y2 J
CDkw = max ( 2p—1— VkVu, 1 • 10-10 ) ; (24)
V J
F2 = tanh(arg2); (25)
( 2Vk 500V .
arg2 = max —-, . (26)
\ в'шу у2ш
Модели рейнольдсовых напряжений. Двухпараметрические модели турбулентности, описанные ранее, предполагают справедливость градиентной гипотезы Буссинеска, линейно связывающей рейнольдсо-вые напряжения с усредненными градиентами скоростей и турбулентной вязкостью. В этом случае главные направления тензора скоростей деформации соответствуют главным направлениям тензора напряжений. Это разумно для весьма простого деформированного состояния, особенно, когда константы модели тщательно откалиброваны для рассматриваемых классов течений, но может быть не применимо при моделировании сложных полей деформаций, появляющихся при наличии закрутки потока или сложной геометрии. При таких условиях необходимо использовать более сложную связь между напряжением и скоростью деформации.
Модели переноса рейнольдсовых напряжений (Reynolds Stress Models) определяют турбулентные напряжения непосредственно по уравнениям переноса для каждого компонента тензора напряжений, что требует решения дополнительных шести взаимозависимых уравнений совместно с уравнением скорости турбулентной диссипации e, по которому определяется пространственный масштаб турбулентности.
Дифференциальное уравнение переноса рейнольдсовых напряжений имеет вид
дри & и ät
+ V • (pU & и & и) =
2k
= P + ф + V • р + -Csp— Vu & и - -8ре, (27)
3
2
3
где Р = -р(и ® и(Уи)т + (Уи)и ® и).
Уравнение переноса рейнольдсовых напряжений в тензорной записи можно представить как
дтц + = Pij + фгз + д
m
дхк
дхк
дт
V-
дхк
+ сг^ ) - е
-ij,
(28)
где Tij = pUiUj — тензор рейнольдсовых напряжений; eij = -5ijpe —
определение диссипации с использованием гипотезы Колмогорова о
2 к2 (дЩщ\
локальной изотропии; cij = -csp— —- — турбулентный перенос
3 e V dxk J
[Launder, Reece and Rodi].
Итак, получаем систему уравнений для модели рейнольдсовых напряжений:
Jt (pUUj) + дх (Uk рии) =
k д = Pij + фц +
дхъ
2 k2 \ ди и
р + ö Csp—
3 е / дх
H^j
- 2 Sijре (29)
2
k
— дифференциальное уравнение переноса для тензора рейнольдсовых напряжений;
д(рс) д с д
+ —(ßUkс) = T(c£lP - ое2ре) +
, Pt р +--
О г
дс dxk
(30)
дЬ Охи к Охи
— уравнение для скорости диссипации турбулентной энергии;
фЦ = Фи! + ф'Ц2 (31)
— корреляционное соотношение для давления и деформации, где
Фт = —рс I Csla + C,
-',s2
aa — i a • aö t t ,
(32)
Фт2 = —Cr\Pa + Cr2pkS — Cr3pkS^a • a+
2
—с 3
2 3(
+ Cripk(aST + SaT — 2a • S^j + Crbpk (aWT + WaT) ; (33)
U ^ U 1 /п77 / Г,,ПТ\
a = —;--- о — тензор анизотропии; S = - ( VU + (VU) —
k
скорость деформации; W = 1 ^VU — (VU— завихренность.
Выше приведена общая форма записи для корреляции. В зависимости от значений соответствующих коэффициентов корреляция может быть линейной и квадратичной. По значениям коэффициентов корреляции различают три стандартных модели рейнольдсовых напряжений.
В моделях LRR-IP и LRR-QI ("IP" — Isotropisation of Production, "QI" — Quasi-Isotropic) [8] корреляция линейная (табл. 1).
SSG модель была разработана в работе [7]; при моделировании использовалась квадратичная корреляция (см. табл. 1).
Таблица 1
Значения констант для LRR и SSG моделей турбулентности
Модель Cs C£1 Cs2 Cs1 Cs2 Cr1 Cr2 Cr3 Cr4 Cr5
LRR-IP 0,22 1,45 1,9 1,8 0 0 0,8 0 0,6 0,6
LRR-QI 0,22 1,45 1,9 1,8 0 0 0,8 0 0,873 0,655
SSG 0,22 1,45 1,83 1,7 -1,05 0,9 0,8 0,65 0,625 0,2
Однако, как показывает практика, учет анизотропии касательных напряжений, т.е. использование полного тензора в моделях высокого порядка, для многих случаев не дает никаких преимуществ по сравнению с двухпараметрическими моделями, хотя требует значительно больших вычислительных ресурсов. К тому же модели рейнольдсовых
напряжений имеют худшую сходимость из-за сильной связи уравнений.
Далее в табл. 2 приведен краткий обзор используемых моделей с точки зрения вышеперечисленных основных преимуществ и недостатков.
Таблица 2
Краткий обзор примененных моделей турбулентности
Модель Преимущества Недостатки
к—ш 1. Применимость данной формулировки к моделированию пристеночных течений при больших продольных перепадах давления 2. Модель не включает в себя сложную нелинейную демпфирующую функцию, учитывающую влияние вязкости на характеристики турбулентности в окрестности стенки, необходимую для к—e модели, следовательно, метод не такой трудоемкий 1. Сильная чувствительность к условиям набегающего потока (в зависимости от параметра ш, принимаемого на входе в расчетную область, может быть получен значительный разброс в результатах моделирования) 2. Низкая точность решения в области свободного течения
BSL 1. Нечувствительность к параметрам свободного потока (в частности, к ш) 2. Сочетает преимущества к—ш и к— e моделей применительно к пристеночным течениям и области свободного течения соответственно 1. Модель плохо предсказывает отрыв потока под действием неблагоприятного градиента давления. 2. Модель дает завышенные значения турбулентной вязкости
SST 1. Обладает всеми преимуществами BSL модели 2. Модель учитывает перенос касательных напряжений 3. Дает очень точное предсказание отрыва потока при неблагоприятных градиентах давления Не изучены
LRR 1. В рамках моделей решается дифференциальное уравнение переноса для каждого компонента тензора рейнольдсовых напряжений 2. Данные модели дают самые лучшие результаты в случае сложных течений, при наличии закрутки потока или сложной геометрии 1. Как показывает практика, учет анизотропии касательных напряжений в моделях высокого порядка часто не дает никаких преимуществ по сравнению с двухпараметрическими моделями 2. Требуется значительно больше вычислительных ресурсов 3. Модели рейнольдсовых напряжений имеют худшую сходимость из-за сильной связи уравнений
SSG
Построение расчетных сеток. При решении задач вычислительной гидродинамики процесс построения сетки обычно наиболее трудоемкий, требующий больших временных затрат. Точность и эффек-
тивность численного исследования существенно зависят от используемой расчетной сетки.
В настоящей работе моделирование трехмерного поля течения около КА "Fire II" проводилось с использованием неструктурированных тетраэдральных сеток.
На рис. 3 приведены различные фрагменты неструктурированной расчетной сетки, применявшейся для моделирования внешнего обтекания модели спускаемого КА "Fire II". Из рисунка хорошо видна степень сгущения узлов сетки к обтекаемой поверхности. Размер ячеек, непосредственно прилегающих к поверхностям КА "Fire II" составляет 10-4 м, что является достаточным для определения местных аэродинамических и тепловых характеристик передней и задней поверхностей КА. Размерность сетки 1188329 тетраэдральных элементов.
Неструктурированные сетки имеют бесспорное преимущество в том, что:
1)для большинства задач, представляющих практический интерес, форма расчетной области очень сложна и построение структурированной расчетной сетки в этом случае становится чрезвычайно трудоемким даже при использовании многооболочного подхода. Таким образом, неструктурированные сетки могут быть приспособлены к любой
Рис. 3. Фрагменты трехмерной расчетной сетки (справа показана структура сетки вблизи обтекаемого модуля и на его передней и задней поверхностях)
желаемой геометрии (сетка имеет большую гибкость применительно к сложным геометриям);
2) процесс построения неструктурированных сеток может быть автоматизирован в высокой степени, что может существенно сократить временные затраты на генерацию сеток. Построение качественной неструктурированной сетки занимает максимум несколько часов. Для создания структурированной сетки иногда требуется несколько недель;
3) неструктурированные сетки проще генерировать.
Неструктурированный метод имеет и ряд недостатков:
1) точность данного метода невысока, так как тетраэдр не имеет параллельных граней. Для получения требуемой точности результатов необходимо гораздо больше расчетных элементов, чем в случае применения структурированной сетки;
2) сетка должна быть в высокой степени изотропна. Соседние элементы могут лишь очень незначительно отличаться по форме и размерам. Это также приводит к существенному увеличению числа необходимых расчетных элементов;
3) неструктурированная сетка зависит от качества самой геометрической модели; даже незначительная ошибка в задании геометрии может сделать построение таких сеток невозможным;
4) некоторые трудности возникают при интерполяции расчетных данных на неструктурированных сетках для наглядного представления результатов расчетов в требуемых плоскостях и на поверхностях;
5) решатель, использующий неструктурированные сетки, требует больше памяти и расчетного времени. Численные алгоритмы, основанные на неструктурированной топологии сетки, более трудоемкие и сложные. Это связано с увеличением числа ячеек и граней по сравнению с шестигранными сетками. Тетраэдральная сетка из N узлов имеет около 6N ячеек и 12N граней, в то время как шестигранная состоит из порядка N ячеек и 3N граней;
6) для записи и хранения неструктурированных сеток требуется больше информации (дополнительная память требуется для хранения информации о связях ячеек сетки).
Производительность современных компьютеров и развитие параллельных вычислений позволяют решать все перечисленные проблемы, возникающие при использовании неструктурированного метода.
Для генерации сетки использовался метод распространяющегося фронта (advancing front method).
Алгоритм включает следующие основные этапы:
1) построение двумерной триангуляции тем же методом на границах расчетной области и поверхности КА;
2) сетка на границах области определяется как начальный фронт;
3) для каждой элементарной площадки фронта определяется оптимальное положение нового узла (из условия непересечения и с учетом расположения всех уже существующих соседних узлов);
4) новый узел и элементарная площадка фронта определяют новый тетраэдрический элемент. По мере образования ячеек сетки новые элементарные площадки добавляются к фронту ("использованные" для построения тетраэдров площадки фронту уже не принадлежат). Таким образом, активный фронт распространяется внутрь расчетной области. Тетраэдры строятся от граней, образующих фронт, и наращиваются внутрь расчетной области.
Результаты расчетов. На рис. 4, 5 (см. 4-ю полосы обложки) приведены результаты численного моделирования в виде полей чисел Маха при различных углах атаки. В каждом случае данные представлены в центральном осевом сечении, в плоскости наклона. Наиболее яркие особенности структуры течения отмечены на рис. 4 (результаты соответствуют нулевому углу атаки а).
Структура течения включает в себя головную ударную волну, отстоящую от передней части летательного аппарата; обширную зону формирования возвратных течений, прилегающую к задней поверхности КА "Fire II" сразу за точкой отрыва пограничного слоя; зону рециркуляции за КА. Внешний невязкий поток над рециркуляционной зоной разворачивается в направлении от оси, образуя ударную волну повторного сжатия, отходящую от горловины потока. Поток, набегающий на лобовую поверхность спускаемого КА, испытывает сильное сжатие и последующее расширение при стекании с кромки переднего теплового щита, образуя область расширения потока.
Сама форма головной ударной волны существенно изменяется при вариации угла атаки (см. рис. 5, 4-я полоса обложки). Ударная волна симметрична относительно оси при нулевом угле атаки. При возрастании угла атаки головная ударная волна все больше прилегает к наветренной стороне (результат сильного сжатия потока) и отходит от подветренной стороны передней части КА. Застойная область вблизи носовой части капсулы также постепенно смещается вниз с ростом угла атаки.
Главный эффект от вариации угла атаки выражается в изменении формы и размера рециркуляционной области за капсулой и изменении структуры течения (см. рис. 5, 4-я полоса обложки). При а = 0 отрыв пограничного слоя происходит сразу после затекания потока на усеченный конус задней поверхности КА, при этом образуется осе-симметричная рециркуляционная зона, прилегающая к задней поверхности. Если угол набегания потока не равен нулю, то зона возвратных течений за КА "Fire II" наклонена и меньше в размерах. При а = 0° набегающий поток плавно обтекает наветренную сторону усеченного
конуса задней поверхности без образования отрывных течений, хотя на подветренной стороне отрывная зона сохраняется и не изменяет расположения. При а = 5° след за обтекаемым телом меньше, чем в случае нулевого угла атаки, главным образом из-за позднего отрыва потока на наветренной поверхности. При дальнейшем возрастании угла атаки эффект, связанный с более поздним отрывом потока, более чем компенсируется благодаря возрастанию наклона пограничного слоя (по отношению к усеченному конусу аппарата). Поэтому при значительных углах падения (а > 10°) ширина следа возрастает по сравнению с нулевым углом атаки. Похожие структуры течений для различных углов атаки были получены в работе [1]. Отличие заключается в том, что по результатам настоящей работы отрыв течения на подветренной поверхности конуса КА "Fire II" происходит несколько позже. Также некоторые расхождения наблюдаются в размерах рециркуляционной зоны.
Значительные изменения в структуре потока при вариации угла атаки также можно проследить по рис. 6 (см. 3-ю полосу обложки). Здесь для нескольких углов атаки структура течения представлена в виде линий тока, нанесенных на поле чисел Маха.
Как видно из рис. 6, при нулевом угле атаки ближний след состоит из осесимметричного тороидального вихря, а также видны две его отдельные петли в центральном осевом сечении. При а = 5° петля вихря наблюдается только на подветренной стороне и ее размер постепенно уменьшается при возрастании угла атаки. Вихревая структура формируется вплоть до значений углов атаки порядка а = 15°. Центр вихревого образования с ростом угла атаки перемещается вверх по потоку вдоль подветренной конусной поверхности.
На рис. 7 с помощью линий тока (линий тензора касательных напряжений) показана структура трехмерных вихревых течений у задней поверхности КА "Fire II" при разных углах атаки.
Расходящиеся по поверхности (не пересекающиеся линии тока) соответствуют безотрывному обтеканию задней поверхности КА. Схождение линий в одной точке соответствует отрыву потока.
При а = 0 вихревая структура осесимметрична и покрывает всю заднюю поверхность летательного аппарата и плоское основание усеченного конуса. Даже при маленьком угле атаки (а = 5°) размеры вихревой зоны значительно уменьшаются, при этом поток имеет безотрывной характер течения уже только на наветренной стороне. Симметрия относительно центрального осевого сечения (плоскости наклона) сохраняется при любых углах атаки.
Местные тепловые характеристики. Далее исследуется влияние различных моделей турбулентности на местные тепловые характеристики (плотность конвективных тепловых потоков) вдоль передней (носовой) и задней поверхностей (рис. 8, 9).
Иллюстрации к статье А.Л. Железняковой, С.Т. Суржикова «Численное моделирование поля течения при входе в атмосферу Земли спускаемого аппарата с аэродинамическим качеством»
Рис. 4. Структура поля течения при нулевом угле атаки
Рис. 5. Распределение числа Маха в поле течения КА «Fire II» в плоскости наклона при разных углах атаки
Иллюстрации к статье А.Л. Железняковой, С.Т. Суржикова «Численное моделирование поля течения при входе в атмосферу Земли спускаемого аппарата с аэродинамическим качеством»
Рис. 6. Структура течений около КА «Fire II» в плоскости наклона при разных углах атаки (линии тока нанесены на поле чисел Маха)
V
J
Рис. 7. Линии тензора касательных напряжений на задней поверхности модуля при разных углах атаки
Рис. 8. Распределения плотности теплового потока вдоль всей поверхности КА "Fire II" от передней критической точки для различных моделей турбулентности
Как и ожидалось, при нулевом угле атаки распределение местных характеристик осесимметрично. Плотность теплового потока достигает максимальных значений в области передней критической точки.
Плотность теплового потока на передней и задней поверхностях летательного аппарата изменяется от 5 Вт/см2 в зоне отрыва потока (соответствует s/D = 0,9, s — расстояние вдоль поверхности аппарата от передней критической точки; D — максимальный диаметр модуля) до 250Вт/см2 (максимальное значение зависит от используемой модели турбулентности) в области передней критической точки (см. рис. 8).
35
30 25
а 20
о CQ
15
10
0
7 я - а = 0 8' - avg а • ~(р = 0 * ~(р = 120( а - <р -240° 1 - SST 2-к-о) 3 - BSL 4- SSG 5 - LRR 6- Lam 3 щ Ш щ
1 U 7. * Я
ч ■ ш 1 8.
/VI1 ? * ■ J д * / * ■ г revisits üVirfu»1 J
■ ■ ■ ■ 1 в ■ ■ glSilМШШШ9* 5 4
0,6 0,8 1,0 1,2 1,4 5/ В
Рис. 9. Распределения плотности теплового потока на задней поверхности КА для различных моделей турбулентности. На графическую область также нанесены данные летного эксперимента для калориметров, расположенных под углами ^ = 0; 120 и 240°, и результаты численного моделирования [1]
Из рис. 8, 9 видно, что значения местных тепловых характеристик на задней поверхности спускаемого КА, полученные с помощью различных моделей турбулентности, отличаются незначительно, в то время как распределения плотности тепловых потоков на передней поверхности очень различны.
Как отмечалось ранее, задняя поверхность спускаемого КА была оснащена двенадцатью калориметрами, расположенными по окружности вдоль трех радиальных направлений, соответствующих углам = 0; 120 и 240° (см. рис.2). Данные летного эксперимента [10], зафиксированные этими измерительными приборами, также помещены в графическую область распределения плотности теплового потока в зависимости от безразмерного параметра s/D (см. рис.9). Значения s/D = 0,64 и 1,43 соответствуют началу и концу задней конической поверхности КА "Fire II" соответственно.
Обратим внимание на значения плотности теплового потока, взятые из летного эксперимента [10] (см. рис.9). Показания датчиков получены в точке траектории спускаемого модуля на высоте 35 км. Угол атаки в данной точке траектории достигал максимального значения порядка 15°. Исходя из структуры течения при ненулевом угле атаки (отрывной характер течения на подветренной стороне и плавное обтекание на наветренной, см. рис. 5,6, 4-я и 3-я полосы обложки) можно было предположить, что показания датчиков на наветренной и подветренной сторонах капсулы должны значительно отличаться в зависимости от угла расположения. Например, датчик, расположенный в наветренной зоне плотного прилегания потока, должен регистрировать большие тепловые нагрузки, чем датчик, расположенный в зоне отрывного течения. Однако данные, полученные с помощью летного эксперимента (см. рис. 9), нанесенные на общий график черными треугольниками, кружками и ромбами, не показывают значительного отличия в зависимости от расположения по углу. Это может быть вызвано тем, что при входе в плотные слои атмосферы спускаемый модуль вращается с угловой скоростью порядка 160 об/мин. Возможно, с помощью датчиков были получены значения тепловых нагрузок, усредненные по углу. В связи с этим для сравнения численных результатов с экспериментальными данными [10] при ненулевом угле атаки были получены усредненные по углу значения плотности теплового потока для каждого рассматриваемого угла атаки. Если учесть, что в процессе спуска, кроме вращения, наблюдаются значительные колебания угла атаки, то для сравнения c экспериментальными данными [10] распределения плотности тепловых потоков были также усреднены по углам атаки. Получившиеся распределения для различных моделей турбулентности приведены на рис. 8, 9. Также на общую диаграмму
(см. рис. 9) нанесены усредненное по углам атаки от 0 до 20° распределения плотности теплового потока и распределение теплового потока для нулевого угла атаки [1].
При нулевом угле атаки в зоне отрыва потока (s/D = 0,9) наблюдаются минимальные значения плотности теплового потока. Результаты расчетов, приведенные в работе [1], указывают на минимальное значение тепловых потоков при s/D = 0,7 и, как следствие, на более ранний отрыв потока на подветренной стороне конуса КА "Fire II". При движении вдоль задней поверхности тепловой поток возрастает незначительно. Распределение, полученное в работе [1], при нулевом угле атаки дает более значительное возрастание плотности теплового потока при движении вдоль конуса задней поверхности КА "Fire II". Усреднение значений плотности теплового потока по углам атаки приводит к снижению этой характеристики приблизительно на 40% для задней поверхности модуля.
Из рис. 9 видно, что полученные результаты хорошо согласуются с экспериментальными данными.
Заключение. В работе приведены результаты численного моделирования обтекания спускаемого КА "Fire II" с использованием уравнений Навье-Стокса в трехмерной постановке при измении угла атаки от 0 до 20° и применении различных моделей турбулентности. Результаты расчетов показали, что структура течения значительно изменяется даже при маленьком угле атаки. Пограничный слой на наветренной стороне плотно прилегает к поверхности летательного аппарата. На подветренной стороне, напротив, при натекании на заднюю коническую поверхность КА "Fire II" происходит отрыв потока с образованием рециркуляционной зоны. Две противоположные тенденции наблюдаются для плотности теплового потока на задней конической поверхности капсулы. Плотность теплового потока возрастает в наветренной зоне плотного примыкания пограничного слоя с ростом угла атаки, в то время как в рециркуляционной зоне отрыва на подветренной стороне усеченного конуса с ростом угла атаки наблюдается снижение тепловых нагрузок. Это приводит к тому, что усредненные по углу (по всей поверхности аппарата) значения плотности теплового потока практически не изменяются при вариации угла атаки. Усреднение по углам атаки приводит к более равномерному распределению плотности теплового потока на задней поверхности КА "Fire II" и, следовательно, к несколько меньшим тепловым нагрузкам на задней поверхности КА по сравнению с результатами моделирования при а = 0.
СПИСОК ЛИТЕРАТУРЫ
1. Sinha K., Vadivelan C. Effect of Angle of Attack on Re-entry Capsule
Afterbody Flowfield // AIAA Paper. - 2008.
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
24
Wilcox D. C. Multiscale Model for Turbulent Flows // AIAA Journal. - 1988. -Vol. 26. - No. 11. - P. 1311-1320.
W i l c o x D. C. Reassessment of the Scale Determining Equation for Advanced Turbulence Models // AIAA Journal. - 1988. - Vol. 26. - No. 11. - P. 1299-1310. M e n t e r F. R. Multiscale model for turbulent flows // In 24th Fluid Dynamics Conference. American Institute of Aeronautics and Astronautics, 1993. M e n t e r F. R. Two-equation eddy-viscosity turbulence models for engineering applications // AIAA Journal. - 1994. - Vol. 32. - No. 8. - P. 269-289. M e n t e r F. R. Zonal Two Equation k—u Turbulence Models for Aerodynamic Flows // AIAA Paper 93-2906. 1993.
SpezialeC. G., SarkarS. and G a t s k i T. B. Modeling the pressure-strain correlation of turbulence: an invariant dynamical systems approach // J. Fluid Mech.
- 1991. - Vol. 227. - P. 245-272.
Launder B.E.,Reece G. J. and R o d i W. Progress in the development of a Reynolds stress turbulence closure // J. Fluid Mech. - 1975. - Vol. 68. - P. 537-566. Cornette E. S. Forebody Temperatures and Calorimeter Heating Rates Measured During Project Fire II Reentry at 11.35 Kilometers Per Second // NASA TM X-1305. Nov. 1966.
S l o c u m b T. H. Project Fire Flight II Afterbody Temperatures and Pressures at 11.35 Kilometers Per Second // NASA TM X-1319. - 1966.
CauchonD. L. Radiative Heating Results from the Fire II Flight Experiment at a Reentry Velocity of 11.4 Kilometers Per Second // NASA TM X-1402. - 1966. GnoffoP. A. A Code Calibration Program in Support of the Aeroassist Flight Experiment // AIAA Paper. - No. 89-1673. - June 1989.
G r e e n d y k e R. B. and H a r t u n g L. C. Convective and Radiative Heat Transfer Analysis for the Fire II Forebody // Journal of Spacecraft and Rockets.
- 1994. - Vol.31. - No.6. - P.986-992.
G o k c e n T. Calculation of Nonequilibrium Radiating Shock Layers // Journal of Thermophysics and Heat Transfer. - 1995. - Vol. 9. - No. 1. - P. 34-40. OlynickD. R. HenlineW. D. Chambers L. H. and Candler G. V. Comparisons of Coupled Radiative Navier-Stokes Flow Solutions with the Project Fire II Flight Data // AIAA Paper. - No. 94-1955. - June 1994. B r o w n J. L. Turbulence model validation for Hypersonic flows // AIAA Paper. -2002-3308. - 2002.
N a n c e R. P., H o 11 i s B. R. H o r v a t h T. J. A 11 e r S. J. & H a s s a n H. A. Computational study of hypersonic transitional wake flow // Journal of thermophysics and heat transfer. - 13(2). - 1999.
S i n h a K. RANS simulation of Fire II re-entry configuration' // International Symposium on Shock Waves. - 2005.
SinhaK., Barnhardt M., and C a n d l e r G. V. Detached Eddy Simulation of Hypersonic Base Flows with Application to Fire II Experiments // AIAA Paper. -2004-2633.-2004.
S i n h a K. Grid Sensitivity of Detached Eddy Simulation of a Mach 16 Re-entry Configuration // AIAA Paper. - 2007-1115. - 2007.
SinhaK. Effect of Reynolds Number on Detached Eddy Simulation of Hypersonic Base Flow // AIAA Paper. - 2007-1457. - 2007.
SinhaK. and R e d d y S. K. Hypersonic Turbulent Reacting Flow Simulation of Fire II Re-entry Vehicle // AIAA Paper. - 2007-0805. - 2007. Wright M. J., P r a b h u D. K. and Martinez E. R. Analysis of Apollo Command Module Afterbody Heating Part I: AS202 // Journal of thermophysics and heat transfer. - 20(1). - 2006.
E d q u i s t K. T. Afterbody Heating Predictions of a Mars Science Laboratory Entry Vehicle // AIAA Paper. - 2005-4817. - 2005.
25. S c a 11 i o n W. I. and L e w i s J. H. Body Motions and Angles of Attack During Project Fire Flight II Re-entry // NASA TM D-4183. - 1967.
26. W i 1 c o x D. C. Turbulence modeling for CFD. - 1994.
Статья поступила в редакцию 23.12.2008
Александра Львовна Железнякова окончила МГТУ им. Н.Э. Баумана в 2008 г. Сотрудник лаборатории "Радиационная газовая динамика" Института проблем механики РАН. Автор 6 научных работ в области теплофизики и вычислительной газовой динамики.
A.L. Zheleznyakova graduated from the Mosckow State Technical University n.a. N.E. Bauman in 2008. Researcher of Laboratiry for Radiative Gas Dinamics of the Institute for Problems in Mechanics RAS. Author of 6 publications in theory of heat and mass transfer and computional fluid dynamics.
Сергей Тимофеевич Суржиков окончил МВТУ им. Н.Э.Баумана в 1975г., МГУ им. М.В. Ломоносова в 1978 г. Чл.-кор. РАН, д-р физ.-мат наук, заведующий лабораторией "Радиационная газодинамика" Института проблем механики РАН, заведующий кафедрой "Физическая и химическая механика" МФТИ, профессор кафедры "Теплофизика" МГТУ им. Н.Э. Баумана. Автор более 350 научных работ в области теплофизики и радиационной газодинамики.
S.T. Surzhikov graduated from the Bauman Moscow Higher Technical School in 1975 and Moscow State University n.a. M.V. Lomonosov in 1978. Dr. Sc(Phis.), Head of the Radiative Gas Dynamics Laboratory of the Institute for Problems in Mechanics Russian Academy of Sciences. Head of the Chair of Physical and Chemical Mechanics, MIPT. Author of more than 350 pub1ications in radiative gas dynamics and theory of heat and mass transfer.
В издательстве МГТУ им. Н.Э. Баумана в 2008 г. вышла в свет книга
Колесников К.С.
Рассказ о моей жизни. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2008. - 360с.
Автобиографическая книга профессора МГТУ им. Н.Э. Баумана, академика РАН Константина Сергеевича Колесникова представляет собой яркое жизнеописание человека интереснейшей судьбы. Перед нами история личности на фоне крупнейших событий двадцатого столетия, пример целеустремленности фронтовика-бауманца, который жаждал учиться и добился максимальной самореализации.
Читатель — студент или выпускник МГТУ им. Н.Э. Баумана — почерпнет из этой книги немало ценной информации о развитии университета во второй половине XX в., воспитании молодежи, замечательных ученых, блестящих педагогах, которыми по праву гордится наша alma mater.
Неподдельная искренность автора, рассказывающего о пройденном им пути, побуждает к серьезному размышлению, поиску ответов на волнующие современника вызовы нынешней эпохи.
По вопросам приобретения обращаться по тел. (499)-263-60-45; e-mail: press@bmstu.ru
