x4e... @th n “rut: 2.. . “adv... . ‘tIDAKr , ~rV-Vy~,. .110 Il‘v‘(‘| . CH?" VI. . ,u ....I, , 5.. ~th 1.4.1.2.»! . 1.4.3.4 PEP. ...... t #5... .22 ..s; ... ) .‘ .. ... . . . . I U n n J .. ;.. . ‘ . .. . .... .. u.» «......i : P‘YH‘WaLvAL .. .. 3.-. . o .1 F ..lu Sun: turn. ...w. .u o I .. . . I .52.... luv... ....w.n!.w-:m. «3%.... . t n a i , {Ifi 9.: bung“... 5.2611 (1.34.. . .‘f, £1; ‘\r v o iv; .....z..-h; .95: sirutfislvflis n J.. . I. ‘ t v t ‘. 4‘ ‘1. w . . ‘. , . . a... . ...- . . . . , .75.:in . . 1 - _ . u ..s $90k. «nah... h x . , V - 4 . .. . ..., ...: : . . YI- .l\: ..7.\. IINIVAI o . unfit?! .54”. 1-7.4: .... , ‘ ‘ . £323.}: 14.). 3.... .. .. z , xiurmr.nu$1.u..iu ....mrek . u...u..L{‘t:h3:.wu: v}... I! ‘ . . I . is .r.,,. 4!? ..... . . . , .. ‘2... )2 .--?!xl... . . . . . . . ...: ...... a zeiifitéfi . . . i BMJI. £1-12! v..u3-y.;g..».‘.l..f&. a . . ....Vf.....i...;! l l,,,£cu.l:...u! ...“..l. .. x ‘. :.\ia.i2..xa!1_4 ......u-l.do$to.rl... ..ofiuaniIllb . 55.33.! -hltll}un, ..l Sit-l..!-l . ...lcltlwxolnuhnu I €5.33S...32.5.1... I ..... $5.11.- .brufir! I! ‘II “44‘: ‘nnrlv. : 5.9114111!“ |plll§£3l§ll¢l|£vIt-'\!J; ‘ ii Iii. u ti-t.f!:.1 -..1(.£12..Li¢|;!11«l.9llhuufl 4 ii . . 11vi$w§l~¢£33 . .. ‘ -58 ...}.av. ‘-:.£a..u...fi.~is..u . a ..o I x... 9,3. fi.xzs.lx...lx.‘ a Iltvlipr‘.’ . . ,. v 5.2.2. . 3?? x . . ‘ .. . ... ... u i: . ......H .. . .. -c 15,..- ao i. . .. ‘1.-- . ...) ,. 2’: .53‘1hbnlhtfl. ......th ..qu\:....1... . t; . i\\ 7‘ (tODXI’bQI‘Il! ' 0" "|V|: II...) n 5...» . .lil: . . . . :1.-. . hunt .... r. . 4.... ...dfl . . furrow: ..is ... . .. , I} «‘11. . . a , 4 illiill Ill-A . innailwflnin... 4 x: . .33.! ...: . 1.. A x . . x . 3 {.ur... . Wr‘hn‘i I x I f . z . .. .....u‘ p u. . . .... n x . sillffrbi ,v v .flléiot. .3 I ...:i. ... . ‘ - rm“ . $5: I). x. ‘ h}. . » .chl .. ail... .. 41.3.... «Xv... «.... Aliit ..-. .... ..- -.. u. . . a w 1. «- .D‘l’-;l'il t3}; . .5 t ‘ 4 >- v ‘c . I ' ...Dno'V'..nnbu f D: HIU‘ vst' ID, .-l\ .17 1c.-IUV.t t‘HJ‘M-thvval.“ Ad Iw‘OID. - |.lfvrh Iiiffluf In -....11 6 ‘ . it... - ..- A. .- (1.3”: -I‘IOAUIP - -3. £uill ‘ w I . L15- ~‘L . i. > ‘ ‘I .. .. Brawl”... ... ,. .. n v t. ‘ I‘ . . . ... . w ... ‘ : , I r \ . A ¢I A 1 ‘ n “V ‘ . . . ...: a 1. v hr: I .~ . ’1' Il’.‘ v : ....il. 6‘22}. .....obhmk. :3: 31-... 0-.A‘v.v. . l’..v|v . u . ‘ .. ‘ tabuif. .....Pwa‘knu....o¥w2ufl. $3.... .mvfl? . U. I§fufl fir..- I IN.” I u; 10. VI.- 9. . . xlnfl..w..u.l::.il. flirtafil: .. O . v I“ o 4. IV . l v ...v . v n ‘ . l . I u I I I v 2 o .. . n 9.71. -99....19. Fr a...-..u..,..u 432...},34 :{31‘2ilxnhh . .. 21.25)} ......... eh..:.....n33 .. . I , .ye, . .--..riil .....t! L». It‘ll '. ~I bl.- 3|." 0 3 ’. 71.an FF. fl.) A 2 .... 7:! llllilllllllllllllllllllli|i|||||lHlllllll|||||||llllllllllll 3 1293 01568844 LIBRARY Mlchlgan State Unlverslty This is to certify that the dissertation entitled AN "OPTIMAL" MODAL REDUCTION OF SYSTEMS WITH FRICTIONAL EXCITATION presented by Ramana Kappagantu V. has been accepted towards fulfillment of the requirements for Ph.D. degree inMechanical Engineering I Major profess}r Date 3/316? MS U is an Affirmative Action/Equal Opportunity Institution 0-12771 PLACE Ii RETURN BOXto mnovothb chookouiirom your record. TO AVOID FINES Mum on or baton duo duo. DATE DUE DATE DUE DATE DUE MSU is An Affirmativo Action/Equal Opportunity Institution . Wan-9.1 AN “OPTIMAL” MODAL REDUCTION OF SYSTEMS WITH FRICTIONAL EXCITATION By Ramana Kappagantu V. A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mechanical Engineering 1997 ABSTRACT AN “OPTIMAL” MODAL REDUCTION OF SYSTEMS WITH FRICTIONAL EXCITATION By Ramana Kappagantu V. We apply a method of order reduction in modeling multi-degree-of-freedom and continuous systems with frictional excitation. We use proper orthogonal decomposi- tion to obtain a basis that has the optimal distribution of energy in the system. The system dynamics are rich, including a variety of periodic, quasi-periodic and chaotic responses. Each of these responses could result in a different set of proper orthogonal modes. We select the proper orthogonal modes that can be used to build the reduced system model by using the chaotic dynamics. The validity of the method of reduc- ing the system and representing the system with this reduced model is studied. In these studies we considered a numerical example and an experimental system. The numerical example is a simple chain of spring-mass-damper system subjected to fric- tional excitation from a moving belt. The experiment consisted of a cantilever beam subjected to friction from an oscillating surface. The results confirmed the validity of the method and exposed the benefits of using proper orthogonal modes rather than linear natural modes of the system. Copyright by Ramana Kappagantu V. 1997 To my parents and family. iv ACKNOWLEDGEMENTS It is a proud feeling to acknowledge with thanks all the guidance and encourage- ment given to me by my advisor Professor Brian Feeny. It is my good fortune to be associated with this wonderful person. His simple thinking, penchant for new ideas, foresight and experimental bent of mind are great sources of inspiration. I express my gratitude to Professor Charles MacCluer for introducing me to the world of nonlinear dynamics and chaos and Professor Alan Haddow for initiating me into the experi- mental research in these areas. I thank Professors Steven Shaw and Clark Radcliffe for their constructive criticism of my work. I am privileged to have them all on my committee. Professor Haddow deserves special thanks for using his good offices in procuring the equipment essential for the experiment. I thank Professor Thompson for his as- sistance with the strain-gage conditioners. I appreciate the help of Professor Kalinath Mukherjee and his student Satyanarayna Kudapa in laser fabrication of some of the equipment needed used in our preliminary research. I acknowledge the help of Roy Bailiff and Bob Rose in the experiment-setup. It is my most humble duty to thank all my teachers who taught me the basics in Mathematics, Sciences and Engineering. Professors Hassan Khalil, Sheldon Axler, Brian Feeny and Charles MacCluer of Michigan State University, Professors Ashitava Ghosal and Udipi Shrinivasa of Indian Institute of Science and Professors G. Rama Murthy, J. S. Ansari and J. V. Subrahmanyam of Osmania University deserves a special mention in this connection. I cherish their excellent teaching. I will be failing in my duty if I do not acknowledge the much needed financial assistance given to me by Professors Stuart Cage and Brian Pijanowski of the De- partment of Entomology. My association with them gave me an exposure to the real world problems on the non-engineering side. Chris Hause, Wei Li, Jin-Wei Liang, Shyh-Leh Chen, Ching-Ming Yuan, Chang- Po Chao, Choong-Min Jung, Yihong Zhang were a great team in the Dynamics lab- oratory whose association was of continuous help in keeping my spirits high. They were instrumental in creating the right atmosphere for research in the laboratory. Mr. Chris Hause deserves special thanks for sharing his skills in the conduct of the experiment. vi Contents 1 Introduction 1.1 1.2 1.3 1.4 1.5 1.6 2.1 2.2 2.3 2.4 2.5 3.1 The Problem Definition .......................... Motivation ................................. Literature Review ............................. 1.3.1 Self-excited oscillations and dry friction ............. 1.3.2 Beams subjected to frictional forces ............... 1.3.3 Phase-space reconstruction and dimension estimates ...... 1.3.4 Review of current research in POD ............... Proposed Research ............................ Contributions ............................... Thesis Organization ............................ Proper Orthogonal Decomposition (POD) Definitions and Theory .......................... The Method ....................... i ......... Relation to Normal Modes — Unforced and Undamped Oscillations . Relation to Normal Modes in General Systems ............. Recent Engineering Applications of POD ................ POD of fiiction Induced Vibrations: A Numerical Experiment System Description and Modeling Issues ................ vii (Dmrprliwi—‘H 12 13 13 14 16 16 19 20 21 22 24 25 3.2 Higher Order Dynamics, Chaos and its Verification .......... 3.2.1 Mutual information ........................ 3.2.2 Embedding dimension ...................... 3.2.3 Lyapunov exponents ....................... 3.3 Finding POMs and System Reduction .................. 3.4 Modal Reduction ............................. 3.4.1 The mathematics ......................... 3.4.2 The order of reduction ...................... 3.5 Validation of Reduced Model ...................... 3.5.1 Qualitative comparison ...................... 3.5.2 Quantitative comparison ..................... 3.5.3 Comparison of bifurcations .................... 3.5.4 Unstable periodic orbits ..................... 3.6 Discussions and Conclusions ....................... Physics and Measurement Methods of the Experiment 4.1 Description ................................ 4.2 System Characteristics .......................... 4.3 Friction Measurement ........................... 4.4 Displacement Measurement ....................... 4.4.1 From strain to displacement ................... 4.4.2 From discrete to continuous ................... 4.5 Measurement of Normal Load Variation ................. Beam Dynamics, Reconstruction and POD 5.1 Introduction ................................ 5.2 Bifurcation Studies ............................ 5.3 Characterization of the Dynamics .................... 5.3.1 Reconstruction Results ...................... viii 31 34 36 39 42 43 43 45 48 48 49 5O 51 53 54 54 58 60 68 69 73 75 77 77 78 82 86 5.4 Proper Orthogonal Modes: Giving a Spatial Description ....... 88 5.4.1 Proper orthogonal modes: using a larger number of sensors . . 92 5.5 Discussion ................................. 94 Mathematical Model, Reduction and Validation 97 6.1 Mathematical Model ........................... 97 6.1.1 Coulomb-friction forcing with infinite contact stiffness ..... 100 6.1.2 Finite contact stiffness with Coulomb friction excitation . . . . 103 6.1.3 Effects of normal vibrations of the beam on normal load . . . 105 6.1.4 Modal reduction using POMs .................. 107 6.2 Numerical Solution Technique ...................... 109 6.3 Reduced Model — Numerical Simulations ............... 110 6.3.1 Friction-force characteristics —— different scenarios ....... 110 6.3.2 Frequency sweeps ......................... 113 6.3.3 Different levels of reduction using POMS ............ 114 6.3.4 Finite and infinite contact-stiffness friction excitation ..... 119 6.3.5 Normal vibrations included ................... 120 6.3.6 Amplitude sweeps ......................... 122 6.3.7 The effect of static friction coefficient .............. 123 6.4 Conclusions ................................ 125 Conclusions 127 7.1 Summary of the Results and Conclusions ................ 127 7.2 Future Work ................................ 130 ix List of Tables 3.1 3.2 5.1 5.2 Statistical deviations at different levels of reduction .......... 47 Quantitative comparison between full and reduced models. AX and AE defined in equations (3.37) and (3.36) respectively ........ 50 Contribution of LNMs to POMFs at 12.0 Hz. ............. 89 Contribution of LNMs to POMFs at 41.1 Hz. ............. 94 List of Figures 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 Configuration of spring mass damper system under frictional excitation 25 Bifurcation diagram of the full model and a few different possible or- bits. The belt speeds for these orbits are 0.20, 0.22 and 0.25 ..... Average mutual information of the 10~SMD system, using velocity time series of the first mass. The numbers of cells used in estimatin prob- ability density are 5, 10 20, 40 and 80 and the corresponding p ots go towards the axes. The first minimum occurs near T = 11. ...... False-nearest-neighbors in the lO-SMD system, using velocity time se- ries of the first mass. Note the first zero occurs near d = 6. ..... A schematic representation of the evolution and replacement rocedure used to estimate Lyapunov exponents from experimental ata. The largest Lyapunov ex onent is computed from the growth of length elements. When the ength of the vector between two points becomes large, a new point is chosen near the reference trajectory, minimizing both the replacement length L and the orientation angle 0 ....... Typical orbits of the 10-SMD system in the first row and corresponding dominant POMs are shown in the bottom row. The first two orbits are obtained at a belt 8 eed of 0.25 with different I.C.’s and the third obtained at a belt spec of 0.20 ...................... Phase portraits of 10—SMD s stem, with different levels of reduction — A qualitative comparison or determining the order of reduction Qualitative comparison between full and reduced models. Top row corres onds to the dynamics of reduced model and bottom row to the ful model. Each column corresponds to a set of identical initial conditions. All other parameters are kept constant. .......... Comparison between full and reduced models: Bifurcations ...... xi 31 36 38 41 43 46 49 51 3.10 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10 4.11 4.12 4.13 4.14 4.15 4.16 Unstable periodic orbits in the chaotic attractor reconstructed from the velocity time series of the first mass. The first row corresponds to data from reduced model and second from_full model. The first column shows part of the attractor. The remaining shows approx1mate unstable periodic orbits of periods 1, 2 and 3 .............. Schematic diagram of the experimental setup ............. A photograph of the experimental setup ................ The effect of tilt on the force transmitted ............... Average power spectrum response to impulses close to the clamped end of the system with zero normal load at the contact. .......... Average power spectrum response to impulses close to the clamped end of the system with high normal load at the contact point. ...... Configuration of friction measurement setup .............. Sample time-series of displacement, velocity, friction force and effective normal load ................................ Friction-displacement relation ...................... Friction-velocity relation ......................... A schematic diagram showing a massless compliant-contact model . . Sample time series of direct (laser) and indirect (strain gage) displace- ment measurements (Y and y) at an operating frequency of 12.0 Hz. Fourier transforms (averaged) of direct (laser) and indirect (strain gages) displacement measurements at an operating frequency of 12.0 Hz. .................................... Sample time series of direct (laser) and indirect (strain gage) displace- ment measurements (Y and y) at an operating frequency of 41.1 Hz. Fourier transforms (averaged) of direct (laser) and indirect (strain gages) displacement measurements at an operating frequency of 41.1 2. .................................... Fourier transforms (averaged) of the noise floor levels for direct (laser) and indirect (strain gages) displacement measurements ........ Relation between normal load and position of the beam tip on the slider. ................................... xii 52 55 56 57 59 60 61 64 65 65 66 72 74 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 5.10 5.11 5.12 5.13 5.14 Amplitude versus frequency of the shaker with other parameters kept constant. A tenth-degree polynomial fit has been used to curve fit the data. Input frequency sweep of the shaker thus takes us along the above curve on the output. ....................... Frequency sweeps (forward and reverse): the beam response is denoted by the strain at the location 6 and measured at uarter phase of the input acceleration. The bifurcation parameter is t e input frequency. Finer fre uency sweeps around 12 and 41.1 Hz : the beam response is denote by the strain at the location 6 and measured at quarter ghase of the input acceleration. The bifurcation parameter is the input equency. ................................. Amplitude sweep at 12.0 Hz: the beam response is denoted by the strain at the location 6 and measured at quarter phase of the input acceleration. The bifurcation parameter is the input amplitude. Amplitude sweep at 41.1 Hz: the beam response is denoted by the strain at the location 6 and measured at quarter phase of the input acceleration. The bifurcation parameter is the input amplitude. Wear sweeps at 12.0 and 41.1 Hz. All other parameters (under our control) are kept constant. ....................... Delay maps at a few selected frequencies on the bifurcation-diagram elucidate the periodicity in the beam. The frequencies are: (a) 6.3, (b) 10.0, (c) 12.0, (d) 26.0, (e) 37.7, (f) 39.0, (g) 41.1 and (h) 47.0 Hz. Dynamics of the contact point of the beam. The top row indicates the base-portraits of the contact point at frequencies 11, 12 and 41.1 Hz. he bottom row indicates corresponding power-spectra. ....... EFT of the strain signal measured at location 6 on the beam at 12.0 2. .................................... EFT of the strain signal measured at location 6 on the beam at 41.1 2. .......................... - .......... FFT of the noise level of the strain signal measured at location 6 on the beam. ................................. Proper ortho onal modes of chaotic dynamics at a driving frequency of 12.0 Hz. . he continuous versions of the discrete modes are not orthonormalized. ............................. POMFs (obtained by orthonormalization of the curves fit on POMS) from chaotic dynamics at a driving frequency of 12.0 Hz. ....... Proper orthogonal modes of chaotic dynamics at a driving frequency of 41.1 Hz. .The continuous versions of the discrete modes are not orthonormalized. ............................. xiii 79 80 81 81 82 83 84 91 5.15 POMFs (obtained by orthonormalization of the curves fit on POMs) of chaotic dynamics at a driving frequency of 41.1 Hz. ........ 5.16 Proper orthogonal modes of chaotic dynamics at a driving frequency of 12.0 Hz, using 41 pseudo-sensors. The continuous versions of the discrete modes are not orthonormalized. ................ 5.17 POMFs of chaotic dynamics at a driving frequency of 12.0 Hz, using 41 pseudo-sensors. ............................ 6.1 Simplified model for normal vibrations ................. 6.2 Friction characteristics of the infinite contact stiffness model. The top row shows the case where p5 = pg and the bottom row shows the case where p5 > pg .............................. 6.3 Friction characteristics of the compliant contact model with p5 = pg. The top row shows the case of macroscopic stick-slip and the bottom row shows the case of micro-stick .................... 6.4 Friction characteristics of the com liant contact model with p5 > pg. The top row shows the case of stic -slip and the bottom row shows the case of micro-stick ............................ 6.5 “Frequency sweeps” of reduced models obtained at different levels of reduction using POMFs (extracted from experiment at 12.0 Hz and the continuous versions orthonormalized) superimposed on the exper- imental sweep. .............................. 6.6 “Frequency sweeps” of reduced models obtained at different levels of reduction using POMFs (extracted from experiment at 12.0 Hz and the continuous versions orthonormalized) superimposed on that of the lO-LNM model (“truth set”). ...................... 6.7 “Frequency swee s” of reduced models obtained at different levels of reduction using NMs of the cantilever superimposed on that of the 10-LNM model (“truth set”). ...................... 6.8 “Frequency sweeps” of reduced models obtained at different levels of reduction - from most dominant mode to all the available POMFs at 41.1 Hz (compliant contact model). .................. 6.9 Phase-portraits (contact point displacement versus velocity) from the three-POMF reduced model at 11, 11.9 and 42.9 Hz. ......... 6.10 “Frequency sweeps” of the system using the reduced model from POMFs (of the three dominant POMs obtained at 12.0 Hz). The top row cor- responds to the zero contact compliance friction excitation and the bottom row corresponds to non-zero contact-compliance. ....... xiv 93 95 96 106 111 112 113 115 117 118 121 6.11 “Frequency sweep” of the system using the reduced model from POMFs (of the three dominant POMs) obtained at 12.0 Hz. Normal vibrations are accounted for in the model. ..................... 6.12 Amplitude sweeps of the system using the reduced model from POMFs (of the three dominant POMs obtained at 12.0 Hz). The top row corresponds to the sweep at 12.5 Hz and the bottom row corresponds to 46.9 Hz. ................................ 6.13 ps/ pg sweeps of the system using the reduced model from POMFs (of the three dominant POMs obtained at 12.0 Hz). The top row corresponds to the sweep at 12.5 Hz and the bottom row corresponds to 43.9 Hz. ................................ XV Chapter 1 Introduction 1.1 The Problem Definition The focus of this research is the application of proper orthogonal decomposition (POD) theory to higher order systems subjected to frictional excitation. The goal is to build a reduced-order model for a continuous system using experimental data. In this process we wish to combine experiment and theory and exploit the rich dy- namics that are inherent of systems with frictional forcing and understand the spatial distribution of “energy” in the system. The key words here are proper orthogonal decomposition theory, friction induced vibrations and higher-order systems. In the next few paragraphs we will present mo- tivation, a brief review of the current-day research in friction dynamics and proper orthogonal decomposition theory followed by proposed research and thesis organiza- tion. 1.2 Motivation Friction is highly complex nonlinear mechanism that dissipates energy. It exists in almost all mechanical systems involving relative motion. When two surfaces in contact slide relative to one another, self-excited vibrations (Nayfeh and Mook [62]) and stick- slip oscillations (Bowden and Leben [8]) are often observed. The energy dissipation from a system undergoing stick-slip oscillations can be in the form of sound which can be heard in daily life as squeak and squeal in automobiles and trains (Irretier [44], Schneider et al., [80], Van Ruiten [77]), music from a violin (Schelleng [79, 78]), etc. Of late manufacturing and designing high—precision equipment is becoming more and more necessary, and a need for a good model for simulation and control of higher order systems has become imperative. Though it is hard to estimate the expenditure involved in combatting friction-related problems, some in the friction community offer figures like 1.6 percent of the gross national product of developed countries, which amounts to $116 billion for the U.S.A. alone in 1995 [50]. Thus friction contributes to the economy both positively as well as negatively. Modeling of the frictional dynamics is therefore important for understanding the above mentioned phenomena. The need becomes more severe when the systems are of high order or continuous. Moreover the usage of dry friction dampers and fasteners to dissipate vibrational energy in turbine blades, aircraft frames and large space structures is becoming more and more common (Ferri and Bindemann [32], Ferri and Heck [33], Dupont and Kasturi [24]). Experimental and mathematical evidence indicate that friction induced vibrations involve interaction of degrees of freedom over broad range of spatial and temporal scales (Van Ruiten [77], Nakai and Yokoi [61, 60], Tworzydlo et al. [90], Arnov et al. [3], Schneider et al. [80], Irretier [44]). When the dynamics of higher-order me- chanical systems involving frictional forces need to be studied in detail, a thorough understanding of modal coupling and a credible model of friction are needed. . We realize that the actual physical processes are too complex to be described by mathematical equations. A great number of factors interacting in the real processes make such descriptions extremely difficult. However, as far as engineering applica- tions are concerned, good models in the form of mathematical equations would be of tremendous use, as they represent, at least approximately, real processes by ex- pressions that are easy to use, understand, communicate and improve. In nonlinear systems new theories are arising for solving problems analytically and sometimes pseudo-analytically. However most of the problems are still being verified by nu- merical methods. This is becoming increasingly feasible and important, thanks to the advances made in the current-day computing technology. There are several fric- tion models available in the literature in various forms. Recent survey papers by Armstrong-Helouvry et al. [2] and Ibrahim [42] gives a comprehensive list of friction models and methods used in tribology, lubrication, dynamics and controls. One can easily understand that one model alone cannot describe the different scenarios. Once when a friction characteristic is identified, the next question that comes up is how to model the other aspects of the system. Typically continuous systems are described by partial differential equations (PDEs). Due to infinite dimensionality, continuous systems are generally difficult to analyze. Analytic solutions to the governing PDE may not be available and boundary condi- tions may not be well specified. The problem gets more complex when the forcing is non-smooth, as is the case with friction and impact. In many cases, finite-difference or finite-element methods are used for numerical analyses. However they are often too complex and computationally intensive for on-line simulations. Different methods have been used in reducing the continuous / higher-order systems to low-order systems, such as finite element analysis, component-mode synthesis, and Galerkin’s methods, to name a few. Galerkin’s method of weighted residuals is the most appealing and treats a continuous system as an n-degree—of-freedom system by assuming the solution in the form of a linear combination of n known functions or assumed approximate modes. Typically these n functions are the linear natural modal functions. These methods become increasingly inconvenient as high-frequency dyanamics come into play because large amount of information (possibly redundant) is to be handled [26]. The dynamics of higher order systems subjected to frictional excitation generally involve high frequencies. As such these methods may not predict the behavior accu- rately. Another issue is the efficacy of using linear natural modes in determining the dynamics of a nonlinear system (which may have nonlinear natural modes [82, 83]). One might rather intuitively want to select the set of functions from a basis in which the energy distribution is optimal. If the major portion of the energy is con- tained in say three or four of the functions in this basis set, one can then pick those functions that have most energy. Towards getting the basis that has the optimally distributed energy we look at the proper modes of the spatial covariance matrix formed of the state time series (Cusumano [17, 18]). This method is analogous to the Karhunen-Loeve expansion which is widely used in pattern recognition (Fuku- naga [37]). The main advantage of proper orthogonal decomposition technique is that it allows the experimental identification of critical coupling information between various degrees-of-freedom. This in turn sheds some light on the noise generation and propagation mechanisms. Also POD “does not do the violence of linearization” as it gives a basis from the energy perspective [5]. Thus the usage of proper orthogonal modes (POMs) preserves the gist of nonlinearity even in the presence of a disconti- nuity [5]. Thus we are motivated to look into the feasibility of applying POD theory (which is a statistical formulation based on spectral theory and gives us a basis that has optimal distribution of “energy”) to systems subjected to frictional excitation. 1 .3 Literature Review 1.3.1 Self-excited oscillations and dry friction Most of the current-day research in the friction-induced dynamics deals with low- order systems (single or sometimes two-degree—of-freedom systems). This is mainly because friction is highly nonlinear and, even in the low-dimensional systems it is yet well modeled. Despite the clear differences in the circumstances (material properties, geometry, scale, etc.) at which the oscillatory phenomena of stick-slip occur it is a fact that most of them are induced by dry friction and that many of them have the common feature of rapidly alternating states of contact and no contact or sticking and sliding. Also the coupling of various degrees of freedom has been shown to play an important role in the generation of dynamic instabilities and various oscillatory phenomena in finite dimensional systems. There are many friction models available in the literature, of which the Coulomb friction law is the foremost and is a great simplification of the recent realistic mod- els. Though it realizes the existence of static and kinetic friction and is reliable in describing the macroscopic frictional behaviors like stick-slip (e.g., Den Hartog [40], Hundal [41], Shaw [81], Kato [57]) it is alone not able to model the subtle frictional features, like the frictional memory in sliding rocks (Ruina [76], Dietrich [20], Ibrahim [42]), rate and time dependencies of the static coefficient (Tworzydlo et al. [90]), the influences of normal degree of freedom in sliding friction (Oden and Martins [63], Tworzydlo et al. [90]). The first question therefore the dynamicist faces in modeling a system with frictional forces would be how to decide on the right law so that the system is appropriately modeled and the appropriate phenomena are noticed. Courtney-Pratt and Eisner [15] and Burwell and Rabinowicz [71] observed the creep sliding or micro-sliding of specimens for long periods of time after each increment of tangential force was applied. The experiments of Johannes, Green and Brockley [45] and Richardson and Nolle [74] indicate that p, decreases with the increase of the rate of application of the tangential force. For sufficiently small load rates, the coefficient of static friction is constant and equal to a value which is the usually quoted coefficient of static friction. For large loading rates the coefficient of static friciton tends to be constant and equal to a value which is usually interpreted as the coefficient of kinetic friction. Sinclair and Manville [84] showed that frictional vibration is due to the rise of the coefficient of friction as the slip speed decreases. While analyzing the experimental results of Bowden and Leben [8] and of Pa- penhuyzen [66], Blok [7] concluded that the essential condition for the occurence of stick-slip motion is a decrease of the frictional force with increasing sliding speeds. In the same paper, Blok also provided the first systematic study of frictional vibrations and established a quantitative criterion for their appearance. Although discrepancies between the results of frictional experiments performed on different apparatuses were noticed by various researchers (e.g., Burwell and Ba- binowicz [12], Rabinowicz et al. [72], Bowden and Tabor [9]) it was Tolstoi [88] who clearly stated for the first time the importance of the dynamic characteristics and vibrations of the apparatus. He considered normal vibrations to explain the sharp reduction in friction (when the sliding begins). He extended the theory developed by Koodinov [48] to explain self-excited vibrations as due to coupling between the normal displacements and tangential velocity of the slider under the conditions of hydrodynamic lubrication. Tolstoi concluded that the interface coefficient of friction does not explicitly depend on sliding velocity, and that the difference between the apparent static and kinetic coefficients of friction is the consequence of microscale vibrations accompanying frictional sliding. This observation was later confirmed by experiments of other researchers. An appropriate model for sliding friciton must incorporate physically reasonable normal contact conditions. In addition, all the major sources of coupling between normal and tangential degrees of freedom should be taken into account in such a model. The normal deformability of the interface is an essential feature of dynamic contact problems involving metallic bodies. Not taking into account this deformability in finite element models of these phenomena (even if frictional effects are negligible) leads to some serious physical inconsistencies. For example, the absence of normal deformability leads to models which can provide oscillations, depending on the mesh used, with frequencies as high as 106 Hz, but these models may be incapable of delivering experimentally observed contact oscillations of frequencies of the order of 103 Hz (Oden and Martins [63]). The ideas of Tolstoi [88] led Oden and Martins [63] to a new approach in the analysis of dynamic friction. They considered a relatively simple constitutive model of the interface, with a power law normal response and the coefficient of friction in- dependent of the velocity. This model was combined with an analysis of motion in particular of normal vibrations of the slider to give an apparent kinetic coefficient of friction different in general from the interface coefficient of friction. In their nu- merical studies they obtained a good qualitative modeling of general experimental observations. Later Tworzydlo and Becker [89] applied it to numerical modeling in the reduction of static friction by vibrations and found a very good agreement with experimental observations of Tolstoi. They modeled a typical pin-on-disk appara- tus using an extended version of the Oden-Martins friction model. They found the coupling between rotational and normal modes to be the primary mechanism of the self-excited oscillations. These oscillations combine with a high-frequency stick-slip motion to produce significant reduction of the apparent kinetic coefficient of friction. Bengisu and Akay [4], have carried out a stability analysis of friction induced vi- brations in multi—degree-of-freedom systems. In their studies they considered a contin- uous friction characteristic. They linearized the nonlinear system (though inherently linear the response is nonlinear due to the nonlinear nature of the friction force at the contact interface) for predicting the system response and qualitatively describ- ing bifurcations associated with the cross-over of the system poles. They obtained the solutions to the nonlinear equations of motion numerically and demonstrated the existence of limit cycles, quasi-periodic and chaotic behavior. New friction models are still being developed. This can be attributed mainly to the need to describe the phenomenon observed using the modern day sensor technology. For example Polycarpou and Soom [67] proposed a two-dimensional friction model which included both tangential and normal dynamics. This model was developed to account for the oscillations during the stick-slip transition observed in a specific system. To account for the spring-like characteristic and hysterisis effects in the sliding regime of a particular experiment, Canudas de Wit et al. [13] suggested a friction model which treats the friction interface as the contact between bristles. A similar model was later verified by Liang and Feeny [52] in their experiments. Marui and Kato [57] noticed oscillations in the friction signal during the stick-slip transition similar to the ones noticed in boundary lubricated system by Polycarpou and Soom [67], and Ko and Brockly [10] in their unidirectional pin-on-disk apparatus. These oscillations were a superposition of high-frequency signal over a low frequency periodic stick-slip signal. To explain this, they introduced a tangential compliance at the contact in a Coulomb friction model. The effect of tangential compliance in case of “pure” sliding motion is to induce an elastic deformation prior to the actual sliding motion. This elastic characteristic is termed “spring-like” sticking behavior by Canudas de Wit et al. [13]. The displacement induced by this elastic deformation is typically small (a few microns) and was referred to as “preliminary displacement”, “micro displacement”, “presliding displacement”, or “micro stick” in the literature (Oden and Martins [63], Canudas de Wit [13], Harnoy et al. [39] and Liang and Feeny [52]). The elastic deformation could belong to either of the contact elements. We have noticed and identified this compliant behavior of the contact and use a slighlty modified model of the same. The contact elements in our experiment were rubber and steel and hence we prefer to attribute the deformation to the former. 1.3.2 Beams subjected to frictional forces In our experiment we consider a cantilever beam subjected to frictional excitation at the free end. In this context we would like to review research done with beams subjected to frictional forces. There is a large catalogue of “vibrating beam” problems available in the literature (as Cusumano [16] states that “sometimes it seems as if two or more such problems appear in each issue of, say, The Journal of Sound and Vibration alone”). We here present a few of those that deal with friction. Bindemann and Ferri [6] have studied the damping characteristics of built-up structures damped predominantly by dry friction. In their studes they considered cantilever free beams subjected to transverse frictional interface. Dowell and Schwartz [21, 22, 23] carried out theoretical and experimental studies of the forced vibration response of a cantilevered beam with Coulomb damping non- linearity. They emphasized on the types of nonlinear behavior and system response spectral characteristics. They used a single mode Galerkin’s approach in the reduc- tion of the PDE to ODE problem and analyzed the ODE using harmonic-balance method. Yagasaki [92, 93] made theoretical and experimental studies of nonlinear vibrations of a straight beam clamped at both ends and forced with two frequencies near the first mode frequency. They used a single mode Galerkin’s approximation in analyzing the beam motion. Finite element approximations of the continuous models can also be developed which feature consistently-derived frictional damping and stiffness matrices. These finite element methods, together with numerical schemes for solving associated sys- tems of nonlinear ordinary differential equations, are capable of modeling stick-slip motion, dynamic sliding, friction damping, and related phenomena in a significant range of practical problems. 1.3.3 Phase-space reconstruction and dimension estimates As stated earlier, we would like to exploit the varied dynamics due to the friction excitation in our experiments. This brings us the issue of classifying the different 10 dynamics in a system into periodic, quasi-periodic, chaotic, etc. Periodic dynamics are easily identifiable. This is not true with quasi-periodic, chaotic and transically choatic dynamics. We use phase-space reconstruction as a foundation for characterizing these dynamics and use the dimensional estimates in deciding the number of POMs that are to be used for model reduction. In this context we here present the recent literature in the field of phase space reconstruction to identify a system using a time series. While carrying out experiments it would be difficult to measure the system vari- ables and more so to get all the system variables. This makes the understanding of the qualitative dynamics of the system difficult. By qualitative system information we mean the phase portrait information, the dimensions, the exponents of the system, etc. In such cases there are methods available in the literature to obtain information of the active states from a single observable. This procedure is called phase space reconstruction. The idea was first suggested by Packard [64]. Later Takens [87] suggested another method of phase portrait reconstruction which is known as method of delays. The work of Takens does not take into account the problems associated with measurement. Broomhead and King [11] developed another method of phase portrait reconstruction based on Takens’ proof and on ideas from the generalized theory of information known as singular systems analysis (SSA). The observation that the dynamics of a system with many degrees of freedom can be investigated using time series of a single scalar observable has broadened the class of experiments in which complex behavior can be interpreted as manifestations of strange attractors. Takens [87] proved that reconstructions are generically topo- logical embeddings of the original dynamics. Fraser [35] optimized the delay-vector phase-space reconstruction with a delay time that satisfies a minimum redundancy criterion, and found a reconstruction better than the one obtained using SSA. The superiority is due to the farmer’s notion of general independence as opposed to the latter’s foundation on the notion of linear independence. Pains and Dvorak [65] in 11 their study confirmed the pre-existing doubts about the reliability of the singular- value decomposition method. They concluded that singular-value decomposition, by nature a linear method, can bring distorted and misleading results when nonlinear structures are studied. So care should be taken in usage of this method in the systems with friction, as they could behave in a nonlinear fashion due to the nonlinear nature of the friction force, though the system inherently could be otherwise linear. Read [73] in his reconstruction attempts made use of multi-variate singular-systems analysis (M-SSA) and characterized a number of different kinds of flow and also de- termined their correlation dimension. In comparision with the more conventional . method of delays using a single variable, they found that the M-SSA offered sig- nificantly improved reconstructions in terms of signal to noise ratio and structural uniformity of the attractor. The method of delays approach suggested by Takens [87] and its variations is the most common method for reconstructing the phase space. Takens’ embedding theo- rem was given for smooth systems. The embedding theorem can be loosely stated as “if basic hypotheses are satisfied, applying the method of delays to an observable pro- duces a trajectory in the reconstructed or pseudo phase space which is diffeomorphic to the trajectory in the reconstructed or pseudo phase space which is diffeomorphic to the trajectory in the real phase space”. The friction—induced vibrations are generally nonsmooth and the presence of the discontinuity could have a profound effect on the analysis of data. Feeny [28] and Popp and Stelter [69] in their attempts at reconstruction using the method of delays, have shown the manifestation of the discontinuity during the analysis of the data. The effect of this problem in a reconstruction is that the consequent analyses for cal- culating the fractal dimensions, Lyapunov exponents and other nonlinear predictions become unavailable. This problem may not be suspected when the reconstruction is attempted from an experimental data where these discontinuties are not clearly identifiable. 12 The question would then be how to identify these problems? Knowing that the problem exists, is there a way of relaxing the smoothness requirements to continue with post-analyses? Some answers to these questions have been ventured for low-order systems (Feeny and Liang [31]). We try to overcome this problem using variables that are least affected by the discontinuity at the source. 1.3.4 Review of current research in POD The POD idea can be traced back to independent investigations by Kosarnbi [49], Loeve [54], Karhunen [46] and Pougachev [70]. Proper orthogonal decomposition is primarily a statistical formulation widely used in pattern recognition and image- processing communities. It is a procedure for extracting a basis for modal decom- position from an ensemble of signals. Its power lies in the mathematical properties that suggest that it is the preferred basis in many circumstances. The attractiveness of the POD lies in the fact that it is a linear procedure. The mathematical theory behind it is the spectral theory of compact, self-adjoint operators. This robustness makes it “a safe haven in the intimidating world of nonlinearity; although this may not do the physical violence of linearization methods” [5], the linear nature of the POD is the source of its limitiations, as will emerge from what follows. However it should be made clear that the POD makes no assumptions about the linearity of the problem to which it is applied. In this respect it is as blind as Fourier analysis [5]. The basic service of this method is to quantify the coherence in an ensemble of data. In the mechanics community, this technique was first exploited by Lumley [55], in 60’s in understanding coherent structures of turbulent flows. Until recently very few dynamicists have evinced interest in exploiting this tool for understanding the spatial distribution of energy in a dynamic system. In the recent past the POD theorem has been applied to estimate number of active states in chaotic attractors (Cusumano and Bai [17, 18]), model distributed systems (FitzSimons and Rui [34]), understand snap-through oscillations of buckled plates (Murphy [59]), investigate fluid structure interaction problems (Sipcic et al. [85]), etc., with different degrees of success. They 13 all have used the technique in quantifying the spatial coherence of the dynamics from observed information at different locations. The foundation of this theory lies on linear concepts like coherence and this could probably be the reason for the nonlinear dynamicists hesitation in using this theory for systems showing nonlinear phenomena. 1.4 Proposed Research The hypothesis for this research is “the proper orthogonal modes obtained from the chaotic dynamics of a higher order system subjected to frictional excitation broadly represent the system dynamics and hence can be used in building a reduced-order model for the system”. In the proposed research we conduct numerical and physical experiments and obtain time-series of displacements at different locations on the system when the system behavior is chaotic. We process this information using POD theory and obtain the dominant modes. We use these modes. in building a reduced order model for the system and validate the same using different qualitative and quantitative comparison techniques. We considered systems subjected to frictional excitation as a feasiblity study for squeak and squeal problems, for which our aim would be to identify the active modes and degrees of freedom that contribute to the system dynamic response. The POD provides modes which can be used in modal reduction in place of linear natural modes which in many situations are unknown. 1 .5 Contributions The chief contribution of this thesis is the exposure of an effective tool for reducing the order of higher order systems subjected to friction forcing. The POD method is simple to implement, is computationally cheap and requires a minimal learning effort in its application, in comparison with usage of traditional-modal-analysis software. In some cases, the number of POMs needed to capture the system dynamics is less than the required number of natural modes. For the systems that we considered, the 14 comparison with a “truth set” indicates that models built of POMS predict better dynamics than models built of an equal number of linear natural modes. The resulting reduced model can be used to predict structural resonant frequen- cies, passive damping levels and thus stratagically enhance the robustness of active control methods for flexible structures. This tool can act as a bridge between high and low order systems and thus transfer the vast amount of research being done in the latter systems to the former. By using POMS as a linear combination of the linear natural modes we throw some light on the degree of modal coupling involved in the dynamics of the system. In this process we also show the redundancy of the information that would have been used had linear natural modes been used. The auxiliary contributions include the experimental characterization of friction between rubber and steel and a validated method to measure displacements in the presence of drift using strain gages. In this process we also reinforce the feasibility of a method for numerically simulating high-order nonsmooth systems, modal projections and apply tools of nonlinear dynamical systems to model validation. 1.6 Thesis Organization In Chapter 2 we discuss the POD theory visa-vis simple linear systems and give an interpretation to proper orthogonal modes (POMS) with reference to linear natural modes. Also we make an attempt to relate the POMS to nonlinear natural modes of a nonlinear system. A feasibility study is described in Chapter 3. In this chapter we consider a numer- ical experiment. The experiment consists of a chain of masses conneected by springs and dampers subjected to friction forcing from a continuously moving belt at one end. The mathematical model for this system is known fully and can be numerically solved without making major assumptions. We carry some bifurcation studies on the full model and identify chaotic regimes using Lyapunov exponents, false nearest 15 neighbors and other phase space reconstruction techniques. We use the chaotic dis- placement information of the individual masses to form a spatial correlation matrix. We determine the POMS and POVs from this matrix and then use the dominant modes in the original mathematical model to obtain a reduced model. We verify this model using different qualitative and quantitative techniques. In Chapter 4 and 5 we extend the method to a physical experiment. These chap- ters include a description and characterization of the physics of the experiment. We present here the different measurement methods and determine the friction char- acteristic. In Chapter 5 we characterize the nonlinear dynamical response of the experiment by presenting frequency and amplitude sweeps and delay maps. We again ' identify the chaotic regime using Lyapunov exponents and F NN. We end the chapter with determining POMs and POVs. In Chapter 6 we describe the mathematical model for the system and discuss the effects of different parameters on the system. We later build and validate a reduced model for the system. The nonlinear dynamic response is characterized to further validate the model. In Chapter 7 we conclude with a summary of the research conducted, lessons learned and directions for future work. Chapter 2 Proper Orthogonal Decomposition (POD) In this chapter, we do a brief review of the ideas behind proper orthogonal decom- position and summarize the relation between proper orthogonal modes as applied in oscillating systems to normal modes of oscillation. 2.1 Definitions and Theory Ensemble is a commonly used term when dealing with statistical phenomena. We present a formal definition for this term in the context of our experiments. Consider a quantity which can be given a numerical value, e.g., displacement of the beam tip u at a particular time t which is being measured in an experiment. In statistical studies we imagine that the same experiment is done repeatedly and each time the value of the same quantity is recorded. At the end of a large number of repetitions of the experiment, we have a collection of the data measured under superficially identical conditions. This collection of data is referred to as ensemble and the members are termed as realizations. The quantity being measured could be the coordinates of a point in a n-dimensional space and the repetition can be done at a uniform sampling rate. Thus the ensemble could consist of displacement information at various locations 16 17 in a dynamical system measured at uniform sampling rate. When there is an ensemble of signals, the proper orthogonal decomposition theo- rem [54] helps in finding a deterministic function which in some average sense has a structure'typical of the members of the ensemble. The functions when found would be the “most similar” to the members of the system represented by the ensemble on the average. Mathematically the notion of “most similar” implies that the deviation of the member functions of the ensemble from the found function is minimum. The advantage of POD to other basis functions lies in the claim that this basis has the optimal distribution in some sense. When talking about fluid-flow systems, the func- tions would be the optimal functions in the sense that they capture more power than any other set of basis functions (Lumley [55]). The notion of “most similar” in mathematics corresponds to seeking a function o such that (Kn, ¢)|’) _ (l(u,¢)l’) max (m ' , (M) ' ‘2'” where () indicates the statistical mean of the quantity under consideration and (2:, y) indicates the inner product of a: and y. That is we are searching for a member «:5 among the dis, which maximizes the normalized inner product with the field u. Let us refer to the maximum value in equation (2.1) by A. This is a classical problem in the calculus of variations. Using kernel theorem and calculus of variations we can show that a necessary condition [55], for equation (2.1) to hold is, 43 Should be an eigenfunction satisfying the eigen-relation Law-(WW. = W). (2.2) The term (u(t)u"(r)) is the two-point correlation tensor, and we define R(t,r) to be this average. Spectral theory [75, 14] guarantees the existence of the maximum in (2.1) and it corresponds to the largest eigenvalue A1 of (2.2). Hilbert-Schmidt theory [75] assures us there is not one, but countably many (the number equals the cardinality of R) of solutions of (2.2), as long as Q is bounded. The non-negative 18 definiteness of R(t, r) assures that A,- 2 0. We also have a diagonal decomposition: RU. T) = Z Ak¢k(t)¢Z(T)- (2-3) I: _ Also almost every member of the ensemble may be reproduced by a modal decompo- sition in the eigenfunctions: ‘U(t) = Zak¢k(t) (2.4) I: The optimality of the basis can be seen by looking at the decomposition of a signal u(:r, t) with respect to an orthonormal basis 1b,: u(r,t) = Z b;(t)¢g($). (2.5) If the w.- have been nondimensionalized and normalized to give («pm/5) = 6.5, then the coefficients b.- carry the dimension of the quantity u. If u(.'c,t) is a displacement measurement and the signal energy is defined as In uu‘dz then, the average signal energy over the experiment is given by [a (uu')d:r = z: b,b‘:. (2.6) Hence, we may say that (bib?) represents the average signal energy in the i-th mode. The optimality of the POD has been summarized by Berkooz et al. [5] in the form of the following proposition. Proposition: Let u(:r,t) be an ensemble member square integrable on the domain It for almost every t and {¢.~,A;} be the POD orthonormal basis set with associated eigenvalues. Let ' u(:c, t) -.= Z a;(t)¢,-(:c) (2.7) be the decomposition with respect to this basis, where equality is almost everywhere. Let w.- be an arbitrary orthonormal set such that "(13, t) = Z b.-(t)¢,-(:c). (2-8) 19 Then the following hold: 1. (a;(t)a;(t)) = 6.3-Ag, i.e. the POD coefi‘icents are uncorrelated. 2. For every n we have )2? (a;(t)a‘;(t)) = E? A; Z 2? (b.-(t)b;'(t)). The claim that POD gives the most optimal basis is based on the above proposi- tion. It implies that among all linear decompositions, POD is the most efficient, in the sense that, for a given number of modes the projection on the subspace used for modeling will contain the most energy possible in an average sense. In addition the time series of the coefficients a.-(t) are uncorrelated. 2.2 The Method To use this method one obtains the covariance of the information, and finds the proper functions and corresponding proper values of the matrix. The proper values indicate the fraction of energy of the system (represented by the ensemble) captured in the corresponding mode. Thus one can determine the modes that dominate the system from the energy sense. These dominant modes can be picked as empirical and be used in building a reduced model for the complex system. To understand how this technique works, let us consider a simple dynamical sys- tem consisting of masses connected by springs and dampers. Let us represent the displacements of the N masses by X(t) = [X1(t),X2(t),....,XN(t)]T. The system equation can then be written as Mi + CX + Kx = F (2.9) where M, C and K are the system’s mass, damping and stiffness matrices. From a real or numerical experiment, one can measure X(t) at different time steps and thus get a T xN matrix 3, where each of the T rows represent the displacements of the masses at a particular time. Note that in matrix 3, we have arranged the dis- placements such that each column represents the displacements of a mass at different 20 instants of time. Using this matrix S we can construct a covariance matrix R = STE. The eigen- vectors of this matrix R (which is an N x N non-negative definite matrix) can easily be determined. As the eigenvalues indicate the energy contained in the correspond- ing mode, we can decide on the modes that dominate the system from the energy perspective using these eigenvalues. 2.3 Relation to Normal Modes — Unforced and Undamped Oscillations Consider an undamped system undergoing oscillations due to initial conditions alone. Let us assume the mass matrix of the system can be reduced to an identity matrix. Let us suppose the system is orbiting with the oscillations as a combination of only two natural modes say V,- and V,-. i.e. X(t) = A; sin(w,°t + ¢;)V.' + 141' sin(wjt + ¢j)Vj (2.10) Let us collect the displacements from the experiment at different time steps to get a matrix of E as shown below. a = av?" + ejvf (2.11) where ' A: sin(w1t1 + 451) l e; = A; sm(w:t2 + 451) (2.12) . Azsin(wjtT + (51) . for l = i, j. After measuring these displacements we can form the covariance matrix, which would be 21 II A T+ V Je-T)(e,-VT + eJVT) = V.-e TegVT+VgeT eJ-VT+VJJ- eJ Te;VT+VJ- -e- TeJ-VT =e Te5VVT+eT eJ-(VVT+VJ- VT)+eT eJ-VJ- VT (2.13) Now if e;, eJ- are sinusoidal vectors with different frequencies, the expected value of e,TeJ- would go to zero, as T goes to 00. Thus, R =e TegV-VT+eT eJ-VJ- VT = aV v?‘ + (N v? (2°14) We can easily verify that V.- and VJ- are the eigenvectors for the above covariance matrix by post-multiplying both sides of the above equation by V.- and VJ- respec- tively. There are only two nonzero eigenvalues far the matrix R in this case and they are in the ratio (Ai/Aj)2. ‘ 2.4 Relation to Normal Modes in General Systems In the case of systems with damping and external forcing, the relation between normal modes and POMS is less obvious and difficult to ascertain analytically. Feeny and Kappagantu [27] used analytical and numerical methods to relate normal modes and POMS in case of simple systems with general damping and forcing. They considered symmetric systems which are typical of structures. In their attempts with damped systems they considered modal damping factors low enough that several cycles can be observed. For proportionally damped systems with synchronous modes [58] they observed the POMS tend to eigenvectors of the system with increasing number of samples. The error decreased with damping. In a simple nonlinear system the POMS were found to lie along the principal axes of inertia of the nonlinear normal modes. In their numerical experiments they assumed the normal modes to be synchronous and the dynamics on the invariant manifold, 22 when projected onto the coordinate space, do not show any hysteretic behavior in the coordinate space during motions. With these assumptions the ensemble accumulates on a nonlinear curve in the coordinate space. They have verified that the POMs represent the principal axes of inertia of the ensemble of data. The POVs represent the mean-squared distances of the ensemble along these axes. They ascertained the fact that though the POVS do not reflect the true energy or power in the system, they represent the mean-squared distance of the ensemble which in mechanical systems is generally associated with energy or power. For the case of fluid flows where velocity measurements form the ensemble, these mean-squared distances represent the kinetic energy consistent with Lumley’s results. We wish to analytically delve into the application of POD to systems with general damping and also nonlinear behavior at a later date. 2.5 Recent Engineering Applications of POD Though Lumley [55] was the first person to report the usage of proper orthogonal decomposition theorem in the study of turbulent flows to understand the coherent structures, POD is better established as a statistical tool for signal analysis and data compression. Almost 95% of the literature found in the application of POD in mechanics is in the study of turbulent flows. Of late a few structural dynamicists have started extending the POD-Galerkin’s projection of a PDE such as the Navier- Stokes equation to understand the behavior of the structure from the resulting finite- dimensional system. Fitzsimons [34] used these ideas in reducing linear and nonlinear smooth systems (consisting of strings). They excited the systems with a random function to generate the ensemble and applied the POD and selected the dominant modes for using in Galerkin’s method to build a reduced order model. Cusumano [17, 18] has applied the POD theorem in studying the spatial complexity of the motions of a 10 degree-of-freedom impact oscillator. He estimated the number of excited degrees of freedom by finding the number of modes needed to exceed 99% 23 of the signal power and gave the distribution of energy in dominant modes. He compared this number with the predictions of fractal dimension theory and found that the dimension of the attractor in the phase space matches with the number of dominant modes. In our preliminary exploratory studies of the POD to friction excited systems we considered a similar system as that of Cusumano but with friction excitation, which we present in detail in the next chapter. Chapter 3 POD of Friction Induced Vibrations: A Numerical Experiment In the end of last chapter we have stated that proper orthogonal decomposition was originally used in understanding coherent structures in turbulent flows, and recently was applied to impacting systems to compare with the phase-space reconstruction results. Here we hypothesize that POD techniques can be applied to the dynamics of systems subjected to frictional forces (which are nonsmooth in general). The applica— tion would be the building of a reduced order model for the continuous or higher-order system. The main question that arises would be how to pick the dominant modes and how to validate the reduced model. The answer for the first question could involve usage of chaotic data that is typical of a friction system. Validation of the reduced model can be done by comparing the time series data, bifurcation results, etc., of the full and reduced systems. To understand the application and validity of applying the POD theorem to the frictional system, we first study a numerical system, which is described in the following section. 24 25 3.1 System Description and Modeling Issues To verify our hypothesis, we conducted some numerical experiments. We have con- sidered a system with a chain of N masses connected by springs and dampers. One end of the chain is anchored. The other end mass is lying over a moving belt with only the friction force acting between the belt and the mass. The configuration is illustrated in Figure 3.1. All the masses, springs and dampers are assumed identical. Proportional viscous damping is assumed. The system equation can be written as l > '—NKA—A§\A *— CJ— m 73— ’“ 73- Figure 3.1: Configuration of spring mass damper system under frictional excitation Mi+cKX+KX=F (3.1) where c is the damping coefficient, M, K are the mass and stiffness matrices and are given by 26 'mOO ooo OmO 000 00m 000 M=EEE'-.EEE (3-2) 000---m00 000 OmO _000 00m ’2k—k0 0001 —k2k-k 000 0—k2k 000 K=:EE'-.EEE (3-3) 000---2k—k0 000---—k2k—-k _000-~-0—kkd where m is the mass of each block and k is the stiffness of each spring. The friction law is an approximation of the Stribeck friction law and is similar to the ones used by Sugimoto [86]. The friction force F R is a function of relative velocity V, and normal load N and can be stated as FR = (11K + (#5 - lime-"MUN V, > 0 (3-4) —p51v < FR < +PsN Vr = o (3.5) FR = -(#K + (#5 - #K)€'°'V")N V, < 0 (3-6) where pg and p5 are the kinetic and static coefficients of friction. Noting that the friction force is acting only on the last mass, (the driving mass), we see that the F vector in the (3.1) is all zeros except the last row. The force F is a column vector given by 27 F = g F,(a,v,,N) (3.7) 0 1 b d The above set of equations are only coupled through the matrix K and can be solved directly by standard numerical ODE solvers. Even when the coupling is through both M and K matrices one can reduce the (3.1) into a non-dimensional form. We assume for simplicity m and k to be unity and the damping to be propor- tional and viscous. Therefore, 5': + ch + xx = F (3.8) Using the relations Y;- = X,- and Y,“ = X;, we can rewrite the above equations in a state—space form: Y = AY + B f(t) (3.9) where the matrices A and B are A = [ i ix 1 (3.10) B = [ g ] (3.11) When the driven mass is slipping over the belt, f (t) in the above equation represents the friction force which is known as a function of relative velocity and normal force. These equations can be solved by any ODE solvers. We have used IMSL ODE solver DI VPRK. For clarity let us denote A and B by A“ and B“. The subscript “s1” stands for slip. 28 However when the driven mass is stuck to belt, f represents the static friction force. We know only the bounds (—psN,+psN) but not the exact value. So the above state space equation cannot be solved in that form directly. Using the ideas of Ferri et al. [32] and P0pp et al. [68] we know that, when the driven mass is stuck to the belt, the acceleration and the velocity of the driven mass would correspond to those of the belt: 2N a(t) = YzN = z A2N,iy£ + BszU), (3.12) i=1 ”(0 = YN = Y2N (3-13) where, a(t) and v(t) are the acceleration and velocity of the belt respectively. The static force f (t) varies with time, but more specifically is dependent on the time]- varying displacements of the system. Denoting the first (2N - 1) elements of the vector Y by Z, we can write the following relation Y=[(I)]Z+[v?t)] (3.14) °’ Y=[;]z.[ag)] (3...) where I and 0 indicate identity and zero matrices of matching sizes. Substituting these relations in (3.9) and writing f (t) in terms of a(t) as can be obtained from the (3.12), we have [31212.1411:1“1t]1+;7w>-F:.Bc11:12+1,1,1] rearranging, I - 0 1 I l 0 l [o]z.[am]#112304012....-EBC)[W,]+§;B.(.) The first 2N —1 equations in the above matrix equation corresponds to the state- space equations for the state Z. By careful observation we can see that the last 29 equation turns out to be an identity relation. One can therefore write the state equations corresponding to the stuck mode as: 2 = AuZ + B.,v(t) (3.18) where A“ is given by the 2N — 1 x 2N — 1 block of the matrix (A — fiBC) and B“ is given by the first 2N -1 elements of the 2N th column of the matrix (A — -B:—NBC). Here subscripts “st” indicate stuck mode. As know the velocity of the belt, we can now integrate (3.18) for Z and evaluate the corresponding Y using the algebraic relation (3.14). The implementation details are as follows. From the friction characteristic, the driving point can be either in the stick mode or in the slip mode. The system is started with some initial conditions and the belt moving at some speed V. From the initial conditions and the belt speed, the system mode whether in slip or stick is determined. The system is in the slip mode when the driving mass velocity is different from that of the belt. If the two speeds are same, then the restraining force 5 on the driving mass (due to the springs and dampers) should be considered. From equation 3.12 the restraining force has the form 5 = 512710) — :1: AgN,.-}’.-. (3.19) If the absolute value of this force 5' is within the static friction force p5, the system is in stick mode. If it is any greater than the static friction force, the system will slip. Let us assume the system is in slip mode. When in the slip mode the instantaneous kinetic friction force is calculated based on the relative velocity between the driving mass and the belt and this is used as an external excitation force. The integration is carried out in small time intervals At, until the crossover of the beltspeed by the speed of end mass (i.e. X N crosses V). When there is a crossover, the program steps back and uses a reduced time step to calculate the state. This is continued till the end mass speed nearly equals the belt speed. At this stage, the restraining force 5 on 30 the end mass is computed. If the restraining force satisfies the static force constraint, i.e. IS I < psN, sticking occurs. If IS I Z psN the system continues to slip and we continue to integrate the slip state equations. When the system gets into the stuck mode, we integrate the state~equations cor- responding to the stuck mode in small time intervals At and at each instant the restraining force 5 is calculated. This restraining force is compared with the its, to check for the slipping criterion. When IS I crosses p5, the program steps back and re-integrates with a reduced time step, until ISI equals ps. The system is then ready to slip. The slip module takes over until it sticks again. These steps are repeated for the specified time duration. For a given initial condition, the system is first checked for the mode in which it is to start off and the integration is continued using the same. The kinetic friction force is to be updated at every time step based on the relative velocities. In the course of time the exact state of transition from one mode to the other is calculated and is given as initial condition for the subsequent mode. To calculate the periodicity we run the system for a specified number of time steps at the end of which we check for any periodicity based on the slip criterion (the absolute difference between the states at successive slips is compared with the last slip recorded). In our simulations we assumed constant belt speed. We have thrown out the first 1000 seconds of data before we check for periodicity. In the absence of periodicity, we continue the simulations for another 5000 seconds before the next check. In the absence of periodicity we continued this process until we reached 100, 000 seconds. 31 12 I I 1' I I I I I I .a O slip displacement a) O 0.3 b §02 -O.7 . .03 4.5 6 5 8.5 5 5 7 1 5 11 x10 x10 x10 Figure 3.2: Bifurcation diagram of the full model and a few different possible orbits. The belt speeds for these orbits are 0.20, 0.22 and 0.25 3.2 Higher Order Dynamics, Chaos and its Verifi- cation Multi-degree-of-freedom systems are known be rich in dynamics. To study this we conducted a sweep of the belt-speed from 0.5 to 0.0 in steps of 0.001. We found that the system gives rise to different kinds of orbits: period 1, period 2, quasi-periodic, chaotic, etc. These are reflected in the bifurcation diagram, where we have the belt- speed on the x-axis and the points of slip on the y-axis. Thus at a given belt-speed, n points indicate the period is 1:. When n is large, i.e., when there is a smear of points it indicates the dynamics are either non-periodic or the system is still in the transient stage. Systems involving frictional forces are known to behave chaotically [31, 42]. This does not obviate the need to check for veracity of the chaotic nature of the seemingly 32 chaotic orbit. It is possible that, for some initial conditions, a nonlinear system can go to a steady-state periodic orbit after a long time, and the observer may collect the data from the transients. Chaos, as a property of orbits z(t), manifests itself as complex time traces, continuous, broadband Fourier spectra, nonperiodic bounded motion and exponential sensitivity to small changes in the orbit [1]. We now have a need to detect and quantify chaos. One way to check for the chaotic nature of the orbit is to look at the spectrum of Lyapunov exponents [91]. Lyapunov exponents are the average exponential rates of divergence or convergence of nearby orbits in phase space. Since nearby orbits correspond to nearly identical states, exponential orbital divergence means that systems whose initial differences we may not be able to visually resolve may soon behave quite differently, i.e. predictive ability is rapidly lost. Any system containing at least one positive Lyapunov exponent is defined to be chaotic, with the magnitude of the exponent reflecting the time scale on which system dynamics become unpredictable. This brings us into the problems faced in reconstruction techniques. Typically most of the theorems that are used in the art of reconstruction from a time series data needs a degree of smoothness. This smoothness is essential in allowing the demonstration that all the invariants of the motion as seen in the reconstructed time delay space are the same as if they were evaluated in the original space. This means we can work in the time delay space and learn as much as we could about the system were we able to make our calculations directly in the true space and other variables. Because of the inherent nature of the friction characteristic we used we have a non-smoothness in the system and this is immediately reflected in the displace- ment / velocity of the driven mass. Some methods have been proposed in the literature to take care of this non-smoothness, by modifying the “now” stande methods [29]. The usage of another state variable measurement such as phase etc., in conjunction with the available time series, is the approach. In our numerical experiment though there is non-smoothness in the driven mass 33 displacement and velocity, we have in hand the time series of the displacement and velocities of even other masses. By taking the velocity of the mass farthest from the driven mass we have a time series that is minimally affected by the the non-smoothness of the forcing. We use this time series in determining the Lyapunov exponents. Let us denote this time series by 3(n). In order to calculate the Lyapunov exponents from a scalar time series one needs to calculate some parameters from the time series. If we keep in mind that we have a time series of only one variable, we realize that we need an equivalent system phase space. We have to unfold the projection back to a multivariate state space that is representative of the original system. This can be obtained by using time delay reconstruction, which involves determining the necessary number of coordinates to unfold the attractor from self overlapping caused by projection. The dimension that unfolds the attractor enough so that we have no overlaps is called the embedding dimension. One has to determine this embedding dimension so that the attractor can be reconstructed before we can attempt to find the Lyapunov exponents. To determine the embedding dimension using time delay reconstruction, we need to determine the time delay. If the time delay is too short, the coordinates 3(a) and s(n + T) will not be independent enough: not enough time will have evolved for the system to have explored enough of its state space to produce, in a practical numerical sense, new information about that state space. On the other hand since chaotic systems are intrinsically unstable, if the time lag is too large, any connection between the measuements s(n) and s(n + T) is numerically tantamount to being random with respect to each other. So we need a time delay which is large enough that s(n) and s(n+T) are rather independent but not so large that they are completely independent in a statistical sense. Towards this goal we make use of chaos as an information source [36, 35]. Fundamental to the idea of information among measurements is Shannon’s idea of mutual informtion [38] between two measurements. This is a measure of the amount of information learned by one measurement from the other. We therefore determine the “mutual information” (auto-correlation is sometimes used instead of 34 mutual information, but it is a linear concept and here we are dealing with nonlinear system). 3.2.1 Mutual information The mutual information between two measurements A = {0;} and B = {bj} is defined as the amount learned by the measurement a,- by the measurement bj and is given by the relation I =1og, [ P‘B(“"b") (3.20) PA(a.-)P3(b,-) where PA3(a, b) is the joint probability density for the measurements A and B re- sulting in values a and b. PA(a) and P3(b) are the individual probability densities for the measurements of A and B. If the measurement of a value from A result- ing in a.- is independent of the measurement of a value from B resulting in ()5, then PAB(a, b) = PA(a)PB(b) and the amount learned is zero, as it should be. The average of this statistic, called the average mutual information between measurements A and B, over all measurements is PAB(aiv bi) PA(a,-)P3(b,-) IAB = Z PAB(a,-, bj) logz [ (3.21) M Now take as the set of measurements A the values of the observable s(n) and for the B measurements, the values of s(n + T). Then the average mutual information between these two measurements, namely the amount learned by measurements of 3(a) through measurements of s(n + T) is P(8(n),3(n + T)) P(8(n))P(3(n + T)) I(T) = Z P(s(n), s(n + T)) logz (3.22) By general arguments [38] I (T) _>_ 0 (as PAB(a.-, b,) is at least PA(a,-)P3(bj) when both the signals a and b are independent of each other). I (T = 0) is directly related to the Kolmogorov-Sinai entropy [35]. When T becomes large, the chaotic behavior 35 of the signal makes the measurements s(n) and s(n + T) become independent in a practical sense, and I (T) will tend to zero. It was the suggestion of Fraser [35] that one use the function I (T) as a kind of nonlinear correlation function to determine when the values of s(n) and s(n + T) are independent enough of each other to be useful as coordinates in a time delay vector but not so independent as to have no connection with each other at all. The actual prescription suggested is to take the T where the first minimum of the average mutual information I (T) occurs as that value to use in time delay reconstruction of phase space. Now that we have a formulae for average mutual information, let us look into the implementation details. The main problem would be determining the probability densities. The probability P(s(n)) can be obtained by taking the ratio of number of points in a 1-D line with the coordinate s(n) to the total number of samples. Similarly joint probability density P(s(n), s(n + T)) is obtained by taking the ratio of number of points in a 2—D plane with coordinates (s(n), s(n + T)) to the square of the total number of samples. However as we have only finitely many samples with finite precision we can divide the 1-D line and 2-D plane into finite grids and keep a count of the number of points falling in the individual cells. Then the ratio of samples in a cell to the total number of samples corresponds to the probabilty of the cell value. In this case, the probability would be a function of the grid size. So one has to decide on the grid size judiciously by playing with the grid size. As mentioned earlier we have taken the time series of velocity of the first mass (which is the farthest from the driven mass) as the sample signal s(n) and applied the above algorithm of finding mutual information. Plots of T versus I (T) are shown in Figure 3.3 for different grid sizes. Clearly we see that the first minimum occurs around T = 11 and hence this would be taken as the time delay in the future sections. 36 0.3 1 0.25 ' 0.2 I of calls - 5. 10. 20. 40. BO Figure 3.3: Average mutual information of the lO-SMD system, using velocity time series of the first mass. The numbers of cells used in estimating probability density are 5, 10, 20, 40 and 80 and the corresponding plots go towards the axes. The first minimum occurs near T = 11. 3.2.2 Embedding dimension We already saw the need to determine the embedding dimension. From the point of view of the mathematics of the embedding process it does not matter whether one uses the minimum embedding dimension (IE or any d 2 d3, since once the attractor is unfolded, it remains unfolded for higher dimensions. For an experimentalist the story is quite different. Working in any dimension larger than the minimum required by the data leads to excessive computation when investigating any subsequent question (Lyapunov exponents, prediction, etc.) one wishes to ask. It also enhances the problem of contamination by round-off or instrumental error since this noise will populate and dominate the additional (1 — d5 dimensions of the embedding space where dynamics are not active. The usual method of choosing the minimum dimension dB is to compute some 37 invariant on the attractor. By increasing the embedding dimension used for the computation one notes when the value of the invariant stops changing. Since these invariants are geometric properties of the attractor, they become independent of d for d 2 (13. The problem with this approach is that it is data intensive and is certainly subjective. Furthermore, the analysis does not indicate the penalty one pays for choosing too low an embedding dimension. Kennel [47] presented the idea of determining the embedding dimension based on false nearest neighbors. In the passage from dimension d to d +1 one can differentiate between points on the orbit y(n) that are “true” neighbors and points on the orbit y(n) which are “false” neighbors. A false neighbor is a point in the data set that is a neighbor solely because we are viewing the orbit in too small an embedding space (d < dE). When we achieve a large enough embedding space, all neighbors of every orbit point in the multivariate phase space will be true neighbors. We used Kennel’s idea in determining the embedding dimension. We follow the following steps in determining the false neighbors. Start with dimension d = 1. For each of the points y(n) determine the closest point yc(n), using the norm defined by the square of the Euclidean distance between ‘ them D(y(n)). This is given by d D(y(n)) = 201001 + iT) - W? + 2T))2 (3-23) i=0 Between each point and its closest neighbor, compute the value of D(y(n)) in higher dimensions. This involves just the addition of the term 61) = (gt/(;(n+kT)—y(n+kT))2 to the previous distance. Here I: is the latest dimension number. A natural criterion for catching embedding errors is to compare the increase in distance between the closeset neighbors y(n) and yc(n) with some threshold. i.e. use the criterion FED > 6101 (3.24) The usual practice in literature is to use c101 = 10 and we used the same. By this we are deciding a neighbor to be false when the proportional change in the distance with 38 an additional dimension is ten times the previous distance. This criterion alone is not sufficient to determine a proper embedding dimension. We have to make sure that the actual distances are also within some threshold. If either of the two are outside the thresholds, the neighbor that we found closest is not really a true neighbor. We have a false nearest neighbor. We keep a count of all these false nearest neighbors and store it across dimension d. Now increase d by 1 and continue the same steps to determine the false nearest neighbors. We continue this till the false nearest neighbors go to zero. When this happens we have the minimal embedding dimension d3. F 38688388 ‘ O r 0° 0M 16 .. l0 0 b an ~10 can 00 dimension 6 Figure 3.4: False-nearest-neighbors in the 10-SMD system, using velocity time series of the first mass. Note the first zero occurs near d = 6. In Figure 3.4 we give a plot of number of false nearest neighbors versus dimension d. Notice that at d = 6 we have zero false-nearest-numbers. Hence we take d5 = 6 in further calculations. 39 3.2.3 Lyapunov exponents We have already seen the meaning and importance of Lyapunov exponents in earlier sections. Here we discuss the method of finding them. Wolf [91] presented the first algorithm to estimate the non-negative Lyapunov exponents from an experimental time series. This method was rooted conceptually in an earlier method which could be applied only to analytically defined systems. It basically monitors the long-term growth rate of small volume elements in an attractor: in an n-dimensional state space, an infinitesimal n-sphere will become an n-ellipsoid due to locally deforming nature of the flow. The ith one-dimensional Lyapunov exponent is defined in terms of the length of the ellipsoidal principal axis p;(t) : P_i(t) Ps"(0) where the A,- are ordered from largest to smallest. Thus the Lyapunov exponents A,- = lim- tl—og2 (3.25) are related to the expanding or contracting nature of different directions in phase space. Since the orientation of the ellipsoid changes continuously as it evolves, the directions associated with a given exponent vary in a complicated way through the attractor. One cannot therefore speak of a well defined direction associated with a given exponent. Notice that the linear extent of the ellipsoid grows as 2"“, the area defined by the first two principal axes grows as 2(*1+’\’)‘, the volume defined by the first three principal axes grows as 2l"1+"’+"3)‘, and so on. This prOperty yields another definition of the spectrum of exponents: the sum of the first j exponents is defined by the long term exponential growth rate of a j-volume element. This alternate definition will provide the basis of our spectral technique for experimental data. Here we have discrete measurements of a single observable. Using the delay coor- dinates (obtained in previous sections) we do a phase-space reconstruction to get an attractor whose Lyapunov spectrum is identical to that of the original attrctor. We define the Lyapunov exponents by the phase space evolution of a sphere of states. At- 40 tempts to apply this definition numerically to equations of motion fail since computer limitations do not allow the initial sphere to be constructed sufficiently small. Assum- ing linear approximation holds at the smallest length scales defined by our data we can work in a reconstructed attractor, examining orbital divergence on length scales that are always as small as possible, using an approximate Gram-Schmidt reorthonor— malization (GSR) procedure in the reconstructed phase space as necessary. To estimate A1 we in effect monitor the long-term evolution of a single pair of nearby orbits. Our reconstructed attractor, though defined by a single trajectory, can provide points that may be considered to lie on different trajectories. We choose points whose temporal separation in the original time series is at least one mean orbital period, because a pair of points with a much smaller temporal separation is characterized by a zero Lyapunov exponent. Two data points may be considered to define the early state of the first principal axis so long as their spatial separation is small. When their separation becomes large we would like to perform GSR on the vector they define (simply normalization for this single vector), which would involve replacing the non-fiducial data point with a point closer to the fiducial point, in the same direction as the original vector. With finite amounts of data, we cannot hope to find a replacement point which falls exactly along a specified line segment in the reconstructed phase space, but we can look for a point that comes close. In eflect, through a simple replacement procedure that attempts to preserve orientation and minimize the size of replacement vectors, we have monitored the long-term behavior of a single principal axis vector. Each replacement vector may be evolved until a problem arises and so on. This leads us to an estimate of A1. The use of a finite amount of experimental data does not allow us to probe the desired infinitesimal length scales of an attractor. These scales are also inaccessible due to the presence of noise on finite length scales and sometimes because the chaos- producing structure of the attractor is of negligible spatial extent. Given the time series s(t), an m-dimensional phase portrait is reconstructed with 41 W Figure 3.5: A schematic representation of the evolution and replacement procedure used to estimate Lyapunov exponents from experimental data. The largest Lyapunov exponent is computed from the growth of length elements. When the length of the vector between two points becomes large, a new point is chosen near the reference trajectory, minimizing both the replacement length L and the orientation angle 0. delay coordinates, i.e., a point on the attractor is given by [s(t),s(t + T), ...,s(t + (d - 1)T)], where T is the time delay that we determined using mutual information and d is the embedding dimension that we determined in the previous section. We locate the nearest neighbor (in the Euclidean sense) to the inital point [s(to), s(to + T), ..., s(to+ (d — 1)T)] and denote the distance between these two points L(to). At a later time t1, the initial length will have evolved to length L’ (t1). The length element is propagated through the attractor for a time short enough so that only small scale attractro structure is likely to be examined. If the evolution time is too large we may see L’ shrink as the two trajectories which define it may pass through a folding region of the attractor. This .would lead to an underestimation of M. We now look for a new data point that satisfies two criteria reasonably well: its separation, L(t1), from the evolved fiducial point is small and replacement elements is small. If an adequate 42 replacement point cannot be found, we retain the points that were being used. This procedure is repeated until the fiducial trajectory has traversed the entire data file, at which point we estimate M LI( )3 I L (1:3) (3.26) =tMl-tok- where M is the total number of replacement steps. In this implementation we have held evolution time step (A = tk+1 — tk) constant. In the limit of no infinite amount of noise-free data this procedure always provides replacement vectors of inifinitesimal magnitude with orientation error and A1 is obtained by definition. This method has been used on the velocity time series of the first mass of the lO-SMD system which resulted in 0.2833 as the most positive Lyapunov exponent, thus verifying the chaotic nature of the data under consideration. 3.3 Finding POMS and System Reduction We collected the displacement data of different orbits and for each of them we de- termined the proper-orthogonal-modes as described in earlier sections and presented them in Figure 3.6. For each of the different orbits we see that the system shows a difl'erent set of POMS. Though for majority of the cases only two modes were found to dominate and these modes happen to be same in the limit, we do find different sets of dominant modes for different kinds of orbits. This gives rise to a problem of finding a single set of POMs that broadly represent the system. We conjecture that the POMs obtained from the chaotic orbit would represent the system in the broadest sense possible. This is because of the fact that a chaotic orbit visits many unstable periodic orbits and hence should contain more information about the system than a single periodic orbit. Also chaos has random like properties. 43 O 3 1 O 3 0 . 2 2° x x -0.3 , —o 3 -0 8 ‘0- 1 s 11 65 05 X10 X10 0 , o o 4 -0.5 —0.5 -0.6 0.5 ' 0:, 0.5 0.0097 0.1 _0.5 -'0.5 —O.5 ... .5 O O .5 O 0 Figure 3.6: Typical orbits of the lO-SMD system in the first row and corresponding dominant POMs are shown in the bottom row. The first two orbits are obtained at a belt speed of 0.25 with different I.C.’s and the third obtained at a belt speed of 0.20. 3.4 Modal Reduction 3.4.1 The mathematics From the chaotic orbit we collected the displacement data of all the N masses and calculated the spatial covariance matrix. The singular value decomposition of this Hermitian matrix gives us the eigenmodes which are the proper orthogonal modes of the system. The eigenvalues are put in a decending order. The dominant modes are then selected. We used these modes in transforming the system coordinates Y, to a new set of coordinates Q, using the relation Y = VQ + Y, where V is the block diagonal matrix, and both the blocks are of size N x L and are formed of the dominant modes (each column in the block represents a modal vector), and Y— is the mean of Y across time. By measuring Y from the points of static deflection of each mass, we can easily show that Y equals zero. Thus we can write Y = VQ. The usage 44 of this transformation in the system’s full model equation results in a reduced system equation. The dimension of the system equation in state space notation changes from 2N to 2L in this reduction. Here N is the number of DOF and L is the number of dominant modes. Equation (3.9) then leads to vi; = AVQ + B fR (3.27) Premultiplying this equation by (VTV)‘1VT results in c“) = (VTV)"VTAVQ + (VTV)“VT B f}; (3.28) Denoting (VTV)‘1VTAV by K.) and (VTV)_1VTB by F“ we can rewrite the above equation as Q = KuQ + FufR (3-29) These 2L equations can be integrated as before to obtain the values of Qgs at each time step during slipping. From Q;s one can get back the value of YES from the relation Y=VQ During stick mode the mathematics is similar to that of the full model. Again we make use of the velocity and acceleration constraints. . 2L v(t) = YN = Y2N = Z V2N,.°Q: (3-30) i=1 0 2L 0 a(t) = Y2N = Z V2N,£Qi ' (3.31) i=1 1 1 2” 2L t = —at -— A ,- Vg- - f() BzN () Burg; 2mg .JQJ l l = - EV 321v BzN Q where E is a row vector corresponding to the 2N th row of matrix A. 45 Now with the constraints in (3.30) and (3.31) we can define a new state R (one dimension less than Q) given by I - .. 0 Q= -1 V ””1".“ -1 n~.1-.]R+[ 1 [v(t) (3.32) V2N,2L “N '1 V2N.2L V2N.2L For simplicity let us denote the coefficients of R and v(t) by S and T respectively. Thus we get Q = SR + Tv(t) (3.33) Substituting for Q, in (3.29) we obtain the following relation: sr’r + Ta(t) = K,,(sn + Tv(t)) + B—l—F,,(a(t) — EV(SR + Tv(t))) 2N (K.. - Lnnwsn + (K,, — —1-F.,EV)Tv(t) + imam 321v 321v 321v Rearranging the terms to bring R to the left side, and premultiplying the equation with (STS)‘IST results in 11 = (sTS)-lsT(K.,—-1—F,,EV)SR 321v +(sTS)-lsT(K., — ernvnvu) + (sTsr‘sT —1—F,, — Ta(t) B2N BZN For simplicity we define the matrices K“, Fm, and Fm to be the coefficients of the terms R, v(t) and a(t) in the above equation. Thus we get the reduced equation for the system when in the stick mode. 1'1 = K.,R + F,,.,v(t) + F..,a(t) (3.34) One can thus build a reduced order model from experimental data which can later be used to study and characterize various phenomenon that occurs in a frictional system. 3.4.2 The order of reduction Now that we know how to construct a reduced order model, the next question is how to decide on the order of reduction? General practice in literature is to pick 46 the first few modes the sum of whose corresponding eigenvlaues exceeds 99% of the sum of all eigenvalues. Cusumano [17, 18] studied the correlation between the 99% criterion and the false nearest neighbours test for dimensionality of the reconstructed phase space. Let us decide on the validity of the 99% criterion in the selection of dominant modes. In Figure 3.7 we show the phase portraits (plots of velocity versus 7. 8, 9, all modes _1 l l l 6 X10 O M m- .5 O .A N Figure 3.7: Phase portraits of 10-SMD system, with different levels of reduction — A qualitative comparison for determining the order of reduction displacement of the driving mass) of different reduced models of the system with different levels of reduction. Note that for the system under consideration, six of the modes have eigenvalues summing up to 99% of the total. We have not shown the one and two mode reduced model phase portraits as they deviate well beyond the range shown in the figure. For quantiative comparison of the full and reduced systems we look at two different measures. One of the measures corresponds in a sense to the total energy in the 47 __Iahl§_3_J;_St_a_ti§tig_al_d§viatignsLfidiflerent levels of reductio # of modes 10 9 8 7 6 5 4 3 2:3;wa 100.0 99.99 99.95 99.90 99.59 98.53 96.50 94.05 Z; -1 ’\' AX 0.00 0.05 0.25 0.39 1.08 2.17 4.67 6.07 AE 0.00 0.01 0.02 0.01 0.06 0.44 0.99 2.07 system. The energy E in the system can be represented by 11 10 - 3:531 .. "—+ —.-x2+ §i(X- 11.-.)” (3.35) i=12i=2 With this relation we can define E f and E, for the energies in the full and reduced system respectively. The deviation of the reduced system from the full system can then be represented by E f - E, E! and this deviation is used as one of the measures. Other measures that we used AE= x 100 (3.36) to compare are the deviations in the displacements and velocities of the individual masses and the sum square of those deviations: 2 2 AX: 1‘1 AX +AX x100 (3.37) _,X3+ X3 This measure would take care of the problems where we have the first measure small in spite of a large phase difference. In the table (3.1) we show the statistical measures of deviation of the models from the full model. From the qualitative comparison we find that with six-mode reduction the orbit is close to that of the full system. From the table we also note that the statistical deviation of the reduced system from the full system is well within the limits. As such the 99% criterion is accepted for selection of number of modes and we continue with the research. We however must mention that for some engineering applications, 48 agreement such as the three-mode case is considered great and settling for two or three-mode cases consideriably simplifies the analysis of the system. 3.5 Validation of Reduced Model Now that we have decided on the order of reduction, the next step is to validate the reduced model. Different methods have been thought of in comparing two non- linear systems. The foremost method is to compare the qualitative behavior of the two systems. The next is to look at the energy content. This gives a quantitative comparison. Another important feature of nonlinear systems is the occurence of bi- furcations and jumps. Towards this we compared the bifurcation diagrams generated from the two systems. To compare the chaotic dynamics, we look at the Lyapunov exponents. Another comparison for chaotic dynamics is the comparison of skeletons of the attractors. 3.5.1 Qualitative comparison For qualitative comparison one can look at the phase portraits of the system obtained from the full and reduced models. Phase portraits of the system are obtained by plotting the displacement versus velocity of the driving mass. With different initial conditions and belt speeds the system gives rise to different dynamics and we plot the corresponding phase portraits from the reduced and full models. In the plots here we show the phase portraits from the reduced model in the top row and those from the full model in the bottom row. One can observe that the phase portraits of the two systems look similar, including the complicated orbits of higher periods, where although there are deviations, the size of the attractor is of the same order. 49 Figure 3.8: Qualitative comparison between full and reduced models. Top row cor- responds to the dynamics of reduced model and bottom row to the full model. Each column corresponds to a set of identical initial conditions. All other parameters are kept constant. 3.5.2 Quantitative comparison For quantitative comparison of the periodic orbits, we can use the same measures that we used in determining the order of reduction. In table (3.2) we list these measures for the three periodic orbits that we used for qualitative comparison. From table 3.2 we realize that the statistical deviations of the reduced model predictions from the full model predictions are rather small, more so from the energy perspective. We do not want to elaborate on this comparison as the tolerances allowed are application dependent. For example, in a computer assembly plant the tolerances could be fractional percentages whereas in a steel plant even 5-6 percent deviation is 50 Table 3.2: Quantitative comparison between full and reduced models. AX and AE defined in equations (3.37) and (3.36) respectively period # AX AE 1 0.614% 0.095% 3 2.527% 1.350% 5 4.873% 0.610% considered acceptable. We also note that in some applications the maximum deviation between the predicted and actual displacements of the system at some particular location would be the deciding factor. An example would be the case of a pick-and- place robot manipulator. The allowable deviation is again dependent on the specific application. 3.5.3 Comparison of bifurcations Another possible method of comparing the two systems is to check for the occurrence of bifurcations. Towards this goal we generated the bifurcation diagrams for both of the models and did a comparison. We varied the belt speed from 0.500 to 0.001 in steps of 0.001. The initial conditions at each belt speed equals the end conditions of the previous speed for the reduced system. Note that the same initial conditions that were used for the reduced model at each belt speed have been used for the full model also. Thus we are comparing the system dynamics in steady state that go close to the subspace spanned by the L modes of the system. Here we present the two bifurcation diagrams. As typical of any two nonlinear systems, we do see differences. The difference is small when the reduced system predicts periodic orbits. However for bifurcations ' into non-periodic orbits, the predictions did not match well with the full model. This could be attributed to a shift in basin boundary and the sensitive dependence to initial conditions of chaotic systems. 51 . o . .0 . . . . . ° . e . ' . . . . . u u. . . . . . 1 s u . o ' ' ' . ‘ ' ' ‘ 6 POM reduced MOI - ........ 'e ......... ' ......... ' .......... ' .......... ' .......... '- n e u n nnnnnnn - . . i . . . . . . . l l l L l 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 2 o 0 ‘ . ‘ . ' I I I 2 . O. I ‘ . . . - - n . . ... 0. .. . . . o . 1 ........ ;. ......... . .......... . .......... . .......... . .......... - ...... original system ...... ; ........ 1 I I o: I I I I Z 2 1 1 1 1 1 '1 1 '1 1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 belt speed Figure 3.9: Comparison between full and reduced models: Bifurcations 3.5.4 Unstable periodic orbits During a chaotic orbit, the system visits countably many unstable periodic orbits [51]. These unstable periodic orbits form the skeleton of the chaotic orbit. One way of comparing the chaotic orbits is to compare these skeletons. Lathrop [51] gives a procedure of extracting the unstable periodic orbits in a reconstructed attractor from a time series. In this method a phase-space attractor is reconstructed from the data using the standard method of time delays as explained in the earlier sections. The attractor is the set of points y(n) = (smsfir, ...,sn+(d_1)1). The time delay T and embedding dimension d are determined as earlier on the basis of mutual information and false nearest neighbors respectively. The orbits are located as follows. Let c > 0, and let y(i) be a point on the reconstructed attractor. We follow the observed images y(i + 1),y(i + 2),... of y(i) until we find the smallest index k > i such that 52 ||y(k) — y(i)|| < e. If such a k exists, we define m = k — i and say that y(i) is an (m, e) recurrent point. In this analysis we fixed 6 = 0.5 x Ay, where Ay is the average distance between two consecutive points on the attractor. In the following we show the skeletons of attractors reconstructed from velocity time series of the first mass, obtained during the sweeps when the belt speed is 0.2, both the systems starting with identcal initial conditions. We find that the approximate unstable orbits have close resemblances qualitatively. —1 o 1 "-1 o 1 '-1 o 1 "-1 o 1 301) s(n) s(n) s(n) Figure 3.10: Unstable periodic orbits in the chaotic attractor reconstructed from the velocity time series of the first mass. The first row corresponds to data from reduced model and second from full model. The first column shows part of the attractor. The remaining shows approximate unstable periodic orbits of periods 1, 2 and 3 53 3.6 Discussions and Conclusions Proper orthogonal decomposition theorem has been used successfully in a Galerkin- based order reduction of a frictionally excited multi-degree-of-freedom numerical sys- tem. This work thus bridges the work of Cusumano [17] which relates the size of a nonsmooth system to spatial coherence of the system and the work of Fitzsimons [34] which deals with modal reduction of simple smooth systems using POD. Also in this chapter we reinforced modal reduction and modal projection in a stick-slip system which is a non-trivial task. Though the reduction is only from an order 20 to an order 12 (this being a discrete system), we suspect a better reduction can be achieved for the same criterion, in case of continuous systems where only a few modes of the system dominate the dynamics. Another important issue is the criterion used for picking the number of modes is application dependent, the same criterion can be better for some applications and can be worse for some others. Detailed analysis of an experimental system with frictional excitation is presented in the following chapter. The reduction in the order helps in applications like building controls viz. the control of robot manipulators where nonsmooth forcing like friction and impact are imperative. Chapter 4 Physics and Measurement Methods of the Experiment In the previous chapter we used POD theory in modal reduction of a discrete nu- merical system. The different validation techniques of the reduced model indicated a great potential for this method of modal reduction as it does not do the violence of linearizaton. In this and following chapters we wish to extend the method to a real continuous system and confirm our hypothesis of using POD technique in a nonlinar and nonsmooth system. 4.1 Description The main ingredients needed to apply the proper orthogonal decomposition theorem to a physical system are the displacement or velocity information at a set of prede- termined locations. Also for using the proper orthogonal projections in a Galerkin’s projection method the physics of the system should be well understood, i.e. the mathematical model for the system should be known in advance. Towards this goal we consider a cantilever beam subjected to frictional excitation. The frictional con- tact is provided by a solid mass of steel connected to a shaker. This system has been selected because it is relatively inexpensive, simple to build and measure, and known 54 55 to show rich dynamics. The displacements of the beam at various locations can be measured using proximity probes, strain gages or accelerometers, and there is abun- dant literature available on the dynamics of beams. The schematic of experimental setup is shown in Figure 4.1. A photograph of this setup is presented in Figure 4.2. ISL—1' 1:. 3315:? 1:. =1 e o 11...... .__A Figure 4.1: Schematic diagram of the experimental setup Before deScribing the setup let us define a coordinate system XYZ. Let X lie along the length of the cantilever, Y the thickness and Z the width of the cantilever. Thus the transverse vibrations of the cantilever is along the Y direction. The cantilever is built by taking a 0.4000 x 0.0128 x 0.00086 m3 mild steel beam (E=128 x 109 N/m“, p=7488 kg/m3) and fixing it at one end. At the free end of the cantilever, we attached a small beam (0.064 x 0.0134 x 0.00056 m“, E = 126 x 109 N/m“, p = 6777 kg/m3) to act as a leaf spring. It is referred to as a loading beam. The mass of the fixture (not shown in figure) is 0.01226 kg. The presence of fixture made the effective length of loading beam 0.076 m. The bending axes of the two beams are perpendicular to each other, i.e. the length of loading beam is along the X 56 Figure 4.2: A photograph of the experimental setup axis, the thickness along the Z axis and the width along the Y axis. The loading beam has low transverse (Z) rigidity and high lateral (Y) rigidity. This configuration allows for the transmission of any force on the tip of loading beam in the Y direction to the cantilever without significant change. A small hemispherical rubber nub is glued to the free end of the loading beam. This nub makes contact with the oscillator. The loading beam is used to provide a spring-loaded normal load. Easy monitoring of the normal load is achieved by attaching strain gages to this loading beam (the amount of strain is directly proportional to the moment which in turn is proportional to the normal force). As the oscillator moves to and fro, the friction force between the rubber nub and the oscillator is transmitted to the beam. We incorporated a tilt 9 in the shaker axis by lifting the rear end. The angle between the rubbing surface and the plane 57 of beam motion allows for a variation in normal load. This would affect the force transmitted to the main beam. In Figure 4.3 we show the force diagram. For ease in understanding we defined a coordinate frame xyz local to the sliding surface as shown in Figure 4.3. Nu a 2N2 ® Figure 4.3: The effect of tilt on the force transmitted From now on, we refer to the force N 2 due to the flexure in loading beam as flezural load and force N2, the force normal to the surface, as normal load. F, (the force tangential to the surface) is referred to as the friction force and Fy (the effective force that excites the main beam) as the excitation force. The relative motion between the slider and the rubber nub results in frictional excitation. The friction force F, is proportional to the normal load N,. The normal load is dependent on flexural load N z. The flexural load is a function of the position of the nub on the surface. Let the distance between positions 1 and 2 be a along Y axis. As the nub moves from position 2 to 1, the flexural load would increase by 3E1 a tan 0/l3. In the figure we have represented the loading beam by a spring of stiffness 3E1 / l3, 1 being the length of the loading beam. Using a simple force balance, we have N2 = N2 c039— F, sin 0. Since F, = pNz, we have N, = Nz/(cos d—p sin 0). 58 Now we can write the resulting excitation force Fy on the beam in direction Y as sin0+pcos0 Fy = N, sin9 + F, cos0 = N2 (4.1) cosd—psinfl Thus the effect of tilt is to change the effective coefficient of friction from p to sinOig c000 coca-usina' During the frictional excitation, the dynamics of the beam comprise of low- frequency drifts. While selecting the displacement measuring devices, care should be taken to see that their low-frequency limit is below this drifting frequency. For example a B&K laser transducer (type 8323 with power supply 2815) when used in the displacement mode has a frequency range of 0.3 Hz - 20 KHZ and therefore cannot measure the“ drifts below 0.3 Hz accurately. Also as displacements at different loca- tions are needed one needs more than 4 - 5 displacement measurement transducers which would therefore be constrained by the resources. Now that we have defined the system and the forces of interest, the problem of conducting the experiment can be divided into a problem of measuring the friction, the normal load variation caused by the tilt, and the displacement. But to know the system more accurately we need to know the system response. In the following section we present the system characteristics in the form of impulse response spectra. 4.2 System Characteristics Neglecting the possibility of microscale loss of contact and impact, the system has two dynamic states. One is the case of the rubber nub slipping over the slider surface, when the system acts a forced cantilever. The other is the case of rubber nub sticking to the slider surface, when the sytem more like a clamped-sliding beam. As such we present two power spectra, one when there is no contact and other when there is contact and the normal load between the contact elements is high. The power spectra are obtained by using a Zonic signal analyzer. The power spectra of the strain signals 59 (strain gage located at 0.05 m from the fixed end) obtained by averaging ten impulses close to the fixed end define the system characteristic. Figure 4.4 shows the power spectrum of the cantilever setup and Figure 4.5 shows the system characteristic when there is a high normal load between the surfaces. The first six modes lie within 300 Hz. I I I I I j _30 . ............................................................... _ 40 . .......................................................................... .1 > %_50b ....................................................................... _. .5 §%. .......................................................................................... -70. ........ f, ......... g .......... g .......... g ......... g ......... ‘ ......... g ..... ......... g ....... - _m|. ................................................................................................ .- _m l l l l M l 1 1 l 0 50 100 150 200 250 300 350 400 450 500 frequencylnHz Figure 4.4: Average power spectrum response to impulses close to the clamped end of the system with zero normal load at the contact. The first six modes of the cantilever setup have frequencies of approximately 3.75, 18.4, 53.75, 102.5, 172.5 and 265 Hz. We have only six strain signal conditioners available for strain measurements and hence are restricted to six modes. In the following sections we present the details of measurement techniques. 60 _55... ............................................................................................. .. .60.. ........................................................................................... -65 1 1 1 1 1 1 l I 0 50 100 150 200 250 300 350 400 450 500 lrequencylnHz Figure 4.5: Average power spectrum response to impulses close to the clamped end of the system with high normal load at the contact point. 4.3 E‘iction Measurement A good friction characteristic is required for the validation of the reduced model. The friction characteristic represents the relation between the friction force, normal force, relative velocity and displacements between the rubbing surfaces. The friction force is most commonly modeled as being directly proportional to normal load. However the relationships between friction force, relative velocity and relative displacement is still a subject of fundamental research and is dependent on many parameters. As such we would like to carry out the studies with the surfaces that are being used. A schematic of the apparatus used to get this friction characteristic is presented in Figure 4.6. The apparatus consists of the friction surface (that has been used in the 61 main eXperiment), two slender beams (similar to the one used in the main experiment) with rubber nubs attached at the end and a carefully designed fixture to connect these beams to the shaker (B&K model 4808) to minimize the moments due to normal loads. Strain gages on the beams are used to sense the normal load. As shown in the figure, the sliding contact surface is now fixed to the steel bracket through a load cell. Thus the axial component of any force on this mass can be measured using the load cell. The two rubber nubs attached to the beams slide over this block, exerting normal (to the surface) and frictional (along the surface) forces onto the block. The fixture is mounted on to the shaker through a roller hearing so that the normal forces balance each other and do not produce a resultant moment which would contaminate the force transducer measurement. The force measured by the transducer is therefore the force that gets transmitted to the oscillator. For determining the friction coefficient we need the normal load acting between the surfaces. This can be measured using the strain gages attached to the beams at appropriate locations. The strains measured are proportional to the moments in the beams at the gage locations which in turn are directly pr0portional to the normal load. Thus one can obtain both the friction and normal loads. Wm) H fl H I _ ‘1 uxsuu PCB we.“ U Pine-woodblock swam Figure 4.6: Configuration of friction measurement setup 62 The normal load has been calibrated by suspending known weights from the free end of loading beam with its fixed end clamped in a vice. For example, with a weight of 0.20 lbs, the strain signal corresponded to 3.82 Volts on the output of the signal conditioner. This way of calibration would automatically take care of the excitation voltage, the conditioner gain and gage factor. Thus for the strain gage voltage we have calibration constant equal to 4.29 Volts/ N. We attached the loading beams to the fixture with the solid block wedged between the two beams. This fixture was attached to the shaker, while the solid block was mounted on to a rigid bracket. The axis of oscillation was parallel to the rubbing surfaces. This we cross-checked by turning on the shaker and examining the varia- tion in the normal load as the rubber nubs rode on the surfaces. The normal load variation was found to be 0.8%. The sum of the two normal loads (taking the sign into account) is the effective normal load between the sliding surfaces. During the friction measurement the mean normal load was found to be 1.86 N. When measurements are taken with different devices there is a scope for phase differences between the various measurements due to the inherent dynamics within the devices, more so when post filtering is done on one of the measurements and not on the other. At 10 Hz frequency the phase difference between measurements of accelerometer PCB 309A and force transducer PCB 208B was found to be 50.4 degrees (0.88 radians) with the accelerometer signal leading. This phase has been adjusted in our measurements. (Incidentally in the above phase calculations we also took into account the phase change that would arise due to a high-pass filtering of the accelerometer signal. We noted that the accelerometer PCB 309A has a low-frequency limit of 5 Hz. Between 5 and 10 Hz there is -3.2 percent deviation in the amplitude. We have accounted for this during our calibration procedure described below.) For measuring the input displacement and velocity we made use of the accelerom- eter. We first obtained the acceleration using accelerometer PCB 309A. As we are measuring the acceleration of the input whose frequency is known, we high-pass fil- 63 tered the acceleration signal to take care of any drifts due to bad ground looping. We took a Fourier transform of this signal and divided this by the frequency compo- nent once to obtain velocity Fourier transform and twice to obtain the displacement Fourier transform. We took the inverse Fourier transform of the two resulting signals to get velocity and displacement. We cross checked this information using a laser transducer B&K 6808 with power amplifier B&K 2815. The laser transducer has a low-frequency limit of 0.3 Hz in the displacement measurement mode. After running the shaker at 10 Hz for about an hour, we recorded the data for 32 seconds at a sampling rate of 1000 Hz. We used Laboratory WorkBench on a Masscomp 5550 for this data acquisition. The four channels of data we collected are the acceleration, both the normal loads and the friction force. In Figure 4.7 we plot a sample of the data after processing for phase correction. A plot of the relative displacement (of the shaker) versus friction force is shown in Figure 4.8 and relative velocity versus friction force is shown in Figure 4.9. The mean slope of the nearly vertical lines of direction reversal in the displacement-friction plots gives the contact stiffness Kc. This slope can be used to predict the velocity with which the rubber nub slips and this has been discussed in the following paragraphs. From the friction-displacement characteristic in Figure 4.8 we find that the effec- tive contact stiffness is —20 kN/ m (this is the average slope of the lines of direction reversal, the standard deviation is 0.1% of the average). A better explanation is given in the following paragraphs. The hysteresis in the friction-velocity plot in Figure 4.9 suggests that the friction characteristic may be close to a Coulomb model with contact compliance. The kinetic friction force limit is close to 1.03 N. Dividing this by the mean normal load of 1.86 N gives us the kinetic friction coefficient of 0.55. The static to kinetic transition if it exists, is generally at high frequencies and gets filtered. As such the static friction limit is difficult to get captured in experiments of this nature and needs a different but simple method. 64 Q E 5 O .. .. 3 . . > I I I . . .02 1 1 1 1 1 1 1.2 1.4 1.6 1.8 2 2.2 2 I I I I I Z 5 S 0 ‘ B 'C "_2 1 1 1 1 1 Z 1 1.2 14 1.6 18 2 2.2 .5 I I I I I I I T I I I '01-865 .................................... .......................... - .3 .16 .. _. 21.855 ..... ......................... .......................................... ....................... .. 3 i 1 1 1 1 1 1 1 1 1 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 ' 2 2.1 limeins Figure 4.7: Sample time-series of displacement, velocity, friction force and effective normal load In a separate experiment we rested a mass of 250 grams underwhich we attached three rubber nubs (similar to the one used in the experimeint) on an inclined plane made of a similar surface (as the one used in the experiment, but not the same one). We increased the tilt of the inclined plane until the mass began to slide. The average tilt at which the mass began to slide was 37.2 degrees (standard deviation across seven different readings was 2.7% of the mean) which amounted to an average coefficient of static friction p5 = 0.75 (standard deviation of 115 is 3.7% of the mean). Thus we have 115 z 1.38pK. In our numerical experiments we have used 1.36. Incidentally we conducted this Madmlorcehhl 65 1.5 —o.s-~- -1--- I I I I I I I I I /’ ' . ~ ......................... L.........,..........,.........-.................... .— .............................. - .... .. .. ...; ................................................... ..1 .1. . . . ... u ‘ I / [K ’ \ ’.\‘ ‘ . . v I ,. . 1 » 0 l 1 l 1 l 1 1 1 "3.5 Figure 4.9: F riction-velocity relation 66 experiment for determining static friction to verify the validity of the estimated value of 1.36. We will describe the method that we adopted for the estimation of static friction later in Chapter 6. Both the friction-displacement and friction-velocity plots indicate “hysteresis” phenomenon. This has been described as “spring-like” sticking behavior in the litera- ture (Canudas de Wit [13]). Experimental evidence of this kind of behavior is reported by Dahl [19] and Courtney-Pratt and Eisner [15]. Both the reports indicate a junc- tion compliance which induces small displacements prior to the actual slip. Eflorts for modeling these elastic contact problems have been put forth by many researchers (Liang and Feeny [52], Oden and Martins [63], Canudas et al. [13], Harnoy [39]). Liang and Feeny aptly described this motion as “microsticking” as during this stick- ing event the elastic displacement is really small for high contact stiffness. Dahl’s model is equivalent to Coulomb friction model incorporated with a lag in friction force change when the direction of motion reverses. Canudas de Wit et al. treated the friction interface as a contact between bristles. Their model was able to describe Stribeck and frictional memory effects in the sliding regime. .Xcl‘) .Xrl') Kc Shaker attachment A sin w t I ' W ' lt-—-—. Rubber “ELK j / é / / Steel block 2 / 2 Figure 4.10: A schematic diagram showing a massless compliant-contact model A schematic of the model similar to the one used by Liang and Feeny is shown in Figure 4.10. This model effectively describes the observed phenomenon in our experiment. YT(t) describes the displacement of the shaker. Yc(t) represents the 67 displacement of the contact point (the point is on the rubber nub). The dynamics of the block and the mass of the nub are negligible. However due to contact compliance there is a relative motion between the contact point (on the rubber nub) and the shaker attachment. During pure sliding motion, the relative motion VT(t) — You) is negligible since the contact is nearly massless and its dynamics are damped out immediately. When Vc(t) remains zero, we have sticking event. This does not neces- sarily mean IQ“) is zero. I’T(t) continues to differ from YOU) which is zero, until the static friction force can no longer sustain the stiffness force exerted by the compliant- contact joint, i.e.IKc(YT(t) — Yc(t))| > IF,|. When the contact elements slip we again have YT(t) = Yc(t). F, is the static friction force bound. Thus the slope in the friction-displacement plots in Figure 4.8 is the effective contact stiffness Kc. From the plot we determined the contact stiffness to be -20 kN/m (the average slope of the lines of direction reversal was found to be 20.2 kN/m with a standard deviation of 10.6% of the average). Now we know from the friction-displacement plot the change in YT(t) during the microsticking event. We also know that micro-sticking begins when the relative velocity is zero, in our case YT = Vc(t) = 0 and ends when stiffness forces overcome the friction forces. Let t1 and t; be the times when microsticking begins and ends respectively. Therefore the relative displacement YT(t2) - YT(t1) is effectively (F, — F1) / Kc. Let us assume harmonic excitation of the shaker attachment, i.e. YT(t) = A0 + Asinwt. Then V701) = Vc(t1) = 0 as prior to t = t1 we have VT(t) = Vc(t). This implies YT(t1) = A0 + A. Also YT(t2) = A0 + Asinwtg. This in turn equals YT(t1) + (F, - Fk)/Kc = A0 + A + (F, — Fig/Kc. Thus we have an estimate of t2. From this we calculate 17702). Going through simple mathematics we get YT(t2) = —wA,/1 — (1 — (F. — F.)/KCA)2. From the friction-displacement plot we have A = 0.0023 m, Kc = —20 KN /m, F. = 1.03 N and F), = —1.0 N. We know the value of w = 201r rad/s. Substituting all these parameters in the equation for V102) we arrive at a slipping velocity equal to 0.044 m/s which is consistent with the the friction-velocity plot, in which a microstick— 68 slip transition occurs at approximately V, = 21:0.045 m/s. 4.4 Displacement Measurement One can obtain the displacements of the beam at various locations directly and in- directly by using proximity probes, laser transducers, accelerometers, strain gages, etc. As we need many probes, the cost of laser transducers preempt their usage. (Though we can use two point correlation measurements using two laser transducers, this would be very difficult to achieve, as we also intend to operate in the chaotic regime.) Moreover the frictional excitation of the beam results in low-frequency drift of the beam over the slider which cannot be measured by the laser transducer due to its low frequency limit. The chaotic regime is usually in the low-frequency range and we find a need to measure the low frequency dynamics of the system. This puts a restraint on the usage of proximity probes and piezo-accelerometers. With the above considerations in mind we are left with the option of strain gages which are relatively inexpensive. The problem then would be to convert the strain data to displacement data and from the discrete displacements to continuous. We used six half-bridge strain-gage circuits to obtain the strains at six uniformly distributed locations. We were limited to the number six for two reasons. The main reason was the availability of signal conditioners and the other was the finite length of the beam being used. We attached the strain gages on either side of the beam to form a half bridge. We spaced the gages at 0.05 m, with the sixth gage located at 0.05 m from the fixed end. Based on numerical simulations we found that this distribution helps us maximize the strains. For easy reference we numbered the gages from 1 to 6, strating from the gage at 0.30 m and ending at 0.05 m. All the strain gages used are of the type Micro-Measurements Precision Strain-Gages model CEA- 06-250UN-350. These have a resistance of 350.0 :1: 0.3%, a gage factor 2.1 :1: 0.5% and a transverse sensitivity of (+0.1 :1: 0.2)% all at 24° C. For minimizing the eflect of wiring, we used Micro-Measurements single conductor wire type l34-AWP. To reduce 69 noise we twisted the three wires coming out of the each half bridge and tried our best to minimize the length. We ran the wires along the beam to the fixed end. We used four strain-gage conditioners of the type 2120 and two of the type 2210, both of which are from the Micro-Measurements Group Inc. We used an excitation voltage of 10 V for both the models. A gain of 2000 was used on model 2120 and gain of 1000 was used on model 2210. The difference in gains have been accounted later in the post-processing programs. For data acquistion we used Laboratory Workbench of Concurrent Computer Cor- poration on a MASSCOMP 5550. We acquired six strain signals on the main beam, one strain signal of the loading beam and one input acceleration signal using an ac- celerometer PCB 309A. We used a sampling frequency of 1000 Hz for the acquisition. 4.4.1 From strain to displacement We know from basic mechanics that the strain 6(1') in a beam subjected to pure bending is related to the transverse displacement y(z) by 62y 633 = 6555 (4.2) where c is half the width of the beam. One of the methods of solving a PDE is to convert it into a ODE and then solve the 'ODE using standard ODE solvers. The conversion is achieved by a Galerkin’s approximation by assuming y(a:, t) = 2;, ¢,($)u,-(t), where 453(3) form a basis satis- fying the geometric conditions. We use the same technique here which enables us to write 6,, in terms of 4).- and 11,-. Thus "I(t) c,,(:c)=c[¢1(x) 1141)] g , (4.3) un(t) where 1b,-(:1:) = 32¢,(z)/02:“. Now by measuring 6 at n diflerent points on the beam 70 we get n such equations. Writing these 11 equations in a matrix form, we get 61 $11 ' °° $111 111 ' = c '- , 5 , (4.4) 6n 1an ° ' ° 1PM ”it where, 6.,- indicates the strain measured at the ith location 2:3, 1b,,- = ¢j(:c,-) is the value of ¢jth function at the ith location and u, is the corresponding weight. The matrix on the right hand side is invertible and therefore we can compute the values of ms. This can be used in the Galerkin’s approximation formula to get back the actual value of y at different locations. Thus we get —1 3’1 4511 ¢1n I(’11 $111 61 fill--| . .. .. . (4.5) yn ¢nl <firm zpnl $1111 611 In our displacement calibration experiment, we applied a high normal load at the contact so that the contact remained in the stick mode at that operating frequency. We operated at 8.45 Hz frequency during this stage. We measured the displacement of the beam at a location of 0.35 m from the fixed end using the laser transducer and simultaneously obtained the strain signals. We then used equation (4.5) to obtain the displacement at the same location. The laser transducer sensitivity is 1 mm/ V. We took Fourier transforms of the laser displacement and the computed displacement and took the ratio of the amplitudes at the operating frequency. This gave us the calibration constant needed in our procedure. This constant should equal half the thickness of the beam. At 8.0 Hz we found this constant to be 99.9% of the thickness. At an operating frequency of 12.0 Hz (with normal load high enough to result in a pure stick motion), we found the constant to be 99.84% of the thickness. We used this calibration constant to verify the displacements at six different locations. Across six different locations, the standard deviation of these constants was found to be 0.71% of the mean of 98.57, with the worst match occuring close to the fixed end where the displacements are very small and we suspect the laser measurement to be not so accurate. 71 The next step was to verify the above procedure during our normal operating conditions where the dynamics would comprise stick-slip oscillations and the dis- placement is not a pure harmonic. We cross-checked our procedure by bringing the normal load to the operating normal load and operating at 12.0 Hz and 41.1 Hz where we found the dynamics to be “chaotic”. During this verification at both the frequencies we noted that the calculated displacement consisted of harmonics of fre- quencies higher than the sixth frequency of the system, which are not present in the laser measured displacement. This might be attributed to the fact that we are using only six modes to approximate the beam dynamics. We use only six modes as we have only six strain signals at any given instant. As such we low pass filtered the calculated displacement, with a cut-off frequency of 300 Hz. We also noted that there is a low-frequency deviation between the two signals. This we attributed to the fact that the laser cannot measure the low-frequency drifts (the low-frequency bound on the laser in displacement mode is 0.3 Hz), while the strain-gages can. As such we high-pass filtered both the laser and calculated displacements. The cut-off frequency at this stage was 0.5 Hz. Now both the signals are ready for comparison. The strain signals have been corrected for by the calibration constant of 0.9857. Figures 4.11 and 4.12 show samples of the direct and indirect measurements (at an operating frequency of 12.0 Hz) in time and frequency domain respectively. At this operating frequency, for the uncorrected signals we found the ratio between the peaks in the frequency domain to be 99.98%. Figures 4.13 and 4.14 represents samples of both the signals obtained at an oper- ating frequency 41.1 Hz. in time and frequency domains. At 41.1 Hz operating frequency, the direct and indirect measurements did not match as well especially in the high frequency regions (above 120 Hz). This could possibly be the result of using a simple pure bending theory for the relation between strain and displacement, without taking into account the effects of shear. Also the ratio between the peak amplitudes of power spectra at this operating frequency is 72 010$ Y&y N L I» time in s Figure 4.11: Sample time series of direct (laser) and indirect (strain gage) displace- ment measurements (Y and y) at an operating frequency of 12.0 Hz. Ill 888° 20"09100110' and Y)) L O -50 Figure 4.12: Fourier transforms (averaged) of direct (laser) and indirect (strain gages) displacement measurements at an operating frequency of 12.0 Hz. 81.11%. Another possibility for this discrepancy is the presence of modes higher than the sixth mode. This can be corrected by using more than six modes during the curve fit. This is a later thought and we have not attempted implementing this idea. We noted that the noise floor in the frequencies above 50 Hz to be 50 dB below zero (see Figure 4.15 where the FFTs of noise floor for the measurements of displacement using laser and strain-gages have been depicted). The signals are found to be of close match, apart from subtle details, and we 73 x10 1 g “I 0.5 ' ‘ . . >~ oi| ‘ i if a I ll 1 1 I l, 1 I >—o.5 . ~1 , p . , —1.5 ; , ,' 1 > ' , ' 1 ~ - oiect- , ,—— l irect - 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 tlme ln s Figure 4.13: Sample time series of direct (laser) and indirect (strain gage) displace- ment measurements (Y and y) at an operating frequency of 41.1 Hz. 20009100110! and y)) 1 50 frequency In Hz Figure 4.14: Fourier transforms (averaged) of direct (laser) and indirect (strain gages) displacement measurements at an operating frequency of 41.1 Hz. adopted this procedure for getting the displacement information in our POD studies. 4.4.2 From discrete to continuous In the previous section we have described a simple way of obtaining displacements of the beam at specified locations by using strain gages. However the beam is continuous and we do not want to restrict our studies to only those few locations where we 74 20'|0910(ffl(Y)) zo'Iogtoauw» s s 2 Figure 4.15: Fourier transforms (averaged) of the noise floor levels for direct (laser) and indirect (strain gages) displacement measurements measure the strains. Here we discuss a way of converting these discrete signals into a continuous signal. Though this is a simple curve fit, we present this for the sake of completion. We can do a simple polynomial fit as suggested by [34] or in this experiment as we know the linear modes of the cantilever, a modal fit can be obtained by using the linear natural modes (LN Ms) of the cantilever. We can later do a simple Gram-Schmidt orthonormalization method to ortho-normalize the functions. Let [ 1:1 - - - 2,. ] represent the 12 locations where we know the displacements. Let the corresponding displacements be [ yl - - - y" ]. Note that we are not taking into account the zero end condition where we have zero displacement all through the motion. Even if the motion is not zero at this end condition we are not taking this point to avoid singularity problems. This problem, if present, can be tackled through a simple translation. Let us try to fit an nth degree polynomial along these 12 points. z.e., y(a:) = tan" (4.6) i=1 where a.- are the coefficients to be determined which can be obtained by a simple 75 least squares fit (“polyfit” of MATLAB would do the job). Once the coefficients a.- are determined, we can substitute into (4.6) to get the value of y(a:) at any given an. Instead of doing a polynomial fit, we can make use of the natural modes of the system to do the curve fitting, i.e., we fall back on our usual Galerkin’s approximation: y(:z:) = Z a;¢,-(:r), (4.7) i=1 where a,- are the coefficients to be determined and ¢,~(a:) are the linear natural modal functions. Writing the above equation for the 72 locations in a matrix notation, we get ¢1($1) '" ¢1($n) [yl y..]=[a1 an] (4-8) ¢n(:c1) ¢n(zn) ¢1($1) ° ' ' ¢1($n) [a1 an]=[y1 ya] (4,9) ¢n(31) ¢n(:rn) --1 We prefer to use natural modes as it would give us a better understanding of the natural modal components in the POMS later. Also, they are natural to the linear part of the system. 4.5 Measurement of Normal Load Variation The case of a single degree-of-freedom oscillator shows us the variation in normal load would give us rich dynamic behavior [28]. As such we introduced a tilt 0 in the shaker axis to allow for normal load variation as the beam vibrates. This normal load variation has to be measured for modeling and we used a quasi-static approach during our measurement. We positioned our beam at various locations on the slider surface and at each location (after waiting for at least fifteen minutes for the system to settle) took the strain measurements at different locations on the main beam and also the strain on the loading beam. The strain on the loading beam represented the 76 normal load. We used the same procedure that we described in the section on friction measurement to convert the strain signal into the normal load. In the previous section we described a way of obtaining the displacements at a specified location on the beam using the strain data. We can use a similar procedure to obtain the displacement of the rubber nub K on the sliding surface. We just need to add I * Y}, to Y1, to get K. Here I is the loading beam length, and YL and YI’, are the transverse displacement and slope of the main beam at the “free” end. 0.38 r u I t I I 20.37 .............. .............. .............. .............. ............ .. 50.36 ............. .............. slopes—JSQN/m ............ d E 3 § i hominalnorinauoad=as45N 0.35. ............. 1 .............. I .............. .0 ..... .............. ............ .- EQM_ ....... ”mammmgmmmms .............. ;mm.”: ............. g ............ _ O : : ; : . : c033,...n.........f .............. 3 .............. 3 .............. i .............. 3. ........... i ............ q ' 2 § i E § 2° 0.3 l 1 1 l l l .02 -0.0‘l 5 -0.01 -0.005 O 0.005 0.01 0.015 displacment in m Figure 4.16: Relation between normal load and position of the beam tip on the slider. In Figure 4.16 we show a plot of normal load versus position on the slider surface. A least-squares fit gives us a slope N gm of —1.54 N / m and the nominal normal load N Zc (at zero position) of 0.345 N. Earlier we have seen that the normal load N z varies as 3E1 a tan 0/13, where a is the change in displacement. By substituting the values of E, I and I we get the value of 0 = 0.312 deg. This value when used in equation (4.1) results in excitation force Fy z 0.557Nz for a coefficient p = 0.55. However N z is still computed using the relation N z = -l.54YT + 0.345 N, where YT is the displacement of the loading beam tip. In terms of the loading beam transverse displacement Z, the above relation can be written as Z = ZmYT + Zc = —5.45e3YT + 1.2263. These numbers will be used later in the mathematical modeling of the system. Chapter 5 Beam Dynamics, Reconstruction and POD 5.1 Introduction Frictional excitation of low-order physical systems have shown a wealth of dynamical behavior [28, 30]. In Chapter 3, we have discussed the rich dynamical behavior of higher-order numerical systems with frictional excitation. In this chapter we try to explore the different possible dynamics of a continuous beam. In our exploration we carried out amplitude and frequency sweeps to understand the different dynamics and the routes. We will first present the bifurcation studies with frequency and amplitude as bifurcation parameters, followed by delay maps of a few interesting cases. Later we will present the phase-space reconstruction results to characterize the chaotic dynamics. We end this chapter with a proper othogonal decomposition of the chaotic dynamics and make an attempt to explain the spatial aspects of the dynamics. 77 78 5.2 Bifurcation Studies The beam in the experimental setup described in the previous chapter is subjected to frictional excitation from the harmonic slider. In the case of pure stick we expect the beam to show harmonic response. At any given frequency as the input amplitude is increased the static friction force is overcome by the restraining force due to the elasticity in the beam and the contact surfaces slip. Similarly for a given input amplitude, as the frequency increases the surfaces slip. However in carrying out a controlled experiment we are constrained by the ranges of the equipment in use. One of the ways of exploring the dynamics of a given system is to carry out a sweep along a well defined curve in the parameter space. The sweep could be an amplitude or frequency sweep with all the other parameters kept constant. The amplitude sweep is easily manageable whereas the pure frequency sweep is difficult to obtain because the shaker amplitude decreases as the frequency increases in order to maintain a constant impedance. However a frequency sweep gives us an indication of the zones of rich dynamics. As such we carried out a sweep of the frequency of the input to the shaker. In Figure 5.1 we show the amplitude versus frequency of the shaker output and a tenth-degree polynomial fit. Thus a frequency sweep takes us along this curve in the parameter space. Incidentally, we later realized that the relation between amplitude and frequency is exponential. During the sweeps, the two signals that we made use of are the acceleration of the slider (input) and the strain at location 6 on the beam. A LabVIEW program on a Macintosh has been used to run the sweep. We sampled both the signals si- multaneously (to keep the phase constant) at a sampling rate of 48 samples/ cycle and omitted the first 500 cycles to account for transients. We collected the next 100 cycles of data and processed the data to obtain 100 “quarter phase strain” signals. By quarter-phase strain, we mean the strain value when the input acceleration is at the peak and the excitation velocity is passing through zero. We carried out this operation at 450 slider frequencies uniformly spaced between 5 and 50 Hz. Figure 5.2 79 driving frequency in Hz Figure 5.1: Amplitude versus frequency of the shaker with other parameters kept constant. A tenth-degree polynomial fit has been used to curve fit the data. Input frequency sweep of the shaker thus takes us along the above curve on the output. shows the resulting bifurcation diagrams in the forward and reverse sweeps. The time for the sweeps took about 6 Hours. In the bifurcation diagrams of Figures 5.2, at each frequency we have 100 points. A single cluster of points (ideally the cluster should be of zero dimension) at any given frequency indicates the response is periodic with period one, multiple clusters indicate higher periods. A smear of points would indicate either a quasi-periodic or chaotic orbit. In a latter section we show the delay maps obtained at a few selected frequencies to elaborate on the same. Notice the hysteresis in start and end of different regimes during the reverse sweep. This is typical of nonlinear systems. A wealth of dynamics is observed near 12.0 and 40.0 Hz. As such we carried finer frequency sweeps in these zones. These are presented in Figure 5.3. The rich dynamics are more pronounced in the latter zone which was accompanied with slight audible chatter. First we suspected normal vibrations to be the cause of this chatter and made an attempt to visualize the loss of contact using a stroboscope. The stroboscopic probe did not show us any loss of contact. We suspect the jump from the static friction to kinetic friction to be the cause of this chatter. 80 x 10" Q n I I I I I I j 3 : iii - 1 c . . s 1 3 ‘5 . . as i i 2 t . . 1 Q 2 I I o‘ 1 1 1 1 1 L 1 1 5 10 15 20 25 30 35 40 45 50 siiderirequency in Hz at 10" o I I I I I I I I § 1 I G E 0.5 g 0 fi-0.5 1 g -1 § slider frequency in Hz Figure 5.2: Frequency sweeps (forward and reverse): the beam response is denoted by the strain at the location 6 and measured at quarter phase of the input acceleration. The bifurcation parameter is the input frequency. The deviation in the plots during finer sweeps might be attributed to the wear on which we did not have any control. To get a better understanding of the dynamics of the beam at frequencies 12.0 and 41.1 Hz where the dynamics seemed chaotic, we carried out amplitude sweeps, similar to the frequency sweep. The only difference is that the shaker-input amplitude sweep implies shaker-output amplitude sweep. We varied the shaker input amplitude between approximately 1.0 to 5.0 mm. All through the sweep we kept the frequency at 12.0 Hz. In Figure 5.4 we present the amplitude sweep plots at 12.0 Hz frequency. A similar amplitude sweep performed at 41.1 Hz is presented in Figure 5.5. 81 ..A C) h) ma: 0 ON L quarter phase strain at loo 6 Q quarter phase strain at loo 6 I M 11,!» dl‘) O 13 45 1 1 12 40 slider frequency in Hz slider frequency in Hz Figure 5.3: Finer frequency sweeps around 12 and 41.1 Hz : the beam response is denoted by the strain at the location 6 and measured at quarter phase of the input acceleration. The bifurcation parameter is the input frequency. .l" 01me .0 m—l quarter phase strain at be 6 C) .5 3 3.5 slider amplitude in mm Figure 5.4: Amplitude sweep at 12.0 Hz: the beam response is denoted by the strain at the location 6 and measured at quarter phase of the input acceleration. The bifurcation parameter is the input amplitude. The amplitude sweeps at both the frequencies indicated that the system behavior is nonperiodic at most of the amplitudes. In the next section we will characterize the dynamics qualitatively. The amplitude sweeps at both the frequencies indicated that the system behavior is nonperiodic. To verify the effects of wear on the broad dynamics of the we carried out sweeps 82 quarter phase strain at be 6 1 1.5 2 2.5 3 3.5 4 4.5 5 slider amplitude in mm Figure 5.5: Amplitude sweep at 41.1 Hz: the beam response is denoted by the strain at the location 6 and measured at quarter phase of the input acceleration. The bifurcation parameter is the input amplitude. over a finite time by keeping all parameters in our control constant. In a way we are attempting to study the wear effects. We sampled the 100 quarter-phase strains after every 500 cycles for about 30,000 cycles. We carried these sweeps at frequencies of 12.0 and 41.1 Hz (the frequencies where the dynamics seemed chaotic). We present these sweeps in Figure 5.6. 5.3 Characterization of the Dynamics In the frequency sweeps, we note that in the vicinity of 12.0 Hz frequency, the system shows a deviation from a periodic orbit and there is a change in sign of the strain. A closer observation in this frequency range indicates as the frequency is increased, the system seems to go from period one to chaotic and back to period one through a quasi-periodic path. This is speculative and has not been rigorously quantified. Similarly there are interesting dynamics in the vicinity of 40 Hz. In this zone as the frequency is varied from 35 Hz to 45 Hz, speculatively the system goes from period one, quasi—periodic, period two, quasi-periodic, chaotic, quasi-periodic and period one in that order. The delay maps obtained at some of the frequencies provide some 83 to § 1 G g? 0.5 8 o .. ............................................................................................... _. §_0.5 .. .................................................................................................. .. g _1 . ................................................................................................ q a : 3 1 i . . O 0 5 1 1.5 2 2 5 3 x 10‘5 x 10“ m l T l T T 3 . if .E E ‘5 ‘43 . g E a 1 . 8' 1 l 1 1 l O 0.5 1 1.5 2 2.5 3 ll number of cycles x 105 Figure 5.6: Wear sweeps at 12.0 and 41.1 Hz. All other parameters (under our control) are kept constant. visualization of these different dynamics. We also observed very low-period motions at these excitation frequencies. We wish to explore these low-periodic motions in our future studies. The delay maps are the projections of phase space reconstructions. The construc- tion of delay maps is straight forward and involves plotting of a time series with the signal on one axis and a delayed version of the same on the other. For the purpose of visualization, we used a time delay of one-fourth period, which converts to 12 sam- ples. In our reconstruction studies we use a different time delay which will be chosen on a more rigorous basis. 84 no" a no“ b no“ d 6 ........i ......... 6' ...... '. ...... 6 6 ...... . ...... 4 4, . . 4 4 ..... 3T 2 2. ..... I ..... 2 2 ._. ...... E 0 fl ....... 0%....” 0 0 ..... O ..... n_2 ............ _2 . ..... _2 _2 .._'. ,. _4 ......... _4 .. ..... 1 ...... ,. _4 _4 ........... _6 ................ _e.~ ...... , ...... ... _6 ............. —‘D 0 -5 O 5 -5 0 5 x10" x1o" x1o" x_1__o" ° no" no" " 6 .. ...... . ...... 6“ 6' e ...... . ...... 4,. ..... ...... . 4, 4. 4. . a 2 , , ..... .. 2 .. 2 2 ...... ..... E o. a 0.. o. o ........ O ..... ”-2. _2 _2. . . _2 ........ ...... _4.. ...... _4. _4_ _4 .............. _6 .. ............. _6. ................ _B. ...... j; ...... _6 ................ O 5 —5 O 0 5 “made" “”11 10" '0”): 10" “"lx 10" Figure 5.7: Delay maps at a few selected frequencies on the bifurcation-diagram elucidate the periodicity in the beam. The frequencies are: (a) 6.3, (b) 10.0, (c) 12.0, (d) 26.0, (e) 37.7, (f) 39.0, (g) 41.1 and (h) 47.0 Hz. Clearly we can identify the period one, two and non-periodic motions in the de- lay map. To determine whether the non-periodic motions are chaotic or not, we can look at the Fourier spectrum for the broad-band characteristic and confirm the same by looking for positive Lyapunov exponents. The latter involves phase-space reconstruction and dimensional studies using time series data. The delay maps presented above were built using the strain signal at location 6 on the beam. To appreciate the end effects of higher-order dyanmics we present the phase-portraits using the displacement and “velocity” of the contact point of the beam, at a few selected frequencies (11, 12 and 41.1 Hz). We measured the displacement using the method described earlier. We used a simple forward-difference scheme to differentiate these displacements to determine the velocity of the beam at the contact point. We have not accounted for any phase deviation that may result due to this process. We present the resulting phase-portraits and the corresponding 85 power-spectra in Figure 5.8. 1 0.5 ........... ..... o 1. . _o.5 ....... ' -0.5 i -1 ; -0.002 0 0.002 —0.01 0 0.01 -0.005 0 0.005 tip displacement in m tip displacement in m tip displacement in m _1.5,.. .. IogtO(abs(lft(YT))) -2 0 501001500 50100150 frequencyinHz frequencyini-iz Figure 5.8: Dynamics of the contact point of the beam. The top row indicates the phase-portraits of the contact point at frequencies 11, 12 and 41.1 Hz. The bottom row indicates corresponding power-spectra. The dynamics at 11 Hz were found to be periodic and is presented in this study as a reference for noise level. The dynamics at 12.0 and 41.1 Hz are non-periodic and are further investigated in the later sections. The effect of higher-order dynamics is to accentuate the chaotic nature of the dynamics and reveal the broad-band nature of the chaos. The usage of a single strain signal from a wrong location could sometimes give misleading results. 86 560- .............. i ............. i ................ i ................ r. ............... 1_ .............. . 3 : 35 .. .. ................................................................................. . 11140 320... ................................................................................... - E 1, i 8 1 1 1 1 1 0 50 100 150 200 250 300 frequencyinHz Figure 5.9: FFT of the strain signal measured at location 6 on the beam at 12.0 Hz. g I l I I I €40 : 1 3 E 3 £20 .. .. .. .................. l ....................... i ........... .. O .' g 15, . . x g of L i i I 1 1 O 50 100 150 200 250 300 frequencyinl-lz Figure 5.10: FFT of the strain signal measured at location 6 on the beam at 41.1 Hz. 5.3.1 Reconstruction Results From the delay map and broad band nature of the FFT of the strain signal and the contact point dynamics at 12.0 Hz and 41.1 Hz, we suspect the dynamics at 12.0 and 41.1 Hz to be chaotic. To verify the suspicion we carried the reconstruction studies detailed in Chapter 3 on the time series consisting of strain at location 6 on the beam. From the mutual information [36] we realized a time delay of 21 sampling intervals. We used this time delay to determine the embedding dimension using 87 20log1 O(abs(fft(96))) Figure 5.11: F F T of the noise level of the strain signal measured at location 6 on the beam. the method of FNN. For determininng the FNN, we used the same method that we described in Chapter 3. As before we decided that a neighbor is false when the proportional change in the distance (with the additional dimension) is greater than 10. The embedding dimension was found to be five. Thus five is the number of delay coordinates necessary to describe the data. There is not a clear relationship between delay coordinates and dimension of the state. But it suggests the number of active modes to be in the ball park of half the number which is close to three in this case, consistent with the 99.99% criterion. We remind the reader that in the measurement of FNN we used some tolerances and these may have some relation to the criterion of 99.99% “energy”. With the embedding dimension of five and the above stated time delay we realized the positive Lyapunov exponent to be 3.73 (using Wolf’s algorithm [91]), thus assuring us of the chaotic nature of signal. A similar study of the dynamics at 41.1 Hz resulted in a time delay of 8 sampling intervals, an embedding dimension of 8 and a positive Lyapunov exponent of 1.91, assuring us of the chaotic nature of the dynamics at 41.1 Hz. A better algorithm to obtain Lyapunov spectrum is given by Eckmann [25]. Our attempts to obtain this spectrum were unsuccessful.- We leave this and the study of 88 the dependence of the predictive ability of the model on the strength of chaos for our future studies. 5.4 Proper Orthogonal Modes: Giving a Spatial Description We have already discussed the significance of proper orthogonal modes and values visa-121's a numerical system. In this section we describe the spatial content of the beam dynamics in terms of proper orthogonal modes and proper orthogonal values. We pick here two particular cases. The first is the dynamics of beam around 12.0 Hz and the other is around 41.1 Hz. In both these regimes we found the dynamics to be chaotic. The selection has been made as we intend to use the dominant modes in these regimes in building a reduced—order model for the beam in the next chapter. We collected the six strain-signal time series at the operating frequencies 12.0 and 41.1 Hz. Using the method described in the previous chapter we have converted the six strain signals into displacement signals at locations 1 — 6. After subtracting the mean displacement from each of these displacement time series, we arranged them to form a matrix X of six columns. We computed the covariance matrix R = XTX. This matrix is a positive definite matrix and the singular values and vectors of this 6 x 6 matrix, yields the proper orthogonal modes. Using the proper values we identified the dominant modes in the system. As we have the values only at six different discrete locations we need to curve-fit the discrete modes to get continuous modal functions. Though the modes are orthonormal, by doing a curve-fit we lose the orthogonality, as the definition of inner-product of the modes is changed from a finite sum to a definite integral. In order to make the continuous functions orthonormal with reference to the distributed mass matrix, we use the Gram-Schmidt reorthonormalization technique, to yield what we call proper orthogonal modal functions (POMFs). We have used the linear natural modes (LNMs) of the cantilever for the curve fit. The motive is to identify the content of the natural modes in the dominant proper orthogonal modes. tributi of LNMs to POMFS at 12.0 89 Table 5.1:Lor 11L POM # LNM 1 LNM 2 LNM 3 LNM 4 LNM 5 LNM 6 1 0.7113 0.2474 0.0296 0.0089 0.0021 0.0007 2 0.2494 0.7265 0.0108 0.0096 0.0035 0.0002 3 0.1978 0.2343 0.5538 0.0053 0.0051 0.0038 In Figure 5.12 we presented the six proper orthogonal modes and their continuous versions, for the dynamics at 12.0 Hz frequency. The energy content (defined as the ratio of corresponding proper value to the sum of all the proper values) of each of the POMS is indicated. The POMS for which we have not indicated any energy content have zero proper values. For example the first three dominant modes account for fractions of 0.9913, 0.0083 and 0.0003 of the total energy. This implies the dimension of the active state is six. We recall that the false-nearest-neighbors have gone to “zero” when the embedding dimension is five, consistent with the number of active modes. In Figure 5.13 we present the POMFs obtained by Gram-Schmidt orthonormal- ization of the continuous versions of the POMs. The effect of orthonormalization is to take out the components that are already present in the previous vectors. As such it would be wrong to assign the same energy content to the POMFs, though the order of dominance is still maintained. We suspect that we can still determine the energy content in the orthonormalized functions, by using the same ideas that we used to orthonormalize the functions. In Table 5.1 we give the percentage contribution of the natural modes in the POMFs. From this table we realize that the POMFs are not purely the LNMs. In fact, the first five LNMs account for 99.99% energy of the first three POMFs. This indicates that we need the first five natural modes to describe the same dynamics (at least from the 99.99% energy perspective) as that of the three POMF model. Definitely this is an advantage favoring POMFs. 90 2 § 20 .......................................... § 0 j ........ 1 E . _20 L L 1 .2m ; L ; 0 0.1 0.2 0.3 0.4 O 0.1 02 0.3 0.4 'sensor' location 'ssnsor' location Figure 5.12: Proper orthogonal modes of chaotic dynamics at a driving frequency of 12.0 Hz. The continuous versions of the discrete modes are not orthonormalized. In Figures 5.14 and 5.15 we present the POMS and the orthonormalized POMFS obtained from the chaotic dynamics observed at the operating frequency of 41.1 Hz. In Table 5.2 we Show the contribution of LNMS of the cantilever to these POMFS. In this case we need the first five LNMS to account for 99.99% energy in the first four POMFS and six LNMs to account for the same proportion of energy content in the fourth POMF. Thus we need six LN MS to describe the same dynamics as that of the four POMF model, again proving that POMFS give a better reduction than LNMS. 91 6 f 1 10 2 § 4 .......................................... 5 ...................................... To '8 2 ........................................... o ........................................ E o ; i i _5 1 1 i O 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4 10 g s a 0 § 3 ; E ; _1 O n 1 1 _1o ; 1 ; 0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4 10 - '9 1-11 0 8 E -10 L L ‘. _1 o l L l 0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4 'sensor' location 'Sensor’ location Figure 5.13: POMFS (obtained by orthonormalization of the curves fit on POMs) from chaotic dynamics at a driving frequency of 12.0 Hz. From this contribution table we note that the higher LNMS have a greater contri- bution to the dominant POMFS. This indicates that a model built from three POMFS would not be similar to the model built from three LNMS. Also as indicated by the F FT of strain signal (Figure 5.10), we suspected more than six modes to be active in these dynamics. The varying degree of contribution of LNMS to the POMFS gives us an indication of the modal coupling involved. This in turn reflects on the better representation of the system dynamics by POMFS rather than LNMS. 92 a—fi —L 4 5 . i a 8 e .2 i L L _10 i L L 0 0.1 02 0.3 0.4 o 0.1 0.2 0.3 0.4 modal coord 0 0.1 0.2 0.3 0.4 0 0.1 02 0.3 0.4 'sansor' location '30an location Figure 5.14: Proper orthogonal modes of chaotic dynamics at a driving frequency of 41.1 Hz. The continuous versions of the discrete modes are not orthonormalized. 5.4.1 Proper orthogonal modes: using a larger number of sensors Until now we have used only six sensors to determine the POMS. Earlier we have converted the displacement information determined at discrete locations into displace- ment at any point on the beam using the first six LNMS. Thus we have an infinite number of pseudo-sensors. However because we used only the first Six LNMS in the curve fit, we are making an implicit assumption that the system dynamics does not 93 ...L U“ 0 model coord O I 0| 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4 8 o modal coord o '3 O _a O modal coord O '2; O ; ; 3 40 1 3 3 0.1 02 0.3 0.4 0 0.1 0.2 0.3 0.4 '50an location 'Sensor' location Figure 5.15: POMFs (obtained by orthonormalization of the curves fit on POMS) of chaotic dynamics at a driving frequency of 41.1 Hz. contain any more than the first six LNMS of the cantilever. Thus if we use any more than six pseudo-sensors in the displacement measurement, the POD would result in a correlation matrix of rank no more than six and hence we can at most have Six non-zero proper values. However the resulting (discrete) POMS could differ. In this section we see the effect of using a larger number of “sensors”. Using the pseudo-sensors we obtained the displacement at points uniformly dis- tributed in steps of 0.01 m all across the beam of length 0.40m. We used this infor- ‘ Table 5.2: Go. 94 tribution of LNMS to POM MAL POM # LNM 1 LNM 2 LNM 3 LNM 4 LNM 5 LNM 6 1 0.2872 0.5748 0.1081 0.0216 0.0060 0.0023 2 0.7764 0.0443 0.1386 0.0326 0.0045 0.0038 3 0.1920 0.0937 0.5582 0.1386 0.0145 0.0030 4 0.1179 0.0819 0.0323 0.6254 0.1093 0.0331 mation and built the covariance matrix and obtained the proper orthogonal modes and modal functions. In Figures 5.16 and 5.17 we present the same. On comparing these modal functions with the ones obtained using the actual sensors we realize that the continuous versions of the dominant modes tend to the orthonormalized functions with increasing number of sensors. This indicates larger number of sensors gives a better representation of the system. However this would render the process more expensive. Thus with an. “optimal” number of sensors we can achieve the POMFS similar to the ones using larger number of pseudo-sensors, using the Gram-Schmidt orthonormalization technique. 5.5 Discussion In this and earlier chapters we presented the outline of the experiment, discussed the dynamics of the system under frictional excitation and described the method of extracting the POMFS from the experimental data. In this process we have seen the ease with which one can obtain the POMS from an experiment compared to the difficulty in obtaining the LNMs from sophisticated and expensive software with a large learning curve. For the systems that we considered we showed that the number of LN Ms needed to span the space of the dominant POMFS is larger than the number of dominant modes. This indicates that we need more LNMS than the POMS to describe the same system dynamics. We will verify these issues later in Chapter 6 with models built of both POMS and LNMS. 95 ..... ....... ....... P N modal coord E a 8 E . . -0.5 A e .05 . 1 o 0.1 0.2 0.3 0.4 0.5 7 : 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4 'sensor' location 'sensor location Figure 5.16: Proper orthogonal modes of chaotic dynamics at a driving frequency of 12.0 Hz, using 41 pseudo-sensors. The continuous versions of the discrete modes are not orthonormalized. 96 0.3 0.4 2 0 'sansor’loca' n 0.1 0 02 0.3 0.4 'sensor' location 1 0 O 9000 .009: Econ .03.: ing frequency of 12.0 Hz, using 41 POMF S of chaotic dynamics at a driv pseudo-sensors. .17: 5 Figure Chapter 6 Mathematical Model, Reduction and Validation 6.1 Mathematical Model The system under consideration (described in the previous chapter) consists of a cantilever with an attachment for convenient variation of normal load. Because of the attachment, we have the friction force F, acting not on the cantilever but at the end of an extension. This force can be transferred to the cantilever tip, by an equivalent force and moment combination (the force is F and the moment is F I , where I is the length of the extension). To take the inertial effects of the extension and the fixture we lump the masses together at the cantilever tip. Hence we need the beam equation where both external forces and moments are present. Assuming that there is no shear, we use the extended Hamilton’s principle [58] to get the standard PDE for a Euler-Bernoulli beam and write our system equation (zero damping is assumed here) as 090; t) 67:1” )62y(:, t)]+ P3( )6_:___ya(t21t )2 P5” t)— where y(z,t) 18613118 deflection at the point a: at time t. E(:I:),I(:ca) and p(z) are the (6.1) elastic modulus, the area moment of inertia and the linear density respectively of the system at location (1:. p(:r, t) and q(:c, t) are the external forces and moments applied 97 98 at location a: on the beam at time t. Note that p and q are implicitly time varying functions of relative velocity, normal load and other parameters in case of frictional excitation. The parameter values are given in Chapter 4. With this assumption of “an Euler-Bernoulli cantilever beam with a lumped mass at the free end” for the system, we numerically estimated frequencies of 3.13, 21.97, 64.72, 130.64, 220.55 and 334.12 Hz. When we compare these frequencies with the experimental ones presented in the Chapter 4 ( 3.75, 18.4, 53.75, 102.5, 172.5 and 265 Hz), we notice that the estimated frequencies are higher, with the exception of the fundamental frequency. Hence we expect to see a shift in the predicted frequency sweeps to the right. With a Timoshenko beam model we could get better astimates of the frequencies [43], however this would increase the order of the system from two to four in the time domain. We defer the Timoshenko model to our future studies. The usual practice in the literature for solving these partial differential equations is to first convert the PDEs into ODEs by using the Galerkin’s approach. In this method, we write the deflection as a weighted sum of functions {so.-(x)w3‘”’(x)dz (6.4) M..- = :pr(x)cps(x)dz + mL‘F,;)‘11>‘K.1U + (Q‘F.1)"1a(t) Now we can substitute for f (t) in (6.9) to get U = K.:U+F.zf(t) = [1 - (FF.:)"F.1‘I’*] K..U + (rm-11111.6) (6.10) (6.11) (6.12) (6.13) (6.14) (6.15) Denoting the coefficients of U and a(t) by C and D, we can write the above equation as 1'} = CU + Da(t) (6.16) 102 Now using (6.10), we can define a new state vector W which is one dimension less than U and is given by . 0 . I2N—l 2N—1 E U = 0 . .x. 0 'WP'] . . . “spar-1F.) W + 0 v(t) ‘PN(3.) «061(1") 1 - ¢N(3.) - Substituting this relation for U in (6.16) results in cw + 116(1) = C(GW + Hv(t)) + Da(t) (6.17) Rearranging the above terms and multiplying the equation by (GTG)'1GT results in w = (GTG)‘1GTCGW + (GTG)‘1GTCHv(t)+ (GTG)-1GT(D - H)a(t) (6.18) By denoting the coefficients of W, v(t) and a(t) by P, Q, R respectively we have w = PW + Qv(t) + 116(1) (6.19) Thus we obtain the state equation for the system in the stick mode. The state is one dimension less than that of the slip mode. This accounts for the temporary constraint associated with the stick. For the sake of consistency/ simplicity in the im- plementation, we wish to use the redundant state U instead of W for the stuck mode also. Hence we use the state equation (6.16) which (for the sake of easy reference) we write as, U = K..U + Fma(t) (6.20) where K“ = C and Feta = D. Now we know the state equations for the system in both the modes: equation (6.9) is for slip and equation (6.20) is for stick. We now need the criteria for the slip-stick transitions. These criteria are based on the relative velocity between the slider and the beam and the elastic force Fy in the beam. When the relative velocity 103 of beam at the point of application 3" is zero and the resultant of the components of Fy and Nz along the direction of sliding is bounded by the limits of static friction [-psN,, +psNz], the system sticks and in all other cases the system slips. Equation (6.15) can be used to obtain Fy. Notice that the usage of equation (6.15) is based on the assumption that the beam were sticking, otherwise the acceleration of the beam would be different from the slider. 6.1.2 Finite contact stiffness with Coulomb friction excita- tion Till now we have assumed the contact stiffness to be infinite. In order to use a finite contact-stiffness model we need to modify the state equations and the criteria for slipping and sticking. The state equation for the slip mode remains the same since f (t) is defined. From a free-body diagram we can easily verify that f (t) is same as the frictional excitation at the massless contact point. Now we determine the frictional excitation during sliding. Let us denote the contact stiffness by Kc. We denote the displacement of the tip by YT(t) and the displacement of contact point as Yc(t). Then as seen in Chapter 4, we have f(t) = KC(YT(t) — Yc(t)). Upon differentiation we have f (t) = KC(YT(t) — Yc(t)). We have assumed Coulomb-friction model. This implies that the friction force f (t) = :1:pr 3 during sliding, the sign decided by the direction of relative motion at the contact. Hence Yc(t) = YT(t) 2]: pKNz(t)/Kc. We note that N z(t) is a function of the contact point displacement and is given by Nz(t) = Nzo + NZmYc(t), where N20 and NZ", are the nominal normal-load and the slope of normal load variation, the measurements of which were described in __£c__ uKNZm+Kc values of Kc,px and NZ", we realize that Yc(t) z YT(t) during sliding. Thus the first criterion for sticking is Yc(t) — v(t) z YT(t) - v(t) = 0. Chapter 4. Using this relation we have Yc(t) = YT(t). Comparing the Earlier for infinite stiffness we compared the effective force on the beam with bounds on the static friction. In that case we used equation (6.15) to get the effective 104 restraining force. When the elements of contact were stuck, we computed the effective force using the acceleration of the slider (which would have been acceleration of the beam tip, had the slider stuck to the beam tip in the infinite stiffness situation). We cannot do the same in the finite stiffness case. With finite stiffness during the stuck state the acceleration of the beam tip is different from the slider acceleration. - However we have a simpler method of determining the force in this case using the finite stiffness Kc. If the beam were to stick, due to the finite contact stiffness, we have a relative displacement between the actual contact point and the tip of the beam. Here Yc(t) = do + d(t), where do is the point on the slider where the contact is made and d(t) is the slider displacement. In case of infinite contact stiffness, Yc(t) = YT(t). However with a finite stiffness Kc, we have f(t) = Kc(Yc(t) - YT(t)), where f(t) is force from the beam. Thus we can replace equation (6.15) with this relation and proceed. Thus we have f(t) = Kc(Yc(t) — a(t)) N = Kc((do + d(t)) — Z ¢.(x')U.-(t))) i=1 = Kc((do + 61(1)) - W‘U) (6-21) where W'=[¢.(z') Mr) 0 0] (6.22) In the above force relation we have the effective force in terms of the state and input displacement which are all known. Thus we know the effective restraining force, which can be compared with the static friction bounds to determine whether the contact elements stick. We can use the force information defined above in the equation (6.9) to obtain the equation for the stuck mode. Thus U K..U — F.1Kc((do + a(t)) - w'U) = (K,( - KcF,(‘I")U + KcF,1(do + (f(t)) (6.23) 105 For easy reference we will rewrite the above equation as U = KmU + 1?...(61o + d(t)) (6.24) where the definitions of Km and Fm are obvious. We now know the state equations for the system with non-zero contact compliance in both the modes (slip is represented by equation (6.9) and stick is represented by equation (6.24)). As before the criteria for slipping and sticking are based on the relative velocity between the slider and the beam and the dynamics of the beam. When the relative velocity of beam at the point of application :6“ is zero and the resultant components of Fy given by equation (6.21) and Nz along the direction of sliding is bounded by the limits of static friction [—p5Nz, +p5Nz], the system sticks. In all other cases the system slips. Equation (6.21) can be used to obtain Fy to determine the sticking criterion. 6.1.3 Effects of normal vibrations of the beam on normal load In the previous sections we have not considered the effects of the normal (i.e. in the Z direction) vibrations of the loading beam and the cantilever. We realize that the coupling of normal, torsional and lateral degrees of freedom may be important in generation of interesting friction-induced dynamics. We have not considered the torsional coupling effects as the geometry of the rubber nub and the torsion of beam are not well known and any assumptions would be arbitrary. The torsional and normal effects could be small but could have a significant impact on the normal load variation. To this end in the following we consider the effect of normal vibrations. We suspect the source of normal vibrations to be the slope of the slider and the effect to be significant in the neighborhood of the fundamental frequency of the normal vibrations (especially when the slope is large) and it could contribute significantly to the bifurcations in that neighborhood. In order to account for this, we used a lumped 106 mass model. We assumed an effective mass mL to be lumped at the Cantilever tip. Though the stiffness of the cantilever in the Z-direction is large compared to the stiffness in Y-direction, it is not 00 (the first frequency of vibration in Z direction is around 70 Hz). We therefore do not neglect the contribution of the cantilever to the normal vibrations. In this case we have the lumped mass m), connected by two springs of stiffnesses kb and Ice. kg, and ke are given by (IE'I)(,z/L3 and (EI),z/I3 respectively. (E1 )b2 and (E1)eZ are the flexural rigidities of the cantilever and extension beam in Z direction respectively. L and I are the lengths of the cantilever and extension beams. The second end of spring k5 is fixed and spring 1:, rides on the slider surface. In Figure 6.1 we show a schematic of this simplified model for normal vibrations. As depicted, Z1, and Ze are the displacements of the cantilever tip and the extension 3 k.= 3(EI)bz/L Zbr mb kc: 3(EI)cZ/13 v Z6]— Figure 6.1: Simplified model for normal vibrations tip in the Z direction. We note that the extension tip is riding on the slider surface and hence is given by Z, = ZmYT + Z, where Zm, Y1 and Zc are the slope of slider, extension tip displacement in the Y direction, and nominal position of the extension tip in the Z direction. We can see that Z, and Ze are connected by mi. + (c. + c.)z'. + (k. + k.)Z.. = 1.2. + c.z'., (6.25) where q, and 6e are damping coefficients added to account for damping in the can- tilever and the rubber nub (both of which are assumed proportional viscous). The 107 damping is not represented in the schematic. The normal load is then given by f~ = k.(Ze - Z.) + c.(Z. — Z'.) + 16.2... (6.26) The normal load can thus be computed and be used in the friction force calculation. We will illustrate the effects of the normal vibrations later in the sections on frequency sweeps. 6.1.4 Modal reduction using POMs In the previous sections we have given a mathematical model for the system under consideration. The information that was not defined was the 905(3). We propose here the usage of proper orthogonal modes for (x)vm.-dx k=1 m=1 n n z=L . = —:vavmi/.=.¢1<2>¢2”>(x)dx, 6.29) k=l m=l 3=L M).- = mL¢j(L)¢a(L)+P/z=o ¢j($)¢e($)d$ = mam/«(m + ,, f: .2} mm.- i ¢m(x)vm.-dx m=1 -_- mL¢,(L)z,b,-(L) + p i in: ”kjvmi/;:L ¢k(x)¢m(x)dz (6.30) Ic=l m=l P.- = I¢;(L)+¢.(L)-I¢.(L) (6.31) Note that ¢,-(a:), are the LNMs of the cantilever section, and due to their orthog- onality, we have simpler forms for the undetermined definite integrals in the above equations. Specifically we note that [L0 ¢j(:r)¢,-(x)dx = 6jg/2 (6.32) I: and j; ¢.(x)¢:i"(z)dx = (666/2 (6.33) 109 Using equations (6.32) and (6.33), in equations (6.29), (6.30) and (6.31), we get for each j,i 6 [1,N], KJ',‘ = —(E1) 2“: i vkjvmgflz/26km = -(EI)1 i ‘Umjvmgflfn/2 (6.34) It=l m=1 m=l M).- = mL¢,(L)¢,(L) + P f: "WM/2 (5535) 16:1 Pi = I1113-(14) + WU!) - MAL) (6-35) The M,K and P matrices thus obtained can now be used in getting the matri- ces required in equations (6.9) and (6.20). Now we have a way to implement the conversion of PDE into a set of ODEs using experimentally obtained POMS. 6.2 Numerical Solution Technique The numerical simulation technique adopted is much similar to the one that we used in Chapter 3. The technique involved the usage of IMSL routine DIVPRK (a Runge- Kutta ODE integrator in double precision). We used a variable time step integration to integrate the slip and stick equations. For a given set of initial conditions we determine the relative velocity and the excitation force. These we use to define a flag for slip or stick. Then we proceed to integrate in a while loop till we reach the specified time of simulation. At each time step we determine the relative velocity and effective force on the beam and compare their values for slip and stick criteria. In case of a cross-over (relative velocity crossing zero or the effective force crossing the static friction force bounds) we step back and use a simple interpolation, to determine a smaller time step and integrate. This we continue until the cross-over is within a predefined tolerance zone. At this point after verifying the respective criterion we move to the following state. In this way we compute the state U at each time step. In turn we use the state U to compute the displacement y(a:,t), velocity g(z,t) or strain e(:c,t) at any given point x on the beam and given time t using the following 110 equations N yet) = Emma-(t) (6.37) i=1 N (KM) = 221(z)U.+~(t) (6.38) I ~ ,. .(.-,t) = 52¢.(2)U;(t) (6.39) (6.40) where h is the beam thickness. 6.3 Reduced Model — Numerical Simulations 6.3.1 Friction-force characteristics — different scenarios In previous sections we described the effects of contact compliance and the transverse inertial effects. In this section we will present the resulting force characteristics, i.e., Fy versus the relative displacement D, and relative velocity V, between the tip of the beam and the slider. In the next section we will present the corresponding frequency and amplitude sweeps. If the contact stiffness is infinite, we do not expect any variation in the displace- ment when there is a direction reversal and we should see a sudden change only in the sign of kinetic friction (in Coulomb’s friction model). This is seen in Figure 6.2. Similar is the case where the static friction coefficient is slightly higher than the ki- netic friction coefficient. The only difference is the beam tip now gets moved to a larger distance in the stick regime. The slope Zm of the shaker axis is present though not perceptible as the displacement is too small. The slope Z,,, is reflected in the top and bottom lines (in the Force-displacement characteristic) when the displacements are higher and is seen as a wedge shape trapezium. In the case of compliant contact model as explained in an earlier chapter, we 111 0.2 . ‘ 0.2 f a 0., ......... ............. ............. 0., ........... .................... ........... g 0 .................................... 0 .......................................... LI. ‘0“ ......................... r. .......... _O1) ........................................... -o.2 3 .0. -1 O 1 2 3.1 -005 0 005 01 x10 0_2..... ....... ............ ........... .. 0.2. .......... .......... _.,m,,..: .......... g o ...... ] ................................... o) .............................. ........... .1 2 _0.2 ......... ............. ............. _0.2. ......... .......... .......... 4 -1 O 1 2 -O.1 -0.05 0 0.05 O 1 '0 1110“ 'V Figure 6.2: Friction characteristics of the infinite contact stiffness model. The top row shows the case where p5 = ,uK and the bottom row shows the case where #8 > [Ix should expect the lines of direction reversal (in the friction-displacement plot) to be slanted and the slope should correspond to the contact stiffness. The finite stiffness results in a non-zero displacement as the the relative velocity changes sign. This would also result in a hysteresis loop in the friction-velocity plot. The friction-velocity plots need a more detailed description in case of compliant contact friction model. This is because now we have the so-called micro-stick. In case of “pure” slip motion, the hysteresis loop becomes more obvious. On direction reversal we have the relative velocity going to zero and hence a micro-stick regime is present which is reflected in the slanted curves of the friction-velocity plots. Notice that the micro-stick event begins at zero relative velocity and ends at a non-zero relative velocity. In the case of stick-slip motion (the top row of Figures 6.3 and 6.4), the sticks are reflected as 112 0.2 . - 0.1 ..... ............. .................................... Friction o 4.1...” ......... ............. i ............. -0.2 0.2 . . . . . . 01.... ...................................... , . . _01... ....................................... . n I . . —0.2 Figure 6.3: Friction characteristics of the compliant contact model with p5 = pg. The top row shows the case of macroscopic stick-slip and the bottom row shows the case of micro-stick loops in the friction-velocity plots, similar to the observations made by Liang [52] and Liang and Feeny [53]. The case of p5 > pg is similar to the one where they are both equal, the only difference being that the sticking continues until a larger force is exerted for a longer time. The effect of transverse inertial mass is not obvious for the low values of mass and slope of the shaker axis under consideration. However for larger slopes, though the friction force characteristic remains same, the effective force Fy characteristic changes considerably. 113 .................... 0.2 ............. ............. .............. 0.2 ........... ......... - V a ?XT‘O)LKOD) . -4 .................... Friction o o 0.] go _02 ............. f ............. f ............. _0.2 ........... ..g .................. .4 -5 0 5 1O -01 -0.05 0 005 01 .4 x 10 0.2.... .. ....... ............. ............. 4 . ,5 E o ...................................... u. _0.2... . ....... ............ .............. -2 0 2 4 'Figure 6.4: Friction characteristics of the compliant contact model with p5 > p g. The top row shows the case of stick-slip and the bottom row shows the case of micro-stick 6.3.2 Frequency sweeps In the friction-characteristics presented above we obtained the relative displacement and velocity between the beam tip and slider velocity using the reduced model of the beam. The model reduction was done using the POMFS at 12.0 Hz. In this section we present the frequency sweeps of the reduced model. We later use these frequency sweeps to compare the reduced models obtained using the LNMs of the cantilever beam and POMs from the experiment. We carried out frequency and amplitude sweeps using the reduced model that was 114 described in the earlier sections. Earlier in Chapter 5 we used a 10 degree polynomial fit to define the relation between the shaker frequency and shaker amplitude. We used the same polynomial in defining the slider amplitude at a given frequency during the numerical simulation of frequency sweep. Each sweep was carried for 600 cycles of which only the last 100 cycles were stored. Thus we are discarding 500 cycles to account for transients. Incidentally we used a one-percent proportional damping and the two-percent settling time is well within the linear damping limits. We however note that the damping effects could be quite different in this nonlinear system. In each of the excitation cycles we measured the strain at location 6 on the beam, whenever v(t) crosses zero with a positive slope. Thus we computed 100 strain values at each frequency step when the slider acceleration crossed its peak. We swept from 5 - 60 Hz in steps of 0.1. We caution the reader that the above plots are not connected to the amplitude-frequency response curves. These plots are sweeps of Poincare sections along the curve described by the 10 degree polynomial in the configuration space. 6.3.3 Different levels of reduction using POMs In this section we present the frequency sweeps for different levels of reduction. For comparison between different reduced models we used the Coulomb friction model with finite contact stiffness. In Figure 6.5 we show the sweeps obtained from the POMFS determined at 12.0 Hz in the physical experiment superimposed on the ex- perimental frequency sweep. From these plots we see that the most dominant mode alone is not able to pre- dict the system’s behavior. With two dominant modes we are able to retrieve the resonant dynamics around 12 Hz range (the POMs are extracted from the chaotic orbits at this frequency), but the high frequency dynamics are not noticed. How- ever with three POMs we are able to see the higher resonant frequencies, though the instabilities (indicated by the smear of points) are not as pronounced as in the ex- periment. These instabilities become more prominent as we include more POMs but 115 .a q.p. st. at be 6 O l _. q.p. 31. at Ice 6 o d l d q.p. st. at loo 6 30 40 50 slider ire on In Hz slider iv an In Hz qu cy —- Experiment, — POMF equ cy 10 20 30 40 50 10 Figure 6.5: “Frequency sweeps” of reduced models obtained at different levels of reduction using POMFS (extracted from experiment at 12.0 Hz and the continuous versions orthonormalized) superimposed on the experimental sweep. the sweep characteristics (represented by the resonant frequencies) remain same. The inclusion of higher modes results in higher computing time and grows exponentially with increasing number of modes. Also the inclusion of higher POMFS could result in deviations from the experiment as the signal noise is generally reflected in higher POMFS whose corresponding proper—values are close to zero. On comparing these predicted sweeps with the ones from the experiment we realize that though the stick-slip instabilities around 12 and 41 Hz of the physical experiment do occur the latter instablities occur at slightly higher frequencies and the resonant levels match well within five percent (the 3-POMF model). The dynamics in the unstable regions in the predicted sweeps seems chaotic but is not as prominent as 116 in the experiment, especially at the higher frequencies. However we note that the periodic orbits and also the zones of pure stick and stick slip are well predicted by the reduced model. We expected the shift from the linear behavior of the numerical model. The re- duction in order though introduces new constraints and in a way stiffens the system, it is not truly the cause for the shift in higher frequency dynamics to higher frequen- cies. This we verified by carrying a frequency sweep using a numerical model with ten LN Ms. The addition of more modes did not significantly improve the dynamic characteristics with reference to the frequency shift. Thus we attribute the shift to the fundamental modeling assumption of “lumped mass at the free end of a Euler- Bernoulli cantilever beam”. We leave the usage of Timoshenko beam model (which is fourth order in time) to our future work. In order to account for the uncertainities involved in our modeling efforts, viz., beam model, damping, friction characteristics, etc. we considered a model similar to that of the POMF model, but with a large number of LNMs (ten) as a “truth set”. We compared the frequency sweeps of reduced models with that of this “truth set” (see Figure 6.6). This comparison reveals that the resonant frequencies predicted by the POMF models are almost same as the ones predicted by the “truth set” model. The convergence of the reduced model frequency sweep with that of the “truth set” is acheived rather earlier just with three POMFs! We suspect that the the contamination of signal noise on higher POMFs is reflected as the deviations in the six POMF model. In this particular system we were fortunate to have the LNMs of the cantilever readily available. As such we built reduced models with the LNMS and carried fre- quency sweeps similar to those of POMF models and in Figure 6.7 we present the same. From Figures 6.6 and 6.7 we realize that in this particular example POMFs have faster convergence than the LNMs. In Figure 6.8 we show the sweeps obtained from the POMFS determined at 41.1 Hz in the physical experiment superimposed on the experimental sweep. 117 q.p.statlocG o d I d ‘ O .8 a «g o a e -1. 1‘0 20 So 40 So 10 20 ab 43 So slider Ir an In Hz slider Ir In Hz equ Cy --10LNM, -POMF ' ’ Figure 6.6: “Frequency sweeps” of reduced models obtained at different levels of reduction using POMFS (extracted from experiment at 12.0 Hz and the continuous versions orthonormalized) superimposed on that of the lO—LNM model (“truth set”). We recall that four of this set of POMFS were dominant (whose combined energy content is at least 99.99%). Thus we expect to reveal the features of the beam dynam- ics with at least four POMS. This is what we observe in the above plots. The models built by using less than four POMFs from the 41.1 Hz dynamics did not reveal the higher resonant frequencies. Also the dynamics in the 12 Hz neighborhood predicted by these models deviated from that of the experimental sweeps. However with four POMF S model the deviation is reduced and also the higher resonant dynamics are revealed. For a better appreciation of the three POMF reduced model (from the 12.0 Hz chaotic dynamics in the experiment) we show the phase portraits (displacement versus 118 a A 1020304050 1020304050 slider hemency In Hz sllderIrequency in Hz ——10LNM,--LNMS Figure 6.7: “Frequency Sweeps” of reduced models obtained at different levels of reduction using LNMS of the cantilever superimposed on that of the 10-LNM model (“truth set”). velocity of the contact point on the beam) in Figure 6.9. These are obtained at 11, 11.9 and 42.9 Hz respectively of the resulting dynamics in the sweeps. On comparison with the experimental phase-portraits (presented in Chapter 5), the qualitative nature and the order of the size of the attractor are found to be well preserved. We also note that the meandering (of the beam tip on the oscillating slider) seen in the experiment could not be effectively realized with this this reduced model. 119 q.p. st. at Ice 6 o ..A l _a _a l —A q.p. st. at loo 6 O —A q.p. st. at Ioc6 O l _a 10 20 30 40 50 10 20 30 40 50 slider fr uen In Hz slider Ir en in Hz eq cy —— Experiment, —- POMF equ cy Figure 6.8: “Frequency sweeps” of reduced models obtained at different levels of re- duction - from most dominant mode to all the available POMFS at 41.1 Hz (compliant contact model). 6.3.4 Finite and infinite contact-stiffness friction excitation In Figure 6.10 we present the plots of frequency sweeps for finite and infinite contact- stiffness model for friction. Both the sweeps were generated by the 3 POM reduced model. Also between the two friction models we note that the compliance model gives a better prediction, in the sense that the occurence of instabilities and resonances are more pronounced in the latter. However both the models Show the shift of stick—Slip instabilities (that are observed around 12 and 40 Hz in the experiment). We also note 120 0.1 0.1 E 0.05 ........................ S. g o .......................... o 2 %_0.05 ........................ -0.1 . ., —o.1 ‘ .0." I -0.002 0 0.002 -0‘.005 0 0.005 -0.001 0 0.001 tip displacement in m tip displacement in m tip displacement in m Figure 6.9: Phase-portraits (contact point displacement versus velocity) from the three-POMF reduced model at 11, 11.9 and 42.9 Hz. that with the infinite stiffness model, we see a better smearing of the quarter-phase strain in the 12.0 Hz neighborhood, indicating the non-periodic nature of the response in this neighborhood. This implies that the infinite contact stiffness friction model is better (compared to finite contact stiffness friction model) in the low frequency region. This leaves us with a suspicion that the friction characteristic is not a simple Coulomb model with a finite contact-stiffness. A possible model could comprise of varying contact stiffness (nonlinear). For example we can replace spring Kc with an hardening Spring such that, the stiffness increases with increase in the amplitude or a decrease in the frequency in these “frequency sweeps”. With this insight we leave more complex friction modeling to future studies. 6.3.5 Normal vibrations included In an earlier section we mentioned that the coupling between torsional, normal and lateral degrees of freedom could significantly contribute to the friction-induced dy- namics of the beam. We described the effects of the normal vibrations caused by the tilt in Slider axis. Here we present a numerical frequency sweep from a reduced model where we accounted for the normal vibrations. For these simulations we lumped one 121 quarter-phase strain at loo 6 quarter-phase strain at Ice 6 5 10 15 20 25 30 35 4O 45 50 55 sliderirequencyInI-lz Figure 6.10: “Frequency sweeps” of the system using the reduced model from POMFS (of the three dominant POMS obtained at 12.0 Hz). The top row corresponds to the zero contact compliance friction excitation and the bottom row corresponds to non- zero contact-compliance. third of the beam mass at the end. The parameters like the lumping factor and damping due to the rubber nub were not determined from the experiment and definitely needs a little fine tuning before for proper predictions. We however note that the bifurcation become more prominent in the 40 Hz region. With this insight and note on the normal vibrations, we leave the studies to future work where we intend to account for higher modes in the normal direction and also better lumping of the masses. 122 S3 m... I .a quenen1uuumHNHnnSkazo 6.. 5 10 15 20 25 30 35 40 45 50 55 afiderflequencyhiHh Figure 6.11: “Frequency sweep” of the system using the reduced model from POMFS (of the three dominant POMs) obtained at 12.0 Hz. Normal vibrations are accounted for in the model. 6.3.6 Amplitude sweeps We already have noted that the unstable zones have moved to higher frequencies and the instabilities are not as prominent as in the experiment. In order to confirm the robustness of the reduced model to parameter changes, we carried our sweeps with amplitude, normal load, static-friction coefficient and slider Slope as sweep param- eters. In this section we present the amplitude sweep. The amplitude sweeps were carried at frequencies of 12.5 and 46.9 Hz rather than 12.0 and 41.1 Hz the frequencies where we conducted the physical experiments, to account for the frequency shift in Euler-Bernoulli model. From the amplitude sweep we note that as the amplitude increases we see an increase in the quarter—phase strain value initially followed by a gradual decrease to zero. By neglecting the higher mode effects we can say that the large strain values indicate larger amplitudes of the beam tip and smaller values indicate small oscillations of the beam tip. At both the frequencies the beam undergoes a steady stick in the low amplitude levels and as the amplitude increases the stick-slip instabilities begin to Show. AS the amplitude increase further the beam tip gets into a pure Slip motion. The pure slip motion is more pronounced at higher frequency. 123 6‘: .. n: ..l .. .1 -( d -i q .5 O O quarter-phase strain at loo 6 01 l om 0.5 1 1 .5 2 2.5 3 3.5 4 4.5 5 1.6 I N I i0 an quarter-phase strain at Ice 6 l 01 I: 0 see i— h p p h - )— b .- 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Slider frequency in Hz x 104 Figure 6.12: Amplitude sweeps of the system using the reduced model from POMF S (of the three dominant POMS obtained at 12.0 Hz). The top row corresponds to the sweep at 12.5 Hz and the bottom row corresponds to 46.9 Hz. The other sweeps involving Slider slope, static friction coefficient and normal load did not reveal any interesting dynamics. As such we with hold their presentation. 6.3.7 The effect of static friction coefficient The static friction coefficient p5 that we used in earlier sections was not experime- natlly determined. We estimated p5 = 1.36p g apriori using sweeps on ratio of static to kinetic friction coefficients, by keeping other parameters constant. We varied the ratio ps/pg from 1.0 to 2.5. We conducted these sweeps at 12.5 and 46.9 Hz, and 124 plotted the quarter-phase strain at location 6 versus this ratio. We presented these plots in Figure 6.13. .s 0| d I . . . . . . . . . . . . . . . . . . . . . . . . . . ..................... .° 01 I quarter—phase strain at be 6 O d 10 ox a i .5 0| quarter-phase strain at loc 6 I i N Figure 6.13: p5 / pg sweeps of the system using the reduced model from POMFS (of the three dominant POMS obtained at 12.0 Hz). The top row corresponds to the sweep at 12.5 Hz and the bottom row corresponds to 43.9 Hz. From these plots we note that the dynamics are nonperiodic in the neighborhood of 1.36 for both the frequencies that we considered. This in a way indicates that p5 is in the ball park of 1.36p g. We later confirmed p5 to be 1.38pg. The higher ratios were not considered in our experiments as we found them to be unrealistic. 125 6.4 Conclusions In this chapter we gave a method for reducing the order of the system under con- sideration and validated the same using frequency sweeps. Though we found some discrepancies between the physical and numerical experiments, we attribute these discrepancies to the original mathematical model of the system, like the assumptions of Euler-Bernoulli model for the cantilever, lumped mass for the loading beam and friction law to be a Simple Coulomb with contact compliance characteristic. Some of these discrepencies have been verified by comparing the natural frequencies de- termined from the physical experiment and numerical estimates. Thus for better predictive ability of the reduced model we need a better mathematical model of the system. We presented two sets of models obtained from the chaotic dynamics at two dif- ferent regions. In the first set (at 12.0 Hz) we found three POMFS to capture 99.99% energy whereas in the second set (at 41.1 Hz) we found four POMFS to capture the same. AS expected we noted that the 99.99% POMFS broadly revealed the dynamics of interest. As the three POMFS model is computationally cheaper we prescribe its usage compared to the four POMFS. We considered the case of a Single POMF dominating a periodic orbit. The model built from this single POMF did not reveal the dynamics clearly and the frequency sweep resembled the sweep obtained from using one POMF of the chaotic dynamics at 12.0 Hz. From the comparison between the LNMS and POMFS we realized that the POMF models converge to both the experimental sweep and the sweep from the “truth set”, rather fast with just three POMFS. Also the experimental nature of the determination of POMFS makes it more advantageous as for complex systems the LNMS may not be readily available. However one should be cautious of the effects of contamination from Signal noise. This noise is usually reflected in higher POMFS and so one should 126 not be tempted to use all the POMFS that are available. With better instrumentation, modeling (usage of Timoshenko beam model) and chaotic dynamics with fewer modes we can possibly obtain a two POMFS (if not one POMF) model that shows us all the dynamics of interest. With this note we conclude this final chapter. Chapter 7 Conclusions 7 .1 Summary of the Results and Conclusions In this study, the reduced order modeling of frictionally excited systems has been con- sidered. The studies were both numerical and experimental. The problem of friction induced vibrations was chosen as it provides us with an insight into understanding the problems of squeak and ”squeal. Our aim was to identify the active modes and degrees of freedom from experiments and the coupling information between differ- ent modes which is critical to understanding noise generation and propagation. In this process we have Shown the ease with which one can obtain the POMS and use them effectively in modal reduction. We also have shown that better reduction is achieved by using POMS compared to using LNMS. This is a major plus point fa- voring POMS as in control applications better reduction means lesser cost per run. The main disadvantage of this method lies in the initial cost involved in obtaining the displacement information at different locations. However this disadvantage pales, when we consider the computational expense and the large learning curve associated with the sophisticated and expensive software involved in determining the LNMS. First we explored a numerical system comprising 10 degrees of freedom. We studied the dynamics of this numerical system under frictional excitation using Sim- 127 128 ulations and dynamical system tools base on phase space reconstruction. Later we used proper othogonal decomposition modes in characterizing the Spatial content of the dynamics. We then built a reduced order model using these POMS (the number of dominant modes is six) and studied the dynamics of the reduced model. We validated the reduced model by comparing the predictions of the reduced model with the full model. In these comparisons we used phase portraits, FFTS, statistical measures, occurrence of bifurcations, Lyapunov exponents and skeletons of the chaotic orbits. Though the reduction in the state dimension is only from 20 to 12, the comparisons were quite promising which led us to apply the ideas to a physical system. In applying the ideas to a physical system we considered a cantilever beam sub- jected to frictional excitation. This seemingly simple system is found to be complex enough considering the limited resources and controlled nature of the experiment. As we needed displacement information at different locations simultaneously we had to consider a cost-effective way of obtaining this information. In this process we made use of the simple bending theory in finding a relation between strains and displacement. AS we were limited to six strain-gage signal conditioners, we made use of six linear natural modes of a cantilever in estimating the displacements. We validated this method using Laser measurements and found the method to be quite effective in the low frequency regions, the estimated displacements were about 99% of the laser measurements. However near 41 Hz, the method was not found to be as effective and the estimates were about 80% of the direct measurements. The second aspect of this experiment was in determining the friction characteristic. We determined the characteristic to be Similar to Coulomb’s model but with a contact compliance. Though we realized the difference between the static and kinetic friction coefficients between the contact elements, we had used a static friction coefficient p5 = 1.36p g. We later experimentally determined the mean static friction coefficient p5 = 1.38pg with a standard deviation of 3.7% of the mean. 129 We assumed a Euler-Bernoulli cantilever (with a lumped mass at the free end) model for our system. We obtained the discrete POMS from the Spatial covariance matrix and curve fit them using linear natural modes of the cantilever, to get the proper orthogonal modal functions. We then built reduced—order models using differ- ent sets of POMFS in a Galerkin’s approximation. With the help of frequency sweeps we validated the reduced-order model. By taking into account the discrepancies that are introduced due to Euler-Bernoulli beam model assumption for the system rather than Timoshenko beam model (where shear effects are considered not negligible) we note that the reduced model’s predic- tive ability is close, as far as zones of stick, slip and stick-Slip and phase reversals are considered. However the chaotic zones and bifurcations are not pronounced and we suspect the neglect of the effects of shear, moment of inertia, normal and torsional modal coupling to be a possible cause. In this study we have attempted to integrate two research areas viz., proper or- thogonal decomposition and friction induced vibrations. In our attempts, we have proposed and validated the use of proper orthogonal decomposition theorem in modal reduction of systems subjected to friction-induced vibrations. The main advantage of this method lies in the fact that the modal functions are obtained from experiments and hence represent the system dynamics in an optimal sense. This however could become a major disadvantage when the spatial information is not easily available. In our preliminary investigations with strings, we encountered this problem. An imme- diate application of the resulting tool from this modal reduction is in the development of models for robust active vibration control of flexible structures where friction acts as a source of passive damping [6]. Another advantage of this method is critical in- formation pertaining to modal coupling comes as a byproduct. This information is needed in some cases to understand noise (squeak / squeal) generation [44, 80, 26] and its propagation. Finally we can use this method in transferring the vast amount of research being done on low order systems subjected to frictional excitation to high order systems which are more common in the real world. 130 7 .2 Future Work During the course of our research we found some improvements and extensions to the current work to be worthy of consideration both by ourselves and future students. Some of these are listed here. 0 We neglected the mass moment of inertia and Shear effects in the beam mod- eling. We expect the perfomance of the model to considerably improve if we take these effects into consideration. To start with one can replace the Euler- Bernoulli beam model with a simple Timoshenko beam model. However we caution that this would result in a two-fold increase in the dimension of the system [43]. e The normal and torsional vibrations of the beam could contribute to the chaos. Though we considered the former vibrations, we treated the source of these to be the slope of the Slider. There are better models in the literature that couple the normal vibrations with the transverse. e The effects of noise on POMS have been considered by some in the literature. We have not considered this in our research and they did not become an issue for us possibly because, we made use of the linear natural modes in curve fitting and in a way filtered out the noise. One can study the effects of noise on POMS and the resulting model predictions. To start with, one can use a polynomial fit rather than modal fit for these studies. 0 In the phase space reconstruction studies we have used Wolf’s algorithm [91] in estimating the Lyapunov exponents. Algorithms which give us the whole spectrum of Lyapunov exponents have been suggested by Eckmann et al. [25]. This spectrum can be used to characterize the strength of chaos. One can then study the effects of the strength of the chaotic dynamics on the predictive ability of the model. 131 o In the experiment, around 41.1 Hz we detected an audible chatter. A strobo- scopic probe did not reveal any loss of contact. We suspect the jump from the static to kinetic friction or micro-scale normal motion to be the source of this chatter. Indeed micro-scale normal motions and p5/ pg may be linked. An experimental confirmation of this suspicion would be of practical use. 0 Another interesting phenomenon that we observed in this neighborhood is the very low periodic motion of the beam. This low periodic motion at these high excitation frequencies further indicates the presence of coupling between other degrees of freedom. This needs a proper study. With this we conclude this thesis work. Bibliography [1] [2] [3] [4] [5] [6i [7] H. D. I. Abarbanel. Analysis of Chaotic Time Series. Springer-Verlag, New York, 1994. Brian Armstrong-Helouvry, Pierre Dupont, and Carlos Canudas de Wit. A survey of models, analysis tools and compensation methods for the control of machines with friction. Automatica, 30(7):]083—1138, 1994. V. Arnov, A. F. D’Souza, S. Kapakjian, and I. Shareef. Experimental investi- gation of the effect of system rigidity on wear and friction induced vibrations. Transactions of the ASME, pages 206-211, April 1983. M. T. Bengisu and A. Akay. Stability of friction induced vibrations in multi- degree-of-freedom system. In Friction-Induced Vibration, Chatter, SqueaI and Chaos - ASME, volume DE—49, pages 59—64, New York, 1992. ASME. G. Berkooz, P. Holmes, and J. L. Lumley. The proper orthogonal decomposi- tion in the analysis of turbulent flows. In Annual Review of Fluid Mechanics, volume 25. Annual Reviews Inc., New York, 1993. A. C. Bindemann and A. A. Ferri. Characteristics of passive damping in built up structures. In Friction-Induced Vibration, Chatter, Squeal and Chaos - ASME, volume DE-49, pages 173—182, New York, 1992. ASME. H. Blok. Fundamental mechanical aspects of boundary lubrication. S.A.E.J., 46(2):54—68, 1940. 132 133 [8] F. P. Bowden and L. Leben. The nature of sliding and analysis of friction. Proceedings Royal Society of London, A169z371—391, 1939. [9] F. P. Bowden and D. Tabor. Friction and Lubrication, volume 24(2). Clarendon Press, Oxford, 1950. [10] C. A. Brockley and P. L. Ko. Quasi-harmonic friction induced vibration. Trans- actions of the ASME: Journal Of Lubrication Technology, 4:550—556, October 1970. [11] D. P. Broomhead and G. P. King. Extracting qualitative dynamics from experi- mental data. Physica D, 20:217-226, 1986. [12] J. T. Burwell and E. Rabinowicz. The nature of the coefficient of friction. Journal of Applied Physics, 24(2):136-139, 1953. [13] C. Canudas de Wit, H. Olsson, K. J. Astrom, and P. Lischinsky. A new model for control of systems with friction. IEEE Transactions on Automatic Control, 40(3):419—425, 1995. [14] J. B. Conway. A Course in Functional Analysis. Springer-Verlag, 1989. [15] J. S. Courtney-Pratt and E. Eisner. The effect of a tangential force on the contact of metallic bodies. Proceedings of The Royal Society of London, A238:529—550, 1957. [16] J. P. Cusumano. Low-Dimensional, Chaotic, Nonplanar Motions of the Elastica. Ph.D. thesis, Applied Mechanics, Cornell University, 1990. [17] J. P. Cusumano and B. Y. Bai. Period-infinity periodic motions, chaos and Spatial coherence in a 10 degree of freedom impact oscillator. Chaos, Solitons and Fractals, 3(5):515—535, 1993. [18] J. P. Cusumano, M. T. Sharkady, and B. W. Kimble. Experimental measure- ments of dimensionality and spatial coherence in the dynamics of a flexible-beam impact oscillator. Philosophical Transactions of Royal Society of London, 1994. 134 [19] P. Dahl. A solid friction model. El Segundo, CA, 1968. Tech. Rep. TOR-0158 (3197-18)-1. ' [20] J. Dietrich. Micro-mechanics of Slip instabilities with rate and state dependent friction. In Eos, Trans. Am. Ceophys. Union, Fall Meeting, volume 324, 1991. [21] E. H. Dowell. The behavior of a linear damped modal system with a non- linear spring mass dry friction damper system attached. Journal of Sound and Vibration, 89(1):65—84, 1983. [22] E. H. Dowell and H. B. Schwartz. Forced response of a cantilever beam with a dry friction damper attached part i: Theory. Journal of Sound and Vibration, 91(2):255—267, 1983. [23] E. H. Dowell and H. B. Schwartz. Forced response of a cantilever beam with a dry friction damper attached part ii: Experiment. Journal of Sound and Vibration, 91(2):269—291, 1983. [24] P. E. Dupont and P. S. Kasturi. Experimental investigation of friction dynam- ics associated with normal load. In Design Engineering Technical Conferences, volume DE-Vol. 84-1, Volume 3, Part A, pages 1109—1115. ASME, 1995. [25] J. P. Eckmann and S. O. Kamphorst. Lyapunov exponents from time series. Physical Review A, 34(6):4971—4979, December 1986. [26] L. Gagliardini F. Bessac and J. L. Guyader. Coupling eigenvalues and eigen- vectors: A tool for investigating the vibroacoustic behavior of coupled vibrating systems. Journal of Sound and Vibration, 191(5):881-899, 1996. [27] B. Feeny and R. Kappagantu. On the physical interpretation of proper orthogonal modes in vibrations. submitted (Journal of Sound and Vibration), 1996. [28] Brian Feeny. Nonlinear dyanamics of stick-slip oscillations. Aerospace Structures ASME, 33:23—33, 1993. 135 [29] Brian Feeny and JinWei Liang. Phase-space reconstructions and stick-slip. Non- linear Dynamics (to appear), 1997. [30] Brian Feeny and F. C. Moon. Bifurcation sequences of a Coulomb friction oscil- [31] [32] [33] [34] [35] [36] [37] [38] lator. Nonlinear Dynamics, 4:25-37, 1993. Brian Feeny and F. C. Moon. Chaos in a forced dry friction oscillator: Experi- ments and numerical modeling. Journal of Sound and Vibration, 170(3):303-323, 1994. A. Ferri and A. Bindemann. Damping and vibration of beams with various types of frictional support conditions. Journal of Vibration and Acoustics, 114:289—296, 1992. A. Ferri and B. Heck. Analytical investigation of damping enhancement using active and passive structural joints. Journal of Guidance, Control and Dynamics, 15(5):1258-1264, 1992. Philip M. Fitzsimons and Chunlei Rui. Determining low dimensional models of distributed systems. Advances in Robust and Nonlinear Control Systems, ASME DSC, 53:9—15, 1993. A. M. Fraser. Reconstructing attractors from scalar time series: A comparison of singular system and redundancy criteria. Physica D, 34:391—404, 1989. A. M. Fraser and H. L. Swinney. Information and entropy in strange attractors. Physical Rev, 33A:1134, 1986. K. Fukunaga. Introduction to Statistical Pattern Recognition. Academic Press, New York, 1972. R. G. Gallager. Information Theory and Reliable Communication. J. Wiley & Sons, New York, 1968. 136 [39] A. Harnoy, B. F riedland, and H. Rachoor. Modeling and Simulations of elastic and friction forces in lubricated bearing for precise motion control. Wear, 172:155- 165, 1994. [40] Den Hartog. Forced vibrations with combined Coulomb and viscous friction. Transactions of the American Society of Mechanical Engineers, 53:107-115, 1930. [41] M. S. Hundal. Response of a base excited system with Coulomb and viscous friction. Journal of Sound and Vibration, 64:371—378, 1979. [42] R. A. Ibrahim. Friction induced vibrations, chatter, squeal and chaos, part i and ii. In Friction-Induced Vibration, Chatter, Squeal and Chaos - ASME, volume DE—49, pages 107—138, New York, 1992. ASME. [43] D. J. Inman. Engineering Vibration. Prentice—Hall, 1996. [44] H. Irretier. The natural and forced vibrations of a wheel disk. Journal of Sound and Vibration, 87:161-177, 1983. [45] Johannes, Green, and Brockley. The role of the rate of application of the tan- gential force in determining the static friction coefficient. Wear, 24:384—385, 1973. [46] K. Karhunen. Zur Spektral Theorie Stochastischer Prozesse. Ann. Acad. Sci. Fennicae, A 1-34, 1946. [47] M. B. Kennel, R. Brown, and H. D. I. Abarbanel. Determining embedding di- mension for phase-space reconstruction using a geometrical construction. Physics Letters A, 45(6):3403-3411, 1992. [48] V. A. Koodinov. An investigation of the vibrations of machine tools. In Proced- ings of 3rd All-Union Conference on Friction and Wear in Machines, page 161, 1960. 137 [49] D. Kosambi. Statistics in function space. Journal of Indian Mathematical Society, 7:76—88, 1943. [50] Jackqueline Krim. Friction at the atomic scale. Scientific American, pages 74-80, October 1996. [51] D. P. Lathrop and E. J. Kostelich. Characterization of an experimental strange attractor by periodic orbits. Physical Review A, 40(7):4028—4031, October 1989. [52] J. W. Liang. Characterizing the Low-Order Friction Dynamics in a Forced Os- cillator. Ph.D. thesis, Mechanical Engineering, Michigan State University, 1996. [53] J. W. Liang and B. F. Feeny. The effects of tangential contact stiffness on a forced oscillator. In Elasto-Impact and Frictionin Dynamic Systems, Proceedings of the ASME 1.996 International Congress and Exposition, November 1996. [54] M. Loeve. Probability Theory. Van Nostrand, Princeton, NJ, 1963. [55] John L. Lumley. Stochastic Tools in Turbulence. Academic Press, New York, 1970. [56] Charles R. MacCluer. Boundary Value Problems for Engineers. John Wiley, New York, 1994. [57] E. Marui and S. Kato. Forced vibration of abase-excited single degree of freedom system with Coulomb friction. Journal of Dynamic Systems, Measurement and Control, 106:280—285, 1984. [58] Leonard Meirovitch. Analytical Methods in Vibrations. Macmillan, New York, 1967. [59] Kevin D. Murphy. Using the Karhunen-Loeve decomposition to examine chaotic snap-through oscillations of a buckled plate. In Proceedings in Sixth Conference on Nonlinear Vibrations, Stability and Dynamics of Structures, Virginia, June 1996. 138 [60] M Nakai and M Yokoi. Generation mechanism of friction noises in dry friction. Japanese Journal of Tribology, 35(5):513—522, 1990. [61] M Nakai, M Yokoi, M. Inoue, and K. Kawakami. Squealing of cylindrical roller bearing. JSME International Journal, 34(1):72—81, 1991. [62] Ali Hasan N ayfeh and Dean T. Mook. Nonlinear Oscillations. John Wiley, New York, 1979. [63] J. T. Oden and J. A. C. Martins. Models and computational methods for dynamic friction phenomena. Computer Methods in Applied Mechanics and Engineering, 52:527—634, September 1985. [64] N. H. Packard, J. P. Crutchfield, J. D. Framer, and R. S. Shaw. Geometry from a time series. Phys. Rev. Lett., 45:712—716, 1980. [65] Milan Palus and Ivan Dvorak. Singular valued decomposition in attractor recon- struction: pitfall and precautions. Physica D, 55:221-234, 1992. [66] P. J. Papenhuyzen. Wrijvnigsproeven in verband methet Slippen van autobanden. De Ingenieur, 53:75-81, 1938. [67] A. Polycarpou and A. Soom. Transitions between sticking and slipping at lubri- cated line contacts. In Friction-Induced Vibration, Chatter, Squeal and Chaos - ASME, volume DE—49, pages 139-148, New York, 1992. ASME. [68] K. Popp, N. Hinrichs, and M. Oestreich. Analysis of a self-excited friction oscil- lator with external excitation. To be published, 1994. [69] K. Popp and P. Stelter. Nonlinear oscillations of structures induced by dry fric- tion. In Proceedings I U TAM Symposium on Nonlinear Dynamics in Engineering Systems, Stuttgart, 1989. [70] V. S. Pougachev. General theory of correlations of random functions. Izv. Akad. Naulc. SSSR, 17:1401-1402, 1953. 139 [71] Ernest Rabinowicz. Autocorrelation analysis of the Sliding process. Journal of Applied Physics, 27(2):131—135, February 1956. [72] Ernest Rabinowicz, B. G. Rightmire, C. E. Tedholm, and R. E. Williams. The statistical nature of friction. Transactions of the ASME, pages 981—984, October 1955. [73] P. L. Read. Phase portrait reconstruction using multivariate singular systems analysis. Physica D, 69:353—365, 1993. [74] R. S. H. Richardson and H. Nolle. Surface friction under time dependent loads. Wear, 37:87—101, 1976. [75] W. Rudin. Functional Analysis. Tata McGraw Hill, 1974. [76] A. L. Ruina. Geometry from a time series. Friction Laws and Instabilities: A Quasistatic Analysis of Some Dry Frictional Behavior, 45:712-716, 1980. [77] C. J. M. Van Ruiten. Mechanism of squeal noise generated by trams. Journal of Sound and Vibration, 120(2):245—253, 1988. [78] J. C. Schelleng. The bowed string and the player. Journal of Acoustic Society of America, 53(1):26—41, 1971. ' ' [79] J. C. Schelleng. The physics of the bowed string. Scientific American, 230:87—95, 1974. [80] E. Schneider, K. Popp, and H. Irretier. Noise generation in a railway wheels due to rail-wheel contact forces. Journal of Sound and Vibration, 120(2):227—244, 1988. [81] S. W. Shaw. On the dynamic response of a system with dry friction. Journal of Sound and Vibration, 108:305-325, 1986. [82] S. W. Shaw and C. Pierre. Normal modes of vibration for nonlinear vibratory systems. Journal of Sound and Vibration, 164(1):85-124, 1993. 140 [83] S. W. Shaw and C. Pierre. Normal modes of vibration for nonlinear continuous systems. Journal of Sound and Vibration, 169(3):319-347, 1994. [84] D. Sinclair and N. J. Manville. Frictional vibration. Journal of Applied Mechan- ics, pages 207—213, 1955. [85] S. R. Sipcic, A. Benguedouar, and A. Pecore. Karhunen-Loeve decomposition in dynamical modeling. In Proceedings in Sixth Conference on Nonlinear Vibra- tions, Stability and Dynamics of Structures, Virginia, June 1996. [86] T. Sugimoto. Stick-Slip friction. Automotive Technology, 35(4):427, 1984. [87] F. Takens. Dynamical Systems and Turbulence, Lecture Notes in Mathematics, volume 898. Springer-Verlag, New York, 1981. [88] D. M. Tolstoi. Significance of the normal degree of freedom and natural normal vibrations in contact friction. Wear, 10:199-213, 1967. [89] W. W. Tworzydlo and E. B. Becker. Influence of forced vibrations on the static coefficient of friction-numerical analysis. Wear, 43:1-22, 1991. [90] W. W. Tworzydlo, E. B. Becker, and J. T. Oden. Numerical modeling of fric- tion induced vibrations and dynamic instabilities. In Friction-Induced Vibration, Chatter, Squeal and Chaos - ASME, volume DE—49, pages 13—32, New York, 1992. ASME. [91] A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano. Determining Lyapunov exponents from a time series. Physica D, 16:285-317, 1985. [92] K. Yagasaki. Chaotic dynamics of a quasi-periodically forced beam. Journal of Applied Mechanics, 59:161-167, 1992. [93] K. Yagasaki. Bifurcations and chaos in quasi-periodically forced beam: Theory, Simulation and experiment. Journal of Sound and Vibration, 183(1):1—31, 1995.