A fully coupled level set-based topology optimization of flexible components in multibody systems

A fully coupled level set-based topology optimization of flexible components in multibody systems is considered. Thereby, using the floating frame of reference approach, the flexible components are efficiently modeled and incorporated in multibody systems. An adjoint sensitivity analysis is utilized to obtain the gradient of the objective function with respect to a set of density-like design variables assigned to elements included in the underlying finite element model. The utilized adjoint sensitivity analysis provides a gradient, which is within numerical limits exact. In this process, the parametrization of material properties of finite elements has a significant influence on the calculated gradient, in particular for poorly filled elements. These influences are studied in detail. As an application example, the compliance minimization problem of a flexible piston rod in a transient slider-crank mechanism is considered. For this model, the influence of different parametrization methods on the obtained gradient is discussed, and a gradient strategy is proposed to overcome numerical issues included in different parametrization laws. Using this gradient strategy within a level set-based algorithm, a topology optimization of the flexible piston rod is performed. The corresponding results are then compared with optimization results provided by the method of moving asymptotes (MMA). Moreover, the computational effort of the sensitivity analysis is high and scales with the number of design variables. In this work, a gradient approximation is introduced using radial basis functions (RBFs). This helps to develop an appropriate gradient for a level set-based topology optimization of the flexible components in multibody systems, where the RBF-based design space reduction decreases the computational effort of the utilized sensitivity analysis. Finally, the efficiency gain obtained by the introduced design space reduction is demonstrated by optimization examples.


Introduction
Light-weight components are widely used in new generations of dynamic systems, which leads to a high material and energy efficiency. In general, the overall motion of such components includes vibrations and deformations. The method of flexible multibody systems is a well-known approach to incorporate the flexible body deformations into modeling and simulation of multibody dynamics, see (Shabana 1997;Wasfy and Noor 2003). In Fig. 1, a flexible multibody system with flexible and rigid bodies is shown. These components are connected to each other, as well as their environment by joints, springs, dampers, and bearings.
For modeling and simulation of flexible multibody systems, different approaches are available. Among others, formalisms based on nonlinear finite element (FE-) methods, such as the absolute nodal coordinate formulation (ANCF), and the floating frame of reference approach are common, see (Shabana 1997;Wasfy and Noor 2003). In case of nonlinear finite element methods, normally large systems of equations of motion with a vast number of degrees of freedom result, whose time integration is connected to a high computational effort. In contrast, assuming that the flexible body deformations remain small and elastic, the floating frame of reference approach provides a more compact formulation of the governing equations of motion, see (Schwertassek and Wallrapp 1999;Seifried 2013;Shabana 2013). Thereby, the body deformations are approximated as elastic coordinates multiplied by a set of global shape functions, which can, for instance, be obtained by a model order reduction from finite element models of flexible bodies. The accuracy and efficiency of the flexible multibody modeling and simulation can be maintained by a careful choice of the global shape functions, see for instance (Held et al. 2016).
The topology optimization of flexible multibody systems is a challenging task. This is among others, due to the complexity and nonlinearity of the equations of motion, as well as time, state and design dependency of dynamic loads. Hence, a gradient-based optimization process is normally required to develop meaningful optimized designs. The gradient calculation can be based on fully and weakly coupled strategies, see (Gufler et al. 2021;Tromme et al. 2018). The fully coupled sensitivity analyses provide a gradient, which is within numerical limits exact. Though, these methods are of a high complexity and computational effort, see (Bestle 1994;Bestle and Seybold 1992;Callejo et al. 2019;Dopico et al. 2018Dopico et al. , 2015Held et al. 2017;Nachbagauer et al. 2015;Sonneville and Brüls 2014) for direct and adjoint sensitivity analyses of multibody systems. On the contrary, a weakly coupled equivalent static load (ESL) method is proposed by Kang et al. (2005). Thereby, the dynamic behavior of the flexible multibody system is replaced by a set of equivalent static load cases. This reduces the required numerical steps, and thus, the computation costs of the gradient calculation. Nevertheless, the developed gradient is an approximation, which lacks the implicit design dependencies of chosen objective functions. For a comparison of topology optimization results of flexible multibody systems using fully and weakly coupled strategies, the interested readers can be referred to Moghadasi et al. (2018).
In literature, the optimization of flexible components in multibody systems is mostly performed by densityor level set-based methods. For instance, in Held et al. (2015Held et al. ( , 2016; Moghadasi et al. (2018) a solid isotropic material with penalization (SIMP) law is used to parametrize the underlying finite element model of flexible components in multibody systems. Fully and weakly coupled topology optimizations are then performed using the method of moving asymptotes (MMA). For details on SIMP approaches, interested readers can be referred to Bendsøe (1989); Du and Olhoff (2007); Olhoff and Du (2005); Pedersen (2000), and for MMA to Svanberg (1987). In contrast, Tromme et al. (2015) presents a fully coupled shape optimization, and Sun et al. (2017) presents a weakly coupled topology optimization of flexible multibody systems based on level set methods (LSMs), where the design boundaries are linked to the zero isosurface of a level set function defined over the design domain, see (van Dijk et al. 2013) for a comprehensive review of level set techniques. Moreover, for detailed comparisons of density-and level set-based optimization approaches, interested readers can be referred to the review paper by Sigmund and Maute (2013).
In the current work, a fully coupled topology optimization of flexible multibody systems is proposed, where the flexible bodies are modeled using the floating frame of reference approach. Thereby, similar to Held et al. (2015); Moghadasi et al. (2018), an adjoint sensitivity analysis is used to obtain the exact gradient of flexible bodies based on Bestle (1994); Held et al. (2017). However, in contrast to Held et al. (2015); Moghadasi et al. (2018), the optimization is performed by a modified level set method, see (Azari Nejat et al. 2022) for static optimization examples. Since in a level set-based optimization algorithm, the borderlines between interior and exterior of the design are clear and the gradient of empty elements in the exterior of the design has no physical meaning, the time-consuming adjoint sensitivity analysis should only be performed for elements in the interior of the design and around the free boundaries, which increases the efficiency of the fully coupled topology optimization of flexible multibody systems compared to the standard density-based cases, where the gradient is in each iteration obtained for all design elements. Moreover, a further design space reduction using radial basis functions (RBFs), see (Wang and Wang 2006a;Wei et al. 2018), is considered. This helps to gain even more efficiency by reducing the number of elements in the interior, for which the exact gradient is calculated by the adjoint sensitivity analysis. A gradient approximation using radial basis functions can in addition help to overcome numerical issues regarding the parametrization of material properties of poorly filled elements included in the underlying finite element model.
The remainder of this work is organized as follows: In Sect. 2, the required basics for the topology optimization of flexible multibody systems are briefly addressed and as rigid body bearing flexible body joints damper spring Fig. 1 Example of a flexible multibody system application example, the compliance minimization problem of a flexible piston rod in a transient slider-crank mechanism is presented. In Sect. 3, the numerical issues regarding different parametrization techniques are shown, and a strategy is proposed to overcome these issues. In Sect. 4, for the application example, the level set-based optimization results are presented and compared with the results from a standard SIMP optimization, which is performed by MMA.
In Sect. 5, the main steps of the implemented adjoint sensitivity analysis are briefly described. The computation times of relevant parts depending on the number of design variables of the flexible piston rod are mentioned. Radial basis functions are then presented for a gradient approximation from sampled gradient data. Subsequently, two sampling patterns are introduced. Based on approximated gradients, optimization examples are performed, and the quality of corresponding optimization results, as well as efficiency gains are compared. In the end, Sect. 6 concludes with a summary of this work.

Topology optimization of flexible multibody systems
In this section, the required basics for a fully coupled level set-based topology optimization of flexible multibody systems are briefly described, which are taken from Azari Nejat et al. (2022); Azari Nejat et al. (2019); Held et al. (2017); Moghadasi et al. (2018). First, the floating frame of reference approach for the incorporation of elastic deformations into the multibody models is addressed. After that, the optimization problem is formulated. The adjoint sensitivity analysis is summarized and the application example is presented. In the end, the utilized level set method for topology optimization is introduced.

Floating frame of reference formulation
In the floating frame of reference method, the deformation of a flexible body, which is assumed to be small and elastic, is described in a reference frame K R . The large nonlinear rigid body motion is captured as the relative movement between the reference frame K R and the inertial frame K I . This subdivision is depicted in Fig. 2. Using a global Ritz approach, the elastic displacement u P and rotation vector P of the coordinate system K P , which is rigidly attached to a point P on the flexible body, can be approximated at each time point t as (1) u P (R RP , t) ≈ P (R RP )q e (t), P (R RP , t) ≈ P (R RP )q e (t).
In this expression, the vector q e ∈ ℝ n e of n e elastic coordinates is time-dependent, whereas the matrix P ∈ ℝ 3×n e of n e global displacement shape functions and the matrix P ∈ ℝ 3×n e of n e global rotation shape functions depend only on the relative position R RP of the point P in the undeformed configuration. The shape functions can, for instance, be obtained by performing a model order reduction of the linear FE model of the flexible body, which reads as Thereby, M , D , and K are the mass, damping and stiffness matrix. Moreover, q ∈ ℝ n no is the vector of n no nodal coordinates and f represents the vector of applied loads. In order to approximate the elastic behavior of flexible bodies with a high accuracy, usually a fine discretization of the FE model is required, which leads to a great number of nodal coordinates. Different model reduction methods can then be used, see Held et al. (2016), to maintain the efficiency of the multibody simulation. In this work, the vector q of nodal coordinates is projected onto the vector q e of elastic coordinates by applying a modal truncation, where n e ≪ n no . It follows In this formulation, ∈ ℝ n no ×n e is the matrix of global shape functions, assembled from n e eigenmodes of the flexible body. The eigenmodes are obtained from the undamped eigenvalue problem of the FE model (2).
Considering the elastic coordinates q e and in addition to the rigid body coordinates as degrees of freedom, the equations of motion of a free flexible body can be formulated by applying, for instance, Jourdain's principle, see Schwertassek and Wallrapp (1999). As shown in Fig. 1, in multibody systems, different elements connect rigid and flexible bodies with each other and their environment and reduce, in this way, the number of degrees of freedom. Performing a coordinate partitioning similar to Held Floating frame of reference approach et al. (2017), the equations of motion of a flexible multibody system with f r redundant coordinates and f degrees of freedom can be derived in minimal form as ordinary differential equations (ODEs) with the following general formulation Thereby, y ∈ ℝ f and z ∈ ℝ f are the generalized position and velocity coordinates, whereas y r ∈ ℝ f r and z r ∈ ℝ f r are the redundant position and velocity coordinates. In the kinematic Eq. (4), Z is the kinematic matrix, which connects the first time derivative of generalized position coordinates y to the redundant velocity coordinates z r . Besides, in the kinetic Eq. (5), M , , and f represent the global mass matrix, the vector of generalized local accelerations, and the vector of generalized applied forces, respectively. In this work, the solver ode15s is utilized to solve the Eqs. (4),(5) of motion in MATLAB.

Optimization problem
In the current work, the compliance minimization problem of flexible bodies in multibody systems is considered. The optimization problem can be defined as Thereby, the vector x ∈ ℝ n contains the continuous densitylike design variables x i of the corresponding design elements from the discretized finite element model, i = 1(1)n . Each design variable x i is connected to the fill volume 0 ≤ v e i ≤ 1 of the design element i as x i = max(v e i , x min ) , where a small lower limit 0 < x min is considered to avoid numerical difficulties regarding ill-conditioned finite element system matrices. Here, x min is set to 0.01.
The objective function is the integral compliance of a flexible body over the simulation time [t 0 , t 1 ] , formulated as where q is the vector of nodal displacements and K the stiffness matrix of the flexible body.
The compliance minimization problem is subject to the Eqs. (4),(5) of motion and the volume constraint h = 0 , which imposes a limit of v max to the volume fraction v = 1 n ∑ n i=1 v e i of the flexible body. It should be noticed that (4) y = Z(y r )z r , (5) M(t, y r )̇z + (t, y r , z r ) = f (t, y r , z r ).
in a level set-based topology optimization, the volume v e i of a design element i depends on the level set function , which is defined over the design domain. Hence, in the optimization problem (6), x can be replaced by . Nevertheless, in this work, the gradient is calculated with respect to the design variables x i and then mapped from elements to nodes such that a nodal update of the level set function can be performed. This procedure is explained in Sect. 2.5.

Adjoint sensitivity analysis
For an exact gradient calculation of flexible multibody systems, direct and adjoint methods are available, see (Bestle 1994;Bestle and Seybold 1992;Callejo et al. 2019;Dopico et al. 2018Dopico et al. , 2015Held et al. 2017;Nachbagauer et al. 2015;Sonneville and Brüls 2014). The key part of these semianalytical approaches is variational calculus, which helps to expose both explicit and implicit dependencies of the selected objective function on the design variables x . In case of optimization problems with a large number of design variables, adjoint methods are normally more efficient. In the current work, an adjoint sensitivity analysis is utilized. Here, a brief introduction of this method is given, see (Held et al. 2017) for a detailed explanation. The variation of the integral objective function from Eq. (7) can be written as In this variational form, both dependent variations y of generalized position coordinates y and independent variations x of design variables x are included. The terms y T 1 and x T 1 represent the derivatives of the integrand F with respect to y and x . If all dependent variations over the entire time interval [t 0 , t 1 ] can be eliminated, the remaining coefficients of the independent variations x provide the desired gradient ∇ , see (Bestle 1994;Bestle and Seybold 1992). For this purpose, in the adjoint sensitivity analysis, additional zero terms are added to the varied objective function from Eq. (8). These zero terms result from the variation of kinematic and kinetic Eqs. (4),(5) multiplied by a set of time-dependent adjoint variables (t), (t) ∈ ℝ f and integrated over the simulation time [t 0 , t 1 ] , see Held et al. (2017). The required adjoint variables (t), (t) for elimination of dependent variations can then be obtained by solving the following adjoint system of ODEs which is a final value problem over the time interval [t 0 , t 1 ] with the final conditions (t 1 ) = 0, M(t 1 ) (t 1 ) = 0 . The terms y T 2 , y T 3 include the derivatives of the kinematic and kinetic Eqs. (4),(5) with respect to the generalized position coordinates y , and z T 2 , z T 3 their derivatives with respect to the generalized velocity coordinates z . The computation of y T 1 , y T 2 , y T 3 , z T 2 , z T 3 can, for instance, be performed using symbolic toolbox in Matlab.
For the numerical solution of the adjoint system (9), the generalized coordinates y and z are required at each time point. These state variables have to be determined by solving the Eqs. (4),(5) of motion in a pre-step. It is worth mentioning that in the current work, both the integration of the Eqs. (4),(5) of motion forward in time and the integration of the adjoint system (9) backward in time are performed using the solver ode15s in MATLAB.
With the obtained adjoint variables, the gradient ∇ can be computed by an evaluation of the following integral function where x T 2 and x T 3 represent the derivatives of the kinematic and kinetic Eqs. (4),(5) with respect to the design variables x . The terms x T 1 , x T 2 and x T 3 are prepared before evaluating Eq. (10). More details about the main numerical steps of the presented adjoint sensitivity analysis and the corresponding computational efforts are given in Sect. 5.1.
In the end, it should also be noticed that in the implemented level set-based topology optimization, the negative of the gradient ∇ is used to update the level set function . For this reason, the notion negative gradient is introduced as ∇ ∶= − ∇ and utilized in the following.

Application example
To demonstrate the functionality and capabilities of the methods presented in the current work, a planar flexible slider-crank mechanism is considered as application example, see Fig. 3, which is a well-known representative for the general class of flexible multibody systems with kinematic loops, see also Held et al. (2017Held et al. ( , 2016; Moghadasi et al. (2018); Sun et al. (2017); Tromme et al. (2015). This mechanism consists of a rigid crank with an effective length of 10cm , a flexible piston rod with dimensions (40 × 4 × 1) cm and a density of 0 = 8750kg/m 3 , as well as a piston with a weight of m = 0.1kg.
The flexible piston rod is discretized by 200 × 20 planar bilinear elements. To ensure an appropriate load transfer around the connections, the interface elements of the flexible piston rod are assumed as rigid, and its design domain consists of 198 × 20 inner elements with a reference Young's modulus of E 0 = 50GPa and a Poisson's ratio of = 0.3 . For a model order reduction of the flexible piston rod, 12 elastic coordinates are used, q e ∈ ℝ 12 . The body deformations are modeled as the multiplication of elastic coordinates q e and the matrix of global displacement shape functions. In this example, no rotation shape functions are necessary.
In the transient analysis, the system motion of the flexible slider-crank mechanism is considered within a time period of T = 2s . Thereby, the rigid crank starts at t 0 = 0s from the rest position with the angle (0) = 0rad , the angular velocity ̇(0) = 0rad/s , and is continuously accelerated with a prescribed acceleration ̈(t) until it reaches the angular position (2) = 12 rad and the angular velocity ̇(2) = 12 rad/s at t 1 = 2s . In Fig. 4, (t) , ̇(t) , and ̈(t) are shown. It is worth mentioning that the slider-crank mechanism is assumed to move horizontally. Hence, the gravity is not involved. Moreover, the inertia moment of the rigid crank has no influence on the dynamic behavior of the modeled system, and thus, it is not specified.
It should also be noticed that all implementations included in this work are carried out using Matlab, and the corresponding calculations are performed on a personal computer equipped with a 3.80 GHz Intel ® Xeon ® E3 -1270 v6 processor and 64 GB RAM. As a first example for the adjoint sensitivity analysis from Sect. 2.3, the negative gradient ∇ of the completely filled flexible piston rod, see Fig. 5, is calculated and shown in Fig. 6 over the design domain. Thereby, each node i on the surface plot corresponds to the gradient value ∇ i with respect to the i-th design variable.
The larger the negative gradient ∇ i , the greater is the role of element i in compliance minimization. The elements around the middle of longitudinal edges have the largest negative gradients. These elements are necessary to dispose of a huge bending of the flexible piston rod. In contrast, the smaller the negative gradient ∇ i , the lower is the contribution of element i to the compliance minimization. For instance, the elements in the middle of the piston rod do not significantly contribute to the structural stiffness. They rather cause bending due to their inertia. These elements possess the lowest negative gradients.
To obtain the presented negative gradient ∇ from Fig. 6, the underlying finite element model of the flexible piston rod is here parametrized using a SIMP law proposed by Olhoff (2005, 2007). More details about the parametrization laws are given in Sect. 3.

Modified level set method
In this work, a modified level set method is utilized, which is mainly adopted from the work proposed by Wei et al. (2018), and further developed in Azari Nejat et al. (2022) for the topology optimization of sparsely filled and slender structures. In the following, first the basics of this method are briefly described. After that, the practical formulation of the normal velocity V n for an appropriate development of the design during the optimization process and the update of the Lagrange multiplier associated with the volume constraint h from Eq. (6b) are addressed. In the end, an algorithmic representation of the optimization process is given.

Basics of the method
In a level set-based topology optimization, a level set function is defined over the design domain D. Depending on the level set function value at each arbitrary point of the design domain with the position vector r , this point can then be assigned to the interior ( (r) > 0 ), exterior ( (r) < 0 ) or free boundaries ( (r) = 0 ) of the structural design. In Fig. 7, this implicit representation is depicted by a 2D example, where the level set function is defined over a rectangular design domain D and creates due to its negative value in the middle region, a design with a central elliptical hole. For the discretization of the level set function in space, a typical sampling on a fixed Cartesian grid is considered. The nodal level set function variables { 1 , … , m } at the m grid nodes with position vectors {r 1 , … , r m } are then evolved in the optimization process, see also (Wang et al. 2003). In other words, the update of the continuous level set function is replaced by an update of the discretized level set function = [ 1 , … , m ] ⊤ . In the modified level set method, is updated by a forward Euler scheme as   where t k is the k-th pseudo-time point, t k+1 is the subsequent one, and t is the corresponding pseudo-time step size. Moreover, B is the augmented velocity vector, which is defined as In this formulation, the normal velocity V n at each node is multiplied by a distribution function ̂ . Depending on the value of , this term gives a smoothing factor between 1, for = 0 , and 0, for large level set function values such as | | ≥ 10 , see (Wang et al. 2003;Wei et al. 2018) for details.
The modification of the normal velocity by the distribution function ̂ hinders, among others, a limitless growth of the level set function in areas with a high value of normal velocity V n and, thus, helps to avoid the corresponding convergence issues.
It should be remarked that the update scheme (11) generally describes a steepest descent method with t as the step size and tB as the movement towards a close-tooptimal design.
To maintain the regularity of the discretized level set function in the optimization process, an approximate reinitialization procedure is used, see (Wei et al. 2018). Thereby, the discretized level set function is divided by the average of the spatial derivative values at selected nodes around the zero-level set. In this way, the spatial derivative value around the zero-level set is kept close to 1. This manipulation is useful and sufficient to retain the regularity of the level set function around the free boundaries. The approximate reinitialization can be considered either after each update or periodically. Here, the reinitialization is performed every 10-th iteration.
The utilized update scheme allows a hole creation in the interior of the design and, thus, is appropriate for a topology optimization. Besides, this scheme alleviates the stability limits included in conventional level set methods, which are, for instance, addressed in Luo et al. (2007); van Dijk et al. (2013). Nevertheless, to avoid issues such as oscillations and loss of the load path, additional adaptations are added to the presented update scheme. These algorithmic steps are adopted from Azari Nejat et al. (2022) and briefly illustrated in Sect. 2.5.4.

Formulation of normal velocity V n
In this work, the fixed uniform Cartesian grid of the level set function coincides with the underlying FE-mesh, see also (Wei et al. 2018(Wei et al. , 2021. Unlike a density-based topology optimization with element-wise design updates, here the nodal variables { 1 , … , m } of the discretized level set function are updated in the optimization process. Hence, the normal velocity V n is required on each node j = 1(1)m with the position vector r j . Consequently, the negative gradient ∇ should be mapped from elements to the nodes. Moreover, to keep the optimization process independent of absolute gradient values, the mapped gradient is here shifted and scaled. The adapted negative gradient ∇ ad j of each node j is constructed as Thereby, ∇ j is the average negative gradient of n su surrounding elements of the j-th node. This type of gradient mapping from elements to nodes has a filtering and smoothing effect, and helps to avoid checkerboard instabilities. However, this smoothing effect decreases with mesh refinement and is not meant to ensure the mesh-independency of optimization results (Sigmund and Maute 2013).
The term ∇ j is shifted by the minimum ∇ min,0 of average negative gradients of the initial design. The result is then divided by the difference between the maximum ∇ max,0 and minimum ∇ min,0 of average negative gradients of the initial design and multiplied by a scaling factor d. Here, d = 10 is chosen. In literature, also other types of gradient scaling, for instance, by the median or average of gradient values are addressed, see (Wei et al. 2018). Depending on the studied examples, the selected gradient scaling can influence the stability and convergence speed of the optimization process. Hence, the scaling strategy should be chosen carefully.
In each pseudo-time point t k , the normal velocity V n (r j ,t k ) of each node j is then formulated by the adapted negative gradient ∇ ad j and the current value k of the Lagrange multiplier as In this formulation, an appropriate selection of the Lagrange multiplier is of great importance. It has a decisive influence on the satisfaction of the volume constraint, the numerical stability, and convergence speed. In literature, the bisectioning algorithm (Sigmund 2001;Wang et al. 2007;Wang and Wang 2006b) and the augmented Lagrangian method (Belytschko et al. 2003;Luo et al. 2008;Wei et al. 2018;Wei and Wang 2006) are mostly used to update the Lagrange multiplier . Unlike bi-sectioning algorithm, in the augmented Lagrangian method, no abrupt correction of the volume fraction is considered, and the constraint violation is decreased gradually. Therefore, in the current work, an augmented Lagrange update scheme is utilized, which is described in the following.

Augmented Lagrange update scheme
Using augmented Lagrange update schemes developed for static optimization problems within the topology optimization of flexible multibody systems, numerical oscillations arise and the convergence is disturbed, see also (Sun et al. 2017). This issue is, among others, due to the design dependency of dynamic loads in flexible multibody systems. For this reason, here the Lagrange update scheme proposed by Wei et al. (2018) is taken and adapted such that the mentioned convergence issues can be avoided. The utilized update scheme consists of three stages: an early stage up to a selected number n E of iterations, a middle stage after that, and in the end a final stage as soon as the current volume fraction v k drops for first time under a certain limit v c = v max + greater but close to the volume limit v max . Thereby, is a prescribed tolerance. For the k-th pseudo-time point, it follows In the early and middle stage, k is selected equal to that adapted negative gradient ∇ ad j , which is greater than (1 − v max ) ⋅ 100% of all gradient values. This choice is based on the fact that, sooner or later, around (1 − v max ) ⋅ 100% of all nodes should be assigned with a normal velocity V n less than 0, such that in the end, a design with the volume fraction v = v max results. Moreover, the process parameter required in the middle and final stage is updated as where n E , , and max are given process parameters. Using the Lagrange update scheme (15) and (16), the optimization can start with an initial volume fraction v 0 far away from the prescribed volume limit v max . The designs in the early stage are developed under a widely relaxed volume constraint. Subsequently, in the middle stage, designs are more strongly moved into the direction of smaller constraint violations. In the final stage, further small volume changes are made such that v gradually converges to v max .
It is worth mentioning that the main difference between the current augmented Lagrange update scheme and the standard ones lies in the additional early and middle stage in Eq. (15). These stages make a meaningful reduction of large volume constraint violations in case of dynamic examples with design-dependent loads possible and help to avoid undesired oscillations of the volume fraction v in the corresponding optimization processes.

Optimization process
The algorithmic representation of the implemented optimization process can be seen in Fig. 8. In this optimization algorithm, first the initial vector 0 of the discretized level set function is defined, which captures the initial design. Based on that, the initial multibody simulation and sensitivity analysis are performed. Subsequently, the optimization loop starts with the update of the Lagrange multiplier . Using ∇ and , the normal velocity V n is constructed. In the next step, the discretized level set function is developed by the update scheme (11). For the resulting design, it is checked whether the load path is lost. In this case, it is reconstructed by an adaptation of Lagrange multiplier (adaptation I), see (Azari Nejat et al. 2022) for details.
If the load path is available, the corresponding design is then passed to the multibody simulation and sensitivity analysis. After that, to avoid a large increase of the compliance from one step to another, and to relax contaminating oscillations, a pseudo-time step size control based on the compliance history is considered (adaptation II), see again (Azari Nejat et al. 2022) for details. After all, the iteration number k is increased by 1 and the optimization loop is repeated until a prescribed iteration number n I is reached. The simple stopping criterion utilized in this work allows to observe the optimization histories and potential stability issues in a desired number of iterations. Nevertheless, in practical terms, it is meaningful to use an improved stopping criterion, which terminates the optimization process based on the history of compliance and volume fraction changes, see, for instance (Wei et al. 2018).
To perform an appropriate level set-based topology optimization by the presented algorithm, the introduced process parameters can, for instance, be set as t = 1 , n E = 10 , n E = 0.2 , = 0.2 , max = 10 , and = 0.02 . Nevertheless, this set of selected values is only a suggestion, and to tune the optimization properties, such as the convergence behavior, the interested readers shall change the process parameters as desired.

Parametrization laws and their influence on the gradient
In topology optimization, the design domain D is typically discretized by a finite element model. Thereby, the parametrization of the material properties of the finite elements is a crucial task, which has a significant influence on the quality of the model order reduction, and thus, on the calculated gradient and corresponding optimization results. Different parametrization laws have been developed and used for structural and topology optimization of eigenvalue, vibration, and transient problems, see, for instance (Allaire and Jouve 2005;Du and Olhoff 2007;Held et al. 2016;Olhoff and Du 2005;Pedersen 2000).
To show the influence of the parametrization on the gradient of the studied transient slider-crank mechanism from Fig. 3, in the following, two approaches, namely a linear parametrization, as well as the solid isotropic material with penalization (SIMP) law proposed by Olhoff and Du (2005) are considered. These methods are utilized to parametrize the underlying finite element model of an example for the flexible piston rod with a central elliptical hole, see Fig. 9. For this example, the numerical artifacts included in the calculated gradients are demonstrated. A strategy is then proposed, which helps to avoid issues regarding the applied SIMP parametrization. This strategy provides an adapted gradient, which can be used within a level set-based optimization process.
Following the level set-based structural optimization, the density i and Young's modulus E i of each design element i can be linearly penalized by the corresponding design variable x i as Thereby, 0 is the density and E 0 is the Young's modulus of the chosen reference material. As mentioned before, to avoid ill-conditioned finite element system matrices, a lower limit x min = 0.01 is introduced for all design variables x i . For completely filled elements, it holds x i = 1 . Due to the high ratio of element density i to the Young's modulus E i in this parametrization, the low density areas with x min ≤ x i < 0.1 affect the model order reduction of the corresponding finite element model and lead to spurious modes, see also (Allaire and Jouve 2005;Pedersen 2000).
To indicate the influence of spurious modes on the gradient calculation, the gradient of the flexible piston rod from Fig. 9 is computed. Thereby, the underlying FE model is parametrized by the linear approach (17). In Fig. 10a, the negative gradient ∇ over the design domain is shown. It turns out that the elements with x i = x min , which form the hole in the middle of the design domain, are assigned with much higher negative gradient values than the filled elements. This is due to the spurious modes and physically not meaningful. This negative gradient is not appropriate for a level set-based topology optimization. Hence, an adapted parametrization concept is required to dispose of mentioned spurious modes and their undesired impact on the gradient.
In Held et al. (2016), the influence of different SIMP laws on the eigenmodes of a test structure is studied. The analyses unveil that the adapted SIMP approach by Olhoff (2005, 2007), defined as is suited to avoid the occurrence of spurious modes in the model order reduction. Thereby, the density of elements with a design variable x i between 0.1 and 1 is penalized linearly, whereas the density of poorly filled elements with a material amount between x min = 0.01 and 0.1 is penalized by a factor c and the power-law x q i . The power q is typically set to 6, and the factor c is selected such that the density penalization around x i = 0.1 is continuous, here c = 10 5 is chosen. Besides, the Young's modulus E i is penalized by the powerlaw x p i , where the power p is typically set to 3, see also (Du and Olhoff 2007) for the set of chosen parameters. In this parametrization approach, the exponentiated penalization of the Young's modulus helps to develop thin connections and struts with a higher stability. Moreover, choosing q = 6 greyscale element-wise representation y x 40 cm 4 cm Fig. 9 Design of the flexible piston rod with a central hole and p = 3 , in the regions with x i = x min , the penalization of the density i compared to the Young's modulus E i is clearly stronger. In this way, the mentioned spurious modes can be avoided. It should be noticed that in literature, also other parametrization methods for dynamic problems are available, which can be used instead of the utilized SIMP approach by Olhoff (2005, 2007), see, for instance, the rational approximation of material properties (RAMP) method proposed by Stolpe and Svanberg (2001).
In the fully coupled adjoint sensitivity analysis, the partial derivatives of the density i and Young's modulus E i with respect to the design variables are included. Based on the SIMP parametrization (18), these partial derivatives can be formulated as Nevertheless, using approach (18) with the partial derivatives from Eq. (19), the obtained gradient includes some irregularities, which are due to the insignificantly small value of i x i in areas with x i = x min and its discontinuity at x i = 0.1 . To indicate these issues clearly, the FE model of the flexible piston rod from Fig. 9 is parametrized using approach (18) and the adjoint sensitivity analysis is repeated. The corresponding negative gradient ∇ over the design domain is shown in Fig. 10b. Here, using a modified ratio of element density i to the Young's modulus E i for elements in the hole area with x i = x min , the spurious modes are avoided, and the negative gradient values of these void elements are insignificantly small and around 0, which is in general meaningful.
However, due to the difference between partial derivatives i x i with respect to x min ≤ x i < 0.1 and 0.1 ≤ x i ≤ 1 , a jump of negative gradient values between the solid part of the design and the hole area can be observed. Besides, some huge peaks are included in the negative gradient ∇ from Fig. 10b. These are caused by the discontinuity of the partial derivative i Thus, considering the mentioned issues, also the negative gradient from Fig. 10b can not be used within a level setbased optimization process. Therefore, an alternative strategy for gradient calculation is considered, where the finite element model is parametrized as before by the SIMP law (18). However, the adjoint sensitivity analysis is used only for gradient calculation of design variables 0.1 ≤ x i ≤ 1 . For the design variables x min < x i < 0.1 , the gradient is approximated using radial basis functions such as multiquadratic (MQ) splines. In Sect. 5.2, this approximation method is presented. Moreover, the negative gradients of void elements with x i = x min have no physical meaning. These are set to a minimum ∇ min , which is selected less than or equal to the minimum of negative gradients of all design variables x min < x i ≤ 1 . Hence, the gradient computation can be described by the following three cases Since, in this strategy, the gradient of poorly filled elements with x min ≤ x i < 0.1 are not obtained by the adjoint sensitivity analysis, numerical artifacts are less likely to occur. Besides, a higher efficiency can be expected than in the case, where the gradient of all elements should be obtained.
To demonstrate the benefits of gradient strategy (21), this procedure is utilized to obtain the gradient of the flexible piston rod from Fig. 9, whereby the flexible piston rod is as before parametrized using approach (18). The obtained negative gradient ∇ over the design domain is shown in Fig. 10c. Compared to the negative gradient from Fig. 10b, the numerical peaks are avoided, and the void elements are provided with the lowest gradient value ∇ min . These changes make the negative gradient ∇ appropriate for a level set-based topology optimization of the flexible piston rod.

Level set-based optimization example and validation
Using the presented level set algorithm, a topology optimization of the flexible piston rod is performed (example 1-A). In addition, for comparison purposes, also a standard SIMP optimization is considered, where the method of moving asymptotes (MMA) is applied to perform the topology optimization (example 1-B methods within a static optimization example can be found in Azari Nejat et al. (2022). It is worth mentioning that in both examples, the underlying FE model of the flexible piston rod is parametrized by the SIMP approach from Eq. (18). In this way, a comparison of results is possible. Nevertheless, as presented in the gradient strategy (21), in the applied level set algorithm, the gradient of elements with a density lower than 0.1 is approximated, whereas in the standard SIMP optimization, this distinction of cases is not included. In each iteration, the gradient with respect to all design variables is calculated by the adjoint sensitivity analysis and directly passed to the MMA algorithm.
In both examples, the volume limit v max is set to 0.4, and the design is developed within a chosen number n I = 100 of optimization iterations. However, if desired, v max and n I can be freely changed. The level set-based optimization in example 1-A starts with a simple initial level set function = 1⋅10 −3 over the entire design domain of the flexible piston rod, which leads to a completely filled initial design with a volume fraction v 0 = 1 . It should be noticed that the completely filled initial design for different level set-based optimizations in this work is the most general case without any preconditioning. Though, there are no limitations on the choice of initial designs with a volume fraction less than 1. Moreover, the SIMP optimization in example 1-B starts from a homogeneous grey initial design with the volume fraction v 0 = v max . Fig. 11 shows the final designs of the examples 1-A and 1-B including their compliance end , volume fraction v end , and normalized greyness ḡ end . As it is noticed in Sigmund and Maute (2013), the qualitative comparison of smooth designs represented by zero-level contours with densitybased greyscale designs is not meaningful. For this reason, in Fig. 11, the final design of example 1-A is represented both as contour plot of the corresponding level set function, created by the function contourf in Matlab, and in addition as density-based greyscale plot. The latter is better comparable to the density-based greyscale plot of the optimization result obtained by the standard SIMP optimization.
The problem definitions and optimization results of the mentioned examples are summarized in Tab. 1. Besides, the corresponding optimization histories are shown in Fig. 12.
The optimization results of both examples are in a good agreement. In both cases, a final design with a compliance end = 0.049 Nmm and a volume fraction v end = 0.400 is reached. Hence, it can be stated that the presented level set algorithm leads to a design with valid quantitative values. Nevertheless, the final designs of both procedures possess clearly different complexities and greynesses. This can be traced back, among others, to the considerable differences of the optimization algorithms, including process parameters and selected initial designs. The optimization histories in As shown in Fig. 12a, due to the chosen grey initial design in example 1-B, the first designs in the corresponding SIMP optimization possess larger compliance values than the first designs in the level set method. Even so, the optimization process in both examples converges to similar compliance values. Moreover, the development of the volume fraction v in example 1-A is different from example 1-B, see Fig. 12b. In example 1-A, the augmented Lagrangian method is used. This method allows to start with a volume fraction v 0 = 1 far away from the volume limit v max and drives the optimization process towards v max . On the contrary, the utilized MMA in example 1-B includes a Lagrange multiplier, which keeps the volume fraction v in the whole optimization process less than or equal to the volume limit v max . Nevertheless, despite different implementations of the volume constraint, here both procedures converge to similar volume fractions.
In the studied examples, the compliance and volume fraction of the standard SIMP optimization converge faster than in case of the level set-based optimization. However, a closeto-optimal design with a high greyness is not useful. Hence, also the convergence of the greyness in both optimizations should be compared.
It is known that the greyness in a level set-based topology optimization is limited to elements intersected by the zero-level set, whereas the SIMP method normally leads to a higher number of grey intermediate elements. This difference can be qualitatively seen in the final designs of examples 1-A and 1-B in Fig. 11, which speaks in favor of less grey final design of example 1-A using the level set algorithm. Besides, for a quantitative comparison, the normalized greyness ḡ can be defined as  see also (Sigmund 2007). In Fig. 12c, the greyness histories of examples 1-A and 1-B are shown. The convergence of the greyness in the studied SIMP optimization is clearly slower than in the level set-based optimization. In example 1-A, the greyness is in the beginning equal to 0, increases slightly with the creation of new boundaries, and reaches a final value ḡ end = 0.12 in the end of the optimization process. On the other hand, in example 1-B, the SIMP optimization starts with a high greyness close to 1. The greyness decreases during the optimization process and reaches a final value ḡ end = 0.21 in the end, which is about 75% higher than the greyness of the final design obtained with the level set algorithm.
The mentioned comparisons with the well-known SIMP optimization show that the proposed level set algorithm provides valid optimization results for the compliance minimization problem of the studied flexible piston rod. Moreover, the possible efficiency gain using the presented level set algorithm is discussed in Sect. 5.

Computational effort and reduction of design space
Compared to weakly coupled topology optimization of flexible multibody systems (Kang et al. 2005), fully coupled topology optimization is computationally more expensive due to the tremendous effort in the adjoint sensitivity analysis, see (Moghadasi et al. 2018;Tromme et al. 2018). Hence, one way to improve the numerical efficiency of fully coupled optimization problems is to reduce the number of design variables, which is typically large in topology optimization. In this section, the main steps of the adjoint sensitivity analysis are briefly reviewed, and the computational effort of the time-consuming steps is shown by test examples. A gradient approximation is presented based on radial basis functions and applied for the gradient calculation of the flexible piston rod. Thereby, different patterns for design space reduction are introduced and compared. In the end, optimization examples are included, which show the efficiency gain and the validity of the obtained results.

Main steps of the adjoint sensitivity analysis
The adjoint sensitivity analysis in the current work consists of a number of computational steps, which are illustrated in Fig. 13 and described in detail in Held et al. (2017). For a given flexible multibody system, the FE models of the flexible bodies are prepared in a preprocessing step. Then, the corresponding shape functions with = 1(1)n e are obtained by modal truncation. In the next step, a number n SID of system matrices V a with a = 1(1)n SID , the so-called standard input data (SID), see Wallrapp (1993), are computed, which contain all necessary information on the flexible body for the multibody simulation. The SID are object-oriented and structured such that the system equations can be evaluated with high efficiency, see (Schwertassek and Wallrapp 1999). Solving the Eqs. (4), (5) of motion forward in time, the generalized coordinates y and z are determined. Using the resulting elastic coordinates q e , the integral compliance of the flexible bodies is then evaluated.
To efficiently calculate the derivatives d i of the shape functions with respect to the design variables x i , Nelson's method (Nelson 1976) is used. After that, the derivatives d i V a of the SID with respect to x i can be determined and gathered in the so-called augmented standard input data (ASID) as suggested in Held (2022), which are later required for the evaluation of the gradient Eq. (10). In the next step, the adjoint variables are obtained by solving the adjoint system (9) backward in time. Finally, the sought negative gradient ∇ can be determined by evaluating the gradient Eq. (10) using the previously determined generalized coordinates y and z , the adjoint variables and , as well as the derivatives d i and ASID. To indicate the connection between the number n of design variables and the computational effort of the utilized  Fig. 5 with different FE-discretizations ranging from 10 × 100 to 38 × 380 is performed. In Fig. 14, the computation times t c of the program steps with a significant computational effort are displayed over the number n of design variables x i .
Using modal truncation for the model order reduction, the computational effort of the multibody simulation (forward integration) and the solution of the adjoint system (backward integration) are independent of the number n of design variables, see Fig. 14a. The computation times of these steps scale with the number n e of selected elastic coordinates q e . Here, 12 elastic coordinates are used, see also the computation times of the adjoint sensitivity analyses with different numbers n e of elastic coordinates included in Azari Nejat et al. (2020).
On the contrary, the computation times for the calculation of the derivatives of the global shape functions, the ASID, and for the evaluation of the gradient Eq. (10) rise with an increased number of design variables as shown in Fig. 14b. While the computation times for the derivatives of the global shape functions and the ASID increase exponentially, the effort to evaluate the gradient Eq. (10) rises linearly, compare also (Moghadasi 2019).

Gradient approximation using radial basis functions
To increase the efficiency of the presented gradient calculation, a design space reduction is of great interest. Here, the aim is to reduce the number of design variables independently of the underlying FE-mesh. To this end, the n design variables x of elements in design domain D are subdivided into two subsets, namely Ŝ = {x 1 , … ,x̂n} and Š = {x 1 , … ,x̌n} , where n =n +ň . The adjoint sensitivity analysis is then performed to obtain the exact value of the negative gradient ∇ only for the n design variables included in subset Ŝ . Based on the obtained gradient values, the negative gradient ∇ with respect to remaining ň design variables included in subset Š is then approximated. For this purpose, radial basis functions (RBFs) are utilized. In the following, this method is briefly described, see (Wang and Wang 2006a;Wang and Kang 2018;Wei et al. 2018) for more details. Radial basis functions are radially symmetric functions positioned at desired points in the space, which can be used for a scattered data interpolation (Morse et al. 2001). In literature, radial basis functions are, among others, applied to reconstruct smooth, manifold surfaces from scattered surface points, and to fill incomplete meshes (Carr et al. 2001). The unique solvability of the interpolation problem, as well as the smoothness and convergence of interpolants, are some of interesting properties of a RBF-based interpolation (Wang and Wang 2006a). There exist various types of radial basis functions, e.g., thin-plate spline, exponential spline, multiquadratic (MQ) spline, as well as Gaussian and compactly supported radial basis functions, see also (Wang and Wang 2006a). In the current work, the well-established multiquadratic spline is used for the gradient approximation.
At an arbitrary point of the 2D design domain D with the position vector r = [r x , r y ] ⊤ , the value g 0 (r) of a multiquadratic spline positioned at a point with the position vector r 0 = [r x,0 , r y,0 ] ⊤ can be calculated as Thereby, || ⋅ || is the Euclidean norm and a shape parameter. In Fig. 15, the shape of a multiquadratic spline positioned at r 0 = [0, 0] ⊤ is shown over a 2D design domain, whereby is set to 10 −4 . Increasing , flatter shapes result, which theoretically lead to more accurate approximations. However, in practical terms, the gained accuracy of the interpolation by an increase of is limited, see (Cheng et al. 2003).
In the considered gradient approximation, radial basis functions are placed at the center of design elements linked to subset Ŝ , with the position vectors r i = [r x,i ,r y,i ] ⊤ , where i = 1(1)n . The negative gradient ∇ (r) at an arbitrary point with the position vector r = [r x , r y ] ⊤ can then be parametrized as Thereby, i is the expansion coefficient of the i-th radial basis function and g i (r) the value of the i-th radial basis function at the point with position vector r . Moreover, P is a first-degree polynomial with coefficients p 0 , p 1 and p 2 , which accounts for the constant and linear parts of ∇ and guarantees the positive-definiteness of the solution; nevertheless, P is not necessary for all types of radial basis functions (Morse et al. 2001).
To guarantee the unique solution of the RBF-based interpolation, the following orthogonality or side conditions are imposed upon the expansion coefficients, see (Carr et al. 2001;Morse et al. 2001). Considering Eq. (24) at the center of elements linked to the subset Ŝ and the side conditions (25), a system of n + 3 linear equations follows, which can be written in matrix formulation as Thereby, = [ 1 …̂n p 0 p 1 p 2 ] ⊤ is the vector of generalized expansion coefficients and G is the collocation matrix gathered as As mentioned in (Kansa et al. 2004), the collocation matrix G is invertible. Solving the system of linear Eqs. (26), the generalized expansion coefficients can be obtained. After that, the unknown negative gradient ∇ with respect to design variables from subset Š can be approximated by evaluating Eq. (24) at the center of corresponding elements.
Thereby, the accuracy of the approximation depends, among others, on the number n and position of the design elements, for which the exact value of the negative gradient ∇ is obtained from the adjoint sensitivity analysis. In the following, two different sampling patterns are utilized for the RBF-based gradient approximation.

Sampling patterns
To indicate the efficiency gain achieved by the described design space reduction, the topology optimization of the flexible piston rod from example 1-A is repeated with two sampling patterns for the RBF-based approximation of the negative gradient ∇ . In Fig. 16a, pattern 1 for the distribution of the chosen design elements over the completely filled design domain is shown, where these sample elements are colored blue. Thereby, from the original regular finite element mesh 100 columns, and in each column, 10 elements are selected such that a reduced design space with n = 1000 instead of n = 3960 design variables results for the sensitivity analysis.
It should be noticed that the chosen elements are not uniformly positioned. Care is only taken that the element distribution is symmetrical with respect to the longitudinal centerline to sustain the symmetry of the evolving design in the optimization process. Besides, the first and last columns of the design elements are completely included in the sampling pattern. In this way, the negative gradient ∇ at the left end and at the right end of the design domain next to the rigid interface elements is obtained with an appropriate accuracy.
In Fig. 16b, pattern 2 for a further element distribution over the completely filled design domain is shown. Thereby 51 columns, and in each column, 6 elements are selected such that a reduced design space with n = 306 elements results. Reducing the number n of the chosen design elements from 1000 in pattern 1 to 306 in pattern 2, the quality of the gradient approximation decreases. For instance, approximating the negative gradient ∇ over the completely filled design domain by multiquadratic splines with = 10 −4 and gradient values of the sample elements in patterns 1 and 2, the approximation results can be compared with the corresponding exact values from Fig. 6. In Fig. 17, a histogram of the absolute errors for the gradient approximation using pattern 1 and 2 over the completely filled design domain is given. It can be seen that the absolute errors using pattern 2 are larger and distributed in a wider range than for pattern 1. As mentioned before, larger shape parameters can theoretically help to decrease the gradient approximation error. Nevertheless, higher values for can also lead to the loss of stability, see (Cheng et al. 2003). Hence, in the following, only = 10 −4 is considered for both pattern 1 and 2.
Moreover, to develop the free boundaries of the design with an appropriate accuracy, it is decided not to reduce the design elements in the boundary areas with a design variable 0.1 ≤ x i < 1 . Hence, in principle, the update of the level set function over the solid areas and, with that, the hole creation are in the following examples based on the presented gradient approximation, whereas the created boundaries are then developed further based on exact gradient values. For instance, in Fig. 16c, the chosen 450 design elements of the flexible piston rod with an elliptical hole in the center are shown, which consist of blue-colored design elements over the solid area as in pattern 2 and, in addition, red-colored elements in boundary areas.

Topology optimization using approximated gradients
In the following, the RBF-based gradient approximation is employed in the topology optimization of the flexible piston rod. In Fig. 18, the optimization result from example 1-A defined in Tab. 1 is shown together with the results using pattern 1 (example 2-A) and pattern 2 (example 2-B) for gradient approximation in the optimization process. and with that, the success of gradient approximation using pattern 1 and 2 in the optimization process.
Nevertheless, it should be noticed that the gradient approximation is here only tested for the compliance minimization problem of the flexible piston rod. Considering other problem types or different application examples, other approximation accuracies, and thus, sampling patterns might be required. The presented examples 2-A and 2-B are only meant to show the potential of design space reduction using radial basis functions.
In Fig. 19, the computation times t c of the calculation of global shape function derivatives with respect to the design variables, ASID terms, and the evaluation of the gradient Eq. (10) for the examples 1-A, 1-B, 2-A, and 2-B, as well as the corresponding number n of considered design variables for the exact gradient calculation are shown over the iteration number k.
The computation time plots in Fig. 19 are jagged, and some of them include also jumps from one iteration to the next. Though, these fluctuations of computation times are in complex numerical processes not uncommon. Apart from that, the presented time courses scale in general with the number n of considered design variables in each iteration k. Unlike in the SIMP optimization example 1-B, where in each iteration, the exact gradient is calculated with respect to all design variables, in the level set-based optimization examples 1-A, 2-A, and 2-B, the number n of considered design variables for the exact gradient calculation decreases, and with that, the computation times of the three mentioned program steps. Hence, in the following, the computation times of example 1-B are taken as reference computation times.
The sum of computation times for the mentioned program steps is listed in Tab. 2. These take in example 1-B 281 min and in example 1-A 144 min, which is an efficiency gain of 49%. Reducing the number of design variables over the solid structure in example 2-A and 2-B, the numerical effort of mentioned program steps is further decreased, and the sum of corresponding computation times drops to 75 min in example 2-A and 50 min in example 2-B. In other words, 73% and 82% less than the sum of mentioned computation times in example 1-B. The presented efficiency gains might be, in particular, in case of large-scale topology optimizations, a worthwhile contribution.
It should be noticed that in Tab. 2, the entire computation times of the optimization process in different examples are listed as well. These include also the computational efforts of other program steps such as multibody simulation and solution of the adjoint sensitivity analysis, which do not depend on the number n of chosen design variables. Hence, a comparison of entire computation times in different examples does not directly reflect the efficiency gain obtained by the presented design space reduction.
In the end, it can be stated that the results discussed in this section show the success and time savings of gradient approximation using pattern 1 and 2 in the level set-based topology optimization of flexible piston rod. Nevertheless, as mentioned before, to develop the free boundaries with an appropriate accuracy, the exact gradient for all elements in boundary areas with a design variable 0.1 ≤ x i < 1 is provided in each iteration k.
In Fig. 20, the ratio n/n of design elements used for exact gradient calculation in example 2-B is shown over the iteration number k, whereby it is distinguished between the design elements distributed in the solid areas and around the boundaries. Starting the optimization process with a completely filled design domain, no free boundaries exist in the beginning. With the creation of new holes, the number of design elements chosen in the solid areas decreases, since the void elements in the hole areas are excluded. On the contrary, new boundaries increase the number of design elements considered in these areas, and it can be seen that in example 2-B, the majority of selected design elements after k = 15 are around the boundaries.
For sure, also an appropriate reduction of design elements around the boundaries can be considered. It is even possible to replace the topology optimization of the moving boundaries by an efficient shape optimization, see for instance (Tromme et al. 2015), whereas coarse sampling patterns such as pattern 1 and 2 are further used for a hole creation. These ideas are not included here, and require comprehensive investigations in a separate work.

Summary and conclusion
In the current work, the fully coupled level set-based topology optimization of flexible bodies in multibody systems modeled with floating frame of reference formulation is addressed. The deformation of the flexible bodies is described by global shape functions, which are obtained from finite element models using model order reduction. The adjoint sensitivity analysis is utilized for an exact gradient calculation, and the optimization is performed using a modified level set method. The presented approach is studied for the compliance minimization of a flexible piston rod within a slider-crank mechanism. At first, the significant influence of parameterizing the material properties of the finite element model on the calculated gradient is shown. In particular, the numerical issues of the gradient calculation with a linear and traditional SIMP parametrization are discussed and shown for the example of a flexible piston rod. In case of a linear penalization, which is typically used in structural level set-based topology optimization, the high ratio of the element density to stiffness leads to spurious modes in void areas of the design domain and distorts in this way the gradient calculation. Using the SIMP parametrization method proposed by Olhoff and Du (2005), these artificial modes can be avoided. Nevertheless, due to the inconsistency between negative gradient values of the solid parts and the hole areas of the design in this approach, the resulting gradient can not be directly used within a level set-based optimization process. Hence, a gradient strategy is included, which replaces the gradient of design elements with a low fill volume by a RBFbased gradient approximation.
Using the proposed gradient strategy, a level set-based topology optimization of the flexible piston rod is performed. Besides, to test the quality of the corresponding results, also a standard SIMP optimization of the flexible piston rod is considered, where the optimization is performed by MMA. Despite clear differences in optimization histories and shapes of final designs, in both cases, similar compliance values are reached. This shows the appropriateness of the optimization result obtained by the level set algorithm.
The adjoint sensitivity analysis is a time-consuming process. To lower the computational effort of program parts, which depend on the number of design variables, the RBFbased gradient approximation is considered. This approach allows a design space reduction independently of the underlying FE-mesh. The topology optimization of the flexible piston rod is performed with two different sampling patterns for the gradient approximation. The validity of corresponding results and the efficiency gains confirm the success of the considered design space reduction.