Алгоритм получения высокой точности в расчетах аксиально-симметричных электростатических полей Текст научной статьи по специальности «Математика»

ISSNÜS6S-5SS6

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2ÜÜ2, том І2, № l, c. 4Ü-45

ОРИГИНАЛЬНЫЕ СТАТЬИ

УДК53.082.72: 621.3.032.26

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

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

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

ВВЕДЕНИЕ

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

ф(го) = гоМ5, (1)

( 5 )

где интегрирование проводится по всей граничной поверхности 5, г е 5, 7(г) — функция плот-

ности поверхностного заряда (ППЗ), С(г, г0) — ядро интегрального уравнения.

Точка г0, в которой ищется потенциал, называется точкой наблюдения (ТН).

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

мр мр Лр

ф(Го) = ^фр(Го) = X I°(1 )С(г, ГоМI, (2)

р =1 р =1 о

где суммирование проводится по всем отрезкам электродов, а интегрирование — по длине контура

каждого отрезка электрода, (рр (го) — вклад

в полный потенциал от одного отрезка, Мр — число отрезков границы, йр — длина каждого отрезка (ниже будем обозначать ее просто d).

Решение задачи нахождения потенциала ф(го) и компонент напряженности электрического поля Е(х), Е(г) осуществляется в два этапа.

1. Соотношение (2) рассматривается как интегральное уравнение относительно неизвестной функции ППЗ 7. Это уравнение решается методом коллокации и в результате находится массив значений 7(1) в точках коллокации.

2. Известная функция ППЗ 7(1) подставляется в выражение (2), проводится численное интегрирование и получается значение потенциала в произвольной точке го .

Как в случае плоской геометрии [2], в рассматриваемом случае аксиально-симметричной геометрии в основе как первого, так и второго этапов решения задачи нахождения потенциала ф(го) и компонент напряженности электрического поля

т~'( х) Т'(г) ^

Е , Е лежит взятие интеграла в правой части

(2) численным методом. Как правило, применяется численное интегрирование по формуле Гаусса.

Для применимости формулы Гаусса N порядка необходимо, чтобы интегрируемая функция была непрерывна вместе с ее производными вплоть до 2N — 1 порядка [3]. Очевидно, что основной путь повышения точности расчетов как на первом, так и на втором этапах решения задачи Дирихле — это увеличение степени гладкости функции, подвергаемой численному интегрированию.

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

ЯДРО ИНТЕГРАЛЬНОГО УРАВНЕНИЯ

Отличие рассматриваемого случая аксиальносимметричной геометрии от случая плоской геометрии заключается в более сложном виде ядра интегрального уравнения G(r, го). Это ядро для задачи Дирихле для уравнения Лапласа имеет вид [4]

( 4

G(r, Го) =

4г

->/(г + Го)2 + (^ - ¿о)2

■K(k):

(3)

где го = (го, го) — координаты точки наблюдения, г = (г, г) — координаты точки интегрирования, К(£) — полный эллиптический интеграл первого рода,

к2 = ■

4гг0

(г + Го) + (г - г о)

( 4

к(к) = 2 аіт1 + 2 ь*т1

і =о

і =о

1п

1

т

(4)

где

т,= 1 - т = 1 - к2 = (Г " Го)2 + (г - г")!

(г + Го) + (г — го)

значения коэффициентов а, и Ъ, приведены в [7].

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

б(г, го) = ¡1 + / 1п[(г — Го)2 + (г — го)2 ], (5)

где

/1 =

4г

УІ(г + Го)2 + (г - ¿о)2

х(( + (21п[( Г + Го)2 + (г - го)2 ])

Л = -

4Г

УІ(г + Го)2 + (г - ¿о)2

( 4

2ат , &2= 2 ьгт'і

, і=о ) , і=о ,

Эта задача Дирихле с ядром вида (3) уже рассматривалась в [4]-[6], где с помощью алгоритма, рассмотренного в общем виде в [3], под интегралом в (2) после некоторых преобразований получали непрерывную функцию с разрывной первой производной, которую затем подвергали численному интегрированию по формуле Гаусса. Ясно, что когда уже первая производная имеет разрыв, то ни о какой приличной точности речи быть не может.

Полный эллиптический интеграл первого рода К(£) допускает интерполяцию [7] (о < т < 1)

& =

Выражение /21п[(г — го)2 + (г — го)2] выделено в отдельный член, т. к. именно оно несет в себе сингулярность и при интегрировании должно быть рассмотрено с особым вниманием. Входящие в (5) функции / и /2 являются гладкими везде в области г > о .

Введем обозначение

Ё(2) = (г — Го)2 + (г — го)2.

Как было показано ранее в [2], функция Ё2) для отрезка границы (ОГ) в виде отрезка прямой линии имеет вид

¿(2) = (I — /о)2 + 52,

где 1о — расстояние от начала рассматриваемого отрезка до точки — проекции ТН на данный отрезок, 5 — расстояние от ТН до отрезка.

Если в качестве ОГ используется некоторая

гладкая линия, то, как показано в [2], функция Ё2) представима в виде

Ё(2) = [(С — Со)2 + 51 ]р(С),

где 51 =5 / d , £ = I / d , £о = 1о / d, d — длина ОГ, функция Р (£) является гладкой на отрезке

С= [о,1].

НАХОЖДЕНИЕ ПЛОТНОСТИ ПОВЕРХНОСТНОГО ЗАРЯДА (ФОРМИРОВАНИЕ И РЕШЕНИЕ МАТРИЧНОГО УРАВНЕНИЯ)

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

Точка наблюдения находится вдали от отрезка интегрирования

В этом случае поведение содержащейся в ядре

функции ЁТ> и всего ядра является гладким. Поэтому интеграл в правой части выражения (2) можно взять численно по формуле Гаусса. Выражение для вклада в потенциал от одного отрезка границы имеет вид

<Р

до=<і 2 лмСк ^к,,

(6)

к=1

N

Точка наблюдения находится вблизи или непосредственно лежит на отрезке интегрирования

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

функции и получим в результате выражение Ф(го) = d17(0/1 (О + /2 (01п|Р (С)|] С +

о

1

+ d 17&)/2&)1п[/ —Со)2 + 512 ]с!С- (7)

о

Все находящиеся под знаком первого интеграла функции являются гладкими. Поэтому к первому интегралу можно применить формулу численного интегрирования Гаусса

Ng

3, = d£Лк7(Ск)() + Л(С* )1п|Р(£к )|) . (8)

к=1

Функция, стоящая под знаком второго интеграла, содержит логарифмическую сингулярность, лежащую на (при 5 ^ о) или (при 5 Ф о) вблизи отрезка интегрирования. Применим к функции, содержащейся под знаком второго интеграла в (7), описанную в [3] методику понижения особенности интегрируемой функции, для чего проведем преобразование подынтегральной функции

з2 = d 17(0/2(0 — {(О/2(0} +{(0Л(0} ][/ — Со)2 +512], (9)

о

где

N1 N1

{(0/2=Х7((,)/2((,)Х т,,,(С—(оУ—'

, =1 ]=1

— интерполяционный полином для функции 7(д) /2(д), 7(Сг-) — значения ППЗ в точках

коллокации, принадлежащих шаблону интерполяции, Nj — число точек в шаблоне интерполирования, т, . — коэффициенты.

Эти коэффициенты т{ ■ можно вычислить для всей геометрии один раз при обработке геометрической информации. Значения ППЗ 7(Сг-) считаются неизвестными при нахождении коэффициентов матрицы.

Далее придерживаемся канвы рассуждений, приведенных в [3]: функция, стоящая под знаком интеграла в правой части (9), не является достаточно гладкой, чтобы численное интегрирование по формуле Гаусса было достаточно обоснованным (требуется гладкость вместе с производными до 2N — 1 порядка). Однако в результате проведенных преобразований степень гладкости подынтегральной функции значительно повысилась. Поэтому можно надеяться, что численное интегрирование этой функции по формуле Гаусса даст результат гораздо более точный, чем интегрирование исходной функции.

В результате проведенных преобразований получаем для интеграла 3 2

N1 N1

32 = d5 о +Х7(Ь )Л(Сг )£ти (Ъ(1) — d5( ^ (Ю)

,=1 ]=1

где

N¡1

5 о = ХЛк7(Ск) /21п[(Ск —Со)2 +52 ],

к Фт(г)

1

Ъ™ = dI(С —Со);—11п[(С—Со)2 +512] ,

о

Ng

5« = X А (Ск —Со);—11п[(Ск —Со)2 +512 ].

к Фт(г)

Видно, что Ъ1( 1> является просто табличным интегралом. Теперь, собирая вместе члены, стоящие при различных значениях ППЗ 7(Сг), мы получим коэффициенты матрицы (матричное уравнение относительно неизвестных значений ППЗ 7(Сг-)). Это матричное уравнение (систему линейных алгебраических уравнений) решаем методом исключения Гаусса с выбором ведущего элемента. В результате получаем значение функции ППЗ 7(Сг-) в точках коллокации.

ВЫЧИСЛЕНИЕ ПОТЕНЦИАЛА ЭЛЕКТРОСТАТИЧЕСКОГО ПОЛЯ

Для вычисления потенциала электростатического поля используем формулу (2), в которой функция ППЗ 7(0 считается уже известной.

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

нованно применение в выражении (2) формулы численного интегрирования Гаусса

ФР (Го) = й£ АМСк){/1 + /21п[2)(Ск)]}.

к=1

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

[3], предварительно выделив особенность, содержащуюся под знаком логарифма

О(г)(г, го) =

4г

л/(г + го)2 + (г - г о)2

г - гп

х-

(г - Го) + (г - го)'

-Е(к),

О(Г)(г, го) =

2г

х

г - г.

уі(г + го)2 + (г - го)

Е(к) - К(к)

г - г,

(г - Го) + (г - го)

где Е(к) — полный эллиптический интеграл второго рода.

После некоторых преобразований и подстановки вместо К(к) и Е(к) интерполяционных выражений из[7] получаем:

О(г )(г, Го) =

= /і(2) + /2(^) 1п[(г - го)2 + (г - го)2 ] +

+ /з( г)

г - гп

(г - го) + (г - г о)2

(11)

Фр (Го) = й £ АМСк )[/1 + лЫЖк )]+

к=1

1

+ й¡а(С)Л1п[(С-Со)2 +512 ] .

о

В стоящем в правой части интеграле применяем процедуру понижения степени особенности интегрируемой функции и получаем формулы (8), (10). Далее суммируем вместе все члены, дающие вклад в потенциал, создаваемый в точке наблюдения различными отрезками границы, и получаем потенциал в точке наблюдения.

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

ПОЛЯ Е(г) И Е(Г)

Известно, что выражения для компонент напряженности электростатического поля Е(г) и

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

где

с,( г) =

4г

д/(г + го)2 + (г - го)2

г - г о

х

(г + го) + (г - го) х(833 + 844 1п[(г + го)2 + (г - го)2 ])

К г ) = -с ( г )_____________г го_____________д

/ \2 / ^ 644 ’

(г + го) + (г - го)

/2( г) = г

/з( г) = 4 г),

833 = £ С‘га1-1 , д 44 = £ ат1-1

г=1

Значения коэффициентов с и й приведены в [7].

Аналогичным способом можно вывести выражения для ядра О(Г) :

О(Г )(г, го) = /1( Г) +

+ /2Г) 1п[(г - Го)2 + (г - г о )2 ]+

+ /3( г)

г - го

(г - го)2 + (г - го)2

(12)

где

/1( г) = ^- х

(г)

х{83 -81) + (84 -82)1п[(г + го)2 + (г-го)2]}

х

х

г

(г)

(г) 1

/2( г) =

го I

/3( г) = 2С1( г),

82 - 844

2 2 2 г - го + (г - го)

(г + го)2 + (г - го)2

г) =

2г

4(г + Го)2 + (г - го)2

ё з = 1 + т1 ё 33, & 4 = т1 ё 44.

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

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

+ а ¡0(0 /3г )(С)

г-г

(С-Со)2 + ¿1 Р(С)

¿С- (14)

Рассмотрим подробно вычисление третьего члена в (14). В этом интеграле функция Р (^) является гладкой на отрезке интегрирования, поэтому ее можно включить в совместную с функциями о (С) и /3(0 интерполяцию:

0(0 /3 г)(Р 1 =

. Р(«)

=2 2,и к-го)-

(15)

1=1

После этого применяем к интегралу процедуру понижения степени особенности интегрируемой функции:

^3 =

= а

х

0(0 /3( г )(Р

Р(С)

¿с

(С-Со)2 +5;

0(0 /3( г )(Р'

Р(С)

х

Е рг,г)(го) = й £ АкО(Ск )О(г,г)(Ск ,Со). (13)

к=1

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

Чтобы обойти эту сложность, используем метод понижения степени особенности интегрируемой функции [3]. Записываем выражение для вклада в компоненту напряженности электроста-

7~’( г)

тического поля Е от одного отрезка границы, предварительно использовав соотношение (5):

Е р )(го) =

мё

= й £ АкО((к)(//г) (к) + /2г) (С к) 1п[Р(Ск )])+

к=1

1

+ й |а(0 /2г )(С)1п[(С-Со)2 +82 ] +

+ а ¡{ОММ 1

о [ Р(С)N

¿с

(С-Со)2 +512

и осуществим численное интегрирование первого интеграла в правой части последнего выражения по формуле Гаусса.

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

^ о(С,)/3(г)(Сг )

Р(Сг)

х

М

х2ча ¡ (^-Со); -1

І =1 о

¿с

(С-Со)2 +51

2

Первые два члена, находящиеся в правой части, мы уже рассмотрели выше: первый — см. (8), а второй — (9).

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

ОБСУЖДЕНИЕ ПОЛУЧЕННОГО АЛГОРИТМА

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

N

+

1

ции N — 1. При этом повышение порядка интерполяции повышает точность вычислений.

Подчеркнем, что в данной работе использовались два разных шаблона интерполяции: один — с числом точек N для интерполяции НПЗ, он определял степень гладкости подынтегральной функции; второй — с числом точек N для интерполяции функции L(2). При этом считалось, что

функция L(2) интерполировалась бесконечно точно. На практике, однако, достаточно было, чтобы выполнялось неравенство N > N + 2 .

В данной работе использовались термины "близко" и "достаточно далеко" от элемента границы, по которому идет интегрирование. Вполне очевидно, и на практике подтвердилось, что формулы для близкого расположения точки наблюдения от рассматриваемого отрезка границы справедливы и на значительном удалении. Однако эти формулы и соответствующие им алгоритмы весьма громоздки и требуют значительного времени для своего выполнения. Поэтому, сравнивая результаты применения тех и других формул, можно найти расстояние, начиная с которого можно перейти к более простым формулам. Например, для N = 3 это расстояние приблизительно равно 8 = 0.2d .

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

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

L Соболев СЛ. Уравнения математической фи- Материал поступил в редакцию 24.12.2001.

зики. М.: Наука, 1966. 443 с.

2. Шевченко С.И. Алгоритм получения предельной точности в электростатических расчетах

AN ALGORITHM FOR HIGH-ACCURACY CALCULATIONS OF AXIALLY-SYMMETRIC ELECTROSTATIC FIELDS

S. I. Shevchenko

Institute for Analytical Instrumentation RAS, Saint-Petersburg

The paper presents an algorithm which helps to achieve high accuracy in the calculation of axially-symmetric electrostatic fields because the reasons leading to the loss of accuracy are taken into account. The singularities of the surface charge density at the electrode ends and bends are ignored.

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

3. Крылов В.И., Шульгина Л.Т. Справочная книга по численному интегрированию. М.: Наука, 1976. 370 с.

4. Антоненко О.Ф. Численное решение задачи Дирихле для незамкнутых поверхностей вращения // Вычислительные системы. Новосибирск: Изд-во ИМ СО АН СССР, 1964. № 12. С.39-47.

5. Иванов В.Я. Методы автоматизированного проектирования приборов электроники. Новосибирск: Изд-во СО АН СССР, 1986. 193 с.

6. Тиунов М.А., Фомель Б.М., Яковлев В.П. SAM — интерактивная программа для расчета электронных пушек на мини-ЭВМ. Новосибирск: Препринт 87-35 Института ядерной физики СО АН СССР, 1987. 63 с.

7. Абрамовиц М., Стиган И. Справочник по специальным функциям. М.: Наука, 1979. 830 с.

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

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