NASS - Finite Element Models with Simscape

We first extract the stiffness and mass matrices.

K = readmatrix('APA300ML_mat_K.CSV');
M = readmatrix('APA300ML_mat_M.CSV');

Then, we extract the coordinates of the interface nodes.

[int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('APA300ML_out_nodes_3D.txt');

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

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

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

The ratio of the developed force to applied voltage is:

where:

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

If we take the numerical values, we obtain:

d33*n*ka % [N/V]

5.64

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

where:

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

If we take the numerical values, we obtain:

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

23.438

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

Let’s say we use two stacks as a force sensor and one stack as an actuator:

• A Relative Motion Sensor block is added between the nodes A and C
• An Internal Force block is added between the remote points E and B

The interface nodes are shown in Figure 1.

One mass is fixed at one end of the piezo-electric stack actuator (remove point F), the other end is fixed to the world frame (remote point G).

We first set the mass to be approximately zero. The dynamics is identified from the applied force to the measured relative displacement. The same dynamics is identified for a payload mass of 10Kg.

m = 10;

The root locus corresponding to Direct Velocity Feedback with a mass of 10kg is shown in Figure 4.

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

For root locus corresponding to IFF is shown in Figure 6.

The goal in this section is to identify the parameters of a simple APA model from the FEM. This can be useful is a lower order model is to be used for simulations.

The presented model is based on (Souleille et al. 2018).

The model represents the Amplified Piezo Actuator (APA) from Cedrat-Technologies (Figure 7). The parameters are shown in the table below.

The goal is to determine \(k_e\), \(k_a\) and \(k_1\) so that the simplified model fits the FEM model.

\[ \alpha = \frac{x_1}{f}(\omega=0) = \frac{\frac{k_e}{k_e + k_a}}{k_1 + \frac{k_e k_a}{k_e + k_a}} \] \[ \beta = \frac{x_1}{F}(\omega=0) = \frac{1}{k_1 + \frac{k_e k_a}{k_e + k_a}} \]

If we can fix \(k_a\), we can determine \(k_e\) and \(k_1\) with: \[ k_e = \frac{k_a}{\frac{\beta}{\alpha} - 1} \] \[ k_1 = \frac{1}{\beta} - \frac{k_e k_a}{k_e + k_a} \]

From the identified dynamics, compute \(\alpha\) and \(\beta\)

alpha = abs(dcgain(G('y', 'Fa')));
beta  = abs(dcgain(G('y', 'Fd')));

\(k_a\) is estimated using the following formula:

ka = 0.8/abs(dcgain(G('y', 'Fa')));

The factor can be adjusted to better match the curves.

Then \(k_e\) and \(k_1\) are computed.

ke = ka/(beta/alpha - 1);
k1 = 1/beta - ke*ka/(ke + ka);

The damping in the system is adjusted to match the FEM model if necessary.

c1 = 1e2;

The analytical model of the simpler system is defined below:

Ga = 1/(m*s^2 + k1 + c1*s + ke*ka/(ke + ka)) * ...
     [ 1 ,              k1 + c1*s + ke*ka/(ke + ka)  , ke/(ke + ka) ;
      -ke*ka/(ke + ka), ke*ka/(ke + ka)*m*s^2 ,       -ke/(ke + ka)*(m*s^2 + c1*s + k1)];

Ga.InputName = {'Fd', 'w', 'Fa'};
Ga.OutputName = {'y', 'Fs'};

And the DC gain is adjusted for the force sensor:

F_gain = dcgain(G('Fs', 'Fd'))/dcgain(Ga('Fs', 'Fd'));

The dynamics of the FEM model and the simpler model are compared in Figure 8.

The simplified model has also been implemented in Simscape.

The dynamics of the Simscape simplified model is identified and compared with the FEM one in Figure 9.

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). The payload mass is set to 10kg.

m = 10;

The obtained dynamics is shown in Figure 10.

The controller is defined below and the loop gain is shown in Figure 11.

Kiff = -1e3/s;

Now the closed-loop system is identified again and compare with the open loop system in Figure 12.

It is the expected behavior as shown in the Figure 13 (from (Souleille et al. 2018)).

We first extract the stiffness and mass matrices.

K = extractMatrix('mat_K_6modes_2MDoF.matrix');
M = extractMatrix('mat_M_6modes_2MDoF.matrix');

Then, we extract the coordinates of the interface nodes.

[int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('out_nodes_3D.txt');

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

The flexor is now imported into Simscape and its parameters are estimated using an identification.

The dynamics is identified from the applied force/torque to the measured displacement/rotation of the flexor. And we find the same parameters as the one estimated from the Stiffness matrix.

Let’s now model the flexible joint with a “perfect” Bushing joint as shown in Figure 15.

The parameters of the Bushing joint (stiffnesses) are estimated from the Stiffness matrix that was computed from the FEM.

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]

The dynamics from the applied force/torque to the measured displacement/rotation of the flexor is identified again for this simpler model. The two obtained dynamics are compared in Figure

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');

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

The flexor is now imported into Simscape and its parameters are estimated using an identification.

The dynamics is identified from the applied force/torque to the measured displacement/rotation of the flexor. And we find the same parameters as the one estimated from the Stiffness matrix.

Let’s now model the flexible joint with a “perfect” Bushing joint as shown in Figure 15.

The parameters of the Bushing joint (stiffnesses) are estimated from the Stiffness matrix that was computed from the FEM.

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]

The dynamics from the applied force/torque to the measured displacement/rotation of the flexor is identified again for this simpler model. The two obtained dynamics are compared in Figure

The stiffness matrix with the flexible joint with a “hinge” size of 0.50mm is loaded.

K_050 = readmatrix('flex050_mat_K.CSV');

Its parameters are compared with the Flexible Joint with a size of 0.25mm in the table below.

Now, the full nano-hexapod strut is modelled using Ansys.

The 3D as well as the interface nodes are shown in Figure 20.

A side view is shown in Figure 21.

The flexible joints used have a 0.25mm width size.

We first extract the stiffness and mass matrices.

K = readmatrix('strut_encoder_mat_K.CSV');
M = readmatrix('strut_encoder_mat_M.CSV');

Then, we extract the coordinates of the interface nodes.

[int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('strut_encoder_out_nodes_3D.txt');

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

Parameters for the APA300ML:

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

na = 2; % Number of stacks used as actuator
ns = 1; % Number of stacks used as force sensor

The dynamics is identified from the applied force to the measured relative displacement. The same dynamics is identified for a payload mass of 10Kg.

m = 10;