Этот пример исследует физику демпфированного гармонического генератора
решение уравнений движения в случае отсутствия движущих сил,
расследование случаев недостаточного, чрезмерного и критического демпфирования

Вывод уравнения движения
Решить уравнение движения (F = 0)
Недостаточно демпфированный случай ()
Перегруженный корпус ()
Критически демпфированный случай ()
Заключение
Рассмотрим принудительный гармонический генератор с демпфированием, показанным ниже. Моделируйте силу сопротивления как пропорциональную скорости, с которой движется генератор.

Задайте уравнение движения где
- масса
- коэффициент демпфирования
является коэффициент упругости
является движущей силой
syms x(t) m c k F(t) eq = m*diff(x,t,t) + c*diff(x,t) + k*x == F
eq(t) =
Переписать уравнение используя и .
syms gamma omega_0 eq = subs(eq, [c k], [m*gamma, m*omega_0^2])
eq(t) =
Разделите массу . Теперь у нас есть уравнение в удобной форме для анализа.
eq = collect(eq, m)/m
eq(t) =
Решить уравнение движения можно используя dsolve в случае отсутствия внешних сил, где . Используйте начальные условия модуля и нулевой скорости.
vel = diff(x,t); cond = [x(0) == 1, vel(0) == 0]; eq = subs(eq,F,0); sol = dsolve(eq, cond)
sol =
Рассмотрим, как упростить решение путем его расширения.
sol = expand(sol)
sol =
Заметьте, что каждый термин имеет коэффициент , или , использовать collect чтобы собрать эти условия
sol = collect(sol, exp(-gamma*t/2))
sol =
Термин появляется в различных частях решения. Перепишите его в более простой форме путем введения коэффициента затухания .
Подстановка
syms zeta; sol = subs(sol, ... sqrt(gamma^2 - 4*omega_0^2), ... 2*omega_0*sqrt(zeta^2-1))
sol =
Дополнительно упростите решение путем замены с точки зрения и ,
sol = subs(sol, gamma, 2*zeta*omega_0)
sol =
Мы вывели общее решение для движения демпфированного гармонического генератора без движущих сил. Далее мы исследуем три особых случая коэффициента затухания где движение принимает более простые формы. Эти случаи называются
underdamped ,
сверхзаглушенный , и
критически демпфирован .
Если , затем чисто мнимый
solUnder = subs(sol, sqrt(zeta^2-1), 1i*sqrt(1-zeta^2))
solUnder =
Заметьте условия в вышеприведенном уравнении и напомнить о тождествах
Перепишите решение с точки зрения .
solUnder = coeffs(solUnder, zeta);
solUnder = solUnder(1);
c = exp(-omega_0 * zeta * t);
solUnder = c * rewrite(solUnder / c, 'cos')solUnder =
solUnder(t, omega_0, zeta) = solUnder
solUnder(t, omega_0, zeta) =
Система колеблется на естественной частоте и распадается с экспоненциальной скоростью .
Постройте график решения с помощью fplot как функцию и .
z = [0 1/4 1/2 3/4];
w = 1;
T = 4*pi;
lineStyle = {'-','--',':k','-.'};
fplot(@(t)solUnder(t, w, z(1)), [0 T], lineStyle{1});
hold on;
for k = 2:numel(z)
fplot(@(t)solUnder(t, w, z(k)), [0 T], lineStyle{k});
end
hold off;
grid on;
xticks(T*linspace(0,1,5));
xticklabels({'0','\pi','2\pi','3\pi','4\pi'});
xlabel('t / \omega_0');
ylabel('amplitude');
lgd = legend('0','1/4','1/2','3/4');
title(lgd,'\zeta');
title('Underdamped');
Если , затем является чисто реальным, и решение может быть переписано как
solOver = sol
solOver =
solOver = coeffs(solOver, zeta); solOver = solOver(1)
solOver =
Заметьте условия и вспомнить тождества .
Переписать выражение в терминах .
c = exp(-omega_0*t*zeta);
solOver = c*rewrite(solOver / c, 'cosh')solOver =
solOver(t, omega_0, zeta) = solOver
solOver(t, omega_0, zeta) =
Постройте график решения, чтобы увидеть, что оно распадается, не колебаясь.
z = 1 + [1/4 1/2 3/4 1];
w = 1;
T = 4*pi;
lineStyle = {'-','--',':k','-.'};
fplot(@(t)solOver(t, w, z(1)), [0 T], lineStyle{1});
hold on;
for k = 2:numel(z)
fplot(@(t)solOver(t, w, z(k)), [0 T], lineStyle{k});
end
hold off;
grid on;
xticks(T*linspace(0,1,5));
xticklabels({'0','\pi','2\pi','3\pi','4\pi'});
xlabel('\omega_0 t');
ylabel('amplitude');
lgd = legend('1+1/4','1+1/2','1+3/4','2');
title(lgd,'\zeta');
title('Overdamped');
Если , затем решение упрощает
solCritical(t, omega_0) = limit(sol, zeta, 1)
solCritical(t, omega_0) =
Постройте график решения для критически демпфированного случая.
w = 1; T = 4*pi; fplot(solCritical(t, w), [0 T]) xlabel('\omega_0 t'); ylabel('x'); title('Critically damped, \zeta = 1'); grid on; xticks(T*linspace(0,1,5)); xticklabels({'0','\pi','2\pi','3\pi','4\pi'});

Мы рассмотрели различные состояния демпфирования для гармонического генератора, решив ОДУ, которая представляет свое движение, используя коэффициент затухания . Постройте все три случая вместе, чтобы сравнить и контрастировать их.
zOver = pi;
zUnder = 1/zOver;
w = 1;
T = 2*pi;
lineStyle = {'-','--',':k'};
fplot(@(t)solOver(t, w, zOver), [0 T], lineStyle{1},'LineWidth',2);
hold on;
fplot(solCritical(t, w), [0 T], lineStyle{2},'LineWidth',2)
fplot(@(t)solUnder(t, w, zUnder), [0 T], lineStyle{3},'LineWidth',2);
hold off;
textColor = lines(3);
text(3*pi/2, 0.3 , 'over-damped' ,'Color',textColor(1,:));
text(pi*3/4, 0.05, 'critically-damped','Color',textColor(2,:));
text(pi/8 , -0.1, 'under-damped');
grid on;
xlabel('\omega_0 t');
ylabel('amplitude');
xticks(T*linspace(0,1,5));
xticklabels({'0','\pi/2','\pi','3\pi/2','2\pi'});
yticks((1/exp(1))*[-1 0 1 2 exp(1)]);
yticklabels({'-1/e','0','1/e','2/e','1'});
lgd = legend('\pi','1','1/\pi');
title(lgd,'\zeta');
title('Damped Harmonic Oscillator');