РАЗРАБОТКА ТРЕХМЕРНОЙ ГИБРИДНОЙ МОДЕЛИ БЕССТОЛКНОВИТЕЛЬНОГО РАСШИРЕНИЯ ПЛАЗМЕННОГО ОБЛАКА В РАЗРЕЖЕННУЮ ИОНИЗОВАННУЮ ЗАМАГНИЧЕННУЮ ОКРУЖАЮЩУЮ СРЕДУ Текст научной статьи по специальности «Физика»

УДК 519.63:533.951

DOI: 10.18698/1812-3368-2021-3-112-132

РАЗРАБОТКА ТРЕХМЕРНОЙ ГИБРИДНОМ МОДЕЛИ БЕССТОЛКНОВИТЕЛЬНОГО РАСШИРЕНИЯ ПЛАЗМЕННОГО ОБЛАКА В РАЗРЕЖЕННУЮ ИОНИЗОВАННУЮ ЗАМАГНИЧЕННУЮ ОКРУЖАЮЩУЮ СРЕДУ

А.С. Дикалюк

ВНИИА, Москва, Российская Федерация ИПМех РАН, Москва, Российская Федерация

aleks.dikalyuk@gmail.com

Аннотация

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

Ключевые слова

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

сти влияние порядка точности аппроксимации чис- Поступила 02.09.2020

ленной схемы решения уравнения индукции на резуль- Принята 13.11.2020

таты численного моделирования © Автор(ы), 2021

Работа выполнена при поддержке РНФ (грант РНФ

№ 16-11-10275-П)

Введение. Некоторые физические явления, возникающие в лабораторных экспериментах или космическом пространстве, реализуются в результате взаимодействия плазменных потоков. Особый интерес представляют плазменные потоки с редкими столкновениями, которые возникают в лабораторных установках для магнитного или инерциального удержания плазмы [1], в лабораторных экспериментах по созданию бес-столкновительных ударных волн [2, 3], в процессе взаимодействия солнечного ветра с магнитосферой планет Солнечной системы [4-8], в результате взрывов и формирования остатков сверхновых [9, 10], в активных экспериментах [11].

Такие взаимодействия плазменных потоков можно моделировать с использованием кинетических, гидрокинетических (гибридных) и магнито-гидродинамических (МГД) подходов [8, 12-15]. МГД-модели используют гидродинамическое описание ионов и электронов, в рамках гибридных моделей ионы описываются кинетически, а электроны — в виде жидкости, полностью кинетические модели моделируют динамику ионов и электронов как отдельных частиц. Каждая модель характеризуется определенными временными и пространственными масштабами. МГД-модели позволяют описывать относительно медленные и крупномасштабные жидкостные процессы. Гибридные модели позволяют изучать процессы на временном масштабе порядка ионной циклотронной частоты и пространственном масштабе порядка ионного радиуса Лармора. Кинетические модели дают возможность моделировать явления с характерным временным масштабом порядка плазменной частоты и электронной циклотронной частоты и характерным пространственным масштабом порядка электронного радиуса Лармора.

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

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

а ионы — как отдельные частицы (т. е. кинетическим образом). Условия применимости гибридной модели можно сформулировать в виде выражений [1]:

л л г г г л 1 тъ±

ь « Л^ > L, с « L, г « L, Х = —, Гс = ——.

ПО \q\B

Здесь , Л, е — длины свободного пробега ионов и электронов; с, тсе — ионный и электронный радиусы Лармора; L — характеристический пространственный масштаб изучаемого явления.

В соответствии с гибридной моделью плазма полагается квазинейтральной:

X q№ - ne|qe| = 0, (1)

к=1

где qk, Пк — заряд и плотность числа ионных компонентов; пе — плотность числа электронов. Суммирование в выражении (1) выполняется по ионным компонентам к. Важно проверять применимость этого предположения для каждой решаемой задачи. Для того чтобы это сделать, необходимо оценить радиус Дебая й и сравнить его с L. Предположение справедливо, если й ^ Ь. Поскольку ионы в рамках гибридной модели рассматриваются кинетически, для них формулируется кинетическое уравнение Власова — Максвелла:

/+б /+1. /=0. дt дх Мк до

Здесь /к — функция распределения для к-й ионной компоненты; Мк — масса ионов сорта к; Р — сила Лоренца. Для решения этого уравнения используется метод частиц в ячейке (particle-in-cell, PIC) [16, 17]. В соответствии с этим методом ионы плазмы моделируются макрочастицами. Для каждой макрочастицы решаются уравнения движения

- "бы, В ^

Е -=

dxki - du kl Zke dt dt Mk

(2)

Здесь Хк1, бы — координаты и скорости макроиона I сорта к; Zк — зарядовое число ионов сорта к; Е, В — электрическое и магнитное поля.

Электромагнитные поля рассчитываются с использованием закона Ампера в приближении Дарвина (неизлучающее приближение, в рамках которого пренебрегается током смещения) и уравнения индукции (закон Фарадея) [1]:

^ 4 ж 7 1 д B - п

rot B = — J;--+ rot E = 0,

c c dt

(3)

где ] — полная плотность тока в системе. Уравнение для расчета полной плотности тока имеет вид

Ns

J = Je + Ji + Jext, Ji = E ZkenkUk, Je =~eneUe. k=1

(4)

Здесь ¡1 — плотность ионного тока; ]е — плотность тока электронов;

¡ех{ — внешние токи в системе; ик — среднемассовая скорость ионов

сорта к; ие — среднемассовая скорость электронов. Суммирование в выражении для плотности тока ионов выполняется по сортам ионов к.

Уравнение для электрического поля получается в результате использования уравнения сохранения ионного импульса, уравнения сохранения импульса для электронов, закона Ампера и закона Фарадея [1]:

V, V, E

^ Г1 + ^

I k = 1 neMk

E +

e

+ —

m

4ж Ns L Zkm \ 7 - 4ж -— El 1 + ^7" I Jk - rot B + — Jext c k=1V Mk J c

x B =

4 жe

f

m

ype +

5xß

( Ns

+ 4 Ж

X Zke div Kk +(vUe )Je k = 1

4 жд J

dt

ext +

m

Ns

Re -X

Zkm

k=1 Mk

Rk

Здесь т — масса электрона; Яе, Як — слагаемые, описывающие среднее изменение импульса частиц данного сорта в результате столкновений

с другими частицами; Кк — тензор кинетической энергии ионов. Это уравнение сформулировано в предположении о конечности массы электрона. Его можно упростить, если учесть, что т ^ Мк. Если умножить это уравнение на т и перейти к пределу т / М ^ 0, то получим

m®peE + e

4л Ns 7 ^ 4л -— X Jk - rot B + — Jext c k = 1 c

f

x B = -4жe

ype +

\

rme ^ dxp

+ 4жeRe.

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

4 ne

С

m®peEJe xB = -4%e

Уре +■

дк

aß

Л

dxp

+ 4KeRe.

У

Учитывая уравнение для шре, окончательное выражение принимает

вид

E = ■

cnee

Je, B

nee

ype +■

дп

aß

dxp

+ — Re, nee

Re = -nem X Ves(Ue ~Ü{).

Здесь ре — давление электронного газа; я^р — тензор вязких напряжений [1]. В простейшем случае это уравнение можно записать следующим образом (в системе единиц СИ с учетом выражения для электронного тока из (4) и пренебрежении вторым, третьим и четвертым слагаемыми в правой части):

(5)

E = -

Üe, B

К системе уравнений (1)-(5) можно добавить уравнение для давления газа электронов в виде [1]:

-2

% + ЦеУре + Уре Це = (у -1) —.

дt О

Однако в настоящей работе это уравнение в модель не включалось.

Таким образом, система уравнений гибридной модели, рассмотренная в работе (в системе единиц СИ), имеет вид:

Ns

= \qe\ X qn;

k = 1

Ns

J = ~\ qe\neÜe + X qknUk5 k = 1

dXki ^ dvkt qki , ..

dt

dt mki

(E + [u kl ,B ]);

(6)

dB

rot E =--, rot B = до J;

dt

E = -

Üe, B

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

- масса электронов полагается равной нулю;

- уравнения Максвелла формулируются в рамках безызлучательного (Дарвина) приближения.

Численный метод. Вычислительный алгоритм решения системы уравнений (6) следующий [18, 19]:

- уравнения движения макроионов решаются с использованием метода Бориса [20];

- на основе решения уравнений движения вычисляются плотность числа ионов и их среднемассовая скорость с использованием взвешивания макрочастиц на сетку [16, 17]; в случае структурированных расчетных сеток используется весовая функция PIC;

- в соответствии с законом Ампера вычисляется полная плотность тока;

- на основе полной плотности тока, плотности числа ионов и сред-немассовой скорости ионов вычисляется плотность тока электронов;

- с использованием условия квазинейтральности и плотности числа ионов рассчитывается плотность числа электронов;

- на основе плотности числа электронов и электронной плотности тока определяется среднемассовая скорость электронов;

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

Гибридная модель реализована с использованием структурированной декартовой сетки. Как указано выше, уравнения движения макроионов решаются с использованием алгоритма Бориса, включающего в себя коррекцию фазового угла [20, 21]. Электрическое и магнитное поля, используемые для вычисления силы Лоренца, интерполируются из узлов сетки к положению частицы на основе выражений [16]:

Здесь гп — вектор координат узла сетки п; Гр — вектор координат частицы р; W — весовая функция [16],

E(rp ) = Х W(rp - rn № );

n

(7)

B(rp) = Х W(rp - rn )B(rn).

n

W (r ) = W (x )W ( 7 )W (z ),

X il

TAT( \ 1 "77, X <H; W (x ) = <! H 1 1

0, XI ^ H.

(8)

Здесь H — длина пространственного шага сетки в соответствующем направлении. Суммирование в выражениях (7) выполняется по узлам, для которых W(rp - rn ) Ф 0.

В соответствии с методом PIC для ионов сорта к плотность числа частиц pk и плотность тока ¡¡к вычисляются по выражениям [16]:

Pk r ) = -1 S qpW(rp - Гп );

. V p (9)

¡¡к (Гп ) =—X qpàpw(rp - Гп ).

Vn p

Здесь Vn — объем, связанный с узлом п. Суммирование выполняется по всем макроионам p сорта к, для которых W(rp - rn ) ^ 0.

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

ÔB

— = rot dt

Ue, B

(10)

В настоящей работе первое требование удовлетворяется путем использования метода ограниченного переноса [22]. В соответствии с этим методом компоненты магнитного поля приписываются граням гекса-эдральных элементов сетки (Вх соответствует х-грани, т. е. грани с нормалью П = (1,0,0), Ву — у-грани, Вг — г-грани, рис. 1). Введем вектор

А =

Ue, B

или покомпонентно:

Ах — UyBz UzBy ;

А7 = UzBx - UxBz; (11)

Az = UxBy — UyBx.

Рис. 1. Иллюстрация к методу ограниченного переноса

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

о

— J BdS = ф A ■ dl.

dt

(12)

ds

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

Необходимо указать способ вычисления компоненты вектора А. Для

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

ие и В внутри ячейки относительно ее центра (в случае численной схе-

мы повышенного порядка точности). Отметим, что необходимо использовать особую форму восстановления компонент вектора магнитного поля для того, чтобы удовлетворить ограничению В = 0. В этой реализации гибридной модели достигается второй порядок точности по пространству путем использования ТУБ-восстановления ие и бездивергентного ТУБ-восстановления магнитного поля В [23].

Для определенности рассмотрим г-грань (верхнюю грань ячейки, ориентированную нормалью п = (0,0,1), см. рис. 1). Остальные грани могут быть рассмотрены аналогичным способом. Для того чтобы рассчитать компоненту Вг, ассоциированную с этой гранью в соответствии с (12), необходимо вычислить компоненты Ах и Ау. Следует отметить, что компонента Ах соотнесена с х-ребром, а компонента Ау — с у-ребром. Для реализации схемы против потока при расчете компонент Ах и Ау можно использовать различные подходы. Один из них заключается в расчете векторов В и ие в центре каждой расчетной ячейки и последующем переносе компонент вектора В в соответствии со знаком определенных компонент скорости ие в центры граней расчетных ячеек. Для расчета значений Ах и Ау на ребрах необходимо выполнить усреднение по соответствующим соседним граням. Этот подход гарантирует отсутствие нефизических осцилляций в численном решении, которое содержит ударные волны и выполнение второго требования, сформулированного выше.

Далее следует обсудить вопрос решения системы уравнений во времени [19]. Предположим, что координаты макроионов хп, электрические Еп и магнитные Вп поля известны в момент времени п, а скорости и""1/2 в момент времени п — 1/2. Сперва скорости и положения частиц рассчитываются в момент времени п +1/2 и п +1 соответственно с использованием метода с перешагиванием, реализованного с помощью алгоритма Бориса [20]:

бп+1/2 = б г1/2 + (Еп + Гбп, Вп ~])At,

у ь л> (13)

хп+1 = хп +ип+1/2 Дt.

Положения макрочастиц в момент времени п +1/2 определяются на основе хп+1 и хп. С использованием выражений (8) и (9) на основе и"+1/2 и х"+1/2 можно рассчитать р"+1/2 и /¿и+1/2. Зная Еп (или ЦЩ) можно

вычислить магнитное поле в момент времени n +1/2 с использованием метода ограниченного переноса и схемы против потока, описанных выше:

Bn+1/2 = Bn-у rot (En ). (14)

На основе рП+1/2, /и+1/2 и Bn+1/2 в соответствии с законом Ампера

и выражениями (4) и (5) рассчитываются иЩ+1/2 и En+1/2. Зная En+1/2, можно рассчитать

B n+1 = Bn -At rot (En+1/2 ).

Последняя величина, которую необходимо вычислить, En+1. Для этого необходимо знание значений Un+1 или, что эквивалентно, р"+1 и /и+1 (величина полной плотности тока J может быть рассчитана на основе Bn+1 с использованием закона Ампера). Величину р"+1 можно определить, так как известны координаты макроионов в момент времени n +1, рассчитанные с использованием (13). Для расчета Jn+1 необходимо вычислить скорости макроионов в момент времени n +1 по формулам

Of = 0,5 (иn+1/2 +и n"1/2 );

иn+1 =иn + (En+1/2 + [uf+1/2, Bn+1/2 ])At. mi

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

Постановка задачи и результаты численного моделирования. Рассмотренная вычислительная задача сформулирована на основе данных, представленных в [12, 24]. В центре расчетной области, которая имеет форму куба с длиной ребра 200 км, помещается плазменное облако, состоящее из однократно ионизованных атомов алюминия со скоростями в диапазоне значений 0...6-105 м/с. Масса облака 1000 кг, начальный радиус 2 км. Плазменное облако расширяется в фоновую плазму, состоящую из однократно ионизованных атомов кислорода. Плотность фоновой плазмы 5,2-10-16 г/см3. Вся система помещена в магнитное поле Bz = = 5-10-5 Тл. Для рассмотренных параметров альфеновское число Маха 9,7. Радиус магнитного торможения 292,6 км, радиус газодинамического торможения 77,1 км. Характеристики взаимодействующих плазменных образований, приведенные выше, приблизительно соответствуют параметрам натурных экспериментов типа Аргус [24].

Для численного моделирования процесса расширения использовалась структурированная гексаэдральная сетка. Число ячеек 100 х 100 х 100. Шаг сетки 2 км. Шаг по времени т = 3,2-10-5 с. Изначально в каждую ячейку помещалось 10 макрочастиц, моделирующих фоновую плазму. Плазменное облако описывалось с использованием ~1,4 млн макрочастиц. Начальная скорость частиц плазменного облака распределялась линейно по радиусу облака.

Результаты численного моделирования представлены ниже. Изменения распределения плотности заряда [Кл/м3] в различные моменты времени приведены на рис. 2. Можно увидеть, что в процессе расширения формируется бесстолкновительная ударная волна. В более поздние моменты времени (рис. 2, в и г) можно заметить формирование струй на полюсах области, возмущенной ударной волной. Направление распространения струй коллинеарно направлению начального магнитного поля.

в г

Рис. 2. Эволюция плотности заряда в процессе расширения облака при г = 1,510-2 с (а), t = 3,0-10-2 с (б), * = 6,0-10-2 с (в), * = 1,210-1 с (г)

Эволюции компоненты Вг [Тл] магнитного поля и компоненты ¡х [А/м2] полной плотности тока приведены на рис. 3, 4. Численные результаты, представленные на рис. 2-4, получены с использованием численной схемы второго порядка точности по пространству с использованием ограничителя вирвтЬвв, который является наименее диссипативным среди всех ТУБ-ограничителей. Следует отметить отсутствие нефизических осцилляций в центральной части области за бесстолкновительной ударной волной. В процессе моделирования также контролировалась точность выполнения условия В = 0 путем вычисления величины В / |В|| в каждой точке расчетной области в каждый момент времени. В рамках выполненных расчетов указанная величина не превышала значения 1,5-10-14.

в г

Рис. 3. Эволюция компоненты Вг [Тл] магнитного поля в процессе расширения облака при г = 1,510-2 с (а), г = 3,0-10-2 с (б), г = 6,0-10-2 с (в), г = 1,210-1 с (г)

в г

Рис. 4. Эволюция компоненты ]х [А/м2] полной плотности тока в процессе расширения облака при г = 1,5-10-2 с (а), г = 3,0-10-2 с (б), г = 6,0-10-2 с (в),

г = 1,2-10-1 с (г)

Одномерные распределения плотности заряда р [Кл/м3] (линии черного цвета) и компоненты В2 [Тл] магнитного поля (линии красного цвета) в различные моменты времени приведены на рис. 5. Распределения построены вдоль линии, проходящей центр расчетной области параллельно оси у. Сплошные кривые на рис. 5 соответствуют моменту времени г = 1,5-10-2 с; штриховые — г = 3,0-10-2 с; штрихпунктирные — г = 6,0-10-2 с; штрихпунктирные с двумя точками — г = 1,2-10-1 с. В наиболее позднюю фазу расширения плазмы, соответствующую моменту времени г = 1,2-10-1 с, амплитуда бесстолкновительной ударной волны (рассматривая значения ^-компоненты магнитного поля В) ниже, чем ее амплитуда в момент времени г = 6,0-10-2 с, что указывает на начало

р, Кл/м

9-10

,-5

6-10

,-5

,-5

3-10

0 50000 100000 150000 у, м

Рис. 5. Одномерные распределения плотности заряда р и ^-компоненты магнитного поля В в различные моменты времени

фазы торможения расширяющейся плазмы. Радиус плазменного облака в момент времени t = 6,0-10-2 с составляет Я ~ 40 км, что согласуется с приведенной выше оценкой газодинамического радиуса торможения.

Сравнение результатов расчетов, полученных с использованием различных численных схем для уравнения индукции, показано на рис. 6. Линиями черного цвета показано распределение плотности заряда р [Кл/м3], красного — В2 [Тл] магнитного поля в момент времени t = 6,0-10-2 с. Распределения построены вдоль линии, проходящей через центр расчетной области параллельно оси у. Сплошные линии соответствуют результатам расчетов, выполненных с использованием схемы первого порядка точности по пространству, штриховые — результатам расчетов, определенных с использованием численной схемы второго порядка точности по пространству и ограничителя выретЬее, штрихпунктирные — результатам, полученным с использованием численной схемы второго порядка точности по пространству и ограничителя Барта — Джесперсена. Использование схемы второго порядка точности по пространству для уравнения индукции предсказывает на ~20 % более высокую амплитуду бесстолкновительной ударной волны, чем в случае схемы первого порядка точности.

р, Кл/м3 Bz, Тл

60000 80000 100000 120000 140000 у, м

Рис. 6. Сравнение одномерных распределений плотности заряда р и ^-компоненты магнитного поля В, полученных с использованием различных численных схем для уравнения индукции

Следует отметить, что в процессе моделирования должен выполняться закон сохранения полной энергии (как сумма кинетической энергии частиц и энергии электромагнитного поля). В наиболее позднюю фазу расширения облака отношение полной энергии в системе к начальной полной энергии составляет 99,4 % для численной схемы первого порядка точности по пространству, 99,91 % для численной схемы второго порядка точности по пространству с ограничителем Барта — Джесперсена и 99,98 % для численной схемы второго порядка точности по пространству с ограничителем вырегЬее.

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

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

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

ЛИТЕРАТУРА

[1] Lipatov A.S. The hybrid multiscale simulation technology. Scientific Computation. Berlin, Heidelberg, Springer, 2002. DOI: https://doi.org/10.1007/978-3-662-05012-5

[2] Захаров Ю.П., Пономаренко А.Г., Тищенко В.Н. и др. Генерация сгустков лазерной плазмы с высокой эффективностью концентрации энергии для лабораторного моделирования бесстолкновительных ударных волн в замагниченной космической плазме. Квантовая электроника, 2016, т. 46, № 5, с. 399-405.

[3] Шайхисламов И.Ф., Захаров Ю.П., Посух В.Г. и др. Экспериментальное исследование бесстолкновительного сверхальфеновского взаимодействия взаимопроникающих плазменных потоков. Физика плазмы, 2015, т. 41, № 5, с. 434-442. DOI: https://doi.org/10.7868/S0367292115050054

[4] Modolo R., Chanteur G.M., Wahlund J.-E., et al. Plasma environment in the wake of Titan from hybrid simulation: a case study. Geophys. Res. Lett., 2007, vol. 34, iss. 24, art. L24S07. DOI: https://doi.org/10.1029/2007GL030489

[5] Modolo R., Chanteur G.M. A global hybrid model for Titan's interaction with the Kronian plasma: application to the Cassini Ta flyby. J. Geophys. Res., 2008, vol. 113, iss. A1, art. A01317. DOI: https://doi.org/10.1029/2007JA012453

[6] Jarvinen R. On ion escape from Venus. Helsinki, Finnish Meteorological Institute, 2011.

[7] Sillanpaa I., Kallio E., Jarvinen R., et al. Oxygen ions at Titan's exobase in a Voyager 1-type interaction from a hybrid simulation. J. Geophys. Res., 2007, vol. 112, iss. A12, art. A12205. DOI: https://doi.org/10.1029/2007JA012348

[8] Kallio E., Chaufray J.-Y., Modolo R., et al. Modeling of Venus, Mars, and Titan. Space Sci. Rev., 2011, vol. 162, pp. 267-307.

DOI: https://doi.org/10.1007/s11214-011-9814-8

[9] Ermishkin M.V., Surzhikov S.T. A three-dimensional numerical study of MHD interaction between supernova remnants and interstellar wind. Proc. 52nd Aerospace Sciences Meeting, 2014, AIAA 2014-0829. DOI: https://doi.org/10.2514/62014-0829

[10] Ermishkin M.V., Surzhikov S.T. A three-dimensional numerical study of supernova remnant type-IA evolution in an inhomogeneous interstellar medium. Proc. 45th AIAA Plasmadynamics and Laser Conf., 2014, AIAA 2014-2238.

DOI: https://doi.org/10.2514Z6.2014-2238

[11] Dyal P. Particle and field measurements of the Starfish diamagnetic cavity. J. Geophys. Res., 2006, vol. 111, iss. A12, art. A12211.

DOI: https://doi.org/10.1029/2006JA011827

[12] Raizer Yu.P., Surhzikov S.T. Magnetohydrodynamic description of collisionless plasma expansion in upper atmosphere. AIAA J., 1995, vol. 33, no. 3, pp. 486-490. DOI: https://doi.org/10.2514/3.12602

[13] Гуськов К.Г., Райзер Ю.П., Суржиков С.Т. Трехмерная вычислительная МГД-модель разлета плазмы в неоднородной ионизированной среде с магнитным полем. Математическое моделирование, 1992, т. 4, № 7, с. 49-66.

[14] Рахманов А.В., Суржиков С.Т. Расширение плазменного облака сложной формы в разреженной плазме с магнитным полем. Математическое моделирование, 1992, т. 4, № 7, с. 67-78.

[15] Суржиков С.Т. Бесстолкновительный разлет двухзарядного плазменного облака в разреженной замагниченной плазме. Физика плазмы, 2000, т. 26, № 9, с. 811-823.

[16] Хокни Р., Иствуд Дж. Численное моделирование методом частиц. М., Мир, 1987.

[17] Бэдсел Ч., Ленгдон А. Физика плазмы и численное моделирование. М., Энер-гоатомиздат, 1989.

[18] Nagy A.F., Balogh A., Cravens T.E. (eds.), et al. Comparative aeronomy. Space Sciences Series of ISSI, vol. 29. New York, Springer, 2008.

DOI: https://doi.org/10.1007/978-0-387-87825-6

[19] Winske D., Yin L., Omidi N., et al. Hybrid simulation codes: past, present, future: a tutorial. In: Space Plasma Simulation. Springer, 2003, pp. 136-165.

[20] Дудникова Г.И., Романов Д.В., Федорук М.П. Об алгоритмах метода частиц на неструктурированных сетках. Журн. вычисл. матем. и матем. физ., 2000, т. 40, № 1, с. 153-165.

[21] Zenitani S., Umeda T. On the Boris solver in particle-in-cell simulation. Phys. Plasmas, 2018, vol. 25, iss. 11, art. 112110. DOI: https://doi.org/10.1063/L5051077

[22] Evans C.R., Hawley J.F. Simulation of magnetohydrodynamic flows: a constrained transport method. Astrophys. J., 1988, vol. 332, pp. 659-677.

DOI: https://doi.org/10.1086/166684

[23] Balsara D.S. Second-order-accurate schemes for magnetohydrodynamics with divergence-free reconstruction. ApJS, 2004, vol. 151, no. 1, pp. 149-184.

DOI: https://doi.org/10.1086/381377

[24] Jones C.B., Doyle M.K., Berkhouse L.H., et al. Operation Argus. Washington, DNA, 1958.

Дикалюк Алексей Сергеевич — канд. физ.-мат. наук, ведущий научный сотрудник ВНИИА (Российская Федерация, 127055, Москва, ул. Сущевская, д. 22); сотрудник ИПМех РАН (Российская Федерация, 119526, Москва, пр-т Вернадского, д. 101, корп. 1).

Просьба ссылаться на эту статью следующим образом:

Дикалюк А.С. Разработка трехмерной гибридной модели бесстолкновительного расширения плазменного облака в разреженную ионизованную замагниченную окружающую среду. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2021, № 3 (96), с. 112-132.

DOI: https://doi.org/10.18698/1812-3368-2021-3-112-132

DEVELOPING A HYBRID THREE-DIMENSIONAL MODEL

OF A PLASMA CLOUD UNDERGOING COLLISIONLESS EXPANSION

INTO A RAREFIED IONISED MAGNETIZED MEDIUM

Dukhov Automatics Research Institute, Moscow, Russian Federation Ishlinsky Institute for Problems in Mechanics, Russian Academy of Sciences, Moscow, Russian Federation

The paper presents the results of developing a hybrid Hybrid model, kinetic ion three-dimensional model of collisionless interaction in description, non-radiative plasma flows. This model considers ions in kinetical approximation, magnetic terms (simulated as a set of individual particles) and field, collisionless shock describes electrons in terms of continuum mechanics waves (simulated as a fluid). We present the system of equations behind the mathematical model and the physical conditions limiting its applicability. The system includes equations describing ion motion in electromagnetic fields, the quasi-neutrality equation, equations for calculating the total current density, non-radiative Maxwell's equations, and the generalised Ohm's law. We outline a numerical method for solving our hybrid model equations and describe an algorithm for solving the system of equations over time. We focus on the numerical method for solving the induction equation, which takes possible discontinuous solutions into account and preserves the divergence-free condition for the magnetic field. The paper discusses the issues of increasing the spatial approximation accuracy for the numerical scheme used to solve the induction equation. We present numerical simulation results for collisionless expansion of a plasma cloud into a rarefied ionised gas in the presence of an external magnetic field. These results were obtained using our computer code that implements the hybrid

A.S. Dikalyuk

aleks.dikalyuk@gmail.com

Abstract

Keywords

model described. The paper demonstrates some numerical properties of the digital simulation developed, specifically, how the order of accuracy for the numerical Received 02.09.2020 scheme approximation designed to solve the induction Accepted 13.11.2020 equation affects numerical simulation results © Author(s), 2021

The study was supported by a Russian Science Foundation (RSF grant no. 16-11-10275-P)

REFERENCES

[1] Lipatov A.S. The hybrid multiscale simulation technology. Scientific Computation. Berlin, Heidelberg, Springer, 2002. DOI: https://doi.org/10.1007/978-3-662-05012-5

[2] Zakharov Yu.P., Ponomarenko A.G., Tishchenko V.N., et al. Generation of laser plasma bunches with a high efficiency of energy concentration for laboratory simulation of collisionless shock waves in magnetised cosmic plasma. Quantum Electron., 2016, vol. 46, no. 5, pp. 399-405. DOI: https://doi.org/10.1070/QEL16067

[3] Shaikhislamov I.F., Zakharov Yu.P., Posukh V.G., et al. Experimental study of collisionless super-Alfvenic interaction of interpenetrating plasma flows. Plasma Phys. Rep., 2015, vol. 41, no. 5, pp. 399-407. DOI: https://doi.org/10.1134/S1063780X15050050

[4] Modolo R., Chanteur G.M., Wahlund J.-E., et al. Plasma environment in the wake of Titan from hybrid simulation: a case study. Geophys. Res. Lett., 2007, vol. 34, iss. 24, art. L24S07. DOI: https://doi.org/10.1029/2007GL030489

[5] Modolo R., Chanteur G.M. A global hybrid model for Titan's interaction with the Kronian plasma: application to the Cassini Ta flyby. J. Geophys. Res., 2008, vol. 113, iss. A1, art. A01317. DOI: https://doi.org/10.1029/2007JA012453

[6] Jarvinen R. On ion escape from Venus. Helsinki, Finnish Meteorological Institute, 2011.

[7] Sillanpää I., Kallio E., Jarvinen R., et al. Oxygen ions at Titan's exobase in a Voyager 1-type interaction from a hybrid simulation. J. Geophys. Res., 2007, vol. 112, iss. A12, art. A12205. DOI: https://doi.org/10.1029/2007JA012348

[8] Kallio E., Chaufray J.-Y., Modolo R., et al. Modeling of Venus, Mars, and Titan. Space Sci. Rev., 2011, vol. 162, pp. 267-307.

DOI: https://doi.org/10.1007/s11214-011-9814-8

[9] Ermishkin M.V., Surzhikov S.T. A three-dimensional numerical study of MHD interaction between supernova remnants and interstellar wind. Proc. 52nd Aerospace Sciences Meeting, 2014, AIAA 2014-0829. DOI: https://doi.org/10.2514Z6.2014-0829

[10] Ermishkin M.V., Surzhikov S.T. A three-dimensional numerical study of supernova remnant type-IA evolution in an inhomogeneous interstellar medium. Proc. 45th AIAA Plasmadynamics and Laser Conf., 2014, AIAA 2014-2238.

DOI: https://doi.org/10.2514/62014-2238

[11] Dyal P. Particle and field measurements of the Starfish diamagnetic cavity. J. Geophys. Res., 2006, vol. 111, iss. A12, art. A12211.

DOI: https://doi.org/10.1029/2006JA011827

[12] Raizer Yu.P., Surhzikov S.T. Magnetohydrodynamic description of collisionless plasma expansion in upper atmosphere. AIAA J., 1995, vol. 33, no. 3, pp. 486-490. DOI: https://doi.org/10.2514Z3.12602

[13] Gus'kov K.G., Rayzer Yu.P., Surzhikov S.T. Three-dimensional computational MHD-model of plasma expansion into non-uniform medium with magnetic field. Ma-tematicheskoe modelirovanie, 1992, vol. 4, no. 7, pp. 49-66 (in Russ.).

[14] Rakhmanov A.V., Surzhikov S.T. Expanding of a plasma cloud of complicated form into rarefied plasma in the magnetic field. Matematicheskoe modelirovanie, 1992, vol. 4, no. 7, pp. 67-78 (in Russ.).

[15] Surzhikov S.T. Collisionless scattering of double charged plasma cloud in attenuated magnetized plasma. Fizika plazmy, 2000, vol. 26, no. 9, pp. 811-823 (in Russ.).

[16] Hockney R.W., Eastwood J.W. Computer simulation using particles. CRC Press, 1989.

[17] Birdsall C.K., Langdon A.B. Plasma physics via computer simulation. CRC Press, 2004.

[18] Nagy A.F., Balogh A., Cravens T.E. (eds.), et al. Comparative aeronomy. Space Sciences Series of ISSI, vol. 29. New York, Springer, 2008.

DOI: https://doi.org/10.1007/978-0-387-87825-6

[19] Winske D., Yin L., Omidi N., et al. Hybrid simulation codes: past, present, future: a tutorial. In: Space Plasma Simulation. Springer, 2003, pp. 136-165.

[20] Dudnikova G.I., Romanov D.V., Fedoruk M.P. On algorithms implementing the particles method on unstructured grids. Comput. Maths. Math. Phys., 2000, vol. 40, no. 1, pp. 147-158.

[21] Zenitani S., Umeda T. On the Boris solver in particle-in-cell simulation. Phys. Plasmas, 2018, vol. 25, iss. 11, art. 112110. DOI: https://doi.org/10.1063/L5051077

[22] Evans C.R., Hawley J.F. Simulation of magnetohydrodynamic flows: a constrained transport method. Astrophys. J., 1988, vol. 332, pp. 659-677.

DOI: https://doi.org/10.1086/166684

[23] Balsara D.S. Second-order-accurate schemes for magnetohydrodynamics with divergence-free reconstruction. ApJS, 2004, vol. 151, no. 1, pp. 149-184.

DOI: https://doi.org/10.1086/381377

[24] Jones C.B., Doyle M.K., Berkhouse L.H., et al. Operation Argus. Washington, DNA, 1958.

Dikalyuk A.S. — Cand. Sc. (Phys.-Math.), Leading Researcher, Dukhov Automatics Research Institute (Suschevskaya ul. 22, Moscow, 127055 Russian Federation); employee, Ishlinsky Institute for Problems in Mechanics, Russian Academy of Sciences (Vernadskogo prospekt 101/1, Moscow, 119526 Russian Federation).

Please cite this article in English as:

Dikalyuk A.S. Developing a hybrid three-dimensional model of a plasma cloud undergoing collisionless expansion into a rarefied ionised magnetized medium. Herald of the Bauman Moscow State Technical University, Series Natural Sciences, 2021, no. 3 (96), pp. 112-132 (in Russ.). DOI: https://doi.org/10.18698/1812-3368-2021-3-112-132