I. Введение

Развитие космической техники предполагает расширяющееся применение электроракетных двигателей (ЭРД) в составе двигательных установок для космических аппаратов. Основными причинами для этого стали значительно возросшие требования к эксплуатационным характеристикам космических аппаратов, в том числе увеличение срока активного существования. От других типов двигателей ЭРД отличаются высоким значением средней скорости истечения рабочего тела и, соответственно, высоким удельным импульсом. Это повышает эффективность использования рабочего тела, что дает заметную экономию массы рабочего тела и двигательной установки при больших сроках активного существования космического аппарата.

Наибольшие успехи в России достигнуты в разработке и практическом применении электроракетных двигательных установок на базе ускорителей плазмы с замкнутым дрейфом электронов - холловских двигателей. За почти полувековую историю разработки холловских двигателей был пройден большой путь от идей до их воплощения в конкретные изделия, многие из которых, например СПД-70, СПД-100 (СПД - стационарный плазменный двигатель) успешно эксплуатируются и имеют высокий уровень эффективности. Тем не менее, при создании холловских двигателей нового поколения возникает ряд проблем по обеспечению их эффективности, которые не удаётся решить только инженерно-экспериментальными методами. Поэтому весьма актуальным является теоретическое исследование холловских двигателей, в том числе при помощи численного моделирования динамики плазмы.

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

II. Принцип работы холловского двигателя

На рис. 1 изображена упрощенная схема холловского двигателя (ХД). Холловский двигатель состоит из анодного блока и катода-компенсатора. Корпус анодного блока выполнен из магнито-проводящей стали. Внутри корпуса расположены электромагнитные катушки. Данная магнитная система, создаёт в канале двигателя радиальное магнитное поле. Между анодом и катодом-компенсатором приложена разность потенциалов, называемая напряжением разряда. Эмитируемые катодом-компенсатором электроны, попадая в область скрещенных электрического и магнитного полей, начинают дрейфовать в канале ХД в азимутальном направлении, ионизируя атомы нейтрального газа, поступающие со стороны анода в канал ХД. Образовавшиеся ионы ускоряются электрическим полем объёмного заряда электронов. Импульс передается на конструкции ХД вследствие взаимодействия азимутального холловского тока электронов с внешним магнитным полем. Электроны в

Численная модель динамики плазмы в холловском

двигателе

Иришков С.В. (serge.irish@mail.ru)

ФГУП Центр Келдыша , Москва, Россия

Нейтральный газ

канале ХД замагничены - параметр Холла для них сот ст >> 1. В тоже время магнитное поле не оказывает существенного влияния на движение ионов.

III. Описание численной модели

A. Описание области моделирования

В качестве моделируемого двигателя был выбран ХД малой мощности1, в котором разрядная камера имеет прианодное расширение и изготовлена как из керамики, так и из металла (см. рис. 2). Современные зондовые измерения параметров плазмы в ХД показывают, что падение потенциала между анодом и катодом имеет место не только в канале ХД, но и в области за срезом ХД2. Также влияние на разряд оказывает плазма, примыкающая к катоду. Положения и параметры анодной и катодной плазмы представляют большой интерес, поскольку позволят определить структуру области ускорения ионов. Таким образом, предполагается, что модель должна описывать помимо области канала ХД и область, примыкающую к срезу ХД, как это показано на рис. 2. В моделируемом ХД диаметр средней линии канала достаточно мал, что требует решать задачу в двумерной осесимметричной постановке.

Магнитные катушки

Рисунок 1. Модель ХД.

B. Описание частиц плазмы

В данной модели рассматриваются следующие компоненты плазмы разряда в ХД - электроны, ионы и нейтральные атомы. Функция распределения каждой компоненты по энергиям (ФРЭ) является двумерной в пространстве координат и

трехмерной в 2D3V.

ФРЭ для

удовлетворяет описывающими ионизации:

f dt

пространстве скоростей

Рисунок 2. Область моделирования.

нейтральных атомов fn

кинетическому уравнению, дополненному модельными слагаемыми, процесс ион-нейтрал рекомбинации и гибель нейтралов вследствие

>f

n + Пп <0 nPion

> f

и

где f on - ФРЭ ионов, nn, n , ne - концентрация, соответственно, нейтралов, ионов электронов, оП - сечение ион-нейтрал перезарядки, о j - сечение ионизации, ue, uion, и п -скорости электронов, ионов и нейтралов.

В модели учитываются только однократно заряженные ионы. Кинетическое уравнение для ионов записывается аналогично уравнению (1):

dt дх}. Мах 1 doimj

где E - вектор напряженности электрического поля, В - вектор индукции магнитного поля, e - заряд электрона, Ма - масса атома нейтрального газа.

ФРЭ электронов fe удовлетворяет кинетическому уравнению, дополненному модельными слагаемыми, описывающими рождение электронов, а также изменение их ФРЭ при столкновениях с нейтральными атомами:

f +Oe,. f (E} X В ]] )-dfe =

где fe - ФРЭ вторичных электронов, штрихами обозначены ФРЭ первичных электронов после соответствующих процессов, oi, Оо, оВ - соответственно, сечение ионизации, сечение упругих столкновений электронов с нейтральными атомами и сечение возбуждения электронами нейтральных атомов, I - интеграл кулоновских столкновений электронов. В уравнениях (1-3) усреднение производится по ФРЭ соответствующих частиц.

C. Граничные условия

В модели различаются пять типов границ области решения: 1) свободная граница, 2) непроницаемая диэлектрическая стенка (керамика), 3) непроницаемая металлическая стенка, 4) анод ХД, 5) ось симметрии. При пересечении оси симметрии все частицы отражаются от нее зеркально, а производная потенциала плазмы на ней равна нулю. Частицы могут покидать область моделирования через свободные границы.

В качестве граничного условия на непроницаемой стенке для нейтралов используется условие диффузного отражения с равновероятно направленной скоростью в полупространстве, определяемой температурой стенки.

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

Для электронов задаются следующие граничные условия. При пересечении границы, совпадающей с анодом, электрон уничтожается, а переносимый заряд усредняется по ионному временному интервалу и составляет ток разряда. По достижении непроницаемых стенок электроны уничтожаются и по коэффициенту вторичной электронной эмиссии

ks = (Ee / E1 )а (где Ee - энергия электрона, E1 - первый порог размножения, а -

экспериментальный коэффициент) рождается требуемое количество вторичных электронов.

На металлической непроницаемой поверхности пристеночный слой падения потенциала определяется самосогласованно. По суммарному заряду, накопленному на поверхности, и её ёмкости определяется плавающий потенциал поверхности, который затем используется в качестве граничного условия при решении уравнения Пуассона.

На диэлектрической непроницаемой стенке считается, что нормальная составляющая электрического поля равна нулю - слой пристеночного падения потенциала напрямую не моделируется. Однако пристеночный скачок потенциала находится в приближении плоского конденсатора по локально накопленному заряду и используется в качестве граничного условия для электронов. Если нормальная к стенке энергия электрона меньше пристеночного скачка потенциала AUw, то он отражается от нее зеркально, в противном случае считается,

что электрон достиг стенки.

Для моделирования катодной плазмы используется модель квазинейтральности. На каждом шаге итераций на границах области моделирования, совпадающих с катодной плазмой, загружается количество электронов необходимое для того, чтобы выполнялось условие Hne * ILn.. Если имеется избыток электронов Hne > ILn., то электроны не загружаются. Положение загружаемых электронов определяется случайно, а их температура Te-cat задается как параметр решения. Считается, что Te-cat * 2 эВ.

электронным

перезарядки

моделируются

D. Моделирование столкновений

Основные параметры плазмы в моделируемом ХД показаны в таблице 1. Это позволяет провести анализ длин свободного пробега частиц и понять, какие типы столкновений учитывать в модели, а какими можно пренебречь (см. таблицу 2).

Столкновения частиц плазмы моделируются методом Монте-Карло. Используется экспериментально измеренные сечения упругого рассеяния электрона на атоме ксенона3,

возбуждения атома ксенона

электроном и ионизации

ударом5. Процесс ионов также методом Монте-Карло. Считается, что в результате перезарядки частицы иона и нейтрала обмениваются

скоростями. Сечение перезарядки ионов взято из6. Кулоновские взаимодействия электронов и ионов моделируются при помощи метода Монте-Карло. При этом модель рассматривает кулоновские взаимодействия как парные столкновения с рассеянием на большой угол.

Таблица 1. Основные параметры плазмы ХД.

IV. Численный алгоритм

В численном алгоритме модели для решения кинетических уравнений 1-3 используется метод частиц в приближении частица-сетка. После установления потока нейтрального газа в область моделирования загружается максвеловская плазма для поджига разряда. Затем алгоритм переходит к основному итерационному циклу. Заряды частиц раздаются в узлы

сетки при помощи метода частиц в ячейках . Используя сеточные плотности зарядов, методом последовательной верхней релаксации по пятиточечному шаблону решается уравнение Пуассона. Затем на сетке вычисляется электрическое поле. На следующем шаге находятся новые ФРЭ частиц. В местоположение каждой частицы интерполируются сеточные значения полей и вычисляется сила, действующая на частицу. Затем частица перемещается и к ней применяется граничные условия. Для интегрирования уравнений

движения используется метод leap-frog в трехмерной осесимметричной реализации Boris7. Затем методом Монте-Карло разыгрываются столкновения частиц. Таким образом, вносятся изменения в ФРЭ частиц, связанные со столкновениями. Далее загружаются частицы, необходимые для подачи нейтрального газа и поддержания катодной плазмы. В модели используется ортогональная сетка с постоянными шагами hz и hr. Считается, что магнитное поле задано в области моделирования и стационарно в процессе решения. Т.е. действием холловского тока, индуцированного замкнутым движением электронов в разряде, пренебрегается по сравнению с внешним магнитным полем. Учитываются две компоненты магнитного поля - радиальная и аксиальная.

A. Ускорение расчета

В модели частица-сетка метода частиц необходимо отслеживать вариации полей на расстояниях не более радиуса Дебая, что при характерных параметрах плазмы для моделируемого ХД (см. табл. 2) потребует сетку порядка 1200 на 600 узлов. Тем не менее, данная сложность обходится путем увеличения константы диэлектрической проницаемости вакуума в т раз е 0 =£ 0т 2. Это приводит, во-первых, к увеличению радиуса Дебая, что

позволяет использовать более редкую сетку. Во-вторых, при введении в расчет т , замедляется динамика плазмы, что позволяет использовать больший шаг интегрирования по времени.

B. Тестирование алгоритма

Для обеспечения корректности вычислений необходимо провести тестирование каждой составляющей алгоритма на простых задачах. Также требуется провести комплексное тестирование алгоритма, включающее проверку общей концепции модели и взаимодействия компонент алгоритма.

Алгоритм численного решения уравнения Пуассона был протестирован на аналитической периодической функции. При этом максимальная ошибка не превышала 4%. Аналогичное тестирование подпрограммы вычисления градиента дало максимальную ошибку 10%. Подпрограмма интегрирования уравнений

11 11 J 1 Ое+0 1е-3 2е-8 Зе-G 4e-G 5e-G Ge-S

движения методом leap-frog проверялась путем моделирования движения одиночного электрона в скрещенных магнитном и электрическом полях. Было получено циклотронное движение с частотой равной теоретической (см. рис.3). При этом радиус орбиты с точностью в третьем знаке соответствовал аналитическому значению радиуса Лармора. Важным с точки зрения моделирования является сохранение магнитного момента при движении электрона в реальном магнитном поле ХД в отсутствии электрического поля. Сохранение магнитного момента лежит в основе такого важного явления, как эффект запирания электрона в магнитной ловушке. В точках разворота электрона продольная составляющая его скорости стремится к нулю, т.е. должны совпадать значения поперечной скорости электрона и, следовательно, значения напряженности магнитного поля. Точкам разворота соответствуют максимумы напряженности магнитного поля на рис. 3.

!..........Ч............Ц......

GE-Q5 ВреМЯ, МС

рис.3 Тестирование подпрограммы интегрирования уравнений движения.

B(t)

Легко видеть, что средние значения напряженности магнитного поля в точках разворота электрона совпадают с хорошей точностью.

Так же было проведено комплексное тестирование численного алгоритма. Распределение параметров при истечении нейтрального газа в вакуум сравнивалось с результатами расчета других программ. При моделировании самосогласованной плазмы в закрытой системе наблюдались плазменные колебания, частота которых соответствует теоретической. Также наблюдался эффект численного разогрева плазмы8. Модель правильно описывает слой пристеночного падения потенциала на металлической стенке. Величина пристеночного падения потенциала для максвелловской плазмы соответствовала теоретической .

V. Апробация модели на холловском двигателе малой мощности

В целях верификации модели были выполнены вычислительные эксперименты по моделированию динамики плазмы в холловском двигателе КМ-37, использующим в качестве рабочего тела ксенон на нескольких режимах с мощностью разряда - от 165 до 360 Вт.

A. Интегральные характеристики ХД

Интегральные характеристики режимов, измеренные экспериментально и полученные численно показаны в таблице 3. Приведены значения удельного импульса и КПД без учета потерей мощности в магнитной системе и расхода газа в катоде. В каждом режиме достигался минимум тока разряда, путем оптимизации параметров магнитной системы. Решения для всех режимов проводились при одинаковых параметрах модели. Здесь и далее приведены значения характеристик, а также распределения параметров плазмы, полученные путем осреднения их мгновенных значений на временном интервале, охватывающем два соседних ионизационных колебания.

С повышением мощности разряда в ХД наблюдается рост КПД и удельного импульса. Данная тенденция наблюдалась и в результате численных расчетов (см. рис. 3, 4). Тем не менее, полученные значения тяги и удельного импульса оказались выше экспериментально измеренных в среднем на 13%. Анодный КПД оказался завышенным в среднем на 18%. Данные расхождения носят систематический характер и не влияют на тенденцию роста КПД. Причины расхождения становятся ясны при анализе структуры анодного КПД.

Таблица 3.

B. Структура анодного КПД

Без учета многозарядных ионов анодный КПД есть произведение трех коэффициентов

эффективности г|а = а-Р-у 2. Первый параметр а= mj /mа - коэффициент использования рабочего тела (КИРТ). Второй параметр р = J +/ J - коэффициент использования тока разряда. Третий параметр у = J m <uj >2 /2eUp - коэффициент потерь скорости ионов (эффективность ускорения). Здесь <и j > - средняя осевая скорость ионов, Up - напряжение

разряда, Jp - ток разряда, J+ - ток пучка ионов, mj, mа - секундный расход ионов и

нейтралов. Коэффициенты эффективности были измерены экспериментально согласно

методике описанной в работе .

Заметное расхождение между экспериментальными и расчетными коэффициентами эффективности видно для в . Его расчетные значения близки к единице, тогда как для ХД КМ-37 значение в обычно лежит в пределах 0.6 - 0.7. Высокие значения в объясняются тем, что в алгоритме используется ёмкостная модель определения плавающего потенциала металлических стенок. В силу использования конечной величины ёмкости имеет место инерционность изменения плавающего потенциала стенки относительно потенциала плазмы. Таким образом, за период одного ионизационного колебания имеют место повышенные потери электронов на стенках, связанные с необходимостью установления слоя пристеночного падения потенциала. Это снижает ток разряда и увеличивает значение в. Следует отметить, что характер зависимости e(Up) не оказывает существенного влияния на структуру анодного КПД. Так же как и в эксперименте, в слабо зависит от напряжения разряда (см. рис. 6). Данная зависимость имеет место и в других ХД10. Завышенные значения в были получены и в работе11, описывающей полностью кинетическую модель применительно к ХД.

Как видно из рис. 5, тенденция изменения а при росте напряжения разряда повторяет экспериментальную зависимость a(Up). Завышенные по сравнению с экспериментом значения а можно отнести на счет заниженных потерь ионов на стенках из-за использования упрощенной модели пристеночного падения потенциала вблизи керамических стенок.

Из рис. 5-7 видно, что рост анодного КПД при увеличении мощности в основном связан с ростом а (зависимость монотонно возрастающая). Данная тенденция вероятней всего связана с ростом температуры электронов при увеличении напряжения разряда (см. рис. 8). Это можно объяснить при помощи следующих оценок. Объемную частоту ионизации ne nn -<ст pe > определяют два конкурирующих фактора - концентрация электронов ne и константа скорости реакции <о р e > . Константа скорости реакции <о р e > при повышении

напряжения разряда и, соответственно температуры электронов в области ионизации с Te 12 эВ (режим 250В) до Te 22 эВ (режим 500В) возрастает более чем в 2 раза, тогда как концентрация электронов с ростом напряжения разряда уменьшается лишь на 15%. Соответственно увеличивается частота ионизации нейтрального газа и а.

Другим механизмом роста анодного КПД при увеличении мощности является монотонный рост у. Это хорошо видно на экспериментальной кривой на рис.7. В численном расчете зависимость уЦр) не была монотонной.

C. Температура электронов

На рис. 11 показано осевое распределение температуры электронов на R = 18.5 мм, что соответствует средней линии канала ХД. В области прианодного расширения средняя температура электронов слабо зависит от напряжения разряда и находится в пределах 2-5 эВ. Для всех моделируемых режимов осевое положение максимума температуры находится в области максимума радиальной компоненты магнитного поля (см. рис. 11, 12). Данный характер распределения температуры электронов вдоль оси ХД подтверждается и результатами зондовых измерений13.

D. Зона ускорения ионов

Как видно из рис. 9 в области, находящейся на расстоянии 3-9 мм от анода ХД и в области находящейся на расстоянии 23-29 мм от анода, аксиальная напряженность электрического поля заметно меньше, чем в области основного падения потенциала Z = 1020 мм. Это подтверждает существование в разряде ХД анодной и катодной плазмы. Таким образом, модель показывает наличие в разряде ХД зоны ускорения, протяженность которой составляет около 12 мм. Оценка протяженности зоны ускорения в приближении

классической проводимости по формуле А.В. Жаринова LycK = jmUjye /eBVt даёт значение 10 мм.

E. Зона ионизации нейтрального газа

Осевые распределения концентрации ионов в целом схожи для всех напряжений разряда (см. рис.10). Можно выделить так называемое ndpo ионизации - область максимальной концентрации ионов. Из рис. 10 и 12 видно, что, как и в работе13 осевое положение ядра ионизации соответствует области с максимальным значением градиента магнитного поля VzB.

Корреляция между положением ядра ионизации и топологией магнитного поля в

двумерном случае выражается через проекцию градиента напряженности магнитного поля на

направление перпендикулярное силовой линии к магнитного поля : gradB\ = B х gradB /B. Из двумерных распределений видно, что положение ядра

ионизации и области максимального градиента магнитного поля совпадают не только в аксиальном, но и в радиальном направлении (см. рис. 13 и 14). Таким образом, вычислительные эксперименты объясняют экспериментально полученную в работе14 зависимость интегральных характеристик ХД от радиального положения максимума градиента напряженности магнитного поля. Наибольшая эффективность ХД достигалась, в случае если область максимального градиента расположена симметрично относительно стенок канала ХД. Поскольку в области максимального градиента магнитного поля происходит основная ионизация, то очевидно, что в этом случае минимизируются потери ионов на стенках канала ХД.

VI. Заключение

Обобщая приведенные выше результаты, можно сделать следующие выводы. Модель адекватно описывает поведение интегральных характеристик на различных режимах работы ХД. Исследование структуры анодного КПД показало, что при росте напряжения разряда имеет место монотонное увеличение коэффициента использования рабочего тела. Вероятней всего это связано с ростом температуры электронов. При этом также наблюдался рост эффективности ускорения, хотя в отличие от эксперимента он не был монотонным.

Пространственные распределения параметров плазмы ХД показали, что в канале ХД существует ядро ионизации, положение которого связано с распределением магнитного поля. Положение максимума температуры электронов в аксиальном направлении совпадает с областью максимальной напряженности магнитного поля. Протяженность зоны ускорения ионов хорошо соотносится с теоретическими оценками.

Литература

1 M.B. Belikov, O.A. Gorshkov, High-performance low power Hall thruster , AIAA-2001-3780, AIAA Joint Propulsion Conference and Exhibit, 37th, Salt Lake City, UT, July 8-11, 2001.

2 Y. Raitses, D. Staak, M. Keidar and N.J. Fisch, Electron-wall Interaction in Hall Thrusters , Report PPPL-4050, Princeton Plasma Physics Laboratory, February 2005. Физические величины , Справочник, ЭнергоАтомИздат, Москва, 1991.

4 Hayashi M. Determination of Electron-Xenon Total Exitation Cross-Sections, From Thresold to 1000-ev, from Experimental Values of Townsends a . J.Phys.D: Appl.Phys, 16:581-589, 1983.

5 Rapp, D. and P. Englander-Golden. Total Cross Sections for Ionization and Attachment in Gases by Electron Impact. I. Positive Ionization. Journal of Chemical Physics, 43(5), 1965.

6 Rapp, D. and W. E. Francis. Charge Exchange Between Gaseous Ions and Atoms. Journal of Chemical Physics, 37(11):2631-2645, 1962.

Ч. Бэдсел, А. Ленгдон, Физика плазмы и численное моделирование , Москва, ЭАИ, 1989. Р. Хокни, Дж. Иствуд, Численное моделирование методом частиц , Москва, Мир, 1987. 9 L.Jolivet, J.-F. Roussel, Effects of the secondary electron emission on the sheath phenomenon in a Hall thruster , IEEE, Vancouver, 16-20 July 2001.

10 Blinov N.V., Gorshkov O.A., Rizakhanov R.N., Shagayda A.A. Hall-Effect Thruster with High Specific Impulse . Proc. 4th Intern. Spacecraft Propulsion Conf. Sardinia, Italy, 2-9 June 2004.

11 V. Blateau, M. Martinez-Sanchez, O. Batishchev, J. Szabo, PIC Simulation of High Specific Impulse Hall Effect Thruster , IEPC-01-037, 27th International Electric Propulsion Conference, Pasadena, CA, 15-19 October, 2001.

12 Noah Warner, James Szabo, Manuel Martinez-Sanchez, Charachterization of a high specific impulse Hall Thruster Using Electrostatic Probes , IEPC-2003-082, 28th International Electric Propulsion Conference, Toulouse, France 2003.

А.И. Морозов, Ю.В. Есипчук, А.М. Капулкин, В. А. Невровский, В. А. Смирнов, Влияние конфигурации магнитного поля на режим работы ускорителя с замкнутым дрейфом электронов (УЗДП) , Журнал Технической Физики, 1972, ТомXLII, в. №3. 14 Blinov N.V., Gorshkov O.A., Shagayda A.A, Experimental Investigation of Magnetic Field Topology Influence on Structure of Accelerating Layer and Performance of Hall Thruster , IEPC-2005-033, 29th International Electric Propulsion Conference, Princeton, 2005.

Удельный импульс

1000 Ф-

Анодный КПД

эксперимент

-расчет

375 Up, [В]

Рис.3 Зависимость удельного импульса от напряжения разряда.

Рис.4 Зависимость анодного КПД от напряжения разряда.

о.э

0.8 0.7 0.6

0.4 0.3 0.2 0.1 0

-эксперимент расчет

375 Up, [В]

Рис.5 Зависимость а от напряжения разряда.

, 0.6

375 Up, [В]

Рис.6 Зависимость в от напряжения разряда.

Рис.7 Зависимость у от напряжения разряда.

Рис.8 Зависимость температуры электронов от напряжения разряда.

Концентрация плазмы

*ЮЕ 1Т

Рис.9 Осевое распределение потенциала плазмы по центру канала ХД.

расстояние от анода, мм

Рис.10 Осевое распределение концентрации ионов по центру канала ХД.

Радиальная компонента магнитного поля

расстояние от анода, мм

Рис.11 Осевое распределение электронов по центру канала ХД.

температуры

расстояние от анода, мм

Рис.12 Осевое распределение напряженности магнитного поля по центру канала ХД.