Ali

Paul L Olivier

Alessandra Alexandre

Abir, Vahan, Gang, Belgasem, Ding, Ramzi, Drezet, Bérenger, Youssef, Ludovic, Gemala, Flaviano Abd El Hamid

Emile Romeo

Aymen, Mohamed Damien Leandro

Karim

Vincent P Dimitri

mother Norma Fawaz

brother Youssef Aline Many

Jinane, Micheline, Mostafa, Ibrahim Anna Maria

Mohamad H Adham

Keywords: Periodic lattices, nonlinear oscillators, collective dynamics, weak coupling, modal interactions, basins of attraction, localization phenomenon, solitons. Methods for solving nonlinear problems

I am thankful for the friendships I developed while at the applied mechanics department

Acknowledgement

Making the commitment to undertake this PhD as a career decision has been a lifechanging and challenging experience that would not have the possibility to end up with a success without the guidance and support that I received from many people. I believe in William Raduchel quote that says: Don't pick a job. Pick a Boss. Your first boss is the biggest factor in your career success. A boss who doesn't trust you won't give you opportunities to grow. I couldn't have chosen more wisely, they inspired me at several scientific levels and ignited in me the love for research and curiosity for further investigations.

For this reason, I would like to express my sincere appreciation to my director, Prof. Noureddine Bouhaddi for his immeasurable support, motivation during hard times and invaluable guidance and instructions. I am grateful for the opportunities he gave me to present and participate to several international conferences, where I gained a large scientific exposure and knowledge. In addition, I will never forget and deeply appreciate his carefulness and patience.

Abstract

Localization phenomenon has attracted significant interest recently due to its important role in the qualification as well as quantification of periodic systems dynamics. The localization can be generated either by the disorder (imperfection) in a lattice or by the interaction between nonlinearities in a discrete system. One of the most frequent localization phenomena is the nonlinear energy localization, called intrinsic localized modes or solitons.

Although the dynamics of periodic nonlinear lattices was thoroughly investigated in the frequency and time-space domains, there is a real need to perform profound analysis of the collective dynamics of such systems in order to identify practical relations with the nonlinear energy localization phenomenon in terms of modal interactions and bifurcation topologies. The principal goal of this thesis consists in exploring the localization phenomenon for modeling the collective dynamics of periodic arrays of weakly coupled nonlinear resonators.

An analytico-numerical model has been developed in order to study the collective dynamics of a periodic coupled Duffing-Van Der Pol oscillators array under simultaneous primary and parametric excitations, where the bifurcation topologies, the modal interactions and the basins of attraction have been analyzed. Two applications were considered: (i) an array of coupled pendulums under harmonic external excitation, where it has been shown that by increasing the number of coupled pendulums, the number of multi-modal solutions increases and the distribution of the basins of attraction of the resonant solution branches becomes larger, (ii) an array of electrostatically coupled nanobeams under parametric excitation, while including the main sources of nonlinearities up to the fifth-order. The model was extended to investigate the collective dynamics of periodic nonlinear two-dimensional arrays, and two applications have been considered: (i) a 2D array of coupled pendulums under harmonic base excitation, where it has been shown that considering a high ratio between the coupling parameters in the two directions leads to additional features distinct from those obtained in 1D arrays, mainly a larger bandwidth and important vibrational amplitudes, (ii) a 2D periodic granular crystals array subjected to compressive loadings and a harmonic horizontal base excitation, where each particle is coupled with fixed-fixed uniform beam while the nonlinear dynamics is governed by the Hertzian contact. A second investigation of this thesis consist in transforming the developed discrete model into a continuous one, where the dynamics is generated

Résumé

La localisation joue un rôle important pour de nombreuses applications en physique car elle permet de quantifier et de qualifier la dynamique d'un système périodique. Elle peut être générée, soit par le désordre (imperfections) dans un réseau ou par l'interaction entre les non-linéarités dans un système discret. L'un des phénomènes de localisation le plus courant est celui de la localisation non-linéaire de l'énergie, connu sous le nom de modes intrinsèques localisés (ILMs) ou Solitons.

Bien que la dynamique des réseaux périodiques non-linéaires ait été investiguée dans les domaines temporel et fréquentiel, il existe un réel besoin d'identifier des relations pratiques avec le phénomène de la localisation d'énergie en termes d'interactions modales et topologies de bifurcation. L'objectif principal de cette thèse consiste à exploiter le phénomène de la localisation pour modéliser la dynamique collective d'un réseau périodique de résonateurs non-linéaires faiblement couplés.

Un modèle analytico-numérique a été développé pour étudier la dynamique collective d'un réseau périodique d'oscillateurs couplés de type Duffing-Van Der Pol sous excitations simultanées primaire et paramétrique, oú les interactions modales, les topologies de bifurcations et les bassins d'attraction ont été analysés. Deux applications ont été investiguées: (i) un réseau de pendules couplés sous excitation extérieure harmonique, oú il a été démontré qu'en augmentant le nombre de pendules, le nombre de solutions multimodales et la distribution des bassins d'attraction des branches résonantes augmentent, (ii) un réseau de nano-poutres couplées électrostatiquement sous excitation paramétrique en incluant les principales sources de non-linéarité jusqu'à l'ordre cinq. Ce modèle a été étendu pour investiguer la dynamique collective d'un réseau périodique non-linéaire bidimensionnel, oú deux applications ont été étudiées: (i) un réseau 2D de pendules couplés sous excitation extérieure harmonique, oú il a été démontré que le choix de différents paramètres de couplage dans les deux directions, enrichit la dynamique collective comparée à celle obtenue dans les réseaux 1D, et conduit principalement à des amplitudes de vibration plus importantes et des bandes passantes plus larges, (ii) un réseau 2D de billes métalliques en compression, guidées élastiquement par des poutres et sous excitation à la base, où la non-linéarité est gouvernée par le contact de Hertz. Une deuxième investigation de cette thèse consiste à transformer le modèle discret développé en modèle continu dont la dynamique est régie par l'équation de Schrödinger. Ce modèle permet d'identifier les solitons associés à la dynamique collective du réseau de Résumé structures périodiques et d'étudier leurs stabilités.

Mots clés: réseaux périodiques, oscillateurs non-linéaires, dynamique collective, couplage faible, interactions modales, bassins d'attraction, localisation d'énergie, solitons.

Introduction

Motivations

Arrangements of periodically repeated united cells, known as periodic structures are ubiquitous. They exist naturally as can be man-made over a large range of length scales, with a wide range of engineering applications, including photonic crystals, elastic plates and shells, honeycombs, civil constructions, antennas, mechanical constructions, MEMS/NEMS and so on. The motivations behind studying similar structures, lie in their ability to exploit periodicity in order to attenuate, isolate or localize vibrations. Their simple fabrication was the main advantage, which motivated researchers in several domains to emerge innovative applications and develop new devices, to explore in deep wave propagation in periodic structures.

In fact, periodic arrangement properties present a particular capacity to attenuate the propagation of elastic waves [Romeo 03], reduce the transmitted vibrations and noise [Li 17] or filter them [ Barbarosie 04] as the dichroic filters which is realized on the concept of photonic crystals [Winnewisser 98]. In addition, interesting electromagnetic properties arise from periodic structures as when white light passes through a diffraction grating (optical component with a periodic structure) or photonic crystal fibers and integrated optical devices based on photonic crystals and periodic arrays of holes [Broeng 99] and frequency-selective surfaces which are periodic surfaces that can be found on a plane, designed to reflect, transmit or absorb electromagnetic fields based on frequency [Mittra 88]. Moreover, it has been shown that mode localization can occur in weakly coupled quasi-periodic linear systems [Hodges 82], resulting to the Anderson localization [ Anderson 58].

Despite the interesting wave properties encountered in periodic linear system as filtering, bandgaps, propagation and attenuation zones, isolate or localize vibrations, periodic nonlinear structures may exhibit variety of rich dynamics properties. For instance, bifurcations, solitary waves and solitons [Cuevas 09], energy dependent nonlinear propagation/attenuation zones [Romeo 06], nonlinear modes localization [Vakakis 97], multitude of nonlinear resonances [Nayfeh 08], wave interaction, spatial and temporal chaos [Waller 84] and energy transfer/harvesting/trapping [ Lydon 15].

In fact, space localization represents an interesting phenomenon in engineering science, which attracted significant interest recently due to its important role in the qualification as well as quantification of periodic systems dynamics. It is based on the fact that the energy does not propagate arbitrary and it is localized in space, while Introduction the wave-function amplitude of the oscillating modal shape decays exponentially in space. This phenomenon has inspired innovative studies in physics and motivated researchers over many years to explore in depth its effects and consequences. The localization can be generated either by the disorder (imperfection) in a lattice or by the interaction between nonlinearities in a discrete system. One of the most frequent localization phenomena is the nonlinear energy localization, called intrinsic localized modes or solitons.

Although the dynamics of periodic nonlinear lattices was thoroughly investigated in the frequency and time-space domains, there is a real need to perform profound analysis of the collective dynamics of such systems in order to identify practical relations with the nonlinear energy localization phenomenon in terms of modal interactions and bifurcation topologies. The principal goal of this thesis consist in exploring the localization phenomenon for modeling the collective dynamics of periodic arrays of weakly coupled nonlinear resonators.

Overview

Chapter 1 presents a state of art of the dynamics of periodic structures, their existence and the characteristics of their assemblies. Then, methods for dealing with wave propagation in multidimensional periodic linear systems are presented, with the highlight of the interesting wave properties that can produce. Variety of rich dynamical nonlinear phenomena motivated the study of wave propagation in periodic nonlinear lattices. In particular, the review of several physical applications that can be modeled as periodic nonlinear lattices is given; in particular, pendulums arrays, MEMS/NEMS arrays and granular crystal lattices. The history of localization phenomenon, the presentation of the techniques employed to investigate localization in discrete and continuous quasi-periodic structure and spatial localization in weakly coupled perfectly periodic and quasi-periodic structures are emphasized.

Chapter 2 is devoted to covering most of the numerical, analytical or combination of both methods employed to study periodic nonlinear lattices. The direct time integration, the asymptotic numerical method (ANM), ANM with the harmonic balanced method (HBM), the multiple scales method with and without standing waves decomposition are presented. An application consisting of periodic two weakly coupled mass-spring system is considered and a comparison between the different used methods is presented. The goal is to validate the capability of the multiple scales method coupled with standing wave decomposition to handle weakly coupled nonlinear periodic oscillators.

In Chapter 3, we consider the general case of an array of weakly coupled Duffing-Van Der Pol oscillators with fifth order nonlinearity under simultaneous primary and external excitations. An important choice of weakly linear coupling parameter between substructures is made enabling the creation of linear closed modes. The coupled nonlinear differential system is solved using an analytico-computational model, based on the method of multiple scales coupled with standing wave modal decomposition, transforming the nonlinear differential system into a set of coupled complex algebraic equations which are numerically solved using the ANM enabling the construction of resonance curves for a large number of dofs. The complexity and the multivaludness of the responses are illustrated by a detailed study of its basins of attraction. The ability of the combination of both primary and parametric excitations to adjust the stability of the system is shown.

Chapter 4 is dedicated to the study of with the collective dynamics of a periodic array of linearly coupled pendulums under primary resonance. A model reduction method is proposed to calculate the dominant dynamics without significant loss of accuracy compared to the full model. The modal interactions and their effects on the nonlinear dynamics and their bifurcation topology are studied. The model is extended to investigate the collective dynamics of 2D coupled pendulums arrays, under harmonic base excitation, where additional feature not obtained in the 1D array are achieved.

Chapter 5 studies the collective dynamics of nonlinear oscillator's arrays with weakly nonlinear coupling, for two different physical applications. An electrostatically coupled nanobeams array under parametric excitation is modeled and investigated, while including the main sources of nonlinearities up to the fifth order. The influence of the fifth order nonlinearity on the complexity of the frequency responses is highlighted. In addition, we investigated, the collective dynamics of a 2D periodic lattice of granular particles coupled with identical beams and subjected to a precompression loads where the nonlinearity is governed by the Hertzian contact law, under horizontal base excitation.

Chapter 6 deals with a different type of behavior in nonlinear periodic arrays, called the intrinsic localized modes or solitons. The multiple scales method is used to transform the differential nonlinear system into an amplitude equation in the form of a Nonlinear Schrödinger (NLS) equation. Exact analytical solitary solutions of the conservative NLS equation are derived, while the dissipative system is solved numerically using the analog Newton's method. Several numerical simulations are performed demonstrating that the stability of localized solutions can benefit from the combination of both parametric and external excitations.

Overview

The current research interests focus on the collective dynamics to explore the localization phenomena in periodic nonlinear lattices. The upcoming section starts introducing periodic systems along with some definitions and examples on both naturally existing and man-made periodic structures. Next, a literature review on the dynamics of multidimensional linear periodic structures is introduced along with the methods devoted to study wave propagation in periodic media such as Floquet, Block, and finite element method and so on. Although linear periodic structures exhibit a number of interesting wave phenomena as bandgaps, propagation and attenuation zones, vibration isolation and localization, variety of rich dynamic properties such as wave localization, solitary waves, energy transfer and frequency-conversion emerged the trend to study the wave propagation in nonlinear periodic systems. Then, a review on the state of art of the different complex dynamical phenomena encountered in periodic nonlinear structures is reported. Three famous physical examples: pendulums, MEMS/NEMS and granular crystals that can be modeled as periodic nonlinear lattices are presented, along with a brief literature review, background study and properties. A historical background, highlighting the existing of the localization phenomena is presented followed by a detailed study on the space localization and its properties in periodic linear and nonlinear arrays. The research motivations, the objectives and the challenges of studying the collective dynamics of periodic nonlinear lattices are addressed and the methodology employed to achieve the present goal is announced. 1.2 Periodic structures

Periodic structures

Periodic structures consist of an arrangement of coupled identical substructures, geometrically repeated and defined by a unit cell at periodic intervals. They can exist naturally (figure 1.1) as they can be man-made (figure 1.2) over a wide range of length scales. For instance, crystal lattices are three-dimensional arrangement of unit cells, made of atoms, ions and molecules. Interactions between these last are characterized by potential energy that result in a series of vibrational waves propagations (Phonons). On the other hand, the origin of the shinning colors of some insects lies in a distinct nanostructure built into their biopolymeric wing scales and exoskeletons which have a high degree of periodicity. Particularly, the propagation of light through the polymeric structures (diffraction and specular reflection) those with nanostructures periodicity produce angle-dependent iridescent colors. Furthermore, honeycombs consist of hexagonal cells built by bees, where their periodicity influence the way waves propagate.

In regards to man-made periodic structures, we can mention many architectural designs as long span bridges and high-rise buildings, aircrafts wing ribs in aerospace, micro-and nano-electromechanical systems, LED arrays in electronics and so on. In addition, they can support many applications as photonic crystals cavities and bandgaps, electromagnetic/optics devices, frequency selective surface or metamaterials. From an engineering viewpoint, the principal advantage of such structures is due to the simplicity of their geometry and manufacturing, the significant impact and high temperature tolerance and their high strength-to-weight ratios. In addition, periodic structures provides interesting wave propagation properties as mechanical filters or waveguides. Moreover, the motivations behind studying similar structures, lie in their ability to exploit periodicity in order to attenuate, isolate or localize vibrations. Therefore, the comprehension of the dynamics of periodic structures is necessary for analyzing and synthetizing the role of periodicity in generating localization phenomena.

Linear systems

Several researches on the dynamics of one-dimensional periodic structures were based on linear model structure. Under the hypothesis of linearity, analytical and numerical methods have been developed in order to analyze their dynamics, where the study of wave propagation in these structures was thoroughly addressed in the literature.

Methods for dealing with linear periodic systems

Brillouin [Brillouin 53] traces a detailed historical background and motivation behind treating wave propagation in the field of solid-state physics, by considering a classical mass-spring system. Mead and his research group [Mead 73,Mead 96] developed methods and techniques based on wave modes to analyze and predict free vibration propagation and forced vibration induced by stationary harmonic loads in continuous periodic engineering structures. In the field of structural mechanics, Langley [Langley 96] studied the response of a two-dimensional linear periodic structure subject to point harmonic forcing using a spring mass model, showing that the location and the extension of the stop bands depend on the considered direction of wave propagation on the structure plane. Later, Langley et al. [Langley 97] presented theoretical and experimental studies to the response of a beam grillage constructed from strips of aluminum with bolted joints (figure 1.3 (a)), using the finite element analysis.

About 50 years ago, the dynamics of wave propagation in periodic linear systems has been investigated by means of transfer matrix method. This method allow reducing computational effort and time of the considered problem to the number of degrees of freedom coupling basic periodic element. Lin et al. [Lin 69] developed an analytical technique applying a transfer matrix for the determination of the frequency response of stiffened plate vibration. A review on the theory of wave propagation in a survey physical applications of one-dimensional locally periodic media using the transfer matrix approach has been reported by Griffiths [Griffiths 01]; proving that as the number of coupled resonators increases, the band structure characteristic of waves in infinite periodic media emerges. Vibrational responses of 1D filter and 2D wave guide modeled as mass-spring structures subjected to periodic loading has been analyzed, while studying the effects of boundaries, viscous damping and imperfections [Jensen 03]. The well know formulation of Floquet [Floquet 83] is a methodology to solve differential equations with periodic coefficients, assuming that any solution of the system can be expressed as a linear combination of functions.

Although Bloch theorem [Bloch 29] was firstly introduced to represent the model of homogenous states of Schrödinger equation with periodic potential, it can be considered as a multidimensional application of the Floquet theorem. The Floquet-Bloch theorem is extensively used for calculating the dispersion properties, wave modes and group velocities of periodic structures. Collet et al. [Collet 09,Collet 11] developed new multifunctional structures integrating electro-mechanical systems in order to optimize their vibro-acoustic behavior over a large frequency band and introduced a Floquet-Bloch decomposition for the computation of dispersion of two-dimensional periodic (figure 1.3 (b)), damped mechanical systems. Recently, an efficient method originally proposed by Zhou et al. [Zhou 15a,Zhou 15b] was illustrated by one and two-dimensional periodic structures; this method is based on a Component Mode Synthesis (CMS) approach to analyze the local behavior of the unit cell using a reduced modal basis in addition to a Wave Finite Element Method (WFEM) that fully exploit the periodic propriety of the structure and extracts directly the propagation parameters. A high number of researches on wave propagation in periodic structures is based on studying linear structural designs. Although these latter produce interesting wave properties as filtering the propagation of waves, bandgaps, propagation and attenuation zones and isolate or localize vibrations, complex dynamical phenomena may be detected in nonlinear periodic structures, which can be exploited to design periodic structure-based smart systems with high performance. In addition, and since most of the structures are nonlinear due to their nature, the interactions with their neighbors, the restoring forces and so on, it's essential to analyze wave propagation in periodically nonlinear systems.

Nonlinearity in periodic systems

Variety of rich dynamic properties such as multi-stability, bifurcations, existence of highly stable localized solutions, solitary waves and solitons, energy transfer, frequency-conversion and variation in wave speeds and propagation direction related to wave amplitude and nonlinearity etc. emerged the trend to devote significant studies to explore these complex phenomena for wave propagation in nonlinear periodic media. For instance, the existence of waves and amplitude dependent frequencies bounding nonlinear propagation and attenuation zones (PZs and AZs) has been found by studying an infinite chain of mono-coupled oscillators with geometric nonlinearities [ King 95]. The existence of solitary wave propagation in strongly nonlinear periodic structures such as uncompressed granular media was first revealed by Nesterenko [Nesterenko 83]. Soliton-defect interactions using simple two-mode models characterized by a mobile soliton and a localized (trapped) mode were first suggested by Goodman et al. [Goodman 04]. Localization and solitary wave phenomenon were observed in a perfectly cyclic, weakly nonlinear periodic system [Vakakis 93]. Nayfeh et al. [Nayfeh 08] explored the multitudes of nonlinear resonances such as subharmonic, superharmonic and combination resonances by applying several perturbation techniques for predicting wave propagation in weakly nonlinear continuous systems. In addition, the synchronization of chaotic circuits has been studied in random network such as secure communication [Zhou 05]. Moreover, it was found that the presence of defects and interfaces in granular media is interesting to control stress propagation and produces energy trapping [ Daraio 06].

Chapter 1 is dedicated to the theoretical framework of methods for solving periodic nonlinear problems, with a short literature review on the methods employed to solve similar systems. In addition, many physical systems in several scientific domains such as materials, acoustics, optics, mechanics, MEMS/NEMS and vibrations can be modeled as periodic nonlinear chains. In the following section, we will be interested in three principal physical application that can be modeled as arrays of coupled nonlinear oscillators. In mechanics, the array of coupled pendulums (figure 1.6 (b), (c)) represents a famous example of periodic nonlinear structures which can be described by the Sine-Gordon model . The latter is suitable for a large variety of physical systems, which explains the vast research area developed around this field. Several theorems on the existence of oscillatory, rotary, and mixed periodic motions of coupled pendulums were proved [ Marlin 68]. In addition, Cai et al. [Cai 93] identified two distinct regimes for a damped Sine-Gordon chain driven by the spatio-temporal periodic potential. The nonlinear dynamics and chaos in a simple pendulum array have been reported experimentally, reproducing the fact that if there was no disorder a state of fully developed spatiotemporal chaos appears, while an optimal amount of disorder returns the dynamics regular [ Shew 99]. In addition, the synchronization of parametrically excited coupled chaotic pendulums has been experimentally observed [Zhang 01]. Vyas et al. [Vyas 01] showed that the bandwidth of absorber effectiveness can be increased substantially by using an array of pendulums with slightly different natural frequencies. For a chain of coupled pendula driven periodically at one end (figure 1.6 (b)), Khomeriki et al. [Khomeriki 08] demonstrated experimentally the existence of a novel regime which produces an output frequency at an odd fraction of the driving frequency. Chen et al. [Chen 07] showed using the Floquet theory that the instability onset in a pendulum array always responds to the forcing frequency subharmonically. Recently, Liang et al. [Liang 15] reported an experiment using a coupled-pendulums chain to demonstrate various aspects of wave motion. Further-more, Jallouli et al. [Jallouli 15] investigated the nonlinear dynamics of a 2D array of coupled pendulums under parametric excitation. In other studies, Intrinsic Localized Modes (ILMs), existing as exact solutions in nonlinear lattices have been widely studied. Particularly, for an array of coupled pendulums [Ikeda 15] which is an ideal playground for showcasing many interesting features of discrete breathers (figure 1.7 (a)), [Cuevas 09] which can be stabilized in the presence of driving and damping [English 14] or propagate energy in a forbidden band gap [Geniet 02]. The nonlinear coupling and decoupling between ILM and third allowed plane-wave mode, as well as both modes interaction with the driver, give rise to pulsating breather instability [Thakur 08]. In addition, it was shown that the mass impurity in a parametrically driven, damped nonlinear coupled pendula has striking influence on the high-frequency modes [Hai- Qing 06]. Recently, the influence of adding external harmonic excitation on the intrinsic localized modes of coupled pendulums chains parametrically excited has been investigated [ Jallouli 17].

MEMS/NEMS arrays

Over the last decade, there has been an increasing attention in the field of MEMS & NEMS dynamics as an important process technology, induced by a set of functional needs which have the potential to affect our lives as well as fundamental questions [ Bao 96]. These small integrated devices are used to create tiny systems that combine electrical and mechanical components, having high frequencies contrary to their microscales. There are plenty of new challenges in this area of researches that results from the request to manufacture new MEMS devises that improve the reliability and performances of the existing ones.

In parallel, there was the effort of miniaturization for designing Nanoelectromechanical systems (NEMS), with a reduced size that is less than one micrometer. This down scaling from MEMS to NEMS has permitted the fabrication of very small devices having faster speed, higher complexity, lower power and cost in the domain beyond electronics, that allow for better performing sensors. While operating at resonance, NEMS display high fundamental frequencies, small masses, tolerable linear stiffness and they are subject to very low damping. These characteristics enable NEMS to be a perfect candidate for several applications like inertial sensors, gas sensing, filters and switches [Craighead 00,Ekinci 05]. The capability to produce arrays of high number of coupled beams in very large scale integration fashion opens new possibilities to study the nonlinear dynamics of coupled MEMS and NEMS. These devices are compatible and co-integrable with standard silicon electronics processes and one can take advantage of their resulting collective dynamics to enhance the performances of MEMS/NEMS sensors and actuators. The collective responses of coupled arrays could be beneficial for signal improvement and noise decrease , as well as sophisticated mechanical signal processing applications. To possess rich dynamic behavior unreachable with uncoupled MEMS resonators, Buks and Roukes [Buks 02] studied the mechanical characteristics of an array of parametrically excited micromechanical doubly clamped gold beams (figure 1.9 (a)). Qualitative comparison with the experiment of Buks and Roukes was done by Lifshitz and Cross [Lifshitz 03], by modeling a parametrically excited MEMS array using a set of coupled discrete nonlinear equations of motion. They used the secular perturbation theory and standing waves decomposition to transform the coupled differential equations of motion into a set of coupled nonlinear complex algebraic equations for the normal modes amplitudes. Although, this procedure helped them to obtain exact solutions for few coupled oscillators, they gave a qualitative understanding of the dynamics of large arrays. Later, Bromberg, Cross and Lifshitz [Bromberg 06] studied the same equations of motion to obtain analytical results for large arrays by approaching the system from the continuous limit of infinitely large degrees of freedom, and obtaining a description of the slow spatiotemporal dynamics of the array of oscillators in terms of an amplitude equation. Moreover, Kenig, Lifshitz and Cross [Kenig 09a] extended the investigations of the amplitude equation to study the problem of pattern selection with the application to MEMS/NEMS systems. They described the transitions between standing-wave patterns of different wave numbers as the drive amplitude is varied either quasistatically, abruptly, or as a linear ramp in time. A review of the nonlinear dynamics and complexity of such nanomechanical and micromechanical resonators was conducted by Lifshitz and Cross [ Lifshitz 10]. The nonlinear dynamics of electrostatically coupled microresonators (figure 1.9 (b)) in several frequency regions have been performed by Gutschmidt et al. [Gutschmidt 12] using a continuum model. We should note that, such arrays may exhibit interesting dynamical phenomena such as the intrinsic localized modes [Sato 03c, Sato 03a], which make them adequate for analyzing the dynamics of large degrees of freedom systems.

Granular crystal lattices

Hertzian chains of spherical particles represent a famous example of nonlinear discrete structures, where the nonlinearity emerges from two characteristics; the first one is the geometry between particles is such that the transmission between two particles in contact is nonlinear and the second one from the applied compressive load. The dynamic behavior of one-dimensional chain of granular crystals emerged by Nesterenko [Nesterenko 83] showing analytical and numerical studies. Later, the concept of "sonic vacuum" and the formation and propagation of highly nonlinear solitary waves in similar arrays was experimentally shown by Lazaridi et al.

[ Lazaridi 85]. The motivation behind studying granular crystals lies from their ability to control the nonlinearity and tune their dynamic responses from linear to weakly and strongly nonlinear by changing the amount of the static precompression loads [ Nesterenko 01]. In addition, they have a great controllability of the assembly as well as the simplicity of their construction and their applicability to engineering devices [Boechler 11,Spadoni 10]. Granular crystals showed the ability to support tunable acoustic band gaps in their linear and weakly nonlinear regime, where a novel perturbation analysis verified by numerical simulations was applied to derive amplitude dependent dispersion relations [Boechler 10a]. In addition, it was observed that the wave speed changes with increase in amplitude at certain frequency leading to small changes in the band edges [Boechler 10b]. Several studies have reported, both experimentally and numerically, energy localization phenomena in Hertzian chains with mass defects considering a one-dimensional diatomic granular crystals [Boechler 10c]. Experimental observation of the modulation instability and discrete breathers in the weakly nonlinear dynamical regime in such arrays was reported [ Theocharis 10]. Heterogeneous, 1D strongly compressed spheres chain has been used to generate solitary waves [ Yang 11]. In the strongly nonlinear regime it has been demonstrated that a sequence of pulses is induced at the interface of two sonic vacua by the deceleration of interfacial particles [ Nesterenko 05]. In addition, Daraio et al. [Daraio 06] demonstrated experimentally and numerically the efficiency of soliton-like and shock-like pulse trapping and disintegration in a composite granular protector and proved that its efficiency depends on the particle's arrangements. Moreover, Job et al. [Job 07] described through experimental, numerical and theoretical studies the formation of solitary wave train by large striker on an alignment of spheres. On the other hand, Theocharis et al.

[ Theocharis 13] reported the dynamic behavior of nonlinear periodic phononic structures, along with how such structures can be utilized to affect the propagation of mechanical waves.

Several studies were devoted to investigate the dynamics of two-dimensional granular lattices, motivated by the richness of the complex dynamical nonlinear phenomena found in one-dimensional ones. Recently, Leonard et al. [Leonard 13]

Localization phenomena

Space localization represents an interesting phenomenon in engineering science, which can occur in periodic structures when the wave-function amplitude of the oscillating modal shape decays exponentially in space. It is based on the fact that the energy does not propagate and it is localized in space. Consequently, researchers were motivated to explore more deeply the effects and consequences of the localization, where they demonstrated analytically and numerically that for linear periodic systems mode localization exists when: (1) the coupling between substructures is weak enough and (2) the periodicity in the structure is perturbed by structural mistuning. However, for perfectly periodic nonlinear systems, strongly localized nonlinear normal modes may occur for sufficiently weak substructures coupling.

History of localization studies

In the previous sections, the dynamics of perfectly periodic structures encountered in many field of engineering and physics have been examined. Manufacturing and assembly defects may produce disorder that can occur in material properties of the system. Spatial localization of normal modes and attenuation of waves in all frequency bands could be the consequence of a disorder in periodic structures, in which energy is confined to the region near disorder where the dynamics behavior of the structure changes.

Starting from the fifties, two notable papers the first one by Dyson [Dyson 53] explored the effect of disorder of a one-dimensional oscillators, coupled by harmonic forces where the inertia of each oscillator and the strength of each coupling were considered as random variables. A vibrating mass-spring system has been considered where a method for calculating the distribution-function of the frequencies of normal modes has been presented. In addition to the frequency spectra, was calculated analytically when the considered distribution law of the oscillator parameter takes an exponential form. The second paper, Schmidt [Schmidt 57] developed a mathematical method which gives fairly generally the density of eigenstates for linear disordered chain system. Although, their results helped explaining in a mathematical way the wave transmission in a disordered periodic systems, they ignored the effect of disorder on the eigenvectors and wave propagation.

(1)

(2) In quantum mechanics, the first scientist to describe eigenstate localization by the absence of diffusion of waves in a disordered potential is the American physicist Philip Warren Anderson [Anderson 58]. He showed that an electron in a semiconductor has the possibility of not being transported, provided that the number of disorder is sufficiently large in a lattice. This phenomenon is named Anderson localization in some cases token of appreciation for his original contribution. Anderson suggestions motivated researchers over many years to well understand localization. This phenomenon has been thoroughly investigated in mistuned periodic structures, inspiring innovative studies in physics. Early studies treated disordered chains of atoms, it has been proved that any infinite disordered chain can be made to have an exponentially localized solution (figure 1.12) for any frequency by modifying the mass of one atom to a suitable real value as a function of frequency [Matsuda 70].

Localization phenomena

In addition, a general proof for the exponential growth of the particular solutions of the stationary equation of motion in the system was given later by Ishii [Ishii 73].

Localization in linear quasi-periodic structures

In the field of structural dynamics and vibrations, Hodges [Hodges 82] was the first to study localization phenomena in disordered periodic structures. He used simple examples such as coupled pendula (figure 1.13) and strings chains to prove that disorder in periodicity could have amazing impact on the structural dynamics and waves attenuation in all frequency bands independently of any dissipation in the system. Later, Hodges et al. [Hodges 83] measured the transmission of energy and the mode shapes in a mass-spring system after moving the masses in a controlled way to provide small degree of irregularities. They satisfied the theoretical predictions by proving that the individual modes making up each passband are localized in specific areas of the disordered structure.

Techniques to investigate mode localization in discrete systems

The possibility of localization in large space structures including the effect of structural imperfection and disorder was investigated theoretically and numerically, proving that mode localization is most likely to occur in structures consisting of a large number of weakly coupled substructures [ Bendiksen 87]. Pierre et al. [Pierre 87b] considered a multispan beams on irregularly spaced supports to study localization using a perturbation technique. In a disordered system consisting of identical nearly coupled component systems, strongly localized modes of vibration may occur if the coupling between the substructures is of the order of or smaller than the spread in natural frequencies of the component system [Pierre 87a,Pierre 88]. Later, Li et al. [Li 92] gave a brief review underlying the effects of near-periodicity on both linear waves and dynamic behavior in the interest of engineering structures. For one-dimensional chain of coupled disordered, Ottarasson et al. [Óttarsson 96] studied mode localization in mistuned bladed disk using random transfer matrices where they used a perturbation technique to obtain an approximation of the localization factor and the spatially averaged rate of exponential attenuation.

Techniques to investigate mode localization in continuous systems

Several researchers thoroughly analyzed the dynamics of continuous disordered systems such as chains of beams, panels, blades, patches or strings using different techniques. Bendiksen et al [Bendiksen 89] investigated mode localization phenomenon in a generic class of large space reflectors, using a Rayleigh-Ritz formulation for few cantilevered beam bending modes and a finite-element formulation using Bernoulli-Euler beam elements. A new perturbation scheme for damped disordered periodic structures, which permits to calculate the localization factor which means the average exponential decay rate of wave transmission with respect to the distance of transmission has been developed [ Cai 91]. In disordered cyclic structures, a localization length scale is proposed as a measure of the spatial extent of a mode using numerical experiments indicating that the particular set of random disorder can have a significant effect on the degree of localization [ Cornwell 92]. Axially vibrating beam continuously restrained by imperfect elastic springs have been analyzed by [Luongo 92] using a deterministic approach. Pierre et al. [Pierre 89] developed a perturbation method which was applied for mistuned assemblies of coupled Mdof component systems, showing that higher modes are more susceptible to localization, which is unavoidable if the mode number is large enough.

More recently, Localization phenomena in one-dimensional imperfect continuous structures have been analyzed in both dynamics and buckling, where Luongo developed a general perturbation method which generalizes the classical WKB approximation (in Honor of Léon Brillouin, Hendrik Anthony Kramers and Gregor Wentzel) which consists in finding approximate solutions to linear differential equations with spatial varying coefficients [ Luongo 01]. Disorder in periodically shunted piezoelectric patches extends the stop bands into two adjacent propagation zones and produces the localization of the vibration energy near the excitation zone [ Thorp 01]. The wave localization in disordered periodic 2-2 piezoelectric composite structures is strengthened due to the piezoelectricity, where a slight disorder can lead to prominent localization phenomenon [Li 05]. Huang [Huang 06] Used the Galerkin method to discretize the localization equation of a continuous system constituting a number of several mistuned blades, revealing that the number of grouped blades and the distribution of multi-disorder in a rotating blade system my markedly affect the localization phenomenon. A review on the literature of mistuned blade disks vibration has been given by Castanier et al. [Castanier 06] enabling better understanding of their forced responses and highlighting important developments and emerging directions in this research area. The Lyapunov exponent has been used in order to calculate the localization factors in a disordered periodically stiffened double-leaf panels, where the numerical results show the significant effect on the pass bands and the localization factor [Yu 10] and the localization characterizing the average exponential rate of decay of the wave amplitude in a disordered piezoelectric axial-bending coupled beams [ Zhu 13].

Recently, an homogenous approach has been applied to highlight the frequency split and localization in disordered rotationally periodic structure [ Bisegna 11]. The frequencies of the wave localization in periodic composite materials with local defects have been analyzed based on the transfer-matrix and the plane-wave expansion methods [Andrianov 14]. Ding et al. [Ding 16] modeled a disordered tunnel with a pipe-beam with periodic joints on elastic foundation, localization factors were calculated showing that the flexural wave is always localized near the moving load when the velocity is less than some critical velocity of a uniform tunnel on elastic foundations. In addition, a numerical approach using the modal data has been suggested in order to quantify vibration mode localization in complex systems using the modal criterion across all the modes sue to changes in some systems parameters [ Chandrashaker 16].

Nonlinear localization

Although, mode localization in linear periodic structures exists only when the coupling is sufficiently weak and the structure is perturbed by small mistuning, It has been shown that for sufficiently weak substructure coupling in perfectly nonlinear periodic structure, strongly localized Nonlinear Normal Modes (NNMs) may occur.

Mode localization in perfectly periodic structures

In the theory of nonlinear waves, spatial nonlinear localization is one of the most important properties encountered in NNMs, where it provides a link between these modes and solitary solutions (solitary waves and solitons). This link was demonstrated by Scott et al. [Scott 85] where they considered a general system of arbitrary large or infinite degrees of freedom. Their analysis showed that in the limit of weak nonlinearity and when the coupling parameters are at the same order, the system possesses 3 n -1 2 NNMs where the majority are spatially localized. Voluminous studies on mode localization in discrete and continuous periodic nonlinear systems exist in the physics of engineering using appropriate analytical and numerical techniques. In 1955, Fermi, Pasta, Ulam (FPU) and Tsingou [Fermi 55] studied numerically the lack of equipartition of energy in one dimensional anharmonic nonlinear lattice of 64 particles (equation (1.1)), where they found that the solution almost recurred and the energy is localized in the lower modes.

M üi = k 2 (u i-1 -2u i + u i+1 ) + k 3 [(u i+1 -u i ) 2 + (u i -u i-1 ) 2 ]
(1.1) After 10 years, Zabusky et al. [Zabusky 65] studied the continuum limit of the corresponding FPU lattice by expressing the displacement as a continuous functions of position and time and expanding the continuous limit of u ±1 into a Taylor series up to the fourth order to obtain the famous integrable partial differential nonlinear Diederik Korteweg and Gustav de Vries (KDV) equation of the form:

∂ t y + y∂ x y + δ 2 ∂ 3 x y = 0, (1.2)
with ∂ x and ∂ t denoting partial derivatives with respect to the space x and time t, and δ a small parameter. They observed a noticeable property of solitary wave pulse named Soliton that they interact elastically with each other and reappear virtually unaffected in size or shape: "interacting localized pulses do not scatter irreversibly".

Researchers paid particular attention to study nonlinear localization in discrete systems. For one dimensional periodic array of nonlinear oscillators, it was numerically discovered that such system undergo time-periodic and spatially localized solutions (Discrete Breathers DBs or Intrinsic Localized Modes ILMs), whose amplitude decays exponentially in space [Takeno 88]. The existence of such strong, stable large amplitude breathers in a discrete Klein-Gordon model [ Kragh 84] has been discovered and calculated numerically [ Dauxois 92] and a qualitative explanation about their importance was given by [ Dauxois 93]. The existence of such spatially localized solutions was rigorously proved for a broad range of weakly coupled oscillators arrays regardless of the lattice dimension [MacKay 94]. Thermal conductivity in one-and two-dimensional randomly disordered nonlinear lattices has been studied when performing numerical simulations of energy transport in these systems [ Payton 68]. Toda [Toda 89] verified by numerical experiments that introducing a nonlinearity of interaction in periodic lattices, only a little energy partition occurred and the state of the systems return periodically to the initial state. In addition, he proved via numerical experiments that the strength of the energy flow in many cases, larger than that in the linear lattice with same impurity concentration and nonlinearity enhance the energy transmission. Nonlinear mode localization in certain classes of discrete periodic oscillators was discussed by Vakakis et al. [Vakakis 96]. In addition, they confirmed for the first time experimentally the existence of mode localization phenomena in a practical flexible structure of coupled beams with actively induced stiffness nonlinearities. Moreover, they presented an analytical methodology to study standing solitary waves of noncyclic continuous oscillators governed by nonlinear partial differential equations. Nonlinear vibrational localization was firstly reported by Brown [Brown 93] confirming experimentally their presence in a rotating turbine-bladed disk assemblies. The main conclusion of these works is that unlike the linear mode localization, nonlinear mode localization may occur in perfectly symmetric periodic nonlinear lattices, where the only constraint is the presence of weak coupling between oscillators.

Conclusion

Variety of rich dynamical phenomena such as propagation and attenuation zones, localization, solitary waves, energy transfer and frequency-conversion emerged the trend to study wave propagation in periodic nonlinear systems. Pendulums, MEMS, NEMS and granular crystal arrays represent famous examples of periodic nonlinear lattices. It has been shown that spatial localization may occur in linear periodic systems when the coupling between substructures is sufficiently weak and the periodicity is perturbed by small structural mistuning. However, perfectly periodic nonlinear structures may exhibit strongly localized nonlinear normal modes for sufficiently weak substructure coupling.

The principal goal of this thesis consists in exploring the localization phenomena by modeling the collective dynamics of periodic arrays of weakly coupled nonlinear resonators. To do so, we propose the study of the collective dynamics of a periodic array of weakly coupled Duffing-Van Der Pol oscillators under simultaneous primary and external excitations. A convenient choice of physical parameters enabling the creation of linear closed modes has to be considered. In addition, one must combine appropriate quasi-analytical and numerical solving techniques to enable the study of the modal interactions and bifurcation topologies. Therefore, the following chapter presents a review of the analytical and numerical methods which are able to calculate periodic responses of weakly coupled nonlinear systems.

Introduction

Many physical problems in engineering can be modeled as periodic nonlinear systems, which explain the considerable advances that have been made in order to compute their periodic solutions. When such systems are subjected to periodic excitations they react with variety of rich dynamical phenomena, therefore analytical, numerical and combination of both techniques have been developed in order to determine the periodic solutions, stability and bifurcations of such nonlinear systems. They can be generally classified into two main categories, namely, time-domain formulations and frequency-domain methods. The first category includes the direct time integration [Subbaraj 89], shooting methods [Nayfeh 95], Poincaré map method [Parker 89] and the orthogonal collocation methods (Gauss-Legendre collocation point) [ Doedel 07]. These methods consist in transforming the original differential system into an algebraic one using a time integration algorithm, generally limited to a single period. Frequency-domain techniques include the harmonic balance method (HBM) [Krylov 49] and perturbation techniques [Nayfeh 81]. Frequency-domain methods consider periodic solutions with a finite number of harmonics or perturbations by replacing them in the given differential system and solving a resulting set of nonlinear algebraic one. This chapter presents a classification of the analytical and numerical methods to calculate periodic responses of nonlinear systems under periodic excitations. Before doing so, a detailed review on the developed numerical techniques, analytical approaches and combination of both to study the dynamics of periodic nonlinear systems has been reported.

Literature review on methods in periodic nonlinear systems

Several techniques have been developed using numerical methods, analytical approaches and combination of both, in order to study the dynamics of nonlinear periodic structures. In the classical multiple scales perturbation method, Asfar et al.

[ Asfar 83] gave a comprehension on the use of this method to derive coupled mode equation for waveguides with periodic systems, while Nayfeh [Nayfeh 83] analyzed the response of a Multi-degrees-of-freedom (Mdofs) system with quadratic nonlinearities to a harmonic parametric excitation. Manktelow et al. [Manktelow 11] developed a bloch wave-based multiple scales perturbation analysis method showing that in the monoatomic chain, the propagation velocity depend on the wavenumber and amplitude. Later, they determined dispersion frequency shifts for a Bloch wave solution containing two determinant harmonic components, for weakly nonlinear Mdofs periodic structures [Manktelow 14]. In continuous nonlinear periodic systems, wave propagation and dispersion relations were studied using a finite-element discretization of a single unit cell followed by a perturbation analysis finding that the band gaps are amplitude dependent and sensitive to systems parameters [ Manktelow 13].

Chakraborty et al. [Chakraborty 01] applied a perturbation approach to investigate the effects of harmonic wave propagation in an infinite, weakly nonlinear periodic chain. A novel perturbation analysis accompanying numerical simulations was applied to analyze the nonlinear dispersion for weakly nonlinear one-dimensional systems, demonstrating the manner in which nonlinearities may be exploited to achieve amplitude dependent dispersion properties [ Narisetti 10]. A general harmonic balance method was used to analytically investigate wave propagation in strongly nonlinear chain of beads under Hertzian contact [ Narisetti 12].For two-dimensional monoatomic lattice of coupled Duffing masses, Narisetti et al. [Narisetti 11] developed a perturbation approach to evaluate the influence of nonlinearities on the location of band gaps, group velocity magnitudes and the direction of energy propagation. Lazarov et al. [Lazarov 07] used the harmonic balanced method to investigate the influence of the nonlinearities on the filtering properties of the chain around the linear natural frequency of the attached oscillators. A combination of harmonic balance method and multiple scales analysis was used to study wave attenuation in nonlinear periodic structures with weak damping [ Marathe 06]. When dealing with cubically nonlinear oscillatory chain, a nonlinear mapping technique has been used to analyze the modification of the boundary of linear propagation/attenuation zones [ Romeo 06]. Romeo et al. [Romeo 15] applied the same technique using the nonlinear propagation region of chain of oscillators with cubic nonlinearity exhibiting existence solutions to identify regions of existence of discrete breathers and to guide their analysis. Recently, Jothimurugan et al. [Jothimurugan 16] analyzed the effect of the coupling strength on the resonance and anti-resonance frequencies and their response amplitudes, for n coupled Duffing oscillators using numerical integration.

The concept of Nonlinear Normal Mode (NNM) was initiated by Rosenberg [Rosenberg 62] for conservative systems as a synchronous periodic oscillation and symmetric nonlinearities. Vakakis [Vakakis 97] gave a comprehensive review on the NNMs which may provide a valuable theoretical tool for analyzing some specificities of nonlinear systems. Later, an efficient Galerkin projection method was developed, which allows the construction of NNMs that are accurate out to large amplitudes of vibration [Pierre 06]. This approach was extended to the generation of nonlinear modes for variety of applications subjected to external resonance. In structural dynamics, the Asymptotic Numerical Method (ANM) which consists in computing power series expansion of solution branches was introduced by Noureddine Damil and Michel Potier-Ferry [Damil 90] to compute the post buckling behavior of elastic plates and shells. Luongo [Luongo 95] was the first to extend the transfer matrix method to post-buckling problems in periodic nonlinear structures. Later, Luongo et al. [Luongo 06] applied a transfer matrix-perturbation approach to the dynamics of chains of continuous nonlinear sliding beams. An accurate approximate analytical solutions for Mdof coupled Van der Pol-Duffing oscillators has been given by Qian et al. [Qian 10] using the Homotopy Analysis Method (HAM). The precise integration method combined with homotopy perturbation method were used to solve a nonlinear dynamic system with Mdof [ Mei 08]. In the upcoming sections, several numerical and analytical techniques will be detailed starting by the numerical direct time integration method.

Direct time integration method

Direct time integration methods are widely used in the analysis of nonlinear dynamic problems. They presents an easy implementation procedure to compute numerical solutions in addition to their phase portrait and frequency spectra. In direct time integration, the resulting ordinary differential equations will be integrated considering a numerical step-by-step procedure.

Consider, the forced nonlinear system described by the following second order nonlinear equations

ẍn = d 2 x n dt 2 = f ( ẋn , x n , t) for n = 1 . . . N (2.1)
Equation ( 2.1) can be written in the form of two first order differential systems by adding an additional variable as follows

     dx n dt = y n dy n dt = f (y n , x n , t) (2.2)
which is equivalent to

dX dt = g(X, t), (2.3)
where X(t) = (x 1 (t), x 2 (t), . . . , x n (t), y 1 (t), y 2 (t), . . . , y n (t)) t . In order to solve the differential system, ones must find a stationary solution of equation (2.3), then solve it numerically for a set of initial conditions using a time integration method as the Rung-Kutta, Newmark, Adams, and so on until the periodic steady state is reached. Based on that, the solution at t + T is obtained and generated by the following equation X(X 0 , t + T ) = X(X 0 , t), (2.4) where X 0 = X(t = 0) describes the initial conditions and T defines the period of the periodic solution which corresponds to the period orbit in the phase plane.

In order to obtain the nonlinear frequency responses, we proceed iteratively to compute the maximum of amplitude a max for a given frequency Ω in the steadystate regime. However, once a multistability domain is reached, jump phenomena may occurs according to the initial conditions that have been taken into account. Nevertheless, these methods does not give enough information about the global dynamics of the system, particularly the bifurcation behavior where they generally fail to capture unstable solutions. In addition, they are time consuming, for the computation of steady-state solutions especially when considering strongly, high order or weakly damped nonlinear models.

Shooting method

When searching for periodic solutions determining the dynamics of nonlinear systems, particularly when the equations are stiff and the period of the responses are long, shooting methods emerges to replace the direct time integration. It consists in assuming a point directly in the periodic regime and then shoots in time for an assumed period T and checks if the periodicity condition is satisfied leading to highly computational time saving. These trajectories represents solutions obtained by using the boundary conditions which need correction as it is off the mark. Considering the following periodic nonlinear system with

   dX dt = g(X, t) X(X 0 , T ) = X 0 (2.5)
To implement the shooting algorithm, one must seek a solution of the nonlinear differential system such that

Q(X 0 , T ) = X(X 0 , T ) -X 0 = 0, (2.6)
Where Q is the boundary conditions residual function. As X 0 is assumed to approximate the initial solution of the system, one must use an iterative Newton-Raphson correction procedure adopted for the present case for locating periodic solutions. The correction phase may be outlined at iteration k as follows

X k+1 0 = X k 0 + δX, (2.7)
where the correction δX satisfies the following system of linear equations

J(X k 0 , T )δX = -G(X k 0 , T ), (2.8)
with J the Jacobian matrix of Q(X 0 , T ) with respect to X 0 which can be written as

J(X k 0 ) = ∂Q ∂X 0 (X 0 ,T ) (2.9)
The above iterative procedure for computing a periodic solution can be summarized as follows. For a given solution X k 0 , we simultaneously integrate the nonlinear differential system (2.5) and the Jacobian matrix defined in (2.9) and solve the linear system (2.8), for X k 0 at each iteration. Once the correction δX is sufficiently low, the iterations are finished and the periodic solution will be achieved ([ Sundararajan 97]). The relative proposed error condition to check the convergence is obtained by

||δX|| ||X k+1 0 || ≤ ε (2.10)
Since the shooting method requires to provide an initial guess for a solution on the orbit and the period of the orbit, the sensitivity of this method is pronounced when the sought periodic solution is highly unstable and when the different sources of errors are highly controlled.

Asymptotic Numerical Method (ANM)

Continuous or discrete nonlinear problems depending on one or more real parameters λ, can be written in the following form

R(U, λ) = 0, (2.11)
Where R is the vector of n equations ∈ R n and U is the unknown vector ∈ R n . The aim of the family of Asymptotic Numerical Method (ANM) is to numerically find the solution curves (or solution branches) U (λ). These ANM are based on different techniques expanding the solution branches into series expansion which is an efficient continuation technique of type prediction-correction based on computing each branch step by step.

Basic concepts of correction-prediction methods

The general concept of continuation techniques of type of type prediction correction is to consider (U 0 , λ 0 ) as a starting solution point of equation (2.11) and generate a sequence of solution points (U 1 , λ 1 ), (U 2 , λ 2 ), . . . , (U j , λ j ), . . . . This series of points constitute a discrete representation of a branch in the (U, λ) space, where each point satisfies the tolerance criterion

||(U j , λ j )|| < ε (2.12)
when ε is a given positive number.

To describe any solution branch, one must introduce a path parameter denoted by a, which is used to identify points along the solution branch. The transition from the solution point (U j , λ j ) to the new one (U j+1 , λ j+1 ) involves two successive steps: the prediction one which consist in calculating an approximation of the desired quantity followed by the correction step which refines the initial approximation by using the predicted value of the function and another method. Prediction-correction methods differ from each other on many criteria such as the predictor, the corrector, the path parameter and the control of the steps length which means the distance between two successive points. In the following, we briefly remind the most classical techniques: the tangent prediction, Newton correction and the pseudo-arclength parameterization.

6 u i (U j ; 6 j ) "6 (0) "U (0) "a P rediction Correction (U j+1 ; 6 j+1 ) ( 7 U j+1 ; 7 6 j+1 ) Figure 2
.1: Stages of Prediction-Correction method.

Tangent Prediction

The prediction point is defined by

Ū j + 1 = U j + ∆U (0) λj+1 = λ j + ∆λ (0) , (2.13)
where (∆U (0) , ∆λ (0) ) is the prediction increment. This later must be such that ( Ū j + 1, λj + 1) is solution of equation (2.11) as

0 = R(U j + ∆U (0) , λ j + ∆λ (0) ) = R(U j , λ j ) + ∂R ∂U j ∆U (0) + ∂R ∂λ j ∆λ (0) + 1 2 ∂ 2 R ∂U 2 j ∆U (0) ∆U (0) + . . . , (2.14)
where ∂R ∂U j

, ∂R ∂λ j , ∂ 2 R ∂U 2 j
, . . . are the respective partial derivatives of the residual vector R with respect to its arguments U and λ at (U j , λ j ).

In the case where R and U are vectors in R n therefore the derivative ∂R ∂U j becomes a n×n matrix, called the Jacobian matrix or the tangent matrix. Using the fact that R(U j , λ j ) 0, we may neglect the higher order terms in equation (2.14) then the prediction increment becomes solution of the following linear problem which consist of n equation and n + 1 unknowns:

∂R ∂U j ∆U (0) + ∂R ∂λ j ∆λ (0) (2.15)
Solutions of this equation form a straight line which is tangent to the required curve at (U j , λ j ) (see figure 2.1). To define the prediction step length we must add the following condition

||∆U (0) || 2 + (∆λ (0) ) 2 = ∆a 2 (2.16)
Where || • || represent the Euclidean norm and ∆a is a given arc length increment.

Newton correction

If the prediction point R( Ū j+1 , λj+1 ) does not verify the tolerance criterion (2.12) we then proceed to one or more corrections. We suppose that (U (i) , λ (i) ) is an approximation at the iteration (i). We define an approximation at the (i + 1) iteration as follows:

U (i+1) = U (i) + ∆U (i) λ (i+1) = λ i + ∆λ (i) , (2.17)
where (∆U (i) , ∆λ (i) ) is the correction at the (i + 1) iteration. We choose to find the point (U (i+1) , λ (i+1) ) in a manner to satisfy equation (2.11), which after canceling the high order terms gives rise to the following equation

∂R ∂U (i) ∆U (i) + ∂R ∂λ (i) ∆λ (i) = -R(U (i) , λ (i) ). (2.18)
To conclude the problem, we add an additional condition which requests that total increment (prediction plus correction) satisfies the given arc-length condition of type (2.16). In practice, it is convenient to linearize this length condition which can be reflected by the fact that the corrections must be orthogonal to the prediction path:

∆U (i) • ∆U (0) + ∆λ (i) ∆λ (0) = 0. (2.19)
If the prediction point is close to the curve, the norm of the R vector goes to 0 during the correction process and we stop the correction once the tolerance criteria is satisfied.

In summary, the new solution point can be written as: 2.20) where the perdition (∆U (0) , ∆λ (0) ) is solution of equations (2.15), (2.16) and the corrections (∆U (i) , ∆λ (i) ) are solutions of equations (2.18), (2.19). This algorithm is called "Newton-Raphson" (see figure 2.1).

U j+1 = U j + ∆U (0) + ∆U (1) + ∆U (2) + • • • + ∆U (N iter ) λ j+1 = λ j + ∆λ (0) + ∆λ (1) + ∆λ (2) + • • • + ∆λ (N iter ) , (
To move from the solution point (U j , λ i ) to (U j+1 , λ i+1 ), we must solve N iter + 1 linear problem, which can require high computational time if the number of equations n is important. Numerous strategies were proposed in order to save computational time. For instance, when the matrix representing the derivative operator ∂R

∂U

(i)

converges relatively slowly over corrections, it is tempting to keep the same matrix

∂R ∂U (1)
for all correction iterations: which is called the modified-Newton algorithm.

However, this gain of time might cause slower convergence in the correction process and the benefits are not guaranteed. Before performing the prediction and the necessary corrections it would be useful to know in advance the optimal path length, small enough to ensure convergence at a reasonable cost and large enough so that the longest path of the response curve does not get too long. The problem becomes more complicated by choosing different predictor-corrector paths as the optimum may vary widely during the progress: it can be large if the curve is almost rectilinear as it gets strongly smaller in areas of high curvature which appears when modeling instabilities.

To conclude, the process of solving a nonlinear problem requires choices, in particular the choice of the iterative algorithm (Newton-Raphson, modified Newton-Raphson or others) in addition to the choice of the path lengths: which is called computational steering strategy. Although, prediction-correction methods manage to construct solution branches, they fail to give a compromise between reliability and efficiency.

Concepts of the ANM

The idea behind the ANM is to numerically determine a solution of equation (2.11) this time not with a sequence of points but a sequence of branches: branch 1 , branch 2 , . . . , branch j , branch j+1 , . . . as represented in figure 2.2.

Each solution branch is represented in a continuous way by a truncated power series expansion of high order comparing to the path parameter a. Posteriori analysis of the truncated series allows an estimation of the branch end a max , beyond which the approximation is no longer valid. u i

(U j ; 6 j ) (U j+1 ; 6 j+1 ) Branch piece j a 2 [0; a max ] Figure 2.2:
The principle of the ANM.

Perturbation technique: series calculation for algebraic equations system

Considering (U j , λ j ) as a starting solution point, we generate a new branch by seeking its power series expansion with regard to the path parameter a as following .21) This branch must define a solution of the equation (2.11), which is reflected by:

U (a) = U j + ∞ p=1 a p U p = U j + aU 1 + a 2 U 2 + . . . λ(a) = λ j + ∞ p=1 a p λ p = λ j + aλ 1 + a 2 λ 2 + . . . ( 2
0 = R(U (a), λ(a)) = R(U j , λ j ) + ∂R ∂U j (aU 1 + a 2 U 2 + . . . ) + ∂R ∂λ j (aλ 1 + a 2 λ 2 + . . . ) (2.22) + 1 2 ∂ 2 R ∂U 2 j (aU 1 + a 2 U 2 + . . . )(aU 1 + a 2 U 2 + . . . ) + . . . .
Using the fact that R(U j , λ j ) = 0 and ordering the terms according to the growing powers of a, equation (2.22) can be written as

0 = a ∂R ∂U j U 1 + ∂R ∂λ j λ 1 + a 2 ∂R ∂U j U 2 + ∂R ∂λ j λ 2 + 1 2 ∂ 2 R ∂U 2 j U 1 U 1 + 1 2 ∂ 2 R ∂λ 2 j λ 2 + 1 2 ∂ 2 R ∂U ∂λ j U 1 λ 1 + a 3 ∂R ∂U j U 3 + ∂R ∂λ j λ 3 + terms depending on U 1 , U 2 , λ 1 andλ 2 + . . . + a p ∂R ∂U j U p + ∂R ∂λ j λ p + terms depending on U 1 , . . . , U p-1 , λ 1 , . . . , λ p-1 Noted -F nl p . . . Or in a condensed form R(U (a), λ(a)) = aR 1 + a 2 R 2 + • • • = 0. (2.23)
This equation must be verified for all values of a, we deduce an infinity of equations which will allow to successively define terms of the initial series (2.21):

R p = 0 ∀ p ≥ 1.
(2.24)

At each p order, the equation R p = 0 is a linear system in term of U p and λ p which can be written as

∂R ∂U j U p + ∂R ∂λ j λ p = F nl p , (2.25)
where the second member of this equation at the order p, depends only on the lower order terms.

Definition of the path parameter

Equation (2.25) is an ill-posed problem as the number of unknowns is greater than the number of equations. As in the prediction-correction problem, we should define an additional condition at each order, however we use a pseudo arc-length parame- (2.26) which corresponds to the projection of the solution increment on the tangent direction (U 1 , λ 1 ). By replacing equations (2.21) into equation ( 2.26) we find the searched additional condition at each order:

terization a = (U -U j ) • U 1 + (λ -λ j )λ 1 ,
||U 1 || 2 + λ 2 1 = 1 (2.27) U p • U 1 + λ p λ 1 = 0 (2.28)
With a series truncation at the order N , computing the branch j requests calculating the N seconds members F nl p and the computation of N linear problems (2.25) and (2.28).

Steering technique

We know that the representation of an analytical function in a power series can't be valid without identifying the radius of convergence of the series a cr where |a| < a cr (2.29) Calculating numerically a finite number N order of terms of the series (2.21) we obtain an approximate representation of the form of a polynomial function of a

U (a) = U j + N order p=1 a p U p λ(a) = λ j + N order p=1 a p λ p (2.30) (2.31)
We define the validity domain of this approximate solution in a manner to satisfy the following tolerance criteria:

||R(U (a), λ(a))|| < ε, (2.32)
where ε defines the precision parameter. This allows to define an interval [0, a max ] in a manner to satisfy the tolerance criteria for an a in this interval. equation (2.31) for an a ∈ [0, a max ] is the representation of the solution branch j.

Once the approximation (2.31) of the branch j is defined, it becomes easy to continue the process and calculate a new truncated series and its validity domain. The stating point of the branch j+1 can be the final point of previous one: (U j+1 , λ j+1 ). It is possible also to add a correction phase between the computation of both series.

To conclude, the determination of solutions in nonlinear problems with the ANM relies on two techniques: first development of solutions in the form of truncated power series expansion to the given order N order (perturbation technique), then an estimation of the validity domain of this approximation and a simple continuation technique.

ANM with Harmonic Balance Method (HBM)

We focus on treating the following periodically forced non-autonomous system of differential equations

Ẋ = f (t, X, λ) (2.33)
where X is the unknown vector, f is a nonlinear vector valued function, periodic in time t, λ is a real parameter and the dot notation is for the time derivative. When λ varies, the considered system is assumed to have branches of periodic solutions. The HBM and a continuation technique (path following procedure) were implemented to follow these branches and finding them.

Principle of the harmonic balance method

The HBM consists basically in writing an approximation of the periodic solutions of this kind of equation as a truncated Fourier series of order H X

(t) = X 0 + H k=1 X c,k cos(kωt) + H k=1 X s,k sin(kωt) (2.34)
We replace this solution into equation ( 2.33) where we expand the nonlinear vector into Fourier series. After developing and linearizing the trigonometric functions, the differential system is then transformed to an algebraic polynomial one of 2H + 1 vector equations and 2H + 1 unknown X i vectors. Solution branches of the algebraic system are followed by a continuation technique, providing only approximate periodic solutions as we consider a finite number H of harmonics.

The quadratic recast

Since nonlinear periodic systems may include nonlinearities of different types, the expansion of the nonlinear vector into Fourier series can be complicated and heavy despite the help of symbolic software. [Cochelin 09] et al. came with the idea of systematically recast the dynamical system (2.33) in quadratic polynomial form before applying the HBM. The new quadratic equation can be written as following

m( Ẏ ) = c(t, λ) + l(Y ) + q(Y, Y ) (2.35)
Y define the unknown vector of size N eq composed of the original components of X and other added variables which belongs to the quadratic transformation. c is a constant vector with respect to the unknowns Y , l(•) and q(•, •) are respectively the linear and the quadratic vectors with respect to the entry vector Y . In addition, m(•) is a linear vector operator with respect to the unknown Y which is null in the case we are treating an algebraic system.

HBM applied to the quadratic system

The next step is to apply the HBM to the quadratic system (2.35), where the unknown vector Y is decomposed into Fourier series with H harmonics:

Y (t) = Y 0 + H k=1 Y c,k cos(kωt) + H k=1 Y s,k sin(kωt) (2.36)
Then we collect all components of the Fourier series into a

(2H + 1)N eq vector U U = [Y t 0 , Y t c,1 , Y t s,1 , Y t c,2 , Y t s,2 , . . . , Y t c,H , Y t s,H ] t (2.37)
After replacing the expression of the unknown vector decomposition (2.36) into the quadratic system (2.35), regrouping the terms having the same harmonic index and canceling the high orders, we will end up having a large system of (2H + 1)N eq equations

ωM (U ) = C + L(U ) + Q(U, U ) (2.38)
The new operators M (•), C, L(•) and Q(•, •) depend only on the operators of equation (2.35), where λ and ω are respectively the continuation parameter and the angular frequency.

The continuation procedure

In order to solve the obtained algebraic system resulting from the quadratic transformation and the HBM using the ANM, we write the system as following

R(U) = C + L(U) + Q(U, U) -ωM (U) = 0 (2.39) where R ∈ R (2H+1)Neq+1 and U = [U t , λ, ω] t ∈ R (2H+1)Neq+2 .
As shown in the previous section, applying the continuation technique requires using a pseudo-arc length parametrization with respect to λ. Therefore, this latter becomes an unknown as U and ω.

Following this logic, we write c and l(•) as in the following form 2.40) where all these operators are independent of λ. Therefore, the final algebraic system (2.39) becomes

C = c 0 + λc 1 L(•) = l 0 (•) + λl 1 (•), (
R(U) = L0 + L(U) + Q(U, U) (2.41)
where

L0 = C 0 L(U) = L 0 (U ) + λC 1 Q(U, U) = Q(U, U ) + λL 1 (U ) -ωM (U ).
(2.42) L0 is a constant vector, L(•) and Q(•, •) are respectively a linear and bilinear vector valued operators.

The main advantages of the perturbation technique ANM is to find solution branches in the form of power series as a function of the pseudo-arclength path parameter a = (U -U 0 ) t U 1 , where U 0 is assumed to be a known starting regular solution point and U 1 is the tangent vector at U 0

U(a) = U 0 + aU 1 + a 2 U 2 + a 3 U 3 + • • • + a n U n (2.43)
Replacing the series expansion (2.43) into equation ( 2.41) and keeping the power terms up to the N th order leads to

R(a) = R 0 + aR 1 + a 2 R 2 + • • • + a N R N = 0.
(2.44)

Equating each power term R i(1≤i≤N ) of a to zero permits to transform the original nonlinear problem (2.41) into a series of N linear systems of N eq equations in order to successively solve U p . As a result, the problem at the p order can be written as

J U 0 U p + p-1 i=1 Q(U i , U p-i ) = 0, (2.45)
where

J U 0 ∈ R N 2 +1
is the Jacobian matrix of R evaluated at U 0 which is defined by J U 0 U 1 = 0 obtained at the first order. We should note that once each U p (p ∈ [1 . . . N ]) has been found and knowing that the range of a truncated series is generally limited it is defined by the value a max such that

∀a ∈ [0a max ], ||R(U (a))|| ≤ ε r (2.46)
where ε r defines the tolerance parameter as mentioned in the section before.

Implementation in MANLAB

MANLAB is a graphical interactive software implemented in Matlab by Karkar et al.

[ Karkar 12] used for the continuation of solution branches based on the alternative ANM continuation method. The latter is used for the continuation and bifurcation analyses of a dynamical algebraic systems where two stability computation methods were proposed depending on the type of the solution under study and on the selected algorithm: frequency domain or time-domain.

The key point of the ANM implemented in MANLAB lies in the quadratic recast in the form of equation (2.35), then providing MANLAB with the vectors m(•), c, l(•) and q(•, •) where the resulting operators are automatically generated by the software.

Multiple scales perturbation method

Introduction to perturbation techniques

The aim of this section is to introduce and highlight the important role of perturbation (asymptotic) techniques in finding approximate analytical solutions of periodic nonlinear problems. According to these techniques, the perturbation is based on expanding the periodic solution in the form of a power series which is usually represented by the first two terms. The development may be performed with respect to a parameter which naturally occurs in the system or which can be conveniently introduced and called "coordinate perturbation".

Coordinate perturbation

All physical problems involving the operator x(t, ε) can be expressed mathematically as a differential equation of the form L(x, t, ε) = 0 and the boundary condition BC(x, ε) = 0, where t is a scalar or an independent vector variable and ε is a small parameter. Perturbation techniques seek to find an approximated analytical solution of the nonlinear problem by writing it as a power series of small ε as

x(t; ε) = x 0 (x) + εx 1 (t) + ε 2 x 2 (t) + . . . (2.47)
x 0 define the solution of the problem when ε = 0 where x n are the high order terms independent of ε. This expansion is substituted into the mathematical differential equation and the boundary condition. Then, they are expanded and the coefficients of all powers of small ε are collected. Therefore, as the sequences of ε are linearly independent, all there coefficients must be vanished. As a result we obtain a system of simple equations governing x n , which can be successively solved.

The method of multiple scales

The method of multiple scales suggest that the solution x(t; ε) is obtained by considering the expansion representing the response to be a function of multiple explicitly independent variables t, εt, ε 2 t, . . . , as well as ε itself. Hence, in order to achieve a truncated expansion of the solution valid for all t up to O(ε -M ), where M is a positive integer, we introduce different time scales T 0 , T 1 , . . . , T M , according to .48) where T m-1 is slower than T m . Therefore, one assumes that the solution can be represented as an expansion of the form

T m = ε m t for m = 0, 1, 2, . . . ( 2
X(t; ε) = x(T 0 , T 1 , T 2 , . . . , T M ; ε) = M -1 m=0 ε m x m (T 0 , T 1 , T 2 , . . . T M ) + O(εT M ).
(2.49)

The initial ordinary differential system is transformed into a partial differential one, where the derivatives with respect to t become expansions in terms of the partial derivatives with respect to the T n according to

d dt = dT 0 dt ∂ ∂T 0 + dT 1 dt ∂ ∂T 1 + dT 2 dt ∂ ∂T 2 + . . . = ∂ ∂T 0 + ε ∂ ∂T 1 + ε 2 ∂ ∂T 2 + . . . (2.50) d 2 dt 2 = ∂ 2 ∂T 2 0 + 2ε ∂ 2 ∂T 0 ∂T 1 + ε 2 ( ∂ 2 ∂T 2 1 + 2ε ∂ 2 ∂T 0 ∂T 2 ) + . . . (2.51)
The above equations define a version of the method of multiple scales called the many-variable version which was first developed by Sturrock [Sturrock 57] and Nayfeh [Nayfeh 65]. Later, Sturrock and Nayfeh called the technique of expanding the derivatives as well as the dependent variables in powers of small parameters the derivative-expansion method. We consider the forced oscillations of a two coupled mass-spring system as shown in figure 2.3 which can be modeled by the following set of coupled equations of motion

Applications

mü n + c un + ku n -k c (u n+1 -2u n + u n+1 ) + k 3 u 3 n = f cos(Ωt) (2.
m (kg) c (kg.s -1 ) k (N.m -1 ) k c (N.m -1 ) k 3 (N.m -3 ) f (N )
0.05 0.01 70 0.001 50 0.1

Time integration

The ordinary second order differential equations (2.52) can be transformed into a system of 4 first order differential system as follows

                     du 1 dt = v 1 du 2 dt = v 2 dv 1 dt = -µv 1 -ω 2 0 u 1 + d(u 2 -2u 1 ) + αu 3 1 -F cos(Ωt) dv 1 dt = -µv 2 -ω 2 0 u 2 + d(-2u 2 + u 1 ) + αu 3 2 -F cos(Ωt) (2.53)
which is equivalent to

dz dt = g(z, t) (2.54) where z(t) = (u 1 (t), u 2 (t), v 1 (t), v 2 (t)) t .
The above system is solved numerically using the time-difference Runge-Kutta method until we reach a steady-state oscillation regime. The process is described in figure 2.4, where we choose the given frequency Ω = 37.4, once the motion settles down onto a steady oscillation, the maximal amplitude gives a point (Ω, a max ) in the response frequency curve as described in figure 2.5. To obtain the overall frequency response, one must repeat the same procedure for different values of Ω.

ANM with HBM

In order to apply the HBM coupled with the ANM we must transform the considered nonlinear differential system (2.52) into its quadratic form by introducing the following auxiliary variables: w 1 = u 2 1 and w 2 = u 2 2 , the system can then be written as following

0 30 -0,2 -0,1 0 0,1 0,2 t u 1 (t)
u1 = +v 1 u2 = +v 2 v1 =f cos(Ωt) -ω 2 0 u 1 -µv 1 + du 2 -2du 1 -αu 1 w 1 v2 =f cos(Ωt) -ω 2 0 u 2 -µv 2 + du 1 -2du 2 -αu 2 w 2 0 = +w 1 -u 2 1 0 = +w 2 -u 2 2 v2 m( Ẋ) f cos(Ωt) c(Ω) -ω 2 0 u 2 -µv 2 + du 1 -2du 2 l(X) -αu 2 w 2 q(X,X)
The unknown vector is using MANLAB, we must define the functions L0, L and Q which are the actual input where the operators m(•), c 0 , c 1 , l(•), l 1 (•) and q(•, •) are programmed once for all as mentioned before in system (2.42). In MANLAB, the unknowns are assembled in the single vector U = [u 1 , u 2 , v 1 , v 2 , w 1 , w 2 , Ω]. The functions allowing the performance of the continuation solutions in MANLAB are given in detail in the Appendix (A.1).

Z = [u 1 , u 2 , v 1 , v 2 ,
In order to analyze the linear stability of dynamical systems we selected the frequency domain algorithm in the lance.m file. The latter has been used to analyze the stability of fixed points in order to identify the stable multi-mode solution branches. This algorithm relies on the computation of the Jacobian matrix JT, of system (2.52).

JT =     0 0 1 0 0 0 0 1 -ω 2 0 -2d 2 -3αu 2 1 d -µ 0 d -ω 2 0 -2d 2 -3αu 2 2 0 -µ     (2.55)
This analyze uses three additional functions J0, JL and JQ as

JT = J0 + JL + JQ (2.

Multiple scales method

Following the concept of the multiple scales analysis, we seek an approximate solution of the system which can vary slowly in time as the form

u 1 = εu 11 (T 0 , T 2 ) + ε 3 u 13 (T 0 , T 2 ) + . . . u 2 = εu 21 (T 0 , T 2 ) + ε 3 u 23 (T 0 , T 2 ) + . . . (2.57)
Where ε is a small dimensionless parameter related to amplitudes, T 0 = t and T 2 = ε 2 t. Note that, the terms ε 2 u 12 , ε 2 u 22 as well as T 1 = εt are missing from equations (2.57) as the effects of nonlinearity appears at O(ε 3 ) Moreover we set µ and d at the order of ε 2 and f of the order of ε 3 so the effect of damping, linear coupling and primary resonances balance the effect of nonlinearity and appear in the same perturbation equation as the nonlinear Duffing terms. Consequently, the periodic nonlinear system can be written as

ü1 + ω 2 0 u 1 = -ε 2 µ u1 + ε 2 du 2 -2ε 2 du 1 -αu 3 1 + ε 3 f cos(ΩT 0 ) ü2 + ω 2 0 u 2 = -ε 2 µ u2 + ε 2 du 1 -2ε 2 du 2 -αu 3 2 + ε 3 f cos(ΩT 0 ) (2.58)
The derivatives with respect to t of solutions u 1 and u 2 can be written as following by expressing them in the form of expansions in terms of the partial derivatives with respect to the T n denoted by

D n = ∂ ∂Tn u1 = εD 0 u 11 + ε 3 (D 2 u 11 + D 0 u 13 ) + . . . ü1 = εD 2 0 u 11 + ε 3 (2D 0 D 2 u 11 + D 2 0 u 13 ) + . . . (2.59) u2 = εD 0 u 21 + ε 3 (D 2 u 21 + D 0 u 23 ) + . . . ü2 = εD 2 0 u 21 + ε 3 (2D 0 D 2 u 21 + D 2 0 u 23 ) + . . . (2.60)
Substituting solutions (2.57) into the system (2.58) and equating coefficients of like powers of ε, we obtain

order ε D 2 0 u 11 + ω 2 0 u 11 = 0 D 2 0 u 21 + ω 2 0 u 21 = 0 (2.61) order ε 3 D 2 0 u 13 + ω 2 0 u 13 = -D 0 (2D 2 u 11 + µu 11 ) + d(u 21 -2u 11 ) -αu 3 11 + f cos(ΩT 0 ) D 2 0 u 23 + ω 2 0 u 23 = -D 0 (2D 2 u 21 + µu 21 ) + d(u 11 -2u 21 ) -αu 3 21 + f cos(ΩT 0 ) (2.62)
Solutions of equations (2.61) can be written in the following form

u 11 = A 1 (T 2 )e iω 0 T 0 + Ā1 (T 2 )e -iω 0 T 0 = A 1 (T 2 )e iω 0 T 0 + c.c. u 21 = A 2 (T 2 )e iω 0 T 0 + Ā2 (T 2 )e -iω 0 T 0 = A 2 (T 2 )e iω 0 T 0 + c.c. (2.63)
Instead of using the primary excitation Ω, we introduce a detuning parameter σ, which characterize quantitatively its approximation to the natural frequency ω 0 , according to

Ω = ω 0 + εσ (2.64)
where A n are complex unknown functions, Ān are there complex conjugates. We choose c.c. to denote the complex conjugate in general. Substituting solutions (2.63) into the third order equations (2.62) and expressing f cos(ω 0 T 0 +σT 2 ) into its complex form yield

D 2 0 u 13 + ω 2 0 u 13 = -iω 0 (2D 2 A 1 + A 1 µ) + d(A 2 -2A 1 ) -3αA 2 1 Ā1 + 1 2 f e iσT 2 e iω 0 T 0 -αA 3 1 e 3iω 0 T 0 + c.c.
(2.65)

D 2 0 u 23 + ω 2 0 u 23 = -iω 0 (2D 2 A 2 + A 2 µ) + d(A 1 -2A 2 ) -3αA 2 2 Ā2 + 1 2 f e iσT 2 e iω 0 T 0 -αA 3 2 e 3iω 0 T 0 + c.c. (2.66)
Any particular solutions of u 13 and u 23 produces secular terms which are proportional to e iω 0 T 0 unless we vanish them to obtain uniform expansions as

-iω 0 (2A 1 + µA 1 ) + d(A 2 -2A 1 ) -3αA 2 1 Ā1 + 1 2 f e iσT 2 = 0 -iω 0 (2A 2 + µA 2 ) + d(A 1 -2A 2 ) -3αA 2 2 Ā2 + 1 2 f e iσT 2 = 0 (2.67)
where the primes denote the derivatives with respect to T 2 . For solving equations (2.67) we introduce the following polar form of the solutions

A 1 = a 1 e iθ 1 A 2 = a 2 e iθ 2 (2.68)
where a n and θ n are real functions of T 2 . Substituting solutions (2.68) into equations (2.77) and separating the results into real and imaginary parts we obtain

ω 0 (2a 1 -µ)a 1 = da 2 sin(γ 1 -γ 2 ) + 1 2 f sin(γ 1 ) (2.69) 2ω 0 a 1 θ 1 = 2da 1 -da 2 cos(γ 1 -γ 2 ) + 3αa 3 1 - 1 2 f cos(γ 1 ) (2.70) ω 0 (2a 2 -µ)a 2 = da 1 sin(γ 2 -γ 1 ) + 1 2 f sin(γ 2 ) (2.71) 2ω 0 a 2 θ 2 = 2da 2 -da 1 cos(γ 2 -γ 1 ) + 3αa 3 2 - 1 2 f cos(γ 2 ) (2.72)
where

γ 1 = σT 2 -θ 1 γ 2 = σT 2 -θ 2 (2.73)
For the steady state responses a 1 = a 2 = γ 1 = γ 2 = 0. Hence it corresponds to the solutions of the following system

ω 0 µa 1 = da 2 sin(γ 1 -γ 2 ) + 1 2 f sin(γ 1 ) (2.74) 2ω 0 a 1 σ = 2da 1 -da 2 cos(γ 1 -γ 2 ) + 3αa 3 1 - 1 2 f cos(γ 1 ) (2.75) ω 0 µa 2 = da 1 sin(γ 2 -γ 1 ) + 1 2 f sin(γ 2 ) (2.76) 2ω 0 a 2 σ = 2da 2 -da 1 cos(γ 2 -γ 1 ) + 3αa 3 2 - 1 2 f cos(γ 2 ) (2.77)
Solving the previous algebraic system numerically, figure 2.8 (a) shows the variation of identical amplitudes a 1 and a 2 as a function of the detuning parameter σ when the frequency Ω is near the natural frequency ω 0 . Solid lines denote stable branches and dashed lines represent the unstable ones.

Multiple scale method coupled with standing waves

As in the multiple scales analysis, we suppose that the damping is weak so it can be written as c m = εµ and a small coupling parameter kc m = εd with ε a small expansion parameter to be a factor of ε smaller than the Duffing nonlinearity. In addition, to ensure that the external excitation g has the ability to cause such weak oscillations by having it enter the system at the same order as the physical effects, we write the amplitude of the drive as f m = ε 3/2 F . The equations of motion can then be written as following

ün + εµ un + ω 2 0 u n -dε(u n+1 -2u n + u n+1 ) + αu 3 n = ε 3/2 F cos(Ωt) (2.78)
We expand the solutions u 1 and u 2 as a sum of standing wave modes decomposition with slowly varying amplitudes [ Buks 02]. The conditions at the end of the system impose the nature of the standing wave modes. As we considered two additional identical immobilized masses at each end, u 0 and u 3 and assuming

u 0 = u 3 = 0.
(2.79)

The standing wave modes are then

u n = sin(nq m ) with q m = mπ N + 1 , n, m = 1 . . . 2.
(2.80)

To solve the coupled nonlinear equations (2.78) analytically, we use a secular perturbation theory combined with a multiple scales analysis, taking the advantages of the natural separation of time scales, governing the slow dynamics. This approach was described in detail by Lifshitz and Cross [ Lifshitz 10]. We expand the displacement of the n th mass far from its equilibrium, with the expectation that it will act at the order ε 1/2 by writing the solution of the form [Lifshitz 03]:

u n (t) = ε 1/2 N m=1 (A m (T ) sin(nq m )e iωmt + c.c.) + ε 3/2 u (1) n (t) + . . . , n = 1, 2, (2.81)
where c.c. describes the complex conjugate and T = εt is a slow time variable that authorizes the complex amplitude A m (T ) to vary slowly in time. Using the following relation

Ȧn = dA n dt = ε dA n dT ≡ εA n (2.82)
and expending the proposed displacement solutions as

u 1 (t) = ε 1/2 √ 3 2 (A 1 e iω 1 t + A 2 e iω 2 t + c.c.) + ε 3/2 u (1) 1 (t) u 2 (t) = ε 1/2 √ 3 2 (A 1 e iω 1 t -A 2 e iω 2 t + c.c.) + ε 3/2 u (1) 2 (t), (2.83)
their first and second order derivatives can be written as

u1 (t) = ε 1/2 √ 3 2 (iω 1 A 1 + εA 1 )e iω 1 t + (iω 1 A 2 + εA 2 )A 2 e iω 2 t + c.c. + ε 3/2 u(1) 1 (t) u2 (t) = ε 1/2 √ 3 2 (iω 1 A 1 + εA 1 )e iω 1 t -(iω 1 A 2 + εA 2 )A 2 e iω 2 t + c.c. + ε 3/2 u(1) 2 (t) (2.84) ü1 (t) = ε 1/2 √ 3 2 (-ω 2 1 A 1 + 2iω 1 εA 1 )e iω 1 t + (-ω 2 2 A 2 + 2iω 2 εA 2 )A 2 e iω 2 t + c.c. + ε 3/2 ü(1) 1 (t) ü2 (t) = ε 1/2 √ 3 2 (-ω 2 1 A 1 + 2iω 1 εA 1 )e iω 1 t -(-ω 2 2 A 2 + 2iω 2 εA 2 )A 2 e iω 2 t + c.c. + ε 3/2 ü(1) 1 (t) (2.85)
Substituting previous equations with the proposed solutions into the nonlinear system term by term, up to the order ε 1/2 we obtain

order ε 1/2 √ 3 2 A 1 (ω 2 0 -ω 2 1 )e iω 1 t + √ 3 2 A 2 (ω 2 0 -ω 2 2 )e iω 2 t = 0 √ 3 2 A 1 (ω 2 0 -ω 2 1 )e iω 1 t - √ 3 2 A 2 (ω 2 0 -ω 2
2 )e iω 2 t = 0 (2.86) this implies that the resonance frequency of the m th mode is equal to the natural frequency as

ω 1 = ω 2 = ω 0 .
(2.87)

Setting the frequency of the primary resonance to be an amount εΩ D away from the natural frequency as

Ω = ω 0 + εΩ D (2.88)
And writing cos(ω 0 + εΩ D ) into its polar form we obtain at the order ε 3/2 the following two equations

order ε 3/2 ü(1) 1 + ω 2 0 u (1) 1 = √ 3 2 -iω 0 2(A 1 + A 2 ) + µ(A 1 + A 2 ) -d(A 1 + 3A 2 ) + 1 √ 3 F e iΩ D T -α 9 4 |A 1 | 2 A 1 + |A 2 | 2 A 2 + Ā1 A 2 2 + Ā2 A 2 1 + 2|A 2 | 2 A 1 + 2|A 1 | 2 A 2 e iω 0 t -α 3 √ 3 8 A 3 1 + A 3 2 + 3A 1 A 2 2 + 3A 2 A 2 1 e 3iω 0 t + c.c. (2.89) ü(1) 2 + ω 2 0 u (1) 2 = √ 3 2 -iω 0 2(A 1 -A 2 ) + µ(A 1 -A 2 ) -d(A 1 -3A 2 ) + 1 √ 3 F e iΩ D T -α 9 4 |A 1 | 2 A 1 -|A 2 | 2 A 2 + Ā1 A 2 2 -Ā2 A 2 1 -2|A 2 | 2 A 1 + 2|A 1 | 2 A 2 e iω 0 t -α 3 √ 3 8 A 3 1 -A 3 2 + 3A 1 A 2 2 -3A 2 A 2 1 e 3iω 0 t + c.c. (2.90)
On the right hand side of both equations we have secular terms that act to drive the oscillators u

(1) 1 and u

(1) 2 at their resonance frequencies. Therefore, we require them to be vanished so that u

(1) 1 and u

(1) 2 remain finite to obtain

                       Eq1: iω 0 2(A 1 + A 2 ) + µ(A 1 + A 2 ) + d(A 1 + 3A 2 ) +α 9 4 |A 1 | 2 A 1 + |A 2 | 2 A 2 + Ā1 A 2 2 + Ā2 A 2 1 + 2|A 2 | 2 A 1 + 2|A 1 | 2 A 2 = 1 √ 3 F e iΩ D T Eq2: iω 0 2(A 1 -A 2 ) + µ(A 1 -A 2 ) + d(A 1 -3A 2 ) +α 9 4 |A 1 | 2 A 1 -|A 2 | 2 A 2 + Ā1 A 2 2 -Ā2 A 2 1 -2|A 2 | 2 A 1 + 2|A 1 | 2 A 2 = 1 √ 3 F e iΩ D T
(2.91) We may use the orthogonality of the modes by multiplying both equations Eq1 and Eq2 respectively by sin(q 1 ) and sin(q 2 ) and summing them as

Eq1 × sin(q 1 ) + Eq2 × sin(2q 1 ) Eq1 × sin(q 2 ) + Eq2 × sin(2q 2 ) (2.92)
After doing so we find the amplitude equation of both amplitudes

2iω 0 dA 1 dT + iω 0 µA 1 + dA 1 + 9 4 α(|A 1 | 2 A 1 + Ā2 A 2 1 + 2|A 2 | 2 A 1 ) = 1 √ 3 F e iΩ D T 2iω 0 dA 2 dT + iω 0 µA 2 + 3dA 2 + 9 4 α(|A 2 | 2 A 2 + Ā1 A 2 2 + 2|A 1 | 2 A 2 ) = 0 (2.93)
We try the steady state solutions of the form

A 1 (T ) = a 1 e iΩ D T A 2 (T ) = a 2 e iΩ D T (2.94)
Where a 1 and a 2 are complex variables in function of T . By replacing the steady state solutions (2.94) into the amplitude equations (2.93), we obtain equations for the time-independent complex amplitudes a 1 and a 2

ω 0 (iµ -2Ω D )a 1 + da 1 + 9 4 α(|a 1 | 2 a 1 + 2|a 2 | 2 a 1 + ā2 a 2 1 ) = F √ 3 ω 0 (iµ -2Ω D )a 2 + 3da 2 + 9 4 α(|a 2 | 2 a 2 + 2|a 1 | 2 a 2 + ā1 a 2 2 ) = 0 (2.95)
In figure 2.9, Mathematica was used to plot the solutions for the response intensity of two forced coupled mass-spring system as a function of the detuning parameter Ω D with the stability analysis. Blue and red colors indicate respectively stable and unstable solutions. Two additional stable solution branches appeared. However, when returning to the proposed solution and plotting the displacement u 1 as a function of the frequency Ω (see figure 2.10), these branches join the resonant and non-resonant branches. In this case and for the chosen physical parameters, the impact of the additional multimodal solutions on the frequency response is not noticeable. In Chapter 4, we present the physical example of weakly coupled pendulums array, for which the additional multimodal solutions are disconnected from the resonant and non-resonant solution branches. The ANM with the HBM represent a powerful analysis tool for nonlinear vibrating systems, however when dealing with weakly coupled nonlinear oscillators it encounters several problems. We can summarize these problems by: (i) convergence problems for dissipative systems, (ii) it fails to capture all branches when they are separated, (iii) enable to provide the solution type in terms of modal interactions and (iv) when having stable solution branches which are completely disconnected from the whole response curve, on must define multiple initial conditions and collected the responses to form the whole response, especially for high number of dofs as in Chapter 3. As the additional multimodal solution branches are not distinguished visually from the resonant and non-resonant branches, we proceed in the following section to use the basins of attraction as a powerful tool to investigate their dynamics.

Basins of attraction

The basins of attraction can be used for qualitative as well as quantitative analysis of the nonlinear dynamics robustness. In a nonlinear nanomechanical resonator, Kozinsky et al. [ Kozinsky 07] experimentally probe the basins of attraction of two fixed points. Moreover, Sliwa et al. [Śliwa 12] investigated the basins of attraction of two coupled Kerr oscillators. Furthermore, Ruzziconi et al. [Ruzziconi 13] studied frequency response curves, behavior charts and attractor-basins phase portraits of a considered NEMS constituted by an electrically actuated carbon nanotube.

In this section, the basins of attraction are used to investigate qualitatively the trajectories of the system response, the robustness of the attractors and their practical implications, for the case of two coupled mass-spring oscillators under primary resonance. They are performed in a classical way where the robustness is only related to the global size of the attractor by solving the nonlinear differential system (2.93) after replacing the steady state solutions (2.94) and the complex values a n into α n + iβ n .

Although the basins of attraction are usually plotted in the phase plane (u n , un ), we choose to represent them for a 1 in Nyquist plane (α 1 , β 1 ) while considering α 1 = α 2 and β 1 = β 2 for several frequencies as shown in figure 2.11. It shows the evolution of the basins of attraction of a 1 between resonant and nonresonant branches, with additional distribution of the multimodal nonresonant solution branch for Ω D = 0.6 and Ω D = 0.7. Although, the distribution of the additional multimodal solutions is narrow in the considered case, Chapter 3 present two different configurations, where the distribution of the basins of attraction of the additional multimodal solutions is important.

Conclusion

In this chapter, we presented a batch of computational methods for solving systems of nonlinear second order differential equations. The classical direct time integration methods define an easy numerical way to solve such system, although they are not able to describe the whole dynamics of the system. They generally fail to reach unstable solutions as the bifurcation points. In case the system is weakly damped, the shooting method came to replace the direct integration ones as the time to reach periodic solution could be extremely high. Shooting methods consist in stating by a solution in the steady state regime ignoring the transient regime while taking into account the boundary conditions. Programing such method is more complicated than the direct ones, however they reduce significantly the computational time. Among the numerical frequency domain methods, the HBM coupled with the ANM is a powerful tool to solve differential nonlinear systems, with the study of stability and bifurcation. However, when dealing with large number of weakly coupled nonlinear oscillators, this method encounters difficulties in convergence, in capturing all solutions specially when having disconnected branches from the entire response. The multiple scales analysis is a very useful method to analytically study weakly nonlinear problems. In addition, combining this method with the standing wave decomposition, enables the capture of additional branches resulting from modal interactions, the identification of their nature and the study of their robustness using the basins of attraction. A comparison between direct time integration method based on the Runge-Kutta algorithm, the ANM with HBM method and the multiple scales method coupled with standing wave modes decomposition is represented in figure 2.12. An excellent agreement between these methods is remarkable while, highlighting the efficiency of the analytical multiple scales method in displaying additional information's in terms of model interactions and bifurcation topology.

Introduction

In this chapter, the collective dynamics of a periodic structure of coupled Duffing-Van Der Pol (D-VDP) oscillators is investigated under simultaneous external and parametric excitations. An analytico-computational model based on a perturbation technique, combined with standing wave decomposition and the asymptotic numerical method is developed for a finite number of coupled oscillators. The frequency responses and the basins of attraction are analyzed for the case of small size arrays, demonstrating the importance of the multi-mode solutions and the robustness in terms of the distribution area of their basins of attraction.

Motivations

Several researches were devoted to the nonlinear dynamics of coupled structures in the presence of a single excitation. For instance, Nayfeh et al. [Nayfeh 83] studied the response of Mdofs systems with quadratic nonlinearities to a harmonic parametric resonance. Perkins et al. [Hikihara 12] illustrated the beneficial effects that noise can produce on the responses of an array of coupled nonlinear oscillators externally excited. The dynamic behavior of a microbeam array subjected to parametric actuation at low, medium and large dc-voltages has been investigated by Gutschmidt et al. [Gutschmidt 12] using a continuum model. To possess rich dynamic behavior unreachable with uncoupled NEMS resonators, Buks and Rukes [ Buks 02] studied the mechanical characteristics of a parametrically excited beam resonators array. Same results were proved qualitatively by Lifshitz et al. [Lifshitz 09] involving a parametrically excited micromechanical resonators array. In addition, Lifshitz et al. [Lifshitz 03] investigated the collective dynamics in periodic coupled resonators in presence of nonlinearities up to the third order under parametric excitation. However few studies have been devoted to simultaneous resonances [Plaut 86,Kacem 11a, Kacem 15,Nguyen 13] and they are mostly limited to single or few dofs systems.

In this chapter, we developed a generic discrete model for the collective dynamics of periodic structures of coupled D-VDP oscillators, with nonlinearities up to the fifth order under simultaneous primary and parametric resonances, which is suitable for several physical applications. The motivations behind considering the nonlinear quintic coupling are detailed in Chapter 5. The main goal is to understand how nonlinearities influence the presence of simultaneous external and parametric excitations and how they may be used to enhance and control modal interactions and bifurcation topology transfer between coupled oscillators. The analytico-computational model is based on the method of multiple scales coupled with standing wave modal decomposition, transforming the nonlinear differential system into a set of coupled complex algebraic equations which are numerically solved using the ANM enabling the construction of resonance curves for a large number of dof. The cases of small size resonator arrays have been analyzed in the frequency domain and it is shown that the multi-mode solutions are stable over a wide frequency-range for a particular set of design parameters. The complexity and the multivaludness of the response were illustrated by a detailed study of its basins of attraction which proves theoretically the robustness of the multi-mode branches.

One-dimensional array of coupled nonlinear oscillators

Equations of motion

The proposed model involves a finite degree of linearly and nonlinearly coupled D-VDP oscillators, under primary and parametric excitations, as shown in Figure 3.1. The corresponding set of coupled Equations of Motion (EoM) can be written in the following form:

Ψ,Φ K d Λ u n-1 C, Δ K, Ξ G(t) Ψ, Φ K d Λ u n K, Ξ G(t) Ψ,Φ K d Λ u n+1 K, Ξ G(t) M H(t) H(t) H(t) M M M M C, Δ C, Δ
M ¨ u n -K d u n+1 -2 u n + u n-1 + H cos(ω p t) u n + C ˙ u n + K u n + Ξ u 3 n +Ψ ( u n -u n+1 ) 3 + ( u n -u n-1 ) 3 + Φ ( u n -u n+1 ) 5 + ( u n -u n-1 ) 5 + ∆ u 2 n ˙ u n +Λ ( u n -u n+1 ) 2 ( ˙ u n -˙ u n+1 ) + ( u n -u n-1 ) 2 ( ˙ u n -˙ u n-1 ) = G cos(ω e t) (3.1)
where u n describes the deviation of the n th resonator from its equilibrium n = 0, • • • , (N + 1), with fixed boundary conditions u 0 = u N +1 = 0. ω p and ω e are respectively the frequencies of the parametric and the external excitations of each oscillator. For the parametric excitation, we consider the largest excitation effect that occurs when the pump frequency ω p is close to twice the resonant frequency of the oscillator. Therefore, we take ω p = 2ω 0 (1 + εΩ D ), where ω 0 is the natural frequency, Ω D is a detuning parameter and ε is a small non-dimensional expansion parameter. In addition we set the external driving frequency of half the pumping frequency as ω e = ωp 2 = ω 0 (1+εΩ D ). M is its effective mass, K = M ω 2 0 is its effective spring constant, Ξ and Ψ are respectively the nonlinear stiffness and the coupling Duffing parameters, Φ the nonlinear quintic coupling parameter, C = M ω 0 Q is the linear damping (Q is the quality factor), K d = ΓK is the coupling spring constant, ∆ and Λ are respectively the VDP damping and the nonlinear dissipative coupling. H and G are the parametric and external excitation amplitudes respectively.

Normalization

For convenience and equation simplicity, we introduce the non-dimensional variables:

u n = u n u D , t = tω 0 (3.2)
where u D = G Cω 0 is the dynamic displacement of the associated linear system while neglecting the linear coupling.

Substituting the non-dimensional variables (3.2) into the EoM (3.1), we obtain after dividing by

M Gω 0 C ün -Γ(u n+1 -2u n + u n-1 ) + 1 Q un + u n + H K cos[2(1 + εΩ D )t]u n + ΞG 2 M C 2 ω 4 0 u 3 n + ΨG 2 M C 2 ω 4 0 [(u n -u n+1 ) 3 + (u n -u n-1 ) 3 ] + ΦG 4 M C 4 ω 6 0 [(u n -u n+1 ) 5 + (u n -u n-1 ) 5 ] + ∆G 2 M C 2 ω 3 0 u 2 n un + ΛG 2 M C 2 ω 3 0 [(u n -u n+1 ) 2 ( un -un+1 ) + (u n -u n-1 ) 2 ( un -un-1 )] = C M ω 0 cos[(1 + εΩ D )t] (3.3)
The aim of this study is to present the influence of the modes localization on the collective dynamics of an array of coupled nonlinear oscillators. Therefore, a study of the linear natural frequency is required.

Linear study

Equation (3.3) describes a non-dimensional system of coupled D-VDP oscillators subjected to simultaneous primary and parametric excitations. To determine the natural frequencies and their eigenvectors, equation (3.3) can be written in matrix form as

M Ü + C U + K L U + H(t) U Linear part + F N L ( U , U ) Nonlinear part = G(t), (3.4)
with the displacement vector U = [u 1 , u 2 , . . . , u N ] T , the excitation vector

G(t) = C M ω 0 cos[(1 + εΩ D )t][1, . . . , 1] T , M = diag(1, 1, . . . , 1), C = 1 Q * diag(1, 1, . . . , 1) , F N L ( U , U
) is a nonlinear restoring force vector and

K L =         1 + 2Γ -Γ -Γ 1 + 2Γ -Γ 0 . . . . . . . . . 0 -Γ 1 + 2Γ -Γ -Γ 1 + 2Γ        
The eigenvalue problem associated to the linear system can be written as:

[K L -ω 2 i M]u i = 0, i = 1, 2, . . . , N. (3.5)
The dimensionless eigenvalues ω i are represented in the following form for few coupled oscillators:

N = 1 : ω 1 = √ 1 + 2Γ N = 2 : ω 1 = √ 1 + Γ ω 2 = √ 1 + 3Γ N = 3 : ω 1 = 1 - √ 2Γ + 2Γ ω 2 = √ 1 + 2Γ ω 3 = 1 + √ 2Γ + 2Γ
We may express all normal frequencies relative to the same non-dimensional natural frequency which is 1, so that

ω n = 1 + λ n Γ (n = 1, . . . , N ), (3.6)
where λ n is the coefficient depending on the mode of vibration (As an example for N = 2, λ 1 = 1 and λ 2 = 3). We suppose that the oscillators are weakly coupled; so that K d << K. Consequently, Γ << 1 and:

ω n ≈ 1 + ∆ n (∆ n = λ n 2 Γ << 1 for n = 1, 2, . . . ) (3.7)
This assumption leads to the creation of linear closed modes which permit to study the effect of the mode localization on the collective dynamics. We proceed in the following subsection to solve these coupled nonlinear differential equations using the multiple time scales method combined with standing waves decomposition. This method has been detailed in the previous Chapter 2, section 2.8.4 for two coupled oscillators.

Derivation of the amplitude equation

The linear damping coefficient is the physical parameter that allows us to use the multiple times scales perturbation approach. We therefore assume it to be small, by expressing it as 1 Q = ε 2 η treating ε as a small expansion parameter. This limit is well verified for targeted applications where Q ≥ 50. We also require that the Duffing and Van-Der-Pol nonlinearities to be a factor of ε smaller than the linear force, or equivalently by taking the leading order in u n to be of the order of ε 1 2 , and expressing Γ = 1 2 ε 2 γ. In addition, we choose to take the parametric excitation amplitude to scale as the damping, i.e. we set H K = ε 2 h. To ensure that the external excitation g has the ability to cause such weak oscillations by having it enter the system at the same order as the physical effects, we write the amplitude of the drive as C

M ω 0 = ε 5/2 g. On the other hand, as u n is proportional to GQ for a regular linear resonance, with u n to be of order ε 1 2 and Q, of order ε -2 thus G has to be of order ε 5/2 . To insure that the nonlinear cubic terms have the ability to cause weak oscillations, we write them with the following corresponding scaling so they inter into the system at the same range as the physical impacts

ΨG 2 M C 2 ω 4 0 = εψ, ΞG 2 M C 2 ω 4 0 = εξ, ∆G 2 M C 2 ω 3 0 = εδ, ΛG 2 M C 2 ω 3 0 = ελ and ΦG 4 M C 4 ω 6 0 = φ. (3.8)
The normalized, scaled differential system becomes

ün -ε 2 1 2 γ(u n+1 -2u n + u n-1 ) + ε 2 h cos[2(1 + εΩ D )t]u n + ε 2 η un + u n + εξu 3 n +εδu 2 n un + εψ[(u n -u n+1 ) 3 + (u n -u n-1 ) 3 ] + φ[(u n -u n+1 ) 5 + (u n -u n-1 ) 5 ] +ελ[(u n -u n+1 ) 2 ( un -un+1 ) + (u n -u n-1 ) 2 ( un -un-1 )] = ε 5 2 g cos[(1 + εΩ D )t]
(3.9)

To calculate the nonlinear response of each simultaneously excited oscillator, we use the method of multiple scales. We express the solution as a sum of standing wave modes with slowly varying amplitudes. As we consider two additional fixed masses labeled 0 and N + 1 at each end u 0 = u n+1 = 0, the standing wave modes are then:

u n = sin(nq m ) with q m = nm N + 1 , m = 1 . . . N (3.10)
To investigate the EoM analytically, we use a multiple scale analysis combined with secular perturbation theory, with the expectation that the displacement of the n th mass far from its equilibrium will act at the order of ε 1 2 , by writing the solution of the form:

u n (t) = ε 1 2 N m=1 (A m (T ) sin(nq m )e iωmt + c.c.) + ε 5/2 u (1) n (t) + • • • , n = 1, ..., N, (3.11)
where T = ε 2 t is a slow time variable, that authorizes the complex amplitude A m (T ) to vary slowly in time and c.c. describes the complex conjugate. (A.2)

Replacing the proposed solution (3.11) into the normalization EoM (3.9) up to the order ε 3/2 as listed in Appendix (A.2). At the order of ε 1/2 , we get the linear dispersion relation, given by

ω 2 m = 1, m = 1 . . . N.
(3.12)

Therefore, at the order of ε To apply the solvability condition allowing to get the system for the slowly changing magnitudes A m (T ), we require to vanish all secular terms that act to drive the coupled oscillators u

(1) n at their resonance frequencies. Thus, we obtain the equations for the slowly varying amplitudes A m (T ), where we benefit from the orthogonality of the modes by multiplying the secular terms by sin(nq m ) and summing over n. Therefore, the m th secular term coefficient can be written as

2i dA m dT + 2γ sin 2 ( q m 2 )A m + h 2 A * m e 2iΩ D T + iηA m + 1 4 j,k,l (iδ + 3ξ)A j A k A * l ∆ (1) jkl;m +4 sin( q m 2 ) j,k,l sin( q j 2 ) sin( q k 2 ) sin( q l 2 )(3ψ + iλ)A j A k A * l ∆ (2) jkl;m +40i φ sin( q m 2 ) j,k,l,o,p sin( q j 2 ) sin( q k 2 ) sin( q l 2 ) sin( q o 2 ) sin( q p 2 )A j A k A l A * o A * p ∆ (3) jklop;m = g (N + 1) e iΩ D T N n=1 sin(nq m ), (3.14)
Ignoring initial transients, and assuming that the nonlinear terms in the equation are sufficient to saturate the growth of the instability, we try a steady-state solution of the form

A m = a m e iΩ D T (3.15)
Substituting equation (3.15) into the differential system (3.14) of amplitude, we obtain the required equation for the fixed complex amplitudes a m .

(iη -2Ω D )a m + 2γ sin 2 ( q m 2 )a m + h 2 a * m + 1 4 j,k,l (iδ + 3ξ)a k a j a * l ∆ (1) jkl;m +4 sin( q m 2 ) j,k,l sin( q j 2 ) sin( q k 2 ) sin( q l 2 )(3ψ + iλ)a k a j a * l ∆ (2) jkl;m +40i φ sin( q m 2 ) j,k,l,o,p sin( q j 2 ) sin( q k 2 ) sin( q l 2 ) sin( q o 2 ) sin( q p 2 )a k a j a l a * o a * p ∆ (3) jklop;m = g (N + 1) N n=1
sin(nq m ), (3.16) Before starting our investigations, we should note that the second member of equations (3.16) is proportional to the sum of standing wave modes, which is null for all even mode modes regardless of the number of coupled pendulums considered. Therefore all even modes a 2n are not excited after modes projection, as they possess null trivial solutions. In addition, if we consider that whenever for a given mode m, ∆

;j = ∆ (2) mmm;j = ∆ (1) mmm
mmmmm;j = 0 for all j = m, then a Single Mode (SM) solution branch can exist with a m = 0 and a j = 0. Thus, the only SM solution can exist for the first mode a 1 in the case of two coupled oscillators, where its amplitude takes the form of a single driven Duffing oscillator response.

Numerical solutions

Equation (3.16) is a complex algebraic system, with a large number of variables and nonlinearities up to the fifth order. Analytical solutions for this type of equations are either too large or simply do not exist. In this case, the only possibility is to solve it numerically, using a conventional method. Mathematica can be used to solve the system for two coupled oscillators, including stability analysis. This method is based on prediction-correction algorithms, such as the Newton-Raphson scheme which is the most popular way to solve a nonlinear structural problem. In general such algorithms are successful for determining nonlinear solution branches in a stepby-step manner, with a load control, a displacement control or an arc-length control but they have two disadvantages: the first one, is that they are time-consuming comparing to a linear problem and the second one is the difficulty of automatization of the continuation process.

Therefore, a graphical interactive software named ManLab [ MAN 10] has been used for the continuation of branches of solutions of equations (3.16) by an alternative method, which is called the Asymptotic Numerical Method (ANM) [Azrar 93,Cochelin 94] is to write the algebraic equations of the form R(U ) = 0, where U is a vector of n + 1 unknowns and R a vector of n smooth equations that must be analytical. The ANM is a perturbation technique which consists in expanding the unknown vector U as a formal power series of a path parameter. It presents several advantages: it provides continuous solutions, the continuation is very robust, and the control of the step length is automatic and always optimal. ManLab provides linear stability analysis for equilibrium points of dynamical systems with an automatic bifurcation detection. This method and its application to our nonlinear differential system (3.14) is detailed here after.

Cartesian transformation

To use the ANM, equation (3.14) is transformed to its cartesian form by defining the amplitude as A m = (α m + iβ m )e iΩ D T . As a result, we obtain the following general equations, for which the unknowns are real: 3.17) and

α m (T ) = - η 2 α m + Ω D β m + h 4 β m -γ sin[ q m 2 ] 2 β m α m (T ) - 1 8 j,k,l δ(α j α k α l + α j β k β l ) + 3ξ(α j α k β l + β j β k β l ) ∆ (1) jkl;m α m (T ) -2 sin ( q m 2 ) j,k,l sin( q j 2 ) sin( q k 2 ) sin( q l 2 ) λ(α j α k α l + α j β k β l ) α m (T ) + 3ψ(α j α k β l + β j β k β l ) ∆ (2) jkl;m α m (T ) -20 φ sin ( q m 2 ) j,k,l,o,p sin( q j 2 ) sin( q k 2 ) sin( q l 2 ) sin( q o 2 ) sin( q p 2 ) α j α k α l α o β p α m (T ) + 2 α j α k β l β o β p + β j β k β l β o β p ∆ (3) jklop;m (
β m (T ) = - g 2(N + 1) N n=1 sin(nq m ) - η 2 β m -Ω D α m + h 4 α m + γ sin[ q m 2 ] 2 α m α m (T ) - 1 8 j,k,l 3ξ(α j α k α l + α j β k β l ) -δ(α j α k β l + β j β k β l ) ∆ (1) jkl;m α m (T ) + 2 sin ( q m 2 ) j,k,l sin( q j 2 ) sin( q k 2 ) sin( q l 2 ) 3ψ(α j α k α l + α j β k β l ) α m (T ) -λ(α j α k β l + β j β k β l ) ∆ (2) jkl;m α m (T ) -20 φ sin ( q m 2 ) j,k,l,o,p sin( q j 2 ) sin( q k 2 ) sin( q l 2 ) sin( q o 2 ) sin( q p 2 ) α j α k α l α o α p α m (T ) + 2 α j α k α l β o β p + α j β k β l β o β p ∆ (3) jklop;m (3.18)
The steady-state motions occur when α m = β m = 0, which corresponds to the singular points of equations (3.17) and (3.18). With these two equations, the stability of fixed points is easy to implement.

Quadratic recast

The key point of the ANM lies in the quadratic recast of equations (3.17) and (3.18) by introducing the following set of auxiliary variables,

α 2 i = c i size N α 4 i = c 2 i = e i -- β 2 i = d i -- β 4 i = d 2 i = f i -- α i α j = g i,j for i = j -N (N -1) 2 β i β j = h i,j for i = j -- c i c j = r i,j for i = j -- d i d j = s i,j for i = j -- c i d j = t i,j -N 2 c k g i,j = v k;i,j -N 2 (N -1) 2 d k g i,j = w k;i,j -- c k h i,j = x k;i,j -- d k h i,j = y k;i,j -- g i,j g k,l = z1 i,j;k,l -N C 4 h i,j h k,l = z2 i,j;k,l -- g i,j h k,l = z3 i,j;k,l -N 2 (N -1) 2 4 (3.19)
These transformations lead to the following quadratic system

R(U ) = L 0 + L(U ) + Q(U , U ) = 0 (3.20)
where R is a vector of 1 3 N 4 + N 3 + 13 6 N 2 + 7 2 N equations, L 0 is a constant vector, L(•) and Q(•) are respectively the linear and quadratic operators with respect to U . U = (α, β, c, d, e, f , g, h, r, s, t, v, w, x, y, z1, z1, z3, Ω D ) T is the vector of

1 3 N 4 + N 3 + 13
6 N 2 + 7 2 N + 1 unknowns, in which α, β, c, d, e, f , g, h, r, s, t, v, w, x, y, z1 , z2 and z3 are vectors, as

α = {α 1 , α 2 , • • • , α i , • • • , α N } β = {β 1 , β 2 , • • • , β i , • • • , β N } c = {c 1 , c 2 , • • • , c i , • • • , c N } d = {d 1 , d 2 , • • • , d i , • • • , d N } e = {e 1 , e 2 , • • • , e i , • • • , e N } f = {f 1 , f 2 , • • • , f i , • • • , f N } g = {g 1,2 , g 1,3 , • • • , g 1,N , • • • , g i,i+1 , • • • , g i,N , • • • • • • , g N -1,N } h = {h 1,2 , h 1,3 , • • • , h 1,N , • • • , h i,i+1 , • • • , h i,N , • • • • • • , h N -1,N } r = {r 1,2 , r 1,3 , • • • , r 1,N , • • • , r i,i+1 , • • • , r i,N , • • • • • • , r N -1,N } s = {s 1,2 , s 1,3 , • • • , s 1,N , • • • , s i,i+1 , • • • , s i,N , • • • • • • , s N -1,N } t = {t 1,1 , t 1,2 , • • • , t 1,N , • • • , t i,1 , • • • , t i,N , • • • , t N -1,1 , • • • , t N,N } v = {v 1;1,2 , v 1;1,3 , • • • , v 1;i,j , • • • , v k;1,2 , v k;1,3 , • • • , v k;i,j , • • • } w = {w 1;1,2 , w 1;1,3 , • • • , w 1;i,j , • • • , w k;1,2 , w k;1,3 , • • • , w k;i,j , • • • } x = {x 1;1,2 , x 1;1,3 , • • • , x 1;i,j , • • • , x k;1,2 , x k;1,3 , • • • , x k;i,j , • • • } y = {y 1;1,2 , y 1;1,3 , • • • , y 1;i,j , • • • , y k;1,2 , y k;1,3 , • • • , y k;i,j , • • • } z1 = {z1 1,2;1,2 , z1 1,2;1,3 , • • • , z1 1;i,j , • • • , z1 k,l;1,2 , z1 k,l;1,3 , • • • , z1 k,l;i,j , • • • } z2 = {z2 1,2;1,2 , z2 1,2;1,3 , • • • , z2 1;i,j , • • • , z2 k,l;1,2 , z2 k,l;1,3 , • • • , z2 k,l;i,j , • • • } z3 = {z3 1,2;1,2 , z3 1,2;1,3 , • • • , z3 1;i,j , • • • , z3 k,l;1,2 , z3 k,l;1,3 , • • • , z3 k,l;i,j , • • • } (3.21)
Substituting the auxiliary variables given in equation (3.19) into the real differential equations and replacing the nonlinear cubic and quintic terms into there quadratic form as given in Appendix A.4, systems (3.17) (3.18) can be rewritten as

α m = 0 -η 2 α m + h 4 β m -γ sin[ qm 2 ] 2 β m +Ω D β m + Q αm β m = -g 2(N +1) N n=1 sin(nq m ) -η 2 β m + h 4 α m + γ sin[ qm 2 ] 2 α m -Ω D α m + Q βm 0 = 0 + c i -a 2 i 0 = 0 + e i -c 2 i 0 = 0 + d i -b 2 i 0 = 0 + f i -d 2 i 0 = 0 + g i,j -a i a j 0 = 0 + h i,j -b i b j 0 = 0 + r i,j -c i c j 0 = 0 + s i,j -d i d j 0 = 0 + t i,j -c i d j 0 = 0 + v k;i,j -c k g i,j 0 = 0 + w k;i,j -d k g i,j 0 = 0 + x k;i,j -c k h i,j 0 = 0 + y k;i,j -d k h i,j 0 = 0 + z1 i,j;k,l -g i,j g k,l 0 = 0 + z2 i,j;k,l -h i,j h k,l 0 = 0 + z3 i,j;k,l -g i,j h k,l M (U ) R(U ) -g2(N + 1)dan = 1 N L0 η2β m + h/4α m + γ sin[q m /2] 2 α m L(U ) -Ω D α m + Q m da Q(U,U ) (3.22)
where Q αm and Q βm are respectively the quadratic transformations of nonlinear cubic and quintic terms in equations (3.17) and (3.18).

Stability analysis

Different algorithms are implemented in ManLab in order to analyze the linear stability of dynamical systems. Two stability computation methods are proposed depending on the type of the solution under study and on the selected algorithm: frequency domain or time-domain. The latter has been used to analyze the stability of fixed points in order to identify the stable multi-mode solution branches. This algorithm relies on the computation of the Jacobian matrix JT , of equations (3.17) and (3.18). The linear stability analysis consists in computing the eigenvalues of the Jacobian matrix at each point of analysis. If any of the eigenvalues has a positive real part, then the current point is unstable. When following a branch that is, at first stable, a bifurcation can be detected when one of the eigenvalues crosses the imaginary axis. This analyze uses three additional functions J0, JL and JQ as

JT = J0 + JL + JQ (3.23)
where J0 is a constant matrix, JL is a linear operator and JQ a quadratic operator on the variables given in U stab = {α, β, c, d, g, h, b, Ω D } which is the variables vector of size N 2 + 3N + 1, where b is a vector of N 2 b i,j = α i β j variables. From a practical point of view, one can easily recast the equations of motion of nonlinear periodic structures with respect to the proposed model and hence one can use the procedure detailed above as a robust solving tool.

Examples and discussions

In order to target various periodic structure-based systems, the numerical simulations have been performed with two sets of design parameters listed in Table 3.1. The first one contains the Duffing coupling term and the nonlinear damping for each oscillator while the second one involves the VDP coupling and the nonlinear cubic stiffness. Indeed, the two configurations can illustrate respectively an array of coupled levitated magnets and a periodic structure of coupled micro or nano-beams.

Moreover, the external excitation G has been chosen sufficiently high in order to generate nonlinear frequency responses for which the displacement at resonance is beyond the critical Duffing amplitude, while the value of the parametric excitation H has been adjusted in order to increase the interactions between both resonances (primary and parametric) and enrich the resulting collective dynamics.

K(N.m -1 ) 40 40 ∆(Kg.m -2 .s -1 ) 50 0 Ψ(N.m -3 ) 5000 0 Ξ(N.m -3 ) 0 800 Λ(Kg.m -2 .s -1 ) 0 6 Φ(N.m -5 ) 0 0 Γ 0.001 0.001

Two coupled nonlinear resonators

As a first step, we started by solving the coupled equations (3.16) with N = 2 for the first design, in order to provide some qualitative and quantitative explanations dealing with the collective nonlinear dynamics of small arrays of periodic structures which can be extended for large periodic lattices. For two coupled oscillators, we have:

q 1 = π 3 and q 2 = 2π 3 .
The coupled algebraic equations have been solved numerically using Mathematica, for several values of Ω D inside the frequency range where the whole dynamic response is represented and the stability of the different branches have been performed based on the eigenvalues of the Jacobian matrix of the differential system (3.14) computed numerically for each point. In this subsection, we are interested in the collective nonlinear dynamics of the coupled oscillators under primary resonance (G = 0.08). By solving numerically the resulting system, the overall collective response of the array can be plotted with respect to the detuning parameter Ω D . With two resonators, there are regions in frequency where three stable solutions can exist. The single mode (SM) of the first intensity corresponds to the forced frequency response of a single Duffing oscillator, and for the second one to a null trivial solution. Remarkably, the double mode (DM) contains two added stable branches as we can see in Figures 3.2 (a) and (b). By zooming over a precise frequency range, we can easily remark that we have a modal interaction and bifurcation topology transfer between these two coupled oscillators. Although, these branches are stable, it is hard to reach them experimentally, because their basins of attraction are very narrow. Consequently we extend the investigations to the case of simultaneous primary and parametric excitations, seeking for additional properties.

Simultaneous primary and parametric resonances

In order to illustrate the complexity of the collective dynamics for the considered periodic structure, the case of simultaneous resonances is numerically investigated for G = 0.1 and H = 7. In addition, with a small change in the excitation amplitudes (G = 0.08 and H = 11), we can obtain up to four solutions for a given frequency for the first intensity response, as shown in Figure 3.4. The frequency response plotted in Figure 3.5, represents the time and space average of the square of the oscillator displacement.

I = 1 N N n=1 < u 2 n >, (3.24)
where the angular brackets denote time average, using the fact that for N = 2,

I = 3 2 (|a 1 | 2 +|a 2 | 2
). The multivaluedness of the response curve due to the nonlinearity has a significance from the physical point of view because it leads to jump phenomena which are localized at the bifurcation points. In addition, the DM has an important amplitude and it is stable over a wide frequency range, which implies significant modal interactions due to nonlinearities.

Basins of attraction

In this subsection, the basins of attraction are used to investigate qualitatively as well as quantitatively the trajectories of the system response, the robustness of the attractors in terms of large basins distribution and their practical implications, for the case of two coupled D-VDP oscillators under simultaneous primary and parametric resonances. The analyzes are performed in a classical way where the robustness is only related to the global size of the attractor. Although the basins of attraction are usually plotted in the phase plane (u n , un ), we choose to represent them in different diagrams for convenience regarding the adopted solving procedure leading to the differential equations (3.17 As a first step, the case Ω D = 0.255 is considered. It corresponds to a multistable response with four attractors for the first intensity and two for the second one. α 1 (0), α 2 (0)) for β 1 (0) = β 2 (0) = 0 and(β 1 (0), β 2 (0)) with α 1 (0) = α 2 (0) = 0. Remarkably, these curves show that the basins of attraction are symmetric with respect to α 2 (0) = 0 and β 2 (0) = 0 and thus, one can investigate their distribution in the Nyquist plane (α 1 (0), β 1 (0)), while setting random positive values of α 2 (0) and β 2 (0). When the first response jumps between SMRB1, SMRB2 and SMNRB, the second one is stabilized on the SM and a similar correspondence exists between the double modes as shown in Figure 3.7. This results in an interesting transfer of basins of attraction topologies between both responses with respect to the type of branches described in Figure 3.4. Thus, one can restrict the analysis to the attractor robustness of the double mode of |a 1 | 2 .

Then, the basins of attraction of the first response are plotted in the Nyquist plane while varying the detuning parameter Ω D to track the evolution of the attractor topology when the oscillators are going from a bistable to a multistable state, as shown in Figure 3.8. For instance, |a 1 | 2 displays two stable resonant solutions for Ω D = 0.07. At the considered detuning parameter, the basins of attraction of SMRB1 are larger than those of SMRB2 due to the difference in the frequency distances separating the corresponding solutions and the bifurcations points B 1 and B 2 . For Ω D = 0.16, another stable solution is added, It is located on the non-resonant branch of the single mode and its basins of attraction take a small domain in the Nyquist plane as shown in The basins of attraction of the double mode take place in Figures 3.8 (c) and (d) and their size increases and becomes significantly large for Ω D = 0.19 with respect to the attractor size of the non-resonant branch which is almost null for Ω D = 0.25. Although these lower branches are very close as shown in Figure 3.4 for |a 1 | 2 , they do not have the same nature and this is proved regarding the frequency response of |a 2 | 2 , which justifies the large differences in their attractor topologies. Figure 3.9 displays the distribution of the basins of attraction of |a 1 | 2 in the Nyquist plane for Ω = 0.255 and three different parametric excitation amplitudes. When H is decreased from 11 to 6.7, the number of stable solutions decreases for the considered detuning parameter at which it is not possible to intercept SMRB2. Indeed, a large part of the basins of attraction of SMRB2 for H = 11 is taken by the double mode for H = 8.5. Moreover, if we decrease the parametric excitation amplitude down to H = 6.7, the oscillator becomes bistable and can switch solution between the double mode and the non-resonant branch of the single mode. Nevertheless, the double mode is more robust, since its basins of attraction represent more than 50% of the whole domain of initial conditions. In practice, the attractor topology can be tuned with respect to the parametric excitation in order to enlarge the basins of attraction of the double mode and obtain an interesting collective dynamics very well adapted for nonlinear energy localization.

Finally, to prove the robustness of the double mode, a quantitative study has been made, based on a random sampling. It consists in solving sequentially the differential equations (3.17) and (3.18) by allowing the initial conditions to be random at each iteration. Of 100 000 sampled initial conditions, the double mode occurs 47%. Interestingly, the basins of attraction of the double mode are large, which demonstrates the attractor robustness. Consequently, once functionalized in term of nonlinearity, the proposed periodic structure can be exploited for jumps-based sensing techniques [Kumar 11,Kumar 12,Nguyen 15]. 3.1. The curves show the squares of the amplitudes of the three different modes, where DM corresponds to the trivial null solution of equation (3.16), for N = 3 and m = 2. Remarkably, all solutions of |a 2 1 | and |a 2 3 | are multimodal. The Double and Triple mode solution branches, result from the coupling between the first and third oscillators, and the three coupled oscillators respectively. They are denoted by: Di j and T k l for i ∈ {1, 3}, j ∈ {1, 2, 3, 4} and k, l ∈ {1, 2, 3}. Thus, in this case the periodic structure is completely driven by the collective dynamics due to the modal interactions between the nonlinear oscillators. This is illustrated in Figure 3.11, where we plotted the average response intensity, defined in equation (3.24), which is In addition, it is noticeable that each bifurcation point due to a multi-modal solution on |a 2 | 2 has a correspondence on the two other intensity responses (for instance, B 1 B 2 and B 3 ), which proves that the bifurcation topology transfer is general for any number of oscillators. Moreover, there are multimode solution branches that are completely disconnected from all other branches. Among them, two Triple Mode solutions localized in the detuning interval [-0.107, 0.007]. These curves were plotted to underline the large number of solution branches even for a small number of coupled oscillators. Indeed, the intensity responses of |a 2 1 | and |a 2 3 | are highly multistable with up to seven stable non-trivial solutions for a given frequency (four double mode and three triple mode solutions). For the configuration, when Ω D = 0.46, the distributions of the basins of attraction for the three intensity responses are plotted in Figures 3.12 (a) and (b) in the Nyquist plane (α 1 (0), β 1 (0)), for random initial conditions α 2 (0) = 0.06, β 2 (0) = 0.48, α 3 (0) = 0.38 and β 3 (0) = 0.39. We remark that |a 1 | 2 and |a 3 | 2 have the same basin distribution and a transfer of basins of attraction topology exists between the three oscillators with respect to the solution branch nature. Furthermore, the probability of reaching the triple mode solution branches is about 32% and their robustness can be adjusted with respect to the design parameters.

Three coupled nonlinear resonators

I = 4 3 (|a 2 1 | + |a 2 2 | + |a 2 3 |) for N = 3.
It is notable in the basins of attraction figures that most fractal parts are negligible with respect to the compact parts of the attractors which implies good agreements with the dynamical integrity of the system [ Soliman 89]. Nevertheless, for large arrays of oscillators, the fractal parts will increase due to the important number of multi-modal solution branches, and therefore, the dynamical integrity [Rega 05, Lenci 08b, Rega 08] must be analyzed by choosing the right definition of safe basin, choosing an appropriate integrity measures to qualify its magnitude and investigating the basins evolution for varying system parameter.

Summary

The collective nonlinear dynamics of periodic nonlinear oscillators was modeled for specific discrete systems of coupled D-VDP oscillators with fifth-order nonlinearity, under simultaneous primary and parametric excitations. A semi-analytical method suitable for nonlinear localized modes problems in perfectly periodic oscillators array has been considered for D-VDP systems. The model is hybrid combining analytical and computational methods and precisely based on the secular perturbation theory with the standing wave decomposition and the asymptotic numerical continuation technique.

The cases of two and three coupled oscillators were investigated in several configurations, where we demonstrate the complexity of the resulting nonlinear frequency curves in terms of modal interactions and bifurcation topology transfer. The complex response curves were shown to emphasize the large number of nontrivial solution branches, even for a small number of coupled oscillators. This can serve as a hint of the important number of multimodal solutions, expected for large number of oscillators. Besides, when the number of coupled oscillators increases, the collective dynamics becomes more complex with a large number of bifurcation points and multimodal interactions over a wide frequency range. Finally, the basins of attraction have been analyzed, precisely in the multistability domain which confirms the robustness of the multi-mode solutions.

In practice, the nonlinearity can be functionalized in such periodic structures in order to generate particular collective dynamics suitable for several applications. Indeed, the model can be used as a design tool in order to increase the number of bifurcations for jump-based multiple mass sensing in micro and nanotechnology or obtaining a large number of stable branches for energy scavenging or trapping applications. Moreover, the stability of these branches can be tuned in the frequency domain for ultra-wide bandwidth filters. For the upcoming chapters, various physical applications can be modeled as coupled nonlinear oscillators array, where they are organized into two chapters, according to their coupling type (linear/nonlinear): 1 and 2D arrays of linearly coupled pendulums under harmonic external excitation in Chapter 4, as well as 1D NEMS array under parametric excitation and 2D granular particles array under harmonic base excitation in Chapter 5.

Introduction

As previously shown in Chapter 3, complex nonlinear phenomena such as the modal interactions, bifurcations and basins of attractions topology transfers have been obtained in weakly coupled nonlinear oscillators. Toward a better understanding of the complex phenomena involved, arrays of axially aligned pendulums with linear torsional coupling are studied as famous example of linearly coupled nonlinear oscillators. The following section is mainly devoted to the study of 1D array of coupled pendulums under harmonic external excitation. The model was extended to investigate the collective dynamics of periodic nonlinear two-dimensional 2D array of coupled pendulums under harmonic base excitation, where the coupled differential equations governing the nonlinear vibrations are solved using an analytical-numerical solving procedure similar to the one employed in 1D arrays.

1D array of coupled pendulums

Introduction

The mechanical pendulum is used to represent classical example to study nonlinear dynamics, which dates back at least five decades. Back to those dates, major contribution were made by Galileo who discovered that the period of swing of a pendulum is independent of its amplitude, Huygens who invented the first pendulum clock was the first to derive the formula for the period of an ideal pendulum and Foucault known for his demonstration of the Foucault pendulum, a device demonstrating the effect of the Earth's rotation. Despite the fact that numerous studies have been made about pendula dynamics, researchers are still active in this research area, including complex dynamics, simple harmonic and subharmonic motions [Skalak 60], coupled pendulums [ Marlin 68], experimental studies [Blackburn 87], chaos [de Paula 06], vibration control [Anh 07] and analytical and numerical studies [Lenci 08a].

Section 1.4.1 in Chapter 1 is devoted to highlight the high number of researches encountered to study the nonlinear dynamics and the complex phenomena in coupled pendulums arrays. Despite the fact that we are interested in studying the nonlinear dynamics in periodic arrays, a system of few dofs, enable analytical solutions for a better understanding of the physical phenomena involved. For instance, Fradkov et al. [Fradkov 07] developed a model for two-pendulum system and obtained from numerical solutions that both in-and anti-phase synchronizations are possible, depending on the initial conditions. The dynamics of the two-pendulum system has been investigated numerically and experimentally, due to the exhibited rich dynamical phenomena [Huynh 13,Witelski 14]. In addition, the synchronization phenomenon of two rotating parametric pendulums attached to common elastic support under harmonic excitation has been studied theoretically and experimentally [ Najdecka 15]. For a three dofs autoparametric system with two connected pendulums, Sado et al. [Sado 10] showed that for this type of system one mode of vibration might excite or damp another one. Although the dynamics of coupled nonlinear pendulums was thoroughly investigated in the frequency and time-space domains, there is a real need to perform profound analysis of the collective dynamics of such systems in order to identify practical relations with the nonlinear energy localization phenomena in terms of modal interactions and bifurcation topologies. .1 shows the theoretical model of an array composed by an horizontal axle, containing N identical pendulums separated by an equal distance. Each pendulum n is attached perpendicularly to the axle, with a mass m and length l, φ n is the rotational dof of the pendulum measured from its stable equilibrium position. The pendulums numbers "n" and "n + 1" are connected together by a linear torsional spring k. Each pendulum is subjected to an external excitation F cos(Ωτ ).

The kinetic energy T of the system can be expressed as

T = 1 2 N n=1 m( ẋ2 n + ẏ2 n ). (4.1)
Where (x n , y n ) represents the Cartesian coordinates system of the n th pendulum and the dot denotes the derivative with respect to time, x n and y n are given as

x n = l sin φ n , y n = l cos φ n . (4.2)
by substituting equation ( 4.2) into equation ( 4.1) the kinetic energy becomes

T = 1 2 N n=1 ml 2 φ2 n . (4.3)
The potential energy V of the system can be written as

V = N n=1 mgl(1 -cos φ n ) + 1 2 N -1 n=1 k(φ n+1 -φ n ) 2 . (4.4)
Where g represents the acceleration due to gravity. The dissipation energy D of the system is given by

D = 1 2 N n=1 αl 2 φ2 n , (4.5)
where α is the viscous damping coefficient. The Lagrange equations can be written in the following form which represents a system of coupled nonlinear equations, subjected to an harmonic external excitation. Figure 4.2 (a) shows the comparison of the response frequency of one pendulum (N=1) using direct time integration method between the full nonlinear problem (4.7) and the one obtained after applying Taylor series (4.11), for the design parameters listen in Table 4.1. equation (4.11) gives a good approximation of the full nonlinear integration for which the error does not exceed 3%. This result is confirmed in Figure 4.2 (b) where we compare the function sin(φ) and its Taylor series as a function of the rotational amplitude φ up to 1.4 rad.

d dτ ( ∂L ∂ φn ) - ∂L ∂φ n = - ∂D ∂ φn + Q n n = 1, . . . , N, (4

Normalized equations

For convenience and equation simplicity, the following nondimensional variables are introduced:

t = ω 0 τ θ n = φ n φ D where φ D = F ω 0 αgl , Q = ω 0 m α and ω 0 = g l (4.8)
Then, we obtain the following nondimensional system of equations

θn + 1 Q θn + 1 φ D sin(φ D θ n ) - k mgl (θ n-1 -2θ n + θ n+1 ) = α ω 0 m cos( Ω ω 0 t) n = 1, 2, . . . ( 4
.9) Equation ( 4.9) represents the normalized nonlinear differential system of a coupled pendulum array, where each pendulum is subjected to an external harmonic excitation. Although this equation can be solved numerically to provide standard reference, the branches resulting from the collective dynamics in the frequency domain can't be identified in terms of modal interactions and bifurcation topologies. To do so, in the following section, we are interested in an analytical solving procedure, enabling the identification of the solution branches type. Therefore, the full model will be replaced by a truncated one valid for small rotational motion.

Solving procedure

As we are interested in small rotational motion, we would expect to get a satisfactory approximation by expanding the nonlinear term sin(φ D θ n ) into its Taylor series up to the third order by writing

sin(φ D θ n ) ≈ φ D θ n - 1 6 φ 3 D θ 3 n . (4.10)
If we replace this approximation in the exact normalized differential equation (4.9), we get the following approximate equation .11) In order to create linear closed modes enabling the study of the collective dynamics, one must choose physical parameters verifying the localization assumption introduced in the previous section. To do so, we suppose that the pendulums are weakly coupled, so that k << mgl (4.12)

θn + 1 Q θn + θ n - k mgl (θ n-1 -2θ n + θ n+1 ) - 1 6 F 2 α 2 gl 3 θ 3 n = α ω 0 m cos( Ω ω 0 t). ( 4
We shall solve equations (4.11) using the multiple-scale perturbation theory, by introducing the appropriate parameters and setting the external frequency an amount εω 0 Ω D away from the resonant frequency, so they can contribute to the amplitudes equations.

1 Q = εc, k mgl = 1 2 εγ, 1 6
F 2 α 2 gl 3 = ξ and α ω 0 m = ε 3 2 f. (4.13)
We express the solution of equations ( 4.9) as a sum of standing-wave modes with slowly varying amplitudes, with fixed boundary conditions (θ 0 = θ N +1 = 0)

θ n (t) = ε 1 2 N m=1 (A m (T ) sin(nq m )e it + c.c.) + ε 3 2 θ (1) n (t) + • • • , n = 1, • • • , N, (4.14)
where T = εt is a slow time variable.Following the same solving procedure employed in the previous chapter for 1D arrays, the complex amplitudes differential system is given by

2i dA m dT + icA m + 2γ sin 2 ( q m 2 )A m - 3 4 ξ j,k,l A j A k A * l ∆ (1) jkl,m = f N + 1 e iΩ D T N n=1
sin(nq m ) (4.15) where the system of complex amplitudes a m is

(ic -2Ω D )a m + 2γ sin 2 ( q m 2 )a m - 3 4 ξ j,k,l a j a k a * l ∆ (1) jkl,m = f N + 1 N n=1 sin(nq m ) (4.16)
The differential system (4.15) of rotational dof φ n (t) has been replaced using perturbation calculations by a time independent mode amplitudes a m system of coupled complex equations. All that remains, in order to study the collective dynamics of an array of coupled pendulums as a function of the original design parameters, is to solve these algebraic coupled complex equations.

Numerical and analytical studies

For the rest of the paper we will consider the design parameters listed in the following Table 4.1, which satisfies the modes localization assumption given in equation (4.12). Note that, since the considered system is periodic and since we expressed the solution as a sum of standing wave modes with slowly varying amplitude (4.14), nearly symmetric responses were obtained. Thereby, for the rest of the paper we choose to plot the intensity responses of the rotational dofs

|φ i | for n ∈ 1, . . . , E( N +1
2 ).

Three coupled pendulums

In order to investigate the modal interactions and their effects on the nonlinear dynamics of a periodic coupled pendulums chain, we consider the case of three coupled pendulums under periodic external harmonic excitation. Before doing so, one must verify qualitatively as well as quantitatively with a brute numerical simulation the validity and the reliability of the proposed analytical-numerical method. 4.1. Comparison between a direct time integration method using Runge-Kutta which is employed to solve the approximate differential system (4.7), and the semi-analytical method.

Concerning numerical simulations, a direct time integration method using Runge-Kutta algorithme has been employed in order to solve the approximate differential system (4.7). However, the resulting algebraic nonlinear system (4.16) of the employed analytical procedure has been solved using a Newton-Raphson algorithm in Mathematica. For several values of the detuning parameter Ω D inside the frequency range, Mathematica represents whole dynamic responses. The stability of the different solution branches has been performed based on the Eigen values of the Jacobian matrix of the differential system (4.15) computed numerically for each point. Figure 4.3 shows the rotational dofs of three coupled pendulums as a function of frequency for the design parameters listed in Table 4.1. The error between both methods rise to 5% on the lower branches. As perturbation methods are well adapted to small damped nonlinear systems; the error decrease by decreasing the damping coefficient.

Direct time integration methods present an easy implementation procedure to plot nonlinear frequency response curves, considering a numerical step-by-step procedure. However, once the multistability domain is reached, jump phenomena may occurs according to the initial conditions that have been taken into account. Consequently, these methods may not give enough information about the global dynamics of the system, particularly the bifurcation behavior where they generally fail to capture unstable solutions and they are not able to identify the nature of multimode branches anyhow. In addition, they are time consuming, since they require highly computational time especially when considering strongly nonlinear high order systems subjected to weak damping.

The employed perturbation technique is based on expanding the periodic solution in the form of a power series in order to obtain an approximate analytical solution of the system. It allows a detailed study regarding the stability, the type of solution branches in term of modal interactions and bifurcation topology transfer. Thus, the proposed solving procedure is robust and efficient to investigate the collective dynamics of weakly coupled nonlinear periodic structures.

Symbolic computations were performed in order to solve the complex algebraic system (4.16) for three coupled pendulums (N = 3) under external excitation, based on prediction-correction algorithms. The natural dimensionless eigenfrequencies of the associated linear system are ω 1 = 1.0017, ω 2 = 1.0059 and ω 3 = 1.01. In addition, the generalized modal forces vector X T F is [0.017, 0, -0.029], which means that the second mode is not excited and does not contribute to the dynamic responses.

In Figures 4.4 (a), we show the solutions for the response intensity of three coupled pendulums under external excitations, as a function of frequency. Branches labeled DM i represent the Double Mode solutions branches involving the excitation of odd modes (a 1 = 0, a 3 = 0 and a 2 = 0), where T M j represent the Triple Mode solutions generated by all modes collectively. For three coupled pendulums all solution branches are multimodal where we can obtain up to six stable solutions for a given frequency. Note that each bifurcation point due to a multimodal solution has a correspondence on the two other intensity responses. These curves were plotted to underline the large number of stable solution branches, even for a small number of coupled pendulums. As the number of multimodal solution branches is important, we analyze each single branch separately to explore the practical appearance of each attractor. Therefore, we make use of the basins of attraction to demonstrate the robustness of these solutions. Figure 4.4 (b) show the basins of attraction in the Nyquist plane (α 1 (0), β 1 (0)) for f = 1.83Hz and a random configuration of initial conditions α 2 (0) = 0.1, β 2 (0) = 0, α 3 (0) = 0.5 and β 3 (0) = -0.3. We can notice that we have basins of attraction topology transfer between odd modes |a 1 | 2 and |a 3 | 2 with respect to the solution branch nature, due to the modal interaction between them. In addition, the distribution of the DM solution branches corresponds to the distribution of the null trivial solution of |a 2 | 2 .

In addition, as the even modes are not excited a Reduced Order Model (ROM) approach considering odd modes only has been proposed in order to reduce the CPU time. Figure 4.5 shows the rotational dofs φ n as a function of frequency comparing results obtained when considering the full model and when applying the ROM approach. For a single pendulum the response intensity corresponds to a forced frequency response of a single Duffing oscillator, however for the three weakly coupled pendulums, we observe additional solution branches resulting from the interaction between localized modes. Starting from the response intensities plotted in Figure 4.4, to the rotational dofs given in Figure 4.5, one can remark that the multimodal solutions were distributed in the multistability domain in a manner to join either resonant or non-resonant branches.

Applying the proposed ROM approach goes to be cancelling the 2nd mode (a 2 = 0), therefore the rotations dofs consist only of DM solutions. The comparison between both approaches shows that the DM solutions capture almost all branches, while the TM (Full model) contributes with additional separated branches for the second rotational dof only. Consequently, as the absence of the 2 nd mode maintains the dominant dynamics without significant loss of accuracy compared to the full model and in order to reduce the CPU time, we propose to exclude it from the projection basis.

Six coupled pendulums

For larger number of coupled pendulums, symbolic computations become timeconsuming, especially when visualizing the whole dynamics of the responses including the unstable branches. We may use an appropriate time integration procedure in order to solve the ordinary differential system (4.15), which is based on the fourth-order Runge Kutta integration method, allowing us to identify the nature of the branches with a low of computational time compared to the case when treating the initial differential system.

For the case of six coupled pendulums (N = 6), the natural dimensionless eigenfrequencies of the associated linear system are: ω 1 = 1.00117, ω 2 = 1.00446, ω 3 = 1.0092, ω 4 = 1.01447, ω 5 = 1.01922 and ω 6 = 1.0225. The generalized modal forces vector X T F is [0.0234, 0, 0.0067, 0, -0.0025, 0], which means that the even modes are not excited after projecting on the standing wave modes. Then, the re-sponses will take the form of a combination of two different types of modes: Triple Mode (T M ) resulting from the interaction between odd excited modes and sextuple mode (6 th M ) solution branches which is driven by the excitation of all modes collectively. As applying the ROM approach consists in writing the odd modes only in the reduction basis, all solution branches in the frequency responses each rotational dof will be of TM type. Therefore, we calculate the square of amplitudes of the odd modes in order to see the contribution of each triple mode solution branch in the rotational dof. Figure 4.6 (b) shows the response intensities of the odd modes as a function of frequency. Dark and light green curves represent the solutions that contribute as a Resonant (R) and Non-Resonant (NR) branches in the rotational dofs respectively. Remarkably, compared to the case of three coupled pendulums, the number of multi-modal solutions increases with up to eight possible solutions for a given frequency, for the case of six coupled pendulums even after applying the ROM approach. In addition, we can see that the multimodal solutions are distributed half-wave between R and NR in the rotational dofs.

N coupled pendulums

For N weakly coupled pendulums, where each of these pendulums is subjected to an external harmonic excitation, the dimensionless eigenfrequencies of the associated linear system can be expressed as: ω m = 1 + ε m and the odd modes are not excited and do not contribute to the dynamic responses. The multi-mode solution branches of the corresponding response intensities are generated either by the excitation of odd modes or by all modes collectively. Then, a ROM approach based on using the even modes only can be applied, reducing the CPU time, while preserving the accuracy of the full model. In addition, the symmetry of the pendulums array leads to identical dynamics behavior of pairs of dofs (φ i , φ N +1-i ). As the number of coupled pendulums increases, the number of multi-modal solution branches increases and the system becomes more complex.

Figure 4.7 shows the rotational dofs of twelve coupled pendulums as a function of frequency, where the ROM approach was applied by only writing the odd modes in the reduction basis. The symmetry of the array of pendulums leads to identical dynamics behavior of six pairs of dof (φ 1 , φ 12 ), (φ 2 , φ 11 ), (φ 3 , φ 10 ), (φ 4 , φ 9 ), (φ 5 , φ 8 ) and (φ 6 , φ 7 ). The response frequencies of N coupled pendulums have similar-trending curves, where the multimodal solutions join either R or NR branches in the rotational dofs. In order to highlight the benefits of treating large number of degrees-of-freedom and extract some interesting features, a study of basins of attraction in the phase portrait can be performed for a given frequency in the multistability domain. Color 4

Basins of attraction analysis

Mode DM DM TM TM QM QM 5M 5M 6M 6M Branch R NR R NR R NR R NR R NR Figure 4
.9 shows the evolution of the basins of attractions in terms of type of branch and nature of mode of |φ 1 | for f = 1.83Hz in the phase portrait (φ 1 (0), φ1 (0)) for three, four, five and six coupled pendulums simultaneously. Table 4.2 shows the color palette, where each color is related to a specific type of branche. DM , T M , QM , 5M and 6M are respectively the Double, Triple, Quadruple, Quintuple and Sextuple Modes.

From the previous diagrams, the basins of attractions of the DM are larger for N = 3 than for N = 4. Moreover the distribution of the TM in the case of three coupled pendulums (N = 3) covers small areas compared to the case of N = 5; this distribution decreases for six coupled pendulums. Remarkably, the distribution of multimodal solution branches generated by all modes simultaneously increases while increasing the number of coupled pendulums. This can serve as a hint of the important distribution of resonant multimodal solutions, expected for large number of coupled pendulums.

2D pendulums array

Introduction

Taking into account the considerable work in the literature devoted to analyze the richness of the nonlinear dynamical effects found in an array of coupled oscillators, 2D systems are expected to present additional dynamics phenomena not realizable in the 1D ones. Recently, the first analysis of a 2D array of coupled pendulums proposed a computational model to study the nonlinear dynamics of such system under parametric excitation [Jallouli 15]. As seen in the previous section, modes localization in periodic nonlinear pendulums array exhibit a rich and complex collective dynamics. Hence, the primary motivation to study the collective dynamics in 2D periodic lattices is to explore the modes localization phenomenon taking into account different coupling linear and nonlinear stiffness parameters in the x and y directions. For a qualitative understanding of their effects on the modal interactions, the bifurcation topology, the responses bandwidth and the vibrational amplitude. The present section investigate the collective dynamics of periodic nonlinear 2D monoatomic lattices under external excitations which are modeled as spring-mass lattices. The coupled differential EOMs are solve using an analytical-numerical approach similar to the one employed in the previous section. Taking into account the considerable work in the literature devoted to analyze the richness of the nonlinear dynamical effects found in an array of coupled pendulums [Vyas 01], 2D systems are expected to present additional dynamics phenomena not realizable in the 1D case. Recently, the first analysis of a 2D array of coupled pendulums proposed a computational model to study the nonlinear dynamics of such system under parametric excitation [Jallouli 15].

Problem Formulation

Figure 4.10 shows the discrete model of 2D coupled N × M -pendulums subjected to harmonic horizontal base excitation. The system is composed of N equidistant axles A n . Each axle has M equidistant pendulum "i, j" attached perpendicularly to the axle and designed as a uniform rod of length l and mass m fixed at the end. φ i,j is its rotational dof from its stable equilibrium position. Pendula "i, j"and "i + 1, j", "i, j" and "i, j+1" are connected respectively by torsional and linear springs k x and k y . The periodic system is subjected to harmonic horizontal base excitation y = y 0 cos(Ωτ ), with fixed boundary conditions φ i,M +1 = φ i,0 = φ N +1,j = φ 0,j = 0. The kinetic energy T of the system can be expressed as

T = 1 2 N i=1 M j=1 m( ẏ2 i,j + ż2 i,j ) (4.17)
where (y i,j , z i,j ) represents the Cartesian coordinates of the pendulum {i, j}, given as

y i,j = y + l sin φ i,j , z i,j = l cos φ i,j , (4.18)
and the dot denotes typically the derivative with respect to time. Substituting Eq. ( 4.18) into Eq. (4.17) the kinetic energy becomes

T = 1 2 N i=1 M j=1 m l 2 φ2 i,j + 2l cos(φ i,j ) ẏ φi,j + ẏ2 . (4.19)
The potential energy of such a problem can be written as:

V = N -1 i=1 M -1 j=1 1 2 k x (φ i+1,j -φ i,j ) 2 + (φ i-1,j -φ i,j ) 2 + 1 2 k y ∆ 2 j,j+1 + ∆ 2 j-1,j + mg l -l cos(φ i,j ) , (4.20)
where g represents the acceleration of gravity and ∆ j,j+1 the distance between pendulums {i, j} and {i, j + 1}, which can be written as

∆ j,j+1 = l(sin(φ i,j+1 ) -sin(φ i,j )).
(4.21)

The Lagrange equations of motion can be written in the following form

d dτ ( ∂L ∂ φi,j ) - ∂L ∂φ i,j = Q i,j for i = 1, . . . , N and j = 1, . . . , M (4.22)
with L the Lagrangian operator defined by L = T -V , Q i,j = αl 2 φi,j represents the generalized forces applied to the pendulum {i, j} and α the viscous damping coefficient. Substituting equations (4.19) and ( 4.20) into the Lagrange equation (4.22), we obtain the equation of motion of each pendulum {i, j} ml 2 φi,j + αl 2 φi,j + mgl sin(φ i,j ) + k x (2φ i,j -φ i-1,j -φ i+1,j ) +k y l 2 cos(φ i,j )(2 sin(φ i,j ) -sin(φ i,j-1 ) -sin(φ i,j+1 )) = mly 0 Ω 2 cos(φ i,j ) cos(Ωτ ) (4.23) By expanding sin(φ i,j ) and cos(φ i,j ) into Taylor series up to the third and first order respectively, the equations (4.23) can be written as:

ml 2 φi,j + αl 2 φi,j + k x (2φ i,j -φ i-1,j -φ i+1,j ) + mgl(φ i,j -1 6 φ 3 i,j ) +k y l 2 (2φ i,j -φ i,j-1 -φ i,j+1 ) -kyl 2 6 (2φ 3 i,j -φ 3 i,j-1 -φ 3 i,j+1 ) = mly 0 Ω 2 cos(Ωτ ) (4.24)
Since we are interested in weakly coupled parameters in both x and y directions, the third order expansion related to the linear coupling parameter in the y direction is considered to be negligible.

Dimensionless equations

For convenience and equation simplicity, the following nondimensional variables are introduced:

θ i,j = φ i,j φ D , t = ω 0 τ, (4.25)
where φ D = my 0 Ω 2 αlω 0 is the dynamic displacement of the associated linear system, Q = ω 0 m α and ω 0 = g l represents the natural frequency of pendulum "i, j". After dropping Eqs. (4.25) into Eq. ( 4.24) and dividing by α m 2 ly 0 ω 0 Ω 2 , we obtain the following nondimensional system of equations

θi,j + 1 Q θi,j + θ i,j + kx m g l (2θ i,j -θ i-1,j -θ i+1,j ) + ky l m g (2θ i,j -θ i,j-1 -θ i,j+1 ) -1 6 ( m y 0 Ω 2 α l ω 0 ) 2 θ 3 i,j = α m ω 0 cos( Ω ω 0 t) for i = 1, . . . , N and j = 1, . . . , M (4.26)
Equation ( 4.26) describes a 2D monoatomoc lattice consisting of an arrangement of identical unit cells assembled together to form the whole structure, where each unit cell adjoins four neighbors via linear springs kx mgl and kyl mg in the x and y directions respectively and subjected to harmonic horizontal base excitation. For boundary condition, we suppose that we have fixed identical systes in both directions, where we define 2(N + M + 2) extra variables and set them to zero as θ 0,j = θ N +1,j = θ i,0 = θ 0,M +1 = 0.

Linear study

To determine the natural frequencies and their eigenvectors, equation (4.26) can be written in matrix form as

Ü + C U + K L U Linear part +K N L (Θ) Nonlinear part = F (t), (4.27)
where

Θ = [θ 1,1 , θ 1,2 , . . . , θ 1,M , . . . , θ i,1 , θ i,2 , . . . , θ i,M , . . . ] is the N×M-dimensional variable vector, C the damping matrix, F (t) = α mω 0 cos[ Ω ω 0 t][1, . . . , 1
] T the excitation vector, K N L (Θ) the nonlinear stiffness vector and

K L =         K 1 -K 2 -K 2 K 1 -K 2 0 . . . . . . . . . 0 -K 2 K 1 -K 2 -K 2 K 1         (4.28)
where

K 1 =           1 + 2kx mgl + 2kyl mg - kyl mg - kyl mg 1 + 2kx mgl + 2kyl mg - kyl mg 0 . . . . . . . . . 0 - kyl mg 1 + 2kx mgl + 2kyl mg - kyl mg - kyl mg 1 + 2kx mgl + 2kyl mg           (4.29)
and K 2 = kx mgl I M with I M the identity matrix. The eigenvalues ω i,j are represented in the following form for few coupled oscillators in contact:

N = 1, M = 1 : ω 1,1 = 1 + 2k x mgl + 2k y l mg N = 2, M = 1 : ω 1,1 = 1 + k x mgl + 2k y l mg ω 2,1 = 1 + 3k x mgl + 2k y l mg (4.30)
We may express all normal frequencies relative to the same non-dimensional reference frequency which is 1, so that

ω i,j = (1 + λ i k x mgl + γ j k y l mg ) (i = 1, . . . , N and j = 1, . . . , M ), (4.31)
where λ i and γ j are the coefficients depending on the mode of vibration (As an example for N = 1 and M = 1, λ 1 = 2 and γ 2 = 2). For weak linear coupling, we suppose that k x << mgl and k y << mg l , consequently, kx mgl << 1 and kyl mg << 1 and:

ω i,j ≈ 1 + 1 2 λ i k x mgl + 1 2 γ j k y l mg (i = 1, . . . , N and j = 1, . . . , M ). (4.32)
This assumption leads to the creation of linear closed modes which permits to the influence of the modes localization on the collective dynamics of a 2D pendulums array. In the following section, we will proceed to solve the normalized equivalent differential system (4.26) using a perturbation technique based on the multiple time scales combined with the standing waves decomposition as presented in Chapter 2, section 2.8.4.

Semi-analytical approach

The multiple scales method coupled with standing wave decomposition is employed as an approximate analytical solution to solve the coupled nonlinear differential systems. The main advantages of the present approach are the capacity of handling weakly coupled nonlinear systems, which permits to visualize all physical responses branches and their properties in terms of modal interactions and bifurcation topology transfer. Therefore, we rescale the variables appearing in the equation (4.26) as follows and set the external frequency an amount εω 0 Ω D away from the resonant frequency

1 Q = εη, k x mgl = 1 2 εγ x , k y l mg = 1 2 εγ y , - 1 6 ( m y 0 ω 2 α l ω 0 ) 2 = ξ, α m ω 0 = ε 3 2 f. (4.33)
In order to calculate the responses of each pendulum, we expand θ i,j (t) as a sum of standing-wave modes with slowly varying amplitude, suitable for weakly coupled 2D nonlinear lattices as described in detail in the previous chapter.

θ i,j (t) = ε 1/2 N r=1 N p=1 (A r,p (T ) sin(nq r ) sin(md p )e it + c.c.) +ε 3/2 θ (1) i,j (t) + • • • , i = 1, ..., N and j = 1, ..., M (4.34)
where T = εt is a slow time variable, that authorizes the complex amplitude A r,p (T ) to vary slowly in time. Since we proposed fixed boundary conditions, the possible wave components q r , d p can be given as

q r = rπ N +1 , r = 1, ..., N d p = pπ M +1 , p = 1, ..., M (4.35)
After replacing the proposed solution (3.11) into the equation of motion, we can get at the order of ε 3 2 , N × M equations of the form:

θ(1) i,j + θ (1)
i,j = r p (r, p th secular term)e it + other terms (4.36)

We must eliminate the N × M secular terms so that θ i,j (t) remains finite which allows us to determine the equations for the slowly varying amplitudes A r,p (T ). To extract the equations for these amplitudes, we make use of the orthogonality of the modes, by multiplying the r, p th secular term by sin(nq r ) sin(dq p ) and summing over n and m. We find that the equation of the r, p th amplitude A r,p (T ) is given by: 2i dAr,p dT + iηA r,p + 3ξ 16

i,k,l j,s,o

A i,j A k,s A * l,o ∆ (1) ikl;r ∆ (1) jso;p +2(γ x sin 2 ( qr 2 ) + γ y sin 2 ( dp 2 ))A r,p = 2f (N +1)(M +1) n,m
sin (nq r ) sin (md p )e iΩ D T (4.37)

Ignoring initial transients, and assuming that the nonlinear terms in the equation are sufficient to saturate the growth of the instability, we try a steady-state solution of the form A j,p = a j,p e iΩ D T (4.38)

Substituting equation ( 4.38) into the differential system (4.37) of amplitude, we obtain the required equation for the fixed complex amplitudes a j,p .

2iΩ D a r,p + iηa r,p + 2(γ x sin( qr 2 ) 2 + γ y sin( dp 2 ) 2 )a r,p + 3ξ 16 i,k,l j,s,o a i,j a k,s a * l,o ∆ (1)
ikl;r ∆

(1)

jso;p = 2f (N + 1)(M + 1) n,m sin (nq r ) sin (md p )(4.39)
The complexity of the nonlinear algebraic system (4.39), limits the possibility to solve it analytically even for few coupled pendulums. We may use the appropriate time integration procedure in order to solve the differential system (4.37) numerically, which implements a version of the fourth-order Runge Kutta integration method. In order to study the collective dynamics of the 2D periodic structure, two different configurations were proposed in Table 4.3, which highlights the influence of linear coupling on the nonlinear behavior in terms of frequency range and vibration amplitude.

Numerical simulations

For the case of 2 × 2 coupled pendulums (N = 2, M = 2) and for the first design parameters listed in Table 4.3, the natural dimensionless eigenfrequencies of the associated linear system are very close to 1. In addition, as we considered k y = kx l 2 , the linear springs are expected to act the same way. Figures 4.11 show the response intensities as a function of the detuning parameter Ω D . The Single Mode (SM) solution branches corresponds to the trivial solution of the system (4.37) with A 1,2 = A 2,1 = A 2,2 = 0. The multimodal solutions DM and QM are generated by the coupling between two and all modes respectively. These curves were ploted to underline: first the effect of modes localization on the high number of solution branches for a given frequency and second the identical contribution of the non-excited modes on the responses. Although, the advantages of modes localization, considering small coupling parameters does not affect the system's rotational dof. Figure 4.12, displays the rotational dof of 2×2 coupled pendulums under harmonic horizontal base excitation as a function of frequency, which corresponds to a single forced Duffing oscillator. In addition, all multimodal solution branches are distributed between resonant and non-resonant branch, hence it's complicated to identify their contribution. We decide to increase linear coupling parameters in such a manner to preserve the localization phenomena. Therefore, we consider the second design parameters listed in Table 4.3, so that the natural dimensionless eigenfrequencies of the associated linear system are ω 1,1 = 1.0331, ω 1,2 = 1.0336, ω 2,1 = 1.0989 andω 2,2 = 1 A 1,1 andA 2,1 or A 2,2 present larger frequency range and their bifurcation point B 2 exceeds the one that corresponds to the non-resonant branch of the SM (B 1 ). In addition, in regards to the system's rotational dof represented in figure 4.14, multimodal solution branches present additional features distinct from those presented by a single forced Duffing oscillator (larger bandwidth and higher vibration amplitude).

Summary

Pendulums arrays were modeled as weak linearly coupled nonlinear oscillators in order to study their collective dynamics. The choice of the physical parameters ensuring sufficiently weak substructure coupling, leads to strongly modes localization for perfectly periodic arrays. A semi-analytical method suitable for nonlinear localized modes problems has been considered, based on a perturbation technique, combined with a standing wave decomposition, transforming the normalized differential systems into a set of coupled complex algebraic ones. In addition, the validity of the proposed semi-analytical method has been verified, and its role in identifying the type of the solution branches has been highlighted. The resulting systems have been numerically solved and plotted as a function of frequency to highlight the complexity and the multivaludness of the responses. High number of multimodal solutions has been obtained for few number of coupled pendulums, resulting from the excitation of several modes collectively.

For 1D array of coupled pendulums under harmonic external excitation, the robustness of the additional multimodal solution branches has been proved with a classical study of their basins of attraction, as well as the bifurcation topology transfer. Identical dynamical behavior between frequency responses, resulting from the symmetry of the array has been identified. In addition, we notice that the multimodal solutions are distributed in a manner to join either resonant or non-resonant branches. A model order reduction approach ROM has been applied to solve large periodic arrays of coupled oscillators under external forces which helps reducing the CPU time. It is based on projecting on odd modes, while maintaining the dominant dynamics of the responses without significant loss of accuracy compared to the results of the full model.

In addition, the model has been extended to study the collective nonlinear dynamics of a 2D periodic array of coupled pendulums under harmonic horizontal base excitation. For few coupled pendulums, the resulting complex system has been numerically solved for two different configurations, highlighting the impact of linear coupling parameters in both directions on the complexity and the multivaludness of the responses.

Remarkably, the distributions of the basins of attraction of the resonant branches and the multimodal solutions governed by all modes simultaneously increase while increasing the number of coupled pendulums. In practice, this model could be used as a design tool to tune the number of multimodal solutions and increase the possibility to reach resonant branches. two physical applications. The first one represents an array of electrostatically coupled nanobeams subjected to parametric excitation. The nonlinear coupling terms arise from the electrostatic and dissipative forces acting between NEMS resonators. The second example is a 2D periodic granular particles array subjected to compressive loadings and a harmonic base excitation, where each particle is coupled with fixed-fixed uniform beam and the nonlinear dynamics is governed by the Hertzian contact.

NEMS array

Introduction

External potentials and geometric effects in MEMS resonators possess strong and easily reachable nonlinearities which requires further investigations of the coupling, stability, nonlinearity and reliability [Zhang 07b]. For instance, nonlinear dynamic characteristics of microbeam-based resonant sensors in MEMS were investigated numerically and experimentally [Zhang 07a, Mestrom 08, Zhang 05, Ghayesh 13]. Recently, Tajaddodianfar et al. [Tajaddodianfar 16] used the homotopy analysis method (HAM) in order to derive approximate analytical solutions for the nonlinear frequency response of MEMS/NEMS resonators characterized by the general Duffing equation having quadratic and cubic nonlinearities. Besides, electrostatically actuated microresonators are adequate for parametric excitation that attracted a considerable attention in the past few years because of its interesting properties of signal amplification [Carr 00] and noise squeezing [ Rugar 91]. For instance, Zhang et al.

[ Zhang 02] presented the effect of the cubic mechanical stiffness and electrostatic stiffness in changing the dynamic behavior of the oscillator response when excited parametrically while Rhoads et al. [Rhoads 06b,Rhoads 06a] investigated the nonlinear response of resonant microbeam systems with purely-parametric electrostatic actuation.

A significant part of Roukes et al. work [Kozinsky 06,Villanueva 13,Matheny 13] involved deeper and detailed understanding of NEMS nonlinearities origins and the resulting limitation for sensing applications. Knowing that, the smaller the structures, the sooner nonlinearities occur which in turns must be taken into account in terms of high order nonlinearities [Kacem 11c, Kacem 10], Kacem et al. [Kacem 09b,Kacem 11b] presented and validated a compact and analytical model of an electrostatically driven clamped-clamped beam, under primary resonance and accounting the principal nonlinearities with cubic-quintic terms. Particularly, they demonstrated the physical meaning of the fifth order nonlinearity, based on the experimental observation of mixed hardening/softening effects in MEMS [Kacem 09a]. Moreover, they studied experimentally the sensitivity of hysteresis suppression to high order nonlinearities in an externally driven MEMS beam at large amplitudes [Kacem 11d]. They also investigated analytically and experimentally the benefits of a superharmonic excitation on the primary resonance of a microresonator by sub-jecting the nonlinear oscillator to simultaneous slow and fast excitations, where they showed how the bifurcation topology of an undesirable unstable behavior is modified when the resonator is simultaneously actuated at its primary and superharmonic resonances. Recently, hysteresis suppression was demonstrated in the case of a nanocantilever with electrostatic actuation and piezoresistive sensing while the sensor dynamics is stabilized thanks to a nonlinear resonance [ Kacem 15].

Section 1.4.2 in chapter Chapter 1 is devoted to study the nonlinear dynamics of coupled MEMS/NEMS arrays. In this study, we model the collective dynamics of an electrostatically actuated, parametrically driven nonbeams array, while including the main sources of nonlinearities up to the fifth order.

Equations of motion

We consider an array of N doubly clamped silicon nanobeams with identical material properties, where the excitation and damping forces are located at the midspan of the resonators, polarized by a DC voltage V dc and subject to an alternative AC voltage V ac at the frequency Ω. b denotes the nanoresonator width, h is its thickness, g is the separation distance between the beams width and L represents the length of the beam (see Fig. 5.1). We may use a system of coupled ODEs that captures the important physical aspects that are required to model one dimensional array of parametrically driven coupled nanomechanical beams. A typical coupled resonators array is described after an appropriate scaling by a dimensionless coupled equations of motion of the form:

ÿn + εµ ẏn + y n + εξy 3 n + ε[Γ + H cos(2(1 + εΩ D )t)](y n-1 -2y n + y n+1 ) +εψ[(y n -y n-1 ) 3 + (y n -y n+1 ) 3 ] + εδ[(y n -y n-1 ) 5 + (y n -y n+1 ) 5 ]
+ελ[(y n -y n-1 ) 2 ( ẏn -ẏn-1 ) + (y n -y n+1 ) 2 ( ẏn -ẏn+1 )] = 0 where y n describes the single degree of freedom of the nth beam, giving it displacement from equilibrium, with n = 1 . . . N , as we consider two additional fixed beams at the ends of the array y 0 = y n+1 = 0. There are three different forces acting on the resonators: the mechanical, the electrostatic and dissipative forces. ξ is the nonlinear mechanical term resulting from midplane stretching and µ is the linear damping coefficient acting on the n th respective nanobeam. The electrostatic coupling between nearest-neighbor beams generates several nonlinearities, as the parametric excitation where H defines its amplitude and Ω D the detuning parameter, the electrostatic cubic and quintic nonlinearities are denoted by ψ and δ respectively. λ denotes the VDP coupling generated by the dissipative coupling forces.

As typical NEMS resonators has high quality factors, we define a small expansion parameter εµ where ε << 1 and µ of order unity. We express the driving amplitude as εH, with H of order unity to be at the same order as the linear damping. To ensure that the nonlinear terms have the ability to cause weak oscillations, we write them with corresponding scaling to get them in the system at the same range as the physical impacts.

Complex amplitude equations

To calculate the nonlinear response of each parametrically excited beam, we use the method of multiple scales, where we express the solution as a sum of standing wave modes with slowly changing magnitude, by writing the displacement of the n th nanobeam far from its equilibrium of the form:

y n (t) = N m=1
(A m (T ) sin(q m )e it + c.c.) + εy (1) n (t) + • • • n = 1, ..., N (5.1) By substituting solution (5.1) term by term into the differential EoM we obtain up to the order ε a system that after applying the solvability condition and vanishing all secular terms, we get the following differential system for the slowly changing magnitudes A m (T )

-2i dAm dT -iηA m + 2H sin 2 ( qm 2 )A * m e 2iΩ D T -3ξ 4 j,k,l A j A k A * l ∆ (1) jkl,m + 4Γ sin 2 ( q m 2 )A m -4(iλ + 3ψ) sin( qm 2 ) j,k,l S j S k S * l ∆ (2)
jklop,m -40δ sin(

q m 2 ) j,k,l,o,p S j S k S l S * o S * p ∆ (3)
jklop,m = 0

(5.2)

with S v = sin( qv 2 )A v . Trying the steady-state solution of the form A m (T ) = a m e iΩ D T which oscillate at half the parametric resonance into the equations (5.2) of time varying amplitudes, we obtain a system of complex time-independent amplitudes

NEMS array

[2Ω D + 4Γ sin 2 ( q m 2 ) -iη]a m - 3ξ 4 j,k,l a j a k a * l ∆ (1) jkl,m + 2H sin 2 ( q m 2 )a * m -4(iλ + 3ψ) sin( q m 2 ) j,k,l T j T k T * l ∆ (2)
jklop,m -40δ sin(

q m 2 ) j,k,l,o,p T j T k T l T * o T * p ∆ (3)
jklop,m = 0

(5.3) with T v = sin( qv 2 )a v .
In order to obtain the overall collective response of an array of electrostatically coupled nanobeams, we proceed in the coming paragraph to investigate the complex algebraic system (5.3). These equations are complex with a high number of parameters and nonlinearities up to the fifth order, which are difficult or even impossible to solve analytically. In order to study the nonlinear dynamics of coupled oscillators, we start by solving the complex algebraic equation analytically. By expressing each amplitude by a m = |a m |e iΦm , taking the magnitude squared of both separated sides and the intensity |a m | 2 of the nontrivial responses as all positive roots, the SM solution branches satisfies the following equation: Beside the SM solutions, additional ones may exist but are hard to find analytically as the complex algebraic system (5.3) has a large number of variables with high nonlinearity up to the fifth order. Consequently, ManLab [ MAN 10] has been employed for iterative continuation and bifurcation investigations of the nonlinear system as described in detail in Chapter 3.

Nonlinear responses to parametric excitation

1 4 sin 4 ( qm 2 ) 3 4 ξ∆ (1) mmm,m + 12ψ sin( q m 2 )∆ (2) mmm,m |a m | 2 + 40δ sin( q m 2 )∆ (3) mmmmm,m |a m | 4 -4Γ sin 2 ( qm 2 ) -2Ω D 2 + η + 4λ sin( qm 2 )∆ (2) mmm,m |a m | 2
We start with the dynamical analysis of two electrostatically coupled nanomechanical beams under parametric excitation. Figures 5.2 (a) (b) show the steady state responses when considering design 1 described in Table 5.1. The dynamic responses were plotted for a gap g = 215 nm while the nanobeam thickness is h = 150 nm. In this case, the beam thickness and the gap width are almost equal, thus the electrostatic coupling parameters are of the same importance as the mechanical Duffing term. (5.4). Note that as |A 1 | 2 is excited after modes projection, its single mode takes the shape of a single Duffing oscillator with cubic and quintic nonlinearities under parametric excitation. The competition between hardening and softening behaviors leads to a mixed nonlinear behavior and the possibility to get range of frequency where four steady-state responses can be obtained. In addition, the excitation of both modes collectively driven by the coupling between the nanobeams and the localization of the normal modes gave birth to three additional stable branches (Labeled Double Mode DM). In Figure 5.2 (c) we zoom to highlight the multimodal interactions and bifurcation topology transfer between both mode amplitudes. We define the time and space average of the square of the beam motion as .5) Chevrons <> denote the average over time, which is

I = 1 N n=1 N < y 2 n > . ( 5
I = 3 2 (|a 1 | 2 + |a 2 | 2
) for two-beams system. Figure 5.2 (d) shows the square of modes amplitudes of twonanobeams system under parametric excitation. Hysteresis effects are observed with high number of steady-state responses for a given frequency with few coupled beams. Contrary to what have been treated before, we choose to decrease the electrostatic coupling parameters by considering a gap width which is greater than the beam thickness as given in for the design 2 in Table 5.1. (5.4), involving the excitation of one mode only when all the others are null. We should underline the fact that, considering a gap greater than the nanobeam thickness, the coupling parameters (Linear, Dissipative and electrostatic third and fifth Duffing terms) are negligible compared to the Duffing parameter generated by the mechanical nonlinearity, which implies a hardening behavior while maintaining the elliptical shape of a single parametrically excited Duffing oscillator.

Moreover, additional stable DM solution branches are obtained. We can distinguish between our study from Lifshitz et al. [Lifshitz 03] work, with the following differences: the consideration of the electrostatic cubic and quintic duffing coupling terms and the linear damping coefficients that are acting of each beam.

With two beams we can obtain regions in frequency where three stable solutions can exist. For instance, for Ω D = 0.06, we plot the basins of attractions in the Nyquist plane (Re(a 1 (0)), Im(a 1 (0))) with Re(a 2 (0)) = 0.2 and Im(a 2 (0)) = -0.2 for the first and the second amplitudes in addition to the distribution of modes type. Figures 5.4 (a) (b), show modal interactions and bifurcation topology transfer between both amplitudes: when the first amplitude is null the second one reaches a SM and vice versa; then, both amplitudes reach a DM branch together. In addition, the distribution of the DM is large and robust so its corresponding branch can be reached experimentally for a specific choice of initial conditions (See yellow distribution in Figure 5.4 (c)).

In term of average, Figure 5.4(d) shows the response frequency for two-nanobeams system under parametric excitation for the 2 nd set of model parameters. The multivaludness of the response curve leads to complex jump phenomena which is located at the bifurcations despite the small number of coupled oscillators. Figure 5.5 shows the frequency curves of three electrostatically coupled nanobeams under parametric excitation, for the 2 th model parameters. We should note that, for three coupled oscillators the only elliptical SM solution branch is related to the second mode (Blue curve), whose exact analytical solution is given by equation (5.4) for m = 2. Red curves refer to three additional DM solution branches which are stable over a wide frequency range. These branches occur when ∆ (1)

mmm,i = ∆ (2) mmm,i = ∆ (3) mmmmm,i = ∆ (1) mmm,j = ∆ (2) mmm,j = ∆ (3)
mmmmm,j = 0, for all i, j = m and i = j. In addition, we find three Triple Mode (TM) responses (Green curves), resulting from the excitation of all three modes collectively. We highlight the high number of nontrivial solution branches obtained with few coupled oscillators, where we can get five stable responses in a specific frequency range. For high number of coupled beams we expect to find enormous number of complex nonlinear solutions.

The following investigations considered the basins of attractions as a qualitative tool for detecting the trajectories of the attractors, their type and their practical involvement for three parametrically driven, electrostatically coupled nanobeams. The analysis were performed in the Nyquist plane (Re(a 1 (0)), Im(a 1 (0))) for a detuning parameter Ω D = 0.16 in the multistability domain with five attractors for the first and third amplitudes, and four for the second one, while setting Re(a 2 (0)) = -0.1, Re(a 3 (0)) = -0.05, Im(a 2 (0)) = Im(a 3 (0)) = 0. Figure 5.7 shows the evolution of the basins of attraction in term of response intensity for a given frequency in the multistability domain of the first amplitude while decreasing the length of the beams. The basins are plotted in the Nyquist plane (Re(a 1 (0)), Im(a 1 (0))), for the case of three coupled beams, the second set of design parameters listed in Table 5.1, Re(a 2 (0)) = -0.1, Re(a 3 ) = -0.05 and Im(a 2 (0)) = Im(a 3 (0)) = 0. For Ω D = 0.18, |a 1 | 2 can reach up to 5 stable solutions according to the chosen initial conditions. In Figure 5.7 (a) we can clearly see that all multimodal solutions are accessible for the chosen initial conditions. Decreasing the length of the beam from 20 µm to 18.5 µm while conserving the same initial conditions, affects the electrostatic coupling terms while conserving the mechanical nonlinearities. We should note that, the red distribution which represents the upper TM solution for the first amplitude disappears for a slight decrease in the beam length, and has been replaced by the distribution of one of the DM solution (green part). Circumstantially, we prove that we should consider the nonlinear electrostatic coupling terms even if they are small, as they play an important role in having complex phenomena.

2D granular particles lattice

Introduction

The study of granular particles as an application of a 2D periodic nonlinear lattice was motivated by several reasons. For instance, we can mention the simplicity of their constructions (multidimensional systems) and their applicability to engineering (passive vibration damping for example), the ability to tune their dynamic response from linear to weakly and strongly nonlinear by changing the amount of static compressive load. In addition, their great controllability of the assembly, where the periodicity leads to dispersion which ends to the creation of solitary waves by the presence of nonlinearity of defects. In Chapter 1 we reported the state of the art of the studies related to multidimensional granular crystal lattices. In this section, we are interested in a novel 2D periodic structure of compressed granular particle and subjected to horizontal base excitation. Before doing so, we will introduce the contact between spherical particles, then present the origin of the nonlinearities in 1D array to end up studying the collective nonlinear dynamics in the 2D lattice. Considering two different spheres 1 and 2 of radii R 1 and R 2 compressed with static forces F, as showed in Fig. 5.8. E 1 and E 2 represent the young's modulus for spheres 1 and 2 and v 1 and v 2 are the Poisson's ratios, respectively. The maximum contact pressure at the center of the circular contact area is given by [Hertz 99]

Contact of spheres

P max = 3F 2πa 2 . (5.6)
a is the radius of the contact area given by:

a = 3F R 4E 1/3
(5.7)

where E represents the contact modulus defined by

1 E = 1 -v 2 1 E 1 + 1 -v 2 2 E 2 (5.8)
And R the effective radius given by

1 R = 1 R 1 + 1 R 2 .
(5.9)

The current approach of both spheres centers δ 0 (The overlap deformation) is related to the maximum contact pressure by

δ 0 = a 2 R = 3 9F 2 16RE 2 .
(5.10) Considering the case of a statically compressed 1D monoatomic granular crystal composed of N identical elastic spherical particles as shown in Figure 5.9. Each sphere has R as radius with E and ν the Young's modulus and the Poisson's ratio respectively. The contact elastic force F i,i+1 between two adjacent particles is given by

1D periodic array of granular particles

F i,i+1 = A i,i+1 [∆ i,i+1 ] n i,i+1 +
( 5.11) where n i,i+1 = 3 2 and A i,i+1 is defined as

A i,i+1 = A = E √ 2R 3(1 -ν 2 )
.

(5.12)

[∆ i,i+1 ] 3 2
+ represents the relative displacement of their centers. [∆ i,i+1 ] + which takes the value of ∆ i,i+1 if ∆ i,i+1 > 0 and 0 when ∆ i,i+1 ≤ 0; which signifies that there is no force between separated particles. Assuming that the chain is subjected to constant static force F 0 applied to both ends and neglecting any dissipation, the dynamics displacement u i of the i th particle can be described by the following differential system

mü i = A(δ 0 + u i-1 -u i ) 3/2 -A(δ 0 + u i -u i+1 ) 3/2
( 5.13) where m is the mass of each sphere and δ 0 = ( F 0 A ) 2/3 representing the initial displacement between neighboring particles.

2D periodic structure

Equations of motion

Consider a 2D periodic structure of identical unit cells subjected to harmonic horizontal base excitation as given in Figures 5.10(a). Each unit cell consists of spherical particle of mass M = 4 3 πR 3 ρ s , with R the radius and ρ s its volumetric mass density. Each particle is coupled with two uniform beams of mass m = bhLρ b with width b, height h, Length L and moment of Inertia I = bh 3 12 . All ends are subjected to compressive forces F 0x and F 0y in both directions x and y respectively 5.10(c). The periodic system is subjected to harmonic horizontal base excitation ÿ(τ ). The main advantage of adding two uniform beams is the ability to control the linear stiffness specific for each mass in a manner to exceed the linear coupling between spheres and ensure the closed linear modes. In addition, the beams helps maintaining the contact between spheres even before applying the pre-compressive loads in both directions. For boundary conditions, we suppose that we have fixed identical systems in both directions, where we define 2(N + M + 2) extra variables and set them to zero as: u 0,j = u N +1,j = u i,0 = u 0,M +1 = 0. The displacement of the i, j th sphere from its equilibrium position in the initially compressed chain is defined by u i,j as m eq üi,j + k eq u i,j + c x ( ui,j -ui-1,j ) + c x ( ui,j -ui+1,j ) +c y ( ui,j -ui,j-1 ) + c y ( ui,j -ui,j+1 ) 5.14) where m eq = M + 13 35 m is the equivalent mass of the two-beams-sphere system, k eq = 192 E b I L 3 is the equivalent beam linear stiffness with E b the Young modulus of the beam. c x = ξ x m eq ω 0 and c y = ξ y m eq ω 0 are the damping coefficients in the x and y directions respectively with ξ x and ξ y the damping ratios. ω 0 = keq meq is the natural frequency. The radius of the elliptical contact area is A = Es √ 2R 3(1-ν 2 s ) where E s and ν s are respectively the Young modulus and the Poisson ration of each spherical particle. δ 0x = F 0x A and δ 0y = F 0y A are the static overlays and Y 0 = m eq Ω 2 y 0 is the harmonic base excitation amplitude (y(τ ) = y 0 cos(Ωτ )).

-A([δ 0x + u i-1,j -u i,j ] 3/2 -[δ 0x + u i,j -u i+1,j ] 3/2 ) -A([δ 0y + u i,j-1 -u i,j ] 3/2 -[δ 0y + u i,j -u i,j+1 ] 3/2 ) = -m eq ÿ (

Weakly nonlinear regime

To approximate the fully nonlinear equations of motion system (5.14), a power series expansion of the forces can be taken. For dynamical displacements with amplitude much less than the static overlay (dynamic displacement have small amplitudes relative to those due to static load), i.e.

| u i-1,j -u i,j δ 0x | << 1 and | u i,j-1 -u i,j δ 0y | << 1 (5.15)
For dynamical displacements with small amplitudes relative to those due to static load, an equivalent spring-mass system can be constructed as shown in Figure 5.11, where a power series expansion was calculated in order to obtain the following i, j th particle equivalent equation of motion. m eq üi,j + k eq u i,j -c x ( ui-1,j -2 ui,j -ui+1,j ) -c y ( ui,j-1 -2 ui,j -ui,j+1 )

-k x (u i-1,j -2u i,j + u i+1,j ) -k y (u i,j+1 -2u i,j + u i,j+1 ) -k 2x [(u i-1,j -u i,j ) 2 -(u i,j -u i+1,j ) 2 ] -k 2y [(u i,j-1 -u i,j ) 2 -(u i,j -u i,j+1 ) 2 ] +k 3x [(u i-1,j -u i,j ) 3 -(u i,j -u i+1,j ) 3 ] + k 3y [(u i,j-1 -u i,j ) 3 -(u i,j -u i,j+1 ) 3 ] = -m eq ÿ = Y 0 cos(Ωτ ) (5.16)
The parameters appearing in equation (5.16) can be written as follows Equation ( 5.16) describes a system of Duffing oscillators linearly coupled by k x and k y in the x and y directions respectively and subjected to harmonic horizontal base excitation.

Normalized equations

It proves convenient to define the following scaled dimensionless variables: 5.17) where u D = Y 0 2(cx+cy)ω 0 is the dynamic displacement of the associated linear system while neglecting the linear coupling with Q = meqω 0 2(cx+cy) . By replacing these variables into the equations of motion and dividing by Y 0 ω 0 meq 2(cx+cy) , we obtain the following nondimensional system: ẍi,j + x i + ξ x (-ẋi-1,j + 2 ẋi,j -ẋi+1,j ) + ξ y (-ẋi,j-1 + 2 ẋi,j -ẋi,j+1 )

τ = ω 0 t x i,j = u i,j u D (
-R x (x i-1,j -2x i,j + x i+1,j ) -R y (x i,j-1 -2x i,j + x i,j+1 ) -β 2x [(x i-1,j -x i,j ) 2 -(x i,j -x i+1,j ) 2 ] -β 2y [(x i,j-1 -x i,j ) 2 -(x i,j -x i,j+1 ) 2 ] +β 3x [(x i-1,j -x i,j ) 3 -(x i,j -x i+1,j ) 3 ] + β 3y [(x i,j-1 -x i,j ) 3 -(x i,j -x i,j+1 ) 3 ] = -2(ξ x + ξ y ) cos( Ω ω 0 t) (5.18)
The parameters appearing in equation ( 5.18) are: In order to create liner closed modes, allowing to study the effects of modes localization on the collective dynamics one must take into account the following assumptions K x << k eq + 2k x + 2k y and K y << k eq + 2k x + 2k y (5.20) In the following section, we will proceed to solve the normalized equivalent differential system (5.18) using a perturbation technique. Before doing so, we should underline the fact that

ξ x = cx meqω 0 , ξ y = cy meqω 0 , R x = kx keq+2kx+2ky , R x = kx keq+2kx+2ky , β 2x = k 2x Y 0 2(cx+cy)ω 3 0 meq , β 2y = k 2y Y 0 2(cx+cy)ω 3 0 meq , β 3x = k 3x Y 2
k 2x k 3x = 6 δ 0x << 1 and k 2y k 3y = 6 δ 0y << 1 (5.21)
Therefore, the quadratic nonlinearities are much smaller than the Duffing coupling parameters, thereby their nonlinear effects are considered negligible. This result is verified in Figure 5.12 where we calculated the displacement of two spheres in contact under compression load, with and without the quadratic nonlinearities.

Analytical-numerical approach

We introduce the parameters listed in equation (5.22) and set the external frequency an amount εω 0 Ω D away from the resonant frequency, whereby they contribute to the equations of amplitude.

ξ x = 1 2 εη x ξ x = 1 2 εη x R x = 1 2 εΓ x R y = 1 2 εΓ y and 2(ξ x + ξ y ) = ε 3/2 g(5.22)
Applying the method of multiple scales combined with the standing waves decomposition as given in the previous in Chapter 3, we obtain the following system of N × M coupled nonlinear equations

2i dAr,p dT + 2i(η x sin 2 ( qr 2 ) + η y sin 2 ( dp 2 ))A r,p + 2(Γ x sin 2 ( qr 2 ) + Γ y sin 2 ( dp 2 ))A r,p +3β 3x sin( qr 2 ) i,k,l j,s,o sin( q i 2 )sin( q k 2 )sin( q l 2 )A i,j A k,s A * l,o ∆ (1) ikl;r ∆ (2) jso;p +3β 3y sin( dp 2 ) i,k,l j,s,o sin( d p 2 )sin( d r 2 )sin( d s 2 )A i,j A k,s A * l,o ∆ (2) ikl;r ∆ (1) jso;p = 2 (N +1)(M +1) g n m
sin (nq r ) sin (md p )e iΩ D T (5.23)

In solving equations (5.23), we write the periodic steady state solution as A r,p = a r,p e iΩ D T , we obtain the following complex amplitudes a r,p .

-2Ω D a r,p + 2i(η x sin 2 ( qr 2 ) + η y sin sin (nq r ) sin (md p ) (5.24)

These complex algebraic equations which define the time independent mode amplitudes are the main result of the perturbation analysis applied to the normalized equivalent differential system. In order to study the collective dynamics of the weakly coupled 2D periodic structure, we start by solving the algebraic system for a set of design parameters listed in Table 5.2 which satisfies the modes localization assumption given by equation (5.20).

Results and discussions

Before starting our investigations we should note the following statements. First, the second member of the equation (5.24) is proportional to the sum of standing waves, which is null for all i,j modes where at least one position is even. For a single direction (N = 1, or M = 1), the general periodic structure is reduced to an array of coupled oscillators. nonlinear coupling parameters in both directions. Figure 5.14 shows the squares of the amplitudes of the different modes as a function of the detuning parameters for the second design parameters in Table 5.2. DM i,j represent the Double Mode solution branches generated by two different modes and QM is the Quadruple Mode involved by the excitation of all modes collectively.

Note that the mode localization has interesting effects on the response intensities with additional multimodal solutions that share the same frequency range with bifurcation topology transfer. In addition, as we choose identical coupling parameters, a 1,2 and a 2,1 become symmetric and interact in the same manner with a 1,1 . Moreover, the Triple Mode (TM) solution branches are missing in this case, as the non-excited modes a 1,2 , a 2,1 or a 2,2 are pairwise symmetric. By choosing different linear and nonlinear coupling parameters, we calculate the response intensities for the third design parameters and we plot them in Figure 5.15.

Remarkably, the DM solution branches are being separated, with additional narrow multimodal solution branches. In addition by decreasing the coupling parameters in the y direction, the DM 1,1;1,2 and QM move away from the resonant branch in |A 1,1 | 2 and their frequency range decreases.

-0. For the fourth set of design parameters, the coupling in the y direction is highly decreased. The responses intensities plotted in Figure 5.16 show that the small solution branches join together in one DM 1,1;1,2 branch heading in the direction of the non-resonant branch of |A 1,1 | 2 , While the response intensity of |A

Summary

In this chapter, we investigated the collective dynamics of nonlinearly coupled oscillators for two applications. The first one is devoted to study electrostatically coupled nanobeams array under parametric excitation, the collective dynamics has been investigated while taking into account the fifth order nonlinear terms. The complexity of the resulting dynamic behavior has been shown, with the appearance of the competition between hardening and softening behaviors for two coupled nanoresonators. In addition, for three coupled nanobeams we can obtain up to five stable solutions for a given frequency, where the basins of attraction have been investigated in the Nyquist plane proving the robustness of the multimodal solution branches. In practice, when the device is used in sensing applications, it is possible to exploit the nonlinearities for the improvement of its sensitivity. The bifurcation topology transfer phenomenon can also be exploited in order to simplify the control of NEMS-based gas sensors.

On the other hand, the study of the collective dynamics of a periodic 2D granular particles array has been studied under a harmonic base excitation, where the nonlinearity is generated by the Hertzian contact law and the applied compressive loads. Complex nonlinear dynamical phenomena with high number of stable multimodal solution branches has been observed, by coupling each particle with two uniform beams allowing the creation of linear closed modes. Unlike the previous physical applications, the choice of the physical parameters can affect the stability of the system, particularly a part of the SM resonant branch.

Several research were devoted to investigate the localization phenomenon whether in pendulums, MEMS/NEMS or granular particles arrays. However, they were limited to the study of a single excitation. Hence, the following Chapter 6 is devoted to study the ILMs in simultaneously driven, damped weakly coupled nonlinear oscillators, where solitary solutions can be found based on the nonlinear Shrödinger equation.

6

Intrinsic localized modes

Introduction

The collective nonlinear dynamics represents an interesting phenomenon that can be encountered in the arrays of weakly coupled nonlinear oscillators. Previous chapters were devoted to the study of the collective dynamics for several physical applications, highlighting the complexity of the responses, the high number of multimodal solutions and other features. The current chapter will focus on a different type of nonlinear behavior, namely Intrinsic Localized Modes (ILMs), also known as Discrete breathers (DBs) or solitons.

In the theory of waves in nonlinear periodic structures, spatial nonlinear localization is one of the most important properties encountered in nonlinear normal modes, providing a link between these modes and solitary solutions. It can be generated either by an extrinsically imposed disorder as in the case of the Anderson localization or by the interaction between the inherent nonlinearities of the resonators. The localization represents an interesting phenomenon in engineering science, which can occur in periodic structures when the wave-function amplitude of the oscillating modal shape is localized in space and decays exponentially. This phenomenon has inspired innovative studies in physics and motivated researchers over many years to explore in depth its effects and consequences. Voluminous studies on mode localization in discrete and continuous periodic nonlinear systems exist in engineering physics using appropriate analytical and numerical techniques. ILMs have been observed in Josephson junctions coupled arrays [Trías 00,Binder 00], antiferromagnet [Schwarz 99,Sato 04], optics [Lederer 08,Sukhorukov 03], photonics [Sievers 88,Christodoulides 04], and atomic lattices [ Bickham 93]. ILMs received a lot of attention in micromechanical resonators arrays [Sato 03c, Sato 03b], Kenig et al. [Kenig 09b] studied the ILMs in arrays of parametrically driven nonlinear oscillators with application to MEMS and NEMS systems. For granular crystals chains, it has been shown that the interplay of periodicity, nonlinearity driving and asymmetry allows the exploration of localization phenomena including solitons and DBs [Huang 93,Theocharis 10,Chong 13].

The following Nonlinear Schrödinger (NLS) equation is used to describe the nonlinear envelop waves as a fundamental model in the previous and other physical applications

iψ t + ψ xx ± β|ψ| 2 ψ = 0 (6.1)
In 1972, Zakharov et al. [Zakharov 72] solved the NLS equation, demonstrating the inerrability and existence of soliton solutions using the Inverse-Scattering Method (ISM). Using the same method, Ablowitz et al. [Ablowitz 73] solved Sine-Gordan (SG) equation given by θ tt -θ xx = sin(θ), (6.2)

and showed that it admits two soliton solutions called kink and anti-kink.

The stability and the bifurcation of the solitons solutions of the damped-driven SG equation in the NLS approximation have been studied analytically by [Terrones 90]. The Parametrically driven, Damped NLS (PDNLS) equation can be written as follows

iψ t + ψ xx + 2ψ|ψ| 2 = he 2it ψ * -iγψ (6.3)
where * denotes the complex conjugate. Exact soliton solutions ψ -and ψ + of the PDNLS equation were found analytically by Barashenkov et al. [Barashenkov 91]. In addition, they demonstrated that in the dissipative case, ψ -solution is unstable for all damping and driving parameters, while ψ + is stable in certain regions. Later, Bondila et al. [Bondila 95] obtained the attractor chart of equation (6.3) on the (γ, h)-plane, identifying stability, period-doubling and quasiperiodic transitions to chaos regions. An approximate analytical expression for the localized solution of the PDNLS equation with positive nonlinear damping, was constructed by Kenig et al. [Kenig 09b].

The externally driven damped NLS equation which can be written as 6.2 Derivation of the amplitude equation

iψ t + ψ xx + 2ψ|ψ| 2 = -iγψ -ge iΩt (6.4)
may describe small amplitude breathers of the long Josephson junction, an acdriven damped SG system or a ferromagnet axis. For zero damping (γ = 0), Barashenkov et al. [Barashenkov 89] derived the exact soliton solutions. For γ > 0, they demonstrated that contrarily to the PDNLS solutions, ψ + is unstable for all γ and h, while ψ -is stable for certain regions. Later, they constructed the existence and stability chart for its soliton solution on the plane of the forcing amplitude g and the dissipation coefficient γ [ Barashenkov 96].

As shown in Chapter 3, section 3.4.1, for two coupled oscillators under primary resonance, we obtained two small stable regions in the additional multimodal solutions, which are hard to find experimentally. However, when adding parametric excitation the system dynamics gains in richness, with additional stable branches over a large frequency range. This remarkable feature motivated the study of ILMs in periodic lattices of nonlinear resonators under simultaneous excitations. The method of multiple scales is employed in order to transform the dimensionless coupled nonlinear equations of motion of the periodic array into a damped driven NLS equation. Exact stationary soliton solutions of the undamped driven NLS equation are derived, while the damped one is numerically solved using the continuous analog of the Newton method. Several numerical simulations have been performed in order to highlight the additional value of employing both external and parametric excitations simultaneously on the stability of localized solutions.

Derivation of the amplitude equation

The normalized equation of motion of an array of coupled D-VDP oscillators up to the third order, under simultaneous external and parametric excitation can be written in the following form

ün -1 2 ε 1 2 D(u n+1 -2u n + u n-1 ) + εγ un + u n + ε ĥ cos 2(1 + εΩ)t u n + λ (u n -u n+1 ) 2 d dt (u n -u n+1 ) + (u n -u n-1 ) 2 d dt (u n -u n-1 ) + ηu 2 n un + δ (u n -u n+1 ) 3 + (u n -u n-1 ) 3 + ξu 3 n = ε 3 2 ĝ cos (1 + εΩ)t (6.5)
The same equation was obtained in Chapter 3, equation (3.9) with different scaling and different coefficients nomenclature. These changes can be resumed as follows

δ -→ η γ -→ D η -→ γ (6.6) ψ -→ δ Ω D -→ Ω ξ -→ ξ (6.7)
As the resonators are collectively oscillating at almost the same frequency, we write the displacement of the n th resonator as [Kenig 09b]:

u n = ε 1 2 [ ψ( Xn , T )e i(t-πn) + c.c.] + ε 3 2 u (1) n (t, T , Xn ) + • • • n = 1, ..., N, (6.8)
with c.c. representing the complex conjugate, T = εt and Xn = ε 1 4 n are slow temporal and spatial variables. Writing the continuous variable X in place of Xn and replacing the displacement solution (6.8) into the normalized equation of motion term by term up to order ε 3 2 we obtain:

un = ε 1 2 ε ∂ ψ ∂ T + i ψ e i(t-πn) + c.c. + ε 3 2 u(1) n , (6.9) ün = ε 1 2 ε 2 ∂ 2 ψ ∂ T 2 + 2iε ∂ ψ ∂ T -ψ e i(t-πn) + c.c. + ε 3 2 ü(1) n , (6.10)
Expanding u n±1 into its Taylor series at X up to the second order it can be written as

u n±1 = -ε 1 2 ψ ± ε 1 4 ∂ ψ ∂ X + ε 1 4 2 ∂ 2 ψ ∂ X2 e i(t-πn) + c.c. + ε 3 2 u (1) n±1 , (6.11)
u n+1 -u n + u n-1 = -ε 3 2 ∂ 2 ψ
∂ X2 e i(t-πn) + c.c., (6.12) -πn) + O(e 3it , e 3iπn ) + c.c, (6.15) -πn) + O(e 3it , e 3iπn ) + c.c, (6. 16) -πn) + O(e 3it , e 3iπn ) + c.c, (6.17) (6.18) where O(e 3it , e 3iπn ) are rapidly oscillating terms with 3 and 3π are the temporal frequency and the spatial wave number respectively. Although, the EoMs (6.5) are trivially satisfied at the order ε 1/2 , one must satisfy a solvability condition at the order ε 3/2 by vanishing all terms proportional to e i(t-πn) so that u

ε cos 2(1 + εΩ)t u n = 1 2 ε 3 2 ψ * e iΩ T e i(t+πn) + O(e 3it ) + c.c, (6.13) cos (1 + εΩ)t = - 1 2 ε 1 2 e i Ω 2 T e i(t-πn) + c.c, (6.14) un u 2 n = ε 3 2 i| ψ| 2 ψe i(t
u 3 n = 3ε 3 2 | ψ| 2 ψe i(t
[(u n -u n+1 ) 3 + (u n -u n-1 ) 3 ] = ε 3 2 48| ψ| 2 ψe i(t
6.3 Soliton solutions [(u n -u n+1 ) 2 ( un -un+1 ) + (u n -u n-1 ) 2 ( un -un-1 )] =ε
(1) n remains finite. Then, we obtain the following partial differential equation (PDE) defining the slow dynamics of the resonators amplitudes

2i ∂ ψ ∂ T + 48 δ +3 ξ +i η +16 λ | ψ| 2 ψ + 1 2 D ∂ 2 ψ ∂X 2 +iγ ψ - ĥ 2 ψ * e 2iΩ T = - ĝ 2 e iΩ T (6.19)
Using the following scaling parameters

ψ = 2Ω 48 δ+3 ξ ψ, X = D 2Ω X, T = 2 Ω T, η + 16 λ = 48 δ+3 ξ 2 η, γ = Ωγ, ĥ = 2Ωh, ĝ = 2Ω 2Ω 48 δ+3 ξ g, (6.20
i ∂ψ ∂T = - ∂ 2 ψ ∂X 2 + (1 -iγ)ψ -(2 + iη)|ψ| 2 ψ + hψ * -g (

Soliton solutions

Exact analytical soliton solutions of the PDNLS equation

Barashenkov et al. [Barashenkov 91] found the exact time-independent soliton solutions of the PDNLS equation (6.22) (for η = 0), which can be written as

Ψ ± (X) = A ± e -iΘ ± sech[A ± (X -X 0 )], (6.23)
where X 0 is an arbitrary position of the soliton, A ± > 0 and

A 2 + = 1 + (h 2 -γ 2 ) 1/2 , 2Θ + = arcsin(γ/h) A 2 -= 1 -(h 2 -γ 2 ) 1/2 , 2Θ -= π -arcsin(γ/h) (6.24)
The pair of solitons solutions ψ + and ψ -exist when h is greater than γ. Barashenkov et al. [ Barashenkov 91] showed that the soliton solution ψ -is unstable for all values of h and γ, while ψ + is stable in certain parameters range. Bondila et al. [Bondila 95] represented the attractor chart of the PDNLS equation solitons solutions on the (γ, h)-plane as shown in figure 6.1. Below the lowest straight line h = γ, the trivial solution ψ ≡ 0 is the only attractor, while above the uppermost line h = 1 + γ 2 it becomes unstable. The Hopf bifurcation is represented by the line 1, where below and to the left of this curve stable stationary solutions exist. [ Bondila 95] on the (γ, h)-plane. Below the lowest straight line h = γ, the trivial solution ψ ≡ 0 is the only attractor, while above the uppermost line h = 1 + γ 2 it becomes unstable. The Hopf bifurcation is represented by the line 1, where below and to the left of this curve stable stationary solutions exist.

Approximate analytical soliton solutions for η > 0

For positive nonlinear damping (η > 0), Kenig et al. [Kenig 09b] constructed an approximate analytical expression for the localized solution of the full damped, parametrically driven NLS equation (6.22) (with g = 0). It is based on writing a function of the same form as Ψ ± given in equation (6.23), with time dependent A and Θ, which is given by .26) which has to be positive and

ψ app (X) = a + e -iθ + sech[a + (X -X 0 )] (6.25) with a 2 ± = 1 -γ η ± h 2 (1 + η2 ) -(γ + η) 2 1 + η2 , ( 6
h cos(2θ ± ) = a 2 ± -1, (6.27) h sin(2θ ± ) = γ + ηa 2 ± .
(6.28)

Where (a + , θ + ) and (a -, θ -) are respectively a stable and a saddle nodes. The saddle-node bifurcation point of these solutions occurs at

h sn = γ + η 1 + η2 , where η = 2 3 η (6.29)
for γ η < 1. An eventual approximation of the n th resonator displacement u n can be obtained after replacing equation (6.25) into equation (6.8)

u n (t) 2 2εΩ 3 a + sech a + 2εΩ D n -X 0 × cos (1 + εΩ)t -πn -θ + (6.30)
In figure 6.2 (a), we plot the scaled analytical approximation (6.25), in the presence of nonlinear damping, while figure 6.2 (b) shows the approximation of the displacements according to the number of dofs n.

Exact analytical solitary solutions of the CSDNLS equation

As a first step, searching for the ground state which has the form of a localized (in space) exact solutions. Given a real solution ψ(x, t) for the following CSDNLS equation (γ = η = 0):

iψ T + ψ XX -Hψ + 2ψ|ψ| 2 = -g, (6.31)

with H = (1 + h). Each ψ(X, T ) = kψ(kX, k 2 T ) is also a solution, corresponding to g = k 3 g and H = k 2 H. Consequently, any solution to equation (6.31) is characterized, up to a simple scaling, by a single combination f = gH -3/2 [ Barashenkov 89].

They found two different soliton solutions of equation (6.31) of the form:

ψ ± (X) = ψ 0 (1 + 2 sinh 2 α 1 ± cosh α cosh(AX)
), (6.32)

Where f is the monotonously decreasing function of h and g defined as [Jallouli 17]:

f = g (1 + h) 3/2 = √ 2 cosh 2 α (1 + 2 cosh 2 α) 3/2 .
(6.33)

f being a monotonically decreasing function , α is uniquely determined by h and g. ψ 0 is the asymptotic value of both ψ -and ψ + solitons: ψ ± (X) → ψ 0 as |X| → ∞ (6.34)

Finally, A has the meaning of "half the area" of both solitons ψ -and ψ + and is equal to (6.35) ψ 0 is real and positive:

A = 2ψ 0 sinh α = 1 2 (ψ ± (X) 2 -ψ 2 0 )dX,
ψ 0 = 1 2(1 + 2 cosh 2 α) = g 4(1 + h) 3/2 cosh 2 α 1/3 (6.36)
We should note that these solutions exist for α ∈ [0, ∞] or for f = g (1 + h) 3/2 ∈ [0, 2/27 0.2722] in terms of excitations amplitudes. Figure 6.3 shows the domain of existence of α according to the different values of the driving forces h and g. Figure 6.4 shows the evolution of the solitary solutions ψ ± of the CSDNLS equation according to α.

Numerical solutions of the DSDNLS equation

For positive linear damping coefficient (γ > 0, η = 0), equation (6.22) is solved numerically using the continuous analog of the Newton method (also known by the variable iteration step Newton method) [Gavurin 58,Galántai 00,Allgower 12]. To solve the PDE (6.22) numerically, we write it on its finite difference form on the discretized domain [-L/2, L/2] as E(ψ) = 0, (6.37) ψ = (ψ 1 , ψ 2 , . . . , ψ N +1 ) discretized solution, with ψ n = ψ(X n ), X n = -L 2 + n∆X, ∆X = L N +1 and E = (E 1 , E 2 , . . . , E N +1 ) is a nonlinear operator defined as follows and .39) This can be obtained using the finite difference approximations for both differential operators ψ XX and ψ X and satisfying the boundary conditions ψ X (±L/2) = 0. The basic concept of the continuous analog Newton's method is to introduce an additional growing variable τ , in such a manner that ψ satisfies the following differential equation d dτ E(ψ(τ )) + E(ψ(τ )) = 0, (6.40)

E n = ψ n+1 -2ψ n + ψ n-1 (
E 0 = -3ψ 0 + 4ψ 1 -ψ 2 2∆X E N +1 = ψ N -1 -4ψ N + 3ψ N +1 2∆X . ( 6
with the initial conditions ψ(0) = ψ (0) . (6.41) ψ (0) is considered to be the exact solitary solutions of the SDNLS equation with no damping. As E(ψ(τ )) -→ τ →+∞ 0, ψ(∞) satisfies equation (6.37).

ψ (k+1) = ψ (k) -∆τ (k+1) ∂E ∂ψ -1 ψ=ψ k E(ψ k )
where k = 1, 2, . . . (6.42) ∆τ (k+1) = τ (k+1) -τ (k) is selected in order to minimize the following residual (6.43) In order to implement the variable iteration step Newton method, we start by choosing a pair (g, h) of parameters in the domain of existence [0, 0.2722]. Equation (6.33) is used to determine the unique α for the chosen (g, h), and we start our continuation using the exact solitary solutions (6.32) as an approximation for γ = 0.01. Then, we calculate the numerical solution (6.42) for the same g, h and γ = 0.01 and keep advancing along the path until convergence is achieved and residual (6.43) is minimized. After, we chose the obtained numerical solution for γ = 0.1 as an approximation for the (g, h) pair with γ = 0.02, and the process repeated until the Newtonian iterations creased to converges for the desired damping parameter γ.

δ (k) = max 1≤n≤N |ReE n (ψ (k) )|, |ImE n (ψ (k) )|
These calculations were performed on the interval (-L/2, L/2), where the solitons solutions decay slowly in space (ψ ± (±L/2) = 0). Finally, for a given ∆x, the obtained ψ is used as a boundary condition for T = 0 to solve the DSDNLS equation (6.22) numerically, using a Runge-Kutta algorithm while setting ψ T = 0. Barashenkov 89] on the plane of both forcing amplitude g and dissipation coefficient γ. They showed that contrarily to the PDNLS soliton solutions, ψ + is unstable for all forcing amplitude g and dissipation coefficient γ. In figure 6.5, blue color indicates the stability region of ψ -, while in the pink region ψ -exist and unstable. Barashenkov et al. [Barashenkov 96] constructed the existence and stability chart for the soliton solutions of the externally driven, damped NLS equation given in figure 6.5. They showed that contrarily to the PDNLS soliton solutions, ψ + is unstable for all forcing amplitude g and dissipation coefficient γ. In figure 6.5, blue color indicates the stability region of ψ -, while in the pink region ψ -exists and it is unstable.

Numerical simulations and interpretations

In this section we are interested in solving numerically the DSDNLS equation (6.22) with zero nonlinear damping (η = 0) for several sets of parameters. According to the existence and stability diagram (6.5) of the externally driven, damped NLS equation constructed by Barashenkov et al [Barashenkov 96], the lower straight line g = 2 π define an approximation for the lower boundary of the domain of existence of ψ ± soliton solutions. Choosing three different pairs of parameters in this region (+ 1), (• 2) and ( 3). For the first two pairs of parameters, ψ ± do not exist for h = 0. However, the effect of adding a certain amount of parametrical driving force, might change the expectations.

Therefore, starting by adding an amount of 0.165 to h for the first point (+ 1) and calculating the numerical solutions of the DSDNLS equation. The top figures in 6.6 show the real and imaginary parts of ψ ± at T = 0, the middle ones show the Temporal response of their imaginary parts over time until convergence, while (e) and (f) show the evolution of their absolute values square |ψ ± | 2 over time and space. Remarkably, both solitons solutions do not only exist but there is one that is stable which is ψ + , while the second one ψ -decays to zero.

In addition for the second pair of parameters (• 2) and for h = 0.05, figure 6.7 shows that similarly both localized solutions exist except that their stability is reversed, whereby ψ + decays to zero while ψ -is stable. For the third point ( 3), we kept the same damping parameter γ = 0.3 while increasing the external parameter g. For these parameters ψ + does not exist while ψ -exists and decay to zero over time as shown in figure 6.8. However, replacing h by the same amount of g = 0.15, numerical simulations given in figure 6.9 show that both localized solutions exists and converge for the same stable solution.

Now choosing the pair of parameters of the fourth point ( * 4) in the blue region of figure 6.5. For these parameters and for h = 0, ψ + does not exist while ψ -is stable as shown in figure 6.10. In contrast, when adding a parametric excitation h = 0.25, both localized solutions ψ ± exist as unexpected interesting results are revealed. Firstly, figure 6.11 shows that ψ -loses its stability and oscillates periodically in time, as shown in the temporal solution and phase portrait curves. In addition, Figures 6.12 (a) and (c) show that the initial transient gives rise to the formation of three solitons ψ + , where the middle one decays to 0 and the two others are not perfectly periodic. Barashenkov et al. [ Barashenkov 96] demonstrated that the stability of solitary solutions is very sensitive to the interval length. Therefore, increasing L to 100 which corresponds to N = 501 dofs, the localized solutions are given in figure 6.12. It is noteworthy that the number of solitons solutions increases up to 6, in addition to their stabilization. Note that before reaching the steady state, several pairs of solitons emerges into one.

The last configuration in figure 6.5 corresponds to the fifth point (• 5) in the pink region where ψ + does not exist and ψ -is unstable and represents a spatio-temporal chaos as shown in figure 6. 13 (a) for h = 0. Figure 6.13 shows the evolution of the localized solution |ψ -| 2 over time for several values of h. We can see that adding a parametric excitation to the system enables to avoid the collapse of spatiotemporal chaos.

Conclusion

We investigated the intrinsic localization in a one dimensional array of weakly coupled D-VDP oscillators under simultaneous parametric and external excitations. The multiple scales method was employed in order to transform the differential system into an amplitude shrödinger equation to describe the spatio-temporal dynamics of the system. The state of art of the soliton solutions for the PDNLS equation with and without damping were given. For zero dissipation, the analytical soliton solutions of the simultaneously driven NLS equation were determined. The DSDNLS equation was numerically solved as a boundary value problem over an interval of length L, using the continuous analog Newton's method for zero nonlinear damping. Several simulations on different pairs of parameters were performed, based on the existence and stability domain of the externally driven, damped NLS equation. It has been shown that adding an amount of parametric driving force enables the mod-ification of the existence and stability chart of the localized solutions and thus avoid spatiotemporel chaos and extend the stability domain according to the linear damping. Therefore, the combination between parametric and the external excitations, helps tuning the existence and stability of solitons solutions.

In practice, and since the damping parameter is often imposed by the system, this study can serve as a predicting numerical tool; allowing the control of the existence and stability of localized solutions by the simultaneous driving forces to localize energy or avoid energy localization.

Conclusions

Summary

This thesis has detailed the importance of periodic systems and their use in engineering devices and physical sciences. The methods employed and the interesting phenomena encountered in linear periodic arrays are given. From linear to nonlinear periodic arrays, complex dynamical phenomena not found in the linear systems may arise. Several physical applications can be modeled as coupled nonlinear oscillators arrays, which can be grouped into two different categories; which are the nonlinear oscillators arrays with weakly linear and nonlinear coupling. Localization is an important phenomena that can be encountered in weakly coupled quasi periodic linear arrays or weakly coupled perfectly periodic nonlinear arrays.

Firstly, various analytical and numerical methods suitable to investigate the nonlinear dynamics of periodic arrays are detailed. The time integration method, is an easy classical numerical solving procedure, generally fails to capture unstable solutions and bifurcation points. In addition, for weakly damped systems, shooting methods succeed to replace the direct time integration method and reduce the computational time. However, this later is much more complicated to program. The ANM coupled with HBM is a powerful tool to solve weakly nonlinear systems, by writing the nonlinear system in its quadratic form. However, this method presents several difficulties when dealing with an important number of weakly coupled dofs. Particularly, it may not converges, as we should provide several initial points and gather the curves specially when having disconnected branches. Finally, the multiple scales perturbation technique coupled with the standing waves decomposition, gives an analytical resolution of the weakly coupled nonlinear system. In addition, it is able to identify the modal interactions and the bifurcation topologies.

The collective dynamics of a general model of weakly coupled nonlinear D-VDP oscillators up to the fifth order nonlinearity has been proposed and modeled under simultaneous primary and parametric resonances. The multiple scales method coupled with standing waves decomposition has been employed, transforming the nonlinear differential system into a complex algebraic one, which was solved numerically using the ANM after applying a quadratic transformation. With nonlinearities up to the third order, the cases of two and three coupled oscillators were investigated for sets of parameters. We have shown that adding a parametric excitation to externally excited system improve its stability with additional stable branches over a large fre-quency range. In addition, for few coupled oscillators, we obtain complex response intensities with a high number of multimodal solutions, which can serve as a hint for the high number of additional solutions for a high number of dofs. Besides, the study of basins of attraction showed a large distribution in the Nyquist plane for the additional multimodal solutions, which validates the robustness of their attractors.

For periodic nonlinear oscillators with weak linear coupling, pendulums arrays represent a famous example in mechanics, which can be described by the SG equation. The collective dynamics of weakly coupled pendulums arrays, has been studied under external excitation. As the coupling is sufficiently weak, the multimodal solutions are distributed on the frequency responses between the resonant and nonresonant branches. In addition and since all even odd modes are excited after modes projection, a ROM approach has been applied in order to solve the EOMs for a high number of rotational dofs. It is based on projecting on odd modes, while maintaining the dominant dynamics of the responses without significant loss of accuracy compared to the results of the full model. In addition, when the number of coupled pendulums increases, the distribution of resonant branches increases.

From 1D to 2D periodic arrays of weakly coupled pendulums, the same solving procedure was employed to study the differential system under a harmonic base excitation. We have shown that the complexity of the responses depends on the linear coupling parameters in both directions. Indeed, when these parameters are different, additional features distinct from those presented in a single array may arise. Particularly, we find multimodal solutions distinguishable from the whole response with larger bandwidth and higher vibration amplitude.

We presented two applications that can be modeled as a periodic nonlinear array with weak nonlinear coupling. An array of electrostatically coupled nanobeams, under parametric excitation where we considered the fifth order electrostatic coupling parameters. The mixed hardening-softening behavior was highlighted, for few coupled nanobeams, regarding to the considered physical parameters.

A 2D periodic granular particles under pre-compression loads in both directions was modeled as a nonlinear oscillators arrays with weak nonlinear coupling governed by the Hertzian contact between particles. In addition, we have ensured the linear closed modes assumption by proposing additional uniform beams coupled to each particle.

Finally, we were interested in determining the ILMs in weakly coupled oscillators arrays under simultaneous parametric and external excitations. The multiple scales method was employed in order to transform the differential nonlinear system into an amplitude NLS equation. The analytical solitary solutions have been found for the undamped system, while the dissipative one was solved numerically using the analog Newton's method with the solitary solution as an initial guess. We have investigated the influence of simultaneous excitations on the stability of localized solutions. Remarkably, the addition of parametric excitation not only gives additional stable multimodal branches as seen in Chapter 3, but it helps also expanding the existence and stability domain of soliton solutions according to the linear damping. magnetic levitation [ Abed 16]. Under appropriate conditions, such process can use solitons to canalize the energy injected into the system following a path defined by an ILM. Each mode determines the type of transduction and its location to convert this canalized vibration energy into electrical one.

Dynamique collective des réseaux périodiques non-linéaires faiblement couplés

Bien que la dynamique des réseaux périodiques non-linéaires ait été investiguée dans les domaines temporel et fréquentiel, il existe un réel besoin d'identifier des relations pratiques avec le phénomène de la localisation d'énergie en termes d'interactions modales et topologies de bifurcation. L'objectif principal de cette thèse consiste à exploiter le phénomène de la localisation pour modéliser la dynamique collective d'un réseau périodique de résonateurs non-linéaires faiblement couplés.

Un modèle analytico-numérique a été développé pour étudier la dynamique collective d'un réseau périodique d'oscillateurs non-linéaires couplés sous excitations simultanées primaire et paramétrique, où les interactions modales, les topologies de bifurcations et les bassins d'attraction ont été analysés. Des réseaux de pendules et de nano-poutres couplés électrostatiquement ont été investigués sous excitation extérieure et paramétrique, respectivement. Il a été démontré qu'en augmentant le nombre d'oscillateurs, le nombre de solutions multimodales et la distribution des bassins d'attraction des branches résonantes augmentent. Ce modèle a été étendu pour investiguer la dynamique collective des réseaux 2D de pendules couplés et de billes sphériques en compression sous excitation à la base, où la dynamique collective est plus riche avec des amplitudes de vibration plus importantes et des bandes passantes plus larges. Une deuxième investigation de cette thèse consiste à identifier les solitons associés à la dynamique collective d'un réseau périodique et d'étudier sa stabilité.

Mots clés: réseaux périodiques, oscillateurs non-linéaires, dynamique collective, couplage faible, interactions modales, bassins d'attraction, localisation d'énergie, solitons.

Collective dynamics of weakly coupled nonlinear periodic structures

Although the dynamics of periodic nonlinear lattices was thoroughly investigated in the frequency and time-space domains, there is a real need to perform profound analysis of the collective dynamics of such systems in order to identify practical relations with the nonlinear energy localization phenomenon in terms of modal interactions and bifurcation topologies. The principal goal of this thesis consists in exploring the localization phenomenon for modeling the collective dynamics of periodic arrays of weakly coupled nonlinear resonators.

An analytico-numerical model has been developed in order to study the collective dynamics of a periodic coupled nonlinear oscillators array under simultaneous primary and parametric excitations, where the bifurcation topologies, the modal interactions and the basins of attraction have been analyzed. Arrays of coupled pendulums and electrostatically coupled nanobeams under external and parametric excitations respectively were considered. It is shown that by increasing the number of coupled oscillators, the number of multimodal solutions and the distribution of the basins of attraction of the resonant solutions increase. The model was extended to investigate the collective dynamics of periodic nonlinear 2D arrays of coupled pendulums and spherical particles under base excitation, leading to additional features, mainly larger bandwidth and important vibrational amplitudes. A second investigation of this thesis consists in identifying the solitons associated to the collective nonlinear dynamics of the considered arrays of periodic structures and the study of their stability.

Keywords: Periodic lattices, nonlinear oscillators, collective dynamics, weak coupling, modal interactions, basins of attraction, localization phenomenon, solitons.

1

  Dynamics of periodic structuresContents 1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.2 Periodic structures . . . . . . . . . . . . . . . . . . . . . . . 23 1.3 Linear systems . . . . . . . . . . . . . . . . . . . . . . . . . 24 1.3.1 Methods for dealing with linear periodic systems . . . . . 24 1.3.2 Multidimensional linear periodic systems . . . . . . . . . . 25 1.4 Nonlinearity in periodic systems . . . . . . . . . . . . . . 27 1.4.1 Pendulums arrays . . . . . . . . . . . . . . . . . . . . . . . 28 1.4.2 MEMS/NEMS arrays . . . . . . . . . . . . . . . . . . . . 29 1.4.3 Granular crystal lattices . . . . . . . . . . . . . . . . . . . 31 1.5 Localization phenomena . . . . . . . . . . . . . . . . . . . 33 1.5.1 History of localization studies . . . . . . . . . . . . . . . . 33 1.5.2 Localization in linear quasi-periodic structures . . . . . . 35 1.5.3 Nonlinear localization . . . . . . . . . . . . . . . . . . . . 37 1.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

Figure 1

 1 Figure 1.1: (a) Crystal lattice of diamond or silicon( Image credit: Dr. Willem Tjerkstra, et al., MESA+ Institute for Nanotechnology, University of Twente). (b) a-Photograph of the Weevil Lamprocyphus Augustus b-Optical micrograph of individual scales attached to the exoskeleton under white-light illumination c-Cross-sectional SEM image of a single scale d-Detailed cross-sectional SEM image of a region of a scale (Reproduced from Phys Rev E 2008;77:050904). (c) Honeycomb. (d) Morpho rhetenor nanostructures (Reproduced from www.asknature.org/media/image/28125#)

  Figure 1.2: Numerous periodic structures among man-made structures.

Figure 1 . 3 :Figure 1 . 4 :

 1314 Figure 1.3: (a) A schematic of a periodic beam grillage system with a single period highlighted The point A represents the reference node [Langley 97]. (b) Generic 3D periodic cells [Collet 11]

Figure 1 . 5 :

 15 Figure 1.5: (a) Bifurcation diagram for the logistic map Jordan Pierce. (b) Spatiotemporal chaos http://www.physics.uci.edu/ dennin/pictures.html. (c) 3D view of the spatio-temporal evolution of a Peregrine soliton Christophe Finot

Figure 1

 1 Figure 1.6: (a) Pictures of three stationary states of the pendula chain obtained for one single given driving amplitude and frequency [Khomeriki 08]. (b)-(c) Experimental setup and schematic diagram of a pendulum array and its frame sits on a platform that can freely move horizontally [Cuevas-Maraver 14]

Figure 1

 1 Figure 1.7: (a) Observation of an inter-site breather [Cuevas 09]. (b) Soliton visualized along a pendulum. Image: pe.soliton.free.fr/.

Figure 1

 1 Figure 1.8: (a) A series of silicon nanowires from [Carr 99]. (b) A scanning electron micrograph showing a fluidics structure fabricated with the sacrificial layer removal technique [Craighead 00]

Figure 1

 1 Figure 1.9: (a) Array of doubly clamped beams forming a diffraction grating [Buks 02]. (b) Micro-electromechanical array; actuation and dissipative forces applied at mid-span of each beam [Gutschmidt 12]

Figure 1 . 10 :

 110 Figure 1.10: One-dimensional chain of particles compressed by a large static force F 0 [Nesterenko 01], causing initial displacement δ 0 between neighboring centers.

  studied the stress wave properties in a basic 2D granular crystals. Coste et al. [Coste 08] investigated the effects of varying the isotropic static compression in linear and nonlinear regimes for 2D granular lattice (figure 1.11 (a)). Experimental and numerical studies have been performed in order to show the impact of the type of excitation imparted on a centered square of highly nonlinear granular system [Leonard 12] (figure 1.11 1.5 Localization phenomena (b)).

Figure 1

 1 Figure 1.11: (a) 2D triangular lattice of beads [Coste 08]. (b) Uncompressed 20 by 20 square array of spherical particles square array of cylindrical particles (intruders) in the interstitial spacings [Leonard 12]

Figure 1 . 12 :

 112 Figure1.12: Normal modes of an isotopically disordered harmonic chain, half-heavy atoms and half•light with a mass ratio of three. The mode numbering is from the lowest to the highest frequency. (1) N=50, (2) N=200[Matsuda 70]

Figure 1 . 13 :

 113 Figure 1.13: Coupled pendula arranged on a chain [Hodges 82].

1. 5 . 3 . 2

 532 Localization in weakly coupled disordered nonlinear systemsKuske et al. [Kuske 93] developed a quasi-degenerate perturbation theory to compute the localization length of the localized wave function while considering the stationary Schrödinger equation in one-dimensional chain of particles with random potentials.

1. 5 . 3 . 3

 533 Nonlinear spatial localization in optics and photonicsThe concept of DBs or ILMs and Discrete Solitons (DS) has been the subject of theoretical, numerical and experimental researchers which attracted increasing attention over the last years. The discovery of such modes and their important role played by nonlinear spatial localization motivated researchers in several domains toward a better understanding of their role. For instance, DBs and ILMs attracted considerable attention in optics due to their novel physics and light-routing applications[Christodoulides 03, Campbell 04]. For a comprehensive review of the theoretical and experimental developments in the area of discrete optical systems which provides new possibilities for optical solitons see [Lederer 08]. In addition, the first experimental observation of a 2D discrete solitons in a biased photorefractive crystals was made by Fleischer et al. [Fleischer 03]. In photonics, Martin et al. [Martin 04] observed 2D discrete solitons in optically induced partially coherent photonic lattices along with a host of new phenomena arising from soliton-lattice interaction (see figure 1.14). Experiments were performed in spatially random nonlinear optical media expanding light propagation, demonstrating that weak interaction tend to restore the coherence in Anderson insulators [Lahini 08]. Maldovan et al. [Maldovan 06] demonstrated theoretically the simultaneous localization of photons and phonons in the same spatial region by introducing lattice defects in a periodic array of dielectric/elastic material that exhibits gaps for both electromagnetic and elastic waves. Schwartz et al. [Schwartz 07] reported experimentally the first observation of Anderson localization of expanding waves packets in spatially random nonlinear 2D optical media. An impressive number of various experiments during the last years observed the existence of ILMs in driven arrays of micromechanical arrays [Sato 06, Kenig 09b].

Figure 1

 1 Figure 1.14: A two-dimensional intrinsic localized mode forms in a photonic lattice [Martin 04].

Figure 2 . 3 :

 23 Figure 2.3: Periodic nonlinear two coupled mass-spring system under primary resonance.

  52) where u n describes the deviation of the nth resonator from its equilibrium n = 1, 2, with fixed boundary conditions u 0 = u 3 = 0. ω 0 = k m define the natural frequency of the uncoupled oscillators, m is its effective mass, k c = d k is the coupling spring constant, µ = c m with c the linear damping, α = k 3 m with k 3 is the nonlinear Duffing parameter and F = f m define the external excitation amplitude.

Figure 2 . 4 :

 24 Figure 2.4: Direct time integration procedure.

Figure 2 . 5 :

 25 Figure 2.5: Forced Duffing frequency response

  Periodic solution for Ω = 37.4

Figure 2 .

 2 Figure 2.6: Amplitudes of the harmonics 1 and 3 of the periodic solution as a function of the frequency Ω.

Figure 2 .

 2 Figure 2.6 shows the frequency amplitudes of the response obtained with the HBM combined with the ANM. The figures show the amplitudes of the harmonics 1 and 3 as a function of the frequency Ω, where solid lines indicate stable solutions and dashed lines indicate unstable ones. We didn't plot the second odd harmonic since it is null.

Figure 2

 2 Figure 2.7: (a) Frequency response curve for Ω ω 0 . (b) Intensity response as a function of frequency.

Figure 2 .

 2 Figure 2.8: (a) Frequency response curve for Ω ω 0 . (b) Intensity response as a function of frequency.

Figure 2 .

 2 Figure 2.9: Response intensity of two resonators as a function of the detuning parameter Ω D .

Figure 2 . 10 :

 210 Figure 2.10: Frequenc response curve.

Figure 2 . 11 :

 211 Figure 2.11: Evolution of the basins of attraction for the first intensity response a 1 with α 1 = α 2 and β 1 = β 2 with respect to several detuning parameter in the Nyquist plane (α 1 , β 1 ).

Figure 3 .

 3 Figure 3.1: A periodic nonlinear system of N+2 coupled D-VDP oscillators, under simultaneous primary and parametric excitations.

5 2

 5 we gets N equations of the form: secular term)e it + other terms(3.13)

Figure 3 . 2 :

 32 Figure 3.2: (a) Response intensity of two resonators as a function of the detuning parameter Ω D , under primary resonance (G = 0.08), for the first design parameters. (b) Zooming in and highlighting the responses areas. Black curves indicate stable solutions. The Single and Double-Mode solution branches are labeled SM and DM respectively.

Figure 3 . 3 :

 33 Figure 3.3: Response intensity of two coupled oscillators as a function of the detuning parameter Ω D , under simultaneous primary and parametric resonances (G = 0.1, H = 7), for the first design parameters. The Single and Double-Mode solution branches are labeled SMi and DMi respectively, with i ∈ {1, 2}

Figure 3 . 4 :Figure 3 . 5 :

 3435 Figure 3.4: Response intensity of two coupled oscillators as a function of the detuning parameter Ω D , under simultaneous primary and parametric resonances (G = 0.08, H = 11), for the first design parameters. SMRB1 and SMRB2 are Single Mode Resonant Branches due respectively to primary and parametric resonances, SMNRB and DMNRB are respectively Single and Double Mode Non Resonant Branches.

Figure 3 Figure 3

 33 Figure 3.6: Distribution of the basins of attraction for the two responses |a 1 | 2 (at the top) and |a 2 | 2 (in the bottom) for Ω D = 0.255 in the planes (α 1 (0), α 2 (0)) and (β 1 (0), β 2 (0)).

Figures

  (a) and (b) display the basins of attraction plotted respectively in the planes (

Figure 3 . 8 :

 38 Figure 3.8: Evolution of the basins of attraction for the first intensity response with α 2 (0) = 0.01, β 2 (0) = 0.02 with respect to the detuning parameter in the Nyquist plane.

Figure 3 Figure 3

 33 Figure 3.9: Evolution of the basins of attraction for the first intensity response with α 2 (0) = 0.05, = β 2 (0) = 0.02 and Ω = 0.255 with respect to the parametric excitation amplitude in the Nyquist plane.

Figure 3 .

 3 Figure 3.10: (a) Response intensity of three coupled oscillators as a function of the detuning parameter Ω D , under simultaneous primary and parametric resonances (G = 0.4, H = 20), for the second design parameters. The Double and Triple Mode solution branches are denoted by: Di j and T k l for i ∈ {1, 3}, j ∈ {1, 2, 3, 4} and k, l ∈ {1, 2, 3} and respectively. DM corresponds to the trivial null solution of equation (3.16), for N = 3 and m = 2. (b) Zooming and highlighting over a region that contains multimode solution branches that are completely disconnected from all other branches.

Figure 3 .

 3 Figure 3.10 displays the response intensity of three coupled oscillators as a function of the detuning parameter Ω D , under simultaneous parametric and external excitations (H = 20 and G = 0.4), for the second set of design parameters given in Table3.1. The curves show the squares of the amplitudes of the three different modes, where DM corresponds to the trivial null solution of equation(3.16), for N = 3 and m = 2. Remarkably, all solutions of |a 2 1 | and |a 2 3 | are multimodal. The Double and Triple mode solution branches, result from the coupling between the first and third oscillators, and the three coupled oscillators respectively. They are denoted by: Di j and T k l for i ∈ {1, 3}, j ∈ {1, 2, 3, 4} and k, l ∈ {1, 2, 3}. Thus, in this case the periodic structure is completely driven by the collective dynamics due to the modal interactions between the nonlinear oscillators. This is illustrated in Figure3.11, where we plotted the average response intensity, defined in equation(3.24), which isI = 4 3 (|a 2 1 | + |a 2 2 | + |a 2 3 |) for N = 3.

Figure 3 .

 3 Figure 3.11: (a) Averaged Response intensity defined in equation(3.24). (b) zooming and highlighting over a region that contains up to five stable solutions. DM and TM denote the branches due to Double and Triple Modes respectively.

  Basins of attraction of |a 1 | 2 and |a 3 | 2 . Basins of attraction of |a 2 | 2 .

Figure 3 . 12 :

 312 Figure 3.12: Distribution of the basins of attraction for the three responses |a 1 | 2 , |a 2 | 2 and |a 3 | 2 for Ω D = 0.46 in the Nyquist plane (α 1 (0), β 1 (0)), with the random initial conditions α 2 (0) = 0.06, β 2 (0) = 0.48, α 3 (0) = 0.38 and β 3 (0) = 0.39.

4. 2 . 2 Figure 4

 224 Figure 4.1: A chain of externally excited equidistant pendulums sharing a common axle.Each pendulum consist of a rigid rod, attached perpendicularly to the axle with a mass at the end. k represents the linear coupling torque constant. Each pendulum is subjected to an external harmonic excitation F cos(Ωτ ) and φ n is the angle between the n th pendulum and the downward vertical.

Figure 4

 4 Figure 4.1 shows the theoretical model of an array composed by an horizontal axle, containing N identical pendulums separated by an equal distance. Each pendulum n is attached perpendicularly to the axle, with a mass m and length l, φ n is the rotational dof of the pendulum measured from its stable equilibrium position. The pendulums numbers "n" and "n + 1" are connected together by a linear torsional spring k. Each pendulum is subjected to an external excitation F cos(Ωτ ).The kinetic energy T of the system can be expressed as

  .6) With L the Lagrangian operator defined by L = T -V and Q n relative to external forces. Substituting equations (4.3), (4.4) and (4.5) into the Lagrange's equation (4.6) leads to the following differential equation of the n th pendulum ml 2 φn + αl 2 φn + mgl sin(φ n ) -k(φ n-1 -2φ n + φ n+1 ) = F cos(Ωτ ) n = 1, . . . , N,

Figure 4 . 2 :

 42 Figure 4.2: (a) Response frequency for one pendulum (N=1) using direct time integration method between the full nonlinear problem (4.7) and the one obtained after applying Taylor series (4.11). (b) Comparison between sin(φ) and its Taylor series φ -1 6 φ 3 as a function of the rotational amplitude φ up to 1.4 Rad.

Table 4 . 1 :

 41 Design parameters for the corresponding periodic structure depicted in Figure 4.1 Parameters m (Kg) g (m.s -2 ) l (m) k (N.m) α (Kg.s -1 ) F (N.m)

Figure 4 . 3 :

 43 Figure 4.3: Rotational dofs of three coupled pendulums as a function of frequency for the design parameters listed in Table4.1. Comparison between a direct time integration method using Runge-Kutta which is employed to solve the approximate differential system (4.7), and the semi-analytical method.

Figure 4 . 4 :

 44 Figure 4.4: (a) Response intensity of three coupled pendulums under external excitation as a function of frequency. Solid curves indicate stable solutions and dashed curves indicate unstable solutions. DM i represent Double Modes solution branches, generated by the first and the thirst mode where the second one is null. T M j are the Triple Mode solution branches involved by exciting all modes respectively. (b) Basins of attractions of these responses for f = 1.83Hz in the Nyquist plane (α 1 (0), β 1 (0)) and random variable α 2 = 0.1, β 2 = 0, α 3 = 0.5 and β 3 = -0.3. Each color reflects the distribution of a specific multimodal solution.

Figure 4 . 5 :

 45 Figure 4.5: Rotational dof as a function of frequency of three coupled pendulums, where all pendulums are excited with external forces, where we zoom and highlight overs a couple of multimodal areas. Black, red and gray curves represent the full model, the reduced order model (The reduction basis contains odd modes only) and the stable branches.

Figure 4 .

 4 Figure 4.6: (a) Rotational dof as a function of frequency of six coupled pendulums, where all pendulums are excited with external forces. Blue curves represent the full model and the Red curves represent the reduced order model (The reduction basis contains odd modes only). (b) Modal intensities of odd modes for the case of six coupled pendulums under external excitation after model order reduction as a function of frequency. Red, blue and gray curves represent the contribution of the TM solution branches in the frequency responses as Resonant (R), Non-Resonant (NR) and stable branches.

Figure 4 .

 4 6 (a) shows the rotational dofs as a function of frequency for six weakly coupled pendulums, where each pendulum is excited with an external force. Blue curves represent the full model (T M +6 th M ) and the Red curves represent the ROM (TM). Results of this model reduction show that the modal reduction preserves the accuracy of the full model in the prediction of all responses. Nevertheless, as the multimodal TM solutions join either resonant or non-resonant branches, further informations about the modal responses are required.

Figure 4 . 7 :

 47 Figure 4.7: Rotational dof as a function of frequency of twelve coupled pendulums, where all pendulums are excited with external forces. Curves represent the reduced order model where the reduction basis contains odd modes only (6 modes).

Figure 4 .

 4 Figure 4.8 displays the basins of attractions of |φ 1 | for f = 1.83Hz in the phase portrait (φ 1 (0), φ1 (0)) for three, four, five and six coupled pendulums simultaneously. Red and blue colors denotes respectively the Resonant (R) and the Non-Resonant (NR) branches of the frequency response of the first dof |φ 1 |. As shown by the sequence of graphs, the distribution of the basins of attraction of the resonant branches (Red areas) increases slightly while increasing the number of coupled pendulums. Although, the study of basins of attraction in term of amplitude (R/NR) is important, we need to investigate the relationship between the type of these amplitudes and the nature of contributed modes.

Figure 4 . 8 :

 48 Figure 4.8: Basins of attraction of |φ 1 | for f = 1.83Hz in the phase portrait (φ 1 (0), φ1 (0)) for three, four, five and six coupled pendulums. Red and blue colors indicate respectively the Resonant (R) and the Non-Resonant (NR) branches in the frequency response of the first dof |φ 1 |.

  Figure 4.9: Basins of attraction of |φ 1 | for f = 1.83Hz in the phase portrait (φ 1 (0), φ1 (0)) for three, four, five and six coupled pendulums. Each color corresponds to different types of mode and branch, they are illustrated in Table 4.2.

Figure 4 .

 4 10(b) shows a front of view of system 4.10.

Figure 4 .

 4 Figure 4.10: (a) 2D array of coupled pendulums subjected to horizontal harmonic base excitation. (b) Front view of the system

Table 4 . 3 :

 43 Design parameters for the corresponding periodic structure depicted in figure 4.10 Parameters m (Kg) g (m.s -2 ) l (m) α (Kg.s -1 ) y 0 (m.s -2 ) k x (N.m) k y (N.m)

Figure 4 . 11 :

 411 Figure 4.11: Response intensity of 2 × 2 coupled pendulums as a function of frequency, under harmonic horizontal base excitation for the first set of design parameters listed in table 4.3.

Figure 4 . 12 : 3 Figure 4 . 13 :

 4123413 Figure 4.12: Rotational dof as a function of frequency of 2 × 2 coupled pendulums, under harmonic horizontal base excitation for the first set of design parameters listed in table 4.3

  .0994.

Figure 4 .

 4 13 shows the response intensities of 2 × 2 coupled pendulums. Unlike to the first case, linear springs k x and k y play different roles. In particular, the multimodal solution branches resulting from the interaction between

Figure 4 .

 4 Figure 4.14: Rotational dof as a function of frequency of 2 × 2 coupled pendulums, under harmonic horizontal base excitation for the second set of design parameters listed in table 4.3

Figure 5 . 1 :

 51 Figure 5.1: Array of N clamped-clamped silicon Micro-or Nano-beams actuated by an electric load S. Gutschmidt et al. Nonlinear Dyn. (2012).

2 = H 2 ( 5 . 4 )

 2254 Ph.D. Thesis -D. Bitar

Figure 5 . 2 :

 52 Figure 5.2: (b) Frequency responses of two coupled nanobeams under parametric excitation for the first set of design parameters. Dashed curves represent unstable solutions. The only elliptical SM solution of the analytical equation (5.4) represented in blue for the second mode amplitude. DM and TM multimodal solution branches are represented in red and green respectively. (a) Zooming and highlighting the multimodal interactions and bifurcation topology transfer between both modes amplitudes. (c) The average response intensity, defined in equation (5.5) for two coupled oscillators and the first design parameters.

  Figures 5.2 (a) (b) show the amplitudes of both symmetric and anti-symmetric modes (|a 1 | 2 and |a 2 | 2 respectively) as a function of Ω D . Full lines indicate stable solutions and dashed curves represent unstable responses, with blue and red lines represent single mode S1 and double mode S2 solution branches respectively, while black straight lines corresponds to the trivial null solution of

Figure 5 . 3 :

 53 Figure 5.3: Square of modes amplitudes of two-nanobeams system under parametric excitation with respect to Ω D for the 2 nd model parameters. The elliptical SM branches (S1 & S2) solutions of the equation (5.4) for m=1 and 2 respectively.

Figure 5 . 4 :

 54 Figure 5.4: (a)-(b) Basins of attraction distribution for both amplitudes |a 1 | 2 and |a 2 | 2 respectively in the Nyquist plane (Re(a 1 (0)), Im(a 1 (0))) for Re(a 2 (0)) = 0.2, Im(a 2 (0)) = -0.2 and Ω D = 0.06. (c) Distribution of the Basins of attractions in terms of modes type: blue and yellow represent SM and DM respectively. (d) The average response intensity, defined in equation (5.5) for two coupled oscillators and the second design parameters.

Figure 5 . 5 :

 55 Figure 5.5: Response intensities of three coupled nanobeams under parametric excitation, for the second design parameters listed in Table 5.1. Solid and dashed lines indicate stable and unstable steady-state solutions respectively. Black, blue, red and green represent trivial, SM, DM and TM solution branches respectively.

Figure 5 .

 5 3 represents the steady state responses, where the top and the bottom graphs show the amplitude of the symmetric and antisymmetric modes respectively with respect to the detuning parameter Ω D . Solid curves indicate stable responses and dashed curves represent unstable solution branches. Particularly, blue curves represent SM responses (S1 and S2), solution of equation

Figure 5 .

 5 Figure 5.6 shows the evolution of the basins of attraction in term of modes type, while decreasing the gap width from 1 µm to 0.9 µm. As shown by the sequence of graphs, the distribution of the basins of attraction of the single and triple modes increase contrary to the trivial and the DM solutions, by slightly decreasing the gap width.

Figure 5

 5 Figure 5.6: Evolution of the basins of attraction in terms of modes type while decreasing the gap width. The basins are plotted in the Nyquist plane (Re(a 1 (0)), Im(a 1 (0))), for the case of three coupled beams, the second set of design parameters listed in Table 5.1, Re(a 2 (0)) = -0.1, Re(a 3 (0)) = -0.05, Im(a 2 (0)) = Im(a 3 (0)) = 0 and the detuning parameter Ω D = 0.18.

Figure 5 . 7 :

 57 Figure 5.7: Evolution of the basins of attraction in terms of response intensity of the first amplitude while decreasing the length of the beams. The basins are traced in the Nyquist plane (Re(a 1 (0)), Im(a 1 (0))), for the case of three coupled beams, the second set of model parameters, Re(a 2 (0)) = -0.1, Re(a 3 (0)) = -0.05, Im(a 2 (0)) = Im(a 3 (0)) = 0 and the detuning parameter Ω D = 0.18.

Figure 5 . 8 :

 58 Figure 5.8: Contact between two spheres

Figure 5 .

 5 Figure 5.9: 1D array of N spheres of radius R in contact, compressed by a static force F 0 at both ends and securing the initial displacement δ 0 between neighboring centers.

Figure 5 .

 5 Figure 5.10: (a) 2D periodic nonlinear lattice of spheres in contact under compression and subjected to horizontal base excitation. (b) A view of the system from the top. (c) Unit cell i, j.

Figure 5 . 11 :

 511 Figure 5.11: An equivalent 2D Mass-Spring system.

Figure 5 .

 5 Figure 5.12: (a) 1D chain of spherical particles under compression. (b) Response frequencies.

Figure 5 .

 5 Figure 5.14: Response intensities of 2 × 2 weakly coupled oscillators under harmonic base external excitation as a function of the detuning parameter Ω D , for the second set of design parameters listed in table 5.2.

Figure 5 . 15 :

 515 Figure 5.15: Response intensities of 2 × 2 weakly coupled oscillators under harmonic base external excitation as a function of the detuning parameter Ω D , for the third set of design parameters listed in table 5.2.

  Contents 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 153 6.2 Derivation of the amplitude equation . . . . . . . . . . . 155 6.3 Soliton solutions . . . . . . . . . . . . . . . . . . . . . . . . 157 6.3.1 Exact analytical soliton solutions of the PDNLS equation 157 6.3.2 Approximate analytical soliton solutions for η > 0 . . . . 158 6.3.3 Exact analytical solitary solutions of the CSDNLS equation 160 6.3.4 Numerical solutions of the DSDNLS equation . . . . . . . 162 6.4 Numerical simulations and interpretations . . . . . . . . 164 6.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170

3 2

 3 16i| ψ| 2 ψe i(t-πn) wwhitewhitewhitehite + O(e 3it , e 3iπn ) + c.c.

  6.22) Equation (6.22) is a generic form of the nonlinear time dependent Schrödinger equation, including both external and parametric driving forces beside the linear and nonlinear damping parameters.With γ = η = 0 and η = 0 we call equation (6.22) the Conservative or Dissipative Simultaneously Driven NonLinear Schrödinger equations (CSDNLS) or (DSDNLS) respectively.

Figure 6 . 1 :

 61 Figure 6.1: Attractor chart of the PDNLS equation as constructed by Bondila et al.[Bondila 95] on the (γ, h)-plane. Below the lowest straight line h = γ, the trivial solution ψ ≡ 0 is the only attractor, while above the uppermost line h = 1 + γ 2 it becomes unstable. The Hopf bifurcation is represented by the line 1, where below and to the left of this curve stable stationary solutions exist.

Figure 6

 6 Figure 6.2: (a) The approximate soliton solution (6.25) of equation (6.22) (with g = 0). (b) Scaled analytical approximation (6.30) for the resonators displacement. The parameters are γ = 0.75, η = 0.25, h = 0.9, D = 0.25, Ω = 46.174, ε = 0.005, X0 = 34.5 and N = 201.

Figure 6 . 3 :

 63 Figure 6.3: Domain of existence of α according to the driving forces h and g.

Figure 6 . 4 :

 64 Figure 6.4: Undamped solitary solution ψ + and ψ -of the SDNLS equation for several values of α.

Figure 6 . 5 :

 65 Figure 6.5: Existence and stability diagram for the soliton solutions of the externally driven damped NLS equation as constructed by Barashenkov et al. [Barashenkov 89] on the plane of both forcing amplitude g and dissipation coefficient γ. They showed that contrarily to the PDNLS soliton solutions, ψ + is unstable for all forcing amplitude g and dissipation coefficient γ. In figure6.5, blue color indicates the stability region of ψ -, while in the pink region ψ -exist and unstable.

  Figure 6.6: (+ 1) γ = 0.1, g = 0.05, h = 0.165 and L = 40.

  (a) Phase portrait of ψ + (b) Phase portrait of ψ - (c) ψ + (X) at T = 200 (d) ψ -(X) at T = 200 (e) Evolution of this |ψ + | 2 over time (f) Evolution of this |ψ -| 2 over time

Figure 6

 6 Figure 6.7: (• 2) γ = 0.3, g = 0.05, h = 0.3 and L = 40.

Figure 6 . 8 :

 68 Figure 6.8: Evolution of |ψ -| 2 over time for ( 3) γ = 0.3, g = 0.15, h = 0 and L = 40.

  Figure 6.9: ( 3) γ = 0.3, g = 0.15, h = 0.15 and L = 60.

Figure 6 .

 6 Figure 6.10: ( * 4) γ = 0.3, g = 0.232, h = 0 and L = 100.

  (a) Temporal response of Im(ψ + ) (b) Phase portrait of ψ - (c) Evolution of |ψ -| 2 over time (d) 3D plot |ψ -| 2

Figure 6 .

 6 Figure 6.11: ( * 4) γ = 0.3, g = 0.232, h = 0.25 and L = 40.

  (a) 3D plot |ψ + | 2 for L = 40 (b) 3D plot |ψ + | 2 for L = 100 (c) Evolution of |ψ + | 2 over time for L = 40 (d) Evolution of |ψ + | 2 over time for L = 100

Figure 6 .

 6 Figure 6.12: ( * 4) γ = 0.3, g = 0.232, h = 0.25.

Figure 6 . 13 :

 613 Figure 6.13: Evolution of |ψ + | 2 solution over time for the pair of parameters (• 5) γ = 0.1, g = 0.232 and L = 40.

n

  N j,k,l,o,p α j α k α l α o β p ∆ α k β l (1 -δ j,k )∆ k α l + 6 j,k,l,o β j α k α l α 2 o (1 -δ k,l )(1 -δ k,o )(1 -δ l,o )∆ k α l + j,k,l,o,l α j α k α l α o α p (1 -δ j,k )(1 -δ j,l )(1 -δ j,o )(1 -δ j,p )(1 -δ k,l ) N j,k,l α j α k α l α o α p + j,k,l,o,p β j α k α l (1 -δ k,o )(1 -δ k,p )(1 -δ l,o )(1 -δ l,p )(1 -δ o,p )∆ β l v j;j,k (1 -δ j,k )∆ k α l + 6 j,k,l,o β j v o;k,l (1 -δ k,l )(1 -δ k,o )(1 -δ l,o )∆ k α l + j,k,l,o,l α j z k,l;o,p (1 -δ j,k )(1 -δ j,l )(1 -δ j,o )(1 -δ j,p )(1 -δ k,l ) N j,k,l α j α k α l α o j,k,l,o,p β j α k (1 -δ k,o )(1 -δ k,p )(1 -δ l,o )(1 -δ l,p )(1 -δ o,p )∆ École doctorale SPIM 1 rue Claude Goudimel F -25030 Besançon cedexn tél. +33 (0)3 81 66 66 02 n ed-spim@univ-fcomte.fr n www.ed-spim.univ-fcomte.fr

Table 2 .

 2 1: Design parameters for the corresponding periodic two coupled mass-spring system under primary resonance depicted in figure 2.3.

  w 1 , w 2 ] t where the external forces are represented in the constant operator c. To compute the periodic solution branches

		0.25
		0.2
	u 1 , u 2	0.15
		0.1
		0.05
		0	36,6	37	37,4	37,8	38,2
					

Table 3 .

 3 1: Design parameters for the corresponding periodic structure depicted in Figure3.1

	Parameters	Design 1 Design 2
	C(Kg.s -1 )	0.001	0.01
	M (Kg)	0.01	0.01

Table 4 .

 4 2.

Table 4 .

 4 2: Color palette to distinguish the solutions in terms of type of mode and branch.

Table 5 .

 5 1: Model parameters of the nanobeams array illustrated in Figure5.1

	Parameters	L (µm)	g (µm)	h (nm)	b (µm)	Q	V ac (mV )	V dc (mV )
	Design 1	11.5	0.215	150	0.2	50	12	0.1
	Design 2	20	1	250	0.5	50	0.75	0.75

Table 5 .

 5 1. Solid and dashed lines indicate stable and unstable steady-state solutions respectively. Black, blue, red and green represent trivial, SM, DM and TM solution branches respectively.

  2 ( dp 2 ))a r,p + 2(Γ x sin 2 ( qr 2 ) + Γ y sin 2 ( dp

	+3β 3x sin( qr 2 )		i,k,l j,s,o	sin(	q i 2	)sin(	q k 2	)sin(	q l 2	2 ))a r,p jso;p (2) ikl;r ∆ l,o ∆ )a i,j a k,s a * (1)
	+3β 3y sin( dp 2 )	i,k,l j,s,o	sin(	d p 2	)sin(	d r 2	)sin(	d s 2	)a i,j a k,s a * l,o ∆ (2) ikl;r ∆ (1) jso;p
	=	2 (N +1)(M +1) g
							n	m

  2,1 | 2 takes an elliptical form. Response intensities of 2 × 2 weakly coupled oscillators under harmonic base external excitation as a function of the detuning parameter Ω D , for the fourth set of design parameters listed in table 5.2.

			1.6 1.8 1.4				Unstable DM 1,1; 2,1 DM 1,1; 2,2 QM		0.07 0.06			Unstable DM 1,1; 1,2 SM QM
			1.2				DM 1,1; 1,2 SM		0.05
	|A 1,1 | 2	0.8 1					|A 1,2 | 2	0.03 0.04
			0.6
										0.02
			0.4
			0.2							0.01
			-0.03 0	-0.02	-0.01	0	0.01	0.02		-0.03 0	-0.02	-0.01	0	0.01	0.02
					 D							 D
		0.18						0.08
							Unstable					Unstable
		0.16				DM 1,1; 2,1	0.07			DM 1,1; 2,2
							SM						QM
		0.14				QM		0.06			SM
		0.12
									0.05
	| 2		0.1					| 2
	|A 2,1	0.08					|A 2,2	0.04
									0.03
		0.06
		0.04						0.02
		0.02						0.01
			-0.03 0	-0.02	-0.01	0	0.01	0.02		-0.03 0	-0.02	-0.01	0	0.01	0.02
					 D							 D
	Figure 5.16:

  ∆X) 2 -ψ n +2|ψ n | 2 ψ n +iγψ n -hψ * n +g; for n = 1 . . . N(6.38)

Ph.D. Thesis -D. Bitar

Ph.D. Thesis -D. Bitar

2.8.5 Basins of attraction . . . . . . . . . . . . . . . . . . . . .

2.9 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

Dr. Najib Kacem for guiding, supporting and encouraging me throughout the doctoral study. I would like to thank him for his unconditional availability, the time and patience for various technical discussions which led to realize major milestones in my research. My warmest thanks go to the jury members, including the reporters Prof. Claude-Henri Lamarque work and the examiners Dr. Cyril Touzé and Dr. Manuel Collet for their participation to the jury. Moreover, I really appreciate the technical support I received from the permanent members of the Laboratory Scott Cogan, Vincent Walter, Emmanuel Ramasso, Emmanuel Foltete and Joseph Lardies. I'm very appreciative of all the support I have received from Prof. Lamine Boubakkar beside his help on various occasions. I greatly acknowledge the funding received towards my PhD from the cooperation with the Labex ACTION program (contract ANR-11-LABX-01-01). Many thanks for the coordinator's assistants of the Labex ACTION Sandrine Chatrenet and Claudia Laou-huen, the SPIM Doctoral School secretariat Alexandrine Vieillard and the Applied Mechanics Department sec-retariat and management team Isabelle Maillotte-Navarro, Christine Froidevaux and Delphine Travaglini.

Configuration (×10 -3 ) η x η y Γ x Γ y β 3x β 3y g 1 8 8 16 16 Particularly, for two coupled complex algebraic equations, we used Mathematica to solve them numerically, where the stability has been performed based on the eigenvalues of the Jacobian matrix of the differential system (5.23) for each point. Figure 5.13 displays the intensity responses of two weakly coupled oscillators (N = 1, M = 2), as a function of the detuning parameter Ω D for the first design parameters listed in Table 5.2. Black, orange and gray curves indicate the Single Mode (SM), Double Mode (DM) and the unstable solution branches respectively. SM solutions are generated by the null trivial solution a 1,2 = 0, where a 1,1 corresponds to a single forced Duffing oscillator. Contrary to the previous results, the SM solutions lose their stability through a pitchfork bifurcation. Similar results were obtained by Touzé et al. [Touzé 02] were they calculated analytically the instability zone of the SM solution branch for the case of two coupled Duffing oscillators under harmonic excitations. The stable solution becomes the coupled one DM, involving the excitation of both modes collectively. With two weakly coupled oscillators we can obtain up to three stable solutions for a given frequency. Hence, for 2D periodic structures we expect additional complex phenomena not obtained for one dimensional array.

Therefore, we study the case of 2 × 2 weakly coupled oscillators, for several design parameters. Starting by the case where we consider the same linear and

Future work

While this thesis developed a general model for studying the collective dynamics of weakly coupled nonlinear oscillators under simultaneous primary and parametric excitations to functionalize the localization phenomena, several research lines arising from this work should be pursued.

Deterministic model

For large arrays, one must find a reduction method different from the one based on the projection on odd modes in order to reduce the computational time. In regards to the basins of attraction, the dynamical integrity must be analyzed by choosing the right definition of safe basin [Rega 08], choosing appropriate integrity measures to qualify its magnitude and investigating the basins evolution for varying system parameters. Nevertheless, when dealing with high number of dofs, one must rely on different criteria allowing to prove the robustness of the additional multimodal solutions. Therefore, probabilistic data analysis tools can be used to explore and visualize high-dimensional data and allow a qualification observation and relationship discovery such as the principal component analysis (PCA) [Jolliffe 02], k-means clustering [Likas 03], Self-organizing maps (SOM) [Kaski 97] or else.

Otherwise, in order to analyze the nonlinear dynamics of an array of electrostatically coupled nanobeams under simultaneous primary and parametric excitations, the proposed model should be extended to investigate the collective dynamics in a diatomic array of weakly coupled nonlinear oscillators.

Uncertainty propagation and robustness analysis Concerning uncertainties, a collaborative work in progress with a Ph.D. student colleague explores the benefits of uncertainties on the collective dynamics of weakly coupled nonlinear oscillators under external excitation [ Chikhaoui 16]. Nonetheless, one should add the parametric excitation to the system and examine the effects of the combination between the simultaneous excitations and uncertainties.

Regarding the 2D periodic granular crystal lattice and since we added in purpose the coupled uniform beams to produce linear closed modes, one must establish rigorous parametric studies for the design and implementation of an experimental prototype. Thus, in order to investigate the behavior of the complex system, a robust parameter optimization must be developed including sensitivity analysis and uncertainties qualifications.

ILMs Despite the results obtained regarding ILMs, several challenges remain, such as constructing a stability chart for the DSDNLS equation, finding the transition between the continuous solutions and the physical displacements, identifying the localized solutions in simultaneously driven, damped, quasi-periodic nonlinear arrays and finally determining the ILMs in a 2D periodic nonlinear arrays. In practice, we aim to implement an experimental setup such as nonlinear oscillator arrays under [0, 0, 0, 0; 0, 0, 0, 0;

Explicit equations of the solution expansion into the EoM

Substituting the proposed solution (3.11) into the normalized, scaled EOM (3.9) term by term. Up to the order ε

εu 2 n un + sin -q j + q k + q l + q o + q p 2 sin[n(-q j + q k + q l + q o + q p )] εu 2 n un + sin

εu 2 n un + sin

εu 2 n un + sin

εu 2 n un + sin

εu 2 n un + sin -q j -q k + q l + q o + q p 2 sin[n(-q j -q k + q l + q o + q p )]

εu 2 n un + sin -q j + q k + q l -q o + q p 2 sin[n(-q j + q k + q l -q o + q p )] εu 2 n un + sin -q j + q k + q l + q o -q p 2 sin[n(-q j + q k + q l + q o -q p )] εu 2 n un + sin q j -q k -q l + q o + q p 2 sin[n(q j -q k -q l + q o + q p )] εu 2 n un + sin q j -q k + q l -q o + q p 2 sin[n(q j -q k + q l -q o + q p )] εu 2 n un + sin q j -q k + q l + q o -q p 2 sin[n(q j -q k + q l + q o -q p )] εu 2 n un + sin

εu 2 n un + sin

εu 2 n un + sin q j + q k + q l -q o -q p 2 sin[n(q j + q k + q l -q o -q p )] εu 2 n un + sin q j + q k + q l + q o + q p 2 sin[n(q j + q k + q l + q o + q p )]

sin(nq j ) sin(nq k ) sin(nq l )(A j e iω j t + c.c.)(A k e iω k t + c.c.)(A l e iω l t + c.c.)

sin[n(-q j + q k + q l )] + sin[n(q j -q k + q l )] + sin[n(q j + q k -q l )]

sin[n(-q j + q k + q l )] + sin[n(q j -q k + q l )] + sin[n(q j + q k -q l )]

εu 2 n un + sin q j -q k + q l 2 sin[n(q j -q k + q l )] + sin q j + q k -q l 2 sin[n(q j + q k -q l )] εu 2 n un + sin q j + q k + q l 2 sin[n(q j + q k + q l )]

εu 2 n un + sin q j -q k + q l 2 sin[n(q j -q k + q l )] + sin q j + q k -q l 2 sin[n(q j + q k -q l )] εu 2 n un + sin

εu 2 n un (A j e iω j t + c.c.)(A k e iω k t + c.c.)(iω l A l e iω l t + c.c.), (A.8)

A.3 ∆ Kronecker functions

The three ∆ functions in Equation ( 5.23) are defined in terms of Kronecker deltas as A.11) where

A.4 Nonlinear cubic and quintic stiffnesses terms into quadratic

jklop;m and N j,k,l,o,p

jklop;m with n ∈ {1, 2}. Struct., vol. 17, no. 1, page 363-372, 2008.

Scientific production

Publications