#+TITLE: Optimization using Finite Element Models
:DRAWER:
#+LANGUAGE: en
#+EMAIL: dehaeze.thomas@gmail.com
#+AUTHOR: Dehaeze Thomas
#+HTML_LINK_HOME: ../index.html
#+HTML_LINK_UP: ../index.html
#+HTML_HEAD:
#+HTML_HEAD:
#+BIND: org-latex-image-default-option "scale=1"
#+BIND: org-latex-image-default-width ""
#+LaTeX_CLASS: scrreprt
#+LaTeX_CLASS_OPTIONS: [a4paper, 10pt, DIV=12, parskip=full, bibliography=totoc]
#+LATEX_HEADER: \input{preamble.tex}
#+LATEX_HEADER_EXTRA: \input{preamble_extra.tex}
#+LATEX_HEADER_EXTRA: \bibliography{nass-fem.bib}
#+BIND: org-latex-bib-compiler "biber"
#+PROPERTY: header-args:matlab :session *MATLAB*
#+PROPERTY: header-args:matlab+ :comments org
#+PROPERTY: header-args:matlab+ :exports none
#+PROPERTY: header-args:matlab+ :results none
#+PROPERTY: header-args:matlab+ :eval no-export
#+PROPERTY: header-args:matlab+ :noweb yes
#+PROPERTY: header-args:matlab+ :mkdirp yes
#+PROPERTY: header-args:matlab+ :output-dir figs
#+PROPERTY: header-args:matlab+ :tangle no
#+PROPERTY: header-args:latex :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/Cloud/tikz/org/}{config.tex}")
#+PROPERTY: header-args:latex+ :imagemagick t :fit yes
#+PROPERTY: header-args:latex+ :iminoptions -scale 100% -density 150
#+PROPERTY: header-args:latex+ :imoutoptions -quality 100
#+PROPERTY: header-args:latex+ :results file raw replace
#+PROPERTY: header-args:latex+ :buffer no
#+PROPERTY: header-args:latex+ :tangle no
#+PROPERTY: header-args:latex+ :eval no-export
#+PROPERTY: header-args:latex+ :exports results
#+PROPERTY: header-args:latex+ :mkdirp yes
#+PROPERTY: header-args:latex+ :output-dir figs
#+PROPERTY: header-args:latex+ :post pdf2svg(file=*this*, ext="png")
:END:
#+latex: \clearpage
* Build :noexport:
#+NAME: startblock
#+BEGIN_SRC emacs-lisp :results none :tangle no
(add-to-list 'org-latex-classes
'("scrreprt"
"\\documentclass{scrreprt}"
("\\chapter{%s}" . "\\chapter*{%s}")
("\\section{%s}" . "\\section*{%s}")
("\\subsection{%s}" . "\\subsection*{%s}")
("\\paragraph{%s}" . "\\paragraph*{%s}")
))
;; Remove automatic org heading labels
(defun my-latex-filter-removeOrgAutoLabels (text backend info)
"Org-mode automatically generates labels for headings despite explicit use of `#+LABEL`. This filter forcibly removes all automatically generated org-labels in headings."
(when (org-export-derived-backend-p backend 'latex)
(replace-regexp-in-string "\\\\label{sec:flexible_org[a-f0-9]+}\n" "" text)))
(add-to-list 'org-export-filter-headline-functions
'my-latex-filter-removeOrgAutoLabels)
;; Remove all org comments in the output LaTeX file
(defun delete-org-comments (backend)
(loop for comment in (reverse (org-element-map (org-element-parse-buffer)
'comment 'identity))
do
(setf (buffer-substring (org-element-property :begin comment)
(org-element-property :end comment))
"")))
(add-hook 'org-export-before-processing-hook 'delete-org-comments)
;; Use no package by default
(setq org-latex-packages-alist nil)
(setq org-latex-default-packages-alist nil)
;; Do not include the subtitle inside the title
(setq org-latex-subtitle-separate t)
(setq org-latex-subtitle-format "\\subtitle{%s}")
(setq org-export-before-parsing-hook '(org-ref-glossary-before-parsing
org-ref-acronyms-before-parsing))
#+END_SRC
* Notes :noexport:
** Notes
Prefix is =detail_fem=
*2DoF and FEM model of the APA is already described in C1 (experimental) report*.
Choice of actuator
- [X] file:~/Cloud/work-projects/ID31-NASS/matlab/nass-simscape/org/optimal_stiffness_disturbances.org
No so interesting
- [X] file:~/Cloud/work-projects/ID31-NASS/documents/state-of-thesis-2020/index.org
- [X] [[file:~/Cloud/meetings/esrf-meetings/2020-09-14-Nano-Hexapod-Design-Review/review-140920.pdf]]
Check the two reports:
- [X] This section should be included here: [[file:~/Cloud/work-projects/ID31-NASS/matlab/test-bench-apa300ml/test-bench-apa300ml.org::*Compare with the FEM/Simscape Model][Compare with the FEM/Simscape Model]]
For instance show that from FEM (/super element/, simscape), we get the characteristics found in the datasheet: stroke, stiffness, resonance frequency, etc...
- [X] Look here: https://gitlab.esrf.fr/dehaeze/nass-fem/-/tree/master?ref_type=heads
** 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>
** DONE [#A] Add APA300ML as actuator type for the NASS model
CLOSED: [2025-02-25 Tue 00:17]
- [X] Should be able to use =ga= and =gs=
- [X] Look at what is done in C6 section
** TODO [#A] Check all figures
- [ ] Legend
- [ ] Units
- [ ] Notations
- [ ] ...
** DONE [#A] Little review of flexible joints for spherical and universal joints
CLOSED: [2025-02-26 Wed 10:13]
For Stewart platform:
- ID16a [[cite:&villar18_nanop_esrf_id16a_nano_imagin_beaml]]
- DCM
** TODO [#B] Talk about limited damping added with IFF when considering flexible joint stiffness
[[file:~/Cloud/work-projects/ID31-NASS/matlab/stewart-simscape/org/control-active-damping.org::*Effect of the Flexible Joint stiffness and Actuator amplification on the Dynamics][Effect of the Flexible Joint stiffness and Actuator amplification on the Dynamics]]
[[file:~/Cloud/work-projects/ID31-NASS/matlab/stewart-simscape/org/control-active-damping.org::*Obtained Damping][Obtained Damping]]
* Introduction :ignore:
- In the detail design phase, one goal is to optimize the design of the nano-hexapod
- Parts are usually optimized using Finite Element Models that are used to estimate the static and dynamical properties of parts
- However, it is important to see how to dynamics of each part combines with the nano-hexapod and with the micro-station.
One option would be to use a FEM of the complete NASS, but that would be very complex and it would be difficult to perform simulations of experiments with real time control implemented.
- The idea is therefore to combine FEM with the multi body model of the NASS.
To do so, Reduced Order Flexible Bodies are used (Section ref:sec:detail_fem_super_element)
- The theory is described
- The method is validated using experimental measurements
- Two main elements of the nano-hexapod are then optimized:
- The actuator (Section ref:sec:detail_fem_actuator)
- The flexible joints (Section ref:sec:detail_fem_joint)
* Reduced order flexible bodies
<>
** Introduction :ignore:
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.
However, a methodological bridge between these two analytical approaches has been established, whereby components whose dynamical properties have been determined through FEA can be successfully integrated into multi-body models [[cite:&hatch00_vibrat_matlab_ansys]].
This combined multibody-FEA modeling approach presents significant advantages, as it enables the selective application of FEA modeling to specific elements while maintaining the computational efficiency of multi-body analysis for the broader system [[cite:&rankers98_machin]].
The investigation of this hybrid modeling approach is structured in three sections.
First, the fundamental principles and methodological approaches of this modeling framework are introduced (Section ref:ssec:detail_fem_super_element_theory).
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).
** 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 :noweb yes
<>
#+end_src
#+begin_src matlab :eval no :noweb yes
<>
#+end_src
#+BEGIN_SRC matlab
%% 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
#+END_SRC
#+begin_src matlab :noweb yes
<>
#+end_src
** Procedure
<>
In this modeling approach, some components within the multi-body framework are represented as /reduced-order flexible bodies/, wherein their modal behavior is characterized through reduced mass and stiffness matrices derived from finite element analysis (FEA) models.
These matrices are generated via modal reduction techniques, specifically through the application of component mode synthesis (CMS), thus establishing this design approach as a combined multibody-FEA methodology.
Standard FEA implementations typically involve thousands or even hundreds of thousands of DoF, rendering direct integration into multi-body simulations computationally prohibitive.
The objective of modal reduction is therefore to substantially decrease the number of DoF while preserving the essential dynamic characteristics of the component.
The procedure for implementing this reduction involves several distinct stages.
Initially, the component is modeled in a finite element software with appropriate material properties and boundary conditions.
Subsequently, interface frames are defined at locations where the multi-body model will establish connections with the component.
These frames serve multiple functions, including connecting to other parts, applying forces and torques, and measuring relative motion between defined frames.
Following the establishment of these interface parameters, modal reduction is performed using the Craig-Bampton method [[cite:&craig68_coupl_subst_dynam_analy]] (also known as the "fixed-interface method"), a technique that transforms the extensive FEA degrees of freedom into a significantly reduced set of retained degrees of freedom.
This transformation typically reduces the model complexity from hundreds of thousands to fewer than 100 DoF.
The number of degrees of freedom in the reduced model is determined by eqref:eq:detail_fem_model_order where $n$ represents the number of defined frames and $p$ denotes the number of additional modes to be modeled.
The outcome of this procedure is an $m \times m$ set of reduced mass and stiffness matrices, which can subsequently be incorporated into the multi-body model to represent the component's dynamic behavior.
\begin{equation}\label{eq:detail_fem_model_order}
m = 6 \times n + p
\end{equation}
** Example with an Amplified Piezoelectric Actuator
<>
**** Introduction :ignore:
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.
Additionally, an Amplified Piezoelectric Actuator (the APA95ML shown in Figure ref:fig:detail_fem_apa95ml_picture) was available in the laboratory for experimental testing.
The APA consists of multiple piezoelectric stacks arranged horizontally (depicted in blue in Figure ref:fig:detail_fem_apa95ml_picture) and of an amplifying shell structure (shown in red) that serves two purposes: the application of pre-stress to the piezoelectric elements and the amplification of their displacement into the vertical direction [[cite:&claeyssen07_amplif_piezoel_actuat]].
The selection of the APA for validation purposes was further justified by its capacity to simultaneously demonstrate multiple aspects of the modeling framework.
The specific design of the APA allows for the simultaneous modeling of a mechanical structure analogous to a flexible joint, piezoelectric actuation, and piezoelectric sensing, thereby encompassing the principal elements requiring validation.
#+attr_latex: :options [b]{0.48\linewidth}
#+begin_minipage
#+name: fig:detail_fem_apa95ml_picture
#+caption: Picture of the APA95ML
#+attr_latex: :float nil :scale 1
[[file:figs/detail_fem_apa95ml_picture.png]]
#+end_minipage
\hfill
#+attr_latex: :options [b]{0.48\linewidth}
#+begin_minipage
#+latex: \centering
#+attr_latex: :environment tabularx :width 0.7\linewidth :placement [b] :align Xc
#+attr_latex: :booktabs t :float nil :center nil
| *Parameter* | *Value* |
|----------------+---------------|
| Nominal Stroke | $100\,\mu m$ |
| Blocked force | $2100\,N$ |
| Stiffness | $21\,N/\mu m$ |
#+latex: \captionof{table}{\label{tab:detail_fem_apa95ml_specs}APA95ML specifications}
#+end_minipage
**** Finite Element Model
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.
#+name: tab:detail_fem_material_properties
#+caption: Material properties used for FEA modal reduction model. $E$ is the Young's modulus, $\nu$ the Poisson ratio and $\rho$ the material density
#+attr_latex: :environment tabularx :width 0.7\linewidth :align lXXX
#+attr_latex: :center t :booktabs t
| | $E$ | $\nu$ | $\rho$ |
|------------------------------+-------------+--------+-----------------------|
| Stainless Steel | $190\,GPa$ | $0.31$ | $7800\,\text{kg}/m^3$ |
| Piezoelectric Ceramics (PZT) | $49.5\,GPa$ | $0.31$ | $7800\,\text{kg}/m^3$ |
The definition of interface frames, or "remote points", constitute a critical aspect of the model preparation.
Seven frames were established: one frame at the two ends of each piezoelectric stack to facilitate strain measurement and force application, and additional frames at the top and bottom of the structure to enable connection with external elements in the multi-body simulation.
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}).
#+attr_latex: :options [htbp]
#+begin_figure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_apa95ml_mesh}Obtained mesh and "remote points"}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :scale 1
[[file:figs/detail_fem_apa95ml_mesh.png]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_apa_model_schematic}Inclusion in multi-body model}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :scale 1
[[file:figs/detail_fem_apa_modal_schematic.png]]
#+end_subfigure
#+end_figure
**** Super Element in the Multi-Body Model
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.
Frame $\{4\}$ was connected to the "world" frame, while frame $\{6\}$ was coupled to a vertically guided payload.
In this example, two piezoelectric stacks were used for actuation while one piezoelectric stack was used as a force sensor.
Therefore, a force source $F_a$ operating between frames $\{3\}$ and $\{2\}$ was used, while a displacement sensor $d_L$ between frames $\{1\}$ and $\{7\}$ was used for the sensor stack.
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.
**** Sensor and Actuator "constants"
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.
From [[cite:&fleming14_desig_model_contr_nanop_system p. 123]], the relation between relative displacement $d_L$ of the sensor stack and generated voltage $V_s$ is given by eqref:eq:detail_fem_dl_to_vs.
\begin{equation}\label{eq:detail_fem_dl_to_vs}
V_s = g_s \cdot d_L, \quad g_s = \frac{d_{33}}{\epsilon^T s^D n}
\end{equation}
From [[cite:&fleming10_integ_strain_force_feedb_high]] the relation between the force $F_a$ and the applied voltage $V_a$ is given by eqref:eq:detail_fem_va_to_fa.
\begin{equation}\label{eq:detail_fem_va_to_fa}
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].
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
#+caption: Stack Parameters
#+attr_latex: :environment tabularx :width 0.4\linewidth :align Xcc
#+attr_latex: :center t :booktabs t
| Parameter | Unit | Value |
|----------------+-----------+------------|
| Nominal Stroke | $\mu m$ | 20 |
| Blocked force | $N$ | 4700 |
| Stiffness | $N/\mu m$ | 235 |
| Voltage Range | $V$ | -20 to 150 |
| Capacitance | $\mu F$ | 4.4 |
| Length | $mm$ | 20 |
| Stack Area | $mm^2$ | 10x10 |
The properties of this "THP5H" material used to compute $g_a$ and $g_s$ are listed in Table ref:tab:test_apa_piezo_properties.
From these parameters, $g_s = 5.1\,V/\mu m$ and $g_a = 26\,N/V$ were obtained.
#+name: tab:test_apa_piezo_properties
#+caption: Piezoelectric properties used for the estimation of the sensor and actuators sensitivities
#+attr_latex: :environment tabularx :width 1\linewidth :align ccX
#+attr_latex: :center t :booktabs t
| *Parameter* | *Value* | *Description* |
|----------------+----------------------------+--------------------------------------------------------------|
| $d_{33}$ | $680 \cdot 10^{-12}\,m/V$ | Piezoelectric constant |
| $\epsilon^{T}$ | $4.0 \cdot 10^{-8}\,F/m$ | Permittivity under constant stress |
| $s^{D}$ | $21 \cdot 10^{-12}\,m^2/N$ | Elastic compliance understand constant electric displacement |
| $c^{E}$ | $48 \cdot 10^{9}\,N/m^2$ | Young's modulus of elasticity |
| $L$ | $20\,mm$ per stack | Length of the stack |
| $A$ | $10^{-4}\,m^2$ | Area of the piezoelectric stack |
| $n$ | $160$ per stack | Number of layers in the piezoelectric stack |
#+begin_src matlab
%% 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]
#+end_src
**** Identification of the APA Characteristics
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.
#+begin_src matlab
%% 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');
#+end_src
The stiffness of the APA95ML was estimated from the multi-body model by computing the axial compliance of the APA95ML (Figure ref:fig:detail_fem_apa95ml_compliance), which corresponds to the transfer function from a vertical force applied between the two interface frames to the relative vertical displacement between these two frames.
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$.
#+begin_src matlab
%% Stiffness estimation
m = 0.0001; % block-free condition, no payload
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]
#+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}$.
#+begin_src matlab :exports none :results none
%% 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]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/detail_fem_apa95ml_compliance.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:detail_fem_apa95ml_compliance
#+caption: Estimated compliance of the APA95ML
#+RESULTS:
[[file:figs/detail_fem_apa95ml_compliance.png]]
In order to estimate the stroke of the APA95ML, first the mechanical amplification factor, defined as the ratio between vertical displacement and horizontal stack displacement, needs to be determined.
This characteristic was quantified through analysis of the transfer function relating horizontal stack motion to vertical actuator displacement, from which an amplification factor of $1.5$ was derived.
#+begin_src matlab :exports none
%% 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]
#+end_src
The piezoelectric stacks, exhibiting a typical strain response of $0.1\,\%$ relative to their length (here equal to $20\,mm$), produce an individual nominal stroke of $20\,\mu m$ (see data-sheet of the piezoelectric stacks on Table ref:tab:detail_fem_stack_parameters, page pageref:tab:detail_fem_stack_parameters).
As three stacks are used, the horizontal displacement is $60\,\mu m$.
Through the established amplification factor of 1.5, this translates to a predicted vertical stroke of $90\,\mu m$ which falls within the manufacturer-specified range of $80\,\mu m$ and $120\,\mu m$.
The high degree of concordance observed across multiple performance metrics provides a first validation of the ability to include FEM into multi-body model.
** Experimental Validation
<>
**** Introduction :ignore:
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.
The test bench illustrated in Figure ref:fig:detail_fem_apa95ml_bench was used, which consists of a $5.7\,kg$ granite suspended on top of the APA95ML.
The granite's motion was vertically guided with an air bearing system, and a fibered interferometer was used to measured its vertical displacement $y$.
A digital-to-analog converter (DAC) was used to generate the control signal $u$, which was subsequently conditioned through a voltage amplification stage providing a gain factor of $20$, ultimately yielding the effective voltage $V_a$ across the two piezoelectric stacks.
Measurement of the sensor stack voltage $V_s$ was performed using an analog-to-digital converter (ADC).
#+name: fig:detail_fem_apa95ml_bench
#+caption: 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}).
#+attr_latex: :options [htbp]
#+begin_figure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_apa95ml_bench_picture}Picture of the test bench}
#+attr_latex: :options {0.34\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.95\linewidth
[[file:figs/detail_fem_apa95ml_bench_picture.png]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_apa95ml_bench_schematic}Schematic with signals}
#+attr_latex: :options {0.72\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.95\linewidth
[[file:figs/detail_fem_apa95ml_bench_schematic.png]]
#+end_subfigure
#+end_figure
**** Comparison of the dynamics
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:&pintelon12_system_ident, chap. 5]].
The most used ones are impulses (particularly suited to modal analysis), steps, random noise signals, and multi-sine excitations.
During all this experimental work, random noise excitation was predominantly employed.
The designed excitation signal is then generated and both input and output signals are synchronously acquired.
From the obtained input and output data, the frequency response functions were derived.
To improve the quality of the obtained frequency domain data, averaging and windowing were used [[cite:&pintelon12_system_ident, chap. 13]]..
The obtained frequency response functions from $V_a$ to $V_s$ and to $y$ are compared with the theoretical predictions derived from the multi-body model in Figure ref:fig:detail_fem_apa95ml_comp_plant.
The difference in phase between the model and the measurements can be attributed to the sampling time of $0.1\,ms$ and to additional delays induced by electronic instrumentation related to the interferometer.
The presence of a non-minimum phase zero in the measured system response (Figure ref:fig:detail_fem_apa95ml_comp_plant_sensor), shall be addressed during the experimental phase.
Regarding the amplitude characteristics, the constants $g_a$ and $g_s$ could be further refined through calibration against the experimental data.
#+begin_src matlab
%% 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]
% 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'};
#+end_src
#+begin_src matlab :exports none :results none
%% Comparison of the identified transfer function from Va to di to the multi-body model
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]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_apa95ml_comp_plant_actuator.pdf', 'width', 'half', 'height', 600);
#+end_src
#+begin_src matlab :exports none :results none
%% 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]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_apa95ml_comp_plant_sensor.pdf', 'width', 'half', 'height', 600);
#+end_src
#+name: fig:detail_fem_apa95ml_comp_plant
#+caption: Comparison of the measured frequency response functions and the identified dynamics from the finite element model of the APA95ML. Both for the dynamics from $V_a$ to $y$ (\subref{fig:detail_fem_apa95ml_comp_plant_actuator}) and from $V_a$ to $V_s$ (\subref{fig:detail_fem_apa95ml_comp_plant_sensor})
#+attr_latex: :options [htbp]
#+begin_figure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_apa95ml_comp_plant_actuator}from $V_a$ to $y$}
#+attr_latex: :options {0.49\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.95\linewidth
[[file:figs/detail_fem_apa95ml_comp_plant_actuator.png]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_apa95ml_comp_plant_sensor}from $V_a$ to $V_s$}
#+attr_latex: :options {0.49\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.95\linewidth
[[file:figs/detail_fem_apa95ml_comp_plant_sensor.png]]
#+end_subfigure
#+end_figure
**** Integral Force Feedback with APA
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.
The IFF controller implementation, defined in equation ref:eq:detail_fem_iff_controller, incorporated a tunable gain parameter $g$ and was designed to provide integral action near the system resonances and to limit the low frequency gain using an high pass filter.
\begin{equation}\label{eq:detail_fem_iff_controller}
K_{\text{IFF}}(s) = \frac{g}{s + 2\cdot 2\pi} \cdot \frac{s}{s + 0.5 \cdot 2\pi}
\end{equation}
The theoretical damped dynamics of the closed-loop system was analyzed through using the model by computed the root locus plot shown in Figure ref:fig:detail_fem_apa95ml_iff_root_locus.
For experimental validation, six gain values were tested: $g = [0,\,10,\,50,\,100,\,500,\,1000]$.
The measured frequency responses for each gain configuration were compared with model predictions, as presented in Figure ref:fig:detail_fem_apa95ml_damped_plants.
The close agreement between experimental measurements and theoretical predictions across all gain configurations demonstrates the model's capability to accurately predict both open-loop and closed-loop system dynamics, thereby validating its utility for control system design and analysis.
#+begin_src matlab :exports none
%% 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
#+end_src
#+begin_src matlab :exports none :results none
%% 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;
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_apa95ml_iff_root_locus.pdf', 'width', 'half', 'height', 600);
#+end_src
#+begin_src matlab :exports none :results none
%% 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]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_apa95ml_damped_plants.pdf', 'width', 'half', 'height', 600);
#+end_src
#+name: fig:detail_fem_apa95ml_iff_results
#+caption: Obtained results using Integral Force Feedback with the APA95ML. Obtained closed-loop poles as a function of the controller gain $g$ are prediction by root Locus plot (\subref{fig:detail_fem_apa95ml_iff_root_locus}). Circles are predictions from the model while crosses are poles estimated from the experimental data. Damped plants estimated from the model (dashed curves) and measured ones (solid curves) are compared in (\subref{fig:detail_fem_apa95ml_damped_plants}) for all tested controller gains.
#+attr_latex: :options [htbp]
#+begin_figure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_apa95ml_iff_root_locus}Root Locus plot}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.95\linewidth
[[file:figs/detail_fem_apa95ml_iff_root_locus.png]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_apa95ml_damped_plants}Damped plants}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.95\linewidth
[[file:figs/detail_fem_apa95ml_damped_plants.png]]
#+end_subfigure
#+end_figure
** Conclusion
:PROPERTIES:
:UNNUMBERED: t
:END:
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.
While this modeling approach provides accurate predictions of component behavior, the resulting model order can become prohibitively high for practical time-domain simulations.
This is exemplified by the nano-hexapod configuration, where the implementation of six Amplified Piezoelectric Actuators, each modeled with 48 degrees of freedom, yields 288 degrees of freedom only for the actuators.
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
<>
** Introduction :ignore:
Goals:
- Based on dynamical models and previous studies, extract specifications for the actuators to be included in the nano-hexapod.
Then choose the most appropriate actuator based on specifications (Section ref:ssec:detail_fem_actuator_specifications)
- Model this actuator accurately using a "reduced order flexible body" to check the dynamics and validate the choice of actuator
and validate this choice with simulations
- Development of a 2DoF model for lower order models (i.e. for simulations)
** 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 :noweb yes
<>
#+end_src
#+begin_src matlab :eval no :noweb yes
<>
#+end_src
#+BEGIN_SRC matlab
%% 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]
#+END_SRC
#+begin_src matlab :noweb yes
<>
#+end_src
** Choice of the Actuator based on Specifications
<>
From previous analysis:
- Actuator stiffness has major impact on the system dynamics and performances due to several factors:
- Spindle rotation: modification of plant dynamics and coupling increase due to Gyroscopic effects
This require to have stiffness above ~
- Limited micro-station compliance / complex dynamics:
The actuator stiffness should be small enough such that the suspension modes of the nano-hexapod are below the problematic modes of the micro-stations.
- There is therefore an intermediate stiffness that is foreseen to give the best compromise, and it is around $1\,N/\mu m$
- HAC-LAC strategy:
Actuator must include a force sensor
Because of the rotation, some stiffness should be present in parallel to the force sensor
- Limited space:
As the maximum height of the nano-hexapod is 95mm, and each strut has a flexible joint at each end, it is estimated that the maximum height of the actuator should be less than 50mm
- Stroke:
The stroke of the each actuator should be large enough such that the nano-hexapod mobility exceed the micro-station positioning errors.
Some margins should be included for mounting errors, and further flexibility of the system (for instance to perform scans with the nano-hexapod, or to align the point of interest with the rotation axis)
Actuator specifications:
- Height (<50mm)
- Stroke (~100um)
- Stiffness (0.1-1 N/um)
- Blocked force?
- Force sensor
Options:
- Two main options: piezoelectric actuators and Lorentz actuator (also known as Voice coil actuators).
Variable reluctance actuators were not considered, even though they have better efficiency than voice coil actuators, they are non linear and induce additional control complexity.
- Voice coil + relatively soft flexible guiding (1N/um):
- required force ~100N for 100um correction
This constant force/current would induce large thermal loads, that may negatively impact system's stability
Advantages of voice coil (longer strokes than piezo + allow for very low stiffness in the direction of actuation, extremely linear for high performance feedforward) are not used here.
- Piezoelectric stack actuators:
- PZT: stroke ~0.1% of its length.
- 50mm length => 50um stroke which is barely enough
- Extremely stiff, in the order of $100\,N/\mu m$, which is not wanted here.
- Amplified Piezoelectric Actuator:
- shell is used to pre-stress the piezoelectric stacks and amplify the motion (roughly by the ratio of the width over the height)
- This also reduce the stiffness in the direction of motion
- This make this design quick compact in the direction of motion (i.e. in height)
- When several stacks are used, one of them can be used as a force sensor, which is therefore very well collocated with the actuators
- Therefore, this actuator is well suited for decentralized IFF, already applied for a Stewart platform with APA [[cite:&hanieh03_activ_stewar]]
#+name: fig:detail_fem_actuator_pictures
#+caption: Example of actuators considered for the nano-hexapod. Voice coil from Sensata Technologies (\subref{fig:detail_fem_voice_coil_picture}). Piezoelectric stack actuator from Physik Instrumente (\subref{fig:detail_fem_piezo_picture}). Amplified Piezoelectric Actuator from DSM (\subref{fig:detail_fem_fpa_picture}).
#+attr_latex: :options [htbp]
#+begin_figure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_voice_coil_picture}Voice Coil}
#+attr_latex: :options {0.25\textwidth}
#+begin_subfigure
#+attr_latex: :height 4.5cm
[[file:figs/detail_fem_voice_coil_picture.jpg]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_piezo_picture}Piezoelectric stack}
#+attr_latex: :options {0.25\textwidth}
#+begin_subfigure
#+attr_latex: :height 4.5cm
[[file:figs/detail_fem_piezo_picture.jpg]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_fpa_picture}Amplified Piezoelectric Actuator}
#+attr_latex: :options {0.45\textwidth}
#+begin_subfigure
#+attr_latex: :height 3.5cm
[[file:figs/detail_fem_fpa_picture.jpg]]
#+end_subfigure
#+end_figure
Based on previous analysis, it was decided to use amplified piezoelectric actuators for the nano-hexapod.
Table ref:tab:detail_fem_piezo_act_models: compares few models that fulfill specifications.
It was decided to go for the APA300ML (shown in Figure ref:fig:detail_fem_apa300ml_picture).
One reason is that we already had experience with APA from Cedrat technologies, and the Finite Element Model was validated experimentally, so we are confident to model the APA300ML with FEA and include it in the NASS model for validation.
- Talk about piezoelectric actuator? bandwidth? noise?
- Resolution: really depends on the electrical noise (induced by DAC and voltage amplifier).
They will be chosen appropriately
#+name: tab:detail_fem_piezo_act_models
#+caption: List of some amplified piezoelectric actuators that could be used for the nano-hexapod
#+attr_latex: :environment tabularx :width 0.9\linewidth :align Xccccc
#+attr_latex: :center t :booktabs t :float t :font \scriptsize
| *Specification* | APA150M | *APA300ML* | APA400MML | FPA-0500E-P | FPA-0300E-S |
|------------------------------------+---------+------------+-----------+-------------+-------------|
| Stroke $> 100\, [\mu m]$ | 187 | 304 | 368 | 432 | 240 |
| Stiffness $\approx 1\, [N/\mu m]$ | 0.7 | 1.8 | 0.55 | 0.87 | 0.58 |
| Resolution $< 2\, [nm]$ | 2 | 3 | 4 | | |
| Blocked Force $> 100\, [N]$ | 127 | 546 | 201 | 376 | 139 |
| Height $< 50\, [mm]$ | 22 | 30 | 24 | 27 | 16 |
** APA300ML - Reduced Order Flexible Body
<>
To validate the choice of the APA300ML (Shown in Figure ref:fig:detail_fem_apa300ml_picture):
- the APA300ML is modeled using a Finite Element Software
- a /super element/ is exported and imported in Simscape where its dynamic is studied
- similarly to what was done with the APA95ML, frames defined for the /super element/ are shown in figure ref:fig:detail_fem_apa300ml_frames
#+name: fig:detail_fem_apa300ml
#+caption: Amplified Piezoelectric Actuator APA300ML. Picture shown in (\subref{fig:detail_fem_apa300ml_picture}). Frames (or "remote points") used for the modal reduction are shown in (\subref{fig:detail_fem_apa300ml_frames}).
#+attr_latex: :options [htbp]
#+begin_figure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_apa300ml_picture}Picture of the APA300ML}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.95\linewidth
[[file:figs/detail_fem_apa300ml_picture.jpg]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_apa300ml_frames}FEM of the APA300ML}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.95\linewidth
[[file:figs/detail_fem_apa300ml_frames.png]]
#+end_subfigure
#+end_figure
- For this reduced order model, 7 frames are defined and 120 additional modes are modelled for a total matrix size of 162.
- This is very large and will not be practical for simulations, but the best model accuracy was wanted for validation
- The blue frames are used to model the force sensor stack: the relative motion between the two frame is measured
- The red frames are used to model the two actuator stacks: /internal force/ are added
- One mass is fixed at one end of the piezo-electric stack actuator (remove point F), the other end is fixed to the world frame (remote point G).
- The link between mechanical properties and electrical properties was discussed in Section ref:ssec:detail_fem_super_element_validation.
As the stacks are the same between the APA300ML and the APA95ML, the values estimated for $g_a$ and $g_s$ are used for the APA300ML.
** Simpler 2DoF Model of the APA300ML
<>
**** Introduction :ignore:
- /super-element/ order is quite large, and therefore not practical for simulations
- the goal here is to develop a low order model, that still represents wanted characteristics of the APA300ML:
- axial stiffness
- actuator and force sensor characteristics
- what is not modelled:
- higher order modes
- the flexibility of the APA in the other directions
- Therefore this model can be useful for simulations as it contains a very limited number of states, but when more complex dynamics of the APA is to be modelled, a flexible model will be used.
**** 2DoF Model
The model is adapted from cite:souleille18_concep_activ_mount_space_applic.
It can be decomposed into three components:
- the shell whose axial properties are represented by $k_1$ and $c_1$
- the actuator stacks whose contribution to the axial stiffness is represented by $k_a$ and $c_a$.
The force source $f$ represents the axial force induced by the force sensor stacks.
The sensitivity $g_a$ (in $N/m$) is used to convert the applied voltage $V_a$ to the axial force $f$
- the sensor stack whose contribution to the axial stiffness is represented by $k_e$ and $c_e$.
A sensor measures the stack strain $d_e$ which is then converted to a voltage $V_s$ using a sensitivity $g_s$ (in $V/m$)
Such a simple model has some limitations:
- it only represents the axial characteristics of the APA as it is modeled as infinitely rigid in the other directions
- some physical insights are lost, such as the amplification factor and the real stress and strain in the piezoelectric stacks
- the creep and hysteresis of the piezoelectric stacks are not modeled as the model is linear
The main advantage is that this model is very simple, only adds 4 states
#+name: fig:detail_fem_apa_2dof_model
#+caption: Schematic of the 2DoF model of the Amplified Piezoelectric Actuator
[[file:figs/detail_fem_apa_2dof_model.png]]
**** Parameter Tuning
9 parameters ($m$, $k_1$, $c_1$, $k_e$, $c_e$, $k_a$, $c_a$, $g_s$ and $g_a$) have to be tuned such that the dynamics of the model (Figure ref:fig:detail_fem_apa_2dof_model) well represents the identified dynamics using the FEM.
- Mass is 5kg (similar to the test bench)
- Tune the parameters:
- From the first zero of the transfer function from Va to Vs, k1 and c1 are tuned
- From the first pole of the transfer function from Va to y, ka, ca, ke, ce are tuned
- because the actuator and sensor stacks are physically the same, we suppose
Then, it is reasonable to assume that the sensor stacks and the two actuator stacks have identical mechanical characteristics.
Therefore, we have $k_e = 2 k_a$ and $c_e = 2 c_a$ as the actuator stack is composed of two stacks in series.
- In the last step, $g_s$ and $g_a$ for the 2DoF motion can be tuned to match the gain of the transfer functions extracted from the FEM
- Found parameters are summarized in Table ref:tab:detail_fem_apa300ml_2dof_parameters
- Comparison of the transfer functions extracted from the high order flexible model with the 4th order (2DoF) model is done in Figure ref:fig:detail_fem_apa300ml_comp_fem_2dof_fem_2dof.
Good match is obtained.
Of course, higher order modes are not represented by the 2DoF model, nor the limited stiffness in the other directions.
#+name: tab:detail_fem_apa300ml_2dof_parameters
#+caption: Summary of the obtained parameters for the 2 DoF APA300ML model
#+attr_latex: :environment tabularx :width 0.3\linewidth :align cc
#+attr_latex: :center t :booktabs t
| *Parameter* | *Value* |
|-------------+-----------------|
| $k_1$ | $0.30\,N/\mu m$ |
| $k_e$ | $4.3\, N/\mu m$ |
| $k_a$ | $2.15\,N/\mu m$ |
| $c_1$ | $18\,Ns/m$ |
| $c_e$ | $0.7\,Ns/m$ |
| $c_a$ | $0.35\,Ns/m$ |
| $g_a$ | $2.7\,N/V$ |
| $g_s$ | $0.53\,V/\mu m$ |
#+begin_src matlab
%% Identify dynamics with "Reduced Order Flexible Body"
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');
#+end_src
#+begin_src matlab :exports none :results none
%% 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]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_apa300ml_comp_fem_2dof_actuator.pdf', 'width', 'half', 'height', 600);
#+end_src
#+begin_src matlab :exports none :results none
%% 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]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_apa300ml_comp_fem_2dof_force_sensor.pdf', 'width', 'half', 'height', 600);
#+end_src
#+name: fig:detail_fem_apa300ml_comp_fem_2dof_fem_2dof
#+caption: Comparison of the transfer functions extracted from the finite element model of the APA300ML and of the 2DoF model. Both for the dynamics from $V_a$ to $d_i$ (\subref{fig:detail_fem_apa300ml_comp_fem_2dof_actuator}) and from $V_a$ to $V_s$ (\subref{fig:detail_fem_apa300ml_comp_fem_2dof_force_sensor})
#+attr_latex: :options [htbp]
#+begin_figure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_apa300ml_comp_fem_2dof_actuator}from $V_a$ to $d_i$}
#+attr_latex: :options {0.49\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.95\linewidth
[[file:figs/detail_fem_apa300ml_comp_fem_2dof_actuator.png]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_apa300ml_comp_fem_2dof_force_sensor}from $V_a$ to $V_s$}
#+attr_latex: :options {0.49\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.95\linewidth
[[file:figs/detail_fem_apa300ml_comp_fem_2dof_force_sensor.png]]
#+end_subfigure
#+end_figure
** Electrical characteristics of the APA
- Mechanical equations and electrical equations are coupled
- This means for instance, that the stiffness of the piezoelectric stack (i.e. the APA) depends on the electrical boundaries of the stacks:
- Short circuited stacks are less stiff than open-circuited ones
- This effect is quite small: example with the APA95ML (Figure ref:fig:detail_fem_apa95ml_effect_electrical_boundaries)
transfer function from Va to di are estimated with the force sensor stack being short circuited or open-circuited.
- In the model used, the electrical phenomena are not modelled.
But as this effect is small, it should be fine
- The electrical characteristics of the APA are very important both from the voltage amplifier side and the ADC measuring the force sensor voltage.
This will be discussed in chapter "instrumentation"
#+begin_src matlab
%% 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)));
#+end_src
#+begin_src matlab :exports none :results none
%% 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'); 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]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/detail_fem_apa95ml_effect_electrical_boundaries.pdf', 'width', 'wide', 'height', 600);
#+end_src
#+name: fig:detail_fem_apa95ml_effect_electrical_boundaries
#+caption: Effect of the electrical bondaries of the force sensor stack on the APA95ML resonance frequency
#+RESULTS:
[[file:figs/detail_fem_apa95ml_effect_electrical_boundaries.png]]
** Validation with the Nano-Hexapod
NASS model + FEM model (or just 2DoF) of APA300ML => validation (based on what?)
- Compare 2DoF model and FEM (Figure ref:fig:detail_fem_actuator_fem_vs_perfect_plants)
- HAC plant
- IFF Plant
- Very similar => can use 2nd order actuator models
- Talk about model order
- 2DoF actuators: 24 states
- FEM actuators:
here matrices have a size of 36
36*6+12 => ~300
#+begin_src matlab
%% 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'};
#+end_src
#+begin_src matlab :exports none :results none
%% 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');
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_actuator_fem_vs_perfect_hac_plant.pdf', 'width', 'half', 'height', 600);
#+end_src
#+begin_src matlab :exports none :results none
%% 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');
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_actuator_fem_vs_perfect_iff_plant.pdf', 'width', 'half', 'height', 600);
#+end_src
#+name: fig:detail_fem_actuator_fem_vs_perfect_plants
#+caption: Comparison of the dynamics obtained between a nano-hexpod having the actuators modeled with FEM and a nano-hexapod having actuators modelled a 2DoF system. Both from actuator force $\bm{f}$ to strut motion measured by external metrology $\bm{\epsilon}_{\mathcal{L}}$ (\subref{fig:detail_fem_actuator_fem_vs_perfect_iff_plant}) and to the force sensors $\bm{f}_m$ (\subref{fig:detail_fem_actuator_fem_vs_perfect_hac_plant}).
#+attr_latex: :options [htbp]
#+begin_figure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_actuator_fem_vs_perfect_hac_plant}$\bm{f}$ to $\bm{\epsilon}_{\mathcal{L}}$}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.9\linewidth
[[file:figs/detail_fem_actuator_fem_vs_perfect_hac_plant.png]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_actuator_fem_vs_perfect_iff_plant}$\bm{f}$ to $\bm{f}_m$}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.9\linewidth
[[file:figs/detail_fem_actuator_fem_vs_perfect_iff_plant.png]]
#+end_subfigure
#+end_figure
* Flexible Joint
<>
** Notes :noexport:
*Here we just look at the wanted stiffness in different directions?*
*Or also design of the flexible joint*?
Because we talk about detail design in last section
*In this section, we talk only about important part of the flexible joint, and not connection to the struts of plates*
Prefix is =fem_joints=
- [ ] *Should I consider rigid micro-station for simplified analysis?*
- [ ] https://research.tdehaeze.xyz/nass-simscape/docs/flexible_joints_study.html
- [X] look at section 5.4 file:/home/thomas/Cloud/work-projects/ID31-NASS/documents/state-of-thesis-2020/index.org
- [X] Show that flexible joints' stiffnesses has little impact on the Stewart platform behavior
Find paper of MrInroy that talks about that with associated conditions: [[cite:&mcinroy02_model_desig_flexur_joint_stewar]]
- [X] Look here: https://gitlab.esrf.fr/dehaeze/nass-fem/-/tree/master?ref_type=heads
- [X] Also check start of this report: file:~/Cloud/work-projects/ID31-NASS/matlab/nass-simscape/org/nano_hexapod.org to show the Simscape model of the Flexible joint
*No this will be for the detail design phase*
Outline:
- Review of flexible joints ?
- Imperfection of the flexible joint: Model
- Study of the effect of limited stiffness in constrain directions and non-null stiffness in other directions
Compare with perfect flexible joint case
- Obtained Specification
- Design optimisation (FEM)
- Implementation of flexible elements in the Simscape model: close to simplified model
** Introduction :ignore:
The flexible joints have few advantages compared to conventional joints such as the *absence of wear, friction and backlash* which allows extremely high-precision (predictable) motion.
The parasitic bending and torsional stiffness of these joints usually induce some *limitation on the control performance*. [[cite:&mcinroy02_model_desig_flexur_joint_stewar]]
In this document is studied the effect of the mechanical behavior of the flexible joints that are located the extremities of each nano-hexapod's legs.
Ideally, we want the x and y rotations to be free and all the translations to be blocked.
However, this is never the case and be have to consider:
- Non-null bending stiffnesses
- Non-null radial compliance
- Axial stiffness in the direction of the legs
This may impose some limitations, also, the goal is to specify the required joints stiffnesses.
Say that for simplicity (reduced number of parts, etc.), we consider the same joints for the fixed based and the top platform.
*Outline*:
- Perfect flexible joint
- Imperfection of the flexible joint: Model
- Study of the effect of limited stiffness in constrain directions and non-null stiffness in other directions
- Obtained Specification
- Design optimisation (FEM)
- Implementation of flexible elements in the Simscape model: close to simplified model
** 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 :noweb yes
<>
#+end_src
#+begin_src matlab :eval no :noweb yes
<>
#+end_src
#+BEGIN_SRC matlab
%% Linearization options
opts = linearizeOptions;
opts.SampleTime = 0;
%% Open Simscape Model
mdl = 'detail_fem_nass'; % Name of the Simulink File
open(mdl); % Open Simscape Model
#+END_SRC
#+begin_src matlab :noweb yes
<>
#+end_src
** Flexible joints for Stewart platforms
Review of different types of flexible joints for Stewart plaftorms (see Figure ref:fig:detail_fem_joints_examples).
Typical specifications:
- Bending stroke (i.e. long life time by staying away from yield stress, even at maximum deflection/load)
- Axial stiffness
- Bending stiffness
- Maximum axial load
- Well defined rotational axes
Typical values?
- $K_{\theta, \phi} = 15\,[Nm/rad]$ stiffness in flexion
- $K_{\psi} = 20\,[Nm/rad]$ stiffness in torsion
- \[ K_a = 60\,[N/\mu m] \] axial stiffness
#+name: fig:detail_fem_joints_examples
#+caption: Example of different flexible joints geometry used for Stewart platforms. (\subref{fig:detail_fem_joints_yang}) [[cite:&yang19_dynam_model_decoup_contr_flexib]]. (\subref{fig:detail_fem_joints_preumont}) [[cite:&preumont07_six_axis_singl_stage_activ]]. (\subref{fig:detail_fem_joints_wire}) [[cite:&du14_piezo_actuat_high_precis_flexib]].
#+attr_latex: :options [htbp]
#+begin_figure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_joints_yang}}
#+attr_latex: :options {0.35\textwidth}
#+begin_subfigure
#+attr_latex: :height 5cm
[[file:figs/detail_fem_joints_yang.png]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_joints_preumont}}
#+attr_latex: :options {0.3\textwidth}
#+begin_subfigure
#+attr_latex: :height 5cm
[[file:figs/detail_fem_joints_preumont.png]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_joints_wire}}
#+attr_latex: :options {0.3\textwidth}
#+begin_subfigure
#+attr_latex: :height 5cm
[[file:figs/detail_fem_joints_wire.png]]
#+end_subfigure
#+end_figure
** Bending and Torsional Stiffness
<>
Because of bending stiffness of the flexible joints, the forces applied by the struts are no longer aligned with the struts (additional forces applied by the "spring force" of the flexible joints).
In this section, we wish to study the effect of the rotation flexibility of the nano-hexapod joints.
- To simplify the analysis, the micro-station is considered rigid, and only the nano-hexapod is considered with:
- 1dof actuators, k=1N/um, without parallel stiffness to the force sensors
- The bending stiffness of all joints are varied and the dynamics is identified
#+begin_src matlab
%% 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
#+end_src
HAC plant (transfer function from f to dL, as measured by the external metrology):
- It increase the coupling at low frequency, but is kept to small values for realistic values of the bending stiffness (Figure ref:fig:detail_fem_joints_bending_stiffness_hac_plant)
- Bending stiffness does not impact significantly the HAC plant.
The added stiffness increases the frequency of the suspension modes
Condition in [[cite:&mcinroy02_model_desig_flexur_joint_stewar]] to have forces aligned with the struts when considering rotational stiffness: kr << k*l^2
For the current nano hexapod configuration, it correspond to << 9000 Nm/rad.
This may be an issue for soft nano-hexapod (for instance k = 1e4 => << 90) => have to design very soft flexible joints.
Here, having relatively stiff actuators render this condition easier to achieve.
IFF Plant:
- Having bending stiffness adds complex conjugate zero at low frequency (Figure ref:fig:detail_fem_joints_bending_stiffness_iff_plant)
- Similar to having a stiffness in parallel to the struts (i.e., to the force sensor).
This can be explained since even if the force sensor is removed (i.e. zero axial stiffness of the strut), the strut will still act as a spring between the mobile and fixed plates because of the bending stiffness of the flexible joints.
The frequency of the zero gives an idea of the stiffness contribution of the flexible joint bending stiffness
- They therefore impose limitation for decentralized IFF, as discussed in [[cite:&preumont07_six_axis_singl_stage_activ]]
- This can be seen in the root locus plot of Figure ref:fig:detail_fem_joints_bending_stiffness_iff_locus_1dof
#+begin_src matlab :exports none :results none
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');
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_joints_bending_stiffness_hac_plant.pdf', 'width', 'half', 'height', 600);
#+end_src
#+begin_src matlab :exports none :results none
%% 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}("fn"+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');
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');
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_joints_bending_stiffness_iff_plant.pdf', 'width', 'half', 'height', 600);
#+end_src
#+name: fig:detail_fem_joints_bending_stiffness_plants
#+caption: Effect of bending stiffness of the flexible joints on the plant dynamics. Both from actuator force $\bm{f}$ to strut motion measured by external metrology $\bm{\epsilon}_{\mathcal{L}}$ (\subref{fig:detail_fem_joints_bending_stiffness_hac_plant}) and to the force sensors $\bm{f}_m$ (\subref{fig:detail_fem_joints_bending_stiffness_iff_plant})
#+attr_latex: :options [h!tbp]
#+begin_figure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_joints_bending_stiffness_hac_plant}$\bm{f}$ to $\bm{\epsilon}_{\mathcal{L}}$}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.9\linewidth
[[file:figs/detail_fem_joints_bending_stiffness_hac_plant.png]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_joints_bending_stiffness_iff_plant}$\bm{f}$ to $\bm{f}_m$}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.9\linewidth
[[file:figs/detail_fem_joints_bending_stiffness_iff_plant.png]]
#+end_subfigure
#+end_figure
#+begin_src matlab :exports none :results none
%% 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'};
#+end_src
#+begin_src matlab :exports none :results none
%% 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;
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_joints_bending_stiffness_iff_locus_1dof.pdf', 'width', 'half', 'height', 500);
#+end_src
However, as the APA300ML was chosen for the actuator, stiffness are already present in parallel to the force sensors:
- The dynamics is computed again for all considered values of the bending stiffnesses with the 2DoF model of the APA300ML
- Root locus for decentralized IFF are shown in Figure ref:fig:detail_fem_joints_bending_stiffness_iff_locus_apa300ml.
Now the effect of bending stiffness has little effect on the attainable damping, as its contribution as "parallel stiffness" is small compared to the parallel stiffness already present in the APA300ML.
#+begin_src matlab
%% 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
#+end_src
#+begin_src matlab :exports none :results none
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;
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_joints_bending_stiffness_iff_locus_apa300ml.pdf', 'width', 'half', 'height', 500);
#+end_src
#+name: fig:detail_fem_joints_bending_stiffness_iff_locus
#+caption: Effect of bending stiffness of the flexible joints on the attainable damping with decentralized IFF. When having an actuator modelled as 1DoF without parallel stiffness to the force sensor (\subref{fig:detail_fem_joints_bending_stiffness_iff_locus_1dof}), and with the 2DoF model of the APA300ML (\subref{fig:detail_fem_joints_bending_stiffness_iff_locus_apa300ml})
#+attr_latex: :options [h!tbp]
#+begin_figure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_joints_bending_stiffness_iff_locus_1dof}1DoF actuators}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.9\linewidth
[[file:figs/detail_fem_joints_bending_stiffness_iff_locus_1dof.png]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_joints_bending_stiffness_iff_locus_apa300ml}APA300ML actuators}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.9\linewidth
[[file:figs/detail_fem_joints_bending_stiffness_iff_locus_apa300ml.png]]
#+end_subfigure
#+end_figure
Conclusion:
- Similar results for torsional stiffness, but less important
- thanks to the use of the APA, the requirements in terms of bending stiffness are less stringent
** Axial Stiffness
<>
- Adding flexibility between the actuation point and the measurement point / point of interest is always detrimental for the control performances.
This is verified, and the goal is to estimate the minimum axial stiffness that the flexible joints should have
- Here, the mass of the strut should be considered.
It is set to 112g as specified in the APA300ML specification sheet.
- Transfer functions are estimated for several axial stiffnesses (Figure ref:fig:detail_fem_joints_axial_stiffness_plants)
- IFF plant is not much affected (Figure ref:fig:detail_fem_joints_axial_stiffness_iff_plant).
Confirmed by the root locus plot of Figure ref:fig:detail_fem_joints_axial_stiffness_iff_locus
- "HAC" plant:
- Additional modes at high frequency corresponding to internal modes of the struts.
It adds coupling to the plant.
This is confirmed by computed the RGA-number for the damped plant (i.e. after applying decentralized IFF) in Figure ref:fig:detail_fem_joints_axial_stiffness_rga_hac_plant
#+begin_src matlab
%% 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
#+end_src
#+begin_src matlab :exports none :results none
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');
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_joints_axial_stiffness_hac_plant.pdf', 'width', 'half', 'height', 600);
#+end_src
#+begin_src matlab :exports none :results none
%% 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');
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_joints_axial_stiffness_iff_plant.pdf', 'width', 'half', 'height', 600);
#+end_src
#+name: fig:detail_fem_joints_axial_stiffness_plants
#+caption: Effect of axial stiffness of the flexible joints on the plant dynamics. Both from actuator force $\bm{f}$ to strut motion measured by external metrology $\bm{\epsilon}_{\mathcal{L}}$ (\subref{fig:detail_fem_joints_axial_stiffness_hac_plant}) and to the force sensors $\bm{f}_m$ (\subref{fig:detail_fem_joints_axial_stiffness_iff_plant})
#+attr_latex: :options [h!tbp]
#+begin_figure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_joints_axial_stiffness_hac_plant}$\bm{f}$ to $\bm{\epsilon}_{\mathcal{L}}$}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.9\linewidth
[[file:figs/detail_fem_joints_axial_stiffness_hac_plant.png]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_joints_axial_stiffness_iff_plant}$\bm{f}$ to $\bm{f}_m$}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.9\linewidth
[[file:figs/detail_fem_joints_axial_stiffness_iff_plant.png]]
#+end_subfigure
#+end_figure
Integral force feedback
- [ ] Maybe show the damped plants instead?
- [ ] Root Locus: not a lot of effect
#+begin_src matlab :exports none :results none
%% 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'};
#+end_src
#+begin_src matlab :exports none :results none
%% 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;
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_joints_axial_stiffness_iff_locus.pdf', 'width', 'half', 'height', 500);
#+end_src
#+begin_src matlab
%% 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
#+end_src
#+begin_src matlab :exports none :results none
%% 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;
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_joints_axial_stiffness_rga_hac_plant.pdf', 'width', 'half', 'height', 500);
#+end_src
#+name: fig:detail_fem_joints_axial_stiffness_iff_results
#+caption: Effect of axial stiffness of the flexible joints on the attainable damping with decentralized IFF (\subref{fig:detail_fem_joints_axial_stiffness_iff_locus}). Estimation of the coupling of the damped plants using the RGA-number (\subref{fig:detail_fem_joints_axial_stiffness_rga_hac_plant})
#+attr_latex: :options [h!tbp]
#+begin_figure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_joints_axial_stiffness_iff_locus}Root Locus}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :scale 0.9
[[file:figs/detail_fem_joints_axial_stiffness_iff_locus.png]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_joints_axial_stiffness_rga_hac_plant}RGA number}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :scale 0.9
[[file:figs/detail_fem_joints_axial_stiffness_rga_hac_plant.png]]
#+end_subfigure
#+end_figure
Conclusion:
- The axial stiffness of the flexible joints should be maximized to limit additional coupling at high frequency that may negatively impact the achievable bandwidth
- It should be much higher than the stiffness of the actuator
- For the nano-hexapod 100N/um is a reasonable axial stiffness specification
- Above the resonance frequency linked to the limited axial stiffness of the flexible joint, the system becomes coupled and impossible to control
- Also, loose control authority at the frequency of the zero
** Obtained design / Specifications
- Summary of specifications (Table ref:tab:detail_fem_joints_specs)
- Explain choice of geometry:
- x and y rotations are coincident
- stiffness can be easily tuned
- high axial stiffness
- Explain how it is optimized:
- Extract stiffnesses from FEM
- Parameterized model in the FE software
- Quick optimization: (few iterations, could probably increase more the axial stiffness)
- There is a trade off between high axial stiffness and low bending/torsion stiffness
- Also check the yield strength
- Show obtained geometry Figure ref:fig:detail_fem_joints_design:
- "neck" size: 0.25mm
- Characteristics of the flexible joints obtained from FEA are summarized in Table ref:tab:detail_fem_joints_specs
#+name: tab:detail_fem_joints_specs
#+caption: Specifications for the flexible joints and estimated characteristics from the Finite Element Model
#+attr_latex: :environment tabularx :width 0.5\linewidth :align Xcc
#+attr_latex: :center t :booktabs t :float t
| | *Specification* | *FEM* |
|-------------------------+------------------------+-------|
| Axial Stiffness $k_a$ | $> 100\,N/\mu m$ | 94 |
| Shear Stiffness $k_s$ | $> 1\,N/\mu m$ | 13 |
| Bending Stiffness $k_f$ | $< 100\,Nm/\text{rad}$ | 5 |
| Torsion Stiffness $k_t$ | $< 500\,Nm/\text{rad}$ | 260 |
| Bending Stroke | $> 1\,\text{mrad}$ | 24.5 |
#+name: fig:detail_fem_joints_design
#+caption: Designed flexible joints.
#+attr_latex: :options [htbp]
#+begin_figure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_joints_3d_view}3D view}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :scale 1
[[file:figs/detail_fem_joints_3d_view.png]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_joint_dimensions}Key dimensions}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :scale 1
[[file:figs/detail_fem_joint_dimensions.png]]
#+end_subfigure
#+end_figure
** Validation with the Nano-Hexapod
To validate the designed flexible joint:
- FEM: modal reduction
two interface frames are defined (Figure ref:fig:detail_fem_joints_frames)
- additional 6 modes are extracted: size of reduced order mass and stiffness matrices: $18 \times 18$
- Imported in the multi-body model
- The transfer functions from forces and torques applied between frames $\{F\}$ and $\{M\}$ to the relative displacement/rotations of the two frames is extracted.
- The stiffness characteristics of the flexible joint is estimated from the low frequency gain of the obtained transfer functions. Same values are obtained with the reduced order model and the FEM.
#+name: fig:detail_fem_joints_frames
#+caption: Defined frames for the reduced order flexible body. The two flat interfaces are considered rigid, and are linked to the two frames $\{F\}$ and $\{M\}$ both located at the center of the rotation.
[[file:figs/detail_fem_joints_frames.png]]
#+begin_src matlab
%% 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]
#+end_src
Depending on which characteristic of the flexible joint is to be modelled, several DoFs can be taken into account:
- 2DoF (universal joint) $k_f$
- 3DoF (spherical joint) taking into account torsion $k_f$, $k_t$
- 2DoF + axial stiffness $k_f$, $k_a$
- 3DoF + axial stiffness $k_f$, $k_t$, $k_a$
- 6DoF ("bushing joint") $k_f$, $k_t$, $k_a$, $k_s$
Adding more degrees of freedom:
- can represent important features
- adds model states that may not be relevant for the dynamics, and may complexity the simulations without adding much information
After testing different configurations, a good compromise was found for the modelling of the nano-hexapod flexible joints:
- bottom joints: $k_f$ and $k_a$
- top joints: $k_f$, $k_t$ and $k_a$
Talk about model order:
- with flexible joints: 252 states:
- 12 for the payload (6 dof)
- 12 for the 2DoF struts
- 216 DoF for the flexible joints (18*6*2)
- 12 states for?
- with 3dof and 4dof: 48 states
- 12 for the payload (6 dof)
- 12 for the 2DoF struts
- 12 states for the bottom joints
- 12 states for the top joints
#+begin_src matlab
%% 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'};
#+end_src
#+begin_src matlab :exports none :results none
%% 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');
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_joints_fem_vs_perfect_hac_plant.pdf', 'width', 'half', 'height', 600);
#+end_src
#+begin_src matlab :exports none :results none
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');
#+end_src
#+begin_src matlab :tangle no :exports results :results file none
exportFig('figs/detail_fem_joints_fem_vs_perfect_iff_plant.pdf', 'width', 'half', 'height', 600);
#+end_src
#+name: fig:detail_fem_joints_fem_vs_perfect_plants
#+caption: Comparison of the dynamics obtained between a nano-hexpod including joints modelled with FEM and a nano-hexapod having bottom joint modelled by bending stiffness $k_f$ and axial stiffness $k_a$ and top joints modelled by bending stiffness $k_f$, torsion stiffness $k_t$ and axial stiffness $k_a$. Both from actuator force $\bm{f}$ to strut motion measured by external metrology $\bm{\epsilon}_{\mathcal{L}}$ (\subref{fig:detail_fem_joints_fem_vs_perfect_iff_plant}) and to the force sensors $\bm{f}_m$ (\subref{fig:detail_fem_joints_fem_vs_perfect_hac_plant}).
#+attr_latex: :options [htbp]
#+begin_figure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_joints_fem_vs_perfect_hac_plant}$\bm{f}$ to $\bm{\epsilon}_{\mathcal{L}}$}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.9\linewidth
[[file:figs/detail_fem_joints_fem_vs_perfect_hac_plant.png]]
#+end_subfigure
#+attr_latex: :caption \subcaption{\label{fig:detail_fem_joints_fem_vs_perfect_iff_plant}$\bm{f}$ to $\bm{f}_m$}
#+attr_latex: :options {0.48\textwidth}
#+begin_subfigure
#+attr_latex: :width 0.9\linewidth
[[file:figs/detail_fem_joints_fem_vs_perfect_iff_plant.png]]
#+end_subfigure
#+end_figure
* Conclusion
:PROPERTIES:
:UNNUMBERED: t
:END:
<>
* Bibliography :ignore:
#+latex: \printbibliography[heading=bibintoc,title={Bibliography}]
* Helping Functions :noexport:
** Initialize Path
#+NAME: m-init-path
#+BEGIN_SRC matlab
%% Path for functions, data and scripts
addpath('./matlab/');
addpath('./matlab/src/'); % Path for scripts
addpath('./matlab/mat/'); % Path for data
addpath('./matlab/STEPS/'); % Path for Simscape Model
addpath('./matlab/subsystems/'); % Path for Subsystems Simulink files
#+END_SRC
#+NAME: m-init-path-tangle
#+BEGIN_SRC matlab
%% 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
#+END_SRC
** Initialize Simscape
** Initialize other elements
#+NAME: m-init-other
#+BEGIN_SRC matlab
%% Colors for the figures
colors = colororder;
freqs = logspace(1,3,500); % Frequency vector [Hz]
#+END_SRC
* Matlab Functions :noexport:
*** =initializeSimplifiedNanoHexapod=: Nano Hexapod
#+begin_src matlab :tangle matlab/src/initializeSimplifiedNanoHexapod.m :comments none :mkdirp yes :eval no
function [nano_hexapod] = initializeSimplifiedNanoHexapod(args)
arguments
args.type char {mustBeMember(args.type,{'none', 'stewart'})} = 'stewart'
%% initializeFramesPositions
args.H (1,1) double {mustBeNumeric, mustBePositive} = 95e-3 % Height of the nano-hexapod [m]
args.MO_B (1,1) double {mustBeNumeric} = 150e-3 % Height of {B} w.r.t. {M} [m]
%% generateGeneralConfiguration
args.FH (1,1) double {mustBeNumeric, mustBePositive} = 15e-3 % Height of fixed joints [m]
args.FR (1,1) double {mustBeNumeric, mustBePositive} = 120e-3 % Radius of fixed joints [m]
args.FTh (6,1) double {mustBeNumeric} = [220, 320, 340, 80, 100, 200]*(pi/180) % Angles of fixed joints [rad]
args.MH (1,1) double {mustBeNumeric, mustBePositive} = 15e-3 % Height of mobile joints [m]
args.MR (1,1) double {mustBeNumeric, mustBePositive} = 110e-3 % Radius of mobile joints [m]
args.MTh (6,1) double {mustBeNumeric} = [255, 285, 15, 45, 135, 165]*(pi/180) % Angles of fixed joints [rad]
%% Actuators
args.actuator_type char {mustBeMember(args.actuator_type,{'1dof', '2dof', 'flexible', 'apa300ml'})} = '1dof'
args.actuator_k (1,1) double {mustBeNumeric, mustBePositive} = 1e6
args.actuator_kp (1,1) double {mustBeNumeric, mustBeNonnegative} = 5e4
args.actuator_ke (1,1) double {mustBeNumeric, mustBePositive} = 4952605
args.actuator_ka (1,1) double {mustBeNumeric, mustBePositive} = 2476302
args.actuator_c (1,1) double {mustBeNumeric, mustBePositive} = 50
args.actuator_cp (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.actuator_ce (1,1) double {mustBeNumeric, mustBePositive} = 100
args.actuator_ca (1,1) double {mustBeNumeric, mustBePositive} = 50
args.actuator_ga (1,1) double {mustBeNumeric} = 1
args.actuator_gs (1,1) double {mustBeNumeric} = 5
%% initializeCylindricalPlatforms
args.Fpm (1,1) double {mustBeNumeric, mustBePositive} = 5 % Mass of the fixed plate [kg]
args.Fph (1,1) double {mustBeNumeric, mustBePositive} = 10e-3 % Thickness of the fixed plate [m]
args.Fpr (1,1) double {mustBeNumeric, mustBePositive} = 150e-3 % Radius of the fixed plate [m]
args.Mpm (1,1) double {mustBeNumeric, mustBePositive} = 5 % Mass of the mobile plate [kg]
args.Mph (1,1) double {mustBeNumeric, mustBePositive} = 10e-3 % Thickness of the mobile plate [m]
args.Mpr (1,1) double {mustBeNumeric, mustBePositive} = 150e-3 % Radius of the mobile plate [m]
%% initializeCylindricalStruts
args.Fsm (1,1) double {mustBeNumeric, mustBePositive} = 1e-3 % Mass of the fixed part of the strut [kg]
args.Fsh (1,1) double {mustBeNumeric, mustBePositive} = 60e-3 % Length of the fixed part of the struts [m]
args.Fsr (1,1) double {mustBeNumeric, mustBePositive} = 5e-3 % Radius of the fixed part of the struts [m]
args.Msm (1,1) double {mustBeNumeric, mustBePositive} = 1e-3 % Mass of the mobile part of the strut [kg]
args.Msh (1,1) double {mustBeNumeric, mustBePositive} = 60e-3 % Length of the mobile part of the struts [m]
args.Msr (1,1) double {mustBeNumeric, mustBePositive} = 5e-3 % Radius of the fixed part of the struts [m]
%% Bottom and Top Flexible Joints
args.flex_type_F char {mustBeMember(args.flex_type_F,{'2dof', '3dof', '4dof', '2dof_axial', 'flexible'})} = '2dof'
args.flex_type_M char {mustBeMember(args.flex_type_M,{'2dof', '3dof', '4dof', '2dof_axial', 'flexible'})} = '3dof'
args.Kf_M (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Cf_M (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Kt_M (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Ct_M (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Kf_F (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Cf_F (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Kt_F (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Ct_F (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Ka_F (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Ca_F (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Kr_F (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Cr_F (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Ka_M (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Ca_M (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Kr_M (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Cr_M (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
%% inverseKinematics
args.AP (3,1) double {mustBeNumeric} = zeros(3,1)
args.ARB (3,3) double {mustBeNumeric} = eye(3)
end
stewart = initializeStewartPlatform();
switch args.type
case 'none'
stewart.type = 0;
case 'stewart'
stewart.type = 1;
end
stewart = initializeFramesPositions(stewart, ...
'H', args.H, ...
'MO_B', args.MO_B);
stewart = generateGeneralConfiguration(stewart, ...
'FH', args.FH, ...
'FR', args.FR, ...
'FTh', args.FTh, ...
'MH', args.MH, ...
'MR', args.MR, ...
'MTh', args.MTh);
stewart = computeJointsPose(stewart);
if strcmp(args.actuator_type, 'apa300ml')
args.actuator_k = 0.3e6;
args.actuator_ke = 4.3e6;
args.actuator_ka = 2.15e6;
args.actuator_c = 18;
args.actuator_ce = 0.7;
args.actuator_ca = 0.35;
args.actuator_ga = -2.7;
args.actuator_gs = 0.53e6;
elseif strcmp(args.actuator_type, 'flexible')
args.actuator_ga = 25.9;
args.actuator_gs = -5.08e6;
end
stewart = initializeStrutDynamics(stewart, ...
'type', args.actuator_type, ...
'k', args.actuator_k, ...
'kp', args.actuator_kp, ...
'ke', args.actuator_ke, ...
'ka', args.actuator_ka, ...
'c', args.actuator_c, ...
'cp', args.actuator_cp, ...
'ce', args.actuator_ce, ...
'ca', args.actuator_ca, ...
'ga', args.actuator_ga, ...
'gs', args.actuator_gs);
stewart = initializeJointDynamics(stewart, ...
'type_F', args.flex_type_F, ...
'type_M', args.flex_type_M, ...
'Kf_M', args.Kf_M, ...
'Cf_M', args.Cf_M, ...
'Kt_M', args.Kt_M, ...
'Ct_M', args.Ct_M, ...
'Kf_F', args.Kf_F, ...
'Cf_F', args.Cf_F, ...
'Kt_F', args.Kt_F, ...
'Ct_F', args.Ct_F, ...
'Ka_F', args.Ka_F, ...
'Ca_F', args.Ca_F, ...
'Kr_F', args.Kr_F, ...
'Cr_F', args.Cr_F, ...
'Ka_M', args.Ka_M, ...
'Ca_M', args.Ca_M, ...
'Kr_M', args.Kr_M, ...
'Cr_M', args.Cr_M);
stewart = initializeCylindricalPlatforms(stewart, ...
'Fpm', args.Fpm, ...
'Fph', args.Fph, ...
'Fpr', args.Fpr, ...
'Mpm', args.Mpm, ...
'Mph', args.Mph, ...
'Mpr', args.Mpr);
stewart = initializeCylindricalStruts(stewart, ...
'Fsm', args.Fsm, ...
'Fsh', args.Fsh, ...
'Fsr', args.Fsr, ...
'Msm', args.Msm, ...
'Msh', args.Msh, ...
'Msr', args.Msr);
stewart = computeJacobian(stewart);
stewart = initializeStewartPose(stewart, ...
'AP', args.AP, ...
'ARB', args.ARB);
nano_hexapod = stewart;
if exist('./mat', 'dir')
if exist('./mat/nass_model_stages.mat', 'file')
save('mat/nass_model_stages.mat', 'nano_hexapod', '-append');
else
save('mat/nass_model_stages.mat', 'nano_hexapod');
end
elseif exist('./matlab', 'dir')
if exist('./matlab/mat/nass_model_stages.mat', 'file')
save('matlab/mat/nass_model_stages.mat', 'nano_hexapod', '-append');
else
save('matlab/mat/nass_model_stages.mat', 'nano_hexapod');
end
end
end
#+end_src
*** =initializeStrutDynamics=: Add Stiffness and Damping properties of each strut
#+begin_src matlab :tangle matlab/src/initializeStrutDynamics.m :comments none :mkdirp yes :eval no
function [stewart] = initializeStrutDynamics(stewart, args)
% initializeStrutDynamics - Add Stiffness and Damping properties of each strut
%
% Syntax: [stewart] = initializeStrutDynamics(args)
%
% Inputs:
% - args - Structure with the following fields:
% - K [6x1] - Stiffness of each strut [N/m]
% - C [6x1] - Damping of each strut [N/(m/s)]
%
% Outputs:
% - stewart - updated Stewart structure with the added fields:
% - actuators.type = 1
% - actuators.K [6x1] - Stiffness of each strut [N/m]
% - actuators.C [6x1] - Damping of each strut [N/(m/s)]
arguments
stewart
args.type char {mustBeMember(args.type,{'1dof', '2dof', 'flexible', 'apa300ml'})} = '1dof'
args.k (1,1) double {mustBeNumeric, mustBeNonnegative} = 20e6
args.kp (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.ke (1,1) double {mustBeNumeric, mustBeNonnegative} = 5e6
args.ka (1,1) double {mustBeNumeric, mustBeNonnegative} = 60e6
args.c (1,1) double {mustBeNumeric, mustBeNonnegative} = 2e1
args.cp (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.ce (1,1) double {mustBeNumeric, mustBeNonnegative} = 1e6
args.ca (1,1) double {mustBeNumeric, mustBeNonnegative} = 10
args.ga (1,1) double {mustBeNumeric} = 1
args.gs (1,1) double {mustBeNumeric} = 1
args.F_gain (1,1) double {mustBeNumeric} = 1
args.me (1,1) double {mustBeNumeric} = 0.01
args.ma (1,1) double {mustBeNumeric} = 0.01
end
if strcmp(args.type, '1dof')
stewart.actuators.type = 1;
elseif strcmp(args.type, '2dof')
stewart.actuators.type = 2;
elseif strcmp(args.type, 'flexible')
stewart.actuators.type = 3;
K = readmatrix('APA300ML_flex_mat_K.CSV');
M = readmatrix('APA300ML_flex_mat_M.CSV');
[int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('APA300ML_flex_out_nodes_3D.txt');
stewart.actuators.M = M;
stewart.actuators.K = K;
stewart.actuators.n_xyz = int_xyz;
stewart.actuators.xi = 0.05;
elseif strcmp(args.type, 'apa300ml')
stewart.actuators.type = 4;
end
stewart.actuators.k = args.k;
stewart.actuators.c = args.c;
% Parallel stiffness
stewart.actuators.kp = args.kp;
stewart.actuators.cp = args.cp;
stewart.actuators.ka = args.ka;
stewart.actuators.ca = args.ca;
stewart.actuators.ke = args.ke;
stewart.actuators.ce = args.ce;
stewart.actuators.ga = args.ga;
stewart.actuators.gs = args.gs;
stewart.actuators.F_gain = args.F_gain;
stewart.actuators.ma = args.ma;
stewart.actuators.me = args.me;
end
#+end_src
*** =initializeJointDynamics=: Add Stiffness and Damping properties for spherical joints
#+begin_src matlab :tangle matlab/src/initializeJointDynamics.m :comments none :mkdirp yes :eval no
function [stewart] = initializeJointDynamics(stewart, args)
% initializeJointDynamics - Add Stiffness and Damping properties for the spherical joints
%
% Syntax: [stewart] = initializeJointDynamics(args)
%
% Inputs:
% - args - Structure with the following fields:
% - type_F - 'universal', 'spherical', 'universal_p', 'spherical_p'
% - type_M - 'universal', 'spherical', 'universal_p', 'spherical_p'
% - Kf_M [6x1] - Bending (Rx, Ry) Stiffness for each top joints [(N.m)/rad]
% - Kt_M [6x1] - Torsion (Rz) Stiffness for each top joints [(N.m)/rad]
% - Cf_M [6x1] - Bending (Rx, Ry) Damping of each top joint [(N.m)/(rad/s)]
% - Ct_M [6x1] - Torsion (Rz) Damping of each top joint [(N.m)/(rad/s)]
% - Kf_F [6x1] - Bending (Rx, Ry) Stiffness for each bottom joints [(N.m)/rad]
% - Kt_F [6x1] - Torsion (Rz) Stiffness for each bottom joints [(N.m)/rad]
% - Cf_F [6x1] - Bending (Rx, Ry) Damping of each bottom joint [(N.m)/(rad/s)]
% - Cf_F [6x1] - Torsion (Rz) Damping of each bottom joint [(N.m)/(rad/s)]
%
% Outputs:
% - stewart - updated Stewart structure with the added fields:
% - stewart.joints_F and stewart.joints_M:
% - type - 1 (universal), 2 (spherical), 3 (universal perfect), 4 (spherical perfect)
% - Kx, Ky, Kz [6x1] - Translation (Tx, Ty, Tz) Stiffness [N/m]
% - Kf [6x1] - Flexion (Rx, Ry) Stiffness [(N.m)/rad]
% - Kt [6x1] - Torsion (Rz) Stiffness [(N.m)/rad]
% - Cx, Cy, Cz [6x1] - Translation (Rx, Ry) Damping [N/(m/s)]
% - Cf [6x1] - Flexion (Rx, Ry) Damping [(N.m)/(rad/s)]
% - Cb [6x1] - Torsion (Rz) Damping [(N.m)/(rad/s)]
arguments
stewart
args.type_F char {mustBeMember(args.type_F,{'2dof', '3dof', '4dof', '2dof_axial', 'flexible'})} = '2dof'
args.type_M char {mustBeMember(args.type_M,{'2dof', '3dof', '4dof', '2dof_axial', 'flexible'})} = '3dof'
args.Kf_M (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Cf_M (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Kt_M (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Ct_M (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Kf_F (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Cf_F (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Kt_F (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Ct_F (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Ka_F (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Ca_F (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Kr_F (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Cr_F (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Ka_M (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Ca_M (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Kr_M (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.Cr_M (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
args.K_M double {mustBeNumeric} = zeros(6,6)
args.M_M double {mustBeNumeric} = zeros(6,6)
args.n_xyz_M double {mustBeNumeric} = zeros(2,3)
args.xi_M double {mustBeNumeric} = 0.1
args.step_file_M char {} = 'flexor_025.STEP'
args.K_F double {mustBeNumeric} = zeros(6,6)
args.M_F double {mustBeNumeric} = zeros(6,6)
args.n_xyz_F double {mustBeNumeric} = zeros(2,3)
args.xi_F double {mustBeNumeric} = 0.1
args.step_file_F char {} = 'flexor_025.STEP'
end
switch args.type_F
case '2dof'
stewart.joints_F.type = 1;
case '3dof'
stewart.joints_F.type = 2;
case '4dof'
stewart.joints_F.type = 3;
case '2dof_axial'
stewart.joints_F.type = 4;
case 'flexible'
stewart.joints_F.type = 5;
K = readmatrix('flex025_mat_K.CSV');
M = readmatrix('flex025_mat_M.CSV');
[int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('flex025_out_nodes_3D.txt');
stewart.joints_F.M = M;
stewart.joints_F.K = K;
stewart.joints_F.n_xyz = int_xyz;
stewart.joints_F.xi = 0.05;
stewart.joints_F.step_file = args.step_file_F;
otherwise
error("joints_F are not correctly defined")
end
switch args.type_M
case '2dof'
stewart.joints_M.type = 1;
case '3dof'
stewart.joints_M.type = 2;
case '4dof'
stewart.joints_M.type = 3;
case '2dof_axial'
stewart.joints_M.type = 4;
case 'flexible'
stewart.joints_M.type = 5;
K = readmatrix('flex025_mat_K.CSV');
M = readmatrix('flex025_mat_M.CSV');
[int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('flex025_out_nodes_3D.txt');
stewart.joints_M.M = M;
stewart.joints_M.K = K;
stewart.joints_M.n_xyz = int_xyz;
stewart.joints_M.xi = 0.05;
stewart.joints_M.step_file = args.step_file_M;
otherwise
error("joints_M are not correctly defined")
end
stewart.joints_M.Ka = args.Ka_M;
stewart.joints_M.Kr = args.Kr_M;
stewart.joints_F.Ka = args.Ka_F;
stewart.joints_F.Kr = args.Kr_F;
stewart.joints_M.Ca = args.Ca_M;
stewart.joints_M.Cr = args.Cr_M;
stewart.joints_F.Ca = args.Ca_F;
stewart.joints_F.Cr = args.Cr_F;
stewart.joints_M.Kf = args.Kf_M;
stewart.joints_M.Kt = args.Kt_M;
stewart.joints_F.Kf = args.Kf_F;
stewart.joints_F.Kt = args.Kt_F;
stewart.joints_M.Cf = args.Cf_M;
stewart.joints_M.Ct = args.Ct_M;
stewart.joints_F.Cf = args.Cf_F;
stewart.joints_F.Ct = args.Ct_F;
end
#+end_src
* Footnotes
[fn:1]The manufacturer of the APA95ML was not willing to share the piezoelectric material properties of the stack.