Физика ослабленного гармонического генератора

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

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

  • исследование случаев под - сверх - и критическое затухание

Содержимое

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

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

  3. Случай Underdamped (ζ<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) = 

m ∂2∂t2 x(t)+c ∂∂t x(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) = 

m ∂2∂t2 x(t)+γ m ∂∂t x(t)+m x(t) ω02=F(t)

Отделите массу m. Теперь у нас есть уравнение в удобной форме, чтобы анализировать.

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

∂2∂t2 x(t)+γ ∂∂t x(t)+x(t) ω02=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

Исследуйте, как упростить решение путем расширения его.

sol = expand(sol)
sol = 

σ1 σ22+σ1 et γ2-4 ω0222-γ σ1 σ22 γ2-4 ω02+γ σ1 et γ2-4 ω0222 γ2-4 ω02where  σ1=e-γ t2  σ2=e-t γ2-4 ω022

Заметьте, что каждый термин имеет фактор σ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 ω02 e-γ t2where  σ1=t γ2-4 ω022

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

Замена ζ в термин выше дает:

γ2−4ω02=2ω0(γ2ω0)2−1=2ω0ζ2−1,

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-ω0 t ζ2-1  σ2=eω0 t ζ2-1

Далее упростите решение путем замены γ в терминах ω0 и ζ,

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

e-ω0 t ζ σ22+σ12+ζ σ22 ζ2-1-ζ σ12 ζ2-1where  σ1=e-ω0 t ζ2-1  σ2=eω0 t ζ2-1

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

  • underdamped (ζ<1),

  • сверхослабленный (ζ>1), и

  • критически ослабленный (ζ=1).

3. Случай Underdamped (ζ<1)

Если ζ<1то ζ2-1=i1-ζ2 является чисто мнимым

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

e-ω0 t ζ σ12+σ22+ζ σ1 i2 1-ζ2-ζ σ2 i2 1-ζ2where  σ1=e-ω0 t 1-ζ2 i  σ2=eω0 t 1-ζ2 i

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

Перепишите решение в терминах cos.

solUnder = coeffs(solUnder, zeta);
solUnder = solUnder(1);
c = exp(-omega_0 * zeta * t);
solUnder = c * rewrite(solUnder / c, 'cos')
solUnder = e-ω0 t ζ cos(ω0 t 1-ζ2)
solUnder(t, omega_0, zeta) = solUnder
solUnder(t, omega_0, zeta) = e-ω0 t ζ cos(ω0 t 1-ζ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 object. The axes object 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-ω0 t ζ σ22+σ12+ζ σ22 ζ2-1-ζ σ12 ζ2-1where  σ1=e-ω0 t ζ2-1  σ2=eω0 t ζ2-1

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

e-ω0 t ζ eω0 t ζ2-12+e-ω0 t ζ2-12

Заметьте термины (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(ω0 t ζ2-1) e-ω0 t ζ
solOver(t, omega_0, zeta) = solOver
solOver(t, omega_0, zeta) = cosh(ω0 t ζ2-1) e-ω0 t ζ

Постройте решение видеть, что оно затухает без колебания.

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 object. The axes object 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-ω0 t ω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 object. The axes object with title C r i t i c a l l y blank d a m p e d , blank zeta blank = blank 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 object. The axes object with title Damped Harmonic Oscillator contains 6 objects of type functionline, text. These objects represent \pi, 1, 1/\pi.