Статья опубликована в №31 (март) 2016
Разделы: Сельское хозяйство
Размещена 20.03.2016. Последняя правка: 31.05.2016.
Просмотров - 2147

ПАРАМЕТРИЧЕСКАЯ ИДЕНТИФИКАЦИЯ ДИНАМИЧЕСКОЙ МОДЕЛИ ПЕРЕХОДНЫХ ПРОЦЕССОВ В ПОРШНЕВЫХ ГИДРОНАСОСАХ И ГИДРОМОТОРАХ

Войтов Антон Викторович

асистент

Харьковский Национальный технический университет сельского хозяйства им. П.Василенка

аспирант, кафедра технических систем и технологий животноводства им. Б.П. Шабельника

Бойко Иван Григорьевич, профессор, кандидат технических наук, Харьковский Национальный технический университет сельского хозяйства им. П.Василенка, Кафедра технических систем и технологий животноводства им. Б.П.Шабельника

УДК 621.891

Актуальность проблемы

Техническое состояние гидростатических трансмиссий (ГСТ), которое выражается в изменении коэффициента подачи насоса и КПД мотора влияет на уменьшение производительности сельскохозяйственных машин при одновременном увеличении расхода топлива [1]. Техническое состояние таких агрегатов гидропривода как аксиально-поршневой насос (НП-90,112) и мотор (МП-90,112), во многом определяет производительность  машин в целом, т.к. влияет на скорость их движения. ДСТУ 2193-93 определяет предельные значения коэффициента подачи для насоса и общего КПД для мотора, при достижении которых эксплуатация гидропривода должна быть прекращена. Указанные параметры не должны снижаться более чем на 20% от начального состояния [2].

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

Анализ публикаций, посвященных данной проблеме 

Автором работы [3] приведен анализ технического состояния ГСТ-90, выявлены закономерности распределения износов и получена математическая модель связи объемного КПД с износами и зазорами в сопряжениях. В работах [4,5] обосновывается одно  из направлений диагностирования гидроприводов - термометрирование, которое позволяет определить техническое состояние по температуре корпуса насоса и мотора во время эксплуатации, где по результатам измерения температуры диагностируются отдельные узлы гидронасоса. Данному подходу посвящены работы [6,7], в которых представлены результаты по исследованию температурного режима и разработана методика графо-аналитичского прогнозирования ресурса гидропривода.

Авторы работы [8] разработали методику определения перерасхода топлива при изменении технического состояния гидропривода. На основании полученных исследований делается вывод о необходимости пересмотра граничных значений КПД в процессе эксплуатации.  

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

В работе [10] выполнена структурная  идентификация объемного гидропривода, как объекта диагностирования. Переходной процесс в насосе и моторе описывается колебательными звеньями второго порядка. На основании передаточных функций получены дифференциальные уравнения переходного процесса, анализ которых позволил установить, что на динамику переходного процесса в насосе влияет скорость отклонения наклонной шайбы насоса, а на динамику переходного процесса в моторе – скорость нарастания давления в контуре нагнетания.

Характер протекания переходного процесса после приложения к системе входного воздействия (отклонения наклонной шайбы на угол α) зависит от коэффициентов усиления К_н и К_м, постоянных времени Т_н и Т_м, а также значений декремента затухания d_н и d_м. Авторами делается вывод, что числовые значения К_i, Т_i, d_i зависят от технического состояния насоса и мотора, а следовательно учитывают утечки рабочей жидкости в сопряжениях. Поиск зависимостей коэффициента подачи насоса и к.п.д. мотора с перечисленными выше параметрами позволит диагностировать техническое состояние насоса и мотора.

В данной работе приводятся решения дифференциальных уравнений переходного процесса для насоса и мотора, которые позволяют моделировать переходный процесс во времени при различных значениях К_i, Т_i, d_i. Для выполнения моделирования необходимо выполнить параметрическую идентификацию гидропривода, как объекта диагностирования. 

Цель исследования. Выполнить параметрическую идентификацию математической модели диагностирования объемного гидропривода ГСТ-90,112, как единой динамической системы насос-мотор-гидравлическая жидкость с учетом утечек в процессе эксплуатации. 

Методический подход в проведении исследований 

В процессе параметрической идентификации необходимо найти значения параметров T_i, d_i, К_i, чтобы при их подстановке в уравнения, которые приведены в работе [10],  правая и левая часть отличались наименьшим образом.

На основании анализа литературных источников, а также статистической информации собранной в эксплуатации, были получены выражения для определения коэффициентов усиления уравнения, которое описывает динамику переходного процесса в поршневых насосах НП-90, 112.

Коэффициент усиления К₁ характеризует степень влияния входного сигнала (скорость  отклонения наклонной шайбы ά) на выходной сигнал (скорость нарастания давления в нагнетательном контуре после насоса Ṗ_н) и определяется как:

(1)

где V_н – рабочий объем насоса, м³;

Р_н – номинальное давление в нагнетательной магистрали за насосом, Па;

Е_ж - модуль упругости рабочей жидкости, Па;

d_ос – диаметр, на котором расположены оси поршней, м.

 Как следует из выражения (1) коэффициент усиления К₁ является безразмерным и учитывает конструктивные особенности поршневого насоса в виде номинального рабочего давления Ṗ_н, эквивалентного зазора d_ос и упругих свойств рабочей жидкости Е_ж.

            В процессе эксплуатации насоса происходит снижение степени влияния входного сигнала ά на выходной сигнал Ṗ по причине изнашивания подвижных сопряжений. Следовательно, коэффициент К₁ является функцией коэффициента подачи насоса Ѱ. Данную зависимость можно представить выражением:

(2)

где ΔѰ – величина уменьшения коэффициента подачи насоса в процессе эксплуатации, выраженная в процентах, %.

На основании выражения (2) были получены зависимости изменения безразмерного коэффициента К₁ при изменении коэффициента подачи насоса – увеличении утечек ΔѰ для
НП-90 и НП-112, которые представлены на рис. 1. Как следует из представленных кривых насос с меньшим рабочим объемом, имеет большие значения коэффициента усиления. Однако при уменьшении коэффициента подачи на 20% различия между насосами НП-90 и НП-112 не проявляются. Полученные зависимости подтверждают, что коэффициент усиления К₁ является переменной величиной, которая зависит от коэффициента подачи насоса и будет учитываться при моделировании переходных процессов в насосе.

Коэффициент усиления К₂, характеризует степень влияния входного сигнала (скорость увеличения давления Ṗ_н) на выходной (утечки жидкости Q_y), которые в процессе эксплуатации насоса увеличиваются. Зависимость коэффициента К₂ при увеличении утечек жидкости, можно представить выражением:

                                                                        (3)

где  µ - динамическая вязкость рабочей жидкости, Па·с;

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

 

где n – рабочие обороты вращения ротора насоса, об/мин.

Коэффициент К₂ в литературных источниках назван как критерий изогональности или обратная величина критерия Зоммерфельда, который пропорционален толщине слоя рабочей жидкости в зазорах насоса [11]. 

Результаты исследований 

Используя выражение (3) были построены зависимости изменения коэффициента К₂ при изменении коэффициента подачи насоса ΔѰ, которые представлены на рис. 2.

Как следует из рис. 2, коэффициент, характеризующий утечки жидкости, увеличивается при уменьшении коэффициента подачи насоса. При этом, значения К₂ для насосов
НП-90 и НП-112 одинаковы, т.к. одинаковы величины номинальных давлений Р_н и угловых скоростей вращения у этих насосов. Следовательно, с уменьшением коэффициента подачи насоса в процессе эксплуатации коэффициент усиления степени влияния скорости нарастания давления на утечки жидкости увеличивается, что необходимо учитывать в процессе моделирования.

 

 

Рис. 1. Зависимости изменения коэффициентов усиления К₁ при изменении коэффициента подачи насосов ΔѰ

 

 

Рис. 2. Зависимости изменения коэффициентов усиления К₂ при изменении коэффициента подачи насосов ΔѰ

Переходный процесс в поршневых насосах характеризуется  постоянными времени Т₁ и Т₂. Согласно структурно-динамической схемы, постоянная времени Т₁ характеризует промежуток времени, за который происходит нарастание давления в нагнетательном контуре за насосом, т.е. инерционность переходного процесса в насосе, которую можно выразить как:

(4)

где J – момент инерции ротора насоса, определяется как:

 кг·м2 ,

где m – масса ротора насоса, кг;

R – радиус ротора, м;

А_п – площадь поршня насоса, м².

Как следует из выражения (4) постоянная времени Т₁ имеет размерность – секунда. Зависимости изменения постоянных времени Т₁ для насосов НП-90 и НП-112 при уменьшении коэффициента подачи ΔѰ, представлены на рис. 3. Из зависимостей следует, что постоянная времени для новых насосов, имеющих максимальное значение коэффициента подачи,  имеет значения 0,106-0,134 с. С увеличением утечек ΔѰ, т.е. с уменьшением коэффициента подачи насоса, увеличивается до значений 3,1-3,9 с. Следовательно, в процессе эксплуатации, по причине изнашивания сопряжений и увеличения утечек жидкости, инерционность поршневых насосов увеличивается. Это приведет к запаздыванию нарастания давления в нагнетательном контуре за насосом при отклонении наклонной шайбы.

Постоянная времени Т₂, характеризует промежуток времени, за который происходит утечка жидкости в зазорах сопряжений насоса, т.е. инерционность утечек, которую можно выразить как:

(5)

Как следует из выражения (5) на инерционность утечек влияет динамическая вязкость рабочей жидкости µ и эквивалентный зазор, который зависит от диаметра, на котором расположены оси поршней.

 

 

Рис. 3. Зависимости изменения постоянной времени Т₁ при уменьшении коэффициентов подачи насосов ΔѰ

 

 

Рис. 4. Зависимости изменения постоянной времени Т₂ при уменьшении коэффициентов подачи насосов ΔѰ

Зависимости изменения постоянной Т₂ при уменьшении коэффициента подачи насоса ΔѰ, представлена на рис. 4. Как следует из представленных зависимостей время, за которое происходят утечки в насосах, не зависит от величины коэффициента подачи и его изменение в процессе эксплуатации. При этом, насос меньшего рабочего объема (НП-90) имеет большие значения времени утечек жидкости, т.к. имеет меньший эквивалентный зазор d_ос, что будет положительно влиять на увеличение ресурса.

Полученные выражения зависимостей изменения Т₁ и Т₂ при изменении коэффициента подачи насоса необходимо учитывать при моделировании переходных процессов в насосах.

На основании формул, которые представлены в работе [10], можно получить общую постоянную времени насоса Т_н и декремента затухания насоса d_н. Зависимости указанных выше параметров, при уменьшении коэффициента подачи насоса, представлены на рис. 5 и 6.

 

 

Рис. 5. Зависимости постоянных времени  Т_н при уменьшении коэффициента подачи насосов ΔѰ

 

Рис. 6. Зависимости декрементов затухания  d_н при уменьшении коэффициента подачи насосов ΔѰ

Как следует из зависимостей представленных на рис. 5 инерционность поршневых насосов в процессе эксплуатации увеличивается. Для новых насосов постоянная времени Т_н составляет 0,1 с и одинакова как для НП-90 так и для НП-112. В процессе эксплуатации, за счет увеличения утечек, постоянная времени увеличивается до значений 1,72-1,85 с. При этом, инерционность насоса НП-90 выше, чем НП-112. Из полученных зависимостей можно сделать вывод, что с уменьшением коэффициента подачи насосов, за счет изнашивания сопряжений и следовательно, увеличения утечек, отклонения наклонной шайбы (входной сигнал) не вызывает быстрого нарастания давления в напорной магистрали (выходной сигнал).

Зависимости изменения декремента затухания d_н при изменении коэффициента подачи насоса, рис. 6, позволяют сделать вывод, что увеличение утечек в насосе уменьшает декремент затухания от значений 0,57-0,62 до значений 0,25-0,27, т.е. в 2 раза. Это приведет к появлению колебания величины давления в нагнетательном контуре за насосом, что может служить одним из диагностических признаков, по которому можно оценивать утечки в насосе, а следовательно и величину коэффициента подачи.

Подставив полученные значения Т_н и d_н в уравнения, полученные в работе [10], можно выполнить моделирование переходных характеристик в поршневых насосах НП-90 и НП-112 с учетом изменения величины утечек (уменьшения коэффициента подачи насосов Δψ).

Используя приведенный выше методический подход можно получить выражения для коэффициентов усиления и постоянных времени, которые входят в дифференциальное уравнение динамики переходного процесса в поршневых моторах МП-90 и МП-112 [10].

Из анализа структурно-динамической схемы гидростатического привода [10], коэффициент усиления К₃ характеризует степень влияния входного сигнала (скорость нарастания давления в нагнетательном контуре перед мотором Ṗ_н) на выходной сигнал (угловая скорость вращения ротора мотора ω_м). Значение коэффициента усиления К₃ и характер его изменения при уменьшении к.п.д. мотора можно представить выражением:

                                                                       (6)

где V_м – рабочий объем мотора, м³;

М_кр – крутящий момент на валу мотора, Н·м;

Δŋ – величина уменьшения к.п.д. мотора в процессе эксплуатации, выраженная в процентах, %.

Из выражения (6) следует, что коэффициент усиления К₃ безразмерная величина, которая характеризует удельную энергоемкость гидромотора.

Используя выражения (6) были получены зависимости изменения коэффициентов усиления К₃ при уменьшении к.п.д. МП-90 и МП-112, которые представлены на рис. 7.

Из представленных зависимостей следует, что при уменьшении к.п.д. гидромоторов коэффициенты усиления К₃ уменьшаются, что приведет к уменьшению степени влияния скорости нарастания давления перед мотором на угловую скорость вращения ротора мотора. При этом, значения коэффициентов К₃ для моторов МП-90 и МП-112 одинаковы. Это следует из выражения (6), при увеличении рабочего объема мотора V_м пропорционально увеличивается крутящий момент М_кр, что не приводит к изменению  К₃. Однако, как следует из зависимостей, рис. 7, значения К₃ у новых моторов составляет 5,8 и уменьшаются в процессе эксплуатации до значений 0,3, что существенно снизит степень влияния входного сигнала Ṗ_н на выходной ω_м. Данные зависимости необходимо учитывать при моделировании динамики переходных процессов в поршневых гидромоторах.

Коэффициент усиления К₄ согласно структурно-динамической схемы [10]   характеризует степень влияния входного сигнала (угловая скорость вращения ротора мотора ω_м) на выходной сигнал (утечки рабочей жидкости в моторе Q_у).  В процессе эксплуатации утечки будут увеличиваться по причине изнашивания сопряжений мотора. Выражение для определения величины К₄  при увеличении утечек в моторе можно представить зависимостью:

                                                                    (7)

где ρ_ж – плотность рабочей жидкости, кг/м³.

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

Используя выражение (7) были получены зависимости изменения коэффициентов усиления К₄ при уменьшении к.п.д. моторов МП-90 и МП-112, которые представлены
на рис. 8. 

 

Рис. 7 Зависимости изменения коэффициентов усиления К₃ при уменьшении к.п.д. моторов ∆ŋ

 

Рис. 8 Зависимости изменения коэффициентов усиления К₄ при уменьшении к.п.д. моторов ∆ŋ

Как следует из представленных зависимостей значения коэффициентов К₄ увеличиваются при уменьшении к.п.д. моторов с значений 1-1,2 до значений 2,7-3,0. На основании данных зависимостей можно сделать вывод, что степень влияния угловой скорости вращения ротора мотора ω_м, на утечки в моторе Q_у, в процессе эксплуатации  усиливается. При этом, степень влияния для моторов с большим рабочим объемом выше. Такая характеристика влияния угловой скорости вращения ротора на утечки жидкости в моторе в процессе эксплуатации будет существенно влиять на динамику переходных процессов, что необходимо учитывать при моделировании.

Согласно полученного дифференциального уравнения динамики переходных процессов в поршневых моторах [10], необходимо определить постоянную времени Т_м, которая зависит от Т₃ и Т₄.

Согласно структурно – динамической схемы гидростатического привода, и передаточной функции мотора, постоянная времени Т₃ характеризует промежуток времени, за который происходит увеличение угловой скорости вращения ротора мотора ω_м при нарастании давления в нагнетательном контуре после насоса Ṗ_н, т.е. инерционность переходного процесса в моторе, которую можно выразить как:

                                                                         (8)

Как следует из выражения (8) постоянная времени Т₃ имеет размерность времени – секунда.

Используя выражение (8) были построены зависимости изменения постоянных времени Т₃ при уменьшении к.п.д. моторов МП-90 и МП-112, которые представлены
на рис. 9. Из представленных зависимостей следует, что постоянная времени Т₃ при уменьшении к.п.д. мотора значительно увеличивается с значений 0,58 с до значений 6,39 с.

На основании этого можно сделать вывод, что в процессе эксплуатации моторов, за счет изнашивания сопряжений и увеличения утечек жидкости, инерционность моторов будет увеличиваться. Это приведет к замедлению увеличения угловой скорости вращения ротора мотора ω_м при увеличении давления в нагнетательном контуре после насоса Ṗ_н. При этом характер изменения зависимостей Т₃ для моторов МП-90 и МП-112 одинаков, т.к. при увеличении рабочего объема мотора пропорционально увеличивается крутящий момент.

Постоянная времени Т₄характеризует промежуток времени, за который происходит утечка жидкости в сопряжениях мотора, т.е. инерционность утечек, которую можно выразить как:

                                                               (9)

Анализ размерности выражения (9) показывает, что Т₄ имеет размерность времени – секунда. При этом, на увеличение времени, за которое происходят утечки, влияет увеличение рабочего объема V_м и динамическая вязкость µ, при одновременном уменьшении крутящего момента М_кр и угловой скорости вращения ротора мотора.

Используя выражение (9) были получены зависимости изменения постоянных времени Т₄ при уменьшении к.п.д. для моторов МП-90 и МП-112, которые представлены на рис. 10. 

 

Рис. 9 Зависимости изменения постоянной времени Т₃ при уменьшении к.п.д. моторов ∆ŋ

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

Полученные выражения для определения К₃, К₄, а также Т₃ и Т₄ и их зависимости при уменьшении к.п.д. моторов, позволяют теоретическим путем получить характер изменения постоянной времени мотора Т_м, и декремента затухания d_м. Полученные теоретические зависимости для моторов МП-90 и МП-112 представлены на рис. 11 и 12.

 

Рис. 10 Зависимости изменения постоянной времени Т₄ при уменьшении к.п.д. моторов ∆ŋ

 

Рис. 11 Зависимости изменения постоянной времени T_м при уменьшении к.п.д. моторов ∆ŋ

 

Рис. 12 Зависимости изменения декремента затухания d_м при уменьшении к.п.д. моторов ∆ŋ

Анализ зависимостей на рис. 11 позволяет сделать вывод, что постоянная времени моторов Т_м в процессе эксплуатации, за счет изнашивания и увеличения утечек, увеличивается от значений 0,9 с до значений 7,05 с. Это вызывает значительное увеличение инерционности переходного процесса в моторе, что будет выражаться в длительном увеличении угловой скорости вращения ротора мотора ω_м при изменении скорости увеличения давления в нагнетательном контуре после насоса Ṗ_н. При этом, разница в значениях Т_м для моторов МП-90 и МП-112 не более 2%, при максимально допустимом уменьшении к.п.д. моторов, равном 20%.

Зависимости изменения декремента затухания колебаний угловой скорости вращения ротора ω_м во время переходного процесса при уменьшении к.п.д. для моторов МП-90 и МП-112 представлены на рис. 12. Из представленных зависимостей следует, что в процессе эксплуатации, при снижении к.п.д. мотора, декремент затухания d_м уменьшается, что будет способствовать возникновению колебаний угловой скорости вращения ротора мотора ω_м. При этом, декремент затухания у новых моторов составляет величины 0,69 – 0,7, а у моторов с минимально допустимыми значениями к.п.д. (∆ŋ=20%), это значение 0,55-0,57. Как и в предыдущих вариантах большой разницы между МП-90 и МП-112 не проявляется.

Полученные выражения Т_м и d_м при подстановке в решения дифференциальных уравнений [10], позволят моделировать переходный процесс в гидромоторах МП-90 и МП-112. 

Выводы 

1. Выполнена параметрическая идентификация динамических моделей переходных процессов в поршневых гидронасосах НП-90,112 и поршневых гидромоторах МП-90,112. На основании метода анализа размерностей получены выражения для определения коэффициентов усиления К₁ и К₄ и постоянных времени Т₁ - Т₄, которые характеризуют переходный процесс в гидростатической трансмиссии ГСТ-90,112.

2. Установлено, что увеличение утечек в насосе и моторе приводит к увеличению их инерционности, что выражается в росте постоянных времени Т_н и Т_м, при одновременном уменьшении декрементов затухания колебаний рабочих параметров насоса d_н и мотора d_м, которые при подстановке в решения дифференциальных уравнений позволяют моделировать динамику переходного процесса в поршневых моторах и насосах.

1. Горбатов В.В. Аналіз впливу зміни технічного стану гідроприводів циклічної дії сільськогосподарських машин на перевитрати палива / Вісник ХНТУСГ, 2009. – Вип. 80. – с. 315-321.
2. ДСТУ 2193-93 Гідроприводи об’ємні. Насоси об’ємні та гідромотори. Загальні технічні умови.
3. Галин Д.А. Анализ технического состояния гидростатической трансмиссии ГСТ-90 // Повышение эффективности функционирования механических и энергетических систем: Сб. науч. тр. Морд. гос. ун-т им. Н.П. Огарева. – Саранск: «Красн. Окт.». – с. 117-120.
4. Башта Т.М. Техническая диагностика гидравлических приводов / Т.М. Башта, Т.В. Алексеева, В.Д. Бабанская. – М.: Машиностроение, 1989. – 264 с.
5. Колганов Е.В. Обґрунтування інформативних діагностичних параметрів технічного стану об’ємного гідроприводі трансмісії ГСТ-90 / Вісник Дніпропетровського ДАУ, 2009. – № 2. – с. 71-74.
6. Романюк А.М., Романюк Ю.А. Тепловой баланс гидроприводов и определение их основных теплотехнических параметров / Науковий вісник Національного аграрного університету. – Київ: - 2004, – вип. 73, частина 2, с. 274-278.
7. Романюк А.М., Романюк Ю.А. Анализ неустановившегося температурного режима гидропривода / Науковий вісник Національного аграрного університету. – Київ: - 2005, – вип. 80, частина 12, с. 227-232.
8. Кравчук В., Горбатов В. Методика визначення перевитрат палива при зміні технічного стану гідроприводів сільськогосподарських машин / MOTROL, LUBLIN, 2009, Vol. 11A, c. 239-246.
9. Динамика транспортно-тяговых колесных и гусеничных машин / Е.Е. Александров, Д.О. Волонцевич, В.А. Карпенко, А.Т. Лебедев, В.А. Перегон, В.Б. Самородов, А.Н. Туренко. – Харьков: Издательство ХГАДТУ (ХАДИ), 2001. – 642 с.
10. Бойко И. Г., Войтов А. В., Войтов В. А. Cтруктурная идентификация математической модели диагностирования объемного гидропривода ГСТ-90,112 / Науковий вісник Національного університету біоресурсів і природокористування України. Серія: техніка та енергетика АПК – К., 2015. – Вип. 226, с. 124-132.
11. Волков А.В., Мищенко В.Ю. Методика испытаний по определению компонентов объемных потерь блоков питания / Электронный журнал “Труды МАИ”, вып.69, c. 100-110.

Рецензии:

11.04.2016, 12:07 Смольникова Фарида Харисовна
Рецензия: В статье рассматривается вопросы выполнения параметрической идентификации математической модели диагностирования объемного гидропривода ГСТ-90,112. Статья написана согласно правилам журнала и рекомендуется к печати



Комментарии пользователей:

Оставить комментарий