Eran Grosu

Isaac Harari
email: harari@eng.tau.ac.il

Stability of semidiscrete formulations for elastodynamics at small time steps

Keywords: Semidiscrete, Small time step oscillation, Rothe method, Spatial stabilization, Principle of causality

Solutions of direct time-integration schemes for elastodynamics that converge in time to conventional semidiscrete formulations may be polluted at small time steps by spurious oscillations that violate the principle of causality, for example by arising before wave fronts. This degradation is not an artifact of the time-marching scheme, but rather a property of the solution of the semidiscrete formulation itself. An analogy to singularly perturbed elliptic problems provides an upper bound on the time step for the onset of these oscillations. A simple procedure of spatial stabilization is proposed to remove this pathology from implicit time-integration schemes, without affecting unconditional temporal stability. Spatially stabilized implicit time-integration methods are free of noncausal oscillations at small time steps.

Introduction

Dynamic behavior of elastic bodies (elastodynamics) is described by the equation of motion, which is continuous in space and time. Finite element methods are usually used to obtain numerical solutions for problems in structural dynamics.

The common approach to transient computation (the method of lines [START_REF] Berezin | Computing Methods, vols. I and II[END_REF]) is time-integration of a semidiscrete formulation, which is obtained by spatial discretization. Thus, the approximation is carried out in two stages. First, spatial approximation, e.g., by standard finite element methods, leads to the semidiscrete formulation (a coupled system of ordinary differential equations in time). Then, temporal discretization by time-integration schemes results in a system of algebraic equations at each time level. Perhaps the most widely used family of direct methods for solving the semidiscrete equation of motion is the Newmark family of methods for time-integration. The unconditionally stable in time, second-order accurate trapezoidal rule is, probably, the most commonly used implicit algorithm among the Newmark family.

It is well known that higher modes, related to higher frequencies of semidiscrete formulations, are approximated poorly [START_REF] Strang | An Analysis of the Finite Element Method[END_REF]Section 6.3]. As a result, algorithmic damping is often introduced in time-integration schemes in order to remove participation of high-frequency modal components [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF][START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: the generalizedmethod[END_REF]. Algorithmic damping cannot be introduced in the Newmark method without degrading the order of accuracy. Hence, the trapezoidal rule lacks algorithmic damping. The implicit, unconditionally stable in time, second-order accurate HHT- [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF] and generalized- [START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: the generalizedmethod[END_REF] algorithms are extensions of the Newmark method with adjustable algorithmic damping introduced to improve this shortcoming.

Conventional wisdom advocates that time step reduction be accompanied by corresponding refinement of the spatial mesh. Indeed, procedures for time step selection often advise against reducing time steps far below the critical values for temporal stability [5, p. 510]. Nevertheless, small time steps are often necessary in practice, for example in problems of radiative transport [START_REF] Ober | Studies on the accuracy of time-integration methods for the radiation-diffusion equations[END_REF] and fluid-structure interaction [START_REF] Masud | Effects of mesh motion on the stability and convergence of ALE based formulations for moving boundary flows[END_REF]. As the time step is reduced with a fixed mesh size, the deleterious effects of higher modes are inevitably admitted into the computation, even with algorithmic damping [START_REF] Krenk | State-space time-integration with energy control and fourthorder accuracy for linear dynamic systems[END_REF]. Thus, spurious spatial oscillations may pollute the solution, at small time steps, of all algorithms that converge in time to consistent semidiscrete formulations.

The principle of causality in classical mechanics states that cause precedes its effect. With relation to computation of discontinuous elastodynamic phenomena, numerical solutions may be polluted by spurious oscillations resulting from the sharp gradient at the front of a stress wave. Spurious oscillations arising before the wave front are noncausal. Specifically, computation of discontinuous phenomena with the trapezoidal rule is accompanied by spurious oscillations, due to its lack of algorithmic damping. At relatively large time steps these oscillations follow the wave front, but for small time steps the oscillations precede the wave front, violating the principle of causality. The HHT- [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF] and generalized- [START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: the generalizedmethod[END_REF] algorithms significantly reduce the spurious oscillations observed in the trapezoidal rule. Nonetheless spurious oscillations precede the wave front at small time steps.

An alternative approach to the method of lines is the Rothe method [START_REF] Rothe | Zweidimensionale parabolische Randwertaufgaben als Grenzfall eindimensionaler Randwertaufgaben[END_REF] (or horizontal method of lines) of first discretizing in time and then in space on each discrete time level. The Rothe method is employed in order to characterize the onset of noncausal oscillations and eliminate them. This method reveals that the time-discrete equation, that is solved at each time level of standard implicit schemes, can be viewed as a Galerkin approximation of a steady reaction-diffusion equation governed by a modified Helmholtz operator. An analogy to such singularly perturbed elliptic problem provides a threshold for the onset of small time step oscillations. Numerical computations obtained for smaller time steps might produce noncausal oscillations.

Any scheme that stabilizes the singularly perturbed equation removes this pathology. We propose a simple procedure of spatial stabilization, which leads to implicit time-integration schemes that are free of noncausal oscillations, without effecting unconditional temporal stability. A stabilized method for the Rothe representation is obtained by appending to its Galerkin equation residuals of its Euler-Lagrange equation in adjoint least-squares form. These residuals are multiplied by a meshdependent stability parameter, which is defined by von Neumann analysis. The stabilized method of adjoint type, originally called the unusual stabilized finite element method [START_REF] Franca | Stabilized finite element methods. I. Application to the advective-diffusive model[END_REF] and also subgrid modeling [START_REF] Cipolla | Subgrid modeling in a Galerkin method for the Helmholtz equation[END_REF] or multiscale stabilization [START_REF] Hughes | Conservation properties for the Gelrekin and stabilized forms of the advection-diffusion and incompressible Navier-Stokes equations[END_REF], may be derived in the variational multiscale framework [START_REF] Hughes | Multiscale phenomena: Green's functions, the Dirichletto-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods[END_REF], and is related to residual-free bubbles [START_REF] Brezzi | A relationship between stabilized finite element methods and the Galerkin method with bubble functions[END_REF][START_REF] Franca | Bubble functions prompt unusual stabilized finite element methods[END_REF]. In practice, for the self-adjoint operators considered herein, this method is equivalent to Galerkin/least-squares, which is often employed for such problems [START_REF] Harari | Improved finite element methods for elastic waves[END_REF].

As mentioned previously, it is well known that effects of higher modes are inevitably admitted into the computation as the time step is reduced with a fixed mesh size. The objective of this paper is to characterize the onset of instability and propose a stabilization procedure for eliminating this pathology from implicit time-integration schemes for elastodynamics. Similar work was done for the diffusion equation, as a model for parabolic problems [START_REF] Harari | Stability of semidiscrete formulations for parabolic problems at small time steps[END_REF]. In this paper, the equation of motion is considered as a model for hyperbolic problems.

The remainder of this paper is organized as follows. In Section 2 we consider the Newmark time-integration method for the semidiscrete equation of motion, and demonstrate the existence of noncausal oscillations at small time steps. Improved methods, such as the generalized-method, concentrate on enhancement of the algorithmic damping capabilities, which can only reduce the nonphysical behavior rather than eliminate it. An alternative approach to the derivation of the time-discrete equations shows, in Section 3, that a steady reaction-diffusion model explains the pathologies observed, characterizes the threshold of instability and leads to a stabilization procedure. By analogy to the simple reaction-diffusion model, these measures are employed to derive a bound for the time step under which noncausal oscillations are expected, and then remove the instabilities from implicit time-integration schemes by spatial stabilization. Conclusions are offered in Section 4.

Formulations for elastodynamics

The equation of motion is considered as a model of a hyperbolic partial differential equation. We examine the Newmark and generalized-families of time-integration schemes, although other algorithms can be considered as well. Consider the following homogeneous, Dirichlet, elastodynamic problem of finding the displacement u(x, t), such that

Semidiscrete formulation of the equation of motion

u ,tt -∇ • (c∇u) = f in × (0, T ), ( 1
)
u = 0 on × (0, T ), (2)
u(x, 0) = u 0 (x), x ∈ , (3)
u ,t (x, 0) = u0 (x), x ∈ . (4)
Here, (x) > 0 is the known density, c(x) is the known fourthrank tensor of elastic coefficients (with the usual symmetry and positive-definiteness properties), f(x, t) is the prescribed body force, u 0 (x) is the specified initial displacement, and u0 (x) is the specified initial velocity. Generalization of the following results to problems with other types of boundary conditions is straightforward (see Section 2.4 for a numerical example with other types of boundary conditions).

For isotropic materials, the tensor c(x) has only two independent parameters, e.g., the Lamé parameters, and . The problem parameters can be combined in terms of the pressure and shear wave velocities, c L and c T , respectively, where c 2 L = ( + 2 )/ and c 2 T = / . Let be partitioned, in the usual way, into nonoverlapping element domains, numbered with index e. The semidiscrete Galerkin approximation is stated in terms of a set of functions that do not depend on time

V h ⊂ H 1 0 ( ). The stan- dard finite element method consists of finding u h (t) ∈ V h , such that ∀ ūh ∈ V h ( ūh , üh ) + a( ūh , u h ) = ( ūh , f), ( 5
) ( ūh , u h (0)) = ( ūh , u 0 ), ( 6
)
( ūh , uh (0)) = ( ūh , u0 ). ( 7
)
Here, x is suppressed as an argument of u, the superposed dot denotes time differentiation, and (•, •) is the L 2 ( ) inner product. The form of the right-hand side assumes sufficiently smooth f. The bilinear operator is

a( ū, u) = (∇ ū, c∇u). ( 8
)
The matrix equations are obtained in the usual way (see, e.g., Ref. [START_REF] Hughes | The Finite Element Method[END_REF]). The functions ūh and u h are expressed in terms of standard finite element shape functions which are not time dependent, whereas the unknown nodal displacements are time dependent.

The semidiscrete form of the initial/boundary-value problem governed by the equation of motion ( 1) is the well-known initial-value problem of finding the vector of unknown nodal displacements, d = d(t), satisfying the coupled system of ordinary differential equations

M d + C ḋ + Kd = F, (9)
and initial conditions

d(0) = d 0 , (10)
ḋ(0) = ḋ0 . ( 11
)
Here, M is the symmetric, positive-definite, mass matrix obtained from the first term on the left-hand side of Eq. ( 5); C is a viscous damping matrix, which is not derived from the equation of motion ( 1), but is often used in structural dynamics; K is the symmetric, positive-semidefinite, stiffness matrix obtained from the second term on the left-hand side of Eq. ( 5); and F = F(t) is the prescribed force vector obtained from the right-hand side of Eq. ( 5). The initial conditions d 0 and ḋ0 are usually taken as nodal values of the given functions u 0 (x) and u0 (x), respectively. The global arrays M, C, K and F are assembled from element arrays m e , c e , k e and f e in the usual way (see, e.g., Ref. [START_REF] Hughes | The Finite Element Method[END_REF]).

The Newmark method

The Newmark method for integrating the semidiscrete equation (9) from time level t n to t n+1 = t n + t is expressed in terms of d n , v n and a n , the approximations to d(t n ), ḋ(t n ) and d(t n ), respectively, as follows:

Ma n+1 + Cv n+1 + Kd n+1 = F n+1 , (12)
d n+1 = dn+1 + 2 t a n+1 , (13)
v n+1 = ṽn+1 + t a n+1 , ( 14
)
where t is the time step and the predictors are

dn+1 = d n + t v n + 2 t 2 (1 -2 )a n , (15)
ṽn+1 = v n + t (1 -)a n . ( 16
)
The Newmark algorithms are a one-step, two-parameter ( and ) family of methods. The second-order accurate, unconditionally stable in time, trapezoidal or average acceleration rule

( = 1 4 , = 1
2 ) is, probably, the most common implicit scheme among the Newmark algorithms. The trapezoidal rule lacks algorithmic damping, which can be obtained by increasing the value of the parameter beyond 1 2 . While this introduces the desired high-frequency dissipation, it also affects the lowfrequency range, resulting in first-order accurate schemes. There is no second-order accurate Newmark method with algorithmic damping.

The following form of implementation is considered:

(M + t C + 2 t K)a n+1 = F n+1 -Cṽ n+1 -K dn+1 . ( 17
)
The initialization of the solution procedure is standard. At the beginning of each time step, the terms on the right-hand side of Eq. ( 17) are known. The equation is solved for a n+1 , and d n+1 and v n+1 are then obtained from the update equations ( 13) and ( 14).

The generalized-method

Consider the integration of the semidiscrete equation ( 9) by the generalized-method for second-order systems [START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: the generalizedmethod[END_REF]. The update equations of the Newmark method ( 13) and ( 14) are retained, whereas the time-discrete equation is modified:

Ma n+1-m + Cv n+1-f + Kd n+1-f = F n+1-f , (18)
d n+1-f = dn+1-f + 2 t (1 -f )a n+1 , (19)
v n+1-f = ṽn+1-f + t (1 -f )a n+1 , (20)
a n+1-m = (1 -m )a n+1 + m a n , (21)
where

t n+1-f = (1 -f )t n+1 + f t n = t n+1 -f t . The predictors are dn+1-f = d n + t (1 -f )v n + 2 t 2 (1 -2 )(1 -f )a n , (22) ṽn+1-f = v n + t (1 -)(1 -f )a n . ( 23
)
The generalized-method is an extension of the Newmark method that exhibits algorithmic damping. It is defined in terms of parameters m and f , in addition to and . Selecting m = f = 0 leads to the familiar Newmark family of methods. Selecting m = 0, f =leads to the HHT-method [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF].

Alternatively, selecting m = 2 ∞ -1 ∞ + 1 , (24) f = ∞ ∞ + 1 , (25) = (1 -m + f ) 2 /4, (26) = (1 -2 m + 2 f )/2, ( 27
)
defines an unconditionally stable in time, second-order accurate, one-parameter (0 ∞ 1), family of methods. The parameter ∞ specifies the dissipation of the high-frequency modal components. Taking ∞ = 0 eliminates the highfrequency response (known as asymptotic annihilation). On the other extreme, setting ∞ = 1 eliminates the algorithmic damping and the trapezoidal time-integration scheme is regained.

The following form of implementation is considered:

((1 -m )M + t (1 -f )C + 2 t (1 -f )K)a n+1 = F n+1-f -m Ma n -Cṽ n+1-f -K dn+1-f . ( 28
)
Again, the initialization of the solution procedure is standard. At the beginning of each time step, the terms on the right-hand side of Eq. ( 28) are known. The equation is solved for a n+1 , and then d n+1 and v n+1 are obtained using the update equations of the Newmark method ( 13) and ( 14)

d n+1 = dn+1-f -f d n + 2 t (1 -f )a n+1 1 -f , ( 29
)
v n+1 = ṽn+1-f -f v n + t (1 -f )a n+1 1 -f . ( 30
)

Motivation: numerical example

Consider a one-dimensional wave propagation problem in a uniform bar of length L and a constant bar wave velocity c 0 = √ E/ , where E is Young's modulus and is the material density. The model lacks physical damping. No body force is applied. The bar is fixed on the left end, and a constant pressure p is applied at t = 0 on the right end (see Fig. 1). The bar is initially at rest. Once the load is applied, a stress wave is generated, propagating towards the fixed end at the bar wave velocity c 0 . Attention is restricted to the time interval before the stress wave reaches the fixed end of the bar (T < L/c 0 ).

The initial/boundary-value problem is

ü -c 2 0 u = 0, 0 < x < L, 0 < t < T , ( 31
) u(0, t) = 0, 0 < t < T , ( 32
) (L, t) = p, 0 < t < T , ( 33
) u(x, 0) = 0, 0 < x < L, ( 34
) u(x, 0) = 0, 0 < x < L, (35)
where = Eu . The exact solution is [18, p. 91]

u = p E c 0 t -(L -x) , ( 36
)
where • is the Macaulay bracket. The bar is at rest before the wave front, since the stress wave has not yet arrived due to the finite velocity of propagation, conforming to the principle of causality. The discretization error inherent in numerical methods invariably leads to precursors before wave fronts that violate the principle of causality. Such precursors are considered acceptable as long as they are monotonic and small. Numerical solutions may also exhibit spurious oscillations on both sides of the wave front resulting from sharp gradients in stress. Noncausal oscillations, appearing before the wave front, are inappropriate. Numerical computations were performed using a uniform mesh of 100 linear rod elements (h = L/100) and three time steps CFL = 1.0, 0.5, 0.1, where CFL = c 0 t /h is the Courant number. Recall, CFL = 1 is the critical time step for temporal stability of the central difference method ( = 0, = 1 4 ), the explicit scheme of the Newmark method (when the mass matrix is lumped).

Snapshots of the numerical solutions for the stress, obtained by the trapezoidal rule and generalized-method ( ∞ = 7 13 ) at t = 0.6L/c 0 , are shown in Fig. 2. The numerical solutions exhibit spurious oscillations resulting from sharp gradients at the wave front. For the largest time step (CFL = 1.0), the solutions exhibit small monotonic precursors but are free of noncausal oscillations, as can be seen in Fig. 2(a). For the trapezoidal rule, spurious oscillations follow the wave front. The damping property of the generalized-method is apparent. For the medium time step (CFL = 0.5), the onset of spurious oscillations before the wave front can be observed in Fig. 2(b), violating the principle of causality. Again, the damping property of the generalized-method is evident, but still noncausal oscillations exist. For the smallest time step (CFL = 0.1), the spurious oscillations before the wave front are no longer localized, as is the case for the medium time step, as can be seen in Fig. 2(c). At this time step, the solutions are essentially converged to the solution of the semidiscrete formulation, hence the two solutions are almost identical.

Spatial stability analysis

The numerical example in the previous section shows that standard time-integration schemes exhibit noncausal phenomena at small time steps. These pathologies can be explained by examining a simple free-space model. Furthermore, this model can be used to formulate stabilized methods that are free of noncausal oscillations at all time steps. Consider the scalar wave equation

ü -c 2 u = f , ( 37
)
where c is the appropriate wave velocity (e.g., c T , c L or c 0 ).

Time-integration by the Newmark method

By the Rothe method, the space-time continuous wave equation ( 37) is discretized in time by the Newmark timeintegration scheme in terms of u n (x), v n (x) and a n (x), the approximations at time t n to u(x, t), u(x, t) and ü(x, t), respectively, as follows:

a n+1 -c 2 u n+1 = f (x, t n+1 ), ( 38
)
u n+1 = ũn+1 + 2 t a n+1 , ( 39
)
v n+1 = ṽn+1 + t a n+1 , ( 40
)
where the predictors are

ũn+1 = u n + t v n + 2 t 2 (1 -2 )a n , (41)
ṽn+1 = v n + t (1 -)a n . ( 42
)
Substitution leads to the following steady reaction-diffusion equation

La n+1 = f , ( 43
)
where Lu = -2 t c 2 u + u is a modified Helmholtz operator and f = f (x, t n+1 ) + c 2 ũn+1 . The degradation of stability for this singularly perturbed elliptic equation at small values of the parameter is well-known.

Stability of spatial Galerkin approximation

The Galerkin approximation of the free-space model (neglecting boundary conditions) consists of finding a h n+1 , such that

( āh n+1 , a h n+1 ) + a( āh n+1 , a h n+1 ) = ( āh n+1 , f ). ( 44
)
In order to recover the usual discrete form, the predictor term on the right-hand side should be integrated by parts. The bilinear operator in this case is

a( ū, u) = 2 t (∇ ū, c 2 ∇u). (45)
Free-space solutions of the constant-coefficient, homogeneous, reaction-diffusion equation exhibit exponential growth and decay with rate 1/( t c), which is related to the Courant number (CFL) in computation. Von Neumann analysis of the Galerkin method of a uniform d-dimensional mesh of linear elements of size h aligned with the direction of growth or decay leads to an expression for the approximate growth or decay rate (for more details, see Ref. [START_REF] Harari | Stability of semidiscrete formulations for parabolic problems at small time steps[END_REF]). Expressing the approximate exponential rate as CFL h , the Galerkin method results in the following relation:

1 CFL h = arccosh 6 CFL 2 + 2 6 CFL 2 -1 . ( 46
)
This expression represents the correct behavior of exponential growth and decay (CFL h ∈ R) when CFL 2 > 1/(6 ). For CFL < 1/ 6 , CFL h is complex valued, which represents spurious oscillations. This is a manifestation of degradation of spatial stability of the time-discrete problem, in the H 1 sense, as the value of the CFL decreases, while maintaining L 2 stability, which is incapable of controlling derivatives. Considerable degradation in accuracy is possible (Fig. 3) even prior to the potential onset of spurious oscillations at CFL = 1/ 6 .

β 1/2 CFL (CFL h -CFL) / CFL Lumped Consistent
Fig. 3. Error in growth and decay rate for the Galerkin method (solid line) and for lumped mass method (dotted line) prior to potential onset of spurious oscillations.

Remark 1. This analysis holds for consistent representation of inertial terms. With lumped mass, the expression for the approximate growth and decay rate becomes

1 CFL h = arccosh 1 + 1 2 CFL 2 . ( 47
)
This relation represents exponential growth and decay for all admissible values of CFL, and hence solutions are free of spurious oscillations (see Fig. 3). However, this may not be the case for lumped mass elements of high-order. For the linear elements considered here, the error of the approximate growth and decay rate, related to the local truncation error, is of the same order as in the consistent representation on uniform meshes. The consistent representation retains second-order spatial accuracy in the local truncation error on non-uniform meshes, whereas the lumped representation can degrade to first order [START_REF] Harari | Accurate finite difference methods for time-harmonic wave propagation[END_REF]. Moreover, the loss of guaranteed bounds, when the variational framework is violated by mass lumping, can be detrimental to the use of techniques such as error estimation.

In summary, the Galerkin solution with linear elements may exhibit spurious spatial oscillations for

CFL < 1 6 , ( 48
)
and considerable degradation in accuracy even prior to the potential onset of the spurious oscillations. Note that the bound for the time step is O(h). For the trapezoidal rule ( = 1 4 , = 1 2 ), the bound for spatial stability is CFL < 2 3 .

Stabilized method

A stabilized finite element method for the Newmark algorithm is obtained by appending to the Galerkin equation (44) terms containing residual-based operators multiplied by method parameters :

( āh n+1 , a h n+1 ) + a( āh n+1 , a h n+1 ) -(L * āh n+1 , La h n+1 ) = ( āh n+1 , f ) -(L * āh n+1 , f ) , (49)
where L * is the adjoint operator of the differential operator L given in Eq. ( 43). The mesh-dependent stabilization parameter for linear elements of size h,

= 1 + 6 CFL 2 1 -cosh(1/( CFL)) 2 + cosh(1/( CFL)) , (50)
is defined on the element level. We note for later reference that 0 < 1 for CFL > 0 (see Fig. 4). Subscripts on inner products denote domains of integration other than , namely, denotes the union of element interiors such that = . Integration of the additional terms is performed over in order to respect regularity requirements of typical, piecewise smooth, finite element functions. In practice, standard finite element procedures which assemble global arrays from element contributions are employed without modification. For linear elements, spatial stabilization is obtained by a simple modification of the element mass matrix and force vector prior to assembly, namely

m e ← (1 -)m e , ( 51
)
f e ← (1 -)f e . ( 52
)
Recall, the predictor terms in the force vector are integrated by parts. Symmetry of the spatially stabilized mass matrix follows directly from the definition of the modified matrix. Positivedefiniteness of the modified mass matrix follows along the lines of the procedure for the original matrix, recalling that 1-> 0.

The analysis of the spatially stabilized Newmark method is similar to that of the conventional method [5, pp. 497-498], indicating that spatial stabilization does not effect unconditional temporal stability.

Repeating the analysis in Section 3.1.1 for the stabilized method (49) shows that setting the method parameter as (50) leads to CFL h = CFL, eliminating errors in the exponential rate in Cartesian meshes of linear elements aligned in the direction of growth and decay, along with the attendant spurious oscillations. Stabilization restores the correct balance between the L 2 and H 1 terms. Furthermore, it reduces the sensitivity of the computation to spurious phenomena that are generated by transitions in mesh size [START_REF] Harari | Reducing spurious dispersion, anisotropy and reflection in finite element analysis of time-harmonic acoustics[END_REF].

Time-integration by the generalized-method

The wave equation ( 37) is now discretized, by the Rothe method, with the generalized-time-integration scheme, as follows:

a n+1-m -c 2 u n+1-f = f (x, t n+1-f ), ( 53
)
u n+1-f = (1 -f )u n+1 + f u n , ( 54
)
v n+1-f = (1 -f )v n+1 + f v n , (55)
a n+1-m = (1 -m )a n+1 + m a n , (56)
where

u n+1 = ũn+1-f -f u n + 2 t (1 -f )a n+1 1 -f , ( 57
)
v n+1 = ṽn+1-f -f v n + t (1 -f )a n+1 1 -f (58) and ũn+1-f = u n + t (1 -f )v n + 2 t 2 (1 -2 )(1 -f )a n , (59) ṽn+1-f = v n + t (1 -)(1 -f )a n . ( 60
)
As in the case of the Newmark method, simple manipulation yields a steady reaction-diffusion equation for a n+1 . Eq. ( 43) is regained with Lu = -((1f )/(1m )) 2 t c 2 u + u, the modified Helmholtz operator, and f = f (x, t n+1-f )m a n + c 2 ũn+1f . Once again, the predictor term of f is integrated by parts in the weak form.

Repeating the analysis in Section 3.1.1 leads to a bound for the time step, under which noncausal spatial oscillations may occur. Specifically, noncausal spatial oscillations may arise in the standard one-parameter ( ∞ ) family of generalizedmethods for

CFL < (1 + ∞ ) 2 -∞ 6 . (61)
Again, the bound for the time step is O(h). Recall, the parameter ∞ specifies the dissipation of the high-frequency modal components. Setting ∞ = 1, eliminating the algorithmic damping, recovers the bound for spatial stability of the trapezoidal time-integration scheme CFL < 2 3 . On the other extreme, taking ∞ = 0, eliminating the highfrequency response, improves the bound for spatial stability CFL < 1 3 compared to the trapezoidal rule. For ∞ = 7 13 (from the numerical example) the bound for spatial stability is CFL < 0.76. A stabilized version of the generalized-algorithm is obtained by the procedure outlined in Section 3.1.2. The mesh-dependent stabilization parameter for linear elements of size h in this case is obtained by substituting in Eq. (50) 1/((2 -∞ )(1 + ∞ ) 2 ) for .

Numerical example revisited

Returning to the numerical example described in Section 2.4, the stability analysis, with c 0 as the appropriate wave velocity, shows that the largest time step (CFL = 1.0) is higher than the bounds for spatial stability (48) and (61) for both algorithms considered ( = 1 4 , ∞ = 7 13 ). Indeed, the solutions for the largest time step are free of noncausal oscillations (see Fig. 2(a)). On the other hand, the medium time step (CFL = 0.5) is a little below and the smallest time step (CFL = 0.1) is well below these bounds, explaining the noncausal oscillations observed in Figs. 2(b) and(c).

Snapshots of the numerical solutions for the stress, obtained by the stabilized trapezoidal rule and stabilized generalizedmethod at t = 0.6L/c 0 , are shown in Fig. 5. The numerical solutions still exhibit spurious oscillations resulting from sharp gradients at the wave front. However, spurious oscillations before the wave front are absent from the solutions of the stabilized methods, even for time steps lower than the bounds for spatial stability (48) and (61), thus conforming to the principle of causality. The solutions of the stabilized trapezoidal rule for all three time steps are similar to the solution of the unstabilized trapezoidal rule for the largest time step. The damping property of the stabilized generalized-method is apparent, providing the best performance.

While the proposed approach successfully removes noncausal oscillations from the computation, considerable spurious oscillations remain. The combination of spatial stabilization and algorithmic damping is probably incapable of addressing this deficiency, which would require combining spatial stabilization with more advanced time-marching schemes. Initialization procedures and alternative forms of implementation should also be examined.

Conclusions

Solutions of semidiscrete formulations for hyperbolic problems may exhibit noncausal oscillations, due to poor approximation of higher modes. Consequently, all conventional direct time-integration schemes that are based on such semidiscrete formulations and are convergent in time may eventually admit these pathologies as the time step is reduced (with a fixed mesh size).

In this paper, a procedure for determining specific values for the potential onset of noncausal spatial oscillations is provided by considering the scalar wave equation as a simple model. Time discretization by the Rothe method leads to a singularly perturbed, steady, reaction-diffusion equation. In particular, bounds (for linear elements) are obtained for implicit algorithms belonging to the Newmark and generalizedmethods.

The time-discrete problem can be stabilized by any procedure that is used for the singularly perturbed equation. For linear elements, spatial stabilization is obtained by a simple modification of the element mass matrix and force vector prior to assembly, in order to improve spatial performance within each time step. This procedure preserves the property of unconditional temporal stability. Employing spatial stabilization, as proposed, eliminates noncausal oscillations observed in implicit computations with small time steps. The best performance is obtained by a combination of spatial stabilization and algorithmic damping, as in the case of the stabilized version of the generalizedmethod.

The onset of noncausal spatial oscillations in implicit timeintegration occurs at time steps of the order of the critical time step for temporal stability of corresponding explicit schemes. Consequently, this pathology may be of particular concern in the application of implicit-explicit mesh partitions. The stabilization procedure proposed in this work is well suited for the implicit part of such algorithms.

  Let ⊂ R d be a d-dimensional, open, bounded domain with smooth boundary . The open time interval is (0, T ).

Fig. 1 .

 1 Fig. 1. Bar fixed on the left end and subject to pressure on the right end.

Fig. 2 .

 2 Fig. 2. Snapshots of stress computed by the trapezoidal rule and generalized-method ( ∞ = 7 13 ) at t = 0.6L/c 0 (compared to the exact solution). (a) Solutions at CFL = 1.0: trapezoidal rule (left) and generalized-method (right). (b) Solutions at CFL = 0.5: trapezoidal rule (left) and generalized-method (right). (c) Solutions at CFL = 0.1: trapezoidal rule (left) and generalized-method (right).

Fig. 4 .

 4 Fig. 4. Stabilization parameter for linear elements.

Fig. 5 .

 5 Fig. 5. Snapshots of stress computed by the stabilized trapezoidal rule and stabilized generalized-method ( ∞ = 7 13 ) at t = 0.6L/c 0 (compared to the exact and unstabilized solutions). (a) Solutions at CFL = 1.0: stabilized trapezoidal rule (left) and stabilized generalized-method (right). (b) Solutions at CFL = 0.5: stabilized trapezoidal rule (left) and stabilized generalized-method (right). (c) Solutions at CFL = 0.1: stabilized trapezoidal rule (left) and stabilized generalized-method (right).


Guanwei Luo

Li Ma

Xiaohong Lv

Dynamic analysis and suppressing chaotic impacts of a two-degree-of-freedom oscillator with a clearance

Keywords: Vibration, Impact, Periodic motion, Bifurcation, Suppressing chaos

Dynamic analysis and suppressing chaotic impacts of a two-degree-of-freedom oscillator with a clearance

Guanwei Luo, Li Ma, Xiaohong

Introduction

Vibrating systems with clearances between the moving parts are frequently encountered in technical applications. Impacts occur when the amplitudes of vibration of some parts of the system exceed critical values. Some examples are linkage machines, heat exchangers, fuel rods in nuclear power plants, rolling railway wheelset, piping systems, gear transmissions and so on. Impacts inevitably occur for the vibrating systems with clearances or gaps, but they are undesirable as they bring about failures, strain, shorter service life and increased noise levels. It is important to be able to accurately model the dynamics of an impacting system, so as to enlarge the profitable effects such as optimum design and high reliability, to minimize adverse effects such as high noise levels, pitting, wear and tear and so on. Compared with single impact, the nonlinear dynamics of vibro-impact systems is more complicated. The trajectories of such systems in phase space have discontinuities caused by the impacts. Consequently, the presence of the nonlinearity and discontinuity complicates the dynamic analysis of such systems considerably, but they can be described theoretically and numerically by discontinuities in good agreement with reality. The broad interest in analyzing and understanding the performance of such systems is reflected by a still increasing amount of investigations devoted to this area. Several methods of the theoretical analysis have been developed and different models of impacts have been assumed in the past several years. Stability and bifurcations of different types of impact oscillators were reported in Refs. [START_REF] Shaw | A periodically forced piecewise linear oscillator[END_REF][START_REF] Luo | Codimension two bifurcations of fixed points in a class of vibratory systems with rigid stops[END_REF][START_REF] Aidanpää | Periodic and chaotic behaviour of a threshold-limited two-degree-of-freedom system[END_REF][START_REF] Quinn | Near-simultaneous impacts[END_REF][START_REF] Knudsen | Dynamic stability of weakly damped oscillators with elastic impacts and wear[END_REF][START_REF] Albert | Period-doubling induced chaotic motion in the LR model of a horizontal impact oscillator[END_REF][START_REF] Dimentberg | Stochastic and/or chaotic response of a vibration system to imperfectly periodic sinusoidal excitation[END_REF][START_REF] De Souza | Basins of attraction changes by amplitude constraining of oscillators with limited power supply[END_REF][START_REF] Wagg | Multiple non-smooth events in multi-degree-of-freedom vibro-impact systems[END_REF][START_REF] Luo | Hopf-flip bifurcations of vibratory systems with impacts[END_REF][START_REF] Luo | Codimension two bifurcations of fixed points in a class of vibratory systems with rigid stops[END_REF]. Blazejczyk-Okolewska et al. [START_REF] Blazejczyk-Okolewska | Determination of geometrical conditions of assembly and impacts in classified types of mechanical systems with impacts[END_REF] provided much information of the fundamental nature that broadens the scope of knowledge on the motion of mechanical systems with impacts. A special feature of impacting systems is the instability caused by grazing bifurcation. The first important work in this area was done by Nordmark [START_REF] Nordmark | Non-periodic motion caused by grazing incidence in an impact oscillator[END_REF]. This work has been further expanded by thorough investigations of two-dimensional maps, where some universal behavior has been found [START_REF] Peterka | Transition to chaotic motion in mechanical systems with impacts[END_REF][START_REF] Whiston | Singularities in virbo-impact dynamics[END_REF][START_REF] Foale | Dynamical complexities of forced impacting systems[END_REF][START_REF] Ivanov | Stabilization of an impact oscillator near grazing incidence owing to resonance[END_REF][START_REF] Budd | The effect of frequency and clearance variations on single-degree-of-freedom impact oscillators[END_REF][START_REF] Ivanov | Bifurcation in impact systems[END_REF][START_REF] Bernardo | Local analysis of C-bifurcations in N-dimensional piecewise-smooth dynamical systems[END_REF][START_REF] Murphy | Grazing instabilities and post-bifurcation behavior in an impacting string[END_REF][START_REF] Halse | C-bifurcations and period-adding in one-dimensional piecewise-smooth maps[END_REF][START_REF] Albert | Periodic motions and grazing in a harmonically forced, piecewise, linear oscillator with impacts[END_REF]. Based on the local singularity theory, Luo and Gegg [START_REF] Albert | Grazing phenomena in a periodically forced, friction-induced, linear oscillator[END_REF] provided a comprehensive investigation of grazing motions in the dry-friction oscillator for a better understanding of the grazing mechanism of such a discontinuous system. Pavlovskaia and Wiercigroch [START_REF] Pavlovskaia | Modeling of ground moling dynamics by an impact oscillator with a frictional slider[END_REF][START_REF] Pavlovskaia | Two-dimensional map for impact oscillator with drift[END_REF][START_REF] Pavlovskaia | Modeling of an impact system with a drift[END_REF][START_REF] Wiercigroch | Material removal rate prediction for ultrasonic drilling of hard materials using an impact oscillator approach[END_REF] developed a series of mathematical models of impact systems with drift and revealed very complex dynamical behavior of such systems. Souza and Caldas [START_REF] De Souza | Calculation of Lyapunov exponents in systems with impacts[END_REF] applied a model-based algorithm for the calculation of the spectrum of the Lyapunov exponents of attractors of mechanical systems with impacts. Peterka [START_REF] Peterka | Explanation of appearance and characteristics of intermittency chaos of the impact oscillator[END_REF] found chaotic motion of an intermittency type of the impact oscillator appearing near segments of saddle-node stability boundaries of subharmonic motions with two different impacts in motion period. Luo [START_REF] Albert | Grazing and chaos in a periodically forced, piecewise linear system[END_REF] presented an idealized, piecewise linear system to model the vibration of gear transmission systems, and predicted analytically the periodic motions in a generalized, piecewise linear oscillator with perfectly plastic impacts. Wagg [START_REF] Wagg | Rising phenomena and the multi-sliding bifurcation in a two-degree of freedom impact oscillator[END_REF][START_REF] Wagg | Periodic sticking motion in a two-degree-of-freedom impact oscillator[END_REF] analyzed the chattering impact and rising phenomena which occur in the stick solutions of impact oscillators. Luo and Gegg [START_REF] Albert | Stick and non-stick periodic motions in periodically forced oscillators with dry friction[END_REF] developed the force criteria for stick and non-stick motions in harmonically forced, friction-induced oscillators from the local theory of non-smooth dynamical systems on connectable domains. Several effective methods of controlling chaos and experimental analysis have been developed in Refs. [START_REF] Hu | Controlling chaos of a periodically forced nonsmooth mechanical system[END_REF][START_REF] Hu | Controlling chaos of a dynamical system with discontinuous vector field[END_REF][START_REF] Lee | Control of impact oscillator[END_REF][START_REF] De Souza | Controlling chaotic orbits in mechanical systems with impacts[END_REF][START_REF] De Souza | Suppressing grazing chaos in impacting system by structural nonlinearity[END_REF][START_REF] Zinjade | Impact damper for controlling friction-driven oscillations[END_REF][START_REF] Zhao | Random Fatigue and Clearance Nonlinearity of Mechanical Vibrating Systems[END_REF][START_REF] Wiercigroch | Experimental study of base excited symmetrically piecewise linear oscillator[END_REF]. Hu [START_REF] Hu | Controlling chaos of a dynamical system with discontinuous vector field[END_REF] presented how to control the chaos of dynamical systems with discontinuous vector field through the paradigm of a harmonically forced oscillator having a set-up elastic stop. Lee and Yan [START_REF] Lee | Control of impact oscillator[END_REF] studied the algorithms of position control of impact oscillator and demonstrated synchronization of two impact oscillators. Souza et al. stabilized desired unstable periodic orbits, embedded in the chaotic invariant sets of mechanical systems with impacts, by applying a small and precise perturbation on an available control parameter [START_REF] De Souza | Controlling chaotic orbits in mechanical systems with impacts[END_REF], and investigated the influence of structural nonlinearity of a simple cantilever beam impacting system on its dynamic responses close to grazing incidence by means of numerical simulation [START_REF] De Souza | Suppressing grazing chaos in impacting system by structural nonlinearity[END_REF]. Zinjade and Mallik [START_REF] Zinjade | Impact damper for controlling friction-driven oscillations[END_REF] designed and fabricated an experimental setup of a single-degree-of-freedom friction-driven oscillator with an attached impact damper, and investigated the control of friction-driven oscillations by using an impact damper. The effects of mass ratio, coefficient of restitution and clearance on the performance of the impact damper have been investigated analytically and the results are verified by numerical integration. Zhao and Wang [START_REF] Zhao | Random Fatigue and Clearance Nonlinearity of Mechanical Vibrating Systems[END_REF] proposed several feedback control methods to suppress chaotic behavior of an impact damper. Experimental study of base excited symmetrically piecewise linear oscillator was performed by Wiercigroch and Sin [START_REF] Wiercigroch | Experimental study of base excited symmetrically piecewise linear oscillator[END_REF]. Along with the basic research into vibro-impact dynamics, a wide range of impacting models have been applied to simulate and analyze engineering systems operating within bounded dynamic responses. For example, in wheel-rail impacts of railway coaches [START_REF] Meijaard | Railway vehicle systems dynamics and chaotic vibrations[END_REF][START_REF] Luo | Periodic-impact motions and bifurcations of the vibratory system with a clearance[END_REF], vibrating hammer [START_REF] Xie | The mathematical model for the impact hammer and global bifurcations[END_REF], pile driver [START_REF] Luo | Dynamics of a small vibro-impact pile driver[END_REF], Jeffcott rotor with bearing clearance [START_REF] Pavlovskaia | Non-linear dynamic interactions of a Jeffcott rotor with preloaded snubber ring[END_REF][START_REF] Karpenko | Piecewise approximate analytical solutions for a Jeffcott rotor with a snubber ring[END_REF][START_REF] Karpenko | Experimental verification of Jeffcott rotor model with preloaded snubber ring[END_REF], excited pendula with impacts [START_REF] Quinn | The dynamics of two parametrically excited pendula with impacts[END_REF], impact dampers [START_REF] Sung | Dynamics of harmonically excited impact damper: Bifurcations and chaotic motion[END_REF][START_REF] Han | Chaotic motion of a horizontal impact pair[END_REF][START_REF] Heiman | Periodic motions and bifurcations in dynamics of an inclined impact pair[END_REF][START_REF] Bapat | The general motion of an inclined impact damper with friction[END_REF][START_REF] Duncan | The damping performance of a single particle impact damper[END_REF], gears [START_REF] Albert | Grazing and chaos in a periodically forced, piecewise linear system[END_REF][START_REF] Kahraman | Non-linear dynamics of a geared rotor-bearing system with multiple clearances[END_REF][START_REF] Kunert | Stochastic model for rattling in gear-boxes[END_REF], etc., impacting models have proved to be useful [START_REF] Blazejczyk-Okolewska | [END_REF].

A two-degree-of-freedom impact oscillator with a clearance is considered in the paper. The maximum displacement of one of the masses is limited to a threshold value by the symmetrical rigid stops. Dynamics of the system is studied with special attention to periodic-impact motions and bifurcations. Period-one double-impact symmetrical motion and transcendental impact map of the system are derived analytically. Stability and local bifurcations of the period-one double-impact symmetrical motions are analyzed by using the impact Poincaré map. The Lyapunov exponents in the impact oscillator with a clearance are calculated by using the transcendental impact map. The influence of the clearance and excitation frequency on the dynamics of the system is analyzed. The bifurcation regularity of periodone double-impact symmetrical motion is revealed. The chattering-impact vibration and the sticking phenomena are found to occur in the region of low forcing frequency. These imply that the dynamics of this system is very complex, varying from different types of periodic motions to chaos, even chattering-impacting vibration and sticking. Chaoticimpact motions are suppressed to minimize the adverse effects by using external driving force, delay feedback and feedback-based method of period pulse.

Mechanical model of a harmonically excited system having symmetrical rigid stops

A two-degree-of-freedom system having symmetrically placed rigid stops and subjected to periodic excitation is shown schematically in Fig. 1. Displacements of the masses M 1 and M 2 are represented by X 1 and X 2 , respectively. The masses are connected to the supporting base by linear springs with stiffness K 1 and K 2 , and linear viscous dashpots with damping constants C 1 and C 2 . Damping in the mechanical model is assumed as proportional damping. The excitations on both masses are harmonic with amplitudes P 1 and P 2 . The excitation frequency Ω and the phase τ are the same for both masses. The masses move only in the horizontal direction. For small forcing amplitudes the system will undergo simple oscillations and behave as a linear system. As the amplitude is increased, the mass M 1 eventually begins to hit the stops and the motion becomes nonlinear (the other mass is not allowed to impact any stop). The impact is described by a coefficient of restitution r .

The motion processes of the system, between consecutive impacts occurring at the stop A, are considered. Between any two consecutive impacts, the time T is always set to zero directly at the starting point (the stop A), and the phase angle τ is used only to make a suitable choice for the origin of time in the calculation. The state of the vibroimpact system, immediately after impact, has become new initial conditions in the subsequent process of the motion. Consequently, the motion of this system is a combination of smooth motions governed by linear differential equation interrupted by a series of non-smooth impacts. The dynamic analysis can be considered in a non-dimensional form, and all parameters and variables are assumed to be non-dimensional, which are introduced by

µ m = M 2 M 1 + M 2 , µ k = K 2 K 1 + K 2 , µ c = C 2 C 1 + C 2 , f 20 = P 2 P 1 + P 2 , ω = Ω M 1 K 1 , ζ = C 1 2 √ K 1 M 1 , δ = B K 1 P 1 + P 2 , t = T K 1 M 1 , x i = X i K 1 P 1 + P 2 . (1)
The governing equations of the system away from the impact can be expressed as

1 0 0 µ m 1 -µ m ẍ1 ẍ2 +   2ζ -2ζ -2ζ 2ζ 1 + µ c   ẋ1 ẋ2 +   1 -1 -1 1 1 -µ k   x 1 x 2 = 1 -f 20 f 20 sin(ωt + τ ), |x 1 | < δ. (2)
An impact occurs wherever |x 1 (t)| = δ. This type of impact may be modelled using an instantaneous coefficient of restitution rule such that ẋ1A+ = -r ẋ1A-, (x 1 = δ), ẋ1 Ā+ = -r ẋ1 Ā-, (x 1 = -δ).

(

) 3
In Eqs. ( 2) and (3), a dot (•) denotes differentiation with respect to the non-dimensional time t, the subscript minus sign denotes the states just before the impact and the subscript plus sign denotes the states just after the impact.

Impact Poincaré map and period-one double-impact symmetrical motion

The impacting system is conveniently studied by the use of an impact Poincaré map derived from the equations of motion, supplemented by transition conditions at the instants of impacts. Here each iterate of this map corresponds to the mass M 1 impacting the stop A once. The state of the vibro-impact system, just immediately after the impact occurring at the stop A, is chosen as the Poincaré section, i.e., σ = (x 1 , ẋ1 , x 2 , ẋ2 , θ) ∈ R 4 × S, x 1 = δ, ẋ1 = ẋ1+ . As a result, the impact Poincaré map is expressed briefly by

X = f (v, X ), (4)
where

X ∈ R 4 , v is a bifurcation parameter, v ∈ R 1 , X = ( ẋ1+ , x 2 , ẋ2 , τ ) T , X = ( ẋ 1+ , x 2 , ẋ 2 , τ ) T .
The periodic motions of the vibro-impact system are characterized by the symbol Npq, where q and p are the numbers of impacts occurring respectively at the stops A and Ā, and N is the number of the forcing cycles. In this section, the periodic motion of the model, with two symmetrical impacts per force cycle, is emphasized, which is called the period-one double-impact symmetrical motion, i.e., 1-1-1 symmetrical motion. The special periodic motion means that if the dimensionless time t is set to zero directly just after an impact occurring at the stopA, it becomes 2π/ω just before the next impact occurring at the same stop. It certainly takes the same time for the mass M 1 to move from the stop A to Ā and from the stop Ā to A, i.e., t 1 = t 2 = π/ω. After the origin of t-coordinate is displaced to an impact point, the determination of period-one double-impact symmetrical motion is based on the fact that it satisfies the following set of periodicity and matching conditions

x(0) ẋ(0) = - x(π/ω) + ẋ(π/ω) + = D x(2π/ω) ẋ(2π/ω) = x 0 ẋ0 , x(π/ω) - ẋ(π/ω) - = - x(2π/ω) ẋ(2π/ω) , (5)
where

D = diag[d i ], (d i = 1, i = 1, 2, 4; d 3 = -r ), x(t) = (x 1 (t), x 2 (t)) T , x 0 = (δ, x 20 ) T , ẋ0 = ( ẋ1+ , ẋ20 ) T .
Combining the general solution of Eq. ( 2) with the set of periodicity and matching conditions, we can find the totalanalytical expressions for period-one double-impact symmetrical motion of the system. Inserting the time t = 0 into the periodic solution, we can obtain the fixed point X * = ( ẋ1+ , x 10 , ẋ2+ , τ 0 ) T , which is the projection of period-one double-impact symmetrical orbit to the Poincaré section σ .

Jacobian matrix of the transcendental map

Eq. ( 2) is amenable to analytical treatment due to proportional damping. Between the n and (n + 1)th impacts, the general solutions of Eq. (2) are give by

x i (t) = 2 j=1 ψ i j (e -η j t (a jn cos ω d j t + b jn sin ω d j t) + A j sin(ωt + τ n ) + B j cos(ωt + τ n )) (6) ẋi (t) = 2 j=1 ψ i j (e -η j t ((b jn ω d j -η j a jn ) cos ω d j t -(η j b jn + a jn ω d j ) sin ω d j t) + A j ω cos(ωt + τ n ) -B j ω sin(ωt + τ n )), (7)
where i = 1, 2, ω j ( j = 1, 2) denote the eigenfrequencies of Eq. ( 2), ψ i j are the components of the canonical model matrix Ψ ; modal damping ratio η j = ζ ω 2 j , damping eigenfrequency ω d j = ω 2 j -η 2 j ; a jn and b jn are the constants of integration determined by model parameters and the state variables X n just after the nth impact, which are expressed by

a 1n = ψ 22 x 1n -ψ 12 x 2n ψ 22 ψ 11 -ψ 12 ψ 21 -A 1 sin τ n -B 1 cos τ n , (8)
a 2n = ψ 21 x 1n -ψ 11 x 2n ψ 21 ψ 12 -ψ 11 ψ 22 -A 2 sin τ n -B 2 cos τ n , (9)
b 1n = ψ 22 ẋ1n+ -ψ 12 ẋ2n ψ 22 ψ 11 -ψ 12 ψ 21 + η 1 a 1n -A 1 ω cos τ n + B 1 ω sin τ n ω d1 , (10)
b 2n = ψ 21 ẋ1n+ -ψ 11 ẋ2n ψ 21 ψ 12 -ψ 11 ψ 22 + η 2 a 2n -A 2 ω cos τ n + B 2 ω sin τ n ω d2 , (11)
where x 1n = δ or -δ, X n = ( ẋ1n+ , x 2n , ẋ2n , τ n ) T denote the state variables just after the impact occurring at the stop A or Ā.

Because of the existence of the impact, there is discontinuity in the velocity ẋ1 . Due to these impacts, we cannot use the algorithms usually applied to compute Lyapunov exponents in smooth systems. Thus, we introduce a transcendental map to develop an adapted algorithm that can be applied to the vibro-impact system described in this section. For the vibro-impact system, the discontinuities in the trajectory have to be supplemented with matching conditions at the impact instants. Such instants are determined by integrating the equations of motion in continuous time. The values of velocities and displacements of two masses and the phase shift at this time are taken as the dynamical variables of the transcendental map.

ẋ1(n+1)+ = -r 2 j=1 ψ 1 j (e -η j t n+1 ((b jn ω d j -η j a jn ) cos ω d j t n+1 -(η j b jn + a jn ω d j ) sin ω d j t n+1 ) + A j ω cos(ωt n+1 + τ n ) -B j ω sin(ωt n+1 + τ n )), (12)
x 2(n+1) = 2 j=1 ψ 2 j (e -η j t n+1 (a jn cos ω d j t n+1 + b jn sin ω d j t n+1 )

+ A j sin(ωt n+1 + τ n ) + B j cos(ωt n+1 + τ n )), (13)
ẋ2(n+1) = 2 j=1 ψ 2 j (e -η j t n+1 ((b jn ω d j -η j a jn ) cos ω d j t n+1 -(η j b jn + a jn ω d j ) sin ω d j t n+1 ) + A j ω cos(ωt n+1 + τ n ) -B j ω sin(ωt n+1 + τ n )), (14)
τ n+1 = ωt n+1 + τ n -2π, (ωt n+1 + τ n ≥ 2π ) ωt n+1 + τ n , (ωt n+1 + τ n < 2π ), (15)
where x 1(n+1) = δ or -δ.

Using the map we can obtain the Jacobian matrix components necessary to calculate the Lyapunov exponents. Using Eqs. ( 12)-( 15) we develop the algorithm used to obtain the Lyapunov exponents. To calculate the Lyapunov exponents, we use the Jacobian matrix J n of the transcendental map

J n =           ∂ ẋ1(n+1)+ ∂ ẋ1n+ ∂ ẋ1(n+1)+ ∂ x 2n ∂ ẋ1(n+1)+ ∂ ẋ2n ∂ ẋ1(n+1)+ ∂τ n ∂ x 2(n+1) ∂ ẋ1n+ ∂ x 2(n+1) ∂ x 2n ∂ x 2(n+1) ∂ ẋ2n ∂ x 2(n+1) ∂τ n ∂ ẋ2(n+1) ∂ ẋ1n+ ∂ ẋ2(n+1) ∂ x 2n ∂ ẋ2(n+1) ∂ ẋ2n ∂ ẋ2(n+1) ∂τ n ∂τ n+1 ∂ ẋ1n+ ∂τ n+1 ∂ x 2n ∂τ n+1 ∂ ẋ2n ∂τ n+1 ∂τ n           . ( 16
)
The derivation of the components in matrix ( 16) is introduced in the Appendix.

Lyapunov exponents

The computation algorithm is described to calculate the Lyapunov exponents in vibro-impact systems using a transcendental map in Ref. [START_REF] De Souza | Calculation of Lyapunov exponents in systems with impacts[END_REF]. Consider m-dimensional discrete dynamical system

X k+1 = F(X k ), X k ∈ R m . ( 17
)
The Lyapunov exponents are expressed by where |Λ N j | are eigenvalues of the Jacobian matrix M given by

λ j = lim N →∞ 1 N ln |Λ N j |, j = 1, 2, . . . . . . , m, (18)
M = N n=1 J n , ( 19
)
where J n is the Jacobian matrix of the nth iteration (impact) of the transcendental map.

As have been mentioned in Ref. [START_REF] De Souza | Calculation of Lyapunov exponents in systems with impacts[END_REF], Eq. ( 18) cannot be used in the computation algorithm due to the matrix M obtained by the product [START_REF] Ivanov | Bifurcation in impact systems[END_REF]. Even for a few iterations, the components of matrix M become very large for chaotic attractors and null for the periodic attractors so that the numerical calculation becomes impracticable. To eliminate this problem, a way is provided to convert the matrix M into a product of the orthogonal matrices and the upper right triangular matrices [START_REF] De Souza | Calculation of Lyapunov exponents in systems with impacts[END_REF], whose eigenvalues are the numbers along the diagonal and the Lyapunov exponents can be expressed by

λ j = lim N →∞ 1 N N n=1 ln |a n j |, j = 1, 2, . . . . . . m, (20)
where a n j are diagonal components of nth iteration of the product of the upper right triangular matrices. Thus, the Lyapunov exponents are obtained from a sum in (20) instead of a product.

Furthermore, we cannot perform exactly the infinite sum in [START_REF] Bernardo | Local analysis of C-bifurcations in N-dimensional piecewise-smooth dynamical systems[END_REF] for numerical calculations, but the convergence of the Lyapunov exponents occurs with good precision after just a few thousand iterations, both for periodic and chaotic attractors.

6. Suppressing chaos in the vibro-impact systems 6.1. Suppressing chaos in the vibro-impact systems by sine periodic force Sine periodic force is a driving one which is easy to apply and control. Zhao and Wang [START_REF] Zhao | Random Fatigue and Clearance Nonlinearity of Mechanical Vibrating Systems[END_REF] presented a method of controlling bifurcation and chaos by adjusting the external driving force. The method is proved to be effective by the presentation of examples of suppressing bifurcation and chaos for an impact damper. Chaotic impacts of vibratory systems with a clearance are suppressed by applying the controlling method in the paper. The schematic of the method for suppressing chaos by using sine periodic force is shown in Fig. 2 [START_REF] Zhao | Random Fatigue and Clearance Nonlinearity of Mechanical Vibrating Systems[END_REF]. The desired target is some measurable state variable while the system is in a stable state, as for example, impact period, displacements or velocities, etc. The desired target is preserved as an ideal output. The practical output in Fig. 2 is the real time test value of measurable state variable. e n+1 is the error between the practical target and the desired target at n +1 time. u is the external driving force. The external sine driving force is expressed by

u = B sin(γ ωt), (21)
where B and γ are adjustable parameters controlled by the adjustment algorithm of external driving force. The error-transfinite decision device satisfies the relationships: when |e n+1 | < ∆, e = 0, there is no control; when |e n+1 | > ∆, e = e n+1 , the control is applied. The symbol ∆ represents the threshold value of allowable deviation of the system.

For Fig. 2, if the error between the practical target and the desired target satisfies the condition |e n+1 | < ∆, e = 0, the adjustment algorithm makes the value of u equal zero. As a result, the external driving force is not applied. If |e n+1 | > ∆, e = e n+1 is produced by the adjustment algorithm of external driving force so that its parameters, B and γ , are dynamically adjusted until the system is in a stable state, i.e., e = 0.

To this end, the adjustable parameters B and γ of driving signal u = B sin(γ ωt) must be adjusted with increase of n until the value of |e n+1 /e n |, from larger than 1 to less than 1, is changed. Consequently, |e n+1 | finally becomes less than the threshold value ∆ of the error-transfinite decision device.

For harmonically excited systems with impacts, one can usually take γ = 1 so as to make the external excitation frequency equal the vibrating frequency of the original system. The practical move to fulfill the control is stated in the ensuing paragraphs.

(1) the initial value B 0 = 0.

(2) if |e n+1 | < ∆ at n + 1 time, no external driving force is applied due to e = 0, B n+1 = B 0 = 0 and u = 0 in Fig. 2 so that the next decision returns to the (2) cycle again; if |e n+1 | ≥ ∆, the decision e = e n+1 brings out an external driving force, and diverts to the (3) process.

(3) u is created by the adjustment algorithm of external driving force and expressed by

B n+1 = B n + α|e n+1 |, u = B n+1 sin(γ ωt), ( 22
)
where α is a weight factor, and 0 < α ≤ 1.

u = B * sin(γ ωt), (23)
This final adjustment brings into effect that the chaos are suppressed by applying a sine periodic force with the amplitude B * to the system.

Suppressing chaos in the vibro-impact systems using period pulse

Consider a general nonlinear mapping system

X n+1 = F(X n ), (24)
where X ∈ R m .

Starting from an initial point X 0 , after pth iteration, the system (24) can be written as

X = F (P) (X 0 ). (25)
The convergence property after pth iteration from the initial points can be known from the Jacobian matrix of the map [START_REF] Pavlovskaia | Modeling of ground moling dynamics by an impact oscillator with a frictional slider[END_REF]. The stability condition of periodic P orbit is all moduli of eigenvalues of Jacobian matrix are less than 1. Assuming that the system exits in stable period-P orbit, from Eq. ( 25), one knows that X P = X 0 . If the variation of the system parameters causes bifurcation or chaos of the system, the map will deviate from the original stable fixed points. That is X P = X 0 .

The idea of the control method is to present the impulse signal Q to the system at the P period point [START_REF] Zhao | Random Fatigue and Clearance Nonlinearity of Mechanical Vibrating Systems[END_REF],

Q = X P -X 0 . (26)
Thus, the orbit, which starts from an initial point X 0 , after pth iteration combined with impulse signal Q, will be

F (P) (X 0 ) -Q = X 0 . (27)
When X P = X 0 , Q = 0 there is no control. When X P = X 0 , Q = 0 will counteract departure. Then the system maps continuously according to Eq. ( 25). Repeating the above process, it is equal to presenting a period pulse to the system.The period pulse can be written by

Q n = k=0 Qδ(n -k P). (28)
k = 0 is the moment when the control is applied. Impulse signal Q makes the orbit end-to-end, thus forms disturbed periodic orbit X 0 , X 1 + ∆X 1 , X 2 + ∆X 2 , . . . , X P-1 + ∆X P-1 and suppresses the chaos eventually. The purpose of this method is not to make the control system return to the original stable periodic orbit but to make it to be in periodic motion all the time so as to suppress the chaos.

As for the vibro-impact system with a clearance, we take one reciprocating impaction as a period, i.e., P = 1. When the system is stable, it satisfies

X 0 = F(X 0 ). (29)
The schematic of the method of periodic impulse suppressing chaos is illustrated in Fig.

3 [START_REF] Zhao | Random Fatigue and Clearance Nonlinearity of Mechanical Vibrating Systems[END_REF]. X i denotes the input at i moment and X 0 denotes the stable output vector (control target). E denotes noise from measurement error, time delay error, and other factors. X i is detected in real time. If X i = X 0 , then Q = 0, no control is applied; if X i = X 0 , then Q = ∆X i , where ∆X i is the derivation of X i from X 0 and the error ∆X i is used as the impulse control signal.

Suppressing chaos of the vibro-impact systems based on delay feedback

The controlling method is to skillfully use certain parts of output signals of the system, through time delaying and subtracting it with original output signals, to get the controlling information to be fed back to the system.

The time delay feedback can be given as the following dynamic forms

X (n + 1) = F(X (n)) + G(X (n -p) -X (n)), (30)
where,

X ∈ R m , G = Diag[g 1 , g 2 , . . . , g i , . . . , g m ], g 1 = g 2 = • • • = g i (i ≤ m)
are the feedback ratio factors. The Fig. 4 shows the schematic of principle of the time delay feedback controlling. When G = 0, the original system exhibits chaotic motion. X (n) and X (np) are the measurements of X at the time n and the time np, respectively. The period of certain unstable period motion can be chosen as the delay time p. When the feedback matrix G has been properly chosen, the system would be driven to the expected orbit.

Dynamic analysis of the vibratory system with a clearance

Numerical analyses are carried out for determining the dynamical responses of the vibro-impact system. According to the non-dimensional parameters and variables (1), we have the ranges of these parameters µ m , µ k and µ c as µ m ∈ (0, 1), µ k ∈ (0, 1) and µ c ∈ (0, 1). The mid-values of parameters µ m , µ k and µ c are considered. The system parameters: µ m = 0.5, µ k = 0.5, µ c = 0.5, ζ = 0.1, r = 0.8 and δ = 0.1-0.6 are chosen for analyses. Clearance is the most important system parameter for the nonlinear system shown in Fig. 1. The influence of parameters ω and δ on symmetrical double-impact periodic motion and bifurcations is depicted by the bifurcation diagram shown in Fig. 5. Representative bifurcation diagrams of before-impact velocity and the corresponding Lyapunov exponents as a function of the excitation frequency ω are shown for δ = 0.1, 0.3, 0.545 and 0.6 in Fig. 6 respectively, which illustrate the bifurcation characteristics of 1-1-1 symmetrical motion. The system exhibits stable 1-1-1 symmetrical motion in the region R 1 of Fig. 5. It undergoes a supercritical Neimark-Sacker bifurcation on the line L 1 , giving rise to a family of attracting invariant circles associated with the fixed point of 1-1-1 symmetrical motion, and a subcritical Neimark-Sacker bifurcation on the line L 2 . There exist quasi-periodic impact motions associated with 1-1-1 symmetrical one, long-period multi-impact motions and chaotic motions in the region R 7 . Crossing the curve L 3 from the region R 1 to the region R 2 , pitchfork bifurcation of 1-1-1 symmetrical motion occurs, and 1-1-1 unsymmetrical motion stabilizes in the region R 2 . Crossing the curve L 4 from the region R 2 to the region R 3 , Neimark-Sacker bifurcation generates a family of invariant circles associated with the fixed point of 1-1-1 unsymmetrical motion. In the region R 3 , the system exhibits quasi-periodic impact motions associated with 1-1-1 unsymmetrical one, long-period multi-impact motions and chaotic motions, etc. Crossing the curve L 5 from the region R 3 to the region R 4 stable 1-2-2 symmetrical motion is generated by a degeneration of chaos. Crossing the curve L 6 from the region R 4 to the region R 5 , pitchfork bifurcation of 1-2-2 symmetrical motion occurs, and 1-2-2 unsymmetrical motion stabilizes. In the region R 5 , the 1-2-2 unsymmetrical motion is of two antisymmetrical forms due to different initial conditions and symmetry of the structure and excitation point. Crossing the curve L 7 , period-doubling bifurcation of 1-2-2 unsymmetrical motion occurs. Thereafter, the motion undergoes a succession of period-doubling bifurcations with decrease in ω, which eventually result in apparently non-periodic or chaotic motions. For smaller values of clearance δ, 1-1-1 symmetrical motion usually undergoes pitchfork bifurcation, and smaller values of clearance generally lead to the enlarged areas of chaotic motions, see Figs. 5 and6(a). For a relatively narrow range of clearance value δ ∈ [0.5352, 0.5602], 1-1-1 symmetrical motion undergoes supercritical Neimark-Sacker bifurcation, advise subcritical Neimark-Sacker bifurcation and pitchfork bifurcation in turn with decrease in the forcing frequency ω; see Figs. 5 and6(c). For large values of clearance δ, 1-1-1 symmetrical motion undergoes supercritical Neimark-Sacker bifurcation with decrease in the forcing frequency ω; see Figs. 5 and6(d). However, for larger values of δ the mass M 1 cannot hit the stops, and the system will undergo simple oscillations and behave as a linear system.

The detailed simulative results are presented in the form of phase plane portraits, time series and projected portraits of Poincaré map for the vibratory system with a clearance δ = 0.1. By looking at the bifurcation diagram in Fig. 6(a), we have a large window of 1-1-1 symmetrical motion for ω in the interval ω ∈ [3.6907, 5]. With the parameter ω reducing, pitchfork bifurcation of 1-1-1 symmetrical motion occurs at ω = 3.6907. Consequently, 1-1-1 symmetrical motion becomes 1-1-1 unsymmetrical one which is of two antisymmetrical forms due to different initial conditions and symmetrical stops. Thereafter, the unsymmetrical double-impact periodic motions undergo Neimark-Sacker bifurcation with further decrease in ω. We give an example of a 1-1-1 symmetrical motion which exists at a forcing frequency of ω = 3.7. A phase portrait of this periodic motion is shown in Fig. 7(a), time series in Fig. 7(b). A phase portrait of 1-1-1 unsymmetrical motions is shown for ω = 3.5 in Fig. 8(a), in which two 1-1-1 antisymmetrical As the forcing frequency ω passes through ω = 2.93193 decreasingly, a grazing contact instability of 2-2-2 motion occurs as seen in Fig. 8(c). Consequently, one impact in the motion period will increase and the 2-2-2 motion transits 6(a), respectively. 1-2-2 symmetrical motion is generated by a degeneration of chaos with decrease in ω. We give an example of a 1-2-2 symmetrical motion which exists at a forcing frequency of ω = 1.72. A phase portrait of this periodic motion is shown in Fig. 10(a), time series in Fig. 10(b). Pitchfork bifurcation of 1-2-2 symmetrical motion occurs as ω passes through ω = 1.7031 decreasingly. Consequently, 1-2-2 unsymmetrical motion is created. A phase portrait of 1-2-2 unsymmetrical motion is shown for ω = 1.65 in Fig. 11(a). With further decrease in ω, the system falls into chaotic-impact motion via period-doubling sequences of 1-2-2 unsymmetrical motion. Fig. 11 shows an important sequence of bifurcations, where the system bifurcates from period one (Fig. 11(a)) to period two (Fig. 11(b)), then from period two to period four (Fig. 11(c)). The chattering-impact motions are found to exist in the region of low forcing frequency. With decrease in ω, non-periodic chattering-impact vibration transits to periodic one. Thereafter, sliding bifurcation of periodic chattering-impact vibration occurs. Consequently, periodic chattering-impact vibration with sticking is created. These imply such motions are very complex, varying from nonsticking types of chattering-impact vibration to sticking types of chattering-impact motion. Time series of non-sticking chattering-impact solutions are shown for ω = 0.65 and 0.6 in Fig. 12(a) and (b), respectively. A time series of sticking chattering-impact vibration is shown for ω = 0.57 in Fig. 12(c), a local zooming of Fig. 12(c) is shown in Fig. 12(c1). This sticking chattering-impact vibration is periodic, which includes regions of free flight, chattering and sticking for the mass M 1 . Sticking in impact systems is mathematically analogous to sliding in electrical systems [START_REF] Di Benardo | Grazing, skipping and sliding: analysis of the non-smooth dynamics of the dc/dc buck converter[END_REF]. The sliding orbits in these electrical systems have been shown to exhibit particular types of sliding bifurcations under parameter variation [START_REF] Bernardo | Self-oscillations and sliding in relay feedback systems: Symmetry and bifurcations[END_REF]. The sliding dynamics on the separation boundary is discussed based on the set-valued vector field theory and the switching bifurcation conditions from the passable boundary to the non-passable boundary are developed in Ref. [START_REF] Albert | On flow switching bifurcations in discontinuous dynamical systems[END_REF]. The normal vector product field function is introduced to determine the switching bifurcation and sliding fragmentation in [START_REF] Albert | On flow switching bifurcations in discontinuous dynamical systems[END_REF]. It is important to note that chattering-impact vibration enlarges the adverse effects such as high noise levels, wear and tear and so on. For δ = 0.545, 1-1-1 symmetrical motion undergoes supercritical Neimark-Sacker bifurcation, advise subcritical Neimark-Sacker bifurcation and pitchfork bifurcation in turn with decrease in the forcing frequency ω; see Fig. 6(c). The sticking and non-sticking types of chattering-impact motions are also found to occur in the region of low forcing frequency; see Fig. 13(a)-(c). It is observed obviously, from Fig. 6(d), that period-one double-impact symmetrical motion of the vibro-impact system undergoes Neimark-Sacker bifurcation, and chaotic-impact motion can occur via phase locking with decrease in the forcing frequency ω. The phase plane portrait and projected Poincaré section near the value of Neimark-Sacker bifurcation are shown in Fig. 14(a) and (b), respectively. The phase plane portrait of period-one double-impact symmetrical motion is plotted for ω = 2.582 in Fig. 14(a), the periodic-impact motion is represented by a stable fixed point in the projected Poincaré section. The attracting invariant circle associated with the fixed point of 1-1-1 symmetrical motion, in the projected Poincaré section, is shown for ω = 2.57 in Fig. 14(b).

8. Suppressing chaotic impacts of the vibratory system with a clearance 8.1. Suppressing chaotic impacts by adjusting the driving force Suppressing chaotic impacts by using sine periodic force is stated in Section 6. The system with non-dimensional parameters: µ m = 0.5, µ k = 0.5, µ c = 0.5, ζ = 0.1 and r = 0.8 are still chosen for analysis. The method effectiveness is analyzed by the presentation of example for the impact oscillator with a clearance δ = 0.1. Bifurcation diagram of before-impact velocity and Lyapunov exponents as a function of the excitation frequency ω are shown for δ = 0.1 in Fig. 6(a) and (a1), respectively. The detailed analysis is presented, which has been illustrated by Figs. 7-12. The system exhibits 1-1-1 symmetrical motion at ω = 3.7 and chaotic-impact motion at ω = 2.83, as seen in Figs. 7(a), 8(d) and 9(c), respectively. The fixed point X 0 = ( ẋ1+ , x 2 , ẋ2 , τ 0 ) = (0.34545790, -0.00634090, 0.02116031, 1.80213070) T of 1-1-1 symmetrical motion at ω = 3.7 is taken as the desired fixed point. The impact period and the time interval between consecutive impacts occurring at the same stop are easy to measure in practical engineering. Consequently, the angle shift is simple to determine. Based on the consideration, we choose the variable τ 0 (phase shift) in the Poincaré section as the desired target. We take the phase shift τ as the practical output, and α = 0.327, ∆ = 10 -5 . Fig. 15 shows the simulative results after applying the control on the chaotic-impact motion which exists at a forcing frequency of ω = 2.83. In order to observe the controlling process clearly, the control acts on the system when n = 2000. Fig. 15(a) shows time response of the amplitude B of the external sine driving signal adjusted by the adjustment algorithm. The initial value of B is zero. After applying the control on the system, we find that the amplitude B swiftly tends to a stable value B * = 0.152618. The controlled system is in a stable state with an external driving force u = 0.152618 sin(ωt) acted on the mass M 1 . The time response of phase shift τ is shown in Fig. 15(b). The time responses of variables ẋ1+ , x 2 and ẋ2 are presented in Fig. 15

Suppressing chaotic impacts by using period pulse

The effect on the periodic impulse signal suppressing the chaotic-impact motion is observed by the presentation of an example. Consider the ideal case of E = 0 in Fig. 3 and the stable fixed point as the control target. Once the state signals deviate from the control target, the system dynamically produces period impulse signals to suppress chaos and bifurcation.

From Fig. 6(a), when no control is applied, the system falls into chaos via Neimark-Sacker bifurcation and phase locking with decrease in the excitation frequency ω.

The system parameters: µ m = 0.5, µ k = 0.5, µ c = 0.5, ζ = 0.1, r = 0.8 and δ = 0.1 are considered. The excitation frequency ω is taken as the bifurcation parameter. As we have known, the system is in stable 1-1-1 symmetric motion for ω ∈ [3.6907, 5]. Hence the projection of the motion to Poincaré section σ is a fixed point. We have given an example of a 1-1-1 symmetrical motion existing at a forcing frequency of ω = 3.7. A phase portrait of this periodic motion is shown in Fig. 7(a), time series in Fig. 7(b). When ω is decreased gradually, pitchfork bifurcation of 1-1-1 symmetrical motion and Neimark-Sacker bifurcation of 1-1-1 unsymmetrical motion occur in turn. Fig. 6(a) and 9 show the simulative result that the system varies from stable fixed point to Neimark-Sacker bifurcation and even chaos.

When ω = 2.83, the system exhibits the chaotic-impact motion, as seen in Figs. 8(d) and 9(c). Fig. 17 shows the simulative results after applying the control on the above chaotic-impact motion while ω = 2.83. In order to observe the controlling process clearly, the control acts on the system when n = 2000. A phase portrait of the system, after applying the control, is shown in Fig. 18(a), time series in Fig. 18(b). The system with ω = 2.83, after applying the control, is found to exhibit 1-1-1 unsymmetrical motion. Comparing Fig. 8(d) with Fig. 18, owing to the control, the original chaotic-impact motion is suppressed. The method needs little time from control action to chaos suppressing and the time response is very fast within a few impact periods, as seen in local zooming of the time response of phase shift τ shown in Fig. 17 with the fixed point of 1-2-1 motion. The simulative result also shows that if the g is too small, the system may converge to period-two fixed points. Fig. 23 shows the system variable's time responses and its phase plane portrait when g = 0.1 and other conditions are unchanged. From Fig. 23, one can observe that the system at ω = 2.83, after controlling with g = 0.1, exhibits 2-4-2 motion.

Conclusions

A two-degree-of-freedom oscillator with a clearance is considered. The maximum displacement of one of the masses is limited to a threshold value by the symmetrical stops. Impacts between the mass and the stops are described by an instantaneous coefficient of restitution. Dynamics of the system is studied with special attention to periodic-impact motions and bifurcations. Period-one double-impact symmetrical motion and transcendental impact map of the system is derived analytically. Stability and local bifurcations of the period-one double-impact symmetrical motions are analyzed by using the impact Poincaré maps. The Lyapunov exponents in vibratory systems with impacts are calculated by using the transcendental impact map. Neimark-Sacker bifurcations of period-one double-impact motions are found to occur in the system. It is well known that no Neimark-Sacker bifurcation of periodic-impact motion occurs in the single-degree-of-freedom oscillators with single stop or symmetrical stops, e.g., see Refs. [START_REF] Shaw | A periodically forced piecewise linear oscillator[END_REF][START_REF] Shaw | The dynamics of a harmonically excited system having rigid amplitude constraints. Part 1, Part 2[END_REF]. However, not only Neimark-Sacker bifurcation associated with period-one doubleimpact symmetrical motion occurs for such two-degree-of-freedom vibratory system with a clearance, but also Neimark-Sacker bifurcation of period-one double-impact unsymmetrical motions does. The influence of the clearance δ and the excitation frequency ω on symmetrical double-impact periodic motion and bifurcations is analyzed. A series of other periodic-impact motions are found and the corresponding bifurcations are analyzed. For smaller values of clearance, period-one double-impact symmetrical motion usually undergoes pitchfork bifurcation with decrease in the forcing frequency. For large values of the clearance, period-one double-impact symmetrical motion undergoes Neimark-Sacker bifurcation with decrease in the forcing frequency. The chattering-impact vibration and the sticking phenomena are found to occur in the region of low forcing frequency, which enlarges the adverse effects such as high noise levels, wear and tear and so on. With decrease in ω, non-periodic chattering-impact vibration transits to periodic one. Thereafter, sliding bifurcation of periodic chattering-impact vibration occurs. Consequently, periodic chattering-impact vibration with sticking is created. These imply that the dynamics of this system is very complex, varying from different types of periodic-impact motions to chaos, even chattering-impacting vibration and sticking, etc.

Impacting vibration arises very naturally in the vibratory systems with clearances. However, chaotic or chattering impacts are undesirable as they bring about the adverse effects such as high noise levels, wear and tear and so on. Based on the impact Poincaré map, we consider some methods of chaos suppression. The state of the vibro-impact system, just immediately after impact, is chosen as the Poincaré section, and the fixed point of 1-1-1 motion is chosen as the desired target. These methods are applied to the two-degree-of-freedom vibratory system with a clearance. Chaotic-impact motions of the system can be suppressed by using sine periodic force, delay feedback and feedbackbased method of period pulse to enlarge the profitable effects and minimize the adverse ones. The results show that the methods can suppress chaotic-impact motions effectively.

Fig. 1 .

 1 Fig. 1. Schematic of a harmonically excited system having symmetrical rigid stops.

Fig. 2 .

 2 Fig.2. Schematic of the method for suppressing chaos by using sine periodic force.

( 4 )

 4 After the driving force is applied, if |e n+1 | ≥ ∆, the adjustment direction of B is modified according to the value of |e n+1 /e n |. If |e n+1 /e n | > 1, B is adjusted in the contrary direction so that let α = -α. Whereas, α is fixed and the (3) process is continued. If new |e n+1 | < ∆, the error-transfinite decision device brings out the output of e = 0, the coefficient adjustment of driving force ends. The final coefficient is expressed by B * , and the external driving force is given by

Fig. 3 .Fig. 4 .

 34 Fig. 3. Schematic of the method of periodic impulse suppressing chaos.

  (a)-(d) and (a1)-(d1),

Fig. 5 .

 5 Fig. 5. Bifurcation diagram associated with parameters (ω, δ), system parameters: µ m = 0.5, µ k = 0.5, µ c = 0.5, ζ = 0.1, r = 0.8.

Fig. 6 .

 6 Fig. 6. Bifurcation diagrams of before-impact velocity and corresponding Lyapunov exponents as a function of the excitation frequency ω: (a) and (a1) δ = 0.1; (b) and (b1) δ = 0.3; (c) and (c1) δ = 0.545; (d) and (d1) δ = 0.6.

Fig. 7 .

 7 Fig. 7. Phase plane portrait and time series: 1-1-1 symmetrical motion, δ = 0.1, ω = 3.7.

Fig. 8 .

 8 Fig. 8. Phase plane portraits: (a) 1-1-1 unsymmetrical motion, δ = 0.1, ω = 3.5; (b) 2-2-2 motion, δ = 0.1, ω = 2.95; (c) 2-2-2 motion with grazing a stop, δ = 0.1, ω = 2.93193; (d) chaotic motion, δ = 0.1, ω = 2.83.

Fig. 9 .

 9 Fig. 9. Projected Poincaré sections (δ = 0.1): (a) quasi-periodic attractor associated with the fixed point of 1-1-1 unsymmetrical motion, ω = 3.03; (b) fixed points of 2-2-2 motion, ω = 2.95; (c) chaos, ω = 2.83.

Fig. 10 .

 10 Fig. 10. Phase plane portrait and time series: 1-2-2 symmetrical motion, δ = 0.1, ω = 1.72.

Fig. 11 .

 11 Fig. 11. Phase plane portraits (δ = 0.1): (a) 1-2-2 unsymmetrical motion, ω = 1.65; (b) 2-4-4 unsymmetrical motion, ω = 1.63; (c) 4-8-8 unsymmetrical motion, ω = 1.628. into chaotic motion via grazing bifurcation; see Figs. 8(d) and 9(c). Looking at the bifurcation diagram in Fig. 6(a), we have a large window of 1-1-1 symmetrical motion for ω in the interval ω ∈ [3.6907, 5], a window of 1-1-1 unsymmetrical motion for ω ∈ [3.0738, 3.6907]. Two windows of 1-2-2 symmetric and unsymmetric motions are also observed from Fig. 6(a), respectively. 1-2-2 symmetrical motion is generated by a degeneration of chaos with decrease in ω. We give an example of a 1-2-2 symmetrical motion which exists at a forcing frequency of ω = 1.72. A phase portrait of this periodic motion is shown in Fig. 10(a), time series in Fig. 10(b). Pitchfork bifurcation of 1-2-2 symmetrical motion occurs as ω passes through ω = 1.7031 decreasingly. Consequently, 1-2-2 unsymmetrical motion is created. A phase portrait of 1-2-2 unsymmetrical motion is shown for ω = 1.65 in Fig. 11(a). With further decrease in ω, the system falls into chaotic-impact motion via period-doubling sequences of 1-2-2 unsymmetrical motion.

Fig. 12 .

 12 Fig. 12. Time series with chattering-impact vibration (δ = 0.1): (a) non-periodic chattering-impact vibration without sticking, ω = 0.65; (b) periodic chattering-impact vibration without sticking, ω = 0.6; (c) chattering-impact vibration with sticking, ω = 0.57; (c1) local zooming of figure (c).

Fig. 13 .

 13 Fig. 13. Time series with chattering-impact vibration (δ = 0.545): (a) chattering-impact vibration without sticking, ω = 0.35; (b) chattering-impact vibration without sticking, ω = 0.3; (c) chattering-impact vibration with sticking, ω = 0.25.

Fig. 14 .

 14 Fig. 14. Phase plane portrait and projected Poincaré section: (a) phase plane portrait of 1-1-1 symmetrical motion, δ = 0.6, ω = 2.582; (b) attracting invariant circle associated with the fixed point of 1-1-1 symmetrical motion, δ = 0.6, ω = 2.57.

Fig. 15 .

 15 Fig. 15. Time responses of some variables (no controlling (n < 2000), after controlling (n ≥ 2000)), ω = 2.83.

Fig. 16 .

 16 Fig. 16. Phase plane portrait and time series after controlling, ω = 2.83.

  (d)-(e), respectively. From the simulative results shown in Fig. 15(b)-(e), we can observe that τ finally trends to τ 0 . However, the other variables ẋ1+ , x 2 and ẋ2 in the Poincaré section, stabilize to new values, i.e., X * 0 = ( ẋ * 1+ , x * 2 , ẋ * 2 , τ * 0 ) T = (-0.419595312, -0.02899704, 0.02903741, 1.80213090) T . A phase portrait of the system, after controlling, is shown in Fig. 16(a), time series in Fig. 16(b). The system at ω = 2.83, after controlling, exhibits 1-2-1 motion. From Figs. 15 and 16, the control can suppress chaos effectively. The method needs little time from control action to chaos suppressing and the time response is very fast within forty-seven impact periods. Local zooming of the time response of phase shift τ is shown in Fig. 15(b1).

Fig. 17 .

 17 Fig. 17. Time responses of some variables (no controlling (n < 2000), after controlling (n ≥ 2000)), ω = 2.83.

  (a1).

Fig. 18 .

 18 Fig. 18. Phase plane portrait and time series after controlling, ω = 2.83.

Fig. 19 .

 19 Fig. 19. Time responses of some variables (no controlling (n < 2000), after controlling (n ≥ 2000)), ω = 2.83.

Fig. 20 .

 20 Fig. 20. Phase plane portrait and time series after controlling, ω = 2.83.

Fig. 21 .

 21 Fig. 21. The simulative results of the time responses of the system transitional process with g: (a) g = 0.8; (b) g = 0.6; (c) g = 0.3; (d) g = 0.25.

Fig. 22 .

 22 Fig. 22. The simulative result with g = 0.97: (a) time response of variable τ (no controlling (n < 2000), after controlling (n ≥ 2000)); (b) and (c) projected portraits of Poincaré map (after controlling); (d) phase plane portrait (after controlling).

Fig. 23 .

 23 Fig. 23. The simulative result with g = 0.1: (a) time response of variable τ (no controlling (n < 2000), after controlling (n ≥ 2000)); (b) phase plane portrait (2-4-2 motion).

Acknowledgements

The authors gratefully acknowledge the support by the National Natural Science Foundation (10572055, 50475109), the Natural Science Foundation of Gansu Province Government of China (3ZS061-A25-043 (key item)) and the 'Qing Lan' Talent Engineering by the Lanzhou Jiaotong University.

Suppressing chaotic impacts based on delay feedback

In the section we study how to suppress chaotic-impact vibration in the vibratory system with clearances when control target cannot be known beforehand. Some state variables of the previous impacts are used to suppress chaotic motion by negative feedback. The purpose of controlling is that when the chaotic-impact motion happens, it continuously feeds back some of the time-delayed state variables until the chaos is suppressed. When the chaos is suppressed, the system will be converged to the fixed point on the Poincaré section. The effect on the method is illustrated with an example of the vibratory system with a clearance.

Set the delay time as p = 1. The non-dimensional parameters: µ m = 0.5, µ k = 0.5, µ c = 0.5, ζ = 0.1, r = 0.8 and δ = 0.1 are still chosen for observing the effect on suppressing chaos based on delay feedback. With changes of the bifurcation parameters ω, the system will deviate from stable periodic orbit and convert to chaotic-impact motion. As we have known, Fig. 6(a) and (a1) respectively show bifurcation diagram of before-impact velocity and Lyapunov exponents as a function of the excitation frequency ω. A chaotic attractor existing at ω = 2.83 is shown in Fig. 9(c).

Fig. 19 shows the suppressing of the above chaotic-impact motion, while G = Diag[g 1 , g 2 , g 3 , g 4 ], g 1 = g 2 = g 3 = g 4 = 0.5 and p = 1. To illustrate the controlling effect clearly, we put the controlling at n = 2000. The chaotic attractor in Fig. 9(c) is suppressed to the fixed point X * 0 = ( ẋ * 1+ , x * 2 , ẋ * 2 , τ * 0 ) T = (-0.41959444, -0.02899695, 0.02903516, 1.80121257) T by the delay feedback control. A phase portrait of the system, after controlling, is shown in Fig. 20(a), time series in Fig. 20(b). The system at ω = 2.83, after controlling, is found to exhibit 1-2-1 motion. It is shown that the time delay feedback controlling method is effective in suppressing the above chaotic-impact motion.

To understand the effect of this method on suppressing such chaotic-impact motions roundly, through numerical simulations, we analyze the relationship between the feedback matrix G = Diag[g 1 , g 2 , . . . , g i , . . .] and the result of controlling, as well as that between the time of controlling and its effect. We discuss the relationship between feedback matrix G = Diag[g 1 , g 2 , . . . , g i , . . .] and its effect at first. While setting G = Diag[g 1 , g 2 , g 3 , g 4 ] and g 1 = g 2 = g 3 = g 4 = g in the region of (0, 1), we decrease g from 0.98 gradually. It is seen that either large or small g will prolong the converging time. Fig. 21 is the simulative results of the responses of the system's transitional process with different g while other conditions are unchanged. The case of phase shift τ is only illustrated in Fig. 21, and other variables are the same. By looking at the time responses of transitional process in Figs. 19(a1) and 21, one can find that the time delay feedback controlling with g = 0.3-0.5 is more effective in suppressing the above chaotic-impact motion.

The simulative result shows that if the g is too large, the system may converge to an attracting invariant circle, rather than converge to a fixed point. For this system, if g > 0.9677, it cannot be converged to the fixed point. Fig. 22 shows the system variable's time responses and its projection on the Poincaré section σ when g = 0.97 and other conditions are unchanged. To clearly show the results on the Poincaré section σ , the first 2000 points at the beginning of the n-axis are not shown in the Poincaré section. It is obviously observed that the system converge to an attracting invariant circle. Comparing with projected Poincaré section (Fig. 22(b) and (c)) and phase portrait (Fig. 22(d)), one can find that the attracting invariant circle on impact Poincaré section represents a quasi-periodic attractor associated Appendix Substituting t = t n+1 and the formulae ( 8)-( 11) into the general solutions ( 6) and ( 7) of Eq. ( 2) between the n and (n + 1)th impacts, we can obtain the state variables X n+1 = ( ẋ1(n+1)+ , x 2(n+1) , ẋ2(n+1) , τ n+1 ) T just after the (n + 1)th impact occurring at the stop A or Ā, which are given by Eqs. ( 12)-( 15).

where x 1(n+1) = δ or -δ. By means of Eq. (A.1), we can define a function h( ẋ1n+ , x 2n , ẋ2n , τ n , t n+1 ) as

Supposing ∂h/∂t n+1 = 0, according to the implicit function theorem, we can solve Eq. (A.2) for

The partial derivative of t n+1 with respect to ẋ1n+ , x 2n , ẋ2n and τ n can be computed by

Inserting the formulae ( 8)- [START_REF] Luo | Codimension two bifurcations of fixed points in a class of vibratory systems with rigid stops[END_REF] and Eq. (A.3) into Eqs. ( 12)-( 15), we get the transcendental impact map derived from the equations of motion, between impacts, supplemented by transition conditions at the instants of impacts, which is expressed briefly by

The components in the matrix ( 16) are obtained by

Although the components in the matrix ( 16) are derived, the state variables ẋ1n+ , x 2n , ẋ2n and τ n , just after the nth impact, are still included in all components in the matrix. It is important to note, that the state variables just after the nth impact cannot be analytically derived by combining Eqs. ( 6) and [START_REF] Dimentberg | Stochastic and/or chaotic response of a vibration system to imperfectly periodic sinusoidal excitation[END_REF] and the matching conditions at the impact instants due to the transcendental equations. Due to the unknown state variables, we cannot directly use the Jacobian matrix to compute its eigenvalues. However, it is easy to obtain the state variables ẋ1n+ , x 2n , ẋ2n and τ n just after the nth impact numerically. Inserting the state variables ẋ1n+ , x 2n , ẋ2n and τ n numerically obtained into the components of the Jacobian matrix [START_REF] Foale | Dynamical complexities of forced impacting systems[END_REF], one can compute the eigenvalues of the matrix.