Научная статья на тему 'Стабилизация динамического состояния станка как основа решения задач повышения точности механической обработки деталей'

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

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

Аннотация научной статьи по механике и машиностроению, автор научной работы — Бржозовский Борис Максович, Мартынов Владимир Васильевич, Янкин Игорь Николаевич, Бровкова Марина Борисовна

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

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

Похожие темы научных работ по механике и машиностроению , автор научной работы — Бржозовский Борис Максович, Мартынов Владимир Васильевич, Янкин Игорь Николаевич, Бровкова Марина Борисовна

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

Results of application of the dynamic processes simulation techniques in the metal-cutting machine tools for it stabilization, status control and increase of the machine working accuracy are presented in this article.

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

УДК 621.9.025.01

Б.М. Бржозовский, В.В. Мартынов, И.Н. Янкин, М.Б. Бровкова СТАБИЛИЗАЦИЯ ДИНАМИЧЕСКОГО СОСТОЯНИЯ СТАНКА КАК ОСНОВА РЕШЕНИЯ ЗАДАЧ ПОВЫШЕНИЯ ТОЧНОСТИ МЕХАНИЧЕСКОЙ ОБРАБОТКИ ДЕТАЛЕЙ

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

B.M. Brzhozovsky, V.V. Martynov, I.N. Yankin, M.B. Brovkova MACHINE TOOL DYNAMIC CONDITION STABILIZATION AS A BASIS OF INCREASE ACCURACY PROBLEMS DECISION OF DETAILS MACHINING

Results of application of the dynamic processes simulation techniques in the metal-cutting machine tools for it stabilization, status control and increase of the machine working accuracy are presented in this article.

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

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

Имитационная модель динамической системы токарного станка, разработанная в среде моделирования 81шиНпк, интегрированного в состав математического пакета МайаЬ, позволяет представить станок как систему взаимосвязанных функциональных модулей в виде упругой подсистемы, блоков задания частоты вращения (скорости резания) и подачи, процесса резания, резца и детали, имеющих блочную внутреннюю структуру (рис. 1), и отличается от разработанных ранее моделей [3, 4] тем, что:

- кроме нелинейностей «зазор» учитывает нелинейности типа «сухое трение» (в блоке «процесс резания»);

- отображает одновременное прохождение вибросигнала как в прямом (от детали к резцу), так и в обратном (от резца к детали) направлениях;

- обладает усовершенствованными (в направлении реализации в процессе моделирования динамического изменения параметров) блоками «зазор»;

- учитывает неравномерность частоты вращения, подачи инструмента и неоднородность материала заготовки;

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

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

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

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

вибросигнала при его прохождении через динамическую систему станка. Тем, следовательно, меньше будет значение НФУ.

Рис. 1. Пример структуры блоков модели: блок «упругая подсистема».

Блоки «Unit Delay», «Buffered-FFT» и «Frame scope» используются, соответственно, для задержки на один квант времени моделирования, получения и просмотра частотной и временной информации

х 10

□ 0.2 0.4 0.Б 0.8 1 1.2 1.4 1.Б 1.8 2

х 104

д

Рис. 2. Реакция НФУ на прохождение вибросигнала через модель: а - до объекта «зазор» между блоками «резец» и «упругая подсистема»; б - после объекта «зазор» между блоками «резец» и «упругая подсистема»; в - после объекта «револьверная головка» в блоке «упругая подсистема»; г - после объекта «суппорт» в блоке «упругая подсистема»; д - после объекта «шпиндель» в блоке «упругая подсистема»; е - средние значения НФУ в перечисленных точках

б

а

г

е

Рис. 3. Пример практического использования НФУ: а - реализация вибросигнала при работе и поломке отрезного резца из сплава Т15К6;

б - результаты мониторинга поломки

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

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

Другой метод стабилизации динамического состояния станков заключается в целенаправленной коррекции их динамических свойств на основе учета динамического взаимодействия формообразующих подсистем станка при резании.

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

Рассмотрим применение данного метода в задаче обеспечения качества обработки при внутреннем шлифовании.

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

Особенностью обработки материалов шлифованием является существенная роль процесса трения в зоне резания. На этой основе построена динамическая модель

процесса шлифования (рис. 4). Относительное движение поверхностей инструмента и изделия в направлении скорости резания вызывает их взаимные смещения в нормальном направлении. В динамической системе станка можно выделить две подсистемы - подсистему инструмента и подсистему изделия. Каждая из подсистем представлена двумя ортогональными парциальными системами, которые отражают некоторые обобщенные параметры динамической системы: mi - приведенные к зоне резания массы; pэi - эквивалентные жесткости групп вдоль соответствующих осей координат с учетом жесткости в области контакта; cэi - эквивалентные коэффициенты диссипативных сил; Qi - обобщенные силы; Xi - обобщенные оси координат.

Рис. 4. Модель динамической системы шлифования: индексы 1, 3 - нормальные парциальные системы; индексы 2, 4 - тангенциальные парциальные системы

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

dT дT dH дT

---=--------+---- , I

dV дH № дV

где обозначено: dT - приращение силы сопротивления движению двух тел; dV - приращение относительной скорости скольжения; dH - приращение нормального сближения поверхностей.

С учетом зависимости (1) определим возбуждающие силы, действующие на площадке трения, в следующем виде: Qn = cэ (х2 + х4) - подъемная (нормальная) сила и

Qt = pэ(л1 + х3) - кулонова (тангенциальная) сила, где cэ и pэ - эквивалентные

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

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

т1Х1 + сэ1Х&1 + рэ1х1 + с(1 - д1л12)( Х2 + Х4) = 0,

т2х2 + Сэ2Х2 + Рэ2х2 + Р(1 - Д2х2) (Х1 + Х3) = 0, (2)

т3Х3 + Сэ3Х3 + РэЗХ3 + С(1 - ^3Х1)(Х2 + •Х4) = 0

т4 Х4 + сэ4 Х4 + рэ4 Х4 + р(1 - д 4 Х^)( Х1 + Х3) = 0, где тг Хг - инерционные силы; ргХг - упругие силы; сгХг - диссипативные силы; с(1 -^гХг2)(Х2 + Х4) - подъемные силы; р (1 -дгХг2)(Х1 + Х3) - тангенциальные силы; -

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

Приведем систему (2) к стандартному виду:

Х1 + ШШ = (®2 - ®21)Х1 - §1Х1 - У1(1 - ^1Х12)(Х2 + Х4) = 0,

&&2 + ®()2 = (©? - © 02 )Х2 - §2Х2 - У2(1 - ^2Х2)(Х1 + Х3) = 0 ,

Х3 + ©23 = (®2 -®23) Х3 -53Х3 -У3(1 -Д3Х32)(Х2 + Х4) = 0, (3)

Х4 + ®24 = (©2 -®о4)Х4 -54Х4 -у4(1 -д4Х^)(Х1 + Х3) = 0,

где ш0г - парциальные частоты; 5г - коэффициенты диссипации; уг - коэффициенты связи парциальных систем через зону контакта.

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

Х = Е (А С08(Ш / -ф; ) , (4)

где А; - амплитуды колебаний в г-й парциальной системе на частоте ш; фу - фазы колебаний (г = 1...4; ] = 1...4).

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

4У сое (ш/ - ф] ) + ф]А вт (ш/ -ф] ) = 0 . (5)

Применение решения в виде (4) к системе (3) с учетом соотношений (5) дает

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

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

Вопрос о статусе подсистем решается из анализа устойчивости движений. Гипотетически устойчивые режимы колебаний возможны на каждой из четырех частот Ш1...Ш4. Решение в виде укороченных уравнений Ван-дер-Поля описывает переходной процесс в динамической системе, который может возникнуть при произвольном внешнем

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

Условия устойчивости колебательного движения на частоте изгиба оправки с кругом имеют вид:

д\{®] -Ш02) + ^2(ш/ -®(м)

)(1) ^

] (\ + ^\)(\ + ^4) (\ + ^2)(\ + СТз)

где коэффициенты о\.. ,о4 определяют взаимовлияние подсистем и имеют вид:

_ _ С08(ф\ -Ф4) ; _ _ 5Ш(ф\ -Ф4) ; _ т

°\ _ Л42 / ч 5 °2 _ Л42 • / ч 5 42 Шт

С0в(ф\-Ф2^ 2 42 вШ (ф\-Ф2^ 42 Т^((0] -®024 )2 + шХ

(Ш2 -Ш02 ) +®22^22 .

_ _ 51п(фз -Ф2) . _ _ С08(фз-Ф2) . _ т

°3 _ Л3\ • , ч ; °4 _ Л3\ , ч ; л3\ _ тп-

(ш2 -®2\)2 + ш12512

7)

^п (Ф\-Ф2^ С0§(Ф\-Ф2^ Пу (ш ( -ш23 )2 + ш|5|

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

(ш°-ш°\)(ш° ш 02) (\ + °\)(\ + ^3)-5\52ш5 (\ + ^2)(\ + °4) _ 0 , (8)

где фазы фг- на частоте Ш/ должны отвечать условиям:

(ш° - ш2\) [ (ф\ - Ф2) + к42 С08 (Ф\ - Ф4)] - °\ш/ (Ф\ - Ф2) + к42 ^ (Ф\ - Ф4)] _ 0 ,

(ш° - ш02) [п(Ф\ - Ф2) + к3\ ^ (Ф3 - Ф2)] - °2Ш/ [с0»(ф\ - Ф2) + к3\ С^(Ф3 - Ф2)] _ 0 ,

(ш2 - Ш03) [[(Ф3 - Ф2) + к42 С0Б (Ф3 - Ф4)] - ^3ш / [8Ш(Ф3 - ф2) + к42 ^^3 - Ф4)] _ 0 , (9)

(ш2/ -ш°;4)[п(ф\ -ф4) + к31Бт(ф3 -ф4)]-о4ш/ [с0б(ф\ -ф4) + к31 С0Б(ф3 -ф4)] _ 0 .

Остановимся подробнее на выражении (6). Здесь в числителях правой части неравенства сгруппированы показатели инерционных, упругих и диссипативных свойств подсистемы инструмента. Коэффициенты а\...а4 отражают связанности между подсистемами инструмента и изделия. Таким образом, в правой части указанного неравенства сосредоточены параметры подсистем инструмента и изделия. Их сочетание характеризует колебательные свойства формообразующих механических систем станка.

Левая часть неравенства (6) имеет вид (у\у2)(\ _ ср/(т\т2). Здесь произведение коэффициентов активных сил ср отражает условия контактирования инструмента с изделием в зоне резания, то есть является «мерой активности» процесса резания.

Таким образом, левая часть данного неравенства отражает действительный уровень возбуждения в подсистеме инструмента, а его правая часть отражает минимальнонеобходимый уровень, при котором в ней может возникнуть устойчивый автоколебательный режим. Обозначим этот уровень, соответствующий границе устойчивости движений на рассматриваемой частоте, через (у\у2)г\). Следует заметить, что это состояние равновесия может быть описано квазилинейной системой, в которую обращается система (2) при дгАу=0.

При уровне возбуждения, превышающем значение (у\у2)^, вступают в действие ограничивающие факторы в виде нелинейных зависимостей возбуждающих сил от амплитуд колебательных смещений в виде \ - А/ и колебательная система

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

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

При исследовании устойчивости практический интерес представляют связи между подсистемами, которые реализуются через зону резания. На рис. 5 приведен пример расчета, показывающий влияние связанности между подсистемами инструмента и изделия на границу устойчивости движений. В качестве варьируемого параметра используется коэффициент отношения парциальных частот ш04 / ©01. Кривая на рис. 5 представляет

собой геометрическое место точек, соответствующих значениям границы устойчивости движений на частоте изгиба оправки с кругом. Расчет выполнен для динамической системы со следующими показателями: парциальные частоты р1=1; р2=2; р3=0,1; р4=0,2...2; коэффициенты динамичности рг=0,01. Параметры указаны в универсальной безразмерной форме путем приведения их к параметрам первой парциальной системы.

Область значений показателя устойчивости (у1у2)Гр , расположенная выше кривой,

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

Вид кривой графика показывает, что изменение связанности двух подсистем оказывает заметное влияние на положение границы устойчивости. В частности, при ©04 / ©01 ~ 1 граница устойчивости снижается до минимальной величины. В этом состоянии система обладает наибольшей склонностью к возбуждению колебаний. В районе ©04 / ©01 ~ 0,8 граница устойчивых движений существенно сужается. В этой области система

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

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

(ШСгрЮ'г

Граница

устойчиво стй

1

} 1

с бл аст ь у ;тойчивых

ое ла :ть [В и жений

неустойчивости Ч *4

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8

Ш01

Рис. 5. Влияние связанности между подсистемами инструмента

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

ЛИТЕРАТУРА

1. Базров Б.М. Устранение автоколебаний при токарной обработке с помощью самоприспособляющихся систем управления / Б.М. Базров, В.И. Горюшкин // Станки и инструмент. 1977. № 4. С. 3-6.

2. Городецкий Ю.И. Создание математических моделей сложных автоколебательных систем в станкостроении / Ю.И. Городецкий // Автоматизация проектирования: сб. статей; под ред. В. А. Трапезникова. М.: Машиностроение, 1986. Вып. 1. С. 203-220.

3. Бржозовский Б. М. Исследование преобразующих свойств динамических систем металлорежущих станков методом математического моделирования / Б.М. Бржозовский, В.В. Мартынов, А.Н. Карпов // Информационные технологии в проектировании и производстве. 1998. № 3. С. 46-50.

4. Бржозовский Б.М. Модель упругой системы станка как основа для разработки динамического мониторинга точности обработки деталей шлифованием / Б.М. Бржозовский, В.В. Мартынов, А.Н. Ворыпаев // Процессы абразивной обработки, абразивные инструменты и материалы: сб. тр. Междунар. науч.-техн. конф. Волжский, 2001. С. 273-277.

5. Толстой Д.Н. К вопросу о роли нормальных перемещений при внешнем трении / Д. Н. Толстой, Р. Л. Каплан // Новое в теории трения: сб. статей. М.: Машиностроение, 1966. С. 76-83.

6. Моисеев Н.Н. Асимптотические методы нелинейной механики / Н.Н. Моисеев. М.: Наука, 1981. 400 с.

Бржозовский Борис Максович -

доктор технических наук, профессор,

заведующий кафедрой «Конструирование и компьютерное моделирование технологического оборудования машино- и приборостроения»

Саратовского государственного технического университета

Мартынов Владимир Васильевич -

доктор технических наук, профессор кафедры «Конструирование и компьютерное моделирование технологического оборудования машино- и приборостроения» Саратовского государственного технического университета

Янкин Игорь Николаевич -

доктор технических наук, профессор кафедры «Конструирование и компьютерное моделирование технологического оборудования машино- и приборостроения» Саратовского государственного технического университета

Бровкова Марина Борисовна -

кандидат технических наук, доцент кафедры «Программное обеспечение вычислительной техники и вычислительных систем»

Саратовского государственного технического университета

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