CC BY
14
-►
Математическое моделирование физических процессов
DOI: 10.18721/JPM.12201 УДК 519.6:533.6.011
ЧИСЛЕННОЕ РЕШЕНИЕ ТРЕХМЕРНОЙ ЗАДАЧИ ОБТЕКАНИЯ УСТАНОВЛЕННОГО НА ПЛАСТИНЕ ЦИЛИНДРИЧЕСКОГО ТЕЛА СВЕРХЗВУКОВЫМ ПОТОКОМ ВЯЗКОГО ГАЗА ПРИ М = 2,95 Е.В. Колесник, Е.М. Смирнов, А.А. Смирновский
Санкт-Петербургский политехнический университет Петра Великого, Санкт-Петербург, Российская Федерация
Представлены результаты численного моделирования ударно-волновых и вихревых структур, формирующихся при обтекании сверхзвуковым потоком удлиненного цилиндрического тела, которое примыкает к пластине и пронизывает развивающийся на пластине пограничный слой. Рассмотрен ламинарный режим течения при числе Маха набегающего потока, равном 2,95, и числе Рейнольдса, равном 4000. Решения получены с использованием двух схем для расчета конвективных потоков (ИЬЬ и ЛИБМ). Проведено сравнение полей течения, рассчитанных с применением указанных численных схем первого и второго порядков точности. Обсуждаются вопросы сходимости численного решения по сетке.
Ключевые слова: высокоскоростное течение, вязко-невязкое взаимодействие, численное моделирование, схема ЛИБМ, схема ИЬЬ
Ссылка при цитировании: Колесник Е.В., Смирнов Е.М., Смирновский А.А. Численное решение трехмерной задачи обтекания установленного на пластине цилиндрического тела сверхзвуковым потоком вязкого газа при М = 2,95// Научно-технические ведомости СПбГПУ. Физико-математические науки. 2019. Т. 12. № 2. С. 7-22. БО1: 10.18721/ 1РМ.12201
NUMERICAL SOLUTION OF A 3D PROBLEM ON A SUPERSONIC VISCOUS GAS FLOW PAST A PLATE-CYLINDRICAL BODY
JUNCTION AT M 2.95 E.V. Kolesnik, E.M. Smirnov, A.A. Smirnovsky
Peter the Great St. Petersburg Polytechnic University, St. Petersburg, Russian Federation
In the paper, results of numerical simulation of a shock-wave pattern and vortex structures forming in supersonic flow past an adjacent-to-plate elongate cylindrical body, which penetrates the developing flat-plate boundary layer, have been presented. The laminar flow regime at Mach number 2.95 was considered, Reynolds number was taken 4000. The solutions were obtained using two schemes for convective flux (HLL and AUSM). Comparison of the flow fields calculated with the mentioned schemes of the first and second orders of accuracy were conducted. Solution grid sensitivity issues were discussed.
Keywords: high-speed flow, viscous-inviscid interaction, numerical simulation, AUSM and HLL schemes
Citation: Kolesnik E.V., Smirnov E.M., Smirnovsky A.A., Numerical solution of a 3D problem on a supersonic viscous gas flow past a plate-cylindrical body junction at M 2.95, St. Petersburg Polytechnical State University Journal. Physics and Mathematics. 12 (2) (2019) 7-22. DOI: 10.18721/JPM.12201
Введение
Необходимость проведения исследований структуры течения при обтекании дозвуковым и сверхзвуковым потоками каких-либо конструктивных элементов, закрепленных на обтекаемой поверхности, возникает во многих практических задачах авиакосмической отрасли и турбомашино-строения. В аэрокосмической отрасли задачи данной направленности встречаются, в частности, при конструировании таких элементов летательных аппаратов, как соединения фюзеляжа и крыла, хвостового оперения и других элементов, а также при оптимизации взаимодействия сверхзвукового потока с вдуваемой газовой струей, используемой в элементах управления.
Особый интерес к задачам такого рода возник с середины прошлого века в связи с развитием сверхзвуковых и гиперзвуковых летательных аппаратов. При сверхзвуковом обтекании вязким газом области сочленения затупленного тела с пластиной возникает трехмерное течение со сложной ударно-волновой структурой, обширной отрывной областью и системой подковообразных вихрей. Другими словами, в данной задаче ярко проявляются эффекты вязко-невязкого взаимодействия в сверхзвуковом течении вязкого газа.
В целом за последние годы была получена достаточно полная информация об особенностях течения, возникающего при взаимодействии отрывающегося пограничного слоя с головным скачком уплотнения. Обзоры литературы на данную тему можно найти, например, в работах [1, 2].
Отметим также, что и случай дозвукового обтекания, где тоже имеют значение особенности формирования отрывной области перед телом, представляет большой интерес. Это относится, в частности, к задачам тур-бомашиностроения, где сведения о сложной вихревой структуре течения в окрестности передней кромки лопатки важны для учета особенностей теплообмена в области сочленения лопатки и торцевой стенки в венцах высокотемпературных газовых турбин. Исследованию структуры течения в условиях его дозвукового режима посвящено множество работ как экспериментальных, так и расчетных (см., например, статьи [3 — 6]).
Предметом настоящего исследования является одна из модельных задач, в которой реализуется сложная структура трехмерного потока. Рассматривается сверхзвуковое
обтекание удлиненного цилиндрического тела, установленного на пластине, вдоль которой развивается пограничный слой (рис. 1). Данной задаче посвящено довольно большое число экспериментальных работ [7 — 13]. Уже первые исследования структуры потока в данной конфигурации [7 — 9] показали, что перед цилиндрическим телом образуются локальные сверхзвуковые области и скачки уплотнения, которые приводят к появлению вторичной зоны отрыва внутри основной зоны.
В случае высокоскоростных летательных аппаратов наиболее важной характеристикой является интенсивность теплообмена в области перед обтекаемым телом; результаты исследований показывают, что в этой области наблюдается резкое увеличение теплового потока к стенке. В частности, опубликованная недавно работа [13] посвящена изучению локального теплообмена на пластине при обтекании цилиндрического тела потоком с числом Маха, равным 5.
В экспериментах весьма проблематично получить объем данных, в деталях определяющих структуру трехмерного сверхзвукового потока во всей области его развития. Соответственно, большой интерес представляют возможности предсказания структуры потока в рассматриваемой модельной конфигурации на основе методов вычислительной гидрогазодинамики. При этом особое внимание следует уделять обеспечению высокого качества разрешения
Рис. 1. Схематическое представление постановки задачи:
обтекание вязким газом удлиненного цилиндрического тела с диаметром затупленной части Б, установленного на пластине; вдоль пластины развивается пограничный слой (5 — его толщина)
областей вязко-невязкого взаимодействия. В первых работах, посвященных численному моделированию ламинарного и турбулентного течений в рассматриваемой конфигурации [14, 15], использовались расчетные сетки, которые не позволяли детально разрешать все характерные особенности течения. В относительно недавней работе [16], содержащей результаты обстоятельных исследований (как экспериментального, так и численного) обтекания затупленного тела при числе Маха M = 6,7, показано, например, что приемлемое по точности разрешение структуры потока достигается при размерности сетки порядка 15 млн. ячеек.
Численное моделирование сверхзвукового течения в рассматриваемых условиях может осложняться проявлением так называемой «карбункул»-неустойчивости [17, 18], которая приводит к сильному искажению фронта ударной волны в численном решении. Для подавления данной неустойчивости, которой подвержены некоторые широко известные численные схемы, разрабатываются различные подходы, в частности гибридные схемы [18, 16] и введение добавочной искусственной вязкости [19]. Вместе с тем имеются схемы, при использовании которых карбункул-неустойчивость обычно не наблюдается. Это схема Хартена — Лакса — Ван Лира (HLL) [20], которая характеризуется повышенной диссипатив-ностью, и семейство схем, основанных на расщеплении вектора потоков, предложенные Лио и Стефаном (метод AUSM, англ. Advection Upstream Splitting Method) [21], относительно которых многими авторами отмечалась их устойчивость к возникновению нефизических осцилляций на фронте ударной волны.
В большинстве опубликованных работ, посвященных численному решению задачи о высокоскоростном обтекании вязким газом затупленного тела, установленного на пластине, рассматривается случай гиперзвукового течения. Особенностям же течений с умеренными значениями числа Маха для набегающего потока и вопросам качества численного воспроизведения эффектов вязко-невязкого взаимодействия в литературе уделялось мало внимания.
В настоящей работе представляются результаты численного решения задачи об обтекании удлиненного затупленного тела, установленного на пластине, при числе Маха набегающего потока, равном
2,95, и числе Рейнольдса, построенном по диаметру затупленной части, равном 4000. Выбор умеренного значения числа Рейнольдса, обеспечивающего ламинарный характер течения в рассматриваемый области, обусловлен нацеленностью работы на получение аккуратного, практически се-точно-сошедшегося решения, с детальным разрешением как сложной структуры течения в вязкой отрывной области с системой подковообразных вихрей, так и газодинамической структуры, характеризуемой отошедшим скачком уплотнения и системой косых волн уплотнения, образующихся при обтекании отрывной области.
Численные решения получены по двум схемам: ИЬЬ и Ли8М; обе позволили избежать проявлений карбункул-неустойчивости. В работе также проведен сравнительный анализ решений.
Численный метод решения задачи
Базовые положения. Представленные далее численные решения для течения вязкого газа получены по методу конечных объемов (МКО), который основан на использовании интегральной формулировки законов сохранения массы, импульса и энергии, записанных для расчетной ячейки (контрольного объема):
я. (!)
d = о,
п ОТ м s
"т.
где О — контрольный объем; М — число его граней; Sm — площадь текущей грани, т = 1,М; вектор потоков на грани контрольного объема; w = [р,ри,ру^,рЕ\ — вектор консервативных переменных (и = и V = и = и3 — компоненты вектора скорости V в декартовой системе координат; Е — полная энергия; р — плотность).
Вектор Б*1 является суммой векторов невязких и вязких потоков Е = Е+Еуис, которые определяются выражениями:
рУ • п рУ • пи + рп • 1
Е = рУ • пу + рп • ] , (2)
рУ • п^ + рп • к рУ • пЯ
0
п-т-1
Г"с = п -т- ] , (2)
п-т-к п-(т-V + я)
где р — давление; Н — полная энтальпия; п — нормаль к грани; 1, ^ к — орты декартовой системы координат (х = ху = х z = х3).
Компоненты тензора вязких напряжений т и вектора плотности теплового потока я записываются в виде
ди,
ди, ^ -+ 1
дх, дх
■(2/3)
ди,,
г У
—5.. дхк Ч
Ч1 =
^(дТ/д* 1),
(3)
(4)
где Т — температура; ц — динамический коэффициент вязкости газа, зависимость которого от температуры определяется по формуле Сазерленда; X — коэффициент теплопроводности газа.
Величины полных энергии и энтальпии определяются выражениями:
Е = сТ + V72, Н = срТ + V72,
где с , с — удельные теплоемкости при постоянном объеме и постоянном давлении соответственно; они принимаются постоянными.
Связь между давлением р, плотностью р и внутренней энергией е определяется уравнением состояния совершенного газа:
Р = (У -1)Ре>
где у — показатель адиабаты (у = с /су).
Схемы расчета конвективных потоков.
При расчете сверхзвуковых течений особую важность имеет способ аппроксимации невязких потоков Г. Схема аппроксимации должна обеспечивать возможность достаточно точного разрешения газодинамических разрывов на небольшом числе внутренних точек при отсутствии осцилля-ций поля течения в окрестности разрывов. Указанными свойствами обладают методы, получившие большую популярность за последние десятилетия, которые используют в том или ином виде характеристические свойства системы уравнений [22]. К ним
относятся схемы расщепления векторов потоков (например, схема расщепления Стегера — Уорминга, метод Ли8М) и схемы, основанные на решении задачи Ри-мана о распаде разрыва (например, схемы Годунова, Роу, ИЬЬ и ИЬЬС). Далее будут более подробно рассмотрены только схемы Ли8М и ИЬЬ, выбранные для проведения настоящих расчетов.
При вычислении потоков широко используется понятие параметров «слева» и «справа» от рассматриваемой грани (далее им присваиваются нижние индексы Ь и К). Если в качестве таких значений используются значения из центров прилегающих ячеек, то численный метод имеет первый порядок точности. Для реализации схем повышенного порядка точности применяются специальные методы определения параметров слева и справа от грани, кратко рассмотренные ниже.
Схема AUSM [21]. Основана на представлении вектора потока Г в виде суммы Г(с) (конвективная составляющая) и Гр (составляющая, связанная с давлением):
pV - п р V - пи + рп -1 рV - пу + рп - ] р V - п^ + рп - к рV - пН
(5)
" рV - п " 0
рV - пи рп -1
рV - пу + рп - ]
рV - п^ рп - к
рV - пН 0
= р(с) + г (р). Соответственно, Г,, (численный поток на
грани) тоже определяется как сумма:
Г, = Г}с) + Г,р).
(6)
При этом для каждой из составляющих используется свой способ расщепления.
Конвективный поток Г(с) выражается через число Маха М, построенное по нормальной составляющей скорости, —
М = V•n/a,
где a -^JyRT — скорость звука (R — газовая постоянная):
F(c) - M
= mF (с).
(7)
pa pau pav paw paH
Расщепление конвективного потока осуществляется в зависимости от знака числа Маха на грани М^ следующим образом:
F/c) -<
M / FF(c) ? L M/ > 0,
M / FF(c) ? R M/ < 0.
(8)
M L, R -
Часть потока, связанная с давлением, определяется как
0
Ff -(p+ + PR)
n • 1 n •j
n • k
0
(10)
Ph R - \
Ml,r| < i;
PL/r (ML/r ± |ML,R |)
(11)
M
L / R
|ml/r| > i.
Схема HLL [20]. Основана на приближенном решении задачи Римана о распаде произвольного разрыва. Предполагается, что решение состоит из двух основных разрывов, которые описывают распространение сильных особенностей типа ударных волн; другие особенности — типа контактных или тангенциальных разрывов — при этом не учитываются.
Скорости основных разрывов (характеристик) SL и SR, из которых состоит решение, определяются следующими выражениями [23]: _
SL = min (VL • n - aL, V • n - ä), (12)
SR — mm
in (VR • n + aR ,V • n + a), (13)
Число Маха на грани определяется как сумма положительной и отрицательной составляющих:
М 7 = М+ + Мд,
где расщепление на составляющие проводится с помощью следующего соотношения:
± 4 (М, / * ± 1)2 ,|М „ < 1;
. 4 (9)
2 (М Ь1 к ± |М Ь1 к\ ),|М £ / д| >1.
где аь, ак — скорости звука, посчитанные по параметрам слева и справа от грани; величины а и V, помеченные волной, — это переменные, которые для текущей грани рассчитываются посредством осреднения по Роу [24].
Согласно приближенному решению задачи Римана [20], характеристики SL и отделяют друг от друга три области с постоянными параметрами газа на диаграмме х - две области с параметрами невозмущенного газа «слева» и «справа» от грани, а третья — область между ними.
Численный поток на грани зависит от конфигурации, отвечающей текущей грани:
F/ —
Fl, 0 <Sl; F*, Sl < 0 < Sr; Fr, 0 > Sr,
(14)
где поток F* определяется формулой Sr Fl - Sl Fr + SlSr (wr - w l )
F*
SR SL
(15)
где расщепление давления на положительную и отрицательную составляющие также проводится в зависимости от числа Маха:
± р^ (МL/л ± 1)2 (2 + МL/л),
Повышение порядка точности. Повышенный порядок точности численной схемы достигается за счет использования техники квазимонотонной интерполяции сеточных решений — подхода MUSCL (Monotonic Upstream-Centered Scheme for Conservation Laws) [25], согласно которому проводится кусочно-полиномиальная реконструкция решения в каждом контрольном объеме и восстановление значений на грани с повышенным порядком. При этом для монотонизации решения в расчетах по схемам вто-
рого порядка точности применяются схемы уменьшения полной вариации (TVD-схемы, англ. Total Variation Diminishing) [26]. Обобщение данного подхода на случай неструктурированных сеток заключается в применении квазиодномерных вычислений при локальном выделении для каждой грани некоего подходящего направления (аналогичного координатному направлению, исходно присутствующему в структурированных сетках). Другими словами, помимо значений переменных в центрах ячеек, расположенных по обе стороны грани, требуется, как минимум, еще две виртуальные точки слева и справа от центральных точек (точки PL, P[ и P , PR на рис. 2), в которых каким-либо образом восстанавливаются значения переменных.
Представленные ниже численные решения получены с использованием подхода, предложенного в работе [27] и подробно описанного в работе [28]. При таком подходе для каждой грани f проводят две прямые (лучи реконструкции), которые проходят через центр грани и центры прилегающих к грани ячеек (прямые реконструкции — I и r на рис. 2); на каждой из этих прямых определяют по две дополнительные точки, значения переменной в которых находят посредством интерполяции по известным значениям переменных в некоторых близ-
лежащих центрах ячеек («интерполирующих тройках»).
Алгоритм нахождения таких ячеек [28] заключается в следующем: вначале определяется множество N1, включающее всех соседей первого уровня для текущей ячейки, т. е. множество ячеек, имеющих хотя бы один общий узел с текущей ячейкой, за вычетом самой ячейки. Далее все ячейки из множества Nl сортируются по возрастанию косинуса угла между лучом реконструкции, выпущенным из центра текущей ячейки, и вектором направления на центр ячейки из этого множества. Далее тройки ячеек перебирают в порядке возрастания суммы индексов этих ячеек в отсортированном массиве и используют первую из найденных «интерполирующих троек». Имеется в виду такая тройка ячеек, центры которых образуют треугольник, а прямая реконструкции его пересекает.
После определения значений в дополнительных точках, значения слева и справа на грани вычисляют в соответствии с линейной односторонней экстраполяцией:
и[ = щ + у (щ-щ)!2, (16)
и* = и(и+-и* )/2, (17)
Рис. 2. Схема к построению алгоритма квазиодномерных вычислений на
примере двумерного поля; темно-серым тоном выделены левая и правая ячейки, прилегающие к рассматриваемой грани / светло-серым — все ячейки шаблона, задействованные при определении значений в дополнительных точках (приведены их обозначения); г, I — прямые реконструкции для грани /
где u — любая из реконструируемых переменных; y(r) — ограничитель, введенный для контроля осцилляций и вычисляемый как функция отношения двух разностей:
Wl = WL ((ul- UL )/(UL - UL-)) > (18)
WR = WR (( UR- UR )/(UR - UR)) . (19)
В качестве функции у в настоящих расчетах использовался TVD-ограничитель Ван Альбады [29].
Постановка задачи и вычислительный инструментарий
На рис. 3 показана расчетная область для рассматриваемой задачи сверхзвукового обтекания удлиненного цилиндрического тела, установленного на адиабатической пластине, вдоль которой развивается ламинарный пограничный слой. Течение предполагается симметричным, и расчетная область охватывает только половину исходной конфигурации. Размеры области: R = 15D, h = 10D, L = 8D, где D - диаметр затупленной части тела, также полагаемого адиабатическим.
Задача определяется следующим набором безразмерных параметров: числом Маха набегающего потока M, числом Рейнольдса ReD, числом Прандтля Pr, показателем адиабаты у и отношением диаметра тела к толщине набегающего пограничного слоя D/S.
В данной работе численные решения получены при М = 2,95, ReD = 4000, Рг = 0,71, у = 1,4, ^/695% = 1. На входной границе расчетной области предписывались профили скорости и температуры для пограничного слоя заданной толщины 695%.
Описанный выше численный метод, в сочетании с неявной схемой «в приращениях», реализован авторами настоящей статьи как одна из опций конечно-объемного неструктурированного программного кода БШР/Р^-Б, разработанного на кафедре «Гидроаэродинамика, горение и теплообмен» Санкт-Петербургского политехнического университета Петра Великого.
Для проведения расчетов использовались вычислительные ресурсы суперкомпьютерного центра того же университета (www.scc.spbstu.ru).
Результаты расчетов и их обсуждение
Структура потока. Численное решение, полученное по схеме ЛИ8М на наиболее подробной из использованных сеток, представлено на рис. 4 — 6. (вопросы сеточной сходимости решения обсуждаются ниже). В целом, структура рассчитанного поля течения аналогична описанной ранее в работах, проведенных для более высоких значений чисел Маха и Рейнольдса [15, 16]. Перед телом образуется головной скачок уплотнения, который, взаимодействуя с погра-
Рис. 3. Расчетная область для задачи сверхзвукового обтекания удлиненного цилиндрического тела (см. рис. 1). Представлены геометрические параметры: Б — диаметр затупленной части тела; Я, Ь, к — размеры расчетной области; стрелкой указано направление потока
ничным слоем, вызывает отрыв последнего. Зона отрыва индуцирует косые волны сжатия; они пересекаются с головным скачком. Внутри отрывной области возникают зоны сверхзвуковых скоростей и местные волны уплотнения, приводящие к повторному отрыву пристенного течения. В результате перед телом формируется протяженная отрывная область с цепочкой вихревых образований, каждое из которых становится «головой» огибающего тело подковообразного вихря.
Поверхностные линии тока, приведенные на рис. 4, показывают области отрыва и присоединения пограничного слоя. На том же рисунке представлено распределение давления (оно отнесено к значению давления в набегающем потоке Рп). Видно, в частности, что максимальное давление в лобовой части обтекаемого тела превышает входное давление примерно в десять раз.
Более детальная визуализация течения в плоскости симметрии дается на рис. 5, где посредством изображения поля модуля градиента плотности иллюстрируется ударно-волновая структура в сочетании с картиной линий тока, свидетельствующей о формировании цепочки вихрей в отрывной области. Анализ рисунка показывает, что вихри, заполняющие отрывную область, индуцируют возникновение косых волн уплотнения, взаимодействующих с головным скачком. Это взаимодействие приво-
дит к изгибу фронта головного скачка в сторону обтекаемого тела, а также к формированию струи газа, которая при нате-кании на тело приводит к появлению зоны локального повышения давления; данная особенность подробнее обсуждается ниже.
На рис. 6 приведено распределение числа Маха в плоскости симметрии, пунктирной линией обозначена звуковая линия (M = 1). В отрывной области, где поток в основном дозвуковой, наблюдаются, однако, две зоны сверхзвукового течения. Поток, который из заторможенной области движется вдоль поверхности обтекаемого тела по направлению к пластине, разгоняется до сверхзвуковых скоростей, а затем сворачивается в вихрь, где также достигается сверхзвуковая скорость.
Сеточная сходимость. Исследование сеточной сходимости проводилось с использованием ряда квазиструктурированных сеток: сетка 1 содержала 0,3 млн. ячеек, сетка 2 — 2,4 млн., сетка 3 — 8,1 млн., сетка 4 — 13,3 млн. Сетка 2 получена путем измельчения сетки 1 в два раза по каждому координатному направлению, сетка 3 — измельчением сетки 2 в 1,5 раза, построение сетки 4 осуществлялось при помощи измельчения сетки 3 только вблизи передней кромки тела (по координате Z/D — от 0 до 2,5 и по координате X/D — от —2 до 0). Все сетки имели одинаковую структуру: со сгущением сеточных линий к обтекаемо-
Рис. 4. Картина поверхностных линий тока и структура течения в продольном (плоскость Х2) и поперечном (ХГ) сечениях потока. Показаны также распределения давления в этих сечениях
и на поверхности обтекаемого тела
Z/DY
IVPl/Pin
+03
7.3x1 0+ 3.9x10 2.1 х1 О+03 1.1x10+03 5.9x1 0+ог 2.8x1 0+ог 1.7х10+ог 9.0x1 0+01 4.8x1 0+°1 2.6x1 0+01 1.8x10+°1 1.4x1 О+01
-4 -3 -2 -1 Х/О
Рис. 5. Поле градиента плотности и линии тока в плоскости симметрии. Приводятся значения градиента плотности, рассчитанные посредством дифференцирования по безразмерным координатам и отнесенные к значению плотности в набегающем потоке р.
-4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 Х/О
Рис. 6. Поле числа Маха в плоскости симметрии; штриховой линией обозначена звуковая
линия М = 1
му телу и к поверхности пластины. Общий вид расчетной сетки 1 иллюстрируется на рис. 7.
Результаты, полученные по схеме Ли8М при использовании разных сеток, представлены на рис. 8, где показано распределение давления вдоль лобовой линии на поверхности тела, а также распределение коэффициента трения на пластине вдоль линии симметрии (давление отнесено к аналитически вычисляемому значению полного давления за прямым скачком уплотнения Ра). Отметим характерное повышение давления в области при ~ 1,5, которое связано с формированием струи газа при взаимодействии косых волн уплотнения с головным скачком. Примечательно также, что схема первого порядка точности не
воспроизводит данный характерный пик давления даже на самой подробной сетке; это обусловлено, прежде всего, недостаточностью разрешения косых волн уплотнения. Более того, по схеме первого порядка предсказывается существенно более простая вихревая структура отрывной области перед телом.
Решение, полученное на сетке 3, весьма близко к решению на сетке 4 как в отношении рассматриваемого распределения давления, так и трения на пластине. Это позволяет заключить, что при использовании схем второго порядка точности сетка 4 является достаточно подробной для разрешения всех деталей вихревой структуры в области сопряжения тела и пластины и обеспечивает также высококачественное разрешение ударно-волновой структуры.
Рис. 7. Расчетная сетка 1, содержащая 0,3 млн. ячеек (см. пояснения в тексте)
Рис. 8. Распределения безразмерного давления вдоль лобовой линии (а) и коэффициента
трения на пластине вдоль линии симметрии (Ь). Даны результаты расчетов по схеме АИ8М второго порядка точности на сетках 1 — 4 (номера кривых совпадают с номерами сеток) и первого порядка точности на сетке 4 (кривая 5)
Сопоставление решений, полученных по разным схемам
Сопоставительные расчеты проводились на самой мелкой сетке (сетка 4) с использованием схем ИЬЬ и Ли8М первого и второго порядков точности. Сравнение результатов, полученных при использовании различных схем, удобно проводить по-
средством анализа расчетных распределений давления вдоль выбранных линий на обтекаемом теле, а также предсказываемых распределений коэффициента трения вдоль отдельных линий на поверхности пластины. Выбранные для анализа линии показаны на рис. 9, а сравниваемые распределения — на рис. 10 и 11.
ю 8 6
40 2 О
-10
15
Рис. 9. Положение линий тока на обтекаемом теле, выбранных для анализа расчетных распределений давления и коэффициента трения (см. сравнительный анализ на рис. 10, 11)
Рис. 10. Распределения безразмерного давления вдоль вертикальных линий на поверхности обтекаемого тела (см. рис. 9): а — на лобовой линии (ф = 0°),
Ь — в конце закругления (ф = 90°), с — ниже по потоку (Х/Б = 4,5). Расчеты выполнены по разным схемам: 1, 2 — ИЬЬ первого и второго порядков точности соответственно; 3, 4 — ЛИ8М первого и второго порядков точности
Рис. 11. Распределения коэффициента трения вдоль нескольких линий на поверхности пластины: а — на линии симметрии (Г/Б = 0); Ь, с — на расстояниях двух (Г/Б = 2)
и четырех (У/Б = 4) диаметров от тела, соответственно. Расчеты выполнены по разным схемам: 1, 2 — ИЬЬ первого и второго порядков точности соответственно; 3, 4 — АИ8М первого и второго порядков точности
Вследствие повышенной диссипатив-ности схем первого порядка точности, решения, полученные по схемам первого и второго порядков, в отдельных частях потока различаются весьма существенно как в случае Аи8М, так и в случае ИЬЬ-схемы. Следует отметить, однако, что при использовании схемы Аи8М первого порядка ударно-волновая структура течения на наиболее подробной сетке воспроизводится в целом удовлетворительно. Рис. 10 показывает, в частности, что распределение давления по поверхности тела, полученное по схеме Аи8М первого порядка, близко к рассчитанному по схеме второго поряд-
ка (за исключением отдельных локальных особенностей); тогда как применение схемы ИЬЬ первого порядка приводит к распределению, которое имеет иной вид.
Что же касается качества разрешения вязких эффектов, то, как уже отмечалось, схема Аи8М первого порядка оказывается излишне диссипативной: даже в случае наиболее подробной сетки ее использование приводит к уменьшению числа разрешаемых вихрей и спаду их интенсивности (см. рис. 11). Применение схемы ИЬЬ первого порядка приводит к еще менее приемлемым результатам: предсказывается совершенно другой уровень трения на большей
части пластины (см. рис. 11). При общей близости результатов, полученных по схемам Ли8М и ИЬЬ второго порядка точности, имеется ряд локальных различий. В частности, из анализа распределения давления вдоль лобовой линии (см. рис. 10) следует, что в решении, полученном по схеме Ли8М, локальный пик давления имеет заметно меньшую ширину. Анализ же распределения трения на пластине (см. рис. 11) показывает, что интенсивность вихревого движения в решении, полученном по ИЬЬ-схеме, несколько ниже, чем в случае Ли8М-схемы.
Заключение
Получено семейство численных решений для модельной задачи о взаимодействии сверхзвукового течения вязкого газа (воздуха) с удлиненным затупленным телом, установленном на пластине, вдоль которой развивается ламинарный пограничный слой. Решения с использованием двух схем для расчета конвективных потоков (ИЬЬ и Ли8М) и на сетках разной размерности получены при числе Маха набегающего потока, равном 2,95, и числе Рейнольдса, равном 4000.
Течение, которое развивается при обтекании рассмотренной конфигурации, является трехмерным, с яркими проявления-
ми эффектов вязко-невязкого взаимодействия. Возникающая перед телом отрывная область имеет сложную вихревую структуру, а вокруг тела формируется семейство подковообразных вихрей, стелющихся вдоль пластины. Обтекание отрывной области сверхзвуковым потоком индуцирует возникновение косых волн уплотнения, взаимодействие которых с головным скачком приводит к формированию струи газа, вызывающей местное повышение давления на теле.
По результатам проведенного исследования сеточной сходимости показано, что в случае использования схем второго порядка точности, сетки, содержащие 13 — 15 млн. гексагональных ячеек, позволяют разрешить все детали вихревой структуры течения в области сопряжения тела и пластины, а также достичь высококачественного разрешения ударно-волновой структуры. Схемы первого порядка точности даже на самой подробной сетке из числа использованных не позволяют воспроизвести существенные особенности течения, обусловленные вязкими эффектами. Решения, полученные по схемам Ли8М и ИЬЬ второго порядка точности, находятся в хорошем согласии, но в целом схема ИЬЬ проявляет себя как более диссипативная.
СПИСОК ЛИТЕРАТУРЫ
1. Dolling D.S. Fifty 2019s of shock-wave/ boundary layer research: What next? // AIAA Journal. 2001. Vol. 39. No. 8. Pp. 1517-1531.
2. Babinsky H., Harvey J.K. Shock wave/ boundary-layer interactions. Cambridge: Cambridge University Press, 2011. 461 p.
3. Devenport W.J., Simpson R.L. Time-dependent and time-averaged turbulence structure near the nose of a wing-body juncture // Journal of Fluid Mechanics. 1990. Vol. 210. Pp. 23-55.
4. Fleming J.L., Simpson R.L., Cowling J.E., Devenport W.J. An experimental study of a turbulent wing-body junction and wake flow // Experiments in Fluids. 1993. Vol. 14. No. 5. Pp. 366-378.
5. Ballio F., Franzetti S. Topological analysis of a junction vortex flow // Proceedings of Advances in Fluid Mechanics. Montreal, Canada. 2000. 24-26 May. WIT Press, Southampton. 2000. Pp. 255-264.
6. Levchenya A.M., Smirnov E.M., Goryachev V.D. RANS-based numerical simulation and vi-
sualization of the horseshoe vortex system in the leading edge endwall region of a symmetric body // International Journal of Heat and Fluid Flow. 2010. Vol. 31. No. 6. Pp. 1107-1112.
7. Войтенко Д.М., Зубков А.И., Панов Ю.А. Обтекание цилиндрического препятствия на пластине сверхзвуковым потоком газа // Известия АН СССР. Механика жидкости и газа. 1966. № 1. С. 121.
8. Адуевский B.C., Медведев К.И. Физические особенности течения в области отрыва при трехмерном взаимодействии пограничного слоя с ударной волной // Известия АН СССР. Механика жидкости и газа. 1967. № 1. С. 25-34.
9. Тетерин М.П. Исследование течения газа в области падения скачка уплотнения на цилиндр, обтекаемый потоком большой сверхзвуковой скорости // Известия АН СССР. Механика жидкости и газа. 1967. № 2. С. 143 -147.
10. Dolling D.S., Bogdonoff S.M. Blunt fin-
induced shock wave/turbulent boundary-layer interaction // AIAA Journal. 1982. Vol. 20. No. 12. Pp. 1674-1680.
11. Ozcan O., Holt M. Supersonic separated flow past a cylindrical obstacle on a flat plate // AIAA Journal. 1984. Vol. 22. No. 5. Pp. 611-617.
12. Houwing A.F.P., Smith D.R., Fox J.S., Danehy D.M., Mudford N.R. Laminar boundary layer separation at a fin-body junction in a hypersonic flow // Shock Waves. 2001. Vol. 11. No. 1. Pp. 31-42.
13. Borovoy V., Mosharov V., Radchenko V., Skuratov A. The shock-waves interference in the flow around a cylinder mounted on a blunted plate //7th European Conference for Aeronautics and Space Sciences (EUCASS). 2017. Report 2017-63. Pp. 1-8.
14. Hung C.-M., Buning P.G. Simulation of blunt-fin-induced shock-wave and turbulent boundary-layer interaction // Journal of Fluid Mechanics. 1985. Vol. 154. Pp. 163-185.
15. Lakshmanan B., Tiwari S.N. Investigation of three-dimensional separation at wing/body junctions in supersonic flows // Journal of Aircraft. 1994. Vol. 31. No. 1. Pp. 64-71.
16. Tutty O.R., Roberts G.T., Schuricht P.H. High-speed laminar flow past a fin-body junction // Journal of Fluid Mechanics. 2013. Vol. 737. Pp. 19-55.
17. Peery K.M., Imlay S.T. Blunt body flow simulations // AIAA. 1988. July. Paper 88-2924.
18. Quirk J.J. A contribution to the great Riemann solver debate // International Journal for Numerical Methods in Fluids. 1994. Vol. 18. No. 6. Pp. 555-574.
19. Rodionov V. Artificial viscosity to cure the carbuncle phenomenon: The three-dimensional case // Journal of Computational Physics. 2018. Vol. 361. Pp. 50-55.
20. Harten A., Lax P.D., van Leer B. On
upstream differencing and Godunov-type schemes for hyperbolic conservation laws // SIAM Review. 1983. Vol. 25. No. 1. Pp. 35-61.
21. Liou M.S., Steffen C.J. A new flux splitting scheme // Journal of Computational Physics. 1993. Vol. 107. No. 1. Pp. 23-39.
22. Toro E.F. Riemann solvers and numerical methods for fluid dynamics. Berlin, Heidelberg: Springer-Verlag, 2009. 724 p.
23. Einfeldt B. On Godunov-type methods for gas dynamics // SIAM Journal on Numerical Analysis. 1988. Vol. 25. No. 2. Pp. 294-318.
24. Roe P.L. Approximate Riemann solvers, parameter vectors, and difference schemes // Journal of Computational Physics. 1981. Vol. 43. No. 2. Pp. 357-372.
25. Van Leer B. Towards the ultimate conservative difference scheme. V. A second order sequel to Godunov's method // Journal of Computational Physics. 1979. Vol. 32. No. 1. Pp. 101-136.
26. Harten A. High resolution schemes for hyperbolic conservation laws // Journal of Computational Physics. 1983. Vol. 49. No. 3. Pp. 357-393.
27. Le Touze С., Murrone A., Guillard H. Mul-tislope MUSCL method for general unstructured meshes // Journal of Computational Physics. 2015. Vol. 284. Pp. 389-418.
28. Бахвалов П.А., Козубская Т.К. Схема с квазиодномерной реконструкцией переменных, определенных в центрах элементов трехмерной неструктурированной сетки // Математическое моделирование. 2016. Т. 28. № 3. C. 79-95.
29. Van Albada G.D., Van Leer B., Roberts
W.W. A comparative study of computational methods in cosmic gas dynamics // Astronomy and Astrophysics. 1982. Vol. 108. No. 1. Pp. 76-84.
Статья поступила в редакцию 26.03.2019, принята к публикации 03.04.2019.
СВЕДЕНИЯ ОБ АВТОРАХ
КОЛЕСНИК Елизавета Владимировна — инженер кафедры «Гидроаэродинамика, горение и теплообмен» Санкт-Петербургского политехнического университета Петра Великого, Санкт-Петербург, Российская Федерация.
195251, Российская Федерация, г. Санкт-Петербург, Политехническая ул., 29 kolesnik_ev@mail.ru
СМИРНОВ Евгений Михайлович — доктор физико-математических наук, заведующий кафедрой «Гидроаэродинамика, горение и теплообмен» Санкт-Петербургского политехнического университета Петра Великого, Санкт-Петербург, Российская Федерация.
195251, Российская Федерация, г. Санкт-Петербург, Политехническая ул., 29 smirnov_em@spbstu.ru
СМИРНОВСКИЙ Александр Андреевич — кандидат физико-математических наук, доцент кафедры «Гидроаэродинамика, горение и теплообмен» Санкт-Петербургского политехнического университета Петра Великого, Санкт-Петербург, Российская Федерация.
195251, Российская Федерация, г. Санкт-Петербург, Политехническая ул., 29 smirta@mail.ru
REFERENCES
1. Dolling D.S., Fifty years of shock-wave/ boundary layer research: What next?, AIAA Journal. 39 (8) (2001) 1517-1531.
2. Babinsky H., Harvey J.K., Shock wave/ boundary-layer interactions, Cambridge University Press, Cambridge (2011).
3. Devenport W.J., Simpson R.L., Time-dependent and time-averaged turbulence structure near the nose of a wing-body juncture, Journal of Fluid Mechanics. 210 (1990) 23-55.
4. Fleming J.L., Simpson R.L., Cowling J.E., Devenport W.J., An experimental study of a turbulent wing-body junction and wake flow, Experiments in Fluids. 14 (5) (1993) 366-378.
5. Ballio F., Franzetti S., Topological analysis of a junction vortex flow, Proceedings of Advances in Fluid Mechanics, Montreal, Canada, 24-26 May, 2000. WIT Press, Southampton (2000) 255-264.
6. Levchenya A.M., Smirnov E.M., Goryachev V.D., RANS-based numerical simulation and visualization of the horseshoe vortex system in the leading edge endwall region of a symmetric body, International Journal of Heat and Fluid Flow. 31 (6) (2010) 1107-1112.
7. Voytenko D.M., Zubkov A.I., Panov Yu.A., Obtekaniye tsilindricheskogo prepyatstviya na plastine sverkhzvukovym potokom gaza [Flow past a cylindrical obstacle on the plate by a supersonic gas stream], Fluid Dynamics. (1) (1966) 121.
8. Aduyevskiy B.C., Medvedev K.I., Fizicheskiye osobennosti techeniya v oblasti otryva pri trekhmernom vzaimodeystvii pogranichnogo sloya s udarnoy volnoy [Physical features of the flow in the separation region in three-dimensional interaction of the boundary layer with a shock wave], Fluid Dynamics. (1) (1967) 25 -34.
9. Teterin M.P., Issledovaniye techeniya gaza v oblasti padeniya skachka uplotneniya na tsilindr, obtekayemyy potokom bolshoy sverkhzvukovoy skorosti [Investigation of gas flow in the region of the shock drop sealing on a cylinder, streamlined by a large supersonic flow], Fluid Dynamics. (2) (1967) 143-147.
10. Dolling D.S., Bogdonoff S.M., Blunt fin-induced shock wave/turbulent boundary-layer interaction, AIAA Journal. 20 (12) (1982) 1674-
1680.
11. Ozcan O., Holt M., Supersonic separated flow past a cylindrical obstacle on a flat plate, AIAA Journal. 22 (5) (1984) 611-617.
12. Houwing A.F.P., Smith D.R., Fox J.S., et al., Laminar boundary layer separation at a fin-body junction in a hypersonic flow, Shock Waves. 11 (1) (2001) 31-42.
13. Borovoy V., Mosharov V., Radchenko V., Skuratov A., The shock-waves interference in the flow around a cylinder mounted on a blunted plate, 7th European Conference for Aeronautics and Space Sciences (EUCASS). Report 2017-63.
(2017) 1-8.
14. Hung C.-M., Buning P.G., Simulation of blunt-fin-induced shock-wave and turbulent boundary-layer interaction, Journal of Fluid Mechanics. 154 (1985) 163-185.
15. Lakshmanan B., Tiwari S.N., Investigation of three-dimensional separation at wing/body junctions in supersonic flows, Journal of Aircraft. 31 (1) (1994) 64-71.
16. Tutty O.R., Roberts G.T., Schuricht P.H., High-speed laminar flow past a fin-body junction, Journal of Fluid Mechanics. 737 (2013) 19-55.
17. Peery K.M., Imlay S.T., Blunt body flow simulations, AIAA, Paper 88-2924 (1988).
18. Quirk J.J., A contribution to the great Riemann solver debate, International Journal for Numerical Methods in Fluids. 18 (6) (1994) 555-574.
19. Rodionov V., Artificial viscosity to cure the carbuncle phenomenon: The three-dimensional case, Journal of Computational Physics. 361
(2018) 50-55.
20. Harten A., Lax P.D., van Leer B., On upstream differencing and Godunov-type schemes for hyperbolic conservation laws, SIAM Review. 25 (1) (1983) 35-61.
21. Liou M.S., Steffen C.J., A new flux splitting scheme, Journal of Computational Physics. 107 (1) (1993) 23-39.
22. Toro E.F., Riemann solvers and numerical methods for fluid dynamics, Springer-Verlag, Berlin, Heidelberg (2009).
23. Einfeldt B., On Godunov-type methods for gas dynamics, SIAM Journal on Numerical
Analysis. 25 (2) (1988) 294-318.
24. Roe P.L., Approximate Riemann solvers, parameter vectors, and difference schemes, Journal of Computational Physics. 43 (2) (1981) 357-372.
25. Van Leer B., Towards the ultimate conservative difference scheme. V. A second order sequel to Godunov's method, Journal of Computational Physics. 32 (1) (1979) 101-136.
26. Harten A., High resolution schemes for hyperbolic conservation laws, Journal of Computational Physics. 49 (3) (1983) 357-393.
27. Le Touze S., Murrone A., Guillard H., Multislope MUSCL method for general
Received 26.03.2019, accepted 03.04.2019.
unstructured meshes, Journal of Computational Physics. 284 (2015) 389-418.
28. Bakhvalov P.A., Kozubskaya T.K., Skhema s kvaziodnomernoy rekonstruktsiyey peremennykh, opredelennykh v tsentrakh elementov trekhmernoy nestrukturirovannoy setki [Cell-centered quasi-one-dimensional reconstruction scheme on 3D hybrid meshes], Matem. modelirovaniye. 28 (3) (2016)79 -95.
29. Van Albada G.D., Van Leer B., Roberts W.W., A comparative study of computational methods in cosmic gas dynamics, Astronomy and Astrophysics. 108 (1) (1982) 76-84.
THE AUTHORS
KOLESNIK Elizaveta V.
Peter the Great St. Petersburg Polytechnic University
29 Politechnicheskaya St., St. Petersburg, 195251, Russian Federation
kolesnik_ev@mail.ru
SMIRNOV Evgueni M.
Peter the Great St. Petersburg Polytechnic University
29 Politechnicheskaya St., St. Petersburg, 195251, Russian Federation
smirnov_em@spbstu.ru
SMIRNOVSKY Alexander A.
Peter the Great St. Petersburg Polytechnic University
29 Politechnicheskaya St., St. Petersburg, 195251, Russian Federation
smirta@mail.ru
© Санкт-Петербургский политехнический университет Петра Великого, 2019
