diff --git a/figs/APA95ML_nodes_1.png b/figs/APA95ML_nodes_1.png new file mode 100644 index 0000000..f7578b6 Binary files /dev/null and b/figs/APA95ML_nodes_1.png differ diff --git a/figs/APA95ML_nodes_2.png b/figs/APA95ML_nodes_2.png new file mode 100644 index 0000000..acc938e Binary files /dev/null and b/figs/APA95ML_nodes_2.png differ diff --git a/figs/compare_Gd_id_simscape.pdf b/figs/compare_Gd_id_simscape.pdf new file mode 100644 index 0000000..ff9c441 Binary files /dev/null and b/figs/compare_Gd_id_simscape.pdf differ diff --git a/figs/compare_Gd_id_simscape.png b/figs/compare_Gd_id_simscape.png new file mode 100644 index 0000000..e65bc5d Binary files /dev/null and b/figs/compare_Gd_id_simscape.png differ diff --git a/figs/compare_Gf_id_simscape.pdf b/figs/compare_Gf_id_simscape.pdf new file mode 100644 index 0000000..acb67bc Binary files /dev/null and b/figs/compare_Gf_id_simscape.pdf differ diff --git a/figs/compare_Gf_id_simscape.png b/figs/compare_Gf_id_simscape.png new file mode 100644 index 0000000..5863767 Binary files /dev/null and b/figs/compare_Gf_id_simscape.png differ diff --git a/figs/dynamics_act_disp_comp_mass.pdf b/figs/dynamics_act_disp_comp_mass.pdf index f0222c3..cfad4c3 100644 Binary files a/figs/dynamics_act_disp_comp_mass.pdf and b/figs/dynamics_act_disp_comp_mass.pdf differ diff --git a/figs/dynamics_act_disp_comp_mass.png b/figs/dynamics_act_disp_comp_mass.png index 6639efa..94202da 100644 Binary files a/figs/dynamics_act_disp_comp_mass.png and b/figs/dynamics_act_disp_comp_mass.png differ diff --git a/figs/gain_Va_to_Vs.pdf b/figs/gain_Va_to_Vs.pdf new file mode 100644 index 0000000..22f2266 Binary files /dev/null and b/figs/gain_Va_to_Vs.pdf differ diff --git a/figs/gain_Va_to_Vs.png b/figs/gain_Va_to_Vs.png new file mode 100644 index 0000000..c6569f2 Binary files /dev/null and b/figs/gain_Va_to_Vs.png differ diff --git a/figs/gain_Va_to_d.pdf b/figs/gain_Va_to_d.pdf new file mode 100644 index 0000000..ca30b2f Binary files /dev/null and b/figs/gain_Va_to_d.pdf differ diff --git a/figs/gain_Va_to_d.png b/figs/gain_Va_to_d.png new file mode 100644 index 0000000..f4a6a58 Binary files /dev/null and b/figs/gain_Va_to_d.png differ diff --git a/index.html b/index.html index 62d027b..17de956 100644 --- a/index.html +++ b/index.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Test Bench - Amplified Piezoelectric Actuator @@ -30,46 +30,52 @@

Table of Contents

-
-

1 Experimental Setup

+
+

1 Experimental Setup

- +

-A schematic of the test-bench is shown in Figure 1. +A schematic of the test-bench is shown in Figure 1.

@@ -105,31 +111,31 @@ The APA95ML has three stacks that can be used as actuator or as sensors.

-Pictures of the test bench are shown in Figure 2 and 3. +Pictures of the test bench are shown in Figure 2 and 3.

-
+

test_bench_apa_schematic.png

Figure 1: Schematic of the Setup

-
+

setup_picture.png

Figure 2: Picture of the Setup

-
+

setup_zoom.png

Figure 3: Zoom on the APA

-
+

Here are the equipment used in the test bench:

@@ -145,11 +151,504 @@ Here are the equipment used in the test bench:
-
-

2 Simscape model of the test-bench

+
+

2 Estimation of electrical/mechanical relationships

- +In order to correctly model the piezoelectric actuator, we need to determine: +

+
    +
  1. \(g_a\): the ratio of the generated force \(F_a\) to the supply voltage \(V_a\) across the piezoelectric stack
  2. +
  3. \(g_s\): the ratio of the generated voltage \(V_s\) across the piezoelectric stack when subject to a strain \(\Delta h\)
  4. +
+ +

+We estimate \(g_a\) and \(g_s\) using different approaches: +

+
    +
  1. Section 2.1: \(g_a\) is estimated from the datasheet of the piezoelectric stack
  2. +
  3. Section 2.2: \(g_a\) and \(g_s\) are estimated using the piezoelectric constants
  4. +
  5. Section 2.3: \(g_a\) and \(g_s\) are estimated experimentally
  6. +
+
+ +
+

2.1 Estimation from Data-sheet

+
+

+ +

+ +

+The stack parameters taken from the data-sheet are shown in Table 1. +

+ + + + +++ ++ ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Table 1: Stack Parameters
ParameterUnitValue
Nominal Stroke\(\mu m\)20
Blocked force\(N\)4700
Stiffness\(N/\mu m\)235
Voltage Range\(V\)-20..150
Capacitance\(\mu F\)4.4
Length\(mm\)20
Stack Area\(mm^2\)10x10
+ +

+Let’s compute the generated force +

+ +

+The stroke is \(L_{\max} = 20\mu m\) for a voltage range of \(V_{\max} = 170 V\). +Furthermore, the stiffness is \(k_a = 235 \cdot 10^6 N/m\). +

+ +

+The relation between the applied voltage and the generated force can be estimated as follows: +

+\begin{equation} + g_a = k_a \frac{L_{\max}}{V_{\max}} +\end{equation} + +
+
ka = 235e6; % [N/m]
+Lmax = 20e-6; % [m]
+Vmax = 170; % [V]
+
+
+ +
+
ka*Lmax/Vmax % [N/V]
+
+
+ +
+27.647
+
+ + +

+From the parameters of the stack, it seems not possible to estimate the relation between the strain and the generated voltage. +

+
+
+ +
+

2.2 Estimation from Piezoelectric parameters

+
+

+ +

+ +

+In order to make the conversion from applied voltage to generated force or from the strain to the generated voltage, we need to defined some parameters corresponding to the piezoelectric material: +

+
+
d33 = 300e-12; % Strain constant [m/V]
+n   = 80;      % Number of layers per stack
+ka  = 235e6;   % Stack stiffness [N/m]
+
+
+ +

+The ratio of the developed force to applied voltage is: +

+\begin{equation} +\label{org8300de8} + F_a = g_a V_a, \quad g_a = d_{33} n k_a +\end{equation} +

+where: +

+
    +
  • \(F_a\): developed force in [N]
  • +
  • \(n\): number of layers of the actuator stack
  • +
  • \(d_{33}\): strain constant in [m/V]
  • +
  • \(k_a\): actuator stack stiffness in [N/m]
  • +
  • \(V_a\): applied voltage in [V]
  • +
+ +

+If we take the numerical values, we obtain: +

+
+
ga = d33*n*ka; % [N/V]
+
+
+ +
+ga = 5.6 [N/V]
+
+ + + +

+From (Fleming and Leang 2014) (page 123), the relation between relative displacement of the sensor stack and generated voltage is: +

+\begin{equation} +\label{org312083f} + V_s = \frac{d_{33}}{\epsilon^T s^D n} \Delta h +\end{equation} +

+where: +

+
    +
  • \(V_s\): measured voltage in [V]
  • +
  • \(d_{33}\): strain constant in [m/V]
  • +
  • \(\epsilon^T\): permittivity under constant stress in [F/m]
  • +
  • \(s^D\): elastic compliance under constant electric displacement in [m^2/N]
  • +
  • \(n\): number of layers of the sensor stack
  • +
  • \(\Delta h\): relative displacement in [m]
  • +
+ +

+If we take the numerical values, we obtain: +

+
+
d33 = 300e-12; % Strain constant [m/V]
+n   = 80;      % Number of layers per stack
+eT  = 5.3e-9;  % Permittivity under constant stress [F/m]
+sD  = 2e-11;   % Compliance under constant electric displacement [m2/N]
+
+
+ +
+
gs = d33/(eT*sD*n); % [V/m]
+
+
+ +
+gs = 35.4 [V/um]
+
+
+
+ +
+

2.3 Estimation from Experiment

+
+

+ +

+

+The idea here is to obtain the parameters \(g_a\) and \(g_s\) from the comparison of an experimental identification and the identification using Simscape. +

+ +

+Using the experimental identification, we can easily obtain the gain from the applied voltage to the generated displacement, but not to the generated force. +However, from the Simscape model, we can easily have the link from the generated force to the displacement, them we can computed \(g_a\). +

+ +

+Similarly, it is fairly easy to experimentally obtain the gain from the stack displacement to the generated voltage across the stack. +To link that to the strain of the sensor stack, the simscape model is used. +

+
+
+

2.3.1 From actuator voltage \(V_a\) to actuator force \(F_a\)

+
+

+The data from the identification test is loaded. +

+
+
load('apa95ml_5kg_Amp_E505.mat', 't', 'um', 'y');
+
+% Any offset value is removed:
+um = detrend(um, 0); % Amplifier Input Voltage [V]
+y  = detrend(y , 0); % Mass displacement [m]
+
+
+ +

+Now we add a factor 10 to take into account the gain of the voltage amplifier and thus obtain the voltage across the piezoelectric stack. +

+
+
um = 10*um; % Stack Actuator Input Voltage [V]
+
+
+ +

+Then, the transfer function from the stack voltage to the vertical displacement is computed. +

+
+
Ts = t(end)/(length(t)-1);
+Fs = 1/Ts;
+
+win = hanning(ceil(1*Fs));
+
+[tf_est, f] = tfestimate(um, y, win, [], [], 1/Ts);
+
+
+ +

+The gain from input voltage of the stack to the vertical displacement is determined: +

+
+
g_d_Va = 4e-7; % [m/V]
+
+
+ + +
+

gain_Va_to_d.png +

+

Figure 4: Transfer function from actuator stack voltage \(V_a\) to vertical displacement of the mass \(d\)

+
+ +

+Then, the transfer function from forces applied by the stack actuator to the vertical displacement of the mass is identified from the Simscape model. +

+
+
m = 5.5;
+
+%% Name of the Simulink File
+mdl = 'piezo_amplified_3d';
+
+%% Input/Output definition
+clear io; io_i = 1;
+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]
+
+Gd = linearize(mdl, io);
+
+
+ +

+The DC gain the the identified dynamics +

+
+
g_d_Fa = abs(dcgain(Gd)); % [m/N]
+
+
+ +
+G_d_Fa = 1.2e-08 [m/N]
+
+ + + +

+And finally, the gain \(g_a\) from the the actuator voltage \(V_a\) to the generated force \(F_a\) can be computed: +

+\begin{equation} + g_a = \frac{F_a}{V_a} = \frac{F_a}{d} \frac{d}{V_a} +\end{equation} + +
+
ga = g_d_Va/g_d_Fa;
+
+
+ +
+ga = 33.7 [N/V]
+
+ + +

+The obtained comparison between the Simscape model and the identified dynamics is shown in Figure 5. +

+ + +
+

compare_Gd_id_simscape.png +

+

Figure 5: Comparison of the identified transfer function between actuator voltage \(V_a\) and vertical mass displacement \(d\)

+
+
+
+ +
+

2.3.2 From stack strain \(\Delta h\) to generated voltage \(V_s\)

+
+

+Now, the gain from the stack strain \(\Delta h\) to the generated voltage \(V_s\) is estimated. +

+ +

+We can determine the gain from actuator voltage \(V_a\) to sensor voltage \(V_s\) thanks to the identification. +Using the simscape model, we can have the transfer function from the actuator voltage \(V_a\) (using the previously estimated gain \(g_a\)) to the sensor stack strain \(\Delta h\). +

+ +

+Finally, using these two values, we can compute the gain \(g_s\) from the stack strain \(\Delta h\) to the generated Voltage \(V_s\). +

+ +

+Identification data is loaded. +

+
+
load('apa95ml_5kg_2a_1s.mat', 't', 'u', 'v');
+
+u = detrend(u, 0); % Input Voltage of the Amplifier [V]
+v = detrend(v, 0); % Voltage accross the stack sensor [V]
+
+
+ +

+Here, an amplifier with a gain of 20 is used. +

+
+
u = 20*u; % Input Voltage of the Amplifier [V]
+
+
+ +

+Then, the transfer function from \(V_a\) to \(V_s\) is identified and its DC gain is estimated (Figure 6). +

+
+
Ts = t(end)/(length(t)-1);
+Fs = 1/Ts;
+
+win = hann(ceil(10/Ts));
+
+[tf_est, f] = tfestimate(u, v, win, [], [], 1/Ts);
+
+
+ +
+
g_Vs_Va = 0.022; % [V/V]
+
+
+ + +
+

gain_Va_to_Vs.png +

+

Figure 6: Transfer function from actuator stack voltage \(V_a\) to sensor stack voltage \(V_s\)

+
+ + +

+Now the transfer function from the actuator stack voltage to the sensor stack strain is estimated using the Simscape model. +

+
+
m = 5.5;
+
+%% Name of the Simulink File
+mdl = 'piezo_amplified_3d';
+
+%% Input/Output definition
+clear io; io_i = 1;
+io(io_i) = linio([mdl, '/Va'], 1, 'openinput');  io_i = io_i + 1; % Actuator Voltage [V]
+io(io_i) = linio([mdl, '/dL'], 1, 'openoutput'); io_i = io_i + 1; % Sensor Stack displacement [m]
+
+Gf = linearize(mdl, io);
+
+
+ +

+The gain from the actuator stack voltage to the sensor stack strain is estimated below. +

+
+
G_dh_Va = abs(dcgain(Gf));
+
+
+ +
+G_dh_Va = 6.2e-09 [m/V]
+
+ + +

+And finally, the gain \(g_s\) from the sensor stack strain to the generated voltage can be estimated: +

+\begin{equation} + g_s = \frac{V_s}{\Delta h} = \frac{V_s}{V_a} \frac{V_a}{\Delta h} +\end{equation} + +
+
gs = g_Vs_Va/G_dh_Va; % [V/m]
+
+
+ +
+gs = 3.5 [V/um]
+
+ + + +
+

compare_Gf_id_simscape.png +

+

Figure 7: Comparison of the identified transfer function between actuator voltage \(V_a\) and sensor stack voltage

+
+
+
+
+ +
+

2.4 Conclusion

+
+

+The obtained parameters \(g_a\) and \(g_s\) are not consistent between the different methods. +

+ +

+The one using the experimental data are saved and further used. +

+
+
save('./matlab/mat/apa95ml_params.mat', 'ga', 'gs');
+
+
+
+
+
+ +
+

3 Simscape model of the test-bench

+
+

+

The idea here is to model the test-bench using Simscape. @@ -160,23 +659,31 @@ Whereas the suspended mass and metrology frame can be considered as rigid bodies

-To model the APA, a Finite Element Model (FEM) is used and imported into Simscape. +To model the APA, a Finite Element Model (FEM) is used (Figure 8) and imported into Simscape.

+ + +
+

APA95ML_FEM.png +

+

Figure 8: Finite Element Model of the APA95ML

-
-

2.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.

-
K = extractMatrix('APA95ML_K.txt');
-M = extractMatrix('APA95ML_M.txt');
+
K = readmatrix('APA95ML_K.CSV');
+M = readmatrix('APA95ML_M.CSV');
 
- +@@ -202,131 +709,131 @@ M = extractMatrix('APA95ML_M.txt'); + - - - - + + + - - - + + + - - - - - - - - - - - - - - - + - - - + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - + + + + - - - - - - - - - - - - - - - - - + + - - + - - - - + + + + - + + - + + + - + + + + + + + + + + + + - - - + + - + + @@ -334,7 +841,7 @@ M = extractMatrix('APA95ML_M.txt');
Table 1: First 10x10 elements of the Stiffness matrixTable 2: First 10x10 elements of the Stiffness matrix
300000000.0-1000.0 -30000.08000.0-200.0-30.0-60000.0-40.070000.0300.0 20000000.0-4000.0500.08-30.0-5000.05
-30000.0100000000.0400.030.0200.0-14000.0-8000000.0800.07
8000.0400.0-1000.0 50000000.0-800000.0-300.0-40.0-7000.0800000.0-20.0 300.0100.05000000.040000.0
-200.030.0-800000.020000.051-10.0-2-40000.0-300.0
-30.0200.0-300.0540000.00.3-4-10.040.00.4
-60000.0-1-40.010.3 3000.07000.00.8-10.00035000000.0400.0-40000.0
-30000.0-7000.0100000000.0-200.0-60.070.03000.03000.0-8000000.0-30.0
-40.0800000.0-200.020000.0-0.4430.040000.07-300.0
70000.0-20.0-60.0-0.43000.01-6000.010.08-0.1
300.0300.070.04140000.0-10.0-10.030.00.1
20000000.04000.0300.03000.03000.030.0-6000.0 -10.0-47000.0 300000000.020000.03000.080.0
-4000.0-8000000.0100.0-2-10.00.820000.0100000000.0-4000.02000.09000.0 -100.0
500.0800.0-30.0 5000000.0-40000.040.0-1 3000.0-4000.040000.010.0-10.02000.0 50000000.0800000.0-3000.0-800000.0
8-5000.0400.0-8000000.0 740000.0830.09000.0-3000.0100000000.0100.0
5-40000.0-30.0 -300.00.40.000380.0-0.10.1 -100.0800000.0-800000.0100.0 20000.0
- +@@ -360,131 +867,131 @@ M = extractMatrix('APA95ML_M.txt'); + - - - - + + + - - - + + + - + - - - - - - - - - - - - - - - - - - - + + + + + - - - - - - - - - - - - + - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - + + + + + + - - + - - - - - - - - - - - - - - - + - + + - - - + + - + + - - - + + + + + + + + + + + + + + + + - - - - - + + + + + @@ -499,6 +1006,10 @@ Then, we extract the coordinates of the interface nodes. +

+The interface nodes are shown in Figure 9 and their coordinates are listed in Table 4. +

+
Table 2: First 10x10 elements of the Mass matrixTable 3: First 10x10 elements of the Mass matrix
0.037e-08 2e-06-2e-071e-082e-080.0002-3e-09-0.0002-6e-08 -0.0012e-07-8e-08-9e-108e-076e-07-8e-09
2e-067e-08 0.02-5e-077e-093e-082e-08-3e-070.0003-1e-081e-10
-2e-07-5e-070.02-9e-054e-09-1e-082e-07-2e-08-1e-069e-05-3e-09-4e-09-1e-06 -0.0006-5e-06
1e-087e-09-9e-051e-066e-114e-10-1e-093e-11-4e-08 5e-063e-08
2e-083e-084e-096e-111e-062e-10-2e-092e-10-7e-09-4e-11
0.00022e-08-1e-084e-102e-10 2e-06-2e-06-1e-09-7e-10-9e-12-1e-060.02-3e-08-4e-081e-081e-07-2e-070.00031e-09
-3e-099e-05-3e-081e-06-3e-11-3e-13-7e-09-5e-06-3e-103e-08
-0.0002-3e-09-4e-08-3e-112e-066e-102e-06-7e-09-2e-097e-11
-6e-08-4e-091e-08-3e-136e-101e-061e-083e-09-2e-092e-13
-0.001-3e-072e-07-1e-09-2e-09-2e-06-1e-061e-07-7e-092e-061e-08 0.034e-08 -2e-06-1e-07-5e-098e-09
2e-070.0003-2e-083e-112e-10-1e-09-2e-060.02-8e-07-1e-08
-8e-08-1e-088e-07 -0.00065e-06-2e-07-5e-06 -7e-09-7e-10-1e-07-8e-073e-094e-08 0.029e-05-9e-07-9e-05
-9e-101e-10-5e-066e-07-4e-080.0003-3e-10-2e-09-2e-09-2e-06-9e-070.022e-08
-8e-095e-061e-09 3e-08-4e-11-9e-12-5e-09-1e-089e-057e-112e-138e-09-9e-052e-08 1e-06
@@ -510,28 +1021,28 @@ Then, we extract the coordinates of the interface nodes. - + - + - + - +
Total number of Nodes1689597
Number of interface Nodes137
Number of Modes306
Size of M and K matrices10848
- - +
Table 3: Coordinates of the interface nodes
+@@ -556,23 +1067,23 @@ Then, we extract the coordinates of the interface nodes. - + - + - + - + - + @@ -580,175 +1091,54 @@ Then, we extract the coordinates of the interface nodes. - - + + - - + + - - + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + - - - - - - - -
Table 4: Coordinates of the interface nodes
1.0168947.040467.0 0.00.03 0.00.029997
2.0168949.040469.0 0.0-0.03 0.0-0.029997
3.0168950.040470.0 -0.035 0.0 0.0
4.0168951.0-0.02840475.0-0.015 0.0 0.0
5.0168952.0-0.02140476.0-0.005 0.0 0.0
6.0168953.0-0.01440477.00.015 0.0 0.0
7.0168954.0-0.0070.00.0
8.0168955.00.00.00.0
9.0168956.00.0070.00.0
10.0168957.00.0140.00.0
11.0168958.00.0210.00.0
12.0168959.040478.0 0.035 0.0 0.0
13.0168960.00.0280.00.0
+ +
+

APA95ML_nodes_1.png +

+

Figure 9: Interface Nodes chosen for the APA95ML

+
+

Using K, M and int_xyz, we can use the Reduced Order Flexible Solid simscape block.

-
-

2.2 Piezoelectric parameters

-
-

-In order to make the conversion from applied voltage to generated force or from the strain to the generated voltage, we need to defined some parameters corresponding to the piezoelectric material: -

-
-
d33 = 300e-12; % Strain constant [m/V]
-n   = 80;      % Number of layers per stack
-eT  = 1.6e-8;  % Permittivity under constant stress [F/m]
-sD  = 1e-11;   % Compliance under constant electric displacement [m2/N]
-ka  = 235e6;   % Stack stiffness [N/m]
-C   = 5e-6;    % Stack capactiance [F]
-
-
- -

-The ratio of the developed force to applied voltage is: -

-\begin{equation} -\label{org1eea105} - F_a = g_a V_a, \quad g_a = d_{33} n k_a -\end{equation} -

-where: -

-
    -
  • \(F_a\): developed force in [N]
  • -
  • \(n\): number of layers of the actuator stack
  • -
  • \(d_{33}\): strain constant in [m/V]
  • -
  • \(k_a\): actuator stack stiffness in [N/m]
  • -
  • \(V_a\): applied voltage in [V]
  • -
- -

-If we take the numerical values, we obtain: -

-
-
d33*n*ka % [N/V]
-
-
- -
-5.64
-
- - -

-From (Fleming and Leang 2014) (page 123), the relation between relative displacement of the sensor stack and generated voltage is: -

-\begin{equation} -\label{org7e93246} - V_s = \frac{d_{33}}{\epsilon^T s^D n} \Delta h -\end{equation} -

-where: -

-
    -
  • \(V_s\): measured voltage in [V]
  • -
  • \(d_{33}\): strain constant in [m/V]
  • -
  • \(\epsilon^T\): permittivity under constant stress in [F/m]
  • -
  • \(s^D\): elastic compliance under constant electric displacement in [m^2/N]
  • -
  • \(n\): number of layers of the sensor stack
  • -
  • \(\Delta h\): relative displacement in [m]
  • -
- -

-If we take the numerical values, we obtain: -

-
-
1e-6*d33/(eT*sD*n) % [V/um]
-
-
- -
-23.438
-
-
-
- -
-

2.3 Simscape Model

-
+
+

3.2 Simscape Model

+

The flexible element is imported using the Reduced Order Flexible Solid Simscape block.

@@ -758,10 +1148,6 @@ To model the actuator, an Internal Force block is added between the A Relative Motion Sensor block is added between the nodes 1 and 2 to measure the displacement and the amplified piezo.

-
    -
  • [ ] Add schematic of the model with interface nodes
  • -
-

One mass is fixed at one end of the piezo-electric stack actuator, the other end is fixed to the world frame.

@@ -772,11 +1158,11 @@ One mass is fixed at one end of the piezo-electric stack actuator, the other end
-
-

2.4 Dynamics from Actuator Voltage to Vertical Mass Displacement

-
+
+

3.3 Dynamics from Actuator Voltage to Vertical Mass Displacement

+

-The identified dynamics is shown in Figure 4. +The identified dynamics is shown in Figure 10.

@@ -793,19 +1179,19 @@ Ghm = -linearize(mdl, io);
-
+

dynamics_act_disp_comp_mass.png

-

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

+

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

-
-

2.5 Dynamics from Actuator Voltage to Force Sensor Voltage

-
+
+

3.4 Dynamics from Actuator Voltage to Force Sensor Voltage

+

-The obtained dynamics is shown in Figure 5. +The obtained dynamics is shown in Figure 11.

@@ -822,17 +1208,17 @@ Gfm = linearize(mdl, io);
-
+

dynamics_force_force_sensor_comp_mass.png

-

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

+

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

-
-

2.6 Save Data for further use

-
+
+

3.5 Save Data for further use

+
save('matlab/mat/fem_simscape_models.mat', 'Ghm', 'Gfm')
 
@@ -846,250 +1232,27 @@ Gfm = linearize(mdl, io);
-
-

3 Estimation of piezoelectric parameters

-
-
-
-

3.1 From actuator voltage to vertical displacement

-
-

-The data from the “noise test” and the identification test are loaded. -

-
-
load('apa95ml_5kg_Amp_E505.mat', 't', 'um', 'y');
-
-
- -

-Any offset value is removed: -

-
-
um = detrend(um, 0); % Amplifier Input Voltage [V]
-y  = detrend(y , 0); % Mass displacement [m]
-
-
- -

-Now we add a factor 10 to take into account the gain of the voltage amplifier. -

-
-
um = 10*um; % Stack Actuator Input Voltage [V]
-
-
- -
-
Ts = t(end)/(length(t)-1);
-Fs = 1/Ts;
-
-win = hanning(ceil(1*Fs));
-
-
- -
-
[tf_est, f] = tfestimate(um, y, win, [], [], 1/Ts);
-
-
- -

-The gain from input voltage of the stack to the vertical displacement is determined: -

-
-
gD = 4e-7; % [m/V]
-
-
- -
-
K = extractMatrix('APA95ML_K.txt');
-M = extractMatrix('APA95ML_M.txt');
-[int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('APA95ML_out_nodes_3D.txt');
-
-
- -

-Define parameters just for the simulation. Should not change any results -

-
-
d33 = 1; % Strain constant [m/V]
-n   = 1; % Number of layers per stack
-eT  = 1; % Permittivity under constant stress [F/m]
-sD  = 1; % Compliance under constant electric displacement [m2/N]
-ka  = 1; % Stack stiffness [N/m]
-C   = 1; % Stack capactiance [F]
-
-
- -
-
m = 5;
-
-%% Name of the Simulink File
-mdl = 'piezo_amplified_3d';
-
-%% Input/Output definition
-clear io; io_i = 1;
-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]
-
-Gd = linearize(mdl, io);
-
-
- -
-
gF = abs(dcgain(Gd)); % [m/N]
-ans = gF
-
-
- -
-6.1695e-09
-
- - -

-\[ g_a = g_D/g_F \] -in [N/V] -

-
-
ga = gD/gF
-ans = ga
-
-
- -
-64.835
-
- - -
-
na = 2; ns = 1; d33 = 300e-12; n = 80; ka = 235e6; gL = 1.7;
-na/(na+ns)*gL*d33*n*ka
-
-
- -
-6.392
-
- - -
    -
  • [ ] Why is there a factor 10 with the “theoretical estimation?”
  • -
-
-
- -
-

3.2 From actuator voltage to sensor Voltage

-
-
-
load('apa95ml_5kg_2a_1s.mat', 't', 'u', 'v');
-
-
- -
-
u = detrend(u, 0); % Input Voltage of the Amplifier [V]
-v = detrend(v, 0); % Voltage accross the stack sensor [V]
-
-
- -
-
u = 20*u; % Input Voltage of the Amplifier [V]
-
-
- -
-
Ts = t(end)/(length(t)-1);
-Fs = 1/Ts;
-
-win = hann(ceil(10/Ts));
-
-[tf_est, f] = tfestimate(u, v, win, [], [], 1/Ts);
-
-
- -
-
gV = 0.022; % [V/V]
-
-
- -
-
m = 5;
-
-%% Name of the Simulink File
-mdl = 'piezo_amplified_3d';
-
-%% Input/Output definition
-clear io; io_i = 1;
-io(io_i) = linio([mdl, '/Fa'], 1, 'openinput');  io_i = io_i + 1; % Actuator Force [N]
-io(io_i) = linio([mdl, '/dL'], 1, 'openoutput'); io_i = io_i + 1; % Sensor Stack displacement [m]
-
-Gf = linearize(mdl, io);
-
-
- -

-\(g_F\) in [m/N] -

-
-
gF = abs(dcgain(Gf));
-ans = gF
-
-
- -
-2.1546e-10
-
- - -

-Finally, we compute the gain from strain of the sensor stack to its generated voltage: -\[ g_S = \frac{g_V}{g_F g_a} \] -Gs in [V/m] -

-
-
gS = gV/gF/ga;
-ans = gS
-
-
- -
-15749000.0
-
- - -
-
d33 = 300e-12; eT = 5.3e-9; sD = 2e-11; n = 80;
-d33/(eT*sD*n)
-
-
- -
-35377000.0
-
-
-
-
- -
-

4 Huddle Test

+
+

4 Huddle Test

- +

-
-

4.1 Time Domain Data

+
+

4.1 Time Domain Data

-
+

huddle_test_time_domain.png

-

Figure 6: Measurement of the Mass displacement during Huddle Test

+

Figure 12: Measurement of the Mass displacement during Huddle Test

-
-

4.2 PSD of Measurement Noise

+
+

4.2 PSD of Measurement Noise

Ts = t(end)/(length(t)-1);
@@ -1105,20 +1268,20 @@ win = hanning(ceil(1*Fs));
 
-
+

huddle_test_pdf.png

-

Figure 7: Amplitude Spectral Density of the Displacement during Huddle Test

+

Figure 13: Amplitude Spectral Density of the Displacement during Huddle Test

-
-

5 Identification of the dynamics from actuator to displacement

+
+

5 Identification of the dynamics from actuator to displacement

- +

  • [ ] List of equipment
  • @@ -1130,8 +1293,8 @@ win = hanning(ceil(1*Fs)); E505 with gain of 10.

-
-

5.1 Load Data

+
+

5.1 Load Data

The data from the “noise test” and the identification test are loaded. @@ -1165,8 +1328,8 @@ ht.u = 10*ht.u;

-
-

5.2 Comparison of the PSD with Huddle Test

+
+

5.2 Comparison of the PSD with Huddle Test

Ts = t(end)/(length(t)-1);
@@ -1183,16 +1346,16 @@ win = hanning(ceil(1*Fs));
 
-
+

apa95ml_5kg_PI_pdf_comp_huddle.png

-

Figure 8: Comparison of the ASD for the identification test and the huddle test

+

Figure 14: Comparison of the ASD for the identification test and the huddle test

-
-

5.3 Compute TF estimate and Coherence

+
+

5.3 Compute TF estimate and Coherence

[tf_est, f] = tfestimate(um, -y, win, [], [], 1/Ts);
@@ -1201,10 +1364,10 @@ win = hanning(ceil(1*Fs));
 
-
+

apa95ml_5kg_PI_coh.png

-

Figure 9: Coherence

+

Figure 15: Coherence

@@ -1216,20 +1379,20 @@ Comparison with the FEM model

-
+

apa95ml_5kg_pi_comp_fem.png

-

Figure 10: Comparison of the identified transfer function and the one estimated from the FE model

+

Figure 16: Comparison of the identified transfer function and the one estimated from the FE model

-
-

6 Identification of the dynamics from actuator to force sensor

+
+

6 Identification of the dynamics from actuator to force sensor

- +

Two measurements are performed: @@ -1291,7 +1454,7 @@ Gfem_a_ss = exp(-s/1

-The transfer function from input voltage to output voltage are computed and shown in Figure 11. +The transfer function from input voltage to output voltage are computed and shown in Figure 17.

@@ -1311,14 +1474,14 @@ win = hann(ceil(10/Ts));
-
+

bode_plot_force_sensor_voltage_comp_fem.png

-

Figure 11: Comparison of the identified dynamics from voltage output to voltage input and the FEM

+

Figure 17: Comparison of the identified dynamics from voltage output to voltage input and the FEM

-
-

6.1 System Identification

+
+

6.1 System Identification

w_z = 2*pi*111; % Zeros frequency [rad/s]
@@ -1332,44 +1495,44 @@ Gi = G_inf*(s^2 
 
 
-
+

iff_plant_identification_apa95ml.png

-

Figure 12: Identification of the IFF plant

+

Figure 18: Identification of the IFF plant

-
-

6.2 Integral Force Feedback

+
+

6.2 Integral Force Feedback

-
+

root_locus_iff_apa95ml_identification.png

-

Figure 13: Root Locus for IFF

+

Figure 19: Root Locus for IFF

-
-

7 Integral Force Feedback

+
+

7 Integral Force Feedback

- +

-
+

test_bench_apa_schematic_iff.png

-

Figure 14: Schematic of the test bench using IFF

+

Figure 20: Schematic of the test bench using IFF

-
-

7.1 First tests with few gains

+
+

7.1 First tests with few gains

iff_g10  = load('apa95ml_iff_g10_res.mat',  'u', 't', 'y', 'v');
@@ -1394,24 +1557,24 @@ win = hann(ceil(10/Ts));
 
-
+

iff_first_test_coherence.png

-

Figure 15: Coherence

+

Figure 21: Coherence

-
+

iff_first_test_bode_plot.png

-

Figure 16: Bode plot for different values of IFF gain

+

Figure 22: Bode plot for different values of IFF gain

-
-

7.2 Second test with many Gains

+
+

7.2 Second test with many Gains

load('apa95ml_iff_test.mat', 'results');
@@ -1440,7 +1603,7 @@ g_iff = [0, 1, 5, 10, 50, 100];
 
-
+

iff_results_bode_plots.png

@@ -1462,13 +1625,13 @@ f_end = 500; % [Hz]
-
+

iff_results_bode_plots_identification.png

-
+

iff_results_root_locus.png

@@ -1483,7 +1646,7 @@ f_end = 500; % [Hz]

Author: Dehaeze Thomas

-

Created: 2020-11-24 mar. 09:17

+

Created: 2020-11-24 mar. 13:54

diff --git a/index.org b/index.org index 1c3670e..27dca3c 100644 --- a/index.org +++ b/index.org @@ -67,18 +67,442 @@ Here are the equipment used in the test bench: - Low Noise Voltage Preamplifier from Ametek ([[file:doc/model_5113.pdf][doc]]) #+end_note +* Estimation of electrical/mechanical relationships +** Introduction :ignore: +In order to correctly model the piezoelectric actuator with Simscape, we need to determine: +1. $g_a$: the ratio of the generated force $F_a$ to the supply voltage $V_a$ across the piezoelectric stack +2. $g_s$: the ratio of the generated voltage $V_s$ across the piezoelectric stack when subject to a strain $\Delta h$ + +We estimate $g_a$ and $g_s$ using different approaches: +1. Section [[sec:estimation_datasheet]]: $g_a$ is estimated from the datasheet of the piezoelectric stack +2. Section [[sec:estimation_piezo_params]]: $g_a$ and $g_s$ are estimated using the piezoelectric constants +3. Section [[sec:estimation_identification]]: $g_a$ and $g_s$ are estimated experimentally + +** Matlab Init :noexport:ignore: +#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) + <> +#+end_src + +#+begin_src matlab :exports none :results silent :noweb yes + <> +#+end_src + +#+begin_src matlab :tangle no + addpath('./matlab/mat/'); + addpath('./matlab/'); +#+end_src + +#+begin_src matlab :eval no + addpath('./mat/'); +#+end_src + +** Estimation from Data-sheet +<> + +The stack parameters taken from the data-sheet are shown in Table [[tab:stack_parameters]]. + +#+name: tab:stack_parameters +#+caption: Stack Parameters +| Parameter | Unit | Value | +|----------------+-----------+----------| +| Nominal Stroke | $\mu m$ | 20 | +| Blocked force | $N$ | 4700 | +| Stiffness | $N/\mu m$ | 235 | +| Voltage Range | $V$ | -20..150 | +| Capacitance | $\mu F$ | 4.4 | +| Length | $mm$ | 20 | +| Stack Area | $mm^2$ | 10x10 | + +Let's compute the generated force + +The stroke is $L_{\max} = 20\mu m$ for a voltage range of $V_{\max} = 170 V$. +Furthermore, the stiffness is $k_a = 235 \cdot 10^6 N/m$. + +The relation between the applied voltage and the generated force can be estimated as follows: +\begin{equation} + g_a = k_a \frac{L_{\max}}{V_{\max}} +\end{equation} + +#+begin_src matlab + ka = 235e6; % [N/m] + Lmax = 20e-6; % [m] + Vmax = 170; % [V] +#+end_src + +#+begin_src matlab :results value replace + ka*Lmax/Vmax % [N/V] +#+end_src + +#+RESULTS: +: 27.647 + +From the parameters of the stack, it seems not possible to estimate the relation between the strain and the generated voltage. + +** Estimation from Piezoelectric parameters +<> + +In order to make the conversion from applied voltage to generated force or from the strain to the generated voltage, we need to defined some parameters corresponding to the piezoelectric material: +#+begin_src matlab + d33 = 300e-12; % Strain constant [m/V] + n = 80; % Number of layers per stack + ka = 235e6; % Stack stiffness [N/m] +#+end_src + +The ratio of the developed force to applied voltage is: +#+name: eq:piezo_voltage_to_force +\begin{equation} + F_a = g_a V_a, \quad g_a = d_{33} n k_a +\end{equation} +where: +- $F_a$: developed force in [N] +- $n$: number of layers of the actuator stack +- $d_{33}$: strain constant in [m/V] +- $k_a$: actuator stack stiffness in [N/m] +- $V_a$: applied voltage in [V] + +If we take the numerical values, we obtain: +#+begin_src matlab + ga = d33*n*ka; % [N/V] +#+end_src + +#+begin_src matlab :results value replace :exports results :tangle no + sprintf('ga = %.1f [N/V]', ga) +#+end_src + +#+RESULTS: +: ga = 5.6 [N/V] + + +From cite:fleming14_desig_model_contr_nanop_system (page 123), the relation between relative displacement of the sensor stack and generated voltage is: +#+name: eq:piezo_strain_to_voltage +\begin{equation} + V_s = \frac{d_{33}}{\epsilon^T s^D n} \Delta h +\end{equation} +where: +- $V_s$: measured voltage in [V] +- $d_{33}$: strain constant in [m/V] +- $\epsilon^T$: permittivity under constant stress in [F/m] +- $s^D$: elastic compliance under constant electric displacement in [m^2/N] +- $n$: number of layers of the sensor stack +- $\Delta h$: relative displacement in [m] + +If we take the numerical values, we obtain: +#+begin_src matlab + d33 = 300e-12; % Strain constant [m/V] + n = 80; % Number of layers per stack + eT = 5.3e-9; % Permittivity under constant stress [F/m] + sD = 2e-11; % Compliance under constant electric displacement [m2/N] +#+end_src + +#+begin_src matlab + gs = d33/(eT*sD*n); % [V/m] +#+end_src + +#+begin_src matlab :results value replace :exports results :tangle no + sprintf('gs = %.1f [V/um]', 1e-6*gs) +#+end_src + +#+RESULTS: +: gs = 35.4 [V/um] + +** Estimation from Experiment +<> +*** Introduction :ignore: +The idea here is to obtain the parameters $g_a$ and $g_s$ from the comparison of an experimental identification and the identification using Simscape. + +Using the experimental identification, we can easily obtain the gain from the applied voltage to the generated displacement, but not to the generated force. +However, from the Simscape model, we can easily have the link from the generated force to the displacement, them we can computed $g_a$. + +Similarly, it is fairly easy to experimentally obtain the gain from the stack displacement to the generated voltage across the stack. +To link that to the strain of the sensor stack, the simscape model is used. + +*** From actuator voltage $V_a$ to actuator force $F_a$ +The data from the identification test is loaded. +#+begin_src matlab + load('apa95ml_5kg_Amp_E505.mat', 't', 'um', 'y'); + + % Any offset value is removed: + um = detrend(um, 0); % Amplifier Input Voltage [V] + y = detrend(y , 0); % Mass displacement [m] +#+end_src + +Now we add a factor 10 to take into account the gain of the voltage amplifier and thus obtain the voltage across the piezoelectric stack. +#+begin_src matlab + um = 10*um; % Stack Actuator Input Voltage [V] +#+end_src + +Then, the transfer function from the stack voltage to the vertical displacement is computed. +#+begin_src matlab + Ts = t(end)/(length(t)-1); + Fs = 1/Ts; + + win = hanning(ceil(1*Fs)); + + [tf_est, f] = tfestimate(um, y, win, [], [], 1/Ts); +#+end_src + +The gain from input voltage of the stack to the vertical displacement is determined: +#+begin_src matlab + g_d_Va = 4e-7; % [m/V] +#+end_src + +#+begin_src matlab :exports none + freqs = logspace(0, 4, 1000); + + figure; + hold on; + plot(f, abs(tf_est), 'DisplayName', 'Identification') + plot([f(2) f(end)], [g_d_Va g_d_Va], '--', 'DisplayName', sprintf('$d/V_a \\approx %.0g\\ m/V$', g_d_Va)) + hold off; + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Amplitude $d/V_a$ [m/V]'); + xlim([10, 5e3]); ylim([1e-8, 1e-5]); + legend('location', 'southwest'); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/gain_Va_to_d.pdf', 'width', 'wide', 'height', 'normal'); +#+end_src + +#+name: fig:gain_Va_to_d +#+caption: Transfer function from actuator stack voltage $V_a$ to vertical displacement of the mass $d$ +#+RESULTS: +[[file:figs/gain_Va_to_d.png]] + +Then, the transfer function from forces applied by the stack actuator to the vertical displacement of the mass is identified from the Simscape model. +#+begin_src matlab :exports none + % Load Parameters + K = readmatrix('APA95ML_K.CSV'); + M = readmatrix('APA95ML_M.CSV'); + [int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('APA95ML_out_nodes_3D.txt'); + + % Define parameters just for the simulation + ga = 1; % [N/V] + gs = 1; % [V/m] +#+end_src + +#+begin_src matlab + m = 5.5; + + %% Name of the Simulink File + mdl = 'piezo_amplified_3d'; + + %% Input/Output definition + clear io; io_i = 1; + 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] + + Gd = linearize(mdl, io); +#+end_src + +The DC gain the the identified dynamics +#+begin_src matlab + g_d_Fa = abs(dcgain(Gd)); % [m/N] +#+end_src + +#+begin_src matlab :results value replace :exports results :tangle no + sprintf('G_d_Fa = %.2g [m/N]', g_d_Fa) +#+end_src + +#+RESULTS: +: G_d_Fa = 1.2e-08 [m/N] + + +And finally, the gain $g_a$ from the the actuator voltage $V_a$ to the generated force $F_a$ can be computed: +\begin{equation} + g_a = \frac{F_a}{V_a} = \frac{F_a}{d} \frac{d}{V_a} +\end{equation} + +#+begin_src matlab + ga = g_d_Va/g_d_Fa; +#+end_src + +#+begin_src matlab :results value replace :exports results :tangle no + sprintf('ga = %.1f [N/V]', ga) +#+end_src + +#+RESULTS: +: ga = 33.7 [N/V] + +The obtained comparison between the Simscape model and the identified dynamics is shown in Figure [[fig:compare_Gd_id_simscape]]. + +#+begin_src matlab :exports none + freqs = logspace(1, 4, 1000); + + figure; + hold on; + plot(f, abs(tf_est), 'DisplayName', 'Identification') + plot(f, ga*abs(squeeze(freqresp(Gd, f, 'Hz'))), 'DisplayName', 'Simscape Model') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); + ylabel('Amplitude $d/V_a$ [m/V]'); xlabel('Frequency [Hz]'); + hold off; + xlim([1, 5e3]); ylim([1e-9, 1e-5]); + legend('location', 'southwest'); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/compare_Gd_id_simscape.pdf', 'width', 'wide', 'height', 'normal'); +#+end_src + +#+name: fig:compare_Gd_id_simscape +#+caption: Comparison of the identified transfer function between actuator voltage $V_a$ and vertical mass displacement $d$ +#+RESULTS: +[[file:figs/compare_Gd_id_simscape.png]] + +*** From stack strain $\Delta h$ to generated voltage $V_s$ +Now, the gain from the stack strain $\Delta h$ to the generated voltage $V_s$ is estimated. + +We can determine the gain from actuator voltage $V_a$ to sensor voltage $V_s$ thanks to the identification. +Using the simscape model, we can have the transfer function from the actuator voltage $V_a$ (using the previously estimated gain $g_a$) to the sensor stack strain $\Delta h$. + +Finally, using these two values, we can compute the gain $g_s$ from the stack strain $\Delta h$ to the generated Voltage $V_s$. + +Identification data is loaded. +#+begin_src matlab + load('apa95ml_5kg_2a_1s.mat', 't', 'u', 'v'); + + u = detrend(u, 0); % Input Voltage of the Amplifier [V] + v = detrend(v, 0); % Voltage accross the stack sensor [V] +#+end_src + +Here, an amplifier with a gain of 20 is used. +#+begin_src matlab + u = 20*u; % Input Voltage of the Amplifier [V] +#+end_src + +Then, the transfer function from $V_a$ to $V_s$ is identified and its DC gain is estimated (Figure [[fig:gain_Va_to_Vs]]). +#+begin_src matlab + Ts = t(end)/(length(t)-1); + Fs = 1/Ts; + + win = hann(ceil(10/Ts)); + + [tf_est, f] = tfestimate(u, v, win, [], [], 1/Ts); +#+end_src + +#+begin_src matlab + g_Vs_Va = 0.022; % [V/V] +#+end_src + +#+begin_src matlab :exports none + freqs = logspace(1, 4, 1000); + + figure; + hold on; + plot(f, abs(tf_est), 'DisplayName', 'Identification') + plot([f(2); f(end)], [g_Vs_Va; g_Vs_Va], '--', 'DisplayName', sprintf('$V_s/V_a \\approx %.0g\\ m/V$', g_Vs_Va)) + hold off; + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); + ylabel('Amplitude $V_s/V_a$ [V/V]'); xlabel('Frequency [Hz]'); + xlim([5e-1 5e3]); ylim([1e-3, 1e1]); + legend('location', 'northwest'); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/gain_Va_to_Vs.pdf', 'width', 'wide', 'height', 'normal'); +#+end_src + +#+name: fig:gain_Va_to_Vs +#+caption: Transfer function from actuator stack voltage $V_a$ to sensor stack voltage $V_s$ +#+RESULTS: +[[file:figs/gain_Va_to_Vs.png]] + + +Now the transfer function from the actuator stack voltage to the sensor stack strain is estimated using the Simscape model. +#+begin_src matlab + m = 5.5; + + %% Name of the Simulink File + mdl = 'piezo_amplified_3d'; + + %% Input/Output definition + clear io; io_i = 1; + io(io_i) = linio([mdl, '/Va'], 1, 'openinput'); io_i = io_i + 1; % Actuator Voltage [V] + io(io_i) = linio([mdl, '/dL'], 1, 'openoutput'); io_i = io_i + 1; % Sensor Stack displacement [m] + + Gf = linearize(mdl, io); +#+end_src + +The gain from the actuator stack voltage to the sensor stack strain is estimated below. +#+begin_src matlab + G_dh_Va = abs(dcgain(Gf)); +#+end_src + +#+begin_src matlab :results value replace :exports results :tangle no + sprintf('G_dh_Va = %.2g [m/V]', G_dh_Va); +#+end_src + +#+RESULTS: +: G_dh_Va = 6.2e-09 [m/V] + +And finally, the gain $g_s$ from the sensor stack strain to the generated voltage can be estimated: +\begin{equation} + g_s = \frac{V_s}{\Delta h} = \frac{V_s}{V_a} \frac{V_a}{\Delta h} +\end{equation} + +#+begin_src matlab + gs = g_Vs_Va/G_dh_Va; % [V/m] +#+end_src + +#+begin_src matlab :results value replace :exports results :tangle no + sprintf('gs = %.1f [V/um]', 1e-6*gs) +#+end_src + +#+RESULTS: +: gs = 3.5 [V/um] + +#+begin_src matlab :exports none + freqs = logspace(1, 4, 1000); + + figure; + hold on; + plot(f, abs(tf_est), 'DisplayName', 'Identification') + plot(f, gs*abs(squeeze(freqresp(Gf, f, 'Hz'))), 'DisplayName', 'Simscape Model') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); + ylabel('Amplitude $V_s/V_a$ [V/V]'); xlabel('Frequency [Hz]'); + hold off; + xlim([1, 5e3]); ylim([1e-3, 1e1]); + legend('location', 'southwest'); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/compare_Gf_id_simscape.pdf', 'width', 'wide', 'height', 'normal'); +#+end_src + +#+name: fig:compare_Gf_id_simscape +#+caption: Comparison of the identified transfer function between actuator voltage $V_a$ and sensor stack voltage +#+RESULTS: +[[file:figs/compare_Gf_id_simscape.png]] + +** Conclusion +The obtained parameters $g_a$ and $g_s$ are not consistent between the different methods. + +The one using the experimental data are saved and further used. +#+begin_src matlab :tangle no + save('./matlab/mat/apa95ml_params.mat', 'ga', 'gs'); +#+end_src + +#+begin_src matlab :exports none :eval no + save('./mat/apa95ml_params.mat', 'ga', 'gs'); +#+end_src + * Simscape model of the test-bench :PROPERTIES: :header-args:matlab+: :tangle matlab/simscape_model.m :header-args:matlab+: :comments org :mkdirp yes :END: <> + ** Introduction :ignore: The idea here is to model the test-bench using Simscape. Whereas the suspended mass and metrology frame can be considered as rigid bodies in the frequency range of interest, the Amplified Piezoelectric Actuator (APA) is flexible. -To model the APA, a Finite Element Model (FEM) is used and imported into Simscape. +To model the APA, a Finite Element Model (FEM) is used (Figure [[fig:APA95ML_FEM]]) and imported into Simscape. + +#+name: fig:APA95ML_FEM +#+caption: Finite Element Model of the APA95ML +[[file:figs/APA95ML_FEM.png]] ** Matlab Init :noexport:ignore: #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) @@ -101,8 +525,8 @@ To model the APA, a Finite Element Model (FEM) is used and imported into Simscap ** Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates We first extract the stiffness and mass matrices. #+begin_src matlab - K = extractMatrix('APA95ML_K.txt'); - M = extractMatrix('APA95ML_M.txt'); + K = readmatrix('APA95ML_K.CSV'); + M = readmatrix('APA95ML_M.CSV'); #+end_src #+begin_src matlab :exports results :results value table replace :tangle no @@ -111,16 +535,16 @@ We first extract the stiffness and mass matrices. #+caption: First 10x10 elements of the Stiffness matrix #+RESULTS: -| 300000000.0 | -30000.0 | 8000.0 | -200.0 | -30.0 | -60000.0 | 20000000.0 | -4000.0 | 500.0 | 8 | -| -30000.0 | 100000000.0 | 400.0 | 30.0 | 200.0 | -1 | 4000.0 | -8000000.0 | 800.0 | 7 | -| 8000.0 | 400.0 | 50000000.0 | -800000.0 | -300.0 | -40.0 | 300.0 | 100.0 | 5000000.0 | 40000.0 | -| -200.0 | 30.0 | -800000.0 | 20000.0 | 5 | 1 | -10.0 | -2 | -40000.0 | -300.0 | -| -30.0 | 200.0 | -300.0 | 5 | 40000.0 | 0.3 | -4 | -10.0 | 40.0 | 0.4 | -| -60000.0 | -1 | -40.0 | 1 | 0.3 | 3000.0 | 7000.0 | 0.8 | -1 | 0.0003 | -| 20000000.0 | 4000.0 | 300.0 | -10.0 | -4 | 7000.0 | 300000000.0 | 20000.0 | 3000.0 | 80.0 | -| -4000.0 | -8000000.0 | 100.0 | -2 | -10.0 | 0.8 | 20000.0 | 100000000.0 | -4000.0 | -100.0 | -| 500.0 | 800.0 | 5000000.0 | -40000.0 | 40.0 | -1 | 3000.0 | -4000.0 | 50000000.0 | 800000.0 | -| 8 | 7 | 40000.0 | -300.0 | 0.4 | 0.0003 | 80.0 | -100.0 | 800000.0 | 20000.0 | +| 300000000.0 | -1000.0 | -30000.0 | -40.0 | 70000.0 | 300.0 | 20000000.0 | -30.0 | -5000.0 | 5 | +| -1000.0 | 50000000.0 | -7000.0 | 800000.0 | -20.0 | 300.0 | 3000.0 | 5000000.0 | 400.0 | -40000.0 | +| -30000.0 | -7000.0 | 100000000.0 | -200.0 | -60.0 | 70.0 | 3000.0 | 3000.0 | -8000000.0 | -30.0 | +| -40.0 | 800000.0 | -200.0 | 20000.0 | -0.4 | 4 | 30.0 | 40000.0 | 7 | -300.0 | +| 70000.0 | -20.0 | -60.0 | -0.4 | 3000.0 | 1 | -6000.0 | 10.0 | 8 | -0.1 | +| 300.0 | 300.0 | 70.0 | 4 | 1 | 40000.0 | -10.0 | -10.0 | 30.0 | 0.1 | +| 20000000.0 | 3000.0 | 3000.0 | 30.0 | -6000.0 | -10.0 | 300000000.0 | 2000.0 | 9000.0 | -100.0 | +| -30.0 | 5000000.0 | 3000.0 | 40000.0 | 10.0 | -10.0 | 2000.0 | 50000000.0 | -3000.0 | -800000.0 | +| -5000.0 | 400.0 | -8000000.0 | 7 | 8 | 30.0 | 9000.0 | -3000.0 | 100000000.0 | 100.0 | +| 5 | -40000.0 | -30.0 | -300.0 | -0.1 | 0.1 | -100.0 | -800000.0 | 100.0 | 20000.0 | #+begin_src matlab :exports results :results value table replace :tangle no @@ -129,16 +553,16 @@ We first extract the stiffness and mass matrices. #+caption: First 10x10 elements of the Mass matrix #+RESULTS: -| 0.03 | 2e-06 | -2e-07 | 1e-08 | 2e-08 | 0.0002 | -0.001 | 2e-07 | -8e-08 | -9e-10 | -| 2e-06 | 0.02 | -5e-07 | 7e-09 | 3e-08 | 2e-08 | -3e-07 | 0.0003 | -1e-08 | 1e-10 | -| -2e-07 | -5e-07 | 0.02 | -9e-05 | 4e-09 | -1e-08 | 2e-07 | -2e-08 | -0.0006 | -5e-06 | -| 1e-08 | 7e-09 | -9e-05 | 1e-06 | 6e-11 | 4e-10 | -1e-09 | 3e-11 | 5e-06 | 3e-08 | -| 2e-08 | 3e-08 | 4e-09 | 6e-11 | 1e-06 | 2e-10 | -2e-09 | 2e-10 | -7e-09 | -4e-11 | -| 0.0002 | 2e-08 | -1e-08 | 4e-10 | 2e-10 | 2e-06 | -2e-06 | -1e-09 | -7e-10 | -9e-12 | -| -0.001 | -3e-07 | 2e-07 | -1e-09 | -2e-09 | -2e-06 | 0.03 | -2e-06 | -1e-07 | -5e-09 | -| 2e-07 | 0.0003 | -2e-08 | 3e-11 | 2e-10 | -1e-09 | -2e-06 | 0.02 | -8e-07 | -1e-08 | -| -8e-08 | -1e-08 | -0.0006 | 5e-06 | -7e-09 | -7e-10 | -1e-07 | -8e-07 | 0.02 | 9e-05 | -| -9e-10 | 1e-10 | -5e-06 | 3e-08 | -4e-11 | -9e-12 | -5e-09 | -1e-08 | 9e-05 | 1e-06 | +| 0.03 | 7e-08 | 2e-06 | -3e-09 | -0.0002 | -6e-08 | -0.001 | 8e-07 | 6e-07 | -8e-09 | +| 7e-08 | 0.02 | -1e-06 | 9e-05 | -3e-09 | -4e-09 | -1e-06 | -0.0006 | -4e-08 | 5e-06 | +| 2e-06 | -1e-06 | 0.02 | -3e-08 | -4e-08 | 1e-08 | 1e-07 | -2e-07 | 0.0003 | 1e-09 | +| -3e-09 | 9e-05 | -3e-08 | 1e-06 | -3e-11 | -3e-13 | -7e-09 | -5e-06 | -3e-10 | 3e-08 | +| -0.0002 | -3e-09 | -4e-08 | -3e-11 | 2e-06 | 6e-10 | 2e-06 | -7e-09 | -2e-09 | 7e-11 | +| -6e-08 | -4e-09 | 1e-08 | -3e-13 | 6e-10 | 1e-06 | 1e-08 | 3e-09 | -2e-09 | 2e-13 | +| -0.001 | -1e-06 | 1e-07 | -7e-09 | 2e-06 | 1e-08 | 0.03 | 4e-08 | -2e-06 | 8e-09 | +| 8e-07 | -0.0006 | -2e-07 | -5e-06 | -7e-09 | 3e-09 | 4e-08 | 0.02 | -9e-07 | -9e-05 | +| 6e-07 | -4e-08 | 0.0003 | -3e-10 | -2e-09 | -2e-09 | -2e-06 | -9e-07 | 0.02 | 2e-08 | +| -8e-09 | 5e-06 | 1e-09 | 3e-08 | 7e-11 | 2e-13 | 8e-09 | -9e-05 | 2e-08 | 1e-06 | Then, we extract the coordinates of the interface nodes. @@ -146,100 +570,47 @@ Then, we extract the coordinates of the interface nodes. [int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('APA95ML_out_nodes_3D.txt'); #+end_src +The interface nodes are shown in Figure [[fig:APA95ML_nodes_1]] and their coordinates are listed in Table [[tab:apa95ml_nodes_coordinates]]. + #+begin_src matlab :exports results :results value table replace :tangle no data2orgtable([length(n_i); length(int_i); size(M,1) - 6*length(int_i); size(M,1)], {'Total number of Nodes', 'Number of interface Nodes', 'Number of Modes', 'Size of M and K matrices'}, {}, ' %.0f '); #+end_src #+RESULTS: -| Total number of Nodes | 168959 | -| Number of interface Nodes | 13 | -| Number of Modes | 30 | -| Size of M and K matrices | 108 | +| Total number of Nodes | 7 | +| Number of interface Nodes | 7 | +| Number of Modes | 6 | +| Size of M and K matrices | 48 | #+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*) data2orgtable([[1:length(int_i)]', int_i, int_xyz], {}, {'Node i', 'Node Number', 'x [m]', 'y [m]', 'z [m]'}, ' %f '); #+end_src +#+name: tab:apa95ml_nodes_coordinates #+caption: Coordinates of the interface nodes #+RESULTS: -| Node i | Node Number | x [m] | y [m] | z [m] | -|--------+-------------+--------+-------+-------| -| 1.0 | 168947.0 | 0.0 | 0.03 | 0.0 | -| 2.0 | 168949.0 | 0.0 | -0.03 | 0.0 | -| 3.0 | 168950.0 | -0.035 | 0.0 | 0.0 | -| 4.0 | 168951.0 | -0.028 | 0.0 | 0.0 | -| 5.0 | 168952.0 | -0.021 | 0.0 | 0.0 | -| 6.0 | 168953.0 | -0.014 | 0.0 | 0.0 | -| 7.0 | 168954.0 | -0.007 | 0.0 | 0.0 | -| 8.0 | 168955.0 | 0.0 | 0.0 | 0.0 | -| 9.0 | 168956.0 | 0.007 | 0.0 | 0.0 | -| 10.0 | 168957.0 | 0.014 | 0.0 | 0.0 | -| 11.0 | 168958.0 | 0.021 | 0.0 | 0.0 | -| 12.0 | 168959.0 | 0.035 | 0.0 | 0.0 | -| 13.0 | 168960.0 | 0.028 | 0.0 | 0.0 | +| Node i | Node Number | x [m] | y [m] | z [m] | +|--------+-------------+--------+-------+-----------| +| 1.0 | 40467.0 | 0.0 | 0.0 | 0.029997 | +| 2.0 | 40469.0 | 0.0 | 0.0 | -0.029997 | +| 3.0 | 40470.0 | -0.035 | 0.0 | 0.0 | +| 4.0 | 40475.0 | -0.015 | 0.0 | 0.0 | +| 5.0 | 40476.0 | -0.005 | 0.0 | 0.0 | +| 6.0 | 40477.0 | 0.015 | 0.0 | 0.0 | +| 7.0 | 40478.0 | 0.035 | 0.0 | 0.0 | + +#+name: fig:APA95ML_nodes_1 +#+caption: Interface Nodes chosen for the APA95ML +[[file:figs/APA95ML_nodes_1.png]] Using =K=, =M= and =int_xyz=, we can use the =Reduced Order Flexible Solid= simscape block. -** Piezoelectric parameters -In order to make the conversion from applied voltage to generated force or from the strain to the generated voltage, we need to defined some parameters corresponding to the piezoelectric material: -#+begin_src matlab - d33 = 300e-12; % Strain constant [m/V] - n = 80; % Number of layers per stack - eT = 1.6e-8; % Permittivity under constant stress [F/m] - sD = 1e-11; % Compliance under constant electric displacement [m2/N] - ka = 235e6; % Stack stiffness [N/m] - C = 5e-6; % Stack capactiance [F] -#+end_src - -The ratio of the developed force to applied voltage is: -#+name: eq:piezo_voltage_to_force -\begin{equation} - F_a = g_a V_a, \quad g_a = d_{33} n k_a -\end{equation} -where: -- $F_a$: developed force in [N] -- $n$: number of layers of the actuator stack -- $d_{33}$: strain constant in [m/V] -- $k_a$: actuator stack stiffness in [N/m] -- $V_a$: applied voltage in [V] - -If we take the numerical values, we obtain: -#+begin_src matlab :results replace value - d33*n*ka % [N/V] -#+end_src - -#+RESULTS: -: 5.64 - -From cite:fleming14_desig_model_contr_nanop_system (page 123), the relation between relative displacement of the sensor stack and generated voltage is: -#+name: eq:piezo_strain_to_voltage -\begin{equation} - V_s = \frac{d_{33}}{\epsilon^T s^D n} \Delta h -\end{equation} -where: -- $V_s$: measured voltage in [V] -- $d_{33}$: strain constant in [m/V] -- $\epsilon^T$: permittivity under constant stress in [F/m] -- $s^D$: elastic compliance under constant electric displacement in [m^2/N] -- $n$: number of layers of the sensor stack -- $\Delta h$: relative displacement in [m] - -If we take the numerical values, we obtain: -#+begin_src matlab :results replace value - 1e-6*d33/(eT*sD*n) % [V/um] -#+end_src - -#+RESULTS: -: 23.438 - ** Simscape Model The flexible element is imported using the =Reduced Order Flexible Solid= Simscape block. To model the actuator, an =Internal Force= block is added between the nodes 3 and 12. A =Relative Motion Sensor= block is added between the nodes 1 and 2 to measure the displacement and the amplified piezo. -- [ ] Add schematic of the model with interface nodes - One mass is fixed at one end of the piezo-electric stack actuator, the other end is fixed to the world frame. #+begin_src matlab m = 5; @@ -354,236 +725,6 @@ The obtained dynamics is shown in Figure [[fig:dynamics_force_force_sensor_comp_ save('mat/fem_simscape_models.mat', 'Ghm', 'Gfm') #+end_src -* TODO Estimation of piezoelectric parameters -** Matlab Init :noexport:ignore: -#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) - <> -#+end_src - -#+begin_src matlab :exports none :results silent :noweb yes - <> -#+end_src - -#+begin_src matlab :tangle no - addpath('./matlab/mat/'); - addpath('./matlab/'); -#+end_src - -#+begin_src matlab :eval no - addpath('./mat/'); -#+end_src - -** From actuator voltage to vertical displacement -The data from the "noise test" and the identification test are loaded. -#+begin_src matlab - load('apa95ml_5kg_Amp_E505.mat', 't', 'um', 'y'); -#+end_src - -Any offset value is removed: -#+begin_src matlab - um = detrend(um, 0); % Amplifier Input Voltage [V] - y = detrend(y , 0); % Mass displacement [m] -#+end_src - -Now we add a factor 10 to take into account the gain of the voltage amplifier. -#+begin_src matlab - um = 10*um; % Stack Actuator Input Voltage [V] -#+end_src - -#+begin_src matlab - Ts = t(end)/(length(t)-1); - Fs = 1/Ts; - - win = hanning(ceil(1*Fs)); -#+end_src - -#+begin_src matlab - [tf_est, f] = tfestimate(um, y, win, [], [], 1/Ts); -#+end_src - -The gain from input voltage of the stack to the vertical displacement is determined: -#+begin_src matlab - gD = 4e-7; % [m/V] -#+end_src - -#+begin_src matlab :exports none - freqs = logspace(0, 4, 1000); - - figure; - hold on; - plot(f, abs(tf_est)) - plot([f(2) f(end)], [gD gD]) - hold off; - set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); - ylabel('Amplitude [m/V]'); set(gca, 'XTickLabel', []); - xlim([10, 5e3]); -#+end_src - -#+begin_src matlab - K = extractMatrix('APA95ML_K.txt'); - M = extractMatrix('APA95ML_M.txt'); - [int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('APA95ML_out_nodes_3D.txt'); -#+end_src - -Define parameters just for the simulation. Should not change any results -#+begin_src matlab - d33 = 1; % Strain constant [m/V] - n = 1; % Number of layers per stack - eT = 1; % Permittivity under constant stress [F/m] - sD = 1; % Compliance under constant electric displacement [m2/N] - ka = 1; % Stack stiffness [N/m] - C = 1; % Stack capactiance [F] -#+end_src - -#+begin_src matlab - m = 5; - - %% Name of the Simulink File - mdl = 'piezo_amplified_3d'; - - %% Input/Output definition - clear io; io_i = 1; - 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] - - Gd = linearize(mdl, io); -#+end_src - -#+begin_src matlab :results value replace - gF = abs(dcgain(Gd)); % [m/N] - ans = gF -#+end_src - -#+RESULTS: -: 6.1695e-09 - -\[ g_a = g_D/g_F \] -in [N/V] -#+begin_src matlab :results value replace - ga = gD/gF - ans = ga -#+end_src - -#+RESULTS: -: 64.835 - -#+begin_src matlab :results value replace - na = 2; ns = 1; d33 = 300e-12; n = 80; ka = 235e6; gL = 1.7; - na/(na+ns)*gL*d33*n*ka -#+end_src - -#+RESULTS: -: 6.392 - -- [ ] Why is there a factor 10 with the "theoretical estimation?" - -#+begin_src matlab :exports none - freqs = logspace(1, 4, 1000); - - figure; - hold on; - plot(f, abs(tf_est)) - plot(f, ga*abs(squeeze(freqresp(Gd, f, 'Hz')))) - set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); - ylabel('Amplitude'); xlabel('Frequency [Hz]'); - hold off; -#+end_src - -** From actuator voltage to sensor Voltage -#+begin_src matlab - load('apa95ml_5kg_2a_1s.mat', 't', 'u', 'v'); -#+end_src - -#+begin_src matlab - u = detrend(u, 0); % Input Voltage of the Amplifier [V] - v = detrend(v, 0); % Voltage accross the stack sensor [V] -#+end_src - -#+begin_src matlab - u = 20*u; % Input Voltage of the Amplifier [V] -#+end_src - -#+begin_src matlab - Ts = t(end)/(length(t)-1); - Fs = 1/Ts; - - win = hann(ceil(10/Ts)); - - [tf_est, f] = tfestimate(u, v, win, [], [], 1/Ts); -#+end_src - -#+begin_src matlab - gV = 0.022; % [V/V] -#+end_src - -#+begin_src matlab :exports none - freqs = logspace(1, 4, 1000); - - figure; - hold on; - plot(f, abs(tf_est)) - plot([f(2); f(end)], [gV; gV]) - hold off; - set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); - ylabel('Amplitude'); xlabel('Frequency [Hz]'); - ylim([1e-3, 1e1]); -#+end_src - -#+begin_src matlab - m = 5; - - %% Name of the Simulink File - mdl = 'piezo_amplified_3d'; - - %% Input/Output definition - clear io; io_i = 1; - io(io_i) = linio([mdl, '/Fa'], 1, 'openinput'); io_i = io_i + 1; % Actuator Force [N] - io(io_i) = linio([mdl, '/dL'], 1, 'openoutput'); io_i = io_i + 1; % Sensor Stack displacement [m] - - Gf = linearize(mdl, io); -#+end_src - -$g_F$ in [m/N] -#+begin_src matlab :results value replace - gF = abs(dcgain(Gf)); - ans = gF -#+end_src - -#+RESULTS: -: 2.1546e-10 - -Finally, we compute the gain from strain of the sensor stack to its generated voltage: -\[ g_S = \frac{g_V}{g_F g_a} \] -Gs in [V/m] -#+begin_src matlab :results value replace - gS = gV/gF/ga; - ans = gS -#+end_src - -#+RESULTS: -: 15749000.0 - -#+begin_src matlab :results value replace - d33 = 300e-12; eT = 5.3e-9; sD = 2e-11; n = 80; - d33/(eT*sD*n) -#+end_src - -#+RESULTS: -: 35377000.0 - -#+begin_src matlab :exports none - freqs = logspace(1, 4, 1000); - - figure; - hold on; - plot(f, abs(tf_est)) - plot(f, ga*gS*abs(squeeze(freqresp(Gf, f, 'Hz')))) - set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); - ylabel('Amplitude'); xlabel('Frequency [Hz]'); - hold off; - ylim([1e-3, 1e1]); -#+end_src - * Huddle Test :PROPERTIES: :header-args:matlab+: :tangle matlab/huddle_test.m