CC BY
71
Вычислительные технологии
Том 11, № 5, 2006
РЕШЕНИЕ САМОСОГЛАСОВАННЫХ ЗАДАЧ ЭЛЕКТРОННОЙ ОПТИКИ МЕТОДОМ ИТЕРАЦИЙ ПО ПОДОБЛАСТЯМ БЕЗ НАЛЕГАНИЯ И СМЕНЫ ТИПА ГРАНИЧНОГО УСЛОВИЯ*
В. М. Свешников Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: [email protected]
Numerical algorithms for the solving of the self-consistent problems of the electronic optics, based on decomposition of calculated domain on the near cathode and the basic subdomains are constructed and realized. Specifics of the present approach is that during the iterative process the matching between subdomains is carried out without intersections, and at the interface boundary the Dirichlet condition for the potential of the electric field is set in both subdomains. Solution of the self-consistent problem is found as the solution of the system of nonlinear equations. Each equation of the system has the difference between normal components of the electric field to the left and to the right of. Iterative algorithms for solving the equations are proposed and their convergence is investigated numerically on a test problem.
Введение. Постановка задачи
Самосогласованные задачи электронной оптики, которые рассматриваются в настоящей статье, заключаются в расчете движения пучка заряженных частиц во внешних и собственных электрических полях. Пучок предполагается интенсивным, при расчете которого нельзя пренебречь объемным зарядом, создаваемым частицами. Такие задачи имеют большое практическое значение, так как возникают при разработке многих электрофизических устройств, например приборов СВЧ-электроники, сильноточных ускорителей, пушек для электронно-лучевой сварки [1]. Согласно монографии [2], математическая постановка нелинейной самосогласованной задачи электронной оптики, которая рассматривается в настоящей статье, состоит в следующем. В замкнутой области G = G U Г, где Г — граница G, требуется найти решение нерелятивистских уравнений движения частиц с зарядом ei и массой mi в электрическом поле с потенциалом ip и напряженностью E:
dri dvi ei - -
и = и = mE E = -gradp- (1)
т. е. определить координаты — и скорости -i частиц, где i = 1, 2,... — номер частицы,
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-01-00487).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.
і — время. Уравнения (1) решаются с начальными условиями
^4=0 = ^ ^4=0 =
(2)
где г0 иг/0 — заданные величины, представляющие собой соответственно координаты и скорости частиц, входящих в расчетную область. Потенциал электрического поля р находится из решения уравнения Пуассона
Ар
с краевыми условиями
Р\
Гі
д,
др
дп
Є0
0,
(3)
(4)
Г 2
где р — плотность объемного заряда, вносимого частицами; е0 — электрическая проницаемость вакуума; д — известная функция координат; Г1, Г2 — куски границы Г такие, что Г и Г2 = Г; п — внутренняя нормаль к Г2. Для плотности тока . потока заряженных частиц выполняются уравнения неразрывности
divJ = 0, <7
ру,
ІІШ
(5)
Ду^ а ’
а на катоде С С Г, представляющем собой поверхность старта заряженных частиц, задаются граничные условия для плотности тока
1
с
і•
(6)
В формуле (5) у — средняя скорость, причем суммы берутся по всем частицам, находящимся в объеме А У.
Наиболее важным для практики вариантом работы катода является режим ограничения тока объемным зарядом (р-режим работы катода). Тогда задача (1)-(5) решается со следующими граничными условиями:
Р\
с
0,
Е = дР
1с дп
0,
(7)
с
где Е — модуль напряженности электрического поля. При этом граничное условие (6) не задается, т. е. плотность тока на катоде является неизвестной величиной. Поясним, что первое условие (7) используется как граничное условие Дирихле при решении уравнения Пуассона, а второе — для нахождения плотности тока на катоде. В дальнейшем модуль вектора плотности тока на катоде (или просто плотность тока на катоде) будем обозначать
через і, т. е. ] = 1 . Отметим, что при решении данной задачи в (2) полагается у0 = 0. В данном варианте самосогласованной задачи наряду со многими проблемами одной из важнейших является проблема нахождения плотности тока і на катоде. Решению именно этого вопроса посвящена данная статья.
Приведенная постановка определяет стационарную задачу. Время і в ней играет роль параметра и, вообще говоря, может быть исключено [1]. Тогда уравнения движения заряженных частиц переходят в уравнения траекторий, которые имеют более сложный вид по сравнению с формулами (1).
Рассматриваются двумерные (плоские или осесимметричные) задачи в декартовой х, у или цилиндрической г, г системах координат. Моделирование пучка заряженных частиц
осуществляется методом трубок тока [2], который представляет собой экономичный вариант широко известного метода больших частиц [3]. Интегрирование уравнений движения и решение уравнения Пуассона проводятся численно [2].
Наиболее часто для отыскания плотности тока і на катоде применяется итерационный метод нижней релаксации, на каждом приближении которого используется закон “3/2” для плоского диода [2]. Такой подход подкупает своей простотой, но может привести к большим ошибкам в расчетах [4, 5], поскольку при этом не учитываются такие существенные факторы, как криволинейность катода и неоднородность плотности тока на нем. Кроме того, самосогласованные задачи при эмиссии в р-режиме имеют особенность, заключающуюся в том, что плотность объемного заряда на катоде обращается в бесконечность. Большинство алгоритмов и программ для их решения строится без выделения прикатод-ной особенности, что также может породить значительные ошибки в расчетах.
Применительно к данным задачам метод итераций по подобластям выполняет две функции: 1) является методом нахождения плотности тока на катоде; 2) служит средством повышения точности расчетов. Он основан на декомпозиции расчетной области О на прикатодную Ос и основную Оъ подобласти. Повышение точности достигается за счет того, что стартовая поверхность входа частиц в область при численных расчетах отодвигается от особой поверхности катода, а решение в прикатодной области находится по аналитическим формулам. Антипараксиальная теория В.А. Сырового [6] наиболее полно и адекватно отражает физическую картину явлений в прикатодной подобласти, и поэтому именно она является основой для построения решения в ней.
В работах автора [7, 8] рассматривался вариант метода итераций по подобластям без налегания для решения данных задач, в котором антипараксиальные разложения в прика-тодной подобласти строились по распределению производной от потенциала, а численные расчеты в основной подобласти проводились с использованием распределения самого потенциала на границе сопряжения подобластей Гсъ, т. е. при переходе из одной подобласти в другую менялся тип граничного условия. В настоящей работе применяется другой подход, который заключается в том, что при проведении итерационного процесса по подобластям без налегания на Гсъ ставится условие Дирихле для потенциала электрического поля в обеих подобластях. При этом задача отыскания плотности тока на катоде формулируется как задача решения системы нелинейных уравнений, в каждое из которых входит разность между нормальными составляющими электрического поля слева и справа от Гсъ. Идея данного подхода опубликована в [9].
В настоящей работе преимущественное внимание уделено исследованию различных видов итерационных процессов по подобластям для решения рассматриваемых задач с целью построения быстро сходящегося итерационного процесса. Основным инструментом проведения исследований при этом является численный эксперимент, так как “для нелинейных задач первостепенное значение для проверки качества численных методов имеют тесты, т. е. численное решение частных задач — типичных представителей класса решаемых задач, для которых известны аналитические решения” [10, с. 457].
1. Методы отыскания плотности тока
Решение нелинейных самосогласованных задач (1)-(5), (7) с неизвестной плотностью тока і на катоде С осуществляется каким-либо итерационным методом, на н-м приближении которого по тем или иным правилам вычисляется очередное распределение іп на С. Затем
решается самосогласованная задача (1)-(6) расчета пучка заряженных частиц, который эмитируется с катода С с заданной плотностью тока, равной іп. Решение данной задачи дает новые значения величин, по которым пересчитывается плотность тока. Описанный процесс повторяется до тех пор, пока не будет достигнута его сходимость по плотности тока с заданной точностью.
Методы отыскания плотности тока і отличаются как видом итерационного процесса, так и способом вычисления і на каждой итерации. Их можно разбить на две группы.
К первой из них (см., например, [2]) относятся методы, использующие закон “3/2”, согласно которому в плоском диоде (0 < х < й) с нулевым потенциалом на катоде и потенциалом на аноде ра, где й, ра заданы, рассматриваемая задача имеет аналитическое решение, в котором плотность тока і и текущий потенциал р определяются по формулам [1]
где в] — константа, зависящая от сорта частиц, например, для электронов о^ ~ 2.33 • 10-6.
Для отыскания плотности тока в двумерной задаче строится итерационный процесс нижней релаксации вида
где н = 0,1,... — номер итерации; 0 < 0 < 1 — заданное число; і0 = 0. На каждой итерации (9) численно решается двумерная самосогласованная задача (1) — (6) с известной плотностью тока іп. Ее решение дает потенциал ра = рП на заданном малом расстоянии й
после чего вновь решается самосогласованная задача с заданной плотностью тока, равной уП+1. Данному подходу можно дать следующую интерпретацию. Прикатодная подобласть представляется в виде совокупности N1 плоских диодов, в каждом из которых потенциал анода постоянен (N4 — число трубок тока).
Такое использование “одномерного” закона “3/2” при решении двумерных задач приводит к ошибкам из-за криволинейности катода и неоднородности плотности тока вдоль него.
Кроме перечисленных ошибок следует отметить, что применение закона “3/2” для плоского диода при решении двумерных задач в итерационном процессе (9) является в принципе непоследовательным. Действительно, основой его использования служит утверждение
о том, что на малом расстоянии от катода распределение потенциала вдоль нормали к поверхности катода ведет себя приблизительно так же, как в плоском диоде. Это действительно так, но данное положение используется только при нахождении плотности тока, а потенциал, напряженность, объемный заряд, фазовые координаты пучка во всей области, включая прикатодную подобласть, вычисляются из решения двумерной задачи. При этом объемный заряд вблизи катода огромен, а потенциал, наоборот, мал по сравнению с их значениями вне прикатодной подобласти, что приводит к большим ошибкам в расчетах вблизи катода и, как следствие, к существенным ошибкам в расчете пучка заряженных частиц в целом. На наличие ошибок указывалось в работе [4]. Их устранению не помогает введение локально-модифицированных сеток и специальных прикатодных аппроксимаций потенциала, как это показано в работе [5].
Ко второй группе методов отыскания плотности тока у можно отнести методы, обеспечивающие непосредственно выполнение второго условия (7), т. е. условия равенства нулю напряженности электрического поля Е на катоде. Здесь в первую очередь необходимо
(8)
і
•п+1
0іп+1/2 + (і _ 0^.п,
(9)
от катода, по которому согласно (8) вычисляется іп+1/2, а затем по формуле (9) — іп+1,
Таблица 1. Погрешности численного расчета плоского диода
16 32 64 128
^ 1.2 1.8 2.5 3.2
отметить работы Г.Т. Головина (см., например, [11]), Н.И. Мокина [12]. Данному вопросу посвящена также работа [13] автора настоящей статьи. В ней предлагается подойти к проблеме нахождения плотности тока как к решению системы нелинейных уравнений
где Ф — заданная функция от напряженности электрического поля на катоде. Такой подход открывает возможность широкого применения известных итерационных методов для решения нелинейных уравнений, что позволяет избежать трудностей, с которыми сталкивался Г.Т. Головин в своих работах.
И в той и в другой группе методов нахождения плотности тока не учитывается особенность решения, природа которой заключается в том, что на катоде плотность объемного заряда обращается в бесконечность.
Для того чтобы убедиться, насколько велико влияние особенности в решении данных задач, автором настоящей статьи были проведены эксперименты по численному решению одномерной задачи (1)-(5), (7) о плоском диоде с параметрами 0 < х < 1, =1, аналити-
ческое решение которой дается формулами (8). Задача решалась на сетках с различным числом узлов N^ = 2Й+3, к = 1, 2, 3,4. Решение самосогласованной задачи проводилось при помощи итерационного процесса (9), в котором на каждом приближении плотность тока вычислялась по формулам (8) на расстоянии = 1/8. В табл. 1 даны значения относительной погрешности ^ расчета потенциала в процентах, которая определялась по формулам
где р* — точные значения потенциала, а р; — его приближенные значения в узлах сетки.
Приведенные результаты свидетельствуют о том, что погрешности не только не падают с увеличением числа узлов, а наоборот, даже возрастают.
2. Формулировка задачи отыскания плотности тока на основе декомпозиции расчетной области
Разобьем расчетную область О линией Гсъ на две подобласти Ос и Оь такие, что О = Ос и
Оь и Гсъ. Их замыкания есть Ос = Ос и Гс, Оь = Оь и Гь, где границы Гс, Гь удовлетворяют равенствам Гс = Г° и Гсь, Гь = Г0 и Гсь, Г° и Г0 = Г, причем Г° включает катод С. В дальнейшем Ос будем называть прикатодной подобластью, а Оь — основной подобластью.
Решение исходной самосогласованной задачи в О с неизвестной плотностью тока і на С будем искать путем построения последовательных приближений в данных подобластях. На каждом приближении проводится решение следующих самосогласованных задач в Ос и Оь. Предположим, что на Гсь задано условие Дирихле для потенциала р вида
Ф(£ (і)) = 0,
(10)
(11)
р|г„ =
(12)
где ги = ги(Т) — известная функция точки Т Є Гсь. Подобласть Ос построим достаточно простой, чтобы стало возможным на основе условия (12) по известным аналитическим формулам рассчитать плотность тока і и пучок заряженных частиц. Данные расчеты осуществляются при помощи антипараксиальной теории В.А. Сырового [6], которая наиболее полно и адекватно описывает явления в прикатодной подобласти, в результате чего определяем координаты и скорости заряженных частиц на Гсь (см. разд. 4 настоящей статьи). Используя величины, полученные в результате расчетов в прикатодной подобласти, и условие (12), решаем самосогласованную задачу в основной подобласти Оь с заданной плотностью тока на поверхности старта заряженных частиц, которой в этом случае является Гсь. Решив данную задачу, по тем или иным правилам, о которых будет идти речь ниже, пересчитываем функцию и(Т) и описанный процесс повторяем на следующем приближении.
Проблема заключается в том, чтобы в результате проведения последовательных приближений найти на границе Гсь такое распределение потенциала, которое давало бы решение исходной задачи во всей области О, включая нахождение плотности тока, обеспечивающей равенство нулю электрического поля на катоде.
“Сшивку” подобластей будем осуществлять без их налегания. Запишем условия сопряжения на Гсь в виде
(р)с = (р)ь; (13)
. І). - (Й). ■ "4)
где п — нормаль к Гсь, а индексы с, Ь означают, что величины в скобках относятся соответственно к подобластям Ос и Оь. Направление нормали в данном случае не важно, но для определенности положим, что П — это внутренняя нормаль по отношению к Оь. Условия (13), (14) должны выполняться на решении исходной задачи в О. Если же на Гсь задан произвольный потенциал и, то условие (13) выполняется, а условие (14) — нет. Производные в формуле (14) являются функциями, зависящими от и. С другой стороны, потенциал ги согласно антипараксиальной теории определяет плотность тока і = і (и). Введем обозначения
* ( ) (дрН А * ( ) (дрН А (15)
/с («>)=( с , Л («>)=( -дпр) ь . (15)
Сформулируем исходную задачу по отысканию плотности тока і следующим образом: найти такое распределение потенциала и, при котором выполняется условие (14). Это требование эквивалентно решению операторного уравнения
Тги - /с (ги) - / (и) = 0, (16)
где Т : Ш ^ Ш — нелинейный оператор, действующий в пространстве Ш функций,
определенных на Гсь.
Уравнение (16) будем решать приближенно. Для этого построим на Гсь сетку ш = {Ті Є Гсь, і =1, 2,... , N} с шагами ^ = І (Ті+1) — І (Ті) > 0, і = 0,1,... , N, где N — заданное целое число, а І (Ті) — длина отрезка [Т0,Т] границы Гсь, отсчитываемая от точки Т0. Пусть и — пространство сеточных функций (векторов), определенных на ш, и = {« — и (Ті) , і = 1, 2,... , N} — вектор пространства и, а Р : Рго (Ті) = ги (Ті) = « — оператор проектирования из Ш в и. Действуя оператором Р на равенства (15), получим
Р/с (и.,7у= /с, - /с(и, Ті) = (др(и)А , Р/ь (и, Ті) = /м - /ь(и, Ті) = (др (и)
Заменим задачу (16) отыскания функции ,ш непрерывного аргумента приближенной задачей нахождения функции и дискретного аргумента, состоящей в решении операторного уравнения
Фи = /с - / = 0, (17)
где Ф : и ^ и — нелинейный оператор; /с = </}, / = {/м} (% = 1, 2,... , N) —
векторы пространства и. Фактически тем самым мы требуем выполнения условий (13), (14) не во всех точках Т € Гсь, а лишь в точках Т € ^. Операторное уравнение (17) можно представить как систему нелинейных уравнений путем введения вектор-функции
Ф(и) = Ши)}
с компонентами
Ф* = Фг(и) = (Фи)* = /с, * - /*, % = 1, 2, . . . , N.
Тогда уравнение (17) можно представить в виде системы нелинейных уравнений относительно вектора и
Ф(и) = 0, (18)
для решения которой можно применять известные итерационные методы (см., например,
[14]). В следующем разделе настоящей статьи предлагается ряд итерационных алгоритмов для решения системы (18), сходимость которых исследуется путем проведения численных экспериментов на тестовой задаче.
3. Итерационные методы решения системы нелинейных уравнений относительно потенциала на границе сопряжения подобластей
Систему нелинейных уравнений (18) будем решать методом простой итерации вида
«п+1/2 = «п + т„Фга«0; (19)
ип+1 = I « , ||ф 11 < ||ф ||, (20)
« 1 „,П ИгЬга+Шп \ I |мм (20)
ип+1/2, ||Фп+1/2|| < ||ФП||,
ип, ||Фп+1/2|| > ||ФП||.
Здесь п = 0, 1,... — номер итерации; Фп = diag {ФП/Ф0} — диагональная матрица; Фп = |Ф™} € и — вектор, содержащий значения Фг на п-й итерации (% = 1, 2,...,N); ||.|| — какая-либо норма в пространстве и, а тп — числовой параметр, который вычисляется как
( Т„, ||Фп+1/2|| < ||ФП||,
ТП+1 = |± Г„, ||ФП+1/2|| > ||ФП||, (21)
где то — задано. В работе [9, с. 51] приведена технология расчета и0, которая обеспечивает Ф0 = 0 для интенсивных пучков, рассматриваемых в настоящей статье. Итерационный процесс прекращается при выполнении хотя бы одного из двух условий
||ФП|| < £, т„ < £, (22)
где є — заданная малая величина.
Наряду с итерационным процессом (19)—(21) будем рассматривать его модификацию, в которой, учитывая, что и = 0 из физических соображений не является решением системы (18), в (19) и0 заменим на ип, в результате чего получим
ига+1/2 = (I + т„Ф га)ига, (23)
где I — единичная матрица.
Одним из наиболее быстрых методов решения системы (18) является метод Ньютона:
ип+1 = ип - [ф'п]-1 фп. (24)
Здесь
( дФп
ф'"={*=а="■2-"ЛГ
— матрица Якоби вектор-функции Фп. Вычисление элементов матрицы Ф/га будем проводить по приближенной формуле
где — заданная малая скалярная величина; ек — вектор, к-я компонента которого равна 1, а остальные — нулю. Вычисление всех N2 элементов дП^, к = 1, 2,... , N, с помощью приведенной формулы требует решения N вспомогательных самосогласованных задач со значениями потенциала на ^, равными ип + ек, к = 1, 2,..., N, которое необходимо
проводить на каждом шаге итерационного процесса (24).
Существенного сокращения числа вспомогательных задач можно достичь, если считать, что матрица Ф/га приближенно является ленточной. Обозначим через 2М — 1 ширину ленты, т. е.
« 0 |і — к|> М. (26)
Покажем, что в данном случае для отыскания дП^ при | і — к| < М достаточно решения 2М — 1 вспомогательной задачи. Действительно, при малом |8и| имеем
Фп (мга + е) « Фп (мга) + Ф'га (мга) е, (27)
где е = (е* = 0,1, і = 1, 2,... , N} — вектор.
Введем векторы ем,к (к = 1, 2,... , 2М — 1), которые определим следующим образом:
ем,к Г ем,к \ ем,к Г 1 і = (2М — 1) v +k, (28)
е = 1е /' е‘ І0 і =(2М — 1)V + к. (28)
Здесь V = 0, 1,... , N1 1, 2,... , N. Решим (2М —
2М----1 (в квадратных скобках — целая часть числа); і =
)-ю вспомогательную задачу со значениями потенциала + ем,к (к = 1, 2,... , 2М — 1) на ^. Тогда из (26)-(28) следует, что, заменив в (25) вектор ек
м к ... 7 ...
на вектор ем,к, в результате решения каждой к-й задачи мы определим элементы матрицы Ф/га, лежащие на столбцах с номерами (2М — 1) V + к, V = 0, 1,... , ^, а решив (2М — 1) задачу, — все ненулевые элементы этой матрицы.
Отметим наиболее важный для практики случай. Полагая М =1, мы получим “одномерный” метод секущих, который запишем в виде
п+1 __ п <¥/пдчп
ип+1 = — ф'п Фп, (29)
?,п — ип-1 \/п _ А^гс ) %
где Ф/п = &ае < — - —і, і = 1, 2,... , N > — диагональная матрица. Применение данно-
Ф п Ф п- 1
го метода не требует решения вспомогательных задач.
4. Решение самосогласованных задач в подобластях
На каждом шаге итерационных процессов, рассмотренных в предыдущем разделе, мы имеем значения потенциала и = (и^, г = 1, 2,... , N} (верхний индекс п опускаем) в узлах сетки ^. Построим по этим значениям приближающую функцию ф(и,Т), которая дает значения потенциала и>(Т) в любой точке Т € Гсь, т. е.
цт ) = д(и,т).
Значения ,ш = эд(Т) будем использовать в качестве граничных условий на Гсь при решении самосогласованных задач в прикатодной Сс и основной Сь подобластях.
Покроем расчетную область С сеткой П^, узлы П которой представим в виде
П = Пс и Пь, (30)
где Пс — узлы прикатодной подобласти, включая узлы, лежащие на Гсь, а Пь — узлы основной подобласти. Потенциал электрического поля и объемный заряд вычисляются в узлах П. Вычисления на Пс проводятся по аналитическим формулам, а на Пь — численно. При расчете потенциала на Пь в узлах, лежащих вблизи Гсь, в качестве граничных условий берутся значения из соседних узлов, входящих в Пс. Разбиение (30) позволяет скрыть алгоритмическую “нефизическую” декомпозицию расчетной области.
4.1. Прикатодная подобласть
Мы будем предполагать, что прикатодная подобласть Сс имеет специальный вид, при котором граница сопряжения Гсь подобластей Сс и Сь представима в виде Гсь = Г^ и Г^, (рис. 1), где Г^ь является геометрическим местом точек, равноудаленных от катода на расстояние й по нормали к нему, а Г^ — это отрезки нормали, ограничивающие пучок заряженных частиц. Такой выбор прикатодной подобласти обеспечивает равное удаление точек старта заряженных частиц от особой поверхности катода при расчетах в Сь. Отметим также, что на практике катод с достаточной точностью может быть аппроксимирован отрезками прямых и дугами окружностей, что в данном приближении обеспечивает существование достаточно гладкой границы Г^.
Введем в подобласти Сс ортогональную систему координат 5, I с точкой начала координат О, расположенной на катоде, где 5 < й — расстояние по нормали к катоду, а
I — длина дуги вдоль катода. В работах В.А. Сырового (см., например, [6]) были получены разложения потенциала электрического поля в прикатодной подобласти в ряд по степеням малого параметра 5. Эти разложения получили название антипараксиальных, чтобы подчеркнуть тот факт, что они применимы не только для описания узких пучков
вблизи оси симметрии, но и для описания широких пучков большого радиуса. Антипарак-сиальные разложения представляют собой по сути дела нелинейные дифференциальные уравнения относительно плотности тока. Приближенное аналитическое решение данных уравнений [15] дает следующие выражения для плотности тока ] на катоде:
w
3/2
3 = С
і2(Т + в?Ь)3/2'
(31)
где
Т = Т(в) = 1 + — Те + — Т2в2 - — Ив2; у ; 15 225 18 ’
2 1 ь = __(и,^ - _ ___^
1 і.
IV =--------------IV
1 а2. .Ж'
(32)
(33)
(34)
Приведем пояснения к формулам (31)—(34). Величины Т, К, К1, входящие в формулы (32), (33), — это кривизны катода, которые выражаются как
Т = + К = К = П' = П" = 2>
Т = К1 + к2, К = К1 к2, К1 =------77) К1 =----777, к2 = ~77,
їх А їх
где г = Д(/), г = А(I) — параметрические уравнения катода. В плоских задачах здесь под г, г понимаются соответственно х,у и полагается к2 = К1 = 0. Отметим, что член Т “отвечает” за криволинейность катода, а Ь — за неоднородность плотности тока 3. Дифференцирование в (34) проводится по переменной /, которая представляет собой длину дуги вдоль кривой Г^, отсчитываемую от точки О на ней , которая является нормальной проекцией О на Г^.
Из последних равенств следует, что потенциал р во всей прикатодной подобласти Сс, в том числе на Г^, рассчитывается по формуле
(в\4/3 Т(в) + в2Ь
р = \і) т (і) + і2ь ,
а нормальная производная от потенциала, входящая в (15), вычисляется на Г^ как
др \ 4u (F (d) + 0.75dF'a (d) + 2.5d2L)
ди)с 3( (Т (() + (2Ь)
Координаты траектории г-й частицы, стартующей из точки (0, 1г,о), согласно антипарак-сиальной теории [6] описываются уравнением
, , _ £ Ф'
г г’° 10^ (г,
а ее скорость уг выражается через потенциал как уг = у/—2рег/тг. Для определения составляющих скорости рассчитываются углы вг наклона траекторий к нормали по формуле
(1г
вг = аге1^-—. Расчет объемного заряда осуществляется по алгоритму, изложенному в [16]. (^г _
Отметим, что размер прикатодной подобласти Сс мы предполагаем таким, что должно
выполняться неравенство
(> к, (35)
где к — максимальный шаг сетки в подобласти Сс.
4.2. Основная подобласть
Решение самосогласованной задачи в основной подобласти Сь проводится численно по алгоритмам, изложенным в [2, 5]. Поверхностью старта заряженных частиц является Г^. В качестве начальных данных для пучка при проведении расчетов в основной подобласти берутся величины, полученные из расчетов в прикатодной подобласти.
Вычисление нормальной производной на Г^ в Сь также осуществляется численно по формуле
др\ др др
— =— cos ах + — cos ay,
on/ b дх ду
где cos ax, cos ay — направляющие косинусы нормали, а x,y — декартовы или цилиндрические координаты (в последнем случае x,y означают соответственно r,z). Производные др/дх, др/ду, а также потенциал в какой-либо точке на Г^, которая в общем случае не совпадает с узлом сетки, вычисляются с помощью билинейной приближающей функции, построенной по значениям соответствующих величин в вершинах сеточного элемента, содержащего расчетную точку, согласно принципам поэлементной технологии [17]. Расчет производных от потенциала (34) также проводится численно по разностным формулам второго порядка точности.
5. Численные эксперименты
Численные эксперименты преследовали две цели: 1) показать повышение точности расчетов с помощью предложенного метода по сравнению с традиционными подходами; 2) исследовать вопросы сходимости итерационных алгоритмов. Итерационные процессы, рассмотренные в разд. 3, останавливались при достижении точности є = 10-4 в критериях (22). Такой выбор є неслучаен: проводились предварительные численные эксперименты с уменьшающимися значениями є, в результате чего установлено, что данное значение
является пороговым и его дальнейшее уменьшение не повышает точность расчетов, так как при этом главенствующую роль начинают играть неустранимые погрешности аппроксимации задачи.
Численные эксперименты показали, что во всех случаях итерационные процессы по подобластям прекращались при выполнении первого неравенства (22). При этом значение итерационного параметра тп, начальное значение которого т0 = 1, уменьшалось не более чем в восемь раз.
В качестве тестовой взята задача о потоке Мельтцера [18] с существенно неоднородной плотностью тока на катоде, имеющая аналитическое решение. Она состоит в следующем.
Электронный пучок движется в поле, аналитическое выражение для которого есть
2р = г-2 вт4/3 ^2Ф^ , г = Vх2 + У2, Ф = аг^ (36)
Катодом является отрезок прямой х = 0, траектории представляют собой концентрические окружности с центром в начале координат, а плотность тока равна
9с, — 5
3 = —У •
8л/2
Для численных расчетов выбиралась область, изображенная на рис. 2, где показаны поведение траекторий (точечные линии). Радиус внешней окружности Я = 4. Катод — это отрезок х = 0, 2 < у < 3.95. На границе области ставилось условие для потенциала согласно формуле (36).
Расчеты проводились на прямоугольных сетках П = П с различным числом узлов и при различных прикатодных расстояниях 4 = 4к. Число узлов в сетках П по х и у соответственно было равным тX = 7 х 2г, тУ = 2г+2, г = 1, 2, 3, что дает приблизительно равные шаги по обеим переменным. Прикатодные расстояния выбирались равными 4к = Я/(7 х 2к—!), к = 1, 2, 3.
Вычислялись следующие погрешности: 5^ — максимальные относительные погрешности расчета потенциала в процентах по формулам, аналогичным формулам (11), 5,, 54 — средние относительные погрешности расчета в процентах соответственно плотности тока
Рис. 2. Расчетная область и траектории заряженных частиц в задаче о потоке Мельтцера.
и траекторий — по формуле
і *
N . 1
г=1
Здесь N — число трубок тока (траекторий); 8 = 8, 8*, а 8г = 8\,8* — относительные погрешности в процентах расчета плотности тока и траектории для г-й трубки тока, причем 8^ определяются аналогично погрешностям расчета потенциала, а 8* вычисляются как отклонение численного значения координаты по х от соответствующего точного значения, вычисляемого по формуле (36) при заданном у.
С целью оценки влияния учета антипараксиальных разложений на точность вычислений прикатодная область рассчитывалась в двух вариантах:
— с учетом производных от плотности тока, т. е. согласно антипараксиальной теории;
— без учета производных от плотности тока, т. е. фактически по закону “3/2” для плоского диода.
В табл. 2, в которой даны погрешности расчета потока Мельтцера, для каждого набора счетных параметров приведены две цифры: верхняя — относится к первому варианту расчета прикатодной области, а нижняя — ко второму; пустые клетки со знаком “—” означают, что для данного набора параметров расчеты не проводились из-за того, что не выполнялось неравенство (35).
Приведенные в данной таблице результаты позволяют сделать следующие выводы.
1. Погрешности в расчетах с применением антипараксиальной теории приблизительно в 2 — 10 раз меньше соответствующих погрешностей в расчетах с использованием закона “3/2”.
2. Поведение погрешностей при изменении прикатодного расстояния и числа узлов сетки не является монотонным. Причиной тому может служить сложный характер взаимодействия погрешностей, вносимых на различных этапах решения рассматриваемой самосогласованной задачи [2].
Итерационный процесс по подобластям проводился в два этапа. На первом этапе сходимость с заданной точностью достигалась без учета неоднородности плотности тока вдоль катода, что в данной задаче означает применение закона “3/2” (Ь = 0 в формуле (31)). Это вызвано тем, что на первых итерациях приближения далеки от решения, поэтому вычисление производных, входящих в Ь, дает большие погрешности. Расчеты второго этапа были продолжением расчетов первого этапа при Ь = 0. В табл. 3 приведено суммарное по двум этапам число итераций, которое потребовалось для достижения заданной точности; в скобках указано число итераций первого этапа.
Таблица 2. Погрешности расчета потока Мельтцера, %
а 5<р 8
й\ ^2 4 й\ ^2 СІ3 й\ ^2 4
14 х 8 3.1 1.4 - 4.1 11 — 3.4 0.63 —
13 1.2 — 23 17 — 9.8 1.7 —
28 х 16 3.5 1.8 1.7 2.8 2.2 8.6 4.2 1.3 0.33
13 4.8 1.5 17 7.8 10 11 3.8 0.99
56 х 32 8.9 0.76 1.1 3.1 0.79 1.3 4.4 1.7 0.51
13 4.4 2.5 15 5.7 3.9 11 4.1 1.3
Таблица 3. Сходимость итераций по подобластям
а (1\ ^2 dз
14 х 8 С М П С М П С М П
6(3) 10(7) 13(10) 6(3) 9(6) 13(10) — — —
28 х 16 6(3) 12(8) 21(18) 6(3) 10(7) 13(10) 5(3) 8(6) 18(16)
56 х 32 16(8) 17(11) 17(11) 7(4) 11(8) 14(11) 6(4) 9(7) 12(10)
Вычисления выполнялись по различным итерационным алгоритмам, рассмотренным в разд. 3. Среди них: П — метод простой итерации (19)—(21), М — модифицированный метод простой итерации (23), (20), (21). Проводились также численные эксперименты по методу секущих (29), которые показали, что он хотя и быстро уменьшает невязку на первых приближениях, но является неустойчивым, зачастую “разбалтывается” при приближении к решению. В связи с этим эксперименты, результаты которых приведены в табл. 3, проводились по методу С — смешанному алгоритму, суть которого состоит в следующем. Первые итерации выполнялись по методу секущих (29). Как только норма невязки начинала возрастать, дальнейшие итерации осуществлялись по модифицированному методу простой итерации. Такая процедура проводилась на обоих этапах расчета.
Из данной таблицы следует, что наиболее быстрым является смешанный метод, который сходится с заданной точностью в большинстве случаев за 5-6 итераций, затем следует модифицированный метод простой итерации (9-12 итераций) и, наконец, метод простой итерации (13-17 итераций).
Список литературы
[1] Молоковский С.И., Сушков А.Д. Интенсивные электронные и ионные пучки. М.: Энер-гоатомиздат, 1991.
[2] Ильин В.П. Численные методы решения задач электрофизики. М.: Наука, 1985.
[3] Григорьев Ю.Н., Вшивков В.А., Федорук М.П. Численное моделирование методами частиц-в-ячейках. Новосибирск: Изд-во СО РАН, 2004.
[4] Акимов П.И„ ОсиповА Г.П., Сыровой В.А. Проблемы повышения точности программ траекторного анализа интенсивных электронных пучков // Журн. вычисл. математики и мат. физики. 1989. Т. 29, № 3. С. 405-422.
[5] Свешников В.М. Численный расчет пучков заряженных частиц на локально модифицированных сетках. Новосибирск, 1997 (Препр. РАН. Сиб. отд-ние. № 1189).
[6] Сыровой В.А. Результаты теории антипараксиальных разложений в оптике плотных электронных пучков // Радиотехника и электроника. 1991. Т. 36, № 3. С. 540-559.
[7] Свешников В.М., Сыровой В.А. О численном расчете пучков заряженных частиц методом итераций по подобластям // Журн. вычисл. математики и мат. физики. 1990. Т. 30, № 11. С. 1675-1688.
[8] Свешников В.М. Повышение точности расчета интенсивных пучков заряженных частиц // Прикл. физика. 2004. № 1. С. 55-65.
[9] Свешников В.М. О расчете интенсивных пучков заряженных частиц методом итераций по подобластям без налегания // Прикл. физика. 2006. № 3. С. 49-57.
[10] Самарский А.А. Теория разностных схем. М.: Наука, 1977.
[11] Головин Г.Т. Комбинированный метод решения двумерных стационарных самосогласованных задач // Журн. вычисл. математики и мат. физики. 1987. Т. 27, № 5. С. 700-710.
[12] Мокин Ю.И. О двух моделях стационарного движения заряженных частиц в вакуумном диоде // Матем. сб. 1978. Т. 106 (148), № 2 (6). С. 234-264.
[13] Свешников В.М. О численном решении самосогласованной задачи расчета стационарных пучков заряженных частиц: Новосибирск, 1988 (Препр. РАН. Сиб. отд-ние. № 789).
[14] Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: БИНОМ. Лаборатория знаний, 2006.
[15] Свешников В.М. Расчет плотности тока эмиссии и объемного заряда в задачах численного моделирования электронно-оптических систем // Матер. Междунар. конф. по вычислительной математике (МКВМ-2004): Новосибирск, 2004. С. 642-647.
[16] Свешников В.М. Расчет прикатодной области в электронно-оптических системах, формирующих интенсивные пучки заряженных частиц // Прикл. физика. 2004. № 1. С. 50-55.
[17] Свешников В.М. Поэлементная технология расчета интенсивных пучков заряженных частиц // Вычисл. технологии. 2004. Т. 9, № 3. С. 90-103.
[18] Meltzer B. Single-component stationary electron flow under space-charge conditions // J. Electron. 1956. Vol. 2, N 2. P. 118-127.
Поступила в редакцию 27 января 2005 г., в переработанном виде — 8 июня 2006 г.
