УДК 517.9:533.9
ЭЛЛИПТИЧЕСКОЕ ПРЕОБРАЗОВАНИЕ УРАВНЕНИЙ ПОЛЯ В НЕЯВНОЙ БЕЗЫЗЛУЧАТЕЛЬНОЙ МОДЕЛИ ПЛАЗМЫ
JI. В. Бородачев
(.кафедра математики) E-mail: [email protected]
Обсуждается методика построения устойчивого численного решения уравнений самосогласованного поля в дискретной модели Власова—Дарвина с неявной схемой расчета динамики заряженных частиц.
Введение
Как известно, широкий круг явлений плазмофи-зики, определяемых коллективными взаимодействиями частиц [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