Modal Analysis - Modal Parameter Extraction

1. Setup
2. Mode extraction and importation
3. Positions of the sensors
4. Solids
5. From local coordinates to global coordinates
6. Modal Matrices
7. Modal Complexity
8. Some notes about constraining the number of degrees of freedom
9. TODO Normalization of mode shapes?
10. Compare Mode Shapes
11. Importation of measured FRF curves
12. Importation of measured FRF curves to global FRF matrix
13. Analysis of some FRFs
14. From local coordinates to global coordinates with the FRFs
15. Analysis of some FRF in the global coordinates
16. Compare global coordinates to local coordinates
17. Verify that we find the original FRF from the FRF in the global coordinates
18. TODO Synthesis of FRF curves

1 Setup

Figure 1: Position and orientation of the accelerometer used

2 Mode extraction and importation

First, we split the big modes.asc files into sub text files using bash.

sed '/^\s*[0-9]*[XYZ][+-]:/!d' modal_analysis_updated/modes.asc > mat/mode_shapes.txt
sed '/freq/!d' modal_analysis_updated/modes.asc | sed 's/.* = \(.*\)Hz/\1/' > mat/mode_freqs.txt
sed '/damp/!d' modal_analysis_updated/modes.asc | sed 's/.* = \(.*\)\%/\1/' > mat/mode_damps.txt
sed '/modal A/!d' modal_analysis_updated/modes.asc | sed 's/.* =\s\+\([-0-9.e]\++[0-9]\+\)\([-+0-9.e]\+\)i/\1 \2/' > mat/mode_modal_a.txt
sed '/modal B/!d' modal_analysis_updated/modes.asc | sed 's/.* =\s\+\([-0-9.e]\++[0-9]\+\)\([-+0-9.e]\+\)i/\1 \2/' > mat/mode_modal_b.txt

Then we import them on Matlab.

shapes = readtable('mat/mode_shapes.txt', 'ReadVariableNames', false); % [Sign / Real / Imag]
freqs = table2array(readtable('mat/mode_freqs.txt', 'ReadVariableNames', false)); % in [Hz]
damps = table2array(readtable('mat/mode_damps.txt', 'ReadVariableNames', false)); % in [%]
modal_a = table2array(readtable('mat/mode_modal_a.txt', 'ReadVariableNames', false)); % [Real / Imag]
modal_a = complex(modal_a(:, 1), modal_a(:, 2));
modal_b = table2array(readtable('mat/mode_modal_b.txt', 'ReadVariableNames', false)); % [Real / Imag]
modal_b = complex(modal_b(:, 1), modal_b(:, 2));

We 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,1)/acc_n/dir_n; % Number of modes

As the 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(:, 1));
T_real = table2array(shapes(:, 2));
T_imag = table2array(shapes(:, 3));

modes = zeros(mod_n, acc_n, dir_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);
      modes(mod_i, acc_i, dir_i) = str2num([T_sign{i}(end-1), '1'])*complex(T_real(i),T_imag(i));
    end
  end
end

The obtained mode frequencies and damping are shown below.

data2orgtable([freqs, damps], {}, {'Frequency [Hz]', 'Damping [%]'}, ' %.1f ');

Frequency [Hz]	Damping [%]
11.4	8.7
18.5	11.8
37.6	6.4
39.4	3.6
54.0	0.2
56.1	2.8
69.7	4.6
71.6	0.6
72.4	1.6
84.9	3.6
90.6	0.3
91.0	2.9
95.8	3.3
105.4	3.3
106.8	1.9
112.6	3.0
116.8	2.7
124.1	0.6
145.4	1.6
150.1	2.2
164.7	1.4

3 Positions of the sensors

We process the file exported from the modal software containing the positions of the sensors using bash.

cat modal_analysis_updated/id31_nanostation_modified.cfg | grep NODES -A 23 | sed '/\s\+[0-9]\+/!d' | sed 's/\(.*\)\s\+0\s\+.\+/\1/' > mat/acc_pos.txt

We then import that on matlab, and sort them.

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);

The positions of the sensors relative to the point of interest are shown below.

data2orgtable(1000*acc_pos, {}, {'x [mm]', 'y [mm]', 'z [mm]'}, ' %.0f ');

x [mm]	y [mm]	z [mm]
-64	-64	-296
-64	64	-296
64	64	-296
64	-64	-296
-385	-300	-417
-420	280	-417
420	280	-417
380	-300	-417
-475	-414	-427
-465	407	-427
475	424	-427
475	-419	-427
-320	-446	-786
-480	534	-786
450	534	-786
295	-481	-786
-730	-526	-951
-735	814	-951
875	799	-951
865	-506	-951
-155	-90	-594
0	180	-594
155	-90	-594

4 Solids

We consider the following solid bodies:

Bottom Granite
Top Granite
Translation Stage
Tilt Stage
Spindle
Hexapod

We create a structure solids that contains the accelerometer number of each solid bodies (as shown on figure 1).

solids = {};
solids.granite_bot = [17, 18, 19, 20];
solids.granite_top = [13, 14, 15, 16];
solids.ty          = [9, 10, 11, 12];
solids.ry          = [5, 6, 7, 8];
solids.rz          = [21, 22, 23];
solids.hexa        = [1, 2, 3, 4];

solid_names = fields(solids);

5 From local coordinates to global coordinates

From the figure above, we can write:

\begin{align*} \vec{v}_1 &= \vec{v} + \Omega \vec{p}_1\\ \vec{v}_2 &= \vec{v} + \Omega \vec{p}_2\\ \vec{v}_3 &= \vec{v} + \Omega \vec{p}_3\\ \vec{v}_4 &= \vec{v} + \Omega \vec{p}_4 \end{align*}

With

\begin{equation} \Omega = \begin{bmatrix} 0 & -\Omega_z & \Omega_y \\ \Omega_z & 0 & -\Omega_x \\ -\Omega_y & \Omega_x & 0 \end{bmatrix} \end{equation}

\(\vec{v}\) and \(\Omega\) represent to velocity and rotation of the solid expressed in the frame \(\{O\}\).

We can rearrange the equations in a matrix form:

\begin{equation} \left[\begin{array}{ccc|ccc} 1 & 0 & 0 & 0 & p_{1z} & -p_{1y} \\ 0 & 1 & 0 & -p_{1z} & 0 & p_{1x} \\ 0 & 0 & 1 & p_{1y} & -p_{1x} & 0 \\ \hline & \vdots & & & \vdots & \\ \hline 1 & 0 & 0 & 0 & p_{4z} & -p_{4y} \\ 0 & 1 & 0 & -p_{4z} & 0 & p_{4x} \\ 0 & 0 & 1 & p_{4y} & -p_{4x} & 0 \end{array}\right] \begin{bmatrix} v_x \\ v_y \\ v_z \\ \hline \Omega_x \\ \Omega_y \\ \Omega_z \end{bmatrix} = \begin{bmatrix} v_{1x} \\ v_{1y} \\ v_{1z} \\\hline \vdots \\\hline v_{4x} \\ v_{4y} \\ v_{4z} \end{bmatrix} \end{equation}

and then we obtain the velocity and rotation of the solid in the wanted frame \(\{O\}\):

\begin{equation} \begin{bmatrix} v_x \\ v_y \\ v_z \\ \hline \Omega_x \\ \Omega_y \\ \Omega_z \end{bmatrix} = \left[\begin{array}{ccc|ccc} 1 & 0 & 0 & 0 & p_{1z} & -p_{1y} \\ 0 & 1 & 0 & -p_{1z} & 0 & p_{1x} \\ 0 & 0 & 1 & p_{1y} & -p_{1x} & 0 \\ \hline & \vdots & & & \vdots & \\ \hline 1 & 0 & 0 & 0 & p_{4z} & -p_{4y} \\ 0 & 1 & 0 & -p_{4z} & 0 & p_{4x} \\ 0 & 0 & 1 & p_{4y} & -p_{4x} & 0 \end{array}\right]^{-1} \begin{bmatrix} v_{1x} \\ v_{1y} \\ v_{1z} \\\hline \vdots \\\hline v_{4x} \\ v_{4y} \\ v_{4z} \end{bmatrix} \end{equation}

This inversion is equivalent to a mean square problem.

mode_shapes_O = zeros(mod_n, length(solid_names), 6);

for mod_i = 1:mod_n
  for solid_i = 1:length(solid_names)
    solids_i = solids.(solid_names{solid_i});

    Y = reshape(squeeze(modes(mod_i, solids_i, :))', [], 1);

    A = zeros(3*length(solids_i), 6);
    for i = 1:length(solids_i)
      A(3*(i-1)+1:3*i, 1:3) = eye(3);

      A(3*(i-1)+1:3*i, 4:6) = [0 acc_pos(i, 3) -acc_pos(i, 2) ; -acc_pos(i, 3) 0 acc_pos(i, 1) ; acc_pos(i, 2) -acc_pos(i, 1) 0];
    end

    mode_shapes_O(mod_i, solid_i, :) = A\Y;
  end
end

6 Modal Matrices

We want to obtain the two following matrices: \[ \Omega = \begin{bmatrix} \omega_1^2 & & 0 \\ & \ddots & \\ 0 & & \omega_n^2 \end{bmatrix}; \quad \Psi = \begin{bmatrix} & & \\ \{\psi_1\} & \dots & \{\psi_n\} \\ & & \end{bmatrix} \]

[ ] How to add damping to the eigen value matrix?

eigen_value_M = diag(freqs*2*pi);
eigen_vector_M = reshape(mode_shapes_O, [mod_n, 6*length(solid_names)])';

\[ \{\psi_1\} = \begin{Bmatrix} \psi_{1_x} & \psi_{2_x} & \dots & \psi_{6_x} & \psi_{1_x} & \dots & \psi_{1\Omega_x} & \dots & \psi_{6\Omega_z} \end{Bmatrix}^T \]

7 Modal Complexity

A method of displaying modal complexity is by plotting the elements of the eigenvector on an Argand diagram, such as the ones shown in figure 3.

To evaluate the complexity of the modes, we plot a polygon around the extremities of the individual vectors. The obtained area of this polygon is then compared with the area of the circle which is based on the length of the largest vector element. The resulting ratio is used as an indication of the complexity of the mode.

A little complex mode is shown on figure 3 whereas an highly complex mode is shown on figure 4. The complexity of all the modes are compared on figure 5.

Figure 3: Modal Complexity of one mode with small complexity

Figure 4: Modal Complexity of one higly complex mode

Figure 5: Modal complexity for each mode

8 Some notes about constraining the number of degrees of freedom

We want to have the two eigen matrices.

They should have the same size \(n \times n\) where \(n\) is the number of modes as well as the number of degrees of freedom. Thus, if we consider 21 modes, we should restrict our system to have only 21 DOFs.

Actually, we are measured 6 DOFs of 6 solids, thus we have 36 DOFs.

From the mode shapes animations, it seems that in the frequency range of interest, the two marbles can be considered as one solid. We thus have 5 solids and 30 DOFs.

In order to determine which DOF can be neglected, two solutions seems possible:

compare the mode shapes
compare the FRFs

The question is: in which base (frame) should be express the modes shapes and FRFs? Is it meaningful to compare mode shapes as they give no information about the amplitudes of vibration?

Stage	Motion DOFs	Parasitic DOF	Total DOF	Description of DOF
Granite	0	3	3
Ty	1	2	3	Ty, Rz
Ry	1	2	3	Ry,
Rz	1	2	3	Rz, Rx, Ry
Hexapod	6	0	6	Txyz, Rxyz
	9	9	18

9 TODO Normalization of mode shapes?

We normalize each column of the eigen vector matrix. Then, each eigenvector as a norm of 1.

eigen_vector_M = eigen_vector_M./vecnorm(eigen_vector_M);

10 Compare Mode Shapes

Let's say we want to see for the first mode which DOFs can be neglected. In order to do so, we should estimate the motion of each stage in particular directions. If we look at the z motion for instance, we will find that we cannot neglect that motion (because of the tilt causing z motion).

mode_i = 3;
dof_i = 6;

mode = eigen_vector_M(dof_i:6:end, mode_i);

figure;
hold on;
for i=1:length(mode)
  plot([0, real(mode(i))], [0, imag(mode(i))], '-', 'DisplayName', solid_names{i});
end
hold off;
legend();

figure;
subplot(2, 1, 1);
hold on;
for i=1:length(mode)
  plot(1, norm(mode(i)), 'o');
end
hold off;
ylabel('Amplitude');

subplot(2, 1, 2);
hold on;
for i=1:length(mode)
  plot(1, 180/pi*angle(mode(i)), 'o', 'DisplayName', solid_names{i});
end
hold off;
ylim([-180, 180]); yticks([-180:90:180]);
ylabel('Phase [deg]');
legend();

test = mode_shapes_O(10, 1, :)/norm(squeeze(mode_shapes_O(10, 1, :)));
test = mode_shapes_O(10, 2, :)/norm(squeeze(mode_shapes_O(10, 2, :)));

11 Importation of measured FRF curves

There are 24 measurements files corresponding to 24 series of impacts:

3 directions, 8 sets of 3 accelerometers

For each measurement file, the FRF and coherence between the impact and the 9 accelerations measured.

In reality: 4 sets of 10 things

a = load('./modal_analysis/frf_coh/Measurement1.mat');

figure;

ax1 = subplot(2, 1, 1);
hold on;
plot(a.FFT1_AvXSpc_2_1_RMS_X_Val, a.FFT1_AvXSpc_2_1_RMS_Y_Mod)
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
set(gca, 'XTickLabel',[]);
ylabel('Amplitude');
title(sprintf('From %s, to %s', FFT1_AvXSpc_2_1_RfName, FFT1_AvXSpc_2_1_RpName))

ax2 = subplot(2, 1, 2);
hold on;
plot(a.FFT1_AvXSpc_2_1_RMS_X_Val, a.FFT1_AvXSpc_2_1_RMS_Y_Phas)
hold off;
ylim([-180, 180]); yticks(-180:90:180);
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'xscale', 'log');

linkaxes([ax1,ax2],'x');
xlim([1, 200]);

12 Importation of measured FRF curves to global FRF matrix

FRF matrix \(n \times p\):

\(n\) is the number of measurements: \(3 \times 24\)
\(p\) is the number of excitation inputs: 3

23 measurements: 3 accelerometers

\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}

n_meas = 24;
n_acc = 23;

dirs = 'XYZ';

% Number of Accelerometer * DOF for each acccelerometer / Number of excitation / frequency points
FRFs = zeros(3*n_acc, 3, 801);
COHs = zeros(3*n_acc, 3, 801);

% Loop through measurements
for i = 1:n_meas
  % Load the measurement file
  meas = load(sprintf('./modal_analysis/frf_coh/Measurement%i.mat', i));

  % First: determine what is the exitation (direction and sign)
  exc_dir = meas.FFT1_AvXSpc_2_1_RMS_RfName(end);
  exc_sign = meas.FFT1_AvXSpc_2_1_RMS_RfName(end-1);
  % Determine what is the correct excitation sign
  exc_factor = str2num([exc_sign, '1']);
  if exc_dir ~= 'Z'
    exc_factor = exc_factor*(-1);
  end

  % Then: loop through the nine measurements and store them at the correct location
  for j = 2:10
    % Determine what is the accelerometer and direction
    [indices_acc_i] = strfind(meas.(sprintf('FFT1_H1_%i_1_RpName', j)), '.');
    acc_i = str2num(meas.(sprintf('FFT1_H1_%i_1_RpName', j))(indices_acc_i(1)+1:indices_acc_i(2)-1));

    meas_dir  = meas.(sprintf('FFT1_H1_%i_1_RpName', j))(end);
    meas_sign = meas.(sprintf('FFT1_H1_%i_1_RpName', j))(end-1);
    % Determine what is the correct measurement sign
    meas_factor = str2num([meas_sign, '1']);
    if meas_dir ~= 'Z'
      meas_factor = meas_factor*(-1);
    end

    % FRFs(acc_i+n_acc*(find(dirs==meas_dir)-1), find(dirs==exc_dir), :) = exc_factor*meas_factor*meas.(sprintf('FFT1_H1_%i_1_Y_ReIm', j));
    % COHs(acc_i+n_acc*(find(dirs==meas_dir)-1), find(dirs==exc_dir), :) = meas.(sprintf('FFT1_Coh_%i_1_RMS_Y_Val', j));

    FRFs(find(dirs==meas_dir)+3*(acc_i-1), find(dirs==exc_dir), :) = exc_factor*meas_factor*meas.(sprintf('FFT1_H1_%i_1_Y_ReIm', j));
    COHs(find(dirs==meas_dir)+3*(acc_i-1), find(dirs==exc_dir), :) = meas.(sprintf('FFT1_Coh_%i_1_RMS_Y_Val', j));
  end
end
freqs = meas.FFT1_Coh_10_1_RMS_X_Val;

13 Analysis of some FRFs

acc_i = 3;
acc_dir = 1;
exc_dir = 1;

figure;

ax1 = subplot(2, 1, 1);
hold on;
plot(freqs, abs(squeeze(FRFs(acc_dir+3*(acc_i-1), exc_dir, :))));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
set(gca, 'XTickLabel',[]);
ylabel('Amplitude');

ax2 = subplot(2, 1, 2);
hold on;
plot(freqs, mod(180+180/pi*phase(squeeze(FRFs(acc_dir+3*(acc_i-1), exc_dir, :))), 360)-180);
hold off;
ylim([-180, 180]); yticks(-180:90:180);
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'xscale', 'log');

linkaxes([ax1,ax2],'x');
xlim([1, 200]);

figure;
hold on;
for i = 1:3*n_acc
  plot(freqs, squeeze(COHs(i, 1, :)), 'color', [0, 0, 0, 0.2]);
end
hold off;
xlabel('Frequency [Hz]');
ylabel('Coherence [\%]');

Composite Response Function.

We here sum the norm instead of the complex numbers.

HHx = squeeze(sum(abs(FRFs(:, 1, :))));
HHy = squeeze(sum(abs(FRFs(:, 2, :))));
HHz = squeeze(sum(abs(FRFs(:, 3, :))));
HH  = squeeze(sum([HHx, HHy, HHz], 2));

exc_dir = 3;

figure;
hold on;
for i = 1:3*n_acc
  plot(freqs, abs(squeeze(FRFs(i, exc_dir, :))), 'color', [0, 0, 0, 0.2]);
end
plot(freqs, abs(HHx));
plot(freqs, abs(HHy));
plot(freqs, abs(HHz));
plot(freqs, abs(HH), 'k');
hold off;
set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Amplitude');
xlim([1, 200]);

14 From local coordinates to global coordinates with the FRFs

% Number of Solids * DOF for each solid / Number of excitation / frequency points
FRFs_O = zeros(length(solid_names)*6, 3, 801);

for exc_dir = 1:3
  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)
      A(3*(i-1)+1:3*i, 1:3) = eye(3);

      A(3*(i-1)+1:3*i, 4:6) = [0 acc_pos(i, 3) -acc_pos(i, 2) ; -acc_pos(i, 3) 0 acc_pos(i, 1) ; acc_pos(i, 2) -acc_pos(i, 1) 0];
    end

    for i = 1:801
      FRFs_O((solid_i-1)*6+1:solid_i*6, exc_dir, i) = A\FRFs((solids_i(1)-1)*3+1:solids_i(end)*3, exc_dir, i);
    end
  end
end

15 Analysis of some FRF in the global coordinates

solid_i = 6;
dir_i = 1;
exc_dir = 1;

figure;

ax1 = subplot(2, 1, 1);
hold on;
plot(freqs, abs(squeeze(FRFs_O((solid_i-1)*6+dir_i, exc_dir, :))));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
set(gca, 'XTickLabel',[]);
ylabel('Amplitude');

ax2 = subplot(2, 1, 2);
hold on;
plot(freqs, mod(180+180/pi*phase(squeeze(FRFs_O((solid_i-1)*6+dir_i, exc_dir, :))), 360)-180);
hold off;
ylim([-180, 180]); yticks(-180:90:180);
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'xscale', 'log');

linkaxes([ax1,ax2],'x');
xlim([1, 200]);

16 Compare global coordinates to local coordinates

solid_i = 1;
acc_dir_O = 6;
acc_dir = 3;
exc_dir = 3;

figure;

ax1 = subplot(2, 1, 1);
hold on;
for i = solids.(solid_names{solid_i})
  plot(freqs, abs(squeeze(FRFs(acc_dir+3*(i-1), exc_dir, :))));
end
plot(freqs, abs(squeeze(FRFs_O((solid_i-1)*6+acc_dir_O, exc_dir, :))), '-k');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
set(gca, 'XTickLabel',[]);
ylabel('Amplitude');

ax2 = subplot(2, 1, 2);
hold on;
for i = solids.(solid_names{solid_i})
  plot(freqs, mod(180+180/pi*phase(squeeze(FRFs(acc_dir+3*(i-1), exc_dir, :))), 360)-180);
end
plot(freqs, mod(180+180/pi*phase(squeeze(FRFs_O((solid_i-1)*6+acc_dir_O, exc_dir, :))), 360)-180, '-k');
hold off;
ylim([-180, 180]); yticks(-180:90:180);
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'xscale', 'log');

linkaxes([ax1,ax2],'x');
xlim([1, 200]);

17 Verify that we find the original FRF from the FRF in the global coordinates

From the computed FRF of the Hexapod in its 6 DOFs, compute the FRF of the accelerometer 1 fixed to the Hexapod during the measurement.

FRF_test = zeros(801, 3);
for i = 1:801
  FRF_test(i, :) = FRFs_O(31:33, 1, i) + cross(FRFs_O(34:36, 1, i), acc_pos(1, :)');
end

Figure 6: Comparison of the original measured FRFs with the recovered FRF from the FRF in the common cartesian frame

The reduction of the number of degrees of freedom from 69 (23 accelerometers with each 3DOF) to 36 (6 solid bodies with 6 DOF) seems to work well.