Физика демпфированного гармонического генератора

Этот пример исследует физику демпфированного гармонического генератора

  • решение уравнений движения в случае отсутствия движущих сил,

  • расследование случаев недостаточного, чрезмерного и критического демпфирования

Содержание

  1. Вывод уравнения движения

  2. Решить уравнение движения (F = 0)

  3. Недостаточно демпфированный случай (ζ<1)

  4. Перегруженный корпус (ζ>1)

  5. Критически демпфированный случай (ζ=1)

  6. Заключение

1. Вывод уравнения движения

Рассмотрим принудительный гармонический генератор с демпфированием, показанным ниже. Моделируйте силу сопротивления как пропорциональную скорости, с которой движется генератор.

Задайте уравнение движения где

  • m - масса

  • c - коэффициент демпфирования

  • k является коэффициент упругости

  • F является движущей силой

syms x(t) m c k F(t)
eq = m*diff(x,t,t) + c*diff(x,t) + k*x == F
eq(t) = 

m2t2 x(t)+ct x(t)+kx(t)=F(t)m * diff (x (t), t, 2) + c * diff (x (t), t) + k * x (t) = = F (t)

Переписать уравнение используя c=mγ и k=mω02.

syms gamma omega_0
eq = subs(eq, [c k], [m*gamma, m*omega_0^2])
eq(t) = 

m2t2 x(t)+γmt x(t)+mx(t)ω02=F(t)m * diff (x (t), t, 2) + gamma * m * diff (x (t), t) + m * x (t) * omega_0^2 = = F (t)

Разделите массу m. Теперь у нас есть уравнение в удобной форме для анализа.

eq = collect(eq, m)/m
eq(t) = 

2t2 x(t)+γt x(t)+x(t)ω02=F(t)mdiff (x (t), t, 2) + gamma * diff (x (t), t) + x (t) * omega_0^2 = = F (t )/m

2. Решить уравнение движения, где F = 0

Решить уравнение движения можно используя dsolve в случае отсутствия внешних сил, где F=0. Используйте начальные условия модуля и нулевой скорости.

vel = diff(x,t);
cond = [x(0) == 1, vel(0) == 0];
eq = subs(eq,F,0);
sol = dsolve(eq, cond)
sol = 

e-tγ2-σ12γ+σ12σ1-e-tγ2+σ12γ-σ12σ1where  σ1=γ-2ω0γ+2ω0(exp ((-t * (gamma/2 - sqrt ((gamma - 2 * omega _ 0) * (gamma + 2 * omega _ 0) )/2)) * (gamma + sqrt ((gamma - 2 * omega _ 0) * (gamma + 2 * omega _ 0))) )/( 2 * sega _)

Рассмотрим, как упростить решение путем его расширения.

sol = expand(sol)
sol = 

σ1σ22+σ1etγ2-4ω0222-γσ1σ22γ2-4ω02+γσ1etγ2-4ω0222γ2-4ω02where  σ1=e-γt2  σ2=e-tγ2-4ω022 (exp (-(gamma*t)/2) *exp (-(t*sqrt (gamma^sym (2) - 4*omega_0^2))/2))/2 + (exp (-(gamma*t)/2) *exp (((t*sqrt (gamma^sym (2) - 4*omega_0^2))/2)))/2 - (gamma*exp (-(gamma*t)/2) *exp (-(t*sqrt (gamma^sym (2) - 4*omega_0^2))/2)) / (2*sqrt (gamma^sym (2) - 4*omega_0^2)) + (gamma*exp (-(gamma*t)/2) *exp (((t*sqrt (gamma^sym (2) - 4*omega_0^2))/2))) / (2*sqrt (gamma^sym (2) - 4*omega_0^2))

Заметьте, что каждый термин имеет коэффициент σ1, или e-γt/2, использовать collect чтобы собрать эти условия

sol = collect(sol, exp(-gamma*t/2))
sol = 

e-σ12+eσ12-γe-σ12γ2-4ω02+γeσ12γ2-4ω02e-γt2where  σ1=tγ2-4ω022(exp ((- (t * sqrt (gamma ^ sym (2) - 4 * omega _ 0 ^ 2) )/2 )/2 + exp ((t * sqrt (gamma ^ sym (2) - 4 * omega _ 0 ^ 2 )/2 )/2 - (gamma * exp ((- (t * sqrt (g

Термин γ2-4ω02 появляется в различных частях решения. Перепишите его в более простой форме путем введения коэффициента затухания ζγ2ω0.

Подстановка

γ24ω02=2ω0(γ2ω0)21=2ω0ζ21,

syms zeta;
sol = subs(sol, ...
    sqrt(gamma^2 - 4*omega_0^2), ...
    2*omega_0*sqrt(zeta^2-1))
sol = 

e-γt2σ22+σ12+γσ24ω0ζ2-1-γσ14ω0ζ2-1where  σ1=e-ω0tζ2-1  σ2=eω0tζ2-1exp ((-(gamma*t)/2)) * (exp (omega_0*t*sqrt (zeta^sym (2) - 1))/2 + exp (-omega_0*t*sqrt (zeta^sym (2) - 1))/2 + (gamma*exp (omega_0*t*sqrt (zeta^sym (2) - 1))) / (4*omega_0*sqrt (zeta^sym (2) - 1)) - (gamma*exp (-omega_0*t*sqrt (zeta^sym (2) - 1))) / (4*omega_0*sqrt (zeta^sym (2) - 1)))

Дополнительно упростите решение путем замены γ с точки зрения ω0 и ζ,

sol = subs(sol, gamma, 2*zeta*omega_0)
sol = 

e-ω0tζσ22+σ12+ζσ22ζ2-1-ζσ12ζ2-1where  σ1=e-ω0tζ2-1  σ2=eω0tζ2-1exp ((-omega_0*t*zeta)) * (exp (omega_0*t*sqrt (zeta^sym (2) - 1))/2 + exp (-omega_0*t*sqrt (zeta^sym (2) - 1))/2 + (zeta*exp (omega_0*t*sqrt (zeta^sym (2) - 1))) / (2*sqrt (zeta^sym (2) - 1)) - (zeta*exp (-omega_0*t*sqrt (zeta^sym (2) - 1))) / (2*sqrt (zeta^sym (2) - 1)))

Мы вывели общее решение для движения демпфированного гармонического генератора без движущих сил. Далее мы исследуем три особых случая коэффициента затухания ζ где движение принимает более простые формы. Эти случаи называются

  • underdamped (ζ<1),

  • сверхзаглушенный (ζ>1), и

  • критически демпфирован (ζ=1).

3. Недостаточно демпфированный случай (ζ<1)

Если ζ<1, затем ζ2-1=i1-ζ2 чисто мнимый

solUnder = subs(sol, sqrt(zeta^2-1), 1i*sqrt(1-zeta^2))
solUnder = 

e-ω0tζσ12+σ22+ζσ1i21-ζ2-ζσ2i21-ζ2where  σ1=e-ω0t1-ζ2i  σ2=eω0t1-ζ2iexp ((-omega_0*t*zeta)) * (exp (-omega_0*t*sqrt (1 - zeta^sym (2)) *sym (1i))/2 + exp (omega_0*t*sqrt (1 - zeta^sym (2)) *sym (1i))/2 + (zeta*exp (-omega_0*t*sqrt (1 - zeta^sym (2)) *sym (1i)) *sym (1i)) / (2*sqrt (1 - zeta^sym (2))) - (zeta*exp (omega_0*t*sqrt (1 - zeta^sym (2)) *sym (1i)) *sym (1i)) / (2*sqrt (1 - zeta^sym (2))))

Заметьте условия eiω0tζ2-1±e-iω0tζ2-1 в вышеприведенном уравнении и напомнить о тождествах eix=cos(x)+isin(x).

Перепишите решение с точки зрения cos.

solUnder = coeffs(solUnder, zeta);
solUnder = solUnder(1);
c = exp(-omega_0 * zeta * t);
solUnder = c * rewrite(solUnder / c, 'cos')
solUnder = e-ω0tζcos(ω0t1-ζ2)exp ((-omega _ 0 * t * zeta)) * cos (omega_0*t*sqrt (1 - zeta ^ sym (2))
solUnder(t, omega_0, zeta) = solUnder
solUnder(t, omega_0, zeta) = e-ω0tζcos(ω0t1-ζ2)exp ((-omega _ 0 * t * zeta)) * cos (omega_0*t*sqrt (1 - zeta ^ sym (2))

Система колеблется на естественной частоте ω01-ζ2 и распадается с экспоненциальной скоростью 1/ω0ζ.

Постройте график решения с помощью fplot как функцию ω0t и ζ.

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');

Figure contains an axes. The axes with title Underdamped contains 4 objects of type functionline. These objects represent 0, 1/4, 1/2, 3/4.

4. Перегруженный корпус (ζ>1)

Если ζ>1, затем ζ2-1 является чисто реальным, и решение может быть переписано как

solOver = sol
solOver = 

e-ω0tζσ22+σ12+ζσ22ζ2-1-ζσ12ζ2-1where  σ1=e-ω0tζ2-1  σ2=eω0tζ2-1exp ((-omega_0*t*zeta)) * (exp (omega_0*t*sqrt (zeta^sym (2) - 1))/2 + exp (-omega_0*t*sqrt (zeta^sym (2) - 1))/2 + (zeta*exp (omega_0*t*sqrt (zeta^sym (2) - 1))) / (2*sqrt (zeta^sym (2) - 1)) - (zeta*exp (-omega_0*t*sqrt (zeta^sym (2) - 1))) / (2*sqrt (zeta^sym (2) - 1)))

solOver = coeffs(solOver, zeta);
solOver = solOver(1)
solOver = 

e-ω0tζeω0tζ2-12+e-ω0tζ2-12exp ((-omega _ 0 * t * zeta)) * (exp (omega_0*t*sqrt (zeta ^ sym (2) - 1) )/2 + exp ((-omega _ 0 * t * sqrt (zeta ^ sym (2) - 1)) )/2)

Заметьте условия (eω0tζ2-1+e-ω0tζ2-1)2 и вспомнить тождества cosh(x)=ex+e-x2.

Переписать выражение в терминах cosh.

c = exp(-omega_0*t*zeta);
solOver = c*rewrite(solOver / c, 'cosh')
solOver = cosh(ω0tζ2-1)e-ω0tζcosh (omega_0*t*sqrt (zeta ^ sym (2) - 1)) * exp ((-omega _ 0 * t * zeta))
solOver(t, omega_0, zeta) = solOver
solOver(t, omega_0, zeta) = cosh(ω0tζ2-1)e-ω0tζcosh (omega_0*t*sqrt (zeta ^ sym (2) - 1)) * exp ((-omega _ 0 * t * 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');

Figure contains an axes. The axes with title Overdamped contains 4 objects of type functionline. These objects represent 1+1/4, 1+1/2, 1+3/4, 2.

5. Критически демпфированный случай (ζ=1)

Если ζ=1, затем решение упрощает

solCritical(t, omega_0) = limit(sol, zeta, 1)
solCritical(t, omega_0) = e-ω0tω0t+1exp ((-omega _ 0 * t)) * (omega_0*t + 1)

Постройте график решения для критически демпфированного случая.

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'});

Figure contains an axes. The axes with title Critically damped, \zeta = 1 contains an object of type functionline.

6. Заключение

Мы рассмотрели различные состояния демпфирования для гармонического генератора, решив ОДУ, которая представляет свое движение, используя коэффициент затухания ζ. Постройте все три случая вместе, чтобы сравнить и контрастировать их.

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');

Figure contains an axes. The axes with title Damped Harmonic Oscillator contains 6 objects of type functionline, text. These objects represent \pi, 1, 1/\pi.