Труды МАИ. Выпуск № 96
http://trudymai.ru/
УДК: 621.455.4,537.5
Оценка параметров прямоточного высокочастотного ионного
двигателя
Гордеев С.В.*, Канев С.В.**, Суворов М.О.***, Хартов С.А.****
Московский авиационный институт (национальный исследовательский университет), МАИ, Волоколамское шоссе, 4, Москва, A-80, ГСП-3,
125993, Россия *e-mail: svyatoslavgordeev@mail. ru
**e-mail: k208mai@mail.ru ***e-mail: maxwell2005@yandex.ru ****e-mail: skhartov@ya.ru
Аннотация
В работе исследована работа высокочастотного ионного двигателя, использующего в качестве рабочего тела смесь газов верхних слоёв атмосферы, который может применяться для компенсации аэродинамического торможения низкоорбитального космического аппарата. Проведено моделирование рабочих процессов в двигателе, проанализированы основные его характеристики и с учетом геометрии космического аппарата показана принципиальная возможность существования аппаратов с прямоточным электроракетным двигателем.
Ключевые слова: низкоорбитальный космический аппарат, аэродинамическое торможение, высокочастотный ионный двигатель.
Введение
В последнее время возрастает интерес к использованию низкоорбитальных космических аппаратов (НОКА) различного назначения. Их достоинствами является низкая цена выведения на орбиту (высоты 160-300 км) и повышенное разрешение целевой аппаратуры. Однако при снижении высоты орбиты космических аппаратов (КА) возрастает концентрация газов остаточной атмосферы Земли, и связанное с этим аэродинамическое торможение аппарата, что приводит к уменьшению срока его активного существования. Для компенсации этого торможения на высотах 200300 км возможно использовать электроракетные двигатели (ЭРД). Один из примеров успешного использования ЭРД на низкоорбитальных аппаратах-спутникООСБ Европейского космического агентства, эксплуатировавшийся более четырех лет на высотах 250-280 км. Аэродинамическое торможение этого аппарата компенсировалось ионным двигателем Т-5 фирмы QinetiQ [1]. КА GOCE выполнял научные задачи по изучению гравитационного поля Земли и океанских течений до тех пор, пока на борту не закончился запас рабочего тела для двигателя -ксенона. Для перспективных низкоорбитальных аппаратов представляет интерес возможность использования в качестве рабочего тела остаточных газов верхних слоев атмосферы, так как в этом случае снимается одно из ограничений для срока активного существования НОКА - запас рабочего тела двигателя.
В мире ведутся работы по созданию прямоточных электроракетных двигателей[2,3]. При разработке таких двигателей необходимо учитывать, что миделевое сечение устройства забора газов определяет не только расход рабочего
тела через двигатель (в конечном итоге его тягу), но также влияет и на торможение аппарата. Площадь солнечных батарей, необходимых для питания целевой аппаратуры, определяет и мощность, которую возможно подвести к двигателю, и также влияет на торможение НОКА[4].
В работе приведена методика расчета и оценены параметры прямоточного ЭРД, использующего в качестве рабочего тела газы остаточной атмосферы, для НОКА, функционирующего на солнечно синхронной орбите высотой около 200 км. Рассматривалась схема ионного двигателя с высокочастотной ионизацией рабочего тела [4].
Методика оценки аэродинамического торможения НОКА на высоте
около 200 км
На высоте около 200 км концентрация окружающих КА частиц имеет
порядок 1016 ^ [1,4]. При этом режим обтекания поверхностей космического
аппарата - свободномолекулярный. Это означает, что соударения молекул друг с другом практически не происходит, и они взаимодействуют только с корпусом космического аппарата. Также известно, что частицы остаточной атмосферы имеют максвелловское распределение по скоростям. Существует модель расчета усилий (торможения), возникающих на поверхностях КА, для таких условий [5], которая использует следующие предположения.
При соударении молекулы с элементом выпуклой поверхности, на нём возникают нормальные и касательные напряжения. Эти напряжения будут определяться из следующих условий:
РТ1,ЛФ)=^АТ^Х1 , (0) , (2)
где: индексом 1 обозначаются функции, которые ищутся на элементе поверхности, обращенном к потоку (направление проекции внутренней нормали к поверхности на направление потока совпадает с этим направлением), индексом 2 обозначаются функции, которые ищутся для теневых по отношению к потоку элементов^ = 7800 м/с - скорость набегающего потока (орбитальная скорость КА на высоте 200 км),
k = 1.38-1 0 _ 2 - постоянная Больцмана, т - молекулярная масса газа, Т5-
температура набегающего потока;ф - угол между вектором скорости набегающего потока и внутренней нормалью к обтекаемой поверхности и - коэффициенты аккомодации нормальной и тангенциальной компоненты импульса молекул соответственно, характеризующие передачу количества движения от молекул
стенке; Тш- температура стенки;р - плотность среды, 55 = и/^- безразмерная
скорость потока, 5П ( ф )=55соБф - проекция безразмерной скорости потока на внутреннюю нормаль к поверхности. Х1 _ 2 ( ф ) -функционал, зависящий от 5П ( ф ) :
Х1 2 ( ф ) = ехр( -5П(ф ) 2) ± ( ф ) ( 1 ± ег^5п(ф) ) ) . (3)
Тогда сила, воздействующая на элемент поверхности dS, если угол между внутренней нормалью к элементу и направлением потока - острый, будет:
dFdl (ф) = dS(0) ((Рщ (ф) cos ф + Рт± (0) sin 0). (4)
А в случае, когда угол между внутренней нормалью к элементу и направлением потока тупой:
dFd2(0) = dS(ф)((Рп2(ф) cos ф + Рт2(ф) sin ф). (5)
Проинтегрировав выражения 4 и 5 по площади, можно найти силу аэродинамического торможения для выпуклого тела любой формы, для компенсации которой необходимо использовать ЭРД.
Оценка параметров высокочастотного ионного двигателя при работе на
атмосферных газах
Использовалась модифицированная модель высокочастотного (ВЧ) разряда в ионизационной камере ЭРД, работающего на атмосферных газах [6].
Для оценки концентраций частиц в ионизационной камере запишем уравнения баланса «прихода» и «ухода» частиц (появление ионов вследствие ионизации, нейтрализация ионов на стенках и вылет частиц через ускоряющую систему двигателя). Примем, что в камере присутствуют только атомарный и молекулярный кислород и азот, так как концентрация других видов частиц в атмосфере на высотах около 200 км существенно меньше [1,4]. При этом будем учитывать только процессы ионизации всех видов частиц и процессы диссоциации молекул. Для
точной оценки необходимо найти концентрацию электронов, а также ионов и нейтральных частиц четырех видов, всего девять неизвестных. Примем, что все ионы, появляющиеся в объеме ионизационной камеры, попадая на ее стенки, нейтрализуются.
Для атомарного кислорода это условие описывается соотношением [7]:
пеп0( о-т;0 ¿) Ккам = 0 . 7п ^ 0 5СТ, (6)
[кТ
где: ^0 = I —-- скорость выпадения ионов атомарного кислорода на стенку (так
МЮ
называемая «бомовская» скорость), Ге - температура электронов, 0 - масса иона кислорода, Ккам - объем ионизационной камеры, 5СТ- площадь стенок, ( бги 0 ¿)-коэффициент скорости ионизации, пе - концентрация электронов и п;0 -концентрация ионов атомарного кислорода в камере, - концентрация нейтральных атомов кислорода в камере.
Аналогичные уравнения можно написать для ионов молекулярного кислорода, а также для атомарного и молекулярного азота.
Концентрацию электронов можно посчитать из условия квазинейтральности плазмы[8]:
П е = П; 0 + П I 0 2 + П ш + 2. (7)
Для заданной конфигурации устройства забора газов известен расход кислорода и азота через двигатель. Соотношения атомарной и молекулярной компонент газов различны на входе и выходе из камеры.
Для кислорода можно записать:
™0_02 = ПЮУЮМЮСТ^5И0С + Щ02УЮ2^Ю2(Т^И0С + + Що^юО- ~
- а) 0 ^ 0 + (По2^ + Щ02^ 0 2 ( 1 - а))^5И 0 ^ 0 2, (8)
где: - суммарный массовый расход атомов и молекул кислорода, -
прозрачность ионно-оптической системы (ИОС) для ионов, а0 - прозрачность ИОС
~ /8 кТу, с
для нейтральных атомов, V 0 = I — — тепловая скорость атомов кислорода, 5И 0 с -
площадь перфорированной части ИОС.
Если сгруппировать слагаемые в правой части, с учетом того, что
п я о п я ^Ю VО
МIо 2 = 2 М; о, V Iо 2 = -щ и Vо 2 = ^, можно получить:
^ 2 + ^ = - (/2"о2 + По) ^^ (9)
Аналогичное уравнение можно написать и для азота.
Рассмотрим уравнение баланса для молекул кислорода в предположении, что на входе в камеру все частицы находятся в состоянии молекул (это справедливо, если используется камера-термолизатор скоростей частиц [9]):
пеп02((ау02+ {(ТУ021)~)Укам - 0.1щ02уЮ2Б^ +
+ (По 2^ + П ;о 2V ;о 2 ( 1 - а) с = ~ (10)
I ¿-(То Мх)2
Если записать аналогичное уравнение для молекул азота, то получим систему из девяти уравнений, решив которую, найдем концентрации всех рассматриваемых видов частиц.
Рассчитаем другие параметры ЭРД с высокочастотным разрядом. Нас будет интересовать потребляемая мощность и тяга двигателя. Потребляемая двигателем
мощность складывается из затрат на ионизацию и ускорение рабочего тела [10]. Тяга двигателя определяется плотностью тока заряженных частиц, извлекаемых из плазмы высокочастотного разряда.
Затраты на ионизацию при известных концентрациях частиц определяются кольцевыми токами в высокочастотной плазме, для расчета которых необходимо знать вихревое электрическое поле, появляющееся в расчетной области вследствие колебаний тока в индукторе двигателя. Зная, что ток в индукторе изменяется по синусоидальному закону, будем искать только амплитудные значения интересующих нас величин.
Амплитудную напряженность поля можно найти, зная векторный магнитный потенциал в каждой точке пространства. Для этого используют уравнение Максвелла [11]:
(11)
Из определения магнитного потенциала:
гоЫ = Б.
(12)
Зная, что:
Б = Б0 5 1 п ( о .
(13)
Примем для амплитудных значений:
Вп = го1:Л
о
0 .
(14)
Тогда с учетом (11):
^ _ ^ . д зт(й)С)
го ЬЕ = -го1:Л0—
И амплитудное значение напряженности поля:
£о = - о4 о . (16)
Таким образом, чтобы найти напряженность поля вокруг индуктора необходимо знать потенциал магнитного поля и частоту токов в нем.
Рассчитаем векторный потенциал магнитного поля индуктора. Магнитный потенциал, создаваемый током I, протекающим по элементу dL (Рис. 1) равен:
сС4 = — ■ —, (17)
471 г ' 4 '
где: ц0 = 4л: ■ 1 0 _ 7 —— магнитная постоянная, ц — магнитная проницаемость среды (в вакууме ц = 1).
С1А
Рис. 1 - Схема расчета векторного магнитного потенциала элементарного
участка проводника с током. Найдем магнитный потенциал в произвольной точке плоскости, содержащей кольцевой виток с током (Рис.2). Магнитный потенциал ищется в точке, находящейся на расстоянии a от центра витка, при известном радиусе витка R и токе I протекающем в нем.
Рис. 2- Схема расчета векторного магнитного потенциала в плоскости кольцевого
витка с током
В этом случае расстояние г от элемента dL до точки, в которой ищется потенциал, будет равно:
г = / а 2 + Д2 — 2 • а • Я - с о 5 р. (18)
Тогда векторный магнитный потенциал запишется в виде:
4=цац^2» !■ с о зф а (19)
4л 0 // а2+К 2- 2-а-Я- с о эф
Если необходимо найти магнитный потенциал витка в точке, лежащей на расстоянии Ь от его плоскости (рис.3) то расстояние от элемента dL до этой точки, будет:
г = / а 2 + Ь2 + Д 2 — 2 • а • Я • с о5 р. (20)
и магнитный потенциал:
4 = Мц^2л 1-с° зф а р . (21)
4л // а 2+ Ь2+ К 2- 2^ а^ с о Бф
пространства вблизи витка с током.
На базе этих соотношений можно рассчитать токи, текущие в плазме. В высокочастотном разряде токи в плазме обусловлены в основном токами электронов, так как их подвижность существенно выше, чем у ионов, поэтому необходимо рассмотреть уравнение движения электронов в плазме. На электрон воздействует вихревое синусоидальное электрическое поле, и он теряет свой импульс в процессе столкновения с атомами и ионами плазмы. Таким образом, уравнение движения электронов можно записать в виде:
7тге1>е = — 5 то t — те1;е V , (22)
где: 7?ге = 9 . 1 • 1 0- 3 1 кг - масса электрона, >е - скорость электрона, V - частота соударений электрона с тяжелыми частицами.
Это уравнение можно решить аналитически:
еЕп
>р =-Г^-7Т ( О С О 5 О { — V 5 1 П О . (23)
е те(й)2+г72) 4 у 47
Если данное выражение умножить и разделить на Vо2 + V 2, можно получить:
еЕ0л/ш2+у2 ( со V . Л /ООЛ
>е =-Т^-ТП —-т===51ПО ) . (28)
CO V
Введя обозначения: , = с о s а, и , = s i п а, после упрощений
va)¿+v¿ va)¿+v¿
получим:
в En
ve = —■==== с о s ( о t + а) . (24)
е meVü)2+v2 v J v 7
Как видно, амплитуда колебаний скорости будет:
v =_eE° = (25)
Vemax (25) Частота соударений так же является функцией от скорости электронов:
У = л/ + <*.„ (26)
где: о е у - сечение рассеяния электронов на атомах ]-го сорта ( может принимать значения О, 02, №N2), ппу - концентрация нейтральных частиц ]-го сорта.
Данное соотношение учитывает только рассеяние электронов на нейтральных атомах, так как их концентрация в ВЧ разряде примерно на порядок больше чем ионов.
Плотность токов в плазме определяется соотношением [12]:
] п е^^еш ах. (27)
Зная концентрации электронов и атомов в плазме можно рассчитать токи, протекающие в ней.
Для расчета процессов в камере ионизации двигателя разобьем исследуемую область (рассматривается осесимметричная постановка задачи) на п элементов по горизонтали и по вертикали. Рассчитаем координаты г и 2 в узлах полученной сетки, и получим, таким образом, две матрицы размерностью п на п. Далее рассмотрим
амплитудные токи в каждом узле сетки, вызываемые магнитным полем токов во всех элементах, где есть плазма и тока в индукторе. Обозначим:
к = — п ее—т==, (28)
где: - площадь сечения элементарной кольцевой ячейки, в которой рассчитывается ток. Тогда в каждой точке:
/ , у = ^4 ¿,у. (29)
В свою очередь, магнитный потенциал зависит от координат точки, в которой ищется, а также радиуса и расположения витка, в котором течет ток, порождающий его. То есть магнитный потенциал в точке кольца с координатами (г,у, г; у), порождаемый витком с током I, с координатами (Я,^) будет:
4 = / (Г ,7,/) . (30)
При этом магнитный потенциал прямо пропорционален току, который его порождает, и можно записать:
4 = / (г,у,^у,Д ,71 ) /. (31)
Для магнитного потенциала, который ищется в точке с координатами (г у, г;,у), и создается кольцевым током с координатами ( ), введем обозначение:
4 = / ( г ,у, г, у, г , ъ г , ъ 1 ) • /, к = /¿к/ , (32)
Магнитный потенциал в точке ( г,у, г^у ), порождаемый индуктором, обозначим:
4 инд¿ у = /инд( г,у,г^ ,уД инд,7инд , 1 ) • / инд. (33)
Тогда можно записать:
/о , О — к , О + /о^о1^ - 1"'" ■ ■ ■ + /оП0П^п,п + А инд О , о) 5
V 1 — к (/одо/о, о + /одЧ) , 1+■ ■ ■ +/()дп/п,п + А индо, О 5 (34)
^п,п — к (/п,п /о, о + /п,п /о, 1+ ■ ■ ■ + /п,п ^п,п + А индп,п) 5
Разделив обе части на к, перенеся слагаемые —, вправо, а Аинд^ у влево
получим систему из п2 уравнений:
С/о, О Ц) ^ , о + /(), о ^ , 1+^ ■ ■ +/о , о — А индо ,о>
/о, 1 + С/о, 1 т) 1 + ■ ■ ■ +Л), 1 — А индо, 15 (35)
/п , п ^ , о + /п , п ^ , 1+ ■ ■ ■ + С/п,п ~ ^ — А индп,п-
Так как частота соударений является функцией скорости, данная система не является линейной. Решим её сначала при 17ет ах — 0, а затем сделаем несколько итераций по нелинейности (решение сходится за 2-3 итерации).
Решив данную систему, найдем токи в каждой элементарной ячейке плазмы. С их учетом можно пересчитать реальный магнитный потенциал в каждом элементе5 расчетной области.
Тогда электрическая мощность, передающаяся от электронов рабочему телу,будет [13]:
^Р 1 — (36)
Электрическая мощность, затрачиваемая на ускорение ионов равна [14]:
Руск = 0 . 7 аХос е //+ £у Пу>у, (37)
где: //+ - потенциал эмиссионного электрода двигателя, п ¿у - концентрация ионов ]-й компоненты (]может принимать значения О, 02, Ки К2), > ¿у-«бомовская» скорость ]-й компоненты.
Тягу двигателя можно найти, зная плотность «бомовского» тока через отверстия ИОС [15]. С учетом наличия в двигателе смеси частиц тяга определяется соотношением:
^ = 0 . 7 О^иос £ у пу> ¿у М ¿уД^, (38)
где: М ¿ у - масса иона ]-й компоненты.
С использованием данных соотношений можно оценить потребляемую мощность и тягу ЭРД.
Пример расчета для полета НОКА по солнечно синхронной орбите
В качестве примера рассмотрим аппарат следующей формы (Рис.4) -цилиндрический корпус диаметром 200-300 мм, по краям которого размещены панели солнечных батарей (СБ) шириной 500-600мм.Длина аппарата варьируется в пределах 0,5-4 м для получения зависимости мощности СБ от величины силы аэродинамического сопротивления КА.
явление полета
СБ —-
Корпус КА с двигателем и целевой аппаратурой
Рис.4 - Форма КА, используемая в расчетах В основном сила сопротивления полету и, тяга ЭРД (расход через двигатель) для ее преодоления будет определяться миделевым сечением корпуса КА. Однако вытянутая форма корпуса аппарата может также влиять на параметры системы аппарат-двигатель - при изменении длины КА будет изменяться мощность, генерируемая солнечными батареями, и также при этом будет возрастать действующая сила торможения.
Прямоточный ЭРД в расчетах имел конфигурацию, представленную на рисунке 5. В данной схеме применено устройство забора газов (УЗАГ ), имеющее входную часть в виде сотовой конструкции [15], которое существенно сокращает длину двигателя. Для уменьшения скоростей частиц между УЗА и камерой ионизации установлена камера термолизатора в форме трубы диаметром 100 мм. Ионизация и ускорение частиц осуществляется с использованием схемы
высокочастотного ионного двигателя с диаметром ИОС 150 мм. Целевая аппаратура НОКА размещена над двигателем в пределах габаритов УЗАГ.
Набегающий поток
УЗАГ
Нейтрализатор
Индуктор Ионизационная камера
ИОС
Рис.5 - Конфигурация ЭРД, принятая в расчетах. Мощность, вырабатываемая солнечными батареями, для солнечно синхронной орбиты КА высотой 200 км равна:
Nb = Г Nyd с о s <р dS ,
^ СБ
(39)
где: Ыуд - мощность получаемая с 1 м солнечных батарей на орбите Земли при
ВтN
стационарном нормальном их освещении (в расчетах принято Nyd = 3 00 —), ср -
м"
угол между внутренней нормалью к поверхности солнечных батарей и вектором
направления от Солнца, S^ солнечными батареями.
площадь поверхности аппарата, покрытой
Результаты оценок параметров ЭРД представлены на рисунке 6 в виде зависимости мощности, вырабатываемой солнечными батареями, от силы торможения, и зависимость мощности, потребляемой двигателем, от силы тяги. В
качестве варьируемого параметра для управления тягой двигателя, выбирался
потенциал эмиссионного электрода. При этом следует помнить, что плотность частиц в атмосфере существенно изменяется на протяжении срока активного существования аппарата. Существуют суточные и сезонные колебания плотности частиц, а кроме того она сильно зависит от солнечной и геомагнитной активности [15]. Изменение плотности частиц влечет за собой изменение силы торможения КА и изменение расхода через двигатель. На рисунке 6 представлен вариант расчета для средней солнечной и геомагнитной активности. Аналогичные расчеты для минимальной и максимальной солнечной и геомагнитной активности представлены на рисунках 7 и 8 соответственно. Как видим, для данной конфигурации КА с прямоточным ЭРД для средней солнечной активности при равенстве тяги двигателя и силы торможения КА5 мощность, вырабатываемая солнечными батареями выше мощности, потребляемой двигателем, то есть возможно стационарное существование аппарата с прямоточным электрореактивным двигателем.
н
рр
А
1.5 х10
2 х10
1 х10
500
0
3
3
3
0
5х10
- 3
0.01
Сила тяги, Н Сила торможения, Н
Рис.6 - Результаты оценок параметров двигателя для средней солнечной
активности для радиуса цилиндрической части КА 150 мм, ширины батареи 600 мм. Длина корпуса - варьируется, при этом изменяется величина торможения и
При минимальной солнечной активности (Рис. 7) расход через двигатель уменьшается настолько, что тяги не хватает для компенсации торможения, даже при достаточно высоких потенциалах эмиссионного электрода ИОС. Следует, однако, заметить, что в этом случае КА будет снижаться под воздействием аэродинамического торможения, концентрация частиц возрастет и реализуется режим работы двигателя аналогичный полету при средней солнечной активности.
мощность солнечных батарей.
н РР
ч
Г)
л н о о к 3
о
(и
ьч ч ю (и о н о
Е
н РР
и
л н о о к 3
о
2.4 Х10
1.8x10
3
1.2 х103 600 0
0
2x10
-3
- 3
4 Х10 3
Сила тяги, Н Сила торможения, Н
Рис. 7 - Результаты оценок параметров двигателя для минимальной солнечной
активности
т
РР
Ч Р
п
ь т с о н 3 о
е я
л ю
е р
т о
Е
т
РР
С
ь т с о н 3 о
3
2.1x10"
3
1.7х10Т
3
1.3Х10Т
900"
500
5x10
- 3
0.01
0.015
Сила тяги, Н Сила торможения, Н
Рис. 8 - Результаты оценок параметров двигателя для максимальной солнечной
активности
При максимальной солнечной активности (Рис.8) практически вся мощность, вырабатываемая солнечными батареями, потребляется двигателем. При этом высота орбиты сохраняется, так как торможение скомпенсировано тягой ЭРД.
Выводы
Была разработана методика расчета параметров прямоточного высокочастотного ионного двигателя, учитывающая геометрические характеристики космического аппарата. С помощью разработанной методики продемонстрирована возможность существования космического аппарата на солнечно-синхронной орбите высотой около 200 км с ЭРД для поддержания орбиты на базе прямоточного высокочастотного ионного двигателя, работающего на остаточных газах верхних слоев атмосферы. При этом, в отличии от известных работ [4, 6, 9], в данной методике более полно учтены геометрические характеристики аппарата, влияющие на его торможение и на выработку электрической энергии.
Разработанная методика расчета может применяться для анализа НОКА различных форм и компоновок прямоточного электроракетного двигателя.
Библиографический список
1. Corbett, M. H., Edwards, C. H. Thrust Control Algorithms for the GOCE Ion Propulsion Assembly // 30th International Electric Propulsion Conference. Florence, Italy, September 17-20, IEPC-2007-210.
2. Nishiyama, K., Air Breathing Ion Engine Concept // Proc. of the 54th Int. Astronautical Congress, Sept. 29-Oct. 3, 2003, Bremen, IAC-03-S4-02.
3. Hruby V., Pote B., Brogan T., et.al. Air Breathing Electrically Powered Hall Effect Thruster // Patent US 6,834,492 B2, Busek Company, Inc., Natick MA, USA, Dec. 2004.
4. Канев С.В., Петухов В.Г., Попов Г.А., Хартов С.А. Прямоточный электрореактивный двигатель для компенсации аэродинамического торможения низкоорбитальных космических аппаратов // Известия Вузов. Авиационная техника. 2015. № 3. С. 35-40.
5. Бондарев Е.Н., Дубасов В.Т., Рыжов Ю.А., и др. Аэрогидромеханика - М.: Машиностроение, 1993. - 608 с.
6. Канев С.В., Кожевников В.В., Хартов С.А. Физико-математическая модель процессов в ионизационной камере электроракетного двигателя, использующего атмосферные газы в качестве рабочего тела // Известия РАН. Энергетика. 2017. № 3. С. 21-30.
7. Белоусов А.П., Мельников А.В., Хартов С.А. Модель динамики электронов в разрядной камере высокочастотного ионного двигателя // Труды МАИ. 2017. № 94. URL: http://trudymai.ru/published.php?ID=80974
8. Godyak V.A., Piejak R.B., Alexandrovich B.M. Electron Energy Distribution Function Measurements and Plasma Parameters in Inductively Coupled Argon // Plasma Sources Science Technology. 2002. Vol. 11, pp. 525-543.
9. Ерофеев А.И., Суворов М.О., Никифоров А.П., Никифоров А.П., Сырин С.А., Попов Г.А., Хартов С.А. Разработка воздушного прямоточного электрореактивного двигателя для компенсации аэродинамического торможения низкоорбитальных космических аппаратов // Вестник НПО им. С.А. Лавочкина. 2016. № 3. С. 104-110.
10. Groh K.H. and Loeb H.W. State-of-the-Art of Radio-Frequency Ion Thrusters // Journal of Propulsion and Power. 1991. Vol. 7, no. 4, pp. 573-579.
11. Канев С.В., Латышев Л.А., Нигматзянов В.В., Хартов С.А. Моделирование рабочего процесса в газоразрядной камере высокоточного ионного двигателя // Труды МАИ. 2012. № 52. URL: https://www. mai.ru/science/trudy/published.php?ID=29483.
12. Кожевников В.В., Хартов С.А. Распределения локальных параметров плазмы в газоразрядной камере высокочастотного двигателя // XLI Академические чтения по космонавтике. Тезисы докладов. Москва, 24-27 января 2017. С. 76.
13. Александров А.Ф., Вавилин К.В., Кралькина Е.А., Неклюдова П.А., Павлов В.Б. Исследование параметров плазмы индуктивного ВЧ-источника плазмы диаметром 46 см. Параметры плазмы в области скин-слоя // Прикладная физика. 2013. № 5. С. 34-37.
14. Алдонин Ф.И., Ахметжанов Р. В. Расчет основных рабочих характеристик ионного двигателя мощностью 20 - 30 Вт. // Труды МАИ. 2015. № 81. URL: https://www.mai.ru/science/trudy/published.php?ID=57827
15. Абрагян В.К., Ахметжанов Р.В., Лёб Х.В., Обухов В.А., Черкасова М.В. Численное моделирование первичного пучка ионов и потока вторичных ионов в ионно-оптической системе ионного двигателя // Труды МАИ. 2013. № 71. URL: http: //trudymai .ru/published.php?ID=46702