Численное моделирование интенсивных стационарных пучков заряженных частиц в диодах

Ма-ю-шан Владислав Витальевич

Численное моделирование интенсивных стационарных пучков заряженных частиц в диодах

Численное моделирование сегодня стало одним из наиболее эффективных средств для изучения плазмы и динамики потоков заряженных частиц в эмитирующих установках, таких как ускорители, электронные и ионные источники. В основе большинства используемых алгоритмов для проведения численных экспериментов в данной области лежат различные методы частиц, в которых пучки представляются как набор большого количества модельных частиц. Существенным недостатком методов данной группы является появление счетных нефизических шумов, приводящих к искажению решения. Причин возникновения шумов много и в настоящее время нет единого подхода к решению данной проблемы. В настоящей работе рассматривается решение самосогласованной задачи моделирования интенсивного потока заряженных частиц в плоском диоде итерационным методом с моделью трубок тока и методом частиц в ячейках. Изучается влияние методов подавления счетного шума, а также способа начальной дискретизации модельных частиц на решения, получаемые с помощью различных методов частиц.

Общественные науки в целом

Дипломы

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 587d36355f1be77c40d58958

UUID: e176fc68-74c0-493f-b055-195f2b4ce1b2

Язык: Русский

Опубликовано: больше 7 лет назад

Просмотры: 9

Ма-ю-шан Владислав Витальевич

Источник: Санкт-Петербургский государственный университет

Комментировать 0

Рецензировать 0

Скачать - 1238641 bytes

Поделиться работой

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ КАФЕДРА ТЕОРИИ СИСТЕМ УПРАВЛЕНИЯ ЭЛЕКТРОФИЗИЧЕСКОЙ АППАРАТУРОЙ Ма-ю-шан Владислав Витальевич Выпускная квалификационная работа бакалавра Численное моделирование интенсивных стационарных пучков заряженных частиц в диодах Направление 010400.62 Прикладная математика и информатика Научный руководитель, доктор физ.-мат. наук, профессор Овсянников Д. А. Санкт-Петербург 2016

Содержание Введение .................................................................................................................. 3 Цель .......................................................................................................................... 5 Глава 1. Модели динамики потоков заряженных частиц ............................ 6 1.1 Уравнение Власова ..................................................................................... 6 1.2 Методы частиц ............................................................................................ 7 1.3 Постановка задачи .................................................................................... 10 Глава 2. Численные алгоритмы модели "частица-сетка" ......................... 13 2.1 Метод частиц в ячейках............................................................................ 13 2.1.1 Эмиссия новых частиц ..................................................................... 13 2.1.2 Раздача заряда на сетку .................................................................... 14 2.1.3 Решение уравнения поля .................................................................. 15 2.1.4 Интерполяция сил ............................................................................. 17 2.1.5 Интегрирования уравнений движения частиц ............................... 17 2.2 Итерационный метод ................................................................................ 20 2.2.1 Схема раздачи заряда в итерационном методе .............................. 20 2.2.2 Подавление вклада погрешности в решение.................................. 21 2.3 Отражение электронов ............................................................................. 23 Глава 3. Сглаживание сеточных величин ..................................................... 24 3.1 Линейные сглаживающие фильтры ........................................................ 25 3.2 Частотная фильтрация .............................................................................. 26 3.3 Сглаживающий кубический сплайн........................................................ 27 Глава 4. Программная реализация численных алгоритмов ...................... 28 Глава 5. Результаты моделирования .............................................................. 29 5.1 Тихий старт. Сравнение методов ............................................................ 29 5.2 Случайный старт. Применение сглаживания ......................................... 32 5.2.1 Случай итерационного метода ........................................................ 32 5.2.2 Случай метода частиц в ячейках ..................................................... 35 Заключение ........................................................................................................... 37 Список литературы ............................................................................................ 38 2

Введение Бурное развитие вычислительной техники и появление суперкомпьютеров открыло новую эпоху в вычислительных науках. Численное моделирование стало наиболее эффективным средством для исследования различных физических явлений. Возможность замены полноразмерного эксперимента компьютерным во много раз сокращает дистанцию между научными гипотезами и их практическим воплощением. Несомненным преимуществом численных экспериментов является и простота их проведения: значительно проще проделать один, два или даже сто расчётов, чем получить удобную для использования формулу. Сегодня вычислительные методы успешно используются для исследования сложнейших процессов во многих областях современной физики: изучения эволюции спиральных структур в галактиках, атмосферных процессов, плазмы в термоядерных установках и ускорителях, турбулентности в жидкостях. Проведение численных экспериментов обеспечивает сокращение сроков исследования и уменьшение его стоимости. Программное обеспечение для численного моделирования непрерывно развивается, возникает необходимость в рассмотрении более сложных моделей, возрастают требования и к точности численных алгоритмов, появляются новые практические и теоретические проблемы. Наиболее широко используемым инструментом изучения плазмы и динамики потоков заряженных частиц в эмитирующих установках, таких как ускорители, электронные и ионные источники, является метод частиц в ячейках [10-12,16,17]. Идея метода заключается в том, что плазма представляется как набор большого количества взаимодействующих частиц. Расчетная область разбивается сеткой, в узлах которой хранится информация об электромагнитном поле, плотности тока и заряда. Метод частиц в ячейках 3

естественным образом моделирует процесс движения частиц, учитывая их собственное влияние друг на друга через сеточные величины. Однако для моделирования стационарного состояния, например, вакуумного диода в режиме насыщения тока, более быстрым и экономичным является итерационный метод с моделью трубок тока [1,4,5], основанный на итерационном повторении процесса расчета лишь для одного поколения частиц, что позволяет заметно сократить время решения задачи. Одним из главных недостатков методов частиц является возникновение счетных нефизических шумов, приводящих к искажению решения. Основной причиной их возникновения является дискретная природа самих методов (представление среды в виде конечного числа частиц). В настоящее время нет единого подхода к решению данной проблемы. В работах [2,11] описаны причины появления вычислительных шумов и методы борьбы с ними, предложены алгоритмы вычитания шумовой добавки и новая схема раздачи заряда в методе частиц в ячейках, рассмотрен вопрос возникновения "самосилы". В [9,15] рассмотрены такие техники подавления шума как схема δf, квазислучайное начальное распределение частиц, применение к сеточным величинам сглаживающих фильтров, основанных на преобразовании Фурье. В работах [1,4,12] предложено использование нижнего релаксационного преобразования плотности заряда (потенциала) для уменьшения вклада погрешности в итерационном методе с моделью трубок тока. В работе [3] рассмотрен способ повышения устойчивости решения самосогласованной задачи итерационным методом, разработан алгоритм изменения параметра релаксации. 4

Цель Целью данной работы является исследование методов решения стационарной самосогласованной задачи о плоском диоде методами частиц. Для достижения данной цели было решено: 1. Составить математическую модель динамики интенсивных стационарных потоков заряженных частиц в плоском диоде 2. Создать программную реализацию численных алгоритмов 3. Провести сравнение сходимости итерационного метода и метода частиц в ячейках на примере плоского диода в условиях случайного и тихого старта. 4. Исследовать влияние применения алгоритмов сглаживания к пространственной плотности заряда для подавления счетного шума на точность решения методов частиц В ходе решения задачи эмиссионная способность катода предполагается достаточно высокой. 5

Глава 1. Модели динамики потоков заряженных частиц 1.1 Уравнение Власова Для описания динамики потока заряженных частиц или плазмы с учетом дальнодействующих кулоновских сил Власовым было предложено использовать кинетическое уравнение с самосогласованным полем [1,2] 𝜕𝑓𝛼 𝜕𝑓𝛼 𝐅𝐚 𝜕𝑓𝛼 +𝐯∙ + =0 𝜕𝑡 𝜕𝐫 𝑚𝑎 𝜕𝑡 (1.1) Здесь 𝑓𝛼 (𝑡, 𝐫, 𝐯) – функция распределения частиц по скоростям, 𝛼 – сорт частиц, 𝐫 – координата, 𝐯 – скорость, 𝐅𝐚 = 𝑞𝛼 (𝐄 + [𝐯𝐁]) – сила Лоренца, 𝐄 – напряженность электрического поля, 𝐁 – индукция магнитного поля, 𝐯 – скорость. Компоненты электромагнитного поля, определяются из системы уравнений Максвелла div 𝐄 = 𝜌 𝜀 div 𝐁 = 0 𝜕𝐁 𝜕𝑡 1 𝜕𝐄 rot 𝐁 = − 2 + 𝜇𝐣 𝑐 𝜕𝑡 rot 𝐄 = − (1.2) где 𝜌 – плотность пространственного заряда, 𝜀 – диэлектрическая проницаемость среды, 𝑐 – скорость света, 𝜇 – магнитная проницаемость среды, 𝐣 – вектор плотности тока. Плотности зарядов и токов вычисляются как моменты функции распределения 6

𝜌(𝑡, 𝐫) = ∑ 𝑞𝛼 ∫ 𝑓𝛼 (𝑡, 𝐫, 𝐯) 𝑑𝑣 𝛼 (1.3) 𝐣(𝑡, 𝐫) = ∑ 𝑞𝛼 ∫ 𝐯 𝑓𝛼 (𝑡, 𝐫, 𝐯) 𝑑𝑣 𝛼 В случае безвихревого электрического поля система уравнений (1.2) сводится к уравнению для потенциала 𝜑 𝐄 = −∇𝜑 (1.4) Так как зачастую аналитическое исследование этой системы затруднено или невозможно, широкое распространение получили численные методы решения уравнения Власова. 1.2 Методы частиц Выделяют три основных подхода для численного решения уравнения Власова: эйлеров, лагранжев и смешанный, эйлеро-лагранжев. Динамика потоков частиц и плазмы описывается эйлеровым способом, если наблюдение за величинами, характеризующих потоки частиц ведется в фиксированных точках пространства, однако нет никакой информации о положении конкретной частицы и ее траектории. Лагранжев подход, напротив, заключается в описании движения отдельных частиц. Смешанные же алгоритмы состоят и двух этапов. На первом этапе вычисляется результат воздействия частиц на моделируемую среду и сохраняется на неподвижной эйлеровой сетке, на втором этапе выполняется очередной временной шаг уравнений динамики частиц, правая часть которых вычислена на первом этапе. Примером применения смешанного и лагранжевого подходов являются методы частиц [8]. Это группа вычислительных методов, в которой потоки заряженных частиц представляются как набор достаточно большого количества модельных частиц. Благодаря своей эффективности, методы частиц нашли свое 7

применение не только в задачах физики плазмы, а также в задачах динамики жидкости и газа, моделировании звездных скоплений и др. Обычно выделяют три основных типа моделей частиц [8]. Модель "частица-частица" является простейшим для понимания и реализации способом учета взаимодействия заряженных частиц. Поток частиц представляется в виде крупных частиц, каждая из которых имеет набор атрибутов (положение, заряд, масса, скорость). Метод основан на непрерывном интегрировании уравнений движения частиц, с учетом их непосредственного взаимодействия. Для расчета силы, действующей на частицу, необходимо рассчитывать влияние на нее всех остальных частиц. Этот метод наиболее трудоемкий, его сложность – O(𝑁 2 ), где N – количество частиц, поэтому схема нашла свое применение только в небольшом числе задач. Данную модель целесообразно применять для систем, где необходим учет дальнодействующих сил (например, при моделировании звездных скоплений). Временной шаг схемы "частица-частица" представляется следующим образом: 1. Вычисление сил непосредственного взаимодействия частиц 2. Проинтегрировать уравнения движения 3. Изменить значение счетчика времени В методах группы "частица-сетка" сила вычисляется значительно быстрее, чем при использовании схемы "частица-частица", но менее точно. В моделях данного типа пространство разбивается сеткой, в узлы которой раздаются характеристики (например, плотность заряда) ближайших частиц. Дифференциальные операторы заменяются конечно-разностными аппроксимациями. Силы, действующие на частицу, вычисляются посредством интерполяции из ближайших узлов сетки. Его сложность – 𝑂(𝑁 𝑙𝑜𝑔 𝑀), где N – количество частиц, 𝑀 – количество узлов сетки. 8

Итерация схемы "частица-сетка" состоит из пяти этапов: 1. Раздача плотности заряда в узлы сетки 2. Решение уравнения Пуассона 3. Вычисление компонент напряженности, как градиент потенциала 4. Интерполяция сил, действующих на частицы. 5. Движение частиц под действием рассчитанных сил Гибридный метод "частица-частица – частица-сетка" сочетает в себе возможность учета близкодействия с точность метода "частица-частица" и быстродействие метода "частица-сетка". Действующая сила на частицу представляется как сумма быстроменяющейся близкодействующей части и медленноменяющейся дальнодействующей. Сложность схемы – 𝑂(𝑁 𝑙𝑜𝑔 𝑀 + 𝑁 𝑁𝑛 ), где N – количество частиц, N𝑛 – среднее количество соседей-частиц, M – количество узлов сетки. 9

1.3 Постановка задачи Рассмотрим математическую постановку задачи расчета динамики интенсивного стационарного потока заряженных частиц в плоском диоде. ̅ = Ω ∩ Γ двумерная расчетная область, где Γ = Γ1 ∪ Γ2 – Пусть Ω граница расчетной области. 𝐿 Г2 Катод U=0 Область эмиссии Г𝑒𝑚 𝐿𝑒𝑚 Анод U = U𝑎 Условия Неймана Г1 Г1 Г2 𝑑 Рисунок 1.1 Модель двумерного плоского диода Так как поток частиц является интенсивным, задача усложняется возникновением пространственного заряда 𝜌: катод работает в режиме ограничения тока, появляется необходимость учитывать кулоновские силы, создаваемые собственным пространственным зарядом пучка. Вычисление потенциала и напряженности электрического поля для стационарного потока заряженных частиц сводится к решению уравнения Пуассона 10

𝛥𝑈(𝐫) = − 𝜌(𝐫) , при 𝐫 𝜖 𝛺 𝜀0 𝐄(𝐫) = −grad(𝑈) 𝑈(𝐫) = 𝑔(𝐫), при 𝐫 𝜖 Г1 𝜕𝑈(𝐫) 𝜕𝐧 (1.5) = 0, при 𝐫 𝜖 Г2 В (1.5) 𝛥 – оператор Лапласа, 𝑈 – потенциал электрического поля, 𝐄 – напряженность электрического поля, 𝐫 – вектор фазовых координат, 𝜌(𝐫) – распределение плотности заряда, 𝑔(𝐫) – функция, описывающая потенциал на электродах, 𝜀0 – диэлектрическая проницаемость среды, 𝐧 – нормаль к границе Γ2 Предполагается, что задача имеет стационарное решение, тогда выполняется уравнения неразрывности плотности тока div 𝐣(𝐫) = 0 𝐣(𝐫) = 𝐣𝐞𝐦 (𝐫), при 𝐫 𝜖 Г𝑒𝑚 (1.6) Движение заряженной частицы в электростатическом поле описывается уравнением Ньютона-Лоренца и уравнением для радиус-вектора [17] 𝑑𝐩𝐢 = 𝑞𝑖 𝐄(𝐫𝐢 ) 𝑑𝑡 (1.7) 𝑑𝐫𝐢 𝐩𝐢 = 𝑑𝑡 𝑚0 В (1.7) 𝐩𝐢 – импульсы частиц, 𝑞𝑖 – заряды частиц, 𝐯𝐢 – скорости частиц, 𝑚0 – масса частиц, 𝐫𝐢 – положения частиц, 𝑖 – номер частицы. 11

Введем независимую переменную 𝜏 = 𝑐𝑡. В релятивистском виде, учитывающем зависимость массы частиц от скорости, уравнения (1.7) примут вид 𝑑 𝐩𝐢 𝑞𝑖 𝐄(𝑟𝑖 ) = 𝑑𝜏 𝑚0 𝑐 2 (1.8) 𝑑 𝐫𝐢 𝐩𝐢 = 𝑑𝜏 𝛾𝑖 Здесь 𝐩𝒊 = 𝛾𝑖 𝛃𝐢 - приведенный импульс частицы, 𝛃𝐢 = 𝐯𝐢 ⁄𝑐 – приведенная скорость частицы, 𝛾𝑖 = 1⁄√1 − 𝛃𝐢 𝟐 – приведенная энергия частицы. Плотность тока эмиссии в задаче считается постоянной и определяется из закона "трех-вторых" Чайлда [6] 3 𝐣𝐞𝐦 4 𝑒 𝑈2 = 𝜀0 √2 9 𝑚 𝑑2 (1.9) В (1.9) 𝑒 и 𝑚 – заряд и масса частицы, 𝑑 – расстояние между катодом и анодом, 𝑈 – напряжение на аноде. Требуется найти решение уравнения движения частиц (1.8) с учетом (1.5)-(1.6), определить распределение электростатического поля и траектории частиц. 12

Глава 2. Численные алгоритмы модели "частица-сетка" 2.1 Метод частиц в ячейках Метод частиц в ячейках относится к группе методов "частица-сетка" и объединяет в себе эйлеров и лагранжев подходы моделирования. Движение частиц в методе описывается подвижной лагранжевой сеткой, при этом плотность распределения пространственного заряда и значения электрических полей определены на неподвижной эйлеровой сетке. Временной шаг метода частиц в ячейках состоит из четырех основных этапов. На первом этапе происходит инжекция новых частиц в расчетную область. На втором этапе происходит расчет движения частиц в электромагнитных полях и раздача плотности заряда, создаваемой частицами, в узлы расчетной сетки. На третьем этапе решаются уравнения поля, определяется новая конфигурация электромагнитных полей. На последнем этапе происходит интерполяция сил, действующих на частицу. На следующем временном шаге все этапы повторяются. 2.1.1 Эмиссия новых частиц Эмиссия новых частиц в методе частиц в ячейках происходят на каждом временном шаге Δ𝑡. Каждой новой частице приписывается  заряд  импульс  координаты Заряд инжектируемой частицы определяется из плотности тока эмиссии. Координаты частицы определяются с помощью равномерного распределения по области эмиссии. Для уменьшения численных шумов возможно использование системы "тихий старт", требующая большего количества модельных частиц. 13

При наличии начального углового разброса y-компонента скорости определяется как величина, подчиняющаяся максвелловскому закону распределения 𝑚 𝑣𝑦 2 𝑚 𝑓(𝑣𝑦 ) = √ exp (− ) 2𝜋𝑇𝑒 2𝑇𝑒 (2.1) где 𝑚, 𝑒 – масса и заряд электрона, 𝑇 – температура эмиссии. 2.1.2 Раздача заряда на сетку В качестве схемы раздачи заряда на сетку выбрана модель cloud-in-cell (облаков в ячейке) [8,10]. Рассмотрим равномерную сетку с прямоугольными ячейками размером [ℎ𝑥 , ℎ𝑦 ]. Пусть одиночная частица с координатами (x𝑝 , y𝑝 ) расположена в ячейке (𝑖, 𝑗). Определим одномерную интерполяционную функцию для плотности заряда W(𝑥 − x𝑝 ) = { 1 − 0 |x − x𝑝 | , |x − x𝑝 | < ℎ𝑥 ℎ𝑥 в других случах (2.2) где (x𝑝 , 𝑦𝑝 ) координаты узла сетки в который распределяется заряд. В двумерном случае (2.2) пример вид W(𝑥 − 𝑥𝑝 , 𝑦 − 𝑦𝑝 ) = W(𝑥 − 𝑥𝑝 )W(𝑦 − 𝑦𝑝 ) (2.3) В двумерном случае заряд каждой частицы распределяется в четыре ближайших узла. В остальных узлах сетки плотность заряда будет равна нулю. Тогда схема распределения заряда для потока из 𝑁𝑝 заряженных частиц примет вид 𝑁𝑝 1 𝜌(𝑥, 𝑦) = ∑ 𝑞𝑖 W(x − x𝑖 , y − y𝑖 ) ℎ𝑥 ℎ𝑦 𝑖=1 где x𝑖 и 𝑞𝑖 – положение и заряд 𝑖-й частицы соответственно. 14 (2.4)

2.1.3 Решение уравнения поля Рассмотрим решение двумерного уравнения Пуассона методом конечных разностей [13] в прямоугольной расчетной области 𝜕2𝑈 𝜕2𝑈 𝜌 + = − 𝜕𝑥 2 𝜕𝑦 2 𝜀0 (2.5) с граничными условиями 𝑈(𝑥, 𝑦) = 𝑔(𝑥, 𝑦), при 𝑥, 𝑦 𝜖 Г1 (2.6) 𝜕𝑈(𝑥, 𝑦) 𝜕𝐧 = 0, при 𝑥, 𝑦 𝜖 Г2 где Γ = Γ1 ∪ Γ2 – граница расчетной области Разобьём расчетную область равномерной сеткой с размером ячейки [ℎ𝑥 , ℎ𝑦 ]. Проведя дискретизацию уравнения (2.5) для внутренних точек сетки с использованием пятиточечного шаблона получим 𝑈(𝑖 + 1, 𝑗) − 2 𝑈(𝑖, 𝑗) + 𝑈(𝑖 − 1, 𝑗) + ℎ𝑥 (2.7) + 𝑈(𝑖, 𝑗 + 1) − 2 𝑈(𝑖, 𝑗) + 𝑈(𝑖, 𝑗 − 1) 𝜌(𝑖, 𝑗) =− ℎ𝑦 𝜀0 Проведем теперь дискретизацию граничных условий 𝑈(𝜉, 𝜂) = 𝑔(𝑥, 𝑦) (2.8) 𝑈(𝜉, 𝜂) − 𝑈(𝑖, 𝑗) =0 ℎ где (ξ, η) номер ячейки, принадлежащей соответствующей границе, (i, j) – номер первой внутренней ячейки. Таким образом, в результате дискретизации получим разреженную СЛАУ размерности [𝑛 ∙ 𝑚 × 𝑛 ∙ 𝑚]. Для решения системы используется 15

стабилизированный метод бисопряженных градиентов с ILU-предобуславливанием. Напряженность электрического поля рассчитывается дифференцированием потенциала электрического поля. Пусть потенциал 𝑈 задан на прямоугольной сетке c размером ячейки [ℎ𝑥 , ℎ𝑦 ], а 𝑈(𝑖, 𝑗) – значение потенциала в (i,j) узле сетки. По определению 𝐄 = −grad(𝑈), тогда в двумерном случае 𝐄(𝐫) = −𝐱 𝜕𝑈(𝑥, 𝑦) 𝜕𝑈(𝑥, 𝑦) −𝐲 𝜕𝑥 𝜕𝑦 (2.9) Обозначим 𝜕𝑈(𝑥, 𝑦) 𝜕𝑥 𝜕𝑈(𝑥, 𝑦) 𝐸𝑦 (𝑥, 𝑦) = − 𝜕𝑦 𝐸𝑥 (𝑥, 𝑦) = − (2.10) Где 𝐸𝑥 (𝑥, 𝑦) – 𝑥-компонента напряженности, 𝐸𝑦 (𝑥, 𝑦) – 𝑦-компонента соответственно. Внутри расчетной области будем использовать центральную разностную схему дифференцирования 𝐸𝑥 (𝑖, 𝑗) = − 𝑈(𝑖 + 1, 𝑗) − 𝑈(𝑖 − 1, 𝑗) 2ℎ𝑥 𝑈(𝑖, 𝑗 + 1) − 𝑈(𝑖, 𝑗 − 1) 𝐸𝑦 (𝑖, 𝑗) = − 2ℎ𝑦 (2.11) на левой границе – левую разностную схему второго порядка 𝐸𝑥 (𝑖, 𝑗) = − −3𝑈(𝑖, 𝑗) + 4𝑈(𝑖 + 1, 𝑗) − 𝑈(𝑖 + 2, 𝑗) 2ℎ𝑥 −3𝑈(𝑖, 𝑗) + 4𝑈(𝑖, 𝑗 + 1) − 𝑈(𝑖, 𝑗 + 2) 𝐸𝑦 (𝑖, 𝑗) = − 2ℎ𝑦 на правой – правую разностную схему соответственно. 16 (2.12)

2.1.4 Интерполяция сил В качестве схемы интерполяции сил в местоположение частицы выбрана модель облаков в ячейке. Использование одинаковой схемы при раздаче заряда на сетку и интерполяции сил продиктовано уменьшением эффекта "самосилы" [2,8], при котором частица действует сама на себя через сетку. Рассмотрим ранее определенную интерполяционную функцию (2.3), тогда сила действия на частицу в точке (𝑥, 𝑦) будет равна 𝐅(𝑥, 𝑦) = ∑ W(x − x𝑝 , y − y𝑝 ) 𝐄(𝑥𝑝 , 𝑦𝑝 ) (2.13) 𝑝 где сумма берется по всем узлам сетки. 2.1.5 Интегрирования уравнений движения частиц Для решения системы движения заряженных частиц (1.7) используется явный метод с перешагиванием [1,10] 𝐩𝐤+𝟏⁄𝟐 − 𝐩𝐤−𝟏⁄𝟐 𝑞𝐄 𝐤 (𝐫𝐢 ) = Δ𝜏 𝑚0 𝑐 2 (2.14) 𝐫 𝐤+𝟏 − 𝐫 𝐤 𝐩𝐢 𝐤+𝟏⁄𝟐 = 𝑘+1⁄2 Δ𝜏 𝛾𝑖 где 𝑘 – номер шага интегрирования, Δ𝜏 – шаг интегрирования по времени В случае декартовых координат (2.14) примет вид 𝑞 p𝑥 𝑘+1⁄2 = p𝑥 𝑘−1⁄2 + Δ𝜏 𝐸𝑘 𝑚0 𝑐 2 𝑥 p𝑦 𝑘+1⁄2 = p𝑦 𝑘−1⁄2 + Δ𝜏 𝑘+1⁄2 𝑥 𝑘+1 = 𝑥 𝑘 + Δ𝜏 𝑦 𝑘+1 p𝑥 𝛾 𝑘+1⁄2 𝑞 𝐸𝑘 𝑚0 𝑐 2 𝑦 (2.15) p𝑦 𝑘+1⁄2 = 𝑦 + Δ𝜏 𝑘+1⁄2 𝛾 𝑘 здесь p𝑥 , p𝑦 – x и y компоненты импульса частицы, 𝑥, 𝑦 – ее декартовы 17

координаты. Рисунок 2.1: Схема метода с перешагиванием. Видоизменено [10]. Метод построен на идее центрирования по времени: вычисление приведенного импульса происходит в момент времени, отстоящий на половину временного шага от вычисления вектора фазовых координат. Метод достаточно прост и устойчив, имеет второй порядок точности. Достоинством алгоритма так же является хорошее сохранение энергии при расчете длинных траекторий [12]. Выбор величина шага интегрирования производится согласно критерию устойчивости Куранта С= |𝐮|∆𝑡 ≤ 𝐶𝑚𝑎𝑥 ℎ (2.16) 𝐶𝑚𝑎𝑥 ≤ 1 где С – число Куранта, ∆𝑡 – временной шаг интегрирования, ℎ - размер ячейки пространственной сетки, |𝐮| – скорость переноса частицы. В двумерном случае (2.16) примет вид С= 𝑢𝑥 ∆𝑡 𝑢𝑦 ∆𝑡 + ≤ 𝐶𝑚𝑎𝑥 ℎ𝑥 ℎ𝑦 (2.17) здесь 𝑢𝑥 , 𝑢𝑦 – x и y компоненты скорость переноса частицы, ℎ𝑥 , ℎ𝑦 – ширина и высота ячеек пространственной сетки. Таким образом, при 𝐶𝑚𝑎𝑥 = 1 критерий означает, что за один временной шаг частица не должна продвинуться больше, чем на величину шага пространственной сетки. 18

Старт Инициализация сетки Расчет плотности тока по закону Чайлда-Ленгмюра Добавление новых частиц в область эмиссии Раздача плотности пространственного заряда на сетку Решение уравнение Пуассона Расчет x-,y-компонент напряженности Интерполяция сил в местоположение частиц Интегрирование уравнений движения Проверка пересечения частицами границ расчетной области нет Проверка условий окончания расчета да Стоп Рисунок 2.2: Схема метода частиц в ячейках 19

2.2 Итерационный метод Итерационный метод [1,4,5] относится к группе методов "частицасетка". В отличие от метода частиц в ячейках, итерационный метод используется только для моделирования стационарных потоков заряженных частиц. Итерация метода состоит из следующих этапов. Сперва происходит инжекция частиц в расчетную область. Затем рассчитываются полные траектории частиц под действием электромагнитных полей, рассчитанных на предыдущей итерации. На каждом временном шаге плотность заряда от частиц накапливается в узлах расчетной сетки. После того как все частицы покинут расчетную электромагнитных область, полей. рассчитывается Процедура новая повторяется до конфигурация сходимости к стационарному решению. Главными структурными отличиями итерационного метода от метода частиц в ячейках являются: 1. Инжекция частиц происходит единожды на каждой итерации (в методе частиц в ячейках эмиссия частиц происходит на каждом временном шаге) 2. 2.2.1 Способ распределения заряда частицы на расчетную сетку Схема раздачи заряда в итерационном методе Рассмотрим равномерную сетку с прямоугольными ячейками размером [ℎ𝑥 , ℎ𝑦 ]. Пусть за временной шаг ∆𝑡 модельная частица переместилась внутри ячейки расчетной сетки из точки (𝑥1 , 𝑦1 ) в точку(𝑥2 , 𝑦2 ). Тогда схема распределения заряда частицы моделью трубок тока [1] на расчетную сетку выглядит следующим образом 20

𝑁𝑝 𝜌(𝑥, 𝑦) = ∑ 𝑖=1 𝐼𝑖 ∆𝑡 [W(x2 − x𝑝 , y2 − y𝑝 ) + W(x1 − x𝑝 , y1 − y𝑝 )] 2 ℎ𝑥 ℎ𝑦 (2.18) где 𝐼𝑖 – ток соответствующей траектории частицы, W – ранее определенная (2.3) интерполяционная функция В случае если за временной шаг частица перешла из одной ячейки расчетной сетки в другую, необходимо: 1. Определить координаты пересечения границ ячеек частицей 2. Распределить заряд по всем пройденным ячейкам с помощью модели трубок тока (2.18) 2.2.2 Подавление вклада погрешности в решение Распределение пространственного заряда определяется с некоторой погрешностью. Для уменьшения вклада погрешности, используется нижнее релаксационное преобразование [1,4,12] 𝜌ℎ𝑛 = (1 − 𝜔)𝜌ℎ𝑛−1 + 𝜔 𝜌ℎ𝑛 (2.19) где 𝜔 𝜖 (0,1] – параметр нижней релаксации. На первой итерации параметр релаксации полагается равным 1.0. Малые значения параметра релаксации усредняют значение плотности заряда последних итераций, что позволяет уменьшить вклад погрешности в вычисление плотности тока, но при этом замедляет сходимость. В качестве критерия остановки [1] итерационного процесса используется 𝜌𝑠𝑛 − 𝜌𝑠𝑛−1 max | |<𝜀 𝑠=1..𝑁 𝜌𝑠𝑛−1 где 𝑁 – количество узлов сетки, 𝜀 – наперед заданная точность. 21 (2.20)

Старт Инициализация сетки Расчет электромагнитных полей без учета пучка Добавление новых частиц в область эмиссии Расчет траекторий частиц Накопление пространственного заряда Расчет электромагнитных полей нет Проверка критерия сходимости да Стоп Рисунок 2.3: Схема итерационного метода 22

2.3 Отражение электронов В случае, когда внутренний проводник анода изготовлен из материала с высоким атомным номером эффект обратного отражения электронов становится существенным. Данное явление заключается в том, что часть электронов, падающих на анод отражается обратно, потеряв часть энергии и заряда. Отраженные электроны попадают в зазор "анод-катод". Затем они пролетают часть зазора в сторону катода, останавливаются и снова попадают на анод, после чего процесс повторяется. В результате это приводит росту пространственного заряда в прианодной области, изменению конфигурации электромагнитных полей. В работе [14] подробно рассмотрен способ моделирования данного эффекта. Пусть 𝛼 – коэффициент потери заряда, 𝛽 – коэффициент потери энергии частицей при отражении. Где 𝛼, 𝛽 зависят от энергии электронов и от материала, из которого изготовлен анод. В случае анода из стали [7] коэффициенты 𝛼, 𝛽 подчиняются нормальным законам распределения с математическими ожиданиями равными E[𝛼] = 0.8 и E[𝛽] = 0.3, среднеквадратическими отклонениями 𝜎𝛼 = 0.8⁄3 и 𝜎𝛽 = 0.1 соответственно. Рисунок 2.4: Отражение электронов в плоском диоде 𝛼 = 0.7 𝛽 = 0.7 23

Глава 3. Сглаживание сеточных величин Одним из главных недостатков методов группы "частица-сетка" является проблема появления счетных "нефизических" шумов. Самым простым способом уменьшения численных шумов является увеличение числа модельных частиц, что ведет к резкому увеличению требуемых вычислительных ресурсов. Для уменьшения счетных шумов изменяют схему раздачи заряда частиц, подбирают оптимальные временной шаг и размер ячейки расчетной сетки. Начальное распределение частиц также играет немаловажную роль. Наиболее распространенными подходами распределения частиц являются системы "тихий" старт и случайный старт. Принцип случайного старта заключается в распределении частиц по области эмиссии методами Монте-Карло. Обычно начальное распределение представляет собой совместное распределение нескольких законов распределения (например, максвелловское распределение по скорости, равномерное по координате и др.) Использование схемы "тихого" старта обеспечивает снижение шума и неопределенности случайного старта. Вся область эмиссии покрывается набором частиц так, чтобы воспроизвести необходимое распределение (например, равномерное распределение может быть воспроизведено помещением набора частиц, расположенных на одинаковом расстоянии друг от друга). Однако широкого распространения схема "тихий старт" не получила [10] по следующим причинам: подход требует большего количества модельных частиц и более сложен в реализации (например, в случае начальной дискретизации по множеству измерений), так же использование регулярно расположенного множества частиц может привести к появлению "нефизических" корреляций среди частиц. Возникновение счетных эффектов приводит к большим флуктуациям рассчитываемых величин, в особенности плотности заряда. В случае же 24

итерационного метода это приводит к плохой сходимости метода. В данной главе рассмотрены различные техники сглаживания распределения пространственной плотности заряда в методах группы "частица-сетка" для подавления счетного шума. 3.1 Линейные сглаживающие фильтры Простейшим способом сглаживания является применение линейных сглаживающих фильтров. Их сущность заключается в усреднении значений элементов по их окрестности. Для вычисления значения используется матрица, называемая ядром свертки, содержащая весовые коэффициенты. Пусть 𝜌(𝑘, 𝑙) – значения распределения пространственного заряда в узлах эйлеровой сетки, 𝑤(𝑠, 𝑡) – ядро свертки, тогда схема фильтрации пространственного заряда выглядит следующим образом 𝑎 𝑏 𝜌(𝑘, 𝑙) = ∑ ∑ 𝑤(𝑠, 𝑡)𝜌(𝑘 + 𝑠, 𝑙 + 𝑡) (3.1) 𝑠=−𝑎 𝑡=−𝑏 где 𝑤(𝑠, 𝑡) – нормированное ядро свертки, суммы вычисляются для всех элементов 𝜌(𝑘, 𝑙) Обычно ядро свертки заполняется по нормальному закону распределения. Ниже приведен пример матрицы размерности 5x5: 2 4 5 4 2 4 9 12 9 4 1 𝑤= × 5 12 15 12 5 159 4 9 12 9 4 5 4 2) (2 4 25 (3.2)

3.2 Частотная фильтрация Частотная фильтрация основана на применении к пространственному заряду преобразования Фурье [15] и состоит из следующих этапов: 1. Прямое двумерное преобразование Фурье для получения спектра пространственного заряда 𝑁−1 𝑀−1 𝜌̅ (𝑤1 , 𝑤2 ) = ∑ ∑ 𝜌(𝑘, 𝑙) 𝑒 −2 𝜋𝑖 𝑘𝑤1 𝑙𝑤 −2 𝜋𝑖 2 𝑁 𝑒 𝑀 (3.3) 𝑘=0 𝑙=0 где 𝑤1 , 𝑤2 – угловые частоты 2. Центрирование Фурье образа и его умножение на частотную функцию фильтра 𝐻(𝑤1 , 𝑤2 ) 3. Обратное двумерное преобразование Фурье 𝑁−1 𝑀−1 𝑘𝑤1 𝑙𝑤2 1 𝜌(𝑘, 𝑙) = ∑ ∑ 𝜌̅ (𝑤1 , 𝑤2 ) 𝑒 −2 𝜋𝑖 𝑁 𝑒 −2 𝜋𝑖 𝑀 𝑁𝑀 (3.4) 𝑤1 =0 𝑤2 =0 Рисунок 3.1: Применение сглаживающей оконной функции Гаусса 𝜎 = 0.1 к спектру плотности пространственного заряда 26

3.3 Сглаживающий кубический сплайн Сглаживающий кубический сплайн 𝑠(𝑥) определяется как сплайн, минимизирующий следующий функционал: 𝑥𝑛 𝑛 𝐼(𝑠, 𝑝) = 𝑝 ∑ 𝑤𝑘 (𝑦𝑘 − 𝑠(𝑥𝑘 𝑘=1 ))2 𝑑 2 𝑠(𝑥) + (1 − 𝑝) ∫ ( ) 𝑑𝑥 𝑑𝑥 2 (3.3) 𝑥1 где (𝑥𝑘 , 𝑦𝑘 ) – приближаемые данные, 𝑤𝑘 – веса данных, 𝑝 – сглаживающий параметр. При 𝑝 = 1 сглаживающий сплайн превращается в обычный кубический сплайн, при 𝑝 = 0 данные аппроксимируются в смысле метода наименьших квадратов. При применении к матрице пространственного заряда сглаживание применяется сначала к каждому столбцу матрицы, затем к ее строкам. Результат сглаживания при 𝑝 = 0.6 на расчетной сетке 128 × 128 отображен на рисунке 3.2. Рисунок 3.2: Распределения плотности заряда (до и после сглаживания кубическими сплайнами) 27

Глава 4. Программная реализация численных алгоритмов Рассмотренные численные алгоритмы реализованы на платформе .NET 4.6 с использованием языка C# 6.0. Многопоточность реализована с помощью библиотеки Task Parallel Library (TPL). Для решения отдельных задач дополнительно использовались Math.Net Numerics, Intel MKL и FFTW. Модули алгоритмов реализованы с помощью интерфейсов, это позволяет легко расширять систему в случае необходимости введения новых схем раздачи заряда, методов решения уравнений поля и др. Для разработки использовалась Тестирование интегрированная алгоритмов среда проводилось на Microsoft Visual Studio 2015. двухъядерном процессоре Intel Core i3 2350M. Графический интерфейс реализован с помощью технологии Windows Presentation Foundation (WPF). Для визуализации результатов расчетов используется библиотека OxyPlot, реализована возможность сохранения результатов в файл. Предусмотрена возможность сохранения и загрузки проектов, содержащих параметры решаемой задачи:  временной шаг интегрирования  геометрия диода  напряжение на аноде  параметры отражения электронов  размеры расчетной сетки  параметр релаксации итерационного метода  количество частиц/траекторий  параметры эмиссии частиц  настройка фильтрации сеточных величин 28

Глава 5. Результаты моделирования Рассмотрим задачу моделирования плоского диода в случае двумерной декартовой геометрии. Значения основных параметров диода приведены в таблице 5.1. Расчеты проводились на сетке размерностью 128 × 128. Выбор шага интегрирования производился согласно критерию устойчивости Куранта (2.16), при 𝐶𝑚𝑎𝑥 = 1⁄2. Плотность тока эмиссии считалась постоянной и определялась из закона Чайлда (1.9). Расстояние от катода до анода, м 𝑑 = 0.1 Длина, м 𝐿 = 0.1 Длина эмиттера, м 𝐿𝑒𝑚 = 0.02 Напряжение на аноде, В 𝑈𝑎 = 100 000 Плотность тока эмиссии, А/м 𝐽𝑒𝑚 − 7368 Таблица 5.1: Параметры тестового диода 5.1 Тихий старт. Сравнение методов При начальном распределении электронов по схеме тихого старта итерационный метод показал хорошую сходимость. Сходимость определялась по критерию (2.20). Была проведена серия расчетов при 𝜀 = 10−10 , количество модельных частиц – 500. Результат сравнения скорости сходимости итерационного метода при разных параметрах релаксации представлен в таблице 5.2. Траектории электронов представлены на рисунке 5.1. Результаты расчета с помощью метода частиц в ячейках и итерационного метода достаточно точно согласуются. На рисунках 5.2, 5.3 представлено сравнение распределений плотности заряда и напряженности электрического поля. Однако, так как метод частиц в ячейках требует пересчет электрических полей на каждом временном шаге, расчет с помощью итерационного метода происходит значительно быстрее. Расчет методом частиц в ячейках потребовал 290 временных шагов и занимал 205 секунд. 29

Параметр релаксации 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Количество итераций 329 155 98 71 52 40 31 Время, с 100 48 32 23 19 17 15 Таблица 5.2: Сходимость итерационного метода с разными параметрами релаксации Рисунок 5.1: Траектории частиц в плоском диоде, рассчитанные итерационным методом с использованием тихого старта 30

Рисунок 5.2: Распределение плотности пространственного заряда в середине диода Рисунок 5.3: Распределение x-компоненты напряженности на катоде 31

5.2 Случайный старт. Применение сглаживания 5.2.1 Случай итерационного метода Использование схемы случайного старта приводит к появлению больших осцилляций в распределении плотности пространственного заряда. В случае итерационного метода сходимости по критерию (2.20) добиться не удалось, сглаживание пространственного заряда также не помогло добиться сходимости (рисунок 5.4). Рисунок 5.4: Большие осцилляции в плотности пространственного заряда приводят к расходимости (критерий 2.1) итерационного метода Однако результаты, рассчитанные с использованием случайного старта, достаточно точно согласуются с результатами, полученными с начальным распределением электронов "тихий старт" (рисунок 5.5). Поэтому для оценки сходимости итерационного метода со случайным стартом будем использовать решение, полученное со схемой тихого старта, как эталонное. Сравнение будем производить по следующему критерию 32

𝑛 𝑛−1 max ∑|𝜌𝑖𝑗 − 𝜌𝑖𝑗 | 𝑖 (5.1) 𝑗 где n – номер итерации, (i,j) – координаты узла сетки. Рисунок 5.5: Распределение плотности пространственного заряда в итерационном методе Примечание: при расчетах использовалось 5000 модельных частиц, параметр релаксации - 0.1 Серией расчетов было определено, что точность вычислений при использовании случайного старта напрямую зависит от количества модельных частиц (рисунок 5.6). Было выявлено, что использование параметра релаксации больше 0.1 нецелесообразно и приводит к большим колебаниям в решении. Также расчеты показали, что критерий (5.1) позволяет сравнивать результаты расчетов, однако не позволяет количественно оценить точность вычислений (рисунок 5.7). Таким образом при расчете итерационным методом случайный старт электронов вносит слишком большие колебания в пространственный заряд, что ухудшает сходимость метода, поэтому предпочтительным способом начального распределения электронов является схема "тихий старт". 33

Рисунок 5.6: Отклонение плотности пространственного заряда в итерационном методе со случайным стартом от эталонного значения Рисунок 5.7: Оценка сходимости итерационного метода по критерию (5.1) 34

Случай метода частиц в ячейках 5.2.2 Рассмотрим теперь влияние случайного старта при расчетах методом частиц в ячейках. Исследуем также эффективность сглаживания пространственной плотности заряда в методе частиц в ячейках. Для оценки точности расчетов будем использовать решение итерационного метода, полученное со схемой "тихий старт" как эталонное. На рисунке 5.8 представлено сравнение отклонений распределений пространственного заряда со сглаживанием кубическими сплайнами и без от эталонного решения. На рисунке 5.9 произведено сравнение эффективности алгоритмов сглаживания (3.1), (3.2), (3.3), результаты представлены как величина среднего отклонения от эталонного значения. Рисунок 5.8: Отклонение распределения плотности заряда со сглаживанием кубическими сплайнами s=0.5 и без от эталонного решения, полученного итерационным методом. 35

Без сглаживания Кубические сплайны Линейная фильтрация Фурье сглаживание 6E-06 5E-06 k=5 σ=0.8 4E-06 3E-06 s=0.5 σ=0.3 k=7 s=0.8 σ=0.5 k=11 s=0.1 s=0.3 σ=0.1 2E-06 1E-06 0E+00 Рисунок 5.9: Сравнение эффективности алгоритмов сглаживания. На диаграмме представлено сравнение среднего отклонения пространственного заряда от эталонного значения, рассчитанного итерационным методом с тихим стартом Примечание: где s – сглаживающий параметр кубических сплайнов, k – размер ядра свертки линейного сглаживания, σ – параметр оконной функции Гаусса при сглаживании. Таким образом, применение сглаживания к пространственной плотности заряда при расчете методом частиц в ячейках, позволяет заметно уменьшить осцилляции, вносимые случайным стартом, тем самым улучшив результаты расчета. 36

Заключение В процессе работы получены следующие результаты: 1. Создана программная реализация алгоритмов для расчета динамики стационарных интенсивных потоков заряженных частиц в диодах. Произведен расчет динамики стационарного интенсивного пучка заряженных частиц в плоском диоде итерационным и методом частиц в ячейках, проведено их сравнение. 2. Итерационный метод с моделью трубок тока является более экономичным и эффективным, по сравнению с методом частиц в ячейках, инструментом моделирования стационарных состояний пучка частиц в случае отсутствия необходимости дискретизации начальных данных по множеству измерений. 3. Использование случайного начального распределения модельных частиц в итерационном методе приводит к большим колебаниям в распределении пространственного заряда и расходимости метода. 4. При расчете методом частиц в ячейках подавление счетного шума сглаживанием пространственной плотности заряда позволяет заметно улучшить точность расчета, снизив осцилляции, вносимые случайным стартом. 37

Список литературы 1. Алцыбеев В. В. Оптимизационный алгоритм расчета плотности тока эмиссии // Вестник Санкт-Петербургского университета. Серия 10. Прикладная математика. Информатика. Процессы управления. 2015. № 4. С. 56–72 2. Месяц, Е. А. Методы оценки и повышения точности решения задач физики плазмы методом частиц в ячейках: Ph.D. thesis / Институт вычислительной математики и математической геофизики СО РАН. — Новосибирск, 2014 — 110 с. 3. Муравьев, А. Г. Математическое моделирование электронных пушек с катодом произвольной формы: Ph.D. thesis / Московский физикотехнический институт. — Москва — 110 с. 4. Свешников, В. М. Численное моделирование интенсивных пучков заряженных частиц. — 2006. 5. César C. Xavier and Cláudio C. Motta. The XMGUN Particle Path FEM Code IEEE Transactions on Magnetics, vol. 46, No. 8, August 2010 6. Child, C. D. Discharge from hot cao / C. D. Child // Phys. Rev. Series I. — 1911. — Vol. 32, № 5. — pp. 255–282. 7. Engelko V., Kuznetsov V., Viazmenova G., Mueller G., and Bluhm H., Influence of electrons reflected from a target on the operation of triode-type electron sources // Journal of Applied Physics 88, 3879–3888 (2000). 8. Hockney, R. Computer Simulation Using Particles / R. Hockney, J. Eastwood. — Francis, 1988. — 540 pp. 9. Jolliet S., Bottino A., Angelino P., et al. A global collisionless PIC code in magnetic coordinates // Computer Physics Communications. — 2007. —Vol. 177. — P. 409-425. 10.Lapenta G., Particle In Cell Method A brief description of the PIC Method // Centrum voor Plasma Astrofysica Katholieke Universiteit Leuven 38

11.Mesyats E.A. A noise-reducing algorithm for Particle-in-Cell plasma simulation // Bull. Nov. Comp. Center, Num. Anal., 14 (2009), pp. 21-30 12.Mudiganti, J. C. An Emission Model for the Particle-in-Cell Method: Ph.D. thesis / Darmstadt: TU Darmstadt. — 2006. — 123 pp. 13.Nagel J. R., Solving the Generalized Poisson Equation Using the FiniteDifference Method (FDM) — February 15, 2012 —18pp 14.Oliver B.V., Genoni T.C., D.V. Rose, D.R. Welch.Space-charge limited currents in coaxial diodes with electron backscatter // 2001 Journal of Applied Physics vol. 90, No. 10, 15 November 2001 15.Sydora R.D. Low-noise electromagnetic and relativistic Particle-in-Cell plasma simulation models // J. Comp. Appl. Math. . — 1999. . — Vol. 109. . —P. 243259. 16.Watrous J.J., Lugisland J.W., and Sasser G.E., An improved space-chargelimited emission algorithm for use in particle-in-cell codes // Phys. Plasmas, vol. 8, pp. 289-296, Jan. 2001. 17.Weisterman T., A particle-in-cell method as a tool for diode simulation // Nuclear Instruments and Methods in Physics Research A 260(1988), pp. 271-279 39

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв