Научная статья на тему 'Численные расчеты характеристик тлеющего разряда в разрядных камерах с поперечным потоком электроотрицательного газа'

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

CC BY
96
24
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
ТЛЕЮЩИЙ РАЗРЯД / ПОТОК ГАЗА / РАЗРЯДНАЯ КАМЕРА / ПРОСТРАНСТВЕННОЕ РАСПРЕДЕЛЕНИЕ ЗАРЯЖЕННЫХ ЧАСТИЦ / ЭЛЕКТРИЧЕСКИЙ ПОТЕНЦИАЛ / GLOW DISCHARGE / GAS FLOW / DISCHARGE CHAMBER / SPATIAL DISTRIBUTION OF CHARGED PARTICLES / ELECTRIC POTENTIAL

Аннотация научной статьи по физике, автор научной работы — Сафиуллин Р. К.

Тлеющий разряд в потоке газов (ТРП) используется для создания активных сред мощных молекулярных лазеров, которые все шире используются для обработки различных стройматериалов и изделий. Важной практической целью исследований ТРП является достижение оптимальных условий превращения электрической энергии разряда в энергию лазерного излучения. В данной работе численно исследуется ТРП в поперечном (относительно взаимного расположения электродов) потоке газа. Рассматривается случай сплошного анода и секционированного катода. Работа является продолжением исследований [1-3].

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

Numerical investigation of glow discharge characteristics in transverse discharge chambers in electronegative gas flow

Transverse glow discharge chambers in electronegative gas flow were investigated numerically by alternating direction method. A case of continuous anode and partitioned cathode was considered. As a result, two-dimensional distributions of charged particles and of electric potential inside such chambers were obtained. These investigations continue those in the papers [1-3]. A system of elliptic equations (4), (7-10) was solved by alternating direction method11 in two-dimensional subspace XY. Calculations were carried out for a rectangle with sides a and b which contained 1 mm width lamellate cathode and continuous anode. The uniform difference grid (100x100) was used.The boundary conditions were taken as in [7]. They were also widely varied with a view to evaluate their influence on the solution. Under given boundary conditions the program elaborated enables to calculate two-dimensional distributions of electric potential and charged particle densities inside discharge chamber. The calculations showed the strong dependence of these distributions on boundary conditions. They showes also that in contrary with the results of one-dimensional analysis in [6], E has a maximum near cathode and decreases towards anode. The same result was also received in [13]. It indicates that in view of strong nonhomogeneity of such GDF the one-dimensional approach evidently turns out to be inadequate. The calculations were carried out for gas velocities in the range 30-250 m/sec. The voltage drop between electrodes was put equal 1400 V, the static pressure p at entry was varied in the range 4-10 kPa, and gas temperature T at entry was put equal 293 K. The qualitative and semi-quantitative agreement with experiment [13] was achieved.The calculations showed also that space distributions of charged particles strongly depend on ambipolar diffusion coefficients Dea, Dpa, Dna. As was proposed in [4] ambipolar diffusion coefficient for negative ions Dna could be even negative. The difference scheme of alternating direction method gives diverging solution at negative values of Dna. As follows from our calculations, the deviation from plasma quasi-neutrality may achieve 3-4 % within the positive column, whereas near electrodes it may be an order of magnitude more (depending on boundary conditions). As an example, the calculated space distributions of electron and ion densities and also of electric potential inside a discharge chamber are presented in figs. 2-14.

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

УДК 537.525

Сафиуллин Р.К. - доктор физико-математических наук, профессор

E-mail: rksaf@mail.ru

Казанский государственный архитектурно-строительный университет

Адрес организации: 420043, Россия, г. Казань, ул. Зелёная, д. 1

Численные расчеты характеристик тлеющего разряда в разрядных камерах с поперечным потоком электроотрицательного газа Аннотация

Тлеющий разряд в потоке газов (ТРП) используется для создания активных сред мощных молекулярных лазеров, которые все шире используются для обработки различных стройматериалов и изделий. Важной практической целью исследований ТРП является достижение оптимальных условий превращения электрической энергии разряда в энергию лазерного излучения. В данной работе численно исследуется ТРП в поперечном (относительно взаимного расположения электродов) потоке газа. Рассматривается случай сплошного анода и секционированного катода. Работа является продолжением исследований [1-3].

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

пространственное распределение заряженных частиц, электрический потенциал.

1. Математическая модель зарядовой кинетики в газоразрядной плазме в электроотрицательном газе

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

Принятая система уравнений имеет вид:

divre = (a - aa)ne + adnn - репепр+ q, (1)

divrp = an - ЬеППр - рцПрПп + q, (2)

divr = aan - adnn - рцпрПп, (3)

Dj = - (np - ne - nn)e/e0, (4)

где

Ге = neU - ne^eE - De Vne, Гр = npu + nPmpE - Dp Vnp, Гп = nnU - nn^nE - Dn Vnn E = -Vj,

(5)

у = в(Гр - Ге - Гп) = е(пр/1Р+пе/1е+пп/1п)Е + в(пр-пе-пп}ы + в(Ое Упе + Оп Уп„- Ор Упр). (6)

Здесь Ге, Гр, Гп - векторы плотности потоков электронов, положительных и отрицательных ионов, соответственно; а , аа , оь - частоты ионизации молекул газа электронным ударом, прилипания электронов к электроотрицательным молекулам и отлипания электронов от отрицательных ионов; Де, Да - коэффициенты электрон-ионной и ион-ионной рекомбинации; q - интенсивность внешнего источника ионизации; пе, пр пп -концентрации электронов, положительных и отрицательных ионов; у - вектор плотности тока; Ве, Ор, Вп, /ле, цр> тп - коэффициенты диффузии и подвижности соответствующих заряженных частиц; и - вектор скорости потока газа в данной точке разрядной камеры (РК); Е - вектор напряженности электрического поля; е0 - диэлектрическая постоянная; V и А -операторы градиента и Лапласа, соответственно; е - абсолютная величина заряда электрона.

В условиях квазинейтральности плазмы тлеющего разряда первый член правой части выражения (6) является доминирующим, поэтому, с учетом того, что /1е » №п,

получаем j ~ enmE. Из уравнений (1-3) следует соотношение divj = 0. Иногда для расчета электрического потенциала это соотношение применяют к приближенному выражению j ~ enemeE . Однако, как будет видно из дальнейшего, более точно распределение электрического потенциала внутри РК можно рассчитать по уравнению Пуассона (4).

Уравнение Больцмана для ФРЭЭ может быть записано в виде:

Ie + Ie-M + In + (SUP) V + (SUP)eiectr = 0. (7)

Здесь Ie - член, описывающий нагрев электронов в электрическом поле; IeM - осредненный по углам интеграл упругих столкновений электронов с атомами и молекулами; Iin - интеграл неупругих столкновений электронов с молекулами, учитывающий возбуждение вращательных, колебательных и электронных уровней энергии; члены (SUP)V и (SUP)el£cr описывают удары второго рода с колебательно- и электронно-возбужденными молекулами, когда электроны получают от них энергию. Конкретный вид членов левой части уравнения (7) приведен в работе [4] с описанием метода его решения. Использование уравнения (7) совместно с уравнениями (1-4) целесообразно ввиду того, что частоты ионизации и прилипания электронов к молекулам о- и о сильно зависят от параметра E/N, где E = \E\, N - полная плотность числа атомов и молекул в газовой смеси. Использование же в зарядовой кинетике имеющихся в литературе аналитических аппроксимаций для коэффициентов о- и О может приводить к серьезным ошибкам ввиду того, что они могут более, чем на порядок, отличаться от данных эксперимента.

Нами была рассмотрена РК с поперечным ТРП, через которую газ прокачивался в направлении оси Х, а электроды располагались, как показано на рис. 1. Использовался сплошной анод и секционированный катод. Катодная плата представляла собой набор узких пластинчатых (ножевых) катодов, каждый из которых подключался к источнику питания через балластное сопротивление. Как известно, секционирование катода применяется для повышения устойчивости тлеющего разряда и увеличения энерговклада в разряд [5, 6].

Y

<-------------------- -Ь

Рис. 1. Схема РК. Межкатодное расстояние а = 4 см, Ь = 3 см, с = 8 см,

Ь = 32 см, длина катодов I = 3,6 см

При переходе в приближение амбиполярной диффузии, которое часто используется для описания процессов в газоразрядной плазме [5-8], уравнения (1-3) для рассматриваемой РК приобретают следующий вид:

ихдпе/дх + иудпе/ду = Оеа (д2пе/дх2 + д2пе/ду2) + (а - аа)пе + аап„ - ре,пепр+ д, (8)

ихдп/дх + иудп/ду = Ора (с^п/дх2 + с^п/ду2) +а пе - ЬеППр - РиПрпп + д, (9)

ихдпп/дх + иудпп/ду = Ош (д2пу/дх2 + д2п/ду2) + а^е - оьпп - Ьппрп№. (10)

Электрическое поле входит в уравнения (8-10) неявным образом через величины а аа, рег и це. Наиболее сильно от параметра Е/И зависит величина а;■. Здесь введены

следующие обозначения: Dea, Dpa, Dna - коэффициенты амбиполярной диффузии для электронов, положительных и отрицательных ионов, соответственно. Общие выражения для них приведены в [8, 9].

Аналогичным образом могут быть записаны уравнения для энергии (энтальпии или температуры) газа, а также уравнения для колебательных температур в случае СО2-лазеров или для населенностей колебательных уровней молекул СО и N2 в случае СО-лазеров.

Следует отметить, что подобная РК рассматривалась в работах [7, 10], в которых считалось, что плазма всюду квазинейтральна (np ~ ne + nn), т.е. рассматривался только положительный столб разряда. Более того, в этих работах использовалось точное условие нейтральности плазмы, так как одно из уравнений зарядовой кинетики, а именно, уравнение для концентрации положительных ионов np, из рассмотрения исключалось. В нашей же модели решаются все уравнения зарядовой кинетики, что позволяет оценить степень отклонения от квазинейтральности как в приэлектродных областях, так и в положительном столбе разряда.

В данной работе мы ограничились рассмотрением системы уравнений (4), (7-10). Течение газа считалось ламинарным. Профиль скорости газа на входе и внутри РК задавался пуазейлевским или полагалось для простоты ux = const, uy = 0.

2. Метод решения системы уравнений зарядовой кинетики и электрического поля

Система эллиптических уравнений (4, 8-10) решалась итерационным методом переменных направлений [11] в плоскости XY. Расчетная область представляла собой прямоугольник со стороной d = 8 см в направлении оси X и стороной b = 3 см, параллельной оси Y. Она содержала один или несколько пластинчатых катодов толщиной 1 мм, расположенных на верхней стороне прямоугольника, и часть сплошного анода (нижняя сторона прямоугольника). Эта область покрывалась равномерной разностной сеткой 100 х 100. Граничные условия для задачи выбирались такими же, как в работах [7, 10]. Помимо этого граничные условия варьировались в широких пределах с целью выяснить зависимость от них характера распределения концентраций заряженных частиц и электрического потенциала. Итерации применялись к каждому из уравнений (4, 8-10) в отдельности и к системе в целом, пока не достигалась относительная погрешность вычислений в пределах 1 % . Следует отметить, что для сравнения с результатами работы [10] распределение потенциала рассчитывалось нами не только согласно уравнению (4), но также из уравнения div(ne V j) = 0, примененного в [10].

3. Результаты расчетов распределения концентраций заряженных частиц и потенциала электрического поля внутри РК

Ниже приведены некоторые результаты ряда расчетов. Они получены при задании граничных условий Дирихле для искомых величин на всей границе рассматриваемой прямоугольной области, причем везде на границе использовалось условие квазинейтральности плазмы. Это означает, что рассматривался положительный столб разряда, который занимает практически весь объем камеры (толщина прикатодного и прианодного слоев в таких РК не превышает 1 мм [10,12]). На входе (х = 0) и выходе из расчетной области (х = d = 8 см) было принято: ne(y) = 106-108 см-3, np(y) = 1,1ne(y), nn(y) = 0,1ne(y), ф(у) = - 1400 y/b (от 0 на аноде до -1400 В при y = b). На аноде (y = 0) полагалось ne(x) = 1010 см-3, np(x)= 1,5.1010 см-3, nn(x) = 0,5.1010 см-3, ф(х) = 0. На катодной плате (y = b, рис. 1) было положено:

ne(x) = ne0(1 - п | sin[n(2xk/d - 1)] |), np(x) = np0(1 - n | sin[n(2xk/d - 1)] |), (11)

nn(x) = nn0(1 - n | sin[n(2xk/d - 1)] |), ф(х) = - 1400 + 500 | sin[n(2xk/d - 1)], (12)

где к полагалось равным 1, 2 и 3, а параметр п - варьировался в диапазоне от 0,1 до 0,9.

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

соответствии с результатами экспериментов, описанными в [13], где было показано, что ne(x) вдоль такой РК имеет ярко выраженный колебательный характер. Для сравнения на

рис. 5 приведено рассчитанное распределение потенциала, соответствующее другим граничным условиям на боковых границах и на катодной плате. Здесь было принято, что при х = 0 и х = d потенциал изменяется как ф(у) = - 900 у/Ь (от 0 на аноде до - 900 В при у = Ь). На катодной плате для концентраций заряженных частиц было принято: пе(х) = 2109 см-3, пр(х)= 1010 см-3, пп(х) = 8109 см-3, ф(х) = изменялось линейно от 900 В при х = 0 и х = 8 см до - 1400 В на пластинчатом катоде.

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

Проведенные расчеты показывают сильную зависимость полученных распределений от характера граничных условий. Расчеты показывают также, что в отличие от результатов одномерного анализа [5, 6] подобных РК, поле Е максимально вблизи катода и убывает в направлении к аноду. Такая же картина была получена ранее в экспериментах [13]. Это свидетельствует о том, что ввиду сильной неоднородности электрического поля в подобных РК одномерный подход здесь является, по-видимому, неадекватным.

Расчеты проводились для скоростей потока в диапазоне (30-250) м/с. Напряжение между анодом и катодом задавалось, как и в эксперименте, равным 1400 В при статическом давлении газа р в диапазоне (4-10) кПа и температуре Т = 293 К. Оказалось, что характер решения сильно зависит не только от вида граничных условий, но и от значений коэффициентов амбиполярной диффузии Беи, Бра, Бпа- Согласно проведенным расчетам, степень отклонения от квазинейтральности в пределах положительного столба достигала 3-4 %, а в приэлектродных областях сильно зависела от характера граничных условий и могла быть на порядок больше. В качестве примера на рис. 2 приведено рассчитанное распределение потенциала внутри РК. Видно, что напряженность электрического поля вблизи катода в несколько раз больше, чем в окрестности анода.

Рис. 2. Распределение концентрации электронов внутри разрядной камеры

Рис. 3. Распределение концентрации положительных ионов внутри РК

Рис. 4. Распределение концентрации отрицательных ионов внутри РК

Рис. 5. Распределение электрического потенциала внутри РК

Рис. 6. Изолинии распределения электрического потенциала внутри РК

Рис. 7. Изолинии распределения концентрации отрицательных ионов в РК

Рис. 8. Распределение концентрации электронов внутри РК

Рис. 9. Распределение концентрации положительных ионов внутри РК

Рис. 10. Распределение отрицательных ионов

Рис. 11. Распределение электрического потенциала

Рис. 12. Электрический потенциал внутри РК

(к=1, п = 0.8)

Рис. 13. Электрический потенциал внутри РК (к = 2, п = 0.8)

Рис. 13. Электрический потенциал внутри РК Рис. 14. Электрический потенциал внутри РК (k = 3, п = 0.8) (k = 1)

Таким образом, методом переменных направлений численно решена двумерная система уравнений зарядовой кинетики и уравнения для электрического потенциала внутри РК с поперечным ТРП. При этом значения коэффициентов ионизации и прилипания электронов к молекулам вычислялись в результате расчета ФРЭЭ. Достигнуто качественное и полуколичественное согласие с имеющимся экспериментом. Численный метод оказался устойчивым при варьировании коэффициентов уравнения (8-10)в пределах: 0 < аг < 1.5- 104 c-1, аа < 1.2104 c-1, q < 5 1020 м3с-1. Данный метод может быть рекомендован для решения более полной системы уравнений, включащей уравнения для газовой и колебательных температур или населенностей колебательных уровней молекул.

Список литературы

1. R. Safioulline. Analytic and numerical calculations of thermal and electric characteristics of glow discharge chambers in gas flow. // International Conference on Atomic and Molecular Pulsed Lasers IV, Tomsk, September 2001. Proceedings of SPIE, Vol. 4747, 2002. - P. 282-285.

2. Сафиуллин Р.К. Распределение концентраций заряженных частиц и потенциала электрического поля в тлеющем разряде в потоке газа. // Известия вузов (Проблемы энергетики), 2002, № 1-2. - С. 69-77.

3. Сафиуллин Р.К. К расчету камер с поперечным тлеющим разрядом в потоке электроотрицательного газа. // Известия вузов (Проблемы энергетики), 2003, № 3-4. - С. 140-145.

4. Сафиуллин Р.К. Математическое моделирование процессов в низкотемпературной плазме тлеющего разряда применительно к СО2- и СО-лазерам. // Диссертация на соискание ученой степени доктора физико-математических наук. - Казань: КГТУ им. А.Н. Туполева, 2006.

5. Голубев В. С., Пашкин С. В. Тлеющий разряд повышенного давления. - М.: Наука, 1990. - 336 с.

6. Велихов Е.П., Ковалев А.С., Рахимов А.Т. Физические явления в газоразрядной плазме. - М.: Наука, 1987. - 160 с.

7. Басыров Р.Ш., Гайсин Ф.М., Миннигулов А.М., Тимеркаев Б. А. Пространственное рапределение параметров тлеющего разряда в потоке электроотрицательного газа. // ТВТ, 1994, 32. - С. 334-338.

8. Rogoff G.L. Ambipolar Diffusion Coefficients for Attachig Gases. // J. Phys.D. Appl. Phys. 18, 1985,№ 8. - P. 1533-1545.

9. Сафиуллин Р.К. Амбиполярная диффузия в газоразрядной плазме. // Известия вузов (Проблемы энергетики), 2000, № 5-6. - С. 110-112.

10. Басыров Р.Ш., Тимеркаев Б.А. Модель тлеющего разряда в поперечном потоке электроотрицательного газа. // ТВТ, 28, 1989, № 1. - С. 30-34.

11. Самарский А.А. Теория разностных схем. - М.: Наука, 1977. - 656 c.

12. Даутов Г.Ю., Тимеркаев Б.А. Генераторы неравновесной газоразрядной плазмы. -Казань: ФЭН, 1996.

13. Миннигулов А.М. Характеристики электронного газа и распределение потенциала электрического поля в тлеющем разряде в потоке воздуха. // Дисс. на соискание ученой степени кандидата физико-математических наук. - Казань: КАИ, 1980. - 186 с.

Safiullin R.K. - doctor of physical-mathematical sciences, professor E-mail: rksaf@mail.ru

Kazan State University of Architecture and Engineering

The organization address: 420043, Russia, Kazan, Zelenaya st., 1

Numerical investigation of glow discharge characteristics in transverse discharge chambers in electronegative gas flow

Resume

Transverse glow discharge chambers in electronegative gas flow were investigated numerically by alternating direction method. A case of continuous anode and partitioned cathode was considered. As a result, two-dimensional distributions of charged particles and of electric potential inside such chambers were obtained. These investigations continue those in the papers [1-3].

A system of elliptic equations (4), (7-10) was solved by alternating direction methodll in two-dimensional subspace XY. Calculations were carried out for a rectangle with sides a and b which contained 1 mm width lamellate cathode and continuous anode. The uniform difference grid (100x100) was used.The boundary conditions were taken as in [7]. They were also widely varied with a view to evaluate their influence on the solution.

Under given boundary conditions the program elaborated enables to calculate twodimensional distributions of electric potential and charged particle densities inside discharge chamber. The calculations showed the strong dependence of these distributions on boundary conditions. They showes also that in contrary with the results of one-dimensional analysis in [6], E has a maximum near cathode and decreases towards anode. The same result was also

received in [13]. It indicates that in view of strong nonhomogeneity of such GDF the one-

dimensional approach evidently turns out to be inadequate.

The calculations were carried out for gas velocities in the range 30-250 m/sec. The

voltage drop between electrodes was put equal 1400 V, the static pressure p at entry was varied in the range 4-10 kPa, and gas temperature T at entry was put equal 293 K. The qualitative and semi-quantitative agreement with experiment [13] was achieved.The calculations showed also that space distributions of charged particles strongly depend on ambipolar diffusion coefficients Dea, Dpa, Dna. As was proposed in [4] ambipolar diffusion coefficient for negative ions Dna could be even negative. The difference scheme of alternating direction method gives diverging solution at negative values of Dna. As follows from our calculations, the deviation from plasma quasi-neutrality may achieve 3-4 % within the positive column, whereas near electrodes it may be an order of magnitude more (depending on boundary conditions). As an example, the calculated space distributions of electron and ion densities and also of electric potential inside a discharge chamber are presented in figs. 2-14.

Keywords: glow discharge, gas flow, discharge chamber, spatial distribution of charged particles, electric potential.

References

1. R. Safioulline. Analytic and numerical calculations of thermal and electric characteristics of glow discharge chambers in gas flow. // International Conference on Atomic and Molecular Pulsed Lasers IV, Tomsk, September 2001. Proceedings of SPIE, Vol. 4747, 2002. - P. 282-285.

2. Safiullin R.K. Charged Particles and Electric Potential Space Distributions in Glow Discharge in Gas Flow // Izvestija Vuzov (Problems of Energetics), 2002, № 1-2. - P. 69-77.

3. Safiullin R.K. Investigation of Discharge Chambers with Transverse Glow Discharge in Electronegative Gas Flow // Izvestija Vuzov (Problems of Energetics), 2003, № 3-4. - P. 140-145.

4. Safiullin R.K. Mathematical Simulation of the Processes in Low Temperature Plasma in Glow Discharge in Connection with CO2 and CO Lasers. // Doctor Dissertation. - Kazan, 2006.

5. Golubev V.S., Pashkin S.V. High Pressure Glow Discharge. - M.: Nauka, 1990. - 336 p.

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

6. Velichov E.P., Kovalev A.S., Rachimov A.T. Physical Phenomena in Gas Discharge Plasma. - M.: Nauka, 1987. - 160 p.

7. Basyrov R.Sh, Gaisin F.M., Minnigulov A.M., Timerkaev B.A. Space Distribution of Glow Discharge Characteristics in Electronegative Gas Flow. // Teplophysica Vysokich Temperatur, 1994, V. 32. - P. 334-338.

8. Rogoff G.L. Ambipolar Diffusion Coefficients for Attaching Gases. // J. Phys.D.- Appl. Phys. 18, 1985, № 8. - P. 1533-1545.

9. Safiullin R.K. Ambipolar Diffusion in Gas Discharge Plasma. // Izvestija Vuzov (Problems of Energetics), 2000, № 5-6. - P. 110-112.

10. Basyrov R.Sh., Timerkaev B.A. The Model of Transverse Glow Discharge in Electronegative Gas Flow. // Teplophysica Vysokich Temperatur, 1989, V. 28, № 1. - P. 30-34.

11. Samarsky A.A. The Theory of Difference Schemes. - M.: Nauka, 1977. - 656 p.

12. Dautov G. Ju, Timerkaev B.A. The Generators of Noneqvilibrium Gas Discharge Plasma. - Kazan: FAN, 1996. - 340 p.

13. Minnigulov A.M. Characteristics of Electron Gas and Distribution of Electric Potential in Glow Discharge in Air Flow. // Candidate Dissertation. - Kazan: KAI, 1980. - 186 p.

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