CC BY
48
Вычислительные технологии
Том 4, № 3, 1999
ПРОБЛЕМА САМОРАЗОГРЕВА МОДЕЛЬНОЙ ПЛАЗМЫ В МЕТОДЕ ЧАСТИЦ *
В. А. Вшивков
Институт вычислительных технологий СО РАН, Новосибирск, Россия
Д. В. РОМАНОВ Новосибирский государственный университет, Россия
В. Н. Снытников
Институт катализа им. Г. К. Борескова СО РАН, Новосибирск, Россия e-mail: [email protected]
The influence of spontaneous self-force a special error of field approximation in the method of large particles on the spontaneous heating of collisionless plasma is considered.
It has been shown that for NGP, PIC and other algorithms of the method of particles self-force does not decrease with the decrease in space grid step. A new VSP class of grid nuclei has been studied free from the above-mentioned drawback and comparing well with PIC algorithm with respect to the consumption of computer resources.
Введение
Модификации метода крупных частиц для бесстолкновительной плазмы отличаются в первую очередь видом ядра частиц. В частности, предложены ядра NGP, PIC (или CIC), TSC, QS и другие [1-3] для расчетов плотности зарядов и тока в узлах сетки по координатам частиц. Из всех известных ядер при решении реальных физических задач наиболее часто используется ядро PIC вследствие хорошо изученных дисперсионных характеристик модельной плазмы, разумного соотношения между затратами и получаемой точностью расчетов вместе с относительной простотой реализации численного алгоритма.
Область применения метода частиц с PIC и другими ядрами ограничена неустранимыми погрешностями метода, к источникам которых необходимо отнести в первую очередь статистические флуктуации и шумы аналогичного происхождения. Неизбежный повышенный уровень флуктуаций плотности зарядов и других параметров модельной плазмы обусловлен малым количеством частиц, используемых в расчетах, по сравнению с реальным числом частиц в дебаевской сфере. Единственный до настоящего времени путь снижения уровня флуктуаций — увеличение количества частиц в модельной плазме.
Другой источник погрешностей в численных моделях крупных частиц связан с нефизическими эффектами, возникающими при введении пространственной сетки, на которой определяются значения электромагнитных полей, токов и зарядов. К ним относятся
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №97-02-18471.
© В. А. Вшивков, Д. В. Романов, В.Н. Снытников, 1999.
изменение потенциала взаимодействия между частицами на масштабе сеточной ячейки, погрешности аппроксимации уравнений Максвелла и погрешности интерполяции силы в местоположение частицы. Взятые вместе, указанные погрешности приводят к нарушениям в численных моделях фундаментальных законов сохранения энергии, импульса и момента импульса и, в частности, к хорошо известному стохастическому нагреву или саморазо-греву [1, 2] ансамбля частиц, который ограничивает допустимый временной интервал, отслеживаемый в расчетах.
Проблему снижения погрешностей и построения алгоритма, имеющего лучшие, в том числе и для практических целей, характеристики, чем PIC, пытались решить путем компромисса между сглаживающими свойствами ядра частицы и числом необходимых операций [1, 3], и поиска методов, подавляющих развитие нефизических численных неустойчивостей [2, 4]. Кроме того, пытались находить точное аналитическое решение уравнений движения в полях, определяемых на сетке, с тем, чтобы максимально снизить стохастический нагрев [5].
Однако существует еще один источник погрешностей, который в методе частиц связан с сеткой, — это появление самосилы, то есть силы, с которой частица действует на саму себя через пространственную сетку. Ранее [1-6] этой погрешности не уделялось должного внимания, так как в расчетах влияние самосилы на динамику отдельной частицы в ансамбле большого числа частиц обычно маскируется действием самосогласованного поля.
1. Самосила в методе частиц
Поясним появление самосилы на примере одной неподвижной частицы в неограниченном трехмерном объеме. Пусть частица с координатой x находится в ячейке между двумя узлами сетки на пространственной оси . Рассмотрим электростатическую задачу, в которой в качестве силы присутствует только электрическое поле , потенциал которого ip определяется из уравнения Пуассона
Др = —4 пр
с нулевым значением потенциала на бесконечности, где р — плотность заряда плазмы. Изменение функции распределения f частиц в ансамбле описывается уравнением Власова
df df „df n
Тй + V Я---+ q^^T =
dt dr dv
Если использовать известные методы NGP и PIC с ядром для PIC
R(x) = { 1 —ljh ' |x|- h
0, | x | > h,
то заряд q частицы переносится в один узел по модели NGP и распределяется между двумя узлами по модели PIC обратно пропорционально расстоянию до узлов. Таким образом, заряды в узлах сетки с номерами i и i + 1 будут равны
q, NGP, 0, NGP,
qth^—£ PIC, Р‘+‘ =( q|, PIC.
Потенциал в модели NGP можно записать в виде
Ю = —Cq/r,
где г — расстояние от точки г. Во втором, Р1С-варианте, потенциал равен
ю = —Ср1/г1 — Срг+1/гг+Ъ
где г и Г{+1 — расстояния соответственно от точек г и г + 1.
Таким образом, на частицу действует сила, направленная вдоль оси :
дю
Е = = —^.
дх
Для NGP-метода она равна Е = q2/х2 и для Р1С-метода
F = aJ^
\xz
pi+1
(h — x)
Из этих формул следует, что сила самодействия частицы тем больше, чем меньше расстояние между узлом сетки и частицей. Благодаря сетке частица рассеивается на самой себе с появлением известной в электродинамике особенности Е ~ 1/г2 в методах с обоими ядрами. К примеру, для Р1С-ядра: если частица с координатой х и зарядом q находится между узлами с номерами г — 1 и г на расстоянии К/3 от узла г — 1, то она передает в этот узел 2/3 своего заряда q, а в узел с номером г — (1/3^. Между частицей и узлами сетки действует сила отталкивания: узел г — 1 отталкивает частицу в 8 раз сильнее узла г. Траектория движения одиночной частицы, у которой была задана начальная скорость вдоль оси У, приведена на рис. 1 на плоскости X, У. Ясно, что такая модельная частица может совершать еще более сложные движения, особенно в трехмерном случае, возбуждать соответствующие волны и вступать с ними в резонансное взаимодействие с возбуждением коллективных степеней свободы.
Рис. 1.
Таким образом, при точном решении на сетке уравнения Пуассона и затем точном определении величины электрического поля в точках нахождения частиц на каждую из этих частиц действует посторонняя и жестко привязанная к самим частицам сила от собственного, распределенного по пространственной сетке заряда. Очевидно, что чем большее число узлов затрагивается при распределении заряда частицы по сетке, последовательно увеличивающееся при переходе от ядра NGP к PIC, к TSC и так далее, тем меньше величина самосилы. Поэтому можно связать увеличение времени саморазогрева в ансамбле частиц с изменением числа точек в шаблоне ядра частицы. При этом и для любых ранее предложенных ядер, и при фильтрации сеточной плотности зарядов, и при подстройке потенциала в k-пространстве [1-5] такая локальная на уровне ячейки несогласованность, как самосила, остается и продолжает влиять на выполнение законов сохранения.
2. Новый класс сеточных ядер частиц
В этой работе нами предлагается новый подход к построению моделей крупных частиц, в которых вклад в нефизические эффекты от самосилы существенно снижен по сравнению с моделями типа NGP или PIC с возможностью использования данных моделей в практических приложениях. Основная идея заключается в построении качественно иного класса сеточных ядер, используемых при определении на сетке плотностей зарядов и токов с тем, чтобы унаследовать простоту и экономичность реализации ядер NGP или PIC. Следуя традиции [1-3], новый класс ядер был нами назван VSP. Идейно VSP-ядра являют собой проекцию на сетку некоторой заряженной оболочки, а не сплошного равномерно заряженного тела, как для NGP- или PIC-ядер. Результатом такого проецирования становится отсутствие особенности ^ ~ q/r для потенциала точечного заряда в уравнениях
Максвелла или для дальнодействующих сил, определяемых по уравнению Пуассона для
любого положения частицы относительно узлов сетки. Другим очевидным следствием выбора частиц в виде оболочки является радикальное изменение потенциала взаимодействия между частицами в сторону ослабления сил на масштабах, меньших размера ячейки, что улучшает бесстолкновительные свойства ансамбля модельных частиц.
Возьмем одно из новых ядер VSP, которое описывается формулой
|x| , ,
2h2, |x| * h,
R(x) = ' h — , h < |x| * 2h, (2'1)
0, |x| > 2h.
Графики всех ядер NGP, PIC и VSP представлены на рис. 2. Для 2D и 3D случаев ядра могут быть получены из (2.1) по соотношениям:
R(x,y) = R(x)R(y),
R(x,y,z) = R(x)R(y)R(z). (2.2)
NGP
PIC
Безусловно, оболочечный вид ядер для плоского или пространственного случаев, в том числе и для неравномерной непрямоугольной сетки может иметь собственный вид, который не получается прямым перемножением одномерных ядер (2.1) по формулам (2.2). Однако для любых ядер необходимо потребовать симметричность расположения вторичных зарядов в узлах сетки относительно точки нахождения частицы при ее произвольном положении в ячейке, что должно привести к отсутствию особенности для локального потенциала p ~ q/r в точке нахождения частицы.
Будем также рассматривать сеточное ядро DSP, которое является линейной комбинацией PIC и VSP с параметром a для 2D и 3D численных кодов:
R(x,y) = aRi(x,y) + (1 — a)R2(x,y), з)
R(x, y, z) = aR1(x, y, z) + (1 — a)R2(x, y, z),
где R1 — ядро PIC-метода, R2 — ядро VSP-метода. Сеточное ядро DSP, как и VSP, обладает четырехточечным шаблоном и близко при a ~ 0.5 по форме к трехточечному ядру TSC [1], отличаясь от последнего экономичностью в вычислительном отношении.
3. Электростатическая одномерная задача
Рассмотрим в первую очередь использование УБР-ядра на примере электростатического случая в сравнении с Р1С-ядром. Для этого рассчитаем зависимость самосилы от параметров сеточного ядра частиц в конкретной численной модели в конечной области пространства. Ограничимся одномерной Ш электростатической задачей, в которой анализ доводится до конца.
Пусть в области решения [0, Ь] введена равномерная сетка с шагом Н. I — количество ячеек сетки (Ь = Н1). В этой области находятся частицы, которые после восстановления плотности заряда в узлах сетки по какому-либо методу дают значения р^, г = 1,..., I — 1. Для определения потенциала в узлах сетки необходимо решить разностное уравнение Пуассона
Р+1 — 2рг + р— = —4пръН2, ро = А, рЬ = В.
Решение этого уравнения можно записать аналитически явным образом:
I — г г . h2
Pi = —^A + iB + 4п—
i-i
(I — г)£ kpk + г Ys(I — k)Pk
k= 1 k=i+1
(3.1)
Рассмотрим одну частицу с координатой x, которая находится между узлами г — 1 и г, xi-1 < x < xi. Тогда 8 = (x — xi-1)/h — относительное расстояние между частицей и узлом г — 1. В методе PIC частица при вычислении плотности в узлах сетки отдает узлу г — 1 часть заряда, равную (1 — 8)q, а узлу г — 8q, где q — заряд частицы. Делением на шаг сетки получаем плотность заряда в узлах
Pi-1 = (1 — 8)q/h, pi = 8q/h.
В остальных узлах плотность равна нулю. Подставляя эти плотности в (3.1), найдем значения потенциала для A = В = 0:
h2
Pi-1 = 4п— [(I — г — 1)(г — 1)Pi-1 + (г — 1)(I — г)Pi] ,
Pi = 4п— [(г — 1)р*-1 + грг] (I — г).
Электрическое поле на отрезке [хг-1,хг] определяется по формуле
Ех = —(рг — рг-1)/Н (3-2)
и равно
Н
Ех = 4пу [(г — 1)рг-1 — (1 — г)рг] =
= 4п| [г — 1 — (I — 1)5].
Точное аналитическое значение электрического поля для частицы в точке х ее нахождения, которое определяется по силе, действующей на частицу
дЕо = д(Ех+о + Ех-о)/2,
равно
Ео = 4п\(х — Ь/2) = 4п| [5 + г — 1 — I/2]. (3.3)
Ь I
Ошибка в определении силы, действующей на частицу, равна
ДЕ = Ех — Ео = 4п^ 2 — 5
При движении частицы к левой границе ячейки с уменьшением 5 ошибка ДЕ положительна, при движении к правой границе 5 ^ 1 и ошибка отрицательна. При этом она имеет конечное значение, не уменьшающееся при измельчении сетки. С точки зрения вычислительной математики, сохранение ошибки расчетов при вычислении силы с уменьшением величины пространственного шага Н представляет собой некорректность традиционных МОР- и Р1С-методов.
Рассмотрим новое ядро УБР. Ширина этого ядра в два раза больше, чем в Р1С-методе, поэтому частица дает вклад в четыре узла:
рг-2 = рг = (1 — 5)д/2Н, р— = рт = 5д/2Н.
Для частицы с координатой х в ячейке [хг-1, хг] электрическое поле равно
Ех = 4п| [5 + г — 1 — I/2],
что соответствует точному аналитическому выражению (3.3) для силы, действующей на заряд. Решение уравнения Пуассона для всего ансамбля частиц находится очевидным суммированием по частицам, исходя из принципа суперпозиции. Ясно, что большего в одномерном случае добиться невозможно, а данный алгоритм обеспечивает аппроксимацию силы — необходимое условие корректности расчетов. Для одномерных физических моделей и методологических задач может быть также полезен и полностью согласованный алгоритм метода частиц на основе соотношений (3.1) и (3.2) с УБР-ядром (2.1).
4. Нестационарная многомерная задача и результаты расчетов
Приведенный выше анализ алгоритма легко проводится только в одномерном случае, когда просто построить и проанализировать аналитическое и разностное решения уравнения Пуассона. Для плоского 2D или для пространственного 3D случаев воспользуемся качественными физическими соображениями и анализом различных алгоритмов [1, 7].
Если в модели крупных частиц используются NGP- или PIC-ядра, то при нахождении заряда около узла сетки практически весь заряд передается в этот узел. В трехточечном ядре TSC в этом же случае заряды узлов оказываются равными 1/4, 1/2, 1/4. Поэтому в NGP-, PIC-, TSC-методах частицы эффективно рассеиваются на своих же сеточных проекциях. Только для TSC предельное значение рассеивающего потенциала оказывается в 2 раза меньшим. При этом фактически сила самодействия не зависит от шага пространственной сетки. Для ядра VSP распределение зарядов по увеличенному числу узлов сетки не просто снижает величину самосилы, но и оказывается связанным с пространственным шагом сетки. При уменьшении размера ячейки оболочка из распределенных на сетке зарядов стягивается к частице. Поэтому следует ожидать ослабления самосилы, хотя ее зависимость от значений распределенных зарядов на сетке достаточно сложна, тем более для непрямоугольных сеток. Кроме того, в отличие от 1D случая в 2D варианте самосила полностью для ядра VSP не устраняется распределением 16 зарядов на сетке с конечным шагом h, так же, как и 64 зарядами для 3D варианта.
В силу громоздкости аналитических выражений для самосилы в 2D и в 3D случаях покажем путем проведения расчетов, как меняется саморазогрев плазмы крупных частиц при смене ядер, но со всеми другими неизменными параметрами и начальными условиями, что позволит оценить влияние самосилы для различных численных моделей крупных частиц. Основной параметр, который контролировался в расчетах, — время саморазогрева. Оно определяется как время, за которое вся энергия ансамбля частиц увеличивается на некоторую фиксированную величину в диапазоне от 5 до 25%. Саморазогрев происходит, в основном за счет возрастания температуры электронов.
Расчеты проводились аналогично вошедшим в монографиях [1,2] численным экспериментам Хокни и Иствуда [7] по изучению саморазогрева плазмы крупных частиц в плоском 2D случае. Созданный нами электростатический 2D код включает NGP-, PIC-, VSP- и DSP-методы. В этом коде уравнение Пуассона решается методом преобразования Фурье. Область прямоугольная с периодическими граничными условиями по обоим направлениям. Отношение массы иона к электрону выбрано равным 100, число узлов 16 х 16, размер ячейки в обоих пространственных направлениях взяты равными 1.5, 1.0 и 0.5 от дебаев-ской длины по рассчитанной температуре распределенных в начальный момент времени электронов. Общее число крупных частиц одного типа заряда равно 4096. Тем самым число частиц в дебаевской сфере увеличивается в 2.25 и в 9 раз. Шаг по времени был постоянным и равным 0.1ц-1. Начальное распределение частиц по скоростям было выбрано максвелловским. Температура ионов равнялась температуре электронов. Распределение ионов по всей пространственной области в начальный момент времени задавалось датчиком случайных чисел, распределение электронов определялось с учетом дебаевского экранирования
[3]. Специальных мер типа задания спокойного старта в этих расчетах не предпринималось. Результаты расчетов представлены на рис. 3 - 5 в виде зависимостей относительного изменения W энергии плазмы для всех четырех моделей W =abs(w(t)/w(t = 0) — 1). Для
тъ о
Рис. 5.
размера пространственной сетки 1.0 зависимости с ядрами VSP и DSP близки друг к другу, поэтому на графике рис. 4 оставлено изменение энергии с ядром VSP. Временные затраты на расчет варианта в относительных единицах ко времени расчета варианта NGP составляют для NGP — 1.000, PIC — 1.0316, VSP — 1.032, DSP — 1.042. Расчеты для ядра DSP выполнены со значением параметра a = 0.6.
Как следует из приведенных графиков, изменение энергии ансамбля частиц от времени можно разбить на три стадии. На начальной (длительностью в один электронный плазменный период) происходит возрастание энергии до нового значения, что отражает процесс подстройки распределения электронов к начальному распределению ионов. На второй стадии с характерной длительностью ионного плазменного периода самосогласуется движение ионов с кратковременным возвращением к более низкому значению энергии по типу “электронного эха”. Третья стадия связана с увеличением энергии по степенному закону с относительно небольшим уровнем флуктуаций. Для ансамбля частиц с NGP-ядрами вторая стадия выражена слабо, а зависимости с ядрами NGP и PIC с h = 1.5 качественно повторяют результаты Р. Хокни и Дж. Иствуда [7]. Энергия плазмы в их расчетах начиная с некоторого момента времени возрастала линейно и за 1000 шагов увеличилась на 70 % для NGP-модели и на 20 % для PIC-модели (к сожалению, ряд параметров, не приведенных в [7], не позволяет произвести точное сопоставление результатов). Подтверждается также и их вывод об увеличении времени саморазогрева при переходе на ядра с большим числом точек в шаблоне или при сглаживании сеточной плотности.
Линейное изменение энергии плазмы W со временем дало основание авторам [7] принять стохастическое накопление погрешностей расчета за причину увеличения энергии.
Однако расчет, представленный на рис. 5, с PIC-ядром и h = 0.5 показал вместо линейного закона степенную зависимость W ~ t1/5. Для нового VSP-ядра фактическое изменение энергии начиная с h ~ 1.0 мало чувствительно к шагу пространственной сетки с зависимостью W ~ ta и показателем степени а ~ 1/3 -1/4. Изменение энергии ансамбля частиц с DSP-ядром также далеко от стохастического линейного закона. Для DSP-ядра с h =1.5 начальная стадия изменения энергии близка к поведению энергии с VSP-ядром, а для h = 0.5 она близка к PIC-ядру. На третьей стадии полная энергия ансамбля частиц с DSP-ядром с h = 0.5 слабо падает при увеличенном уровне флуктуаций.
Полученные зависимости не соответствуют стохастическому рассеянию частиц на своих проекциях на сетку. Так как ядро DSP изменяет для однородной плазмы величину самосилы с непрерывным переходом от PIC-ядра к VSP при неизменных физических параметрах плазмы, то эти зависимости следует интерпретировать как коллективное следствие действия самосилы с возбуждением колебаний плазмы при нескомпенсированном взаимодействии частиц с пространственной сеткой. Это взаимодействие наиболее сильно проявляется при h > Ad, когда оно не экранируется флуктуациями заряда внутри дебаев-ской сферы. В PIC-алгоритме самосила (т. е. взаимодействие частиц с пространственной сеткой) максимально компенсируется флуктуациями заряда при h ~ Ad, как было найдено опытным путем в многочисленных работах. Поэтому расчеты обычно проводились при последнем условии. Следовательно, для 2D расчетов с изменением локальной плотности плазмы и, тем самым, с изменением локального отношения шага пространственной сетки к дебаевской длине, наиболее предпочтительными характеристиками обладает алгоритм с VSP-ядром.
На рис. 6 приведены аналогичные результаты для изменения полной энергии ансамбля крупных частиц во времени для 3D электростатического кода с полностью периодическими граничными условиями. Расчеты проведены для отношения масс иона к электрону 100 на сетке 16 х 16 х 16 с шагом, равном дебаевской длине. Шаг по времени был постоянным и равным 0.1ц-1, общее число отслеживаемых частиц достигало 221184. В схемах DSP параметр а = 0.5 и 0.25. В начальный момент времени частицы были распределены по условию спокойного старта с учетом дебаевского экранирования. Расчеты проводились на Pentium-100 с оперативной памятью 16 Мб и 64-разрядной арифметикой. Временные затраты на расчеты отличались от варианта к варианту не более 0.1%. Из анализа этого графика видно, что переход к более точной аппроксимации силы, действующей на частицу, с использованием нового типа ядра VSP приводит к увеличению времени саморазогрева.
1000.00
щ%
0.01
Увеличение самосилы при переходе к ядру DSP с ухудшенной симметрией распределения заряда по узлам сетки уменьшает время саморазогрева.
5. Заключение
Таким образом, в данной работе показано, что саморазогрев ансамбля крупных частиц связан, в первую очередь, не со стохастическими свойствами системы, а с недостатками определения функции плотности зарядов на сетке, с нефизическим вкладом в действующую силу на отдельную частицу в методах типа NGP и PIC. Это дает основание для дальнейшего развития метода частиц с контролируемым или точным выполнением фундаментальных законов сохранения в численных моделях. Предложен новый класс VSP-ядер для решения уравнения Власова, которое позволяет получать решения с уменьшением численной ошибки при уменьшении пространственной ячейки. Новый метод найдет свое применение для задач пылевой плазмы, в задачах многофазных реагирующих сред и в других приложениях, где особенно важно точное восстановление траекторий частиц в самосогласованных полях.
Список литературы
[1] Хокни Р., Иствуд Дж. Численное моделирование методом частиц. Мир, М., 1987.
[2] БэдсЕл Ч., ЛЕнгдон А. Физика плазмы и численное моделирование. Энергоатом-издат, М., 1989.
[3] Березин Ю. А., Вшивков В. А. Метод частиц в динамике разреженной плазмы. Наука, Новосибирск, 1980.
[4] Сигов Ю. С. Численные методы кинетической теории плазмы. Изд-во МИФИ, М., 1984.
[5] Снытников В. Н. Численное моделирование нестационарных возмущений в плазменных потоках: Дис. ... канд. физ.-матем. наук. ИТПМ СО АН СССР, Новосибирск, 1990.
[6] Вшивков В. А. Аппроксимационные свойства метода частиц в ячейках. Журн. вы-числ. матем. и матем. физ., 36, №4, 1996, 106-113.
[7] HocKNEY R. W. Measurements of collision and heating times in a two-dimensional thermal computer plasma. J. Comput. Phys., 8, 1971, 19-44.
Поступила в редакцию 24 августа 1998 г., в переработанном виде 25 сентября 1998 г.
