tdehaeze/phd-micro-station-modal-analysis
Archived
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Initial Commit

Browse Source
This commit is contained in:
Thomas Dehaeze
2024-03-19 15:11:22 +01:00
commit 1fe38c4ba3
172 changed files with 11389 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

236
matlab/modal_1_meas_setup.m Normal file
View File

@@ -0,0 +1,236 @@
%% Clear Workspace and Close figures
clear; close all; clc;
%% Intialize Laplace variable
s = zpk('s');
%% Path for functions, data and scripts
addpath('./mat/'); % Path for data
%% Colors for the figures
colors = colororder;
% Location of the Accelerometers
% <<ssec:modal_accelerometers>>
% 4 tri-axis accelerometers are used for each solid body.
% Only 2 could have been used as only 6DOF have to be measured, however, we have chosen to have some *redundancy*.
% This could also help us identify measurement problems or flexible modes is present.
% The position of the accelerometers are:
% - 4 on the first granite
% - 4 on the second granite
% - 4 on top of the translation stage (figure ref:fig:accelerometers_ty)
% - 4 on top of the tilt stage
% - 3 on top of the spindle
% - 4 on top of the hexapod (figure ref:fig:accelerometers_hexapod)
% In total, 23 accelerometers are used: *69 DOFs are thus measured*.
% The precise determination of the position of each accelerometer is done using the SolidWorks model (shown on figure ref:fig:location_accelerometers).
% #+name: fig:accelerometer_pictures
% #+caption: Accelerometers fixed on the micro-station
% #+begin_figure
% #+attr_latex: :caption \subcaption{\label{fig:accelerometers_ty}$T_y$ stage}
% #+attr_latex: :options {0.49\textwidth}
% #+begin_subfigure
% #+attr_latex: :height 6cm
% [[file:figs/accelerometers_ty.jpg]]
% #+end_subfigure
% #+attr_latex: :caption \subcaption{\label{fig:accelerometers_hexapod}Micro-Hexapod}
% #+attr_latex: :options {0.49\textwidth}
% #+begin_subfigure
% #+attr_latex: :height 6cm
% [[file:figs/accelerometers_hexapod.jpg]]
% #+end_subfigure
% #+end_figure
% #+name: fig:location_accelerometers
% #+caption: Position of the accelerometers using SolidWorks
% #+attr_latex: :width \linewidth
% [[file:figs/location_accelerometers.png]]
% The precise position of all the 23 accelerometer with respect to a frame located at the point of interest (located 270mm above the top platform of the hexapod) are shown in table ref:tab:position_accelerometers.
%% Load Accelerometer positions
acc_pos = readtable('mat/acc_pos.txt', 'ReadVariableNames', false);
acc_pos = table2array(acc_pos(:, 1:4));
[~, i] = sort(acc_pos(:, 1));
acc_pos = acc_pos(i, 2:4);
% Signal Processing :noexport:
% <<ssec:modal_signal_processing>>
% The measurements are averaged 10 times corresponding to 10 hammer impacts in order to reduce the effect of random noise.
% Windowing is also used on the force and response signals.
% A boxcar window (figure ref:fig:modal_windowing_force_signal) is used for the force signal as once the impact on the structure is done, the measured signal is meaningless.
% The parameters are:
% - *Start*: $3\%$
% - *Stop*: $7\%$
%% Boxcar window used for the force signal
figure;
plot(100*[0, 0.03, 0.03, 0.07, 0.07, 1], [0, 0, 1, 1, 0, 0]);
xlabel('Time [\%]'); ylabel('Amplitude');
xlim([0, 100]); ylim([0, 1]);
% #+name: fig:modal_windowing_force_signal
% #+caption: Boxcar window used for the force signal
% #+RESULTS:
% [[file:figs/modal_windowing_force_signal.png]]
% An exponential window (figure ref:fig:modal_windowing_acc_signal) is used for the response signal as we are measuring transient signals and most of the information is located at the beginning of the signal.
% The parameters are:
% - FlatTop:
% - *Start*: $3\%$
% - *Stop*: $2.96\%$
% - Decreasing point:
% - *X*: $60.4\%$
% - *Y*: $14.7\%$
%% Exponential window used for acceleration signal
x0 = 0.296;
xd = 0.604;
yd = 0.147;
alpha = log(yd)/(x0 - xd);
t = x0:0.01:1.01;
y = exp(-alpha*(t-x0));
figure;
plot(100*[0, 0.03, 0.03, x0, t], [0, 0, 1, 1, y]);
xlabel('Time [\%]'); ylabel('Amplitude');
xlim([0, 100]); ylim([0, 1]);
% Force and Response signals
% <<ssec:modal_measured_signals>>
%% Load raw data
meas1_raw = load('mat/meas_raw_1.mat');
% Sampling Frequency [Hz]
Fs = 1/meas1_raw.Track1_X_Resolution;
% Time just before the impact occurs [s]
impacts = [5.937, 11.228, 16.681, 22.205, 27.350, 32.714, 38.115, 43.888, 50.407]-0.01;
% Time vector [s]
time = linspace(0, meas1_raw.Track1_X_Resolution*length(meas1_raw.Track1), length(meas1_raw.Track1));
% The force sensor and the accelerometers signals are shown in the time domain in Figure ref:fig:modal_raw_meas.
% Sharp "impacts" can be seen for the force sensor, indicating wide frequency band excitation.
% For the accelerometer, many resonances can be seen on the right, indicating complex dynamics
%% Raw measurement of the Accelerometer
figure;
tiledlayout(1, 3, 'TileSpacing', 'Compact', 'Padding', 'None');
ax1 = nexttile([1,2]);
hold on;
plot(time, meas1_raw.Track2, 'DisplayName', 'Acceleration [$m/s^2$]');
plot(time, 1e-3*meas1_raw.Track1, 'DisplayName', 'Force [kN]');
hold off;
xlabel('Time [s]');
ylabel('Amplitude');
xlim([0, time(end)]);
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 1);
ax2 = nexttile();
hold on;
plot(time, meas1_raw.Track2);
plot(time, 1e-3*meas1_raw.Track1);
hold off;
xlabel('Time [s]');
set(gca, 'YTickLabel',[]);
xlim([22.19, 22.4]);
linkaxes([ax1,ax2],'y');
ylim([-2, 2]);
% #+name: fig:modal_raw_meas
% #+caption: Raw measurement of the acceleromter (blue) and of the force sensor at the Hammer tip (red). Zoom on one impact is shown on the right.
% #+RESULTS:
% [[file:figs/modal_raw_meas.png]]
%% Frequency Analysis
Nfft = floor(5.0*Fs); % Number of frequency points
win = hanning(Nfft); % Windowing
Noverlap = floor(Nfft/2); % Overlap for frequency analysis
%% Comnpute the power spectral density of the force and acceleration
[pxx_force, f] = pwelch(meas1_raw.Track1, win, Noverlap, Nfft, Fs);
[pxx_acc, ~] = pwelch(meas1_raw.Track2, win, Noverlap, Nfft, Fs);
% The "normalized" amplitude spectral density of the two signals are computed and shown in Figure ref:fig:modal_asd_acc_force.
% Conclusions based on the time domain signals can be clearly seen in the frequency domain (wide frequency content for the force signal and complex dynamics for the accelerometer).
%% Normalized Amplitude Spectral Density of the measured force and acceleration
figure;
hold on;
plot(f, sqrt(pxx_force./max(pxx_force(f<200))), 'DisplayName', 'Force');
plot(f, sqrt(pxx_acc./max(pxx_acc(f<200))), 'DisplayName', 'Acceleration');
hold off;
set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Normalized Spectral Density');
xlim([0, 200]);
xticks([0:20:200]);
ylim([0, 1])
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 1);
% #+name: fig:modal_asd_acc_force
% #+caption: Normalized Amplitude Spectral Density of the measured force and acceleration
% #+RESULTS:
% [[file:figs/modal_asd_acc_force.png]]
% The frequency response function from the applied force to the measured acceleration can then be computed (Figure ref:fig:modal_frf_acc_force).
%% Compute the transfer function and Coherence
[G1, f] = tfestimate(meas1_raw.Track1, meas1_raw.Track2, win, Noverlap, Nfft, Fs);
[coh1, ~] = mscohere( meas1_raw.Track1, meas1_raw.Track2, win, Noverlap, Nfft, Fs);
%% Frequency Response Function between the force and the acceleration
figure;
plot(f, abs(G1));
xlabel('Frequency [Hz]'); ylabel('FRF [$m/s^2/N$]')
set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'log');
xlim([0, 200]);
xticks([0:20:200]);
% #+name: fig:modal_frf_acc_force
% #+caption: Frequency Response Function between the measured force and acceleration
% #+RESULTS:
% [[file:figs/modal_frf_acc_force.png]]
% The coherence between the input and output signals is also computed and found to be good between 20 and 200Hz (Figure ref:fig:modal_coh_acc_force).
%% Frequency Response Function between the force and the acceleration
figure;
plot(f, coh1);
xlabel('Frequency [Hz]'); ylabel('Coherence [-]')
set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'lin');
xlim([0, 200]); ylim([0,1]);
xticks([0:20:200]);

378
matlab/modal_2_frf_processing.m Normal file
View File

@@ -0,0 +1,378 @@
%% Clear Workspace and Close figures
clear; close all; clc;
%% Intialize Laplace variable
s = zpk('s');
%% Path for functions, data and scripts
addpath('./mat/'); % Path for data
%% Colors for the figures
colors = colororder;
% Frequency Response Matrix
% <<ssec:modal_frf_matrix>>
% All the Frequency Response Functions measured are combined into one big array called the Frequency Response Matrix.
% The frequency response matrix is an $n \times p \times q$ matrix with:
% - $n$ the number of measurements: $23 \times 3$ (23 accelerometers measuring 3 directions each)
% - $p$ the number of excitation inputs: $3$
% - $q$ the number of frequency points $\omega_i$
% Thus, the FRF matrix is an $69 \times 3 \times 801$ matrix.
% For each frequency point $\omega_i$, a 2D matrix is obtained:
% \begin{equation}
% \text{FRF}(\omega_i) = \begin{bmatrix}
% \frac{D_{1_x}}{F_x}(\omega_i) & \frac{D_{1_x}}{F_y}(\omega_i) & \frac{D_{1_x}}{F_z}(\omega_i) \\
% \frac{D_{1_y}}{F_x}(\omega_i) & \frac{D_{1_y}}{F_y}(\omega_i) & \frac{D_{1_y}}{F_z}(\omega_i) \\
% \frac{D_{1_z}}{F_x}(\omega_i) & \frac{D_{1_z}}{F_y}(\omega_i) & \frac{D_{1_z}}{F_z}(\omega_i) \\
% \frac{D_{2_x}}{F_x}(\omega_i) & \frac{D_{2_x}}{F_y}(\omega_i) & \frac{D_{2_x}}{F_z}(\omega_i) \\
% \vdots & \vdots & \vdots \\
% \frac{D_{23_z}}{F_x}(\omega_i) & \frac{D_{23_z}}{F_y}(\omega_i) & \frac{D_{23_z}}{F_z}(\omega_i) \\
% \end{bmatrix}
% \end{equation}
%% Load frequency response matrix
load('frf_matrix.mat', 'freqs', 'frf');
%% Accelerometers ID connected to each solid body
solids = {};
solids.gbot = [17, 18, 19, 20]; % bottom granite
solids.gtop = [13, 14, 15, 16]; % top granite
solids.ty = [9, 10, 11, 12]; % Ty stage
solids.ry = [5, 6, 7, 8]; % Ry stage
solids.rz = [21, 22, 23]; % Rz stage
solids.hexa = [1, 2, 3, 4]; % Hexapod
% Names of the solid bodies
solid_names = fields(solids);
%% Save the acceleromter positions are well as the solid bodies
save('mat/geometry.mat', 'solids', 'solid_names', 'acc_pos');
% #+name: fig:aligned_accelerometers
% #+caption: Aligned measurement of the motion of a solid body
% #+RESULTS:
% [[file:figs/aligned_accelerometers.png]]
% The motion of the rigid body of figure ref:fig:aligned_accelerometers is defined by its displacement $\delta p$ and rotation $\vec{\Omega}$ with respect to the reference frame $\{O\}$.
% The motions at points $1$ and $2$ are:
% \begin{align*}
% \delta p_1 &= \delta p + \Omega \times p_1 \\
% \delta p_2 &= \delta p + \Omega \times p_2
% \end{align*}
% Taking only the $x$ direction:
% \begin{align*}
% \delta p_{x1} &= \delta p_x + \Omega_y p_{z1} - \Omega_z p_{y1} \\
% \delta p_{x2} &= \delta p_x + \Omega_y p_{z2} - \Omega_z p_{y2}
% \end{align*}
% However, we have $p_{1y} = p_{2y}$ and $p_{1z} = p_{2z}$ because of the co-linearity of the two sensors in the $x$ direction, and thus we obtain
% \begin{equation}
% \delta p_{x1} = \delta p_{x2}
% \end{equation}
% #+begin_important
% Two sensors that are measuring the motion of a rigid body in the direction of the line linking the two sensors should measure the same quantity.
% #+end_important
% We can verify that the rigid body assumption is correct by comparing the measurement of the sensors.
% From the table ref:tab:position_accelerometers, we can guess which sensors will give the same results in the X and Y directions.
% Comparison of such measurements in the X direction is shown on figure ref:fig:modal_solid_body_comp_x_dir.
% Similar result is obtained for the Y direction.
meas_dir = 1; % X
exc_dir = 1; % X
% Pair of accelerometers aligned in the X direction
acc_i = [1 , 4 ;
2 , 3 ;
5 , 8 ;
6 , 7 ;
9 , 12;
10, 11;
14, 15;
18, 19;
21, 23];
%% Comparaison of measured frequency response function for in the X directions for accelerometers aligned along X
figure;
tiledlayout(3, 3, 'TileSpacing', 'Compact', 'Padding', 'None');
for i = 1:size(acc_i, 1)
nexttile();
hold on;
plot(freqs, abs(squeeze(frf(meas_dir+3*(acc_i(i, 1)-1), exc_dir, :))), ...
'DisplayName', sprintf('%i', acc_i(i, 1)))
plot(freqs, abs(squeeze(frf(meas_dir+3*(acc_i(i, 2)-1), exc_dir, :))), ...
'DisplayName', sprintf('%i', acc_i(i, 2)))
hold off;
set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'log');
if i > 6
xlabel('Frequency [Hz]');
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 1);
else
set(gca, 'XTickLabel',[]);
legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 1);
end
if rem(i, 3) == 1
ylabel('Amplitude');
else
set(gca, 'YTickLabel',[]);
end
xlim([0, 200]); ylim([1e-5, 2e-2]);
end
% TODO Verification of the principle of reciprocity :noexport:
% <<ssec:modal_reciprocity_principle>>
% Because we expect our system to follow the principle of reciprocity.
% That is to say the response $X_j$ at some degree of freedom $j$ due to a force $F_k$ applied on DOF $k$ should be the same as the response $X_k$ due to a force $F_j$:
% \[ H_{jk} = \frac{X_j}{F_k} = \frac{X_k}{F_j} = H_{kj} \]
% This comes from the fact that we expect to have symmetric mass, stiffness and damping matrices.
% In order to access the quality of the data and the validity of the measured FRF, we then check that the reciprocity between $H_{jk}$ and $H_{kj}$ is of an acceptable level.
% We can verify this reciprocity using 3 different pairs of response/force.
dir_names = {'X', 'Y', 'Z'};
figure;
for i = 1:3
subplot(3, 1, i)
a = mod(i, 3)+1;
b = mod(i-2, 3)+1;
hold on;
plot(freqs, abs(squeeze(frf(3*(11-1)+a, b, :))), 'DisplayName', sprintf('$\\frac{F_%s}{D_%s}$', dir_names{a}, dir_names{b}));
plot(freqs, abs(squeeze(frf(3*(11-1)+b, a, :))), 'DisplayName', sprintf('$\\frac{F_%s}{D_%s}$', dir_names{b}, dir_names{a}));
hold off;
set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'log');
if i == 3
xlabel('Frequency [Hz]');
else
set(gca, 'XTickLabel',[]);
end
if i == 2
ylabel('Amplitude [$\frac{m/s^2}{N}$]');
end
xlim([0, 200]);
legend('location', 'northwest');
end
% From accelerometer DOFs to solid body DOFs - Matlab Implementation
% First, we initialize a new FRF matrix =FRFs_O= which is an $n \times p \times q$ with:
% - $n$ is the number of DOFs of the considered 6 solid-bodies: $6 \times 6 = 36$
% - $p$ is the number of excitation inputs: $3$
% - $q$ is the number of frequency points $\omega_i$
% #+begin_important
% For each frequency point $\omega_i$, the FRF matrix =FRFs_O= is a $n\times p$ matrix:
% \begin{equation}
% \text{FRF}_\text{O}(\omega_i) = \begin{bmatrix}
% \frac{D_{1,T_x}}{F_x}(\omega_i) & \frac{D_{1,T_x}}{F_y}(\omega_i) & \frac{D_{1,T_x}}{F_z}(\omega_i) \\
% \frac{D_{1,T_y}}{F_x}(\omega_i) & \frac{D_{1,T_y}}{F_y}(\omega_i) & \frac{D_{1,T_y}}{F_z}(\omega_i) \\
% \frac{D_{1,T_z}}{F_x}(\omega_i) & \frac{D_{1,T_z}}{F_y}(\omega_i) & \frac{D_{1,T_z}}{F_z}(\omega_i) \\
% \frac{D_{1,R_x}}{F_x}(\omega_i) & \frac{D_{1,R_x}}{F_y}(\omega_i) & \frac{D_{1,R_x}}{F_z}(\omega_i) \\
% \frac{D_{1,R_y}}{F_x}(\omega_i) & \frac{D_{1,R_y}}{F_y}(\omega_i) & \frac{D_{1,R_y}}{F_z}(\omega_i) \\
% \frac{D_{1,R_z}}{F_x}(\omega_i) & \frac{D_{1,R_z}}{F_y}(\omega_i) & \frac{D_{1,R_z}}{F_z}(\omega_i) \\
% \frac{D_{2,T_x}}{F_x}(\omega_i) & \frac{D_{2,T_x}}{F_y}(\omega_i) & \frac{D_{2,T_x}}{F_z}(\omega_i) \\
% \vdots & \vdots & \vdots \\
% \frac{D_{6,R_z}}{F_x}(\omega_i) & \frac{D_{6,R_z}}{F_y}(\omega_i) & \frac{D_{6,R_z}}{F_z}(\omega_i)
% \end{bmatrix}
% \end{equation}
% where $D_i$ corresponds to the solid body number i.
% #+end_important
% Then, as we know the positions of the accelerometers on each solid body, and we have the response of those accelerometers, we can use the equations derived in the previous section to determine the response of each solid body expressed in the frame $\{O\}$.
FRFs_O = zeros(length(solid_names)*6, 3, 801);
for solid_i = 1:length(solid_names)
solids_i = solids.(solid_names{solid_i});
A = zeros(3*length(solids_i), 6);
for i = 1:length(solids_i)
acc_i = solids_i(i);
A(3*(i-1)+1:3*i, 1:3) = eye(3);
A(3*(i-1)+1:3*i, 4:6) = [ 0 acc_pos(acc_i, 3) -acc_pos(acc_i, 2) ;
-acc_pos(acc_i, 3) 0 acc_pos(acc_i, 1) ;
acc_pos(acc_i, 2) -acc_pos(acc_i, 1) 0];
end
for exc_dir = 1:3
FRFs_O((solid_i-1)*6+1:solid_i*6, exc_dir, :) = A\squeeze(FRFs((solids_i(1)-1)*3+1:solids_i(end)*3, exc_dir, :));
end
end
% Center of Mass of each solid body
% From solidworks, we can export the position of the center of mass of each solid body considered.
% These are summarized in Table ref:tab:modal_com_solid_bodies
%% Extract the CoM of considered solid bodies
model_com = reshape(table2array(readtable('mat/model_solidworks_com.txt', 'ReadVariableNames', false)), [3, 6]);
% From accelerometer DOFs to solid body DOFs - Expressed at the CoM
% First, we initialize a new FRF matrix which is an $n \times p \times q$ with:
% - $n$ is the number of DOFs of the considered 6 solid-bodies: $6 \times 6 = 36$
% - $p$ is the number of excitation inputs: $3$
% - $q$ is the number of frequency points $\omega_i$
% #+begin_important
% For each frequency point $\omega_i$, the FRF matrix is a $n\times p$ matrix:
% \begin{equation}
% \text{FRF}_\text{CoM}(\omega_i) = \begin{bmatrix}
% \frac{D_{1,T_x}}{F_x}(\omega_i) & \frac{D_{1,T_x}}{F_y}(\omega_i) & \frac{D_{1,T_x}}{F_z}(\omega_i) \\
% \frac{D_{1,T_y}}{F_x}(\omega_i) & \frac{D_{1,T_y}}{F_y}(\omega_i) & \frac{D_{1,T_y}}{F_z}(\omega_i) \\
% \frac{D_{1,T_z}}{F_x}(\omega_i) & \frac{D_{1,T_z}}{F_y}(\omega_i) & \frac{D_{1,T_z}}{F_z}(\omega_i) \\
% \frac{D_{1,R_x}}{F_x}(\omega_i) & \frac{D_{1,R_x}}{F_y}(\omega_i) & \frac{D_{1,R_x}}{F_z}(\omega_i) \\
% \frac{D_{1,R_y}}{F_x}(\omega_i) & \frac{D_{1,R_y}}{F_y}(\omega_i) & \frac{D_{1,R_y}}{F_z}(\omega_i) \\
% \frac{D_{1,R_z}}{F_x}(\omega_i) & \frac{D_{1,R_z}}{F_y}(\omega_i) & \frac{D_{1,R_z}}{F_z}(\omega_i) \\
% \frac{D_{2,T_x}}{F_x}(\omega_i) & \frac{D_{2,T_x}}{F_y}(\omega_i) & \frac{D_{2,T_x}}{F_z}(\omega_i) \\
% \vdots & \vdots & \vdots \\
% \frac{D_{6,R_z}}{F_x}(\omega_i) & \frac{D_{6,R_z}}{F_y}(\omega_i) & \frac{D_{6,R_z}}{F_z}(\omega_i)
% \end{bmatrix}
% \end{equation}
% where 1, 2, ..., 6 corresponds to the 6 solid bodies.
% #+end_important
% Then, as we know the positions of the accelerometers on each solid body, and we have the response of those accelerometers, we can use the equations derived in the previous section to determine the response of each solid body expressed in their center of mass.
%% Frequency Response Matrix - Response expressed at the CoM of the solid bodies
FRFs_CoM = zeros(length(solid_names)*6, 3, 801);
for solid_i = 1:length(solid_names)
% Number of accelerometers fixed to this solid body
solids_i = solids.(solid_names{solid_i});
% "Jacobian" matrix to go from accelerometer frame to CoM frame
A = zeros(3*length(solids_i), 6);
for i = 1:length(solids_i)
acc_i = solids_i(i);
acc_pos_com = acc_pos(acc_i, :).' - model_com(:, solid_i);
A(3*(i-1)+1:3*i, 1:3) = eye(3);
A(3*(i-1)+1:3*i, 4:6) = [ 0 acc_pos_com(3) -acc_pos_com(2) ;
-acc_pos_com(3) 0 acc_pos_com(1) ;
acc_pos_com(2) -acc_pos_com(1) 0];
end
for exc_dir = 1:3
FRFs_CoM((solid_i-1)*6+1:solid_i*6, exc_dir, :) = A\squeeze(frf((solids_i(1)-1)*3+1:solids_i(end)*3, exc_dir, :));
end
end
%% Save the computed FRF at the CoM
save('mat/frf_com.mat', 'FRFs_CoM');
% Verify that we find the original FRF from the FRF in the global coordinates
% We have computed the Frequency Response Functions Matrix representing the response of the 6 solid bodies in their 6 DOFs with respect to their center of mass.
% From the response of one body in its 6 DOFs, we should be able to compute the FRF of each of its accelerometer fixed to it during the measurement, supposing that this stage is a solid body.
% We can then compare the result with the original measurements.
% This will help us to determine if:
% - the previous inversion used is correct
% - the solid body assumption is correct in the frequency band of interest
% From the translation $\delta p$ and rotation $\delta \Omega$ of a solid body and the positions $p_i$ of the accelerometers attached to it, we can compute the response that would have been measured by the accelerometers using the following formula:
% \begin{align*}
% \delta p_1 &= \delta p + \delta\Omega p_1\\
% \delta p_2 &= \delta p + \delta\Omega p_2\\
% \delta p_3 &= \delta p + \delta\Omega p_3\\
% \delta p_4 &= \delta p + \delta\Omega p_4
% \end{align*}
% Thus, we can obtain the FRF matrix =FRFs_A= that gives the responses of the accelerometers to the forces applied by the hammer.
% It is implemented in matlab as follow:
FRFs_A = zeros(size(frf));
% For each excitation direction
for exc_dir = 1:3
% For each solid
for solid_i = 1:length(solid_names)
v0 = squeeze(FRFs_CoM((solid_i-1)*6+1:(solid_i-1)*6+3, exc_dir, :));
W0 = squeeze(FRFs_CoM((solid_i-1)*6+4:(solid_i-1)*6+6, exc_dir, :));
% For each accelerometer attached to the current solid
for acc_i = solids.(solid_names{solid_i})
% We get the position of the accelerometer expressed in frame O
pos = acc_pos(acc_i, :).' - model_com(:, solid_i);
% pos = acc_pos(acc_i, :).';
posX = [0 pos(3) -pos(2); -pos(3) 0 pos(1) ; pos(2) -pos(1) 0];
FRFs_A(3*(acc_i-1)+1:3*(acc_i-1)+3, exc_dir, :) = v0 + posX*W0;
end
end
end
% We then compare the original FRF measured for each accelerometer =FRFs= with the "recovered" FRF =FRFs_A= from the global FRF matrix in the common frame.
% The FRF for the 4 accelerometers on the Hexapod are compared on figure ref:fig:recovered_frf_comparison_hexa.
% All the FRF are matching very well in all the frequency range displayed.
% The FRF for accelerometers located on the translation stage are compared on figure ref:fig:recovered_frf_comparison_ty.
% The FRF are matching well until 100Hz.
%% Comparaison of the original accelerometer response and reconstructed response from the solid body response
exc_names = {'$F_x$', '$F_y$', '$F_z$'};
DOFs = {'x', 'y', 'z', '\theta_x', '\theta_y', '\theta_z'};
solid_i = 6; % Considered solid body
exc_dir = 1; % Excited direction
accs_i = solids.(solid_names{solid_i}); % Accelerometers fixed to this solid body
figure;
tiledlayout(2, 2, 'TileSpacing', 'Compact', 'Padding', 'None');
for i = 1:length(accs_i)
acc_i = accs_i(i);
nexttile();
hold on;
for dir_i = 1:3
plot(freqs, abs(squeeze(frf(3*(acc_i-1)+dir_i, exc_dir, :))), '-', 'DisplayName', sprintf('$a_{%i,%s}$', acc_i, DOFs{dir_i}));
end
set(gca,'ColorOrderIndex',1)
for dir_i = 1:3
plot(freqs, abs(squeeze(FRFs_A(3*(acc_i-1)+dir_i, exc_dir, :))), '--', 'HandleVisibility', 'off');
end
hold off;
if i > 2
xlabel('Frequency [Hz]');
else
set(gca, 'XTickLabel',[]);
end
if rem(i, 2) == 1
ylabel('Amplitude [$\frac{m/s^2}{N}$]');
else
set(gca, 'YTickLabel',[]);
end
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlim([5, 200]); ylim([1e-6, 1e-1]);
legend('location', 'northwest', 'FontSize', 8, 'NumColumns', 1);
end

271
matlab/modal_3_analysis.m Normal file
View File

@@ -0,0 +1,271 @@
%% Clear Workspace and Close figures
clear; close all; clc;
%% Intialize Laplace variable
s = zpk('s');
%% Path for functions, data and scripts
addpath('./mat/'); % Path for data
%% Colors for the figures
colors = colororder;
% Singular Value Decomposition - Modal Indication Function
% The Mode Indicator Functions are usually used on $n\times p$ FRF matrix where $n$ is a relatively large number of measurement DOFs and $p$ is the number of excitation DOFs, typically 3 or 4.
% In these methods, the frequency dependent FRF matrix is subjected to a singular value decomposition analysis which thus yields a small number (3 or 4) of singular values, these also being frequency dependent.
% These methods are used to *determine the number of modes* present in a given frequency range, to *identify repeated natural frequencies* and to pre-process the FRF data prior to modal analysis.
% From the documentation of the modal software:
% #+begin_quote
% The MIF consist of the singular values of the Frequency response function matrix. The number of MIFs equals the number of excitations.
% By the powerful singular value decomposition, the real signal space is separated from the noise space. Therefore, the MIFs exhibit the modes effectively.
% A peak in the MIFs plot usually indicate the existence of a structural mode, and two peaks at the same frequency point means the existence of two repeated modes.
% Moreover, the magnitude of the MIFs implies the strength of the a mode.
% #+end_quote
% #+begin_important
% The *Complex Mode Indicator Function* is defined simply by the SVD of the FRF (sub) matrix:
% \begin{align*}
% [H(\omega)]_{n\times p} &= [U(\omega)]_{n\times n} [\Sigma(\omega)]_{n\times p} [V(\omega)]_{p\times p}^H\\
% [CMIF(\omega)]_{p\times p} &= [\Sigma(\omega)]_{p\times n}^T [\Sigma(\omega)]_{n\times p}
% \end{align*}
% #+end_important
% We compute the Complex Mode Indicator Function.
% The result is shown on Figure ref:fig:modal_indication_function.
%% Computation of the modal indication function
MIF = zeros(size(frf, 2), size(frf, 2), size(frf, 3));
for i = 1:length(freqs)
[~,S,~] = svd(frf(:, :, i));
MIF(:, :, i) = S'*S;
end
%% Modal Indication Function
figure;
hold on;
for i = 1:size(MIF, 1)
plot(freqs, squeeze(MIF(i, i, :)));
end
hold off;
set(gca, 'Xscale', 'lin'); set(gca, 'Yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('CMIF Amplitude');
xticks([0:20:200]);
xlim([0, 200]);
ylim([1e-6, 2e-2]);
% Composite Response Function
% An alternative is the Composite Response Function $HH(\omega)$ defined as the sum of all the measured FRF:
% \begin{equation}
% HH(\omega) = \sum_j\sum_kH_{jk}(\omega)
% \end{equation}
% Instead, we choose here to use the sum of the norms of the measured frf:
% \begin{equation}
% HH(\omega) = \sum_j\sum_k \left|H_{jk}(\omega) \right|
% \end{equation}
% The result is shown on figure ref:fig:modal_composite_reponse_function.
%% Composite Response Function
figure;
hold on;
plot(freqs, squeeze(sum(sum(abs(frf)))), '-k');
hold off;
xlabel('Frequency [Hz]'); ylabel('Amplitude');
xlim([0, 200]);
xticks([0:20:200]);
% Importation of the modal parameters on Matlab
% The obtained modal parameters are:
% - Resonance frequencies in Hertz
% - Modal damping ratio in percentage
% - (complex) Modes shapes for each measured DoF
% - Modal A and modal B which are parameters important for further normalization
%% Load modal parameters
shapes_m = readtable('mat/mode_shapes.txt', 'ReadVariableNames', false); % [Sign / Real / Imag]
freqs_m = table2array(readtable('mat/mode_freqs.txt', 'ReadVariableNames', false)); % in [Hz]
damps_m = table2array(readtable('mat/mode_damps.txt', 'ReadVariableNames', false)); % in [%]
modal_a = table2array(readtable('mat/mode_modal_a.txt', 'ReadVariableNames', false)); % [Real / Imag]
modal_b = table2array(readtable('mat/mode_modal_b.txt', 'ReadVariableNames', false)); % [Real / Imag]
%% Guess the number of modes identified from the length of the imported data.
acc_n = 23; % Number of accelerometers
dir_n = 3; % Number of directions
dirs = 'XYZ';
mod_n = size(shapes_m,1)/acc_n/dir_n; % Number of modes
%% Mode shapes are split into 3 parts (direction plus sign, real part and imaginary part)
% we aggregate them into one array of complex numbers
T_sign = table2array(shapes_m(:, 1));
T_real = table2array(shapes_m(:, 2));
T_imag = table2array(shapes_m(:, 3));
mode_shapes = zeros(mod_n, dir_n, acc_n);
for mod_i = 1:mod_n
for acc_i = 1:acc_n
% Get the correct section of the signs
T = T_sign(acc_n*dir_n*(mod_i-1)+1:acc_n*dir_n*mod_i);
for dir_i = 1:dir_n
% Get the line corresponding to the sensor
i = find(contains(T, sprintf('%i%s',acc_i, dirs(dir_i))), 1, 'first')+acc_n*dir_n*(mod_i-1);
mode_shapes(mod_i, dir_i, acc_i) = str2num([T_sign{i}(end-1), '1'])*complex(T_real(i),T_imag(i));
end
end
end
% Modal Matrices
% We would like to arrange the obtained modal parameters into two modal matrices:
% \[ \Lambda = \begin{bmatrix}
% s_1 & & 0 \\
% & \ddots & \\
% 0 & & s_N
% \end{bmatrix}_{N \times N}; \quad \Psi = \begin{bmatrix}
% & & \\
% \{\psi_1\} & \dots & \{\psi_N\} \\
% & &
% \end{bmatrix}_{M \times N} \]
% \[ \{\psi_i\} = \begin{Bmatrix} \psi_{i, 1_x} & \psi_{i, 1_y} & \psi_{i, 1_z} & \psi_{i, 2_x} & \dots & \psi_{i, 23_z} \end{Bmatrix}^T \]
% $M$ is the number of DoF: here it is $23 \times 3 = 69$.
% $N$ is the number of mode
eigen_val_M = diag(2*pi*freqs_m.*(-damps_m/100 + j*sqrt(1 - (damps_m/100).^2)));
eigen_vec_M = reshape(mode_shapes, [mod_n, acc_n*dir_n]).';
% Each eigen vector is normalized: $\| \{\psi_i\} \|_2 = 1$
% However, the eigen values and eigen vectors appears as complex conjugates:
% \[ s_r, s_r^*, \{\psi\}_r, \{\psi\}_r^*, \quad r = 1, N \]
% In the end, they are $2N$ eigen values.
% We then build two extended eigen matrices as follow:
% \[ \mathcal{S} = \begin{bmatrix}
% s_1 & & & & & \\
% & \ddots & & & 0 & \\
% & & s_N & & & \\
% & & & s_1^* & & \\
% & 0 & & & \ddots & \\
% & & & & & s_N^*
% \end{bmatrix}_{2N \times 2N}; \quad \Phi = \begin{bmatrix}
% & & & & &\\
% \{\psi_1\} & \dots & \{\psi_N\} & \{\psi_1^*\} & \dots & \{\psi_N^*\} \\
% & & & & &
% \end{bmatrix}_{M \times 2N} \]
eigen_val_ext_M = blkdiag(eigen_val_M, conj(eigen_val_M));
eigen_vec_ext_M = [eigen_vec_M, conj(eigen_vec_M)];
% We also build the Modal A and Modal B matrices:
% \begin{equation}
% A = \begin{bmatrix}
% a_1 & & 0 \\
% & \ddots & \\
% 0 & & a_N
% \end{bmatrix}_{N \times N}; \quad B = \begin{bmatrix}
% b_1 & & 0 \\
% & \ddots & \\
% 0 & & b_N
% \end{bmatrix}_{N \times N}
% \end{equation}
% With $a_i$ is the "Modal A" parameter linked to mode i.
modal_a_M = diag(complex(modal_a(:, 1), modal_a(:, 2)));
modal_b_M = diag(complex(modal_b(:, 1), modal_b(:, 2)));
modal_a_ext_M = blkdiag(modal_a_M, conj(modal_a_M));
modal_b_ext_M = blkdiag(modal_b_M, conj(modal_b_M));
% Matlab Implementation
Hsyn = zeros(acc_n*dir_n, acc_n*dir_n, length(freqs));
for i = 1:length(freqs)
Hsyn(:, :, i) = eigen_vec_ext_M*((j*2*pi*freqs(i)).^2*inv(modal_a_ext_M)/(diag(j*2*pi*freqs(i) - diag(eigen_val_ext_M))))*eigen_vec_ext_M.';
end
% Because the synthesize frequency response functions are representing the displacement response in $[m/N]$, we multiply each element of the FRF matrix by $(j \omega)^2$ in order to obtain the acceleration response in $[m/s^2/N]$.
for i = 1:size(Hsyn, 1)
Hsyn(i, :, :) = squeeze(Hsyn(i, :, :)).*(j*2*pi*freqs).^2;
end
% Original and Synthesize FRF matrix comparison
acc_o = 1; dir_o = 1; dir_i = 1;
figure;
ax1 = subplot(2, 1, 1);
hold on;
plot(freqs, abs(squeeze(frf(3*(acc_o-1)+dir_o, dir_i, :))), 'DisplayName', 'Original');
plot(freqs, abs(squeeze(Hsyn(3*(acc_o-1)+dir_o, 3*(11-1)+dir_i, :))), 'DisplayName', 'Synthesize');
hold off;
set(gca, 'yscale', 'log');
set(gca, 'XTickLabel',[]);
ylabel('Magnitude [$\frac{m/s^2}{N}$]');
title(sprintf('From acc %i %s to acc %i %s', 11, dirs(dir_i), acc_o, dirs(dir_o)))
legend('location', 'northwest');
ax2 = subplot(2, 1, 2);
hold on;
plot(freqs, mod(180/pi*phase(squeeze(frf(3*(acc_o-1)+dir_o, dir_i, :)))+180, 360)-180);
plot(freqs, mod(180/pi*phase(squeeze(Hsyn(3*(acc_o-1)+dir_o, 3*(11-1)+dir_i, :)))+180, 360)-180);
hold off;
yticks(-360:90:360); ylim([-180, 180]);
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
linkaxes([ax1,ax2],'x');
xlim([1, 200]);
% Synthesize FRF that has not yet been measured
accs = [1]; dirs = [1:3];
figure;
ax1 = subplot(2, 1, 1);
hold on;
for acc_i = accs
for dir_i = dirs
plot(freqs, abs((1./(j*2*pi*freqs').^2).*squeeze(Hsyn(3*(acc_i-1)+dir_i, 3*(acc_i-1)+dir_i, :))), 'DisplayName', sprintf('Acc %i - %s', acc_i, dirs(dir_i)));
end
end
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
set(gca, 'XTickLabel',[]);
ylabel('Magnitude [$\frac{m}{N}$]');
legend('location', 'southwest');
ax2 = subplot(2, 1, 2);
hold on;
for acc_i = accs
for dir_i = dirs
plot(freqs, mod(180/pi*phase((1./(j*2*pi*freqs').^2).*squeeze(Hsyn(3*(acc_i-1)+dir_i, 3*(acc_i-1)+dir_i, :)))+180, 360)-180);
end
end
hold off;
yticks(-360:90:360); ylim([-180, 180]);
set(gca, 'xscale', 'log');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
linkaxes([ax1,ax2],'x');
xlim([1, 200]);