tdehaeze/phd-nass-fem
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
phd-nass-fem/nass-fem.org

3872 lines
173 KiB
Org Mode
Raw Blame History

#+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: <link rel="stylesheet" type="text/css" href="https://research.tdehaeze.xyz/css/style.css"/>
#+HTML_HEAD: <script type="text/javascript" src="https://research.tdehaeze.xyz/js/script.js"></script>
#+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:
<<sec:flexible_strut_encoder>>
**** 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)
<<matlab-dir>>
#+end_src
#+begin_src matlab :exports none :results silent :noweb yes
<<matlab-init>>
#+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:
<<sec:flexible_flexor_ID16>>
**** 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)
<<matlab-dir>>
#+end_src
#+begin_src matlab :exports none :results silent :noweb yes
<<matlab-init>>
#+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
<<sec:detail_fem_super_element>>
** 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)
<<matlab-dir>>
#+end_src
#+begin_src matlab :exports none :results silent :noweb yes
<<matlab-init>>
#+end_src
#+begin_src matlab :tangle no :noweb yes
<<m-init-path>>
#+end_src
#+begin_src matlab :eval no :noweb yes
<<m-init-path-tangle>>
#+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
<<m-init-other>>
#+end_src
** Procedure
<<ssec:detail_fem_super_element_theory>>
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
<<ssec:detail_fem_super_element_example>>
**** 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
<<ssec:detail_fem_super_element_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(f<f_start):length(f)-sum(f>f_end));
f_id = f(sum(f<f_start):length(f)-sum(f>f_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 Selection
<<sec:detail_fem_actuator>>
** Introduction :ignore:
The selection and modeling of actuators constitutes a critical step in the development of the nano-hexapod.
This chapter presents the approach to actuator selection and modeling.
First, specifications for the nano-hexapod actuators are derived from previous analyses, leading to the selection of the actuator type and ultimately to a specific model (Section ref:ssec:detail_fem_actuator_specifications).
Then, the chosen actuator is modeled using the reduced-order flexible body approach developed in the previous section, enabling validation of this selection through detailed dynamical analysis (Section ref:ssec:detail_fem_actuator_apa300ml).
Finally, a simplified two-degree-of-freedom model is developed to facilitate time-domain simulations while maintaining accurate representation of the actuator's essential characteristics (Section ref:ssec:detail_fem_actuator_apa300ml_2dof).
# TODO Add link to other sections
** Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
<<matlab-dir>>
#+end_src
#+begin_src matlab :exports none :results silent :noweb yes
<<matlab-init>>
#+end_src
#+begin_src matlab :tangle no :noweb yes
<<m-init-path>>
#+end_src
#+begin_src matlab :eval no :noweb yes
<<m-init-path-tangle>>
#+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
<<m-init-other>>
#+end_src
** Choice of the Actuator based on Specifications
<<ssec:detail_fem_actuator_specifications>>
The actuator selection process was driven by several critical requirements derived from previous dynamic analyses.
A primary consideration is the actuator stiffness, which significantly impacts system dynamics through multiple mechanisms.
The spindle rotation induces gyroscopic effects that modify plant dynamics and increase coupling, necessitating sufficient stiffness.
Conversely, the actuator stiffness must be carefully limited to ensure the nano-hexapod's suspension modes remain below the problematic modes of the micro-stations to limit the coupling between the two structures.
These competing requirements suggest an optimal stiffness of approximately $1\,N/\mu m$.
Additional specifications arise from the control strategy and physical constraints.
The implementation of a HAC-LAC (High Authority Control-Low Authority Control) architecture necessitates integrated force sensing capability.
The system's geometric constraints limit the actuator height to 50mm, given the nano-hexapod's maximum height of 95mm and the presence of flexible joints at each strut extremity.
Furthermore, the actuator stroke must exceed the micro-station positioning errors while providing additional margin for mounting adjustments and operational flexibility, which is estimated at $\approx 100\,\mu m$.
Three actuator technologies were evaluated (examples are shown in Figure ref:fig:detail_fem_actuator_pictures): voice coil actuators, piezoelectric stack actuators, and amplified piezoelectric actuators.
Variable reluctance actuators were not considered despite their superior efficiency compared to voice coil actuators, as their inherent nonlinearity would introduce unnecessary control complexity.
#+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
Voice coil actuators (shown in Figure ref:fig:detail_fem_voice_coil_picture), when combined with flexure guides of wanted stiffness $\approx 1\,N/\mu m$, would require forces above $100\,N$ to achieve the specified $100\,\mu m$ displacement.
While these actuators offer excellent linearity and long strokes, the constant force requirement would result in significant steady-state current, leading to thermal loads that could compromise system stability.
Their advantages were not considered adapted for this application, diminishing their benefits relative to piezoelectric solutions.
Conventional piezoelectric stack actuators (shown in Figure ref:fig:detail_fem_piezo_picture) present two significant limitations for the current application.
Their stroke is inherently limited to approximately $0.1\,\%$ of their length, meaning that even with the maximum allowable height of $50\,mm$, the achievable stroke would only be $50\,\mu m$, insufficient for the application.
Additionally, their extremely high stiffness, typically around $100\,N/\mu m$, exceeds the desired specifications by two orders of magnitude.
Amplified Piezoelectric Actuators (APAs) emerged as the optimal solution by addressing these limitations through an specific mechanical design.
The incorporation of a shell structure serves multiple purposes: it provides mechanical amplification of the piezoelectric displacement, reduces the effective axial stiffness to more suitable levels for the application, and creates a compact vertical profile.
Furthermore, the multi-stack configuration enables one stack to be dedicated to force sensing, ensuring excellent collocation with the actuator stacks, a critical feature for implementing robust decentralized control strategies.
Moreover, using APA for active damping has been successfully demonstrated in similar applications [[cite:&hanieh03_activ_stewar]].
Several specific APA models were evaluated against the established specifications (Table ref:tab:detail_fem_piezo_act_models).
The APA300ML emerged as the optimal choice.
This selection was further reinforced by previous experience with APAs from the same manufacturer[fn:2], and particularly by the successful validation of the modeling methodology with a similar actuator (Section ref:ssec:detail_fem_super_element_example).
The demonstrated accuracy of the modeling approach for the APA95ML provides confidence in the reliable prediction of the APA300ML's dynamic characteristics, thereby supporting both the selection decision and subsequent dynamical analyses.
#+name: tab:detail_fem_piezo_act_models
#+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
<<ssec:detail_fem_actuator_apa300ml>>
The validation of the APA300ML started by incorporating a "reduced order flexible body" into the multi-body model as explained in Section ref:sec:detail_fem_super_element.
The FEA model was developed with particular attention to the placement of reference frames, as illustrated in Figure ref:fig:detail_fem_apa300ml_frames.
Seven distinct frames were defined, with blue frames designating the force sensor stack interfaces for strain measurement, red frames denoting the actuator stack interfaces for force application and green frames for connecting to other elements.
120 additional modes were added during the modal reduction for a total order of 162.
While this high order provides excellent accuracy for validation purposes, it proves computationally intensive for simulations.
#+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
The sensor and actuator "constants" ($g_s$ and $g_a$) derived in Section ref:ssec:detail_fem_super_element_example for the APA95ML were used for the APA300ML model, as both actuators employ identical piezoelectric stacks.
** Simpler 2DoF Model of the APA300ML
<<ssec:detail_fem_actuator_apa300ml_2dof>>
To facilitate efficient time-domain simulations while maintaining essential dynamic characteristics, a simplified two-degree-of-freedom model was developed, adapted from [[cite:&souleille18_concep_activ_mount_space_applic]].
This model, illustrated in Figure ref:fig:detail_fem_apa_2dof_model, comprises three components.
The mechanical shell is characterized by its axial stiffness $k_1$ and damping $c_1$.
The actuator is modelled with stiffness $k_a$ and damping $c_a$, incorporating a force source $f$.
This force is related to the applied voltage $V_a$ through the actuator constant $g_a$.
The sensor stack is modeled with stiffness $k_e$ and damping $c_e$, with its deformation $d_L$ being converted to the output voltage $V_s$ through the sensor sensitivity $g_s$.
#+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]]
While providing computational efficiency, this simplified model has inherent limitations.
It considers only axial behavior, treating the actuator as infinitely rigid in other directions.
Several physical characteristics are not explicitly represented, including the mechanical amplification factor and the actual stress the piezoelectric stacks.
Nevertheless, the model's primary advantage lies in its simplicity, adding only four states to the system model.
The model requires tuning of 8 parameters ($k_1$, $c_1$, $k_e$, $c_e$, $k_a$, $c_a$, $g_s$, and $g_a$) to match the dynamics extracted from the finite element analysis.
The shell parameters $k_1$ and $c_1$ were determined first through analysis of the zero in the $V_a$ to $V_s$ transfer function.
The physical interpretation of this zero can be understood through Root Locus analysis: as controller gain increases, the poles of a closed-loop system converge to the open-loop zeros.
In this context, the zero corresponds to the poles of the system with a theoretical infinite-gain controller that ensures zero force in the sensor stack.
This condition effectively represents the dynamics of an APA without the force sensor stack.
This physical interpretation enables straightforward parameter tuning: $k_1$ determines the frequency of the zero, while $c_1$ defines its damping characteristic.
The stack parameters ($k_a$, $c_a$, $k_e$, $c_e$) were then derived from the first pole of the $V_a$ to $y$ response.
Given that identical piezoelectric stacks are used for both sensing and actuation, the relationships $k_e = 2k_a$ and $c_e = 2c_a$ were enforced, reflecting the series configuration of the dual actuator stacks.
Finally, the sensitivities $g_s$ and $g_a$ were adjusted to match the DC gains of the respective transfer functions.
The resulting parameters, documented in Table ref:tab:detail_fem_apa300ml_2dof_parameters, yield dynamic behavior that closely matches the high-order finite element model, as demonstrated in Figure ref:fig:detail_fem_apa300ml_comp_fem_2dof_fem_2dof.
While higher-order modes and non-axial flexibility are not captured, the model accurately represents the fundamental dynamics within the operational frequency range.
#+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
<<ssec:detail_fem_actuator_apa300ml_electrical>>
The behavior of piezoelectric actuators is characterized by coupled constitutive equations that establish relationships between electrical properties (charges, voltages) and mechanical properties (stress, strain) [[cite:&schmidt20_desig_high_perfor_mechat_third_revis_edition, chapter 5.5]].
To evaluate the impact of electrical boundary conditions on the system dynamics, experimental measurements were conducted using the APA95ML, comparing the transfer function from $V_a$ to $y$ under two distinct configurations.
With the force sensor stack in open-circuit condition (analogous to voltage measurement with high input impedance) and in short-circuit condition (similar to charge measurement with low output impedance).
As demonstrated in Figure ref:fig:detail_fem_apa95ml_effect_electrical_boundaries, short-circuiting the force sensor stack results in a minor decrease in resonance frequency.
This relatively modest effect validates the simplifying assumption made in the model of the APA.
#+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 [m/N]'); set(gca, 'XTickLabel',[]);
hold off;
ylim([1e-6, 1e-4]);
leg = legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 1);
leg.ItemTokenSize(1) = 15;
ax2 = nexttile;
hold on;
plot(f, 180/pi*angle(G_oc), '-')
plot(f, 180/pi*angle(G_sc), '-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
yticks(-360:45:360);
ylim([-20, 200]);
axis padded 'auto x'
linkaxes([ax1,ax2], 'x');
xlim([100, 300]);
#+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]]
However, the electrical characteristics of the APA remain crucial for instrumentation design.
Proper consideration must be given to voltage amplifier specifications and force sensor signal conditioning requirements.
These aspects, being fundamental to system implementation, will be addressed in the instrumentation chapter.
** Validation with the Nano-Hexapod
<<ssec:detail_fem_actuator_apa300ml_validation>>
The integration of the APA300ML model within the nano-hexapod simulation framework served two validation objectives: to validate the APA300ML choice through analysis of system dynamics with APA modelled as flexible bodies, and to validate the simplified 2DoF model through comparative analysis with the full FEM implementation.
The dynamic characteristics predicted using the flexible body model align well with the design requirements established during the conceptual phase.
The dynamics from $\bm{u}$ to $\bm{V}_s$ exhibits the desired alternating pole-zero pattern (Figure ref:fig:detail_fem_actuator_fem_vs_perfect_hac_plant), a critical characteristic for implementing robust decentralized Integral Force Feedback.
Additionally, the model predicts no problematic high-frequency modes in the dynamics from $\bm{u}$ to $\bm{\epsilon}_{\mathcal{L}}$ (Figure ref:fig:detail_fem_actuator_fem_vs_perfect_iff_plant), maintaining consistency with earlier conceptual simulations.
These findings suggest that the control performance targets established during the conceptual phase remain achievable with the selected actuator.
Comparative analysis between the high-order FEM implementation and the simplified 2DoF model (Figure ref:fig:detail_fem_actuator_fem_vs_perfect_plants) demonstrates remarkable agreement in the frequency range of interest.
This validates the use of the simplified model for time-domain simulations, where computational efficiency is paramount.
The reduction in model order is substantial: while the FEM implementation results in approximately 300 states (36 states per actuator plus 12 additional states), the 2DoF model requires only 24 states for the complete nano-hexapod.
These results validate both the selection of the APA300ML and the effectiveness of the simplified modeling approach for the nano-hexapod.
#+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 Design
<<sec:detail_fem_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)
<<matlab-dir>>
#+end_src
#+begin_src matlab :exports none :results silent :noweb yes
<<matlab-init>>
#+end_src
#+begin_src matlab :tangle no :noweb yes
<<m-init-path>>
#+end_src
#+begin_src matlab :eval no :noweb yes
<<m-init-path-tangle>>
#+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
<<m-init-other>>
#+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
<<sec:joints_rot_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
<<sec:joints_trans_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:
<<sec:detail_fem_conclusion>>
* 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:2]Cedrat technologies
[fn:1]The manufacturer of the APA95ML was not willing to share the piezoelectric material properties of the stack.
Reference in New Issue View Git Blame Copy Permalink