tdehaeze/nass-micro-station-measurements
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
8ca2d8ed20
nass-micro-station-measurem.../modal-analysis/modes_analysis.org

828 lines
28 KiB
Org Mode
Raw Normal View History

Update the modal analysis
2019-07-02 17:48:34 +02:00
#+TITLE: Modal Analysis
:DRAWER:
#+STARTUP: overview
#+LANGUAGE: en
#+EMAIL: dehaeze.thomas@gmail.com
#+AUTHOR: Dehaeze Thomas
#+HTML_LINK_HOME: ../index.html
#+HTML_LINK_UP: ../index.html
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/htmlize.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/readtheorg.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/zenburn.css"/>
#+HTML_HEAD: <script type="text/javascript" src="../js/jquery.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="../js/bootstrap.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="../js/jquery.stickytableheaders.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="../js/readtheorg.js"></script>
#+HTML_MATHJAX: align: center tagside: right font: TeX
#+PROPERTY: header-args:matlab :session *MATLAB*
#+PROPERTY: header-args:matlab+ :comments org
#+PROPERTY: header-args:matlab+ :results none
#+PROPERTY: header-args:matlab+ :exports both
#+PROPERTY: header-args:matlab+ :eval no-export
#+PROPERTY: header-args:matlab+ :output-dir figs
#+PROPERTY: header-args:shell :eval no-export
#+PROPERTY: header-args:latex :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/MEGA/These/LaTeX/}{config.tex}")
#+PROPERTY: header-args:latex+ :imagemagick t :fit yes
#+PROPERTY: header-args:latex+ :iminoptions -scale 100% -density 150
#+PROPERTY: header-args:latex+ :imoutoptions -quality 100
#+PROPERTY: header-args:latex+ :results raw replace :buffer no
#+PROPERTY: header-args:latex+ :eval no-export
#+PROPERTY: header-args:latex+ :exports both
#+PROPERTY: header-args:latex+ :mkdirp yes
#+PROPERTY: header-args:latex+ :output-dir figs
:END:
* 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
* Setup
#+name: fig:nass-modal-test
#+caption: Position and orientation of the accelerometer used
[[file:figs/nass-modal-test.png]]
* Mode extraction and importation
First, we split the big =modes.asc= files into sub text files using =bash=.
#+begin_src bash :results none
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
#+end_src
Then we import them on Matlab.
#+begin_src 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));
#+end_src
We guess the number of modes identified from the length of the imported data.
#+begin_src matlab
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
#+end_src
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.
#+begin_src matlab
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
#+end_src
The obtained mode frequencies and damping are shown below.
#+begin_src matlab :exports both :results value table replace :post addhdr(*this*)
data2orgtable([freqs, damps], {}, {'Frequency [Hz]', 'Damping [%]'}, ' %.1f ');
#+end_src
#+RESULTS:
| 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 |
* Positions of the sensors
We process the file exported from the =modal= software containing the positions of the sensors using =bash=.
#+begin_src bash :results none
Add few notes about modal complexity
2019-07-03 09:08:35 +02:00
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
Update the modal analysis
2019-07-02 17:48:34 +02:00
#+end_src
We then import that on =matlab=, and sort them.
#+begin_src matlab
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);
#+end_src
The positions of the sensors relative to the point of interest are shown below.
#+begin_src matlab :exports both :results value table replace :post addhdr(*this*)
data2orgtable(1000*acc_pos, {}, {'x [mm]', 'y [mm]', 'z [mm]'}, ' %.0f ');
#+end_src
#+RESULTS:
| 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 |
* 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 [[fig:nass-modal-test]]).
#+begin_src matlab
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);
#+end_src
* From local coordinates to global coordinates
#+begin_src latex :file local_to_global_coordinates.pdf :post pdf2svg(file=*this*, ext="png") :exports results
\newcommand\irregularcircle[2]{% radius, irregularity
\pgfextra {\pgfmathsetmacro\len{(#1)+rand*(#2)}}
+(0:\len pt)
\foreach \a in {10,20,...,350}{
\pgfextra {\pgfmathsetmacro\len{(#1)+rand*(#2)}}
-- +(\a:\len pt)
} -- cycle
}
\begin{tikzpicture}
\draw[rounded corners=1mm] (0, 0) \irregularcircle{3cm}{1mm};
\node[] (origin) at (4, -1) {$\bullet$};
\begin{scope}[shift={(origin)}]
\def\axissize{0.8cm}
\draw[->] (0, 0) -- ++(\axissize, 0) node[above left]{$x$};
\draw[->] (0, 0) -- ++(0, \axissize) node[below right]{$y$};
\draw[fill, color=black] (0, 0) circle (0.05*\axissize);
\node[draw, circle, inner sep=0pt, minimum size=0.4*\axissize, label=left:$z$] (yaxis) at (0, 0){};
\node[below right] at (0, 0){$\{O\}$};
\end{scope}
\coordinate[] (p1) at (-1.5, -1.5);
\coordinate[] (p2) at (-1.5, 1.5);
\coordinate[] (p3) at ( 1.5, 1.5);
\coordinate[] (p4) at ( 1.5, -1.5);
\draw[->] (p1)node[]{$\bullet$}node[above]{$p_1$} -- ++(1, 0.5)node[right]{$v_1$};
\draw[->] (p2)node[]{$\bullet$}node[above]{$p_2$} -- ++(-0.5, 1)node[right]{$v_2$};
\draw[->] (p3)node[]{$\bullet$}node[above]{$p_3$} -- ++(1, 0.5)node[right]{$v_3$};
\draw[->] (p4)node[]{$\bullet$}node[above]{$p_4$} -- ++(0.5, 1)node[right]{$v_4$};
\end{tikzpicture}
#+end_src
#+RESULTS:
[[file:figs/local_to_global_coordinates.png]]
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.
#+begin_src matlab
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
#+end_src
* 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?
#+begin_src matlab
eigen_value_M = diag(freqs*2*pi);
eigen_vector_M = reshape(mode_shapes_O, [mod_n, 6*length(solid_names)])';
#+end_src
\[ \{\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 \]
* Modal Complexity
Add few notes about modal complexity
2019-07-03 09:08:35 +02:00
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 [[fig:modal_complexity_small]].
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 [[fig:modal_complexity_small]] whereas an highly complex mode is shown on figure [[fig:modal_complexity_high]].
The complexity of all the modes are compared on figure [[fig:modal_complexities]].
Add importation of FRF and COH
2019-07-03 11:57:59 +02:00
#+begin_src matlab :exports none
Update the modal analysis
2019-07-02 17:48:34 +02:00
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');
Add few notes about modal complexity
2019-07-03 09:08:35 +02:00
title(sprintf('Mode %i', mod_i));
Update the modal analysis
2019-07-02 17:48:34 +02:00
axis manual equal
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
Add few notes about modal complexity
2019-07-03 09:08:35 +02:00
#+begin_src matlab :var filepath="figs/modal_complexity_small.pdf" :var figsize="normal-normal" :post pdf2svg(file=*this*, ext="png")
Update the modal analysis
2019-07-02 17:48:34 +02:00
<<plt-matlab>>
#+end_src
Add few notes about modal complexity
2019-07-03 09:08:35 +02:00
#+NAME: fig:modal_complexity_small
#+CAPTION: Modal Complexity of one mode with small complexity
[[file:figs/modal_complexity_small.png]]
Update the modal analysis
2019-07-02 17:48:34 +02:00
Add importation of FRF and COH
2019-07-03 11:57:59 +02:00
#+begin_src matlab :exports none
Add few notes about modal complexity
2019-07-03 09:08:35 +02:00
mod_i = 8;
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');
title(sprintf('Mode %i', mod_i));
axis manual equal
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/modal_complexity_high.pdf" :var figsize="normal-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:modal_complexity_high
#+CAPTION: Modal Complexity of one higly complex mode
[[file:figs/modal_complexity_high.png]]
Add importation of FRF and COH
2019-07-03 11:57:59 +02:00
#+begin_src matlab :exports none
Update the modal analysis
2019-07-02 17:48:34 +02:00
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');
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/modal_complexities.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:modal_complexities
#+CAPTION: Modal complexity for each mode
[[file:figs/modal_complexities.png]]
* 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 | |
#+TBLFM: $4=vsum($2..$3)
#+TBLFM: @>$2..$>=vsum(@I..@II)
* TODO Normalization of mode shapes?
We normalize each column of the eigen vector matrix.
Then, each eigenvector as a norm of 1.
#+begin_src matlab
eigen_vector_M = eigen_vector_M./vecnorm(eigen_vector_M);
#+end_src
* 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).
#+begin_src matlab
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();
#+end_src
#+begin_src matlab
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();
#+end_src
#+begin_src matlab
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, :)));
#+end_src
Add importation of FRF and COH
2019-07-03 11:57:59 +02:00
* 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
#+begin_src matlab
a = load('./modal_analysis/frf_coh/Measurement1.mat');
#+end_src
#+begin_src matlab
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]);
#+end_src
* 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}
#+begin_src matlab
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;
#+end_src
* Analysis of some FRFs
#+begin_src matlab
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]);
#+end_src
#+begin_src matlab
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 [\%]');
#+end_src
Composite Response Function.
We here sum the norm instead of the complex numbers.
#+begin_src matlab
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));
#+end_src
#+begin_src matlab
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]);
#+end_src
* From local coordinates to global coordinates with the FRFs
#+begin_src matlab
% 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
#+end_src
* Analysis of some FRF in the global coordinates
#+begin_src matlab
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]);
#+end_src
* Compare global coordinates to local coordinates
#+begin_src matlab
solid_i = 1;
FRF from accelerometers to global cartesian and comparison
2019-07-03 13:40:02 +02:00
acc_dir_O = 6;
Add importation of FRF and COH
2019-07-03 11:57:59 +02:00
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
FRF from accelerometers to global cartesian and comparison
2019-07-03 13:40:02 +02:00
plot(freqs, abs(squeeze(FRFs_O((solid_i-1)*6+acc_dir_O, exc_dir, :))), '-k');
Add importation of FRF and COH
2019-07-03 11:57:59 +02:00
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
FRF from accelerometers to global cartesian and comparison
2019-07-03 13:40:02 +02:00
plot(freqs, mod(180+180/pi*phase(squeeze(FRFs_O((solid_i-1)*6+acc_dir_O, exc_dir, :))), 360)-180, '-k');
Add importation of FRF and COH
2019-07-03 11:57:59 +02:00
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]);
#+end_src
FRF from accelerometers to global cartesian and comparison
2019-07-03 13:40:02 +02:00
* 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.
#+begin_src matlab
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
#+end_src
#+begin_src matlab :exports none
figure;
ax1 = subplot(3, 1, 1);
hold on;
plot(freqs, abs(squeeze(FRFs(1, 1, :))));
plot(freqs, abs(squeeze(FRF_test(:, 1))), '--k');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
set(gca, 'XTickLabel',[]);
xlim([1, 200]);
title('FRF $\frac{D_{1x}}{F_x}$');
ax2 = subplot(3, 1, 2);
hold on;
plot(freqs, abs(squeeze(FRFs(2, 1, :))));
plot(freqs, abs(squeeze(FRF_test(:, 2))), '--k');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
set(gca, 'XTickLabel',[]);
ylabel('Amplitude');
xlim([1, 200]);
title('FRF $\frac{D_{1y}}{F_x}$');
ax3 = subplot(3, 1, 3);
hold on;
plot(freqs, abs(squeeze(FRFs(3, 1, :))));
plot(freqs, abs(squeeze(FRF_test(:, 3))), '--k');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]');
xlim([1, 200]);
legend({'Original Measurement', 'Recovered Measurement'}, 'Location', 'southeast');
title('FRF $\frac{D_{1z}}{F_x}$');
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/compare_original_meas_with_recovered.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:compare_original_meas_with_recovered
#+CAPTION: Comparison of the original measured FRFs with the recovered FRF from the FRF in the common cartesian frame
[[file:figs/compare_original_meas_with_recovered.png]]
#+begin_important
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.
#+end_important
Add importation of FRF and COH
2019-07-03 11:57:59 +02:00
Update the modal analysis
2019-07-02 17:48:34 +02:00
* TODO Synthesis of FRF curves
Reference in New Issue Copy Permalink