M C

R R

 S M

MODAL ANALYSIS OF A CRACKED BEAM

This paper addresses the problem of vibrations of a cracked beam. In general, the motion of such a beam can be very complex. This phenomenon can be attributed to the presence of the non-linearity due to the opening and closing of cracks. The focus of this paper is the modal analysis of a cantilever beam with a transverse edge crack. The non-linearity mentioned above has been modelled as a piecewise-linear system. In an attempt to define effective natural frequencies for this piecewise-linear system, the idea of a ''bilinear frequency'' is utilized. The bilinear frequency is obtained by computing the associated frequencies of each of the linear pieces of the piecewise-linear system. The finite element method is used to obtain the natural frequencies in each linear region.

In order better to understand the essential non-linear dynamics of the cracked beam, a piecewise-linear two-degree-of-freedom model is studied. Perturbation methods are used to obtain the non-linear normal modes of vibration and the associated period of the motion. Results of this piecewise-linear model problem are shown to justify the definition of the bilinear frequency as the effective natural frequency. It is therefore expected that calculating piecewise mode shapes and bilinear frequencies is useful for understanding the dynamics of the infinite degree of freedom cracked beam.

INTRODUCTION

Cracks are present in structures due to various reasons. The presence of a crack could not only cause a local variation in the stiffness, but it could affect the mechanical behavior of the entire structure to a considerable extent. Cracks present in vibrating/rotating components could lead to catastrophic failure. For these reasons, there is a need to understand the dynamics of cracked structures. The vibration characteristics of cracked structures can be useful for an on-line detection of cracks (non-destructive testing) without actually dismantling the structure. In particular, the natural frequencies and mode shapes of cracked beams can provide an insight into the extent of damage. As a representative model of a cracked structure, a cantilever beam with a transverse edge crack is studied in this paper.

In the past there have been considerable attempts to understand the dynamics of a cracked beam [1][2][START_REF] S | Vibrations of beams with a fatigue crack[END_REF][START_REF] C | [END_REF][5][6][7][8][9]. Christides and Barr [1] developed a one-dimensional cracked beam theory at the same level of approximation as the Bernoulli-Euler beam theory. Several assumptions on the displacement, velocity and stress fields are built into this model. The pair of symmetric cracks is always assumed to remain open as the beam is vibrating, so as to avoid the non-linear characteristics of an opening and closing crack. An approximate Galerkin solution to the one-dimensional cracked beam theory was obtained by Shen and Pierre [2]. The comparison functions used in the Galerkin procedure consisted of mode shapes of an uncracked beam. Shen and Chu [START_REF] S | Vibrations of beams with a fatigue crack[END_REF] extended the cracked beam theory to account for the opening and closing of the crack-the so-called ''breathing crack'' model.

A Galerkin procedure was used to obtain the bilinear equation for each vibration mode. The non-linear dynamic response of the bilinear equation to a forcing excitation was calculated through a numerical analysis.

Chu and Shen [START_REF] C | [END_REF] obtained a closed form solution for a forced single-degree-of-freedom bilinear oscillator under low frequency excitation. They extended the procedure in order to study the dynamics of cracked beams with binlinear forcing functions. Yuen [5] proposed that the change in the stiffness of the cracked beam at the location of the crack can be modelled as a change in the modulus of elasticity of the cracked location. The finite element method was used to carry out the analysis. Shen and Taylor [6] developed an identification procedure for an on-line detection of the size and location of cracks. A mean square difference and a mini-max criterion were used to demonstrate the reliability of the identification procedure.

Ostachowicz and Krawczuk [7] replaced the crack section with a spring and then carried out modal analysis for each part of the beam using appropriate matching conditions at the location of the spring. The equivalent stiffness of the spring was calculated using the stress intensity factor at the crack location. Qian et al. [8] derived an element stiffness matrix of a beam with a crack, based on the integration of stress intensity factors. The finite element method was used to study the vibration response of the beam. Abraham and Brandon [9] modelled the opening and closing of a crack using a substructuring approach. Lagrange multipliers and time varying connection matrices were used to represent the interaction forces between the two segments of the cantilever beam separated by the crack. The effect of dry friction when the crack is closed has been accounted for in this model.

A cracked vibrating beam has an inherent non-linearity present in it due to the opening and closing of the crack. Not all of the papers mentioned above model this non-linearity. For example, some authors assume that a crack always remains open. In this paper, we explicitly account for this non-linear behaviour of the beam by not allowing interpenetration of the crack faces. Some of the distinctive features of this paper are: (i) modelling of the opening and closing of the crack as a multi-degree-of-freedom piecewise-linear system; (ii) definition of an effective natural frequency, called the ''bilinear frequency''; and (iii) a systematic study of a two-degree-of-freedom piecewise-linear model to understand the essential non-linear dynamics.

The paper is organized as follows. In section 2, we describe the modal analysis procedure for obtaining the effective natural frequencies for the cracked beam. In section 3, we give the numerical results of the finite element method used to carry out the analysis. In section 4, we begin studying the bilinear two-degree-of-freedom oscillator, and in section 5, we present a brief introduction to non-linear normal modes. In section 6, a perturbation method is used to obtain expressions for the non-linear normal modes of the bilinear system, and in section 7 our results are compared with the results of numerical integration. Our conclusions are presented in section 8.

MODAL ANALYSIS

Consider a uniform cantilever beam, of length L, having a transverse crack at location b of depth a (Figure 1). One can qualitatively see that the global stiffness of the cracked beam would depend on whether the crack is open or closed. In the past (see references [START_REF] S | Vibrations of beams with a fatigue crack[END_REF][START_REF] C | [END_REF]), the cracked beam has been modelled as a bilinear single-degree-of-freedom system (Figure 2). In what follows, we shall briefly describe the procedure for obtaining the ''bilinear frequency'' of this piecewise-linear single-degree-of-freedom system.

--- 

The system shown in Figure 2 consists of a mass m which is attached to a linear spring k 1 and which makes contact with the linear spring k 2 only when x Q 0. The equations of motion for this system are

mx¨+ k 1 x = 0 for x q 0, ( 1
)
mx¨+ (k 1 + k 2 )x = 0 for x Q 0, (2)
where the dots indicate differentiation with respect to time. The free vibration problem can be easily solved analytically (see, for example, reference [START_REF] C | [END_REF]). However, the solution is presented here for the sake of completeness. Before we solve the free vibration problem, let us solve for the two natural frequencies in the two linear regions, x q 0 and x Q 0. The two frequencies are

v 1 = X k 1 m for x q 0, v 2 = X k 1 + k 2 m for x Q 0. (3)
To solve equations ( 1) and ( 2), let the initial conditions be x(0) = x 0 and x˙(0) = 0. Without loss of generality, let us assume that x 0 q 0. The solution to equation ( 1) is

x(t) = x 0 cos v 1 t for t Q t 1 (here x q 0). ( 4
)
At time t 1 = p/2v 1 , x(t 1 ) = 0 and the mass m enters the region x Q 0 with the new initial conditions x(t 1 ) = 0 and x˙(t 1 ) = -x 0 v 1 . The solution to the equation ( 2) is given as

x(t) = - x 0 v 1 v 2 sin v 2 (t -t 1 ) for t 1 E t E t 2 (here x Q 0). ( 5
)
At time t 2 = t 1 + p/v 2 , x(t 2 ) = 0 and the mass m enters the region x q 0 with the initial conditions x(t 2 ) = 0 and x˙(t 2 ) = x 0 v 1 . Note that the velocity with which the mass m enters the region x Q 0 is exactly the same in magnitude as the velocity with which it left the region x Q 0. The solution to equation (1) now is

x(t) = x 0 sin v 1 (t -t 2 ) for t 2 E t E t 3 (here x q 0). ( 6
)
Figure 2. The one-degree-of-freedom system.

At the time t 3 = t 2 + p/2v 1 , x(t 3 ) = x 0 and x˙(t 3 ) = 0; i.e., the mass m, comes back exactly to the same point in the state space (x, x˙) from where the motion started at t = 0. Thus, the motion is periodic and the effective period (T 0 ) of this piecewise-linear system is

T 0 = t 3 = p v 1 + p v 2 = 2p v 0 , (7)
where the effective natural frequency v 0 of this piecewise-linear single-degree-of-freedom system is defined as

v 0 = 2v 1 v 2 v 1 + v 2 . ( 8
)
The effective natural frequency v 0 is called the ''bilinear frequency''. The central idea of this paper is to extend the definition of bilinear frequency from the single-degree-of-freedom bilinear system to the infinite-degree-of-freedom cracked beam problem. Just as the bilinear system (Figure 2) has two linear regions (x q 0 and x Q 0), the cracked beam problem can be modelled as having two configurations: crack open and crack closed. In each of the individual linear regions, one can define a set of natural frequencies by solving an eigenvalue problem. The two sets of natural frequencies can then be combined, using equation (8), to obtain the effective natural frequencies for the cracked beam. In what follows we shall describe the procedure for calculating the frequencies for the two configurations of the cracked beam.

 

The eigenvalues of the cracked beam in the two configurations (crack open and crack closed) can be extracted using a numerical method. Among the various numerical methods present, the most popular are the finite difference method (FDM), the finite element method (FEM) and the boundary element method (BEM). The eigenvalue problem of the cracked beam (Figure 1) is solved using the FEM in this paper.

The basic idea of the finite element method (see reference [START_REF] B | Finite Element Procedures[END_REF]) is to discretize a body into an assemblage of discrete finite elements which are interconnected at the nodal points on element boundaries. The displacement field is approximated over each finite element, in terms of the nodal displacements. The equations of motion can then be formulated using the displacement-based formulation in conjunction with the principle of virtual displacements. These equations of motion can be written as

MU + KU = 0, (9)
where overdots indicate differentiation with respect to time, M and K are the mass and stiffness matrix respectively, and

U (=[u 1 v 1 u 2 v 2 . . . u n v n ] T )
is the column vector of nodal displacements. Each node i has two degrees of freedom, u i and v i , in the horizontal (x) and vertical (y) directions, respectively.

To explain the idea of the two configurations used to calculate the natural frequencies, consider the two crack faces shown in Figure 3. The numbers on the left and right faces of the crack are the node numbers of the interacting nodes. For the first configuration, assume that the crack faces remain open which basically implies that the degrees of freedom associated with the nodes on the crack surfaces are independent. This configuration will be referred to as the ''unconstrained'' configuration. For the second configuration, assume that the crack is closed. Now, to prevent a physically unrealistic situation with interpenetration of crack faces, the horizontal degrees of freedom (u i ) of the nodes on the left and right faces of the crack, which are initially at the same geometric location, need to be restricted. This condition can be stated as

u leftface l = u rightface r , ( 10
)
where node numbers l and r are the corresponding interacting nodes on the left and right face of the crack. This implies that the imposed virtual displacements must also obey the constraint [START_REF] B | Finite Element Procedures[END_REF]. This configuration will be referred to as the ''constrained'' configuration.

The effect of this constraint can be taken into account by an appropriate addition of the corresponding rows and columns of the mass and stiffness matrices. Note that there is no constraint in the vertical direction (v i ) on the crack faces. Therefore, there could be frictionless sliding of the crack faces, and the closed crack configuration need not correspond to an uncracked beam, as assumed in some previous research (see, for example, reference [START_REF] S | Vibrations of beams with a fatigue crack[END_REF]).

Let us assume a time harmonic solution for the nodal displacements:

U = A sin vt, ( 11
)
where A is the amplitude of the nodal displacements. Substituting equation [START_REF] R | [END_REF] into equation ( 9) leads to

[-v 2 M + K]{A} = 0. (12)
For non-trivial solutions,

det (K -v 2 M) = 0, (13)
where ''det'' denotes the determinant and v is the natural frequency. Using equation ( 13) for the two configurations of the cracked beam, two sets of natural frequencies can be obtained. Let us denote the mass and stiffness matrices corresponding to the unconstrained configuration with the subscript 1 and that of the constrained configuration with subscript 2. For the unconstrained configuration, the eigenvalue problem takes the form

det (K 1 -v 2 1i M 1 ) = 0, ( 14
)
where v 1i is the set of natural frequencies corresponding to the unconstrained configuration arranged in order of increasing magnitude. Similarly, for the constrained configuration,

det (K 2 -v 2 2i M 2 ) = 0, ( 15
)
where v 2i is the set of natural frequencies corresponding to the constrained configuration arranged in the order of increasing magnitude.

In this work, we define an effective natural frequency, called the bilinear frequency, which permits us to obtain an approximation for the frequency associated with the ith mode of a piecewise-linear vibrating system without computing the actual mode shape. The bilinear frequency is obtained by computing the associated frequencies of each of the linear pieces of the piecewise-linear system. The important question that we want to address is the relationship between the ith mode shape and natural frequency to those obtained from a ''simplified'' finite element analysis, where one computes the mode shapes and frequencies of the individual linear segments of the piecewise-linear system. We call the latter the ''bilinear mode shapes''. It is expected that if the crack is sufficiently shallow, the two results, namely the actual and bilinear ith mode shapes, will be close to each other. In such a case, we would also expect that the overall frequency of the ith mode will be well approximated by the bilinear frequency as we define it. According to this line of reasoning, calculation of these bilinear mode shapes and frequencies for a vibrating cracked beam appears meaningful. This is done later in this paper. Also, the actual and bilinear mode shapes are computed and compared for a simplified two-degree-of-freedom model of the cracked beam. It is seen that the bilinear mode shapes are comparable to the actual ones. Details are given later in the paper.

The two sets of frequencies (v 1i and v 2i ) for the two configurations of the cracked beam can be combined, using equation ( 8), to define effective natural frequencies (V i ) of the cracked beam as

V i = 2v 1i v 2i v 1i + v 2i . (16)

NUMERICAL RESULTS AND DISCUSSION

Due to the presence of the crack, the beam cannot be modelled using beam elements; but, rather, we need to use plate/shell elements. The cracked cantilever beam has been modelled as a problem in plane stress using eight-noded quadrilateral plate elements. The idea of using beam elements and plate elements in conjunction with transition elements was considered, but since the amount of computational saving was limited the idea was not pursued further. The commercially available finite element package ABAQUS was used to obtain the numerical results. The material and geometric parameters used to carry out the analysis are: E = 210 GPa, r = 7850 kg/m 3 , n = 0•3, L = 10 m and A = 1 m 2 . The location b and the length a of the crack have been varied to study their effect on the natural frequencies of the cracked beam.

First, a convergence study is presented in Tables 1 and2 for a crack location of b = 5 m and a = 0•5 m. Table 1 contains the frequencies for the unconstrained configuration and Table 2 contains the frequencies for the constrained configuration. It can be seen that the effect of increasing the number of elements from 192 to 704 does not have a significant effect on the natural frequencies for the first few modes. All further calculations have been carried out using 704 elements. The undeformed mesh consisting of 704 elements, with b = 5 m and a = 0•5 m, is shown in Figure 4.

Due to the presence of the crack, there is a 1/zr singularity in the stress field at the crack tip. Usually, quarter-point elements are used to model the singular behavior in the stress field. However, in this particular analysis, global changes in the dynamics of the beam are desired, rather than the local variation in the stress field due to the presence of the crack. A comparison in natural frequencies when the singular elements are used and when they are not used is presented in In Figure 7 is shown the variation in the normalized first bilinear effective frequency with the normalized crack depth a/d. The variation shown in Figure 7 can be easily understood if we look at the corresponding mode shapes associated with the first mode (see Figures 5(a In order to justify the intuitive assertion that the bilinear frequency is a good approximation for the effective natural frequency of the multi-degree-of-freedom piecewise-linear system we next investigate a two-degree-of-freedom piecewise-linear system using perturbation methods. It is shown that the bilinear frequency a good approximation to the frequencies of the (non-linear) normal modes in this case.

TWO-DEGREE-OF-FREEDOM SYSTEM

In order better to understand the essential nonlinear dynamics of the cracked beam, we study in this section a simplified piecewise-linear two-degree-of-freedom (2-DOF) system. The system is displayed in Figure 10 and consists of two masses constrained to a straight line and restrained by three springs. Note that spring k 2 applies no force to mass m 1 unless co-ordinate x 1 Q 0. This feature represents, in this simplified model, the essential character of the cracked beam. The equations of motion for this system are (17,[START_REF] L | Regular and Stochastic Motion[END_REF] where x 1 and x 2 are the displacements of the masses m 1 and m 2 respectively, dots indicate differentiation with respect to time, and

m 1 x¨1 + (k 1 + k 3 + k 2 H(x 1 ))x 1 -k 3 x 2 = 0, m 2 x¨2 -k 3 x 1 + k 3 x 2 = 0,
H(x 1 ) = 6 1, x 1 Q 0, 0, x 1 q 0.
As the mass m 1 goes through x 1 = 0, this system switches between each of the two linear problems, one for x 1 q 0 and another for x 1 Q 0. We now obtain the corresponding two sets of eigenvalues and eigenvectors, which will be used in what follows. For x 1 q 0, the characteristic equation is

l 2 -$ k 1 + k 3 m 1 + k 3 m 2% l + k 1 k 3 m 1 m 2 = 0, (19)
where l = v 2 . Equation ( 19) can be solved for the two natural frequencies (v 11 and v 12 ) and then the corresponding eigenvectors can be obtained. These frequencies would actually be the fundamental frequencies of the system in the absence of the spring k 2 .

Carrying out a similar analysis for the case when x 1 Q 0, the characteristic equation takes the form

l 2 -$ k 1 + k 2 + k 3 m 1 + k 3 m 2 % l + (k 1 + k 2 )k 3 m 1 m 2 = 0 (20)
which can be solved for the two natural frequncies (v 21 and v 22 ) and the corresponding mode shapes. The eigenfrequencies in the two regions can be combined in pairs to form two effective bilinear frequencies for this non-linear system. These are given by the same formula used to define the effective bilinear frequency for a one-DOF bilinear system. The bilinear frequencies are defined as

V 1 = 2v 11 v 21 v 11 + v 21 , V 2 = 2v 12 v 22 v 12 + v 22 . ( 21, 22)
In the next two sections we will use perturbation methods and numerical integration to show that the frequencies for the non-linear normal modes associated with the system shown in Figure 10 are well approximated by the bilinear frequencies (21) and (22).

NON-LINEAR NORMAL MODES

Non-linear normal modes are a class of special solutions of n-DOF conservative systems. These could be considered as a non-linear analog of the classical normal modes in a linear system. They were shown by Rosenberg [START_REF] R | [END_REF] to be important because resonance in the forced system typically occurs in the neighborhood of the non-linear normal modes of the free system. Over the years since Rosenberg's original work, there have been many studies of non-linear normal modes. An outstanding recent book by Vakakis et al. [START_REF] V | Normal Modes and Localization in Nonlinear Systems[END_REF] summarizes much of the research in this area to date. We mention here a few works of non-linear normal modes, but the reader is referred to reference [START_REF] V | Normal Modes and Localization in Nonlinear Systems[END_REF] for a complete bibliography.

Rand [START_REF] R | [END_REF] obtained explicit approximate expressions for normal modal curves in non-linear two-degree-of-freedom systems by applying a perturbation method to the modal equation. Rand et al. [14] studied the stability and bifurcation of non-linear normal modes in a class of two-degree-of-freedom systems. King and Vakakis [15] developed an energy-based formulation for computing non-linear normal modes in undamped continuous systems. A partial differential equation is constructed using conservation of energy, which is solved asymptotically using a perturbation methodology. Nayfeh et al. [16] used two approaches to obtain the non-linear normal modes of a cantilever beam. The same non-linear normal mode is obtained by the application of the method of multiple scales to the discretized problem and the original boundary value problem. Shaw and Pierre [17] extended the concept of non-linear normal modes to a multi-degree-of-freedom system with damping. They used the theory of invariant manifolds in conjunction with the center manifold reduction technique.

In the problem at hand, the concept of non-linear normal modes must be applied to a piecewise-linear system. This topic has been studied by Chen and Shaw [START_REF] C | Normal modes for piecewise linear vibratory systems[END_REF], who presented a general discussion of non-linear normal modes in piecewise-linear systems and an example of a 2-DOF system which is similar (but not identical) to the 2-DOF system studied in our paper. They obtained the non-linear normal modes using an invariant manifold approach in the four-dimensional space. In contrast to their approach, we present an analysis of a 2-DOF piecewise-linear system based on a modal equation in the two-dimensional configuration space.

We will define non-linear normal modes as vibrations in unison; i.e., synchronous periodic motions in which all generalized coordinates vibrate equiperiodically, achieving instantaneous zero velocity at the same instant of time. (This is a slight generalization of Rosenberg's definition, suitable for the piecewise-linear system at hand.)

In order to obtain the non-linear modes for the system in Figure 10, we will use Rosenberg's modal equation, derived as follows. Consider a conservative system for which the potential energy V is of the form

V = 1 2 (ax 2 + 2bxy + cy 2 ), ( 23
)
where x and y are the generalized co-ordinates. The constants a, b and c are such that V is positive definite, and kinetic energy T is of the form

T = 1 2 (x˙2 + y˙2). ( 24
)
Then

x¨= -V x , y¨= -V y , T + V = h, (25-27)
where h is a constant equal to the total energy of the system. In what follows, we will eliminate time t by considering y as a function of x only. Using equations (25-27) and y˙= y'x˙, y¨= y0x˙2 + y'x¨, (28, 29)

we obtain the governing equation for the normal modes of the system,

2(h -V)y0 + (1 + y' 2 )(V y -y'V x ) = 0, ( 30
)
where the primes denote differentiation with respect to x. We will define a non-linear normal mode to a be a solution to equation (30) which is single valued and which has two rest points; i.e., two points at which x˙and y˙simultaneously vanish giving T = 0, or, from equation ( 27), V = h. From equation ( 30), at such points where V = h, we find that

V y -y'V x = 0 at V = h. ( 31
)
Geometrically, this implies that a non-linear normal mode intersects the bounding curve V = h orthogonally at two points, as shown in Figure 11. In what follows we shall use a perturbation approach to obtain the nonlinear normal modes of the bilinear two degree of freedom system (Figure 10).

PERTURBATION METHOD

We begin by expressing the equations ( 17) and ( 18) in normal co-ordinates x and y for the x 1 q 0 system. The modal matrix consisting of orthonormal eigenvectors obtained for the region x 1 q 0 is used to achieve the transfer to normal co-ordinates. This transformation can be written in the abbreviated form as

6 x 1 x 27 = $ a 1 a 3 a 2 a 4%6 x y 7 . (32)
The coefficients a i are the components of the eigenvectors and are not listed here for brevity. Note that we will be using the same modal matrix to affect the transformation to the normal co-ordinates for both x 1 q 0 and x 1 Q 0. The equations of motion for the region (x 1 q 0) in the normal co-ordinates then take the form x¨= -v 2 11 x, y¨= -v 2 12 y.

(33, 34)

The potential energy is of the form

V = 1 2 (v 2 11 x 2 + v 2 12 y 2 ) ( 35
)
and the kinetic energy has the same form as equation ( 24). The equations of motion for the region x 1 Q 0 in the normal co-ordinates take the form

x¨= -(v 2 11 + a 2 1 k 2 )x -a 1 a 2 k 2 y, y¨= -a 1 a 2 k 2 x -(v 2 12 + a 2 2 k 2 )y.
(36, 37)

The potential energy is of the form

V = 1 2 (v 2 11 x 2 + v 2 12 y 2 ) + k 2 2 (a 2 1 x 2 + 2a 1 a 2 xy + a 2 2 y 2 ) (38)
and the kinetic energy has the same form as equation ( 24). Note that the presence of an xy term in the equation ( 38) is due to the fact that while co-ordinates x and y are normal co-ordinates for x 1 q 0, they are not normal co-ordinates for x 1 Q 0.

A perturbation method will be used to solve equation (30) approximately. We will assume that the spring k 2 is small:

k 2 = ko, o W 1, ( 39
)
where e is a small parameter. The non-linear normal mode that we desire to solve for can be taken to be a perturbation to the linear mode (y = 0). To solve for a non-linear normal mode y(x), we assume

y(x) = 6 of 1 (x), x 1 q 0, of 2 (x), x 1 Q 0 . ( 40
)
Substituting equation (40) into the modal equation (30), i.e., with V as in equation ( 35), we obtain, neglecting terms of order o, the differential equation governing the non-linear normal mode for the region x 1 q 0:

(2h -v 2 11 x 2 )f0 1 -v 2 11 xf' 1 + v 2 12 f 1 = 0 + O(o). ( 41
)
Similarly, for the region x 1 Q 0 we obtain

(2h -v 2 11 x 2 )f0 2 -v 2 11 xf' 2 + v 2 12 f 2 = -ka 1 a 2 x + O(o). ( 42
)
The solution to the linear ordinary differential equation ( 41) is of the form (see reference [START_REF] R | [END_REF])

f 1 (x) = C 1 sin (D arcsin zpx) + C 2 cos (D arcsin zpx), (43)
where C 1 and C 2 are arbitrary constants,

D = v 11 /v 12 (44)
and

p = v 2 11 /2h. (45)
Equation ( 42) is a non-homogeneous version of equation ( 41). The complementary solution to equation ( 42) is of the same form as equation (41),

f 2,comp (x) = C 3 sin (D arcsin zpx) + C 4 cos (D arcsin zpx), (46)
where C 3 and C 4 are arbitrary constants. A particular solution to equation ( 42) is

f 2,part (x) = ka 1 a 2 v 2 11 -v 2 12
x.

(47)

Combining equations ( 46) and (47), the solution to equation (42) becomes

f 2 (x) = f 2,comp (x) + f 2,part (x). (48)
In order to solve for the constants C 1 , C 2 , C 2 and C 4 , the following conditions will be employed. Comparing equation (41) with equation (30), the boundary condition, equation (31), becomes

f' 1 = v 2 12 f 1 v 2 11 X 1 at X 1 = X 2h v 2 11
q 0 (49) at a rest point (X 1 , Y 1 ) for which X 1 q 0, where Y 1 = of 1 (X 1 ). Similarly,

f' 2 = v 2 12 f 2 v 2 11 X 2 -ka 1 a 2 at X 2 = - X 2h v 2 11 Q 0 (50) at a rest point (X 2 , Y 2 ) with X 2 Q 0, where Y 2 = of 2 (X 2 ). The modal curve runs from (X 1 , Y 1 ) to (X 2 , Y 2 ).
Using the normality conditions (equations (49) and (50)) in equations (43-48), we obtain

C 1 /C 2 = tan (Dp/2) 0 Q for x 1 q 0, (51)
C 3 /C 4 = -tan (Dp/2) = -Q for x 1 Q 0. ( 52
)
Two more conditions are obtained by requiring that the modal curve be continuous and have a continuous first derivative at x 1 = 0. Let (x*, y*) be the point at which x 1 = 0. The matching conditions at the point x 1 = 0 gives that

y* = - a 1 a 2 x* = of 1 (x*), y* = - a 1 a 2 x* = of 2 (x
x* 6 [1 -koa 2 2 /(v 2 11 -v 2 12 )](Q cos z + sin z) -Q sin z + cos z + Q cos z -sin z Q sin z + cos z 7 - koa 2 2 z1 -px* 2 (v 2 11 -v 2 12 )Dzp = 0, (56)
where z = D sin zpx*. Equation ( 56) can be solved using a Taylor series approximation. Since, from equation (56), x* = 0 when o = 0, we write

x* = ox* 1 + O(o 2 ).
(57) Substituting equation (57) into equation (56), and using the computer algebra package MACSYMA to Taylor expand in o, we may collect terms and solve for x* 1 . Once x* is obtained all the constants (C 1 , C 2 , C 3 and C 4 ) can be calculated. This procedure results in an algebraic approximation for the associated non-linear normal mode. An example calculation will be presented in the next section.

To obtain the period associated with the non-linear normal mode, the equations of motion need to be integrated. For the region x 1 q 0, equation (33) can be integrated from x = X 1 to x = x* to give

t 1 = 1 v 11 arcos 0 x* X 1 1 (58)
and, similarly, equation (36) can be integrated from x = X 2 to x = x* to give

t 2 = 1 zv 2 11 + a 2 1 ko arccos 0 x* X 2 1 . ( 59
)
Thus the associated period of the motion can be written as

T p = 2(t 1 + t 2 ). ( 60
)
Note that t 1 and t 2 depend on x*, which is obtained as a Taylor series approximation (57) to equation (56). This period could then be compared with the ''bilinear period'' defined for the 2-DOF piecewise-linear system. In order to obtain an approximation for the other non-linear normal mode in this 2-DOF piecewise-linear system, we may interchange the roles of x and y throughout, and utilize a similar perturbation calculation for x = og(y).

NUMERICAL RESULTS

The results obtained using the perturbation method will be presented for the 2-DOF piecewise-linear system. The values for the parameters in the 2-DOF model (Figure 10 In order to verify the results obtained by the perturbation method, the equations of motion are integrated numerically using the fourth order Runge-Kutta method. For any given initial conditions (x 0 , x˙0 , y 0 , y˙0 ), the equations of motion can be integrated numerically. To obtain the modal curve, we need to choose an approximate starting point on the V = h curve. However, the process of choosing the appropriate initial conditions for numerical simulation is that of trial and error. If any other initial condition is chosen the motion of the phase point is quasi-periodic, as shown in Figure 12. In Figures 13 and15 are shown the non-linear normal modes, obtained using numerical integration, and the bilinear mode shapes for k 2 = 1•0. In Figures 13 and15, the eigenvectors of the individual linear modes (i.e., the modes for x 1 q 0 and x 1 Q 0) are shown as dotted lines. The computation of the bilinear frequency is based upon approximating actual non-linear normal modes by these bilinear shapes. Note that these non-linear normal modes do not pass through the origin x 1 = x 2 = 0. This could be attributed to the fact, that, due to the piecewise-linearity, the bounding curve V = h and associated vector field is unsymmetric about the origin. The corresponding motions of x 1 and x 2 with time are shown in Figures 14 and16.

A comparison of the non-linear normal modes obtained by the perturbation method, and by numerical integration, is presented in Figure 17 for k 2 = 0•1. We see that the resultsagree well. The mismatch between the perturbation method and the numerical results could be attributed to the fact that the results of the perturbation method are based on an order o approximation.

In Table 4 we present a comparison of the period associated with the first non-linear normal mode obtained using the bilinear formula, numerical integration and the perturbation method. It is clearly seen that the bilinear formula provides a good approximation for the period associated with the non-linear normal mode. The perturbation method is seen to do well, as expected, for small values of k 2 .

Figure 19. A special kind of non-linear normal mode corresponding to points C, D and E in Figure 18. V is the potential energy of the system and h is the total energy of the system. Figure 20. The Poincare´map for k2 = 0•0, h = 1•0 and surface of section x1 = 0, x˙1 q 0. This figure corresponds to a linear system; in contrast to Figure 18, which corresponds to a piecewise-linear system: h is the total energy of the system.

Our understanding of the non-linear dynamics of this piecewise-linear system may be further enhanced by examining a numerically generated Poincare´map corresponding to total energy h = 1•0 and a surface of section x 1 = 0, x˙1 q 0 (see reference [START_REF] L | Regular and Stochastic Motion[END_REF]). The non-linear normal modes that we have investigated in the foregoing correspond to the fixed points A and B in the Poincare´map of Figure 18. We note that, in addition to these priodic motions, we expect to see other periodic motions in this piecewise-linear system. As an example, see Figure 18, where C, D and E lie on a 3-cycle. This motion corresponds to a kind of non-linear normal mode; see Figure 19. However, the associated linear system with k 2 = 0 would exhibit no such periodic motion; see Figure 20.

Generalizing this phenomenon from the 2-DOF system to the infinite-dimensional cracked beam, we expect in the case of the latter that there will exist additional non-linear Figure 21. The Poincare´map for k2 = 10•0, h = 1•0 and surface of section x1 = 0, x˙1 q 0: h is the total energy of the system. normal modes to those discussed in this paper (i.e., in addition to those associated with the bilinear frequencies, V i ; equation ( 16)). Their significance is that resonance is expected to occur in their neighborhood, i.e., at forcing frequencies close to their natural frequencies. In addition, we expect that chaotic behavior will occur for large crack depths, in analogy to the behavior of the 2-DOF system for large k 2 ; cf., Figure 21.

CONCLUSIONS

The main goal of this paper is to understand cracked beam dynamics. An eigenvalue problem cannot be formulated for the cracked beam, due to the non-linear boundary conditions associated with the crack being open or closed. However, by modelling the cracked beam as two linear configurations (crack open and crack closed), a definition for the effective natural frequency, called the ''bilinear frequency'', is developed.

The bilinear period is obtained by adding the half-periods of the individual linear modes in the constrained and unconstrained configurations. The reason that this is an approximation rather than an exact result is as follows: if initial conditions were chosen corresponding to the ith linear normal mode of the unconstrained configuration then, when the crack first closes, this would give rise to a motion in the constrained configuration which would in general consist of a superposition of all associated linear modes rather than just the ith constrained mode. However, we have shown that if the initial conditions were appropriately chosen, a periodic motion would result. If the crack depth is not too large, this initial condition would be close to that of the ith unconstrained mode, and the resulting frequency would be close to the ith bilinear frequency. Furthermore, the use of the bilinear frequency as the effective natural frequency was justified through a systematic study of a simplified model of the cracked beam.

A two-degree-of-freedom piecewise-linear system was studied in order to model and understand the essential non-linear dynamics associated with the cracked beam. The non-linear normal modes associated with the model system and also the associated period are predicted using the perturbation method. It is seen that the bilinear formula is indeed a good approximation for the effective natural frequency, especially when the difference between the linear regions is small, i.e., for small values of k 2 . For large values of k 2 , simulation of the 2-DOF system shows that more complicated motions are possible including chaos.

Figure 1 .

 1 Figure 1. A cantilever beam with a transverse edge crack, where a is the length of the crack, b is the location of the crack, L is the length of the beam and d is the depth of the beam.

Figure 3 .

 3 Figure 3. Crack faces.

a

  = 0•5 m. The frequencies are quite close. Therefore, we decided not to use singular quarter-point elements in further calculations.The mode shapes of the beam, with b = 5 m and a = 0•5 m, in the two configurations are shown in Figures5 and 6. In Figures5(a), (b) and (c) are shown the first, second and third modes, respectively, for the unconstrained configuration of the beam, while in Figures6(a), (b) and (c) are shown the first, second and third modes for the constrained configuration of the beam. These mode shapes will be utilized to explain the trends in the effective bilinear frequencies with the variation of parameters such as crack location b and crack length a.

Figure 4 .

 4 Figure 4. The undeformed mesh.

Figure 5 .

 5 Figure 5. Unconstrained mode shapes for the (a) first, (b) second and (c) third mode shapes.

  ) and 6(a)). Since in the first mode of vibration the major contribution to the stiffness comes from the clamped end of the beam, a crack present at b/L = 0•2 results in a substantial decrease in the bilinear frequencies, as compared to b/L = 0•5 or 0•7. The variation in the second and third effective bilinear frequencies with crack depth is shown in Figures8 and 9

  .

Figure 6 .

 6 Figure 6. Constrained mode shapes for the (a) first, (b) second and (c) third mode shapes.

Figure 7 .

 7 Figure 7. The variation of the first bilinear frequency with crack depth, where b is the location of the crack and L is the length of the beam. *, b/L = 0•2; +, 0•5; ×, 0•7.

Figure 8 .

 8 Figure 8. The variation of the second bilinear frequency with crack depth, where b is the location of the crack and L is the length of the beam. *, b/L = 0•2; +, 0•5; ×, 0•7.

Figure 9 .

 9 Figure 9. The variation of the third bilinear frequency with crack depth, where b is the location of the crack and L is the length of the beam. *, b/L = 0•2; +, 0•5; ×, 0•7.

Figure 10 .

 10 Figure 10. The two-degree-of-freedom model.

Figure 11 .

 11 Figure 11. The non-linear normal mode.

Figure 12 .

 12 Figure 12. The motion of a phase point in the x1-x2 plane for a general motion, not a non-linear normal mode.

  ) are chosen as k 1 = k 3 = 1•0 and m 1 = m 2 = 1•0, and the spring constant k 2 is varied. The equation for the non-linear normal mode is computed to be (cf. equation (40)): f 1 (x) = -0•0874032zhk sin (D arcsin zpx)-0•0597945zhk cos (D arcsin zpx), f 2 (x) = 0•0874032zhk sin (D arcsin zpx)-0•0597945zhk cos (D arcsin zpx) -0•2kx, where D = 2•618034 and p = 0•190983/h.

Figure 13 .

 13 Figure 13. The first mode: x1 versus x2 . The dashed curve represents the eigenvectors in each of the linear regions. Note that the bilinear frequency is a better approximation than the mode shape. ---, Non-linear normal mode; • • •, bilinear mode shape.

Figure 14 .

 14 Figure 14. The first mode: ---. x1 ; ---, x2 ; versus time.

Figure 15 .

 15 Figure 15. The second mode: x1 versus x2 . The dashed curve represents the eigenvectors in each of the linear regions. Note that the bilinear frequency is a better approximation than the mode shape. ---, Non-linear normal mode, • • •, bilinear mode shape.

Figure 16 .

 16 Figure 16. The second mode: ---, x1 ; ---, x2 ; versus time.

Figure 17 .

 17 Figure 17. A comparison of the numerical integration and perturbation methods: k2 = 0•1. ---, Numerical integration; ---, perturbation method.

Figure 18 .

 18 Figure 18. The Poincare´map for k2 = 1•0, h = 1•0 and surface of section x1 = 0, x˙1 q 0; h is the total energy of the conservative system.

Table 3 .

 3 The results presented in this table are the natural frequencies for the unconstrained configuration of the beam with b = 2 m and T 1A convergence study for natural frequencies (rad/s) in the unconstrained configuration

			Number of elements
		ZXXXXXXXXXXXXXCXXXXXXXXXXXXXV
		4	20	40	192	704
	First transverse mode	53•40	49•83	49•63	48•33	48•24
	Second transverse mode	385•17	266•51	263•79	247•24	246•40
	First longitudinal mode	789•27	768•75	765•41	746•61	745•88
	Third transverse mode	1165•65	827•83	826•25	825•07	820•98
	Fourth transverse mode	2331•80	1369•1	1361•9	1334•60	1324•40
			T 2
	A convergence study for natural frequencies (rad/s) in the constrained configuration
			Number of elements
		ZXXXXXXXXXXXXXCXXXXXXXXXXXXXV
		4	20	40	192	704
	First transverse mode	54•31	52•09	52•18	52•12	52•03
	Second transverse mode	424•21	309•76	312•34	311•08	310•13
	First longitudinal mode	815•97	806•13	811•84	809•82	809•41
	Third transverse mode	1177•6	827•94	826•38	825•18	821•07
	Fourth transverse mode	2492•70	1496•00	1503•30	1501•60	1487•10
			T 3
	A comparison of natural frequencies (rad/s) for singular versus non-singular elements
			Singular		Non-singular
	First transverse mode		39•62		39•47
	Second transverse mode		311•72		311•68
	First longitudinal mode		688•57		666•96
	Third transverse mode		818•76		818•67
	Fourth transverse mode		1381•10		1380•50

  Solving equation (51) for C 1 , and substituting into equation (53), permits us to obtain the equation for C 2 in terms of x*. Similarly, solving equation (52) for C 3 and substituting into equation (54) gives an equation for C 4 in terms of x*. When these expressions are substituted into equation (55), we obtain a resulting equation for x* as

	*)	(53, 54)
	and
	f' 1 (x*) = f' 2 (x*).	(55)

ACKNOWLEDGMENTS

This research has been supported by contract number ECS-9321508 of the US National Science Foundation with Cornell University. The computing for this research was performed at the Cornell National Supercomputer Facility.


O V Gendelman

Exact solutions for discrete breathers in a forced-damped chain

Exact solutions for symmetric on-site discrete breathers (DBs) are obtained in a forced-damped linear chain with on-site vibro-impact constraints. The damping in the system is caused by inelastic impacts; the forcing functions should satisfy conditions of periodicity and antisymmetry. Global conditions for existence and stability of the DBs are established by a combination of analytic and numeric methods. The DB can lose its stability through either pitchfork, or Neimark-Sacker bifurcations. The pitchfork bifurcation is related to the internal dynamics of each individual oscillator. It is revealed that the coupling can suppress this type of instability. To the contrary, the Neimark-Sacker bifurcation occurs for relatively large values of the coupling, presumably due to closeness of the excitation frequency to a boundary of the propagation zone of the chain. Both bifurcation mechanisms seem to be generic for the considered type of forced-damped lattices. Some unusual phenomena, like nonmonotonous dependence of the stability boundary on the forcing amplitude, are revealed analytically for the initial system and illustrated numerically for small periodic lattices.

I. INTRODUCTION

Discrete breathers (DBs) or intrinsic localized modes (ILMs) are well known in many nonlinear lattices [1,2]. Generally, they are exponentially localized (if a coupling between the neighbors in the lattice has a linear component) and can demonstrate remarkable stability also in two-and three-dimensional lattices [2]. Numerous systems which exhibit the DBs include chains of mechanical oscillators [3], superconducting Josephson junctions [4], nonlinear magnetic metamaterials [5], electrical lattices [6], michromechanical cantilever arrays [7], antiferromagnets [8], and many other physical systems [9].

In the majority of theoretical studies related to the DBs, Hamiltonian models are considered. Still, in many applications the damping cannot be neglected; in order to maintain the DB, one should compensate it by some kind of direct or parametric external forcing [2]. Experimentally, many of the DBs observed in the experiments were in fact created and maintained in the presence of damping and under homogeneous forcing from the external sources.

It is more or less easy to "explain" the nonlinear localization in a chain of forced-damped oscillators on a qualitative level. It is well known that single forced-damped nonlinear oscillators can exhibit stable steady-state responses with different amplitudes, depending on initial conditions [START_REF] Nayfeh | Nonlinear Oscillations[END_REF]. If one or few oscillators in the lattice are excited at high amplitude, and if all the others are at low amplitude, and coupling between the oscillators is weak enough to preserve this structure-one obtains an example of strongly localized excitation in conditions of a homogeneous forcing. However, quantitative description of such breathers remains a major challenge. Lack of Hamiltonian structure changes the properties of the DBs. Instead of a continuous family of localized solutions, one expects to obtain a discrete set of attractors. Consequently, many of the methods devised for computation and analysis of the Hamiltonian DBs are not applicable in forced-damped systems. Even more important, the DBs in the forced-damped systems can have some qualitative properties absent in their Hamiltonian counterparts. For instance, one can observe stable moving DBs in forced-damped systems [START_REF] Marin | [END_REF]. In Hamiltonian lattices, the existence of such moving stable DBs is generally denied [2], with the known exception of the integrable Ablowitz-Ladik model [12]. Also, there exist some special nonintegrable Hamiltonian discrete models where the DBs can be computed exactly [13]. To the best of the author's knowledge, no such exact solutions are known in forced-damped discrete lattices. Even in continuous nonlinear models with forcing and damping, some exact solutions for the breathers were obtained only for a handful of special cases [14]. This paper is devoted exactly to this problem and suggests a model, which allows one to derive exact solutions for the forced-damped DBs and to study some of their properties.

We are going to demonstrate that the exact solutions can be derived for the DBs in a vibro-impact (VI) chain. Dynamical systems involving impacts have an important peculiaritythey exhibit extreme (actually, the strongest possible) nonlinearity, but the latter reveals itself only at the moments of impacts; between the impacts the system obeys linear equations of motion, if other sources of nonlinearity are absent. Consequently, the VI models can offer a relatively simple description of complicated nonlinear phenomena. Celebrated examples of this sort are problems of bouncing ball [START_REF] Tufillaro | An Experimental Approach to Nonlinear Dynamics and Chaos[END_REF] and kicked rotor [START_REF] Casati | Stochastic Behavior in Classical and Quantum Hamiltonian Systems[END_REF], but many other tractable VI models are known and used for description of realistic physical systems, such as Bose-Einstein condensates [17]. One of the few exact solutions for DBs in Hamiltonian systems also uses the VI model [18]; in the current paper, we extend this result for the forced-damped case.

II. DESCRIPTION OF THE MODEL AND ANALYTIC SOLUTION

The model used here is a homogeneous chain with linear nearest-neighbor interactions; in addition, each particle can move only between inelastic impacts constraints. Such model is topologically equivalent to smooth models with linear interaction and on-site potential, widely used in investigations devoted to the DBs [7]. A sketch of the model system is presented in Fig. 1.

Each particle is subject to the same external forcing. The forcing function is considered to be 2π -periodic and antisymmetric, as specified below in the equations of motion. Between the impacts, displacement u n of the nth particle is described as

ün + γ (2u n -u n-1 -u n+1 ) = F (t), F (t) = F (t + 2π), (1)
F (t + π ) =-F (t).
Without affecting the generality, the clearance and the particle mass are set to unity and the only parameter which characterizes the chain is a coupling coefficient γ . Equation ( 1) is valid between successive impacts, i.e., as |u n | < 1 for all particles. If the particle with number n impacts the constraint at time instance t i,n , the following conditions are satisfied:

t = t i,n ; u n (t i,n ) =±1;
(2) un (t i,n + 0) =-k un (t i,n -0), 0 k 1.

Here k is a restitution coefficient. The notation ± 0i su s e d to denote the particle's velocities immediately after ( + ) and before (-) the impact. Values of k<1 correspond to inelastic impacts; it is the only source of damping in the model. The model is not solvable in a general form. Therefore, we are going to look for partial solutions for the DBs. In this paper, the simplest on-site symmetric DB is considered. Our motivation for this choice is that the DBs of this type are the most explored ones, both experimentally and numerically. Hence, we adopt that responses of all particles are periodic and symmetric, and only a particle with n = 0 impacts each of its constraints symmetrically with period 2π :

u n (t) = u n (t + 2πm) =-u n [t + π (2m + 1)] , t i,0 = ϕ + πm, u 0 (ϕ) = 1, (3) u0 (ϕ + 0) =-k u0 (ϕ -0),u 0 (ϕ + π ) =-1, u0 (ϕ + π + 0) =-k u0 (ϕ + π -0).
Here ϕ is a "phase lag" of the impacts with respect to the external forcing. The effect of each impact is equivalent to a transfer of a certain amount of momentum to the impacting particle. Consequently, the solution we look for should satisfy the following modification of Eq. (1) (in terms of distributions):

ün + γ (2u n -u n-1 -u n+1 ) = F (t) + 2pδ n0 ∞ j =-∞ {δ [t -ϕ + π (2j + 1)] -δ(t -ϕ + 2πj )} , (4)
where δ(t) is the Dirac delta function, δ nm is the Kronecker symbol. 2p is the yet unknown momentum transferred to the "central" particle in the course of each impact. For further analysis, the variables are changed as follows:

u n = v n + f (t), f (t) = F (t). ( 5
)
DuetotheantisymmetryofF (t), it is always possible to find a unique function f (t) satisfying Eq. ( 5) and the antisymmetry condition f (t + π ) =-f (t). In terms of the new variables, the external forcing term on the right-hand side of Eq. ( 4) disappears:

vn + γ (2v n -v n-1 -v n+1 ) = 2pδ n0 ∞ j =-∞ {δ [t -ϕ + π (2j + 1)] -δ(t -ϕ + 2πj )}. (6)
In Eq. ( 6) the only remaining forcing term is related to the impacts. The stationary regime is possible only if the impact is symmetric in terms of variables v n (t):

v0 (ϕ -0) =-v 0 (ϕ + 0) = p.
With the account of this relationship, Eq. (2) for the inelastic impact is rewritten as

v 0 (ϕ) + f (ϕ) = 1; u0 (ϕ + 0) =-p + ḟ (ϕ) =-k[p + ḟ (ϕ)] ⇒ ḟ (ϕ) = p 1 -k 1 + k . (7)
Equation ( 7) has three unknowns, p, ϕ, and v 0 (ϕ). To define the solution completely, an additional relationship between p, ϕ, and the value of v 0 in time instance of the impact v 0 (ϕ)is required. One obtains it by rewriting Eq. ( 6) with the help of a formal Fourier series:

vn + γ (2v n -v n-1 -v n+1 ) =- 4pδ n0 π ∞ l=0 cos [(2l + 1)(t -ϕ)]. (8)
Equation ( 8) represents a common linear chain with an external forcing applied to the particle n = 0. The steady-state solution will be spatially localized if all frequencies in the forcing term lie in an attenuation zone of the chain. The latter requirement will be satisfied if γ ∈ [0,0.25). In this case the steady-state solution of Eq. ( 8) may be presented in the following form [17]:

v n (t) = 4p π -1 2γ |n| ∞ l=0 [(2l + 1) 2 -2γ -(2l + 1) 4 -4γ (2l + 1) 2 ] |n| (2l + 1) 4 -4γ (2l + 1) 2 cos [(2l + 1)(t -ϕ)]. (9)
For the sake of completeness, the derivation details of Eq. ( 9) are presented in Appendix A.

Fourier series ( 9) converges for all n and for γ ∈ [0,0.25), i.e., for all values of coupling, for which the forcing frequency stays in the attenuation zone. An explicit expression for v 0 (ϕ) is written as

v 0 (ϕ) = pχ(γ ),χ (γ ) = 4 π ∞ l=0 1 (2l + 1) 4 -4γ (2l + 1) 2 . ( 10
)
Taking into account Eq. ( 10), the system (7) can be closed:

pχ(γ ) + f (ϕ) = 1, ḟ (ϕ) = p 1 -k 1 + k . (11)
Then, the exact solution for the symmetric one-site DB is expressed as

u n (t) = f (t) + 4p π -1 2γ |n| ∞ l=0 [(2l + 1) 2 -2γ -(2l + 1) 4 -4γ (2l + 1) 2 ] |n| cos [(2l + 1)(t -ϕ)] (2l + 1) 4 -4γ (2l + 1) 2 . ( 12
)
Parameters p and ϕ can be determined from Eq. ( 11) for each specific choice of F (t).

In order to study the DB described by Eqs. ( 11) and ( 12) in more depth, the case of a simple harmonic forcing is considered. In this case, one can solve Eq. ( 11) explicitly:

F (t) =-a cos t ⇒ f (t) = a cos t; -a cos ϕ = pχ(γ ) -1, -a sin ϕ = pq, q = 1 -k 1 + k , (13)
p = χ (γ ) ± χ 2 (γ )a 2 -q 2 (1 -a 2 ) χ 2 (γ ) + q 2 ,ϕ =-arcsin pq a .
Here a is the amplitude of the harmonic forcing. A stable response corresponds to a positive sign in expression (13) for p.

From ( 13) one obtains the condition for minimum forcing amplitude necessary to support the DB solution, as a function of the restitution coefficient and the coupling:

a q χ 2 (γ ) + q 2 . ( 14
)
To ensure self-consistency of the solution starting from Eq. ( 4), two conditions should be satisfied:

(1) The particle with n = 0 should impact the barriers only at the time instances t = ϕ, ϕ + π ;

(2) The particles with n = 0 should never reach the impact constraints. Condition (1) is satisfied for solution (12)-( 13) for not-too-small values of the restitution coefficient. It is demonstrated in Appendix B that for the values of k 0.1 the solution is self-consistent for all possible values of the other parameters. To satisfy Condition (2), it is enough to take care about this condition for particles with n =±1, since the amplitudes of all particles with |n| > 1 would be even smaller (a proof of this fact and more detailed exploration of the solution self-consistency are also presented in Appendix B). Then, for any time instance, one should obtain

|u 1 (t)| = a cos t - 2p γπ ∞ l=0 (2l + 1) 2 -2γ -(2l + 1) 4 -4γ (2l + 1) 2 (2l + 1) 4 -4γ (2l + 1) 2 cos [(2l + 1)(t -ϕ)] < 1. (15)
Conditions ( 14) and (15) determine the zone of existence of the DB in the space of parameters. This zone is illustrated in Fig. 4, together with the stability thresholds for the DBs.

III. STABILITY ANALYSIS

The DB solution ( 12)-( 13) is 2π -periodic by construction, and thus its stability may be established by analysis of eigenvalues of a monodromy matrix computed on a single time period [START_REF] Landau | Mechanics[END_REF]. Such matrix cannot be computed for an infinite chain; therefore, it is commonly accepted to compute it for a finite chain and to check whether the stability properties depend on the number of particles [START_REF] Marin | [END_REF]. In this paper we adopt the same approach and consider the system with N particles and periodic boundary conditions. For systems with smooth nonlinearity, such monodromy matrices are usually computed numerically, by solving a complete set of N ordinary differential equations (ODEs) for 2N independent sets of initial conditions [START_REF] Marin | [END_REF][START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF]. However, for the vibro-impact chain considered here, in addition to availability of exact solutions, one can compute the monodromy matrix in a much simpler manner. A state vector of the system is defined as We adopt w k = uk ,k= 0,...,N. The periodic boundary conditions are imposed as

V = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ u 0 . . . u N-1 w 0 . . . w N-1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (16)
u 0 = u N ,w 0 = w N . (17)
In the exact solution ( 12)-( 13) the impacts occur only at the site n = 0 in time instances t = ϕ, ϕ + π . Consequently, between impacts, the time evolution of the state vector is described by a linear equation. With account of (1) and ( 16), this equation is written as

V = AV + F, A = 0 N×N I N×N γL N×N 0 N×N , F = F (t) 0 N×1 I N×1 ϕ<t<φ+ π, ϕ + π<t<φ+ 2π. ( 18
)
Here

L N×N = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -21 0•••
10 ••• 01 -2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
is the discrete Laplacian adjacency matrix accounting for the linear coupling.

Due to the linearity of Eq. ( 18), evolution of small perturbations of the state vector V between the impacts is described by the exponent of matrix A. At the moments of impact, the evolution of small perturbations of the state vector is described by the so-called saltation matrix [START_REF] Fredriksson | [END_REF], which takes into account both the perturbation of the state vector at the point of discontinuity and the infinitesimal variation of the impact time instance. For the system under consideration, this matrix at time instance t = ϕ in the lowest-order (linearized) approximation [START_REF] Fredriksson | [END_REF] may be written as

S = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -k 0 ••• 0 01 . . . . . . . . . . . . . . . 0 0 ••• 01 0 N×N ω(1+k) w 0 (ϕ-0) 0 ••• 0 00 . . . . . . . . . . . . . . . . . . 0 ••• ••• 0 -k 0 ••• 0 01 . . . . . . . . . . . . . . . 0 0 ••• 01 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (19)
where ω = ü0 (ϕ -0). The complete monodromy matrix at the period 2π is thus expressed as

M = (QS) 2 , Q = exp(π A). (20)
Values of velocity and acceleration of the impacting particle in the moment of impact are taken from the exact solution ( 12)- (13). Expression [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF] still requires the numeric computation of the exponent matrix and the multiplication of the matrices, but all of these tasks are much easier from a computational viewpoint than the solution of N second-order ODEs with 2N independent sets of the initial conditions. The loss of stability has been recorded numerically in the space of parameters, when the eigenvalue of the monodromy matrix (or couple of the complex conjugate eigenvalues) with the largest absolute value crossed a threshold of 1.005. For all explored examples, it was found that for N 400 the critical parameter values do not depend on N within a relative margin of about 0.1%. Thus, it is possible to suggest that these results on stability of the DBs are also relevant for the infinite chain. Two generic scenarios of the loss of stability were revealed. The first one corresponds to the transition of two complex conjugate eigenvalues through the unit circle and is related to the Neimark-Sacker bifurcation. The other scenario corresponds to the passage of a single eigenvalue through unity. This scenario corresponds to a pitchfork bifurcation; the latter results in appearance of a pair of stable asymmetric DBs. Typical structure of the eigenvalues of the monodromy matrix for both of these bifurcation scenarios is presented in Figs. 2 and3.

Therefore, for instance, it is possible to infer from Fig. 2 that for the set of parameters corresponding to case (a) the symmetric DB will be stable, and for case (b) it will be unstable.

Examples for zones of existence and stability for the DBs are presented in Fig. 4. The lower boundary of the zones of existence is described by Eq. ( 14); the upper one is described by Eq. [START_REF] Tufillaro | An Experimental Approach to Nonlinear Dynamics and Chaos[END_REF].

One can see that for lower values of the restitution coefficient, the DB exists for a narrower range of amplitudes of the external forcing, but it is stable in a wider range of the values of coupling coefficient. The pitchfork bifurcation in the upper left corner is similar to the loss of symmetry observed in a single vibro-impact oscillator [START_REF] Shaw | [END_REF]; naturally, this case is equivalent to zero coupling. A similar pitchfork bifurcation, related to a "loss of symmetry," is also known in the regular forced-damped Duffing oscillator [START_REF] Nayfeh | Nonlinear Oscillations[END_REF]. Thus, this scenario of the loss of stability is caused by the internal dynamics of each individual oscillator. Our computations reveal an interesting and, it seems, previously unknown fact: The coupling between the oscillators can suppress this scenario of the loss of stability and stabilize the symmetric DBs.

It is also possible to conjecture that the Neimark-Sacker bifurcation in the lower right corner is related to four-wave interaction with the boundary of the propagation zone [7]. One can also see that the size of the interval of existence of the DBs shrinks and approaches zero as the frequency approaches the boundary of the propagation zone. This observation also correlates with numeric and experimental findings presented in Ref. [7], where the amplitude of the DB was seen to decrease as the forcing frequency gets closer to an optical band. These results allow us to conjecture that both bifurcation mechanisms described above are generic for forced lattices with nonlinear on-site potentials, and not only relevant for the VI chain studied here.

In Fig. 4(a) one can see that for relatively low damping (k = 0.8) the stability threshold on the parameter plane is not monotonic. It is not clear whether this property is generic. In the next section, this finding is illustrated by direct numeric simulations.

IV. NUMERIC SIMULATIONS

Usually, experiments with the DBs are performed with a relatively small number of particles or nodes. Therefore, in order to illustrate the analytic findings, numeric simulations are performed for a system with N = 20 particles and periodic boundary conditions. Additionally, the inelastic impacts were replaced by strongly nonlinear smooth force and damping, which exactly correspond to the Newtonian model of the inelastic impact in the limit case specified below [23]. Also, small on-site linear damping with coefficient μ = 0.005 was included. Thus, the following smooth system has been simulated:

ün + γ (2u n -u n-1 -u n+1 ) + μ un + λ(2m + 1)u 2m n un + (2m + 1)u 4m+1 n = a cos t, u 0 = u N . ( 21
)
Here m is a natural number. In Ref. [23] it was proved that the two last terms on the left-hand side of Eq. ( 21) in the limit m →∞ faithfully model the inelastic impact with a velocity-independent restitution coefficient. The value of this coefficient is related to parameter λ as follows [23]:

k = exp -πλ 2 1 -λ 2 /4 . ( 22
)
The system has been simulated numerically by the Rosenbrock method suitable for stiff systems. It was observed that for values of m>5 the simulation results ceased to depend on the exact choice of m.

Convergence to the exact DB solution has been simulated for a = 0.1, k = 0.8, and γ = 0.08. The system has been initialized by initial velocity u0 (0) = 3; remaining initial conditions were set to zero. The simulated response of the system after an initial transient is compared to the exact solution ( 12)-( 13) for the vibro-impact counterpart in Figs. 5(a)-5(c).

One can conclude that the DB for the given set of parameters is a genuine attractor of the dynamics and that the exact solution ( 12)-( 13) is robust in a sense that it nicely approximates similar excitation in a smooth system with a relatively small number of degrees of freedom.

In order to illustrate the nonmonotonic dependence of stability on the coupling revealed in Fig. 4(a), system [START_REF] Shaw | [END_REF] was simulated for initial conditions exactly corresponding to solution ( 12)-( 13), the restitution coefficient k = 0.8, coupling γ = 0.104, and three different values of the external forcing a = 0.1, 0.2, 0.7. From Fig. 4(a) it is clear that all these values of the forcing and coupling are rather close to the boundary of stability. The results of this simulation are presented in Fig. 6.

For cases (a) and (c), the system settles in the attractor corresponding to the symmetric DB after some initial transient. For case (b), the system preserves the localization for some time, but then the response becomes delocalized and the DB structure is eventually lost. The results of simulations completely correspond to predictions in Fig. 4(a), including the nonmonotonous behavior of the stability threshold.

One can mention that in the simulation presented in Fig. 6(b) the DB is destroyed. It does not turn out to be quasiperiodic, as one would expect from the generic Neimark-Sacker bifurcation scenario described above. This result suggests more careful study of the DB structure transformation when the stability boundary is crossed. For this sake, the symmetric DB was created according to the exact solution in the region of stability. Then, the amplitude of the forcing was varied (at the rate of about 10 -3 with respect to the forcing frequency), in a way that the system has crossed the stability threshold. It was revealed that the localized DBs in such conditions could persist for a rather large variation of the forcing amplitude; however, their structure undergoes significant changes. The simulations were performed for both stability boundaries (corresponding to the pitchfork and the Neimark-Sacker bifurcations). The time series for the final DB shape are presented in Figs. 7 and8, respectively.

The time series obtained in both simulations correspond to the DBs, but the latter are no more symmetric. Analytic solutions for the unstable symmetric DBs are depicted in these figures by solid lines for the sake of comparison. In Fig. 7 one can observe that the numeric solution is still periodic, but asymmetric. The solution presented in Fig. 8 is no more 2π-periodic, but still spatially localized. One can suppose that this solution might be either quasiperiodic or phase locked with some large period. Both results presented in Figs. requires many additional analytic and numeric efforts and will be presented elsewhere.

V. CONCLUDING REMARKS

As was demonstrated above, the vibro-impact model allows derivation of exact analytic solutions for discrete breathers in forced-damped essentially nonlinear chain, without any simplifications or approximations. The zone of existence and stability of the DB solution in the space of parameters has been computed. Two mechanisms of the loss of stability (pitchfork and Neimark-Sacker bifurcations) were revealed. Existence of the asymmetric and quasiperiodic DBs was predicted analytically and verified numerically. It is interesting to mention that the increase of the coupling can both stabilize and destabilize the symmetric DB in different zones on the plane of parameters and for different bifurcation mechanisms. Also, nonmonotonic behavior of the stability boundary in the space of parameters was revealed analytically for the infinite chain and illustrated numerically for small smooth systems with periodic boundary conditions.

One should note that the numeric simulations were performed for the model with smooth, rather that vibroimpact, on-site potentials. The results were very similar to the theoretical counterparts obtained for the vibro-impact model. This coincidence allows one to suggest that the results concerning the behavior of the on-site DBs (generic bifurcations, "topography" of the existence and stability zones on the plane of parameters) are in fact relevant for a wider class of systems with topologically similar potentials of interaction. This suggestion will be explored in a future work.

= p + 16γp π sin(t -ϕ) √ 1 -4γ [1 + √ 1 -4γ ] + ∞ l=1 sin [(2l + 1)(t -ϕ)] (2l + 1)[2l + 1 + (2l + 1) 2 -4γ ] (2l + 1) 2 -4γ p + 16γp sin(t -ϕ) π 1 √ 1 -4γ (1 + √ 1 -4γ ) - ∞ l=1 1 [2l + 1 + (2l + 1) 2 -4γ ] (2l + 1) 2 -4γ p + 16γp sin(t -ϕ) π 1 2 - 1 8 ∞ l=1 1 l 2 = p + 16γp sin(t -ϕ) π 1 2 - π 2 48 . ( B5
)
The first inequality in (B5) explores the fact that sin [(m + 1)x] (m + 1) sin x for all non-negative whole numbers m and for x ∈ [0,π ]. This fact easily follows from recurrence relationships for the second-order Chebyshev polynomials [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]. In the last transformation we took into account that γ ∈ [0,0.25). Using estimation (B5), one obtains from Eq. (B4)

V 0 (t) =-a sin t - 4p π ∞ l=0 sin [(2l + 1)(t -ϕ)] (2l + 1) 2 -4γ -a sin t -p - 16γp sin(t -ϕ) π 1 2 - π 2 48 -a sin t -p. (B6)
The last transformation holds since sin(tϕ) is non-negative for ϕ t ϕ + π .I tf o l l o w sf r o m(B6) that u 0 (t) will monotonously decrease in the interval ϕ<t<ϕ+ π for all values of γ ∈ [0,0.25) if it monotonously decreases for γ = 0. It is clear from (B6) that the latter condition will be fulfilled if p>ain the interval 1 a q √ χ 2 (0)+q 2 ; from expression (13), it is easy to derive that the inequality holds for all k 0.1. So, one can conclude that Statement 1 is correct for not-too-small values of the restitution coefficient. It should be mentioned that for very small values of k Statement 1 really turns out to be wrong-that happens, for instance, for the set of parameters k = 0.005,a= 0.99,γ= 0. Therefore, although it might be possible to tighten the estimations used in inequalities (B5) and (B6), still the consistency of solution ( 12)-( 13) will be violated for small values of the restitution coefficient k.

The proof of Statement 2 essentially relies on the fact that the harmonic component with unit frequency is dominant in all sums u n (t)forn = 0. Solution ( 12)-( 13) may be transformed as follows: Note that we consider only n>0; the opposite case is similar due to the symmetry of the exact solution. Using the former notations, Eq. ( 15) can be transformed in the following way:

u n (t
max t [u 1 (t)] 1 ⇔ max t [A 1 (t) + R 1 (t)] 1 ⇒ max t [A 1 (t)] 1 + max t [R 1 (t)] . (B9)
Further, the following notations will be used:

Ân = max t A n (t), Rn = max t R n (t).
As for Rn , the following estimation may be obtained:

Rn = max

FIG. 1 .

 1 FIG.1. Sketch of the model system.

FIG. 2 .

 2 FIG. 2. (Color online) Zones of existence and stability for the symmetric DB on the (γ , a) parameter plane for two different values of the restitution coefficient: (a) k = 0.8; (b) k = 0.4. Crosses denote the pitchfork bifurcation line and diamonds denote the line of Neimark-Sacker bifurcation. Transitions corresponding to illustrations in Figs. 2 and 3 are marked by arrows in (a). "S" and "US" denote zones of stability and instability, respectively.

FIG. 3 .

 3 FIG. 3. Evolution of eigenvalues of the monodromy matrix corresponding to Neimark-Sacker bifurcation. N = 400, a = 0.2, k = 0.8; (a) γ = 0.09, (b) γ = 0.11.

FIG. 4 .

 4 FIG. 4. Evolution of eigenvalues of the monodromy matrix corresponding to pitchfork bifurcation. N = 400, a = 0.8, k = 0.8; (a) γ = 0.05, (b) γ = 0.045.

FIG. 5 .

 5 FIG. 5. (Color online) Comparison between exact solution for infinite chain (12)-(13) and numeric simulation of system (22) for N = 20; periodic boundary conditions, a = 0.1, k = 0.8, γ = 0.08, and μ = 0.005. The plots present almost indistinguishable exact solution (red, solid line) and numeric result (black, circles) for (a) u 0 (t); (b) u 1 (t); (c) u 10 (t).

FIG. 6 .

 6 FIG. 6. (Color online) Numeric simulation of the DBs in small system with smoothened potentials and damping, N = 20, k = 0.8, γ = 0.104; (a) a = 0.1, (b) a = 0.2, (c) a = 0.7. Time series for five central particles (-2 n 2) are presented.

FIG. 7 .

 7 FIG. 7. (Color online) Comparison between exact solution for the symmetric DB in the infinite chain (12)-(13) (red, solid line) and numeric simulation of system (22) (black, circles) for N = 20; periodic boundary conditions, a = 0.82, k = 0.8, γ = 0.04, and μ = 0.005. The time series are depicted for (a) u 0 (t); (b) u 1 (t); (c) u 10 (t).

  7 and 8 completely conform to the predictions of the stability analysis presented above. A study of precise structure and stability properties of the asymmetric and quasiperiodic DBs

FIG. 8 .

 8 FIG. 8. (Color online) Comparison between exact solution for the symmetric DB in the infinite chain (12)-(13) (red, solid line) and numeric simulation of system (22) (black, circles) for N = 20; periodic boundary conditions, a = 0.15, k = 0.8, γ = 0.104, and μ = 0.005. The time series are depicted for (a) u 0 (t); (b) u 1 (t); (c) u 10 (t).

t[+ 2 )

 2 R n (t)] = max t (-1) n 4p π ∞ l=1 ξ n l cos [(2l + 1)(tϕ)] (2l + 1) 4 -4γ (2l + 1) = pγ n μ n ,μ n = ζ (2n + 2) 4 n π . (B10)

  In this expression A n (t) is the component with the unit frequency, and R n (t) denotes the sum of all other harmonics,

	) = a cos t +	4p π	2γ -1	|n| ∞
	A n (t) = a cos t + R n (t) = (-1) n 4p π	(-1) n 4pξ n 0 π √ 1 -4γ ∞ l=1 ξ n l cos [(2l + 1)(t -ϕ)] cos(t -ϕ). (2l + 1) 4 -4γ (2l + 1) 2	.	(B7)
	ξ l =	(2l + 1) 2 -2γ -(2l + 1) 4 -4γ (2l + 1) 2 2γ
	=	2γ (2l + 1) 2 -2γ + (2l + 1) 4 -4γ (2l + 1) 2	=	2l + 1 -(2l + 1) 2 -4γ 2l + 1 + (2l + 1) 2 -4γ	.	(B8)

l=0 ((2l + 1) 2 -2γ -(2l + 1) 4 -4γ (2l + 1) 2 ) |n| cos [(2l + 1)(tϕ)] (2l + 1) 4 -4γ (2l + 1) 2 = A n (t) + R n (t).

ACKNOWLEDGMENTS

The author is very grateful to Professor Yuli Starosvetsky for useful discussions and suggestions and to the anonymous reviewer for great help in improving the paper.

APPENDIX A

In this Appendix we present the detailed derivation of Eq. (9). It is natural to look for the steady-state solution of Eq. (8) in the form of the Fourier series:

Substituting this expression into Eq. ( 8), one obtains -(2l + 1) 2 v n,l + γ (2v n,l -v n-1,l -v n+1,l ) =-4pδ n0 π .

(A2)

For n = 0Eq.(A2) is homogeneous. If one requires vanishing of the solution as n →±∞, then the appropriate solution of (A2) is expressed as

For n = 0 one substitutes Eq. (A3) into Eq. (A2) and obtains the following expression for v 0,l :

Substituting Eqs. (A3) and (A4) into Eq. (A1), one obtains Eq. (9).

APPENDIX B

In this Appendix the proof of the self-consistency of analytic solution (12)-( 13) for not-too-small values of the restitution coefficient k is presented for the case of simple harmonic forcing. The proof is based on two statements, which will be analyzed separately:

(1) The particle with n = 0 undergoes impacts with the constraints at time instances t = ϕ, ϕ + π and does not impact the constraints at any other time instance;

(2) The upper boundary for existence of the single-site DB is limited by the curve corresponding to the condition

for all values of |n| > 1.

If Statements 1 and 2 are correct, they guarantee that no particle besides the one with n = 0 is engaged in the impacts with the constraints, and the central particle meets the constraints only at t = ϕ, ϕ + π . Thus, the transition from Eqs. ( 1) and (2) to Eq. ( 4) and all analytic treatment that follows would be justified.

Let us first consider Statement 1. It is equivalent to the statement that

Due to the symmetry of the considered solution, it is enough to prove (B3) for ϕ<t<ϕ+ π . If one recalls that u 0 (ϕ) = 1,u 0 (ϕ + π ) =-1, it is clear that it is sufficient to demonstrate a monotonous decrease of u 0 (t)f o rϕ<t< ϕ + π .

From Eq. ( 12), velocity of the particle with n = 0 is written as

Due to a discontinuity of the velocity at the impacts, series (B4) does not converge uniformly. In order to obtain suitable estimations of the velocity, this series should be presented in somewhat modified shape:

To derive estimation (B11), we used expression (B8) and the fact that γ ∈ [0,0.25). ζ (x) denotes the Riemann zeta function.

For values of Ân one can obtain exact expressions from definition (B7) with the help of conditions (13) as follows:

From inequality (B9) it follows that

If one considers inequality (B2) in a similar manner, it is possible to derive the following estimation:

In order to prove the latter inequality [and, consequently, inequality (B2)] it is enough to demonstrate that 2(r 1 + γμ 1 )

Then, (B13) would follow from (B12). Inequality (B13) can be rewritten as

From (B11) one can derive the following estimations:

Using these estimations, one can easily demonstrate the correctness of inequality (B15) for all γ ∈ [0,0.25),q ∈ [0,1]. It should be mentioned that the equality in (B10) and (B15) holds only for γ = 0 and in (B13), and, consequently, in (B2)-only for a degenerate case with a = 1,γ= 0.