diff --git a/figs/apa300ml_comp_simpler_model.pdf b/figs/apa300ml_comp_simpler_model.pdf index 059ce03..1141940 100644 Binary files a/figs/apa300ml_comp_simpler_model.pdf and b/figs/apa300ml_comp_simpler_model.pdf differ diff --git a/figs/apa300ml_comp_simpler_model.png b/figs/apa300ml_comp_simpler_model.png index a13fd9f..493b5d5 100644 Binary files a/figs/apa300ml_comp_simpler_model.png and b/figs/apa300ml_comp_simpler_model.png differ diff --git a/index.html b/index.html index 3029f89..0e7e6ce 100644 --- a/index.html +++ b/index.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Finite Element Model with Simscape @@ -34,60 +34,61 @@

Table of Contents

-
-

1 Amplified Piezoelectric Actuator - 3D elements

+
+

1 Amplified Piezoelectric Actuator - 3D elements

The idea here is to: @@ -101,8 +102,8 @@ The idea here is to:

-
-

1.1 Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates

+
+

1.1 Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates

We first extract the stiffness and mass matrices. @@ -128,8 +129,8 @@ Then, we extract the coordinates of the interface nodes.

-
-

1.2 Output parameters

+
+

1.2 Output parameters

load('./mat/piezo_amplified_3d.mat', 'int_xyz', 'int_i', 'n_xyz', 'n_i', 'nodes', 'M', 'K');
@@ -168,7 +169,7 @@ Then, we extract the coordinates of the interface nodes.
 
 
 
-

+

 
amplified_piezo_interface_nodes.png
 

 
Figure 1: Interface Nodes for the Amplified Piezo Actuator

@@ -626,8 +627,8 @@ Using K, M and int_xyz, we can use the 
 
 
-

-
1.3 Piezoelectric parameters

+

+
1.3 Piezoelectric parameters

 

 

 Parameters for the APA95ML:
@@ -688,8 +689,8 @@ where:
 

 

 
-

-
1.4 Identification of the Dynamics

+

+
1.4 Identification of the Dynamics

 

 

 The flexible element is imported using the Reduced Order Flexible Solid simscape block.
@@ -741,7 +742,7 @@ And the dynamics is identified.
 

 
 

-The two identified dynamics are compared in Figure 2.
+The two identified dynamics are compared in Figure 2.
 

 

 
%% Name of the Simulink File
@@ -757,7 +758,7 @@ Ghm = -linearize(mdl, io);
 

 
 
-

+

 
dynamics_act_disp_comp_mass.png
 

 
Figure 2: Dynamics from \(F\) to \(d\) without a payload and with a 10kg payload

@@ -765,8 +766,8 @@ Ghm = -linearize(mdl, io);
 

 

 
-

-
1.5 Comparison with Ansys

+

+
1.5 Comparison with Ansys

 

 

 Let’s import the results from an Harmonic response analysis in Ansys.
@@ -778,11 +779,11 @@ Gresp10 = readtable('FEA_HarmResponse_10kg.txt');
 

 
 

-The obtained dynamics from the Simscape model and from the Ansys analysis are compare in Figure 3.
+The obtained dynamics from the Simscape model and from the Ansys analysis are compare in Figure 3.
 

 
 
-

+

 
dynamics_force_disp_comp_anasys.png
 

 
Figure 3: Comparison of the obtained dynamics using Simscape with the harmonic response analysis using Ansys

@@ -790,15 +791,15 @@ The obtained dynamics from the Simscape model and from the Ansys analysis are co
 

 

 
-

-
1.6 Force Sensor

+

+
1.6 Force Sensor

 

 

 The dynamics is identified from internal forces applied between nodes 3 and 11 to the relative displacement of nodes 11 and 13.
 

 
 

-The obtained dynamics is shown in Figure 4.
+The obtained dynamics is shown in Figure 4.
 

 
 

@@ -838,7 +839,7 @@ Gfm = linearize(mdl, io);
 

 
 
-

+

 
dynamics_force_force_sensor_comp_mass.png
 

 
Figure 4: Dynamics from \(F\) to \(F_m\) for \(m=0\) and \(m = 10kg\)

@@ -846,8 +847,8 @@ Gfm = linearize(mdl, io);
 

 

 
-

-
1.7 Distributed Actuator

+

+
1.7 Distributed Actuator

 

 

 
m = 0;
@@ -896,8 +897,8 @@ Gdm = linearize(mdl, io);
 

 

 
-

-
1.8 Distributed Actuator and Force Sensor

+

+
1.8 Distributed Actuator and Force Sensor

 

 

 
m = 0;
@@ -937,8 +938,8 @@ Gfdm = linearize(mdl, io);
 

 

 
-

-
1.9 Dynamics from input voltage to displacement

+

+
1.9 Dynamics from input voltage to displacement

 

 

 
m = 5;
@@ -950,7 +951,7 @@ And the dynamics is identified.
 

 
 

-The two identified dynamics are compared in Figure 2.
+The two identified dynamics are compared in Figure 2.
 

 

 
%% Name of the Simulink File
@@ -972,8 +973,8 @@ G = -linearize(mdl, io);
 

 

 
-

-
1.10 Dynamics from input voltage to output voltage

+

+
1.10 Dynamics from input voltage to output voltage

 

 

 
m = 5;
@@ -996,19 +997,19 @@ G = -linearize(mdl, io);
 

 

 
-

-
2 APA300ML

+

+
2 APA300ML

 

 
-

+

 
apa300ml_ansys.jpg
 

 
Figure 5: Ansys FEM of the APA300ML

 

 

 
-

-
2.1 Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates

+

+
2.1 Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates

 

 

 We first extract the stiffness and mass matrices.
@@ -1040,8 +1041,8 @@ Then, we extract the coordinates of the interface nodes.
 

 

 
-

-
2.2 Output parameters

+

+
2.2 Output parameters

 

 

 
load('./mat/APA300ML.mat', 'int_xyz', 'int_i', 'n_xyz', 'n_i', 'nodes', 'M', 'K');
@@ -1482,11 +1483,11 @@ Using K, M and int_xyz, we can use the 
 

 
-

-
2.3 Piezoelectric parameters

+

+
2.3 Piezoelectric parameters

 

 

-Parameters for the APA95ML:
+Parameters for the APA300ML:
 

 
 

@@ -1500,8 +1501,8 @@ C = 5e-6; % Stack capactiance [F]
 

 
 

-
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
 

 

 
@@ -1516,7 +1517,7 @@ We denote this constant by \(g_a\) and:
 

 
 
-1.88
+0.42941
 

 
 
@@ -1544,12 +1545,12 @@ where:
 

 

 
-

-
2.4 Identification of the APA Characteristics

+

+
2.4 Identification of the APA Characteristics

 

 

-

-
2.4.1 Stiffness

+

+
2.4.1 Stiffness

 

 

 The transfer function from vertical external force to the relative vertical displacement is identified.
@@ -1574,16 +1575,16 @@ The specified stiffness in the datasheet is \(k = 1.8\, [N/\mu m]\).
 

 

 
-

-
2.4.2 Resonance Frequency

+

+
2.4.2 Resonance Frequency

 

 

 The resonance frequency is specified to be between 650Hz and 840Hz.
-This is also the case for the FEM model (Figure 6).
+This is also the case for the FEM model (Figure 6).
 

 
 
-

+

 
apa300ml_resonance.png
 

 
Figure 6: First resonance is around 800Hz

@@ -1591,8 +1592,8 @@ This is also the case for the FEM model (Figure 6).
 

 

 
-

-
2.4.3 Amplification factor

+

+
2.4.3 Amplification factor

 

 

 The amplification factor is the ratio of the axial displacement to the stack displacement.
@@ -1625,8 +1626,8 @@ If we take the ratio of the piezo height and length (approximation of the amplif
 

 

 
-

-
2.4.4 Stroke

+

+
2.4.4 Stroke

 

 

 Estimation of the actuator stroke:
@@ -1657,8 +1658,8 @@ This is exactly the specified stroke in the data-sheet.
 

 

 
-

-
2.5 Identification of the Dynamics

+

+
2.5 Identification of the Dynamics

 

 

 The flexible element is imported using the Reduced Order Flexible Solid simscape block.
@@ -1684,7 +1685,7 @@ The same dynamics is identified for a payload mass of 10Kg.
 

 
 
-

+

 
apa300ml_plant_dynamics.png
 

 
Figure 7: Transfer function from forces applied by the stack to the axial displacement of the APA

@@ -1692,29 +1693,29 @@ The same dynamics is identified for a payload mass of 10Kg.
 

 

 
-

-
2.6 IFF

+

+
2.6 IFF

 

 

 Let’s use 2 stacks as actuators and 1 stack as force sensor.
 

 
 

-The transfer function from actuator to sensors is identified and shown in Figure 8.
+The transfer function from actuator to sensors is identified and shown in Figure 8.
 

 
-

+

 
apa300ml_iff_plant.png
 

 
Figure 8: Transfer function from actuator to force sensor

 

 
 

-For root locus corresponding to IFF is shown in Figure 9.
+For root locus corresponding to IFF is shown in Figure 9.
 

 
 
-

+

 
apa300ml_iff_root_locus.png
 

 
Figure 9: Root Locus for IFF

@@ -1722,46 +1723,221 @@ For root locus corresponding to IFF is shown in Figure 9
 

 
-

-
2.7 DVF

+

+
2.7 DVF

 

 

-Now the dynamics from the stack actuator to the relative motion sensor is identified and shown in Figure 10.
+Now the dynamics from the stack actuator to the relative motion sensor is identified and shown in Figure 10.
 

 
-

+

 
apa300ml_dvf_plant.png
 

 
Figure 10: Transfer function from stack actuator to relative motion sensor

 

 
 

-The root locus for DVF is shown in Figure 11.
+The root locus for DVF is shown in Figure 11.
 

 
 
-

+

 
apa300ml_dvf_root_locus.png
 

 
Figure 11: Root Locus for Direct Velocity Feedback

 

 

 

+
+

+
2.8 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 (Souleille et al. 2018).
+

+
+

+The model represents the Amplified Piezo Actuator (APA) from Cedrat-Technologies (Figure 12).
+The parameters are shown in the table below.
+

+
+
+

+
souleille18_model_piezo.png
+

+
Figure 12: Picture of an APA100M from Cedrat Technologies. Simplified model of a one DoF payload mounted on such isolator

 

 
-

-
3 Flexible Joint

+
+
+
+++
++
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Table 7: 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.
+

+
+

+\[ \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}} \]
+

+
+

+If we can fix \(k_a\), we can determine \(k_e\) and \(k_1\) with:
+\[ k_e = \frac{k_a}{\frac{\beta}{\alpha} - 1} \]
+\[ k_1 = \frac{1}{\beta} - \frac{k_e k_a}{k_e + k_a} \]
+

+
+

+From the identified dynamics, compute \(\alpha\) and \(\beta\)
+

+

+
alpha = abs(dcgain(G('y', 'Fa')));
+beta  = abs(dcgain(G('y', 'Fd')));
+

+

+
+

+\(k_a\) is estimated using the following formula:
+

+

+
ka = 0.8/abs(dcgain(G('y', 'Fa')));
+

+

+

+The factor can be adjusted to better match the curves.
+

+
+

+Then \(k_e\) and \(k_1\) are computed.
+

+

+
ke = ka/(beta/alpha - 1);
+k1 = 1/beta - ke*ka/(ke + ka);
+

+

+
+
+
+
+++
++
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+




 Value [N/um]
ka42.9
ke1.5
k10.4

+
+

+The damping in the system is adjusted to match the FEM model if necessary.
+

+

+
c1 = 1e2;
+

+

+
+

+Analytical model of the simpler system:
+

+

+
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'};
+

+

+
+

+Adjust the DC gain for the force sensor:
+

+

+
lambda = dcgain(Ga('Fs', 'Fd'))/dcgain(G('Fs', 'Fd'));
+

+

+
+
+

+
apa300ml_comp_simpler_model.png
+

+
Figure 13: Comparison of the Dynamics between the FEM model and the simplified one

+

+

+

+

+
+

+
3 Flexible Joint

 

 
-

+

 
flexor_id16_screenshot.png
 

-
Figure 12: Flexor studied

+
Figure 14: Flexor studied

 

 

 
-

-
3.1 Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates

+

+
3.1 Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates

 

 

 We first extract the stiffness and mass matrices.
@@ -1787,8 +1963,8 @@ Then, we extract the coordinates of the interface nodes.
 

 

 
-

-
3.2 Output parameters

+

+
3.2 Output parameters

 

 

 
load('./mat/flexor_ID16.mat', 'int_xyz', 'int_i', 'n_xyz', 'n_i', 'nodes', 'M', 'K');
@@ -1827,7 +2003,7 @@ Then, we extract the coordinates of the interface nodes.
 
 
 
-
+@@ -1869,7 +2045,7 @@ Then, we extract the coordinates of the interface nodes.
 
Table 7: Coordinates of the interface nodesTable 8: Coordinates of the interface nodes
 
 
 


 
 
-
+@@ -2027,7 +2203,7 @@ Then, we extract the coordinates of the interface nodes.
 
 
 
Table 8: First 10x10 elements of the Stiffness matrixTable 9: First 10x10 elements of the Stiffness matrix
 
 
 


-
+@@ -2189,8 +2365,8 @@ Using K, M and int_xyz, we can use the 
 
 
-

-
3.3 Flexible Joint Characteristics

+

+
3.3 Flexible Joint Characteristics

 

 
Table 9: First 10x10 elements of the Mass matrixTable 10: First 10x10 elements of the Mass matrix
 
 
 


 
@@ -2238,8 +2414,8 @@ Using K, M and int_xyz, we can use the 
 
 
-

-
3.4 Identification

+

+
3.4 Identification

 

 

 
m = 10;
@@ -2309,16 +2485,16 @@ G = linearize(mdl, io);
 

 

 
-

-
4 Integral Force Feedback with Amplified Piezo

+

+
4 Integral Force Feedback with Amplified Piezo

 

 

 In this section, we try to replicate the results obtained in (Souleille et al. 2018).
 

 

 
-

-
4.1 Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates

+

+
4.1 Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates

 

 

 We first extract the stiffness and mass matrices.
@@ -2339,11 +2515,11 @@ Then, we extract the coordinates of the interface nodes.
 

 

 
-

-
4.2 IFF Plant

+

+
4.2 IFF Plant

 

 

-The transfer function from the force actuator to the force sensor is identified and shown in Figure 13.
+The transfer function from the force actuator to the force sensor is identified and shown in Figure 15.
 

 
 

@@ -2380,19 +2556,19 @@ Gf = linearize(mdl, io);
 

 
 
-

+

 
piezo_amplified_iff_plant.png
 

-
Figure 13: IFF Plant

+
Figure 15: IFF Plant

 

 

 

 
-

-
4.3 IFF controller

+

+
4.3 IFF controller

 

 

-The controller is defined and the loop gain is shown in Figure 14.
+The controller is defined and the loop gain is shown in Figure 16.
 

 

 
Kiff = -1e12/s;
@@ -2400,16 +2576,16 @@ The controller is defined and the loop gain is shown in Figure 
+

 
piezo_amplified_iff_loop_gain.png
 

-
Figure 14: IFF Loop Gain

+
Figure 16: IFF Loop Gain

 

 

 

 
-

-
4.4 Closed Loop System

+

+
4.4 Closed Loop System

 

 

 
m = 10;
@@ -2452,18 +2628,18 @@ G.OutputName  = {'x1', 'Fs'};
 

 
 
-
@@ -2481,7 +2657,7 @@ G.OutputName  = {'x1', 'Fs'};
 

 

 
Author: Dehaeze Thomas

-
Created: 2020-08-03 lun. 15:42

+
Created: 2020-08-04 mar. 11:53

 

 
 
diff --git a/index.org b/index.org
index 8a4f9a8..12d71b1 100644
--- a/index.org
+++ b/index.org
@@ -1380,10 +1380,12 @@ The parameters are shown in the table below.
 
 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}} \]
+
+If we can fix $k_a$, we can determine $k_e$ and $k_1$ with:
+\[ k_e = \frac{k_a}{\frac{\beta}{\alpha} - 1} \]
+\[ k_1 = \frac{1}{\beta} - \frac{k_e k_a}{k_e + k_a} \]
 
 #+begin_src matlab :exports none
   m = 10;
@@ -1408,24 +1410,20 @@ Three unknowns and three equations:
   G.OutputName = {'y', 'Fs', 'd'};
 #+end_src
 
+From the identified dynamics, compute $\alpha$ and $\beta$
 #+begin_src matlab
   alpha = abs(dcgain(G('y', 'Fa')));
   beta  = abs(dcgain(G('y', 'Fd')));
-  gamma = abs(dcgain(G('d', 'Fa')));
 #+end_src
 
+$k_a$ is estimated using the following formula:
 #+begin_src matlab
-  na = 1;
-  ns = 2;
+  ka = 0.8/abs(dcgain(G('y', 'Fa')));
 #+end_src
+The factor can be adjusted to better match the curves.
 
-Amplification Factor
+Then $k_e$ and $k_1$ are computed.
 #+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
@@ -1441,7 +1439,7 @@ Amplification Factor
 | ke |          1.5 |
 | k1 |          0.4 |
 
-Adjust the damping in the system.
+The damping in the system is adjusted to match the FEM model if necessary.
 #+begin_src matlab
   c1 = 1e2;
 #+end_src
@@ -1456,7 +1454,7 @@ Analytical model of the simpler system:
   Ga.OutputName = {'y', 'Fs'};
 #+end_src
 
-Adjust the DC gain for the force sensor;
+Adjust the DC gain for the force sensor:
 #+begin_src matlab
   lambda = dcgain(Ga('Fs', 'Fd'))/dcgain(G('Fs', 'Fd'));
 #+end_src