-
Notifications
You must be signed in to change notification settings - Fork 0
/
Homework_8.m
169 lines (156 loc) · 6.74 KB
/
Homework_8.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
%% MECHANICAL VIBRATIONS (2021/1) - HOMEWORK 8
% Docente: Michael John Brennan
% Discente: Estevao Fuzaro de Almeida
% Data: 06/05/2021
% INICIALIZACAO
clc; clear all; close all; format long; %#ok<*CLALL>
set(groot,'defaultAxesTickLabelInterpreter','latex');
set(groot,'defaultLegendInterpreter','latex');
set(groot,'defaultTextInterpreter','latex');
txtsize = 26;
lgndsize = 18;
%% VARIAVEIS
ms = 1; % Massa principal [kg]
ks = 1e4; % Rigidez principal [N/m]
zs = 0.001; % Zeta principal [adimensional]
mu = 0.1; % mu = ma/ms [adimensional]
Fs = 1000; % Freq. de Amostragem [Hz]
T = 60; % Periodo [s]
dt = 1/Fs; % Incremento de Tempo [s]
df = 1/T; % Incremento de Frequencia [Hz]
t = 0:dt:T; % Vetor de Tempo [s]
f = 0:df:Fs; % Vetor de Frequencia [Hz]
w = 2*pi*f; % Velocidade Angular [rad/s]
%% PARAMETROS DO SISTEMA
wn = sqrt(ks/ms); % Freq. Natural [rad/s]
cs = 2*zs*sqrt(ks*ms); % Amortecimento principal [N.s/m]
wa = wn/(1+mu); % Freq. Nat. do absorvedor [Hz]
ma = ms*mu; % Massa do absorvedor [kg]
ka = ma*wa^2; % Rigidez do absorvedor [N/m]
za = sqrt((3/8)*(mu/(1+mu)^3)); % Zeta do absorvedor [adimensional]
ca = 2*za*sqrt(ka*ma); % Amortecimento do absorvedor [N.s/m]
%% CALCULO DAS FRF's Hs E Ha
M = [ms 0; 0 ma];
K = [ks+ka -ka; -ka ka]; % Modelo
C = [cs+ca -ca; -ca ca]; % matricial
F = [1; 0];
st=0;
for fAux=0:df:Fs
st=st+1;
wAux = 2*pi*fAux;
D = K - wAux.^2*M + 1i*wAux*C;
H = D\F;
Hs(st) = H(1); %#ok<*SAGROW>
Ha(st) = H(2);
end
% PLOTANDO Hs E Ha
figure
set(gcf,'Units','Normalized','OuterPosition',[0 0 1 0.6])
loglog(f,abs(Hs),'r','linewidth', 2), hold on
loglog(f,abs(Ha),'k--','linewidth', 2), hold on
xlabel('$f$ [Hz]')
ylabel('Receptance [m/N]')
legend({'$H_s$','$H_a$'},'Location','northeast','fontsize',lgndsize)
grid on, grid minor
axis([1e-1 Fs/2 1e-7 2e-3])
set(gca,'fontsize',txtsize,'Xtick',[1e-2 1e-1 1e0 1e1 1e2 1e3],'Ytick',[1-7 1e-6 1e-5 1e-4 1e-3],'XColor','k','YColor','k','ZColor','k','GridColor','k')
axes('Position',[.25 .25 .3 .4]); box on
loglog(f,abs(Hs),'r','linewidth',2), hold on, grid on
loglog(f,abs(Ha),'k--','linewidth', 2), hold on
xlim([5 25]); ylim([1e-4 2e-3]);
grid on, grid minor
set(gca,'fontsize',txtsize-8,'Xtick',[1e-2 1e-1 1e0 1e1 1e2 1e3],'Ytick',[1-7 1e-6 1e-5 1e-4 1e-3],'XColor','k','YColor','k','ZColor','k','GridColor','k')
figure
set(gcf,'Units','Normalized','OuterPosition',[0 0 1 0.6])
semilogx(f,rad2deg(unwrap(angle(Hs))),'r','linewidth', 2), hold on
semilogx(f,rad2deg(unwrap(angle(Ha))),'k--','linewidth', 2), hold on
xlabel('$f$ [Hz]')
ylabel('Receptance $\phi$ [$^{\circ}$]')
legend({'$H_s$','$H_a$'},'Location','northeast','fontsize',lgndsize)
grid on, grid minor
axis([1e-1 Fs/2 -315 5])
set(gca,'fontsize',txtsize,'Xtick',[1e-2 1e-1 1e0 1e1 1e2 1e3],'XColor','k','YColor','k','ZColor','k','GridColor','k')
axes('Position',[.2 .3 .3 .4]); box on
semilogx(f,rad2deg(unwrap(angle(Hs))),'r','linewidth',2), hold on, grid on
semilogx(f,rad2deg(unwrap(angle(Ha))),'k--','linewidth', 2), hold on
xlim([5 25]); ylim([-315 5]);
grid on, grid minor
set(gca,'fontsize',txtsize-8,'Xtick',[1e-2 1e-1 1e0 1e1 1e2 1e3],'XColor','k','YColor','k','ZColor','k','GridColor','k')
%% EXCITACAO RANDOMICA x(t)
x = randn(1,length(t));
%% CALCULO DE xs E xa
X = fft(x)*dt;
Xs = Hs.*X;
Xa = Ha.*X;
xs = ifft(Xs)*Fs;
xa = ifft(Xa)*Fs;
% PLOTANDO xs, xa
figure
set(gcf,'Units','Normalized','OuterPosition',[0 0 1 0.6])
plot(t,xa,'r','linewidth', 2), hold on
plot(t,xs,'k','linewidth', 2), hold on
xlabel('$t$ [s]','Interpreter','Latex')
ylabel('$x(t)$ [m]','Interpreter','Latex')
legend({'$x_s(t)$','$x_a(t)$'},'location','northeast','fontsize',lgndsize)
xlim([0 50]); ylim([-2.5e-4 2.5e-4])
grid on, grid minor
set(gca,'fontsize',txtsize,'XColor','k','YColor','k','ZColor','k','GridColor','k')
%% ESTIMANDO Ha e Hs
N = round(length(t));
[HsEst, fsEst] = tfestimate(x,xs,[],[],N,Fs);
[HaEst, faEst] = tfestimate(x,xa,[],[],N,Fs);
%% PLOTANDO E COMPARANDO AS FRF's
% ABS(X/F) PELA FREQUENCIA
figure
set(gcf,'Units','Normalized','OuterPosition',[0 0 1 0.6])
subplot(1,2,1) %Hs
loglog(f,abs(Hs),'r','linewidth', 2), hold on
loglog(fsEst,abs(HsEst),'--k','linewidth', 2), hold on
xlabel('$f$ [Hz]')
ylabel('Receptance [m/N]')
legend({'$H_s(j\omega)$ Analytic','$H_s(j\omega)$ Estimated'},'Location','southwest','fontsize',lgndsize)
grid on, grid minor
xlim([1e-1 Fs/2])
set(gca,'fontsize',txtsize,'Ytick',[1e-8 1e-7 1e-6 1e-5 1e-4 1e-3],'XColor','k','YColor','k','ZColor','k','GridColor','k')
axes('Position',[.21 .39 .11 .33]); box on
semilogx(f,abs(Hs),'r','linewidth',2), hold on, grid on
semilogx(fsEst,abs(HsEst),'k--','linewidth', 2), hold on
xlim([10 20]); ylim([2 5]*1e-4);
grid on, grid minor
set(gca,'fontsize',txtsize-8,'Xtick',[1e-2 1e-1 1e0 1e1 1e2 1e3],'XColor','k','YColor','k','ZColor','k','GridColor','k')
subplot(1,2,2) %Ha
loglog(f,abs(Ha),'r','linewidth', 2), hold on
loglog(faEst,abs(HaEst),'--k','linewidth', 2), hold on
xlabel('$f$ [Hz]')
ylabel('Receptance [m/N]')
legend({'$H_a(j\omega)$ Analytic','$H_a(j\omega)$ Estimated'},'Location','southwest','fontsize',lgndsize)
grid on, grid minor
xlim([1e-1 Fs/2])
set(gca,'fontsize',txtsize,'Ytick',[1e-10 1e-9 1e-8 1e-7 1e-6 1e-5 1e-4 1e-3],'XColor','k','YColor','k','ZColor','k','GridColor','k')
axes('Position',[.63 .38 .11 .33]); box on
semilogx(f,abs(Ha),'r','linewidth',2), hold on, grid on
semilogx(faEst,abs(HaEst),'k--','linewidth', 2), hold on
xlim([10 20]); ylim([2 15]*1e-4);
grid on, grid minor
set(gca,'fontsize',txtsize-8,'Xtick',[1e-2 1e-1 1e0 1e1 1e2 1e3],'XColor','k','YColor','k','ZColor','k','GridColor','k')
% ANGULO(X/F) PELA FREQUENCIA
figure
set(gcf,'Units','Normalized','OuterPosition',[0 0 1 0.6])
subplot(1,2,1) %Hs
semilogx(f,rad2deg(unwrap(angle(Hs))),'r','linewidth', 2), hold on
semilogx(fsEst,rad2deg(unwrap(angle(HsEst))),'--k','linewidth', 2), hold on
xlabel('$f$ [Hz]')
ylabel('Receptance $\phi$ [$^{\circ}$]')
legend({'$H_s(j\omega)$ Analytic','$H_s(j\omega)$ Estimated'},'Location','southwest','fontsize',lgndsize)
grid on, grid minor
axis([1e-1 1e2 -200 10])
set(gca,'fontsize',txtsize,'XColor','k','YColor','k','ZColor','k','GridColor','k')
subplot(1,2,2) %Ha
semilogx(f,rad2deg(unwrap(angle(Ha))),'r','linewidth', 2), hold on
semilogx(faEst,rad2deg(unwrap(angle(HaEst))),'--k','linewidth', 2), hold on
xlabel('$f$ [Hz]')
ylabel('Receptance $\phi$ [$^{\circ}$]')
legend({'$H_a(j\omega)$ Analytic','$H_a(j\omega)$ Estimated'},'Location','southwest','fontsize',lgndsize)
grid on, grid minor
axis([1e-1 80 -350 10])
set(gca,'fontsize',txtsize,'XColor','k','YColor','k','ZColor','k','GridColor','k')