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

ISSN Ü868-5886

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2002, том 12, № 2, c. 23-29

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

УДК53.082.72: 621.3.032.26

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

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

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

ВВЕДЕНИЕ

Данная работа является продолжением работы

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

Аф = -

Q

(І)

где £0 = 8.8542 • 10 Ф/м — электрическая постоянная, Q — функция плотности объемного заряда (ПОЗ).

Уравнение (1), имеет эквивалентное ему интегральное уравнение [3], записанное для случая аксиально-симметричной геометрии:

<p(ro) = \а G dL + -—f QG dS,

(L) 4П—0 (S)

(2)

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

Задача Дирихле для аксиально-симметричного уравнения Лапласа рассматривалась в работах [47], где после некоторых преобразований под интегралом в (2) получали непрерывную функцию

с разрывной первой производной, которую затем подвергали численному интегрированию по формуле Гаусса.

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

G(r, ro) =

4r

:K(k ), (3)

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

к2 =

4rr0

(r + ro)2 + (z - z0)

2

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

с точностью не хуже 2 -10-8 [8]

4 4 і

K(k ) = У Ö,.m1 + (У bm1 )ln—

i=o to mi

где

m1(z, r) = 1 - m = 1 - к =

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

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

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

£

o

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

6(г, г0) =

= Мг,г) + Л(г, г)1п[(г - г0)2 + (г - г0)2], (5)

где

/1( г, г) =

4г

-\/(г + Г0)2 + (г - г0)2

:Х

X {^1 (г, г) + g2 (г, г)1п[(г + г,)2 + (г - г0)2]},

Л( гг) = -

4г

у1(г + Г0)2 + (г - ^)2

gl( г г) = 2 °^т1

g2 (^ г) = 2Ьг‘

г=0

Выражение,

содержащее

функцию

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

Как обычно (см., например, [1]), решение уравнения Пуассона осуществляется следующим образом: полный потенциал ф разбивается на две части: вклад от ОЗ ф2 и некоторый потенциал (р1. Если вклад ф2 от ОЗ, удовлетворяющий уравнению Пуассона (1), представить в виде

Ф2 =

(6)

то для оставшейся части потенциала (р1 получаются уравнения (граничная задача)

Аф1 = 0,

^ к- = иг ,

(7)

значение потенциала на этом электроде

чение функции ф2 на поверхности этого электрода.

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

1. По формуле (6) находится вклад ф2 от ОЗ в точках коллокации задачи (7), и в результате подстановки этого вклада во второе равенство (7) находится функция плотности поверхностного заряда (ППЗ) ст(г).

2. Известная функция ППЗ (7 (г) подставляется в

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

Особенности, которые появлялись в первом интеграле, были ранее рассмотрены в работе [7]. В данной работе рассмотрим особенности второго интеграла в (2), т.е. интеграла (6).

Чтобы не проводить сразу интегрирование по всей площади, занимаемой ОЗ, эту площадь покрывают некоторой прямоугольной сеткой (возможно, неравномерной). Тогда интеграл в (6) сводится к совокупности интегралов по каждой отдельной ячейке

Ф2(го , Го ) = 22^ (гг , гг Ь

-г =1 г г =1

(8)

где гг, гГ — номера ячеек (номер левой нижней вершины) вдоль направлений Z и Я соответственно, МГ — количества ячеек сетки в направлениях Z и Я соответственно, Цг, 1Г) — вклад от

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

1 ч гк

(-г,гг) = 4^ 1 ¿г 1 drQ(z,г) в(г,г,г0,Го), (9)

где гп, гк — левый и правый, а гп, Гк — нижний и верхний пределы рассматриваемой ячейки, г0, го — координаты точки наблюдения.

Вполне очевидно, что простейшим способом взятия двойного интеграла в (9) является повторное численное интегрирование. Если, например, применить повторное численное интегрирование методом Гаусса [9], то получаем формулу

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

Здесь Si — поверхность электрода номер г, иг —

— зна-

4п^0

N N

Х 2 Аг 2 А Q(гг , Г1 ) °(гг , Г1 , *о , Го ^

г =1 1 =1

(10)

2„ Г

П 'п

где Q(zi, г.) — значения плотности объемного заряда в узлах интегрирования (хі, у.); 8г = гк - гп

— длина (протяженность) рассматриваемой ячейки в направлении X; 8Г = гк - гп — длина (протяженность) рассматриваемой ячейки в направлении Я; Аі, А. — коэффициенты квадратурной формулы Гаусса.

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

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

В известных нам работах [5, 6] в пределах каждой ячейки функцию ПОЗ считали константой. Очевидно, что такая аппроксимация не может быть удовлетворительной для неоднородных функций распределения ПОЗ.

Если ядро интегрального уравнения О, преобразованное к виду (5), подставить в формулу (9), то получим в последней под знаком интеграла (помимо всего прочего) произведение О•Q. Т. к. функции /1(г, г), _/2(г, г) и Q мы считаем гладкими, то имеет смысл проводить интерполяцию не одной функции ППЗ, а произведений

/(г, гШг,г) и Яг, гШ.г, г) .

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

1 zk rk

-------I dz I dr[d00 + d10 z + d01 r + d11 zr] +

4ne„ j j

fiQ = d 00 + d 10 z + d 01Г + d„ zr, f2Q = C00 + C10z + C01r + C11zr-

(11)

В этой формуле уже проведено приведение координат к ТН (см. [1]).

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

Из предыдущего материала следует, что рассматриваемый вклад в потенциал от прямоугольной ячейки имеет вид

+--------X

4п^0

zk rk

X |dz|dr[c00 + c10 z + c01 r + c11 zr]ln[z2 + r2]. (12)

“n n

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

Первый интеграл равен сумме

JI — — \dgoР00 + рю + doiPoi + dnPn],

где P00, P10, P01, P11 — элементарно берущиеся интегралы:

P00 =Az • Ar .

P01 = Az

2 2 rk - rn

P10 1 z 2 - z2 • Ar ,

2

P11 = zk2 - z2 r2 -

2 2

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

ПОЛЯ Е(г) ОТ ПРЯМОУГОЛЬНОЙ ЯЧЕЙКИ, ЗАПОЛНЕННОЙ ОБЪЕМНЫМ ЗАРЯДОМ

Вклад в компоненту напряженности электро-

7~’( г) и и

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

E(q z) = | QG(z }dV

(13)

где ядро [6]

G(z)( z, r) =

4r

(z - z0)

-E(k).

•^(г + г0)2 + (г - г0)2 (г - г0)2 + (г - г0)

Б(£) — полный эллиптический интеграл второго рода в области 0 < т1 < 1 допускает интерполяцию с точностью не хуже 2 -10 8 [8]:

nn

2

4 4 1

е(К) = 2 cm1 + (2 dm\ )ln—• 1=0 7=T m

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

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

О(г) =

= /')(г, г) + /2=>( 2, г )1п[( г - 2 0)2 + (г - г,)2] +

/2 ') = -с( ')

/3 ') =

(r + r0) + (z - z0)2

+ /з( z )( z, r )

z - z

(z - zo) + (r - Г0)2

(14)

где

ciz ) =

4r

V(r + Г0)2 + (z - zo)2

z - z 0

-X

f ( z ) = c ( z ) ___________________

J1 1 / \2 Í \2 ( r + r0) + (z - z0)

X(g33 + g44 ln[(r + r0)2 + (z - z0)2 ])

g 33 = 2 cm1-1, g 44 = 2 dm1-1 •

i=0

i=1

Обратим внимание, что функции

/1( z), /2(z ), /3( z) и Q являются гладкими, поэтому

имеет смысл проводить совместную интерполяцию следующих функций

/1( z )Q = d 00) + d^) z + d 0z) r + d(¡) zr,

f2 z)Q = c00} + c1(0)z + c0? r + СП )zr,

/( z )Q = p00) + PÍ0) z + p0z ) r + ) °r-

(H этих формулах уже сделано приведение координат к ТН)

При таком подходе вклад в компоненту El'0) принимает вид

л “П 'k

E[z) =-í dz í dr[d 00) + d1(0) z + d 0°) r + d1(f) zr ] +

4 4Я£0 J J

0 zn rn

1 zn rk

-----í dz í dr[c00) + c|0)z + c0z)r + с!!)zr]ln[z2 + r2] +

0 zn rn

1 zn rk

+-— Jdz Jdr[ p00 )+P1(0z) z+p0z ) r+píf zr]

22 z + r

(15)

°n ’n

Первый из стоящих в правой части последнего выражения интегралов рассмотрен нами в предыдущем параграфе. Второй и третий интегралы были нами рассмотрены в работе [1].

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

ПОЛЯ Е(г) ОТ ПРЯМОУГОЛЬНОЙ ЯЧЕЙКИ, ЗАПОЛНЕННОЙ ПРОСТРАНСТВЕННЫМ ЗАРЯДОМ

G( r)(z, r) =

2r / r0

V(r + r0)2 + (z - z0)2

'r2 - r02+ ( z - z 0)2

(r - r0)2 + (z - z0)2

X

X

E(k) - K(k)

Если провести с ядром G(r) такие же преобра-

/'"’(z)

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

Г^(г) ~ ~

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

К г) =

E(r = JQG(r) dV,

следующее выражение G(r) =

= /1(r)(z,r) + /(r) (z, r)ln[(z - z0)2 + (r - r0)2] +

(16)

+ /3 r)( z, r)

где ядро [6]

r - r0

(z - z0) + (r - r0)

(17)

где

fír) = с ( r )

— Í&3 - (О + ({4 - g2)ln[(r + r0)2 + (z - Z0)2 ] ro

(r)

g 2 - g 44

2 2 2 r - r0 + (z - zo)

(r + ro)2 + (z - Zo)2

f3r) = 2ci(

C (r) = 1

(r)

2r

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

g3 = 1 + mi g33. g4 = mi g44.

Обратим внимание, что функции f1(r), f2( r), f3(r) и Q являются гладкими, поэтому имеет смысл проводить совместную интерполяцию следующих функций:

fi( r )Q = d oo} + di(or) z + d or) r + d(iri) zr,

f2 r)Q = coo) + C(o)z + Co()r + C(()zr , f((r)Q = ^oo) + Рш)z + po() r + Pn)zr .

(В этих формулах уже выполнено приведение координат к ТН.)

При таком подходе вклад в компоненту E(r) принимает вид

zn rk

Е(r) = q 4ле

л “П 'К

— í dz í dr[doo) + dio)z + dor) r + d(((r) zr] +

i£o J J

o zn rn

( zn rk

— í dz í dr[co o) + c(o) z + co r) r + с((Г) zr]ln[z2 + r2] +

i£o J J

zn rn

( zn rk

+ -— Idz Idr[ po o) + p(or) z + po () r + p(r) zr]

4neo

22 z + r

(18)

Входящие в последнее выражение интегралы были нами рассмотрены ранее в данной работе или в работе [1].

ЗАКЛЮЧЕНИЕ

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

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

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

дом Гаусса с числом узлов гауссовой квадратуры

** = N. =10.

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

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

r

o

nn

Табл. 1

Номер ячейки гу Относительная ошибка

1 7.462025е-003

20 5.613663е-005

40 4.254564е-005

Табл. 2

Номер ячейки гу Относительная ошибка

1 3.813173е-003

20 4.349788е-005

40 2.596540е-005

Табл. 3

Номер ТН Относительная ошибка

1 5.600178е-004

2 1.207571е-004

3 7.608428е-005

Табл. 4

Отклонение по У Относительная ошибка

0.03 6.692424е-004

0.01 4.834191е-003

0.001 1.350306е-002

Видно (табл. 1), что алгоритм данной работы уступает в точности алгоритму, основанному на формуле (10), который можно принять за идеальный. Причем наиболее заметно это вблизи оси, где функции /1, /2 меняются вдоль оси Я наиболее резко. Поэтому увеличение числа ячеек (более частое деление) вдоль направления Я должно уменьшить погрешность. Это демонстрируется результатами, приведенными в табл. 2, полученными для приведенного выше облака ОЗ, поделенного вдоль направления Я на 80 ячеек.

Вполне очевидно, что при нахождении вклада от всего облака ОЗ, наибольшая относительная ошибка вычисления будет уменьшаться за счет усреднения ошибок от отдельных ячеек. Это демонстрируется результатами табл. 3, где приведены значения вклада в потенциал от всего облака ОЗ для трех ТН: 1 — находящейся посредине самой нижней, 2 — средней и 3 — самой верхней ячеек.

Видно, что хотя результаты алгоритма данной статьи и уступают результатам алгоритма форму-

лы (10), однако могут считаться вполне допустимыми по точности.

Особое внимание следует обратить на интегрирование в случае, когда ТН принадлежит ячейке, по площади которой осуществляется интегрирование. Очевидно, что при близости ТН к одному из узлов численного интегрирования, точность интегрирования может ухудшаться. Для демонстрации этого приведем в табл. 4 значения вклада в потенциал от одной ячейки ОЗ при приближении ТН к узлу гауссовой квадратуры. Величина отклонения по У от одного из узлов гауссовой квадратуры дается в длинах стороны ячейки.

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

— первую стадию решения уравнения Пуассона (нахождение плотности поверхностного за-

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

— расчет компонент электростатического поля на сетке для последующего вычисления траекторий.

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

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

1. Иванов В.Я., Шевченко С.И. О расчете плоских электростатических полей в приборах, имеющих области, заполненные объемным зарядом // Научное приборостроение. 1999. Т. 9, № 4. С.88-94.

2. Ландау Л.Д. и Лифшиц Е.М. Теория поля. М.: Наука, 1967. 460 с.

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

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

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

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

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

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

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

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

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

ON THE CALCULATION OF AXIALLY-SYMMETRIC ELECTROSTATIC FIELDS IN THE SPACE-CHARGE REGIONS

S. I. Shevchenko

Institute for Analytical Instrumentation RAS, Saint-Petersburg

The electrostatic field potential in the space-charge region is found as a solution to the Dirichlet problem for the axially-symmetric Poisson equation. The problem is solved by the collocation method for the respective integral boundary equation. The contribution of the space-charge cloud to the potential at each point of the space is defined as the sum of contributions from the square grid cells computed by taking respective integrals wherein the space-charge density function is interpolated by a bilinear distribution.