Add simplier model identification and comparison

index.org

@@ -778,7 +778,7 @@ Then, we extract the coordinates of the interface nodes.
 Using =K=, =M= and =int_xyz=, we can use the =Reduced Order Flexible Solid= simscape block.
 
 ** Piezoelectric parameters
-Parameters for the APA95ML:
+Parameters for the APA300ML:
 
 #+begin_src matlab
   d33 = 3e-10; % Strain constant [m/V]
@@ -790,8 +790,8 @@ Parameters for the APA95ML:
 #+end_src
 
 #+begin_src matlab
-  na = 3; % Number of stacks used as actuator
-  ns = 0; % Number of stacks used as force sensor
+  na = 2; % Number of stacks used as actuator
+  ns = 1; % Number of stacks used as force sensor
 #+end_src
 
 The ratio of the developed force to applied voltage is $d_{33} n k_a$ in [N/V].
@@ -802,7 +802,7 @@ We denote this constant by $g_a$ and:
 #+end_src
 
 #+RESULTS:
-: 1.88
+: 0.42941
 
 From cite:fleming14_desig_model_contr_nanop_system (page 123), the relation between relative displacement and generated voltage is:
 \[ V_s = \frac{d_{33}}{\epsilon^T s^D n} \Delta h \]
@@ -1357,6 +1357,181 @@ Root locus
   xlabel('Real Part'); ylabel('Imaginary Part');
 #+end_src
 
+** Identification for a simpler model
+The goal in this section is to identify the parameters of a simple APA model from the FEM.
+This can be useful is a lower order model is to be used for simulations.
+
+The presented model is based on cite:souleille18_concep_activ_mount_space_applic.
+
+The model represents the Amplified Piezo Actuator (APA) from Cedrat-Technologies (Figure [[fig:souleille18_model_piezo]]).
+The parameters are shown in the table below.
+
+#+name: fig:souleille18_model_piezo
+#+caption: Picture of an APA100M from Cedrat Technologies. Simplified model of a one DoF payload mounted on such isolator
+[[file:./figs/souleille18_model_piezo.png]]
+
+#+caption: Parameters used for the model of the APA 100M
+|       | Meaning                                                        |
+|-------+----------------------------------------------------------------|
+| $k_e$ | Stiffness used to adjust the pole of the isolator              |
+| $k_1$ | Stiffness of the metallic suspension when the stack is removed |
+| $k_a$ | Stiffness of the actuator                                      |
+| $c_1$ | Added viscous damping                                          |
+
+The goal is to determine $k_e$, $k_a$ and $k_1$ so that the simplified model fits the FEM model.
+
+Three unknowns and three equations:
+\[ \alpha = \frac{x_1}{f}(\omega=0) = \frac{\frac{k_e}{k_e + k_a}}{k_1 + \frac{k_e k_a}{k_e + k_a}} \]
+\[ \beta = \frac{x_1}{F}(\omega=0) = \frac{1}{k_1 + \frac{k_e k_a}{k_e + k_a}} \]
+\[ \gamma = \frac{dL}{f}(\omega=0) = \frac{1}{k_a + \frac{k_e k_1}{k_e + k_1}} \]
+
+#+begin_src matlab :exports none
+  m = 10;
+#+end_src
+
+#+begin_src matlab :exports none
+  %% Name of the Simulink File
+  mdl = 'APA300ML_test_bench';
+
+  %% Input/Output definition
+  clear io; io_i = 1;
+  io(io_i) = linio([mdl, '/Fd'], 1, 'openinput');  io_i = io_i + 1; % External Vertical Force [N]
+  io(io_i) = linio([mdl, '/w'],  1, 'openinput');  io_i = io_i + 1; % Base Motion [m]
+  io(io_i) = linio([mdl, '/Fa'], 1, 'openinput');  io_i = io_i + 1; % Actuator Force [N]
+  io(io_i) = linio([mdl, '/y'],  1, 'openoutput'); io_i = io_i + 1; % Vertical Displacement [m]
+  io(io_i) = linio([mdl, '/Vs'], 1, 'openoutput'); io_i = io_i + 1; % Force Sensor [V]
+  io(io_i) = linio([mdl, '/d'],  1, 'openoutput'); io_i = io_i + 1; % Stack Displacement [m]
+
+  G = linearize(mdl, io);
+
+  G.InputName = {'Fd', 'w', 'Fa'};
+  G.OutputName = {'y', 'Fs', 'd'};
+#+end_src
+
+#+begin_src matlab
+  alpha = abs(dcgain(G('y', 'Fa')));
+  beta  = abs(dcgain(G('y', 'Fd')));
+  gamma = abs(dcgain(G('d', 'Fa')));
+#+end_src
+
+#+begin_src matlab
+  na = 1;
+  ns = 2;
+#+end_src
+
+Amplification Factor
+#+begin_src matlab
+  A = abs(dcgain(G('y', 'Fa'))/dcgain(G('d', 'Fa')));
+#+end_src
+
+#+begin_src matlab
+  ka = 0.8*(1/A)/gamma;
+  ke = ka/(beta/alpha - 1);
+  k1 = 1/beta - ke*ka/(ke + ka);
+#+end_src
+
+#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
+  data2orgtable(1e-6*[ka; ke; k1], {'ka', 'ke', 'k1'}, {'Value [N/um]'}, ' %.1f ');
+#+end_src
+
+#+RESULTS:
+|    | Value [N/um] |
+|----+--------------|
+| ka |         42.9 |
+| ke |          1.5 |
+| k1 |          0.4 |
+
+Adjust the damping in the system.
+#+begin_src matlab
+  c1 = 1e2;
+#+end_src
+
+Analytical model of the simpler system:
+#+begin_src matlab
+  Ga = 1/(m*s^2 + k1 + c1*s + ke*ka/(ke + ka)) * ...
+       [ 1 ,              k1 + c1*s + ke*ka/(ke + ka)  , ke/(ke + ka) ;
+        -ke*ka/(ke + ka), ke*ka/(ke + ka)*m*s^2 ,       -ke/(ke + ka)*(m*s^2 + c1*s + k1)];
+
+  Ga.InputName = {'Fd', 'w', 'Fa'};
+  Ga.OutputName = {'y', 'Fs'};
+#+end_src
+
+Adjust the DC gain for the force sensor;
+#+begin_src matlab
+  lambda = dcgain(Ga('Fs', 'Fd'))/dcgain(G('Fs', 'Fd'));
+#+end_src
+
+#+begin_src matlab :exports none
+  freqs = logspace(0, 5, 1000);
+
+  figure;
+
+  ay = subplot(2, 3, 1);
+  title('$\displaystyle \frac{x_1}{w}$')
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(G( 'y', 'w'), freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(Ga('y', 'w'), freqs, 'Hz'))));
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Amplitude [m/m]');xlabel('Frequency [Hz]');
+
+  ax2 = subplot(2, 3, 2);
+  title('$\displaystyle \frac{x_1}{f}$')
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(G( 'y', 'Fa'), freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(Ga('y', 'Fa'), freqs, 'Hz'))));
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');
+
+  ax3 = subplot(2, 3, 3);
+  title('$\displaystyle \frac{x_1}{F}$')
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(G( 'y', 'Fd'), freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(Ga('y', 'Fd'), freqs, 'Hz'))));
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');
+
+  ax4 = subplot(2, 3, 4);
+  title('$\displaystyle \frac{F_s}{w}$')
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(lambda*G( 'Fs', 'w'), freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(Ga('Fs', 'w'), freqs, 'Hz'))));
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Amplitude [m/m]');xlabel('Frequency [Hz]');
+
+  ax5 = subplot(2, 3, 5);
+  title('$\displaystyle \frac{F_s}{f}$')
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(lambda*G( 'Fs', 'Fa'), freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(Ga('Fs', 'Fa'), freqs, 'Hz'))));
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');
+
+  ax6 = subplot(2, 3, 6);
+  title('$\displaystyle \frac{F_s}{F}$')
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(lambda*G( 'Fs', 'Fd'), freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(Ga('Fs', 'Fd'), freqs, 'Hz'))));
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
+
+  linkaxes([ax1,ax2,ax3,ax4,ax5,ax6],'x');
+#+end_src
+
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/apa300ml_comp_simpler_model.pdf', 'width', 'full', 'height', 'full');
+#+end_src
+
+#+name: fig:apa300ml_comp_simpler_model
+#+caption: Comparison of the Dynamics between the FEM model and the simplified one
+#+RESULTS:
+[[file:figs/apa300ml_comp_simpler_model.png]]
+
 * Flexible Joint
 ** Introduction                                                      :ignore: