Modal Analysis
+Table of Contents
+-
+
- 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. 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
++Complexity of one mode +
+mod_i = 1; +i_max = convhull(real(eigen_vector_M(:, mod_i)), imag(eigen_vector_M(:, mod_i))); +radius = max(abs(eigen_vector_M(:, mod_i))); +theta = linspace(0, 2*pi, 100); + +figure; +hold on; +plot(radius*cos(theta), radius*sin(theta), '-'); +plot(real(eigen_vector_M(i_max, mod_i)), imag(eigen_vector_M(i_max, mod_i)), '-'); +plot(real(eigen_vector_M(:, mod_i)), imag(eigen_vector_M(:, mod_i)), 'ko'); +hold off; +xlabel('Real Part'); ylabel('Imaginary Part'); +axis manual equal ++
+
Figure 3: Modal Complexity of one mode
++Complexity function of the mode order. +
+modes_complexity = zeros(mod_n, 1); +for mod_i = 1:mod_n + i = convhull(real(eigen_vector_M(:, mod_i)), imag(eigen_vector_M(:, mod_i))); + area_complex = polyarea(real(eigen_vector_M(i, mod_i)), imag(eigen_vector_M(i, mod_i))); + area_circle = pi*max(abs(eigen_vector_M(:, mod_i)))^2; + modes_complexity(mod_i) = area_complex/area_circle; +end + +figure; +plot(1:mod_n, modes_complexity, 'ok'); +ylim([0, 1]); +xlabel('Mode Number'); ylabel('Modal Complexity'); ++
+
Figure 4: 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, :))); ++