ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 8, с. 1363-1379
УДК 519.634
Посвящается светлой памяти А.П. Фаворского
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТРЕХМЕРНЫХ ТЕЧЕНИЙ КВАЗИНЕЙТРАЛЬНОГО ГАЗА НА ОСНОВЕ СГЛАЖЕННЫХ УРАВНЕНИЙ МАГНИТНОЙ ГИДРОДИНАМИКИ^
© 2015 г. Т. Г. Елизарова*, М. В. Попов**
(*125047Москва, Миусская пл., 4, ИПМРАН;
**École Normale Supérieure de Lyon, CRAL, UMR CNRS 5574, Université de Lyon 1, 46 allée d'Italie 69007Lyon, France) e-mail: telizar@mail.ru; mikhail.v.popov@gmail.com Поступила в редакцию 26.01.2015 г.
Представлен новый конечно-разностный метод для численного моделирования сжимаемых МГД-течений, применимый к широкому классу задач. Метод состоит в использовании магнитных квазигазодинамических уравнений (КМГД-уравнений), которые, по сути, являются системой уравнений Навье—Стокса и уравнений Фарадея, к которым была применена процедура усреднения на малом интервале по времени. КМГД-уравнения дискретизируются на расчетной сетке с помощью центральных разностей. Усреднение позволяет стабилизировать численное решение и не использовать дополнительные ограничивающие процедуры (лимитеры и пр.). Бездивергентность магнитного поля обеспечивается применением теоремы Стокса. Представлены результаты расчетов тестовых 3Б-задач: центральный взрыв в магнитном поле, взаимодействие ударной волны с облаком и трехмерный тест Орсзага—Танга. Также продемонстрированы предварительные расчеты плазменного пинча, удерживаемого магнитным полем в ловушке. Библ. 16. Фиг. 1З.
Ключевые слова: магнитная квазигазодинамика, КМГД, МГД-течения, конечно-разностный алгоритм, центрально-разностная аппроксимация.
DOI: 10.7868/S0044466915080086
1. ВВЕДЕНИЕ
В статье представлен численный алгоритм и результаты моделирования нестационарных пространственных течений идеальной квазинейтральной плазмы, находящейся в поле электромагнитных сил. Численный алгоритм является расширением построенных ранее конечно-разностных схем, основанных на КГД-уравнениях для описания течений вязкого сжимаемого газа (см. [1]—[3]). Работа тесно связана с кругом интересов Антона Павловича Фаворского, который, являясь признанным специалистом в области вычислительной гидродинамики, в свое время активно поддержал еще молодое в то время научное направление, развитие которого привело к формированию КГД-подхода. Антон Павлович Фаворский как руководитель научной школы в области вычислительной магнитной гидродинамики положительно оценил бы полученные здесь результаты.
Квазигазодинамические уравнения (КГД-уравнения) выражают законы сохранения для усредненных по малому временному интервалу газодинамических переменных — плотности, компонентов импульса и энергии. При этом делается предположение, что усредненная величина является гладкой функцией времени, которую можно разложить в ряд Тейлора в окрестности каждого момента времени t, т.е.
t+Т
f(r,t) = 1 f f(r, f)df « f(r,t) + + ••• . (1)
T J д t
t
1) Работа выполнена при финансовой поддержке РФФИ (код проекта 13-01-00773_а).
1363
6*
Закон сохранения для усредненной величины /(г, г) должен содержать слагаемые, отражающие как сохранение / (г, г), так и сохранение поправки, пропорциональной малому параметру т, имеющему размерность времени. Таким образом, КГД-уравнения представляют собой уравнения Навье—Стокса, содержащие дополнительные диссипативные слагаемые. Эти слагаемые играют стабилизирующую роль при численном решении уравнений, так как являются причиной дополнительной диссипации.
КГД-уравнения, в которых учтено влияние магнитных полей, впервые были рассмотрены в [3], [4] для описания течений вязкого газа и жидкости. При этом влияние поля учитывалось в виде магнитных сил и диссипативных т -поправок в части газодинамики, а само поле описывалось уравнениями Максвелла без коррекции с помощью т -поправок. Для этой системы было построено уравнение баланса энтропии, выписано точное решение задачи Гартмана, проведен расчет течения электропроводного расплава в безындукционном приближении.
Однако описанную выше процедуру усреднения можно применить и к уравнениям для магнитной индукции, записанным в рамках единой системы магнитной гидродинамики. Этот подход позволит описывать магнитные вязкие течения с помощью квазигазодинамических уравнений для магнитной гидродинамики (КМГД-уравнений) самосогласованным образом. Такие уравнения впервые рассмотрены в работах [5], [6], где выполнено их исследование на некоторых стандартных Ш- и 2Э-тестах: задаче Римана, задаче о распространении магнитных волн, задаче о диссипации и распаде альфвеновской волны, задаче о взрывной волне в замагниченной среде, задаче о вихре Орсзага—Танга и задаче о взаимодействии ударной волны с облаком. Во всех случаях была продемонстрирована хорошая сходимость решения к точному при сгущении разностной сетки.
В [7] КМГД-уравнения были обобщены на случай уравнения состояния неидеального газа при наличии внешних сил и источника тепла. Было получено уравнение теплового баланса и исследованы энтропийные свойства КМГД-уравнений.
В данной работе мы приведем систему КМГД-уравнений для 3Э-случая, записанную покомпонентно и представим результаты расчетов тестовых 3Э-задач: центральный взрыв в магнитном поле, взаимодействие ударной волны с облаком и трехмерный тест Орсзага—Танга. Также будут продемонстрированы предварительные расчеты плазменного пинча, удерживаемого продольным магнитным полем. Проблема устойчивости пинча является важнейшей в технологической задаче удержания плазмы в магнитных ловушках.
2. КМГД-УРАВНЕНИЯ
Запишем КМГД-уравнения в декартовой системе координат с использованием стандартных обозначений для независимых переменных: р — плотность, их, иу, и, — компоненты скорости, Вх, Ву, В, — компоненты магнитной индукции, Е — полная энергия единицы объема.
Будем также использовать краткие обозначения для квадратов модулей векторов скорости и магнитной индукции:
и2 = и2 + и2у + иI, В2 = ВX + В2 + В2.
Множитель >/1/4п мы включили в определение вектора магнитной индукции В. В этих обозначениях полная энергия единицы объема записывается в виде
Е = ре + Ри! + в!, 2 2
где е — удельная внутренняя энергия. Для замыкания системы уравнений потребуется уравнение состояния. Для случая идеального газа оно имеет вид р = (у - 1)рб, где р — газодинамическое давление, у — показатель адиабаты. Уравнение состояния, выраженное через температуру, имеет вид р = рЯТ/п, где Я — универсальная газовая постоянная, п — средний молекулярный вес газа. Отсюда температура Т выражается по формуле Т = рц/рЯ.
Комбинация Е, р и р дает полную удельную энтальпию:
Н = (Е + р)/р.
Малый параметр, имеющий размерность времени, по которому проводится усреднение, обозначим через т (см. (1)). Для удобства второе слагаемое в разложении Тейлора в (1) обозначим через приращение А:
7 = 7 + А/.
Выпишем приращения всех величин, которые будут использованы в дальнейшем. Данные выражения получены из МГД-уравнений для невязкой и нетеплопроводной квазинейтральной плазмы (см. [8]):
А1 = —т| ^гаё- —1 Шуи |, Р Р
Аб = -т| иgгad6 + р ё1уи I, Ар = -т (ugгadp + урё1уи), Р
Аи;- = -т
и^аёи;. + - 1( р+в!] -1 {дВЛ+дМ+§ВА
рдIV 2 у ру дх ду д-
АВ( = т -д (Ви - ихВ,) + -д (Вущ - иуВ1) + 5 (Ви - иВ) ,
дх ду д- _
где I = х, у, z, и — вектор скорости. Скалярное произведение скорости на градиент величины / и дивергенция скорости имеют вид
л, д/ д/ д/
ugгad7 = и^ + и^ + , дх ду д,. дих диу ди ё1уи = —х + —- + —-.
дх ду д-
Закон сохранения массы в рамках системы КМГД записывается в виде
др + дк + дк + дк = о, дг дх ду д-
где в качестве потоков выступают скорректированные компоненты плотности потока массы
}1 = Р(и. - ^), I = х, у,
Величина корректировки пропорциональна т и записывается через пространственные производные, являющиеся, по сути, производными от потоков в уравнении Эйлера для МГД:
(2)
т I 5
р [дх т Гд
2 , В „2 рих + р + — - Вх
д
Л Л
+ — (риуих - ВуВх) + — (ри-их - В-Вх) ду д-
Ыу =- ^ — (риуих - ВуВх) + —
р [дх ду
2 , В „2
Риу + р + — - Ву
+д- (ри-иу - в^у )Р
Р||-(ри-их - В,ВХ) + -д(ри-иу - ВВу) + ^
р [ох ду д-
Законы сохранения компонент импульса имеют вид
2 , В „2
Ри2 + р + — - В2
дри ^ дтх1 ^ дТу, дт._ дп* ^ дпу- дп
• + —+
+ ■
• + ■
дг дх ду д- дх ду д-
где I = х, у, z■ Компоненты тензора Ту (/, ] = х, у, z) выражают действие силы, связанной с потоком скорректированного импульса, газодинамическим и магнитным давлением вдоль каждого направления:
Т = 1ч
12 2
}%и% + р+2 В - ¡уих- вуВх
¡1их - В,Вх
ухиу -ВуВх ухи1 - ВВ
¡уиу + р + 2 В2 - В2у ¡и - ВВу ]уи, - ВВу ¡и, + р + 2 В2 - В2
Здесь и далее первый индекс обозначает номер столбца, второй — номер строки. Тензор Пу
(/, ] = х, у, z) включает в себя тензор вязких напряжений Навье—Стокса ПЩ, пропорциональный
коэффициенту динамической вязкости ц, и тензор П ч™м, связанный с поправками КМГД, пропорциональными коэффициенту т:
-Г-Г _ ТТК5
(4)
где
П П =
г4 ди_х _ 2 _ 2 ди^ 3 дх 3 ду 3 д1
дх ду
+ дих дх д,
дь+дих
дх ду
4 ди^ _ 2 дих _ 2 ди^ 3 ду 3 дх 3 д1
диу + ди,
д1 ду
ди1 +
дх д1 диу ди1 д1 ду
4 ди^ _ 2 ди* _ 2 диу.
3 д1 3 дх 3 ду
п -Г-Г П5 т-г П5
Заметим, что Пу = П н , а диагональные члены:
П ? = Ц (2
^ - 2 а1уц), 81 3 !
I = х, у, 1.
Тензор ПУпНё в (4) записывается в виде
П
qmhd
-рихАих -Ар - ^Р + Д(Вх2) -рихАиу + А(ВуВх)
-р их А и1 + А(ВВх)
-риуАих + А(ВуВх)
л . А(В 2) Л/„2Ч
-риуАиу - Ар - + А(Ву)
-р иу А и1 + А(ВВ)
-ри1 Аих + А(ВВх) -ригАиу + А(ВВу) -р и1А и1 -Ар-А<р + А(В2)
К приращению А нужно применять правила дифференцирования, в соответствии с которыми имеем А аЬ = аА Ь + ЬА а■
Уравнения для магнитной индукции имеют вид
дВ д_п 3_Т1 ел __дтг _дТТ ,_ ХУ7 + + + _ , 1 _ л> у1
(5)
д1 дх ду д1 дх ду д1 где тензор, содержащий компоненты электрического поля, имеет вид
Ту = и]В1 - щВу, (6)
тензор, выражающий КМГД-поправку, определяется через комбинацию приращений:
гртн _
=
) = В1Аи] + и]АВ1 - В3Ьщ + и1АВу.
Уравнение для полной энергии имеет вид
дЕ дГх дЕу дГ, д0х дОу дО, д /Т1 ^ ^ ч — + —х + —- + —- + + —- + - — (П ххих + П хуиу + П х-и-) +
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.