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

Кудряшов И.Ю.

Численное моделирование трансзвуковых

отрывных течений

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

программ

АВТОРЕФЕРАТ

диссертации на соискание ученой степени

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

Москва – 2013

Работа выполнена в ИПМ им. М.В. Келдыша РАН.

Научный руководитель: д. ф.-м. н., Луцкий Александр Евгеньевич

Официальные оппоненты: д. ф.-м. н., проф., Елизарова Татьяна Геннадьевна, ИПМ им. М.В. Келдыша РАН, главный научный сотрудник к. ф.-м. н., Георгиевский Павел Юрьевич НИИ Механики МГУ им. М.В. Ломоно­ сова, ведущий научный сотрудник

Ведущая организация: ФГУП ЦНИИмаш

Защита состоится 21 марта 2013г. в часов на заседании диссертаци­ онного совета Д 002.024.03 при Институте прикладной математики им.

М.В. Келдыша РАН, расположенном по адресу: 125047, Москва, Миусская пл., д.

С диссертацией можно ознакомиться в библиотеке ИПМ им. М.В. Келдыша РАН.

Автореферат разослан « » 2013 г.

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

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

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

Актуальность работы Данная работа посвящена математическому мо­ делированию и разработке высокоэффективного параллельного комплекса программ для проведения вычислительных экспериментов и прикладных рас­ четов задач транс- и сверхзвукового обтекания летательных аппаратов с уче­ том турбулентности.

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

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

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

Для достижения поставленных целей решены следующие задачи:

Разработка параллельной(MPI) программы на языке Fortran 90 для ре­ шения двумерных задач в рамках уравнений Навье-Стокса и Рейнольд­ са (RANS).

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

Проведение двумерных расчетов обтекания головной части ракеты-носи­ теля. Сравнение результатов расчетов с экспериментом.

Создание на имеющейся базе программы для расчета трёхмерных тече­ ний в рамках уравнений Навье-Стокса, RANS и метода отсоединенных вихрей (DES) [8].

Реализация трёхмерной версии программы для гибридных вычисли­ тельных систем на базе CUDA-MPI на языке CUDA Fortran.

Проведение двумерных расчетов методом RANS и трёхмерных расчетов методом DES обтекания крыла. Сравнение с данными эксперимента.

Разработка генератора синтетической турбулентности для задания па­ раметров турбулентности набегающего потока при расчетах методом DES.

Проведение трёхмерных расчетов обтекания крыла с генератором син­ тетической турбулентности. Сравнение с экспериментом.

Исследование ускорения и эффективности параллельных алгоритмов в режимах MPI и CUDA-MPI.

Научная новизна диссертации отражена следующими элементами:

1. Разработан и реализован параллельный комплекс программ для рас­ чета вязких сжимаемых течений в рамках уравнений Навье-Стокса, RANS и DES на гибридных вычислительных системах с графически­ ми процессорами NVIDIA. В процессе расчетов подтверждена высокая работоспособность и эффективность комплекса.

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

3. Проведено моделирование процесса разгона ракеты из состояния покоя до скорости M=1.3 с использованием данных телеметрии по ускорению.

4. Показано, что наличие разрешаемых турбулентных пульсаций в набега­ ющем на крыло потоке может значительно улучшить качество 3D DES расчета на больших углах атаки.

Практическая значимость. Результаты, изложенные в диссертации, использованы для проведения численных экспериментов и решения инже­ нерных задач. В ходе работы над диссертацией получен опыт по переносу Fortan приложений на гибридную архитектуру и освоен новый компилятор PGI CUDA Fortran; выявлены его возможности, недостатки и ограничения.

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

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

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

2. В рамках разработанного комплекса реализованы и распараллелены вспомогательные алгоритмы — неявное сглаживание невязки (residual smoothing) и генератор синтетической турбулентности.

3. При помощи созданного программного комплекса:

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

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

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

Молодежная конференция Устойчивость и турбулентность течений го­ могенных и гетерогенных жидкостей (г. Новосибирск, 2010).

XXII Юбилейный семинар с международным участием Струйные, от­ рывные и нестационарные течения (г. Санкт-Петербург, 2010).

Международная конференция The 8th Pacific Symposium on Flow Visualization and Image Processing (г. Москва, 2011).

Международная конференция 16th International Conference on Aerophysics Research Methods (г. Казань, 2012).

Публикации. Материалы диссертации опубликованы в 2 статьях в ре­ цензируемых журналах [1, 2] из перечня ВАК и 4 статьях в сборниках трудов конференций [3–6] [7].

Структура и объем диссертации Диссертация состоит из введения, 4 глав, заключения и библиографии. Общий объем диссертации 103 страниц, из них 95 страницы текста, включая 46 рисунков. Библиография включает 68 наименований.

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

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

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

Интегрирование по времени, в зависимости от настроек расчета осуществ­ ляется либо при помощи явной схемы Эйлера, либо пятистадийным методом Рунге-Кутта.

Вычисление потоков на границах ячеек осуществляется с помощью ре­ шения задачи о распаде разрыва. Для этого используется один из двух ре­ ализованных алгоритмов: нахождение точного решения нелинейной задачи, методом Ньютона [10] и нахождение решения линеаризованной задачи мето­ дом HLLC [11].

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

Однако, явная процедура осреднения не позволяет заметно увеличить максимально допустимый шаг по времени [12]. В работе [13] показано, что при неявном осреднении со значением параметра, лежащем в определенном интервале, ограничение на число Куранта практически снимается. В резуль­ тате мы приходим к системе линейных алгебраических уравнений(СЛАУ):

Большое количество научных и инженерных задач требует в процессе численного моделирования создания искусственной турбулентности. В число указанных задач входит моделирование взаимодействия летательных аппа­ ратов со свободной турбулентностью, анализ влияния атмосферной турбу­ лентности на распространение электромагнитных волн радио- и оптического диапазона и многие другие. Широкое применение получил также данный под­ ход при моделировании течений вихреразрешающими методами DNS и LES, которые для корректной работы требуют задания на входной границе не толь­ ко параметров осредненного течения, но и параметров турбулентности(а они, как правило, априори не известны).

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

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

Суть метода, применённого нами, описана в работе [14]. Предполагается, что корреляция компонент скорости не зависит от координаты и времени, а зависит только от дистанции, и единственного линейного масштаба.

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

Для наложения пульсаций на основное течение нами применен подход актуатора, предложенный в статьях [15, 16], однако существенным отли­ чием является то, что в указанных статьях рассматривался несжимаемый случай, а численное решение осуществлялось методом конечных разностей с помощью алгоритмов PISO или SIMPLE. Для реализации актуатора при этом использовались добавки к давлению, величина которых вычислялась из уравнения Бернулли в приближении стационарного одномерного течения.

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

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

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

Разбиение области на подобласти осуществляется автоматически. Алго­ ритм разбиения, разработанный при участии автора, изначально был пред­ назначен для симулятора МКТ, моделирующего процессы фильтрации в пла­ сте [17]. Он многократно опробован и подтвердил свою эффективность в ходе промышленных расчетов.

Для алгоритма сглаживания невязки (residual smoothing) в программе реализована параллельная версия метода Якоби. Расчеты показывают, что на реальных задачах на операцию сглаживания тратиться в 5-6 раз меньше времени, чем на вычисление потоков или вычисление производных на гранях.

Ниже представлены результаты замеров ускорения полученные при рас­ чете задачи двойного маховского отражения (double mach reflection) в трёх­ мерной постановке. Расчетная область состояла из 15827000 ячеек. Расчеты проводились на кластере МВС-100к МСЦ. Из представленных на рис. 1 и графиков видно, что разработанный алгоритм обладает хорошими показате­ лями эффективности и масштабируемости.

Рис. 1. Зависимость ускорение от числа ядер (из-за размера задачи, ускорение замерялось относительно расчета на 16 ядрах) Рис. 2. Зависимость эффективности от числа ядер В ИПМ им. М.В. Келдыша действует вычислительная система гибрид­ ной архитектуры К-100. В числе прочего математического обеспечения на ней установлены компиляторы Portland Group, что дало возможность дора­ ботать программу под использование технологии CUDA без переноса её на C/C++ и определить характеристики полученного параллельного кода.

Условно процесс переработки программы можно разделить на три основ­ ных этапа:

1. Создание необходимых дублирующих структур данных на видеокарте (массивы основных переменных, данные о сетке, константы и т.д.). Все расчеты будут проводиться в памяти граф. плат, однако есть три про­ цесса, которые принципиально осуществляются только через обычную загрузка данных с жесткого диска выгрузка данных на жесткий диск обмен данными по MPI на границах областей различных процессов Массивы, участвующие в данных операциях нуждаются в дублирова­ 2. Замена основных циклов по сеточным массивам на вызов специальных процедур предназначенных для выполнения на GPU (kernel-процедур).

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

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

Ускорение работы, по сравнению со счетом на CPU, составило для фи­ нальной версии от 10 до 16 раз, в зависимости от конкретной задачи и числа используемых устройств. Так, в задаче по обтеканию крыла, на сетке из млн. ячеек при использовании 128 CPU скорость счета составляла 4 шага в секунду. При счете на 32 GPU скорость счета составила 13.5 шагов в секунду, на 64 — 24.5, а на 128 — 40 шагов в секунду.

В третьей главе описано применение разработанного программного обеспечения для решения задачи об обтекании головной части ракеты носите­ ля. Представлено сравнение численных результатов, полученных в рамках мо­ делей уравнений Эйлера, Навье-Стокса и осредненных уравнений Рейнольдса Рис. 3. Основные элементы поверхности модели.

с моделью турбулентности Спаларта-Аллмараса с экспериментальными дан­ ными [18–21]. Также приведены предварительные результаты моделирования разгона ракеты и сравнение их с результатами расчетов на установление.

Расчеты обтекания модели (рис. 3) выполнены для чисел Маха набегаю­ щего потока M = 0.82, 0.90, 0.95, 1.1, 1.3 при числах Рейнольдса Re = 3 * 106 (D — диаметр модели).

Одной важных характеристик является наличие или отсутствие отрыва пограничного слоя по действием замыкающего скачка. Анализ полученных в расчетах данных позволяет утверждать, что при M 0.95 отрыв не возника­ ет. Отрыв наблюдается при M < 0.90. Область отрыва достигает наибольших размеров при M = 0. Решение уравнений Навье-Стокса имеет нестационарный характер. Над цилиндрической поверхностью образуются вихри (области отрыва), которые сносятся вниз по потоку рис. 4а. На рис. 4б представлена зависимость от вре­ мени давления в точке = 1.78. Осредненная по всему промежутку времени величина существенно отличается, как от экспериментального значения, так и от решения уравнений Рейнольдса.

Течения около обратных уступов также является классической областью исследований (см. напр. [22] и приведенный там список литературы). Извест­ но, что положительный градиент давления, как правило, приводит к отрыву и образованию обширной области обратного течения. Модель уравнений Эй­ Рис. 4. Мгновенное распределение давления на поверхности (a). Зависимость от времени давления в точке и усредненное значение (б).

лера, естественно, не описывает такие явления. Сопоставление численных и экспериментальных данных по распределению давления на поверхности моде­ ли (рис. 5) свидетельствует о неадекватности модели Эйлера в этой области.

Указанный факт существенно влияет на получаемую в расчете величину со­ противления. Наличие для решений уравнений Эйлера области низкого дав­ ления на обратном конусе приводит к росту суммарного сопротивления на 30% по сравнению с решением уравнений Рейнольдса.

Приведенные данные свидетельствуют, в первую очередь, о необходимо­ сти при численном моделировании рассматриваемых течений учитывать вли­ яние вязкости и турбулентности. Турбулентный характер течения в рассмат­ риваемом случае удовлетворительно описывается моделью Спаларта-Аллма­ раса.

Во второй части третьей главы приведены результаты расчетов, модели­ рующих процесс разгона ракеты. Получена серия мгновенных картин течения для различных чисел Маха в интересующем диапазоне. Поскольку нами так­ же проводились расчеты на установление по данной задаче, было выполнено Рис. 5. Распределение давления на поверхности модели (M=0.95). Решения уравнений Эйлера и Рейнольдса, экспериментальные данные.

сравнение результатов обоих методов. Установлено, что мгновенная картина течения совпадает с полученной в расчете на установление для соответству­ ющего числа Маха.

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

Полученные в одном расчете данные можно использовать для построения зависимостей интегральных характеристик от числа Маха (рис. 6). Кроме то­ го, данный подход позволяет наблюдать за процессами происходящими при трансзвуковой перестройке течения, которые не воспроизводятся в расчетах стационарных режимов.

Четвёртая глава посвящена численному моделированию течения во­ круг крыла с симметричным профилем с использованием подхода RANS до­ полненного моделью турбулентности Спаларта-Алмараса, приводится сравне­ ние с экспериментальными данными, устанавливается область применимости используемых методик.

В начале главы идет общее описание задачи, а также особенностей и ха­ Рис. 6. Зависимость сопротивления от числа Маха.

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

Формулируется постановка задачи об отекании профиля NACA0012, опи­ сываются начальные и граничные условия и сетка. Приводятся результаты расчетов в двумерной постановке по методу RANS(рис. 7). В целом получен­ ные результаты вполне адекватны и соотетствуют тому, что может дать метод RANS. Для коэффициента подъемной силы расчеты дают хорошее совпаде­ ние с экспериментом на линейном участке. При этом момент начала свалива­ ния (резкого падения подъемной силы) оказывается затянут. Получающиеся течения на всех режимах является стационарным.

Третья часть главы посвящена трёхмерным расчетам обтекания крыла с профилем NACA0012 методом DES на дозвуковых режимах. Рассказыва­ ется об особенностях построения сетки и постановке задачи. Обосновывается необходимость применения генератора синтетической турбулентности.

На первом этапе исследований возмущения вносились в установившееся двумерное RANS решение для угла атаки = 14. Из эксперимента извест­ но [23], что при данном угле атаки коэффициент примерно равен 0.65 и крыло находится в режиме сваливания. В результате однократно внесенных возмущений образовалось новое устойчивое периодическое решение с отры­ Рис. 7. Зависимость коэффициента подъемной силы от угла атаки [23].