On design sensitivities in the structural analysis and optimization of flexible multibody systems

The structural analysis and optimization of flexible multibody systems become more and more popular due to the ability to efficiently compute gradients using sophisticated approaches such as the adjoint variable method and the adoption of powerful methods from static structural optimization. To drive the improvement of the optimization process, this work addresses the computation of design sensitivities for multibody systems with arbitrarily parameterized rigid and flexible bodies that are modeled using the floating frame of reference formulation. It is shown that it is useful to augment the body describing standard input data files by their design derivatives. In this way, a clear separation can be achieved between the body modeling and parameterization and the system simulation and analysis.


Introduction
Analysis and optimization of flexible multibody systems are important steps in the computed-aided design and dimensioning process of dynamic mechanisms. In [14,20], typical examples from optimal control and structural optimization of flexible multibody systems are given and reviewed. In the optimization, gradient-based optimization algorithms usually work most efficiently, but they require sensitivity information of the objective and constraint functions [9].
One way to determine these gradients is by numerical differentiation. It is a simple and easy-to-implement method, but it suffers from various deficiencies. For example, the gradient can only be approximated, the perturbations of the design variables are not known a priori, and the computational effort increases proportionally with the number of design variables. However, in structural optimization, which is the focus of this work, the number of design variables is often large. For instance, the topology optimization of a flexible slidercrank mechanism presented in [13]   Due to these disadvantages, semianalytical methods, such as direct differentiation and the adjoint variable method, are nowadays often preferred over numerical differentiation, despite the fact that their derivation and implementation effort is usually more demanding and significantly higher. The key idea of semianalytical methods is to derive a set of additional sensitivity equations from which the gradient can be computed. Thereby the structure of the sensitivity equations depends on the formulation of the system equations and the type of the criterion function defined.
Detailed derivations of the direct differentiation method and the adjoint variable method for the sensitivity analysis of rigid multibody systems are presented in [3][4][5]12]. The methods have also been applied to flexible multibody systems. In [7], structural sensitivity analysis is performed with the direct differentiation method for a flexible slider-crank mechanism and a vehicle chassis. For the description of the flexible bodies, the floating frame of reference approach is used, whereby the body deformations are approximated with a set of eigenmodes. The application of the direct differentiation method to flexible multibody systems modeled with nonlinear finite element methods is shown, for example, in [6], where topology optimization of truss structures assembled from beams is performed. Also, both the direct differentiation and the adjoint variable method are employed for sensitivity analysis of beam structures modeled with the absolute nodal coordinate formulation; see [16].
The current work focuses on flexible multibody systems modeled with the floating frame of reference approach. It is further assumed that the material properties and the geometry of the bodies are arbitrarily parameterized by the design variables p ∈ R h , which is typically the case in structural optimization [11].
The general structure of the considered multibody system and the notation used are given in Fig. 1. It is assembled from rigid and flexible bodies, spring and damper elements, and actuators. The bodies are connected with each other and the foundation by ideal joints and bearings. Formulating the system equations in implicit descriptor form yields The first two equations (1a) and (1b) are the initial conditions at position and velocity level. They are followed by the kinematic equations (1c) and the kinetic equations (1d). The latter are often denoted as equations of motion. The system equations are completed by the constraint equations (1e).
To analyze the dynamic system (1) in the time domain t ∈ [t 0 , t 1 ] a scalar criterion function ψ ∈ R has to be defined first. If the transient system behavior is of interest, then ψ is typically an integral-type function of the form If the adjoint variable method is employed to compute the gradient ∇ψ , then two steps are necessary. At first, a system of adjoint differential equations has to be derived. Then the gradient can be calculated using the results of the transient analysis and the adjoint analysis. For criterion functions of the form (2) and system equations of the form (1), the derivation of the adjoint system and the gradient equation is comprehensively described, for instance, in [1].
For system (1), the adjoint variables μ ∈ R f and ν ∈ R f can be obtained from the solution of the differential equation whereby ξ ∈ R f and γ ∈ R nc are auxiliary variables. They are necessary to take into account the accelerations and constraint equations in the adjoint analysis and can be determined from the linear equation After solving the adjoint system of equations (3), only the adjoint variables η 0 ∈ R f and ζ 0 ∈ R f at the initial time t 0 remain to be determined: A. Held Then the gradient can finally be computed solving It can be observed from the Eqs. (3)-(6) that two types of derivatives of the system equations are required in the sensitivity analysis. On the one hand, the system equations must be differentiated with respect to the position and velocity variables and, on the other hand, with respect to the design variables. The former derivatives can be determined systematically for the different multibody system formulations. The latter, which are denoted as design sensitivities throughout the paper, are, however, strongly problem-specific. Therefore setting up Eq. (6) is an individual and time-consuming process.
In the literature, symbolic algebra systems are often used to derive the sensitivity equations. For instance, the numerical examples in [4,7] are created with MAPLE. Alternatively, it is recommended in [8] to apply automatic differentiation [10] to obtain the derivatives of the system equations with respect to the state and design variables, in particular, for complex multibody systems. A third way is presented in [22], where the state sensitivities are computed analytically and the design sensitivities by numerical differentiation.
Aside from the numerical advantages and disadvantages of the three ways, they all suffer from mixing up two fundamental problems: the derivative of body-specific properties, such as the center of gravity, with respect to the design variables, and the derivative of the system equations with respect to body-specific properties. To improve structural optimization procedures, this work aims to remove this mixing. Therefore the kinematic and kinetic equations are systematically differentiated with respect to arbitrary geometry and material parameters. Further, it is recommended to adapt the body describing object-oriented standard input data [21] and augment them by their design derivatives. In this way the modeling and parameterization of the bodies can be considered separately from the dynamic simulation and adjoint sensitivity analysis.
The rest of this paper is organized in the following way. Section 2 addresses the partial derivatives of the system equations with respect to the design variables. After a brief review of the body kinematics and kinetics, their design sensitivities are systematically derived. In Sect. 3 the standard input data are summarized, and an augmentation of the object-oriented data set is suggested to allow a general and systematic parameterization of rigid and flexible bodies. The procedure is then demonstrated in Sect. 4 using an example from the topology optimization of a flexible multibody system. Finally, Sect. 5 concludes with a brief summary and discussion.

Design sensitivities of system equations
To evaluate the gradient equation (6), the partial derivatives of the implicit initial conditions, of the implicit kinetic equations, and of the constraints at acceleration level with respect to the design variables are required. In the following section, these derivatives are systematically derived. At first, the marker kinematics and relative marker kinematics of flexible bodies are briefly summarized and differentiated, whereby a similar notation as in [18,19]  is used. After that, the kinetic equations are presented, and the quantities necessary for their derivatives with respect to the structural design parameters are discussed.

Marker kinematics
In the floating frame of reference formulation, the deformation of a flexible body i is described with regard to a reference frame R i , which undergoes large translational and rotational motions; see Fig. 2. Thus the absolute position ρ k,i of a marker k at point P on the flexible body i can be displayed as where ρ i is the position of the reference frame, R k,i is the position of marker k with respect to R i in the undeformed configuration, and u k,i is the space-and time-dependent displacement due to the body deformation. It is worth noting that all vectors in Eq. (7) are represented in the reference frame R i . The absolute orientation of a frame that is fixed in P and described by the Cartesian basis e k,i is represented by the rotation matrix S k,i , which results from two successive rotations as Thereby the first rotation matrix S i is here parameterized by the parameters β i ∈ R 3 and describes the rotation from the inertial into the reference frame, whereas D k,i represents the rotation from the reference into the marker system. Provided that there is no initial rotation and assuming the deformations to be small and linear elastic, D k,i can be approximated as where ϑ k,i ∈ R 3 are small rotation angles, which are converted by the tilde operator(..) into a skew-symmetric matrix as follows: A. Held Both the elastic displacement u k,i and rotation ϑ k,i depend on space and time. Using a Ritz approach, space-and time-dependency can be separated as where k,i = i (R k,i ) and k,i = i (R k,i ) are the matrices of the shape functions defined over the entire body, and q i ∈ R n i q are time-dependent weighting factors, which are often denoted as elastic coordinates.
The position and orientation of an arbitrary marker k on the flexible body i can be expressed by the position variables Next to the position and orientation, also the velocity and angular velocity of marker k are required. Again, all kinematic quantities are represented in the body reference frame R i . For the absolute velocity and angular velocity of marker k represented in R i it holds, and where v i and ω i are the velocity and angular velocity of the reference frame R i . The velocities (13) and (14) can be written in a more compact form as and using the velocity coordinates of body i, an auxiliary matrix for the translation and an auxiliary matrix for the rotation Finally, the absolute acceleration and angular acceleration of marker k represented in R i are and respectively. Like the velocities, they can be displayed in compact notation as and where and

Derivatives of marker kinematics
From the Eqs. (7) and (8) it can be seen that the position and orientation of marker k depends both implicitly and explicitly on the design variables p. On the one hand, the position variables z i I depend implicitly via the system equations on time t and the design variables p. On the other hand, explicit dependencies of the position R k,i in the undeformed configuration or the shape functions k,i and k,i can occur depending on the parameterization of the body.
To provide a systematic way to compute the partial derivatives of the system equations with respect to the design variables, the derivatives of the marker kinematics are derived first. They are required later to compute the derivatives of the constraint equations, the reaction forces, and the applied discrete forces and torques, which act at marker k.
Differentiating the marker position (7) and rotation (8) with respect to the lth component of p yields and To compute the partial derivatives of the velocity, angular velocity, acceleration, and angular acceleration, only the partial derivatives of the auxiliary matrices T k,i t and T k,i r and of the auxiliary vectors ζ k,i t and ζ k,i r have to be provided. They are A. Held It is worth mentioning that the derivatives (28), (30), and (31) are not constant but depend on the position and velocity coordinates of body i.

Relative marker kinematics
Kinematic constraints lock 1 to 6 degrees of freedom between markers to model joints or bearings. Since they are of particular importance, their formulation and their derivatives are presented in the following. To this end, at first, the relative marker kinematics at joint s defined between marker k on body i and marker m on body j is briefly presented; see Fig. 3. A detailed description of the relative marker kinematics can be found in [18]. For readability, the upper index (. . . ) k,i is replaced by (. . . ) t to denote the "to" marker, and (. . . ) m,j is replaced by (. . . ) f for the "from" marker. Thus the absolute positions of the "to" and "from" markers are The relative position d s of the "from" and "to" markers represented in the "to" marker frame P t is The relative rotation between the "from" and "to" markers is described by the rotation matrix B s (β s ), β s ∈ R 3 , which can be computed by To derive the relative velocity and relative angular velocity of the "from" and "to" markers, at first, the absolute time derivative of the relative joint position d s with respect to the inertial frame is computed. Thereby V s is the relative time derivative of the relative marker position, and t ω t is the angular velocity of the "to" marker frame represented in P t . The latter can be obtained by transforming Eq. (16) into the "to" marker frame as The absolute time derivative of the relative joint position (35) can alternatively be expressed by the absolute velocities v f and v t of the "from" and "to" markers. Transforming the absolute velocities into the "to" marker frame it holds, Rearranging Eq. (37) for the relative joint velocity V s and using the compact notation (15) for the absolute marker velocities yield The relative angular velocity between the "from" and "to" markers represented in P t is given by which can be expressed in terms of the velocity coordinates of the "from" and "to" bodies as In the same way as the relative velocities, the relative acceleration and angular acceleration between the "from" and "to" markers can be determined aṡ respectively.
A. Held

Derivatives of relative marker kinematics
The relative marker kinematics is the basis to formulate the constraint equations. Hence their derivatives with respect to the design variables are required to compute the partial derivatives of the constraint equations with respect to p. In the following, the partial derivatives of the relative marker kinematics with respect to the design variables are presented. Thereby it is assumed that the "from" body is parameterized. The partial derivative of the relative marker position d s with respect to the lth component of p is where ∂ρ f /∂p l can be determined from Eq. (26). In contrast, the partial derivatives of the relative marker orientation represented by β s cannot be directly computed. The rotation parameters β s are only auxiliary variables, which depend on the position variables of the "from" and "to" bodies, as can be seen from Eq. (34). However, to find their derivatives with respect to the design variables p, the implicit equation can be rewritten in vector form as The total derivative of Eq. (45) yields ∂χ s ∂β s β s + ∂χ s ∂p p = 0.
Using either three independent equations from (46) or the pseudo-inverse ∂χ s /∂β s T ∂χ s /∂β s −1 , the partial derivatives of the relative orientation parameters with respect to the design variables can be determined from Eq. (46) as The partial derivatives of the relative marker velocity V s and relative angular marker velocity s with respect to p l can be computed by Eqs. (28), (29), and (43) as and On design sensitivities in the structural analysis and optimization. . .
Finally, the derivatives of the relative marker accelerationV s and relative angular ac-celeration˙ s with respect to the design parameters p can be determined using the derivatives of the marker kinematics (28), (29), (30), and (31) and of the relative marker kinematics (43), (48), and (49) as and respectively.

Constraint equations
Using the relative marker kinematics presented in Sect. 2.3, implicit kinematical constraint equations can be formulated. For example, if all degrees of freedom are locked at joint s it holds Differentiating Eq. (52) with respect to time t yields the implicit constraint equations at velocity level for joint s in terms of the relative joint coordinateṡ It can be seen that to compute the time derivative of the rotation parametersβ s , the kinematic relation X s r is required, which depends in turn on the rotation parameters β s ; see [17]. With Eqs. (38) and (40) the relative joint velocity and angular velocity can be expressed in terms of the velocity coordinates of the "from" and "to" bodies as Plugging Eq. (54) into Eq. (53) yields the constraints at velocity level in terms of the local Jacobian matrix C s and the velocity coordinates z i II and z j II of the "to" and "from" bodies. Representation (55) is very useful as it shows the structure of C s , whose derivatives will be required in the differentiation of the reaction forces.
Differentiating the constraint equations (53) once more with respect to time t , the constraint equations at acceleration level are obtained as with the relative joint velocities x II .

Derivatives of constraint equations
In this final step of the differentiation of the kinematics with respect to the design variables, the partial derivatives of the constraint equations with respect to p are given. With Eq. (43) and (47), the partial derivatives of the constraint equations at position level with respect to p l are simply Derivatives (57) can be used in the evaluation of gradient (6) as part of the derivatives of the implicit initial conditions at position level φ 0 .
To represent the partial derivatives of the constraint equations at velocity level with respect to p, there are two ways. On the one hand, ∂ċ s /∂p l can be written in terms of the joint coordinates as whereby the partial derivatives of the velocity and of the angular velocity are given in Eq. (48) and (49), respectively. In Eq. (58), the partial derivatives of the kinematic matrix are needed. They can be determined from the derivatives of X s with respect to the rotation parameters β s and the derivatives of the rotation parameters with respect to the design variables. Whereas the former depend on the choice of the rotation parameters, the latter are given by Eq. (47).
On the other hand, ∂ċ s /∂p l can be expressed in terms of the velocity coordinates of the "to" and "from" bodies as where Derivatives (58) or (60) can be used in the evaluation of gradient (6) as part of the derivatives of the implicit initial conditions at velocity levelφ 0 . Moreover, the partial derivatives of the local Jacobian (61) are required for the differentiation of the reaction forces in the kinetic equations.
Finally, the derivatives of the constraint equations at the acceleration level, which are required in the solution of the adjoint problem (3) In Eq. (64), all quantities are known.

Kinetic equations
For the evaluation of the gradient equation (6) With the variation of the velocity coordinates δv k,i = T k,i t δz i II and Eq. (22) for the acceleration of point P, the virtual power of the inertias yields The integral in the first term represents the local mass matrix M i . It depends on the elastic coordinates and can be written with Eq. (18) as In the second term of Eq. (67), h i ω contains the inertias from the acceleration of the reference frame. In the literature, h i ω is often separated into translational, rotational, and elastic parts as follows: Thereby the three parts are defined as and For the evaluation of the virtual power of inner forces δP i e , first, the virtual distortion velocity δε k,i and the stress σ k,i in point P have to be determined. Assuming that the deformations are small and that the material behavior is linearly elastic, it holds for the distortion velocity and its variationε and, for the stress, with the linear differential operator L L , the material matrix E i , and the distortion matrix B k,i . The virtual power of the inner forces can then be written in terms of the velocity coordinates of body i as with the vector h i e of the inner forces Incorporating the virtual velocities δv k,i and angular velocities δω k,i = T k,i r δz i II into the virtual power of the applied body forces and discrete loads yields Using the principle of the independent variation, the local equations of motion for body i follow from the Eqs. (67), (73), and (75) as It can be seen from this very brief derivation of the kinetic equations that for the evaluation of Eq. (77), a set of body integrals, which are not constant but dependent on the elastic coordinates has to be solved. A summary of the body integrals is given in Table 1.
The body integrals play an important role in the calculation of the design sensitivities because differentiating the local equations of motion (77) with respect to the design variables, their derivatives must essentially be computed. However, since the parameterization of the flexible body is assumed to be arbitrary, it can be everything from the body domain i 0 (p), the parameters of the material matrix E i (p) and the shape functions (p) and (p). In most proprietary and academic multibody programs, neither the required information on the body domain, material properties, nor shape functions are available. An integrated design modeler and finite element program would be required for this task. Therefore it is suggested to provide the partial derivatives of the body integrals via a standardized data interface. There already exists a standardized data set for flexible multibody systems. These so-called standard input data and their suggested augmentation are presented in the next section.

Augmentation of object-oriented standard input data
The original standard input data have been proposed in [21] and pursue essentially two goals. On the one hand, they represent a general and unified dataset, which completely describes flexible bodies of multibody systems. This includes all necessary information to evaluate the marker kinematics and the body integrals. On the other hand, they are meant to accelerate the evaluation of the equations of motion. To this end, the body integrals in Table 1, most of which depend on the elastic coordinates, are expressed by Taylor series. Since the coefficient matrices of the Taylor series can be evaluated in a preprocessing step, the repeated quadrature of the body integrals during the time integration is avoided.
In [21] the standard input data are stored in an object-oriented dataset using the four classes modal, refmod, node, and taylor. A detailed list of their properties is given in Tables 2, 3, 4, 5. It can be seen that the class modal contains all properties of the flexible body. This includes general body properties, which are stored in the variable refmod. The latter is of class refmod and defines, for instance, the global properties mass or the number of elastic coordinates nelastq.
The quantities required to describe all markers on the body and their kinematics are given in an array of objects of class node; see Table 4. In these objects, for example, the position and orientation of the markers are stored. Since these kinematic quantities depend on the elastic coordinates, they are of type taylor; see Table 5. In the class taylor the coefficient matrices of the Taylor series and their dimensions are saved. The remaining properties of the class modal are all of type taylor because they contain the coefficient matrices of the body integrals.
It is worth mentioning that the standard input data can be directly computed from finite element models without quadrature if these are assembled from continuum elements only; see [18]. In case there are additional structural elements in the finite element model, an approximate calculation is still possible.
Due to the object-oriented structure of the dataset, its augmentation is comparatively easy. All necessary changes in the class definitions are marked in gray in the tables. In fact, there are only two changes. Firstly, the property mass in class refmod is now of type taylor instead of double as in the original definition. Secondly, the taylor class has to be augmented by the property nd, which describes the number of design variables, and dM0, dM1, and dMn, which contain the partial derivatives of the coefficient matrices. By these simple changes all necessary information to compute the design sensitivities is provided. On the one hand, the design sensitivities of the marker kinematics in Eqs. (26)-(31) and the subsequent relative marker kinematics and constraints can be computed. On the other hand, the design sensitivities of the kinetic equations (65), which are based on the derivatives of the local equations of motion and hence on the derivatives of the body integrals in Table 1, can be determined.
Depending on the number of design variables, the amount of data can rise significantly. Therefore a careful examination of the influence of the derivatives of the individual coefficient matrices with respect to the overall gradient should be carried out as it has been done for large-scale topology optimization problems in [13]. In this way, it is possible to drop derivatives of body integrals with a low impact and limit the total size of the standard input data.

Application example: flexible slider-crank mechanism
In this section, an example from topology optimization, which closely resembles that in [1], is presented to demonstrate the practical use of the augmented standard input data in the structural analysis and optimization of multibody systems. Thereby the optimal design for a flexible crank of a slider-crank mechanism shall be found; see Fig. 4. To this end, the material properties of an underlying finite element model, from which the shape functions for the crank are obtained by modal truncation are parameterized by the so-called solid isotropic material with penalization (SIMP) approach [2]. In this method the density ρ i and stiffness E i of each finite element i = 1(1)n e are penalized by design variables p i . The finite element model of the flexible crank is parameterized with the SIMP formulation ρ i = cp q i ρ 0 for p min = 0.01 ≤ p i < 0.1, p i ρ 0 for 0.1 ≤ p i ≤ 1, A. Held   where c = 10 5 , q = 6, and s = 3. This SIMP approach works well for finite element models that are modally reduced; see [15]. The material properties and geometry of the flexible crank and the rigid body data of the piston rod are summarized in Table 6.
Determining the augmented standard input data for this problem is not trivial. The following procedure is used in this work: a parameterized finite element model is established, and modal analysis is performed in the commercial finite element solver ANSYS. Then the finite element system matrices and the results of the modal analysis are imported in Matlab. In the next step the standard input data are computed as described in [18], and their derivatives are determined. It is worth noting that for the current example, all body integrals in Table 1 depend on the design variables p. The type of dependency is irrelevant for the evaluation of the gradient equation. Therefore it is sufficient to extend the standard input data by the derivatives of the coefficient matrices to the design variables.  With the standard input data at hand, the transient analysis of the flexible slider-crank mechanism can be performed. The crank rotation is prescribed by the rheonomic constraint A. Held is evaluated as objective function. Thereby q and K are the elastic coordinates and the modally reduced stiffness matrix of the flexible crank, respectively. After solving the adjoint problem (3)-(4), the gradient ∇ψ can be computed according to Eq. (6) using the design sensitivities of the system equations discussed in Sect. 2. The results of the gradient computation for p i = 1 are shown in Fig. 5. The sensitivities of the single design variables are plotted over the discretized design domain of the flexible crank. The results of the adjoint variable method are validated using the finite difference method (forward differences, perturbations 5 · 10 −2 ). It turns out that both results are in good agreement.
Even though the application example is chosen from topology optimization, the augmented standard input data concept can be directly transferred to any structural analysis and optimization of both rigid and flexible bodies in multibody systems.

Summary and conclusion
The adjoint variable method can be used to derive a set of equations for multibody systems from which the gradient of arbitrary criterion functions can be computed. In this method, among others, the partial derivatives of the system equations with respect to the design variables are required. These so-called design sensitivities are presented in detail for both the kinematical constraint equations and the equations of motion of multibody systems modeled with the floating frame of reference formulation. Thereby it has been shown that for the systematic differentiation of the constraints equations, the marker kinematics derivatives and the relative marker kinematics must be determined. When computing the partial derivatives of the kinetic equations, it turns out that the design sensitivities of the body integrals are required. Since they cannot be determined for arbitrarily parameterized rigid and flexible bodies without an integrated design modeler and finite element program, augmenting the standard input data by their derivatives has been suggested. With only minor changes in the original definition of the standard input data, it is possible to clearly separate the two tasks of modeling and parameterization of bodies and the system simulation and analysis. To demonstrate the capability of the suggested augmentation of the standard input data, the sensitivity analysis required in the topology optimization of a flexible slider-crank mechanism has been presented.