In this paper, the direct differentiation of generalized- αtime integration is derived, 1 equations are introduced and results are shown. Although generalized- αtime integration has found 2 usage, the derivation and the resulting equations for the analytical sensitivity analysis via direct 3 differentiation are missing. Thus here, the sensitivity equations of generalized- αtime integration via 4 direct differentiation are provided. Results with generalized- αare compared with Newmark- βand 5 their sensitivities with numerical sensitivities via forward finite differencing in terms of accuracy 6 and performance. An example is shown for each linear structural dynamics and flexible multibody 7 dynamics. 8 Keywords: sensitivity analysis; time integration; structural dynamics; flexible multibody dynamics 9 1. Introduction 10 The analytical design sensitivity analysis is developed in this paper for the generalized- α 11 family of time integration. After deriving and introducing the design sensitivity using direct 12 differentiation, results are compared with numerical differentiation via forward differencing 13 in terms of accuracy and computational effort. 14 Efficient sensitivity analysis is the centerpiece of efficient and effective gradient-based 15 design optimization. Other uses of sensitivity analysis include uncertainty analysis and 16 quantification. The adjoint variable method applied to the sensitivity analysis with respect 17 to uncertain and random variables is developed in [ 1]. Overviews of design sensitivity 18 analysis are provided by [ 2–6]. In transient problems, sensitivity analysis must be extended 19 to time integration [7,8], which is handled here. 20 A number of studies have been carried out for sensitivity analysis in transient structural 21 dynamics and multibody dynamics. Early works on sensitivity analysis with transient 22 structural dynamics include [ 9], which bridges sensitivities in static and dynamic analysis via 23 the adjoint variable method, including maximum-value objectives and constraints. In [ 10], 24 sensitivity analysis with adjoint variable method and direct differentiation is introduced for 25 time-average and time-specific constraints. The German-language monograph [ 11] clearly 26 demonstrates analytic sensitivity analysis of transient linear dynamics with Newmark time 27 integration, providing all necessary equations for implementation, which also includes the 28 sensitivity of the duration of a dynamic operation. The adjoint variable method for transient 29 structural dynamics with Newmark time integration is derived by [ 12]. The adjoint variable 30 method and direct differentiation of nonlinear and transient problems is handled by [ 13] in 31 which the crucial aspects of the nonlinear solver for sensitivity analysis are introduced. 32 The sensitivity analysis of multibody systems is often limited to the governing equations 33 as for direct differentiation with the divide-and-conquer based modal superposition method 34 in [14], for direct differentiation and the adjoint variable method with general formulations 35 Version December 4, 2023 submitted to Machines https://www.mdpi.com/journal/machines Version December 4, 2023 submitted to Machines 2 of 35 of differential–algebraic equations in [ 15], for direct differentiation and the adjoint variable 36 method of ordered differential equations in [ 16] and for the adjoint variable method with the 37 floating frame of reference formulation in [ 17]. The adjoint variable method of the governing 38 equations in descriptor form and for generic optimization functions is shown by [ 18] for 39 several time integration methods including backward difference formula, Newmark, Adams– 40 Bashforth–Moulton and Diagonally Implicit Runge–Kutta. The derivation of the adjoint 41 variable method based on generalized- αtime integration and applied to the finite-element 42 based formulation is provided by [19]. 43 Higher-order sensitivities can also be of interest, especially in cases of low number of 44 degrees of freedom and design variables. Second-order senstivities via the adjoint variable 45 method with the Runge–Kutta–Fehlberg method are shown in [ 20]. Second-order sensitivities 46 for multibody systems with the adjoint variable method and differential–algebraic equations 47 are shown in [ 21]. Higher-order sensitivites are used in [ 22,23] to develop a surrogate model 48 for optimization. This work is limited to first-order sensitivities as needed for gradient- 49 based optimization, though the equations provided can be readily extended to higher-order 50 derivatives. 51 The generalized- αfamily of time integration methods offers several time integrators via 52 a single implementation and the generalized- αmethod is therefore chosen here. The α-family 53 of methods add numerical damping to the system to avoid numerical stability problems [ 24]. 54 Generalized- αtime integration method has found usages in both nonlinear structural 55 dynamics and multibody dynamics. Multibody dynamics addresses the dynamic response 56 with kinematic constraints. Rigid multibody dynamics considers the bodies rigid, while 57 flexible multibody dynamics allows for deformation in the bodies of the system. Applications 58 are shown including both linear structural dynamics and flexible multibody dynamics. 59 In the dynamic systems addressed, the transient response (response in time domain) 60 is of interest and therefore a temporal discretization is applied. The calculation of the 61 system responses, here position, velocity and acceleration, can be carried out with explicit 62 or implicit time integration methods. Implicit time integration allows for larger time steps 63 than explicit methods, which is bounded by the Courant criterion to remain stable [25]. 64 Newmark introduced an implicit integration scheme, now referred to as Newmark- β, 65 in 1959 with two integration parameters γandβ[26]. The Hilber–Hughes–Taylor (HHT) 66 method was introduced in 1977 to improve upon Newmark with the additional integration 67 parameter αffor numerical damping [ 27], while the Wood–Bossak–Zienkiewicz (WBZ, also 68 known as simply Bossak or Wood) introduced the parameter αm[28]. Generalized- αwas 69 introduced by Chung and Hulbert in 1993 [ 24] to generalize all these methods with the 70 integration parameters αm,αf,βandγ. These methods can be generalized via a single 71 implementation and is therefore referred to as the generalized- αfamily of time integrators. 72 In the generalized- αmethod, for which the family is named, all parameters can be controlled 73 by a single value, namely the spectral radius at infinity1ϱ∞ranging from zero, being the 74 highest amount of damping, to unity, being none [ 24]. The mapping of ϱ∞toαm,αf,β 75 andγis shown below. The choice of this value can affect the response for a given time step, 76 see e.g. [29]. 77 This paper is concerned with the sensitivity analysis, which is used when one or more 78 parameters effect the system responses. The calculation of the system responses is referred 79 to here as primal analysis to distinguish from the sensitivity analysis. This is of interest 80 for efficient design optimization via gradient-based optimization algorithms in addition to 81 investigating the sensitivities themselves or for uncertainty analysis. Methods for sensitivity 82 analysis include analytical, numerical, complex step and algorithmic approaches [5]. 83 Analytical sensitivities can be derived either manually (by hand) or via automatic 84 differentiation of the implemented code. Automatic differentiation and manual differentiation 85 of rigid multibody systems are compared in [ 30]. Although prone to human error and not 86 1It should be noted that we differentiate between variants of the Greek letter rho: ϱis used here for the numerical damping term, while ρis used below for density. Version December 4, 2023 submitted to Machines 3 of 35 automatic, the former represents an exact and one of the fastest methods for sensitivity 87 analysis [ 31]. The high implementation effort of the manual sensitivity analysis though 88 results in a typically highly reusable implementation and must therefore only carried out once. 89 Due to its physics- and knowledge-based nature in addition to the increased transparency, 90 hand differentiation is preferred by the authors. 91 This study limits itself to analytical sensitivity analysis, which is further divided into 92 direct differentiation (also known as forward sensitivity analysis) and the adjoint variable 93 method (also known as backward sensitivity analysis). Generally, it can be said that direct 94 differentiation is more efficient than the adjoint variable method if the number of design 95 variables is less than the number of functions to differentiate, i.e. in optimization the sum 96 of the number of objective functions (one in case of single objective optimization) and the 97 number of constraints. Due to its nature, the adjoint variable method must be performed 98 after the time integration of the primal analysis and integrated backwards in time. Thus, the 99 adjoint variable method, in contrast to direct differentiation, needs dedicated solver routines 100 and cannot reuse these implementations [ 7,8]. Further, for nonlinear analysis the system 101 matrices must be stored or recalculated when integrating backwards in time. Therefore, 102 direct differentiation is chosen in this work. The reader is referred to the following literature 103 for the adjoint variable method: [17,18,32–34]. 104 The scope here is thus the investigation of direct differentiation of the generalized- α 105 method. A method in this regard is introduced in [ 35], yet the needed equations are not 106 derived and shown. The present work is an extension of [ 7,36] where Newmark- βsensitivities 107 are developed, which are in turn based on [ 11,12]. In this paper, Newmark- βserves as a 108 comparison. 109 This work first generally introduces generalized- αtime integration and its equations. 110 Afterwards, the sensitivity equations are introduced for structural dynamics and flexible 111 multibody dynamics, respectively. A numerical example is shown for each and results 112 are provided of both primal and sensitivity analysis in comparison with Newmark- βtime 113 integration. We close with a conclusion, including discussion of this work. Further, the 114 detailed derivation of the effective system matrices is provided in the appendix. 115 2. Sensitivity analysis with generalized- αtime integration 116 In this section, we introduce the equations for an acceleration-form ( a-form) of 117 generalized- αtime integration using a predictor–corrector scheme. The acceleration-form 118 refers to the formulation, where the acceleration is the variable for which is being solved, 119 and is contrasted by the displacement-form ( d-form), where the displacement is the primary 120 variable. The predictor–corrector scheme allows for nonlinear terms in the governing equa- 121 tion in which the predictor values are used to determine these values in the nonlinear solver, 122 e.g. Newton–Raphson. 123 In this work, we extend this well-established time integration method to include design 124 sensitivity analysis, which is defined by the differentiation with respect to the design variables 125 (but could include other parameters, e.g. uncertain parameters). Sensitivity analysis is 126 carried out here by directly differentiating both the equations of the time integration as 127 well as the governing equations. We consider both linear structural dynamics as well 128 as flexible multibody dynamics, which is intrinsically nonlinear. After introducing the 129 general equations for generalized- αtime integration and its sensitivities, the domain specific 130 governing equations for linear structural dynamics and flexible multibody dynamics are given. 131 Flexible multibody dynamics is shown using the specific approach developed in [ 7,37,38] 132 with an index-1 approach and the floating frame of reference formulation. 133 Version December 4, 2023 submitted to Machines 4 of 35 tn−1 tn tn+1−αf tn+1−αmtn+1 tn+2¨qn−1¨qn¨qn+1¨qn+2 αfαm¨qn ¨qn+1 time t Figure 1. Time integration with —acceleration, —γ=0.5,—γ=0.75 2.1. Generalized- αtime integration 134 The method for primal and sensitivity analysis with generalized- αtime integration is introduced in the following. The basis of which are Newmark’s equations for the updated position qand velocity vectors ˙qwith the calculated value of the acceleration vector ¨q, qn+1=qn+∆t˙qn+1 2−β ∆t2¨qn | {z } predictor qpred+β∆t2¨qn+1|{z} corrector, (1) ˙qn+1=˙qn+(1−γ)∆t¨qn| {z } predictor ˙qpred+γ∆t¨qn+1|{z} corrector, (2) where γandβare the Newmark integration parameters. Though the simulation is evaluated at the times tnandtn+1, values at intermediate points are interpolated. The position q and velocity ˙quse time tn+1−αfand the acceleration ¨quses time tn+1−αm, giving the approximations by linear interpolation (see fig. 1), qn+1−αf=(1−αf)qn+1+αfqn, (3) ˙qn+1−αf=(1−αf)˙qn+1+αf˙qn, (4) ¨qn+1−αm=(1−αm)¨qn+1+αm¨qn, (5) where αmandαfare integration parameters. A predictor–corrector scheme is used in order to accommodate nonlinear calculation and state-dependency of the system parameters, i.e. system matrices with respect to position and velocity. The predicted state variables with Newmark’s equations are qpred=qn+∆t˙qn+1 2−β ∆t2¨qn, (6) ˙qpred=˙qn+(1−γ)∆t¨qn, (7) and we calculate ¨qn+1using the effective system equations (see below for explanation of structural dynamics and flexible multibody dynamics, respectively). The predicted values are then corrected with Newmark’s equations, qn+1=qpred+β∆t2¨qn+1, (8) ˙qn+1=˙qpred+γ∆t¨qn+1. (9) Version December 4, 2023 submitted to Machines 5 of 35 For the sensitivity analysis within the generalized- αtime integration, we use direct differen- tiation of the approximations, ∇qn+1−αf=(1−αf)∇qn+1+αf∇qn, (10) ∇˙qn+1−αf=(1−αf)∇˙qn+1+αf∇˙qn, (11) ∇¨qn+1−αm=(1−αm)∇¨qn+1+αm∇¨qn, (12) where the nabla operator ∇denotes the design sensitivity, i.e. total derivative with respect 135 to the design variables, e.g. for the sensitivity scalar response with respect to the design 136 variable vector, 137 ∇(·)=d(·) dx. (13) The same direct differentiation is applied to the predictors, which are then defined by ∇qpred=∇qn+∆t∇˙qn+1 2−β ∆t2∇¨qn, (14) ∇˙qpred=∇˙qn+(1−γ)∆t∇¨qn. (15) Aftercalculating ∇¨qn+1byuseoftheeffectivesensitivitysystem(seebelowforexplanationof structural dynamics and flexible multibody dynamics, respectively), the following correctors are used: ∇qn+1=∇qpred+β∆t2∇¨qn+1, (16) ∇˙qn+1=∇˙qpred+γ∆t∇¨qn+1. (17) For nonlinear analysis, a Newton solver can be integrated into the solving routine. This is 138 explained with flexible multibody systems below in Sec. 2.3. 139 The time integration constants of generalized- αareαm,αf,βandγand through their choice can lead to a number of established time integration methods, including the Newmark family of integrators, Hilber–Hughes–Taylor and Wood–Bossak–Zienkiewicz. Therefore, this can then be referred to as the α-family of time integrators. The equations introduced above in this section can therefore be used generally within this family. Specifically for the generalized- αtime integration used here, the choice of integration constants can be tied to a single user-defined parameter, the spectral radius at infinity ϱ∞[24], αm=2ϱ∞−1 ϱ∞+1, (18) αf=ϱ∞ ϱ∞+1, (19) γ=1 2−αm+αf, (20) β=1 4(1−αm+αf)2. (21) The spectral radius at infinity is sometimes referred to as the (numerical) damping and 140 is limited to a range of [0, 1]. The value ϱ∞=1 returns Newmark- βretaining frequencies 141 are preserved, while ϱ∞=0 annihilates all frequencies above∆t Tafter one time step in 142 which Tis the period of the highest frequency of interest. The results converge to those 143 without numerical damping when the time steps goes to zero. Thus, the level of numerical 144 damping introduced to generalized- αtime integration depends on both the step size as well 145 as the spectral radius at infinity and the choice is critical for proper results. A parameter 146 study with fixed time step and varying the values of the spectral radius at infinity is shown 147 in [29]. As the present study develops a general framework for this methodology, this is not 148 addressed here and the reader is referred to [ 24,39–41]. In this paper, we use a moderate 149 Version December 4, 2023 submitted to Machines 6 of 35 value of ϱ∞=0.55, though the sensitivity analysis methodology developed is general and 150 independent of this parameter value. 151 2.2. Structural dynamics 152 The governing equation for linear structural dynamics is defined by the following second-order equation of motion: m¨q(t)+d˙q(t)+kq(t)−F(t)=0, (22) where mis the mass matrix, dis the damping matrix, kis the stiffness matrix and Fis 153 the external force vector. For linear structural dynamics the system matrices are constant 154 and independent of position and time. For the solving scheme, we add the residual force 155 vector R(t), which is to be zero (to an allowable value) at all times, to the right-hand side. 156 The α-family of integration methods use intermediate point approximations between tnand 157 tn+1, namely tn+1−αfandtn+1−αm, revealing 158 Rn+1−αf=m¨qn+1−αm+d˙qn+1−αf+kqn+1−αf−Fn+1−αf, (23) where the linear approximation of the intermediate external force is 159 Fn+1−αf=(1−αf)Fn+1+αfFn. (24) The effective equation of motion is defined by 160 meff¨qn+1=Feff, (25) where meff=(1−αm)m+(1−αf)γ∆td +(1−αf)β∆t2k, (26) Feff=(1−αf)Fn+1+αfFn+ −αmm¨qn−d (1−αf)˙qpred+αf˙qn −k (1−αf)qpred+αfqn . (27) For the derivation of the effective mass matrix meffand effective force vector Fefffor 161 generalized- αmethod, see Sec. A.1. 162 The sensitivity analysis is derived using the direct differentiation which is applied to the governing equation (23). This reveals the sensitivity of the force balance for the α-family of integration methods as ∇Rn+1−αf=∇m¨qn+1−αm+m∇¨qn+1−αm+∇d˙qn+1−αf+d∇˙qn+1−αf+ +∇kqn+1−αf+k∇qn+1−αf−∇Fn+1−αf, (28) where the intermediate external force sensitivity is found with the approximation 163 ∇Fn+1−αf=(1−αf)∇Fn+1+αf∇Fn. (29) The design sensitivities of the system matrices are three-dimensional matrices and therefore 164 efficient handling is needed. 165 The effective equation of motion for sensitivities is defined by the differentiation of (25), meff∇¨qn+1=Fpseudo, (30) where 166 Fpseudo =∇Feff−∇meff¨qn+1, (31) Version December 4, 2023 submitted to Machines 7 of 35 and fully expressed by Fpseudo =(1−αf)∇Fn+1+αf∇Fn−∇m (1−αm)¨qn+1+αm¨qn −αmm∇¨qn+ −∇d (1−αf)˙qn+1+αf˙qn −d (1−αf)∇˙qpred+αf∇˙qn + −∇k (1−αf)qn+1+αfqn −k (1−αf)∇qpred+αf∇qn . (32) For the derivation of the effective mass matrix meffand effective pseudo force Fpseudofor 167 generalized- αmethod, see Sec. A.1 and Sec. A.2. 168 2.3. Flexible multibody dynamics 169 Flexible multibody dynamics is the analysis of mechanical systems consisting of flexible 170 bodies constrained by kinematic joints modeled with differential–algebraic equations. This 171 work builds upon previous work of analytical sensitivity analysis with multibody dynamics, 172 including [ 15,42,43]. In addition to the kinematic constraints w.r.t. structural dynamics, 173 multibody dynamics also must integrate a nonlinear solver in order to consider dependencies 174 of the system matrices and vectors on the state variables. This is carried out with Newton 175 steps. 176 The flexible multibody system is modeled with the floating frame of reference formula- tion [44–46]. The generic index-3 differential–algebraic equation for a flexible multibody system with holonomic constraints is written as follows: m q,t¨q(t)+d q,t˙q(t)+k q,t q(t)+JT Φ q,t λ(t)=Fext q,˙q,t +Fv q,˙q,t , (33) Φ q,t =0, (34) where q,˙qand¨qare the generalized position, velocity and acceleration vectors, λis the 177 vector of Lagrange multipliers of kinematic constraints, the matrices m,d, and kare the 178 mass, damping and stiffness matrices, JΦis the Jacobian matrix of kinematic constraints, 179 Fextis the vector of external forces, Fvis the quadratic velocity vector and Φis the vector 180 of kinematic constraints at position level. 181 Numerical problems related to the solution of index-3 differential–algebraic equa- tions [47] can be avoided by a number of methods including index reduction, coordinate reduction, Augmented Lagrangian formulation and preconditioning [ 48,49]. Although generalized- αtime integration is suitable for solving index-3 or index-2 differential–algebraic equations directly, in this work we address the widely used index reduction leading to an index-1 differential–algebraic equation in which the kinematic constraint function is differentiated twice with respect to time, m¨q+d˙q+kq+JT Φλ=Fext+Fv, (35) JΦ¨q=Fc, (36) where Fcis the right-hand side of acceleration constraints. As index-1 differential–algebraic 182 equations are susceptible to drift in the kinematic constraints, Baumgarte stabilization [ 50] 183 is used. Kinematic constraint drift is stabilized in this method by adding damping and 184 restoring terms to the acceleration constraint equations that are proportional to the velocity 185 and position constraints. Details of its implementation, including the sensitivity analysis 186 are described in [51,52]. 187 The generalized- αtime integration is based on 188  m¨q n+1−αm+ d˙q n+1−αf+ kq n+1−αf+ JT Φλ n+1−αm=Fext,n+1−αf+Fv,n+1−αf, (37)  JΦ¨q n+1−αm=Fc,n+1−αf, (38) Version December 4, 2023 submitted to Machines 8 of 35 with the intermediate approximations of all force terms of the equations of motion. These intermediate approximations are shown in Sec. B.1. Inserting the intermediate approximations of the force terms acting on the flexible multibody system leads to the equations of motion in residual form given by R1,n+1−αf=(1−αm)mn+1¨qn+1+αmmn¨qn+(1−αf)dn+1˙qn+1+αfdn˙qn+ +(1−αf)kn+1qn+1+αfknqn+(1−αm)JT Φ,n+1λn+1+αmJT Φ,nλn+ −(1−αf)Fext,n+1−αfFext,n−(1−αf)Fv,n+1−αfFv,n, (39) R2,n+1−αf=(1−αm)JΦ,n+1¨qn+1+αmJΦ,n¨qn−(1−αf)Fc,n+1−αfFc,n. (40) In the case of a nonlinear mass matrix, as is the case with the floating frame of reference 189 formulation, this weighted formulation of the residual equations are first-order accurate. 190 For a formulation of second-order accuracy see e.g. [ 53], where the dynamic equilibrium 191 is enforced exactly at the time steps of the simulation. The derived time integration and 192 its differentiation with the weighted residual formulation can be extended and applied to 193 further formulations. 194 The effective equation of motion is defined by 195 meffJT Φ,eff JΦ,eff0¨qn+1 λn+1 =Fa,eff Fc,eff , (41) where meff=(1−αm)mn+1+(1−αf)γ∆tdn+1+(1−αf)β∆t2kn+1, (42) JΦ,eff=(1−αm)JΦ,n+1, (43) Fa,eff=(1−αf)Fext,n+1+αfFext,n+(1−αf)Fv,n+1+αfFv,n+ −αmmn¨qn−(1−αf)dn+1˙qpred−αfdn˙qn+ −(1−αf)kn+1qpred−αfknqn−αmJT Φ,nλn, (44) Fc,eff=(1−αf)Fc,n+1+αfFc,n−αmJΦ,n¨qn. (45) The derivation of the effective equation of motion and the terms meff,JΦ,eff,Fa,effandFc,eff196 is shown in Sec. B.1. 197 With flexible multibody dynamics, the system parameters (matrix- and vector-valued 198 forndof>1) including mass, damping, stiffness, constraints and forces are generally 199 dependent on state variables. Therefore, a nonlinear solver is needed to solve the system 200 and the Newton–Raphson method is used here of the following form: 201  ∂R1,n+1−αf ∂¨qn+1∂R1,n+1−αf ∂λn+1 ∂R2,n+1−αf ∂¨qn+1∂R2,n+1−αf ∂λn+1 ∆¨qn+1 ∆λn+1 +R1,n+1−αf R2,n+1−αf =0. (46) Version December 4, 2023 submitted to Machines 9 of 35 With the dependencies shown in (33) and (34), the terms of the Jacobian matrix (tangent operator) of the residuals with respect to ¨qn+1andλn+1are given by, ∂R1,n+1−αf ∂¨qn+1=(1−αm) β∆t2∂mn+1 ∂qn+1¨qn+1+mn+1! + +(1−αf) β∆t2∂dn+1 ∂qn+1˙qn+1+γ∆tdn+1! + +(1−αf) β∆t2∂kn+1 ∂qn+1qn+1+β∆t2kn+1! + +(1−αm)β∆t2∂JT Φ,n+1 ∂qn+1λn+1+ −(1−αf) β∆t2∂Fext,n+1 ∂qn+1+γ∆t∂Fext,n+1 ∂˙qn+1! + −(1−αf) β∆t2∂Fv,n+1 ∂qn+1+γ∆t∂Fv,n+1 ∂˙qn+1! , (47) ∂R2,n+1−αf ∂¨qn+1=(1−αm) β∆t2∂JΦ,n+1 ∂qn+1¨qn+1+JΦ,n+1! + −(1−αf) β∆t2∂Fc,n+1 ∂qn+1+γ∆t∂Fc,n+1 ∂˙qn+1! , (48) ∂R1,n+1−αf ∂λn+1=(1−αm)JT Φ,n+1, (49) ∂R2,n+1−αf ∂λn+1=0. (50) The derivation of these terms is shown in Sec. B.2. 202 The direct differentiation of the equations of motion (39) and (40) leads to the governing equation of the sensitivity analysis, ∇R1,n+1−αf=(1−αm)∇m n+1¨qn+1+αm∇m n¨qn+ +(1−αm)mn+1∇¨q n+1+αmmn∇¨q n+ +(1−αf)∇d n+1˙qn+1+αf∇d n˙qn+ +(1−αf)dn+1∇˙q n+1+αfdn∇˙q n+ +(1−αf)∇k n+1qn+1+αf∇k nqn+ +(1−αf)kn+1∇q n+1+αfkn∇q n+ +(1−αm)∇JT Φ,n+1λn+1+αm∇JT Φ,nλn+ +(1−αm)JT Φ,n+1∇λn+1+αmJT Φ,n∇λn+ −(1−αf)∇Fext,n+1−αf∇Fext,n−(1−αf)∇Fv,n+1−αf∇Fv,n, (51) ∇R2,n+1−αf=(1−αm)∇JΦ,n+1¨qn+1+αm∇JΦ,n¨qn +(1−αm)JΦ,n+1∇¨q n+1+αmJΦ,n∇¨q n+ −(1−αf)∇Fc,n+1−αf∇Fc,n. (52) Version December 4, 2023 submitted to Machines 10 of 35 The governing sensitivity equation of the effective system is defined by the differentiation of 203 (41), 204meffJT Φ,eff JΦ,eff0" ∇¨q n+1 ∇λn+1# =" Fa,pseudo Fc,pseudo# , (53) where 205 " Fa,pseudo Fc,pseudo# =" ∇Fa,eff ∇Fc,eff# −" ∇meff∇JT Φ,eff ∇JΦ,eff0#¨qn+1 λn+1 , (54) leading to Fa,pseudo=(1−αf)∇Fext,n+1+αf∇Fext,n+(1−αf)∇Fv,n+1+αf∇Fv,n+ −(1−αm)∇m n+1¨qn+1−αm∇m n¨qn−αmmn∇¨q n+ −(1−αf)∇d n+1˙qn+1−αf∇d n˙qn−(1−αf)dn+1∇˙q pred−αfdn∇˙q n+ −(1−αf)∇k n+1qn+1−αf∇k nqn−(1−αf)kn+1∇q pred−αfkn∇q n+ −(1−αm)∇JT Φ,n+1λn+1−αm∇JT Φ,nλn−αmJT Φ,n∇λn, (55) Fc,pseudo=(1−αf)∇Fc,n+1+αf∇Fc,n+ −(1−αm)∇JΦ,n+1¨qn+1−αm∇JΦ,n¨qn−αmJΦ,n∇¨q n. (56) The derivation of the effective mass matrix and the effective constraint Jacobian matrix is 206 shown in Sec. B.1 and the derivation of the effective pseudo force terms is shown in Sec. B.3. 207 The Newton–Raphson step for the sensitivity analysis is given by 208  ∂∇R1,n+1−αf ∂∇¨qn+1∂∇R1,n+1−αf ∂∇λn+1 ∂∇R2,n+1−αf ∂∇¨qn+1∂∇R2,n+1−αf ∂∇λn+1 " ∆∇¨q n+1 ∆∇λn+1# +" ∇R1,n+1−αf ∇R2,n+1−αf# =0. (57) As shown in [ 7] the Jacobian of the residual sensitivity with respect to ∇¨q n+1and 209 ∇λn+1reduces to 210  ∂∇R1,n+1−αf ∂∇¨qn+1∂∇R1,n+1−αf ∂∇λn+1 ∂∇R2,n+1−αf ∂∇¨qn+1∂∇R2,n+1−αf ∂∇λn+1 = ∂R1,n+1−αf ∂¨qn+1∂R1,n+1−αf ∂λn+1 ∂R2,n+1−αf ∂¨qn+1∂R2,n+1−αf ∂λn+1 e, (58) wherefore a conventional solver can be used for the solution of the following: 211  ∂R1,n+1−αf ∂¨qn+1∂R1,n+1−αf ∂λn+1 ∂R2,n+1−αf ∂¨qn+1∂R2,n+1−αf ∂λn+1 " ∆∇¨q n+1 ∆∇λn+1# +" ∇R1,n+1−αf ∇R2,n+1−αf# =0. (59) 3. Numerical results 212 In this section, numerical results are shown comparing the introduced method for the 213 sensitivity analysis for generalized- αtime integration (primal analysis introduced in Sec. 2) 214 with that of Newmark- β. The methods are compared with applications to both structural 215 dynamics as well as flexible multibody dynamics. 216 Version December 4, 2023 submitted to Machines 11 of 35 k1 k2 q1 q2 q3F1 m1 m2 m3 Figure 2. Structural dynamics example 3.1. Structural dynamics 217 The linear structural dynamic problem used here as a demonstration example is based 218 on [54], which has found usage as a benchmark in [ 55–61], and shown schematically in fig. 2. 219 This problem is defined by 220 mz }| { m10 0 0m20 0 0 m3 ¨qz}|{ ¨q1 ¨q2 ¨q3 +kz }| { k1 −k10 −k1k1+k2−k2 0 −k2 k2 qz}|{ q1 q2 q3 −Fz}|{ F1 0 0 =Rz}|{ 0 0 0 , (60) with the prescribed displacement, 221 q1=sinωt, (61) and the unitless values m1=0, m2=1, m3=1, k1=107, k2=1, ω=1.2. The masses m2andm3as well as the stiffnesses k1andk2are considered design variables 222 with respect to which the design sensitivities are calculated. The problem consists of stiff 223 (k1)and flexible ( k2)parts typical of complex, practical structural-mechanical systems [ 54]. 224 Considering the prescribed displacement, this problem can be reduced to two degrees 225 of freedom, 226 mz}|{m20 0m3¨qz}|{¨q2 ¨q3 +kz}| {k1+k2−k2 −k2 k2qz}|{q2 q3 −Fz}|{k1q1 0 =Rz}|{0 0 , (62) which we implement in our in-house, Python-based solver framework SiMuLi in the time t=[0, 10]with the initial conditions q2=0, q3=0, ˙q2=0, ˙q3=0. Version December 4, 2023 submitted to Machines 12 of 35 0246810−101 q2 0246810−202 ˙q2 0246810−5000500 ¨q2 0246810−2024 time tq3 0246810−4−2024 time t˙q3 0246810−4−202 time t¨q3 Figure 3. Responses of the structural dynamic example calculated with —Newmark- βand —generalized- αcompared with —exact solution (undermost line) As reported in [ 54], time steps of ∆t=0.2618 are used, giving the following values for the ratio step size to fundamental periods of the system: ∆t TF=0.05, ∆t T1=0.041667, ∆t T2=131.76, where TFis the fundamental period of the forcing function and T1,2for the first and second 227 modes, T1=6.283 and T2=0.001987. 228 The results for both Newmark- βand generalized- αtime integration methods are shown 229 in fig. 3 and compared with the exact solution. The exact solution is calculated via mode 230 superposition per [ 62]. The numerical damping parameter of generalized- αtime integration 231 is set as ϱ∞=0.55. As expected with generalized- αtime integration, the unphysical, 232 high-frequent oscillatory behavior of the response can be avoided. As such, the beneficial 233 effects of numerical damping with generalized- αtime integration can be clearly seen in the 234 velocity and acceleration of node 2, while maintaining the integrity of the other responses. 235 The sensitivity analysis developed above is applied and these results can be seen in 236 fig. 4–6. Analytical sensitivity analysis with Newmark- βand generalized- αtime integration 237 methods are compared and also with their numerical sensitivity implementations via finite 238 differencing using a relative perturbation of 1 ×10−6. 239 In cases where Newmark- βtime integration shows instability or noise, the numerical 240 damping allows for usable sensitivities. In some cases, noise — albeit at a much lower 241 level — persists in the first iterations before being damped out: after approximately 2.5s 242 (or approximately 10 iterations). Attention is drawn to those sensitivities susceptible to 243 stability issues, namely: 244 •sensitivities of position 2 w.r.t. stiffness 1 and mass 2 (fig. 4), 245 •sensitivities of velocity 2 w.r.t to stiffness 1 and mass 2 (fig. 5), 246 •sensitivities of acceleration 2 w.r.t. stiffness 1, stiffness 2 and mass 2 (fig. 6), 247 •sensitivities of acceleration 3 w.r.t. stiffness 1 and mass 2 (fig. 6). 248 Further, generalized- αtime integration eliminates the noise in those sensitivities present 249 when using Newmark- β. This is seen with the following sensitivity responses: 250 •sensitivities of velocity 2 w.r.t. stiffness 2 (fig. 5), 251 •sensitivities of velocity 3 w.r.t. stiffness 1 and mass 2 (fig. 5). 252 Version December 4, 2023 submitted to Machines 13 of 35 0246810−505·10−12 ∇q2w.r.t. k1 0246810−101·10−6w.r.t. k2 0246810−505·10−5w.r.t. m2 0246810−0.500.51·10−6w.r.t. m3 0246810−1−0.500.5·10−13 time t∇q3 0246810−10−505 time t0246810−2024·10−7 time t0246810−50510 time t Figure 4. Position sensitivities ∇qof the structural dynamic example w.r.t. k1,k2,m2andm3 calculated analyically with —Newmark- βand—generalized- αcompared with those calculated numerically via finite differencing with —Newmark- βand—generalized- α(undermost lines). Note: no visible difference between finite differencing and analytical sensitivities. The analytical and numerical results are in the most plots nearly exact. There are cases 253 though where the negative consequences of the numerical instability with Newmark- βare 254 increased. The readers attention is called to especially 255 •sensitivity of velocity 2 w.r.t. stiffness 2 (fig. 5), 256 •sensitivity of velocity 2 w.r.t. mass 3 (fig. 5), 257 •sensitivity of acceleration 2 w.r.t. stiffness 2 (fig. 6), 258 •sensitivity of acceleration 2 w.r.t. mass 3 (fig. 6). 259 It can be seen that both analytical and numerical sensitivities with generalized- αtime 260 integration are not negatively effected by numerical instabilities as the case with Newmark- 261 β. 262 The number of time steps needed for the damping to take effect can be reduced by 263 reducing the numerical damping parameter ϱ∞towards zero (i.e. more damping), though 264 also effecting the physical response to a greater degree. As time continues, the instability 265 grows with Newmark- βtime integration and therefore a clear advantage can be seen with 266 generalized- α. It can be clearly seen how gradient-based solvers would struggle to use the 267 gradient information provided with Newmark- βtime integration. 268 For our modest case of four design variables, we show a speedup of two, where the 269 speedup Sis defined by the ratio of computational time for finite differencing tFDto that of 270 analytical sensitivities tana, 271 S=tFD tana. (63) The computational time required for analytical sensitivity analysis represents approximately 272 the ideal computational effort twice that of the primal analysis, i.e. analysis without 273 sensitivity analysis. Computation reduction is the same for both Newmark- βand generalized- 274 αtime integration methods, see tab. 1. 275 3.2. Flexible slider–crank mechanism 276 The slider–crank mechanism as shown in fig. 7 is modeled with flexible multibody 277 dynamics using the floating frame of reference formulation. The mechanism consists of four 278 bodies that are connected by three revolute joints and one prismatic joint. The ground 279 and the slider are considered to be rigid, while the crank and the rod are considered to be 280 flexible and are modeled with Euler–Bernoulli beams. Fig. 8 shows the flexible multibody 281 Version December 4, 2023 submitted to Machines 14 of 35 0246810−202·10−9 ∇˙q2w.r.t. k1 0246810−101·10−6w.r.t. k2 0246810−202·10−2w.r.t. m2 0246810−0.500.51·10−6w.r.t. m3 0246810−1−0.500.5·10−13 time t∇˙q3 0246810−10−505 time t0246810−505·10−7 time t0246810−50510 time t Figure 5. Velocity sensitivities ∇˙qof the structural dynamic example w.r.t. k1,k2,m2andm3 calculated analyically with —Newmark- βand—generalized- αcompared with those calculated numerically via finite differencing with —Newmark- βand—generalized- α(undermost lines). Note: no visible difference between finite differencing and analytical sensitivities. 0246810−4−2024·10−7 ∇¨q2w.r.t. k1 0246810−505·10−5w.r.t. k2 0246810−4−2024w.r.t. m2 0246810−101·10−6w.r.t. m3 0246810−505·10−12 time t∇¨q3 0246810−10010 time t0246810−505·10−5 time t0246810−10010 time t Figure 6. Acceleration sensitivities ∇¨qof the structural dynamic example w.r.t. k1,k2,m2and m3calculated analyically with —Newmark- βand—generalized- αcompared with those calculated numerically via finite differencing with —Newmark- βand—generalized- α(undermost lines). Note: no visible difference between finite differencing and analytical sensitivities. Table 1. Computational effort for the structural dynamics example time integration sensitivity analysis calculation time [s]1speedup2 Newmark- β none 0.0555 – Newmark- β finite differencing 0.2088 – Newmark- β analytical 0.1023 2.0410 generalized- α none 0.0458 – generalized- αfinite differencing 0.2122 – generalized- α analytical 0.0993 2.1370 1Intel Core i5-7200U CPU 2.50GHz, 2speedup of analytical sensitivity analysis w.r.t. finite differencing Version December 4, 2023 submitted to Machines 15 of 35 OA B groundslidercrank rod xy xCyCxRyR xSyS (a)Schematic diagramcrank wchcrod wrhr (b)Cross sections Figure 7. Slider–crank mechanism xyxC yCxRyR xSyS xSθCxC,f yC,fθC,f Figure 8. Flexible multibody system of the slider–crank mechanism with the undeformed bodies in gray, the deformed bodies in black, the inertial reference frame in blue, the body reference frames in red and the coordinates that are investigated in green system of the mechanism. The inertial reference frame coinciding with the reference frame 282 of the ground is shown in blue. The body reference frames of the crank, the rod and the 283 slider are shown in red. The undeformed bodies of the crank and the rod are shown in gray 284 and the deformed bodies are shown in black, where the points represent the nodes of the 285 finite elements. The coordinates that are investigated in detail for the primal analysis and 286 the sensitivity analysis are shown in green. The dimensions of the bodies are listed in tab. 2. 287 The bodies are made of steel and tab. 3 shows the material properties. In this example, the 288 MPa unit system is used with length in mm, mass in t, time in s, force in N, angle in rad 289 and stress in MPa. 290 The generalized coordinates of a flexible multibody system modeled with the floating frame of reference formulation and Euler–Bernoulli beams are given by q=rTθqT fT, (64) where ris the position of the body reference frame, θis the orientation of the body reference frame and qfare the nodal displacements in body coordinates of each node on the flexible Table 2. Geometry of slider–crank mechanism property symbol crank rod slider units width w 20 20 30 mm height h 30 30 30 mm length ℓ120 180 40 mm Version December 4, 2023 submitted to Machines 16 of 35 Table 3. Material properties property symbol value units density ρ7.85×10−9t mm3 elastic modulus E 210000 MPa Poisson ratio ν 0.3 − body. Since no damping is considered, the equations of motion of a flexible body are given by m¨q+kq+JT Φλ=Fext+Fv, (65) JΦ¨q=Fc, (66) where the system matrices of the flexible bodies are given by m= mttmtrmtf mrrmrf sym. mff , (67) k= 0 0 0 0 0 sym. kff , (68) and the system vectors are given by Fext= Fext,t Fext,r Fext,f , (69) Fv= Fv,t Fv,r Fv,f , (70) with the index denoting translational terms t, rotational terms rand flexible terms f. 291 For the dynamic simulation of the mechanism, the above-mentioned solver framework 292 SiMuLi is used. The simulation is performed in the time interval t=[0s, 0.025s ]with time 293 steps ∆t=1·10−4s. The resulting values for the ratio step size to fundamental periods 294 of the system for the first ten non-zero eigenfrequencies are found in fig. 9. The values for 295 ratio step size to fundamental periods range from 0.2119 for first non-zero eigenfrequency 296 to 3.9919 for the tenth non-zero eigenfrequency and to 17.1157 for the twentieth non-zero 297 eigenfrequency (not pictured in plot). 298 For the initial conditions, the angular position and velocity of the crank are defined with θC=0rad, ˙θC=150πrad s, and the initial conditions of the other bodies are derived through the kinematic analysis of 299 the system. The dimensions of the cross sections of the crank hC,wCand the rod hR,wR300 are used as design variables in the sensitivity analysis. 301 The results of the simulation of the coordinates that are highlighted in fig. 8 are shown 302 in fig. 10. With Newmark- β, vibrations at high frequencies are observed. With generalized- α, 303 vibrations are observed in the beginning of the simulation that are numerically damped with 304 ongoing simulation time. In addition to the generalized coordinates, the maximum stresses 305 on the upper and lower fiber of the beam elements of the flexible bodies including the 306 Version December 4, 2023 submitted to Machines 17 of 35 0501001502002503003501234 crank orientation θC[deg]ratio ∆t/T Figure 9. Ratio time step to fundamental periods for first ten non-zero eigenfrequencies of the slider–crank mechanism crank and the rod are computed. The Kreisselmeier–Steinhauser function [ 63,64] is used to 307 approximate the maximum stress value of the bodies at each time step with a differentiable 308 function and the results are shown in fig. 11. With Newmark- β, the vibrating behavior of 309 the system causes higher flexible deformation and therefore unrealistic high stress in the 310 components. The numerical damping with generalized- αleads to lower maximum stress 311 values and a more realistic stress distribution over time. 312 The results of the sensitivity analysis of the coordinates that are highlighted in fig. 8 313 are shown in figs. 12–14. The analyzed coordinates include the rotation of the floating frame 314 of the crank θC, the translation in x-direction of the floating frame of the slider xSand the 315 three flexible coordinates of the node in A on the crank xC,f,yC,fandθC,f. The design 316 variables used for the sensitivity analysis are the cross-sectional dimensions of the crank hC, 317 wCand the rod hR,wR. 318 Similar to above, with Newmark- βthe results show a vibrating behavior that would lead 319 to significant problems when using the sensitivities in gradient-based design optimization. 320 With generalized- α, the numerical damping leads to usable sensitivity values after a few 321 time steps. Fig. 12 shows that the position sensitivities of the position and orientation 322 coordinates of the floating frames are less noisy compared to the flexible coordinates. Position 323 sensitivities that suffer little from noise include the following: 324 •sensitivities of θCw.r.t. wCandwR, 325 •sensitivities of xSw.r.t. hC,wC, and wR. 326 Position sensitivities with significant noise include the following: 327 •sensitivities of θCw.r.t. hCandhR, 328 •sensitivities of xSw.r.t. hR, 329 •sensitivities of xC,fw.r.t. . hC,wC,hRandwR, 330 •sensitivities of yC,fw.r.t. . hC,wC,hRandwR, 331 •sensitivities of θC,fw.r.t. . hC,wC,hRandwR. 332 The sensitivities of velocities as shown in fig. 13 and the sensitivities of accelerations as 333 shown in fig. 14 have noisy responses with Newmark that are damped with generalized- αfor 334 all investigated coordinates and all used design variables. The same behavior also continue 335 with the sensitivity values for the maximum stresses in the crank and the rod shown in 336 fig. 15. 337 To check the correctness of the results, the sensitivity values with the introduced 338 method are compared to numerical sensitivity analysis with forward differencing. The 339 results are in good agreement, with only small numerical differences. In fact there are no 340 visible differences in figs. 12–15 with the finite differencing values lying under the respective 341 analytical values. Tab. 4 shows a comparison of the computational effort with different 342 methods. The computational effort with generalized- αis about 30% lower compared to 343 Version December 4, 2023 submitted to Machines 18 of 35 00.010.0202468 θC 00.010.02200300400500 ˙θC 00.010.02−202·106 ¨θC 00.010.02100200300 xS 00.010.02−4−2024·104 ˙xS 00.010.02−505·107 ¨xS 00.010.0205·10−2 xC,f 00.010.02−1,000−50005001,000 ˙xC,f 00.010.02−202·107 ¨xC,f 00.010.02−0.200.2 yC,f 00.010.02−5,00005,000 ˙yC,f 00.010.02−202·108 ¨yC,f 00.010.02−4−2024·10−3 time tθC,f 00.010.02−1000100 time t˙θC,f 00.010.02−4−2024·106 time t¨θC,f Figure 10. Reponses of the slider–crank mechanism (values in MPa units, i.e. mm,mm s,mm s2for translations and rad,rad s,rad s2for rotations) calculated with —Newmark- βand—generalized- α 00.010.020100200 time tσU C 00.010.020100200 time tσL C 00.010.020100200300 time tσU R 00.010.020100200 time tσL R Figure 11. Maximum stress in MPa on the upper and lower fiber of the crank and the rod calculated with —Newmark- βand—generalized- α Version December 4, 2023 submitted to Machines 19 of 35 00.010.020123·10−2 ∇θCw.r.t. hC 00.010.02024·10−2w.r.t. wC 00.010.02−8−6−4−20·10−3w.r.t. hR 00.010.02−6−4−20·10−3w.r.t. wR 00.010.02−202 ∇xS 00.010.02−4−202 00.010.0200.5 00.010.02−0.500.5 00.010.02−0.100.1 ∇xC,f 00.010.02−202·10−2 00.010.02−505·10−2 00.010.02−505·10−3 00.010.02−1−0.500.5 ∇yC,f 00.010.02−0.2−0.100.10.2 00.010.02−0.200.2 00.010.02−202·10−2 00.010.02−101·10−2 time t∇θC,f 00.010.02−202·10−3 time t00.010.02−505·10−3 time t00.010.02−505·10−4 time t Figure 12. Position sensitivities ∇qof the slider–crank mechanism w.r.t. hC,wC,hRandwR (values in MPa system of units, i.e. mm,mm s,mm s2for translations and rad,rad s,rad s2for rotations) calculated with —Newmark- βand—generalized- αcompared with those calculated numerically via finite differencing with —Newmark- βand—generalized- α(undermost lines). Note: no visible difference between finite differencing and analytical sensitivities. Version December 4, 2023 submitted to Machines 20 of 35 00.010.02−2000200 ∇˙θCw.r.t. hC 00.010.02−50050w.r.t. wC 00.010.02−100−50050w.r.t. hR 00.010.02−10−50510w.r.t. wR 00.010.02−5,00005,000 ∇˙xS 00.010.02−4,000−2,0000 00.010.02−4,000−2,00002,0004,000 00.010.020500 00.010.02−2,00002,000 ∇˙xC,f 00.010.02−5000500 00.010.02−1,00001,000 00.010.02−1000100 00.010.02−202·104 ∇˙yC,f 00.010.02−5,00005,000 00.010.02−5,00005,000 00.010.02−1,00001,000 00.010.02−400−2000200400 time t∇˙θC,f 00.010.02−100−50050100 time t00.010.02−1000100 time t00.010.02−20−1001020 time t Figure 13. Velocity sensitivities ∇˙qof the slider–crank mechanism w.r.t. hC,wC,hRandwR (values in MPa system of units, i.e. mm,mm s,mm s2for translations and rad,rad s,rad s2for rotations) calculated with —Newmark- βand—generalized- αcompared with those calculated numerically via finite differencing with —Newmark- βand—generalized- α(undermost lines). Note: no visible difference between finite differencing and analytical sensitivities. Version December 4, 2023 submitted to Machines 21 of 35 00.010.02−101·107 ∇¨θCw.r.t. hC 00.010.02−202·106w.r.t. wC 00.010.02−505·106w.r.t. hR 00.010.02−2024·105w.r.t. wR 00.010.02−101·108 ∇¨xS 00.010.02−202·107 00.010.02−505·107 00.010.02−505·106 00.010.02−0.500.51·108 ∇¨xC,f 00.010.02−1012·107 00.010.02−2024·107 00.010.02−505·106 00.010.02−101·109 ∇¨yC,f 00.010.02−202·108 00.010.02−505·108 00.010.02−505·107 00.010.02−2−1012·107 time t∇¨θC,f 00.010.02−4−2024·106 time t00.010.02−1−0.500.51·107 time t00.010.02−505·105 time t Figure 14. Acceleration sensitivities ∇¨qof the slider–crank mechanism w.r.t. hC,wC,hRandwR (values in MPa system of units, i.e. mm,mm s,mm s2for translations and rad,rad s,rad s2for rotations) calculated with —Newmark- βand—generalized- αcompared with those calculated numerically via finite differencing with —Newmark- βand—generalized- α(undermost lines). Note: no visible difference between finite differencing and analytical sensitivities. Version December 4, 2023 submitted to Machines 22 of 35 00.010.02−5000500 ∇σU Cw.r.t. hC 00.010.02−1000100w.r.t. wC 00.010.02−1000100200w.r.t. hR 00.010.02−20020w.r.t. wR 00.010.02−5000500 ∇σL C 00.010.02−1000100 00.010.02−1000100 00.010.02−20−1001020 00.010.02−2000200400 ∇σU R 00.010.02−50050 00.010.02−5000500 00.010.02−20020 00.010.02−400−2000200400 time t∇σL R 00.010.02−50050 time t00.010.02−5000500 time t00.010.02−20020 time t Figure 15. Maximum stress sensitivities on the upper and lower fiber of the crank and the rod (values in MPa system of units, i.e. mm,mm s,mm s2for translations and rad,rad s,rad s2for rotations) calculated with —Newmark- βand—generalized- αcompared with those calculated numerically via finite differencing with —Newmark- βand—generalized- α(undermost lines). Note: no visible difference between finite differencing and analytical sensitivities. Version December 4, 2023 submitted to Machines 23 of 35 Table 4. Computational effort for the slider–crank mechanism time integration sensitivity analysis calculation time [mm:ss]1speedup2 Newmark- β none 03:44 – Newmark- β finite differencing 18:12 – Newmark- β analytical 04:29 4.0595 generalized- α none 02:42 – generalized- αfinite differencing 13:22 – generalized- α analytical 03:30 3.8190 1Intel Core i7-8700 CPU 3.20GHz 2speedup of analytical sensitivity analysis w.r.t. finite differencing the computational effort with Newmark- β. The main reason for that is the continuing 344 vibrating behavior with Newmark- β, requiring more iteration steps of the nonlinear solver. 345 The comparison of the sensitivity methods with numerical finite differencing and analytical 346 direct differentiation shows the high efficiency of the analytical method. With the four 347 shown design variables, a speedup as shown in (63) of about S≈4 of the analytical method 348 compared to the numerical method is obtained. The additional computational effort of 349 analytical sensitivity analysis is just a fraction of the computational effort of the primal 350 analysis. This is because of the simplification shown in (59), which allows to use the Jacobian 351 matrix of the residual built in the primal analysis also for the sensitivity analysis. This 352 allows for a further reduced computational effort. 353 4. Conclusion 354 Analytical sensitivity analysis is crucial for efficient design optimization, parameter 355 fitting and uncertainty analysis. In this work, the direct differentiation of the generalized- α 356 time integration is derived and implemented for linear structural dynamics and flexible 357 multibody dynamics, the latter of which requiring a nonlinear solve routine. All equations 358 needed for implementation are provided and extendable to further governing equations. The 359 beneficial effects of the step-size dependent numerical damping provided by the generalized- α 360 method can be clearly seen in the comparison of the results with Newmark- β. 361 The objective of this work is to develop the analytical sensitivity analysis for the 362 α-family of time integration, which has a single implementation for several common time 363 integration schemes. This thus has consequences for gradient-based optimization frameworks 364 for dynamic systems. 365 If the functions of interest are nonsmooth, discontinuous and bifurcated, non-gradient- 366 based methods can be successfully used, e.g. design cases with crash analysis. Such cases 367 can be solved with gradient-free optimization algorithms (e.g. [ 65–68]) and approximation 368 methods (e.g. [69–71]). 369 The system responses with respect to the design variables are indeed often smooth, 370 continuous and unbifurcated and thus enabling gradient-based methods where efficient and 371 accurate sensitivity analysis is critical for both effectiveness and efficiency. This is shown 372 above with linear structural dynamics and flexible multibody dynamics. The method and 373 equations developed are of general nature and are extendable to further formulations and 374 applicable to larger problems, both of larger number of degrees of freedom and number of 375 design variables. 376 With the above study and introduction of the direct differentiation of the α-family of 377 time integrators, the authors hope to provide the readers as design engineers a tool in the 378 pursuit of optimal structures and mechanisms under dynamic loading. 379 Author Contributions: Conceptualization, E.W. and V.G.; methodology, E.W. and V.G.; 380 software, E.W. and V.G.; validation, E.W. and V.G.; formal analysis, E.W. and V.G.; investigation, 381 E.W. and V.G.; resources, E.W. and V.G.; data curation, E.W. and V.G.; writing—original draft 382 Version December 4, 2023 submitted to Machines 24 of 35 preparation, E.W. and V.G.; writing—review and editing, E.W. and V.G.; visualization, E.W. and 383 V.G.; supervision, E.W.; project administration, E.W.; funding acquisition, E.W. All authors have 384 read and agreed to the published version of the manuscript. 385 Funding: This work is supported by the project CRC 2017 TN2091 doloMULTI Design of 386 Lightweight Optimized structures and systems under MULTI disciplinary considerations through 387 integration of MULTI body dynamics in a MULTI physics framework, funded by the Free University 388 of Bozen-Bolzano. 389 Data Availability Statement: Data available upon request. 390 Acknowledgments: This work is supported by the project CRC 2017 TN2091 doloMULTI Design 391 ofLightweight Optimized structures and systems under MULTI disciplinary considerations through 392 integration of MULTI body dynamics in a MULTI physics framework, funded by the Free University 393 of Bozen-Bolzano. 394 Conflicts of Interest: The authors declare no conflict of interest 395 Appendix A. Derivation of effective equations of motion for structural dynamics 396 Appendix A.1. Primal analysis 397 The base for the α-family of integration methods is the equation of motion which for for linear structural dynamics is 398 given by 399 m¨qn+1−αm+d˙qn+1−αf+kqn+1−αf=Fn+1−αf, (A1) and the intermediate values shown in eqs. (3)–(5). Inserting the intermediate values into the equation of motion leads to m (1−αm)¨qn+1+αm¨qn +d (1−αf)˙qn+1+αf˙qn +k (1−αf)qn+1+αfqn = (1−αf)Fn+1+αfFn. (A2) With Newmark’s eqs, (1) and (2), this becomes m (1−αm)¨qn+1+αm¨qn +d (1−αf) ˙qpred+γ∆t¨qn+1 +αf˙qn + +k (1−αf) qpred+β∆t2¨qn+1 +αfqn =(1−αf)Fn+1+αfFn. (A3) This can be rearranged to the following form  (1−αm)m+(1−αf)γ∆td +(1−αf)β∆t2k¨qn+1= (1−αf)Fn+1+αfFn−αmm¨qn−d (1−αf)˙qpred+αf˙qn −k (1−αf)qpred+αfqn . (A4) Therefore, the effective equation of motion is defined by 400 meff¨qn+1=Feff, (A5) where meff=(1−αm)m+(1−αf)γ∆td +(1−αf)β∆t2k, (A6) Feff=(1−αf)Fn+1+αfFn−αmm¨qn−d (1−αf)˙qpred+αf˙qn −k (1−αf)qpred+αfqn . (A7) Appendix A.2. Sensitivity analysis 401 The governing equation for sensitivities is the first derivative of the equation of motion, ∇m¨qn+1−αm+m∇¨qn+1−αm+∇d˙qn+1−αf+d∇˙qn+1−αf+∇kqn+1−αf+k∇qn+1−αf= ∇Fn+1−αf, (A8) Version December 4, 2023 submitted to Machines 25 of 35 where we introduce the intermediate sensitivity values (10)–(12), ∇m (1−αm)¨qn+1+αm¨qn +m (1−αm)∇¨qn+1+αm∇¨qn + +∇d (1−αf)˙qn+1+αf˙qn +d (1−αf)∇˙qn+1+αf∇˙qn + +∇k (1−αf)qn+1+αfqn +k (1−αf)∇qn+1+αf∇qn = (1−αf)∇Fn+1+αf∇Fn, (A9) and which we update with Newmark’s sensitivity eqs. (16) and (17) ∇m (1−αm)¨qn+1+αm¨qn +m (1−αm)∇¨qn+1+αm∇¨qn + +∇d (1−αf)˙qn+1+αf˙qn +d (1−αf) ∇˙qpred+γ∆t∇¨qn+1 +αf∇˙qn + +∇k (1−αf)qn+1+αfqn +k (1−αf) ∇qpred+β∆t2∇¨qn+1 +αf∇qn = (1−αf)∇Fn+1+αf∇Fn. (A10) This equation can be rearranged to the following form,  (1−αm)m+(1−αf)γ∆td +(1−αf)β∆t2k ∇¨qn+1= (1−αf)∇Fn+1+αf∇Fn−∇m (1−αm)¨qn+1+αm¨qn −αmm∇¨qn+ −∇d (1−αf)˙qn+1+αf˙qn −d (1−αf)∇˙qpred+αf∇˙qn + −∇k (1−αf)qn+1+αfqn −k (1−αf)∇qpred+αf∇qn . (A11) Therefore, the effective equation of motion for sensitivities is defined by 402 meff∇¨qn+1=Fpseudo, (A12) where meffis given by eq. (A6) and the pseudo load is defined by Fpseudo =(1−αf)∇Fn+1+αf∇Fn−∇m (1−αm)¨qn+1+αm¨qn −αmm∇¨qn+ −∇d (1−αf)˙qn+1+αf˙qn −d (1−αf)∇˙qpred+αf∇˙qn + −∇k (1−αf)qn+1+αfqn −k (1−αf)∇qpred+αf∇qn . (A13) Appendix B. Derivation of effective equations of motion for flexible multibody dynamics 403 Appendix B.1. Primal analysis 404 The base of the generalized- αmethod for index-1 differential algebraic equations of a flexible multibody system is given by  m¨q n+1−αm+ d˙q n+1−αf+ kq n+1−αf+ JT Φλ n+1−αm=Fext,n+1−αf+Fv,n+1−αf, (A14)  JΦ¨q n+1−αm=Fc,n+1−αf. (A15) Version December 4, 2023 submitted to Machines 26 of 35 The intermediate point approximations of the force terms acting on the flexible multibody system are given by  m¨q n+1−αm=(1−αm)mn+1¨qn+1+αmmn¨qn, (A16)  d˙q n+1−αf=(1−αf)dn+1˙qn+1+αfdn˙qn, (A17)  kq n+1−αf=(1−αf)kn+1qn+1+αfknqn, (A18)  JT Φλ n+1−αm=(1−αm)JT Φ,n+1λn+1+αmJT Φ,nλn, (A19)  JΦ¨q n+1−αm=(1−αm)JΦ,n+1¨qn+1+αmJΦ,n¨qn, (A20) Fext,n+1−αf=(1−αf)Fext,n+1+αfFext,n, (A21) Fv,n+1−αf=(1−αf)Fv,n+1+αfFv,n, (A22) Fc,n+1−αf=(1−αf)Fc,n+1+αfFc,n. (A23) Applying these linear approximations of the generalized- αmethod to the equations of motion leads to (1−αm)mn+1¨qn+1+αmmn¨qn+ +(1−αf)dn+1˙qn+1+αfdn˙qn+ +(1−αf)kn+1qn+1+αfknqn+ +(1−αm)JT Φ,n+1λn+1+αmJT Φ,nλn=(1−αf)Fext,n+1+αfFext,n+ +(1−αf)Fv,n+1+αfFv,n, (A24) (1−αm)JΦ,n+1¨qn+1+αmJΦ,n¨qn=(1−αf)Fc,n+1+αfFc,n. (A25) Inserting Newmark’s eqs. (1) and (2) gives (1−αm)mn+1¨qn+1+αmmn¨qn+ +(1−αf)dn+1 ˙qpred+γ∆t¨qn+1 +αfdn˙qn+ (1−αf)kn+1 qpred+β∆t2¨qn+1 +αfknqn+ (1−αm)JT Φ,n+1λn+1+αmJT Φ,nλn=(1−αf)Fext,n+1+αfFext,n+ +(1−αf)Fv,n+1+αfFv,n, (A26) (1−αm)JΦ,n+1¨qn+1+αmJΦ,n¨qn=(1−αf)Fc,n+1+αfFc,n. (A27) which can be rewritten to  (1−αm)m+(1−αf)γ∆td +(1−αf)β∆t2k¨qn+1+ +(1−αm)JT Φ,n+1λn+1=(1−αf)Fext,n+1+αfFext,n+ +(1−αf)Fv,n+1+αfFv,n−αmmn¨qn+ −(1−αf)dn+1˙qpred−αfdn˙qn+ −(1−αf)kn+1qpred−αfknqn+ −αmJT Φ,nλn, (A28) (1−αm)JΦ,n+1¨qn+1=(1−αf)Fc,n+1+αfFc,n−αmJΦ,n¨qn. (A29) The equations of motion of the effective system is given by 405 meffJT Φ,eff JΦ,eff0¨qn+1 λn+1 =Fa,eff Fc,eff , (A30) Version December 4, 2023 submitted to Machines 27 of 35 where meff=(1−αm)mn+1+(1−αf)γ∆tdn+1+(1−αf)β∆t2kn+1, (A31) JΦ,eff=(1−αm)JΦ,n+1, (A32) Fa,eff=(1−αf)Fext,n+1+αfFext,n+(1−αf)Fv,n+1+αfFv,n−αmmn¨qn+ −(1−αf)dn+1˙qpred−αfdn˙qn−(1−αf)kn+1qpred−αfknqn−αmJT Φ,nλn, (A33) Fc,eff=(1−αf)Fc,n+1+αfFc,n−αmJΦ,n¨qn. (A34) Appendix B.2. Nonlinear solver 406 The base for the nonlinear solver are the residual equations of the flexible multibody system, given by R1,n+1−αf=(1−αm)mn+1¨qn+1+αmmn¨qn+(1−αf)dn+1˙qn+1+αfdn˙qn+ +(1−αf)kn+1qn+1+αfknqn+(1−αm)JT Φ,n+1λn+1+αmJT Φ,nλn+ −(1−αf)Fext,n+1−αfFext,n−(1−αf)Fv,n+1−αfFv,n, (A35) R2,n+1−αf=(1−αm)JΦ,n+1¨qn+1+αmJΦ,n¨qn−(1−αf)Fc,n+1−αfFc,n. (A36) The Newton–Raphson method is used here 407  ∂R1,n+1−αf ∂¨qn+1∂R1,n+1−αf ∂λn+1 ∂R2,n+1−αf ∂¨qn+1∂R2,n+1−αf ∂λn+1 ∆¨qn+1 ∆λn+1 +R1,n+1−αf R2,n+1−αf =0, (A37) where the terms of the Jacobian matrix are given by ∂R1,n+1−αf ∂¨qn+1=(1−αm)∂mn+1 ∂¨qn+1¨qn+1+αm∂mn ∂¨qn+1¨qn+(1−αm)mn+1∂¨qn+1 ∂¨qn+1+αmmn∂¨qn ∂¨qn+1+ +(1−αf)∂dn+1 ∂¨qn+1˙qn+1+αf∂dn ∂¨qn+1˙qn+(1−αf)dn+1∂˙qn+1 ∂¨qn+1+αfdn∂˙qn ∂¨qn+1+ +(1−αf)∂kn+1 ∂¨qn+1qn+1+αf∂kn ∂¨qn+1qn+(1−αf)kn+1∂qn+1 ∂¨qn+1+αfkn∂qn ∂¨qn+1+ +(1−αm)∂JT Φ,n+1 ∂¨qn+1λn+1+αm∂JT Φ,n ∂¨qn+1λn+(1−αm)JT Φ,n+1∂λn+1 ∂¨qn+1+αmJT Φ,n∂λn ∂¨qn+1+ −(1−αf)∂Fext,n+1 ∂¨qn+1−αf∂Fext,n ∂¨qn+1−(1−αf)∂Fv,n+1 ∂¨qn+1−αf∂Fv,n ∂¨qn+1, (A38) ∂R2,n+1−αf ∂¨qn+1=(1−αm)∂JΦ,n+1 ∂¨qn+1¨qn+1+αm∂JΦ,n ∂¨qn+1¨qn+(1−αm)JΦ,n+1∂¨qn+1 ∂¨qn+1+αmJΦ,n∂¨qn ∂¨qn+1+ −(1−αf)∂Fc,n+1 ∂¨qn+1−αf∂Fc,n ∂¨qn+1, (A39) Version December 4, 2023 submitted to Machines 28 of 35 ∂R1,n+1−αf ∂λn+1=(1−αm)∂mn+1 ∂λn+1¨qn+1+αm∂mn ∂λn+1¨qn+(1−αm)mn+1∂¨qn+1 ∂λn+1+αmmn∂¨qn ∂λn+1+ +(1−αf)∂dn+1 ∂λn+1˙qn+1+αf∂dn ∂λn+1˙qn+(1−αf)dn+1∂˙qn+1 ∂λn+1+αfdn∂˙qn ∂λn+1+ +(1−αf)∂kn+1 ∂λn+1qn+1+αf∂kn ∂λn+1qn+(1−αf)kn+1∂qn+1 ∂λn+1+αfkn∂qn ∂λn+1+ +(1−αm)∂JT Φ,n+1 ∂λn+1λn+1+αm∂JT Φ,n ∂λn+1λn+(1−αm)JT Φ,n+1∂λn+1 ∂λn+1+αmJT Φ,n∂λn ∂λn+1+ −(1−αf)∂Fext,n+1 ∂λn+1−αf∂Fext,n ∂λn+1−(1−αf)∂Fv,n+1 ∂λn+1−αf∂Fv,n ∂λn+1, (A40) ∂R2,n+1−αf ∂λn+1=(1−αm)∂JΦ,n+1 ∂λn+1¨qn+1+αm∂JΦ,n ∂λn+1¨qn+(1−αm)JΦ,n+1∂¨qn+1 ∂λn+1+αmJΦ,n∂¨qn ∂λn+1+ −(1−αf)∂Fc,n+1 ∂λn+1−αf∂Fc,n ∂λn+1. (A41) The derivatives of the system parameters of the previous time step with respect to state variables of the next time step disappear, ∂#n ∂ n+1=0,∂#n ∂ n+1=0. (A42) The dependencies of the system parameters are shown in (33) and (34) and applying the chain rule of differentiation leads to ∂m ∂¨q=∂m ∂q∂q ∂¨q, (A43) ∂d ∂¨q=∂d ∂q∂q ∂¨q, (A44) ∂k ∂¨q=∂k ∂q∂q ∂¨q, (A45) ∂JΦ ∂¨q=∂JΦ ∂q∂q ∂¨q, (A46) ∂Fext ∂¨q=∂Fext ∂q∂q ∂¨q+∂Fext ∂˙q∂˙q ∂¨q, (A47) ∂Fv ∂¨q=∂Fv ∂q∂q ∂¨q+∂Fv ∂˙q∂˙q ∂¨q, (A48) ∂# ∂λ=0. (A49) ∂# ∂λ=0. (A50) Version December 4, 2023 submitted to Machines 29 of 35 The derivatives of the state variables with respect to ¨qn+1andλn+1are given by ∂qn+1 ∂¨qn+1=β∆t2e,∂qn+1 ∂λn+1=0, ∂˙qn+1 ∂¨qn+1=γ∆te,∂˙qn+1 ∂λn+1=0, ∂¨qn+1 ∂¨qn+1=e,∂¨qn+1 ∂λn+1=0, ∂λn+1 ∂¨qn+1=0,∂λn+1 ∂λn+1=e. (A51) Applying eq. (A42), eqs. (A43)–(A50) and eq. (A51) leads to the terms of the Jacobian matrix given by ∂R1,n+1−αf ∂¨qn+1=(1−αm) β∆t2∂mn+1 ∂qn+1¨qn+1+mn+1! + +(1−αf) β∆t2∂dn+1 ∂qn+1˙qn+1+γ∆tdn+1! + +(1−αf) β∆t2∂kn+1 ∂qn+1qn+1+β∆t2kn+1! + +(1−αm)β∆t2∂JT Φ,n+1 ∂qn+1λn+1+ −(1−αf) β∆t2∂Fext,n+1 ∂qn+1+γ∆t∂Fext,n+1 ∂˙qn+1! + −(1−αf) β∆t2∂Fv,n+1 ∂qn+1+γ∆t∂Fv,n+1 ∂˙qn+1! , (A52) ∂R2,n+1−αf ∂¨qn+1=(1−αm) β∆t2∂JΦ,n+1 ∂qn+1¨qn+1+JΦ,n+1! + −(1−αf) β∆t2∂Fc,n+1 ∂qn+1+γ∆t∂Fc,n+1 ∂˙qn+1! , (A53) 408∂R1,n+1−αf ∂λn+1=(1−αm)JT Φ,n+1, (A54) 409∂R2,n+1−αf ∂λn+1=0. (A55) Version December 4, 2023 submitted to Machines 30 of 35 Appendix B.3. Sensitivity analysis 410 The derivative of the equations of motion leads to the governing equation of the sensitivities given by (1−αm)∇m n+1¨qn+1+αm∇m n¨qn+ +(1−αm)mn+1∇¨q n+1+αmmn∇¨q n+ +(1−αf)∇d n+1˙qn+1+αf∇d n˙qn+ +(1−αf)dn+1∇˙q n+1+αfdn∇˙q n+ +(1−αf)∇k n+1qn+1+αf∇k nqn+ +(1−αf)kn+1∇q n+1+αfkn∇q n+ +(1−αm)∇JT Φ,n+1λn+1+αm∇JT Φ,nλn+ +(1−αm)JT Φ,n+1∇λn+1+αmJT Φ,n∇λn=(1−αf)∇Fext,n+1+αf∇Fext,n+ +(1−αf)∇Fv,n+1+αf∇Fv,n, (A56) and (1−αm)∇JΦ,n+1¨qn+1+αm∇JΦ,n¨qn+ +(1−αm)JΦ,n+1∇¨q n+1+αmJΦ,n∇¨q n=(1−αf)∇Fc,n+1+αf∇Fc,n. (A57) Including Newmark’s sensitivity eqs. (16) and (17) leads to (1−αm)∇m n+1¨qn+1+αm∇m n¨qn+ +(1−αm)mn+1∇¨q n+1+αmmn∇¨q n+ +(1−αf)∇d n+1˙qn+1+αf∇d n˙qn+ +(1−αf)dn+1 ∇˙q pred+γ∆t∇¨q n+1 +αfdn∇˙q n+ +(1−αf)∇k n+1qn+1+αf∇k nqn+ +(1−αf)kn+1 ∇q pred+β∆t2∇¨q n+1 +αfkn∇q n+ +(1−αm)∇JT Φ,n+1λn+1+αm∇JT Φ,nλn+ +(1−αm)JT Φ,n+1∇λn+1+αmJT Φ,n∇λn=(1−αf)∇Fext,n+1+αf∇Fext,n+ +(1−αf)∇Fv,n+1+αf∇Fv,n, (A58) and (1−αm)∇JΦ,n+1¨qn+1+αm∇JΦ,n¨qn+ +(1−αm)JΦ,n+1∇¨q n+1+αmJΦ,n∇¨q n=(1−αf)∇Fc,n+1+αf∇Fc,n. (A59) Version December 4, 2023 submitted to Machines 31 of 35 which can be rewritten to  (1−αm)m+(1−αf)γ∆td +(1−αf)β∆t2k ∇¨qn+1+ +(1−αm)JT Φ,n+1∇λn+1=(1−αf)∇Fext,n+1+αf∇Fext,n+ +(1−αf)∇Fv,n+1+αf∇Fv,n+ −(1−αm)∇m n+1¨qn+1+ −αm∇m n¨qn−αmmn∇¨q n+ −(1−αf)∇d n+1˙qn+1−αf∇d n˙qn+ −(1−αf)dn+1∇˙q pred−αfdn∇˙q n+ −(1−αf)∇k n+1qn+1−αf∇k nqn+ −(1−αf)kn+1∇q pred−αfkn∇q n+ −(1−αm)∇JT Φ,n+1λn+1+ −αm∇JT Φ,nλn−αmJT Φ,n∇λn, (A60) and (1−αm)JΦ,n+1∇¨q n+1=(1−αf)∇Fc,n+1+αf∇Fc,n+ −(1−αm)∇JΦ,n+1¨qn+1−αm∇JΦ,n¨qn−αmJΦ,n∇¨q n. (A61) The governing equation of sensitivity analysis of the effective system is given by 411 meffJT Φ,eff JΦ,eff0" ∇¨q n+1 ∇λn+1# =" Fa,pseudo Fc,pseudo# , (A62) where Fa,pseudo=(1−αf)∇Fext,n+1+αf∇Fext,n+(1−αf)∇Fv,n+1+αf∇Fv,n+ −(1−αm)∇m n+1¨qn+1−αm∇m n¨qn−αmmn∇¨q n+ −(1−αf)∇d n+1˙qn+1−αf∇d n˙qn−(1−αf)dn+1∇˙q pred−αfdn∇˙q n+ −(1−αf)∇k n+1qn+1−αf∇k nqn−(1−αf)kn+1∇q pred−αfkn∇q n+ −(1−αm)∇JT Φ,n+1λn+1−αm∇JT Φ,nλn−αmJT Φ,n∇λn, (A63) Fc,pseudo=(1−αf)∇Fc,n+1+αf∇Fc,n+ −(1−αm)∇JΦ,n+1¨qn+1−αm∇JΦ,n¨qn−αmJΦ,n∇¨q n. (A64)