CC BY
135
Вычислительные технологии
Том 18, № 6, 2013
О выборе числа частиц в методе частиц-в-ячейках для моделирования задач физики плазмы*
Е.А. Месяц1 , А. В. Снытников1, К. В. Лотов2,3 1 Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия 2Институт ядерной физики им. Г.И. Будкера СО РАН, Новосибирск, Россия 3Новосибирский государственный университет, Россия e-mail: [email protected], [email protected], [email protected]
Исследована точность решения созданного параллельного комплекса программ, основанного на методе частиц-в-ячейках. Исследование проводилось на примере задачи эволюции тёплого электронного пучка малой плотности в плазме в трёх разных режимах неустойчивости. Проведён сравнительный анализ двух методов диагностики неустойчивости, из которых выбран оптимальный, позволяющий оценить уровень ошибки вычислений. Найдена зависимость точности полученных результатов от количества модельных частиц. Определено минимально необходимое количество модельных частиц для корректного воспроизведения инкремента неустойчивости.
Ключевые слова: метод частиц-в-ячейках, пучковая неустойчивость, уравнение Власова, уравнения Максвелла.
Введение
В настоящее время для моделирования нестационарных задач физики плазмы широко используется метод частиц-в-ячейках. Данный метод сочетает в себе ряд преимуществ, что ставит его в этих задачах на первое место. Метод физически нагляден, довольно прост в реализации, менее затратен по сравнению с конечно-разностными методами вычисления функции распределения частиц по скоростям в трёхмерном пространстве и методами частица-частица, когда рассчитывается непосредственное взаимодействие всех частиц друг с другом. Существует ряд пакетов, например, российские KARAT [1] и Mandor [2], построенных на этом методе и позволяющих проводить моделирование методом частиц в различных системах координат в одно-, двух- и трёхмерной постановках. С развитием вычислительной техники и появлением высокопроизводительных вычислительных систем появилась возможность решать всё более сложные задачи и моделировать всё более тонкие физические эффекты. В связи с этим вопрос точности метода частиц-в-ячейках вновь становится актуальным, так как наличие численных шумов всегда было недостатком данного метода. Источником нефизических шумов является дискретность модели: наличие сетки, на которой вычисляются плотности, токи и поля, разбиение среды на модельные частицы, дискретный шаг по времени. Для уменьшения
* Работа выполнена при поддержке РФФИ (гранты № 11-01-00249а и 12-07-00065а), Министерства образования и науки Российской Федерации (соглашение 14.B37.21.0784) и Интеграционного проекта СО РАН № 130.
численных шумов используются тихий старт, ядра более высоких порядков для частиц [3], фильтры Фурье для обрезания шумовых гармоник [4] и др. Разрабатываются также модификации алгоритма, например, так называемый 8/-метод [5, 6] или метод с периодическим восстановлением функции распределения частиц по скоростям [7]. Все эти методы, используемые по отдельности или в сочетании друг с другом, имеют как свои преимущества, так и свои области применимости. Вместе с тем они ведут к усложнению алгоритма. Самым простым и наиболее часто используемым методом уменьшения численных шумов является увеличение количества модельных частиц. Однако в силу ограниченности ресурсов ЭВМ это не всегда возможно. В работе [3] довольно подробно освещён вопрос точности метода частиц-в-ячейках, даны оценки и порядок аппроксимации для применяемых схем. В настоящей работе приводятся конкретные оценки необходимого числа частиц на примере задачи эволюции тёплого электронного пучка в плазме. Известно, что шумовая погрешность убывает медленно, как N-1/2, что ведёт к резкому увеличению времени расчётов, а также объёмов хранимых и передаваемых данных. Поэтому важно ответить на вопрос, сколько частиц достаточно для качественно и количественно удовлетворительного решения конкретной задачи. Ответ на него также позволит понять, насколько в том или ином случае необходимо использовать модификации алгоритма метода частиц-в-ячейках.
В рамках работы по совместным проектам с ИЯФ СО РАН авторами была создана трёхмерная численная модель взаимодействия пучка электронов с плазмой, основанная на методе частиц и ориентированная на исследование устойчивости и нагрева плазмы электронным пучком. Несмотря на то что задача о релаксации электронного пучка в плазме является классической проблемой физики плазмы [8, 9] и существует целый ряд теоретических моделей, описывающих различные режимы пучково-плазменного взаимодействия, исследования в этой области остаются актуальными [10, 11]. Максимально приближенная к условиям лабораторных экспериментов постановка задачи зачастую требует отказа от привычных для теории идеализаций, таких как слабое или сильное магнитное поле, гидродинамический или кинетический характер пучковой неустойчивости, приближение случайных фаз возбуждаемых в плазме турбулентных пульсаций. Кроме того, при длительной инжекции пучка эволюция пучково-плазменной системы может проходить через целую последовательность стадий, определяемых совершенно разными нелинейными процессами. Проводимая работа была направлена на изучение сценариев возбуждения плазменной турбулентности электронным пучком. При этом наибольший интерес представляют параметры пучка и плазмы, характерные для экспериментов по нагреву плазмы в открытой ловушке ГОЛ-3 (ИЯФ СО РАН) [12]. Для исследования влияния нелинейностей на поведение неустойчивости в условиях развитой турбулентности необходимо численное моделирование, которое, с одной стороны, способно на больших временах отслеживать эволюцию возбуждаемой пучком турбулентности, а с другой — позволяет обеспечить достаточно подробное описание кинетических эффектов, связанных с захватом пучка.
Основная цель настоящей работы — исследование точности метода частиц-в-ячейках в зависимости от количества частиц. Проводится моделирование возбуждения электронным пучком плазменной волны. Формулируются критерии оценки точности решения. На примере задачи о взаимодействии в разных режимах электронного пучка с плазмой даётся ответ на вопрос, сколько частиц достаточно не только для получения качественно правильных результатов, но и для совпадения значений физических величин с результатами теоретических оценок. В разделах 2, 3 приводятся физическая
модель и алгоритмы численного решения, в разделе 4 дано краткое описание параллельной реализации алгоритма, в разделе 5 представлены два подхода к вычислению инкремента неустойчивости, по которым далее оценивается точность решения, в разделах 6 и 7 подробно описаны результаты расчётов для трёх режимов развития неустойчивости при разном количестве частиц, оцениваются точность полученных результатов и качественное поведение электронов пучка на фазовой плоскости. Показано, что вычисление инкремента по амплитуде главной неустойчивой моды предпочтительнее, поскольку данный способ свободен от эффекта саморазогрева модельной плазмы и даёт большую точность при меньшем числе модельных частиц.
1. Основные уравнения
Модель высокотемпературной бесстолкновительной плазмы представляется кинетическим уравнением Власова и системой уравнений Максвелла [13], которые в безразмерной форме имеют вид
д/\е д/^е /-,-, г т->тчд/^е п ПЧ
~дГ + ""дГ + ^(Е + "д"" = 0' (1)
д Е „ .
ж = ™»в -(2)
дВ
ж =- го( Е' (3)
а1у Е = р, (4)
в = о. (5)
Здесь индексами г и е помечены величины, относящиеся соответственно к ионам и электронам; qi = ше/ш^ ^е = -1; /г,е(£, г, р) — функция распределения частиц; ш^е, р1е, г1,е — масса, импульс, положение иона или электрона; Е, В — напряжённости электрического и магнитного полей. Для перехода к безразмерному виду используются следующие базовые величины: скорость света с = 3-1010 см/с; масса электрона ше = 9.1 • 10-28 г; плотность плазмы п0 = 1014 см-3; время £ = '}, где плазменная электронная частота шре = 56 • 1011 с-1.
В начальный момент времени в трёхмерной области решения, имеющей форму прямоугольного параллелепипеда
х е [0,Ь], у, г е [0,Ь±],
находятся плазма, состоящая из электронов и ионов водорода, и пучок электронов. Задаются плотности пучка пь и электронов плазмы пе = 1 — пь, температура пучка ТЬ и электронов плазмы Те; температура ионов считается нулевой Т = 0. Начальное распределение частиц по скоростям максвелловское с плотностью распределения
, 1 ( (V — ^о)2 / (v) = -т-ехр —-
А^л/2^ V 2Av2
где Av — разброс частиц по скоростям, v0 — средняя скорость пучка. Средняя скорость ионов и электронов фона нулевая. Все частицы (электроны и ионы плазмы, электроны пучка) распределены по области равномерно, начальная средняя скорость пучка
направлена по оси х и равна го = 0.2. Граничные условия для всех функций по всем направлениям периодические.
В расчётах нас интересует развитие отдельно взятой неустойчивой моды, поэтому длина области в направлении х выбрана равной одной длине исследуемой плазменной волны: Ьу = Ьх = = 4Нх. Так как граничные условия для всех функций по всем направлениям периодические, то из непрерывного спектра волновых чисел к вырезаются дискретные значения к3 = 2пз/Ь3 (в = х,у,г), из которых только основная мода с кх = 2п/Ь, ку = кх = 0 попадает в область большого инкремента. Эта волна и характеризующий её инкремент и будет рассматриваться в расчётах. В зависимости от разброса электронов пучка по скоростям различают несколько режимов неустойчивости: гидродинамический (кДг ^ 7), переходный (кДг ~ 7) и кинетический (кДг ^ 7). Здесь Дг — разброс электронов пучка по скоростям, 7 — инкремент развивающейся неустойчивости. Требуется найти распределение ионов и электронов по энергиям и инкременты плазменных неустойчивостей для физических параметров, соответствующих этим трём режимам.
2. Методы и алгоритмы решения
Решение уравнения Власова проводится методом частиц-в-ячейках [3]. Плазма представляется набором модельных частиц, траекториями движения которых являются характеристики уравнения Власова
(Рм
(И
дг,е(Е + [V, В]),
(Г1,6
(г
VI,е, Р1,е
VI,.
1 - V?
Для решения системы уравнений (6) используется схема с перешагиванием
т+1/2 т—1/2
Ра - Ра
%е[ Ет +
т+1/2 . т—1/2
Vа + Vа
Вт
„т+1 _ „т
Га_„а = vm+l/2
Т
Здесь а — номер частицы.
Плотность заряда вычисляется по положениям частиц га = (ха,уа, га) с использованием Р1С-ядра
Рт+1/2 ,.7+1/2,1+1/2 = X/ Ча*
г+1/2 ха )П(Уу+1/2 - Уа)Щ^+1/2 - Zа)),
где да — заряд частицы а, а Я(х) имеет вид
адН К1 "1Г1 ■ 'х'-Н
0,
' х ' > Н.
2
Для решения системы уравнений Максвелла используется метод Лэнгдона — Лазинского, описанный в работе [14], в котором поля определяются из разностных аналогов законов Фарадея и Ампера:
B m+1/2 _ Bm-1/2
-= - roth Em
т
pm+1 _ pm
E-— = roth Bm+1/2 - jm+1/2.
т
Значения компонент напряжённостей электрического и магнитного полей вычисляются в узлах сеток, смещённых относительно друг друга по времени и пространству:
-r>m+1/2_/R1 r2 r3 \m+1/2
B — VBi+1/2,j,l, Bi,j+1/2,1 ,Bi,j,1+1/2) ,
pm _ / тр1 тр2 тр3 \m
E — VEi,j+1/2,Z+1/2, Ei+1/2,j,1+1/2, Ei+1/2,j+1/2,iy , m+1/2 — ( -1 • 2 • 3 )m+1
•m+1/2 — ( -1 • 2 • 3 )
j — \Ji,j+1/2,l+1/2,Ji+1/2,j,l+1/2,Ji+1/2,j+1/2,l)
m — m
P — pi+1/2,j+1/2,l+1/2
Разностные операторы roth и divh на такой сетке имеют следующий вид (в индексах опустим 1/2):
roth H
u3 _ u3 тт2 _ тт2
Hi,j,l Hi,j-1,l Hi,j,l Hi,j,l-1
hy hz
тт1 _ тт1 тт3 _ тт3 tt2 _ тт2 тт1 _ тт1
Hi,j,l Hi,j,l-1 Hi,j,l Hi-1,j,l Hi,j,l Hi-1,j,l Hi,j,l Hi,j-1,l
hz hx hx hy
тт1 _ 7т1 7Т2 _ 71Г2 гтЗ _ гтЗ
Hi+1,j,l Hi,j-1,l , Hi,j+1,l Hi,j,l , Hi,j,l+1 Hi,j,l
т XT i+1,/,i i,'/- 1,l , i,'/ + 1,t i,/,l ,
divh H — —+ 7--—- +---— +
кХ ку к2:
Данная схема имеет порядок аппроксимации 0(т2 + к2). Если в начальный момент времени divh В0 = 0, а токи вычисляются таким образом, что выполнено уравнение неразрывности
0т+1 _ пт
Р-^ + divh jm+1/2 — 0,
т
то разностные аналоги уравнений (4) и (5) выполняются автоматически. Чтобы удовлетворить этим условиям, токи вычисляются по схеме, приведённой в работах [3, 15], и В0 = 0.
Большой объём данных, связанных с трёхмерностью задачи, требует использования распределённых вычислительных систем. Распараллеливание выполнено методом декомпозиции расчётной области по направлению у, перпендикулярному движению электронного пучка х. Применена смешанная эйлерово-лагранжева декомпозиция (рис. 1). Сетка, на которой решаются уравнения Максвелла, разделена на одинаковые подобласти по одной из координат. С каждой подобластью связана группа процессоров. Далее модельные частицы каждой из подобластей разделяются между процессорами связанной с этой подобластью группы равномерно, вне зависимости от координат. Каждый из процессоров группы решает уравнения Максвелла во всей подобласти. Затем решаются уравнения движения модельных частиц. После этого происходит суммирование значений тока по всей области. Один из процессоров группы производит обмен граничными
Рис. 1. Декомпозиция по области и частицам. Разбиение области — декомпозиция эйлерова этапа метода частиц-в-ячейках. Разбиение частиц — декомпозиция лагранжева этапа. Каждый тип маркеров обозначает частицы, принадлежащие одному процессору в данной группе
значениями тока и полей с соседними подобластями и далее рассылает полученные граничные значения всем процессорам своей группы. Подробности параллельной реализации описаны в работе [16].
3. Вычисление инкремента неустойчивости
Для оценки точности получаемого решения проводится сравнение инкремента развивающейся неустойчивости с его аналитическими оценками для той же задачи в одномерной постановке [17]. Известно, что на начальной стадии развития неустойчивости энергия электрического поля Ш нарастает по экспоненте: Ш х е27^. Поэтому инкремент неустойчивости можно вычислять как производную от логарифма энергии электрического поля во всей счётной области:
1 д
и = , дг )• (9)
Расчёты, однако, показали, что данный способ вычисления инкремента неустойчивости не подходит для вычислений в кинетическом режиме. Сложность перехода от гидродинамического режима к кинетическому состоит в следующем. Когда разброс по скоростям электронов пучка увеличивается, а плотность пучка становится меньше, во взаимодействие с возбуждаемой волной вступает все меньшая часть электронов пучка, и время насыщения неустойчивости увеличивается. При этом энергия электрического поля волны уменьшается и при недостаточном количестве частиц становится по величине сравнима с численными шумами. Такие расчёты требуют большого количества частиц, что не всегда возможно обеспечить. Поэтому важно, во-первых, подобрать менее чувствительный к численным шумам способ вычисления инкремента неустойчивости, а во-вторых, ответить на вопрос, сколько частиц в ячейке можно считать достаточным для вычислений.
В качестве такой диагностики был выбран расчёт инкремента неустойчивости по амплитуде отдельной неустойчивой моды Е0:
д
72 = дг 1п(Ео). (10)
Формулы для вычисления инкремента имеют следующий вид:
1 д ( 1 Пх Пу п* \ * = 5дг|п(Е Е• (11)
Пх Пу Пг / \ \ 2
-=IчпПлтИ £ ^.мЧ^^+
¿,¿7—1
I ПХ ПуV V 7 — 1 \ LХ Ьу Ьг
♦О^^ц^(12)
Здесь щ — число узлов сетки в направлении 5 (в = х, у, £), к = Ь/щ.
4. Результаты расчётов
На графиках (см. ниже) представлено сравнение результатов расчётов инкремента неустойчивости по энергии электрического поля (11) и по амплитуде главной моды (3) с решением, описанным в [17] (для одномерной задачи). Расчёты проводились при параметрах, соответствующих трём режимам развития пучковой неустойчивости. Счётные параметры, общие для всех расчётов, были следующими: сетка по пространству 100 х 4 х 4, шаг по времени т = 0.001. Обозначим через 1р число частиц в одной ячейке сетки и рассмотрим вопрос точности получаемого решения в разных режимах в зависимости от 1р.
4.1. Гидродинамический режим (kДv < 7)
Основные параметры: отношение плотности пучка к плотности плазмы щ/п0 = 2 • 10-3, разброс электронов пучка по скоростям = 0.035, длина области Ь = 1.2566.
При заданных параметрах величина инкремента неустойчивости должна быть 7 = 0.0722 [17].
На рис. 2 приведены темп роста энергии электрического поля и значение инкремента 71, вычисленное по формуле (11). На рисунках 2, а —6, а прямой линией показан
а б
Рис. 2. Временная зависимость логарифма энергии электрического поля (а) и производной логарифма энергии электрического поля (инкремент 71) (б) при разном числе частиц в ячейке. Гидродинамический режим
темп роста волны с аналитическим инкрементом, положение этой линии снесено по координате £, максимумы графиков совмещены по времени. На рис. 2, б — 6, б прямые линии — точное значение инкремента. На рис. 2 видна одна из особенностей моделирования холодной плазмы методом частиц — саморазогрев модельной плазмы, который выражается в том, что при малых температурах плазмы в самом начале вычислений энергия возрастает до определённой величины [18], зависящей от счётных параметров. Как видно из рисунка, уровень саморазогрева уменьшается с увеличением числа частиц. Когда амплитуда развивающейся неустойчивой волны превышает уровень шумов, связанных с саморазогревом, на графиках наблюдается экспоненциальная стадия роста волны (см. кривые на рис. 2, а). Далее решение выходит из линейного режима, и темп роста энергии спадает. Наступает стадия насыщения неустойчивости. В работе рассматривается только линейная стадия развития неустойчивости, на которой и вычисляется значение инкремента. На графиках, представленных на рис. 2, хорошо видна сходимость по количеству частиц к требуемому значению в линейной стадии.
Второй способ вычисления инкремента (по амплитуде главной волны электрического поля, формула (10)) представлен на рис. 3. Как видно из графиков, здесь уже не играют роли шумы от саморазогрева модельной плазмы, что позволяет практически
а б
Рис. 3. Временная зависимость логарифма амплитуды электрического поля (а) и производной логарифма амплитуды (инкремент 72) (б) при разном числе частиц в ячейке. Гидродинамический режим
Таблица 1. Инкремент и ошибка (%) в зависимости от количества частиц в ячейке. Гидродинамический режим
1р 71 Ошибка, % 72 Ошибка, %
50 0.058 20 0.066 9
100 0.062 14 0.068 6
500 0.067 7 0.069 4
1000 0.068 6 0.07 3
5000 0.07 3 0.071 2
10000 0.07 3 0.071 2
сразу наблюдать экспоненциальный рост. Более того, из представленных графиков следует, что для той же точности вычисления инкремента можно использовать меньшее количество частиц.
В табл. 1 приведены значения инкремента, вычисленные по энергии и по амплитуде, а также характерная величина ошибки (в %). В гидродинамическом случае ошибка по энергии примерно в 2 раза больше, чем по амплитуде. Из приведенных графиков, а также из таблицы видно, что в гидродинамическом режиме достаточным можно считать число частиц в ячейке, равное 50-100, если использовать вычисление инкремента по амплитуде. Ошибка в этом случае не превышает 10%.
4.2. Переходный режим (kДv ~ 7)
Основные параметры: отношение плотности пучка к плотности плазмы пь/п0 = 2 • 10-3, разброс электронов пучка по скоростям Аю/ю0 = 0.14, длина области Ь = 1.1424. При этих параметрах величина инкремента неустойчивости должна быть 7 = 0.0232, т. е. по сравнению с таковой гидродинамического режима уменьшилась в три раза. Таким образом, энергия в рассматриваемом случае нарастает медленнее, насыщение неустойчивости происходит позже по времени. Как и при гидродинамическом режиме, на графике логарифма энергии электрического поля (рис. 4) виден саморазогрев модельной плазмы. На графике амплитуды основной моды (рис. 5) этого эффекта нет. Из приведённых рисунков следует, что в переходном режиме по сравнению с гидродинамическим необходимо использовать гораздо большее число частиц в ячейке. При увеличении числа частиц значение инкремента приближается к точному. При недостаточном числе частиц (1р = 100) график инкремента 72 (рис. 5, б) не имеет области, близкой к горизонтальной прямой, которая чётко наблюдалась в гидродинамическом случае и которая появляется в переходном режиме при большем числе частиц. При 1р = 500 область
а
б
0,04
• 1р-Л 00 ♦ //7 = 500 1р= 1000 1р = 5000 -/7? =10000 -----Ут = 0.0232
0,03-
0,01
0,02
0
-12
У-12.9
-0,01
0
50 100 150 200 250
50 100 150 200 250
г
Рис. 4. Временная зависимость логарифма энергии электрического поля (а) и производной логарифма амплитуды (инкремент 71) (б) при разном числе частиц в ячейке. Переходный режим
а б
Рис. 5. Временная зависимость логарифма амплитуды электрического поля (а) и производной логарифма амплитуды (инкремент 72) (б) при разном числе частиц в ячейке. Переходный режим
Таблица 2. Инкремент и ошибка (%) в зависимости от количества частиц в ячейке. Переходный режим
1р 71 Ошибка (%) 72 Ошибка (%)
500 0.017 27 0.02 14
1000 0.019 18 0.02 14
5000 0.02 14 0.021 9
10000 0.021 9 0.022 5
экспоненциального роста видна уже лучше. Далее с ростом числа частиц кривая 1п(Е0)/2 в интервале [0,150] приближается к прямой линии, а 72(£) выходит на уровень, примерно соответствующий точному значению инкремента. Таким образом, в данном случае для того, чтобы вычислить инкремент, необходимо использовать как минимум 500 частиц в ячейке.
По ресурсам вычисления становятся более затратными, считать приходится на большие времена с большим числом частиц, чем в гидродинамическом случае. Из табл. 2 видно, что ошибка по энергии также в 1.5-2 раза больше, чем по амплитуде. Чтобы ошибка вычисления инкремента по амплитуде была в пределах 10%, число частиц в одной ячейке должно быть 1р ~ 5000. Для достижения той же точности при вычислении инкремента по энергии электрического поля потребуется 10 000 частиц в одной ячейке.
4.3. Кинетический режим (кАу ^ 7)
Основные параметры следующие: отношение плотности пучка к плотности плазмы Щ/п0 = 2 • 10-4, разброс электронов пучка по скоростям Ау/у0 = 0.14, длина области Ь = 1.1424. Аналитическое значение инкремента неустойчивости 7 = 0.0027.
Для вычисления инкремента в кинетическом режиме требуется гораздо больше частиц, чем в переходном. Причина этого станет понятной далее при рассмотрении картины поведения электронов пучка на фазовой плоскости. Величина инкремента по срав-
а б
Рис. 6. Временная зависимость логарифма амплитуды электрического поля (а) и производной логарифма амплитуды (инкремент 72) (б) при разном числе частиц в ячейке. Кинетический режим
нению с таковой переходного режима уменьшилась почти в 9 раз, время насыщения неустойчивости увеличилось примерно в 6 раз, из-за численных шумов вычисление инкремента по энергии электрического поля стало невозможным. Поэтому здесь показан только график амплитуды электрического поля (рис. 6). Несмотря на перечисленные сложности, вычисление инкремента по амплитуде даёт удовлетворительный результат и показывает, что развивается именно та неустойчивая волна, которая соответствует выбранным начальным и граничным значениям функций и значениям параметров.
Как и в переходном режиме, при недостаточном числе частиц график инкремента 72 при 1р = 1000 не имеет области, близкой к горизонтальной прямой, и соответственно выбрать точку для вычисления инкремента невозможно. С увеличением числа частиц область экспоненциального роста становится шире. Для вычисления значения инкремента в кинетическом режиме с ошибкой в 25 % требуется около 5000 частиц в ячейке. Отметим, что в отличие от гидродинамического режима графики в переходном и кинетическом режимах приведены со сглаживанием.
5. Фазовые плоскости
Помимо инкремента, признаком развития неустойчивости является характерное поведение электронов пучка на фазовой плоскости (х, ух). Для визуализации этого поведения отслеживалась эволюция выбранных частиц, начальное положение которых представлено на рис. 7 (верхний ряд). В нижнем ряду показано положение электронов пучка в моменты времени, соответствующие насыщению неустойчивости. Для кинетического режима (рис. 7, в) приведена только узкая область по скорости 0.17 < Ух < 0.2, но ширина начального распределения частиц здесь такая же, как в переходном режиме (рис. 7, б, £ = 0). Видно, во-первых, что в гидродинамическом режиме возбуждаемой волной захватывается весь пучок, а в переходном и кинетическом режимах — всё более узкая его часть, остальные же электроны во взаимодействие с волной не вступают.
Этим, в частности, и обусловливается то, что в кинетическом режиме требуется гораздо большее число частиц. Несмотря на малое (с точки зрения точности вычисления инкремента в кинетическом режиме) количество частиц, наблюдается захват именно тех частиц, скорость которых близка к расчётной фазовой скорости волны V = 0.9^0 = 0.18. Это соответствует аналитическим оценкам для одномерной задачи [19].
Обозначим через 5vx ширину по vx области захваченных частиц. Для переходного
7
и кинетического режимов Svx = С0 —, где значение константы С0 для всех расчётов
к
в переходном и кинетическом режимах одинаково и равно С0 ~ 14, что соответствует теоретическому значению 13.6 [19]. Следовательно, можно заключить, что построенная модель позволяет не только качественно, но и количественно правильно моделировать пучковую неустойчивость.
Рис. 7. Электроны пучка на фазовой плоскости (x, vx) для а — гидродинамического (Av = 0.007, lp = 100), б — переходного (Av = 0.02, lp = 500) и в — кинетического (Av = 0.02, lp = 5000) режимов
Рис. 8. Электроны пучка на фазовой плоскости (х, ух) для а — переходного (А^ = 0.02, 1р = 100) и б — кинетического (А^ = 0.02, 1р = 1000) режимов
Для сравнения на рис. 8 приведены фазовые портреты пучка при £ = 480 и 2100 для переходного и кинетического режимов при недостаточном числе частиц (1р = 100 и 1000 соответственно), когда вычислить инкремент невозможно. При недостаточном числе частиц область захваченных частиц формируется гораздо раньше и разрушается быстрее. Видно, что в этом случае картина не столь чёткая, как на рис. 7. Помимо основной волны здесь возникает ряд волн сравнимой с ней амплитуды. Эти волны не просто маскируют решение, а разрушают его. Такие волны, имеющие чисто вычислительную природу, называются в методе частиц-в-ячейках численными шумами. Опасность их заключается в том, что они не просто накладываются сверху как помеха, которую легко отделить от решения, но взаимодействуют друг с другом и с решением. В таком случае назвать полученные картины решением исходной задачи нельзя. Можно заключить, что невозможность вычислить инкремент по формуле (3) вызвана ошибками в моделировании явления, а не малой чувствительностью диагностики. Следовательно, совпадение рассчитанного по формуле (3) значения инкремента с теоретически предсказанным его значением может служить критерием соответствия результатов моделирования физическому явлению, а также количественной мерой его точности.
Таким образом, по начальному количеству частиц можно оценить, эффекты какой величины при моделировании будут воспроизводиться правильно. В среднем при рассмотренных параметрах в трёхмерной постановке для того, чтобы иметь возможность вычислять в переходном и кинетическом режимах инкремент с точностью 25 %, необходимо иметь около 200 000 модельных частиц пучка в интервале скоростей 8ух, в котором происходит взаимодействие частиц с неустойчивой волной. В случае, когда характеристики ЭВМ не позволяют использовать такое количество частиц, необходимо искать алгоритмы решения, позволяющие для получения более точных результатов понизить уровень численных шумов при моделировании.
Заключение
На основе метода частиц-в-ячейках создан комплекс программ для моделирования в трёхмерной постановке взаимодействия пучка электронов с плазмой. При различных счётных параметрах исследована точность получаемого решения. С помощью данного комплекса программ можно моделировать эволюцию и неустойчивость тёплого электронного пучка в плазме в трёх разных режимах неустойчивости. Проведён сравнительный анализ двух методов диагностики неустойчивости. Показано, что вычисление инкремента по амплитуде электрического поля отдельной неустойчивой моды в сравнении с вычислением по энергии электрического поля является более чувствительной диагностикой. Предложенный метод позволяет надёжно идентифицировать экспоненциальный рост плазменной волны на начальной стадии развития неустойчивости, когда полная энергия волны ещё мала по сравнению с энергией тепловых шумов. Показано: чтобы помимо качественного воспроизведения картины развития неустойчивости было ещё и совпадение значений физических величин с результатами теоретических оценок, необходимо иметь не менее 200 000 модельных частиц пучка в интервале скоростей, в котором происходит взаимодействие частиц с неустойчивой волной.
Авторы выражают благодарность Г.И. Дудниковой и В.А. Вшивкову за полезные обсуждения представленного в статье материала.
Список литературы
[1] Tarakanov V.P. User's Manual for Code KARAT. Springfield, VA: Berkley Res., 1992.
[2] http://mandor.ilc.edu.ru/mandor3
[3] Григорьев Ю.Н., Вшивков В.А., Федорук М.П. Численное моделирование методами частиц-в-ячейках. Новосибирск: Изд-во СО РАН, 2004.
[4] Sydora R.D. Low-noise electromagnetic and relativistic particle-in-cell plasma simulation models // J. of Comput. and Appl. Math. 1999. Vol. 109. P. 243-259.
[5] Parker S.E., Lee W.W. A fully nonlinear characteristic method for gyrokinetic simulation // Phys. of Fluids. 1993. Vol. B5. P. 77-86.
[6] Lee W.W., Jenkins T.G., Ethier S. A generalized weight-based particle-in-cell simulation scheme // Comput. Phys. Communicat. 2011. Vol. 182. P. 564-569.
[7] Wang B., Miller G.H., Collela P. A particle-in-cell method with adaptive phase-space remapping for kinetic plasmas // SIAM J. on Sci. Comput. 2011. Vol. 33, No. 6. P. 3509-3537.
[8] МихАйловский А.Б. Теория плазменных неустойчивостей. Т. 1. М.: Атомиздат, 1975.
[9] Цывович В.Н. Теория турбулентной плазмы. М.: Атомиздат, 1971.
[10] Timofeev I.V., Lotov K.V. Relaxation of a relativistic electron beam in plasma in the trapping regime // Phys. of Plasmas. 2006. Vol. 13, No. 6. P. 062312.
[11] Burdakov A.V., Erofeev V.I., Kotelnikov I.A. Explanation of turbulent suppression of electron conductivity in the GOL-3 facility at the stage of relativistic electron beam injection // Fusion Sci. and Technol. 2005. Vol. 47, No. 1T. P. 74-77.
[12] Koidan V.S., Arzhannikov A.V., Astrelin V.T. et al. Progress in multimirror trap GOL-3 // Fusion Eng. and Design. 2005. Vol. 47, No. 1T. P. 35-42.
[13] Кролл Н., Трайвелпис А. Основы физики плазмы. М.: Мир, 1975.
[14] Langdon A.B., Lasinski B.F. Electromagnetic and relativistic plasma simulation models // Methods in Comput. Phys. 1976. Vol. 16. P. 327-366.
[15] Вшивков В.А., Вшивков К.В., ДудниковА Г.И. Алгоритмы решения задачи взаимодействия лазерного импульса с плазмой // Вычисл. технологии. 2001. Т. 6, № 2. С. 47-63.
[16] Snytnikov A.V. Supercomputer simulation of plasma electron heat conductivity decrease due to relativistic electron beam relaxation // Procedia Comput. Sci. 2010. Vol. 1, iss. 1. P. 607-615.
[17] Лотов К.В., Тимофеев И.В. Переходный режим одномерной двухпотоковой неустойчивости // Вестник НГУ. Физика. 2008. Т. 3, вып. 1. С. 62-65.
[18] Хокни Р., Иствуд Дж. Численное моделирование методом частиц. М.: Мир, 1987.
[19] Лотов К.В., Терехов А.В., Тимофеев И.В. О насыщении двухпотоковой неустойчивости электронного пучка в плазме // Физика плазмы. 2009. Т. 35, № 6. С. 567-574.
Поступила в 'редакцию 29 мая 2013 г., с доработки —10 октября 2013 г.
