Thomas Dehaeze
6e3677eb29
Folder name is changed, rework the html templates Change the organisation.
758 lines
26 KiB
Org Mode
758 lines
26 KiB
Org Mode
#+TITLE: Spindle Analysis
|
|
#+SETUPFILE: ../config.org
|
|
|
|
The report made by the PEL is accessible [[file:documents/Spindle_report_test.pdf][here]].
|
|
|
|
* Notes
|
|
#+name: fig:setup_spindle
|
|
#+caption: Measurement setup at the PEL lab
|
|
[[file:./img/setup_spindle.png]]
|
|
|
|
* Data Processing
|
|
:PROPERTIES:
|
|
:header-args:matlab+: :tangle matlab/spindle_data_processing.m
|
|
:header-args:matlab+: :comments org :mkdirp yes
|
|
:END:
|
|
<<sec:spindle_data_processing>>
|
|
|
|
#+begin_src bash :exports none :results none
|
|
if [ matlab/spindle_data_processing.m -nt data/spindle_data_processing.zip ]; then
|
|
cp matlab/spindle_data_processing.m spindle_data_processing.m;
|
|
zip data/spindle_data_processing \
|
|
mat/10turns_1rpm_icepap.txt \
|
|
mat/10turns_60rpm_IcepapFIR.txt \
|
|
src/getAsynchronousError.m \
|
|
spindle_data_processing.m
|
|
rm spindle_data_processing.m;
|
|
fi
|
|
#+end_src
|
|
|
|
#+begin_note
|
|
All the files (data and Matlab scripts) are accessible [[file:data/spindle_data_processing.zip][here]].
|
|
#+end_note
|
|
|
|
** Matlab Init :noexport:ignore:
|
|
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
|
|
<<matlab-dir>>
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none :results silent :noweb yes
|
|
<<matlab-init>>
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none :results silent :noweb yes
|
|
%% Add the getAsynchronousError to path
|
|
addpath('./src/');
|
|
#+end_src
|
|
|
|
** Load Measurement Data
|
|
#+begin_src matlab :results none :exports code
|
|
spindle_1rpm_table = readtable('./mat/10turns_1rpm_icepap.txt');
|
|
spindle_60rpm_table = readtable('./mat/10turns_60rpm_IcepapFIR.txt');
|
|
#+end_src
|
|
|
|
#+begin_src matlab :results output :exports code
|
|
spindle_1rpm_table(1, :)
|
|
#+end_src
|
|
|
|
#+begin_src matlab :results none :exports code
|
|
spindle_1rpm = table2array(spindle_1rpm_table);
|
|
spindle_60rpm = table2array(spindle_60rpm_table);
|
|
#+end_src
|
|
|
|
** Convert Signals from [deg] to [sec]
|
|
#+begin_src matlab :results none :exports code
|
|
speed_1rpm = 360/60; % [deg/sec]
|
|
spindle_1rpm(:, 1) = spindle_1rpm(:, 2)/speed_1rpm; % From position [deg] to time [s]
|
|
|
|
speed_60rpm = 360/1; % [deg/sec]
|
|
spindle_60rpm(:, 1) = spindle_60rpm(:, 2)/speed_60rpm; % From position [deg] to time [s]
|
|
#+end_src
|
|
|
|
** Convert Signals
|
|
#+begin_src matlab :results none :exports code
|
|
% scaling = 1/80000; % 80 mV/um
|
|
scaling = 1e-6; % [um] to [m]
|
|
|
|
spindle_1rpm(:, 3:end) = scaling*spindle_1rpm(:, 3:end); % [V] to [m]
|
|
spindle_1rpm(:, 3:end) = spindle_1rpm(:, 3:end)-mean(spindle_1rpm(:, 3:end)); % Remove mean
|
|
|
|
spindle_60rpm(:, 3:end) = scaling*spindle_60rpm(:, 3:end); % [V] to [m]
|
|
spindle_60rpm(:, 3:end) = spindle_60rpm(:, 3:end)-mean(spindle_60rpm(:, 3:end)); % Remove mean
|
|
#+end_src
|
|
|
|
** Ts and Fs for both measurements
|
|
#+begin_src matlab :results none :exports code
|
|
Ts_1rpm = spindle_1rpm(end, 1)/(length(spindle_1rpm(:, 1))-1);
|
|
Fs_1rpm = 1/Ts_1rpm;
|
|
|
|
Ts_60rpm = spindle_60rpm(end, 1)/(length(spindle_60rpm(:, 1))-1);
|
|
Fs_60rpm = 1/Ts_60rpm;
|
|
#+end_src
|
|
|
|
** Find Noise of the ADC [$\frac{m}{\sqrt{Hz}}$]
|
|
#+begin_src matlab :results none :exports code
|
|
data = spindle_1rpm(:, 5);
|
|
dV_1rpm = min(abs(data(1) - data(data ~= data(1))));
|
|
noise_1rpm = dV_1rpm/sqrt(12*Fs_1rpm/2);
|
|
|
|
data = spindle_60rpm(:, 5);
|
|
dV_60rpm = min(abs(data(50) - data(data ~= data(50))));
|
|
noise_60rpm = dV_60rpm/sqrt(12*Fs_60rpm/2);
|
|
#+end_src
|
|
|
|
** Save all the data under spindle struct
|
|
#+begin_src matlab :results none :exports code
|
|
spindle.rpm1.time = spindle_1rpm(:, 1);
|
|
spindle.rpm1.deg = spindle_1rpm(:, 2);
|
|
spindle.rpm1.Ts = Ts_1rpm;
|
|
spindle.rpm1.Fs = 1/Ts_1rpm;
|
|
spindle.rpm1.x = spindle_1rpm(:, 3);
|
|
spindle.rpm1.y = spindle_1rpm(:, 4);
|
|
spindle.rpm1.z = spindle_1rpm(:, 5);
|
|
spindle.rpm1.adcn = noise_1rpm;
|
|
|
|
spindle.rpm60.time = spindle_60rpm(:, 1);
|
|
spindle.rpm60.deg = spindle_60rpm(:, 2);
|
|
spindle.rpm60.Ts = Ts_60rpm;
|
|
spindle.rpm60.Fs = 1/Ts_60rpm;
|
|
spindle.rpm60.x = spindle_60rpm(:, 3);
|
|
spindle.rpm60.y = spindle_60rpm(:, 4);
|
|
spindle.rpm60.z = spindle_60rpm(:, 5);
|
|
spindle.rpm60.adcn = noise_60rpm;
|
|
#+end_src
|
|
|
|
** Compute Asynchronous data
|
|
#+begin_src matlab :results none :exports code
|
|
for direction = {'x', 'y', 'z'}
|
|
spindle.rpm1.([direction{1}, 'async']) = getAsynchronousError(spindle.rpm1.(direction{1}), 10);
|
|
spindle.rpm60.([direction{1}, 'async']) = getAsynchronousError(spindle.rpm60.(direction{1}), 10);
|
|
end
|
|
#+end_src
|
|
|
|
** Save data
|
|
#+begin_src matlab
|
|
save('./mat/spindle_data.mat', 'spindle');
|
|
#+end_src
|
|
|
|
* Time Domain Data
|
|
:PROPERTIES:
|
|
:header-args:matlab+: :tangle matlab/spindle_time_domain.m
|
|
:header-args:matlab+: :comments org :mkdirp yes
|
|
:END:
|
|
<<sec:spindle_time_domain>>
|
|
|
|
#+begin_src bash :exports none :results none
|
|
if [ matlab/spindle_time_domain.m -nt data/spindle_time_domain.zip ]; then
|
|
cp matlab/spindle_time_domain.m spindle_time_domain.m;
|
|
zip data/spindle_time_domain \
|
|
mat/spindle_data.mat \
|
|
spindle_time_domain.m
|
|
rm spindle_time_domain.m;
|
|
fi
|
|
#+end_src
|
|
|
|
#+begin_note
|
|
All the files (data and Matlab scripts) are accessible [[file:data/spindle_time_domain.zip][here]].
|
|
#+end_note
|
|
|
|
** Matlab Init :noexport:ignore:
|
|
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
|
|
<<matlab-dir>>
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none :results none :noweb yes
|
|
<<matlab-init>>
|
|
#+end_src
|
|
|
|
** Load the processed data
|
|
#+begin_src matlab
|
|
load('./mat/spindle_data.mat', 'spindle');
|
|
#+end_src
|
|
|
|
** Plot X-Y-Z position with respect to Time - 1rpm
|
|
#+begin_src matlab :results none :exports code
|
|
figure;
|
|
hold on;
|
|
plot(spindle.rpm1.time, spindle.rpm1.x);
|
|
plot(spindle.rpm1.time, spindle.rpm1.y);
|
|
plot(spindle.rpm1.time, spindle.rpm1.z);
|
|
hold off;
|
|
xlabel('Time [s]'); ylabel('Amplitude [m]');
|
|
legend({'tx - 1rpm', 'ty - 1rpm', 'tz - 1rpm'});
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_xyz_1rpm
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/spindle_xyz_1rpm.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_xyz_1rpm
|
|
#+CAPTION: Raw time domain translation - 1rpm
|
|
#+RESULTS: fig:spindle_xyz_1rpm
|
|
[[file:figs/spindle_xyz_1rpm.png]]
|
|
|
|
** Plot X-Y-Z position with respect to Time - 60rpm
|
|
The measurements for the spindle turning at 60rpm are shown figure [[fig:spindle_xyz_60rpm]].
|
|
|
|
#+begin_src matlab :results none :exports code
|
|
figure;
|
|
hold on;
|
|
plot(spindle.rpm60.time, spindle.rpm60.x);
|
|
plot(spindle.rpm60.time, spindle.rpm60.y);
|
|
plot(spindle.rpm60.time, spindle.rpm60.z);
|
|
hold off;
|
|
xlabel('Time [s]'); ylabel('Amplitude [m]');
|
|
legend({'tx - 60rpm', 'ty - 60rpm', 'tz - 60rpm'});
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_xyz_60rpm
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/spindle_xyz_60rpm.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_xyz_60rpm
|
|
#+CAPTION: Raw time domain translation - 60rpm
|
|
#+RESULTS: fig:spindle_xyz_60rpm
|
|
[[file:figs/spindle_xyz_60rpm.png]]
|
|
|
|
** Plot Synchronous and Asynchronous - 1rpm
|
|
#+begin_src matlab :results none :exports code
|
|
figure;
|
|
hold on;
|
|
plot(spindle.rpm1.time, spindle.rpm1.x);
|
|
plot(spindle.rpm1.time, spindle.rpm1.xasync);
|
|
hold off;
|
|
xlabel('Time [s]'); ylabel('Amplitude [m]');
|
|
legend({'tx - 1rpm - Sync', 'tx - 1rpm - Async'});
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_1rpm_sync_async
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/spindle_1rpm_sync_async.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_1rpm_sync_async
|
|
#+CAPTION: Comparison of the synchronous and asynchronous displacements - 1rpm
|
|
#+RESULTS: fig:spindle_1rpm_sync_async
|
|
[[file:figs/spindle_1rpm_sync_async.png]]
|
|
|
|
** Plot Synchronous and Asynchronous - 60rpm
|
|
The data is split into its Synchronous and Asynchronous part (figure [[fig:spindle_60rpm_sync_async]]). We then use the Asynchronous part for the analysis in the following sections as we suppose that we can deal with the synchronous part with feedforward control.
|
|
|
|
#+begin_src matlab :results none :exports code
|
|
figure;
|
|
hold on;
|
|
plot(spindle.rpm60.time, spindle.rpm60.x);
|
|
plot(spindle.rpm60.time, spindle.rpm60.xasync);
|
|
hold off;
|
|
xlabel('Time [s]'); ylabel('Amplitude [m]');
|
|
legend({'tx - 60rpm - Sync', 'tx - 60rpm - Async'});
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_60rpm_sync_async
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/spindle_60rpm_sync_async.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_60rpm_sync_async
|
|
#+CAPTION: Comparison of the synchronous and asynchronous displacements - 60rpm
|
|
#+RESULTS: fig:spindle_60rpm_sync_async
|
|
[[file:figs/spindle_60rpm_sync_async.png]]
|
|
|
|
** Plot X against Y
|
|
#+begin_src matlab :results none :exports code
|
|
figure;
|
|
hold on;
|
|
plot(spindle.rpm1.x, spindle.rpm1.y);
|
|
plot(spindle.rpm60.x, spindle.rpm60.y);
|
|
hold off;
|
|
xlabel('X Amplitude [m]'); ylabel('Y Amplitude [m]');
|
|
legend({'1rpm', '60rpm'});
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_xy_1_60rpm
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/spindle_xy_1_60rpm.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_xy_1_60rpm
|
|
#+CAPTION: Synchronous x-y displacement
|
|
#+RESULTS: fig:spindle_xy_1_60rpm
|
|
[[file:figs/spindle_xy_1_60rpm.png]]
|
|
|
|
** Plot X against Y - Asynchronous
|
|
#+begin_src matlab :results none :exports code
|
|
figure;
|
|
hold on;
|
|
plot(spindle.rpm1.xasync, spindle.rpm1.yasync);
|
|
plot(spindle.rpm60.xasync, spindle.rpm60.yasync);
|
|
hold off;
|
|
xlabel('X Amplitude [m]'); ylabel('Y Amplitude [m]');
|
|
legend({'1rpm', '60rpm'});
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_xy_1_60_rpm_async
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/spindle_xy_1_60_rpm_async.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_xy_1_60_rpm_async
|
|
#+CAPTION: Asynchronous x-y displacement
|
|
#+RESULTS: fig:spindle_xy_1_60_rpm_async
|
|
[[file:figs/spindle_xy_1_60_rpm_async.png]]
|
|
|
|
* Model of the spindle
|
|
:PROPERTIES:
|
|
:header-args:matlab+: :tangle matlab/spindle_model.m
|
|
:header-args:matlab+: :comments org :mkdirp yes
|
|
:END:
|
|
<<sec:spindle_model>>
|
|
|
|
#+begin_src bash :exports none :results none
|
|
if [ matlab/spindle_model.m -nt data/spindle_model.zip ]; then
|
|
cp matlab/spindle_model.m spindle_model.m;
|
|
zip data/spindle_model \
|
|
spindle_model.m
|
|
rm spindle_model.m;
|
|
fi
|
|
#+end_src
|
|
|
|
#+begin_note
|
|
All the files (data and Matlab scripts) are accessible [[file:data/spindle_model.zip][here]].
|
|
#+end_note
|
|
|
|
** Schematic of the model
|
|
The model of the spindle used is shown figure [[fig:model_spindle]].
|
|
|
|
$f$ is the perturbation force of the spindle and $d$ is the measured displacement.
|
|
|
|
#+name: fig:model_spindle
|
|
#+caption: Model of the Spindle
|
|
#+attr_latex: :float t
|
|
[[./figs/uniaxial-model-spindle.png]]
|
|
|
|
** Matlab Init :noexport:ignore:
|
|
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
|
|
<<matlab-dir>>
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none :results silent :noweb yes
|
|
<<matlab-init>>
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none :results silent :noweb yes
|
|
%% Add the getAsynchronousError to path
|
|
addpath('./src/');
|
|
#+end_src
|
|
|
|
** Parameters
|
|
#+begin_src matlab :exports code :results none
|
|
mg = 3000; % Mass of granite [kg]
|
|
ms = 50; % Mass of Spindle [kg]
|
|
|
|
kg = 1e8; % Stiffness of granite [N/m]
|
|
ks = 5e7; % Stiffness of spindle [N/m]
|
|
#+end_src
|
|
|
|
** Compute Mass and Stiffness Matrices
|
|
#+begin_src matlab :exports code :results none
|
|
Mm = diag([ms, mg]);
|
|
Km = diag([ks, ks+kg]) - diag(ks, -1) - diag(ks, 1);
|
|
#+end_src
|
|
|
|
** Compute resonance frequencies
|
|
#+begin_src matlab :exports code :results none
|
|
A = [zeros(size(Mm)) eye(size(Mm)) ; -Mm\Km zeros(size(Mm))];
|
|
eigA = imag(eigs(A))/2/pi;
|
|
eigA = eigA(eigA>0);
|
|
eigA = eigA(1:2);
|
|
#+end_src
|
|
|
|
** From model_damping compute the Damping Matrix
|
|
#+begin_src matlab :exports code :results none
|
|
modal_damping = 1e-5;
|
|
|
|
ab = [0.5*eigA(1) 0.5/eigA(1) ; 0.5*eigA(2) 0.5/eigA(2)]\[modal_damping ; modal_damping];
|
|
|
|
Cm = ab(1)*Mm +ab(2)*Km;
|
|
#+end_src
|
|
|
|
** Define inputs, outputs and state names
|
|
#+begin_src matlab :exports code :results none
|
|
StateName = {...
|
|
'xs', ... % Displacement of Spindle [m]
|
|
'xg', ... % Displacement of Granite [m]
|
|
'vs', ... % Velocity of Spindle [m]
|
|
'vg', ... % Velocity of Granite [m]
|
|
};
|
|
StateUnit = {'m', 'm', 'm/s', 'm/s'};
|
|
|
|
InputName = {...
|
|
'f' ... % Spindle Force [N]
|
|
};
|
|
InputUnit = {'N'};
|
|
|
|
OutputName = {...
|
|
'd' ... % Displacement [m]
|
|
};
|
|
OutputUnit = {'m'};
|
|
#+end_src
|
|
|
|
** Define A, B and C matrices
|
|
#+begin_src matlab :exports code :results none
|
|
% A Matrix
|
|
A = [zeros(size(Mm)) eye(size(Mm)) ; ...
|
|
-Mm\Km -Mm\Cm];
|
|
|
|
% B Matrix
|
|
B_low = zeros(length(StateName), length(InputName));
|
|
B_low(strcmp(StateName,'vs'), strcmp(InputName,'f')) = 1;
|
|
B_low(strcmp(StateName,'vg'), strcmp(InputName,'f')) = -1;
|
|
B = blkdiag(zeros(length(StateName)/2), pinv(Mm))*B_low;
|
|
|
|
% C Matrix
|
|
C = zeros(length(OutputName), length(StateName));
|
|
C(strcmp(OutputName,'d'), strcmp(StateName,'xs')) = 1;
|
|
C(strcmp(OutputName,'d'), strcmp(StateName,'xg')) = -1;
|
|
|
|
% D Matrix
|
|
D = zeros(length(OutputName), length(InputName));
|
|
#+end_src
|
|
|
|
** Generate the State Space Model
|
|
#+begin_src matlab :exports code :results none
|
|
sys = ss(A, B, C, D);
|
|
sys.StateName = StateName;
|
|
sys.StateUnit = StateUnit;
|
|
|
|
sys.InputName = InputName;
|
|
sys.InputUnit = InputUnit;
|
|
|
|
sys.OutputName = OutputName;
|
|
sys.OutputUnit = OutputUnit;
|
|
#+end_src
|
|
|
|
** Bode Plot
|
|
The transfer function from a disturbance force $f$ to the measured displacement $d$ is shown figure [[fig:spindle_f_to_d]].
|
|
|
|
#+begin_src matlab :exports code :results none
|
|
freqs = logspace(-1, 3, 1000);
|
|
|
|
figure;
|
|
plot(freqs, abs(squeeze(freqresp(sys('d', 'f'), freqs, 'Hz'))));
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('Amplitude [m/N]');
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_f_to_d
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/spindle_f_to_d.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_f_to_d
|
|
#+CAPTION: Bode plot of the transfer function from $f$ to $d$
|
|
#+RESULTS: fig:spindle_f_to_d
|
|
[[file:figs/spindle_f_to_d.png]]
|
|
|
|
** Save the model
|
|
#+begin_src matlab
|
|
save('./mat/spindle_model.mat', 'sys');
|
|
#+end_src
|
|
|
|
* Frequency Domain Data
|
|
:PROPERTIES:
|
|
:header-args:matlab+: :tangle matlab/spindle_psd.m
|
|
:header-args:matlab+: :comments org :mkdirp yes
|
|
:END:
|
|
<<sec:spindle_psd>>
|
|
|
|
#+begin_src bash :exports none :results none
|
|
if [ matlab/spindle_psd.m -nt data/spindle_psd.zip ]; then
|
|
cp matlab/spindle_psd.m spindle_psd.m;
|
|
zip data/spindle_psd \
|
|
mat/spindle_data.mat \
|
|
mat/spindle_model.mat \
|
|
spindle_psd.m
|
|
rm spindle_psd.m;
|
|
fi
|
|
#+end_src
|
|
|
|
#+begin_note
|
|
All the files (data and Matlab scripts) are accessible [[file:data/spindle_psd.zip][here]].
|
|
#+end_note
|
|
|
|
** Matlab Init :noexport:ignore:
|
|
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
|
|
<<matlab-dir>>
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none :results silent :noweb yes
|
|
<<matlab-init>>
|
|
#+end_src
|
|
|
|
** Load the processed data and the model
|
|
#+begin_src matlab
|
|
load('./mat/spindle_data.mat', 'spindle');
|
|
load('./mat/spindle_model.mat', 'sys');
|
|
#+end_src
|
|
|
|
** Compute the PSD
|
|
#+begin_src matlab :exports code :results none
|
|
n_av = 4; % Number of average
|
|
|
|
[pxx_1rpm, f_1rpm] = pwelch(spindle.rpm1.xasync, hanning(ceil(length(spindle.rpm1.xasync)/n_av)), [], [], spindle.rpm1.Fs);
|
|
[pyy_1rpm, ~] = pwelch(spindle.rpm1.yasync, hanning(ceil(length(spindle.rpm1.yasync)/n_av)), [], [], spindle.rpm1.Fs);
|
|
[pzz_1rpm, ~] = pwelch(spindle.rpm1.zasync, hanning(ceil(length(spindle.rpm1.zasync)/n_av)), [], [], spindle.rpm1.Fs);
|
|
|
|
[pxx_60rpm, f_60rpm] = pwelch(spindle.rpm60.xasync, hanning(ceil(length(spindle.rpm60.xasync)/n_av)), [], [], spindle.rpm60.Fs);
|
|
[pyy_60rpm, ~] = pwelch(spindle.rpm60.yasync, hanning(ceil(length(spindle.rpm60.yasync)/n_av)), [], [], spindle.rpm60.Fs);
|
|
[pzz_60rpm, ~] = pwelch(spindle.rpm60.zasync, hanning(ceil(length(spindle.rpm60.zasync)/n_av)), [], [], spindle.rpm60.Fs);
|
|
#+end_src
|
|
|
|
** Plot the computed PSD
|
|
The Amplitude Spectral Densities of the displacement of the spindle for the $x$, $y$ and $z$ directions are shown figure [[fig:spindle_psd_xyz_60rpm]]. They correspond to the Asynchronous part shown figure [[fig:spindle_60rpm_sync_async]].
|
|
|
|
#+begin_src matlab :exports code :results none
|
|
figure;
|
|
hold on;
|
|
plot(f_1rpm, (pxx_1rpm).^.5, 'DisplayName', '$P_{xx}$ - 1rpm');
|
|
plot(f_1rpm, (pyy_1rpm).^.5, 'DisplayName', '$P_{yy}$ - 1rpm');
|
|
plot(f_1rpm, (pzz_1rpm).^.5, 'DisplayName', '$P_{zz}$ - 1rpm');
|
|
% plot(f_1rpm, spindle.rpm1.adcn*ones(size(f_1rpm)), '--k', 'DisplayName', 'ADC - 1rpm');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
|
|
legend('Location', 'northeast');
|
|
xlim([f_1rpm(2), f_1rpm(end)]);
|
|
#+end_src
|
|
#+NAME: fig:spindle_psd_xyz_1rpm
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/spindle_psd_xyz_1rpm.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_psd_xyz_1rpm
|
|
#+CAPTION: Power spectral density of the Asynchronous displacement - 1rpm
|
|
#+RESULTS: fig:spindle_psd_xyz_1rpm
|
|
[[file:figs/spindle_psd_xyz_1rpm.png]]
|
|
|
|
#+begin_src matlab :exports code :results none
|
|
figure;
|
|
hold on;
|
|
plot(f_60rpm, (pxx_60rpm).^.5, 'DisplayName', '$P_{xx}$ - 60rpm');
|
|
plot(f_60rpm, (pyy_60rpm).^.5, 'DisplayName', '$P_{yy}$ - 60rpm');
|
|
plot(f_60rpm, (pzz_60rpm).^.5, 'DisplayName', '$P_{zz}$ - 60rpm');
|
|
% plot(f_60rpm, spindle.rpm60.adcn*ones(size(f_60rpm)), '--k', 'DisplayName', 'ADC - 60rpm');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
|
|
legend('Location', 'northeast');
|
|
xlim([f_60rpm(2), f_60rpm(end)]);
|
|
#+end_src
|
|
#+NAME: fig:spindle_psd_xyz_60rpm
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/spindle_psd_xyz_60rpm.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_psd_xyz_60rpm
|
|
#+CAPTION: Power spectral density of the Asynchronous displacement - 60rpm
|
|
#+RESULTS: fig:spindle_psd_xyz_60rpm
|
|
[[file:figs/spindle_psd_xyz_60rpm.png]]
|
|
|
|
** Compute the response of the model
|
|
#+begin_src matlab :exports code :results none
|
|
Tfd = abs(squeeze(freqresp(sys('d', 'f'), f_60rpm, 'Hz')));
|
|
#+end_src
|
|
|
|
** Plot the PSD of the Force using the model
|
|
#+begin_src matlab :exports code :results none
|
|
figure;
|
|
plot(f_60rpm, (pxx_60rpm.^.5)./Tfd, 'DisplayName', '$P_{xx}$');
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$N/\sqrt{Hz}$]');
|
|
xlim([f_60rpm(2), f_60rpm(end)]);
|
|
#+end_src
|
|
#+NAME: fig:spindle_psd_f_60rpm
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/spindle_psd_f_60rpm.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_psd_f_60rpm
|
|
#+CAPTION: Power spectral density of the force - 60rpm
|
|
#+RESULTS: fig:spindle_psd_f_60rpm
|
|
[[file:figs/spindle_psd_f_60rpm.png]]
|
|
|
|
** Estimated Shape of the PSD of the force
|
|
#+begin_src matlab :exports code :results none
|
|
s = tf('s');
|
|
|
|
Wd_simple = 1e-8/(1+s/2/pi/0.5)/(1+s/2/pi/100);
|
|
Wf_simple = Wd_simple/tf(sys('d', 'f'));
|
|
TWf_simple = abs(squeeze(freqresp(Wf_simple, f_60rpm, 'Hz')));
|
|
|
|
% Wf = 0.48902*(s+327.9)*(s^2 + 109.6*s + 1.687e04)/((s^2 + 30.59*s + 8541)*(s^2 + 29.11*s + 3.268e04));
|
|
% Wf = 0.15788*(s+418.6)*(s+1697)^2*(s^2 + 124.3*s + 2.529e04)*(s^2 + 681.3*s + 9.018e05)/((s^2 + 23.03*s + 8916)*(s^2 + 33.85*s + 6.559e04)*(s^2 + 71.43*s + 4.283e05)*(s^2 + 40.64*s + 1.789e06));
|
|
|
|
Wf = (s+1697)^2*(s^2 + 114.5*s + 2.278e04)*(s^2 + 205.1*s + 1.627e05)*(s^2 + 285.8*s + 8.624e05)*(s+100)/((s+0.5)*3012*(s^2 + 23.03*s + 8916)*(s^2 + 17.07*s + 4.798e04)*(s^2 + 41.17*s + 4.347e05)*(s^2 + 78.99*s + 1.789e06));
|
|
|
|
TWf = abs(squeeze(freqresp(Wf, f_60rpm, 'Hz')));
|
|
#+end_src
|
|
|
|
** PSD in [N]
|
|
Above $200Hz$, the Amplitude Spectral Density seems dominated by noise coming from the electronics (charge amplifier, ADC, ...). So we don't know what is the frequency content of the force above that frequency. However, we assume that $P_{xx}$ is decreasing with $1/f$ as it seems so be the case below $100Hz$ (figure [[fig:spindle_psd_xyz_60rpm]]).
|
|
|
|
We then fit the PSD of the displacement with a transfer function (figure [[fig:spindle_psd_d_comp_60rpm]]).
|
|
|
|
#+begin_src matlab :exports code :results none
|
|
figure;
|
|
hold on;
|
|
plot(f_60rpm, (pxx_60rpm.^.5)./Tfd, 'DisplayName', '$\sqrt{P_{xx}}/|T_{d/f}|$');
|
|
plot(f_60rpm, TWf, 'DisplayName', 'Wf');
|
|
plot(f_60rpm, TWf_simple, '-k', 'DisplayName', 'Wfs');
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$N/\sqrt{Hz}$]');
|
|
xlim([f_60rpm(2), f_60rpm(end)]);
|
|
legend('Location', 'northeast');
|
|
#+end_src
|
|
#+NAME: fig:spindle_psd_f_comp_60rpm
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/spindle_psd_f_comp_60rpm.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_psd_f_comp_60rpm
|
|
#+CAPTION: Power spectral density of the force - 60rpm
|
|
#+RESULTS: fig:spindle_psd_f_comp_60rpm
|
|
[[file:figs/spindle_psd_f_comp_60rpm.png]]
|
|
|
|
** PSD in [m]
|
|
To obtain the PSD of the force $f$ that induce such displacement, we use the following formula:
|
|
\[ \sqrt{PSD(d)} = |T_{d/f}| \sqrt{PSD(f)} \]
|
|
|
|
And so we have:
|
|
\[ \sqrt{PSD(f)} = |T_{d/f}|^{-1} \sqrt{PSD(d)} \]
|
|
|
|
The obtain Power Spectral Density of the force is displayed figure [[fig:spindle_psd_f_comp_60rpm]].
|
|
|
|
#+begin_src matlab :exports code :results none
|
|
figure;
|
|
hold on;
|
|
plot(f_60rpm, pxx_60rpm.^.5, 'DisplayName', '$\sqrt{P_{xx}}$');
|
|
plot(f_60rpm, TWf.*Tfd, 'DisplayName', '$|W_f|*|T_{d/f}|$');
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
|
|
xlim([f_60rpm(2), f_60rpm(end)]);
|
|
legend('Location', 'northeast');
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_psd_d_comp_60rpm
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/spindle_psd_d_comp_60rpm.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_psd_d_comp_60rpm
|
|
#+CAPTION: Comparison of the power spectral density of the measured displacement and of the model
|
|
#+RESULTS: fig:spindle_psd_d_comp_60rpm
|
|
[[file:figs/spindle_psd_d_comp_60rpm.png]]
|
|
|
|
** Compute the resulting RMS value [m]
|
|
#+begin_src matlab :exports code :results none
|
|
figure;
|
|
hold on;
|
|
plot(f_60rpm, 1e9*cumtrapz(f_60rpm, (pxx_60rpm)).^.5, '--', 'DisplayName', 'Exp. Data');
|
|
plot(f_60rpm, 1e9*cumtrapz(f_60rpm, ((TWf.*Tfd).^2)).^.5, '--', 'DisplayName', 'Estimated');
|
|
hold off;
|
|
set(gca, 'XScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('CPS [$nm$ rms]');
|
|
xlim([f_60rpm(2), f_60rpm(end)]);
|
|
legend('Location', 'southeast');
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_cps_d_comp_60rpm
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/spindle_cps_d_comp_60rpm.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_cps_d_comp_60rpm
|
|
#+CAPTION: Cumulative Power Spectrum - 60rpm
|
|
#+RESULTS: fig:spindle_cps_d_comp_60rpm
|
|
[[file:figs/spindle_cps_d_comp_60rpm.png]]
|
|
|
|
** Compute the resulting RMS value [m]
|
|
#+begin_src matlab :exports code :results none
|
|
figure;
|
|
hold on;
|
|
plot(f_1rpm, 1e9*cumtrapz(f_1rpm, (pxx_1rpm)), '--', 'DisplayName', 'Exp. Data');
|
|
plot(f_1rpm, 1e9*(f_1rpm(end)-f_1rpm(1))/(length(f_1rpm)-1)*cumsum(pxx_1rpm), '--', 'DisplayName', 'Exp. Data');
|
|
hold off;
|
|
set(gca, 'XScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('CPS [$nm$ rms]');
|
|
xlim([f_1rpm(2), f_1rpm(end)]);
|
|
legend('Location', 'southeast');
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_cps_d_comp_1rpm
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/spindle_cps_d_comp_1rpm.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+NAME: fig:spindle_cps_d_comp_1rpm
|
|
#+CAPTION: Cumulative Power Spectrum - 1rpm
|
|
#+RESULTS: fig:spindle_cps_d_comp_1rpm
|
|
[[file:figs/spindle_cps_d_comp_1rpm.png]]
|
|
|
|
* Functions
|
|
** getAsynchronousError
|
|
:PROPERTIES:
|
|
:header-args:matlab+: :tangle src/getAsynchronousError.m
|
|
:header-args:matlab+: :comments org :mkdirp yes
|
|
:END:
|
|
<<sec:getAsynchronousError>>
|
|
|
|
This Matlab function is accessible [[file:src/getAsynchronousError.m][here]].
|
|
|
|
#+begin_src matlab :eval no :results none
|
|
function [Wxdec] = getAsynchronousError(data, NbTurn)
|
|
%%
|
|
L = length(data);
|
|
res_per_rev = L/NbTurn;
|
|
|
|
P = 0:(res_per_rev*NbTurn-1);
|
|
Pos = P' * 360/res_per_rev;
|
|
|
|
% Temperature correction
|
|
x1 = myfit2(Pos, data);
|
|
|
|
% Convert data to frequency domain and scale accordingly
|
|
X2 = 2/(res_per_rev*NbTurn)*fft(x1);
|
|
f2 = (0:L-1)./NbTurn; %upr -> once per revolution
|
|
|
|
%%
|
|
X2dec = zeros(size(X2));
|
|
% Get only the non integer data
|
|
X2dec(mod(f2(:), 1) ~= 0) = X2(mod(f2(:), 1) ~= 0);
|
|
|
|
Wxdec = real((res_per_rev*NbTurn)/2 * ifft(X2dec));
|
|
|
|
%%
|
|
function Y = myfit2(x,y)
|
|
A = [x ones(size(x))]\y;
|
|
a = A(1); b = A(2);
|
|
Y = y - (a*x + b);
|
|
end
|
|
end
|
|
#+end_src
|