Поэлементная технология расчета интенсивных пучков заряженных частиц Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 9, № 3, 2004

ПОЭЛЕМЕНТНАЯ ТЕХНОЛОГИЯ РАСЧЕТА ИНТЕНСИВНЫХ ПУЧКОВ ЗАРЯЖЕННЫХ

ЧАСТИЦ*

В. М. Свешников Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: victor@lapasrv.sscc.ru

The element-oriented technology for calculation of intensive charged particle beams is suggested. Algorithms for the integration of the equations of the charged particle motion, calculation of the electric field intensity and calculation of the space charge, which are developed by virtue of the element-oriented technology are described. The algorithms have the second order of accuracy. Numerical solution of a test problem is presented.

Введение

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

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

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

*Работа поддержана грантом "Университеты России" (№ УР.03.01.017).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.

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

— интегрирование уравнений движения заряженных частиц;

— расчет напряженности электрического поля;

— распределение объемного заряда.

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

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

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

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

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

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

Решение исходной задачи в общем случае ищется на квазиструктурированных сет-

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

В настоящей работе излагаются алгоритмы интегрирования уравнений движения заряженных частиц, расчета напряженности электрического поля и распределения объемного заряда. Предполагается, что потенциал электрического поля вычислен в узлах сетки П. Алгоритмы расчета потенциала приведены в [1], а алгоритмы решения самосогласованной задачи в целом изложены, например, в [3, 5].

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

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

Ду = -

Р_

£о

с граничными условиями

Iу = 0,

уравнения движения заряженных частиц

¿г

V,

¿V

т

ПЕ, Е = — grad у

с начальными условиями

г , V

4=0

V

4=0

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

div j = 0, j = pv.

(1) (2)

(3)

(4)

(5)

Здесь р — плотность объемного заряда; е0 — диэлектрическая проницаемость вакуума; I — оператор граничных условий; г = (х,у) — радиус-вектор заряженной частицы, где х,у — декартовы или цилиндрические координаты; V — ее скорость; п — отношение заряда к массе; Ь — время; Е — напряженность электрического поля; j — плотность тока. Область О предполагается двумерной (плоской или осесимметричной). В осесимметричных областях расчет проводится в меридиональной плоскости, причем х означает радиальную, а у — осевую координаты.

Задача (1)-(5) аппроксимируется на сетке П, покрывающей область О. В общем случае П может быть непрямоугольной и неструктурированной сеткой. Одним из примеров непрямоугольных сеток являются квазиструктурированные сетки, предложенные в работе [4]. Их суть заключается в том, что расчетная область О покрывается набором N<0 подобластей О к, к = 1, 2,... , N0, так, что О = У О к. В каждой О к строится структурированная

к= 1

ма

сетка Пк, а результирующая квазиструктурированная сетка есть П = У Пк. Отметим, что

к=1

в Ок возможно введение локальной системы координат, удобной для проведения расчетов

г

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

2. Интегрирование уравнений движения

Уравнения движения (3) с начальными условиями (4) будем интегрировать численно с шагом тп по следующей неявной схеме второго порядка точности:

r«+1 _ r«. + V«

vn+l _ v«

2

E«+1 + E™

(6)

v n-2-' (7)

' n ^

где индекс n = 0,1,... означает номер временного шага.

Данная схема обладает двумя важными физическими свойствами. Во-первых, она обратима по времени. Это означает, что, заменив в ней т« на _ т« и предполагая известной (n + 1)-ю точку, мы получим исходную n-ю точку, как это непосредственно следует из формул (6) и (7). Во-вторых, рассматриваемая схема аппроксимирует закон сохранения энергии, который имеет вид

«+112 i « |2 _ о ( .«+1

|^+1|2 _|^|2 = - <^п), (8)

п. Здесь И2 = (^)2 + « бы это показать, запишем (6), (7) покомпонентно и умножим уравнения (6) на уравнения (7) соответственно для компонент х и у. В результате получим соотношения

с точностью о(т3), где т = тахтп. Здесь |уп|2 = (^П)2 + (^П)2, п = ^(хп,уп). Для того что-

п

(v«+1)2 _ (v«)2 = n(x«+1 _ x«)(E«+1 + E«),

(v«+1)2 — (v«)2 = n(y«+1 _ y«)(E«+1 + E«).

Так как

(9)

(E«+1 + E«) = e«+1/2 + o(t2), r«+1 _ r« = 0(т),

где E«+1/2 — значение напряженности в средней точке,

n+1 n

r«+1/2 = £-, (10)

2

равенства (9) можно записать в виде

(v«+1)2 _ (v«)2 = 2n(x«+1 _ x«)E«+1/2 + 0(т3), (v«+1)2 _ (v«)2 = 2n(y«+1 _ y«)E«+1/2 + 0(т3).

n

Здесь и в дальнейшем 0(тр) означает вектор с компонентами 0(тр). Выразим напряженность в средней точке при помощи соотношений

Еп+1/2

ЕГ =

_ (Еп'° + Б"'1) = 2 ' (^(жп+1'Уп+) - <^(жп

Еп,г = ) _ ^(хп+\уп) . = 0 1

2 уП+1 _ уП ' '

Подстановка данных выражений в равенства (11), а затем сложение последних дают закон сохранения (8) с точностью 0(т3).

Исключая из уравнений (6), (7) скорость уп+1, получим уравнение для нахождения координат

гп+1 = гп + TnVn + тП 4(Еп+1 + Еп). (12)

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

Пусть точка гп принадлежит сеточному элементу е, уравнения сторон которого есть

/¿(ж,у) = 0' г = 1, 2,...'П3'

где п — число сторон. Решение уравнения (12) будем искать в пределах элемента е, для чего необходимо выбрать тп равным времени пролета частицы через е. Обозначим через гп+1 точку выхода частицы из рассматриваемого элемента. Пусть гп+1 € з^, где ^ — одна из сторон е, т. е.

/г(гп+1) = 0. (13)

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

Так как гп+1 лежит строго на одной из сторон элемента, напряженность Еп+1 в ней с точностью 0(к2) можно определить при помощи линейной интерполяции по значениям напряженности на концах стороны по формуле

-п+1

Еп+1 = Е1 + (Е2 _ Е1). ¿21

Здесь Е1' Е2 — напряженность в точках 1, 2, являющихся концами -21 — длина Si, а -п+1 — длина дуги между точкой гп+1 и узлом 1 (рис. 1).

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

?п+1 _ А

Еп+1 = р(Ап+1) = Е1 + А (Е2 _ Е1)' (14)

где Р = (Рс Ру). Здесь

= А _ А А Г х при ^ > |-2У1|' I У пРи |-21| < |-21|.

1

Рис. 1.

Подставляя (14) в (12) и исключая из получившихся соотношений координаты (п + 1)-й точки с помощью (13), получим уравнение четвертого порядка для определения тп — времени пролета частицы через е. Решение этого сложного уравнения требует значительных вычислительных затрат. В связи с этим найдем более простые соотношения, которые дают значение тп с приемлемой точностью.

Из разложения гп+1 в ряд Тейлора следует линейное уравнение для вычисления тп

тпvn + rn - Г+1 =0. (15)

Тогда, если уравнение Si есть

Лг x + Biy + Ci = 0,

то искомое время вычисляется как

= - Лг ХП + ВгУП + Сг (ifi)

тп,г л п , тэ п . (16)

Лг vn + Вг vn

Так как сторона s,i, которую пересекает траектория, заранее не известна, то окончательное время пролета частицы через элемент рассчитывается как

Тп = min тп,г, тп,г > 0, i = 1, 2, ...,ns. (17)

г

Определив таким образом тп, координаты гп+1 будем вычислять по формуле (12). Покажем, что и в этом случае Еп+1 можно рассчитывать при помощи линейной интерполяции. Из (12), (15) следует, что точки гп+1 и гп+1 близки друг к другу, а именно

Гп+1 - гп+1 = 0(т2).

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

Еп+1 = Е(гп+1) = Е(гп+1) + 0(т2) = Р(Г+1) + 0(т2) = Р(£п+1) + 0(т2). Запишем функцию Р в виде

Р(£ ) = а + Ъ£, где векторы а = (аж, ау), Ъ = (Ьх, ) выражаются как

£2 Е1 — £1Е^ Е2 — Е1 а =-т-, Ъ =-т-.

621 /21

Тогда координаты (n + 1)-й точки рассчитываются по формулам

_ Х + TnVx + ^n(Ex + ax)

1 - (18) уп+1 = уп + + + + Ьу хп+1)

при проведении интерполяции по х (£ = х) и по формулам

п+1 = Уп + Тп^П + жп(ЕП + ау)

У 1 - ЖпЬу ' (19)

хп+1 = хп + Тп< + Жп(Еп + аж + ЬхУп+1)

при проведении интерполяции по У (£ = у). Здесь жп = 0.25^^.

Итак, расчет (п + 1)-й точки траектории заряженной частицы в элементе е с прямолинейными сторонами состоит из следующих шагов:

— определение временного шага по формулам (16), (17);

— вычисление координат по формулам (18) или (19);

— расчет компонент скорости по формулам (6).

При вычислении времени пролета определяется сторона зе элемента, которую пересекает траектория. Следующим расчетным элементом будет элемент, имеющий зе своей границей. В этом случае при наличии соответствующей структуры данных изложенный подход не требует выполнения трудоемкой процедуры поиска элемента, которому принадлежит рассматриваемая точка. Исключение составляет начальная точка траектории. Если окажется, что гп+1 с точностью 0(т2) лежит вблизи одного из узлов г5, являющегося концом зе, то анализ положения точки необходим, но он проводится только среди элементов, имеющих г5 своей вершиной.

Отметим, что при поэлементной технологии т = О (Л), что дает возможность в приводимых оценках заменить т на Л.

3. Расчет напряженности электрического поля

Компоненты напряженности электрического поля

E , ,

dx dy

вычисляются в узлах сетки П до расчета траекторий по следующему алгоритму.

На рис. 2 изображен фрагмент сеточной области. Цифрами обозначены узлы сетки, а буквами s, 1 — линии сетки, пересекающиеся в узле 0, в котором требуется вычислить напряженность.

Тогда производные по направлениям s и I выражаются как

— _ — cos ax + — cos a,, ds ax dy y

"тгт _ t— cos ax + — cos a,, dl dx x dy y

s s i i где cos ax, cos ay — направляющие косинусы вектора s, а cos ax, cos ay — направляющие

косинусы вектора l в точке 0. Отсюда найдем производные по x и y:

1

рХ = 6 (pScos аУ- picos аУ), 1

рУ = 6 (Picos аХ- pS cos аХ )5

6 = cos aX cos aly — cos ay cos aX.

Значения производных p's, p\ рассчитываются по трехточечным схемам. Например,

dp 1

pS вычисляется как

ds soi + S03

— (po — Pi) + — (рз — Po) Soi S03

где вц — расстояние между узлами г,].

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

Приведенный алгоритм обеспечивает вычисление напряженности электрического поля с точностью 0(Н2).

4. Распределение объемного заряда

В общем случае при проведении расчетов с произвольным шагом по времени отрезок к-й траектории (к = 1, 2,..., , где — суммарное число траекторий) ёЩ = [гЩ, гЩ+1 ] может пересекать ^ > 1 сеточных элементов в{,г = 1, 2,..., N. Известно несколько способов распределения заряда с данного отрезка в сеточные узлы. Например, один из них заключается в том, что ёЩ разбивается на подынтервалы в предположении равноускоренного движения и заряд с каждого такого подынтервала относится к узлу, ближайшему к его центру [3]. Другой способ состоит в том, что ёЩ разбивается на такие подынтервалы, каждый из которых полностью лежит внутри ячейки Дирихле — Вороного, а заряд с него относится к центральному узлу ячейки [5]. Так или иначе, данные алгоритмы требуют проведения трудоемких действий по отысканию сеточных узлов, ближайших к заданной точке.

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

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

Требуется рассчитать распределение заряда д по узлам сетки П, которые являются вершинами е. В основу решения данной задачи положен алгоритм распределения единичного точечного заряда, изложенный в монографии Р. Хокни, Дж. Иствуда [6]. Его суть заключается в следующем.

Обозначим через Жр долю единичного заряда, находящегося в точке г, которая относится к точке гр — вершине е. Точный потенциал в точке г' £ е есть

¥(г') = С(г' - г), (20)

а приближенный потенциал равен

¥>(0 = £ ЖС(г' - гр), (21)

р

где С(г) — функция Грина, а суммирование проводится по всем узлам, которые являются вершинами е. Разлагая С(г' — гр) = С((г' — г) + Др) в ряд Тейлора по параметру Др = г — гр = (Д^, Др), получим

¥(г') = с 5>р+^Ежд + ¿¡тЕ жд+^рдхду+°(^2)- (22)

р р р р

Потребуем, чтобы точный потенциал (20) совпадал с приближенным потенциалом (21) с точностью О (Л,2), в результате чего получим уравнения для определения Жр. Ниже рассматриваются алгоритмы распределения объемного заряда для треугольного и четырехугольного сеточных элементов.

4.1. Треугольный элемент

Пусть е = е4 — треугольный элемент с вершинами, которые обозначим цифрами 1, 2, 3. В этом случае, взяв в (22) три первых члена (порядок точности при этом сохраняется), получим следующую систему уравнений:

3 3 3

рр

р=1 р=1 р=1

¿Жр = 1, ¿Хд = 0, = 0. (23)

Ее решение можно выразить как

Ж = АД^ — В1ДЗ, Ж2 = — В2Д1, Жз = АзД^ — ВзДУ,

где

/V /X

IX IV

А = /У /23 ¿1 ' В = /X /23 ¿1 ' ¿1 /X /V = /13/23

А = /У /31 ¿2 ' В = ¡X /31 ¿2 ' ¿2 _ /X 7У = /21 /31

со /У /21 ¿3 ' В = /X /21 ¿3 ' ¿3 /X 7 V = /31 /21

¡х /У

X

Здесь /Х = xi — х], /у = yi — У], где у — координаты ¿-го узла сетки.

Перейдем теперь к решению задачи о распределении заряда д, вносимого отрезком 8 £ et между и-й и (и + 1)-й точками траектории. Разобьем 8 на щ подынтервалов равной длины

/

dX

nn

i=iiI=jx+Ti.

Здесь 1 = (/Х, /У), где /Х = хп+1 — хп, /У = уп+1 — уп. Координаты центров этих подынтервалов Г определяются по формуле

Г, = rn , i =1,2,

2n

,n„

q

а заряд, вносимый каждым подынтервалом, равен д = д/щ. Сосредоточим заряд д в центре ¿-го подынтервала. Его вклад др в узлы сетки, являющиеся вершинами е4, выражается как

qi = Wz

Чр ' ' р

q

n„

где Шр — доля единичного заряда, находящегося в точке г^ отнесенного к р-му узлу.

Для Шр справедливы уравнения, аналогичные уравнениям (23). Их решение можно записать в виде

ш; = АрД^ — Врдy/)i, р,р' =1, 2, 3, г = 1, 2,..., щ, р = р'.

Заряд др, отнесенный к р-му узлу, со всего участка 8 равен сумме зарядов с каждого подынтервала, т. е.

др = Е & = уЕШ^-

i i

Окончательно заряд в р-м узле определим как

где

Так как по определению

АХ

^p' ,i

qp = Jim & = W pq>

l

Wp = If W;dx. 0

Xi - Xp', Ap'i = y, - yp'

(24)

а из геометрических соображений

X, = xn + X cos a, y, = yn + X sin a,

где a — угол наклона рассматриваемого отрезка траектории к оси x, а X — длина отрезка траектории, ограниченного точками rn и г,, то (24) можно записать в виде

i i Ai B Г

Wp = -y (xn - xp' + X cos a)dX--^ (yn - yp' + X sin a)dX.

00

Проводя интегрирование и преобразования для Жр, получим окончательное выражение

Ж = Лр(жга+1/2 - жр) - Вр(уп+1/2 - ур),

где жга+1/2, уп+1/2 — координаты центра отрезка траектории, определяемые по формуле (10).

Таким образом, доказано следующее

Утверждение 1. Вклад объемного заряда q, вносимого отрезком траектории 5, принадлежащим треугольному элементу е^ в его вершины равен вкладу точечного заряда q, сосредоточенного в центре (жга+1/2,уп+1/2) отрезка 5.

Из построения qp следует, что при данном распределении заряда погрешность £ расчета потенциала удовлетворяет равенству

£ = - ^ = 0(Л,2).

4.2. Четырехугольный элемент

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

= 1, АрХ = 0, ХУрДР = 0, Ар

р=1

р=1

р=1

р=1

Ее решение можно выразить как

Ж = С^АУ - ^1А|АУ,

Ж = С2А^А| - ^А!АУ,

СзАХАУ - ДзА|АУ,

(25)

где

С = С2 =

]х ¡У

/23/24

¿1 /х /У

113114

¿2

/х

24

Сз = ТГ" С2 ,

/Х

1з

Ж4 = С4АХАУ - Д^АУ,

£1 = £ =

/х /У

/24/23

¿1 /х /У

114113

„ _ 7Х 7Х 7У 7У 7Х 7Х 7У 7У

¿1 = /14/2.3/1.3/24 /13/24/14/23,

¿2 У

„ _ 7Х 7Х 7У 7У 7Х 7Х 7У 7У

^2 = /24/13/23/14 /23/14/24/1:

£ = /|4 £>2,

113

/У

11 о

23 14 24 13

/!

Ь1

С4 = тУ3 С1, £4 = 711 А-

/У

124

/Х

124

(26)

Отметим, что в случае прямоугольного сеточного элемента формулы (25), (26) существенно упрощаются и переходят в формулы, приведенные в монографии Р. Хокни, Дж. Ист-вуда [6].

Выражения для Жр можно записать в виде

Ж,

СрАХАт - ДрА^Ак, р,к,т = 1,2, 3,4, р = к = т.

^р^т^к'

4

4

4

4

0

Чтобы получить выражения для Жр, определяющие вклад заряда д со всего отрезка траектории, так же, как и в случае треугольного элемента, разобьем данный отрезок на равные подынтервалы и перейдем к пределу при ^Л ^ 0. При этом в Жр войдут интегралы вида

i

0

Проводя интегрирование и преобразования, получим

I* =(хП+1/2 - Хк)(уП+1/2 - Ут) + ^.

Отсюда следует

Утверждение 2. Вклад объемного заряда д, вносимого отрезком траектории 5, принадлежащим четырехугольному элементу е*, в его вершины определяется по формулам

qP = Wpq, Wp = + ^(Cp - Dp), p =1, 2, 3, 4,

где Жр,* — вклад от единичного заряда, находящегося в центре 5.

При этом так же, как в случае треугольного элемента, обеспечивается второй порядок точности расчета потенциала.

5. Численный эксперимент

Рассматривалось движение электрона в поле, потенциал которого определяется по формуле

cosh 2y — cos 2x , .

= ---, (27)

cosh 2y + cos 2x

где — декартовы координаты (см., например, [7]). Координаты траектории удовлетворяют соотношению

cosh 2y + cos 2x = Ca = const, (28)

а скорости вычисляются по формулам

sinh 2y sin 2x

cosh 2y + cos 2x cosh 2y + cos 2x

Рассчитывалась траектория, выходящая из точки x0 = п/8, y0 = 0. Расчет проводился в области, изображенной на рис. 3 на квазиструктурированной сетке П = П У П2 U Пэ, состоящей из равномерных прямоугольных сеток Hi = N х N, Пэ = N х 2N и равномерной треугольной сетки П2. Здесь N — число разбиений по одной переменной. Треугольная сетка имела N разбиений по каждой стороне. В точках траектории вычислялась относительная ошибка по формуле

¿n% =

сn — с ca

100%, n = 1, 2,

к к я х

¥ 7 ~2

Рис. 3.

где величина Сп рассчитывалась по формуле (28), в которую подставлялись значения координат, полученные в результате численного интегрирования уравнений движения по изложенным выше алгоритмам.

В таблице приведены значения максимальной ошибки в процентах 8% = тах 8п% на сетках П с различным числом узлов N. При этом 81% означает максимальные ошибки в расчетах траектории, в которых напряженность поля вычислялась по аналитическим формулам, полученным дифференцированием (27), а 82% — максимальные ошибки в расчетах траектории с численным расчетом напряженности по рассмотренным выше алгоритмам.

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

Ошибки в расчете траектории

N ¿1% ¿2%

8 1.1 1.2

16 2.4 ■ 10-1 3.5 ■ 10- 1

32 6.4 ■ 10-2 1.1 ■ 10- 1

64 1.6 ■ 10-2 3.0 ■ 10- 2

128 4.0 ■ 10-3 8.0 ■ 10- 3

256 1.0 ■ 10-3 2.1 ■ 10- 3

512 2.6 ■ 10-4 5.5 ■ 10- 4

1024 6.4 ■ 10-5 1.4 ■ 10- 4

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

[1] Ильин В.П., Павлов М.В., Свешников В.М. Решение двумерных краевых задач на квазиструктурированных сетках // Тр. Междунар. конф. RDAMM-2001. Новосибирск, 2001. Т. 6, Ч. 2, спецвыпуск. С. 51-59.

[2] Молоковский С.И., Сушков А.Д. Интенсивные электронные и ионные пучки. М.: Энергоатомиздат, 1991.

[3] Ильин В.П. Численные методы решения задач электрофизики. М.: Наука, 1985.

[4] Ильин В.П., Свешников В.М., Сынах В.С. О сеточных технологиях для двумерных краевых задач // Сиб. журн. индустр. мат. 2000. Т. 3, № 1. С. 124-136.

[5] Свешников В.М. Численный расчет пучков заряженных частиц на локально-модифицированных сетках. Новосибирск, 1997 (Препр. РАН Сиб. отд-ние. № 1109).

[6] Хокни Р., Иствуд Дж. Численное моделирование методом частиц: Пер. с англ. М.: Мир, 1987.

[7] Сыровой В.А. Проблема адекватности математических моделей в оптике плотных релятивистских электронных пучков // Радиотехника и электроника. 2003. Т. 48, № 4. С. 467.

Поступила в редакцию 9 января 2003 г., в переработанном виде — 2 марта 2004 г.