CC BY
20
УДК 533.9
КИНЕТИЧЕСКИЕ ХАРАКТЕРИСТИКИ ЗАРЯДКИ МАКРОЧАСТИЦ В ПЛАЗМЕ
С. В. Владимиров1, С. А. Майоров2
В работе на основе численного моделирования методом частиц в ячейке исследуются временные характеристики зарядки пылинок и кинетические характеристики стационарного состояния. Проведено сравнение с некоторыми теоретическими результатами (в основном диффузионно-дрейфовой моделью) и расчетами из перво-принципов - методом молекулярной динамики.
Пылевая плазма является объектом интенсивных экспериментальных и теоретических исследований последних лет [1 - 3]. В настоящей работе, продолжающей работы [4-6], на основе вычислительного эксперимента исследованы различные кинетические характеристики пылинок микронных размеров (макрочастиц), помещенных в плазму. Обычно в реальных экспериментах температура электронов определяется энерговкла дом от внешнего источника (СВЧ излучение, разряд, фотоионизация и т.д.), ионы же находятся в тепловом равновесии с холодными атомами буферного газа. В приэлектрод-ном слое дополнительно формируется также ионный поток, направленный к электроду, и макрочастица при определенных условиях может находиться в равновесии (левитировать) за счет равенства электрических и гравитационной сил. В работе рассмотрены случаи двухтемпературной, однотемпературной плазм и движущейся плазмы.
В настоящей работе на основе численного моделирования рассчитывалась временная зависимость процесса зарядки первоначально электрически нейтральной пылинки, поглощающей все попадающие в нее ионы и электроны, а также кинетические характеристики стационарного (установившегося) режима.
1 Университет Сиднея, г. Сидней, Австралия.
28. Maiorov@physics.usyd.edu.au
В основном моделирование проводилось методом частиц в ячейке. В центр вычислительной ячейки помещалась тяжелая макрочастица и решались уравнения Ньютона для системы точечных, заряженных частиц. Полагалось, что границы куба отражают частицы плазмы в соответствии с их функцией распределения по скоростям на большом удалении от пылинки, которая полагалась известной. Обычно выбиралось распределение Максвелла (температуры электронов и ионов могут быть различны). Такие граничные условия будем называть термостатирующими стенками. При моделировании движущейся плазмы для ионов использовались граничные условия другого типа, имитирующие входящий поток (сдвинутое максвелловское распределение) и поглощающая граница на выходе.
Большая часть расчетов выполнена по упрощенной модели, в которой при вычислении траекторий частиц учитывалось только взаимодействие макрочастицы с частицами плазмы и влияние экранирования. Несмотря на грубость такой модели, она позволяет достаточно точно рассчитывать кинетические характеристики макрочастицы в неподвижной плазме, что подтверждается сравнением получаемых результатов, как с имеющимися теоретическими результатами, так и с более точными расчетами методом молекулярной динамики (МД). Но для движущейся плазмы такая постановка приводит к более значительным ошибкам, и для ее проверки использовались более трудоемкие с вычислительной точки зрения МД-расчеты.
Общая постановка задачи и методы решения. Моделирование из первых принципов (ab initio) широко используется при решении разнообразных плазменных задач [7]. Основное преимущество этого подхода заключается в том, что в основу вычислительной модели закладывается только одно предположение - вид потенциала взаимодействия между частицами.
В данной работе применялось моделирование как методом МД, так и методом частиц в ячейке (PIC - particles in cell). Исследовались системы, состоящие лишь из некольких десятков тысяч частиц. Рассмотрим полностью ионизованную плазму, состоящую из ионов с массой m,, положительным зарядом е и электронов с массой те, зарядом —е. Рассматривается временная эволюция системы из пе электронов и п, ионов, заключенных в куб, в центре которого находится поглощающее сферическое тело радиуса R с зарядом Q — Zoe < 0, обладающее большой массой. Число ионов и электронов выбирались такими, чтобы в целом система являлась электронейтральной: п,- — пе + Zq = 0.
Траектории частиц определяются путем решения уравнений Ньютона:
d2rk/dt2 = Fk/mk, k = 1,2,...,а
(1)
В методе МД сила, действующая на любую из частиц в системе, определяется как сумма сил со стороны всех прочих частиц:
Здесь Гуь(£) - радиус-вектор /с-той частицы с массой тк и зарядом qk, г3(£) - радиус-вектор макрочастицы с зарядом С^, пр - общее число частиц. Сила кулоновского взаимодействия {к1 на близких расстояниях между частицами для устранения особенности заменялась силой кулоновского взаимодействия однородно заряженных взаимно проницаемых сфер очень маленького радиуса [7]. Уравнения (1), (2) решались методом Рунге-Кутты четвертого порядка точности. Макрочастица в настоящей работе полагалась бесконечно массивной - неподвижной, хотя в работе [4], при исследовании раскачки макрочастицы в ионном потоке, учитывалось ее движение.
Метод частиц в ячейке отличается от описанного выше метода МД в том, что используется более экономный способ вычисления сил, действующих на частицы (естественно, за счет их огрубления). В настоящей работе траектории электронов и ионов определялись путем решения уравнений Ньютона (1), как и в методе МД. Но действующая на частицы плазмы сила вычислялась другим способом:
Здесь гкд = гк — гд. В (3) при вычислении силы учитывается экранирующее влияние ион-электронного облака, окружающего макрочастицу. Суммирование в (3) ведется по частицам, расположенным ближе к макрочастице, чем рассматриваемая ¿-тая частица. Сила, действующая на макрочастицу, определялась по третьему закону Ньютона. Для определения степени адекватности модели проводились расчеты и сравнение результатов, получаемых методом МД (1), (2) и частиц в ячейке (1), (3). Как показывает анализ результатов, подход (2) для неподвижной плазмы совершенно оправдан. Вследствие этого появляется возможность проводить расчеты с большим числом частиц при реальном отношении масс электронов и ионов в течение длительного промежутка времени, что необходимо для вычисления флуктуационных характеристик зарядки макрочастицы.
(3)
Для движущейся плазмы этот подход приводит к несколько большим ошибкам, хотя и позволяет получать достаточно точные результаты.
При моделировании плазмы методом частиц необходимым условием адекватности модели (не всегда достаточным) является то, чтобы число частиц в моделируемом объеме было значительно больше, чем число частиц в сфере Дебая. В экспериментах по пылевой плазме характерными параметрами являются плотность ионов Л^ = 109 см~3, температура электронов Те = 1 эВ, ионов Г, = 0.025 эВ. При таких параметрах плазмы в дебаевской сфере находится примерно 106 частиц. Поскольку число арифметических операций на один шаг по времени в методе МД пропорционально квадрату числа частиц, то даже для современных параллельных компьютеров такой объем вычислений нереален. Поэтому приходится прибегать к масштабированию физических параметров и применять упрощения вычислительной модели, уменьшающие объем вычислений.
В пылевой плазме важную роль играют процессы зарядки макрочастиц, флуктуации их заряда. Одной из общих вычислительных проблем при моделировании плазмы является сильное различие масс ионов и электронов, что вызывает необходимость расчета медленных процессов, обусловленных движением ионов, с шагами разностной сетки, которые определяются быстрыми электронами. Часто возможно использование моделей, в которых отношение масс электронов и ионов увеличено для сокращения разницы в ионных и электронных характерных временах [5, 7]. Однако такой прием не годится для вычисления флуктуаций заряда макрочастиц.
В качестве основного упрощающего предположения, позволяющего рассчитать флуктуационные характеристики заряда пылинки, в настоящей работе используется сферическая симметрия плотности заряда вокруг пылинки при вычислении силы (3). Фактически при таком подходе частицы плазмы взаимодействуют между собой только посредством усредненных по радиусу флуктуаций заряда и через накопленный пылинкой заряд. Следует отметить, что такой подход уменьшает влияние слабой неидеальности ионной компоненты при рассмотрении двухтемпературной плазмы с холодными ионами.
Зарядка макрочастицы в плазме - некоторые формулы. В идеальном газе среднее число атомов испытавших столкновение с телом за время пропорционально площади поверхности, плотности потока и времени наблюдения. Соответственно, для шара в газе с максвелловским распределением атомов:
AN(t) = AvR2JMt,
(4)
где Jm — (Те/2тгт)1/2п - плотность потока частиц на поверхность в максвелловском газе, п - числовая плотность атомов. Зависимость (4) будет определять начальный участок зарядки макрочастицы, когда ее заряд мал по сравнению со средним зарядом.
Если макрочастица в плазме имеет среднее значение заряда в стационарном режиме Q, то эволюция отклонения заряда макрочастицы от среднего значения АQ(t) = Q(t)—Q обычно описывается уравнением:
dAQ/dt = -А Q/tj, (5)
где характерное время флуктуаций заряда ту для максвелловской плазмы определяется следующим выражением:
т = (irk) 1 .г,
5 \viRj 1+Т{/Те + <р1Те'
Здесь Vi = föTi/irM)1/2 - средняя тепловая скорость ионов, (р = - величина энергетического барьера. Выражение (6) справедливо для небольших отклонений от среднего значения. Обычно полагается, что это же время является характерным временем затухания в автокорреляционной функции (AQ(t -f r)AQ(t)) = (AQ2) ехр(—т/ту).
В приближении Фоккера- Планка средний квадрат флуктуаций заряда определяется отношением времени флуктуации к среднему времени между столкновениями тс частиц плазмы с макрочастицей:
а2 = <AQ2)/e2 = тс/2у. (7)
В работе [10] получена следующая аппроксимация для зависимости среднеквадратич ного отклонения заряда от среднего значения:
a^CrZ1'2, с, = 1/2, (8)
где Z - заряд макрочастицы, выраженный в единицах заряда электрона.
В работе [8] получена следующая аппроксимация для зависимости среднеквадратичного отклонения заряда от среднего значения:
-fL = l__I__(9)
t'RT,. l + T,/T„ + :f,/r,_
При аппроксимации процесса зарядки функцией:
Q(t) = Q + дд(0) ехрН/ту) (10)
выбирался промежуток времени, соответствующий не всему расчету, а его части. Из рассмотрения выбрасывался некоторый начальный участок, на котором не выполняется условие линейности (5), но в то же время начальное отклонение было значительно больше величины флуктуаций.
Однотемпературная плазма. Рассмотрим вначале случай однотемпературной плазмы гелия с z = 1, N{ — 2 ■ 1012 см~3, температурой электронов и ионов Те = Г, = 1 эВ. Электронный дебаевский радиус гре = (Те/4тт выраженный в единицах харак-
терного межионного расстояния, равен roeN— 6.6, число электронов в дебаевской сфере равно 1220. Условие идеальности плазмы хорошо выполняется как для электронов, так и для ионов.
Начальное распределение электронов и ионов по координатам выбиралось равновероятным в объеме куба. Распределение по скоростям соответствовало распределению Максвелла на бесконечности. В зависимости от начального расстояния до макрочастицы распределение Максвелла сдвигалось по модулю скорости на величину энергии взаимодействия с макрочастицей. Направление скорости выбиралось изотропным. Таким образом, формировалось начальное распределение без связанных частиц, которые при определенных условиях могут сильно искажать кинетические характеристики.
В начальный момент времени в центр куба помещалась нейтральная макрочастица заданного размера с бесконечной массой, поглощающая все падающие на нее частицы плазмы. В предыдущих работах число ионов и число электронов в системе было строго фиксировано. Поэтому из-за флуктуаций заряда макрочастицы система была электронейтральной лишь в среднем по времени (и то лишь в лучшем случае - при правильном выборе числа ионов и электронов в системе на основе предварительных расчетов). В настоящей работе впервые использован алгоритм вбрасывания частиц (взамен поглощенных), точно сохраняющий нейтральность всей системы. Взамен поглощенной частицы со случайно выбранной точки на поверхности куба вбрасывались частицы таким образом, чтобы обеспечить нейтральность всей системы, включая заряд макрочастицы. Эта процедура состояла в том, что при поглощении электрона их число в системе уменьшалось, а при поглощении иона в систему с поверхности куба вбрасывалась электрон-ионная пара. Таким образом, число ионов в системе поддерживалось постоянным, число же электронов было переменным, обеспечивающим нейтральность системы. При отражении от стенок куба использовались термостатирующие граничные условия. Такая постановка задачи позволяет самосогласованно учитывать процесс
зарядки и флуктуации заряда макрочастицы.
Рис. 1. Расчеты процесса зарядки неподвижных поглощающих сфер с радиусами Я = 0.25,0.5,1,2 м к м в однотемпературной плазме с Т{ = 1 эВ.
Рис. 2. Зависимости от времени числа поглощенных макрочастицами электронов и ионов для сфер с радиусами Л = ®.2Ь,2 мкм в однотемпературной плазме с Г, = 1 эВ.
На рис. 1 и 2 приведены зависимости от времени заряда макрочастиц и числа поглощенных макрочастицей электронов и ионов для сфер с различными радиусами. Время на всех рисунках нормировано на ионную плазменную частоту = (47ге2Л'1/М)1/'2, Т{ = ш^1 = 1.076 не, М - масса иона. Время расчета было равно ¿о = 1.3 - Ю-7 с, что значительно больше как времени ионного плазменного периода 6.7 • Ю-9 с, так и характерного времени зарядки макрочастицы. На графиках приведены расчеты только начального периода времени, сравнимого со временем флуктуации заряда.
На рис. 1 представлены результаты четырех расчетов зависимостей от времени заряда первоначально нейтральных макрочастиц - сфер с радиусами Д = 0.25,0.5,1,2 мкм. Начальное число ионов в системе полагалось равным 15000 (оно во время расчета не менялось), начальное число электронов в системе было равно также 15000 (во время расчета оно менялось), но система оставалась нейтральной.
На рис. 2 приведены зависимости числа поглощенных макрочастицей электронов и ионов от времени для сфер с радиусами И — 0.25,2 мкм. Там же штрихованными прямыми нанесены зависимости £) = SJмt^
5000 4000 3000 2000 1000 0
/ ' 2 мкм
г 0.5 - 1 О 25
0
6 8 10 ^
Рис. 3. Отклонение заряда от среднего значения при зарядке неподвижной поглощающей сферы с радиусом Я = 0.25 мкм в однотемпературной плазме с Т, = 1 эВ. а - отклонение заряда от среднего; б - аппроксимация зависимостью (7); в - отклонение от аппроксимации.
Рис. 4. Расчеты процесса зарядки неподвижных поглощающих сфер с радиусами Д = 0.25,0.5,1,2 м к м в двухтемпературной плазме с Т, = 1/40 эВ, Те = 1 эВ.
На рис. 3 приведены для сферы с радиусом Я = 0.25 мкм следующие зависимости от времени: отклонение заряда макрочастицы от среднего значения, ее аппроксимация (7) и отклонение заряда от аппроксимирующей зависимости (7). Этот график иллюстрирует способ определения среднего значения заряда макрочастицы, флуктуаций заряда и характерного времени флуктуаций. Все эти данные для однотемпературной плазмы приведены в табл. 1.
Таблица 1
Характеристики зарядки пылинки различных размеров в однотемпературной плазме 1 эВ. Приведены: время флуктуации, полученное в результате расчета и по теории для максвелловской плазмы, потенциальный барьер, средний заряд и его среднеквадратичное отклонение от аппроксимирующей кривой
В., мкм 0.25 0.5 1 2
т,нс, расчет 9.5 6.3 4.8 2.2
г, не, (4) 11.6 5.3 2.5 1.2
пот. барьер 2.9 3.3 3.6 4.0
ф, ср. заряд 501 1145 2520 5587
<т, фл. заряда 12.8 14.1 26.2 36.6
а/заряд1/2 0.57 0.42 0.52 0.49
Двухтемпературная плазма. Рассмотрим случай двухтемпературной плазмы с холодными ионами Г, = 0.025 эВ, все остальные физические и вычислительные параметры оставались прежними. Ионный дебаевский радиус гд,- = (Гг747ге2Лг2)1/2 при таких параметрах плазмы, выраженный в единицах характерного межионного расстояния, равен
1/3
= 1.06, число ионов в дебаевской сфере равно 5. Условие идеальности плазмы для ионов выполнено.
На рис. 4, как и на рис. 1, представлены результаты четырех расчетов зависимостей от времени заряда первоначально нейтральных макрочастиц - сфер с радиусами Я = 0.25,0.5,1,2 мкм.
Рис. 5. Зависимости от времени числа поглощенных макрочастицами электронов и ионов для сфер с радиусами Я = 0.25,2 мкм в двухтемпературной плазме с Г, = 1/40 эВ, Те = 1 эВ.
На рис. 5 приведены зависимости числа поглощенных макрочастицей электронов
и ионов от времени для сфер с радиусами II = 0.25,2 мкм. Там же штрихованными прямыми нанесены зависимости (¿(Ь) = БЗм^-
Интересен эффект немонотонного увеличения заряда у макрочастицы наибольшего размера. Анализ зависимости ионного тока от времени показывает, что он обусловлен задержкой в установлении экранирующего облака ионов из-за большого заряда макро частицы. Соответственно, в наибольшей мере он проявляется для больших частиц, у которых больше заряд.
Таблица 2
То же самое для двухтемпературной плазмы Т, = 0.025 эВ, Те = 1 эВ
Л, мкм 0.25 0.5 1 2
г, и с, расчет 30.8 16.8 14.9 12.2
г, не, (4) 2.7 1.2 0.57 0.28
<р, пот. барьер 2.3 2.6 2.9 3.0
<5, ср. заряд 391 909 2016 4123
¡г, фл. заряда 9.7 15.9 24.0 37.2
с / заряд1!2 0.49 0.53 0.54 0.58
Аналогично случаю однотемпературной плазмы в табл. 2 приведены результаты обработки расчетных данных. Отметим, что имеется значительное отклонение от теории (в сторону увеличения) времени флуктуаций.
Движущаяся плазма. Рассмотрим наиболее близкий большей части экспериментов случай движущейся плазмы. Были выполнены расчеты в аналогичной постановке, только ионам задавалась направленная скорость с кинетической энергией А, = 2 эВ. Выбор значения кинетической энергии ионов обусловлен тем фактом, что в этом случае кинетическая энергия ионов равна средней кинетической энергии максвелловского потока с Г, = 1 эВ:
оо оо
(К) = I 1V= 2Т.
о о
Здесь /м(") = и2ехр(—ти2/2Т), К [у) = ти2/2. Все остальные физические и вычислительные параметры были такими же, как при расчетах однотемпературной плазмы.
Были проведены расчеты зависимости заряда от времени для первоначально нейтральных макрочастиц - сфер с радиусами Н, = 0.25,0.5,1 и 2 мкм. Из сравнения результатов этих расчетов с результатами для однотемпературной плазмы (рис. 1),
процессы зарядки в однотемпературной плазме и ионном потоке практически совпадают.
Таблица 3 То же самое для движущейся плазмы с К = 2 эВ
R, MKM 0.25 0.5 1 2
t, не, расчет 10.6 7.3 3.7 2.3
<p, пот. барьер 3.0 3.4 3.7 4.1
Q,ср. заряд 516 1191 2602 5707
а, фл. заряда 7.4 19.6 21.6 34.8
а/заряд1/2 0.33 0.57 0.43 0.46
В табл. 3 приведены результаты обработки расчетных данных. В этой таблице значения флуктуационного времени по теории не приводятся, так как плазма не максвел-ловская.
Отметим, что имеется хорошее совпадение характеристик движущейся плазмы и однотемпературной. Причем для макрочастиц большего размера оно лучше.
Выполненное численное моделирование является попыткой применения методов численного моделирования для изучения кинетических характеристик макрочастиц в плазме. Разработана методика численного моделирования методом частиц, результаты находятся в согласии с теоретическими и экспериментальными данными, получены временные характеристики установления заряда.
Авторы благодарят Австралийский исследовательский совет (ARC) за финансовую поддержку работы.
ЛИТЕРАТУРА
[1] Т о m a s Н. and М о г f i 1 1 G. Е. Nature (London), 379, 806 (1996).
[2] Ц ы т о в и ч В. Н. УФН, 167, N 1, 57 (1997).
[3] И г н а т о в А. М. Физика плазмы, PLASMA PHYS. REP., 24(8), 677 (1998).
[4] В л а д и м и р о в С. В., Крамер Н. Ф., Майоров С. А. Краткие сообщения по физике ФИАН, N 7, 7 (2000).
[5] М a i о г о V S. А., V 1 a d i m i г o v S. V., and Cramer N. F. Phys. Rev. E, 63, 17401 (2001).
[6] V 1 a d i m i г о V S. V., М a i о г о v S. A., and Cramer N. F. Phys. Rev. E, 63, 45401 (2001).
[7] X о к н и Р., Иствуд Дж. Численное моделирование методом частиц. М., Мир, 1987.
[8] Mutsoukas Т. and Rüssel М. Phys. Rev. Е, 55, N 1, 55 (1997).
[9] К h г а р a k S.A., Nefedov А. P., Pet го v О. F., and Vau Ii na О. S. Phys. Rev. E, 59, 6017 (1999).
[10] С u i С. and G o r e e G. IEEE Trans. Plasma Sei., 22, 151 (1994).
Институт общей физики РАН Поступила в редакцию 4 февраля 2002 г.
