5w. re... , «git. an. #th , t v. t ...)«.‘47. 11+). p i??? . 13.. r. ”.33. 9?? an .- . . I » ‘ ‘ , . y. \Znflfifiiv‘ t , . , . . , , ‘ l ‘ 12:1? I. I l. . v VJ... I134 . LIBRARY Michigan State University This is to certify that the dissertation entitled MODEL REDUCTION OF NONLINEAR STRUCTURAL SYSTEMS USING NONLINEAR NORMAL MODES AND COMPONENT MODE SYNTHESIS presented by Polarit Apiwattanalunggarn has been accepted towards fulfillment of the requirements for the PhD. degree in Mechanical Ermineering ii” wit ' U Major Rififessor’s Signature K/ZI/oz I I Date MSU is an Affirmative Action/Equal Opportunity Institution PLACE IN RETURN BOX to remove this checkout from your record. To AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 6/01 c:/ClRC/DateDue.065-p. 15 MODEL REDUCTION OF NONLINEAR STRUCTURAL SYSTEMS USING NONLINEAR NORMAL MODES AND COMPONENT MODE SYNTHESIS By Polarit Apiwattanalunggarn A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mechanical Engineering 2003 ABSTRACT MODEL REDUCTION OF NONLINEAR STRUCTURAL SYSTEMS USING NONLINEAR NORMAL MODES AND COMPONENT MODE SYNTHESIS By Polarit Apiwattanalunggarn This work addresses the general problem of model size reduction for describing the nonlinear vibration of structural elements and systems. The aim is to provide computational tools that allow one to accurately capture nonlinear dynamic behavior using a minimal number of degrees of freedom. In typical applications the finite element (FE) method is used to generate structural dynamic models, and model size reduction is carried out using linear modal analysis with truncation. However, in some cases one must retain many modes in order to accurately capture essential nonlinear coupling between the linear modes. In this work we utilize nonlinear normal modes (NNMs) defined in terms of invariant manifolds for the purposes of model size reduction, since it directly addresses modal coupling. This approach, which makes use of master and slave modes, along with the concept of dynamic invariance, allows one to generate accurate reduced order models (ROM) with only a few DOF, while capturing the effects of all modeled linear modes without directly simulating them. There are three main contributions of the present effort: (1) Two new numerical approaches for solving the invariant manifold equations are introduced. Both approaches employ master modal displacement and velocity coordinates and are based on weighted-residual techniques. When compared with previous methods that utilize amplitude and phase variables, the new methods are found to be superior in terms of computational time but inferior in terms of accuracy. (2) A specific application is considered: the finite amplitude vibrations of a rotating beam, which is a crude model for a rotorcraft blade. This system is known to possess essential nonlinear coupling between axial and transverse displacements, thereby leading to slow modal convergence. The proposed method systematically captures this coupling and provides an accurate single degree of freedom ROM. These results demonstrate the utility of NNM-based ROM, since they combine the versatility of the finite element method with the accurate NNM model reduction technique. (3) A model reduction technique suitable for structures that can be partitioned into substructures is developed. This allows one to build ROMs using NNMs at the substructure level and to assemble these using a component mode synthesis (CMS) technique. It is found that the proposed nonlinear CMS technique generally provides an accurate model only when the couplings between substructures are weak. To my father and mother ACKNOWLEDGMENTS I would like to express my gratitude and appreciation to many people who have contributed to this work. Without them, it could not be finished. I would first like to thank my advisor, Professor Steve Shaw, for his generous support, patience and guidance throughout my research work. I would also like to thank my co- advisor, Professor Christophe Pierre, from the University of Michigan for his helpful insight on Component Mode Synthesis. I would also like to thank Professors Alejandro Diaz, Ronald Averill, Milan Miklavcic, and Brian Feeny for serving in my dissertation committee. I would like to acknowledge the contributions of professors who taught me the knowledge in engineering and mathematics. Professor Steve Shaw deserves thanks for teaching me solid foundations in system dynamics, vibrations and control and most importantly training me to use my physical insight to cross check with mathematical results. Professors Brian Feeny, Alan Haddow, Ronald Rosenberg, Clark Radcliffe, Hassan Khalil, N ing Xi, and Fathi Salem also deserve thanks for their lectures in system dynamics, vibrations and control. I would also like to thank Professors Ronald Averill, Patrick Kwon, Hungyu Tsai, and Darren Ma- son for their lectures in finite element method and mechanics. I wish to thank Professor Abraham Engeda for allowing me to visit his class ME802 Advanced Classical Thermody- namics. I also wish to thank Professors Milan Miklavcic, Jerry Schuur, John McCarthy, Patricia Lamm, and Michael Hazier for expanding my knowledge in mathematics. I would like to thank the colleagues from the University of Michigan, Dr. Eric Pesheck, a former Ph.D. student, and Mr. Dongying Jiang, a Ph.D. student, for their suggestions, supplying some results, and allowing me to use their computational programs to generate some results in this work. Thanks also goes to the colleagues in the System Dynamics and Vibration Laboratory, Dr. Choong-Min Jung, a former Ph.D. student, Brian Olson, Tyler Nester, Peter Schmitz, Jeffrey Quinby, and Jeff Rhoads, for their suggestions, comments, and friendship. I wish to thank the Royal Thai Government and Kasetsart University for giving me the opportunity and support to pursue higher education at l\x’lichigan State University. I would also like to thank a friend of mine from Thailand, Mr. l\-’Ianoch Srinangyam, who I met during my freshman year at Chulalongkorn University and travelled with me on the first trip to the U.S., for being a good friend and giving helpful suggestions during my study at Michigan State University. I wish to express my thanks to the special person, Dr. Kunlakarn Lekskul, who has walked in to my life for seven years. It has been the joyful time to get to know her during all these years. I could not possibly go through many difficulties without her encouragement and support. Finally, I would like to express my special thanks and love to my mother who always values the importance of education. Her farsightedness is the basis for my success today. I also thank my sister and brothers for supplying all my needs from Thailand. vi TABLE OF CONTENTS LIST OF TABLES ix LIST OF FIGURES x 1 Introduction 1 1.1 Motivation .................................... 1 1.2 Nonlinear Normal Modes ............................ 2 1.3 Beam (Blade) Dynamics ............................ 7 1.4 Substructure Synthesis ............................. 9 1.4.1 Component Mode Synthesis ...................... 9 1.5 Dissertation Organization ........................... 11 2 Comparison of Methods for Constructing Invariant Manifolds for Non- linear Normal Modes 13 2.1 Introduction ................................... 13 2.2 Formulations of N NM Invariant Manifold ................... 15 2.2.1 Invariant Manifold Equations in Terms of Modal Position and Velocity 16 2.2.2 Invariant Manifold Equations in Terms of Modal Amplitude and Phase 22 2.3 Solution of the Invariant Manifold Equations ................ 26 2.3.1 Solution by Asymptotic Series Expansion ............... 27 2.3.2 Solution by Galerkin and Collocation Methods ............ 31 2.3.3 Comparison of Computational Efforts ................. 47 2.4 Examples .................................... 53 2.4.1 A Two-DOF Nonlinear Spring-Mass System ............. 54 2.4.2 A Finite-Element-Based Model of a Rotating Beam ......... 58 2.5 Conclusions ................................... 62 2.6 Tables ...................................... 63 2.7 Figures ...................................... 64 2.8 Appendix 2A .................................. 80 2.9 Appendix 28 .................................. 81 3 Finite-Element-Based Nonlinear Modal Reduction of a Rotating Beam with Large-Amplitude Motion 85 3.1 Introduction ................................... 85 3.2 Rotating Blade Formulation .......................... 88 vii 3.3 Modal Reduction ................................ 94 3.3.1 Galerkin-based Invariant Manifold Approach ............. 94 3.3.2 Collocation-based Invariant Manifold Approach ........... 99 3.4 Results ...................................... 100 3.4.1 Linear System Convergence ...................... 100 3.4.2 Nonlinear Model Development ..................... 101 3.4.3 System Response ............................ 103 3.4.4 Computational Considerations ..................... 106 3.5 Conclusions ................................... 107 3.6 Tables ...................................... 108 3.7 Figures ...................................... 109 3.8 Appendix 3 ................................... 121 4 Component Mode Synthesis Using Nonlinear Normal Modes 128 4.1 Introduction ................................... 128 4.2 Formulation of the NNM Invariant Manifolds ................. 130 4.2.1 Invariant-Manifold Governing Equations ............... 130 4.2.2 Invariant Manifold Solution ...................... 133 4.3 Fixed-Interface Nonlinear Component Mode Synthesis ............ 136 4.3.1 System Representation in Physical Coordinates ........... 137 4.3.2 System Representation in Linear Modal Coordinates ........ 138 4.3.3 Synthesis with Nonlinear Modal Reduction .............. 143 4.4 A Five-DOF Nonlinear Spring-Mass System ................. 147 4.5 A Class of Systems ............................... 151 4.5.1 Model Development ........................... 152 4.5.2 Example ................................. 154 4.6 Conclusions ................................... 165 4.7 Tables ...................................... 166 4.8 Figures ...................................... 168 4.9 Appendix 4 ................................... 198 5 Conclusions and Erture Works 205 BIBLIOGRAPHY 211 viii 2.1 2.2 3.1 3.2 4.1 4.2 4.3 4.4 LIST OF TABLES Comparison of the total number of unknown coefficients and the number of evaluation points needed for the integrands in order to compute the integrals for the global domain approaches ........................ Comparison of the total number of unknown coefficients and the number of evaluation points needed for the integrands in order to compute the integrals for the local domain approaches ......................... Dimensionless linear natural frequencies. (.01 and tag were computed using results from [1]. (123 and (124 were computed using results from [2] ....... Comparison of the linear natural frequencies obtained by power series and finite element methods. The analytical results for an and tag were computed using results from [1]. The analytical results for rug and M] were computed using results from [2] ............................... Fixed—interface linear natural frequencies of substructure [3 using the second— trial parameters. ................................ Fixed-interface linear natural frequencies of substructure 7 .......... Comparison of the linear natural frequencies of the three—DOF linear-CMS nonlinear model and the three—DOF nonlinear-CMS model .......... Comparison of the first-sixteen linear-natural frequencies of the forty-one- DOF linear-CMS nonlinear model and the forty-one-DOF original model. . 63 63 108 108 166 166 167 167 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 2.10 2.11 2.12 LIST OF FIGURES A two-DOF nonlinear spring-mass system. .................. The contribution of the second linear mode position 712 = P2(a1, m) to the first nonlinear mode manifold obtained by the global a — (b approach. The contribution of the second linear mode position 02 : X2071, in) to the first nonlinear mode. manifold obtained by the global u — v approach with Ub =1.5 and Vb =16. ............................. The contribution of the second linear mode position 772 = X2(n1(a1,¢1),171(a1, Q51» to the first nonlinear mode manifold obtained by the global u - v approach with Ub = 1.5 and Vb = 1.5. ........... The contribution of the second linear mode position 772 = X 2071,51) to the first nonlinear mode manifold obtained by the global u — v approach with Ub = 1.7 and Vb = 1.17 .............................. The contribution of the second linear mode position 172 = X2(7)1(a1, $1),7'71(a1,gb1)) to the first nonlinear mode manifold obtained by the global u — v approach with Ub = 1.7 and Vb = 1.17 ............ The contribution of the second linear mode position 712 = X2(7)1,1'71) to the first nonlinear mode manifold obtained by the global u — v approach with Ub = 2.6 and Vb = 1.79 .............................. The contribution of the second linear mode position 772 = X2(7)1(a1, r151),r)1(a1, 951)) to the first nonlinear mode manifold obtained by the global u — v approach with Ub = 2.6 and Vb = 1.79 ............ The time responses of the first nonlinear mode displacement for the ROM of the global u — 12 approach with 111(0) 2 1.3, 7'}1(0) = 0.0, the ROM of the global (L — (15 approach with (11(0) 2 1.3, ¢1(0) = 0.0, and the original model with initial conditions from the shooting algorithm. ............. The time responses of the second linear mode displacement n2 on the first nonlinear mode manifold for the ROM of the global u — 21 approach, the ROM of the global a — d) approach, and the original model .......... The time responses of the displacement $1 from Figure 2.1 for the ROM of the global u — v approach, the ROM of the global a — d) approach, and the original model. ................................. The contribution of the second linear mode position 02 = P2(al,el) to the first nonlinear mode manifold obtained by the local a — go approach. . . . . 64 64 65 65 66 66 67 67 68 68 69 69 2.13 2.14 2.15 2.16 2.17 2.18 2.19 2.20 2.21 2.22 2.23 2.24 2.25 2.26 2.27 2.28 The contribution of the second linear mode position 772 = X2(r}1,7')1) to the first nonlinear mode manifold obtained by the local u — v approach ..... The time responses of the first nonlinear mode displacement for the ROM of the local u — 22 approach with 771(0) = 1.5,7'}1(0) = 0.0, the ROM of the local 0. — ¢ approach with a1(0) = 1.5, (1)1(0) : 0.0, and the original model with initial conditions from the shooting algorithm. ............. The time responses of the second linear mode displacement 712 on the first nonlinear mode manifold for the ROM of the local u — '11 approach, the ROM of the local (I — (15 approach, and the original model .............. The time responses of the displacement $1 from Figure 2.1 for the ROM of the local u — v approach, the ROM of the local 0. — (i) approach, and the original model. ................................. Finite element representation of a rotating beam with Q = constant. . . . . The contribution of the second linear flapping mode 772 = P2(a1,¢1) to the first nonlinear mode manifold obtained by the local a —- d) approach ..... The contribution of the fifth linear flapping mode 02 = P5(a1,¢1) to the first nonlinear mode manifold obtained by the local (1. — qb approach ..... The contribution of the first linear axial mode 7710 = P10(a1,¢)1) to the first nonlinear mode manifold obtained by the local (1 — ab approach. ...... The contribution of the sixth linear axial mode 015 = P15(a1,¢1) to the first nonlinear mode manifold obtained by the local (2 — 45 approach. ...... The contribution of the second linear flapping mode 1);; = X2(n1,r')1) to the first nonlinear mode manifold obtained by the local u — 11 approach ..... The contribution of the fifth linear flapping mode '02 2 X5(n1,7°)1) to the first nonlinear mode manifold obtained by the local u — v approach. . . . . The contribution of the first linear axial mode 7710 = X10(n1,1'71) to the first nonlinear mode manifold obtained by the local u — v approach ........ The contribution of the sixth linear axial mode 7715 = X15071, 7'71) to the first nonlinear mode manifold obtained by the local u — v approach ........ The time responses of the first nonlinear flapping mode displacement m for the ROM of the local u — v approach with 171(0) = 2.375, {71(0) 2 0.0, the ROM of the local 0. — (i) approach with (11(0) = 2.375, ¢1(0) = 0.0, and the original model with initial conditions from the shooting algorithm ...... The time responses of the second linear flapping mode displacement in on the first nonlinear mode manifold for the ROM of the local u - 12 approach, the ROM of the local a — (23 approach, and the original model. ....... The time responses of the fifth linear flapping mode displacement 775 on the first nonlinear mode manifold for the ROM of the local u — v approach, the ROM of the local (1 — g5 approach, and the original model. ......... xi 70 70 71 71 72 72 73 73 74 74 76 77 2.29 2.30 2.31 2.32 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 The time responses of the first linear axial mode displacement 7710 on the first nonlinear mode manifold for the ROM of the local u — 12 approach, the ROM of the local a — d) approach, and the original model. ......... 78 The time responses of the sixth linear axial mode displacement 7715 on the first nonlinear mode manifold for the ROM of the local u — 12 approach, the ROM of the local a — (35 approach, and the original model. ......... 78 The time dependence of the beam—tip flap deflection, w(L, t), for the ROM of the local u — 21 approach, the ROM of the local a — qt) approach, and the original model. ................................. 79 The time dependence of the beam—tip axial deflection, u(L, t), for the ROM of the local u — v approach, the ROM of the local a — d) approach, and the original model. ................................. 79 Finite element representation of a rotating beam with Q = constant. . . . . 109 Comparison of the static axial blade deflection, us, obtained from an ana- lytical solution [3] and by the finite element method. ............ 110 Convergence of the first four linear natural frequencies (flapping modes) of the finite element model ............................. 111 The first nonlinear flapping natural frequency as a function of initial mode amplitude a(0) with (b = 0 obtaining using three methods: the FE—based ROM, the Rayleigh-Ritz based ROM, and the reference model (182-element and 45-kept—modes). .............................. 112 The projection of the motion occurring on the first nonlinear flapping mode manifold onto the linear eigenspace of the first linear flapping mode, in am- plitude and phase coordinates .......................... 112 The time responses of the first nonlinear flapping mode displacement 171 for the F E—based ROM, the Rayleigh-Ritz based ROM, and the reference model. 113 The contribution of the second linear flapping mode 772 = P2(a, (15) to the first nonlinear mode manifold, solved by the collocation method, as a function of the first nonlinear mode amplitude and phase. ................ 113 The time response of the second linear flapping mode deflection 772 = P2(a(t), ¢(t)) on the first nonlinear mode manifold for the FE—based ROM, the Rayleigh-Ritz based ROM, and the reference model. .......... 114 The contribution of the tenth linear flapping mode 1710 = P10(a, (ii) to the first nonlinear mode manifold, solved by the collocation method, as a function of the first nonlinear mode amplitude and phase. ................ 114 The time response of the tenth linear flapping mode deflection 7710 = P10(a.(t), ¢(t)) on the first nonlinear mode manifold for the F E-based ROM, the Rayleigh-Ritz based ROM, and the reference model. .......... 115 The contribution of the first linear axial mode 1722 = P2201, Q5) to the first nonlinear mode manifold, solved by the collocation method, as a. function of the first nonlinear mode amplitude and phase. ................ 115 xii 3.12 3.13 3.14 3.15 3.16 3.17 3.18 4.1 4.2 4.3 4.4 4.5 The time response of the first linear axial mode deflection 1722 = P22(a(t), ¢(t)) on the first nonlinear mode manifold for the FE—based ROM, the Rayleigh-Ritz based ROM, and the reference model. .......... The contribution of the twelfth linear axial mode n33 = P33(a, (b) to the first nonlinear mode manifold, solved by the collocation method, as a function of the first nonlinear mode amplitude and phase. ................ The time response of the twelfth linear axial mode deflection 7733 = P33(a(t), ¢(t)) on the first nonlinear mode manifold for the F E—based ROM, the Rayleigh-Ritz based ROM, and the reference model. .......... The time dependence of the beam-tip transverse deflection, w(L, t), for the FE—based ROM, the Rayleigh-Ritz based ROM, and the reference model. . Transverse deflection w(:v, t) for a quarter-period of motion on the first non- linear mode manifold. The motion starts at the maximum deflection (the bottom curve) and moves as shown to the zero deflection at a quarter-period, after which it moves upward in a symmetric manner about w = 0 and then repeats. ..................................... Angle deflection 0(1v, t) for a quarter-period of motion on the first nonlinear mode manifold. The motion starts at the maximum deflection (the top curve) and moves as shown to the zero deflection at a quarter-period, after which it moves downward in a symmetric manner about 0 = 0 and then repeats. ..................................... Axial deflection u.(;1:, t) for a quarter-period of motion on the first nonlinear mode manifold. The dashed line denotes the static deflection, u3(:r). The bottom curve corresponds to the initial condition. The top curve from Figure 3.16 and the bottom curve from Figure 3.17 correspond to the top u(:r) curve here, all at one quarter-period. Note that the axial motion occurs at twice the frequency of the transverse and angle motions, due to symmetry. As time progresses beyond the quarter-period, the axial motion moves downward again, and oscillates between the top and bottom curves shown. ...... A five-DOF nonlinear spring-mass system. .................. System substructures. ............................. The contribution of the second fixed-interface linear mode amplitude 779% = Xzfimiv’fi, 7.7{Vfl) to the first fixed—interface nonlinear mode manifold of sub— structure 6. ................................... The contribution of the second fixed-interface linear mode amplitude "$7.7 = X3079”, 779,”) to the first fixed-interface nonlinear mode manifold of sub- structure 7. ................................... The time response of the displacement $1 of mass ml from Figure 4.1, which corresponds to 513?, :13? , and :517 from Figure 4.2, for the nonlinear-CMS model, the linear—CMS nonlinear model, the linear-CMS linear model, and the orig- inal model ..................................... xiii 116 116 117 117 118 119 120 168 168 169 169 170 4.6 4.7 4.8 4.9 4.10 4.11 4.12 4.13 4.14 4.15 4.16 4.17 4.18 4.19 4.20 The time responses of the displacement $2 from Figure 4.1, which corre- sponds to :13? from Figure 4.2, for the nonlinear—CMS model, the linear-CMS nonlinear model, the linear-CMS linear model, and the original model. . . 170 The time responses of the displacement $3 from Figure 4.1, which corre- sponds to I? from Figure 4.2, for the nonlinear-CMS model, the linear-CMS nonlinear model, the linear-CMS linear model, and the original model. . . . 171 The time responses of the displacement :54 from Figure 4.1, which corre- sponds to it); from Figure 4.2, for the nonlinear-CMS model, the linear-CMS nonlinear model, the linear-CMS linear model, and the original model. . . . 171 The time responses of the displacement 11:5 from Figure 4.1, which corre sponds to 173 from Figure 4.2, for the nonlinear—CMS model, the linear-CMS nonlinear model, the linear-CMS linear model, and the original model. . . . 172 A class of nonlinear spring-mass system. ................... 172 System substructures. ............................. 173 The first four fixed-interface LNMs of substructure [3 using the first-trial parameters. ................................... 173 The coefficients of the pqr cubic nonlinear terms expressed in fixed-interface linear modal coordinates and associated with the first mode (master mode) of substructure B using the first-trial parameters. .............. 174 The contribution of the second fixed-interface linear mode amplitude 77? ”[3 = X2fi(n]V,fi, niv’fi) to the first fixed-interface nonlinear mode manifold of sub- structure 5 using the first-trial parameters ................... 174 5: to the first fixed-interface nonlinear mode manifold of sub- The contribution of the third fixed-interface linear mode amplitude 1797’ 3 (771 v 771 ) structure 6 using the first-trial parameters ................... 1 The contribution of the fourth fixed-interface linear mode amplitude 1751,13 = X f (nfv ’fl , nil/fl ) to the first fixed-interface nonlinear mode manifold of sub— "‘1 CI! structure 6 using the first-trial parameters ................... 175 The time responses of the first fixed-interface mode displacement nil/’3 for the ROM and the original model of substructure 6 using the first-trial parameters. 176 The time responses of the second fixed-interface linear mode displacement név’fi = Xgmiv’fi, nil/fl) for the ROM and the original model of substructure B using the first-trial parameters. ....................... 176 The time responses of the third fixed-interface linear mode displacement név’fi = Xgmiv’fl, niv’fi) for the ROM and the original model of substructure 6 using the first-trial parameters. ....................... 177 The time responses of the fourth fixed-interface linear mode displacement nil/fl = X f (77]V ’5 , 7-)]Vfl ) for the ROM and the original model of substructure L3 using the first-trial parameters. ....................... 177 xiv 4.21 4.22 4.23 4.24 4.25 4.26 4.27 4.28 4.29 4.30 4.31 4.32 The time responses of the first fixed-interface mode displacement 779% for the linearized one-mode model, the one-mode model, and the ROM of sub- structure 6 using the first-trial parameters ................... The first four fixed-interface LNMs of substructure 3 using the second-trial parameters. ................................... The coefficients of the pqr cubic nonlinear terms expressed in fixed-interface linear modal coordinates and associated with the first mode (master mode) of substructures ,6 using the second-trial parameters. ............ The contribution of the second fixed-interface linear mode amplitude 779513 = Xgmiv’fi, niv’fi) to the first fixed-interface nonlinear mode manifold of sub- structure 6 using the second-trial parameters. ................ The contribution of the third fixed-interface linear mode amplitude név’fi = X g (7)]V ’3 , niv’fi ) to the first fixed-interface nonlinear mode manifold of sub- structure 6 using the second-trial parameters. ................ The contribution of the fourth fixed-interface linear mode amplitude nil/(3 = X f (nil/fl , rift/6 ) to the first fixed—interface nonlinear mode manifold of sub- structure 6 using the second-trial parameters. ................ The time responses of the first fixed-interface mode displacement 779MB for the linearized one-mode model, the one-mode model, the ROM, and the original model of substructure 3 using the second-trial parameters. The time responses of the second fixed-interface linear mode displacement név’fi = Xgmiv’fi, niv’fi) on the first fixed-interface nonlinear mode manifold for the ROM and the original model of substructure ,8 using the second—trial parameters. ................................... The time responses of the third fixed-interface linear mode displacement N, N, .N, 773 fi=X§0Il fin] 5) for the ROM and the original model of substructure (3 using the second-trial on the first fixed-interface nonlinear mode manifold parameters. ................................... The time responses of the fourth fixed-interface linear mode displacement nf’fl = X f (77?“3 , 7.7]Vfi ) on the first fixed-interface nonlinear mode manifold for the ROM and the original model of substructure B using the second-trial parameters. ................................... The time responses of the displacement 2:? from Figure 4.11 with the inter- face fixed for the linearized one—mode model, the one-mode model, the ROM, and the original model of substructure 6 using the second-trial parameters. The time responses of the displacement mg] from Figure 4.11 with the inter— face fixed for the linearized one—mode model, the one-mode model, the ROM, and the original model of substructure 6 using the second—trial parameters. XV 178 178 179 179 180 180 181 181 182 182 183 183 4.33 4.34 4.35 4.36 4.37 4.38 4.39 4.40 4.41 4.42 4.43 4.44 i3 i fixed-interface linear mode manifold (linear eigenspace) and on the first fixed- Comparison of deflections 2: (t) for a quarter-period of motion on the first interface nonlinear mode manifold of substructure B using the second-trial parameters. The motion starts at the maximum deflection (the bottom curve) and moves as shown to the zero deflection at a quarter-period. Note that linear mode shape is normalized such that r? is the same as 1’13 of the nonlinear mode shape at each time ....................... The first four fixed-interface LNMs of substructure 7 ............. The coefficients of the pqr cubic nonlinear terms expressed in fixed-interface linear modal coordinates and associated with the first mode (master mode) of substructure 7 ................................. The contribution of the second fixed-interface linear mode amplitude 17?,” = X g (nil/’7, 179,”) to the first fixed-interface nonlinear mode manifold of sub- structure 7. ................................... The contribution of the third fixed~interface linear mode amplitude 77:]; {’7 = X g (179,”, 17¢”) to the first fixed-interface nonlinear mode manifold of sub— structure 7. ................................... The contribution of the fourth fixed-interface linear mode amplitude 041‘,” = XZ(77[V’7,1'7[V’7) to the first fixed-interface nonlinear mode manifold of sub- structure 7. ................................... The time responses of the first fixed-interface mode displacement nix,” for the linearized one-mode model, the one-mode model, the ROM, and the original model of substructure 7. ....................... The time responses of the second fixed-interface linear mode displacement 179,” = X 3 (779m, 7'79”) on the first fixed-interface nonlinear mode manifold for the ROM and the original model of substructure 7. ........... The time responses of the third fixed-interface linear mode displacement n?” = X g (7)]V ’7, nil/’7) on the first fixed-interface nonlinear mode manifold for the ROM and the original model of substructure 7. ........... The time responses of the fourth fixed-interface linear mode displacement nil/’7 = X 207?”, fly”) on the first fixed-interface nonlinear mode manifold for the ROM and the original model of substructure 7. ........... The time responses of the displacement £133 from Figure 4.11 with the inter- face fixed ( as? (t) = 0 ) for the linearized one-mode model, the one-mode model, the ROM, and the original model of substructure 7 .......... The time responses of the displacement $31 from Figure 4.11 with the in- terface fixed for the linearized one-mode model, the one-mode model, the ROM, and the original model of substructure 7. ............... xvi 184 184 185 185 186 186 187 187 188 188 189 189 4.45 4.46 4.47 4.48 4.49 4.50 4.51 4.52 4.53 4.55 Comparison of deflections 2:270) for a quarter-period of motion on the first fixed-interface linear mode manifold (linear eigenspace) and on the first fixed- interface nonlinear mode manifold of substructure 7. The motion starts at the maximum deflection (the bottom curve) and moves as shown to the zero deflection at a quarter-period. ......................... The time responses of the first synthesized fixed-interface mode displacement 779”" for the three-DOF linear-CMS linear model, the three-DOF linear-CMS nonlinear model, the three-DOF nonlinear-CMS model, and the forty-one- DOF linear-CMS nonlinear model of substructure B. ............ The time responses of the second synthesized fixed—interface linear mode displacement 779% for the three-DOF nonlinear-CMS model and the forty- one—DOF linear-CMS nonlinear model of substructure fl. .......... The time responses of the third synthesized fixed-interface linear mode dis- placement 17f;v ’fl for the three-DOF nonlinear-CMS model and the forty-one- DOF linear-CMS nonlinear model of substructure fl. ............ The time responses of the fourth synthesized fixed-interface linear mode displacement niv’fi for the three-DOF nonlinear-CMS model and the forty- one-DOF linear-CMS nonlinear model of substructure fl. .......... The time responses of the first synthesized fixed-interface mode displacement nil/37 for the three-DOF linear-CMS linear model, the three-DOF linear-CMS nonlinear model, the three-DOF nonlinear-CMS model, and the forty-one— DOF linear-CMS nonlinear model of substructure 7 .............. The time responses of the second synthesized fixed-interface linear mode displacement 1);] ’7 for the three-DOF nonlinear-CMS model and the forty- one-DOF linear-CMS nonlinear model of substructure 7. .......... The time responses of the third synthesized fixed-interface linear mode dis- placement "97,7 for the three—DOF nonlinear-CMS model and the forty-one- DOF linear-CMS nonlinear model of substructure 7 .............. The time responses of the fourth synthesized fixed-interface linear mode displacement nil/’7 for the three-DOF nonlinear-CMS model and the forty— one-DOF linear-CMS nonlinear model of substructure 7. .......... The time response of the displacement $1 of mass m1 from Figure 4.10, which corresponds to r?,x?, and 2517 from Figure 4.11, for the three-DOF linear-CMS linear model, the three—DOF linear-CMS nonlinear model, the three-DOF nonlinear-CMS model, and the forty-one—DOF original model. . The time response of the displacement 51:2 of mass m2 from Figure 4.10, which corresponds to as? from Figure 4.11, for the three-DOF linear-CMS linear model, the three-DOF linear-CMS nonlinear model, the three-DOF nonlinear-CMS model, and the forty-one-DOF original model. ....... xvii 190 190 191 191 192 192 193 193 194 194 4.56 4.57 4.58 4.59 4.60 The time response of the displacement 1:21 of mass mm from Figure 4.10, which corresponds to 22% from Figure 4.11, for the three-DOF linear-CMS linear model, the three-DOF linear-CMS nonlinear model, the three-DOF nonlinear-CMS model, and the forty-one—DOF original model. ....... The time response of the displacement 1:22 of mass mm from Figure 4.10, which corresponds to 2:; from Figure 4.11, for the three-DOF linear-CMS linear model, the three-DOF linear-CMS nonlinear model, the three-DOF nonlinear-CMS model, and the forty-one—DOF original model. ....... The time response of the displacement $41 of mass m41 from Figure 4.10, which corresponds to 2:31 from Figure 4.11, for the three-DOF linear-CMS linear model, the three—DOF linear-CMS nonlinear model, the three-DOF nonlinear-CMS model, and the forty-one—DOF original model. ....... The time responses of the first synthesized fixed-interface mode displacement 179w for the three-DOF nonlinear-CMS model and the forty-one—DOF linear- CMS nonlinear model of substructure 6 with ”955(0) 2 1.3, "{qu) = 1.3, nC(0) = 0.04,Nfl(0) = 0, fifWO) = 0,1')C(0) = 0. .......... The time responses of the first synthesized fixed-interface mode displacement nil/’7 for the three-DOF nonlinear-CMS model and the forty-one-DOF linear- CMS nonlinear model of substructure 7 with nil/(3(0) = 1.3, n[v’7(0) = 1.3, 770(0) = 0.04,[V’5(0) = 0, 77]“(0) = 0. r;C(0) = 0. .......... xviii 195 196 196 197 197 CHAPTER 1 Introduction 1.1 Motivation h‘lodern engineers are sometimes required to accurately model the motions of complex structural systems, i.e., rotorcraft, turbo machines (turbines), vehicles, aircraft, bridges, off-shore platforms, robotic arms, etc., so that they can function properly and safely under certain operating conditions. Nowadays such structural systems are designed to be lighter, have more complex geometry, and operate at high speeds in order to enhance performance. These conditions can cause structural systems to experience responses in the nonlinear regime. To mathematically describe such motions, one must formulate nonlinear models, either discrete or continuous in nature. If they are continuous, in practice engineers typ- cially discretize them by using either direct modal discretization (Rayleigh-Ritz method) or the finite element (FE) method. The FE method is typically prefered in practice since it is a versatile tool for modeling very complex structures. However, it has a disadvantage in that it often requires a very large number of degrees of freedom (DOF) to accurately model the structure ([4]). The dynamic analysis (e.g., determination of natural frequen- cies, simulations, etc.) of nonlinear structures is typically more easily performed in modal coordinates than in physical coordinates, since it requires a smaller number of DOF. How- ever, many modal coordinates are still needed for some nonlinear structures, typically those with some type of essential coupling between the linear modes. Therefore, this can result in considerable computational effort when ones try to analyze their dynamics ([5]). The main objective of this research is to develop a class of techniques for generating ac- curate reduced order models for nonlinear structural systems. Our starting point is gener- ally a nonlinear model, from either modal or FE analysis, whose dynamics are expressed by equations of motion expressed in terms of the linear modal coordinates. First, this work demonstrates a general framework on how to obtain nonlinear reduced order mod- els (ROMS) using nonlinear normal modes (NNMS) that are defined in terms of invariant manifolds. Next, a methodology for obtaining nonlinear ROMS from FE based descriptions of structural elements is demonstrated through an example, the nonlinear rotating beam, which is a crude model for a helicopter rotor blade. Lastly, this work demonstrates how one can synthesize nonlinear ROMS using NNMs of relatively simple substructures that are assembled to form a complex structural system. In this first chapter we offer brief overviews and literature surveys of the main themes considered in the thesis, Specifically: nonlinear normal modes, rotating beams, and component mode synthesis. These are considered sequentially, and are followed by a description of the outline of the remainder of the thesis. 1.2 Nonlinear Normal Modes Linear modal analysis is a fundamental tool used in linear vibration theory. It allows one to decompose a general motion into a linear combination of the fundamental motions (modal motions) that take place on the eigenspaces in the system phase space; this is the essence of superposition. These modal responses are invariant, since, if one starts with a purely modal response, the system stays in that mode for all time. Similarly, if the system is started with energy in a subset of modes, only those modes will be active during the ensuing response. For general responses, the individual modal responses remain uncoupled from one another during the motion. Of course, these facts do not generally hold for nonlinear systems. Certainly superposition does not extend to nonlinear systems. But, the concept of fundamental invariant responses does carry over, even if one cannot construct general responses by combining them, whether in a linear or nonlinear manner. Rosenberg and Atkinson ([6] and [7]) provided the pio- neering ideas which tried to extend the idea of fundamental motions to nonlinear systems. Rosenberg and coworkers carried out several studies along these lines ([6], [8], [7], [9], [10], [11]) and summarized these results in a classical review paper [12], in which Rosenberg offered a break-through definition of nonlinear normal modes (N NMS). In that work, the NNMS of conservative systems with 71 degrees of freedom were defined as vibration in uni- son, i.e., vibrations such that all degrees of freedom reach their extrema at the same time and pass through zero at the same time. By this definition, one can express all generalized displacements as functions of a chosen generalized displacement. When these systems pos- sess certain types of symmetries, the constraint relations are linear in the configurations space (just as they always are in the linear case) and are referred to as “similar” NNMS. For more general, but still conservative, systems the constraint relations are nonlinear and the nonlinear modes are called “non-similar” N N Ms, which can be depicted by curved lines in the configuration Space. This definition was picked up and used by subsequent researchers who used it to construct NNMS and study their stability and bifurcations that occur due to changes in system parameters and system energy levels. These were typically conservative systems with two degrees of freedom (DOF) system ([13], [14], [15]). This definition, to- gether with the conservation of energy, were employed by Vakakis to construct non-similar NNMS of a two DOF-conservative system in [16]. It was also used to study the steady state motions of two DOF systems subjected to periodic forcing; see [17] for Similar N NMS and [18] for nonsimilar NNMS. Vakakis also discovered that the important phenomenon of mode localization can occur for NNMS, even for perfectly tuned subsystems [16]. This is in contrast to linear systems, wherein localized modes exist only when the subsystems are Slightly mistuned. A good account of NNMS, their stability and bifurcations in unforced nonlinear systems, their existence in forced nonlinear systems, and their application to engineering systems is given in the monograph by Vakakis [19]. The N NM concept has been generalized to a wide class of systems, which includes discrete systems with dissipation and gyroscopic terms, as well as quite general continuous systems by Shaw and Pierre in ([20] and [21]). Therein, a definition of NNMS is given in terms of invariant manifolds, which are a natural way to define and construct fundamental motions of nonlinear systems. The procedures used to obtain the NNMS by this approach are closely related to center manifold theory, which is used for bifurcation analysis ([22]), and inertial manifold theory, which is used to study the long time behavior of dissipative partial diflereiitial equations (PDES) ([23], [24]). Using this definition, the nonlinear system equations are restricted to a two-dimensional invariant manifold that describes the NN M of interest. The behavior on the N N M manifold is governed by a single second order differential equation of motion, which corresponds to the equation of motion for a nonlinear single mode reduced order model. In Shaw and Pierre ([20], [21]) approximate solutions of the invariant manifold equation were obtained by asymptotic expansions using polynomials expressed in terms of a generalized position-velocity pair of state variables. King and Vakakis ([25]) developed an energy—based N N M approach based on [21] to investi— gate N N Ms of one dimensional, conservative, continuous systems. N ayfeh and N ayfeh ([26]) computed NNMS of continuous systems based on a complex amplitude / phase formulation and the method of multiple scales. Nayfeh ([27]) compared the various methods for con- structing N NMS of continuous systems as developed by Shaw and Pierre, King and Vakakis, Nayfeh and Nayfeh, and a new approach that employed normal form theory. They con- cluded that all expansion methods yielded the same results, but claimed that the method of multiple scales with complex variables was the simplest to implement. All of these approaches are in some manner equivalent, but differ in terms of formulation, solution, and range of applicability. However. the invariant manifold definition is the most general, since it covers the widest range of systems and responses. It can also be generalized to the case of multiple modes, as required for the case of internal resonances, or if one fecmires a model that is valid over a wide frequency range. By defining the NNMS in terms of two-dimensional invariant manifolds, the individual NNMS can be constructed, but they can not interact with each other once a motion is initiated on any one of them. Therefore, a motion involving multiple modes cannot be captured by this definition. In addition, the concept of mode superposition cannot be applied to construct non-linear multi-mode models using individual NNMS, since the essential interaction between the NNMS will be missing. Boivin et al. ([28]) generalized the individual N NM concept by defining a motion involving M N NMS as a motion that takes place on a 2M -dimensional invariant manifold in the system’s phase space. By this definition, the non-linear ODES are restricted to the M-NNMS invariant manifold of interest. The behavior on the manifold is governed by M second order, coupled equations of motion (a nonlinear M -degree-of-freedom system). Boivin et al. also described how to detect the case of internal resonances from the multi- mode invariant manifold formulation [29]. N NMS were also constructed for discrete systems with internal resonances, based on the complex formulation, by Nayfeh et al. [30]. King et al. ([31]) also extended the energy-based NN M approach to cover the case of internal resonances. These latter two studies were for the case of M = 2 N N Ms. Slater ([32]) developed a numerical method for determining individual NNMS based on numerical searching techniques for periodic solutions of conservative non-linear systems. Using this approach, individual solutions are found, but the entire family of motions on the NNM manifold cannot be obtained, therefore one cannot develop ROMS using this ap- proach. A similar shooting technique has been used in this work and by previous resarchers in the author’s research group ([3], [33], [34], [35] and [36]). In this thesis it is employed for comparing ROMS with simulations of the original system, which has the full number of DOF. Other works, related to the use of invariant manifolds for the generation of ROMS, include the use of Karhunen-Loeve (K-L) decomposition to develop accurate low-order models, by using data obtained from transient simulations of large-scale systems. See [37], [38], [39], and [40], for example. In references [20], [21], [28], [41], and [29], the N NM invariant manifolds were approximated by a polynomial expansion (asymptotic series) of position-velocity pairs of chosen (master) modes, which provides a solution that is locally valid. The approximate invariant manifolds obtained by such an asymptotic series will diverge from the actual invariant manifold in some amplitude regime. Typically, the domain of validity of the approximation is not known a priori. Hence, the reduced equations of motion will generate inaccurate time responses when the amplitudes of the modes are “too large”. Pesheck et al. ([33]) improved the approximation of the invariant manifolds by expressing the invariant manifolds as an expansion of basis functions defined over a specified domain and numerically solving the invariant manifold equations. The expansion of basis functions is introduced into the PDEs governing the NNM invariant manifold, and, using a Galerkin projection, the nonlinear equations for the expansion coefficients were obtained and then solved numerically. By this approach, the domain of approximation can be selected by the user, and the error of approximation can be minimized over the chosen domain by selection of the number and type of basis functions. Since the computational cost associated with the Galerkin projection can be quite expensive, the collocation method has recently been adopted in [35] to minimize computational efforts associated with this method. In this approach, instead of projecting each manifold governing equation onto each basis function, each manifold governing equation is projected onto a set of Dirac delta functions in the master coordinates, thereby providing a solution that minimizes an error that is measured in a point-wise manner. Research on the construction and use of NNMS continues; here we outline recent work by others in the MSU/UM NNM research group, in particular by Mr. Dongying Jiang, a Ph.D. student at UM. Jiang et al. ([36]) have applied the work of Pesheck et al. ([33]) to construct NNMS of piecewise linear systems, and are working on systems with friction elements. They have also extended the work of Pesheck et al. ( [33]) to numerically construct multi-NNM models which can capture internal resonances among participating modes, and are valid over a large range of amplitudes ([42]). They have also numerically constructed NNMS of nonlinear systems under periodic excitation ([43]). In this case the manifolds are time-periodic in nature, and responses on them represent the steady-state responses of the full system. Current work is aimed at distilling ROMS from detailed FE models for rotor blades. A detailed summary of the invariant manifold approach to NNMS, and the computational issues associated with its solution, are presented in Chapter 2 of this thesis. 1.3 Beam (Blade) Dynamics The dynamic analyses of helicopter blades, turbopropeller blades, wind-turbine blades and robotic arms has provided motivation for investigations of the vibration of rotating beams. To predict the dynamic characteristics of rotating flexible structures, the kinematics must be carefully modeled, which leads to nonlinear coupling effects between degrees of freedom (DOF) in different directions. These coupling effects can cause slow modal convergence, thereby often requiring large system models for accurate dynamic representation. Simula- tion of such large—scale models is a time consuming process, which slows parametric studies and design cycles. Much work has been done using finite elements (FE) to model the nonlinear, large amplitude vibrations of rotating beams, including [44], [45], [46], [47], [48], and [49]. These models are typically complex in nature due to their geometry, degrees of freedom (flap, lead-lag, axial, and torsion), and nonlinear coupling effects. Furthermore, because of the nature of the finite element approach, many elements are required in order to obtain an accurate model. A common approach is to use linearization of the finite element model about the nonlinear static equilibrium position and solve the eigenvalue problem of the resulting linearized model to obtain the linear natural frequencies of the system ([44] and [49]). Bauchau and Hong ([45]) also utilized finite elements in time to obtain nonlinear responses and stability results of the rotating beam undergoing large deflections. However, the computational time associated with obtaining the equilibrium solution was expensive, because all of the spatial degrees of freedom are coupled at all time steps. Perturbation modes ([50] and [51]) were applied to the finite element model of a helicopter rotor blade in order to obtain a reduced order model ([46]). Bauchau and Bottasso ([52]) applied the perturbation modes to the Space-time finite element model of a beam subjected to a sinusoidal load in order to obtain a reduced order model. Crespo da Silva ([53]) utilized a truncated set of eigenfunctions or eigenvectors obtained from the linearized system of PDES or the linearized finite element model about the nonlinear static equilibrium position in order to obtain a reduced order model of a beam in planar motion. Crespo da Silva ([54]) also extended his work to handle multi-beam structures in planar motion. In most nonlinear structures, there is no simple expansion of basis vectors which decouples the DOF (i.e., modes) in the frequency range of interest from those outside that range. For the rotating-beam problem, this is evident in the nonlinear coupling between transverse and axial motions. Therefore, some (potentially important) nonlinear effects may be ignored in the truncation process, unless one is careful. Generally, many linear modes must be retained in the nonlinear model in order to minimize these effects. NNMS is a natural approach for handling this issue. Over the past decade, systematic procedures have been developed to obtain ROMS via NNMS that are based on invariant manifolds in the state space of nonlinear systems, as described above. These procedures initially used asymptotic series to approximate the geometry of the invariant manifold and have been used to study the nonlinear rotating Euler-Bernoulli Beam ([3]). More recent work has employed a numerically-based Galerkin approach to obtain the geometry of the NNM invariant manifolds out to large amplitudes ([33]). These procedures can be applied to more general nonlinearities over wider amplitude ranges, and have been recently applied to study the vibrations of a rotating Euler-Bernoulli beam ([34]). These approaches have provided accurate models for the fundamental nonlinear flapping mode, by systematically capturing the essential dynamic coupling that exists between the linear modes of the system. Chapter 3 of this thesis considers this problem in detail, by investigating the N NMS of a rotating beam that is modeled using a nonlinear finite element formulation. 1.4 Substructure Synthesis Many complex structures are composed of several relatively simple substructures that are assembled together. This occurs in trusses, bladed disk assemblies in turbine rotors, aerospace and ground vehicles, and other applications. In such cases it is convenient to develop a dynamic model for the overall structure by taking advantage of the dynamic properties of the substructures. Methods for doing this for linear structural models are well developed and have been used extensively, especially in the aerospace industry ([55], [56], [57], [58], etc). These techniques construct ROMS of the overall structure by making use of modal-based ROM descriptions of the substructures and combining these using a technique known as Component Mode Synthesis (CMS). In this section we offer an overview of CMS for linear systems and describe the two main CMS approaches, fixed- and free- interface CMS. A thorough review of substructuring and CMS can be found in [59]. 1.4.1 Component Mode Synthesis CMS was developed to synthesize models that are described in terms of substructures, and to take advantage of model size reduction carried out at the substructure level ([60], [61], [62]). In CMS, the dynamics of each substructure is described by a set of dynamic (normal) modes and a set of static (constraint or residual attachment) modes that are used to describe the interfaces between the substructures. A set of component normal modes is selected from each substructure and are chosen and truncated in such a way that the modes lie-in the frequency range of interest. A set of static (so-called constraint) modes is a key component for the low frequency response of the structure ([63]) and are used to couple the substructures together. There are two general types of CMS methods; they are known as the fixed-interface and the free-interface approaches, as briefly described below. 1.4.1.1 Fixed—Interface Linear CMS The fixed-interface CMS technique, developed by [62], is widely used, since the procedures are straightforward. Moreover, it produces very accurate models with very few component modes ([64]). The dynamics of each substructure is described by its modal description, which is composed of substructure normal modes (component modes) and constraint modes. The component normal modes are the normal modes of the substructure derived for the case when the interface coordinates between the substructures are held fixed. A constraint mode (static mode) of the substructure is the deflection obtained by imposing a unit displacement on one of the interface coordinates and holding the remaining interface coordinates fixed. To obtain all of the constraint modes, the process is applied in turn to each of the interface coordinates. By applying displacement and force compatibility conditions at the interface coordinates, one can obtain the reduced synthesized system, which is described by the component normal-mode coordinates and the generalized constrained coordinates. This approach iS known as the Craig-Bampton method. 1.4.1.2 Free-Interface Linear CMS Ree-interface CMS methods are more attractive than fixed-interface CMS methods when the component modes are obtained from modal testing or when an experimental verification of the component modes is required ([65]). The free-interface CMS technique developed by Craig and Chang ([66], [67] and [63]) is the most accurate among the free-interface CMS techniques. It is a modified version of Rubin’s method [68] and MacNeal’s method [69]. It is superior to the CMS of Craig-Bampton in terms of accuracy, but is more difficult to implement ([64] ). 10 In this synthesis technique, the dynamics of each substructure is described by a set of substructure component normal modes and residual attachment modes. The component normal modes are the normal modes of the substructures for the case when the interface coordinates between the substructures are free. The residual attachment modes of the substructure are a Special type of static modes that are used to couple the substructures together, and they also account (at least partially) for the static deflections of the truncated normal modes of the substructures ([66], [67] and [63]). By applying the displacement and force compatibility conditions at the interface coordinates, and neglecting the inertial effects associated with the generalized residual attachment coordinates, one can obtain the reduced synthesized system. It is described in terms of only the component normal mode coordinates since, in this approach, the residual attachment modes can be condensed out of the equations of motion ([66], [67] and [63]). This approach is known as the Craig-Chang method. Chapter 4 of this dissertation describes a CMS methodology developed for nonlinear systems, wherein the substructure ROMS are developed using NNMS, and these are assembled using a newly developed technique that is a extension of the fixed-interface CMS method. 1.5 Dissertation Organization The Dissertation is organized as follows. In Chapter 2, the formulation and solution of NN M invariant manifold equations are de- scribed. Both previous and original work is described. Asymptotic series expansions and weighted-residual type methods, i.e., Galerkin and collocation, are the methods employed to solve for the manifolds. In this chapter, two new methods, similar to those described in [33], are developed and implemented to solve for the single-mode manifold solutions. They 11 are both formulated in terms of modal position and velocity. The first alternative is to solve for the manifold solution using globally defined basis functions using the Galerkin method. The second is to obtain the manifold solution over several small patches using locally de- fined basis functions, and solving these smaller problems using a collocation method. A rough comparison of the relative computational efforts of these approaches is provided. Also, a two-DOF nonlinear spring-mass system and a FE model of a rotating beam are used for demonstrating these approaches and comparing their accuracies, by direct time Simulations of the corresponding ROMS and the original system models (with full DOF). In Chapter 3, the generation of nonlinear ROMS from FE based descriptions is demon- strated through the application to a nonlinear rotating beam, which is a highly idealized model for a helicopter rotor blade. First, the nonlinear FE model of a rotating beam is gen- erated and converted into a truncated (but still large—scale) modal form that is convenient for the NNM analysis. The invariant manifold equations are formulated, and a numerical collocation method is used to obtain the solution of the NN M invariant manifold for the fundamental flapping mode of the beam. This invariant manifold is used to construct a nonlinear single DOF ROM, which is subsequently used for a Simulation study. Note that the results of Chapter 3 have been recently published [35]. In Chapter 4, the development of N N Ms, as needed for the individual substructures is first reviewed. The associated invariant manifold equations, parameterized by modal position and velocity, are formulated, and the numerical collocation method is reviewed, since it allows one to obtain the solution of the NNM invariant manifold. Then, the procedures for the fixed-interface nonlinear CMS are described. A five—DOF spring-mass system and a forty-one-DOF spring-mass system are used to demonstrate the effectiveness of the method, via direct time—simulation comparisons of different ROMS. Chapter 5 offers a summary of the contributions of this research to the field of nonlinear structural dynamics, along with a discussion of potential areas for further work in this area. 12 CHAPTER 2 Comparison of Methods for Constructing Invariant Manifolds for Nonlinear Normal Modes 2.1 Introduction In linear vibration theory one can decompose a general free vibration response as a linear combination of the fundamental motions (modal motions) that take place on the eigenspaces of the system. These motions are invariant, that is, if one initiates a purely modal response, the system remains in that mode for all time. The works of Rosenberg and Atkinson ([6] and [7]) were the pioneering works which extended the idea of fundamental motions to nonlinear systems. In those works, nonlinear spring-mass systems with two degrees of freedom (DOF) were studied. Due to special symmetries possessed by the systems, they were able to construct linear modal-constraint relations, the attendant nonlinear natural frequencies, and the stability of each nonlinear mode in terms of its amplitude. Rosenberg summarized his works ( [6], [8], [7], [9], [10], [11]) on normal mode vibrations of nonlinear systems in [12], in which a fundamental definition of nonlinear normal modes (NNMS) for multi-DOF systems was given. In that work the concepts of similar and non-similar modes were introduced. Caughey et al. ([14]) studied similar normal modes and their bifurcations in terms of system parameters for two-DOF nonlinear conservative system, 13 based on Atkinson’s work ([7]). Rand et al. ([15]) studied the general class of two-DOF nonlinear conservative systems. In that study, they were able to determine the number and stability of periodic motions of which the non-similar modes were special cases and also studied the bifurcation of the periodic motions in terms of system parameters. Vakakis ([16]) studied the non-similar modes of two-DOF nonlinear conservative systems, where the construction of the non-similar modes exploited the conservation of energy. Shaw and Pierre ([70] and [20]) generalized the definition of NNMS by using concepts from invariant manifold theory, which encompassed a wide class of quite general nonlinear systems. The technique was constructive and allowed for traveling wave (that is, complex) and damped NNMS. Boivin and co-workers et al. ([28], [29], [71]) extended the construction of the Single—mode NNM manifold developed by Shaw and Pierre to the multi—mode case, which allowed for the dynamic interactions among a subset of nonlinear modes of interest. A summary of these works by Shaw, Pierre, and coworkers ([72], [73], [41], [28], [29], and [71]) can be found in [74]. In these papers, the solution of the invariant manifold equatiosn were sought in terms of locally valid asymptotic power series expansions. Pesheck et al. ([33]) developed a numerical method for solving the invariant manifold equations, such that the solution, while not analytically available, was valid over a relatively large amplitude range. An added feature of this approach is that the nonlinear terms are not limited to being smooth functions, which allows for extensions to more general classes of problems ([36]). Even though the Galerkin method gives accurate solutions, the computational cost of acquiring solutions is quite expensive, especially when one attempts to generate multi-mode models ([75]). This motivated the application of the collocation method, which significantly reduces computational costs, but does not sacrifice much in the way of accuracy ( [76], [77], and [78]). This method has been applied in [35] to study the dynamics of a finite-element model of a rotating beam. In the forthcoming work of Jiang et al. ([42]), the Galerkin method developed by Pesheck et al. ([33]) has been generalized to construct accurate multi-mode manifolds that are able to capture internal resonances, and this method is applied to a model of a rotating beam. 14 In this chapter, we develop two alternative methods for the numerical solution of the single-mode invariant manifold equations, and offer a comparison of the various methods that have been used in this and previous studies. The new methods are Galerkin-based, in the same spirit of that developed in [33], except that the equations are formulated in terms of modal position and velocity, instead of modal amplitude and phase. The first new method uses global basis functions and obtains the unknown coefficients for the manifold solution using the Galerkin method; this method is described in section 2.3.2.2. The second new method uses a local patchwork of basis functions, and the unknown coefficients for the manifold solution are obtained using the collocation method; this method is described in section 2.3.2.4. In terms of computational time, both new methods have advantages over the approaches used in [33] and [35]; this is summarized in section 2.3.3. This chapter is outlined as follows. The formulations for the NNM invariant manifolds are first introduced. The power series expansions methods and the weighted-residual type methods, i.e., Galerkin and collocation methods, are described in detail, for both amplitude—phase and displacement-velocity formulations. The computational effort associated with each approach is then discussed, and the accuracy of each approach is considered through direct comparisons in general terms and for calculations carried out for two example problems. Some conclusions are drawn at the end of the chapter. 2.2 Formulations of NNM Invariant Manifold We begin with a general discrete representation of the vibrations of a nonlinear structural system obtained either by a finite element model followed by linear modal expansion, or by a Rayleigh-Ritz approach. We assume that the system at linear order is undamped. (This assumption greatly simplifies the problem, but can be relaxed, in principle.) In this case, the equations of motion for a Q-DOF system are uncoupled at linear order and can 15 be expressed in the form: Ifi+An=f(n,ii) (2.1) where I is the identity matrix, A the diagonal matrix of squared linear natural frequencies, f (17, 1'7) a vector of nonlinear forces, 77 the modal position vector, and 7'] the modal velocity vector. The component form of equation (2.1) is given by 777i+wi2 7h = fiMjflij) (2-2) for i,j=1,2,3,...,Q where w,- is the linear natural frequency of mode 2' and Q is the number of retained linear modes. 2.2.1 Invariant Manifold Equations in Terms of Modal Position and Velocity The fundamental concepts for nonlinear normal modes of discrete, conservative, nonlinear systems were laid out by Rosenberg ([12]). Shaw and Pierre ([70], [20], [21]) used the theory of invariant manifolds for dynamical systems to generalize the concept of a nonlinear normal mode to a wide class of systems, including continuous systems and systems with dissipation and gyroscopic terms. Boivin et al. ([28], [29], [71]) extended the single-nonlinear-mode construction developed by Shaw and Pierre to the multi-mode 16 case, which is able to capture the dynamic interactions among a set of nonlinear modes of interest. A summary on the initial series of works by Shaw, Pierre, and coworkers ([72], [73], [41], [28], [29], and [71]) can be found in [74]. Below we provide a short account of the geometric definition of NNM invariant manifolds and methods of constructing single- and multi-mode versions of N N M models. 2.2.1.1 Singe-Mode Manifold Formulation For present purposes one can consider an invariant manifold to be a low dimensional surface living in the system state space, such that motions initiated on the surface will remain on it for all time. This concept is the key to obtaining a reduced order NNM model, that is, a lower-dimensional dynamical system, Specifically by restricting the equations of motion to this surface. In an N-degree-of-freedom, undamped, nongyroscopic, linear vibratory system, individual modal motions are synchronous responses that take place on a two- dimensional plane (linear eigenspace) in the 2N -dimensional system state space. In this case the response of the system is a time-harmonic standing wave in which all degrees of freedom reach their extrema and pass through zero simultaneously. Such motions are governed by two first-order linear differential equations or, equivalently, by one second-order linear differential equation, which in this case is a Simple undamped oscillator. For linear systems with gyroscopic and/or general dissipative terms this concept is also applicable, however the motions taking place on the invariant planes are generally nonsynchronous, and represent traveling wave responses of the system. This concept can be extended to nonlinear systems as well, however, the motions generally take place on non-planar surfaces, and the motions are governed by two first-order nonlinear differential equations, or by one second—order nonlinear differential equation. For smooth nonlinear systems, the manifolds are curved surfaces that are tangent to the manifolds (eigenplanes) of the linearized systems. In [74], a definition for NNMS based on invariant manifolds was offered, as follows: A normal mode for a nonlinear system is a motion that takes place on a two- 17 dimensional invariant manifold in the system’s phase space. This manifold passes through the stable equilibrium point of interest and, at that point, is tan- gent to a two-dimensional eigenspace of the system linearized about that equi- librium. On this manifold, the system dynamics are governed by an equation of motion involving a pair of state variables; that is, it behaves like a single- degr‘ee-of-freedom system. Based on the above definition, the construction of the nonlinear-mode manifold is straight- forward. In order to search for a particular individual NNM, it is assumed that the NNM manifold is parameterized by a single modal position-velocity pair corresponding to the mode of interest, referred to as the master mode. This is accomplished by using the fact that for a NN M response all of the remaining modal positions and velocities are slaved (constrained) to this master mode. For the kth nonlinear mode, we take uk 2 71k, and vk = 71k as the master states. The remaining slave states are expressed as 772' = Xi('“ks'vk) = Xif’lkflik) 7h = Witt-Wk.) = 33(77ka'7ik) for i = 1,2,3,...,Q, i 74 k. Equations (2.3) and (2.4) constitute a set of constraint equations that are to be determined. The constraint functions in equations (2.3) and (2.4) are obtained by an invariant manifold procedure that generates equations that can be solved for the unknown constraint relations. The process begins by taking a time derivative of the constraint equations, yielding 18 . 8X, , 0X, , - = ' . — i. 2.5 "I auk uh + aUk 1A ( ) .. (913 . BY . r), — 0a]. uk + —ka vk (2.6) for i:1,2,3,...,Q, 2'75 k. The time dependence in these equations is eliminated by using the following relations: ti), = vk, v'k = 7)), = —w,%, 11k+fk(nj,7'7j), and 77',- = ‘wi2 n,+f,-(r)j,7'7j). Then, the constraints (equations (2.3) and (2.4)) are substituted in the resulting expression everywhere in place of the Slave state variables, resulting in a set of partial differential equations for the functions (X,(uk, vk), Y,(uk, vk)). This set of 2Q - 2, time-independent, partial differential equations govern the geometry of the kth manifold, and are given by 8X,- 8X, , 2 Buk ,le + 81—ik (—wk uk +fk(X',Yj, ltk,'Uk)) (2.7) 81" BY' w? X.- + Miner... vi.) = 57212;. + 5—7 (we? ”is +fk(Xj,YJ-,uk,vk)) (2.8) for i,j:1,2,3,...,Q, i,j;£ k. These equations are not solvable in closed form (except in very Special cases) and meth- ods for obtaining approximate solutions for X .1- and Y, are described in detail in this chapter. 19 2.2.1.2 Multi-Mode Manifold Formulation Single-mode motions are not general motions, since they take place on two dimensional manifolds in the system state space. In structural dynamics, we are also interested in more general motions in which more than a single mode participates. For linear systems, multi—mode motions take place in the Space spanned by a set of two-dimensional planes that are the linear eigenspaces of the modes of interest. Furthermore, using superposition, multi-mode motions of linear systems can be obtained through a linear combination of responses of the individual modes of interest. For nonlinear systems, multi-mode motions cannot be obtained through the superposition of individual NNM responses, since the coupling between the NNMS cannot be accounted for in this manner. However, the idea of a nonlinear single-mode manifold can be generalized by defining a nonlinear multi-mode manifold to be a multi- dimensional surface that is tangent to the multi-dimensional linear eigenspace of the corresponding linear modes of interest. Based on these ideas, in [74] the following definition of a nonlinear-multi mode was offered: A nonlinear ll/I-mode invariant motion of a system is a response that takes place on a 2M -dimensional invariant manifold in the system’s phase space; the manifold passes through the stable equilibrium point of interest, and at that point it is tangent to a 2M -dimensi0nal eigenspace of the system linearized about that equilibrium (representing M linear modes). On this manifold, the system dynamics are governed by ll”! pairs of state variables; that is, it behaves like an M’ -degree-0f— freedom system. The procedure for constructing the multi-mode manifold is Similar to the procedure given in section 2.2.1.1, except that the manifolds are now parameterized by 2M variables. Note that the Single-mode manifold is a special case of the multi-mode manifold with M = 1. Specifically, we take uk 2 71k: and vk 2 77),, for h 6 SM. where SM is a set of indices that describes the modes of interest (master modes). For example, if one is interested in 20 constructing a model that captures the interactions between the second, third and fifth nonlinear modes, then SM = {2,3,5}. The sets {uk} and {vk}, k 6 SM, are denoted by u M and vM respectively. The remaining slave states are expressed as m = Xi(UM, vM) (2-9) le'umwm (2.10) H 7h for i=1,2,3,...,Q, re SM. The procedure described in section 2.2.1.1 is again applied here by taking a time derivative of the constraint equations, eliminating the time dependence by using the default rela- tions, the equations of motion of the master modes, and the equations of motion of the slave modes, and substituting the constraints every where in place of the Slave state vari- ables, yielding the following set of 2Q — 2M partial differential equations for the constraint functions: Yi = Z [351(ka + % (_w% “k keSM +fk(Xj»}/ja"-Ika ”kl” (2.11) “W? X1 + fr(XjalG~‘llkvl"k~) I Z [3:ka + 37):: (“9f “k keSM , +fk(Xj1Yj»“k~1/’l.~.))] (2.12) for i,j=1,2,3,...,Q, i,j ¢ SA]. AS in section 2.2.1.1, closed-form solutions do not generally exist. The method for 21 obtaining approximate solutions for X,- and Y,- in terms of power series expansions is reviewed in section 2.3.1.2. 2.2.2 Invariant Manifold Equations in Terms of Modal Ampli- tude and Phase Pesheck et al. ([33]) developed an alternative form of the single-mode invariant manifold equations that are expressed in terms of modal amplitude and phase. The procedure was similar to the formulation in modal position and velocity, and yields results that are, of course, equivalent. Recently Jiang et al. ([42]) have extended this formulation to the multi-mode case. These formulations are described in the following sections. 2.2.2.1 Singe-Mode Manifold Formulation The coordinate transformations used here are typical of those used in the method of averag- ing ([79]), and are similar to those used in the method of variation of parameters, since they employ a time-varying amplitude and phase. They relate the master modal displacement and velocity (nk, 77k) to the master modal amplitude and phase (ak, (bk) via the following invertible transformation: 77k : ak cos((.bk) (2.13) ifk = —ak wk 8111(Ok). (2.14) Hence, the constraints for the slave modes, equations (2.3) and (2.4), can be expressed in term of ak and 3,, as 22 771' = PAM-.951.) = XiUIt-(a-k,$k)fl7'i.-(ak,¢k)) (2-15) 7h = Qi(aka¢kl = K(72k(ak»¢k)ifik(ak,¢kll (2.16) for i: 1,2,3,...,Q, 1'72 k. Using equations (2.13) and (2.14), the equation of motion governing the master mode can be expressed as two coupled 1st order ODES in ak and risk, as follows, — .P', -, , . - (ii.- 2 f1.( 1le :1: <31.) 81110151.), (2.17) k .- .P', ',a.,ct). coscb. @k : wk_f1(JQ] i. A) (A), (2.18) “k. wk for j=1,2,3,...,Q, jaé k. In order to determine the slave constraints, the procedure is identical to that for the (uk, vk) formulation, adapted to these coordinates. Using steps similar to those in equations (2.5) and (2.6), we obtain , 3P,- , 0P,- ,- 4 ,- = — .. — . 2.19 7)! aak (1k + 00k Qbh ( ) _ 8Qi . 3Qz‘ ' r), — Oak "k + 00k¢k (2.20) for i=1,2,3,...,Q, 1:75 k. 23 The time dependence in these equations is eliminated by using equations (2.17), (2.18), and r); 2: wot-2 771+ fi(7lj, 17]). Then, the constraints (equations (2.15) and (2.16)) are substituted in the resulting expression everywhere in place of the slave state variables, resulting in a set of partial differential equations for the functions (Pi(ak, (pk), Q,(ak, pk». This set of 2Q - 2, time-independent, partial differential equations govern the geometry of the kth manifold, and are given by 3Pi(-fk(Pinj~ak~C/bkl Sinf‘r’lkl) Qi = 8 (1k wk 8P, fk(Pj9Qj‘a'k3¢k) COS +86% (wk - ak wk ) (221) , (9 ‘ —fk(P*Qaa“¢) 5111(0) _w’12Pi+fi(Pjanaakv@k) = a—Q—I( J J k A A ) ak wk ' .P'. ',,..‘..‘I., argue _ f1.( 3 Q3111. €151) (08051)) (2.22) 00h. Gk Wk for i,j:1,2,3....,Q, i,j;£ k. These equations are more complicated than the invariant-manifold governing equations represented in modal position and velocity, and have potential singularities at zero modal amplitude. However, as mentioned in [33], the constraint equations of the slave state variables are periodic in ct], with period 2rr, and therefore a Fourier basis can be used for shape functions in the (bk direction. The methods for obtaining approximate solutions for P,- and Q,- are described in section 4.2.2. 2.2.2.2 Multi-Mode Manifold Formulation The multi-mode manifold formulation in terms of modal amplitudes and phase angles can 24 be developed in the same manner as was done for the case of modal displacements and velocities. To avoid redundancy, only the important steps are given here. The trans- formations relating the master modal position and velocity (77k, in.) to the master modal amplitude and phase (ak, (bk) are given in equations (2.13) and (2.14); these are employed with k E S M, where SM has the same definition as before. The sets {ak} and {m} are denoted by a M and d) M respectively. The remaining slaved states are expressed as 7h = PilaMs (PM) (223) 7h = Qi(aMa (PM) (224) for i=1,2,3,...,Q, i¢ SA]. Equations (2.17) and (2.18) are still valid with j if S M. Next, we take a time derivative of the constraint equations, eliminate the time dependence by using the first—order equations of motion for the master and slave modes, and substitute the constraints every where in place of the slave state variables. This yields the following 2Q — 2M equations for the mnlti—mode invariant manifolds: 25 Qt _ Z [911 (-fi.~(Pj»Qj.ak.ak) sums“) kESAI (90k wk +11 (Wk _ f1.( JQJ a}. <91.) C05(C’A))] (225) 30k ak Wk —w,-2 P7; +fi(Pj,Qjaak,¢k) : Z [991(-fk(Pinjiak~¢L-) Sin(¢kl) k E S A! 80 k wk +6121 (wk _ fi( JQJ 0]. Oi.) 009090)] (2.26) dcfik ak wk. for i,j=1,2,3,...,Q, i,j¢ SM. As before, closed-form solutions do not generally exist. A method for obtaining accurate approximate solutions for the P,- and Q,- is described and implemented on example systems in [42]. 2.3 Solution of the Invariant Manifold Equations In general, solutions of the invariant-manifold equations cannot be determined in closed form, since such a solution represents a family of solutions to the original equations of motion. However, in some special cases involving symmetries, closed-form solutions, which typically represent flat manifolds, can be determined ([14]). More generally, though, the manifolds are not flat, and other approaches are required (such as those described in sections 2.3.1 and 2.3.2). In the early works by Shaw and Pierre ([70], [20], [21]), approximate solutions for the NNM manifolds of smooth systems were obtained in the. form of asymptotic 26 series expansions. Shaw, Pierre, and co—workers ([72], [73], [41], [28], [29], and [71]) used this approach to solve a variety of problems, as did others, notably [26], [27], [30], [18], [16], and [19]. However, the asymptotic series approximation of the manifold is only locally valid; therefore, if one tries to Simulate a NNM model at large amplitudes, the time response of the NNM model will deviate from the time response of the original model. The primary difference is in the frequency of oscillation, which results in a phase drift. As mentioned in section 2.1, the work by Pesheck et al. ([33]) employs a Galerkin method to solve the invariant manifold equations, which allows one to obtain accurate approximate individual NNM manifolds up to large amplitude ranges. This approach is also not limited to smooth systems, so that the nonlinear terms can be more general functions. However, the computational cost of acquiring such solutions can be quite expensive. Therefore, efficient means of computing coefficients in the Galerkin expansions were sought and, in particular, the collocation method has been found to provide accurate solutions with much greater computational efficiency. The series expansion method is reviewed briefly in section 2.3.1, for both individual and multi-mode NNM models. The Galerkin and collocation methods for individual NNM models are reviewed in section 2.3.2. However, the Galerkin method used to solve multi-mode NNM manifolds, developed by Jiang et al. ([42]) is not reviewed here. Interested readers are referred to [42]. 2.3.1 Solution by Asymptotic Series Expansion A good summary of the series expansion method for individual and multi-mode NNM manifolds was given in [74]. A brief version of the method is given here. To use the series expansion method, the nonlinear terms in equation (2.2) are assumed to be second and third degree polynomials in the modal displacements. This assumption can be generalized 27 such that velocity-dependent and mixed nonlinear terms are included, but the solutions will be more complicated ([20]). Therefore, the form of the forces f,- in equation (2.2) is given by Q Q Q Q Q f1 = —ZZa1pq ups—2:23am 11111., 11. (2.27) p=1q=p p=lq=PT=q 2.3.1.1 Solutions for Individual NNM Manifolds To construct the single—mode manifold, approximate local solutions of X,- and Y,- can be found using a power series expansions in terms of uk and vk as X, = afiuk + 03,171}; + aging, + allaukvk +a§7ivg + aging + afiuzvk + ogiukvg + agivg + . .. (2.28) Y1 = b11111: + birl’k + biaui + bi,1’uk'l’k +b§ge§ + 115,112 + 11%,.11711, + 115,111,112, + 1153,11}: + . .. (2.29) The a’s and b’s are the unknown coefficients to be determined, and there are 9 x (262 — 2) of them. These coefficients can be determined by substituting equations (2.28) and (2.29) into equations (2.7) and (2.8) and collecting like powers in uk and 11),. This yields a set of linear equations for the unknown coefficients, which can be solved successively. The first-order (linear) coefficients are solved for first, and they are all found to be zero. This is because equation (2.2) is given in terms of the linear modal coordinates, and therefore the linear parts are already uncoupled, and must remain so since the linear model satisfies invariance for its modes. Next, the 3 x (262 — 2) second—order coefficients are obtained, and they are required in the solution of the 4 x (2Q — 2) third-order coefficients (as is typical in expansion solutions). General formulas for the as and b’s can be found in [74]. 28 Once all of the expansion coefficients are obtained, the X’s and Y’s, which describe the slaved modes, are known functions of the master states (11k, vk). For i = k, these known functions are used to express fk in equation (2.2) in terms of only uk and 11),, rendering a single-DOF oscillator as the reduced-order-model for the kth NNM. The coeflicients of the polynomials in the expansions become Singular when there exists low order resonances between modes. For the quadratic and cubic terms, 1:1, 2:1, and/or 321 resonances ([79]) between the linear natural frequencies of the master and slave modes will give rise to these Singularities. In these cases there exists nonlinear coupling such that energy is exchanged among the resonant modes, and therefore these modes cannot be dynamically separated from one another, as is required for the construction of an individual NNM model. We should note that in this case the response of an individual N NM is not stable, and any perturbation will cause the response to leave the neighborhood of the NNM manifold, as the energy flows from that mode to others (and in most cases, back to that mode again at a later time). In these situations, the multi-mode formulation described in section 2.2.1.2 is required, with the necessary inclusion of the resonant modes. The series solution of these multi-mode manifolds is described in the next section. 2.3.1.2 Solution of Multi-NNM Manifolds Here there are M pairs of master displacements and velocities. The assumed forms of the slave constraints X ,- and Y,- in terms of 11k and U}, are given by 29 k k k.l k,l k,l Xi(uM,vM) = 2 (anti), + a2,,ei1k)+ Z Z ((23,,l-Itkll11‘l' a4’i’kaUl + 05,11‘kvl) kESAy ICES)” [65]” + Z Z Z (algj‘qukujuq + aéj‘qukujvq + agj’qukvlvq ICES)” lESM (1681” + ugli’qvkvlvq) + . .. (2.30) 1(uM. vM) — Z ( 1,114 + 2,114.) + Z Z ( 3,111,11111 + 4,1141% + 5,2-vkv1) kESM ICES!” [63M k 1, k1, '.l. + Z: Z Z (b6:1‘ qukuluq + biz-qukulvq + bggi'qukvlvq kESAJ 165)” (1651” + bgjg‘qvktqvq) + . .. (2.31) The as and b’s are the unknown coefficients to be determined, of which there are (2M + 3M2 + 4M3) x (26.2 -— 2M ). These coefficients can be determined by substituting equations (2.30) and (2.31) into equations (2.7) and (2.8) and collecting like powers in the master displacements and velocities. This yields a set of (2M + 3M 2 + 4M3) x (2Q - 2M) linear equations for the (2M +3M2 +4M3) x (262- 2M ) unknown coefficients, which can be solved at successive orders. Since the multi—mode invariant—manifold equations are formulated in terms of the linear modal coordinates, the linear terms in the above expansions are again zero. The formulas for the remaining a’s and b’s can be found in [74]. The coefficients a’s and b’s in equations (2.28) and (2.29) of the individual NN M manifold formulation become Singular when w,- : wk and/or to,- = 2112;, and/or w,- = 3wk. This is a flag indicating that there exists 1:1, 2:1, and/or 3:1 internal resonance(s) among slave modes 2' 75 k and master mode is: ([79]). In this case there exists energy exchange among the resonant modes and the motions necessarily take place on higher-dimensional (dimension greater than two) manifolds. Once all of the participating modes in the internal resonance are included in the set S M, the singularities in the coefficients a’s and b’s will be removed (at least up to polynomial degree 3) in uk and 11),. Then we obtain the multi-mode invariant manifold that captures the internal resonance, and allows energy exchange among the 30 master modes. Sin'iilarly, there are possibilities that the coefficient a’s and b’s in equations (2.30) and (2.31) of the multi-mode manifold formulation become singular when 112, 2 wk and/or to, = 2wk and/or w,- = 3a)];C and/or 112,- : [wk ital] and/or to,- = [2112), ital] and/or to,- 2 [wk in); :Ewml. This is another flag indicating that there exists 1:1, 2:1, and/or 3:1 internal resonances among a Slave mode i d S M and a master mode k, l, m E S M- The singularities in the coefficients a’s and b’s can be removed at least up to polynomial degree 3 (third order) in uk and 11],, by redefining the new set SM such that it includes those extra modes. Once all of the expansion coefficients are obtained, the X ’s and Y’s, which describe the slaved modes, are known functions of the master states (11 k1 vk), where k E S M- For i = k, these known functions are used to express fk in equation (2.2) in terms of only the uk and vk, rendering a M -DOF oscillator. which is the reduced—order—model that includes the desired NNMS. 2.3.2 Solution by Galerkin and Collocation Methods The Galerkin and Collocation methods are techniques that employ a. variational formulation to solve for discrete approximate solutions of continuous problems. Both methods belong to the general class of techniques known as weighted—residual methods. A good account of the theoretical aspects of these methods, and examples of how to apply them to engineering systems, can be found in [77]. Many useful guidelines and rules of thumb on how to use the methods effectively can be found in [78]. A brief summary of the methods is given here. Consider the operator equation of the form 31 A(u(a:)) = f(;1:) (2.32) in an open domain Q with apprOpriate homogeneous boundary conditions, where :1: E Q is the independent variable, A is an operator, u(r) is the dependent variable to be determined such that it satisfies the operator equation and the boundary conditions, and f (:13) is a known function that represents a driving term. In the weighted-residual method, u(:r) is approximated by uN(a:) in the expanded form N 1e) 2111(1) = Z ea) (2.33) i=1 where {dz-(12)} are elements of a complete set of basis functions that are differentiable up to the order of differential operator A and satisfy all homogeneous boundary conditions, and the c, are constants to be determined. Substituting equation (2.33) in (2.32), we obtain the residual, R N RNl$1 C1) Al'urvtv» — f(~’r), (1334) which is generally not zero, since uN(r) is an approximation of u(.r). Note that uN(a:) lies in a finite-dimensional space since uN(;1:) is represented by a finite linear combination of elements of the basis. However, u(a?) lies in infinite—dimensional space, since to represent u(.r) all elements of the basis are generally required. If {tb,(:r)} is a complete set of basis functions (possibly different from {gb,-(.r)}) in the space 32 that A(u(;r)) and f(:r) live in, then by projecting RN(a:,c,-) onto 111,,(1) using the usual L2(Q) inner product, and equating the resulting terms to zero, we obtain = fflRniiricz'W’kffldx = 0 (235) If the set {'1l',(:r)} is an orthonormal basis, 7‘}, = < R N(:r, oi), u’1k(:1:) > is the representation of RN(:r, Ci) in that basis. We then have N algebraic equations with N unknowns, the {C1}. Solving for the {C1} that satisfy equation (2.35) would make RN(a:,c,-) small over $1. In the limit N —1 00, R N(r, ci) must go to zero. This method is called the weighted-residual method. If one chooses different basis functions for expansion and projection, 1141A. # pk, the approach is known as the Petrov-Galerkin method. If one chooses the same basis functions wk 2 (pk, it is known as the Bubnov-Galerkin method, or, more commonly in the West, the Galerkin method. If one chooses "l’k : 6(r —:rk), where 6 (:13) is the Dirac delta function and it], are N selected points (called collocation points), this becomes the collocation method. The collocation method will force the residual R N(a:, c2) to be zero at points 313;, in (I. This greatly speeds up the calculations, due to the simplification of the integrals. In the following sections, the Galerkin method will be applied to equations (2.21) and (2.22) in a domain ak E [0,a0] and dk E [0,27r], and also to equations (2.7) and (2.8) in a domain uk 6 [—Ub, U5] and vk E [—Vb, Vb]. The collocation method will be applied to equations (2.21) and (2.22) in a local domain defined later, and also to equations (2.7) and 33 (2.8) in another local domain defined later. 2.3.2.1 Galerkin Method with Global-Basis Functions in Modal Amplitude and Phase Angle Pesheck et al. ([33]) applied the Galerkin method to solve equations (2.21) and (2.22) for P2: and Q,- in the domain (1;, E [0, no] and 1,5,, 6 [0, 2r] using an expansion of basis functions. A brief summary of the approach in [33] is reviewed here. The unknown position and velocity constraints (P,, (2,) are expanded as a double series in the master modal amplitude ark and phase (bk as Na N92 , l.- P1<111e1f> = [Zamnmeeea (2.36) l=1m=1 Na N11 1 Q1(a~i.-1¢1~.) = ZEDi’mUz,1—n(ak1¢k) (2-37) l=1rn=1 for i=1,2.3,...,Q,i71- k, where sz‘m and Dg’m are coefficients to be determined, T1,",(abqbk) and Ul,m(ak~ m) are known basis functions. The basis functions T1,",(ak. m.) and Ul,111(aka¢k) are products of selected basis functions in the (1k and (bk directions, and they are defined over the domain 61k E [0, a0] and 1.9;, E [0, 27r]. Na and Ng are the number of shape functions used in ak and Qk respectively. Since. (pk has period 2rr, the basis functions in the (1);, direction are chosen to be the Fourier basis, {1, cos(ncf)), sin(no)}. The basis functions in the (1;, direction are Chosen to be polynomial functions, L1(ak). Therefore the form of T17m(ak,q)k) and Ul,m(ak1¢k) are given by 34 71111011119511) = LAG/1) 608((771-1)¢k) (2.38) Ul,m(akv¢>k) = L1(ak)sin(m¢k). (2.39) It can be shown that for a conservative non-gyroscopic system to posses synchronous mo— tions only the cosine functions are needed for basis functions in 111,, for Tim» and only sine functions are needed in 115k for Ul,m- The polynomial functions L,(ak) are chosen to be a set of polynomials defined over the domain ak E [0. a0] with zero value and zero slope at ak = 0, which will satisfy the conditions that the invariant manifold pass through the ori- gin (ak, pk) = (0,0) (the stable-equilibrium point) and be tangent to the two-dimensional eigenspace at that point . Also, for convenience, they are chosen to satisfy the orthogonality condition “0 (1.", 1 for 'l :j / —2- Li(ak) Lj(ak) dak = (2.40) 0 (10 O for i 75 j, which can be achieved by starting with a set of polynomial functions “k 2 0k 3 “k {(—) ,(--) .(—— a0 a0 00 of the L1(ak) can be found in [33]. Equations (2.36) and (2.37) are substituted into the )4, }, and then applying the Gram—Schmidt process ([80]). The form invariant-manifold governing equations, equations (2.21) and (2.22), and each of these is projected onto each basis function using the Lg-inner product over the entire domain. This leads to 35 or, _ 0 ___ / (1_2 )Up.q ‘ak EDI, n'liU 1,111+ +20%, 111 l. n( fk a}; SiIl(¢1)) a k ‘f’k 031,111 l,m 8(1k(wk 0T . +2: Cl 171 lm( l. 1.1.111. - M) dak (lgbk (2.41) [771 wk 1 1 8U1 —f a sin (b ) 0 = / <—2> T11 .. aizC’ ",—T1m 111 “2080111 U; (k) “A ilk (’0 l, m l. m k k 0U . 1 + 2: D1 m [m( 1. , wk — ————f]‘ (02%)) dak 111,11. (2.42) l m. wk for i: 1,2,3...-,Q, 1% k; The integrations in equations (2.41) and (2.42) are computed numerically. Therefore, equa- tions (2.41) and (2.42) result in a set of 2(Q — 1)N1,N11, nonlinear algebraic equations in the C ’s and D’s. For simplicity in referring to the method subsequently, this method is named “the global a —— ct approach”. Note that for the special case of a conservative non-gyroscopic system with only cubic nonlinearities in the modal positions, all harmonic terms in the expansion of basis functions for each state are needed. Therefore, the number of nontrivial coefficients remains 2(Q — 1)NaN¢. This fact is observed by performing numerical experiments. In contrast, for the global u — 11 approach, this special case results in a reduction in the number of coefficients, as described in section 2.3.2.2. Once all of the expansion coefficients are obtained, the PS and Q’s, which describe the slaved modes, are known functions of the master states ((11,31). These known functions are used to express f1. in equations (2.17) and (2.18) in terms of only a1. and (111., rendering 36 a single-DOF oscillator as the reduced-order-model for the km NN M. 2.3.2.2 Galerkin Method with Global-Basis Functions in Modal Position and Velocity In this section, we apply the Galerkin method to solve equations (2.7) and (2.8) for X,- and Y,- in the domain uk 6 [—U1,,U1,] and v1, 6 [—V1,, V1,] using an expansion of polynomial basis functions. For nonlinear structures, the form of nonlinear terms is often a polynomial function, and solving the invariant manifold equations in modal position and velocity would have advantages over those in modal amplitude and phase angle. This is so Since the terms that are necessarily zero from symmetry , i.e., the trivial terms, are easily identified and can be removed from the expansion of basis functions. Moreover, the integration of polynomial functions can be computed exactly using Gaussian integration, while the integration of sine and cosine functions in equations (2.41) and (2.42) cannot be computed explicitly. Therefore, the integrands in equations (2.41) and (2.42) have to be evaluated at many points, which in turn increases the computational effort of the procedure. The unknown position and velocity constraints (X1, Y1) are expanded as a double series in the master modal position 111,. and velocity vk as Np,u Npgv l X1('tik,vk) = Z Z Czi’lel/I1,m(uk,111.) (2.43) [:1 m=1 Nan Nprv 1’. 3301111011) = Z ZDiijlIlflnhl'kivk) (244) 1:1 "1:1 for i = 1,2,3, ...,Q, i 729 k, l,m l,-m . . where C,- and DI. are to-be-determmed coeffic1ents, and the L1,",(1.1.k.vk) are known 37 basis functions. The basis functions L1,",(uk,vk) are chosen to be two—dimensional polynomial-basis functions in uk and 11k such that L1g.,,,(0,0) = 0, (31m (0 0): 0, 11k and #(0,0) = 0 , which insure that the invariant manifold passes through the ori- Uk gin (11.1., v11) = (0, 0) (the stable equilibrium point) and is tangent to the two-dimensional eigenspace at that point. The construction of LM1,m(u1.,vk) is obtained by performing the tensor product of two sets of polynomial functions given by 1_11 11.2111. 3 11.1. 1_1 L = 1 —— — ,... 2.45 (I’kb (L7; 2 “A: 3 ”k m—l [I : Q 0.. 0.. O I v . u . The terms 1, (vi), and (F2) are eliminated from the set {Lh11’m(uk,vk)}, since they do not satisfy the boundary conditions at the origin. The Gram-Schmidt process is applied to the reduced set {LML,,,(111.,v1.)}, such that the following orthogonality conditions are satisfied: Vb Ub 1 1 1 for [771 2 pg / / (U)(V) L’ll1m(u1.,11)L'llpq(u1.,1k)duk (111: (2.47) Vb Ub b b 0 for lrn. sé pq. Samples of the {LA/1,",(uk, 111)} can be found in Appendix 2A. Equations (2.43) and (2.44) are substituted into the invariant-manifold governing equations, equations (2.7) and (2.8), and each of these is projected onto each basis function using the Lg—inner product over the entire domain. This leads to: 38 1 1 lm 0 = f —, — LAI, — 0’ LM 11k;11k( lb) (V1,) pq Z i l,m (l-l-m.)yé2,3 1.171 aLi’l/Mn l.m CLAIM" 2 't' 2 Ci ’Uk W 'i- 2 Ci W(—wk ltk + fk) dllk dUk (l+m)¢2,3 (l+m)7é2.3 _ Yrq : Fz‘ (2.48) _ 1 1 , 2 1,171 0 _ (-U—) (V) LMW .1,- 2 C, LMW — f,- aka/'1. b b (l+m)7é2,3 l.m 814])!th [.m aLjubm 2 ‘i‘ Z _ Di l'k T + 2 Di w(—wk 11k + fk) d111, d‘l’k (l+m)¢2,3 (l+-m);£2,3 s Ff“ (2.49) for i = 1,2,3,...,Q, i 79 k; p :1""’1N1)‘u; 7 q :1,...,1‘\lp.11; and (1) + (1) 7f 2, 3 . The integration in equations (2.48) and (2.49) is computed numerically. Therefore, equa- tions (2.48) and (2.49) represent a set of 2(Q — 1)(1’Vp,u1’Vp,v — 3) nonlinear-algebraic equa- tions in the C ’s and D’s. For simplicity in referring to the method subsequently, this method is named “the global u — 11 approach”. For a conservative non-gyroscopic system to posses a synchronous modal response, some of the C ’s and D’s must be zero. The basis functions corresponding to such C ’s are the ones that are generated from 39 “Exit we). (“433(3), (“—kW—h, (fire). It, 31:3 36:23 “$233 $323 E5433 35533 I) . 11 ’U . U. . U . ll . ‘U . ll, . ”U . ”(1. U. —‘— 5, <—" —* 5 i 2 —")5 —"—>3 —">5 —‘— 4 ——‘—>5 i— 5 —")5, ...}, (2.50) Vb Ub Vb ’ Ub Vb ’ Ub Vb ’ Ub Vb ’ Ub Vb . ‘ , . I . 'Uk 1 NA. 3 'l’k 5 ‘3‘ . where the elements in the set are linear in (—) , (—) , (———) , The basis functions Vb Vb Vb corresponding to such D’s are the ones that are generated from on, {(302, (E93, (”—w, (“—k>5,. Ub Ub Ub Ub 952 a a2 3.23:2 a3a2 £5432 32.532 1&4 ”_k 21:4 2521;“ 113114 E434 “4524 251 Vb ’(Ub Vb ’Ub Vb ’Ub) Vb ’(Ub) Vb ’Ub) vb) "'"} (' ) . . , v . v. v . where the elements in the set are linear 1n (~15)0 —A 2 a ($7504 b Vb ’Vb the number of nonlinear algebraic equations to be solved for the remaining 0’5 and D’s is (Q _ 1)(J’Vp‘u]\’rpyv "' 3). , Therefore, in this case For a conservative non-gyroscopic system with cubic nonlinearities, the additional C ”s and D’s that correspond to l +m 2 an even number, would all be zero. Therefore, the number of nonlinear-algebraic equations in the C ”s and D’s to be solved is about 15(6) - 1)(i’Vp,uNp,.,/. — 3). Once all of the expansion coefficients are obtained, the X ‘s and Y’s, which describe the slaved modes, are known functions of the master states (uk, vk). For i = k, these known functions are used to express fk in equation (2.2) in terms of only uk and ck, rendering a single—DOF oscillator as the reduced-order—model for the km N NM. 40 2.3.2.3 Collocation Method with Local-Basis Functions in Modal Amplitude and Phase Angle (Annular Strip Representation) In [33], Pesheck et al. mentioned that as the number of retained linear modes Q increases, the computational time associated with solving equations (2.41) and (2.42) increases drasti- cally. Therefore, they developed a new approach to solve the invariant manifold equations. In this approach, the desired domain in the (IA. direction is divided into K small sections, producing annular subdomains given by (35}: E [0, 27r] and ak E [Cl/(J, (Maj + Aak]. For each subdomain, as before the shape functions in the qbk direction are chosen to be the Fourier basis. However, the shape functions in the ak direction are chosen to be piecewise linear segments, which are accurate, providing the subdomains are sufficiently small. Therefore Pi(ak, (bk) and Q,(ak, (bk) can be expressed over the jth interval as: . 1,- ,. Paar-9 99k) = Z Ci‘mfl-mUl/c» at) 1,171 NO 0) (1k - l,m ‘ _ “,j = c, __ z [ . < > m=l , a . — (1.,- +C.2""(1 — i—Ei) cos((m — n.3,.) (2.52) 1' Auk . l, Qliak’gok) : ZDimUlmlmke‘rl’k) Lm Né 1 m 0k — Gk], 2 m 0k "' akj = 2 [DIN (——'——) + Dz? (1 — —) sin(mc;')) (2.53) 77121 Aak Auk for i:1,2,3,...,Q,i7$k. The Galerkin method is then applied over each subdomain. Therefore, equations (2.41) 1 . Y and (2.42) without the scaling factor (7) are valid for each subdomam as well. Note that. “'0 41 there are K sets of 0’5 and D’s, one for each Aak interval. These individual solutions are assembled to construct the invariant manifold. Note that this approach was applied to construct the accurate ROM of a nonlinear rotating beam in [34]. The computational effort associated with evaluating and solving equations (2.41) and (2.42) can be quite high. To reduce it, Apiwattanalunggarn et al. ([35]) implemented the collo- cation method to construct an accurate ROM for a nonlinear finite element model of the rotating beam. For this implementation of the collocation method, the nonlinear-algebraic equations (similar to equations (2.41) and (2.42)) are given by - _, . l,m () = / ()(a;C -— ”1.31)st — (Dhql) “'01; 2 Di Ul,m (‘ki‘f’k 1,771 +2 Cl, m 0T! III (—fk a}; Si“(ok)) (In aak< wk + Z Cl md _l_0m( kwk — M) dak dqbk (2.54) (In wk - l. . 0 = / , 001k - ak,p~¢k - 65m?) w? 6% ZCi‘mTIm - 0k fi ah ¢k l,m [JuaUl m (“f/,- ak 5111(6)” +2 0 a . . > I In ”k wk . 0U . co.‘ '. + 2 D9," 0;:1 (aka) wk— iii—M) . dak dcfik (2.55) w . [.In I" for i=1,2,3,...,Q,i# 1:; 13:12. ()2 =1 "NC“ where (ak.p~¢k,q1) E [(Lk‘j,(lk‘j + Auk] x[0,27r] are collocation points associated with 42 Qi(ak,¢k) and (ak,p~¢k.q2) E lak,jaak,j + Auk] ><[0,27r] are collocation points associated with [Dz-(alt, an). The collocation method greatly simplifies the integrals. Equations (2.54) and (2.55) are a set of 4(Q —- 1)N¢ nonlinear-algebraic equations in the C ’s and D’s. Recall, that one must solve K such sets of equations to obtain the manifold over the entire domain. For simplicity in referring to the method subsequently, this method is named “the local a — qb approach”. For a conservative non-gyroscopic system, only half of the harmonic terms in the expan— sion of the basis functions for each state are needed. Therefore the number of nontrivial coefficients is 2(Q — 1)N¢,. This fact is known from the symmetry of such solutions, and is also confirmed by numerical calculations. Once all of the expansion coefficients are obtained, the P’s and Q’s, which describe the slaved modes, are known functions of the master states (ak, (bk). These known functions are used to express fk in equations (2.17) and (2.18) in terms of only ak and (bk, rendering a single-DOF oscillator as the reduced-order—model for the kth NNM. 2.3.2.4 Collocation Method with Local-Basis Functions in Modal Position and Velocity (Patch Representation) The approach developed in section 2.3.2.2 has the same disadvantage as that described in section 2.3.2.1 when the number of retained linear modes Q increases. Hence, in this section we develop an approach similar to the approach given in section 2.3.2.3. However, we take the slightly different approach from that of section 2.3.2.3, wherein we use rectangular subdomains defined in the modal displacement and velocity. The solution is expanded over each such subdomain using two-dimensional polynomial basis functions. Since these subdomains are small, we can use low degree polynomial functions to describe the invariant manifold over the subdomain, which significantly reduces the computational time of solving for the invariant manifold. 43 Here we describe the development of the method for a local domain described by u): E [—ub, ub] and v}: E [-'ub, vb]. Once this procedure has been developed, it can be applied to each patch, and the final result is obtained by collecting the results for all patches such that the entire domain is covered. The relations between the coordinates in the global domain (uk, ck) and those in the local domain (’uz, vi) are given by (I). = (1” +112, (2.56) L’k : (11} + vi, (257) 2U where (1“ : —Ub + ( N:)(eu — 1) + (1),, (2.58) u r 2"?) '— and (iv = —l’b + (WW2), — 1) + vb, (2.09) L! where eu and 6,, are the patch indices in the uk and ck directions respectively, and N5 and N5 are the number of patches used in the uk and ck directions, respectively, and (in and ~ - - , e .,e e (11, represent the shift from the or1g1n of (uk, wk) to (uk, bk). Here cu runs from 1 to Nu and ev runs from 1 to NS. Substituting equations (2.56) and (2.57) into equations (2.7) and (2.8), the partial differen- tial equations governing the geometry of the kth manifold in the local coordinates (iii, vi) are given by: 44 —w,-2 Xi + fiinanadu + afidv + 12f.) for 6X- 3 - ’(d. + vi.) + ’(—w£ (d. + at) 57; av; + fk(XJ-. Yjadu + ui..dv + vi.» (2.60) ar- 01/- ggéw. + vi.) + 51%(_°°’12c (.1, + at) + fk(XJ-, Yj,du + 11):, dv + 1%)) (2.61) 2'.I=1,2,3.....Q. zyjsék. The solution of equations (2.60) and (2.61) on the local domain is obtained by expanding (X i: Y,) in terms of basis functions as: 1Vpau Npav x,(uf;,u,f) = Z ZCf'mLML7,,(ui,vg) (2.62) [:1 m=1 NW NW Yi(ui.,v£) = Z Z Dg’mLA/Il’mhtz,vi), (2.63) [:1 m=1 where Lfl>ll,,.,,(u‘):;,'v;) = Tl_1(ui./ub)me_1(v,‘;/vb), (2.64) where T1_1 (:L‘) and Tm_1 (.1) are standard Chebyshev polynomials defined over :17 E [—1, +1], and the C ’s and D’s are the to-be-determined expansion coefficients. Equations (2.62) and (2.63) are substituted into the local manifold-governing equations, equations (2.60) and (2.61). Normally, each of the resulting equations is projected onto the basis functions, but here we employ a collocation method, which is computationally more efficient, yet retains very good accuracy. This is carried out by projection of the equations onto Dirac delta functions in the local master coordinates over the local domain, as follows: 45 for , 1.171 / e 6(112 — 1124,. 1'; - 'vi‘q) — 2 D, LMI,rIz '“Z’Uk l,m e BLAIIJH "i”:Cz- (d1) +Uk)—_8—(;Z—- l,m ‘ Fri/J”! l l,m BLAIM, Lm 2 + Z Ci w(—wk(du ‘i‘ HZ) + fk) duck (iv/i: 1,171 (2.65) , _. . 2 [.771 , / e 5(112 _ {1.24), v; — vixq) “’2' 2C1 Ljubm _ f2. 6 .,,. 11k,1k 8 0141111.," 1m +2 DY d»,+1..’.—,—.— l (i 1‘) 011;; lm RX 1%] Z i=1,2,3,...,Q,i7€k; q :19 ..., AGLU, 1,171 BLI’ULm Lm . 2 , € . e ,e ’i‘ E: Di W(—wk(du + 11k) + fk) duk dl'k ,m (2.66) where (ui‘pwg‘q) E [—ub,ub] x[—vb,vb] are the collocation points, which are zeroes of TNp‘u(uf;, / 11b) and er,L,(Ui)~/vb)~ respectively ([78], [81]). Equations (2.65) and (2.66) con- stitute a set of 2(Q — 1)} nuNpm nonlinear equations in the C ’s and D’s. Note that there are N5 x N5 sets of Cs and D’s. However, if the system is conservative, non-gyroscopic, and the nonlinear terms are functions of solely the modal positions, then only :14 x N5 x N5 sets of C’s and D‘s need to be obtained, and the remaining coefficients can be generated using symmetries in the solution. For simplicity in referring to the method subsequently, this method is named “the local 21 -— 11 approach”. Once all of the expansion coefficients are obtained, the X ’s and Y’s, which describe the slaved modes, are known functions of the master states (1123,, iii). For 2' = k, these known 46 functions are used to express fk in equation (2.2) in terms of only 11% and vi, rendering a single-DOF oscillator as the reduced-order-model for the km NN M. 2.3.3 Comparison of Computational Efforts To compare the computational efforts of the approaches described above, there are three main aspects that we must consider: the total number of unknown coefficients, the number of evaluation points needed for the integrands in order to compute the integrals, and the methods used to solve the nonlinear-algebraic equations. First, the reduction in the number of unknown coefficients C’s and D’s for a specific class of systems of all four approaches (the global a — 66, the global u — v, the local (1 — <25, and the local 11 — v approaches) is described. Next, a. comparison of the number of unknown coefficients 0’8 and D’s and a comparison of the number of evaluation points needed for the integrands of the global approaches (the global a — (b and the global u — v approaches) are described. Then, a comparison of the number of unknown coefficients 0’3 and D’s of the local approaches (the local (1 — (b and the local 11 — v approaches) are described. This section is closed with a comparison of the methods used to solve the nonlinear-algebraic equations in this study. Summaries for the total number of unknown coefficients and the number of evaluation points for the global-domain and local-domain approaches are shown in Tables 2.1 and 2.2, respectively. 2.3.3.1 Coefficient Reduction for a Specific Class of Systems For the most general class of systems, that is, those with general (i.e., non-Caughey) damping, gyroscopic effects, etc., one must obtain the entire set of C and D coefficients. However, in some commonly encountered cases the number of coefficients can be reduced simply due to symmetries in the equations of motion and the resulting responses. These symmetries may be in terms of temporal or spatial dependence. Specifically, consider the 47 forms of equations (2.41), (2.48), (2.54), and (2.65). If the nonlinear terms are functions of the modal positions only, it can be seen that the D’s can be expressed as functions of the C’s. By using this result in equations (2.42), (2.49), (2.55), and (2.66), the number of equations and the number of unknowns in all four approaches can be reduced by half, thereby reducing the computational time accordingly. We now turn to the details of the four approaches. 2.3.3.2 Comparison of Coefficients and Evaluation Points for the Global- Domain Approaches Here a comparison of the number of unknown coefficients C ’s and D’s of the global a — (b and the global u — v approaches is described first. This comparison is based on the systems having nonlinear forces of polynomial types. The numbers of unknown coefficients 0’3 and D’s of the global a — (b and the global u — v approaches for this class of systems have been described in sections 2.3.2.1 and 2.3.2.2. Second, a comparison on the reuse of previously-computed integrals, which is based on how the slave states are represented by the expansion of basis functions, is described. Then comparisons of the number of evaluation points needed for the integrands in order to compute the integrals of the global (L — 45 and the global u -— v approaches are described. For the global-domain approaches, solving the invariant—manifold equations in terms of modal position and velocity has three advantages over solving them in modal amplitude and phase angle. First, for systems having specific polynomial nonlinear terms, either quadratic or cubic, the trivial basis functions of the global u — v approach (mentioned in section 2.3.2.2) can be easily determined, and therefore the number of coefficients can be significantly reduced. However for the global a. — (I) approach (mentioned in section 2.3.2.1), we can not take such an advantage from this specific class of systems. 48 Second, the basis functions of the slave positions and slave velocities, expressed in terms of the master position and velocity, are the same, therefore equations (2.48) and (2.49) share the same integrals. However, this is not true for equations (2.41) and (2.42) for the global a — 65 approach, since the basis functions of the slave positions and slave velocities expressed in terms of the master amplitude and phase angle are different. Third, as mentioned in [82], the number of evaluation points in the (bk direction of the global a — (I) approach is taken to be at least 10N¢ to provide good results, and typically M, 2 12. The number of evaluation points in the 1);, direction of the global u — 11 approach, for which gaussian integration ([80]) is used, is (ngxv — 2). Assuming that the number of terms used for the modal displacement is comparable to that used for the modal amplitude, i.e. Npru % Na, therefore the number of evaluation points in the uk and the 0.), directions, for which gaussian integration are used, are comparable, i.e., (ngm — 2) z (3N0 + 3). Hence, one can make a comparison by considering the number of terms needed for the modal velocity versus those needed for the modal phase. If NW, is assumed to be 12, which is the minimum number required for N¢, it can be seen that the number of evaluation points in the 63;, direction is much higher than in ii), direction, i.e., (10M), = 10 x 12) >> (314,,“ — 2 = 3 x 12 —- 2). 2.3.3.3 Comparison of Coefficients for the Local-Domain Approaches Here, a comparison on the reuse of already-computed integrals, which is based on how the slave states are represented by the expansion of basis functions, is described first. Second, a comparison of computational times, which is based on how many terms in the expansion of basis functions are needed to sufficiently describe the slave states, is described. Note that the number of terms in the expansion of basis functions is determined from experience gained from numerical experiments. Note that a comparison of the number of evaluation points needed for the integrands in order to compute the integrals for the local a.—— and the local 'u — v approaches is not considered since both approaches use the collocation method. 49 First, the basis functions for the slave positions and slave velocities, expressed in terms of master position and velocity, are the same, in the local 21 — I) approach, therefore equations (2.65) and (2.66) share the same integrals. However, this is not true for equations (2.54) and (2.55) for the local 11 - (f) approach. Second, the numbers of basis functions in the a), and uk directions for the local a — 45 and the local 11 — v approaches, respectively, are comparable to each other since the manifold is described in the “k and U), directions using local basis functions and typically Na = 2 and typically NW) 2 3. Many harmonic terms, typically N¢ Z 12 (as previously mentioned), are needed to describe the manifold globally in the 45k direction. In contrast, only a few Chebyshev-polynomial basis terms, typically Np,” : 2, are needed to describe the manifold locally in the 2);, direction. Therefore, the dimension of each subproblem of the local 11 — v approach is much smaller than the local a — 65 approach, which implies that the computational time for each subproblem of the former method is much less than that for the latter method. Even though for each grid (piece) along the uk direction many pieces along the 1);, direction are needed, the total computational time for each grid (piece) along the u), direction is just the sum of the computational times for the pieces along the 1);, direction. However the computational time for each grid (piece) along the a), direction does not grow proportionally with the number of harmonic terms N¢, i.e., it grows nonlinearly with Ngb' Hence, there is a tendency that the local 21. — 12 approach uses less computational time than the local 0. — (25 approach. Note that this is only an educated estimate. In section 2.4.2.1, the computational times for the local a — q!) and the local 11 — v approaches will be measured and reported. Note that for both local-domain approaches applied to conservative non-gyroscopic systems, the conservative condition and the associated symmetry in the manifold solution are fully employed. 5O 2.3.3.4 Comparison of Methods Used to Solve the Nonlinear-Algebraic Equa- tions In this study, three methods have been implemented to solve the system of nonlinear- algebraic equations for the unknown coefficients. They are Powell’s Hybrid method ([83]), the Newton-Raphson method ([80]), and the Secant method ([80]). All three methods have been used in this study, and by previous researchers in this line of research. Note that this is a cooperative research program between the Department of Mechanical Engineering at Michigan State University and the Department of Mechanical Engineering at the University of Michigan. Powell’s Hybrid method, as implemented via the NAG (Numerical Algorithms Group) routines, is available at the University of Michigan, and was used in the studies by Pesheck et al. ([33], [34], [71], and [82]), and in the study by Apiwattanalunggarn et al. ([35]). However the NAG routines are not available at Michigan State University, and therefore the Newton-Raphson method was implemented via the Numerical-Recipes routines([84]), and the Secant method was developed in-house and used in other parts of this study. Powell’s Hybrid method is a method used to find a minimum solution of an objective function, which in our case it is a summation of nonlinear-algebraic equations squared, whose global minimum point is the best available solution of the system of nonlinear- algebraic equations. The method requires no information on derivatives of the objective function. However it maintains a set of independent directions, where the number of directions is equal to the number of dimensions of the solution, and it performs successive line searches along the set in a cyclical manner. The set is also updated at each iteration. As mentioned in [84] and [83], for a problem with N unknowns, the quadratic-objective function would take N iterations of the basic procedure, which results in N 2 + 0(N) exact line searches, i.e., the minima is obtained along each line search. This estimate can be used as a guide for a general function since locally it behaves like a quadratic function at an arbitrary point. As mentioned in [82], it takes on average about 1.25N iterations to 51 reach a solution to equations (2.41) and (2.42) (with the D’s expressed as functions of the C’s), for either the global domain approach or the local domain approach (for one single strip) using this approach. Also, as mentioned in [82], the computational effort of the local (1 — 45 approach using Powell’s Hybrid method is only £2926 of the global a - (I) approach using Powell’s Hybrid method. This is because for large :mplitude motions many terms of polynomial basis (Na) are needed to describe the manifold globally along the ak direction for the global (L — (p approach. This causes the number of iterations to reach a solution to equations (2.41) and (2.42) to increase as Na increases, and the number of evaluation points needed for the integrands increases as well. However this is not the case for the local a — (I) approach since Na is fixed at 2, the number of evaluation points needed for the integrands is fixed, and the total computational effort for the whole domain is just the sum of the computational efforts for each annular strip. The Newton-Raphson method uses gradient information (the J acobian matrix) of the sys- tem of nonlinear-algebraic equations, evaluated at a current step, to determine the direction of the subsequent step, which advances the process towards the solution. The gradient in— formation of the Newton-Raphson method is determined in closed form from the system of nonlinear-algebraic equations. For a single nonlinear equation, the order of convergence ([80]) of the method is determined from lim M = ' w = c, where 7‘ is the n—+oo [enlp 71—400 [Tn — TIP th iteration, rn+1 is the approximate exact solution, Tn is the approximate solution at the 71 solution at the (n + 1)th iteration, c is a constant called the asymptotic error (c 75 0 ), and p is a constant called the order of convergence, where p 2 1. For a system of nonlinear equations, similar definitions can be formed using norms instead of absolute values ([85]). For the Newton-Raphson method, the order of convergence is 2 (quadratic convergence) ([80], [85]). There is no such estimate for the number of iterations needed in order to reach the solution from an initial guess, since in practice the solution is not known beforehand, so the error at each iteration is not known. Also, the order of convergence and the asymptotic error constant hold only when the iterates are close to the solution. The Secant method is similar to the Newton-Raphson method. However, the gradient 52 information is determined approximately using finite difference methods. The order of convergence of the Secant method is 1.6 (between linear and quadratic convergence) ([80], [85]). Just as for the Newton-Raphson method, there is no estimate for the number of iterations needed in order to reach the solution from an initial guess. Among these three methods, Powell’s Hybrid method has an advantage over the Newton- Raphson and the Secant methods. Since both N ewton-Raphson and Secant methods require the Jacobian matrix, and the cost of computing the J acobian grows rapidly as the dimension of the problem increases. The Newton-Raphson method would have an advantage over the Secant method since it generally requires fewer steps to reach the solution, due to the order of convergence ([80]). However, when the system of nonlinear-algebraic equations is complicated, the Secant method is efficient for computing the Jacobian matrix, since each entry of the J acobian cannot be determined in closed form, or is complicated and therefore it takes significant time to compute. In the next section, direct comparisons of the accuracy of the four approaches, both global-demain and local—domain, each using modal displacement/velocity and modal amplitude / phase, are made by comparing time simulations from two example systems. 2.4 Examples A two-DOF nonlinear spring-mass system and a finite element (FE) model of a rotating beam are used as examples to demonstrate the accuracy of the four approaches used to solve the invariant manifold equations. The first example is used as a “proof of concept,” while the second demonstrates the utility of the approach for systems with several DOF. For simplicity in identifying the various approaches, they are labeled as follows: global (L — (6, global u — 1), local (1 — (f), and local 11 — v, as described in sections 2.3.2.1, 2.3.2.2, 2.3.2.3, and 2.3.2.4, respectively. For the two—DOF nonlinear spring-mass system, all four 53 methods are considered. For the rotating beam, only the local (I — a) and local 11 — v approaches are considered. 2.4.1 A Two-DOF Nonlinear Spring-Mass System The two-DOF nonlinear spring-mass system considered here is the same as the system studied in [33]. The system has nonlinear springs, described by linear spring stiffnesses In and 163 and nonlinear spring parameters kg and k4; it is schematically depicted in Figure 2.1. The system parameters are m1 = n12 2: 1, and k1 = 1,162 = 2,k3 = 5, and k4 = 1. A comparison of a — gb and u — v global methods is considered first, followed by a comparison of local methods in a -— (f) and u — v. 2.4.1.1 Comparison of Global Methods The accuracy of the global methods is studied through a comparison of motions taking place on the NNM manifold of the first mode of the system. This manifold can be depicted as a pair of constraint surfaces which depend on either the first modal amplitude-phase or the first modal position—velocity. The two constraint surfaces restrict the second mode displacement and velocity, respectively, and describe motions that take place on the 2- dimensional NNM manifold in the 4-dimensional state space. The C coefficients describing the manifold solution in equation (2.36), with the D’s con- densed out, were solved using Powell’s Hybrid method implemented in the C language on a work-station computer (this was carried out by Mr. Dongying J iang, a graduate student at the University of Michigan, and a collaborator on this project). A sample of such surfaces is shown in Figure 2.2, which corresponds to the contribution from the second linear mode displacement to the NNM manifold. The boundary of all constraint surfaces in terms of modal amplitude is given by a0 : 2.6. 54 These surfaces are obtained by using Na = 5 polynomials and N¢ = 12 harmonics, which are the same as used in [33]. Note that Na 2 5 and Nd) = 12 result in 60 coefficients for each state when the conservative (symmetry) condition is used. Here the range of validity of the surfaces is based on comparing simulations of the original model and the N NM ROM. If the motion is initiated beyond a certain amplitude of the master mode, the time responses of the original model begin to have high frequency components. We define this amplitude to be the boundary of the valid domain. As mentioned in [33], the surfaces are valid up to a = 1.5 with (10 = 2.22. Here with (10 = 2.6, the surfaces are valid up to a = 1.7. Therefore, we suspect that an instability of the first NNM manifold may exist at near this amplitude limit. This suspicion is based on the fact that instabilities of NNMS were observed in 2-DOF systems with symmetry by Rosenberg et al. ( [6], [15], and [19]). The C coefficients describing the manifold solution in equation (2.43), with the D’s con- densed out, are solved using the secant method implemented with the symbolic processor TA! Mathematica on a personal computer. To yield a comparable number of C’s for this case, the number of polynomials along uk and 1),, would need to be NW, 2 11 and Np,” = 11. TM , and are therefore computa- Here the manifold-solver codes are written in Mathematica tionally slower. Thus, for acceptable computational time we used Npm = 8 and NW, = 8. Note that Np,“ = 8 and Np,” = 8 result in 15 coefficients for each state when the con- servative and odd-polynomial conditions are imposed. The C coefficients are solved with three different sets of boundaries, given by: (Ub, Vb) = (15,15), (Ub, Vb) = (1.7,1.7 X an), and (Ub, Vb) = (26,26 x wl), where w] = 0.69 rad/sec is the linear natural frequency of the first mode. The contributions from the second linear mode position to the first NN M manifold for the three different sets of boundaries are depicted in Figures 2.3, 2.5, and 2.7, respectively. Using the transformations given in equations (2.13) and (2.14), the previous contributions can be expressed as functions of modal amplitude and phase, as depicted in Figures 2.4, 2.6, and 2.8. It can be seen from Figure 2.5 that the surface bends sharply at around 11 = 1.3. Also, from Figure 2.8, which is solved for the larger amplitude, the solution surface shape looks completely different from the lower amplitude surfaces. These two observations are a sign of some type of instability, either dynamic or numerical, in the 55 first N N M manifold or its solution. Therefore, in order to compare the accuracy of the global a. — d) and global u. — v approaches, the surface in Figure 2.4 is used to generate the time responses. Figure 2.9 depicts the time responses of the first nonlinear mode displacement 771 obtained from the ROM of the global u —- 22 approach, the ROM of the global a — 62 approach, and the original model. The dotted line is the time response of the ROM defined at the end of section 2.3.2.2 with initial modal displacement 771(0) 2 1.3 and velocity 771(0) = 0.0. The dashed line is the time response of the ROM defined at the end of section 2.3.2.1 with initial modal amplitude (11(0) = 1.3 and phase angle 631(0) 2 0.0. The solid line is the time response of the original model (2 DOFS) initiated quite precisely on the first-NNM manifold (as determined by a shooting algorithm that finds periodic responses). As these responses all appear to be quite close, an examination of the slaved modes is used to offer a more refined look at the methods. Figure 2.10 depicts the time responses of the second linear mode displacement obtained from the ROM of the global u - v approach, the ROM of the global a — 4) method, and the original model. The dotted line is obtained by applying the master—slave constraint (depicted in Figure 2.3) to the numerical solution of the ROM of the global u — v approach (the dotted line in Figure 2.9). The dashed line is obtained by applying the master-slave constraint (depicted in Figure 2.2) to the numerical solution of the ROM of the global (L — (b approach (the dashed line in Figure 2.9). The solid line is obtained from the 2-DOF original model initiated on the NNM manifold. Figure 2.11 depicts the time responses of the physical displacement 2:1 from Figure 2.1 for the ROM of the global u — 22 approach, the ROM of the global a - qb approach, and the original model. The line types definitions are the same as in the previous figure. Horn Figures 2.9, 2.10, and 2.11, it can be seen that the time responses from the global a — 63 approach are closer to the original model than the global u — v approach. However, 56 this is not surprising, since the number of basis functions used for the global (L — ab approach was twice the number of basis functions used for the global u — v approach. 2.4.1.2 Comparison of Local Methods Similarly, the accuracy of the two local methods is studied through a comparison of motions taking place on the first NNM manifold. For the local 0. — (.75 method the C coefficients describing the manifold solution in equation (2.52), with the D’s condensed out, were solved using Powell’s Hybrid method implemented with the C language on a computer work-station (again, by Mr. D. Jiang). A sample of such a surface is shown in Figure 2.12, which depicts the contribution from the second linear mode displacement to the first NNM manifold. This N N M manifold is obtained over an amplitude range of a E [0, 2.6], which is generated using 52 piecewise linear manifold segments in a of width Aa = 0.05, and N4, = 12 harmonics in ¢- Note that Na = 2 and N9), = 12 result in 24 coefficients for each state within each strip when the conservative symmetry condition is used. The total number of coefficients for the local a — 65 approach is (Q — 1) x 2N¢ x Nsm‘p = 1248. Here the surfaces are valid up to a = 1.7. For the local 11 - v method the C coefficients describing the manifold solution in equation (2.62), with the D’s condensed out, are solved using the secant method implemented with the C language on a personal computer. The contribution from the second linear mode displacement to the first NNM manifold is depicted in Figure 2.13. The boundaries of all surfaces along the modal position and velocity are Ub = 2.6 and Vb = 1.79, respectively. All surfaces are obtained by using Chebyshev polynomials along 11; and v]: with Np," = 3 and NW) = 2, respectively, and the number of pieces along 11;, and vk to be N5 = 100 and N5 = 100, respectively. The total number of coefficients for the local 11 — v approach is I (Q _ 1) X 2N1),uNp,v X ENSIVS = 30,000 . 57 Figure 2.14 depicts the time responses of the first nonlinear mode displacement m obtained from the ROM of the local 11 — v approach, the ROM of the local a - (15 approach, and the original model. The dotted line is the time response of the ROM defined at the end of section 2.3.2.4 with initial modal displacement 711(0) 2 1.5 and velocity 1')1(0) = 0.0. The dashed line is the time response of the ROM defined at the end of section 2.3.2.3 with initial modal amplitude (11(0) = 1.5 and phase angle (61(0) = 0.0. The solid line is the time response of the original model (2 DOFS) initiated quite precisely on the first-NNM manifold (again, found using a shooting algorithm). Again, since these responses all appear to be quite close, an examination of the slaved modes is in order. Figure 2.15 depicts the time responses of the second linear mode displacement obtained from the ROM of the local 21 — v approach, the ROM of the local a — 45 approach, and the original model, using the previously defined line types. Figure 2.16 depicts the time responses of the physical displacement 2:1 from Figure 2.1 for the ROM of the local 11 — v approach, the ROM of the local (I — (15 approach, and the original model. From Figures 2.14, 2.15, and 2.16, it can be seen that the time responses from the local a — 45 approach are generally closer to the original model than those from the local 11 — v approach. Hence, for this example, and considering the number of coefficients used for each method, we conclude that the local a — qb approach provides more accurate manifold solutions than the local 11 — v approach. 2.4.2 A Finite-Element-Based Model of a Rotating Beam The finite element (FE) model of a rotating beam used here is the same as that described in [35], and analyzed in detail in the next chapter. The system consists of a uniform, unloaded Euler-Bernoulli beam attached to a hub that is rotating at a constant rate. The beam is restricted to vibrate in a plane that is parallel to the axis of rotation and rotating 58 with the hub, so that only transverse vibrations in one direction and axial vibrations are allowed. The corresponding transverse vibrations are commonly referred to as “flapping” motions, while transverse vibrations out of this plane are known as “lead-lag” motions. In this model, torsion, lead-lag, and other motions are neglected. The system is schematically depicted in Figure 2.17. The rotating beam parameters are L = 9 m, m = 10 kg/m, E1 = 3.99 x 105 N -m2, EA = 2.23 x 108 N, o = 30 rad/s, and h = 0.5 m. The objective here is to compare the accuracy of the methods for generating NNM manifolds (this model is discussed in more detail in the next chapter). The beam is modeled using 182 elements, which are used to construct linear modes and associated nonlinear coupling terms. This model is truncated to Q = 21 DOF, expressed in linear modal coordinates, and this is used as a starting point for developing the N N M ROM. This model has 9 transverse modes and 12 axial modes, selected according to the information described in the next chapter. For this system only the local methods in a — d) and u —- v are considered. As mentioned in section 2.3.3.4 the computational effort for the local a — d) approach with Powell’s Hybrid method implemented is only %,6—§% of the global a -— 65 approach with Powell’s Hybrid method implemented (as estimated by Pesheck in his dissertation [82]). Based on this estimate, it is reasonable to assume that the local approach is generally much faster than the global approach. Therefore the global approaches are not considered in this example. 2.4.2.1 Comparison of Local Methods The accuracy of the local methods is studied by direct comparison of motions taking place on the computed first N NM manifold. In collaboration with Mr. D. Jiang, the C coefficients describing the manifold solutions in equation (2.52) (with the D’s condensed out) are obtained using Newton’s method, implemented via the Numerical-Recipes routines ([84]) with the C language on a personal 59 computer. Note that for a single NN M there are 40 slave functions for this problem, that is, 40 constraint surfaces. Samples of these surfaces obtained by the local a — 66 approach are shown in Figures 2.18, 2.19, 2.20, and 2.21. These correspond to the contributions from the second and fifth linear flapping mode displacements and the first and sixth linear axial mode displacements to the first NNM manifold, respectively. The first NN M manifold is obtained over an amplitude range of a 6 [0.6.0], which is generated using 120 piecewise linear manifold segments in a with width AG. 2 0.05, and M), = 16 harmonics in 4). Note that Na = 2 and Ne = 16 result in 16 coefficients for each state within one strip, when the conservative symmetry condition is used. Thus, there are a total of 120 problems, each with 16 unknowns, that must be solved, and it took about 12 hours to solve for the manifold solution using a 2.0 GHz personal computer. We will consider simulations of this model after describing the solution obtained using the local 11 — v method. The C coefficients describing the manifold solution in equation (2.62) (again, with the D’s condensed out) were next solved using the secant method implemented with the C language on a 2.0 GHz personal computer. To yield a comparable number of the C ’s over the whole domain, the number of segments along 11;, and vk are taken to be approximately N5 = 40 and N5 = 40 with NW, = 3 and Np, = 2 to cover 11 E [-6.0,6.0] and v E [—6.0 x w1,6.0 x 621], where 6.21 = 34.03 rad/sec is the linear natural frequency of the first flapping mode. It takes about 40 minutes for the manifold solution with N5 = 40 and N5 = 40 to converge. However, the residuals of the manifold-governing equations are large. To obtain the appropriate boundary Uh. we fix Ub to be 6.0 and then increase the number of patches, for instance, N5 = 100 and N5 = 100. Then we simply observe the residuals for pieces along the 1)), direction for a few of the first pieces along the uk direction. From these observations, we can determine the range in the vk direction that gives acceptable residuals. This allows one to compute the appropriate Vb and the appropriate Ub. We then choose the same N5 and NS for these reduced values of Ub and Vb, i.e., we choose N5 = 100 and N5 = 100. Therefore, we take N5 = 100 and N5 = 100 to cover the smaller domain, 11 E [—4.0,4.0] and v E [-4.0 x 621,40 x 1.11]. This procedure yields small residuals over the entire domain. In this manner it takes about 4 hours to obtain the manifold solution. 60 Samples of constraint surfaces obtained by the local 11 — v approach are shown in Figures 2.22, 2.23, 2.24, and 2.25, which correspond to the contributions from the second and fifth linear flapping mode displacements and the first and sixth linear axial mode displacements to the first NNM manifold, respectively. Figure 2.26 depicts the time responses of the first nonlinear mode displacement 771 obtained from the ROM of the local 11 — v approach, the ROM of the local (1 — cf) approach, and the original model with the full 21 DOF. The dotted line is the time response of the ROM of the local 11 —- v approach with initial modal displacement 771(0) = 2.375 and velocity 171(0) = 0.0. The dashed line is the time response of the ROM of the local (1 — 63 approach with initial modal amplitude (11(0) = 2.375 and phase angle 961(0) 2 0.0. The solid line is the time response of the original model (21 DOF) initiated on first NNM manifold as obtained from a shooting algorithm. Figures 2.27, 2.28, 2.29, and 2.30 depict the time responses of the second and fifth linear flapping mode displacements and the first and sixth linear axial mode displacements obtained from the ROM of the local 11 — v approach, the ROM of the local (1 — (25 approach, and the original model, respectively. Figures 2.31 and 2.32 depict the time responses of the beam-tip flapping deflection and the time responses of the beam—tip axial deflection, respectively. From the above figures, it is obvious that the time responses from the local (1 — 65 approach are closer to the original model than the local u—v approach. We conclude that the local (1- ¢ approach yields more accurate manifold solutions than the local 11 —— v approach. However the local 11 — v approach yields the manifold solutions with significantly less computational time. Therefore, the local 11. — v approach can be used as a preliminary study or for a study in which very high accuracy is not required. The drawbacks of using 11 —— v coordinates might be eliminated by first applying the invertible van der Pol transformation ([22]) to the master u -— v coordinates. This would results in a new rectangular coordinate system that rotates with the linear response, which would uncouple the governing equations of the master states up to linear order. In 61 other words, we would parameterize the invariant manifold using a dynamically motivated rotating-coordinate system. Once the invariant-manifold governing equations are reformu- lated, then the patch approach can be applied again with low-degree polynomials. Note that this idea has yet to be tested. 2.5 Conclusions In this chapter, we have proposed two new methods for solving for individual NN M in— variant manifolds. For both of these methods the manifold is parameterized by the master displacement and velocity and the manifold is obtained using Galerkin approaches. For the first method, the manifold is solved using global polynomial basis functions over the domain of interest. For the second method, the domain of interest is subdivided into small pieces and the manifold is solved using the collocation method over these sub domains, and then pieced together. These methods are compared, in terms of computational times and accuracy, to similar methods that employ amplitude and phase master coordinates, as developed by Pesheck et al. ([33]). In terms of computational time, the proposed methods have advantages over the. methods of Pesheck et al.. In terms of accuracy, the methods of Pesheck et al. are superior to the methods proposed herein. In general, a combination of these approaches, which uses rotating rectangular coordinates, as described above, may offer good efficiency and accuracy. 62 2.6 Tables Approaches No. of Coefficients No. of Eval. Points Consv. Syst. Consv.-Cubic-NL Syst. per integrand Section 2.3.2.1 (Q — 1) (Q — 1) (3N0 + 3) Global (1 — Q X2NaN¢ XZNQNQ, x10N¢ Section 2.3.2.2 (Q — 1) (Q - 1) (57%,, — 2) Global 11 — 6 x (Np,,,Np,., — 3) x §(NP,.,NP,, - 3) x (336,, — 2) Table 2.1. Comparison of the total number of unknown coefficients and the number of evaluation points needed for the integrands in order to compute the integrals for the global domain approaches. Approaches No. of Coefficients for Consv. Syst. No. of Eval. Points per Integrand Section 2.3.2.3 (Q —— 1) 1 Local (1 — ¢ x2N¢, x Ngmp Section 2.3.2.4 (Q - 1) 1 Local 11 — v x2NEuM x iNfiNg Table 2.2. Comparison of the total number of unknown coefficients and the number of evaluation points needed for the integrands in order to compute the integrals for the local domain approaches. 63 2.7 Figures X1 X2 —> —> ; k3’ k4 /, ,/ m, m, 2; s\\ \\\\\ x \ \ Figure 2.1. A two-DOF nonlinear spring-mass system. ‘5 c ”. '0 Q o 0 0. a o \ o c 3 1’... 9:44'0, ’l’ Figure 2.2. The contribution of the second linear mode position 712 = P2(a1, Q1) to the first nonlinear mode manifold obtained by the global a — cf) approach. 64 : ~ s~§~ - . s ¢ :V J a. a... .- - - - Q s Figure 2.3. The contribution of the second linear mode position in = X2(n1,I'71) to the 1.5 and first nonlinear mode manifold obtained by the global u — v approach with Ub V6 1.5. 3.952.954? s The contribution of the second linear mode position 772 X2(171(a1, $1), 7')1(a1, (91)) to the first nonlinear mode manifold obtained by the global u - v approach with Ub = 1.5 and Vb = 1.5. Figure 2.4. 65 ~ . ~ 4 w . ¢~¢-~- . t .3 ow . .041.“ v.6 . . c .1 ~“QQJQQO‘. . :3 cf \ h ~0§§\\ .4 y l n. ..:...... . . 0 . . - . . ~ ~ ~ ~ . . ~ ..........:.s~ . 3 -~4 33ft: Figure 2.5. The contribution of the second linear mode position 112 = X2(n1,1'71) to the 1.7 and first nonlinear mode manifold obtained by the global u — v approach with U), V, =1.17. agiééeééax n 3 linear mode position 1n = second The contribution of the X2(n1(a1,¢1),1)1(a1,¢1)) to the first nonlinear mode manifold obtained by the global u— v approach with Ub = 1.7 and Vb Figure 2.6. 1.17. 66 o ~ ~ ~ ~ .-...~... .1... I 1...: .1 I -~ . Jo-~- -Q 42... . 1 .3... .~ 0 O Q O O r - QOQa. A 99 .fi .0 O A X2(171,7'71) t0 the Figure 2.7. The contribution of the second linear mode position 712 2.6 and first nonlinear mode manifold obtained by the global u — v approach with Ub Vb 1.79. O _ 3.96.295?” a mode position 02 linear second The contribution of the Figure 2.8. 1.79. 2.6 and Vb = X2(n1(a1, 1151), 1')1(a1, ¢1)) to the first nonlinear mode manifold obtained by the global u — v approach with Ub 67 2 GloLbal u—v. Nth. 1 DOF [ — — Global a—(p. NNM. 1 DOF -~— Original Model, 2 DOF O N)- A "0,, co .1 O .1 N Figure 2.9. The time responses of the first nonlinear mode displacement for the ROM of the global u — v approach with 771(0) = 1.3171(0) = 0.0, the ROM of the global a — 12 approach with 111(0) = 1.3, 651(0) 2 0.0, and the original model with initial conditions from the shooting algorithm. 0.2 ~ ' Gldbal u—v, Nth, 1 DOF _ — - Global a—¢. NNM, 1 DOF —~ Original Model. 2 DOF 0.15" 0.1 r 0.05[ [ :N0 -0.05 * —O.1 r —0.15* 1 Figure 2.10. The time responses of the second linear mode displacel'nent 772 on the first nonlinear mode manifold for the ROM of the global u — v approach, the ROM of the global a - (25 approach, and the original model. 68 Gldbal u-v, NfiM, 1 DOF — - Global a—Q, NNM, 1 DOF 1 )- ———— Original Model, 2 DOF H bx 0.5— \ x" 0i —0.5 r ..1 i- " O 2 4 6 8 10 12 Figure 2.11. The time responses of the displacement 51:1 from Figure 2.1 for the ROM of the global u — v approach, the ROM of the global a — Q approach, and the original model. Figure 2.12. The contribution of the second linear mode position 772 = P2(a1,Q1) to the first nonlinear mode manifold obtained by the local a — Q approach. 69 xx .x xexx xxxxxx xx xx xx xx xxx xx . xx x xx x n xxx xx xxx\\\xxx\ x H xx xxx xxx xx x xx xx xx xx x xx xx xx xxxxxx xx xx xx xx x xxxxxxx \ x x x x x x. 8.. xxx xxxx \ xx x x x xx xxxxxx xxxxxxxxxx x xxxxxxxxxxxxxx xx xx xxx x xx xxxxxxx xxxxx xx x x x xxxxxx xxxx xxxxxx x x x x x x x xx xxxx . xxxxxxxx xx xx xx xxxxxxx. x x x x.x xxxxxxxxxx xx x x x x x x x x x xxxxxx..x xxxxx x x xxxx xxx x xxxxxxxx x x x x xxxxxxx x x \ ‘ x x x‘zx“¢‘ x‘ 8‘ ‘ K ‘\ x x x“ x ‘x x‘. lilim. 4202 . . . 00 00 __ 3.5;" 9 xx x x x xxx xxx xxx xxx xxx xx x xx xx x x x x x x x xx xxx xxx xx xxxxxx x x xx xx xx x \\ xx xx xx xx xxxxxxxx xxxxxxxxxx xxxxxx xxxxx x xx xx xx .xx» xxx xx x. x x x. xx xx xx xxxxxxx xx xx xx xx xx xx xxxxxxxxx xx x x x x xx xx .x .x i x x x xxxxxxxxx xxxxx xxx xxx xx Figure 2.13. The contribution of the second linear mode position n2 = X2(n1,1'71) to the first nonlinear mode manifold obtained by the local 11 — v approach. Local u—v. NNM, 1 DOF Local a—Q. NNM, 1 DOF — Original Model. 2 DOF l 2.. 1.5“ 12 10 t Figure 2.14. The time responses of the first nonlinear mode displacement for the ROM the ROM of the local 11 — Q 1.5, Q1(0) = 0.0, and the original model with initial conditions from 7 0.0 ) O = 1.5,7)1( 0) ( of the local 11 — v approach with 171 approach with a1(0) the shooting algorithm. 70 0.3% I ~ Lgal u—v, NNrM, 1 DOF l. — — Local a—Q, NNM, 1 DOF ——— Original Model. 2 DOF [ Figure 2.15. The time responses of the second linear mode displacement W on the first nonlinear mode manifold for the ROM of the local 11 — v approach, the ROM of the local a — Q approach, and the original model. Local u-v, NfiM, 1 DOF _- ~ Local a-o, NNM, 1 DOF 1 _ —— Original Model. 2 DOF Figure 2.16. The time responses of the displacement 3:1 from Figure 2.1 for the ROM of the local 11 — v approach, the ROM of the local (1 — Q approach, and the original model. 71 We(xe: 1) ®e(xe,t) Ue(xe.t) e=1. 2. 3...,, M l l l l l e1 I | | l l l J x9 xe+1 \E, A, I, L, m Figure 2.17. Finite element representation of a rotating beam with Q = constant. P2(a1,¢1) (2"d Flapping Mode) Figure 2.18. The contribution of the second linear flapping mode ”()2 = P2(a1,Q1) to the first nonlinear mode manifold obtained by the local a — Q approach. 72 ) (5th Flapping Mode) .3: 8 P5(a|1’¢1 .0 03 Figure 2.19. The contribution of the fifth linear flapping mode 712 = P5(a1,Q1) to the first nonlinear mode manifold obtained by the local (1 — Q approach. O '_. O 51 . P1o(a1.¢1) (1 Ax1a| Mode) Figure 2.20. The contribution of the first linear axial mode 1710 = P10(a1,Q1) to the first nonlinear mode manifold obtained by the local a — Q approach. 73 $8.2 53. 52 333 m— n. la] mode 7715 = P15(a1,Q1) t0 the first lnear ax Figure 2.21. The contribution of the sixth l nonlinear mode manifold obtained by the local (1 — Q approach. . x xxx . x 1 xxxxxxxxx x xxx xxx xx x ..xxxx Xxx .o xxx xxx x xx . x .x xx .xx. xx x .. xxxxxxexxxxxxxxxxxxxxxxxx. .xxx. waxes......s.x........................ ii. xxxxxxxxxxxxxxxxxxxx xxx. xx x xxx. xx x x xxx. x x xx xxxx xx x x x x x. xx xx. x xx ...x .11.. 1.2.... xx x xx x x xx x x xxxxxxxxxx x xx xxnxxxxxxxxx “xx 4......,......,.....N........a..n.....a. . \ x 9 x x s x x x Q x x x x x x . x x x xxxx xxxx xx xxxx xx x xxxx x x xxxxxxxx xxxxxxx xxxx xxxx x xxx xxxxx xx ..xx x. xxxx xxcxxxxxxxxxxxxx x xx x x x x xx x x xxxxxxx x x x x Figure 2.22. The contribution of the second linear flapping mode 712 = X2(n1,r']1) to the first nonlinear mode manifold obtained by the local 11 — v approach. 74 X5(T}1,fil) t0 the first x xxxxxxxx xxxx xx x x \\\\ x xx xxx xxxxx xx x x x x x x\xxxx\\\\ \\ .xxxxx xxxx xx x x «xxx \i\\~\\ \ \\ xx xxxxx x xx xxx x x \ xx xxx\\ x xx xx xx x x x xxx xx \ K x xxxxxxxx xx x x xx \ \ \ . x x x x x xx xx xx .xxxx xxxxx xxx xx xxxx \ . . x . x x x x x x x x xx...xx xxxxxx x . x x . x x x x x x x xx xx xx x x x xxxxxxxxxxxxxxxx xx ...x... .x xxx xx .xxxxxxx . . ..xxxxx xx xxxxxx xxxxxxx x x x xx x x \ x xx xxxx . x xx . xx. xxxxx xxxx xx xxxxx x . x xxxx. x xxxx xxxx xx.xxxx xx xx x x x xx .xxxxx. xxxx xxxxxxxxxxxxxxxxxnxxxxx . x.x xxxxx xxx xxxxx. xxxxxxxxxxxxxxxxxxx x x x xxxx. xxx. xxxx xxxxxxx xxxx xxxx x x xx x xxxxxx xxxxxx. x x x x x . x x x x xxxxxxxxxx xxxxxxxxxx x x xxxxxxxxxxx xxx xx xxx xxxxxx x ion of the fifth linear flapping mode 712 nonlinear mode manifold obtained by the local u — 11 approach. Figure 2.23. The contribut ial mode 7710 = X10(n1,i)1) to the first nonlinear mode manifold obtained by the local u -— 11 approach. inear ax 75 ion of the first 1' but 1 Figure 2.24. The contr \\ \ \ \ \ \ \ \ x \ x x x“ \x\\\\\\\\\x\ \\x\ xx\\xxxx\x xx ‘ ‘ x ‘x\‘ \‘ x‘ x“x~\ x x \\\‘x \‘x‘x‘ x \ \\ \\\\\\ Figure 2.25. The contribution of the sixth linear axial mode m5 = X 150113771) to the first nonlinear mode manifold obtained by the local u — 11 approach. - —- Local 8—4). NNM. 1 DOF Y Local u—v. NNM. 1 DOF 3' — Orlginal Model, 21 DOF .1 r Mode 1 (1St Flapping Mode, Master Mode) Deflection O 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 t Figure 2.26. The time responses of the first nonlinear flapping mode displacement n1 for the ROM of the local u — v approach with n1(()) 2 2.375,7‘]1(0) = 0.0, the ROM of the local (1 — ¢ approach with (11(0) = 2.375, (151(0) = 0.0, and the original model with initial conditions from the shooting algorithm. 76 ~ Local u-v, NNM, 1 DOF 0.15% ~ Local a—o, NNM, 1 DOF - : —— Original Model. 21 DOF 0 § 0.1 l‘ \ "“5 o 0.05 R, a) \ é o~ / . a K \ .s / Q-O.O5 ~ 3 / LI. ‘2 -O 1 r ' N —0.1 5 ~ j 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Figure 2.27. The time responses of the second linear flapping mode displacement in on the first nonlinear mode manifold for the ROM of the local u —— v approach, the ROM of the local a — Q5 approach, and the original model. —3 6 x 10 Local u—v. NNM. 1 DOF Local a-o. NNM. 1 DOF 4 _ Original Model. 21 DOF 2 h i ".1va.. 5th Flapping Mode Deflection 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Figure 2.28. The time responses of the fifth linear flapping mode displacement 775 on the first nonlinear mode manifold for the ROM of the local u — v approach, the ROM of the local a — qb approach, and the original model. 77 Local u-v. NiuM, 1 06F 0.01 i — - Local a-¢, NNM, 1 DOF F, ~-— Original Model, 21 DOF O'- /\ / fl / —o.02L ’ —o.03 —o.04~ 4.057 —o.06~ \ ‘ —0.07 4 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 t T I 1St Axial Mode Deflection Figure 2.29. The time responses of the first linear axial mode displacement 7710 on the first nonlinear mode manifold for the ROM of the local u — 11 approach, the ROM of the local or — 96 approach, and the original model. x 10"4 9 . Local u-v. NNM. 1 DOF —- - Local a-c. NNM, 1 DOF 3. , .. Original Model. 21 DOF 7 .8 7 3 6 8 5 § 4 2 E 3 5 2 (D 1 0 l ..A 0.25 0.3 0.35 o o o 01 .0. A o d m o m Figure 2.30. The time responses of the sixth linear axial mode displacement 7715 on the first nonlinear mode manifold for the ROM of the local u — 11 approach, the ROM of the local a — qb approach, and the original model. 78 0.8 fl 1 T I I I f , Local u—v, NNM. 1 DOF — — Local a-o, NNM, 1 DOF E 0.6 r - — Original Model, 21 DOF g 0.4 ~ C 8 0.2 ~ :3 ‘55 o 0 r \ 8' L? -O.2 — - .9 l / —O.6 r ‘ ’0'80 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Figure 2.31. The time dependence of the beam-tip flap deflection, w(L, t), for the ROM of the local u — v approach, the ROM of the local a -— o approach, and the original model. Local u~—v. NNM. 1 06’: a n Local a—o, NNM. 1 DOF 4- —-— Original Model. 21 DOF O —2 -4 -6— —8 Tip-Axial Deflection (m) —10 ~12 -14 Figure 2.32. The time dependence of the beam-tip axial deflection, u(L, t), for the ROM of the local u — v approach, the ROM of the local (1 — ab approach, and the original model. 79 2.8 Appendix 2A The samples of {L1’lIL.,n(uk. 1.115)} for NW, 2 3 and NIH) = 3 in section 2.3.2.2 are as follows: Lilli,1(“ka“k) Lrl11.2(ur~~ vk) “11.30%: 'l’k) L1112.1(‘llk,llk) L1\[2,2(“k~1’kl LAI‘ZBULk, PA.) LAI3.1(“kavk) LM3.2(UA~»'U1) LA'13,3(UA~, We) 0 0 1.11803(Lf‘)2 LI) 0 _ ilk U}: 1.3 — —— “-1- ‘Uk 2 1.93649 — — (U1, )< vb) 1.55005(%/‘i)2 — 65.7441(%)10 +172.187(-:’/i)8 b b (J —163.124(%)6 + 67.1689(%~)4 —11.1948(;7k-)2 b b b we 9 Uk 7 ’Uk 5 17.4667 — — 40.4492 —— 32.1214 — — 9.99334(:’/—")3 + 2.88269(%5)2(:—5) b b ’b —1.47321(;‘j—"')2 + 62.4849(:)/—k)10 —163.651(%)8 12 b b +155.03s(5’f—‘)6 — 63.839(”—{“‘)"1 + 8.91533(%)2 + 5.17352( ‘ b Vb Vb 80 1152352 2.9 Appendix 2B Equations (2.48) and (2.49) in section 2.3.2.2 are solved using the Secant method. The method implemented is described here by assuming that the nonlinear term fz- in equation (2.2) is a function of the modal positions only. The orthogonality conditions from equation (2.47) are used in equations (2.48) and (2.49), the simplified equations are given respectively by 0 = Ff“ : _Dg’~q+ Z A'gq‘l“"‘(C)Cf“"' (2.67) (l+m)#2,3 0 = Fix“ = ..ch’~q—f§’*q(0)+ Z Kgq’l’m(C)D:’m, (2.68) (l+m);é2,3 , where Kgq‘l‘m, Kgq‘l’m, and fip‘q are given by r1).q.l.m _ ,p.q.l.m AC _ AD 1 1 81/11”). = / (—) (7) LAIDJ] “k TE “k’vk b b llk 6LAI, + (“-1.01% “k + fk) —()i—lr£:l (Ink dvk (2.69) ’k , l l '“kd’k L'b b Note that if either the conservative conditions (equations (2.50) and (2.51)) or both the 81 conservative conditions (equations (2.50) and (2.51)) and the odd-polynomial conditions are used, then Kgq‘l’m 76 Kgq'l’m. From equation (2.67), it can be seen that we can express the D’s as functions of the C’s. Those expressions are used in equation (2.68). Hence Fix’p’q will be functions of the C’s only and they are given by F.X"’“’(C) = 0 . (2.71) If C’“ is a solution of equation (2.71) in vector form, then equation (2.71) in vector form FX ’1’"? (C) can be approximated using the Taylor series expansion about 00)) , which is close to C", up to the first-order terms by 81?" 11M! 0 = FX’P’WC) z FX’qu(C(°)) + T(C(°))(C — 0(0)) . (2.72) anip’q 0 Defining the Jacobian matrix J E —8—C_(C( )), the approximate C* is given by 0* a: —J-1 vapv4(c(0)) + 0(0) . (2.73) If J is determined in closed form, then the method is called the N ewton-Raphson method. If J is determined approximately using the finite difference method and its components are given by 82 X,l.m thsm IMHO) ‘ _ X,l,m p,q,(0) J’~"W~Q=——05 (C(0))zFi (Ci “:2 F2- (Cj ), (2.74) where AC is a chosen small number, then the method is called the Secant method. For each component Ji’JI-nip’q, it is required to recompute Kgq‘l’m and K $71,137:. at each step due to the perturbed C's. Therefore. the computation of each component can be very expensive. The following algorithm is implemented instead: (i) Initialize n. = 0, tol 2 small number, norm > to], and 0“». (ii) \Nhile norm > tol (a) Compute Kg” = {Ix’gq‘l‘m}(”) and Kg) = {A’lqu’l‘m}(") using equation (2.69). (b) Compute FX’p’q(C(n)) using equation (2.71). (c) Compute J(") : J(")(C(n),K(Cn),Kgl),FX’p’q(C(n))) using equation (2.74). Note that here K gt) and K g) from step (a) are used instead of recom- puting them with the perturbed C ’s. (d) Compute 6C(") = —(J("))_1Fx’p’q(C(n)). (e) Compute = || 60“” ”2 (f) Compute C("+1) = C(") + 600‘). (g) Compute n = n + 1_ Similarly, equations (2.65) and (2.66) in section 2.3.2.4 are solved using the Secant method. Therefore, equations similar to equations (2.67) and (2.68) are given by 83 0 = FY‘M 2 .,Im D1771 (,I,(m = 2K6? +ZI<2P [.771 I m 0 = FXPPPQ 7 :. Z(_w 12)A,pr11”C.I771 1,777 1,771 (1 I 771 A4). .1,,qI m A.p,(1,I,-m A,]).(1.I.In ,rp $111,771 _ [pq./.772 I‘CJ _ ADJ : —LIl’[[’m (11:31” 17;”) A,p.q.l,m x1) (12. I. 771 C2 AD 8L1“ : (d1; + U z)_ I. 771 ). f: q = fll(u’. _le)~(qC BLIl-I 0 + a, éfll ”k bk C)CIP’” (2.75) 2 )+ Z Kp’ .7127 "’( C)Df.P'"‘, (2.76) 1,772 , and ff“) are given by (2.77) 63(2) +112.) + ml (2.78) J (ui2u‘z’p,vi=v£’q) (2.79) Note that the algorithm used to solve equations (2.75) and (2.76) can be developed in the same manner as the one described above. 84 CHAPTER 3 Finite-Element-Based Nonlinear Modal Reduction of a Rotating Beam with Large-Amplitude Motion 3.1 Introduction The dynamic analyses of helicopter blades, turbopropeller blades, wind-turbine blades and robotic arms has provided motivation for investigations of the vibration of rotating beams. To predict the dynamic characteristics of rotating flexible structures, the kinematics must be carefully modeled, which leads to nonlinear coupling effects between degrees of freedom (DOF) in different directions. These coupling effects can cause slow modal convergence, thereby often requiring large system models for accurate dynamic representation. Simula- tion of such large—scale models is a time consuming process, which slows parametric studies and design cycles. Much work has been done using finite elements (FE) to model the nonlinear, large amplitude vibrations of rotating beams, including [44], [45], [46], [47], [48], and [49]. These models are typically quite complicated due to their geometry. degrees of freedom (flap, lead-lag, axial, and torsion), and nonlinear coupling effects. Furthermore, because of the nature of the finite element approach, many elements are required in order to obtain an accurate model. 85 A common approach is to use linearization of the finite element model about the nonlinear static equilibrium position and solve the eigenvalue problem of the resulting linearized model to obtain the linear natural frequencies of the system ([44] and [49]). Bauchau and Hong ([45]) also utilized finite elements in time to obtain nonlinear responses and stability results of the rotating beam undergoing large deflections. However, the computational time associated with obtaining the equilibrium solution was expensive, because all the spatial degrees of freedom at all time steps are coupled. Perturbation modes ([50] and [51]) were applied to the finite element model of a helicopter rotor blade in order to obtain a reduced order model ([46]). Bauchau and Bottasso ([52]) applied perturbation modes to the space- time finite element model of a beam subjected to a sinusoidal load in order to obtain a reduced order model. Crespo da Silva ([53]) utilized a truncated set of eigenfunctions or eigenvectors obtained from the linearized system of partial differential equations (PDEs) or the linearized finite element model about the nonlinear static equilibrium position in order to obtain a reduced order model of a beam in planar motion. Crespo da Silva ([54]) also extended his work to handle multi-beam structures in planar motion. In most nonlinear structures, there is no simple expansion of basis vectors which decouples the DOF (i.e., modes) in the frequency range of interest from those outside that range. Therefore, some (potentially important) nonlinear effects will be ignored in the truncation process. Generally, many linear modes must be retained in the nonlinear model in order to minimize these effects. Over the past decade, systematic procedures have been developed to obtain reduced order models (ROM) via nonlinear normal modes (NN M) that are based on invariant manifolds in the state space of nonlinear systems ([20], [21], and [74]). These procedures initially used asymptotic series to approximate the geometry of the invariant manifold and have been used to study the nonlinear rotating Euler-Bernoulli Beam ([3]). More recent work has employed a numerically-based Galerkin approach to obtain the geometry of the N NM invariant manifolds out to large amplitudes ([33]). These procedures can be applied to more general nonlinearities over wider amplitude ranges, and have been recently applied to study the vibrations of a rotating Euler-Bernoulli beam ([34]). In that study, the system 86 of PDEs governing the axial and transverse motions of the nonlinear rotating beam are derived by Hamilton’s principle and discretized by a Rayleigh-Ritz method. These PDEs are similar to those in [86], which were derived by Newtonian methods. The study presented here uses the same PDE model considered in [34]. However, in this study we use Hamilton’s principle as a basis for deriving a nonlinear finite element model of the rotating beam. The primary goal of this effort is to demonstrate that the invariant manifold NN M approach can interface with nonlinear FE models, thereby Opening the door to the application of the approach to systems with more complex geometries and additional degrees of freedom (e.g., lead-lag and torsion). This chapter is outlined as follows. The nonlinear FE model is first generated and converted into a truncated (but still large-scale) modal form that is convenient for the NNM analysis. The invariant manifold equations are formulated, and a numerical collocation method is used to obtain the solution of the NN M invariant manifold for the fundamental flapping mode of the beam. This invariant manifold is used to construct a nonlinear single DOF ROM, which is subsequently used for a simulation study. The linear natural frequencies of the transverse motion, obtained by solving the eigenvalue problem of the linearized FE model, are verified against analytical solutions from [1] and [2]. A study is also carried out to determine the size of an appropriate reference model for the full nonlinear system, which is used for comparisons with the ROM. Detailed results for the fundamental nonlinear flapping mode form the bulk of the numerical results, including the amplitude-frequency relationship, the nature of the NNM invariant manifold, the associated master-slave modal relations, and the actual beam dynamics. In all cases, the results obtained from the ROM are checked against the full reference model, and against the collocation-based ROM determined in [34]. Some conclusions are drawn at the end of the chapter. 87 3.2 Rotating Blade Formulation A dynamic model for the vibrations of a uniform, cantilevered, Euler-Bernoulli beam, attached to a rigid rotating hub is considered. The system is schematically shown in Figure 3.1. In the formulation, the following assumptions are made: rotary inertia is neglected; the motion is restricted to axial and transverse vibrations in a rotating plane, i.e., twist and lead-lag motions are not considered; and, the nonlinear axial strain due to element extension is included *. A finite element approach is adopted, which is based on the element. kinetic energy, T, and the element potential energy, U which are expressed as follows, 1 “+1 -2 -2 2 2 T = —2- m(u. + w )+ 7720 (h + :L‘ + u.) (11: (3.1) Jig 1 ‘1'P+1 1 U = 2 / ‘ 131(wm)2 + 132107,, + §(w,3;)2)2d:r (3.2) 1'8 where u(;r, t) and w(;r, t) are the axial and transverse displacements, respectively, (-),x is a derivative with respect to the spatial variable :13, an overdot represents a time derivative, h is the hub radius, 9 is the constant angular velocity of the hub/ beam, m is the beam mass per unit length, E, A, I are the usual beam material and geometric parameters, and are and :66“ are the global coordinates at nodes 6 and e + 1, respectively, of a typical element. By applying Hamilton’s principle to these two expressions, the weak formulation for the equations of motion is obtained: ‘Note that other sources of nonlinearity, such as finite curvature effects, while potentially significant in magnitude, are not considered here, since our primary interest is the axial/transverse coupling. 88 t2 .1“€+1 .. 2 .. f / [—-mu + 7779 (h + :r +11)](5u + [—mw]6w + II fife [——EA(u,_.,; + %(Pu7,r)2)](6uvr + 16,751.31.) + [—EI‘lU,1"r]6u/I9Ix (1.7:(11: : 0. (3.3) A finite element model of the simplified rotorcraft blade is developed by introducing fol- lowing expansions into the weak formulation “gill t) R5 2 ¢1(I)ui(tl iue(:17,t) % Z lI’ihrlsift) (3-4) . 721 721 where ue and we are u and in restricted to the spatial domain of element 6, (352(2) are standard linear shape functions, 11,-(t) are nodal variables for axial displacements, 16,-(33) are standard cubic shape functions, and 3,-(t) are nodal variables for transverse displacements and rotation angles. By substituting these approximations into the weak formulation, one obtains the linear element mass and stiffness matrices, the force vector due to rotational effects, and the quadratic and cubic nonlinear force vectors, for both axial and transverse deflections. These are given by the following expressions: 89 ~71 PTe+1 . ,. 777,-]- : / mo,(a:)r3j(:r) dx, (3.5) ’l C , Ic+l ([0 (1G) 6:3 = < 7,— (,1 n Meme» dx (36) Te ' u re+1 2 ‘ f, = mil (I7+:1')(,9.,P(.r) (1.7 (3 7) 1‘e 4 4 9:1 : 2211]],157377» (3.8) [=1 71:1 ”38+1 1 (1(5- du”) d773,, Pl . I" = —EA—‘—£—'—'—— (in, . “ 1e1e all” [.8 2 (11‘ (II (1.1: I (3 9) . 1.6+] and 777.21 : / 777.72’)k(:r)7.5',(.r) (1.7:, (3.10) Te :7) 6’ +1 (1121,71. (i211!) I I78: = / E1 ' (1;, 3.11 A" 1'7: 611:2 (11:2 E ( ) 4 2 g“: = Z Z ailmslum, (3.12) [:1 77721 $e+l dwk (It?) d- 1.9777 VVllCI‘C (121," = A EAE—IZE— (11' 1, (3.13) --e . . 4 4 4 and hi. 2 ZZZbilmslsnsr, (3.14) [=1 n:1r=l xe+l 1 (1717;, d717, (1112,, (IL/’7 where bZ'Inr : [L -2- fi-E—d—T— dT (11‘, (315) e .. .. .. where the indices take on values as follows {7' = 1,2}, {j = 1,2}, {k = 1,2,3,4}, {I = 1,2,3,4}, {777 = 1,2}, {77, = 1,2,3,4}, and {r = 1,2,3,4}. (-)“‘ denotes quantities in the u direction. 777?]- is the ij component of the linear axial mass matrix, It]: is the ij component of the linear axial stiffness matrix, f,” is the 7' component of the axial force vector, y? is the 2' component of the quadratic nonlinear axial force vector, and ai‘m is the In quadratic nonlinear coefficient of 913‘. ()3 denotes quantities in the .9 direction. mil is the kl component of the linear transverse mass matrix, k2] is the kl component of the linear transverse stiffness matrix, 9;: is the k component of the quadratic nonlinear transverse force vector, (721m is the [777 quadratic nonlinear coefficient of 9):, I72, is the I: component of 90 the cubic nonlinear transverse force vector, and bfdm is the [777 cubic nonlinear coefficient of hi. M one-dimensional beam elements are used and assembled using standard procedures ([87]). Since there are three degrees of freedom at each node, and accounting for the boundary conditions, one obtains 3M nonlinear equations of motion in the nodal coordinates, which can be expressed as follows: MUU + KUU + GU(S) : Fa (3.16) M3s+KSS+CS(U, S) +HS(S) = o (3.17) where U is the global axial displacement vector and S contains the global transverse displacement W and rotation angle vectors 6). Quantities in the U direction are denoted (-)U, as follows: M U and K U are the linear axial global mass and stiffness matrices, respectively, and GU( S ) is the quadratic nonlinear axial global force vector. F9 is the axial global force vector due to rotational effects. Quantities in the 8' direction are denoted by (0)5 , as follows: M S and K S are the linear transverse global mass and stiffness matrices, respectively, and GS (U, S) and H S (S ) are the quadratic and cubic nonlinear transverse global force vectors, respectively. To include rotational effects in the transverse stiffness K S , thus accounting for centrifugal stiffening, the axial displacement U is separated into static and dynamic components, as follows: 91 where Us is the “static” axial blade deflection due to F9 and Ud is the dynamic axial displacement relative to U3. U3 is obtained by solving for the axial deformation in the case of no vibration in either the transverse or axial directions, i.e., KUUS = F9, (3.19) where GU(0) = O has been used. Figure 3.2 depicts the comparison between Us obtained by a 182-element model and an analytical solution ([3]), indicating very close agreement. (In fact, good agreement can be achieved with as few as five elements for this simple shape.) Substituting equation (3.18) into equations (3.16) and (3.17) offers a convenient form of the equations of motion. A linear transverse stiffness term in S, Kq(U3), which arises from rotational effects, can be separated out from G5 (U, S). This represents the linear transverse centrifugal stiffening due to rotation, which is combined with K S to provide the total transverse stiffness. Also, the rotational force vector is cancelled by using equation (3.19). The finite element model then takes the form: MUU‘d + KUUd + GU(S) = 0 (3.20) M53 + (KS + Kq)S + GS(Ud, S) + H5(S) = 0. (3.21) These equations have the form desired for the application of modal analysis. In particular, they have a zero solution, (Ud, S) = (0, 0). Equations (3.20) and (3.21) are then rear- ranged such that the nodal variables at each node are grouped together, resulting in the following equations of motion in (finite element) physical coordinates: 92 M)? + KX + G(X) + H(X) = 0 (3.22) where X = lUd,2W292l---|Ud.11+1I'l”.i1+19M+1lT- (323) These equations of motion are next transformed into linear modal coordinates and trun- cated using the following coordinate transformation: X 2 (P77 (3.24) where the 3M x Q matrix <1) is the collection of transverse and axial normal modes of interest of the linearized rotating beam. Here Q = Nt + N0 is the number of kept modes, where N)- is the number of kept transverse (flapping) modes and Na is the number of kept axial modes. This truncation is used to remove unreliable linear modes from the full finite element model. The model is still relatively large and is to be ultimately reduced to a single DOF. The con'iponent form of the equations of motion in these linear modal coordinates is given by Q Q Q 7771+ 4'72. 7In + Z Z anjk 77)- 'le + Z Z Z (371)-17 '77)- Wk 777 = 0 (325) j=l k2] j=l k=l [=1 for 77=1,2,3,...,Q 93 where am are the natural frequencies associated with the 77th mode, anjk are the coefficients of the jk quadratic nonlinear terms associated with the 77th mode, and ,8”ij are the coefficients of the jkl cubic nonlinear terms associated with the nth mode. The formulae for anjk and finjkl can be found in Appendix 3. The determination of these nonlinear coefficients has been computationally automated using the finite element formulation and the coordinate transformations outlined above. 3.3 Modal Reduction 3.3.1 Galerkin-based Invariant Manifold Approach Typically, one is interested in the dynamics of modes that exhibit large amplitudes and/or lie in a certain frequency range, but not in the response of all system modes. To obtain an accurate ROM of a system of coupled nonlinear ODES, one should not simply truncate the linear modes that lie outside the frequency range of interest. The contribution of the truncated modes can have important effects and should be accounted for in the ROM in order to accurately represent the dynamics of the system. The nonlinear model reduction approach based on N N M is a systematic procedure that accounts for the contributions of all linear modes to the NNM of interest, without the direct dynamic simulation of these linear modes. It does so by slaving the linear modes to the dynamics of the NNM of interest in a particular manner. The ROM obtained by this procedure can be made to be very accurate over a large amplitude range if one can accurately construct the invariant manifold corresponding to the NNM of interest. This has been accomplished by Pesheck et al. ([34]), by making use of a Galerkin-based numerical approach to solve the invariant manifold equations. This approach was applied to the rotating beam problem, where a Raleyigh-Ritz approach was used to create a discretized dynamic model of the blade ([34]). Some results from that paper are used for comparisons in the present work. 94 A brief summary of the invariant manifold procedure is reviewed here. It is followed by a description of a solution approach for the invariant manifold equations that is based on collocation methods, which offers some advantages in terms of computational efficiency. From equation (3.25), the mode of interest, i.e., the master mode, is taken to be the mth mode. This linear mode is to be extended in order to generate the corresponding NNM. The generalized master modal position and velocity, 77m and 77m, are expressed in terms of a master mode amplitude, a, and phase, 45, using the coordinate transformation 77m : 0 005(6)) (3.26) Tim : —(I, Wyn Sin(¢). (3.27) Using these dynamic variables, the equation of motion for the master mode can be expressed as two coupled lst order ODEs in a and (b, as follows, — fm Sin (Q) a = ————-——. (3.28) Wm <5 2 Wm _ fm COSfQ), (3.29) 0 Wm Q Q Q Q Where fm Z —( Z: amjk 773' III; + Z Z :78ij 7Ij "It 771 ) (3'30) j:1k=1 j=1k=11=1 corresponds to the nonlinear terms in the 777th equation of motion, which contain coupling effects from all linear modes. The remaining Q - 1 linear modal positions and velocities are slaved to a and at through functions that are to be determined. These are expressed as 95 772' = Pifaa 45) (3.31) 77'.- = QM, <7) (3.32) for 7:1,2,3,...,Q,i# m. Equations (3.31) and (3.32) are substituted back into the components of the equations of motion (equation (3.25)) associated with the slaved coordinates. The functions P,- and Q,- are then obtained using a standard invariant manifold approach. The time derivatives are removed from the equations of motion by the chain rule, resulting in the following nonlinear, time-independent, partial differential equations which govern the geometry of the manifold through the functions F, and Q,: 6P, —- m S. 8P, m 0.- = Eth-fw) + at m — L331?) (333) Jaw.- = aagi(‘f";:“(f)>+gii(wn.—f—”;{f—:@) (3.34) for i=1,2,3,...,Q, 274 m. Equations (3.33) and (3.34) are PDEs in the PS and Q’s and are solved using an expansion of basis functions. Via a Galerkin projection, each PDE, after the expansion is introduced, is projected onto each basis function. These basis functions are products of selected shape functions in the a. and (6 directions. Since 96 has period 277, the shape functions in the (,0 direction are chosen to be a Fourier basis, {1, cos(77¢), sin(776>)}. The shape functions in the a. direction are chosen to be piecewise linear segments. The desired domain in the a direction is divided into K small sections, producing annular subdomains given by (n 96 E [0, 277] and a. E [a], aj + An]. Therefore P,(a, <75) and Q,-(a, qt) can be expressed over the jth interval as: P7(a.¢>) = :Cl’n7lm(a.¢) 1,77. a—aj 7v. ‘ G9 ‘ (l — (1' 23., . , : Z[Ci1'n(—El—J)+ Ci n(1— -—A—a—)] cos((n "1)QD) (3.35) 7721 ' Q7(a.e) = ZDfiP’iUmw) 1,77 Na ' a — a- n a — a- . = Z]D:’"(_£a_3)+Di2’ (1— (13)] Sln(n(,f>) (3.36) 7721 for i=1,2,3,...,Q,-i7ém where the C ’s and D’s are to-be-determined coefficients. (Note that only a particular subset of harmonic functions in (7') are required for this conservative, gyroscopic system, since it possesses synchronous modes [34].) Equations (3.35) and (3.36) are substituted into the manifold-governing equations, equations (3.33) and (3.34), and each of these is projected onto each basis function using a Galerkin projection over the domain of interest. This leads to: 97 , 8T -— , ' ’ 0 = /¢Up.q -aZDi-‘"Uz,n+ZCf’" gift fmasmmh a. lvn 1,71 Wm . 6T. 1, Ln f7, cos(q§) + €01. n a¢ ((1. Wm, "' —”—’u}rm—) d0 dgb (337) -3U7- —f a sin(¢) 0 z: T _ . 2 . CART. __ _ DI,n ,n m -/a.¢> 17.71 617,012”: 2 13,, af2+l§g 2 8a( wm ) 8U , , . + Z DéJl I I." (a w,” _ M) da d¢ (338) In 696 Wm for i=1,2,3,...,Q,77£ m; Equations (3.37) and (3.38) are a set of 2(Q— 1) x 2 X My nonlinear equations in the C ’s and D’s. Note that there are K sets of C’s and D’s, one for each Aa interval. These individual solutions are assembled to construct the invariant manifold over the domain of interest. Once all the expansion coefficients are obtained, the PS and Q’s , which are the slaved linear modal positions and velocities, are known functions in terms of the master variables (a, (1)). These functions dictate how the slaved linear modes must follow the master mode such that the equations of motion are satisfied and the overall motion is invariant. These known functions are used to express fm in equation (3.28) and equation (3.29) in terms of only a and (f7, rendering the ROM single degree of freedom oscillator. The key to the invariant manifold approach is that the solutions of this oscillator represent solutions of the full equations of motion for a particular NNM, and they systematically account for the dynamics of the slaved linear modes. The numerical solution of the invariant manifold equations allows one to obtain NNMS that are accurate over a large amplitude range. 98 3.3.2 Collocation-based Invariant Manifold Approach The computational cost associated with evaluating and solving equations (3.37) and (3.38) can be quite expensive. To reduce these computational costs, the collocation method has been adopted here ([77]). Instead of projecting each manifold-governing equation onto each basis function, it is instead projected onto a set of Dirac delta functions in the master coordinates. This amounts to an approximation of the integrals over the domain, and it leads to: _ . ,. Ln, 21,1187) n( “fm a Sin(¢) 99') [an 8T : + Z 0’; n 0’" (a w,” — N) da da (3.39) Wm 1,71 12:: n X 1,871 Ln. fm a fink?) 1,11 ) +2040 ” —,—’—"( (a cum — f$94—$65 da dd) (3.40) Wm l. n for i=1,2,3,...,Q, 3% m; 19:12, 611—1 ,Na: (12 :11 'aNd)» where (a.p,o,‘)ql) E [(Lj,(lj + An] x[0,27r] are collocation points associated with Qi(a,gb) and ((lpszqoZ) E [an (1]- + An] x [0, 2n] are collocation points associated with 192-((1.46). This projection, while an approximation, greatly simplifies the integrals. 99 3.4 Results The rotating beam parameters used here are the same as in [34]: L = 9 m, m = 10 kg/m, E1 = 3.99 x 105 N -m2, EA 2 2.23 x 108 N, Q = 30 rad/s, and h = 0.5 m, which are feasible for a rotor blade with sub-sonic tip velocity. First, the convergence of the linear natural frequencies is considered as the number of elements is varied, and the results are verified with analytical solutions from [1] and [2]. Then the convergence of the nonlinear amplitude-dependent fundamental flapping natural frequency is studied and verified against the reference model and the collocation-based ROM of the model from [34]. Finally, the time response of the ROM of the fundamental flapping N NM is studied and compared with the reference model and the ROM based on results from [34]. 3.4.1 Linear System Convergence The finite—element model of the rotating beam is obtained by the procedures outlined in section 2. There are 8 cases considered: 5, 10, 15, 20, 45, 90, 136 and 182 elements. The linear natural frequencies of the rotating beam are calculated by solving the eigenvalue problem of the linear part. of equation (3.22). The first four flapping linear natural fre- quencies are shown in Figure 3.3 as a function of the number of elements. These natural frequencies are compared with the natural frequencies from table 7 of [1] and from Table 3 of [2]. In [1], the following dimensionless parameters are defined: ‘1 ‘2 2 mL (2 El ’ z ‘ z __ t 3.41 7’ E1 T mL4 ( ) where a is the dimensionless hub radius, 7} is the dimensionless rotating speed, and T is the dimensionless time. Since the parameters used here do not fall in the parameter range used in [1] and [2], the frequencies must be extrapolated in 7} and interpolated in a. The 100 first four linear natural (flapping) frequencies in the T time scale are presented in Table 3.1. These frequencies are then converted to frequencies in the t time scale and compared with the natural frequencies from the finite-element model in Table 3.2. From Figure 3.3, it is clear that these four linear natural frequencies quickly converge as the number of elements is increased. In addition, from Table 3.2, the maximum difference be- tween the reference natural frequencies and the natural frequencies from the finite—element model is 0.2%. 3.4.2 Nonlinear Model Development The finite-element model of the rotating blade generated by M one-dimensional beam elements has 3M degrees of freedom (with boundary conditions included). Due to the fact that transverse displacement and element rotation degrees of freedom have the same natural frequencies and mode shapes, the natural frequencies that lie in the frequency range of modes 1.5M — 3M may be unreliable ([88]). To be conservative, we consider a frequency range that includes only the first M natural modes (that is, the M with the lowest natural frequencies), and we desire to keep Q g M retained modes for our reference model. For this system it has been found that a reliable approach is to maintain an equal ratio of transverse modes to axial modes when reducing from M to Q. If Nt, M and Na, M represent the number of transverse and axial modes among the first M modes, respectively, then we take Nt = Nt.M x (Q/M) and Na, 2 Na,M x (Q/M) I. These Q modes are used in equation (3.24) to arrive at equation (3.25). The procedure outlined in section 3.3.2 is then applied to the model in linear modal coordinates to obtain the single DOF ROM. In particular, the first linear modal coordinate (first flapping mode) is chosen to be the master mode of interest. The dynamics for this NNM, equations (3.28) and (3.29), are numerically integrated with initial conditions for mode 1 (flapping) amplitude and phase, lSince NM! + NM” 2 M, then N, + N,, = Q, as desired. 101 a(0) = 4.0 and (,b(0) = 0 I. For this amplitude, the NNM frequency is studied as the number of linear modes and the number of elements are varied. For 90, 136, and 182 elements, the number of kept modes (Q) used to generate the frequency data are 15, 21, 27, 33, 39 and 45, respectively §. From the frequency data, it found that as the number of elements is increased, the frequency is nearly converged for 182 elements, over a wide range of values of Q. As the the number of kept modes, Q, increases, the frequency decreases. For 182 elements and 45 modes, the frequency differs from the 182 element, 39 mode case by only 0.0375 rad /s. However, it appears that more than 45 kept modes are needed to demonstrate convergence in a convincing manner. At this point, we employ the 45-mode model generated from 182 elements for use as the reference model. The NNM manifold is obtained from this model over an amplitude range of a E [0,60], which is generated using 120 piecewise linear manifold segments in a with width Aa = 0.05, and N¢ = 16 harmonics in (p'. This amplitude range, in physical terms, corresponds to the beam tip moving with a peak—to—peak amplitude of about 2.4 meters (recall that the beam is 9 meters in length). The single DOF ROM is obtained by restricting the dynamics of the 45 DOF reference model to this invariant manifold. This amounts to restricting the dynamics, which occur in a 90-dimensional state space, to a two-dimensional invariant manifold that is described by the master-slave relations given by the Pis and Qi’s. Note that there are 88 such relations in this case. Figure 3.4 depicts the natural frequency of the first nonlinear flapping mode as a function of the initial mode amplitude, a(0), with ¢(0) = 0. The + line is a result of the ROM (1 DOF) which is generated from the finite element reference model. The 0 line is a result of the ROM (1 DOF) which is generated from the Rayleigh-Ritz model from [34], using the same N3 and Na as used in the reference model. The solid line is a result of the full reference model (45 DOF), which is obtained by searching for periodic solutions using a shooting algorithm 1] The results from the three models are virtually identical over this 1T his amplitude corresponds to the beam having a peak-to—peak tip vibration amplitude of 1.59 meters. §Due to computer memory limitations, the maximum number of kept modes we can achieve is 45. 1Note that all modal solutions based on the reference model require use of this shooting algorithm. 102 amplitude range. 3.4.3 System Response We now have a single DOF ROM (equations (3.28) and (3.29)), which is generated from 182 elements and 45 kept modes. This ROM is simulated with initial modal amplitude (1(0) 2 5.9 and phase 95(0) = 0, and the result is shown as the dashed line in Figure 3.5. This Figure can be geometrically interpreted as the projection of the motion that occurs on the NNM manifold onto the linear eigenspace of the first linear flapping mode. This numerical solution is verified against the numerical solution of the ROM from [34] with the same initial conditions, which is shown as a dotted line. These two solutions are nearly identical. Figure 3.6 depicts the corresponding time responses of the N NM displacement 771 obtained from three different models. The dashed line is the time response of the ROM obtained by applying equation (3.26) to the numerical solution shown in Figure 3.5. Obtained by applying the same procedures, the dotted line is the time response of the ROM from [34], while the solid line is the time response of the reference model (45 DOF) initiated on the NNM manifold. These three responses are nearly identical. Figure 3.7 depicts one of the attendant master-s]ave-constraint relations, specifically the second linear flapping mode displacement as a function of the NNM amplitude (a) and phase ((1)). Note that solutions of the ROM are used to obtain (curb), after which the second linear flapping mode displacement is obtained by the corresponding trace on this surface. Of course, there are 87 other such surfaces, 43 more slaved displacements and 44 slaved velocities, for the linear modes. Note also that the contribution of this (and all) linear mode(s) start at zero at a = 0 and increase in amplitude as 0. increases. Figure 3.8 depicts the time response for the second linear flapping mode displacement 103 obtained from the ROM, the ROM from [34], and the reference model. The dashed line is obtained by applying the master-slave constraint (depicted in Figure 3.7) to the numerical solution of the ROM (the dashed line in Figure 3.5). The dotted line is obtained by applying the master-slave constraint from [34] (which is similar to Figure 3.7 but is not shown here) to the numerical solution of the ROM from [34] (the dotted line in Figure 3.5). The solid line is obtained from the 45 DOF reference model. Again, the results are virtually identical, demonstrating that the invariant manifold accurately captures the effects of higher linear modes. Further sample results are shown in Figure 3.9, Figure 3.11, and Figure 3.13, which depict the master-slave constraint relations for the tenth flapping mode, the first axial mode, and the twelfth axial mode, respectively. Associated with these are Figure 3.10, Figure 3.12, and Figure 3.14, which are similar in nature to Figure 3.8, and depict the slaved time responses for the tenth flapping mode, the first axial mode, and the twelfth axial mode, re- spectively. In all cases the N NM manifold accurately captures the contributions from these linear modes. It is interesting to note that the first linear axial mode makes a significant contribution to the fundamental flapping N N M, thus demonstrating the importance of the nonlinear axial / transverse coupling. The 88 master—slave—constraint relations will generate the NN M manifold for the first non- linear flapping mode in the 90-dimensiona1 state space. From Figure 3.7 through Figure 3.14, we see that the responses of the slaved linear flapping modes have the same fundamen- tal frequency, as that of the N N M. However, the axial modes move at twice the frequency of the NNM. This arises from the physics of the beam motion, and is discussed below. Figure 3.15 depicts time responses of the beam-tip transverse deflection obtained from the present ROM, the ROM from [34], and the reference model. The dashed line is obtained by applying equation (3.24) to the time responses of the linear flapping and axial modes of the ROM, and plotting the result for the transverse deflection of the end node of the beam. Obtained by applying the same procedures, the dotted line is the result from the ROM from 104 [34], and the solid line is the result from the reference model. Note that three responses are very close and correspond to a tip-to-tip deflection of 2.34 meters. The discrepancy in amplitude is, we believe, related to convergence difficulties with the Rayleigh-Ritz approach near the beam tip [3]. Figure 3.16, Figure 3.17, and Figure 3.18 depict the spatial nature of the transverse, angle, and axial motions of the beam in the NNM response for the same initial conditions as the previous Figures. They represent the motion starting from t = 0 up to a quarter of the period of the NNM. The time intervals between each deflection curve are equal. The bottom curves of Figure 3.16 and Figure 3.18 and the top curve of Figure 3.17 are at the initial time, corresponding to the peak NNM amplitude. The top curves of Figure 3.16 and Figure 3.18 and the bottom curve of Figure 3.17 are one quarter-period later. The dashed line in Figure 3.18 (close to the top curve in the Figure) is the static axial extension us due to rotation of the beam. At the initial time, the beam experiences axial foreshortening due to nonlinear effects through the transverse deflection, as can be seen by the negative axial displacements in Figure 3.18. At the quarter period, the transverse deflection is at zero, and the beam does not experience axial foreshortening, but deforms axially due to centrifugal loads. Note that during the first quarter-period, the transverse and angle responses go through a quarter-period of the periodic oscillation, but the axial response goes through one half of its full oscillation. This is due to the symmetry of the axial response about the zero transverse and angle deflections condition, that is, the axial foreshortening is the same for —w and —6 as for +11) and +6. This is why the axial motion has twice the frequency of the transverse and angle motions. By combining the transverse deflection from Figure 3.16 and the axial deflection from Figure 3.18, the deformed shape of the rotating beam on the first nonlinear mode manifold can be generated. Due to the small scale of the axial motions, the beam shape will be very close to the shapes shown for the transverse deflection (this is not true at all amplitudes, but is a very good approximation here). It is important to note that the shape of the beam changes in a periodic manner through the course of the NNM motion. This is due to the 105 contributions from the slaved linear modes, which are not fixed, but vary in time according to the NNM invariant manifold. This is in distinct contrast to more typical reduction methods in which the equations of motion are projected onto a single mode shape, resulting in a response during which the shape of the beam remains the same throughout the course of the motion. 3.4.4 Computational Considerations An important issue related to this procedure is that of computational time. As a compar- ison, we consider here the relative times used to generate the invariant manifold solution and run simulations on it with those of generating the manifold using a shooting technique. It takes about 1,600 seconds to generate one strip of the manifold, yielding a total of about 192,000 seconds to generate the manifold. This is not trivial, but with the manifold in hand, a simulation for a set of initial conditions takes under 500 seconds for about two periods of response. In contrast, if one generates the manifold directly by using a shoot- ing technique to find the correct initial conditions, each strip takes nearly 8,400 seconds to generate, a factor of five higher. And, note that this approach uses the master-slave relations as starting points for the shooting technique, which are not generally available, making this result quite conservative. When parametric studies are required, it is possible that the ROM using invariant manifold techniques can be extended to handle the situation. A proposed approach is to augment the system by considering the parameters of interest as elements of the master mode co- ordinates. This is the so—called “suspension trick” ([89]). By this approach, the single nonlinear mode manifold of interest will also be a function of the parameters of interest. Hence, the computational cost of solving the invariant manifold does not become a recur- ring cost when the parameters of interest are changed, however the initial cost of generating the “suspended” manifold will be higher. This should be feasible and efficient for a single parameter. 106 3.5 Conclusions In this chapter, the problem of obtaining nonlinear single—mode ROMS from a finite element formulation of a rotating beam was addressed using N NMs. Accurate simulations of this system typically require many DOF, due to the nonlinear coupling effects between axial and transverse deflections. The problem of obtaining accurate ROM was solved using invariant manifold-based NN Ms. The ROMS obtained govern the dynamics of the NNM of interest (the master mode), and the approach allows one to accurately capture the dynamics of the slaved linear modes. Hence, only a single nonlinear oscillator is required for accurate simulation of the NNM out to large amplitudes. The numerical solution of the NNM invariant manifold equations was facilitated by employing a collocation method. Excellent agreement was found for the dynamics of the fundamental flapping mode of the beam, by comparing results from the single DOF ROM obtained here with the single DOF ROM based on a Rayleigh-Ritz discretization ([34]) and a full 45 DOF reference model. The method can be similarly applied to any mode of interest, and can be generalized to include dissipation effects, both linear and nonlinear. The results presented demonstrate significant promise, since they combine procedures from the versatile finite element method with the accurate reduction procedures via numerically generated NN Ms. This will allow the NN M method to be applied to the large amplitude vibrations of beams with more complex geometry and more nodal DOF, as required for more realistic models of rotor blades and other structural members. Also, current work on multi-mode manifolds shows promise for developing nonlinear ROM with two DOF in a similar manner. 107 3.6 Tables Table 3.1. Dimensionless linear natural frequencies. wl and tag were computed using results a = 0.06 77 w1 w2 w3 w4 12.17 13.8248 38.9163 81.37 142.606 from [1]. (.123 and M4 were computed using results from [2]. h :2 0.5 in Q w1 w2 w3 w4 (rad / s) Analy. FEM Analy. FEM Analy. FEM Analy. FEM 30 34.09 34.03 95.97 95.84 200.67 200.49 351.68 351.49 Table 3.2. Comparison of the linear natural frequencies obtained by power series and finite element methods. The analytical results for an and UJ2 were computed using results from [1]. The analytical results for w3 and (124 were computed using results from [2]. 108 3.7 Figures Q We(xe:t) D I fees.» Vii/Ci/r/ Ungt) 8:1, 2, 3,...,M V///, J l l lLlell I I ll I] 1 U _ x. x...” \E,A,I,L,m h 3: Figure 3.1. Finite element representation of a rotating beam with Q = constant. 109 t I j j T I —l— AnaIytical Sfilution q 0012 + 182 Elements 0.01 r . ' . 0008* 4 :5. 3 0.003 - U) 3 0.004- a 0.002[ , G L l i 1 1 l l 1 0 1 2 3 4 5 6 7 8 9 X (m) Figure 3.2. Comparison of the static axial blade deflection, 'Us, obtained from an analytical solution [3] and by the finite element method. 110 53 1st Linear (Flap) Natural Freq (rad/s) 8 . 2nd Linear (Flap) Natural Freq (rad/s) co .‘1 or 1 $3 8 96.5 ’ K 2 i 1 1 I r r 2% 50 100 150 200 0 50 100 150 200 Number of Elements Number of Elements 23 o m co 0: '01 O 203 (rad/s) c.) or x: q - 202.5 1 1 00 01 O) 202 f 3 201.5 1 201 [ + 0) 01 Q) 200.5 ' v if + (a) (II N fin- l l l I F I I #A 4A 50 100 150 200 50 100 150 200 Number of Elements Number of Elements § 4th Linear (Flap) Natural Freq. (rad/s) 8 a 3'“ Linear (Flap) Natural Freq 0 Figure 3.3. Convergence of the first four linear natural frequencies (flapping modes) of the finite element model. 111 34.4 . , , . . ’0? + ROM from FE ‘/4>—4> g 34.35» 0 ROM from Rayleigh-Ritz ,-a/ - b —- Reference Model / >‘ 34.3 » - 8 fl 8 ,/ ,ar g, 34.25 » // _ u. [,1 E 34.2 » ,/ - 2 ,ar :5 / E 34.1 5 l' ,1"! _, 8 ,/ g 34.1 r ,/ - C ,/ 0 /gr .2 34.05% ,/ q 0’ _,_ -49 34 g ‘ . 1 . O 1 2 3 4 5 6 Mode (1St Flapping Mode) Amplitude Figure 3.4. The first nonlinear flapping natural frequency as a function of initial mode am- plitude a(0) with ¢ = 0 obtaining using three methods: the FE—based ROM, the Rayleigh- Ritz based ROM, and the reference model (182-element and 45-kept-modes). — — ROM from FE 5.945 F ROM from Rayleigh—Ritz ‘ 5.94» ,1 ‘\ / “\ ~ A / \\ / \ g 5.935- x \\ ,/ \\ _ Q) / ‘x / \ / I :3 5.93» ,’ ‘, ,’ \. 3 z: \ ‘ Q. E 5.925r ,’ ‘, ,-’ ‘, - < I \ I \ - . 2» ' ‘ ’ ‘ 4 a 5 9 1’ \ I’ \ g 5.915» ,’ l, I \\ ~ ’I ‘ I, l 5.91 r- r ‘\ I ‘\ " [I \ I' \ 5.905— r \ I \ ~ I '\ / \ 5 9 ' 1 1 “ 4' ‘ Q Mode 1 Phase Angle ((1) Figure 3.5. The projection of the motion occurring on the first nonlinear flapping mode manifold onto the linear eigenspace of the first linear flapping mode, in amplitude and phase coordinates. 112 - - ROM from FE E .g 8 g) 8— ROM from Rayleigh— —thz '07 6‘» _ E I- 4" ' 2'3 3 2— - 2 05 or 8 2 _21 0" .E ca. _4r 3 LL -6‘ 257- : —8— ~ <1: . . . . n A E O 0.05 0.1 0.15 0.2 0.25 0.3 Time Figure 3.6. The time responses of the first nonlinear flapping mode displacement m for the FE-based ROM, the Rayleigh-Ritz based ROM, and the reference model. P2(a,¢) (2nd Flapping Mode) Figure 3.7. The contribution of the second linear flapping mode 772 = P2(a,¢) to the first nonlinear mode manifold, solved by the collocation method, as a function of the first nonlinear mode amplitude and phase. 113 - - ROM from FE ROM from Raylelgh-thz 4 ——~ ~ Reference Model .0 o: .0 a .0 m I l .o .o A M 2"d Flapping Mode Deflection O l .0 a: I J 0.05 0.1 0.15 0.2 0.25 0.3 Time 0 Figure 3.8. The time response of the second linear flapping mode deflection 172 = P2(a(t), ¢(t)) on the first nonlinear mode manifold for the FE—based ROM, the Rayleigh- Ritz based ROM, and the reference model. P10(a,¢) (1 0th Flapping Mode) Figure 3.9. The contribution of the tenth linear flapping mode 7710 = P10(a,¢) to the first nonlinear mode manifold, solved by the collocation method, as a function of the first nonlinear mode amplitude and phase. 114 — — ROM from FE ' ROM from Rayleigh—thz ‘ — Reference Model ,5 1.5 » - 13 - a) _ 8 2 ~~ _ U) C — ‘a. E - LI. 5 o -l .l 0.1 0.1 5 0.2 0.25 0.3 Time Figure 3.10. The time response of the tenth linear flapping mode deflection rho = P10(a(t), ¢(t)) on the first nonlinear mode manifold for the FE—based ROM, the Rayleigh- Ritz based ROM, and the reference model. P22(a,¢) (1“ Axial Mode) Figure 3.11. The contribution of the first linear axial mode 7122 = P22(a,¢) to the first nonlinear mode manifold, solved by the collocation method, as a function of the first nonlinear mode amplitude and phase. 115 - - ROM from FE - ROM from Rayleigh-Ritz -— Reference Model 0.05 -0.05 —0.1 » / —O.1 5 ~ —O.2 * -O.25r / \ -O.3 r / l \ -335.) . (Z —o.4 “ o I 3‘ Axial Mode Deflection 1 0.15 0.2 0.25 0.3 Figure 3.12. The time response of the first linear axial mode deflection 7722 = P22(a(t), ¢(t)) on the first nonlinear mode manifold for the F E—based ROM, the Rayleigh-Ritz based ROM and the reference model. 1 x10 6\ P33(a,¢) (12th Axial Mode) Figure 3.13. The contribution of the twelfth linear axial mode n33 = P33(a,¢) to the first nonlinear mode manifold, solved by the collocation method, as a function of the first nonlinear mode amplitude and phase. 116 I [ - - ROM from FE -l ROM from Rayleigh-Ritz ~—— Reference Model 8 7 5 z..- 6* .33 \ E‘S’ 5* 5 4* 2 l .5 3“ , E 2i- \ E N 'r- 1,_ \ o» \v \«' -1 l 1 J A A I O 0.05 0.1 0.15 0.2 0.25 0.3 Figure 3.14. The time response of the twelfth linear axial mode deflection 7733 = P33(a(t), ¢(t)) on the first nonlinear mode manifold for the FE—based ROM, the Rayleigh- Ritz based ROM, and the reference model. 1.52 — - ROM from FE _ ROM from Rayleigh-Ritz “\ —— Reference Model , [ \ 1 \ a \ E, 0.5[ . c _O E o. .9- -o.sl . l— / -1 >— / / -l —1 .5 r 3 0 0.05 0.1 0.1 5 0.2 0.25 0.3 Time Figure 3.15. The time dependence of the beam-tip transverse deflection, w(L, t), for the FE-based ROM, the Rayleigh-Ritz based ROM, and the reference model. 117 l- .— 1 I 1 0123456789 X(m) r— p— —1.4 m Figure 3.16. Transverse deflection ’lU(LlI, t) for a quarter-period of motion on the first non- linear mode manifold. The motion starts at the maximum deflection (the bottom curve) and moves as shown to the zero deflection at a quarter-period, after which it moves upward in a symmetric manner about 11) = 0 and then repeats. 118 0.18 0.16“ 0.14~ 0.12’ 0.11 9(X) (m) 0.08 - 0.06 - 0.041 0.02[ ' x (m) Figure 3.17. Angle deflection 6(23, t) for a quarter-period of motion on the first nonlinear mode manifold. The motion starts at the maximum deflection (the top curve) and moves as shown to the zero deflection at a quarter-period, after which it moves downward in a symmetric manner about 0 = 0 and then repeats. 119 0.02 l T I I I I I I 0.01 —0.01 —0.02 -0.03 UlX) (m) -0.04 -0.05 -0.06 —0.07 l l l l n— l— 41080123456789 x (m) Figure 3.18. Axial deflection u(:1:, t) for a quarter-period of motion on the first nonlinear mode manifold. The dashed line denotes the static deflection, 113(23). The bottom curve corresponds to the initial condition. The top curve from Figure 3.16 and the bottom curve from Figure 3.17 correspond to the top u(:1:) curve here, all at one quarter-period. Note that the axial motion occurs at twice the frequency of the transverse and angle motions, due to symmetry. As time progresses beyond the quarter—period, the axial motion moves downward again, and oscillates between the top and bottom curves shown. 120 3 .8 Appendix 3 The coefficients of the Jk quadratic nonlinear terms associated with the nth mode are as follows: Ill-F1 . LI , ‘1? , e 0an = Z (0(3e—2ln G'ejk + 99(3e—1)n Oejk + @315)” 063-1.) (342) 821 {n,J', It 2 1,2,3, Q} (1W eJk’ aejk and \vhere 90,-]- is the Ij component of the modal matrix in equation (24), and 0U ,9 , . 061k are. 4 4 U U aejk : :1 "2 Gel Injk (3°43) 1141112 W’ W aeJk 21 "1:10 aelrnjk (3'44) (41 m2 8 - “BI/r :1 mzlae Qelrnjk (3.40) N m1,2,3, ...,M +1}; {J',k =1,2,3,...,Q} PM rs ll ,9 . . . elnjk’ 06:1ka and aelmjk (“9‘ 121 ,U aelnjk a 11' eI rnJ'k (aglnIe ¢(38+I—5)j ¢(3e+n—5)k + (allInIE ¢(3€+I—2)j ¢(38+n—2)Ic {1:12}; {n :12}; {e = 1,2,3,...,!v1 +1}; {13k =1»2131~-1Q} u e—l . u e (0211.) ¢(3e+l—5)j ¢(3e+n—4)l~ + (Gun) ¢(3e+l—2)j ¢(3e+n—1ll.~ {[21’2}; {71: 3’4}1 {8 :112131""A/I +1}; {jak : 1,2131°"1Q} (a’é‘inle‘l ¢(3e+I—4)J' ¢(36+n—5)k + (aflrrle ¢(Be+I—1)J' ¢(3e+n—2)k (3-45) {1: 3,4}; {71:12}; {6 =1,2,3,...,M +1}; {J',k =1,2,3,...,Q} )e—l e . U . ' ' (“2m C29(3e+l—4)j ¢(36+n—-4)k + (“(11111) ¢(3e+l—1)j ¢(3e+n—1)k {1: 3,4}, {71: 3,4}, {8 : 1,2,3,-~1A1 +1}; {jk =112131'°'1Q} (”firmly—1 ¢(3e+l—5)j ¢(3e—5)k + (aII~nz)e ¢(3e+l—2)j ¢(3e—2)k {1:192}? {771:1}? {6 : 19213IH'1AI +1}’ {jfk :1’2’3’°"’Q} . 3—1 ' ~ (”filmy ¢(:3e—+l—5)j @(3e—2lk + (ail-"116 ¢(3e+l—2)j ¢(3e+1)k {121,2}; {m = 2}; {e =1,2,3,...,M +1}; {j,k = 1,2,3,...,Q} (”153111118—1 ¢(3e+l—4)j ¢(3r—5)l~ + (aIImIe ¢(3e+l—1)j 9(3e—2)k {I : 3,4}; {77121}; {e = 1,2,3,...“71/1 +1}; {j,k =1,2,3,...,Q} I s e—l s e (”31:71) @(3e+I—4)J CTb(3(3—2)Ic + (allrn) ¢(3€+I‘~I)j Q75(3e+1)k {I = 3,4}; {m = 2}; {e =1,2.3,...,.M +1}; {Ik :112I3I”'1Q} 122 Q ._ 3 0-1 3 e ' 0(1ka — ((141171) ¢(3e+1—5)j ¢(3e—5)k+(“21m) ¢(3e+l—2)j ¢(36—2)k {121,2}; {77121}; {e 21,2,3,...,]\II +1}; {j,k =1,2,3,...,Q} ._ —1 4 . . , = (“3sz $(3e+z—5)j ¢(3e—2)k + (aglmf ¢(3e+1—2)j ¢(3e+1)k {1:1’2}; {771: 2}; {e =1,2,3,...,AI +1}; {],k =1,2,3,...,Q} v' c —1 . = (031,”)6 ¢(3e+1—4)j @(3e—5)k+(031m)e ¢(3e+1—1)j¢(3e—2)k (3-48) {I = 3,4}; {m =1}; {e =1,2,3,...,A'I +1}; {J'Jf- =1,2,3,---,Q} ‘3 ‘_1 . , . = (“-317er ¢(3e+1—4)j ¢(36—2)k + (03177;)6 ¢(3e+z—1)j ¢(3e+1)k {1: 3,4}; {m = 2}; {e =1,2,3,...,M+1}; {j,k =1,2,3,..-,Q} vvhere (-)e—l denotes quantities associated with element 6 — 1, (~)" denotes quantities asso- ciated with element 6?, and of“ ”n and (121m are quantities defined by equations (9) and (13), 1‘ es pect ively. The coefficients of the jkl cubic nonlinear terms associated with the nth mode are as fEDIIOWS Al+l , , ,/ w .~ ,1 e flnjkz = Z (0(3e—l)n flejkz + ¢’(3e)n fjejkl) (3-49) 621 {71,j,k,l = 1,2,3» --:Q} n a o / ’1, / w}1ere (bl-j is the ij component of the modal matrix In equation (24), and c.1523“ and .52“ flare : 123 4 4 33kt = Z 22:33,:li vvhere [3le and [[33 e nrrjkl 3."er are: 124 (3.50) (3.51) ,/ W "361*nrjkl __1 , [( 31*721')e ¢(3e+l*—5)j ¢(3€+n-—5)k ¢(36+r—5)l +(b‘fp.,,,.)e ¢(3e+z*—2)j ¢(3e+n—2)k¢(3e+r—2)z} {1* 21,2}; {n =1,2}; {r =1,2}; {e =1,2,3,. .—1 [( gm")? ¢(3e+z*—5)j ¢(3e+n—5)k ¢(3e+r—4)1 +( ‘iz*n1~)e ¢(3e+1*—2)j ¢(3e+n—2)k¢(3e+r—1)z} {1* 21,2};{71 = 1,2}; {1" = 3,4}; {6 =1,2,3,. s .—1 - }( 3m”)? ¢(3e+1*—5)j ¢(3e+n—4)k ¢(3e+r—5)1 +(bf1*n,.)e ¢(3e+z*—2)j ¢(3e+n—1)k¢(3e+r-2)t} {1* 21,2}; {n = 3,4}; {7‘=1,2}; {e =1,2,3,. [( 31“an ¢(3€+l*—5)j ¢(38+n—4)k ¢(3e+r—4)l +(’)f[*,,r)e ¢(3e+l*—2)j ¢(3e+n—1)k¢(3e+r—1)1} {1* 21,2}; {71 = 3,4}; {7‘ = 3,4}; {6 =1,2,3,. , ,—l , [(bgfinr)? ¢(36?+l*-4)j ¢(3€+Tl*5)k gD(38-+-1"-—5)l +(b}[*m‘)e ¢(3e+l*—1)j ¢(3e+n—2)k¢(3e+r—2)l} {1* = 3,4}; {n = 1,2}; {7" 21,2}; {6 =1,2,3,. . —1 1 . [( §l*-nr)€ ¢(3e+1*—4)j $(3e+n—5)k ¢(3e+r—4)1 +1 'iz*m~)e $(3e+1*—1)j ¢’(3e+n—2)k¢(3e+r—1)l} {1* = 3,4}; {7121,2}; {1‘ = 3,4}; {(3 =1,2,3,. .' —1 . 1 [( Sunf 0(3e+1*—4)j l +( 21‘2”)? ¢(3()+l*—l)j ¢(3e+n—1)k¢(36+7‘—1)1} {1* = 3,4}; {n = 3,4}; {7~ = 3,4}; {6 = 1,2,3, 126 M +1}; {j,k,l=1,2,3,... A] +1}; {j,k,l=1,2,3,... M+1}; {3,3,1 = 1,2,3,... M+1}; {j,k,l=1,2,3,... M+ 1}; {3,132 = 1,2,3,... M+1}; {j,k,l=1,2,3,... ,1” +1}; {j,k,l=1,2,3,... ,3! +1}; {Jk7,l=1,2,3,... (3.53) ,Q} ,Q} ,Q} ,Q} ,Q} ,Q} ,Q} ,Q} w11ere (-)€-1 denotes for quantities associated with element 6 — 1, (-)e denotes quantities zlssociated with element 6, and bZl*nr are quantities defined by equation (15). 127 CHAPTER 4 Component Mode Synthesis Using Nonlinear Normal Modes 4 - 1 Introduction 1\-I any complex structures are composed of several relatively simple substructures that are assembled together. This occurs in trusses, bladed-disk assemblies in turbine rotors, aerospace and ground vehicles, and many other applications. In such cases it is convenient to develop a dynamic model for the overall structure by taking advantage of the dynamic properties of the substructures. Component Mode Synthesis (CMS) was developed using these ideas, in order to synthesize models that are described in terms of the component Structures and to take advantage of model size reduction carried out at the substructure leVel. There are two general types of CMS methods; they are known as the fixed-interface and the free-interface approaches. The fixed-interface CMS technique, developed by Hurty ( [60] ,[61]) and simplified by Craig and Bampton ([62]), is widely used, since the reduction procedure is straightforward and typically produces highly accurate models with relatively few component modes ([64]). Free-interface CMS methods are more attractive than fixed- interface CMS methods when the component modes are obtained from modal testing or ‘Vhen an experimental verification of the component modes is required ([65]). The free- interface CMS technique developed by Craig and Chang ([66], [67] and [63]) is the most accurate among the free-interface CMS techniques. It is a modified version of Rubin’s 128 Inethod [68] and MacNeal’s method [69]. It is superior to the CMS of Craig-Bampton in terms of accuracy, but is more difficult to implement ([64]). An extensive review on CMS can be found in [59]. CMS-type methods are well developed for linear structural models alld they have been used extensively, especially in the aerospace industry, see ([55], [56], [57] , and [58], for example). T11e finite element (FE) method is used in CMS to model and characterize the behavior of the individual substructures. Often, accurate models of the component structures require a large number of FE degrees of freedom (DOF). For linear systems, model reduction can be achieved using linear modal analysis, wherein one simply truncates the higher modes of vi- bration. In fact, this is how one carries out model size reduction at the substructure level in linear CMS. However, the dynamic analysis (e.g., the determination of natural frequencies, rnode shapes, time responses) of such structures can require considerable computational effort, and this is especially true for nonlinear structures ([5]). The difficult issue of model size reduction for nonlinear structures continues to be a ma- jor Challenge for computational vibration analysis. For relatively simple systems one can use nonlinear normal modes (NNMS) for this task. For example, Shaw and Pierre have developed systematic procedures to obtain reduced-order models (ROM) via NNMS using invariant manifolds ([20], [21], and [74]). Asymptotic series were initially used to approxi- mate the geometry of the invariant manifold, and this approach was applied to the study Of a variety of systems, such as a nonlinear rotating Euler-Bernoulli Beam ([3]). Pesheck e t al. developed a numerically-based Galerkin approach to calculate the geometry of the NNlVI invariant manifolds out to large amplitudes of vibration ([33]). These procedures Can be used for quite general nonlinearities over a wide range of amplitudes, and they have been applied to many systems, including the rotating beam ([34]). Recently, these methods h ave been shown to be applicable in conjunction with nonlinear FE models ([35]), which Opens a new frontier for their application to more complex nonlinear structural systems. Jiang et al. have recently extended the work of Pesheck et al. to calculate the geometry 3f Inulti-Inode invariant manifolds that are able to capture the dynamics for Situations 129 where more than one mode is active, for example, in the case of internal resonance ([42]). However the computational cost involved in generating the multi-mode invariant manifold is still quite high. The present chapter describes recent research that is aimed at extending the fixed-interface linear CMS method to nonlinear structures by making use of the fixed-interface NNMS of the component structures. By synthesizing the reduced-order representations of the substructures, one can obtain accurate low-order models of structures that are composed of assemblies of nonlinear substructures, with significantly less computational cost than by computing the ROM directly from the fully coupled system. It is found that such an approach is valid, so long as the coupling between substructures is relatively weak. It is also shown how this method relates to the usual linear CMS, and reduces to it under linearization of the model. The chapter is outlined as follows. We first review the development of N NMs, as needed for the individual substructures. The associated invariant manifold equations, parameterized by modal position and velocity, are formulated, and a numerical collocation method is presented that allows one to obtain the solution of the NNM invariant manifold out to moderate amplitudes. Then, the procedures for fixed-interface nonlinear component mode synthesis are described. A five-DOF system and a forty-one—DOF system are used as examples to demonstrate the method, and some conclusions are drawn to close the chapter. 4.2 Formulation of the N NM Invariant Manifolds 4.2.1 Invariant-Manifold Governing Equations We begin with a general discrete representation of the vibrations of a nonlinear structural system (in our case, a subsystem component), obtained either by a. finite element model 130 followed by linear modal expansion, or by a Rayleigh-Ritz approach using the linear mode shapes. We assume that the system at linear order is undamped. (This assumption greatly simplifies the problem, but can be relaxed, in principle.) In this case, the equations of motion for a Q-DOF system are uncoupled at linear order and can be expressed in the form: Iii+An=f(n,7i) (H) where I is the identity matrix, A is the diagonal matrix of squared linear natural frequen- cies, f (17, 1'7) a vector of nonlinear forces, 7) the modal position vector, and 1'7 the modal velocity vector. The component form of equation (4.1) is given by 772+w12 m = fina’lj) (4-2) for i,j=1,2,3,...,Q where w,- is the linear free vibration natural frequency of mode 2' and Q is the number of retained linear modes. In order to search for a particular individual NN M, it is assumed that the N NM manifold is parameterized by a single modal position-velocity pair corresponding to the mode of interest, referred to as the master mode. This is accomplished by using the fact that for a N NM response all of the remaining modal positions and velocities are slaved (constrained) to this master mode. For the kth nonlinear mode, we take uk 2 77k» and Uk 2 77)C as the master states. The remaining slave states are expressed as 131 7h“ = Xifuka‘l’k) = Xi(77k)71k) (4.3) '77: = W‘llka’l-‘H = Wnkflik) (4-4) for 1: = 1,2,3,...,Q; 2' 7f k. Equations (4.3) and (4.4) constitute a set of constraint equations that are to be determined. The constraint functions in equations (4.3) and (4.4) are obtained by an invariant manifold procedure that generates equations that can be solved for the unknown constraint relations. The process begins by taking a time derivative of the constraint equations, yielding ax,- ax,- : ,. 4. 77! (91% “It avk I). ( 5) .. BY, . 81’,- , - = —' . —— .’. 4. 7]? 011}: ”k + 8L1.“ ( 6) for i=1,2,3,...,Q; 3% k. The time dependence in these equations is eliminated by substitutions that make use of the the equations of motion, as follows: 1'1). 2 Uk, bk 2 ijk 2 flag m. + fk(nj,1'7j), 7),- = —w,2 7},- + f,(-r)j,'r'}J-). Then, the constraints (equations (4.3) and (4.4)) are substituted in the resulting expression everywhere in place of the slave state variables, resulting in a set of partial differential equations for the functions (X z(21k, ek), Yi(uk, 'vk)). This set of 2Q - 2, time-independent, partial differential equations govern the geometry of the kth manifold, and are given by 132 2 o —wz X2+f2(Xja)/Jiuk9vk) : Uk+—' for i,j=1,2,3,...,Q; 12¢ k. These equations are not solvable in closed form (except in very special cases). Methods for obtaining approximate solutions for X ,j and Y,- are described in the next section. 4.2.2 Invariant Manifold Solution Asymptotic expansions can be used to obtain approximate solutions of the invariant man- ifold equations, equations (4.7) and (4.8), for smooth nonlinearities, but such solutions are only locally valid ([20], [21]). Numerical Galerkin-based solutions have also been developed, wherein one can utilize either local or global basis functions to describe the invariant man- ifold ([33]). In this work we employ a Galerkin-based procedure that relies on a patchwork of local solutions to obtain a solution that is valid over a given domain: uk 6 [—Ub, Ub] and “k E l-le Vbl- We describe the development of the method for a local domain described by ui E [—ub, ub] and v}: E [—v),,vb]. Once this procedure has been developed, it can be applied to each patch, and the final result is obtained by collecting the results for all patches such that the entire domain is covered. The relations between the coordinates in the global domain , , 4 ° , ' . f? . (.9 ' ~ (uk, wk.) and those in the local domain (uh, Uk) are given by 133 u k 'U k where du and dU du + “Z, (4.9) (iv + 22):, (4.10) 2U ‘Ub + (TV—5X6“ — l) + ab, (4.11) u 2V —V:. + (7,?er — 1) + vb, (4.12) where eu and en are the patch indices in the uk and U), directions, respectively, and N5 and N5 are the number of patches used in the uk and vk directions, and du and dv represent the shift from the ori 'in of u ., v. to u‘five, , res ectivel I. Here 6,, runs from 1 to N e 8 k A k 1, P l u and ev runs from 1 to N5. Substituting equations (4.9) and (4.10) into equations (4.7) and (4.8), the partial differential equations governing the geometry of the kth manifold in the local coordinates (uz, vi.) are given by: BX' . (9X; Y2. : 875% + vi) + 8125 (flag ((1,, + ui.) +fk(vaYjad'u+uiadv+v2)) (4.13) 2 X' , _ , e e _ BY; 6 Yi . 2 , e “’00)“ z+fz(X]aYJadu+Ukadv+vk) — 51"?de +vk)+51_'5(—wk (d‘u‘l‘uk) k ’k + fk(X)-, Yjadu + ui..dv + 222)) (4-14) for i,j=1,2,3,...,Q; 2',j 34 k. The solution of equations (4.13) and (4.14) on the local domain is obtained by expanding (X1, Y2) in terms of basis functions as: 134 Npm NW! . l, X,(u‘,:.,vf;) = Z: Z CimLMl,m(uZ,v}:) (4.15) [:1 "1:1 Npm NW 1 WM) = Z Z Di‘mLA/Iz,m(uti,vi), (4.16) 1:1 m=1 where LAI),.,n(u‘Z,v;) = T1_1(uZ/ub)XTm_1(vf./vb), (4.17) where T)_1(:L‘) and Tm_1 (1') are standard Chebyshev polynomials defined over :1: E [- 1, +1], and the C ”s and D’s are the to-be-determined expansion coefficients. Equations (4.15) and (4.16) are substituted into the local manifold-governing equations, equations (4.13) and (4.14). Normally, each of the resulting equations is projected onto the basis functions, but here we employ a collocation method, which is computationally more efficient, yet retains very good accuracy. This is carried out by projection of the equations onto Dirac delta functions in the local master coordinates over the local domain, as follows: 1, 1,, BLM), 0 : ff’ 6 6(212—ui‘pwi—vgq) — 2 DZ. mLIWLm + Z Cz- m(dv + ”fl—Efifl U ”U; k k It I m 1,771 1,. 3L1”) , + 2 Ci m ave m (-wf.(du + ui) + fk.) d“): dvi (4.18) l,m k . .. , ,, 2 1. 1, (”Mam 0 = [6 e 6(‘112—‘ltz‘pyvz —- “1:21) a), Z CimLIl/Il’m — f,- + 2: Di "‘(dv + Up—Bue U‘k‘vk l,m l,m k 8L3! + Z Dg'ln—Elél(—wg(du +212.) + fk) (in?C (it); (4.19) l,m for i=1,2,3,...,Q; 1'75 k; _ T . P— 13'“) 1).“) q :1) "'9 Arne) 135 where (ufwwz‘q) E [—ub,ub] x[—vb,vb] are the collocation points, which are zeroes of TNmei/“bl and TNpJ,(vf:,/vb), respectively (see [78], [81]). Equations (4.18) and (4.19) constitute a set of 2(Q — 1) x Np,” X Np‘v nonlinear equations in the C’s and D’s. Note that there are N; x NS sets of Cs and D’s, that is, a set for each patch. However, if the system is conservative and the nonlinear terms are functions of solely the modal positions, then only 0.25 x N3 x N5 sets of C’s and D’s need to be solved for, and the remaining coefficients can be generated using symmetries in the NNM manifold. Once all of the expansion coefficients are obtained, the X ‘s and Y’s, which describe the slaved modes, are known functions of the master states ((12.1%). For 2' = k, these known functions are used to express fk in equation (4.2) in terms of only 21E. and vi, rendering a single-DOF oscillator as the reduced-order—model for the km N N M. Once all of the coefficients have been determined for each subdomain, the entire manifold is known, in a. patchwork form. It should be noted that a motion started on this manifold, that is, one that satisfies the master-slave relations captured by the X’s and Y’s, remains on that manifold for all time. by dynamic invariance. However, a given response will visit many subdomains, and thus the solution moves across this two dimensional surface like a curve tracing across a tiled floor. It should also be noted that these “tiles” are obtained in a patchwork manner, without any continuity or smoothness constraints, so that the slave constraints may exhibit jumps as the response moves across the tile boundaries. If the subdomains are taken to be. sufficiently small, these jumps will likewise be small. Examples of nonlinear modes constructed from a similar approach can be found in [33] and [34]. 4.3 Fixed-Interface Nonlinear Component Mode Syn- thesis In this section, we extend the concept of the fixed—interface CMS of Craig and Bamp- ton ([62]) by making use of the fixed-interface NNMS instead of the fixed-interface linear 136 normal modes (LNMs). The aim here is to account for nonlinearities at the substructure level without having to resort to the retention of several linear modes. First, the nonlin- ear structural system is partitioned into component substructures. Then, fixed-interface dynamic representations of the substructures are constructed, where the nonlinear models are expressed in terms of fixed-interface LNMS, which are coupled through nonlinear terms. These substructure models are then reduced, each to a single nonlinear mode, using the NNM constraint relations. and subsequently synthesized with linear constraint modes to produce a ROM that describes the dynamics of the combined structure. 4.3.1 System Representation in Physical Coordinates Consider a structure that is partitioned into two substructures, denoted by a and B. The equations of motion of the substructures in physical coordinates can be written separately, along with appropriate conditions on the common junction coordinates, as follows: 0 .- O’ O 0 M11 MIJ XI K11 KIJ XI MJI MJJ Jh KJI KJJ XJ (4.20) 0 ‘]Q 0100 _ 0 GJ(X) FJ ,1} .. ,3 l3 )3 M11 MIJ XI + - K11 KIJ XI MJI MJJ XJ [KJI KJJ XJ (4.21) 13 ,3 GI(X) _‘ 0 GJ(X) FJ 137 with junction boundary conditions on displacements and forces, as follows, X3 X’B, (4.22) F9+Ff : o, (4.23) where ()0 denotes quantities associated with substructure a, and all terms defined below have ()3 counterparts. The vector X 0: contains all coordinates for the a substructure and has dimension N a, 3’ is the vector of interior (non-interface) physical coordinates of dimension 71?, and X“; is the vector of junction (interface) physical coordinates of dimension nC. G? is a vector of the nonlinear forces associated with the X? equations, 3‘ is the vector of nonlinear forces associated with the X3" equations, and F? is the vector of reaction forces from substructure ,8 acting on substructure (1. Mg], Kf’I, etc., are mass and stiffness matrices, respectively, defined in an obvious manner. 4.3.2 System Representation in Linear Modal Coordinates In this section, the (typically) large number of physical coordinates for each substructure is first reduced to a smaller (yet still possibly large) number of coordinates by truncating using the substructures’ fixed-interface linear modes. The physical coordinates of substructures a and x3 are transformed into a truncated set of fixed—interface linear modal coordinates using the following coordinate transformations: 0: Xa : we: "a : l: (pa ‘1’0 :l "N (4.24) N C a 17C 138 i0 nN a- a a- B a X _\11 n _ [‘I’N QC] "[3 , (4.25) where (Pg, is the. 1 ’0 x n0 matrix of retained fixed-interface linear normal modes, 1710‘; is the vector of retained fixed-interface linear modal coordinates (of dimension n“ ), \I'g. is the N a X "C matrix of linear constraint modes, 778. is the vector (of dimension TLC ) of linear constraint. modal coordinates, and similar terms are defined for substructure 3. The procedures for calculating the retained fixed-interface linear normal modes and the linear constraint modes are described below. Note that at this step, the dimension of the system of equations for substructure a has been reduced from (N0 = n? + 72C) to (n0 + mg) with n” < 71‘}, and similarly for substructure 3 4.3.2.1 Fixed—Interface Linear and Nonlinear Normal Modes The concept of fixed-interface linear CMS of Craig and Bampton ([62]) is extended by making use of the fixed-interface NNMS in place of the fixed-interface LNMS. Since the procedures for substructures a and ,3 are the same, these superscripts will be omitted from the development here. Consider equations (4.20) and (4.21) with XJ = 0, that is, with the interface fixed. This yields MIIXI + KIIXI + GI(XIs XJ = 0) = 0. (426) Using standard linear modal analysis (with the modes obtained for G I = 0), we obtain the 139 normalized fixed~interface linear modal matrix <1) I N associated with the interior physical coordinates (of dimension n] X n, where n << 71,). Then, the retained fixed-interface linear modal matrix N is defined as ‘I’N = . (4.27) The nonlinear equations of motion (4.26) are then transformed into fixed-interface linear modal coordinates using the following coordinate transformation: XI = ‘I’IN 771v. (428) This yields the fixed interface nonlinear equations of motion expressed in fixed—interface linear modal coordinates. as follows: INN fiN+ANN 77N+GN(77N) = 0 (429) where INN = @fNMuem (4.30) ANN = @fNKHem (4.31) éNWN) = @fNGfl‘PINnN’XJzol- (432) This fixed-interface model is still relatively large and is to be ultimately reduced to a single DOF using the nonlinear modal reduction procedure described in sections 4.2 and 4.2.2. Note that in terms of matching notation from those sections, the function G N and 140 dimension n are — f and Q, respectively. The NNM reduction provides the constraint relations between the slaved states and the master states of the kth fixed-interface NNM for a substructure, and can be expressed as follows: n.” X471? .77.? ) (4.33) , r N , r n? = Yzlnk .721?) (4.34) for i=1,2,3,...,n; iaé k. After obtaining the constraints X,- and Yi, and doing so for each substructure, we have single mode nonlinear models for the substructures, based on fixed interface models. We now turn to the procedure for coupling these component structure models to one another through the interface. 4.3.2.2 Linear Constraint Modes The definition of a linear constraint mode from Craig and Bampton [62] is adopted here. It is defined by statically imposing a unit displacement on one physical coordinate of the set of junction coordinates and zero displacements on the remaining coordinates of the set. This procedure is applied consecutively to all junction coordinates. Again, this process is the same for the a and t3 substructures, and therefore these superscripts are omitted. The collection of linear constraint modes is obtained from we 2 = , (4.35) 141 where ‘I’IC has dimension n I X TLC and I JC has dimension nC X 720. Note that nonlinear constraint modes are not considered in the present study, for a few rea- sons. First. to do so would require additional steps involving the derivation of another set of master-slave relations to describe nonlinear constraint modes (even if they are static). Also, if one were to use nonlinear constraint modes in the synthesis procedure, the coefficients of the nonlinear terms would be functions of the constraint mode amplitudes, which would make the synthesis process cumbersome to the point of being impractical. Furthermore, and perhaps most importantly, since fixed-interface methods work best for assemblies of weakly coupled subsystems, for which fixed—interface (linear or nonlinear, as the case may be) modes are well-suited to describe component motions accurately, it is appropriate to use linear constraint modes to represent the small interface-induced motions in the com- ponents. This restriction must be kept in mind, however, since it limits the method to a class of systems with weak coupling. Equations (4.20) and (4.21) are transformed into linear modal coordinates using equations (4.24) and (4.25). The equations of motion for both substructures expressed in terms of these coordinates have the form (1 .0 C! a INN MNC 771v ' + ANN 0 UN ' MCN MCC fie 0 K00 ’10 (4.36) ~ 0 O GN(77Na770) __ 0 (30(77ch) ' F0 142 INN MNC fiN + ANN 0 7m MON MCC fie 0 K00 no (4.37) ~ B X3 GN(77N, 77C) _ 0 éC(’7Na 770) PC with the following junction boundary conditions on the displacements and forces: X31 = ng 2 n3 -_- Xfi, (4.38) Fg = Ff; and pg = Ff,’ (4.39) Fg + pg 2 0. (4.40) Note that the stiffness matrices are block diagonal due to the orthogonality between <1) N and WC with respect to K, which results from how (I) N and ‘110 are defined in equations (4.27) and (4.35), respectively. We now turn to the synthesis procedure. 4.3.3 Synthesis with Nonlinear Modal Reduction In this section, the model size of the substructures is reduced using the NNM constraint relations. The substructure ROMS are then synthesized to obtain the final ROM that describes the dynamics of the combined structure. In summary, this nonlinear CMS process involves two steps of model reduction, one of which is standard and uses linear modes, while the other involves the NNMS. First, we 143 use the linear modal reduction described in section 4.3.2 to reduce the system DOF from (n? + n? + ac) to (n0 + n3 + 71.0). Note that this reduction is not central to the technique, and may not be used in some systems (in fact, it is not done in the examples considered in the next section). Then, the key step makes use of the nonlinear modal reduction based on the approach described in sections 4.2 and 4.2.2, which reduces the system size from (n0 + 72:3 + 72.0) to (2 + 720) DOFs. Again, the superscripts a and ,3 are omitted in the derivation. The component forms of equations (4.36) and (4.37) are given by n 1 ..N .NC--C 2 N ~N N N C 77,- + mis 713 +w, 77,: +9,- (77! ,7]k ,177.) = O (4.41) 3:1 ”C ~N ’VCuC 2 N ~N N N C m, +2774... Us +wkm, +91. (721 .77.. .771) = 0 (4-42) 321 n "C ”C CV~V CC~C CC C ~C N N C J 7”er ”3' +anrs Us +Zkrs "ls +9r (771 ’77k inr) = Fri (4'43) j=1 3-1 3:1 for i,l = 1,2,3,...,n; i,l;ék; r = 1, 2, 3, ..., nC, where k is the kth fixed—interface master mode defined in section 4.3.2.1 (typically taken to be the fundamental mode). From equation (4.41), the expression for 77;“, is obtained , and it is used in equation (4.43). Rearranging and grouping terms in equation (4.43) yields 144 ”C CC I I .. ~C 1 t N N mfiN ijk+ 3: 771,3 “pm e 7),?' + 821 kCC ”SC + 9 up” 8(7), ,nk ,0?) = F,:](4.44) where mgc‘wdate = 771C? + mac , (4.45) Tl 777g.C * = Z —777,C.;-N m 79:0, (4.46) j=ld¢k g?"“”"““’ = a? + a? (4.47) n and 99* = Z 77123] (‘sz 7)]N — 9?), (4.48) j=ld¥k for r : 1,2,3,...,n The modal constraints given in equation (4.33) are used in equations (4.42) and (4.44) to reduce the equations of motion associated with the synthesized fixed-interface linear modal coordinates (which are coupled to the constraint modes) onto the km fixed-interface NN M invariant manifold (which is uncoupled from the constraint modes). This approximation will yield an accurate ROM if substructures a and [3 are weakly coupled, that is, the actual NN M invariant manifold of dimension 2 X (2+nC) of the combined structure has invariant 2- kth dimensional submanifolds that are close to the fixed-interface NNM invariant manifolds of the substructures, which are each of dimension 2. With the NNM reduction of substructures a and ,3, equations (4.42) and (4.44) can be expressed in matrix form as O (1 (1 r O 1 MKC ”iv “1% 0 n? d t .. MCK M31}; 0 e "C 0 K00 ”C (4.49) (JIM/[l]. ‘I'ij. 770) _ O .. d t 7V — GE? (1 30)}? ,77iV "0) F] 145 1 MKC 77,13’ + ...,: 0 77,1," d t .. MCK Mg? 0 e no 0 K00 no (4.50) ~ N V 3 fl Gk(nk 777; an) __ 0 d t '— ’ Cup ”(77,, 77,. We) FJ along with the junction boundary conditions from equations (4.38) and (4.40). The final form of the synthesized ROM is obtained by imposing the junction conditions, equations (4.38) and (4.40), on equations (4.49) and (4.50), resulting in r 1 0 Mgc - ‘77]“0‘1 0 1 M13“; fifcv’fi MCK MgK M32510“, L "C _ , (4.51) "(4,2 0 0 q Fm?“- [ 53 - ’0- + 0 (of)? 0 77,1)” + (7]: = 0 _ 0 0 K00] _nc _ _c':‘g’d“‘8_ _0J where Mupdate = (Mupdatea+Mupdate,fi) (4-52) KCC = (KgC+KgC), (4.53) and ézépdate : Gupdate,a(nllcV,a,77i:V,a,nC) +G“pd“‘e’f’> “’1 , “’1 . Applying the procedures from section 4.3.2.1, the fixed-interface nonlinear mode manifolds for the fundamental modes of substructures fl and 7 are obtained. In this simple case the manifolds can be expressed as constraint functions in which the second linear modal position and velocity for a substructure depend on the first modal position and velocity of that substructure. The results obtained numerically for the modal position functions for the two multi-DOF substructures are shown in figures (4.3) and (4.4). For substructure 6, the boundaries of the surface along modal position and velocity are Ub = 2.6 and Vb = 1.79, respectively, the number of Chebyshev polynomials along each u": and U]: are NP,“ 2 3 and NP,” 2 2 respectively, and the number of pieces along uk and vk are N5 = 100 and Nf,’ = 100, respectively. For substructure 7, the boundaries of the surface along modal position and velocity are Ub = 12.0, Vb = 16.78, respectively, the number of Chebyshev polynomials along 21 k and along 17], are Np,“ 2: 3 and NW) = 2, respectively, and the number of pieces along uk and U), are N5 = 100 and N5 = 100, respectively. Note that similar surfaces exist for the slaved modal velocities, although they are not shown here. Applying the procedures from section 4.3.3, a nonlinear three-DOF ROM is obtained. In this case it describes the motions of the coupling mass and the first NNM for each of the attached substructures. This model is simulated with initial conditions initiated on the first fixed-interface nonlinear mode manifolds of substructures B and 7, as follows: N , 73 7)] (0) 2: 1.5,77‘1N’7(0) = 0.4, along with a. constraint amplitude of 770(0) = 0.1, with initial modal velocities taken to be zero. For comparison purposes we use four models constructed in various manners. They are: the original five-DOF nonlinear model given in equation (4.55), a three-DOF nonlinear 150 model obtained using linear CMS (that is, by projecting the nonlinear subsystems onto their fixed-interface linear modes and the linear constraint mode), a three-DOF linear model obtained using linear CMS (this is the linearized version of the previous model), and the threeDOF nonlinear CMS model developed by the proposed approach. Figures (4.5) - (4.9) depict the comparison of the time responses of the various physical degrees of freedom obtained by simulations of these four models. In these figures the dashed line labeled NLCMS corresponds to the three-DOF nonlinear CMS model, the dash-dot line labeled LCMS corresponds to the three-DOF nonlinear model obtained using linear CMS, the dotted line labeled LCMSLPrb corresponds to the three-DOF linear model obtained using linear CMS, and the solid line labeled Original Model corresponds to the five-DOF original model. From the simulation results it is clear that the fixed-interface nonlinear CMS model outperforms the other low-order models, including the linear CMS model, which is the approach most commonly used in this type of problems. In particular, the NLCMS correctly captures the frequencies of the substructures, as well as the effects of their coupling. It should be noted that there is a range of system parameters where this method gives satisfactory results. One particular trend worth noting is obtained by varying the stiffness to ground k1. As k1 is reduced and approaches the same order as the stiffnesses of the substructures, the fixed-interface approximation begins to break down and the accuracy of the NLCMS model deteriorates. On the other hand, as [$1 is increased, the NLCMS model becomes more accurate, but the amount of substructure interaction is reduced, leading to a system which consists of essentially decoupled motions for the two attached substructures. 4.5 A Class of Systems In order to illustrate the general applicability of the proposed method, a class of nonlinear spring-mass systems is examined. This class of systems is schematically depicted in Figure 4.10, and naturally partitions into Na appendage substructures, as shown in Figure 4.11. 151 This system has the same general features as the system in section 4.4 except that each substructure has (Nmasg + 1) DOFs. (Note that this is easily generalized to the case where the substructures have different number of DOF). In the following sections, the model of this class of systems is first developed and then a system having Na 2 2 and Nmass = 20 is used as an example. 4.5. 1 Model Development The component forms of the equations of motion of the N = (Na Nmass + 1)—DOF system are given by N N Z "LU .I‘j + Z kij .Tj ‘l' gi(X) : 0 (4.62) j=l j=1 for '= 1, 2, 3, ..., N. The formulae for mij. Isl-j and 9,- can be found in the Appendix. The component forms of the equations of motion for the individual substructures and the junction boundary conditions are given by -r . . 5 m1 .1)? + k1 1'? = 2 F? (4.63) SI;’3.‘7.5,... (Nmass +1) (Nnmss +1) 615.65 + Z 165. 1:53 + g;9(XS)=f.S (4.64) j=1 i=1 for i: 1, 2, 3, ..., (Nmass +1), and S 2,13, 7, 6, with junction boundary conditions on displacements and forces: 1‘? = .17? = :17? = 51:31; 2 (4.65) F?” +Ff“ = Ff,” +F}"’ 2 F35 +F§a = = 0. (4.66) The formulae for mg, kg, 925‘, and fl-5 can be found in the Appendix. For bookkeeping of the indices, the procedures from section 4.3.2.1 are slightly modified. Instead of applying the coordinate transformation from equation (4.28) to equation (4.26), the coordinate transformation X = q’N'IN is applied to equations (4.20) and (4.21) with the interface fixed (X J = 0); this step also lead to equation (4.29). Therefore the nonlinear vectors of the fixed-interface nonlinear equations of motion, ex- pressed in fixed—interface linear modal coordinates of substructures S, are given in their component forms by 715 ”S n5 - E :j :E : N,S, N,S N,S N,S p=1q=pr=q for 2': 1, 2, 3, 723 and S = i3, 7, 6, ..., where the formulae for (2327:” can be found in the Appendix. Then the fixed-interface non— linear mode manifolds for the fundamental modes of substructures ,B, 7, are computed. Next, the linear constraint modes of the substructures S are computed. Equation (4.64) is then transformed into linear modal coordinates using the coordinate transformation defined in section 4.3.2. The nonlinear vectors of the nonlinear equations of motion (expressed in linear modal coordinates of substructures S) are given in their component forms by S (115+nC) (nS+nC) (11 +110) - Z Z Z 5,... ])=1 (1:!) 7:0 for 1:1,293a”°3(n5 + 71C) and S = ,3, ’7, 6, The formulae for big], can be found in the Appendix. 4.5.2 Example As an example we consider a system with Na 2 2 and Nnmss = 20. In the following sections, the fundamental fixed-interface NNMS of substructures [3 and 7 are first 154 developed. Then, by applying the procedures described in section 4.3.3, a three-DOF nonlinear CMS model is obtained which describes the motion of the coupling mass and the first NNM for each of the attached substructures. For comparison purposes, four other models are also developed. They are: the forty-one—DOF nonlinear original model given in equation (4.62), a forty-one—DOF nonlinear model obtained by using linear CMS that retains all the fixed—interface LNMs from individual appendage substructures fl and 7 (71/3 = 20, n7 = 20), a three-DOF nonlinear model obtained by using linear CMS that retains only the first fixed-interface LN M from individual appendage substructures fl and 7 (n5 = 1, 127’ = 1), and a three-DOF linear model obtained by using linear CMS (this is the linearized version of the previous model). 4.5.2.1 Substructure ,3 In this study, the system parameters used to construct the fixed-interface NNM of sub- structure 6 are first chosen such that the system is similar to the equal—length FE model of a rod under axial vibration, i.e., each block has the same mass and each spring has the same spring constant. For the discrete system with equal parameters and under axial vi- bration, the ratios of the higher linear natural frequencies to the fundamental linear natural frequency are closed to 3:1 (mode 2 to mode 1), 5:1 (mode 3 to mode 1), 7:1 (mode 7 to mode 1), etc., which are the ratios of the continuous system under axial vibration ( [90]). Since the nonlinearities present in the system are cubic polynomials, the system will expe- rience internal resonances. Moreover, the system parameters chosen in this way will result in weak nonlinear couplings between the fundamental linear mode and the higher modes, which causes the fundamental nonlinear mode shape to be close to the fundamental linear mode shape. To remove internal resonances from the system, and to generate stronger nonlinear coupling and modal distortion, the linear spring that connects the subsystem to ground is taken to be softer than the other springs. This causes the fundamental nonlin- ear mode shape to be very different from the fundamental linear mode shape at the large amplitudes. Some details of this result follow. 155 For the first trial, the system parameters used to construct the fixed-interface NNM of substructure 6 are 7712 2 m3 2 m4 = : ”7.21 = 1, k2 2 k3 2 k4 = 2 k4] = 1000. The number of retained fixed-interface LN Ms (n3) is chosen to be 20. The first four fixed- interface (linear) natural frequencies are a)? = 2.42 rad/s, tag 2 7.25 rad/s, w? = 12.04 rad/s, cuff = 16.76 rad/s. The first four fixed-interface LN Ms are shown in Figure 4.12. Figure 4.13 depicts coefficients of the pqr cubic nonlinear terms expressed in fixed-interface linear modal coordinates and associated with the first mode (master mode). Note that there are 1540 nonlinear terms ( (n!3 x (n5 + 1) x (nfl + 2)) / 6 ) for each retained fixed-interface LNM. The horizontal axis corresponds to the nth term (i.e. the first term is (17?“fi )3 , the second term is (wy‘lw’fi)2 n‘zw’fl,” the last term is (n%’fi)3). Note that there are 20 envelopes in the figure, and each envelope corresponds to a value of the first index p. We can see that most of the coefficients are zero, which means that the first mode is not strongly coupled with higher retained fixed—interface LNMs, and therefore they will not be reflected in the N NM. The first fixed-interface NN M manifold can be depicted as constraint surfaces depending on the first modal position and velocity. A total of thirty-eight such constraint surfaces are needed to describe motions on the first fixed-interface NNM manifold in the 40-dimensional state space. Samples of such surfaces are shown in Figures 4.14, 4.15, and 4.16, which cor- respond to the contribution from the second, third, and fourth fixed-interface linear mode amplitudes to the first fixed-interface NN M manifold. The boundaries of these surfaces along modal position and velocity are Ub = 40.0 and Vb = 96.9, respectively. The sur- faces are obtained by using Chebyshev polynomials along ui and v}: with Np,” = 3 and NW, 2 2, respectively, and the number of pieces along uk and vk to be N5 = 50 and N5 = 50, respectively. Figure 4.17 depicts the time responses of the first fixed-interface mode displacement ob— 156 tained from two different models. The dashed line is the time response of the ROM defined N fl at the end of section 4.2.2 with initial modal displacement 771 (0) = 33.9 and velocity 01“3 (0) = 0. The solid line is the time response of the original model (20 DOFs) initiated on the fixed-interface N NM manifold. Figures 4.18, 4.19, and 4.20 depict the time responses of the second, third, and fourth, respectively, fixed-interface linear mode displacement ob- tained from the ROM and the original model. The dashed line is obtained by applying the master-slave constraints (depicted in Figures 4.14, 4.15, and 4.16) to the numerical solution of the ROM (the dashed line in Figure 4.17). The solid line is obtained from the 20-DOF original model. The time responses of the ROM and the original model are almost identi- cal at the beginning, and as time progress they begin to diverge. From the linear natural frequencies, we can observe that rug z 3014}, w? x 2h)? —wf , and a)? z a)? +0)? —w1fi, which indicates that there are internal resonances (to the first-order perturbation approximation ([79])) among these modes in the fixed-interface system of substructure 3. Therefore, en- ergy is being exchanged among the modes involved in the internal resonances, which cause the motion of the original model initiated on (in practice, very near) the first fixed-interface NNM to eventually leave the manifold. In order to obtain the ROM that embeds the in- ternal resonances, a multi-mode invariant manifold formulation is needed, as described in [29], [71], and [42]. The multi-mode invariant manifold formulation is beyond the scope of the current study. From Figures 4.17, 4.18, 4.19, and 4.20, we can see that the contribution of the higher modes to the first fixed-interface NNM is less than 3% of the initial modal displacement nifV’/3(0) 2: 33.9. Hence we expect that the one-mode nonlinear model, which is obtained by projecting the fixed-interface subsystem onto the first fixed—interface LNM, will yield a good time response, at least for a limited time interval. Figure 4.21 depicts the time responses of the first fixed-interface modal displacement ob- tained from the linearized one-mode model, the one-mode nonlinear model, and the ROM. The dotted line is obtained from the linearized one-mode model, which is the linearized version of the one-mode nonlinear model. The dashed line is obtained from the one-mode nonlinear model, which is acquired by projecting the fixed-interface subsystem onto the first fixed-interface linear mode. The solid line is obtained from the NNM ROM. We can 157 see that the time response of the one mode nonlinear model is close to the time response of the ROM for many periods. This is due to the fact that the first fixed-interface linear mode is not strongly coupled with the higher retained fixed-interface LNMs, as mentioned above. To eliminate the internal resonances from the fixed-interface subsystem and to yield strong coupling between the first fixed-interface linear mode and higher linear modes, a new set of system parameters are introduced. These will make enhance the nonlinear aspects of the problem to more fully demonstrate the accuracy of the proposed approach. By observ- f A ’S and bN’S mm. ,1qu from the Appendix and Figure 4.12, we can see that the in g coefficients b displacements of DOFs near the fixed-interface of the low-frequency linear mode shapes N,x3,* 1qu coefficients to be small, or zero. In order to are small, which cause many of b yield stronger coupling between the first fixed-interface linear mode and higher modes, a new set of system parameters are considered: m2 = m3 2 m4 = = mm = 1, k2 = 10, k3 = 164 =2 = k4] = 1000, which allows DOF near the interface to have significant displacements, in particular for the fundamental mode. Twenty fixed-interface linear natural frequencies are presented in Table 4.1, and it is seen that they do not possess low—order internal resonances among them. Figure 4.22 depicts the first four fixed-interface LN Ms obtained from using the second set of parameters. Figure 4.23 depicts the coefficients of the pqr cubic nonlinear terms expressed in fixed—interface linear modal coordinates and associated with the first mode (master mode), obtained using these parameters. Comparing Figure 4.23 with Figure 4.13 indicates that the first fixed- interface linear mode strongly couples with higher modes for this set of parameters. Note that for these parameters the first fixed-interface linear mode is similar to a rigid-body mode. Therefore one might argue that the spring next to the wall is soft and the other springs are very stiff, which would imply that the fixed-interface subsystem would behave like a one-DOF lumped-mass system. Therefore the one-mode nonlinear model obtained by projecting the fixed-interface subsystem onto the first fixed-interface linear mode would be 158 adequate to capture the dynamics of the subsystem. However, Figure 4.23 indicates that the first fixed-interface linear mode is strongly coupled with higher retained fixed-interface LNMs, which indicates that such an argument is not justified. Hence, the concept of NN M is needed in order to capture the correct dynamics in the large amplitude regime, which is illustrated in following results. Figures 4.24, 4.25, and 4.26 depict the contributions from the second, third, and fourth fixed-interface linear mode amplitudes to the first fixed-interface NNM manifold. The surfaces are obtained by using Ub = 3.0, V5 = 2.06, NW, = 3, Np,” = 2, N5 = 40 and N5 = 40. Note that the maximum contribution from the second linear mode amplitudes to the first NN M manifold is around 10%, which indicates significant coupling between the first and second modes. Figure 4.27 depicts the time responses of the first fixed-interface mode displacement ob- tained from four different models. The dashed-dotted line, the dotted line, and the dashed line are the time responses of the linearized one-mode model, the one-mode nonlinear model, and the ROM (based on invariant manifolds), respectively, with initial modal displacement TAM/3(0) = 1.20 and velocity fiiN’fiW) = 0. The solid line is the time response of the original model with initial conditions initiated on the fixed-interface NN M manifold. The time responses from the ROM and the original model are nearly identical. Figures 4.28, 4.29, and 4.30 depict the time responses of the second, third, and fourth, respectively, fixed-interface linear mode displacements obtained from the ROM and the original model. The dashed and solid lines are the time responses of the ROM and the original model, respectively. The time responses from the ROM are close to the original model, and the errors are very small when compared with the amplitude of the overall response (> 1). (From these figures, there exist discontinuities in the time responses since continuity between patches of the constraint surfaces is not enforced. However if the subdomains are taken to be sufficiently small, these jumps will be small.) 159 Figures 4.31 and 4.32 depict the time responses of the physical displacements 3’26 and 56%, respectively, (the leftmost masses of each substructure) with the interface fixed for the lin- earized one-mode model, the one-mode model, the ROM, and the original model. The time responses from the ROM and the original model are nearly identical. The time responses from the linearized one-mode model and the one—mode nonlinear model are significantly different from the original model, since the contribution from the higher fixed-interface linear modes is not accounted for in those two models. In other words, the motions of those two models occur on a flat approximation of the invariant manifold (i.e., the linear eigenspace), which is inadequate in the present case. Figure 4.33 depicts the comparison of deflections of all physical coordinates as? (t) at three different times: t = 0, the initial time, where the system is at its peak configuration; t = 7% , one eighth of a period of the nonlinear system; and t = 1:; , a quarter of the period, at which the system is at zero displacement (with nonzero velocity, of course). Two cases are shown: the first fixed-interface LNM and the first fixed-interface NNM. Note that the shape of deflection on the linear mode manifold (linear mode shape) does not change with time, however the deflection shape on the nonlinear mode manifold (that is, the nonlinear mode shape) does change with time, due to the time-dependent contributions from the other retained linear modes, which is significant for this case. The distortion of the NNM is evident. Note that the computational time required to obtain the manifold for this case (with N5 = 40 and N5 = 40)is about 30 minutes. The computational time associated with the simulation of the twenty-one—DOF original model is about 13 seconds. The computational time associated with the simulation of the ROM is about 2 seconds. Thus, if one is to carry out many simulations, the generation of the ROM is worthwhile. 4.5.2.2 Substructure 7 The system parameters used for substructure 7 are mg? 2 m23 = 777.24 = = 77141 = 1, 160 [£42 2 42.5, 1643 = [€44 = = 1681 = 4250. This yields strong nonlinear coupling between the linear modes, just as was done for substructure 6. The results that are analogous to Table 4.1 and Figures 4.22 to 4.33 for substructure 6 are shown in Table 4.2 and Figures 4.34 to 4.45. To avoid redundancy, detailed descriptions are not repeated here. Note that since the ratio among the linear and nonlinear spring coefficients are the same as for substructure 3, the surfaces shown in Figures 4.36 to 4.38, obtained by using Ub = 3.0, Vb = 4.24, NW, 2 3, Np,” = 2, N5 = 40, N5 = 40, are similar to those shown in Figures 4.24 to 4.26. However, the fixed-interface linear-mode natural frequencies of substructure 7 are higher than those of substructure 6. 4.5.2.3 Synthesized Structure The system parameters for the base structure (substructure a) are ml 2 1,500 and k1 = 24,000, which yields Lac 2 4.0, which is close to the fundamental frequencies of the substructures. However, since this coupling structure is very stiff, its response will remain small, consistent with the assumptions needed to employ linear constraint modes. A comparison of the linear natural frequencies of the three-DOF linear-CMS nonlinear model and the three-DOF nonlinear-CMS model is shown in Table 4.3. Comparison of the first sixteen Linear natural frequencies of the forty-one-DOF linear-CMS nonlinear model and the forty-one-DOF original model are shown in Table 4.4. These two comparisons are used as the first check to verify that the synthesized models are correct. The linear natural frequencies from Table 4.3 are in agreement with those from Table 4.4, at least up to the third decimal digit. Note that the forty—one—DOF linear-CMS nonlinear model is equivalent to the forty-one-DOF original model since the modal matrices from substructures fl and 7 together represent a complete set of forty-one bases that are used to describe the original model in the synthesized linear-modal coordinates. Figure 4.46 depicts the time responses of the first synthesized fixed-interface mode dis- 161 placement "(V213 obtained from four different models. The dashed-dotted line, the dotted line, and the dashed line are the time responses of the three-DOF linear-CMS linear model, the three-DOF linear-CMS nonlinear model, and the three-DOF nonlinear-CMS model, re- spectively, with initial modal displacements niw’fim) = 1.2, nil/“7(0) = 1.2, 770(0) = 0.03 and modal velocities fjiw’fi (0) = 0, nil/”(0) = 0, 770(0) 2 0. The solid line is the time response of the forty-one-DOF linear-CMS nonlinear model with initial conditions initiated on the fixed-interface N NM manifolds of substructures 6 and 7 (this is the full model). We can see that the time responses of the three-DOF nonlinear-CMS model and the forty-one— DOF linear-CMS nonlinear model are nearly identical over the time interval shown, and the other two models show significant errors, especially the fully linear model. Figures 4.47, 4.48, and 4.49 depict the time responses of the second, third, and fourth, respectively, synthesized fixed-interface linear mode displacements (77$V ’fl , "11;! ’fi , and nfiv’fi), respectively, obtained from the three-DOF nonlinear-CMS model and the forty—one-DOF linear-CMS nonlinear model. The dashed and solid lines are the time responses of the three-DOF nonlinear-CMS model and the forty-one—DOF linear—CMS nonlinear model, re- spectively. They are nearly identical except that the forty-one—DOF linear-CMS nonlinear model has some high frequency “ringing”. This is because the initial conditions for the forty-one-DOF linear-CMS nonlinear model are initiated on the fixed-interface NNM man- ifolds of substructures .13 and 7, which are only an approximation of the actual 6-dimensional N NM invariant manifold of the whole structure. Figures 4.50, 4.51, 4.52, and 4.53 are analogous to Figures 4.46, 4.47, 4.48, and 4.49. Therefore, the descriptions similar to those given for substructure 13 are not repeated here. The comments related to those analogous figures are also applicable to substructure 7. Note that since the first fixed-interface nonlinear-mode natural frequency of substructure 7 is higher than that of substructure 6, the time windows in Figures 4.50, 4.51, 4.52, and 4.53 are adjusted in order to depict clear comparisons. Figures 4.54, 4.55, 4.56, 4.57, and 4.58 depict the time responses of masses 7721, 171-2, 77121, 162 77122, and 777.41, respectively, obtained from four different models ( these are the coupling mass and the end masses of each substructure). The dashed-dotted lines and the dotted lines are the time responses of the three-DOF linear-CMS linear model and the three-DOF linear-CMS nonlinear model, respectively, obtained by applying the coordinate transforma- tions defined in section 4.3.2 to the time responses of the synthesized fixed—interface linear modal coordinates of both models with nff = 1, n”7 = 1, and 110 = 1. The dashed lines are the time responses of the three-DOF nonlinear-CMS model obtained by applying the coordinate transformations defined in section 4.3.2 to the time responses of the synthesized fixed-interface linear modal coordinates having n3 = 20, n7 = 20, and no = 1. The solid lines are the time responses of the forty-one-DOF original model with initial condi- tions initiated on the fixed-interface NNM manifolds of substructures fl and 7. It is seen that the time responses of the three—DOF nonlinear-CMS model and the forty-one—DOF original model are nearly identical. However, the time responses of the three-DOF linear- CMS linear model and the three-DOF linear-CMS nonlinear model are quite different from the time responses of the forty-one-DOF original model. This is due to the fact that the contributions from the higher fixed-interface linear modes are not accounted for in those two models. In other words, the motions of those two models occur on the flat approx- imation of the invariant manifold of dimension 6, which is not a good approximation in this case. (Note that such an approximation would be quite good for the original set of parameters. for which the modal coupling was small.) Figures 4.59 and 4.60 depict the time responses of the first synthesized fixed—interface fl and "4171M, respectively, obtained from the three-DOF nonlinear- mode displacements niv’ CMS model and the forty-one—DOF linear-CMS nonlinear model. The dashed lines are the time responses of the three-DOF nonlinear-CMS model with initial modal displace- ments niv’flm) = 1.3, ni\r’7(0) = 1.3, nC(0) = 0.04 and modal velocities Trim/3(0) :- 0, 7'11V’7(0) = 0, 7")C(0) : 0. The solid lines are the time response of the forty-one—DOF linear-CMS nonlinear model with initial conditions initiated on the fixed-interface NNM manifolds of substructures B and 7. It is seen that the time responses of the forty-one-DOF linear-CMS nonlinear model starts to deviate from the three-DOF i‘ionlinear-CMS model 163 after several oscillations. This indicates that at these initial conditions, with the base- structure parameters m1 = 1,500 and k1 = 24, 000, the actual NNM invariant manifold of the whole structure deviates from the fixed-interface NN M manifolds of substructures ,8 and 7. Note that we also performed a numerical experiment by fixing the ratio between m1 and k1 and reducing their values. It was found that the valid time domains of the fixed- interface N NM manifolds of substructures ,8 and 7 are reduced as well. This is expected to occur, since this parameter trend allows more interaction between substructures fl and 7. This causes the actual N N M invariant manifold to deviate further from the fixed-interface NN M manifolds of the individual substructures, especially for large-amplitude motions. Finally, the computational time associated with the simulation of the three-DOF nonlinear- CMS model is compared with that of the forty-one—DOF linear-CMS nonlinear model. It takes about 180 seconds to simulate the forty—one-DOF linear-CMS nonlinear model. However it takes only 60 seconds to simulate the three-DOF nonlinear-CMS model, which is a factor of three lower. The primary reason for this is that a smaller time step (a factor of 2.5 smaller) is needed for the ode-solver to integrate the stiff forty-one-DOF linear-CMS nonlinear model. Therefore, in terms of computational time associated with the simulation, the three-DOF nonlinear-CMS model does not have much advantage over the forty-one- DOF linear—CMS nonlinear model, as we might expect. This is mainly due to the fact that one is required to compute the nonlinear update terms gg’uPdate in equation (4.44) of both substructures, which represent nonlinear terms that arise from the slave fixed- interface linear-modal-coordinate equations of motion. This is the same scenario as for the forty-one-DOF linear-CMS nonlinear model, since all nonlinear terms in this model need to be computed. This situation is in contrast to the individual fixed-interface substructures. From section 4.5.2.1, it can be seen that the simulation time of the ROM is a factor of 6.5 smaller than the simulation time of the twenty-DOF original model. This is because the nonlinear terms of the master mode EOM are only computed in order to simulate the ROM, but all nonlinear terms in the twenty—DOF original model must be computed in order to simulate the original model. 164 4.6 Conclusions In this chapter we have shown how to extend the fixed-interface linear CMS technique of Craig and Bampton ([62]) to nonlinear structures by making use of fixed-interface NNMS in place of fixed-interface LN Ms. This approach allows one to build nonlinear reduced-order models with improved accuracy for systems that are composed of assemblies of substruc- tures. The first example system presented here is quite simple, and the model reduction is not very significant (from five to three DOFs). A second system with many more DOFs is also used, and it shows a significant reduction in model size (from forty-one to three DOFs). In fact, the N NM method has been developed to the point where it can be applied to systems with thousands of DOFs. The roadblock to pushing this method further is that it is very computationally demanding to generate substructure ROM with more than one DOF using a NN M approach ([42]). Current work is focusing on systems with many more DOFs at the unreduced substructure level, and it includes substructures that are modeled by nonlinear FE techniques. 165 4.7 Tables 4th Mode of 1th Mode of 1 0.68596232102650 11 44.72697600017000 2 5.06072722592600 12 48.09670186700000 3 9.94306090168100 13 51.17011958090000 4 14.79648130027000 14 53.92824229550000 5 19.56717359004000 15 56.35403805246000 6 24.22073003072000 16 58.43253131062000 7 28.72676467269000 17 60.15089310481000 8 33.05675986752000 18 61.49851874983000 9 37.18364689122000 19 62.46709231770000 10 41.08177243723000 20 63.05063732260000 Table 4.1. F ixed-interface linear natural frequencies of substructure B using the second—trial parameters. 7th Mode 0217 7th Mode w? 1 1.41414755239200 11 92.20702318159000 2 10.43295644747000 12 99.15389102074001 3 20.49814516979000 13 105.4899039538000 4 30.50372764424000 14 111.1759195941000 5 40.33876175326000 15 116. 1768256602000 6 49.93231412311000 16 120.4617492830000 7 59.22174251389000 17 124.0042428732000 8 6814825628723000 18 126.7824443123000 9 76.65605183909000 19 128.7792098755000 10 84.69224352313999 20 129.9822187218000 Table 4.2. Fixed-interface linear natural frequencies of substructure 7. 166 4th Mode 47,1401” 5 4th Mode waCMS 1 0.68582383924810 1 0.68582383933910 2 1.41280504607300 2 1.41280505021300 3 4.00456599773800 3 4.00460940350200 Table 4.3. Comparison of the linear natural frequencies of the three-DOF linear-CMS nonlinear model and the three—DOF nonlinear-CMS model. 7th Mode ijCMSAlDOF 1th Mode 018M 1 0.685823839247 1 0.68582383924720 2 1.41280504587300 2 1.41280504587300 3 4.00452879688600 3 4.0045287 9688500 4 5.06079466571000 4 5.06079466571000 .5 9.94306486412500 5 9.94306486412500 6 10.43301741537000 6 10.43301741538000 7 14. 79648235748000 7 14.79648235748000 8 19.56717401354000 8 19.56717401354000 9 20.49815228625000 9 20.49815228625000 10 24.22073023840000 10 24.2207 3023840000 11 28.72676478757000 11 28.72676478757000 12 30.50372969952000 12 30.50372969952000 13 33.05675993624000 13 33.05675993624000 14 37.18364693456000 14 37.18364693456000 15 40.33876259811000 15 40.33876259810000 16 41 .08177246556000 16 41 .08177246556000 Table 4.4. Comparison of the first-sixteen linear-natural frequencies of the forty-one-DOF linear—CIVIS nonlinear model and the forty-one—DOF original model. 167 4.8 Figures é \‘ \\\\\\\\\\\ g g a Figure 4.1. A five-DOF nonlinear spring-mass system. [i substructure X‘l r—> 2 an .// .4va m. , '7, ———> F: ///, 2,2 l ./ —’ Y x1 x5 X3 or substructure y substructure Figure 4.2. System substructures. 168 . . . ...”...tan. . 4 4 . . .... . . e 44 x 33.: : 3.... \\ $¢s¢~ -~-~-~-~ J \ é $5 .. i .. ... see 4 . e . ... \ ss¢¢¢¢¢ ... ‘\\\ S 5‘ ~. Q ~ s 4. 4 ~ ~¢ \ \ x“ \“ \ \“0‘ s ‘s x ‘\ .-¢“s . \x . ‘s ‘ .r‘s “ 0‘ ‘s ‘\ -- s . ¢ .8 ~ S... : X‘se‘ekscxtz. ... . . \s‘ss‘s s~ - - .3 - .o u o e $3333... ,.... . .. . 4.4.4... §ss§§s QQQQ .\ S “ Q‘ §Q ‘5 Q‘ Q \\ s § 4 ~§ Q \éfite ‘ ‘3 \ \‘a N fl itude 772 ' nonl h first fixed-interface Xfi( N’fi,7'/[V‘fi) tot e 2 I Figure 4 \‘A‘ 5 \ s .63 $66 . “ 8.84.44“... ... ~-N-WM~§SM§~WO h. - QQ §§ -~ :h -~ - 1.3"...4.....~4~.N.Nu¢ . ~ ~ ._ in g--~--~~0~§~§~fifi~ . ......w... - 44.............:...: 45:”... ~ ~ ~ é ex: ¢ $sss-s-n~nwwu ~ -~ ssssss§~s9-- ss ‘ ss ss - ¢ xxx 9 ss s 3‘33 s s o‘§‘\‘\“‘\§§$§§ .§‘&“‘~§§~ Q‘.“‘§“‘§‘ ‘ §$ “ o. “ \ 349‘s s .... v .... ...: “ 0 ~.' #ué ss s “$44 . é ~ Q Q. ~Q § ‘ ..N.....:$- ~ . 4 ~ sswsswos #3» .00. ss \\ ss§ ‘5‘““\‘ ....s § sss ......» QMCx 6 ss ~ s §§§W§§S .. ~u~ . ins : s . ...... - Na : m e 7. bstructur 'fold of su ' de am face nonlinear mo d-inter first fixe 2 1 ’ 169 0.25 I I r 1 1 I L - - NLCMS. 3 DOF (12. --LCNEL3[M3F a LCMSLPrb. 3 DOF Original Model. 5 DOF 0.151 ‘ H 1‘ I‘ ’ 1 o 1 r " '1 .1 '1 ‘ '1 « o ‘ . \ I 1 ‘ I /‘ I 1;“ " 1" I1 I ‘1‘,“ I 1 If. 1 1 I” ‘ '1‘” ,6" 1.1 ‘1 ‘I‘."‘ I [,11 0‘05 b' i {I ' ‘ 1| "1“ | ,.>‘ f \1\’ f .l 1 1 ' V '11 '1 ' '11“ 11.1 “ I)1” " 1| I . lr"1\’|\ I, ”h If.“ [1‘ I‘1!"‘1 ’ f :\1=.1I, [ ’1’ I“ ‘ t l’ 3 I" "' 7 . ,‘ .‘ ‘ I . xv— OL \' I, 1‘ “| ;’.‘\ I I\ [ll 1 ’1‘ '1 I “ ‘ '1' \II “(7 i“ V" ,fl I l H I J I I 1‘ .- N l "‘l " “‘11 ' I" ' ‘ ' ‘ f1 -‘ 4‘ 1 ‘I 1 I" il I ._\‘ I .1 ‘1‘]11 1' ‘1‘12'1 1‘), 1' {1: ‘1‘} [i '11,,” ‘ \11” ,1f -0‘05 ’- ‘ ' l I ' l ‘ | 1 J I." ‘ ) [I J . I 1‘1!» 1‘ 1 I 1 1’ ’7' \I 1 . 1 \ 1 ‘ I 11 .1 ‘ \ 1 I l , I u _ Pl‘ ‘ II | 1’ x" 11‘; ‘ '1 " 1' -001 [- 1! ‘1 I" " 1' \l_, l) f] 1’ [I I —O.15 h ...1 —O.2 - a _O.25 EL 1 1 1 1 1 1 O 5 1O 15 2O 25 3O 35 40 Figure 4.5. The time response of the displacement $1 of mass m1 from Figure 4.1, which corresponds to x?, :13? , and x] from Figure 4.2, for the nonlinear-CMS model, the linear- CMS nonlinear model, the linear-CMS linear model, and the original model. 2 r - - NLCMS. 3 DOF ~ - - LCMS. 3 DOF 1 5, LCMSLPrb, 3 DOF _ ' Original Model. 5 DOF 1‘ 1:}. I,“ 1,.‘\’\ ' ”“1 A" 17‘1 'A d 1 :I 1‘ ’-’ 11 I‘ ‘ I ’l \ I 1 ‘ , .., \ 1 II 1‘, ' 1 (I ‘ \ 1 ' ‘ 1 ’ I ‘ 1‘ ,‘ 05"I 1' I, I 1! 1‘ 'I1\ I '1‘ ‘ I 11‘ "‘ O ‘s H‘ [I ‘ ‘ , 1 \ ' I ‘ 1 I , I ‘ I ‘1 ‘ 1;. ‘ I " I, 1 ‘ f i ‘ ‘ ' ’I ‘ ‘ I I ‘ “ ’ ," ‘ 1| N I ’ I ‘ 1 ‘ 1 x O ’- ‘i I ‘ [I ‘1 \ I 1' ‘ l ’ I ‘ ‘1 ,I ,l ‘ 1 ’ ; ‘ 1'1 I 1 I1 ‘ ‘ ‘ , 1 \ 1‘ ‘ .‘ I . . , 1 I ‘ . f . I I, (I I," 1‘ I ‘ 1 1‘ ' 1 ‘ ‘ I ' \ \ I i \ ‘ -O-5 _ \ 1 \\ I; \ \ I ‘ \ 1‘ ’ I \ ‘I \ f \ I I‘. 1/ ‘1 1 I ‘ ‘ I I ‘ "I l ‘ ' ‘ \ d I ' I/ ‘ \f I \ ’ / ‘ H ’I \ _1__ ‘J ‘\_g {\1 \,\/ ,\ \/ \ [L —1.5 7 a _2 >- _. 1 1 1 1 4 L 1 O 5 1O 15 20 25 3O 35 40 t Figure 4.6. The time responses of the displacement 3:2 from Figure 4.1, which corresponds to 32/23 from Figure 4.2, for the nonlinear-CMS model, the linear-CMS nonlinear model, the linear-CMS linear model, and the original model. 170 2 *- fi r Y I' L If L -5 - - NLCMS. 3 DOF - - LCMS. 3 DOF 2 * LCMSLPrb. 3 DOF ~ 1 5 Original Model. 5 DOF . f" r ,\ 1‘ II/\ ‘ ’6" ll f" \ : ‘ ‘2 ‘ 1" If "1 /’\\ 'I ‘ I ‘1 l I V \ I \7 f ‘ . , H I \ \ \ I y/\ . A I 7 i ‘ I! I‘ ‘ ’ I 1 1 l 1 \ ’ I 1 I 3 '1 05*) H u ,, “ "f i I 1‘. "\ ‘\ l :1 I... \ 1, ‘ 1' H I . l l I 1 ‘ | , i l .\ ‘ "i m [l N ‘) '5 ‘i 'f “ " 1] 1' ‘ l ' ’1 ‘ \l '.\ x O ' ,1 1, .t ‘1‘ , I 1 1‘ I '1 ‘ 1‘ 1 ’1 1 1 I I ‘ 1‘— ' 1 '. I - [ 1 q! H. II 1‘. ,’ I ‘ l , I I \ ’ 1' \ i , .5 1 I —0.5 " \ [I \\ I ' \ ‘ I ‘ ‘ I ' \ ‘ (q ‘ 1 1 \i ,l \ i 1 ‘ I ’ ‘ I ' ‘ I, I ‘ 1 L f 1 \l I \l f ' ‘ U ’ \ V I \ } l 4 _ \ '1 \ I \\/ ’ \ ) I \ /\ \/§ I \ VI \ ’ ., \/ .\’ \Jl -’ -1.5 ~ ~ _2 - -4 —2'5 >- 1 4 1 41* P J 2_L ‘ O 5 1O 15 20 25 30 35 40 t Figure 4.7. The time responses of the displacement :173 from Figure 4.1, which corresponds 5 to $3 from Figure 4.2, for the nonlinear-CMS model, the linear-CMS nonlinear model, the linear-CMS linear model, and the original model. - ~ NLCMS. 3 DOF 06* — — LCMS.3DOF _ ' . LCMSLPrb. 3 DOF Original Model, 5 DOF 0.4 F .1 q I} I \ I“ I . . I ‘ 1, 11 I ‘ 'I 1 I \ ‘l \. 0.2 ~. ,’ .\ v I ‘. 1r 1‘ I \ 1 l ‘ /'\_ .1 ‘1 ‘ 1 1! \1 l.\\ ’ \ g ‘ f l I’ \ j \l ' ‘ i 11 i I ‘ j / ‘ 1 1‘. Vt ' ‘ II \- A 1! ‘1 ' ill If f 1 ’f : A“ X 0’" I f ‘1 I, \\ / l I; i- ‘ 1‘ ’1 1 ’7 K I] ‘1 -[ . ‘. I I , 1 1 1 1 1 I ’ (ll ‘\ I 1; ‘1 1 \. I ‘ ‘I ‘ If H a f 1 '/ 1 . ,t ‘1 ,' 1‘ I «(I ' I, 11 -—O 2 k- : , s/ \./ I ‘l 1' \ I \\ . ‘\ / I 1‘ .4 l '1 I 11 1 \-l I ‘ 1’ 1“ L- f l ,J' 1 ’ \l ‘J \‘l 1 I \ —O.4 ' a -0.6 ~ - O 5 1O 15 20 25 30 35 40 t Figure 4.8. The time responses of the (:lisplacen‘ient .174 from Figure 4.1, which corresponds to 1'3 from Figure 4.2, for the nonlinear-CMS model, the linear—CMS nonlinear model, the linear—CMS linear model, and the original model. 171 — - NLCMS. 3 0015 0.6 t -- - LCMS. 3 DOF F‘ LCMSLPrb, 3 DOF Original Model. 5 DOF 0.4 r . . .1 I" - 9, I“. I 1 ‘7 " \1" II I" I] ‘1 ‘1 i '1‘ 5' IA‘ : {Tl 0.2 r I A j .\ 1, 1, 1/1‘ _, ' ;~ \1 5 '1- r . , I 11 1 i \l I \\ ' ' .‘ ‘, I l ‘ \‘ . )l 1 i / 3 I \l l I II \ I { R. I ‘I '1’ ‘ 1 . I 1 to ‘ , I, , I , 1‘, l 1“ I l , II ‘ 1“ x 0 ” ‘1 f ,, ,1 I ,, 1; I 1‘ ' I i I! ,I 11 ‘ I I I . 1. 1 \ I \ 1 \. I l . ‘ I \|\ ., l’ ll / .1 / I \ \ I } ’I ii. '1 ,I 1 , 1x, \ ,2 1 , I 1 1 1 «1 - 1 11 11 _O 2 1- ‘ 1 ‘1 \J’ .’ i, ’ ll 4 \ I ‘\ l I \l -4 ' -. I \ II. “‘ ’ . \ I V l I ‘1‘“ '1 ,1 1 ’ "J ” " 1 I w J ‘1‘) ‘i,’ ‘u‘ ‘ A —o.4 - , « —O.6 - J O 5 1 O 15 20 25 3O 35 40 t Figure 4.9. The time responses of the displacement 9:5 from Figure 4.1, which corresponds to :rg from Figure 4.2, for the nonlinear-CMS model, the linear-CMS nonlinear model, the linear—CMS linear model, and the original model. XNmass meass+1 I“? -~+ 'kz Nmass’lfz Nmass +1 X 2Nmass x2Nmass+1 ,—> i? k 4 Nmass’ , 4 Nmass + 1 xNa Nmass XNa Nmass + 1 1’ 1 b l t k IkzNaNmassI 2NaNmass+1 WIND Figure 4.10. A class of nonlinear spring-mass system. 172 +22 x9x2x9xfi x2 4 mass xamassfl [—> kg Nmass’k2 Nmass ¢ 1 If} h. :ke. ‘ m kiwi 1 u 7 ' F043 F30: I m m a £7 7 W Bsubstructure J J / , /, » / Xlrtmass meassH + r” ngmass'k; Nmasu1 , / / m? t + F3” L7 L4 7 substructure // / p i D b j" + v i i 1 ‘ ‘ j 1 ,I / __’ ‘ A ‘ /7l 7 z / / a substructure Figure 4.11. System substructures. —~—— Fixed-Interfage Linear Mode Shape 1 0-5 “ — — — Fixed—Interface Linear Mode Shape 2 fl - Fixed—Interface Linear Mode Shape 3 0 4 _ — Fixed—Interface Linear Mode Shape 4 0.2L s ax... O i‘ J _0.2 _ - -O.4 ~ ‘ —0.6 ~ - O 5 10 15 20 25 im Degree of freedom Figure 4.12. The first four fixed-interface LNMs of substructure ,3 using the first-trial parameters. 173 A 1 500 n 1 000 60 40* -40- 500 —60 m term n Figure 4.13. The coefficients of the pqr cubic nonlinear terms expressed in fixed-interface linear modal coordinates and associated with the first mode (master mode) of substructure 3 using the first-trial parameters. s s é é sé ss \ssxss ésss‘ ~ é s ss 0. ss é ss s 0 \s\ x s s~¢s~ss s ~§~s~s§s§ \\\§ us 3‘ s~ x x. ¢ 5 s. s~ s ~s - é ys L é ~ ss ~ s s fé -- s ‘ \~ ~ - #Néhé sssHMsHVQ ss‘sVQ s~ ¢s-~ $‘ é o : é ~ s~ ‘5 A.“~%--h~§§ §§~¢§§ “~§‘ 33:: *av ~ss ss ost§bp . 3..: ~ s ~ 9? ., ......:..&$ sewex ss s ...~ ~ ~ s s V I ~ ~ ~ ~ § . . ...w---“-“s~“~\s .....u.....:..:u¢$¢ws ‘ . 2????4is ‘1» 1111111~és s11. .1111..111 ‘1 . .1115: .....u::. . 0 .~ . ~ -~V§ ~ s ~ essfiss ss \ x \ \ \So SS s 1 I ¢ ¢ s ~. -~ ‘ . s ¢ ¢ ~ ~ : . s s®¢~¢~¢$$¢ ...: s @033 ¢¢~$¢¢~$- $3 s 3 ¢ ¢ ~ .. .. . s \ ~§ ss s~ Q ~Q ‘0‘ ‘§ Q s a: eoo»e s» s o8 é 3. ~u¢ ~ ... o $ $x~$.:.....i.. , ¢ cox ~.~.-~ .~ ~ ~ s 3% ~ - o I \ s 9 ss s s ~¢ ¢ ¢ s s 11 s ‘ s s ‘ ssss~ ~ ~ 0 reevewgvfi... ~ s ss es Q QQQNQQ s Q o s ss‘ssossQ‘sé‘s‘QQQ s s e ~ s s s s s s s. s s so s O s ésés essavsss ssss NMB Figure 4.14. The contribution of the second fixed-interface linear mode amplitude 712 X 5 ( {V46 , niv‘fi) to the first fixed-interface nonlinear mode manifold of substructure B using the first-trial parameters. 174 Figure 4.15. The contribution of the third fixed-interface linear mode amplitude "91,5 = X g (”{Vfl , fiiv‘fi ) to the first fixed-interface nonlinear mode manifold of substructure fl using the first—trial parameters. Figure 4.16. The contribution of the fourth fixed—interface linear mode amplitude rig/fl = X :3 (niv fl , fiiv’fl) to the first fixed-interface nonlinear mode manifold of substructure B using the first-trial parameters. 175 50 rLtlLllllldl ' a 4 a Inn-UH.” lllll NNM. 1 DOF —-—- Orlglnal Model. 20 DOF A 50 40 30 20 1O Figure 4.17. The time responses of the first fixed-interface mode displacement "{Vfi for the ROM and the original model of substructure {3 using the first-trial parameters. iflwfl) for the ROM and the original model of substructure fl using the 176 The time responses of the second fixed-interface linear mode displacement H N J3 first-trial parameters. ure 4.18. 5.12 1 772 Fi 0.5- 1 «a 111 111 I 111111 1 11111111 1111111 1'1'11111 °‘* 11111 11'111111111111111"" 1111111. 0 1O 20 t 50 4O 50 Figure 4.19 The time responses of the third fixed-interface linear mode displacement n3 N1fl=3X (nN’ 5,119”) for the ROM and the original model of substructure fl using the first- trial parameters. " ' ' — - N.M 1oor= —— OI'l lnal Model, 20 DOF 1.5“ " 1 0.5‘ 1 1 11 1' , 1 a; O " ‘51.1'1’11"1‘1""1 : 1 1-1'1'1 ' '1‘1'1111'1 1 :1 "',""'11,"1 ' 111‘ S—‘v ‘1"1'1'1' "1 ‘ ' ' ""1" "1" " ' '11'1111"""'1 "‘1' —o.5- 11 1 —1.5r -2 - 0 1O 20 30 40 50 Figure 4. 20. The time responses of the fourth fixed-interface linear mode displacement N—fl X f ( N fl, N’fl ) for the ROM and the original model of substructure 3 using the first- trial parameters. 177 60 , , r LNMiPrb. 1 DOF - — LNM. 1 DOF —— NNM. 1 DOF 0 1 O 20 30 40 50 Figure 4.21. The time responses of the first fixed-interface mode displacement 721V ’16 for the linearized one-mode model, the one-mode model, and the ROM of substructure 5 using the first-trial parameters. 0.6 —~— Fixed-Interface Linear Mode Shape 1 ~ — -— Fixed-Interface Linear Mode Shape 2 - 1 Fixed-Interface Linear Mode Shape 3 0.4 _ ~ Fixed—Interface Linear Mode Shape 4 H l I I 1 f .— r0 0.2 '- / / f I >\:’ .1 tax" 0 » , » ’ '1 ~ ' / ’ / I. \ \ \ .l /' ..0.2 ~ ‘\ j , " A / 4 1% 1 1 _. ’ T 1 \ —O.4 ~ - —-O.6 r 1 J 1 . ‘ O 5 10 15 20 25 im Degree of freedom Figure 4.22. The first four fixed-interface LNMs of substructure 3 using the second-trial parameters. 178 600 800 1 000 1 200 1 400 1 600 nm term 400 4G —4c 0 Figure 4.23. The coefficients of the pqr cubic nonlinear terms expressed in fixed-interface linear modal coordinates and associated with the first mode (master mode) of substructures 5 using the second-trial parameters. Nfl = Figure 4.24. The contribution of the second fixed-interface linear mode amplitude n2 X g (”{Vfl’ nil/fl) to the first fixed-interface nonlinear mode manifold of substructure B using the second-trial parameters. 179 Figure 4.25. The contribution of the third fixed-interface linear mode amplitude név’fl = X g (niv ”[3 , 7'7?” ) to the first fixed-interface nonlinear mode manifold of substructure B using the second-trial parameters. Figure 4.26. The contribution of the fourth fixed-interface linear mode amplitude nil/fl = X fluff} , nil/'5 ) to the first fixed-interface nonlinear mode manifold of substructure fl using the second-trial parameters. 180 2.5 l w 1 - - LNMLPrb. 1 DOF 2 _ LNM. 1 DOF - - NNM. 1 DOF 1 5 _ Original Model, 20 DOF 1:3. \ \ \\ r I // 0.5F \ \'\ / / .1 c: \ \ \ / I // Z 1— 0’. \ \\ / / 'l \ / \ I / \ I /\ / “0.5 h \ \ / / / T 1‘ I. \. l -1 i- \ l/ \ \ I I I \\ l/ " —1.5~ « -2 _ q _2.5 1 4 A 4 0 2 4 6 8 10 t Figure 4.27. The time responses of the first fixed-interface mode displacement nil/fl for the linearized one-mode model, the one-mode model, the ROM, and the original model of substructure ,3 using the second-trial parameters. V Y I — — NNM. 1 06F 0,15 _ *—* Original Model, 20 DOF H I l —O.15 )- p— .- .— Figure 4.28. The time responses of the second fixed-interface linear mode displacement 77‘?)V ”3 = Xgmiv’fl, 1'79“} ) on the first fixed-interface nonlinear mode manifold for the ROM and the original model of substructure ,8 using the second-trial parameters. 181 YY)—T .0 O .5 l .0 o w L — — NNM. 1 DOF —— Original Model. 20 DOF Figure 4.29. The time responses of the third fixed-interface linear mode displacement név’fi = X :43 (nfV ’fi , fiiv’fi ) on the first fixed-interface nonlinear mode manifold for the ROM and the original model of substructure fl using the second-trial parameters. 10 0.015 5— NNM. 1 05F ——- Orkflml Model. 20 DOF 0.01 M\\ I f A \\ II I, 0.005l \\ _ n \ 2“ v o~ \ ‘ p _ : —o.oosi ‘ \ ‘\ —o.01 L ; ‘7’ ‘ —o.o15 ‘ ‘ l o 2 4 6 8 Figure 4.30. The time responses of the fourth fixed-interface linear mode displacement ) on the first fixed—interface nonlinear mode manifold for the ROM Nfi .Nfl N, ' 774 H=Xfil(nl .721 and the original model of substructure [3 using the second—trial parameters. 182 10 0.5 T T L L '— ‘ - LNMLPrb. 1 DOF 0.4L- LNM, 1 DOF H - - NNM, 1 DOF 0 3 - Original Model. 20 DOF 0.2 " l/I, w ~- \\I I I .\>\ ‘ Ill/,... “\\ . -1 0.1 L (I / \\ \'\ I \'\. / " 1‘ /\ \. -l / / I r" // '\ \. CORN O T j/ / _ \y / '\’ \\ 1 I I \. / \ \ \ -O.1 r /,I l.’ \y [,1 \ \\ _ If I \I\ II" \\ ‘ \.\ —o.2~r ,’ “" . "t : —0.3[- - ‘0'4i a —O.5 1 L 1 L O 2 4 6 8 10 t 5 Figure 4.31. The time responses of the displacement $2 from Figure 4.11 with the interface fixed for the linearized onemode model, the one-mode model, the ROM, and the original model of substructure fl using the second-trial parameters. 0.6% ' —-- LNMLPrb, 1 DOF ~ - LNM, 1 DOF - — NNM. 1 DOF 0,4 r- Original Model, 20 DOF H II/ \ 3 I , _ § ‘ \ I, / \ \ 02* // ,‘1/ \\ x/ \\ J I” 1.] \ \I/ \l\ v— ,/ I \ f \ \ \‘ QXN o l" / ’ / \\ \ \\ _, / I \ / \ \ ,/ x/ \ // \ \ _02P I, / I / \\ ’1’ \\‘\ \\ : J4 I \~\> ‘ I ‘ — ‘1 —0.4 " .4 —O.6 r L l l l a O 2 4 6 8 10 t Figure 4.32. The time responses of the displacement xgl from Figure 4.11 with the interface fixed for the linearized one-mode model, the one-mode model, the ROM, and the original model of substructure ,3 using the second-trial parameters. 183 O ‘, ' ; ' —e— Fixed-Interface Linear Mode Shape at t = 0 _+_ Fixed-Interface Nonlinear Mode Shape at t: 0 —O 05 _ _O Fixed—Interface Linear Mode Shape at t = 0.125 Tm _ - , Fixed—Interface Nonlinear Mode Shape at t = 0.125 T"l O Fixed-Interface Linear Mode Shape at t = 0.25 Tnl —0.1 r ( Fixed-Interface Noninear Mode Shape at t = 0.25 Tn. J —0.1 5 - _, g . cn.x._ 0 Q4}? 956 0‘9 O—GQ-O 0—9- 0-0- O-GO —o.2~ ‘ " l \K\ —o.25 ~ *\~\ J \‘K*.\ —0.3 " *‘fi‘w‘ *Hfi ' -—0.35 ‘ 4 t ‘ 0 5 10 15 20 25 ith Degree of freedom Figure 4.33. Comparison of deflections :L‘fU) for a quarter-period of motion on the first fixed-interface linear mode manifold (linear eigenspace) and on the first fixed-interface nonlinear mode manifold of substructure B using the second—trial parameters. The motion starts at the maximum deflection (the bottom curve) and moves as shown to the zero 5 deflection at a quarter-period. Note that linear mode shape is normalized such that $1 is ’3 . . the same as ar’l of the nonlmear mode shape at each time. 0.6 r —-— Fixed-Interface Linear Mode Shape 1 ~ — - Fixed—Interface Linear Mode Shape 2 Fixed—Interface Linear Mode Shape 3 0.4} Fixed-Interface Linear Mode Shape 4 . / a g , , - ’ '1'“ 0.2 - , , ’ _ ~ / / /' \ >‘x—O __ i I / . \ J I / 7" \ —o.2 - ‘, , - ’ ” - _ j J -O.4 r ~ -0.6 l- 1 l l l d O 20 25 , 10 1 5 Im Degree of freedom Figure 4.34. The first four fixed-interface LNMs of substructure 7. 184 600 800 1000 1200 1400 1600 nlh term 00 200 200 150’ -150- —200 Figure 4.35. The coeflicients of the pqr cubic nonlinear terms expressed in fixed-interface linear modal coordinates and associated with the first mode (master mode) of substructure Na = Figure 4.36. The contribution of the second fixed-interface linear mode amplitude 172 '7) to the first fixed-interface nonlinear mode manifold of substructure 7. .N 37]] 185 Figure 4.37. The contribution of the third fixed-interface linear mode amplitude 17?,” = X ; (niv’l, 779]”) to the first fixed-interface nonlinear mode manifold of substructure 7. Figure 4.38. The contribution of the fourth fixed-interface linear mode amplitude nil/'7 = XZ(n{V’l,17{V‘7) to the first fixed—interface nonlinear mode manifold of substructure 7. 186 — - LNMLPrb, 1 DOF 2" LNM. 1 DOF ~ - - NNM, 1 DOF 1.5 P Original Model, 20 DOF \\ I \ 1" ~ / N \ ‘ , I I T \ ‘l 1 ” \ \ ‘/ \\_ ’ C \ I ./ \ \ f / 0.5- ‘ \ \ I / J \ / I )1. \ \ _’ / z r- — \ \ I I _ c. 0 \ z \\ , I I/ \ / /\ / —0.5— X ’ \ / ‘ I \~ L /" \ / I \\ // -1 \ // \ ‘ _ I / \\ x z/ T -1.5* 1 -2» 7 O 1 2 3 4 5 t Figure 4.39. The time responses of the first fixed-interface mode displacement "{Vn for the linearized one-mode model, the one-mode model, the ROM, and the original model of substructure 7. 7 I - — NLNM, 1 DOF —- Original Model, 20 DOF H 0.1 0.05 >- p— 1.. 1- Figure 4.40. The time responses of the second fixed-interface linear mode displacement ; 7 n3}7 2 X; (17fV ’7, 17?“) on the first fixed-interface nonlinear mode manifold for the ROM and the original model of substructure 7. 187 0.04 1 fi - — NLNM. 1 DOF ~—- Original Model. 20 DOF 0.03 t _ "\ 0.02r \‘ 0.01 i >: Z or) Ol- 1;: -0.01 ~ —o.02~ \A ' —0.03 r a —0.04 4 L ‘ g 0 1 2 3 4 5 Figure 4.41. The time responses of the third fixed—interface linear mode displacement I)?” 2 X37 (niN’l, nil/‘7) on the first fixed-interface nonlinear mode manifold for the ROM and the original model of substructure 7. I - — NLNM, 1 D]OF 0.015 r — Original Model, 20 DOF r\ C: 3—3 \‘\ I’ I, 0.005 r ‘\ \ 1 N"! 4 O ..,_,u "s c \ / —0.005[ \\ \ n -o.o1 l \..// .. / . —0.015r 7 Figure 4.42. The time responses of the fourth fixed-interface linear mode displacement ; r . . 1):,“7 = X](niv’7, 7'19”) on the first fixed-interface nonlinear mode mamfold for the ROM and the original model of substructure 7. 188 . Y -—I — LNMLPrb, 1 DOF 0_4L ‘ LNM. 1 DOF a -‘- — NNM, 1 DOF Original Model, 20 DOF 0.3 r 0.2 r , , , A \ a I \\/ ,./ ‘\ I .\ \ / 0.1 L \/ \»\ ‘1 l ,, ‘x \ >$‘>< 0 ’- 1 / / I \ / \ \ \'\ d / / \- / .\ l\ ’ I \\ \ \ / I l [I \ \ \\ —O.2 ’- "I I f \\ 1’ \ \ '\ z 4’ I \ J ‘ _ ,.\’ ./ —O.4t - —O.6 l“ 4 J 1 _ 0 1 2 3 4 5 Figure 4.44. The time responses of the displacement $31 from Figure 4.11 with the interface fixed for the linearized one-mode model, the one-mode model, the ROM, and the original model of substructure 7. 189 O \ ' l —-e— Fixed-Interface Linear Mode Shape at t = O l -+ Fixed-Interface Nonlinear Mode Shape at t: O —0.05 _ l O Fixed—Interface Linear Mode Shape at t = 0.125 Tn. , l , Fixed-Interface Nonlinear Mode Shape at t = 0.125 Tnl \ O Fixed—Interface Linear Mode Shape at t = 0.25 Tm —O.1 r \ Fixed-Interface Nonlinear Mode Shape at t = 0.25 Tnl r —O.1 5 ~ \3 _Q a fig- l @ ? TH? 0‘0“0‘0-G“O»e—o—e~o—e-o~e—o .412? - —O.25 r - -413- * -O.35 ‘ ‘ ‘ ‘ O 5 10 15 20 25 ith Degree of freedom Figure 4.45. Comparison of deflections $170) for a quarter-period of motion on the first fixed-interface linear mode manifold (linear eigenspace) and on the first fixed-interface nonlinear mode manifold of substructure '7. The motion starts at the maximum deflection (the bottom curve) and moves as shown to the zero deflection at a quarter-period. 2.5 7 I I I '- - LCMSLbe, 3 DOF 2 " ' LCMS. 3 DOF .., ' ’ NLCMS. 3 DOF \ ’ f ‘ fl IA - \. 1’; 1‘, " “ll [\ Fu‘ 1 ”,4, , 1" i’\ 1,", - \\ ‘4 I ;I “ll ‘ 1‘ ill! ‘. I’ll .: N l’ \ , ‘il 0.5 _ l 1‘ l \ l , )I \f l‘ l 1» ‘ y l‘ II I l.‘ , i “I , f? l. .i l l ’I ll “ ‘ l I. \ l I i , $ ’ . n. \ l l l ’ ‘\ I l‘ l [’4 \ ‘ 1 ‘l f ’ . I H I, | x 2 v- 0 ~ ' ,I l |\ z , l l x, l I ’ \\ l '7 \ .-4 :- \ I ’l l ‘l I i ‘ ‘ ‘ ‘l l" , ll I, ill I \ \ II \ I ‘ I l l l “ l g _ . ‘ . ‘ \ _0 5 ’- l. ( l x l \ l l' l I l ‘l I;‘ l A b I l‘ (r! a l K ll, '|-.i '. (l \l l l A} \ ‘1' ‘ / '- ‘z‘ ’ \ l A} j ‘ (I ‘ r l 'I ’1 ‘ f,- I ‘1’ \ _’l' I -‘ I \ I l‘ , “j" ’ V \ —1~ . , , , ,Jg, ”1,, a ‘, , \j i, p v —1 .5 r . —2 d —2 5 4 . . x 0 10 20 30 40 50 t Figure 4.46. The time responses of the first synthesized fixed-interface mode displacement 179/fl for the three-DOF linear-CMS linear model, the three-DOF linear-CMS nonlinear model, the three-DOF nonlinear-CMS model, and the forty-one-DOF linear-CMS nonlinear model of substructure B. 190 0.05 — - NLEMS. 3 DOF —— LCMS, 41 DOF 0 10 Figure 4.47. The time responses of the second synthesized fixed-interface linear mode displacement név‘fi for the three-DOF nonlinear—CMS model and the forty-one—DOF linear- CMS nonlinear model of substructure fi. 30 4O 50 0.04 1 - - NLCMS. 3 DOF —— LCMS. 41 DOF 0.03 h '- 0.02 ' - —0.02 —-0.03 —0.04 i i . l ll (3:030- ll -o.o1{ W K Mr} : l 1 0 1O 20 3O 4O 50 Figure 4.48. The time responses of the third synthesized fixed-interface linear mode dis- placement név’fi for the three-DOF nonlinear-CMS model and the forty-one-DOF linear- CMS nonlinear model of substructure 5. 191 0.015 . r , - - NLéMS. 3 DOF —— LCMS. 41 DOF 0.01 n [l .’ ~ 0.005 [\f NB _1 -o oos~ l '[ -ool V V .. - —o.01 5 . . g % Figure 4.49. The time responses of the fourth synthesized fixed-interface linear mode displacement niv’fi for the three-DOF nonlinear-CMS model and the forty-one—DOF linear- CMS nonlinear model of substructure i3- — - LCMSLPf‘b. 3 DOF 2 * LCMS. 3 DOF t '— NLCMS. 3 DOF 1-5t LCMS. 41 DOF g l‘ \ / f 1»\ ’1‘ i, ,x ’l ’8 I,“ A Mn '. I y l ‘ f ' l- \ I ’ \\ I \ il i ll | / I x‘ [ , l I‘ ‘ ‘ I \ . I ‘ \l ‘I i 'I“ \ ‘ l \1\ j ‘ \ J 0-55‘ ’ I" ,4 i g l ‘ l l I! ‘ ‘ ' l ' ' (l l . \ ’ ‘“ I I H " 1" l ‘I’' ‘ U \‘ I l ' ’H I >: \ l l . '- l ‘l \ , 'l~ \. ' \l ( :- \ l \ i \ l I, l l‘ l 'p \ l l y ’ 4 ; I l | , l ' l l l l l ‘ r , -0.5 '— \ I \\ I l \ l l, l l \ 4' l if I‘ll I ‘ l., ‘ I l \ "l l H I i I ‘ I \ \ I \ \ l I y ‘ l ‘ . L ‘ I ( I ‘ ‘ \ j \ \/ \_ I \l I \VI ‘1‘, —1 \’ ll 4 “I U ‘W \’ {LN -1.5F fl _2>- _‘ 0 5 10 15 20 25 t Fi ure 4.50. The time responses of the first synthesized fixed-interface mode displacement r); ’7 for the three-DOF linear-CMS linear model, the three—DOF linear-CMS nonlinear model, the three-DOF nonlinear—CMS model, and the forty-one-DOF linear-CMS nonlinear model of substructure '7'. 192 0.2 r t — - NLCIMS.3DOF ‘ —-—- LCMS. 41 DOF Figure 4.51. The time responses of the second synthesized fixed-interface linear mode displacement 779/” for the three-DOF nonlinear-CMS model and the forty-one-DOF linear- CMS nonlinear model of substructure 7. 0.04 — - NLCrMS. a DOF w— LCMS. 41 DOF All '. 0.03 0.02 0.01 N91 —0.01 f _.__I___...-—+—r r M‘ ‘2. ~52: v} 1" ‘7' -o.02L ‘ -0.03 - —0.04 Figure 4.52. The time responses of the third synthesized fixed-interface linear mode dis- placement 179’” for the three-DOF nonlinear-CMS model and the forty-one-DOF linear- CMS nonlinear model of substructure 7'. 193 0.015 . . 1 - - NLéMs. 3 DOF —— LCMS. 41 DOF 0.01 l ' 0.005 ‘ ‘ My 114 o Kr... "/ -o.oos ~ l -o.o1 - l4 ~' “0'0150 5 1o 15 20 25 Figure 4.53. The time responses of the fourth synthesized fixed-interface linear mode displacement n?” for the three-DOF nonlinear-CMS model and the forty-one-DOF linear- CMS nonlinear model of substructure 'y. - - LCMSLPrb, 3 DOF 0.06 ~ LCMS, 3 DOF W a- - NLCMS, 3 DOF O 04 , Original Model, 41 DOF . A 1 1 1‘ ' I \ D , .\ ll ' \l‘ l . {l r" ‘3 PE ll fl P 1"; [la-I rlll [A (m (ill, Ill-.1 If 0.02". '1 11. [l 1" l‘ .13 ll '? 1" .l '1 '1‘! a? ',‘ 3" 5' l, Ml. IT I 1’ '1 .’ l f l : l ,.' i j l f l I} l ,' 3 l l‘ ‘4 ll 1. I' If” “ l, “ l' "'I -. 1'. 91 n1; 111. s 1*{“1H.11H‘:.tl".l'w — o~"’i=ifi'!'i!;i s!ir-mu..."!.w".,«"..1"~u'u " 1’ «w fluvial-1:: szwzwwm"w * l l I l l l ’ | " l i i l l ”I i 'I ’5 L “ '1 ‘ u ‘ l l l. ‘- l' l l‘ ‘1 ~: I' 'l *‘ " ‘r’l‘l’l'l‘m‘ L4! 'fll'li-‘J ' 1 l I l ‘ , | .3 P 1 - J l ‘l l l. l il l 1‘ ‘I ‘l -002 _ 3‘ ‘1 l l l 1' l l, l. 1 '1 l .., 'l 3‘ 1" l {:1 l" l! '3 11“ 5 1‘ u ‘11" ‘11:" “14‘ a if “v >4 if v a" u" 4 M ‘3’ “1’ ll ‘. ’ iv” ‘o" -0.04~ ‘ -0.06>— 1 O 5 1O 15 20 25 t Figure 4.54. The time response of the displacement 11:1 of mass m1 from Figure 4.10, which corresponds to x?,xlfi, and :13? from Figure 4.11, for the three-DOF linear-CMS linear model, the three-DOF linear-CMS nonlinear model, the three—DOF nonlinear-CMS model, and the forty-oneeDOF original model. 194 T 0.5 ' - - LCMSLPrb. 3 DOF LCMS. 3 DOF 0.4 *‘ - - NLCMS. 3 DOF H O 3T Original Model. 41 DOF 0,2 I‘ . \ I A \ / I/ \ r _\ / f I I ’ I I‘ I ’l‘ \ ‘ 1‘ I, \ 3’ , I O 1 P I H \\ I ,I' " I. ,l \l‘ I/ ll ‘9‘ \\ ’I ‘ \ I, \ I1 ' I I l \ , ‘ I 1 I‘ ’ ‘l I I I ‘ ‘I t \ I l I ’ \ I , II ‘ l ' I‘ ’\ . ' l ‘ ‘ . N O P l I ‘ I I l I I u I ‘ -‘ II \ \ [ \ .I I x "I I " I ' ‘ ’ ‘-’ I" V“ I I, ‘ ‘ ‘ l ‘1 l: I ’ 'H I I .' l I l ‘1' I,‘ I ' '\ I \ ‘ f \ , II ’A -01 I." I " ‘ \I I I I l .7 \l I! I, ‘ I] ‘ l / l I! ‘ ‘ I " J ’ x ’ " I x ,I' "x. I ‘w ‘ ‘\ \, \ I ,2 “0.2 F, \ , ’ \ x, \.\ I \ I _, _O.3 *- .4 "'O.4L -1 —o.5 » . 1 A A . 0 1O 20 30 40 50 t Figure 4.55. The time response of the displacement 1:2 of mass m2 from Figure 4.10, from Figure 4.11, for the three-DOF linear-CMS linear model, the three-DOF linear-CMS nonlinear model, the three-DOF nonlinear-CMS model, and the which corresponds to a: forty-one-DOF original model. — - LCMSLPrb. 3 DOF 0.5 * LCMS. 3 DOF ‘ —- - NLCMS. 3 DOF (Ml Original Model, 41 DOF . 0.3 ~ - 1.4 ,‘1 J ’1“ ’- \ I ‘ ‘ 1 '~ ’ \ I ’ - . 0.2 " 1., ’ i 1‘ I’ K. 'I ‘ I ‘ I ‘, I ‘\ ‘\ v\ 'H I \ I ‘ , J l I 0 1 _ I ’. 1 I ‘ ’ \‘l ' ‘I ‘1 'l l \ ' | I I . I \f I I I y 3 l I \ I | I I I \ ‘I I I. , , \\ , ’I. 1 .1 ‘ ‘ I ' I & O_I, ‘ ’\ ‘ " I ll ’ \l ’ 1‘ I” I " l ‘I .1 ' I" x ’ l I ‘ I I , f I \ I ' ‘ ‘ I 5 5 ', I I I ’ I ‘ \ I H I j I ‘ I ‘ l 1 'I -0.1_” I ll! ’l l ‘ I l I \l ,— II ’ \ ' 5 A ' .‘ “ I" I ‘ " h 1 ‘1 I, I I‘ I , , C ' \ I l 1 ‘ I I \ I . —o-2 I \ /‘\ ‘ I “\I / ‘J \ l I" \,’r \ ,I ‘ ’ -< 4. / I f A " I ’ \ U \ \ -O.3 r ‘ —OI4 _ - —O.5 * ‘ I l I l 0 1O 20 3O 4O 50 t Figure 4.56. The time response of the displacement 3:21 of mass "121 from Figure 4.10, which corresponds to 33% from Figure 4.11, for the three-DOF linear-CMS linear model, the three-DOF linear-CMS nonlinear model, the three-DOF nonlinear-CMS model, and the forty—one—DOF original model. 195 o 4. - - LCMSLPrb. 3 DOF I LCMS, 3 DOF -‘ - NLCMS, 3 DOF p 0'3 b Original Model, 41 DOF _ ‘N I /. -_‘> '\-‘ 0.2 I, I \ \ A U‘ I \\ 1“ I \\ f I \ / .I I. \I \ I r . _l \ r I ‘ I ‘I‘\ V ‘ \ ’N ' ' I ’ ’ ‘ ' \ I ‘ A l I \ , f l l 01—” “‘~ ~* M,‘ "I, u, u, 1,, . . [I \ I 1 | ‘ I I I II ‘I I .I I ll ‘ II I . I I ‘l ‘ I I .I I h I ‘ l I, . 1 I I ‘ II \ I I I ‘ I. I I l, I I ' I4 I I I, i . g . I l I l I l 1‘ I I f I . ' I I l I; l 0 — ’ l , I - .I ' ’ I I ‘ I H X I I s i H l ' . '\ I \‘ I I \ I l I I ‘7 l I. I I II I I I I I I I. . I ‘ I I , I‘ I II l. Io I II 5 II I I‘ l I I. I \ I‘ II 5. ‘0'1 "7 I i I ‘ ‘ III " A 7 ‘I- I ’ I1 I n‘ I '~ ‘ i i- '. ’I 4+ "I \ I \ \I ,, i I “I I ‘1' I.‘ I" \I \ \II I I l ‘ I - I k ' \ . \ / I K \ ‘ l \ I] i I J \\,/l 1 ‘ L -0.2I7 \ I \I I \ I \‘l' \ ‘1 _4 -0-3I I —o.4L _ O 5 1O 15 20 25 t Figure 4.57. The time response of the displacement 3:22 of mass mm from Figure 4.10, which corresponds to at; from Figure 4.11, for the three-DOF linear-CMS linear model, the three-DOF linear-CMS nonlinear model, the three-DOF nonlinear-CMS model, and the forty-one-DOF original model. J J 0.5 r- - - LCMSLPrb, 3 DOF F LCMS. 3 DOF 0.4 ~ - - NLCMS, 3 DOF ~ Original Model. 41 DOF . -J O . 3 r- \ I" / I~ I‘ ! {‘3 F & I‘r‘ 02L- " ’ \ III yo‘ ‘i IfiI‘I II“ ”II III" \ (II III II ' ,, "I l i f ’. I ' ‘ l - \ ‘. I 4 I I I ‘ I I ' I! l I ‘ I Q ' ‘ ‘ I II I‘ \ yr 1 I l ,‘ I .‘ I I ‘ I \ I I ‘ I, In \‘ {I I O.1~ . ,. I g , x | , e I , ,. \. I I l: l ,3 r I ‘I l '1' i ‘ I . I « ‘ ' I ' I l I I ‘, I I f l ._ ,~ . . ,, . \ z I I I \ t I . . I '4 , V o I' I I I I ' I. ‘ I I l\ 1 I ' l l I I ‘ I‘_ l II LI >< ‘ \ l I \ \ l I , l I \ I I i I I ,I I 1. ,I‘ l I! . \l I I V II , ‘ l \I 'I’ ‘ ‘I‘ 3‘ I: ‘ "0-1 ”I I ‘ \ ‘. I l I I ll ‘ ,' l I \I i . 1 ‘ “a, I . 1 h I . 1 I, V. : I I l I l I ‘ I — HI, ‘ I I‘ J I I l $ l. I \ I , ‘ vI- _ 0.2} I) \ gyI 1‘ _ \\_l " ‘ \I ‘I v] ‘I \l \ i . \ ‘1 k‘ l_- ll \’ l; ‘I J . 4 —O.3 r _ -0.4 r ~ —O.5 * -* L 1 1 l 0 5 10 1 5 20 25 Figure 4.58. The time response of the displacement 11:41 of mass m4] from Figure 4.10, which corresponds to 11:31 from Figure 4.11, for the three—DOF linear-CMS linear model, the three-DOF linear-CMS nonlinear model, the three-DOF nonlinear-CMS model, and the forty-one-DOF original model. 196 — — NLICMS, 3 DOF 1.5— — LCMS. 41 DOF I .A 0‘ O 1 O 20 3O 4O 50 Figure 4.59. The time responses of the first synthesized fixed-interface mode displacement 77;”) for the three-DOF nonlinear-CMS model and the forty-one-DOF linear-CMS nonlinear model of substructure 5 with niv‘fim) = 1.3, ”{mm = 1.3, 170(0) = 0.04, 7796(0) = .N, . 0, 721 7(0) = 0, Wm = 0. - — NLéMS. 3 DOF 1.5~ -—~— LCMS, 41 DOF f l 1 , ' \ >- \ 2 v- O*-\ \ ‘ l \ —o.5r \‘ ‘ \ _1_ J —1.5~ 1 o 5 1o 15 20 25 Figure 4.60. The time responses of the first synthesized fixed-interface mode displacement nix/‘7 for the three-DOF nonlinear-CMS model and the forty-one-DOF linear-CMS nonlinear model of substructure 'y with niv’flm) = 1.3, ni’vq(0) = 1.3, 77C(O) = 0.04, 77?,”8(O) = 0, 7'2{V”(0) = 0. THO) = 0. 197 4.9 Appendix 4 The component forms of the mass and stiffness matrices and the nonlinear force vector of the class of spring-mass systems from section 4.5 are given by m- if j = 2' 7710' "—2 2 (4.69) .1 0 otherwise ' 198 (Na-1) k1 + Z k(23Nrnass+2) if i: 1’ j: 1 3:0 —k(2-9Nmass+2) If 2: 1, j : SNma33 + 2, S : 0,1,2,...,(Na *1) —k(231lvfll(188) If I: S Nrnass +1, J :2"1, 3 :1,2,...,Na k(291v'l7l(188) if 7: 2 S anass + 1, j Z 7;, S :1,2,...,Na _k(2sifvn1088+2) if I = S NTTIOSS + 2, j = 1, s =1,2,...,(Na,— 1) k(2anmss+2) + k(2staSS+4) if i = SNmass + 2’ j = i, s = 1,2, ..., (Na — 1) (4.70) 412an....” +4, if i = s Nnms + 2, j = 2' + 1, s =1,2,...,(Na — 1) _k-‘(2z'—2) if (2 g 2' S Nnmss) or (s Nmass + 3 S i S (s + 1) Nmass), j=i— 1, s: 1,2,...,(Na— 1) k(2,_2, + km, if (2 S i S Nmass) 0r (8 Nmass + 3 S 2' S (8 + 1) Nmass), j: 2', s =1,2,...,(Na — 1) —k'(2,-) if (2 S i _<_ Nmass) or (s Nmass + 3 gig (8 +1)N7TIGSS)7 j =i+1, s =1,2,...,(Na — 1) 0 otherwise 199 f (NW—1) 4 ,. _ 3 — Z: A(2 8 N77103:; +3) (’L(«5' 4N‘TI‘L(£.S’8+2) $1) 320 if i: 1 [C(23Nmass+1)($i - $(i—1))3 if i : SNmass +1, 5 =1,2,...,Na ”2‘ MID , 3 3 k(231\77nass+3)(xi — I1) — k(2sta33+5)(~T(i+l) — 12,) if l28N771088+23 S 21,2,...,(Na_1) k(2i—1)($i — 1*(i—1))3 “ k(2i+1)($(i+1) _ $223 if i = 2, “'3 Nmass, 01' i = SNmass + 3, ..., (3 +1)Nmass, s = 1,2, ..., (Na — 1) L for i,j I 1, 2, 3, ..., (Ara [\r'NULSS +1). The component forms of the mass and stiffness matrices, the nonlinear force vector, and the reaction force vector for substructures S are given by 0 fiizszl my: 7% iszi M73 0 otherwise 200 1:25 if 2f: 1, j =1 —k§ if i: 1, j =2 _kig2Nmass) if i : Nmass + 1’ j : i - 1 k5 nizN..+Lj=i 1?} = 1 (”ma“) (4.73) kéi—Q) + kg?) if (2 S 1 S Nrnass), j :2 ‘kén if (2 S i S Afmass): .7: i +1 k 0 otherwise r —k§(1‘§— 1‘?)3 if i=1 5' ' . . {12' = klséNmass-l-l)(rig _ xii-UP 1f 1 = Nmass + 1 (4'74) .5 ,.S S 3 S .S' S 3 ' - __ \ A‘(2i—1)(lz’ — 10.4)) - k(22‘+1)(l(i+1)— 51:2.) 1f 2 —— 2, ..., Nmass ( F5" 1f 2: J f? = (4.75) 0 otherwise for 77,1 -_—. 1, 2, 3 (NW, + 1), and s = a, 7, 5, . Note that mff,’ :- . 1’3 __ , .13 _ f3 _ ,. .13 _ ., 5 ... mg, 7713 — m3, mleassd'l) —— m(2’Vnm.ss+l)’ k2 — k2, k3 —— k3, ..., k(2Nmass+l) _ 4 '7 _ .7 _ ,. '7 _ A’(2N7nass+1)’ "L2 — "2(Nmass+2)’ m3 _ mmeass—f3)’ m’(N-mass+1) _ ”1(2Nma33‘f1)’ 7__, 3_. .7 _ , . ,.7 _ , k2 *— k(2Nmass+2)‘ A“ k(2Nm(133+3)’ "" k(2Nmass+l) _ k(4Nmass+1)’ and SO on. The coefficients of the pqr cubic nonlinear terms expressed in terms of fixed-interface linear modal coordinates and associated with the ith mode of substructures S are given by 201 bfV,S ipqr N,S N,S N,S N,S.xr __ bipqr + biprq + birpq I‘M" _ N,S N,S N,S bipqr + biqpr + biqrp N,S N,S N,S N,S N,S N,S bipq’r + biprq + biqp'r + biqrp + birpq + birqp N,S ‘ where bipqr are N S (Ninass+1) N S N S ipfzr Z (phi, bh-Pb" h=2 for S = 13, 7, 6, ..., where (1)295 is the hi component of the modal matrix defined in section 4.3.2.1, and b are given by 202 if p = q = 7‘ if p = q (4.76) if q = r otherwise, (4.77) N,S hpqr 7N.S h pqr — The coefficients of the pqr cubic nonlinear terms expressed in terms of the linear N,S N,S .N,S f ‘5 .N,S ,N.S N,S 5 «NS .N,S [NS .3 - [A C93p @311 993,, + 3 k5 $31) Cb3q CZ)21' 3 02 ¢2q ¢2r — 4‘5 .S .N.S N. .N,S _S (N,S N,S N,S -3 A5 993,, go'Zq (927' + [‘5 ¢2p cb2q d)27‘ if h=2 N,S N,S N,S ,3 .N,S .N,S N,S ,5 [k )Cbhp Cbhq Qbhr —3A'(2h—1)¢hp Cbhq ¢(h—1)1‘ (2h—l .S .N,S .N,S .N.5 _ ,S .N,S N,S N,S +3 A'(2h-—l) Cohp ¢(h—l)q GUI—UT A(2h—1) "leirovitch, Computational Methods in Structural Dynamics. Sijthofl Noordhoff, 1980. [89] J. Carr, Applications of Centre Manifold Theory. Springer-Verlag, 1981. [90] L. Meirovitch, Analytical Methods in Vibrations. Macmillan Publishing, New York, 1967. [91] M. P. Castanier, Y. C. Tan, and C. Pierre, “Characteristic constraint modes for com- ponent mode synthesis,” in Proceedings of Design Engineering Technical Conference 1909, no. DETC99/VIB-8187, September 1999. 217 IIIIIIIIIIIII UNIVERSIT [[1]] [[l]]]]]]] [[1]] 3 1293 0250