Математическое моделирование методом декартовых сеток задачи о взаимодействии ударной волны с системой тел тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат наук Сидоренко, Дмитрий Алексеевич

• Сидоренко, Дмитрий Алексеевич
• кандидат наук
• 2018

Оглавление диссертации кандидат наук Сидоренко, Дмитрий Алексеевич

ВВЕДЕНИЕ.................................................................................................................4

ГЛАВА 1. МЕТОДЫ ПОГРУЖЕННОЙ ГРАНИЦЫ И ДЕКАРТОВЫХ СЕТОК ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ УРАВНЕНИЙ ЭЙЛЕРА В ОБЛАСТЯХ СЛОЖНОЙ ФОРМЫ.........................................................................13

1.1. Проблема «малой» ячейки и метод h-ячеек.................................................15

1.2. Проблема «малой» ячейки и метод слияния ячеек......................................17

1.3. Проблема «малой» ячейки и метод распределения потоков......................20

1.4. Метод Forrer, Berger.......................................................................................22

1.5. Метод погруженной границы........................................................................25

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

ОБЛАКОМ ЧАСТИЦ...............................................................................................28

ГЛАВА 3. МОДЕЛИРОВАНИЕ ВЗАИМОДЕЙСТВИЯ УДАРНОЙ ВОЛНЫ С СИСТЕМОЙ НЕПОДВИЖНЫХ ЦИЛИНДРОВ (ЗАСЫПКОЙ ГРАНУЛИРОВАННОЙ СРЕДЫ)............................................................................39

3.1. Постановка задачи..........................................................................................39

3.2. Математическая модель.................................................................................42

3.3. Вычислительный алгоритм метода h-ячеек.................................................43

3.3.1. Расчет численного потока через регулярные ребра.................................45

3.3.2. Расчет численного потока через усеченные ребра..................................47

3.4. Верификация...................................................................................................51

3.4.1. Нормальное отражение ударной волны от стенки...................................52

3.4.2. Отражение ударной волны от клина.........................................................54

3.4.3. Взаимодействие ударной волны с цилиндром.........................................59

3.4.4. Взаимодействие ударной волны с регулярной решеткой цилиндров .... 61

3.5. Результаты и обобщение вычислительных экспериментов.......................63

3.5.1. Взаимодействие УВ с засыпками гранулированных сред......................63

3.5.2. Комплексный подход в проблеме взаимодействия ударной волны с

облаком частиц......................................................................................................70

ГЛАВА 4. ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ РЕЛАКСАЦИИ ТЕЛ ЗА ПРОХОДЯЩЕЙ УДАРНОЙ ВОЛНОЙ..................................................................74

4.1. Постановки задач............................................................................................74

4.1.1. Постановка задачи о взаимодействии ударной волны с подвижным цилиндром .............................................................................................................. 74

4.1.2. Постановка задачи о динамике движения двух частиц в сверхзвуковом потоке газа за ударной волной.............................................................................75

4.2. Математическая модель и вычислительный алгоритм...............................76

4.3. Верификация...................................................................................................83

4.4. Результаты вычислительных экспериментов..............................................92

4.4.1. Движение цилиндров различной массы....................................................92

4.4.2. Движение двух сфер в сверхзвуковом потоке..........................................96

ЗАКЛЮЧЕНИЕ.........................................................................................................98

СПИСОК СОКРАЩЕНИЙ И УСЛОВНЫХ ОБОЗНАЧЕНИЙ...........................100

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

Введение диссертации (часть автореферата) на тему «Математическое моделирование методом декартовых сеток задачи о взаимодействии ударной волны с системой тел»

ВВЕДЕНИЕ

Многомерное газодинамическое моделирование может использоваться для уточнения качественных и количественных характеристик ударно-волновых процессов в двухфазных средах [1]. Так, в канонической задаче о взаимодействии ударной волны (УВ) с облаком частиц [2] рассмотрение дифракции УВ на системе тел (сферической [3], Рис. 1 или более сложной формы [4], Рис. 2) позволяет описать влияние коллективных эффектов аэродинамической интерференции на коэффициент сопротивления облака. Двумя основными параметрами, от которых зависит характер протекания процесса, являются интенсивность падающей волны, определяемая ее числом Маха М, и пористость или проницаемость системы тел £, то есть доля объема, занятая газом.

Рис. 1 - Пространственные распределения давления в последовательные моменты времени в задаче о взаимодействии УВ с системой сфер различного размера [3]. М = 3.0, £ = 0.85.

Важным аспектом, влияющим на возможность подобного многомерного газодинамического моделирования (в литературе также употребляют термины «прямое численное моделирование» [5] или «моделирование на микроуровне»

[6], которые также будут использованы в работе), является вычислительный алгоритм, поскольку речь идет о решении уравнений газовой динамики в многосвязной области, которую нельзя покрыть структурированной расчетной сеткой. Так, в [7, 8] расчеты осуществлялись пакетом ANSYS Fluent; о расчетных сетках в двумерном случае известно, что они были с четырехугольными ячейками, детали устройства сетки в многосвязной области, как и какие-либо особенности вычислительного алгоритма с этим связанные, не приводятся.

Рис. 2 - Численная шлирен-визуализация в задаче о падении УВ на матрицу из

тел различной формы, М = 1.4 [4].

Рис. 3 - Поля статического давления для сферы, М = 1.5 [7]. Переход от коллективного обтекания (а) к индивидуальному с Маховским взаимодействием УВ (б) и с регулярным взаимодействием (в).

В [7] была рассмотрена решетка из неподвижных поперечно расположенных цилиндров, классифицированы типы взаимодействия отошедших УВ (см. Рис. 3), а также продемонстрировано, что волновые структуры, возникающие при обтекании цилиндра и сферы, существенно различаются. В [8] был рассмотрен случай двух неподвижных сфер (двумерные расчеты в осесимметричной постановке, см. Рис. 4), расположенных вдоль направления движения УВ; получены кривые изменения коэффициента сопротивления для каждой из сфер во времени при варьировании расстояния между сферами.

Рис. 4 - Поля плотности для случая взаимодействия УВ с двумя сферами в последовательные моменты времени, М = 3.0 [8].

Рис. 5 - Пространственные распределения давления в задаче о взаимодействии УВ с системой цилиндров, моделирующей облако частиц [9]. М = 1.67, £ = 0.85.

В [9] двумерные газодинамические расчеты использовались для того, чтобы качественно описать многомерные эффекты при взаимодействии УВ с облаком частиц в трехмерном случае, см. Рис. 5. Полученное решение исполь-

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

Актуальность темы исследования и степень ее разработанности. Различные аспекты взаимодействия УВ с телами исследуются экспериментально, аналитически и численно уже очень давно [10, 11]. Значительно позже стали изучаться процессы взаимодействия УВ с системой тел применительно к уточнению описания течения газа в гранулированной или пористой среде.

Рис. 6 - Теневые снимки процесса взаимодействия УВ с системой цилиндров

[12]. М = 1.67, 5 = 0.72.

Вероятно, первая работа подобного рода - исследование [12], в котором экспериментально рассматривалась задача о взаимодействии УВ с системой неподвижных цилиндров для получения качественного представления о характере процесса, см. Рис. 6. Также решались одномерные уравнения Эйлера с источниковыми членами в правой части для описания течения газа в пористой среде. Актуальной и до конца не решенной задачей, к которой, так или иначе, приближаются рассмотренные выше работы [3 - 9, 12], как и настоящая работа - получение, вообще говоря, нестационарного коэффициента сопротивления засыпки гранулированной среды или пористого тела в результате детального многомерного газодинамического моделирования. Данный коэффициент сопротивления необходим в моделях механики гетерогенных сред, например, [13, 14], а также в Лагранжевых моделях, в качестве замыкающего соотношения в выражении для силы межфазного трения. Используемые в настоящее время замыкающие соотношения - широко распространенной является, так называемая, формула Эргана [15] - получены, как правило, в опытах с медленными течениями газа или жидкости без УВ через загроможденное пространство. Математические проблемы, возникающие на этапе стыковок моделей на микроуровне (взаимодействие УВ с отдельными частицами) и макроуровне (течение двухфазной среды), обсуждаются в [6].

Цели диссертационной работы:

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

2. Численное исследование задачи о взаимодействии УВ с одной и двумя подвижными сферами для описания релаксации частиц за УВ. Задачи диссертационной работы:

1. Разработка двух вычислительных алгоритмов метода декартовых сеток для численного моделирования распространения УВ в областях сложной формы с подвижными и неподвижными границами.

2. Разработка комплекса программ для численного исследования распространения УВ в областях сложной формы, в том числе с подвижными границами.

3. Проведение параметрических вычислительных экспериментов по взаимодействию УВ с системой цилиндров при варьировании длины системы тел, ее проницаемости и диаметра цилиндров.

4. Проведение вычислительных экспериментов по взаимодействию проходящей УВ с одной и двумя подвижными сферами.

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

Теоретическая и практическая значимость работы. Задача взаимодействия УВ с системой тел представляет интерес с нескольких точек зрения. Во-первых, она связана с рядом практических приложений, таких как ослабление взрывных волн при авариях в добывающей отрасли [16]. Результаты вычислительных экспериментов, выполненных в работе, дают возможность оценок ослабления УВ засыпками гранулированных сред. Во-вторых, подобная постановка соответствует фундаментальной задаче Римана для систем уравнений, описывающих течения двухфазных сред. Формулировка математических моделей и вычислительных алгоритмов, позволяющих рассчитывать широкий

спектр режимов течений двухфазных сред - от разреженных до плотных - в рамках подходов механики гетерогенных сред, является сегодня очень актуальным и не до конца решенным вопросом [14, 17]. Использование многомерного моделирования позволяет получать «точные» решения, которые могут быть использованы при разработке моделей механики гетерогенных сред.

Методология и методы исследования. В работе использовались методы математического моделирования и проведения вычислительных экспериментов, в том числе на многопроцессорных ЭВМ, а именно: метод декартовых сеток, конечно-объемный метод Годунова расчета потоков, явная схема Эйлера интегрирования по времени, метод декомпозиции расчетной области.

На защиту выносятся следующие положения:

1. Два вычислительных алгоритма метода декартовых сеток для численного моделирования распространения УВ в областях сложной формы с неподвижными и подвижными границами.

2. Соответствующие программные комплексы для численного моделирования распространения УВ в областях сложной формы с неподвижными и подвижными границами.

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

4. Результаты численных исследований взаимодействию проходящей УВ с одной и двумя подвижными сферами.

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

Основные результаты диссертационной работы докладывались и обсуждались на следующих конференциях:

- 4-ая международная научная школа молодых ученых «Волны и вихри в сложных средах» (Москва, 26 - 29 ноября 2013 г.);

- 40-ая международная молодежная научная конференция «Гагаринские чтения» (Москва, 7 - 11 апреля 2014 г.);

- 10-ый Международный коллоквиум по пульсирующей и непрерывной детонации (Санкт-Петербург, 4 - 8 июля 2016 г.);

- 7-ой Международный симпозиум по неравновесным процессам, плазме, горению и атмосферным явлениям (Сочи, 2 - 7 октября 2016 г.);

- Ежегодные научные конференции отдела Горения и взрыва Института химической физики им. Н.Н. Семенова РАН (Москва, 8 - 10 февраля 2017 г., 7 - 9 февраля 2018 г.);

- 31-ый Международный симпозиум по ударным волнам (Нагоя, Япония, 9 - 14 июля 2017 г.);

- 5-ый Минский международный коллоквиум по физике ударных волн, горению и детонации (Минск, Беларусь, 13 - 16 ноября 2017 г.);

- Квазилинейные уравнения, обратные задачи и их приложения (Долгопрудный, 5 - 7 декабря 2017 г.);

- 56-ая, 57-ая, 58-ая, 59-ая и 60-ая Научные конференции МФТИ (Долгопрудный, 25 - 30 ноября 2013 г., 24 - 29 ноября 2014 г., 23 - 28 ноября 2015 г., 21 - 26 ноября 2016 г., 20 - 26 ноября 2017 г.);

- 33-ья Международная конференция по уравнениям состояния вещества (Эльбрус, 1 - 6 марта 2018 г.);

- Международная конференция «50 лет развития сеточно-характеристического метода» (Долгопрудный, 31 марта - 3 апреля 2018 г.);

- Международная конференция по методам аэрофизических исследований (Новосибирск, 13 - 19 августа 2018 г.).

Публикации автора по теме диссертационного исследования. Основные результаты опубликованы в 18 работах, включая 3 статьи в изданиях из Перечня ВАК ведущих рецензируемых научных журналов и изданий, в которых должны быть опубликованы основные результаты диссертации на соискание степеней кандидата наук [18 - 20], 1 в издании из системы Web of Science [21], 2 статьи в журнале, включенном в систему РИНЦ [22, 23], 2 публикации в сборниках статей [24, 25], 10 публикаций в трудах конференций [26 - 35].

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

Благодарности. Автор выражает глубокую признательность д.ф.-м.н. В.М. Бойко (ИТПМ СО РАН) за пояснения в постановке задачи о динамике движения частиц в сверхзвуковом потоке за проходящей ударной волной. Расчеты проводились на суперкомпьютере МВС-10П Межведомственного суперкомпьютерного центра РАН.

Глава 1. Методы погруженной границы и декартовых сеток для численного решения уравнений Эйлера в областях сложной формы.

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

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

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

криволинейными границами расчетной области. Отметим также подход, связанный с использованием системы сеток — сетки в ядре течения, которая не согласована с границами тел, и сеток в окрестности криволинейных границ, согласованных с ними [41, 42]. Согласование решений на перекрывающихся участках сеток достигается за счет процедур интерполяции. В работе [42] метод дополнен локальной адаптацией сеток и подвижностью сетки, построенной вблизи криволинейной границы, что дает возможность вести расчёт в областях с подвижными границами. Также адаптация сеток наряду с отдельной процедурой обработки усеченных ячеек успешно применяется на полностью декартовых сетках [43, 44].

Сами методы погруженной границы можно условно разделить на два класса. В методах первого класса в определяющую систему уравнений вводятся фиктивные члены, которые моделируют действие стенки. Расчет при этом ведется единообразно сквозным образом во всей расширенной расчетной области, включающей в себя ячейки, пересекаемые криволинейными границами непроницаемых объектов или целиком накрываемые этими объектами. Новый метод из этого класса для расчета течений сжимаемых сред был предложен в недавних работах [45, 46]. В [45] особое внимание уделяется виду «компенсационного» потока, который фигурирует в правых частях уравнений и который обеспечивает эквивалентность модифицированной постановки со сквозным счетом в расширенной области исходной краевой задаче. Таким образом, среди основных достоинств обсуждаемого класса методов погруженной границы — их алгоритмическая однородность, то есть единый вид разностной схемы во всех расчётных ячейках. К недостаткам можно отнести отсутствие однозначного физического смысла параметров непосредственно в ячейках, пересекаемых границами объектов, особенно в задачах физико-химической гидродинамики, а также необходимость контроля консервативности вычислительного алгоритма в условиях наличия в схеме источниковых членов.

Второй класс объединяет в себе методы декартовых сеток или методы «cut-cell», когда наряду с внутренними регулярными ячейками дискретизация исходной системы уравнений производится особым образом в граничных ячейках, пересекаемых границами расчетной области. Cut-cell методы и их развитие широко используются как в зарубежных коммерческих пакетах для численного моделирования задач механики сплошной среды, таких как Ansys CFX и Ansys Fluent, так и в российских (FlowVision, см., например, [47]). Детальное рассмотрение и учет того, каким образом криволинейная граница расчетной области пересекает прямоугольную ячейку декартовой сетки, приводит, однако, к другой проблеме, известной в литературе как «проблема малой ячейки». Рассмотрим ее далее на примере одномерного линейного уравнения переноса.

1.1. Проблема «малой» ячейки и метод h-ячеек

Рассмотрим одномерное скалярное нелинейное уравнение переноса:

Э? Эх

Будем искать его решение методом конечных объемов на сетке {х{} с «малой» ячейкой с индексом к (см. Рис. 1.1):

где и" и сеточные аналоги функций и и /, / - пространственный индекс,

п - временной, ^ - длина / -ой ячейки, т" - шаг интегрирования по времени. При использовании метода Годунова поток через ребро I -1/2 определяется выражением:

{Xi j : {Хк-1/2 = 0, Хк+1/2 = h , Xk±1/2±i = Xk±12 ± 1 ' h, 1 = 1,2,Kj

к-12

к+12

Разностная схема запишется в виде:

Куг = /(и* [<_!,<]),

* I п п \ п п

где и (и-1,и ) - решение задачи о распаде разрыва с параметрами иг-1 и и1 сле-

(п п\

и-1> и )

ва и справа от него.

Н

и

.п

к-г

Н п Щ п ик+1

п ик-1

, Н ,

ик+1/г ик+1/г

Н -^

ь ик-1/г К ик-1/г

и

п

к+г

Рис. 1.1 - К расчету потоков через ребра малой ячейки.

Поскольку вычислительный алгоритм на основе метода И-ячеек в двумерном случае будет изложен ниже (см. Главу 3), здесь остановимся на схематическом описании его в одномерном случае. В методе И-ячеек [48] при расчете потока через ребро I -1/2 в качестве аргументов задачи о распаде разрыва предлагается выбирать не ип_х и иги, а значения

К

х 1 г х1_1/г+Н п/\7 ь

Щ-Цг =Т I и(Х) ЛХ и иг_Чг ,, Н

Н'!х• -12 Н'!х'-12-Н

^_1 г — I и (х) Лх,

- -1Г

"1г Н ^

которые являются средними взвешенными функции и по отрезку длиной Н .

Рассмотрим случай линейного уравнения переноса / (и) - и , для которого численный поток в методе Годунова имеет вид:

Аг - / (и * Кг, иК/г ])

и

-1/г-

ь

Тогда, в предположении кусочно-постоянной интерполяции сеточной функции и условии к' £ к, переменная в малой ячейке методом И-ячеек найдется в виде:

т"

к

и"--1 к

ки" + (к - к) и"_х Л " (11)

и +(к - к ) и , "

, и 1-1 V к У

и" -т-(и" -и"_1),

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

т"

"+1 " " ( " " \ и = и -к(и -и-1 )•

Критерий устойчивости построенной схемы для рассматриваемого линейного уравнения переноса записывается в виде т"/к £ 1, что в случае расчета потока стандартным конечно-объемным методом приводит к сильному уменьшению шага по времени при к' □ к, тогда как использование схемы (1.1) лишает необходимости учета размера малой ячейки при расчете шага по времени.

1.2. Проблема «малой» ячейки и метод слияния ячеек

Другой метод решения проблемы малой ячейки заключается в слиянии малой усеченной ячейки с одной или несколькими соседними с целью увеличения площади расчетной ячейки как альтернативы уменьшения шага по времени в ней, и, соответственно, шага интегрирования во всей расчетной области [49, 50]. Этот подход также успешно применяется к задачам с подвижными телами [51] и задачам химической гидродинамики [52].

Рассмотрим широко используемый метод слияния ячеек второго порядка аппроксимации в двумерном случае [53]. Рассмотрим участок расчетной области, покрытый декартовой сеткой и содержащий пересекаемые ячейки. Алгоритм поиска ячеек, требующих слияния, основывается на критерии «малости» ячейки. Очевидно, что не все усеченные ячейки, а только достаточно малые из них могут являться причиной неустойчивости. Авторы работы [49] предлагают

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

А / С

<

в

Рис. 1.2 - Поиск соседних ячеек для слияния.

Итак, для каждой малой усеченной ячейки необходимо найти большую усеченную или регулярную ячейку, соседствующую с малой по ребру, и параметры газа в малой ячейке впредь рассчитывать при помощи потоков газа через ребра ячейки, получившейся в результате слияния. Например, ячейки А и В на Рис. 1.2 сливаются с ячейками С и Б соответственно.

Уравнения Эйлера записываются в интегральной форме:

Э_

£ Ш О + ^ Р • = 0,

" р'

и = ри , Р • п =

Ру

е

ри • п

рии • п + рп х руи • п + рп ^

(е + р )• и • п

и

р

7~ 1

+ Р-

2 , 2 и + V

2

Здесь обозначения стандартные, О - тело, ограниченное кривой £. Разностная схема строится следующим образом:

е

т т

(Ш =(ки )П-I £ Р (и к )■ 8;. (1.2)

2 к=1

Выписанная схема является первым неконсервативным полушагом в используемой схеме предиктор-корректор. Здесь I и у - индексы при естественной

нумерации ячеек декартовой сетки, V - площадь ячейки, 8 - вектор, нормальный к ребру ячейки. Суммирование в (1.2) ведется по ребрам ячейки, которых у регулярных ячеек 4, а у усеченных от 3 до 7. Вектор потока Р(и к) определяется

по реконструированному значению вектора консервативных переменных

ик = иг" + гк ■ Уи; на соответствующем ребре к, гк - радиус вектор середины

ребра в системе, связанной с центром масс рассматриваемой ячейки, а Уи; -

градиент консервативных переменных.

Градиенты в регулярных ячейках вычисляются тривиально. Но для обеспечения второго порядка аппроксимации в ячейках вблизи границ требуется отдельная процедура вычисления градиентов и применения ограничителей, например, [37, 54].

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

т

(ки у =(ки)"-г£ Р (и? )■ 8^2. (1.3)

к=1

Здесь Р(и£, и^) - поток вектора консервативных переменных, полученный в

результате решения задачи Римана с параметрами слева и справа, определяемыми формулами:

и £ = и ;+1/2 + гк -у и;,

и кк = и;;12 + г* -уит,

где I и т - индексы ячейки, граничащей с ячейкой /, у по ребру к.

Для решения задачи Римана может использоваться любой монотонный численный поток. Авторами [53] используется ИЬЬ метод [55].

Авторы работы [52] использовали аналогичный метод слияния ячеек, дополненный локальной адаптацией сетки для решения задач химической гидродинамики.

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

1.3. Проблема «малой» ячейки и метод распределения потоков

Еще одним методом, решающим проблему малой ячейки, является метод распределения потоков, когда сумма потоков, входящая в правую часть конечно объемной схемы вида (1.3) в малой ячейке уменьшается до значения, пропорционального ее объему за счет распределения потоков в другие ячейки [56].

На двумерной декартовой сетке с шагом h расположено плоское тело О. Пусть V - площадь рассматриваемой усеченной ячейки. Введем обозначение для доли тела в ячейке k = V/к2. Будем решать систему уравнений газовой динамики

и

р

ри ру

рЕ

эи

+ У- Р = 0,

Р = [рх Р У ]

ри

2

ри + р риу и (рЕ + р)

ру

риу

2

ру + р V (рЕ + р)

(1.4)

Рассматриваются две аппроксимации V- Р через граничное ребро. Первая является консервативной:

(V • Р )С=Ж (^ " ^+аЛ - аРУ+аР • пв)

(1.5)

где ц - отношение длины рассматриваемого ребра к И, пв - вектор внешней

нормали к граничному ребру. Вторая аппроксимация является неконсервативной, но устойчивой:

(У ■ РГ = И(Ре - Р^ + РК - Р^). (1.6)

В (1.5) и (1.6) индексы Е, ^^ К, Б и В обозначают правое, левое, верхнее, нижнее и граничное ребра соответственно. Использование консервативной аппроксимации приводит к схеме:

и"+1 = и; -М(У ■ Р)с,

которая, однако, не является устойчивой. Поэтому предлагается использовать линейную комбинацию аппроксимаций (1.5) и (1.6) для построения разностной схемы для (1.4):

и;+1 = и; -М [^(У ■ Р)с + (1 - 77)(У ■ Р)Кс ]. (1.7)

Рис. 1.3 - К методу распределения потоков.

Таким образом, при использовании схемы (1.7) возникает следующая ошибка в расчете вектора консервативных переменных, связанная с использованием неконсервативной аппроксимации:

dM = -k (1 -h) [(V ■ F )C + (V • F )NC '

Если выбрать 77 = к, то малый знаменатель в формуле (1.5) сократится, и получится устойчивый метод. Но при этом он будет неконсервативным. Идея метода распределения потоков заключается в том, чтобы разделить величину между соседними с текущей ячейками (см. Рис. 1.3):

и«+1 ® и"+1 + м>Ш, м>. > 0, V м>к. = 1,

1 1 I * I * ^^^ 11 ^

г

где / - индексы окружающих ячеек, на которые указывают стрелки на Рис. 1.3. Величины можно выбирать равными:

f л-1

w =

V i У

I k,

В работе [57] для расчета wi предлагается следующая формула:

Р™

w =———

r ki

i

где rNC - это неконсервативная оценка плотности на следующем шаге по времени в текущей ячейке:

rNC =r-D (V^ F )ГС.

Аналогичный метод используется в более поздней работе [58].

1.4. Метод Forrer, Berger

Ещё один подход к решению проблемы малой ячейки предложен авторами работы [59], а в работе [60] несложно обобщен на случай подвижных границ.

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

Список литературы диссертационного исследования кандидат наук Сидоренко, Дмитрий Алексеевич, 2018 год

// -

* ■..................

V;

■ / аГ

У- Е В -

г \ \ \ 1111 1 1

0

250 500 750 Число процессоров

1000

Рис. 4.2 - Оценка ускорения работы алгоритма. Сплошная линия - идеальное ускорение. Пунктирная линия и точки - ускорение при расчете 1000 шагов по времени. Пунктирная с точкой - ускорение при расчете 100 шагов по времени.

4.3. Верификация

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

Алгоритм тестировался на задаче о колебании поршня. Расчетная область представляет собой отрезок, на левой границе которого выставлено условие непротекания. В области покоится поршень, на который действует постоянная сила. Поршень и газ между ним и левой границей приходит в движение. Координата левой границы расчетной области 0, начальная координата поршня 1.0. Уравнение движения поршня:

э2 Хр АРр (Хр)- Р

дг2 т

где хр - координата левой границы поршня, А - площадь поршня, рр - давление газа на его левую границу, Р - сила, приложенная к поршню в направлении левой границы расчетной области, т - масса поршня. В расчете приняты значения А/т = 2 и Р/т = 4, в начальный момент времени поршень и газ покоятся, плотность и давление в газе равны 1.0. На Рис. 4.3 приведены кривые зависимости координаты левой границы поршня от времени в расчете [91] и в расчете автора на двух сетках - 50 и 250 ячеек.

Рис. 4.3 - Зависимость координаты левой границы поршня от времени. Красная пунктирная линия - расчет [91], зеленая - расчет автора на сетке 250 ячеек, синяя - расчет автора на сетке 50 ячеек.

Методика была протестирована на предмет консервативности по массе. Очевидно, что в рассматриваемом тесте масса газа, заключенного между левой стенкой и поршнем, должна оставаться постоянной. На трех сетках - 50, 250 и

1250 ячеек - на Рис. 4.4 приведены графики погрешности массы. Наблюдается уменьшение этой погрешности с измельчением сетки.

Рис. 4.4 - Зависимость погрешности полной массы от времени. Зеленый цвет сетка 50 ячеек, красный цвет - 250 ячеек, синий цвет - 1250 ячеек.

Рис. 4.5 - Изолинии плотности на момент времени 0.1. Слева - расчет [91] на сетке 800 х 800 ячеек, справа - расчет автора на сетке 400 х 400.

В той же работе [91] рассмотрен двумерный тест о течении в окрестности круглого тела диаметром 0.2 с центром в точке с координатами (0.5, 0.5), колеблющемся в покоящемся идеальном газе с показателем адиабаты 1.4 и плотностью и давлением 1.0. Движение тела описывается уравнением:

X (*) = Ус (*) = (10Р), X (0) = у (0) = 0.5.

На Рис. 4.5 приведены изолинии давления в расчете автора и из [91].

Также был рассмотрен часто встречающийся в литературе тест о подъеме цилиндра за проходящей УВ под действием сил давления. В начальный момент времени в прямоугольном канале длиной 1.0 и высотой 0.2 на расстоянии 0.15 от левой границы на нижней границе покоится цилиндр диаметром 0.1 и плотностью 10.77. Канал заполнен покоящимся идеальным газом с показателем адиабаты 1.4 с плотностью 1.0 при давлении 1.0. Все величины приведены в безразмерном виде. На цилиндр слева направо набегает плоская УВ с числом Маха 3.0. Начальная координата фронта 0.08. На левой границе расчетной области задано условие втекания газа с параметрами за УВ, на всех остальных -условия непротекания. Время расчета 0.33. Расчетная сетка 2 000 х 400 ячеек. Рис. 4.6 иллюстрирует динамику изменения положения центра цилиндра при его движении, вызванном силами давления. Там же отмечены расчетные данные из работ [60, 92 - 96]. Стоит отметить расхождение всех расчетных результатов в динамике изменения абсциссы и ординаты цилиндра. Представленные результаты говорят о том, что методики расчета подвижных тел нуждаются в дальнейшем развитии и анализе на предмет определения границ их применимости. Различия в результатах расчетов с использованием различных методик наблюдается и на поле течения. Сравнение Рис. 4.7 с аналогичными пространственными распределениями из [60, 92 - 96] показывает, что ударно-волновая картина в левой части, связанная с отражением падающей УВ от цилиндра, воспроизводится практически одинаково во всех расчетах.

(а)

(б)

Рис. 4.6 - Динамика движения центра цилиндра - (а) абсциссы и (б) ординаты -в тестовой задаче о подъеме цилиндра за УВ: 1 - [92], 2 - [93] (в расчете, выполненном в [93], было использовано три различных подхода), 3 - [60], 4 - [94], 5 -[95], 6 - [96].

0.4 0.5 0.6 0.7 0.8 0.9 Рис. 4.7 - Распределение модуля градиента плотности газа в тестовой задаче о подъеме цилиндра на финальной стадии процесса.

В то же время вихревые структуры в окрестности цилиндра сильно варьируются от методики к методике, что, вероятно, связано с использованием в данной работе и в исследованиях [60, 92 - 95] модели невязкого газа. Эта же задача была рассмотрена в работе [96] в рамках уравнений Навье-Стокса с использованием метода из класса погруженной границы. Число Рейнольдса выбиралось равным 240. Рис. 4.8 иллюстрирует рассчитанные пространственные распределения модуля градиента плотности газа в расчетах в рамках моделей невязкого и вязкого газа.

Рис. 4.8 - Шлирен-визуализация «невязкого» (а, Ь, с) и «вязкого» (ё, е, : расчетов задачи о подъеме цилиндра [96] за проходящей УВ в три момента времени: 0.0 (а, ё), 0.1 (Ь, е) и 0.3 (с, ф.

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

(а) (б)

Рис. 4.9 - Изолинии плотности на момент времени 11.6 мс в тестовой задаче о взаимодействии двух цилиндров за проходящей УВ, полученные в (а) работе [71] и (б) в настоящей работе.

В модели рассматривается неупругое нецентральное соударение двух тел произвольной формы с коэффициентами трения и восстановления при ударе. На основе законов сохранения импульса и момента импульса в приближении мгновенного удара и отсутствия деформаций, находятся поступательные и угловые скорости тел после соударения. Реализованная модель была протестирована на задаче из [71] о взаимодействии УВ с числом Маха 5.0 с парой изначально покоящихся цилиндров. В результате цилиндры начинают движение, после чего сталкиваются и разлетаются. На Рис. 4.9 представлены рассчитанные изолинии плотности, а также данные [71].

Давление

Угол, в градусах

Рис. 4.10 - Распределение давления по поверхности цилиндра, угол 0 градусов соответствует точке с координатами (2;0) на Рис. 3.10а, сплошная линия - расчет авторов, пунктирная - расчет [64]. Давление нормировано на величину фонового давления перед УВ. Безразмерный момент времени 1г = 2.4.

Интерес также представляет сравнение результатов расчетов с использованием описанного вычислительного алгоритма с данными, которые получаются при использовании расчетных сеток, согласованных с границами расчетной области. Рассмотрена задача о взаимодействии УВ с числом Маха 1.3 с неподвижным круговым цилиндром радиуса Я = 1.0 из раздела 3.4.3. Расчет проводился на сетке 2 000 х 1 200. Количественное сравнение распределения давления по поверхности цилиндра (см. Рис. 4.10) с результатами расчета из [64] демонстрирует совпадение в пределах погрешности 2%. В [64] расчеты выполнялись в рамках уравнений Навье-Стокса на сетке, согласованной с границами расчетной области, с разрешением пограничного слоя и динамической адаптацией. Поскольку в расчете методом декартовых сеток решение в пересекаемых границами тел ячейках не строится, давление бралось в ближайших к ним внутренних ячейках. Это, в том числе, приводит к наличию высокочастотных пульсаций с амплитудой в доли процента от основного решения на сплошной кри-

вой на Рис. 4.10. Максимальная погрешность наблюдается в направлении 90°, то есть у верхней точки цилиндра, и связана с тем, что в методике отсутствует «подсеточное» разрешение геометрии границы, и решение не строится в пересекаемых границами тел ячейках. Таким образом, в верхней точке цилиндра присутствует одиночная выступающая внешняя ячейка, что заметно на кривой давления.

Наконец, верифицируем методику на основе сравнения с экспериментальными данными. Рассматривается динамика движения уединенной сферической бронзовой частицы диаметра 170 мкм в сверхзвуковом потоке за УВ с числом Маха 2.6. Постановка задачи соответствует натурному эксперименту [98, 99]. Время расчета - 500 мкс. Ошибка в определении координаты частицы в сравнении с экспериментальными данными не превышает 2% (см. Рис. 4.11).

Рис. 4.11 - Динамика движения частиц в сверхзвуковом потоке за УВ. Линии -расчет (сплошная - координата, пунктирная - скорость). Символы - эксперимент [98, 99].

4.4. Результаты вычислительных экспериментов 4.4.1. Движение цилиндров различной массы

Перейдем к задаче о взаимодействии УВ с цилиндрами разной массы. В случае тяжелого цилиндра (т = 1 - 2), его ускорение будет малым, а картина течения на начальном этапе взаимодействия будет мало отличаться от картины течения вокруг неподвижного цилиндра. Давление по поверхности распределено таким образом, что, сила, действующая на цилиндр, возрастает, пока волна не начнет огибать его, после чего монотонно знакопостоянно убывает, см., например, [64]. Такой характер изменения силы приводит к монотонному увеличению скорости цилиндра (см. Рис. 4.12). Рис. 4.12 иллюстрирует динамику изменения скорости цилиндров различной массы. Характер кривых существенно меняется при варьировании массы.

Рис. 4.12 - Динамика изменения скорости цилиндров различной массы (массы подписаны у каждой кривой) от времени. Пунктирная линия - скорость газа за падающей УВ.

В случае легкого цилиндра (т = 0.005) сила давления, которую создаст сжатый УВ газ, сообщит за время огибания волной цилиндра последнему скорость, значительно превосходящую скорость окружающего газа, что создаст область повышенного давления справа от цилиндра и разрежения слева. Сила давления изменит знак, а скорость начет убывать, что с отсрочкой по времени приведет к росту давления слева от цилиндра и разрежению справа. Таким образом, реализуется волнообразное убывание скорости цилиндра до скорости окружающего газа. То есть реализуется колебательный режим изменения силы давления.

1.4 1.3 1.2 1.1 1

0.9 0.8 0.7 0.6

Рис. 4.13 - Ударно-волновая конфигурация при обтекании (а) подвижного (т = 0.2, момент времени 0.23) и (б) неподвижного (момент времени 0.18) цилиндра на начальном этапе. Распределение модуля градиента плотности.

Интерес также представляет ситуация промежуточной массы цилиндра, когда сила давления, достигнув максимума, убывает до нуля, после чего монотонно возрастает (случай т = 0.2 на Рис. 4.12). Опишем картину течения, реализующуюся в данном случае. В работах [42, 100, 101] описаны ударно-волновые конфигурации, возникающие на начальном этапе взаимодействия УВ с подвижным цилиндром, на момент времени, когда фронт УВ еще не опередил

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

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

Рис. 4.14 - Ударно-волновая конфигурация на этапе формирования второй пары тройных точек. Распределение модуля градиента плотности, т = 0.2, момент времени 0.29.

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

на левую - разреженный в результате движения цилиндра. Анализ поля модуля градиентов давления и плотности позволяет наблюдать эффект рождения второй пары тройных точек - точек пересечения волн Маха с новой УВ, разбивающих первую на два участка - СТ2 (от цилиндра до тройной точки Т2) и Т2Т1 (от тройной точки Т2 до тройной точки Т1). Из тройной точки Т2 также выходит контактный разрыв, заканчивающийся на контактном разрыве, выходящим из точки Т1 (см. Рис. 4.14).

Рис. 4.15 - Ударно-волновая конфигурация на этапе удаления друг от друга второй пары тройных точек. Распределение модуля градиента плотности, т = 0.2, момент времени 0.58.

Далее волны Маха встречаются в окрестности точки С. Происходит взаимодействие «верхней» волны Маха с контактным разрывом, выходящим из точки Т2, а также взаимодействие «нижней» волны Маха с контактным разрывом, выходящим из симметрично расположенной тройной точки. Приблизительно когда точки контакта верхней и нижней волн Маха с поверхностью цилиндра достигают соответственно нижней и верхней точки цилиндра (см. Рис. 4.15), сила давления достигает наименьшего значения, а скорость цилиндра - точки перегиба.

4.4.2. Движение двух сфер в сверхзвуковом потоке

Далее представлены результаты по задаче о релаксации пары частиц за проходящей УВ. При решении данной задачи возникла вычислительная сложность, связанная с, так называемым, «эффектом карбункула» [112], который возникал у оси симметрии расчетной области после прохождения УВ второй частицы. Для ее устранения в (4.1), (4.2) был использован численный поток Steger-Warming [103], вместо метода Годунова.

Течение газа за УВ инициирует движение сначала одной, потом другой частицы под действием силы давления газа. Рис. 4.16 иллюстрирует пространственное распределение модуля градиента плотности (так называемая, численная шлирен-визуализация) и несколько мгновенных линий тока. На Рис. 4.16а и Рис. 4.16б, соответствующих моментам времени 30 и 60 мкс изображены основные элементы структуры течения на начальном этапе. Когда первая частица увлекается сверхзвуковым потоком, формируются отошедшая УВ В1, волна, отраженная от оси симметрии Я1 и волна Т1 (см. Рис. 4.16а). Схожая система волн наблюдается и за второй частицей. Лидирующая УВ обозначена через I, волна Маха через М (Маховская конфигурация реализуется за второй частицей). В более поздние моменты времени, конфигурации содержат три косых ударных волны В1, В2 и Я2, а также две волны Т1 и Т2, сходящихся с В1 и В2 в тройных точках Р1 и Р2. Отошедшая от второй частицы волна В2 исчезает приблизительно в момент времени 120 мкс. Через некоторое время происходит столкновение частиц. Все перечисленные наблюдения соответствуют экспериментальным данным [68].

Картины течений на рисунках 30 мкс и 60 мкс качественно соответствуют известному в литературе режиму обтекания двух тел с двумя отошедшими волнами [104]. Картины течений в более поздние моменты времени соответствуют режиму с отрывным течением.

Материал, представленный в Главе 4, опубликован в [20, 21, 23, 34, 35].

Я (мм)

(а) 30 мкс

Я (мм )

(б) 60 мкс

Я ( мм )

(в) 540 мкс

Рис. 4.16 - Численная шлирен-визуализация и несколько мгновенных линий тока

в последовательные моменты времени в задаче о движении пары частиц за УВ.

97

ЗАКЛЮЧЕНИЕ

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

2. Разработан параллельный вычислительный алгоритм метода декартовых сеток второго порядка аппроксимации для математического моделирования распространения ударных волн в областях с изменяющимися границами. Разработан соответствующий программный комплекс. Алгоритм и его программная реализация протестированы на задаче о подъеме цилиндра за проходящей ударной волной и о взаимодействии двух цилиндров за проходящей ударной волной.

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

женных, так и плотных засыпок. Результаты вычислительных экспериментов в терминах избыточного давления за отраженной и прошедшей волнами обобщены с использованием безразмерного критерия в = ( 1 - £ ) • I / ( £ • В ), где £ -проницаемость системы цилиндров, I - ее длина, В - диаметр цилиндров. Для отраженной волны во всех сериях расчетов получено увеличение интенсивности отраженной волны с увеличением в, то есть с увеличением длины засыпки, уменьшением диаметра тел и с уменьшением проницаемости. Относительная погрешность расчетных и доступных экспериментальных данных в диапазоне 2.5 < в < 5.0 находится в пределах 5 - 7%. Расхождение расчетных и экспериментальных данных для случая прошедшей волны существенно больше, чем для отраженной, и достигает десятков процентов при в > 7.5. Наилучшее соответствие между расчетными, экспериментальными и теоретическими результатами наблюдается в том же диапазоне 2.5 < в < 5.0, что и для отраженной волны.

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

СПИСОК СОКРАЩЕНИЙ И УСЛОВНЫХ ОБОЗНАЧЕНИЙ

Латинские символы

с - скорость звука В - диаметр цилиндра Е - полная энергия единицы объема газа е - внутренняя удельная энергия газа / - функция потока в уравнении переноса Н - ширина канала Ь - длина канала I - длина системы тел М - число Маха падающей волны

т - масса; индекс для формулы численного интегрирования в разделе N - число ячеек расчетной области р - давление

£ - проницаемость (пористость) системы тел ^ - время

и - скорость по оси х V - скорость по оси у (х, у) - декартова система координат

Греческие символы

у - показатель адиабаты р - плотность О - безразмерный параметр

Векторы

{ - вектор потока в дифференциальной задаче и - вектор скорости

и - вектор консервативных переменных

Нижние индексы

0 - начальные параметры перед фронтом волны / - номер расчетной ячейки

Верхние индексы

п - номер шага по времени

Иные символы А? - шаг интегрирования по времени

Ах - шаг интегрирования по пространственной координате

Сокращения

УВ - ударная волна

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

1. Ben-Dor G., Igra O., Elperin T. Handbook of Shock Waves. V. 2. - Academic Press. - 2001.

2. Boiko V.M., Kiselev V.P., Kiselev S.P., Papyrin A.N., Poplavsky S.V., Fomin V.M. Shock wave interaction with a cloud of particles // Shock Waves. - 1997. - V. 7. - P. 275 - 285.

3. Mehta Y., Neal C., Salari K., Jackson T.L., Balachandar S., Thakur S. Propagation of a strong shock over a random bed of spherical particles // Journal of Fluid Mechanics. - 2018. - V. 839. - P. 157 - 197.

4. Chaudhuri A., Hadjadj A., Sadot O., Ben-Dor G. Numerical study of shockwave mitigation through matrices of solid obstacles // Shock Wave. - 2013. - V. 23.

- P. 91 - 101.

5. Бедарев И.А., Федоров А.В. Прямое моделирование релаксации нескольких частиц за проходящими ударными волнами // Инженерно-физический журнал. - 2017. - Т. 90, № 2. - С. 450 - 457.

6. Sen O., Davis S., Jacobs G., Udaykumar H.S. Evaluation of convergence behavior of metamodeling techniques for bridging scales in multi-scale multimaterial simulation // Journal of Computational Physics. - 2015. - V. 294. - P. 585 - 604.

7. Бедарев И.А., Федоров А.В., Фомин В.М. Численный анализ течения около системы тел за проходящей ударной волной // Физика горения и взрыва. - 2012.

- Т. 48, № 4. - С. 83 - 92.

8. Бедарев И.А., Федоров А.В. Расчет волновой интерференции и релаксации частиц при прохождении ударной волны // Прикладная механика и техническая физика. - 2015. - Т. 56, № 5. - С. 18 - 29.

9. Regele J.D., Rabinovitch J., Colonius T., Blanquart G. Unsteady effects in dense, high speed, particle laden flows // International Journal of Multiphase Flow. -2014. - V. 61. - P. 1 - 13.

10. Witham G. A new approach to problems of shock dynamics: Part I: Two-dimensional problems // Journal of Fluid Mechanics. - 1957. - V. 2. - P. 145 - 171.

11. Bryson A.E., Gross R.W. Diffraction of strong shocks by cones, cylinders, and spheres // Journal of Fluid Mechanics. - 1961. - V. 10. - P. 1 - 16.

12. Rogg B., Hermann D., Adomeit G. Shock-induced flow in regular arrays of cylinders and packed beds // International Journal of Heat and Mass Transfer. - 1985. -V. 28, No. 12. - P. 2285 - 2298.

13. Saur el R., Abgrall R. A multiphase Godunov method for compressible multifluid and multiphase flows // Journal of Computational Physics. - 1999. - V. 150, No. 2. - P. 425 - 467.

14. Houim R. W., Oran E.S. A multiphase model for compressible granular-gaseous flows: formulation and initial tests // Journal of Fluid Mechanics. - 2016. - V. 789. -P. 166 - 220.

15. Ergun S. Fluid flow through packed columns // Chemical Engineering Progress. - 1952. - V. 48. - P. 89 - 94.

16. Медведев С.П., Фролов С.М., Гельфанд Б.Е. Ослабление ударных волн насадками из гранулированных материалов // Инженерно-физический журнал. -1990. - Т. 58, № 6. - С. 924 - 928.

17. McGrath T.P., Clair J.G.St., Balachandar S. A compressible two-phase model for dispersed particle flows with application from dense to dilute regimes // Journal of Applied Physics. - 2016. - V. 119. - Paper No. 174903.

18. Сидоренко Д.А., Уткин П.С. Метод декартовых сеток для численного моделирования распространения ударных волн в областях сложной формы // Вычислительные методы и программирование. - 2016. - Т. 17, № 4. - С. 353 - 364.

19. Сидоренко Д.А., Уткин П.С. Двумерное газодинамическое моделирование взаимодействия ударной волны с засыпками гранулированных сред // Химическая физика. - 2018. - Т. 37, № 9. - С. 43 - 49.

20. Сидоренко Д.А., Уткин П.С. Численное моделирование релаксации тела за проходящей ударной волной // Математическое моделирование. - 2018. - Т. 30, № 11. - С. 91 - 104.

21. Sidorenko D.A., Utkin P.S. Two-dimensional gas dynamics modeling of the relaxation of particles behind the transmitted shock wave // AIP Conference Proceedings. - 2018. - V. 2027. - Paper 030068. - 6 P.

22. Сидоренко Д.А., Уткин П.С. Комплексный подход к проблеме численного исследования взаимодействия ударной волны с плотным облаком частиц // Горение и взрыв. - 2017. - Т. 10, № 2. - С. 47 - 51.

23. Сидоренко Д.А., Уткин П.С. Численное моделирование взаимодействия ударной волны с подвижным цилиндром // Горение и взрыв. - 2018. - Т. 11, № 3. - С. 79 - 86.

24. Sidorenko D.A., Utkin P.S. Cartesian grid approach for the modeling of shock waves propagation in complex-shape domains // Progress in Detonation Physics / [Edited by S.M. Frolov and G.D. Roy]. - Moscow: Torus Press, 2016. - P. 62 - 70.

25. Sidorenko D.A., Utkin P.S. Unsteady effects in the problem of shock wave -array of cylinders interaction // Nonequilibrium processes in physics and chemistry. Vol. 2. Combustion and detonation / [Edited by A.M. Starik and S.M. Frolov]. -Moscow: Torus Press, 2016. - P. 288 - 295.

26. Sidorenko D., Utkin P. To the complex approach to the numerical investigation of the shock wave - dense particles bed interaction // Proceedings of 31st International Symposium on Shock Waves, 9-14 July 2017, Nagoya, Japan. - Paper SBM000098. - 6 P.

27. Сидоренко Д.А., Уткин П.С. Многомерные эффекты при взаимодействии ударной волны с плотной засыпкой частиц // Сборник докладов V Минского международного коллоквиума по физике ударных волн, горения и детонации, Минск, Беларусь, 13 - 16 ноября 2017 г. - Минск, Беларусь: Ин-т тепло- и мас-сообмена им. А.В. Лыкова. - С. 159 - 165.

28. Sidorenko D.A., Utkin P.S. One Cartesian grids approach for the numerical solution of gas dynamics equations // Waves and Vortices in Complex Media. Proceedings of 4-th International Scientific School of Young Scientists. Moscow, November 26 - 29, 2013. - M.: MAKS Press, 2013. - P. 39 - 41.

29. Сидоренко Д.А., Уткин П.С. Метод декартовых сеток для решения двумерных задач газовой динамики // Научные труды Международной научной конференции XL Гагаринские чтения. Москва, 7 - 11 апреля 2014 г. - Т. 5. - М.: МАТИ, 2014. - С. 174 - 175.

30. Сидоренко Д.А., Уткин П.С. Об одном методе декартовых сеток для решения уравнений газовой динамики // Труды 56-й научной конференции МФТИ. Аэрофизика и космические исследования. Том 2. - Москва - Долгопрудный - Жуковский: МФТИ, 2013. - С. 14 - 15.

31. Сидоренко Д.А., Уткин П.С. Вычислительный алгоритм метода h-box для интегрирования уравнений Эйлера в областях с криволинейными границами // Труды 57-й научной конференции МФТИ. Аэрофизика и космические исследования. - М.: МФТИ, 2014. - С. 138 - 140.

32. Сидоренко Д.А., Уткин П.С. Математическое моделирование взаимодействия ударной волны с клином методом декартовых сеток // Труды 58-ой научной конференции МФТИ. Аэрофизика и космические исследования. Под общ. ред. к.т.н. С.С. Негодяева. - М.: МФТИ, 2015. - С. 167 - 168.

33. Сидоренко Д.А., Уткин П.С. Математическое моделирование взаимодействия ударной волны с системой тел методом декартовых сеток // Труды 59-ой научной конференции МФТИ. Аэрофизика и космические исследования. - М.: МФТИ, 2016. - С. 179 - 180.

34. Сидоренко Д.А., Уткин П.С., Максимов Ф.А. Математическое моделирование взаимодействия ударной волны с подвижным цилиндром // Труды 60-ой научной конференции МФТИ. Аэрокосмические технологии. - М.: МФТИ, 2017. - С. 139 - 141.

35. Sidorenko D.A., Utkin P.S. Two-dimensional gas dynamics modeling of the relaxation of particles behind the transmitted shock wave // Abstracts of International conference on the methods of aerophysical research, August 13 - 19, 2018, Novosibirsk, Russia. Part I. - Novosibirsk: Parallel, 2018. - P. 259 - 260.

36. Barth T., Jespersen D. The design and application of upwind schemes on unstructured meshes // Proceedings of 27th Aerospace Sciences Meeting. USA, Reno. January 9 - 12 1989. - Paper AIAA-89-0366.

37. Barth T., Frederickson P. Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction // Proceedings of 28th Aerospace Sciences Meeting. USA, Reno. January 8 - 11 1990. - Paper AIAA-90-0013.

38. Hu C., Shu C.-W. WENO schemes on triangular meshes // Journal of Computational Physics. - 1999. - V. 150, No. 1. - P. 97 - 127.

39. Mittal R., Iaccarino G. Immersed boundary methods // Annual Review of Fluid Mechanics. - 2005. - V. 37 - P. 239 - 261.

40. Gorsse Y., Iollo A., Telib H., Weynans L. A simple second order Cartesian scheme for compressible Euler flows // Journal of Computational Physics. - 2012. -V. 231, No. 23. - P. 7780 - 7794.

41. Максимов Ф.А. Сверхзвуковое обтекание системы тел // Компьютерные исследования и моделирование. - 2013. - Т. 5, № 6. - С. 969 - 980.

42. Henshaw W.D., Schwendeman D.W. Moving overlapping grids with adaptive mesh refinement for high-speed reactive and non-reactive flow // Journal of Computational Physics. - 2006. - V. 216, No. 2. - P. 744 - 779.

43. Quirk J.J. An alternative to unstructured grids for computing gas dynamic flows around arbitrarily complex two-dimensional bodies // Computers & Fluids. -1994. - V. 23, No. 1. - P. 125 - 142.

44. Zeeuw D., Kenneth G. Powell K.G. An adaptively-refined Cartesian mesh solver for the Euler equations // Journal of Computational Physics. - 1993. - V. 104, No. 1. - P. 56 - 68.

45. Меньшов И.С., Корнев М.А. Метод свободной границы для численного решения уравнений газовой динамики в областях с изменяющейся геометрией // Математическое моделирование. - 2014. - Т. 26, № 5. - С. 99 - 12.

46. Меньшов И.С., Павлухин П.В. Эффективный параллельный метод сквозного счета задач аэродинамики на несвязанных декартовых сетках // Журнал вычислительной математики и математической физики. - 2016. - Т. 56, № 9. -С.1677 - 1691.

47. Жлуктов С.В., Аксенов А.А., Карасев П.И. Моделирование байпасного ламинарно-турбулентного перехода в рамках k-e подхода // Компьютерные исследования и моделирование. - 2014. - Т. 6, № 6. - С. 879 - 888.

48. Berger M., Helzel C. A simplified h-box method for embedded boundary grids // Journal of Scientific Computing. - 2012. - V. 34, No. 2. - P. 861 - 888.

49. William J.C., Powell K.G. An accuracy assessment of Cartesian-mesh approaches for the Euler equations // Journal of Computational Physics. - 1995. - V. 117, No. 1. - P. 121 - 131.

50. Chiang Y., Bram V.L., Powell K.G. Simulation of unsteady inviscid flow on an adaptively refined Cartesian grid // Proceedings of 30th Aerospace Sciences Meeting and Exhibit. USA, Reno. - January 6 - 9 1992. - Paper AIAA-92-443.

51. Yang G., Causon D., Ingram D., Saunders R., Battent P. A Cartesian cut cell method for compressible flows. Part B: moving body problems // The Aeronautical Journal. - 1997. - V. 101, Issue 1002. - P. 57 - 65.

52. Ji H., Lien F., Yee E. Numerical simulation of detonation using an adaptive Cartesian cut-cell method combined with a cell-merging technique // Computers & Fluids. - 2010. - V. 39, No. 6. - P. 1041 - 1057.

53. Yang G., Causon. D., Ingram D., Saunders R., Battent P. A Cartesian cut cell method for compressible flows Part A: Static body problems // The Aeronautical Journal. - 1997 - V. 101, No. 1002. - P. 47 - 56.

54. May S., Berger M. Two-dimensional slope limiters for finite volume schemes on non-coordinate-aligned meshes // Journal on Scientific Computing. - 2013 - V. 35, No. 5. - P. 2163 - 2187.

55. Harten A., Lax P.D., Van Leer. B On upstream differencing and Godunov-type schemes for hyperbolic conservation laws // SIAM Review. - 1983 - V. 25, No. 1. -P. 35 - 61.

56. Colella P., Graves D.T., Keen B.J., Modiano D. A Cartesian grid embedded boundary method for hyperbolic conservation laws // Journal of Computational Physics. - 2006 -V. 211, No. 1. - P. 347 - 366.

57. Pember R.B., Bell J.B., Colella P., Crutchfield W.Y., Welcome M.L. An adaptive Cartesian grid method for unsteady compressible flow in irregular regions // Journal of Computational Physics. - 1995. - V. 120, No. 2. - P. 278 - 304.

58. Hu X.Y., Khoo B.C., Adams N.A., Huang F.L. A conservative interface method for compressible flows // Journal of Computational Physics. - 2006. - V. 219, No. 2. - P. 553 - 578.

59. Forrer H., Jeltsch R. A higher-order boundary treatment for Cartesian-grid methods // Journal of Computational Physics. - 1998. - V. 140, No. 2. - P. 259 -277.

60. Forrer H., Berger M. Flow simulations on Cartesian grids involving complex moving geometries // Proc. of 7th International Conference «Hyperbolic Problems: Theory, Numerics, Applications». Switzerland, Zurich, 1999. - V. 1. - P. 315 - 324.

61. Годунов С.К., Забродин А.В., Иванов М.Я., Крайко А.Н., Прокопов Г.П. Численное решение многомерных задач газовой динамики. - М.: Наука, 1976.

62. Gueyffier D., Li J.,Nadim A., Scardovelli R., Zaleski S. Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows // Journal of Computational Physics. - 1999. - V. 152, No. 2. - P. 423 - 456.

63. Forrest E.G., Edward W.P. Drag of circular cylinders of a wide range of Reynolds numbers and Mach numbers // National Advisory Committee for Aeronautics. Technical Note 2960. - USA, Washington, 1953.

64. Drikakis D., Ofengeim D., Timofeev E., Voionovich P. Computation of non-stationary shock wave/cylinder interaction using adaptive-grid methods // Journal of Fluids and Structures. - 1997. - V. 11, No. 6. - P. 665 - 692.

65. Igra O., Takayama K. Shock tube study of the drag coefficient of a sphere in a non-stationary flow // Proceedings of the Royal Society of London A. - 1993. - V. 442, No. 1915. - P. 231 - 247.

66. Jourdan G., Houas L., Igra O., Estivalezes J.-L., Devals C., Meshkov E.E. Drag coefficient of a sphere in a non-stationary flow: new results // Proceedings of the Royal Society of London A. - 2007. - V. 463. - P. 3323 - 3345.

67. Бойко В.М., Клинков К.В., Поплавский С.В. Коллективный головной скачок перед поперечной системой сфер в сверхзвуковом потоке за проходящей ударной волной // Известия Российской академии наук. Механика жидкости и газа. - 2004. - № 2. - С. 183.

68. Boiko V.M., Klinkov K. V., Poplavski S. V. On a mechanism of interphase interaction in non-relaxing two-phase flow // Proceedings of 11th International Conference on Methods of Aerophysical Research. Novosibirsk, Russia, 1 - 7 July 2002. -P. 24 - 27.

69. Фролов С.М. Эффективность ослабления ударных волн в каналах различными способами // Журнал прикладной механики и технической физики. -1993. - № 1. - С. 5 - 39.

70. Abe A., Takayama K., Itoh K. Experimental and numerical study of shock wave propagation over cylinders and spheres // Transactions on Modelling and Simulation. - 2001. - V. 30 - P. 209 - 218.

71. Nourgaliev R.R., Dinh T.N., Theofanous T.G., Koning J.M., Greenman R.M., Nakafuji G.T. Direct numerical simulation of disperse multiphase high-speed flows // Proceedings of 42nd AIAA Aerospace Sci. Meet. and Exhibit, Reno, Nevada, USA, January 5 - 8 2004, AIAA paper 2004-1284.

72. Chaudhuri A., Hadjadj A., Sadot O., Glazer E. Computational study of shockwave interaction with solid obstacles using immersed boundary methods // Interna-

tional Journal for Numerical Methods in Engineering. - 2012. - V. 89. - P. 975 -990.

73. Wan Q., Eliasson V. Numerical study of shock wave attenuation in two-dimensional ducts using solid obstacles: how to utilize shock focusing techniques to attenuate shock waves // Aerospace. - 2015. - V. 2. - P. 203 - 221.

74. Wagner, J.L., Beresh, S.J., Kearney, S.P. A multiphase shock tube for shock wave interactions with dense particle fields // Experiments in Fluids. - 2012 - V. 5 2(6) - P. 1507 - 1517.

75. Clemins A. Representation of two-phase flows by volume averaging // International Journal of Multiphase Flow. - 1988 - V. 14, No. 1. - P. 81 - 90.

76. Metha Y., Neal C., Jackson T., Balachandar S., Thakur S. Shock interaction with three-dimensional face centered cubic array of particles // Physical Review Fluids. - 2016 - V. 1 - 054202.

77. Barth T., Ohlberger M. Finite Volume Methods: Foundation and Analysis // Encyclopedia of Computational Mechanics. - 2004 - V. 1 - P. 439 - 470.

78. Toro E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction - Springer, 2009.

79. BergerM.J., Helzel C., LeVeque R. J. H-box methods for the approximation of hyperbolic conservation laws on irregular grids // SIAM Journal on Numerical Analysis. - 2003. - V. 41, No. 3. - P. 893 - 918.

80. Sutherland I.E., Hodgman G.W. Reentrant polygon clipping // Communications of the ACM. - 1974. - V. 17, No. 1. - P. 32 - 42.

81. Levy D.W., Powell K.G., van Leer B. Use of a rotated Riemann solver for the two-dimensional Euler equations // Journal of Computational Physics. - 1993. - V. 106, No. 2. - P. 201 - 214.

82. Ren Y.X. A robust shock-capturing scheme based on rotated Riemann solvers // Computers & Fluids. - 2003. - V. 32, No. 10. - P. 1379 - 1403.

83. Glaz H.M., Colella P., Glass I.I., Deschambault R.L. A numerical study of oblique shock-wave reflections with experimental comparisons // Proceedings of the Royal Society of London A. - 1985. - V. 398, No. 1814. - P. 117 - 140.

84. Баженова Т.В., Гвоздева Л.Г. Нестационарное взаимодействие ударных волн. - М.: Наука, 1977. - 274 С.

85. Takayama K., Itoh K. Unsteady Drag over Cylinders and Aerofoils in Transonic Shock Tube Flows // Proc. 15th International Symposium on Shock Waves and Shock Tubes. - 1985 - P. 439 - 485.

86. Crowe C. T. Multiphase Flow Handbook. - CRC Press. - 2005.

87. Уткин П.С. Математическое моделирование взаимодействия ударной волны с плотной засыпкой частиц в рамках двухжидкостного подхода // Химическая физика. - 2017. - Т. 36, № 11. - С. 61 - 71.

88. Tanino Y., Nepf H.M. Laboratory investigation of mean drag in a random array of rigid, emergent cylinders // Journal of Hydraulic Engineering. - 2008. - V. 134, No. 1. - P. 34 - 41.

89. Абалкин И.В., Жданова Н.С., Козубская Т.К. Реализация метода погруженных границ для моделирования задач внешнего обтекания на неструктурированных сетках // Математическое моделирование. - 2015. - Т. 27, № 10. - С. 5 - 20.

90. Колган В.П. Применение принципа минимальных значений производной к построению конечноразностных схем для расчета разрывных решений газовой динамики // Ученые записки ЦАГИ. - 1972. - Т. 3, № 6. - С. 68 - 77.

91. Chertock A., Kurganov A. A simple Eulerian finite-volume method for compressible fluids in domains with moving boundaries // Communications in Mathematical Sciences. - 2008. - V. 6, No. 3. - P. 531 - 556.

92. Sambasivan S.K., Udaykumar H.S. Ghost fluid method for strong shock interactions. Part 2: Immersed solid boundaries // AIAA Journal. - 2009. - V. 47, No. 12. - P.2923 - 2937.

93. Arienti M., Hung P., Morano E., Shepherd J.E. A level set approach to Eulerian - Lagrangian coupling // Journal of Computational Physics. - 2003. - V. 185, No. 1. - P. 213 - 251.

94. Tan S., Shu C.-W. A high order moving boundary treatment for compressible inviscid flows // Journal of Computational Physics. - 2011. - V. 230, No. 15. - P. 6023 - 6036.

95. Shyue K.M. A moving-boundary tracking algorithm for inviscid compressible flow // Proceedings of 11th International Conference «Hyperbolic Problems: Theory, Numerics, Applications». France, Lyon, 2006. - Springer, 2008. - P. 989 - 996.

96. Das P., Sen. O., Jacobs, G., Udaykumar H.S. A sharp interface Cartesian grid method for viscous simulation of shocked particle-laden flows // International Journal of Computational Fluid Dynamics. - 2017. - V. 31. - P. 269 - 291.

97. Гольдсмит. В. Удар. Теория и физические свойства соударяемых тел. -М.: Издательство литературы по строительству. - 1965 - 448 с.

98. Boiko V.M., Fedorov A.V., Fomin V.M., Papyrin A.N., Soloukhin R.I. Ignition of small particles behind shock waves // Shock Waves, Explosions and Detonations. Progress in Astronautics and Aeronautics. AIAA. - 1983. - V. 87. - P. 71 - 87.

99. Бойко В.М. Исследование динамики ускорения, разрушения и воспламенения частиц за ударными волнами методами лазерной визуализации // Диссертация к.ф.-м.н. - Новосибирск, 1984. - 205 с.

100. Muralidharan B., Menon S. Simulations of unsteady shocks and detonation interactions with structures // Proc. 49th AIAA/ASME/SAE/ASEE Joint Propulsion Conference. - 2013. - AIAA Paper 3655.

101. Luo K., Luo Y., Jin T, Fan J. Studies on shock interactions with moving cylinders using immersed boundary method // Physical Review Fluids. - 2017 - V. 2. -Paper 064302.

102. Quirk J.J. A contribution to the great Riemann solver debate // International Journal for Numerical Methods in Fluids. - 1994. - V. 18. - P. 555 - 574.

103. Steger J.L., Warming R.F. Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods // Journal of Computational Physics. - 1981. - V. 40. - P. 263 - 293.

104. ChangP.K. Separation of Flow. - London: Pergamon Press. - 1970. - 760 p.