Введение
Данная статья являет собой логическое продолжение темы, начатой в предыдущей публикации . Как и было обещано в комментариях, рассмотрим применимость метода избыточных координат к динамическому анализу механических систем движущихся под действием сил сухого кулоновского трения. В качестве иллюстративного примера решим следующую задачу
Тонкий однородный стержень массы m = 2 кг, длины AB = 2l = 1 м в точке A шарнирно прикреплен к невесомому ползуну, перемещающемуся в горизонтальных шероховатых направляющих. В начальный момент времени стержень расположен вертикально, затем его отклоняют от вертикали на ничтожно малый угол и отпускают без начальной скорости. Необходимо составить уравнения движения данной механической системы и найти закон её движения. Коэффициент трения между ползуном и направляющими равен f = 0,1.
Прежде чем приступить к решению задачи предлагаемым автором методом, рассмотрим немножко элементарной теории, касающейся сухого трения.
1. Что может быть «проще» трения?
Нет более страшного наказания для механика, чем сила трения. Появляясь в задаче, эта сила сразу делает её существенно нелинейной, ибо ведет себя достаточно интересным образом.
Рассмотрим довольно простой пример. Пусть на шероховатой поверхности покоится горизонтальный брусок.
Пусть в начале к нему не прилагают никаких сил (кроме силы тяжести и нормальной реакции). В этом случае и сила трения между бруском и плоскостью будет равна нулю.
Теперь приложим к бруску небольшую горизонтальную силу. Брусок не сдвинется с места, так как в ответ на наше воздействие со стороны поверхности на него станет действовать сила трения, которая будет удовлетворять условию
Будем постепенно увеличивать силу 
и, согласно (1) расти будет и сила трения, которая в этом случае называется силой трения покоя. Так будет продолжаться до тех пор, пока сила трения покоя не достигнет величины
называемой предельной величиной силы трения покоя. Здесь f — коэффициент сухого трения между бруском и плоскостью; N — нормальная реакция со стороны плоскости. После этого сила трения расти перестанет, а при дальнейшем увеличении горизонтальной силы начнется скольжение бруска. Сила трения перейдет в силу трения скольжения, равную
где 
— скорость бруска.
Пример весьма тривиальный, однако он раскрывает суть поведения силы сухого трения. Таким образом, получаем следующий алгоритм расчета силы трения:
Если точка, где приложена сила трения неподвижна:
- Расчитываем силу трения покоя и нормальную реакцию
- Проверяем условие
при нарушении которого принимаем силу трения равной предельной силе трения покоя
Если точка приложения силы трения движется:
- Вычисляем нормальную реакцию
- Вычисляем силу трения скольжения, согласно выражению (2)
2. Моделирование движения системы с трением
Теперь решим нашу задачу. Рассматриваемая нами система имеет две степени свободы, однако из-за необходимости определения нормальной реакции расширяем число степеней свободы до трех и получаем следующую расчетную схему
Здесь в качестве обобщенных координат берем вектор
где x,y — координаты точки A; 
— угол наклона стержня к вертикали. Вооружаемся Maple'ом
# Чистим память
restart;
# Подключаем линейную алгебру
with(LinearAlgebra):
# Подключаем лагранжеву механику
read `/home/maisvendoo/work/maplelibs/mechanics/lagrange.m`:
Определяемся с кинематикой системы
# Обобщенные координаты системы
q := [x(t), y(t), phi(t)]:
# Координаты точки A
xA := q[1]:
yA := q[2]:
# Координаты центра масс стержня
xC := q[1] - L*sin(q[3]):
yC := q[2] + L*cos(q[3]):
# Радиус-векторы точек A и C
rA := Vector([xA, yA]):
rC := Vector([xC, yC]):
# Вычисляем скорость ценра масс стержня
VectorCalculus[BasisFormat](false):
vC := VectorCalculus[diff](rC, t):
Вычисляем её кинетическую энергию
# Момент инерции стержня относительно центра масс
J := m*(2*L)^2/12:
# Кинетическая энергия системы
T := simplify(J*diff(q[3], t)^2/2 + m*DotProduct(vC, vC, conjugate=false)/2);
Maple выдает такой результат
Довольно громоздко, но «ковырять» нам не руками на листочке, поэтому двигаемся дальше. Задаем векторы и точки приложения сил
# Векторы сил, приоложенных к системе
Mg := Vector([0, -m*g]): # Сила тяжести
F_A := Vector([-F, 0]): # Сила трения
N_A := Vector([0, N]): # Нормальная реакция
# Формируем массив сил
Fk := [Mg, F_A, N_A]:
# Формируем массив точек приложения сил
rk := [rC, rA, rA]:
Получаем уравнения движения системы в форме Лагранжа 2 рода
# Составляем уравнения Лагранжа 2 рода
EQs := LagrangeEQs(T, q, rk, Fk):
Получаем трёх «крокодилов»
Эти уравнения пришлось вбить в статью руками, ибо «копипаста» LaTeX-вывода Maple приводит к неприглядному виду результата. Но даже так видно — уравнения сложны и с учетом того что F — это сила трения, аналитически не интегрируемы.
Теперь введем уравнения связей. Во-первых, ползун движется по горизонтальным направляющим, поэтому
Кроме того, в том случае когда ползун неподвижен, а сила трения покоя не превысила предельного значения, включается ещё одна связь
где 
— некоторая горизонтальная координата ползуна. Теперь преобразуем систему (4) — (6) с учетом уравнений (7) и (8) и найдем силу трения покоя и нормальную реакцию
# Уравнения связей
link_eq1 := q[1] = A: # Ползун неподвижен, если сила трения - сила трения покоя
link_eq2 := q[2] = 0: # Ползун движется вдоль оси X
Link_EQs := {link_eq1, link_eq2}:
# Трансформируем систему для поиска силы трения покоя и нормальной реакции
Reduced_EQs := ReduceSystem(EQs, Link_EQs, q):
solv_reduced := SolveAccelsReacts(Reduced_EQs, [q[3]], [F, N]):
# Выделяем силы из решения
for i from 1 to numelems(solv_reduced) do
if has(solv_reduced[i], F) then F1 := rhs(solv_reduced[i]); end if;
if has(solv_reduced[i], N) then N1 := rhs(solv_reduced[i]); end if;
end do:
Тут приведу результат непосредственно выданный Maple
Если посмотреть на полученные выражения, то они вполне соответствуют логике процесса. Теперь получим выражение для расчета нормальной реакции в случае, когда ползун скользит по направляющим, учитывая, что в этом случае сила трения будет определятся выражением
(минус уже имеется в уравнениях движения)
# Трансформируем систему для поиска нормальной реакции при скольжении ползуна
Slade_EQs := ReduceSystem(EQs, {link_eq2}, q):
# Силу трения принимаем равной силе трения скольжения
for i from 1 to numelems(Slade_EQs) do
Slade_EQs[i] := subs(F = f*N*signum(diff(q[1], t)), Slade_EQs[i]);
end do:
solv_slade := SolveAccelsReacts(Slade_EQs, [q[1], q[3]], [N]):
# Выделяем нормальную реакцию из решения
for i from 1 to numelems(solv_reduced) do
if has(solv_slade[i], N) then N2 := rhs(solv_slade[i]); end if;
end do:
М-да, «крокодилище» вышел ещё тот, особенно с учетом что Maple таки довольно избыточно генерирует LaTeX-код
Все необходимые нам выражения получены, теперь можно переходить к моделированию. В отличие от задачи с маятником, о которой я уже писал, тут мы честно трансформируем наши уравнения Maple-средствами для вида пригодного к численному решению. Прежде всего решим уравнения (4) — (6) относительно обобщенных ускорений
# Разрешаем основную систему уравнений относительно обобщенных ускорений
Main_EQs := SolveAccelsReacts(EQs, q, []):
# Число обобщенных координат
s := numelems(q):
Результат уже не буду приводить — он тоже довольно громоздкий. Перейдем к фазовым координатам
# Воводим фазовые координаты
# Y[1] -> x(t)
# Y[2] -> y(t)
# Y[3] -> phi(t)
# Y[4] -> vx(t) - горизонтальная проекция скорости ползуна
# Y[5] -> vy(t) - вертикальная проекция скорости ползуна
# Y[6] -> omega(t) - угловая скорость стержня
# Переходим к фазовым координатам в выражениях для реакций и трения покоя
for i from 1 to s do
N2 := subs(diff(q[i], t) = y[i+s], N2);
N2 := subs(q[i] = y[i], N2);
N1 := subs(diff(q[i], t) = y[i+s], N1);
N1 := subs(q[i] = y[i], N1);
F1 := subs(diff(q[i], t) = y[i+s], F1);
F1 := subs(q[i] = y[i], F1);
end do:
# Переходим к фазовым координатам в уравнения движения
for i from 1 to s do
for j from 1 to s do
eq := Main_EQs[j];
if has(eq, diff(diff(q[i], t), t)) then accel[i] := rhs(eq); end if;
end do;
for j from 1 to s do
accel[i] := subs(diff(q[j], t) = y[j+s], accel[i]);
accel[i] := subs(q[j] = y[j], accel[i]);
end do:
end do:
Формируем функции вычисления необходимых нам сил
# Вычисление нормальной реакции при покоящемся ползуне
GetN1 := proc(mass, length, grav_accel, fric_coeff, Y)
local react := subs(m = mass,
L = length,
g = grav_accel,
f = fric_coeff, N1);
local i;
for i from 1 to numelems(Y) do
react := subs(y[i] = Y[i], react);
end do:
return evalf(react);
end proc:
# Вычисление нормальной реакции при скользящем ползуне
GetN2 := proc(mass, length, grav_accel, fric_coeff, Y)
local react := subs(m = mass,
L = length,
g = grav_accel,
f = fric_coeff, N2);
local i;
for i from 1 to numelems(Y) do
react := subs(y[i] = Y[i], react);
end do:
return evalf(react);
end proc:
# Вычисление силы трения покоя
GetF1 := proc(mass, length, grav_accel, fric_coeff, Y)
local react := subs(m = mass,
L = length,
g = grav_accel,
f = fric_coeff, F1);
local i;
for i from 1 to numelems(Y) do
react := subs(y[i] = Y[i], react);
end do:
return evalf(react);
end proc:
Приведенный код хоть и объемный, но довольно прост — выполняется подстановка численных параметров в соответствующие выражения и их вычисление. Такую же функцию формируем и для вычисления обобщенных ускорений
# Вычисление вектора обобщенных ускорений
GetAccel := proc(mass, length, grav_accel, fric_force, normal_react, Y)
local acc;
local i, j;
for i from 1 to numelems(Y)/2 do
acc[i] := subs(m = mass,
L = length,
g = grav_accel,
F = fric_force,
N = normal_react, accel[i]);
end do:
for i from 1 to numelems(Y)/2 do
for j from 1 to numelems(Y) do
acc[i] := evalf(subs(y[j] = Y[j], acc[i]));
end do:
end do:
return [seq(acc[i], i=1..numelems(Y)/2)];
end proc:
Задаем параметры, данные нам в условии задачи
# Задаем параметры системы
m1 := 2.0;
L1 := 0.5;
f1 := 0.1;
g1 := 9.81;
Время задать основную логику модели, которая определяет расчет силы трения. При этом задаемся погрешностью скорости ползуна, при которой будем считать её равной нулю.
# Вычисление силы трения и нормальной реакции
GetFricNormal := proc(mass, length, grav_accel, fric_coeff, Y)
local F1, N1; # Сила трения и нормальная реакция
local eps_v := 1e-6; # Точность нуля скорости ползуна
# Если ползун неподвижен
if abs(Y[4]) < eps_v then
# Вычисляем силу трения покоя и нормальную реакцию
F1 := GetF1(mass, length, grav_accel, fric_coeff, Y);
N1 := GetN1(mass, length, grav_accel, fric_coeff, Y);
# Если трение покоя превышает максимально возможное значение,
# то сила трения равна силе трения скольжения
if abs(F1) > fric_coeff*abs(N1) then F1 := fric_coeff*abs(N1)*signum(F1); end if;
else
# Если ползун уже движется, считаем нормальную реакцию и трение скольжения
N1 := GetN2(mass, length, grav_accel, fric_coeff, Y);
F1 := fric_coeff*abs(N1)*signum(Y[4]);
end if;
return [F1, N1];
end proc:
Определяем callback для решателя
# Вычисление правой части ОДУ в форме Коши для численного решателя (собственно математическая модель)
EQs_func := proc(N, t, Y, dYdt)
local F1, N1; # Сила трения и нормальная реакция
local acc; # Обобщенные ускорения
local ret;
# Вычисляем силу трения и нормальную реакцию
ret := GetFricNormal(m1, L1, g1, f1, Y);
F1 := ret[1];
N1 := ret[2];
# Вычисляем производные фазовых координат
acc := GetAccel(m1, L1, g1, F1, N1, Y);
dYdt[1] := Y[4];
dYdt[2] := Y[5];
dYdt[3] := Y[6];
dYdt[4] := acc[1];
dYdt[5] := acc[2];
dYdt[6] := acc[3];
end proc:
Формируем для решателя список фазовых координат и начальные условия (угол отклонения стержня от вертикали делаем малым) и выполняем численное интегрирование (на самом деле последний вызов dsolve() лишь обозначает наши намерения по численному решению — оно будет поизведено при вычислении конкретных значений фазовых координат).
# Список фазовых координат
vars := [X(t), Y(t), Phi(t), Vx(t), Vy(t), Omega(t)];
# Начальные условия
initc := Array([0.0, 0.0, 1e-4, 0.0, 0.0, 0.0]);
# Интгрируем уравнения
dsolv := dsolve(numeric, number = 6, procedure = EQs_func, start = 0, initial = initc, procvars = vars, output=listprocedure);
Выполняем некоторые подготовительные операции
# Выделяем нужную нам часть решения
x := eval(X(t), dsolv);
y := eval(Y(t), dsolv);
phi := eval(Phi(t), dsolv);
vx := eval(Vx(t), dsolv);
vy := eval(Vy(t), dsolv);
omega := eval(Omega(t), dsolv);
Далее просчитаем движение системы в течение некоторого интервала времени и сформируем массивы данных для вывода на графики
# Рассчитываем движение системы на интересующем нас интервале времени
t0 := 0.0:
t1 := 10.0:
num_plots := 1000:
dt := (t1 - t0)/num_plots:
t := t0:
i := 1:
while t <= t1 do
Time[i] := t;
Y := [x(t), y(t), phi(t), vx(t), vy(t), omega(t)];
x1[i] := Y[1];
phi1[i] := Y[3];
fric[i] := GetFricNormal(m1, L1, g1, f1, Y)[1];
norm_react[i] := GetFricNormal(m1, L1, g1, f1, Y)[2];
lim_fric[i] := f1*abs(norm_react[i])*fric[i]/abs(fric[i]);
t := t + dt;
i := i + 1;
end do:
Ну вот, у нас практически всё готово
# Формируем графики
G_x := [ [Time[k],x1[k]] $k=1..num_plots]:
G_phi := [ [Time[k],phi1[k]] $k=1..num_plots]:
G_fric := [ [Time[k],fric[k]] $k=1..num_plots]:
G_norm_react := [ [Time[k],norm_react[k]] $k=1..num_plots]:
G_lim_fric := [ [Time[k],lim_fric[k]] $k=1..num_plots]:
# Выводим их на экран
gr_opts := captionfont=['ROMAN', 16], axesfont=['ROMAN','ROMAN', 12],titlefont=['ROMAN', 14],gridlines=true:
plot(G_x, gr_opts, view=[t0..t1, -1..1.0]);
plot(G_phi, gr_opts, view=[t0..t1, 0.0..7.0]);
plot({G_fric, G_lim_fric}, gr_opts, view=[t0..t1, -20..20]);
plot(G_norm_react, gr_opts, view=[t0..t1, 0.0..200.0]);
Получаем графики. Красоты ради, графики были конвертированы из Maple в *.eps и немножко обработаны в inkscape.
Перемещение ползуна
Угол отклонения стержня от вертикали
Сила трения
Здесь синей линией показано предельное значение трения покоя, а красной — фактическое значение силы трения.
Нормальная реакция со стороны направляющих
Видно, что ползун покоится в течение чуть более двух секунд, а затем, после преодоления трения покоя приходит в движение, которое постепенно затухает и полностью прекращается через 6,5 секунд после начала движения. После этого сила трения никогда не превышает предельного для покоя значения, ползун остается на месте, а стержень совершает гармонические колебания около устойчивого положения равновесия.
Заключение
Рассмотрено применение метода избыточных координат и уравнений Лагранжа 2 рода к анализу движения систем с сухим трением. Видно, что при внешней громоздкости получаемых результатов, процес синтеза уравнений движения может быть автоматизирован средствами символьной математики, а это существенно для современных технических задач.
Благодарю за внимание к моему труду.
Автор: maisvendoo
