Que Para

Obtener El

DOCTOR EN Grado De

Ingeniería Mecánica

Presenta Juan

Emmanuel Ayala

Dr J Jesús Sébastien Briot

Sánchez Cervantes

MÉXICO Gto Salamanca

Junio

UNIVERSIDAD DE GUANAJUATO CAMPUS IRAPUATO-SALAMANCA DIVISIÓN DE INGENIERÍAS TOPOLOGY OPTIMIZATION OF REACTIONLESS ROBOTS TESIS

Keywords: dynamic balancing, reactionless robots, multibody topology optimization, reactionless four-bar linkage, partial dynamic balancing of five-bar mechanism ABS Acrylonitrile butadiene styrene. AM Additive manufacturing. API Application program interface. AAN Aritficial neural networks. BESO Bidirectional evolutionary structural optimization CPU Central processing unit. MUBOTO Computational platform for multibody topology optimization. CAD Computer-aided design. CAM Computer-aided manufacturing. CNC Computerized numerical control

Dynamic balancing is an important field of study in high-speed robotics and spatial robots. Taking into account robot dynamic balancing performance for robot design leads to low base vibrations, high precision and short cycle times. With the aim to develop a comprehensive robot design for dynamic balancing, structural topology optimization is studied in this research work as a tool for designing dynamically balanced robots, also called reactionless robots.

Classical dynamic balancing techniques do not consider the linkage elastic performance, which is of the utmost importance. This research deals with a novel approach that overcomes this shortcoming through topology optimization, which is used to perform an optimal robot link design.

Topology optimization is a mathematical method that aims to redistribute the material into an initial domain taking into account design specifications. In our research, the links are treated as three-dimensional flexible bodies, and the optimization process is performed for all the bodies simultaneously, while the dynamic balancing conditions are fulfilled.

Considering the robot as a flexible multibody system allows to obtain an optimized design based on the structural performance indices while dynamic balancing conditions are considered. These structural performance indices can be the natural vibration response, static deformations and robot stiffness, among others. Topology optimization of multibody systems, in some cases, leads to face large-scale analysis requiring huge computational effort. Therefore a parallel computational platform was specifically developed in order to optimize flexible multibody systems under dynamic balancing constraints.

The suitability of the proposed methodology is confirmed by accomplishing an optimized design of a reactionless four-bar linkage and the partial dynamic balancing of five-bar robotic mechanism. The significance of the dynamically balanced four-bar linkage is related to the possibility to exploit this optimized linkage as a special leg for building reactionless robots. Besides, the five-bar robot is very important due to its industrial applications, where it is typically used in pick-and-place operations. Even if the five-bar mechanism is partially balanced, the shaking force elimination and the shaking moment reduction can impact positively in the mechanism diminishing the vibrations and increasing the position accuracy. i Numerical validations of the optimized linkage properties were carried out using commercial software. The dynamic balancing performance of the optimized mechanisms was numerically validated using ADAMS. Besides, ANSYS software was used to evaluate the optimized structural properties.

Finally, it is important to mention that a prototype of the optimized reactionless four-bar linkage was built in order to verify the feasibility of the proposed methodology for dynamic balancing. The prototype was used to conduct experimental tests to evaluate its dynamic balancing performance.

Résumé

L'équilibrage dynamique est un domaine d'étude important en robotique à grande vitesse et pour la robotique spatiale. La prise en compte des performances d'équilibrage dynamique des robots lors de leur conception permet de réduire les vibrations du bâti, une meilleure précision et des temps de cycle réduits. Dans le but de développer une méthode de conception de robot équilibrés dynamiquement, l'optimisation topologique structurelle est étudiée dans ce travail de recherche en tant qu'outil pour concevoir des robots équilibrés dynamiquement, c'est à dire un robot ne transmettant pas de réaction sur le châssis.

Les techniques classiques d'équilibrage dynamique ne tiennent pas compte de l'élasticité des corps, ce qui est crucial pour les performances en robotique. Cette recherche porte sur une nouvelle approche qui permet de surmonter cette lacune grâce à l'optimisation topologique, qui est utilisée pour concevoir de manière optimale les corps du robot. L'optimisation topologique est une méthode qui vise à redistribuer la matière dans un espace de recherche donné en tenant compte des spécifications de conception. Dans notre étude, les corps du robot sont traités comme des corps flexibles tridimensionnels, et le processus d'optimisation est effectué pour tous les corps simultanément, tout en respectant les conditions d'équilibre dynamique.

Considérer le robot comme un système multi-corps flexible permet de réaliser une conception optimisée basée sur des indices de performance structurelle tout en respectant les conditions d'équilibre dynamique. Ces indices de performance structurelle peuvent être, entre autres, la réponse aux vibrations, ou les déformations statiques du robot. L'optimisation topologique des systèmes multi-corps nécessite souvent un effort de calcul considérable. C'est pourquoi une plateforme de calcul parallèle a été spécifiquement développée pour réaliser ces optimisations de la conception.

La pertinence de la méthodologie proposée est confirmée par la réalisation d'une conception optimisée d'un mécanisme à quatre barres équilibré dynamiquement et l'équilibrage dynamique partiel d'un mécanisme à cinq barres. L'importance de la conception d'un mécanisme à quatre barres équilibré dynamiquement est liée à la possibilité d'exploiter ce mécanisme optimisé comme un module utile à la conception de robots équilibrés dynamiquement. Le mécanisme à cinq barres a iii une importance toute particulière en raison de ses applications industrielles, où il est généralement utilisé dans des opérations de prise et de dépose. Même si le mécanisme à cinq barres est partiellement équilibré, l'élimination des forces et moments dynamiques peuvent avoir un impact positif sur le mécanisme car ils permettent la diminution des vibrations et l'augmentation de la précision.

Des validations numériques des propriétés des corps optimisés ont été effectuées à l'aide de logiciels de simulations. La performance d'équilibrage dynamique des mécanismes optimisés a été validée numériquement en utilisant ADAMS. En outre, le logiciel ANSYS a été utilisé pour évaluer les propriétés élastiques.

Enfin, il est important de mentionner qu'un prototype du mécanisme optimisé à quatre barres équilibré dynamiquement a été construit afin de vérifier la faisabilité de la méthodologie proposée pour l'équilibrage dynamique. Le prototype a été utilisé pour effectuer des tests expérimentaux afin d'évaluer ses performances d'équilibrage dynamique.

Mots-clés: équilibrage dynamique, robots équilibré dynamiquement, optimisation topologique multi-corps, mécanisme à quatre barres équilibré dynamiquement, équilibrage dynamique partiel d'un mécanisme à cinq barres.

Resumen

El balanceo dinámico es un campo de estudio muy importante en aplicaciones de robots de alta velocidad y robots espaciales. Considerar el balanceo dinámico en el diseño de los robots nos conduce a reducir las vibraciones en la base, incrementar la precisión y reducir los ciclos de trabajo. Con el objetivo de realizar un diseño integral para el balanceo dinámico de robots, en este trabajo de investigación la optimización topológica estructural es estudiada como una herramienta para el diseño de robots dinámicamente balanceados, los cuales no transmiten ninguna carga reactiva desbalanceada sobre su base.

Las técnicas clásicas de balanceo dinámico no consideran el desempeño elástico de los mecanismos, lo que es de gran importancia. Esta investigación trata con un enfoque novedoso que atiende este problema a través de la optimización topológica, la cual es utilizada para realizar el diseño óptimo de los eslabones del robot. La optimización topológica es un método matemático cuyo objetivo es la distribución óptima de material dentro de un espacio inicial, tomando en cuenta especificaciones de diseño. En nuestra investigación, los eslabones son tratados como cuerpos flexibles tridimensionales y el proceso de optimización se realiza de manera simultánea en todos los cuerpos, mientras las condiciones de balanceo dinámico se cumplen.

Considerar al robot como un sistema flexible multicuerpo permite obtener un diseño optimizado con base en índices de desempeño estructural, tomando en cuenta la condiciones de balanceo dinámico. Estos índices de desempeño estructural pueden ser la frecuencias naturales, las deformaciones estáticas y la rigidez del robot, entre otras más. La optimización topológica de sistemas multicuerpo, en algunos casos, nos lleva a enfrentar análisis de gran escala los cuales requieren de un enorme esfuerzo computacional. Por esta razón, se desarrolló una plataforma computacional basada en cómputo paralelo, específicamente diseñada para la optimización de sistemas flexibles multicuerpo considerando restricciones de balanceo dinámico.

La idoneidad de la metodología propuesta es confirmada mediante la el diseño óptimo del mecanismo de cuatro barras sin reacción y el balanceo dinámico parcial del mecanismo de cinco v barras. La importancia del mecanismo de cuatro barras reside en la posibilidad de utilizar este mecanismo como una cadena cinemática especial para construir robots sin reacción. El mecanismo de cinco barras es muy importante debido a sus aplicaciones industriales, el cual típicamente es utilizado para operaciones de pick-and-place. A pesar de que el mecanismo de cinco barras se balancea parcialmente, la elimininación de las fuerzas resultantes desbalanceadas y la reducción de los momentos resultantes desbalanceados tienen un efecto positivo en el mecanismo, reduciendo las vibraciones e incrementando la precisión de posicionamiento.

Se realizaron validaciones numéricas de los mecanismos optimizados utilizando software comercial. El desempeño del balanceo dinámico de los mecanismo optimizados fue numéricamente validado utilizando ADAMS. Por otro lado, ANSYS fue empleado para evaluar las propiedades estructurales que fueron optimizadas.

Finalmente, es importante mencionar que se construyó un prototipo del mecanismo de cuatro barras dinámicamente balanceado y optimizado, ésto con el objetivo de probar la aplicación práctica de la metodología propuesta para el balanceo dinámico. El prototipo fue utilizado para realizar pruebas experimentales dónde se evaluó el desempeño de este mecanismo respecto al balanceo dinámico.

Palabras clave: balanceo dinámico, robots dinámicamente balanceados, optimización topológica multicuerpo, mecanismo de cuatro barras sin reacción, balanceo dinámico parcial del mecanismo de cinco barras.

Dedicated to the memory of my beloved mother María de Lourdes Hernández-Zaragoza, woman of infinite kindness and unbreakable strength.

vii

Introduction 1
The phases in design [START_REF] Budynas | Mechanical Engineering Design[END_REF] [START_REF] Herder | A counter-rotary counterweight (CRCW) for light-weight dynamic balancing[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.5 Full dynamic balancing of a 3-DOF planar parallel robot using counterweights and inertial flywheel [START_REF] Arakelian | Design of planar 3-DOF 3-RRR reactionless parallel manipulators[END_REF]. . . . . . . . . . . . . . . . . . . . . 10

1.6 Reactionless planar 3-DOF parallel robot [START_REF] Gosselin | Synthesis and design of reactionless three-degree-of-freedom parallel mechanisms[END_REF]. . . . . . . . . . . 11 1.7 Dynamic optimization of the planar four-bar linkage [START_REF] Feng | A new optimization method for dynamic design of planar linkage with clearances at joints-optimizing the mass distribution of links to reduce the change of joint forces[END_REF]. . . . . . . 12 1.8 Shape optimization of planar four-bar and five-bar linkages for dynamic balancing (Chaudhary and Chaudhary, 2015a). Initial linkage is drawn by a dashed line and the optimized linkage by solid line. . . . . . . . . . . . . . . . . . . . . . . . . . . 13

1.9 Typical process of structural topology optimization. . . . . . . . . . . . . . . . . . 15 1.10 Optimized MBB beam for minimum compliance constrained by the 40% of material. 17 1.11 Topology optimization of robotic manipulator for tracking error reduction [START_REF] Seifried | Optimal design of lightweight machines using flexible multibody system dynamics[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 1.12 Topology optimization of multibody system for compliance minimization [START_REF] Ghandriz | Structural topology optimization of multibody systems[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.13 Topology optimization of the upper frame of a painting robot [START_REF] Kim | Topology optimization of industrial robots for system-level stiffness maximization by using part-level metamodels[END_REF]. . . 20 1.14 Topology optimization of five-bar mechanism considering technology-oriented performance indices [START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF]. . . . . . . . . . . . . . . . . . . . 21 1.15 Structural topology optimization process for flexible multibody systems. . . . . . . 21 1.16 Topology optimization of four-bar linkage under balancing constraints (Briot and Goldsztejn, 2018a) [START_REF] Wu | Synthesis of reactionless spatial 3-DoF and 6-DoF mechanisms without separate counter-rotations[END_REF]. The reactionless spatial 3-DOF mechanism is used as a building leg in order to design the 6-DOF parallel robot. . . . . . . . . . . . . . . .

Symbols and abbreviations

Throughout the manuscript, vectors are represented by bold lowercase symbols and matrices by bold uppercase symbols. Scalar variables and names of points are in italic styles.

Symbols B i

Body i.

N i

Number of elements in the body mesh.

E 0 Material Young's modulus.

E min Young's modulus for very weak material.

E i j

Young's Modulus of the element j corresponding to body i.

ρ i j Density of the element i j.

p Penalization parameter.

q Partial penalization parameter.

ρ ρ ρ

Vector of design variables (densities).

K i j

Stiffness matrix of the element i j.

m i j Mass of the element i j.

m i Mass of body B i .

M i j Mass matrix of the element i j.

u i j Displacement vector of the element i j.

U e i j

Potential elastic energy of the element i j.

T e i j Kinetic energy of the element i j.

f Force vector.

xviii I Identity matrix.

Φ Φ Φ si Matrix of static modes of body B i .

Φ Φ Φ di

Matrix of dynamic modes of body B i .

R

Rotation matrix.

Q

Block-diagonal matrix staking on its diagonal rotation matrices.

I xxi , I yyi , I zzi The axial moments of inertia around x i , y i and z i axes, respectively, for body B i . These axes define a local frame attached to body B i with origin at O i .

I xyi , I yzi , I zxi The inertia products for body B i , expressed in its local frame defined by the x i , y i and z i axes and the origin O i .

I xi , I yi , I zi

The first moments of inertia (also known as static moments) around x i , y i and z i axes, respectively, for body B i . These axes define a local frame attached to body B i with origin at O i .

L

The Lagrangian of a system.

r P/Q Position vector of point P with respect to point Q.

v P/Q Velocity of point P with respect to point Q.

( * ) Time-derivative.

p Linear momentum.

h A Angular momentum with respect to point A.

h A x , h A y , h A z Components of the angular momentum h A around x, y and z axes, respectively.

f sh Shaking force vector.

m sh Shaking moment vector.

θ i Orientation of body B i with respect to X 0 -axis (for the plane case).

ℓ i Length of the body B i .

S i

Center of mass of the body B i .

r i Norm of the position vector of S i . The position vector of S i has been defined in a local frame rigidly attached to the body B i .

ψ i Orientation of the position vector of S i with respect to a local frame rigidly attached to the body B i .

n Unit vector.

f ( * ) Objective function.

h( * ) Equality constraint.

g( * )

Inequality constraint.

xix

Introduction Context of the thesis

Robots have become standard elements in the automated industrial processes. The robot mechanical structure (mechanical system) is constituted by several bodies interconnected between them with the aim to transmit motion and loads. The human nature pursues to improve the performance of any type of systems, as much as possible or as far as the frontiers of knowledge allow. Hence, the search of new improvements on the robots performance have been a very active research area in the last decades. These areas of study in the field of robotics are vast, but among the diverse topics which concern to robotics, the unbalanced inertial loads generated at the base of the robot during its operation, and the way to reduce or eliminate them, is the problem that concerns this thesis.

The mechanical design of the robot structure and its optimization is part of our interest. As can be seen in Fig. 1, the design process starts with the identification of the need, and culminated with the satisfactory solution, obtained by an exhaustive iterative process. In the definition of problem, the design specifications are established. Then, in the synthesis, the invention of the concept or concept design is how the need will be satisfied [START_REF] Budynas | Mechanical Engineering Design[END_REF]. The synthesis stage involves to fulfill requirements or performance criteria in terms of the geometric performance, kinematics, dynamics and structural behavior. After that, the proposed design is assessed in order to determine if its performance is satisfactory; subsequently the potential solutions are optimized.

The optimization stage can be a complex labor which implies, in some cases, to fulfill antagonistic requirements. Indeed, performing design optimization considering all the performance criteria simultaneously is a very complex task. Therefore, it is usual to conduct design optimization focusing in a subset ot the performance criteria. Thereby, this thesis is focused on the dynamic balancing of robotic mechanisms taking into account their structural optimization. In our proposal, mechanism links are optimized based on the structural response of the flexible multibody system in addition with the dynamic balancing properties. At its early days, the structural design was driven by trial and error, but it was time consuming, [START_REF] Budynas | Mechanical Engineering Design[END_REF].

expensive and prone to errors. Moreover, the design success depended on the experience of the designer, without the possibility to systematize the process. Nowadays, virtual prototyping enhance the evaluation stage, easing the hard labor of physical realization until the proposed design is well assessed. Once the successful design is obtained, the physical prototype can be built for a further evaluation. The technological development made possible the access to computer programs which allow to develop virtual prototypes where the definition of problem, synthesis, analysis and optimization can be conducted in the same software environment. Although optimization has been a complex task, currently structural optimization is commonly used in high-tech aerospace and automotive industry [START_REF] Cavazzuti | High performance automotive chassis design: a topology optimization based approach[END_REF][START_REF] Zhu | Topology optimization in aircraft and aerospace structures design[END_REF] to design functional systems obtained as result of the rigorous optimization process. Nevertheless, the general purpose software cannot be useful in all cases, and it happens when a novel methodology is under study. In that case, it is necessary to develop its own computational framework with the aim to evaluate the proposed methodology. Then, in order to verify the correctness of the obtained results, these should be validated by the specialized software if possible, this is usually known as numerical validation.

The dynamic balancing of robots is a fundamental topic in high-speed robots and also for freefloating space robots. Hence, the dynamic balancing is desirable in order to reduce vibrations, noise and wear. Moreover, positive effects of dynamic balancing are the increment of manipulator accuracy and cycle time reduction. When a mechanism is in motion, it produces inertia forces (shaking forces) and inertial moments (shaking moments), which generates dynamic reactions at its base due to the mass unbalance of the moving elements. When these dynamic reactions are canceled by a suitable dynamic balancing technique, the mechanism does not exert any unbalanced reaction at its base, and it is called a reactionless mechanism.

The design of reactionless mechanisms has focused on traditional directions, such as the use of counterweights and counter-rotations, as well as the addition of an auxiliary mechanism with the aim to produce dynamic balancing. Other researchers have proposed innovative alternatives, which include dynamic balancing via optimal motions of robot moving links, and the design of reactionless robots based on the dynamically balanced four-bar linkage (without counter-rotations). The major drawbacks of the classical balancing techniques are: (a) the total mass in the system is increased, (b) when auxiliary structures are used, the resultant mechanism is formed by a complex arrangements of bodies, and (c) all dynamics balancing design procedures lead to the design of links which do not take into account the elastic behavior of the links or the mechanism. Indeed, balancing a mechanism under elastic behavior constraints is a very complicated task.

In order to take into account the elastic behavior of the links, in this thesis we propose to use structural Topology Optimization (TO) as a tool for robot dynamic balancing. The proposed approach makes it possible to perform an optimal link design under balancing conditions. Consequently, we addressed the design of reactionless mechanisms as flexible mechanical systems, without attaching additional masses or using external mechanical elements. The methodology allows to perform an optimal design since the initial stage based on performance indices, which can be static deformations or structural vibration response, among others.

Scope of the doctoral thesis

The goal of this thesis is the study of topology optimization as a tool for the dynamic balancing of robotic mechanisms.

The specific objectives are the following:

1. Development of a computational platform capable of performing multibody topology optimization for dynamic balancing.

2. Carry out dynamic balancing of four-bar linkage using topology optimization.

Experimental validation of the dynamically balanced four-bar linkage.

4. To conduct the partial dynamic balancing of five-bar linkage using topology optimization.

Numerical validation of the optimized linkages.

In the following section a brief summary of thesis constributions is addressed.

Contributions of the thesis

The main contributions of this research work are described below:

• Development of parallel computational platform for topology optimization of multibody systems. Based on the formulation for topology optimization of multibody systems, a parallel computational platform was developed. It is specifically designed for the optimal design of flexible multibody systems, and it is used to perform dynamic balancing of mechanisms.

• Optimal design of mechanisms using topology optimization as a tool for dynamic balancing.

In order to evaluate the feasibility of the proposed methodology, two linkages were optimized taking into account their balancing conditions. The first case study is the four-bar linkage, which has one degree of freedom (DOF), but once it is dynamically balanced, it can be used as a module to build dynamically balanced planar and spatial parallel robots. In the second case it is the planar five-bar mechanism, which is very important due to its industrial applications: it is typically used in pick-and-place operations. Even if the mechanism is partially balanced, the shaking force elimination and shaking moment reduction can impact positively in reducing vibrations and increasing the position accuracy.

• Numerical validation of the optimized mechanisms by commercial software. In order to prove the accuracy of the proposed methodology, numerical validations of the optimized mechanisms were carried out. The optimization of the four-bar and five-bar mechanisms was focused on different performance indices alongside the balancing conditions. These optimized properties and the dynamic balancing were validated numerically using commercial software. The structural properties were validated with ANSYS and the dynamic balancing performance with ADAMS. In all cases the results shows a low error between the optimized values and the values obtained with commercial package.

• Experimental validation of the optimized four-bar linkage designed by the proposed approach. In the case of the four-bar linkage, a prototype was created in order to accomplish an experimental validation of its dynamic balancing. For this experiment, the base of the optimized linkage is suspended from a fixed structure with four cables, thus the linkage can moves freely in the case of an unbalanced reactive force is generated at its base. A direct current motor was used to provide oscillating movement. The unbalanced effects are embodied as translations and rotations of the suspended platform, these displacements are evaluated using video processing.

The results of the experiment shows a satisfactory dynamic balancing behavior.

Thesis structure

The manuscript is structure as follows. The thesis begins with a review of the classical balancing methods for shaking force and moment balancing, where the advantage and drawbacks of these methods are enlisted. Moreover topology optimization and its applications to flexible multibody systems and dynamic balancing were introduced in Chapter 1. Chapter 2 presents the complete mathematical formulation for modeling flexible multibody systems using topology optimization.

In addition, it is described the computational platform where this methodology is implemented.

Besides, the proposed approach for the solution of the dynamic balancing problem is then applied in the design of the four-bar and five-bar mechanisms. Chapter 3 presents an optimized design of a reactionless four-bar linkage, where the linkage stiffness is maximized constrained by the dynamic balancing conditions. The optimized properties of the mechanism are validated numerically and experimentally. The optimization of five-bar mechanism for partial dynamic balancing is studied in Chapter 4. In this case, the fundamental frequency of the mechanism is optimized while the shaking force balancing conditions are fulfilled. Furthermore the mechanism compliance is optimized, and shaking moment reduction is obtained by constraining the mechanism inertial properties. The numerical validation of these optimized properties is also reported. Finally, in Chapter 5 the conclusions of this research work are presented, highlighting the most important achievements.

Additionally, the perspectives of future research work are pointed out.

With the aim of enriching the work presented in the chapters, two appendices are included.

Appendix A is devoted to the dynamic balancing analysis of the five-bar linkage, and Appendix B consists in some comments on the manufacturing issues that can be faced during prototyping of dynamically balanced systems.

Chapter 1

State of the art

El trabajo aleja el vicio, el fastidio y la miseria.

"El amigo Pancho", Radio announcer.

D

YNAMIC balancing of robotic mechanism is of utmost importance for a comprehensive mechanical design. Mechanism balancing has been studied for a long time, and despite this there are still new contributions to the field. It happens largely, due to new challenges arising from the extensive use of robots for high-speed tasks. These new contributions for dynamic balancing are ranging from control approaches to optimum design and structural topology optimization. In such way that this topic continues being a very active field of research. In this chapter the state-of-the-art of dynamic balancing techniques for robotic mechanisms are described and analyzed. In addition, topology optimization is introduced, as well as its applications to the optimization of multibody systems and dynamic balancing of mechanisms.

Introduction

During high-speed tasks, large accelerations of mechanisms and robots lead to the generation of high inertial forces and moments, which are transmitted to the ground, causing the so-called shaking forces and shaking moments. These fluctuating loads are a significant source of vibration excitation, and lead to noise, fatigue and wear [START_REF] Lowen | Survey of investigations into the balancing of linkages[END_REF]). However, they can be canceled or reduced by suitable dynamic balancing techniques. Thus, if a mechanism does not exert any unbalanced reaction at its base it is called a reactionless or a dynamically balanced mechanism. Dynamic balancing is desirable to reduce vibrations, noise and wear [START_REF] Lowen | Survey of investigations into the balancing of linkages[END_REF]. Moreover, positive effects of dynamic balancing are the increment of manipulator accuracy and cycle time reduction [START_REF] Raaijmakers | Besi zoekt snelheidslimiet pakken en plaatsen op (Besi attacks the speedlimit for pick and place motion)[END_REF].

Classical dynamic balancing techniques

Shaking forces and moments can be reduced or eliminated using specific mechanism design, and proper inertia parameters. Moreover, there are different self-balanced systems generated by duplicating the original mechanism and other types of mechanical arrangements, examples of these auxiliary structures can be found in [START_REF] Arakelian | Balancing of Linkages and Robot Manipulators -Advanced Methods with Illustrative Examples[END_REF][START_REF] Bagci | Complete shaking force and shaking moment balancing of link mechanisms using balancing idler loops[END_REF][START_REF] Hilpert | Weight balancing of precision mechanical instruments[END_REF], however those mechanism are not considered in this work. There are several literature reviews which present the achievements and trends on classical dynamic balancing techniques [START_REF] Arakelian | Inertia forces and moments balancing in robot manipulators: a review[END_REF][START_REF] Arakelian | Shaking force and shaking moment balancing fo mechanisms: a historical review with new examples[END_REF][START_REF] Lowen | Survey of investigations into the balancing of linkages[END_REF][START_REF] Lowen | Balancing of linkages-an update[END_REF][START_REF] Wei | A review of dynamic balancing for robotic mechanisms[END_REF].

In addition, a detailed review of the classical balancing techniques can be found in the textbook [START_REF] Arakelian | Balancing of Linkages and Robot Manipulators -Advanced Methods with Illustrative Examples[END_REF]. In general, the dynamic balancing techniques are divided in two main groups: shaking force balancing (SFB) and shaking moment balancing (SMB), either to eliminate or to reduce the unbalanced reaction loads. The following sections are dedicated to describe these methods.

Shaking force balancing techniques

The methods for shaking force balancing relies on keeping the mechanism center of mass (COM) stationary, and they can be classified in the following subgroups [START_REF] Arakelian | Inertia forces and moments balancing in robot manipulators: a review[END_REF]:

(a) Shaking force balancing by adding counterweights [START_REF] Berkof | A new method for completely force balancing simple linkages[END_REF][START_REF] Fisher | Über die reduzierten systeme und die hauptpunkte der glieder eines gelenkmechanismus[END_REF][START_REF] Crossley | Dynamics in Machines[END_REF]. These methods use counterweights, attached to the movable links of the mechanism, in order to keep the linkage COM stationary. The drawback is the addition of mass, which decreases the mechanism's natural frequencies and enlarges the efforts in the (a) SFB by auxiliary structures [START_REF] Agrawal | Reactionless space and ground robots -novel designs and concept studies[END_REF].

(b) SFB by auxiliary structures [START_REF] Fattah | On the design of reactionless 3-DOF planar parallel mechanisms[END_REF]. joints. Figure 1.1 shows a schematic diagram of a four-bar linkage which has been balanced by adding counterweights. (b) Shaking force balancing by adding auxiliary structures [START_REF] Agrawal | Reactionless space and ground robots -novel designs and concept studies[END_REF][START_REF] Arakelian | Shaking moment cancellation of self-balanced slider-crank mechanical systems by means of optimum mass redistribution[END_REF][START_REF] Kamenskii | On the question of the balancing of plane linkages[END_REF][START_REF] Van Der Wijk | Dynamic balancing of Claver's Delta robot[END_REF]. These methods are based on the use of an axially symmetric duplicate mechanism in order to make the combined center of mass stationary. Theses approaches imply to have a more complex linkage arrangement, with an increase in the total mass of the system. Figure 1.2 present two examples of this approach. In Fig. 1.2a an auxiliary parallelogram was added in order to suitably locate the center of mass. In Fig. 1.2b a 3-DOF planar parallel mechanism is force balanced by using an auxiliary structure.

(c) Shaking force balancing by adjustment of kinematic parameters (AKP) [START_REF] Ouyang | Integrated design of robotic mechanisms for force balancing and trajectory tracking[END_REF][START_REF] Ouyang | Force balancing of robotic mechanisms based on adjustment of kinematic parameters[END_REF][START_REF] Yu | Partial shaking moment balancing of spherical parallel robots by a combined counterweight and adjusting kinematic parameters approach[END_REF]. For this method, the mass of the links is determined a priori, and the length and mass center of the links are computed in order to perform (Briot et al., 2012, 2010[START_REF] Geng | Shaking force balancing of the delta robot[END_REF]. This alternative does not require mass redistribution of the moving links.

The technique consists in planning the displacements of the total COM of the moving links, in order to minimize the center of mass acceleration. However, this implies to follow some specific trajectories which may not be the desired ones by the user.

Besides, strictly speaking, the force balancing by using elastic components [START_REF] Alici | Optimum force balancing of a planar parallel manipulator[END_REF] cannot be considered as a technique for shaking force balancing. The elastic elements connected to the driving links only affect the gravitational forces and the input torques, but they do not have any influence on the shaking force.

Shaking moment balancing techniques

On the other hand, the methods for shaking moment balancing can be classified in the next subgroups:

(e) Shaking moment balancing by counter-rotation [START_REF] Arakelian | Design of planar 3-DOF 3-RRR reactionless parallel manipulators[END_REF][START_REF] Berestov | Comparative analysis of the reactions in the kinematic pairs of the fourbar linkages for the different balancing methods[END_REF][START_REF] Kochev | General theory of complete shaking moment balancing of planar linkages: a critical review[END_REF][START_REF] Van Der Wijk | Comparative analysis for low-mass and low-inertia dynamic balancing of mechanisms[END_REF]. This technique enables the balancing of planar linkages by adding a system which generates a counter-rotation with proportional angular acceleration. Thus the system will supply an equal and opposite balancing moment. This technique has the disadvantage of increasing the linkage total mass, which leads to higher [START_REF] Herder | A counter-rotary counterweight (CRCW) for light-weight dynamic balancing[END_REF].

Figure 1.5: Full dynamic balancing of a 3-DOF planar parallel robot using counterweights and inertial flywheel [START_REF] Arakelian | Design of planar 3-DOF 3-RRR reactionless parallel manipulators[END_REF].

input torque requirements. The most common way to generate counter-rotations is by using gears, as can be seen in Fig. 1.4.

(f) Shaking moment balancing by adding an inertia flywheel rotating with a prescribed angular velocity [START_REF] Arakelian | Design of planar 3-DOF 3-RRR reactionless parallel manipulators[END_REF][START_REF] Van Der Wijk | Dynamic balancing of Claver's Delta robot[END_REF][START_REF] Van Der Wijk | Active dynamic balancing unit for controlled shaking force and shaking moment balancing[END_REF]. This method is based on the fact that, after performing shaking force balancing, the resultant shaking moment for a given linkage configuration has the same value for any location on the ground [START_REF] Berkof | Theory of shaking moment optimization of force-balanced four-bar linkage[END_REF]. Thus, any planar linkage can be balanced adding an inertia flywheel rotating with a prescribed angular velocity. Nevertheless, the main issue with this technique is how to generate precisely the required input motion. An illustrative scheme of this technique appears in Fig. 1.5, where the a 3-DOF planar parallel robot is shaking moment balanced using an inertial flywheel. (g) Shaking moment balancing by generating optimal trajectories of moving links (He and Lu, Planar 2-dof mechanism [Gosselin et al., (a) Reactionless planar 2-DOF mechanism based on the dynamically balanced four-bar linkage.

2004] a b Planar 3-dof mechanism [Gosselin et al., 2004] (b) Planar 3-DOF mechanism built with four-bar dynamic balanced modules.

Figure 1.6: Reactionless planar 3-DOF parallel robot [START_REF] Gosselin | Synthesis and design of reactionless three-degree-of-freedom parallel mechanisms[END_REF].

2006, [START_REF] Papadopoulos | Design and motion planning for a zero-reaction manipulator[END_REF]. Once the shaking force is solved, then the shaking moment balancing can be achieved by planning the linkage optimal motions. Again, this implies to follow some specific trajectories which may be different from the desired ones.

In addition to the described solutions for shaking force and moment balancing, there are other approaches which use a dynamically balanced four-bar linkage as a leg to building reactionless robots. The dynamic balancing conditions to design a reactionless four-bar linkage without counterrotations were described in (Ricard and [START_REF] Gosselin | On the development of reactionless parallel manipulators[END_REF], where the authors pointed out the viability of building reactionless robots using this mechanism. Subsequently, in [START_REF] Gosselin | Synthesis and design of reactionless three-degree-of-freedom parallel mechanisms[END_REF]Gosselin, 2004) the design of planar and spatial robots was reported, proving the viability of the proposed strategy. In Fig. 1.6 is shown the designed 3-DOF reactionless planar parallel robot based on the reactionless four-bar linkage. This approach is very attractive because of the simple design of the reactionless four-bar mechanism. In [START_REF] Foucault | Synthesis, design and prototyping of a planar three degeeof-freedom reactionless parallel mechanism[END_REF] the authors presented a similar approach addressing the dynamic balancing of a planar 3-DOF parallel mechanism, but the building leg is based on the five-bar parallelogram linkage, which require counter-rotations in order to be fully dynamically balanced, thus resulting in a more complex linkage arrangement.

Alternatively, there are some works which explore the optimization of the external links shape in order to minimize the inertial forces and moments. The resultant shape is generated by means of optimal mass distribution which minimizes the shaking force and moment. Some of these works are described in the following section. [START_REF] Feng | A new optimization method for dynamic design of planar linkage with clearances at joints-optimizing the mass distribution of links to reduce the change of joint forces[END_REF].

Shape optimization for inertial loads minimization

A shape optimization procedure for dynamic design of planar linkages is reported in [START_REF] Feng | A new optimization method for dynamic design of planar linkage with clearances at joints-optimizing the mass distribution of links to reduce the change of joint forces[END_REF], where the authors performed dynamic optimization of the planar four-bar linkage controlling the changes in the joint forces while the joint clearances are considered. In addition a technique called "small element superposing method" was proposed in the same paper, in order to generate the external link shape, based on the optimized link parameters. This technique does not include the elasticity of the links in its formulation. The resultant optimized four-bar linkage is presented in Fig. 1.7.

Another optimization framework approach to optimize linkages was proposed in (Chaudhary and Chaudhary, 2015a,b), where the authors used shape optimization based on the rigid body dynamics for optimizing the external link shape (only) of some linkages. Figure 1.8 presents the results of this methodology, where the initial linkage is drawn by a dashed line, and the optimized linkage is represented by the solid line. Since the approach is based on rigid body formulation, the deformations of the link cannot be considered during the optimization process. Moreover, there is no possibility of modifying the internal link shape (including some voids, for instance).

After having exposed all the aforementioned works for dynamic balancing, it becomes clear that the linkage elastic performances are never optimized, despite them being of the utmost importance in robot design. Therefore, it is concluded that all the methods presented above have some major drawbacks: (1) the total mass in the system is increased, (2) when auxiliary structures / counterrotations are used, the resultant mechanism is formed by a complex arrangements of bodies, and (3) typical dynamic balancing design procedures lead to the design of links which do not take into account the elastic behavior of the links or the mechanism. Although the linkage dynamic balancing under elastic behavior constraints is very complicated task, it is an important matter in high-speed robot design.

There are few research works which studied the dynamic balancing conditions when flexible members are considered. In [START_REF] Kalas | Shaking force balance in parallelmanipulators with flexible links[END_REF] the author considered the inclusion of flexible links in the five-bar linkage, thus the flexible links are modeled as pseudorigid bodies [START_REF] Edwards | A pseudo-rigid-body model for initiallycurved pinned-pinned segments used in compliant mechanisms[END_REF][START_REF] Howell | Compliant mechanisms[END_REF] and then the author presents some strategies in order to improve the shaking force balancing under these conditions. More recently Meijaard and van der Wijk (2022) have presented the dynamic balancing of flexible mechanisms considering two approaches, the symilarity and modal balancing. However, considering the direct effects on the dynamic balancing conditions due to the flexibility of the links are out of the scope of this research work.

On that account, in order to consider the elastic behavior of the links, the structural topology optimization (TO) [START_REF] Bendsøe | Topology Optimization -Theory, Methods and Applications[END_REF]) is introduced as a potential tool that can be used to overcome the aforementioned problems. The next section is devoted to present a topology optimization overview.

Topology optimization

Since [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF] introduced a method for material distribution, topology optimization has been a very active research field, ranging from structural topology optimization, microstructure design, acoustics, and fluids, among others. Topology optimization is a powerful mathematical method which aims to redistribute the material into an initial domain taking into account design specifications [START_REF] Bendsøe | Topology Optimization -Theory, Methods and Applications[END_REF]. The optimization process is based on the structural response, typically computed by the Finite Element Method (FEM). Topology optimization considers the whole domain of the structure, instead of the shape (external) or size of the structural members.

Topology optimization can be formulated as a problem with continuous or discrete design variables. Topology optimization formulated with continuous design variables can be solved by different approaches such as density-based (Bendsøe, 1989, Zhou and[START_REF] Zhou | The COC algorithm, Part II: Topological, geometrical and generalized shape optimization[END_REF], topological derivatives [START_REF] Novotny | Topological Derivatives in Shape Optimization[END_REF]Sokołowski, 2013, Sokolowski and[START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF], level set [START_REF] Allaire | A level-set method for shape optimization[END_REF][START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF] or phase field method [START_REF] Bourdin | Design-dependent loads in topology optimization[END_REF]. Despite of the variety of topology optimization approaches, there are small differences between them [START_REF] Sigmund | Topology optimization approaches[END_REF], but the most well known and mature is the density-based approach.

The density-based approach can be modeled using a material interpolation scheme called Simplified Isotropic Material with Penalization (SIMP), which was initially proposed in [START_REF] Bendsøe | Optimal shape design as a material distribution problem[END_REF], and later on a modified SIMP scheme was defined in [START_REF] Sigmund | Morphology-based black and white filteres for topology optimization[END_REF]. However, there are other material interpolation methods such as RAMP [START_REF] Stolpe | An alternative interpolation scheme for minimum compliance optimization[END_REF], but SIMP method is one of the approaches used in industrial software [START_REF] Rozvany | A critical review of established methods of structural topology optimization[END_REF].

The SIMP scheme is a method of topology optimization based on the Finite Element (FE) formulation, thereby, the design domain is divided (discretized) into many small pieces (finite elements), then each element is associated with an artificial density ρ e , which is the design variable bounded to take values between 0 and 1. Using the density of each element its Young's modulus E e can be defined, and the modified SIMP scheme [START_REF] Sigmund | Morphology-based black and white filteres for topology optimization[END_REF] is used to penalize intermediate density values, thus it is expressed as follows:

E e = E min + ρ p e (E 0 -E min ) , with ρ e ∈ [0, 1]. (1.1)
where p is the penalization power (usually p = 3, for elasticity problems), E min is the stiffness at ρ e = 0 (E min is different from 0 in order to avoid singularity of the stiffness matrix), and E 0 is the Young's modulus of the material. Hence, the stiffness matrix of each element is given by:

K e (ρ e ) = E e (ρ e )K 0 e (1.2)
where K 0 e is a constant stiffness matrix for an element with Young's modulus equal to one. Besides, the global stiffness matrix can be obtained by a standard assembly procedure. Consequently, the optimization problem modeled by the SIMP scheme has the standard form of a nonlinear program: where f (ρ ρ ρ) is the objective function, while g (ρ ρ ρ) and h (ρ ρ ρ) are the inequality and equality constraints, respectively. The vector ρ ρ ρ contains decision variables, and each component of this vector represents the density of a given element. The objective function typically models a performance criteria, and the constraints describe restrictions on the structure.

min ρ ρ ρ∈[0, 1] : f (ρ ρ ρ) subject to : g (ρ ρ ρ) ≤ 0 (1.3) : h (ρ ρ ρ) = 0 Initialize FEM FEM Filtering Optimization solver Converged? YES K(ρ)u = F F K 0 i K(ρ) = N i=1 E i (ρ i ) K 0 i SIMP E i (ρ i ) min ρ∈[0, 1] n : f (ρ) s.t.: g (ρ) ≤ 0, h (ρ) = 0 MMA, OC
The formulation given in (1.3) is the so-called nested approach, where the equilibrium equations coming from the FE model must be satisfied for each optimization step. Therefore a suitable linear solver for the solution of the FE model is necessary. The schematic process of topology optimization for a single body is described in Fig. 1.9. The optimizer, gradients, and filtering are described in the following sections.

Other popular TO formulations are developed for discrete design variables, they are the evolutionary approaches usually known as Evolutionary Structural Optimization (ESO) [START_REF] Xie | A simple evolutionary procedure for structural optimization[END_REF] and Bidirectional Evolutionary Structural Optimization (BESO) [START_REF] Querin | Evolutionary structural optimisation (ESO) using a bidirectional algorithm[END_REF][START_REF] Yang | Bidirectional evolutionary method for stiffness optimization[END_REF]. In their early days these techniques were criticized for their heuristic approach, and the difficulty for assess convergence qualities [START_REF] Rozvany | A critical review of established methods of structural topology optimization[END_REF]. Afterwards, they were considered as the discrete version of the standard SIMP scheme [START_REF] Sigmund | Topology optimization approaches[END_REF]. Nowadays evolutionary approaches are still being improved [START_REF] Huang | Evolutionary Topology Optimization of Continuum Structures: Methods and Applications[END_REF] and are part of several commercial software due to their simplicity to be utilized with finite element software [START_REF] Deaton | A survey of structural and multidisciplinary continuum topology optimization: post 2000[END_REF].

Optimizer

The optimization problem described in Eq. (1.3) represents a smooth and differentiable problem which can be efficiently solved by well-establish gradient-based optimization methods [START_REF] Sigmund | Topology optimization approaches[END_REF] such as Optimality Criteria (OC) and the Method of Moving Asymptotes (MMA) [START_REF] Svanberg | The method of moving asymptotes -a new method for structural optimization[END_REF]. These gradient-based methods require the analytical description of the first order derivatives of the objective function and of the constraints with respect to the design variables. The derivatives are commonly called gradients and also sensitivities, thus it is usual to designate the computation of the derivatives as a sensitivity analysis.

The gradients provide useful information to the gradient-based optimizer, which is used for the algorithm to converge more efficiently. Hence, the computation of the gradients is a fundamental part of these types of optimizers. In some cases since the gradients cannot be computed analytically, some approximations to the gradient should be computed. Among these approximation techniques are the finite-difference methods, which are popular because of their simplicity, but these techniques are neither accurate nor computationally efficient (Martins and[START_REF] Martins | Review and unification of methods for computing derivatives of multidisciplinary computational models[END_REF][START_REF] Martins | Engineering Design Optimization[END_REF]. Indeed, when possible, analytical gradients are preferred instead of numerical approximations, because the analytical gradients are evaluated instead of computed in each optimization step.

Besides, in topology optimization literature some optimizer known as Non-Gradient Topology Optimization schemes (NGTO) appears as an alternative to the Gradient-based Topology Optimizations (GTO) methods. In [START_REF] Sigmund | On the usefulness of non-gradient approaches in topology optimization[END_REF] the strengths of GTO methods were described, compared with those NGTO alternative methods, such as Genetic Algorithms (GA), simulated annealing, etc. The author pointed out that those NGTO methods should be used in very special cases due to their higher computational time required to find a satisfactory solution for small-scale problems. In general, the GTO methods are much better option than NGTO methods, because the gradient provides important information on where the material or void is needed [START_REF] Sigmund | On the usefulness of non-gradient approaches in topology optimization[END_REF].

On the other side, Neural Networks (NN) or Machine Learning (ML) techniques are commonly used as accelerators in topology optimization by data-driven training [START_REF] Banga | 3D Topology optimization using convolutional neural networks[END_REF][START_REF] Sosnovik | Neural networks for topology optimization[END_REF][START_REF] Ulu | A data-driven investigation and estimation of optimal topologies under variable loading configurations[END_REF], obtaining interesting results in terms of computational time reduction, but these papers discuss simple benchmarks problems, mostly with coarse mesh. However, in latest studies conducted by [START_REF] Xue | Efficient, high-resolution topology optimization method based on convolutional neural networks[END_REF], more complex problems are analyzed and finer meshes are used, obtaining promising results. In addition, [START_REF] Chandrasekhar | TOuNN: Topology optimization using neural networks[END_REF] propose to execute directly (instead of using as a training or accelerator tool) topology optimization using NN. The results show that the proposed methodology generates optimized designs by twice more the computational cost of the typical optimization formulation given in [START_REF] Andreassen | Efficient topology optimization in MAT-LAB using 88 lines of code[END_REF].

Therefore, an efficient approach for topology optimization relies on a gradient-based optimization procedure alongside the analytic derivation of the gradients. For structural design the most natural choice are the gradient-based optimization methods, and the problem faced in this thesis is well suited for these methods. In addition, the problem we treat is characterized by a high number of variables. Hence, in order to solve the optimization problem it is necessary to use an optimizer, able to handle large number of variables, as well as multiple constraints. Thus, the method of moving asymptotes introduced in [START_REF] Svanberg | The method of moving asymptotes -a new method for structural optimization[END_REF] is the optimizer used for this research work.

It is able to handle multiple inequality and equality constraints, and is the most accepted optimizer in the structural optimization community because of its excellent convergence properties [START_REF] Bendsøe | Topology Optimization -Theory, Methods and Applications[END_REF].

Numerical problems in topology optimization

For density-based topology optimization formulation there are some issues that affect the computational results. The common numerical problems that appears in TO when no regularization scheme is applied are: a) checkerboard pattern, b) mesh dependency, and c) local minima [START_REF] Sigmund | Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards mesh-dependencies and local minima[END_REF]. In order to overcome checkerboard patterns and mesh dependency it is necessary to use restriction methods [START_REF] Sigmund | Morphology-based black and white filteres for topology optimization[END_REF], which operate as filters. The classical approach is the application of a filter to modify either the sensitivities or densities. Local minima are mainly generated by the optimizer used to solve the problem. Topology optimization problems have many local minima, even more for those problems with many intermediate density values. Therefore, by using a method to avoid mesh-dependency it is possible to prevent local minima [START_REF] Sigmund | On the design of compliant mechanisms using topology optimization[END_REF].

The checkerboard patterns appear as regions with alternating solid and void elements resembling a checkerboard. It is due to a wrong numerical solution of the FE model, where there appears artificially high stiffness in the elements [START_REF] Diaz | Checkerboard patterns in layout optimization[END_REF]. Checkerboard problem can be avoided partially by using higher order elements or it can be avoided by using restriction methods (Lazarov and[START_REF] Sigmund | On the usefulness of non-gradient approaches in topology optimization[END_REF][START_REF] Sigmund | Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards mesh-dependencies and local minima[END_REF]. Among these restriction methods are the filtering schemes that ensure mesh independence. Mesh dependence refers to the problem of obtaining qualitatively different solutions for different mesh discretizations [START_REF] Bendsøe | Topology Optimization -Theory, Methods and Applications[END_REF]. Thus, when a mesh is refined it is expected to obtain qualitatively the same optimal structure not a different structure, but when mesh-dependency appears, the mesh density controls the solution, which is a undesirable effect. Therefore, in order to prevent mesh-dependence and checkboard patterns some restriction methods should be used, where the most extended approach is the use of filters applied to the sensitivities [START_REF] Sigmund | Design of Material Structures Using Topology Optimization[END_REF][START_REF] Sigmund | On the design of compliant mechanisms using topology optimization[END_REF] or the densities (Bourdin, 2001, Bruns and[START_REF] Bruns | Topology optimization of non-linear elastic structures and compliant mechanisms[END_REF]. The sensitivity filter modifies the sensitivities of the objective function, while the density filter modifies the objective function and the constraints. Figure 1.10 represents the solution of a typical benchmark problem (MBB beam) for TO, where the compliance is minimized and it is constrained by the amount of material. The optimized result without use a filter is shown in Fig. 1.10a, where the checkboard patter is clearly present. When the sensitivity filter is used, a more regular design is obtained, as it can be appreciated in Fig. 1.10b. Some applications of topology optimization in the field of Multibody Systems (MBS) are described in the following section. TO has been used to generate optimized single bodies, nevertheless there are relatively few research works which study Multibody Topology Optimization (MTO).

Topology optimization of flexible multibody systems

The design of flexible multibody systems based on structural optimization has attracted the attention of some researchers. One of the first attempts to face the structural optimization of multibody systems was described in [START_REF] Brüls | Multibody Dynamics, chapter Opti-mization of multibody systems and their structural components[END_REF], where the authors considered the bodies of the system like truss structures (defined a priori), and minimized the compliance of the truss members constrained by a predefined amount of material. Nevertheless this approach is very restrictive because the general shapes of the bodies cannot be taken into account. The inclusion of structural topology optimization in multibody systems begins with the optimization of isolated parts which are components of the whole system [START_REF] Albers | Upper body of a new humanoid robot -the design of ARMAR III[END_REF], 2007[START_REF] Lohmeier | Leg design for a humanoid walking robot[END_REF][START_REF] Rozvany | A critical review of established methods of structural topology optimization[END_REF]. After that, some studies were conducted considering the whole multibody system for structural topology optimization (density-based approach). [START_REF] Seifried | Optimal design of lightweight machines using flexible multibody system dynamics[END_REF] performed shape optimization and topology optimization of a 3-DOF planar robotic manipulator, based on the dynamic simulation of the flexible multibody system for a set of displacements fields. They conducted TO for compliance minimization while tracking errors were reduced for only one arm of the mechanism. The robotic Figure 1.11: Topology optimization of robotic manipulator for tracking error reduction [START_REF] Seifried | Optimal design of lightweight machines using flexible multibody system dynamics[END_REF].

mechanism and the optimized arm can be seen in Fig. 1.11. On the other side, in [START_REF] Ghandriz | Structural topology optimization of multibody systems[END_REF] the authors presented an algorithm for topology optimization of multibody systems, based on the dynamic (time dependent) behavior of the entire multibody system, which is focused on planar systems with large rotational and translational motions. The initial mechanism and the optimized one are shown in Fig. 1.12. The approaches described in (Ghandriz et al., 2017, Seifried and[START_REF] Seifried | Optimal design of lightweight machines using flexible multibody system dynamics[END_REF] deal with the relevant loads on the flexible multibody system, trying to capture the global dynamic behavior. Nevertheless both works reported numerical difficulties for the sensitivity analysis. These issues were overcome assuming independence between the compliance (objective function) and the dynamic loads, which resembles the static problem.

As a matter of fact, it is still a challenge to develop an efficient methodology capable to capture the complete dynamic behavior of the multibody system while it is optimized using topology optimization. For that reason alternative methodologies have been proposed to solve this problem assuming specific conditions, but with satisfactory results [START_REF] Jang | Dynamic response topology optimization in the time domain using equivalent static loads[END_REF][START_REF] Kang | Optimization of flexible multibody dynamic systems using the equivalent static load method[END_REF][START_REF] Moghadasi | Topology optimization of members of flexible multibody systems under dominant inertia loading[END_REF]. For instance, in [START_REF] Kim | Topology optimization of industrial robots for system-level stiffness maximization by using part-level metamodels[END_REF] the authors performed topology optimization of industrial robots using part-level metamodels on a commercial software. It implies to solve the problem by isolating each body from the rest of the system, but the metamodel incorporates the corresponding information of the entire system (system-level). The researchers minimized the strain energy constraining the robot mass and also solved the problem by considering dynamic loads, but they obtained only a close-to-optimum due to part-level metamodels. Figure 1.13 presents the results of this optimization procedure.

Typically, structural topology optimization is employed for energy efficiency purposes. Hence, Figure 1.12: Topology optimization of multibody system for compliance minimization [START_REF] Ghandriz | Structural topology optimization of multibody systems[END_REF].

Figure 1.13: Topology optimization of the upper frame of a painting robot [START_REF] Kim | Topology optimization of industrial robots for system-level stiffness maximization by using part-level metamodels[END_REF].

almost all the research works deal with the objective to generate lightweight design by means of mass reduction while the stiffness is maximized. There are few research works which attempt to solve more general problems in the field of multibody topology optimization. Additionally, as it has been mentioned, due to the high computational cost and the complexity of this type of problems, it is imperative to find an adequate trade-off between the complexity in the problem formulation and the accuracy of the results.

As a consequence, a detailed methodology to design industrial robots using topology optimization was presented in [START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF]. The authors described the mathematical formulation of the optimization problem as a flexible multibody system. Besides, they proposed to use model reduction techniques in order to decrease the computational effort in the analysis. The authors presented the optimal design of a planar 2-DOF robot, and they analyzed two cases: in the first case the links are treated as planar bodies and for the second case they are considered as threedimensional bodies. In both cases the optimization was performed for all bodies simultaneously, and it was based on the full system response for a set of critical trajectories in order to enforce the global optimization. The results of the two cases are given in Fig. 1.14.

Typical problems in multibody topology optimization deal with minimum compliance and mass reduction in order to generate lightweight design and energy-efficient robots. Nevertheless, it is difficult to use the compliance value (solely) as a design parameter, and also there are other criteria which can affect the performance of the robots. Therefore, [START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF] proposed some technology-oriented performance indices which are easy to understand when the structure of the robot is designed, as it is the case of the static deformations, energy consumption and the natural frequencies. An schematic description of the proposed topology optimization for multibody systems is shown in Fig. 1.15, where unlike the TO for a single body, in this case it is proposed to model the entire multibody system and apply model reduction techniques in order to reduce the computational effort.

Dynamic balancing using topology optimization

The current literature shows that there are just a few research works which attempt to solve the dynamic balancing of flexible multibody systems. The optimal design based on TO for a reactionless four-bar linkage was reported in (Briot and Goldsztejn, 2018a), which is based on the methodology described in [START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF]. In (Briot and Goldsztejn, 2018a) the links are treated as a two-dimensional (2D) bodies, the first natural frequency of the four-bar linkage is maximized, and the mechanism compliance is constrained while ensuring the balancing conditions for rigid-body mechanisms. These dynamic balancing conditions were obtained from (Ricard and [START_REF] Gosselin | On the development of reactionless parallel manipulators[END_REF], whose advantage is that they not require the use of counter-rotations in order to achieve the full dynamic balancing. In addition, these conditions were rewritten in terms of the design variables in order to include them as constraints in the topology optimization formulation. These optimized links of the four-bar linkage are presented in Fig. 1.16. In contrast to typical dynamic balancing, the resultant optimized bodies can have a different shape (internal/external) depending on the load system, but the dynamic balancing conditions must be fulfilled.

The results obtained by Briot and Goldsztejn (2018a) reveal the potential of topology optimization as a valuable tool for linkage design under dynamic balancing constraints. Nevertheless, the optimized properties were not validated, and the links of four-bar linkage were considered as a planar bodies. This avoids the possibility of considering the link bending in all directions, and restricts the optimal design for planar loads systems.

The design approach reported in the present thesis introduces several advantages with respect to previously published researches, related to dynamic balancing of flexible multibody systems using topology optimization. With respect to (Briot and Goldsztejn, 2018a), at first, the links are treated as three-dimensional (3D) flexible bodies, being thus possible to optimize the links for a general load system. Secondly, in order to validate the reliability of our solution proposal, numerical validation of the optimized properties were realized using specialized software (ADAMS and ANSYS). In third place, the optimization is performed considering the experimental evaluation of the dynamic balancing, thus the optimization includes the inertia of axes and bearings. Thereby, a prototype was built in order to carry out an experimental evaluation of the linkage dynamic balancing behavior.

Additionally, in this thesis the topology optimization of the five-bar mechanism is conducted in order to accomplish its partial dynamic balancing. To the best of our knowledge, this is the first time that TO is used for this purpose. The five-bar mechanism is a robotic architecture of 2-DOF used typically for pick-and-place industrial operations, because this mechanism has the possibility to move its end-effector on a plane.

The formulation for multibody topology optimization used in this thesis is based on [START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF]. Nonetheless, even if the model reduction techniques were successfully applied, the computational cost for solving the problem was high, even more for the three dimensional case, where in some cases the use of a fine discretization becomes prohibitive. Thus, in this thesis a computational platform is specifically developed for the efficient computation of multibody topology optimization using parallel processing and the scientific programming language C++. The used formulation for multibody topology optimization allows to optimize the system by considering the response of the entire system. In addition, the optimization can be conducted for 22 a set of trajectories. Nevertheless, in order to reduce the computational complexity and due to the numerical instabilities generated during the trajectory optimization, for this research work the mechanisms are analyzed for one configuration. This assumption does not affect the main objective of the thesis, because the shaking force balancing is independent of the mechanism configuration.

Regarding to the shaking moment balancing, it is minimized using the inertial mechanism parameters, which also do not depend on the configuration. Indeed, the structural performance of the mechanism is optimized only for the selected configuration.

In addition to the dynamic balancing techniques, there are some important numerical and experimental investigations related to the effects of mass-balancing on the elastodynaic performance of the balanced mechanisms [START_REF] Martini | Modaland kineto-elastodynamic analyses of balanced four-bar linkages[END_REF][START_REF] Martini | Elastodynamic effects of mass-balancing: experimental investigation of a four-bar linkage[END_REF]. For instance [START_REF] Martini | Elastodynamic effects of mass-balancing: experimental investigation of a four-bar linkage[END_REF] conducted an experimental evaluation of the elastodynamic effects of a mass-balanced four-bar linkage, where the authors concluded that the use of counter-weights on the four-bar linkage generates a deteriorate dynamic operation at low/medium speed. Besides the authors suggest that these conclusions should considered in general for closed-loop mechanisms. In our studies we do not considered the dynamic effects of the mass-redistribution in the optimized linkages, however we optimize the stiffness of the whole linkage (while dynamic balancing conditions are fulfilled) assuming that this feature have a positive impact on the mechanism dynamic performance [START_REF] Thompson | A survey of finite element techniques for mechanism design[END_REF]). Nevertheless we envisage that further investigations should be conducted in this research direction.

Summary

This Chapter presented an overview of the classical balancing methods related to mechanisms and robots. These methods have the main disadvantage of being formulated for the analysis of rigid bodies, therefore they do not consider the elastic behavior of the mechanisms, which is an important consideration for comprehensive design, mainly in high-speed robot design.

Classical methods propose to include additional components in order to mitigate or eliminate the effects of the inertial loads, but the use of counter-weights and counter-rotations, and even more 1.6. SUMMARY auxiliary structures increase the total mass of the system and its complexity. Thereby, it is more appealing to conceive solutions which do not need additional components, as it is the case of the self-balanced four-bar linkage, which does not require counter-rotations. Therefore, for this thesis the inclusion of additional components is not considered in order to perform dynamic balancing.

There are some works which solve the partial dynamic balancing problem by means of shape optimization. These works are based on the optimization of the mass distribution, wherewith the external shape is defined. Nevertheless, these formulations do not take into account the elastic behavior of the links and the shape optimization makes impossible to modify the internal link shape. A more general formulation is the TO, where the material is redistributed within a design domain subjected to specific constrains.

A brief introduction in the field of topology optimization was presented. For the problem formulation the density-based approach and the SIMP scheme were selected, which implies to model the system by the FEM. Because of the nature of the problem faced in this thesis, the best approach to solve the optimization problem are the gradient-based methods. Additionally, due to the high number of variables and the use of several constraints, the most suitable method for this endeavor is the MMA. One important stage in topology optimization is the solution of the equilibrium equations belonging to the FE model. Due to the large-scale model generated in multibody analysis, the solution of these equations should be carried out using an efficient iterative solver.

In order to model the mechanisms as flexible multibody systems it is proposed to use multibody topology optimization. The proposed methodology makes possible to take into account the linkage elastic behavior, while dynamic balancing conditions are satisfied within an optimization process, based on the whole system response. Besides, the high computational cost of the MTO forced us to propose the use of parallel computing in order to achieve results in a reasonable amount of time by efficiently exploiting the resources of the hardware. Therefore, in this thesis a general methodology for reactionless robot design based on multibody topology optimization is proposed. Then, this methodology is applied to the design of a reactionless four-bar linkage and the partial dynamic balancing of the five-bar mechanism. The relevance of the reactionless four-bar linkage is based on the possibility of using this linkage as modular leg to building reactionless robots. Regarding to the five-bar mechanism, its industrial application gives it enough importance. In addition, the proposed optimization process allows to include structural performance indices in order to conduct an integrated optimization at an early stage of the design process.

Chapter 2

Problem formulation for multibody topology optimization

There is nothing noble in being superior to your fellow man; true nobility is being superior to your former self.

Ernest Hemingway.

T HE mathematical model for topology optimization of multibody systems is described in this chapter. Topology optimization uses the finite element method to modeling the body elastic behavior. The optimization problem of any multibody system involves the complete finite element model of n bodies, connected by kinematic joints under boundary conditions [START_REF] Géradin | Flexible Multibody Dynamics: A finite Element Approach[END_REF]Cardona, 2001, Shabana, 2005). This problem generates a large model, which requires an efficient problem formulation and suitable computational tools in order to find a solution in a reasonable amount of time.

X 0 Y 0 Z 0 Body B i element i, j O 0 O i x i y i z i

Introduction

In the next sections, the elastostatic and elastodynamic models of single body and multibody systems are exposed. The proposed methodology uses finite element formalism to model the individual bodies which are connected by ideal rigid joints. Additionally, model reduction techniques are employed in order to reduce the computational cost of solving the multibody model. This approach is based on [START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF], where the authors address the multibody topology optimization problem.

Elastostatic model

Elastostatic model of a single body

Usually the topology optimization problem is formulated for a single body, but in our case it is necessary to model a system composed for multiple bodies connected by kinematic joints. In order to describe the multibody system, we use subscripts to identify the element and the body. Thus the element j of the body i is defined by the subscript i j. Figure 2.1 is used to represent a general body B i , from where we can notice that the body has its own reference frame O i x i y i z i , and it is discretized in a regular mesh with N i elements.

Standard topology optimization is based on the finite element method, and it uses finite elements of regular shape, where each element is associated with an artificial density, which is the design variable. The design variable is bound to take values between 0 and 1, however in order to avoid optimization results with several intermediate material densities, it is necessary to use a material interpolation scheme [START_REF] Bendsøe | Topology Optimization -Theory, Methods and Applications[END_REF]. In this case, we use the modified simplified isotropic material with penalization (SIMP) interpolation scheme [START_REF] Sigmund | Morphology-based black and white filteres for topology optimization[END_REF]. This method consists in assigning to each finite element a density ρ i j that determines its Young's modulus E i j .

Therefore the modulus of elasticity for element j belonging to body B i is given by:

E i j = E min + ρ p i j (E 0 -E min ) , with ρ i j ∈ [0, 1]. (2.1)
where p is the penalization power (usually p = 3, for elasticity problems), E min is the stiffness at ρ i j = 0 (E min is different from 0 in order to avoid singularity of the stiffness matrix), and E 0 is the Young's modulus of the material.

By using finite element discretization, and resorting to the theory of linear elasticity, the stiffness matrix K i j (ρ i j ) of a single element is written under the SIMP scheme as follows:

K i j = E i j K (0) i j = E min + ρ p i j (E 0 -E min ) K (0) i j (2.2)
where

K (0)
i j is a constant stiffness matrix for an element with Young's modulus equal to one. Therefore, the potential elastic energy (strain energy) of the element i j is:

U e i j = 1 2 u T i j K i j u i j (2.3)
where u i j is the vector of the element i j nodal displacements. Hence, the total potential elastic energy of the body B i is equal to:

U e i = N i ∑ j=1 U e i j = 1 2 N i ∑ j=1 u T i j K i j u i j = 1 2 u T itot K itot u itot (2.4)
where:

•

u itot = u T i1 . . . u T iN i
T , and it is the vector of nodal displacements of the body B i .

• K itot is a block-diagonal matrix stacking on its diagonal all elementary stiffness matrices as follows:

K itot =        K i1 0 . . . 0 K iN i        (2.5)
Accordingly, the stiffness matrix K i of the body B i can be obtained by taking into account the fact that the nodal displacements of the element i j are equal to the nodal displacements of its adjacent elements, as it is usual in standard methods for assembling finite elements [START_REF] Zienkiewicz | The Finite Element Method: Its Basis and Fundamentals[END_REF]. The expression of the vector u itot can be obtained from a reduced set of independent coordinates u i :

u itot = A i u i (2.6)
where A i is a constant matrix, and it is the matrix which generates the assembled stiffness matrix K i . Introducing (2.6) into (2.4) we obtain:

U e i = 1 2 u T i K i u i (2.7)
with matrix K i defined as follows:

K i = A T i K itot A i (2.8)
Thus, K i is the body stiffness matrix which relates the nodal displacements u i to the forces f i exerted on the nodes by the relation [START_REF] Shabana | Dynamics of Multibody Systems[END_REF]:

f i = ∂U e i ∂ u i = K i u i (2.9)
Equation (2.9) shows an explicit relation between the external loads acting in a single body and its nodal displacements. Normally, the nodal displacements are the variables to be computed, and the accuracy of the displacements relies upon body discretization. Nonetheless, the size of the stiffness matrix K i depends on the number of elements in which the body is discretized: the larger the size, the longer the computation time. Thereby it is important to define the best trade-off between the mesh resolution and the computational cost.

The complete elastostatic model of the MBS can be derived using the body stiffness model given in Eq. (2.9). Nevertheless, if the formulation is used in this basic form, it generates a very large MBS stiffness matrix (in our case studies the dimension of K i is typically greater than 10 5 × 10 5 ). Therefore, the use of a classical approach is not adequate for efficient computations, and it is necessary to apply a model reduction technique.

Model reduction technique: static condensation

Using a model reduction technique [START_REF] Qu | Model order reduction techniques with applications in finite element analysis[END_REF] is possible to reduce the computational cost of solving the elastostatic model. The reduction technique uses interface nodes: they are artificial nodes which control the surrounding nodes, in our case by a rigid connection. Basically, the rigid

X 0 Y 0 element i, j Z 0 Body B i Body B (i+1) interface nodes O 0 element (i + 1), j
connection implies that the translations and rotations of the interface node are transmitted directly to its surrounding nodes.

The use of interface nodes allows to connect several bodies. In Fig. 2.2 it is shown the schematic representation of the interface nodes used for bodies connection. Each interface node has 3 DOFs in the planar case (two Cartesian translations and one rotation), and 6 DOFs for the three-dimensional case (three Cartesian translations and three rotations). Moreover, since the MBS under study has bodies connected only in two locations, each body will have two interface nodes. The vector of interface nodes for body B i can be expressed as a concatenation of each individual joint:

u il =    u (1) il u (2) il    (2.10)
where the superscript indicates the body joint. For the planar case u il has size of (6 × 1), and its size is (12 × 1) for the spatial case.

In order to include the interface nodes vector u il as part of the independent set of coordinates, u i , the augmented set of independent coordinates vector is defined as follows:

u a i =    u i u il    (2.11)
Hence, the relation between the interface nodes and its surrounding nodes, as well as the independent set of coordinates, is defined by the following expression:

u a i = J i u c i (2.12)
where vector u c i is the condensed independent coordinates vector, and it is given by:

u c i =    u i f u il    (2.13)
where u i f is the vector of nodes on which no force or displacement is imposed and it has a size of (a × 1). Vector u i f is a subset of u i , and it excludes the nodes which surround the interface nodes.

The size of interface nodes vector u il is defined as (b × 1), and usually a ≫ b.

The matrix J i in Eq. (2.12), is a transformation matrix given by:

J i =    J (11) i J (12) i 0 I    (2.14)
where:

• J (11) i
is the matrix which relates the nodes u i f with the nodes from u i . Its size is

(d × a),
where d is the size of the independent nodes vector u i .

• J

(12) i

is the matrix which relates the interface nodes u il with their connected nodes from u i .

Its size is (d × b).

• 0 is the matrix zero, because the interface nodes u il are independent from u i f , this matrix has a size of (b × a).

• I is an identity matrix. It represents a direct relation between the interface nodes, and its size is (b × b).

In order to compute the potential elastic energy of body B i based the augmented set of independent coordinates vector given in Eq. (2.11), we have to rewrite the body stiffness matrix K i coming from Eq. (2.9). For this purpose we define matrix K a i , which is the matrix K i rewritten as follows:

K a i =    K i(d×d) 0 (d×b) 0 (b×d) 0 (b×b)    (2.15)
Hence, the potential elastic energy of body B i based on the model reduction technique is given as:

U e i = 1 2 u T a i K a i u a i (2.16)
Introducing (2.11) in (2.16) we obtain the following expression:

U e i = 1 2 u T c i K c i u c i (2.17)
where matrix K c i is the condensed body stiffness matrix, and it is given by K c i = J T i K a i J i . This matrix is obtained by the following computations, and it can be rewritten as a block matrix:

K c i =          J (11) i T K i J (11) i (a×a) J (11) i T K i J (12) i (a×b) J (12) i T K i J (11) i (b×a) J (12) i T K i J (12) i (b×b)          =    K (11) c i K (12) c i K (21) c i K (22) c i    (2.18)
On the other hand, the arrangement of the condensed body stiffness matrix affects the vector f i .

Thus, considering the absence of external loads and displacement on the nodal coordinates u i f , the loads vector is rewritten as:

f c i =    0 f il    (2.19)
in which 0 is a zero vector of size (a × 1), which represents the absence of interaction between the environment and the nodes parameterized by u i f . The size of f il is (b × 1).

Taking into account the model reduction technique that we have described, the relation between the nodal displacements u c i to the forces f c i , is given by:

f c i = ∂U e i ∂ u c i = K c i u c i (2.20)
In order to solve the condensed body stiffness model, the linear system

f c i = K c i u c i given in Eq.
(2.20), it can be written as:

   0 f il    =    K (11) c i K (12) c i K (21) c i K (22) c i       u i f u il    (2.21)
from where two equations can be obtained:

0 a×1 = K (11) c i u i f + K (12) c i u il (2.22) f il = K (21) c i u i f + K (22) c i u il (2.23)
Solving Eq. (2.22) for u i f we get:

u i f = Φ Φ Φ si u il (2.24)
where, Φ Φ Φ si is a matrix whose size is (a × b) and, it is the matrix of the static modes [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF] given by:

Φ Φ Φ si = -K (11) c i -1 K (12) c i (2.25)
The force exerted on the interface nodes is obtained by combining Eq. (2.23) and (2.24), thus:

f il = K red i u il (2.26)
where K red i is the reduced body stiffness matrix with size of (b × b) and it is given by:

K red i = K (21) c i Φ Φ Φ si + K (22) c i (2.27)
Besides, the body potential elastic energy from Eq. (2.7) can be rewritten as:

U e i = 1 2 u il T K red i u il .
(2.28) from which can be realized that:

• K red i is the reduced stiffness matrix associated with the displacements of the interface nodes of the body B i .

• K red i has dimensions of (6 × 6) for 2D and (12 × 12) for 3D, and it does not depend on the MBS configuration.

• The most time consuming process is related to the computation of matrix Φ Φ Φ si , Eq. (2.25):

Φ Φ Φ si = -K (11) c i -1 K (12) c i
which requires either the inversion of the matrix K (11) c i or solving the linear system of equations

K (11) c i Φ Φ Φ si = -K (12) c i . If the size of K (11) c i
is small, the matrix of the static modes can be computed using a direct solver. When K (11) c i is large, as is usual in topology optimization of three-dimensional bodies, the computation of this matrix must be done using a iterative solver alongside a suitable preconditioner [START_REF] Borrvall | Large-scale topology optimization in 3D using parallel computing[END_REF].

The advantage of using the model reduction technique is quite remarkable. The typical body stiffness matrix usually has a considerable size which depends on the discretization mesh, then it is compacted into the reduced body stiffness matrix K red i with standard size of (6 × 6) for 2D and (12 × 12) for 3D, thus leading to a small stiffness matrix for the assembled MBS. Moreover, the computation of the matrices Φ Φ Φ si , which is intensive (computationally speaking), can be made in parallel in order to save computational time.

Elastostatic model of the MBS

Once the elastostatic model for single body is computed it is possible to compute the MBS elastostatic multibody model. This multibody model considers the MBS configuration, and the boundary conditions.

In order to take into account the MBS configuration, the orientation of each body should be considered. Since the reduced body stiffness matrix K red i is expressed in its own local reference frame, it must be expressed in the global frame, which is done using the coordinates transformation.

Coordinate transformation is performed by using a block-diagonal matrix Q i :

Q i =        R i 0 . . . 0 R i        (2.29)
where R i is a rotation matrix with size of (3 × 3), which represents the body orientation, and it comes from the solution of MBS kinematics. Certainly R i depends on the MBS configuration q, but not on nodal displacements (assuming small perturbations). Usually the matrix Q i has size of (6 × 6) when the bodies are considered as planar ones, and (12 × 12) when the bodies are modeled as three-dimensional (if two interface nodes per body are used).

Therefore, the reduced body stiffness matrix is expressed in the global coordinate system as follows:

K red i 0 = Q i K red i Q T i .
(2.30)

Considering a MBS made of n bodies, as seen in the Fig. 2.2, the full potential elastic energy of the system is given by:

U e = n ∑ i=1 U e i = 1 2 n ∑ i=1 (u il ) T 0 K red i 0 (u il ) 0 = 1 2 u red tot T K red tot u red tot (2.31)
where:

• u red tot = (u il ) T 0 . . . (u nl ) T 0
T is the vector of interface nodes displacements for all n bodies in the global frame.

• K red tot is a block-diagonal matrix stacking on its diagonal all bodies stiffness matrices as follows:

K red tot =        K red i 0 0 . . . 0 K red n 0       
.

(2.32)

The MBS stiffness matrix K r can be obtained by taking into account the fact that the bodies are connected altogether through the interface nodes. As a result, the expression of the vector u red tot can be obtained from a reduced set of independent coordinates u r , as follows:

u red tot = J r u r (2.33)
where J r depends on how the bodies of the MBS are connected, namely their kinematic pairs.

The kinematic relations between the interface nodes depend on the joint type. In the case of the revolute joint, and considering a rigid connection, all displacements are constrained to be the same, and only the rotation about the joint axis is independent. Thus, considering two adjacent bodies, B i and B (i+1) as we see in Fig. 2.2, we have:

u red tot = u il u (i+1)l T = u (1) il u (2) il u (1) (i+1)l u (2) (i+1)l T (2.34)
If these two bodies are connected by the joints u

(2)

il and u

(1) (i+1)l , and choosing u

(1) (i+1)l as the dependent coordinates set, the reduced set of independent coordinates is given by:

u r = u (1) il u (2) il θ (i+1)r u (2) (i+1)l T (2.35)
where θ (i+1)r is the independent coordinate associated with body B (i+1) rotation. Furthermore, the boundary conditions related to the MBS supports can be defined using the remaining interface nodes in the reduced set of independent coordinates vector.

Introducing (2.33) into (2.31), the potential elastic energy of the MBS can be computed by:

U e = 1 2 u T r K r u r (2.36)
where the MBS stiffness matrix K r is:

K r = J T r K red tot J r (2.37)
Moreover, the relation between the nodal displacements u r and external forces f r exerted on the considered nodes is given by:

f r = ∂U e ∂ u r = K r u r (2.38)
For a MBS with n bodies, the size of K r is lower than (6n × 6n) in the 2D case, and for 3D case the size is lower than (12n × 12n). Equation (2.38) describes the relation between the external loads acting on the MBS and its joint displacements. By solving this model it is possible to compute the joint displacements in the MBS when it is subjected to external loads. As it is usual, this linear system is solved by applying the boundary conditions defined by the supports of the multibody system. This structural response is used in the optimization process.

Computing the nodal displacements for the elastostatic model

Equation (2.38) represents a linear system from where the vector u r is computed. Then, it is necessary to compute the complete vector of the reduced set of interface nodes using the relation given in (2.33), namely u red tot = J r u r .

Since the vector of the reduced set of interface nodes

u red tot = (u il ) T 0 . . . (u nl ) T 0 T , is formed by
all the interface nodes of each body B i , it is necessary to split the vector on its corresponding interface nodes vectors (u il ) 0 . These vectors are described in the global reference frame, thus it is necessary to express them in their local frames. This is done using the block-diagonal matrix Q i :

u il = Q T i (u il ) 0 (2.39)
Then, with the body interface displacements u il defined in its local reference frame, it is possible to compute the complete nodal displacements of the body by means of Eq. (2.24), namely

u i f = Φ Φ Φ si u il .
Thus, the condensed independent coordinates vector u c i is given as follows:

u c i =    u i f u il    =    Φ Φ Φ s u il u il    (2.40)
Besides, the augmented set of independent coordinates vector is given by the Eq. (2.11), u a i = J i u c i , therefore:

   u i u il    =    J (11) i J (12) i 0 1       Φ Φ Φ s u il u il    =    J (11) i Φ Φ Φ s u il + J (12) i u il u il    (2.41)
Finally, the independent body nodal displacements for the elastostatic model is given by:

u i = J (11) i Φ Φ Φ s u il + J (12) i u il .
(2.42)

The vector displacements given in Eq. (2.42) is a computationally efficient expression based on the model reduction technique.

Elastodynamic model 2.3.1 Elastodynamic model of a single body

The body dynamic response related to natural frequencies is important because these frequencies are associated with the highest level of energy due to vibrations. The following sections describe the methodology to compute natural frequencies.

For the elastodynamic model the mass matrix associated with the element i j is given by:

M i j = ρ i j M 0 i j (2.43)
where M 0 i j is the consistent mass matrix of a single element computed for a density equal to 1. Thus, the kinetic energy of the element due to elastic oscillations is:

T e i j = 1 2 uT i j M i j ui j (2.44)
being ui j the vector of the nodal velocities. Hence, the total kinetic energy of the body B i is equal to:

T e i = N i ∑ j=1 T e i j = 1 2 N i ∑ j=1 uT i j M i j ui j = 1 2 uT itot M itot uitot (2.45)
where M itot is a block-diagonal matrix stacking on its diagonal all elementary mass matrices as follows:

M itot =        M i1 0 . . . 0 M iN i        (2.46)
Besides, differentiating (2.6) with respect to time (matrix A i is constant) the following expression is obtained:

uitot = A i ui (2.47)
Introducing (2.47) into (2.45), we obtain:

T e i = 1 2 uT i M i ui (2.48)
where M i = A T i M itot A i is the mass matrix of body B i . As can be noticed in Eq. (2.43), the material interpolation scheme used for the mass matrix is a linear approach, which is the natural physical option for the mass matrix.

The elastodynamic model for the MBS can be generated directly form the results already described, but this model can be of considerable size. In consequence a model reduction technique is required in order to reduce the computational effort associated to the solution of the complete elastodynamic model for the MBS.

Model reduction technique: dynamic condensation

In order to reduce the computational cost associated with the MBS natural frequencies, and because of deformations remain small, the dynamic substructuring scheme or Craig-Bampton model reduction technique is used. The Craig-Bampton model reduction technique is based on the assumption that body nodal coordinates u i f can be expressed as a function of the nodal coordinates u il , corresponding to the interface nodes plus a term characterizing their vibratory free behavior, as follows:

u i f = Φ Φ Φ si u il + Φ Φ Φ di µ µ µ i (2.49)
where the term Φ Φ Φ si u il comes from (2.24) and it characterizes the node static displacements, and the term Φ Φ Φ di µ µ µ i is an additional term characterizing the body oscillatory behavior.

Classically, matrix Φ Φ Φ di has the following form:

Φ Φ Φ di = u (1) i f • • • u (s) i f (2.50)
where the vector u

(k)
i f is the k-th eigenmode associated with the equation:

M (11) c i üi f + K (11) c i u i f = 0 (2.51) in which matrix K (11) c i is defined in (2.18). Matrix M (11) c i
is obtained following a procedure similar that the used to obtain the condensed body stiffness matrix. With J i given in Eq. (2.14), the condensed body mass matrix is given by:

M c i =          J (11) i T M i J (11) i (a×a) J (11) i T M i J (12) i (a×b) J (12) i T M i J (11) i (b×a) J (12) i T M i J (12) i (b×b)          =    M (11) c i M (12) c i M (21) c i M (22) c i    (2.52) The dimension of matrix M (11) c i is (a × a).
The matrix Φ Φ Φ di is formed using s vectors, chosen among the a eigenmodes associated with (2.51). Usually, s < a (in our examples, s is lower than 10).

Similarly as in the elastostatic mode, most of the computational cost is due to the computation of the matrix Φ Φ Φ di . Nevertheless, for a multibody system, the computation of this matrix is independent for each body. Now, using (2.49) the condensed nodal displacement vector, u f req c i , for the oscillatory free vi-bration is given by:

u f req c i =    u i f u il    = B i d i , B i =    Φ Φ Φ si Φ Φ Φ di I b 0 b×s    , d i =    u il µ µ µ i    .
(2.53) in which I b is the identity matrix of dimension (b × b). Taking into account once again that matrix B i does not depend on the MBS configuration or the link deformation (hypothesis of small perturbations), the derivative of Eq. (2.53) with respect to time leads to:

u f req c i = B i ḋi (2.54)
By using these transformation of coordinates, the body potential and kinetic elastic energies given at (2.7) and (2.48) becomes:

U e i = 1 2 d T i K f req i d i (2.55
)

T e i = 1 2 ḋT i M f req i ḋi (2.56)
where the reduced body stiffness matrix is:

K f req i = B T i (K c i ) B i , (2.57)
and the reduced body mass matrix is:

M f req i = B T i (M c i ) B i .
(2.58)

The reduced body stiffness and mass matrices are usually of small dimension. Typically their sizes are {(6 + s) × (6 + s)} for the planar case, and {(12 + s) × (12 + s)} for spatial case, where s is the number of chosen eigenmodes. These matrices not depend on the MBS configuration, which is a great advantage because they must be computed only once at each step of the optimization algorithm, whatever the number of tested configurations for the MBS. Moreover, after simplifying the computations of the reduced body stiffness matrix given by (2.57), we obtain [START_REF] Géradin | Flexible Multibody Dynamics: A finite Element Approach[END_REF]:

K f req i =    K (21) c i Φ Φ Φ si + K (22) c i 0 0 Φ Φ Φ T di K (11) c i Φ Φ Φ di    (2.59)
Conversely, the reduced body mass matrix given by (2.58), can be simplified as:

M f req i =    M (21) c i + Φ Φ Φ T si M (11) c i Φ Φ Φ si + M (22) c i + Φ Φ Φ T si M (12) c i M (21) c i + Φ Φ Φ T si M (11) c i Φ Φ Φ di Φ Φ Φ T di M (12) c i + M (11) c i Φ Φ Φ T si Φ Φ Φ T di M (11) c i Φ Φ Φ di    (2.60)
where the block reduced body mass matrix M f req i 21) is the transpose of the block matrix M f req i 12) .

(
(
Equations (2.59) and (2.60) are explicit expression, useful for an efficient computation of the reduced body matrices.

Elastodynamic model of the MBS

In order to define the multibody elastodynamic model of a system, it is necessary to take into account the orientation of each body, it means the multibody system configuration. It can be performed using the matrix Q

f req i
, which is a block-diagonal matrix generated by staking on its diagonal the rotation matrix R i(3×3) which represents the body B i orientation, and it comes from the solution of MBS kinematics:

Q f req i =            R i . . . 0 . . . . . . . . . R i 0 . . . I (s×s)            (2.61)
where I (s×s) is an identity matrix.

Usually the matrix Q f req i has size of {(6 + s) × (6 + s)} for the planar case, and {(12 + s) × (12 + s)} for spatial case, where s is the number of chosen eigenmodes. For a single configuration there is only one matrix Q

f req i
for each body, but for a set of trajectories there will be as many as configurations the MBS reaches. Thus, Q f req i depends on the MBS configuration q but not on the nodal displacements u f req i

(assumption of small perturbations).

In order to take into account the configuration of the multibody system, the orientation of each body must be considered. In case of a link with a planar motion, this is accomplished by the following matrix transformations:

K f req i 0 = Q f req i K f req i Q f req i T (2.62) M f req i 0 = Q f req i M f req i Q f req i T (2.63)
Considering now a MBS composed of n bodies, the full potential and kinetic energies of the system are given by the summation of individual energies from Eq. (2.55) and Eq. (2.56), then:

U e = n ∑ i=1 U e i = 1 2 n ∑ i=1 (d i ) T 0 K f req i 0 (d i ) 0 = 1 2 u f req tot T K f req tot u f req tot
(2.64)

T e = n ∑ i=1 T e i = 1 2 n ∑ i=1 ḋi T 0 M f req i 0 ḋi 0 = 1 2 u f req tot T M f req tot u f req tot (2.65)
where:

• u

f req tot = d T 1 0 . . . d T n 0
T is the vector composed of vectors (d i ) 0 for all n MBS bodies, defined in the global frame.

• K f req tot and M f req tot are block-diagonal matrices stacking on their diagonal all bodies stiffness and mass matrices as follows:

K f req tot =        K f req 1 0 0 . . . 0 K f req n 0        , M f req tot =        M f req 1 0 0 . . . 0 M f req n 0        . (2.66)
The reduced elastodynamic model of the MBS can be obtained by taking into account the fact that the bodies are connected through the interface nodes. As a result, the expression of the vector u f req tot can be obtained from the MBS reduced set of independent coordinates u d as follows

u f req tot = J d u d (2.67)
Moreover, for the analysis of the oscillatory free behavior, the matrix J d is constant and it depends on the MBS kinematic relations, thus:

u f req tot = J d ud (2.68)
Introducing (2.67) and (2.68) into (2.64) and (2.65), we get

U e = 1 2 u T d K d u d , T e = 1 2 uT d M d ud (2.69)
where the reduced MBS stiffness matrix, K d , and the reduced MBS mass matrix, M d , are computed as follows:

K d = J T d K f req tot J d , M d = J T d M f req tot J d (2.70)
Using the Lagrange equations, in absence of external efforts, we can express:

d dt ∂ L ∂ u - ∂ L ∂ u = 0 (2.71)
where L = T e -U e , the reduced dynamic equation characterizing the MBS free oscillations is:

M d üd + K d u d = 0 (2.72)
A solution u dk of this equation satisfies:

ω 2 k M d -K d u d k = 0 (2.73)
where the pulsation ω k of the MBS is equal to ω k = 2π f k . Besides, f k is the natural frequency associated with the k-th natural mode of vibrations, and u d k is its associated eigenvector. Indeed, for the solution of the Eq. (2.73) the boundary conditions must be applied, affecting both matrices

M d and K d .
For a 2D problem, the typical size of matrices K d and M d is lower than {(6 + s)n × (6 + s)n}, being n the number of bodies in the MBS, and s the number of chosen eigenmodes. Besides, for 3D problems the size is given by: {(12 + s)n × (12 + s)n}. The resolution of Eq. (2.73) is considerably more efficient than the classical solutions due to the applied model reduction techniques.

Computing the nodal displacements for the elastodynamic model

Once the associated vibration mode is computed, the resultant eigenvector u d,1 is used in order to compute the displacements due to the first natural vibration.

u f req tot = J d u dk
(2.74)

The vector u f req tot , defined in Eq. (2.64), contains all the joints displacements of the multibody system. Therefore, this vector must be splitted in their corresponding displacements bodies, thus each vector (d i ) 0 corresponds to the joint displacement of the body B i . Moreover, since the vector (d i ) 0 is expressed in the global reference frame, it is necessary to express it in the local frame, this is performed using the block-diagonal matrix Q f req i :

d i = Q f req i T (d i ) 0
(2.75) thus the vector d i is expressed in its local reference frame. Besides, the independent condensed displacement vector is computed as u

f req c i = B i d i , hence: u f req c i =    Φ Φ Φ si Φ Φ Φ di I b 0 b×s       u il µ µ µ i    =    Φ Φ Φ si u il + Φ Φ Φ di µ µ µ i u il    (2.76)
The vector u f req c i corresponds to the independent nodes in the body B i , then the augmented nodal displacements it is computed as:

u f req a i = J i u f req c i (2.77) =    J (11) i J (12) i 0 I       Φ Φ Φ si u il + Φ Φ Φ di µ µ µ i u il    (2.78) =    J (11) i (Φ Φ Φ si u il + Φ Φ Φ di µ µ µ i ) + J (12) i u il u il    (2.79)
The body nodal displacements due to the mode of vibration u dk is given by:

u f req i = J (11) i (Φ Φ Φ si u il + Φ Φ Φ di µ µ µ i ) + J (12) i u il (2.80)
Usually the first natural frequency is used as a performance index, and the eigenvector u d1 , associated to this frequency is used as initial displacement for the modal analysis.

Topology optimization based on SIMP scheme is prone to numerical problems such as meshdependence and checkerboard patterns, which are well documented and successful solutions have been applied [START_REF] Diaz | Checkerboard patterns in layout optimization[END_REF][START_REF] Sigmund | Morphology-based black and white filteres for topology optimization[END_REF][START_REF] Sigmund | Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards mesh-dependencies and local minima[END_REF]. In order to overcome these undesired effects it is essential to use restriction methods, which are usually applied as a filters. Thereby, among the broad options of filters, we decided to use one of the most used types of density filters, which is described below.

Density filter

Density filtering was introduced in [START_REF] Bruns | Topology optimization of non-linear elastic structures and compliant mechanisms[END_REF] and the mathematical proof of the solution existence was reported in [START_REF] Bourdin | Filters in topology optimization[END_REF]. This filter transforms the original densities ρ i j as follows:

ρi j = ∑ k∈N e H ek ρ ik ∑ k∈N e H ek (2.81)
where N e is the set of elements k for which the center-to-center distance D(i j, k) to element i j is smaller than the filter radius r min and H ek is a weight factor defined as:

H ek = max (0, r min -D(i j, k)) (2.82)
where the value of the filter radius r min in our problems is from 2 to 3 times the element size.

Besides, usually the original densities ρ i j are referred to as the design variables, and the filtered densities ρi j are referred to as the physical densities.

When a density filter is applied, the original densities ρ i j do not have any physical meaning [START_REF] Sigmund | Morphology-based black and white filteres for topology optimization[END_REF]. Therefore, the filtered density field ρi j should be always presented rather than the original density field ρ i j as the solution to the optimization problem. Additionally, the sensitivities of the objective function f (ρ ρ ρ) and constraints functions must be modified replacing the variable ρ i j with ρi j . The sensitivities with respect to the design variables ρ i j are obtained by means of the chain rule. Moreover, as it was pointed out by [START_REF] Sigmund | Morphology-based black and white filteres for topology optimization[END_REF], it is necessary to have special treatment for the boundaries of the structure. Therefore the supported boundaries are specified as solid elements (region) in order to obtain adequate filtering.

In the previous sections the formulation for multibody topology optimization was developed. The multibody system was described in terms of the design variables, but the inertial properties of the individual bodies must be also rewritten in terms of the design variables. Hence, those inertial parameters can be properly incorporated into the topology optimization formulation as it is shown in the next section.

Inertial parameters of an individual body

Each body of the MBS can be characterized by its inertial parameters such as its mass, the first moments of inertia (also known as the static moments), the second moments of inertia, and the products of inertia. The inertia parameters are defined in a local frame rigidly attached to a body,

O 0 Y 0 X 0 y i M ij m ij Body B i O i Z 0 z i x i r i
and they do not depend on the MBS configuration.

From [START_REF] Khalil | Modeling, Identification and Control of Robot[END_REF], we have the formulas for the computation of the body inertial parameters. Nevertheless, these formulas must be rewritten as a function of the design variables, i.e., the element density ρ i j . Because of the bodies are modeled using finite elements, each body B i is discretized in a total of N i elements, and the element j of the body i is defined by the subscript i j. The mass of the element i j is defined as m i j , and x i j , y i j and z i j are the positon of the origin M i j of this element in its local frame (O i , x i , y i , z i ) attached to the body B i . In Fig. 2.3 the schematic description of the finite element i j in the body B i is represented. As a consequence, considering a linear relation between the mass and the design variables, the following expressions for the inertial parameters of body B i are obtained: Body's mass:

m i = N i ∑ j=1 m i j ρ i j
(2.83a)

Static moments:

I xi = N i ∑ j=1 m i j x i j ρ i j , I y i = N i ∑ j=1 m i j y i j ρ i j , I z i = N i ∑ j=1 m i j z i j ρ i j (2.83b)
Moments of inertia:

I xxi = N i ∑ j=1 m i j y 2 i j + z 2 i j ρ i j , I yy i = N i ∑ j=1 m i j x 2 i j + z 2 i j ρ i j , I zz i = N i ∑ j=1 m i j x 2 i j + y 2 i j ρ i j (2.83c)
Products of inertia:

I xy i = N i ∑ j=1 m i j x i j y i j ρ i j , I xz i = N i ∑ j=1 m i j x i j z i j ρ i j , I yz i = N i ∑ j=1 m i j y i j z i j ρ i j (2.83d)
where Eq. (2.83a) gives the total mass of the body, while Eqs. (2.83b), (2.83c) and (2.83d) are the corresponding static moments, the moments of inertia and the products of inertia, respectively.

Once the formulation of the multibody topology optimization problem has been properly de-rived, the next stage is its resolution. In this thesis, the optimization problem is solved by the MMA optimizer in all the study cases. Considering the MBS compliance and its natural frequency as important performance indices, the analytical expression of the gradients of these functions are derived in the following sections.

Optimizer and sensitivity analysis

The problem we face is characterized by a high number of variables. In order to solve the optimization problem it is necessary to use an optimization solver, able to handle large number of variables, as well as multiple constraints. Thus, the MMA algorithm introduced in [START_REF] Svanberg | The method of moving asymptotes -a new method for structural optimization[END_REF] is the optimizer used for this research work. It is able to handle multiple inequalities and equalities constraints, and it is the most accepted optimizer in the structural optimization community because of its excellent convergence properties [START_REF] Bendsøe | Topology Optimization -Theory, Methods and Applications[END_REF]. A special attention is paid to the equality constraints, which can be converted into two inequality constraints as it is shown in the section below. In addition, this method is a first-order gradient-based algorithm which requires the derivatives of the objective function and of the constraints with respect to the design variables, therefore a sensitivity analysis is conducted.

Equality constraints

As it is usual in the constrained optimization problem, the constraints can be defined as inequality or equality constraints. Equality constraints can be transformed into inequalities without loss of generality, in order to avoid a very restrictive optimization scheme. It is done by replacing the equality constraint h i (ρ ρ ρ) = 0, with two inequality constraints, h i (ρ ρ ρ) -ε ≤ 0 and -h i (ρ ρ ρ) + ε ≤ 0, where ε is usually a small number called relaxation parameter.

The parameter ε can be defined as a constant value, or it can be considered as a variable parameter able to change throughout the optimization process. The way to define the value of this parameter, is typically based on a heuristic approach. In our case we decided to use ε as a decreasing parameter in order to be closer to the equality constraint at each iteration step. In order to update the relaxation parameter we use an asymptotic function defined by:

ε = ε 0 (κ + 10) β (2.84)
where κ is the iteration number, and β is a parameter defined by the user, as well as ε 0 . In our case we define ε 0 = 0.9 and β = 3.5 This relaxation parameter is used for all the cases analyzed in the thesis.

Sensitivity analysis

Sensitivity analysis or gradient computation is the centerpiece of the gradient-based optimization methods. Hence, in the ideal case the analytical expression of the gradients are preferred instead of the numerical approximations. Commonly, the numerical methods, such as finite differences (or more robust methods) are used when the analytical expressions are not available. The problem with the numerical methods is the high computational cost. Nevertheless, for the problems we face in this thesis, the gradients of the performance indices, such as the MBS compliance and its natural frequency can be computed analytically and they are derived in the sections below.

With respect to the body's inertial parameters, the gradients computation is simple because of the linear relation between the inertial parameters and the design variables. Thus, the gradients can be easily computed analytically.

Compliance

A standard objective function or performance index in topology optimization is the compliance [START_REF] Bendsøe | Topology Optimization -Theory, Methods and Applications[END_REF], which usually is considered as the inverse of the stiffness. The compliance as it is defined in [START_REF] Sigmund | Design of Material Structures Using Topology Optimization[END_REF] can be considered as twice the strain energy (potential elastic energy), or also as the dot product of the nodal wrenches by the nodal displacements.

This performance index will be used in this work and its gradient will be derived. Because of the displacements computed using the static condensation are completely equivalent to those values computed using the full elastic model, the compliance gradient can be derived from the general elastostatic model of the MBS, without involving the model reduction technique. Consequently, the compliance function based on the full model for the MBS is given by:

c(ρ ρ ρ) = u T mbs K mbs u mbs (2.85)
where:

• u mbs = u T 1 . . . u T n is the vector stacking all nodal displacements for all n bodies, and each vector u i can be computed using the expression shown in (2.42).

• K mbs is a block-diagonal matrix stacking on its diagonal all bodies stiffness matrices:

K mbs =        K 1 0 . . . 0 K n        (2.86)
Therefore, the compliance gradient is given by:

∂ c(ρ ρ ρ) ∂ ρ ρ ρ = u T mbs ∂ K mbs ∂ ρ i j u mbs (2.87)
where the partial derivative of K mbs is:

∂ K mbs ∂ ρ i j =                0 . . . ∂ K i ∂ ρ i j . . . 0               
, and

∂ K i ∂ ρ i j = -pρ p-1 i j E 0 A T i                0 . . . K (0) i j . . . 0                A i .
(2.88) where the expression for K i is defined in (2.8).

Natural frequency

Other important performance indices to be optimized in multibody systems are the natural frequencies, which are significant in order to avoid serious vibrations when the system is affected by dynamic loads that are close to these natural frequencies. Therefore, usually the first natural frequency or fundamental frequency is maximized in order to prevent this undesirable situations. In fact, optimize the fundamental frequency wont reduce the dynamic unbalance in a MBS by itself, but optimizing this frequency when shaking force balancing by mass redistribution is applied, will be beneficial [START_REF] De Jong | The influence of flexibility on the force balance quality: a frequency domain approach[END_REF].

The first natural frequency of the MBS is given as follows:

f (ρ ρ ρ) = f 1 = ω 1 2π (2.89)
where ω 1 is the first natural pulsation obtained from the solution of Eq.

(2.73):

ω 2 1 M d -K d u d 1 = 0.
Moreover, the natural frequency and the pulsation are related by ω 1 = 2π f 1 . [START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF] derived the gradient of the MBS natural frequency, and it is given as follows:

∂ f 1 ∂ ρ i j = 1 4πω 1 u T d 1 M d u d 1 u T all ∂ K all ∂ ρ i j u all -ω 2 1 u T all ∂ M all ∂ ρ i j u all (2.90)
where the involved terms are:

∂ K all ∂ ρ i j =                0 . . . ∂ K itot ∂ ρ i j . . . 0               
, and

∂ M all ∂ ρ i j =                0 . . . ∂ M itot ∂ ρ i j . . . 0                (2.91) ∂ K itot ∂ ρ i j = -pρ p-1 i j E 0                0 . . . K (0) i j . . . 0               
, and

∂ M itot ∂ ρ i j =                0 . . . M (0) i j . . . 0                (2.92) with K (0)
i j and M (0) i j defined in (2.2) and (2.43), respectively. In addition, matrix M d was derived in (2.70). Furthermore, the vector u all stacks all nodal displacements for all elements in the multibody system, and it is given by:

u all = u T 1tot . . . u T ntot T (2.93) being u itot = A i u f req i
, with A i defined in (2.6), and u

f req i
given in (2.80).

Topology optimization for three dimensional multibody systems leads to large-scale analysis, which makes imperative to solve the problem using a suitable computational tool. Indeed, an efficient modeling formulation is also necessary to reduce the computational effort. For that reason model reduction techniques have been proposed in Section 2.2.2 and 2.3.2. These strategies by themselves are not enough to reduce the high computational cost associated to the solution of the problem. Therefore, in order to solve the problem at hand in a reasonable amount of time, a computational platform for fully parallel processing is developed. In the following Section is presented an overview of the computational platform.

Computational platform

The structural topology optimization is based on finite element method. The solution of equilibrium equations coming from the finite element model is the most time-consuming stage. [START_REF] Borrvall | Large-scale topology optimization in 3D using parallel computing[END_REF] report that 97% of the computational time is used in the solution of the resulting linear system.

Topology optimization problems involves a large amount of design variables, and when a three dimensional multibody topology optimization problem is considered, the computational effort required to solve the problem increases drastically. On the other hand, the mesh size has an important role on the accuracy of results, usually a fine mesh is desired to get more realistic designs. Due to these conditions, the optimization problem to be solved implies a large-scale analysis. For reactionless robot design, in order to fulfill balancing conditions accurately, a fine mesh is always desired, since the results can be affected by the mesh discretization.

Even though MatLab is an excellent software for develop and test algorithms, the problems of medium or large size require specialized programming languages for a better management of the computational resources. A suitable programming language using with the appropriate computing paradigm can exploit the full potential of the given computational system. As a consequence, a computational platform was developed specifically for the solution of multibody topology optimization problems. It is programmed in C++ language, using the object oriented paradigm and the distributed memory model for parallel computing. It is based on the Portable, Extensible Toolkit for Scientific Computation (PETSc) [START_REF] Balay | PETSc Web page[END_REF], and the platform incorporates some classes provided in [START_REF] Aage | Topology optimization using PETSc: An easy-to-use, fully parallel, open source topology optimization framework[END_REF]. The computational platform for multibody topology optimization (MUBOTO) is not software dependent, and it is designed to solve large scale problems. The platform was tested in a desktop computer with a processor AMD Ryzen 9 3900X, using the OS Ubuntu 20.04 64 bits, and PETSc 3.16.0.

Codes for topology optimization

There are several computational codes for topology optimization. Some of the most popular are those written in MatLab [START_REF] Andreassen | Efficient topology optimization in MAT-LAB using 88 lines of code[END_REF][START_REF] Liu | An efficient 3D topology optimization code written in matlab[END_REF][START_REF] Sigmund | A 99 line topology optimization code written in matlab[END_REF]. These programs have an educational purpose, are written for a single body analysis, and they solve classical benchmarks for TO. With respect to multibody topology optimization [START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF] describe the optimization of five-bar mechanism using MatLab. The authors report some computational issues when the number of elements in the FE model start to increase, mainly for a 3D model. In this case the computational cost is increased, and they refer the necessity of specific solver (preconditioned conjugate gradient) for a better performance. Even more, the solution of the eigenvalue problem for the 3D model took several days for a relatively small size problem.

Thus, for this scenario MatLab is not practical, since it limits the optimization size, and becomes prohibitive for systems of bigger size.

The preferred software to develop educational codes for TO is MatLab [START_REF] Zhu | Design of compliant mechanisms using continuum topology optimization: A review[END_REF]. It is because MatLab has an easy to use programming language, and has highly optimized algorithms and executes some of them in parallel1 . Nevertheless, when the problem turns into large-scale analysis MatLab is not the best option, and then it is necessary to use standard scientific programming languages such as Fortran or C/C++, in conjunction with parallel processing techniques.

An overview of parallel computing

Parallel computing is the simultaneous use of multiple processors resources to solve a computational problem [START_REF] Pacheco | An introduction to parallel programming[END_REF], it implies to execute multiples computations at the same time. In contrast, in the serial execution the operations are executed one-at-a-time. Today, most of the desktop computers are build using processors of multiple cores, thus basically we have parallel computers in our homes. In contrast, the first personal computers were packaged with a single Central Processing Unit (CPU) or core.

Traditionally, software has been written for serial computation, designed for a single processor computer. In a typical serial program, instructions are executed sequentially one after another and only one instruction may execute at any moment in time. In fact, nowadays serial programming is still very common, even if the programs are running in multiprocessor desktops. It happens because of the complexity of parallel programming.

In the last years, multiprocessor computers have been available, from personal computers to High Performance Computing (HPC) centers [START_REF] Eijkhout | Introduction to High Performance Scientific Computing[END_REF]. Hence, in order to exploit the hardware resources, the serial programming has been replaced by parallel programming. Parallel computing paradigm is a natural response for the efficient use of multiprocessor hardware. Generally in a parallel program, the problem is broken into discrete parts that can be solved concurrently. Each part is further broken down to a series of instructions, these instructions from each part execute simultaneously on different processors.

Parallel computing requires some additional features in order to exploit the parallel hardware capabilities. In this case there are some parallel programming models, but the most used are the share-memory and distributed-memory models. Parallel programming models exist as an abstraction above hardware and memory architectures. Commonly the share-memory approach is programmed using the standard Open Multiprocessing (OpenMP). For the distributed-memory model, the standard is Message Passing Interface (MPI). It is common to find many different implementations of these standards, which depends on the hardware vendors. OpenMP is an Application Program Interface (API), which is an explicit (not automatic) programming model, and supports C/C++ and Fortran on a wide variety of architectures. This standard is included in many compilers and provides a portable and scalable model for developers of shared memory parallel applications, it uses the multi-thread concept under the fork-join model. Sometimes this API is called masterworker approach because the thread-based coding.

Despite of the different models, the industrial standard for high performance computing is MPI. This library allows to use the real parallel hardware capacities, because it is designed specifically for distribute-memory systems, additionally it is possible to use multi-thread implementations in the same code, alike it supports C/C++ and Fortran programming languages. Additionally, this standard can be executed in a share memory system, such as desktop computer.

Some researcher choose the OpenMP to develop their computational platforms, and usually it is because its simplicity at implementation level. Sometimes they need to use sequential code previously written, and it can be easily modified using OpenMP directives [START_REF] Paris | Parallel computing in topology optimization of structures with stress constraints[END_REF]. Nevertheless, the codes developed based on MPI offers the best scalability and true high performance, but their use requires a bigger effort in order to develop the complete code in parallel.

Message passing interface

Originally, the MPI library (not a language) was designed for distributed memory architectures.

Today, MPI runs on virtually any hardware platform, either distributed-memory, shared-memory or hybrid hardware. The real power of MPI is in being able to use multiple computers. The MPI library is the de facto tool for large scale parallelism as it is used in engineering sciences.

The Message-passing is a communication paradigm in which processes communicate by ex-changing messages via communication channels. The processes is a task executing on a given processor at a given time. The Operating System (OS) treats each process as an independent entity and schedules it to run on system resources. Each process maintains its own virtual address space, which the OS maps into physical memory. In MPI model there are always multiple processes active.

It is possible to start a large number of MPI processes, even on a laptop. In MPI all parallelism is explicit. The programmer is responsible for correctly identifying parallelism and implementing parallel algorithms using MPI constructs. MPI defines their own data-type to produce a portable application. Thus, depending on the system, MPI cast the data to fitting the datatype of the host system. The basic model of MPI is Single Program Multiple Data (SPMD). It means that each process is an instance of the same program, each process operates on its own data. The synchronization between the MPI processes is done through explicit send and receive calls. In MPI applications there is no master process, all processes are equal, start and end at the same time, this property is called symmetry.

Libraries for parallel programming

Use a parallel programming paradigm implies that all tools must be compatible within them. Hence, there are some libraries written for parallel processing which provides capabilities for linear algebra computations and data management. Some functions of these libraries are developed on OpenMP or MPI, and in some cases for both. Commonly, for economical reasons and flexibility, the opensource or free libraries are preferred, and many of them have an excellent performance. Some classical libraries for linear algebra operations, which also provide parallel functionality based on share-memory model are Eigen [START_REF] Guennebaud | Eigen[END_REF] and Armadillo [START_REF] Sanderson | Armadillo: a template-based C++ library for linear algebra[END_REF]. Examples of libraries which support distribute-memory (and also share-memory) are Elemental [START_REF] Poulson | Elemental: A new framework for distributed memory dense matrix computations[END_REF], PETSc [START_REF] Balay | PETSc Web page[END_REF] and Trillinos [START_REF] Heroux | An overview of the Trilinos project[END_REF]. Additionally, there are libraries specifically focused on the solution of large linear systems of equations (linear solvers), as is the case of PARDISO [START_REF] Kourounis | Towards the next generation of multiperiod optimal power flow solvers[END_REF], MUMPS [START_REF] Amestoy | MUMPS: A general purpose distributed memory sparse solver[END_REF], and HYPRE [START_REF] Falgout | hypre: A library of high performance preconditioners[END_REF]. On the other hand, in many cases the results from large-scale parallel processing must be treated with suitable tools. These tools must be able to display or process huge amount of data. Some packages which provides these benefits are VisIt [START_REF] Childs | High Performance Visualization-Enabling Extreme-Scale Scientific Insight, chapter VisIt: An End-User Tool For Visualizing and Analyzing Very Large Data[END_REF], TecPlot Software and ParaView [START_REF] Ahrens | The Visualization Handbook, chapter 36[END_REF]. Besides, because of the use of the FEM, the need for a meshing tool capable to have a good intregration with the processing libraries is important. It is possible to find some libraries such as LibMesh [START_REF] Kirk | libMesh: A C++ Library for Parallel Adaptive Mesh Refinement/Coarsening Simulations[END_REF] and Gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF] which provide several tools for mesh generation and usually are compatible with the most used solvers.

In recent years, the HPC has been improved by the use of accelerator-based architectures like Graphics Processing Units (GPU) and Field Programmable Gate Arrays (FPGA), among others.

These architectures has the purpose to accelerate the computations, but their use is not straightforward and depends on the type of application. In general, because of the differences in their architectures each device uses its own programming language, producing a complicated integration between them.

Improving the performance for multibody topology optimization analysis

In order to increase the performance in large-scale topology optimization problems, the researchers typically focus on the numerical algorithms and computational tools. These two approaches have a strong interaction and dependence.

The algorithms can be improved, or new ones can be developed in order to reach a better performance. Some approaches attempt to apply a combination of different algorithms with the aim of increasing the computational efficiency [START_REF] Liao | A triple acceleration method for topology optimization[END_REF]. Other promising strategies reside on the time reduction by solving the nested analysis problem using multigrid methods [START_REF] Amir | On multigrid-CG for efficient topology optimization[END_REF]. A recent survey of the different strategies for speed-up the solution of the topology optimization problems is presented by [START_REF] Mukherjee | Accelerating large-scale topology optimization: State-of-the-art and challenges[END_REF].

In the case of computational tools, the key aspects are the programming models and hardware features. When the algorithms are adequate and well designed, the computational time depends on the programming model and hardware capability. For these reasons, parallel computing opens the gates to increase the computational efficiency for the problem under study.

In the last years, several works have been focused on topology optimization for a single body using parallel computing on CPUs such as [START_REF] Aage | Topology optimization using PETSc: An easy-to-use, fully parallel, open source topology optimization framework[END_REF][START_REF] Aage | Parallel framework for topology optimization using the method of moving asymptotes[END_REF][START_REF] Aage | Topology optimization of large scale stokes flow problems[END_REF][START_REF] Borrvall | Large-scale topology optimization in 3D using parallel computing[END_REF][START_REF] Evgrafov | Large-scale parallel topology optimization using a dual-primal substructuring solver[END_REF][START_REF] Kim | Parallelized structural topology optimization for eigenvalue problems[END_REF][START_REF] Liu | Fully parallel level set method for large-scale structural topology optimization[END_REF][START_REF] Mahdavi | Topology optimization of 2D continua for minimum compliance using parallel computing[END_REF][START_REF] Vemaganti | Parallel methods for optimality criteria-based topology optimization[END_REF]. Moreover, parallel processing on GPUs is another vibrant area of interest [START_REF] Challis | High resolution topology optimization using graphics processing units (GPUs)[END_REF][START_REF] Martínez-Frutos | Efficient topology optimization using GPU computing with multilevel granularity[END_REF][START_REF] Schmidt | A 2589 line topology optimization code written for the graphics card[END_REF][START_REF] Wadbro | Megapixel topology optimization on a graphics processing unit[END_REF][START_REF] Zegard | Toward GPU accelerated topology optimization on unstructured meshes[END_REF]. Among these trends, parallel computing on CPU remains as the dominant one. Regarding to multibody topology optimization for large-scale, there are few works on this topic, an example of this is the work done in the doctoral thesis presented by Moghadasi (2019).

The PETSc library for large-scale topology optimization

In [START_REF] Aage | Topology optimization using PETSc: An easy-to-use, fully parallel, open source topology optimization framework[END_REF] the authors present a public and flexible framework for parallel topology optimization in CPU. The code is written in C++ using object oriented programming, and it is based on PETSc. The default version of this program solves the minimum compliance of cantilever beam problem, and it includes a parallel version of important algorithms such as the method of moving asymptotes, and filters. Besides, in [START_REF] Aage | Giga-voxel computational morphogenesis for structural design[END_REF] the authors report a giga-voxel analysis for morphogenesis in structural desing, where the methodology was implemented using the PETSc library. These examples evidence the suitability of PETSc for large-scale problems in topology optimization.

The use of libraries like PETSc eliminates the need to write from scratch all the low-level math libraries and specialized algorithms for parallel processing. PETSc is a collection of parallelized libraries for linear algebra operations, which includes several types of solvers, methods for partition domain, and it is freely available. Its implementation is in parallel and scalable to thousands of cores, and it is portable to the most popular operating systems.

In order to solve the problem at hand in a reasonable amount of time, a suitable computational platform for multibody topology optimization is developed for parallel processing in this thesis. The platform it is written in C++ programming language, using the object oriented paradigm, it is not software dependent, and it is designed to solve large scale problems. The computational platform is based on PETSc. The next sections are dedicated to provide an overview of the parallel framework for multibody topology optimization.

MUBOTO framework

The code of the computational framework is written based on the object oriented paradigm, therefore the MUBOTO structure relies on the class concept, which are conceived as modular units allowing their operation independently. The mathematical formulation to modeling multibody systems based on topology optimization is described in the Sections 2.2 and 2.3. These algorithms are encapsulate in three classes: Body, Multibody and Optimization.

The layout of MUBOTO framework is presented in Fig. 2.4. The first level represents the most general classes or libraries which do not need modifications when the problem changes (different mechanism). The second level presents those classes which are general and they are independent of the application, but need some modifications when the problem is different. In the upper level there are classes which must be adapted in order to analyze the multibody system under study. In the main file is where the entire process of MUBOTO is managed. which contains sensitivity and density filters. The remaining classes comprising the framework are:

• Body: In this class the geometry, mesh, material properties, joints and the interface nodes are defined. The class contains methods for the computation of the matrix of static modes and also to solve the natural frequency analysis. The proposed methodology is based on the finite element formulation, thus the first step is the discretization of bodies. The mesh of each link is generated using structured grids handled by distributed arrays DMDA provided by PETSc.

The DMDA object manages the parallel communication required while working with data stored in regular arrays.

• Multibody: It contains information related to the multibody model, such as the type of multibody system, number of bodies, and boundary conditions. This class solves the reduced multibody FE models by means of the SUPERLU_DIST solver. In this class the kinematic connections of the multibody system are considered.

• Optimization: This class works as a manager between the Body and Multibody classes. It includes the definition of the objective function, constraints and their gradients. In our case the gradients are computed analytically and they are evaluated instead of being approximated.

The results obtained with MUBOTO can be visualized in ParaView [START_REF] Ahrens | The Visualization Handbook, chapter 36[END_REF], which is a specialized package for display results from large-scale analysis. Besides, the computation of the matrix of static modes (Eq. (2.25)) and the solution of the natural frequency analysis (Eq.

(2.50)) are the more demanding computations in MUBOTO, therefore they are solved using specific solvers, which are described in the following section.

Linear solvers and eigenvalue solvers

Computing the matrix of static modes implies to solve a large sparse linear system, which are the equilibrium equations of the FE model. This computation is accomplished for each body of the mechanical system, generating a high computational effort. In order to solve these linear systems there are two types of solvers: direct solvers (factorization of the stiffness matrix), and iterative solvers. Iterative solvers are the best option for large-scale problems because they are faster than direct solvers, use less memory, and they are well suited for parallelization, since these solvers requires less computations and memory resources [START_REF] Bendsøe | Topology Optimization -Theory, Methods and Applications[END_REF].

Since MUBOTO is based on PETSc, several linear solvers and preconditioners are available throughout the object KSP (Krylov methods). In addition, PETSc interfaces to other specialized high-performance libraries of linear solvers such as HYPRE [START_REF] Falgout | hypre: A library of high performance preconditioners[END_REF], and Su-perLU_DIST [START_REF] Li | SuperLU_DIST: A scalable distributed-memory sparse direct solver for unsymmetric linear systems[END_REF], among others. The static modes computations are performed by a iterative solver due to the size of the problem. Thus, we use the KSPCR solver (provided by PETSc), which implements the conjugate residual method. It is an iterative Krylov subspace numeric method used for solving systems of linear equations, and it is very similar to the popular conjugate gradient method. The preconditioner used is PCJACOBI, which is a diagonal scaling preconditioning. The solution of the reduced linear systems whose size is small (reduced models), is performed by a direct solver provided by SuperLU_DIST, which is a parallel direct solver package for LU factorization.

On the other hand, in order to solve the free vibration problem, we use the Scalable Library for Eigenvalue Problem Computations (SLEPc) [START_REF] Hernandez | SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems[END_REF]. This library is a software for the solution of large-scale sparse eigenvalue problems, it is an extension of PETSc. In our case the Generalized Davidson method is selected to solve the eigenvalue problem. This is an iterative method which can be useful to solve any type of problem (Generalized Hermitian, Generalized Non-Hermitian, etc.), and this method has the possibility to computing different number of eigenvalues.

Summary

In this Chapter the general formulation for multibody topology optimization was described. The methodology is based on the density-based SIMP approach, thereby the finite element method is used in order to model the multibody system. The elastostatic and elastodynamic models of the MBS were formulated based on the SIMP scheme in order to capture the MBS elastic behavior.

In addition, the optimization is performed simultaneously in all the bodies based on the structural response of the entire system. Besides, with the purpose of reduce the computational cost associated to the solution of MBS models, the static and dynamic condensation were exploited. These model reduction techniques leads to reduced matrices which are used to build computationally efficient elastic models of the MBS.

In topology optimization typically appears numerical problems such as mesh-dependence and checkerboard pattern. These issues are solved by means of restriction methods which operate as a filters, in our case it is done by filtering the densities. Regarding the optimization problem resolution, it was concluded that due to the type of problem faced on this thesis, the best option for solving the TO problem is a gradient-based algorithm. Moreover in the view of the optimization problem involves a high number of variables, multiple and different types of constraints, thus the MMA optimizer is the most suitable tool for this situation. With respect to the equality constraints, it is proposed to convert them into inequality constraints in order to avoid a very restrictive optimization scheme. On the other hand, the gradient-based optimizer requires the availability of the gradient of each function involved in the optimization process. These gradients or sensitivities can be computed analytically or numerically, being preferred the analytical expression. Hence, the analytical gradients of the compliance and natural frequency of the MBS were computed.

The mathematical model described in this Chapter was implemented in a computational platform specifically developed for this research work, which is called MUBOTO. Topology optimization for three dimensional multibody systems leads to large models. Therefore, in order to solve the model in a reasonable amount of time, a computational platform for fully parallel processing was developed in C++ programming language. It uses the object oriented paradigm and the distributed memory model for parallel computing, which implies that all algorithms must be written for parallel processing. Moreover, since the platform is developed based on PETSc library, and its complementary package SLEPc, there are available high-performance iterative solvers for linear systems and for eigenvalue problems. As a result, these libraries make possible take the advantage of parallel computing and exploit efficiently the hardware resources.

The computational platform MUBOTO can be used to model the MBS as a flexible systems. Therefore, using this formulation it can be possible to perform structural topology optimization of the whole multibody system while the other constraints are considered.

Chapter 3

Topology optimization for dynamic balancing of four-bar linkage ... la derrota tiene una dignidad que la ruidosa victoria no merece ... Jorge Luis Borges, La Cifra 1981.

H

ISTORICALLY the four-bar linkage has been an object of study for a long time. This linkage is a common component in many mechanical systems, being the simplest 1-DOF linkage. Therefore, the research works which treat different engineering topics on this mechanism are vast. The dynamic balancing is no exception, and in this chapter it is presented as a new approach for the optimum design of a reactionless four-bar linkage relied on multibody topology optimization. The reliability of the proposed approach is validated numerically with commercial software, and a prototype was built for experimental validation of its dynamic balancing performance.

X 0 O 3 O 1 ≡ O 0 O 2 Y 0 S 1 S 2 S 3 O ′ 2 ψ 3 ψ 1 ψ 2 θ 1 B 0 B 1 B 2 B 3 r 1 r 3 r 2 ℓ 1 ℓ 2 ℓ 3 ℓ 4

Introduction

A general scheme of the four-bar linkage is given in Fig. 3.1. This mechanism has 1-DOF, and it is composed of three moving bodies B i , i = 1, . . .3, and one fixed body, B 0 . The mass of body B i is m i , and its length is defined as ℓ i . The center of mass S i of the moving bodies is represented by distances r 1 , r 2 , r 3 and constant angles ψ 1 , ψ 2 , and ψ 3 , which have been defined in a local frame rigidly attached to the respective moving body. Furthermore, θ 1 is the angular position of body B 1 with respect to the X 0 axis.

The dynamic balancing of the four-bar linkage has been the object of study in many researches. These works can be classified as those which perform (i) shaking force balancing [START_REF] Artobolevskii | Methods of inertia calculation for mechanisms of agricultural machines[END_REF][START_REF] Berkof | A new method for completely force balancing simple linkages[END_REF][START_REF] Fisher | Über die reduzierten systeme und die hauptpunkte der glieder eines gelenkmechanismus[END_REF], F.R.E. [START_REF] Crossley | Dynamics in Machines[END_REF][START_REF] Talbourdet | Mathematical solution of 4-bar linkages -IV, balancing of linkages[END_REF][START_REF] Hilpert | Weight balancing of precision mechanical instruments[END_REF][START_REF] Smith | Inertia forces in a four-bar linkage[END_REF], (ii) the complete shaking force and partial shaking moment balancing (Arakelian and Dahan, 2001, Arakelian et al., 2001, Berkof and Lowen, 1971[START_REF] Carson | Feasible parameter design spaces for force and root-mean-square moment balancing an inline 4R 4-bar synthesized for kinematic criteria[END_REF][START_REF] Elliott | The theory of torque shaking force and shaking moment balancing of four link mechanisms[END_REF][START_REF] Haines | shaking moment or driving torque of a force-balanced 4-bar linkage using feasible counterweights[END_REF], Lowen and Berkof, 1971[START_REF] Shchepetilnikov | The determination of the mass centers of mechanisms in connection with the problem of mechanism balancing[END_REF][START_REF] Wiederrich | Momentum balancing of four-bar linkages[END_REF], and the (iii) complete shaking force and shaking moment balancing [START_REF] Arakelian | Complete shaking force and shaking moment balancing of linkages[END_REF][START_REF] Shabana | Dynamics of Multibody Systems[END_REF][START_REF] Bagci | Complete shaking force and shaking moment balancing of link mechanisms using balancing idler loops[END_REF][START_REF] Berestov | Comparative analysis of the reactions in the kinematic pairs of the fourbar linkages for the different balancing methods[END_REF][START_REF] Berkof | Complete force and moment balancing of inline four-bar linkages[END_REF], Briot and Arakelian, 2012[START_REF] Dresig | Vollständiger und harmonischer Ausgleich ebener Mechanismen[END_REF][START_REF] Esat | A theory of complete force and moment balancing of planar linkage mechanisms[END_REF][START_REF] Feng | Complete shaking force and shaking moment balancing of 17 types of eight-bar linkages only with revolute pairs[END_REF][START_REF] Kamenskii | On the question of the balancing of plane linkages[END_REF][START_REF] Ye | Complete balancing of planar linkages by an equivalence method[END_REF], where the complete shaking moment balancing is achieved using additional components. As can be concluded, the shaking force balancing of the four-bar linkages is usually addressed by employing counterweights, and there are well known solutions with the disadvantage of increasing the total mass. However, the shaking moment balancing of the four-bar linkage is more challenging, and it is generally carried out using counter-rotations. Figure 3.2: Design of reactionless robot using a reactionless four-bar linkage as a building leg [START_REF] Wu | Synthesis of reactionless spatial 3-DoF and 6-DoF mechanisms without separate counter-rotations[END_REF]. The reactionless spatial 3-DOF mechanism is used as a building leg in order to design the 6-DOF parallel robot.

On the other hand, to avoid the increase of the moving masses and of the input torques, several works use optimization techniques in order to minimize some unbalanced reaction loads of the four-bar linkage [START_REF] Berkof | Theory of shaking moment optimization of force-balanced four-bar linkage[END_REF][START_REF] Chaudhary | Balancing of four-bar linkages using maximum recursive dynamic algorithm[END_REF][START_REF] Demeulenaere | Ultimate limits for counterweight balancing of crank-rocker four-bar linkages[END_REF][START_REF] Farmani | Multiobjective optimization for force and moment balance of a four-bar linkage using evolutionary algorithms[END_REF][START_REF] Orvañanos-Guerrero | Gradient descent-based optimization method of a four-bar mechanism using fully cartesian coordinates[END_REF][START_REF] Qi | Optimum balancing of four-bar linkages[END_REF][START_REF] Zhang | The optimum balance of shaking force and shaking moment of linkages[END_REF]. Usually these works solves the dynamic balancing problem as an optimization problem, in which the shaking force and moment are used as functions to be minimized or constrained.

In order to avoid the addition of counter-rotations for the shaking moment cancellation, it was shown by Ricard and [START_REF] Gosselin | On the development of reactionless parallel manipulators[END_REF] that four-bar linkages can be fully balanced without the addition of counter-rotations or any auxiliary linkages by forcing a combination of geometric relationships and proper mass distributions. Three different types of reactionless four-bar linkages, characterized by their link lengths have been found.

Additionally, it is possible to exploit this self-balanced linkage as a special module for building reactionless robots. For instance, in [START_REF] Gosselin | Synthesis and design of reactionless three-degree-of-freedom parallel mechanisms[END_REF], the authors present the synthesis of a planar and spatial 3-degrees-of-freedom reactionless mechanism, using the dynamically balanced four-bar linkage without counter-rotations, while in [START_REF] Wu | Synthesis of reactionless spatial 3-DoF and 6-DoF mechanisms without separate counter-rotations[END_REF] this approach is extended to the synthesis of spatial 3-DOF and 6-DOF reactionless mechanism, as can be seen in Fig. 3.2. This strategy represents one of the best alternatives in order to obtain reactionless robots while reducing their design complexity, because it does not require supplementary components for shaking moment balancing.

For robot design purposes, obtaining reactionless four-bar linkage is of interest, but is usually necessary to optimize other performances in parallel. In particular, in [START_REF] Jiang | Dynamic optimization of reactionless four-bar linkages[END_REF], it is shown that the linkage input torque is mainly affected by the mass of the links to be balanced. Therefore, in order to optimize the self-balanced four-bar linkage, and to take into account the elastic behavior of the links, in this research work we propose to use structural topology optimization as a tool for the design of a reactionless four-bar linkage. This approach allows to perform a complete shaking force and moment balancing, because the mass redistribution is constrained by the dynamic balancing conditions included in the optimization problem, and at the same time the elastic behavior of the links is taken into account.

Dynamic balancing conditions of the four-bar linkage

The balancing conditions to achieve a full dynamic balancing of the four-bar linkage without counter-rotations were settled down in the seminal paper of Ricard and [START_REF] Gosselin | On the development of reactionless parallel manipulators[END_REF]. These balancing conditions are based on a set of geometric relations and constraints on the inertial parameters of the links. Indeed, there are three families of this reactionless four-bar linkage, characterized by their links lengths, which are: S 1 : ℓ 1 = ℓ 4 and ℓ 2 = ℓ 3 , S 2 : ℓ 1 = ℓ 3 and ℓ 2 = ℓ 4 , and S 3 : ℓ 1 = ℓ 2 and ℓ 3 = ℓ 4 (Ricard and [START_REF] Gosselin | On the development of reactionless parallel manipulators[END_REF], which lead to the full dynamic balancing without the use of any counter-rotations. These families are schematically represented in Fig. 3.3, where the parameter x i is the location of S i along the straight line between the joints of the body B i (see Fig.

3.1).

The second set of link lengths has proven to be an effective option for design reactionless parallel mechanisms [START_REF] Gosselin | Synthesis and design of reactionless three-degree-of-freedom parallel mechanisms[END_REF]Gosselin, 2004). Therefore, we decided to focus our work on this particular linkage, which is depicted in Fig. 3.3b.

The dynamic balancing conditions, given in (Ricard and [START_REF] Gosselin | On the development of reactionless parallel manipulators[END_REF] for the second set of link lengths, are expressed in terms of bodies inertial parameters (Section 2.5) in order to take into account the elastic model for topology optimization. Following (Briot and Goldsztejn, 2018a), the dynamic balancing conditions are given by: Equations (3.1a) are included in order to simplify the link design and keep its center of mass along the longitudinal link axis. Based on Fig. 3.1, these equations are formulated when ψ 1 = ψ 3 = π and ψ 2 = 0.

I y 1 = 0, I y 2 = 0, I y 3 = 0 (3.1a) I x1 ℓ 1 + m 2 - I x2 ℓ 2 = 0 (3.1b) I x3 ℓ 3 + I x2 ℓ 2 = 0 (3.1c) I zz 1 -I x1 ℓ 1 + I zz 2 -I x2 ℓ 2 = 0 (3.1d) I zz 3 -I x3 ℓ 3 + I zz 2 -I x2 ℓ 2 = 0 (3.1e) O 1 O 2 O ′ 2 O 3 ℓ 1 ℓ 2 ℓ 3 ℓ 4 θ 1 S 1 S 2 S 3 x 1 x 2 x 3 B 1 X Y B 2 B 3 B 0 (a) Set S 1 : ℓ 1 = ℓ 4 and ℓ 2 = ℓ 3 . O 1 O 2 O ′ 2 O 3 ℓ 1 ℓ 2 ℓ 3 ℓ 4 θ 1 S 1 S 2 S 3 x 1 x 2 x 3 B 1 Y X B 3 B 2 B 0 (b) Set S 2 : ℓ 1 = ℓ 3 and ℓ 2 = ℓ 4 . O 1 O 2 O ′ 2 O 3 ℓ 1 ℓ 2 ℓ 3 ℓ 4 θ 1 S 1 S 2 S 3 x 2 x 3 x 1 B 2 B 3 B 0 B 1 X Y (c) Set S 3 : ℓ 1 = ℓ 2 and ℓ 3 = ℓ 4 .
The following sections present a formulation intended to optimize the four-bar linkage taking into account its elastic behavior, while dynamic balancing conditions are satisfied.

Topology optimization for dynamic balancing

In this section we report the optimal design of the reactionless four-bar linkage using multibody topology optimization. In previous sections we described the methodology for linkage modeling, and we presented the dynamic balancing conditions. The following section is dedicated to the definition of the optimization problem and its resolution. All the following analyzes were carried out using the MUBOTO framework running on a computer with a processor AMD Ryzen 9 3900X, using the OS Ubuntu 20.04 64 bits, and PETSc 3.16.0, while the results are visualized with ParaView.

Optimization problem

The four-bar linkage under study was presented in Fig. 3.3b, and the following assumptions are referred to this scheme. The link lengths selected for our studies are: ℓ 1 = ℓ 3 = 60 mm, and ℓ 2 = ℓ 4 = 200 mm, which are proposed in (Briot and Goldsztejn, 2018a). The initial design domain for the movable links is presented in Fig. 3.4. Each link has two joints, which are represented by holes of 4 mm of diameter. Bodies B 1 and B 3 have the same dimensions and thickness of 32 mm, while body B 2 has 10 mm of thickness. The initial dimensions of the links were defined after several simulations, choosing the design/solution with less material between the joints, but with a meaningful design. This implies to have links with bigger thickness. Performing specific analysis in order to select the dimensions of these initial guesses can be a particular topic of study, which is out of the scope of our research work.

The assembled linkage is subject to external loads which are applied, with respect of the global frame, as follows (see Fig. 3.3b): at O 2 and O ′ 2 a force equal to f = [10, 10, 10] T N, applied along the x, y and z global axes. Additionally a moment at O 3 of 1 Nm is applied around z-axis. These loads were defined arbitrarily. Besides, the linkage was optimized with the joint O 1 fixed at θ 1 = π/2.

The approach consists of generating a mesh for each body, and then analyze the bodies as a multibody system under a set of loads and boundary conditions, such as it was described in Chapter 2. The mesh is generated as a structured 3D grid using 8-node linear hexahedral elements. The finite element (FE) analysis is assumed linear elastic, and a penalization factor for the SIMP scheme equal to p = 3. The element size is equal to 2 mm for all bodies and each node has three degree of freedom for Cartesian displacements. Hence, the four-bar linkage is modeled using a total of 116,789 elements and it generates a FE model with 362,700 degrees of freedom.

In order to model a link we use passive and active elements. Both types of elements are considered in the finite element analysis, but only the active ones are included in the optimization process. It is important to mention that the density of passive elements never changes, and their initial state can be zero or one (void or solid).

In our case, we use passive elements with the purpose of modeling the link joints. A void region is used to create the hole in the joint, and the solid region represents the material required to create the joint. This solid region is considered in the computation of inertial parameters.

As usual in topology optimization problems, we apply a filtering procedure in order to obtain a layout without checkerboard problem. This filter modifies the density variables based on the density of their neighborhoods, it is know as density filter and was proposed in [START_REF] Bourdin | Filters in topology optimization[END_REF].

Additionally, we study how materials combination with different densities affects the mechanism footprint. For this purpose, we consider to include in links B 1 and B 3 a cylinder made of steel, while the link is made of material with lower density. Thus, using the cylinder of steel the link footprint will be smaller compared whit a link made of the same material but without the cylinder.

Hence, we performed four optimization processes for different links materials, and they are as follow:

• Case I: All the links are made of nylon.

• Case II: All the links are made of aluminium.

• Case III: All the links are made of nylon and bodies B 1 and B 3 includes a cylinder of steel.

• Case IV: All the links are made of aluminium and bodies B 1 and B 3 includes a cylinder of steel.

For case III and IV, the bodies B 1 and B 3 have a hole to contain a cylinder, as can be seen in All dimensions are in millimeters. Fig. 3.5. The cylinder has diameter of 18 mm and it length is 32 mm, thus its mass is 63.5 g.

The cylinder of steel only affects the inertial parameters and the total mass of the links, and it is not included in the optimization process. Thus, we propose a practical approach where the links are analyzed with a hole instead of multimaterial model. The material properties used for the optimization process are summarized in Table 3.1.

Definition of the optimization problem

In order to prevent negative effects due to excessive deformations, such as wear, low accuracy, and even vibrations, we decide to optimize the linkage stiffness. Stiffness can be defined as a linkage property to sustain loads without excessive deformations, and its reciprocal is know as compliance. Therefore, compliance is established as the objective function, and it is minimized while constant load is applied, under the dynamic balancing constraints of the four-bar linkage. In the following sections the objective function and its constraints are defined as well as their gradients.

Objective function

With the purpose to design a stiff linkage, we decide to minimize the compliance which is twice the strain energy stored in the structure under a given loading. In essence, minimize the compliance is equivalent to maximize the stiffness. Hence, the compliance of the four-bar linkage can be defined by means of the strain energy stored in the linkage (Eq. (2.36)), and it is defined as follows:

f (ρ ρ ρ) = u T r K r u r (3.2)
where the compliance is computed using the linkage stiffness matrix K r defined in Eq. (2.37), and u r is the reduced set of independent coordinates, Eq. (2.33). The gradient of this objective function is computed analytically and it was presented in Section 2.6.2.1.

Constraints

The set of balancing conditions given in Eq. (3.1) are essentially equality constraints, but they are transformed into inequalities as it was described in Section 2.6.1. Hence, the inequality balancing conditions written in terms of the body inertial parameters (Section 2.5) are expressed in the following form:

g 1, 2 (ρ ρ ρ) = ± N 1 ∑ j=1 m 1 j y 1 j ρ 1 j ∓ ε ≤ 0 (3.3a) g 3, 4 (ρ ρ ρ) = ± N 2 ∑ j=1 m 2 j y 2 j ρ 2 j ∓ ε ≤ 0 (3.3b) g 5, 6 (ρ ρ ρ) = ± N 3 ∑ j=1 m 3 j y 3 j ρ 3 j ∓ ε ≤ 0 (3.3c) g 7, 8 (ρ ρ ρ) = ± 1 ℓ 1 N 1 ∑ j=1 m 1 j x 1 j ρ 1 j + N 2 ∑ j=1 m 2 j ρ 2 j - 1 ℓ 2 N 2 ∑ j=1 m 2 j x 2 j ρ 2 j ∓ ε ≤ 0 (3.3d) g 9, 10 (ρ ρ ρ) = ± 1 ℓ 3 N 3 ∑ j=1 m 3 j x 3 j ρ 3 j + 1 ℓ 2 N 2 ∑ j=1 m 2 j x 2 j ρ 2 j ∓ ε ≤ 0 (3.3e) g 11, 12 (ρ ρ ρ) = ± N 1 ∑ j=1 m 1 j x 2 1 j + y 2 1 j ρ 1 j -ℓ 1 N 1 ∑ j=1 m 1 j x 1 j ρ 1 j + N 2 ∑ j=1 m 2 j x 2 2 j + y 2 2 j ρ 2 j -ℓ 2 N 2 ∑ j=1 m 2 j x 2 j ρ 2 j ∓ ε ≤ 0 (3.3f) g 13, 14 (ρ ρ ρ) = ± N 3 ∑ j=1 m 3 j x 2 3 j + y 2 3 j ρ 3 j -ℓ 3 N 3 ∑ j=1 m 3 j x 3 j ρ 3 j + N 2 ∑ j=1 m 2 j x 2 2 j + y 2 2 j ρ 2 j -ℓ 2 N 2 ∑ j=1 m 2 j x 2 j ρ 2 j ∓ ε ≤ 0 (3.3g)
where ε is the relaxation parameter defined in Eq. (2.84), and after generating the inequalities we have twice number of restriction equations. If we consider only the dynamic balancing conditions, then will have fourteen constraints. Besides, the corresponding gradients of the constraint functions are computed assuming a linear relation between the design variables and the inertial properties, thus:

∂ g 1, 2 (ρ ρ ρ) ∂ ρ ρ ρ = ± N 1 ∑ j=1 m 1 j y 1 j (3.4a) ∂ g 3, 4 (ρ ρ ρ) ∂ ρ ρ ρ = ± N 2 ∑ j=1 m 2 j y 2 j (3.4b) ∂ g 5, 6 (ρ ρ ρ) ∂ ρ ρ ρ = ± N 3 ∑ j=1 m 3 j y 3 j (3.4c) ∂ g 7, 8 (ρ ρ ρ) ∂ ρ ρ ρ = ± 1 ℓ 1 N 1 ∑ j=1 m 1 j x 1 j + N 2 ∑ j=1 m 2 j - 1 ℓ 2 N 2 ∑ j=1 m 2 j x 2 j (3.4d) ∂ g 9, 10 (ρ ρ ρ) ∂ ρ ρ ρ = ± 1 ℓ 3 N 3 ∑ j=1 m 3 j x 3 j + 1 ℓ 2 N 2 ∑ j=1 m 2 j x 2 j (3.4e) ∂ g 11, 12 (ρ ρ ρ) ∂ ρ ρ ρ = ± N 1 ∑ j=1 m 1 j x 2 1 j + y 2 1 j -ℓ 1 N 1 ∑ j=1 m 1 j x 1 j + N 2 ∑ j=1 m 2 j x 2 2 j + y 2 2 j -ℓ 2 N 2 ∑ j=1 m 2 j x 2 j (3.4f) ∂ g 13, 14 (ρ ρ ρ) ∂ ρ ρ ρ = ± N 3 ∑ j=1 m 3 j x 2 3 j + y 2 3 j -ℓ 3 N 3 ∑ j=1 m 3 j x 3 j + N 2 ∑ j=1 m 2 j x 2 2 j + y 2 2 j -ℓ 2 N 2 ∑ j=1 m 2 j x 2 j (3.4g)
Consequently, the mathematical formulation of the optimization problem for the optimum design of the reactionless four-bar linkage is expressed as: subject to : g(ρ ρ ρ) ≤ 0,

min ρ ρ ρ∈[0, 1] : f (ρ ρ ρ) = u T r K r u r (3.5)
where g(ρ ρ ρ) is the set of the dynamic balancing constraints given in (3.3). Besides, the objective function and the constraints are normalized using their values computed for the first iteration, excluding Eq. (3.1a) whose initial values are null.

Furthermore, the optimization problem presented in this Chapeter is solved using the gradientbased MMA optimizer [START_REF] Svanberg | The method of moving asymptotes -a new method for structural optimization[END_REF], because it is able to handle large number of variables, as well as multiple constraints, and has excellent convergence properties [START_REF] Bendsøe | Topology Optimization -Theory, Methods and Applications[END_REF]. Since MMA is a first-order optimizer, the analytical gradients of the objective function and constrains are derived in order to improve the computational performance of the optimizer. The dynamic balancing conditions, used as a constraints, are linear with respect to the decision variable, which simplify the gradient computation.

Numerical results

In our case, constraint violation is defined as an index intended to evaluate the fulfillment of the dynamic balancing conditions. Hence, because the optimization problem behavior is monotonic when approaching convergence, the stopping criteria (convergence criterion) is defined based on the number of iterations. Whereby the optimization process was run for each case until it reached 1000 iterations, in such a way that it is ensured a lowest constraints violation. The constraints are normalized using their values for the initial design, thus the constraint violation percentage is computed based on the initial and final values of the constraints.

The optimization results are summarized in Table 3.2, where the computation time, maximum constraint violation and objective function values are displayed for each case. The values of the constraint violation shows that the constraints are properly satisfied. Besides, the objective function evolution is shown in Fig. 3.6 for the four cases. It worth to mention that, because of the optimization process removes material, the strain energy (compliance) increases. Then the optimization algorithm find the optimum material distribution based on the given constraints. This is a typical behavior in structural topology optimization problems.

The next section provides the results after perform the multibody topology optimization for the dynamic balancing of the four-bar linkage. The finite element model, loads, boundary conditions, objective function, and constraints are the same for all study cases, the only difference is the material used in each case.

Results and discussion

The results from our topology optimization procedures are essentially structured meshes, where all the elements are hexahedral of the same size, as shown in Fig. 3.7. These type of results are know as a voxel-based results. Figure 3.7a shows a density field which ranging from 0 to 1, and its corresponding colors are blue and red, respectively. A section view of the link is presented in order to visualize the internal part of the link, which can be seen in Fig. 3.7b. Finally, Fig. 3.7c was generated after applying a filter to remove the elements with density below 0.95.

In the following sections we examine the linkage footprint reduction comparing the four described cases. Then, we describe and analyze the linkage optimized properties focusing on the results of Case I. Numerical validations corresponding to the linkage compliance and dynamic balancing were carried out using commercial specialized software.

Considering the results from Case I, we built a prototype in order to evaluate by an experiment the dynamic balancing of the optimized linkage. Details on the manufacturing process and experi- (c) Elements with a density higher than 0.95. mental setup are described at the end of the section. Besides, the audiovisual material is provided alongside the manuscript.

It is worth mentioning that the links shape obtained for the optimization problems that we define are valid only for those particular problems. While the optimized dynamic balancing properties are valid for the general linkage configuration, the structural properties are optimized for the given linkage configuration.

Linkage footprint reduction

Based on the optimization results, we can consider, for a practical point of view, that the optimized bodies can generate a bulky mechanism. Thus, in order to reduce the footprint of the mechanism, we carried out the optimization process for different material combinations, described as a Case I, Case II, Case III, and Case IV.

We decide to use the volume of the bodies as an index to evaluate the linkage footprint reduction. The volume of each body is computed from the voxel-based results of the optimization process, it means that we evaluate the volume of the body taking into account the elements with a density higher than 0.95. The results of the optimization cases are presented in Fig. 3.8. In these figures we show only body B 1 for all optimization cases, but the volumetric reduction is computed for the entire linkage.

Figure 3.8a shows the optimized body B 1 from Case I. If we analyze the same problem, but including now a cylinder of steel (Case III), then we obtain the body shown in Fig. 3.8b. In addition, in order to compare the differences of using a material with higher density, the overlapping meshes are presented in Fig. 3.8c. The volumetric reduction when Case I and Case III (including the cylinder hole) are compared is 11.21%.

With respect to Fig. 3.8d, the results presented correspond to the optimized body B 1 when the From the described results, we can conclude that the cylinder works as a counterweight, and at the same time, due to its high density less material is required to achieve dynamic balancing in bodies B 1 and B 3 . Indeed, the location of the cylinder should be optimized for an optimum footprint reduction. For our examples, we define this location performing several simulations, choosing the location that produced best results.

Numerical validations

The optimization results that have been presented, can not be treated as a solid object in standard computer-aided design (CAD) software. In order to have an editable CAD file, it is mandatory to perform the body-fitted post-processing, i.e., convert the voxel-based results to a smooth design suitable for a CAD software. Using Paraview we apply some filters in order to generate a triangular mesh which can be exported to STL file. The STL format usually can be read by any standard CAD package, and then the CAD program converts the STL file into its native format. In Fig. 3.9, the conversion process of the optimized body B 1 , from voxel-based results to CAD format is depicted.

Editable CAD files are necessary to make changes in the optimized bodies, as it it shown in Fig. 3.9c, where the link joints where modified. These changes are necessary in order to create an adequate CAD model for the physical prototype and numerical validations.

The optimized properties of the four-bar linkage are the compliance an dynamic balancing. Hence, in order to validate these properties, the compliance is evaluated using ANSYS, and the dynamic balancing of the linkage is analyzed with ADAMS. For this purpose the results from Case I are considered. In Table 3.3 the properties of the optimized reactionless four-bar linkage are presented, as well as the mass, link length and the location of the center of mass given by x i .

Compliance validation

The objective function was defined as the compliance, and it was minimized in order to ensure the stiffest linkage design. With the aim to verify the results obtained with the proposed approach, we perform a numerical validation of the linkage compliance using ANSYS.

The numerical validation with ANSYS is based on the static structural analysis of the linkage. It was carried out using the conditions described in the optimization problem definition. Since the compliance is twice the strain energy, considering the strain energy computed with ANSYS, the compliance can be obtained by adding up the strain energy of all elements. In the optimization process we normalize the objective function, but in order to compare the values, a non-normalized compliance value is used. The compliance value from the optimization procedure is 1.8346 • 10 -3 J, and the ANSYS compliance value is 1.8469 • 10 -3 J, based on these results we obtain a small error of 0.6704%.

Even though the aim of this work is not to perform an intensive linkage structural analysis, it is necessary to verify the structural integrity of the linkage for the described conditions. Thus, the total deformations were analyzed in ANSYS, having a maximum total deformation of 5.37 • 10 -5 m. The total deformation plot is shown in Fig. 3.10, where deformations are scale excessively for display purposes.

Besides, Table 3.4 summarizes the computed joint deformations using ANSYS and the corresponding deformations obtained with the optimization framework. From Table 3.4, and considering the ANSYS results, the maximum joint displacement occurs at the joint O ′ 2 . The linkage has a less stiff behavior when is analyzed in ANSYS, nevertheless the comparative error for the joints deformation is very low, having a value of 2.96 % for joint O 2 , and 1.94 % for joint O ′ 2 . In general, these results agree with the compliance validation, the higher strain energy is due to larger deformations when a constant load is applied.

Differences in the comparative analyses arise form the particularities in each finite element model. In ANSYS, we use the CAD model obtained after file format conversions, and the mesh is generated by using the hex dominant method with the option of low order elements (linear hexahedral), but due to the complex shape of the links, in some regions tetrahedral elements appears. In our optimization program we use structured mesh with hexahedral elements. Thus, we do not expect to have equal values, but very similar ones as they are shown.

Validation of dynamic balancing

Numerical validation of the dynamic balancing was carryout using the ADAMS software. In order to analyze the results, the optimized model is compared with an unbalanced four-bar linkage. This unbalanced linkage has been modeled with the same links length and mass as the optimized links, but with the center of mass placed between joints for each link.

The reactions loads on the base, for the balanced and unbalanced linkages, have been simulated for arbitrary motions. The input angular velocity in the body B 3 is a sinusoidal function equal to 0.37 sin(πt), and the gravity is not considered. In Fig. 3.11 the corresponding results of balanced and unbalanced models are presented. Figure 3.11a shows the corresponding reaction force of balanced and unbalanced linkages. With respect to the shaking moments, the reaction moment for the dynamic balanced linkage and unbalanced one is shown in Fig. 3.11b. The computed reaction forces and moments clearly demonstrate the dynamic balancing of the four-bar linkage.

In addition, in Table 3.5 we present a comparison between the balanced and unbalanced linkages, taking into account the RMS of the reaction loads. The format conversion from structured mesh to CAD format produces some numerical errors, which generate this very small unbalance in the linkage. Using a finer mesh can be reduced these numerical errors, but finer mesh leads to a huge computational cost, which is not justified in our case due to the small error.

Experimental validation

Because of the complex shapes generated with topology optimization, the best option to fabricate the optimized bodies is the Additive Manufacturing (AM) technique. Nevertheless, there are some situations where the optimization parameters generates some particular shapes, which can be fabricated with traditional manufacturing techniques.

In our optimization problems, the applied loads affect considerably the resultant shapes. If the applied load in Case I is replaced with a force equal to f = [1, 1, 1] T N (ten times smaller than the original), some cavities begin to emerge in the optimized bodies, what makes impossible to fabricate the links with traditional manufacturing techniques. In Fig. 3.12, the optimization results based on this modified Case I are presented, and as we can appreciate in Fig. 3.12b, the body has cavities. For the load conditions that we define in our optimization problems, the resultant shapes do not present internal hollows. Thus, the prototype presented in this work was fabricated with a CNC milling machine, using the optimal solution of the Case I. The CAD model and the physical prototype are shown in Fig. 3.13. For the physical realization, the body B 2 was divided in two parts (not for numerical validations) in order to allow the motion without interference with the other bodies, and at the same time the bodies disposition keeps the center of mass of all the bodies in the same plane.

Experimental setup

In order to verify the dynamic balancing of the optimized linkage, experimental tests were carried out. The purpose is to evaluate the effects of the residual shaking force and moment in the case of the dynamically balanced linkage. The usual experimental setup for the performance evaluation of a dynamically balanced mechanisms includes force/torque sensors [START_REF] Laliberté | Synthesis, optimization and experimental validation of reactionless two-dof parallel mechanisms using counter-mechanisms[END_REF][START_REF] Martini | Elastodynamic effects of mass-balancing: experimental investigation of a four-bar linkage[END_REF][START_REF] Van Der Wijk | Design and experimental evaluation of a dynamically balanced redundant planar 4-rrr parallel manipulator[END_REF][START_REF] Zomerdijk | Structural design and experiments of a dynamically balanced inverted four-bar linkage as manipulator arm for high acceleration applications[END_REF], moreover the prototype and the remaining electronic components are industrial grade devices, which is totally reasonable and they are the best option to conduct this type of experimental evaluations. Nevertheless, we are proposing an affordable and simple alternative technique to evaluate the dynamic balancing of the optimized four-bar linkage by means of video-processing.

The experimental setup consists in a fixed structure, the linkage base and the cables. Hence, the base of the linkage was suspended from the fixed structure with four cables, as it can be seen in Fig. 3.14. This arrangement allows the linkage to move freely in the presence of unbalanced reaction loads. We use a direct current motor to provides oscillating movement in the linkage, which is attached to body B 3 . With the aim of studying the prototype performance, we record with a camera the location of its suspended based during its motion, and the resultant oscillatory motion is analyzed using video processing 1 . Indeed, the amplitude of the oscillatory motion of the platform is directly related to the amplitude of the shaking force and shaking moment exerted by the moving linkage: the smaller the motion, the smaller the unbalance effects. For this task, we use the open source software called Tracker [START_REF] Brown | Video modeling: Combining dynamic model simulations with traditional video analysis[END_REF], which is an image and video analysis package. Besides, in order to compute the linkage base displacements, it is fundamental to define an inertial reference frame, calibrate the video scale, and designate the feature to be tracked. The video scale is the ratio of a real dimension to an image dimension in pixels between two points, and it is computed providing the real dimension of an object in the video.

As can been seen in Fig. 3.14, the linkage is recorded from the top-view. The home position is defined by the reference frame XY shown in Fig. 3.15, and it remains fixed while the mechanism is in motion. We define the bearing on the joint O 1 as the feature of interest, this feature will be tracked in the video during the motion of the linkage. Therefore, using the position of the tracked feature based on the reference frame, the linkage translations were determined. Besides, the platform rotation is computed around the Z-axis, using the translation of two features belonging to the linkage.

The experiments were performed for two scenarios, as we can see in Fig. 3.15. In the first scenario, we use the optimized linkage to observe its behavior, Fig. 3.15a. Then, for the second scenario we attach an extra mass (46.5 g) on the body B 2 in order to induce an unbalance. This mass is placed near to the joint O 2 , see Fig. 3.15b. For the two scenarios we use the same input motion equals to cos(1.6πt), which is generated by the DC motor attached to body B 3 .

Prototype performance

In Fig. 3.16 we present a comparative evolution of the computed displacements for the first and second scenarios. In Fig. 3.16a, the magnitude of the displacements for the balanced linkage are presented, and the corresponding displacements for the unbalanced case are given in Fig. 3.16b. The platform rotation is shown in Fig. 3.16c for the balanced linkage, and the Fig. 3.16d shows the respective platform rotation for the unbalanced linkage.

From the experimental results, we can observe small translations and rotations in the case of the balanced linkage, which implies that there are residual unbalanced reaction loads. These unbalanced reaction loads can be the consequence of numerical errors, generated when the file format conversions were made, and of course due to the errors generated for the manufacturing process, and because of the quality of the mechanical components. Nonetheless, if we compare the first and second scenarios, the transnational displacement reduction based on their RMS values is 85.48%, having a maximum translation of 1.16 • 10 -3 m for balanced case, and 6.77 • 10 -3 m for unbalanced one. Regarding to the platform rotation, we can find a maximum absolute value of 0.42 deg when the mechanism is dynamically balanced, and a maximum rotation of 3.36 deg for the unbalanced case, these results lead to a reduction of 89.26 % with respect of their RMS values. All these results show that, by using TO, we were able to design a reactionless four-bar linkage while its elastic behavior is taken into account during the design process. Note that a movie showing the frame's motion for the balanced and unbalanced cases is provided in the attached multimedia content.

Summary

The dynamic balancing of the four-bar linkage was addressed by means of topology optimization. Through this approach the mass of the links is redistributed, based on a optimization procedure, in order to fulfill the dynamic balancing conditions (without counter-rotations). The proposed methodology makes possible to model the four-bar linkage as flexible MBS in order to take into account the elastic behavior of the linkage. Even more, the linkage optimization can be performed based on the entire linkage response and considering structural performance indices. Thus, dynamic balancing conditions can be fulfilled while the link elasticity is considered. The proposed methodology was implemented in the computational framework MUBOTO, which has been specially designed for the large-scale multibody topology optimization analysis.

The dynamically balanced four-bar linkages was numerically validated using ADAMS software. This optimized linkage evinced an excellent dynamic balancing performance since the unbalanced reaction loads were almost canceled. The residual shaking forces and moments which appear in the ADAMS simulation were generated because of the body-fitted post-processing, but this unbalance was almost negligible. Besides, the optimized four-bar linkage was compared with an unbalanced linkage. The comparison shown a shaking force and moment reduction of 97.26% and 98.98%, respectively.

Even if our main objective was not to perform an intensive structural design, in order to take into account the link elasticity, we optimized the linkage stiffness through minimizing the compliance. In this way, we avoided excessive deformations, which are source of low accuracy, wear and even vibrations. The linkage compliance was validated with the commercial software ANSYS. The comparison between the optimized compliance value and the compliance obtained with ANSYS by structural analysis has an error of 0.6704%. This small error indicated that the optimization was successfully performed.

Furthermore, with the aim to reduce the linkage footprint, it was proposed a practical approach.

Using a cylinder of steel in the bodies B 1 and B 3 it acts as a counterweight and because of its high density less material is required in the links. This method was applied for two different material combinations, obtaining the maximum linkage footprint reduction of 11.21%.

On the other hand, it was proposed to use video processing as an affordable technique for experimental dynamic balancing evaluation. For the experimental setup the linkage was mounted on a platform which is hanged by four cables attached to a fixed base. In this way free motion of the platform is allowed when an unbalanced load is generated during the linkage operation (driven by a DC motor). A prototype was built, and its dynamic balancing performance was analyzed by means of this experimental approach using video processing, wherewith the effects of the residual shaking force and moment were evaluated. These effects are embodied as translations and rotations on the base of the optimized linkage, which were computed and measured with an unbalanced linkage. The experimental displacements of the optimized linkage, compared with the unbalanced one, are clearly smaller. Hence, comparing the balanced and unbalanced scenarios, the translational displacement reduction based on their RMS values is 85.48%, while the reduction in the platform rotation is 89.26%. This simple, but enlightening experiment, allowed to evaluate the dynamic balance of the optimized linkage beyond a numerical validation.

B C(x, y) D E θ 2 θ 3 Y 0 X 0 θ 4 θ 1 O ℓ 1 ℓ 2 ℓ 3 ℓ 4 d

Introduction

The five-bar mechanism is a two degree of freedom parallel robotic mechanism, composed of five revolute joints with parallel axes. Its schematic representation is shown in Fig. 4.1. This mechanism is able to position its end-effector, located at point C with coordinates (x, y), in a plane. The joints located in A and E are usually the active ones, and they are in gray color to represent the actuated joints.

It is usual that the bodies between joints A -B, and E -D are denoted as proximal links, and the bodies between joints B -C and D -C are called distal links. Besides, the five-bar linkage can be analyzed as a mechanism embodied by two "legs", being the left leg the kinematic chain formed by the bodies between the joints A -B -C, and the right leg that formed by joints E -D -C.

Dynamic balancing conditions of five-bar mechanism

The partial and complete dynamic balancing of the five-bar linkage has been studied by several researchers such as [START_REF] Alici | Optimum dynamic balancing of planar parallel manipulators based on sensitivity analysis[END_REF], [START_REF] Ilia | A novel formulation of the dynamic balancing of five-bar linkages with applications to link optimization[END_REF], [START_REF] Jean | Static balancing of planar parallel manipulators[END_REF], [START_REF] Lecours | Reactionless two-degree-of-freedom planar parallel mechanism with variable payload[END_REF]Gosselin (2010), van der Wijk (2014) and [START_REF] Acevedo | Conditions for dynamic balancing of planar parallel manipualtors using natural coordinates and their application[END_REF], among others. All these papers agree that for this mechanism, the shaking force balancing can be fully accomplished by relocating the center of mass of the links. Nevertheless, the complete shaking moment balancing is only possible with the addition of external components.

A B C D E Y 0 X 0 O S 1 S 2 S 3 S 4 r 1 r 2 r 3 r 4 ψ 1 ψ 2 ψ 3 ψ 4 B 2 B 1 B 0 B 3 B 4

Shaking force balancing conditions

For the sake of completeness, in Appendix A are derived the shaking force balancing conditions. These conditions allow to cancel the unbalanced forces at the base of the five-bar mechanism, keeping the linkage center of mass stationary. Hence, based on the mechanism parametrization shown in Fig. 4.2, the general shaking force balancing conditions are written as follows:

m 1 r 1 cos ψ 1 + ℓ 1 ℓ 3 m 3 r 3 cos ψ 3 + m 2 ℓ 1 = 0 (4.1a) m 1 r 1 sin ψ 1 + ℓ 1 ℓ 3 m 3 r 3 sin ψ 3 = 0 (4.1b) m 2 r 2 cos ψ 2 + ℓ 2 ℓ 3 m 3 r 3 cos ψ 3 = 0 (4.1c) m 2 r 2 sin ψ 2 + ℓ 2 ℓ 3 m 3 r 3 sin ψ 3 = 0 (4.1d) m 4 r 4 cos ψ 4 - ℓ 4 ℓ 3 m 3 r 3 cos ψ 3 + m 3 ℓ 4 = 0 (4.1e) m 4 r 4 sin ψ 4 - ℓ 4 ℓ 3 m 3 r 3 sin ψ 3 = 0 (4.1f)
where m i and ℓ i are the total mass and length of body B i , respectively. Besides, the location of the center of mass S i of body B i , is described by the distance r i and the constant angle ψ i .

For practical purposes, the balancing conditions given in Eq. (4.1) can be simplified. Thus, the following considerations are assumed:

• angles ψ 1 = π and ψ 4 = π,

• the center of mass S 2 is located on B, and S 3 on D, which implies that r 2 = r 3 = 0,

• all links are of equal lengths:

ℓ 1 = ℓ 2 = ℓ 3 = ℓ 4 .
Therefore, the simplified shaking force balancing conditions are given by: Considering links of equal lengths makes possible to impose that the leg A -B -C is symmetrical to the leg E -D -C. It means that, both proximal links and both distal links, should have the same shape and the optimization problem can be solved using symmetry, i.e., the mechanism is analyzed for only one leg (but the whole system is considered). When the problem is analyzed using symmetry, the optimization problem is solved only for one leg of the mechanism, but the structural response takes into account the whole system. It is possible because the symmetrical leg is obtained by reflecting the optimized leg, and thus the results involve the entire linkage. Thus, by applying conditions of symmetry is possible to reduce the computational time in the optimization process because of the considered design variables are the half of the total.

-m 1 r 1 + m 2 ℓ 1 = 0 (4.
In addition, defining the same length for all links makes possible to avoid holes on the five-bar mechanism workspace [START_REF] Figielski | Towards development of a 2-DOF planar parallel robot with optimal workspace use[END_REF]. Furthermore, this condition allows a symmetrical distribution of the mechanism performance in the dextrous workspace1 with respect to the Y 0 -axis [START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF]. After these considerations, the resultant shaking force balancing conditions described in terms of their inertial parameters, are given by the following expressions:

I y 1 = 0, I y 4 = 0 (4.3a)
I y 2 = 0, I y 3 = 0 (4.3b)

I x2 = 0, I x3 = 0 (4.3c) -I x1 + m 2 ℓ 1 = 0, -I x4 + m 3 ℓ 4 = 0 (4.3d)
where Eq. (4.3a) and (4.3b) are necessary in order to constraint the center of mass of each link along the longitudinal axis. Using these balancing conditions it will be generated a five-bar mechanism with the characteristics shown in Fig. 4.3.

A B, S 2 C(x, y) D, S 3 E Y 0 X 0 θ 4 θ 1 O ℓ 1 ℓ 2 ℓ 3 ℓ 4 S 1 S 4

Inertial properties of the five-bar mechanism

In this section the angular momentum of the linakge is studied in order to identify the inertial parameters which produce a lower shaking moment at the base of the linkage. In Appendix A, the general expression of the angular momentum was derived, and it is replicated here:

h A z = I S 1 + m 1 r 2 1 + m 2 ℓ 2 1 + m 2 ℓ 1 r 2 cos(θ 1 -θ 2 -ψ 2 ) θ1 + I S 2 + m 2 r 2 2 + m 2 ℓ 1 r 2 cos(θ 1 -θ 2 -ψ 2 ) θ2 + I S 3 + m 3 r 2 3 + m 3 ℓ 4 r 3 cos(θ 3 -θ 4 -ψ 3 ) + m 3 dr 3 cos(θ 3 -ψ 3 ) θ3 + I S 4 + m 4 r 2 4 + m 3 ℓ 2 4 + m 3 ℓ 4 r 3 cos(θ 3 -θ 4 -ψ 3 ) + m 4 dr 4 cos(θ 4 -ψ 4 ) + m 3 dℓ 4 cos(θ 4 ) θ4 (4.4)
where the parameters m i , and ℓ i are the mass and length of the link B i , respectively. The constant angle ψ i and the distance r i describes the location of the COM of the body B i , while the angle θ i defines the orientation of the i-th body with respect to an horizontal axis, and θi is its timederivative. I S i is the moment of inertia of the moving body B i around an axis perpendicular to the plane of motion.

Analyzing Eq. ( 4.4), it can be noticed that the angular momentum of the five-bar linkage is not constant, since it depends on the angular velocities of the bodies. Therefore, in order to eliminate the shaking moment it is mandatory to include external components. Nevertheless, instead of eliminating the shaking moment by external components, the constant factors can be constrained to be lower than a given threshold, with the aim to reduce the angular momentum effects.

Considering a fully force balanced linkage, the angular momentum given in Eq. (4.4) can be simplified applying the shaking force balancing conditions, namely ψ 1 = ψ 4 = π, r 2 = r 3 = 0, and equal links lengths. In consequence, the angular momentum of a force balanced linkage is defined as follows:

h A z = I S 1 + m 1 r 2 1 + m 2 ℓ 2 1 θ1 + I S 2 + m 2 r 2 2 θ2 + I S 3 + m 3 r 2 3 θ3 + I S 4 + m 4 r 2 4 + m 3 ℓ 2 4 θ4 (4.5)
which revels the well known fact: the angular momentum of a force balanced linkage is independent of the reference point [START_REF] Berkof | Theory of shaking moment optimization of force-balanced four-bar linkage[END_REF].

Because in our analyses we use symmetry, the angular momentum is considered as:

h A z =2 I S 1 + m 1 r 2 1 + m 2 ℓ 2 1 θ1 + 2 I S 2 + m 2 r 2 2 θ2 (4.6)
The angular momentum given in Eq. (4.6), can be expressed in terms of the moment of inertia of the bodies with respect to their axis of rotation, i.e., I zz1 = I S 1 + m 1 r 2 1 and I zz2 = I S 2 + m 2 r 2 2 , thereby it is given as follows:

h A z =2 I zz1 + m 2 ℓ 2 1 θ1 + 2I zz2 θ2 (4.7)
From Eq. (4.7) can be visualized that, in order to reduce the shaking moment without use external components, the best approach is to constrain or minimize the moment of inertia of body B 1 , and also to constrain or minimize the mass and moment of inertia of body B 2 .

In the following section the optimization procedure for the partial dynamic balancing is described. The shaking force balancing conditions give in Eq. ( 4.3) are used as a constraints in the optimization problem, along with the corresponding constraints on the moment of inertia of body B 1 , and the mass and moment of inertia of body B 2 .

Topology optimization for partial dynamic balancing

In this section the optimization problem for partial dynamic balancing of the five-bar mechanism is defined. The optimization solver used in this work requires the computation of the gradients of all the involved functions, and they are computed using an analytical approach.

The numerical solution of the optimization problem was carried out using the computational platform MUBOTO. All the following analyzes were performed on a computer with a processor AMD Ryzen 9 3900X, using the OS Ubuntu 20.04 64 bits, and PETSc 3.16.0, while the results are visualized with ParaView.

Optimization problem

The link lengths defined for our studies are ℓ i = 90 mm, with i = 1 . . . The problem formulation includes external loads applied to the joint C. These loads are equal to f = [10, 10, 10] T N, applied along the x, y and z global axes, in addition a moment of 0.1 Nm is applied at the same point around z-axis. The end-effector of the linkage is placed in the center of the regular dextrous workspace. This workspace is defined as a rectangle centered in (0, 110) mm of width equal to 120 mm and height equal to 40 mm. The material of the bodies is defined as steel, with a density of 7800 kg/m 3 , Young's modulus of 2.1 • 10 11 Pa and a Poisson ration of 0.3. The bodies are meshed with a structured 3D grid using 8-node linear hexahedral elements. The finite element model is considered as linear elastic model, and the penalization factor for the SIMP scheme is equal to p = 3. The element size is equal to 2.5 mm for all bodies and each node has three degree of freedom for Cartesian displacement. In order to modeling the bodies, passive and active elements are used. The passive elements are used to model the link joints. Moreover, since the analysis relies on the symmetry approach, the size of the optimization problem is reduced in a half. Therefore the five-bar mechanism is modeled using a total of 162,560 elements and it generates a FE model with 532,926 degrees of freedom. Besides, in order to avoid the checkerboard and mesh-dependence problems, the density filter was applied.

The shaking force balancing of the five-bar linkage is accomplished keeping stationary the center of mass of the entire mechanism. Thus, because the shaking force balancing conditions are derived for generic mechanism configuration, the optimization solution for one configuration guarantees a shaking force balancing for any mechanism configuration, unlike the shaking moment which depends on the joint variables. Indeed, the structural performance can be different depending on the mechanism configuration.

For multibody systems subjected to dynamic loads, frequency optimization is of great importance in order to avoid severe vibration. This excessive vibration occurs when the frequency of the dynamic loading is close to the natural frequencies of the MBS, therefore is important to maximize the first natural frequency. In this case study the objective function is defined as the first natural frequency and it is maximized. Additionally, in order to perform partial dynamic balancing, the optimization problem is constrained by: (a) shaking force balancing conditions, (b) mechanism inertial properties, and (c) mechanism compliance. The latter ensures the mechanism capability to support loads adequately. In the following sections the objective function and the constraints are defined as well as their gradients.

Objective function

The first natural frequency or fundamental frequency is associated with the highest level of energy due to vibrations and at this frequency appears the displacements with the highest amplitude [START_REF] Briot | On the optimal design of parallel robots taking into account their deformations and natural frequencies[END_REF]. Maximizing the fundamental frequency allows a better control of the mechanism vibration response. Consequently, the objective function, defined as the fundamental fre-quency of the mechanism, is given as follows:

f (ρ ρ ρ) = f 1 = ω 1 2π (4.8)
where ω 1 is the first natural pulsation obtained from the solution of Eq. (2.73): ω 2 1 M d -K d u d 1 = 0. Moreover, the natural frequency and the pulsation are related by ω 1 = 2π f 1 . The gradient of this objective function is computed analytically and it is presented in Section 2.6.2.2. Topology optimization for frequency optimization is more computationally expensive than compliance analysis and it is especially challenging to apply to multibody topology optimization problems. Therefore, the methodology takes advantage of dynamic condensation as was described in Section 2.3.2. In our problem formulation we define s = 6, therefore the first six eigenmodes associated with the Eq. (2.51) are computed.

Constraints

The optimization problem is constrained by three different types of functions. These constraints and their gradients are derived in the following sections.

Shaking force balancing conditions

The shaking force balancing conditions are given in Eq. (4.3), nevertheless because of the use of symmetry, only the set of balancing conditions corresponding to body B 1 and B 2 are considered.

These equations are equality constraints, but they can be written as inequality constrains without loss of generality as it was reported in Section 2.6.1. Writing the balancing conditions in terms of the body inertial parameters (Section 2.5), and converting the equality conditions into inequalities, the following constraints are obtained:

g 1 (ρ ρ ρ) = N 1 ∑ j=1 y 1 j m 1 j ρ 1 j -ε ≤ 0 (4.9a) g 2 (ρ ρ ρ) = - N 1 ∑ j=1 y 1 j m 1 j ρ 1 j + ε ≤ 0 (4.9b) g 3 (ρ ρ ρ) = N 2 ∑ j=1 y 2 j m 2 j ρ 2 j -ε ≤ 0 (4.9c) g 4 (ρ ρ ρ) = - N 2 ∑ j=1 y 2 j m 2 j ρ 2 j + ε ≤ 0 (4.9d) g 5 (ρ ρ ρ) = N 2 ∑ j=1
x 2 j m 2 j ρ 2 j -ε ≤ 0 (4.9e)

g 6 (ρ ρ ρ) = - N 2 ∑ j=1
x 2 j m 2 j ρ 2 j + ε ≤ 0 (4.9f)

g 7 (ρ ρ ρ) = N 1 ∑ j=1 x 1 j m 1 j ρ 1 j + ℓ 1 N 2 ∑ j=1 m 2 j ρ 2 j -ε ≤ 0 (4.9g) g 8 (ρ ρ ρ) = - N 1 ∑ j=1 x 1 j m 1 j ρ 1 j + ℓ 1 N 2 ∑ j=1 m 2 j ρ 2 j + ε ≤ 0 (4.9h)
where ε is the relaxation parameter given in Eq. (2.84).

Considering a linear relation between the decision variables and the inertia parameters, the gradients of these constraints are computed as follows:

∂ g 1 (ρ ρ ρ) ∂ ρ ρ ρ = N 1 ∑ j=1 y 1 j m 1 j (4.10a) ∂ g 2 (ρ ρ ρ) ∂ ρ ρ ρ = - N 1 ∑ j=1 y 1 j m 1 j (4.10b) ∂ g 3 (ρ ρ ρ) ∂ ρ ρ ρ = N 2 ∑ j=1 y 2 j m 2 j (4.10c) ∂ g 4 (ρ ρ ρ) ∂ ρ ρ ρ = - N 2 ∑ j=1 y 2 j m 2 j (4.10d) ∂ g 5 (ρ ρ ρ) ∂ ρ ρ ρ = N 2 ∑ j=1
x 2 j m 2 j (4.10e)

∂ g 6 (ρ ρ ρ) ∂ ρ ρ ρ = - N 2 ∑ j=1 x 2 j m 2 j (4.10f) ∂ g 7 (ρ ρ ρ) ∂ ρ ρ ρ = N 1 ∑ j=1 x 1 j m 1 j + ℓ 1 N 2 ∑ j=1 m 2 j (4.10g) ∂ g 8 (ρ ρ ρ) ∂ ρ ρ ρ = - N 1 ∑ j=1 x 1 j m 1 j + ℓ 1 N 2 ∑ j=1 m 2 j (4.10h)
In this case, the analytical gradient computation is straightforward due to the linear relation between the design variables and the functions.

Compliance

Compliance is twice the strain energy stored in a mechanical system when it is subjected to external loads, and minimizing the compliance is equivalent to maximizing the stiffness of the system. In order to prevent excessive deformations in the mechanism, the compliance is defined as a constraint function, and it is obtained from the strain energy of the MBS given in Eq. (2.36). Therefore, the constraint function based on the mechanism compliance must be lower than a given threshold, and it is expressed as:

g 9 (ρ ρ ρ) = u T r K r u r -c 0 ≤ 0 (4.11)
where mechanism stiffness matrix K r is given in Eq. (2.37), and the vector of reduced set of independent coordinates u r is defined in Eq. (2.33). The constant parameter c 0 is computed as a fraction of the compliance obtained at the first iteration, and its values is c 0 = 1.1 in our case study. It should be noted that because the optimization process removes material, the strain energy (compliance) increases.

Mechanism inertial properties

In Section 4.2.2 the angular momentum of the five-bar mechanism was studied. There it is concluded that in order to reduce the shaking moment, the constant terms of the angular momentum (I zz1 + m 2 ℓ 2 1 ) and (I zz2 ) should be constrained to a given threshold. Considering the first constant term of the angular momentum (I zz1 + m 2 ℓ 2 1 ), it can be written in the following form:

I zz1 + m 2 ℓ 2 1 -I th1 ≤ 0 (4.12)
where I th1 is the threshold based on the percentage of the initial function value. In our case study this threshold is defined as I th1 = 0.15.

On the other hand, the body B 2 is divided in two parts, the right-side and left-side, taking as a reference the vertical line which passing by the interior link joint in Fig. 4.4b (joint B in Fig. 4.1). Hence, the second constant term of the angular momentum, I zz2 , can be written in the following form:

I zz2 = I zz2 L + I zz2 R (4.13)
with I zz2 L as the moment of inertial generated by the corresponding elements on the left-side, while the I zz2 R is the moment of inertial of the right-side. Because of I zz2 R is generated for the elements in the right-side of the link, which is the side that supports the loads, thus the most simple approach is to define the following condition:

I zz2 L = I zz2 R (4.14)
where it should be considered that the left-side of the links is used for balancing purposes. Therefore, the constraints are given as follows:

I zz1 + m 2 ℓ 2 1 -I th1 ≤ 0, (4.15) I zz2 L -I zz2 R = 0. (4.16)
Besides, because of Eq. (4.16) is a equality constraint, this must be written as inequality. Therefore the constraints related to the mechanism inertial parameters can be written in the following form:

g 10 (ρ ρ ρ) = N 1 ∑ j=1 m 1 j x 2 1 j + y 2 1 j ρ 1 j + ℓ 2 1 N 2 ∑ j=1 m 2 j ρ 2 j -I th1 ≤ 0 (4.17a) g 11 (ρ ρ ρ) = N 2 ∑ j=1 m 2 j x 2 2 j + y 2 2 j ρ 2 j L - N 2 ∑ j=1 m 2 j x 2 2 j + y 2 2 j ρ 2 j R + ε ≤ 0 (4.17b) g 12 (ρ ρ ρ) = - N 2 ∑ j=1 m 2 j x 2 2 j + y 2 2 j ρ 2 j L + N 2 ∑ j=1 m 2 j x 2 2 j + y 2 2 j ρ 2 j R -ε ≤ 0 (4.17c)
where ε is the relaxation parameter given in Eq. (2.84). Consequently, the gradients of these constraints are given as follows:

∂ g 10 (ρ ρ ρ) ∂ ρ ρ ρ = N 1 ∑ j=1 m 1 j x 2 1 j + y 2 1 j + ℓ 2 1 N 2 ∑ j=1 m 2 j (4.18a) ∂ g 11 (ρ ρ ρ) ∂ ρ ρ ρ = N 2 ∑ j=1 m 2 j x 2 2 j + y 2 2 j , ∂ g 12 (ρ ρ ρ) ∂ ρ ρ ρ = - N 2 ∑ j=1 m 2 j x 2 2 j + y 2 2 j (4.18b)
notice that the last two equations are written for the whole body B 2 .

Definition of the optimization problem

The optimization problem is defined as maximizing the first natural frequency of the five-bar mechanism, constrained by the shaking force balancing conditions, mechanism compliance and mecha- g i (ρ ρ ρ) ≤ 0, i = 10 . . . 12. Mechanism inertial parameters (4.19d)

Solving the optimization problem in the manner that was defined in Eq. (4.19) produces results that show a high number of elements with intermediate density, which is a undesirable behavior. Figure 4.5 shows an example of this issue, which occurs mainly in the left-side of the body.

Therefore, it is necessary to apply partial penalization to some constraints in order to reduce the elements with intermediate density (Briot and Goldsztejn, 2018a). Partial penalization consists in divided the link in two parts, the left-side and the right-side, where the division is done in the first joint, such as it is shown in Fig. 4.6. Then, the densities of the elements on the left-side are penalized in order to force a binary density field (elements with 0 or 1 density values). This procedure is done for those bodies which present high number of elements with intermediate density, in our case it happens in the left-side of all the links. Hence, the constraints which are penalized are the shaking force balancing conditions (Eqs. 4.9) and the inertial mechanism parameters (Eqs. 4.17). These functions are linear with respect to the decision variable, therefore taking into account all the bodies, these constraints can be written in the following matrix form:

Pρ ρ ρ = P L ρ ρ ρ L + P R ρ ρ ρ R = 0 (4.20)
where Pρ ρ ρ is the expression of the shaking force balancing conditions and the inertial mechanism parameters under a matrix form. Besides ρ ρ ρ is the full design vector. In consequence, P L is a sub-matrix of P, which corresponds to the design variables in the left-side of the links, thereby ρ ρ ρ L contains the design variables corresponding to the elements on the left-side of the links. Regarding to P R and ρ ρ ρ R , they are arranged in the same way already described, corresponding to the right-side of the links.

With the aim to reduce the elements with intermediary density values in the left-side of the links, it is necessary to force the values of ρ ρ ρ L to be only 0 or 1, thus the densities in the left-side of Eq. (4.20) are penalized as follows:

P L ρ ρ ρ q L + P R ρ ρ ρ R = 0 (4.21)
where q is the partial penalization parameter, which should be q ≥ 1. In our case study, the value of the partial penalization parameter is defined as q = 2 in order to obtain binary results.

Numerical results

The convergence criterion was defined based on the change of the objective function. The optimization is performed while the computed value in the k-th iteration is greater than the previous one. Thus, while f (ρ ρ ρ k ) > f (ρ ρ ρ k-1 ) the optimization is running, otherwise the algorithm stops. The optimization process achieve an optimum solution in 94 iterations and it takes 8.78 h. The value of the objective function was 681.28 Hz, and the maximum constraint violation was 1 • 10 -3 , which happens for the constraint (4.3d). The evolution of the objective function is shown in Fig. 4.7.

In the following section, the numerical results are transformed into structured mesh which are post-processed in order to generate a proper CAD solids. Using the post-processed bodies, the optimized properties of the linkage are validated using commercial software.

Results and discussion

The results from topology optimization problems are post-processed in order to generate a editable CAD file. Then, the reactive loads at the base of the optimized five-bar mechanism are verified using ADAMS software, while the natural frequency, compliance and joint displacements are analyzed through ANSYS.

Post-processing

In Fig. 4.8 the voxel-based results of the optimized bodies B 1 and B 2 are presented. From these results, one may observe that the partial penalization allows to generate binary results. Besides, the body B 1 in Fig. 4.8b shows a hollow, which can be better appreciated in the Fig. 4.9. This cavity has a volume of 34.36 cm 3 (3.436 • 10 -5 m 3 ), representing the 4.32% of the total body volume.

The voxel-based results are post-processing by a body-fitted procedure, the process for proximal and distal links is depicted in Fig. 4.10. Therefore, the complete CAD model of the optimized fivebar linkage is shown in Fig. 4.11, where the right-leg is generated reflecting the left-leg with respect to the Y -axis (see Fig. 4.1). In Table 4.1, the mass and the COM (r i ) of each body of the optimized five-bar mechanism are presented. (f) Elements with a density higher than 0.95 in body B 2 .

Fundamental frequency validation

The first property to be validated is the mechanism fundamental frequency. Figure 4.12 shows the solution of the modal analysis performed in ANSYS. For this case, we are interested in the value of the first natural frequency, which is computed as f 1 ANSY S = 666.77 Hz. Besides, the value of the first natural frequency obtained from the optimization process is f 1 = 681.28 Hz.

Therefore, by comparing the values from ANSYS and the computational platform, the error is 2.17%. This error can be attributed to the post-processing and also due to the differences with the mesh type. In the computational platform all the elements are linear hexahedral, while in ANSYS due to the geometry some elements are thetrahedral, making less stiff the bodies.

Compliance validation

The compliance is a index of how stiff is a system. In this case, the compliance of the linkage is computed using ANSYS performing a static structural analysis, and this value is compared with that obtained in the optimization process.

The total compliance of the whole system computed with ANSYS is 6.21 • 10 -6 J, and the optimized compliance is 6.58 • 10 -6 J, therefore the error is 5.62%. It is an expected result because of the mesh type used in ANSYS. In ANSYS some regions are modeled using thetrahedrons, therefore these elements reduce slightly the total compliance due to they are less stiff than the hexahedrons used in MUBOTO. In addition, the displacements at the joint C are computed using ANSYS, and they are compared with the values obtained with the computational platform. Table 4.2 shows the displacements on the joint C evaluated with ANSYS and those computed by the computational platform. From this table can be observed that the joint deformations are very similar, and comparing the 2-norm of these displacements we obtain a 2.04% of error.

Dynamic validation

In order to evaluate the dynamic performance of the optimized five-bar linkage, the shaking force of the mechanism is evaluated at its base using ADAMS. For this endeavor, a zigzag trajectory was selected in order to move the end-effector inside of its regular workspace. The trajectory is defined in such way that it sweeps the regular workspace of the five-bar linkage as it is shown in The shaking force of the opitmized linkage is presented in Fig. 4.17. This plot shows a very small values of shaking force acting at the base of the linkage, hence the shaking force balancing is validated satisfactorily.

The optimized linkage is compared with an unbalanced linkage, which has links with the same mass as the optimized links, but with their COM located at the geometric center of the link. The resultant reaction loads of the optimized linkage and unbalanced linkage are plotted in Fig. 4.18.

Using the RMS values of the shaking force for the optimal and unbalanced five-bar mechanism presented in Fig. 4.18a, the shaking force reduction is 99.26%.

Regarding to the shaking moment, a considerable difference between the optimized and unbalanced linkages is obtained, as it can be seen in Fig. 4.18b. Even if the methodology described in this work does not deal directly with the shaking moment minimization, the proposed strategy of constrain the mechanism inertia parameters reduce the shaking moment by 28.35%, based on their RMS values.

Moreover, because of the largest values of the mechanism acceleration take place at the beginning and at the end of the motion, higher values of shaking force and moment arise there. Hence, considering the highest mechanism acceleration which occur at the beginning of the motion, the peak value of the shaking force is 0.0209 N, while for the unbalanced case the peak values is 2.7 N, thus it means a reduction of 99.22% for the shaking force. Likewise, the peak value of the shaking moment for the optimized mechanism is 0.239 Nm, and the peak value for the unbalanced case is 0.321 Nm obtaining a reduction of 25.54% with respect to the shaking moment.

Summary

In this chapter the partial dynamic balancing of five-bar robotic mechanism is conducted using multibody topology optimization. The methodology used allows to take into account the entire mechanism elastic behavior in order to conduct its dynamic balancing. The optimization problem is defined as maximizing the mechanism fundamental frequency, constrained by the shaking force balancing conditions, mechanism inertial parameters and the mechanism compliance. The fulfillment of the shaking force balancing conditions ensures a complete force balanced mechanism. In addition, the shaking moment is reduced by constraining the mechanism inertial parameters. On the other hand, the problem formulation requires to use a partial penalization scheme in order to obtain meaningful designs in the case of the natural frequency optimization. This problem formulation was implemented and solved in the MUBOTO computational framework.

In order to prove the feasibility of the proposed approach, numerical validations of the optimized properties were carried out. The results of these numerical validations showed a low error for all the optimized properties. These results demonstrate that the proposed approach is suitable for a comprehensive mechanical design of robotic mechanisms, taking into account its elastic behavior and dynamic balancing. The structural properties were evaluated using ANSYS and the dynamic balancing performance was assessed by ADAMS.

The computed value of the fundamental frequency using ANSYS is 666.77 Hz, while the optimized value obtained with MUBOTO is 681.28 Hz, thereby the percentage error is 2.17%. With respect to the compliance, the obtained error between the model analyzed in ANSYS and the values obtained by MUBOTO is 5.62%. As it was pointed out, these differences can be generated by the post-processing and due to the element type used in the ANSYS analysis. These errors are very low and are acceptable in both cases. In addition, the end-effector displacement was evaluated using ANSYS and the 2-norm displacement value obtained is 4.89 • 10 -7 m, which compared with the computed one with MUBOTO it is obtained a 2.04% of error. Hence, the displacement value is considered low and the error is small. The dynamic balancing of the five-bar mechanism was evaluated using ADAMS. Thus, the mechanism is evaluated moving its end-effector with a predefined trajectory inside of its regular workspace. Then, the shaking force and moment are computed and compared with an unbalanced mechanism. The latter is formed by links with the same mass and inertial parameters like the optimized ones, but with the center of mass between its joints. Therefore, the RMS value of the shaking force magnitude is computed for the optimized mechanism and compared with the unbalanced one, thus obtaining a reduction of 99.56%. Moreover, since the optimization problem deals with the shaking force elimination by proper mass redistribution, the expected value is zero, and the computed RMS value is 0.0089 N. These former results are acceptable (and can be neglected) since the post-processing step can generate numerical errors during the body-fitting operation. Regarding to the shaking moment, the corresponding shaking moment reduction via RMS values is 28.35%.

Based on the results described above, the proposed approach for dynamic balancing using multibody topology optimization has shown to be a feasible approach in order to design dynamically balanced robotic mechanism, while the mechanism elastic behavior is considered.

Chapter 5

Conclusions

Summary

The Chapter 1 presents the state-of-the-art related to this thesis. The main topics exposed are the classical balancing techniques, topology optimization and its applications to structural design of flexible multibody system and dynamic balancing of mechanisms. Analyzing the main drawbacks of the classical balancing technique we found that the elastic behavior of the links were never taken into account, which is an important consideration for advanced robot design. Therefore, structural topology optimization for multibody systems arise as an option to overcome this challenging problem. In the literature review, very few works that faced the dynamic balancing of mechanism by means of a multibody topology optimization approach were found. Therefore in this thesis this problem is studied in depth.

The mathematical formulation for multibody topology optimization was described in Chapter 2. The topology optimization formulation is based on SIMP (density-based) scheme, thus the FE method is used to model the flexibility of the multibody system. Therefore, in order to reduce the computational cost of solving the complete model, static and dynamic condensation techniques were included. Regardless, when the bodies are considered as three-dimensional ones, the computational cost increases drastically, leading us to develop a computational platform (MUBOTO) based on parallel computing where the proposed methodology was implemented.

In order to verify the suitability of the proposed methodology, two mechanism were optimized.

In Chapter 3 a reactionless four-bar linkage was optimized: the linkage compliance was minimized while the dynamic balancing constraints where fulfilled. The numerical validation with commercial software of the optimized four-bar linkage exhibits an excellent performance in terms of dynamic balancing and linkage stiffness. The linkage compliance was properly optimized and the dynamic unbalanced loads at the base of the linkages were almost null. Besides, the experimental validation of the dynamic performance was conducted by mounting the linkage prototype in a platform which is hanged by four cables attached to a fixed frame. In this way, the platform can moves freely, and the unbalanced effects are visible as translations and rotations of the linkage platform when the mechanism is operating: the higher the unbalance, the largest the displacements. The displacements of the linkage platform were computed by video processing and compared with an unbalanced case. After analyzing these results it was concluded that the dynamic performance of the prototype was satisfactory.

The partial dynamic balancing of the five-bar robotic mechanism was conducted in Chapter 4. This mechanism was optimized maximizing its fundamental frequency while it was constrained by the dynamic balancing conditions, the mechanism compliance and its inertial parameters. In this case, the corresponding numerical validations with ADAMS and ANSYS were carried out, obtaining very low error percentages when the results of the proposed approach are compared against the results obtained with the corresponding simulation. The five-bar mechanism can be fully force balanced, but the shaking moment can only be reduced. Otherwise, in order to eliminate the shaking moment additional components are necessary, which in our case is not desirable. Thereby the resultant shaking force in the optimized mechanism was almost negligible, and the shaking moment was reduced by constraining the mechanism inertial parameters. With respect to the structural optimization, the fundamental frequency was increased more than twice the initial value, and the compliance was properly constrained. Moreover, when these optimized properties are compared with the ANSYS values, the resultant percentage of error is very low, which reveals the correctness of the problem formulation.

It is worth to clarify that the presented results were obtained in a computational framework specially designed for the multibody topology optimization problem, which we named MUBOTO. This framework is software independent, written in C++ programming language using the object oriented paradigm, and it is based on the distributed memory model for parallel computing. In consequence, this computational platform has the potential to solve large-scale problems, despite the fact that the results presented here were obtained on a desktop computer. The computational platform is designed to optimize multibody systems in general by means of topology optimization.

Therefore this platform can be used not only for dynamic balancing, but it can be used in order to explore the best design (shape and material distribution) of multibody systems taking into account structural performance indices.

In this research work we studied the possibility to use topology optimization as a tool for reactionless robot design. The optimum designs of the four-bar and five-bar mechanisms were conducted considering them as flexible multibody systems, optimizing structural performance indices while their dynamic balancing conditions were constrained. Comparisons with numerical cross validations showed excellent performances with respect to the shaking force and moment cancellation/reduction, as well as with respect to the defined elastic constraints and objectives. Moreover, the experimental results showed the efficiency of the approach. Therefore, based on the results presented in this work it can be concluded that topology optimization is a suitable and reliable tool for robot design under dynamic balancing constraints.

As it is usual in the research activity, the answers to the main question generate new paths to be explored. Several research directions appear during these years of scientific work, and they are worth to be mentioned in this chapter. The second section is devoted to them.

Perspectives of future research work

The dynamic balancing of MBS using topology optimization is a novel approach, which has several potential directions to be explored. The research work presented in this report covers four fields of study which can be expanded and/or improved. These fields are the followings:

• Extend the formulation of multibody topology optimization.

• Dynamic balancing of different robots: planar and spatial.

• Computational platform and efficient methods.

• Post-processing and physical realization.

These research directions are detailed in the following sections.

Extend the formulation of multibody topology optimization

Although the shaking force balancing conditions are completely independent of the robot configuration, the shaking moment can be a function of the joint variables in some cases, but even more, the structural performance depends on the robot configuration in general. With the aim of conducting a comprehensive structural optimization, the mathematical formulation of multibody topology optimization must be expanded to consider multiple robot configurations.

Besides, the optimization problems described in this thesis are focused on the use of two phases of the same material, i.e., solid and void. Nevertheless in some circumstances it should be helpful to have the possibility to use different materials, namely with the aim to reduce the mechanism footprint by considering materials of different properties (higher density or strength). A first practical attempt was conducted in Chapter 3, where a cylinder of steel was incorporated into the link 2004] a b Planar 3-dof mechanism [Gosselin et al., 2004] Figure 5.1: Reactionless planar parallel robot of 3-DOFs [START_REF] Gosselin | Synthesis and design of reactionless three-degree-of-freedom parallel mechanisms[END_REF]. in order to reduce the link footprint. This cylinder worked as a counterweight, acting as a separated body and its elasticity was not considered. Thereby, a more general formulation is required to incorporate different materials (for SIMP scheme) in multibody topology optimization for dynamic balancing. Additionally, functionally-graded materials (FGMs) could be used to achieve the desired structural performance, without forgetting the challenges that result from this formulation.

Considering different performance indices for a MBS optimization could be an important approach. Thus, optimize a mechanisms for a target stiffness could be an interesting research activity. Furthermore a multi-objective optimization could be tackle in order to optimize the total material cost or even to include some manufacturing constraints.

Dynamic balancing of different robots: planar and spatial

For a comprehensive study of topology optimization as a tool for dynamic balancing of robots, it is necessary to extend the analysis towards planar mechanisms with more than 2-DOFs and spatial mechanism.

Dynamic balancing of planar robots with 3-DOFs

An remarkable case study is the design of the dynamically balanced 3R planar parallel robot with 3-DOFs. For this case, the optimized reactionless four-bar linkage can be used as a leg for building reactionless robots as it was described in (Gosselin et al., 2004, Ricard and[START_REF] Gosselin | On the development of reactionless parallel manipulators[END_REF]. In Fig. 5.1 a schematic representation of the reactionless robot designed with modules of reactionless four-bar linkages is depicted.

Following the same premise, in [START_REF] Wu | Synthesis of reactionless spatial 3-DoF and 6-DoF mechanisms without separate counter-rotations[END_REF] the authors present the synthesis of reactionless spatial 3-DOF and 6-DOF mechanisms using the dynamically balanced four-bar link- ages. Exploring the design of these mechanisms using topology optimization could be an attractive complementary work.

Dynamic balancing of spatial robots

In case of spatial robots, an interesting case study is the Delta robot [START_REF] Clavel | Device for the movement and positioning of an element in space[END_REF], which is typically used for high-speed operations mainly pick-and-place tasks. The kinematic architecture depicted in Fig. 5.2 is a simplified version of the Delta robot, which allows to perform spatial motions due to the three translational DOFs. This figure represents a proposal for force balancing of the Delta robot (Van der [START_REF] Van Der Wijk | Dynamic balancing of Claver's Delta robot[END_REF], with the pink balls representing counterweights.

Using the proposed methodology, an optimized design of the Delta robot under shaking force balancing conditions can be conducted. Because of the three legs are equal, the optimization should be performed for only one leg, but taking into account the whole structure.

Computational platform and efficient methods

The computational platform developed for this project is implement for large-scale analysis. Nevertheless, there was no opportunity to test the platform in a High Performance Computing (HPC)

center. Analyze large-scale models allows to generate high-quality realistic designs, owing to the mesh size that can be use. Therefore, conduct a large-scale analysis in a HPC center, of any of the study cases described in the thesis is an appealing task.

Regarding the computational methods, it is worth exploring the possibility of including the multigrid method (MG) in the solution of equilibrium equations arising from the FEM. The MG method should reduce the computational time [START_REF] Amir | On multigrid-CG for efficient topology optimization[END_REF], but it must be adequate to the model reduction techniques (static/dynamic condensation) used in this work.

Besides, it would be interesting to examine the Linearization Method (LM), used as optimizer for the topology optimization problem in [START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF]. This motivation arises because of [START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF] indicated that the linearization method is less complex in its implementation than the MMA optimizer, and the results obtained with LM are competitive with respect to MMA results, in terms of convergence time. In this case, precise studies of the results obtained with LM and MMA should be conducted with the aim to perform a proper comparison between these optimizers, mainly in parallel computing applications.

It will be significant to analyze which are the most time consuming processes of the proposed methodology in order to define which is the best approach to speed up the optimization. In this manner it is possible to define if it is worth to use a different hardware (GPU) or alternative algorithms. Other important direction to be explored is the use of methods based on data-driven training to accelerate the optimization process, mainly Neural Networks and more broadly Machine Learning.

Post-processing and physical realization

Post-processing implies to convert the voxel-based results obtained with topology optimization to a smooth design suitable for a CAD software. Because the optimized shape is defined in voxel-based the smooth design is not exactly the same as the one optimized. In [START_REF] Liu | A realization method for transforming a topology optimization design into additive manufacturing structures[END_REF] the authors propose a method for post-processing and physical realization of structures coming from TO, but they do no evaluate the optimized properties and the examples are for 2D problems. Thereby, even if we already conduct successful post-processing procedure, comprehensive studies should be carried out in order to define the most appropriated scheme to generate the lowest error in the mesh conversion process. These studies should be done with the aim of preserving the optimized properties within an acceptable error range, at least for the simulation stage.

As it is described in Appendix B, when the aim is to build a functional prototype of dynamically balanced linkage, the designer must choose printing solid1 3D parts. Otherwise, the center of mass is modified due to the material distribution during the printing process. Therefore, it would be interesting to study the possibility to develop a tool for 3D printing which allows to print non-solid (infill percentage lower than 100%) objects, while the total mass and the location of the COM of the designed/optimized body is preserved. Of course, the structural performance will be different, but the dynamic balancing will be maintained. where the position vector r i/ j represents the location of point i with respect to point j. Then, using this expression the following two scalar equations derived:

ℓ 1 cos θ 1 + ℓ 2 cos θ 2 -ℓ 3 cos θ 3 -ℓ 4 cos θ 4 -d = 0 (A.2) ℓ 1 sin θ 1 + ℓ 2 sin θ 2 -ℓ 3 sin θ 3 -ℓ 4 sin θ 4 = 0 (A.3)

With the aim to express the angular variables as a function of the independent parameters Eq. (A.3) is solved for sin θ 3 , and Eq. (A.2) for cos θ 3 , therefore their corresponding equations are obtained: cos θ 3 = ℓ 1 cos θ 1 + ℓ 2 cos θ 2 -ℓ 4 cos θ 4 -d ℓ 3 (A.4) sin θ 3 = ℓ 1 sin θ 1 + ℓ 2 sin θ 2 -ℓ 4 sin θ 4 ℓ 3 (A.5) Equations (A.4) and (A.5) are function of the dependent variable θ 2 . In order to find a solution for cos θ 2 and sin θ 2 , it is necessary to use the trigonometric identity sin θ 2 3 + cos θ 2 3 = 1, which leads to an equation of the form ∆ 1 cos θ 2 + ∆ 2 sin θ 2 = ∆ 3 . The analytic solution of the latter equation is Solving the latter equations for θ2 and θ3 , the following equations are obtained: θ2 = -sin(θ 3 -θ 1 ) sin(θ 3 -θ 2 ) ℓ 1 ℓ 2 θ1 -sin(θ 4 -θ 3 ) sin(θ 3 -θ 2 ) ℓ 4 ℓ 2 θ4 (A.14) θ3 = -sin(θ 2 -θ 1 ) sin(θ 3 -θ 2 )

ℓ 1 ℓ 3 θ1 -sin(θ 4 -θ 2 ) sin(θ 3 -θ 2 ) .15) where the dependent angular velocities are expressed as function of the velocities of the actuated joints, namely θ1 and θ4 .

ℓ 4 ℓ 3 θ4 (A
On the other hand, in order to compute the velocities of the COM of each body, Eqs. (A.5) and (A.4) are substituted into the vector r S 3 /A (Eq. A.11) with the purpose to reduce one dependent variable. Hence, computing the time derivatives of the position vectors given in the resultant Eq. (A.16)

A.2 Shaking force formulation

The shaking force balancing conditions require to keep the linkage center of mass stationary during the linkage operation. The shaking force f sh generated at the base of the mechanism is computed as the time rate of its linear momentum, thus:

f sh = dp dt (A.17)
When the linear momentum of the system is constant for any trajectory, the resultant reaction forces at the linkage base sum up to zero. The linkage is composed of four moving bodies, thus its linear momentum computed with respect to point A is given by: .18) where v S i /A is the velocity of the link's COM given by S i , and it is defined with respect to the origin of the reference system, point O. Therefore, the time-derivative of the linear momentum is expressed as: with vS i /A as the acceleration of the link's COM. Therefore, in order to obtain the shaking force balancing conditions, the following equation must be fulfilled:

p = m 1 v S 1 /A + m 2 v S 2 /A + m 3 v S 3 /A + m 4 v S 4 /A (A
f sh = dp dt = 0 (A.20)
In the next sections, a detailed computation of the linear momentum is described, as well as the derivation of the shaking force balancing conditions.

A.2.1 Shaking force balancing conditions

The analytical expression of the linear momentum of the five-bar linkage is given as follows: In order to obtain the shaking force balancing conditions without superfluous constraints, the expression for the linear momentum must include only independent variables and their derivatives.

Thus, the angular velocity θ2 computed in Eq. (A.14) is substituted in Eq. (A.21). Afterwards, a compact expression for the linear momentum is formed by grouping similar terms and substituting the trigonometric identities in Eq. (A.21). The resultant expression is written in the following form: with λ 12 = ℓ 1 /ℓ 2 and λ 42 = ℓ 4 /ℓ 2 . Besides, the vector Ω Ω Ω is composed as:

Ω Ω Ω = θ1 sin θ 1 θ1 cos θ 1 θ1 k 1 sin θ 2 θ1 k 1 cos θ 2 θ4 k 2 sin θ 2 θ4 k 2 cos θ 2 θ4 sin θ 4 θ4 cos θ 4 T , (A.23) and the parameters k i are given by: k 1 = sin(θ 3 -θ 1 ) sin(θ 3 -θ 2 ) , k 2 = sin(θ 4 -θ 3 ) sin(θ 3 -θ 2 ) . (A.24)

Therefore, the vector Ω Ω Ω is expressed as a function of the input angular velocities θ1 and θ4 , and the dependent angular variables are already defined in Section A.1.1. As a consequence, the terms a i and b i which appear in Eq. (A.22) are constant coefficients and they are give by: a 1 = -m 1 r 1 cos ψ 1 + m 2 ℓ 1 + ℓ 1 ℓ 3 m 3 r 3 cos ψ 3 , b 1 = -m 1 r 1 sin ψ 1 + ℓ 1 ℓ 3 m 3 r 3 sin ψ 3 a 2 = m 2 r 2 cos ψ 2 + ℓ 2 ℓ 3 m 3 r 3 cos ψ 3 , b 2 = m 2 r 2 sin ψ 2 + ℓ 2 ℓ 3 m 3 r 3 sin ψ 3 a 3 = -m 4 r 4 cos ψ 4 - ℓ 4 ℓ 3 m 3 r 3 cos ψ 3 + m 3 ℓ 4 , b 3 = -m 4 r 4 sin ψ 4 - ℓ 4 ℓ 3 m 3 r 3 sin ψ 3 (A.25)

Consequently, because of vector Ω Ω Ω (Eq. A.23) depends on the angular velocities and joint variables, the linear momentum given by Eq. (A.22) cannot generate a constant value for any arbitrary mechanism motion. Hence, in order to fulfill the requirement for shaking force balancing dp/dt = 0, the constant coefficients of Eq. (A.22) must be zero (with λ 1 = 0 and λ 2 = 0). Thereby, taking into account this condition, the general shaking force balancing conditions are obtained: m 1 r 1 cos ψ 1 + ℓ 1 ℓ 3 m 3 r 3 cos ψ 3 + m 2 ℓ 1 = 0 (A.26) m 1 r 1 sin ψ 1 + ℓ 1 ℓ 3 m 3 r 3 sin ψ 3 = 0 (A.27) using a grid pattern, with its internal material distribution depicted in Fig. B.1. As a result the material distribution is not homogeneous and it generates a COM in the printed part different from the optimized one.

Therefore, the main problem is the impossibility to generate a uniform material distribution using an infill percentage different from 100%. The algorithms are designed to generate proper external shape, while the internal shape is generated by specific patterns with the aim of influence in the strength of the printed part. Nevertheless, usually these algorithms do not preserve the inertial properties of the designed object. Obviously, the most straightforward solution is printing with a full infill, but it can no be a practical solution.

Consequently, one solution relies on the algorithms for path planing of the 3D printing machine.

It will be necessary to include a new constraints on the material distribution algorithm, forcing to keep the total mass and the center of mass in the printed part equal to the designed part. This constraint must ensure the printability of the optimum mechanism. [START_REF] Prévost | Make it stand: Balancing shapes for 3D fabrication[END_REF] developed a semi-automatic software for balancing printed objects by modifying the internal shape and the total volume, while keeping some features of the external shape. In our case, the external shape should be preserved or some metrics should be defined in order to produce a functional prototype.

B.3 Traditional manufacturing process

When the shapes of the optimized bodies do not have hollows, and for specific curvatures, the traditional manufacturing process are useful. The viability of the traditional Computerized Numerical Control (CNC) manufacturing should be verified by Computer-Aided Manufacturing (CAM) software. Nevertheless, because some intricate features on the part can be complex but able to be machined, the realization of the part will depend on the expertise of the workshop technician, and the features of the CNC machine.

In the case of the four-bar linkage, due to the loading conditions, the optimized bodies do not present hollows, thus the they were able to manufacture using a CNC machine. Besides, the parts created with the traditional manufacturing process using isotropic materials, has the advantage to be able for structural validation, since the FE model and the machined part have the same mechanical properties. The suitability of the proposed methodol-ogy is confirmed by accomplishing an optimized design of a reactionless four-bar linkage and the partial dynamic balancing of five-bar robotic mechanism. The significance of the dynamically balanced four-bar linkage is related to the possibility to exploit this optimized linkage as a special leg for building reactionless robots. Besides, the five-bar robot is very important due to its industrial applications, where it is typically used in pick-and-place operations.

  3.3 The three kinds of dynamically balanced four-bar linkage. . . . . . . . . . . . . . 3.4 Initial design domain of four-bar links. All dimensions are in millimeters. . . . . . 3.5 Initial design domain of bodies B 1 and B 3 (thickness of 32 mm) when the cylinder is included. All dimensions are in millimeters. . . . . . . . . . . . . . . . . . . . . 3.6 Compliance evolution for the four cases. . . . . . . . . . . . . . . . . . . . . . . . 3.7 Voxel-based results of body B 1 for Case I. . . . . . . . . . . . . . . . . . . . . . . 3.8 Linkage footprint reduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Post-processing of the optimized body B 1 for Case I. . . . . . . . . . . . . . . . . 3.10 Numerical validation of total deformations using ANSYS. . . . . . . . . . . . . . 3.11 Numerical validation of dynamic balancing using ADAMS. . . . . . . . . . . . . . 3.12 Cavities in an optimized body. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Prototype of the optimized reactionless four-bar linkage. . . . . . . . . . . . . . . 3.14 Experimental setup for the dynamically balanced four-bar linkage. . . . . . . . . . 3.15 The four-bar linkage in its home position for the two scenarios. . . . . . . . . . . . 3.16 Experimental platform displacements for the balanced and unbalanced four-bar linkages. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Scheme of the five-bar mechanism. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Five-bar mechanism with bodies of generic shape. . . . . . . . . . . . . . . . . . . 4.3 Shaking force balanced five-bar mechanism. . . . . . . . . . . . . . . . . . . . . . 4.4 Initial design domain of five-bar links. The figures have different scales, and all dimensions are in millimeters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Full density field of body B 1 without partial penalization. . . . . . . . . . . . . . . 4.6 Body division for partial penalization. . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Evolution of the fundamental frequency for the optimization of five-bar mechanism. 4.8 Voxel-based results of bodies B 1 and B 2 . . . . . . . . . . . . . . . . . . . . . . . 4.9 Description of the cavity in body B 1 . The body is represented by a wireframe for visual purposes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv 4.10 Post-processing of the optimized bodies B 1 and B 2 . . . . . . . . . . . . . . . . . 4.11 CAD model of optimized five-bar linkage. . . . . . . . . . . . . . . . . . . . . . . 4.12 Total deformation for the first natural frequency. . . . . . . . . . . . . . . . . . . . 4.13 Total deformation for static structural analysis. . . . . . . . . . . . . . . . . . . . . 4.14 Trajectory for the dynamic balancing evaluation. The workspace of the five-bar mechanism is a rectangle centered in (0, 110) mm of width equal to 120 mm and height equal to 40 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Translational velocity and acceleration of the end-effector. . . . . . . . . . . . . . 4.16 Motors angular velocity and acceleration for the zig-zag trajectory. . . . . . . . . . 4.17 Shaking force at the base of the optimized five-bar mechanism. . . . . . . . . . . . 4.18 Comparative shaking force and moment at the base of the optimized and unbalanced five-bar mechanisms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Reactionless planar parallel robot of 3-DOFs (Gosselin et al., 2004). . . . . . . . . 5.2 Force balanced Delta robot (Van der Wijk and Herder, 2009). . . . . . . . . . . . . A.1 Schematics of the five-bar linkage. . . . . . . . . . . . . . . . . . . . . . . . . . . B.1 Material distribution for 3D printed part. . . . . . . . . . . . . . . . . . . . . . . . xvi 3.1 Material properties. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Optimization results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Properties of the optimized reactionless four-bar linkage (Case I). . . . . . . . . . . 3.4 Joint deformations. All quantities are in meters. . . . . . . . . . . . . . . . . . . . . 3.5 Reaction loads comparison using their RMS values. . . . . . . . . . . . . . . . . . . 4.1 Properties of the optimized five-bar mechanism. . . . . . . . . . . . . . . . . . . . . 4.2 Joint C deformations. All quantities are in meters. . . . . . . . . . . . . . . . . . . xvii

Figure 1 :

 1 Figure1: The phases in design[START_REF] Budynas | Mechanical Engineering Design[END_REF].

( a )

 a Unbalanced four-bar linkage. (b) Balanced four-bar linkage.

Figure 1 . 1 :

 11 Figure 1.1: Shaking force balancing by adding counterweights.

Figure 1 . 2 :

 12 Figure 1.2: Shaking force balancing by adding auxiliary structures.

( a )

 a Unbalanced four-bar linkage. (b) SFB by adjusting kinematic parameters.

Figure 1 . 3 :

 13 Figure 1.3: Shaking force balancing by adjustment of kinematic parameters.

  SMB by Gears with internal teeth .

Figure 1 . 4 :

 14 Figure 1.4: Shaking moment balancing by counter-rotation (and counterweight)[START_REF] Herder | A counter-rotary counterweight (CRCW) for light-weight dynamic balancing[END_REF]).

  Figure 1.7: Dynamic optimization of the planar four-bar linkage[START_REF] Feng | A new optimization method for dynamic design of planar linkage with clearances at joints-optimizing the mass distribution of links to reduce the change of joint forces[END_REF].

( a )

 a Optimized shapes of four-bar linkage. (b) Optimized shapes of five-bar linkage.

Figure 1 . 8 :

 18 Figure 1.8: Shape optimization of planar four-bar and five-bar linkages for dynamic balancing (Chaudhary and Chaudhary, 2015a). Initial linkage is drawn by a dashed line and the optimized linkage by solid line.

Figure 1 . 9 :

 19 Figure 1.9: Typical process of structural topology optimization.

( a )

 a Optimized beam without filter. (b) Optimized beam filtered with sensitivity filter.

Figure 1 .

 1 Figure 1.10: Optimized MBB beam for minimum compliance constrained by the 40% of material.

( a )

 a Schematic diagram of the robotic manipulator. (b) Optimized arm.

( a )

 a Initial mechanism. (b) Optimized mechanism.

( a )

 a Two dimensional model. (b) Three dimensional model.

Figure 1 . 14 :

 114 Figure 1.14: Topology optimization of five-bar mechanism considering technology-oriented performance indices (Briot and Goldsztejn, 2018b).

Figure 1 . 15 :

 115 Figure 1.15: Structural topology optimization process for flexible multibody systems.

( a )

 a Optimal design of body 1. (b) Optimal design of body 2.(c) Optimal design of body 3.

Figure 1 . 16 :

 116 Figure 1.16: Topology optimization of four-bar linkage under balancing constraints(Briot and Goldsztejn, 2018a).

  Figure 2.1: General scheme of the body B i .

Figure 2 . 2 :

 22 Figure 2.2: Interface nodes for bodies connection.

Figure 2 . 3 :

 23 Figure 2.3: The finite element i j of the body B i .

The

  Figure 2.4: MUBOTO framework.

Figure 3 . 1 :

 31 Figure 3.1: Scheme of the four-bar linkage.

  Reactionless spatial 3-DOF mechanism. (b) Reactionless spatial 6-DOF parallel robot.

Figure 3 . 3 :

 33 Figure 3.3: The three kinds of dynamically balanced four-bar linkage.

  Initial design domain of bodies B 1 and B 3 . Thickness of 32 mm. Initial design domain of body B 2 . Thickness of 10 mm.

Figure 3 . 4 :

 34 Figure 3.4: Initial design domain of four-bar links. All dimensions are in millimeters.

Figure 3 . 5 :

 35 Figure 3.5: Initial design domain of bodies B 1 and B 3 (thickness of 32 mm) when the cylinder is included.All dimensions are in millimeters.

Figure 3 . 6 :

 36 Figure 3.6: Compliance evolution for the four cases.

( a )

 a Full density field of body B 1 , ranging from 0 to 1.(b) Internal section view.

Figure 3 . 7 :

 37 Figure 3.7: Voxel-based results of body B 1 for Case I.

  (a) Case I: Optimized body B 1 . (b) Case III: Optimized body B 1 with cylinder. (c) Overlapping of results: Case I and III. (d) Case II: Optimized body B 1 .(e) Case IV: Optimized body B 1 with cylinder.(f) Overlapping of results: Case II and IV.

Figure 3 . 8 :

 38 Figure 3.8: Linkage footprint reduction.

( a )

 a Optimized body represented by structured mesh. (b) Optimized body represented by triangular mesh. (c) Optimized body in CAD software.

Figure 3 . 9 :

 39 Figure 3.9: Post-processing of the optimized body B 1 for Case I.

Figure 3 .

 3 Figure 3.10: Numerical validation of total deformations using ANSYS.

  Comparative shaking moment.

Figure 3 . 11 :

 311 Figure 3.11: Numerical validation of dynamic balancing using ADAMS.

  Section view with cavities.

Figure 3 . 12 :

 312 Figure 3.12: Cavities in an optimized body.

Figure 3 . 13 :

 313 Figure 3.13: Prototype of the optimized reactionless four-bar linkage.

  Figure 3.14: Experimental setup for the dynamically balanced four-bar linkage.Origin O 1 Y

  Four-bar linkage with additional mass.

Figure 3 .

 3 Figure 3.15: The four-bar linkage in its home position for the two scenarios.

  Balanced four-bar linkage rotation around the Z axis.

  Unbalanced four-bar linkage rotation around the Z axis.

Figure 3 . 16 :

 316 Figure 3.16: Experimental platform displacements for the balanced and unbalanced four-bar linkages.

Figure 4 . 1 :

 41 Figure 4.1: Scheme of the five-bar mechanism.

A

Figure 4 . 2 :

 42 Figure 4.2: Five-bar mechanism with bodies of generic shape.

  2c) and (4.2d) are the necessary conditions to locate the COM of bodies B 2 and B 3 on joints B and D respectively, and they are derived from Eq. (4.1c).

Figure 4 . 3 :

 43 Figure 4.3: Shaking force balanced five-bar mechanism.

;

  (a) Initial design domain of body B 1 . Thickness of 50 mm. ; (b) Initial design domain of body B 2 . Thickness of 20 mm.

Figure 4 . 4 :

 44 Figure 4.4: Initial design domain of five-bar links. The figures have different scales, and all dimensions are in millimeters.

  4, and the dimension between the fixed joints A and E is d = 112 mm. The initial design domain for the links is show in Fig. 4.4. Each link has two joints denoted by the holes of 10 mm of diameter.

Figure 4 . 5 :Figure 4 . 6 :

 4546 Figure 4.5: Full density field of body B 1 without partial penalization.

  Figure 4.7: Evolution of the fundamental frequency for the optimization of five-bar mechanism.

( a )

 a Full density field of body B 1 ranging from 0 to 1. (b) Internal section view of body B 1 .(c) Elements with a density higher than 0.95 in body B 1 .(d) Full density field of body B 2 ranging from 0 to 1. (e) Internal section view of body B 2 .

Figure 4 . 8 :

 48 Figure 4.8: Voxel-based results of bodies B 1 and B 2 .

( a )

 a Top view of body B 1 . (b) Perspective view of body B 1 . (c) Front view of body B 1 . (d) Right view of body B 1 .

Figure 4 . 9 :

 49 Figure 4.9: Description of the cavity in body B 1 . The body is represented by a wireframe for visual purposes.

( a )

 a Optimized body B 1 represented by structured mesh. (b) Optimized body B 1 represented by triangular mesh. (c) Optimized body B 1 in CAD software. (d) Optimized body B 2 represented by structured mesh. (e) Optimized body B 2 represented by triangular mesh. (f) Optimized body B 2 in CAD software.

Figure 4 .

 4 Figure 4.10: Post-processing of the optimized bodies B 1 and B 2 .

Figure 4 .

 4 Figure 4.11: CAD model of optimized five-bar linkage.

Figure 4 .

 4 Figure 4.12: Total deformation for the first natural frequency.

Figure 4 .

 4 Figure 4.13: Total deformation for static structural analysis.

  Fig. 4.14. Besides, Fig. 4.15 shows the magnitude of translational velocity and acceleration of the end-effector, while Fig.4.16 presents the corresponding angular velocity and acceleration of the motors.

Figure 4 .

 4 Figure 4.14: Trajectory for the dynamic balancing evaluation. The workspace of the five-bar mechanism is a rectangle centered in (0, 110) mm of width equal to 120 mm and height equal to 40 mm.

  Magnitude of the translational acceleration of the endeffector.

Figure 4 . 15 :

 415 Figure 4.15: Translational velocity and acceleration of the end-effector.

Figure 4 . 16 :

 416 Figure 4.16: Motors angular velocity and acceleration for the zig-zag trajectory.

Figure 4 . 17 :

 417 Figure 4.17: Shaking force at the base of the optimized five-bar mechanism.

Figure 4 . 18 :

 418 Figure 4.18: Comparative shaking force and moment at the base of the optimized and unbalanced five-bar mechanisms.

Figure 5 . 2 :

 52 Figure 5.2: Force balanced Delta robot (Van der Wijk and Herder, 2009).

  Five-bar mechanism with bodies of general shape.

Figure A. 1 :

 1 Figure A.1: Schematics of the five-bar linkage.

-

  (A.11), we obtain the velocities of the COM of each link:v S 1 /A ≡ ṙS 1 /A =    -r 1 sin(θ 1 + ψ 1 ) θ1 r 1 cos(θ 1 + ψ 1 ) θ1    , v S 2 /A ≡ ṙS 2 /A =    -θ1 ℓ 1 sin θ 1 -θ2 r 2 sin(θ 2 + ψ 2 )θ1 ℓ 1 cos θ 1 + θ2 r 2 cos(θ 2 + ψ 2 ) θ 4 + r 3 sin(θ 4 + ψ 3 ) ℓ 3 cos θ 4 -r 3 cos(θ 4 + ψ 3 ) θ4 r 4 sin(θ 4 + ψ 4 ) θ4 r 4 cos(θ 4 + ψ 4 )

  dp dt = m 1 vS 1 /A + m 2 vS 2 /A + m 3 vS 3 /A + m 4 vS 4 /A (A.19)

  θ 4 + r 3 sin(θ 4 + ψ 3 ) ℓ 3 cos θ 4 -r 3 cos(θ 4 + ψ 3 )

  1 a 2 λ 12 b 2 λ 12 a 2 λ 42 b 2 λ 42 a 3 b 3 b 1 -a 1 b 2 λ 12 -a 2 λ 12 b 2 λ 42 -a 2 λ 42 b 3 -a 3

Titre:

  Optimisation de la topologie des robots dynamiquement équilibrés Mot clés : équilibrage dynamique, robots équilibré dynamiquement, optimisation topologique multi-corps, mécanisme à quatre barres équilibré dynamiquement, équilibrage dynamique partiel d'un mécanisme à cinq barres Résumé : L'équilibrage dynamique est un domaine d'étude important en robotique à grande vitesse et pour la robotique spatiale. La prise en compte des performances d'équilibrage dynamique des robots lors de leur conception permet de réduire les vibrations du bâti, une meilleure précision et des temps de cycle réduits. Dans le but de développer une méthode de conception de robot équilibrés dynamiquement, l'optimisation topologique structurelle est étudiée dans ce travail de recherche en tant qu'outil pour concevoir des robots équilibrés dynamiquement, c'est à dire un robot ne transmettant pas de réaction sur le châssis. La pertinence de la méthodologie proposée est confirmée par la réalisation d'une conception optimisée d'un mécanisme à quatre barres équilibré dynamiquement et l'équilibrage dynamique partiel d'un mécanisme à cinq barres. L'importance de la conception d'un mécanisme à quatre barres équilibré dynamiquement est liée à la possibilité d'exploiter ce mécanisme optimisé comme un module utile à la conception de robots équilibrés dynamiquement. Le mécanisme à cinq barres a une importance toute particulière en raison de ses applications industrielles, où il est généralement utilisé dans des opérations de prise et de dépose. Title: Topology optimization of reactionless robots Keywords: dynamic balancing, reactionless robots, multibody topology optimization, reactionless four-bar linkage, partial dynamic balancing of five-bar mechanism Abstract: Dynamic balancing is an important field of study in high-speed robotics and spatial robots. Taking into account robot dynamic balancing performance for robot design leads to low base vibrations, high precision and short cycle times. With the aim to develop a comprehensive robot design for dynamic balancing, structural topology optimization is studied in this research work as a tool for designing dynamically balanced robots, also called reactionless robots.

  . . . . . . . . . . . . . . . . . . 2 1.1 Shaking force balancing by adding counterweights. . . . . . . . . . . . . . . . . . 8 1.2 Shaking force balancing by adding auxiliary structures. . . . . . . . . . . . . . . . 8 1.3 Shaking force balancing by adjustment of kinematic parameters. . . . . . . . . . . 9 1.4 Shaking moment balancing by counter-rotation (and counterweight)

  General scheme of the body B i . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Interface nodes for bodies connection. . . . . . . . . . . . . . . . . . . . . . . . . 2.3 The finite element i j of the body B i . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 MUBOTO framework. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Scheme of the four-bar linkage. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Design of reactionless robot using a reactionless four-bar linkage as a building leg

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 xiv 2.1

Table 3 .

 3 1: Material properties.

	Material	Density [kg/m 3 ]	Young's modulus [GPa]	Poisson ratio
	Nylon (PA66)	1150	2.76	0.39
	Aluminium	2700	68.9	0.33
	Steel	7800	210	0.3

Table 3 .

 3 2: Optimization results.

	Case	Time [h]	Max. const. violation [%]	f (ρ ρ ρ)
	I	23.63	6.64 • 10 -4	1.24
	II	20.91	1.09 • 10 -4	1.16
	III	23.05	1.69 • 10 -4	1.13
	IV	20.94	1.13 • 10 -4	1.13

Table 3 . 3 :

 33 Properties of the optimized reactionless four-bar linkage (Case I).

	Body	mass (kg)	length (m)	x i (m)
	1	0.134808	0.06	-0.012414
	2	0.044840	0.2	0.102785
	3	0.133670	0.06	-0.012179

Table 3 .

 3 4: Joint deformations. All quantities are in meters.Table3.5: Reaction loads comparison using their RMS values.

	Joint	X	Y	ANSYS	Z	Norm
	O 2	4.3685 • 10 -5	2.5500 • 10 -9	8.8538 • 10 -6	4.4573 • 10 -5
	O ′ 2	3.6929 • 10 -5	-2.5008 • 10 -5	4.4942 • 10 -6	4.4826 • 10 -5
				MUBOTO
		X	Y			Z	Norm
	O 2	4.3097 • 10 -5	2.4891 • 10 -9	4.1009 • 10 -6	4.3292 • 10 -5
	O ′ 2	3.6673 • 10 -5	-2.3787 • 10 -5	4.7843 • 10 -6	4.3973 • 10 -5
	Reaction load (magnitude)	Unbalanced linkage	Reactionless linkage	Reduction (%)
	Shaking force [N]	0.0951			0.0026	97.26
	Shaking moment [Nm]	0.0083			8.46 • 10 -5	98.98

Table 4 . 1 :

 41 Properties of the optimized five-bar mechanism.

	Body	mass [kg]	length [m]	r i [m]
	1	6.1923	0.09	-0.01931
	2	1.2912	0.09	-0.0004
	3	1.2912	0.09	-0.0004
	4	6.1923	0.09	-0.01931

Table 4 . 2 :

 42 Joint C deformations. All quantities are in meters.

1.2. CLASSICAL DYNAMIC BALANCING TECHNIQUES

In specific cases, the function mldivide automatically execute multiple threads, without requiring the user specification (MatLab,

2022).

A movie of the camera records is available here: https://uncloud.univ-nantes.fr/index.php/s/xfzrJMipy6zJAYN

The dextrous workspace is the set of locations of the end-effector for which any orientations can be reached(Merlet,

2006).

Usually for 3D printing by fused deposition modeling the amount of material of the printed part is controlled by an infill percentage, where 100% generates a solid object. Commonly this parameter is used below of 60%.

Acknowledgments

I would like to express my gratitude to all those who were involved during these years of my PhD studies.

Contents

I N this chapter partial dynamic balancing of the five-bar robotic mechanism is conducted using topology optimization. The optimization problem is addressed by maximizing the mechanism first natural frequency, while it is constrained by: (a) shaking force balancing conditions, (b) mechanism inertial properties, and (c) mechanism compliance. Constraining the mechanism inertial properties, the shaking moment generated at the frame is reduced. Furthermore, in order to obtain meaningful designs for the natural frequency optimization, partial penalization scheme was applied to the shaking force balancing constraints and also to the mechanism inertial parameters. This Appendix presents a formulation to obtain the shaking force and moment balancing conditions of the five-bar mechanism. The shaking force balancing conditions make zero the resultant force on the base of the linkage for any arbitrary trajectory. Nevertheless, for the case of shaking moment it is not possible to obtain balancing conditions which produce zero unbalanced moment on the linkage base. Thus, shaking moment balancing is only possible using external components.

The five-bar linkage is a two degree of freedom mechanism, composed of five revolute joints of parallel axes. A general schematic representation of the five-bar mechanism is depicted in Fig.

A.1a, where angle θ i describes the orientation of each body B i with respect to the X 0 -axis of the fixed frame. The parameters required to describe the center of mass of each link are shown in Fig.

A.1b. The points S i represent the COM of each link and m i is the total mass of the link. The location of the center of mass S i is defined by the distance r i and the constant angle ψ i .

A.1 Kinematic analysis

In order to compute the linear momentum and angular momentum, the kinematic analysis of the five-bar linkage is conducted, and it involves the position and velocity analysis. The complete kinematic solution of the five-bar linkage is expressed as a function of independent variables, namely θ 1 , θ 2 and its angular velocities θ1 and θ2 .

given as follows:

where:

It should be noted that the terms ∆ i are function only of the independent variables and constant parameters. The sign of the square root will define the assembly mode of the linkage.

Besides, once the position analysis is completed, from Fig.

A.1 we can define the position vectors of the center of mass of each link:

which can be readily computed.

A.1.2 Velocity analysis

The velocity analysis is carried out by differentiating with respect to time the loop-closure equation defined in Eq. (A.2) and (A.3). Thereby the following is obtained:

When the five-bar linkage is considered as a "in-line" linkage (i.e., the angles ψ i only takes values of 0 or π), the resultant shaking force balancing conditions are:

Moreover, if we choose a linkage with equal link lengths, we obtain:

m 2 r 2 cos ψ 2 + m 3 r 3 cos ψ 3 = 0 (A.36) m 4 r 4 cos ψ 4 -m 3 r 3 cos ψ 3 + m 3 ℓ 4 = 0 (A.37)

Using the shaking force balancing conditions derived in this section, different types of force balanced linkages can be designed.

A.3 Shaking moment formulation

The shaking moment is computed as the time-derivative of the total angular momentum of the system. The angular momentum is computed with respect to the reference point A, which coincides with the origin of the reference frame. Hence, the angular momentum h A of the five-bar linkage is given by: .38) where h i/ j defines the angular momentum of each moving body, identified by the body number i and its reference point j. The angular momentum for each body is expressed as follows:

Where the angular momentum of each body, h S i , is defined with respect to its center of mass, and it is given as follows:

where I S i is the moment of inertia of the moving body B i , around its centroidal axis, which is perpendicular to the plane of motion. Besides, the absolute angular velocity of each body is given by:

where vector n is a unit vector, normal to the plane of motion.

The shaking moment generated at the base of the mechanism, m sh , is computed as the timederivative of the angular momentum, thus it is give by:

Hence, in order to design a shaking moment balanced mechanism, the requirement is to have a constant angular momentum or zero shaking moment at the base, thus the following equation must be fulfilled:

In the following section the angular momentum is studied in order to determine if a suitable set of shaking moment balancing conditions can be derived from this equation.

A.3.1 Angular momentum of the five-bar mechanism

Using the results of the kinematic analysis of the five-bar linkage, and considering a planar motion normal to the z-axis, the total angular momentum is computed, resulting in the following scalar equation:

where can be notice that the angular velocities θ2 and θ3 appears, and these are functions of the independent velocities θ1 and θ4 , as can be seen in Eqs. (A.14) and (A.15).

Inspecting the angular momentum of the linkage given in Eq. (A.47), it is reveled that it cannot be constant, because it depends on the linkage angular velocities. The only possibility is to operate the linkage with a constant angular velocities, and at the same time have a constant or zero value for the remaining variable terms. These conditions are almost impossible for practical applications. Therefore, the complete shaking moment balancing cannot be achieved by mass redistribution, it is necessary to use external components like counter-rotations.

Appendix B Notes on the manufacture of dynamically balanced mechanisms

This Appendix has the purpose to describe some insights for the manufacturing process of the optimized bodies coming from the dynamically balanced linkages using topology optimization.

These notes are result of the prototyping phase, when several questions arise and they are described here.

As it is natural for topology optimization, the first manufacturing technology to be considered is Additive Manufacturing (AM). In our case, we describe Fused Deposition Modeling (FDM), because it is the most widely used and affordable technique for AM. Indeed, a prototype built using AM based on FDM technology has the aim to validate only the dynamic balancing performance. The structural performance is more challenging to evaluate using 3D printed parts. The FE model for the optimization problem is considered as an isotropic elastic model, while the 3D printing process generates orthotropic properties in the printed part. In general, due to AM is generated by printing layer by layer, the mechanical properties of the printed parts are affected, presenting anisotropy, porosity and residual stress [START_REF] Meng | From topology optimization design to additive manufacturing: Today's success and tomorrow's roadmap[END_REF].

B.1 Additive manufacturing based on FDM

There are different AM technologies, among them is the fused deposition modeling. The FDM process consist in generating a part layer by layer until a three-dimensional object is created. These layers are produced by heating thermoplastic in filament form and then deposited onto a surface.

The most common thermoplastic filament are Acrylonitrile Butadiene Styrene (ABS) and Polylactic Acid (PLA). In the FDM process the exterior surface is generated ensuring the dimensionality of the part, but the interior of part can be printed completely solid or with a specific patterns which basically generate void spaces. The control of this property is called infill percentage, and it is representative of the FDM technology. The infill is used to reduce the printing time, the used material and also to reduce the weight of the printed part.

High density region

Low density region

Even if there are more AM technologies, the most accessible is the FDM due to wide availability of 3D printers at low cost for personal use. In addition, the materials used for these machines are also affordable.

B.2 Limitations of the AM based on FDM for dynamically balanced bodies

The 3D printers for personal use are not generally designed to printing for a long time, therefore print a solid object can be prohibitive or impractical. Thus, printing time can be one of the most limiting factors, and it can be the main source of the printed part cost. The printing time depends on the size of the part, but for instance, the time required for printing one link of the four-bar linkage as a solid object was 24 h, which generate overheating and malfunction in a personal 3D printer. Thus, if a personal 3D printer based on FDM technology is used, the infill percentage should be lower than 100%.

Considering the optimized four-bar linkage, where the objective function is the compliance and it is constrained by the balancing conditions. If the links are printed with an infill percentage different from 100%, the optimized compliance will be affected, and the center of mass of the links cannot be preserved. Thus, a link of the four-bar linkage is studied with an infill of 40%