test-bench-apa/matlab/piezo_parameters.m

318 lines
8.8 KiB
Matlab

%% Clear Workspace and Close figures
clear; close all; clc;
%% Intialize Laplace variable
s = zpk('s');
addpath('./mat/');
% Estimation from Data-sheet
% <<sec:estimation_datasheet>>
% 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}
ka = 235e6; % [N/m]
Lmax = 20e-6; % [m]
Vmax = 170; % [V]
ka*Lmax/Vmax % [N/V]
% Estimation from Piezoelectric parameters
% <<sec:estimation_piezo_params>>
% 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:
% #+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:
ga = d33*n*ka; % [N/V]
% #+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:
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]
% 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]
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');
% #+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.
% 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]
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]
% #+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}
ga = g_d_Va/g_d_Fa;
% #+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]].
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');
% 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 [[fig:gain_Va_to_Vs]]).
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]
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');
% #+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.
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));
% #+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}
gs = g_Vs_Va/G_dh_Va; % [V/m]
% #+RESULTS:
% : gs = 3.5 [V/um]
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');
save('./mat/apa95ml_params.mat', 'ga', 'gs');