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

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

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

ГАЛАНИНА Анна Михайловна

МЕТОД ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ

ГАЗОДИНАМИЧЕСКИХ ТЕЧЕНИЙ И ЕГО ПРИМЕНЕНИЕ В ЗАДАЧЕ О Т-СЛОЕ

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

численные методы и комплексы программ

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

14 НОЯ 2013

005538559

Москва — 2013

005538559

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

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

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

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

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

фессор Фаворский Антон Павлович |

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

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

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

Автореферат разослан 1 ноября 2013 г.

Ученый секретарь диссертационного совета Д 501.001.43 доктор физико-математических наук, профессор

Е.В. Захаров

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

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

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

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

ложили схему для расчета течений с разрывами.

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

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

Проблемой решения гиперболических уравнений, в частности, уравнений газовой динамики, занимались такие учёные, как Courant R., Einfeldt В., Friedrichs К.О., Harten A., Isaacson Е., Lax P.D., Osher S., Rees M., Roe P.L., von Neumann J., Wendroff В., Белоцерковский О.M., Годунов C.K., Головизнин В.M., Кузнецов O.A., Попов Ю.П., Рождественский Б.Л., Самарский A.A., Тишкин В.Ф., Фаворский А.П., Яненко H.H. Их работы послужили важным и ценным источником для выполнения данной диссертации.

Правильное решение задач гиперболического типа позволяет получить практически важные результаты в прикладных задачах. К таким задачам относится, в частности, моделирование эффектов Т-слоя. Эффект Т-слоя обнаружен в 60-х годах прошлого века в результате чисто теоретических исследований и вычислительного эксперимента группой ученой во главе с А.Н. Тихоновым (Самарский A.A., Заклязьминский JI.A., Волосевич П.П., Дегтярев JIM., Курдюмов С.П., Попов Ю.П., Соколов B.C., Фаворский А.П.). Позднее появились факты, экспериментально подтверждающие существование эффекта Т-

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

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

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

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

2. построение и исследование консервативной однородной монотонной разностной схемы для решения квазилинейного уравнения переноса;

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

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

5. построение численного метода решения двумерной системы МГД уравне-

ний;

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

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

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

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

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

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

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

3. МГД математическая модель образования и развития Т-слоя в потоке сла-бопроводящего газа в присутствии магнитного поля; условия образования и развития Т-слоя в зависимости от параметров магнитного поля, скорости течения и геометрических характеристик возмущения.

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

в переменных Лагранжа и в переменных Эйлера, в двумерном случае в эйлеровых координатах, а также уравнений магнитной гидродинамики в двумерном случае.

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

Проведено сравнение численных результатов, полученных с помощью предложенного алгоритма и с помощью известных методов решения соответствующих задач. Так, для линейного уравнения переноса проведено сравнение результатов численного решения с помощью сконструированной сплайн-схемы с численным решением, полученным по схеме Лакса-Вендроффа. Для квазилинейного уравнения переноса сравнение проведено с монотонизированной схемой К.И. Вабенко (схемой „квадрат"). Квазиакустическая схема для решения уравнений газовой динамики сравнивалась со схемой Роу-Эйнфельдта-Ошера повышенного порядка аппроксимации.

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

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

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

1. Международной конференции „Современные проблемы вычислительной математики и математической физики", секция „Математическое моделирование" (г. Москва, 2009 г.);

2. Тихоновских чтениях в Московском государственном университете им. М.В. Ломоносова, секция вычислительной математики и кибернетики (г. Москва, 2009 г.);

3. Международной молодежной конференции-школе „Современные проблемы прикладной математике и информатике" (г. Дубна, 2012 г.);

4. Ломоносовских чтениях в Московском государственном университете им. М.В. Ломоносова, секция вычислительной математики и кибернетики (г. Москва, 2013 г.);

5. Одиннадцатом международном семинаре „Mathematical Models к Modeling in Laser-Plasma Processes к Advanced Science Technologies" (г. Будва, Черногория, 2013 г.).

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

Публикации. Положения диссертации отражены в 11 печатных работах: 6 статьях [1-6] (5 из которых опубликованы в изданиях, включенных в Перечень ведущих научных журналов и изданий для опубликования основных научных результатов диссертации, рекомендованных ВАК РФ [1-5]), 5 тезисах докладов [7-11].

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

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

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

Первая глава посвящена построению и исследованию сплайн-схемы для линейного и квазилинейного уравнения переноса. Для линейного уравнения переноса поставлена задача Коши

!«!=„, -о°<*<+~, «»о. (1)

/(я,0) = /о(®), а = const > 0.

Решение построено численно на равномерной сетке ш =

7

{(хк, Ъ), к = 1,..., ЛГ, j = 0,..., Г} с шагами Лиг.

Построение схемы проведено интегро-интерполяционным методом. В результате интегрирования уравнения переноса по пространственно-временно ячейке \xk-xi2,2^+1/2] х ^+1] имеем интегральное тождество:

К?к+1-Гк) + Ж1к+у-Ш1к_и = о, . (2)

где

я = 1 I /(*• = I а/(хкЦ, №. (3)

Соотношение (2) описывает семейство разностных схем. Задача поиска численного решения тем самым сведена к аппроксимации интегральных потоков И/7. Геометрический смысл интегрального потока Ш1к+х/2^ - площадь криволинейной трапеции с основанием [2^+1/2 ~ ат> хк+1/2]-

Для вычисления интегральных потоков \У1 проинтерполируем сеточную функцию у1 — у^{хк) на текущем временном слое, получив функцию непрерывного аргумента ${х). Для интерполяции воспользуемся простейшими линейными сплайнами, которые на каждом отрезке [2^-1/21 ®*+1/2] в момент времени ^ имеют вид:

&(х)=й + 0{(х-хк),

п ЫУх + \Ух\Ух ^

к Ш + Ы '

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

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

нения переноса:

2/Г = 2/£ - 7(24 " - " ~ 7 = Ъ (5)

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

Следующий подраздел первой главы посвящен численному решению задачи Коши для квазилинейного уравнения переноса:

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

Слайн-схема для решения квазилинейного уравнения переноса (в том числе со знакопеременными начальными данными) имеет вид:

^-+/^- = 0, -ос Сс<+оо, £ > О, ся ах

/(*,()) = /„(х).

д}

дх

(6)

_ р _

1к ~ ■> к

к-1

н

уЬ{ = = у{ - уЯ{ = У*Ыь) = у1 +

,0, уЯ1< о, уЯ{ = У:(чц) = У>к + -21-

Рис. 1. Геометрическая интерпретация интегрального потока для квазилинейного

уравнения переноса.

В работе приведены результаты тестовых расчетов и проведено сравнение решений, полученных с помощью сплайн-схемы и с помощью монотонизиро-ванной схемы К.И. Бабенко („квадрат"). Показано, что обе схемы дают практически неотличимые результаты при значительных различиях в сложности алгоритма.

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

д_ = ди

дь \ Р ) а?

д(Р°(Ои) дР_

Ы ' (8)

к

дх _

здесь ро(£) ~ значение плотности в начальный момент времени. Численное реше-

ние задачи строится на равномерной сетке с узлами где & = (г+кИ, к =

0,1,..., ТУ, ^ = з = 0,1,..., J, N - фиксированное число ячеек пространственной сетки с шагом к — «2 - СО/Л^, т - постоянный шаг по времени.

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

< Нч?1 ~ <&) + МСЛЦ - = 0, , (9)

-е1) + ШЕ1Ц-ШЕ11г О,

где Ах{ - объём одномерной частицы, - объёмный импульс, а - средняя по объёму величина количества энергии в ячейке разностной сетки. Слагаемые 1/2> ^ЯРк+1/2 и ШЕРк+1/2 представляют собой интегральные потоки соответствующих величин через границы ячейки расчётной сетки за время т. Задача построения разностной схемы тем самым состоит в аппроксимации потоков.

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

+ = « + = + (Ю)

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

скоростями:

' ™и1{+1/2 = йт + Ш{+у2 + Ш{+1/2, ' 1/2 = йрт + (рШ1+1/2 + йШ&{+1/2) +

+ (рШ{+1/2+й^:+1/2).

Для разложения функций на суммы вида (11) предложен следующий алгоритм. В каждой ячейке [Сь—г /г, Сь-ы /2] используем локальную сплайн-реконструкцию опорных функций р, и и р в виде (4). На отрезке про-

водим расслоение (разбиение) линейных сплайн-функций (под f пони-

мается одна из функции р, и или р) и /¿+1(£, tj) на совокупность горизонтальных слоёв (с номерами т = 1,..., М), параллельных оси абсцисс (рис. 2 (а)). Количество слоёв М выбираем таким образом, чтобы слои были достаточно тонкими и их можно было считать малыми возмущениями. Горизонтальное сечение, ближайшее к среднему между значениями сплайн-функций 1/2, и /1+1^+1/2, на границе раздела £ = двух соседних ячеек (рис. 2 (а)), принимаем общим постоянным фоном для функции / между этими ячейками.

Каждую из сплайн-функций и /¿+1(£, tj) на полуячейках ^+1/2],

[&+1/2! заменяем изображённой на рис. 2 (а) конструкцией, состоящей из общего постоянного фона />,*;+1/2, на котором последовательно располагаются „кирпичи" малых ступенчатых возмущений (с учётом знака). Существенно при этом, что каждый „кирпич" имеет своим индивидуальным постоянным фоном либо общий постоянный фон либо другой „кирпич", примыкающий к

нему со стороны общего постоянного фона.

Таким образом, в соответствии с геометрической интерпретацией интегральных потоков их вычисление сводится к суммированию площадей тех частей „кирпичей", которые при движении в сторону границы £ = &+1/2 со своей скоростью звука пересекли её за время г (рис. 2 (б)).

Дк+1/2

к+1

Рис. 2. (а) Расслоение сплайн-функции на „кирпичи" малых возмущений, определение фонового значения; (б) движение „кирпичей" и вычисление интегральных потоков.

Из описанного алгоритма непосредственно вытекает условие устойчивости данного метода: величина шага по времени т должна быть такой, чтобы ни один „кирпич" не вышел за пределы рассматриваемого отрезка:

срт . Ъ,

< то Ф = тахскРк'

р (О 2 кл

(12)

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

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

0.0401 11Р"Р«1к

0.00 0.02 0.04 0.06 0.08 0.10 К х=0.1И

Рис. 3. Решение задачи с гладким начальным профилем. График функции плотности на последовательные моменты времени (а) и график зависимости погрешности решения от шагов сетки при постоянном числе Кураната (6).

полученного решения с учётом конвективных потоков. Алгоритм второго этапа основан на методе крупных частиц О.М. Белоцерковского. Основной величиной, используемой на этом этапе, является поток массы через границу ячейки, то есть вычисленный на первом этапе. Такая модификация оригиналь-

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

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

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

1.0

0.5

0.0 0.2 0.4 0.6

(а)

1.0

0.8 ' 1.0 X

0.5

0.0 0.2 0.4

~06~

(б)

0.8

1.0

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

ограничения на шаги сетки.

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

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

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

Электропроводность газа, зависящая от его температуры, недостаточна для обеспечения взаимодействия газа с магнитным полем. В начальный момент времени магнитное поле однородно и имеет только одну отличную от нуля компоненту, направленную вдоль оси г. Невозмущённая начальная температура равна То. В начальный момент времени £ = 0 в поток газа вносится локальное возмущение температуры до величины Тв > Т0. При этом электропроводность в возмущённой зоне возрастает и становится достаточной для взаимодействия газа с магнитным полем.

Для поставленной задачи использован метод расщепления по физическим процессам, в результате чего численное решение системы распадается на три этапа:

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

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

3. пересчет газодинамических параметров с учетом электромагнитного взаимодействия.

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

На рис. 5 (а) представлены результаты численного решения поставленной задачи. После небольшого падения температуры в первые моменты времени

4)8 л Т401, К

4,6 4'4 1

4,2 - /' \! '

4,0 3,8 3,6" 3,4 " 3,2 -3,0

' \

т

л- ✓

0

I

10

15

20

1,0 — 0,80,6" 0,40,2-

о,о —

с о-

о

Л

г>

Л

00 л

10

(а)

(б)

Рис. 5. (а) Профиль температуры вдоль оси канала в различные моменты времени (¿о = 0, = 2 мсек, ¿2 = 4 мсек, 43 = 6 мсек, ¿4=8 мсек); (б) линии уровня температуры для возмущения на различные моменты времени.

начинается ее интенсивный рост, и вскоре максимальное значение превышает начальную температуру возмущения. На рис. 5 (б) представлены линии уровня температуры. Первоначально двумерный профиль с течением времени эволюционирует в одномерную структуру.

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

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

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

Основные результаты диссертационной работы

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

2. метод обобщен для решения уравнений газовой динамики и уравнений магнитной гидродинамики в двумерной постановке на регулярной равномерной прямоугольной сетке;

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

4. построенные алгоритмы использованы для решения МГД математической модели образования и развития Т-слоя в потоке слабопроводящего газа в присутствии магнитного поля; в вычислительных экспериментах исследованы условия образования и развития Т-слоя в зависимости от начальных параметров магнитного поля, скорости течения и геометрических характеристик возмущения.

Основные результаты диссертации изложены в работах

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

2. Фаворский А.П., Тыглиян М.А., Тюрина H.H., Галанина A.M., Исаков В.А. Численное моделирование распространения акусти-

ческих импульсов в гемодинамике // Дифф.уравнения. 2009. Т. 45, №8. С. 1179-1187.

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

4. Галанина A.M., Фаворский А.Н. Численное решение уравнений газовой динамики в лагранжевых переменных // Матем. моделирование. 2012. Т. 24, №12. С. 119-123.

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

6. Galanina A., Favorskiy A. Local two-dimensional perturbations evolution in smail-conductivity gas flow in magnetic field // Mathematica Montisnigri. 2013. V. XXVII. P. 53-64.

7. Фаворский А.П., Тыглиян M.А., Тюрина H.H., Галанина A.M., Исаков В.А. Численное моделирование распространения гемодинамиче-ских импульсов // Современные проблемы вычислительной математики и математической физики: тез. докл. международной конференции памяти академика A.A. Самарского. М.:, 2009. С. 376.

8. Галанина A.M. Численное решение уравнений газовой динамики в лагранжевых переменных // Сборник тезисов лучших дипломных работ 2010 года. МГУ им. М.В. Ломоносова. Ф-т ВМК, 2010. С. 32-34.

9. Фаворский А.П., Галанина A.M., Исаков В.А. Конструктивный подход к численному решению квазилинейных уравнений переноса //

Современные проблемы прикладной математики и информатики: тез. до-кл.международной молодежной конф.-школы, 2012. С. 34.

10. Фаворский А.П., Галанина A.M. Развитие двумерных локальных возмущений в потоке слабопроводящего газа в магнитном поле // Ломоносовские чтения 2013: тез. докл. конференции, 2013. С. 88-90.

11. Galanina A., Favorskiy A. Local two-dimensional perturbations evolution in small-conductivity gas flow in magnetic field // Mathematical models modeling in laser-plasma processes advanced science technologies: abstract. Eleventh international seminar. Budva, 2013. P. 43.

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

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

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

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

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

04201365262

ГАЛАНИНА Анна Михайловна

МЕТОД ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ

ГАЗОДИНАМИЧЕСКИХ ТЕЧЕНИЙ И ЕГО ПРИМЕНЕНИЕ В ЗАДАЧЕ О Т-СЛОЕ

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

численные методы и комплексы программ

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

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

д. ф. - м. н., проф. А. П. Фаворский

д ф. - м. н. С. И. Мухин

Москва — 2013

Оглавление

Введение 4

1 Сплайн-схема для решения линейного и квазилинейного уравнений переноса 15

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

1.2 Построение и исследование сплайн-схемы для линейного уравнения переноса............................................................16

1.3 Результаты расчётов. Перенос различных профилей................25

1.4 Сплайн-схема для квазилинейного уравнения переноса............30

1.4.1 Построение точного решения. Геометрическая интерпретация ................................................................30

1.4.2 Построение разностной сплайн-схемы. Учет нелинейности

при аппроксимации потоков....................................33

1.5 Примеры численного решения. Сравнение с монотонизированной схемой К.И. Бабенко....................................................37

2 Квазиакустическая схема 43

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

2.2 Квазиакустическая схема для одномерной системы уравнений газовой динамики в переменных Лагранжа............................47

2.2.1 Уравнение неразрывности......................................47

2.2.2 Уравнение движения............................................49

2.2.3 Уравнение энергии..............................................50

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

2.3.1 Линеаризованная система......................................52

2.3.2 Аппроксимация потоков........................................54

2.3.3 Интерполяция сеточных функций. Построение фоновых значений..........................................................57

2.3.4 Алгоритм вычисления интегральных потоков................61

2.4 Результаты расчётов для одномерной системы в переменных Лагранжа................................................................63

2.5 Переход от переменных Лагранжа к переменным Эйлера..........72

2.6 Результаты расчётов для одномерной системы в переменных Эйлера. Сравнение со схемой Роу-Эйнфельдта-Ошера................75

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

2.8 Результаты численного решения двумерной системы уравнений газовой динамики ......................................................84

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

газа в магнитном поле 86

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

3.2 Алгоритм численного решения........................................91

3.2.1 Уравнение диффузии магнитного поля ......................93

3.2.2 Джоулев нагрев и сила Лоренца. Пересчет решения с учетом электромагнитного взаимодействия......................98

3.3 Результаты вычислительных экспериментов ........................99

Введение

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

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

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

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

сколь угодно гладких начальных данных решение может оставаться гладким и, соответственно, классическим лишь в течение ограниченного отрезка времени. Построение методов расчета разрывных решений представляет собой особую проблему при решении гиперболических уравнений [37]. Важнейший шаг в этом направлении сделали фон Нейман и Рихтмайер, которые впервые в 1950 г. предложили схему для расчета течений с разрывами [92].

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

Достаточно полный обзор методов решения гиперболических уравнений можно найти в [36]. Основные методы на настоящий момент времени можно условно разделить на несколько типов, при этом большинство методов существенно использует характеристические свойства гиперболических систем [16, 36, 40] (достаточно, например, отметить метод характеристик [49]). Исторически одними из первых схем, используемыми для решения, в частности, уравнений газовой динамики, стали однородные разностные схемы [51, 85, 86]. Использование однородных разностных схем удобно, т.к. в этом случае не требуются специальные методы для выявления разрывов и специальные алгоритмы для их расчетов. Тем не менее существует ряд существенных ограничений, накладываемых на подобные схемы. Во-первых, как показано в [60], разностные схемы для гиперболических уравнений должны обладать свойством консервативности, т.к. нарушение этого принципа при расчете разрывов может приводить к серьезным ошибкам, лишающим разностную схему важнейшего свойства сходимости. Во-вторых, в соответствии с теоремой С.К. Годунова [15] из трех свойств разностной схемы для решения линейной системы гипербо-

лических уравнений - линейности, аппроксимации с порядком выше первого и монотонности - одновременно могут иметь место только два. В результате при использовании схем повышенного порядка аппроксимации данного типа (например, схемы Лакса-Фридрихса [79, 85], схемы Лакса-Вендроффа [86], схемы Куранта-Изаксона-Риса [77]) в окрестностях разрыва возникают сильные осцилляции, что связано с недостаточной диссипацией подобных разностных схем.

Для решения этой проблемы предложено два принципиальных решения. Один из наиболее распространенных вариантов - введение так называемой искусственной вязкости [52, 70]. В этом случае в уравнение вводятся дополнительные слагаемые, формально аппроксимирующие диссипативные процессы и стремящиеся к нулю при измельчении шагов разностной сетки. Дополнительная вязкость обеспечивает монотонность решения в окрестности скачков, однако это приводит к сглаживанию разрывной функции. Другой способ преодоления противоречия, вытекающего из теоремы С.К. Годунова, - это схемы линейной гибридизации, представляющие собой схемы повышенного порядка аппроксимации в областях гладкости решения, и схемы первого порядка аппроксимации с достаточной вязкостью в окрестностях разрыва. При этом разностный алгоритм часто перестает быть однородным. В дальнейшем появились схемы с ограничителями антидиффузионных потоков [8. 9, 32, 91], среди которых важно выделить ТУБ-схемы [74, 82], ЕШ-схемы [83], \¥ЕШ-схемы [87]. Основными недостатками данного класса схем являются неоднородность алгоритма, необходимость использования анализаторов гладкости, а также введение дополнительных слагаемых, не существующих в исходной дифференциальной физической модели. Кроме того, ЕМО- и \УЕГ,'Ю-схемы достаточно трудоемки.

В настоящий момент времени все большую популярность получают компактные [42, 90] и бикомпактные схемы [45, 46, 47, 48]. Эти схемы характеризуются абсолютной устойчивостью, монотонностью, высоким порядком аппроксимации и экономичностью, т.к. имеют компактный шаблон и могут быть

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

Важным шагом в развитии методов решения гиперболических уравнений стал переход от схем, предполагающих гладкость решения, к схемам, базирующимся на взаимодействии разрывов. Так как системы гиперболического типа допускают существование разрывных решений [38], одним из важнейших критериев качества предлагаемого численного алгоритма является аккуратный расчет разрывов. Первоначально С.К. Годунов предложил для вычисления потоков через границу расчетной сетки решать задачу о распаде произвольного разрыва (задачу Римана) [15, 17, 18]. В этом случае на каждой границе ячейки решается система нелинейных алгебраических уравнений, что требует значительных временных и вычислительных затрат. Поэтому дальнейшее развитие этой идеи привело к появлению методов годуновского типа с приближенным решением задачи о распаде разрыва [78, 84, 88, 89].

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

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

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

нать с решения уравнения переноса - одного из фундаментальных уравнений математической физики [61]. Важнейшим свойством точного решения данного уравнения является перенос начального профиля без искажения [61]. Однако устойчивые схемы первого порядка аппроксимации, как правило, дают „расплываюдщееся" со временем по пространству решение [57], а схемы повешенного порядка аапроксимации в соответствии с теоремой С.К. Годунова [31] зачастую сильно искажают решение вследствие нормальной или аномальной дисперсии [20, 30, 44, 51]. Существуют различные способы повышения качества получаемого решения. Один из наиболее распространенных основан на использовании однородных схем повышенного порядка аппроксимации с коррекцией решения путем введения „лимитеров" таким образом, чтобы разностная схема удовлетворяла принципу максимума [54]. Число вариантов монотонизации разностных схем очень велико: одни из первых представлены в работах [5, 25, 53, 70], более новые, основанные на знании области зависимости точного решения, представлены в [20, 21]. Методом иного класса для решения уравнения переноса (в том числе и нелинейного) является разрывный метод Галер-кина (ШШС-метод) [75]. Этот метод дает более высокую точность решения, а также меньшее „размазывание" разрывов по сравнению со многими конечно-разностными схемами [3]. Однако этот метод оказывается более трудоемким в реализации.

Решение системы уравнение газовой динамики представляет особый интерес. Традиционно для описания движения сплошной среды используются два подхода - Эйлера и Лагранжа, связанных с выбором системы координат [38]. Интерес представляют собой лагранжевы массовые координаты [57] и ряд разностных схем [57], применимых в этом случае. Однако лагранжевы массовые координаты успешно удается ввести только в одномерном случае, поэтому такой подход не является универсальным. Традиционные лагранжевы и эйлеровы координаты играют важную роль в современных прикладных задачах. Лагранжевы координаты чаще используются в схемах, основанных на вариа-

ционном подходе [24]. Этот подход универсален, допускает использование любых пространственных координат (декартовых, цилиндрических, сферических и т.д.) и любое количество пространственных измерений [22]. Другим серьезным преимуществом лагранжевых координат является тот факт, что контактные разрывы в этом случае выделяются фактически в явной форме и остаются неподвижными, что значительно облегчает, например, выбор коэффициентов искусственной вязкости [19]. Идея сочетания лагранжева и эйлерова подходов лежит в основе Р1С-метода [43] и метода ,.крупных частиц" [4]. На настоящий момент построение качественных схем для уравнений газовой динамики в обеих системах координат, а также в смешанных эйлерово-лагранжевых координатах [23, 71] является актуальной и важной задачей.

Настоящая работа посвящена построению и исследованию разностной схемы для решения гиперболических уравнений с основным упором на решение системы уравнений газовой динамики. В качестве прикладной задачи, для решения которой применяется предложенный алгоритм, выбрана магнитогидро-динамическая задача о развитии и устойчивости Т-слоя [63]. Эффект Т-слоя обнаружен в 60-х годах прошлого века в результате чисто теоретических исследований и вычислительного эксперимента [63]. Позднее появились факты, экспериментально подтверждающие существование эффекта Т-слоя [27]. Т-слой представляет собой высокотемпературное самоподдерживающееся образование, возникающее и развивающееся в плазме при определенных условиях в процессе ее взаимодействия с магнитным полем [62]. Физическая сущность явления состоит в том, что при нестационарном взаимодействии плазмы и магнитного поля при определенных условиях в плазме образуется самоподдерживающийся высокотемпературный электропроводный слой (Т-слой), усиливающий взаимодействие плазмы с магнитным полем [6, 55].

После открытия явления Т-слоя появилось большое количество работ, посвященных исследованию данного процесса. Пространственно-одномерная модель магнитогидродинамического течения рассмотрена в [62], двумерная мо-

дель изучена методом вычислительного эксперимента позднее в работе [6|. Теоретическое исследование магнитогидродинамической структуры, содержащей Т-слой, проведено в работе [59]. Рассмотрению дополнительных факторов, влияющих на возникновение и развитие Т-слоя (таких, как излучение и теплопроводность), посвящена работа [7]. Суммарное достаточно полное исследование, результаты вычислительных и практических экспериментов, возможности практического применения явления Т-слоя описаны в [26].

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

Цель и задачи исследования

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

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

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

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

2. построение и исследование консервативной однородной монотонной раз-

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

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

4. построение и исследован