CC BY
11
УДК 517.9:533.9
ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЯ ДЛЯ СОЛЕНОИДАЛЬНОГО ЭЛЕКТРИЧЕСКОГО ПОЛЯ В ДАРВИНСКОЙ МОДЕЛИ ПЛАЗМЫ
JI. В. Бородачев, Д. О. Коломиец, В. В. Литвинюк
(.кафедра математики) E-mail: boroda@afrodita.phys.msu.su; litvinyuk@bk.ru
Работа посвящена вопросам построения экономичной процедуры нахождения соленоидальной компоненты электрического поля в дискретной дарвинской модели самосогласованной плазмы. Предложенная методика позволяет снизить вычислительные затраты, связанные с использованием операции векторного разложения полей, и упростить выбор граничных условий.
Введение
Приближение Дарвина [1] является, как известно, одним из наиболее эффективных для представления внутренних электромагнитных полей низкочастотной плазмы. Формально оно может быть получено из полной системы уравнений Максвелла при отбрасывании соленоидальной компоненты тока смещения, что физически приводит к описанию мгновенного дальнодействия полей без учета излучения в виде свободных электромагнитных волн [2, 3].
Дарвинский формализм хотя и беднее полного максвелловского по описательным возможностям, тем не менее труднее поддается численной интерпретации: незапаздывающая природа дарвинского приближения обусловливает необходимость исключения из аналитического представления производных по времени, т.е. эллиптической переформулировки соответствующих уравнений, и прежде всего уравнения для соленоидальной компоненты электрического поля [2, 4] (решение уравнений для продольного (потенциального) электрического и магнитного полей не представляет принципиальных трудностей и уже было рассмотрено в ряде работ (см., напр., [5]). В результате мы получаем сложную нелинейную эллиптическую задачу с физически неочевидными граничными условиями, которая может быть численно решена лишь с помощью достаточно громоздкого итерационного процесса. В настоящей работе на примере одного из вариантов такого процесса обсуждаются проблемы практической реализации данной методики и предлагается ее развитие на основе дальнейшего преобразования уравнения для соленоидальной компоненты электрического поля. Последнее, с одной стороны, позволяет упростить формулировку граничных условий соответствующей краевой задачи, а с другой — построить простой и экономичный алгоритм ее численного решения.
Задача для соленоидального электрического поля
В соответствии с подходом, предложенным в работах [3, 5], соленоидальное электрическое поле Еа(г) (далее для краткости будем писать просто Еи) в заданной области П может быть найдено из следующей краевой задачи:
ЛЕг - <5; + (рЕг).,. г€<». ЧЕ„ = О, (1)
Е» ъгл — (Е Еп
где величины в правой части определяются следующим образом:
в = £ + цЕр + <;х В,
4тг qs з 4тг^
mscz
где s = t,e
'' msc6
s 6
сорт частиц, Ep — потенциальная составляющая поля, £ — дивергенция тензора переноса тока, fs — функция распределения частиц сорта s, определяемая уравнением Власова, ps — плотность заряда частиц сорта s.
Решение уравнения (1) прямыми методами невозможно, поскольку правая часть явно содержит величину pEv, из которой должна быть выделена соленоидальная компонента. Естественным выходом в данной ситуации представляется итерационный процесс, предложенный в работе [4].
Представим текущее значение плотности макрочастиц сорта s в виде ns = по + Sns, где Sns — флуктуация плотности частиц около среднего значения «о,г = п0,е = п0 = const.
Таким образом,
р
sr Ажсй „
/ j ... -<>ns
msc
msc
Учитывая полученное выражение для /л, приходим к итерационному процессу решения краевой задачи (1):
(АЕ% - (ю ЕЧ = С« + геП,
^"Ьп =
При всей естественности данного подхода его практическая реализация сопряжена с рядом неудобств. Во-первых, нам нужно достаточно аккуратно произвести векторное разложение правой части й на потенциальную и соленоидальную компоненты и решить вопрос задания соответствующих граничных условий. Во-вторых, на каждом шаге требуется производить аналогичное разложение величины 6(мЕи. Кроме того, следует заметить, что сама конструкция уравнения (1) допускает некоторую неоднозначность определения решения. Даже при небольшой погрешности векторного разложения правой части возможно возникновение и нарастание паразитной добавки к искомому соленоидальному электрическому полю (причины появления этой добавки подробно исследованы, например, в работе Хьюитта [6]). Способ гарантированно избавиться от нее неочевиден. По-видимому, требуется формирование определенного набора граничных условий в эллиптических задачах, решаемых при векторном разложении й на нулевой итерации и 8ц,Еи на каждом шаге итерационного процесса.
Таким образом, полученный алгоритм оказывается трудоемким как в плане постановки сопутствующих граничных условий, так и с точки зрения вычислительных затрат. Кроме того, уместно предположить, что скорость его сходимости при сильных градиентах ц и сложной структуре поля может быть достаточно невысокой. Однако существует способ обойти указанные трудности, произведя дальнейшее преобразование исходного уравнения для Еи.
Перепишем (1) в следующем виде:
АЕу = (<? + //£..) - (<? + //£.),,.
Введем функцию ф следующим образом:
-Ар7ф) = ф + цЕ0)р.
С учетом этого обозначения уравнение (1) принимает вид
А(Еи-Чф) = 0 + [лЕи.
В полученном выражении слева появилось слагаемое, которое компенсирует продольную компоненту правой части, избавляя нас от необходимости ее векторного разложения. Обозначив К = Е11 — 'Чф, получим для соленоидального поля систему
(АК - цК = в + ц^ф, ген
Аф = —УЖ", (3)
К\дй=Би\дй
Поле Еу определяется из соотношения Еу = К + Ч<!>. При этом второе уравнение системы (3) обеспечивает его соленоидальность.
Заметим, что для нас не принципиально, каково значение на границе переменных К и ф по отдельности, поскольку они связаны указанным выше соотношением. Таким образом, граничные условия для соленоидального поля будут удовлетворены при любом удобном для нас выборе граничного значения ф.
Следует отметить, что изложенная выше методология решения имеет хороший потенциал для развития. Определенные преобразования уравнений исходной дарвинской системы, идеологически сходные с описанными выше, позволяют получить краевую задачу для полного электрического поля Е, по виду близкую к (3). Таким образом, появляется возможность полностью избавиться от разделения поля на потенциальную и соленоидальную части. Однако корректность этих преобразований и нюансы численного решения полученной в результате задачи являются предметом дальнейших исследований.
Численное решение
Решение уравнений (3) прямыми методами возможно, но требует значительных вычислительных затрат, поэтому представляется целесообразным использовать итерационный процесс, причем оба уравнения по возможности следует решать совместно на каждой итерации. В случае 2.5-мерной постановки (х,у,их,иу,иг) для прямоугольной области Ьх х Ьу мы предлагаем использовать метод переменных направлений (см., напр., [7]).
В рамках данного метода система (3) сводится к следующей разностной задаче:
вдоль координаты х для узла (г,/)
(2/т - Лх + ¡¡¿)КТХ'2 + фхфп+т =
= {2/т-{\-!)ц, + Лу)Кпх-Ох, (4)
(2/т — Л * + ¡ц)Ку+1 ^ = = (2/т - (1 - /)/* + Лу)Ку - [мОуфп - Су, (5)
(2/т-Лх + /^ЖГ1/2 =
= (2/т — (1 — + Лу)К" — Сг, (6)
(2/т-Ах)фп+1^-ОхКпх+1/2 =
= (2/т + Лу)фп + ОуКу\ (7)
вдоль координаты у
= (2/т -(1-Пц + ЛХ)Ц+Х'2 - фхфп+]^ - Сх,
(8)
(2/т -Ау+ !ц)Кпу+х + цЩфп+х =
= (2/т - (1 -Пц + ЛЖ+>/2 - Су, (9)
(2/т — Лу + [ц)К"+1 =
= (2/т - (1 - /)уи + Лх)К?т - Ог, (10)
(2/т — Лу)фп+1 — БуКу+1 =
= (2/т+Лх)фп+]/2 + ОхКпх+]/2. (11)
Здесь величина / — параметр неявности слагаемого цКХгУ, Лх, Лу и Их, Оу — центрированные разностные операторы второй и первой производных по переменным х и у соответственно (см., напр., [7]).
Системы разностных уравнений (5), (6), (8) и (10) имеют обычную трехдиагональную матрицу, и поэтому неизвестные могут быть найдены методом скалярной прогонки [7]. В то же время (4) и (9) оказываются связаны соответственно с (7) и (11) и должны решаться совместно. Объединим две искомые величины (Кт и ф, т = х,у) в двухкомпо-нентный вектор, тогда получим систему векторных уравнений с матрицами 2 х 2 в качестве коэффициентов. Матрица этой системы, составленная из блоков 2x2, будет блочно-трехдиагональной, поэтому решение может быть найдено методом матричной прогонки (см., напр., [7]).
В качестве методического обоснования использования описанного итерационного метода нами были проведены экспериментальные исследования его сходимости для различных значений источников и коэффициента ц с периодическими граничными условиями и условиями Дирихле. Расчеты проводились на прямоугольной сетке 81x81 узлов в области Ьх х Ьу, Ьх = Ьу = 20. Правая часть первого уравнения системы (3) вычислялась по заранее заданному полю в соответствии с формулой й = —V х V х Е — цЕу. В качестве критерия достижения необходимой точности решения 6 мы использовали отношение Со нормы невязки системы (3) к норме ее правой части \\АКп-(1Кп-(1^фп-в\\/Щ\ и \\Афп + ЧКп\\/Щ\.
Основной вопрос, возникающий в случае применения метода переменных направлений, — выбор постоянного значения или последовательности значений итерационного параметра т. В случае использования априорно заданного постоянного т его оптимальная величина в основном определяется значениями коэффициента ц. Форма зависимости имеет сложный характер, поэтому представляется затруднительным выработать на основе эмпирических данных универсальные рекомендации по подбору параметра. С другой стороны, теоретический анализ операторов разностной системы с целью априорного формирования оптимальной последовательности значений т\ также является достаточно сложной задачей. В связи с этим нами был использован алгоритм автоматического подбора параметра непосредственно во время итерационного процесса. Его основная идея заключается в том, что на каждом шаге производится оценка отношения ¿2
норм ||0п+2 — ип\\/\\ип+2 — ип|| и в соответствии с этой величиной определяется необходимое изменение итерационного параметра. Здесь 1]п — текущее значение искомой сеточной функции, ип+2 — ее значение через две итерации с текущим значением т, а 0п+2 — ее значение в случае, если вместо этих двух итераций производится одна, но с удвоенной величиной т (подробнее см. [8]).
Как показали дальнейшие исследования с использованием автоматического алгоритма выбора шага, основными факторами, определяющими скорость сходимости, являются величина ц и параметр неявности /, традиционно принимающий значения 0.5 или 1. Нами были проведены тесты для различных вариантов переменного ц\ в периодическом случае ц = /¿о(2 + ъ\п(2жу/Ьу) ъ\п(2жх/Ьх)) и в случае условий Дирихле ц = ц$(х/Ьх+\ — ъ\п(Ажу/Ьу)), Но = 10—1 — 103. Основные результаты вычислений для 6 = Ю-4 представлены в таблице. Отличие полученного из расчетов поля Еи от априорно заданного составило ||(Щ,||с0 ~ 10^3.
Несмотря на то что в определенном диапазоне значений но вариант алгоритма с / = 0.5 не уступал по скорости сходимости алгоритму с / = 1, в целом метод переменных направлений на основе разностной схемы с / = 1 продемонстрировал лучшие показатели сходимости для широкого спектра значений ц,. При дальнейшем уменьшении величины 6 и усложнении пространственной структуры искомого поля разрыв в скорости сходимости между итерационными схемами с / = 0.5 и / = 1 существенно возрос. Тем не менее следует отметить, что использование е < 10^4 не привело к уменьшению ошибки определения поля Еи. Это объясняется ограниченной точностью разностной аппроксимации выражения Еи= К + 'Чф для данной величины шага сетки.
Скорость сходимости итерационного метода решения уравнения для Е0 в зависимости от граничных условий, коэффициента /л
и параметра /. Требуемая точность е = Ю-'1
МО Число итераций для граничных условий
периодических Дирихле
/ = 0.5 / = 1 / = 0.5 / = 1
0.1 28 28 168 102
1 118 32 252 182
10 96 60 122 80
100 670 18 212 50
1000 8 4 20 48
Ввиду перспективы применения данного алгоритма в составе самосогласованного кода для моделирования безызлучательной плазмы в дарвинском приближении был также проведен ряд тестов с наложением шума на значение параметра ц и правую часть й. Алгоритм показал хорошую устойчивость
к значениям шума (5-10% от базовых значений источников).
По результатам методических расчетов можно констатировать, что метод переменных направлений в сочетании с алгоритмом автоматического выбора итерационного параметра продемонстрировал хорошие показатели сходимости в широком диапазоне значений источников и коэффициентов разностных уравнений.
В дальнейшем представляется интересным, во-первых, все же попытаться сформировать априорно заданную последовательность значений итерационного параметра и произвести сравнение скорости сходимости построенного таким образом итерационного процесса и варианта с автоматическим изменением т и, во-вторых, изучить возможности улучшения используемого автоматического алгоритма выбора параметра.
Литература
1. Darwin С.G. // Phil. Mag. 1920. 39. P. 537.
2. Нильсон К., Льюис Г. // Управляемый термоядерный синтез. М„ 1980. С. 395.
3. Бородачев Л.В. Препринт. Физический факультет МГУ. № 19/2000. М., 2000.
4. Бородачев Л.В. Матем. моделирование. 2005. 17, №9. С. 53.
5. Бородачев Л.В., Мингалев И.В., Мингалев О.В. Матем. моделирование. 2006. 18, № 11. С. 117.
6. Hewett D. W. // Computer Physics Communications. 1994. 84. P. 247.
7. Самарский A.A., Николаев E.C. Методы решения сеточных уравнений. M., 1978.
8. Doss S., Miller K. // SIAM J. Numer. Anal. 1979. 16, N 5. P. 837
Поступила в редакцию 19.12.05
