Математическая модель переноса вещества в винтовом магнитном поле с использованием граничных условий на бесконечности

Полный текст

1. Введение Для проектировщиков установок для удержания плазмы основным источником информации о динамике вещества является вычислительный эксперимент. Разнообразие и сложность параметров течения плазмы определяют необходимость создания специализированных моделей различных процессов [5]. Исследования течения плазмы в магнитном поле представляют существенный интерес для управляемого термоядерного синтеза [9], изучения стойкости материалов под воздействием мощных тепловых нагрузок [18], лабораторного моделирования астрофизических процессов [16] и ряда других фундаментальных и прикладных научных задач. Наибольшие успехи математики достигнуты в области моделирования систем с тороидальной топологией магнитного поля [1, 4, 15]. Альтернативным подходом является удержание плазмы в открытых магнитных системах, где поле близко к осесимметричному, а его силовые линии пересекают границу области удержания в двух точках [2, 12]. Большой вклад в развитие этого направления внесло математическое моделирование установок, основанных на принципе многопробочного удержания. В качестве развития метода было предложено [13] удержание плазмы магнитным полем с винтовой симметрией. Установка СМОЛА (Спиральная Магнитная Открытая Ловушка) разработана и построена в 2017 году в Институте ядерной физики СО РАН им. Г.И. Будкера для экспериментальной проверки этой идеи [19, 22]. Область, в которой удерживается плазма в установке СМОЛА, с одной стороны ограничена классической пробкой, с другой - многопробочной секцией с винтовым магнитным полем. Установка СМОЛА предназначена для моделирования эффектов винтового удержания при низкой (а потому легко достижимой) температуре плазмы. Для масштабирования винтового удержания на системы термоядерного класса необходимо детальное сравнение экспериментально наблюдаемых потоков вещества с модельными и дальнейший расчёт эффективности системы большего масштаба на основе математической модели. В настоящее время показано соответствие наблюдаемых результатов приближённым теоретическим оценкам. Точного аналитического решения для теории винтового удержания не построено, поэтому востребованы результаты численного решения уравнения переноса плазмы. Рассматриваемая в статье математическая модель переноса вещества в винтовом магнитном поле построена на основе уравнений из работ [14, 19, 22] и параметров установки СМОЛА. Проведена валидация математической модели [17], определено влияние на эффект удержания глубины гофрировки магнитного поля, диффузии и потенциала плазмы. Расчёты показали эффект пинчевания (уменьшения среднего радиуса) плазменной струи, проявляющийся и в эксперименте. Для расширения области параметров, при которых модель будет обладать достаточной предсказательной силой, потребовалось изменения граничных условий. В статье представлены результаты расчётов для условия нулевой концентрации плазмы на большом расстоянии от её источника. Цель работы состоит в предсказании параметров удержания установок на стадии их проектирования. 2. Постановка задачи Рассмотрим движение плазмы в центральной части ловушки, имеющей форму цилиндра. Вещество входит в область удержания из источника плазмы через левый торец цилиндра и выходит в расширитель через правую границу. Рассмотрим поперечное сечение цилиндра в плоскости (r,z). В статье [14] получены выражения для компонентов радиального и продольного переноса частиц в винтовом магнитном поле. Система уравнений описывает динамику плазмы в МГД-приближении в аксиально-симметричной постановке. Различия в движении запертых и пролётных ионов учитываются в виде эффективной силы трения, зависящей от взаимной скорости компонент и доли запертых частиц. Продольная сила, действующая на плазму, возникает в результате взаимодействия радиального электрического тока захваченных ионов с азимутальной компонентой винтового магнитного поля. Учитывается диффузия плазмы поперёк магнитного поля. Исключение зависимых переменных сводит систему уравнений к уравнению неразрывности потока: где u - концентрация вещества, T = Ti +Te, Ti = 4 Эв и Te = 30(1-(r/r0)2) -ионная и электронная температуры, Λ - отношение длины системы к длине свободного пробега λ иона, κ(r,Rm) - доля запертых частиц, l = 216 см - длина системы вдоль силовой линии, Z - среднее зарядовое число одного иона, D - коэффициент диффузии в поперечном поле. Доля запертых частиц κ(r,Rm) = 1 - 1/R(r,Rm), R(r,Rm) = 2(Rm - 1)(r/a)2 + 1 где Rm = 1,52 -глубина гофрировки. Параметр ζ = c/Vz является отношением скорости звука cs = (Te/M)1/2 и продольной скорости Vz движения магнитных возмущений при вращении плазмы в собственном амбиполярном электрическом поле. В уравнении (2.1) физические величины обезразмерены на: r0 = a, z0 = l, φ0 = Te/e, u0 = umax, T0 = Te, где a = 8 см - граница камеры, в которой может существовать плазма. Потенциал электрического поля φ(r) введен в виде полинома, интерполирующего экспериментальные данные [22]: φ(r) = -2,21776 + 1,31r - 7,79r2 + 31,18r3 - 33,43r4 + 10,94r5. В моделируемом эксперименте происходит установление параметров плазмы и электромагнитного поля в течении 40 мс, и затем наступает фаза, когда процесс стационарен в течении 120 мс, после чего разряд отключается. Основная задача экспериментов и математического моделирования состоит в исследовании удержания плазмы, при которых все параметры постоянны. Измеряемый с помощью зондов потенциал плазмы в экспериментах на установке СМОЛА зависит от параметров эксперимента. Максимум φ варьируется в пределах от 2Te/e до 3Te/e. Максимальное значение безразмерного потенциала φ для установок следующего поколения, в которых возможно использование принципа винтового удержания, также лежит в этих же пределах [10, 21]. Экспериментально наблюдаемое распределение потенциала в центральной области плазмы (при значениях радиуса, меньших 0,6 в безразмерных величинах) близко к квадратичному, далее в периферийной области плазмы производная потенциала по радиусу снижается. Погрешность измерения потенциала в эксперименте составляет около 5%. Степень и коэффициенты аппроксимирующего полинома подобраны таким образом, чтобы его отклонение от экспериментально измеренных в референсном эксперименте значений было сопоставимым с экспериментальной погрешностью. Для задания влияния поля в радиальном направлении в модели используется производная абсолютной величины электрического поля: r r2 r3 . Известно, что величина потенциала снижается с увеличением за счёт наличия поперечной проводимости плазмы. Его распределение в пространстве задается следующим образом: . (а) (б) Рис. 1: Зависимости от радиуса установки модуля потенциала электрического поля и его производной (а) и граничного распределения плазмы (б). Fig. 1: Dependences on the radius of the plant: electric field potential module and its derivative (a) and the boundary plasma distribution (b). Рассмотрим область [0,rmax] × [0,zmax] в поперечном сечении центральной части установки. В безразмерных величинах область центральной части установки является единичным квадратом. Рассмотрим расширенную область, имеющую длину, которая в несколько раз превышает длину ловушки. Это позволит использовать условие равенства нулю концентрации плазмы на бесконечности. Теоретические оценки и расчёты показали, что при троекратном и выше превышении длины центральной части установки результаты расчётов практически не изменяются. Считаем, что вещество не достигает стенок ловушки, на оси z ставим условие симметрии, на входе вещества в расчётную область задаем граничное распределение концентрации вещества u(r,0) = uL(r), на выходе u(r,zmax) = 0, zmax = 4. Для расчётов использованы экспериментальные данные [15] для распределения концентрации на границе. На рис. 1(б) приведено экспериментальное значение концентрации плазмы на входе и интерполяция вида: uL(r) = 1,03 + 0,46r - 1,52r2 + 14,48r3 - 44,17r4 + 43,77r5 - 14,05r6. Таким образом, стационарная задача имеет вид: ∂2u(r,z) ∂(Φ(r,z)u(r,z)) ∂ ∂u(r,z) ∂ ∂u(r,z) ⎪⎪⎪ ∂z2 ∂z ∂u L ⎪(0,z) = 0, u(r,0) = u (r), ⎪⎪⎩ ∂r где коэффициенты уравнения: ∂r ∂z ∂u (1,z) = 0, u(r,4) = 0, ∂r ∂r ∂r (2.2) C1(r) = lrT(r)/Λ, C2(r) = lrκ(r)(1 + κ(r))ζ(r), C3(r) = rκ(r), C4(r) = lζ(r)T(r)/Z, C5(r) = κ(r), C6(r) = lrD. ⎧ C1(r)+C2(r)+C3(r) C4(r) +C5(r) C6(r) = 0, 3. Метод решения Задача (2.2) была реализована методом установления [6] и более экономичным методом Зейделяи [11]. Зададим на сетке с узлами ri = ih, i = 1,...,Nr, zk = kh, k = 1,...,Nz сеточные функции, Ci5 = C5(ri), Ci6 = C6(ri), uLi = uL(ri). Критерий сходимости . Для реализации задачи (2.2) методом установления при τ = 10-4 использована безусловно устойчивая схема стабилизирующей поправки и метод прогонки [8]. Новый вариант разностной схемы имеет вид: τ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧⎪⎪⎪⎪⎨ u+-i,k∗ C- ui,kn = C-i1(-ui,k∗ -+1-- 2-ui,k∗ + ui,k∗-1)/hn2 +) +Ci2(Φi,k+1u∗i,k+1 - Φi,k--1u∗i,k--1)/2h + + Cii4-(1C/i24+1(u/i,kn2(+1uin+1,kui,kn ui,knuni-1u,kin+1+1,k-1 n+ ui,k-1 2 + Ci5(Ci6+1/2(uin+1,k ui,kn ) 3 + ui-1,k))/2h Ci6-1/2(uni,k uni-1,k)/h2, (3.1) i = 1,...,Nr, k = 2,...,Nz - 1, u∗1,k = u2∗,k, u∗N r,k = u∗Nr-1,k, ⎪⎪ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪ iu= 2i,kn+1,...,N- ui,k∗ =-C11,(uink+1+1= 1,k -,...,N2uni,k+1 +, uuni-+1n1+1,k)/h= 2uLi-, Ci1u(ni,Nu+1in+1z = 0,k -. 2ui,kn + uni-1,k)/h2, τ i r z i,1 Решение уравнения (2.2) методом Зейделя на каждом шаге по времени позволяет построить более экономичный алгоритм. Разностная схема, основанная на методе Зейделя, имеет вид: ⎪⎪⎨⎪⎪⎪⎧⎪⎪⎪⎪⎪ βuni,k= 2+1 = (Ci1 β-+1Cuini5+1(C,ki6+1+/β2-1 uin-+161,k +) +β+2 Cui,kni3(+1Ci4+1+/β2-2+uni,kC+1i4--11/+2)β/+32,uni+1,k-1 + β-3 uin-1,k+1)/β, + Ci-1/2 β±1 = Ci5Ci6±1/2 + (Ci3Ci4±1/2)/2, β±2 = Ci1 ± (hCi2Φi,k±1 + Ci3Ci4∓1/2)/2, (3.2) ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪ βi ±= 23 =,...,N-(Ci3Ci4-±11/2, k)/2= 2, ,... ,N - 1, un+1 = un2,k+1, r z 1,k unN+1r,k = unN+1r-1,k, uni,+11 = uLi , uni,N+1z = 0. Можно заметить, что задача (3.1) и алгоритм решения (3.2) содержит параметр ζ(r) = 1/Ar, полученный аппроксимацией экспериментальных данных, здесь A = 20. На оси при r = 0 параметр ζ ограничен значением при r = ρB, где ρB = VTi· мс/еВ - ларморовский радиус. Для параметров установки СМОЛА ρB ≈ 0,3-0,4 см. Это обусловлено тем, что ион движется по ларморовской орбите (вращается в магнитном поле), поэтому его радиальная координата осциллирует. В рассматриваемой модели все воздействия осредняются, иону приписывается координата центра окружности, по которой он движется в магнитном поле. То есть, координата равна нулю для ионов, которые облетают ось вокруг и находятся на расстоянии ларморовского радиуса от неё. Поэтому для исключения особенностей решения при расчётах в окрестности оси симметрии безразмерный параметр ζ задается следующим образом: (3.3) . Наряду с другими преимуществами, метод Зейделя интересен удобством использования в цилиндрических координатах. Принцип выражения искомого элемента через соседние по схеме точек вида «крест» универсален и не зависит от выбора системы координат. В ходе расчётов выяснилось, что выбор шаблона (рис. 2) при аппроксимации смешанной производной влияет на решение. Рис. 2: Шаблоны для аппроксимации смешанных производных. Fig. 2: Stencils for approximation of mixed derivatives. Для рассмотрены следующие шаблоны: 1. Девятиточечный шаблон [20]: Λrz = Λzr = C3(r)((C4-(r)uz¯)r + (C4-(r)uz)r¯ + (C4+(r)uz)r + (C4+(r)uz¯)r¯)/4, C4±(r) = C4(r) ± |C4(r)|, 4 Λ+++rzuu=i+1C,ki3+1((-u((Ci-C+1C,k(-| - ||2||C|i4-|1/24|) +| |) +ui-1u,ki,k(| - |--12(|CCii444+1-1//22| - ||||) +) +|Cuui,kii4--4 11(2,k/2|+1C| -|))i4(-C/C14i/4+1ih42+1|2/.+ 2/22| - || - ||CiC4+1Ci4+1i4/+12/|/) +22|) +|) + (3.4) i+1,k-1 i-1/2 + Ci-1/2 ui,k+1 i4-4 1/2 + CCi4i-4-11//22+) +Ci4u+1i-/12,k-1(CCii+14+1//22 + Ci+1/2 i-1/2 Для тестовой задачи: Λrz = C3(r)((C4-(r)uz¯)r + (C4-(r)uz)r¯ + (C4+(r)uz)r + (C4+(r)uz¯)r¯)/4, C4±(r) = C4(r) ± |C4(r)|, (3.5) Λrz = Ci3(-2ui+1,k - 2ui-1,k + 4ui,k - 2ui,k-1 - 2ui,k+1 + 2ui+1,k+1 + 2ui-1,k-1)/4h2. 2. Четырехточечный шаблон [17]: Λrz = C3(r)((C4(r)uz¯)r), Λrz = Ci3(Ci4(ui,k - ui,k-1)r)/h = Ci3(Ci4+1/2(ui+1,k - ui,k) - Ci4-1/2(ui+1,k-1 - ui,k-1))/h2. (3.6) 3. Шеститочечный шаблон [3]: (3.7) 4. Семиточечный левый шаблон [7]: (3.8) 5. Семиточечный правый шаблон [7]: (3.9) В тестовой задаче все коэффициенты равны единице. Результаты применения различных способов взятия смешанной производной можно разделить на три группы: сокращение конечноразностной аппроксимации смешанных производных, возникновение осцилляций, и сходимость к точному решению (см. таб. 1). В случае, когда конечно-разностная аппроксимация смешанных производных для тестовой задачи сокращается, расчёт сводится к решению уравнения Лапласа. Наилучший результат дает использование семиточечных левого для смешанной производной Λrz и правого шаблонов для Λzr. Номер шаблона 1 2 3 4 5 4-5 Сокращаются смешанные производные Осцилляция решения Сходится к точному решению + - - -+ - -+ - + - - + - - Таб. 1: Результат использования шаблонов для аппроксимации смешанной производной для тестовой задачи. Stencil No. 1 2 3 4 5 4-5 Mixed derivatives are reduced Oscillation of solution Converges to the exact solution + - - -+ - -+ - + - - + - - Tab. 1: Results of using stencils to approximate the mixed derivative for the test problem. В случае расчёта задачи о переносе плазмы сокращения симметричной конечно-разностной аппроксимации смешанных производных не происходит, так как в уравнении (2.2) есть только одна смешанная производная и одна первая производная по z. При использовании сеточного оператора для смешанной производной на девятиточечном, четырехточечном и шеститочечном шаблоне возникали осцилляции при коэффициенте диффузии D меньшем, чем 0,001 (см. таб. 2). При использовании семиточечных шаблонов осцилляций не возникает, был выбран левый шаблон. Номер шаблона 1 2 3 4 5 Осцилляция решения Согласуется с экспериментальными данными - - + - + - - + - + Таб. 2: Результат использования шаблонов для аппроксимации смешанной производной для расчёта концентрации плазмы. Stencil No. 1 2 3 4 5 Oscillation of solution Consistent with experimental data - - + - + - - + - + Tab. 2: Results of using stencils to approximate the mixed derivative in calculation of plasma concentration. 4. Результаты расчётов В качестве тестовой задачи рассмотрено уравнение: (4.1) , Задача (4.1) имеет аналитическое решение (рис. 4(а)): . Рис. 3: Решение тестовой задачи (а) и точность сходимости к решению (б) при использовании метода Зейделя (2) и метода установления (1). Fig. 3: Solution of the test problem (a) and the accuracy of convergence to the solution (b) using the Seidel method (2) and the establishment method (1). Методы решения (метод установления и метод Зейделя) потребовали тщательного подбора шаблона для аппроксимации смешанной производной конечными разностями. Получены графики (рис. 4(б)) относительной погрешности для достижения критерия сходимости итерационного процесса εn = 10-8 на последовательности измельчающихся сеток. Метод Зейделя более экономичен и прост в реализации, чем метод установления. Необходимое число итераций метода установления превышает 4000. Метод Зейделя требует до 1000 итераций, причём в случае предлагаемого шаблона для аппроксимации смешанной производной сходимость происходит в два раза быстрее. Распределение концентрации плазмы показано в расчётной области, имеющей в безразмерном виде форму прямоугольника (рис. 4(а)). Для дальнейшего анализа результатов расчётов планируется сравнение с данными экспериментов в сечениях z = 0,4 и z = 1. Расчёты показали (рис. 4(б)), что норма разности решения δ(r) = |u(r,1)|z=10 - u(r,1)|z=k|, k = 2,3,4, меньше 10-4 при размере расчётной области по оси z более, чем z = 4. Для прогноза результатов работы проектируемых установок для удержания плазмы в винтовом магнитном поле необходимо определить более точно коэффициент диффузии. Диффузию можно определить, как ту, при которой расчётные и экспериментальные распределения концентрации плазмы при z = 0,4 и z = 1 совпадают с экспериментальными данными при малых значениях потенциала плазмы. Тестовая задача выбрана с граничными условиями, близкими к условиям в задаче (2.2). На входе в область задается близкое к экспериментальным данным условие (см. рис. 5(а)), в результате получена близкая картина убывания искомой функции вдоль оси (см. рис. 5(а)). Цель работы - это прогноз результатов экспериментов на проектируемых установках для удержания плазмы в винтовом магнитном поле. (а) (б) Рис. 4: Распределение концентрации плазмы в винтовом магнитном поле (а) и зависимость нормы разности решения (б) от размера расчётной области по оси z. Fig. 4: Distribution of plasma concentration in a helical magnetic field (a) and dependence of the norm of the solution difference (b) on the size of the computational domain along the z axis. Рис. 5: Граничное условие (а) и значение искомой функции на оси z (б) для тестовой задачи (пунктир) и задачи (2.2) (сплошная линия). Fig. 5: Boundary condition (a) and the value of the unknown function on the z axis (b) for the test problem (dashed line) and problem (2.2) (solid line). 5. Выводы В работе представлена математическая модель переноса плазмы в спиральном магнитном поле с использованием условия равенства нулю концентрации плазмы на бесконечности. Стационарное уравнение переноса вещества в аксиально-симметричной постановке содержит вторые производные по пространству. Выбран оптимальный шаблон для аппроксимации смешанной производной на основе тестовой задачи с близкими граничными условиями. Используется безусловно устойчивая схема стабилизирующей поправки и более экономичный метод Зейделя, который интересен удобством использования в цилиндрических координатах. Принцип выражения искомого элемента через соседние по схеме точек вида «крест» универсален и не зависит от выбора системы координат. Такая постановка задачи (магнитное поле новой сложной конфигурации) и математическая модель (учёт параметров удержания через задание коэффициентов) используется впервые, возможно в дальнейшем её сравнение с другими моделями (МГД, кинетические модели), которые, возможно, будут развиваться для этой задачи.

Об авторах

Г. Г. Лазарева

Российский университет дружбы народов

Автор, ответственный за переписку.
Email: lazareva-gg@rudn.ru
Москва, Россия

И. П. Оксогоева

Российский университет дружбы народов

Email: oksogi@mail.ru
Москва, Россия

А. В. Судников

Институт ядерной физики имени Г.И. Будкера СО РАН

Email: a.v.sudnikov@inp.nsk.su
Новосибирск, Россия

Список литературы