МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ НЕЛИНЕЙНЫХ ДИНАМИЧЕСКИХ ЭФФЕКТОВ ПРИ МЕДЛЕННОМ ДВИЖЕНИИ С СУХИМ ТРЕНИЕМ
ПРИБОРЫ ТОЧНОЙ МЕХАНИКИ
УДК 621.01
Г. Б. ЗАМОРУЕВ, А. Л. ТКАЧЕВ
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ НЕЛИНЕЙНЫХ ДИНАМИЧЕСКИХ ЭФФЕКТОВ ПРИ МЕДЛЕННОМ ДВИЖЕНИИ С СУХИМ ТРЕНИЕМ
Предложена математическая модель трибологического взаимодействия, рассмотрены варианты взаимодействия объектов с сухим трением. Результаты расчетов приведены в графической форме.
Ключевые слова: эффект сухого трения, моделирование, расчеты, графическое представление результатов.
Известно, что коэффициент трения скольжения по своей природе не является постоян-
ным и зависит, в частности, от относительной скорости. При скорости, равной нулю, коэф-
фициент трения выше, чем при движении. Тело „прилипает“ к месту и для его сдвига тре-
буется большее, чем при равномерном движении, усилие. В момент начала движения коэф-
фициент трения практически мгновенно падает на некоторую величину от f0 до f00 . Далее происходит изменение (обычно уменьшение) коэффициента трения в зависимости от скоро-
сти относительного движения до более или менее стабильного значения f , которое и приво-
дят в таблицах в качестве коэффициента трения. Если скорость движения очень мала, неста-
бильность коэффициента трения (особенно его конечный мгновенный скачок в момент тро-
гания) приводит к скачкообразному характеру движения с остановками. Описанное явление
затрудняет тонкое позиционирование объекта, так как при попытке сдвинуть его на малое
расстояние происходит сначала накопление потенциальной энергии упругости связи, потом —
скачкообразная „разрядка“ с нерегулируемым смещением на некоторое расстояние, которое и
оказывается интервалом неопределенности. Эти особенности силы трения скольжения доста-
точно подробно на качественном уровне описаны Н. И. Колчиным в его работе [см. лит.].
Ситуация, описанная выше, проиллюстрирована на рис. 1. Объект с массой m располо-
жен на горизонтальной поверхности и находится под действием движущей силы и силы тре-
ния. Усилие от кинематического привода, поддерживающего постоянную горизонтальную
скорость v , передается объекту посредством упругой связи с коэффициентом c . Сила трения
является функцией скорости, т.е. Fтр = f (x) . Коэффициент трения при этом — переменная
величина — f (x) ≡ f (v) .
Для моделирования зависимости f (x) выберем гиперболоидальную зависимость в виде
f
( x )
=
x
d +
a
+ b,
f
′( x )
=
−
d (x + a)2
.
(1)
Коэффициенты a, b и d найдем из следующих условий:
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 5
Математическое моделирование нелинейных динамических эффектов
— при x = 0 коэффициент трения равен f00 ; — при некоторой „установившейся“ скорости vуст коэффициент трения равен f ; — производная f ′(x) в момент трогания равна некоторой величине f ′(0) .
v = const
с
37
Fтр x
0 mg x vt
Рис. 1
Представим систему уравнений для нахождения коэффициентов в выражениях (1):
f00
=
d a
+
b,
f
=
d vуст +
a
+
b,
f
′(0)
=
−
d a2
.
(2)
Окончательный вид зависимость приобретает в момент трогания, когда значение коэф-
фициента трения мгновенно падает от f0 до f00 .
f
( x )
=
⎛ ⎝⎜
x
d +
a
+
b
⎞ ⎟⎠
(x
>
0)
+
f0
( x
≤
0)
.
(3)
Зависимость, представленная выражением (3), приведена на рис. 2 (1 — f, 2 — f (v) ).
f, f(v)
0,13
0,12
0,11 2
0,1 1
0,09 0
0,5
1
1,5 v, м/с
Рис. 2
Рассмотрим силы, действующие на объект. Роль движущей силы выполняет упругая связь.
Действие упругой связи в данной модели будем считать „односторонним“, т.е. связь формирует
силу, только когда она „растянута“ ( vt > x ). Выражение для упругой связи будет следующим:
c (v t − x) (v t > x) .
(4)
Действие силы трения в различных ситуациях:
— Fтр = 0 , если vt ≤ x и x = 0 — связь „провисает“, объект неподвижен;
— Fтр = c (vt − x) , если vt > x , x = 0 и c (vt − x) < m g f0 — объект неподвижен, связь
„натянута“, но сила трения покоя не достигнута;
— Fтр = m g f0 — в момент, когда связь „натянулась“ до значения силы трения покоя;
— Fтр = m g f (x) — в процессе движения объекта ( x > 0 ).
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 5
38 Г. Б. Заморуев, А. Л. Ткачев
Исходя из указанных предпосылок составим динамическую модель в форме дифферен-
циального уравнения второго порядка:
mx − c(vt − x)(vt ≥ x) + c(vt − x)[x ≤ 0, 001 ∧ vt ≥ x ∧ c(vt − x) < mgf0 ] +
+mgf (x)(x > 0) = 0, x(0) = 0, x(0) = 0.
(5)
Интегрирование выражения (5) затруднено вследствие сильной нелинейности, особенно
в момент перехода от движения объекта к неподвижности под действием силы трения и при
отрицательном ускорении. С помощью численных методов интегрирования сложно опреде-
лить момент, когда скорость и ускорение становятся равными нулю. При отрицательном ус-
корении скорость переходит через нуль, только затем ускорение принимается равным нулю,
таким образом, на участке неподвижности объекта может сохраняться небольшая (а иногда и
значительная) отрицательная скорость. Результат интегрирования выражения (5) приведен на
рис. 3, а. На втором неподвижном участке можно заметить очень малый обратный наклон ли-
нии графика, свидетельствующий об описанном явлении.
а) x(t), м vt, м
v(t), м/с
v, м/с
0,3 2
0,2 1
0,1
0
0 б) x(t), м
vt, м v(t), м/с v, м/с
0,3
4
34 6 8 10 t, с
2
0,2 1
3 0,1
4 0
05
10 15 t, с
Рис. 3
Кривые 2 и 4 на рис. 3 отображают перемещение vt и скорость v кинематического при-
вода соответственно. Пересечение графиков перемещения объекта и кинематического привода на рис. 3, а объясняется тем, что в какой-то момент объект „обгоняет“ привод, и участки графика, лежащие выше линии движения привода, соответствуют ненапряженной упругой связи. Последнее явление хорошо иллюстрирует график силы (F0), развиваемой упругой связью и показанный на рис. 4, в кривой 1. На рис. 3 скорость объекта и его кинематического привода приведены кривыми 1 и 3. Графики ускорения объекта приведены на рис. 4, а, б.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 5
Математическое моделирование нелинейных динамических эффектов
39
Для сглаживания описанных выше трудностей возможно использовать другие варианты
модели рассматриваемого явления. Так, можно редуцировать порядок дифференциального
уравнения и построить модель в виде системы из двух уравнений первого порядка. Если обо-
значить x = x1, то система будет иметь следующий вид:
x − x1 = 0,
m x1 − c (vt − x) (vt ≥ x) +
+ c (vt − x)[ x1 < 0, 001∧ vt ≥ x ∧ c (vt − x) < m g f0 ] +
(6)
+m g f (x1) (x1 > 0, 001) = 0,
x(0) = 0, x1(0) = 0.
Проинтегрировав систему (6), получим значения перемещения и скорости объекта.
а) а(t), м/с2 2
0
–2 0
б) а, м/с2 1
0
5
10 15 t, с
–1 0
в) F0, Н Fi, Н 0,2
0,1
5 10
1 2
15 t, с
0
05
10 15 t, с
Рис. 4
Также может быть использован конечно-разностный метод. В данной задаче можно использовать принцип импульс силы (количество движения). Введем обозначения: число уча-
стков для расчета — n ; время моделирования — T ; интервал времени ∆t = T n ; перемен-
ная — счетчик цикла — i = 0, ..., n ; начальные условия следующие x0 = 0, v0 = 0, a0 = 0 ( v и
a — скорость и ускорение объекта). Суммарная сила на i-м интервале:
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 5
40 Г. Б. Заморуев, А. Л. Ткачев
Fi = c (vti − xi ) (vti ≥ xi ) −
− c (vti − xi )[ vvi < 0, 001 ∧ vti ≥ xi ∧ c (vti − xi ) < m g f0 ] −
−m g fi (vi > 0, 001).
(7)
Скорость объекта на интервале i+1:
vi+1
=
vi
+
Fi m
∆t
.
(8)
Перемещение и ускорение объекта на интервале i+1:
xi+1 = xi + vi ∆t,
ai+1
=
vi+1 − vi ∆t
.
(9)
Все три варианта динамической модели (5)—(9) фактически являются расчетными фор-
мулами и используются для расчетов при различных исходных параметрах в среде математи-
ческого моделирования MathCAD. В частности, в приведенном на рис. 3, б варианте разница
между коэффициентами трения покоя f0 и движения f не так велика, как в случае, показан-
ном на рис. 3, а. Остановки объекта происходят при „напряженной“ упругой связи, и коорди-
наты объекта постоянно смещены относительно уровня перемещения кинематического при-
вода. На рис. 4, в представлен график состояния упругой связи (кривая 2), соответствующего
режиму движения (см. рис. 3, б). Ускорение объекта при этом режиме приведено на рис. 4, б.
В заключение следует отметить, что все представленные варианты модели дают одина-
ковые числовые результаты. Модель подтвердила высокую нелинейность описываемого яв-
ления, особенно проявляющуюся на чрезвычайно коротком временном промежутке перехода
объекта от движения к покою. Решить указанную проблему удается за счет значительного
увеличения количества точек (участков) интегрирования. Число участков доводится до
1500—2500 при 10—20 с модельного времени, причем изменение числа участков даже на
единицу может сдвинуть границы интервала и изменить результат.
ЛИТЕРАТУРА
Колчин Н. И. Механика машин. Л.: Машиностроение, 1972. Т. 2.
Сведения об авторах Георгий Борисович Заморуев — канд. техн. наук, доцент; Санкт-Петербургский государственный универси-
тет информационных технологий, механики и оптики, кафедра мехатроники; E-mail: georgyz09@gmail.com Алексей Леонидович Ткачев — аспирант; Санкт-Петербургский государственный университет информационных технологий, механики и оптики, кафедра мехатроники
Рекомендована кафедрой мехатроники
Поступила в редакцию 25.12.09 г.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 5
УДК 621.01
Г. Б. ЗАМОРУЕВ, А. Л. ТКАЧЕВ
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ НЕЛИНЕЙНЫХ ДИНАМИЧЕСКИХ ЭФФЕКТОВ ПРИ МЕДЛЕННОМ ДВИЖЕНИИ С СУХИМ ТРЕНИЕМ
Предложена математическая модель трибологического взаимодействия, рассмотрены варианты взаимодействия объектов с сухим трением. Результаты расчетов приведены в графической форме.
Ключевые слова: эффект сухого трения, моделирование, расчеты, графическое представление результатов.
Известно, что коэффициент трения скольжения по своей природе не является постоян-
ным и зависит, в частности, от относительной скорости. При скорости, равной нулю, коэф-
фициент трения выше, чем при движении. Тело „прилипает“ к месту и для его сдвига тре-
буется большее, чем при равномерном движении, усилие. В момент начала движения коэф-
фициент трения практически мгновенно падает на некоторую величину от f0 до f00 . Далее происходит изменение (обычно уменьшение) коэффициента трения в зависимости от скоро-
сти относительного движения до более или менее стабильного значения f , которое и приво-
дят в таблицах в качестве коэффициента трения. Если скорость движения очень мала, неста-
бильность коэффициента трения (особенно его конечный мгновенный скачок в момент тро-
гания) приводит к скачкообразному характеру движения с остановками. Описанное явление
затрудняет тонкое позиционирование объекта, так как при попытке сдвинуть его на малое
расстояние происходит сначала накопление потенциальной энергии упругости связи, потом —
скачкообразная „разрядка“ с нерегулируемым смещением на некоторое расстояние, которое и
оказывается интервалом неопределенности. Эти особенности силы трения скольжения доста-
точно подробно на качественном уровне описаны Н. И. Колчиным в его работе [см. лит.].
Ситуация, описанная выше, проиллюстрирована на рис. 1. Объект с массой m располо-
жен на горизонтальной поверхности и находится под действием движущей силы и силы тре-
ния. Усилие от кинематического привода, поддерживающего постоянную горизонтальную
скорость v , передается объекту посредством упругой связи с коэффициентом c . Сила трения
является функцией скорости, т.е. Fтр = f (x) . Коэффициент трения при этом — переменная
величина — f (x) ≡ f (v) .
Для моделирования зависимости f (x) выберем гиперболоидальную зависимость в виде
f
( x )
=
x
d +
a
+ b,
f
′( x )
=
−
d (x + a)2
.
(1)
Коэффициенты a, b и d найдем из следующих условий:
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 5
Математическое моделирование нелинейных динамических эффектов
— при x = 0 коэффициент трения равен f00 ; — при некоторой „установившейся“ скорости vуст коэффициент трения равен f ; — производная f ′(x) в момент трогания равна некоторой величине f ′(0) .
v = const
с
37
Fтр x
0 mg x vt
Рис. 1
Представим систему уравнений для нахождения коэффициентов в выражениях (1):
f00
=
d a
+
b,
f
=
d vуст +
a
+
b,
f
′(0)
=
−
d a2
.
(2)
Окончательный вид зависимость приобретает в момент трогания, когда значение коэф-
фициента трения мгновенно падает от f0 до f00 .
f
( x )
=
⎛ ⎝⎜
x
d +
a
+
b
⎞ ⎟⎠
(x
>
0)
+
f0
( x
≤
0)
.
(3)
Зависимость, представленная выражением (3), приведена на рис. 2 (1 — f, 2 — f (v) ).
f, f(v)
0,13
0,12
0,11 2
0,1 1
0,09 0
0,5
1
1,5 v, м/с
Рис. 2
Рассмотрим силы, действующие на объект. Роль движущей силы выполняет упругая связь.
Действие упругой связи в данной модели будем считать „односторонним“, т.е. связь формирует
силу, только когда она „растянута“ ( vt > x ). Выражение для упругой связи будет следующим:
c (v t − x) (v t > x) .
(4)
Действие силы трения в различных ситуациях:
— Fтр = 0 , если vt ≤ x и x = 0 — связь „провисает“, объект неподвижен;
— Fтр = c (vt − x) , если vt > x , x = 0 и c (vt − x) < m g f0 — объект неподвижен, связь
„натянута“, но сила трения покоя не достигнута;
— Fтр = m g f0 — в момент, когда связь „натянулась“ до значения силы трения покоя;
— Fтр = m g f (x) — в процессе движения объекта ( x > 0 ).
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 5
38 Г. Б. Заморуев, А. Л. Ткачев
Исходя из указанных предпосылок составим динамическую модель в форме дифферен-
циального уравнения второго порядка:
mx − c(vt − x)(vt ≥ x) + c(vt − x)[x ≤ 0, 001 ∧ vt ≥ x ∧ c(vt − x) < mgf0 ] +
+mgf (x)(x > 0) = 0, x(0) = 0, x(0) = 0.
(5)
Интегрирование выражения (5) затруднено вследствие сильной нелинейности, особенно
в момент перехода от движения объекта к неподвижности под действием силы трения и при
отрицательном ускорении. С помощью численных методов интегрирования сложно опреде-
лить момент, когда скорость и ускорение становятся равными нулю. При отрицательном ус-
корении скорость переходит через нуль, только затем ускорение принимается равным нулю,
таким образом, на участке неподвижности объекта может сохраняться небольшая (а иногда и
значительная) отрицательная скорость. Результат интегрирования выражения (5) приведен на
рис. 3, а. На втором неподвижном участке можно заметить очень малый обратный наклон ли-
нии графика, свидетельствующий об описанном явлении.
а) x(t), м vt, м
v(t), м/с
v, м/с
0,3 2
0,2 1
0,1
0
0 б) x(t), м
vt, м v(t), м/с v, м/с
0,3
4
34 6 8 10 t, с
2
0,2 1
3 0,1
4 0
05
10 15 t, с
Рис. 3
Кривые 2 и 4 на рис. 3 отображают перемещение vt и скорость v кинематического при-
вода соответственно. Пересечение графиков перемещения объекта и кинематического привода на рис. 3, а объясняется тем, что в какой-то момент объект „обгоняет“ привод, и участки графика, лежащие выше линии движения привода, соответствуют ненапряженной упругой связи. Последнее явление хорошо иллюстрирует график силы (F0), развиваемой упругой связью и показанный на рис. 4, в кривой 1. На рис. 3 скорость объекта и его кинематического привода приведены кривыми 1 и 3. Графики ускорения объекта приведены на рис. 4, а, б.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 5
Математическое моделирование нелинейных динамических эффектов
39
Для сглаживания описанных выше трудностей возможно использовать другие варианты
модели рассматриваемого явления. Так, можно редуцировать порядок дифференциального
уравнения и построить модель в виде системы из двух уравнений первого порядка. Если обо-
значить x = x1, то система будет иметь следующий вид:
x − x1 = 0,
m x1 − c (vt − x) (vt ≥ x) +
+ c (vt − x)[ x1 < 0, 001∧ vt ≥ x ∧ c (vt − x) < m g f0 ] +
(6)
+m g f (x1) (x1 > 0, 001) = 0,
x(0) = 0, x1(0) = 0.
Проинтегрировав систему (6), получим значения перемещения и скорости объекта.
а) а(t), м/с2 2
0
–2 0
б) а, м/с2 1
0
5
10 15 t, с
–1 0
в) F0, Н Fi, Н 0,2
0,1
5 10
1 2
15 t, с
0
05
10 15 t, с
Рис. 4
Также может быть использован конечно-разностный метод. В данной задаче можно использовать принцип импульс силы (количество движения). Введем обозначения: число уча-
стков для расчета — n ; время моделирования — T ; интервал времени ∆t = T n ; перемен-
ная — счетчик цикла — i = 0, ..., n ; начальные условия следующие x0 = 0, v0 = 0, a0 = 0 ( v и
a — скорость и ускорение объекта). Суммарная сила на i-м интервале:
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 5
40 Г. Б. Заморуев, А. Л. Ткачев
Fi = c (vti − xi ) (vti ≥ xi ) −
− c (vti − xi )[ vvi < 0, 001 ∧ vti ≥ xi ∧ c (vti − xi ) < m g f0 ] −
−m g fi (vi > 0, 001).
(7)
Скорость объекта на интервале i+1:
vi+1
=
vi
+
Fi m
∆t
.
(8)
Перемещение и ускорение объекта на интервале i+1:
xi+1 = xi + vi ∆t,
ai+1
=
vi+1 − vi ∆t
.
(9)
Все три варианта динамической модели (5)—(9) фактически являются расчетными фор-
мулами и используются для расчетов при различных исходных параметрах в среде математи-
ческого моделирования MathCAD. В частности, в приведенном на рис. 3, б варианте разница
между коэффициентами трения покоя f0 и движения f не так велика, как в случае, показан-
ном на рис. 3, а. Остановки объекта происходят при „напряженной“ упругой связи, и коорди-
наты объекта постоянно смещены относительно уровня перемещения кинематического при-
вода. На рис. 4, в представлен график состояния упругой связи (кривая 2), соответствующего
режиму движения (см. рис. 3, б). Ускорение объекта при этом режиме приведено на рис. 4, б.
В заключение следует отметить, что все представленные варианты модели дают одина-
ковые числовые результаты. Модель подтвердила высокую нелинейность описываемого яв-
ления, особенно проявляющуюся на чрезвычайно коротком временном промежутке перехода
объекта от движения к покою. Решить указанную проблему удается за счет значительного
увеличения количества точек (участков) интегрирования. Число участков доводится до
1500—2500 при 10—20 с модельного времени, причем изменение числа участков даже на
единицу может сдвинуть границы интервала и изменить результат.
ЛИТЕРАТУРА
Колчин Н. И. Механика машин. Л.: Машиностроение, 1972. Т. 2.
Сведения об авторах Георгий Борисович Заморуев — канд. техн. наук, доцент; Санкт-Петербургский государственный универси-
тет информационных технологий, механики и оптики, кафедра мехатроники; E-mail: georgyz09@gmail.com Алексей Леонидович Ткачев — аспирант; Санкт-Петербургский государственный университет информационных технологий, механики и оптики, кафедра мехатроники
Рекомендована кафедрой мехатроники
Поступила в редакцию 25.12.09 г.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 5