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

Московский государственный университет имени М.В. Ломоносова Факультет вычислительной математики и кибернетики

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

Исаков Виктор Александрович

Квазиакустическая схема для уравнений Эйлера газовой динамики

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

АВТОРЕФЕРАТ

диссертации на соискание учёной степени кандидата физико-математических наук

14 НОЯ 2013

005538449

Москва — 2013

005538449

Работа выполнена в федеральном государственном образовательном учреждении высшего профессионального образования «Московский Государственный Университет имени М.В. Ломоносова»

Научный руководитель:

Официальные оппоненты:

Ведущая организация:

Фаворский Антон Павлович ,

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

Головизнин Василий Михайлович,

доктор физико-математических наук, профессор кафедры вычислительных методов, федеральное государственное образовательное учреждение высшего профессионального образования «Московский Государственный Университет имени М.В. Ломоносова» Жуков Виктор Тимофеевич, доктор физико-математических наук, заведующий отделом, федеральное государственное бюджетное учреждение науки «Институт прикладной математики им. М.В. Келдыша Российской академии наук» Котеров Владимир Николаевич, кандидат физико-математических наук, ведущий научный сотрудник, федеральное государственное бюджетное учреждение науки «Вычислительный центр им. A.A. Дородницына Российской академии наук» федеральное государственное бюджетное учреждение науки «Институт автоматизации проектирования Российской академии наук»

Защита состоится «4» декабря 2013 г. в 15 час. 30 мин. на заседании диссертационного совета Д 501.001.43 при Московском Государственном Университете имени М.В. Ломоносова по адресу: 119991, ГСП-1, Москва, Ленинские горы, МГУ, 2-й учебный корпус, факультет ВМК, аудитория 685.

С диссертацией можно ознакомиться в научной библиотеке Московского Государственного Университета имени М.В. Ломоносова по адресу: 119192, Москва, Ломоносовский проспект, д. 27.

Автореферат разослан *_/_> 2013 г.

Ученый секретарь диссертационного

совета Д 501.001.43 доктор физико - ма- . Е.В. Захаров

тематических наук, профессор

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

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

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

В работах [1,2] описывается предложенный А.П. Фаворским новый алгоритм решения линейного и простейшего квазилинейного уравнения переноса. В основе алгоритма лежит кусочно-линейная реконструкция численного решения, опирающаяся на характеристические свойства уравнений переноса. Для случая линейного уравнения переноса проведено достаточно подробное исследование предлагаемой методики. Показано, что построенная на основе данной методики явная, консервативная, однородная разностная схема обладает вторым порядком аппроксимации в областях гладкости решения и позволяет обеспечить монотонность решения без использования искусственных регуляриз аторов.

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

Целью диссертационной работы является

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

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

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

1. разработан алгоритм решения уравнений газовой динамики в одномерном случае;

2. проведено обобщение алгоритма на случай двух и трёх пространственных измерений;

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

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

5. разработан параллельный алгоритм, с использованием которого проведено математическое моделирование прикладной задачи о взаимодействии газовых струй, истекающих из двигателей самолёта с отражающим экраном-отбойником на суперкомпьютере IBM eServer pSeries 690 (Regatta).

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

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

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

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

Основные положения, выносимые на защиту:

1. алгоритм для решения уравнений газовой динамики в одномерном случае;

2. обобщение предлагаемого алгоритма на случай двух и трёх пространственных измерений;

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

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

1. научной конференции „Тихоновские чтения" (г. Москва, 2009);

2. международной молодёжной конференции-школе „Современные проблемы прикладной математики и информатики" (МРАМСЭ' 2012) (г. Дубна, 22 - 27 августа 2012);

3. ХХ-ой международной научной конференции студентов, аспирантов и молодых учёных ,Ломоносов - 2013" , секция „Вычислительная математика и информатика" (г. Москва, 9-12 апреля 2013);

4. научной конференции ,Домоносовские чтения" (г. Москва, 24 апреля 2013);

5. научно-исследовательском семинаре в Институте прикладной математике им. М.В. Келдыша РАН под руководством В.Т. Жукова (г. Москва, февраль 2013).

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

Публикации. Основные результаты диссертационной работы опубликованы в 9 печатных работах: 5 статей [1-5] в рецензируемых журналах, рекомендованных ВАК, 4 тезиса докладов [6-9].

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

Диссертационная работа была выполнена при поддержке ФЦП „Научные и научно-педагогические кадры инновационной России" на 2009-2013 годы.

Содержание работы

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

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

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

где / = {р, д — ри, е — р (е + и2/2) }Т - вектор-столбец консервативных переменных системы, ^ (/) = {ир, щ+р,и(е + р)}т - вектор-столбец потоков, р - плотность газа, д = ри - объёмная плотность импульса газа, р - давление

6

газа, е = р/ (7 — 1) + ри2/2 - полная энергия газа единицы объёма (использовано уравнение состояния идеального газа в форме р = (7 — 1) ре, где г -удельная внутренняя энергия газа), и - скорость газа, 7 > 1 - показатель адиабаты. Все газодинамические функции системы зависят от времени Ь и декартовой координаты х.

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

В третьем разделе описан процесс построения квазиакустической схемы. Построение схемы проведено на прямоугольной равномерной сетке, состоящей из пространственно-временных ячеек = [хк-1/2,Хк+1/2\ х [£ги Ьп + Д£], с помощью интегро-интерполяционного метода. Значения газодинамических функций системы определены в центре ячейки £1к = [^-1/2)^+1/2] и представляют собой средние интегральные значения непрерывно распределённых величин по данной ячейке на момент времени 4 = Ьп. Получены основные балансные соотношения, которые выражают законы сохранения массы, импульса и полной энергии

(/Г1 - Ю Ах + [щ+1/2 - Щ_1/3] = О,

где /£, - средние интегральные значения функции / по ячейке на моменты времени г = и Ь = + а величины ~ интегральные

потоки функции Р через границы х = хк±\/2 ячейки

Построение схемы сводится к аппроксимации потоков состоя-

щей из следующих этапов.

На первом шаге в пределах каждой ячейки Пк расчётной сетки каждая из опорных функций системы / = {р, и, р}т на момент времени £ = ¿п приближается линейной функцией

/ (я; хк,Ъ) =г Д" + (х - хк) Щ, х <Е Пк = [хк-1/2,хк+1/2]

Величина задаёт угол наклона линейной функции f{x;xk,tn) и определяется по формуле

If" I I fn I

туп \Jx,k\ m I__l^g.fel tn

k~m+\fik\'Jx* i/ij+i/^i'h*'

гДе flk = (fk ~ fk-i) /btJlk = (Л+i ~ fk) /АТакой способ определения величины позволяет обеспечить монотонность решения на соответствующих участках монотонности сеточной функции {/£}.

Далее на каждом отрезке [zb^fc+i] выбирается горизонтальное сечение, которое назовём общим постоянным фоном функции / между ячейками и fifc+i. Значение опорной функции на фоне определяется как среднее арифметическое значений линейных функций / (x;xk,tn) и f (x;xk+i,tn) на общей границе х — x^+i/2 и обозначается fk+\/2- На общем постоянном фоне определены значения трёх опорных функций системы: рк+1/2,Щ+1/2,Рк+1/2-

На следующем шаге на полуотрезках [хь 2^+1/2] и \xk+i/2,Xk+i] заменяем линейные функции конструкцией, состоящей из общего постоянного фона и расположенных на нём М горизонтальных слоёв разбиения, параллельных координатному направлению х. Построенные горизонтальные слои разбиения отождествляются с малыми возмущениями опорных функций системы. Фоном для распространения каждого такого горизонтального слоя является либо общий постоянный фон, либо поверхность другого горизонтального слоя, примыкающего со стороны общего постоянного фона.

Такой способ приближения функций / (х; Xk, tn) и / (х; Xk+i, tn) приводит к тому, что интегральный поток через границу ячейки складывается из потока, обусловленного фоновыми значениями опорных функций системы, и трёх интегральных компонент потока, соответствующих характеристическим скоростям и,и + с и и — с, с которыми по соответствующим характеристикам распространяются бегущие волны малых возмущений. Бегущие волны представляют собой решение линеаризованных уравнений газовой динамики.

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

ложенных на полуотрезках [я^я^-и/г] и [хк+1/2, '¿к+1}, которые за время Д£, двигаясь каждый со своей скоростью по своему фону, пересекли границу ячейки Пк-

В четвёртом разделе представлены результаты численных расчётов ряда тестовых задач одномерной газовой динамики. Проведено сравнение численного решения, полученного с помощью квазиакустической схемы, с решением, полученным по схеме Годунова первого порядка точности, схемам годуновского типа второго порядка точности: МийСЬ-Напсоск схеме и \¥АЕ схеме (см. рис. 2), а также схеме Роу-Ошера с модификацией Эйнфельдта второго порядка точности (см. рис. 1).

О.ОСЯ

от о.«м 0005 | СШ04 сига

* !

0и«амсоив1с кЬеак -Со&ют тсОюй -Мимх-Нагсоск твЬой -ШАГтс*о<1 - 1 1 :

Рис. 1: Расчёт задачи с гладким начальным профилем скорости.

Рис. 2: Расчёт задачи о распаде разрыва с сильным перепадом давления (с 1000 до 0.01).

Разработанная квазиакустическая схема обладает следующими свойствами:

• схема является явной, консервативной, однородной;

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

• физически корректно передаёт разрывные решения;

• правильно воспроизводит звуковую точку и поведение решения в её окрестности;

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

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

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

где / = {р,дЫ=ри,дЮ=ру,е = р(е+(и2+ь2)/ 2)}Т -

вектор-столбец консервативных переменных системы, Р (/) = {ри,идт+р,щ{2\и(е+р)}Т = {ру, V1', +р, V (е +р)} -

вектор-столбцы потоков.

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

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

Численное решение системы уравнений (2) строится на равномерной сетке, состоящей из пространственно-временных ячеек Пу = [а^-1/2, £¿+1/2] х [2/7-1/2,2/7+1/2] х [¿ти ¿п + имеющих форму прямоугольного параллелепипеда. Значения газодинамических функций системы определены в центре ячейки Пу = [2:^1/2, х!+1/2] х [г/,-1/2,2/7+1/2] и представляют собой средние интегральные значения функций по данной ячейке на момент времени4 =

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

(ЯГ1 - ® ЬхЬу + [щ+1/2 -и?_1/2] + [щ+и2 - Ю]_1/2] = О,

где /у, - средние интегральные значения функции / по ячейке Пу на моменты времени £ = и Ь = величины представляют

собой интегральные потоки ^ и (? через соответствующие границы ячейки Пу.

Процесс аппроксимации интегральных потоков состоит из следующих этапов.

На первом этапе в пределах каждой ячейки расчётной сетки каждая из опорных функций системы / = {р, и, V, р}т на момент времени Ь = приближается линейной плоскостью, проходящей через значение опорной функции в центре ячейки /¿". Уравнение плоскости имеет вид

/ (х, у;хи у, г„) ~ Щ + (х - ян) + (у - Ц^,

Величины и у задают вектор нормали к линейной плоскости и определяются по формуле

!/■- I 1/™ I

7)1 _ ЫЫ1 1П , 1-ЧЛЛ ХП

где £ = {х, у}т, - производная назад и производная вперёд функции

/ вдоль координатного направления £ в узле в момент времени f = Ьп.

Такой способ определения величин £>г™ж и у позволяет обеспечить монотонность решения на соответствующих участках монотонности сеточной функции {/,"} вдоль каждого координатного направления.

Далее производим разбиение каждой ячейки расчётной сетки на х прямоугольных подъячеек. В пределах каждой подъячейки линейная функция f {x,y■,Xi,yj,tn) заменяется средним значением, определённым в центре подъячейки, которое назовём вертикальным столбцом. Данная процедура проводится для каждой опорной функции системы в отдельности.

На следующем шаге для каждой ячейки Пу на расчётной сетке выделяется четыре фоновые подобласти. Каждая фоновая подобласть представляет собой прямоугольник, вершинами которого являются центры рассматриваемой Пу и трёх соседних с ней ячеек. Значение опорной функции на фоновой подобласти /г-+а^+/з, где а, ¡3 = ±1/2 определяется как среднее по объёму, составленному из суммы объёмов вертикальных столбцов из ячеек, пересека-

емых фоновой подобластью. На каждой фоновой подобласти определены значения четырёх опорных функций системы: йi+aj+Jз,i;i+aJ■+/3,pг•+aJ•+lз.

Далее в пределах каждой фоновой подобласти представим линейную функцию f {х,у, XI, уj,tTl) в виде композиции фонового значения /г+ад-+/з и усечённых вертикальных столбцов, высота которых определяется как разность высоты вертикального столбца и фонового значения опорной функции. Каждый усечённый вертикальный столбец, в свою очередь, разделим на Л?ыосЬ блоков одинаковой высоты, которые будем отождествлять с малыми возмущениями опорной функции системы. Будем считать, что фоном для распространения таких малых возмущений является либо значение опорной функции системы на соответствующей фоновой подобласти, либо значение опорной функции на границе, разделяющей данный блок и блок, примыкающий к нему со стороны фонового значения.

Такой способ приближения линейной функции / (х, у; у^ ¿п) в пределах ячейки расчётной сетки приводит к тому, что интегральный поток через границу ячейки Яу складывается из потока, обусловленного значениями опорных функций на фоновых подобластях, пересекающих границу ячейки, и двух интегральных добавок потока, которые соответствуют фоновым подобластям, пересекающим границу ячейки. Каждая интегральная добавка потока состоит из трёх компонент, отвечающих характеристическим скоростям, с которыми вдоль соответствующего координатного направления распространяются одномерные бегущие волны малых возмущений. Бегущие волны являются решением линеаризованных уравнений газовой динамики.

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

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

что численное решение, полученное по квазиакустической схеме, соответству-

12

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

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

Рис. 4: Расчёт цилиндрического взрыва: сравнение численного решения, полученного с помощью квазиакустической схемы, с решением, полученным по ШАР схеме.

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

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

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

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

газовой динамики

df | dFjJ) | dG(f) | dH (/) _ ^ 8t дх dy dz

где / = {р.д^ = ри,д^ = ру,д^ = рю,е = р(е+(и2+ v2+ w2)/2)}Т

- вектор-столбец консервативных переменных системы, а величины — {ри, ид№ ид^2\ид^\и(е +р)}Т,С (/) = {рг>,ид^,-^2' + р,г></3\V(е + р)} , Я(/) = {рш,адд^,гид(2\юд^ + р,ги(е+р)}

- вектор-столбцы потоков.

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

Рис. 5: Расчётная область П.

На границах расчётной области Г2 поставлены следующие условия:

® стенки экрана-отбойника и нижняя граница расчётной области Л считаются непроницаемыми (reflective boundaries).

® на других границах расчётной области Q ставятся условия свободного „вытекания" газа через границу (transmissive boundaries).

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

14

Цель математического моделирования данной задачи состоит в том,

чтобы

» выяснить качественное поведение течения во внутренней и внешней областях экрана-отбойника.

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

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

Рис. 6: Линии тока во внутренней области экран а-отбойника.

Рис. 7: Линии тока с внешней стороны экранаготбойника.

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

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

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

струями, на стенки экрана-отбойника.

15

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

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

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

Разработан параллельный алгоритм, с использованием которого проведено математическое моделирование прикладной задачи на суперкомпьютере IBM eServer pSeries 690 (Regatta, 16 процессоров POWER 4 1,1 GHz). Использование параллельного алгоритма для решения прикладной задачи на сетке размером 42 х 42 х 42 на 16 процессорах позволило сократить время счёта одного шага по времени почти в 16 раз по сравнению с последовательной версией. Расчёт данной задачи также проводился на сетке размером 102 х 102 х 102.

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

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

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

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

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

и Василию Михайловичу Головизнину за И.К. Ермолаеву за любезно предоставленные данные натурных экспериментов, а также С.И. Мухину, М.В. Абакумову, А.Ю. Мокину за ценные советы по подготовке диссертационной работы.

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

[1] Фаворский А.П., Тыглиян М.А., Тюрина H.H., Галанина A.M., Исаков В.А. Численное моделирование распространения акустических импульсов в гемодинамике // Дифф. уравнения. 2009. Т. 45, № 8. С. 1179 - 1187.

[2] Фаворский А.П., Тыглиян М.А., Тюрина H.H., Галанина A.M., Исаков В.А. Численное моделирование распространения гемодинамических импульсов // Матем. моделирование. 2009. Т. 21, № 12. С. 21 - 34.

[3] Абакумов М.В., Галанина A.M., Исаков В.А., Тюрина H.H., Фаворский А.П., Хруленко А.Б. Квазиакустическая схема для уравнений Эйлера газовой динамики // Дифф. уравнения. 2011. Т. 47, № 8. С. 1092 - 1098.

[4] Исаков В.А., Фаворский А.П. Квазиакустическая схема для уравнений Эйлера газовой динамики в случае двух пространственных измерений // Матем. моделирование. 2012. Т. 24, № 12. С. 55 - 59.

Антону Павловичу Фаворскому постоянное внимание к работе,

[5] Галанина A.M., Исаков В.А., Тюрина H.H., Фаворский А.Я.Конструктивный подход к численному решению квазилинейных уравнений переноса // Матем. моделирование. 2013. Т. 25, № 8. С. 80 - 88.

[6] Исаков В.А. Решение уравнений газовой динамики в эйлеровых переменных с использованием сплайн-функций // Сборник тезисов лучших дипломных работ факультета ВМК МГУ имени М.В. Ломоносова. 2010. С. 34 - 36.

[7] Фаворский А.П., Галанина A.M., Исаков В.А. Конструктивный подход к численному решению квазилинейных уравнений переноса // Тезисы международной молодёжной конференции-школы „Современные проблемы прикладной математики и информатики" .г. Дубна. 22 - 27 августа 2012. С. 34.

[8] Исаков В. А. Квазиакустическая схема для уравнений Эйлера газовой динамики // Сборник тезисов XX Международной научной конференции студентов, аспирантов и молодых учёных „Ломоносов - 2013". Москва. МГУ имени М.В. Ломоносова. Ф-т ВМК. 9-12 апреля 2013. Секция „Численные методы и математическое моделирование" . С. 134 - 135.

[9] Фаворский А.П., Исаков В.А. Квазиакустическая схема для уравнений Эйлера газовой динамики // Сборник тезисов конференции „Ломоносовские чтения" 2013. Москва. МГУ имени М.В. Ломоносова. Ф-т ВМК. 24 апреля 2013. Секция XI. С. 87 - 88.

Напечатано с готового оригинал-макета

Подписано в печать 30.10.2013 г. Формат 60x90 1/16. Усл.печ.л. 1,0. Тираж 100 экз. Заказ 343.

Издательство ООО "МАКС Пресс" Лицензия ИД N 00510 от 01.12.99 г. 119992, ГСП-2, Москва, Ленинские горы, МГУ им. М.В. Ломоносова, 2-й учебный корпус, 527 к. Тел. 8(495)939-3890/91. Тел./факс 8(495)939-3891.

Московский государственный университет имени М.В. Ломоносова Факультет вычислительной математики и кибернетики

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

04201364699

Исаков Виктор Александрович

Квазиакустическая схема для уравнений Эйлера газовой динамики

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

программ

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

Научный руководитель:

д. ф. - м.н., проф. А. П. Фаворский д.ф. - м.н., проф. В.М. Головизнин

Москва — 2013

Оглавление

Введение 5

1 Квазиакустическая схема для уравнений Эйлера газовой динамики 20

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

1.2 Линеаризация уравнений газовой динамики..........................21

1.2.1 Решение линеаризованной системы уравнений..............23

1.3 Построение разностной схемы..........................................28

1.4 Аппроксимация интегральных потоков ..............................33

1.4.1 Линейная реконструкция опорных функций системы ... 34

1.4.2 Выбор фона и определение фоновых значений опорных функций системы................................................35

1.4.3 Расслоение (разбиение) опорных функций системы .... 37

1.4.4 Аппроксимация интегральных потоков ......................39

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

1.5 Результаты численных расчётов......................................45

1.5.1 Расчёт задачи с гладким начальным профилем скорости . 45

1.5.2 Расчёт задач с разрывными начальными данными..........51

1.5.2.1 Задача об ударной трубе ............................52

1.5.2.2 Задача о распаде разрыва с сильным перепадом давления ..............................................54

1.5.3 Сравнение квазиакустической схемы со схемами годунов-ского типа второго порядка точности ........................57

1.5.4 Влияние количества горизонтальных слоёв разбиения на качество численного решения..................................60

1.5.5 Зависимость качества численного решения от числа Куранта 61

2 Обобщение алгоритма на случай двух пространственных измерений 64

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

2.2 Линеаризация двумерных уравнений Эйлера газовой динамики . 65

2.3 Построение разностной схемы..........................................69

2.4 Аппроксимация интегральных потоков ..............................75

2.4.1 Линейная реконструкция опорных функций системы ... 75

2.4.2 Замена линейной функции суперпозицией вертикальных столбцов..........................................................77

2.4.3 Определение фоновых подобластей и фоновых значений опорных функций системы....................................78

2.4.4 Разбиение вертикальных столбцов на блоки малых возмущений ............................................................81

2.4.5 Аппроксимация интегральных потоков ......................83

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

2.5 Результаты численных расчётов......................................89

2.5.1 Сравнение результатов расчётов с горизонтальным и вертикальным разбиением опорных функций системы..........90

2.5.2 Расчёт цилиндрического взрыва..............................93

3 Применение квазиакустической схемы к решению прикладной задачи 97

3.1 Обобщение алгоритма на случай трёх пространственных измерений 97

3.1.1 Построение разностной схемы..................................97

3.1.2 Аппроксимация интегральных потоков ........... 99

3.1.2.1 Линейная реконструкция опорных функций системы ......................... 99

3.1.2.2 Замена линейной функции суперпозицией „вертикальных столбцов"..................100

3.1.2.3 Определение фоновых подобластей и фоновых значений опорных функций.............101

3.1.2.4 Разбиение „вертикальных столбцов" на блоки малых возмущений...................101

3.1.2.5 Аппроксимация интегральных потоков......103

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

3.1.3 Результаты расчётов......................104

3.2 Применение квазиакустической схемы к решению прикладной задачи ...................................106

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

3.2.2 Математическое моделирование взаимодействия газовых

струй с экраном-отбойником .................107

3.2.2.1 Картина течения во внутренней и внешней областях экрана-отбойника................110

3.2.2.2 Давление газовых струй на стенки экрана-отбойника .......................111

3.2.2.3 Параллельная версия алгоритма..........114

Заключение 116

Литература 117

Введение

Актуальность темы

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

На сегодняшний день существует большое количество численных методов решения уравнений газовой динамики. Их число продолжает расти, пополняясь всё новыми разработками. С одной стороны, это свидетельствует о важности численного решения уравнений газовой динамики в различных приложениях, а с другой, такое обилие методов численного решения говорит о том, что пока не существует универсального метода, удовлетворяющего всем предъявляемым к нему требованиям. Обзор наиболее распространённых методов численного решения уравнений газовой динамики представлен в работах [45],[41],[39].

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

беппостей которых является возможность возникновения и распространения разрывных решений. При этом разрывы могут появляться даже в случае гладких начальных данных [45].

Во времена, когда производительность ЭВМ была сравнительно небольшой, для расчёта разрывных решений газовой динамики применялись методы, основанные на явном выделении разрывов [30, 1, 17, 41]. Явное выделение разрывов обеспечивало возможность использования небольшого числа сеточных узлов для получения приемлемой точности численного решения и позволяло даже на слабой вычислительной базе решать достаточно сложные классы задач. Методы данного направления показали свою эффективность при расчёте непрерывных, достаточно гладких решений. Однако методы с выделением разрывов зачастую требовали знания априорной информации об особенностях решения, что приводило к усложнению расчётного алгоритма.

С ростом вычислительных мощностей широкое распространение получили методы сквозного счёта (однородные разностные схемы) [51, 47], которые при расчёте разрывных решений не требуют явного выделения разрывов. Исторически схемы сквозного счёта первого порядка точности сильно „размазывали" разрывы (например, схема Лакса-Фридрихса [74, 66],КИР [64]). В тоже время такие схемы позволяли сохранить монотонность профиля решения. Среди схем сквозного счёта первого порядка точности широкое распространение получил метод, предложенный С.К. Годуновым [15, 16]. В классическом методе Годунова для получения решения па новом временном слое применяется кусочно-постоянная реконструкция решения в пределах ячейки расчётной сетки на текущем временном слое и используется точное решение задачи о распаде разрыва (задачи Римана).

Для повышения точности численного решения стали применять схемы сквозного счёта, обладающие вторым порядком аппроксимации в областях гладкости решения (например, схема Лакса-Веидроффа [76], МакКормака [79]). Такие схемы позволили уменьшить область „размазывания" разрывов. Однако

при проведении расчётов в областях, характеризующихся сильным перепадом значений сеточных величин, схемы приводили к появлению осцилляций [48, 44], не имеющих физического характера, а являющихся следствием разностной схемы. Это часто приводило к невозможности проведения адекватного вычислительного эксперимента. Для борьбы с нефизичсскими осцилляциями в разностные схемы с порядком аппроксимации выше первого стали вводить искусственную вязкость [96, 40, 46, 43], которая представляет собой дополнительные члены в разностных уравнениях, имеющие вид формальных аппроксимаций вязких и других диссипативных процессов с коэффициентами, стремящимися к нулю при измельчении сетки. Также применялись линейные и нелинейные фильтры высокочастотных осцилляций [71, 35].

Таким образом, возникла проблема, как повысить порядок точности схемы и при этом обеспечить монотонность численного решения при наличии слабых и сильных разрывов. В работах С.К. Годунова [15, 16[ показано, что не существует линейных монотонных разностных схем выше первого порядка точности (теорема С.К. Годунова). Выходом из противоречия меж/1у получением монотонного решения и повышением порядка аппроксимации стало появление нелинейных разностных схем. К таким нелинейным монотонным схемам относятся гибридные схемы, которые могут локально менять свои свойства (например, порядок аппроксимации). В частности, гибридные схемы позволяют применять схемы с повышенным порядком аппроксимации в областях гладкости решения, а в узлах, в которых решение имеет разрывы, использовать монотонные разностные схемы первого порядка точности. В работе Р.П. Федоренко [58] предложена первая гибридная схема для линейного уравнения переноса, в которой переключение между базовыми схемами первого и второго порядка точности происходит на основе анализа отношения второй конечной разности решения к её первой разиости. Гольдин, Калиткин, Шишова [25] построили несколько гибридных схем для линейного и квазилинейного уравнений переноса, применив гладкое переключение между схемами первого и второго порядков. В их

подходе коэффициент переключения между базовыми схемаим зависел от градиента решения. Harten и Zwas [71, 72] предложили первые гибридные схемы для гиперболической системы уравнений общего вида. Так в работе [72] использовалась комбинация схемы Лакса-Фридрихса [74] первого порядка точности и схемы Лакса-Вендроффа второго порядка точности [75, 76]. Краткий обзор подходов к созданию первых гибридных разностных схем может быть найден в работах [36, 37].

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

Дальнейшее развитие идей, заложенных в работах Бориса-Бука, привело к созданию TVD (Total Variation Diminishing) схем. В работах [67, 68] Харте-ном для численного решения квазилинейного уравнения переноса предложен класс нелинейных разностных схем, получены условия, которыми должны обладать коэффициенты, чтобы схема удовлетворяла условию TVD, которое со-

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

Для того, чтобы численное решение удовлетворяло условию ТУБ, в настоящее время развита и широко применяется процедура кусочно-линейной (кусочно-полиномиальной) реконструкции сеточных функций. При этом наклоны кусочно-линейных распределений сеточных функций для удовлетворения условию ТУБ в пределах ячеек расчётной сетки ограничиваются специальными ограничителями, которые получили название лимитеры (Нт^егв). Среди первых работ в данном направлении следует отметить работы В.П. Колгана [34, 35]. В схеме Колгана [34] для повышения порядка точности классического метода С.К. Годунова [15] решению задачи о распаде разрыва предшествует этап кусочно-линейной реконструкции сеточных величин. Для подготовки начальных данных для задачи о распаде разрыва применяется процедура интерполяции (или экстраполяции) сеточных величии на текущем временном слое в граничные точки ячеек расчётной сетки, которая основана на принципе минимальных значений производных [34]. Применение принципа минимальных значений производных позволяет обеспечить монотонность решения. Построенная в работе [34] разностная схема аппроксимирует нестационарные решения со вторым порядком точности по координатному направлению и с первым порядком точности по временному направлению. Схема Колгана сохраняет монотонность решения и способствует уменьшению размазывания контактных разрывов и слабых скачков уплотнения, а также позволяет достичь большей точности в

областях гладкости решения.

Позднее Брам Ван Лир опубликовал серию работ [93, 94, 95[, в которых предложил получившую в дальнейшем широкое распространение схему MUSCL (Monotone Upstream-centered Schemes for Conservation Laws). Схема MUSCL, также как и схема Колгана [34], использует кусочно-линейную реконструкцию сеточных величин и обладает вторым порядком точности. Однако, в отличие от схемы Колгана, построение схемы MUSCL [95] для уравнений газовой динамики проводится в лагранжевых переменных, а отдельным этапом реализуется пересчёт полученного решения на эйлерову сетку. Для обеспечения монотонности решения применяются более сложные, чем в схеме Колгана, алгоритмы монотонизации [94], которые основаны на применении лимитерных функций, ограничивающих наклон линейных восполнений сеточных величин для борьбы с нефизическими осцилляциями. На сегодняшний день известно большое количество лимитерных функций, которые обладают различными свойствами и выбираются в зависимости от специфики конкретно решаемой задачи. Обзор наиболее распространённых лимитеров можно найти в работах [89, 97, 73, 90].

Использование кусочно-параболической реконструкции сеточных величин позволяет повысить порядок точности схемы до третьего по пространственной переменной. Среди методов данного направления широкое распространение получили PPM (Piecewise Parabolic Method) метод [63], метод Чакраварти-Ошера [62].

Среди методов, использующих реконструкцию сеточных величин, можно также отметить ENO метод (Essentially Non-Oscillatory), основанный на автоматическом анализе гладкости решения [70]. Процедура автоматического анализа гладкости решения состоит в определении ячеек сетки, в которых находятся разрывы сеточной функции, разрывы её первой производной (слабые разрывы) и разрывы производных более высокого порядка в соответствии с требуемым порядком аппроксимации. В зависимости от локального уровня гладкости, для реконструкции сеточных величин применяются интерполяционные многочле-

ны Ыыотона соответствующего порядка, степень которого определяет порядок аппроксимации. Это приводит к тому, что в методе ENO для интерполяции значений сеточных величин в ячейках расчётной сетки используется несколько шаблонов и среди них выбирается тот, на котором решение является наиболее гладким. Дальнейшее развитие метод ENO получил в работах [87, 88].

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

Интенсивное развитие TVD-алгоритмов привело к значительному повышению качества получаемых численных ре