Введение

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

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

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

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

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

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

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

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

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

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

где — скорость бруска.

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

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

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

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

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

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