Новый метод численного решения граничной задачи для системы стационарных уравнений власоваи его применение для моделирования тонкого токового слоя в хвосте магнитосферы Земли Текст научной статьи по специальности «Физика»

УДК 533.95

М. Н. Мельник, О. В. Мингалев, И. В. Мингалев, П. В. Сецко, Т. Г. Когай

НОВЫЙ МЕТОД ЧИСЛЕННОГО РЕШЕНИЯ ГРАНИЧНОЙ ЗАДАЧИ ДЛЯ СИСТЕМЫ СТАЦИОНАРНЫХ УРАВНЕНИЙ ВЛАСОВА И ЕГО ПРИМЕНЕНИЕ ДЛЯ МОДЕЛИРОВАНИЯ ТОНКОГО ТОКОВОГО СЛОЯ В ХВОСТЕ МАГНИТОСФЕРЫ ЗЕМЛИ

Аннотация

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

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

уравнения Власова, граничная задача, тонкий токовый слой, магнитосфера Земли.

M. N. Melnik, O. V. Mingalev, I. V. Mingalev, P. V. Setsko, T. G. Kogai

A NEW METHOD OF THE NUMERICAL SOLUTION OF THE BOUNDARY PROBLEM FOR THE SYSTEM OF STATIONARY EQUATIONS OF VLASOV AND ITS APPLICATION FOR SIMULATION OF A THIN CURRENT LAYER IN THE TAIL OF THE EARTH MAGNETOSPHERE

Abstract

A new method for the numerical solution of the boundary problem for the system of stationary Vlasov equations in a bounded region of space is developed. Typical examples of such problems are thin current sheets. The new method has a rigorous and clear mathematical justification, and is also convenient for efficiently organizing massively parallel computations on cluster supercomputers with the execution of the bulk of the computations on GPUs. In the new method, a fixed regular grid in the coordinate space and an adaptive (with the possibility of orientation over the magnetic field) regular grid in the velocity space with a fixed size and step, and with a center at the local hydrodynamic speed are used to approximate the distribution function.

Keywords:

Vlasov equations, boundary problem, thin current sheet, Earth magnetosphere. Введение

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

Типичными примерами таких актуальных задач являются тонкие токовые слои размерности 1D3V и 2D3V (1- или 2-мерные по пространству и 3-мерные по скоростям). Ввиду сложности системы уравнений и наблюдаемой по спутниковым данным экспериментальной картины, тонкие токовые слои в солнечном ветре и в магнитосфере Земли, а также в магнитосфере Юпитера не описываются известными точными или приближенными аналитическими решениями. Достаточно детальное воспроизведение их равновесной конфигурации может быть получено только с помощью численного решения граничной (краевой) задачи для системы уравнений, которая состоит из стационарных уравнений Власова для каждой компоненты плазмы, уравнения Ампера для магнитного поля и уравнения для стационарного электрического поля. При этом, в зависимости от специфики задачи и ее размерности, возможны различные варианты уравнения для стационарного электрического поля. В настоящее время в основном используются два типа методов численного решения нестационарной системы уравнений Власова: 1) метод крупных частиц (метод PIC в англоязычной литературе) (см., например, [1-2]); 2) сеточные методы (см., например, [3-5]). Каждый из этих методов имеет свой набор достоинств и недостатков. Но для численного решения стационарной граничной задачи применение обоих этих методов встречает ряд труднопреодолимых препятствий. Поэтому нами разработан новый метод численного решения этой задачи, который, по сути, является методом характеристик. Новый метод имеет строгое математическое и ясное физическое обоснование, а также позволяет при помощи ограниченных вычислительных ресурсов получать численные решения сложных и крупномасштабных задач, решение которых указанными выше методами либо затруднительно, либо требует на порядки больших вычислительных ресурсов.

1. Постановка граничной задачи для стационарной системы уравнений Власова

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

kxDkvV, где 1 < kx < kv < 3, где kx — пространственная размерность задачи,

kv - размерность пространства скоростей, рассматриваются на основе общей

постановки и упрощений, вытекающих из специфики конкретной задачи.

Рассмотрим плазму из K сортов ионов и электронов в ограниченной

области Q пространства M . Будем обозначать через (а;б) и [ахб] скалярное

и векторное произведения векторов в пространстве M , а через а®Ь — образованный этими векторами диадный тензор с компонентами

(a®b)kl =akb{ . Для каждой из плазменных компонент а = 1,...

обозначим через fa ( t, x, v ) — функцию распределения, где

т т

х = (x^x^xj) еК — пространственная координата, v = (v,,v>2,et —

скорость, через еа обозначим заряд частицы ( е — заряд протона, для электронов ее = —е ), через та — массу частицы, через па(х, ^) — концентрацию, через х, ?) — плотность тока. Через £0 и ¡Л0 обозначим

электрическую и магнитные постоянные, а через с = 1 ^ £0^0 — скорость света в вакууме.

Краевая задача для нерелятивистской стационарной системы уравнений Власова состоит из: 1) стационарных уравнений Власова для одночастичных функций распределения fa (х, V) каждой из плазменных компонент а;

2) формул, определяющих концентрации, плотности тока и тензор напряжений каждой из плазменных компонент а; 3) выражений для суммарных плотности заряда и тока; 4) вытекающих из уравнений Власова уравнений для моментов 0-го и 1-го порядка для каждой из плазменных компонент а; 5) системы стационарных уравнений Максвелла для плазмы; 6) граничных условий для полей и для функций распределения. В системе СИ для плазмы, состоящей из К-компонент, система стационарных уравнений Власова и формулы для моментов имеют вид

Гу; С qafw-Л.Г ...„/„М.С

' СХ

E ( x ) + [ vxB ( x )] ;-

= 0, а = \,... (1.1)

dv

(*)=J/a(*>VMV P,(X) = J]eana.(X), Р(х) = Р,(х)~ еПе(Х) (L2)

а V J

K

a=1

K

A(*) = eaJv/a(*>vKv> j(X)=j.(X)+Je(X)» (13)

К a=1

Па(х) = та J v®vfa{x,v)d\ . (1.4)

Магнитное поле определяется уравнениями Гаусса и Ампера:

divB(x) = 0, rotB(x) = fj,0 j(x) . (1.5)

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

E (x ) = -Vp(x), x ей. (1.6)

В случае пространственного разрешения модели на уровне характерного дебаевского расстояния электронов XDe = VTej($pe (здесь и далее для каждой

компоненты плазмы а через VTa =*JeTa/ma и o)pa =| ea naJ(e0ma) обозначены тепловая скорость и плазменная частота, через Ta обозначена температура в эВ), то есть при Ах ~ , выполнены материальное

уравнение вакуума

D( x ) = s0 E( x ), (1.7)

а также уравнение Пуассона

divD( x ) = p( x ), (1.8)

которое с учетом уравнений (1.7) и (1.6) принимает вид

A^(x) = -divE(x) = -р(x)/e0, x e^ . (1.9)

Наибольшим кинетическим масштабом в бесстолкновительной космической плазме является характерный гирорадиус тепловых протонов rcp = VTpj(0СР , где G)cp = eB jmp — гирочастота протонов. В плазме магнитосферы

Земли и солнечного ветра этот масштаб не менее чем на три порядка превосходит характерное дебаевское расстояние электронов, т. е. выполнена оценка

гф>1000^. (1.10)

В случае грубого пространственного разрешения Ах ~ ]Х/А:

уравнения (1.7) и (1.8) неприменимы и вектор электрической индукции D( x)

исключается из рассмотрения. В этом случае вместо уравнений (1.9) нужно рассматривать условие квазинейтральности, которое с учетом обозначения в (1.2) для плотности заряда pi( x, t) ионов примет вид

K

ene ( x ) = #( x ) = ^ ea na ( x ) , (111)

a=1

а также вытекающие из уравнений Власова (1.1) уравнения непрерывности

divja( x ) = 0, (1.12) и уравнения силового баланса (потока импульса)

eana(x)E(x)-[B(x)xja(x)]-dwIla(x) = 0 (1.13)

для каждой компоненты плазмы а = 1,... . Каждое из этих уравнений приводит к равенству

еаПа\Х)

которое определяет напряженность электрического поля как функцию от B (x) , (x) , ja(x) и П«(x) . Отметим, что полный тензор напряжений Па(x) принято представлять в виде суммы тензора инерции mana( x) ua( x)® ua( x) и

тензора давления Pa( x ):

Пa(x) = mana(x)ua(x)® ua(x) + Pa(x) , (115)

где гидродинамическая скорость ua(x) и тензор давления определяются следующими формулами:

Ua(x ) = = j V fa(x, V ) d\ , (1.16)

eana(x) na(x)

= (1.17)

R

Тогда с учетом уравнения (1.12) и формулы для дивергенции диадного тензора

div ( a ® b ) = ( a; V) b + b diva,

получим следующие формулы

divn«(x) = divpa+madiv(na ua® ua) = divpa+mana( ua;V)ua . (1.18)

Подстановка этих формул в уравнения силового баланса (1.13) приводит их к виду

eanaE^BXJaY&vYa+mana(ua'V)ua , (Х = 1,... . (1.19)

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

в (x, t\a = Bb{x, t), E (x, t)\da= Eb{x, t), x едП, (1.20)

где Bb (x, t) и Eb (x, t) — заданные функции. Отметим, что в задачах космической плазмы возможно присутствие заданных полей от внешних источников B (x) и E0 (x) , которые должны удовлетворять в области Ü уравнениям

divE0 (x) = 0, rot E0 (x) = 0, divB0 (x) = 0, rotB0 (x) = 0, x e П, (1.21) Учет этих полей можно ввести в граничные условия (1.20). Обозначим через n (x) единичную внешнюю нормаль к границе дП. Граничные условия для функций распределения состоят в том, что для каждого сорта частиц на множестве Г(+) = |(x,v) , x едП, (n (x) ;v)< 0 j задана

функция влета f^ ^ ( x, v) .

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

5= |B|, ь = B, v£ = ( •4 „■ Ч|. (1.22)

В нулевом приближении дрейфовой теории для замагниченных электронов их функция распределения f (x, v) переходит в функцию

распределения их ведущих центров Fe{r,vц 4 (см. [6]), которая зависит

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

полем B(r) декартовой системе координат. В этом случае пространство

скоростей для электронов является 2-мерным, что приводит к существенной экономии вычислительных ресурсов при численном моделировании. Стационарное уравнение Власова в дрейфовом приближении имеет вид [6]:

г,

„

дРЛ

О,

(1.23)

У

где скорость ведущего центра / \

ускорение вдоль магнитного поля

а

и модуль «ларморовского» ускорения ^ в нулевом

приближении дрейфовой теории определяются формулами, которые удобно представить в форме

„

"*е ^ ^

V, Г (Ь;УВ )

(^ ;v) ь)

V

, (1.24)

Отметим, что характеристической системой для уравнения (1.23) является 0-е приближение к системе дрейфовых уравнений Морозова — Соловьева (см. [6]), которое с учетом обозначений (1.19) и (1.20) выше имеет вид

йг ( г) _

= ие( Г (г), 1

/ л\

¿/V,,

/ / л /л

II ? II

/л\

Ш = (0).

(1.25)

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

В этом случае, согласно дрейфовой теории, плотность тока электронов определяется следующей формулой (см. [6]):

Ч / , , [ЬУ^РеЛ (1.96)

¡е(г) = -епе(и^Ъ

В

где концентрация пе (х), гидродинамическая скорость вдоль магнитного поля ие\' 4 , а также продольное 4 и поперечное ре±_(х) давление

определяются через функцию распределения ведущих центров Ее(х, 1 следующими формулами (см. [6]):

,.). 1

е{Г)=\Рг('>Л\

м

1 Г

"еУ /т

\\ ^ ~ Г 1 Г „/

-00^ 0

\ ^ 2я Г [ Г

II

(1.27)

) J

\

/

\

а

\

;

/ Л

п

/ \

/\ / \ I / \|2 Г 9 „ /

Ре\\ г г II г II г II ^

? I

(1.28)

~ Г ? I Г ^ / \ I ,

= 2лте\ V,, г „ ц

-00 V о )

( \ те Г 2Г/ \ ^ Г I Г \/ чЗ , I , (Л опл

„ „ (1.29)

м -А о )

Согласно дрейфовой теории тензор давлений электронов Ре(X) и их

гидродинамическая скорость в нулевом приближении определяются формулами (см. [6]):

Ре=ре±1 + (реп , (1-30)

где I — единичный тензор. Отметим, что подстановка в уравнение силового баланса для электронов (1.19) формул (1.30) с учетом формулы

сНу(Ре11+(Ре[] V ± Ъ'ЪЧ' „ ^

1

V

В

у

приводит уравнение силового баланса для электронов к следующему виду:

еапа(Ь;Е)Ь=-(Ь;?ре ' „ + ^

+™епе\ие\\7 'Ь У)Ь~теПе(ие'У)ие

Отметим, что в общем случае для вычисления в узлах сетки моментов функции распределения электронов, указанных в (1.2)—(1.4), нужно вычислить 10 интегралов по 3-мерному пространству скоростей. В случае описания электронов в дрейфовом приближении для вычисления этих моментов нужно вычислить 4 интеграла (1.27)—(1.29) по 2-мерному пространству скоростей. Поэтому применение дрейфового приближения для электронов в случае их замагниченности приводит к существенной экономии вычислительных ресурсов.

2. Теоретическое обоснование численного метода

Допустим, имеется классическое решение описанной выше краевой

задачи для системы уравнений Власова (1.1). Будем считать, что поля Е (X) и В (X) известны в некоторой ¿-окрестности области О (например,

продолжены постоянными вдоль нормали п (х) к границе дО). Тогда через

каждую точку ^ = е /2 х М К. фазового пространства одной

частицы проходит характеристика каждого уравнения Власова (1.1), которая является фазовой траекторией автономной системы обыкновенных дифференциальных уравнений Ньютона — Лоренца. В системе единиц СИ в

нерелятивистском случае система Ньютона — Лоренца вместе с начальными условиями в момент времени t = 0 имеет вид:

1Г='(') • тг=тг(Е(х ('))+[ V (') *в (- ('))]) •

х ( 0) = X , V ( 0) = V.

(2.1)

Решение задачи Коши (2.1) для последующего изложения удобно обозначить как

х (t ) = = ^ • XV) • V (t ) = иа( t • XV) •

(2.2)

то есть эти векторные функции удовлетворяют следующим уравнениям и

начальным условиям:

дЯ( t, X •V ) , ч ^ ) = иа( t • X •V),

= Е(Яа{^X,V)) + [и*(t• X,V) х В(Яа{t• X,V))] ),

дt та

Яа( 0, X •V ) = X • и а 0, X •V ) = V.

(2.3)

Как известно, функция распределения /а( х, V) постоянна вдоль фазовой траектории (2.2), то есть верно тождество

/а (Яа (*, X, V), и а (I, X, V)) - / ( X, V) - пош1. (2.4)

Отметим, что функция распределения вылетающих из области О частиц ^ (х,V) задана на множестве Г(_)=|(XV) , х ЕдО, (п(х)¡V)> 01,

то есть множества Г(-) и являются разложением границы фазового

пространства: д£2хК ^ ^ Для каждой точки

е и ^ рассмотрим характеристику (2.2), которая при ¿<0

уходит «в прошлое». Возможны только два случая:

1) при некотором t = —^ (XV )< 0 фазовая траектория

(Я (и XV), иа (и XV)) попадает в область влета Г(+) , на которой задана

функция влета +(х, V), а при —^ (X, V) < t < 0 она не выходит из области й:

Яа(г, X •V )еО 1 бе -1(+) (X•V)< t < 0 Яа (—1(+) (XV), X, V ) Е дО,

а(-^+) ( X, V), X, V) Е дО,

(иа (—^ ( X, V), X, V); п (Яа (—V) ( X, V), X, V))) < 0 ^

2) при любом t < 0 фазовая траектория не выходит из О :

(2.4)

Яа( Г, X, V ) е П , то есть *(+) (X, V) = .

На фазовой траектории 1-го типа в силу тождества (2.4) функция распределения равна функции влета в точке

(Яа (-^ (X, V), X, V) , иа (-(X, V), X, V) ) пересечения границы дП:

/а (X, V) = /« (Яа (-^ (X, V), X, V), иа (-^ (X, V), X, V) ) (2.5)

Для фазовой траектории 2-го типа нет общих явных условий, которые бы определяли на ней функцию распределения. В ходе численного моделирования фазовые траектории вычисляются не точно, а с некоторой погрешностью. Кроме того, стационарное решение является лишь абстрактной идеализацией, с помощью которой аппроксимируют достаточно медленно изменяющееся состояние реальной системы. Поэтому у реальной системы всегда имеется характерное время 0С существования квазистационарного решения. Отсюда следует, что в ходе численного решения достаточно рассчитывать фазовые траектории (2.2) в прошлое на время &^ ~ . Если за это время фазовая

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

Из приведенных рассуждений вытекает схема построения итерационного процесса, которая будет изложена ниже в пункте 4.

3. Дискретизация в координатном пространстве и в фазовом пространстве

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

Обозначим через {еъе2,е3} декартов базис в пространстве М. , то есть х = (х, х2, х ) = х • ех + х2 • е2 + X • е3. Будем использовать регулярную

пространственную сетку с шагом Дх и узлами

з

г (к ) = г0 + Д х • к = > ( г^ + Д х • )• е, к = *з)е2 . (3.1)

г=\

Пусть область моделирования /ЗсМ является прямоугольником с центром г 0 :

П = {х: |хг -г0г| <Ы^ • Дх = Ц, Е = 1,2,3 }.

Будем использовать сетку, покрывающую область П: ^ ={г(к), \к\ < Ыза + 2, Е = 1,2,3 }.

Шаг сетки Дуа в пространстве скоростей для частиц сорта а следует определять через характерную температуру этих частиц Тао

(в электронвольтах) по формуле А \-а = у • УТа0 , где у ~ и

УТа0 =^1 ЦаТао/ та — характерная тепловая скорость.

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

магнитным полем В(г(к)) . Декартов базис связанной с магнитным полем системы координат в узле пространственной сетки г(к) обозначим через к), к2 (к ), к3( к )}, причем будем считать, что вектор Н3(к) = В(г(к|В(г(к)) | направлен вдоль местного магнитного поля.

Для каждой компоненты плазмы а ее функцию распределения /а(х, V)

в любой момент времени ^ в каждом узле координатной сетки г(к) будем считать финитной по скоростям. Будем использовать в пространстве скоростей регулярную прямоугольную сетку фиксированного размера 2по каждому измерению , = 1,2,3 , с центром в местной гидродинамической скорости иа( г (к)), которая связана с введенными в (1.2) и (1.3) концентрацией и током этой компоненты формулой

иа(Г (к)) = ]а(Г (к))/(еаПа(Г (к))) . (3.2)

Обозначим эту сетку через

Уы{к)={Уа{Кч), -Куа1 + \<Ч1<Куа1, / = 1,2,3},

где д = {(¡I, ¿/2 7 % ) е ^ целый индекс узла сетки в пространстве скоростей.

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

Уа( к, д) = иа( Г (к) ) + ДУаХ ( Ч, - , |^ , (3.3)

а в случае ориентации по магнитному полю определяется формулой

Уа(к, д) = иа(г (*)) + ДУаХ( Ч, - ^к,(к) . (3 4)

Эта сетка должна «с запасом» содержать носитель /а(г (к), V) по скоростям, то есть на граничных и приграничных (по скоростям) слоях этой сетки /а(г(к)Уа(к,д) ) = 0 .

Функция распределения /а (х, V) заменяется конечным набором своих значений в точках

(г(к)Уа(к,д) ): /а(X V) ^ {/ап(к,д) } ,

где

fh(k,q) = fa(r (k)Уа{k,q) ), ^ < ,

+ 1 < q < ^vai, i = 1,2,3 Дискретизация в координатном пространстве состоит в замене полей и моментов функций распределения конечным набором их значений в узлах r ( k )

координатной сетки Q h :

F(х) (k) = F(r (k)), F = E, B, na, ja, Па . (3.6)

При этом интегралы по пространству скоростей в формулах (1.2), (1.3) и (1.4) заменяются соответствующими суммами:

nah (k ) = ( Ava)3 X fah(k, q ) , jah ( k ) = ea( Avaf ^Va(k, q )-fah(k, q ) ,

nah ( k ) = ma( AVa)3 ^a( k, q )® K( k, q Ш k, q ) .

q

Значение полей и всех функций, указанных в (3.6), в произвольной точке х е £2, как и в методе частиц, вычисляются по значениям в ближайших к точке х узлах сетки Qh в координатном пространстве при помощи линейного взвешивания:

Л

(3.7)

F(х)=Ss!^ )• F.(k) ,

k V J (3.8)

3

Sx( х )=П*1( X ), = max {l-|^;0}

i=1

Значение функции распределения /а( х, V) в произвольной точке (х, V) фазового пространства вычисляются по значениям (3.5) в ближайших узлах сетки (г (к) ,Уа( к, q)) в фазовом пространстве также при помощи линейного взвешивания:

а, к, ^^ ]'/а.(к, , (39)

к q V /V а /

где Бу(к,п>) = 8Х(п>) для сетки (3.3), а для сетки (3.4) взвешивание в пространстве скоростей происходит вдоль осей связанной с магнитным полем

системы координат, то есть к,) =Ц;(к))) , где (ж;(к)) —

г=1

скалярное произведение векторов ж и к1 (к) .

4. Схема итерационного процесса

Рассмотрим кратко основные детали схемы итерационного процесса

в терминах перехода от текущего приближения,

обозначенного верхним индексом , к следующему, которое обозначено

верхним индексом (к +1) , для функций распределения и функций, перечисленных в (3.6). Пусть известно приближение для полей Е(к^(к) и В^( к) . Для каждого узла (г (к) ,Уа( к, д)) сетки на множестве

(п х! и из фазового пространства выпускаем в прошлое фазовую траекторию (2.2) (которая является решением задачи Коши (2.1) ):

(Яа(t, г (к) ,Уа(к, д)) , иа(t, г (к) ,Уа(к, д) ) ) на время не более указанного в пункте 2 времени . При этом поля

в произвольной точке из О определяются через значения текущего приближения полей на сетке по формулам (3.8). В случае, когда траектория вышла из О и попала в приграничную область Г^+ ) в некоторый момент времени

= (г(к)уа(к,д)) <0, находим /(г(к),Уа(к,д)) по ф°рмуда (2.5):

/( г ( к )Уа( к, д )) = /+>(, г ( к )Уа( к, д ) ), иа(-+, г (к ) ,Уа( к, д ) ) ) .

В результате таких действий мы найдем значения функций распределения /(г (к) ,Уа( к, д)) в тех узлах сетки на множестве

(п х! и из фазового пространства, через которые проходят траектории

из области влета Г^+). Затем на основе формул (3.7) численно находим вклады найденных частей функций распределения в их моменты в узлах пространственной сетки г (к) . Как уже указывалось выше, либо значения

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

В результате находим полные значения фигурирующих в (3.7) моментов функций распределения каждой компоненты плазмы. Затем по формулам в (1.3) находим плотность полного тока (к) , по которой в результате решения краевой задачи для системы уравнений (1.5) находим в очередном приближении магнитное поле В^1 (к) . Далее из соответствующей специфике задачи

системы уравнений находим электрическое поле Е(к+1)(к) . Для контроля

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

4к+1} = 2||^(к+1) -Р(к\/(||Р^Ц + ) , Р = Е, В, где ||Р||А — какая-либо сеточная норма.

5. Применение нового метода к моделированию тонкого токового слоя в хвосте магнитосферы

Новый метод был применен к моделированию тонких токовых слоев в ближнем и среднем хвосте магнитосферы Земли. Постановка задачи, а также основные теоретические вопросы и ключевые детали моделирования изложены в работах [7-9].

На основе нового метода были созданы два варианта программы: в первом расчеты выполнялись на 8 нитях 4-ядерного процессора Intel i7 с практически 100 % распараллеливанием при помощи системы OpenMP, а во втором варианте основной объем вычислений — расчет траекторий, выполнялся на графическом процессоре (GPU) Titan 1080. Второй вариант программы продемонстрировал примерно в 20 раз более высокое быстродействие.

Результаты расчетов по сравнению с ранее выполненными расчетами на основе метода частиц подтвердили существенно лучшие свойства нового метода. В частности, он оказался априорно свободен от ряда модельных эффектов, присущих методу. Также равновесная конфигурация в новом методе достигается всего за несколько итераций, в пределах 10, в то время как в расчетах по методу частиц требовалось несколько десятков итераций. Таким образом, использование нового метода позволяет на хорошем персональном компьютере с одним или двумя современными графическими процессорами (GPU) Titan 1080 позволяет моделировать задачи, для которых в случае использования метода частиц требовался бы достаточно мощный кластерный суперкомпьютер.

6. Выводы

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

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

Благодарность. Работа выполнена при финансовой поддержке РФФИ, проект 17-01-00100, и программы ОФН РАН V.15 «Динамика разреженной плазмы в космосе и лаборатории».

Литература

1. Сигов Ю. С. Численные методы кинетической теории плазмы. М.: МФТИ, 1984. 94 с.

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

3. Eliasson B. Numerical simulations of the Fourier transformed Vlasov-Maxwell system in higher dimensions — Theory and applications // Transport Theory and Statistical Physics. 2011. Vol. 39, Issue 5-7.

4. Rieke M., Trost T., Grauer R. Coupled Vlasov and two-fluid codes on GPUs // J. Computational Physics. 2015. Vol. 283. P. 436-452.

5. Elkina N. V., Buchner J. A new conservative unsplit method for the solution of the Vlasov equation // J. Computational Phys. 2006. Vol. 213. P. 862-875.

6. Волков Т. Ф. Вопросы теории плазмы / под ред. М. А. Леонтовича. М.: Госатомиздат, 1964. Вып. 4. С. 3-19.

7. Быков А. А., Зеленый Л. М., Малова Х. В. Тройное расщепление тонкого токового слоя: новый тип плазменного равновесия // Физика плазмы. 2008. Т. 34, № 2, С. 148-155.

8. Тонкие токовые слои в бесстолкновительной плазме: равновесная структура, плазменные неустойчивости и ускорение частиц / Л. М. Зеленый и др. // Физика плазмы. 2011. Т. 37, № 2. С. 137-182.

9. Кинетические модели токовых слоев с широм магнитного поля / О. В. Мингалев // Физика плазмы. 2012. Т. 38, № 4. С. 329-344.

Сведения об авторах

Мельник Михаил Николаевич

м. н. с., Полярный геофизический институт, Апатиты

E-mail: melnik@pgia.ru

Мингалев Олег Викторович

к. ф.-м. н., зав. cектором, Полярный геофизический институт, Апатиты

E-mail: mingalev_o@pgia.ru

Мингалев Игорь Викторович,

д. ф.-м. н., зам. директора, Полярный геофизический институт, Апатиты

E-mail: mingalev_i@pgia.ru

Сецко Павел Владимирович

Программист, Полярный геофизический институт, Апатиты

E-mail: setsko@pgia.ru

Когай Татьяна Герасимовна

стажер-исследователь, Полярный геофизический институт, Апатиты