Формализм Лагранжа в задачах с сухим трением

в 22:36, , рубрики: maple, аналитическая механика, математика, математическое моделирование, Программирование, сухое трение

Введение

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

Формализм Лагранжа в задачах с сухим трением - 1

Тонкий однородный стержень массы m = 2 кг, длины AB = 2l = 1 м в точке A шарнирно прикреплен к невесомому ползуну, перемещающемуся в горизонтальных шероховатых направляющих. В начальный момент времени стержень расположен вертикально, затем его отклоняют от вертикали на ничтожно малый угол и отпускают без начальной скорости. Необходимо составить уравнения движения данной механической системы и найти закон её движения. Коэффициент трения между ползуном и направляющими равен f = 0,1.

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

1. Что может быть «проще» трения?

Нет более страшного наказания для механика, чем сила трения. Появляясь в задаче, эта сила сразу делает её существенно нелинейной, ибо ведет себя достаточно интересным образом.

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

Формализм Лагранжа в задачах с сухим трением - 2

Пусть в начале к нему не прилагают никаких сил (кроме силы тяжести и нормальной реакции). В этом случае и сила трения между бруском и плоскостью будет равна нулю.

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

Формализм Лагранжа в задачах с сухим трением - 3

Будем постепенно увеличивать силу Формализм Лагранжа в задачах с сухим трением - 4 и, согласно (1) расти будет и сила трения, которая в этом случае называется силой трения покоя. Так будет продолжаться до тех пор, пока сила трения покоя не достигнет величины

Формализм Лагранжа в задачах с сухим трением - 5

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

Формализм Лагранжа в задачах с сухим трением - 6

где Формализм Лагранжа в задачах с сухим трением - 7 — скорость бруска.

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

Если точка, где приложена сила трения неподвижна:

  1. Расчитываем силу трения покоя и нормальную реакцию
  2. Проверяем условие
    Формализм Лагранжа в задачах с сухим трением - 8
    при нарушении которого принимаем силу трения равной предельной силе трения покоя

Если точка приложения силы трения движется:

  1. Вычисляем нормальную реакцию
  2. Вычисляем силу трения скольжения, согласно выражению (2)

2. Моделирование движения системы с трением

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

Формализм Лагранжа в задачах с сухим трением - 9

Здесь в качестве обобщенных координат берем вектор

Формализм Лагранжа в задачах с сухим трением - 10

где x,y — координаты точки A; Формализм Лагранжа в задачах с сухим трением - 11 — угол наклона стержня к вертикали. Вооружаемся 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 выдает такой результат

Формализм Лагранжа в задачах с сухим трением - 12

Довольно громоздко, но «ковырять» нам не руками на листочке, поэтому двигаемся дальше. Задаем векторы и точки приложения сил

# Векторы сил, приоложенных к системе 
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):

Получаем трёх «крокодилов»

Формализм Лагранжа в задачах с сухим трением - 13

Формализм Лагранжа в задачах с сухим трением - 14

Формализм Лагранжа в задачах с сухим трением - 15

Эти уравнения пришлось вбить в статью руками, ибо «копипаста» LaTeX-вывода Maple приводит к неприглядному виду результата. Но даже так видно — уравнения сложны и с учетом того что F — это сила трения, аналитически не интегрируемы.

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

Формализм Лагранжа в задачах с сухим трением - 16

Кроме того, в том случае когда ползун неподвижен, а сила трения покоя не превысила предельного значения, включается ещё одна связь

Формализм Лагранжа в задачах с сухим трением - 17

где Формализм Лагранжа в задачах с сухим трением - 18 — некоторая горизонтальная координата ползуна. Теперь преобразуем систему (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

Формализм Лагранжа в задачах с сухим трением - 19

Формализм Лагранжа в задачах с сухим трением - 20

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

Формализм Лагранжа в задачах с сухим трением - 21

(минус уже имеется в уравнениях движения)

# Трансформируем систему для поиска нормальной реакции при скольжении ползуна
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-код

Формализм Лагранжа в задачах с сухим трением - 22

Все необходимые нам выражения получены, теперь можно переходить к моделированию. В отличие от задачи с маятником, о которой я уже писал, тут мы честно трансформируем наши уравнения 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.

Перемещение ползуна
Формализм Лагранжа в задачах с сухим трением - 23
Угол отклонения стержня от вертикали
Формализм Лагранжа в задачах с сухим трением - 24
Сила трения
Формализм Лагранжа в задачах с сухим трением - 25
Здесь синей линией показано предельное значение трения покоя, а красной — фактическое значение силы трения.

Нормальная реакция со стороны направляющих
Формализм Лагранжа в задачах с сухим трением - 26

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

Заключение

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

Благодарю за внимание к моему труду.

Автор: maisvendoo

Источник

* - обязательные к заполнению поля


https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js