Диссертация (1149802), страница 5

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 5)

Ng .(1.34)Здесь функция f c представляет собой последовательность следующих операций:1. решение уравнений поля (1.17), (1.18) с использованием сеточных распределений заряда ρg , jg ;2. расчет плотноcти тока эмиссии потоков исходя из условия (1.26);3. решение уравнений движения макрочастиц (1.21);4. расчет распределения пространственного заряда в соответствии с (1.28),(1.29).В том случае, если оператор f c является сжимающим, для решения уравнения (1.34) можно применять метод простой итерацииρnh = f c (ρn−1h ),(1.35)28где n — номер итерации. Функция f c при этом определяется с некоторой погрешностью, т.

е. ее реальное значение будет равноf c (ρh ) = f (ρh )(1 + ϵ(ρh )).(1.36)Здесь ϵ(ρh ) — относительная погрешность вычисления точного значения f (ρh ).При расчете ρnh , согласно (1.35), значение f (ρh )ϵ(ρh ) будет вносить свой вкладв величины ρnh , что может привести к существенным колебаниям решения. Дляподавления вклада погрешности f (ρh )ϵ(ρh ) в решение можно использовать такойподход: уравнение (1.34) будет эквивалентно следующемуρh = f c (ρh ) + (1 − ω) (ρh − f c (ρh )) = (1 − ω)ρh + ωf c (ρh ),(1.37)где ω ∈ (0, 1] — вещественное число. Для решения уравнения (1.37) используетсяметод простой итерацииρnh = (1 − ω)ρn−1+ ωf c (ρn−1hh ).(1.38)Учитывая (1.36), (1.38) перепишется какn−1ρnh = (1 − ω)ρn−1+ ωf (ρn−1hh ) + ϵ(ρh )ωf (ρh ).(1.39)Уменьшая коэффициент ω, можно уменьшить вклад погрешности ϵ(ρh )ωf (ρn−1h )в вычисление ρnh , при этом замедляя сходимость.

Для остановки итерационногопроцесса (1.39) применяется критерий n ρs − ρn−1s < εiter .max s=1...NhρnsЗдесь ρns — плотность заряда на n итерации в s узле расчетной сетки, Nh — числоузлов сетки, εiter — заранее определенная точность.Схема алгоритма итерационного метода представлена на рис. 1.4.29Риcунок 1.4: Схема итерационного метода30Глава 2Численные алгоритмы в методах частица-сеткаРассмотрим в данной главе детальное описание численных методов, которыебудут применяться в дальнейшем при решении задач с помощью итерационногометода и метода частиц в ячейках.2.1 Расчет траекторий частицИнтегрирование траекторий макрочастиц как правило является наиболее трудоемким этапом в методах частица-сетка, поскольку число макрочастиц можетдостигать порядков > 105 , при этом каждую траекторию необходимо рассчитывать отдельно от остальных. Среди множества явных методов интегрированиятраекторий частиц [33, 42, 74], стандартом де-фактов в программах для моделирования динамики пучков и плазмы стали метод с перешагиванием в случаеэлектрических сил [75] и схема Бориса [76–79] с случае наличия магнитногополя.

Перечислим основные причины популярности этих методов [78, 79]:• Второй порядок точности при необходимости только одного вычислениясилы, действующей на частицы на каждом шаге. При этом, другие методы,такие как метод Рунге-Кутты второго порядка требуют двух вычисленийсилы на каждом временном шаге.• Метод с перешагиванием и схема Бориса обратимы по времени.

Ошибкиаппроксимации законов сохранения полной энергии и момента импульса,возникающие при расчете данными методами, ограничены. Данные методыдостаточно устойчивы при расчете длинных траекторий.31Схема Бориса для интегрирования релятивистских уравнений (1.21) в векторномвиде в потоке с номером α имеет вид()k+1/2k−1/2p+ppk+1/2 − pk−1/2qα=Ek +× Bk ,∆τm0α c22γ k(2.2)rk+1 − rkpk+1/2= k+1/2 .∆τγЗдесь k — номер шага интегрирования, ∆τ — шаг интегрирования по времени,Ek = E(rk ), Bk = B(rk ).

Последовательность вычисления новых координат иимпульсов по этой схеме схематично приведена на рис. 2.2.Риcунок 2.2: Последовательность вычисления новых координат и импульсов частиц2.1.1 Случай декартовых координатВ случае двумерной задачи в декартовых координатах без магнитного поля,схема (2.2) для интегрирования уравненений движения макрочастиц в потоке сномером α (1.21) примет видqαExk ,2m0α cqαEyk ,= p k−1/2+ ∆τy2m0α cp xk+1/2 = p xk−1/2 + ∆τp yk+1/2k+1/2k+1xpx= x + ∆τ k+1/2 ,γkk+1/2yk+1py= y + ∆τ k+1/2 .γk(2.3)32Здесь px = vx γ/c и py = vy γ/c — компоненты импульса частицы, x и y — еёдекартовые координаты.2.1.2 Случай цилиндрических координатДля расчета траекторий в цилиндрических координатах (r, φ, z) будем использовать приведенные импульсы частиц pr = vr γ/c, pφ = rvφ γ/c, pz = vz γ/c.Тогда схема (2.2) в цилиндрических координатах с азимутальным магнитным полем для интегрирования уравнений движения макрочастиц в потоке с номеромα (1.21) примет вид [79][p k+1/2rqα= p rk−1/2 + ∆τm0α c2(k+1/2Erk −k+1/2+(p φc(p zk−1/2 2+ pφ4(rk )3 γ k)k−1/2+ pz2γ k]))Bφk +,qαEφk rk ,2m0α c[]k+1/2k−1/2qαc(p r+ pr) kk= p zk−1/2 + ∆τE+Bφ ,zm0α c22γ kp k+1/2= p φk−1/2 + ∆τφp k+1/2zk+1/2rpr= r + ∆τ k+1/2 ,γk+1kk+1/2k+1φ4p φ= φ + ∆τ k+1/2 k+1,γ(r+ rk )2kk+1/2zk+1pz= z + ∆τ k+1/2 .γkПри вычислении импульсов γ k определяется с помощью подхода Борисаkγ =√√1+(p k )2=1 + (p kr )2 + (p kz )2 + (p kφ /rk )2 .(2.4)33Для того, чтобы вычислить p k , необходимо сделать половинный шаг с отключенным магнитным полем∆τ qαEk ,22 m0α c()k−1/2 2∆τqα(p φ)Erk + k 3 k−1/2 ,p kr = p rk−1/2 +22m0α c(r ) γ(p k ) = pk−1/2 +∆τ qαEφk rk ,22 m0α c∆τ qαp kz = p zk−1/2 +Ezk .22 m0α cp kφ = p φk−1/2 +Использую данные выражения, из схемы (2.4) можно явным образом выразитьk+1/2значения p rk+1/2, pzk+1/2, pφ.Задачи расчета траекторий в двумерной задаче в полярных координатах безмагнитного поля и в двумерной осесимметричной задаче с азимутальным магнитным полем являются частными случаями задачи расчета траекторий в рассмотренной трехмерной задаче.

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

Вне зависимости от используемой системы координат (q1 , . . . , qNq ) на каждом шаге интерирования kдля каждого потока c номером α для каждой частицы с номером i выбираетсяиндивидуальный шаг ∆ταki из условия устойчивости Куранта — Фридрихса —Леви34∆ταki = min hkqs α /(Kβqks α )).s=1...Nqi(2.5)iЗдесь hkqs α — размер ячейки эйлеровой сетки по координате qs , в которой нахоiдится частица с номером i из потока c номером α на k-м шаге интегрирования;βqks α — компонента qs приведенной скорости данной частицы; K > 1 — задаваеiмый коэффициент.Для метода частиц в ячейках шаг интегрирования основной процедуры расчета ∆τ для всех потоков частиц выбирается исходя из соотношения∆τ = min min min ∆ταki .αik(2.6)Для итерационного метода для каждого потока α выбирается индивидуальныйшаг интегрирования основной процедуры расчета ∆τα исходя из соотношения∆τα = min min ∆ταki .ik(2.7)Выбор шага согласно предложенному алгоритму обеспечивает прохождение всеми частицами не более чем K-й части ячейки эйлеровой сетки.

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

В таких областях со сгущенной сеткой необходимо использовать значительно меньший шаг интегрирования, чем в остальной области. Для35сокращения объема вычислений в данных задачах предлагается алгоритм ступенчатого изменения шага.Как и в случае выбора постоянного шага, перед основным этапом расчетазапускается тестовая процедура и вычисляются ∆τα исходя из (2.7). На каждомшаге основной процедуры расчета вычисляется ∆ταki исходя из (2.5). Вычисляется соотношениеRαk∆ταki= min.i∆ταШаг интегрирования выбирается исходя из соотношения∆ταk = rj ∆τα ,rj : Rαk ∈ (rj , rj+1 ).(2.8)Здесь r — вектор, содержащий коэффициенты увеличения шага. Например выбор r = (1, 10, 100, 1000, .

. .) обеспечивает ступенчатое десятикратное изменениешага. Предлагаемый алгоритм является компромиссом между необходимостьюсократить объем вычислений и сохранить точность расчета траекторий макрочастиц.2.2 Форма макрочастиц, расчет пространственного заряда исил2.2.1 Случай метода частиц в ячейкахВыбор функции формы макрочастиц Sr (r, r(t)) в значительной степени характеризует методы частица-сетка. Например, широко известными формамимакрочастиц являются: точечная частица (NGP), облако в ячейке (CIC), облакотреугольной формы (TSC).

Рассмотрим далее более подробно форму CIC, которая в дальнейшем будет использоваться при всех расчетах. В одномерном случае36в декартовых координатах функция Sr (x, x) выбирается следующим образом1/hx , |x − x| ≤ hx /2;Sr (x, x) =0, |x − x| > hx /2.Здесь hx — размер макрочастицы. Положим его равным шагу эйлеровой сеткиhx = xg+1 − xg . Тогда, используя (1.27), получим весовую функцию SxSx (xg , x) =(hx − |xg − x|)/hx ,0,|xg − x| ≤ hx ;(2.9)|xg − x| > hx .Чтобы свойства гладкости одномерных весовых функций сохранились, для двумерных или трехмерных случаев, весовые функции нужно выбирать в виде произведения соответствующих одномерных функций, напримерS(xg , yg , zg , x, y, z) = Sx (xg , x)Sy (yg , y)Sz (zg , z).(2.10)Используя (2.10), запишем расчетную форму выражений (1.28), (1.31) для трехмерных декартовых координат:ρg (t) =Eαi =1 ∑∑qαi Sx (xg , xαi (t))Sy (yg , yαi (t))Sz (zg , zαi (t)).hx hy hz α i∑Eg (xg , yg , zg , t)Sx (xg , xαi (t))Sy (yg , yαi (t))Sz (zg , zαi (t)).(2.11)(2.12)Аналогичным образом запишутся эти выражения в случае двумерных осесимметричных координат:∑∑ρg (t) =Eαi =∑αqαi Sr (rg , rαi (t))Sz (zg , zαi (t))iπ((rg + hr )2 − (rg − hr )2 )hz.Eg (rg , zg , t)Sr (rg , rαi (t))Sz (zg , zαi (t)).(2.13)(2.14)37В методе частиц в ячейках выражения (2.11)-(2.14) необходимо применять накаждом шаге дискретизации задачи, обрабатывая все макрочастицы и узлы эйлеровой сетки.2.2.2 Случай итерационного методаВ случае применения итерационного метода предполагается, что у задачи существует стационарное решение, и соответственно, плотности тока частиц j(r)и плотности тока отдельных потоков jα (r) не зависят от времени и для них выполняются уравнения неразрывностиdivj(r) = 0,divjα (r) = 0,α = 1 .

Характеристики

Список файлов диссертации