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

УДК 621.313

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

А. Б. МЕНЖИНСКИЙ, А. Н. МАЛАШИН, Ю. Г. КОВАЛЬ

Учреждение образования «Военная академия Республики Беларусь», г. Минск

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

Введение

Система энергоснабжения современных робототехнических комплексов требует разработки электромеханических преобразователей энергии с высокими энергетическими и минимальными массогабаритными показателями [1]. В связи с этим в про-мышленно развитых странах (США, России, Великобритании, Японии, ФРГ, Швеции, Нидерландах, Китае, Израиле и др.) в качестве перспективной энергоустановки рассматривается свободнопоршневой двигатель (СПД) с генератором [2]-[5].

Подобными энергоустановками на базе СПД в настоящее время занимается множество фирм и научных университетов, таких как: General Motors, Toyota Central, Sandia National Laboratories (P.V. Blarigan), NASA, Stirling Technology Company (США), Национальный университет науки и технологии Тайваня, университет Тянжина, прикладной институт науки и технологии Кореи, Стэндфордский университет, университет технологий Петронас, университет Ньюкасла и др. [5].

Основной особенностью такой системы является отсутствие кривошипно-шатунного механизма в конструкции двигателя. Это позволяет увеличить его КПД до 50-60 %, в 2,5-3 раза увеличить габаритную мощность, уменьшить удельную массу, металлоемкость СПД [2] и удельный расход топлива приводного двигателя на 30 %, реализовать модульную структуру, увеличить ресурс до капитального ремонта до 50 тыс. ч [3], [4].

В энергоустановках на базе СПД в качестве электрической машины чаще всего применяются возвратно-поступательные электрические генераторы с поперечным приращением магнитного потока (ВПЭГ с ПМП) [2]-[5]. Основным недостатком этих генераторов является отсутствие согласования электрической и механической подсистем в крайних точках рабочего цикла, что ограничивает эффективность использования СПД и снижает надежность энергоустановки.

В [6], [7] для решения проблемы согласования электрической и механической подсистем энергоустановки на базе СПД в крайних точках рабочего цикла было предложено использовать электромеханический преобразователь энергии с поперечным и продольным приращением магнитного потока (комбинированный генератор). Одна из возможных конструкций комбинированного генератора представлена на рис. 1.

Вид в плоскости 20Х Вид в плоскости У02

1 ч /

<1 5 г-

ЦШШ1 4 ■к-:-:-:- У А кап

\ \ / 1_/ < -<

2 3

Рис. 1. Генератор комбинированной конструкции: 1 - рабочая обмотка; 2 - магнитопровод ВПЭГ поперечного типа;

3 - магнитопровод ВПЭГ продольного типа; 4 - изолятор;

5 - общий магнитопровод

Электромагнитная сила генератора комбинированной конструкции (ГКК) принимает максимальное значение в крайних точках рабочего цикла СПД. Благодаря этому появляется возможность обеспечения согласования электрической и механической подсистем энергоустановки на базе СПД на всем рабочем цикле.

Поэтому целью работы является математическое описание ГКК для его дальнейшего исследования и оценки эффективности применения в энергоустановке на базе СПД.

Основная часть

Математическому описанию ВПЭГ поперечного типа посвящено достаточное количество работ [2]-[5] в отличие от ВПЭГ продольного типа [6]-[8], анализ которых показал, что все они основываются на теории цепей (уравнениях Кирхгофа). Основное преимущество цепных методов заключается в том, что построенные на их основе математические модели (ММ) электрических машин позволяют получать ключевые характеристики за малый промежуток времени, поэтому они применяются в задачах оптимизации, позволяя перебирать множество вариантов за ограниченное время. Недостаточная точность таких моделей требует применения более сложных моделей для последующего уточнения полученного результата.

Однако в отличие от электромеханических преобразователей энергии (ЭМПЭ) вращательного типа, возвратно-поступательные преобразователи обладают рядом особенностей: неравномерностью распределения магнитного поля в воздушном зазоре зубцово-пазовой зоны; переменным характером воздушного зазора между подвижной и неподвижной частью генератора и в некоторых случаях разомкнутостью магнитопровода, учет которых имеет важное значение при проектировании возвратно-поступательных преобразователей энергии. Поэтому принятие некоторых упрощений (допущений), характерных для теории цепей, при математическом описании возвратно-поступательных преобразователей энергии может оказаться достаточно грубым приближением, что повлечет за собой неточности в вычислениях их характеристик. В связи с этим для исследования подобных ЭМПЭ целесообразно применять ММ на основе теории поля, использующие численные методы [9]. Это позволяет учесть специфику геометрии машины, насыщение участков магнитопровода, различие магнитных свойств среды, неравномерность воздушного зазора и другие особенности распределения магнитного поля [10], что позволяет описывать процессы, протекающие в возвратно-поступательных преобразователях с высокой достоверностью.

Электромагнитные процессы в ГКК описываются известными уравнениями Максвелла в дифференциальной форме [11], [12]:

(1)

div B = 0;

Г, ^ dD

rot H = у + -;

t E dB rot E=" m;

div D = p,

где В - вектор магнитной индукции; Н - векторы напряженности магнитного поля; ] - вектор плотности тока; Е - вектор напряженности электрического поля; Б -вектор электрической индукции; р - объемная плотность электрических зарядов.

Основные четыре уравнения Максвелла необходимо дополнить системой уравнений, описывающих свойство материалов [11]:

B = Д o^H = д а H ;

D = s 0 sE = s а E; j = УE,

(2)

где д 0 - абсолютная магнитная проницаемость вакуума; д - относительная магнитная проницаемость среды; да - абсолютная магнитная проницаемость среды; в 0 -абсолютная диэлектрическая проницаемость вакуума; в - относительная диэлектрическая проницаемость среды; ва - абсолютная диэлектрическая проницаемость среды; у - удельная электрическая проводимость среды.

В общем случае абсолютная магнитная проницаемость среды, связывающая между собой векторы напряженности и индукции магнитного поля, имеет различные значения по каждой оси [11].

Векторы электромагнитного поля (В, Н, Е и ]) являются функциями не только пространственных координат, но и времени. Причем во времени эти векторы изменяются по произвольному периодическому закону. При такой постановке задачи не удается получить из систем уравнений (1) и (2) аналитическое выражение, которое можно было бы использовать для дальнейших расчетов электромагнитного поля даже численными методами [11].

С учетом этого был принят ряд допущений относительно свойств магнитных материалов и характера протекания электромагнитных процессов.

Первым является допущение о стационарном характере поля. Источниками такого поля являются постоянные токи или постоянные магниты [13]. Стационарное магнитное поле можно рассматривать независимо от стационарного электрического и наоборот (они не влияют друг на друга) [10].

Второе допущение - ферромагнитные сердечники представляются средами с линейными или нелинейными, но изотропными свойствами (д х = д у = д). Это допущение свидетельствует о том, что свойства магнитопровода по различным осям одинаковы. В нелинейной постановке задачи свойства активных материалов ВПЭГ задаются зависимостью В = / (Н) [14].

Третье допущение - магнитная проницаемость постоянна по всей длине магни-топровода (д = const).

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

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

целесообразно записывать через вектор остаточной индукции Br [10]. При этом учитываются следующие допущения [11]:

- вектор остаточной индукции постоянного магнита отличается от нуля только по главной оси намагничивания;

- вектор остаточной индукции зависит только от напряженности магнитного поля по главной оси намагничивания;

- применительно к высококоэрцитивным постоянным магнитам (ПМ) вектор остаточной индукции принимается постоянным в пределах изменения напряженности магнитного поля от нуля до значения, равного коэрцитивной силе по индукции;

- магнитная проницаемость ПМ по всем координатам одинакова и равна магнитной проницаемости по главной оси намагничивания (для высококоэрцитивных ПМ принимается равной д0).

Всю магнитную систему ГКК рассмотрим в виде совокупности следующих областей: область рабочего воздушного зазора; область проводников с током; область магнитопровода; область постоянных магнитов.

С учетом принятых допущений магнитостатическая векторная модель магнитного поля (МП) ГКК на основе уравнений Максвелла приобретает вид:

rot H = j; (3)

div B = 0; (4)

уравнение материальной связи:

B = д a H + Br. (5)

При построении модели на внутренних и внешних границах области задаются нижеперечисленные граничные условия [12].

Условие Неймана. Составляющие магнитного поля B и H можно найти при соблюдении граничных условий неразрывности нормальных и тангенциальных составляющих магнитного поля на границах раздела сред (условие Неймана) с различными магнитными проницаемостями и :

- граничные условия неразрывности нормальных составляющих вектора индукции магнитного поля B+n = Bn (однородное условие Неймана);

- граничные условия неразрывности тангенциальных составляющих вектора напряженности магнитного поля H^ = H~ однородное условие Неймана).

Условие Дирихле позволяет задать на внешней границе модели наперед известный векторный магнитный потенциал A. Это граничное условие характеризует поведение нормальной составляющей вектора индукции на границе модели Bn. В данной задаче зададим нулевое граничное условие Дирихле Bn = 0, для указания полного затухания поля ( A = 0) на удаленной от источников границе.

Возьмем ротор левой и правой части уравнения (5) и получим:

rot8 = rot (| а H ) + rot Br. (6)

Условие непрерывности магнитных силовых линий (div B = 0) позволяет ввести некоторую векторную функцию A (векторный магнитный потенциал) такую, что [10]:

B = rot A. (7)

Подставим уравнение rot H из (3) и B из (7) в уравнение (6) и получим основное уравнение для расчета магнитостатического поля:

rot (rot A) = |а j + rot Br. (8)

Уравнение (8) можно решить численным методом. С учетом того, что [10]:

rot (rot A) = grad(divA)- V2 A (9)

можно записать

grad(divA)- V2 A = / j + rot Br. (10)

Так как для магнитостатического поля линии вектора A замкнуты сами на себя, то divA = 0. [12].

Тогда уравнение Пуассона для векторного магнитного потенциала примет вид:

V2 A = -|J - rot Br. (11)

В двумерной плоскопараллельной задаче вектор индукции B всегда ориентирован в плоскости модели (x, y), а вектор плотности тока j и векторный потенциал A перпендикулярны к ней [13]. Это значит, что отличны от нуля только компоненты jz и Az. Таким образом, в декартовой системе координат уравнение (11) относительно векторного магнитного потенциала для магнитной системы ГКК примет вид:

( д2 Az + д2 Az л

dx2 dy2 )

= -|0jz --^rot Br. (12)

(j,

Уравнения (12) и граничные условия представляют двухмерную магнитостатиче-скую векторную модель ГКК в плоскопараллельной постановке.

Для решения уравнения (12) при такой постановке задачи был выбран метод конечных элементов (МКЭ), реализованный в программном продукте Elcut 5.1. Программный продукт Elcut 5.1 может применяться для решения линейных и нелинейных двумерных задач магнитостатики в плоскопараллельной и осесимметричной постановке, при этом используется формулировка задачи относительно векторного магнитного потенциала.

Имитируя движения подвижной части генератора с некоторым шагом ёхш и интерполируя в программе Matlab/Simulink или Mathcad полученные значения потокос-цепления на один виток ^го контура, определим мгновенные значения:

- ЭДС движения к-го контура:

d ilAzkdsЛ

(13)

f f

^ d l Azkds

E = ~Wkdlt

S

y

где wk - количество витков к-го контура; 1/$ £А2кёя - потокосцепление на один виток к-го контура.

Интегрирование в данной формуле ведется по поперечному сечению обмотки, а $ обозначает площадь этого поперечного сечения;

- собственной индуктивности рабочей обмотки:

I» = - ^, (14)

гк $

где 1/ $ £ Ак^ ё - потокосцепление на один виток к-го контура, созданное током к-го контура /к;

- взаимной индуктивности:

w

| Azk (,, )ds

L,k = y"" , (15)

n

где 1/$ £ Агк^ ё - потокосцепление на один виток к-го контура, созданное током п-го контура гп.

Разработанная магнитостатическая векторная модель ГКК на основе уравнений Максвелла позволяет получить мгновенные значения основных параметров (Е, Ь0, Ьпк) ГКК с учетом принятых допущений и граничных условий, а также специфики геометрии магнитной системы генератора, нелинейности кривой намагничивания материалов, насыщения участков магнитопровода, различия магнитных свойств сред и неравномерности распределения магнитного потока в воздушном зазоре, что способствует повышению точности полученных результатов.

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

Уравнения электрического равновесия для к-го контура магнитоэлектрического ГКК (рис. 1) можно записать в виде:

,=1

= 0, (16)

к

,

где 1н - индуктивность нагрузки; Яок - активное сопротивление к-го контура; Ян -

- ЭДС движения к-го контура.

активное сопротивление нагрузки; > -——

1 СII

,=1 ы,

Уравнению электрического равновесия (16) может быть поставлена в соответствие эквивалентная электрическая схема (рис. 2).

ik

Rn LH

Unk

Рис. 2. Эквивалентная электрическая схема ГКК при работе на линейную нагрузку

Напряжение нагрузки ^го контура можно записать в виде:

U* = ДА + Lн Щ-.

ш

(17)

Подставляя в (16) ЭДС движения и собственную индуктивность, полученные по выражениям (13) и (14), ММ ГКК при линейной нагрузке на основе уравнений Кирхгофа и Максвелла может быть представлена системой уравнений вида:

w,

0k

f Azk(4 d

S

+L

di,

d

+ik (( + Rn Кwokdft

■ r

f Azk (ik )ds

S

- w,

0k

d_ dt

f )ds

S

= 0;

(18)

Unk = Kh + L §,

где

1S f Azk (

Mn=L .fMs

)ds - потокосцеплепие па один виток k-го контура, созданное

п-м ПМ (п = 1...^ ) и учитывающее продольное, поперечное или комбинированное приращения магнитного потока.

Для расчета собственной индуктивности ^го контура и потокосцепления на один виток ^го контура, созданного п-м ПМ (п = 1... необходимо решить уравнение Пуассона (12) относительно компоненты Л2 векторного магнитного потенциала Л с заданными граничными условиями для магнитной системы ГКК.

Вся магнитная система ГКК состоит из следующих областей: область рабочего воздушного зазора (Овозд); область проводников с током (Оток); область магнитопровода

(Ост); в магнитоэлектрических генераторах область постоянных магнитов (О магн ).

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

Геометрия двухмерной модели обобщенной магнитной системы магнитоэлектрического ГКК в плоскопараллельной постановке с граничными условиями представлена на рис. 3.

k

4 = о

А =о

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

На основании (12) получим уравнения для каждой области магнитной системы ГКК: - область постоянных магнитов (Омагн) :

1 Г д2 А, д2 А

До

2 Л \

дх2

ду2

У

= ---го1 Бг; До

- область проводников с током (О ток):

-м

ГдЧ

дх2

ду2

= -л - о;

- область рабочего воздушного зазора (О ):

1 Г д2 А, д2 А

- в

2

ах2

ду2

= 0;

- область магнитопровода (Ост) :

1 Г д2А, д2А, л

-с

дх2

ду2

= о.

(19)

(2о)

(21)

(22)

Уравнения (19)-(22) и граничные условия представляют двухмерную магнито-статическую векторную модель ГКК в плоскопараллельной постановке, решение которых позволит определить компоненту Аг векторного магнитного потенциала А с заданными граничными условиями для каждой области магнитной системы магнитоэлектрического генератора комбинированного типа.

1

В программном продукте Е1еШ 5.1 построена двухмерная конечно-элементная модель МП магнитоэлектрического генератора комбинированного типа (рис. 4).

а) б)

Рис. 4. Двухмерная конечно-элементная модель МП ГКК

Посредством имитации перемещения подвижной части генератора с некоторым шагом Лх двухмерная конечно-элементная модель МП магнитоэлектрического ГКК, представленная на рис. 4, а, позволяет получить дискретную функцию собственной индуктивности к-го контура, а на рис. 4, б - потокосцепления на один виток k-го контура, созданного n-м ПМ, а также значения коэффициентов рассеяния и выпучивания магнитного потока в магнитной системе генератора, в зависимости от координаты перемещения подвижной части генератора и геометрических размеров магнитной системы.

Интерполяция полученных значений собственной индуктивности и потокосцеп-ление на один виток рабочей обмотки в программе Elcut 5.1 была проведена в программе Matlab/Simulink с помощью кубического сплайна одномерной таблицы Look-Up Table [14], что позволило получить закон изменения во времени собственной индуктивности и потокосцепление на один виток рабочей обмотки при изменении координаты положения подвижной части генератора.

На основании полученных результатов и системы уравнений (18) разработана имитационная модель ГКК, структурная схема которой представлена на рис. 5.

15 $ Azk Fm-Fu )ds

Рис. 5. Структурная схема имитационной модели ГКК

На рис. 6 представлены временные диаграммы мгновенных значений мощности, тока и напряжения на выходе ГКК, а также мгновенные значения мощности, тока и напряжения на выходе ГКК, полученные посредством ММ на основе уравнений Кирхгофа (16) и (17).

P, Вт 800 еоо

400 200 о

I, А 5

0

-5

и, В 200

0

-200

0,05 0,052 0,054 0,056 0,058 0,06 0,062 0,064 0,066 0,068 ^ С

Рис. 6. Временные диаграммы мгновенной мощности, тока, напряжения на выходе ГКК при линейной нагрузке: 1 - линейная ММ на основе уравнений Кирхгофа с расчетом магнитных проводимостей воздушных зазоров по методу Ротерса; 2 - линейная ММ на основе уравнений Кирхгофа с уточняющими коэффициентами;

3 - нелинейная ММ на основе уравнений Кирхгофа и Максвелла

Под номером 1 представлены временные диаграммы, полученные посредством линейной ММ на основе уравнений Кирхгофа с расчетом магнитных проводимостей воздушных зазоров по методу Ротерса. Под номером 2 даны временные диаграммы, полученные посредством линейной ММ на основе уравнений Кирхгофа с расчетом магнитных проводимостей воздушных зазоров по методу Ротерса, учитывая при этом геометрию магнитной системы генератора, потоки рассеяния и выпучивания магнитной системы, посредством соответствующих коэффициентов рассчитанных МКЭ. Под номером 3 преведены временные диаграммы, полученные посредством нелинейной ММ на основе уравнений Кирхгофа и Максвелла (18).

Приняв за истинные значения результаты, полученные посредством нелинейной ММ на основе уравнений Кирхгофа и Максвелла (18), относительная погрешность расчетов активной мощности ГКК по линейной ММ без уточняющих коэффициентов и линейной ММ, учитывающей геометрию, потоки рассеяния и выпучивания магнитной системы генератора, представлена на рис. 7.

др.0 о 31

АР, % 5

Рис. 7. Относительная погрешность расчетов активной мощности генератора: а - по линейной ММ без уточняющих коэффициентов; б - по линейной ММ с уточняющими коэффициентами

Из рис. 7, а видно, что, учитывая в линейной ММ на основе уравнений Кирхгофа геометрию, магнитные потоки рассеяния и выпучивания магнитной системы генератора, посредством соответствующих коэффициентов рассеяния и выпучивания, точность полученных результатов возрастает не менее чем на 15 % по сравнению с традиционными [15] линейными ММ на основе уравнений Кирхгофа, в которых расчет магнитных проводимостей воздушных зазоров осуществляется по методу Ротерса. Кроме того, из рис. 7, б видно, что расхождение результатов, полученных по нелинейной и линейной ММ, учитывающей геометрию, потоки рассеяния и выпучивания магнитной системы генератора, не превышает 5 %, что является приемлемой для большинства инженерных расчетов точностью. Основным недостатком нелинейной ММ является сложность ее применения для решения задач оптимизации и управления, поэтому ее целесообразно использовать на завершающих этапах проектирования с целью уточнения полученных результатов.

Заключение

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

Литература

1. Военно-патриотический сайт «Отвага» Российской Федерации. - Режим доступа: www.otvaga2004.ru/na-zemle/na-zemle-11/modern_land_robots_1/. - Дата доступа: 28.10.2017.

2. Пинский, Ф. И. Энергоустановки со свободнопоршневыми двигатель-генераторами / Ф. И. Пинский // Бортовая энергетика. - 2004. - № 2. - С. 13-17.

3. Cawthorne, W. R. Optimization of a Brushless Permanent Magnet Linear Alternator for Use with a Linear Internal Combustion Engine: Diss. College Eng. and Mineral Resources / W. R. Cawthorne. Morgantown, 1999. - 113 р.

4. Темнов, Э. С. Разработка теоретических основ расчета и конструирования малоразмерных двигатель-генераторных установок как единой динамической системы : дис. ... канд. техн. наук : 05.04.02 / Э. С. Темнов. - Тула, 2005. - 134 л.

5. Hanipah, M. R. Recent commercial free-piston engine developments for automotive applications. Applied Thermal Engineering / M. R. Hanipah, R. Mikalsen, A. P. Roskilly. -2015. - Р. 493-503.

6. Применение возвратно-поступательного генератора комбинированной конструкции для повышения КПД и уменьшения удельной массы энергоустановок автономных образцов вооружения / А. Б. Менжинский [и др.] // Вестн. Воен. акад. Респ. Беларусь. - 2017. - № 4 (57). - С. 62-72.

7. Применение активного выпрямителя с возвратно-поступа-тельными генераторами комбинированной конструкции для повышения эффективности энергоустановок автономных объектов / А. Б. Менжинский [и др.] // Магист. Вестн. Нац. акад. наук Респ. Беларусь. - 2017. - С. 40-50.

8. Использование возвратно-поступательной схемы электрического генератора для повышения эффективности энергоустановок автономных образцов вооружения / А. Б. Менжинский [и др.] // Вестн. Воен. акад. Респ. Беларусь. - 2016. - № 4 (53). -С.108-114.

9. Копылов, И. П. Математическое моделирование электрических машин / И. П. Копылов. - М. : Высш. шк., 2001. - 327 с.

10. Буль, О. Б. Методы расчета магнитных систем электрических аппаратов: Магнитные цепи, поля и программа FEMM : учеб. пособие для студентов высш. учеб. заведений / О. Б. Буль. - М. : Академия, 2005. - 336 с.

11. Ледовский, А. Н. Электрические машины с высококоэрцитивными постоянными магнитами / А. Н. Ледовский. - М. : Энергоатомиздат, 1985. - 169 с.

12. Кулон, Ж.-Л. САПР в электротехнике ; пер. с фр. / Ж.-Л. Кулон, Ж.-К. Сабоннадь-ер. - М. : Мир, 1988. - 203 с. : ил.

13. Сочава, М. В. Решение полевых задач с помощью программы ELCUT 6.0. Задачи магнитостатики и магнитного поля переменных токов : учеб. пособие / М. В. Со-чава. - СПб., 2014. - 38 с.

14. Нейман, Л. А. Решение задачи учета нелинейных свойств динамической модели электромагнитного привода / Л. А. Нейман, А. С. Шабанов, В. Ю. Нейман // Материалы XIX Международной научно-практической конференции, Москва, 7-8 окт. 2015 г. - М., 2015. - С. 58-62.

15. Иванов-Смоленский, А. В. Электрические машины : учеб. для вузов / А. В. Иванов-Смоленский. - М. : Энергия, 1980. - 928 с.

Получено 19.03.2018 г.