A semi-analytical approach to sensitivity analysis with flexible multibody dynamics of a morphing forward wing section

In the present paper, semi-analytical design sensitivity analysis of a morphing forward wing section modeled as a flexible multibody system is developed. The flexible multibody dynamics model consists of a flexible external wing skin and an actuation mechanism with rigid bodies. The floating frame of reference formulation is used for the formulation of the equations of motion, while time integration is carried out with the generalized-α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\alpha$\end{document} method and Newton–Raphson iterations for the nonlinear analysis. These steps are considered in both the primal analysis and the developed efficient design sensitivity analysis. The calculation of the design sensitivities is based on direct differentiation and carried out with a semi-analytical approach. The responses of interest include the deviations of the wing profile from the targeted morphed geometry and the resulting stresses in the wing skin. The design variables considered in this work are geometric parameters, material parameters and loading parameters. The calculations with the shown method provide reliable and accurate results combined with high computational efficiency. The visualization of the computed sensitivity values gives an easy interpretation of the sensitivities and facilitates the understanding of the design engineer on how to change the design variables to improve the design of the system. The introduced sensitivity analysis enables future investigations, including the application of the method to gradient-based design optimization and uncertainty analysis.

here are the stress and the geometric deviations from the target geometry. These are differentiated with respect to the design variables, which include geometric parameters, material parameters and loading parameters. The method is general and not limited by those design variables explored here. Future application fields of the developed sensitivity method include gradient-based design optimization and uncertainty analysis.

Flexible multibody dynamics
Flexible multibody simulation is comprised of three components [8,11]: governing equations, time integration and nonlinear solver. The main simulation to calculate the responses is referred as primal analysis to differentiate from the design sensitivity analysis, in which the sensitivity of the responses with respect to design variables are calculated. The simulation routine is shown as a flowchart in Fig. 1 and includes the components, which are introduced in § 2 for the primal analysis and in § 3 for the design sensitivity analysis.

Governing equation
The governing equations of flexible multibody systems are given by the equations of motion and the kinematic constraint equations, which lead to a set of index-3 differential-algebraic equations given by where q is the vector of generalized positions, λ is the vector of Lagrangian multipliers of the kinematic constraints, m is the mass matrix, d is the damping matrix, k is the stiffness Fig. 2 Floating frame of reference formulation (FFRF) for a two-dimensional beam with the global reference frame in blue, the body reference frame in red, the undeformed beam in gray and the deformed beam in black. (Color figure online) matrix, Φ is the vector of kinematic constraints, J Φ is the Jacobian matrix of the constraints (i.e. the partial derivative of the constraints with respect to position), F ext is the external force vector and F v is the quadratic velocity force vector. Overdots represent the first( ·) and second( ·) temporal derivatives. Symbols with single underline (·) denote vectors, symbols with double underline (·) denote two-dimensional matrices, symbols with triple underline denote three-dimensional terms, symbols with quadruple underline denote four-dimensional terms and those without underlines are scalars. Further, symbols with an overline (·) are expressed in floating coordinates of the body reference frame. To reduce numerical problems related to index-3 differential-algebraic equations [12,13], the kinematic constraint equations are differentiated twice with respect to time to obtain a system of index-1 differential-algebraic equations expressed in matrix form by where R is the residual vector and F c is the right-hand side of the kinematic constraints in acceleration form. As index-1 differential-algebraic equations are affected by drift, Baumgarte stabilization [14] is used as described by [15]. FMBD is reviewed by [13,[16][17][18]. Available flexible multibody formulations include the absolute nodal coordinate formulation (ANCF) [19], the absolute coordinate formulation (ACF) [20], the equivalent rigid link system (ERLS) [21,22], geometrically exact beam formulation [23,24] and the floating frame of reference formulation (FFRF) [25,26]. The lattermost is used here, which is suitable for beam, shell and solid FE elements and can be used with a linear-elastic material model as well as with geometric or material nonlinearities [13,27].
In the following derivation of the equations of motion for FFRF, the index of each body is omitted for a simpler notation. The equations are shown for the general three-dimensional case. The generalized positions in terms of FFRF are paramterized by the position r and orientation β of the reference frame and flexible deformations in local coordinates q f , see Fig. 2, The continuous position of a material point P on a flexible body can be seen in Fig. 2 and is described in global coordinates as with the transformation matrix A and the local position vector u P that is decomposed by the undeformed term u P ,o and the deformed term u P ,f . In this implementation, FFRF is used in combination with finite-element analysis. Each flexible body is discretized and the flexible deformations in local coordinates are given by with the matrix of shape functions S expressed in floating coordinates and the nodal deformations q f of the finite-element mesh. A linear-elastic material model is used, which leads to a FFRF stiffness matrix k that is linear, with the finite-element stiffness matrix k ff . In contrast, the FFRF mass matrix m is highly nonlinear and defined by with the volume V , the density ρ, the identity matrix e and the matrix B. It should be noted that in the system matrices of FFRF, t denotes the translatory degrees of freedom of the body frame, r the rotational degrees of freedom of the body frame and f the flexible degrees of freedom. Accordingly, the index ff represents the system matrices from the finite-element model. The matrix B is defined by whereũ P is the skew-symmetric matrix of positions in local coordinates and the matrix G relates the angular velocity vector of the floating frame ω and the time derivatives of the orientation parametersβ, The quadratic velocity vector is highly nonlinear and defined by In FFRF, the external forces are assembled in a vector as The vectors and matrices of each body are assembled to obtain the system vectors and matrices. The governing equations are to be fulfilled and hence are solved at all time steps.

Time integration
The system of differential-algebraic equations shown in § 2 is solved in time with the generalized-α time integration method originally introduced by [28] and implemented as a predictor-corrector scheme, as described in [8,11]. Generalized-α time integration is based on Newmark's equations [29], with intermediate approximations of the generalized-α method, The resulting effective system of equations for FMBD is where the components are defined as The position and velocity values are updated with the solved accelerations and continue to the next time step until the termination time is reached.

Nonlinear solver
The residual equations with generalised-α time integration are Newton-Raphson iterations (i.e. exact Newton method) are used to consider nonlinearity of the system by solving with the Jacobians w.r.t. accelerationJ and Lagrangian multipliers λ J giving the expressions for the Jacobians of the residual byJ The components of the Jacobian of the residual arë For a complete derivation of the terms above using generalized-α time integration including time integration, the readers are referred to [11].

Design sensitivity analysis
Although referred to as design sensitivities, the sensitivities of the system responses with respect to parameters are useful in a wide range of applications. Such applications include design optimization (including sizing, shape and topology optimization), uncertainty analysis as well as the direct application of the sensitivities to assess robustness and parameter influence. This paper investigates discrete sensitivity analysis, also referred as discretize-thendifferentiate, where the differentiation w.r.t. the design variables is carried out on discretized governing equations [30]. These include temporal discretization into time steps of numerical time integration and spatial discretization into finite elements of the flexible bodies. This is in contrast to continuous sensitivity analysis, also known as the variational approach or differentiate-then-discretize, where continuum governing equations are differentiated before discretization.
To avoid the high computational effort and lack of precision of numerical sensitivity analysis, analytical methods including the direct differentiation method [13,31] and the adjoint variable method [32][33][34] are preferred for the sensitivity analysis of flexible multibody systems. When the number of design variables is higher than the number of optimization responses including objective and constraint functions, the adjoint variable method is generally more efficient. In the reverse case, where the number of design variables is less than the number of optimization responses, the direct differentiation method is more suitable. The direct differential method is chosen here and integrated into the primal solve routine to avoid the solving backwards in time as needed with an adjoint variable methodology (therefore often referred to the backward method). Related challenges to the adjoint variable method include the handling of time-dependent system matrices (save and reload or recalculate), especially when using variable time steps and a required special solver for the adjoint system.
The design sensitivity analysis is performed here with direct differentiation using a semianalytical approach. The differentiation is carried out through all three steps of the calculation routine: the governing equation, the sensitivity analysis and the nonlinear solver (cf. Fig. 1). The introduced method is similar to the staggered corrector method introduced by [35] with the backward difference formula. This method originates from the staggered direct scheme, where an iteration loop of the primal analysis is followed by the iteration loop of the sensitivities for each time step [36]. Efficiency improvements compared to the staggered direct scheme have been obtained by the simultaneous corrector scheme, where primal and sensitivity analysis are solved simultaneously with one iteration loop [37]. With the staggered corrector method as used here, further efficiency improvements are obtained using two sequential iteration loops for primal and sensitivity analysis and by reducing the updates of the residual Jacobians [35].

Governing equation
The direct differentiation of the governing equations of the primal analysis (3) results in turn with the governing equations for the sensitivity analysis, where ∇ is the partial derivative with respect to the design variables, It should be noted that (38) is of the same form as (3) and accordingly the same solving routine can be used. The pseudo load F pseudo is comprised of the partial derivatives of the system parameters with respect to the design variables, This method is shown in [9], where the partial derivatives of (40) are approximated numerically with forward differencing, This thus results in a semi-analytical approach, where the governing equations are derived analytically, but the partial derivatives of these equations are computed numerically. In contrast, the partial derivatives of those parameters dependent on the finite-element model m, k and F v are differentiated analytically by a further layer based on an invariant approach, see [10]. This is done to avoid loading the finite-element model at each time step and leads to a further decrease in computation effort compared to the method shown in [9].

Time integration
The sensitivity analysis with time integration is analogous to the time integration of the primal analysis. For the generalised-α method, we utilize a predictor-corrector scheme with Newmark's equations and intermediate approximations [8,11]. The resulting effective sensitivity system is where the pseudo load case is The sensitivity values of position and velocity are carried out from the sensitivity values of acceleration analogously to the primal analysis.

Nonlinear solver
The Jacobian matrix of the sensitivity analysis simplifies to the Jacobian matrix of the primal analysis [8,11], allowing to reuse the already computed Jacobian of the primal analysis for an efficient computation method of the sensitivity analysis, with the Jacobian with respect to acceleration sensitivity ∇J and with respect to the Lagrangian multipliers sensitivity ∇ λ J. With this simplification, the Jacobian of the sensitivity analysis is reduced from four dimensions to two. This saves excessive memory usage and thus represents the key to efficient sensitivity analysis of multibody systems with direct differentiation.

Morphing wing model
In this work, the developed method is applied to the design of a morphing wing concept developed in [2,3,38,39]. Specifically, the geometry of the morphing wing sailplane for both  undeformed (high-speed) and morphed (low-speed) configurations is introduced in [2,40]. The airfoil shapes are designed using a numerical optimization approach where aerodynamic lift and drag are considered [2,40]. The geometries of the airfoils along with the aircraft for which they were designed are shown in Fig. 3. A morphing concept shows a significant performance advantage over conventional sailplanes with a camber changing flap: 20% higher lift coefficients are achieved at equal drag coefficients compared to conventional camber changing trailing edge flap laminar airfoils, thus allowing to reduce the wing area and decrease profile drag by the same ratio. This results in 20 km h higher interthermal dash speed.
In [5], the morphing actuation of this concept is achieved using compliant mechanisms designed with topology optimization. This concept uses a stack of six individual compliant mechanisms and is shown in Fig. 4. In contrast to a compliant concept, we will address a traditional hinged mechanism to drive the deformation of the wing in this work.
The forward wing section is modeled as flexible multibody system using FFRF. Figure 5a shows the planar system actuated with a conventional hinged mechanism and the two-dimensional case here is the simplification of the three-dimensional case shown in § 2. The flexible outer shell (body 1), i.e. wing skin, is modeled as flexible body using planar Euler-Bernoulli beams with a FE mesh consisting of 31 nodes and 30 elements. For this system, a wing section with the length of 1000 mm and a wall thickness of 1.5 mm is considered. On the upper end of the forward wing section (point F), the wing skin is fully constrained in all three degrees of freedom, i.e. x, y and θ . On the lower end of the forward wing section (point G), the wing skin is constrained in two degrees of freedom y and θ , allowing for free motion in x.
The morphing forward wing section is actuated by a hinged mechanism consisting of four rigid bodies 2, 3, 4 and 5. The motion of these bodies is constrained by six joints including five revolute joints in the points A, B, C, D and E and one rigid joint between the bodies 2 and 4 in the point B. For actuation, a torque is applied to body 2 on point A. The torque value is zero at the initial time and increases with the curve in Fig. 6 to a maximum value of 4.2 N · m.
A generic aluminum alloy is used with a density ρ of 2.7×10 −9 t/mm 3 , a Young modulus E of 70 000 MPa and a Poisson ratio ν of 0.35 [−]. The material is assumed to be isotropic and to remain in the linear elastic regime.
The transient solution of the flexible multibody system is obtained with generalised-α time integration and Newton-Raphson iterations as nonlinear solver. The system is solved at each time step Δt = 0.05 s. The deformation of the external wing skin and the motion of the rigid bodies of the multibody system from the low-lift (high-speed, undeformed) configuration to the high-lift (low-speed, morphed) configuration is shown in Fig. 5b.
The key for a proper design of the morphing forward wing section is the interaction between the actuation mechanism and the flexible external wing skin. Firstly, it is important that the deformation of the morphed forward wing section approximates the target wing profile as closely as possible. Secondly, the material limits of the wing skin must not be exceeded. Flexible multibody dynamics is applied to the morphing forward wing section to assess both morphed shape as well as the stresses. Figure 7 shows the stress distribution caused by tension or compression and bending in the wing skin during the motion of the mechanism from the low-lift (high-speed, undeformed) configuration to the high-lift (lowspeed, morphed) configuration. The stress values increase during the motion and the highest values are found in morphed configuration close to the constraints of the wing skin (points F and G). The design is limited by the maximum stress of all elements and all time steps. To achieve this with a continuous, differentiable function, we approximate the maximum value via constraint aggregation (also called constraint lumping) with the Kreisselmeier-Steinhauser function F KS [41], The modified Kreisselmeier-Steinhauser function is used here to avoid numerical overflow and is defined by [42], where n f is the number of function values and ρ is the aggregation parameter (also referred to as penalty), for which a value of 200 is used here. Specifically in this case, we approximate the maximum stress by with the stress σ , the number of time steps n t and the number of elements n e . Table 1 shows a comparison of the primal analysis of the design described above, where the maximum stress is given by 111.586975 MPa. When using a small value for the aggregation parameter, there is a small error between the true maximum stress and the approximated value. Increasing the value of the aggregation parameter, the error reduces and with large aggregation parameters with ρ ≥ 50, the approximated value is exactly the true maximum stress and no relative error can be detected. The undeformed high-speed configuration, the target morphed low-speed configuration and the achieved deformation is seen in Fig. 8a. The error between the target morphed geometry and the deformation field achieved is calculated for each node. For the kth node, the distance between the simulated node position and the target profile in morphed configuration is where x k,sim and y k,sim are the simulated node coordinates, x k,tar and y k,tar are the target node coordinates and n N is the number of nodes. Figure 8b shows the error between the simulated profile and the target profile in morphed configuration. The largest error is 1.185 mm and is located at approximated the front quarter point on the lower side of the wing profile.
We quantify the deviation of the achieved deformation from the target shape via rootmean-square error, which is given by 0.6348 mm for the given design. The deviations from the target wing profile with a maximum value of 1.124 mm and a root-mean-square error of 0.6225 mm are deemed to be acceptable accuracy, considering the chord length of 550 mm, the wing height of 67.5 mm and the chord length of the forward wing section of 137.5 mm. As the shape of the wing profile is crucial for the aerodynamics and therefore the performance of sailplanes, it is desired to minimise the deviations. The design variables for this investigation include the thickness of the external wing skin t w , the torque value applied to actuate the mechanism M t , the coordinates x A and y A of the torque application point A and the Young's modulus E. The method shown in § 3 is applied to compute the sensitivity values and the chain rule is applied to all operations. Figure 9 shows the design sensitivities of the stress during the motion of the flexible multibody system. 1 During the motion, the highest design sensitivities for all design variables are found at the end of the simulation in morphed configuration. The location of the highest values on the wing profile for the sensitivities with respect to the x-coordinate x A of the torque application point are found on the tip of the wing profile and at the revolute joint in point C, where the actuation mechanism is coupled to the external wing skin, while for all other design variables the location of the highest sensitivity values is close to the constraints of the external wing skin (points F and G). While the maximum stress values in the primal analysis are found on the lower constraint (point G), the highest stress sensitivity values in the sensitivity analysis are found on the upper constraint (point F). Comparing the Figs. 7 and 9 shows that a reduction of the stress σ can be obtained by increasing the thickness of the external wing skin t w , the x-coordinate x A of the torque application point and the Young's modulus E or by reducing the actuation torque M t and the y-coordinate y A of the torque application point.  Figure 9 shows the design sensitivities of the error between the simulated profile and the target profile in morphed configuration. Thereby stands out that especially in the nose of the forward wing section, the sign of the design sensitivities on the lower side is opposite to the sign of the design sensitivities on the upper side. The comparison of the error values on Fig. 8b with the error sensitivities on Fig. 10 shows how the design variables should be modified for reducing the error. The maximum values of the geometric deviation between the simulated profile and the target profile are found on the lower side of the forward wing section and a reduction of the error is therefore obtained by increasing the thickness of the external wing skin t w , the x-coordinate x A of the torque application point and the Young's modulus E or by reducing the actuation torque M t and the y-coordinate y A of the torque application point. Table 2 gives the results of the sensitivity analysis of the maximum stress and the rootmean-square error. The differentiation of the Kreisselmeier-Steinhauser function is defined by where the function of interest f is the stress σ . The direct differentiation of the error is ∇ k = x k,sim − x k,tar ∇x k,sim + y k,sim − y k,tar ∇y k,sim and the direct differentiation of the root-mean-square error is It can be seen that Eq. (53) is infinity when the error at one point is exactly zero. Where this can have adverse effects, e.g. gradient-based optimization, mean square error, which does not have this issue, should be used.
The results in Table 2 show consistency with the results in Fig. 9 and Fig. 10. For the maximum stress, the highest value is found on the constraint at the lower end of the forward wing section (point G). The design sensitivities of Fig. 9 in this location correspond exactly to the approximated Kreisselmeier-Steinhauser sensitivities reported in Table 2. In the computation of the root-mean-square error, all nodes are considered by the function and therefore it is not possible to perform a similar comparison. However, the order of magnitude of the values from Fig. 10 and Table 2 is the same.
The introduced semi-analytical method for the sensitivity analysis of flexible multibody dynamics is validated with a numerical sensitivity analysis with forward differences using a relative perturbation equal to 10 −6 × the value of the design variable itself. The values coincide with the relative difference shown in Table 2. Table 3 shows a comparison of the computational effort of both methods. 2 This comparison highlights the high computational efficiency of the introduced semi-analytical sensitivity  With numerical sensitivity analysis, the computational effort of the numerical sensitivity analysis is given by the sum of n x + 1 primal evaluations with a mean computation time of 62 s. With the semi-analytical sensitivity analysis, the primal analysis and the sensitivity analysis result in computational effort that is less than twice (1.403×) the computational effort of one primal evaluation. The high efficiency of the semi-analytical method compared to the numerical method is already seen with the five design variables used here. This difference becomes even more evident for increasing numbers of design variables. The main reasons for the speed-up are the decoupling of the sensitivity analysis from the FE model explained in § 3.1 and the simplification shown in § 3.3 that allows to reuse the two-dimensional Jacobian of the primal analysis in the sensitivity analysis.

Conclusion
The present work introduces an efficient method for the computation of design sensitivities of flexible multibody dynamics using floating frame of reference formulation with a semi-analytical direct differentiation approach with generalized-α time integration and Newton-Raphson iterations for a morphing wing concept. This concept for a morphing forward wing section in combination with a trailing flap for low-lift (high-speed) and highlift (low-speed) configurations increases the performance in the wide speed range of highperformance sailplanes.
The morphing forward wing section is modeled with flexible elements for the external wing skin and rigid bodies for the actuation mechanism. As system responses, the error between the simulated forward wing profile and the target profile in morphed configuration and the stress values in the wing skin are investigated. The maximum stress is considered as material limit and knockdown factor, though the consideration of the fatigue (including sensitivities) would be of great practical use in future work. Special interest is devoted to the design sensitivities of the system responses.
A semi-analytical sensitivity approach is implemented for its computational efficiency and accuracy. A direct differentiation is chosen to avoid the reverse time integration needed for an adjoint variable methodology. That said, the adjoint variable method exhibits great computational effort advantages with high numbers of design variables, e.g. in topology optimization, and further investigation and comparison in the regard is warranted.
The high efficiency of the method is guaranteed by the reuse of the primal Jacobian for the sensitivity analysis and the decoupling of the sensitivity analysis from the FE model by the analytical differentiation of the system parameters w.r.t. the design variables to a further layer with an invariant-based approach. The design sensitivities calculated include those with respect to geometric properties, material properties and the position and value of the actuator load.
The present study is an integral part of a design optimization framework to optimally design a morphing wing and its actuation mechanism, which is being developed. In this context, the shown sensitivity analysis is the basis for efficient design optimization and uncertainty analysis specifically of the morphing wing and its actuation mechanism and flexible multibody dynamics in general.