Sprawozdanie - PSIWMUM, Lab6

Zadanie 1.1

Dany jest system mechaniczny o jednym stopniu swobody, dla którego równanie dynamiczne ruchu:

$$\begin{array}{lr}M\ddot{x}+C\dot{x}+Kx=F& (1)\end{array}$$

Mając zadane parametry systemu:

• Masa: $M=1$ [$\mathrm{kg}$],
• Tłumienie: $C=\pi$ [$\frac{\mathrm{N}}{\mathrm{m}}\cdot \mathrm{s}$],
• Sztywność sprężyny: $K=1089{\pi }^2$ [$\frac{\mathrm{N}}{\mathrm{m}}$].

Należało wyznaczyć:

• Częstotliwość drgań własnych: $f$ [$\mathrm{Hz}$],
• Bezwymiarowy współczynnik tłumienia: $\xi$,
• Stosunek amplitud wierzchołków kolejnych oscylacji: $\frac{x_0}{x_1}$.

W tym celu skorzystałem z poniższych zależności:

$$\begin{array}{lr}f=\frac{1}{2\pi }\cdot \sqrt{\frac{K}{M}}& (2)\\ \xi =\frac{C}{C_c}& (3)\end{array}$$

gdzie $C_c$ to tłumienie krytyczne, dane jako: $C_c=2M\cdot \sqrt{\frac{K}{M}}$ Dodatkowo, dla systemów oscylujących ($\xi \ll 1$):

$$\begin{array}{lr}\xi =\frac{\ln (\frac{x_0}{x_1})}{2\pi }& (4)\end{array}$$

Kod programu:

M = 1;                              % Masa
C = pi;                             % Tłumienie
K = 1089*(pi^2);                    % Sztywność sprężyny
Cc=2*M*sqrt(K/M);                   % Tłumienie krytyczne
f=sqrt(K/M)/(2*pi);                 % Częstotliwość własna
Xi=C/Cc;                            % Bezwymiarowy współczynnik tłumienia
X01=exp(2*pi*Xi);                   % Stosunek kolejnych amplitud

disp('f = '); disp(f)
disp('ξ = '); disp(Xi)
disp('x0/x1 = '); disp(X01)

Wynik działania programu:

f = 
   16.5000

ξ = 
    0.0152

x0/x1 = 
    1.0999

Zadanie 1.2

Następnie należało sprawdzić wyniki uzyskane w zadaniu 1.1 metodą numeryczną. W tym celu utworzyłem model w Simulinku układu mechanicznego opisanego równaniem (1), wymuszonego impulsem o amplitudzie 1, przy częstotliwości próbkowania $f_s=1000\mathrm{Hz}$.

Schemat modelu:

Kod programu:

Fs = 1000;                          % Częstotliwość próbkowania
M = 1;                              % Masa
C = pi;                             % Tłumienie
K = 1089*(pi^2);                    % Sztywność sprężyny

% Symulacja układu:
out = sim('Lab6_model_1_2023');
t=out.tout;                         % Wektor czasu
x=out.Displacement.Data;            % Przemieszczenie

% Szukanie pierwszych dwóch wierzchołków:
[x_peaks,location]=findpeaks([0,x'],'NPEAKS',2);

X01=x_peaks(1)/x_peaks(2);          % Stosunek kolejnych amplitud

Xi=log(X01)/(2*pi);                 % Bezwymiarowy współczynnik tłumienia

% Widmo amplitudowe:
N = length(x);                      % Ilość próbek
X = fft(x);                         % FFT z sygnału (x)
df = Fs/N;                          % Rozdzielczość
f_vec = (0:(N-1))*df;               % Wektor częstotliwości
X_amp = abs(X);                     % Amplituda

% Znormalizowana amplituda:
X_amp_norm = X_amp/N*2; X_amp_norm(1)=X_amp_norm(1)/2;

% Szukanie częstotliwości własnej:
[A_peaks,Location]=findpeaks([0,X_amp_norm'],'NPEAKS',1);

f=f_vec(Location);                  % Częstotliwość własna

figure(1)
subplot(2,1,1)
plot(t,x)
xlabel('Czas [s]');
ylabel('Przemieszczenie [m]');
title("Przebieg czasowy")
subplot(2,1,2)
plot(f_vec,X_amp_norm)
xlabel('Częstotliwość [Hz]');
ylabel('Amplituda');
title("Jednostronne widmo amplitudowe")
xlim([0 Fs/20])

disp('f = '); disp(f)
disp('ξ = '); disp(Xi)
disp('x0/x1 = '); disp(X01)

Wynik działania programu:

f = 
   16.5983

ξ = 
    0.0150

x0/x1 = 
    1.0989

Przy pomocy funkcji findpeaks jestem w stanie znaleźć wartości amplitud pierwszych dwóch wierzchołków. Następnie tworzę jednostronne widmo amplitudowe, z którego odczytuję częstotliwość własną.

Warto zauważyć, że identyczny wynik można było uzyskać korzystając z Toolbox’a SimScape, przy czym należy pamiętać o odwróceniu kierunku wymuszenia lub znaku zbieranych odpowiedzi.

Schemat modelu:

Porównanie wyników:

Simulink:
SimScape:

Zadanie 2

W tym zadaniu należało znaleźć taką wartość parametru tłumienia C, aby stosunek amplitud wierzchołków kolejnych oscylacji $\frac{x_0}{x_1}$ był równy 2, dla systemu mechanicznego z zadania 1.1.

Przekształcając wzór (4): $\ln (\frac{x_0}{x_1})=2\pi \cdot \xi$ Korzystam z własności logarytmu: ${\log }_{e}x=y\iff x=e^y$ A więc: $\frac{x_0}{x_1}=e^{2\pi \cdot \xi }$ Wstawiając wartość 2 za $\frac{x_0}{x_1}$ oraz wzór (3) za $\xi$: $2=e^{2\pi \cdot \frac{C}{C_c}}$ Ostatecznie: $e^{\frac{C}{33}}=2\Rightarrow C=33\ln (2)$ Kod programu:

Fs = 1000;                          % Częstotliwość próbkowania
M = 1;                              % Masa
C = 33*log(2);                      % Tłumienie
K = 1089*(pi^2);                    % Sztywność sprężyny

out = sim('Lab6_model_1_2023');
t=out.tout;                         % Wektor czasu
x=out.Displacement.Data;            % Przemieszczenie

[x_peaks,location]=findpeaks([0,x'],'NPEAKS',2);

% Położenie znaczników:
t1=t(location(1));
t2=t(location(2));
x1=x_peaks(1);
x2=x_peaks(2);

figure(1)
plot(t,x)
hold on
% Dodaj znaczniki:
plot(t1,x1,'_','MarkerSize',30,'MarkerEdgeColor','red','MarkerFaceColor',[1 .6 .6])
plot(t2,x2,'_','MarkerSize',30,'MarkerEdgeColor','red','MarkerFaceColor',[1 .6 .6])
xlabel('Czas [s]');
ylabel('Przemieszczenie [m]');
title('Przemieszczenie w czasie')
xlim([0,1])
hold off

Wynik działania programu:

Jak widać, dobrana wartość parametru C jest odpowiednia, gdyż stosunek $\frac{x_0}{x_1}$ jest równy w przybliżeniu 2.

Zadanie 3

W tym zadaniu należało wykorzystać model z zadania 1.2 do uzyskania funkcji odpowiedzi impulsowej (ang. Impulse Response Function - IRF), po czym powtórzyć symulację dla wymuszenia białym szumem. Na koniec trzeba było porównać wynik splotu IRF z białym szumem oraz odpowiedź systemu na wymuszenie białym szumem.

Schemat modelu:

Kod programu:

Fs = 1000;                          % Częstotliwość próbkowania
M = 1;                              % Masa
C = pi;                             % Tłumienie
K = 1089*(pi^2);                    % Sztywność sprężyny

out = sim('Lab6_model_1_2023');
t=out.tout;                         % Wektor czasu
IRF=out.Displacement.Data;          % Funkcja odpowiedzi impulsowej [1]

out = sim('Lab6_model_3_2023');
x2=out.Wymuszenie.Data;             % Wektor wymuszenia [2]
x3=out.Odp_Systemu.Data;            % Odpowiedź systemu [3]

g=conv(IRF,x2);                     % Splot IRF z wymuszeniem

figure(1)
subplot(2,1,1)
plot(t,g(1:length(x3)))
xlabel('Czas [s]');
ylabel('Amplituda [-]');
title('Splot IRF z wymuszeniem drugiej symulacji')
subplot(2,1,2)
plot(t,x3)
xlabel('Czas [s]');
ylabel('Amplituda [-]');
title('Odpowiedź systemu w drugiej symulacji')

Wynik działania programu:

Jak widać, uzyskane wykresy są identyczne, co się zgadza z teorią, gdyż odpowiedź systemu $y(t)$ na wymuszenie $x(t)$ można uzyskać na dwa sposoby:

• Poprzez symulację układu:
• Poprzez splot wymuszenia z IRF: $y(t)=x(t)\ast h(t)$

Zadanie 4

W tym zadaniu należało utworzyć model systemu mechanicznego o dwóch stopniach swobody, opisanego równaniami:

$$\{\begin{array}{l}{M}_1{\ddot{x}}_1+(C_1+C_2){\dot{x}}_1-C_2{\dot{x}}_2+(K_1+K_2)x_1-K_2{x}_2=F_1\\ M_2{\ddot{x}}_2+(C_2+C_3){\dot{x}}_2-C_2{\dot{x}}_1+(K_2+K_3)x_2-K_2{x}_1=F_2\end{array}$$

Zadane parametry:

• $M_1=1$ [$\mathrm{kg}$],
• $M_2=1$ [$\mathrm{kg}$],
• $K_1=10^4$ [$\frac{\mathrm{N}}{\mathrm{m}}$],
• $K_2=10^5$ [$\frac{\mathrm{N}}{\mathrm{m}}$],
• $K_3=10^4$ [$\frac{\mathrm{N}}{\mathrm{m}}$],
• $C_1=10$ [$\frac{\mathrm{N}}{\mathrm{m}}\cdot \mathrm{s}$],
• $C_2=10$ [$\frac{\mathrm{N}}{\mathrm{m}}\cdot \mathrm{s}$],
• $C_3=10$ [$\frac{\mathrm{N}}{\mathrm{m}}\cdot \mathrm{s}$].

Schemat modelu:

Zadanie 5

W tym zadaniu należało utworzyć model systemu mechanicznego o 10 stopniach swobody. Parametry identyczne jak w zadaniu 1.2.

Schemat modelu:

Zadanie 6

W tym zadaniu należało wykonać symulacje modelu z zadania 5, po czym wykreślić przebiegi czasowe i ich widma amplitudowo-częstotliwościowe dla wybranej masy. Na koniec należało odczytać częstotliwości rezonansowe badanego układu.

Dane:

• Czas pierwszej symulacji: $t_1=100\mathrm{s}$,
• Czas drugiej symulacji: $t_2=1\mathrm{s}$,
• Częstotliwość próbkowania: $f_s=1000\mathrm{Hz}$.

Kod programu:

A = 10000;                          % Amplituda wymuszenia
Fs = 1000;                          % Częstotliwość próbkowania
M = 1;                              % Masa
C = pi;                             % Tłumienie
K = 1089*(pi^2);                    % Sztywność sprężyny
simTime1=100;                       % Czas symulacji
simTime2=1;

out = sim('Lab6_model_10DOF_2023', "StopTime",'simTime1');
t1=out.tout;                        % Wektor czasu
x1=out.M10;                         % Masa 10

N1 = length(x1);                    % Ilość próbek
X1 = fft(x1);                       % FFT z sygnału
df1 = Fs/N1;                        % Rozdzielczość
f_vec1 = (0:(N1-1))*df1;            % Wektor częstotliwości
X1_amp = abs(X1);                   % Amplituda
X1_amp_norm = X1_amp/N1*2; X1_amp_norm(1)=X1_amp_norm(1)/2;

out = sim('Lab6_model_10DOF_2023', "StopTime",'simTime2');
t2=out.tout;                        % Wektor czasu
x2=out.M10;                         % Masa 10

N2 = length(x2);                    % Ilość próbek
X2 = fft(x2);                       % FFT z sygnału
df2 = Fs/N2;                        % Rozdzielczość
f_vec2 = (0:(N2-1))*df2;            % Wektor częstotliwości
X2_amp = abs(X2);                   % Amplituda
X2_amp_norm = X2_amp/N2*2; X2_amp_norm(1)=X2_amp_norm(1)/2;

figure(1)
subplot(2,1,1)
plot(t1,x1)
xlabel('Czas [s]');
ylabel('Przemieszczenie masy 1 [m]');
title('t=100 s')
subplot(2,1,2)
semilogy(f_vec1,X1_amp_norm)
xlim([0 50])
xlabel('Częstotliwość [Hz]')
ylabel('Amplituda')
title('Jednostronne widmo amplitudowe w skali logarytmicznej dla t=100 s')

figure(2)
subplot(2,1,1)
plot(t2,x2)
xlabel('Czas [s]');
ylabel('Przemieszczenie masy 1 [m]');
title('t=1 s')
subplot(2,1,2)
semilogy(f_vec2,X2_amp_norm)
xlim([0 50])
xlabel('Częstotliwość [Hz]')
ylabel('Amplituda')
title('Jednostronne widmo amplitudowe w skali logarytmicznej dla t=1 s')

[X1_peakValue,X1_location]=findpeaks([0,X1_amp_norm'],'NPEAKS',8);
ResonanceFreq_1=f_vec1(X1_location(1:end));

disp('Częstotliwości rezonansowe dla t=100 s:'); disp(ResonanceFreq_1)

[X2_peakValue,X2_location]=findpeaks([0,X2_amp_norm'],'NPEAKS',8);
ResonanceFreq_2=f_vec2(X2_location(1:end));

disp('Częstotliwości rezonansowe dla t=1 s:'); disp(ResonanceFreq_2)

Wynik działania programu:

Częstotliwości rezonansowe dla t=100 s:
    2.4800    7.3499   12.0699   16.5298   20.6498   24.3598   27.5397   30.0497

Częstotliwości rezonansowe dla t=1 s:
    2.9970    7.9920   12.9870   17.9820   21.9780   25.9740   28.9710   30.9690

Jak widać, dla czasu symulacji $t=1\mathrm{s}$, rozdzielczość jest zbyt mała, aby częstotliwości rezonansowe zostały odczytane poprawnie.