nass-micro-station-measurements/dynamical-meas-granite/index.org

#+TITLE: Vibrations induced by the Slip-Ring and the Spindle
: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="https://research.tdehaeze.xyz/css/style.css"/>
#+HTML_HEAD: <script type="text/javascript" src="https://research.tdehaeze.xyz/js/script.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
:END:

* Experimental Setup
*Setup*:
One geophone is located on the marble, the other one on the floor (see figure [[fig:experimental_setup]]).


Two geophone are use:
- One on the floor (corresponding to the first column in the data)
- One at the marble location (corresponding to the second column in the data)

Two voltage amplifiers are used, their setup is:
- gain of 60dB
- AC/DC switch on AC
- Low pass filter at 1kHz

A first order low pass filter is also added at the input of the voltage amplifiers.

*Goal*:
- Identify the marble dynamics in all the directions


*Measurements*:
Three measurements are done:
| Measurement File   | Description |
|--------------------+-------------|
| =mat/data_037.mat= | Z direction |
| =mat/data_038.mat= | N direction |
| =mat/data_039.mat= | E direction |

Each of the measurement =mat= file contains one =data= array with 3 columns:
| Column number | Description       |
|---------------+-------------------|
|             1 | Geophone - Floor  |
|             2 | Geophone - Marble |
|             3 | Time              |

#+name: fig:experimental_setup
#+caption: Picture of the experimental setup
#+attr_html: :width 500px
[[file:./img/IMG_20190513_161729.jpg]]

#+name: fig:experimental_setup_bix
#+caption: Picture of the experimental setup
#+attr_html: :width 500px
[[file:./img/IMG_20190513_161718.jpg]]

* Data Analysis
  :PROPERTIES:
  :header-args:matlab+: :tangle matlab/marble_dynamics.m
  :header-args:matlab+: :comments org :mkdirp yes
  :END:
  <<sec:marble_dynamics>>

** ZIP file containing the data and matlab files                     :ignore:
#+begin_src bash :exports none :results none
  if [ matlab/marble_dynamics.m -nt data/marble_dynamics.zip ]; then
    cp matlab/marble_dynamics.m marble_dynamics.m;
    zip data/marble_dynamics \
        mat/data_037.mat \
        mat/data_038.mat \
        mat/data_039.mat \
        marble_dynamics.m
    rm marble_dynamics.m;
  fi
#+end_src

#+begin_note
  All the files (data and Matlab scripts) are accessible [[file:data/marble_dynamics.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 data
#+begin_src matlab
  m_z = load('mat/data_037.mat', 'data'); m_z = m_z.data;
  m_n = load('mat/data_038.mat', 'data'); m_n = m_n.data;
  m_e = load('mat/data_039.mat', 'data'); m_e = m_e.data;
#+end_src

** Time domain plots
#+begin_src matlab
  figure;
  subplot(1, 3, 1);
  hold on;
  plot(m_z(:, 3), m_z(:, 2), 'DisplayName', 'Marble - Z');
  plot(m_z(:, 3), m_z(:, 1), 'DisplayName', 'Floor - Z');
  hold off;
  xlabel('Time [s]'); ylabel('Voltage [V]');
  xlim([0, 100]); ylim([-2 2]);
  legend('Location', 'northeast');

  subplot(1, 3, 2);
  hold on;
  plot(m_n(:, 3), m_n(:, 2), 'DisplayName', 'Marble - N');
  plot(m_n(:, 3), m_n(:, 1), 'DisplayName', 'Floor - N');
  hold off;
  xlabel('Time [s]'); ylabel('Voltage [V]');
  xlim([0, 100]); ylim([-2 2]);
  legend('Location', 'northeast');

  subplot(1, 3, 3);
  hold on;
  plot(m_e(:, 3), m_e(:, 2), 'DisplayName', 'Marble - E');
  plot(m_e(:, 3), m_e(:, 1), 'DisplayName', 'Floor - E');
  hold off;
  xlabel('Time [s]'); ylabel('Voltage [V]');
  xlim([0, 100]); ylim([-2 2]);
  legend('Location', 'northeast');
#+end_src

#+NAME: fig:marble_floor_motion_time
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/marble_floor_motion_time.pdf" :var figsize="full-normal" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:marble_floor_motion_time
#+CAPTION: Floor and ground motion
#+RESULTS: fig:marble_floor_motion_time
[[file:figs/marble_floor_motion_time.png]]

** Compute the power spectral densities
We first compute some parameters that will be used for the PSD computation.
#+begin_src matlab :results none
  dt = m_z(2, 3)-m_z(1, 3);

  Fs = 1/dt; % [Hz]

  win = hanning(ceil(10*Fs));
#+end_src

Then we compute the Power Spectral Density using =pwelch= function.
#+begin_src matlab
  [px_fz, f] = pwelch(m_z(:, 1), win, [], [], Fs);
  [px_gz, ~] = pwelch(m_z(:, 2), win, [], [], Fs);

  [px_fn, ~] = pwelch(m_n(:, 1), win, [], [], Fs);
  [px_gn, ~] = pwelch(m_n(:, 2), win, [], [], Fs);

  [px_fe, ~] = pwelch(m_e(:, 1), win, [], [], Fs);
  [px_ge, ~] = pwelch(m_e(:, 2), win, [], [], Fs);
#+end_src


The results are shown on figure [[fig:floor_marble_psd_z]] for the Z direction, figure [[fig:floor_marble_psd_n]] for the north direction, and figure [[fig:floor_marble_psd_e]] for the east direction.

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(f, sqrt(px_fz), 'DisplayName', 'Floor - Z');
  plot(f, sqrt(px_gz), 'DisplayName', 'Granite - Z');
  hold off;
  set(gca, 'xscale', 'log');
  set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('ASD of the measured Voltage $\left[\frac{V}{\sqrt{Hz}}\right]$')
  legend('Location', 'southwest');
  xlim([0.1, 500]);
#+end_src

#+NAME: fig:floor_marble_psd_z
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/floor_marble_psd_z.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:floor_marble_psd_z
#+CAPTION: Amplitude Spectral Density of the measured voltage corresponding to the geophone located on the floor and on the marble - Z direction
#+RESULTS: fig:floor_marble_psd_z
[[file:figs/floor_marble_psd_z.png]]

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(f, sqrt(px_fn), 'DisplayName', 'Floor - N');
  plot(f, sqrt(px_gn), 'DisplayName', 'Granite - N');
  hold off;
  set(gca, 'xscale', 'log');
  set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('ASD of the measured Voltage $\left[\frac{V}{\sqrt{Hz}}\right]$')
  legend('Location', 'southwest');
  xlim([0.1, 500]);
#+end_src

#+NAME: fig:floor_marble_psd_n
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/floor_marble_psd_n.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:floor_marble_psd_n
#+CAPTION: Amplitude Spectral Density of the measured voltage corresponding to the geophone located on the floor and on the marble - N direction
#+RESULTS: fig:floor_marble_psd_n
[[file:figs/floor_marble_psd_n.png]]

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(f, sqrt(px_fe), 'DisplayName', 'Floor - E');
  plot(f, sqrt(px_ge), 'DisplayName', 'Granite - E');
  hold off;
  set(gca, 'xscale', 'log');
  set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('ASD of the measured Voltage $\left[\frac{V}{\sqrt{Hz}}\right]$')
  legend('Location', 'southwest');
  xlim([0.1, 500]);
#+end_src

#+NAME: fig:floor_marble_psd_e
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/floor_marble_psd_e.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:floor_marble_psd_e
#+CAPTION: Amplitude Spectral Density of the measured voltage corresponding to the geophone located on the floor and on the marble - E direction
#+RESULTS: fig:floor_marble_psd_e
[[file:figs/floor_marble_psd_e.png]]

** Compute the transfer function from floor motion to ground motion
We now compute the transfer function from the floor motion to the granite motion.

The result is shown on figure [[fig:tf_granite]].
#+begin_src matlab :results none
  [TZ, f] = tfestimate(m_z(:, 1), -m_z(:, 2), win, [], [], Fs);
  [TN, ~] = tfestimate(m_n(:, 1), -m_n(:, 2), win, [], [], Fs);
  [TE, ~] = tfestimate(m_e(:, 1), -m_e(:, 2), win, [], [], Fs);
#+end_src

#+begin_src matlab :results none :exports none
  figure;
  ax1 = subplot(2, 1, 1);
  hold on;
  plot(f, abs(TZ), 'DisplayName', 'Z');
  plot(f, abs(TN), 'DisplayName', 'N');
  plot(f, abs(TE), 'DisplayName', 'E');
  hold off;
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
  set(gca, 'XTickLabel',[]);
  ylabel('Magnitude');
  legend('Location', 'southwest');

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(f, mod(180+180/pi*phase(TZ), 360)-180);
  plot(f, mod(180+180/pi*phase(TN), 360)-180);
  plot(f, mod(180+180/pi*phase(TE), 360)-180);
  hold off;
  set(gca, 'xscale', 'log');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');

  linkaxes([ax1,ax2],'x');
  xlim([10, 100]);
#+end_src

#+NAME: fig:tf_granite
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/tf_granite.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:tf_granite
#+CAPTION: Transfer function from floor motion to granite motion
#+RESULTS: fig:tf_granite
[[file:figs/tf_granite.png]]

** Conclusion
#+begin_important
  - We see resonance of the granite at 33Hz in the horizontal directions
  - We see two resonances for the z direction: at 60Hz and 75Hz
#+end_important