CC BY
46
Вычислительные технологии
Том 3, № 6, 1998
О МОДЕЛЯХ ЧАСТИЦ НА НЕСТРУКТУРИРОВАННЫХ СЕТКАХ
Г. И. ДудниковА, Д. В. Романов, М. П. Федорук Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: dudn@adm.ict.nsc.ru, fed@net.ict.nsc.ru
The algorithms for plasma simulation by the particle in cell method on the unstructured grid are considered. In particular, the algorithm based on the finite-size particle conception is developed. This algorithm provides the charge conservation automatically. Results of the test calculations are presented.
Введение
Известно, что хорошей исходной моделью для описания кинетических процессов в бес-столкновительной плазме является система уравнений Власова — Максвелла [1-3]:
dfa + vdfa + fa f = 0, (1)
dt dr djj
F 1 дЕ
rot B = + c.w (2)
rot E = -(3)
с dt
div E = 4np, div B = 0, (4)
p=E qa f Udp-J= E * Vafadp. (5)
a a J
Здесь fa(f,p,t) — функции распределения электронов (ионов), описывающие плотность частиц на единицу объема в фазовом пространстве в точке с координатами F и импульсами p; Fa = qa (Е +—F х B]) — сила Лоренца, действующая на частицу, qa — заряд частицы с
сорта а, Е — напряженность электрического поля, B — индукция магнитного поля, p — плотность пространственного заряда, j — плотность тока. В общем виде система уравнений (1)-(5) представляет собой очень сложную нелинейную интегродифференциальную модель, которая в реальных постановках решается численно, с привлечением ЭВМ.
Следует отметить, что наиболее универсальным и широко применяемым методом для численного решения исходной системы уравнений (1) - (5) является метод частиц. В этом
© Г. И. Дудникова, Д. В. Романов, М.П. Федорук, 1998.
методе плазма представляется набором некоторого числа модельных заряженных частиц, зависящим от мощности имеющихся в распоряжении исследователя компьютеров (как правило, в настоящее время в расчетах обычно используется 105 — 107 модельных частиц). Частицы движутся в соответствии с законами классической (или релятивистской) механики в самосогласованном электромагнитном поле, и их траектории являются характеристиками кинетического уравнения. Самосогласованное электромагнитное поле определяется из уравнений Максвелла с использованием зарядов и токов в качестве источников, а плотности этих зарядов и токов вычисляются, в свою очередь, по координатам и скоростям частиц с помощью какой-либо процедуры.
Первоначально метод частиц был развит для численного моделирования бесстолкно-вительных плазменных процессов в областях простой формы, когда в моделируемой области можно ввести ортогональную сетку: декартову, цилиндрическую или сферическую. В этих случаях можно реализовать классические варианты метода частиц, которые описаны в широко известных монографиях, посвященных моделированию плазмы этим методом [1-3]. Общепринятая классификация различных вариантов метода основана на способе определения плотности зарядов и токов по характеристикам модельных частиц. Для уменьшения уровня флуктуаций плотности заряда в узлах конечно-разностной сетки с частицей связывают облако определенной формы, отвечающее данной схеме распределения заряда. Площадь перекрытия облака с ячейкой определяет долю заряда, приписываемую сеточному узлу в центре этой ячейки. Классификацию моделей частиц можно произвести также на основе формы функции распределения заряда [2, 4] (или сеточного ядра [1]). В этом случае с частицей связывают функцию распределения и значение этой функции в каждом узле сетки дает долю заряда, приписываемую этому узлу.
Следует подчеркнуть, что в практических расчетах используются лишь несколько наиболее простых вариантов метода частиц. Так, например, в случае двумерной, равномерной декартовой сетки наиболее употребительным является PIC-метод [5]. В этом методе частицы представляют собой прямоугольники со сторонами hi, h2 (здесь hi,h2 — шаги пространственной сетки по x,y соответственно). Алгоритм распределения заряда сводится к тому, что заряд каждой частицы распределяется с помощью обратной билинейной интерполяции между четырьмя ближайшими узлами сетки. Процедуры, аналогичные схемам раздачи заряда, используются для интерполяции электромагнитных сил в местоположение частиц. Для осуществления процедур интерполяции силы и раздачи заряда необходимо определить адрес ячейки, в которой находится частица с номером к. Для равномерной сетки номер ячейки, содержащей частицу, находится следующим образом:
I = INT[(xfc — Xo)/hi] + 1, J = INT[(yfc — yo)/h2] + 1. (6)
Здесь операция INT означает взятие целой части от выражения, заключенного в квадратные скобки, а x0, y0 — координаты левого нижнего угла расчетной сетки. Для решения уравнений движения частиц, как правило, используется метод с перешагиванием [1-3]. Наконец, для решения уравнений Максвелла методом конечных разностей используется очень удобная в численной реализации схема с перешагиванием [1-3].
Вместе с тем потребности фундаментальной науки и практических приложений привели к необходимости моделирования плазменных процессов посредством метода частиц в областях сложной геометрии. Такая ситуация имеет место, в частности, в ионных диодах — установках для получения мощных ионных пучков [6]. Поэтому в последнее время получили развитие различные варианты метода частиц для моделирования динамики плазмы в областях сложной формы [7-10, 14-16]. Все эти варианты основаны на использовании в
расчетах сеток (структурированных или неструктурированных), адекватно описывающих границы моделируемой области ( в зарубежной литературе этот подход получил название boundary — fitted coordinates).
Необходимо отметить, что в современной научной литературе нет устоявшегося определения понятий структурированной и неструктурированной сеток. Наиболее общее представление о расчетной сетке соответствует понятию размерно однородного связного геометрического (топологического) комплекса [17]. При этом сетка рассматривается как набор взаимосвязанных и взаимоподчиненных геометрических элементов: узлов, ребер, граней, ячеек и т. п. Под структурированными сетками обычно подразумеваются такие сетки, для которых топологическая структура соответствующего комплекса однозначно характеризуется набором индексов узловых точек. Неструктурированной сеткой называют всякую сетку, которая не является структурированной.
Были развиты подходы, основанные на использовании структурированных сеток: треугольных [7], четырехугольных [8-10, 14]. Алгоритмы расчета основаны на отображении моделируемой (физической) области на логическую. При этом границы моделируемой области в плоскости переменных (x, y) отображаются на границы прямоугольника в плоскости переменных (£,п). В прямоугольной области вводится равномерная сетка, которая отображается на физическую область. Алгоритмы расчета, описанные в литературе, различаются по способам решения полевых уравнений и уравнений движения частиц. Так, в работе [9] был реализован подход, при котором полевые уравнения, как и уравнения движения частиц, интегрируются на логической сетке. В работе [10] для численного моделирования стационарной системы уравнений Максвелла — Лоренца был реализован алгоритм, основанный на решении уравнения Пуассона на логической сетке, а уравнений движения частиц — на физической сетке. Такой подход является более эффективным, поскольку для решения уравнений движения частиц можно использовать стандартную (простую и экономичную) схему интегрирования [11]. Однако в этом случае возникает проблема локализации частиц, интерполяции сил в их местоположение и раздачи заряда частиц в узлы неравномерной сетки, состоящей из выпуклых четырехугольников. Построению этих процедур посвящен ряд работ [12, 13]. В [14] численное интегрирование двумерных нестационарных уравнений Максвелла производилось в физическом пространстве с использованием MUSCL-схемы второго порядка аппроксимации, построенной на основе метода конечных объемов. Для локализации частиц, интерполяции сил и раздачи заряда использовались процедуры, аналогичные [13].
Однако, по нашему мнению, более предпочтительным является использование в расчетах неструктурированных сеток. Хотя структура данных при использовании таких сеток весьма сложна, они обладают рядом преимуществ по сравнению с сетками регулярной структуры. Это касается, например, численного моделирования процессов в неодносвяз-ных областях, где применение сетки регулярной структуры попросту невозможно; решения широкого класса проблем, связанных с адаптацией разностной сетки, которая, как правило, осуществляется локально, что приводит к необходимости изменения структуры сетки и т.д. Построению подобных алгоритмов посвящены работы [15, 16]. Для решения уравнений Максвелла на неструктурированной сетке используется метод конечных элементов, или метод конечных объемов [18-20]. Эффективные процедуры для локализации частиц на такой сетке, состоящей из треугольников в случае двух пространственных измерений и тетраэдров в случае трех, были построены в [15, 21-22]. Для осуществления процедур интерполяции сил и раздачи заряда частиц в сеточные узлы используются базисные функции (функции формы), определенные на каждом элементе неструктурированной
сетки [15, 21].
Ниже, в значительной мере следуя работам [15, 21], опишем алгоритмы локализации частиц, интерполяции сил в их местоположение и раздачи заряда и представим результаты тестовых расчетов в случае двух и трех пространственных измерений. Хотя эти алгоритмы раздачи заряда и тока весьма эффективны и просты в реализации, они обладают одним существенным недостатком, связанным с нарушением дискретного аналога закона сохранения заряда
(!) к+*=°- (7)
Способы преодоления этого принципиального недостатка электромагнитных кодов описаны в литературе, в основном, применительно к алгоритмам метода частиц на декартовой сетке [11, 23-26]. Так, в работе [11] для выполнения уравнения (7) было предложено на каждом временном шаге подправлять продольную часть электрического поля. В работе [23] был предложен очень простой способ коррекции разностного закона Гаусса, основанный на введении в уравнение (2) члена, пропорционального величине VГ, который получил название "псевдотока". Здесь
Г = Е - 4пр.
Улучшенная модификация этого метода и ее сравнение с "классическим"методом [11] представлены в работе [24]. В работах [25, 26] разработан электромагнитный код, в котором для выполнения разностного аналога уравнения непрерывности было предложено находить вклад в плотность тока от каждой макрочастицы, непосредственно подсчитывая потоки через стороны (грани) прямоугольных ячеек. Поскольку этот закон автоматически выполняется для отдельной макрочастицы, он, естественно, выполняется и для всей системы в целом.
В разделе 3 нами предлагается алгоритм, основанный на концепции частиц конечного размера. Для того, чтобы удовлетворить дискретному аналогу уравнения непрерывности (7), мы использовали подход, аналогичный [25], который был обобщен здесь на случай неструктурированной сетки. Представлены также результаты тестовых вычислений для случая двух пространственных измерений.
1. Алгоритм "точечная частица — неструктурированная сетка"
Приступая к описанию алгоритмов метода частиц на неструктурированной сетке, рассмо трим вначале схему интегрирования уравнений движения частиц
(гк (г)
(г
4(г), (1.1)
Ри (г), (1.2)
(г
которая была предложена в работе [11]. Далее мы рассмотрим только одну частицу и поэтому опустим индекс к. Запишем схему с перешагиванием в следующем виде (рассматривая для общности релятивистский случай):
Рг+1 - Р1 = тьп+1/2, (1.3)
(£п+1/2 - аЕп) - (ип-1/2 + аЕп) = — [£п х Вп], (1.4)
7 пс
9Т Г~Г - - ТТ
где а =-, и = — релятивистская скорость. Для алгоритмическом реализации схемы
2т
(1.3) — (1.4) введем величины
и- = ип-1/2 + аЕп, (1.5),
и+ = ип+1/2 - аЕп, (1.6),
заменим ип в правой части (1.4) ее средней величиной, а также аппроксимируем 7п, используя величину и-:
7п « 7- = 1 + ^. (1.8)
с2
Тогда формулу (1.4) можно записать в виде
и+ — и- = (и+ + и-) х Ь, (1.9)
^п и и
где Ь = —-—Ьп, Ьп = . Для определения и+ из (1.9) сначала найдем величину и
7- |-°п|
а|В п|
Ъ, Ъ = |ж\. Для Шфе/
ь! = и- + и- х Ь (1.10)
и определим и+ из формулы
2
и+ = и +--2 их £ (1.11)
1 + |Ь|2
Наконец, величина релятивистской скорости £п+1/2 определяется после второго "полуускорения" под действием электрического поля из выражения
£п+1/2 = и+ + аЕп. (1.12)
Затем из (1.3) определим положение частицы в момент времени Ь = Ьп+1. Данная схема интегрирования уравнений движения имеет второй порядок точности 0(т2) и обратима во времени.
Перейдем теперь к описанию алгоритмов локализации частиц на неструктурированной сетке. Алгоритм локализации частиц, предложенный в [21, 22], основан на вычислении линейных базисных функций [27] (функций формы) Si,a = Si,a(Xcn+1) для макрочастицы, локализованной в момент времени Ь = Ьп внутри элемента С (рис. 1). Здесь этот алгоритм описан для двумерного случая, однако он легко обобщается и на трехмерный случай. Итак, в двумерном случае Ci — треугольник с вершинами и координатами вершин а,,,,а = (аа, Ъа)Т, где 1 < а < 3. Для того, чтобы определить, будет ли частица Х?п+1 = (хп+1 ,уп+1) принадлежать элементу С,,,, используется следующий алгоритм. 1) Вычисляем величины Si)1, Si,2, Si)з:
з з
^ ^ Si,aai,a = х + , ^ ^ Si,a = 1. (1.13)
а=1 а=1
Из формул (1.13) можно получить выражения для функций формы ^,2 = 1[(Ъз - Ъ1)(хп+1 - а1) - (аз - а^у^1 - Ъ^],
Рис. 1. К алгоритму локализации частиц в двумерном случае.
1
Si,3 = - bi)(xn+1 - ai) - (a2 - a1)(yn+1 - bi)]
n+1
J
Si,1 — 1 - Si,2 - S^
4,3-
(1.14)
Здесь
3 = (Ъз - Ъ\)(а<2 - аг) - (62 - 61 )(аз - «1). (1.15)
2) Если шт(Бгд, Бг,2, Бг,3) > 0 и шах(Бгд, Бг,2, Бг,3) < 1, то в момент времени г = Г+1 частица лежит внутри элемента С г.
3) Если шт(Бгд, Бг,2, Бг,3) < 0, то хп+1 (// Сг и поиск частицы необходимо продолжить в элементе, противолежащем вершине Мг,а с наименьшей величиной функции формы Бг,а.
Алгоритм локализации частиц в трехмерном случае был предложен в работе [15]. Рассмотрим элемент (в данном случае тетраэдр) неструктурированной сетки Сг с вершинами
—* т
Мг,а, имеющими координаты аг,а = (аа,Ъа, са) , 1 < а < 4 (рис. 2, а). Вычислим четыре детерминанта Ога, описывающих положение частицы хп+1 относительно этого элемента:
где Ac
ai
Ад — [A2 x A3 ] A4, D i,2 — [A1 x A4] A3,
Di,3 — [A4 x A1 ]A2, D i,4 — [A3 x A2] A1, - xn+1. Далее возможны следующие варианты:
(1.16)
Di,k > е; Vk, 1 < k < 4, (1.17)
где е — "машинный нуль", частица в момент времени t — tn+1 находится в элементе Ci;
Dik < е, k е [1, 2, 3,4], (1.18)
xn+1 е Ci, тогда поиск частицы необходимо продолжить в элементе, противолежащем вершине Miyk;
Dik < е, Du < е; k,l е [1,2,3,4], k — l. (1.19)
В этом случае два или даже три детерминанта становятся отрицательными. Тогда мы должны определить одну из граней текущего тетраэдра Ci, через которую проходит траектория частицы. Уравнение траектории частицы имеет, в параметрической форме, следующий вид:
+ Avn+1/2. (1.20)
n
rp — x
Рис. 2. К алгоритму локализации частиц в трехмерном случае.
Для определенности предположим, что < е, Ог2 < £■ Рассмотрим треугольник
(М 2, иг>3, ЫгА)
(рис. 2, б). Введем локальный базис (а, Ь) данного треугольника и разложим вектор Гр — аг,2 по этому базису. Для нахождения локальных координат п точки гр получим уравнение
Гр — аг, 2 = £ (аг,з — аг, 2) + п(аг, 4 — а*,2), (1.21)
откуда можно определить локальные координаты точки гр. Если п £ [0,1] и £ + п < 1, то частица пересекает плоскость (Мг,2, Мг,3, Мг,4) и ее поиск необходимо продолжить в элементе Сг1; в противном случае ее поиск необходимо продолжить в элементе Сг2■
Рассмотрим далее процедуры интерполяции сил в местоположение частиц и раздачи заряда. Пусть частица с фазовыми координатами (Хк, щ) находится в элементе Сг неструктурированной сетки. Определим плотность заряда и тока, вносимые данной макрочастицей в узел г с координатами аг. Этот узел окружен конечным числом элементов иг, формирующих область Ог данного элемента (рис. 3). Как известно [1], плотность заряда и тока точечных частиц определяется по формулам
N
р(х, I) = 6[Х — Хк(£)], (1.22)
к=1
N
2 (х,г) = ^ Як$к(Рк)^[Х — Хк(£)]■ (1.23)
к=1
В этом случае заряд макрочастицы с номером к приписывается ближайшему узлу неструктурированной сетки. В случае ограниченного числа макрочастиц это, как известно, приводит к большому уровню флуктуаций плотности заряда и тока. Для уменьшения уровня
п
Рис. 3. Область, связанная с узлом г неструктурированной сетки.
флуктуаций при вычислении зарядовой и токовой плотности в узлах % используем линейные базисные функции ¿¿(ж), определенные на каждом элементе неструктурированной сетки [27]. Пусть частица локализована в некотором злементе неструктурированной сетки. Будем определять заряд и ток, вносимый частицей в узел % данного элемента, по формулам
Из определения линейных базисных функций [27] следует, что ^ г Бг (Хк) = 1. Поэтому полный заряд макрочастицы при таком способе распределения заряда сохраняется. Суммарный заряд и ток, внесенные макрочастицами в данный узел, есть
к=1
Здесь сумма берется по всем макрочастицам. Тогда плотность заряда и тока в узле % определим как
0%к — 0,к Бг(хк ), ■гк — 0к$к Бг(хк).
к=1
(1.24)
(1.25)
Здесь У обозначает связанный с узлом % объем:
(1.26)
В формуле (1.26) сумма берется по всем элементам сетки (треугольникам или тетраэдрам), примыкающим к данному узлу % (см. рис. 3), У — объем элемента с номером I, б — размерность пространства.
Наконец, силы действующие на частицу Xk (t), определяются по формулам
Ог
F (X,t) = Y, Si,a(Xk (t))F (aha,t), (1.27)
,a V
a=1
где сг = 3, 4 соответственно для двух- или трехмерного случая. Здесь Г(аг,а,£) — силы, заданные в вершинах элемента Сг.
Рассмотрим результаты некоторых тестовых расчетов, проведенных с использованием описанных алгоритмов. Сначала рассмотрим процедуры локализации и раздачи заряда в двух измерениях. Область, в которой проводились тестовые вычисления, представляет собой единичный квадрат. Примеры неструктурированных сеток, введенных в этой области, показаны на рис. 4.
Рис. 4. Примеры неструктурированных сеток в тестовых расчетах: a — 130 узлов, 222 элемента; б — 701 узел, 1320 элементов.
Для тестирования процедур локализации и раздачи заряда в вычислительной области равномерно распределялось определенное число макрочастиц Нр £ [400 — 409600]. Предположим, что плотность плазмы в расчетной области постоянна: р0 = 1. На первом этапе тестовых вычислений была проведена локализация частиц в расчетной области и вычислены функции формы (1.14) каждой макрочастицы по отношению к вершинам элементов Сг, в которых находятся частицы. Затем, в соответствии с (1.24) - (1.25), производилась раздача заряда макрочастицы в узлы и барицентры неструктурированной сетки. Барицентрические координаты гт элемента с номером т определялись следующим образом:
гт = ^ гг/б, (1.28)
где сумма берется по координатам узлов принадлежащих данному элементу т, б — размерность пространства. Результаты численных экспериментов показаны на рис. 5. На рис. 5, а представлена максимальная абсолютная величина |р0 — рг | = |1 — рг | отклонений плотности в узлах и барицентрах неструктурированной сетки в зависимости от числа частиц. Видно, что это отклонение для данной сетки составляет менее 1 %, если число частиц в
Рис. 5. Результаты тестовых расчетов по восстановлению равномерной плотности в двумерном случае (130 узлов, 222 элемента): максимальное отклонение (в процентах) от тестовой плотности в зависимости от числа частиц (а), гистограммы отклонений расчетной плотности от тестовой в узлах неструктурированной сетки (б) и в барицентрах неструктурированной сетки (в); число
макрочастиц N = 6400.
области N > 6400. Гистограммы отклонений расчетной плотности от теоретического значения для N — 6400 представлены на рис. 5, б, в. Аналогичные гистограммы для случая трех пространственных измерений и неструктурированной сетки, содержащей 200 узлов и 744 элемента, приведены на рис. 6.
Для тестирования алгоритмов локализации частиц, интерполяции сил в их местоположение и интегрирования уравнений движения была рассмотрена задача о движении пробной частицы в центрально-симметричном поле [28] вида
— , г^^——ч з — , г^^——\ з, (1-29)
-0.1гх Р — —0.1гу
(у/гх+^)3' у (УГХТ7|У
где гх — х — 0.5, гу — у — 0.5. В начальный момент времени Ь — 0 ж(0) — 0.2, у(0) —
0.5, г>х(0) — 0, уу(0) — —=. Траектория частицы в таком поле представляет собой окруж-
V 3
ность радиуса г — 0.3. Траектории частицы на фазовых плоскостях (ух,х), (уу,у) также
представляют собой окружности с радиусом гу — —р. Результаты численных расчетов для
3
неструктурированных сеток, приведенных на рис. 4, показаны на рис. 7. На них также
Рис. 6. Результаты тестовых расчетов по восстановлению равномерной плотности в трехмерном случае (200 узлов, 744 элемента): гистограммы отклонений расчетной плотности от тестовой в узлах (а) и в барицентрах неструктурированной сетки (б); число макрочастиц Нр = 64000.
представлены (для сравнения) результаты аналитического решения задачи.
2. Алгоритм "частица конечного размера — неструктурированная сетка"
Построим теперь алгоритм, в котором для восстановления плотности заряда и тока предлагается использовать макрочастицы конечного размера. В этом случае вклад заряда макрочастицы в барицентр соответствующего элемента неструктурированной сетки Ст будет пропорционален площади пересечения "частицы-облака" и данного элемента Ст. Для простоты здесь ограничимся рассмотрением только двумерной ситуации. В предлагаемом нами алгоритме частице можно придать форму любого выпуклого многоугольника. Однако структура данных наиболее проста для частицы прямоугольной формы. В этом случае каждая макрочастица определяется заданием координат ее левого нижнего угла и сторонами к\,к2. Основным элементом алгоритма является разработанная нами процедура определения площади пересечения "частицы-выпуклого многоугольника" с элементом неструктурированной сетки. Последовательность шагов, реализованная в этой процедуре, изображена на рис. 8. Образно говоря, мы берем "ножницы" и последовательно режем ими вдоль каждой стороны "элемента-треугольника" и отбрасываем ту часть "частицы-прямоугольника", которая остается снаружи от линии разреза. В итоге (после трех разрезов) получаем выпуклый многоугольник, который является результатом пересечения макрочастицы с элементом неструктурированной сетки.
Вклад частицы с номером р в плотность заряда элемента с номером т определяется следующим образом:
Б ■
Арт = -т Рр, (2Л)
Бт
где Бт,г — площадь пересечения частицы и элемента, Бт — площадь элемента с номером т, рр — плотность заряда частицы (здесь она предполагается равномерной).
Рис. 7. Движение "пробной частицы" в центрально-симметричном поле: I — 130 узлов, 222 элемента, II — 701 узел, 1320 элементов. а — расчет, б — аналитическое решение.
Рис. 8. Последовательность шагов при определении площади пересечения частицы конечного размера и элемента неструктурированной сетки.
Силу, действующую на частицу с номером p, определим равенством
j = ^ m Sm,jFm (2 2)
Sp
В (2.2) сумма берется по всем элементам, имеющим ненулевую площадь пересечения с "частицей-облаком". Электромагнитные силы Рт определяются в барицентрах элементов Cm-
Для осуществления процедуры локализации частиц в расчетной области вводится прямоугольная равномерная сетка (см. например [21]). Каждой ячейке этой сетки сопоставляется элемент вспомогательного массива данных, который содержит глобальные номера элементов неструктурированной сетки, принадлежащих данной ячейке. Процедура локализации частицы производится в два этапа. Сначала мы определяем с какими ячейками декартовой сетки пересекается данная частица. Затем в каждой из этих ячеек устанавливаем номера элементов неструктурированной сетки, с которыми пересекается частица, и вычисляем площадь пересечения частицы с каждым из таких элементов.
Для вычисления плотности тока определим "токовую сетку" на серединах сторон каждого треугольника (рис. 9, а). Заряд, переносимый частицей через какую-либо сторону элемента неструктурированной сетки в течение временного шага т, очевидно равен
SQ = PpSi, (2.3)
где Si — площадь пересечения "частицы-прямоугольника" с параллелограммом aba b , образованным векторами ab и a b' = ab — vt (рис. 9, б). Плотность тока через данную сторону элемента можно определить из условий j\\v, 5Q = (j,n)LT, где n — вектор внешней нормали, L = |ab|, т — временной шаг, а
6Q = ppSi = |j |(eV ,u)Lt,
j = pPSi = pSi\j\ (24)
171 Lt(ev ,n) lt(v, n)' 1 " ;
ppSi ->
j = t^vt v-(v, n)LT
Рассмотрим произвольный элемент неструктурированной сетки с номером m. Изменение заряда в этом элементе, вызванное движением отдельной макрочастицы с номером p, в течение временного шага т есть
3
$Qmp т ^ ^(jp,n%)Lm,i-
i=1
Здесь плотность тока jp определяется в соответствии с формулой (2.4). Взяв сумму по всем макрочастицам, получим, что изменение заряда в элементе с номером m, 5Qm = (pm1 — pnm)Sm определяется соотношением
(pm+1 —fm)Sm + ¿ ¿(jp, n)Lmi = 0. (2.5)
p i=1
Соотношение (2.5) представляет собой дискретный аналог уравнения непрерывности заряда.
Рассмотрим теперь результаты некоторых тестовых расчетов. Прежде всего был выполнен тест на восстановление равномерной плотности р0 = 1. Для этого в расчетной области размещалось определенное количество макрочастиц. Макрочастицы представляли собой квадраты определенного размера. Заряд всех макрочастиц предполагался одинаковым и равным Q = 1 /Ыр. Затем, в соответствии с формулой (2.1), определялась плотность в барицентрах неструктурированной сетки. Максимальные отклонения плотности от ее теоретического значения (в процентах) в зависимости от размеров и числа частиц представлены в таблице. Как показывают результаты численных экспериментов, если размер макрочастицы выбирался равным Норх = 1/\f~Np, отклонение плотности от теоретического значения не превышало величины "машинного нуля".
Число частиц в области Размер частицы Отклонение, %
2500 0.0199990 3.6 Е-3
2500 0.0180 7.11
2500 0.0150 16.7
9000 0.0333330 1.09 Е-3
9000 0.0320 4.23
9000 0.030 10.2
Для проверки дискретного аналога уравнения непрерывности мы рассмотрели движение пробной частицы на неструктурированной сетке, содержащей 130 узлов, 222 элемента (см. рис. 4, а) и 351 "токовый" узел. Частица представляла собой квадрат размером 0.1x0.1. В процессе движения данной макрочастицы вычислялась величина
6 = тах
Е(рт+1 - рпт^т+т]г в,я)ьг.
(2.6)
которая дает максимальное отклонение от дискретного аналога уравнения непрерывности (2.5). Вычислительные эксперименты показали, что абсолютное значение 6 не превышало величины "машинного нуля".
т
Рис. 9. К вычислению плотности тока от частиц конечного размера: а — элемент сетки токов: 1 , 2 , 3 — узлы "токовой" сетки, б — "токовый" параллелограмм.
Рис. 10. Движение "пробной частицы" конечного размера в центрально-симметричном поле: I — 130 узлов, 222 элемента, II — 701 узел, 1320 элементов. а — расчет, б — аналитическое решение.
Решение задачи Кеплера (для сеток, изображенных на рис. 4) представлено на рис. 10. Рассматривалось движение "пробной частицы" конечного размера в форме квадрата со стороной Л = 0.05. Для "пробной частицы" на грубой сетке (см. рис. 4, а) следует иметь в виду, что размер частицы сравним с размером типичного элемента сетки. В этом случае частица в процессе своего движения имеет ненулевую площадь пересечения (как правило) с одним элементом сетки, и сила, действующая на частицу, определяется по значению силы в барицентре этого элемента. Здесь обнаруживается аналогия с моделью "ближайшего узла" (МОР) для регулярной сетки. Поэтому имеет место весьма большое отличие от аналитического решения. При расчетах на более подробной сетке (см. рис. 4, б) интерполяция силы в местоположение частицы осуществляется из нескольких барицентров сетки (модель С1С для регулярной сетки). В этом случае соответствие между результатами численного и аналитического решений достаточно хорошее.
3. Заключение
В настоящей работе были рассмотрены численные алгоритмы, реализующие метод частиц на неструктурированной сетке. Предъявлены процедуры, обеспечивающие передвижение
частиц, их локализацию на неструктурированной сетке (двумерной или трехмерной), интерполяцию сил в их местоположение, а также раздачу зарядов в узлы (барицентры) неструктурированной сетки. Рассмотрены алгоритмы, основанные как на концепции точечных частиц, так и частиц конечного размера. Использование последней приводит к выполнению дискретного аналога уравнения непрерывности. Результаты тестовых расчетов свидетельствуют о работоспособности этих алгоритмов.
Список литературы
[1] БЕРЕЗИН Ю. А., Вшивков В. А. Метод частиц в динамике разреженной плазмы. Наука, Новосибирск, 1980.
[2] Хокни Р., Иствуд Дж. Численное моделирование плазмы методом частиц плазмы. Мир, М., 1987.
[3] БЭДСЕЛ Ч., ЛЕНГДОН А. Физика плазмы и численное моделирование плазмы методом частиц плазмы. Энергоатомиздат, М., 1989.
[4] Сигов Ю. С., ХОДЫРЕВ Ю. В. К теории дискретных моделей разреженной плазмы. Численные методы механики сплошной среды, 7, 1976, 109-117.
[5] MORSE R. L., NlELSON C. W. Numerical simulation of warm two beam plasma. Phys. Fluids, 12, 1969, 2418-2425.
[6] Bluhm H., Hoppe P., BACHMANN H. ET.AL. Progress in the development of a high power focusing B-applied extractortype ion diode for 1.5 TW pulse generator KALIF. Proc. 8th Int. Conf. on High-Power Particle Beams, Novosibirsk, 2-5 July 1990, World Scientific Singapore, 1990, 451-456.
[7] Matsumoto M., Kawata S. TRIPIC: triangular-mesh PIC code. J. Comp. Phys., 87, 1990, 488-493.
[8] Thopmson J.F., Warsi Z.U.A., Mastin C.W. Boundary-fitted coordinate systems for numerical solution of partial differential equations — a review. Ibid., 47, 1982, 1-108.
[9] JONES M. E. Electromagnetic PIC codes with body-fitted coordinates. Proc. 12th Int. Conf. on the Numerical Simulation of Plasmas, San Francisco, 20-24 September 1984, IM3, 1987.
[10] WESTERMANN T. Numerical modelling of the stationary Maxwell — Lorentz system in technical devices. Int. J. Num. Mod., Electronic Network, Devices and Fields, 7, 1994, 43-67.
[11] BORIS J. P. Relativistic plasma simulation — optimization of a hybrid code coordinates. Proc. 4th Int. Conf. on the Numerical Simulation of Plasmas, Washington, 20-24 September 1970, 3-6.
[12] SELDNER D., WESTERMANN T. Algorithms for interpolation and localization in irregular 2D meshes. J. Comp. Phys., 79, 1988, 1-11.
[13] Westermann T. Localization shemes in 2D boundary-fitted grids. Ibid., 101, 1992, 307313.
[14] Halter E., KrauSS M., C.-D. Munz ET. AL. A concept for numerical solution of the Maxwell — Vlasov system. Forshungszentrum Karlsruhe — Umwelt und Technik, FZKA 5564, 1995.
[15] Assous F., Degond P., Segre J. A particle-tracking method for 3D electromagnetic PIC codes on unstructured meshes. Comput. Phys. Commun., 72, 1992, 105-114.
[16] Sonnedrucker E., Ambrosiano J. J., Brandon S. T. A finite element formulation of the Darwin PIC model on unstructured grids. J. Comp. Phys., 121, 1995, 281-297.
[17] Кругляковл Л. В., Неледова А. В., Тишкин В.Ф., Филатов А. Ю. Неструктурированные адаптивные сетки для задач математической физики (обзор). Математическое моделирование, 10, 1998, 93-116.
[18] Assous F., Degond P., Heintze E. ET. AL. On a finite — element method for solving the three-dimensional Maxwell equations. J. Comp. Phys., 109, 1993, 222-237.
[19] Ambrosiano J. J., Brandon S.T., Lohner R., DeVore C.R. Electromagnetics via the Taylor — Galerkin finite element method on unstructured grids. Ibid., 110, 1994, 310-319.
[20] Shang J. S., Fithen R. M. A comparative study of characteristics-based algorithms for the Maxwell equations. Ibid., 125, 1996, 378-394.
[21] Lohner R., Ambrosiano J. A vectorized particle tracer for unstructured grids. Ibid., 91, 1990, 22-31.
[22] Lohner R. Robust, vectorized search algorithms for interpolation on unstructured grids. Ibid., 118, 1995, 380-387.
[23] Marder B. A method for incorporating Gauss' law into electromagnetic PIC codes. Ibid., 68, 1987, 48-51.
[24] Langdon A. B. On enforcing Gauss' law in electromagnetic particle-in-cell codes. Comp. Phys. Commun., 70, 1992, 447-450.
[25] Villasenor J., Buneman O. Rigorous charge conservation for local electromagnetic field solvers. Ibid., 69, 1992, 306-361.
[26] Березин Ю. А., Брейзман Б.Н., Вшивков В. А. Численное моделирование ин-жекции мощного электронного пучка в вакуумную камеру с сильным магнитным полем. Журн. прикл. мех. и техн. физики, №1, 1981, 3-9.
[27] Зенкевич О., Морган К. Конечные элементы и аппроксимация. Мир, М., 1986.
[28] Ландау Л. Д., Лифшиц е. М. Механика. Наука, М., 1988.
Поступила в редакцию 6 марта 1998 г., в переработанном виде 4 июня 1998 г.
