Nonlinear Dynamics, Volume 1

River Rapids Conference Proceedings of the Society for Experimental Mechanics Series Nonlinear Dynamics, Volume 1 Gaetan Kerschen Proceedings of the 36th IMAC, A Conference and Exposition on Structural Dynamics 2018 River Publishers

Conference Proceedings of the Society for Experimental Mechanics Series Series Editor Kristin B. Zimmerman, Ph.D. Society for Experimental Mechanics, Inc., Bethel, CT, USA

River Publishers Nonlinear Dynamics, Volume 1 Proceedings of the 36th IMAC, A Conference and Exposition on Structural Dynamics 2018 Gaetan Kerschen Editor

Published, sold and distributed by: River Publishers Broagervej 10 9260 Gistrup Denmark www.riverpublishers.com ISBN 978-87-7004-965-8 (eBook) Conference Proceedings of the Society for Experimental Mechanics An imprint of River Publishers © The Society for Experimental Mechanics, Inc. 2019 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 Nonlinear Dynamics represents one of nine volumes of technical papers presented at the 36th IMAC, A Conference and Exposition on Structural Dynamics, organized by the Society for Experimental Mechanics, and held in Orlando, Florida, February 12–15, 2018. The full proceedings also include volumes on Dynamics of Civil Structures; Model Validation and Uncertainty Quantification; Dynamics of Coupled Structures; Special Topics in Structural Dynamics; Structural Health Monitoring, Photogrammetry and DIC; Rotating Machinery, Vibro-Acoustics and Laser Vibrometry; Sensors and Instrumentation, Aircraft/Aerospace and Energy Harvesting; and Topics in Modal Analysis and Testing. Each collection presents early findings from experimental and computational investigations on an important area within structural dynamics. Nonlinearity is one of these areas. The vast majority of real engineering structures behave nonlinearly. Therefore, it is necessary to include nonlinear effects in all the steps of the engineering design: in the experimental analysis tools (so that the nonlinear parameters can be correctly identified) and in the mathematical and numerical models of the structure (in order to run accurate simulations). In so doing, it will be possible to create a model representative of the reality which, once validated, can be used for better predictions. Several nonlinear papers address theoretical and numerical aspects of nonlinear dynamics (covering rigorous theoretical formulations and robust computational algorithms) as well as experimental techniques and analysis methods. There are also papers dedicated to nonlinearity in practice where real-life examples of nonlinear structures will be discussed. The organizers would like to thank the authors, presenters, session organizers, and session chairs for their participation in this track. Liege, Belgium Gaetan Kerschen v

Contents 1 Interface Reduction on Hurty/Craig-Bampton Substructures with Frictionless Contact....................... 1 Patrick J. Hughes, Wesley Scott, Wensi Wu, Robert J. Kuether, Matthew S. Allen, and Paolo Tiso 2 Experimental Path Following of Unstable Static Equilibria for Snap-Through Buckling ...................... 17 T. van Iderstein and R. Wiebe 3 Direct Detection of Nonlinear Modal Interactions and Model Updating Using Measured Time Series ....... 23 Keegan Moore, Mehmet Kurt, Melih Eriten, D. Michael McFarland, Lawrence A. Bergman, and Alexander F. Vakakis 4 Pareto Optimization of a Nonlinear Tuned Mass Damper to Control Vibrations in Hand Held Impact Machines.................................................................................................................. 27 Seyed Milad Mousavi Bideleh and Viktor Berbyuk 5 Inverse Methods for Characterization of Contact Areas in Mechanical Systems ................................ 45 Matthew Fronk, Kevin Eschen, Kyle Starkey, Robert J. Kuether, Adam Brink, Timothy Walsh, Wilkins Aquino, and Matthew Brake 6 Experimental Characterization of a New Benchmark Structure for Prediction of Damping Nonlinearity ... 57 Aabhas Singh, Matteo Scapolan, Yuta Saito, Matthew S. Allen, Daniel Roettgen, Ben Pacini, and Robert J. Kuether 7 Nonlinear System Identification for Joints Including Modal Interactions ........................................ 79 Alexander H. Haslam, Gaurav Chauda, Nalik Kenia, Esther S. Rufat-Meix, Matthew S. Allen, Robert M. Lacayo, Malte Krack, and Matthew R. W. Brake 8 Performance of Nonlinear Modal Model in Predicting Complex Bilinear Stiffness ............................. 101 Benjamin R. Pacini, Wil A. Holzmann, and Randall L. Mayes 9 Low Order Nonlinear Dynamic Modelling of Fuel Supply Pipes .................................................. 113 Alberto Sanchez and Christoph W. Schwingshackl 10 System Identification to Estimate the Nonlinear Modes of a Gong ................................................ 121 Daniel Piombino, Matthew S. Allen, David Ehrhardt, Tim Beberniss, and Joseph J. Hollkamp 11 An Enhanced Static Reduction Algorithm for Predictive Modeling of Bolted Joints............................ 137 Iman Zare, Matthew S. Allen, and Emily Jewell 12 Time-varying Spectral Submanifolds: Analytic Calculation of Backbone Curves and Forced Response ..... 141 Thomas Breunung and George Haller 13 Operational Modal Analysis Based Stress Estimation in Friction Systems ....................................... 143 Marius Tarpø, Tobias Friis, Bruna Nabuco, Sandro Amador, Evangelos Katsanos, and Rune Brincker 14 Damping Estimation of Friction Systems in Random Vibrations .................................................. 155 Tobias Friis, Evangelos Katsanos, Sandro Amador, and Rune Brincker vii

viii Contents 15 System Identification of Jointed Structures: Nonlinear Modal Testing Vs. State-Space Model Identification............................................................................................................. 159 Maren Scheel, Gleb Kleyman, Ali Tatar, Matthew R. W. Brake, Simon Peter, Jean-Philippe Noël, Matthew S. Allen, and Malte Krack 16 Effect of Boundary Conditions on Finite Element Submodeling................................................... 163 Michael W. Sracic and William J. Elke 17 On Euler Buckling and Snap-Through ............................................................................... 171 Richard Wiebe, Mihaela Nistor, and Ilinca Stanciulescu 18 Solitons in Cyclic and Symmetric Structures......................................................................... 175 Filipe Fontanela, Aurelien Grolet, Loic Salles, and Norbert Hoffmann 19 Experimental and Numerical Nonlinear Modal Analysis of a Beam with Impact: Part I – Numerical Investigation.............................................................................................................. 179 F. Schreyer, S. Peter, and R. I. Leine 20 Experimental and Numerical Nonlinear Modal Analysis of a Beam with Impact: Part II – Experimental Investigation.................................................................................. 183 S. Peter, F. Schreyer, and R.I. Leine 21 The Effect of Non-Flat Interfaces On System Dynamics............................................................ 187 I. Lawal, S. Shah, M. Gonzalez-Madrid, T. Hu, C. W. Schwingshackl, and M. R. W. Brake 22 Investigating Modal Contributions Using a Galerkin Model ....................................................... 199 A. J. Elliott, A. Cammarano, and S. A. Neild 23 Acoustic Excitation of a Flanged Joint ................................................................................ 211 Trevor W. Jerome, Micah R. Shepherd, and Stephen A. Hambric 24 In Situ Measurements of Interfacial Contact Pressure During Impact Hammer Tests.......................... 225 B. Seeger, P. Butaud, M. V. Baloglu, F. Du, M. R. W. Brake, and C. W. Schwingshackl 25 An Improved Shape Reconstruction Methodology for Long Rod Like Structures Using Cosserat Kinematics- Including the Poisson’s Effect ........................................................................... 237 Mayank Chadha and Michael D. Todd 26 Computing Nonlinear Normal Modes of Aerospace Structures Using the Multi-harmonic Balance Method.................................................................................................................... 247 Christopher I. VanDamme, Ben Moldenhauer, Matthew S. Allen, and Joseph J. Hollkamp 27 Nonlinear Identification of an Aero-Engine Component Using Polynomial Nonlinear State Space Model ... 261 Samson B. Cooper, Koen Tiels, and Dario DiMaio 28 Curved Structures That Can Elastically Snap-Through............................................................ 275 Lawrence N. Virgin, Yue (Susan) Guan, Raymond H. Plaut 29 Experimental Analysis of Non-Linear Damping Performance in Composites Materials Thanks to Local Transduction-Dissipation Phenomenon..................................................................... 279 Mathieu Chevalier, Camille Bessaguet, Luis Quiroga-Cortès, Eric Dantras, and Guilhem Michon 30 Subspace-Based Identification of a Distributed Nonlinearity in Time and Frequency Domains ............... 283 D. Anastasio, S. Marchesiello, J. P. Noël, and G. Kerschen 31 Reduced Order Modelling for Non-linear Rotating Systems in ALE Formulation with Contact .............. 287 Tim Weidauer and Kai Willner 32 Damage Precursor Indicator for Aluminum 7075-T6 Based on Nonlinear Dynamics........................... 303 Robert A. Haynes, Ed Habtour, Todd C. Henry, Daniel P. Cole, Volker Weiss, Antonios Kontsos, and Brian Wisner 33 Application of Control-Based Continuation to a Nonlinear System with Harmonically Coupled Modes ..... 315 L. Renson, D. A. W. Barton, and S. A. Neild

Contents ix 34 Numerically Assessing the Relative Significance of Nonlinear Normal Modes to Forced Responses .......... 317 T. L. Hill, S. A. Neild, and A. Cammarano 35 Direct Frequency Domain Identification of Time Varying Systems................................................ 327 Lee Mazurek and Richard Christenson 36 Reduced-Order Modelling for Investigating Nonlinear FEM Systems ............................................ 335 I. Tartaruga, S. A. Neild, T. L. Hill, and A. Cammarano 37 Nonlinear Forced Response of a Composite Fan Blade Actuated by Piezoelectric Patches: Simulation and Testing ............................................................................................................... 351 Antoine Mabilia, Claude Gibert, Fabrice Thouverez, and Edouard De Jaeghere 38 Locating Nonlinearity in Mechanical Systems: A Dynamic Network Perspective ............................... 363 J. P. Noël, M. Schoukens, and P. M. J. Van den Hof 39 Tracing a Prescribed Force-Displacement Curve Using Topology Optimization................................. 369 Jongsuh Lee, Thibaut Detroux, and Gaetan Kerschen 40 Model Updating of a Wing-Engine Structure with Nonlinear Connections ...................................... 373 Mingming Song, Ludovic Renson, Jean-Philippe Noël, Babak Moaveni, and Gaetan Kerschen 41 Modal Analysis of Axially Deforming Rods with Isolated Lap Joints ............................................. 375 D. Dane Quinn 42 Identification of Nonlinear Viscoelastic Parameters Based on an Enhanced Oberst Beam Method........... 379 Kévin Jaboviste, Emeline Sadoulet Reboul, Nicolas Peyret, Gaël Chevallier, C. Arnould, and E. Collard 43 A General Framework for Time Domain Finite Element Analysis of Viscoelastically Damped Structures... 383 J.-F. Deü and L. Rouleau 44 Nonlinear Structural, Inertial and Damping Effects in an Oscillating Cantilever Beam........................ 387 Michal Raviv Sayag and Earl H. Dowell 45 Reduced Order Modeling of Structures with Preloaded Bolted Joints by the Use of Trial Vector Derivatives................................................................................................................ 401 Wolfgang Witteveen and Florian Pichler 46 Towards the Development of a Model for Nonlinear Elements in Machine Tools................................ 405 Steven M. Whitican, Charles Van Karsen, and Jason Blough 47 Nonlinear Characterization of a Machine Tool Energy Absorber ................................................. 419 Steven M. Whitican, Charles Van Karsen, and Jason Blough 48 Tutorial: Bolted Joints and Tribomechadynamics ................................................................... 427 M. R. W. Brake

Chapter 1 Interface Reduction on Hurty/Craig-Bampton Substructures with Frictionless Contact Patrick J. Hughes, Wesley Scott, Wensi Wu, Robert J. Kuether, Matthew S. Allen, and Paolo Tiso Abstract Contact in structures with mechanical interfaces has the ability to significantly influence the system dynamics, such that the energy dissipation and resonant frequencies vary as a function of the response amplitude. Finite element analysis is commonly used to study the physics of such problems, particularly when examining the local behavior at the interfaces. These high fidelity, nonlinear models are computationally expensive to run with time-stepping solvers due to their large mesh densities at the interface, and because of the high expense required to update the tangent operators. Hurty/Craig-Bampton substructuring and interface reduction techniques are commonly utilized to reduce computation time for jointed structures. In the past, these methods have only been applied to substructures rigidly attached to one another, resulting in a linear model. The present work explores the performance of a particular interface reduction technique (system-level characteristic constraint modes) on a nonlinear model with node-to-node contact for a benchmark structure consisting of two c-shape beams bolted together at each end. Keywords Component mode synthesis · Substructuring · Hurty/Craig-Bampton Method · Interface reduction · Mechanical joints · Contact 1.1 Introduction Contact in structures with mechanical interfaces can have a significant influence on the dynamic response to time varying loads. The magnitude of clamping forces near the contact due to bolt preloads influences the overall stiffness in the system, while frictional effects in joints typically contribute to the overall structural damping [1]. A common approach in interface modeling is to approximate contact areas with linear springs and dashpots, or with nonlinear elements such as Jenkins elements [2] or Iwan elements [3]. These approaches are capable of simplifying the interface model at the expense of losing local kinematics and stresses at these locations. When higher interface fidelity is desired, the interfaces must be modeled in full detail with a fine mesh resolution to adequately resolve nonlinear contact and friction effects. Sophisticated models of this type are computationally expensive to run with transient solvers, due to the high number of degrees-of-freedom (DOF) required to accurately capture the local contact forces at the interface. Furthermore, the contact state can change in time, requiring continuous updating of internal force vectors and Jacobian operators within implicit integration schemes. This work seeks to address this issue by exploring model order reduction techniques to speed up transient simulations for structures with nonlinear contact. The current study focuses on the frictionless case such that there are no tangential loads due to friction. P. J. Hughes Department of Structural Engineering, University of California – San Diego, San Diego, CA, USA W. Scott · M. S. Allen Department of Engineering Physics, University of Wisconsin-Madison, Madison, WI, USA W.Wu School of Civil and Environmental Engineering, Cornell University, Ithaca, NY, USA R. J. Kuether ( ) Sandia National Laboratories, Albuquerque, NM, USA e-mail: rjkueth@sandia.gov P. Tiso Institute for Mechanical Systems, Zurich, Switzerland © The Society for Experimental Mechanics, Inc. 2019 G. Kerschen (ed.), Nonlinear Dynamics, Volume 1, Conference Proceedings of the Society for Experimental Mechanics Series, https://doi.org/10.1007/978-3-319-74280-9_1 1

2 P. J. Hughes et al. Fig. 1.1 (a) Finite element model of bolted C-beam assembly with coordinate axes. (b) Close-up view of interface surface with bolt DOF spider Component mode synthesis (CMS) methods in structural dynamics are used to reduce the linear portion of a model, while preserving the physical DOF at interfaces containing nonlinear elements [4, 5]. One commonly used approach is the Hurty/Craig-Bampton (HCB) method, which was originally proposed by Hurty [6], and later simplified by Craig and Bampton [7]. The HCB method represents each subcomponent by a truncated set of dynamic fixed-interface modes, augmented with a static constraint mode for every physical interface DOF. For finite element models with a fine interfacial mesh, the number of constraint modes needed for the HCB method becomes excessively high and prohibits the effectiveness for applications involving contact. Furthermore, the critical time step for HCB models with explicit time solvers is dictated by the static constraint modes, which are localized in shape and therefore associated with high frequencies. As such, model reduction does little if anything to improve this [4]. The objective of this research is to further decrease the HCB model order by performing a secondary reduction on the interface DOF using a set of mode shapes describing the interface kinematics. Some of the authors have recently reviewed interface reduction methods for HCB models in [8], but these techniques were only applied to linear substructures with rigid compatibility enforced (i.e. linear assembly models). The review paper reveals that the system-level characteristic constraint mode method originally developed by Craig and Chang [9], and later elaborated by Castanier et al. [10] is most accurate for rigidly connected boundary DOF. This modal basis is slightly modified in this work to preserve some of the physical boundary DOF in the subspace. A few other works have explored interface reduction techniques on contacting interfaces. Becker and Gaul [5] applied component mode synthesis to structure with bolted joints, but the number and type of interface modes needed to resolve the local and global response remains an open research question. In this research, the HCB method with a modified version of system-level interface reduction [9, 10] is used to model a system consisting of two C-shape beams bolted together at each end (see Fig. 1.1 in Sect. 1.3.1). Interface contact is considered frictionless here, modeled using node-to-node penalty springs [11] in the normal direction only. This work examines the necessary interface basis vectors to resolve the nonlinear kinematics at the frictionless contacting surfaces. Furthermore, this research examines the effect of interface reduction on transient solution accuracy and simulation time, relative to the full-interface HCB model. 1.2 Theory 1.2.1 Nonlinear Hurty/Craig-Bampton Method Modern finite element models tend to have extremely fine meshes, with potentially hundreds of thousands of elements and millions of DOF. In models with substructures connected to one another via contacting boundary conditions, the interface DOF may control the dynamic response of the system, while the numerous interior (non-interface) DOF provide unnecessary model redundancy. As such, the interior DOF can be reduced through CMS methods, which approximate the substructure interior with a relatively small set of mode shapes, while leaving the interface DOF unchanged. One such method is the Hurty/Craig-Bampton technique, summarized below for the case when substructures are connected through nonlinear contact elements. Consider an arbitrary substructure in an assembly with nonlinear elements at the interface. In discretized form, the equations of motion for the substructure are written as MRuCKuCfN Dfext (1.1)

1 Interface Reduction on Hurty/Craig-Bampton Substructures with Frictionless Contact 3 where Mand Kare respectively the mass and stiffness matrices, u is the relative displacement vector, fN is the vector of displacement-dependent normal contact forces (discussed in Sect. 1.2.3), andfextis the externally applied loading vector. The matrices are n n and the vectors are n 1. Each dot placed over u represents its derivative with respect to time. These equations can be partitioned into interior i and interface b DOFas Mii Mib Mbi Mbb Rui Rub C Kii Kib Kbi Kbb ui ub C 0 fN;b D fext;i fext;b (1.2) The mass and stiffness from the interior degrees of freedom are used to form the fixed-interface (FI) modes of the system through a partitioned eigenanalysis, i.e. hKii ¨ FI s 2Miii® FI s D0 (1.3) The eigenvectors ®FI s (s D1, 2, : : : , nFI) form the FI modal matrix ˆ FI D ® FI 1 ® FI 2 ® FI nFI (1.4) wherenFI is an integer much smaller than the original number of DOF in set i. The truncation of modes in this step introduces the model reduction desired by the HCB process. The eigenshapes are combined with a set of static constraint modes that describe how the interior DOF respond when each interface DOF is moved independently. This is defined by the static condensation ‰ HCB D K 1 ii Kib (1.5) The combination of FI modes and constraint modes form a basis for the interior DOF of the substructure, and are arranged to form the HCB transformation matrix as THCB D ˆ FI ‰ HCB 0 I (1.6) which converts the interior DOF of the substructure to generalized HCB coordinates through the transformation u DTHCBv (1.7a) ui ub D ˆ FI ‰ HCB 0 I qi ub (1.7b) The vector v is the HCB generalized coordinate and qi is the modal representation of the interior partition of the substructure. The reduced mass and stiffness matrices, as well as the reduced load vectors, are computed via MHCB D THCB TMTHCB ; KHCB D THCB TKTHCB (1.8) fHCB N D THCB Tf N; f HCB ext D THCB Tf ext where T is the transpose operator. The dimension of the problem after HCB reduction is n FI Cnb, where nFI is the number of retained fixed-interface modes, and nb is the total number of interface DOF. The equations of motion for a substructure in HCB coordinates are MHCBRvCKHCBvCfHCB N Df HCB ext (1.9a) Iii MHCB ib MHCB bi MHCB bb Rqi Rub C ƒFI ii 0 0 KHCB bb qi ub C 0 fHCB N;b D ( fHCB ext;i fHCB ext;b ) (1.9b) where Iii and ƒFI ii are the identity matrix and diagonal matrix of fixed-interface eigenvalues, respectively. This form arises from the orthogonality properties of the modal matrixˆFI in the HCB coordinate transformation.

4 P. J. Hughes et al. 1.2.2 Interface Modal Basis Any interface reduction basis must satisfy a number of criteria: (1) allow for realistic deformations in the contact area, (2) reproduce the distribution of contact forces, (3) match the overall dynamic response of the HCB model within a reasonable margin of error, and (4) be efficient enough to provide overall computational savings. This research attempts to meet each of these criteria by utilizing an interface reduction technique that is a modification of the system-level characteristic constraint modes [9, 10]. In addition to these mode shapes, the basis is augmented with the static deformation shape obtained from a nonlinear static preload analysis, as it is inefficient to resolve the interface kinematics with dynamic modes alone. 1.2.2.1 Extended System-Level Characteristic Constraint Modes About Preloaded State The traditional method of system-level characteristic constraint modes, or SCC method, transforms all physical interface DOF in the HCB model into truncated modal DOF, such that the resulting model is defined entirely in the modal domain. In some cases, it is desirable to retain some portion of the interface as physical DOF. For example, it may be necessary to retain physical DOF when loads are applied near substructure boundaries (e.g. bolt preload force). A novel method to achieve this is introduced here, and referred to as the extended system-level characteristic constraint mode method, or SCCe method. To reduce some of the interface DOF to modal DOF, and retain the rest as physical DOF, a given substructure must now be partitioned into three sets: the interior DOFi (i.e. fixed-interface mode partition), physical interface DOFp, and reduced interface DOFr. The HCB system mass and stiffness matrices, according to these partitions, are MHCB D 2 64 Iii MHCB ir MHCB ip MHCB ri MHCB rr MHCB rp MHCB pi MHCB pr MHCB pp 3 75; KHCB D 2 64 ƒFI ii 0 0 0 KHCB rr KHCB rp 0 KHCB pr KHCB pp 3 75 (1.10) The contact state of the model about its preloaded state is important to capturing the correct dynamic response, and hence the SCCe modal basis must account for this. This is achieved by linearizing the preloaded HCB model to include stiffness contributions from the contact areas that form upon fastening/preloading substructures coming into contact. The static deformation due to the preload force (denoted as ve in HCB coordinates) is solved using a nonlinear static solver (see Sect. 1.2.2.2). A linearized stiffness matrix (denoted with –c due to the added constraints) is formed by summing the linear HCB stiffness matrix and the Jacobian of the normal contact force, evaluated about the deformed state, KHCB c DKHCB C @fHCB N @v v Dve (1.11) A set of SCCe mode shapes is computed from a second partitioned eigenanalysis, employing the linearized stiffness matrix and isolating the reduced interface DOFr hKHCB c rr ¨ SCCe c 2MHCB rr i® SCCe c s D0 (1.12) A limited set of eigenvectors ®SCCe c s D(s D1, 2, : : : , nS) are assembled into the columns of the constrained SCCe modal matrix ˆ SCCe c D ®SCCe c 1 ®SCCe c 2 ®SCCe c ns (1.13) where nS is a number smaller than the original number of DOF in set r. The modes contained in ˆSCCe c provide a good representation of the local interface dynamics when the interface is in contact due to the additional term in the contact stiffness, KHCB c. If the interface loses contact during a vibration response, then these modes may be insufficient. In order to address this, the eigenanalysis can be repeated using the original stiffness matrixKHCB, such that a new set of unconstrained modes (denoted with –u) can be added to the reduction basis. hKHCB rr ¨ SCCe u 2MHCB rr i® SCCe u s D0 (1.14)

1 Interface Reduction on Hurty/Craig-Bampton Substructures with Frictionless Contact 5 The first nS unconstrained eigenvectors are used to build the unconstrained SCCe modal matrix as ˆ SCCe u D ®SCCe u 1 ®SCCe u 2 ®SCCe u ns (1.15) Both sets of eigenmodes described in Eqs. (1.13) and (1.15) use ther-r partition of the system matrices, so they inherently fix the motion of the physical interface DOF (set p) to be zero. To alleviate this issue, static constraint modes are added to the interface reduction set, similar to those in original the HCB method ‰ SCCe D KHCB c rr 1KHCB c rp (1.16) Combining the dual set of eigenmodes with the constraint modes, along with appropriately-sized zero and identity partitions, the SCCe transformation matrix is formed as TSCCe D2 4 I 0 0 0 0 ˆ SCCe c ˆ SCCe u ‰ SCCe 0 0 0 I 3 5 (1.17) This transformation uses a combination of dynamic (ˆSCCe c, ˆSCCe u) and static (‰SCCe) shapes to convert the r portion of the model to modal DOF, while retaining the p portion as physical DOF. The transformation between HCB coordinates and SCCe coordinates becomes v DTSCCew (1.18a) 8< : qi ur up 9= ; D2 4 I 0 0 0 0 ˆ SCCe c ˆ SCCe u ‰ SCCe 0 0 0 I 3 5 8ˆˆ< ˆ: qi qc qu up 9>>= >; (1.18b) where wis the SCCe generalized coordinate vector, and qc and qu are the constrained and unconstrained modal DOF in the r partition, respectively. 1.2.2.2 Static Preload Deformation The inclusion of the static preload shape ve in the reduction basis ensures that the SCCe model is exact for the initial preload analysis performed on the nonlinear HCB model in Eqs. (1.9a) and (1.9b). This shape is determined by solving the following equations for static equilibrium, stated in HCB coordinates KHCBve CfHCB N Db efe (1.19) where be is a mapping vector that positions and orients the scalar preload force amplitude fe at the correct physical DOF. In partitioned form, the equations are 2 64 ƒFI ii 0 0 0 KHCB rr KHCB rp 0 KHCB pr KHCB pp 3 75 8< : qe i ue r ue p 9= ; C 8ˆ< ˆ: 0 fHCB N;r fHCB N;p 9>= >; D8< : 0 0 be pf e 9= ; (1.20) Due to the nonlinearity introduced by the contact force fHCB N , the displacement vector ve D hqe i ue r ue p i T must be determined using an iterative solver. In this research, a damped Newton Raphson method was employed.

6 P. J. Hughes et al. 1.2.2.3 Total SCCe Transformation Matrix The static deformation vector is concatenated at the end of TSCCe to build the final SCCe transformation matrix TSCCeC TSCCeC D TSCCe ve (1.21a) TSCCeC D2 4 I 0 0 0 qe i 0 ˆ SCCe c ˆ SCCe u ‰ SCCe ue r 0 0 0 I ue p 3 5 (1.21b) This matrix performs the same transformation as Eqs. (1.18a) and (1.18b), but adds one more DOF from the preload vector. The HCB coordinate v is transformed to the new SCCe coordinates wCvia v DTSCCeCwC (1.22a) 8< : qi ur up 9= ; D2 4 I 0 0 0 qe i 0 ˆ SCCe c ˆ SCCe u ‰ SCCe ue r 0 0 0 I ue p 3 5 8ˆˆˆ < ˆˆ ˆ: qi qc qu up qe 9>>> = >> >; (1.22b) where qe is the scalar DOF introduced by the inclusion of ve. The variety of shapes that make upTSCCeCcan lead to a basis set that is not necessarily linearly independent. The singular value decomposition (SVD) reforms TSCCeCas TSCCeC DUSCCeC† SCCeCVSCCeC (1.23) where USCCeC and VSCCeC are the left and right singular vectors, respectively, and †SCCeC is a diagonal matrix of singular values. The vectors in USCCeC form an orthonormal basis for the column space of TSCCeC, so they are taken as the final transformation matrix to convert between HCB and SCCe coordinates. Thus, the SCCe system matrices and load vectors are computed as MSCCe D USCCeC TMHCBUSCCeC; KSCCe D USCCeC TKUSCCeC (1.24) fSCCe N D USCCeC TfHCB N ; f SCCe ext D USCCeC TfHCB ext The dimension of the problem is now nFI C2nS Cnp C1, where nFI is the number of retained fixed-interface modes, nS is the number of retained SCCe modes (one set of constrained, one set of unconstrained), np is the number of physical interface DOF, and the 1 accounts for the static preload shape. The equations of motion for the system in SCCe coordinates are MSCCe RwCCKSCCewCCfSCCe N Df SCCe ext (1.25a) 2 66 66 64 Iii MSCCe ic MSCCe iu MSCCe ip MSCCe ie MSCCe ci MSCCe cc MSCCe cu MSCCe cp MSCCe ce MSCCe ui MSCCe uc MSCCe uu MSCCe up MSCCe ue MSCCe pi MSCCe pc MSCCe pu MSCCe pp MSCCe pe MSCCe ei MSCCe ec MSCCe eu MSCCe ep MSCCe ee 3 77 77 75 8ˆˆˆ < ˆˆ ˆ: Rqi Rqc Rqu Rup Rqe 9>>> = >> >; C 2 66 66 64 ƒ FI ii 0 0 0 0 0 KSCCe cc KSCCe cu KSCCe cp KSCCe ce 0 KSCCe uc KSCCe uu KSCCe up KSCCe ue 0 KSCCe pc KSCCe pu KSCCe pp KSCCe pe 0 KSCCe ec KSCCe eu KSCCe ep KSCCe ee 3 77 77 75 8ˆˆˆ < ˆˆ ˆ: qi qc qu up qe 9>>> = >> >; C 8ˆˆˆ < ˆˆ ˆ: 0 fSCCe N;c fSCCe N;u fSCCe N;p fSCCe N;e 9>>> = >> >; D 8ˆˆˆ < ˆˆ ˆ: fSCCe ext;i fSCCe ext;c fSCCe ext;u fSCCe ext;p fSCCe ext;e 9>>> = >> >; (1.25b) Note that, although the normal contact force vector fSCCe N has been reduced to SCCe coordinates, it must still be computed in HCB coordinates, where all interface DOF are physically represented. In nonlinear dynamic simulations, this requires a

1 Interface Reduction on Hurty/Craig-Bampton Substructures with Frictionless Contact 7 transformation from SCCe to HCB coordinates at every time step to evaluate the nonlinear forces. The advantage that comes from this reduction is twofold: (1) the size of the equations to solve is smaller compared to the HCB model and (2) the high frequency content has been truncated and thus improves the critical time step in explicit solvers. 1.2.3 Normal Contact Model Contact in the normal direction is modeled using a node-to-node approach with the penalty method [11]. In this method, the normal gap g between a contacting node pair can be either negative (in contact) or positive (out of contact). For a particular node pair j, the normal gap is gj D Yj 1 C Yj 1 Yj 2 C Yj 2 (1.26) where the Y’s are the undeformed normal coordinates, and the Y’s are the physical relative displacements in the normal direction. The subscripts 1 and 2 refer to the two nodes in node pair j. The penalty method states that the normal contact force at node j fj N can be written as fj N D kpengj 0 ; gj <0 ; gj 0 (1.27) wherekpen is the penalty stiffness value. This defines the penalty force as proportional to the gap when the nodes are in contact (gj < 0), and zero when the nodes are out of contact (gj 0). The node-to-node contact model assumes that relative tangential displacements are small, and thus node pairs do not change during dynamic simulation. Construction of the normal contact force vector fN is a matter of positioning and orienting each nodal force f j N using a corresponding mapping vector bj N as fN D nn XjD1 bj Nf j N (1.28) where nn is the number of node pairs in all interfaces of the system. 1.2.4 Time Integration Schemes The response time histories of the reduced order models are computed using one of two different 2nd order accurate methods: an explicit scheme from Chung and Lee [12], and the implicit Hilber-Hughes-Taylor ’(HHT-’) method extended to nonlinear systems [13]. The implicit method is implemented with ’D0, yielding the constant average acceleration method, which is unconditionally stable for linear systems. The Chung-Lee scheme is conditionally stable, with a critical time step requirement that is dependent on its single free parameter, “CL, tcr D 2 ¨maxp4“CL 3 (1.29) where tcr is the critical time step length, and ¨max is the largest natural circular frequency in the system. The Chung-Lee parameter “CL is chosen to be 28/27, the maximum permissible value in the stability limit of 1 “CL 28/27. This selection improves stability by introducing slight numerical dissipation of high-frequency responses, which mitigates the amplification of numerical errors.

8 P. J. Hughes et al. Table 1.1 Summary of material properties used in finite element model of C-beams Symbol Description Numerical value E Young’s modulus 194 GPa Poisson’s ratio 0.290 ¡ Mass density 8000 Kg/m3 1.3 Numerical Studies: Two-Beam Assembly with Frictionless Contact 1.3.1 Model Description Figure 1.1a shows the benchmark finite element model under consideration, which consists of two identical C-beams bolted together at each end. Each individual C-beam is 50.8 cm long, 3.2 cm wide, 1.27 cm thick at its ends, and 0.95 cm thick at the mid-section. A total of 24,944 first order hexahedral elements are used to model each individual C-beam, and they are made of linear elastic structural steel. Table 1.1 provides a summary of material properties used for the C-beams. These are connected by 5/16-inch diameter steel bolts at each end that are modeled with 25 discretized beam elements along the y-axis, as depicted in Fig. 1.1b. The bolt ends are connected rigidly to a circular patch of nodes on the C-beams to represent clamping boundary conditions through a nut-washer-bolt assembly. The C-beam assembly has simply supported boundary conditions at all DOF along the z-directional line on the lower beam ends. The full fidelity finite element model has a total of 94,244 physical DOF. The model has two areas where the C-beams come into contact; each is 3.2 cm wide by 5.0 cm long, as indicated in Fig. 1.1b. The interface surfaces are meshed in such a way that their nodal locations are coincident and unmerged in the undeformed state, allowing for direct implementation of node-to-node contact elements. The Hurty/Craig-Bampton model is defined such that the boundary DOF include the nodes along these interface surfaces and the four nodes at each end of the beam elements (to apply the preload force); the remaining DOF are partitioned to interior DOF. Three different reduced order models are generated with fixed-interface mode frequencies cut off at 2, 3, and 4 kHz. The DOF count in the reduced model is dominated by the 3660 boundary DOF that account for the contacting surfaces, while 24 DOF pertain to the bolt DOF and 8– 16 DOF for the fixed-interface modes, depending on the cutoff frequency (see Sect. 1.3.3). The HCB reduction is performed using the Sierra Structural Dynamics (Sierra/SD) [14] finite element code developed at Sandia National Laboratories. Contact at the interface surfaces is considered through node-to-node penalty elements in the normal direction [11]. Selection of the penalty stiffness kpen is critical in the accuracy and stability of the penalty method - using a stiffness value that is too small allows excessive and non-physical nodal interpenetration, while a value that is too high can cause numerical ill-conditioning or instability. Preliminary studies indicated that a value of 3.78 109 N/m was sufficient to avoid excessive node overlap, while retaining stability and convergence. 1.3.2 Static Preloading of Bolt Elements A static preload force fe is applied to the ends of the bolt elements, to simulate the effect of bolt torque clamping the C-beams together. The preload force is obtained via an artificial strain"art, which leads to the following expression. fe DEA©art (1.30) The cross-sectional area A of each bolt is 4.95 10 5 m2, and the Young’s modulus E is the same as that of the C-beams (194 GPa). To ensure that a physically reasonable preload force is applied, a target internal bolt force is obtained using the DIN 946 standard [15] equation to convert applied torque T to transmitted axial force ftr. ftr D T 0:159PC0:578d T C0:5Df H (1.31) where P is the bolt thread pitch, d is the nominal bolt diameter, Df is the effective diameter in contact, T is the thread friction coefficient, and H is the head friction coefficient. Inputs for these variables are summarized in Table 1.2. For a typical applied torque of 18.5 ft-lb (25.1 N-m), the transmitted axial force is computed as 2.91 kN.

1 Interface Reduction on Hurty/Craig-Bampton Substructures with Frictionless Contact 9 Table 1.2 Summary of bolt preload parameters Symbol Description Numerical value T Applied bolt torque 25.1 N-m P Bolt thread pitch 0.00106 m d Nominal bolt diameter 0.00794 m Df Effective diameter in contact 0.0191 m T Thread friction coefficient 0.600 H Head friction coefficient 0.600 ftr Transmitted axial force 2.91 kN "art Artificial preload strain 3.32E-04 fe Preload force 3.19 kN Fig. 1.2 Amplified deformation after static preload, with close-up view of interface Table 1.3 HCB model comparison for different cutoff frequencies Cutoff FI frequency [Hz] FI modes per substructure Max FI frequency [Hz] Maxmodel frequency [Hz] Model size [DOF] tcr [sec] ta [sec] 2000 4 1544 9.14EC06 3692 3.25E-08 2.93E-08 3000 7 2534 9.13EC06 3698 3.25E-08 2.93E-08 4000 8 3763 9.13EC06 3700 3.25E-08 2.93E-08 The correct "art (and therefore f e) is obtained through an iterative procedure that incrementally alters "art, and applies the corresponding fe, until the resulting internal bolt axial force matches the theoretical transmitted force f tr. This procedure indicates that an artificial strain of 3.32 10 4 (corresponding to a preload force of 3.19 kN) is appropriate to produce an internal axial force equal to 2.91 kN. A visualization of the preloaded deformation mapped to the full-field model is displayed in Fig. 1.2, along with a close-up view of the interface, showing how the preload force indents the C-beam surfaces and causes receding of the contact area. 1.3.3 Results: Full-Interface Model A critical step in the Hurty/Craig-Bampton (HCB) reduction is determining the number of fixed-interface (FI) modes to retain, which is established by discarding FI modes with frequencies above a defined threshold. For a given loading bandwidth and pattern, the dynamic response of the system will converge to the “true” solution by increasing the frequency cut-off limit. The number of FI modes required for convergence is considered adequate to capture the important dynamic characteristics of the system. For this research, so-called “cutoff” frequencies of 2, 3, and 4 kHz are examined and summarized in Table 1.3. The lowest fixed-interface mode frequency in all cases is 288 Hz. Critical time step lengths are computed for the Chung-Lee Scheme [12] using Eq. (1.29). The maximum natural circular frequency ¨max was determined using the regular HCB mass matrix and the linearized HCB stiffness matrix computed in Eq. (1.11). This method of computing the critical time step assumes a linear model, so the analysis time step ta is taken as 90%of tcr to ensure stability in the nonlinear case. Table 1.3 demonstrates that the critical time step length is controlled by the HCB constraint modes computed in Eq. (1.5), and remains unchanged as fixed-interface modes are truncated. This suggests that HCB substructuring with nonlinear

10 P. J. Hughes et al. Fig. 1.3 (a) Fundamental mode shape of the HCB model. (b) Dynamic loading pattern interface DOF does not gain any computational advantage directly from explicit time integration, except for the lower order of equations solved. In other words, the critical time step is determined by the constraint modes. An externally applied loading pattern fext(t) was designed to excite the fundamental mode of the C-beam model, which has the displacement field shown in Fig. 1.3a. This mode exhibits bending about the z-axis, with most of the deformation concentrated in the top beam. To excite this mode, a vertical haversine impulse was applied to a strip of nodes at center of the top C-beam. The impulse has a duration of 1 ms and was applied with various amplitudes. The loading pattern can be seen in Fig. 1.3b. Dynamic simulation results are presented for the three HCB models described in Table 1.3, and for three loading amplitudes of 15, 100, and 500 N. These loading amplitudes are selected to elicit linear, moderately nonlinear, and strongly nonlinear responses, respectively. It is important to study different loading amplitudes to ensure that the HCB models retain accuracy for varying degrees of nonlinear response. A viscous damping matrix with 1% damping in all modes is included in the system to simulate a more realistic system. Response convergence is determined based on system-level and interface-level outputs: the drive point displacement (vertical displacement at the loading point) is used to evaluate the global response convergence, while contact area and summation of normal contact forces determine local interface-level convergence. Time histories are computed out to 10 ms to adequately show the transient vibration response of the system subjected to a 1 ms haversine impulse. A comparison of the different loading amplitudes and HCB models in Fig. 1.4 reveals that increasing the cutoff frequency from 2 to 3 kHz significantly alters the response at the system-level and interface-level. Drive point displacements, as well as contact area and interface forces, change dramatically in magnitude and frequency between the 2 and 3 kHz models. Furthermore, a visual comparison shows strong similarity between the 3 and 4 kHz results, indicating the reduced order models have converged. Similar to results from [4], the 2 kHz model is not accurate for this loading because the 2/£ impulse frequency is 2000 Hz, but the model only kept frequencies up to 1544 Hz. Thus, the 2 kHz model cannot produce response frequencies in the range generated by the 1 ms impulse. It is interesting to observe that for the 3 and 4 kHz models, the qualitative behavior of the drive point displacement does not change between different input levels (although the amplitude is scaled proportional to force). The contact area changes significantly with force level, and the total force at the interface oscillates considerably about the 2.91 kN preload force level. A seemingly linear global response in fact contains nonlinearities that are hidden in the interface-level outputs. This implies that, should friction be included, the damping characteristics of the system would change dramatically in time as the interfaces come in and out of contact. The 3 and 4 kHz models provide satisfactory starting points for interface reduction, and both would be suitable for the dynamic loading considered here. Nonetheless, the 4 kHz model provides roughly 1200 Hz of increased frequency bandwidth over the 3 kHz model at the cost of only two more DOF, so the 4 kHz model is used for the following interface reduction studies. 1.3.4 Results: Interface-Reduced Model The dynamic analysis from the full-interface studies in Sect. 1.3.3 is repeated here, but now using the interface-reduced models described by Eqs. (1.25a) and (1.25b) The base HCB model includes all fixed-interface frequencies up to 4 kHz,

1 Interface Reduction on Hurty/Craig-Bampton Substructures with Frictionless Contact 11 (a) (b) (c) 2 kHz cutoff 3 kHz cutoff 4 kHz cutoff Fig. 1.4 Time history of (a) drive point displacement; (b) contact area (percentage of nodes in contact); (c) total normal contact force on one interface. Comparison between different HCB cutoff frequencies Fig. 1.5 Amplified SCCe mode shapes. Constrained modes are formulated in Eqs. (1.12) and (1.13). Unconstrained modes are formulated in Eqs. (1.14) and (1.15) which equates to 16 FI modes. The effect of interface fidelity is considered through several reduction cases that retain between 20 and 3659 total SCCe modes (constrained Cunconstrained). A few SCCe mode shapes with their corresponding frequencies are plotted in Fig. 1.5.

12 P. J. Hughes et al. Table 1.4 Model comparison for various levels of SCCe interface reduction SCCemodes retained Max constrained SCCe frequency [Hz] Max unconstrained SCCe frequency [Hz] Maxmodel frequency [Hz] Model size [DOF] tcr [sec] ta [sec] 20 1.95EC03 1.90EC03 1.87EC05 61 2.51E-07 2.26E-07 100 5.01EC04 4.18EC04 4.28EC05 141 1.37E-07 1.23E-07 200 1.12EC05 1.06EC05 5.88EC05 241 9.54E-08 8.59E-08 500 2.46EC05 2.23EC05 8.91EC05 541 6.46E-08 5.81E-08 993* 3.77EC05 3.48EC05 1.19EC06 1034 4.37E-08 3.93E-08 3659* 8.48EC05 7.97EC05 1.63EC06 3700 3.26E-08 2.93E-08 Cases marked with * had some modes removed during the orthonormalization process The same process shown in Sect. 1.3.3 is employed to determine ¨max and tcr, but now the eigenvalue problem is formulated in the SCCe system using MSCCe and a linearized SCCe stiffness matrix KSCCe c, obtained via KSCCe c D USCCeC TKHCB cUSCCeC (1.32) Similar to Sect. 1.3.3, the analysis time step ta is obtained by multiplying the critical time step tcr by a factor of 0.9. Table 1.4 summarizes the different SCCe models in terms of maximum frequency, number of DOF, and critical time step length. Note that the analysis time steps of the SCCe models are significantly larger than that of the HCB model (2.93E-08), which results in nearly an order of magnitude speed-up in explicit time solvers. Drive point displacement, contact area, and normal contact force time histories for the different SCCe models are plotted against those of the 4 kHz HCB model, shown in Fig. 1.6. Visual comparison of the drive point displacement histories in Fig. 1.6a shows that all SCCe solutions lie on top of the HCB solution, with no observable deviations. For 500 or more SCCe modes, the normal contact force in Fig. 1.6c also shows near-perfect agreement. The contact area (percentage of nodes in contact) in Fig. 1.6b is the slowest to converge to the HCB solution, with significant errors observed for the case of 20 and 100 SCCe modes. All results converge to the HCB solution when all SCCe modes (3659 C1 augmentation vector) are retained and the transformation basis is not truncated. Drive point displacements from the HCB and SCCe models match well because the loading excites vibrational modes with deformations concentrated away from the interfaces, regardless of loading amplitude. The contact area and contact force, however, show a strong dependence on the loading amplitude. Errors in these quantities are negligible for the linear excitation of 15 N, but become more significant as the loading increases to the nonlinear excitations of 100 and 500 N. This implies that the accuracy of the SCCe method not only depends on the number of modes in the reduction basis, but also on the degree of nonlinearity (and therefore the loading amplitude). If the only output of interest is some system-level displacement at a location sufficiently distant from the interfaces, then only a small number of SCCe modes is required. Still, the rapid changes in contact area (Fig. 1.6b) would have a significant effect on the system-level damping if friction was included at the interface. In that case, more SCCe modes would be necessary to accurately capture system-level response. It is important to note that at certain levels of interface reduction, the total contact force (Fig. 1.6c) may be accurate, while the contact area (Fig. 1.6b) has significant errors. Interface forces in the SCCe model with 500 modes show near-perfect agreement with the HCB model, but the contact areas are in error by up to 20%. A relatively small number of SCCe modes are needed for convergence in overall contact forces, but more modes must be added to achieve suitably accurate interfacial stresses, which require accurate representation of the contact area. This is confirmed by the distribution of contact forces plotted in Fig. 1.7. The spread of contact forces over one interface surface is plotted at single instant in time (t D3.10 ms), and under the most severe loading (500 N). At this selected time, the structure is transitioning from low to high contact area, a phenomenon that seemingly requires a high number of included SCCe modes. The sum of instantaneous contact forces quickly converges to around 3500 N, but doesn’t accurately represent the distribution of those forces within the interface until many more modes are added. To the extent previously described, the interface reduction basis provides an acceptable level of accuracy relative to the full-interface HCB model. The final and arguably most important metric for gauging the effectiveness of the interface reduction is its potential for computational savings. The previous analyses were repeated using the implicit and explicit time integration schemes described in Sect. 1.2.4. Simulation run times for each of these cases are plotted in Figs. 1.8 and 1.9. The solve time reduction factor is the ratio of HCB solve time to SCCe solve time. Explicit time integration is carried out using 90% of the critical time step as determined by Eq. (1.29). This resulted in a time step length of 2.93E-08 sec for the HCB model, and ranged from 2.26E-07 sec to 2.93E-08 sec for the SCCe models (see Table 1.4). The implicit scheme is unconditionally stable, so a time step length is chosen provide adequate resolution

RkJQdWJsaXNoZXIy MTMzNzEzMQ==