CC BY
42
2009
НАУЧНЫЙ ВЕСТНИК МГТУ ГА серия Математика и физика
№ 140
УДК 533.9:51-7
МОДЕЛИРОВАНИЕ ДИНАМИКИ СТОЛКНОВИТЕЛЬНОЙ ПЛАЗМЫ ВБЛИЗИ ЦИЛИНДРИЧЕСКОГО ЗОНДА
И.А. КУДРЯВЦЕВА, А.В. ПАНТЕЛЕЕВ
В работе предложена модель, описывающая динамику столкновительной двухкомпонентной плазмы вблизи бесконечного заряженного до некоторого потенциала цилиндра [1]. Предполагается, что плазма состоит из электронов и однозарядных ионов, поэтому учитывались лишь кулоновские взаимодействия. Модель включает уравнение Фоккера-Планка, описывающее данный тип взаимодействия, и уравнение Пуассона. Предложен метод решения, базирующийся на методе Монте-Карло. Приведены и проанализированы результаты численного моделирования.
Ключевые слова: столкновительная плазма, динамика, уравнение Фоккера - Планка, численное моделирование.
Постановка задачи
Рассматривается следующая физическая постановка зондовой задачи [1,2]. В невозмущенную бесконечно-протяженную плазму, состоящую из электронов и однозарядных ионов, внесен бесконечный цилиндр радиусом гр, заряженный до потенциала рр. Поверхность цилиндра является идеально поглощающей для электронов, ионы при ударе о поверхность нейтрализуются. Концентрации ионов п¥ и электронов пе¥, а также температуры данных частиц Т¥ Т¥ в невозмущенной плазме заданы. Каждая частица движется в макроскопическом электрическом поле, магнитное поле отсутствует. Начальные распределения для всех типов частиц предполагаются максвелловскими.
Требуется с учетом столкновений между заряженными частицами найти самосогласованное электрическое поле Е(г, *), плотности токов на зонд ионов (г,*) = ц■ | /(г,V,*)у ёг и электронов ]е (г,*) = е ■ | /е(г,V,*)у ёу, где ц = 2{е , = 1 - заряд иона, е - заряд электрона, а также концентрации вблизи зонда ионов п{ (г,*) = | /(г,V,*)ёу и электронов
пе (г,*) = | /е(г,V,*)ёу , где /е (г,V,*), / (г,V,*) - функции распределения электронов и однозарядных ионов соответственно. Поведение частиц во времени * характеризуется радиус-вектором г и вектором скорости V .
Математическая модель, соответствующая данной физической постановке, имеет вид:
Э/а(г,?,*) . V Э/„(г,V,*) + Ра(г,*) Э/Лг,V,*)= (д/а(г,V,*)^
- + V-------------+
а
Э* Эг та ЭV
е
Э*
+ ^( Г, V, * ),
(1)
Ар(г,*) =-------------------(п (г, *) -пе (г,*)), Е(г,*) = -Vр(г,*),
ео
где первое из уравнений системы (1) - уравнение Фоккера-Планка для частиц сорта а (а = /, е); второе - уравнение Пуассона для самосогласованного электрического поля;
Г Э/а( ^ ^ 1
----------- - оператор столкновения
I Э )с
Фоккера-Планка; функция 8а( г, V, *) описывает источники или стоки частиц;
/а(г,V,*) - функция распределения частиц сорта а;
с
- - \ЧЕ(-,t), а = г, - -
Еа(г, I ) = < - , где Е (г, I) - напряженность самосогласованного электрического
еЕ (г, t), а = е,
поля; ср(г,t) - потенциал самосогласованного электрического поля; па(г,t), а = г,е - концентрация частиц сорта а; та - масса частицы сорта а; £0 - электрическая постоянная.
Оператор столкновения Фоккера-Планка имеет вид [3,4]:
Т1- 001 = 2 : (/а ^,&,(t ))-У. • (/а ^).
где VVVVga (г, V, t) - ковариантная тензорная производная второго ранга, символ двоеточия (:) обозначает операцию двойного суммирования:
Г
In D„
D _ ПжЄокТ^
¿У
f Є0 kT л 1/2
0_____Є¥
V ПЄ¥Є J
Za_ 1, a_i,e,
Z і
ga (^ V t )_ Z If f ( P, ^ ’ 1) V - V 1 ^ ', a_ ^ e,
b_i,e V Za J
¿a(F,v 't)_ Z f f^lf f,( ^ V', t) dv', a_ i, e .
да,
V Z a ,
V - V
b=i,e "*b
К системе уравнений (1) необходимо добавить начальные и краевые условия:
t = 0: fa(r,F.0) = /Г”. a = i,e,
r eOp: f„(r.v.t)| = 0. a = i.e.
p Ja\ ’ 3 ' IreO p ’ ’ ’
Ft )l NO, = П. r eO. : /„ (F. v. t )| = /„—. a = i. e.
j( Л t)
?єй_
_ 0,
(2)
maksv
■n„
ma
V 2kpTa¥ J
exp
ma
V 2kTa¥
V - V.
, v. _ 0, a_ i, e; W p, W¥ - множество ра-
J
диус-векторов частиц, концы которых принадлежат поверхности зонда и границе возмущенной зоны соответственно.
Метод решения задачи
Для решения поставленной задачи перейдем к цилиндрической системе координат, а цилиндр разместим таким образом, чтобы его ось симметрии совпадала с осью 0і. Тогда положение частицы в пространстве будет определяться координатами г, в, 2, где г = г (і) - величина радиус-вектора, в = в(і) - полярный угол, а скорость определится координатами уг, ув, у2 . Функции распределения частиц и потенциал самосогласованного электрического поля инвариантны относительно сдвига по 2 | = 0, = 0 |. В результате в цилиндрической системе
^ ді ді )
координат система (1) - (2) будет иметь вид:
2
V З/а + -вЛ + Зі г дг г дв
( 2 Е ^
+д«Ег
г т„
л
(Ь’,
+
ЧаЕв -г-в
т„
^ГД/«, а = і, е,
д-в
2
д2Л д^а + д'/а дУа + 2_ д2/а дV«
д [-г ]2 д [-г ]2 д [-в]2 д [-в]2 З- д-вд-г д-в
Эу
г V
Эу
З ( / Э*а^
./с
г У
Эу,
V д-в У
д2р + 1 др + 1 д2р _ е
дг2 г дг г2 дв2
_-------------( Пг - Пе ) ,
Е _-Зр Е _-1 Зр
г дг ’ в г дв’
г
і _ 0: /а (г, в, V, , Ув,0)_ /
так&-
г _ г^ : Ла( г^,в -г , в І )_ 0
Р( г,, І )_Р,,
а_ г, е, а _ і, е,
(3)
таЬ-
а_ г, е,
г _ г¥ : /а ( г¥,в, -г , в і )_ Л -
р(г¥, і) _ 0.
Перейдем к безразмерным величинам, воспользовавшись системой масштабов [5]. Тогда задача (3) в безразмерном виде будет выглядеть следующим образом:
Л
Зі
+
К
( Л ( Л
- д/а+ -в Л + -в + ^ аЕг д/а + гаЕв -г-. д/а
г дг г дв г 2Єа З-г 2Єа г д-в
V У V У __
К/а = 2 АС
( Л
З2 ? д ^
¿с
+-
д[V]2 V д[V]21 д[,-в]21 -д[V.]
-да
д'?.
л
д/г„
V д-г У
_д_
д-в
( ? З/а ^
J а -ч ^
З-в У
■>?а
З2
+ 2-
д
_ ГаК/а, а_ г, е,
2 ( 32л Л
д-г д-в
Ла
V д-г д-в у
(4)
2
д
2
д
д р 1 др 1 д р ,Л Л ч ? др ? 1 др
^ + г?эТ+?в=-("'-п}, Е =-э7’ е?._-г?эв-
і _ 0: /а (г, в, -г, -в,0)_ , а_ і, е,
г _ і : Ла ( і,в, - , ^ ? )_ 0, а_ К e,
Ф ( ? , ?)_ р?р,
г _ г- : /а (і ,в, -г, -в, і)_ /ґ", а_ і, е,
Ф (г-, ?)_ 0,
Л
так&і
-ехр (-'; - % )
а = г,
р
1
р
ехр
-т^( V;-+V;)А
, т Т.
\ е га
а = е,
у
где 8а=^, та
Т
са
Т
Ма
т Мам,
."'а да_ г , " , А<? / \4
т г /- а4
(ма)
Аа = мам,
(ма)
ма, маа - коэффициен-
ты, определяемые аналогично задаче диагностики столкновительнои плазмы плоским зондом
[5].
Для решения поставленной задачи применим метод статистических испытании Монте-Карло [5]. Для этого от системы (4), содержащей уравнение Фоккера-Планка, перейдем к системе уравнений, включающей стохастическое дифференциальное уравнение (СДУ) Ито:
^а(ї) = аа(^ ^а(,Л)) + Sа(ї, ^а(,)¥Ж(,) С = К е
а \ / а \ ’ а'
Э;р 1 Эр 1 Э2 р „ ч
3^+7 37+? ;=-( "г- "е) •
Е = -ЭР Е =-1ЭР 7 эя’ ; , э;’
(5)
7 = 0: [V,(?) ',(()] ~/Г*”’, а = г,е,
, = гр : р (гр,ї ) = р .г = ,а : р (,а, ї ) = 0,
где
^ = [,(?) ;(/) V, (?) V; (?)]Т , аа(7, ^а(Г)):
А
■■ас
Эа
Э'г
э ; #а
-+
да
а:
Эа
-+
Да
*«е, 'Л
, С, «О; (?,^„(?)) =
Э' Э^Р„
5,q = ,,;, а = г,е, Ж(?)- стандартный винеровский случайный процесс.
Решение системы уравнений (5) начинается с решения граничной задачи для уравнения Пуассона:
Э ;р 1 Эр 1 Э ;р ,
Э2+7 *+? ;=*( 7,;),
р| ,=ї =рр, р ?=,_ =0,
(6)
(7)
где 7/(г,в) = -(ц -ц) является периодической четной функцией по в: г](г,в)=л(г,в + 2р), г}(г,в) = ц(г, -в). По найденным значениям потенциала р вычисляются значения элементов
Е, ^ = г, в вектора напряженности электрического поля с использованием аппроксимации второго порядка точности на трехточечном шаблоне [6].
Уравнение Пуассона, записанное в цилиндрической системе координат, предлагается решать методом разделения переменных Фурье. Этот метод позволяет перейти от граничной задачи для уравнения Пуассона к краевой задаче для обыкновенных дифференциальных уравнений.
Т
Согласно методу разделения переменных Фурье решение задачи (6) - (7) в силу четности функции )(r,q) ищется в виде разложения в ряд по системе косинусов:
j( r,q) =¿b (r)cos (qq) •
q=0
Для перехода к системе обыкновенных дифференциальных уравнений необходимо разложить в ряд по косинусам также функцию ) (r, в):
1 г„(„ 2cos (чвУ
h( r,q)=—|h( r,q) dq+Z—■—jh( r, q) cos (nq) dq.
p 0 9=1 p 0
В результате приходим к следующей системе уравнений:
і 1 Р
К + тД'= — fh(r,q)dq,
Г Р о
1 п ^ 2 Р
К+7К -^Д=p Иrq)cos(nq)dq,
Д (rp ) = j, Д (r~) = °
Д (rp) = 0 Д (r~) = 0, п
Р,
Решение данной системы ищем путем перехода к ее конечно-разностному аналогу и применения метода прогонки [6].
Далее ищется решение СДУ Ито помощью стохастического метода Эйлера [7] согласно следующей явной разностной схеме:
^ =П + Vаа&,^) + оа(1к,^)АЖк, к = 0,...,К,, а = /,е, (8)
где - количество шагов по времени, АЖк = Ж(?+1) - Ж(?) , - величина приращения винеров-
ского процесса Ж (?) на отрезке [?к, ?к+1 ], по определению независимая от V?, АЖ0,..., АЖк-1:
АЖк ~ N(0,Ит) , кТ - шаг интегрирования по времени, то есть АЖк представляют собой гауссовские случайные величины с нулевыми математическими ожиданиями и дисперсиями, рав-
Г" ~\Т ^ ^
ными шагу интегрирования; У0 = Уг (?0), Ув (?0) , У0 ~ £аЛь'у
эа
Для решения СДУ Ито (8) требуется вычислить значения элементов ——, 5 = г,в вектора
к Э2 І? л к Т ~ к
аа(?к,^а) и элементов -------а, 5,q = г,в матрицы (7а(,к,¥а)са(,к,¥а). Частные производные
Э' Э?
—а, 5 = г, в находятся с использованием аппроксимационных формул численного дифференцирования второго порядка точности для трехточечного шаблона по найденным значениям
? Э ;г
функции Иа. Частные производные ----------а, 5, q = г,в ищутся с помощью аппроксимационных
Э' Э')q
формул первого порядка точности для аналогичного трехточечного шаблона [6] по найденным значениям ?а. Для вычисления интегралов, входящих в выражения для функций Ъа и ?а, применяется метод Монте-Карло [6]:
ц =м (V [¿.г') ), ([^=
V ~ 1,ї = [V, V. (9)
Из решения СДУ Ито находятся значения *¥а вектора состояния *¥а(ї) в дискретные моменты времени ік = і0 + ккг, к = 0,..., Ы(, на основании которых вычисляются методом Монте-Карло интегралы (9), а также концентрации частиц пі, пе и токи частиц на зонд ]і, ]е.
Результаты численного моделирования
В данной работе согласно описанному выше методу решения задачи составлен алгоритм и разработана программа, позволяющая по известным в начальный момент времени потенциалу зонда (рр, радиусу зонда гр, а также концентрациям частиц п¥, пе¥ в невозмущенной плазме,
гг! гг! гтакяу гтак^у ^
их температурам 1і¥ 1е¥ и распределениям )і , )е получить зависимости плотностей токов
частиц на зонд ]і, ]е, концентраций частиц пі, пе и компонент вектора напряженности самосогласованного электрического поля Ег, Ев от времени и пространственных координат.
В начале проведения вычислительного эксперимента невозмущенные концентрации задавались равными пі¥ = пе¥ = 1012 -*101б[-^у]. Для данных выбранных начальных значений концен-°° м
траций правая часть уравнения Фоккера-Планка близка к нулю и имеет место бесстолкнови-тельный режим обтекания цилиндрического зонда плазмой. Данный факт позволил сравнить полученные результаты с результатами, представленными в [2]. При заданных невозмущенных концентрациях исследовалась зависимость выходных данных от изменения значений потенциала, подаваемого на зонд, и радиуса зонда. Это позволило выявить совпадение результатов решения данной задачи с результатами, представленными в [2].
Для второй серии вычислительных экспериментов невозмущенные концентрации варьировались в диапазоне п¥ = пе¥ = 1020 ^1022[^т]. В данном случае правая часть уравнения Фоккера-
м
Планка становится отличной от нуля, в результате на характерные параметры задачи начинают влиять кулоновские столкновения, имеет место переходный режим обтекания зонда плазмой.
На рис.1-2 представлены входящие в оператор столкновений Фоккера-Планка функции £е, ке,
называемые потенциалами Розенблюта, для фиксированного момента времени при условии выбора начальных максвелловских функций распределения.
г»
Рис. 1. Вид функции 7е при і = 0.1 Рис. 2. Вид функции Ъе при і = 0.1
На рис. 3-4 приведены графики, иллюстрирующие динамику изменения значений концентрации ионов и электронов от расстояния до зонда для некоторого фиксированного момента времени и значения полярного угла, а также эволюцию плотности тока ионов во времени для двух различных начальных значений концентраций в невозмущенной плазме, относящихся к случаям бесстолкновительного и столкновительного режимов обтекания.
Рис. 3. График изменения зависимости концентраций частиц от расстояния до зонда: а - ионов; б - электронов без учета столкновений; в - ионов; г - электронов с учетом столкновений
Рис. 4. Эволюция плотности тока ионов во времени: а - без учета столкновений;
б - с учетом столкновений
Результаты численного моделирования, представленные на рис. 3-4, получены для фр = -10, гр = 2. В силу того, что на зонд подается отрицательный потенциал, отталкивающий электроны, вблизи зонда образуется объемный положительный заряд. О данном факте свидетельствуют равные нулю концентрации электронов и отличные от нуля концентрации ионов в окрестности г = 2. Учет столкновений влияет на профиль графиков изменения концентраций и величину объемного положительного заряда вблизи зонда. Кроме того, из результатов, приведенных на рис.4, вытекает, что при наличии столкновений процесс установления идет быстрее. Плотность тока ионов с учетом столкновений меньше, чем в случае отсутствия столкновений.
В данной работе рассмотрена задача обтекания заряженного цилиндрического тела столкновительной плазмой. Сформирована математическая модель, представляющая собой систему самосогласованных уравнений на основе уравнения Фоккера-Планка и уравнения Пуассона. Предложен метод решения поставленной задачи, разработанный на базе метода Монте-Карло, и сформировано программное обеспечение для реализации данного метода. Получены результаты численного моделирования, и проведено исследование зависимости данных результатов от характерных параметров задачи.
ЛИТЕРАТУРА
1. Пантелеев А. В., Кудрявцева И.А. Формирование математической модели двухкомпонентной плазмы с учетом столкновений заряженных частиц в случае плоского зонда//Теоретические вопросы вычислительной техники и программного обеспечения: Межвуз. сб. науч. тр. - М.: МИРЭА, 2006, с. 11 - 21.
2. Алексеев Б.В., Котельников В.А. Зондовый метод диагностики плазмы. - М.: Энергоатомиздат, 1989.
3. Олдер Б. Вычислительные методы в физике плазмы. - М.: Мир, 1974.
4. Montgomery D.C., Tidman D.A. Plasma kinetic theory. - New York, 1964.
5. Кудрявцева И.А., Пантелеев А.В. Применение метода Монте-Карло для анализа поведения двухкомпонентной плазмы с учетом столкновений между заряженными частицами//Теоретические вопросы вычислительной техники и программного обеспечения: Межвуз. сб. науч. тр. - М.: МИРЭА, 2008, с. 122-128.
6. Киреев В.И., Пантелеев А.В. Численные методы в примерах и задачах. - М.: Высшая школа, 2006.
7. Семенов В. В., Пантелеев А. В., Руденко Е.А., Бортаковский А. С. Методы описания, анализа и синтеза нелинейных систем управления. - М.: МАИ, 1993.
THE MODELING OF COLLISIONAL PLASMA DYNAMICS NEAR THE CYLINDRICAL PROBE
Kudryavtseva I.A., Panteleyev A.V.
The mathematical model describing the behavior of two-component plasma with Coulomb interactions is proposed. The model is based on Poisson and Fokker-Planck equations. The results of modeling of the two-component plasma with collisions of charged particles are given. These results agree with previous research.
Сведения об авторах
Кудрявцева Ирина Анатольевна, окончила МАИ (2005), аспирантка МАИ им. С. Орджоникидзе, автор 16 научных работ, область научных интересов - физика плазмы, зондовые методы исследования.
Пантелеев Андрей Владимирович, 1955 г.р., окончил МВТУ им. Н.Э. Баумана (1978), доктор физико-математических наук, профессор, заведующий кафедрой математической кибернетики МАИ им. С. Орджоникидзе, автор более 110 научных работ, область научных интересов - синтез оптимальных систем управления, методы оптимизации, численные методы.
