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

00500"^

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

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

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

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

программ

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

Р МАР ?0|3

Москва - 2013

005050346

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

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

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

проф.,

Елизарова Татьяна Геннадьевна к. ф.-м. и.,

Георгиевский Павел Юрьевич Ведущая организация: ФГУП ЦНИИмаш

Защита состоится 21 марта 2013г. в _ часов на заседании диссертационного совета Д 002.024.03 при Институте прикладной математики им. М.В. Келдыша РАН, расположенном по адресу: 125047, Москва, Миусская пл., д.4

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

Автореферат разослан «_»_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. Проведено моделирование процесса разгона ракеты из состояния покоя до скорости М=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 Visualiza 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] и нахождение решения линеаризованной задачи методом НЬЬС [11].

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

4

щет = щы + - ДО ¿=1

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

интервале, ограничение на число Куранта практически снимается. В результате мы приходим к системе линейных алгебраических уравнений(СЛАУ):

(1 + пе)Щет RTW = Rf

i

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

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

u = U + и'

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

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

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

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

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

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

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

Для алгоритма сглаживания невязки (residual smoothing) в программе

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

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

Таблица 1. Замеры времени счета

Число ядер Время счета (мин.)

16 189

32 96

64 50

128 27

Ускорение

1 "

>

16 22 48 64 80 Число ядер 96 1 2 128 1

Рис. 1. Зависимость ускорение от числа ядер (из-за размера задачи, ускорение замерялось относительно расчета на 16 ядрах)

Эффективность

IDO 00% ------ —--.___

г:0 00%--------- - --

30 00%--------

70 00%---------

в 60 00%--------

50 00%--------

30 00%--------i

П 10 00%--------

20 00%--------

10 00%--------

0 00% ---------1

0 20 40 60 80 100 120 МО

Число ядер

Рис. 2. Зависимость эффективности от числа ядер

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

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

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

• загрузка данных с жесткого диска

• выгрузка данных на жесткий диск

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

Массивы, участвующие в данных операциях нуждаются в дублировании.

2. Замена основных циклов по сеточным массивам на вызов специальных процедур предназначенных для выполнения на GPU (kernel-процедур).

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

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

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

Таблица 2. Результаты замеров скорости счета.

Вычислители 128xCPU 32xGPU 64xGPU 128xGPU

Скорость (шаг/сек.) 4 13.5 24.5 40

Ускорение (отн. 128х CPU) 1 3.375 6.125 10

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

Рис. 3. Основные элементы поверхности модели.

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

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

Одной важных характеристик является наличие или отсутствие отрыва пограничного слоя по действием замыкающего скачка. Анализ полученных в расчетах данных позволяет утверждать, что при М > 0.95 отрыв не возникает. Отрыв наблюдается при М < 0.90. Область отрыва достигает наибольших размеров при М = 0.85

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

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

Рис. 4. Мгновенное распределение давления на поверхности (а). Зависимость от времени давления в точке и усредненное значение (б).

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

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

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

Рис. 5. Распределение давления на поверхности модели (М=0.95). Решения уравнений Эйлера и Рейнольдса, экспериментальные данные.

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

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

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

В начале главы идет общее описание задачи, а также особенностей и ха-

Рис. 6. Зависимость сопротивления от числа Маха.

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

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

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

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

Рис. 7. Зависимость коэффициента подъемной силы от угла атаки [23].

вом. Усредненный коэффициент Су при этом равен 0.75, что значительно ближе к эксперименту, чем результат первоначального расчета на установление (рис. 7).

Для трёхмерной задачи возмущения генерировались специальном образом, подробно описанном в первой главе. Получившееся в результате решение, как и в двумерном случае является отрывным с наличием периодических структур. Однако, в отличие от двумерного случая оно менее регулярно и после осреднения дает значение Су приблизительно 0.65, что совпадает с экспериментом (рис. 8).

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

Выполнены расчеты обтекания профиля КАСА0012 при числе Маха набегающего потока М = 0.7 и различных углах атаки. В режиме ИА^ для всех рассмотренных углов атаки может быть получено численное решение, близкое к стационарному.

Рис. 8. Слева:Зависимость коэффициента подъемной силы от времени для возмущенного 3D DES решения. Угол атаки а = 14°

Справа.-Обтекание крыла с симметричным профилем, визуализация вихревых структур, с раскраской по модулю скорости (изоповерхность Q = (||П||2 - ||5||2)/2).

Основная особенность обтекания рассматриваемого профиля при M = 0.7 состоит в формировании по мере увеличения угла атаки местной сверхзвуковой области над подветренной стороной. Сверхзвуковая область замыкается внутренней ударной волной.

Расчеты показали, что при достаточно малых углах атаки (а < 1.5°) размеры сверхзвуковой области весьма невелики. При увеличении угла атаки размеры сверхзвуковой области быстро увеличиваются, замыкающий скачок вызывает отрыв пограничного слоя.

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

К настоящему времени в ряде экспериментальных и численных работ (см., напр. (24, 25]) показано, что на трансзвуковых режимах стационарное решение может терять устойчивость, что проявляется как в неединственности, так и нестационарности. Согласно [24] такой эффект может наблюдаться при обтекании профиля NACA 0012 для М = 0.7а = 4.8°.

Рис. 9. Распределение давления в различные моменты времени для нестационарного решения.

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

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

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

I проводить в том числе и трёхмерные расчеты вязких сжимаемых тече-

ний методом DES. Показана эффективность созданных параллельных алгоритмов

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

3. Реализованный программный комплекс использован при моделировании задач обтекания головной части ракеты-носителя и обтекания крыла с симметричным профилем. Достигнуты следующие результаты:

• для головной части в расчетах получено хорошее совпадение с экс-

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

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

Список публикаций

1. Кудряшов И. Ю., Луцкий А. Е. Моделирование турбулентного отрывного трансзвукового обтекания тел вращения // Математическое моделирование. 2011. Т. 23, № 5. С. 71-80.

2. Кудряшов И. Ю., Луцкий А. Е. Адаптация кода для расчета течений вязких жидкостей под гибридные вычислительные системы на базе СиБА-МР1 // Математическое моделирование. 2012. Т. 24, № 7. С. 31-44.

3. Кудряшов И. Ю., Луцкий А. Е. Численное исследование осесимметрично-го трансзвукового обтекания модели с развитыми отрывными зонами // Сборник трудов конференции "Устойчивость и турбулентность течений гомогенных и гетерогенных жидкостей : Доклады Всероссийской молодежной конференции. Вып. XII". 2010. С. 189-192.

4. Кудряшов И. Ю., Луцкий А. Е. Численное моделирование задач трансзвукового отрывного обтекания тел вращения // Струйные, отрывные и нестационарные течения. 2010. С. 283-285.

5. Кудряшов И. Ю., Луцкий А. Е. Vortex shedding from an airfoil visualization // The 8th Pacific Symposium on Flow Visualization and Image Processing. No. 038/1. 2011.

6. Kudryashov I. Y., Lutsky A. E. Numerical simulation of axisymmetrical transonic flow with developed separation zones // XVI International conference on the methods of aerophysical research. Abstracts Pt. II. Kazan: 2012. P. 169-170.

7. Кудряшов И. Ю., Луцкий А. Е., Северин А. В. Численное исследование отрывного трансзвукового обтекания моделей с сужением хвостовой части. Препринт ИПМ №7 2010.

Цитированная литература

8. Волков К. Н., Емельянов В. Н. Моделирование крупных вихрей в расчетах турбулентных течений. ФИЗМАТЛИТ, 2008. ISBN: 978-5-9221-0920-8.

9. Edwards J. R., Chandra S. Comparison of Eddy Viscosity-Transport Turbulence Models for Three-Dimensional, Shock-Separated Flowfields // AIAA Journal. 1996. Vol. 34, no. 4. P. 756-763.

10. Годунов С. К., Забродин А. В. Численное решение многомерных задач газовой динамики. Москва: Издательство «Наука», 1976.

11. Того Е. F., Spruce М., Speares W. Restoration of the contact surface in the HLbRiemann solver // Shock Waves. 1994. Vol. 4. P. 25-34.

12. Jameson A., Mavriplis D. Finite Volume Solution of the Two-Dimensiona I Euler Equations on a Regular Triangular Mesh // AIAA 23rd aerospace sciences meeting. No. AIAA-85-0435. 1985.

13. Jameson A. Transonic Flow Calculations: Tech. Rep. 1651: Princeton University, 2003.

14. Bin Mohamad Badry A. B. Synthetic Turbulence Generation for LES on Unstructured Cartesian Grids: Ph. D. thesis / Cranfield University, School of Mechanical Engineering. 2008.

15. Troldborg N. Actuator LineModelling ofWind TurbineWakes: Ph. D. thesis / Technical University of Denmark. 2008.

16. Gilling L., Sorensen N. N., Rethore P.-E. M. Imposing Resolved Turbulence by an Actuator in a Detached Eddy Simulation of an Airfoil // EWEC 2009 Proceedings online. EWEC, 2009.

17. Кудряшов И. Ю., Максимов Д. Ю. Моделирование задач многофазной многокомпонентной фильтрации на многопроцессорных вычислительных комплексах // Препринт ИПМ. 2009. № 68.

18. Даньков Б. Н. и др. Особенности трансзвукового обтекания конусоци-линдрического тела при большом угле излома образующей на передней угловой кромке // Изв. РАН. МЖГ. 2006. » 2. С. 46-60.

19. Даньков Б. Н. и др. Особенности трансзвукового обтекания конусоцилин-дрического тела при малом угле излома образующей на передней угловой кромке // Изв. РАН. МЖГ. 2006. № 3. С. 140-154.

20. Даньков Б. Н. и др. Волновые возмущения в трансзвуковых отрывных течениях //' Изв. РАН. МЖГ. 2006. № 6. С. 153-165.

21. Даньков Б. Н. и др. Особенности трансзвукового течения за задней угловой кромкой надкалиберного конусоцилиндрического тела // Изв. РАН. МЖГ. 2007. 3.

22. Петров К. П. Аэродинамика тел простейших форм. Москва: Факториал, 1998. С. 432.

23. Sheldahl R. Е., Klimas Р. С. Aerodynamic Characteristics of Seven Airfoil Sections Through 180 Degrees Angle of Attack for Use in Aerodynamic Anal-

ysis of Vertical Axis Wind Turbines: Tech. rep.: Sandia National Laboratories, 1981.

24. Crouch J. D., Garbaruk A., Magidov D. Predicting the onset of flow unsteadiness based on global instability // Journal of Computational Physics. 2007. Vol. 224. P. 924-940.

25. Barakos G., Drikakis D. Numerical simulation of transonic buffet flows using various turbulence closures // International Journal of Heat and Fluid Flow. 2000. Vol. 21. P. 620 - 626.

О ИПМ им .М.В.Келдыша РАН, 2013

в печать 12.02.2013. Формат 60x84/16. Усл. печ. л. 1,4. Тираж 75 экз. Заказ 2. ИПМ им.М.В.Келдыша РАН. 125047, Москва, Миусская пл., 4

Институт прикладной математики им. М.В. Келдыша Российской академии наук

04201354910

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

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

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

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

программ

ДИССЕРТАЦИЯ на соискание ученой степени кандидата физико-математических наук

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

Москва - 2013

Содержание

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

Глава 1. Алгоритмы...........................18

1.1. Основные уравнения и их аппроксимация ............18

1.2. Особенности трёхмерной реализации...............23

1.3. Модель турбулентности Спаларта-Аллмареса..........24

1.4. Интегрирование по времени ....................28

1.5. Сглаживание невязки........................29

1.6. Синтетическая турбулентность...................31

Глава 2. Параллельная реализация .................38

2.1. Эффективность, ускорение, масштабируемость .........38

2.2. Обзор вычислительных систем...................41

2.3. Программная реализация......................42

2.4. Параллельная реализация сглаживания невязки.........44

2.5. Использование графических процессоров и гибридных вычислительных систем..........................48

Глава 3. Численное исследование обтекания головной части ракеты носителя..............................55

3.1. Введение...............................55

3.2. Постановка задачи..........................57

3.3. Постановка задачи при моделировании разгона ракеты.....58

3.4. Результаты расчетов ........................59

3.5. Результаты моделирования процесса разгона ракеты......67

3.6. Выводы................................71

Глава 4. Численное моделирование течения вокруг крыла с симметричным профилем.......................73

4.1. Введение...............................73

4.2. Основные особенности обтекания профилей крыла........74

4.3. Численное моделирование обтекания крыла в рамках подхода RANS SA................................76

4.4. Трёхмерные расчеты методом DES на дозвуковых режимах . . 82

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

4.6. Расчет обтекания профиля NACA 0012 на трансзвуковом режиме 88

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

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

Введение

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

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

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

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

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

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

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

Re3/4° = Re9/4 > Re2. В расчетах для инженерных приложений число Рейнольдса по порядку величины, как правило не меньше 106. Это соответствует расчетной сетке из более чем 1012 ячеек. Если добавить к этому тот факт, что с измельчением сетки мы вынуждены уменьшать и шаг по времени, становится ясно, что при текущем уровне и темпах развития вычислительной техники, полный расчет турбулентного течения для реальной задачи в рамках чистых уравнений Навье-Стокса не представляется возможным.

Возникает ситуация, когда по техническим причинам выполнить расчет невозможно, однако проводить его надо. Для выхода из этого затруднительного положения, начиная с работы Буссинеска 1877 года [6] и на протяжении всего XX века, разрабатывались специальные математические методы, объединяемые под общим названием «моделирование турбулентности». До 1940-ых годов разрабатывались и вводились основные понятия. На данном этапе, помимо упомянутой работы Буссинеска, в которой высказывается гипотеза о линейной связи между тензором Рейнольдсовых напряжений и тензором скоростей деформаций осредненного течения, можно выделить работы Рейнольдса(1895, вводится осреднение по Рейнольдсу) [7], Прандтля(1925) [8] и Кармана(1930) [9] (понятие пути перемешивания). Начиная с 1940 года

началась разработка математической базы и теоретических основ большинства моделей турбулентности. В этот период, существенный вклад внесли такие ученые как Колмогоров(1942, формула Колмогорова, первая модель к —и) [10], Ротта(1951, первая модель Рейнольдсовых напряжений) [11], Ван-Дрист(1956, демпфирующий множитель Ван-Дриста) [12] и многие другие. Начиная с 60-ых годов XX века началось использование моделей турбулентности для замыкания системы уравнений Рейнольдса. Однако, несмотря на большое число работ в данной области, единого метода, подходящего без существенных изменений и настройки к любой задаче найдено не было. Работы в данном направлении по поиску новых подходов к описанию турбулентности продолжаются до сих пор.

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

При классификации подходов к моделирования турбулентности можно выделить две основные черты. Первая — это способ осреднения исходных непрерывных величин, которая определяет физическую интерпретацию получаемых в расчете полей дискретных значений. По данному признаку подходы делятся на DNS (Direct Numerical Simulation), RANS (Reynolds Averaged Navier Stokes), LES(Large Eddy Simulation) и гибридные RANS-LES методы, наиболее распространенным и проработанным из которых является DES (Detached

Eddy Simulation).

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

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

Метод LES был предложен Смагоринским в 1963 [13] году для моделирования атмосферных течений. Суть его заключается в том, что дискретное поле величин рассматривается как сглаженное(осредненное по пространству) исходное поле непрерывных величин. Получающееся решение является нестационарным. При этом крупные вихри разрешаются явно, а для мелкомасштабных вихрей, как и в методе RANS используется та или иная модель. В основе этого лежит предположение о том, что крупные элементы течения напрямую зависят от граничных условий (в том числе формы обтекаемого

тела), а мелкомасштабная турбулентность практически изотропна и её вид слабо зависит от конкретной задачи. Как правило, предполагается, что размер носителя сглаживающего фильтра в каждой точке равен размеру соответствующей ячейки. Уже из этого видно, что вид получающегося решения сильно связан с тем, как построена сетка. На практике шаг сетки выбирается так, чтобы разрешаемые пульсации попадали в так называемый инерциаль-ный интервал турбулентного спектра. Хотя данный метод и не требует такого высокого разрешения как DNS, тем не менее он по прежнему остается вычислительно сложным. Это вызвано высокими требованиями к разрешению сетки в пограничном слое. Как компромиссный вариант между RANS и LES был разработан метод DES [14]. Суть его заключается в том, что в области пограничного слоя строится стационарное RANS решение, а вне его, на удалении от тела работает метод LES.

По мнению специалистов в данной области [15], метод RANS будет оставаться самым восстребованым в инженерных расчетах еще на протяжении многих лет. Однако для некоторых задач, в том числе рассмотренных в данной диссертации, предпочтительным выбором является LES/DES. Причина этого в том, что подход RANS не обеспечивает достаточной точности при описании принципиально нестационарных течений или наличии интенсивных отрывов. Модели подсеточных напряжений, удовлетворительно описыващие влияние мелкомасштабной изотропной турбулентности, плохо подходят для крупных вихревых структур, возникающих при отрыве, в виду их сильной зависимости от граничных условий. Метод LES менее требователен к подсеточ-ным моделям и опирается на них только в той области, где они действительно могут дать правдоподобный результат.

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

зданных моделей выходит за рамки данной работы. Достаточно сказать, что разработано их большое количество. По этому признаку модели делятся на линейные, использующие гипотезу Буссинеска (они в свою очередь делятся на алгебраические модели, модели с одним уравнением, с двумя уравнениями и т.д.) и на модели Рейнольдсовых напряжений (нелинейные, делятся на дифференциальные, алгебраические, явные алгебраические). На данный момент одними из наиболее распространенных являются две модели: модель Спаларта-Алмараса (линейная, с одним дополнительным уравнением) и модель SST к —и (Menter) (линейная модель с двумя уравнениями). Каждая из них имеет множество различных модификаций, отличающихся сложностью и областью применения.

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

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

Уровень развития вычислительной техники сегодня во многом определяет темпы технического прогресса и успехи в решении фундаментальных научных задач, среди которых общепризнанным является класс проблем Grand challenges: это фундаментальные научные или инженерные задачи с широкой областью применения, эффективное решение которых возможно только с использованием мощных (суперкомпьютерных) вычислительных ресурсов,

с производительностью сотен Gflops 1012 операций в секунду) и выше. Одной из таких задач является расчет нестационарных турбулентных течений методами DES, LES, DNS.

Хорошо известно, что в пристеночной области, составляющей 20% от толщины пограничного слоя генерируется до 80% турбулентной энергии. Поэтому, даже при проведении RANS расчетов стационарных течений, для корректного учета турбулентных эффектов сохраняется необходимость хорошо разрешать пограничный слой. Это приводит к известному ограничению на толщину ячеек у стенки в нормальном к ней направлении — безразмерная величина у+ должна составлять порядка 1, что равносильно тому, что на область пограничного слоя должно приходиться не менее 10 ячеек.

Для корректной работы численных схем и уверенности в адекватности получаемых результатов расчетная сетка должна удовлетворять определенным требованиям. Так, вне области пограниченого слоя, где нет явно выделенных направлений, соотношение производных по различным координатным направлениям определяет допустимое соотношение сторон ячеек сетки (aspect ratio). В области LES ячейки должны быть по возможности более кубическими и отношение сторон не должно слишком сильно отличаться от 1. Неоднородность сетки, быстрое растяжение её в каком-либо направлении также негативно сказывается на точности получаемых результатов, а в случае конечно-разностных схем это может привести к потере свойства консервативности.