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

ISSN08б8-588б

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2GG2, том 12, № 1, с. 4б-52

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

УДК53.082.72: 621.3.032.26 © С. И. Шевченко

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

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

ВВЕДЕНИЕ

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

ФМ = Jff(r)G(r, ro )dS .

(1)

(S)

где интегрирование проводится по всей граничной поверхности « = ^ « , « — поверхность 1 отрезка границы, 7(г) — функция плотности поверхностного заряда (ППЗ), С(г, г0) — ядро интегрального уравнения.

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

Как было показано в [3], если граница области (электроды) имеет концы или изломы, то функция ППЗ 7 представима в виде

a(l) = ■

~(l)

(2)

7І I1 -d

где 7 (I) — некоторая гладкая функция от I, называемая иногда модифицированной функцией плотности поверхностного заряда (МФППЗ), I — расстояние от начала отрезка до точки интегрирования, й — длина отрезка границы (ОГ), вдоль которого проводится интегрирование, к1 и к2 — показатели сингулярности поверхностного заряда в точках концов и изломов электродов, определяемые по формуле

к=

п-т

2п-т

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

Подобная проблема уже была решена в [4] для случая плоской симметрии геометрии.

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

Ф(го) = ^Рр (гоК

р

где рр (г0) — вклад в потенциал в ТН г0 от отрезка границы с номером р .

Из (1) и (2) легко получаем, что вклад рр (г0) определяется выражением:

Ф

і

p(ro) = dJ<~(Z) G(ro,r)

dZ

ZK1(1 -Z)

к2

(3)

где £ = I / й .

В рассматриваемом случае аксиальной симметрии ядро интегрального уравнения С(г, г0) имеет вид [3]

G(r, ro) =

4r

V(r + ro)2 + (z - zo)2

:K(k ),

(4)

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

Следуя работе [1], можно получить для ядра С(г, г0) выражение

в(г, О = / + Л1п[(2)], где / и /2 — гладкие при (г, г0) > 0, г = г0 = 0, г Ф г0 функции, имеющие следующий вид:

/ =

4г

/2 = -

Здесь

л/(г + Г))2 + (г - ^)2 х(1 + (2 1п[( Г + Г))2 + ( 4г

7( г + Г))2 + (г - г 0 )2

С 4 Л

&1 = £ агт1 , § 2 =

1г=0 ;

Г 4 £».

г =0

до точки интегрирования.

(2)

X

СК1(1 -С)"1

і

+ й |<г(С) /2(С)1и[(С-Со)2 +5і2 ]х

0

X-------------.

£4(1 -£)К

(5)

Видно, что в первом интеграле содержатся две точки сингулярности £ = 0,1, а во втором:

£ = 0, С0, 1, если 8 = 0, т.е. точка наблюдения лежит на отрезке интегрирования (ОИ), или две точки сингулярности и одна негладкость (резкий пик), если 8 Ф 0, но 8 << 1 (ТН лежит вблизи ОИ).

Как и в работе [4], перейдем от численного интегрирования по формуле Гаусса к более общему случаю численного интегрирования по формуле типа Гаусса, когда сингулярности концов ОИ включаются в весовую функцию

Р =

1

ГЧ1 -о

К2

(6)

а Ь(2) = д/(г - г0)2 + (г - г0)2 — расстояние от ТН

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

Ь(2) = (I -10)2 +52,

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

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

гладкая линия, то, как показано в [4], функция Ьт> может быть представлена в виде

Ь(2) =[(С-С0)2 +52 ]р(С),

где 51 =5 / й, ^0 = 10 / й , й —длина ОГ. Функция Р (£) является гладкой на отрезке £ = [0,1].

Полученные выражения подставляем в (3) и после замены переменной интегрирования £ = I / й получаем для вклада в потенциал:

Фр (г0) =

1

= й |<г(С) [/1(0 + /2Іп( Р(д))]

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

первой — формирование матрицы и нахождение функции ППЗ, или на второй — вычисление искомой функции (потенциала) или ее производных

Е(г), Е(г) (компонент напряженности электростатического поля).

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

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

X

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

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

Для вклада в потенциал от одного ОИ имеем выражение (5), в котором функции, ответственные за сингулярности ППЗ, связанные с концами отрезка интегрирования, включены в "вес" правила численного интегрирования (6). Поэтому первый интеграл в выражении (5), содержащий гладкую интегрируемую функцию, просто представляем в виде квадратурной суммы типа Гаусса, а второй требует еще некоторых преобразований, описанных в [7] и названных процедурой понижения степени особенности интегрируемой функции:

Фр (Г0) =

=й £‘к<~(Ск ) х

к=1

х/[) + МСк )1п| Р((к)]+

1

+ й |<г(С) /2(С)1п[(С-С0)2 +512 ]х

х-

-, (7)

Ск1(1 -С)к2

где Лк — коэффициенты квадратурной формулы, ^к — абсциссы квадратурной формулы.

Со вторым интегралом в (7), следуя рекомендациям [4], проведем следующее преобразование:

|[)/2 - {{(С)/2 } ]х

х 1п[-£0? + 512 ]

СК1(1 -о

К2

х

х 1п[-С0)2 +512 ]

(8)

т 41-СГ2’

где {?(£) /2 } — интерполяционный полином от функции <~(С)/2(0, построенный на шаблоне с N.. точками

{~(С) /2(0} =

£а(С, ) /2(Сг )£ (£-£,)

І -1

(9)

1 =1 ] =1

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

всех точках шаблона интерполяции. Хотя интегрируемая функция и не удовлетворяет условию применимости правила численного интегрирования типа Гаусса [7] (гладкость порядка 2^ -1),

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

В результате вычисления первого интеграла в (8) получаем:

^ 21 =

N2

= Е Лк в((к М(, )1п[к-С0)! +82 ]-

кет(^0)

N1 N1

(10)

1=1 ]=1

где т(^0) — массив (набор) номеров точек кол-

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

N г т

«2л = Е Лк (Ск-С0)'-Чп[(Ск-С0)2 +82 ].

к<£т(^0)

Второй интеграл в (8), после подстановки в него (9), принимает вид

Мг

Мг

3 22 = £ ~(Сг) МСг )£ 22 ( І),

г=1 ] =1

(11)

где

1

І (С-С0); -11п[(С-С0)2 + 52 ]■

ас

СК1(1 -о

К2

Стоящий в выражении (11) интеграл 322(І) не относится к числу аналитически берущихся, а прямому применению формул численного интегрирования Гаусса или типа Гаусса препятствует присутствующая под знаком интеграла логарифмическая сингулярность (при 5 = 0, т.е. когда ТН лежит на ОИ) или негладкость (острый пик, при 5 << й , т.е. когда ТН лежит вблизи ОИ).

Методика обращения с этим интегралом была рассмотрена в работе [4], посвященной решению уравнения Лапласа в случае плоской симметрии. Поэтому мы на нем останавливаться не будем.

Таким образом, для вклада в потенциал от одного отрезка границы (7) получаем:

Фр (Го) =

= й £лк~({к)[ (Ск) + / (Ск)]р(Ск)| +

кет(Со)

+ а £лк~(Ск)вк'„

к<£т(С0)

N N г ,

- а £ ~(£,) ^ (С, )£ с,у [ 2 ([) - з22 (])]. (12)

1=1 У=1

Если в правой части соотношения (12) выделить величины ст(С,), то стоящие при них выражения будут являться коэффициентами матрицы. Вид этих коэффициентов в значительной мере зависит от того, принадлежит ли точка — проекция ТН (или сама ТН, если 8 = 0 ) шаблону интерполяции или нет:

1. к е т(Со):Ст(к,п) = аАкОк п,

2. к е т(£о): Ст (к,п) =

= йАк [(£к) + f2{Zk)]1пР(Ск) +

N Г 1

,2<Ск )£ ^ [у202 ( І) - 5 2 ( Л\.

І =1

Полученное матричное уравнение (систему линейных уравнений) решаем методом исключения Гаусса с выбором ведущего элемента. В результате получаем значения "модифицированной" плотности поверхностного заряда (МППЗ) ~(С,)

в точках коллокации.

ВЫЧИСЛЕНИЕ ПОТЕНЦИАЛА ЭЛЕКТРОСТАТИЧЕСКОГО ПОЛЯ ф(г0)

В этом разделе мы приведем описание второй стадии (этапа) решения уравнения Лапласа методом коллокации для соответствующего граничного интегрального уравнения: функция МППЗ

~(С,) уже известна, и нам требуется вычислить значение потенциала в произвольной точке наблюдения ф(г0). Для этого используем выражение (3). Качество функции (гладкость или сингуляр-ность/негладкость), находящейся под знаком интеграла, зависит от того, далеко находится ТН от ОИ или близко от него (или на нем).

ТН находится далеко от ОИ

В этом случае логарифмическую функцию, стоящую в (5) под знаком интеграла, можно считать достаточно гладкой для проведения численного интегрирования по формуле типа Гаусса [7]:

Ф(гс) = а £ Ак <~(Ск )ок'„, (13)

к=1

где °к, п = °(гк , ^ Г0 = Гк .

ТН расположена близко от ОИ

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

Фр (Го) =

ро 1

ас

= й I <7(0 [/ + /24 рЩК1 (1 _0'2 + й}[<7(С)/2(0 _ {<7(0ЛЮ}, ]х

х 1и[(С_Со)2 +5і2І

ас

СК1(1 _С)

к2

і

/{(О /2(0}

х 4с-Со)2+512)

ас

СК1(1 _СГ2

(14)

В отличие от предыдущего параграфа в данном функцию МППЗ ~(С,) считаем уже известной, поэтому все квадратурные суммы вычисляем до конца.

Под знаком первого интеграла в правой части (14) стоит гладкая интегрируемая функция (считаем весовую функцию применяемого правила численного интегрирования типа Гаусса, равной (6)). Поэтому реализация для этой функции численного интегрирования по указанному правилу вполне законна (оправдана).

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

N

+

+

2Ыг — 1. Т.е. применение правила численного интегрирования типа Гаусса ко второму интегралу не является вполне обоснованным, однако значительно повышает точность по сравнению с его применением в выражении (5).

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

В результате всего получим для вклада в потенциал:

ФР (Го) = ^0 + dSl +

Nі Ы,

+ d£&(С,)МС, )£{(7) — ^(;)}, (15)

7=1

где

s 0 =

= £ лк д(Ск)[/ (С к) + / (Ск )]1п| Р(Ск)| +

кеш((0)

+ £ Лк~(Ск )Ок'„,

к«т(Со)

а величины $2( у) и Т2( у) определены выше.

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

ПОЛЯ Е(2), Е(Г)

При вычислении компонент напряженности электростатического поля используем выражение (3), в которое подставлено соответствующее ядро О(2) или О(Г):

О(г )(г, Го) =

4г

д/(г + Го)2 + (^ — г о)2

:Х

г — гп

О( г)(г, Го) =

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

2г 1

-Е(к),

Х

г — г

л/(г + Го)2 + (г — го)

Е(к) — К(к)

г — г,

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

после некоторых преобразований (см. [1]) дают для ядер следующие выражения:

( ^ )

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

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

г — гп

3 (г — Го)2 + (г — го)2’

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

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

(16)

. г(Г) .

Г—Г

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

(17)

где все коэффициенты выписаны в работе [1].

Рассмотрим подробно получение компоненты Е(2). С ядром (16) выражение для вклада в Е(2) от одного отрезка границы имеет вид:

Е «(Го) = d | <7(Ск )[/1( г) + /2г) 1п|р(С)|]х

Х

ас

СК1 (1 — С)

К2

1

+ d | <7(0 /2г )(С)1и[(С—Со)2 +5,2 ]х

Х

ас

СК1(1 — СГ2

+ d | <7(0 /з( г )(С)

г —

1

(С —Со)2 +5,2 Р(С)

ас

-Х

Х-

СК1(1 — СГ2

(18)

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

Как обычно, для функций К(к) и Е(к) используем интерполяционные представления [8], которые

Как обращаться с первым и вторым интегралами мы уже разобрали выше. Интеграл, похожий на третий, был рассмотрен в работе [4] (вместо функции <г(0(х - х0)/Р(0 в данной работе стоит

функция <7(0/3(2)(0(2 - г0)/Р(О). Все остальные действия не отличаются от действий, проведенных выше, или от действий, проведенных в работе [4]. В результате получим интеграл, вычисление которого уже разобрано в [1]. Обозначим его

Т (2 )

3

Таким образом, при вычислении компоненты Е(2)(2) у нас не остается неразрешимых проблем.

Для получения компоненты Е(г) напряженности электростатического поля подход полностью аналогичен.

+

+

Х

Х

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

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

— на первом этапе решения интеграл заменяем квадратурной суммой по правилу численного интегрирования типа Гаусса; коэффициенты в полученной сумме объявляем коэффициентами матрицы; матричное уравнение решаем методом исключения Гаусса с выбором ведущего элемента;

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

(<р(Го), Е(2),Е(г)).

Из вышесказанного очевидно, что неточности вычислений могут возникать при задании геометрии системы, при замене интегралов на квадратурные суммы, особенно, когда ТН близка к ОИ, или при решении матричного уравнения.

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

Решение матричного уравнения методом Гаусса имеет много дополнительных методов уточнения результатов [9], которые позволяют значительно повысить точность решения матричного уравнения. Это дает право считать этап решения матричного уравнения более точным, чем все остальные.

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

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

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

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

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

1. Шевченко С.И. Алгоритм получения высокой точности в расчетах аксиально-симметричных электростатических полей // Здесь. С. 35-39.

2. Соболев С.Л. Уравнения математической физики. М.: Наука, 1966. 443 с.

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

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

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

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

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

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

9. Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических вычислений. М.: Мир, 1980. 279 с.

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

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

Материал поступил в редакцию 24.12.2001.

AXIALLY-SYMMETRIC ELECTROSTATIC FIELD CALCULATIONS TAKING INTO ACCOUNT SURFACE CHARGE DENSITY SINGULARIES

S. I. Shevchenko

Institute for Analytical Instrumentation RAS, Saint-Petersburg

The sources of accuracy losses in electrostatic calculations of electron- and ion-optical instruments with axial symmetry and their parts are analyzed. An algorithm which takes into account these sources and reduces their contribution to computational errors is offered. The singularities of the surface charge density at the electrode ends and bends are taken into account.