Эллиптическое преобразование уравнений поля в неявной безызлучательной модели плазмы Текст научной статьи по специальности «Физика»

УДК 517.9:533.9

ЭЛЛИПТИЧЕСКОЕ ПРЕОБРАЗОВАНИЕ УРАВНЕНИЙ ПОЛЯ В НЕЯВНОЙ БЕЗЫЗЛУЧАТЕЛЬНОЙ МОДЕЛИ ПЛАЗМЫ

JI. В. Бородачев

(.кафедра математики) E-mail: boroda@afrodita.phys.msu.su

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

Введение

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

В этом классе приближений наиболее интересным представляется дарвинское (магнитоиндукци-онное) [2], которое, являясь по сути «незапаздываю-щим», описывает тем не менее ряд не свойственных «мгновенным» системам индукционных эффектов, связанных с законом Фарадея [3].

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

Однако практически эта методика успешно реализована лишь при использовании в самосогласованном алгоритме традиционной явной схемы («с перешагиванием») численного интегрирования уравнений движения частиц [4]. Применение же неявных схем интегрирования (в частности, оптимизированной схемы из работы [8]) вносит, как будет показано ниже, определенную специфику в решение проблемы актуальной численной устойчивости безызлучатель-ных моделей на пути упомянутой выше эллиптической редакции полевого описания (эффективность использования в этих целях неявной разностной ап-

проксимации исходных уравнений поля обсуждается в работе [9]).

1. Система уравнений и вопросы устойчивости

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

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

V х а; = 0, Vat = 0, а; + а^ = а.

Во-вторых, отбросить в законе Ампера поперечную составляющую тока смещения ¿:<ЭЕt/dt.

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

\7Е = 4тг р, VEt = 0, (1)

1 (9F5

V х Е =---—, VxEl = 0, (2)

с at

VxB = -j + -^, VB = 0, (3)

с с at

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

p(r,t) = Y,4jS(r-rj(t)),

Л (4)

j

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

| + Vj, = 0. (5)

При этом кинетические уравнения Власова, описывающие эволюцию распределения частиц, как

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

дг j ~dt

dvj

~at

3L

ГПп

E ■

1

(vj x В)

(6)

Здесь и ниже mj, qj, — соответственно масса, заряд и скорость частицы с индексом .

Легко видеть, что от точного полевого описания магнитоиндукционное фактически отличается лишь опущенной производной по времени от поперечного электрического поля, что, однако, приводит к весьма существенным для корректной численной интерпретации последствиям: оставаясь формально гиперболическими, дарвинские уравнения отвечают полям с мгновенным дальнодействием, т.е. соответствуют системам без запаздывания. (Доказательство этого неочевидного, на первый взгляд, утверждения можно найти в работе [5].)

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

Теоретически этот эффект легко объяснить, если воспользоваться критерием устойчивости Куранта [10], который в нашем случае имеет вид

V/ < = 1г/т,

где V} — скорость распространения электромагнитных возмущений; Л, г — величины сеточных шагов, определяемые характерными пространственно-временными масштабами рассматриваемого процесса. Очевидно, что в условиях мгновенного дальнодействия сил (г;/ = оо) требуемое отношение не может быть удовлетворено ни при каких реальных значениях Л, т. Иными словами, передача достоверной информации через пространственную ячейку с конечной скоростью по сути не возможна в системах без запаздывания.

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

2. Эллиптическая переформулировка

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

^ „ 15В VxEi = --;

с dt

„ ^ 4тг. V х В = —jt,

с

4тг)i

m

dt

= 0.

(8) (9)

Уравнение (8) уже является эллиптическим, поэтому в контексте рассматриваемой переформулировки нас будут интересовать уравнения (7) и (9). Последнее является следствием уравнения непрерывности заряда (5) и носит поверочный характер: его невязка, определяющая внутреннюю согласованность магнитоиндукционного алгоритма, должна равняться нулю (дивергентно-свободной функции в общем случае). При этом численное решение уравнения (9) подразумевает лишь необходимую корректировку текущего значения dEi/dt или j;. Для большинства алгоритмов, использующих схему «с перешагиванием» (или ее модификации) при интегрировании уравнений движения частиц, указанная корректировка носит чисто технический характер и не вызывает проблем, поскольку текущие значения полей и токов уже известны и разнесены во времени на г/2, а уравнение (9) рассматривается на полуцелых (токовых) временных слоях. Ситуация резко меняется при использовании для численного расчета динамики частиц неявной схемы, определяющей фазовые координаты, а стало быть и величины j и Е, на одном временном слое. Теперь необходимость разностного представления dEi/dt потребует экстраполяции текущих значений продольного электрического поля, провоцирующей неустойчивость того же типа, что и прямое разностное интегрирование уравнения (7), т. е. непосредственная корректировка численного решения уравнения непрерывности заряда (точнее его следствия) оказывается невозможной.

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

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

1.

дв

V(VEt)^AEt = ^-Vx ( —

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

AEt =

4тг <9j 1 д2Вi

dt

Ôt2

(Ю)

Отметим, что разностное представление уравнения (10) потребует экстраполяции текущих значений ] и Е| и, как следствие, развитие уже известной

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

Для этого распишем частную производную по времени от плотности тока (4), опираясь на уравнение движения заряда в электромагнитном поле (воспользуемся для наглядности простейшей модификацией метода макрочастиц [4]):

И <и

Л

т л к "

г,)

(Н)

3

"Ч

г,).

Заметим, что первый и второй члены справа представляют собой свертки электрического и магнитного полей соответственно с плотностями заряда и тока.

Далее представим в разностном виде <Э2Ег/<%2 для момента времени полагая, что уравнению (9) численно удовлетворено. Воспользовавшись центральными разностями [11], получим

Е;

п+1

Е?

2т

_2Е"

яП™ 1

• Е,

(12)

где члены справа уже известны (Е; из решения уравнения (1), а — из решения стандартной краевой задачи [9, Приложение]).

Обозначив (Е"^Е"-1)/т = Б_Е"/Бт, подставим выражения (11) и (12) в разностный аналог уравнения (10). Имея в виду замечание выше о физическом смысле членов выражения (11), окончательно получим в текущий момент времени

ДЕр

4 7Г

а <?а

СГПг

а»хв»ьу>ау)л

а

(13)

1 Б^Е" 4тг

Бт

Здесь правая часть имеет смысл источника поперечного электрического поля, поэтому в компактном виде ее можно записать как ^4я-£Г2С" (г).

Формальное решение уравнения (13) допускает использование прямых методов, однако экономически это оправдано лишь в одномерной пространственной геометрии (Е^Е^), где оно сводится к системе обычных скалярных уравнений [12].

Альтернативой алгебраического решения в многомерном случае является организация итерационного процесса. Представим его в форме, удобной для физической интерпретации. Для этого запишем текущее значение плотности макрочастиц сорта а (рассмотрим для простоты двухкомпонент-ную ионно-электронную плазму: а = г,е) в виде па(г) = Щ+$па(г), где 8п — флуктуация плотности частиц около среднего значения (напомним, что ПО УСЛОВИЮ КВаЗИНеЙТраЛЬНОСТИ П()г = «0е = щ)-Вводя ленгмюровскую частоту шре [3], перепишем уравнение (13) следующим образом: ,2

Д — ^

4 7Г

1 + — пц

С

Е

(14)

5П

(г) - (г) Щ Здесь ¿¿(г) = {Ча/па) г), V — номер итера-

а

ции, уточняющей текущее значение Е". Отметим, что правая часть, имеющая сложную структуру, носит нелинейный характер и учитывает взаимо-зависящие флуктуации поперечных электрических полей и токов.

Физический смысл проведенной редакции полевого описания усматривается из сравнения уравнения (14) с уравнением для экранированного (дебаев-ского) потенциала [3]: ,2

1)гр

те т

V = ~ 17Г/Л,

(15)

где р3 — нелинейный источник, обусловленный колебаниями плотности заряда около среднего значения. Тогда по аналогии с уравнением (15) полученное выше уравнение (14) описывает эффективную в масштабах скин-слоя сыр} экранировку Щ токами, обратными по отношению к порождающим магнитное поле и связанное с ним поперечное электрическое. Иными словами, предложенный подход исключает паразитную взаимоиндукцию полей и токов, определяющую актуальную численную неустойчивость традиционного дарвинского формализма.

Наконец обсудим характер и степень влияния возможной невязки разностного решения уравнения (9) на внутреннюю согласованность дискретной модели. И здесь прежде всего отметим, что основная причина невязки кроется в независимом от полевых уравнений и приближенном (интерполяционном) характере операций численного расчета плотностей тока и заряда (а следовательно, и продольного электрического поля), используемых в методе макрочастиц [4]. При этом точность указанных расчетов определяется в конечном итоге плотностью и формой

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

Вместе с тем этот эффект можно дополнительно минимизировать, потребовав (по аналогии со схемой «с перешагиванием») тождественного удовлетворения дискретного аналога уравнения непрерывности заряда (точнее, его следствия) на полуцелых временных слоях. Действительно, обозначив через ф невязку решения уравнения (9), записанного в разностном виде для момента времени

лп— 1

' 2тг (jf + j"-1) = Ф,

легко получить величину требуемой для этого поправки текущего значения плотности тока 5}? = ф/ 2тг.

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

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

Литература

9 10

И 12

Власов A.A. Теория многих частиц. М.; Л., 1950. Darwin C.G. // Phil. Mag. 1920. 39. P. 537. Франк-Каменецкий Д.А. Лекции по физике плазмы. М., 1968.

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

Нильсон К., Льюис Г. // Управляемый термоядерный синтез. М., 1980. С. 395.

Hewett D. W. // Space Sc. Rev. 1985. 42. P. 29.

DiPeso G., Hewett D. W., Simonson G.F. // J. Comp. Phys.

1994. 111. P. 237.

Бородачев Л.В. I ! ЖВМ и МФ. 1991. 30, №6. С. 934.

Бородачев Л.В. // Матем. моделир. 2005. 17, № 9. С. 53. Рихтмайер Р.Д. Разностные методы решения краевых задач. М„ 1960.

Самарский A.A. Теория разностных схем. М., 1977. Бородачев Л.В. // Вестн. Моск. ун-та. Физ. Астрон. 1993. 34, №3. С. 87.

Поступила в редакцию 01.12.04