Analytical Sensitivity Analysis of Dynamic Problems with Direct Differentiation of Generalized-α Time Integration

Abstract

In this paper, the direct differentiation of generalized-$\mathsf{\alpha }$ time integration is derived, equations are introduced and results are shown. Although generalized-$\mathsf{\alpha }$ time integration has found usage, the derivation and the resulting equations for the analytical sensitivity analysis via direct differentiation are missing. Thus, here, the sensitivity equations of generalized-$\mathsf{\alpha }$ time integration via direct differentiation are provided. Results with generalized-$\mathsf{\alpha }$ are compared with Newmark-$\mathsf{\beta }$ time integration and their sensitivities with numerical sensitivities via forward finite differencing in terms of accuracy and performance. An example is shown for each linear structural dynamics and flexible multibody dynamics.

1. Introduction

Analytical design sensitivity analysis is developed in this paper for the generalized-α family of time integration. After deriving and introducing design sensitivity using direct differentiation, results are compared with numerical differentiation via forward differencing in terms of accuracy and computational effort.

Efficient sensitivity analysis is the centerpiece of efficient and effective gradient-based design optimization. Other uses of sensitivity analysis include uncertainty analysis and quantification, e.g., the adjoint variable method applied to the sensitivity analysis with respect to uncertain and random variables was developed in [1]. Overviews of design sensitivity analysis are provided by [2,3,4,5,6]. In transient problems, sensitivity analysis must be extended to time integration [7,8], which is handled here by introducing and deriving the equations for the direct differentiation of the generalized-α family of time integration.

In the dynamic systems addressed, the transient response (response in time domain) is of interest and therefore a temporal discretization is applied. The calculation of the system responses, here position, velocity and acceleration, can be carried out with explicit or implicit time integration methods. Implicit time integration allows for larger time steps than explicit methods, which is bounded by the Courant criterion to remain stable [9].

Newmark introduced an implicit integration scheme, now referred to as Newmark-β, in 1959 with two integration parameters $\gamma$ and $\beta$ [10]. The Hilber–Hughes–Taylor (HHT) method was introduced in 1977 to improve upon Newmark with the additional integration parameter ${\alpha }_f$ for numerical damping [11], while the Wood–Bossak–Zienkiewicz (WBZ, also known as simply Bossak or Wood) introduced the parameter ${\alpha }_m$ [12]. Generalized-α was introduced by Chung and Hulbert in 1993 [13] to generalize all these methods with the integration parameters ${\alpha }_m$, ${\alpha }_f$, $\beta$ and $\gamma$. These methods can be generalized via a single implementation which is therefore referred to as the generalized-α family of time integrators. In the generalized-α method, for which the family is named, all parameters can be controlled by a single value, namely the spectral radius at infinity ${\varrho }_{\infty }$ ranging from zero, being the highest amount of damping, to unity, being none [13]. It should be noted that we differentiate between the variants of the Greek letter rho: $\varrho$ is used here for the spectral radius at infinity, while $\rho$ is used below for density. The mapping of ${\varrho }_{\infty }$ to ${\alpha }_m$, ${\alpha }_f$, $\beta$ and $\gamma$ is shown below. The choice of this value can affect the response for a given time step, see, e.g., [14].

Generalized-α time integration has found usage in both nonlinear structural dynamics and multibody dynamics. Multibody dynamics addresses the dynamic response with kinematic constraints. Rigid multibody dynamics considers the bodies rigid, while flexible multibody dynamics allows for the deformation in the bodies of the system. Applications are shown including both linear structural dynamics and flexible multibody dynamics.

This paper is concerned with the sensitivity analysis of dynamic problems, which is used when one or more parameters affect the system responses. The calculation of the system responses is referred to here as primal analysis so that it may be distinguished from the sensitivity analysis. This is of interest for the efficient design optimization via gradient-based optimization algorithms in addition to investigating the sensitivities themselves or for uncertainty analyses. Methods for sensitivity analysis include analytical, numerical, complex step and algorithmic approaches [5].

Analytical sensitivities can be derived either manually (by hand) or via the automatic differentiation [15,16] of the implemented code. The automatic differentiation and manual differentiation of rigid multibody systems are compared in [17]. Although prone to human error and not automatic, the former represents an exact and one of the fastest methods for sensitivity analysis [18]. Manual sensitivity analysis requires high implementation effort, but the resulting code is typically highly reusable and must therefore only be derived and implemented once. Due to its physics- and knowledge-based nature, in addition to its increased transparency, hand differentiation is preferred by the authors.

This study limits itself to analytical sensitivity analysis, which is further divided into direct differentiation (also known as forward sensitivity analysis) and the adjoint variable method (also known as backward sensitivity analysis). Generally, it can be said that direct differentiation is more efficient than the adjoint variable method if the number of design variables is less than the number of functions to differentiate, i.e., in optimization, the sum of the number of objective functions (one in case of single objective optimization) and the number of constraints. Due to its nature, the adjoint variable method must be performed after the time integration of the primal analysis and integrated backwards in time. Thus, the adjoint variable method, in contrast to direct differentiation, needs dedicated solver routines to integrate backwards in time and cannot reuse primal forward time integration [7,8]. Furthermore, for nonlinear analysis, the system matrices must be stored or recalculated when integrating backwards in time. Therefore, direct differentiation is chosen in this work. The reader is referred to the following literature for the adjoint variable method: [19,20,21,22,23,24,25,26,27,28,29,30].

A number of studies have been carried out for sensitivity analysis in transient structural dynamics and multibody dynamics. Early works on sensitivity analysis with transient structural dynamics include [31], which bridges sensitivities in static and dynamic analysis via the adjoint variable method, including maximum-value objectives and constraints. In [32], sensitivity analysis with the adjoint variable method and direct differentiation is introduced for time-averaged and time-specific constraints. The German-language monograph [33] clearly demonstrates the analytic sensitivity analysis of transient linear dynamics with Newmark-β time integration, providing all necessary equations for implementation, which also includes the sensitivity of the duration of a dynamic operation. The adjoint variable method for transient structural dynamics with Newmark-β time integration is derived by [34]. The adjoint variable method and direct differentiation of nonlinear and transient problems is handled by [35] in which the crucial aspects of the nonlinear solver for sensitivity analysis are introduced.

The sensitivity analysis of multibody systems is often limited to the governing equations and often not extended to the nonlinear solver or time integration. Exceptions to this can be found for direct differentiation with the divide-and-conquer-based modal superposition method in [36]; for direct differentiation and the adjoint variable method with general formulations of differential–algebraic equations in [37]; for direct differentiation and the adjoint variable method of ordered differential equations in [38]; and for the adjoint variable method with the floating frame of reference formulation in [20]. The adjoint variable method of the governing equations in descriptor form and for generic optimization functions is shown by [21] for several time integration methods, including the backward difference formula, Newmark-, Adams–Bashforth–Moulton and diagonally implicit Runge–Kutta. The derivation of the adjoint variable method based on generalized-α time integration and applied to the finite-element-based formulation is provided by [39]. The transient response sensitivity analysis of localized nonlinear behavior using Newmark-β [40].

Higher-order sensitivities can also be of interest, especially in cases of a low number of degrees of freedom and design variables. Second-order sensitivities via the adjoint variable method with the Runge–Kutta–Fehlberg method are shown in [41]. Second-order sensitivities for multibody systems with the adjoint variable method and differential–algebraic equations are shown in [42]. Higher-order sensitivities are used in [43,44] to develop a surrogate model for optimization. This work is limited to first-order sensitivities, as required by gradient-based optimization, although the equations provided can be readily extended to higher-order derivatives.

The scope here is thus the investigation of the direct differentiation of the generalized-α method. A method in this regard is introduced in [45], yet the needed equations are not derived and shown. The present work builds upon and extends the works to generalized-α. The authors developed Newmark-β sensitivities in [7,46], which are in turn based on [33,34]. In this paper, Newmark-β serves as a comparison to the developed method for generalized-α.

This work first generally introduces generalized-α time integration and its equations. Afterwards, the sensitivity equations are introduced for structural dynamics and flexible multibody dynamics, respectively. A numerical example is shown for each and the results are provided for both the primal and sensitivity analyses in comparison with Newmark-β time integration. We close with a conclusion, including a discussion of this work. Furthermore, a detailed derivation of the effective system matrices is provided in the Appendix A and Appendix B.

2. Sensitivity Analysis with Generalized-α Time Integration

In this section, we introduce the equations for an acceleration-form (a-form) of generalized-α time integration using a predictor–corrector scheme and extend this for sensitivity analysis. The acceleration-form refers to the formulation, where the acceleration is the primary variable (i.e., the variable for which it is being solved), and contrasts with the displacement-form (d-form) in which the displacement is the primary variable. The predictor–corrector scheme allows for nonlinear terms in the governing equation in which the predicted values are updated in the nonlinear solver, e.g., Newton–Raphson, and then corrected.

In this work, we extend the well-established generalized-α time integration method to include design sensitivity analysis, which is defined by the differentiation with respect to the design variables (but could include other parameters, e.g., uncertain parameters). Sensitivity analysis is carried out here by directly differentiating both the equations of the time integration as well as the governing equations. We consider both linear structural dynamics as well as flexible multibody dynamics, which is intrinsically nonlinear. After introducing the general equations for generalized-α time integration and its sensitivities, the domain-specific governing equations for linear structural dynamics and flexible multibody dynamics are given. Flexible multibody dynamics is shown using the specific approach developed in [7,47,48] with an index-1 approach and the floating frame of reference formulation.

2.1. Generalized-$\alpha$ Time Integration

The method for primal and sensitivity analysis with generalized-α time integration is introduced in the following. This is based on Newmark’s equations for the updated position $\underline{q}$ and velocity vectors $\underline{\dot{q}}$ with the calculated value of the acceleration vector $\underline{\ddot{q}}$,

$$\begin{array}{cc}\hfill {\underline{q}}_{n+1}& =\underset{\mathrm{predictor}{\underline{q}}_{\mathrm{pred}}}{\underbrace{{\underline{q}}_n+\Delta t{\underline{\dot{q}}}_n+\left(\frac{1}{2}-\beta \right)\Delta t^2{\underline{\ddot{q}}}_n}}+\underset{\mathrm{corrector}}{\underbrace{\beta \Delta t^2{\underline{\ddot{q}}}_{n+1}}},\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\underline{\dot{q}}}_{n+1}& =\underset{\mathrm{predictor}{\underline{\dot{q}}}_{\mathrm{pred}}}{\underbrace{{\underline{\dot{q}}}_n+\left(1-\gamma \right)\Delta t{\underline{\ddot{q}}}_n}}+\underset{\mathrm{corrector}}{\underbrace{\gamma \Delta t{\underline{\ddot{q}}}_{n+1}}},\hfill \end{array}$$

(2)

where $\gamma$ and $\beta$ are the Newmark integration parameters. Though the simulation is evaluated at the times $t_n$ and $t_{n+1}$, values at intermediate points are interpolated. The position q and velocity $\dot{q}$ use time $t_{n+1-{\alpha }_f}$ and the acceleration $\ddot{q}$ uses time $t_{n+1-{\alpha }_m}$, giving the approximations by linear interpolation (see Figure 1),

Figure 1. Time integration with — acceleration, — $\gamma =0.5$, — $\gamma =0.75$.

Figure 1. Time integration with — acceleration, — $\gamma =0.5$, — $\gamma =0.75$.

$$\begin{array}{cc}\hfill {\underline{q}}_{n+1-{\alpha }_f}& =\left(1-{\alpha }_f\right){\underline{q}}_{n+1}+{\alpha }_f{\underline{q}}_n,\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\underline{\dot{q}}}_{n+1-{\alpha }_f}& =\left(1-{\alpha }_f\right){\underline{\dot{q}}}_{n+1}+{\alpha }_f{\underline{\dot{q}}}_n,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill {\underline{\ddot{q}}}_{n+1-{\alpha }_m}& =\left(1-{\alpha }_m\right){\underline{\ddot{q}}}_{n+1}+{\alpha }_m{\underline{\ddot{q}}}_n,\hfill \end{array}$$

(5)

where ${\alpha }_m$ and ${\alpha }_f$ are integration parameters. A predictor–corrector scheme is used in order to accommodate the nonlinear calculation and state-dependency of the system parameters, i.e., the system matrices with respect to position and velocity. The predicted state variables with Newmark’s equations are

$$\begin{array}{cc}\hfill \underline{q}{ }_{\mathrm{pred}}& ={\underline{q}}_n+\Delta t{\underline{\dot{q}}}_n+\left(\frac{1}{2}-\beta \right)\Delta t^2{\underline{\ddot{q}}}_n,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \underline{\dot{q}}{ }_{\mathrm{pred}}& ={\underline{\dot{q}}}_n+\left(1-\gamma \right)\Delta t{\underline{\ddot{q}}}_n,\hfill \end{array}$$

(7)

and we calculate ${\underline{\ddot{q}}}_{n+1}$ using the effective system equations (see below for an explanation of structural dynamics and flexible multibody dynamics, respectively). The predicted values are then corrected with Newmark’s equations,

$$\begin{array}{cc}\hfill {\underline{q}}_{n+1}& ={\underline{q}}_{\mathrm{pred}}+\beta \Delta t^2{\underline{\ddot{q}}}_{n+1},\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill {\underline{\dot{q}}}_{n+1}& ={\underline{\dot{q}}}_{\mathrm{pred}}+\gamma \Delta t{\underline{\ddot{q}}}_{n+1}.\hfill \end{array}$$

(9)

For the sensitivity analysis within the generalized-α time integration, we use the direct differentiation of the approximations

$$\begin{array}{cc}\hfill \underline{\underline{\nabla q}}{ }_{n+1-{\alpha }_f}& =\left(1-{\alpha }_f\right)\underline{\underline{\nabla q}}{ }_{n+1}+{\alpha }_f\underline{\underline{\nabla q}}{ }_n,\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \underline{\underline{\nabla \dot{q}}}{ }_{n+1-{\alpha }_f}& =\left(1-{\alpha }_f\right)\underline{\underline{\nabla \dot{q}}}{ }_{n+1}+{\alpha }_f\underline{\underline{\nabla \dot{q}}}{ }_n,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \underline{\underline{\nabla \ddot{q}}}{ }_{n+1-{\alpha }_m}& =\left(1-{\alpha }_m\right)\underline{\underline{\nabla \ddot{q}}}{ }_{n+1}+{\alpha }_m\underline{\underline{\nabla \ddot{q}}}{ }_n,\hfill \end{array}$$

(12)

where the nabla operator ∇ denotes the design sensitivity, i.e., the total derivative with respect to the design variables, e.g., for the sensitivity scalar response with respect to the design variable vector,

$$\underline{\nabla \left(·\right)}=\frac{\mathrm{d}\left(·\right)}{\mathrm{d}\underline{x}}.$$

(13)

The same direct differentiation is applied to the predictors, which are then defined by

$$\begin{array}{cc}\hfill \underline{\underline{\nabla q}}{ }_{\mathrm{pred}}& =\underline{\underline{\nabla q}}{ }_n+\Delta t\underline{\underline{\nabla \dot{q}}}{ }_n+\left(\frac{1}{2}-\beta \right)\Delta t^2\underline{\underline{\nabla \ddot{q}}}{ }_n,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \underline{\underline{\nabla \dot{q}}}{ }_{\mathrm{pred}}& =\underline{\underline{\nabla \dot{q}}}{ }_n+\left(1-\gamma \right)\Delta t\underline{\underline{\nabla \ddot{q}}}{ }_n.\hfill \end{array}$$

(15)

After calculating $\underline{\underline{\nabla \ddot{q}}}{ }_{n+1}$ by use of the effective sensitivity system (see below for explanation of structural dynamics and flexible multibody dynamics, respectively), the following correctors are used:

$$\begin{array}{cc}\hfill \underline{\underline{\nabla q}}{ }_{n+1}& =\underline{\underline{\nabla q}}{ }_{\mathrm{pred}}+\beta \Delta t^2\underline{\underline{\nabla \ddot{q}}}{ }_{n+1},\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill \underline{\underline{\nabla \dot{q}}}{ }_{n+1}& =\underline{\underline{\nabla \dot{q}}}{ }_{\mathrm{pred}}+\gamma \Delta t\underline{\underline{\nabla \ddot{q}}}{ }_{n+1}.\hfill \end{array}$$

(17)

For a nonlinear analysis, a Newton solver can be integrated into the solving routine. This is explained with flexible multibody systems below in Section 2.3.

The time integration constants of generalized-α are ${\alpha }_m$, ${\alpha }_f$, $\beta$ and $\gamma$ and 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 generally used 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, and the spectral radius at infinity ${\varrho }_{\infty }$ [13],

$$\begin{array}{cc}\hfill {\alpha }_m& =\frac{2{\varrho }_{\infty }-1}{{\varrho }_{\infty }+1},\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\alpha }_f& =\frac{{\varrho }_{\infty }}{{\varrho }_{\infty }+1},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \gamma & =\frac{1}{2}-{\alpha }_m+{\alpha }_f,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \beta & =\frac{1}{4}{\left(1-{\alpha }_m+{\alpha }_f\right)}^2.\hfill \end{array}$$

(21)

The spectral radius at infinity is sometimes referred to as the (numerical) damping and is limited to a range of $\left[0,1\right]$. The value ${\varrho }_{\infty }=1$ returns the Newmark-β retaining frequencies, which are preserved, while ${\varrho }_{\infty }=0$ annihilates all frequencies above $\frac{\Delta t}{T}$ after one time step in which T is the period of the highest frequency of interest. The results converge to those without numerical damping when the time step goes to zero. Thus, the level of numerical damping introduced to generalized-α time integration depends on both the step size as well as the spectral radius at infinity and the choice is critical for proper results. A parameter study with a fixed-time step and varying the values of the spectral radius at infinity is shown in [14]. As the present study develops a general framework for this methodology, this is not addressed here and the reader is referred to [13,49,50,51]. In this paper, we use a moderate value of ${\varrho }_{\infty }=0.55$, though the sensitivity analysis methodology developed is general and independent of this parameter value.

2.2. Structural Dynamics

The governing equation for linear structural dynamics is defined by the following second-order equation of motion:

$$\begin{array}{cc}\hfill \underline{\underline{m}}\underline{\ddot{q}}\left(t\right)+\underline{\underline{d}}\underline{\dot{q}}\left(t\right)+\underline{\underline{k}}\underline{q}\left(t\right)-\underline{F}\left(t\right)& =\underline{0},\hfill \end{array}$$

(22)

where $\underline{\underline{m}}$ is the mass matrix, $\underline{\underline{d}}$ is the damping matrix, $\underline{\underline{k}}$ is the stiffness matrix and $\underline{F}$ is the external force vector. For linear structural dynamics, the system matrices are constant and independent of position and time. For the solving scheme, we add the residual force vector $\underline{R}\left(t\right)$, which is to be zero (to an allowable value) at all times, to the right-hand side. The α-family of integration methods uses intermediate-point approximations between $t_n$ and $t_{n+1}$, namely $t_{n+1-{\alpha }_f}$ and $t_{n+1-{\alpha }_m}$, revealing

$${\underline{R}}_{n+1-{\alpha }_f}=\underline{\underline{m}}{\underline{\ddot{q}}}_{n+1-{\alpha }_m}+\underline{\underline{d}}{\underline{\dot{q}}}_{n+1-{\alpha }_f}+\underline{\underline{k}}{\underline{q}}_{n+1-{\alpha }_f}-{\underline{F}}_{n+1-{\alpha }_f},$$

(23)

where the linear approximation of the intermediate external force is

$${\underline{F}}_{n+1-{\alpha }_f}=\left(1-{\alpha }_f\right){\underline{F}}_{n+1}+{\alpha }_f{\underline{F}}_n.$$

(24)

The effective equation of motion is defined by

$$\underline{\underline{m}}{ }_{\mathrm{eff}}\ddot{\underline{q}}{ }_{n+1}=\underline{F}{ }_{\mathrm{eff}},$$

(25)

where

$$\begin{array}{cc}\hfill {\underline{\underline{m}}}_{\mathrm{eff}}=& \left(1-{\alpha }_m\right)\underline{\underline{m}}+\left(1-{\alpha }_f\right)\gamma \Delta t\underline{\underline{d}}+\left(1-{\alpha }_f\right)\beta \Delta t^2\underline{\underline{k}},\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill {\underline{F}}_{\mathrm{eff}}=& \left(1-{\alpha }_f\right){\underline{F}}_{n+1}+{\alpha }_f{\underline{F}}_n+\hfill \\ \hfill & -{\alpha }_m\underline{\underline{m}}{\underline{\ddot{q}}}_n-\underline{\underline{d}}\left(\left(1-{\alpha }_f\right){\dot{\underline{q}}}_{\mathrm{pred}}+{\alpha }_f{\dot{\underline{q}}}_n\right)-\underline{\underline{k}}\left(\left(1-{\alpha }_f\right){\underline{q}}_{\mathrm{pred}}+{\alpha }_f{\underline{q}}_n\right).\hfill \end{array}$$

(27)

For the derivation of the effective mass matrix ${\underline{\underline{m}}}_{\mathrm{eff}}$ and effective force vector ${\underline{F}}_{\mathrm{eff}}$ for generalized-α method, see Appendix A.1.

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

$$\begin{array}{cc}\hfill \underline{\underline{\nabla R}}{ }_{n+1-{\alpha }_f}=& \underline{\underline{\underline{\nabla m}}}\underline{\ddot{q}}{ }_{n+1-{\alpha }_m}+\underline{\underline{m}}\underline{\underline{\nabla \ddot{q}}}{ }_{n+1-{\alpha }_m}+\underline{\underline{\underline{\nabla d}}}{\underline{\dot{q}}}_{n+1-{\alpha }_f}+\underline{\underline{d}}\underline{\underline{\nabla \dot{q}}}{ }_{n+1-{\alpha }_f}+\hfill \\ \hfill & +\underline{\underline{\underline{\nabla k}}}\underline{q}{ }_{n+1-{\alpha }_f}+\underline{\underline{k}}\underline{\underline{\nabla q}}{ }_{n+1-{\alpha }_f}-\underline{\underline{\nabla F}}{ }_{n+1-{\alpha }_f},\hfill \end{array}$$

(28)

where the intermediate external force sensitivity is found with the approximation

$$\underline{\underline{\nabla F}}{ }_{n+1-{\alpha }_f}=\left(1-{\alpha }_f\right)\underline{\underline{\nabla F}}{ }_{n+1}+{\alpha }_f\underline{\underline{\nabla F}}{ }_n.$$

(29)

The design sensitivities of the system matrices are three-dimensional matrices and therefore efficient handling is needed.

The effective equation of motion for sensitivities is defined by the differentiation of (25),

$$\begin{array}{cc}\hfill \underline{\underline{m}}{ }_{\mathrm{eff}}\underline{\underline{\nabla \ddot{q}}}{ }_{n+1}& =\underline{\underline{F}}{ }_{\mathrm{pseudo}},\hfill \end{array}$$

(30)

where

$$\underline{\underline{F}}{ }_{\mathrm{pseudo}}=\underline{\underline{\nabla F}}{ }_{\mathrm{eff}}-\underline{\underline{\underline{\nabla m}}}{ }_{\mathrm{eff}}\ddot{\underline{q}}{ }_{n+1},$$

(31)

and fully expressed by

$$\begin{array}{cc}\hfill \underline{\underline{F}}{ }_{\mathrm{pseudo}}=& \left(1-{\alpha }_f\right)\underline{\underline{\nabla F}}{ }_{n+1}+{\alpha }_f\underline{\underline{\nabla F}}{ }_n-\underline{\underline{\underline{\nabla m}}}\left(\left(1-{\alpha }_m\right){\underline{\ddot{q}}}_{n+1}+{\alpha }_m{\underline{\ddot{q}}}_n\right)-{\alpha }_m\underline{\underline{m}}\underline{\underline{\nabla \ddot{q}}}{ }_n+\hfill \\ \hfill & -\underline{\underline{\underline{\nabla d}}}\left(\left(1-{\alpha }_f\right){\underline{\dot{q}}}_{n+1}+{\alpha }_f{\underline{\dot{q}}}_n\right)-\underline{\underline{d}}\left(\left(1-{\alpha }_f\right)\underline{\underline{\nabla \dot{q}}}{ }_{\mathrm{pred}}+{\alpha }_f\underline{\underline{\nabla \dot{q}}}{ }_n\right)+\hfill \\ \hfill & -\underline{\underline{\underline{\nabla k}}}\left(\left(1-{\alpha }_f\right){\underline{q}}_{n+1}+{\alpha }_f{\underline{q}}_n\right)-\underline{\underline{k}}\left(\left(1-{\alpha }_f\right)\underline{\underline{\nabla q}}{ }_{\mathrm{pred}}+{\alpha }_f\underline{\underline{\nabla q}}{ }_n\right).\hfill \end{array}$$

(32)

For the derivation of the effective mass matrix ${\underline{\underline{m}}}_{\mathrm{eff}}$ and effective pseudo force $\underline{\underline{F}}{ }_{\mathrm{pseudo}}$ for the generalized-α method, see Appendix A.1 and Appendix A.2.

2.3. Flexible Multibody Dynamics

Flexible multibody dynamics is the analysis of mechanical systems consisting of flexible bodies constrained by kinematic joints modeled with differential–algebraic equations. This work builds upon the previous work of analytical sensitivity analysis with multibody dynamics, including [37,52,53]. In addition to the kinematic constraints with respect to structural dynamics, multibody dynamics also must integrate a nonlinear solver in order to consider the dependencies of the system matrices and vectors on the state variables. This is carried out with Newton steps.

The flexible multibody system is modeled with the floating frame of the reference formulation [54,55,56]. The generic index-3 differential–algebraic equation for a flexible multibody system with holonomic constraints is written as follows:

$$\begin{array}{cc}\hfill \underline{\underline{m}}\left(\underline{q},t\right)\underline{\ddot{q}}\left(t\right)+\underline{\underline{d}}\left(\underline{q},t\right)\underline{\dot{q}}\left(t\right)+\underline{\underline{k}}\left(\underline{q},t\right)\underline{q}\left(t\right)+\underline{\underline{\mathrm{J}}}{ }_{\Phi }^T\left(\underline{q},t\right)\underline{\lambda }\left(t\right)& ={\underline{F}}_{\mathrm{ext}}\left(\underline{q},\underline{\dot{q}},t\right)+{\underline{F}}_v\left(\underline{q},\underline{\dot{q}},t\right),\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill \underline{\Phi }\left(\underline{q},t\right)& =\underline{0},\hfill \end{array}$$

(34)

where $\underline{q}$, $\underline{\dot{q}}$ and $\underline{\ddot{q}}$ are the generalized position, velocity and acceleration vectors; $\underline{\lambda }$ is the vector of Lagrange multipliers of kinematic constraints; the matrices $\underline{\underline{m}}$, $\underline{\underline{d}}$, and $\underline{\underline{k}}$ are the mass, damping and stiffness matrices; $\underline{\underline{\mathrm{J}}}{ }_{\Phi }$ is the Jacobian matrix of kinematic constraints; ${\underline{F}}_{\mathrm{ext}}$ is the vector of external forces; ${\underline{F}}_v$ is the quadratic velocity vector; and $\underline{\Phi }$ is the vector of kinematic constraints at the position level.

Numerical problems related to the solution of index-3 differential–algebraic equations [57] can be avoided by a number of methods including index reduction, coordinate reduction, augmented Lagrangian formulation and preconditioning [58,59,60]. Although generalized-α time integration is suitable for directly solving index-3 or index-2 differential–algebraic equations, in this work, we address the problem using index reduction leading to an index-1 differential–algebraic equation in which the kinematic constraint function is differentiated twice with respect to time,

$$\begin{array}{cc}\hfill \underline{\underline{m}}\underline{\ddot{q}}+\underline{\underline{d}}\underline{\dot{q}}+\underline{\underline{k}}\underline{q}+\underline{\underline{\mathrm{J}}}{ }_{\Phi }^T\underline{\lambda }& ={\underline{F}}_{\mathrm{ext}}+{\underline{F}}_v,\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill \underline{\underline{\mathrm{J}}}{ }_{\Phi }\underline{\ddot{q}}& ={\underline{F}}_c,\hfill \end{array}$$

(36)

where ${\underline{F}}_c$ is the right-hand side of acceleration constraints. As index-1 differential–algebraic equations are susceptible to drift in the kinematic constraints, Baumgarte stabilization [61] is used. Kinematic constraint drift is stabilized in this method by adding damping and restoring terms to the acceleration constraint equations that are proportional to the velocity and position constraints. Details of its implementation, including the sensitivity analysis, are described in [14,62].

The generalized-α time integration is based on

$$\begin{array}{cc}\hfill {\left(\underline{\underline{m}}\underline{\ddot{q}}\right)}_{n+1-{\alpha }_m}+{\left(\underline{\underline{d}}\underline{\dot{q}}\right)}_{n+1-{\alpha }_f}+{\left(\underline{\underline{k}}\underline{q}\right)}_{n+1-{\alpha }_f}+{\left(\underline{\underline{\mathrm{J}}}{ }_{\Phi }^T\underline{\lambda }\right)}_{n+1-{\alpha }_m}& ={\underline{F}}_{\mathrm{ext},n+1-{\alpha }_f}+{\underline{F}}_{v,n+1-{\alpha }_f},\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill {\left(\underline{\underline{\mathrm{J}}}{ }_{\Phi }\underline{\ddot{q}}\right)}_{n+1-{\alpha }_m}& ={\underline{F}}_{c,n+1-{\alpha }_f},\hfill \end{array}$$

(38)

with the intermediate approximations of all force terms of the equations of motion. These intermediate approximations are shown in Appendix 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

$$\begin{array}{cc}\hfill {\underline{R}}_{1,n+1-{\alpha }_f}=& \left(1-{\alpha }_m\right){\underline{\underline{m}}}_{n+1}{\underline{\ddot{q}}}_{n+1}+{\alpha }_m{\underline{\underline{m}}}_n{\underline{\ddot{q}}}_n+\left(1-{\alpha }_f\right){\underline{\underline{d}}}_{n+1}{\underline{\dot{q}}}_{n+1}+{\alpha }_f{\underline{\underline{d}}}_n{\underline{\dot{q}}}_n+\hfill \\ \hfill & +\left(1-{\alpha }_f\right){\underline{\underline{k}}}_{n+1}{\underline{q}}_{n+1}+{\alpha }_f{\underline{\underline{k}}}_n{\underline{q}}_n+\left(1-{\alpha }_m\right)\underline{\underline{\mathrm{J}}}{ }_{\Phi ,n+1}^T{\underline{\lambda }}_{n+1}+{\alpha }_m\underline{\underline{\mathrm{J}}}{ }_{\Phi ,n}^T{\underline{\lambda }}_n+\hfill \\ \hfill & -\left(1-{\alpha }_f\right){\underline{F}}_{\mathrm{ext},n+1}-{\alpha }_f{\underline{F}}_{\mathrm{ext},n}-\left(1-{\alpha }_f\right){\underline{F}}_{v,n+1}-{\alpha }_f{\underline{F}}_{v,n},\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill {\underline{R}}_{2,n+1-{\alpha }_f}=& \left(1-{\alpha }_m\right)\underline{\underline{\mathrm{J}}}{ }_{\Phi ,n+1}{\underline{\ddot{q}}}_{n+1}+{\alpha }_m\underline{\underline{\mathrm{J}}}{ }_{\Phi ,n}{\underline{\ddot{q}}}_n-\left(1-{\alpha }_f\right){\underline{F}}_{c,n+1}-{\alpha }_f{\underline{F}}_{c,n}.\hfill \end{array}$$

(40)

In the case of a nonlinear mass matrix, as is the case with the floating frame of reference formulation, this weighted formulation of the residual equations are first-order accurate. For a formulation of second-order accuracy, see, e.g., [63], where the dynamic equilibrium is enforced exactly at the time steps of the simulation. The derived time integration and its differentiation with the weighted residual formulation can be extended and applied to further formulations.

The effective equation of motion is defined by

$$\left[\begin{array}{cc}{\underline{\underline{m}}}_{\mathrm{eff}}& \underline{\underline{\mathrm{J}}}{ }_{\Phi ,\mathrm{eff}}^T\\ \underline{\underline{\mathrm{J}}}{ }_{\Phi ,\mathrm{eff}}& \underline{\underline{0}}\end{array}\right]\left[\begin{array}{c}{\underline{\ddot{q}}}_{n+1}\\ {\underline{\lambda }}_{n+1}\end{array}\right]=\left[\begin{array}{c}{\underline{F}}_{a,\mathrm{eff}}\\ {\underline{F}}_{c,\mathrm{eff}}\end{array}\right],$$

(41)

where

$$\begin{array}{cc}\hfill {\underline{\underline{m}}}_{\mathrm{eff}}=& \left(1-{\alpha }_m\right){\underline{\underline{m}}}_{n+1}+\left(1-{\alpha }_f\right)\gamma \Delta t{\underline{\underline{d}}}_{n+1}+\left(1-{\alpha }_f\right)\beta \Delta t^2{\underline{\underline{k}}}_{n+1},\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill \underline{\underline{\mathrm{J}}}{ }_{\Phi ,\mathrm{eff}}=& \left(1-{\alpha }_m\right)\underline{\underline{\mathrm{J}}}{ }_{\Phi ,n+1},\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill {\underline{F}}_{a,\mathrm{eff}}=& \left(1-{\alpha }_f\right){\underline{F}}_{\mathrm{ext},n+1}+{\alpha }_f{\underline{F}}_{\mathrm{ext},n}+\left(1-{\alpha }_f\right){\underline{F}}_{v,n+1}+{\alpha }_f{\underline{F}}_{v,n}+\hfill \\ \hfill & -{\alpha }_m{\underline{\underline{m}}}_n{\underline{\ddot{q}}}_n-\left(1-{\alpha }_f\right){\underline{\underline{d}}}_{n+1}{\underline{\dot{q}}}_{\mathrm{pred}}-{\alpha }_f{\underline{\underline{d}}}_n{\underline{\dot{q}}}_n+\hfill \\ \hfill & -\left(1-{\alpha }_f\right){\underline{\underline{k}}}_{n+1}{\underline{q}}_{\mathrm{pred}}-{\alpha }_f{\underline{\underline{k}}}_n{\underline{q}}_n-{\alpha }_m\underline{\underline{\mathrm{J}}}{ }_{\Phi ,n}^T{\underline{\lambda }}_n,\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill {\underline{F}}_{c,\mathrm{eff}}=& \left(1-{\alpha }_f\right){\underline{F}}_{c,n+1}+{\alpha }_f{\underline{F}}_{c,n}-{\alpha }_m\underline{\underline{\mathrm{J}}}{ }_{\Phi ,n}{\underline{\ddot{q}}}_n.\hfill \end{array}$$

(45)

The derivations of the effective equation of motion and the terms ${\underline{\underline{m}}}_{\mathrm{eff}}$, $\underline{\underline{\mathrm{J}}}{ }_{\Phi ,\mathrm{eff}}$, ${\underline{F}}_{a,\mathrm{eff}}$ and ${\underline{F}}_{c,\mathrm{eff}}$ are shown in Appendix B.1.

With flexible multibody dynamics, the system parameters (matrix- and vector-valued for $n_{\mathrm{dof}}>1$) including mass, damping, stiffness, constraints and forces are generally dependent on state variables. Therefore, a nonlinear solver is needed to solve the system and the Newton–Raphson method is used here of the following form:

$$\left[\begin{array}{cc}\frac{\partial {\underline{R}}_{1,n+1-{\alpha }_f}}{\partial {\underline{\ddot{q}}}_{n+1}}& \frac{\partial {\underline{R}}_{1,n+1-{\alpha }_f}}{\partial {\underline{\lambda }}_{n+1}}\\ \frac{\partial {\underline{R}}_{2,n+1-{\alpha }_f}}{\partial {\underline{\ddot{q}}}_{n+1}}& \frac{\partial {\underline{R}}_{2,n+1-{\alpha }_f}}{\partial {\underline{\lambda }}_{n+1}}\end{array}\right]\left[\begin{array}{c}\Delta {\underline{\ddot{q}}}_{n+1}\\ \Delta {\underline{\lambda }}_{n+1}\end{array}\right]+\left[\begin{array}{c}{\underline{R}}_{1,n+1-{\alpha }_f}\\ {\underline{R}}_{2,n+1-{\alpha }_f}\end{array}\right]=\underline{0}.$$

(46)

With the dependencies shown in (33) and (34), and the terms of the Jacobian matrix (tangent operator) of the residuals with respect to ${\underline{\ddot{q}}}_{n+1}$ and ${\underline{\lambda }}_{n+1}$ are given by

$$\begin{array}{cc}\hfill \frac{\partial {\underline{R}}_{1,n+1-{\alpha }_f}}{\partial {\underline{\ddot{q}}}_{n+1}}=& \left(1-{\alpha }_m\right)\left(\beta \Delta t^2\frac{\partial {\underline{\underline{m}}}_{n+1}}{\partial {\underline{q}}_{n+1}}{\underline{\ddot{q}}}_{n+1}+{\underline{\underline{m}}}_{n+1}\right)+\hfill \\ \hfill & +\left(1-{\alpha }_f\right)\left(\beta \Delta t^2\frac{\partial {\underline{\underline{d}}}_{n+1}}{\partial {\underline{q}}_{n+1}}{\underline{\dot{q}}}_{n+1}+\gamma \Delta t{\underline{\underline{d}}}_{n+1}\right)+\hfill \\ \hfill & +\left(1-{\alpha }_f\right)\left(\beta \Delta t^2\frac{\partial {\underline{\underline{k}}}_{n+1}}{\partial {\underline{q}}_{n+1}}{\underline{q}}_{n+1}+\beta \Delta t^2{\underline{\underline{k}}}_{n+1}\right)+\hfill \\ \hfill & +\left(1-{\alpha }_m\right)\beta \Delta t^2\frac{\partial \underline{\underline{\mathrm{J}}}{ }_{\Phi ,n+1}^T}{\partial {\underline{q}}_{n+1}}{\underline{\lambda }}_{n+1}+\hfill \\ \hfill & -\left(1-{\alpha }_f\right)\left(\beta \Delta t^2\frac{\partial {\underline{F}}_{\mathrm{ext},n+1}}{\partial {\underline{q}}_{n+1}}+\gamma \Delta t\frac{\partial {\underline{F}}_{\mathrm{ext},n+1}}{\partial {\underline{\dot{q}}}_{n+1}}\right)+\hfill \\ \hfill & -\left(1-{\alpha }_f\right)\left(\beta \Delta t^2\frac{\partial {\underline{F}}_{v,n+1}}{\partial {\underline{q}}_{n+1}}+\gamma \Delta t\frac{\partial {\underline{F}}_{v,n+1}}{\partial {\underline{\dot{q}}}_{n+1}}\right),\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill \frac{\partial {\underline{R}}_{2,n+1-{\alpha }_f}}{\partial {\underline{\ddot{q}}}_{n+1}}=& \left(1-{\alpha }_m\right)\left(\beta \Delta t^2\frac{\partial \underline{\underline{\mathrm{J}}}{ }_{\Phi ,n+1}}{\partial {\underline{q}}_{n+1}}{\underline{\ddot{q}}}_{n+1}+\underline{\underline{\mathrm{J}}}{ }_{\Phi ,n+1}\right)+\hfill \\ \hfill & -\left(1-{\alpha }_f\right)\left(\beta \Delta t^2\frac{\partial {\underline{F}}_{c,n+1}}{\partial {\underline{q}}_{n+1}}+\gamma \Delta t\frac{\partial {\underline{F}}_{c,n+1}}{\partial {\underline{\dot{q}}}_{n+1}}\right),\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill \frac{\partial {\underline{R}}_{1,n+1-{\alpha }_f}}{\partial {\underline{\lambda }}_{n+1}}=& \left(1-{\alpha }_m\right)\underline{\underline{\mathrm{J}}}{ }_{\Phi ,n+1}^T,\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill \frac{\partial {\underline{R}}_{2,n+1-{\alpha }_f}}{\partial {\underline{\lambda }}_{n+1}}=& \underline{\underline{0}}.\hfill \end{array}$$

(50)

The derivation of these terms is shown in Appendix B.2.

The direct differentiation of the equations of motion (39) and (40) leads to the governing equation of the sensitivity analysis,

$$\begin{array}{cc}\hfill \underline{\underline{\nabla R}}{ }_{1,n+1-{\alpha }_f}=& \left(1-{\alpha }_m\right){\underline{\underline{\underline{\nabla m}}}}_{n+1}{\underline{\ddot{q}}}_{n+1}+{\alpha }_m{\underline{\underline{\underline{\nabla m}}}}_n{\underline{\ddot{q}}}_n+\hfill \\ \hfill & +\left(1-{\alpha }_m\right){\underline{\underline{m}}}_{n+1}{\underline{\underline{\nabla \ddot{q}}}}_{n+1}+{\alpha }_m{\underline{\underline{m}}}_n{\underline{\underline{\nabla \ddot{q}}}}_n+\hfill \\ \hfill & +\left(1-{\alpha }_f\right){\underline{\underline{\underline{\nabla d}}}}_{n+1}{\underline{\dot{q}}}_{n+1}+{\alpha }_f{\underline{\underline{\underline{\nabla d}}}}_n{\underline{\dot{q}}}_n+\hfill \\ \hfill & +\left(1-{\alpha }_f\right){\underline{\underline{d}}}_{n+1}{\underline{\underline{\nabla \dot{q}}}}_{n+1}+{\alpha }_f{\underline{\underline{d}}}_n{\underline{\underline{\nabla \dot{q}}}}_n+\hfill \\ \hfill & +\left(1-{\alpha }_f\right){\underline{\underline{\underline{\nabla k}}}}_{n+1}{\underline{q}}_{n+1}+{\alpha }_f{\underline{\underline{\underline{\nabla k}}}}_n{\underline{q}}_n+\hfill \\ \hfill & +\left(1-{\alpha }_f\right){\underline{\underline{k}}}_{n+1}{\underline{\underline{\nabla q}}}_{n+1}+{\alpha }_f{\underline{\underline{k}}}_n{\underline{\underline{\nabla q}}}_n+\hfill \\ \hfill & +\left(1-{\alpha }_m\right)\underline{\underline{\underline{\nabla \mathrm{J}}}}{ }_{\Phi ,n+1}^T{\underline{\lambda }}_{n+1}+{\alpha }_m\underline{\underline{\underline{\nabla \mathrm{J}}}}{ }_{\Phi ,n}^T{\underline{\lambda }}_n+\hfill \\ \hfill & +\left(1-{\alpha }_m\right)\underline{\underline{\mathrm{J}}}{ }_{\Phi ,n+1}^T{\underline{\underline{\nabla \lambda }}}_{n+1}+{\alpha }_m\underline{\underline{\mathrm{J}}}{ }_{\Phi ,n}^T{\underline{\underline{\nabla \lambda }}}_n+\hfill \\ \hfill & -\left(1-{\alpha }_f\right){\underline{\underline{\nabla F}}}_{\mathrm{ext},n+1}-{\alpha }_f{\underline{\underline{\nabla F}}}_{\mathrm{ext},n}-\left(1-{\alpha }_f\right){\underline{\underline{\nabla F}}}_{v,n+1}-{\alpha }_f{\underline{\underline{\nabla F}}}_{v,n},\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill \underline{\underline{\nabla R}}{ }_{2,n+1-{\alpha }_f}=& \left(1-{\alpha }_m\right)\underline{\underline{\underline{\nabla \mathrm{J}}}}{ }_{\Phi ,n+1}{\underline{\ddot{q}}}_{n+1}+{\alpha }_m\underline{\underline{\underline{\nabla \mathrm{J}}}}{ }_{\Phi ,n}{\underline{\ddot{q}}}_n\hfill \\ \hfill & +\left(1-{\alpha }_m\right)\underline{\underline{\mathrm{J}}}{ }_{\Phi ,n+1}{\underline{\underline{\nabla \ddot{q}}}}_{n+1}+{\alpha }_m\underline{\underline{\mathrm{J}}}{ }_{\Phi ,n}{\underline{\underline{\nabla \ddot{q}}}}_n+\hfill \\ \hfill & -\left(1-{\alpha }_f\right){\underline{\underline{\nabla F}}}_{c,n+1}-{\alpha }_f{\underline{\underline{\nabla F}}}_{c,n}.\hfill \end{array}$$

(52)

The governing sensitivity equation of the effective system is defined by the differentiation of (41),

$$\left[\begin{array}{cc}{\underline{\underline{m}}}_{\mathrm{eff}}& \underline{\underline{\mathrm{J}}}{ }_{\Phi ,\mathrm{eff}}^T\\ \underline{\underline{\mathrm{J}}}{ }_{\Phi ,\mathrm{eff}}& \underline{\underline{0}}\end{array}\right]\left[\begin{array}{c}{\underline{\underline{\nabla \ddot{q}}}}_{n+1}\\ {\underline{\underline{\nabla \lambda }}}_{n+1}\end{array}\right]=\left[\begin{array}{c}{\underline{\underline{F}}}_{a,\mathrm{pseudo}}\\ {\underline{\underline{F}}}_{c,\mathrm{pseudo}}\end{array}\right],$$

(53)

where

$$\left[\begin{array}{c}{\underline{\underline{F}}}_{a,\mathrm{pseudo}}\\ {\underline{\underline{F}}}_{c,\mathrm{pseudo}}\end{array}\right]=\left[\begin{array}{c}{\underline{\underline{\nabla F}}}_{a,\mathrm{eff}}\\ {\underline{\underline{\nabla F}}}_{c,\mathrm{eff}}\end{array}\right]-\left[\begin{array}{cc}\underline{\underline{\underline{\nabla m}}}{ }_{\mathrm{eff}}& \underline{\underline{\underline{\nabla \mathrm{J}}}}{ }_{\Phi ,\mathrm{eff}}^T\\ \underline{\underline{\underline{\nabla \mathrm{J}}}}{ }_{\Phi ,\mathrm{eff}}& \underline{\underline{\underline{0}}}\end{array}\right]\left[\begin{array}{c}{\underline{\ddot{q}}}_{n+1}\\ {\underline{\lambda }}_{n+1}\end{array}\right],$$

(54)

leading to

$$\begin{array}{cc}\hfill {\underline{\underline{F}}}_{a,\mathrm{pseudo}}=& \left(1-{\alpha }_f\right){\underline{\underline{\nabla F}}}_{\mathrm{ext},n+1}+{\alpha }_f{\underline{\underline{\nabla F}}}_{\mathrm{ext},n}+\left(1-{\alpha }_f\right){\underline{\underline{\nabla F}}}_{v,n+1}+{\alpha }_f{\underline{\underline{\nabla F}}}_{v,n}+\hfill \\ \hfill & -\left(1-{\alpha }_m\right){\underline{\underline{\underline{\nabla m}}}}_{n+1}{\underline{\ddot{q}}}_{n+1}-{\alpha }_m{\underline{\underline{\underline{\nabla m}}}}_n{\underline{\ddot{q}}}_n-{\alpha }_m{\underline{\underline{m}}}_n{\underline{\underline{\nabla \ddot{q}}}}_n+\hfill \\ \hfill & -\left(1-{\alpha }_f\right){\underline{\underline{\underline{\nabla d}}}}_{n+1}{\underline{\dot{q}}}_{n+1}-{\alpha }_f{\underline{\underline{\underline{\nabla d}}}}_n{\underline{\dot{q}}}_n-\left(1-{\alpha }_f\right){\underline{\underline{d}}}_{n+1}{\underline{\underline{\nabla \dot{q}}}}_{\mathrm{pred}}-{\alpha }_f{\underline{\underline{d}}}_n{\underline{\underline{\nabla \dot{q}}}}_n+\hfill \\ \hfill & -\left(1-{\alpha }_f\right){\underline{\underline{\underline{\nabla k}}}}_{n+1}{\underline{q}}_{n+1}-{\alpha }_f{\underline{\underline{\underline{\nabla k}}}}_n{\underline{q}}_n-\left(1-{\alpha }_f\right){\underline{\underline{k}}}_{n+1}{\underline{\underline{\nabla q}}}_{\mathrm{pred}}-{\alpha }_f{\underline{\underline{k}}}_n{\underline{\underline{\nabla q}}}_n+\hfill \\ \hfill & -\left(1-{\alpha }_m\right)\underline{\underline{\underline{\nabla \mathrm{J}}}}{ }_{\Phi ,n+1}^T{\underline{\lambda }}_{n+1}-{\alpha }_m\underline{\underline{\underline{\nabla \mathrm{J}}}}{ }_{\Phi ,n}^T{\underline{\lambda }}_n-{\alpha }_m\underline{\underline{\mathrm{J}}}{ }_{\Phi ,n}^T{\underline{\underline{\nabla \lambda }}}_n,\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill {\underline{\underline{F}}}_{c,\mathrm{pseudo}}=& \left(1-{\alpha }_f\right){\underline{\underline{\nabla F}}}_{c,n+1}+{\alpha }_f{\underline{\underline{\nabla F}}}_{c,n}+\hfill \\ \hfill & -\left(1-{\alpha }_m\right)\underline{\underline{\underline{\nabla \mathrm{J}}}}{ }_{\Phi ,n+1}{\underline{\ddot{q}}}_{n+1}-{\alpha }_m\underline{\underline{\underline{\nabla \mathrm{J}}}}{ }_{\Phi ,n}{\underline{\ddot{q}}}_n-{\alpha }_m\underline{\underline{\mathrm{J}}}{ }_{\Phi ,n}{\underline{\underline{\nabla \ddot{q}}}}_n.\hfill \end{array}$$

(56)

The derivation of the effective mass matrix and the effective constraint Jacobian matrix is shown in Appendix B.1 and the derivation of the effective pseudo force terms is shown in Appendix B.3.

The Newton–Raphson step for the sensitivity analysis is given by

$$\left[\begin{array}{cc}\frac{\partial {\underline{\underline{\nabla R}}}_{1,n+1-{\alpha }_f}}{\partial {\underline{\underline{\nabla \ddot{q}}}}_{n+1}}& \frac{\partial {\underline{\underline{\nabla R}}}_{1,n+1-{\alpha }_f}}{\partial {\underline{\underline{\nabla \lambda }}}_{n+1}}\\ \frac{\partial {\underline{\underline{\nabla R}}}_{2,n+1-{\alpha }_f}}{\partial {\underline{\underline{\nabla \ddot{q}}}}_{n+1}}& \frac{\partial {\underline{\underline{\nabla R}}}_{2,n+1-{\alpha }_f}}{\partial {\underline{\underline{\nabla \lambda }}}_{n+1}}\end{array}\right]\left[\begin{array}{c}\Delta {\underline{\underline{\nabla \ddot{q}}}}_{n+1}\\ \Delta {\underline{\underline{\nabla \lambda }}}_{n+1}\end{array}\right]+\left[\begin{array}{c}{\underline{\underline{\nabla R}}}_{1,n+1-{\alpha }_f}\\ {\underline{\underline{\nabla R}}}_{2,n+1-{\alpha }_f}\end{array}\right]=\underline{\underline{0}}.$$

(57)

As shown in [7], the Jacobian of the residual sensitivity with respect to ${\underline{\underline{\nabla \ddot{q}}}}_{n+1}$ and ${\underline{\underline{\nabla \lambda }}}_{n+1}$ reduces to

$$\left[\begin{array}{cc}\frac{\partial {\underline{\underline{\nabla R}}}_{1,n+1-{\alpha }_f}}{\partial {\underline{\underline{\nabla \ddot{q}}}}_{n+1}}& \frac{\partial {\underline{\underline{\nabla R}}}_{1,n+1-{\alpha }_f}}{\partial {\underline{\underline{\nabla \lambda }}}_{n+1}}\\ \frac{\partial {\underline{\underline{\nabla R}}}_{2,n+1-{\alpha }_f}}{\partial {\underline{\underline{\nabla \ddot{q}}}}_{n+1}}& \frac{\partial {\underline{\underline{\nabla R}}}_{2,n+1-{\alpha }_f}}{\partial {\underline{\underline{\nabla \lambda }}}_{n+1}}\end{array}\right]=\left[\begin{array}{cc}\frac{\partial {\underline{R}}_{1,n+1-{\alpha }_f}}{\partial {\underline{\ddot{q}}}_{n+1}}& \frac{\partial {\underline{R}}_{1,n+1-{\alpha }_f}}{\partial {\underline{\lambda }}_{n+1}}\\ \frac{\partial {\underline{R}}_{2,n+1-{\alpha }_f}}{\partial {\underline{\ddot{q}}}_{n+1}}& \frac{\partial {\underline{R}}_{2,n+1-{\alpha }_f}}{\partial {\underline{\lambda }}_{n+1}}\end{array}\right]\underline{\underline{\underline{\underline{e}}}},$$

(58)

which reduces the dimensions of the coefficient matrix from $\left(n_q+n_{\Phi }\right)\times n_x\times \left(n_q+n_{\Phi }\right)\times n_x$ to $\left(n_q+n_{\Phi }\right)\times \left(n_q+n_{\Phi }\right)$. The solution of the resulting equation,

$$\left[\begin{array}{cc}\frac{\partial {\underline{R}}_{1,n+1-{\alpha }_f}}{\partial {\underline{\ddot{q}}}_{n+1}}& \frac{\partial {\underline{R}}_{1,n+1-{\alpha }_f}}{\partial {\underline{\lambda }}_{n+1}}\\ \frac{\partial {\underline{R}}_{2,n+1-{\alpha }_f}}{\partial {\underline{\ddot{q}}}_{n+1}}& \frac{\partial {\underline{R}}_{2,n+1-{\alpha }_f}}{\partial {\underline{\lambda }}_{n+1}}\end{array}\right]\left[\begin{array}{c}\Delta {\underline{\underline{\nabla \ddot{q}}}}_{n+1}\\ \Delta {\underline{\underline{\nabla \lambda }}}_{n+1}\end{array}\right]+\left[\begin{array}{c}{\underline{\underline{\nabla R}}}_{1,n+1-{\alpha }_f}\\ {\underline{\underline{\nabla R}}}_{2,n+1-{\alpha }_f}\end{array}\right]=\underline{\underline{0}},$$

(59)

is found using the same solving scheme for the nonlinear solver of the primal analysis (46). Furthermore, both the primal analysis (46) and sensitivity analysis (59) have the same coefficient matrix, which also allows one to reuse the coefficient matrix of the primal analysis in the sensitivity analysis without the need for further evaluations.

3. Numerical Results

In this section, numerical results are shown comparing the introduced method for the sensitivity analysis for generalized-α time integration (primal analysis introduced in Section 2) with that of Newmark-β. The methods are compared with applications to both structural dynamics as well as flexible multibody dynamics.

3.1. Structural Dynamics

The linear structural dynamic problem used here as a demonstration example is based on [64], which has found usage as a benchmark in [65,66,67,68,69,70,71]. This unitless problem is schematically shown in Figure 2. This problem is defined by

$$\stackrel{\underline{\underline{m}}}{\overbrace{\left[\begin{array}{ccc}{m}_1& 0& 0\\ 0& m_2& 0\\ 0& 0& m_3\end{array}\right]}}\stackrel{\underline{\ddot{q}}}{\overbrace{\left[\begin{array}{c}{\ddot{q}}_1\\ {\ddot{q}}_2\\ {\ddot{q}}_3\end{array}\right]}}+\stackrel{\underline{\underline{k}}}{\overbrace{\left[\begin{array}{ccc}{k}_1& -k_1& 0\\ -k_1& k_1+k_2& -k_2\\ 0& -k_2& k_2\end{array}\right]}}\stackrel{\underline{q}}{\overbrace{\left[\begin{array}{c}{q}_1\\ q_2\\ q_3\end{array}\right]}}-\stackrel{\underline{F}}{\overbrace{\left[\begin{array}{c}{F}_1\\ 0\\ 0\end{array}\right]}}=\stackrel{\underline{R}}{\overbrace{\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]}},$$

(60)

with the prescribed displacement,

$$q_1=sin\omega t,$$

(61)

and the unitless values

$$\begin{array}{cc}\hfill m_1& =0,\hfill \\ \hfill m_2& =1,\hfill \\ \hfill m_3& =1,\hfill \\ \hfill k_1& ={10}^7,\hfill \\ \hfill k_2& =1,\hfill \\ \hfill \omega & =1.2.\hfill \end{array}$$

The masses $m_2$ and $m_3$ as well as the stiffnesses $k_1$ and $k_2$ are considered design variables with respect to which the design sensitivities are calculated. The problem consists of stiff ($k_1)$ and flexible ($k_2)$ parts that are typical of complex, practical structural–mechanical systems [64].

Considering the prescribed displacement, this problem can be reduced to two degrees of freedom,

$$\stackrel{\underline{\underline{m}}}{\overbrace{\left[\begin{array}{cc}{m}_2& 0\\ 0& m_3\end{array}\right]}}\stackrel{\underline{\ddot{q}}}{\overbrace{\left[\begin{array}{c}{\ddot{q}}_2\\ {\ddot{q}}_3\end{array}\right]}}+\stackrel{\underline{\underline{k}}}{\overbrace{\left[\begin{array}{cc}{k}_1+k_2& -k_2\\ -k_2& k_2\end{array}\right]}}\stackrel{\underline{q}}{\overbrace{\left[\begin{array}{c}{q}_2\\ q_3\end{array}\right]}}-\stackrel{\underline{F}}{\overbrace{\left[\begin{array}{c}{k}_1{q}_1\\ 0\end{array}\right]}}=\stackrel{\underline{R}}{\overbrace{\left[\begin{array}{c}0\\ 0\end{array}\right]}},$$

(62)

which we implement in our in-house, Python-based solver framework SiMuLi in the time $t=\left[0,10\right]$ with the initial conditions

$$\begin{array}{cc}\hfill q_2& =0,\hfill \\ \hfill q_3& =0,\hfill \\ \hfill {\dot{q}}_2& =0,\hfill \\ \hfill {\dot{q}}_3& =0.\hfill \end{array}$$

As reported in [64], unitless time steps of $\Delta t=0.2618$ are used, giving the following values for the ratio step size to the fundamental periods of the system:

$$\begin{array}{cc}\hfill \frac{\Delta t}{T_F}& =0.05,\hfill \\ \hfill \frac{\Delta t}{T_1}& =0.041667,\hfill \\ \hfill \frac{\Delta t}{T_2}& =131.76,\hfill \end{array}$$

where $T_F$ is the fundamental period of the forcing function and $T_{1,2}$, for the first and second modes, $T_1=6.283$ and $T_2=0.001987$.

The results for both Newmark-β and generalized-α time integration methods are shown in Figure 3 and compared with the exact solution. The exact solution is calculated via mode superposition per [72]. The numerical damping parameter of generalized-α time integration is set as ${\varrho }_{\infty }=0.55$. As expected with generalized-α time integration, the unphysical, high-frequent oscillatory behavior of the response can be avoided. As such, the beneficial effects of numerical damping with generalized-α time integration can be clearly seen in the velocity and acceleration of node 2, while maintaining the integrity of the other responses.

The sensitivity analysis developed above is applied and these results can be seen in Figure 4, Figure 5 and Figure 6. Analytical sensitivity analysis with Newmark-β and generalized-α time integration methods are also compared with their numerical sensitivity implementations via finite differencing using a relative perturbation of ${10}^{-6}$.

In cases where the Newmark-β time integration shows instability or noise, the numerical damping allows for usable sensitivities. In some cases, noise—albeit at a much lower level—persists in the first iterations before being damped out—after approximately $2.5\mathrm{s}$ (or approximately 10 iterations). Attention is drawn to those sensitivities susceptible to stability issues, namely

• sensitivities of position 2 with respect to stiffness 1 and mass 2 (Figure 4);
• sensitivities of velocity 2 with respect to stiffness 1 and mass 2 (Figure 5);
• sensitivities of acceleration 2 with respect to stiffness 1, stiffness 2 and mass 2 (Figure 6);
• sensitivities of acceleration 3 with respect to stiffness 1 and mass 2 (Figure 6).

Generalized-α time integration is not afflicted by this instability, resulting in sensitivities that can be successfully used in gradient-based algorithms.

Furthermore, generalized-α time integration eliminates the noise in those sensitivities present when using Newmark-β. This is seen with the following sensitivity responses:

• sensitivities of velocity 2 with respect to stiffness 2 (Figure 5);
• sensitivities of velocity 3 with respect to stiffness 1 and mass 2 (Figure 5).

The analytical and numerical results are, in most plots, nearly exact. There are cases, however, where the negative consequences of the numerical instability with Newmark-β are increased. The reader’s attention is especially called to

• sensitivity of velocity 2 with respect to stiffness 2 (Figure 5);
• sensitivity of velocity 2 with respect to mass 3 (Figure 5);
• sensitivity of acceleration 2 with respect to stiffness 2 (Figure 6);
• sensitivity of acceleration 2 with respect to mass 3 (Figure 6).

It can be seen that both analytical and numerical sensitivities with generalized-α time integration are not negatively affected by numerical instabilities, as is the case with Newmark-β.

The number of time steps needed for the damping to take effect can be reduced by reducing the numerical damping parameter ${\varrho }_{\infty }$ towards zero (i.e., more damping), though also affecting the physical response to a greater degree. As time continues, the instability grows, with Newmark-β time integration and therefore a clear advantage can be seen with generalized-α. It can be clearly seen how gradient-based solvers would struggle to use the gradient information provided with Newmark-β time integration.

For our modest case of four design variables, we show an acceleration of two, where the speedup S is defined by the ratio of computational time for finite differencing $t_{\mathrm{FD}}$ to that of analytical sensitivities $t_{\mathrm{ana}}$,

$$S=\frac{t_{\mathrm{FD}}}{t_{\mathrm{ana}}}.$$

(63)

The computational time required for analytical sensitivity analysis represents approximately the ideal computational effort twice that of the primal analysis, i.e., analysis without sensitivity analysis. Computation reduction is the same for both Newmark-β and generalized-α time integration methods, as can be seen in

Table 1. Computational effort for the structural dynamics example.

Time Integration	Sensitivity Analysis	Calculation Time $\left(\mathbf{s}\right){ }^1$	Speedup ${ }^2$
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

${ }^1$ Intel Core i5-7200U CPU 2.50 GHz, ${ }^2$ acceleration of analytical sensitivity analysis with respect to finite differencing.

.

3.2. Flexible Slider–Crank Mechanism

The slider–crank mechanism, as shown in Figure 7, is modeled with flexible multibody dynamics using the floating frame of reference formulation. The mechanism consists of four bodies that are connected by three revolute joints and one prismatic joint, modeled as ideal joints without friction. The ground and the slider are considered to be rigid, while the crank and the rod are considered to be flexible and are modeled with Euler–Bernoulli beams. Figure 8 shows the flexible multibody system of the mechanism. The inertial reference frame coinciding with the reference frame of the ground is shown in blue. The body reference frames of the crank, the rod and the slider are shown in red. The undeformed bodies of the crank and the rod are shown in gray and the deformed bodies are shown in black, where the points represent the nodes of the finite elements for the geometric discretization. The coordinates that are investigated in detail for the primal analysis and the sensitivity analysis are shown in green. The dimensions of the bodies are listed in

Table 2. Geometry of the 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

. The bodies are made of steel and

Table 3. Material properties.

Property	Symbol	Value	Units
Density	$\rho$	$7.85\times {10}^{-9}$	$\frac{\mathrm{t}}{\mathrm{m}{\mathrm{m}}^3}$
Elastic modulus	E	210,000	$MPa$
Poisson ratio	$\nu$	0.3	−

shows the material properties. A linear elastic material model with isotropic behavior is applied. In this example, the SI-derived MPa unit system is used with length in millimeters (mm), mass in tonnes (t), time in second (s), force in Newtons (N), Angle in radian (rad) and stress in megapascals (MPa).

The generalized coordinates of a flexible multibody system modeled with the floating frame of the reference formulation and Euler–Bernoulli beams are given by

$$\begin{array}{cc}\hfill \underline{q}& ={\left[\begin{array}{ccc}{\underline{r}}^T& \theta & {\overline{\underline{q}}}_f^T\end{array}\right]}^T,\hfill \end{array}$$

(64)

where $\underline{r}$ is the position of the body reference frame, $\theta$ is the orientation of the body reference frame and ${\overline{\underline{q}}}_f$ are the nodal displacements in the body coordinates of each node on the flexible body. Since no damping is considered, the equations of motion of a flexible body are given by

$$\begin{array}{cc}\hfill \underline{\underline{m}}\underline{\ddot{q}}+\underline{\underline{k}}\underline{q}+\underline{\underline{\mathrm{J}}}{ }_{\Phi }^T\underline{\lambda }& ={\underline{F}}_{\mathrm{ext}}+{\underline{F}}_v,\hfill \end{array}$$

(65)

$$\begin{array}{cc}\hfill \underline{\underline{\mathrm{J}}}{ }_{\Phi }\underline{\ddot{q}}& ={\underline{F}}_c,\hfill \end{array}$$

(66)

where the system matrices of the flexible bodies are given by

$$\begin{array}{cc}\hfill \underline{\underline{m}}& =\left[\begin{array}{ccc}{\underline{\underline{m}}}_{tt}& {\underline{\underline{m}}}_{tr}& {\underline{\underline{m}}}_{tf}\\ & {\underline{\underline{m}}}_{rr}& {\underline{\underline{m}}}_{rf}\\ \mathrm{sym}.& & {\underline{\underline{m}}}_{ff}\end{array}\right],\hfill \end{array}$$

(67)

$$\begin{array}{cc}\hfill \underline{\underline{k}}& =\left[\begin{array}{ccc}0& 0& 0\\ & 0& 0\\ \mathrm{sym}.& & {\underline{\underline{k}}}_{ff}\end{array}\right],\hfill \end{array}$$

(68)

and the system vectors are given by

$$\begin{array}{cc}\hfill {\underline{F}}_{\mathrm{ext}}& =\left[\begin{array}{c}{\underline{F}}_{\mathrm{ext},t}\\ {\underline{F}}_{\mathrm{ext},r}\\ {\underline{F}}_{\mathrm{ext},f}\end{array}\right],\hfill \end{array}$$

(69)

$$\begin{array}{cc}\hfill {\underline{F}}_v& =\left[\begin{array}{c}{\underline{F}}_{v,t}\\ {\underline{F}}_{v,r}\\ {\underline{F}}_{v,f}\end{array}\right],\hfill \end{array}$$

(70)

with the index denoting translational terms t, rotational terms r and flexible terms f.

For the dynamic simulation of the mechanism, the above-mentioned solver framework SiMuLi is used. The simulation is performed in the time interval $t=\left[0\mathrm{s},0.025\mathrm{s}\right]$. For both generalized-α and Newmark-β time integration, a time step of $\Delta t=1\times {10}^{-4}s$ is chosen for comparison. The resulting values for the ratio step size to fundamental periods of the system for the first ten non-zero eigenfrequencies are found in Figure 9. The values for the ratio step size to fundamental periods range from $0.2119\mathrm{s}$ for first non-zero eigenfrequency to $3.9919\mathrm{s}$ for the tenth non-zero eigenfrequency and to $17.1157\mathrm{s}$ for the twentieth non-zero eigenfrequency (not pictured in the plot).

For the initial conditions, the angular position and velocity of the crank are defined with

$$\begin{array}{cc}\hfill {\theta }_C& =0\mathrm{r}\mathrm{a}\mathrm{d},\hfill \\ \hfill {\dot{\theta }}_C& =150\pi \frac{\mathrm{r}\mathrm{a}\mathrm{d}}{\mathrm{s}},\hfill \end{array}$$

and the initial conditions of the other bodies are derived through the kinematic analysis of the system. The dimensions of the cross-sections of the crank $h_C$, $w_C$ and the rod $h_R$, $w_R$ are used as design variables in the sensitivity analysis.

The results of the simulation of the coordinates that are highlighted in Figure 8 are shown in Figure 10. With Newmark-β, vibrations at high frequencies are observed. With generalized-α, vibrations are observed in the beginning of the simulation that are numerically damped with ongoing simulation time. In addition to the generalized coordinates, the maximum stresses on the upper and lower fiber of the beam elements of the flexible bodies including the crank and the rod are computed. The Kreisselmeier–Steinhauser function [73,74] is used to approximate the maximum stress value of the bodies at each time step with a differentiable function and the results are shown in Figure 11. With Newmark-β, the vibrating behavior of the system causes higher flexible deformation and therefore unrealistic high stress in the components. The numerical damping with generalized-α leads to lower maximum stress values and a more realistic stress distribution over time.

The results of the sensitivity analysis of the coordinates that are highlighted in Figure 8 are shown in Figure 12, Figure 13 and Figure 14. The analyzed coordinates include the rotation of the floating frame of the crank ${\theta }_C$, the translation in the x direction of the floating frame of the slider $x_S$ and the three flexible coordinates of the node in A on the crank ${\overline{x}}_{C,f}$, ${\overline{y}}_{C,f}$ and ${\overline{\theta }}_{C,f}$. The design variables used for the sensitivity analysis are the cross-sectional dimensions of the crank $h_C$, $w_C$ and the rod $h_R$, $w_R$.

Similarly to the above, with Newmark-β, the results show a vibrating behavior that would lead to significant problems when using the sensitivities in gradient-based design optimization. With generalized-α, the numerical damping leads to usable sensitivity values after a few time steps. Figure 12 shows that the position sensitivities of the position and orientation coordinates of the floating frames are less noisy compared to the flexible coordinates.

Position sensitivities that suffer little from noise are limited to the position and orientation coordinates of the floating frames with respect to several design variables and include the following:

• sensitivities of ${\theta }_C$ with respect to $w_C$ and $w_R$;
• sensitivities of $x_S$ with respect to $h_C$, $w_C$, and $w_R$.

Position sensitivities with significant noise include the position and orientation coordinates of the floating frames with respect to a few specific design variables as well as the flexible coordinates with respect to all design variables. The list of noisy responses of position sensitivities is given in the following:

• sensitivities of ${\theta }_C$ with respect to $h_C$ and $h_R$;
• sensitivities of $x_S$ with respect to $h_R$;
• sensitivities of ${\overline{x}}_{C,f}$ with respect to $h_C$ $w_C$, $h_R$ and $w_R$;
• sensitivities of ${\overline{y}}_{C,f}$ with respect to $h_C$; $w_C$, $h_R$ and $w_R$;
• sensitivities of ${\overline{\theta }}_{C,f}$ with respect to $h_C$, $w_C$, $h_R$ and $w_R$.

The sensitivities of the velocities as shown in Figure 13 and the sensitivities of the accelerations as shown in Figure 14 have noisy responses with Newmark that are damped with generalized-α for all investigated coordinates and all used design variables. The same behavior also continues with the sensitivity values for the maximum stresses in the crank and the rod shown in Figure 15.

To check the correctness of the results, the sensitivity values with the introduced method are compared to numerical sensitivity analysis with forward differencing. The results are in good agreement, with only small numerical differences. In fact, there are no visible differences in Figure 12, Figure 13, Figure 14 and Figure 15 with the finite differencing values lying under the respective analytical values.

Table 4. Computational effort for the slider–crank mechanism.

Time Integration	Sensitivity Analysis	Calculation Time [mm:ss] ${ }^1$	Speedup ${ }^2$
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

${ }^1$ Intel Core i7-8700 CPU 3.20 GHz, ² speedup of analytical sensitivity analysis with respect to finite differencing.

shows a comparison of the computational effort with different methods. The computational effort with generalized-α is about $30\%$ lower compared to the computational effort with Newmark-β. The main reason for that is the continuing vibrating behavior with Newmark-β requires more iteration steps of the nonlinear solver. The comparison of the sensitivity methods with numerical finite differencing and analytical direct differentiation shows the high efficiency of the analytical method. With the four shown design variables, as shown in (63), an acceleration of about $S\approx 4$ of the analytical method compared to the numerical method is obtained. The additional computational effort of analytical sensitivity analysis is just a fraction of the computational effort of the primal analysis. This is because of the simplification shown in (59), which also allows one to reuse the Jacobian matrix of the residual built in the primal analysis for the sensitivity analysis.

4. Conclusions

This work develops the analytical sensitivity analysis for the α-family of time integration, which has a single implementation for several common time integration schemes. This therefore has consequences for gradient-based optimization frameworks for dynamic systems. Analytical sensitivity analysis is crucial for efficient design optimization, parameter fitting and uncertainty analysis. In this work, the direct differentiation of the generalized-α time integration is derived and implemented for linear structural dynamics and flexible multibody dynamics, the latter of which requires a nonlinear solution routine. All equations needed for implementation are provided and extendable to further governing equations. The beneficial effects of the step-size-dependent numerical damping provided by the generalized-α method can clearly be seen in the comparison of the results with Newmark-β. The numerical noise of the responses and sensitivities can lead to convergence problems when using gradient-based optimization algorithms. This can be avoided when using the generalized-α time integration and its sensitivities as developed herein.

If the functions of interest are nonsmooth, discontinuous and bifurcated, non-gradient-based methods can be successfully used, e.g., design cases with crash analysis. Such cases can be solved with gradient-free optimization algorithms (e.g., [75,76,77]) and approximation methods (e.g., [78,79,80]).

The system responses with respect to the design variables are indeed often smooth, continuous and unbifurcated and thus enabling gradient-based methods where efficient and accurate sensitivity analysis is critical for both effectiveness and efficiency. This is shown above with linear structural dynamics and flexible multibody dynamics. The method and equations developed herein are of general nature and are extendable to further formulations and applicable to larger problems, both of larger number of degrees of freedom and number of design variables. Future investigations are needed on these large-scale problems to assess the acceleration and include a comparison with the adjoint variable method.

With the above study and introduction of the direct differentiation of the α-family of time integrators, the authors hope to provide the readers as design engineers a tool in the pursuit of optimal structures and mechanisms under dynamic loading.