Разработка параллельного неявного метода решения задач динамики вязкого сжимаемого газа тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат наук Борисов, Виталий Евгеньевич

• Борисов, Виталий Евгеньевич
• кандидат наук
• 2016

Оглавление диссертации кандидат наук Борисов, Виталий Евгеньевич

Оглавление

Введение

1. Математическая модель и вычислительные алгоритмы

1.1 Основные уравнения и определяющие соотношения

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

1.3 Дискретизация уравнений

1.3.1 Аппроксимация основных уравнений

1.3.2 Аппроксимация уравнений модели турбулентности

1.3.3 Постановка начальных и граничных условий

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

1.4.1 Метод Ньютона

1.4.2 Метод 1Л"-808

1.4.3 Метод BiCGStab

1.4.4 Общий алгоритм решения

2. Программный комплекс

2.1 Особенности программной реализации

2.2 Валидационные расчёты

2.2.1 Течение вокруг крыла с симметричным профилем

2.2.2 Переход между регулярным и маховским отражением ударных волн в области двойных решений

2.3 Анализ эффективности параллельной реализации

2.4 Выводы

3. Моделирование течений в воздухозаборнике прямоточного воздушно-реактивного двигателя

3.1 Введение

3.2 Описание модели

3.3 Моделирование пограничного слоя заданной толщины

3.4 Расчёты стационарного режима течения в воздухозаборнике

3.4.1 Расчёт с дополнительной областью

в двумерной постановке

3.4.2 Расчёт с заданным пограничным слоем

в двумерной постановке

3.4.3 Расчёт в трёхмерной постановке

3.5 Расчёты нестационарного режима течения в воздухозаборнике

3.5.1 Расчёт с вложением энергии в двумерной постановке

3.5.2 Расчёт с вложением энергии в трёхмерной постановке

3.6 Выводы

Заключение

Приложение

Литература

Введение диссертации (часть автореферата) на тему «Разработка параллельного неявного метода решения задач динамики вязкого сжимаемого газа»

Введение

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

История вычислительной гидрогазодинамики (в западной литературе CFD, Computational Fluid Dynamics) насчитывает чуть более 60-и лет, за которые она прошла значительный путь. За это время были разработаны принципы построения математических моделей, основы теории разностных схем, фундаментальные положения вычислительной математики. Существенный вклад в развитие численных методов для задач гидрогазодинамики внесли научные группы, работающие в Институте прикладной математики им. М.В. Келдыша РАН, Вычислительном центре им. А.А. Дородницына РАН, Институте вычислительных технологий СО РАН, Центральном институте авиационного моторостроения им. П.И. Баранова, Центральном аэрогидродинамическом институте имени профессора Н.Е. Жуковского, Институте теоретической и прикладной механики им. С.А. Христиановича СО РАН, Институте прикладной механики УрО РАН, РФЯЦ-ВНИИЭФ, Московском государственном университете имени М.В. Ломоносова, Санкт-Петербургском политехническом уни-

верситете, и многие другие. Из зарубежных научных центров можно отметить Los Alamos National Laboratory, University of Stanford (США), Imperial College, University of Manchester, University of Cambridge (Великобритания), University of Tolouse, ONERA (Франция) и т.д.

Основными уравнениями гидрогазодинамики являются уравнения Эйлера и Нивье Отокси для несжимаемой или сжимаемой среды. Анализ этих уравнений, способы их замыкания, построение численных методов их решения (в том числе параллельных) представлены в работах О.М. Белоцерковского [1], С.К. Годунова [2], A.A. Самарского, Ю.П. Попова [3], К.И. Бабенко, В.В. Русанова [4], А.Г. Куликовского, Н.В. Погорелова, А.Ю. Семёнова [5], Б.Н. Четверушкина [6], Т.Г. Елизаровой [7], Р.П. Федоренко [8], A.M. Липанова, Ю.Ф. Кисарова, И.Г. Ключникова [9], И.А. Белова, С.А. Исаева [10], P. Lax, В. van Leer, A. Harten, S. Osher [11-14] и многих других. Один из наиболее полных обзоров методов вычислительной гидрогазодинамики можно найти в работе С. Hirsch [15].

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

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

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

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

Один из возможных подходов к построению вычислительно эффективной неявной схемы заключается в использовании метода Ш-БСБ, впервые предложенного в работе [26] и далее модифицированного в работах [27-30]. Метод Ш-БСБ относится к классу безматричных методов решения систем линейных алгебраических уравнений, т.е. решение системы данным методом не требует отдельного хранения коэффициентов матрицы системы целиком, позволяя строить компоненты решения локально по мере прохода циклов по ячейкам сетки через произведения матриц-векторов [31]. В этом методе блочно-диагональная матрица сводится к диагональной за счёт использования аппроксимаций специального вида при линеаризации правой части системы уравнений. Внедиагональные элементы матрицы при этом учитываются в прямом и обратном циклах итерационного процесса Гаусса-Зейделя. Это позволяет существенно сократить объемы используемой памяти и исключить операции с большой матрицей системы, что часто об-

легчает параллельную реализацию методов данного класса [24,32].

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

\ t/ "Л t/ t/

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

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

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

1) Разработка алгоритмов неявной схемы на основе методов LU-SGS и BiCGStab для решения нестационарных осреднённых по Рейнольдсу уравнений Нивье Отокси URANS. В качестве модели турбулентности используется дивергентная форма однопараметрической модели турбулентности Спаларта-Аллмараса для сжимаемых течений с модификацией Эдвардса.

2) Программная реализация разработанных алгоритмов на языках программирования С/С++ с использованием интерфейса параллельного программирования MPI для расчёта трёхмерных турбулентных течений вязкого сжимаемого газа, в том числе течений с энерговложением в поток.

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

4) Проведение расчётов трансзвукового обтекания крыла 0NERA Мб с симметричным профилем. Сравнение результатов расчётов с экспериментом.

5) Проведение расчётов перехода между регулярным и маховским отражением RR^MR ударных волн в области двойных решений, инициируемого вложением энергии в сверхзвуковой поток газа. Сравнение результатов расчётов с известными данными.

6) Моделирование различных режимов работы воздухозаборника прямоточного воздушно-реактивного двигателя. Сравнение результатов расчётов с экспериментом и известными данными.

7) Исследование масштабируемости и эффективности разработанных параллельных алгоритмов.

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

1) Разработаны алгоритмы неявной схемы на основе методов LU-SGS и BiCGStab для решения нестационарных осреднённых по Рейнольд-су уравнений Нииье Стокси URANS совместно с дивергентной формой уравнения однопараметрической модели турбулентности Спаларта-Аллмараса для сжимаемых течений с модификацией Эдвардса и возможностью энерговложения в поток.

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

3) Проведено моделирование стационарных и нестационарных течений в псевдоскачке в тракте воздухозаборника прямоточного воздушно-реактивного двигателя. Показано существенное влияние толщины пограничного слоя на положение и характер ударно-волновой структу-

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

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

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

1) Разработаны алгоритмы неявной схемы на основе методов LU-SGS и BiCGStab для решения нестационарных осреднённых по Рейнольд-су уравнений Нииье Отокси URANS совместно с дивергентной формой уравнения однопараметрической модели турбулентности Опилирти Аллмараса для сжимаемых течений с модификацией Эдвардса и возможностью энерговложения в поток.

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

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

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

Апробация работы. Основные результаты диссертационной работы докладывались на:

• научной конференции «Ломоносовские чтения» (г. Москва, 2014);

• семинарах кафедры вычислительной механики механико-математического факультета МГУ имени М.В. Ломоносова (г. Москва, 2014, 2015, 2016);

ческой и прикладной механики (г. Казань, 2015);

• 9-й Всероссийской школе-семинаре «Аэротермодинамика и физическая механика классических и квантовых систем» (г. Москва, 2015);

•

(г. Москва, 2016);

следований «ICMAR 2016» (г. Пермь, 2016);

• семинаре 11-го отдела ИПМ им. М.В. Келдыша РАН «Вычислитель»

•«

»

2016).

Публикации. Основные результаты диссертационной работы опубликованы в 5 печатных работах [33-37], в том числе в 4 печатных работах [33,34,36,37] в изданиях из перечня ВАК.

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

Структура и объем диссертации. Диссертация состоит из введения, трёх глав, заключения, приложения и списка литературы. Работа представлена на 93 страницах, содержит 35 иллюстраций и 5 таблиц. Список литературы содержит 99 наименований.

Глава 1

Математическая модель и вычислительные

алгоритмы

1.1 Основные уравнения и определяющие соотношения

Для описания трёхмерных течений совершенного вязкого сжимаемого газа используется система нестационарных осреднённых по Рейнольдсу уравнений Навье-Стокса (URANS), в декартовой системе координат (х\ = х7 х2 = У-, х3 = z) имеющая следующий вид:

dq + f = dg_

dt дх4 дх4

+ w,

ал)

где суммирование ведется по ] = 1, 2, 3, Ч _ вектор консервативных переменных, ^ и gj - векторы невязких и вязких потоков соответственно, w - вектор внешних источников (в настоящей работе используется только источник в уравнении баланса энергии):

q

( р ^

рщ

pU'i риз

\рЕ V

/

pUj

\

pUj Щ + öi,j р pUj U'i + 02 ,jP pUj из + Ö3,jP

\ Uj (pE* + p)

, g =

0

TlJ T2j T33

\ nj иг + hj у

w=

0 0 0 0

\W J

Здесь р - плотность, и^ - компоненты вектора скорости и _ компоненты тензора вязких напряжений, ^ - тензор Кронекера, Е* - полная энергия турбулентного течения:

Е * = Е + к, (1.2)

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

j

нию состояния совершенного газа:

1

(Е - 2 ±

р = (7 - 1 )р[Е -->и1] . (1-3)

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

(диг диЛ 2 дщ 2

= МеЧдх; + дх7- з ^ ах^"- зр•

дТ

Ьп = АейР

3 'е11 дх^'

Т

ляется из формулы

Т = р

р Л'

где Л - газовая постоянная, для воздуха Л = 287 Дж/(кг^К).

«Эффективное» значение коэффициентов вязкости и теплопроводности определяется как

Мей = М + Мъ,

и = Ц£ + .

Здесь - коэффициент удельной теплоемкости газа при постоянном давлении, Рг и Ргъ - ламинарное и турбулентное числа Прандтля, для воздуха принимается Рг = 0.7, Ргъ = 0.9. Молекулярная вязкость определяется следующей степенной зависимостью:

Мо

(£ )"•

где в = 0.76 Мо _ молекулярная вязкость при опорной температуре Т).

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

1.2 Модель турбулентности Спаларта^Аллмараса

В работе используется вариант однопараметрической модели турбулентности Спилирти Алл.мирен (SA) [38] для сжимаемых течений [39] с модификацией Эдвардса [40]. В рамках этой модели осредненная величина кинетической энергии турбулентных пульсаций не может быть найдена напрямую, поэтому в уравнении (1.2) полагается Е* = Е. Турбулентная вязкость задается соотношением

з ^

М = pvfvu fvl = 3 , ^з , X = Р -, (1.4)

X + М

где v - модельная величина, которая определяется из основного уравнения модели

+ V • (Pi>u) = p(Pt -D, + T,} + — V • [(М + pv}VV] at a,

+ ^p(Vv)2 -—(}i + pv}Vp^VV. (1.5) a, а,р

Величины Pz и Dz, отвечающие соответственно за производство и диссипацию турбулентности, и T, - за определение ламинарно-турбулентного перехода в пограничном слое, записываются в виде

Сы 4 г ~1 2

Pz = C61(1 — /i2}v, Dz = ( Cw\ fw--р-ft2

, T, = fti(AU }2,

fw = 9

1 + С6

w 3

£6 + ^

w 3

1/6

, g = r + CW2(r6 - r).

Здесь А и - модуль разности между скоростями в потоке и ближайшей точ-

Модификация Эдвардса модели Спилирти Аллмириси [40] предназначена для улучшения поведения решения в пограничном слое. Соответству-

ющие величины £ и г записываются в виде:

£ = л/Я

1 .

- + л1

.X

5 =

/ дщ диЛ дщ 2 / \

V дх* дх{) дxj 3 V дхк )

%,2 = 1 47 / о \ /

г = 1апЬ

V

к ^/С2 /

) / 1апЬ(1.0).

Остальные величины являются константами модели БА и представлены в таб. 1.1.

2

^Ь'р Сы Съ2 СVI С-ш1 Ст2 Сы3

2 3 0.41 0.1335 0.622 7.1 П , = СЬ1 1 1+^Ь2 ^■Ю1 = к2 + а 0.3 2.0

Таблица 1.1. Константы модели Б А

При моделировании полностью турбулентного пограничного слоя учёт /и и Л2 не вносит существенных изменений в решение, поэтому обычно ими пренебрегают [39].

1.3 Дискретизация уравнений

Уравнения (1.1) и (1.5) дискретизируются по пространственным переменным методом конечных объемов [41]. Этот метод имеет ряд преимуществ перед альтернативным методом конечных разностей, так как применим к любой геометрии, оперирует любыми сетками и позволяет избежать проблем с метрическими сингулярностями обобщённых координат. Предполагая, что расчётная область покрыта сеткой, состоящей из неперекрывающихся многогранных ячеек, метод конечных объемов реализуется интегрированием уравнений (1.1) и (1.5) по каждой счётной ячейке с последующим преобразованием объемных интегралов от потоков в интегралы по граням ячейки.

Для дискретизации уравнений по времени используются явная и неявная схемы Эйлера [41].

1.3.1 Аппроксимация основных уравнений

После аппроксимации по пространству методом конечных объемов полудискретный вариант системы уравнений Навье-Стокса (1.1) имеет вид:

Уг ^ = - ^ (Ъ - ga) +

(1.6)

^ = ({к ,па), ga = ,па),

где г = 1... Щ N - число счётных ячеек, ^ - объем ячейки, ва - площадь боковой грани ячейки, па - единичный вектор внешней нормали к грани ячейки. Нижний индекс а указывает на грань, суммирование ведется во всем граням, составляющим боковую поверхность рассматриваемой счётной ячейки г.

Далее па каждой грани а вводится ортонормированный базис (п,к,1), где к = (к1,к2,к3) и I = (11,12,13) - касательные векторы. Матрица преобразования Та к локальному базису при этом определяется координатами базисных ортов:

/10 0 0 0\

0 п1 п2 п3 0

та = 0 к1 к2 к3 0

0 11 12 13 0

0 0 0 0 1

Невязкий поток на грани ячейки представляется через локально-одномерный поток Еа в нормальном для грани направлении:

Ъ = Т-1 ¥а, ¥а =

Здесь поток f1 имеет такое же определение, как в уравнении (1.1), но аргументом его являются контрвариантные компоненты вектора q в локальном

базисе (п, к, I) боковой грани ячейки.

В настоящей работе для вычисления невязкого потока Еа на грани счётной ячейки используется обобщенная схема Годунова [2] с интерполяционной схемой УЕКГОЗ [42]. Детальное описание данных методов представлено ниже.

Аппроксимация вязкого потока gcr па грани счётной ячейки в правой части уравнения (1.6) вычисляется интерполяцией сеточных значений и их производных:

gа = g (да, ) ,

_ УгЧа(г) + Уа(г)Чг = ^Ча^г) + Уа(г) Уч,

Ча — Т7 , Тг , V Ча —

К + К(г) ^ К + У а (г)

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

Обозначим в качестве Иг правую часть уравнения (1.6), во введённых обозначениях имеющую вид:

Иг _ - ^ {Т-1¥а - gа) «а +

В конечном итоге после дискретизации по времени явным методом Эйлера система уравнений (1.6) примет вид:

„п+1 Чп

КЧ-ААГЧ- _ В?. (1-7)

Аналогично проводится дискретизация неявным методом Эйлера:

КЧ Аг Ч _ В?+1. (1.8)

Обобщенная схема Годунова решения задачи Римана

Для того чтобы определить численное значение невязкого потока Еа па грани счётной ячейки, используется обобщенная схема Годунова [2]. В данной схеме поток вычисляется на решении задачи Римана о распаде про-

извольного разрыва на грани ячейки между состояниями «слева» q^ и

TD

«справа» q^, которые вычисляются по той или иной интерполяционной схеме, т.е. Fa = Fa (q^, q^).

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

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

Схема реконструкции WEN03

В настоящей работе для реконструкции значений q^ и q^ слева и справа от расчётной грани а используется интерполяционная схема WEN03 [43-47], адаптированная для структурированных неравномерных сеток [42]. Данный метод проводит реконструкцию величин вдоль одного координатного направления, при этом значения по каждой координате вычисляются по очереди.

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

Пусть грань а = i +1 /2 разделяет ячейки с индексами гиг +1 как показано на рис. 1.1. Зеленым и синим цветами обозначены линейные размеры ячеек в рассчитываемом направлении, используемые для реконструкции значений слева и справа от грани соответственно. Под Xj будем подразуме-

вать координаты центра г-й ячейки, под ха - координаты центра грани а.

Рис. 1.1. Интерполяция переменных на границе ячеек

Рассмотрим процедуру построения значения согласно схеме

ЫЕКТОЗ [42]. Заметим, что реконструкция каждой компоненты^, ] — 1, 5 вектора ч проводится отдельно от остальных. Интерполяционные полиномы Р^ (х) и Р^ (х) для компоненты д^ будут иметь вид

*<х)_<*)г + (х - хг)

РЛх)_(Ч3)„ + 2(1 -('<')г (х - хг)

Соответствующие индикаторы гладкости строятся как

(л {(Яз)г - (Оз)г-1

(1 + р)

ь _ л, + оч2 {(Яз)г - (Яз)г-0 ,

1г ((Ъ 1 - (Ъ )г )2 .

(1 + )2 г 1 г

Оптимальные линейные веса имеют вид:

ЮI юг

(¿1 = -——:г-, = -——:г-,

Щ + иг и I + и г

СО 1 = -тг , Юг = -тг ,

(е + II)2' (е + /г)2'

Г + = 1 Г+= 1+Р 1 1 + Р + 7' г 1 + Р + 7'

где £ = 0.5 к2 - зависящий от шага сетки регуляризующий параметр, предохраняющий от деления на ноль в области равномерного течения.

После вычисления указанных величин реконструированное значение (^ строится по формуле

(<Ц )Ьа = и1 Р^Ха) + Рг1(ха).

Для вычисления компонент вектора q? используются аналогичные формулы с заменой соответствующих величин к, ¡3 и 7:

С01 Сиг С} Сг

Ш1 = -, иг = -, ю [ = -2 ' = -2 '

и I + (Ог + юг (е + 11 )2 (е + 1Г )2

СГ =1+ + + 7' 11 = (ТГ^ ^)г+1 - ^)г)2 '

= 1 + ^ + 7' 1г = (1+^ )г+2 - (Ъ'

? м = (* )„(х - х,1:

^(х) = (ъ)+1 + 2(Ц^^^ (Х - хт) ' )? = Р^(Ха )+ Р? (Ха ).

1.3.2 Аппроксимация уравнений модели турбулентности

Дискретизация уравнения модели турбулентности Опилирти Аллмараса (1.5) проводится аналогично дискретизации основных уравнений. Перепишем его в следующем виде (без учёта величин /ц

и УУ:

^ = -V • (рип) + -V • [(М + ри(1 + Сю)) Vй] - ^йУ • [рУй] ел

+ Аф + А2й2 -—(ц + рй )Ур^Уй, (1.9) а^р

где А\ = рСы§, А2 = -рСШ1 ¡и}/(I2.

Для произвольного £ положим = 0.5 (£ ± |£|). Обозначим в качестве Щ правую часть уравнения (1.9). Тогда после дискретизации уравнения Спаларта-Аллмараса по пространству методом конечных объемов с учётом аппроксимации конвективного потока методом сноса против потока во введённых обозначениях для г-й ячейки она примет вид:

Щ = - ^ + и-А1)8аРа + Г^ (еепУ)а

а

+ У (А\Щ + А2й2)--— (цг + ргйг)Ург • V

а„рг

Г^ = —[ра + (1 + СЪ2)райа - СЪ2 Ра^г] -

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

В конечном итоге после дискретизации по времени явным методом Эйлера уравнение модели турбулентности Опилирти Алмириси примет вид:

й'п+1 _ -п

РгУг ^ = Щ- (1.Ю)

Аналогично проводится дискретизация неявным методом Эйлера:

~П+1 _ ~п

Р,У, = Щ+1. (1.11)

1.3.3 Постановка начальных и граничных условий

Уравнения (1.7), (1.8), (1.10), (1.11) должны быть дополнены соответствующими граничными и начальными условиями, конкретный вид которых зависит от задачи. Для постановки граничных условий в настоящей работе применена схема с использованием двух (согласно шаблону схемы УЕКГОЗ) рядов фиктивных ячеек за каждой границей расчётной области. При этом допускаются следующие типы границ:

• граница типа «вход»: задаются все 5 значений вектора примитивных переменных: 3 компоненты вектора скорости и, давление р и плотность р; модельная переменная V = С^/р, где для полностью турбулентного пограничного слоя полагается С = 1 ^ 5;

• граница типа «выход»: значения всех переменных экстраполируются в фиктивные ячейки изнутри расчётной области;

• граница типа «твердая стенка» и условия симметрии: для твердой стенки со скольжением потока и условия симметрии в фиктивных ячейках значения давления р, плотности р и касательной компоненты скорости задаются симметричным отображением значений во внутренних точках относительно границы, а значение нормальной компоненты скорости и модельной переменной V - антисимметричным отображением, т.е. знак меняется на обратный. В случае твердой стенки без скольжения потока (т.е. условия прилипания) в отличие от преды-

Список литературы диссертационного исследования кандидат наук Борисов, Виталий Евгеньевич, 2016 год

Литература

[1] Белоцерковский О.М. Численное моделирование в механике сплошных сред. М.: Физматлит, 1994, 444 с.

[2] Годунов С.К., Забродин A.B., Иванов М.Я., Крайко А.Н., Прокопов Г.П. Численное решение многомерных задач газовой динамики. М.: Наука, 1976, 400 с.

[3] Самарский A.A., Попов Ю.П. Разностные методы решения задач газовой динамики. М.: Наука, 1992, 424 с.

[4] Бабенко К.П., Воскресенский Т.П., Любимов А.Н., Русанов В.В. Пространственное обтекание гладких тел идеальным газом. М.: Наука, 1964, 505 с.

[5] Куликовский А.Г., Погорелов Н.В., Семёнов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001, 608 с.

[6] Четверушкин Б.Н. Кинетически-согласованные схемы в газовой динамике. М.: Изд-во МГУ, 1999, 232 с.

[7] Елизарова Т.Г. Математические модели и численные методы в механике жидкости и газа. M.: Изд-во МГУ, 2005, 224 с.

[8] Федоренко Р.П. Введение в вычислительную физику. 2-е изд, испр. и доп. Долгопрудный: Издательский дом «Интеллект», 2008, 503 с.

[9] Липанов A.M., Кисаров Ю.Ф., Ключников И.Г. Численный эксперимент в классической гидромеханике турбулентных потоков. Екатеринбург: Изд-во УрО РАН, 2001, 160 с.

[10] Белов И.А., Исаев С.А. Моделирование турбулентных течений. СПб: Изд-во БГТУ, 2001, 108 с.

[11] Harten A., Lax P., van Leer B. On upstream differencing and Godunov type methods for hyperbolic conservation laws / / SI AM review. 1983. V. 25. P. 35-61.

[12] Harten A., Osher S. Uniformly high-order accurate non-oscillatory schemes // Int. J. of Numerical Analysis. 1987.

[13] Lax P.D., Wendroff B. Difference schemes for hyperbolic equations with high order of accuracy // Comm. Pure Appl. Math. 1964. V. 17, P. 381-398.

[14] Lax P.D. Hyperbolic systems of conservation laws and mathematical theory of shock waves. I. Philadelphia: SIAM, 1973.

[15] Hirsch C. Numerical computation of internal and external flows. Vol 1. Elsevier, 2nd edition, 2007, 680 p.

[16] Beam R.M., Warming R.F. An implicit factored scheme for the compressible Navier-Stokes equations // AIAA Journal. 1978. V. 16. P. 393-402.

[17] Ковеня B.M., Тарнавский Г.А., Яненко H.H. Неявная разностная схема для численного решения пространственных уравнений газовой динамики // Журн. вычисл. математики и матем. физики. 1980. Т. 20. № 6. С. 1466-1482.

[18] Ковеня В.М., Тарнавский Г.А., Чёрный С.Г. Применение метода расщепления в задачах аэродинамики. Новосибирск: Наука, 1990, 247 с.

[19] Haga Takanori, Sawada Keisuke and Wang Z.J. An Implicit LU-SGS Scheme for the Spectral Volume Method on Unstructured Tetrahedral Grids. Commun. Com/put. Phys. 2009. V. 6. № 5. P. 978-996.

[20] Balay S., Gropp W.D., Mclnnes L.C., Smith B.F. PETSc 2.0 users manual. Technical Report ANL-95/11- Revision 2.0.24. Argonne National Laboratory, 1999.

[21] Balay S., Buschelman K., Gropp W.D., Kaushik D., Knepley M.G., Mclnnes L.C., Smith B.F., Zhang H. PETSc Web page. 2001.

http://www.mcs.anl.gov/petsc

[22] Balay S., Gropp W.D., Mclnnes L.C., Smith B.F. Efficiency management of parallelism in object oriented numerical software libraries // Modern Software Tools in Scientific Commuting. Birkhauser Press, 1997. P. 163202.

[23] Rasetarinera P., Hussaini M.Y. An efficient implicit discontinuous Galerkin method // Journal of Computational Physics. 2001. V. 172. № 2. P. 718-738.

[24] Luo H., Sharov D., Baum J.D., Lohner R. Parallel unstructured grid GMRES+LU-SGS method for turbulent flows. American Institute of Aeronautics and Astronautics, 2003.

[25] Копьев В.Ф., Титарев В.А., Беляев И.В. Разработка методологии расчёта шума винтов с использованием суперкомпьютеров // Ученые записки ЦАРИ. 2014. Т. XLV. № 2. С. 78-106.

[26] Yoon S., Jameson A. Lower-upper symmetric-Gauss-Seidel method for the Euler and Navier-Stokes equations // AIAA Journal. 1988. V. 26. № 9. P. 1025-1026.

[27] Chen R.F., Wang Z.J. Fast, block lower-upper symmetric Gauss-Seidel scheme for arbitrary grids // AIAA Journal. 2000. V. 38. № 12. P. 22382245.

[28] Jameson A., Caughey D.A. How many steps are required to solve the Euler equations of steady compressible flow: in search of a fast solution algorithm. In: 15th AIAA computational fluid dynamics conference, Anaheim, CA. June 2001.

[29] Li Dong, Men'shov Igor, Nakamura Yoshiaki. Detached-Eddy Simulation of Three Airfoils with Different Stall Onset Mechanisms // Journal of Aircraft. 2006. V. 43. № 4. P. 1014-1021.

[30] Sun Y., Wang Z.J., Liu Y., Chen C.L. Efficient implicit LU-SGS algorithm for high-order spectral difference method on unstructured hexahedral grids. AIAA Paper № 2007-0313, 2007.

[31] Brown P.N., Hindmarsh A.C. Matrix-Free Methods for Stiff Systems of ODE's // SI AM Journal on Numerical Analysis. 1986. V. 23. № 3. P. 610-638.

[32] Павлухин П.В., Меньшов И.С. Эффективная реализация метода LU-SGS для задач газовой динамики // Научный вестник МГТУ ГА. 2011. Т. 165. С. 46-55.

[33] Борисов В.Е., Давыдов А.А., Кудряшов И.Ю., Луцкий А.Е., Меньшов И.С. Параллельная реализация неявной схемы на основе метода LU-SGS для моделирования трёхмерных турбулентных течений // Математическое моделирование. 2014. Т.26. № 10. С. 64-78.

[34] Борисов В.Е. Моделирование обтекания крыла ONERA Мб с помощью параллельной реализации неявной схемы // Вестник МГУ. Серия 1. Математика. Механика. 2015. № 4. С. 65-68.

[35] Борисов В.Е., Луцкий А.Е. Моделирование течений в воздухозаборнике ПВРД // Физико-химическая кинетика в газовой динамике. 2015. Т. 16, вып. 1. 12 с. http://chemphys.edu.ru/issues/ 2015-16-1/articles/530/

[36] Борисов В.Е., Кудряшов И.Ю., Луцкий А.Е. Численное исследование формирования псевдоскачка в канале // Препринты ИПМ им. М.В. Келдыша. 2016. № 2. 24 с.

http://library.keldysh.ru/preprint.asp?id=2016-2

[37] Борисов B.E., Луцкий А.Е. Моделирование перехода между регулярным и маховским отражением ударных волн с помощью неявной схемы на основе методов LU-SGS и BiCGStab // Препринты ИПМ им. М.В. Келдыша. 2016. № 68. 36 с.

http://library.keldysh.ru/preprint.asp?id=2016-68

[38] Spalart P.R., Allmaras S.R. A One-equations Turbulence Model for Aerodynamics Flows. AIAA Paper № 92-0439, 1992.

[39] Allmaras S.R., Johnson F.T., Spalart P.R. Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model. In: 7th International Conferecnce on CFD, Big Island, Hawaii, July 2012.

[40] 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.

[41] Blazek J. Computational Fluid Dynamics: Principles and Applications. New York: Elsevier, 2001.

[42] Cravero I., Semplice M. On the accuracy of WENO and CWENO reconstructions of third order on nonuniform meshes // Journal of Scientific Computing. 2016. V. 67. № 3. P. 1219-1246.

[43] Liu X.D., Osher S., Chan T. Weighted essentially non-oscillatory schemes // Journal of Computational Physics. 1994. V. 115. P. 200212.

[44] Shu C.W. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. Chapter in: Advanced numerical approximation of nonlinear hyperbolic equations. Volume 1697 of the series Lecture Notes in Mathematics, p. 325-432. Springer, Berlin, 1998.

[45] Hu C., Shu C.W. Weighted essentially non-oscillatory schemes on triangular meshes // Journal of Computational Physics. 1999. V. 150. P. 97-127.

[46] Shu C.W. High order weighted essentially non-oscillatory schemes for convection dominated problems / / SI AM Review. 2009. V. 51. P. 82-12.

[47] Arandiga F.. Baeza A., Belda A.M., Mulet P. Analysis of WENO schemes for full and global accuracy // SIAM J. Numer. Anal. 2011. V. 49. № 2. P. 893-915.

[48] Jameson A. Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings. AIAA paper № 911596, 1991.

[49] Рябенький B.C. Введение в вычислительную математику: Учеб. пособие. 2-е изд., исправ. М.: ФИЗМАЛИТ, 2000, 296 с.

[50] Русанов В.В. Расчёт нестационарных ударных волн с препятствиями // ЖВМ и МФ. 1961. Т. 1. № 2. С. 267-279.

[51] Men'shov I., Nakamura Y. On implicit Godunov's method with exactly linearized numerical flux // Computers and Fluids. 2000. V. 29. P. 595616.

[52] Men'shov I., Nakamura Y. An implicit advection upwind splitting scheme for hypersonic air flows in thermochemical nonequilibrium. Collection of technical papers of 6th Int. Symp. on CFD, Lake Tahoe, Nevada, 1995. V. 2. P. 815-821.

[53] Saad Y. Iterative Methods for Sparse Linear Systems. Second Edition. SIAM, 2003.

[54] Knoll D.A., Keyes D.E. Jacobian-free Newton-Krylov methods: a survey of approaches and applications // Journal of Computational Physics. 2004. V. 193. № 2. P. 357-397.

[55] Krasnopolsky B. The reordered bicgstab method for distributed memory computer systems // Proceedia Computer Science. 2010. V. 1. P. 213218.

[56] MPI: A message passing interface standard. University of Tennessee, Knoxville, Tennessee: 1994, 228 p.

[57] Chaco: Software for Partitioning Graphs.

https://cfwebprod.sandia.gov/cfdocs/CompResearch/ templates/insert/softwre.cfm?sw=36

[58] Cuthill E., McKee J. Reducing the Bandwidth of Sparse Symmetrical Matrices // Proc. ACM Nat. Conf. New York, 1969. P. 157-172.

[59] Волков К.H., Емельянов В.H. Вычислительные технологии в задачах механики жидкости и газа. М.: ФИЗМАЛИТ, 2012. 468 с.

[60] Гергель В.П. Теория и практика параллельных вычислений: учебное пособие. М.: Бином. Лаборатория знаний; Интернет-университет открытых технологий, 2007. 423 с.

[61] Семёнов И.В., Ахмедьянов И.Ф., Уткин П.С. Разработка вычислительного комплекса для решения двух и трёхмерных задач газодинамики реагирующих течений на многопроцессорных ЭВМ // Высокопроизводительные параллельные вычисления на кластерных системах: Мат-лы 6 Междупарод, науч.-практ. семинара. 12-17 декабря 2006 г. СПб.: Изд-во Санкт-Петербургского гос.университета, 2007. Т. 2. С. 138-145.

[62] Семёнов И.В., Ахмедьянов И.Ф. Разработка параллельного алгоритма LU-SGS для решения многомерных задач вычислительной газодинамики // Материалы Четвертой Сибирской школы-семинара по параллельным и высокопроизводительным вычислениям Изд-во Томск: Дельтаплан, 2008. С. 122-129.

[63] Intel® Compilers.

http ://software.intel.com/ru-ru/intel-compilers

[64] Гибридный вычислительный кластер К-100. http ://www.kiam.ru/MVS/resourses/klOO.html

[65] NPARC Alliance Validation Archive.

https://www.grс.nasa.gov/WWW/wind/valid/m6wing/m6wing.html

[66] Schmitt V., Charpin F. Pressure Distributions on the ONERA-M6-Wing at Transonic Mach Numbers. Experimental Data Base for Computer Program Assessment. Report of the Fluid Dynamics Panel Working Group 04, AGARD AR 138, 1979.

[67] Gaultier Sylvain, ONERA research lab.

http://www.onera.fг/en/actualites/image-du-mois/ onera-m6-wing-star-cfd

[68] Кудрявцев A.H. Вычислительная аэродинамика сверхзвуковых течений с сильными ударными волнами: дис. ... д-ра ф.-м. наук. ИТПМ СО РАН, Новосибирск, 2014.

[69] Иванов M.C., Кудрявцев А.Н., Никифоров С.В., Хотяновский Д.В. Переход между регулярным и маховским отражением ударных волн: новые численные и экспериментальные результаты // Аэромеханика и газовая динамика. 2002. № 3, С. 3-12.

[70] Ivanov M.S., Kudryavtsev A.N., Nikiforov S.B., Khotyanovsky D.V., Pavlov A.A. Experiments on shock wave reflection transition and hysteresis in low-noise wind tunnel // Physics of Fluids. 2003. V. 15. P. 1807-1810.

[71] Khotyanovsky D.V., Knight D.D., Kudryavtsev A.N., Ivanov M.S. Numerical study on laser-induced shock wave reflection transition. In: The 5th Int. Workshop on Shock-Vortex Interaction, Kaohsiung, Taiwan, Oct. 27-31, 2003. P. 84-90.

[72] Khotyanovsky D.V., Kudryavtsev A.N., Ivanov M.S. Effects of a Singlepulse Energy Deposition on Steady Shock Wave Reflection // Shock Waves. 2006. V. 15. № 5. P. 353-362.

[73] Бормотова Т.А., Володин В.В., Голуб В.В., Ласкин И.Н. Тепловая коррекция входного диффузора гиперзвукового прямоточного воздушно-реактивного двигателя // Теплофизика высоких температур. 2003. Т. 41. № 3. С. 472-477.

[74] Замураев В.П., Калинина А.П. Численно-аналитическое моделирование структуры сверхзвукового течения газа в канале переменного сечения с подводом энергии // Инженерно-физический журнал. 2015. Т. 88. № 1. С. 210-219.

[75] Латыпов А.Ф. Численное моделирование течения в канале переменной площади сечения при импульсно-периодическом подводе энергии // Прикладная механика и техническая физика,. 2009. Т. 50. № 1. С. 3-11.

[76] Третьяков П.К. Организация пульсирующего режима горения высокоскоростных ПВРД // Физика горения и взрыва. 2012. Т. 48. № 6. С. 21-27.

[77] Забайкин В.А. Управление псевдоскачком нестационарным воздействием // Физико-химическая кинетика в газовой динамике. 2011. Т. 12. 7 с. http://chemphys.edu.ru/issues/2011-12/articles/ 353/

[78] Quirk J.J. AM RITA - A computational facility (for CFD modeling) // VKI 29th CFD Series. 1998. P. 23-27.

[79] Hadiadi A., Kudryavtsev A. Computation and flow visualization in highspeed aerodynamics // Journal of Turbulence. 2005. V. 6. № 16.

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

[81] Гуськов О.В., Копчёнов В.И., Липатов И.И., Острась В.Н., Старухин В.П. Процессы торможения сверхзвуковых течений в каналах. М.: ФИЗМАЛИТ, 2008.

[82] Anderson J. Hypersonic and high-temperature gas dynamics, Second edition. AIAA, 2006.

[83] Andreadis D. Scramjet Engines Enabling the Seamless Integration of

Air & Space Operations // The Industrial Physicist. August/September 2004.

[84] Mutzman R.C., Murphy J.S. X-51 Development: A Chief Engineer's Perspective. In: 17th AIAA International Space Planes and Hypersonic Systems and Technologies Conference, 13 April 2011.

[85] Peebles C. Road to Mach 10: Lessons learned from the X-34A flight research program. AIAA, 2008.

[86] Железнякова А.Л., Суржиков С.Т. На пути к созданию модели виртуального ГЛА. I. М.: НПМех РАН, 2013, 160 с.

[87] Быков Л.В., Молчанов A.M., Щербаков М.А., Янышев Д.С. Вычислительная механика сплошных сред в задачах авиационной и космической техники. М.: ЛЕНАНД, 2015, 688 с.

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

[89] Кудряшов И.Ю., Луцкий А.Е. Численное моделирование эффектов турбулизации и реламинаризации потока в результате активных внешних воздействий // Математическое моделирование. 2014. Т. 26. № 3. С. 3-13.

[90] Menter F.R., Kuntz М., Langtry R. Ten Years of Industrial Experience with the SST Turbulence Model // Turbulence, Heat and Mass Transfer 4: Proceedings of the Fourth International Symposium on Turbulence, Heat and Mass. Begell House, Inc., 2003. P. 625-632.

[91] Wagner J.L. Experimental studies of unstart dynamics in inlet/isolator configurations in a Mach 5 flow. PhD Thesis. The University of Texas at Austin, 2009.

[92] Koo H., Raman V. Large-eddy simulation of a supersonic inlet-isolator // AIAA Journal. 2012. V. 50. № 7. P. 1596-1613.

[93] Xu S., Martin M.P. Assessment of inflow boundary conditions for compressible turbulent boundary layers // Physics of Fluids. 2004. V. 16. № 7. P. 2623-2639.

[94] Do H., Im S., Mungal M.G. et al. The influence of boundary layers on supersonic inlet flow unstart induced by mass injection // Experiments in Fluids. 2011. V. 51. № 3. P. 679-691.

[95] Emami S., Trexler C.A., Auslender A.H., Weidner J.P. Experimental Investigation of Inlet-Combustor Isolators for a Dual-Mode. Technical Report 3502. NASA Langley Technical Report Server, 1994.

[96] Wilcox D.C. Turbulence Modeling for CFD (Third Edition). DCW Industries, Inc., 2006.

[97] Jang I., Pecnik R., Moin P. A numerical study of the unstart event in an inlet/isolator model // Center for Turbulence Research Annual Research Briefs. 2010. P. 93-103.

[98] Жуков В.Т. [и др.] Исследование картины течения в модельном тракте двигателя высокоскоростного летательного аппарата // Препринты ИПМ им. М.В. Келдыша. 2015. № 5. 23 с.

[99] Топлива, смазочные материалы, технические жидкости. Ассортимент и применение: Справочник. 2-е изд. Под ред. В.М. Школьни-кова. М.: Химия, 1999.