Численное моделирование трансзвуковых отрывных течений тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат физико-математических наук Кудряшов, И. Ю.

• Кудряшов, И. Ю.
• кандидат физико-математических наук
• 2013

Оглавление диссертации кандидат физико-математических наук Кудряшов, И. Ю.

Содержание

Введение

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1. Введение

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

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

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

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

3.6. Выводы

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

4.1. Введение

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

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

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

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

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

Заключение

Литература

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

Введение

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

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

Отдельно необходимо упомянуть отмеченное в экспериментах по трансзвуковой перестройке течения [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. Неоднородность сетки, быстрое растяжение её в каком-либо направлении также негативно сказывается на точности получаемых результатов, а в случае конечно-разностных схем это может привести к потере свойства консервативности.

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

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

Развитие традиционных процессоров в сторону усложнения их внутренней структуры привело к хорошо известному специалистам в области высокопроизводительных вычислений феномену доминирования вспомогательных операций во времени исполнения программ вычислительного характера (см. напр. [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 наименований.

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

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

Уравнения Навье-Стокеа и Рейнольдеа с 2 пространственными переменными (плоские или осесимметричные течения) при использовании гипотезы Буссинеска могут быть записаны в единообразной форме.

ди дР дв тт

т

и = (р, ри, ру, е)

У

= (ри, ри2 + р, риу, (е + р) и)Т ,

Р — (О, Тхх, ТХу, ^Тхх иТху Цх) ;

X1

Сг = [ру, риу, рь + р,(е + р)у) ,

= (0) ^жг/5 Туу, иТХу ИТуу С[у) ,

Нг = (—ру, —риу, —ру2, — (е + р) у)Т ,

гр

Н — (0, 7"®?/) Туу ^(/3; Т^Тху ~Ь иТуу Цу)

а; = ( 0 — в плоском случае, 1 — в осесимметричном) Компоненты тензора вязких напряжений задаются соотношениями

2 2

Тжа; = - (^ + 1н) (2г4ж - , Туу = - (р + р^ (2уу - их) ,

у 1

тжу = = (/2 + рь) (иу + г^) , г^ = -2(^ + + ~(Иу{У}).

У 3

Рис. 1.1. Схема расчета ячейки в двумерном алгоритме.

Турбулентная вязкость необходима для замыкания осредненных по Рейнольдсу уравнений Навье-Стокса, с моделью турбулентности, использующей гипотезу Бусинеска. Молекулярная вязкость вычисляется по степенной формуле ц = цо(Т/Т0)0-76

Для вычисления давления данная система должна быть дополнена уравнением состояния. В нашем случае в качестве уравнения состояния выбрано уравнение состояния идеального газа: р = рє{7 — 1)

Уравнение для турбулентной вязкости, как правило, может быть записано в форме, аналогичной уравнениям системы (1.1), поэтому расчет производится единообразно для всех уравнений. Численный алгоритм строится методом конечных объемов.

7771+1 _ ттп 4

П^+№+1/2 '^+1/2,т/2 + + сщ)г = ш {12)

г=\

Здесь — площадь ячейки, (Л^, Иу) — внешняя нормаль. Рассмотрим аппроксимацию потоков і*1, Є на примере ребра с номером к+1,1+1/2 — см. рис. 1.1 . Пусть С4+і/2,г+і/2 ~~ величины, отнесенные к центрам ячеек,

ик,1 = (£4+1/2,2+1/2 + £4+1/2,2+1/2 + £4+1/2,2+1/2 + £4+1/2,2+1/2)/4- величины в узлах сетки.

Производные (11х, иу)к+1: ¿+1/2, входящие в выражения потоков вы-

числяются через разности (£4+3/2,2+1/2 ~ £4+1/2,2+1/2) и (£4+1,2+1 - £4+1,1)

Значения функций £4+1,2+1/2 на ребре определяются из решения задачи Римана с начальными данными

£4+1 = £4+1/2,/+1/2 + 1/2Аж£/а.)/к+1/21г+1/2 + 1/2Ду£/2/;/с+1/2,г+1/2 £4+1 = £4+3/2,2+1/2 - 1/2Да;^4,/с+з/2,г+1/2 - 1/2А|/£/г/^+з/2,г+1/2

Для обеспечения монотонности разностной схемы производные в ячейках определяются в соответствии с принципом минимума модуля производных на противоположных ребрах

£4,&+1/2,¿+1/2 = тттос1(£4,/с,г+1/2) £4^+1,2+1/2) иу>к+1/2,1+1/2 = тттос!(£/у5/С!г+1/2, £4^+1,2+1/2) Таким образом, аппроксимация конвективных членом аналогична разностной схеме В.П.Колгана[27].

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

выглядит следующим образом:

ЯГ Г ЯТ? ЯП -ц- иь ил иу '

и = (р,ри,ру,е)т,

Р = рг + ^ С = Сг +

Н = ~(Нг + Нь) + #а, Я° = (0, ра, 0, рпа)т.

'г;

<у

(1.3)

В дополнительной правой части На, а обозначает величину ускорения. В работе [28] показано, что модель турбулентности Спаларта-Алмараса в данной системе отсчета остается неизменной.

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

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

В инерциальной системе отсчета явная схема будет устойчивой если возмущение. распространяющееся со скоростью (и + с) за время шага пройдет расстояние меньшее либо равное шагу сетки (рассматриваем только самую быструю волну). Т.е. должно быть выполнено условие (u+c)At < К. В неинерциальной системе отсчета это же возмущение за время А1 пройдет расстояние (■и + с)Д£ + аД£2/2. В результате мы приходим к видоизмененному условию

КФЛ:

(и + с)At + a(A¿)2/2 ^ 1

G P L = -:- < 1, или же

h

-0 + c) + + с)2 - 2ah

At < —---—---

a

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

В нашей программе реализовано два алгоритма решения задачи Римана. Первый и основной метод — это нахождение точного решения нелинейной задачи, методом Ньютона [29]. Второй метод — Harten-Lax-van Leer (Contact) или HLLC использует линеаризацию задачи Римана, является вспомогательным и применятся на первых стадиях метода Рунге-Кутта для ускорения работы программы.

Метод HLLC является развитием двухволнового метода HLL[30] и впервые был предложен в работе [31] . Отличием от оригинального метода является учет в решении контактного разрыва.

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

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

Рис. 1.2. Иллюстрация к методу НЬЬ.

и вихрей. Метод НЬЬС, вводя в рассмотрение контактные разрывы, позволяет избавится от этих недостатков и дает решение по качеству сравнимое с точным.

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

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

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

Для вычисления градиента в ячейке используется теорема Грина:

дгас1 (/) ¿V = О п до.

или переходя к дискретной форме записи:

г

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

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

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

* ~ ( X3 Р"

\н = иР", и = , , гз , Х = С1-4)

Величина фигурирующая в определения турбулентной вязкости (1.4), находится из решения дифференциального уравнения

дрй дрйщ , ,

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

1

х 2 У \

(1.6)

Gu — производство турбулентности

Су = Сы&й, § =

5= УЩЩ (1.7)

с _ 1 (диг . ди3 \ 2 (дик \ °г3 2 удх3 "т" дхг) 3 \дхк) '

— диссипацию турбулентности

1

Ь = 9 Я-г + С^-т),

Величина ¿и, в вышеприведенных формулах есть расстояние от твердой стенки, что в данной модели характеризует линейный масштаб турбулентности. Величины Сы, Сь2,Си\, С№х, С.^, С^з ~~ константы модели ЗА.

0,1 0V Сы С}>2 Cw 1 Cw 2 Cw з К

7.1 2 3 0.1335 0.622 г* _ Сы 1 1+Сь2 - «2 + ffi/ 0.3 2.0 0.41

Здесь представлена одна из модификаций модели Спаларта-Аллмара-са — S A-Edwards [32]. В качестве граничных условий используется 0 на стенке и значение молекулярной вязкости на бесконечности. В начальный момент времени во всей области значение модельной переменной равно значению молекулярной вязкости на бесконечности.

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

Трудности, возникающие при моделировании пристеночных течений с помощью LES можно преодолеть, если рассматривать гибридные модели, такие как метод отсоединенных вихрей (DES) [33]. Анализ возможностей и ограничений двух традиционных подходов к моделированию турбулентности, таких как RANS и LES, показал, что проблема расчета турбулентных отрывных течений может быть успешно решена путем создания «гибридного» подхода. Этот подход сочетает в себе сильные стороны этих двух методов — высокой точности существующих полуэмпирических RANS моделей в областях присоединенного пограничного слоя и универсальности и приемлемых вычислительных затрат LES в отрывных областях потока.

Метод моделирования отсоединенных вихрей представляет собой гибридную модель RANS/LES, что является разумным компромиссом RANS и LES методиками. Моделирование крупных вихрей (LES) включается только в областях течения, где размер (шаг) сетки Д достаточен для разрешения турбулентности с масштабом LtUrb (А LtUrb)- RANS работает в остальной области течения (Д > Lturb). Например, для течений с обширной отрывной зоной это означает следующее:

1. область присоединенного пограничного слоя считается по RANS подходу, который обеспечивает хорошую точность в подобной ситуации;

2. отрывная область, с относительно крупными вихрями, рассчитывается в LES моде.

Основная идея DES состоит в использовании одной и тоже модели турбулентности в RANS и LES модах. LES мода представляет собой подсеточную версию модели Спаларта-Аллмараса и получается заменой линейного масштаба турбулентности dw из RANS на масштаб турбулентности С desд из модели LES в (1.8). В результате получаем подсеточную модель из одного уравнения — версия модели SA с моделированием вихрей подсеточного масштаба

(SGS — SubGrid Scale). В «равновесии», когда генерация равна диссипации, эта A4 одел ь эквивалентна модели Смагоринского — ¡it = р(С&А) у/2 <SV,SV,, где Cs = С des ■ Сшивка RANS и LES мод состоит во введение масштаба

lDES = min (dw, CDESA), CDes = O(l)1, (1.9)

где А — максимальный диаметр расчетной пристеночной ячейки. Величина Ides заменит масштаб турбулентности dw в уравнениях (1.8), определяющих величину вихревой (турбулентной) вязкости.

1.3.1. Аппроксимация уравнения SA

Для решения дифференциального уравнения модели Спаларта-Аллмара-са используется метод конечных объемов. Роль переносимой величины играет

у/рй.

Легко заметить, что конвективные потоки для уравнений" неразрывности и SA имеют схожий вид: (рип)гз и (рипь>) соответственно (индексы ij означают, что поток вычисляется на грани между ячейками inj). Поэтому, если известен поток плотности между двумя ячейками, конвективный поток величины V может быть вычислен простым домножением на Uk, где индекс к определяется по принципу сноса по потоку:

(pUnXj еСЛИ (Рип)г3 >

KPWhj = S

[ (PUn)ij Vji если {pun)t3 < 0;

Для вычисления величины Cb2 {^d^j входящей в состав диффузионного члена необходимо аппроксимировать градиент модельной переменной в центре ячейки. Для этого используется теорема Грина:

1 В большинстве случаев величина берется равной Cd es = 0 65

дгас1 (у/ру) (IV = о {у/рй) ^с1Б п дп

Переходя к значениям в ячейках:

Г Л( г~\ 1 1 Ур»]%+Ур*\3. \grad {у/ри)]г = ^г -2-Пу

1 з

где суммирование проводится по всем граням ячейки 1, 5г обозначает площадь 1 (или объем в трёхмерном случае), а п^ обозначает внешнюю по отношению к ячейке 1 нормаль, к их общей грани с ячейкой ^

Для аппроксимации члена ^ {(м + у/р») воспользуемся форму-

лой (Ну {к ■ дга(1{/)} ¿V = к ■ . То есть на каждой грани нам

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

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

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

+ - от] = £ X ЯНЗ^),

где <5г обозначает площадь ячейки, Р{у/г) представляет собой сумму конвективных потоков, 0(У/г) — сумму вязких диффузионных потоков, а ЛНЗ(\Уг) — работу источников. В более компактном виде эта система может быть записана как

=

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

ууп+1 = ууп +

Второй — это пятистадийный метод Рунге-Кутты. Он может быть записан следующим образом:

ж(0) = жп

ж(1) =

Ж; + а! —

цгп+1 =

113 1

ос\ = а2 = = -, а4 = а5 = 1

При использовании метода Рунге-Кутты на первых четырёх стадиях невязкие потоки вычисляются из решения линеаризованной задачи Римана методом НЬЬС. На заключительной пятой стадии потоки вычисляются по точном решению.

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

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

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

4

Я?™ = Я?* + е^Я, - Яг)

7=1

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

(1 + ПЕ)Я™Ш ПТШ = ПОгМ

3

В наших расчетах применяется значение параметра е = 0.5. Параметр п обозначает количество соседних ячеек и равен 4 в двумерном случае и 6 в трёхмерном. Для ячеек лежащих вне расчётной области значение Я^еги принимается равным 0.

Несмотря на то, что размерность такой системы может быть достаточно велика ^ х N. где N — это число ячеек расчетной области), а решать её необходимо на каждом шаге по времени (а в случае Рунге-Кутта и на каждой стадии), это не слишком замедляет расчет, так как получившаяся матрица не зависит от времени, симметрична и обладает свойством сильного диагонального преобладания. СЛАУ может быть приближенно решена достаточно дешевыми методами, а информация о факторизации матрицы, получившаяся при первом её решении может быть использована в дальнейшем, т.к. меняется только правая часть системы.

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

Применение сглаживания невязки позволяет заметно сократить скорость сходимости к стационарному решению. На тестовых расчетах задачи обтекания профиля в трёхмерной постановке использование метода Рунге-Кутта совместно со сглаживанием позволило во-первых за одно и тоже процессорное время сосчитать в 1.5 раза больший интервал псевдо-времени, и во-вторых для одних и тех же величин псевдо-времени получать решение более близкое к стационарному. На рис. 1.3 приведены поля плотности во всей расчетной области, полученные методом Эйлера (слева) и Рунге-Кутта со сглаживанием (справа) на один и тот же момент псевдо-времени. Из рисунка видно, что при расчете вторым методом возмущение от крыла распространилось дальше. В совокупности, экономия процессорного времени составляет в среднем около 2 раз.

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

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

%

%

-40 -20

О 20

X

0 20 X

Рис. 1.3. Поле плотности: слева метод Эйлера, справа Рунге-Кутта со сглаживанием (моменты псевдо-времени совпадают).

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

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

u = U + и'

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

числе и крупномасштабных энергосодержащих, что приводит к быстрому его распаду из-за процессов диссипации [37].

Для того, чтобы избавиться от упомянутых недостатков на генерируемые случайные сигналы накладывают дополнительные ограничения, определяемые известными статистическими параметрами турбулентного течения. Этими параметрами могут быть пространственная и временная корреляции, характерные линейный и временной масштабы, энергетический спектр пульсаций, тензор напряжений Рейнольдса. Последний, в частности, является по сути тензором ковариации его компонент (и'г,и'3) и необходим для описания анизотропии турбулентности (например при моделировании турбулентного погранслоя). Для его учета принят подход впервые описанный в статье [38]. Для того, чтобы компоненты поля пульсаций обладали необходимой корреляцией используется разложение тензора Г1 по Холецкому К = АТА, где

А = КЛ =

л/7й~і О О

— у/#22 - О

о11 о

#31 ІХ32 — «21 «31

V#33 - «31 - «32

(1.10)

«11 «22

После этого, компоненты искомого и' (г, ¿) можно представить как и'г = аг]ь3, где компоненты векторного поля у(г, £) являются независимыми случайными сигналами с нулевым математическим ожиданием и единичной дисперсией, т.е. (ьг) = 0 и (уг,у3) = 6гз. В результате задача свелась к определению трёх случайных сигналов, обладающих указанными свойствами. Кроме того, в поле у(г, ¿) могут быть учтены другие статистические свойства турбулентности.

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

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

Другим популярным классом методов является фильтрация сгенерированного сигнала у(г, £). Наиболее популярным является Гауссовский фильтр. К данному классу относится и метод, примененный нами. Суть данного метода описана в работе [39]. Приведем здесь краткое его описание. Предполагается, что корреляция компонент скорости не зависит от координаты и времени, а зависит только от дистанции, и единственного линейного масштаба Ь. На основании сделанных предположений строится модельная функция корреляции:

Д««(г) = ехр (-^Г) (1-11)

Или, переходя к дискретным величинам

= (1.12) (ит,ит) 4(пД)2у

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

n

ит = ^ ^ ЬпГт+п, (1-13)

п

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

(ит, ит+к) _ Ь3Ь3-к ^

(ит,ит) Т^=-мЬз2

В результате мы приходим к уравнению для коэффициентов фильтра:

Егі-м+кЬзЬз-к ( ттк2

-ехр{-^) (1Л5)

Приближенным его решением является

Ьк=ехр (-Й?) (1Л6)

Процесс генерации турбулентности на входной границе х = 0 с размерам по у иг Му и Мг соответственно можно разделить на следующие шаги:

1. Из априорных соображений выбирается линейный масштаб Ь.

2. В прямоугольной области размерами [—Ых : —А^+1 : Му+Му, —

1 : Мг + генерируются три случайных сигнала га, а = 1,2,3 удовлетворяющие условиям (га) = 0, (гг,г3) = 510. А^, Ыу и А"2] выбираются так, чтобы размеры области по каждому из направлений были не меньше 2 Ь.

3. С помощью (1.16) вычисляются коэффициенты фильтра.

4. Производится фильтрация сгенерированного поля га

Щ N.

Уа

(¿*0= Е Е Е Ъ(г'.,э\к')га(г\з' + з,к' + к) (1.17)

г'=~Мх з'=-Ыу к'=-Ыг

5. Производится корреляция компонент вектора с помощью тензора (1.10).

и г

Оценка для компонент тензора напряжений Рейнольдса на входной границе была предложена в работе [40]. Яц = —(10_4й2), где й — это средняя ско-

О

рость потока. Поскольку турбулентность предполагается изотропной, остальные диагональные элементы тензора совпадают с Яц. Внедиагональные элементы принимаются равными 0.1 от Яц.

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

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

Альтернативный метод наложения пульсаций на усредненное течение был предложен в статьях [42, 43]. Суть его заключается в том, что плоскость наложения располагается не на входной границе, а значительно ниже по потоку, на некотором удалении от обтекаемого тела. Данный подход позволяет разредить сетку вдали от тела и сэкономить машинные ресурсы. Технически это осуществляется с помощью так называемого «актуатора» — плоскости, при прохождении которой поток получает соответствующее ускорение.

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

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

Заключение диссертации по теме «Математическое моделирование, численные методы и комплексы программ», Кудряшов, И. Ю.

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. Реализованный программный комплекс использован при моделировании задач обтекания головной части ракеты-носителя и обтекания крыла с симметричным профилем. Достигнуты следующие результаты:

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

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

Список литературы диссертационного исследования кандидат физико-математических наук Кудряшов, И. Ю., 2013 год

Литература

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

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

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

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

5. Goldstein S. Modern Developments in Fluid Dynamics. Oxford University Press, 1938. Vol. 2. P. 331.

6. Boussinesq J. Theorie de l'Ecoulement Tourbillant.: Tech. Rep. 23: Mem. Presentes par Divers Savants Acad. Sei. Inst. Fr., 1877.

7. Reynolds O. On the Dynamical Theory of Incompressible Viscous Fluids and the Determination of the Criterion // Philosophical Transactions of the Royal Society of London. 1895. Vol. 186. P. 123.

8. Prandtl L. Uber die ausgebildete turbulenz // ZAMM. 1925. Vol. 5. P. 136-139.

9. Von Karman T. Mechanische Ähnlichkeit und Turbulenz // Proc. Int. Congr. Appl. Mech. 1930. P. 85-105.

10. Колмогоров А. Н. Уравнения турбулентного движения несжимаемой жидкости // Изв. АН СССР. 1942. Т. 6, № 1,2. С. 56-58.

11. Rotta J. С. Statistische Theorie Nichthomogener Turbulenz // Zeitschrift fur. Physik. 1951. Vol. 129. P. 547-572.

12. Van Driest E. R. On Turbulent Flow Near a Wall // Journal of the Aeronautical Sciences. 2012. Vol. 23. P. 1007.

13. Smagorinsky J. General Circulation Experiments with the Primitive Equations // Monthly Weather Review. 1963. Vol. 91. P. 99.

14. Spalart P. R. et al. Comments on the feasibility of LES for wing and on a hybrid RANS/LES approach // 1st ASOSR Conference on DNS/LES. 1997.

15. Spalart P. Reflections on RANS Modelling // Progress in Hybrid RANS-LES Modelling. 2010. Vol. Ill of Notes on Numerical Fluid Mechanics and Mul-tidisciplinary Design. P. 7-24.

16. Корнеев В. В. Подход к программированию суперкомпьютеров на базе многоядерных мультитредовых кристаллов // Вычислительные методы и программирование. 2009. Т. 10. С. 123-128.

17. Каляев И. А., Левин И. И. Семейство реконфигурируемых вычислительных систем с высокой реальной производительностью // Вычислительные методы и программирование. 2009. Т. 10.

18. Рейтинг Тор500. URL: http://www.top500.org.

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

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

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

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

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

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

25. 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.

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

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

28. Speziale С. G. Turbulence Modeling in Noninertial Frames of Reference // Theoretical and Computational Fluid Dynamics. 1989. Vol. 1, no. 1. P. 3-19.

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

30. Harten A., Lax P., van Leer В. On upstream differencing and Godunov type methods for hyperbolic conservation laws // SIAM review. 1983. Vol. 25. P. 35-61.

31. Того E. F., Spruce M., Speares W. Restoration of the contact surface in the HLL-Riemann solver // Shock Waves. 1994. Vol. 4. P. 25-34.

32. 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.

33. Breuer M., Jovicic N., Mazaev K. Comparison of DES, RANS and LES for the separated flow around a flat plate at high incidence // Int. J. Numer. Meth. Fluids. 2003. no. 41. P. 357-388.

34. 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.

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

36. Eli Turkel, Vatsa N. V. Local preconditioners for steady and unsteady flow applications // ESAIM: M2AN. 2005. Vol. 39, no. 3. P. 515-535.

37. Jarrin N. et al. Synthetic turbulent inflow conditions for large eddy simulation. // 4th International Turbulence, Heat and Mass Transfer Conference. 2003.

38. Lund T., Wu X., Squires D. Generation of turbulent inflow data for spatially-developing boundary layer simulations. // Journal of Computational Physics. 1998. no. 140. P. 233-258.

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

40. Ferziger J. H., Peric M. Computational Methods for Fluid Dynamics. Third Edition. Berlin, Germany: Springer, 2002.

41. Taylor G. I. The spectrum of turbulence // Proceedings of the Royal Society. 1938. Vol. A164. P. 476-490.

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

43. 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.

44. Gratien J. M. et al. Scalability and load balancing problems in parallel reservoir simulation // Proceedings of 10th ECMOR. 2006. B028.

45. Gustafson J. L. Reevaluating Amdahl's law // Communications of the ACM. 1988. Vol. 31, no. 5. P. 532-533.

46. Amdahl G. Validity of the single processor approach to achieving large-scale computing capabilities // AFIPS Conf. Proc. 1967. P. 483-485.

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

48. Barrett R. et al. Templates for the solution of linear systems: building blocks for iterative methods. Philadelphia: SIAM, 1994.

49. Portland Group. URL: http://www.pgroup.com.

50. Страница вычислительного комплекса K-100. URL: http://www.kiam. ru/MVS/resourses/klOO.html.

51. Spalart P. R., Allmaras S. R. A one-equation turbulence model for aerodynamic flows //La Recherche Aerospatiale. 1994. no. I. P. 5-21.

52. Боровой В. Я. Течение газа и теплообмен в зонах взаимодействия ударных волн с пограничным слоем. Машиностроение, 1983. С. 144.

53. Delery J., Dussauge J.-P. Some physical aspects of shock wave/boundary layer interaction // Shock Waves. 2009. Vol. 19, no. 6. P. 453-468.

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

1998. С. 432.

55. Lovely D., Haimes R. Shock detection from computational fluid dynamics results // 1999 AIAA Computational Fluid Dynamics Conference. No. 99-3288.

1999.-June.

56. Stockdill B. et al. Simulation of unsteady turbulent flow over a stalled airfoil // Computational Fluid Dynamics. 2006. Vol. 14(4), no. 42. P. 359-378.

57. Shur M. et al. Detached-Eddy Simulation of an Airfoil at High Angle of Attack // Engineering Turbulence Modelling and Experiments. 1999. Vol. 13, no. 4. P. 669-678.

58. Торэнбик Э. Проектирование дозвуковых самолетов. Москва: Машиностроение, 1983.

59. Sheldahl R. Е., Klimas Р. С. Aerodynamic Characteristics of Seven Airfoil Sections Through 180 Degrees Angle of Attack for Use in Aerodynamic Analysis of Vertical Axis Wind Turbines: Tech. rep.: Sandia National Laboratories, 1981.

60. Reichardt H. Volstaendige Darstellung der Turbulenten Geschwindigkeitsverteilung in glatten Leitungen // Zeitschrift fuer angewaelende Mathematik und Mechanik. 1957. no. 7. P. 132.

61. McDevitt J. В., Okuno A. F. Static and Dynamic Pressure Measurements on a NACA 0012 Airfoil in the Ames High Reynolds Number Facility: Technical Paper 2485: NASA, 1985.

62. Gilling L., Sorensen N. N., Davidson L. Detached Eddy Simulations of an Airfoil in Turbulent Inflow // 47th AIAA aerospace sciences meeting. No. AIAA 2009-270. 2009.

63. Wilcox D. C. Turbulence Modeling for CFD. DCW Industries, Inc., 1994. ISBN: 0-9636051-0-0.

64. Mann J. Wind Field Simulation // Probabilistic Engineering Mechanics. 1998. Vol. 13, no. 4. P. 269-282.

65. di Mare D. et al. Synthetic turbulence inflow conditions for large eddy simulation // Physics of Fluids. 2006.

66. 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.

67. 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.

68. Crouch J. D. et al. Global Structure of Buffeting Flow on Transonic Airfoils // IUTAM Bookseries / IUTAM Symposium on Unsteady Separated Flows and their Control. Vol. 14. Springer Science+Business Media B.V., 2009.

i