Диссертация (1149802), страница 6
Текст из файла (страница 6)
. . Nα .(2.15)Представим jα (r) в виде суммыjα (r) =∑jαi (r).(2.16)iЗдесь jαi (r) — вектор плотности тока, вычисленный вдоль траектории макрочастицы с номером αi .Рассмотрим производную по направлению1 ∂|jαi (r)| grad|jαi (r)|vαi (r) ∂jαix ∂jαiy ∂jαiz==++= divjαi (r). (2.17)|vαi (r)| ∂vαi (r)|vαi (r)|∂x∂y∂zИз (2.15)-(2.17) следует, что производная по направлению траектории любоймакрочастицы от модуля вектора плотности тока вдоль этой траектории равнанулю. Соответственно, плотность тока вдоль траектории любой макрочастицысохраняется. Используем далее модель "токовых трубок".
Предполагается, чтофизические частицы, входящие в расчетную область через некоторую площадку∆s образуют некоторую токовую трубку постоянного поперечного сечения ∆s.При необходимости также проводится дискретизация по начальным скоростямна интервалы ∆v. Тогда, исходя из распределения плотности тока эмиссии, не38зависящего от времени, траектории каждой моделируемой макрочастицы, эмиттируемой в центре интервалов ∆s и ∆v можно поставить в соответствие некоторый ток Iαi :Iαi = jemα (r, v).(2.18)Ток Iαi , также будет сохраняться вдоль каждой траектории, как и модуль плотности тока.Для того, чтобы вычислить заряд, вносимый траекторией каждой макрочастицы в узлы эйлеровой сетки, рассмотрим участок траектории между точкойk+1/2r αik−1/2и r αik+1/2, которые соответствуют временным моментам tαk−1/2и tα, по-лученными в процессе интегрирования уравнений движения макрочастицы спомощью схемы с перешагиванием, или Бориса [44].
В соответствии с предположением о непрерывности и стационарности потока частиц, данный участок траектории можно рассмотреть как участок, непрерывно заполненный Np макрочаn−1/2n+1/2стицами с зарядами qp,αi = Iαi ∆tp расположенными в точках rp ∈ (r αi , r αi )∑n+1/2n−1/2∆tp = tα− tα, Iαi ток, соответству(см. рис.
2.3, здесь p = 1 . . . Np ,pющий рассматриваемой траектории). В соответствии с (1.28) заряд, вносимыйкаждой этой макрочастицей в узел rg определяется следующим выражениемgqp,α= qp,αi S(rg , rp ) = Iαi ∆tp S(rg , rp ).iТаким образом, общий заряд, вносимый в узел rg участком траекторииk−1/2(r αik+1/2, r αi) будет определяться какgqk,αi=Np∑p=0gqp,αi=Np∑∆Iαi ∆tp S(rg , rp ).p=0Данное выражение представляет собой сумму Римана. Используя предельныйпереход ∆tp → 0, мы можем записать его в виде интеграла39Риcунок 2.3: Схема метода раздачи пространственного заряда.t∫αk+1/2gqk,α= ∆IαiiS(rg , r(t))dt.k−1/2tαСледуя концепции метода с перешагиванием, мы рассматриваем движениеk−1/2частицы на каждом интервале (r αik+1/2, r αi) как движение с постоянной ско-ростью. Таким образом, t = ξ/v αi , где ξ — расстояние от точки r k−1/2 до точкиr(t).
И (2.2.2) может быть записано какk+1/2gqk,α= ∆Iαiitαk−1/2− tαl∫lS(rg , r(ξ))dξ,0где l — длина рассматриваемого участка траектории. На практике данный интеграл может будем вычислять с использованием квадратурной формулы][k+1/2k−1/2S(r,r)+S(r,r)αiαigggqk,α≈ ∆Iαi (tαk+1/2 − tk−1/2.)αi2(2.19)Данное выражение имеет дополнительное преимущество. оба значения весовойфункции на данной стадии уже должны быть вычислены на стадии интерполя-40ции поляE(r αk−1/2)i=∑W (rg , r k−1/2)Eg .αiТакже данное выражение является достаточно общим, таким образом, что использование различных функци S в (2.19) позволяет простым образом построитьразличные схемы раздачи заряда (CIC, TSC, PQS) [33]. В случае, когда участоктраектории пересекает границы ячейки сетки, мы должны разбить этот участокна подынтервалы, каждый из которых полностью лежит внутри некоторой ячейки и вычислять выражение (2.19) и время прохождения частицей для каждого изних.
Также в данном случае необходимо производить дополнительные вычисления функции S(rg , r) для точек пересечения участка траектории и границ ячеексетки. Общий заряд, вносимый в узел сетки rg будет определяться следующимвыражениемqg =∑∑∑αkgqk,α.ii2.3 Решение уравнения ПуассонаРассмотрим задачу решения уравнения Пуассона в осесимметричных координатах (r, z) (в случае двумерных и трехмерных декартовых координат решениепроводится аналогичным образом). В этом случае уравнение имеет вид()∂U∂ 2Uρr+=− ,2∂r∂zε0∂U∂UEr = −, Ez = −.∂r∂z1 ∂r ∂r(2.20)(2.21)Для решения используется метод конечных разностей на пятиточечном шаблоне "крест". Для более точной аппроксимаци граничных условий, эйлерова сетка строится таким образом, что ее ячейки полностью покрывают41расчетную область, а дляU’upприграничных узлов шаблонUupразностной схемы модифицируется с учетом действитель-h upU lefth righth leftUright U’rightU0h downU downных расстояний от узла до границы [81].
На рис. 2.4 приведен пример шаблона разностной схемы, когда границарасчетной области пересекаетячейки эйлеровой сетки.В общем случае, конечно-Риcунок 2.4: Шаблон разностной схемы.Сплошными линиями обозначены ячейкиэйлеровой сетки, жирными сплошными —границы расчетной областиразностное уравнение для любого внутреннего узла сеткиможно представить в следующем виде2Ulef t2Uright2Uup+++hlef t (hlef t + hright ) hright (hlef t + hright ) hup (hup + hdown )Uright − U02U02U0ρ2Udown+−−= .+hdown (hup + hdown )r × hrighthup hdown hup hdownϵ0(2.22)Для решения системы, составленной из уравнений вида (2.22) реализованы метод верхней релаксации и стабилизированный метод бисопряженных градиентов.При расчете динамики потока частиц с вблизи криволинейной поверхности, необходимо точно вычислять напряженность электрического поля вблизинеё [38, 53].
Для того, чтобы избежать роста погрешности, возникающей причисленном дифференцировании потенциала, когда шаги hlef t , hright , hup , hdownдостаточно малы, применяется экстраполяция потенциала за границу расчетнойобласти. Например, в случае, приведенном на рис. 2.4 вычисляются значения′′Urightи Uupс помощью линейной экстраполяции42′Uright= U0 +hhright′Uup= U0 +(Uright − U0 ),h(Uup − U0 ).hupПосле чего напряженности электрического поля вычисляются в точках, находящихся за границей, с использованием левой или правой разностной производнойпо трем значениям потенциала в узлах эйлеровой сетки.2.4 Методы расчета тока, ограниченного пространственнымзарядом2.4.1 Модель эмиссии Чайлда-ЛенгмюраАналитические решения для тока, ограниченного пространственным зарядом, были получены Чайлдом [82] для одномерного диода и Ленгмюром с Блоджетт для одномерного цилиндрического [83] и сферического диодов [84].1.
Плоский случай4ε0J=9√2q V 3/2.m 0 d2(2.23)2q V 3/2.m0 ra rc β 2 (γ)(2.24)2. Цилиндрический случай4ε0J=9√3. Сферический случай4ε0J=9√2q V 3/2.m0 α2 (γ)rc2(2.25)Здесь V — напряжение между электродами диода, ε0 - диэлектрическая постоянная, rc — радиус катода, ra — радиус анода, d — расстояние между катодом ианодом, γ = log(rc /ra ) функции β и α вычисляются следующим образом:β(γ) = γ − 0.4γ 2 + 0.091667γ 3 − 0.014242γ 4 + 0.001679γ 5 ,43α(γ) = γ − 0.3γ 2 + 0.075γ 3 − 0.014318γ 4 + 0.002161γ 5 .Наиболее широко используемым способом описания эмиссии в двумерных илитрехмерных случаях с помощью итерационного метода является использованиеэтих аналитических решений [42,43,45,47,50]. Следуя данному подходу, необходимо выбрать параметр d (в цилиндрическом (2.24) и сферическом (2.25) случаях параметр ra ) и множество точек, принадлежащих эмиссионной поверхности.Таким образом, вблизи эмиттера будет выделено множество виртуальных диодов.
На каждой итерации метода применяется одно из выражений (2.23)-(2.25)для каждого виртуального диода с целью пересчета плотности тока после каждого пересчета электрического поля. При этом геометрическая форма границыэмиссии определяет выбор используемого закона. Так, в плоской геометрии очевидным будет выбор закона (2.23), а в случае цилиндрического диода — закона(2.24).
Однако в случае сложной геометрии эмиттера ни один из приведенныхтрех аналитических законов может не дать приемлимый результат (например,случай эллиптического эмиттера, который будет приведен ниже). Также необходимо отметить, что выражения (2.23)-(2.25) необходимо модифицировать принеобходимости учета начальной энергии частиц и при наличии начального разброса пучка по скоростям [83].Метод виртуального катода также применяется для описания эмиссии в методе частиц в ячейках [40, 85]. В таком случае заряд новых частиц, входящихв расчетную область на каждом временном шаге дискретизации, вычисляется спомощью одного из выражений (2.23)-(2.25).2.4.2 Оптимизационный алгоритм в итерационном методеВ некоторых работах, например [42, 46] предлагается подход, который заключается в прямом нахождении функции плотности тока эмиссии из условияравенства нулю электрического поля на катоде.
В предложенном в работе [42]подходе предлагается аппроксимировать плотность тока эмиссии с помощью44кусочно-постоянной функции, заданной с помощью набора значений в точкахстарта траекторий с последующим применением многомерного метода Ньютона. При этом задача определения тока, ограниченного пространственным зарядом, сводится к набору задач с известным током эмиссии. Такой метод требуетзначительного объема вычислений — в общем случае для расчета производныхнеобходимо решить число задач с известным током эмиссии, равное числу точек,задающих функцию распределения плотности тока.Однако объем вычислений, требующийся в данном методе, можно существенно сократить с помощью следующего подхода. Для двумерных и осесимметричных задач в случае эмиттера, представляющего собой некоторую кривую,представим плотность тока в виде полинома [53, 56]J(l) =Np∑cp l p .(2.26)p=0В (2.26) l — длина кривой эмиттера от некоторой начальной точки до точкина эмиттере, в которой определяется плотность тока, c1 , .
. . , cNp — параметры,подлежащие определению. В таком случае напряженность электрического поляна эмиттере будет являться функцией от параметров c1 , . . . , cNp и условие (1.26)будет эквивалентно задаче минимизации вектор-функцииE 2 (Γ0 ) av e F (c1 , . . . , cNp ) = . . . → 0.Np2Eav (Γe )(2.27)Тогда для нахождения параметров c1 , . . . , cNp , определяющих распределение токаэмиссии, будем использовать модифицированный метод Ньютона [86]k−1)−1(ckc11 1 1 ∂F(c,...,c)1Np k−1F (ck−1. . . = .















