Об особенностях нахождения аксиальных электростатических полей вблизи оси. Ii. Метод аналитической замены Текст научной статьи по специальности «Математика»

ISSN 0868-5886

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2009, том 19, № 1, c. 58-66

МАТЕМАТИЧЕСКИЕ МОДЕЛИ

УДК517.956.225: 621.319.7

© С. И. Шевченко

ОБ ОСОБЕННОСТЯХ НАХОЖДЕНИЯ АКСИАЛЬНЫХ ЭЛЕКТРОСТАТИЧЕСКИХ ПОЛЕЙ ВБЛИЗИ ОСИ. II. МЕТОД АНАЛИТИЧЕСКОЙ ЗАМЕНЫ

Приведен алгоритм решения проблемы оси в решении аксиального уравнения Лапласа применительно к методу граничных интегральных уравнений, в котором интегральное уравнение с помощью метода коллока-ции и метода аналитической замены и последующего аналитического интегрирования трансформируется в матричное уравнение. Разработан алгоритм вычисления потенциала и его производных на оси до четвертого порядка. (1-я статья серии — ж. "Научное приборостроение", 2007 г., т. 17, № 1, с. 83-90).

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

ВВЕДЕНИЕ

В процессе работы с пакетами прикладных программ нахождения электростатического поля в аксиально-симметричных системах (электродов) было замечено [1, 2], что точность вычислений вблизи оси и на оси значительно ниже, чем в остальном далеком от оси пространстве. Это, по-видимому, связано с сингулярностью аксиального оператора (уравнения) Лапласа на оси.

Алгоритм решения уравнения Лапласа методом граничных интегральных уравнений описан, например, в [1, 2]. В основе этого алгоритма лежит переход от дифференциального уравнения Лапласа к эквивалентному ему интегральному уравнению [3]:

0(го) = | ^(г ) G(г ,Го)^, (1)

(«)

где <г(г) — функция плотности поверхностного заряда (ППЗ); G(r, г0) — ядро интегрального уравнения; г0 — точка, в которой ищется потенциал (точка наблюдения, ТН); г — точка интегрирования (ТИ); интеграл в правой части берется по площади поверхности всех электродов (границы); коэффициент 1/4п0, который должен быть перед интегралом в правой части, включен в функцию <г(г).

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

В результате выполнения первого этапа решения вблизи оси получается скачок ППЗ, достигающий 1.5-2 % от теоретических значений. Эти результаты были получены при вычислениях с помощью известного комплекса программ "Ро18-80п2" [1], а также с помощью программы, описанной в работе [2]. Ясно, что этот скачок является следствием вычислительных, а не физических причин. Это означает, что если в системе электродов присутствуют части, близкие или касающиеся оси, то без решения проблемы оси на первом этапе решения получается неточность в определении ППЗ, что не позволяет и на втором этапе решения получить приличную точность в вычислении требуемых величин как вблизи, так и вдали от оси. Поэтому существует необходимость внести в алгоритм решения уравнения Лапласа некоторые поправки, учитывающие влияние особенности оператора Лапласа (или соответствующего ему ядра G интегрального уравнения) вблизи оси.

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

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

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

гральных уравнений. В стоящем в правой части (1) интеграле проводится замена гладкой части подынтегральной функции на ее интерполяционное представление и далее осуществляется аналитическое интегрирование.

ОСОБЕННОСТИ АЛГОРИТМА УЧЕТА ЭФФЕКТА ОСИ

Ядро для задачи Дирихле для уравнения Лапласа в случае аксиальной геометрии может быть представлено в виде [6]

О (г, г0) = Л( 2, г) +

+ /2(2, Г)1П[(Г - Го)2 + (2 - 2о)2],

(2)

где функции /, /2 являются гладкими везде в области г Ф 0.

Ясно, что сингулярность дифференциального уравнения при переходе к интегральному уравнению трансформируется в некоторую особенность ядра О(г, г0) этого интегрального уравнения.

Входящие в ядро параметры с, т1 имеют некоторые особенности вблизи оси. Для изучения этих особенностей рассмотрим частный случай геометрии — прямолинейный электрод, перпендикулярный оси и опирающийся на нее своим началом. Если ТН принадлежит поверхности электрода, по контуру которого идет интегрирование, то для параметров с, т можно записать простые выражения:

особенности (функция с — полюс первого порядка в точке г = —г0, а функция т — полюс второго порядка в точке г = —г0).

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

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

1) на первом этапе функция ППЗ еще не известна, а на втором уже известна;

2) ТН на первом этапе располагается в ТК, а на втором этапе — в любой точке оси.

Существует и общее в этих двух этапах. Это главным образом то, что на обоих этапах приходится вычислять вклад в потенциал от каждого элемента границы (отрезка электрода).

В заключение данного пункта отметим, что когда ТН лежит на оси, то функция Грина для потенциала и его производных до четвертого порядка определяется выражениями [1]:

О(0) = 2пг-

1

Е'

(3)

4г

(-

с = -

г + г0

т

у

V г + г0 ]

Графики функций с, т1 при г0 = 0 представлены на рис. 1.2. работы [5]. Существуют некоторые не очень существенные отличия от того, что рассматривается в данной работе, которые происходят из того факта, что в [5] точки коллокации (ТК) расположены в точках-абсциссах гауссового интегрирования. В данной работе ТК расположены равномерно вдоль длины рассматриваемого отрезка границы. Поэтому функция с(г) для первой точки коллокации (для случая, когда начало отрезка границы лежит на оси, ядро имеет упрощенный вид, и представление типа (2) для него не требуется) является константой (г0 = 0, с = 4), а графики функции с(г) для остальных ТК несколько сдвинуты по отношению к графикам работы [5]. То же самое можно сказать относительно графиков функции т (г), которые опираются на ось в других точках, расположенных равномерно. Но эти отличия не меняют главного — эти функции имеют

О(1) = 2пг- 20

Е

3'2

О(2) = 2пг

О(3) = 2пг

,(2 — 20)2

Е

5'2

Е

3'2

15(2 — 20)3 — 9 2 — 20

Е

7'2

Е

5'2

О(4) = = 2пг

105 (2 — 1р)4 — 90-(2 — 20)2

Е

9'2

Е

7'2

+ 9

(2 — 20)

Е5'2

(4)

(5)

(6)

(7)

Здесь Е2 = г2 + (2 — 20 )2 , 20 — 2-координата ТН. Из этих выражений видно, что для вычислений с использованием этих выражений необходимо получить выражения для величин г, 2 — 20 и Е.

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

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2009, том 19, № 1

РЕАЛИЗАЦИЯ АЛГОРИТМА УЧЕТА БЛИЗОСТИ ОСИ

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

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

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

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

Выходом из этой ситуации может послужить идея использовать разные шаблоны интерполяции для разных функций [7], стоящих под знаком интеграла: один для интерполяции функции а с N точками, каждая из которых является точкой кол-локации, второй с ^ = N ■ п точками — шаблон для интерполяции функций /1, /. Конечно, применение такого более частого шаблона на том же самом отрезке должно привести (и приводит) к значительному увеличению времени расчетов, но только для тех точек, которые лежат на или вблизи оси. А таких точек во всей электронно-

оптической системе (ЭОС) единицы. Т. е. на полное время расчетов при формировании матрицы это увеличение числа точек в шаблоне не должно сказаться заметно.

На указанных шаблонах для функции а и для функций /1, /2 формируем сплайны третьего порядка [8] (на каждом отрезке соответствующего шаблона д{ <= д <= д+1 с 7 = 0 + N - 2 или

д <=д<=д+1 с ] = 0 ^ N2 - 2):

Г) = X V (8)

к=0

/1(') = X 41 )дк, (9)

к=0

/') = X Ак )дк, (10)

к=0

где д = I -10; I — расстояние вдоль контура границы от начала отрезка границы до точки интегрирования; 10 — расстояние вдоль контура границы от начала отрезка границы до точки коллока-ции (ТК). Верхний индекс 7 в скобках у функции Г (и далее у других функций) указывает на то, что равенство (8) справедливо на отрезке шаблона с номером 7, а верхний индекс ] — на то, что равенства (9, 10) справедливы на отрезке более частого шаблона с номером ].

Первый этап: формирование и решение матричного уравнения

Особенностью реализации первого этапа решения уравнения Лапласа (формирование и решение матричного уравнения) является то, что ТН помещаются поочередно в ТК, которые расположены на границе (контуре электродов). Эта особенность, возникающая только вблизи оси, состоит в негладком характере параметров с, т1. В этом случае особенности близости к оси могут проявиться (см. предыдущий раздел), только когда ТН находится на отрезке электрода, по контуру которого производится интегрирование и который или опирается одним своим концом на ось или подходит к оси весьма близко. Причем, если ТН помещена в ТК, которая лежит на оси (это первая ТК на данном электроде и электрод опирается своим началом на ось), то ядро имеет вид (3). Если же или ТН помещена во вторую ТК, или отрезок границы не опирается на ось, то ядро имеет общий вид (2), и трудности наибольшие.

Численные эксперименты показали, что учет близости к оси для остальных ТН = ТК рассматриваемого отрезка не является необходимым, однако в данной работе использовалось рассмотрение первых трех ТК от оси.

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2008, том 18, № 4

Если ТН находится на отрезке, не совпадающем с отрезком интегрирования (ОИ), то всегда можно дополнительным делением отрезка интегрирования добиться того, чтобы ТН находилась на достаточном расстоянии от ОИ, сравнимом или большем длины отрезка интегрирования. Поэтому всегда можно применить метод разных шаблонов и осуществить сплайн-представление, и далее провести интегрирование.

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

Отрезок границы — отрезок прямой линии

Если отрезок границы опирается своим началом на ось (в случае конца отрезка все выводы вполне аналогичны) и наклонен к оси под углом }, то для ТН, расположенной в первой ТК, 1о = 0,

г = 181п}, Е = l и ядро (3) имеет вид о(0) = = 2п81п}. Т. е. является константой, и вклад в потенциал от этого отрезка равен

N—2 3

Ф = 2П81п}^^) 4')

1=0 к=0

(11)

где

= 1 | <*=

Л+1 + гк +1 Ь1+1 + Ь'

к +1

А6 = В3 С3 ,

А5 — В3 С 2 + В 2 Сз ,

А4 = В3 С1 + В2 С2 + В1С3, А = В3 С 0 + В2 С1 + В1С2 + В0 С3, А2 = В 2 С 0 + В1 С1 + В00 С 2 , А1 = В1 С0 + В0 С1 , А0 = В0 С0 .

Для получения коэффициентов ^ соотношения (13) подставить

(13)

Ак =

следует

(1)

Вк = ^

Если отрезок границы опирается своим началом на ось, и ТН расположена во второй или третьей ТК, то для ядра уже требуется применять общее выражение (2). В этом случае, как уже описывалось выше, применяем два разных шаблона интерполяции, подставляем выражения (8-10) в правую часть уравнения (1), записанную для одного элементарного отрезка границы, и после некоторых преобразований получаем

d N2 —2 6

Ф = 1О-ад! = ХХС + Ви ),(12)

0 1=0 к=0

где ) = 1-дк 1п||<|, а для коэффициентов

Л( 1) г>( 1)

Ак и Вк существуют простые выражения, определяемые формулами перемножения двух кубичных полиномов:

■,(1) Ск = Ак). А для коэффициентов В^' следует в формулах (13) подставить Ак = В^ ),

В = е (0 С = А 1)

пк - °к > ^к - Л2к ■

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

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

Выписываем коэффициенты матрицы (коэффициенты, стоящие при а и mi в матричной форме записи уравнения (1)). Полученное матричное уравнение решаем методом Гаусса с выделением ведущего элемента.

Отрезок границы — часть дуги окружности

Для ТН, помещенной в первую ТК, принадлежащую отрезку границы, который опирается на ось, ядро (3) легко привести к виду

О(0) = 2п

(

Яа . А} а А}

— 81П--+ 008 «О 008-

Я 2 О 2

л

(14)

где Яс — радиус окружности, описывающей границу; ЯО — г-координата центра этой окружности; }О — угол, под которым из центра окружности видна ТН; « = } +А}— угол, под которым из центра окружности видна ТИ.

Видно, что ядро не имеет особенностей. При подстановке ядра (14) в правую часть (1), записанную для одного отрезка границы, легко получить

N—2 3

Ф = ЬЬА

1=0 к=0

+ 008}

(1)

(15)

к

в

к

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2009, том 19, № 1

где

= I д *

81П-

2Д

¿д, А? = I д

к ео^-^ ¿д.

2Д

Т. е. получили совокупность табличных интегралов.

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

Отрезок границы — гладкая кривая

Используем подход работы [6] и представим функцию Ь в виде Ь = д ■ ¥(д), где ¥(д) — гладкая функция, близкая к единице в окрестности точки д = 0.

В этом случае вклад в потенциал имеет вид

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

Когда ТН расположена на оси и вдали от всех отрезков границы (электродов), то ядро имеет гладкий характер. В этом случае формируем для каждого ядра (3-7) сплайн-представление

О (т) = X О^т) дк

к=0

где т — порядок производной потенциала по переменной г.

Подставляем ядро в выражение для вклада в потенциал от одного отрезка границы:

фТ=X X А"*К

(т,7)

(18)

где коэффициенты А^'1 определяются соотно-

фр = Х((д){/1(д) + /2(д)1п¥(д) + /2(д)1п|д|}, шениями 7(13), в ^которые ^ледует подставить

и целесообразно дополнительно строить сплайн-представление на более частом шаблоне для функции /1 (д) + /2 (д) 1п ¥ (д):

/1(д) + /2(д)1п¥(д) = XАЛ .

(16)

Подставляя (8, 9, 16) в (1), после весьма несложных преобразований получим для вклада в потенциал

N2 -2 6

фр = XX(() #)+Ак) А3)), (17)

j=0 к=0

. (3) . (3)

где интегралы Гк3 , Ак и коэффициенты А1к , А2к уже определены выше.

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

Второй этап: вычисление потенциала и его производных на оси

Вычисление потенциала и его производных на оси имеет свои особенности:

- все ядра имеют относительно простой вид;

- в этих ядрах отсутствуют те особенности вблизи оси, которые принесли нам немалые трудности на первом этапе;

- в отличие от первого этапа решения все ТН лежат на оси и могут совпадать с ТК (или ТИ), только когда отрезок границы, вдоль которого производится интегрирование, опирается одним своим концом на ось.

Ак = АГ', Вк = ^, Ск = Ок .

Когда ТН расположена на оси и вблизи от одного из отрезков границы, то выражение для вклада в потенциал зависит от вида (геометрии) этого отрезка.

Отрезок границы — отрезок прямой линии

Методика подхода к данному виду отрезка границы подобна рассмотренной в [6]. На основе этого подхода легко найти выражения для требуемых величин г, г - г0 и Ь2:

Ь2 = д2 + £2,

г = ад + Ь ,

г - г0=сд+ё,

где ¿у = Ау[72] - Ау[71]; ¿х = Ах[72] - Ах[71]; а = =ду\й; Ь = Ау[71 ] + &х/ё ■10; ё=Ах[71] + ■10 -г0; Ау[71], Ау[72] — у-координаты, а Ах[71], Ах[72]— х-координаты начала и конца отрезка границы.

Подставляем эти величины в выражения для функций Грина (3-7), а последние в правую часть (1), записанную для вклада от одного отрезка границы:

фГ = X X а(0,,)К

(0,<) к

7=0 к=0 N-2 6

фТ = XX а^ ) Ки ),

7= 0 к=0

и так далее до четвертого порядка. Здесь

К

(0,7)

Ъ1+\

= I

д

КГ) = |

д

(д2+*2)

-ад.

к-0

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2008, том 18, № 4

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

Отрезок электрода — часть дуги окружности

Для случая, когда центр окружности не лежит на оси, а ТН лежит на оси, легко выводятся следующие соотношения: r = R0 + Rc sin ú, L =

= 2R sin

Aú

€} = 2^114

f

1=0 k=0

Ro R

iCk) + cosúo Jk

(i )

Л

где cosúo = . 1 -

R

N-2 6

|(m,1) j¿r(m,i)

N2 -2 6

фт) = Z Z4mj) *kmj) •

j=0 k=0

L = Jg2 +S2 F (g),

где функция ^(|) — гладкая функция, равная единице в точке | = 0 .

Строим на шаблоне, в котором точками шаблона являются ТК, сплайн-представления следующих функций, входящих в ядра (3, 4):

ú = ú0 +Aú, Rc sin ú = R0, и функ-

2п

2п

VFg)

= Z Fk(0'¿ )gk

r (2 - Z0)

F 3/2(g)

ция Грина (3) принимает вид (14). Вклад в потенциал от одного отрезка границы

= Z F/V

k-0

Подставляем эти сплайн-представления в ядра (3, 4) и далее в (1) и получаем для вклада в потенциал и его первую производную (т = 0, 1) от одного отрезка границы

N -2 6

yi(m) _ ^^ 17(m,i)m,i)

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

Отрезок электрода — гладкая кривая

Когда ТН находится на оси и вдали от электрода, вдоль контура которого проводится интегрирование, то все ядра (3-7) являются гладкими функциями и для них можно сформировать сплайн-представление, подставить его в (1) и, как в (18), получить

с=цаг к .

1=0 к=0

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

Если расстояние до одного из отрезков границы становится совсем малым и для более частого шаблона число отрезков превышает несколько тысяч, то действуем так, как действовали при формировании матрицы. Применяем методику, описанную в [6], когда функция Е представляется в виде

ФТ = хх^

1=0 к=0

Для производных потенциала порядка выше первого, как уже упоминалось, реального улучшения этим методом достичь не удалось. Поэтому производные высших порядков в пределе 8 ^ 0 рассматривать не будем.

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

ТЕСТОВЫЙ ПРИМЕР И ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

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

1) плоский конденсатор в аксиальной геометрии с электродами, расположенными перпендикулярно оси;

2) сферический конденсатор;

3) гиперболический конденсатор.

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

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

r

k -0

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2009, том 19, № 1

Результаты вычисления потенциала Р0 и его первых четырех производных Р1-Р4 на оси сферического конденсатора в сравнении с теоретическими значениями (каждая вторая строка значений)

х Р0 Р1 Р2 Р3 Р4

0.0000 0.999999994 -13.04472 -446.76468 -16190.993 -622150.84

1.000000000 -23.33333 -466.66666 -14000.000 -560000.0

0.0003 0.992978931 -23.4739659 -470.89111 -14080.896 -4.4419+09

0.992978936 -23.4739658 -470.89199 -14169.267 -568476.132

0.0006 0.985915487 -23.6158735 -475.16802 -14339.500 -1.5570е+07

0.985915493 -23.6158736 -475.16848 -14341.1011 -577106.68

0.0009 0.607476637 -31.8441219 -479.49662 -14515.563 -733500.72

0.607476635 -31.8441203 -479.49691 -14515.5473 -585894.94

0.0144 0.607476637 -31.844121 -744.02147 -26075.483 -1218490.4

0.607476635 -31.844120 -744.02150 -26075.519 -1218482.2

0.0147 0.597889802 -32.068507 -751.89928 -26444.254 -1240068.9

0.597889800 -32.068505 -751.89931 -26444.289 -1240060.4

0.0150 0.588235296 -32.295272 -759.88869 -26819.567 -1262107.8

0.588235294 -32.295271 -759.88873 -26819.602 -1262098.9

0.0153 0.578512398 -32.524452 -767.99170 -27201.562 -1284618.8

0.578512396 -32.524450 -767.99174 -27201.596 -1284609.0

0.0156 0.568720381 -32.756079 -776.21033 -27590.383 -1307613.7

0.568720379 -32.756077 -776.21037 -27590.416 -1307602.6

относительный скачок ППЗ на оси уменьшился до величины, не превышающей 10-6.

При вычислении потенциала электростатического поля и его производных до 4-го порядка были получены результаты, приведенные в таблице. Видно, что точность вычисления потенциала вплоть до поверхности и на поверхности электрода, опирающегося на ось, лежит в пределах 10-7 +10-8.

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

Вторая производная вдали от электродов, опирающихся на ось, имеет точность 6-7 знаков, третья — тоже 6-7 знаков, а четвертая — 5 знаков.

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

Еще остается проблема с выбором пороговых значений параметров перехода от одного алгоритма расчета к другому. При вычислении как на первом, так и на втором этапах решения уравне-

ния Лапласа нами применялось или сплайн-представление всего ядра или выделение в ядре логарифмического члена. Порог этого перехода был определен численными экспериментами. Было показано, что в диапазоне расстояний до отрезка границы, вклад от которого вычисляется, За = (0.01 + 0.2)^ оба алгоритма дают одинаковый

вклад с относительной точностью не хуже 10-8.

Для алгоритма учета близости к оси граничные значения расстояния до оси, при котором следует переходить от неучета этих эффектов к их учету, следует выбирать из вида функций с(х), т1 (х) с точки зрения возможности их представления (и представления функций /1, /2) сплайном третьего порядка. Если рассмотреть случай, когда отрезок границы опирается на ось, то для первой ТК проблем с функциями с(х), т1 (х) не существует, т. к. с = 4, а т1 = 1. Для второй ТК гладкость функции т (х) явно недостаточна для ее интерполирования сплайном третьего порядка на шаблоне точек кол-локации. Однако для более частого шаблона с п = 10 точек шаблона уже явно достаточно для хорошей интерполяции функции т (х).

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2008, том 18, № 4

ЗАКЛЮЧЕНИЕ

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

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

Тестовые вычисления, проведенные для нескольких геометрий электродов, показали, что описанный алгоритм является работоспособным и позволяет получить в результате решения уравнения Лапласа относительную величину пика функции ППЗ на оси не более 10-6, а относительную точность вычисления потенциала электростатического поля и его первой производной на оси получить в диапазоне 10-7 +10-8.

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

Приведенный в данной работе алгоритм был реализован в виде программ, включенных в очередную версию пакета прикладных программ "Shift" [9].

СПИСОК ЛИТЕРАТУРЫ

сибирск: Изд-во СО АН СССР, 1986. 193 с.

2. Шевченко С.И. Алгоритм получения предельной точности в электростатических расчетах элементов электронно- и ионно-оптических приборов, имеющих плоскую симметрию // Научное приборостроение. 1997. Т. 7, № 1-2. С. 45-53.

3. Соболев С.Л. Уравнения математической физики. М., 1966. 400 с.

4. Монастырский М.А., Иванов В.Я., Куликов Ю.В. Устранение особенностей в методе интегральных уравнений при расчете осевых распределений потенциала и его производных вблизи границы // Новые методы расчета электронно-оптических систем. М.: Наука, 1983. С.187-191.

5. Шевченко С.И. Об особенностях нахождения аксиальных электростатических полей вблизи оси. I. Метод прямого интегрирования // Научное приборостроение. 2007. Т. 17, № 1. С. 8390.

6. Шевченко С.И. Алгоритм получения высокой точности в расчетах аксиально-симметричных электростатических полей // Научное приборостроение. 2002. Т. 12, № 1. С. 40-45.

7. Шевченко С.И. I. Метод аналитической замены в задаче вычисления аксиальных электростатических полей // Здесь. С. 00-00.

8. Стечкин С.Б., Субботин Ю.Н. Сплайны в вычислительной математике. М.: "Наука", 1976. 248 с.

9. Шевченко С.И. Пакет прикладных программ "Shift" для решения уравнений Лапласа и Пуассона // Материалы Всес. семинара "Методы расчета электронно-оптических систем". Алма-Ата, 1992. С. 11.

Институт аналитического приборостроения РАН, Санкт-Петербург

1 Ивтж КЯ Методы автоматизир°ванног° Материал поступил в редакцию 13.10.2008. проектирования приборов электроники. Ново-

ON SOME PECULIARITIES ARISING IN CALCULATING NEAR-AXIS AXIAL ELECTROSTATIC FIELDS. II. ANALYTICAL REPLACEMENT METHOD

S. I. Shevchenko

Institute for Analytical Instrumentation RAS, Saint-Petersburg

This paper presents the algorithm of resolving the axial problem in axial Laplace equation solving by boundary integral equation method where the integral equation is transformed into matrix equation by collocation,

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2009, том 19, № 1

analytical replacement, and analytical integrating. The method of potential and its derivatives before forth order at axis calculation have been developed.

Keywords: electrostatic field potential, Laplace equation, axial problem, boundary integral equation method, collocation method, method of analytical replacement.

HAYHHQE nPHBQPQCTPQEHHE, 2008, tom 18, № 4