Моделирование трансзвуковых течений является одной из актуальных задач вычислительной аэродинамики. На трансзвуковом режиме проводит большую часть полета гражданская и военная авиация, течения внутри многих реактивных двигателей также являются трансзвуковыми. Вычислительные эксперименты в данной области имеют высокую научную и практическую значимость. Вместе с тем, в связи с рядом физических особенностей трансзвуковые течения крайне сложны для моделирования.
В первую очередь следует отметить, что в указанных задачах в большей части расчетной области течение дозвуковое, а значит возмущения в нём могут распространяться вверх по потоку. Это накладывает требования на удаленность входной границы от тела, что приводит к увеличению размера сетки. Вместе с тем, характерной чертой трансзвуковых течений являются локальные сверхзвуковые области, замыкающиеся скачком уплотнения. Наличие в течении ударных волн требует рассматривать задачу с учетом сжимаемости.
Отдельно необходимо упомянуть отмеченное в экспериментах по трансзвуковой перестройке течения [1−4] явление аэродинамического гистерезиса. Суть данного явления заключается в том, что конкретная конфигурация течения при одном и том же числе Маха может зависеть от того, в каком направлении оно изменяется (переход от дозвукового обтекания к сверхзвуковому или обратно). Моделирование таких эффектов требует нестационарной постановки задачи с зависящими от времени граничными условиями.
В большинстве задач практической важности рассматриваемое течение является отрывным. Отрыв пограничного слоя индуцируется неблагоприятным градиентом давления, который обусловлен либо геометрией обтекаемого тела, либо взаимодействием замыкающего сверхзвуковую область скачка с пограничным слоем. Корректное моделирование указанных явлений требует учета вязкости. Независимо от конкретного механизма его образования, для правильного предсказания отрыва, а также определения положения отрыва и его влияния на исследуемые параметры необходимо определить, какую роль в данном течении играет явление турбулентности.
Согласно определению, данному в 1937 году Тейлором и фон Карманом турбулентность это «нерегулярное движение, которое обычно происходит в жидкостях и газах при обтекании твердой поверхности, или даже при движении слоев одной жидкости друг относительно друга» [5]. Она характеризуется наличием в течении большого диапазона пространственных и временных масштабов. Без преувеличения можно сказать, что все реальные течения, имеющие инженерное приложение, турбулентны. Детальный анализ решений уравнений Навье-Стокса, или их приближения для пограничного слоя, показывает, что турбулентность развивается как неустойчивость ламинарного течения при достижении числом Рейнольдса некоторого критического значения. Для исследования устойчивости ламинарного течения, как правило рассматривают линеаризованные уравнения. И хотя такой подход позволяет получить некоторое представление о начальном этапе развития турбулентности, принципиальная нелинейность уравнений Навье-Стокса не позволяет получить полное аналитическое описание процесса перехода, не говоря уже о самом полностью турбулентном режиме. В реальных (т.е. вязких) течениях неустойчивость возникает из взаимодействия нелинейных инерциальных членов и вязких членов. Сложный характер этого взаимодействия связан с сильной завихренностью течения, его трёхмерностью и нестационарностью. Трёхмерность является неотъемлемой частью любого турбулентного течения, так как без неё не работает механизм растяжения вихрей и каскадной передачи энергии от крупномасштабных энергосодержащих к мелким диссипативным вихрям.
С уверенностью можно сказать, что даже наименьшие из вихрей являются макроскопическими структурами и по своему размеру значительно превосходят любой молекулярный масштаб. Из этого следует, что теоретически, вся физика турбулентных течений должна описываться теми же трёхмерными нестационарными уравнениями Навье-Стокса. Однако, из соотношений, полученных Колмогоровым, следует оценка, согласно которой отношение размера наибольших вихрей к размеру наименьших пропорционально Re3//4. Для оценки снизу размера расчетной сетки можно предположить, что наибольшие вихри занимают всю расчетную область, а линейный размер наименьших соответствует одной ячейке. В таком случае число ячеек трёхмерной сетки составит.
Re¾° = 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. Неоднородность сетки, быстрое растяжение её в каком-либо направлении также негативно сказывается на точности получаемых результатов, а в случае конечно-разностных схем это может привести к потере свойства консервативности.
Перечисленные ограничения приводят к тому, что для расчета турбулентных течений, даже при простой геометрии расчетной области и двумерной постановке должны использоваться сетки содержащие сотни тысяч ячеек. Высокое пространственное разрешение требует соответствующего временного разрешения. Все это делает подобные задачи крайне дорогими с точки зрения потребляемых машинных ресурсов. Если учесть, что для анализа одной задачи, как правило требуется неоднократное проведение расчета (возможно с различными начальными или граничными условиями), становится ясно, что программный код для расчета турбулентных течений должен быть ориентирован на массивно-параллельные вычислительные системы.
Развитие традиционных процессоров в сторону усложнения их внутренней структуры привело к хорошо известному специалистам в области высокопроизводительных вычислений феномену доминирования вспомогательных операций во времени исполнения программ вычислительного характера (см. напр. [16, 17]). Современные процессоры выполняют «полезные» (с точки зрения пользователя-математика) операции с плавающими числами так же быстро, как «вспомогательные, ненужные» операции по вычислению адресов этих чисел, и в десятки, сотни раз быстрее, чем эти числа выбираются из памяти. В итоге «вспомогательные» действия доминируют во времени исполнения программы. Рост скорости выполнения операций над вещественными числами перестает приводить к значительному ускорению счета, поскольку реальные программы тратят львиную долю времени вовсе не на эти операции, а на «вспомогательные» действия.
При наращивании числа традиционных универсальных процессоров в суперкомпьютерных системах разработчик также сталкивается с проблемой масштабируемости приложения. Заключается она в том, что относительная доля операций, необходимых для организации совместной работы большого числа процессоров, быстро растет при разбиении задачи на все меньшие подзадачи. Кроме того, на пути наращивания числа процессоров в системе есть и чисто технические ограничения, связанные с объемом потребляемой энергии количеством выделяемого тепла. Такой компьютер представляет собой сложную и дорогую инженерную конструкцию, включающую в себя систему бесперебойного питания достаточной мощности, систему кондиционирования и систему пожаротушения. Спрос на мощные вычислительные системы заметно опережает предложение и реальных альтернатив нетрадиционным архитектурам нет уже сегодня.
На данный момент, три из пяти самых мощных суперкомпьютеров в рейтинге top500 являются гибридными вычислительными системами, в основе которых лежат графические процессоры [18]. Начинают появляться такие комплексы и в нашей стране. В качестве примера можно привести систему К-100 в ИПМ им. М. В. Келдыша РАН и супперкомпьютер «Ломоносов» в МГУ им. М. В. Ломоносва, занимающий 26 место в рейтинге ТорбОО по сосо-тоянию на ноябрь 2012. Подобные решения существуют не только в виде дорогих промышленных систем — графический процессор с поддержкой CUDA можно найти сегодня практически в любом настольном ПК. В связи с этим встает вопрос о необходимости переноса старых проверенных программ, многие из которых написаны на языке Fortran, на новую архитектуру.
Актуальность работы Данная работа посвящена математическому моделированию и разработке высокоэффективного параллельного комплекса программ для проведения вычислительных экспериментов и прикладных расчетов задач транси сверхзвукового обтекания летательных аппаратов с учетом турбулентности.
Высокая востребованность для инженерных приложений упомянутых расчетов, а также большое количество нерешенных проблем в данной области, обуславливает научную и практическую значимость поставленной задачи.
Таким образом, разработка программного обеспечения для расчета трёхмерных нестационарных сжимаемых вязких течений с учетом турбулентности является актуальной. Не менее важно и эффективное распараллеливание программного обеспечения, позволяющее проводить сложные расчеты на современных высокопроизводительных вычислительных системах, в том числе гибридной архитектуры.
Цель диссертационной работы состоит в создании параллельного программного комплекса для вычислительных систем с гибридной архитектурой, предназначенного для решения инженерных задач. В качестве математической модели выбраны уравнения Рейнольдса, с моделью турбулентности Спаларта-Алмараса. Второй целью является применение комплекса для решения ряда задач, имеющих большое практическое значение.
Для достижения поставленных целей решены следующие задачи:
• Разработка параллельной (MPI) программы на языке Fortran 90 для решения двумерных задач в рамках уравнений Навье-Стокса и Рейнольдса (RANS).
• Разработка вспомогательных модулей и приложений (в том числе эллиптического сеточного генератора, модуля для моделирования условий реального старта и разгона обтекаемого тела, модуля реализующего неявное сглаживание невязки (residual smoothing)).
• Проведение двумерных расчетов обтекания головной части ракеты-носителя. Сравнение результатов расчетов с экспериментом.
• Создание на имеющейся базе программы для расчета трёхмерных течений в рамках уравнений Навье-Стокса, RANS и метода отсоединенных вихрей (DES) [19].
• Реализация трёхмерной версии программы для гибридных вычислительных систем на базе 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 Visualizati and Image Processing» (г. Москва, 2011).
• Международная конференция «16th International Conference on Aerophysics Research Methods» (г. Казань, 2012).
Публикации. Материалы диссертации опубликованы в 2 статьях в рецензируемых журналах [20, 21] из перечня ВАК и 4 статьях в сборниках трудов конференций [22−25] [26].
Структура и объем диссертации
Диссертация состоит из введения, 4 глав, заключения и библиографии. Общий объем диссертации 103 страниц, из них 95 страницы текста, включая 46 рисунков. Библиография включает 68 наименований.
4.6.2. Выводы результатам расчетов трансзвукового обтекания.
1. Выполнены расчет трансзвукового обтекания профиля NACA 0012 M = 0.7, 0.5 < а < 4.8.
2. Для всех указанных углов атаки в режиме RANS получены стационарные решения. Наблюдается вполне удовлетворительное согласие с экспериментальными данными по коэффициентам давления и подъемной силы.
Для M = 0.7 а = 4.8° в режиме DES получено нестационарное решение.
Заключение
.
В процессе работы над диссертацией достигнуты следующие результаты:
1. Разработан и реализован комплекс программ под высокопроизводительные вычислительные системы гибридной архитектуры, позволяющий проводить в том числе и трёхмерные расчеты вязких сжимаемых течений методом DES. Показана эффективность созданных параллельных алгоритмов.
2. В рамках разработанного комплекса реализованы и распараллелены вспомогательные алгоритмы — неявное сглаживание невязки (residual smoothing) и генератор синтетической турбулентности.
3. Реализованный программный комплекс использован при моделировании задач обтекания головной части ракеты-носителя и обтекания крыла с симметричным профилем. Достигнуты следующие результаты:
• для головной части в расчетах получено хорошее совпадение с экспериментом и воспроизведены основные элементы трансзвуковой перестройки течения при увеличении числа Маха набегающего потокаустановлено, что для аккуратного численного моделирования задач трансзвукового обтекания рассматриваемых тел принципиальное значение имеет учет факторов вязкости и турбулентности.
• в двумерных расчетах обтекания крыла удалось удовлетворительно воспроизвести линейный участок экспериментальной кривой зависимости подъемной силы от угла атаки. В трёхмерном расчете, с помощью синтетической турбулентности в набегающем потоке удалось получить совпадение с экспериментом на больших углах атаки.
