Use of RecurDyn for the analysis of a compliant mechanism by modeling with techniques for the flexible multibody.

The following article is an extract from a project carried out for the course of Advanced Techniques for the design of Prosthetic Devices. Held for Medical Engineering at the University of Rome Tor Vergata .

Authors: Mastrofini Alessandro & Muscedere Erica


Table of Contents

Flexible multibody

The techniques for flexible multibody analysis make it possible to take into account the flexibility of the body within dynamic analyzes. In particular, this is essential when this flexibility acts on the inertial characteristics or is indispensable to motion, as in flexible hinges. Within RecurDyn we can use different formulations to take into account the flexibility of the bodies.

Discrete flexible multibody

The discrete flexible techniques take flexibility into account by segmenting the element into different segments. These segments are connected by elastic forces expressed by means of a stiffness matrix correlated to the elastic properties of the body.

These methods are more than valid for the simulations of beam elements as it is clear how to segment and the elastic forces can be deduced from the properties of the beam and its section.

Full flexible

In full flex techniques the entire element is made flexible and therefore all degrees of freedom of finite element discretization are included within the equation of motion.

A corotational approach is used within the RecuDyn simulation environment. Therefore, the nodal displacements are calculated by superimposing a deformation deriving from neighboring elements on the rigid motion. Therefore the internal elastic energy, due to tensions and deformations, is taken into account in the overall equation of motion.

Structure under consideration. You see the flexible element (right element) represented as a line in the technical drawing (a); as a discrete set of bodies of type Beam Group (b); flexible body intermingled with shell elements (c). The centers of mass that identify the different solid bodies and the rotoidal and interlocking torque constraints are also evident. In the figure there are no graphic markers for the elastic forces between the elements of the Beam Group.
Fig. 1: Structure under consideration. You see the flexible element (right element) represented as a line in the technical drawing (a); as a discrete set of bodies of type Beam Group (b); flexible body intermingled with shell elements (c). The centers of mass that identify the different solid bodies and the rotoidal and interlocking torque constraints are also evident. In the figure there are no graphic markers for the elastic forces between the elements of the Beam Group.

This approach allows to combine a not excessive computational burden and to manage a formulation similar to the finite elements that allows to directly calculate stresses and deformations.

Modeling with RecurDyn

To model the structure under analysis (fig. 1a) we chose to model the rigid bodies using CAD and import the parasolid file into RecurDyn. The flexible body was then modeled and constrained, making sure to keep the motion of the flat structure as much as possible.

For the discrete flexible we used a BeamG element with 40 segments (fig. 1b). The section properties were then set taking into account both the geometric and material properties for harmonic steel EN 10270-1 class SM. For the full-flexible, a mesh with triangular plate-like elements of the PShell1 type is used, taking into account, within the element, the material properties of harmonic steel. 5982 elements are considered (fig. 1c).

Clearly the full flex formulation requires a greater computational burden which, as can be seen from the following results, is not always justified but often the results are similar to the discrete flexible case.

Moving torque

Considering the rigid steel bodies, these can act more or less intensely on the force field depending on how the structure is positioned. We therefore report a case with the presence of gravity and the structure vertically and the analogous case in the absence of gravity for which, as can be seen from the following results, the forces acting on the flexible hinge are greatly reduced.

Torque applied with respect to the angle of rotation of the crank in the range [0โ—ฆ; 360โ—ฆ] with a rotation speed of โˆ’๐œ‹ โ‹… 10โˆ’4 rad / s, in the presence of gravity. Comparison between discrete flexible and full flex.
Fig. 2: Torque applied with respect to the angle of rotation of the crank in the range [0โ—ฆ; 360โ—ฆ] with a rotation speed of โˆ’๐œ‹ โ‹… 10โˆ’4 rad / s, in the presence of gravity. Comparison between discrete flexible and full flex.

To analyze the torque applied in static conditions, a speed was identified low enough to neglect the dynamic effects. In particular, starting from a speed of โˆ’๐œ‹ rad / s, this has been scaled by dividing by an increasing scale factor: ๐™ฟ๐š…๐šœ๐šŒ๐šŠ๐š•๐šŽ โˆˆ {1, 2; 5; 10; 102; 103; 104; 105}. We therefore chose the case in which the speed is equal to ๐™ฟ๐š…๐šŸ๐šŽ๐š• = โˆ’10โˆ’4 โ‹… ๐œ‹ rad / s in which the trend of the applied torque is constant and equivalent to the following case.
To consider the torque acting in the crank rotation range in [0โ—ฆ; 360โ—ฆ] a simulation was set on a time equal to ๐‘ก = 2 โ‹… ๐™ฟ๐š…_๐šœ๐šŒ๐šŠ๐š•๐šŽ so as to have a complete revolution (2๐œ‹), by appropriately scaling the integration step.

Gravity

In the case in which the force of gravity is considered, the weight of the rigid members, in steel, tends to crush the flexible hinge. Furthermore, the presence of these forces tends to generate a motion in the opposite direction to that imposed by the (negative) speed and initially a discordant torque is observed (fig. 2).

Figura 3: Discrete flex - g - Andamento della coppia applicata al membro movente nel range [0; 2๐œ‹] rad al variare del parametro di scala (e quindi della velocitร ) (a); zoom sul caso ๐™ฟ๐š…_๐šŸ๐šŽ๐š• โˆˆ {๐œ‹โˆ•2; ๐œ‹ โ‹… 10โˆ’2; ๐œ‹ โ‹… 10โˆ’5} con limite a ยฑ2000 N/mm, in presenza della forza di gravitร . Si osservi come nei casi a velocitร  piรน alta sono presenti oscillazioni piรน o meno grandi mentre riducendo la velocitร  si arriva ad una curve sempre piรน stabile e costante. In particolare, le curve relative a velocitร  di ๐œ‹ โ‹… 10โˆ’4 e ๐œ‹ โ‹… 10โˆ’5 rad/s coincidono.
Fig. 3Discrete flex – g – Trend of the torque applied to the driving member in the range [0; 2๐œ‹] rad as the scale parameter (and therefore the speed) (a) varies; zoom on the case ๐™ฟ๐š…_๐šŸ๐šŽ๐š• โˆˆ {๐œ‹ โˆ• 2; ๐œ‹ โ‹… 10โˆ’2; ๐œ‹ โ‹… 10โˆ’5} with a limit of ยฑ 2000 N / mm, in the presence of gravity. It should be observed that in the cases at higher speed there are more or less large oscillations while reducing the speed one arrives at an increasingly stable and constant curve. In particular, the curves relating to velocities of ๐œ‹ โ‹… 10โˆ’4 and ๐œ‹ โ‹… 10โˆ’5 rad / s coincide.

Inertial contribution

Comparing the case of discrete with the full flex case, it is possible to see how the maximum values, in modulus, are similar and in both cases a crossing of zero occurs. However, the two trends show a phase shift even though the speed set is the same (โˆ’๐‘๐‘– โ‹… 10โˆ’4 rad / s) and is such that the system can be considered as static. Being a quasi-static analysis we are not considering the inertial contribution of the bodies (neither rigid nor flexible) but we can see, even by comparing with the case in which gravity is not present (fig. 7), how the static contribution of the flexible hinge is very small compared to the weight of the rigid members.

Figura 4: Full flex - g - Andamento della coppia applicata al membro movente nel range [0; 2๐œ‹] rad al variare del parametro di scala (e quindi della velocitร ) (a); zoom sul caso ๐™ฟ๐š…_๐šŸ๐šŽ๐š• โˆˆ {๐œ‹โˆ•2; ๐œ‹ โ‹… 10โˆ’1; ๐œ‹ โ‹… 10โˆ’4} con limite a ยฑ2000 N/mm (b), in presenza della forza di gravitร . Si osservi come nei casi a velocitร  piรน alta sono presenti oscillazioni piรน o meno grandi mentre riducendo la velocitร  si arriva ad una curve sempre piรน stabile e costante.
Fig. 4Full flex – g – Trend of the torque applied to the moving member in the range [0; 2๐œ‹] rad as the scale parameter (and therefore the speed) (a) varies; zoom on the case ๐™ฟ๐š…_๐šŸ๐šŽ๐š• โˆˆ {๐œ‹ โˆ• 2; ๐œ‹ โ‹… 10โˆ’1; ๐œ‹ โ‹… 10โˆ’4} with a limit of ยฑ 2000 N / mm (b), in the presence of the force of gravity. It should be observed that in the cases at higher speed there are more or less large oscillations while reducing the speed one arrives at an increasingly stable and constant curve.

No gravity

If the force of gravity is not considered, it is possible to have a better estimate of the effect of the flexible hinge.

Figura 5: Discrete flex - no g - Andamento della coppia applicata al membro movente nel range [0; 2๐œ‹] rad al variare del parametro di scala (e quindi della velocitร ) (a); zoom sul caso ๐™ฟ๐š…_๐šŸ๐šŽ๐š• โˆˆ {๐œ‹โˆ•2; ๐œ‹ โ‹… 10โˆ’2; ๐œ‹ โ‹… 10โˆ’5} con limite a ยฑ2000 N/mm (b), in assenza della forza di gravitร . Si osservi come nei casi a velocitร  piรน alta sono presenti oscillazioni piรน o meno grandi mentre riducendo la velocitร  si arriva ad una curve sempre piรน stabile e costante. In particolare, le curve relative a velocitร  di ๐œ‹ โ‹… 10โˆ’4 e ๐œ‹ โ‹… 10โˆ’5 rad/s coincidono.
Fig. 5Discrete flex – no g – Trend of the torque applied to the moving member in the range [0; 2๐œ‹] rad as the scale parameter (and therefore the speed) (a) varies; zoom on the case ๐™ฟ๐š…_๐šŸ๐šŽ๐š• โˆˆ {๐œ‹ โˆ• 2; ๐œ‹ โ‹… 10โˆ’2; ๐œ‹ โ‹… 10โˆ’5} with a limit of ยฑ 2000 N / mm (b), in the absence of the force of gravity. It should be observed that in the cases at higher speed there are more or less large oscillations while reducing the speed one arrives at an increasingly stable and constant curve. In particular, the curves relating to velocities of ๐œ‹ โ‹… 10โˆ’4 and ๐œ‹ โ‹… 10โˆ’5 rad / s coincide.
Figura 6: Full flex - no g - Andamento della coppia applicata al membro movente nel range [0; 2๐œ‹] rad al variare del parametro di scala (e quindi della velocitร ) (a); zoom sul caso ๐™ฟ๐š…_๐šŸ๐šŽ๐š• โˆˆ {๐œ‹โˆ•2; ๐œ‹ โ‹… 10โˆ’1; ๐œ‹ โ‹… 10โˆ’4} con limite a ยฑ2000 N/mm (b), in assenza della forza di gravitร . Si osservi come nei casi a velocitร  piรน alta sono presenti oscillazioni piรน o meno grandi mentre riducendo la velocitร  si arriva ad una curve sempre piรน stabile e costante.
Fig 6Full flex – no g – Trend of the torque applied to the moving member in the range [0; 2๐œ‹] rad as the scale parameter (and therefore the speed) (a) varies; zoom on the case ๐™ฟ๐š…_๐šŸ๐šŽ๐š• โˆˆ {๐œ‹ โˆ• 2; ๐œ‹ โ‹… 10โˆ’1; ๐œ‹ โ‹… 10โˆ’4} with a limit of ยฑ 2000 N / mm (b), in the absence of the force of gravity. It should be observed that in the cases at higher speed there are more or less large oscillations while reducing the speed one arrives at an increasingly stable and constant curve.

Torque analysis

In particular, placing ourselves in a quasi-static condition we observe (fig. 7) an initial torque in agreement with the imposed speed. Subsequently, the torque changes sign, indicating that the static contribution of the hinge tends to provide a more intense elastic return than would be sufficient for the imposed motion. The maximum and minimum values are considerably lower than the case in which gravity is present.

Torque applied with respect to the angle of rotation of the crank in the range [0โ—ฆ; 360โ—ฆ] with a rotation speed of โˆ’๐œ‹ โ‹… 10โˆ’4 rad / s, in the absence of gravity. Comparison between discrete flexible and full flex.
Fig. 7: Torque applied with respect to the angle of rotation of the crank in the range [0โ—ฆ; 360โ—ฆ] with a rotation speed of โˆ’๐œ‹ โ‹… 10โˆ’4 rad / s, in the absence of gravity. Comparison between discrete flexible and full flex.

Driving torque

The trend of the drive torque increases as the speed set increases. So, the numerical results are present in table 1 and have been read from the graph of the driving torque acting on the crank, as soon as the motion starts.

By considering the traces of the center of the rotoidal pair between the intermediate member and the flexible member in fig. 8 and comparing the full flex cases with the discrete flex. The trends are very similar. In all cases there is at least one oscillation around the starting point, more or less pronounced.

Path of the center of the rotoidal pair between the intermediate member and the flexible element. Observe how the scale is double on the y axis and centered with respect to the origin of the joint, to give a better representation of what is happening. The cases with speed of 90 ยฐ / s are shown (a); 60 ยฐ / s (b); 40 ยฐ (s) (c) and 20 ยฐ / s (d). Also shown is a part of the intermediate element in the initial configuration.
Fig. 8: Path of the center of the rotoidal pair between the intermediate member and the flexible element. Observe how the scale is double on the y axis and centered with respect to the origin of the joint, to give a better representation of what is happening. The cases with speed of 90 ยฐ / s are shown (a); 60 ยฐ / s (b); 40 ยฐ (s) (c) and 20 ยฐ / s (d). Also shown is a part of the intermediate element in the initial configuration.
Deformation of the full flex structure without gravity in the initial configuration (a); maximum flexion at about 135 ยฐ (b) and maximum extension at about 315 ยฐ (c)circa 315ยฐ (c)
Fig. 9: Deformation of the full flex structure without gravity in the initial configuration (a); maximum flexion at about 135 ยฐ (b) and maximum extension at about 315 ยฐ (c)

Stress

For the full flex model it is easy to read the stress simply by activating the contour function and calculating the maximum value. A Von Mises criterion was chosen by reading the maximum values โ€‹โ€‹for the 4 different speeds considered. These values, in table 2, are reached in the moment of maximum bending, that is in the rotation range between 90 ยฐ and 180 ยฐ. An example has been reported for the case of speed 20 ยฐ / s in fig. 10.

For the discrete flexible case, on the other hand, it is necessary to do some more reasoning. First of all we are talking about rigid bodies so it is impossible to go to read deformations and / or stresses from the simulation. However, it is possible to read the applied loads and report them on a beam-like structure, referring to the single element of the beam group. To do this it is necessary to make some practical considerations.

The most loaded element is the wedged one so we can take the reactions and bring them back to the element, considering it as a deformable material characterized by a certain section ๐ด with its own moment of inertia ๐ผ๐‘–๐‘– , on the relative references.

Eulero-Bernoulli beam

Therefore, looking at the element as an elementary segment of a macroscopically slender beam, we can focus on the numerically significant contributions, i.e. the main stresses along the beam axis. Such stresses along ๐‘ฆ in the global reference system will hear of two calculable contributions using the Navier formula.

A contribution made by the tensile load:

\begin{equation}
\sigma_{y}=\frac{N}{A}
\end{equation}

And another one given by the acting couple:

\begin{equation}
\sigma_{y}=\frac{M_{z}}{I_{z z}} l_{x}
\end{equation}

Where ๐‘™๐‘ฅ is half the length in the third direction (i.e. half the thickness).

These two contributions are also the most relevant in the calculation of Von Mises stresses. This is due to the particular geometry of the motion which essentially, due to the way the structure is built, can be said to be flat. Furthermore, given the geometry of the beam and the distribution of the acting loads, we can consider a plane state of tension. It results in fact that the calculated values are comparable with what obtained with the full flex formulation (table 2), with the exception of the 90 ยฐ / s case where the gap is considerably wider.

Maximum stresses

For comparison, the cases in which the flexion is maximum are also reported in the case in which gravity is considered in fig. 10b. From these deformations we have excluded this case from these latest analyzes as some deformations are really excessive and probably attributable only to the numerical analysis.

Plot of the Von Mises stresses for the full flex structure with speed of 20 ยฐ / s in the configuration where the stress is maximum, about 110 ยฐ of rotation (a); deformed with the presence of gravity in the full flex simulation (b) with 90 ยฐ rotation (above) and 320 ยฐ (below)rotazione (a); deformate con presenza di forza di gravitร  nella simulazione full flex (b) con 90ยฐ di rotazione (sopra) e 320ยฐ (sotto)
Fig. 10: Plot of the Von Mises stresses for the full flex structure with speed of 20 ยฐ / s in the configuration where the stress is maximum, about 110 ยฐ of rotation (a); deformed with the presence of gravity in the full flex simulation (b) with 90 ยฐ rotation (above) and 320 ยฐ (below)
fulldiscrete
90 ยฐ/s-61.56-1155.940
60 ยฐ/s-27.36-756.94
40 ยฐ/s-12.16-498.54
20 ยฐ/s-3.04-246.23
Tab 1: In N mm. Drive torques tabulated with respect to the speed set and the formulation used.
20 ยฐ/s40 ยฐ/s60 ยฐ/s90 ยฐ/s
discrete76.6078.2580.30115.90
full78.8481.4183.47298.23
Tab 2: Resulting stresses in the discrete case (calculated) and in the full flex case (Von Mises). Values in MPa, case no ๐‘”โƒ—.

Frequency of the mechanism

For the analysis of the natural frequency of the mechanism, the RFlex analysis can be used within RecurDyn. In particular, once the translation modes have been removed, the first frequency results at 6.50 Hz and the second at 17.92 Hz. Further details are present in fig. 11.

Figura 11: Primi 5 modi di buckling (dopo i moti rigidi) a frequenza 6.50 Hz (a), 17.92 Hz (b), 35.13 Hz (c), 58.09 Hz (d) e 86.81 (e)
Fig 11: Top 5 buckling modes (after rigid motions) at frequency 6.50 Hz (a), 17.92 Hz (b), 35.13 Hz (c), 58.09 Hz (d) and 86.81 (e)

For the free oscillation frequency of the mechanism, which clearly cannot be analyzed using methods that make use of linear compositions in the regime of small displacements and small displacement gradients, it is necessary to make a higher level discussion. Probably this frequency is linked to the type of system considered, i.e. we are considering a set of masses connected by forces (stiffness matrices or internal elastic energy).

This is nothing more than a damped mass-spring system which, however distributed the masses and individual springs may be, could characterize a harmonic oscillator. Within the simulations, however, a damping factor is introduced, through the damping ratio considered equal to 10โˆ’3. The introduction of this damping tends to stabilize the system but tends to strongly reduce the oscillatory contributions, especially if it is raised.

This frequency can be analyzed by setting a speed condition with a frequency sweep and going to see, through the FFT, the frequency where there seems to be resonance.