Содержание к диссертации
Введение
1 История разработки и исследования ионных двигателей 9
1.1 История применения ионных двигателей 9
1.2 Ионный двигатель и процессы, протекающие в его газоразрядной камере 9
1.3 Моделирование как метод исследования 12
1.4 Методы моделирования динамики плазмы 13
1.5 Моделирование плазмы в газоразрядной камере ионного двигателя 15
1.6 Моделирование близких по параметрам плазменных устройств 20
1.7 Основные результаты обзора работ по моделированию ИД 21
2 Постановка модели 24
2.1 Выбор методики моделирования 24
2.2 Область моделирования 25
2.3 Моделируемые процессы 26
2.4 Плазма в электростатическом приближении 26
2.5 Столкновения между частицами 28
2.6 Кинетические уравнения компонент плазмы 29
2.7 Полная система уравнений 31
2.8 Граничные условия 31
2.9 Методика задания начальных условий 32
2.10 Методика упрощения кинетического моделирования плазмы 35
2.11 Обобщение постановки модели 37
3 Численная реализация модели 38
3.1 Область устойчивости решения 38
3.2 Общий алгоритм моделирования 40
3.3 Задание параметров моделирования 40
3.4 Построение расчетной области 41
3.5 Алгоритм шага по времени 42
3.6 Взаимосвязь локальных и сеточных величин 44
3.7 Расчет потенциала электрического поля 46
3.8 Интегрирование движения макрочастиц 49
3.9 Методика индивидуального подбора шага по времени для частиц 51
3.10 Моделирование столкновений 53
3.11 Взаимодействие частиц с границами области моделирования 57
3.12 Моделирование катода 58
3.13 Учет влияния источника разряда 61
3.14 Обобщение численной реализации модели 61
4 Тестирование численных алгоритмов 63
4.1 Тестирование алгоритмов сеточной интерполяции 63
4.2 Тестирование алгоритмов расчета электрического поля 64
4.3 Тестирование алгоритмов движения частиц 67
4.4 Тестирование алгоритмов взаимодействия частиц с границами области моделирования 69
4.5 Тестирование алгоритмов взаимодействия между частицами 71
4.6 Обобщение результатов тестирования вычислительных алгоритмов 72
5 Результаты моделирования 73
5.1 Ионный двигатель ИД-50 73
5.2 Параметры и схема области моделирования для ГРК ИД-50 76
5.3 Результаты моделирования ГРК ИД-50 77
5.4 Распределения параметров плазмы в ГРК ИД 83
5.5 Распределение электронов по энергиям 90
5.6 Динамика ионной компоненты 98
5.7 Обобщение результатов моделирования 1 6 Заключение 103
7 Список литературы 106
- Моделирование как метод исследования
- Кинетические уравнения компонент плазмы
- Построение расчетной области
- Тестирование алгоритмов взаимодействия частиц с границами области моделирования
Моделирование как метод исследования
В ионных двигателях такого типа ионизация рабочего тела осуществляется электронным ударом. Электроны, эмитированные катодом ГРК и ускорившиеся в приложенной разности потенциалов, обладают достаточной энергией для ионизации нейтральных атомов ксенона. Повышение эффективности ионизации осуществляется путем удержания электронов внутри ГРК с помощью приложенного магнитного поля, для поддержания которого используются либо магнитные катушки, либо постоянные магниты. Конфигурация магнитного поля обычно выбирается таким образом, чтобы наиболее эффективно препятствовать транспорту электронов из середины газоразрядной камеры на аноды.
Для получения тяги ионы, образованные в ГРК, ускоряются электрическим полем между перфорированными электродами ионно-оптической системы. Внутренний электрод, называемый эмиссионным, поддерживается под катодным потенциалом. Это обеспечивает отражение большей части электронов обратно в плазму и вытягивание положительно заряженных ионов. Второй электрод (ускоряющий) имеет высокий потенциал. В электрическом поле между двумя электродами осуществляется ускорение ионов до высоких скоростей.
Трех-электродная конфигурация ионно-оптической системы предусматривает дополнительный замедляющий электрод. Функция этого электрода заключается в предотвращении распыления ускоряющего электрода обратным потоком ионов, образованных в результате перезарядки первичных ионов на нейтральных атомах [10]. Несмотря на это, эрозия электродов остается основным фактором, ограничивающим ресурс ионного двигателя. Одним из методов увеличения ресурса двигателя наряду со снижением ускоряющего напряжения и использованием замедляющего электрода является выполнение элементов ионно-оптической системы из устойчивых к распылению материалов: молибдена, углерод-углеродных композитов или пиролитического графита.
В ионных двигателях с разрядом постоянного тока в качестве источников электронов, инициирующих и поддерживающих разряд, используется, как правило, полый катод [11]. Это обусловлено его эффективностью и надежностью. Полый катод позволяет создать наибольшую плотность тока при меньших энергетических затратах. Недостаток полых катодов, связанный с дополнительным расходом рабочего тела, в данном случае не играет роли. Ксенон, затраченный на работу катода ГРК, затем оказывается внутри газоразрядной камеры и вносит свой вклад в тягу двигателя. Второй катод (катод-нейтрализатор) располагается снаружи ГРК и служит для компенсации положительного заряда ионного пучка.
Численная модель – это компьютерная программа, работающая на отдельном компьютере, суперкомпьютере или множестве взаимодействующих компьютеров (вычислительных узлов), реализующая свойства объекта или понятия в форме, отличной от реальной, но приближенной к алгоритмическому описанию, включающей набор данных, характеризующих свойства системы и их динамику во времени [12].
В последнее время компьютерное моделирование стало обычным инструментом исследования слишком сложных для полного аналитического описания систем во всех отраслях науки и техники. Вычислительные эксперименты, являющиеся заключительным этапом подобных работ, позволяют ясно увидеть все детали и особенности процессов, протекающих в объекте исследования. Особенно значимым это оказывается в тех случаях, когда постановка и проведение эксперимента являются трудоемкими, дорогостоящими или вовсе невозможными. Так вычислительный эксперимент позволяет заглянуть в физику процессов, не поддающихся прямому измерению, либо проверить новаторские идеи, существующие лишь на бумаге.
Плазма как совокупность частиц, проявляющих коллективные свойства, требует особых подходов к моделированию. При этом важно понимать, что выбор методики моделирования той или иной системы основывается как на свойствах самой системы, так и на тех процессах и явлениях, которые требуется исследовать. Симуляция, создаваемая в процессе моделирования, не является полной копией исследуемой системы. Однако некоторые ее параметры могут совпадать или коррелировать с аналогичными параметрами в реальности. Наша задача заключается в том, чтобы построить модель таким образом, при котором создаваемые симуляции, будучи неизбежно несовершенными в некоторых, а то и во многих аспектах, тем не менее позволят исследовать интересующие нас свойства подлинной системы.
В общем случае различают аналитическое и имитационное моделирование. При аналитическом моделировании изучаются математические (абстрактные) модели реального объекта в виде алгебраических, дифференциальных и других уравнений, предусматривающих осуществление однозначной вычислительной процедуры, приводящей к их точному решению. При имитационном моделировании исследуются математические модели в виде алгоритмов, воспроизводящих функционирование исследуемой системы путем последовательного выполнения большого количества элементарных операций. Забегая вперед, можно отметить, что модель, которой посвящена данная работа, относится ко второму типу.
Для описания плазмы как совокупности большого числа частиц применяются уравнения Больцмана, дополненные уравнениями Максвелла. Однако прямое полностью кинетическое моделирование подобной системы потребовало бы невероятно большой вычислительной мощности. Поэтому для моделирования плазмы прибегают к различным допущениям и упрощениям, позволяющим строить модели, приемлемые для реализации.
Для нашего класса исследовательских задач, как правило, целесообразно не рассматривать электромагнитные волны. Это справедливо, поскольку для нас важны процессы, протекающие на временах, существенно превышающих время прохождения света через систему. Иными словами, предполагается, что электрическое и магнитное поля мгновенно подстраиваются под текущее распределение зарядов и токов в системе.
Различают три основных подхода к моделированию плазмы: магнитогидродинамические модели, электростатические модели и гибридные модели. В магнитогидродинамическом приближении плазма полагается квазинейтральной, а все пространственные масштабы существенно превышают дебаевский радиус экранирования. Зарядовая квазинейтральность исключает из рассмотрения механизм, ответственный за возбуждение плазменных колебаний и высокочастотного движения частиц.
Таким образом, описание компонент плазмы кинетическими уравнениями заменяется на уравнения магнитной гидродинамики. Основными уравнениями магнитогидродинамических моделей являются уравнения сохранения массы, импульса и энергии для всех рассматриваемых компонент плазмы, а также уравнения Максвелла. Для решения этих уравнений зачастую применяются стандартные конечно-разностные методы [13].
Электростатическое приближение, в отличие от магнитной гидродинамики, рассматривает разделение зарядов. Модели плазмы в этом приближении отражают процессы, протекающие на малых пространственных и временных масштабах. Эти модели могут рассматривать изменение параметров на длинах порядка радиуса Дебая и с частотами порядка плазменной частоты.
Обычно предполагается, что токи, связанные с разделением зарядов, малы. Это позволяет исключить нестационарные правые части уравнений Максвелла. В моделях с электростатическим приближением плазма описывается системой кинетических уравнений. Существует несколько наиболее распространенных методов для их решения.
Конечно-разностный метод решения кинетических уравнений предполагает замену уравнений на конечно-разностные, которые решаются как система алгебраических уравнений. Однако существенным недостатком этого метода является то, что переход к конечно-разностной форме приводит к сеточной диффузии в пространстве скоростей. Это ограничивает время, в течение которого результаты моделирования остаются корректными.
Кинетические уравнения компонент плазмы
Процессы ионизации, рекомбинации ионов, кулоновские взаимодействия в плазме описываются при помощи столкновений между частицами. В плазме газоразрядной камеры ионного двигателя длины пробега частиц достаточно велики по сравнению с размерами области моделирования. Другими словами, действие на частицу суммарной силы со стороны окружающих ее компонент плазмы в среднем превосходит взаимодействие этой частицы с ближайшими соседями. Это дает возможность решать кинетические уравнения в приближении Власова. Столкновения между частицами можно моделировать отдельно, при помощи одного из методов статистических испытаний. В данной модели используется метод Монте-Карло. Этот метод весьма эффективен при сложной геометрии и наличии внешних полей, что является особенностью данной моделируемой системы.
В применении к методу частиц метод Монте-Карло заключается в вероятностном моделировании отдельных актов взаимодействий по известным сечениям реакций. С помощью псевдослучайных последовательностей производится розыгрыш типов столкновений, углов рассеяния, обмен энергией и т.п. Подробно реализация этого метода будет описана в разделе 3.10.
Анализ работ, посвященных моделированию плазмы в газоразрядной камере ионного двигателя, позволяет выделить типы взаимодействий, оказывающих наибольшее влияние на формирование и динамику разряда. В данной модели были учтены следующие взаимодействия: 1. Кулоновские столкновения электронов; 2. Упругие столкновения электронов с нейтральными атомами; 3. Неупругие столкновения электронов с нейтральными атомами; 4. Реакции ионизации нейтральных атомов электронным ударом; 5. Реакции перезарядки ионов с нейтральными атомами. Как показывают оценки, известные механизмы электронной проводимости, основанные на столкновениях с тяжелыми частицами (ионами и атомами нейтрального газа), неудовлетворительно описывают электронный ток в плазме с характерными для данной задачи концентрациями компонент и величинами магнитного поля. Как правило, величина электронного тока оказывается существенно выше того значения, которое предсказывает классическая теория столкновительной диффузии электронов поперек магнитного поля. Такая проблема носит название проблемы аномального транспорта и имеет место не только в плазме газоразрядной камеры ионного двигателя, но и во многих других типах устройств, в которых тем или иным способом ограничивается подвижность электронов при помощи магнитного поля. Анализ этого явления впервые проводился в классической работе [63] и в последствии получил название бомовской диффузии. Моделирование аномальной электронной проводимости осуществлялось введением дополнительных упругих столкновений электронов так, чтобы с помощью столкновительной диффузии воспроизвести эффект бомовской диффузии, которую напрямую учесть в двухмерной модели невозможно. Как известно, выражение для коэффициента бомовской диффузии имеет вид [63]: кТе еВ где коэффициент g определен эмпирически как 1/16 [6]. Коэффициент классической диффузии поперек магнитного поля имеет вид: 1 2 1 3 со] где vTe - тепловая скорость электронов, veN - частота столкновений электронов с нейтралами, а сос - циклотронная частота. Заменяя в модели бомовскую диффузию дополнительной столкновительной, можно выразить величину добавочной частоты столкновений vB : кТе 1 v2yn g А TeVeN еВ 3 со] Подставив известные выражения для средней тепловой скорости электронов и циклотронной частоты, получим: еВ vB = g— . me
При моделировании плазмы рассматривались три компоненты: нейтральные атомы, электроны и ионы. В модели считалось, что функции распределения всех компонент являются двухмерными в пространстве координат и трехмерными в пространстве скоростей. Нейтральный газ поступает в область моделирования через газораспределитель и через катод. Распределение этих частиц по скоростям в модели полагалось максвелловским. Температуры и интенсивность каждого из потоков задавались как параметры модели. Атомы нейтрального газа могут покидать область моделирования через ионную оптику с некоторой вероятностью.
Поскольку в данной модели используется приближенный интеграл столкновений, соответствующий релаксационному приближению, то с учетом ионизационных столкновений и столкновений перезарядки можно записать следующее кинетическое уравнение для нейтральной компоненты: dfn dfn 1Г + vnj = -ne( TiVe)fn + nn{arVi)fi - щ{огУі)Гп Ot OX; где (ТІ - сечение ионизации, ar - сечение перезарядки, ve - скорость электронов, v i скорость ионов, vn - скорость нейтральных атомов. Усреднение производится по функции распределения соответствующих частиц. Для описания функции распределения ионов также необходимо записать кинетическое уравнение. Однако, в отличие от нейтралов, на ионы оказывает воздействие электрическая сила. В данной модели влиянием магнитного поля на движение ионной компоненты пренебрегалось, поскольку в характерном для ГРК ИД магнитном поле радиус Лармора для ионов на порядки превышает размеры системы. Кинетическое уравнение для ионов с учетом ионизационных столкновений и столкновений перезарядки имеет следующий вид:
Построение расчетной области
Схема общего алгоритма моделирования представлена на Рис. 6. На первом этапе задаются параметры моделирования и области, строится геометрия ГРК, загружается распределение магнитного поля. Подготовка к моделированию предполагает операции по построению расчетной сетки, предварительному расчету необходимых величин и другие операции. Задание начальных условий, приближенных к ожидаемому стационарному решению, позволяет сократить время моделирования и избежать разрушения решения вследствие стартового заброса параметров плазмы в области моделирования. Основная фаза расчета предполагает непосредственное моделирование плазмы в ГРК и продолжается до получения стационарного решения. Подробная реализация каждого из этапов и отдельных их элементов будет описана в следующих разделах.
Перед началом моделирования необходимо задать все требуемые параметры, которые делятся на две группы: параметры исследуемого двигателя и параметры математической модели. К важным для данной задачи параметрам двигателя относятся: 1. Электрические потенциалы различных границ области моделирования; 2. Коэффициенты прозрачности эмиссионного электрода для ионов и нейтралов; 3. Величины расхода нейтрального ксенона в катод и газораспределитель. Параметры моделирования необходимо задавать с учетом условий обеспечения устойчивости решения. В данной работе применялась следующая схема.
Шаг пространственной сетки определялся исходя из возможностей вычислительной машины. Разумеется, более мелкая сетка предпочтительна с точки зрения точности решения, однако увеличение числа ячеек и частиц требует больших вычислительных ресурсов.
Остальные параметры (длительность шага по времени и величина коэффициента г ) рассчитывались на основании ожидаемых величин концентрации и температуры. Коэффициент г выбирался таким, чтобы удовлетворить первому критерию устойчивости. Здесь вместо коэффициента взята величина 2, чтобы отойти от границы устойчивости и иметь некоторый запас [69]:
Для удовлетворения второго критерия устойчивости необходимо выбрать величину шага по времени: „ 2 Т . со , Однако третий критерий (разрешимость траекторий быстрых частиц) накладывает на шаг по времени более сильное условие: Т - Н . (3.2) 2Ude me Здесь Ud - напряжение разряда. Это условие обеспечивает то, что наиболее быстрые электроны (ускоренные полным разрядным напряжением) за один шаг по времени будут смещаться не более чем на половину ячейки.
Область моделирования задается замкнутым набором отрезков, не имеющих пересечений. Каждый элемент границы характеризуется координатами своих концов и типом, к которому относятся границы. Также на этом этапе необходимо произвести загрузку заранее рассчитанной конфигурации магнитного поля.
После задания геометрии области выполняется подготовка ее к процедуре моделирования. Осуществляется анализ взаимного расположения узлов расчетной сетки и границ области моделирования. В результате определяется тип каждого узла: внутренний, внешний либо граничный. Также осуществляется создание дополнительных вспомогательных узлов сетки во всех точках пересечения границ области моделирования с плечами сетки между основными узлами. В дальнейшем эти вспомогательные узлы будут использоваться для расчета распределения потенциала электрического поля в области моделирования. На Рис. 7 показана схема участка области моделирования после процедуры анализа.
На этом рисунке желтым цветом обозначена область моделирования, зеленым цветом – внутренние узлы, фиолетовым – граничные узлы, черным – наружные узлы. Красным цветом обозначены дополнительные вспомогательные граничные узлы. Тип каждого узла закрепляется за ним в начале процесса моделирования и используется в дальнейшем при расчете распределения потенциала электрического поля, плотности компонент плазмы и других сточных величин.
Основной этап моделирования представляет собой последовательное циклическое выполнение шагов по времени. На каждом шаге осуществляется расчет малых перемещений частиц в самосогласованном электромагнитном поле с учетом взаимодействий со стенками и другими частицами. Схема численного алгоритма представлена на Рис. 8. Рис. 8 Алгоритм шага по времени
Для данного алгоритма полагается, что в начальный момент времени текущее распределение концентраций частиц в области моделирования известно. Для первого шага оно определяется в процессе загрузки стартового распределения частиц в область моделирования. Для последующих шагов определяется в конце предыдущего шага.
На схеме под «активными» частицами подразумеваются те, перемещение которых будет рассчитываться на текущем шаге. Более подробно данная методика описана в разделе 3.9. Зеленым цветом выделены операции, которые реализуются в моделирующей программе с помощью параллельных вычислений. 3.6 Взаимосвязь локальных и сеточных величин
При моделировании методом частица-сетка распределения всех пространственных величин определяются в узлах двухмерной сетки. В то же время сами макрочастицы располагаются в произвольных точках. Это, с одной стороны, требует распределения вклада от каждой макрочастицы в величины плотности в узлах, а с другой стороны – определения сил, действующих на каждую частицу в зависимости от ее расположения.
Существуют различные методы интерполяции вклада макрочастицы в значения плотности в узлах [69], [73]. Наиболее простой метод заключается в передаче всего веса макрочастицы в ближайший к ней узел. Более сложные методы предусматривают вклад частицы в несколько окрестных узлов. При этом предполагается, что облако реальных частиц, которое олицетворяет макрочастица, имеет некоторую форму. Использование более сложных форм позволяет добиться большей гладкости получаемых распределений, однако это требует и больших вычислительных ресурсов.
В данной модели в качестве формфактора макрочастицы берется квадрат с шириной стороны, равной ширине ячейки пространственной сетки. Учитывая то, что моделирование осуществляется в осесимметричной геометрии, форма макрочастицы представляет собой кольцо квадратного сечения. На Рис. 9 показано представление макрочастицы в области моделирования.
На этом рисунке красная прямая – ось симметрии. Черные точки – узлы сетки. Узлы ниже оси симметрии показаны для наглядности, они являются отражением реальных узлов. Красной точкой представлена макрочастица (и ее отражение от оси симметрии), а зеленым – ее форма. Синими прямыми ограничены окрестности каждого узла сетки. В процессе раздачи плотности от частицы узлам каждый узел получает ту часть, которая приходится на его окрестность. Фиолетовым цветом выделены те узлы, вклад в которые внесет частица.
Определение концентрации каждой из компонент плазмы осуществляется в два этапа. Сперва в узлах сетки определяется суммарный вклад веса от всех макрочастиц. Затем по известному объему окрестности каждого узла определяется концентрация. На Рис. 10 показана схема раздачи веса макрочастицы в узлы.
Тестирование алгоритмов взаимодействия частиц с границами области моделирования
Алгоритмы сеточной интерполяции обеспечивают связь между локальными и сеточными величинами. То есть они, с одной стороны, осуществляют распределение вклада каждой частицы в значения плотности, определенные в узлах сетки, а с другой -определяют величину электромагнитной силы, приложенной к каждой частице, исходя из величин полей, определенных в узлах. Ошибки в данных алгоритмах могут привести к неправильному подсчету распределений концентраций компонент плазмы или к неправильному движению частиц. И то и другое нарушит самосогласованный расчет динамики плазмы.
Тестирование данных алгоритмов осуществлялось следующим образом. В область расчета загружалась равномерная невозмущенная плазма с заданной плотностью, после этого проводились подсчет концентрации в узлах сетки и сравнение с исходным значением. Загрузка плазмы выполнялась весьма грубым способом. Исходя из объема области Va, заданной плотности п 1 и размера частицы s находилось число частиц, которые требовалось расположить в области моделирования:
После этого для каждой из новых частиц выбиралось случайно положение внутри области моделирования. Аксиальная и радиальная координаты определялись с помощью псевдослучайных чисел г1 иг2: z = 2maxЛ г = r Jr max v 2 В том случае, если частица из-за особенностей геометрии не попала в область моделирования, выбираются новые координаты. Полученное в результате распределение плотности в узлах должно быть равномерным и соответствовать заданной величине.
На Рис. 19 приведено распределение плотности вдоль радиальной координаты после загрузки однородной плазмы плотностью 1,01018 м-3 на сетку с шагом 0,5 мм при размере частицы 1,0108. Яркой красной линией показано осредненное значение по десяти ячейкам.
Отклонение величины плотности в узлах от заданного значения связано с размером макрочастиц и их количеством. При измельчении частиц в пределе отклонение стремится к нулю. Нарастание отклонения по мере приближения к оси связано с особенностями осесимметричной геометрии: в ячейках меньшего объема оказывается меньше частиц.
По результатам тестирования можно заключить, что алгоритмы сеточной интерполяции работают правильно. Дополнительно это проверялось в процессе тестирования движения частиц, поскольку величины электромагнитных полей всегда определяются в узлах.
Алгоритмы расчета электрического поля выполняют решение уравнения Пуассона с использованием распределений концентраций заряженных частиц в узлах сетки и граничных условий. Ошибки в данных алгоритмах могут привести к неправильному расчету электромагнитных сил, действующих на частицы, и нарушению самосогласованного движения компонент плазмы. Тестирование этих алгоритмов осуществлялось в два этапа. На первом этапе решалось уравнение Пуассона с нулевым пространственным зарядом. Это позволило проверить правильность учета граничных условий. В области моделирования была задана область, представляющая собой цилиндр. На основаниях этого цилиндра были заданы граничные условия с фиксированной величиной электрического потенциала. Результат тестирования показан на Рис. 20. Рис. 20 Расчет распределения электрического потенциала в цилиндрической области без пространственного заряда Полученное в результате линейное распределение соответствует аналитическому решению и подтверждает правильность работы алгоритмов.
Для тестирования работы алгоритма решения уравнения Пуассона при наличии пространственного заряда в области применялась задача расчета распределения потенциала внутри шара с однородным распределением объемного заряда и фиксированным потенциалом на поверхности.
Для удобства сравнения результата работы алгоритма с точным аналитическим решением условия тестовой задачи были выбраны следующим образом. Радиус шара равен 1 м. Плотность объемного заряда внутри шара вє0Кл/м\ Потенциал электрического поля в каждой точке поверхности шара задавался равным нулю.
Поскольку моделирующая программа позволяет задавать область моделирования в виде набора замкнутых отрезков, неизбежное присутствие углов на границе моделирования вносило свою погрешность в решение. Однако, как показали вычисления, задание контура полукруглой области моделирования двенадцатью отрезками достаточно для получения решения приемлемой точности. На Рис. 21 приведен результат работы программы для данной тестовой задачи.