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

УДК 519-711-3

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

OBTAINING PROPERTY MANAGEMENT MODELS IN TRAFFIC UNDER RAW

SLURRY VIBRATION EXPOSURE

С.В. Андрущак, П.В. Беседин S.V. Andruschak, P.V. Besedin

Белгородский государственный технологический университет им. В.Г. Шухова, Россия, 308012, Белгород, ул. Костюкова, д. 46

Belgorod state technological university after V.G. Shukhov, 56 Kostikova St, Belgorod, 308012. Russia

e-mail: [email protected], [email protected]

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

Resume. Currently, an important area of research is the development of energy-saving techniques for the cement industry. Methods of intensification of transporting sludge through mechanical action effectively influence increase mobility and improve the rheological properties of cement slurries. However, these effects are drifting extreme dependence, for which the study was developed a mathematical model in the form of a nonlinear dynamical system, which includes a model of the actuator, an extreme static characteristics, as well as a measuring device. The resulting model can be applied for designing adaptive in the presence of extreme static characteristic drift control system napimer, in case of change of sludge moisture.

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

Keywords: cement residue, a vibration control object, an extreme model.

Введение

Снижение влажности шлама является существенным резервом повышения производительности вращающихся печей и снижения расхода топлива. Однако, при снижении влажности шлама необходимо сохранить его подвижность на прежнем уровне [Беседин и др., 2011]. Для достижения данного результата существует множество различных факторов, влияющих на структурно-механические свойства шламов, одним из которых является вибрационное воздействие. В результате экспериментальных исследований [Беседин и др., 2015] было установлено, что вибрационное воздействие позволяет снизить влажность транспортируемых шламов на 5...9% при воздействии в диапазоне частот от 10 до 35 Гц, т.е. данная зависимость вязкости шлама имеет экстремальный характер (рис. 3). При этом заранее неизвестно точное оптимальное значение частоты, так как цементные заводы имеют сырьевые материалы различной минералогической природы, а также это значение может «смещаться» на конкретном цементном заводе из-за изменений реологических свойств шламов.

Для решения данной задачи рассмотрим участок шламового питателя (рис. 1), на котором шлам из шламового бассейна 1 поступает во всасывающий патрубок 2 центробежного шламового насоса 3, при этом вибратор 4 целесообразно установить на соответствующий всасывающий патрубок, так как такая компоновка позволяет снизить сопротивление движению жидкости к шламовому насосу и повысить его производительность [Беседин и др., 2015].

Рис. 1. Схема установки вибратора в шламовый питатель Fig. 1. The scheme of installation of the vibrator in a slurry feeder

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

Целью исследования является разработка математической модели объекта управления, включающей экстремальную статическую характеристику при движении сырьевого шлама под вибрационным воздействием, а также модель исполнительного механизма и измерительного устройства. В общем случае рассматриваемый объект экстремального управления можно представить в виде нелинейной динамической системы, структуру которой можно изобразить в виде, показанном на рис. 2 [Рубанов, Бушуев, 2012].

w(t) ^

u ш AP z

ИМ [ эх\ ДЧ

Рис. 2 Структура о&ьекта экстремального управления Fig. 2. The structure of the object of extreme control

На вход исполнительного механизма (ИМ) поступает управляющее воздействие, которое передается на объект управления, представленный в виде нелинейной экстремальной статической характеристики (ЭХ), которая при дрейфе является также функцией времени [Казакевич, Родов, 1977]. Измерительное устройство представлено датчиком (ДЧ).

В общем случае вязкость шлама п будет зависеть от частоты вибрации ш, влажности шлама w, плотности шлама р и его химического состава X:

ц = f (œ, w, р, X ). (1)

Однако, если рассматривать решение задачи регулирования для конкретного цементного завода, то с некоторым допущением можно утверждать, что физико-химические свойства перекачиваемого шлама будут постоянными. Кроме того, влажность и плотность шлама взаимосвязаны между собой некоторой характеристикой, в этой связи для оценки вязкости шлама, достаточно учитывать только влажность сырьевой смеси w(t), которая может вносить дрейф экстремума статической характеристики [Александров, 1989].

В итоге уравнение для вязкости шлама имеет вид:

р = f (œ, w), (2)

где a(t) - частота вибрации вибратора, w(t) - влажность шлама, которая зависит от времени, так как в процессе подачи сырьевой смеси возможно ее изменение.

Перевод объемного расхода шлама в перепад давления

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

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

Потеря напора в трубопроводе в общем случае зависит от сопротивления трения и местных сопротивлений [Касаткин, 1971].

к = к + к , (3)

п тр м.с.'

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

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

Для ламинарного движения по прямой трубе потери напора на трение ктр можно определить на основании уравнения Пуазейля [Романков, Курочкина, 1982]:

Г , (4)

128«

где ц - вязкость жидкости, I - длина круглой прямой трубы, ё - диаметр трубы, АР - перепад давления на участке трубы.

Согласно уравнению Бернулли для горизонтального трубопровода постоянного сечения напор, теряемый на трение:

Р1 - Р2 АР ,

--2 =— = Нтр , (5)

Рё Рё

где р - плотность жидкости, д- ускорение свободного падения.

При подстановке АР=рдктр в уравнение (4) и замене объемного расхода произведением средней скорости потока ш на площадь поперечного сечения трубы пё2/4 получим:

лй2 лй^РёКр

(О-=--. (6)

4 128«

Отсюда находим потерю напора:

32о1«

Умножая числитель и знаменатель правой части на 2ш и группируя величины, получим:

64« I (О2 64 I о2

к = ——----=-----. (8)

тр ойр й 2ё Яе й 2ё

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

, . I о2

к = Л----, (9)

тр й 2ё' У где Х=б4/Яв - коэффициент гидравлического трения для ламинарного потока. Для перепада давления АРтр (с учетом АР=рдктр):

АР = Л- 1- РО(10) тр й 2

Уравнение (10), в первом приближении, может быть использовано для определения перепада давления при турбулентном движении жидкости. После обобщения результатов экспериментов методом теории подобия была получена зависимость коэффициента X для турбулентного потока [7]:

. 0,316 _ ойр

Л = ^=, где Яе =-. (11)

ктр =-1Г • (7)

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

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

J— = Мд - Mc, (12)

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

Уравнение движения трехфазного асинхронного электродвигателя при нулевых начальных условиях, а с учетом только электромеханических переходных процессов [Невраев, Петелин, 1964], можно записать в виде:

da dt

где J - момент инерции, приведенный к валу двигателя, ш - угловая скорость двигателя; Мд - вращающий момент двигателя, Мс -момент сопротивления на валу двигателя.

В инженерной практике широко используется запись уравнения движения в виде:

GD2 dn . .

--= МД - МС, (13)

375 dt Д

где GD2 - маховой момент, приведенный к валу двигателя, n - скорость вращения двигателя. Скорость асинхронного двигателя связана со скольжением s:

( = ( (1 - s) или n = nc (1 - s). (14)

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

ds

- J^Ht = МД -Мс; (15)

GD2nc ds ,, ,,

---— = М^ - Мг. (16)

375 dt Д

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

Jd(= МД (у,() - Мс (). (17)

В установившемся режиме отклонения координат равны нулю и уравнение статики примет

вид:

Мд0 - Мсо = 0. (18)

Подставляя уравнение (18) в (17), получим уравнение асинхронно электродвигателя в отклонениях:

(Тдp + 1)А( = кдAy, где p = d; (19)

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

Т = Тд

дМс дМД

J_

i \ .

(20)

Коэффициент усиления двигателя:

да да

ГдМ Д }

У=Уо

У=Уо

ду

к =_(=(о (21)

Д г дМс дМД ^ ■ ( )

v да да

У=Уо

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

Аа к„

ЖД(р) = — =-. (22)

Д ^ Ау Тдр +1 ( )

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

С0-Ю

о

_ GD 2sHnr

^ ыъ. (24)

Получение экстремальной статической характеристики перепада давления

от частоты вибрации

Для построения статической экстремальной характеристики (ЭХ) будем использовать экспериментальные зависимости расхода шлама (табл. 1) при истечении сырьевой смеси Балаклейского цементного завода от частоты вибрационного воздействия с рабочим органом поршневого типа через патрубок постоянного сечения 15 мм.

Таблица 1 Table 1

Исходные данные зависимостей объемного расхода шлама от частоты вибрации Baseline volume of slurry flow dependency from the vibration frequency

№ w = 37,1% w = 38,1% w = 4 ^1,6%

Ы1 V1 Ы2 V2 Ы3 V3

1 3 0,25 3 0,7625 3 1,2375

2 7 0,33 6 0,95 6 1,425

3 11 0,4 9,5 1,1 9,5 1,675

4 16 0,46 14 1,1875 14 1,8375

5 21 0,49 19,5 1,2625 19,5 1,9375

6 26 0,46 24 1,2375 24 1,9125

7 31 0,39 31,5 1,1375 31,5 1,8625

8 40 0,38 40 0,95 40 1,675

9 49 0,24 49 0,65 49 1,4125

Аппроксимируем данные с помощью функции МаЙаЪ polyfit(ш,V,n), которая находит коэффициенты полинома степени п и аппроксимирует функцию V(x) в смысле метода наименьших квадратов [Беседин и др., 2015]:

V (а) = ап ап +... + а3 а3 + а2 а2 + ах а + а0, (25)

где ш - частота вибрации, Гц; a - коэффициенты аппроксимирующего полинома; V - объемный расход шлама, л/мин.

На рис. 3 приведена графическая иллюстрация, на которой точками обозначены экспериментальные данные согласно табл. 1.

Частота. Гц Частота. Гц Частота. Гц

а б в

Рис. 3. Графическая иллюстрация аппроксимации данных полиномами для трех значений влажности:

а - 3 степени, б - 4 степени, в - 5 степени Fig. 3. Graphic illustration of the data polynomial approximation for the three humidity values: and - 3 degrees, b - 4 degree to - 5 degrees

Для оценки точности аппроксимации экспериментальных данных рассчитаем среднеквадратические погрешности для каждого полинома:

n

^-Х (V - ))2, (26)

i-1

Таблица 2 Table 2

Оценка среднеквадратической погрешности для каждого полинома Estimation of the mean square error for each polynomial

Степень полинома w = 37,1% w = 38,1% w = 41,6%

3 0,0035 0,0028 0,0029

4 0,0032 0,0024 0,0022

5 0,0007 0,0020 0,0019

Анализируя табл. 2 можно утверждать, что минимальная ошибка аппроксимации достигается для полинома 5 степени, однако, форма получившейся кривой для влажности 37,1% (рис. 2 б) имеет несколько точек перегиба. Поэтому выберем 4 порядок аппроксимирующего полинома, так как он наиболее точно описывает экспериментальные зависимости. В итоге, аппроксимирующие полиномы имеют вид:

V = 0,0000001-ю4 -0,0000105ю3 -0,0000893ю2 + 0,0165721-® + 0,2053385;

V = -0,0000003 ю4 + 0,0000514 ю3 -0,0036931-ю2 + 0,0963513 ю + 0,4262834; (27)

V = -0,0000004 ю4 + 0,0000658 ю3 - 0,0043459 ю2 + 0,1144684 ю + 0,918743 5;

где V полином для зависимости объемного расхода от частоты вибрации с влажностью 37,1%; V2 для 38,1%; и V3 для 41,6%, соответственно.

Далее рассмотрим реализацию модели экстремальной статической характеристики в среде Matlab Simulink (рис. 4) на основе приведенных зависимостей [Черных, 2003]. В данной модели входными величинами являются частота вибрации (frequency), Гц и влажность шлама (humidity),%, а выходной перепад давления на участке трубопровода с вибратором (pressure), Па. В блоках (Extreme w) заложены аппроксимированные экспериментальные характеристики, приведенные в табл. 1. Дрейф экстремума осуществляется изменением влажности шлама. С помощью блоков (Prelookup) и (Interpolation Using Prelookup) реализована линейная интерполяция для расчета перепада давления при промежуточных значениях влажности.

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

и влажности шлама в среде Matlab Simulink Fig. 4. The static model for pressure drop characteristics depending on the frequency of vibration and moisture in the

slurry medium Matlab Simulink

После выполнения моделирования при входных параметрах частоты вибрации от о до 50 Гц и влажности шлама от 37,1 до 41,6% представим выходные данные давления в виде трехмерной поверхности (рис. 5).

Рис. 5. Экстремальная зависимость перепада давления от частоты вибрации и влажности шлама Fig. 5. Extreme dependence of the pressure drop on the frequency of vibration and moisture sludge

Получение модели движения сырьевого шлама и построение графика переходного процесса

Далее получим математическую модель (рис. 6) исполнительного механизма, основным элементом которой является трехфазный асинхронный электродвигатель 5АИ 56 А2 IM со следующими техническими характеристиками (Engine): скорость вращения n = 2700 об/мин; число пар полюсов p = 1; мощность P = 0,18 кВт; напряжение питания U = 380 В; частота сети переменного тока f = 50 Гц. Вращающий момент при номинальном режиме работы:

p 018

М„ = 9550 • — = 975—,— = 0,065 Н ■ м. (28)

Н п 2700

Скорость циклического изменения магнитного потока статора, а именно синхронная скорость двигателя:

П = 60-f = 60-50 = 3000 об/мин. (29)

Р 1

Скольжение электродвигателя имеет значение:

п -n 3000- 2700 Л1

SH = -=-2ППП =0Д. (30)

пс 3000

Значение коэффициента чувствительности двигателя по моменту к изменению частоты питающего напряжения приближено вычисляется по формуле:

кч =

(дМ„ Л

Д

М„ 0,065

Н _ 3

= 1,3-10-3 Н-м/Гц. (31)

/ Л / 50

Маховой момент нагрузки (вибронаконечник) составляет GD = 0,055 Нм.

В итоге, значение постоянной времени и коэффициент усиления двигателя вычисляются в соответствии с выражениями (23) и (24):

= (1 - ОД).»,0552-3000 =

Д 375-0,065

(1 - 0,1) - 0,17 -10-3 -3000 ^

кД -'-¿—^-= 54.

Д 0,065

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

Передаточная функция датчика обратной связи (индукционного дифференциального манометра) имеет вид:

к 1

lMs +1 0,05s +1

172

НАУЧНЫЕ ВЕДОМОСТИ

Ниже представлен пример модели (рис. 6), созданной в среде Simulink, для исследования динамики экстремального управления подвижностью шлама при вибрационном воздействии [Кривилев, 2005].

Рис. 6. Структурная модель движения сырьевого шлама при вибрационном воздействии как объект управления

в среде Matlab Simulink

Fig. 6. The structural model of the movement of raw slurry with vibration exposure as a control object in the environment

Matlab Simulink

Экстремальная статическая характеристика, модель которой изображена на рис. 5, реализована с помощью блока (Extreme), на вход которого подаются сигналы влажности (humidity) и частоты вибрации вибратора (frequency). В рассматриваемой модели влажность в качестве примера задана константой 40% и не меняется на протяжении процесса моделирования. Реакция модели на единичное ступенчатое воздействие в момент времени моделирования 1 секунда приведен на рис. 7. По результатам моделирвоания видно, что движение системы началось еще до подачи входного ступенчатого воздействия и установилось около значения 1,38 Па. Это говорит о том, что сырьевой шлам начал свое движение по патрубку без вибрации.

Рис. 7. Построение графика переходного процесса по результатам моделирования в среде Matlab Simulink Fig. 7. Construction of the schedule of the transition process on the results of the simulation in the environment Matlab

Simulink

Выводы

Механическое воздействие вибрирующего органа с частотой вибрации 10-35 Гц позволяет значительно увеличить подвижность цементного шлама, и дает возможность снизить его влажность на 5-9%. Однако, для обеспечения максимальной подвижности шлама необходимо разработать систему управления вибрирующим устройством, которое установлено на всасывающем патрубке шламнасоса. При этом система управления должна учитывать дрейф экстремальной статической характеристики, которая обусловлена изменением влажности шлама и особенностями его движения в гладкой трубе всасывающего патрубка.

Предложена математическая модель экстремального объекта управления, которая включает вибрирующее устройство, экстремальную статическую характеристику, которая описывает зависимость перепада давления от частоты вибрации и влажности шлама, а также дифференциальный манометр. Полученная математическая модель реализована в среде МаЙаЬ Simulink и отражает динамику объекта управления.

НАУЧНЫЕ ВЕДОМОСТИ

Серия Экономика. Информатика. 2016. № В (230). Выпуск 3S

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

Александров А. Г. l989. Оптимальные и адаптивные системы: Учеб. пособие для вузов по спец. Автоматика и упр. в техн. системах. М.: Высш. шк, c. 263.

Aleksandrov A. G. 1989. Optimal'nye i adaptivnye sistemy: Ucheb. posobie dlja vuzov po spec. Avtomatika i upr. v tehn. sistemah. M.: Vyssh. shk, c. 263.

Беседин П.В., Андрущак С.В., Козлов В.К. 20l5. Исследование и моделирование процесса движения цементного шлама в лабораторных условиях. Вестник БГТУ им. В.Г. Шухова, 4: ll3-ll9.

Besedin P.V., Andrushhak S.V., Kozlov V.K. 2015. Issledovanie i modelirovanie processa dvizhenija cementnogo shlama v laboratornyh uslovijah. Vestnik BGTU im. V.G. Shuhova, 4: 113-119.

Беседин П.В., Андрущак С.В., Козлов В.К. 2015. Методы физико-механического регулирования реологических свойств сырьевых шламов. Цемент и его применение, 3: l27-l30.

Besedin P.V., Andrushhak S.V., Kozlov V.K. 2015. Metody fiziko-mehanicheskogo regulirovanija reologicheskih svojstv syr'evyh shlamov. Cement i ego primenenie, 3: l27-l30.

Беседин П.В., Трубаев П.А., Панова О.А., Гришко Б.М. 20ll. Некоторые направления энергосбережения в технологии цемента. Цемент и его применение. 2: l30-l34.

Besedin P.V., Trubaev P.A., Panova O.A., Grishko B.M. 2011. Nekotorye napravlenija jenergosberezhenija v tehnologii cementa. Cement i ego primenenie. 2: l30-l34.

Казакевич, В.В., Родов А.Б. 1977. Системы автоматической оптимизации. М.: Энергия. с. 288.

Kazakevich, V.V., Rodov A.B. 1977. Sistemy avtomaticheskoj optimizacii. M.: Jenergija. s. 288.

Касаткин А.Г. l97l. Основные процессы и аппараты химической технологии. 8-ое изд. М.: Химиздат, с. 783.

Kasatkin A.G. 1971. Osnovnye processy i apparaty himicheskoj tehnologii. 8-oe izd. M.: Himizdat, s. 783.

Кривилев А. 2005. Основы компьютерной математики с использованием системы MATLAB. Лекс-Книга,

с. 496.

Krivilev A. 2005. Osnovy komp'juternoj matematiki s ispol'zovaniem sistemy MATLAB. Leks-Kniga, s. 496.

Невраев В.Ю., Петелин Д.П. l964. Системы автоматизированного электропривода переменного тока. М: Издательство «Энергия». с. l04.

Nevraev V.Ju., Petelin D.P. 1964. Sistemy avtomatizirovannogo jelektroprivoda peremennogo toka. M: Iz-datel'stvo «Jenergija». s. l04.

Ромаиков П.Г., Курочкина М.И. 1982. Гидромеханические процессы химической технологии. Л.: Химия.

c. 288.

Romankov P.G., Kurochkina M.I. 1982. Gidromehanicheskie processy himicheskoj tehnologii. L.: Himija. s. 288.

Рубанов В.Г., Бушуев Д.А. 20l2. Моделирование экстремальных систем управления в среде matlab и simulink как средство анализа динамики. Научные ведомости Белгородского государственного университета. Серия: История. Политология. Экономика. Информатика. l9-l (l38): l69-l75.

Rubanov V.G., Bushuev D.A. 20l2. Modelirovanie jekstremal'nyh sistem upravlenija v srede matlab i simulink kak sredstvo analiza dinamiki. Nauchnye vedomosti Belgorodskogo gosudarstvennogo universiteta. Serija: Istorija. Politologija. Jekonomika. Informatika. l9-l (l38): l69-l75.

Черных, И.В. 2003. SIMULINK: среда создания инженерных приложений. Под общ. Ред. к.т.н. В.Г. Потемкина. М.: ДИАЛОГ-МИФИ, с. 496.

Chernyh, I.V. 2003. SIMULINK: sreda sozdanija inzhenernyh prilozhenij. Pod obshh. Red. k.t.n. V.G. Potem-kina. M.: DIALOG-MIFI, s. 496.

Пат. 20l4ll2882. Российская Федерация, МПК F 27 B 7/32. Способ регулирования подачи шлама в цементную печь. Беседии П.В., Аидрущак С.В., Шарапов Р.Р., Козлов В.К., Новичков С.Г.; заявитель и патентообладатель Белгородский государственный университет; заявл. 02.04.20l4; опубл. 20l5.06.30.

Pat. 2014112882. Rossijskaja Federacija, MPK F 27 B 7/32. Sposob regulirovanija podachi shlama v cementnuju pech'. Besedin P.V., Andrushhak S.V., Sharapov R.R., Kozlov V.K., Novichkov S.G.; zajavitel' i patentoobladatel' Belgo-rodskij gosudarstvennyj universitet; zajavl. 02.04.2014; opubl. 2015.06.30.