CC BY
25
УДК 621.95.08:51-74
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ И ПАРАМЕТРИЧЕСКАЯ ИДЕНТИФИКАЦИЯ ДИНАМИЧЕСКИХ СВОЙСТВ ПОДСИСТЕМ ИНСТРУМЕНТА И ЗАГОТОВКИ ПРИ ТОЧЕНИИ
© 2011 г. В.Л. Заковоротный, Фам Динь Тунг, Нгуен Суан Тьем
Донской государственный технический Don State Technical
университет University
Рассматривается математическое моделирование динамики упругих деформационных смещений вершины режущего инструмента и обрабатываемой заготовки в точке контакта с ней вершины режущего инструмента. Приводятся алгоритмы и результаты идентификации параметров модели для случая продольного точения
Ключевые слова: металлорежущий станок; математическое моделирование; процесс резания; динамические подсистемы.
The mathematical modeling of the dynamics of elastic deformation displacements of the corner of the cutting tool and workpiece at the contact point of the top of the cutting tool with workpiece is considered. The algorithms and the results of the model parameters for the case of linear turning are resulted.
Keywords: metal cutting machine tools; mathematical modeling; process turning; dynamic subsystems.
Введение. Одна из проблем, стоящих при изучении динамики процесса резания, связана с моделированием деформационных свойств подсистем режущего инструмента и заготовки [1—4]. В работах [1, 2] эти модели рассматриваются в скалярном представлении, в работах [3, 4] — в плоскости, нормальной к оси вращения заготовки, без раскрытия механизмов формирования деформационных смещений. В этом случае не удаётся изучить многие важные механизмы потери устойчивости равновесия, а также выяснить особенности преобразования траекторий движения исполнительных элементов станка в траектории формообразующих движений инструмента относительно заготовки. Статья посвящена построению пространственных динамических моделей и их параметрической идентификации, что существенно дополняет известные свойства преобразующей системы процесса резания.
Обоснование математической модели и её свойства. Для анализа динамики процесса резания используются пространственные конечномерные модели, приводящие к необходимости анализа следующего дифференциального уравнения:
—2 X —Х
т(Х,Бр,гр) —рр + ЫХ^р^р)^ + с(Х,БР,гР)Х = ¿г ¿г
= F(X,Sp,tp) + т,
где
F = {ад Sp, гР), ад Sp, гр),..., ад Sp, гР )}Т —
вектор-функция динамической характеристики процесса резания, раскрывающая зависимость сил резания от упругих деформационных смещений инструмента и заготовки, а также от технологических режимов: величины подачи на оборот Sp и глубины резания гр при заданной скорости; X = {Х1, Х2, Х3, Х4, Х5, Х6}Т — вектор упругих деформационных смещений вершины инструмента (первые три координаты) и заготовки в точке контакта с ней режущего инструмента (последние три координаты);
I (г) = {т т т т т л(тТ — изменяющиеся во времени составляющие сил резания, не объяснимые в координатах упругих деформационных смещений, которые интерпретируются как шум; т(Х, Sp, гР) = [ш5,к (Х, Sp, гР)], к(Х, Бр, гр) = [н5,к (Х, Бр, гр)],
с^, Sp, гР) = |са (X, Sp, гР)], к = 1,2,...6 — соответственно функциональные матрицы инерционных и диссипативных коэффициентов, а также функциональная матрица формирования упругой составляющей сил в зависимости от вектора деформационных смещений и технологических режимов. Во всех случаях здесь и ниже символ {...}т есть операция транспонирования.
Вначале проанализируем свойства матрицы с(Х, Sp, (р) = [с^к (X, Sp, (р)]. Рассмотрим подсистему режущего инструмента. Моделирование упругих деформационных смещений в линеаризованном виде, т. е. в виде матрицы жёсткости, позволяет представить деформационные смещения в новой системе координат {у^72,Уз} е Y (рис. 1), которая может быть получена на основе вращения системы координат {X!,Х2,Х3} е X с помощью углов Эйлера. В пространстве У матрица жёсткости является диагональной. Для перехода от координат пространства X к координатам Y выполним три последовательных положительных вращений на углы Эйлера ф, 0 , у . В результате получим матрицу преобразования (или матрицу перехода) в виде
пространственного случая представляет большие сложности, так как полученные недиагональные элементы матрицы с^ ^ являются сложными функциями от синусов и косинусов углов Эйлера.
Л =
cos ф cos у - sin ф cos 0 sin у - cos ф sin у - sin ф cos 0 cos у sin Ф sin 0 sin Ф cos у + cos Ф cos 0 sin у - sin ф sin у + cos Ф cos 0 cos у - cos Ф sin 0
Рис. 1. Система координат, в которой отчитывается упругое деформационное смещение и внешние силы
Подчеркнём, что матрица жёсткости с в статике является вещественной, положительной и
симметричной. Отметим два фундаментальных
sin 0 sin у
sin 0 cos у
0
Отметим некоторые свойства матрицы Л :
1. Т = Л-1 — матрица перехода координат пространства Y к координатам пространства X;
2. ЛТЛ = Е — единичная матрица, поэтому матрица Л ортогональна.
Аналогичное свойство справедливо и для матрицы Т. Тогда
X = Лу; АГ = ЛАГ1 или у = ТХ ;
АГ ? = ТАГ. (4)
Здесь АГ, АГт — вариации сил. Следовательно, с« . у = АГ« ; с№ = Л-1сЛ . (5)
Если для малых деформационных смещений оси (У1, У2, Уз)т можно считать главными, то сила, имеющая направление, совпадающее с какой-то из осей (У1, У2, Уз)т , вызывает упругое деформационное смещение только в этом направлении. В этом случае матрица упругих деформационных смещений с1 является диагональной. Таким образом, необходимо выбрать углы Эйлера так, чтобы матрица с1 была диагональной. Это можно сделать на основе решения уравнений сЦу(ф, 0, у) = 0,1 ^ у . Однако эта процедура для
свойства симметричных матриц:
1. Собственные числа вещественной, положительной, симметричной матрицы положительны, вещественны.
2. Собственные векторы, соответствующие различным собственным числам, ортогональны.
Известно [5, 6], что для вещественной симметричной матрицы с существует такая ортогональная матрица Л , что ЛтсЛ = с^^ , где ^ ^ = diag (сХх, с%2, сХз) — диагональная матрица, элементами которой являются собственные числа матрицы с. Таким образом, можно привести симметричную матрицу жёсткости с к диагональной форме, и вычислить матрицу преобразования Л на основе определения её собственных чисел и собственных векторов. Матрица жёсткости с(т) в новой системе координат (Ух, 72, Уз)Т имеет диагональные элементы, являющиеся собственными числами матрицы с. Наибольшее значение имеет случай, когда матрица с состоит из собственных чисел. Тогда справедливо
det(c - ХЕ) = 0.
Для каждого собственного числа определим собственный вектор Я1 = (ягд,gi,2,gi,з)T , который удовлетворяет системе уравнений
(с -Х,-Е)^ = 0.
Собственные векторы £ , £ , £ , соответствующие различным собственным числам Я-2, ^з, взаимно ортогональны. Найдем по векторам ^ , £2 , £3 нормированные собственные векторы в1 ° 5 г = 1,2,3
-1 в2, в3 при условии (в1, в1) = 1,
е = g'
Тогда ортогональная матрица преобразования Л , столбцами которой являются нормированные собственные векторы вг,
Л = (е1, е2, е3) =
е1,1 e2,1 e3,1 e1,2 e2,2 e3,2 e1,3 e2,3 e3,3
(1)
Матрицу (1) можно считать матрицей угловых коэффициентов преобразования координат пространства X в Y. Теперь можно определить углы Эйлера
0 = arccos е3 3; ф = arcsin(
e3,1
д/1 - (ез,з)
О;
у = arcsin(
е1,3
y¡1 - (ез,з)2
(2)
Аналогичный алгоритм можно использовать и для определения деформационных свойств подсистемы обрабатываемой заготовки. В этом случае можно считать, что деформационные смещения заготовки в направлении скорости подачи отсутствуют. Тогда деформационную модель можно рассматривать в плоскости, нормальной к оси вращения заготовки. В этом случае эллипсоид жёсткости для пространственного случая преобразуется в эллипс жёсткости в плоскости (рис. 2). Тогда можно сразу получить выражение для ориентации эллипса жёсткости и значений его диагональных элементов. Если известны матрицы упругих деформационных смещений с в простран-
стве X, то их выражения
с(Л)
в пространстве Л
получаем на основе следующего преобразования:
лл) _ дТ са с _ ах,лсах,л
Для этой матрицы выполним преобразование поворота координат
(л) =
с4 ' =
cosa - an а sin а cosa
с1,1 c1,2
С2,1 с2,2
cosa - sin a
a
cosa
Потребуем, чтобы недиагональные элементы в матрице (3) были равны нулю
[(c2,2 - cj,j) sin 2a + 2c12 cos 2a] = 0 , тогда получаем
1 . / 2ci 2 . a = — arctg(--!--)
2 (c1,1 - c2,2)
Матрица является диагональной, т.е.
с(Л) = diag (со,1, со,2), (4)
где
со 1 = (c1^os2a + С2 2 sin2 a - 2c12 sin a cos a) ;
c0 2 = (C11 sin2 a + C2 2со^a + 2c2 sin a cos a).
Выполняя обратное преобразование координат {Y1, y2)t в x = {^1, ^2)^ , получаем выражение для определения матрицы деформационных смещений вершины инструмента при заданной матрице (4)
с =
cosa - sin a
a
cosa
c0,1 0
0
c0,2
a - a sin a cosa
(c01cos2a + c0,2 sin2 a) [2(c0,2 - C0,1)sin2a]
[2(c0,2 - C0,1)sin2a] (c01 sin2 a + c0 2 cos2 a)
При с01 _ с0,2 матрица является диагональной.
Этот случай характерен для подсистемы обрабатываемой заготовки. Для этого случая любая ортогональная система координат будет характеризоваться тем, что сила любого направления вызывает деформационные смещения этого же направления. Это естественно, так как подсистема заготовки обладает симметричными деформационными свойствами. Поэтому при моделировании деформационных свойств подсистемы заготовки можно использовать скалярное представление.
Изучение статических и динамических деформационных свойств подсистемы режущего инструмента, в том числе корреляционно- спектральный анализ реакций системы на внешние дельтообразные возмущения, показывает, что в частотном диапазоне, ограниченном первыми
формами колеба-
(c11cos2a + С2 2 sin2 a - 2q2 sin a cos a) [(c2 2 - C11) sin 2a + 2c12 cos 2a]
[(c2 2 - C11) sin 2a + 2c12 cos 2a] (c11 sin2 a + C2 2 cos2 a + 2q 2 sin a cos a)
(3
нии, справедливо следующее свойство. Колебательные реакции подсистемы инструмента, вызванные внешними дельтообразными возмущениями, как и статическими си-
лами, обладают свойством пространственной избирательности. Причём для динамических и для статических реакций направления осей коллине-арного направления совпадают. Поэтому для подсистемы инструмента справедлива схематизация, приведённая на рис. 2 в плоскости. Следовательно, справедливы следующие свойства преобразования координат и параметров матриц т, h, с: 1. ЛТсЛ = с(7), ЛThЛ = с(т), ЛттЛ = с(т). 2. Преобразование пространства X в Y или Y в X характеризуется линейной операцией. Следовательно, характеристические полиномы в этих пространствах будут одни и те же, т. е. А(р) = А^(р).
Параметрическая идентификация. В пространстве X реакции упругих смещений, скоростей и ускорений при неизменных по модулю силовых возбуждениях описывают эллипсоид, главные оси которого остаются одними и теми же. Поэтому идентификации подлежат диагональные элементы матриц т, h, с пространства Y, а также матриц угловых коэффициентов (1) или углов Эйлера (2). Так как матрицы угловых коэффициентов для т, h, с одни и те же, то идентификация выполняется в два этапа.
На первом этапе осуществляется идентификация с. Придавая силам следующие измеримые приращения
АГ(1) = {АГ1,0,0}т , АГ(2) = {0, АГ2,0}т , АГ(3) = {0,0, АГ3}Т
и соответственно измеряя вариации х(1) = {х(1) х(1) х(1)}т х(1) = {х(1) х(1) х(1)}т и
X — \Х"1 5 Х2 5 Х3 ) X — \Х"1 5 Х2 5 Х3 |
х® = {х{^, х2р, х3^}Т деформационных смещений, получаем системы для оценивания элементов матрицы )
с^*)х() = АГ(>, i = 1, 2, 3. (5)
Система (5) позволяет оценить параметры матрицы жесткости ). Отметим, что при проведении исследований деформации измеряются приборами, база которых находится на несущей системе станка.
Для фиксированного значения точки равновесия из (5) получаем массив значений сг-,у , i, у = 1,2,3, сг-,у = сл при i ф у , по которому получаем оценки параметров жесткости по математическому ожиданию и значения их дисперсии. Математические ожидания сг-,у являются оценками матрицы жёсткости с в системе координат 0X1X2X3 для заданной точки равновесия. Для них также вычисляются дисперсии дi, у. По полученным данным определяются матрицы угловых коэффициентов и диагональных элементов матрицы жёсткости, проведённые к нормальным координатам.
На втором этапе оцениваются параметры матриц т и h. При этом в качестве исходной используется информация об оценках весовых функций системы (колебательных реакциях системы
Рис. 2. Система координат, в которой отсчитывается упругое деформационное смещение вершины режущего инструмента и внешние силы в плоскости
на дельтообразные возмущения), и их частотных образах. При определении обобщённых масс в основу положены значения резонансных частот, близких при малых коэффициентах затухания к собственным частотам системы. В свою очередь, резонансные частоты определяются из условия минимума характеристического полинома системы А( р) = А^ (р) после подстановки р = уш и вычисления квадрата модуля, т. е.
А( ую)А(-уш) = А«( уш)А«(-уш).
При этом А(т)(р) = /7(с0/ + %р + т01р2). В основу оценивания обобщённых масс ту положены легко измеримые автоспектры колебательных реакций системы. Так как су определены на первом этапе проведения исследований, то обобщённые массы оцениваются по очевидным зависимостям т0г = сщ /(^0 г')2 , где — резонансная частота соответствующей консервативной системы.
Так как колебательные контуры добротные, то коэффициент затухания можно определить также на основе анализа уширения спектральной линии каждой из резонансных частот, т. е. ^ = (^ - o0":í>)("o,i)-1, i = 1,2, где (о0+/ -^0/)— ширина резонансной кривой на высоте амплитуд 0,707 от резонансной амплитуды. Тогда
^ = 2^ 0.0^ = 2^7т<ус07, i = 1,2.
Затем уточняются собственные частоты для неконсервативной системы
о(рез-) = (0.0,1 )[1 - )2]2, i = 1,2,
где £ — коэффициент затухания рассматриваемого колебательного контура. Коэффициенты затухания в механических системах металлорежущих станков обычно лежат в диапазоне (0,02—0,06). Для определения moJ не обязательно измерять возмущающие силы, так как значения обобщённых жёстко-стей с0,; являются известными. Частоты 00,1 и О0 2 определяются на основе автоспектров в безразмерном виде. Обратим внимание на правило, связывающее соотношения собственных частот системы со значениями диагональных элементов матрицы инерции. Если справедливо
(00,3)2 /(00Д)2 = с0,1/¿0,3
и (О0 3)2 /(О0 2)2 = с0,2 / с0,3, то Щц = const .
Тогда в пространстве X матрица инерции также диагональна, и все диагональные элементы рав-
ны = const. В этом случае массу подвески инструмента можно считать неизменной и независимой от матрицы угловых коэффициентов.
Главная сложность параметрической идентификации h и m заключается в необходимости оценивания параметров для различных точек равновесия системы. Для этого необходимо создавать внешнее силовое поле, смещающее точку равновесия (X3). Для этого введем в систему дополнительный упругий элемент, который вносит погрешности, прежде всего, при определении собственных частот. Главное требование к упругому элементу — линейность упругих свойств, высокая добротность изгибных колебаний и ориентация по направлению осей Yi, i = 1,2,3. Тем самым его диссипацией, по сравнению с диссипацией подсистемы режущего инструмента, можно пренебречь. Тогда введенная связь изменяет один из диагональных элементов жесткости, что, не изменяя общей методики, позволяет оценить изменения параметров системы в зависимости от точки равновесия.
Анализ особенностей динамических подсистем. Проанализируем особенности динамических подсистем инструмента и заготовки на примере станка 1К62. Ограничимся случаем, когда подсистема инструмента рассматривается в плоскости, нормальной к оси вращения заготовки. На рис. 3 представлена экспериментально определенная диаграмма деформационных смещений вершины инструмента при двух постоянных по модулю силах, имеющих различные в плоскости направления. На диаграмме выделены траектории деформационных смещений вершины инструмента (две жирных точечных кривых, соответствующих двум значениям модуля силы 100 кг и 200 кг). В данном случае вылет инструмента из зажимного приспособления равен 35 мм. Пунктирными прямыми показаны некоторые направления действия силы. Кроме этого рассматриваются направления внешней силы, соответствующие осям X1 (отжимающее направление смещений) и X2 (направление скорости резания). Стрелками показаны направления деформационных смещений, ортогональные направлению действия силы. На иллюстрации углам а = 19° и а = 109° соответствуют только коллинеарные деформационные смещения. Всем остальным направлениям соответствуют как коллинеар-ные, так и ортогональные деформационные смещения. Важно подчеркнуть, что направления только коллинеарных деформационных смещений не зависят от величины модуля силы в диапазоне его изменения (20,0 кг—200,0 кг).
X 2, мкм 70,0
60,0 50,0
40,0 30,0
20,0 10,0
20,0 _ 40,0 4.- -. v . . i
100,0 Xj, мкм
Рис. 3. Типичная диаграмма суммарных деформационных смещений вершины режущего инструмента в зависимости от направления и модуля внешней силы для суппортной группы станка 1К62
По мере увеличения внешних сил изменяются значения матриц жёсткости. При этом ориентация эллипса жёсткости практически не меняется. На рис. 4 показаны зависимости жёсткостей от {У1, у2}Т е Л. Они хорошо аппроксимируются экспоненциальными функциями, т. е.
MYi) = $ + c$[1 - exp(-*)]; C2,o(Yi) = 4Ю + c220 [1 - exp(- tY2))]
где с«, с20 — значения жёсткостей на начальном этапе приложения внешней нагрузки;
(2) (2)
С0'с2о _ значения приращения жесткостеи по мере увеличения деформационных смещений; Ту(1),Ту(2) — параметры, характеризующие скорость увеличения параметров жёсткости. Функции жёсткостей имеют размерность [кг/мм], а Т^ имеют размерность [мм].
Приведём также данные об основных параметрах, характеризующих деформационные свойства вершины режущего инструмента для различных станков токарной группы (таблица). Поперечное сечение инструмента определяется следующими раз-
мерами: ширина основания — 30 мм, высота — 26 мм. Длина приведена в таблице.
Данные таблицы показывают, что параметры, характеризующие деформационные смещения, зависят не только от типа станка, но и от параметров инструмента, в частности, от его вылета из зажимного приспособления. Не трудно показать, что при этом упругие деформационные смещения вершины режущего инструмента в направлении Х2 возрастают существенно больше, чем деформации стержня инструмента. Здесь приходится считаться с тем, что за счёт увеличения вылета возрастает момент вращения, поворачивающий всю суппортную группу, а также с перераспределением контактных напряжений и деформаций в точках сопряжения отдельных деталей, входящих в суппортную группу. Именно этим объясняется зависимость ориентации осей эллипса жёсткости в зависимости от внешних сил при больших вылетах инструмента из зажимного приспособления, задающих точку равновесия системы инструмента. Что касается матрицы демпфирования, то она обладает такими свойствами: ориентация осей коллинеарных направлений эллипса диссипации совпадает с ориентацией осей эллипса жёсткости.
с01, кг/мм • 103 4,5
3,0
1,5
0,0
0,0 20,0 40,0 60,0 80,0 Yi, мкм
с02,кг/мм • 103 6,0
4,5
3,0
1,5
0,0
0,0 10,0 20,0 30,0 40,0 У 2,мкм
Рис. 4. Зависимость элементов матриц жёсткости деформационных смещений вершины инструмента по направлениям ориентации эллипсов жёсткости в зависимости от деформационных смещений по этим же направлениям
Параметры матриц жёсткости деформационных смещений для различных станков и вылета инструмента из зажимного приспособления
Тип станка с0 1 [кг/мм] с0 2 [кг/мм] а0 В^глет инструмента, 1 [мм]
1К62 2,9 103 5,4 103 100 35
1К62 0,7 102 5,1 103 (100 -180) 120
УТ16Ф3 4,3 103 6,7 103 120 35
16К20 с УЧПУ №210 3,7 103 6,5 103 140 54
1В340Ф30 2,8 103 4,5 103 190 54
Однако коэффициент демпфирования по направлению большей жёсткости меньше, чем по направлению меньшей жёсткости. Кроме этого обнаружена всегда проявляющаяся тенденция. По мере увеличения смещения точки равновесия, зависящего от внешней силы, уменьшается дисси-пативное влияние суппортной группы на колебания. Здесь также возможна аппроксимация, аналогичная (6):
Мто = И,1) + л$ехр(- А);
МУ1) = А$ + А220)ехр(- Т22)),
Ту
где (/г(о + Иц^ЛИу) + й220)) — значения коэффициентов демпфирования на начальном этапе приложения внешней нагрузки; й^с?, ИЦ2) — значения приращения коэффициентов демпфирования по мере увеличения деформационных смещений; Тт(1),Ту(2) — параметры, характеризующие скорость уменьшения коэффициентов демпфирования. Коэффициенты демпфирования имеют размерность [кг/с • м^, а Т^ 1\Т^2) имеют размерность [мм].
Что касается подсистемы заготовки, то, как уже отмечено, её свойства представимы в рассматриваемом частотном диапазоне в виде обобщённой массы, подвешенной в пространстве на упруго-дисси-пативных подвесках, имеющих равные значения жёсткостей и коэффициентов демпфирования.
Упругие свойства подсистемы заготовки зависят от координаты перемещения точки контакта
инструмента относительно заготовки. При этом собственные частоты колебаний подсистемы заготовки остаются неизменными. Для согласования этого противоречия необходимо рассматривать дополнительно изменение приведённой массы вдоль траектории движения инструмента относительно заготовки (рис. 5).
Выводы. Приведённые математические модели упругих деформационных свойств подсистем режущего инструмента и обрабатываемой заготовки позволяют учитывать главные особенности динамики деформационных смещений вершины режущего инструмента и заготовки в точке контакта с ней режущего инструмента. При этом в подсистеме режущего инструмента в частотном диапазоне до 5,0 кГц можно выделить три ортогональных направления, деформации по которым имеют только коллинеарные составляющие. Для этих осей матрицы динамической жёсткости и демпфирования и являются диагональными. В любой другой ортогональной системе координат эти матрицы симметричны и положительно определённы. Элементы этих матриц зависят от матриц угловых коэффициентов или углов Эйлера. Для анализа деформационных смещений подсистемы заготовки достаточно иметь информацию о распределении скалярной величины жёсткости по длине, дополненную данными об изменении инерционного коэффициента.
Основные результаты в статье получены при финансовой поддержке РФФИ по проекту 07-09-90000.
Рис. 5. Изменения приведенных коэффициентов жёсткости и массы в зависимости от координаты вдоль оси вращения заготовки
Литература
4. Заковоротный В. Л., Флек М. Б. Динамика процесса резания. Синергетический подход. Ростов н/Д., 2006. 876 с.
1. Соколовский А. П. Научные основы технологии
машиностроения. М., 1955. 514 с.
5. Хорн Р., Джонсон Ч. Матричный анализ. М., 1989. 655 с.
2. Адаптивное управление станками / под ред. Б. С.
Балакшина. М., 1973. 688 с.
6. Гантмахер Ф. Р. Лекции по аналитической механике. М., 1966. 299 с.
3. Кудинов В. А. Динамика станков. М., 1967. 360 с.
Поступила в редакцию
22 ноября 2010 г.
Заковоротный Вилор Лаврентьевич — д-р техн. наук, профессор, Донской государственный технический университет. Тел. 2738510. E-mail: vzakovozotny@dstu.edu
Фам Динь Тунг — канд. техн. наук, докторант, Донской государственный технический университет. Тел. 89604540318. E-mail: phamdinhtung@mail.ru
Нгуен Суан Тьем — аспирант, Донской государственный технический университет. Тел. 89508558472. E-mail: cyclone_rus0309@yahoo.com
Zakovorotny Vilor Lavrentevich — Doctor of Technical Sciences, professor, Don State Technical University. Tel. 2738510. E-mail: vzakovozotny@dstu.edu
Pham Dinh Tung — Candidate of Technical Sciences, рost-doctoral student, Don State Technical University. Tel. 89604540318. E-mail: phamdinhtung@mail.ru
Nguyen Xuan Chiem — post-graduate student, Don State Technical University. Tel. 89508558472. E-mail: cyclone_rus0309@yahoo .com
