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

Возможности гибридного метода аппроксимации конвективных потоков при моделировании течений сжимаемых

сред

М.В. Крапошин, <[email protected]> Институт системного программирования РАН, Россия, Москва, ул. Солженицына, д. 25

Аннотация. Для моделирования течений в широком диапазоне чисел Маха предложен гибридный метод аппроксимации конвективных слагаемых, основанный на схеме Курганова-Тадмора и разновидности метода проекций PISO (Pressure Implicit With Splitting Operators). Особенность данного метода состоит в неявной выражении конвективных потоков из схемы Курганова-Тадмора и введении специальной функции-переключателя, обеспечивающей в зависимости от локальных характеристик потока переход от «сжимаемой» схемы (Курганова-Тадмора) к «несжимаемой» схеме (стандартная аппроксимация, используемая в методе PISO). Использование такой гибридной схемы позволяет получить следующие преимущества: а) за счёт неявного учёта диффузионных слагаемых шаг по времени не ограничен скоростью распространения волн диффузионным механизмом; б) за счёт аппроксимации конвективных слагаемых неявным способом и перехода к стандартным схемам PISO шаг по времени ограничивается только потоковым числом Куранта; в) при необходимости разрешения акустических волн достаточно снижения шага по времени до достижения акустическим числом Куранта значений меньше 1 во всей области; г) использование схемы Курганова-Тадмора позволяет получить неосциллирующее решение в задачах с распространением акустических сигналов или при М > 0.3. В данной работе выполнено тестирование реализации гибридного метода для широкого класса задач с известным аналитическим решением и экспериментальными данными: а) сжимаемые течения — распространение волны в прямом канале (Задача Сода), обтекание плоского клина, обтекание обратного уступа сверхзвуковым потоком, обтекание прямого уступа сверхзвуковым потоком, течение в сверхзвуковом сопле при наличии прямого скачка уплотнения в закритической части; б) несжимаемые течения — дозвуковое течение ламинарного вязкого потока в канале круглого сечения, обтекание цилиндра в ламинарном режиме; обтекание цилиндра турбулентным потоком, течение струй газов со смешением; в) промышленные и академические верификационные задачи — истечение струи газа из сверхзвукового сопла, истечение квазиравновесной расширяющейся струи плазмы в вакуум; г) качественное исследование адекватности модели для задач промышленного масштаба — моделирование течения в высокоскоростном компрессоре, модель гидродинамики

водокольцевого насоса. Все материалы работы и исходный код свободно доступны через проект GiHub https:// github. com/unicfdlab.

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

DOI: 10.15514/ISPRAS-2016-28(3)-16

Для цитирования: Крапошин М.В. Возможности гибридного метода аппроксимации конвективных потоков при моделировании течений сжимаемых сред. Труды ИСП РАН, том 28, вып, 3, 2016, стр. 267-326. DOI: 10.15514/ISPRAS-2016-28(3)-16

1. Список сокращений и обозначений

NASA — National Aero Space Agency

PISO — Pressure Implicit With Splitting Operators

SIMPLE — Semi-Implicit Method for Pressure Linked Equations

U — Скорость

R — Индивидуальная газовая постоянная f p — удельная изобарная теплоёмкость Су — удельная изо хорная теплоёмкость 1 — показатель адиабаты Р — давление среды Т — температура среды Р — плотность среды

Уг — массовая доля компоненты или фазы в объёме среды М — число Маха а — локальная скорость звука D — коэффициент диффузии

2. Введение

Современный рост вклада численного моделирования в процессы проектирования и эксплуатации инженерных сооружений связан в первую очередь с ростом вычислительных мощностей. В то же время одним из основных факторов, тормозящих применение численных методов в промышленности, является жёсткая привязка последних к выбору математической модели (даже в рамках единого физического представления) и как результат — к конкретной прикладной области. Примером такого разделения численных методов может служить такой параметр сплошной среды как сжимаемость, разбивающая механику жидкости и газа на два больших класса — сжимаемые и несжимаемые течения. Изменение класса 268

решаемой прикладной задачи приводит к изменению выбора численного метода, программных средства его реализации и набора исходных данных, вплоть до смены технологического процесса решения. Численные методы, предназначенные для решения дозвуковых задач (в первую очередь в смысле несжимаемости) малопригодны для решения транс- и сверхзвуковых задач в силу их особенностей, обусловленных построением соответствующих численных схем. В первом случае чаще всего используются методы проекций (такие как PISO/SIMPLE), а во втором случае — методы на основе приближённого решения задачи Римана. Метод проекций, хорошо зарекомендовавший себя в дозвуковых течениях (M < 0.3, [27]), обладает высокой устойчивостью решения, обеспечивающей возможность интегрирования с большим шагом по времени. При этом переход в «сжимаемую область» сопровождается паразитными осцилляциями, которые ставят под сомнение целесообразность использования этого инструмента для случаев с M > 0.3 или для разрешения акустических волн. Ещё одним важным достоинством метода проекций является его расширяемость и относительная простота разработки сопряжённых моделей благодаря стандартной процедуре связывания плотности, скорости и давления [28].

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

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

Метод реализован в качестве самостоятельного приложения-«решателя» на основе открытой библиотеки OpenFOAM [2].

3. Краткое описание гибридного метода

Основная идея ранее предложенного гибридного метода [1] состоит во введении функции переключателя, которая в зависимости от состояния среды в двух соседних ячейках «смешивает» конвективные потоки, аппроксимированные с использованием метода Курганова-Тадмора, с потоками, вычисленными в соответствии со стандартным методом проекций, адаптированным к сжимаемым течениям. Процедура вычисления потоков гибридным методом включает в себя следующие этапы.

1. Выражение потоков (энергии, импульса, массовых долей) в неявном виде относительно интенсивных переменных в соответствии с процедурой KT/KNP (Kurganov-Tadmor/Kurganov-Noelle-Petrova) и в соответствии с методом проекций.

2. Введение функции-переключателя для «смешивания» конвективных потоков, вычисляемых двумя разными методами.

3. Формирование системы линейных алгебраических уравнений для каждого балансового соотношения и решение этой системы.

4. Формирование алгебраического уравнения для давления на основе дискретизированного уравнения неразрывности. Поиск нового поля давления, удовлетворяющего условию неразрывности, коррекция массовых потоков в соответствии с новым давлением среды.

4. Реализация гибридного метода

Предложенный гибридный метод и его модификации для случаев течения многокомпонентных и двухфазных гомогенных сред был реализован с использованием открытой библиотеки OpenFOAM в виде самостоятельных приложений-«решателей», см. рис. 1. Реализованные приложения находятся в открытом доступе и хранятся в архиве веб-сервиса GitHub по адресу: https ://github. co m/unicfdlab.

Многофазные

течения (multiphase)

reactingCentralFoam

Горение/Течение реагирующих сред

| (combustion) Теплофизические модели (thermophysicalM odels)

Движение модели расчетной ^ ^ g- 11 турбулентности области и .. JIT , ... . . . (dynamicFvMesh) (turbulentModels) Методы аппроксимации уравнение (finiteVolume)

Рис. 1. Приложения-решатели, реализующие гибридный метод и расширяющие стандартный набор моделей пакета OpenFOAM

Fig. 1. Solver applications implementing hybrid method and extending the standard model set

of the OpenFOAM package

Разработанные приложения расширяют существующий функционал стандартных моделей OpenFOAM и разделяются на три группы в соответствии с принятым в OpenFOAM способом классификации (рис. 1), полное описание возможностей приложений приведено в табл. 1.

• Сжимаемые и несжимаемые течения — решатели pisoCentralFoam/rhoPisoCentralFoam/pisoCentralDyMFoam.

• Сжимаемые течения реагирующих сред (совершенных газов) — решатель reactingCentralFoam.

• Гомогенная двухфазная модель течения двух сжимаемых сред — twoPhaseMixingCentralFoam/twoPhaseMixingCentralDyMFoam.

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

Табл. 1. Список приложений- «решателей», разработанных на основе библиотеки OpenFOAM и реализующих гибридный метод моделирования сжимаемых течений.

Table 1. List of applications-"solvers" developed on basis of the OpenFOAM library and implementing a hybrid method for modeling compressible flows.

№№ Наименование Описание

Сжимаемые и несжимаемые течения однофазные течения

1. pisoCentralFoam Модель ламинарного или турбулентного течения сжимаемой среды (совершенный газ) при числах Маха от 0 до 6 с возможностью переключения между стационарной или нестационарной численной схемой.

2. rhoPisoCentralFoam Модель ламинарного или турбулентного течения сжимаемой среды с уравнением состояния реального газа при числах Маха от 0 до 6 и с возможностью переключения между стационарной или нестационарной численной схемой.

3. pisoCentralDyMFoam Нестационарная модель ламинарного или турбулентного течения сжимаемой среды (совершенный газ) при числах Маха от 0 до 6 с возможностью моделирования случаев в условиях подвижной расчётной области.

Сжимаемые течения реагирующих сред (совершенных газов)

4. reactingCentralFoam Модель ламинарного или турбулентного течения сжимаемой многокомпонентной среды (совершенный газ) при числах Маха от 0 до 6 и учётом кинетики химических превращений составляющих потока.

Двухфазные течения

5. twoPhaseMixingCentralFoam Модель ламинарного или турбулентного течения гомогенной двухфазной сжимаемой смеси.

6. twoPhaseMixingCentralDyMFoam Модель ламинарного или турбулентного течения гомогенной двухфазной сжимаемой смеси моделирования случаев в условиях подвижной расчётной области.

5. Рассмотренные задачи

С целью тестирования гибридного метода были рассмотрены следующие группы задач.

1. Валидационные задачи для случая сжимаемого течения. В эту группу входят задачи с относительно простой геометрией области течения, имеющие либо аналитическое решение, либо эталонные данные из эксперимента, либо результаты численного моделирования с помощью других численных методов или пакетов. Всего в эту группу входят следующие задачи: а) распространение волны в прямом канале (Задача Сода); б) обтекание плоского клина; в) обтекание обратного уступа сверхзвуковым потоком; г) обтекание прямого уступа сверхзвуковым потоком; д) течение в сверхзвуковом сопле при наличии прямого скачка уплотнения в закритической части.

2. Валидационные задачи для случая несжимаемого течения: а) дозвуковое течение ламинарного вязкого потока в канале круглого сечения; б) обтекание цилиндра в ламинарном режиме; в) обтекание цилиндра турбулентным потоком; г) течение струй газов со смешением.

3. Промышленные верификационные задачи: а) истечение струи газа из сверхзвукового сопла; б) истечение квазиравновесной расширяющейся струи плазмы в вакуум.

4. Задачи промышленного масштаба: а) моделирование течения в высокоскоростном компрессоре; б) моделирование гидродинамики водокольцевого насоса.

6. Результаты тестирования

Проведённое тестирование охватило широкий круг задач — сжимаемые и несжимаемые течения, течения в подвижных областях, многокомпонентные и двухфазные течения.

6.1 Распространение волны в прямом канале (Задача Сода)

Рассматривается случай распространения ударной волны в цилиндрическом канале (ударной трубе). Волна создаётся расширением сжатого воздуха с высоким давлением и температурой, в область, заполненную газом с более низкими давлением и температурой. На рис. 2 представлена схема рассматриваемой задачи. Газ слева и справа от перегородки — воздух при разных температурах и давлениях. Задача является одномерной и имеет аналитическое решение [3]. В начальный момент времени области с разным давлением разделены диафрагмой, в момент разрыва диафрагмы начинается распространение волны сжатия в сторону газа низкого давления, волны 274

разрежения в сторону газа высокого давления, а также контактного разрыва. В зависимости от соотношения давлений в правой и левой частях ударной трубы перетекание газа может происходить со звуковой и дозвуковой скоростями. При тестировании решателя pisoCentralFoam моделировались оба этих случая (см. табл. 2).

Табл. 2.Начальныеусловиявзадачераспадаразрыва.

Table 2.The initial conditions for the Sod's problem.

Область Докритическоеистечение Критическоеистечение

Давление, Па Температура, К Давление, Па Температура, К

Слеваотдиафрагмы P4=6,897*104 T4=288,89 P4=6,897*104 T4=288,89

Справа от диафрагмы pi=5,897*104 T1=288,89 P1=6,897*103 ti=23 1,11

Diaphragm

Pi = 6.897 ■ 104 Pa © ® PI = 6.897 103 Pa

ТА = 288.89 К Driver Driven Ti = 231.11 К

0cm 15.42cm 30.48cm

Рис. 2. Схема расчетной области для случая распространения волны в канале

Fig. 2. Computational domain for the case of wave propagation in the channel (Sod's

problem)

Тестирование решателя pisoCentralFoam проводилось в одномерной (количество ячеек - 100), двумерной (1000 ячеек), осесимметричной двумерной (1000 ячеек) и полностью трёхмерной постановке (30000 ячеек) с целью определения влияния размерности задачи на результат. На рис. 3 показаны расчётные сетки для двумерного (2а) и трёхмерного (2б) случаев.

Рис. 3. Варианты расчетной сетки (а — двумерная, б — трехмерная) для задачи Сода

Fig. 3. Types of computational mesh (а — two-dimensional, and б — three-dimensional) used

for the Sod problem

Расчёты проводились до момента времени t=0.00025c. Как показали результаты, распределение давления по оси трубы не зависит от размерности задачи. На рис. 4 и 5 приводится сравнение результатов, полученных по одномерному расчёту (т. к. было установлено, что размерность задачи не влияет на результат) с аналитическим решением. Из представленных графиков видно отсутствие осцилляций и хорошая согласованность результатов, полученных с помощью разработанного метода, с аналитическим решением. Вторым важным результатом является совпадение решения, полученного с использованием гибридного метода и исходной схемы Курганова-Тадмора.

Рис. 4. Сравнение расчётного и аналитического распределения давления и скорости вдоль оси трубы для случая докритического течения. Обозначения: «analytical» -аналитическое решение, «pisoCentralFoam» — численное решение, полученное с помощью реализация гибридного метода в OpenFOAM.

Fig. 4. Comparison of calculated and analytical distributions ofpressure and velocity along

the tube axis in the case of subcritical flow. Legend: «analytical» - analytical solution, «pisoCentralFoam» - the numerical solution obtained with the implemented hybrid method

Рис. 5. Сравнение расчётного и аналитического распределения давления и скорости вдоль оси трубы для случая критического течения. Обозначения: «analytical» -аналитическое решение, «pisoCentralFoam» — численное решение полученное с помощью реализация гибридного метода в OpenFOAM.

Fig. 5. Comparison of calculated and analytical distributions ofpressure and velocity along the tube axis in the case of critical flow. Legend: «analytical» - analytical solution, «pisoCentralFoam» - numerical solution obtained with the implemented hybrid method

6.2 Обтекание плоского клина

Рассматривалось сверхзвуковое обтекание плоского клина (рис. 6). В качестве исходных данных использовались материалы, представленные в работе [4]. Число Маха набегающего потока М=2.5. Рабочая среда - сухой воздух, молярная масса - 28.96 г/моль, удельная газовая постоянная - 287.05 Дж/кг/К. Принималась гипотеза о возможности моделирования рабочей среды в качестве идеального газа. Давление и температура в набегающем потоке -101350 Па и 288.9 К.

Изобарная теплоёмкость принималась Cp=1004 Дж/кг/К, показатель адиабаты Y = 1.4. Скорость звука среды - a = 340.73 м/с. Скорость набегающего потока U = M • а = 851.84 м/с. Динамическая вязкость среды принималась равной 18.4 мкПа-с. Число Прандтля Pr = 1.

Как известно, для данной задачи существует приближенное аналитическое решение в рамках теории косых скачков уплотнений (см. [4]). Таким образом, с помощью данной задачи можно проверить возможности схемы по воспроизведению скачков уплотнений: их положения и «размытости», что позволяет в целом оценить близость численного решения к аналитическому. В качестве параметра, характеризующего близость численного решения к аналитическому, рассматривалось число Маха, скачкообразно меняющееся при прохождении через скачок уплотнения. С этой целью были расставлены точки отбора расчётных значений параметров потока, отстоящие на 0.05 м по координате Y от твёрдой стенки.

а)

б)-

Рис. 6. Расчетная схема для случая набегания потока на клин (а) и сетка в расчётной области для рассматриваемого случая (б)

Fig. 6. Geometry and settings for the case of flow attack on the inclined wedge (а) and the mesh in the computational domain for this case (б)

Для данного случая строилась двухблочная двумерная сетка (рис. 6б). Первый блок представлял собой прямоугольник размерами 0.1522x0.3048 м, второй

(над клином) - прямоугольную трапецию высотой 0.3048 м, равной нижнему основанию, и с боковой стороной, наклонённой к оси ОХ на 15°. Для проверки сходимости проводилось моделирование для трёх уровней сгущения сетки. Первоначальное разбиение блоков на ячейки: 75х50 ячеек (первый блок -25х50; второй блок - 50х50). Для получения более грубой и более точной сетки количество ячеек на каждый блок уменьшалось и увеличивалось в 1.5 раза от базового соответственно. Сеточная сходимость проиллюстрирована на рис. 8.

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

Рис. 7. Поле давления при обтекании косого уступа сверхзвуковым потоком Fig. 7. The field of pressure at the supersonic flow attack on the inclined wedge

2,55 2,45 2.35

а

£ 2,15 | 2,05 1,95 1,85 1,75

Чл

ч

\

V

к i 'Л

0,02

0,04

0,06

0,0В 0,1 X', m

-Normal-----------Ma finer- ■

0,12 0,14 0,16 Marougher-theoretical

0,18

0,2

Рис. 8. Сеточная сходимость и сравнение численного и точного решения для случая

набегания потока на клин

Fig. 8. Mesh convergence and comparison of numerical and exact solutions for the case of attack of stream inclined wedge

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

6.3 Обтекание обратного уступа сверхзвуковым потоком

В качестве третьего тестового примера рассматривалась классическая задача из теории отрывных течений - плоское сверхзвуковое обтекание обратного уступа. На рис. 9 представлена схематичная картина течения. Поток при прохождении кромки уступа расширяется, образуя веер волн разрежения. Наличие преграды в виде горизонтальной поверхности за уступом обуславливает отрыв вязкого потока. Присоединение потока ведёт к образованию Х-образного скачка уплотнения.

Рис. 9. Схема течения при сверхзвуковом обтекании обратного уступа

Fig. 9. Geometry and seetings for case of supersonic flow over backward step

В качестве исходных данных использовались материалы, представленные в [6]. Число Маха набегающего потока М=2.5. Рабочая среда - сухой воздух, молярная масса - 28.96 г/моль, удельная газовая постоянная - 287.05 Дж/кг/К. Принималась гипотеза о возможности моделирования рабочей среды в качестве совершенного газа. Использовалась стандартная k-ю SST модель турбулентности. Задание параметров турбулентности k, т на входе расчетной области производилась посредством известных зависимостей: k=3/2*(U I)2, где U - средняя скорость течения, а I - интенсивность турбулентности, принимаемая равной 5%; т =s/(k^Cp), где кинетическая энергия турбулентности е= C3/4k3/2/l, l=0,07L, где L- характерный размер, а коэффициент С принимается равным 0.09.

Статическое давление в набегающем потоке - 13316.6Па, давление торможения — 227527Па, температура набегающего потока — 153.04К, температура торможения - 344.44К.

Изобарная теплоёмкость принималась равной Cp=1005 Дж/кг/К, показатель адиабаты у=1.4. Скорость звука среды a = 248м/с. Скорость набегающего потока U = M * a = 620м/с. Динамическая вязкость среды принималась равной 18.27мкПа*с. Число Прандтля Pr =0.7.

Высота уступа 0.01125м, расстояние от уступа до входного сечения 0.1016м, до выходного сечения 0.3048м, расстояние до верхней границы расчётной области 0.1475м. Диапазон характерных для течения чисел Рейнольдса: 7*103 — 5*106 [5].

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

и насчитывал 240x40 ячеек, второй располагался над уступом (104х112 ячеек), третий располагался над первым блоком и замыкал расчётную область (240x112 ячеек). На всех твёрдых поверхностях задавалось граничное условие прилипания, k и ю аппроксимировались при помощи пристеночных функций. Результаты моделирования сравниваются с опытными данными по обтеканию обратного уступа той же геометрии, представленными в работе [5], а также с расчётными данными, полученными в кодах PARC,WIND и ANSYSFluidDynamics. Расчёт проводился до 0.06c, что соответствовало установившемуся режиму течения.

На рис. 10 представлены графики распределения давления за уступом в расчёте, в эксперименте и в расчётах других авторов. Давление отнесено к статическому давлению перед уступом, а горизонтальная координата отсчитывается от стенки уступа в дюймах.

Как видно из графика, модель, реализованная в pisoCentralFoam, приводит к результатам, наиболее близким к полученным с помощью ANSYS [6] и в эксперименте, указанном в [6]. В то же время положение кривой, полученной с помощью разработанной модели, несколько отличается от экспериментальных данных из источника [5] и результатов кодов WIND и PARC [7], а предсказываемое давление в отрывной зоне несколько завышено. Следует отметить, что распределение давления, которое приводится в руководстве пользователя ANSYS [6], в отчёте [5] найдено не было. Для наглядной оценки адекватности воспроизведения пространственных характеристик среды при сверхзвуковом обтекании со скачками уплотнения выполнено сравнение с имеющимися экспериментальными и расчётными данными. На рис. 11 представлена картина течения, полученная с помощью разработанной модели, реализованной pisoCentralFoam, и сравнение с экспериментальным положением скачка уплотнения [5] и расчётным положением скачка уплотнения, полученным в коде PARC [7]. Как видно из рисунка, положение скачка уплотнения лежит достаточно близко как к экспериментальным, так и к расчётным данным.

Результат проведённого исследования сеточной сходимости представлен на рис. 12. Схема имеет второй порядок точности по пространству и решение сходится к точному; при измельчении сетки в области отрыва по каждому направлению в 4 и более раз погрешность не превышает 1%.

0 0,5 1 1,5 2 2,5 3 3,5 4 4,5 5

x, inches

Рис. 10. Сравнение распределения давления за обратным уступом Fig. 10. Comparison of the pressure distribution behind backward step

Рис. 11. Сравнение положений скачка уплотнения, вызванного присоединением потока. Серыми точками показано экспериментальное положение скачка уплотнения, синий линией — результат расчёта кодом PARC[7]

Fig. 11. Comparison of the shock position caused by the flow reattachment. Grey dots show the experimental position of the shock wave, the blue line - the result of the calculation with

PARC code [7]

Mesh relative resolution

Рис. 12. Сеточная сходимость разработанной численной схемы на примере случая обтекания обратного уступа

Fig. 12. Mesh convergence of the developed numerical scheme by the example of the case of

a flow over backward step

6.4 Обтекание прямого уступа сверхзвуковым потоком

Рассматривается сверхзвуковое течение идеального газа в канале с резким сужением (рис. 13). Данная задача является классическим тестом методов моделирования сверхзвуковых течений (см. [8]). В начальный момент времени по всему пространству канала скорость, давление и температура распределены равномерно. Граничные условия следующие.

• На уступе (заштрихован) - условие непротекания для скорости, условие адиабатичности для температуры (нулевая нормальная производная), условие непроницаемости для давления (нулевая нормальная производная).

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

• На входе в расчётную область (левая вертикальная граница) -фиксированные значения скорости (3м/с), давления (1 Па), температуры (1 К).

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

Физические свойства газа подобраны так, чтобы в начальный момент времени число Маха потока в горизонтальном направлении равнялось 3, а показатель адиабаты (Ср/^) - 1.4.

• Молярная масса: 11640.3 г/моль.

• Адиабатная теплоёмкость 2.5 Дж/кг/К.

• Число Прандтля 1.

Параметры разбиения приведены в табл. 3. Расчёт проводился до момента времени t=4с. Сравнивалось положение «ножки» А-скачка и наличие/отсутствие неустойчивости Кельвина-Гельмгольца, которая наблюдалась во многих работах при воспроизведении данного случая. При корректной дискретизации в момент времени t=4с положение ножки А-скачка должно приходиться на передний край уступа (Х=60см).

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

зоо cm

60 cm

NXl ■*-----------------------— NX? * T

L X -*■- ;:>> В £! 1= IK. R z i p=lPa, У/ у i ! U=3m/s

О

Рис. 13. Схема расчётной области и разбиения на блоки.

Fig. 13. Geometry of computational domain and settings for block mesh for the case of supersonic flow over forward step.

Табл. 3. Параметры разбиения расчётной области

Tab. 3. The parameters of the block mesh of the computational domain

№ п.п. Отрезок Число отрезков разбиения

грубая сетка средняя сетка улучшенная сетка мелкая сетка

1 NX1 96 192 384 768

2 NX2 384 768 1536 3072

3 NY1 32 64 128 256

4 NY2 128 256 512 1024

Число ячеек в расчётной сетке

60 тыс 250 тыс 1 млн 4 млн

Из расчётов видно, что вне зависимости от степени измельчения расчётной сетки, положение «ножки» Л— скачка соответствует положению уступа (рис. 14). При этом в случае использования схемы первого порядка по времени, неустойчивость Кельвина-Гельмгольца начинает проявляться только на самой мелкой сетке (4 млн ячеек), в то время как использование схемы второго порядка аппроксимации по времени позволяет воспроизводить этот эффект и на заметно более грубой сетке (1 млн ячеек) — рис. 15. Исследование масштабируемости (рис. 16) решателя показало удовлетворительное ускорение на сетке 1 млн ячеек и сверхлинейное ускорение для сетки 4 млн. ячеек. Сверхлинейное ускорение связано с невозможностью размещения всех данных численной модели (4 млн ячеек) в кэше процессора в однопроцессорном режиме расчета, что привело к завышенному относительному ускорению.

Moderate mesh,

1st order

time scheme

Fine mesh, 1st order time scheme

i.500e+00

Very fine mesh,

1st order

time scheme

Рис. 14. Сравнение картин течения (поле плотности) для «средней» сетки, улучшенной сетки и мелкой сетки в момент времени 1=4с. 1-й порядок аппроксимации по времени и 2-й порядок аппроксимации по пространству.

Fig. 14. A comparison of flow visualisation (field density) for the moderate grid, fine grid and very fine grid at time t = 4s. First order approximation with respect to time and the second order approximation in space were used.

rho

1 Very fine mesh, 1st order time scheme

1 Fine mesh, 2nd order time scheme

i.500e+00 .5

3

1.5

2.000е-02

Рис. 15. Сравнение картин течения (плотности) для случая, посчитанного с 1-м порядком аппроксимации по времени на самой мелкой сетке и случая, посчитанного со 2-м порядком аппроксимации по времени на более грубой сетке

Fig. 15. Comparison of the flow visulisation (density) for the case calculated with the 1st order approximation in time on the finest mesh and the case calculated with the 2nd order approximation with respect to time on a coarser mesh

Kraposhin M.V. Study of capabilities of hybrid scheme for advection terms approximation in mathematical models of compressible flows. Trudy ISP RAN / Proc. ISP RAS, vol. 28, issue 3, 2016, pp. 267-326.

80 70 60 % 50

"O CD

40

Ul

30 20 10 0

-Linear speedup * Speedup, fine mesh X Speedup, very fine mesh > <

X

X /

/С

Xу / hi;.

/ *

0 10 20 30 40 50 60 70 N0 с^ СРи

Рис. 16. Сравнение масштабируемости разработанного решателя для сеток 1 млн

ячеек и 4 млн ячеек

Fig. 16. Scalability of the developed solver for 1 million cells and 4 million cells grids

6.5 Течение в сверхзвуковом сопле при наличии прямого скачка уплотнения в закритической части

Рассматривалась задача о течении в простейшем одномерном сверхзвуковом сопле, геометрия которого задавалась комбинацией двух усеченных конусов. Начальные данные соответствовали расчетному случаю из [6]. Результаты сравнивались с приближённым аналитическим решением, основанным на законах изоэнтропического течения идеального газа и теории прямых скачков уплотнений [4,6], и с расчетом в АК8У8Р1ш®упат1с8.

Схема сопла и структура установившегося течения изображены на рис.17. Отношение площадей на входе и на выходе к критическому сечению принималось равным 3, длина сопла равной 2м (для удобства оперирования с обезразмеренной координатой скачка).

Граничные условия определялись давлениями на входе и на выходе, которые принимались равными 300 и 175 кПа соответственно, на стенках ставилось ГУ проскальзывания. Течение принималось идеальным.

Поскольку существующее аналитическое решение справедливо лишь для одномерного случая (идентичные параметры течения по всему поперечному сечению сопла), была выбрана одномерная расчетная сетка (по одной ячейке в направлениях ОУ, 02; ось ОХ расположена по оси симметрии сопла). Количество ячеек по X - 100.

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

м

Рис. 17. Картина установившегося течения в сопле

Fig. 17. Visualization of the stationary supersonic flow in the nozzle

Сравнение распределения числа Маха по длине сопла с аналитическим решением представлено на рис. 18. Графики практически полностью совпадают, исключая небольшие расхождений в области скачка уплотнения, которые могут быть объяснены схемной диффузией.

2.5

1.5

iE i e

ъ 1

0,5

-Analytical ■ ■ ■ pisoCentralFoam У А

/

/

/

/

/

У *

-1 -0,8 -0,6 -0,4 -0,2 0 0,2 0,4 0,6 0 8 1

x,m

Рис. 18. Сравнение распределения числа Маха по длине сопла. Fig. 18. Comparison of the distribution of the Mach number along the length of the nozzle.

6.6 Дозвуковое течение ламинарного вязкого потока в канале круглого сечения (течение Пуазейля)

С помощью данного расчётного случая проверяется корректность воспроизведения диффузионных слагаемых в уравнении сохранения импульса (тензора напряжений) при малых числах Маха. Поскольку для данного случая известно аналитическое решение, то можно количественно оценить разницу между точным и приближённым решениями. Для постановки задачи принимается, что число Re = 200, вязкость вычисляется из физических данных, задаваемых при постановке начальных и краевых условий. Принимается, что граничные условия соответствуют нормальным условиям: скорость на входе Uвх=0,68369м/с; давление на выходе рвых=101325Па; температура — 25 0С; газ — воздух. Профиль скорости на входе в исследуемую область равномерный.

В соответствии с указанными параметрами среды число Прандтля Pr=0.73, динамическая вязкость д=1.85Т0-5 Па-с, теплоемкость Ср=1007 Дж/кг/К, молярная масса M=28.96г/моль.

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

Для получения заданного значения критерия Re, диаметр расчётной области выбирается равным 4.6мм. Длина полагалась равной 161мм. Для получения равномерной сетки расчётная область разбивается на 23 отрезка по радиусу и 1610 отрезков по длине.

Результаты сравнения численного решения, полученного с помощью настоящей модели на выходе из расчётной области, с аналитическим решением (см., например, [4,6]) для ламинарного профиля представлены на рис. 19. Число Маха составляло 0.002. Расчёт вёлся с шагом по времени 30-40 мкс, что соответствует акустическому критерию Куранта порядка 1300.

1.6

1.4

1.2

0.8

0;6

0,4

0.2

1 17 Дп' ilytica solution CentralFoam

■ ■ »Liz, pise

0.5

1,5

2.5

r.mm

Рис. 19. Сравнение аналитического и расчётного распределений Uz по радиусу канала

круглого сечения расчета.

Fig. 19. Comparison of analytical and calculated distributions of Uz along radius of channel with circular cross section

6.7 Обтекание цилиндра в ламинарном режиме

С помощью данного теста исследуется пригодность реализованной модели для моделирования дозвуковых течений для такого широко известно случая, как течение вокруг плохо обтекаемых тел. Исследование проводится для двух случаев — ламинарное течение и турбулентное течение. Исследование последнего случая особенно важно в свете поставленной в работе цели обеспечения использования уже имеющейся в ОрепБОЛМ библиотеки моделей турбулентности без её изменения.

За основу берутся результаты расчётов, полученные в работе [9]. Сравнение проводится для двух значений числа Маха - меньше 0.1 и 0.3, которые соответствуют двум предельно допустимым случаям - "глубокий" дозвук (полностью несжимаемое течение) и сжимаемое течение. Число Re равно 100. В качестве давления и температуры выбираются близкие к нормальным условия - 101325 Па и 300К соответственно, рабочая среда - воздух, молярная масса - 28.9 г/моль.

Таким образом, плотность среды при этих условиях будет равна 1.17404 кг/м3 Изобарная теплоемкость принимается равной 1004 Дж/кг/К, следовательно, показатель адиабаты равен 1.4 Скорость звука среды - a=^(yRT) = 347.6 м/с Динамическая вязкость среды взята равной 18.5 мкПа-с Число Прандтля Рг = 0.73. Диаметр цилиндр определяется по заданной скорости на входе и числу Re.

Приняв и = 10 м/с, что соответствует числу Маха 0.029, получаем значение диаметра цилиндра 0.000157м (0.157мм). Приняв и = 100 м/с, что соответствует числу Маха 0.29, получаем значение диаметра цилиндра 0.0157мм.

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

Табл. 4. Сравнение коэффициента сопротивления цилиндра полученных разными методами

Table 4. Comparison of the cylinder drag coefficient obtained by different methods

pisoCentralFoam ACL 2008-4 Sharman 05 Mene-01 Kang (2003) Ding 07

Cd 1.37 1.365 1.33 1.37 1.33 1.356

6.8 Обтекание цилиндра турбулентным потоком

Моделируется обтекание одиночного цилиндра (рис. 20), результаты сопоставляются с исследованием [10]. Помимо сопоставления с экспериментом, было выполнено сравнение с расчетом по несжимаемой и по сжимаемой дозвуковой моделям, имеющимся в OpenFOAM. Рабочая среда — воздух, условия — близкие к нормальным ( давление 101325 Па и температура 300К), плотность — 1.18 кг/м3, кинематическая вязкость — 1.5*10-5 м2/^ скорость звука — около 330 м/с. Скорость набегающего потока в экспериментах принималась равной 10 м/с, т.е. число Маха было меньше 0.1

— глубоко дозвуковое течение. В ходе экспериментов измерялись значения силы лобового сопротивления Бд и по её значению рассчитывался коэффициент лобового сопротивления Сд .

На первом этапе моделирования для отладки модели расчёты проводились в несжимаемом приближении. Расчёт проводился до наступления установившегося режима течения. Турбулентность учитывалась с использованием модели k-omegaSSТ [11]. Были рассмотрены варианты расчётных сеток с низким разрешением вблизи поверхности цилиндра (у+ ~ 100, расчёты проводились с использованием пристеночных функций), так и сеток с высоким разрешением вблизи поверхности цилиндра (у+ ~1, пристеночные функции не использовались). По полученным результатам для дальнейшего исследования был выбран второй вариант сетки (рис. 21). На границах расчётной области задавались следующие граничные условия:

1) на входе — значение скорости, интенсивности кинетической энергии и температуры (для сжимаемых случаев);

2) на выходе — значение давления;

3) на поверхности цилиндра — нулевое значение скорости (условие прилипания), пристеночная функция или нулевое значение для кинетической энергии турбулентности;

4) на верхней и нижней границах расчётной области — условие проскальзывания.

Внутри расчётной области в начальный момент температура и скорость задавались равными температуре и скорости на входе, давление — равным давлению на выходе.

При наступлении установившегося режима течения в расчётной области наблюдался периодический отрыв вихрей от поверхности цилиндра и образование за цилиндром вихревой дорожки (т. н. дорожки Кармана). Вследствие отрыва вихрей значение коэффициента лобового сопротивления Сд колеблется во времени. Поэтому в качестве результатов рассматривалось усреднённое по времени значение Сдна интервале времени после наступления установившегося режима течения.

Моделирование проводилось для случая обтекания цилиндра с диаметром 12.22 см при числе Рейнольдса ^е), равном 7.3*104 (отсюда скорость набегающего потока и= 9м/с). Интенсивность кинетической энергии турбулентности в эксперименте изменялась путем установки на входе решеток с разной формой и размером ячеек и составляла в рассмотренном случае 0.7% (рис. 22). Результирующее значение Сд в эксперименте составило 1.22.

Рис. 20. Расчётнаяобласть Fig. 20. Geometry and initial settings for the case offlow over cylinder

Рис. 21. Расчётнаясетка Fig. 21. Computational mesh

в " -. _

»

— -lu=5.4% lu=3.S%

-Iu=0.7%

О.ОЕ+ОО 2.0Е+04 4.0Е+04 6.0Е+04 8.0Е+04 1.0Е+05

Re

Рис. 22. Экспериментальные данные по обтеканию одиночного цилиндра

Fig. 22. Experimental data (drag coefficient) for the flow around a single cylinder

На втором этапе были проведены расчеты с помощью гибридного метода. Расчеты проводились на сетке с y+ ~0.8, пристеночные функции не использовались. Результаты сравнения эксперимента, расчёта в несжимаемом солвере pimpleFoam, стандартном сжимаемом решателе rhoPimpleFoam, а также в исследуемом в данной работе решателе pisoCentralFoam приведены в табл. 5.

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

Table. 5. The results of calculations of the cylinder's drag coefficient in a turbulent flow through a variety of numerical models.

№ Описание Используемый солвер Значение Cd

1 Эксперимент - 1,22

2 Расчет в несжимаемом решателе ( у+~0.8, 1шЫ=0о/о) pimpleFoam 1,15

3 Расчёт на мелкой сетке в сжимаемом решателе pisoCentralFoam 1,07

4 Расчёт на мелкой сетке в сжимаемом решателе rhoPimpleFoam 1,12

6.9 Течение струй газов со смешением

Предложенный в [1] гибридный метод может быть расширен на случай движения многокомпонентной среды совершенных газов. В этом случае к общей системе уравнений из [1] добавляются уравнения переноса (баланса) массы каждой из компонент. С учётом диффузионного приближения [16] уравнение переноса i-й компоненты смеси приобретает вид:

Данный подход был реализован в виде отдельного решателя ОрепРОЛМ и протестирован на простейшей задаче ламинарного смешения двух газов. Для сравнения были взяты результаты расчётов, полученные с помощью схемы ЛШЫ [12].

В данном случае рассматривается смешение двух различных газов текущих в плоском канале (рис. 23). Свойства газа в верхней струе соответствуют азоту (N2), текущему со скоростью 0.1м/с, нижней струи — водороду (Н2), текущему со скоростью 0.3м/с. Длина расчётной области составляет 1.2м, высота — 0.16м. Давление на выходе — 100кПа, температура обеих струй на входе — 300К. Моделирование производится до достижения стационарного состояния. Физические свойства сред были взяты следующие:

• Азот N2: молярная масса 28 г/моль, удельная изобарная теплоёмкость 1040 Дж/кг/К, число Прандтля 0.73, динамическая вязкость 17 мкПа-с, коэффициент теплопроводности 0.026 Вт/м/К.

• Водород Н2: молярная масса 2 г/моль, удельная изобарная теплоёмкость 15000 Дж/кг/К, число Прандтля 0.73, динамическая вязкость 8.9 мкПа-с, коэффициент теплопроводности 0.172Вт/м/К.

Описание граничных условий и их типы приведены в табл. 6. Сравнение расчётов с результатами моделирования [12]. Распределения плотности смеси и скорости смеси построены на расстоянии 0.7м от входа в расчётную область и представлены на рис. 24.

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

Ь\ = ».1 m/s. и2 = 0.3 m/s

Рис. 23. Схема течения в рассматриваемой задаче о смешении двух газов Fig. 23. Geometry and settings for the considered problem offlow with mixing of two gases

Таблица 6. Граничные условия

Table 6. Boundary conditions

Поле Вход азота (1п1еИ) Вход водорода (ш1е12) Выход смеси (outlet) Стенки (walls)

Давление, Pa Условие постоянного потока Условие постоянного потока Полное давление 100кПа Условие непроницаемости

Массовая доля, N2 Фиксированное значение, 1 Фиксированное значение, 0 Условие свободного выхода Условие непроницаемости компоненты

Массовая доля, H2 Фиксированное значение, 0 Фиксированное значение, 1 Условие свободного выхода Условие непроницаемости компоненты

Скорость среды/компоненты, m/s Фиксированное значение, 0.1 Фиксированное значение, 0.3 Условие свободного выхода Условие проскальзывания

Температуры среды/компоненты, K Фиксированное значение, 300 Фиксированное значение, 300 Условие свободного выхода Условие адиабатичности

U_OpenFOAM-3.0.1 U_ref

rh o_0 pen FOAM-3.0.1

■U_OpenFOAM-2.3.0

rho_ref ■ Rho_OpenFOAM-2_3.0

0.4 0,35 0,3 0,25 0,2 0,15 0,1 0,05 0

-J -к

3 Ч /

i^rift 1 и а . в

У

0,4

0,2

-0,06 -0,04

0 V. m

0,04 0,06

Рис. 24. Сравнение распределения плотности (ось справа) и скорости смеси (ось слева), полученные для задачи смешения струй различными методами и в различных пакетах (U_OpenFOAM-3.0.1 — скорость смеси, гибридный метод, OFv3.0.1, U_OpenFOAM-2.3M — скорость смеси, гибридный метод, OFv2.3.0, Uref— скорость смеси, AUSM, rho_OpenFOAM-3.0.1 — плотность смеси, гибридный метод, OFv3.0.1, rho_OpenFOAM-2.3.0 — плотность смеси, гибридный метод OFv2.3.0, rhoref— плотность смеси, AUSM)

Fig. 24. Comparison of density distribution (axis on the right) and speed (axis on the left) of

mixture for the problem of mixing jets, computed using different methods and packages (U_OpenFOAM-3.0.1 — speed of mixture, hybrid method, OFv3.0.1, U_OpenFOAM-2.3.0

— speed of mixture, hybrid method, OFv2.3.0, Uref— speed of mixture, AUSM, rho_OpenFOAM-3.0.1 — density of mixture, hybrid method, OFv3.0.1, rhoOpenFOAM-2.3.0 — density of mixture, hybrid method OFv2.3.0, rhoref— density of mixture, AUSM)

6.10 Моделирование распространения акустических волн

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

Тестирование гибридного метода на задачах распространения акустических возмущений позволяет определить степень его пригодности для решения подобных задач, или иными словами - «близость» свойств метода к «исходному» методу Курганова-Тадмора в акустическом диапазоне. Кроме того, тестирование позволит качественно «оценить» дисперсионные и диссипативные свойства схемы [13].

Для тестирования были отобраны два случая с аналитическим решения однородного волнового уравнения, соответствующие монополю и диполю [14].

6.10.1 Моделирование акустических волн, производимых «дышащей» сферой

В первом случае рассматриваются колебания плотности (давления) среды, производимые изменением радиуса сферы по гармоническому закону, рис. 25а . Акустическое давление, возникающее в этом случае вычисляется согласно следующему закону [15]:

А

рМ)

^j(ujt-kr) -jka

где амплитуда вычисляется как А — росГ'оое при условии совпадения направления излучения с радиус-вектором наблюдателя.

Рис. 25а. Схема излучения акустических

волн монополем

Fig. 25а. Diagram of acoustic wave emission by monopole

Рис. 25б. Расчетная схема для моделирования монополя.

Fig. 256.Computational domain of monopole

Для построения расчётной модели, была сформирована область, включающая в себя сегмент поверхности сферы. На «узкой» стороне сегмента задавалась скорость движения поверхности сферы, на противоположной стороне — условие свободного распространения волны. На остальных гранях — условие симметрии, рис. 25б. Далее в расчётной области выбиралась точка в которой сравнивалась динамика пульсации давления с аналитическим выражением. Константы среды и параметры движения выбирались следующие:

Необходимо понимать, что аналитическое решение было получено для распространения акустических волн в стоячей среде (однородное волновое уравнение), в то время как численное решение получено для уравнений Эйлера. Это означает, что между двумя решениями должен присутствовать сдвиг по времени и, значит, для их корректного сравнения требуется совместить оба временных ряда данных. Различие между аналитической и расчётной временными зависимостями обусловлено первым периодом колебаний, на котором численная модель вынуждена преодолевать «начальную инерцию» среды.

Было определено, что диссипативные и дисперсионные свойства схемы становятся пренебрежимо малыми при использовании более чем 20 ячеек на длину волны и акустическом числе Куранта меньшем 1/2. В этих условиях решение стремится к аналитическому, рис. 26.

Рис. 26. Сравнение аналитического и численного решения в различных точках и с различным сеточным разрешением для случая пульсирующей сферы.

Fig. 26. Comparison of analytical and numerical solutions at different points and for different mesh resolution for the case ofpulsating sphere

6.10.2 Моделирование акустических волн, производимых колеблющейся сферой

Второй тест является расширением предыдущего на трёхмерный случай (рис 27а). Расчётная область представляет собой «дольку» шара из которой вырезан шар малого диаметра.

На поверхности последней задан гармонический закон колебаний вдоль выбранной оси, на внешней поверхности — условие свободного распространения волны. Результаты сравнения представлены на рис. 28, качественная картина распространения волн дана на рис. 27б. Как видно из рисунка, присутствует расхождение по фазе между численным и аналитическим решением, которое как и в первом случае объясняется начальными условиями и соответствующим им начальным возмущениям. Параметры расчётной сетки (разрешение) были подобраны на основе предыдущего теста. Хорошее совпадение с аналитическим решением при задании соответствующего сеточного разрешения, выбранного на основе выводов из предыдущего раздела, позволяет говорить о принципиальной возможности прогнозирования погрешности численной схемы в зависимости от шага по пространству и времени. Данное свойство является критическим при использовании что в промышленных приложениях.

>1

Рис. 27а. Схема излучения акустических волн диполем - «дрожащей» сферой

Fig. 27а.Diagram of acoustic wave emission by dipole («trembling» sphere)

Рис. 27б. Расчетная область и результат для моделирования диполя — дрожащей сферы.

Fig. 276.Computational domain and result for modelling of the trembling sphere

Prriie5 pressure

(.n-jLi'e б v 1=й, '1 v priiit I ir=ii,j:

.......................J..................... г"ч (r*E]»4)

sq (г=б.Ч)

о a.iii o.ul: u.ih li.a-i u,aa

TJl№t С

Рис. 28. Сравнение аналитического и численного решения в различных точках для

случая дрожащей сферы.

Fig. 28. Comparison of analytical and numerical solutions at different points for the case of

trembling sphere

6.11 Истечение струи газа из сверхзвукового сопла

Объектом исследования в данном расчётном случае является плоское двумерное сверхзвуковое сопло [17,18], истечение через которое детально исследовалось как экспериментальным, так и расчётным способом в NASA, рис. 29а и 29б. Коэффициент расширения сопла (отношение площади выходного сечения к критическому) составляет 1.797, проектное отношение давления - 8.78. Эксперименты показали, что при отношении давлений меньше проектного, поток перерасширен и подвержен неустойчивым срывам пограничного слоя, сопровождаемым образованием ударных волн.

Исследуемой средой является сухой воздух при температуре Т=294.44К. Физические свойства следующие - число Прандтля 0.7, динамическая вязкость ^=18.27 мкПа-с, удельная изобарная теплоёмкость Ср=1005 Дж/кг/К, молярная масса 28.96кг/кмоль, коэффициент теплопроводности 1=0.024 Вт/м/К.

Были построены сетки трёх типов:

• Двумерная четырёхугольная сетка, порядка 200 тысяч элементов.

• Двумерная треугольная сетка, порядка 800 тысяч элементов.

• Трёхмерная тетраэдральная сетка, порядка 5500 тысяч элементов. Расчётная область включала в себя верхнюю часть сопла (рассматривалось симметричное течение) и объём на выходе. Сеточное расширение выбиралось равномерным. Турбулентные свойства течения воспроизводились с помощью к-ю-$$Т [11] модели и пристеночных функций, величина у+ не превышала 3. В качестве граничных условий задавались полное давление и температура на входе в расчётную область и полное давление на выходе. Для моделирования выбран режим с отношением давлений 5.02.

Для учёта смены режима на выходной границе применялось смешанное граничное условие - если число Маха в прилегающих ячейках было равно или превышало 1, то задавалось условие свободного выхода, иначе - полное давление среды на бесконечности. Расположение и обозначение граничных условий показано на рис. 30, описание граничных условий дано в табл. 7. Качественная картина течения показана на рис. 31 для двух различных типов сетки (четырёхугольная и треугольная). Сравнение расчётных данных [17] с гибридным методом показано на рис. 32, сравнение расчётных и экспериментальных данных ([18]) показано на рис. 33.

Рисунок 29а. Сопло NASA Fig. 29a.NASAnozzle

Рисунок 29б. Схематическое изображение перерасширенного сопла с отрывом

Fig. 29b. Schematic representation of over-expanded nozzle with separation

Сопоставление экспериментальных и расчётных данных показывает их хорошую согласованность для сеток трёх разных типов, что позволяет говорить о наличии сеточной сходимости. Положительным аспектом является «схожесть» картин течения, полученных на четрырёхугольной и треугольной сетках с примерно одинаковой длиной рёбер, позволяя надеяться на независимость получаемого результата от топологии.

Рис. 30. Схема расположения граничных условий и их наименования в модели истечения сверхзвуковой струи в сопле NASA

Fig. 30. Geometry, settings and boundary conditions in the model of flow expansion in a supersonic jet in a NASA nozzle

Таблица 7. Список граничных условий в задаче истечения газа из сверхзвукового сопла

Table 7. List of boundary conditions for the problem of expansion of gas from a supersonic nozzle

Граница/ Поле P, Pa U, m/s T, K k omega

inlet Полное давление, 101325Па Условие свободного входа 294.44 Интенсивность турбулентности Турбулентная длина смешения

nozzle-walls Условие непроницаемости Условие прилипания Условие адиаба- тичности Пристеночная функция Пристеночная функция

symmetry-plane Условие симметрии

vert-outlet Условие свободного выхода или полное давление Условие свободного выхода Условие свободного выхода Условие свободного выхода Условие свободного выхода

horiz-outlet Условие свободного выхода или полное давление Условие свободного выхода Условие свободного выхода Условие свободного выхода Условие свободного выхода

left-walls Условие непроницаемости Условие прилипания Условие адиаба- тичности Пристеночная функция Пристеночная функция

Рис. 31а. Мгновенное поле числа Маха и линии тока, раскрашенные по полю плотности при истечении струи газа из сопла NASA, гексаэдральная сетка

Fig. 3^. Instantaneous field of Mach number and stream lines colored according to field of density during jet expansion from a NASA nozzle, hexahedral mesh

Рис. 31б. Мгновенное поле числа Маха и линии тока, раскрашенные по полю плотности при истечении струи газа из сопла NASA,тетраэдральная сетка

Fig. 31б. Instantaneous field of Mach number and stream lines colored according to field of density during jet expansion from a NASA nozzle, tetrahedral mesh

Рис. 32. Сравнение визуализации числа Маха в сопле, полученной численным методом в [17] (слева), и в настоящей работе (справа)

Fig. 32. Comparison of visualizations ofMach number in a nozzle, calculated by the numerical method in [17] (left) and in this paper (right)

V

Рисунок 33. Сравнение экспериментального распределения давления по стене сопла (оранжевые ромбы) и расчётного (синяя линяя).

Fig 33. Comparison of experimental (orange diamonds) and calculated (blue line) distribution of pressure on the nozzle wall

6.12 Истечение квазиравновесной расширяющейся струи плазмы в вакуум

Одним из важных приложений, в которых гибридный метод может оказаться эффективным, является моделирование механики плазмы. Этот тип течений характеризуется большими градиентами плотности, скорости, давления и температуры как в пространстве, так и во времени (пульсации). Впрочем, даже в случае стационарных задач эти течения характеризуются экстремально широкими диапазонами изменения термодинамических параметров — в сотни и тысячи раз, что несомненно является «вызовом» для численных методов. Для тестирования гибридного метода был взяты результаты эксперимента и соответствующих расчётов [19,20], выполненных средствами коммерческого пакета Fluent. Данный эксперимент проводился в университете г. Эйндховен и позволяет оценить применимость приближения сплошной среды для указанного класса задач.

Опуская подробности формирования высокотемпературной плазмы, расчётная область представляет собой цилиндрический канал малого сечения, который расширяется в цилиндрический канал широкого сечения, заканчивающийся цилиндрической щелью, так что в разрезе стенки широкой части образуют «уступ» (рис. 34).

В канал малого сечения поступает горячая плазма Ar при M =1, которая, переходя в канал большого сечения, расширяется до давления 20-100Па и, огибая уступ, проходит к выходу. Для сравнения с расчётом из эксперимента доступны данные по распределению скоростей и температуры на оси канала. Диаметр последнего участка перед расширением составляет 6мм, его длина -10 мм, из которых последние 5 мм приходятся на конический раструб с углом раствора 450, так что диаметр канала перед областью свободного течения -16мм Согласно постановке задачи [19] известны температура газа плазмы (9283К, или 0.8эВ), число Маха (M=1) и расход — 3 SLM. Исходя из этих параметров вычисляется скорость потока на входе и давление.

Рис. 34. Схематическое описание установки для генерации плазмы.

Fig. 34. Schematic representation of the plasma generation device

Качественная картина течения приведена на рис. 35. Сравнение экспериментальных и расчётных картин течения показано на рисунке 36. Данная задача решалась в стационарной постановке, что позволило сэкономить время вычислений в несколько раз даже для такого относительно сложного случая с изменением температуры и давления в 10 и 100 раз. Результаты расчётов показывают в первую очередь хорошее совпадение с математической моделью, заложенной в коммерческий пакет Fluent. При этом расхождение с экспериментальными данными можно объяснить не только некорректностью выбора предположения о сплошности среды, но и граничными условиями, неопределённость которых связана с погрешностью экспериментальных данных. Последняя же составляет порядка 10% по заявлению самих авторов эксперимента.

Рис. 35. Поле плотности (в логарифмическом масштабе) плазмы Ar при стационарном расширении в вакуум 20Па

Fig. 35. Field of density (in log scale) of Ar plasma during stationary expansion into 20Pa

vacuum

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

Fig. 36a. Distribution of axial speed ofplasma, computed using hybrid method. Plotted in log

scale.

Рис. 36б. Распределение осевой скорости плазмы, полученное экспериментальным способом в работе [19]. Построеновлогарифмическоммасштабе.

Fig. 36б. Distribution of axial speed ofplasma, measured in experiment in [19]. Plotted in

log scale.

Рабочая область

О Q.

О ^

(J

- H 1

\

\ - H

г

......\ .........[.........

f ......... \ ......... .........

г

----

.........1.........

- - -- - -

0.1 0.2 Расстояние от входа в сопло, м

Чт. авг. 18 13:59:57 2016

Рис. 36в. Распределение осевой скорости плазмы, полученное расчётным способом с помощью гибридного метода. Построеновлинейноммасштабе.

Fig. 36в. Distribution of axial speed ofplasma, computed using hybrid method. Plotted in

linear scale.

Рис. 36г. Сравнение распределений осевой скорости плазмы, полученных экспериментальным способом в работе [19] и расчётно-теоретическими средствами в работе [20]. Построеновлинейноммасштабе

Fig. 36г. Comparison of axial plasma velocity between experiment [19] and numerical simulation [20]. Plottedinlinearscale.

6.13 Моделирование течения в высокоскоростном компрессоре

Гибридный метод был также доработан для моделирования устройств с подвижными частями. Использование стандартного алгоритма PISO позволило выполнить процедуру интегрирования уравнений на подвижной сетке аналогично имеющимся в OpenFOAM моделям[21]. Реализованный метод был валидирован на стандартном тесте ERCOFTAC, см. [22, 23]. На рис. 37, представлены результаты сравнения разработанного метода, стандартной модели OpenFOAM и эксперимента для обезразмеренного давления в радиальном зазоре компрессора (рис. 37б).

Рис. 37а. Схема компрессора ERCOFTAC Fig. 37a. Diagram of an ERCOFTAC compressor

B=s

as S 0.5

OA

0.3 0.2

VTI-O.o

¡sípBFlmenai -¡U--oDyMFrram (Ratit-Msson)-plseCeirtralFeam -

▼ -r-

T )». ''.........

w

■ ■ ■

1

we¡

Рис. 37б. Сравнение профилей давления в зазоре компрессора ERCOFTAC, полученных расчётным и экспериментальным способом

Fig. 37b. Comparison of computed and experimental pressure profile in an ERCOFTAC

compressor opening

Далее с использованием гибридного метода и алгоритма POD были получены и исследованы характерные моды и частоты потока в типичном высокоскоростном микрокомпрессоре [24], конструкция которого показана на рис. 38. В настоящее время статистические методы анализа, схожие с POD, являются аналогами модальных методов анализа конструкций, применяемых для выявления собственных форм и частот колебаний конструкций. Совместное использование таких статистических методов анализа для конструкций и проточных частей позволит избежать возникновения резонансных явлений и связанных с ними разрушений.

С другой стороны, эти методы могут использоваться для тестирования численных схем — если для заданного течения заранее известны моды и соответствующие им частоты, то их отсутствие может вполне сигнализировать о недостаточном качестве численной схемы. Например, анализ мод течения в модели компрессора показал наличие характерных частот, соответствующих лопаточным частотам ротора и статора — моды №0, №5 и №7 или оборотным частотам, моды №1 и №3 (рис. 39). Визуализация соответствующих мод показана на рис. 40 для мод №0 и №3. Статистический анализ показывает, что даже на относительно грубой сетке основные частоты, характерные для данного типа машин, — оборотная, лопаточная импеллера и лопаточная диффузора могут быть разрешены с использованием гибридного метода.

Рис. 38. Схематическое изображение проточной части рассматриваемого компрессора (слева) и построенной блочно-структурированной сетки (справа).

Fig. 38. Schematic presentation of flow channel in the considered compressor (left) and consequent block-structured mesh (right)

10 10 10 . 10 10 ; ю 10 10 10 10

m -

— mode 0 - mode 1 - mode 3 mode 5 — mode 7 mode 9

л /

\J

- 4 v 4

V

50000

150000

200000

100000 Frequency (HzI

Рис. 39. Характерные моды течения в микрокомпрессоре. Синими вертикальными линиями отмечены частоты относящиеся к ротору, зелёными — относящиеся к статору. Частотавращения 108000об/мин

Fig. 39. Characteristic modes of flow in a microcompressor. Frequencies related to rotor are marked with blue vertical lines, frequencies related to stator with green vertical lines. Rotation frequency is 108000rpm.

Рис. 40а. Мгновенное (слева) поле давления в компрессоре и его 0-я мода (справа) в

области импеллера

Fig. 40a. Instantaneous field ofpressure in the compressor (left) and its zeroth mode (right)

in the impeller area

Рис. 40б. Мгновенное (слева) поле давления в компрессоре и его 3-я мода (справа) в

области диффузора

Fig. 40b. Instantaneous field ofpressure in the compressor (left) and its third mode (right) in

the diffusor area

Тестирование модели на задаче ERCOFTAC позволяет говорить о качественно удовлетворительном совпадении результатов расчётов и эксперимента. При этом наблюдающиеся различия могут быть связаны с рядом причин, которые требуют рассмотрения в отдельном исследовании - «трёхмерность» течения, влияние «модельного» числа Маха на расчёт, выбор модели турбулентности и пр.

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

6.14 Модель гидродинамики водокольцевого насоса

Ещё одним важным направлением, заслуживающим интереса, является возможность применения гибридного метода для моделирования двухфазных течений в гомогенном сжимаемом приближении. Такие модели могут быть полезными для первичной оценки интегральных характеристик устройств со смешением потоков сред с большим отношением плотностей (например, вода и воздух). Кроме того, сжимаемые модели позволяют оценить пульсации давления и следовательно, уровень шума, что также является актуальной инженерной задачей в настоящее время. В случае течения гомогенной двухфазной смеси система уравнений баланса массы, импульса и энергии смеси дополняется уравнением переноса массы одной из компонент, например воды (1^), выраженной в массовой доли этой фазы ¿у: дрУид ,

dt

(UpYLiq) =

= 0

алгебраическим уравнением сохранения массовых долей всех фаз:

и алгебраическим выражением для плотности смеси:

¥Ьга . ^Ъаз \

Р ~

+

Pbiq PGas )

а также выражением для изоэнтропийной скорости звука среды (смеси):

Приведённая функция является нелинейной относительно давления, что создаёт дополнительные трудности, поскольку даже при относительно небольших скоростях ( порядка 10м/с) течение может становиться «около-» или даже «сверхзвуковым», приводя к появлению скачков плотности.

Важным приложением таких гомогенных моделей может являться моделирование водокольцевых насосов [26], которые используются в энергетике для создания разрежения высокой степени.

Принцип работы насоса (рис. 41) относительно прост и базируется на двух законах — законе сохранения массы и законе сохранения импульса. Ротор машины размещён с эксцентриситетом относительно статорной части, имеющей цилиндрическую форму. При вращении ротора жидкость в рабочей части за счёт центробежных сил «разбрасывается» к периферии, образуя между валом ротора и межфазной поверхностью жидкий кольцевой канал переменного сечения. При проталкивании прокачиваемой среды (газа) лопастями ротора через расширяющуюся часть жидкого кольцевого канала, происходит расширение среды и как следствие — создаётся разрежение на всасе. Затем кольцевой канал сужается и проталкивание газа через него приводит к росту давления на выходе из насоса.

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

входящий поток среды, красными точками — уходящий поток среды. 1 — лопасти ротора, приводящие в движение жидкость 2, создающую канал переменного сечения между поверхностью ротора и межфазной границей. 3 — статорнаячастьнасоса.

Fig. 41. Diagram of operation and flow structure in a liquid ring pump. Liquid forming the converging-diverging channel is shown with blue, passage for pumped gas with white, incoming flow of gas with green points, outgoing flow of gas with red points. 1 — rotor blades which move the liquid 2, which forms variable cross-section channel between the rotor surface and phase boundary. 3 — stator part of the pump.

В работе в качестве модели для тестирования кода была выбрана конструкция близкая к реальной, переданная Dr.-Ing. Jörn Beilke [25]. Для соединения подвижных и неподвижных частей модели использовались поверхности интерполяции данных. Для этого между соответствующими частями создавался зазор, который выбирался либо исходя из конструкторской документации, либо из соображений снижения времени расчёта (чем тоньше слой, тем больше время счёта). В начальный момент времени расчётная область была «залита» жидкостью согласно её предполагаемому положению при работе на номинальной мощности. Скорость вращения вала увеличивалась ступенчатой функцией от 0 до 200 рад/с, давление на всасе снижалось с 100кПа до 60кПа. В результате расчёта были получены распределения полей давления, скорости и объёмной и массовых долей в водокольцевом насосе (рис. 41). Сделана оценка подачи при скорости вала 200 рад/с, перепаде давлений 40 кПа — 5 м3/ч. Сравнивая данную оценку с экспериментальной величиной для перепада 40 кПа — 16м3/ч при скорости вращения вала 298 рад/с (см [25]), можно сделать вывод о качественно правильном воспроизведении явлений с помощью данной модели, поскольку:

• переход с частоты вращения 298рад/с до 200рад/с при сохранении перепада должен снизить подачу по крайней мере на 1/3;

• величина зазора между вращающимся ротором и подводящими/отводящими патрубками значительно меньше чем заданная в модели (в несколько раз), что очевидно сказывается на увеличении модельных протечек.

Качественный анализ течения (рис. 42) показывает сжатие-расширение среды в водокольцевом зазоре, наличие трансзвуковых зон в областях с объёмной долей воздуха около 50%.

Исходя из предварительного анализа результатов расчётов можно утверждать применимость модели для решения задач подобного класса.

Рис. 42. Качественная картина течения в водокольцевом насосе. Вверху слева: поле статического давления. Вверху справа: поле объёмной доли жидкости. Внизуслева: эффективное число Маха. Внизусправа:полемодуляскорости.

Fig. 42. Qualitative analysis offlow in a liquid ring pump. Top left: field of static pressure. Top right: field of volumetric fraction of liquid. Bottom left: effective Mach number. Bottom

right: magnitude of speed.

7. Заключение

1. Сформирован список задач для оценки адекватности гибридного метода моделирования сжимаемых течений. Список включает в себя задачи с аналитическим решением или экспериментальными данными и относящиеся к следующим прикладным направлениям: сжимаемые транс- и сверхзвуковые течения, распространение акустических волн, несжимаемые течения, течения газа и плазмы с высоким отношением давления (5), струйные течения, многокомпонентные и двухфазные течения.

2. Тестирование на задачах с известным аналитическим решением показало, что по крайней мере в областях сверхзвукового или околозвукового течения свойства гибридного метода слабо отличаются от метода Курганова-Тадмора, взятого за основу. В области несжимаемых течений поведение гибридного метода полностью соответствует стандартным схемам типа PISO/SIMPLE.

3. Метод позволяет решать задачи с распространением акустические колебаний, разрешать сжимаемые двухфазные или многокомпонентные течения, исследовать движение плазмы. Продемонстрирована возможность использования метода для решения задач с подвижной расчётной сеткой (в неинерциальной системе отсчёта).

Список литературы

[1]. M. Kraposhin, A. Bovtrikova, S. Strijhak. Adaptation of Kurganov-Tadmor Numerical Scheme for Applying in Combination with the PISO Method in Numerical Simulation of Flows in a Wide Range of Mach Numbers. Procedia Computer Science, 66:43-52, 2015

[2]. OpenFOAM: http://openfoam.org/

[3]. J.D. Anderson, Jr. Modern Compressible Flow: With Historical Perspective. New York: McGraw-Hill, third edition, 2003

[4]. F.M. White. Fluid Mechanics. McGraw-Hill Book Co., New York, NY, third Edition, 1994

[5]. H.E. Smith. The Flow Field and Heat Transfer Downstream of a Rearward Facing Step in Supersonic Flow. Technical report ARL 67-0056. Aerospace Research Laboratories, Ohio, 1967 (Mar.)

[6]. ANSYS Fluid Dynamics Verification Manual, Release 15.0, 2013

[7]. G.D. Garrard, W.J. Phares. Calibration of the PARC Program for Propulsion-Type flows. AEDC-TR-90-7, July, 1990

[8]. М.П. Галанин, Е.Б. Савенков. Методы численного анализа математических моделей. М.: ИздательствоМГТУим. Н.Э. Баумана, 2010

[9]. C. Liang. High-order accurate simulation of low-Mach laminar flow past two side-by-side cylinders with Spectral Difference method. Report ACL 2008-4 Aerospace Computing Laboratory, Aeronautics and Astronautics , Stanford University , May 2008

[10]. X. Liu. Wind loads on multiple cylinders arranged in tandem with effects of turbulence and surface roughness. Master thesis, Department of Civil and Environmental Engineering, Louisiana State University, 2003

[11]. F.R. Menter, M. Kuntz, R. Langtry. Ten Years of Industrial Experience with the SST Turbulence Model. Turbulence, Heat and Mass Transfer 4: Proceedings of the Fourth International Symposium on Turbulence, Heat and Mass Transfer, Antalya, Turkey, 1217 October, 2003. Publisher: 2003 Begell House, Inc.

[12]. J.R. Edwards, M. Ling. Low-Diffusion Flux-Splitting Methods for Flows at All Speeds. AIAA Journal 1998

[13]. R.A.C. Germanos, L.F. de Souza. Analysis of Dispersion Errors in Acoustic Wave Simulations. Thermal Engineering, Vol. 5 - No 01 - July 2006

[14]. Y.-H. Kim. Sound Propagation. An Impedance Based Approach. John Wiley Sons, first edition, 2010

[15]. L.E. Kinsler. Fundamentals of acoustics. Wiley, New York, 2000

[16]. Л.Г. Лойцянский. Механика жидкости и газа. М.: Дрофа, 2003

[17]. K.S. Abdol-Hamid et al. Numerical Investigation of Flow in an Over-expanded Nozzle with Porous Surfaces. 41st AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, 2005

[18]. S.C. Asbury, C.A. Hunter. Static Performance of a Fixed-Geometry Exhaust Nozzle Incorporating Porous Cavities for Shock-Boundary Layer Interaction Control. NASA Langley Research Cente, 1999

[19]. R. Engeln, S. Mazouffre, P. Vankan, D.C. Schram, N. Sadeghi. Flow dynamics and invasion by background gas of a supersonically expanding thermal plasma. Plasma Sources Sci. Technol. 10 (2001) 595-605

[20]. S.E. Selezneva, M.I. Boulos, M.C.M. van de Sanden, R. Engeln, D.C. Schram. Stationary supersonic plasma expansion: continuum fluid mechanics versus direct simulation Monte Carlo method. Journal of Physics D: Applied Physics, Volume 35, Number 12, http://dx.doi.org/10.1088/0022-3727/35/12/312

[21]. H. Jasak, Z. Tukovic. Dynamic mesh handling in OpenFOAM applied to fluid-structure interaction simulations. V European Conference on Computational Fluid Dynamics, ECCOMAS CFD 2010, Lisbon, Portugal,14-17 June 2010

[22]. O. Petit, H. Nilson, M. Page, M. Beaudoin. The ERCOFTAC Centrifugal Pump OpenFOAM Case-Study. In Proceedings of the 3rd IAHR International Meeting of the Workgroup on Cavitation and Dynamic Problem in Hydraulic Machinery and Systems, Brno, Czech Republic, 2009

[23]. J.F. Combes. Test Case U3: Centrifugal Pump with a Vaned Diffuser. ERCOFTAC Seminar and Workshop on Turbomachinery Flow Prediction VII, Aussois, jan 4-7, 1999

[24]. K. Wittig. Konstruktion einer Gastubine fuer Modellflugzeuge und Dokumentation der Auslegungsrechnungen. Muenchen, 24 September 1993

[25]. Strömungssimulation Flüssigkeitsringpumpe. Projekt 1. Ingenieurburo beilke, 28.09.2015, Dresden

[26]. H. Ding, Y. Jiang, H. Wu, J. Wang. Two Phase Flow Simulation of Water Ring Vacuum Pump Using VOF Model. ASME/JSME/KSME 2015 Joint Fluids Engineering Conference, Volume 1: Symposia, Seoul, South Korea, July 26-31, 2015

[27]. M. Kraposhin, C. Brouzet, T. Dauxois, E. Ermanyuk, S. Joubaud, I. Sibgatullin. Direct numerical simulation of internal gravity wave attractor in trapezoidal domain with oscillating vertical wall. Proceedings of ISP RAS, 26(5):117-142, 2014. DOI: 10.15514/ISPRAS-2014-26(5)-6

[28]. В.А. Васильев, М.В. Крапошин, А.Ю. Ницкий, А.В. Юскин. Применение HPC-технологий для решения пространственных задач мультифизики. Вычислительныеметодыипрограммирование, 12(1):160-169, 2011

Study of capabilities of hybrid scheme for advection terms approximation in mathematical models of compressible

flows

M.V. Kraposhin <[email protected]> Institute for System Programming of Russian Academy of Sciences, Russia, Moscow, Solzhenitsyna str., 25

Annotation. The hybrid method for approximation of advective terms is proposed in order to resolve flows in the wide Mach numbers region. This hybrid method is based on the Kurganov-Tadmor (KT) scheme and projection method PISO (Pressure Implicit with Splitting Operators). To construct this method Kurganov-Tadmor scheme for convective fluxes was formulated in implicit manner together with introduced blending function which switches between compressible regime (KT) and incompressible regime (PISO) depending on local characteristics of the flow. Such hybrid scheme gives next benefits: a) implicit treatment of diffusive terms allows to remove time step restrictions imposed by this kind of processes when approximated with explicit scheme; b) implicit formulation of convective terms together with mixing between PISO and KT produces better stability relied only on the flow Courant number, removing acoustic Courant number restrictions at low Mach number flows;

c) however, acoustic flows can can also be reproduced — in this case, local acoustic Courant number must be decreased to values less the 1 in the whole computational domain; d) utilization of KT scheme as the basis for approximation of convection terms allows to achieve non-oscillating solution for both acoustic and compressible cases (Mach number larger then 0.3). In order to study hybrid method properties a set of cases with analytical solution or experimental data for different classes of flows was considered: a) compressible flows — propagation of the wave in straight channel (Sod's Problem), supersonic flow over flat wedge, supersonic flow over backward step, flow over forward step with supersonic velocities, flow in supersonic converging-diverging nozzle with shock wave; b) incompressible flows — subsonic flow of laminar viscous fluid in the channel with circle cross section, flow around cylinder in laminar and turbulent regimes, mixing of two gases in 2D flat channel; c) industrial and academic verification tests — superisonic flow of air in NASA nozzle for pressure ratio 5, expansion of stationary equilibrium hot plasma in vacuum;

d) qualitative assessment of the hybrid method adequacy for industrial cases — numerical simulation of flows in high speed micro-compressor, simulation of two-phase flow in liquid ring pump. All materials are available for public access through GitHub project https://github.com/unicfdlab.

Keywords: mathematical models, numerical simulation, numerical schemes, compressible flows, acoustics, computational fluid dynamics, open source software.

DOI: 10.15514/ISPRAS-2016-28(3)-16

For citation: Kraposhin M.V. [Study of capabilities of hybrid scheme for advection terms

approximation in mathematical models of compressible flows]. Trudy ISP RAN / Proc. ISP

RAS, vol. 28, issue 3, 2016, pp. 267-326 (in Russian). DOI: 10.15514/ISPRAS-2016-28(3)-16

References

[1]. M. Kraposhin, A. Bovtrikova, S. Strijhak. Adaptation of Kurganov-Tadmor Numerical Scheme for Applying in Combination with the PISO Method in Numerical Simulation of Flows in a Wide Range of Mach Numbers. Procedia Computer Science, 66:43-52, 2015

[2]. OpenFOAM: http://openfoam.org/

[3]. J.D. Anderson, Jr. Modern Compressible Flow: With Historical Perspective. New York: McGraw-Hill, third edition, 2003

[4]. F.M. White. Fluid Mechanics. McGraw-Hill Book Co., New York, NY, third Edition, 1994

[5]. H.E. Smith. The Flow Field and Heat Transfer Downstream of a Rearward Facing Step in Supersonic Flow. Technical report ARL 67-0056. Aerospace Research Laboratories, Ohio, 1967 (Mar.)

[6]. ANSYS Fluid Dynamics Verification Manual, Release 15.0, 2013

[7]. G.D. Garrard, W.J. Phares. Calibration of the PARC Program for Propulsion-Type flows. AEDC-TR-90-7, July, 1990

[8]. M.P.Galanin, E.B. Savenkov. [Methods of numerical analysis of mathematical models]. BMSTU, Moscow, Russia, 2010 (in Russian)

[9]. C. Liang. High-order accurate simulation of low-Mach laminar flow past two side-by-side cylinders with Spectral Difference method. Report ACL 2008-4 Aerospace Computing Laboratory, Aeronautics and Astronautics , Stanford University , May 2008

[10]. X. Liu. Wind loads on multiple cylinders arranged in tandem with effects of turbulence and surface roughness. Master thesis, Department of Civil and Environmental Engineering, Louisiana State University, 2003

[11]. F.R. Menter, M. Kuntz, R. Langtry. Ten Years of Industrial Experience with the SST Turbulence Model. Turbulence, Heat and Mass Transfer 4: Proceedings of the Fourth International Symposium on Turbulence, Heat and Mass Transfer, Antalya, Turkey, 1217 October, 2003. Publisher: 2003 Begell House, Inc.

[12]. J.R. Edwards, M. Ling. Low-Diffusion Flux-Splitting Methods for Flows at All Speeds. AIAA Journal 1998

[13]. R.A.C. Germanos, L.F. de Souza. Analysis of Dispersion Errors in Acoustic Wave Simulations. Thermal Engineering, Vol. 5 - No 01 - July 2006

[14]. Y.-H. Kim. Sound Propagation. An Impedance Based Approach. John Wiley Sons, first edition, 2010

[15]. L.E. Kinsler. Fundamentals of acoustics. Wiley, New York, 2000

[16]. L.G. Loitsiansky. [Fluid and Gas Mechanics]. Drofa, Moscow, Russia, 2003 (in Russian)

[17]. K.S. Abdol-Hamid et al. Numerical Investigation of Flow in an Over-expanded Nozzle with Porous Surfaces. 41st AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, 2005

[18]. S.C. Asbury, C.A. Hunter. Static Performance of a Fixed-Geometry Exhaust Nozzle Incorporating Porous Cavities for Shock-Boundary Layer Interaction Control. NASA Langley Research Cente, 1999

[19]. R. Engeln, S. Mazouffre, P. Vankan, D.C. Schram, N. Sadeghi. Flow dynamics and invasion by background gas of a supersonically expanding thermal plasma. Plasma Sources Sci. Technol. 10 (2001) 595-605

[20]. S.E. Selezneva, M.I. Boulos, M.C.M. van de Sanden, R. Engeln, D.C. Schram. Stationary supersonic plasma expansion: continuum fluid mechanics versus direct simulation Monte Carlo method. Journal of Physics D: Applied Physics, Volume 35, Number 12, http://dx.doi.org/10.1088/0022-3727/35/12/312

[21]. H. Jasak, Z. Tukovic. Dynamic mesh handling in OpenFOAM applied to fluid-structure interaction simulations. V European Conference on Computational Fluid Dynamics, ECCOMAS CFD 2010, Lisbon, Portugal,14-17 June 2010

[22]. O. Petit, H. Nilson, M. Page, M. Beaudoin. The ERCOFTAC Centrifugal Pump OpenFOAM Case-Study. In Proceedings of the 3rd IAHR International Meeting of the Workgroup on Cavitation and Dynamic Problem in Hydraulic Machinery and Systems, Brno, Czech Republic, 2009

[23]. J.F. Combes. Test Case U3: Centrifugal Pump with a Vaned Diffuser. ERCOFTAC Seminar and Workshop on Turbomachinery Flow Prediction VII, Aussois, jan 4-7, 1999

[24]. K. Wittig. Konstruktion einer Gastubine fuer Modellflugzeuge und Dokumentation der Auslegungsrechnungen. Muenchen, 24 September 1993

[25]. Strömungssimulation Flüssigkeitsringpumpe. Projekt 1. Ingenieurburo beilke, 28.09.2015, Dresden

[26]. H. Ding, Y. Jiang, H. Wu, J. Wang. Two Phase Flow Simulation of Water Ring Vacuum Pump Using VOF Model. ASME/JSME/KSME 2015 Joint Fluids Engineering Conference, Volume 1: Symposia, Seoul, South Korea, July 26-31, 2015

[27]. M. Kraposhin, C. Brouzet, T. Dauxois, E. Ermanyuk, S. Joubaud, I. Sibgatullin. [Direct numerical simulation of internal gravity wave attractor in trapezoidal domain with oscillating vertical wall]. Trudy ISP RAN/Proc. ISP RAS, 26(5):117-142, 2014 (in Russian). DOI: 10.15514/ISPRAS-2014-26(5)-6

[28]. V.A. Vasiliev, M.V. Kraposhin, A.Yu. Nitzkiy, A.V. Yuskin. [Application of HPC-technologies to solving spatial multiphysics problems]. Vychislitel'nye metody i programmirovanie [Numerical Methods and Programming], 12(1):160-169, 2011 (in Russian)