diff --git a/figs/detail_fem_apa95ml_compliance.pdf b/figs/detail_fem_apa95ml_compliance.pdf index 6ad228f..8ef119f 100644 Binary files a/figs/detail_fem_apa95ml_compliance.pdf and b/figs/detail_fem_apa95ml_compliance.pdf differ diff --git a/figs/detail_fem_apa95ml_compliance.png b/figs/detail_fem_apa95ml_compliance.png index 3322569..db0adca 100644 Binary files a/figs/detail_fem_apa95ml_compliance.png and b/figs/detail_fem_apa95ml_compliance.png differ diff --git a/matlab/detail_fem_1_flexible_body.m b/matlab/detail_fem_1_flexible_body.m new file mode 100644 index 0000000..3935bbd --- /dev/null +++ b/matlab/detail_fem_1_flexible_body.m @@ -0,0 +1,307 @@ +%% Clear Workspace and Close figures +clear; close all; clc; + +%% Intialize Laplace variable +s = zpk('s'); + +%% Path for functions, data and scripts +addpath('./src/'); % Path for scripts +addpath('./mat/'); % Path for data +addpath('./STEPS/'); % Path for Simscape Model +addpath('./subsystems/'); % Path for Subsystems Simulink files + +%% Linearization options +opts = linearizeOptions; +opts.SampleTime = 0; + +%% Open Simscape Model +mdl = 'detail_fem_super_element'; % Name of the Simulink File +open(mdl); % Open Simscape Model + +%% Colors for the figures +colors = colororder; +freqs = logspace(1,3,500); % Frequency vector [Hz] + +%% Estimate "Sensor Constant" - (THP5H) +d33 = 680e-12; % Strain constant [m/V] +n = 160; % Number of layers per stack +eT = 4500*8.854e-12; % Permittivity under constant stress [F/m] +sD = 21e-12; % Compliance under constant electric displacement [m2/N] + +gs = d33/(eT*sD*n); % Sensor Constant [V/m] + +%% Estimate "Actuator Constant" - (THP5H) +d33 = 680e-12; % Strain constant [m/V] +n = 320; % Number of layers + +cE = 1/sD; % Youngs modulus [N/m^2] +A = (10e-3)^2; % Area of the stacks [m^2] +L = 40e-3; % Length of the two stacks [m] +ka = cE*A/L; % Stiffness of the two stacks [N/m] + +ga = d33*n*ka; % Actuator Constant [N/V] + +%% Load reduced order model +K = readmatrix('APA95ML_K.CSV'); % order: 48 +M = readmatrix('APA95ML_M.CSV'); +[int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('APA95ML_out_nodes_3D.txt'); + +%% Stiffness estimation +m = 0.0001; % block-free condition, no payload +k_support = 1e9; +c_support = 1e3; + +clear io; io_i = 1; +io(io_i) = linio([mdl, '/Fd'], 1, 'openinput'); io_i = io_i + 1; +io(io_i) = linio([mdl, '/y'], 1, 'openoutput'); io_i = io_i + 1; + +G = linearize(mdl, io); + +% The inverse of the DC gain of the transfer function +% from vertical force to vertical displacement is the axial stiffness of the APA +k_est = 1/dcgain(G); % [N/m] + +%% Estimated compliance of the APA95ML +freqs = logspace(2, log10(5000), 1000); + +% Get first resonance indice +i_max = find(abs(squeeze(freqresp(G, freqs(2:end), 'Hz'))) - abs(squeeze(freqresp(G, freqs(1:end-1), 'Hz'))) < 0, 1); + +figure; +hold on; +plot(freqs, abs(squeeze(freqresp(G, freqs, 'Hz'))), 'DisplayName', 'Compliance'); +plot([freqs(1), freqs(end)], [1/k_est, 1/k_est], 'k--', 'DisplayName', sprintf('$1/k$ ($k = %.0f N/\\mu m$)', 1e-6*k_est)) +xline(freqs(i_max), '--', 'linewidth', 1, 'color', [0,0,0], 'DisplayName', sprintf('$f_0 = %.0f$ Hz', freqs(i_max))) +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Amplitude [m/N]'); +leg = legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; +xlim([100, 5000]); + +%% Estimation of the amplification factor and Stroke +clear io; io_i = 1; +io(io_i) = linio([mdl, '/Fa'], 1, 'openinput'); io_i = io_i + 1; +io(io_i) = linio([mdl, '/y'], 1, 'openoutput'); io_i = io_i + 1; +io(io_i) = linio([mdl, '/d'], 1, 'openoutput'); io_i = io_i + 1; + +G = linearize(mdl, io); + +% Estimated amplification factor +ampl_factor = abs(dcgain(G(1,1))./dcgain(G(2,1))); + +% Estimated stroke +apa_stroke = ampl_factor * 3 * 20e-6; % [m] + +%% Experimental plant identification +% with PD200 amplifier (gain of 20) - 2 stacks as an actuator, 1 as a sensor +load('apa95ml_5kg_2a_1s.mat') + +Va = 20*u; % Voltage amplifier gain: 20 + +% Spectral Analysis parameters +Ts = t(end)/(length(t)-1); +Nfft = floor(1/Ts); +win = hanning(Nfft); +Noverlap = floor(Nfft/2); + +% Identification of the transfer function from Va to di +[G_y, f] = tfestimate(detrend(Va, 0), detrend(y, 0), win, Noverlap, Nfft, 1/Ts); +[G_Vs, ~] = tfestimate(detrend(Va, 0), detrend(v, 0), win, Noverlap, Nfft, 1/Ts); + +%% Plant Identification from Multi-Body model +% Load Reduced Order Matrices +K = readmatrix('APA95ML_K.CSV'); % order: 48 +M = readmatrix('APA95ML_M.CSV'); +[int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('APA95ML_out_nodes_3D.txt'); + +m = 5.5; % Mass of the suspended granite [kg] +k_support = 4e7; +c_support = 3e2; + +% Compute transfer functions +clear io; io_i = 1; +io(io_i) = linio([mdl, '/Va'], 1, 'openinput'); io_i = io_i + 1; % Voltage accros piezo stacks [V] +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; % Sensor stack voltage [V] + +Gm = linearize(mdl, io); +Gm.InputName = {'Va'}; +Gm.OutputName = {'y', 'Vs'}; + +%% Comparison of the identified transfer function from Va to di to the multi-body model +freqs = logspace(1, 3, 500); +figure; +tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +plot(f, abs(G_y), '-', 'color', [colors(2,:), 0.5], 'linewidth', 2.5, 'DisplayName', 'Measured FRF'); +plot(freqs, abs(squeeze(freqresp(Gm('y', 'Va'), freqs, 'Hz'))), '--', 'color', colors(2,:), 'DisplayName', 'Model') +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude $y/V_a$ [m/V]'); set(gca, 'XTickLabel',[]); +hold off; +ylim([1e-8, 1e-5]); +leg = legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; + +ax2 = nexttile; +hold on; +plot(f, 180/pi*angle(G_y), '-', 'color' , [colors(2,:), 0.5], 'linewidth', 2.5); +plot(freqs, 180/pi*angle(squeeze(freqresp(Gm('y', 'Va'), freqs, 'Hz'))), '--', 'color', colors(2,:)) +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +hold off; +yticks(-360:45:360); +ylim([-45, 180]); + +linkaxes([ax1,ax2],'x'); +xlim([10, 1e3]); + +%% Comparison of the identified transfer function from Va to Vs to the multi-body model +figure; +tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +plot(f, abs(G_Vs), '-', 'color', [colors(1,:), 0.5], 'linewidth', 2.5, 'DisplayName', 'Measured FRF'); +plot(freqs, abs(squeeze(freqresp(Gm('Vs', 'Va'), freqs, 'Hz'))), '--', 'color', colors(1,:), 'DisplayName', 'Model') +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]); +hold off; +ylim([1e-3, 1e1]); +leg = legend('location', 'northwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; + +ax2 = nexttile; +hold on; +plot(f, 180/pi*angle(G_Vs), '-', 'color', [colors(1,:), 0.5], 'linewidth', 2.5); +plot(freqs, 180/pi*angle(squeeze(freqresp(Gm('Vs', 'Va'), freqs, 'Hz'))), '--', 'color', colors(1,:)) +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +hold off; +yticks(-360:90:360); ylim([-180, 180]); + +linkaxes([ax1,ax2],'x'); +xlim([10, 1e3]); + +%% Integral Force Feedback Controller +K_iff = (1/(s + 2*2*pi))*(s/(s + 0.5*2*pi)); +K_iff.inputname = {'Vs'}; +K_iff.outputname = {'u_iff'}; + +% New damped plant input +S1 = sumblk("u = u_iff + u_damp"); + +% Voltage amplifier with gain of 20 +voltage_amplifier = tf(20); +voltage_amplifier.inputname = {'u'}; +voltage_amplifier.outputname = {'Va'}; + +%% Load experimental data with IFF implemented with different gains +load('apa95ml_iff_test.mat', 'results'); + +% Tested gains +g_iff = [0, 10, 50, 100, 500, 1000]; + +% Spectral Analysis parameters +Ts = t(end)/(length(t)-1); +Nfft = floor(1/Ts); +win = hanning(Nfft); +Noverlap = floor(Nfft/2); + +%% Computed the identified FRF of the damped plants +tf_iff = {zeros(1, length(g_iff))}; +for i=1:length(g_iff) + [tf_est, f] = tfestimate(results{i}.u, results{i}.y, win, Noverlap, Nfft, 1/Ts); + tf_iff(i) = {tf_est}; +end + +%% Estimate the damped plants from the multi-body model +Gm_iff = {zeros(1, length(g_iff))}; + +for i=1:length(g_iff) + K_iff_g = -K_iff*g_iff(i); K_iff_g.inputname = {'Vs'}; K_iff_g.outputname = {'u_iff'}; + Gm_iff(i) = {connect(Gm, K_iff_g, S1, voltage_amplifier, {'u_damp'}, {'y'})}; +end + +%% Identify second order plants from the experimental data +% This is mandatory to estimate the experimental "poles" +% an place them in the root-locus plot +G_id = {zeros(1,length(results))}; + +f_start = 70; % [Hz] +f_end = 500; % [Hz] + +for i = 1:length(results) + tf_id = tf_iff{i}(sum(ff_end)); + f_id = f(sum(ff_end)); + + gfr = idfrd(tf_id, 2*pi*f_id, Ts); + G_id(i) = {procest(gfr,'P2UDZ')}; +end + +%% Comparison of the Root-Locus computed from the multi-body model and the identified closed-loop poles +gains = logspace(0, 5, 1000); +figure; +hold on; +plot(real( pole(Gm('Vs', 'Va'))), imag( pole(Gm('Vs', 'Va'))), 'kx', 'HandleVisibility', 'off'); +plot(real(tzero(Gm('Vs', 'Va'))), imag(tzero(Gm('Vs', 'Va'))), 'ko', 'HandleVisibility', 'off'); +for i = 1:length(gains) + cl_poles = pole(feedback(Gm('Vs', 'Va'), gains(i)*K_iff)); + plot(real(cl_poles(imag(cl_poles)>100)), imag(cl_poles(imag(cl_poles)>100)), 'k.', 'HandleVisibility', 'off'); +end +for i = 1:length(g_iff) + cl_poles = pole(Gm_iff{i}); + plot(real(cl_poles(imag(cl_poles)>100)), imag(cl_poles(imag(cl_poles)>100)), '.', 'MarkerSize', 20, 'color', colors(i,:), 'HandleVisibility', 'off'); + plot(real(pole(G_id{i})), imag(pole(G_id{i})), 'x', 'color', colors(i,:), 'DisplayName', sprintf('g = %0.f', g_iff(i)), 'DisplayName', sprintf('$g = %.0f$', g_iff(i))); +end +xlabel('Real Part'); +ylabel('Imaginary Part'); +axis equal; +ylim([-100, 2100]); +xlim([-2100,100]); +leg = legend('location', 'northwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; + +%% Experimental damped plant for several IFF gains and estimated damped plants from the model +figure; +tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); + +ax1 = nexttile([2, 1]); +hold on; +plot(f, abs(tf_iff{1}), '-', 'DisplayName', '$g = 0$', 'color', [0,0,0, 0.5], 'linewidth', 2.5) +plot(f, abs(squeeze(freqresp(Gm_iff{1}, f, 'Hz'))), 'k--', 'HandleVisibility', 'off') +for i = 2:length(results) + plot(f, abs(tf_iff{i}), '-', 'DisplayName', sprintf('g = %0.f', g_iff(i)), 'color', [colors(i-1,:), 0.5], 'linewidth', 2.5) +end +for i = 2:length(results) + plot(f, abs(squeeze(freqresp(Gm_iff{i}, f, 'Hz'))), '--', 'color', colors(i-1,:), 'HandleVisibility', 'off') +end +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); +ylabel('Amplitude $y/V_a$ [m/N]'); set(gca, 'XTickLabel',[]); +hold off; +ylim([1e-6, 2e-4]); +leg = legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; + +ax2 = nexttile; +hold on; +plot(f, 180/pi*angle(tf_iff{1}./squeeze(freqresp(exp(-s*2e-4), f, 'Hz'))), '-', 'color', [0,0,0, 0.5], 'linewidth', 2.5) +plot(f, 180/pi*angle(squeeze(freqresp(Gm_iff{1}, f, 'Hz'))), 'k--') +for i = 2:length(results) + plot(f, 180/pi*angle(tf_iff{i}./squeeze(freqresp(exp(-s*2e-4), f, 'Hz'))), '-', 'color', [colors(i-1,:), 0.5], 'linewidth', 2.5) + plot(f, 180/pi*angle(squeeze(freqresp(Gm_iff{i}, f, 'Hz'))), '--', 'color', colors(i-1,:)) +end +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); +ylabel('Phase [deg]'); xlabel('Frequency [Hz]'); +hold off; +yticks(-360:45:360); +ylim([-10, 190]); + +linkaxes([ax1,ax2], 'x'); +xlim([150, 500]); diff --git a/matlab/detail_fem_2_actuators.m b/matlab/detail_fem_2_actuators.m new file mode 100644 index 0000000..ab6f873 --- /dev/null +++ b/matlab/detail_fem_2_actuators.m @@ -0,0 +1,318 @@ +%% Clear Workspace and Close figures +clear; close all; clc; + +%% Intialize Laplace variable +s = zpk('s'); + +%% Path for functions, data and scripts +addpath('./src/'); % Path for scripts +addpath('./mat/'); % Path for data +addpath('./STEPS/'); % Path for Simscape Model +addpath('./subsystems/'); % Path for Subsystems Simulink files + +%% Linearization options +opts = linearizeOptions; +opts.SampleTime = 0; + +%% Open Simscape Model +mdl = 'detail_fem_APA300ML'; % Name of the Simulink File +open(mdl); % Open Simscape Model + +% Piezoelectric parameters +ga = -25.9; % [N/V] +gs = -5.08e6; % [V/m] + +%% Colors for the figures +colors = colororder; +freqs = logspace(1,3,500); % Frequency vector [Hz] + +%% Identify dynamics with "Reduced Order Flexible Body" +K = readmatrix('APA300ML_mat_K.CSV'); +M = readmatrix('APA300ML_mat_M.CSV'); +[int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('APA300ML_out_nodes_3D.txt'); + +m = 5; % [kg] +ga = 25.9; % [N/V] +gs = 5.08e6; % [V/m] + +clear io; io_i = 1; +io(io_i) = linio([mdl, '/Va'], 1, 'openinput'); io_i = io_i + 1; +io(io_i) = linio([mdl, '/Fd'], 1, 'openinput'); io_i = io_i + 1; +io(io_i) = linio([mdl, '/z'], 1, 'openoutput'); io_i = io_i + 1; +io(io_i) = linio([mdl, '/Vs'], 1, 'openoutput'); io_i = io_i + 1; + +G_fem = linearize(mdl, io); +G_fem_z = G_fem('z','Va'); +G_fem_Vs = G_fem('Vs', 'Va'); +G_fem_comp = G_fem('z', 'Fd'); + +%% Determine c1 and k1 from the zero +G_zeros = zero(minreal(G_fem_Vs)); +G_zeros = G_zeros(imag(G_zeros)>0); +[~, i_sort] = sort(imag(G_zeros)); +G_zeros = G_zeros(i_sort); +G_zero = G_zeros(1); + +% Solving 2nd order equations +c1 = -2*m*real(G_zero); +k1 = m*(imag(G_zero)^2 + real(G_zero)^2); + +%% Determine ka, ke, ca, ce from the first pole +G_poles = pole(minreal(G_fem_z)); +G_poles = G_poles(imag(G_poles)>0); +[~, i_sort] = sort(imag(G_poles)); +G_poles = G_poles(i_sort); +G_pole = G_poles(1); + +% Solving 2nd order equations +ce = 3*(-2*m*real(G_pole(1)) - c1); +ca = 1/2*ce; + +ke = 3*(m*(imag(G_pole)^2 + real(G_pole)^2) - k1); +ka = 1/2*ke; + +%% Matching sensor/actuator constants +% ga = dcgain(G_fem_z) / (1/(ka + k1*ke/(k1 + ke))); +clear io; io_i = 1; +io(io_i) = linio([mdl, '_2dof', '/Fa'], 1, 'openinput'); io_i = io_i + 1; +io(io_i) = linio([mdl, '_2dof', '/z'], 1, 'openoutput'); io_i = io_i + 1; +ga = dcgain(G_fem_z)/dcgain(linearize([mdl, '_2dof'], io)); + +clear io; io_i = 1; +io(io_i) = linio([mdl, '_2dof', '/Va'], 1, 'openinput'); io_i = io_i + 1; +io(io_i) = linio([mdl, '_2dof', '/dL'], 1, 'openoutput'); io_i = io_i + 1; +gs = dcgain(G_fem_Vs)/dcgain(linearize([mdl, '_2dof'], io)); + +%% Identify dynamics with tuned 2DoF model +clear io; io_i = 1; +io(io_i) = linio([mdl, '_2dof', '/Va'], 1, 'openinput'); io_i = io_i + 1; +io(io_i) = linio([mdl, '_2dof', '/Fd'], 1, 'openinput'); io_i = io_i + 1; +io(io_i) = linio([mdl, '_2dof', '/z'], 1, 'openoutput'); io_i = io_i + 1; +io(io_i) = linio([mdl, '_2dof', '/Vs'], 1, 'openoutput'); io_i = io_i + 1; + +G_2dof = linearize([mdl, '_2dof'], io); +G_2dof_z = G_2dof('z','Va'); +G_2dof_Vs = G_2dof('Vs', 'Va'); +G_2dof_comp = G_2dof('z', 'Fd'); + +%% Comparison of the transfer functions from Va to vertical motion - FEM vs 2DoF +freqs = logspace(1, 3, 500); +figure; +tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +plot(freqs, abs(squeeze(freqresp(G_fem_z, freqs, 'Hz'))), '-', 'color', [colors(2,:), 0.5], 'linewidth', 2.5, 'DisplayName', 'FEM') +plot(freqs, abs(squeeze(freqresp(G_2dof_z, freqs, 'Hz'))), '--', 'color', colors(2,:), 'DisplayName', '2DoF Model') +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude $y/V_a$ [m/V]'); set(gca, 'XTickLabel',[]); +hold off; +ylim([1e-8, 2e-4]); +leg = legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; + +ax2 = nexttile; +hold on; +plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_fem_z, freqs, 'Hz')))), '-', 'color', [colors(2,:), 0.5], 'linewidth', 2.5); +plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_2dof_z, freqs, 'Hz')))), '--', 'color', colors(2,:)) +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +hold off; +yticks(-360:45:360); ylim([-20, 200]); + +linkaxes([ax1,ax2],'x'); +xlim([10, 1e3]); + +%% Comparison of the transfer functions from Va to Vs - FEM vs 2DoF +figure; +tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +plot(freqs, abs(squeeze(freqresp(G_fem_Vs, freqs, 'Hz'))), '-', 'color', [colors(1,:), 0.5], 'linewidth', 2.5, 'DisplayName', 'FEM'); +plot(freqs, abs(squeeze(freqresp(G_2dof_Vs, freqs, 'Hz'))), '--', 'color', colors(1,:), 'DisplayName', '2DoF Model') +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]); +hold off; +ylim([6e-4, 3e1]); +leg = legend('location', 'northwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; + +ax2 = nexttile; +hold on; +plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_fem_Vs, freqs, 'Hz')))), '-', 'color', [colors(1,:), 0.5], 'linewidth', 2.5); +plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_2dof_Vs, freqs, 'Hz')))), '--', 'color', colors(1,:)) +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +hold off; +yticks(-360:45:360); ylim([-20, 200]); + +linkaxes([ax1,ax2],'x'); +xlim([10, 1e3]); + +%% Effect of electrical boundaries on the +oc = load('detail_fem_apa95ml_open_circuit.mat', 't', 'encoder', 'u'); +sc = load('detail_fem_apa95ml_short_circuit.mat', 't', 'encoder', 'u'); + +% Spectral Analysis parameters +Ts = sc.t(end)/(length(sc.t)-1); +Nfft = floor(2/Ts); +win = hanning(Nfft); +Noverlap = floor(Nfft/2); + +% Identification of the transfer function from Va to di +[G_oc, f] = tfestimate(detrend(oc.u, 0), detrend(oc.encoder, 0), win, Noverlap, Nfft, 1/Ts); +[G_sc, f] = tfestimate(detrend(sc.u, 0), detrend(sc.encoder, 0), win, Noverlap, Nfft, 1/Ts); + +% Find resonance frequencies +[~, i_oc] = max(abs(G_oc(f<300))); +[~, i_sc] = max(abs(G_sc(f<300))); + +%% Effect of the electrical bondaries of the force sensor stack on the APA95ML resonance frequency +figure; +tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +plot(f, abs(G_oc), '-', 'DisplayName', sprintf('Open-Circuit - $f_0 = %.1f Hz$', f(i_oc))) +plot(f, abs(G_sc), '-', 'DisplayName', sprintf('Short-Circuit - $f_0 = %.1f Hz$', f(i_sc))) +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); +ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]); +hold off; +ylim([1e-6, 1e-4]); +leg = legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; + +ax2 = nexttile; +hold on; +plot(f, 180/pi*angle(G_oc), '-') +plot(f, 180/pi*angle(G_sc), '-') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); +ylabel('Phase'); xlabel('Frequency [Hz]'); +hold off; +yticks(-360:45:360); +ylim([-20, 200]); +axis padded 'auto x' + +linkaxes([ax1,ax2], 'x'); +xlim([100, 300]); + +%% Compare Dynamics between "Reduced Order" flexible joints and "2-dof and 3-dof" joints +% Let's initialize all the stages with default parameters. +initializeGround('type', 'rigid'); +initializeGranite('type', 'rigid'); +initializeTy('type', 'rigid'); +initializeRy('type', 'rigid'); +initializeRz('type', 'rigid'); +initializeMicroHexapod('type', 'rigid'); +initializeSample('m', 50); + +initializeSimscapeConfiguration(); +initializeDisturbances('enable', false); +initializeLoggingConfiguration('log', 'none'); +initializeController('type', 'open-loop'); +initializeReferences(); + +mdl = 'detail_fem_nass'; + +% Input/Output definition +clear io; io_i = 1; +io(io_i) = linio([mdl, '/Controller'], 1, 'openinput'); io_i = io_i + 1; % Actuator Inputs +io(io_i) = linio([mdl, '/Tracking Error'], 1, 'openoutput', [], 'EdL'); io_i = io_i + 1; % Errors in the frame of the struts +io(io_i) = linio([mdl, '/NASS'], 3, 'openoutput', [], 'fn'); io_i = io_i + 1; % Force Sensors + +% Flexible actuators +initializeSimplifiedNanoHexapod('actuator_type', 'flexible', ... + 'flex_type_F', '2dof', ... + 'flex_type_M', '3dof'); + +G_flex = linearize(mdl, io); +G_flex.InputName = {'f1', 'f2', 'f3', 'f4', 'f5', 'f6'}; +G_flex.OutputName = {'l1', 'l2', 'l3', 'l4', 'l5', 'l6', 'fm1', 'fm2', 'fm3', 'fm4', 'fm5', 'fm6'}; + +% Actuators modeled as 2DoF system +initializeSimplifiedNanoHexapod('actuator_type', 'apa300ml', ... + 'flex_type_F', '2dof', ... + 'flex_type_M', '3dof'); + +G_ideal = linearize(mdl, io); +G_ideal.InputName = {'f1', 'f2', 'f3', 'f4', 'f5', 'f6'}; +G_ideal.OutputName = {'l1', 'l2', 'l3', 'l4', 'l5', 'l6', 'fm1', 'fm2', 'fm3', 'fm4', 'fm5', 'fm6'}; + +%% Comparison of the dynamics for actuators modeled using "reduced order flexible body" and using 2DoF system - HAC plant +freqs = logspace(1, 4, 1000); + +figure; +tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +for j = 1:5 + for k = j+1:6 + plot(freqs, abs(squeeze(freqresp(G_flex("l"+k,"f"+j), freqs, 'Hz'))), 'color', [colors(1,:), 0.1], ... + 'HandleVisibility', 'off'); + plot(freqs, abs(squeeze(freqresp(G_ideal("l"+k,"f"+j), freqs, 'Hz'))), 'color', [colors(2,:), 0.1], ... + 'HandleVisibility', 'off'); + end +end +plot(freqs, abs(squeeze(freqresp(G_flex("l1","f1"), freqs, 'Hz'))), 'color', colors(1,:), 'DisplayName', 'FEM'); +plot(freqs, abs(squeeze(freqresp(G_ideal("l1","f1"), freqs, 'Hz'))), 'color', colors(2,:), 'DisplayName', '2-DoF'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]); +leg = legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; +ylim([1e-10, 1e-4]); + +ax2 = nexttile; +hold on; +plot(freqs, 180/pi*angle(squeeze(freqresp(G_flex("l1","f1"), freqs, 'Hz')))); +plot(freqs, 180/pi*angle(squeeze(freqresp(G_ideal("l1","f1"), freqs, 'Hz')))); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +ylabel('Phase [deg]'); xlabel('Frequency [Hz]'); +ylim([-20, 200]); +yticks([-360:45:360]); + +linkaxes([ax1,ax2],'x'); + +%% Comparison of the dynamics for actuators modeled using "reduced order flexible body" and using 2DoF system - IFF plant +freqs = logspace(0, 3, 1000); + +figure; +tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +for j = 1:5 + for k = j+1:6 + plot(freqs, abs(squeeze(freqresp(G_flex("fm"+k,"f"+j), freqs, 'Hz'))), 'color', [colors(1,:), 0.1], ... + 'HandleVisibility', 'off'); + plot(freqs, abs(squeeze(freqresp(G_ideal("fm"+k,"f"+j), freqs, 'Hz'))), 'color', [colors(2,:), 0.1], ... + 'HandleVisibility', 'off'); + end +end +plot(freqs, abs(squeeze(freqresp(G_flex("fm1","f1"), freqs, 'Hz'))), 'color', colors(1,:), 'DisplayName', 'FEM'); +plot(freqs, abs(squeeze(freqresp(G_ideal("fm1","f1"), freqs, 'Hz'))), 'color', colors(2,:), 'DisplayName', '2-DoF'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude [N/N]'); set(gca, 'XTickLabel',[]); +leg = legend('location', 'northwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; +ylim([1e-5, 1e1]); + +ax2 = nexttile; +hold on; +plot(freqs, 180/pi*angle(squeeze(freqresp(G_flex("fm1","f1"), freqs, 'Hz')))); +plot(freqs, 180/pi*angle(squeeze(freqresp(G_ideal("fm1","f1"), freqs, 'Hz')))); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +ylabel('Phase [deg]'); xlabel('Frequency [Hz]'); +ylim([-20, 200]); +yticks([-360:45:360]); + +linkaxes([ax1,ax2],'x'); diff --git a/matlab/detail_fem_3_flexible_joints.m b/matlab/detail_fem_3_flexible_joints.m new file mode 100644 index 0000000..1cd8176 --- /dev/null +++ b/matlab/detail_fem_3_flexible_joints.m @@ -0,0 +1,623 @@ +%% Clear Workspace and Close figures +clear; close all; clc; + +%% Intialize Laplace variable +s = zpk('s'); + +%% Path for functions, data and scripts +addpath('./src/'); % Path for scripts +addpath('./mat/'); % Path for data +addpath('./STEPS/'); % Path for Simscape Model +addpath('./subsystems/'); % Path for Subsystems Simulink files + +%% Linearization options +opts = linearizeOptions; +opts.SampleTime = 0; + +%% Open Simscape Model +mdl = 'detail_fem_nass'; % Name of the Simulink File +open(mdl); % Open Simscape Model + +%% Colors for the figures +colors = colororder; +freqs = logspace(1,3,500); % Frequency vector [Hz] + +%% Identify the dynamics for several considered bending stiffnesses +% Let's initialize all the stages with default parameters. +initializeGround('type', 'rigid'); +initializeGranite('type', 'rigid'); +initializeTy('type', 'rigid'); +initializeRy('type', 'rigid'); +initializeRz('type', 'rigid'); +initializeMicroHexapod('type', 'rigid'); +initializeSample('m', 50); + +initializeSimscapeConfiguration(); +initializeDisturbances('enable', false); +initializeLoggingConfiguration('log', 'none'); +initializeController('type', 'open-loop'); +initializeReferences(); + +% Input/Output definition +clear io; io_i = 1; +io(io_i) = linio([mdl, '/Controller'], 1, 'openinput'); io_i = io_i + 1; % Actuator Inputs +io(io_i) = linio([mdl, '/Tracking Error'], 1, 'openoutput', [], 'EdL'); io_i = io_i + 1; % Errors in the frame of the struts +io(io_i) = linio([mdl, '/NASS'], 3, 'openoutput', [], 'fn'); io_i = io_i + 1; % Force Sensors + +% Effect of bending stiffness +Kf = [0, 50, 100, 500]; % [Nm/rad] +G_Kf = {zeros(length(Kf), 1)}; + +for i = 1:length(Kf) + % Limited joint axial compliance + initializeSimplifiedNanoHexapod('actuator_type', '1dof', ... + 'flex_type_F', '2dof', ... + 'flex_type_M', '3dof', ... + 'actuator_k', 1e6, ... + 'actuator_c', 1e1, ... + 'actuator_kp', 0, ... + 'actuator_cp', 0, ... + 'Fsm', 56e-3, ... % APA300ML weight 112g + 'Msm', 56e-3, ... + 'Kf_F', Kf(i), ... + 'Kf_M', Kf(i)); + + G_Kf(i) = {linearize(mdl, io)}; + G_Kf{i}.InputName = {'f1', 'f2', 'f3', 'f4', 'f5', 'f6'}; + G_Kf{i}.OutputName = {'l1', 'l2', 'l3', 'l4', 'l5', 'l6', 'fm1', 'fm2', 'fm3', 'fm4', 'fm5', 'fm6'}; +end + +freqs = logspace(0, 3, 1000); + +%% Effect of the flexible joint bending stiffness on the HAC-plant +figure; +tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +for i = 1:length(Kf) + for j = 1:5 + for k = j+1:6 + plot(freqs, abs(squeeze(freqresp(G_Kf{i}("l"+k,"f"+j), freqs, 'Hz'))), 'color', [colors(i,:), 0.1], ... + 'HandleVisibility', 'off'); + end + end + plot(freqs, abs(squeeze(freqresp(G_Kf{i}("l1","f1"), freqs, 'Hz'))), 'color', colors(i,:), 'DisplayName', sprintf('$k_f = %.0f $ [Nm/rad]', Kf(i))); + for j = 2:6 + plot(freqs, abs(squeeze(freqresp(G_Kf{i}("l"+j,"f"+j), freqs, 'Hz'))), 'color', colors(i,:), 'HandleVisibility', 'off'); + end +end +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]); +leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; +ylim([1e-10, 1e-4]); + +ax2 = nexttile; +hold on; +for i = 1:length(Kf) + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_Kf{i}(1, 1), freqs, 'Hz'))))); +end +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +ylabel('Phase [deg]'); xlabel('Frequency [Hz]'); +ylim([-200, 20]); +yticks([-360:45:360]); + +linkaxes([ax1,ax2],'x'); + +%% Effect of the flexible joint bending stiffness on the IFF plant +figure; +tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +for i = 1:length(Kf) + for j = 1:5 + for k = j+1:6 + plot(freqs, abs(squeeze(freqresp(G_Kf{i}("fm"+k,"f"+j), freqs, 'Hz'))), 'color', [colors(i,:), 0.1], ... + 'HandleVisibility', 'off'); + end + end + plot(freqs, abs(squeeze(freqresp(G_Kf{i}("fm1","f1"), freqs, 'Hz'))), 'color', colors(i,:), 'DisplayName', sprintf('$k_f = %.0f $ [Nm/rad]', Kf(i))); + for j = 2:6 + plot(freqs, abs(squeeze(freqresp(G_Kf{i}("fm"+j,"f"+j), freqs, 'Hz'))), 'color', colors(i,:), 'HandleVisibility', 'off'); + end +end +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude [N/N]'); set(gca, 'XTickLabel',[]); +leg = legend('location', 'northwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; +ylim([1e-4, 1e2]); + +ax2 = nexttile(); +hold on; +for i = 1:length(Kf) + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_Kf{i}("fm1", "f1"), freqs, 'Hz'))))); +end +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +ylabel('Phase [deg]'); xlabel('Frequency [Hz]'); +ylim([-20, 200]); +yticks([-360:45:360]); + +linkaxes([ax1,ax2],'x'); + +%% Decentalized IFF +Kiff = -200 * ... % Gain + 1/s * ... % LPF: provides integral action + eye(6); % Diagonal 6x6 controller (i.e. decentralized) + +Kiff.InputName = {'fm1', 'fm2', 'fm3', 'fm4', 'fm5', 'fm6'}; +Kiff.OutputName = {'f1', 'f2', 'f3', 'f4', 'f5', 'f6'}; + +%% Root Locus for decentralized IFF - 1dof actuator - Effect of joint bending stiffness +gains = logspace(-1, 2, 400); + +figure; +tiledlayout(1, 1, 'TileSpacing', 'compact', 'Padding', 'None'); +nexttile(); +hold on; + +for i = 1:length(Kf) + plot(real(pole(G_Kf{i}({"fm1", "fm2", "fm3", "fm4", "fm5", "fm6"}, {"f1", "f2", "f3", "f4", "f5", "f6"}))), imag(pole(G_Kf{i}({"fm1", "fm2", "fm3", "fm4", "fm5", "fm6"}, {"f1", "f2", "f3", "f4", "f5", "f6"}))), 'x', 'color', colors(i,:), ... + 'DisplayName', sprintf('$k_f = %.0f$ Nm/rad', Kf(i))); + plot(real(tzero(G_Kf{i}({"fm1", "fm2", "fm3", "fm4", "fm5", "fm6"}, {"f1", "f2", "f3", "f4", "f5", "f6"}))), imag(tzero(G_Kf{i}({"fm1", "fm2", "fm3", "fm4", "fm5", "fm6"}, {"f1", "f2", "f3", "f4", "f5", "f6"}))), 'o', 'color', colors(i,:), ... + 'HandleVisibility', 'off'); + + for g = gains + clpoles = pole(feedback(G_Kf{i}({"fm1", "fm2", "fm3", "fm4", "fm5", "fm6"}, {"f1", "f2", "f3", "f4", "f5", "f6"}), g*Kiff, +1)); + plot(real(clpoles), imag(clpoles), '.', 'color', colors(i,:), ... + 'HandleVisibility', 'off'); + end + +end + +xline(0, 'HandleVisibility', 'off'); yline(0, 'HandleVisibility', 'off'); +hold off; +axis equal; +xlim(1.1*[-900, 100]); ylim(1.1*[-100, 900]); +xticks(1.1*[-900:100:0]); +yticks(1.1*[0:100:900]); +set(gca, 'XTickLabel',[]); set(gca, 'YTickLabel',[]); +xlabel('Real part'); ylabel('Imaginary part'); +leg = legend('location', 'northwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; + +%% Identify the dynamics for several considered bending stiffnesses - APA300ML +G_Kf_apa300ml = {zeros(length(Kf), 1)}; + +for i = 1:length(Kf) + % Limited joint axial compliance + initializeSimplifiedNanoHexapod('actuator_type', 'apa300ml', ... + 'flex_type_F', '2dof', ... + 'flex_type_M', '3dof', ... + 'Fsm', 56e-3, ... % APA300ML weight 112g + 'Msm', 56e-3, ... + 'Kf_F', Kf(i), ... + 'Kf_M', Kf(i)); + + G_Kf_apa300ml(i) = {linearize(mdl, io)}; + G_Kf_apa300ml{i}.InputName = {'f1', 'f2', 'f3', 'f4', 'f5', 'f6'}; + G_Kf_apa300ml{i}.OutputName = {'l1', 'l2', 'l3', 'l4', 'l5', 'l6', 'fm1', 'fm2', 'fm3', 'fm4', 'fm5', 'fm6'}; +end + +Kiff = -1000 * ... % Gain + 1/(s) * ... % LPF: provides integral action + eye(6); % Diagonal 6x6 controller (i.e. decentralized) + +Kiff.InputName = {'fm1', 'fm2', 'fm3', 'fm4', 'fm5', 'fm6'}; +Kiff.OutputName = {'f1', 'f2', 'f3', 'f4', 'f5', 'f6'}; + +%% Root Locus for decentralized IFF - APA300ML actuator - Effect of joint bending stiffness +gains = logspace(-1, 2, 300); + +figure; +tiledlayout(1, 1, 'TileSpacing', 'compact', 'Padding', 'None'); +nexttile(); +hold on; + +for i = 1:length(Kf) + plot(real(pole(G_Kf_apa300ml{i}({"fm1", "fm2", "fm3", "fm4", "fm5", "fm6"}, {"f1", "f2", "f3", "f4", "f5", "f6"}))), imag(pole(G_Kf_apa300ml{i}({"fm1", "fm2", "fm3", "fm4", "fm5", "fm6"}, {"f1", "f2", "f3", "f4", "f5", "f6"}))), 'x', 'color', colors(i,:), ... + 'DisplayName', sprintf('$k_f = %.0f$ [Nm/rad]', Kf(i))); + plot(real(tzero(G_Kf_apa300ml{i}({"fm1", "fm2", "fm3", "fm4", "fm5", "fm6"}, {"f1", "f2", "f3", "f4", "f5", "f6"}))), imag(tzero(G_Kf_apa300ml{i}({"fm1", "fm2", "fm3", "fm4", "fm5", "fm6"}, {"f1", "f2", "f3", "f4", "f5", "f6"}))), 'o', 'color', colors(i,:), ... + 'HandleVisibility', 'off'); + + for g = gains + clpoles = pole(feedback(G_Kf_apa300ml{i}({"fm1", "fm2", "fm3", "fm4", "fm5", "fm6"}, {"f1", "f2", "f3", "f4", "f5", "f6"}), g*Kiff, +1)); + plot(real(clpoles), imag(clpoles), '.', 'color', colors(i,:), ... + 'HandleVisibility', 'off'); + end + +end + +xline(0, 'HandleVisibility', 'off'); yline(0, 'HandleVisibility', 'off'); +hold off; +axis equal; +xlim(1.4*[-900, 100]); ylim(1.4*[-100, 900]); +xticks(1.4*[-900:100:0]); +yticks(1.4*[0:100:900]); +set(gca, 'XTickLabel',[]); set(gca, 'YTickLabel',[]); +xlabel('Real part'); ylabel('Imaginary part'); +leg = legend('location', 'northwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; + +%% Identify the dynamics for several considered axial stiffnesses +% Let's initialize all the stages with default parameters. +initializeGround('type', 'rigid'); +initializeGranite('type', 'rigid'); +initializeTy('type', 'rigid'); +initializeRy('type', 'rigid'); +initializeRz('type', 'rigid'); +initializeMicroHexapod('type', 'rigid'); +initializeSample('m', 50); + +initializeSimscapeConfiguration(); +initializeDisturbances('enable', false); +initializeLoggingConfiguration('log', 'none'); +initializeController('type', 'open-loop'); +initializeReferences(); + +% Input/Output definition +clear io; io_i = 1; +io(io_i) = linio([mdl, '/Controller'], 1, 'openinput'); io_i = io_i + 1; % Actuator Inputs +io(io_i) = linio([mdl, '/Tracking Error'], 1, 'openoutput', [], 'EdL'); io_i = io_i + 1; % Errors in the frame of the struts +io(io_i) = linio([mdl, '/NASS'], 3, 'openoutput', [], 'fn'); io_i = io_i + 1; % Force Sensors + +% Effect of bending stiffness +Ka = 1e6*[1000, 100, 10, 1]; % [Nm/rad] +G_Ka = {zeros(length(Ka), 1)}; + +for i = 1:length(Ka) + % Limited joint axial compliance + initializeSimplifiedNanoHexapod('actuator_type', '1dof', ... + 'flex_type_F', '2dof_axial', ... + 'flex_type_M', '4dof', ... + 'actuator_k', 1e6, ... + 'actuator_c', 1e1, ... + 'actuator_kp', 0, ... + 'actuator_cp', 0, ... + 'Fsm', 56e-3, ... % APA300ML weight 112g + 'Msm', 56e-3, ... + 'Ca_F', 1, ... + 'Ca_M', 1, ... + 'Ka_F', Ka(i), ... + 'Ka_M', Ka(i)); + + G_Ka(i) = {linearize(mdl, io)}; + G_Ka{i}.InputName = {'f1', 'f2', 'f3', 'f4', 'f5', 'f6'}; + G_Ka{i}.OutputName = {'l1', 'l2', 'l3', 'l4', 'l5', 'l6', 'fm1', 'fm2', 'fm3', 'fm4', 'fm5', 'fm6'}; +end + +freqs = logspace(1, 4, 1000); + +%% Effect of the flexible joint axial stiffness on the HAC-plant +figure; +tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +for i = 1:length(Ka) + for j = 1:5 + for k = j+1:6 + plot(freqs, abs(squeeze(freqresp(G_Ka{i}("l"+k,"f"+j), freqs, 'Hz'))), 'color', [colors(i,:), 0.1], ... + 'HandleVisibility', 'off'); + end + end +end +for i = 1:length(Ka) + plot(freqs, abs(squeeze(freqresp(G_Ka{i}("l1","f1"), freqs, 'Hz'))), 'color', colors(i,:), 'DisplayName', sprintf('$k_a = %.0f$ [N/$\\mu$m]', 1e-6*Ka(i))); + % for j = 2:6 + % plot(freqs, abs(squeeze(freqresp(G_Ka{i}("l"+j,"f"+j), freqs, 'Hz'))), 'color', colors(i,:), 'HandleVisibility', 'off'); + % end +end +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]); +leg = legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; +ylim([1e-10, 1e-4]); + +ax2 = nexttile; +hold on; +for i = 1:length(Ka) + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_Ka{i}(1, 1), freqs, 'Hz'))))); +end +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +ylabel('Phase [deg]'); xlabel('Frequency [Hz]'); +ylim([-200, 20]); +yticks([-360:45:360]); + +linkaxes([ax1,ax2],'x'); + +%% Effect of the flexible joint axial stiffness on the IFF plant +figure; +tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +for i = 1:length(Ka) + for j = 1:5 + for k = j+1:6 + plot(freqs, abs(squeeze(freqresp(G_Ka{i}("fm"+k,"f"+j), freqs, 'Hz'))), 'color', [colors(i,:), 0.1], ... + 'HandleVisibility', 'off'); + end + end +end +for i = 1:length(Ka) + plot(freqs, abs(squeeze(freqresp(G_Ka{i}("fm1","f1"), freqs, 'Hz'))), 'color', colors(i,:), 'DisplayName', sprintf('$k_a = %.0f$ [N/$\\mu$m]', 1e-6*Ka(i))); + % for j = 2:6 + % plot(freqs, abs(squeeze(freqresp(G_Ka{i}("fm"+j,"f"+j), freqs, 'Hz'))), 'color', colors(i,:), 'HandleVisibility', 'off'); + % end +end +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude [N/N]'); set(gca, 'XTickLabel',[]); +leg = legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; +ylim([1e-4, 1e2]); + +ax2 = nexttile(); +hold on; +for i = 1:length(Ka) + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_Ka{i}("fm1", "f1"), freqs, 'Hz'))))); +end +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +ylabel('Phase [deg]'); xlabel('Frequency [Hz]'); +ylim([-20, 200]); +yticks([-360:45:360]); + +linkaxes([ax1,ax2],'x'); + +%% Decentalized IFF +Kiff = -200 * ... % Gain + 1/s * ... % LPF: provides integral action + eye(6); % Diagonal 6x6 controller (i.e. decentralized) + +Kiff.InputName = {'fm1', 'fm2', 'fm3', 'fm4', 'fm5', 'fm6'}; +Kiff.OutputName = {'f1', 'f2', 'f3', 'f4', 'f5', 'f6'}; + +%% Root Locus for decentralized IFF - 1dof actuator - Effect of joint bending stiffness +gains = logspace(-1, 2, 400); + +figure; +tiledlayout(1, 1, 'TileSpacing', 'compact', 'Padding', 'None'); +nexttile(); +hold on; + +for i = 1:length(Ka) + plot(real(pole(G_Ka{i}({"fm1", "fm2", "fm3", "fm4", "fm5", "fm6"}, {"f1", "f2", "f3", "f4", "f5", "f6"}))), imag(pole(G_Ka{i}({"fm1", "fm2", "fm3", "fm4", "fm5", "fm6"}, {"f1", "f2", "f3", "f4", "f5", "f6"}))), 'x', 'color', colors(i,:), ... + 'DisplayName', sprintf('$k_a = %.0f$ N/$\\mu$m', 1e-6*Ka(i))); + plot(real(tzero(G_Ka{i}({"fm1", "fm2", "fm3", "fm4", "fm5", "fm6"}, {"f1", "f2", "f3", "f4", "f5", "f6"}))), imag(tzero(G_Ka{i}({"fm1", "fm2", "fm3", "fm4", "fm5", "fm6"}, {"f1", "f2", "f3", "f4", "f5", "f6"}))), 'o', 'color', colors(i,:), ... + 'HandleVisibility', 'off'); + + for g = gains + clpoles = pole(feedback(G_Ka{i}({"fm1", "fm2", "fm3", "fm4", "fm5", "fm6"}, {"f1", "f2", "f3", "f4", "f5", "f6"}), g*Kiff, +1)); + plot(real(clpoles), imag(clpoles), '.', 'color', colors(i,:), ... + 'HandleVisibility', 'off'); + end + +end + +xline(0, 'HandleVisibility', 'off'); yline(0, 'HandleVisibility', 'off'); +hold off; +axis equal; +xlim(1.1*[-900, 100]); ylim(1.1*[-100, 900]); +xticks(1.1*[-900:100:0]); +yticks(1.1*[0:100:900]); +set(gca, 'XTickLabel',[]); set(gca, 'YTickLabel',[]); +xlabel('Real part'); ylabel('Imaginary part'); +leg = legend('location', 'northwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; + +%% Compute the damped plants +Kiff = -500 * ... % Gain + 1/(s + 2*pi*0.1) * ... % LPF: provides integral action + eye(6); % Diagonal 6x6 controller (i.e. decentralized) + +Kiff.InputName = {'fm1', 'fm2', 'fm3', 'fm4', 'fm5', 'fm6'}; +Kiff.OutputName = {'u1iff', 'u2iff', 'u3iff', 'u4iff', 'u5iff', 'u6iff'}; + +% New damped plant input +S1 = sumblk("f1 = u1iff + u1"); +S2 = sumblk("f2 = u2iff + u2"); +S3 = sumblk("f3 = u3iff + u3"); +S4 = sumblk("f4 = u4iff + u4"); +S5 = sumblk("f5 = u5iff + u5"); +S6 = sumblk("f6 = u6iff + u6"); + +G_Ka_iff = {zeros(1,length(Ka))}; +for i=1:length(Ka) + G_Ka_iff(i) = {connect(G_Ka{i}, Kiff, S1, S2, S3, S4, S5, S6, {'u1', 'u2', 'u3', 'u4', 'u5', 'u6'}, {'l1', 'l2', 'l3', 'l4', 'l5', 'l6'})}; +end + +%% Interaction Analysis - RGA Number +rga = zeros(length(Ka), length(freqs)); +for i = 1:length(Ka) + for j = 1:length(freqs) + rga(i,j) = sum(sum(abs(inv(evalfr(G_Ka_iff{i}({"l1", "l2", "l3", "l4", "l5", "l6"}, {"u1", "u2", "u3", "u4", "u5", "u6"}), 1j*2*pi*freqs(j)).').*evalfr(G_Ka_iff{i}({"l1", "l2", "l3", "l4", "l5", "l6"}, {"u1", "u2", "u3", "u4", "u5", "u6"}), 1j*2*pi*freqs(j)) - eye(6)))); + end +end + +%% RGA number for the damped plants - Effect of the flexible joint axial stiffness +figure; +hold on; +for i = 1:length(Ka) + plot(freqs, rga(i,:), 'DisplayName', sprintf('$k_a = %.0f$ N/$\\mu$m', 1e-6*Ka(i))) +end +hold off; +xlabel('Frequency [Hz]'); ylabel('RGA number'); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +ylim([0, 10]); xlim([10, 5e3]); +leg = legend('location', 'northwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; + +%% Extract stiffness of the joint from the reduced order model +% We first extract the stiffness and mass matrices. +K = readmatrix('flex025_mat_K.CSV'); +M = readmatrix('flex025_mat_M.CSV'); +% Then, we extract the coordinates of the interface nodes. +[int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('flex025_out_nodes_3D.txt'); + +m = 1; + +%% Name of the Simulink File +mdl = 'detail_fem_joint'; + +%% Input/Output definition +clear io; io_i = 1; +io(io_i) = linio([mdl, '/T'], 1, 'openinput'); io_i = io_i + 1; % Forces and Torques +io(io_i) = linio([mdl, '/D'], 1, 'openoutput'); io_i = io_i + 1; % Translations and Rotations + +G = linearize(mdl, io); + +% Stiffness extracted from the Simscape model +k_a = 1/dcgain(G(3,3)); % Axial stiffness [N/m] +k_f = 1/dcgain(G(4,4)); % Bending stiffness [N/m] +k_t = 1/dcgain(G(6,6)); % Torsion stiffness [N/m] + + +% Stiffness extracted from the Stiffness matrix +k_s = K(1,1); % shear [N/m] +% k_s = K(2,2); % shear [N/m] +k_a = K(3,3); % axial [N/m] +k_f = K(4,4); % bending [Nm/rad] +% k_f = K(5,5); % bending [Nm/rad] +k_t = K(6,6); % torsion [Nm/rad] + +%% Compare Dynamics between "Reduced Order" flexible joints and "2-dof and 3-dof" joints +% Let's initialize all the stages with default parameters. +initializeGround('type', 'rigid'); +initializeGranite('type', 'rigid'); +initializeTy('type', 'rigid'); +initializeRy('type', 'rigid'); +initializeRz('type', 'rigid'); +initializeMicroHexapod('type', 'rigid'); +initializeSample('m', 50); + +initializeSimscapeConfiguration(); +initializeDisturbances('enable', false); +initializeLoggingConfiguration('log', 'none'); +initializeController('type', 'open-loop'); +initializeReferences(); + +mdl = 'detail_fem_nass'; + +% Input/Output definition +clear io; io_i = 1; +io(io_i) = linio([mdl, '/Controller'], 1, 'openinput'); io_i = io_i + 1; % Actuator Inputs +io(io_i) = linio([mdl, '/Tracking Error'], 1, 'openoutput', [], 'EdL'); io_i = io_i + 1; % Errors in the frame of the struts +io(io_i) = linio([mdl, '/NASS'], 3, 'openoutput', [], 'fn'); io_i = io_i + 1; % Force Sensors + +% Fully flexible joints +initializeSimplifiedNanoHexapod('actuator_type', 'apa300ml', ... + 'flex_type_F', 'flexible', ... + 'flex_type_M', 'flexible', ... + 'Fsm', 56e-3, ... % APA300ML weight 112g + 'Msm', 56e-3); + +G_flex = linearize(mdl, io); +G_flex.InputName = {'f1', 'f2', 'f3', 'f4', 'f5', 'f6'}; +G_flex.OutputName = {'l1', 'l2', 'l3', 'l4', 'l5', 'l6', 'fm1', 'fm2', 'fm3', 'fm4', 'fm5', 'fm6'}; + +% Flexible joints modelled by 2DoF and 3DoF joints +initializeSimplifiedNanoHexapod('actuator_type', 'apa300ml', ... + 'flex_type_F', '2dof_axial', ... + 'flex_type_M', '4dof', ... + 'Kf_F', k_f, ... + 'Kt_F', k_t, ... + 'Ka_F', k_a, ... + 'Kf_M', k_f, ... + 'Kt_M', k_t, ... + 'Ka_M', k_a, ... + 'Cf_F', 1e-2, ... + 'Ct_F', 1e-2, ... + 'Ca_F', 1e-2, ... + 'Cf_M', 1e-2, ... + 'Ct_M', 1e-2, ... + 'Ca_M', 1e-2, ... + 'Fsm', 56e-3, ... % APA300ML weight 112g + 'Msm', 56e-3); + +G_ideal = linearize(mdl, io); +G_ideal.InputName = {'f1', 'f2', 'f3', 'f4', 'f5', 'f6'}; +G_ideal.OutputName = {'l1', 'l2', 'l3', 'l4', 'l5', 'l6', 'fm1', 'fm2', 'fm3', 'fm4', 'fm5', 'fm6'}; + +%% Comparison of the dynamics with joints modelled with FEM and modelled with "ideal joints" - HAC plant +freqs = logspace(1, 4, 1000); + +%% Effect of the flexible joint axial stiffness on the HAC-plant +figure; +tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +for j = 1:5 + for k = j+1:6 + plot(freqs, abs(squeeze(freqresp(G_flex("l"+k,"f"+j), freqs, 'Hz'))), 'color', [colors(1,:), 0.1], ... + 'HandleVisibility', 'off'); + plot(freqs, abs(squeeze(freqresp(G_ideal("l"+k,"f"+j), freqs, 'Hz'))), 'color', [colors(2,:), 0.1], ... + 'HandleVisibility', 'off'); + end +end +plot(freqs, abs(squeeze(freqresp(G_flex("l1","f1"), freqs, 'Hz'))), 'color', colors(1,:), 'DisplayName', 'Reduced Order Flexible Joints'); +plot(freqs, abs(squeeze(freqresp(G_ideal("l1","f1"), freqs, 'Hz'))), 'color', colors(2,:), 'DisplayName', 'Bot: $k_f$, $k_a$, Top: $k_f$, $k_t$, $k_a$'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]); +leg = legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; +ylim([1e-10, 1e-4]); + +ax2 = nexttile; +hold on; +plot(freqs, 180/pi*angle(squeeze(freqresp(G_flex("l1","f1"), freqs, 'Hz')))); +plot(freqs, 180/pi*angle(squeeze(freqresp(G_ideal("l1","f1"), freqs, 'Hz')))); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +ylabel('Phase [deg]'); xlabel('Frequency [Hz]'); +ylim([-20, 200]); +yticks([-360:45:360]); + +linkaxes([ax1,ax2],'x'); + +freqs = logspace(0, 3, 1000); + +%% Effect of the flexible joint axial stiffness on the HAC-plant +figure; +tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +for j = 1:5 + for k = j+1:6 + plot(freqs, abs(squeeze(freqresp(G_flex("fm"+k,"f"+j), freqs, 'Hz'))), 'color', [colors(1,:), 0.1], ... + 'HandleVisibility', 'off'); + plot(freqs, abs(squeeze(freqresp(G_ideal("fm"+k,"f"+j), freqs, 'Hz'))), 'color', [colors(2,:), 0.1], ... + 'HandleVisibility', 'off'); + end +end +plot(freqs, abs(squeeze(freqresp(G_flex("fm1","f1"), freqs, 'Hz'))), 'color', colors(1,:), 'DisplayName', 'Reduced Order Flexible Joints'); +plot(freqs, abs(squeeze(freqresp(G_ideal("fm1","f1"), freqs, 'Hz'))), 'color', colors(2,:), 'DisplayName', 'Bot: $k_f$, $k_a$, Top: $k_f$, $k_t$, $k_a$'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude [N/N]'); set(gca, 'XTickLabel',[]); +leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; +ylim([1e-5, 1e1]); + +ax2 = nexttile; +hold on; +plot(freqs, 180/pi*angle(squeeze(freqresp(G_flex("fm1","f1"), freqs, 'Hz')))); +plot(freqs, 180/pi*angle(squeeze(freqresp(G_ideal("fm1","f1"), freqs, 'Hz')))); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +ylabel('Phase [deg]'); xlabel('Frequency [Hz]'); +ylim([-20, 200]); +yticks([-360:45:360]); + +linkaxes([ax1,ax2],'x'); diff --git a/matlab/detail_fem_APA300ML.slx b/matlab/detail_fem_APA300ML.slx index 3fd8111..1262ef9 100644 Binary files a/matlab/detail_fem_APA300ML.slx and b/matlab/detail_fem_APA300ML.slx differ diff --git a/matlab/detail_fem_super_element.slx b/matlab/detail_fem_super_element.slx index 44c1df4..f2ca5b9 100644 Binary files a/matlab/detail_fem_super_element.slx and b/matlab/detail_fem_super_element.slx differ diff --git a/matlab/mat/detail_fem_apa95ml_open_circuit.mat b/matlab/mat/detail_fem_apa95ml_open_circuit.mat new file mode 100644 index 0000000..f0910bb Binary files /dev/null and b/matlab/mat/detail_fem_apa95ml_open_circuit.mat differ diff --git a/matlab/mat/detail_fem_apa95ml_short_circuit.mat b/matlab/mat/detail_fem_apa95ml_short_circuit.mat new file mode 100644 index 0000000..6cacaef Binary files /dev/null and b/matlab/mat/detail_fem_apa95ml_short_circuit.mat differ diff --git a/matlab/mat/nass_model_controller.mat b/matlab/mat/nass_model_controller.mat index 9c2b12b..49d2b09 100644 Binary files a/matlab/mat/nass_model_controller.mat and b/matlab/mat/nass_model_controller.mat differ diff --git a/matlab/mat/nass_model_disturbances.mat b/matlab/mat/nass_model_disturbances.mat index edf5852..9d9104b 100644 Binary files a/matlab/mat/nass_model_disturbances.mat and b/matlab/mat/nass_model_disturbances.mat differ diff --git a/matlab/mat/nass_model_stages.mat b/matlab/mat/nass_model_stages.mat index dee7f0e..5c3890d 100644 Binary files a/matlab/mat/nass_model_stages.mat and b/matlab/mat/nass_model_stages.mat differ diff --git a/nass-fem.bib b/nass-fem.bib index c9045d9..cd54d16 100644 --- a/nass-fem.bib +++ b/nass-fem.bib @@ -1,15 +1,15 @@ -@article{souleille18_concep_activ_mount_space_applic, - author = {Souleille, Adrien and Lampert, Thibault and Lafarga, V and - Hellegouarch, Sylvain and Rondineau, Alan and Rodrigues, - Gon{\c{c}}alo and Collette, Christophe}, - title = {A Concept of Active Mount for Space Applications}, - journal = {CEAS Space Journal}, - volume = 10, - number = 2, - pages = {157--165}, - year = 2018, - publisher = {Springer}, - keywords = {parallel robot, iff}, +@article{mcinroy02_model_desig_flexur_joint_stewar, + author = {J.E. McInroy}, + title = {Modeling and Design of Flexure Jointed Stewart Platforms + for Control Purposes}, + journal = {IEEE/ASME Transactions on Mechatronics}, + volume = 7, + number = 1, + pages = {95-99}, + year = 2002, + doi = {10.1109/3516.990892}, + url = {https://doi.org/10.1109/3516.990892}, + keywords = {parallel robot, flexure}, } @@ -109,6 +109,19 @@ +@book{pintelon12_system_ident, + author = {Rik Pintelon and Johan Schoukens}, + title = {System Identification : a Frequency Domain Approach}, + year = 2012, + publisher = {Wiley IEEE Press}, + url = {https://doi.org/10.1002/9781118287422}, + address = {Hoboken, N.J. Piscataway, NJ}, + doi = {10.1002/9781118287422}, + isbn = 9780470640371, +} + + + @phdthesis{hanieh03_activ_stewar, author = {Hanieh, Ahmed Abu}, keywords = {parallel robot}, @@ -120,18 +133,29 @@ -@article{mcinroy02_model_desig_flexur_joint_stewar, - author = {J.E. McInroy}, - title = {Modeling and Design of Flexure Jointed Stewart Platforms - for Control Purposes}, - journal = {IEEE/ASME Transactions on Mechatronics}, - volume = 7, - number = 1, - pages = {95-99}, - year = 2002, - doi = {10.1109/3516.990892}, - url = {https://doi.org/10.1109/3516.990892}, - keywords = {parallel robot, flexure}, +@article{souleille18_concep_activ_mount_space_applic, + author = {Souleille, Adrien and Lampert, Thibault and Lafarga, V and + Hellegouarch, Sylvain and Rondineau, Alan and Rodrigues, + Gon{\c{c}}alo and Collette, Christophe}, + title = {A Concept of Active Mount for Space Applications}, + journal = {CEAS Space Journal}, + volume = 10, + number = 2, + pages = {157--165}, + year = 2018, + publisher = {Springer}, + keywords = {parallel robot, iff}, +} + + + +@book{schmidt20_desig_high_perfor_mechat_third_revis_edition, + author = {Schmidt, R Munnig and Schitter, Georg and Rankers, Adrian}, + title = {The Design of High Performance Mechatronics - Third Revised + Edition}, + year = 2020, + publisher = {Ios Press}, + keywords = {favorite}, } @@ -187,3 +211,16 @@ keywords = {parallel robot}, } + + +@book{preumont18_vibrat_contr_activ_struc_fourt_edition, + author = {Andre Preumont}, + title = {Vibration Control of Active Structures - Fourth Edition}, + year = 2018, + publisher = {Springer International Publishing}, + url = {https://doi.org/10.1007/978-3-319-72296-2}, + doi = {10.1007/978-3-319-72296-2}, + keywords = {favorite, parallel robot}, + series = {Solid Mechanics and Its Applications}, +} + diff --git a/nass-fem.org b/nass-fem.org index 4e03421..17a2249 100644 --- a/nass-fem.org +++ b/nass-fem.org @@ -22,7 +22,7 @@ #+BIND: org-latex-bib-compiler "biber" #+PROPERTY: header-args:matlab :session *MATLAB* -#+PROPERTY: header-args:matlab+ :comments org +#+PROPERTY: header-args:matlab+ :comments none #+PROPERTY: header-args:matlab+ :exports none #+PROPERTY: header-args:matlab+ :results none #+PROPERTY: header-args:matlab+ :eval no-export @@ -139,993 +139,6 @@ Outline: -** Not used -*** Complete Strut with Encoder -:PROPERTIES: -:header-args:matlab+: :tangle matlab/strut_encoder.m -:END: -<> -**** Introduction - -Now, the full nano-hexapod strut is modelled using Ansys. - -The 3D as well as the interface nodes are shown in Figure ref:fig:strut_encoder_points3. - -#+name: fig:strut_encoder_points3 -#+caption: Interface points -#+attr_latex: :width 1.0\linewidth -[[file:figs/strut_fem_nodes.jpg]] - -A side view is shown in Figure ref:fig:strut_encoder_nodes_side. - -#+name: fig:strut_encoder_nodes_side -#+caption: Interface points - Side view -#+attr_latex: :width 1.0\linewidth -[[file:figs/strut_fem_nodes_side.jpg]] - -The flexible joints used have a 0.25mm width size. - -**** 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/'); - addpath('matlab/strut_encoder/'); -#+end_src - -#+begin_src matlab :eval no - addpath('strut_encoder/'); -#+end_src - -#+begin_src matlab - open('strut_encoder.slx'); -#+end_src - -**** Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates -We first extract the stiffness and mass matrices. -#+begin_src matlab - K = readmatrix('strut_encoder_mat_K.CSV'); - M = readmatrix('strut_encoder_mat_M.CSV'); -#+end_src - -#+begin_src matlab :exports results :results value table replace :tangle no - data2orgtable(K(1:10, 1:10), {}, {}, ' %.1g '); -#+end_src - -#+name: tab:strut_mat_K -#+caption: First 10x10 elements of the Stiffness matrix -#+attr_latex: :environment tabularx :width \linewidth :align XXXXXXXXXX -#+attr_latex: :center t :booktabs t :font \tiny -#+RESULTS: -| 2000000 | 1000000 | -3000000 | -400 | 300 | 200 | -30 | 2000 | -10000 | 0.3 | -| 1000000 | 4000000 | -8000000 | -900 | 400 | -50 | -6000 | 10000 | -20000 | 3 | -| -3000000 | -8000000 | 20000000 | 2000 | -900 | 200 | -10000 | 20000 | -300000 | 7 | -| -400 | -900 | 2000 | 5 | -0.1 | 05 | 1 | -3 | 6 | -0007 | -| 300 | 400 | -900 | -0.1 | 5 | 04 | -0.1 | 0.5 | -3 | 0001 | -| 200 | -50 | 200 | 05 | 04 | 300 | 4 | -01 | -1 | 3e-05 | -| -30 | -6000 | -10000 | 1 | -0.1 | 4 | 3000000 | -1000000 | -2000000 | -300 | -| 2000 | 10000 | 20000 | -3 | 0.5 | -01 | -1000000 | 6000000 | 7000000 | 1000 | -| -10000 | -20000 | -300000 | 6 | -3 | -1 | -2000000 | 7000000 | 20000000 | 2000 | -| 0.3 | 3 | 7 | -0007 | 0001 | 3e-05 | -300 | 1000 | 2000 | 5 | - - -#+begin_src matlab :exports results :results value table replace :tangle no - data2orgtable(M(1:10, 1:10), {}, {}, ' %.1g '); -#+end_src - -#+name: tab:strut_mat_M -#+caption: First 10x10 elements of the Mass matrix -#+attr_latex: :environment tabularx :width \linewidth :align XXXXXXXXXX -#+attr_latex: :center t :booktabs t :font \tiny -#+RESULTS: -| 0.04 | -0.005 | 0.007 | 2e-06 | 0.0001 | -5e-07 | -1e-05 | -9e-07 | 8e-05 | -5e-10 | -| -0.005 | 0.03 | 0.02 | -0.0001 | 1e-06 | -3e-07 | 3e-05 | -0.0001 | 8e-05 | -3e-08 | -| 0.007 | 0.02 | 0.08 | -6e-06 | -5e-06 | -7e-07 | 4e-05 | -0.0001 | 0.0005 | -3e-08 | -| 2e-06 | -0.0001 | -6e-06 | 2e-06 | -4e-10 | 2e-11 | -8e-09 | 3e-08 | -2e-08 | 6e-12 | -| 0.0001 | 1e-06 | -5e-06 | -4e-10 | 3e-06 | 2e-10 | -3e-09 | 3e-09 | -7e-09 | 6e-13 | -| -5e-07 | -3e-07 | -7e-07 | 2e-11 | 2e-10 | 5e-07 | -2e-08 | 5e-09 | -5e-09 | 1e-12 | -| -1e-05 | 3e-05 | 4e-05 | -8e-09 | -3e-09 | -2e-08 | 0.04 | 0.004 | 0.003 | 1e-06 | -| -9e-07 | -0.0001 | -0.0001 | 3e-08 | 3e-09 | 5e-09 | 0.004 | 0.02 | -0.02 | 0.0001 | -| 8e-05 | 8e-05 | 0.0005 | -2e-08 | -7e-09 | -5e-09 | 0.003 | -0.02 | 0.08 | -5e-06 | -| -5e-10 | -3e-08 | -3e-08 | 6e-12 | 6e-13 | 1e-12 | 1e-06 | 0.0001 | -5e-06 | 2e-06 | - - -Then, we extract the coordinates of the interface nodes. -#+begin_src matlab - [int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('strut_encoder_out_nodes_3D.txt'); -#+end_src - -#+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 - -#+name: tab:parameters_fem_strut -#+caption: Some extracted parameters of the FEM -#+attr_latex: :environment tabularx :width 0.4\linewidth :align lc -#+attr_latex: :center t :booktabs t -#+RESULTS: -| Total number of Nodes | 8 | -| Number of interface Nodes | 8 | -| Number of Modes | 6 | -| Size of M and K matrices | 54 | - -#+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:coordinate_nodes_strut -#+caption: Coordinates of the interface nodes -#+attr_latex: :environment tabularx :width 0.6\linewidth :align ccccc -#+attr_latex: :center t :booktabs t -#+RESULTS: -| Node i | Node Number | x [m] | y [m] | z [m] | -|--------+-------------+---------+--------+----------| -| 1.0 | 504411.0 | 0.0 | 0.0 | 0.0405 | -| 2.0 | 504412.0 | 0.0 | 0.0 | -0.0405 | -| 3.0 | 504413.0 | -0.0325 | 0.0 | 0.0 | -| 4.0 | 504414.0 | -0.0125 | 0.0 | 0.0 | -| 5.0 | 504415.0 | -0.0075 | 0.0 | 0.0 | -| 6.0 | 504416.0 | 0.0325 | 0.0 | 0.0 | -| 7.0 | 504417.0 | 0.004 | 0.0145 | -0.00175 | -| 8.0 | 504418.0 | 0.004 | 0.0166 | -0.00175 | - -Using =K=, =M= and =int_xyz=, we can use the =Reduced Order Flexible Solid= simscape block. - -**** Piezoelectric parameters -Parameters for the APA300ML: - -#+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 - -#+begin_src matlab - na = 2; % Number of stacks used as actuator - ns = 1; % Number of stacks used as force sensor -#+end_src - -**** Identification of the Dynamics -#+begin_src matlab :exports none - m = 0.01; -#+end_src - -The dynamics is identified from the applied force to the measured relative displacement. -#+begin_src matlab :exports none - %% Name of the Simulink File - mdl = 'strut_encoder'; - - %% Input/Output definition - clear io; io_i = 1; - io(io_i) = linio([mdl, '/Fa'], 1, 'openinput'); io_i = io_i + 1; - io(io_i) = linio([mdl, '/L'], 1, 'openoutput'); io_i = io_i + 1; - - Gh = -linearize(mdl, io); -#+end_src - -The same dynamics is identified for a payload mass of 10Kg. -#+begin_src matlab - m = 10; -#+end_src - -#+begin_src matlab :exports none - Ghm = -linearize(mdl, io); -#+end_src - -#+begin_src matlab :exports none - freqs = logspace(0, 4, 5000); - - figure; - tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); - - ax1 = nexttile([2,1]); - hold on; - plot(freqs, abs(squeeze(freqresp(Gh, freqs, 'Hz'))), '-'); - plot(freqs, abs(squeeze(freqresp(Ghm, freqs, 'Hz'))), '-'); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - ylabel('Amplitude'); set(gca, 'XTickLabel',[]); - axis padded 'auto x' - hold off; - - ax2 = nexttile; - hold on; - plot(freqs, 180/pi*angle(squeeze(freqresp(Gh, freqs, 'Hz'))), '-', ... - 'DisplayName', '$m = 0kg$'); - plot(freqs, 180/pi*angle(squeeze(freqresp(Ghm, freqs, 'Hz'))), '-', ... - 'DisplayName', '$m = 10kg$'); - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); - yticks(-360:90:360); - ylim([-180 180]); - xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); - hold off; - axis padded 'auto x' -#+end_src - -#+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/dynamics_encoder_full_strut.pdf', 'width', 'wide', 'height', 'tall'); -#+end_src - -#+name: fig:dynamics_encoder_full_strut -#+caption: Dynamics from the force actuator to the measured motion by the encoder -#+RESULTS: -[[file:figs/dynamics_encoder_full_strut.png]] - -*** Identification of the stiffness properties of the APA300ML -**** APA Alone -#+begin_src matlab :exports none - m = 10; -#+end_src - -#+begin_src matlab :exports none - %% Name of the Simulink File - mdl = 'APA300ML_characterisation'; - - %% 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, '/D'], 1, 'openoutput'); io_i = io_i + 1; % Displacement/Rotation [m] - - G = linearize(mdl, io); -#+end_src - -#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*) - data2orgtable([1e-6./dcgain(G(1,1)), 1e-6./dcgain(G(2,2)), 1e-6./dcgain(G(3,3)), 1./dcgain(G(4,4)), 1./dcgain(G(5,5)), 1./dcgain(G(6,6))]', {'Kx [N/um]', 'Ky [N/um]', 'Kz [N/um]', 'Rx [Nm/rad]', 'Ry [Nm/rad]', 'Rz [Nm/rad]'}, {'*Caracteristics*', '*Value*'}, ' %.1f '); -#+end_src - -#+RESULTS: -| *Caracteristics* | *Value* | -|------------------+---------| -| Kx [N/um] | 0.8 | -| Ky [N/um] | 1.6 | -| Kz [N/um] | 1.8 | -| Rx [Nm/rad] | 71.4 | -| Ry [Nm/rad] | 148.2 | -| Rz [Nm/rad] | 4241.8 | - -**** See how the global stiffness is changing with the flexible joints -#+begin_src matlab - flex = load('./mat/flexor_ID16.mat', 'int_xyz', 'int_i', 'n_xyz', 'n_i', 'nodes', 'M', 'K'); -#+end_src - -#+begin_src matlab :exports none - %% Name of the Simulink File - mdl = 'APA300ML_characterisation_with_joints'; - - %% 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, '/D'], 1, 'openoutput'); io_i = io_i + 1; % Displacement/Rotation [m] - - Gf = linearize(mdl, io); -#+end_src - -#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*) - data2orgtable([1e-6./dcgain(Gf(1,1)), 1e-6./dcgain(Gf(2,2)), 1e-6./dcgain(Gf(3,3)), 1./dcgain(Gf(4,4)), 1./dcgain(Gf(5,5)), 1./dcgain(Gf(6,6))]', {'Kx [N/um]', 'Ky [N/um]', 'Kz [N/um]', 'Rx [Nm/rad]', 'Ry [Nm/rad]', 'Rz [Nm/rad]'}, {'*Caracteristic*', '*Value*'}, ' %.1f '); -#+end_src - -#+RESULTS: -| *Caracteristic* | *Value* | -|-----------------+---------| -| Kx [N/um] | 0.0 | -| Ky [N/um] | 0.0 | -| Kz [N/um] | 1.8 | -| Rx [Nm/rad] | 722.9 | -| Ry [Nm/rad] | 129.6 | -| Rz [Nm/rad] | 115.3 | - - -#+begin_src matlab - freqs = logspace(-2, 5, 1000); - - figure; - hold on; - plot(freqs, abs(squeeze(freqresp(G(2,2), freqs, 'Hz'))), '-', 'DisplayName', 'APA'); - plot(freqs, abs(squeeze(freqresp(Gf(2,2), freqs, 'Hz'))), '-', 'DisplayName', 'Flex'); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - xlabel('Frequency [Hz]'); ylabel('Amplitude [m/N]'); - hold off; - legend('location', 'northeast'); -#+end_src - -#+begin_src matlab - freqs = logspace(-2, 5, 1000); - - figure; - hold on; - plot(freqs, abs(squeeze(freqresp(G(3,3), freqs, 'Hz'))), '-', 'DisplayName', 'APA'); - plot(freqs, abs(squeeze(freqresp(Gf(3,3), freqs, 'Hz'))), '-', 'DisplayName', 'Flex'); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - xlabel('Frequency [Hz]'); ylabel('Amplitude [m/N]'); - hold off; - legend('location', 'northeast'); -#+end_src - -*** Effect of APA300ML in the flexibility of the leg -#+begin_src matlab :exports none - m = 10; -#+end_src - -#+begin_src matlab :exports none - %% Name of the Simulink File - mdl = 'APA300ML_flex_joints'; - - %% 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, '/D'], 1, 'openoutput'); io_i = io_i + 1; % Displacement/Rotation [m] -#+end_src - -#+begin_src matlab :exports none - G = linearize(mdl, io); -#+end_src - -#+begin_src matlab :exports none - Gf = linearize(mdl, io); -#+end_src - -#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*) - data2orgtable([[1e-6./dcgain(G(1,1)), 1e-6./dcgain(G(2,2)), 1e-6./dcgain(G(3,3)), 1./dcgain(G(4,4)), 1./dcgain(G(5,5)), 1./dcgain(G(6,6))]', [1e-6./dcgain(Gf(1,1)), 1e-6./dcgain(Gf(2,2)), 1e-6./dcgain(Gf(3,3)), 1./dcgain(Gf(4,4)), 1./dcgain(Gf(5,5)), 1./dcgain(Gf(6,6))]'], {'Kx [N/um]', 'Ky [N/um]', 'Kz [N/um]', 'Rx [Nm/rad]', 'Ry [Nm/rad]', 'Rz [Nm/rad]'}, {'*Caracteristic*', '*Rigid APA*', '*Flexible APA*'}, ' %.3f '); -#+end_src - - -#+RESULTS: -| *Caracteristic* | *Rigid APA* | *Flexible APA* | -|-----------------+-------------+----------------| -| Kx [N/um] | 0.018 | 0.019 | -| Ky [N/um] | 0.018 | 0.018 | -| Kz [N/um] | 60.0 | 2.647 | -| Rx [Nm/rad] | 16.705 | 557.682 | -| Ry [Nm/rad] | 16.535 | 185.939 | -| Rz [Nm/rad] | 118.0 | 114.803 | - -*** First Flexible Joint Geometry -:PROPERTIES: -:header-args:matlab+: :tangle matlab/flexor_ID16.m -:END: -<> -**** Introduction :ignore: -The studied flexor is shown in Figure ref:fig:flexor_id16_screenshot. - -The stiffness and mass matrices representing the dynamics of the flexor are exported from a FEM. -It is then imported into Simscape. - -A simplified model of the flexor is then developped. - -#+name: fig:flexor_id16_screenshot -#+caption: Flexor studied -#+attr_latex: :width 0.3\linewidth -[[file:figs/flexor_id16_screenshot.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) - <> -#+end_src - -#+begin_src matlab :exports none :results silent :noweb yes - <> -#+end_src - -#+begin_src matlab :tangle no - addpath('matlab/'); - addpath('matlab/flexor_ID16/'); -#+end_src - -#+begin_src matlab :eval no - addpath('flexor_ID16/'); -#+end_src - -#+begin_src matlab - open('flexor_ID16.slx'); -#+end_src - -**** Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates -We first extract the stiffness and mass matrices. -#+begin_src matlab - K = extractMatrix('mat_K_6modes_2MDoF.matrix'); - M = extractMatrix('mat_M_6modes_2MDoF.matrix'); -#+end_src - -#+begin_src matlab :exports results :results value table replace :tangle no - data2orgtable(K(1:10, 1:10), {}, {}, ' %.2e '); -#+end_src - -#+name: tab:APA300ML_mat_K_2MDoF -#+caption: First 10x10 elements of the Stiffness matrix -#+attr_latex: :environment tabularx :width \linewidth :align XXXXXXXXXX -#+attr_latex: :center t :booktabs t :font \tiny -#+RESULTS: -| 11200000 | 195 | 2220 | -0.719 | -265 | 1.59 | -11200000 | -213 | -2220 | 0.147 | | 195 | 11400000 | 1290 | -148 | -0.188 | 2.41 | -212 | -11400000 | -1290 | 148 | -| 2220 | 1290 | 119000000 | 1.31 | 1.49 | 1.79 | -2220 | -1290 | -119000000 | -1.31 | | | | | | | | | | | | -| -0.719 | -148 | 1.31 | 33 | 000488 | -000977 | 0.141 | 148 | -1.31 | -33 | | | | | | | | | | | | -| -265 | -0.188 | 1.49 | 000488 | 33 | 00293 | 266 | 0.154 | -1.49 | 00026 | | | | | | | | | | | | -| 1.59 | 2.41 | 1.79 | -000977 | 00293 | 236 | -1.32 | -2.55 | -1.79 | 000379 | | | | | | | | | | | | -| -11200000 | -212 | -2220 | 0.141 | 266 | -1.32 | 11400000 | 24600 | 1640 | 120 | | | | | | | | | | | | -| -213 | -11400000 | -1290 | 148 | 0.154 | -2.55 | 24600 | 11400000 | 1290 | -72 | | | | | | | | | | | | -| -2220 | -1290 | -119000000 | -1.31 | -1.49 | -1.79 | 1640 | 1290 | 119000000 | 1.32 | | | | | | | | | | | | -| 0.147 | 148 | -1.31 | -33 | 00026 | 000379 | 120 | -72 | 1.32 | 34.7 | | | | | | | | | | | | - -#+begin_src matlab :exports results :results value table replace :tangle no - data2orgtable(M(1:10, 1:10), {}, {}, ' %.1g '); -#+end_src - -#+name: tab:APA300ML_mat_M_2MDoF -#+caption: First 10x10 elements of the Mass matrix -#+attr_latex: :environment tabularx :width \linewidth :align XXXXXXXXXX -#+attr_latex: :center t :booktabs t :font \tiny -#+RESULTS: -| 0.02 | 1e-09 | -4e-08 | -1e-10 | 0.0002 | -3e-11 | 0.004 | 5e-08 | 7e-08 | 1e-10 | -| 1e-09 | 0.02 | -3e-07 | -0.0002 | -1e-10 | -2e-09 | 2e-08 | 0.004 | 3e-07 | 1e-05 | -| -4e-08 | -3e-07 | 0.02 | 7e-10 | -2e-09 | 1e-09 | 3e-07 | 7e-08 | 0.003 | 1e-09 | -| -1e-10 | -0.0002 | 7e-10 | 4e-06 | -1e-12 | -6e-13 | 2e-10 | -7e-06 | -8e-10 | -1e-09 | -| 0.0002 | -1e-10 | -2e-09 | -1e-12 | 3e-06 | 2e-13 | 9e-06 | 4e-11 | 2e-09 | -3e-13 | -| -3e-11 | -2e-09 | 1e-09 | -6e-13 | 2e-13 | 4e-07 | 8e-11 | 9e-10 | -1e-09 | 2e-12 | -| 0.004 | 2e-08 | 3e-07 | 2e-10 | 9e-06 | 8e-11 | 0.02 | -7e-08 | -3e-07 | -2e-10 | -| 5e-08 | 0.004 | 7e-08 | -7e-06 | 4e-11 | 9e-10 | -7e-08 | 0.01 | -4e-08 | 0.0002 | -| 7e-08 | 3e-07 | 0.003 | -8e-10 | 2e-09 | -1e-09 | -3e-07 | -4e-08 | 0.02 | -1e-09 | -| 1e-10 | 1e-05 | 1e-09 | -1e-09 | -3e-13 | 2e-12 | -2e-10 | 0.0002 | -1e-09 | 2e-06 | - -#+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 - -#+name: tab:parameters_fem_2MDoF -#+caption: Some extracted parameters of the FEM -#+attr_latex: :environment tabularx :width 0.4\linewidth :align lc -#+attr_latex: :center t :booktabs t -#+RESULTS: -| Total number of Nodes | 2 | -| Number of interface Nodes | 2 | -| Number of Modes | 6 | -| Size of M and K matrices | 18 | - -Then, we extract the coordinates of the interface nodes. -#+begin_src matlab - [int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('out_nodes_3D.txt'); -#+end_src - -#+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:coordinate_nodes_2MDoF -#+caption: Coordinates of the interface nodes -#+attr_latex: :environment tabularx :width 0.6\linewidth :align ccccc -#+attr_latex: :center t :booktabs t -#+RESULTS: -| Node i | Node Number | x [m] | y [m] | z [m] | -|--------+-------------+-------+-------+-------| -| 1.0 | 181278.0 | 0.0 | 0.0 | 0.0 | -| 2.0 | 181279.0 | 0.0 | 0.0 | -0.0 | - -Using =K=, =M= and =int_xyz=, we can use the =Reduced Order Flexible Solid= simscape block. - -**** Identification of the parameters using Simscape and looking at the Stiffness Matrix -The flexor is now imported into Simscape and its parameters are estimated using an identification. - -#+begin_src matlab :exports none - m = 1; -#+end_src - -The dynamics is identified from the applied force/torque to the measured displacement/rotation of the flexor. -#+begin_src matlab :exports none - %% Name of the Simulink File - mdl = 'flexor_ID16'; - - %% Input/Output definition - clear io; io_i = 1; - io(io_i) = linio([mdl, '/T'], 1, 'openinput'); io_i = io_i + 1; - io(io_i) = linio([mdl, '/D'], 1, 'openoutput'); io_i = io_i + 1; - - G = linearize(mdl, io); -#+end_src - -And we find the same parameters as the one estimated from the Stiffness matrix. - -#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*) - data2orgtable([1e-6*K(3,3), K(4,4), K(5,5), K(6,6) ; 1e-6./dcgain(G(3,3)), 1./dcgain(G(4,4)), 1./dcgain(G(5,5)), 1./dcgain(G(6,6))]', {'Axial Stiffness Dz [N/um]', 'Bending Stiffness Rx [Nm/rad]', 'Bending Stiffness Ry [Nm/rad]', 'Torsion Stiffness Rz [Nm/rad]'}, {'*Characteristic*', '*Value*', '*Identification*'}, ' %0.f '); -#+end_src - -#+name: tab:comp_identified_stiffnesses -#+caption: Comparison of identified and FEM stiffnesses -#+attr_latex: :environment tabularx :width 0.7\linewidth :align lcc -#+attr_latex: :center t :booktabs t -#+RESULTS: -| *Characteristic* | *Value* | *Identification* | -|-------------------------------+---------+------------------| -| Axial Stiffness Dz [N/um] | 119 | 119 | -| Bending Stiffness Rx [Nm/rad] | 33 | 33 | -| Bending Stiffness Ry [Nm/rad] | 33 | 33 | -| Torsion Stiffness Rz [Nm/rad] | 236 | 236 | - -**** Simpler Model -Let's now model the flexible joint with a "perfect" Bushing joint as shown in Figure ref:fig:flexible_joint_simscape. - -#+name: fig:flexible_joint_simscape -#+caption: Bushing Joint used to model the flexible joint -#+attr_latex: :width 0.3\linewidth -[[file:figs/flexible_joint_simscape.png]] - -The parameters of the Bushing joint (stiffnesses) are estimated from the Stiffness matrix that was computed from the FEM. -#+begin_src matlab - Kx = K(1,1); % [N/m] - Ky = K(2,2); % [N/m] - Kz = K(3,3); % [N/m] - Krx = K(4,4); % [Nm/rad] - Kry = K(5,5); % [Nm/rad] - Krz = K(6,6); % [Nm/rad] -#+end_src - -The dynamics from the applied force/torque to the measured displacement/rotation of the flexor is identified again for this simpler model. -#+begin_src matlab :exports none - %% Name of the Simulink File - mdl = 'flexor_ID16_simplified'; - - %% Input/Output definition - clear io; io_i = 1; - io(io_i) = linio([mdl, '/T'], 1, 'openinput'); io_i = io_i + 1; - io(io_i) = linio([mdl, '/D'], 1, 'openoutput'); io_i = io_i + 1; - - Gs = linearize(mdl, io); -#+end_src - -The two obtained dynamics are compared in Figure - -#+begin_src matlab :exports none - freqs = logspace(0, 5, 1000); - - figure; - tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None'); - - ax1 = nexttile; - hold on; - set(gca,'ColorOrderIndex',1) - plot(freqs, abs(squeeze(freqresp(G(1,1), freqs, 'Hz'))), '-', 'DisplayName', '$D_x/F_x$'); - set(gca,'ColorOrderIndex',1) - plot(freqs, abs(squeeze(freqresp(Gs(1,1), freqs, 'Hz'))), '--', 'HandleVisibility', 'off'); - set(gca,'ColorOrderIndex',2) - plot(freqs, abs(squeeze(freqresp(G(2,2), freqs, 'Hz'))), '-', 'DisplayName', '$D_y/F_y$'); - set(gca,'ColorOrderIndex',2) - plot(freqs, abs(squeeze(freqresp(Gs(2,2), freqs, 'Hz'))), '--', 'HandleVisibility', 'off'); - set(gca,'ColorOrderIndex',3) - plot(freqs, abs(squeeze(freqresp(G(3,3), freqs, 'Hz'))), '-', 'DisplayName', '$D_z/F_z$'); - set(gca,'ColorOrderIndex',3) - plot(freqs, abs(squeeze(freqresp(Gs(3,3), freqs, 'Hz'))), '--', 'HandleVisibility', 'off'); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - xlabel('Frequency [Hz]'); ylabel('Amplitude [m/N]'); - hold off; - legend('location', 'southwest'); - - ax2 = nexttile; - hold on; - set(gca,'ColorOrderIndex',1) - plot(freqs, abs(squeeze(freqresp(G(4,4), freqs, 'Hz'))), '-', 'DisplayName', '$R_x/M_x$'); - set(gca,'ColorOrderIndex',1) - plot(freqs, abs(squeeze(freqresp(Gs(4,4), freqs, 'Hz'))), '--', 'HandleVisibility', 'off'); - set(gca,'ColorOrderIndex',2) - plot(freqs, abs(squeeze(freqresp(G(5,5), freqs, 'Hz'))), '-', 'DisplayName', '$R_y/M_y$'); - set(gca,'ColorOrderIndex',2) - plot(freqs, abs(squeeze(freqresp(Gs(5,5), freqs, 'Hz'))), '--', 'HandleVisibility', 'off'); - set(gca,'ColorOrderIndex',3) - plot(freqs, abs(squeeze(freqresp(G(6,6), freqs, 'Hz'))), '-', 'DisplayName', '$R_z/M_z$'); - set(gca,'ColorOrderIndex',3) - plot(freqs, abs(squeeze(freqresp(Gs(6,6), freqs, 'Hz'))), '--', 'HandleVisibility', 'off'); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - xlabel('Frequency [Hz]'); ylabel('Amplitude [rad/Nm]'); - hold off; - legend('location', 'southwest'); -#+end_src - -#+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/flexor_ID16_compare_bushing_joint.pdf', 'width', 'wide', 'height', 'normal'); -#+end_src - -#+name: fig:flexor_ID16_compare_bushing_joint -#+caption: Comparison of the Joint compliance between the FEM model and the simpler model -#+RESULTS: -[[file:figs/flexor_ID16_compare_bushing_joint.png]] - -*** APA300ML - transfer functions - -**** Identification of the Dynamics from actuator to replace displacement -We first set the mass to be approximately zero. -#+begin_src matlab :exports none - m = 0.01; -#+end_src - -The dynamics is identified from the applied force to the measured relative displacement. -#+begin_src matlab :exports none - %% Name of the Simulink File - mdl = 'APA300ML'; - - %% Input/Output definition - clear io; io_i = 1; - io(io_i) = linio([mdl, '/F'], 1, 'openinput'); io_i = io_i + 1; - io(io_i) = linio([mdl, '/z'], 1, 'openoutput'); io_i = io_i + 1; - - Gh = -linearize(mdl, io); -#+end_src - -The same dynamics is identified for a payload mass of 10Kg. -#+begin_src matlab - m = 10; -#+end_src - -#+begin_src matlab :exports none - %% Name of the Simulink File - mdl = 'APA300ML'; - - %% Input/Output definition - clear io; io_i = 1; - io(io_i) = linio([mdl, '/F'], 1, 'openinput'); io_i = io_i + 1; - io(io_i) = linio([mdl, '/z'], 1, 'openoutput'); io_i = io_i + 1; - - Ghm = -linearize(mdl, io); -#+end_src - -#+begin_src matlab :exports none - freqs = logspace(0, 4, 5000); - - figure; - tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); - - ax1 = nexttile([2,1]); - hold on; - plot(freqs, abs(squeeze(freqresp(Gh, freqs, 'Hz'))), '-'); - plot(freqs, abs(squeeze(freqresp(Ghm, freqs, 'Hz'))), '-'); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - ylabel('Amplitude'); set(gca, 'XTickLabel',[]); - hold off; - - ax2 = nexttile; - hold on; - plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gh, freqs, 'Hz')))), '-', ... - 'DisplayName', '$m = 0kg$'); - plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Ghm, freqs, 'Hz')))), '-', ... - 'DisplayName', '$m = 10kg$'); - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); - yticks(-360:90:360); - ylim([-360 0]); - xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); - hold off; - linkaxes([ax1,ax2],'x'); - xlim([freqs(1), freqs(end)]); - legend('location', 'southwest'); -#+end_src - -#+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/apa300ml_plant_dynamics.pdf', 'width', 'wide', 'height', 'tall'); -#+end_src - -#+name: fig:apa300ml_plant_dynamics -#+caption: Transfer function from forces applied by the stack to the axial displacement of the APA -#+RESULTS: -[[file:figs/apa300ml_plant_dynamics.png]] - -The root locus corresponding to Direct Velocity Feedback with a mass of 10kg is shown in Figure ref:fig:apa300ml_dvf_root_locus. -#+begin_src matlab :exports none - figure; - - gains = logspace(0, 5, 500); - - hold on; - plot(real(pole(Ghm)), imag(pole(G)), 'kx'); - plot(real(tzero(Ghm)), imag(tzero(G)), 'ko'); - for k = 1:length(gains) - cl_poles = pole(feedback(Ghm, gains(k)*s)); - plot(real(cl_poles), imag(cl_poles), 'k.'); - end - hold off; - axis square; - xlim([-500, 10]); ylim([0, 510]); - - xlabel('Real Part'); ylabel('Imaginary Part'); -#+end_src - -#+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/apa300ml_dvf_root_locus.pdf', 'width', 'wide', 'height', 'tall'); -#+end_src - -#+name: fig:apa300ml_dvf_root_locus -#+caption: Root Locus for Direct Velocity Feedback -#+RESULTS: -[[file:figs/apa300ml_dvf_root_locus.png]] - -**** Identification of the Dynamics from actuator to force sensor -Let's use 2 stacks as a force sensor and 1 stack as force actuator. - -The transfer function from actuator voltage to sensor voltage is identified and shown in Figure ref:fig:apa300ml_iff_plant. -#+begin_src matlab :exports none - m = 10; -#+end_src - -#+begin_src matlab :exports none - %% Name of the Simulink File - mdl = 'APA300ML'; - - %% Input/Output definition - clear io; io_i = 1; - io(io_i) = linio([mdl, '/Va'], 1, 'openinput'); io_i = io_i + 1; - io(io_i) = linio([mdl, '/Vs'], 1, 'openoutput'); io_i = io_i + 1; - - Giff = -linearize(mdl, io); -#+end_src - -#+begin_src matlab :exports none - freqs = logspace(0, 4, 5000); - - figure; - tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); - - ax1 = nexttile([2,1]); - hold on; - plot(freqs, abs(squeeze(freqresp(Giff, freqs, 'Hz'))), '-'); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - ylabel('Amplitude'); set(gca, 'XTickLabel',[]); - hold off; - - ax2 = nexttile; - hold on; - plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Giff, freqs, 'Hz')))), '-'); - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); - yticks(-360:90:360); - ylim([-180 180]); - xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); - hold off; - linkaxes([ax1,ax2],'x'); - xlim([freqs(1), freqs(end)]); -#+end_src - -#+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/apa300ml_iff_plant.pdf', 'width', 'wide', 'height', 'tall'); -#+end_src - -#+name: fig:apa300ml_iff_plant -#+caption: Transfer function from actuator to force sensor -#+RESULTS: -[[file:figs/apa300ml_iff_plant.png]] - -For root locus corresponding to IFF is shown in Figure ref:fig:apa300ml_iff_root_locus. -#+begin_src matlab :exports none - figure; - - gains = logspace(0, 5, 500); - - hold on; - plot(real(pole(Giff)), imag(pole(Giff)), 'kx'); - plot(real(tzero(Giff)), imag(tzero(Giff)), 'ko'); - for k = 1:length(gains) - cl_poles = pole(feedback(Giff, gains(k)/s)); - plot(real(cl_poles), imag(cl_poles), 'k.'); - end - hold off; - axis square; - xlim([-500, 10]); ylim([0, 510]); - - xlabel('Real Part'); ylabel('Imaginary Part'); -#+end_src - -#+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/apa300ml_iff_root_locus.pdf', 'width', 'wide', 'height', 'tall'); -#+end_src - -#+name: fig:apa300ml_iff_root_locus -#+caption: Root Locus for IFF -#+RESULTS: -[[file:figs/apa300ml_iff_root_locus.png]] - -**** Integral Force Feedback -In this section, Integral Force Feedback control architecture is applied on the APA300ML. - -First, the plant (dynamics from voltage actuator to voltage sensor is identified). -#+begin_src matlab :exports none - Kiff = tf(0); -#+end_src - -The payload mass is set to 10kg. -#+begin_src matlab - m = 10; -#+end_src - -#+begin_src matlab :exports none - %% Name of the Simulink File - mdl = 'APA300ML_IFF'; - - %% Input/Output definition - clear io; io_i = 1; - io(io_i) = linio([mdl, '/w'], 1, 'openinput'); io_i = io_i + 1; - io(io_i) = linio([mdl, '/F'], 1, 'openinput'); io_i = io_i + 1; - io(io_i) = linio([mdl, '/Fd'], 1, 'openinput'); io_i = io_i + 1; - io(io_i) = linio([mdl, '/z'], 1, 'openoutput'); io_i = io_i + 1; - io(io_i) = linio([mdl, '/APA300ML'], 1, 'openoutput'); io_i = io_i + 1; - - G_ol = linearize(mdl, io); - G_ol.InputName = {'w', 'f', 'F'}; - G_ol.OutputName = {'x1', 'Fs'}; - - G = G_ol({'Fs'}, {'f'}); -#+end_src - -The obtained dynamics is shown in Figure ref:fig:piezo_amplified_iff_plant. - -#+begin_src matlab :exports none - freqs = logspace(1, 5, 1000); - - figure; - tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); - - ax1 = nexttile([2,1]); - hold on; - plot(freqs, abs(squeeze(freqresp(G, freqs, 'Hz')))); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - ylabel('Amplitude'); set(gca, 'XTickLabel',[]); - hold off; - - ax2 = nexttile; - hold on; - plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G, freqs, 'Hz'))))); - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); - yticks(-360:90:360); - ylim([-390 30]); - xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); - hold off; - linkaxes([ax1,ax2],'x'); - xlim([freqs(1), freqs(end)]); -#+end_src - -#+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/piezo_amplified_iff_plant.pdf', 'width', 'wide', 'height', 'tall'); -#+end_src - -#+name: fig:piezo_amplified_iff_plant -#+caption: IFF Plant -#+RESULTS: -[[file:figs/piezo_amplified_iff_plant.png]] - -The controller is defined below and the loop gain is shown in Figure ref:fig:piezo_amplified_iff_loop_gain. -#+begin_src matlab - Kiff = -1e3/s; -#+end_src - -#+begin_src matlab :exports none - freqs = logspace(1, 5, 1000); - - figure; - tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); - - ax1 = nexttile([2,1]); - hold on; - plot(freqs, abs(squeeze(freqresp(G*Kiff, freqs, 'Hz')))); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - ylabel('Amplitude'); set(gca, 'XTickLabel',[]); - hold off; - - ax2 = nexttile; - hold on; - plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G*Kiff, freqs, 'Hz'))))); - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); - yticks(-360:90:360); - ylim([-180 180]); - xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); - hold off; - linkaxes([ax1,ax2],'x'); - xlim([freqs(1), freqs(end)]); -#+end_src - -#+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/piezo_amplified_iff_loop_gain.pdf', 'width', 'wide', 'height', 'tall'); -#+end_src - -#+name: fig:piezo_amplified_iff_loop_gain -#+caption: IFF Loop Gain -#+RESULTS: -[[file:figs/piezo_amplified_iff_loop_gain.png]] - -Now the closed-loop system is identified again and compare with the open loop system in Figure ref:fig:piezo_amplified_iff_comp. - -It is the expected behavior as shown in the Figure ref:fig:souleille18_results (from cite:souleille18_concep_activ_mount_space_applic). - -#+begin_src matlab :exports none - clear io; io_i = 1; - io(io_i) = linio([mdl, '/w'], 1, 'openinput'); io_i = io_i + 1; - io(io_i) = linio([mdl, '/F'], 1, 'openinput'); io_i = io_i + 1; - io(io_i) = linio([mdl, '/Fd'], 1, 'openinput'); io_i = io_i + 1; - io(io_i) = linio([mdl, '/z'], 1, 'openoutput'); io_i = io_i + 1; - io(io_i) = linio([mdl, '/APA300ML'], 1, 'output'); io_i = io_i + 1; - - Giff = linearize(mdl, io); - Giff.InputName = {'w', 'f', 'F'}; - Giff.OutputName = {'x1', 'Fs'}; -#+end_src - -#+begin_src matlab :exports none - freqs = logspace(0, 3, 1000); - - figure; - tiledlayout(2, 3, 'TileSpacing', 'None', 'Padding', 'None'); - - ax1 = nexttile; - hold on; - plot(freqs, abs(squeeze(freqresp(G_ol('x1', 'w'), freqs, 'Hz')))); - plot(freqs, abs(squeeze(freqresp(Giff('x1', 'w'), freqs, 'Hz')))); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - set(gca, 'XTickLabel',[]); ylabel('$x_1/w$ [m/m]') - - ax2 = nexttile; - hold on; - plot(freqs, abs(squeeze(freqresp(G_ol('x1', 'f'), freqs, 'Hz')))); - plot(freqs, abs(squeeze(freqresp(Giff('x1', 'f'), freqs, 'Hz')))); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - set(gca, 'XTickLabel',[]); ylabel('$x_1/f$ [m/N]'); - - ax3 = nexttile; - hold on; - plot(freqs, abs(squeeze(freqresp(G_ol('x1', 'F'), freqs, 'Hz')))); - plot(freqs, abs(squeeze(freqresp(Giff('x1', 'F'), freqs, 'Hz')))); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - set(gca, 'XTickLabel',[]); ylabel('$x_1/F$ [m/N]'); - - ax4 = nexttile; - hold on; - plot(freqs, abs(squeeze(freqresp(G_ol('Fs', 'w'), freqs, 'Hz')))); - plot(freqs, abs(squeeze(freqresp(Giff('Fs', 'w'), freqs, 'Hz')))); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - xlabel('Frequency [Hz]'); ylabel('$F_s/w$ [N/m]'); - - ax5 = nexttile; - hold on; - plot(freqs, abs(squeeze(freqresp(G_ol('Fs', 'f'), freqs, 'Hz')))); - plot(freqs, abs(squeeze(freqresp(Giff('Fs', 'f'), freqs, 'Hz')))); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - xlabel('Frequency [Hz]'); ylabel('$F_s/f$ [N/N]'); - - ax6 = nexttile; - hold on; - plot(freqs, abs(squeeze(freqresp(G_ol('Fs', 'F'), freqs, 'Hz')))); - plot(freqs, abs(squeeze(freqresp(Giff('Fs', 'F'), freqs, 'Hz')))); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - xlabel('Frequency [Hz]'); ylabel('$F_s/F$ [N/N]'); - - linkaxes([ax1,ax2,ax3,ax4,ax5,ax6],'x'); -#+end_src - -#+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/piezo_amplified_iff_comp.pdf', 'width', 'full', 'height', 'full'); -#+end_src - -#+name: fig:piezo_amplified_iff_comp -#+caption: OL and CL transfer functions -#+RESULTS: -[[file:figs/piezo_amplified_iff_comp.png]] - -#+name: fig:souleille18_results -#+caption: Results obtained in cite:souleille18_concep_activ_mount_space_applic -#+attr_latex: :width 1.0\linewidth -[[file:figs/souleille18_results.png]] - ** DONE [#B] Add Matlab files related to APA95ML and compute all the figures CLOSED: [2025-02-21 Fri 15:09] SCHEDULED: <2025-02-21 Fri> @@ -1170,6 +183,9 @@ Section ref:sec:detail_fem_joint addresses the design of flexible joints, where In both cases, the hybrid modeling approach enables detailed component optimization while maintaining the ability to predict system-level dynamic behavior, particularly under closed-loop control conditions. * Reduced order flexible bodies +:PROPERTIES: +:HEADER-ARGS:matlab+: :tangle matlab/detail_fem_1_flexible_body.m +:END: <> ** Introduction :ignore: @@ -1291,7 +307,7 @@ Six additional modes were considered, resulting in total model order of $48$. The modal reduction procedure was then executed, yielding the reduced mass and stiffness matrices that form the foundation of the component's representation in the multi-body simulation environment. #+name: fig:detail_fem_apa95ml_model -#+caption: Obtained mesh and defined interface frames (or "remote points") in the finite element model of the APA95ML (\subref{fig:detail_fem_apa95ml_mesh}). Interface with the multi-body model is shown in (\subref{fig:detail_fem_apa_modal_schematic}). +#+caption: Obtained mesh and defined interface frames (or "remote points") in the finite element model of the APA95ML (\subref{fig:detail_fem_apa95ml_mesh}). Interface with the multi-body model is shown in (\subref{fig:detail_fem_apa_model_schematic}). #+attr_latex: :options [htbp] #+begin_figure #+attr_latex: :caption \subcaption{\label{fig:detail_fem_apa95ml_mesh}Obtained mesh and "remote points"} @@ -1335,7 +351,7 @@ From [[cite:&fleming10_integ_strain_force_feedb_high]] the relation between the F_a = g_a \cdot V_a, \quad g_a = d_{33} n k_a, \quad k_a = \frac{c^{E} A}{L} \end{equation} -Unfortunately, it is difficult to know exactly which material is used in the amplified piezoelectric actuator[fn:1]. +Unfortunately, it is difficult to know exactly which material is used in the amplified piezoelectric actuator[fn:detail_fem_1]. However, based on the available properties of the stacks in the data-sheet (summarized in Table ref:tab:detail_fem_stack_parameters), the soft Lead Zirconate Titanate "THP5H" from Thorlabs seemed to match quite well the observed properties. #+name: tab:detail_fem_stack_parameters @@ -1408,6 +424,8 @@ A value of $23\,N/\mu m$ was found which is close to the specified stiffness in #+begin_src matlab %% Stiffness estimation m = 0.0001; % block-free condition, no payload +k_support = 1e9; +c_support = 1e3; clear io; io_i = 1; io(io_i) = linio([mdl, '/Fd'], 1, 'openinput'); io_i = io_i + 1; @@ -1420,7 +438,7 @@ G = linearize(mdl, io); k_est = 1/dcgain(G); % [N/m] #+end_src -The multi-body model predicted a resonant frequency under block-free conditions of $2024\,\text{Hz}$ (Figure ref:fig:detail_fem_apa95ml_compliance), which is in agreement with the nominal specification of $2000\,\text{Hz}$. +The multi-body model predicted a resonant frequency under block-free conditions of $\approx 2\,\text{kHz}$ (Figure ref:fig:detail_fem_apa95ml_compliance), which is in agreement with the nominal specification of $2\,\text{kHz}$. #+begin_src matlab :exports none :results none %% Estimated compliance of the APA95ML @@ -1548,6 +566,8 @@ M = readmatrix('APA95ML_M.CSV'); [int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('APA95ML_out_nodes_3D.txt'); m = 5.5; % Mass of the suspended granite [kg] +k_support = 4e7; +c_support = 3e2; % Compute transfer functions clear io; io_i = 1; @@ -1562,6 +582,7 @@ Gm.OutputName = {'y', 'Vs'}; #+begin_src matlab :exports none :results none %% Comparison of the identified transfer function from Va to di to the multi-body model +freqs = logspace(1, 3, 500); figure; tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); @@ -1829,6 +850,9 @@ This is exemplified by the nano-hexapod configuration, where the implementation However, the methodology remains valuable for system analysis, as the extraction of frequency domain characteristics can be efficiently performed even with such high-order models. * Actuator Selection +:PROPERTIES: +:HEADER-ARGS:matlab+: :tangle matlab/detail_fem_2_actuators.m +:END: <> ** Introduction :ignore: @@ -1931,7 +955,7 @@ Moreover, using APA for active damping has been successfully demonstrated in sim Several specific APA models were evaluated against the established specifications (Table ref:tab:detail_fem_piezo_act_models). The APA300ML emerged as the optimal choice. -This selection was further reinforced by previous experience with APAs from the same manufacturer[fn:2], and particularly by the successful validation of the modeling methodology with a similar actuator (Section ref:ssec:detail_fem_super_element_example). +This selection was further reinforced by previous experience with APAs from the same manufacturer[fn:detail_fem_2], and particularly by the successful validation of the modeling methodology with a similar actuator (Section ref:ssec:detail_fem_super_element_example). The demonstrated accuracy of the modeling approach for the APA95ML provides confidence in the reliable prediction of the APA300ML's dynamic characteristics, thereby supporting both the selection decision and subsequent dynamical analyses. #+name: tab:detail_fem_piezo_act_models @@ -2027,6 +1051,10 @@ While higher-order modes and non-axial flexibility are not captured, the model a #+begin_src matlab %% Identify dynamics with "Reduced Order Flexible Body" +K = readmatrix('APA300ML_mat_K.CSV'); +M = readmatrix('APA300ML_mat_M.CSV'); +[int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('APA300ML_out_nodes_3D.txt'); + m = 5; % [kg] ga = 25.9; % [N/V] gs = 5.08e6; % [V/m] @@ -2420,6 +1448,9 @@ exportFig('figs/detail_fem_actuator_fem_vs_perfect_iff_plant.pdf', 'width', 'hal #+end_figure * Flexible Joint Design +:PROPERTIES: +:HEADER-ARGS:matlab+: :tangle matlab/detail_fem_3_flexible_joints.m +:END: <> ** Introduction :ignore: @@ -2617,13 +1648,13 @@ hold on; for i = 1:length(Kf) for j = 1:5 for k = j+1:6 - plot(freqs, abs(squeeze(freqresp(G_Kf{i}("fn"+k,"f"+j), freqs, 'Hz'))), 'color', [colors(i,:), 0.1], ... + plot(freqs, abs(squeeze(freqresp(G_Kf{i}("fm"+k,"f"+j), freqs, 'Hz'))), 'color', [colors(i,:), 0.1], ... 'HandleVisibility', 'off'); end end plot(freqs, abs(squeeze(freqresp(G_Kf{i}("fm1","f1"), freqs, 'Hz'))), 'color', colors(i,:), 'DisplayName', sprintf('$k_f = %.0f $ [Nm/rad]', Kf(i))); for j = 2:6 - plot(freqs, abs(squeeze(freqresp(G_Kf{i}("fn"+j,"f"+j), freqs, 'Hz'))), 'color', colors(i,:), 'HandleVisibility', 'off'); + plot(freqs, abs(squeeze(freqresp(G_Kf{i}("fm"+j,"f"+j), freqs, 'Hz'))), 'color', colors(i,:), 'HandleVisibility', 'off'); end end hold off; @@ -3820,5 +2851,5 @@ end #+end_src * Footnotes -[fn:2]Cedrat technologies -[fn:1]The manufacturer of the APA95ML was not willing to share the piezoelectric material properties of the stack. +[fn:detail_fem_2]Cedrat technologies +[fn:detail_fem_1]The manufacturer of the APA95ML was not willing to share the piezoelectric material properties of the stack. diff --git a/nass-fem.pdf b/nass-fem.pdf index 6ca9218..2174975 100644 Binary files a/nass-fem.pdf and b/nass-fem.pdf differ diff --git a/nass-fem.tex b/nass-fem.tex index be5d844..bee0bd9 100644 --- a/nass-fem.tex +++ b/nass-fem.tex @@ -1,4 +1,4 @@ -% Created 2025-02-27 Thu 01:24 +% Created 2025-02-27 Thu 10:38 % Intended LaTeX compiler: pdflatex \documentclass[a4paper, 10pt, DIV=12, parskip=full, bibliography=totoc]{scrreprt} @@ -38,7 +38,6 @@ Section \ref{sec:detail_fem_joint} addresses the design of flexible joints, wher In both cases, the hybrid modeling approach enables detailed component optimization while maintaining the ability to predict system-level dynamic behavior, particularly under closed-loop control conditions. \chapter{Reduced order flexible bodies} -\label{sec:orgfd22661} \label{sec:detail_fem_super_element} Components exhibiting complex dynamical behavior are frequently found to be unsuitable for direct implementation within multi-body models. These components are traditionally analyzed using Finite Element Analysis (FEA) software. @@ -50,7 +49,6 @@ First, the fundamental principles and methodological approaches of this modeling It is then illustrated through its practical application to the modelling of an Amplified Piezoelectric Actuator (APA) (Section \ref{ssec:detail_fem_super_element_example}). Finally, the validity of this modeling approach is demonstrated through experimental validation, wherein the obtained dynamics from the hybrid modelling approach is compared with measurements (Section \ref{ssec:detail_fem_super_element_validation}). \section{Procedure} -\label{sec:org93ab665} \label{ssec:detail_fem_super_element_theory} In this modeling approach, some components within the multi-body framework are represented as \emph{reduced-order flexible bodies}, wherein their modal behavior is characterized through reduced mass and stiffness matrices derived from finite element analysis (FEA) models. @@ -74,7 +72,6 @@ m = 6 \times n + p \end{equation} \section{Example with an Amplified Piezoelectric Actuator} -\label{sec:org1e66a5f} \label{ssec:detail_fem_super_element_example} The presented modeling framework was first applied to an Amplified Piezoelectric Actuator (APA) for several reasons. Primarily, this actuator represents an excellent candidate for implementation within the nano-hexapod, as will be elaborated in Section \ref{sec:detail_fem_actuator}. @@ -105,7 +102,6 @@ Stiffness & \(21\,N/\mu m\)\\ \captionof{table}{\label{tab:detail_fem_apa95ml_specs}APA95ML specifications} \end{minipage} \paragraph{Finite Element Model} -\label{sec:orgce5afdb} The development of the finite element model for the APA95ML necessitated the specification of appropriate material properties, as summarized in Table \ref{tab:detail_fem_material_properties}. The finite element mesh, shown in Figure \ref{fig:detail_fem_apa95ml_mesh}, was then generated. @@ -142,11 +138,10 @@ The modal reduction procedure was then executed, yielding the reduced mass and s \end{center} \subcaption{\label{fig:detail_fem_apa_model_schematic}Inclusion in multi-body model} \end{subfigure} -\caption{\label{fig:detail_fem_apa95ml_model}Obtained mesh and defined interface frames (or ``remote points'') in the finite element model of the APA95ML (\subref{fig:detail_fem_apa95ml_mesh}). Interface with the multi-body model is shown in (\subref{fig:detail_fem_apa_modal_schematic}).} +\caption{\label{fig:detail_fem_apa95ml_model}Obtained mesh and defined interface frames (or ``remote points'') in the finite element model of the APA95ML (\subref{fig:detail_fem_apa95ml_mesh}). Interface with the multi-body model is shown in (\subref{fig:detail_fem_apa_model_schematic}).} \end{figure} \paragraph{Super Element in the Multi-Body Model} -\label{sec:org809b5e7} Previously computed reduced order mass and stiffness matrices were imported in a multi-body model block called ``Reduced Order Flexible Solid''. This block has several interface frames corresponding to the ones defined in the FEA software. @@ -158,7 +153,6 @@ This is illustrated in Figure \ref{fig:detail_fem_apa_model_schematic}. However, to have access to the physical voltage input of the actuators stacks \(V_a\) and to the generated voltage by the force sensor \(V_s\), conversion between the electrical and mechanical domains need to be determined. \paragraph{Sensor and Actuator ``constants''} -\label{sec:orgb3a075f} To link the electrical domain to the mechanical domain, an ``actuator constant'' \(g_a\) and a ``sensor constant'' \(g_s\) were introduced as shown in Figure \ref{fig:detail_fem_apa_model_schematic}. @@ -217,7 +211,6 @@ From these parameters, \(g_s = 5.1\,V/\mu m\) and \(g_a = 26\,N/V\) were obtaine \end{table} \paragraph{Identification of the APA Characteristics} -\label{sec:orge041867} Initial validation of the finite element model and its integration as a reduced-order flexible model within the multi-body model was accomplished through comparative analysis of key actuator characteristics against manufacturer specifications. @@ -225,7 +218,7 @@ The stiffness of the APA95ML was estimated from the multi-body model by computin The inverse of the DC gain this transfer function corresponds to the axial stiffness of the APA95ML. A value of \(23\,N/\mu m\) was found which is close to the specified stiffness in the datasheet of \(k = 21\,N/\mu m\). -The multi-body model predicted a resonant frequency under block-free conditions of \(2024\,\text{Hz}\) (Figure \ref{fig:detail_fem_apa95ml_compliance}), which is in agreement with the nominal specification of \(2000\,\text{Hz}\). +The multi-body model predicted a resonant frequency under block-free conditions of \(\approx 2\,\text{kHz}\) (Figure \ref{fig:detail_fem_apa95ml_compliance}), which is in agreement with the nominal specification of \(2\,\text{kHz}\). \begin{figure}[htbp] \centering @@ -243,7 +236,6 @@ Through the established amplification factor of 1.5, this translates to a predic The high degree of concordance observed across multiple performance metrics provides a first validation of the ability to include FEM into multi-body model. \section{Experimental Validation} -\label{sec:org354cea4} \label{ssec:detail_fem_super_element_validation} Further validation of the reduced-order flexible body methodology was undertaken through experimental investigation. The goal is to measure the dynamics of the APA95ML and compared it with predictions derived from the multi-body model incorporating the actuator as a flexible element. @@ -269,7 +261,6 @@ Measurement of the sensor stack voltage \(V_s\) was performed using an analog-to \caption{\label{fig:detail_fem_apa95ml_bench}Test bench used to validate ``reduced order solid bodies'' using an APA95ML. Picture of the bench is shown in (\subref{fig:detail_fem_apa95ml_bench_picture}). Schematic is shown in (\subref{fig:detail_fem_apa95ml_bench_schematic}).} \end{figure} \paragraph{Comparison of the dynamics} -\label{sec:orgc0ac08b} Frequency domain system identification techniques were used to characterize the dynamic behavior of the APA95ML. The identification procedure necessitated careful choice of the excitation signal \cite[, chap. 5]{pintelon12_system_ident}. @@ -304,7 +295,6 @@ Regarding the amplitude characteristics, the constants \(g_a\) and \(g_s\) could \end{figure} \paragraph{Integral Force Feedback with APA} -\label{sec:org4b65d74} To further validate this modeling methodology, its ability to predict closed-loop behavior was verified experimentally. Integral Force Feedback (IFF) was implemented using the force sensor stack, and the measured dynamics of the damped system were compared with model predictions across multiple feedback gains. @@ -338,7 +328,6 @@ The close agreement between experimental measurements and theoretical prediction \end{figure} \section*{Conclusion} -\label{sec:org105aef7} The modeling procedure presented in this section will demonstrate significant utility for the optimization of complex mechanical components within multi-body systems, particularly in the design of actuators (Section \ref{sec:detail_fem_actuator}) and flexible joints (Section \ref{sec:detail_fem_joint}). Through experimental validation using an Amplified Piezoelectric Actuator, the methodology has been shown to accurately predict both open-loop and closed-loop dynamic behavior, thereby establishing its reliability for component design and system analysis. @@ -348,7 +337,6 @@ This is exemplified by the nano-hexapod configuration, where the implementation However, the methodology remains valuable for system analysis, as the extraction of frequency domain characteristics can be efficiently performed even with such high-order models. \chapter{Actuator Selection} -\label{sec:org3ec4809} \label{sec:detail_fem_actuator} The selection and modeling of actuators constitutes a critical step in the development of the nano-hexapod. This chapter presents the approach to actuator selection and modeling. @@ -356,7 +344,6 @@ First, specifications for the nano-hexapod actuators are derived from previous a Then, the chosen actuator is modeled using the reduced-order flexible body approach developed in the previous section, enabling validation of this selection through detailed dynamical analysis (Section \ref{ssec:detail_fem_actuator_apa300ml}). Finally, a simplified two-degree-of-freedom model is developed to facilitate time-domain simulations while maintaining accurate representation of the actuator's essential characteristics (Section \ref{ssec:detail_fem_actuator_apa300ml_2dof}). \section{Choice of the Actuator based on Specifications} -\label{sec:org929a34b} \label{ssec:detail_fem_actuator_specifications} The actuator selection process was driven by several critical requirements derived from previous dynamic analyses. @@ -431,7 +418,6 @@ Height \(< 50\, [mm]\) & 22 & 30 & 24 & 27 & 16\\ \end{table} \section{APA300ML - Reduced Order Flexible Body} -\label{sec:orgeff9b1b} \label{ssec:detail_fem_actuator_apa300ml} The validation of the APA300ML started by incorporating a ``reduced order flexible body'' into the multi-body model as explained in Section \ref{sec:detail_fem_super_element}. @@ -459,7 +445,6 @@ While this high order provides excellent accuracy for validation purposes, it pr The sensor and actuator ``constants'' (\(g_s\) and \(g_a\)) derived in Section \ref{ssec:detail_fem_super_element_example} for the APA95ML were used for the APA300ML model, as both actuators employ identical piezoelectric stacks. \section{Simpler 2DoF Model of the APA300ML} -\label{sec:org120c274} \label{ssec:detail_fem_actuator_apa300ml_2dof} To facilitate efficient time-domain simulations while maintaining essential dynamic characteristics, a simplified two-degree-of-freedom model was developed, adapted from \cite{souleille18_concep_activ_mount_space_applic}. @@ -532,7 +517,6 @@ While higher-order modes and non-axial flexibility are not captured, the model a \end{figure} \section{Electrical characteristics of the APA} -\label{sec:orga7af5a1} \label{ssec:detail_fem_actuator_apa300ml_electrical} The behavior of piezoelectric actuators is characterized by coupled constitutive equations that establish relationships between electrical properties (charges, voltages) and mechanical properties (stress, strain) \cite[, chapter 5.5]{schmidt20_desig_high_perfor_mechat_third_revis_edition}. @@ -553,7 +537,6 @@ Proper consideration must be given to voltage amplifier specifications and force These aspects, being fundamental to system implementation, will be addressed in the instrumentation chapter. \section{Validation with the Nano-Hexapod} -\label{sec:orgba951c9} \label{ssec:detail_fem_actuator_apa300ml_validation} The integration of the APA300ML model within the nano-hexapod simulation framework served two validation objectives: to validate the APA300ML choice through analysis of system dynamics with APA modelled as flexible bodies, and to validate the simplified 2DoF model through comparative analysis with the full FEM implementation. @@ -586,7 +569,6 @@ These results validate both the selection of the APA300ML and the effectiveness \end{figure} \chapter{Flexible Joint Design} -\label{sec:org93c9b2c} \label{sec:detail_fem_joint} High-precision position control at the nanometer scale requires systems to be free from friction and backlash, as these nonlinear phenomena severely limit achievable positioning accuracy. This fundamental requirement prevents the use of conventional joints, necessitating instead the implementation of flexible joints that achieve motion through elastic deformation. @@ -622,7 +604,6 @@ The analysis of bending and axial stiffness effects enables the establishment of These specifications guide the development and optimization of a flexible joint design through finite element analysis (Section \ref{ssec:detail_fem_joint_specs}). The validation process, detailed in Section \ref{ssec:detail_fem_joint_validation}, begins with the integration of the joints as ``reduced order flexible bodies'' in the nano-hexapod model, followed by the development of computationally efficient lower-order models that preserve the essential dynamic characteristics. \section{Bending and Torsional Stiffness} -\label{sec:org582c93a} \label{ssec:detail_fem_joint_bending} The presence of bending stiffness in flexible joints causes the forces applied by the struts to deviate from the strut direction. @@ -679,7 +660,6 @@ A parallel analysis of torsional stiffness revealed similar dynamic effects, tho \end{figure} \section{Axial Stiffness} -\label{sec:org05206a1} \label{ssec:detail_fem_joint_axial} The limited axial stiffness (\(k_a\)) of flexible joints introduces an additional compliance between the actuation point and the measurement point. @@ -735,7 +715,6 @@ Based on this analysis, an axial stiffness specification of \(100\,N/\mu m\) was \end{figure} \section{Specifications and Design flexible joints} -\label{sec:org2310ef0} \label{ssec:detail_fem_joint_specs} The design of flexible joints for precision applications requires careful consideration of multiple mechanical characteristics. @@ -784,7 +763,6 @@ The final design, featuring a neck dimension of 0.25mm, achieves mechanical prop \end{figure} \section{Validation with the Nano-Hexapod} -\label{sec:org4e972fc} \label{ssec:detail_fem_joint_validation} The designed flexible joint was first validated through integration into the nano-hexapod model using reduced-order flexible bodies derived from finite element analysis. @@ -822,7 +800,6 @@ While additional degrees of freedom could potentially capture more dynamic featu \end{figure} \chapter*{Conclusion} -\label{sec:org5382fe7} \label{sec:detail_fem_conclusion} In this chapter, the methodology of combining finite element analysis with multi-body modeling has been demonstrated and validated, proving particularly valuable for the detailed design phase of the nano-hexapod.