Моделирование и оптимизация динамики интенсивных пучков заряженных частиц тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат наук Алцыбеев Владислав Владимирович

• Алцыбеев Владислав Владимирович
• кандидат наук
• 2016

Оглавление диссертации кандидат наук Алцыбеев Владислав Владимирович

Введение

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

1.1 Система уравнений Власова

1.2 Модели частиц

1.3 Макрочастицы в моделях частица-сетка

1.4 Постановка задачи моделирования потоков заряженных частиц в электростатическом приближении

1.5 Дискретизация задачи на расчетной сетке. Весовые функции

1.6 Концепция метода частиц в ячейках

1.7 Концепция итерационного метода

Глава 2. Численные алгоритмы в методах частица-сетка

2.1 Расчет траекторий частиц

2.1.1 Случай декартовых координат

2.1.2 Случай цилиндрических координат

2.1.3 Алгоритмы выбора шага интегрирования

2.2 Форма макрочастиц, расчет пространственного заряда и сил

2.2.1 Случай метода частиц в ячейках

2.2.2 Случай итерационного метода

2.3 Решение уравнения Пуассона

2.4 Методы расчета тока, ограниченного пространственным зарядом

2.4.1 Модель эмиссии Чайлда-Ленгмюра

2.4.2 Оптимизационный алгоритм в итерационном методе

2.4.3 Модель эмиссии Гаусса в итерационном методе

Глава 3. Комплекс программ для моделирования динамики пучков

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

3.1 Основные характеристики комплекса программ

3.1.1 Формат входных данных

3.1.2 Варианты использования

3.2 Общее описание архитектуры

3.3 Алгоритмы хранения, обработки и удаления частиц

3.3.1 Структура данных хранения частиц

3.3.2 Параллельная реализация методов частица-сетка с использованием ОрепМР

Глава 4. Результаты моделирования эмиссионных устройств

4.1 Цилиндрический диод

4.1.1 Случай цилиндрических координат

4.1.2 Случай декартовых координат

4.2 Диод с эллиптическим эмиттером

4.3 Источник радиально сходящегося пучка электронов триодного типа

4.3.1 Параметры источника

4.3.2 Исследование дефокусирующего эффекта ускоряющей сетки

4.3.3 Траектории и вольт-амперные характеристики

Глава

Моделирование и оптимизация ускорителя с переменно-фазовой фокусировкой

5.1 Ускоритель с ПОКФ

5.2 Расчет геометрии резонатора с ПФФ

5.2.1 Модель динамики пучка

5.2.2 Расчет последовательности синхронных фаз

5.2.3 Генератор трехмерной компьютерной модели резонатора

5.2.4 Методика выбора геометрии периодов

5.3 Результат расчета резонатора на 60 периодов

5.4 Задача оптимизации

5.4.1 Постановка задачи

5.4.2 Метод оптимизации

5.4.3 Результаты численной оптимизации

Заключение

Список литературы

Введение диссертации (часть автореферата) на тему «Моделирование и оптимизация динамики интенсивных пучков заряженных частиц»

Введение

Задачи моделирования и оптимизации динамики плазмы, а также пучков заряженных частиц представляют собой широкую и важную область вычислительной физики, поскольку приблизительно 99 % материи в галактике, включая все звезды, межзвездную и межпланетную среду, верхние слои планетных атмосфер, находится в плазменном состоянии. А также, процессы, протекающие в установках, предназначенных для генерации и ускорения потоков заряженных частиц, описываются теми же законами, что и динамика плазмы. Такие установки (ускорители, токамаки, электронные и ионные источники и т. д.) интенсивно развиваются в последние годы и применяются в широком спектре прикладных задач.

В частности, большое значение приобретают разработки в области построения ускорителей с фокусировкой ускоряющим полем. Подобная фокусировка осуществляется, например, в линейных ускорителях с пространственно-однородной квадрупольной фокусировкой (ПОКФ) и переменно-фазовой фокусировкой (ПФФ). Структуры с ПОКФ давно уже стали основной начальной частью ускорителей на большие энергии. Также они используются и самостоятельно для различных прикладных целей. Однако их использование ограничивается невысоким темпом ускорения в случаях, когда необходимо получать высокую мощность высокочастотного поля для ускорения интенсивного пучка. Поскольку дальнейшее ускорение в таких структурах не столь выгодно, то разрабатывается идея создания ускорителей, где начальной частью служит ускоритель ПОКФ, а далее используется структура с ПФФ. Структура с ПФФ обладает высоким темпом ускорения, что в совокупности с отсутствием фокусировки магнитными полями, делает её эффективнее и дешевле ускорителей типа Альвареца. Проблемам расчета этих структур посвящено много работ [1-24]. Отметим, что при создании и проектировании современных ускорительных комплексов и структур прикладного назначения растут требования к качеству получаемого пучка.

Один из способов обеспечения желаемого результата — применение специального математического аппарата теории управления [25], позволяющего строить эффективные направленные методы оптимизации параметров ускорителя. Задачи управления динамикой пучков заряженных частиц, впервые были поставлены и изучены в работах [10,11]. В частности, различные математические модели оптимизации динамики пучков в структурах с ПОКФ и ПФФ были предложены в [1,2,6,13-16,26,27]. Например, предлагаются модели, позволяющие проводить оптимизацию в поле эквивалентной бегущей волны. Тем не менее, несмотря на значительный прогресс в вопросах разработки ускорителей с ПФФ, при ускорении пучка с высоким значением тока в импульсе, может наблюдаться существенное ухудшение его характеристик по сравнению со слаботочным пучком. Общих универсальных методик выбора параметров резонаторов ускорителей с ПФФ, обеспечивающих минимизацию эффектов этого явления в настоящее время не существует. Поэтому актуальными являются проблемы разработки данных ускорителей и улучшения качества пучков в них.

Различные источники электронов могут применяться, например, для облучения различных мишеней с целью обработки их поверхностей [28-31]. В таких случаях электронный поток может вызывать эмиссию ионов с поверхности мишени [32]. При этом, зачастую, ток эмиссии электронов и ионов в таких установках ограничен пространственным зарядом, таким образом повышение электронного тока вызывает увеличение ионного тока, что в свою очередь вызывает увеличение электронного тока и так далее. При нормальных режимах работы, после стадии колебаний токов и пространственного заряда, источник достигает стационарного состояния, которое продолжается в течение длительности импульса генератора напряжения. Однако различные физические эффекты (ионная эмиссия, отражение электронов, собственное магнитное поле пучка) могу приводить к существенным колебаниям токов и мощностей на облучаемых объектах в течение импульса и нарушению нормального режима работы. Таким образом, при разработке конструкции источника, необходимо проводить численный

анализ динамики потоков заряженных частиц с целью нахождения параметров, обеспечивающих стабильные режимы работы.

Как в случае линейных ускорителей, так и источников электронов, динамика пучков заряженных частиц описывается с помощью системы уравнений Власова-Максвелла или Власова-Пуассона в случае безвихревого электрического поля. В обоих случаях необходимо корректно учитывать сложную геометрию установок, что ведет к необходимости использования различных систем координат при проведении расчетов (декартовы, цилиндрические, полярные координаты), уравнения Максвелла и уравнения движения в которых имеют раз-

VJ ТЛ VJ

ный вид. В настоящее время численные методы частиц, в частности эйлерово-лагранжевые методы частица-сетка является наиболее эффективным средством для решения уравнения Власова [33-35]. Число макрочастиц в процессе моделирования методами частиц может достигать порядков 109, а интервалы расчета порядков микросекунд, что налагает серьезные требования к производительности прикладных программ. Таким образом, стоит отметить, что актуальной и одной из важнейших задач, стоящей перед российскими исследователями и разработчиками является разработка отечественных программ моделирования динамики заряженных пучков и плазмы, которые могут составить конкуренцию коммерческим зарубежным аналогам в производительности и широте применения, таким как VORPAL, CST Particle Studio, Magic.

Среди методов частица-сетка наиболее широко используется метод частиц в ячейках [36-41]. В случае независимости от времени рассматриваемого процесса (например, стационарное состояние источника), возможно применять более экономичный итерационный метод [42-49]. Данный подход основан на повторяющемся процессе расчета траекторий частиц, заряда, вносимого ими в пространство и расчете электрического поля до достижения сходимости итераций процесса. На данный момент существуют несколько программных пакетов, использующих данный подход, например CST particle studio [50] и POISSON [51].

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

Требуемое число модельных частиц для метода частиц в ячейках значительно больше (может достигать порядков 109), чем требуемое число частиц для итерационного метода, поскольку в случае итерационного метода частицы инжектируются в расчетную область единственный раз на каждой итерации, в отличии от инжекции частиц на каждом временном шаге в случае метода частиц в ячейках.

Обычно в задачах с эмиссией, ограниченной пространственным зарядом совместно с итерационным методом применяются модели, основанные на одномерных аналитических решениях Чайлда и Ленгмюра [42,43,45,46,50,51]. Но, случае криволинейной эмиттирующей поверхности, применение этого подхода может привести к возникновению существенных ошибок в решении [42, 45]. Этот факт составляет основной недостаток итерационного метода.

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

Цель диссертационной работы заключается в развитии методов математического моделирования и оптимизации динамики пучков заряженных частиц. Для достижения поставленной цели были решены следующие задачи:

1. Разработан метод расчета плотности тока эмиссии в режиме ограничения пространственным зарядом при использовании итерационного метода решения стационарных самосогласованных задач. Предлагаемое улучшение, основанное на применении модели эмиссии Гаусса, позволяет более точно определять распределение плотности тока эмиссии в случае криволинейной эмиттирующей поверхности по сравнению с моделью эмиссии Чайлда-Ленгмюра. Также, в случае наличия различных физических эффектов, таких как дополнительные потоки и отражение электронов, предлагаемый подход позволяет получить корректные результаты.

2. Произведены расчеты динамики пучков в различных эмиссионных устройствах, иллюстрирующие эффективность разработанного метода расчета плотности тока эмиссии. Получены условия достижения необходимой плотности мощности для модификации цилиндрических мишеней в источнике электронов три-одного типа ГЕЗА 4м.

3. Рассмотрена новая постановка задачи оптимизации динамики пучка траекторий. Особенностью данной задачи является то, что программное управление является кусочно-постоянной функцией, в которой точки переключения управления зависят от самих значений функции. Для оценки качества пучка вводится интегральный функционал. Получено аналитическое представление градиента функционала качества по управляющим параметрам и формулируются необходимые условия оптимальности.

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

тическом приближении в ускоряющих структурах до и после оптимизации при разной интенсивности пучка.

5. Разработан комплекс программ для моделирования динамики пучков заряженных частиц в электростатическом приближении а также для расчета и оптимизации параметров ускорителей с переменно-фазовой фокусировкой, который позволяет поддерживать возможность использования различных систем координат при расчетах, таких как декартовые двумерная и трехмерная, цилиндрическая с учетом осевой симметрии, полярная. Комплекс программ создан с использованием языка С++. Процесс вычислений оптимизирован под работу на системах с общей памятью с использованием стандарта ОрепМР, проведены исследования масштабируемости. Разработанный комплекс программ реализует разработанную математическую модель оптимизации, позволяет производить расчет геометрии резонатора с учетом его настройки на расчетную частоту, а также автоматически создавать его трехмерную компьютерную модель.

Методы исследования. Для решения задач, рассматриваемых в диссертации используются методы математического моделирования, вычислительной математики, математической физики, дифференциальных уравнений, физики пучков заряженных частиц, теории управления и оптимизации,объектно-ориентированного программирования, высокопроизводительных вычислений.

Основные положения и результаты, выносимые на защиту:

1. Метод расчета плотности тока эмиссии, ограниченного пространственным зарядом, для итерационного метода решения самосогласованных электростатических задач.

2. Результаты расчетов динамики пучков в различных эмиссионных устройствах (в том числе в источнике электронов триодного типа ГЕЗА 4м).

3. Математическая модель и алгоритм оптимизации динамики пучка траекторий, в которой программное управление является кусочно-постоянной функцией с точками переключения, зависящими от значений функции.

4. Результаты оптимизации параметров ускорителя дейтронов с переменно-фазовой фокусировкой.

5. Комплекс программ моделирования динамики пучков заряженных частиц в электростатическом приближении и для расчета параметров ускорителей с переменно-фазовой фокусировкой.

Научная новизна.

Разработан метод расчета плотности тока эмиссии в режиме ограничения пространственным зарядом при использовании итерационного метода решения стационарной самосогласованной задачи. Эффективность метода (возможность получать более точные решения в случае криволинейной эмиттирующей поверхности) по сравнению со стандартным подходом (модель Чайлда-Ленгмюра) подтверждена на примерах моделирования различных эмиссионных устройств.

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

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

Все результаты, изложенные в оригинальной части диссертационной работы, получены впервые и являются новыми.

Научная и практическая ценность работы состоит в создании эффективного метода расчета плотности тока эмиссии с криволинейной эмиссионной поверхности при использовании итерационного метода, разработке математической модели и алгоритма оптимизации динамики пучка в случае функции программного управления специального вида, а также в создании необходимого для расчетов программного обеспечения. Разработанные модели, методы и программы могут быть использованы при расчете различных эмиссионных устройств и ускорителей, а также при проектировании и оптимизации ускорителей с переменно-фазовой фокусировкой.

Степень достоверности и апробация результатов. Основные положения и результаты диссертационной работы докладывались и обсуждались на следующих конференциях и семинарах: IVESC-2014 (International Vacuum Electron Sources Conference, 2014, Санкт-Петербург); BDO-2014 (Beam Dynamics and Optimization, 2014, Санкт-Петербург); LINAC14 (27th Linear Accelerator Conference, 2014, Женева); ВСПУ-2014 (XII Всероссийское совещание по проблемам управления, 2014, Москва); XI Международный семинар по проблемам ускорителей заряженных частиц памяти В. П. Саранцева (2015, Алушта); SCP-2015 (III International conference stability and control processes, 2015, Санкт-Петербург).

Публикации.

Основные положения диссертации изложены в 15 опубликованных в печати работах [19,21,22,26,27,38,52-60], из них 4 статьи в рецензируемых журналах, входящих в перечень ВАК [55-58]. Работы [19,22,26,27,38,52,53] опубликованы в изданиях, индексируемых в Scopus.

Личный вклад автора. Все положения, выносимые на защиту, получены лично автором.

Структура диссертационной работы.

Диссертация состоит из введения, пяти глав, заключения и списка литературы. Объем работы составляет 130 страниц, среди них 5 таблиц и 59 рисунков. Список литературы включает 109 наименований.

Глава 1

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

1.1 Система уравнений Власова

Рассмотрим расчетную область О = О и Г, где Г — граница расчетной области. Потоком частиц в дальнейшем будем называть множество частиц одного типа, эмиттируемых в расчетную область О с одной непрерывной поверхности или кривой старта (эмиттера). Массу, массу покоя и заряд частиц в потоке с номером а будем обозначать далее как та, т0а и да соответственно. Рассмотрим модели описания динамики пучка заряженных частиц, который может состоять из нескольких потоков частиц разного типа.

Пусть па(1,х,у,х,рх,ру,рг) = иа(1, г, р) — функция распределения плотно-

сти частиц в потоке с номером а = 1... N в пространстве координат г и импульсов р. Эти функции удовлетворяют уравнению Власова [61-63]

П = дпа + дпа v + дпа р = 0 (12)

дЛ, дЬ дг др а

Здесь Еа = да (Е + [уБ]) — сила Лоренца; E — напряженность электрического поля; B — индукция магнитного поля; v — скорость.

Уравнение Власова является математической моделью, описывающей динамику пучка взаимодействующих заряженных частиц или плазмы с учётом даль-нодействующих кулоновских сил. Компоненты электромагнитного поля опреде-

ляются из уравнений Максвелла

&уЕ = р/е,

а1уБ = о,

дВ (1.3)

го1Е = ——,

дг

д Е

го1Б = —— + до. с дг

Здесь е и ц — диэлектрическая и магнитная проницаемости материала среды; р и j — объемная плотность пространвенного заряда и вектор плотности тока, определяемые через фазовую плотность частиц

р(г, г) = £ д^ Па (г, г, р)ф, (1.4)

j(г, г) = £ да VПа(г, г, (1.5)

а

Электромагнитное поле в системе уравнений (1.2), (1.3) является самосогласованным. Для решения уравнения Власова (1.2) необходимо определить распределения компонент электромагнитного поля из уравнений Максвелла (1.3), которые, в свою очередь, зависят от решения уравнения Власова (1.2) согласно (1.4), (1.5). В общем случае, система уравнений Власова-Максвелла является нелинейной, при этом полное её аналитическое исследование невозможно [61,62,64-67]. В связи с этим, интенсивное развитие получили различные численные методы решения уравнения Власова — методы частиц, конечно-разностные методы, метод конечных элементов.

1.2 Модели частиц

Для численного решения системы уравнений Власова (1.2)-(1.5) применяются три типа вычислительных алгоритмов — эйлеровы алгоритмы, лагранжевы и смешанные, эйлерово-лагранжевы алгоритмы.

В подходе Эйлера ведется наблюдение за изменением характеристик моделируемой среды в фиксированных точках пространства, например в узлах непо-

^ ^У / \ Т/*

движной расчетной (эйлеровой) сетки. К данной группе методов, например относятся методы конечных разностней и конечных элементов [34,68,69]. Однако, широкого распространения в задачах моделирования динамики пучков и плазмы данные методы не получили, причиной чего стало большое количество недостатков данных методов, например большой требуемый объем вычислений, рекуррентные явления, неточность работы для областей с большими градиентами функции распределения.

Более широкое распространение при практических расчетах получил подход Лагранжа, при котором необходимо следить за движением различных выбранных частиц, и смешанный подход. Примерами применения данных подходов к исследованию динамики пучков заряженных частиц и плазмы являются модели частиц [33,35,39]. Выделяют три основных типа вычислительных моделей частиц: модель частица-частица, частица-сетка и модель частица-частица — частица-сетка [33].

Метод частица-частица является простейшим способом учета взаимодействия заряженных частиц с понятийной точки зрения. В таком случае пучок частиц представляется в виде макрочастиц определенной формы, зависящей от размерности задачи и свойств симметрии рассматриваемой системы. Например к такой группе моделей относятся модель плоских листов, дисковые и кольцевые модели [12,70-72]. Частицы могут обладать нулевым объемом (например модель тонких дисков), или иметь конечный объем (напрмиер модель толстых дисков). Использование частиц нулевого размера приводит к завышению влияния ближних взаимодействий, или короткодействующих сил. В то время как для частиц конечного объема сила взаимодействия при их сближении сначала возрастает, а затем, при проникновении частиц друг в друга снижается до нуля, что означает снижение влияния ближних взаимодействий. Схематично сравнение силы взаимодействия частиц конечного и нулевого размеров представлено

на рис. 1.2. Существенным недостатком данного метода является его высокая вычислительная сложность — 0(Ы2), где N — число модельных частиц. Однако, в ряде случаев (например кольца или набор дисков в полой металлической трубе), с помощью данного метода удается получить аналитические выражения для кулоновского поля пучка взаимодействующих частиц. Данные модели нашли свое применение в задачах оптимизации динамики пучков в линейных ускорителях заряженных частиц [10,11,19,26,27,73]. Метод частица-частица является лагранжевым алгоритмом.

В методе частица-сетка за-__ ряженная среда моделируется

усредненными по некоторому правилу значениями плотностей заряда и тока в узлах эйлеровой сетки, получаемых из расчетов траекторий макрочастиц. Распределения электромагнитных полей вычисляют-

-Сила Кулона

-----Сила взаимодействия

частиц конечного размера

* \ / /

/ /

5 г/А 6 г

Рисунок 1.2: Взаимодействие частиц. Здесь г — ся в узлах сетки с п°м°щью расстояние между частицами, Дг — размер решения уравнения Пуассона частицы или полной системы уравне-

ний Максвелла. Сила вычисляется значительно быстрее, чем в методе частица-частица. Ситуация в данном случае будет аналогична случаю, представленному на рис. 1.2. Однако для систем с преобладающими дальнодействующими силами, этот недостаток не вносит существенных ошибок в расчеты. Дополнительно, модели частиц позволяют естественным образом учесть многие физические эффекты, такие как эмиссия и поглощение частиц, их рассеяние или отражение, возникающие при изучении источников заряженных частиц или ускорителей. Таким образом, методы частица-сетка относятся к смешанным эйлерово-лагранжевым подходам. К та-

кой группе методов относятся метод частиц в ячейках и итерационный метод решения стационарных задач.

Смешанный метод частица-частица — частица-сетка может применяться при необходимости учета ближних взаимодействий, при этом сохраняя простоту вычисления дальнего взаимодействия в методе частица-сетка. Действующая между частицами сила расщепляется на две части — быстроменяющаяся близкодействующая часть, отличная от нуля только на нескольких межчастичных расстояниях и медленноменяющаяся дальнодействующая часть, достаточно гладкая для точного представления на расчетной сетке. Метод частица-частица используется для нахождения близкодействующей силы для каждой частицы, а метод частица-сетка используется для вычисления суммарной дальнодействующей составляющей силы.

1.3 Макрочастицы в моделях частица-сетка

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

na(t, r, p) = nai(t, r, p). (1.6)

i

Каждая функция nai представляет собой функцию распределения плотности макрочастицы, совокупности из NPai физических частиц, близких друг к другу в фазовом пространстве. В общем случае, функции na. представляются в виде

nai(t, r, p) = NpatSr(r, r«i(t))Sp(p, pai(t)). (1.7)

Здесь Sr(r, ra.(t)) и SP(p, pai(t)) — функции распределения распределения плотности в макрочастице по координатам и импульсам (функции формы), а ra. (t) и pai (t) — координаты, определяющие положения макрочастиц в фазовом пространстве. Функции формы симметричные, положительно определенные и для них выполняется условие нормировки

J S(х,у)йх = 1. (1.8)

х

В стандартных методах частиц функция Sp(p, р^ (г)) берется в виде дельта-функции Дирака

Sp(p, Раг (г)) = 5(р - Ра (г)). (1.9)

Такой выбор означает, что в состав макрочастицы входят физические частицы с одинаковой скорость. Выбор функции Sr(г, га.(г)) также в виде дельта-функции приведет к модели частиц нулевого размера. Однако, такой выбор не используется на практике, поскольку он порождает возникновение значительных нефизических шумов в решении и флуктуаций плотности.

Подставив уравнения (1.6)-(1.9) в уравнение Власова и рассмотрев уравнения сохранения его первых трех моментов, можно получить уравнения

=0, (1.10)

dt

drai pai

dt y а^ mo

dpai

(1.11)

dt

= qa (Ea + [vai X Bai]) , (1.12)

Eai = J E(r,t)Sr(r, rai)dr, (1.13)

Bai =f B(r,t)Sr(r, rai)dr. (1.14)

Уравнение (1.10) означает постоянство во времени числа физических частиц, входящих в состав макрочастицы. Уравнения (1.11)-(1.14) представляют собой уравнения движения макрочастицы. Подставив уравнения (1.6), (1.7), (1.9) в уравнения (1.4), (1.5), получим выражения для вычисления плотности пространственного заряда и плотности тока с использованием макрочастиц

а

p(t' r) = S S qa-i Sr (r' rai (t))

a i

(1.15)

r) = EE Sr (r, r«i (t)) = E r)- (!-16)

a i a

Здесь qa. = qaNPa. — заряд г-й макрочастицы в потоке с номером a, ja — плотность тока отдельно взятого потока с номером a.

1.4 Постановка задачи моделирования потоков заряженных частиц в электростатическом приближении

С учетом вышесказанного, рассмотрим постановку задачи расчета самосогласованных стационарных электромагнитных полей и динамики стационарных потоков макрочастиц.

Рассмотрим расчетную область Q = Q U Г, где Г = Г U Г2 — граница расчетной области. Предположение о том, что электрическое поле в рассматриваемой системе является безвихревым, позволяет свести систему уравнений Максвелла (1.3) к уравнениям Ампера и Пуассона для вычисления индукции магнитного поля B, электрического потенциала U и напряженности электрического поля E.

Список литературы диссертационного исследования кандидат наук Алцыбеев Владислав Владимирович, 2016 год

3.1.2 Варианты использования

Комплекс прикладных программ, предназначенный для моделирования динамики потоков заряженных частиц должен удовлетворять следующим требованиям:

• Поддерживать возможность использования различных систем координат при расчетах, таких как декартовые двумерная и трехмерная, цилиндрическая с учетом осевой симметрии и цилиндрическая трехмерная, полярная. Данная возможность позволяет получать максимально точные решения прикладных задач с учетом особенностей их геометрии. Например, в цилиндрической геометрии выбор полярных или цилиндрических координат, структурированные сетки в которых более точно аппроксимируют границы расчетной области, позволяет сократить объем вычислений и приводит к более точным и устойчивым решениям, чем при использовании декартовых координат.

• Возможность выбора пользователем метода решения задачи. В зависимости от наличия или отсутствия стационарного решения рассматриваемой задачи, пользователю стоит использовать метод частиц в ячейках или итерационный метод. В некоторых же задачах, таких как моделирование динамики пучка в линейном резонансном ускорителе, заранее допустимо использование только метода частиц в ячейках.

• Возможность точного описания пользователем параметров моделируемого устройства: сложной геометрии расчетной области, граничных условий на электродах, начальных данных для потоков частиц.

• Возможность настройки параметров используемых численных методов: генерация неоднородной расчетной сетки с её визуализацией, настройка шага интегрирования уравнений движения, точности решения уравнений поля и т. д.

• Визуализация результатов расчета с возможностью сохранения результатов и вывода в файл.

• Загрузка и сохранение проектов — файлов, содержащих все введенные пользователем параметры задачи, геометрию и некоторые результаты расчетов (структуру расчетной сетки, положения частиц, распределения полей).

• Возможность расширения комплекса программ на другие задачи, например двумерная задача с экструзией сетки.

UML (unified modeling language) диаграмма вариантов использования разработанного комплекса программ, отвечающего перечисленным выше требованиям представлена на рис. 3.2.

Рисунок 3.2: иМЬ диаграмма вариантов использования

3.2 Общее описание архитектуры

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

Одной из основных задач при разработке иерархии классов сервера была минимизация использования виртуальных функций для увеличения производительности и легкости дальнейшего изменения кода с целью использования технологии CUDA (Compute Unified Device Architecture). Технология виртуальных функций и полиморфизма времени выполнения была эффективно заменена шаблонами С++ и полиморфизмом времени компиляции. В итоге, виртуальные функции используются только в некритичной части взаимодействия клиента и сервера.

При описании архитектуры комплекса, под моделью (Model) будем понимать объединение параметров моделируемого устройства (Device) и набора применяемых численных методов с параметрами (Solver).

В зависимости от типа задачи и типа применяемого метода, класс DeviceStatus, описывающий состояние моделируемого устройства, будет агрегировать в себе данные различных типов и являться шаблонным. Такая структура позволяет использовать единую реализацию численного метода на высоком уровне для всех задач (общая схема метода) и индивидуальную реализацию методов для задачи каждого типа на низком уровне (решение уравнений поля, расчет траекторий частиц и т. д.).

Диаграмма классов объекты которых составляют класс DeviceStatus представлена на рис. 3.3.

I dataType j

ParticlesBase

Н>

-particlesNumber : int

-cellsNumbers : int

-charges : std::vector: dataType

dataType

Particles2dpolar

coordinates : std::vector: dataType momentums : std::vector: dataType

dataType

coordinates : std::vector: momentums : std::vector:

dataType dataType

dataType

particles2d

coordinates : std::vector: momentums : std::vector:

dataType dataType

coordinates : std::vector: momentums : std::vector:

dataType dataType

dataType

ParticlesSource2d

currentDensity : std::vector: double sourceSurface : BoundaryContainer2d

о

dataType

ParticlesSource3d

currentDensity : std::vector: double sourceSurface : BoundaryContainer3d

dataType

EmitterDeviceBase

-¿y

dataType

EmitterDevice2d

particlesSource : ParticlesSource2d

+generateParticles(p : particles2d *) +generateParticles(p : Particles2dpolar *) +generateParticles(p : Particles3d *) +generateParticles(p : Particles3dcyl *)

dataType

l<>

EmitterDevice2daxs

particlesSource : ParticlesSource2d

+generateParticles(p : particles2daxs *)

dataType

EmitterDevice3d

particlesSource : ParticlesSource3d

+generateParticles(p : Particles3d *) +generateParticles(p : Particles3dcyl *)

particlesDataType

emitterDeviceType

dataType

ParticlesFlow

dynamicsData : std::vector<particlesDataType mass : double charge : double

emitterDeviceStatus : emitterDeviceType boundaryConditions : BoundaryConditions

generateParticles()

dataType

GridData2d

-nodesCoordinates : std::vector: dataType -fieldComponents : std::vector: dataType -chargeDensity : std::vector: dataType

dataType

GridData2daxs

-nodesCoordinates : std::vector: dataType -fieldComponents : std::vector: dataType -chargeDensity : std::vector: dataType

dataType

GridData2dpolar

-nodesCoordinates : std::vector: dataType -fieldComponents : std::vector: dataType -chargeDensity : std::vector: dataType

dataType

GridData3d

nodesCoordinates : std::vector: dataType fieldComponents : std::vector: dataType chargeDensity : std::vector: dataType

k>

г-----------

particlesDataType gridDataType emitterDeviceType boundaryContainerType meshDataType dataType

DeviceStatus

particlesFlowsStatus : std::vector: particlesFlow domainBoundaries : boundaryContainerType mesh : meshDataType gridData : gridDataType boundaryConditions : BoundaryConditions electrodes : std::vector: Electrode

BoundaryConditions

-boundariesList:: std::vector : int -conditionType : string -conditionProperties : std::vector: double

dataType

Point

-x : dataType

-y : dataType

-z : dataType

dataType Edge

pointl : Point point2 : Point

dataType

Triangle

-pointl : Point -point2 : Point -point3 : Point

dataType

BoundaryContainer2d

geomData : std::vector: Edge

dataType

BoundaryContainer3d

geomData : std::vector: Triangle

— MeshContainer2d

-meshData : std::vector: Point -cellsPoints : std::vector: int +convertToGridData2d() +convertToGridData2daxs() +convertToGridData2dpolar()

MeshContainer3d

meshData : std::vector: Point cellsPoints : std::vector: int

+convertToGri d D ata3d()

electrodeEdges : std::vector: Edges currenrAtEdges : std::vector: double

Рисунок 3.3: UML диаграмма классов части device

Классы с именами Particles и GridData представляют собой структуры для хранения данных о частицах и значениях полей и плотностей в узлах эйлеровых сеток, а также содержат некоторые функции обработки данных с учетом особенностей используемой системы координат. Классы EmitterDevice содержат дан-

ные об эмиттере и свойствах начального распределения частиц, а также имеют методы генерации начальных данных. При этом класс ParticlesFlow агрегирует общую информацию о каждом потоке частиц. Классы MeshContainer2d и MeshContainer3d содержат более детальные описания структур эйлеровых сеток без значений полей, функции экспорта в формат VTK для визуализации, а также методы конвертирования в форматы GridData. Классы Point, Edge, Triangle, BoundaryContainer предназначены для работы с геометрией. Для описания граничных условий на электродах и для потоков заряженных частиц используется класс BoundaryConditions. Класс Electrode предназначен для расчета плотности тока, протекающего по электродам моделируемого устройства и плотности их облучения.

Дополнительно, все крупные структуры данных являются шаблонными по использумому типу числового представления dataType для поддержания возможности проводить расчеты с плавающей точкой двойной или одинарной точности.

Данная архитектура позволяет удобным образом, правильно скомбинировав параметры шаблона класса DeviceStatus, создавать объект нужного типа для решения той или иной задачи. Например, для решения двумерной задачи в декартовых координатах с использованием типа dataType можно ввести тип device2d.

typedef DeviceStatus<Particles2d, GridData2d, EmitterDevice2d, BoundaryContainer2d, MeshContainer2d, dataType> device2d

По аналогии введем типы device2d, device2dpolar, device2daxs, device3d.

Диаграмма классов SolverPIC и SolverPTI, описывающих применяемые численные методы, представлена на рис. 3.4. Класс SolverPIC реализует метод частиц в ячейках, а класс SolverPTI — итерационный метод. Класс PaticleGridInterface реализует функции интерполяции сил, раздачи заряда

[dataTypel

Рисунок 3.4: UML диаграмма классов части device

по узлам сетки, проверки прохождения частицами границ области и т. д. Класс ParticlesMover содержит функции интегрирования уравнений движения макрочастиц. Класс MeshGenerator предназначен генерирования расчетной сетки, FieldSolver — для расчета распределений электромагнитного поля. Классы EmissionCurrentSolverPTI и EmissionCurrentSolverPIC реализуют модели и алгоритмы нахождения тока эмиссии в режиме ограничения тока пространственным зарядом. Функции-члены данных классов могут являться перегруженными (например функции UpdateMomentums() и UpdatePositions()) или иметь единственную реализацию для нескольких типов входных данных.

На рис. 3.5 представлена общая UML диаграмма основных классов серверной части с учетом описанных выше составных частей DeviceStatus и Solver. Класс ModelInterface является абстрактным и предоставляет клиентской части интерфейс доступа к данным сервера. Его наследник, класс ModelBase, содержит реализацию большинства виртуальных функций интерфейса ModelInterface.

particlesDatatype

gridDataType

emitterDeviceType

boundaryContainerType

meshDataType

dataType

Modelinterfaoe

TV

DeviceStatus

->

device2d, solverPIC

Model2dPIC

device2d, solverPTI

ModelBase

deviceType | solvertype j dataType

Й-

dataType

->

-deviceStatus : deviceType -solver : SolverType

+Simulate() +FieldSimulate()

£-Л-7Х

SolverPIC

dataType

<b

SolverPTI

device2daxs, solverPIC

Model2daxsPIC

device2d, solverPTI

Model2dPTI

Model2daxsPTI

device2dpolar, solverPIC

Model2dpolarPIC

device2dpolar, solverPTI

device3d, solverPIC

Model3dPIC

device3d, solverPTI

Model2dpolarPTI

Model3dPTI

Рисунок 3.5: Общая UML диаграмма классов

Он агрегирует в себе объекты типа DeviceStatus и Solver и является шаблонным.

Таким образом, вызов общей функции Simulate() класса ModelBase будет приводить к использованию необходимых методов над выбранными структурами данных, определяемых параметрами шаблонов класса DeviceStatus.

Класс ModelBase имеет несколько наследников, в зависимости от типа задачи и типа решателя, определяющих спецификацию шаблона:

Model2dPIC — двумерная задача с использованием метода частиц в ячейках;

Model2dPTI — двумерная задача с использованием итерационного метода;

Model2daxsPIC — задача в цилиндрических координатах с учетом осевой симметрии с использованием метода частиц в ячейках;

Model2daxsPTI — задача в цилиндрических координатах с учетом осевой симметрии с использованием итерационного метода;

Model2dpolarPIC — двумерная задача в полярных координатах с использованием метода частиц в ячейках;

Model2dpolarPTI — двумерная задача в полярных координатах с использованием итерационного метода;

Model3dPIC — трехмерная задача с использованием метода частиц в ячейках;

Model3dPIC — трехмерная задача с использованием итерационного метода.

При создании пользователем нового проекта, в зависимости от его выбора используемых системы координат и метода, создается объект одного из данных классов-наследников, аттрибуты которого затем меняются пользователем. Процесс сохранения проекта заключается в сериализации всех аттрибутов класса ModelBase с использованием библиотеки Boost.

При загрузке существующего проекта, в зависимости от его типа, также создается объект одного из наследников ModelBase, в который затем происходит десериализация данных из файла проекта.

3.3 Алгоритмы хранения, обработки и удаления частиц

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

3.3.1 Структура данных хранения частиц

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

над крупным массивом данных выполнялось множество однотипных операций. Классическим в этом плане является вопрос о выборе между массивом объектов, или одном объекте, содержащем массивы или указатели на массивы, или контейнеры ("structure of array vs array of structure"). Согласно многим источникам, например [97, 98], необходимо делать выбор в пользу одного объекта, содержащего массивы. Таким образом, в разработанном комплексе программ используется следующая структура данных для хранения положений, импульсов и зарядов макрочастиц (пример для 2d задач):

template <class DataType>

class Particles2d {

std::vector<DataType> x, y, px, py, q;

std::vector<int> cellNumber; }

Кроме того, для каждой частицы запоминается последний номер ячейки, в которой она находилась, что позволяет значительно сократить время выполнения процедуры локализации частицы на сетке.

С точки зрения быстроты выполнения кода в каждом отдельном пункте алгоритмов частица-сетка (рис. 1.3, 1.4) выгодно обрабатывать макрочастицы блоками, записывать результаты во временные переменные и затем использовать эти результаты на следующем этапе. Данный подход позволяет повысить эффективность использования кэш-памяти, при этом объем дополнительной памяти будет пропорционален размеру блока. Аналогичный подход используется, например, в алгоритмах перемножения матриц [99] (алгоритм Фокса). Проводился тест зависимости скорости расчета от размера блока. Использовался процессор AMD FX8320, один вычислительный поток, 367000 макрочастиц рассчитывались в течение 7 нс. Результаты расчета показывают, что производительность кода будет иметь некоторую зависимость от размера блока — существует некоторое оптимальное значение, обеспечивающее наилучшее использование кэш-

памяти и скорость расчета (рис. 3.6). Данное оптимальное значение, по всей видимости, будет зависеть от конкретных характеристик задачи и архитектуры вычислительной системы.

Также, важным аспектом, влияющим на производительность, является алгоритм удаления частиц, покинувших расчетную область. Простая реализация, подразумевающая удаление частиц на каждом шаге будет приводить к необходимости часто дефрагментиро-вать вектора, хранящие информацию о частицах, что приведет к существенному замедлению расчета. Для минимизации этих затрат, на места тех частиц, что должны быть удалены, перезаписываются новые, помещаемые на каждом шаге в расчетную область. Если число частиц, которые необходимо удалить, превышает число новых частиц, остаток помечается специальным образом и не участвует в дальнейших расчетах. Раз в 20 шагов все подлежащие удалению частицы удаляются и вектора дефрагментируются.

3.3.2 Параллельная реализация методов частица-сетка с использованием OpenMP

Методы частица-сетка успешно параллелятся как для систем с общей памятью, так и для распределенных систем, или GPU, и обладают хорошей масштабируемостью.

Существует два основных подхода к распараллеливанию методов частиц в ячейках, описанные в литературе [37,100-102].

Рисунок 3.6: Зависимость производительности кода от размера блока

Первый основан на декомпозиции расчетной области на подобласти между вычислителями (потоками или процессами). Каждый вычислитель при этом ответственен за обработку тех частиц, которые находятся в его подобласти. В таком случае, при операции расчета пространственного заряда не будет возникать конфликтов чтения, однако нагрузка на вычислители при этом может быть не сбалансирована (например, когда большая часть частиц локализована в малой области). Дополнительно, вычислительный алгоритм в таком случае значительно усложняется и требует дополнительных балансировок нагрузки. Такая стратегия, однако, подходит в тех случаях, когда объем памяти, который возможно выделить для хранения сеточных значений пространственного заряда ограничен, например, при распараллеливания для GPU.

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

Для реализации в разработанном комплексе программ был выбран второй подход, как наиболее подходящий для систем с общей памятью и позволяющий добиться наилучшей масштабируемости.

3.2.1 Исследование масштабируемости метода частиц в ячейках

Реализованную последовательность вычислений в методе частиц в ячейках с учетом разбиения частиц на потоки и на блоки внутри потока можно представить с помощью следующего алгоритма:

Зарезервировать память для хранения временных переменных Выделить для каждого потока локальный массив для хранения плотности заряда

Распределеить начальные данные между потоками

#pragma omp parallel num_threads(numThreads) {

while время меньше заданного {

#pragma omp barrier

#pragma omp single {

Просуммировать заряд, накопленный во всех потоках Решить уравнения поля

Вычислить ток эмиссии }

thread = omp_get_thread_num();

Разбить частицы, обрабатываемые потоком thread на блоки

for все блоки частиц {

Вычислить весовые коэффициенты для всех частиц в блоке Вычислить поле для всех частиц в блоке Сохранить текущие положения частиц в блоке

Вычислить новые положения и импульсы для всех частиц в блоке Проверить взамодействие частиц блока с границами расчетной области Локализовать частицы в блоке на сетке

Вычислить весовые коэффициенты для всех частиц в блоке

Вычислить плотность заряда }

Сгенерировать новые частицы

Удалить частицы, покинувшие расчетную область }

}

Таким образом, на каждом шаге интегрирования в методе частиц в ячейках необходимо проводить минимум две синхронизации потоков выполнения — перед расчетом поля и после него (неявная синхронизация потоков присутствует также после выполнения блока omp single), а вычислительные потоки порождаются только один раз — на входе в процедуру расчета. Расчет поля производится с использованием только одного потока. На рис. 3.7 представлены результаты исследования масштабируемости реализованного алгоритма в зависимости от числа макрочастиц. Для расчетов использовался восьмиядерный процессор Intel Xeon X5675. При увеличении числа макрочастиц, алгоритм достигает значительных показателей ускорения (более 6 раз), чего удалось добиться за счет хорошей балансировки нагрузки между потоками и минимизации времени выполнения последовательной части программы.

Рисунок 3.7: Зависимость ускорения от числа вычислительных потоков и числа макрочастиц. Здесь Ыр — число макрочастиц

3.2.2 Исследование масштабируемости итерационного метода

Реализованную последовательность вычислений в итерационном методе с учетом разбиения частиц на потоки и на блоки внутри потока можно представить с помощью следующего алгоритма:

Зарезервировать память для хранения временных переменных Выделить для каждого потока локальный массив для хранения плотности заряда

Распределеить начальные данные между потоками

#pragma omp parallel num_threads(numThreads) {

while достигнута сходимость {

#pragma omp barrier

#pragma omp single {

Просуммировать заряд, накопленный во всех потоках Решить уравнения поля Вычислить ток эмиссии

Проверить условие сходимости }

#pragma omp barrier

thread = omp_get_thread_num();

Сгенерировать новые частицы

Разбить частицы, обрабатываемые потоком thread на блоки

for все блоки частиц {

Расчитать динамику частиц в блоке до выхода их из расчетной области,

накапливая пространственный заряд в узлах сетки }

} }

На рис. 3.8 представлены результаты исследования масштабируемости реализованного алгоритма в зависимости от числа макрочастиц. Для расчетов использовался восьмиядерный процессор Intel Xeon X5675. При увеличении числа макрочастиц, алгоритм достигает значительных показателей ускорения (более 7 раз).

Ускорение

Рисунок 3.8: Зависимость ускорения от числа вычислительных потоков и числа

макрочастиц. Здесь N — число макрочастиц

Глава 4

Результаты моделирования эмиссионных устройств

4.1 Цилиндрический диод

В качестве теста рассмотрим задачу моделирования цилиндрического диода с эмиттером конечной длины Lem — 0.012 m. Ионы H + также могут эмиттиро-ваться с поверхности длиной Lem — 0.012 m. Диод имеет длины L — 0.06 тс Rc — 0.0055 т радиусом катода and Ra — 0.0005 m радиусом анода. Напряжение Ua приложено к аноду заземленным катодом.

4.1.1 Случай цилиндрических координат

т-ч VJ VJ

В силу наличия аксиальной и центральной симметрии, можно рассмотреть двумерную задачу в RZ координатах на половине диода (рис. 4.2).

1.1.1 Электронный монопоток

Будем исследовать сходимость предложенного подхода вычисления плотности тока в итерационном методе (модель эмиссии Гауса) с помощью сравнения полученного решения с аналитическим решением Ленгмюра (2.24) для цилиндрического диода. С уменьшение шага расчетной сетки, значение плотности тока эмиссии в центре диода должно приближаться к аналитическому решению при достаточно невысоких (нерелятивистских) значениях анодного напряжения. Была проведена серия расчетов с шагами расчетной сетки вдоль оси R hr — 0.0002 m, hr — 0.0001 m, hr — 0.00005 m, hr — 0.000025 m, hr — 0.00002 m. Шаг сетки вдоль оси Z использовался равным hz — 0.0001 m для всех расчетов. Параметры эмиттирующих ячеек Lem — hz and Hem — hr для всех расчетов, число траекторий макрочастиц 800, параметр выбора шага интегрирования (коэффициент прохождения ячейки сетки) по времени K — 2.

Z( м

0.03

0.025

0.02

0.015

0.01

0.005 "

Анод U = U

Поверхность

эмиссии

ионов

L /2

em

Катод U = 0

L/2

Поверхность

эмиссии

электронов

L /2

em

R' м , Л -3 X 10

Рисунок 4.2: Модель цилиндрического диода. Граничные условия Неймана используются на пунктирных линиях

Для всех расчетов итерационный метод показал хорошую сходимость (см рис. 4.3 а, значение тока эмиссии достигалось примерно за 20 итераций). На рис. 4.3 б относительное отклонение Я = \(Зо - Зь)/Зь\ (здесь ,1ь — аналитическое решение Ленгмюра и — численное значение плотнсоти тока эмиссии, полученное в центре диода) представлено в зависимости от значения 1/Нг. Распределения плотностей токов эмисиии для различных значений кг с представлены на рис. 4.4. На рис. 4.5 представлены траектории электронов. Из полученных результа-

тов можно сделать вывод, что предлагаемый подход обладает хорошей сходимостью и согласованностью с аналитическим решением.

1.1.2 Биполярный поток с наличием отраженных электронов

Ионы Н+ эмиттируются с анода, при этом, некоторая часть падающих на анод электронов может отражаться назад, в зазор анод-катод, потеряв при этом часть своей энергии. Этот процесс повторяется до тех пор, пока электрон не поглотится. Таким образом, пространственный заряд в прианодной области растет, что приводит к росту электрического поля на аноде и росту тока эмиссии ионов, что, в свою очередь увеличивает ток эмиссии электронов.

a

К

c

К

а

0

0

2

4

а б

Рисунок 4.3: Исследование использования модели эмиссии Гаусса в итерационном методе для цилиндрического диода. Сходимость итераций метода для Нг = 0.00002 м (а) и сходимость численного решения к аналитическому при использовании различных шагов сетки (б)

Рисунок 4.4: Распределния плотности тока эмиссии вдоль половины эмиттера

для различных шагов расчетной сетки

Следуя работе [32] параметризуем отраженные электроны с помощью коэффициента а, jscat = Щ{иЫс1епь и в, Е3а = вЕпЫсе**. В общем случае, коэффициенты а и в будут зависеть от материала мишени, её диаметра и энергии

Щ ы

0.035 0.03 0.025 0.02 0.015 0.01 0.005 0

0.005 0.01 0.015 0.02 0.025 0.03 0.035 2, м

Рисунок 4.5: Траектории частиц в цилиндрическом диоде

падающих элетронов. Проводился ряд расчетов с различными коэффициентами а и в для 2 МВ напряжения на аноде с использованием шага расчетной сетки Нг = 0.0001 т. На рис. 4.6 представлено отношение I/1с1 (здесь 1с1 — ток эмиссии электронов, полученный с помощью решения Чайлда (2.23) с й = гс — га) в зависимости от значений в для различных значений а. Мы можем сравнить численное решение с зависимостями, полученными в работе [32] с использованием одномерного аналитического решения. Результаты достаточно точно согласуются, за исключением того, что в аналитическом решении существует диапазон значений а и в, при которых не существует стационарного решения. Этот факт может быть объяснен тем, что в двумерном случае частицы могут двигаться вдоль оси 2, таким образом электрическое поле может уменьшаться.

1,5

1

0.75

!

0,5

а

отношение 1/1 с|

-^-а=0.2

. ---а=0.4

-*-а=0.5

----- - ~~ ** —V-V—

0.2

0.4

0.6

°.8 р 1

б

1.5

а)

О а=0.2 □ а=0.4 х а=0.5

ё 1-0

-О

о 0.5

а)

ш

0.0

0.0 0.2 0.4 0.6 0.8

и

1.0

Рисунок 4.6: Отношение тока эмиссии к решению Чайлда для биполярного потока и наличия обратного рассеяния электонов. Численное решение, полученное в двумерном осе симметричном случае (а) и аналитическое решение (источник — работа [32]) (б)

0

0

4.1.2 Случай декартовых координат

Рассмотрим задачу моделирования цилиндрического диода в случае двумерной декартовых координат. В данном разделе будем считать диод бесконечно длинным. Данный пример хорошо иллюстрирует работу алгоритма определения тока эмиссии, основанного на применении закона Гаусса, поскольку поверхность эмиссии является криволинейной, а аналитическое решение Ленгмюра известно. Была проведена серия расчетов с различными шагами расчетной сетки Нх = Ну. Длина и ширина эмиттирующих ячеек выбрана следующей Ьет = 0.8 мм для эмиттера элеткронов и Ьет = 0.05 мм мм для эмиттера ионов, Нет = Нх. Число траекторий макрочастиц 5000 для электронного и ионного потоков, параметр выбора шага интегрирования (коэффициент прохождения ячейки сетки) по времени К = 2. Выбор достаточно длинных эмиттирующих ячеек позволяет дополнительно сгладить распределение плотности тока эмиссии.

1.2.1 Результаты, полученные с использование модели эмиссии Гаусса

В случае электронного монопотока предлагаемый метод показал хорошую сходимость. Численное решение приближается к аналитическому при сгущении эйлеровой сетки. На рис. 4.7 (а) относительное отклонение Я = |(3о — 3ь)/3ь\ (здесь 3 — аналитическое решение Ленгмюра и — численное значение плотности тока эмиссии) представлено в зависимости от значения 1/Нх в случае электронного монопотока. В случае биполярного потока метод также сходится. На рис. 4.7 (б) представлены полученные численно распределения плотности тока эмиссии в случае монопотока и биполярного потока, а также аналитическое решение Ленгмюра. Аналитическое и численное решения достаточно точно согласуются. На рис. 4.8 представлено распределение модуля напряженности электрического поля в расчетной области диода.

а

Относительное отклонение, К

б

2.2 2 1.8 1.6 1.4 1.2 1

0.8 0.6 0.4 0.2 0

х 105 Плотность тока, А/м

-Численное решение в случае монопотока

---Численное решение в случае бипотока

X Аналитическое решение Ленгмюра

-*хххххххх-х-ххх-хх-хххххх

0.005 0.01 0.015 0.02 0.025 0.03 0.035

Длина эмиттера, м

Рисунок 4.7: Исследование использования модели эмиссии Гаусса в итерационном методе для цилиндрического диода в декартовой геометри.

Сходимость численного решения к аналитическому при использовании различных шагов сетки (а) и распределния плотностей тока эмиссии вдоль

эмиттера (б)

о

X

1.2.2 Результаты, полученные с использование модели эмиссии Чайлда-Ленгмюра

Дополнительно проводи-

лись попытки получить решения рассматриваемой задачи с помощью модели эмиссии 4.81е+006 Чайлда-Ленгмюра. В случае элек-

,01е+006

тронного монопотока модель

3.21е+006

удалось удачно применить, однако в случае биполярного потока, 2.52е+ооз метод разошелся (см. рис. 4.9). Таким образом, предлагаемый в диссертации метод показал

Рисунок 4.8: Распределение напряженности

преимущество над стандартным

электрического поля в цилиндрическом

подходом в данной задаче.

диоде в декартовой геометрии

а

12000 10000 8000 6000 4000 2000 0

Ток эмиссии, А/м

0

10 15

Номер итерации

б

х 104 Ток эмиссии, А/м

10 15 20 25 Номер итерации

Рисунок 4.9: Исследование сходимости метода нахождения тока эмиссии, основанного на модели Чайлда-Ленгмюра в случае биполярного потока в цилиндрическом диоде в декартовой геометрии. Сходимость в случае электронного монопотока (а) и в случае биполярного потока (б)

5

0

5

4.2 Диод с эллиптическим эмиттером

В качестве второй тестовой задачи будем рассматривать задачу моделирования динамики пучка в диоде с эмиттером в виде полуэллипсоида. Осевое сечение геометрии диода представлено на рис. 4.10. Малая и большая полуоси эллипсоида Ь = 0.01 м, а = 0.02 м, радиус диода ( = 0.15 м, расстояние Ь = 0.07 м, высота цилиндрического основания эмиттера К = 0.02 м. Напряжение иа = 10 000 В приложено к аноду с заземленным катодом. В силу наличия аксиальной и центральной симметрии, задачу можно решать в двумерной геометрии в ЯЕ координатах на половине диода. Расчетная модель задачи представлена на Fig. 4.11.

Ось

симметрии

с1

Рисунок 4.10: Геометрия задачи с эллиптическим эмиттером

Была проведена серия расчетов с шагами расчетной сетки Кг = кг = 0.001 м, Кг = кг = 0.0005 м, Кг = кг = 0.00025 м. Параметры эмиттирующих ячеек Ьет = 0.001 т и Нет = Кг для всех расчетов,число траекторий макрочастиц 1200, шаг интегрирования по времени выбирался из условия прохождения макрочастицей не более чем половины ячейки расчетной сетки за один шаг. Распределения плотностей токов эмиссии для различных значений Кг, полученные

с использованием модели эмиссии, основанной на законе Гаусса представлены рис. 4.12 (а). Также представлены траектории электронов на рис. 4.13.

Также для сравнения проводились расчеты с использование модели эмиссии Чайлда-Ленгмюра. Для расчета плотности тока использовались выражения (2.24) или (2.25) для различных участков эмиттера (выбор выражения зависит от формы участка). Расстояние между катодом и анодом каждого виртуального диода выбиралось равным кг. Распределение плотности тока эмиссии для различных

Рисунок 4.11: Расчетная модель задачи с эллиптическим эмиттером

значений Нг полученные с использованием модели эмиссии Чайлда-Ленгмюра представлены на рис. 4.12 (б). В результате сравнения рис. 4.12 (а) и (б) можно заключить, что решение, полученное с помощью разработанного метода обладает большей степенью гладкости, что можно трактовать, как снижение уровня счетного шума вызванного дискретизацией задачи и погрешностью методов.

аб

Рисунок 4.12: Распределения плотности тока эмиссии вдоль поверхности эмиттера, полученные с помощью модели эмиссии, основанной на законе Гаусса (а) и Чайлда-Ленгмюра (б) для различных шагов расчетной сетки

Рисунок 4.13: Траектории частиц в плоскости RZ в диоде с эллиптическим

эмиттером

4.3 Источник радиально сходящегося пучка электронов триодного типа

4.3.1 Параметры источника

В работе [28] предлагается конструкция источника электронов триодного типа, генерирующего радиально сходящийся пучок электронов с целью обработки поверхностей цилиндрических мишеней (например, оболочек тепловыделяющих элементов реакторов).

Источник состоит из цилиндрических катода, анода и ускоряющей сетки с внешним эмиттером электронов (рис. 4.14, 4.15). Ускоряющая сетка состоит из 108 проволок диаметров 0.2 мм. Значения основных параметров источника приведены в табл. 4.1.

Таблица 4.1: Параметры источника

Длина катода, м Ь = 0.98

Радиус катода, м тс = 0.15

Радиус ускоряющей сетки, м тд = 0.1

Радиус анода/мишени, м Га = 0.005

Напряжение катод-сетка, кВ ид = 20 - 35

Напряжение катод-анод, кВ иа = 120

0.2 0.15 0.1 0.05 0

-0.05 -0.1 -0.15

-0-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2

X, m