398 lines
11 KiB
Matlab
398 lines
11 KiB
Matlab
%% Clear Workspace and Close figures
|
|
clear; close all; clc;
|
|
|
|
%% Intialize Laplace variable
|
|
s = zpk('s');
|
|
|
|
% Numerical Values for the NASS
|
|
% Let's define the parameters for the NASS.
|
|
|
|
mlight = 35; % Mass for light sample [kg]
|
|
mheavy = 85; % Mass for heavy sample [kg]
|
|
|
|
wlight = 2*pi; % Max rot. speed for light sample [rad/s]
|
|
wheavy = 2*pi/60; % Max rot. speed for heavy sample [rad/s]
|
|
|
|
kvc = 1e3; % Voice Coil Stiffness [N/m]
|
|
kpz = 1e8; % Piezo Stiffness [N/m]
|
|
|
|
wdot = 1; % Maximum rotation acceleration [rad/s2]
|
|
|
|
d = 0.01; % Maximum excentricity from rotational axis [m]
|
|
ddot = 0.2; % Maximum Horizontal speed [m/s]
|
|
|
|
save('./mat/parameters.mat');
|
|
|
|
labels = {'Light sample mass [kg]', ...
|
|
'Heavy sample mass [kg]', ...
|
|
'Max rot. speed - light [rpm]', ...
|
|
'Max rot. speed - heavy [rpm]', ...
|
|
'Voice Coil Stiffness [N/m]', ...
|
|
'Piezo Stiffness [N/m]', ...
|
|
'Max rot. acceleration [rad/s2]', ...
|
|
'Max mass excentricity [m]', ...
|
|
'Max Horizontal speed [m/s]'};
|
|
data = [mlight, mheavy, 60*wlight/2/pi, 60*wheavy/2/pi, kvc, kpz, wdot, d, ddot];
|
|
data2orgtable(data', labels, {}, ' %.1e ')
|
|
|
|
% Euler and Coriolis forces - Numerical Result
|
|
% First we will determine the value for Euler and Coriolis forces during regular experiment.
|
|
% - *Euler forces*: $m d_v \ddot{\theta}$
|
|
% - *Coriolis forces*: $2 m \dot{d_v} \dot{\theta}$
|
|
|
|
|
|
Felight = mlight*d*wdot;
|
|
Feheavy = mheavy*d*wdot;
|
|
|
|
Fclight = 2*mlight*ddot*wlight;
|
|
Fcheavy = 2*mheavy*ddot*wheavy;
|
|
|
|
|
|
|
|
% The obtained values are displayed in table [[tab:euler_coriolis]].
|
|
|
|
|
|
data = [Fclight, Fcheavy ;
|
|
Felight, Feheavy];
|
|
data2orgtable(data, {'Coriolis', 'Euler'}, {'Light', 'Heavy'}, ' %.1fN ')
|
|
|
|
% Negative Spring Effect - Numerical Result
|
|
% The negative stiffness due to the rotation is equal to $-m{\omega_0}^2$.
|
|
|
|
|
|
Klight = mlight*wlight^2;
|
|
Kheavy = mheavy*wheavy^2;
|
|
|
|
|
|
|
|
% The values for the negative spring effect are displayed in table [[tab:negative_spring]].
|
|
|
|
% This is definitely negligible when using piezoelectric actuators. It may not be the case when using voice coil actuators.
|
|
|
|
|
|
data = [Klight, Kheavy];
|
|
data2orgtable(data, {'Neg. Spring'}, {'Light', 'Heavy'}, ' %.1f[N/m] ')
|
|
|
|
% Numerical Analysis
|
|
% We plot on the same graph $\frac{|-m \omega^2 + (k - m {\omega_0}^2)|}{|2 m \omega_0 \omega|}$ for the voice coil and the piezo:
|
|
% - with the light sample (figure [[fig:coupling_light]]).
|
|
% - with the heavy sample (figure [[fig:coupling_heavy]]).
|
|
|
|
|
|
f = logspace(-1, 3, 1000);
|
|
|
|
figure;
|
|
hold on;
|
|
plot(f, abs(-mlight*(2*pi*f).^2 + kvc - mlight * wlight^2)./abs(2*mlight*wlight*2*pi*f), 'DisplayName', 'Voice Coil')
|
|
plot(f, abs(-mlight*(2*pi*f).^2 + kpz - mlight * wlight^2)./abs(2*mlight*wlight*2*pi*f), 'DisplayName', 'Piezo')
|
|
plot(f, ones(1, length(f)), 'k--', 'HandleVisibility', 'off')
|
|
hold off;
|
|
xlim([f(1), f(end)]);
|
|
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
|
|
xlabel('Frequency [Hz]');
|
|
legend('Location', 'northeast');
|
|
|
|
|
|
|
|
% #+LABEL: fig:coupling_light
|
|
% #+CAPTION: Relative Coupling for light mass and high rotation speed
|
|
% [[file:./figs/coupling_light.png]]
|
|
|
|
|
|
figure;
|
|
hold on;
|
|
plot(f, abs(-mheavy*(2*pi*f).^2 + kvc - mheavy * wheavy^2)./abs(2*mheavy*wheavy*2*pi*f), 'DisplayName', 'Voice Coil')
|
|
plot(f, abs(-mheavy*(2*pi*f).^2 + kpz - mheavy * wheavy^2)./abs(2*mheavy*wheavy*2*pi*f), 'DisplayName', 'Piezo')
|
|
plot(f, ones(1, length(f)), 'k--', 'HandleVisibility', 'off')
|
|
hold off;
|
|
xlim([f(1), f(end)]);
|
|
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
|
|
xlabel('Frequency [Hz]');
|
|
legend('Location', 'northeast');
|
|
|
|
% Limitations due to negative stiffness effect
|
|
% If $\max{\dot{\theta}} \ll \sqrt{\frac{k}{m}}$, then the negative spring effect is negligible and $k - m\dot{\theta}^2 \approx k$.
|
|
|
|
% Let's estimate what is the maximum rotation speed for which the negative stiffness effect is still negligible ($\omega_\text{max} = 0.1 \sqrt{\frac{k}{m}}$). Results are shown table [[tab:negative_stiffness]].
|
|
|
|
data = 0.1*60*(1/2/pi)*[sqrt(kvc/mlight), sqrt(kpz/mlight);
|
|
sqrt(kvc/mheavy), sqrt(kpz/mheavy)];
|
|
data2orgtable(data, {'Light', 'Heavy'}, {'Voice Coil', 'Piezo'}, ' %.0f[rpm] ')
|
|
|
|
|
|
|
|
% #+NAME: tab:negative_stiffness
|
|
% #+CAPTION: Maximum rotation speed at which negative stiffness is negligible ($0.1\sqrt{\frac{k}{m}}$)
|
|
% #+RESULTS:
|
|
% | | Voice Coil | Piezo |
|
|
% |-------+------------+-----------|
|
|
% | Light | 5[rpm] | 1614[rpm] |
|
|
% | Heavy | 3[rpm] | 1036[rpm] |
|
|
|
|
% The negative spring effect is proportional to the rotational speed $\omega$.
|
|
% The system dynamics will be much more affected when using soft actuator.
|
|
|
|
% #+begin_important
|
|
% Negative stiffness effect has very important effect when using soft actuators.
|
|
% #+end_important
|
|
|
|
% The system can even goes unstable when $m \omega^2 > k$, that is when the centrifugal forces are higher than the forces due to stiffness.
|
|
|
|
% From this analysis, we can determine the lowest practical stiffness that is possible to use: $k_\text{min} = 10 m \omega^2$ (table sec:tab:min_k)
|
|
|
|
|
|
data = 10*[mlight*2*pi, mheavy*2*pi/60]
|
|
data2orgtable(data, {'k min [N/m]'}, {'Light', 'Heavy'}, ' %.0f ')
|
|
|
|
% Voice coil actuator
|
|
|
|
m = mlight;
|
|
k = kvc;
|
|
|
|
ws = linspace(0, 2*pi, 5); % Rotation speed vector [rad/s]
|
|
|
|
Gs = {zeros(1, length(ws))};
|
|
Gcs = {zeros(1, length(ws))};
|
|
|
|
for i = 1:length(ws)
|
|
w = ws(i);
|
|
Gs(i) = {(m*s^2 + (k-m*w^2))/((m*s^2 + (k - m*w^2))^2 + (2*m*w*s)^2)};
|
|
Gcs(i) = {(2*m*w*s)/((m*s^2 + (k - m*w^2))^2 + (2*m*w*s)^2)};
|
|
end
|
|
|
|
freqs = logspace(-2, 1, 1000);
|
|
|
|
figure;
|
|
ax1 = subplot(2,1,1);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, abs(squeeze(freqresp(Gs{i}, freqs, 'Hz'))));
|
|
end
|
|
hold off;
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
set(gca, 'XTickLabel',[]);
|
|
ylabel('Magnitude [m/N]');
|
|
|
|
ax2 = subplot(2,1,2);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
|
|
end
|
|
hold off;
|
|
yticks(-180:90:180);
|
|
ylim([-180 180]);
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
legend('Location', 'southwest');
|
|
linkaxes([ax1,ax2],'x');
|
|
|
|
|
|
|
|
% #+LABEL: fig:G_ws_vc
|
|
% #+CAPTION: Bode plot of the direct transfer function term (from $F_u$ to $D_u$) for multiple rotation speed - Voice coil
|
|
% [[file:figs/G_ws_vc.png]]
|
|
|
|
|
|
|
|
freqs = logspace(-2, 1, 1000);
|
|
|
|
figure;
|
|
ax1 = subplot(2,1,1);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, abs(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))));
|
|
end
|
|
hold off;
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
set(gca, 'XTickLabel',[]);
|
|
ylabel('Magnitude [m/N]');
|
|
|
|
ax2 = subplot(2,1,2);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, 180/pi*angle(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
|
|
end
|
|
hold off;
|
|
yticks(-180:90:180);
|
|
ylim([-180 180]);
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
legend('Location', 'southwest');
|
|
linkaxes([ax1,ax2],'x');
|
|
|
|
% Piezoelectric actuator
|
|
|
|
m = mlight;
|
|
k = kpz;
|
|
|
|
ws = linspace(0, 2*pi, 5); % Rotation speed vector [rad/s]
|
|
|
|
Gs = {zeros(1, length(ws))};
|
|
Gcs = {zeros(1, length(ws))};
|
|
|
|
for i = 1:length(ws)
|
|
w = ws(i);
|
|
Gs(i) = {(m*s^2 + (k-m*w^2))/((m*s^2 + (k - m*w^2))^2 + (2*m*w*s)^2)};
|
|
Gcs(i) = {(2*m*w*s)/((m*s^2 + (k - m*w^2))^2 + (2*m*w*s)^2)};
|
|
end
|
|
|
|
freqs = logspace(2, 3, 1000);
|
|
|
|
figure;
|
|
ax1 = subplot(2,1,1);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, abs(squeeze(freqresp(Gs{i}, freqs, 'Hz'))));
|
|
end
|
|
hold off;
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
set(gca, 'XTickLabel',[]);
|
|
ylabel('Magnitude [m/N]');
|
|
|
|
ax2 = subplot(2,1,2);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
|
|
end
|
|
hold off;
|
|
yticks(-180:90:180);
|
|
ylim([-180 180]);
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
legend('Location', 'southwest');
|
|
linkaxes([ax1,ax2],'x');
|
|
|
|
|
|
|
|
% #+LABEL: fig:G_ws_pz
|
|
% #+CAPTION: Bode plot of the direct transfer function term (from $F_u$ to $D_u$) for multiple rotation speed - Piezoelectric actuator
|
|
% [[file:figs/G_ws_pz.png]]
|
|
|
|
|
|
figure;
|
|
ax1 = subplot(2,1,1);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, abs(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))));
|
|
end
|
|
hold off;
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
set(gca, 'XTickLabel',[]);
|
|
ylabel('Magnitude [m/N]');
|
|
|
|
ax2 = subplot(2,1,2);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, 180/pi*angle(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
|
|
end
|
|
hold off;
|
|
yticks(-180:90:180);
|
|
ylim([-180 180]);
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
legend('Location', 'southwest');
|
|
linkaxes([ax1,ax2],'x');
|
|
|
|
% Campbell diagram
|
|
% The poles of the system are computed for multiple values of the rotation frequency. To simplify the computation of the poles, we add some damping to the system.
|
|
|
|
|
|
m = mlight;
|
|
k = kvc;
|
|
c = 0.1*sqrt(k*m);
|
|
|
|
wsvc = linspace(0, 10, 100); % [rad/s]
|
|
|
|
polesvc = zeros(2, length(wsvc));
|
|
|
|
for i = 1:length(wsvc)
|
|
polei = pole(1/((m*s^2 + c*s + (k - m*wsvc(i)^2))^2 + (2*m*wsvc(i)*s)^2));
|
|
polesvc(:, i) = sort(polei(imag(polei) > 0));
|
|
end
|
|
|
|
m = mlight;
|
|
k = kpz;
|
|
c = 0.1*sqrt(k*m);
|
|
|
|
wspz = linspace(0, 1000, 100); % [rad/s]
|
|
|
|
polespz = zeros(2, length(wspz));
|
|
|
|
for i = 1:length(wspz)
|
|
polei = pole(1/((m*s^2 + c*s + (k - m*wspz(i)^2))^2 + (2*m*wspz(i)*s)^2));
|
|
polespz(:, i) = sort(polei(imag(polei) > 0));
|
|
end
|
|
|
|
|
|
|
|
% We then plot the real and imaginary part of the poles as a function of the rotation frequency (figures [[fig:poles_w_vc]] and [[fig:poles_w_pz]]).
|
|
|
|
% When the real part of one pole becomes positive, the system goes unstable.
|
|
|
|
% For the voice coil (figure [[fig:poles_w_vc]]), the system is unstable when the rotation speed is above 5 rad/s. The real and imaginary part of the poles of the system with piezoelectric actuators are changing much less (figure [[fig:poles_w_pz]]).
|
|
|
|
|
|
figure;
|
|
% Amplitude
|
|
ax1 = subplot(1,2,1);
|
|
hold on;
|
|
plot(wsvc, real(polesvc(1, :)))
|
|
set(gca,'ColorOrderIndex',1)
|
|
plot(wsvc, real(polesvc(2, :)))
|
|
plot(wsvc, zeros(size(wsvc)), 'k--')
|
|
hold off;
|
|
xlabel('Rotation Frequency [rad/s]');
|
|
ylabel('Pole Real Part');
|
|
ax2 = subplot(1,2,2);
|
|
hold on;
|
|
plot(wsvc, imag(polesvc(1, :)))
|
|
set(gca,'ColorOrderIndex',1)
|
|
plot(wsvc, -imag(polesvc(1, :)))
|
|
set(gca,'ColorOrderIndex',1)
|
|
plot(wsvc, imag(polesvc(2, :)))
|
|
set(gca,'ColorOrderIndex',1)
|
|
plot(wsvc, -imag(polesvc(2, :)))
|
|
hold off;
|
|
xlabel('Rotation Frequency [rad/s]');
|
|
ylabel('Pole Imaginary Part');
|
|
|
|
|
|
|
|
% #+LABEL: fig:poles_w_vc
|
|
% #+CAPTION: Real and Imaginary part of the poles of the system as a function of the rotation speed - Voice Coil and light sample
|
|
% [[file:figs/poles_w_vc.png]]
|
|
|
|
|
|
|
|
figure;
|
|
% Amplitude
|
|
ax1 = subplot(1,2,1);
|
|
hold on;
|
|
plot(wspz, real(polespz(1, :)))
|
|
set(gca,'ColorOrderIndex',1)
|
|
plot(wspz, real(polespz(2, :)))
|
|
plot(wspz, zeros(size(wspz)), 'k--')
|
|
hold off;
|
|
xlabel('Rotation Frequency [rad/s]');
|
|
ylabel('Pole Real Part');
|
|
ax2 = subplot(1,2,2);
|
|
hold on;
|
|
plot(wspz, imag(polespz(1, :)))
|
|
set(gca,'ColorOrderIndex',1)
|
|
plot(wspz, -imag(polespz(1, :)))
|
|
set(gca,'ColorOrderIndex',1)
|
|
plot(wspz, imag(polespz(2, :)))
|
|
set(gca,'ColorOrderIndex',1)
|
|
plot(wspz, -imag(polespz(2, :)))
|
|
hold off;
|
|
xlabel('Rotation Frequency [rad/s]');
|
|
ylabel('Pole Imaginary Part');
|