автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Моделирование взаимодействий ударных волн с использованием неструктурированных расчётных сеток

На правах^рукописи

005531838

Эпштейн Дмитрий Борисович

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

05.13.18 - Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук - 1 ДВГ

тг ш I -

Новосибирск - 2013

005531838

Работа выполнена в Федеральном государственном бюджетном учреждении науки Институте теоретической и прикладной механики им. С.А. Христиано-вича Сибирского отделения Российской академии наук

Научный руководитель: кандидат физико-математических наук

Кудрявцев Алексей Николаевич

Официальные оппоненты: Ковеня Виктор Михайлович,

доктор физико-математических наук, профессор,

ИВТ СО РАН, главный научный сотрудник

Зудов Владимир Николаевич, доктор физико-математических наук, ИТПМ СО РАН, ведущий научный сотрудник

Ведущая организация: Федеральное государственное бюджетное

образовательное учреждение высшего профессионального образования «Балтийский государственный технический университет «ВОЕНМЕХ» им. Д.Ф. Устинова»

Защита состоится 23 октября 2013 г. в 15:30 на заседании диссертационного совета ДООЗ .046.01 в Институте вычислительных технологий СО РАН по адресу: 630090, г. Новосибирск, пр. Ак. Лаврентьева 6.

С диссертацией можно ознакомиться в специализированном читальном зале вычислительной математики и информатики ГПНТБ СО РАН.

Автореферат разослан июля 2013 г.

Ученый секретарь диссертационного совета, к.ф.-м.н., доцент

А.С. Лебедев

Общая характеристика работы

Актуальность темы. Численное решение уравнений газовой динамики является в настоящее время одним из наиболее разработанных разделов вычислительной механики сплошных сред. Различные подходы к численному моделированию газодинамических течений описаны в целом ряде известных монографий (Годунов С.К. и др, 1976, Ковеня В.М., Яненко H.H., 1981, Куликовский А.Г. и др., 2001). Быстрый прогресс в области вычислительной техники привел к тому, что в наши дни стал возможным расчёт весьма сложных трехмерных нестационарных течений. Это в полной мере относится и к задачам, связанным с распространением и взаимодействием ударных волн. Современные схемы сквозного счета позволяют с высокой точностью рассчитывать сверхзвуковые течения газа, избегая при этом возникновения нефизических осцилляций вблизи газодинамических разрывов. Это открывает широкие возможности для применения численного моделирования как при решении научных задач, так и в различных практических приложениях.

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

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

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

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

Цели диссертационной работы. Целями данной работы, являются:

— Разработка численного алгоритма и расчетного кода, основанного на использовании неструктурированных расчётных сеток и позволяющего проводить моделирование на многопроцессорных ЭВМ двумерных и трёхмерных течений сжимаемого газа;

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

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

1. Алгоритм решения уравнений газовой динамики на неструктурированных расчётных сетках на основе конечно-объёмной TVD (Total Variation Diminishing) схемы 2-го порядка точности;

2. Расчетный код, реализующий численный алгоритм, и способ его паралле-лизации с помощью метода декомпозиции расчётной области и библиотеки MPI;

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

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

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

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

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

— Международный симпозиум по ударным волнам (ISSW), Гёттинген, Германия, 2007, Санкт-Петербург, 2009, Манчестер, Великобритания, 2011

— Международная студенческая конференция: «Студент и научно-технический прогресс», НГУ, Новосибирск, 2007;

— Всероссийская конференция молодых учёных «Проблемы механики: теория, эксперимент и новые технологии», Новосибирск, 2007

— Международная конференция по методам аэрофизических исследований (ICMAR), Новосибирск, 2007, 2008, 2010;

— Международный симпозиум по взаимодействию ударных волн (ISIS), Ру-ан, Франция, 2008

— Международная конференция по военным аспектам ударных и взрывных волн (MABS), Иерусалим, Израиль, 2010.

— Семинар «Математическое моделирование в механике» ИТПМ им. С.А. Христиановича СО РАН (Новосибирск) под руководством академика В.М. Фомина и проф. A.B. Федорова.

— Семинар «Информационно-вычислительные технологии» ИВТ СО РАН (Новосибирск) под руководством академика Ю.И. Шокина и проф. В.М. Ковени.

Публикации. По материалам диссертации опубликовано 10 работ, в том числе 2 статьи в рецензируемых журналах, рекомендованных ВАК для для представления основных результатов диссертации. Список работ приведен в конце автореферата.

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

Структура и объем диссертации. Диссертационная работа состоит из введения, четырех глав, заключения и списка литературы из 101 наименования. Общий объем диссертации составляет 113 страниц, включая 45 рисунков и одну таблицу.

Основное содержание работы

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

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

<9(3 д¥ дО + Ж = 0 дг дх ду дг где <2 — вектор консервативных переменных:

т

(} = (р,ри,ру,ры,Е) ,

с, Н — векторы потоков соответственно вдоль осей х, у и г:

Е= (ри,ри2 + р,риу,рию,и(Е + р))Т, в = (ру,рм,ру2 + р,руы,у(Е + р)) , Н = (ри>, рии>, рун>, ры2 + р, \\>(Е + р)) ,

здесь р — плотность газа, и, V, уу — компоненты вектора скорости, р — давление газа, Е — полная энергия на единицу объёма.

д_ д1

///<2 *** = -III (Тх + д£ + д1) =

п п

- £(¥пх + Спу + Нлг) ¿8, (2)

где - контрольный объем, ограниченный замкнутой поверхностью Г = д£1, пх, пу и пг - составляющие вектора внешней нормали к поверхности Г, с13 — элемент площади поверхности. Раздел 1.2 посвящен обсуждению метода сквозного счета для расчёта течений, включающих газодинамические разрывы. В данной работе используется одна из схем сквозного счёта — конечнообъемная МШСЬ ТУЕ) схема второго порядка точности. Искомыми величинами являются значения газодинамических переменных, осреднённые по объемам ячеек сетки, состоящей из треугольников (в двумерном случае) или тетраэдров (в трёхмерном). Процедура численного решения состоит из трех последовательных этапов:

1. Реконструкция переменных на гранях, разделяющих ячейки.

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

3. Интегрирование по времени и нахождение значения газодинамических переменных на новом временном слое.

Эти этапы подробно обсуждаются в следующих трех разделах. Раздел 1.3 посвящен процедуре реконструкции переменных. При реконструкции первого порядка решение считается постоянным внутри каждой ячейки и соответствующие значения присваиваются и переменным в центре ограничивающих ячейку граней. Реконструкция второго порядка основана на кусочно-линейном представлении решения. Линейная функция строится исходя из величины градиента, вычисленного из средних по объему самой ячейки и соседних ячеек (двух в двумерном случае, см. рис. 1, трёх — в трёхмерном). В соответствии с построенным кусочно-линейным представлением вычисляются переменные в центрах граней (в двумерном случае — сторон). Для того, чтобы избежать возникновения численных осцилляций вблизи газодинамических разрывов и получить схему, удовлетворяющую локальному принципу максимума, что исключает появление новых локальных экстремумов при реконструкции, расчитанные значения ограничиваются с помощью нелинейного ограничителя (лимитера) тттос!.

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

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

метод Роу, метод Хартена-Лакса-ван Леера (HLL) и его модификации, метод расщепления потоков по ван Лееру.

После нахождения потоков уравнения интегрируются по времени и находятся значения переменных на новом временном слое (раздел 1.5). Для интегрирования использовалась двухстадийная схема Рунге-Кутта второго порядка точности. Шаг по времени выбирается так, чтобы удовлетворять условию Куранта-Фридрихса-Леви.

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

Раздел 2.1 посвящен вопросам, связанным с построением неструктурированных расчетных сеток. Для построения используются свободно распространяемые генераторы сеток: Triangle (в двумерном случае), Tetgen и Netgen (в трехмерном). Для задания геометрии области и преобразования данных, выдаваемых генераторами сеток, в формат, используемый расчётной программой, разработан ряд вспомогательных программ-утилит. Строящиеся сетки могут быть легко сгущены в областях, где требуется более высокое пространственное разрешение. На рис. 2 показан пример такого локального сгущения сетки при моделирования маховского отражения ударной волны.

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

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

Процедура параллелизации обсуждается в разделе 2.4. Производится геометрическая декомпозиция расчётной области — разбиение ее на подобласти, присваиваемые отдельным процессорам. Для декомпозиции используется про-

Рис. 2: Пример сгущения расчётной неструктурированной сетки.

грамма METIS, стремящаяся минимизировать объем данных, которые необходимо передавать между процессорами. Для самого обмена данных используется библиотека MPI.

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

Как видно из рис. 3, эффективность масштабированного расчёта выше — около 93%, в то время как при немасштабированном расчёте эффективность параллелизации составила около 76%.

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

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

Рис. 3: Немасштабированное (линия 1) и масштабированное (линия 2) параллельное ускорение при проведении расчётов на нескольких процессорах при различных значениях загрузки процессоров. Сплошная линия показывает идеальное ускорение.

клиньями, помещенными в сверхзвуковой поток. Эта задача, в которой наблюдается достаточно тонкий эффект гистерезиса перехода, была ранее подробно изучена М.С. Ивановым и соавторами (2001) с помощью кода на структурированной сетке. Результаты наших расчётов на неструктурированной сетке полностью согласуются с их данными.

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

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

Далее исследуется, как меняются возникающие ударно-волновые конфигурации с изменением числа Маха набегающего потока М (раздел 3.2) или относительного расстояния между цилиндрами H/R (раздел 3.3, R — радиус цилиндров, Я — расстояние между их центрами). При изменении параметров задачи каждый новый расчет стартовал со стационарного поля, полученного при предыдущем значении параметра. Впервые показано, что в обоих рассмотренных случаях наблюдается гистерезис перехода между регулярным и маховским режимами взаимодействия головных ударных волн соседних цилиндров. Именно, переход к маховскому отражению при уменьшении М (или H/R) и обратный переход при увеличении М (соответственно, H/R) происходят при разных значениях числа Маха (соответственно, H/R). Это иллюстрируется рис. 4, на котором показаны поля течения (численные шлирен-визуализации) полученные при изменении относительного расстояния. Для хорошего разрешения ножки Маха малых размеров использовалось локальное сгущение сетки, при котором мини-

мальный линейный размер ячейки достигал 1/40 от размера ячейки основной сетки, пример сгущения был показан на рис. 2

На рис. 5а и 56 показано как изменяется высота ножки Маха при изменении соответственно числа Маха набегающего потока и относительного расстояния между цилиндрами. Очевиден гистерезисный характер зависимости.

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

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

(а) (б)

Рис. 5: Зависимость высоты ножки Маха от числа Маха набегающего потока при H/R = 5 (а) и относительного расстояния между цилиндрами при М = 4 (б).

ных данных выбирался равномерный поток, всегда наблюдалось регулярное отражение, в экспериментах — маховское (последнее объясняется влиянием присутствующих в экспериментальной установке возмущений набегающего потока, способных при ccj < а < ад? инициировать переход к маховскому отражению).

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

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

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

В разделе 4.2 прохождение ударной волны через систему цилиндров моделировалось в точности в той постановке, как оно изучалось в эксперименте Suzuki et al., 2000, выполненном в ударной трубе. Расчёты были выполнены для

регулярного и «шахматного» расположения цилиндров. Характерная картина течения в момент времени, когда уже сформировались единые отраженная и прошедшая ударные волны, показана на рис. 6.

1S 10 Б 0 -5 -10

0 -30 у ' ; -20 -10 0 х 10 20 ЗО 40 і . BO

Рис. 6: Численная шлирен-визуализация поля течения в момент времени t = 60 при прохождении ударной волны через систему цилиндров при их регулярном расположении.

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

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

В разделе 4.4 приведено исследование взаимодействия взрывной волны с системой цилиндров, образующейся при выделении большого количества энергии в цилиндрической зоне малого радиуса, расположенной на некотором расстоянии от барьера. На рис. 8а показана сложная картина течения, возникающего при нестационарном взаимодействии взрывной волны с проницаемой преградой. Распределения давления вдоль оси X при различных значениях H/R приведены на рис. 86 (они даны для момента, когда прошедшая ударная волна достигает точки X = 39). При уменьшении H/R ослабевает интенсивность прошедшей через барьер ударной волны и, соответственно, уменьшается скорость ее фронта. Однако ослабляющее действие барьера становится практически неощутимым при H/R > 5.

Заключение содержит основные результаты и выводы по работе:

1. Разработан и верифицирован численный метод на основе конечно-объёмной TVD (Total Variation Diminishing) схемы 2-го порядка точности, позволяющей проводить сквозной счет течений, содержащих газодинамиче-

Рис. 7: Пример построенной сетки (1 458 239 ячейки) в случае H/R = 3 (а) и изоповерхности давления при взаимодействии ударной волны с системой сфер при H/R = 2,5 и М = 3,0 в момент времени t = 15(6).

(а) (б)

Рис. 8: Шлирен-визуализация взаимодействия взрывной волны с системой цилиндров, Я//? = 3, г = 18,8 (а) и распределения давления вдоль оси X для различных Я/Я в момент времени, когда прошедшая ударная волна достигает точки X = 39 (б).

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

2. Произведена эффективная параллелизация разработанного алгоритма для расчётов на многопроцессорных ЭВМ. Для параллелизации в работе использован метод геометрической декомпозиции — разбиение расчётной области на подобласти. Обмен данными между процессорами производился с помощью библиотеки MPI. В тестовых случаях при использовании

64 процессоров эффективность масштабированного ускорения составила около 93%, а немасштабированного — около 76%.

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

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

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

Автор выражает глубокую признательность своему научному руководителю А.Н. Кудрявцеву за его всестороннюю поддержку. Также автор хотел бы поблагодарить Д.В. Хотяновского за неоднократные полезные обсуждения полученных результатов.

Список основных работ по теме диссертации

Публикации в периодических изданиях, рекомендованных ВАК:

1. Kudryavtsev, A.N. Hysteresis phenomenon at interaction of shock waves generated by a cylinder array. / A.N. Kudryavtsev, D.B Epstein // Shock Waves. — 2012. — V. 22, № 4. — P. 341-349.

2. Кудрявцев, А.Н. Явление гистерезиса при обтекании системы цилиндров сверхзвуковым потоком. / А.Н. Кудрявцев, Д.Б. Эпштейн // Изв. РАН. МЖГ. — 2012. — № 3. — С. 122-131.

Публикации в трудах международных и всеросскийских конференций:

3. Kudryavtsev, A.N. Investigation of interaction between shock waves and flow disturbances with different shock-capturing schemes / A.N. Kudryavtsev, D.V. Khotyanovsky, D.B. Epstein // Proceedings of 26th International Symposium on Shock Waves. — Springer: Berlin, Heidelberg, 2009. — Vol. 2 —P. 1023-1028.

4. Эпштейн, Д.Б. Моделирование взаимодействий ударных волн с использованием неструктурированных расчётных сеток / Д.Б. Эпштейн // Труды XLV международной студенческой конференции «Студент и научно-технический прогресс» — 2007. — С. 26-31

5. Эпштейн, Д.Б. Моделирование взаимодействий ударных волн с использованием неструктурированных расчётных сеток / Д.Б. Эпштейн, А.Н. Кудрявцев // IV Всероссийская конференция молодых учёных «Проблемы механики: теория, эксперимент и новые технологии». — 2007. — С. 67-68.

6. Kudryavtsev, A.N. Numerical simulations of interaction of shock waves with small disturbances / A.N. Kudryavtsev, D.V. Khotyanovsky, D.B. Epstein, A.Yu. Ovsyannikov // Proceedings of 18th International Shock Interaction Symposium. — 2008. — P. 87-90

7. Kudryavtsev, A.N. Interaction of steady and unsteady shock waves generated by a cylinder array / A.N. Kudryavtsev, D.B. Epstein // Proceedings of 18th International Shock Interaction Symposium. — 2008. — P. 135-137.

8. Epstein, D.B. Numerical simulation of shock-wave interaction with a system of cylinders / D.B. Epstein // XIV International Conference on the Method of Aerophysical Research. — 2008. — P. 11-12

9. Epstein, D.B. Hysteresis phenomenon at the interaction of bow shock waves from cylinder bodies /D.B. Epstein, A.N. Kudryavtsev // Proceedings of 27th International Symposium on Shock Waves — 2009.— Paper № 30576, — P. 371.

10. Kudryavtsev, A.N. Numerical simulation of shock and blast wave propagation over a porous barrier / A.N. Kudryavtsev, D.B. Epstein // XXI International Symposium on Military Aspects of Blast and Shock — 2010. — P. 25

11. Epstein, D.B. Hysteresis phenomenon at aerodynamic interaction of cylinder bodies in supersonic flow / D.B. Epstein, A.N. Kudryavtsev // XV International Conference on the Method of Aerophysical Research. — 2010. — P. 57-58

12. Epstein, D. Shock and blast propagating throw a porous barrier / D. Epstein // 28th International Symposium on Shock Waves — 2011. — Paper№ 2551, 6 p.

Подписано в печать 05.07.2013 Тираж 135 экз., Заказ № 038 Отпечатано «Документ-Сервис» 630090, Новосибирск, ул. Институтская, 4/1, тел. (383) 335-66-00

Федеральное государственное бюджетное учреждение науки Институт теоретической и прикладной механики им. С.А. Христиановича Сибирского отделения Российской академии наук

На правах рукописи

04201364221

Эпштейн Дмитрий Борисович

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

РАСЧЁТНЫХ СЕТОК

05.13.18 - Математическое моделирование, численные методы

и комплексы программ

Диссертация на соискание ученой степени кандидата физико-математических наук

Научный руководитель

кандидат физико-математических наук

Кудрявцев Алексей Николаевич

Новосибирск - 2013

Содержание

Введение............................................................................4

Глава 1. Численный алгоритм решения уравнений газовой динамики на

неструктурированных расчетных сетках..................................10

1.1. Основные уравнения....................................................10

1.2. Метод сквозного счёта для решения уравнений газовой динамики . 12

1.3. Реконструкция решения на гранях ячеек..............................17

1.4. Нахождение потоков через грани......................................21

1.5. Интегрирование по времени.......................33

Глава 2. Программная реализация алгоритма для многопроцессорных ЭВМ, верификация разработанного расчетного кода..........35

2.1. Построение неструктурированной сетки...............35

2.2. Структуры данных...........................39

2.3. Граничные условия...........................42

2.4. Параллелизация.............................43

2.5. Тестовые расчёты............................46

Глава 3. Взаимодействие ударных волн при обтекании системы цилиндрических тел сверхзвуковым потоком..................62

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

3.2. Гистерезис при изменении числа Маха набегающего потока.....66

3.3. Гистерезис при изменении относительного расстояния между цилиндрами .................................70

3.4. Образование коллективной головной ударной волны ........73

Глава 4. Численное моделирование прохождения ударных и взрывных волн через системы цилиндров и сфер ..................

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

4.2. Численное моделирование эксперимента в ударной трубе .....83

4.3. Взаимодействие ударной волны с системой сфер...........86

4.4. Взаимодействие взрывной волны с системой цилиндров.......92

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

Литература....................................101

Введение

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

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

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

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

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

Разумеется, за данные преимущества приходится и платить. Необходимы специальные алгоритмы, позволяющие решать уравнения газовой динамики на таких сетках. Вообще говоря, при этом обычно приходится ограничиться численными схемами второго порядка точности, поскольку разработка методов более высокого порядка на неструктурированных сетках оказывается намного более сложной задачей, чем на структурированных (отметим, однако, что в последние годы в этом направлении были достигнуты важные успехи, связанные с развитием т. н. разрывных методов Галеркина — см. [10]). Далее, структуры данных, необходимых для программной реализации алгоритмов решения уравнений с частными производными на неструктурированных сетках, оказываются заметно более громоздкими и требуют для своего хранения большего объема машинной памяти.

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

Неудивительно поэтому, что разработка таких подходов для решения уравнений газовой динамики началась достаточно давно. Детальный обзор зарубежных работ в этом направлении содержится в [11, 12]. На русском языке следует назвать обзорную статью [13], содержащую подробную библиографию как российских, так и зарубежных публикаций. Отметим, что несмотря на наличие целого ряда важных работ по применению неструктурированных сеток в вычислительной газовой динамики, выполненных отечественными авторами (см. [14]-[20]), в целом методы, основанные на таких сетках, используются у нас недостаточно широко. Это особенно относится к схемам для решения трехмерных уравнений Эйлера. Возможно, это является одной из причин, затрудняющих развитие отечественных универсальных вычислительных пакетов.

Расчет достаточно сложных трехмерных течений требует использования многопроцессорных ЭВМ. Это означает, что расчетный код должен быть распараллелен. Параллелизация расчетных программ, использующих неструктурированные сетки, заметно сложнее, чем кодов, основанных на применении структурированных сеток. Из отечественных работ, посвященных этой теме, можно указать разве что работу [20].

Целями настоящей диссертационной работы являются:

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

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

Диссертация состоит из введения, четырех глав, заключения и списка литературы. Она содержит 45 рисунков и одну таблицу. Объем работы 115 страниц.

Библиографический список включает 101 наименование.

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

Первая глава посвящена изложению численного алгоритма решения уравнений газовой динамики на неструктурированных расчетных сетках. В ней формулируются решаемые уравнения и описываются основные составляющие части используемой нами TVD (Total Variation Diminishing) схемы сквозного счета: реконструкция со вторым порядком газодинамических переменных на гранях расчетных ячеек из известных величин, средних по объемам ячеек, нахождение потоков через грани путем приближенного решения задачи о распаде разрыва и интегрирование уравнений по времени.

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

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

В четвертой главе описываются результаты численного моделирования про-

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

Заключение содержит основные результаты и выводы по работе.

Результаты, представленные в диссертации, докладывались на следующих семинарах и конференциях:

— Международный симпозиум по ударным волнам (ISSW), Гёттинген, Германия, 2007, Санкт-Петербург, 2009, Манчестер, Великобритания, 2011;

— Международная студенческая конференция: «Студент и научно-технический прогресс», НГУ, Новосибирск, 2007;

— Всероссийская конференция молодых учёных «Проблемы механики: теория, эксперимент и новые технологии», Новосибирск, 2007;

— Международная конференция по методам аэрофизических исследований (ICMAR), Новосибирск, 2007, 2008, 2010;

— Международный симпозиум по взаимодействию ударных волн (ISIS), Руан, Франция, 2008;

— Международная конференция по военным аспектам ударных и взрывных волн (MABS), Иерусалим, Израиль, 2010;

— Семинар «Математическое моделирование в механике» ИТПМ им. С.А. Христиановича СО РАН (Новосибирск) под руководством академика В.М. Фомина и проф. A.B. Федорова;

— Семинар «Информационно-вычислительные технологии» ИВТ СО РАН (Новосибирск) под руководством академика Ю.И. Шокина и проф. В.М. Кове-ни.

По теме диссертации опубликованы работы [21]-[32].

Глава 1

Численный алгоритм решения уравнений газовой динамики на неструктурированных расчетных

сетках

1.1. Основные уравнения

Настоящая работа посвящена численному моделированию взаимодействия ударных волн в областях сложной геометрической формы. Рассматриваемые явления описываются нестационарными уравнениями Эйлера. Полная система уравнений Эйлера для сжимаемого невязкого и нетеплопроводного газа записывается в безразмерном виде следующим образом [33].

Уравнение сохранения массы:

где г — время, х/ — компоненты вектора пространственных координат, р — плотность газа, и] — компоненты вектора скорости. Уравнение сохранения импульса:

Здесь Е — полная энергия на единицу объёма, которая равняется сумме внутренней и кинетической энергий, т.е.

(1.1)

(1.2)

где р — давление газа.

Уравнение сохранения энергии:

дЕ д ч ч

(1.3)

(1.4)

где у — показатель адиабаты, равный отношению удельных теплоемкостей при постоянном давлении и постоянном объеме. Для двухатомных газов (в том числе воздуха) при не очень высоких температурах показатель адиабаты 7=1,4.

В выражениях (1.1) — (1.4) = 1,2,3, по повторяющимся индексам предполагается суммирование; — символ Кронекера,

О, хф)

Предполагается, что газ является совершенным, так что выполняется уравнение Клапейрона:

Р = РЯТ, (1.5)

где Т — температура, Я — удельная газовая постоянная, равная универсальной газовой постоянной, отнесенной к молярной массе газа.

При записи системы уравнений (1.1) — (1.5) предполагается, что координаты отнесены к некоторому характерному масштабу длины Ь, плотность и скорость — к их значениям в некотором состоянии, например, в набегающем потоке, р«,, время — к величине Ь/и.«,, давление — к р^и2.

В вычислительной аэродинамике система уравнений Эйлера часто записывается в следующем виде:

+ ^ + ^ + (16) <9? дх ду дг '

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

Q = (p,pu,pv,pw,E)т, Г, О, Н — векторы потоков вдоль, соответственно, осейх,у иг:

¥ = (ри,ри2 + р,риу,рию,и(Е + р)) ,

т

(ру,риу,ру2 + р,руц>,у(Е+р)) ,

т

И= (рм?,рим?,рум',р\у2 + р,ы(Е + р)) .

При численном решении уравнения Эйлера обычно используются в интегральной форме. Интегрируя (1.6) по некоторой области (контрольному объему) £2, ограниченному замкнутой поверхностью Г = дО., получаем

где при переходе к поверхностному интегралу использована теорема Гаусса—Остроградского. В (1.7) через пх, пу и пг обозначены составляющие вектора внешней нормали к поверхности Г, ёБ — элемент площади поверхности.

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

1.2. Метод сквозного счёта для решения уравнений газовой динамики

Численное решение системы уравнений Эйлера при высоких скоростях потока сопряжено со значительными трудностями, связанными с наличием в течении сильных ударных волн и других газодинамических разрывов. Причём, разрывные течения могут возникать из гладких начальных данных. Такие особенности течения предъявляют к используемым численным алгоритмам до некоторой степени противоречивые требования. С одной стороны, численный метод должен сохранять монотонность в тех областях, где решения имеют сильные разрывы. С другой стороны, желательно, чтобы метод был возможно более высокого (как минимум второго) порядка точности там, где решение является гладким. Теорема Годунова [34] показывает, что в рамках линейных разностных схем эти два требования невозможно удовлетворить одновременно.

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

йхйу(1г —

(1-7)

г

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

Другим решением этой проблемы является использование схем сквозного счёта, в которых происходит размазывание разрыва на отрезке, величина которого определяется численной диссипацией. Первые схемы сквозного счёта, появившиеся в начале 50-х годах прошлого века [35], основывались на явном введении в правую часть уравнения Эйлера искусственной вязкости. В дальнейшем существенно большее распространение получил иной подход, при котором численная диссипация появляется как результат используемой разностной аппроксимации решаемых уравнений. В таких случаях обычно говорят о «схемной» диссипации.

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