River Rapids Conference Proceedings of the Society for Experimental Mechanics Series Topics in Modal Analysis, Volume 10 Michael Mains Proceedings of the 33rd IMAC, A Conference and Exposition on Structural Dynamics, 2015 River Publishers
Conference Proceedings of the Society for Experimental Mechanics Series Series Editor TomProulx Society for Experimental Mechanics, Inc. Bethel, CT, USA
River Publishers Michael Mains Editor Topics in Modal Analysis, Volume 10 Proceedings of the 33rd IMAC, A Conference and Exposition on Structural Dynamics, 2015
Published, sold and distributed by: River Publishers Broagervej 10 9260 Gistrup Denmark www.riverpublishers.com ISBN 978-87-7004-914-6 (eBook) Conference Proceedings of the Society for Experimental Mechanics An imprint of River Publishers © The Society for Experimental Mechanics, Inc. 2015 This work is subject to copyright. All rights are solely and exclusively licensed by the Publisher, whether the whole or part of the material is concerned, specifically the rights of translation, reprinting, reuse of illustrations, recitation, broadcasting, or reproduction in any other physical way, and transmission or information storage and retrieval, electronic adaptation, computer software, or by similar or dissimilar methodology now known or hereafter developed. The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective laws and regulations and therefore free for general use. The publisher, the authors, and the editors are safe to assume that the advice and information in this book are believed to be true and accurate at the date of publication. Neither the publisher nor the authors or the editors give a warranty, expressed or implied, with respect to the material contained herein or for any errors or omissions that may have been made. The publisher remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Preface v Topics in Modal Analysis represents one of ten volumes of technical papers presented at the 33rd IMAC, A Conference and Exposition on Structural Dynamics, 2015, organized by the Society for Experimental Mechanics, and held in Orlando, Florida, February 2–5, 2015. The full proceedings also include volumes on Nonlinear Dynamics; Dynamics of Civil Structures; Model Validation and Uncertainty Quantification; Dynamics of Coupled Structures; Sensors and Instrumentation; Special Topics in Structural Dynamics; Structural Health Monitoring & Damage Detection; Experimental Techniques, Rotating Machinery & Acoustics; and Shock & Vibration, Aircraft/Aerospace, Energy Harvesting. Each collection presents early findings from experimental and computational investigations on an important area within Structural Dynamics. Topics in Modal Analysis I represents papers on enabling technologies for Modal Analysis measurements and applications of Modal Analysis in specific application areas. Topics in this volume include: Experimental Techniques Processing Modal Data Rotating Machinery The organizers would like to thank the authors, presenters, session organizers, and session chairs for their participation in this track. Cincinnati, Ohio, USA Michael Mains Bruel & Kjaer North America
Contents 1 An Efficient Treatment of Parameter Identification in the Context of Multibody System Dynamics Using the Adjoint Method.................................................................................. 1 Karim Sherif, Karin Nachbagauer, Stefan Oberpeilsteiner, and Wolfgang Steiner 2 Calibration and Validation of a Car Subframe Finite Element Model Using Frequency Responses .......... 9 T.J.S. Abrahamsson, F. Bartholdsson, M. Hallqvist, K.H.A. Olsson, M. Olsson, and Å. Sällström 3 Subspace Identification of a 5 MW Reference Wind Turbine ...................................................... 23 Jennifer M. Rinker and Henri P. Gavin 4 Dynamic Analysis of Complex Mechanical Structures Using a Combination of Numerical and Experimental Techniques.......................................................................................... 31 Dimitrios Giagopoulos and Sotirios Natsiavas 5 Predicting Approximate Governing Formula from Experimental Observations ...................................................................................... 41 Sagar Agarwal and R.P. Shimpi 6 Experimental Characterization and Simulation of Vibration Environmental Test .............................. 45 Washington J. DeLima and Melanie N. Ambrose 7 System Identification of an MDOF Experimental Structure with a View Towards Validation and Verification.......................................................................................................... 57 K. Worden, O.D. Tiboaca, I. Antoniadou, and R.J. Barthorpe 8 Some Non-conventional Boundary Conditions (From Marshmallows to Plungers: Who Would Have Guessed) ........................................................................................................... 67 Tina Dardeno, Patrick Logan, and Peter Avitabile 9 Robustness of Disc Brake Systems Regarding Squeal ............................................................... 79 Philippe Stegmann, Sebastian Kruse, and Klaus Augsburg 10 Measurement of Vibration Resulting from Non-contact Ultrasound Radiation Force .......................... 87 Thomas M. Huber, Spencer M. Batalden, and William J. Doebler 11 The Use of Fiber Bragg Grating Sensors for Strain Modal Analysis .............................................. 93 Fábio Luis Marques dos Santos, Bart Peeters, Ludo Gielen, Wim Desmet, and Luiz Carlos Sandoval Góes 12 Using Mode Shapes for Real Time ODS Animation................................................................. 103 Brian Schwarz, Shawn Richardson, and Mark Richardson 13 Removing Unwanted Noise from Operational Modal Analysis Data .............................................. 115 William K. Bonness and David M. Jenkins 14 Adaptive-Like Vibration Control in a Three-Story Building-Like Structure with a PZT Stack Actuator.... 123 G. Silva-Navarro, F. Beltrán-Carbajal, and L.G. Trujillo-Franco vii
viii Contents 15 A Fast Maximum Likelihood-Based Estimation of a Modal Model................................................ 133 Mahmoud El-kafafy, Giampiero Accardo, Bart Peeters, Karl Janssens, Tim De Troyer, and Patrick Guillaume 16 An Improved Implementation of the Orthogonal Polynomial Modal Parameter Estimation Algorithm Using the Orthogonal Complement....................................................................... 157 William Fladung and Håvard Vold 17 An Orthogonal View of the Polyreference Least-Squares Complex Frequency Modal Parameter Estimation Algorithm................................................................................................... 171 William Fladung and Håvard Vold 18 Operational Modal Parameter Estimation from Short-Time Data Series ........................................ 183 R. Arora, A. Phillips, and R. Allemang 19 Order Based Modal Analysis Versus Standard Techniques to Extract Modal Parameters of Operational Wind Turbine Gearboxes ............................................................................. 199 S. Manzato, E. Di Lorenzo, A. Medici, F. Vanhollebeke, B. Peeters, and W. Desmet 20 Modal Identification Results of Quasi-statically Tested RC Frames at Different Damage Levels.............. 215 Ozgur Ozcelik, I. Serkan Misir, Carmen Amaddeo, Umut Yucel, and Erkan Durmazgezer 21 An Innovative Tool for Simulation and Control of a Flying-Cutting Machine ................................... 227 S. Cinquemani and H. Giberti 22 Modal Analysis and Testing of Honeycomb Sandwich Composites................................................ 237 Adarsh Kumar and Ramesh S. Sharma 23 Towards an Automatic Modal Parameter Estimation Framework: Mode Clustering........................... 243 Majid Khorsand Vakilzadeh, Vahid Yaghoubi, Anders T. Johansson, and Thomas J.S. Abrahamsson
Chapter 1 An Efficient Treatment of Parameter Identification in the Context of Multibody System Dynamics Using the Adjoint Method Karim Sherif, Karin Nachbagauer, Stefan Oberpeilsteiner, and Wolfgang Steiner Abstract In multibody system dynamics, a wide range of parameters can occur where some of them may not be known a priori. Therefore, this work presents an efficient adjoint method for parameter identification that can be utilized in multibody simulation software. Compared to standard system sensitivity based approaches the adjoint method has the major advantage of being independent on the number of parameters to identify. Especially when dealing with large and probably flexible multibody systems this characteristic is crucial. Formulating parameter identification as an automatable procedure, of course, leads to a complicated structure of the involved matrices and equations. However, during a forward simulation of the system many of the matrices needed for solving the so called “adjoint system equations” are already evaluated. Adopting the functionality of the forward solver for the adjoint system solver therefore results in little additional effort. In order to illustrate the performance of the adjoint method two examples are presented. A planar example shows the possibility of identifying non-linear control parameters and a three-dimensional example is presented for identifying time-invariant inertia parameters. Keywords Adjoint method • Parameter identification • Optimal control • Inverse dynamics • Multibody systems 1.1 Introduction The adjoint method is probably the most efficient way to solve a variety of optimization problems in engineering sciences. Much attention to this approach has been paid recently in the context of continuous systems for sensitivity analysis (see, e.g., [1–4]). The class of dynamic programming methods for the computation of gradients in an optimization problem includes the adjoint method, which has a long history in optimal control theory [5]. In the explanations of Giles and Pierce [6], the adjoint method is seen as a special case of linear duality, in which the dual problem has to solve only a single linear system. Obviously, huge benefits can be achieved by solving the dual formulation. The adjoint method utilizes this powerful aspect of duality to dramatically improve the efficiency of the computation. Previous work on the adjoint method in multibody dynamics can be found in the work of Bottasso et al. [7], where the solution of inverse dynamics and trajectory optimization problems for multibody systems is reflected by an indirect approach combining optimal control theory with control and adjoint equations and transversality conditions. The design of the indirect method for solving optimal control problems for multibody systems presented by Bertolazzi et al. [8] seems to be familiar with the idea of the adjoint method. The work by Schaffer [9] presents a numerical algorithm, the piecewise adjoint method, which formulates the coordinate partitioning underlying ordinary differential equations as a boundary value problem, which is solved by multiple shooting methods. The group around Petzold, Cao, Li and Serban [10–12] describe forward and adjoint methods for sensitivity analysis for differential-algebraic equations and partial differential equations, and state that the results of sensitivity analysis have wide-ranging applications in science and engineering, including model development, optimization, parameter estimation, model simplification, data assimilation, optimal control, uncertainty analysis and experimental design [10]. In the work of Eberhard [4], the adjoint method is used for sensitivity analysis in multibody systems interpreted as a continuous, hybrid form of automatic differentiation. Despite of the great potential of the adjoint method, in multibody dynamics the adjoint method is rarely applied, since the structure of the equations of motion is usually extremely complicated, in particular if flexible bodies are included, and the effort to obtain the set of adjoint equations seems tremendously high. Hence, dealing with the adjoint method is obviously unattractive for most developers of multibody simulation software, as far as they are familiar with it at all. The main goal K. Sherif ( ) • K. Nachbagauer • S. Oberpeilsteiner • W. Steiner University of Applied Sciences Upper Austria, Wels Campus, Stelzhamerstr. 23, 4600 Wels, Austria e-mail: karim.sherif@fh-wels.at; karin.nachbagauer@fh-wels.at; stefan.oberpeilsteiner@fh-wels.at; wolfgang.steiner@fh-wels.at © The Society for Experimental Mechanics, Inc. 2015 M. Mains (ed.), Topics in Modal Analysis, Volume 10, Conference Proceedings of the Society for Experimental Mechanics Series, DOI 10.1007/978-3-319-15251-6_1 1
2 K. Sherif et al. of the present paper, which is a shortened form of [13], is to show how the adjoint method can be embedded efficiently to a multibody system described by a system of differential-algebraic equations of index 3 for optimal control problems or parameter identification applications. The present paper shows the potential of the adjoint method for solving classical optimization problems in multibody dynamics and presents applications for optimal control problems and a parameter identification. 1.2 The Adjoint Method in Multibody Dynamics In this section we discuss the application of the adjoint gradient computation to multibody dynamics. The adjoint equations for the equations of motion of a multibody system will be derived and a flowchart for embedding the adjoint method in multibody dynamics is illustrated. 1.2.1 Adjoint Equations In the augmented formulation the dynamics of a multibody system is described by a set of differential-algebraic equations in the following form: Pq Dv (1.1) MPv Df.q; v; u/ CT q (1.2) C.q/ D0 (1.3) where qand v denote the generalized coordinates and velocities. Moreover, umay either describe a vector of time-invariant parameters or a vector of time-dependent control signals actuating the system. For simplicity we assume that uappears only in the vector of generalized forces and gyroscopic forces f, which is e.g. the case if u is a stiffness or damping parameter or an actuating force. The matrix Mrepresents the symmetric mass matrix of the system. The algebraic constraint equation C.q/ D0 influences the equations of motion via the Jacobian Cq and the vector of Lagrange multipliers . The key idea of the adjoint method (see [4, 10–12] for example) is to determine the parameter/control u such that a functional of the form J.u/ DZ T 0 h.q; v; u/dt CS.q; v/ˇ ˇ ˇtDT (1.4) is minimized, in which the function h.q; v; u/ and the end point termS.q; v/jT have to be defined appropriately. The end point termS.q; v/jT is often also denoted as scrap function. Numerous methods are available to compute the argument for which a function or a functional attains a minimum. We just refer to the method of the steepest descent, the conjugate gradient method, the Gauss-Newton method or quasi Newton methods like the BFGS algorithm when estimating the Hessian or the DFP algorithm when estimating the inverse of the Hessian. Some authors embed these methods in a homotopy continuation to obtain a global minimum [14]. In any cases the gradient of J.u/ must be determined. For this purpose several strategies can be pursued again. On the one hand, if uis a vector of Nu parameters, the sensitivity equations for xu D@x=@uare usually considered (see [15, 16] or [14] for example). The computational effort for this approach is equal to solvingNu linear sets of equations with the same dimension as Eqs. (1.1)–(1.3). On the other hand, if urepresents a vector of time-dependent control signals, these signals are often discretized in order to transform the problem into a finite dimensional one. The adjoint method is a powerful alternative to compute the gradient of J.u/ in both cases. First, we notice that J.u/ does not change if we incorporate Eqs. (1.1)–(1.3) to the integrand J.u/ DZ T 0 h.q; v; u; t/ Cp T .Pq v/ CwT .MPv f CCT q / C TC dt CS.q; v/ˇ ˇ ˇtDT (1.5) no matter how the adjoint variables p.t/, w.t/ and .t/ are chosen.
1 An Efficient Treatment of Parameter Identification in the Context of Multibody System Dynamics Using the Adjoint Method 3 Let us now consider a forward solution q.t/, v.t/ and .t/ of the system Eqs. (1.1)–(1.3) for a set of parameters or control variables u. A variation of uwill result in variations of the variables q, v and , and, moreover, in a variation of the functional J. Thus, the variation of the functional J can be written in the form ıJ DZ T 0 (hqıqChvıvChuıuCp T .ıPq ıv/C CwT h .MPv/q ıqCMıPv fqıq fvıv fuıuC CT q q ıqCCT qı iC TCqıq) dt CSqıqˇ ˇ ˇtDT CSvıvˇ ˇ ˇtDT (1.6) where hq, hv and hu denote the partial derivatives of the function hwith respect to the components of q, v and u. In order to get rid of the variations of Pqand Pv we integrate by parts Z T 0 pT ıPqdt D Z T 0 P pT ıqdt Cp T ıqˇ ˇ ˇtDT Z T 0 wTMıPvdt D Z T 0 d dt wTM ıvdt CwTMıvˇ ˇ ˇtDT (1.7) where the fact has been used that ıq.0/ D0 and ıv.0/ D0 since the initial conditions for the variables q and v are given. Using Eq. (1.7) the variation of the functional J given by Eq. (1.6) can be rewritten as ıJ DZ T 0 ( hhq Pp T CwT .MPv/q fq C CT q q C TCqiıq C hv p T wTf v d dt wTM ıvC wTCT q ı C hu wTf u ıu) dt C Sq Cp T ıqˇ ˇ ˇtDT C Sv CwTM ıvˇ ˇ ˇtDT (1.8) The key question now is how the variation of J is related to the variation of u. In order to obtain this relationship we choose our adjoint variables p.t/, w.t/ and .t/ such that the terms involving ıq, ıv and ı vanish from the integrand. Thus, we end up with the following system of adjoint equations: dp dt D h T q CAwCCT q (1.9) d dt .Mw/ Dh T v p Bw (1.10) 0 DCqw (1.11) 0 DS T q Cp.T/ (1.12) 0 DS T v CM.q.T//w.T/ (1.13) with the following abbreviations AD.MPv/ T q f T q C.CT q / T q (1.14) BDfT v (1.15)
4 K. Sherif et al. Using this set of adjoint equations, Eq. (1.8) reduces to ıJ DZ T 0 hu wTf u ıudt (1.16) which directly relates the independent variation ıu to the variation of the objective function. In order to find the minimum of J, we walk along its negative gradient for a finite distance and update the gradient form time to time until we reach a point where the gradient gets zero. If u is a control signal, it is easy to show that the largest infinitesimal increase of ıJ is obtained, if ıu.t/ D hu wTf u T .control optimization/ (1.17) where is an infinitesimal positive factor. Hence Eq. (1.17) is the gradient formula for the cost functional we are looking for. Moreover, if uis a vector of time-invariant parameters, Eq. (1.16) can be written as ıJ D Z T 0 hu wTf u dt ıu DrJ T ıu (1.18) where the gradient of the functional J is defined as rJ DR T 0 hu wTf u T dt. Thus, the variation of ucan be written for the case of time-invariant parameters as ıu D rJ D Z T 0 hu wTf u T dt .parameter optimization/ (1.19) If is sufficiently small, the updated control/parameter uCıuwill always reduce J. The adjoint equations (1.9)–(1.13) have to be solved backwards in the physical time. Therefore it is advantageous to introduce a new time coordinate by the following transformation DT t; 2Œ0;T ; d dt D d d d dt D d d For the numerical solution of the adjoint equations we propose a backward differentiation scheme which approximates the derivative of a function F. / at a time instant n by using the function values at n; n 1; : : : ; n k. The backward differentiation formula (BDF) reads dF d ˇ ˇ ˇ ˇ n 1 k Xi D0 ˛i F. n i / (1.20) The coefficients ˛i result from differentiating an interpolation polynomial throughF. n/; : : : ;F. n k/ and are chosen as the standard coefficients presented, e.g., in [17, p. 349]. A closer look at the adjoint equations (1.9)–(1.13) shows that the boundary condition (1.13) of the variable wis in general incompatible with the adjoint constraint equation (1.11). Only when Sv D0, i. e. when the scrap function does not depend on v, all equations are satisfied by setting p.T/ D ST q and w.T/ D0. For the case that Sv ¤0 a consistent boundary condition approach has to be applied, see [13] for more details on this issue. 1.2.2 Flowchart of the Adjoint Method Embedded in Multibody Systems At this point it should be mentioned that only two systems of DAEs must be integrated for computing the direction of the gradient, that is Eqs. (1.1)–(1.3) and (1.9)–(1.13). Subsequently, we give a flowchart of those steps, which are needed for the computation of the gradient of the objective function J with the adjoint method. Let u denote a given vector of controls or parameters of a multibody system. The flowchart describes how one can successively decrease the cost function J until a user-defined limit is achieved by varying u.
1 An Efficient Treatment of Parameter Identification in the Context of Multibody System Dynamics Using the Adjoint Method 5 1. Solve the equations of motion Eqs. (1.1)–(1.3) forward in time t 2 Œ0;T yielding q.t/, v.t/ and .t/. This may be done e.g. by choosing the Hilbert-Hughes-Taylor (HHT) integration scheme, as proposed in [18] and its application for a differential algebraic system given in an index three formulation in [19]. 2. Compute the objective function J by inserting q.t/, v.t/ and u.t/ into Eq. (1.4). Note, that the integration must be done numerically. If J is smaller than a user defined minimum then stop. 3. Along the forward simulation of the equations of motion compute the mass matrixM, the constraint JacobianCq and the Jacobian matrices A, Bfrom Eqs. (1.14) and (1.15). 4. Determine the consistent boundary conditions at D0 for the adjoint variables w, pas suggested in [13]. 5. Solve the adjoint equations (1.9)–(1.13) for p. /, w. / and . /, where DT t. 6. Compute the adjoint variables as functions of the original time by settingp.t/ Dp. DT t/ andw.t/ Dw. DT t/. Moreover, determine hu and fu along the forward simulation. 7. Compute u DuCıu. Use Eq. (1.17) for the computation of ıu if u is a control signal. Otherwise, if u is a set of timeinvariant parameters, use Eq. (1.19) for the update. A sufficiently small number >0 has to be chosen for that purpose. Go to step 1. For the sake of overview, it should be mentioned at this point that the Jacobian matrices.MPv/q, fq, fv and.CT q /q whichare needed in the adjoint equations may be required already for the simulation of the multibody system if an implicit integration scheme such as, e.g., the HHT-algorithm [18, 19] is applied. 1.3 Numerical Examples Two numerical examples will be presented to illustrate the application of the adjoint method in typical multibody systems. As a first example, a planar overhead crane is considered as an example of an underactuated mechanical system which follows a given trajectory and the optimization process identifies the control force and the control torque in a specific time domain. The second example incorporates a single rigid body which is parametrized by the four redundant Euler parameters. A point of the body is excited in order to follow a specific trajectory and the optimization process identifies all entries of the inertia tensor. A time history of the identified forces and moments as well as convergence analyses for the cost functionals and the inertia parameters are presented. 1.3.1 Planar Overhead Crane The overhead crane presented in [20, 21] is a classical example of an underactuated mechanical system consisting of the cart A, the hoisting drumB, thecableC andamass Dmounted at the end of the cable, see Fig. 1.1a. The generalized coordinates of this two-dimensional problem are chosen as the position of the cart in x-direction xc, the length of the cable l and the position of the mass inx- andy-direction, xm andym, respectively, see Fig. 1.1b. The goal of this example is to compute the drive signals F and M, depicted in Fig. 1.1b, such that the point mass Dfollows a given trajectory defined by a linear path from a given starting point .xm0;ym0/ D.0=4/ to a given end point .xmf ;ymf / D.5=1/ within 3 s, see Fig. 1.2. The control force F and control torque Mare identified within the time periodt 2Œ0;3 . For the computation we set the mass of the cart and the hoisting drumACBto 10kg, the mass of Dis set to 100kg, the moment of inertia of the hoisting drum is defined as 0:1kgm2 and the radius of the hoisting drum is given by0:1m. As starting point of our drive signals we choose the static solution, i.e. M0 D100 9:81 0:1 D98:1Nmand F0 D0N. The results in Fig. 1.3 are computed with a constant step size and show the initial settings for the first iteration and the identified control force F and control torque Mafter 300 iterations. The optimization process reduces the costs to a factor 10 7 within 300 iterations, see Fig. 1.4. 1.3.2 Single Rigid Body Parametrized with Euler Parameters The second example fall in the category of parameter identification, i.e. uis time-invariant. The goal of the present example is to identify the components of the inertia tensor (I11;I22, I33, I12;I13 and I23) of a single three-dimensional rigid body. The position of the rigid body is described by the coordinates of the center of mass x D.x;y; z/ and for the description of
6 K. Sherif et al. Fig. 1.1 Geometry description of the planar, rigidly modeled overhead crane (a) Planar overhead crane (b) Parameters of planar overhead crane A B C D xc M F xm ym l a b Fig. 1.2 The point mass has to follow a linear trajectory from a specified starting point to a fixed end point Fig. 1.3 Time history of identified force F and torque M after 300 iterations and initial settings for the first iteration for Sect. 1.3.1. The initial input for the force is set toF0 D0Nfor the first iteration. The initial input for the torque is defined as the static torque M0 D98:1Nm −600 600 Force [N] time [s] F0 F Torque [Nm] time [s] M0 M 0 0 0 0.5 0.5 1 1 1.5 1.5 2 2 2.5 2.5 3 3 -400 -200 200 400 140 120 100 80 the orientation of the body Euler parameters D. 0; 1; 2; 3/ are used. Figure 1.5 shows the arbitrarily shaped rigid body for which the inertia parameters describing the inertia tensor are not known in advance. The point S of the body follows a prescribed motion realized by a constraint and the velocity of the point P is measured (in global coordinates) in order to identify the entries of the inertia tensor. The geometry used for generating the inertia data is a cuboid with one diagonal congruent with the z-axis. The excitation for the identification process is applied in the point S situated at the origin at time t D0. Velocity measurements are taken in the point P situated in the opposite direction toS along thez-axis. In order to start the optimization, initial values for the parameters to identify have to be defined. Table 1.1 shows the correct, the initial values for the first iteration and the final identified values for the moments of inertia after 126 iterations, for which the values are assumed to be converged, see Fig. 1.6. The convergence analysis of the cost functional shows that the optimization process
1 An Efficient Treatment of Parameter Identification in the Context of Multibody System Dynamics Using the Adjoint Method 7 Fig. 1.4 The convergence analysis of the cost functional shows that the optimization process reduces the costs to a factor of 10 7 within300 iterations iteration costs 100 10−1 10−2 10−3 10−4 10−5 10−6 10−7 0 50 100 150 200 250 300 S P x Fig. 1.5 A single rigid body is studied for which the moments of inertia parameters describing the inertia tensor are not known. Point S follows a specified motion and the velocity of point P is measured in order to identify the entries of the inertia tensor Table 1.1 Moments of inertia: the correct values, the initial values for the first iteration and the final identified values after 126 iterations, for which the values are assumed to be converged Parameter Correct value Initial value Final identified value Unit I11 132.285 1.0 133.441 kgm2 I22 162.054 1.0 162.073 kgm2 I33 92.809 1.0 96.815 kgm2 I12 22.382 0.0 22.198 kgm2 I13 42.342 0.0 44.593 kgm2 I23 2.689 0.0 3.011 kgm2 Fig. 1.6 The convergence analysis of the cost functional shows that the optimization process reduces the costs already tremendously within the first 100 iterations. It has to be mentioned that only the costs of every second iteration are depicted here iteration costs 10−3 10−6 10−9 0 20 40 60 80 100 120 140 reduces the costs already within the first 100 iterations to a factor 10 8, see Fig. 1.6. It has to be mentioned that in case of identifying only the diagonal entries of the inertia tensor, I11;I22 and I33, the results show higher accordance to the correct values as in case of identifying the whole inertia tensor also incorporating the deviation moments of inertia.
8 K. Sherif et al. 1.4 Conclusions The adjoint method has been embedded to multibody system dynamics. The presented method can be used for inverse dynamic problems as well as for parameter identification in multibody systems. The computation of the direction of the gradient of the cost functional is based on solving two DAEs. On the one hand, the equations of motion of the multibody system have to be integrated forward in time, and, on the other hand, the adjoint equations have to be integrated backwards. Depending on the chosen forward time integration scheme, almost all necessary matrices for the backward time integration can be reused from the forward time integration and, therefore, a time- and memory-efficient simulation tool for inverse dynamics in the field of multibody dynamics can be constructed. Therefore, the adjoint method shows an efficient way to incorporate inverse dynamics to flexible multibody system applications arising from modern engineering problems. Acknowledgements This research has been funded by the European Regional Development Fund and the government of the Upper Austria via a Regio 13 project. References 1. Haug E, Ehle P (1982) Second-order design sensitivity analysis of mechanical system dynamics. Int J Numer Methods Eng 18(11):1699–1717 2. Haug E, Wehage R, Mani N (1984) Design sensitivity analysis of large-scaled constrained dynamic mechanical systems. Trans ASME 106: 156–162 3. Bestle D, Eberhard P (1992) Analyzing and optimizing multibody systems. Mech Struct Mach 20:67–92 4. Eberhard P (1996) Adjoint variable method for sensitivity analysis of multibody systems interpreted as a continuous, hybrid form of automatic differentiation. In: Proceedings of the 2nd international workshop on computational differentiation, Santa Fe. SIAM, Philadelphia, pp 319–328 5. Lions J (1971) Optimal control of systems governed by partial differential equations. Springer, New York 6. Giles M, Pierce N (2000) An introduction to the adjoint approach to design. Flow Turbul Combust 65:393–415 7. Bottasso C, Croce A, Ghezzi L, Faure P (2004) On the solution of inverse dynamics and trajectory optimization problems for multibody systems. Multibody Sys Dyn 11:1–22 8. Bertolazzi E, Biral F, Lio MD (2005) Symbolic-numeric indirect method for solving optimal control problems for large multibody systems. Multibody Sys Dyn 13:233–252 9. Schaffer A (2005) On the adjoint formulation of design sensitivity analysis of multibody dynamics. Dissertation, University of Iowa. http://ir. uiowa.edu/etd/93 10. Petzold L, Li S, Cao Y, Serban R (2006) Sensitivity analysis for differential-algebraic equations and partial differential equations. Comput Chem Eng 30:1553–1559 11. Cao Y, Li S, Petzold L (2002) Adjoint sensitivity analysis for differential-algebraic equations: algorithms and software. J Comput Appl Math 149:171–191 12. Cao Y, Li S, Petzold L, Serban R (2003) Adjoint sensitivity analysis for differential-algebraic equations: the adjoint DAE system and its numerical solution. SIAM J Sci Comput 24(3):1076–1089 13. Nachbagauer K, Oberpeilsteiner S, Sherif K, Steiner W (2014) The use of the adjoint method for solving typical optimization problems in multibody dynamics. J Comput Nonlinear Dyn. doi: 10.1115/1.4028417 14. Vyasarayani CP, Uchida T, McPhee J (2011) Parameter identification in multibody systems using lie series solutions and symbolic computation. J Comput Nonlinear Dyn 6(4):041011. doi: 10.1115/1.4003686 15. Ding JY, Pan ZK, Chen LQ (2012) Parameter identification of multibody systems based on second order sensitivity analysis. Int J Non-Linear Mech 47:1105–1110 16. Özyurt DB, Barton PI (2005) Cheap second order directional derivatives of stiff ODE embedded functionals. SIAM J Sci Comput 26(2): 1725–1743 17. Süli E, Mayers D (2003) An introduction to numerical analysis. Cambridge University Press, Cambridge 18. Hilbert H, Hughes T, Taylor R (1977) Improved numerical dissipation for time integration algorithms in structural dynamics. Earthq Eng Struct Dyn 5:283–292 19. Negrut D, Rampalli R, Ottarsson G, Sajdak A (2005) On the use of the HHT method in the context of index 3 differential algebraic equations of multibody dynamics. In: Goicolea JM, Cuadrado J, Garcia Orden JC (eds) Proceedings of the ECCOMAS conference on advances in computational multibody dynamics, Madrid 20. Blajer W, Kołodziejczyk K (2004) A geometric approach to solving problems of control constraints: theory and a DAE framework. Multibody Sys Dyn 11(4):343–364 21. Betsch P, Uhlar S, Quasem M (2009) Numerical integration of mechanical systems with mixed holonomic and control constraints. In: Arczewski K, Fra¸czek J, Wojtyra M (eds) Proceedings of the ECCOMAS thematic conference on multibody dynamics, Warsaw University of Technology, 29th June–2nd July
Chapter 2 Calibration and Validation of a Car Subframe Finite Element Model Using Frequency Responses T.J.S. Abrahamsson, F. Bartholdsson, M. Hallqvist, K.H.A. Olsson, M. Olsson, and Å. Sällström Abstract A finite element model of a car front subframe has been calibrated against test data. Stepped-sine testing has been used to give frequency response function estimates on an ensemble of seemingly identical subframes. Therefore, the deviation between test data and simulation results can be compared in a meaningful way by the outcome of model calibration and cross-validation. Emphasis has been put on the preparation of the test pieces for high fidelity testing and on bettering the chances of getting a calibration outcome that provides insight into the physical processes that govern the subframe dynamics. The front subframe model has more than 200,000 degrees-of-freedom and 17 model calibration parameters. The efficiency of the calibration procedure under these conditions is reported. To achieve efficiency, a calibration with a smooth deviation metric is used together with a damping equalization method that eliminates the need for matching of experimental and analytical eigenmodes. The method is combined with surrogate model frequency response evaluation based on model reduction for increased speed. The Matlab based open-domain software tool FEMcali that employs the Levenberg-Marquardt minimizer with randomized starts has been used for calibration and an unregularized Gauss-Newton minimizer has been used in the cross-validation. Keywords FEM model calibration • FEM model validation • Large model calibration • Surrogate model • Damping equalization 2.1 Introduction For car industry, as well as other advanced industry, it is important to shorten the development time and reduce cost for testing. This has increased the need for attribute predictions that are based on computational models with high credibility. Good modeling guidelines are essential to achieve trust in model predictions that are made prior to the availability of physical test objects. Such modeling guidelines should ideally be created based on the outcome of model correlation and validation studies made on test data from already existing components, subsystems or full systems. Since industry-size models tend to be very large, efficient methods for calibration, correlation and validation are needed in this process of increasing the model credibility. This paper concerns the calibration and validation of a subframe model for a Volvo V40 car. This part of the car was selected since it is represented by a reasonable large computational model and was believed to behave linearly under mild loading conditions. A Nastran model with 207,912 degrees-of-freedom is the target for calibration. The subframe and its physical location in the car is indicated in Fig. 2.1. Many FEM calibration practices are based on modal analysis and modal comparisons. The accuracy of the FEM is determined by comparing modal data from test and analysis. The eigenfrequencies are compared directly, while the corresponding mode shapes are compared with various metrics. The metrics most often employed are based on the orthogonality and cross-orthogonality of the modes with respect to the mass distribution quantified by a FEM mass matrix, see [1]. A significant drawback of these metrics is the fact that there is not an explicit connection between the metric values and a corresponding deviation in a predicted response. Another drawback is that measurement noise, limited spatial resolution in the vibration sensing, and shortcomings that affect the test data quality and processed data often magnify the problems associated with closely spaced modes that produce a spurious coupling sensitivity between the test and FEM mode shapes. This sensitivity makes a modal orthogonality metric problematic for matching test and FEM mode shapes T.J.S. Abrahamsson ( ) Department of Applied Mechanics, Chalmers University of Technology, Hörsalsvägen 7A, S-41296 Göteborg, Sweden e-mail: thomas.abrahamsson@chalmers.se F. Bartholdsson • M. Hallqvist • K.H.A. Olsson • M. Olsson • Å. Sällström Volvo Car Corporation, S-40531 Göteborg, Sweden © The Society for Experimental Mechanics, Inc. 2015 M. Mains (ed.), Topics in Modal Analysis, Volume 10, Conference Proceedings of the Society for Experimental Mechanics Series, DOI 10.1007/978-3-319-15251-6_2 9
10 T.J.S. Abrahamsson et al. Fig. 2.1 Front subframe seen as placed in the car between the front wheels. The subframe is a vital part in the front wheel suspension and connects the suspension link arms to the car body and engine which makes it an important part regarding transmitted vibrations from road and engine in the process of FEM calibration for systems with high modal density. Also, as the vibrational wavelengths approach the scale of structural geometry variations, the mode shapes become very sensitive to modelling errors and uncertainties [2–4]. Alternatively, modal approaches could be avoided by using frequency responses for FEM validation. Many researchers, e.g. Grafe [5], have studied model calibration using frequency domain data. The methods developed can loosely be categorized based on their criterion deviation metric as equation error methods or output error methods, or on their parameterization strategy, as direct matrix entity parameterization or physical property parameterization. In frequency domain equation error methods, such as described in [6–12], sensitivities of the FEM impedance matrix with respect to the design variables must be computed. One observed drawback of equation error approaches is that the comparably small number of responses measured in comparison to the large number of dofs in the FE model can produce non-smooth reduced impedance functions that have sharp peaks at the resonances that occur when the sensor dofs are fixed, see [6]. This can lead to a complicated calibration process in which a numerical calibration minimizer can easily be stuck in a local minimum. Output error methods on the other hand, such as described in References [13–18], attempt to minimize the difference between the measured data and the analytical prediction of that data. Balmès [13] minimized the least-squares or log-least-squares deviation in the frequency response functions. Herendeen et al. [14] presented a general multidisciplinary optimization framework for minimizing the deviation in any type of structural response. Dascotte and Strobbe [15] minimized the error in a frequency response function correlation metric that is based on the modal assurance criterion adapted to the frequency domain. The drawback of this approach is that the deviation metric, like many of the frequency response based deviation metrics in use [18], are not physically well motivated. There is no direct connection between a metric deviation and the deviation in the FEM prediction of responses due to loading. The advantage of the output approach is that no model dof reduction is necessary. However, these techniques require the minimization of an objective function that is non-linear in the parameters that can lead to minimizer convergence problems and a large computational effort. Early developed FEM calibration strategies directly targeted the matrix entities of the FE mass and stiffness matrices by adjusting these to achieve best fit to test data. Recent papers, of more practical use to an engineer in the product development process, deal with model parameters that have clear physical interpretations, such as Youngs’ moduli, material densities and thicknesses of shell-like parts. In a paper by Abrahamsson and Kammer [19] a discrete-frequency response based deviation metric was presented. It was shown that, together with a method for damping equalization that makes the calibration problem regularized, the metric gives good convergence properties for the calibration optimizer. The smoothness of the calibration metric was demonstrated to increase with increased system damping. In the damping equalization procedure the system damping can be arbitrarily set to any fictitious level for both the experimental data and the FEM. By assigning the same amount of modal damping to all system modes, modal correlation analysis for mapping of experimentally found damping to FEM damping can be avoided. By considering the FE coefficient matrices to be linear in the parameters around a linearization point, the reduced order model can be easily formed without expensive evaluation of the FEM equation as the model parameter settings vary as the iteration process proceeds to achieve the calibration minimum.
2 Calibration and Validation of a Car Subframe Finite Element Model Using Frequency Responses 11 In this paper, the calibration and validation problem of a large scale FE model is treated by the method set out in [19]. The measurement procedure, involving tests of five seemingly identical subframes, is described. The outcome of the calibration, validation and cross-validation is reported and observations made during the process are discussed. 2.2 Theory The theory presented in [19] is reiterated in condensed form here for reader convenience. We set out the theory from the FE representation of the linear time-invariant mechanical system as the set of ordinary differential equations MRq CV Pq CKq DQ.t/ (2.1) withM;V;K2<n n being the mass, viscous damping and stiffness matrices respectively. The displacement vector is here q, and the dot notation is used for time differentiation. The vectorially associated loadings to q is Q, for which just but a few entities are normally non-zero. These non-zero entities constitute the stimuli vector u 2<nu from which the full load vector is given by the Boolean transformation matrix Pu as QDPuu. It is assumed that only some mass and stiffness properties are subjected to calibration and therefore the FE problem is parameterized, with physical property parameters collected in data vector P 2 <np, such that the mass and stiffness matrices depend on these, i.e. K DK.P/ and M DM.P/. The physical parameter setting P, with entities that can be of any order, is related to a normalized parameter vector p and some fixed non-zero nominal parameter setting P0 such that P DP0 : .1Cp/. The normalization of the parameters is made to better suit numerical minimizing schemes that often works better with scaled unknowns. Using the dummy equation I Rq–I Rq D0, the first-order differential equation counterpart to Eq. (2.1) can be formed as the state-space equations Px DAxCBu y DCxCDu (2.2) with the state vector xT D˚ qT PqT . It can easily be verified that the coefficient matrices Aand Bof Eq. (2.2) relate to the mass, damping and stiffness matrices such that AD 0 I M 1K M 1V BD 0 M 1Pu (2.3) and the matrices Cand Dare the appropriate matrices forming the system output y 2 <ny from linear combination of the states in x 2 <2n and the stimuli in u. Using the state-space representation (2.2), the system transfer function, relating the input to the output, can be expressed as H.p/ DC.i!I A/ 1BCD (2.4) The transfer functionHestablished from analytical FE data will here be denotedHA. The experimentally determined system transfer function (denoted HX) is often obtained from vibration testing and the use of signal processing and state-space system identification. The model calibration problem is thus about finding that parameter setting p Dp (called the oracle parameter setting) that minimizes the deviation betweenHA(p) andHX, i.e. p Dargmin f HA.p/ HX withf being a calibration metric. 2.2.1 A Frequency Response Calibration Metric A calibration scheme that uses a gradient based minimizer, needs to work with a smooth deviation metric for high likelihood of success. That is to obtain convergence in the search for parameter optimum from multiple start settings of the parameters. A well calibrated model should give high accuracy in simulation of test output quantities, and ideally predictions with high credibility of other output quantities not tested. In a frequency domain context, this often translates to that model which accurately captures the structural resonances and possibly also its anti-resonances. A metric that does not discriminate against deviations at frequencies where the structural response is small is the quadratic functional ı D"H"=N (2.5)
12 T.J.S. Abrahamsson et al. which is called the square deviation (SD). In this the deviation vector "is ".p/ Dlog10 vect H A .p/ :=vect H X (2.6) where the ./ operator denotes the element-by-element division and vect(.) is the vectorizing operation that makes all frequency response function elements of the ny nu transfer function, at a dense set of nf discrete frequencies used for evaluation, into an nynunf 1column vector. Nis the number of elements of that vector. Since finite element model calibration tends to be computationally demanding, calibration criteria that lead to computational efficiency are strongly of the essence. One such efficiency concern targets the sampling strategy for the discrete frequencies that are selected for frequency response function evaluation and the selection is based on the half-band width of the eigenmodes. The half-bandwidth !i of a damped structural resonance at frequency !i is given by !i D 2 i !i with i being the relative modal damping of the i:th mode. One observes that the half-bandwidth increases linearly with increasing resonance frequency. It has been found to be a good frequency sampling strategy to utilize frequency steps that increase linearly with frequency to give a balanced contribution to the deviation metric from various system modes. Such sampling keeps the number of samples over one half-bandwidth constant over the range. That sampling strategy seems reasonable, provided that relative damping of all modes in the range is equal, which rarely happens for experimentally found eigenmodes. However, the damping can be equalized by a procedure that is described next. 2.2.2 Damping Equalization A central issue for FRF based model calibration is that of model damping. Since in general, damping has been found to be very difficult to model using first principles, it is most often assigned a simple representation for modeling convenience. For modal damping modelling, the model’s damping is usually set using the outcome of experimental modal analysis of test data from the structure under investigation. The modal dampings found in experiments are normally used for FE simulation without further attempts to understand their physical background. The nature of the damping mechanisms is normally such that the modal damping varies from mode to mode. That makes a mapping of experimentally obtained modal damping into modal damping of FE modes cumbersome. The difficulty arises since the mapping of modal damping relies on mode shape pairing, meaning that the same amount of modal damping should be assigned to modes that are in some sense similar by their deformation pattern. Such pairing is normally difficult, especially for systems with a high modal density. To overcome the problem of mode pairing, a method of damping equalization has been suggested [19]. The damping equalization is achieved by imposing the same modal damping on all experimentally found system modes by perturbation of a mathematical model of the experimental data found from system identification using raw frequency response function data HX raw. The experimentally found system transfer function HX raw can then be represented by an identified state-space system Px D Ax CBu, y D CxCDusuch that HXraw C.i!I A/ 1BCD (2.7) The experimental state-space system can be brought to diagonal form by a similarity transformation. Using the mode matrix X, pertinent to the eigenvalue problemAXDXƒ, for transformation we have for the diagonal realization that x DXz Pz DAz CX 1Bu y DCXz CDu (2.8) with ADX 1AXDdiag. i / (2.9) forwhich i are the complex-valued system poles as given by the experimental data. For small damping, the relative modal damping i, obtained from these poles are i D <. i /=j=. i /j (2.10) In the process of damping equalization, the real parts of the poles are perturbed such that the damping is made equal for all modes. The modal dampings are then set to a single fixed value 0, i.e. i D 0 8i. The effect of such damping equalization
2 Calibration and Validation of a Car Subframe Finite Element Model Using Frequency Responses 13 is that the oscillatory imaginary part of the poles are preserved and the real damping part is modified such that the perturbed system poles Q i are now Q i D=. i /. 0 Ci/ ; 8=. i />0; Q i D=. i /. 0 Ci/ ; 8=. i /<0 (2.11) and the modified state-space realization is then Pz D QAz CX 1Bu y DCXz CDu (2.12) with QA Ddiag Q i . This in turn gives us a modified transfer function for the experimental model, such that the transfer function used for calibration with damping equalization is HX CX i!I QA 1 X 1BCDE (2.13) At this stage it should be obvious that the application of a system identification procedure on the raw test data HX raw has led us to a mathematical model which we can evaluate for any frequency. We are also able to make fictitious modifications of the system under test. A particularly useful such modification is that we can adjust the system damping level, leaving stiffness and inertia properties intact, such that all system modal damping are set equal. The model calibration of the FE model can then be made towards this fictitious experimental model for calibration of parameters that relate to mass and stiffness only. For the FE based system representation, the modal damping allows for a simple representation. For a system with given mass and stiffness matrices Mand Kwe have the viscous damping matrix Vtobe [20] V DMXdiag.1=mi /diag.2 0mi !i /diag.1=mi /X T M (2.14) with eigenfrequencies !i, modal masses mi, and the modal matrix Xgiven by the undamped system’s eigenvalue problem KXDMXdiag !2 i diag.mi / DX T MX (2.15) In a calibration procedure we are then able to search for the mass and stiffness related parameters p of the FE model fK(p), M(p)g that render the transfer function HA given by Eq. (2.4) and that let the criterion function of Eq. (2.5) to be minimal. The discrete frequencies used to evaluate Eq. (2.5) do not have to match the discrete frequencies used in testing. After a successful calibration, the FE model is in better agreement with test data. The eigenmodes from test and analysis are therefore normally well suited for mode-pair matching by correlation analysis. After such matching, the modal damping found in the experiment can easily be mapped to the finite element model’s viscous damping matrix using Eq. (2.14). 2.2.3 Surrogate Modeling To be practical, an FRF based model calibration requires rapid calculations of the frequency responses of the model as the parameter settings change in the iterative search for minimum deviation to test data. For big models this is undoable without model reduction. To obtain efficiency, a modal reduction scheme is here applied to create a surrogate model used by the calibration procedure. The eigenmodes of the corresponding undamped system at the nominal parameter configuration, belonging to all eigenvalues in a frequency range that significantly overlaps the frequency range of interest are then used for reduction. Let that range be ! D !low;!high . To save computational effort, that reduction basis is kept constant within one full calibration cycle and thus not modified as the parameter settings vary during minimization procedure iterations. Let the eigenvalue problem formulated at the nominal parameter setting p0 be K.p0/T DM.p0/T Ddiag !2 i 8!D !low;!high (2.16) Then the mass and stiffness matrices at any parameter settingp of the reduced model are M.p/ DT T M.p/T K.p/ DT T K.p/T (2.17)
RkJQdWJsaXNoZXIy MTMzNzEzMQ==