Расчет вязкого ударного слоя около поверхности затупленных тел с использованием алгебраической модели турбулентности Текст научной статьи по специальности «Физика»

УДК 533.6.011.5:532.582.33

Д. А. З а б а р к о, В. П. К о т е н е в, И. А. Шлякова

РАСЧЕТ ВЯЗКОГО УДАРНОГО СЛОЯ ОКОЛО ПОВЕРХНОСТИ ЗАТУПЛЕННЫХ ТЕЛ С ИСПОЛЬЗОВАНИЕМ АЛГЕБРАИЧЕСКОЙ МОДЕЛИ ТУРБУЛЕНТНОСТИ

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

E-mail: d9454@mail.ru

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

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

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

Постановка задачи. Для моделирования турбулентного обтекания используются уравнения Рейнольдса осредненного турбулентного движения. Применяя гипотезу Буссинеска, согласно которой (турбулентные сдвиговые напряжения для плоских течений, например, связаны со скоростью средней деформации через скалярную турбулентную вязкость (т = ит—, где иТ — коэффициент турбулентной (виду

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

движение, можно записать выражение полного касательного напряжения трения, понимая под последним как ламинарное (молекулярное),

так и турбулентное трение, в виде т = (лл + лт) —. Таким образом,

дУ

коэффициент динамической вязкости л представляет собой сумму молекулярной составляющей вязкости ¡л и турбулентной (вихревой) вязкости ¡т [5].

Для получения стационарного решения задачи обтекания тела равномерным сверхзвуковым потоком используется метод установления решения по времени для нестационарной системы осредненных уравнений Навье-Стокса [5]

р— = — grad р — ^ grad (^¡л divVj + 2Б1у (^рБ^ . Уравнение неразрывности для смеси имеет вид

др + (р9) =0.

Уравнение энергии записывается как

- Л. V2 )

р-Т \Н + т] =

др ( ( V2 1 ^ -Л 2 = —- + div I 2л I grad —---V х го^ I — -¡V div V — д

дТ \ \ 2 2 /О

В приведенных уравнениях Н — энтальпия газа; р — давление; р — плотность; V — вектор скорости; р = ¡л + ¡т — динамическая вязкость; V — вектор плотности теплового потока; Б — тензор скоростей деформации; Бгу — дивергенция тензора.

Тепловой поток к поверхности тела обусловлен молекулярным и турбулентным переносом теплоты:

—V = (Р + £) ^ н,

ЛСр ¡¡т Ср _

где Рт = —— = 0,72 — число Прандтля; Ргт = —— = 0,9 — турбу-

Л Лт

лентное число Прандтля [1]; ср — удельная теплоемкость смеси при

постоянном давлении; Л, Лт — коэффициенты молекулярной и турбулентной теплопроводности соответственно.

_ л р р

Введем безразмерные величины: ¡безр = —; рбезр = —; рбезр = —;

Рж рж рж

т = _т_. V = V • = х • = у • = 1 •

т безр гр ; —безр у---; хбезр т ; убезр т ; ^безр т ;

тж VPж/Pж ь ь ь

где L — характерный линейный размер. Далее индекс

"безр" будет опущен.

Полученную систему уравнений (здесь не приводится) дополняет безразмерное уравнение состояния в виде

Сформулируем основные положения, при которых решается задача:

1) набегающий на затупленное тело поток газа сверхзвуковой, однородный и невозмущенный с показателем адиабаты 7 = 1,4;

2) течение во всей возмущенной телом области симметрично относительно вертикальной плоскости, проходящей через ось симметрии тела, условия симметричности течения имеют вид

где и,у,т — компоненты вектора скорости V, переменная р соответствует координатной линии, перпендикулярной плоскости симметрии;

3) осредненные уравнения Навье-Стокса справедливы для описания течения во всей возмущенной телом области;

4) из вязкостных эффектов учитываются молекулярная и турбулентная вязкость.

В качестве граничных условий на поверхности обтекаемого тела для компонент скорости задаются условия прилипания ии = = = = 0 и температура стенки Ти. (Индекс ы соответствует параметрам течения газа, задаваемым на стенке.)

Для обеспечения устойчивости расчета налагается условие постоянства давления поперек ударного слоя в непосредственной близости

При задании тонкой головной ударной волны ("внешней" границы области интегрирования) необходимы априорные данные о структуре поля течения. Для больших чисел Рейнольдса (Яе^ > 104), характерных для рассматриваемых турбулентных режимов течений, можно пренебречь влиянием структуры тонкой головной ударной волны на течение вниз по потоку и принять, что ударная волна является поверхностью разрыва газодинамических параметров, на которой выполняются нестационарные условия Рэнкина-Гюгонио, записываемые в безразмерной форме как

p = рТ.

du dv др дТ дф дф дф дф

= w = 0,

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

0.

р (Vn - D) = - D); p + р (Vn - D)2 = 1 + (Vn,TO - D)2 ;

h (p,p) + 2(Vn - D)2 =

2 y — I

к = V™,

12

+ — D)2 ;

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

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

В турбулентном пограничном слое выделяется по меньшей мере пять областей: вязкий подслой, переходная (или буферная) область, область логарифмического профиля скорости, область закона следа, область перемежаемости [2]. Первые три принято объединять в одну внутреннюю область или область закона стенки (рис. 1). Внутренняя область пограничного слоя занимает примерно 15-20 % толщины всего слоя. Согласно измерениям, в ней генерируется до 80% энергии турбулентности, причем первые 5 % толщины дают более половины вклада в полное производство турбулентной энергии. Область закона следа и область перемежаемости обычно объединяют во внешнюю область турбулентного пограничного слоя, которая занимает порядка 80% толщины всего слоя.

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

¡т = р12 5, (1)

'/ди\ 2 [ ду \ 2 где для плоских течений [2] 5 = \1 ( — ) + ( — ) , а для полно-

\ду) \дх)

стью трехмерных сдвиговых слоев величина 5 эквивалентна вектору

Рис. 1. Схематическая структура ударного слоя:

1, 2 — внешняя и внутренняя области пограничного слоя соответственно

завихренности:

д^ ди\2 /ды дv\2 /дм дых 2 \дх ду) \ду дг) \дг дх

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

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

ди

можно интерпретировать как характерную

[1]. Произведение I 0

дУ

скорость турбулентности V..

При вычислении I для внутренней области, расположенной в непосредственной близости от стенки, хорошие результаты дает оценка согласно формуле [1]

1 = ХУ

1 - e a*

(2)

Здесь 1 — е а* — демпфирующая функция Ван-Дриста, необходимая для связи между полностью развитым пограничным слоем (I = ху) и вязким подслоем (I ^ 0); х = 0,41 — постоянная Кармана; Л* = 26 — демпфирующая константа; у+ = ит— сеточный параметр, характеризующий качество сетки в пристеночной области течения и

являющийся безразмерной величиной; ит = \ /1—^1 — постоянная для

V Р'Ш

луча динамическая скорость; = — — кинематическая вязкость;

Р'Ш

— динамическая вязкость; рп — плотность на поверхности тела;

ди

= — — касательное напряжение трения; индекс ы соответ-

V ду) п

ствует параметрам на стенке.

Для расчета турбулентной вязкости во внешней области пограничного слоя используется формула Клаузера [1]

= ар8*ие/1,

где ие — скорость на границе пограничного слоя; параметр а ~ 0,0168

(его модификация позволяет учитывать эффекты при малых числах

5

Рейнольдса); = (1--) ду — толщина вытеснения; I = 1 +

\ ие]

0

+ 5,5(у/8)6 — аппроксимационная величина; 8 — толщина пограничного слоя; индекс е соответствует параметрам на внешней границе пограничного слоя.

Для окончательного вычисления коэффициента турбулентной вязкости из двух коэффициентов, полученных для внутренней и внешней областей турбулентного пограничного слоя, берется минимальной значение ¡т = шт(¡т(внутр), ¡т(внеш)) .

При зарождении турбулентного пограничного слоя размеры как внутренней, так и внешней областей стремятся к нулю и применение модели турбулентности, использующей выражение (2), требует обоснования. Трудности связаны с тем, что малые значения 8 при зарождении турбулентного пограничного слоя, вызывают переключение на модель внешней области, прежде чем демпфирующий эффект позволит развиться полностью турбулентной области закона стенки. Это приводит к тому, что разностные схемы, использующие эту модель, дают уменьшенные значения касательного напряжения на стенке. Такое уменьшение мало для течений несжимаемой жидкости, оно гораздо больше для сжимаемой жидкости и становится все более заметным при увеличении числа Рейнольдса, так как число Маха растет из-за утолщения вязкого подслоя вследствие тепловых эффектов, на что влияет интенсивность охлаждения стенки в случае сжимаемых течений [1].

Можно добиться хорошего соответствия результатов расчета с экспериментальными данными при малых числах Рейнольдса, вычисленных по толщине потери импульса, путем простого запаздывания переключения с модели внутренней области слоя на модель внешней области. Если 1/8 ^ 0,09, то в модификации нет необходимости. Если соотношение (2) дает 1/8 > 0,09, то длина пути смешения I ограничивается искусственно (/тах = 0,098) и становится постоянной для вихревой вязкости, рассчитываемой по соотношению (1), во внешней области турбулентного пограничного слоя. Такая простая модификация дает линейно-логарифмический закон скорости, что согласуется с экспериментом.

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

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

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

Определение д — толщины пограничного слоя. Если скорость в невязкой части ударного слоя сверхзвукового течения не постоянна, то обычное определение границы пограничного слоя теряет смысл. В невязкой части ударного слоя должна сохраняться постоянной полная энтальпия. В случае охлаждаемой стенки, когда основное изменение полной энтальпии имеет место в пограничном слое, его граница определяется из условия [7]

\Н — Нте\

I Hqq Hw i

= 0,005.

Полная энтальпия увеличивается от значения Hw на стенке до постоянного значения Нж в невязком потоке и невозмущенном течении. Вследствие погрешности аппроксимации профиль полной энтальпии может иметь локальный максимум в той области, где определяется граница пограничного слоя. Было установлено [7], что при мелком (вдоль радиального направления) шаге сетки такой локальный максимум Н не возникает. Таким образом, нахождение 5 на основе приведенного соотношения ограничивалось областью, расположенной ниже сечения, где располагался этот максимум.

Математическая модель. Для моделирования обтекания затупленных тел используем сферическую систему координат (ССК) R, р, 9 (рис. 2). Переход от декартовой системы координат x,y,z в исходную сферическую задается преобразованием координат

x = R sin 9 sin р, y = — R sin 9 cos р, z = z0 — R cos 9,

где z0 — смещение по оси z центра ССК относительно начала декартовой системы координат.

Введем новую переменную £ в следующем виде

£ = R (R-RT (R,PL ), (0 < £ < 1, 0 < р < п, 0 < 9 < п/2), Rb (t, 9, р) — Rt (9, р)

где R = Rt (9, р) + £ (Rb М,р) — Rt (9,р)), а RT = Rt (9,р), RB = RB (t, 9, р) — уравнения, задающие поверхности тела и ударной

Рис. 2. Исходная сферическая система координат

Рис. 3. Нормированная сферическая система координат волны соответственно, причем Яв определяется в процессе решения задачи обтекания.

Расчет течений около затупленных тел будем проводить в нормированной сферической системе координат £, р, 9 (рис. 3), полученной из исходной ССК.

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

В локальной системе координат компоненты скорости определяются соотношениями V = V • п для нормальной составляющей и и = \/V2 — V2 для касательной составляющей.

Уравнение для нормали к поверхности Я — Ят + £ (Яв — Ят) =0 в нормированной сферической системе координат записывается в виде

Яту + (Яту — Яту) £ Яте + (Яве — Яте)

n =

R sin в

R

1 +

R

'■Tip + (rtp — RTp) £\ 2 + f RTe + (Rße — RTe) £

R sin в

R

нижние индексы 9, р обозначают дифференцирование по соответствующей переменной.

2

Рис. 4. Локальная система координат, связанная с обтекаемым телом

Для поверхности тела, где £ = 0, имеем

Нт^ Нте

ey = n =

R sin в R

1 +

R

Т v

R sin в

+

R

Т в

R

еу ,ех — единичные орты в локальной системе координат; еу = п — нормаль к поверхности тела.

Сам вектор скорости определен в обеих системах координат — исходной ССК и локальной системе координат:

V = ивх + vey = ивев + и^е^ + щ ее,

где ев,ее,еф — единичные орты в сферической системе координат; ив,ие,иф — компоненты скорости в исходной ССК.

Запишем выражения для орта в касательном направлении

V — vn

ex

v- [V-fDn

V — vn

V- V-n)n

и градиента функции f в сферической системе координат

1

f = dv e

f = öReR + R sin в dp

dv ^ 1 dv ^ •- ev + R • двев ■

Производная компонент скорости в локальной системе координат вычисляется по формулам:

1

dv dv

dx = dR (VR ^ ex) + R sin в dp

dv 1 dv

—JeV • ex) + R • de (ee • ex)

2

2

по касательному направлению и du du 1 du

dy dR по нормали.

du = ^• ) + R—e • •ey) +- ^ ® • )

1 du

R • de

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

Этап 1. Интегрируется "невязкая" часть системы, включающая в себя уравнения движения Эйлера, неразрывности и энергии только с невязкими членами:

dU ^ = к •

где U =

Р

V , к =

, V2

h +

2

—p2div V —grad p

dp dt

"Невязкая" часть газодинамической системы уравнений решается с использованием явной двухшаговой конечно-разностной схемы Мак-Кормака.

Введем в пространстве t, р, в сетку с шагами At, Д£, Др = п/К, Дв и обозначим координаты узлов сетки как V = ]Дt, вт = тДв, = кДр, а любую функцию f в точке как ¡П,тк. При этом 0 < п < N, 0 < т < М, 0 < к < К, где Ы,Ы,К'— число интервалов в направлениях в, р соответственно.

Этап 2. Решение уточняется с использованием "вязкой" части уравнений Навье-Стокса и энергии:

dU ^ q

Р^тт = K, где U = dt

Р V

V2 h + Т

к =

— - grad (p div V) + 2Div (p S

V2 1 ^ 2 ^ ^ J

div ( 2p ( grad —--- V x rot Vi — - pV div V — q

0

"Вязкая" подсистема интегрируется при помощи метода прогон-

п,ш,к п,ш,к тЗ-'1+2 тт / ■ , 1 \

ки: -—- = КПараметры на (] + 1)-м слое считаются

известными из решения "невязкой" части уравнений на первом этапе.

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

дГ л д / дГ\ л д / дГ\ л д ( дГ\

Иъ = Ллд£ ) + лфд£ {^ж) + Лед£ +

+ Вм ^ + вдГ + сг + В,

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

Для получения решения на (] + 2)-м слое при заданных граничных условиях интегрирование одномерных уравнений сводится к последовательному решению отдельных разностных уравнений с трехдиаго-нальной матрицей методом прогонки. Эти уравнения имеют стандартный вид

Лп/п+1 — Вп/п + Сп/п-1 = Dn, 0 <п < N,

где в качестве целевой функции / рассматривается любая из неизвестных функций; /0 и /^ заданы или определяются из граничных условий.

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

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

Рис. 5. Разделение физической области интегрирования

£ (n) =

сеток. Для адекватного учета поведения параметров в областях больших градиентов вводится сгущение сетки в радиальном направлении вида

—, если п> N

11 11

+ ^^^ -П(2П ~ , если п < N„

2 2 V 2N, у' *'

где Л — произвольное целое число, задающее отношение наибольшего шага к наименьшему; N — номер сеточного узла, с которого вводится неравномерная сетка.

Для исследования течений около удлиненных тел используется разделение всей области интегрирования на ряд взаимно перекрывающихся подобластей (рис. 5) и проводится последовательный расчет в каждой из них. Сначала задача решается в окрестности затупления. Затем центр сферической системы координат переносится по оси тела. Выстраивается новая расчетная область, где на левой границе задаются "жесткие" граничные условия из уже рассчитанной области, на выходных границах — "мягкие" граничные условия вида линейной экстраполяции искомых функций; на теле задаются условия прилипания и температура стенки, на ударной волне — нестационарные соотношения Рэнкина-Гюгонио. Решение в полученной области устанавливается, и процедура построения новой расчетной сетки повторяется. Описанное разделение области возможно в силу слабой передачи возмущений вверх по потоку при обтекании тел сверхзвуковым набегающим потоком вязкого газа [4] и позволяет проводить расчет длинных затупленных тел.

Из-за сложной структуры уравнений Навье-Стокса невозможно получить аналитическое выражение условия устойчивости для схемы Мак-Кормака. Однако практическое применение численного метода показало возможность без потери устойчивости счета использовать

эмпирическую формулу

АЬ = а Д^кфл ,

где а — коэффициент запаса (а = 0,8); Аьькфл определяется по критерию Куранта-Фридрихса-Леви для линейных гиперболических уравнений в частных производных [1]. Перед очередным шагом по времени для всех точек сетки рассчитывается шаг интегрирования АЬ. Затем наименьшее значение шага используется для получения решения на следующем временном слое.

Анализ результатов. Для проверки адекватности алгебраической модели турбулентности реальным процессам было проведено сравнение результатов моделирования осесимметричного обтекания сферически затупленных конусов с экспериментальными данными [6, 7]. Значения параметров для рассмотренных режимов приведены в таблице. Выбор данных режимов обусловлен наличием доступных экспериментальных данных.

Параметры экспериментальных режимов

№ Т.Х,; K а, град ©, град Re^, м-1 Rn, м Tw, к

1 5 74 0 9 60-106 0,0635 103

2 8 54,5 0 7 12,1 • 106 0,0127 300

Первый режим соответствует осесимметричному обтеканию сферически затупленного конуса с углом полураствора в = 9°, радиусом носка Яп = 6,35 см и общей длиной 10 калибров. Параметры набегающего потока заданы числом Маха М^ = 5, числом Рейнольдса Яе^ = 60 • 106 м-1 и углом атаки а = 0°. Температура поверхности тела принималась постоянной и равной 103 К.

Для указанного режима при моделировании ламинарных и турбулентных течений использовались различные расчетные сетки. Для ламинарного режима число узлов поперек ударного слоя принималось N = 80, сгущение сетки вводилось в обоих направлениях (к телу и к ударной волне) с соответствующими параметрами А1 = 500, N = 60, А2 = 5; параметры сетки по координате в: М = 46, Ав = 2°. Для турбулентного режима первая сетка задавалась значениями параметров N = 160, А1 = 500, N = 160, А2 = 1; вторая — N = 160, А1 = 1000, N = 160, А2 = 1. Число итераций варьировалось в пределах 50 000— 150 000 до достижения сходимости. При этом затраты машинного времени (для ПЭВМ с тактовой частотой процессора ~ 2 ГГц) на моделирование обтекания тел удлинением до 100 калибров составляли от 10 до 30 ч на один вариант расчета в зависимости от параметров

6JE+C5

/

3

О

1

2

3

4

5

6

7

S

s/R,

Рис. 6. Распределение размерного теплового потока по длине образующей тела:

= 5, = 60 • 106 м-1; 1, 2 — расчет для случая турбулентного обтекания, соответственно сетка 1 и 2; 3 — расчет для случая ламинарного обтекания; 4 — экспериментальные данные для случая турбулентного обтекания [6, 7]

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

Анализируя результаты расчетов (рис. 6), можно отметить, что развитая турбулентность наступает уже на сфере на расстоянии примерно 0,6 калибра по длине образующей 8 (расстояние вдоль поверхности тела), в этом месте заметно резкое увеличение теплового потока (до 2,5 раз) по сравнению с критической точкой. В целом можно отметить удовлетворительное согласование (различия не превышают 30%) расчетных и экспериментальных данных, учитывая погрешность эксперимента в пределах 10-20%. Здесь же приведены результаты расчета теплового потока в ламинарном режиме (нижняя кривая). Сравнение показывает, что турбулентный тепловой поток намного превышает ламинарный, а его максимум смещается по потоку от точки торможения.

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

Для второго рассмотренного режима число Маха задавалось равным 8, длина тела — 70 калибрам, угол полураствора затупленного

конуса равен 7°, радиус затупления 0,0127 м, температура поверхности тела 300К. Число Рейнольдса для данного режима 12,1 • 106 м-1.

График зависимости числа Стантона (безразмерный тепловой поток) от координаты, отсчитываемой от носка тела вдоль оси конуса приведен на рис. 7, причем число Стантона определяется по формуле

St =

Р»У» (Ь (Тг) - К)' где — тепловой поток; Ь,ш — энтальпия газа на стенке; Ь (Тг) — энтальпия газа в пограничном слое, соответствующая температуре

7 — 1

восстановления Тг = Тх ( 1 + г—

М2^ . Коэффициент восстановления г для ламинарного слоя принимается равным 0,85, для турбулентного — 0,89.

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

На рис. 7 показано качественное различие в распределении тепловых потоков вдоль поверхности тела. Так, для турбулентного погра-

З.ОЕ 03

2.5Е-03

2.0Е-03

ß 1.5Е-03

5.0Е-04

0.0Е+00

! U 4 /

* т / 1 2

ч • ••/• ' • »* • « • ■ • •

-г Ч V 3

10

20

30 40

7/Rn

50

60

70

Рис. 7. Распределение числа Стантона (безразмерный тепловой поток) в зависимости от координаты, отсчитываемой от носка тела вдоль оси конуса:

= 8, = 12,1 • 106 м-1; 1 — расчет для случая турбулентного обтекания по модели настоящей работы; 2 — расчет для случая ламинарного обтекания; 3 — расчет для случая турбулентного обтекания по модели Себечи-Смита [7]; 4 — экспериментальные данные для случая турбулентного обтекания [7]

ничного слоя характерно увеличение числа Стантона с увеличением расстояния от носка, для ламинарного течения — уменьшение. Здесь же приведены экспериментальные данные. Общее сравнение с экспериментальными данными для турбулентного режима обтекания показывает удовлетворительное согласование тепловых нагрузок с различием не более 20 %. Причем значения тепловых потоков, рассчитанные в работе [7] с использованием модели Себечи-Смита (кривая 3), равно как и с использованием предложенной модели (кривая 1), оказываются заниженными по сравнению с экспериментальными данными.

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

Далее рассматривается осесимметричное обтекание типичной "об-гарной" формы, имеющей место при уносе массы, которая представляет собой затупление с выступом на критической линии. Режим обтекания соответствует первому тестируемому случаю с параметрами: число Маха М^ = 5, угол полураствора затупленного конуса 9°, радиус затупления 0,0635 м, температура поверхности тела 103 К. Число Рейнольдса для данного режима соответствует 60 • 106 м-1. Для данной геометрии тела наиболее интересно течение около выступа обгарной формы, поэтому рассматривается только область затупления тела в районе критической точки и сопряжения выступа с "обгарной" формой.

На рис. 8 приведено распределение числа Маха и линии тока, описывающее характер течения в области затупления. На картине течения видно образование вихревой зоны в месте сопряжения носка с "об-гарной" формой. При сопоставлении рис. 8 и 9 можно отметить, что образование вихря сопровождается увеличением теплового потока в месте сопряжения носового выступа с затуплением (здесь безразмерная координата вдоль оси конуса г/Яп ~ 0,18).

На рис. 9 приведена соответствующая зависимость относительного теплового потока в ламинарном и турбулентном режимах обтекания. Турбулентный тепловой поток в районе своего максимума превышает значение ламинарного потока на величину от 10 до 50%, и влияние турбулентности приводит к более медленному уменьшению теплового потока ниже по течению по сравнению с ламинарным режимом обтекания.

Выводы. Удовлетворительное согласование расчетных и экспериментальных данных, продемонстрированное на примере затупленного

Рис. 8. Распределение числа Маха и линий тока в ударном слое при обтекании "обгарной" формы:

М™ = 5, Яе™ = 60 • 106 м-1

1.4

1.2

0.8

0.6

1 %

4 \

Л

\ \ \ Л \ 2

0.4

0.2

г/Д.

Рис. 9. Распределение относительного теплового потока:

= 5, = 60 • 106 м-1; 1 и 2 — результаты расчета для случаев турбулентного и ламинарного обтекания

конуса (различие не более 20 %), подтверждает корректность применения алгебраической модели турбулентности для моделирования пристенных турбулентных течений при внешнем обтекании ЛА.

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

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

Работа выполнена по Госконтракту № П608 от 06.06.2009 г. в рамках Федеральной целевой программы "Научные и научно-педагогические кадры инновационной России".

СПИСОК ЛИТЕРАТУРЫ

1.Андерсон Д., ТаннехилДж., ПлетчерР. Вычислительная гидромеханика и теплообмен. - М.: Мир, 1990. Т. 1. - 384 с., Т. 2. - 392 с.

2. Б е л о в И. А. Моделирование турбулентных течений: Учеб. пособие / И.А. Белов, С.А. Исаев. - СПб.: Балт. гос. техн. ун-т, 2001. - 108 с.

3. ЗабаркоД. А., КотеневВ. П. Численное исследование ламинарных течений вязкого химически реагирующего газа около затупленных тел // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. - 2006. - № 1. - С. 77-95.

4. Кокошинская Н. С., Павлов Б. М., Пасконов В. М. Численное исследование сверхзвукового обтекания тел вязким газом. - М.: МГУ, 1980. - 248 с.

5. Лойцянский Л. Г. Механика жидкости и газа: Учебник для вузов. - 7-е изд., испр. - М.: Дрофа, 2003. - 840 с.

6. У и д х о п ф Д ж. Ф., Холл Р. Измерение теплопередачи на затупленном конусе под углом атаки при переходном и турбулентном режимах течения // Ракетная техника и космонавтика. - 1972. - Т. 10. - № 10. - С. 71-79.

7. Ш и р а х и С. А., Т р у м е н К. Р. Сравнение алгебраических моделей турбулентности на примере расчета с помощью параболизованных уравнений Навье-Стокса сверхзвукового обтекания конуса со сферическим носком // Аэрокосмическая техника. 1990. - № 10. - С. 69-81.

Статья поступила в редакцию 30.09.2010

Владимир Пантелеевич Котенев окончил МГУ им. М.В. Ломоносова в 1978 г. Д-р техн. наук, начальник отдела ОАО "ВПК "НПО машиностроения". Автор более 30 научных работ в области прикладной математики, численных и аналитических методов исследования течения газа при обтекании поверхности летательных аппаратов.

V.P. Kotenev graduated from the Lomonosov Moscow State University in 1978. D. Sci. (Eng.), head of department of JSC "VPK "NPO Mashinostroeniya". Author of more than 30 publications in the field of applied mathematics, numerical and analytical methods of investigating gas streams in flows about the surface of blunt flying vehicles.

Дмитрий Александрович Забарко окончил МГТУ им. Н.Э. Баумана в 1995 г. Канд. физ.-мат. наук, старший научный сотрудник ОАО "ВПК "НПО машиностроения". Автор более 10 научных работ в области прикладной математики, разработки численных методов и математического программного обеспечения аэрогазодинамических расчетов.

D.A. Zabarko graduated from the Bauman Moscow State Technical University in 1995. Ph. D. (Phys.-Math), senior researcher of JSC "VPK "NPO Mashinostroeniya". Author of more than 10 publications in the field of applied mathematics, development of numerical methods and mathematical software for aero-gasdynamical of computations.

Ирина Андреевна Шлякова окончила МГТУ им. Н.Э. Баумана в 2010 г. Инженер ОАО "ВПК "НПО машиностроения".

I.A. Shlyakova graduated from the Bauman Moscow State Technical University in 2010. Engineer of JSC "VPK "NPO Mashinostroeniya".