CC BY
28
Серия «Математика»
2010. Т. 3, № 3. С. 67-79
Онлайн-доступ к журналу: http://isu.ru/izvestia
УДК 536.6.01
Численное моделирование турбулентных течений при сверхзвуковом обтекании тел *
Г. С. Журавлева
Иркутский государственный университет
Аннотация. Рассматривается гиперзвуковое обтекание затупленных осесимметричных тел потоком вязкого сжимаемого газа. Решение задачи получено конечноразностным методом для случая переменного вдува по обводу пористого затупления тела в широком диапазоне чисел Рейнольдса.
Ключевые слова: уравнения Навье-Стокса; вязкий ударный слой, турбулентное течение; условия Рэнкина-Гюгонио; ударная волна; численные методы.
Летательный аппарат при движении в атмосфере со сверхзвуковой скоростью испытывает интенсивное силовое и тепловое воздействие, которое может привести к разрушению его теплозащитного покрытия. Для снижения температуры поверхности используют вдув газа через пористое затупление тела. Наличие вдува газа приводит к изменению условий обтекания и, следовательно, к изменению теплового и силового воздействия атмосферы на движущийся летательный аппарат.
При сверхзвуковом обтекании пары тел, расположенных друг за другом, перед вторым телом при его достаточном удалении от первого реализуется стационарное течение с головным скачком уплотнения. Имеющиеся теоретические и экспериментальные данные свидетельствуют
о том, что в зависимости от расстояния между телами, возможны два типа обтекания второго тела: обтекание с возвратно-циркуляционными течениями в ударном слое и безотрывное обтекание. Численное исследование обтекания с образованием возвратно-циркуляционных течений требует привлечения полных уравнений Навье-Стокса. Во втором случае оказывается возможным применение более простых моделей вязкого ударного слоя. Следует заметить, что проведение лабораторных и натурных экспериментов с телами, находящимися на больших расстояниях друг от друга, представляет значительные трудности.
* Работа выполнена при частичной финансовой поддержке РФФИ, гранта 07-01-
Экспериментальному исследованию теплообмена при вдуве газа в пограничный слой на пластине и конусе посвящено большое количество работ. Однако для затупленных тел количество экспериментальных работ весьма ограничено. Поэтому важное значение приобретает численное исследование влияния вдува газа различной интенсивности вдоль образующей затупления на характеристики теплообмена и сопротивления тела. С увеличением скорости движения и давления в газе большую роль начинает играть аэродинамический нагрев и трение на поверхности тела при турбулентном течении в ударном слое.
В работе численно исследуется неравномерное сверхзвуковое обтекание затупленных осесимметричных тел. Для защиты поверхности тела от высоких тепловых потоков в пограничный слой через пористое затупление подается охладитель - другой газ (гелий или азот). Рассматривается процесс бинарной диффузии вдуваемого газа при условии, что в смеси не происходят химические реакции.
Система уравнений вязкого ударного слоя получается из осреднен-ных уравнений Навье-Стокса с помощью тех же предположений, что и в [2]. В безразмерных переменных в системе ортогональных координат х, у, связанных с поверхностью тела, уравнения имеют вид [4]
1. Постановка задачи
(1.1)
дх К \ду
V + д ( № дИЛ
) ду \а^К ду)
ду [К \ат, Зет,
Р = рТ,
гпШ TS Т) Т) Р°0 ^^О ^^0
ц = Т , K = eRe, Re =---------------------,
Ц0
e = ^, Y = ^, Цо = ц(То)= Т%,
2y cv
V2 Ц0ЦЕср с ЦоЦЕ
Т0 = -Z—, -Е = ---------7---, ScE = п ,
2cp Ae pDE
h = (ccp2 + (1 - c)cpi)T, cv = CCv2 + (1 - c)Cv1.
Здесь x определяет расстояние вдоль образующей тела, измеренное от критической точки; у - расстояние по нормали к поверхности тела, uVtx, evV0 - компоненты вектора скорости, соответствующие осям X, у; pooV2Р, е-1рор, Т0Т, h0h, cp, cv - соответственно давление, плотность, температура, энтальпия, удельные теплоемкости газовой смеси при постоянном давлении и объеме; Ц0ЦЕ, As, De - коэффициенты полных вязкости, теплопроводности и диффузии; с - массовая концентрация вдуваемого газа; к - продольная кривизна поверхности тела. Все линейные размеры отнесены к характерному линейному размеру R0, нормальная координата - к eR0, в качестве R0 выбирался радиус кривизны затупления тела при x = 0. Индекс w относится к величинам на поверхности тела, индекс £ обозначает суммарные коэффициенты, обусловленные молекулярным и турбулентным переносом. Индекс 1 относится к компоненте внешнего потока, 2 - к вдуваемой компоненте. Влиянием термодиффузии пренебрегаем, так как рассматриваются случаи умеренных градиентов температуры.
На ударной волне зададим обобщенные условия Рэнкина-Гюгонио, которые в сверхзвуковом приближении при неравномерном обтекании имеют вид [5]
p{v - U<ddx^j = PiViV о Р = Pi Vi2v ~ ’
PiViVoo (и - Viuo) = KдУу, (1.2)
т/ _ ЦЕ дс
piViV00C -- Q т,- О ,
ScsK ду
Т/ 2 т/2\ _ ЦЕ дТ 2ЦЕ ( ди
piViv°v +u -V0 = + -к{иду
Здесь иоо = cos a, v0 = - sin а, а - угол между касательной к
поверхности тела и осью симметрии.
На поверхности тела зададим условие прилипания для продольной составляющей скорости, расход газа, температуру стенки и условие для концентрации вдуваемого газа [7]
и = 0, рь = С(х), Т = Тт(х),
/1 \ де ,, „ч
рь (1 - е) = - зетк дУ (1.3)
Распределение газодинамических параметров в набегающем потоке
имеет вид [3]
У (г) = 1 — а ехр(—Ьг2), рг(т) = В [1 + С (1 — йУ?)]-\ (1.4)
pi(r) = const, ci(r) = 0,
B = 1 + C(1 - d), d = (1 - a)-2.
Здесь G (x) — pw vw/ pooVoo; VooV1, poop1, pooVoo p1, ci - соответственно скорость, плотность, давление и концентрация в набегающем потоке; R0r - расстояние до оси симметрии потока; a, b, C - параметры, характеризующие неравномерность набегающего потока; Voo, poo -скорость и плотность газа при r ^ то. Обтеканию тела равномерным потоком газа соответствует a = 0.
Модель тонкого вязкого ударного слоя (ТВУС), несмотря на широкое распространение, имеет меньшую точность по сравнению с другой известной моделью - уравнениями полного вязкого ударного слоя (ПВУС). Сравнение расчетов газодинамических параметров по обеим моделям для равномерного сверхзвукового ламинарного обтекания конусов приведено в [11], для неравномерного обтекания сферы - в [9]. Установлено, что максимальное отличие теплового потока для обеих моделей в случае неравномерного обтекания сферы при Re > 104 составляет около 10 %. Модель ПВУС дает более точные результаты на боковой поверхности удлиненных тел [11] при умеренных Re и в случае сильной неравномерности набегающего потока [9] при условиях, близких к отрыву.
Применимость модели тонкого вязкого ударного слоя с обобщенными условиями Рэнкина-Гюгонио на ударной волне для описания гиперзву-кового равномерного обтекания гладких затупленных тел обоснована в [10]. Указанная модель обобщена на случай неравномерного обтекания
в [3]. Система (1.1) правильно описывает течение в ударном слое при умеренном вдуве газа с поверхности тела [1].
В работе используется модель ТВУС для расчета параметров течения около гиперболоида при условиях, далеких от отрыва, и при Ев > 104, как более простая при численной реализации и достаточно точная в рассматриваемой области.
2. Модель турбулентности
Для замыкания системы уравнений (1.1)—(1.4) применяется алгебраическая модель [8] для полных коэффициентов переноса, согласно которой коэффициенты переноса равны
цт = Ц-(тг,Пк),
Ат = А-(п,щ)
1 + (- — 1) - (г—*
\(Гт ) - \1 — тг
(2.1)
Бт = Б-(п, Пк)
1+
Бе
Бс
Е
- V1 — тг
-(тг,Пк) =
(т2 — п1) + \/ (пк— т2)2 +
2т г
Цт ,2
тг = —, цт = р12 Ц
ди
ду
Здесь ц, цт - коэффициенты молекулярной и турбулентной вязкости; А, Ат - коэффициенты молекулярной и турбулентной теплопроводности; а = ц0цеР/А, ат = Ц0Цтер/Ат - ламинарное и турбулентное числа Прандтля; Бе = ц0ц/рБ, Бет = Ц0Цт/рБт - ламинарное и турбулентное числа Шмидта. Ламинарные, турбулентные числа Прандтля и Шмидта полагались в расчетах постоянными.
Молекулярная вязкость смеси ц в общем случае зависит от состава смеси и от температуры. Для каждого компонента ] использована приближенная зависимость ц- ~ Тш(и = 0, 5 ^ 1, 0). Для вязкости смеси использована формула Манна.
Длина пути смешения I в (2.1) вычисляется по формуле
1_ = 0 11 — ехр(—20ку/6)
6 ’ 1+ехр(—15ку/5)'
(2.2)
Здесь 6 - толщина пограничного слоя. Толщина пограничного слоя определялась из условия Н = 0, 995Н1.
Значение параметра Кармана к определяется с учетом поправки на кривизну обтекаемой поверхности. Аппроксимация опытных данных дает для выпуклой поверхности следующую зависимость к от радиуса кривизны Е
6
к = 0, 4 — 1, 97
Е
Для определения параметра Пк применяется аппроксимация по локальному числу Рейнольдса Евв = реиев/це, построенному по толщине потери импульса
Пк = 10 + Евв — 4) + ^(^ Евв — а1)2, а1 = 3, 2; Ь1 = 57; А1 = 3;
Ь1 =0 (^ Евв > а1), А1 = 0 (^ Евв > 4).
Анализ этой формулы показал, что она дает уточнение лишь при малых Евв. В общем случае параметры а1,Ъ1,А1 зависят от величин, определяющих переход пограничного слоя, и могут быть связаны со значениями Евв в точке перехода, найденными из эксперимента [8].
3. Метод решения
Решение уравнений (1.1)—(1.4) проводится в преобразованных переменных типа Дородницына
1 ГУ Г Уа
х = £, п = д Уо рАу, А = 0 рАу, ф = т£Ди*/ь
и = и. «) Щ.е = дп,т = т. К) М (3Л)
Функцию тока ф, поделенную на 2п, определим из системы уравнений, заменяющих первое уравнение системы (1.1)
дф дф
дХ = —рт£у’ ду = рт£и-
Функции и.(О, Т.(£) в расчетах принимались равными соответственно иоо, V^, при этом особенности, возникающие в коэффициентах уравнений при £ = 0, разрешаются.
Вдув газа с поверхности тела является эффективным средством снижения теплового потока, поэтому возникает задача определения наилучшего способа организации вдува. Для моделирования различных
способов вдувания газа проведено исследование тепло- и массообмена на поверхности гиперболоида. Расход вдуваемого газа задается в виде [6]
G = Go exp ((ж - xm)/xo)2
(3.2)
Здесь x - продольная координата, xm - продольная координата точки поверхности гиперболоида, в которой тепловой поток принимает максимальное значение, xo - фиксированное значение.
Так как на практике величина суммарного расхода вдуваемого газа всегда ограничена
Г хк
Gs = 2п / G (x) rw (x) dx,
Jo
наилучший закон распределения вдува (3.2), обеспечивающий снижение максимального значения теплового потока к телу, целесообразно искать при фиксированном значении Gs. В процессе вычислений условие Gs = const после задания xo обеспечивается выбором величины Go, определяемой по формуле
Go = GsI ^2n Jo exp ((x - xm) /xo)2 rw (x) dx
Данная зависимость моделирует различные способы организации вдува. G (x) ^ Go = const при больших xo, что соответствует постоянному вдоль поверхности тела расходу вдуваемого газа. При малых xo вдув сосредоточен в небольшой окрестности точки xm, вне которой он практически равен нулю. Случай Go = 0 соответствует непроницаемой поверхности.
При численном интегрировании краевой задачи (1.1)—(1.3) использован метод, основанной на применении конечно-разностной схемы четвертого порядка точности по нормальной координате к поверхности тела и первого по продольной, разработанный для численного решения уравнений параболического типа.
Полученная краевая разностная задача решалась методом прогонки с итерациями. Линеаризация разностных уравнений проводилась по значениям функций, полученным на предыдущей итерации. Уравнения
(1.1) решались поочередно. Вначале определялся градиент давления Р2 = Ц— Щ, входящий во второе уравнение системы (1.1), из уравнения, которое получается дифференцированием по переменной x уравнения импульсов в проекции на нормаль. В переменных Дородницына
(3.1) для функций р и р2 получаются обыкновенные дифференциальные уравнения 1-го порядка
дР = кА и* (! V
дп дп )
др2 1 д кАи0 ( д!
дп -иу д 8 ' ^ ^ОО 1 дп
(3.3)
Уравнения (3.3) интегрировались от ударной волны до тела с использованием квадратурной формулы Симпсона. Граничное условие для функции р2 получается дифференцированием условия для функции р на ударной волне.
Последовательность численного интегрирования системы (1.1) на каждой итерации была следующей. Задавалось линейное начальное приближение на критической линии. Затем определялись профили давления р и градиента давления р2 из (3.3). Плотность р определялась из уравнения состояния. Полная вязкость определялась по формулам
(2.1)—(2.2), молекулярная вязкость у - из закона Сазерленда. Затем интегрировались последовательно уравнения для приведенной функции тока f, концентрации вдуваемого газа р и температуры в, подсчитывалось новое значение отхода ударной волны в переменных Дородницына А. Для нахождения А использовалось первое граничное условие в системе (1.2). Итерации на каждом луче продолжались до тех пор, пока максимальное отличие всех профилей и параметра А на данной итерации от предыдущей не становилось меньше заданной точности. Начальное приближение на последующих лучах задавалось с предыдущего слоя.
После решения разностных уравнений вычислялись распределения теплового потока дш, коэффициента трения с/, коэффициента массо-обмена в (тепловой поток отнесен к р00У3, сила трения - к рооУ*, коэффициент массообмена - к росУос ) по формулам
^ да/ = дс.
Для уравнений движения тела по траектории и его полного аэродинамического нагрева важную роль играют интегральные характеристики - коэффициент полного теплового потока
1 дТ
сн =-^ I ЯшйБ, ^ Е дутт2 , (3.4)
Бы -]я О.Ьроои2
коэффициент волнового сопротивления
со = 7^ [ рш айБ, рш = Рш , (3.5)
Ьы ЛЗ О.Ьроои2
коэффициент сопротивления трения
1 f duu
Ст = -^ Tw cos adS, Tw = ——, (3.6)
SM Js Q.OpooU2
а также коэффициент полного массообмена
1 г D dc
Cm =-^ l^wdS, pw = S dy2 . (3.7)
Sm Js 0.5pooU2
Здесь Sm и S - площади миделя тела и его боковой поверхности. Поверхностные интегралы в формулах (3.4)-(3.7) после несложных преобразований вычислялись по формулам, записанным на примере преобразования формулы (3.6)
2п С1
ст = —— rw (ж) tw (ж) cos adx. (3.8)
Sm Jо
Криволинейные интегралы типа (3.8) вычислялись по формуле трапеций.
4. Результаты численного решения
Численные расчеты проводились на неравномерной сетке со сгущением узлов в поперечном направлении к поверхности тела, в продольном - в области затупления.
Уравнение образующей гиперболоида в цилиндрической системе координат, связанной с носком тела, имеет вид
г2 г2
2 2 1,
Р1 я{
где г, г - координаты вдоль оси симметрии и по нормали к поверхности, Р1, - параметры. Все длины отнесены к радиусу носового затупле-
ния Но. Изменением параметров р1, Я1 можно менять эффективный раствор гиперболоида.
Безразмерные параметры а, Ь, с, характеризующие неравномерность профилей скорости и температуры набегающего потока, по оценкам [3], лежат в диапазоне 0 < а < аи, Ь > 0, 3 < с < 5. При любых а < аи не будет наблюдаться возвратно-циркуляционного течения у лобовой поверхности затупленного тела. Рассматривалось обтекание серии гиперболоидов неравномерным потоком вязкого газа (М0 ^ то).
Решение получено при следующих значениях параметров задачи: Не = 4 ■ 105 ^ 107; а = 0; 0, 02; 0, 04; Ь = 2, 6; С = 4; а = 0, 7; ат = 0, 9; Бс = 0, 7; Бст = 0, 9; Хш = 0,15; и = 0, 5; 72 = 1, 67 (вдув гелия); р1 = 1; д1 = 0, 4 ^ 2.
Результаты расчетов обтекания гиперболоидов показывают, что при турбулентном течении профили продольной компоненты скорости становятся более наполненными как в критической точке, так и на обводе гиперболоида. При этом температура газа и концентрация вдуваемого газа внутри ударного слоя повышаются. Максимумы теплового потока, коэффициентов трения и массообмена расположены на боковой поверхности гиперболоида и существенно превышают аналогичные величины при ламинарном режиме течения.
Для снижения максимальной температуры поверхности гиперболоида осуществлялся неравномерный вдув газа с поверхности пористого затупления. На рис. 1-3 приведены распределения теплового потока, коэффициентов трения и массообмена для Не = 2 ■ 106; а = 0, 02 при вду-ве газа в ударный слой по закону (3.2) при фиксированном суммарном расходе вдуваемого газа Сц = 0, 0029.
3
4*
1,5
0
2 $ 4
Рис. 1. Распределение теплового потока вдоль поверхности гиперболоида
Кривая 1 соответствует равномерному вдуву газа при Со = 0, 00027, кривые 2-4 - значениям параметра Хо = 0, 3; 0, 4; 1, 2. Заданным хо соответствуют значения С0 = 0, 00717; 0, 00543; 0, 00186.
Как видно из рисунков, с уменьшением параметра Хо вдув газа в точках максимального значения теплового потока становится существенным. С увеличением хо тепловой поток, коэффициенты трения и мас-сообмена возрастают на боковой поверхности гиперболоида.
О -I1-----------------------------------------------
1,5 С 3
Рис. 2. Распределение коэфф. трения вдоль поверхности гиперболоида
Рис. 3. Распределение коэфф. массообмена вдоль поверхности гиперболоида
Анализ результатов показал,что зависимость максимального значения теплового потока от параметра Хо является немонотонной. Следовательно, существует распределение вдува газа (3.2) с фиксированным суммарным расходом, при котором распределение теплового потока вдоль поверхности гиперболоида является наилучшим. С точки зрения теплозащиты вдув газа по закону (3.2) является более эффективным, чем вдув с равномерным распределением расхода при фиксированном значении суммарного расхода газа.
Исследовалось влияние числа Рейнольдса, неравномерности набегающего потока и эффективного раствора гиперболоида на аэродинамические характеристики. С ростом числа Рейнольдса увеличиваются тепловой поток, коэффициенты трения и массообмена как в области затупления, так и на боковой поверхности гиперболоида. При этом с
ростом неравномерности набегающего потока наблюдается уменьшение локальных аэродинамических коэффициентов. При увеличении эффективного раствора гиперболоида значительно уменьшаются тепловой поток, коэффициенты трения и массообмена в области затупления.
При исследовании движения тела по траектории и его аэродинамического нагрева важную роль играют интегральные аэродинамические характеристики. Из анализа расчетов следует, что с ростом неравномерности набегающего потока наблюдается уменьшение коэффициентов полного теплового потока, волнового сопротивления, сопротивления трения и полного массообмена. При росте числа Рейнольдса наблюдается уменьшение интегральных коэффициентов. Следует отметить, что зависимость коэффициентов полного теплового потока, сопротивления трения и полного массообмена от числа Рейнольдса носит монотонный характер. При этом зависимость коэффициента волнового сопротивления от числа Рейнольдса монотонно снижается и выходит на значение, соответствующее расчету коэффициента волнового сопротивления по модели невязкого газа.
Список литературы
1. Гершбейн Э. А. Теория гиперзвукового вязкого ударного слоя при больших числах Рейнольдса и при сильном вдуве инородных газов / Э. А. Гершбейн // Приклад. математика и механика. - 1974. - Т. 38, вып. 6. - С. 99-105.
2. Головачев Ю. П. Численное моделирование течений вязкого газа в ударном слое / Ю. П. Головачев. - М. : Наука. Физматлит, 1996. - 376 с.
3. Исследование аэродинамических характеристик и теплообмена тел в неравномерном сверхзвуковом потоке газа / И. Г. Еремейцев, Н. Н. Пилюгин, И. С. Хлебников, С. А. Юницкий. - М. : Изд-во МГУ, 1988. - 106 с.
4. Еремейцев И. Г. Исследование турбулентного течения в вязком ударном слое при обтекании газом затупленных удлиненных тел / И. Г. Еремейцев, Г. С. Журавлева, Н. Н. Пилюгин // Приклад. механика и техн. физика. - 1993. - № 1. - С.69-75.
5. Журавлева Г. С. Влияние вдува газа с поверхности сферы на трение и теплообмен при неравномерном турбулентном гиперзвуковом обтекании / Г. С. Журавлева, Н. Н. Пилюгин // Теплофизика высоких температур. -1999. - Т. 37, № 3. - С.427-433.
6. Журавлева Г. С. Влияние вращения тела на теплообмен и сопротивление удлиненного затупленного тела в неравномерном сверхзвуковом потоке газа / Г. С. Журавлева, Н. Н. Пилюгин // Внутрикамерные процессы, горение и газовая динамика дисперсных систем : 5-я междунар. школа-семинар : сб. материалов. - СПб. : Изд-во БГТУ, 2007. - Т. 2. - С. 62-69.
7. Журавлева Г. С. Турбулентный ударный слой на вращающихся затупленных телах, обтекаемых неравномерным сверхзвуковым потоком вязкого газа при наличии вдува / Г. С. Журавлева, Н. Н. Пилюгин // Внутрикамерные процессы, горение и газовая динамика дисперсных систем : 4-я междунар. школа-семинар : сб. материалов. — СПб. : Изд-во БГТУ, 2005. — Т. 1. - С. 119-128.
8. Котляр Я. М. Методы и задачи теплообмена / Я. М. Котляр, В. Д. Совершенный, Д. С. Стриженов. - М. : Машиностроение, 1987. - 302 с.
9. Пилюгин Н. Н. Численное исследование неравномерного обтекания сферы в рамках модели вязкого ударного слоя / Н. Н. Пилюгин, Р. Ф. Талипов // Приклад. механика и техн. физика. — 1991. - № 5. - С. 62-67.
10. Тирский Г. А. К теории гиперзвукового обтекания плоских и осесимметричных затупленных тел вязким химически реагирующим многокомпонентным потоком газа при наличии вдува / Г. А. Тирский // Науч. тр. ИМ МГУ. - М. : Изд-во МГУ, 1975. - № 39. - С. 5-38.
11. Тирский Г. А. Сравнение моделей тонкого и полного вязкого ударного слоя в задаче сверхзвукового обтекания притупленных конусов вязким газом / Г. А. Тирский, С. А. Утюжников // Прикладная математика и механика. -1989. - Т. 53, вып. 6. - С. 963-970.
G. S. Zhuravleva
Numerical model operation turbulent flows at hyperacoustic flowstreamlining of skew fields
Abstract. Hyperacoustic flow of the blunt axisymmetrical skew fields by a stream of a viscous compressible gas is considered. The solution of a problem is obtained by a finite-difference method for a case of a variable injection on contour of spherical bluntness of a skew field in a broad band of Reynold’s numbers.
Keywords: equations of Navier-Stokes; sinuous flow; shock wave; numerical methods.
Журавлева Галина Степановна, кандидат физико-математических наук, доцент, Институт математики, экономики и информатики, Иркутский государственный университет, 664003, Иркутск, ул. К. Маркса, 1, тел.: (3952) 24-22-14 (zhur@math.isu.ru)
Zhuravleva Galina, Irkutsk State University, 1, K. Marks St., Irkutsk, 664003, associated professor, Phone: (3952) 24-22-14 (zhur@math.isu.ru)
