Научная статья на тему 'Оптимизация проектно-баллистических параметров космического аппарата с солнечной тепловой двигательной установкой при формировании и коррекции низких околоземных орбит'

Оптимизация проектно-баллистических параметров космического аппарата с солнечной тепловой двигательной установкой при формировании и коррекции низких околоземных орбит Текст научной статьи по специальности «Механика и машиностроение»

CC BY
264
70
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
КОСМИЧЕСКИЙ АППАРАТ / ТЕПЛОВОЙ АККУМУЛЯТОР / ПРОЕКТНАЯ МОДЕЛЬ / МЕТОД УСРЕДНЕНИЯ / ОПТИМАЛЬНОЕ УПРАВЛЕНИЕ / SPACECRAFT / THERMAL STORAGE / DESIGN MODEL / AVERAGING METHOD / OPTIMAL CONTROL

Аннотация научной статьи по механике и машиностроению, автор научной работы — Храмов Андрей Александрович

Рассматривается задача оптимизации проектно-баллистических параметров космического аппарата с солнечной тепловой двигательной установкой при межорбитальных манёврах в нецентральном гравитационном поле Земли. Используется модель космического аппарата с нерегулируемым двигателем ограниченной мощности с аккумулятором энергии. Общая задача разделяется на динамическую и параметрическую части. С использованием усреднённой модели движения получен алгоритм определения приближённо-оптимального управления. Оптимизация проектных параметров сводится к итерационной процедуре. Показана эффективность использования двигательной установки для формирования и коррекции низких околоземных орбит.

i Надоели баннеры? Вы всегда можете отключить рекламу.

Похожие темы научных работ по механике и машиностроению , автор научной работы — Храмов Андрей Александрович

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

OPTIMIZATION DESIGN AND BALLISTIC PARAMETERS OF SPACECRAFT WITH SOLAR THERMAL PROPULSION SYSTEM FOR THE FORMATION AND CORRECTION LOW-EARTH ORBIT

The problem of optimizing the design and ballistic parameters of the spacecraft with solar thermal propulsion system for interorbital maneuvers in non-central gravitational field of the Earth is considered. A model of a spacecraft with a fixed engine of limited power with power storage is used. The overall objective is divided into dynamic and parametric parts. Using the averaged model of the motion detection algorithm obtained approximately-optimal control. Optimization of design parameters is reduced to an iterative procedure. The efficiency of the propulsion system for the formation and correction of low-Earth orbit is shown.

Текст научной работы на тему «Оптимизация проектно-баллистических параметров космического аппарата с солнечной тепловой двигательной установкой при формировании и коррекции низких околоземных орбит»

УДК 629.78

ОПТИМИЗАЦИЯ ПРОЕКТНО-БАЛЛИСТИЧЕСКИХ ПАРАМЕТРОВ КОСМИЧЕСКОГО АППАРАТА С СОЛНЕЧНОЙ ТЕПЛОВОЙ ДВИГАТЕЛЬНОЙ УСТАНОВКОЙ ПРИ ФОРМИРОВАНИИ И КОРРЕКЦИИ НИЗКИХ ОКОЛОЗЕМНЫХ ОРБИТ

© 2013 А.А. Храмов

Самарский государственный аэрокосмический университет имени академика С.П. Королёва (национальный исследовательский университет)

Поступила в редакцию 10.07.2013

Рассматривается задача оптимизации проектно-баллистических параметров космического аппарата с солнечной тепловой двигательной установкой при межорбитальных манёврах в нецентральном гравитационном поле Земли. Используется модель космического аппарата с нерегулируемым двигателем ограниченной мощности с аккумулятором энергии. Общая задача разделяется на динамическую и параметрическую части. С использованием усреднённой модели движения получен алгоритм определения приближённо-оптимального управления. Оптимизация проектных параметров сводится к итерационной процедуре. Показана эффективность использования двигательной установки для формирования и коррекции низких околоземных орбит.

Космический аппарат, тепловой аккумулятор, проектная модель, метод усреднения, оптимальное управление.

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

В настоящее время существуют разработки перспективных двигательных систем с накоплением энергии, которые по уровню тяги (порядка десятков Ньютон) и удельного импульса занимают промежуточное положение между ЭРД и жидкостными ракетными двигателями (ЖРД). В исследовательском центре имени Келдыша разработана солнечная тепловая двигательная установка (СТДУ) с подогревом рабочего тела, принцип работы которой заключается в следующем. В течение освещённого участка орбитального полёта с помощью солнечных батарей (СБ) происходит преобразование световой энергии в

Храмов Андрей Александрович, инженер-программист кафедры космического машиностроения. E-mail: hramovaa76@rambler.ru

тепловую, которая накапливается в тепловом аккумуляторе (ТА). Перед подачей рабочего тела (водород) в двигатель осуществляется его подогрев в ТА, что повышает удельный импульс тяги двигателя. СТДУ может работать в режиме жидкостного ракетного двигателя (ЖРД) на холодных компонентах топлива (водород-кислород). Характерной особенностью рассматриваемой двигательной системы является ограниченное время работы за одно включение, определяемое рабочей ёмкостью ТА, а также достаточно высокий уровень тяги, что позволяет сократить время манёвра при некотором уменьшении массы полезной нагрузки по сравнению с ЭРД.

В состав СТДУ могут входить следующие основные системы [1]:

- система электроснабжения КА, включающая солнечные батареи, аппаратуру преобразования и регулирования;

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

- система хранения и подачи компонентов топлива.

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

1 ^ пн

заданную орбиту при заданных граничных условиях перелёта, начальной массе КА т0 и времени перелёта Т. Для решения задачи будем использовать проектную модель КА с нерегулируемым двигателем ограниченной мощности с аккумулятором энергии [2]. Предполагается, что ДУ работает в двух режимах: двигатель включён, тяга р и ско-

рость истечения рабочего тела с постоянны, происходит разрядка теплового аккумулятора с постоянной мощностью N; двигатель выключен -параметры р и с равны нулю, производится зарядка аккумулятора от солнечных батарей с постоянной мощностью NS. Для анализа вводится массовая модель КА, представляемая в виде суммы масс отдельных компонентов аппарата:

т0 = т + т + тЕ + тЬ + то + т + т , (1)

О пн к Е Ь о рт спх ' V /

где тпн - масса полезной нагрузки, тк - масса элементов конструкции, системы управления и пр., тЕ - масса ТА, тЬ - масса СБ, шд - масса двигателя, ш - масса рабочего тела, ш - масса

рт г спх

системы подачи и хранения рабочего тела.

Массы отдельных компонентов КА принимаются в виде функций:

тк = Мкт0 , тЕ = УеЕ0 , т = УSNS ,

m

P

= у P, m =— T , m = ß

<д ' рт м ' спх f с

c

m

спх рт '

(2)

где Е0 - максимальная накапливаемая аккумулятором энергия (рабочая ёмкость аккумулятора), NS - мощность солнечных батарей, р - тяга двигательной установки, с - эффективная скорость истечения, Тм - моторное время, Мк, УЕ , УЬ , У0у , РСпх - соответствующие удельные массовые характеристики. Так как при проектировании КА удельные характеристики аппарата, как правило, известны, то принимаем их постоянными заданными величинами.

Для массы полезной нагрузки можно записать:

m

= m0 (1 -uk )~yee0 ~ysns -

-УдP -(1 + ßCnx )m0 1 - exp

Vx_

c

(3)

u

opt

(t )= arg min

К (u, x0, XK , T )

T = fixe, x о — fixe, x к — fixe

Параметрическая часть задачи формализуется следующим образом:

mnn = arg max mnH (( Vx = Vxmin ( x0 , xk , T^

Для решения динамической задачи используется модель движения в оскулирующих элементах [3]. Вследствие ограниченного времени работы ДУ (до 20 мин.) и относительно небольшой тяги (до 100 Н) приращения элементов орбиты за виток незначительны, а траектория перелёта является многовитковой. Это позволяет упростить систему уравнений движения на основе метода усреднения.

Используя оптимальное управление на витке, полученное в работе [3], модель движения на длительных интервалах времени с учётом второй зональной гармоники в разложении геопотенциала Земли представляется в виде:

dA „ A3 w0 г „ чП

— = 2--[- (п-а)\

dt \ и п

dq „ ¡A w0 . f ai fa, . .

— = 4--Lsin\ g + — Icos\ — |cos(f) -bk,

dt \ и п ^ 2) ^ 2

dk „ ¡A w0 . f ai ,

— = 4--sinl g + — I cosl — | sin(f) + bq,

dt V и п ^ 2) ^ 2

(4)

V dt

= l1 -f 1W0-

Здесь А - большая полуось орбиты, q и & -компоненты вектора Лапласа, £ - половина ширины разгонного участка с аргументом широты его центра а - длительность пассивного участка по аргументу широты, Ух - характеристическая скорость перелёта, М - гравитацион-

7 В В А + Лк

ный параметр, Ь =-«-, Л = —-—,

В состав вектора оптимизируемых проектных параметров входят тяга ДУ, мощность СБ и рабочая ёмкость ТА: у = (Р, NS, Е0). Если зафиксировать проектные параметры, то максимум полезной нагрузки тпн соответствует минимуму затрат характеристической скорости Ух. При этом общая вариационная задача разделяется на динамическую и параметрическую части.

Динамическая задача заключается в поиске оптимального управления, доставляющего минимум затратам характеристической скорости при ограничениях на векторы управления и и состояния х:

B =g(5c°sIi -1) , , = 2,6333-1010

2т] М

наклонение орбиты.

Для поиска оптимальной по быстродействию программы управления воспользуемся принципом максимума Понтрягина. Гамильтониан для системы (4) имеет вид:

H = 4 A — \ U п

. , _ п - a ,

ava\ i +

aa + sin(c + — )cos — X 2 2

X

t

cos f] + wk sin ^

> +

+ b(qVk -kfq)-1.

ср

Здесь У л, У1 ,ук - сопряжённые множители, уравнения для которых запишутся в виде:

H = 4.

Aw0_ U п

(

АУа £

п-а

2

+

. а а

+ sin(£ +— )cos—х 22

х( cos) + у/к sin)

+

b{qvk - k¥q)-1-

> +

(5)

Законы изменения параметров управления определяются из необходимого условия максимума гамильтониана:

dH л -= 4

d£ V

Awo. /и п

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

АУ А +

+ cosí £ + аа 1 cosí аа I х

х

((

cos r + ук sin

in))

= o,

dH

-= 4

dr V U п

A w0 . í а sin I £ +

а

cosí — I х

2 ) I 2

х (- у sinr + ук cosr)= 0.

После преобразований с учётом ограничения на величину активных участков на витке ( 0 < £<Ж — а) длительность разгонного участка равна:

~ п-а

£opt = —^+ arcsin

a¥a

V2 2 \ + У.

q 0

при

A\¥a

\ (i+g 2r)

< cos

а

£„, = (п.

а

) = const

при

a\\a

\ (i+tg ))

> cos

а

(6)

В общем случае программа изменения £{t) может включать три этапа: первый содержит один разгонный участок на витке: £ = Ж — X = const, второй - разгонный и тормозной участки: 0 < £ < Ж — X = var, третий -один тормозной участок: £ = 0.

Из условия максимума для параметра ^ получаем:

ш

opt

Уq

= tg)o + bt), )opt =)o + bt. (7)

Таким образом, оптимальное расположение центра активных участков относительно линии узлов Т] линейно зависит от времени.

Системы (4) и (5) с полученным управлением (6) и (7) при заданных граничных условиях образуют краевую задачу. Неизвестными являются начальные значения сопряжённых множителей Ул0, У10, Ук0 и время перелёта Г. В качестве невязок используются отклонения конечных значений фазовых параметров от требуемых и значение гамильтониана Н, которое должно быть равно нулю на оптимальной траектории.

Сформулированная краевая задача может быть решена приближённо. Для первого и третьего этапов уравнения движения (4) интегрируются аналитически:

л = Ло

(1 -?t )2

q = e0 cos(®0 + bt)-yin(i-pt)cos(r0 + bt), (8) к = e0 sin(®0 + bt) - у in (l - (pt)sin(r0 + bt).

Здесь e0 = Jk0 + q0 , co0 = arctg —— экс-

q o

центриситет и аргумент перигея начальной орбиты соответственно, р = ±w0 [i — — | ,

V м У я J

sin —

у = ±2-, знак «плюс» соответствует пер-

я — а

вому этапу, знак «минус» - третьему.

Для описания движения КА на втором этапе необходимо получить в явном виде зависимость параметра £{t). Система (4) и сопряжённая с ней система (5) допускает интеграл вида:

d (A уА) = H + 1 - b (q У к - кУ q)

= const

Ж 2

Так как для задачи на быстродействие справедливо тождество: Н = 0 , то

a\a =

Ao\ao -

1 - b(qоУкo - \o) 2

Г~2 Г х . л/у^0 + у 10С08 у

При этом закон изменения длительности разгонного участка записывается в виде:

. п — а . / \

ь=—2— + агс§1П —^21), (9)

2

t

где ( =

AVa

г~2 г а

V +V cos"

(2 =

1 ~ b{qoVkо - kvqo)

„ Г~2 Т а 2vVko +Vqo cos у

При полученной программе управления параметром £ аналитическое решение можно получить только для изменения большой полуоси орбиты:

л

А =

(

A w f {{)

2

1 + 2

"" ¡ л (р2 где f {t)= {( - (21)arcsin {( - (21)-

( arcsin (

j + Vi - t(j -(21)2 -V 1 -(j2.

,(Ю)

Приближённый закон изменения компонентов вектора Лапласа можно получить, полагая, что длительность активных участков на витке мала и большая полуось орбиты меняется незначительно, и можно принять: A(t) » Аср = const. Учитывая (9), запишем:

sin \£-

л - а

J = sin y = sin {( - (2t). Пр

и ма-

лых значениях ( y < 50 ): sin y ~ y = { )•

Тогда искомое решение имеет вид:

q = e0 cos(ffl„ + bt) - — x

2b

((2 I t sin {r0 + bt )-

sin r¡0 sin bt b

((

2b

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

- (t2 sin {r0 + bt)

1 -

4b

2

( y

t cos {r0 + bt )-

cos r0 sin bt

b

k = e0 sin(®0 + bt) + — x

(1(2 ^- t costr0 + bt) +

x

cosr0 sin bt b

(( 4b2>

2b

1 -■

(2

t sin {r0 + bt )-

V

sinr0 sin bt b

2

+ 12 costr0 + bt)

, (11)

где — = 4

A ср W0

л л

а

cos — .

2

С использованием полученных решений (10) и (11) приближённо-оптимальное управление и траектория перелёта определяются посредством следующего алгоритма:

1) задаются начальные значения коэффициента р?, начальное значение аргумента широты линии переключения и время манёвра Т;

2) с учётом Н = 0 вычисляется коэффициент (р1 : ср2 = 2 IА Л(1 + р2 );

у / л

3) вычисляются начальное и конечное значения длительности разгонного участка:

л - а

„ л - а

е0 = —z—+(u

4 =

+ {( -(2 T );

22

4) с учётом ограничения на длительность активных участков на витке ( 0 < £ < л - а ) определяется время каждого из этапов:

t =

^0 n

kn

(2

где п = 1..3 - номер этапа, £кп, £0п - конечное и начальное значения продолжительности разгонного участка для п-го этапа;

5) по соотношениям (10), (11) последовательно определяются значения фазовых параметров в конце каждого этапа Акп, qkn, ккп, при этом конечное значение параметра для предыдущего этапа является начальным для последующего;

6) определяются отклонения полученных конечных значений параметров орбиты от требуемых, по которым определяется следующее приближение коэффициента р?, аргумент широты линии переключения и времени манёвра Т.

Затраты характеристической скорости определятся как: V = I 1--. Полученный ал' I л)

горитм может использоваться при незначительных изменениях большой полуоси орбиты А и малых (менее 200) длительностях активных участков на витке.

Результаты могут использоваться для решения краевой задачи в исходной постановке, где начальное приближение сопряжённых множителей определяются как:

V = (1

A° %0&Г - k0 )sign{cosr0)

/ A0

2(2 +-

■>¡1 + tg2r0 cosa/2

2

x

+

AVa

Ф

iV^+Ä)

. I а —sign\ cos—

а I 2.

cos— 2

¥к о =¥ч&щ.

При наличии пассивных участков на витке полученное управление обеспечивает ограниченную область существования оптимальных перелётов. Границу области можно определить по управлению, состоящему из одного активного участка (разгонного или тормозного) на витке максимальной продолжительности:

л-а

sin а

>

И.

(12)

где И =

ln(Ak / A)

Vе2 - 2eke0CosАт,

■ + ео

а 0

s= 1 + - A0) (л-а),

Vi =V0+bt, V2 =V0+л + b(( +1).

(13)

Здесь Г]0 = arctg

A~

- bT

аргумент ши-

V™* J

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

T 1

мент времени. T = — Ф

(

1 -

время манев-

ра. - программы изменения расположения

центров активных участков до и после переключения, Ак = ек sin сок - е0 sm(сy0 + ЬТ), А~~ = ек cosсk - е0 cos(сy0 + ЬТ).

Момент переключения ¿1 равен:

ti =-

1-

exp

-ja ~2 + a ~2

V

. (14)

Затраты характеристической скорости в данном случае зависят только от требуемого изменения большой полуоси орбиты:

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

(

V. =

L A

Л

slgn(Ak - A0 ).

к J

Ащ = щ -Щ) - bT

- изменение аргумента перигея под действием тяги ДУ. При малых длительностях активных участков ( а ^ я ) условие (12) определяет множество пересекающихся начальных и конечных орбит. С увеличением активных участков область существования перелётов расширяется и захватывает часть множества непересекающихся орбит. В пределе, при управлении без пассивных участков ( а = 0 ) эта область включает в себя любые замкнутые орбиты.

Если условие (12) не выполняется, то оптимальное управление содержит максимальный разгонный или тормозной участок на витке с переключением аргумента широты центра его приложения на 1800 в определённый момент времени ¿1. При этом параметры управления определятся как:

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

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

1) Задаются тяга двигательной установки ( Pmin < P < Pmax ), начальная масса КА m0, граничные условия x0, xk и время перелёта T. Минимальная величина тяги P ограничивается

min

временем перелёта без выключения ДУ; максимальная P = 100 Н - техническими ограни-

max 1

чениями.

2) Определяется оптимальное управление и соответствующая траектория перелёта. Рассчитываются затраты характеристической скорости Vx, максимальное значение большой полуоси орбиты Amax и длительность пассивного участка на витке а .

3) Из условия эффективного использования ТА, предполагающего полный его разряд после активного участка максимальной длительности, рассчитывается рабочая ёмкость теплового аккумулятора:

E0 = At

0 а

x Ne =(л-а).

í3

max

Pc

I Vt

максимальное время работы ДУ; NE

где At

^ a max

мощность, отводимая от теплового аккумулятора к двигателю.

4) Мощность СБ рассчитывается из условия полного заряда ТА от СБ после цикла работы двигательной установки, включающего в общем случае два активных и два пассивных участка:

Ns =

At,.

T б -At

об а

-Ne =

л-а Pc а 2vt

где Тоб - период обращения, - тяговый КПД.

5) Рассчитывается масса полезной нагрузки по соотношению (3).

На рис. 1 приведены результаты параметрической оптимизации при формировании рабочей орбиты для исходных данных, представленных в табл. 1, 2. Формирование рабочей ор-

Таблица 1. Параметры КА с СТДУ

М Ув ,кг/Мдж У5 , кг/кВт Уд, кг /Н 5 1 спх Лт,% Ш0 ,кг Р Н жрд с, м/с

0,1 1,18 50 (СТДУ-1) 20 (СТДУ-2) 0,39 0,5 80 7700 234 4532 (режим 1) 7456 (режим 2)

Таблица 2. Параметры орбиты выведения и рабочей орбиты

Орбита выведения Рабочая орбита

А0, км е0 а>0, град А, , км к 5 ек а,, град

6578 0,00205 0 7858 0,00248 0

-СТДУ-1 -В-СТДУ-2 -О-ЖРД

20 25

Время, сут.

Рис. 1. Зависимость массы полезной нагрузки от времени манёвра

биты КА осуществляется в два этапа: переход с орбиты выведения на круговую орбиту высотой 300 км на холодных компонентах топлива в режиме ЖРД (режим 1) и многовитковый переход на рабочую орбиту в режиме без дожигания водорода, нагретого в тепловом аккумуляторе (режим 2).

Результаты параметрической оптимизации при коррекции орбиты в режиме работы СТДУ без дожигания водорода приведены на рис. 2 для исходных данных, представленных в табл. 1, 3. Начальная масса КА принималась равной 6570 кг, наклонение орбит при коррек-

ции 69,9 град.

Результаты расчётов показывают, что при определённом увеличении времени межорбитального манёвра СТДУ может обеспечить выведение полезной нагрузки большей массы по сравнению с кислород-водородным ЖРД. Эффективность рассматриваемой двигательной установки в значительной мере определяется энергомассовым совершенством системы электроснабжения. Приведённая методика определения проектных параметров может использоваться на начальных этапах проектирования КА с СТДУ.

Таблица 3. Параметры начальной и конечной орбиты при коррекции

Начальная орбита Конечная орбита

А0, км е0 ю0, град А, , км ек а,, град

6750 0,0149 120 6896 0,0239 0

о ЖРД -&-СТДУ-1 -В- СТДУ-2

I 5500

о 5400

15 20

Время, сут.

Рис. 2. Зависимость массы полезной нагрузки от времени манёвра

2.

СПИСОК ЛИТЕРАТУРЫ

Солнечный тепловой ракетный двигатель [Электронный ресурс]: Патент Российской Федерации. URL: http://ru-patent.info/21/25-29/2126493.html (дата обращения 29.05.2013).

Механика космического полета: Проблемы оптимизации / Г.Л. Гродзовский, Ю.Н. Иванов, В.В. Токарев.

М. : Наука, 1975. 704 с.

3. Храмов А.А. Оптимальные программы коррекции слабоэллиптических и круговых орбит космических аппаратов с двигателем ограниченной тяги // Вестник Самарского государственного аэрокосмического университета имени академика С. П. Королева (национального исследовательского университета). Самара, 2011. Вып. 2. С. 112 - 122.

5100

OPTIMIZATION DESIGN AND BALLISTIC PARAMETERS OF SPACECRAFT WITH SOLAR THERMAL PROPULSION SYSTEM FOR THE FORMATION AND CORRECTION LOW-EARTH ORBIT

© 2013 A.A. Khramov

Samara State Aerospace University named after Academician S.P. Korolyov (National Research University)

The problem of optimizing the design and ballistic parameters of the spacecraft with solar thermal propulsion system for interorbital maneuvers in non-central gravitational field of the Earth is considered. A model of a spacecraft with a fixed engine of limited power with power storage is used. The overall objective is divided into dynamic and parametric parts. Using the averaged model of the motion detection algorithm obtained approximately-optimal control. Optimization of design parameters is reduced to an iterative procedure. The efficiency of the propulsion system for the formation and correction of low-Earth orbit is shown. Spacecraft, thermal storage, design model, averaging method, optimal control.

Andrey Khramov, Software Engineer at the Space Engineering Department. E-mail: hramovaa76@rambler.ru

i Надоели баннеры? Вы всегда можете отключить рекламу.