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

-►

РАДИОФИЗИКА

УДК 537.86

А.Ю. Рыкшин, П.М. Зацепин, П.Н. Уланов

Алтайский государственный университет

моделирование рассеяния короткого

электромагнитного импульса на двухслойной структуре при помощи вейвлет-преобразования

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

ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЕ, ИМПУЛЬС, РАССЕЯНИЕ, ДВУХСЛОЙНАЯ СТРУКТУРА.

Введение

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

В связи с этим значительный интерес представляют задачи распространения и рассеяния коротких и сверхкоротких электромагнитных импульсов в различных структурах [1 — 3].

При теоретическом исследовании задач распространения и рассеяния коротких и сверхкоротких электромагнитных импульсов возникают трудности, связанные с выбором метода решения таких задач. В данном случае главное требование к методу — это возможность учета импульсного характера источника первичного поля. Традиционные методы решения, в основе которых лежат преобразования Фурье или Лапласа, во многих случаях непригодны, так как приводят либо к трудностям в ин-

терпретации полученного решения, либо к сложности численных расчетов. Решение этих проблем возможно, если использовать подход, который бы учитывал локальный характер первичного поля. Указанным требованиям может удовлетворять метод, в котором вместо преобразования Фурье по временной координате используется интегральное вейвлет-преобразование.

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

Постановка задачи

Бесконечно протяженный вдоль оси у линейный источник с плотностью электрического поверхностного тока / имеет координаты (х0, £0) в декартовой системе. Источник излучает первичный импульс, имеющий горизонтальную поляризацию, и его поле рассеивается на двухслойной структуре с толщинами слоев ¿2 и ¿у Диэлектрические проницаемости верхнего и нижнего полупространств е1, е4, а также слоев исследуемой структуры е2, е3 могут

Рис. 1. Геометрия задачи:

й?2, й?3 — толщины слоев структуры; — относительные диэлектрические проницаемости сред; (х0, х0) — координаты источника излучения

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

Геометрия задачи представлена на рис. 1.

Источник тока является импульсным с временной зависимостью вида:

j с ) = - ^

exp

(t - Q2

2т2

(1)

где т — размерный коэффициент, задающий полуширину импульса.

Поле источника падает на диэлектрический слой, частично отражается и проходит внутрь этого слоя; там оно испытывает множественные переотражения и затем выходит из слоя (см. рис. 1) [4].

Для решения поставленной задачи нами применяется преобразование Фурье по пространственным координатам и вейвлет-преобразование по временной координате. Вейвлет-преобразование имеет следующий вид:

Ж (О, Т) = | / () |П|1/2 - Т ))А;

1

/а) = — Г -Ж(О, Т) |о| 1/2 - Т))с!ШТ; С -

|| -ад

здесь С^ — нормировочная константа, которая рассчитывается следующим образом:

У(ю) 12

с. = =

œ

-dœ,

(у(ю) — Фурье-образ вейвлетной функции

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

\|/ПТ а) = -оа - Т ) ехр (' 1)

Электромагнитные поля и граничные условия

В силу двумерности задачи система уравнений Максвелла распадается на две независимых подсистемы с компонентами (Еу, Нх, Ни (Ну, Ех, Е) — горизонтальная и вертикальная поляризации соответственно.

Далее, рассмотрим решение для случая горизонтальной поляризации. Система уравнений Максвелла для выбранной поляризации примет следующий вид:

dE„

dz dE„

= hr

dHr

dt

dH

dx

= -hr

dt

(2)

dH dH

dE

d dx

Z = so8 —TT + j

dt

где jey — плотность тока импульса (ток проходит по нити источника).

Для вывода компонент в разных областях пространства воспользуемся техникой, описанной в работах [5, 6]. Представление компонент электромагнитного поля через вейвлет-образы по временным координатам и Фурье-образы по пространственным координатам в полупространстве z > 0 записывается следующим образом:

1

Ey ( x, z, t ) =

С. (2п)2

x exp[/^x + inz]

1/2

Hx ( x, z, t ) =

1

(2n)

j y, П, T) x

—да

. QT (t )d £,d n d QdT ; j x, n, n,T) x

(3)

I 11/2

х ехр[/Е,х + ¿пг ] О .пТ (^ )Т ^Т пТ ОТТ; 1

Н (X, г, t) =

С. (2п)2

I к; & п, О, Г) х (3)

|1/2

х ехр[/Е,х + ¿пг] О .пТ(t)ТпТОТТ.

Для полного описания решения доста- где t =

точно найти компоненты (Е , Нх).

Подставляя выражения для компонент полей (3) в систему уравнений Максвелла (2), получим следующие соотношения для спектральных плотностей компонент падающего поля:

еУ&п,П,Т) = -ЦоIе&п,П,Т) х

.1ОТ (t)

п2 + ^2 + 8оМ"08.ШТ (t)

.ОТ (t)

ых & п, п, Т) = ^

.От (t) _

/ / п.ОТ 0 )

^0.1ОТ () '

(4)

где

Г & п, О, Т) = I у- (х', г', t >

-¡Чх '-щг' х

х ^|П|.пТ (t ')Тх' й;' Т', /(х', г', t') = 8(х'- Х0)8(;'- го)Г(');

.1пг (t) = ■

(5)

Выражения для компонент электромагнитного поля, с учетом применения теоремы о вычетах при интегрировании по координате п и подстановки спектральной плотности тока (5), примут следующий вид:

-I +да

Еу(х,;, t) = I --^ехрК(х - ±

У 4пС. -1 Ж

± щ(г - го Ж(О, Т) |п|1/2 .ют(tОТТ;

1 +да

Нх (х, г, t) = —— I ± ехр [/%(х - Х0) ± (6)

± ¿ж1(г - г0)]/(О, Т) |о| 1/2 .пТ ^)Т^ТОТТ;

(

Ч2 -

80 ^0 8п .ЮТ (t)

.От ( )

1/2

V ТОТГ/ у

Для упрощения записи введем оператор запаздывания по времени Ъ(t0) : /(t) ^ /(t - t0). Тогда граничное

условие на границе диэлектриков г = 0 запишется в следующем виде:

7Т • "У ;

|Ъ (^0)Еу + Еу = ЕТ [Ъ (tо)нX + ну = НГ,

(7)

2 У/2 Х0 + г0

ч с 2Ке 81,

Представляя прошедшее и отраженное поля излучения аналогично компонентам падающего поля (6) и используя систему (7), получим для отраженного и прошедшего полей на границе г = 0 следующие выражения:

Еу (х, г, t) = Ъ (tо) х

4пС.

+да

; | е1уЯи ехр[/^х + щг]ТОТТ;

-да

1 ъ,. /г Ж1.пТ(t)

(8а)

Щ (х, г,t) =

4пС,

■Ъ(tо) I

^0.1ОТ (t) (8б)

е'уЯ12 ехр[/+ /ж1г]ТОТТ;

1 +да

У(х, г, t) = -— Ъ(0)| еу/ х

4ПС. -да

х ехр[/£,х - /ж2г]Т\сТОТТ;

(8в)

НГ (х, г, t) =

1

4пС,

■Ъ (tо) I

'№2.ОТ 0) х

М'о .1ОТ (0 (8г)

еуТ12 ехр[/Е,х - /ж2г]ТОТТ,

где

еу = е;(£,, О, Г, t) = -^°ехр[-/- /^1 го]

ж

/ (о, г) |о| 1/2.1от (*);

м> -м> 2ж

п = п т гг = п

пт1 ' пт

+ м> т

п т

тп + м> т

п т

Используя тот же способ, получим с помощью формул (8в), (8г) выражения для напряженностей прошедшей и отраженной компонент поля, образованных в результате падения излучения на границу г = — Т2:

1 +да

(х, г, t) = (tо + t2) I еу

Я23Т12еЪж2йг ехр[/^х + /ж2г]ТОТТ;

жп =

Е1 (х, г, t) =

1

+с

Ъ(tn + Г в'т23т,2

(п 2)-| ^ 23 12 (9)

х ехр[г'Я2^2 - ^3)] ехр[/^х - тъг]ЯОЯТ,

где tn = Яп /сп — время задержки при прохождении среды с номером п, сп — скорость света в среде с номером п.

В результате падения излучения с компонентой напряженности электрического поля Еу из области первого слоя структуры (см. формулы (9)) на границу между слоями г = П компоненты поля излучения принимают следующий вид:

Еу(х, г, ^) = -—=-Ъ(^ + Ъ2) Г вЦ х

4пЧ -с

X в^^ * --2гй^ШТ-

л +с

Ету (х, г, ^) = —— Щп + Ъ2) { вУу^ Т12 х

4пЧ -с

х exp[2iw2d2]T21 ехр[/£,х + iw1г]d ОЯТ.

Рассмотрим переотражения от границ в среде 2 и в результате получим конечное выражение для компонент полей излучения в этой среде:

Е2+ (х, г, t) =

1

4пСх

( ( Ъ

V V

+?+■»( ^ ( Я + -Л

I" Е П , . ^ . , Я2 + г

tn + (2п + 1)^ + ■

вуТ12С23 Х

2 У

ехр)[/Е+х + г>2 г]Я О ОЯТ;

Е2 (х, г, t) =

1

(1П)

4пС..

(^ ( Ъ

V V

tn + 2Ш2 - ■

е'Т Яп

12Л1

2

х ехр^х - iw2г]ЯОЯТ;

Я1 = С^Сд ехр[2/^2Я2].

Верхние плюсовый и минусовый индексы при величине Е2 означают движение поля излучения по направлению оси г и в противоположном направлении соответственно (нижний индекс указывает на номер среды). Результирующее поле Е2( х, г, /) во второй среде (первый слой структуры) движется против оси г и попадает на границу г = —Я2, где частично отражается и частич-

но проходит в третью среду. Это прошедшее поле создает в третьей среде (второй слой структуры) множественные внутренние переотражения:

Е+ (х, г, t) =

1

4тсС

{ЕЕ (Ъ (О:

Ъ ([2 п + Щ + [2 т + Щ) Ъ

( Я2 + Я3 + г ^

х вуТ12Я23Т23Я34ЯпЩ ехр[2/^2Я2] х ехр[/'^2 - w3)d2] ехр[2/^3(Я2 + Я3)]) х

х ехр[/'^х + iw3 г ]Я ОЯТ; (11)

Е- (х, г, t) = -1— {ЕЕ (Ъ (tn) х

Ъ

[2п + 1К + 2mt3 -

Я2 + г

вУТ12Т23Я1пЯТ

3 У

х ехр[/^2 - w3)d2]) ехр[/'^х - iw3г]ЯОЯТ; С = Я,2Я34 ехр[2/^3Я3].

Далее поле Е3- (х, г, ?) падает на границу сред 3 и 4 в плоскости г = —Я3 и частично проходит в среду 4:

1 +с

Е-, (х, г, 0 = ^ {ЕЕ (Ъ (О:

-с0 е=П и=П

4кС„

х Д[2п + 1]t2 + [2т + 1]^) х

х Ъ ^ +Я3 + г 1 в2'' х (12)

V С4 Уу

х ехр^'^ - w3)d2 + i(w3 - w-)(d2 + Я3)]) х х ехр^х - iw-г]ЯОЯТ.

Поле Е3+ (х, г, t) в свою очередь падает на границу сред 3 и 2 в плоскости г = — и частично проходит в среду 2:

+с +с +с +с

Е2+ (х, г, t) = ^ {ЕЦ (Ъ (О х

4тсС

1=п т=П п=п

(

х Ъ

^ -с

[2п + 1]t2 + [2т + 2]^ +

+ г

л

2У

Ъ (21^ ')вуТ12С23Т23С34Т32С\+^ С2т

(13а)

х ехр^^^] ехр^г'^ - w3)d2] х х ехр[2г>3(Я2 + Я3)]) ехр^х + iw2 г]Я ОЯТ;

Е2 (х, г, ^) =

1

+с +с +с

4тсС

{ЕЕЕ(П(tn)х (13б) 139

1=п т= п п=п

х Ъ

[2п + 1]! + [2т + 2К +

Т2 + г ^

С2

\п+1 тут

х Ъ(.и^е'уТцУ23Т23У,4Т32У21У1П+'Ут х х ехр[2/ж2Т2] ехр[2/(ж2 - ж3)Т2] х х ехр[2/ж3(Т2 + Т3)]) ехр[/;х + /ж2 г ]Т ОТТ.

Поле Е2+ (х, г, t) (см. формулы (13а)) падает на границу сред 2 и 1, в результате чего появляется компонента Е :

1 +да

(х, г, t) = —— I е1, О, Г, г, ^):

4пС. -да

х ехр[/;х + /ж1г ]Т ОТТ;

(16а)

1 +да

Е- (х, г, t) = —— I е~— О, Г, г, t):

4пС. -да

х ехр[/£,х - /ж1г]Т\ТТОТТ;

Е1+ (х, г, t) =

1

+да +да +да

4пС. -да

ЦИ (Ъ (О:

/=0 т=0 п=0

I)([2п + Щ + [2т + 2]^ + [2/ + 1]t2):

х Ъ

' г ^

V с1 у

(14)

e■yГ12R13Г23R34Г32Г21R1 У х

х ехр[2/ж2Т2] ехр[2/(ж2 - ж3)Т2] х х ехр[2/ж3(Т2 + Т3)]) ехр[/;х + /ж2 г]Т ОТТ.

Затем движение поля излучения повторяется: Е2 (х, г, t) проходит в среду 3, испытывает переотражения, проходит в среду 4, а также возвращается в среду 2, где вновь испытывает переотражения, проходит в среду 1 и падает на границу сред 2 и 3, а затем все вновь повторяется.

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

Е -2 = Ъ (2п<2)Г12Я;;

Е +2 = Ъ([2п + 1]t2)Т12R23R1n ехр[2Ж2Т2];

Е32 = Ъ(t3 + 2п<2)Г32Л1п ехр[/(ж2 - ж3)Т2];

Е -2 = Ъ (<2)В21 Е+2;

Е+1 = Ъ (t2)Т21 Е+2;

Е23 = Ъ(t2 + 2п<3 )Т23У2п ехр[/(ж2 - ж3)Т2];

Е+3 = Ъ(<3)Я,4 ехр[2/ж3(Т2 + Т3)]Е-3;

Е-4 = Ъ(t3)Т34 ехр[/(ж3 - ж4)(Т2 + Т3)]Е-3;

М = Е32 Е23, X = Е32 Е23.

Тогда выражения для напряженностей электрического поля примут следующий вид:

(15)

1 +да

(х, г, t) = I е2+, п, Г, г, t):

. -да

х ехр[...]Т ОТТ;

1 +да

(х, г, t) = —— I е--, О, Г, г, t):

4пС. --да

х ехр[/;х - /ж2г]Т\ТТОТТ;

(16б)

2

+да

1 +да

(х, г, t) = —— I е3+(;, п, Г, г, t) х 4пС -1

. -да

х ехр[/;х + ^г ]Т ОТТ;

1 +да

Е3 (х, г, t) = —— I п, Г, г, t) х

--да (16в)

х ехр[/;х - /ж3г]ТОТТ;

1 +да

Е4- (х, г, t) = —— I е4, О, Г, г, t) х

4пС. --да

х ехр[/;х - /ж4г]ТОТТ,

где спектральные плотности полей будут иметь следующий вид:

еу;

е- (;, О,Т, г, t) = Ъ

е+ (;, п,г, г, t) = Ъ

+да п

адо;, п, Г, t) + X Е 31Е 23 ГЕ 12еу

' г ^

V с1 У

^ г >

V с1 У

(17а)

е2 (;, О,Т, г, t) = Ъ

( - \ +да

V с2 У

е2+ (;, О,Т, г, t) = Ъ

X X Е 12еу;

п=0

Т2 + г

п=0

^ Т2 ^^

<о+ 2

(17б)

2 У

Е +2еу +Х МХПВ -2еу

е32 (;, О, Т, г, t) = Ъ

<о-

тТ2 + г ^

-3 у

(17в)

х

Е П-3 -2ву (О, Т, t);

в3+ & О, Т, г, t) = Ъ

( Я2 + Я3 + г ^

t +—2-3-

»п т

Е П23 П 12ву ;

(17в)

е4- & О, Т, г, t) = Ъ

( + Я3 + г ^

t -г'2 1 "3

'п

1-4 у

Е П24 П 12ву .

Результаты численных расчетов

Рассмотрим некоторые численные результаты, полученные с использованием приведенных выше выражений для на-пряженностей электрической компоненты поля в верхнем полупространстве (среда 1). На рис. 2 приведены три графика для на-пряженностей полей, нормированных на максимум, во временной развертке. Они построены для относительных диэлектрических проницаемостей е14 = 1, е3 = 6, е2 = {2; 1п; 8п} и толщины слоев = = Я3 = п,1 м. Заметно смещение относительных максимумов амплитуды напряжен-ностей для различных комбинаций диэлектрических проницаемостей.

Рис. 2. Временная развертка нормированных

напряженностей электрического поля при изменяющейся величине е2: 2 (кривая 1), 1п (2), 8п (3); = 1, е3 = 6, е4 = 1; толщина слоев структуры Я2 = = п,1 м

£>Л£>\, 1,0 -0,8 -0.6 -0.4 -0,2 -0,0 --0.2 --0.4 --0.6 --0,В --1,0 -

-10 -5 0 5 10 1 5 20 25 30 35

{/ - Гц/с,)/1,

Рис. 3. Временная развертка нормированных

напряженностей электрического поля при изменяющейся толщине слоя структуры Я2, м: п,1 (1), п,2 (2), п,3 (3); йъ = п,1 м;

= е4 = 1, е2 = 4, е3 = 8п

На рис. 3 приведены еще три графика для нормированных на максимум на-пряженностей электрического поля во временной развертке. Они построены для одинаковых значений относительных диэлектрических проницаемостей е14 = 1, е2 = 4, е3 = 8п, толщины второго слоя структуры Я3 = п,1 м и различных значений толщины ее первого слоя Яг Как и на рис. 2, заметно смещение относительных максимумов напряженности поля для различных толщин слоев.

Заключение

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

Работа выполнена при поддержке РФФИ конкурс «р_сибирь_а» проект № 13-п5-98п37.

п=п

х

п=п

х

п=п

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

1. Рыкшин А.Ю., Зацепин П.М. Моделирование рассеяния электромагнитного импульса на диэлектрическом слое // Научно-технические ведомости СПбГПУ. Физико-математические науки. 2009. № 4(88). С. 14-18.

2. Zatsepin P.M., Ulanov P.N., Rykshin A.Yu. Wavelet-transform usage for solving the problem of short electromagnetic pulse scattering by the perfect conductive strip // International Siberian Conference on Control and Communications Proceedings. Krasnoyarsk, 2013.

3. Ulanov P.N., Zatsepin P.M. Applying of the wavelet transform to solve the problem of coaxial waveguide pulsed excitation // International Siberian Conference on Control and Communications Proceedings. Krasnoyarsk, 2013.

4. Бреховских Л.М., Годин О.А. Акустика слоистых сред. М.: Наука. Гл. ред. физ.-мат. лит-ры, 1989. 416 с.

5. Rykshin A., Zatsepin P., Komarov S. Scattering of electro-magnetic pulse on the impedance strip into the dielectric layer // XIIIth International Conference on Mathematical Methods in Electromagnetic Theory (MMET'2010) Proceedings. 2010. Kiev, Ukraine. 2010.

6. Zatsepin P.M., Rykshin A.Y. Electromagnetic pulse scattering by impedance cylinder // Proceedings of 2008 IEEE Region 8 International Conference on Computational Technologies in Electrical and Electronics Engineering, Sibircon 2008. Novosibirsk, Russia. 2008. Pp. 350-354.

СВЕДЕНИЯ ОБ АВТОРАХ

РыКШИН Алексей Юрьевич — кандидат физико-математических наук, доцент кафедры вычислительной техники и электроники Алтайского государственного университета. 656049, Россия, г. Барнаул, пр. Ленина, 61 а1ехеу-г@уа^ех.га

ЗАЦЕПИН Павел Михайлович — кандидат физико-математических наук, завидующий кафедрой вычислительной техники и электроники Алтайского государственного университета. 656049, Россия, г. Барнаул, пр. Ленина, 61 zpm@phys.asu.ru

УЛАНОВ Петр Николаевич — аспирант кафедры вычислительной техники и электроники Алтайского государственного университета.

656049, Россия, г. Барнаул, пр. Ленина, 61 u1anovpn@gmai1.com

Rykshin A.Yu., Zatsepin P.M., Ulanov P.N. THE MODELING OF A SHORT ELECTROMAGNETIC PULSE SCATTERING BY DUAL LAYERS STRUCTURE USING WAVELET TRANSFORM.

In the paper the problem of a short electromagnetic pulse scattering by dual dielectric layers structure has been considered. An incident field was generated by infinitely thread-like source of electric or magnetic current located parallel to the boundaries of layers. The solution was written through the Fourier transform with respect to the spatial coordinates and wavelet transform with respect to the time coordinate. Electromagnetic fields characteristics within the given structure and in the surrounding space were studied.

WAVELET TRANSFORM, PULSE, SCATTERING, TWO-LAYER STRUCTURE.

REFERENCES

1. Rykshin A.Yu., Zatsepin P.M. Modelirovanie rasseyaniya elektromagnitnogo impulsa na dielektricheskom sloe// St. Petersburg State Polytechnical University Journal : Physics and mathematics, 2009, No. 4(88), pp. 14-18. (rus)

2. Zatsepin P.M., Ulanov P.N., Rykshin A.Yu. Wavelet-transform usage for solving the problem of short electromagnetic pulse scattering by the perfect

conductive strip // 2013 International Siberian Conference on Control and Communications Proceedings. Krasnoyarsk, 2013.

3. Ulanov P.N., Zatsepin P.M. Applying of the wavelet transform to solve the problem of coaxial waveguide pulsed excitation // International Siberian Conference on Control and Communications Proceedings. Krasnoyarsk, 2013.

4. Brekhovskikh L.M., Godin O.A. Acoustics of layered media. Moscow: Nauka. Head Ed. Sci. Lit. 1989. 416 p. (rus)

5. Rykshin A., Zatsepin P., Komarov S.

Scattering of electro-magnetic pulse on the impedance strip into the dielectric layer // Xlllth International Conference on Mathematical Methods in Electromagnetic Theory

(MMET'2010) Proceedings. 2010. Kiev, Ukraine. 2010.

6. Zatsepin P.M., Rykshin A.Yu. Electromagnetic pulse scattering by impedance cylinder // Proceedings of 2008 IEEE Region 8 International Conference on Computational Technologies in Electrical and Electronics Engineering, Sibircon 2008. Novosibirsk, Russia. 2008. pp. 350-354.

THE AUTHORS

RYKSHIN Alexei Yu.

Altay State University,

61 Lenin St., Barnaul, 656049, Russia

alexey-r@yandex.ru

ZATSEPIN Pavel M.

Altay State University,

61 Lenin St., Barnaul, 656049, Russia

zpm@phys.asu.ru

ULANOV Petr N.

Altay State University,

61 Lenin St., Barnaul, 656049, Russia

ulanovpn@gmail.com

© Санкт-Петербургский государственный политехнический университет, 2014