Nicolas Van Goethem
email: vangoeth@fc.ul.pt

DIRECT EXPRESSION OF INCOMPATIBILITY IN CURVILINEAR SYSTEMS

Keywords: 2010 Mathematics Subject Classification. 35J48, 35J58, 53A05, 74A45, 74B99 incompatibility, dislocations, elasticity, differential geometry, curvilinear basis, spherical system

With this note, we would like present a method to compute the incompatibility operator in any system of curvilinear coordinates/components. The procedure is independent of the metric in the sense that the expression can be expressed by means of the basis vectors only, which are first defined as normal and tangent to the domain boundary and then extended to the whole domain. It is in some sense an intrinsic method, since the chosen curvilinear system depends solely on the geometry of the domain boundary. As an application, the in extenso expression of incompatibility in a spherical system is given.

Introduction

Let us consider a smooth body Ω ⊂ R 3 . The incompatibility is a well-known operator in elasticity, since as applied to the linearized strain tensor , it determines whether the strain derives from a displacement field. Specifically, let the elastic strain be obtained by a constitutive law from the stress tensor σ, i.e., = Cσ, with C the compliance fourth-rank tensor, then inc = 0 if and only if = ∇ S u for some displacement field u.

If on the contrary it turns out not to vanish, then Kröner's works [START_REF] Kröner | Continuum theory of defects[END_REF] tell us that dislocations are present, preventing the existence of a well-defined displacement field defined in the whole body. In general a dislocation is a three-dimensional line singularity for the strain field, reducing to a straight line in some simplified cases where a 2D treatment is sufficient, cf., e.g., [START_REF] Van Goethem | A distributional approach to 2D Volterra dislocations at the continuum scale[END_REF]). Specifically, Kröner's relation reads

Curl κ = inc , (1.1)
where the contortion tensor κ is related to the tensor-valued dislocation density Λ by κ = Λ -1 2 tr ΛI 2 . At the mesoscopic scale the dislocation density reads Λ = Λ L = τ ⊗ bH 1 L , where H 1 L stands for the one-dimensional Hausdorff measure concentrated in the dislocation loop L. At the mesoscale, Kröner's relation also holds, as proved in [START_REF] Van Goethem | Fields of bounded deformation for mesoscopic dislocations[END_REF][START_REF] Van Goethem | Incompatibility-governed singularities in linear elasticity with dislocations[END_REF]. At the macro (or continuous) scale (which is the scale considered in the present work), Λ is a smooth tensor obtained from its mesoscopic counterpart by some regularization. The fact that at the mesoscale dislocations are closed loops or end at the boundary implies that div Λ = 0. However div κ = 0 except in particular cases, for instance if one considers pure edge dislocation loops in 3D, i.e., satisfying tr Λ = 0, and therefore the knowledge of the right-hand side of (1.1) is in general not sufficient to uniquely determine the field κ. Note that in this case, The Frank tensor Curl t and the dislocation density are univoquely related, since (1.1) reduces to Curl t = κ in Ω with ×N = 0 on ∂Ω, by virtue of (1.4) and a uniqueness result as proved in [START_REF] Scala | Dislocations at the continuum scale: functional setting and variational properties[END_REF].

Being a symmetric tensor, the elastic strain satisfies by Beltrami decomposition [START_REF] Maggiani | A compatible-incompatible decomposition of symmetric tensors in L p with application to elasticity[END_REF],

= ∇ S u + 0 (1.2)
with u a vector field and where 0 = inc F represents in a Cartesian system the incompatible part of the strain, for some symmetric and solenoidal second-rank tensor F . Thus, the field F , related to the presence of dislocations satisfies by (1.1) and (1.3) an equation of the following form:

inc inc F = Curl κ, (1.3)
which is proven in [START_REF] Amstutz | Analysis of the incompatibility operator and application in intrinsic elasticity with dislocations[END_REF] to be well posed (with appropriate boundary conditions on F and Curl t F × N ), provided the dislocation density is given (here we assume that κ is known, for instance by solving a transport-reactiondiffusion equation, as done for point-defects in [START_REF] Van Goethem | Dynamic prediction of point defects in Czochralski silicon growth. An attempt to reconcile experimental defect diffusion coefficients with the V /G criterion[END_REF]).

The incompatibility operator on second-rank tensors is classically defined as

inc T := Curl Curl t T, (1.4)
meaning (in a Cartesian system) that the curl is taken over the rows and the columns of a second-rank tensor T , consecutively. In the present work, our concern is to compute the incompatibility in a subset of Ω, say an inclusion, whose shape might be arbitrary, or solution of a geometric optimization problem. For this reason, there is a need to express incompatibility in curvilinear systems, chosen to fit the inclusion geometry.

In a general curvilinear system, the same definition (1.4) holds, but care must be taken, since the covariant derivatives do not commute in general because the basis vectors depend on the position, and hence must also be differentiated in (1.4). We will write the general second-rank tensor as T in the Cartesian basis as T = T CART ij e i ⊗ e j , and in the curvilinear basis as T = T CURV ij g i ⊗ g j . Our point of view is that the basis

{g i } i = {N, τ R } for R = A, B
, where τ R are tangent to the boundary, and N is its unit outwards normal, is first defined on the domain boundary ∂Ω and then extended in Ω, where its differentials can be computed. This latter operation gives rise to 5 numbers: κ R , the 2 surface curvatures, γ R , the 2 divergences of τ R , and ξ, the deviation with respect to the principal directions. From these five numbers, the Christoffel symbols can be found and hence the covariant derivatives, thus the expression of the curvilinear curl, and eventually of the incompatibility. An important preliminary step is to express the differentials of the basis vectors in terms of κ R , γ R , ξ, which is given in Theorem 3.2 whose proof can be found in [START_REF] Amstutz | Analysis of the incompatibility operator and application in intrinsic elasticity with dislocations[END_REF]. The curvilinear coordinates are simply the abscissae of the curves with tangent vector τ R , and the radial coordinate r associated to N .

Expressions of the incompatibility in a general curvilinear system are rarely found in extenso in the literature. Let us mention L. E. Malvern's textbook [START_REF] Malvern | Introduction to the mechanics of a continuous medium[END_REF], where an expression can be found in Appendix II, expressed in terms of the metric factors h i defining the intrinsic metric. Our approach can be considered as an metric-free alternative, since we base our method on the sole geometry of the domain boundary and on the natural orthogonal basis that we may define on it. In Section 5, we apply our method to the spherical system, and provide explict expressions of all six components of inc T .

Applications of this method can be found in dislocation theory, where the incompatibility in an inclusion is to be found, in order to determine its dislocation content, and design an optimization method where inclusions are inserted as to obtain a maximal increase (or decrease) of certain functionals.

Notations and conventions.

Let Ω be a bounded domain of R 3 with smooth boundary ∂Ω. By smooth we mean C ∞ , but this assumption could be considerably weakened. Curl, incompatibility and cross product with 2nd order tensors are defined componentwise as follows with the summation convention on repeated indices.

Here, E represents 2nd rank tensors, N is a vector, and is the Levi-Civita 3rd rank tensor. One has:

( Curl E) ij := (∇ × E) ij = jkm ∂ k E im , ( inc E) ij := ( Curl Curl t E) ij = ikm jln ∂ k ∂ l E mn .
Note that the transpose of Curl E will be denoted by Curl t E. Moreover the tensorial product of two vectors a and b will be denoted by a ⊗ b, while a b = 1 2 (a ⊗ b + b ⊗ a).

Some physical motivations

In this section, we provide two examples of models in which the incompatibility plays a crucial role and must be expressed in a curvilinear system.

2.1. The incompatibility operator in linearized elasticity with dislocations. The strain energy density in small-strain elasticity of an isotropic material reads

W e ( ) = 1 2 A • , (2.1)
where is the linearized elastic strain tensor. The stress tensor is classically defined as σ := ∂We ∂ = A . Furthermore, by the symmetry property of , Beltrami decomposition (1.3) holds. The potential energy is defined as

W( ) = Ω (W e ( ) -f • u -G • F ) dx, (2.2)
which, in the absence of dislocations, i.e., as F = 0, yields by minimization the standard Equilibrium equation,

-div (A ) = -div A∇ S u = f , (2.3)
with f the body force and u the displacement field.

= δW( ) δu , ũ = Ω σ • ∇ S ũ -f • ũ dx, 0 = δW( ) δF , F = Ω σ • inc F -G • F dx, (2.5)
providing after some easy part integrations, the strong forms

-div σ = f, inc σ = G,
which appear clearly as a generalization of (2.3). Recalling (1.3), the complete problem consists in solving the coupled problem with unknowns u and F :

-div A∇ S u = f + div (A inc F ) , inc (A inc F ) = G -inc A∇ S u . (2.6)
Material isotropy yields A = µI 4 + λI 2 ⊗ I 2 and hence (2.8) rewrites as

-div A∇ S u = f + div λ tr 0 I 2 , inc (A inc F ) = G -inc (λ div uI 2 ) . (2.7)
Note that the decoupled problem is found as soon as either λ = 0, or incompressibility is assumed, i.e., tr = tr 0 = div u = 0, and reads

   -div Ã∇ S u = f inc à inc F = G in Ω,
where the special form of à due to incompressibility can be found in [START_REF] Hoger | Linear elasticity for constrained materials: Incompressibility[END_REF].

The point here is to observe that according to the geometry of Ω it is usefull to have an expression of incompatibility in curvilinear coordinates/components systems according to the geometry of Ω.

Dislocation-induced dissipation.

Let us define the system Helmholtz free energy density as

Ψ := W e ( ) + W defect ( Curl κ), (2.8)
whose elastic part is simply the strain energy of previous section, and whose defect part is assumed to depend on the curl of the contortion tensor κ. Therefore, the free energy Ψ is partially of second-order in the sense that the defect internal variable appears in the form of its derivatives (here its curl). For simplicity, let us assume a quadratic law in the higher-order terms, viz., W defect ( Curl κ) = 1 2 M Curl κ • Curl κ, with M a positive-definite fourth-rank tensor. By Kröner's relation (1.1), the energy of an inclusion ω ⊂ Ω reads

W defect := ω 1 2 M inc 0 • inc 0 dx.
(2.9)

Therefore, minimizing this energy will again lead one to evaluate or express the incompatibility in local basis appropriate to the geometry of ω.

Note that a full second-order energy density would read for instance Ψ := W e ( ) + Ŵe ( Curl t , div ) + Ŵdefect (κ, Curl κ, div κ) + Wdefect ( 0 ), where Curl t is recognized as the Frank tensor, i.e. the gradient of the rotation field.

Extension and differentiation of the normal and tangent vectors to a surface

The aim here is to construct a curvilinear basis on the boundary which should be smooth and also orthonormal, starting from the vector N ∂Ω normal to the boundary and defining two tangent vectors perpendicular to N ∂Ω . This basis is then extended to the whole body. The natural moving frame sought is close in spirit to the Darboux frame of surfaces, though in principle the latter may only be defined at non-umbilical points. As a matter of fact, in order to achieve a certain level of generality, we will not consider principal lines of curvature with their associated principal curvatures, and hence the gradient of the normal vector will be given by a symmetric matrix with possibly non-zero extradiagonal components. Detail on this section, and in particular, the proofs, can be found in [START_REF] Amstutz | Analysis of the incompatibility operator and application in intrinsic elasticity with dislocations[END_REF].

3.1. Signed distance function and extended unit normal. We denote by N ∂Ω the outward unit normal to ∂Ω, and by b the signed distance to ∂Ω, i.e.,

b(x) = dist(x, ∂Ω) if x / ∈ Ω, -dist(x, ∂Ω) if x ∈ Ω.
We recall the following results ([3], Chap. 5, Thms 3.1 and 4.3).

Theorem 3.1.

There exists an open neighborhood W of ∂Ω such that

(1) b is smooth in W ;

(2) every x ∈ W admits a unique projection p ∂Ω (x) onto ∂Ω;

(3) this projection satisfies

p ∂Ω (x) = x - 1 2 ∇b 2 (x), x ∈ W ; (3.1) (4) it holds ∇b(x) = N ∂Ω (p ∂Ω (x)), x ∈ W.
In particular, this latter property shows that ∇b(x) = N ∂Ω (x) for all x ∈ ∂Ω and |∇b(x)| = 1 for all x ∈ W . Therefore, we define the extended unit normal by

N (x) := ∇b(x) = N ∂Ω (p ∂Ω (x)), x ∈ W. (3.2)
3.2. Tangent vectors on ∂Ω. For all x ∈ ∂Ω, we denote by T ∂Ω (x) the tangent plane to ∂Ω at x, that is, the orthogonal complement of N ∂Ω (x). As ∂Ω is smooth, there exists a covering of ∂Ω by open balls B 1 , ..., B M of R 3 such that, for each index k, two smooth vector fields τ A ∂Ω , τ B ∂Ω can be constructed on ∂Ω ∩ B k where, for all x ∈ ∂Ω ∩ B k , (τ A ∂Ω (x), τ B ∂Ω (x)) is an orthonormal basis of T ∂Ω (x). In all the sequel, the index k will be implicitly considered as fixed and the restriction to B k will be omitted. In fact, for our needs, global properties and constructions will be easily obtained from local ones through a partition of unity subordinate to the covering.

Using that the Jacobian matrix

DN (x) = D 2 b(x) of N (x) is symmetric, differentiating the equality |N (x)| 2 = 1 entails ∂ N N (x) = DN (x)N (x) = 0, x ∈ W. (3.3)
In other words, N (x) is an eigenvector of DN (x) for the eigenvalue 0. For all x ∈ ∂Ω, the system (τ A ∂Ω (x), τ B ∂Ω (x), N ∂Ω (x)) is an orthonormal basis of R 3 . In this basis, DN (x) takes the form

DN (x) =   κ A ∂Ω (x) ξ ∂Ω (x) 0 ξ ∂Ω (x) κ B ∂Ω (x) 0 0 0 0   , x ∈ ∂Ω, (3.4)
where κ A ∂Ω , κ B ∂Ω and ξ are smooth scalar fields defined on ∂Ω.

If R ∈ {A, B}, we denote by R * the complementary index of R, that is, R * = B if R = A and R * = A if R = B.
d = 1 + b κ A ∂Ω • p ∂Ω 1 + b κ B ∂Ω • p ∂Ω -(b ξ ∂Ω • p ∂Ω ) 2 .
Possibly adjusting W so that d(x) > 0 for all x ∈ W , we define in W , and for R = A, B,

τ R = τ R ∂Ω • p ∂Ω , κ R = d -1 (1 + b κ R * ∂Ω • p ∂Ω )(κ R ∂Ω • p ∂Ω ) -b (ξ ∂Ω • p ∂Ω ) 2 , (3.5)
ξ = d -1 ξ ∂Ω • p ∂Ω , κ = κ A + κ B , γ R = div τ R . (3.6)
Obviously, for each x ∈ W , the triple τ A (x), τ B (x), N (x) forms an orthonormal basis of R 3 . Next, we compute the normal and tangential derivatives of these vectors. We denote the tangential derivative ∂ τ R by ∂ R for simplicity, i.e., ∂ R u := Duτ R , where Du stands for the differential of u, and ∂ R u its value in the direction τ R .

Theorem 3.2 ([1]

). It holds in W :

∂ N τ R = 0, ∂ R N = κ R τ R + ξτ R * , ∂ R τ R = -κ R N -γ R * τ R * , ∂ R * τ R = γ R τ R * -ξN, div N = tr DN = ∆b = κ. (3.7) Corollary 3.3 ([1]). If f is twice differentiable in Ω it holds ∂ R ∂ N f = ∂ N ∂ R f + κ R ∂ R f + ξ∂ R * f. (3.8)

Differential geometry on the boundary with curvinormal basis

At each point x ∈ ∂Ω the curvilinear basis g i (x) i=A,B,N := τ A (x), τ B (x), N ∂Ω (x) is orthonormal and differentiable by Theorem 3.2. Therefore it will be called curvinormal. Remark that indices P, Q, R will stand for A or B, and denote one of the two orthogonal tangent vectors on the boundary, whereas index N will always be associated to the normal N ∂Ω . In some sense, the chosen curvilinear basis is a generalization to general surfaces of the spherical or cylindrical basis. We recall that ∂ i means the differential in the direction g i . Let u be a scalar. Then,

∂ i u = ∂ R u = τ R • ∇u = Duτ R for R = A, B, or ∂ N u = N • ∇u = DuN for i = N , with ∇ = e i ∂ xi
the Cartesian gradient operator, where e i stands for the ith Cartesian basis vector. For instance, the gradient in spherical coordinates reads ∇u = ∂ r ue r + 1 r ∂ φ ue Φ + 1 r sin φ ∂ θ ue θ and hence

∂ A = 1 r ∂ φ and ∂ B = 1 r sin φ ∂ θ .
Recall that partial curvilinear derivatives do not commute, as shown in Corollary 3.3.

4.1.

Christoffel symbols and Riemannian curvature. A general vector field will be written as v = v i g i with v i its covariant components. Moreover, the extrinsic metric is Euclidean, since g ij := g i • g j = δ ij . Let g i := g ij g j be the dual of the basis vector. The second Christoffel symbol Γ p ij is defined as the linear operator such that [START_REF] Ciarlet | An introduction to differential geometry with applications to elasticity[END_REF]

∂ j g p = -Γ p ij g i . (4.1)
In other words, Γ p ij := -g i • ∂ j g p . Note also that since Ω is embedded in a Euclidean space, one has ∂ j g i = ∂ j g i . Connection. As a consequence, for vector v it holds

∂ j v = ∂ j (v i g i ) = (∂ j v i -Γ p ij v p )g i = v i j g i , (4.2)
where the covariant derivative of the covariant component of v reads

v i j := ∂ j v i -Γ p ij v p . (4.3)
Thus, for vector v = v i g i = vj e j , one has (∇v) mn = ∂ xm u n , and hence

gradv := (∇v) mn e m ⊗ e n = v i j g i ⊗ g j = ∂ j v ⊗ g j . (4.4)
Accordingly, the curl of a vector in the curvinormal basis writes as

Curl v := ( Curl v) k g k = kij v i j g k . (4.5) Curvilinear coordinates. Let q R ∈ ω R be the curvilinear coordinate associated to g R in the sense that g R = ∂q R x G R , with G R := ∂ q R x
, and where x stands for the position vector of a point. Otherwise said, q R is the curvilinear abcissa of the curve with tangent vector τ R . In general one has

∂ q R u = ∂ xi u ∂x i ∂ q R = g R G R i ∂ xi u = G R ∂ R u. (4.6)
Hence the gradient of scalar u reads

gradu = ∂ i ug i = 1 G i ∂ qi ug i , (4.7)
and of vector v = v i g i as

gradv = ∂ j v ⊗ g j = 1 G j ∂ qj v ⊗ g j . (4.8)
We call the curvilinear expression of the gradient the operator ∇ CURV (•) := h j ∂ qj (•)g j with h j := 1 Gj , the jth metric factor. Remark that the ∂ R derivatives do not commute, contrarily to ∂ q R , because of the factors G R . As an example, consider the spherical base system, where

G N = G r = 1, G A = G φ = 1 r , G B = G θ = 1
r sin φ , and q A = φ (polar angle), q B = θ (azimuthal angle); it holds,

(∂ A ∂ r -∂ r ∂ A ) = 1 r 2 ∂ φ , (∂ B ∂ r -∂ r ∂ B ) = 1 r 2 sin φ ∂ θ , (∂ B ∂ A -∂ A ∂ B ) = 1 r 2 sin φ tan φ ∂ θ . (4.9)
Christoffel symbols in the curvinormal basis. By Theorem 3.2, it is easily deduced by identification with (4.1) that the only nonvanishing components of Γ p ij read (with no sum on repeated indices)

Γ N RR * = -ξ, Γ N RR = -κ R , Γ R * N R = ξ, Γ R R * R = γ R * , Γ R N R = κ R , Γ R R * R * = -γ R . (4.10)
Moreover, it is observed that Γ p ij is not symmetric, i.e., Γ p ij = Γ p ij . Therefore, the torsion is nonvanishing, and reads T p ij := Γ p ij -Γ p ji . In the curvinormal basis, it is easily computed that the only nonvanishing components of

T p ij are T R ij = κ R δ iN δ jR + ξδ iN δ jR * + (γ R * -γ R )δ iR * δ jR .
Note that the Riemann curvature tensor is defined as [START_REF] Dubrovin | Modern geometry -methods and applications, Part 1[END_REF] Riem

q ijk := ∂ k Γ q ij -∂ j Γ q ik + Γ p ij Γ q pk -Γ p ik Γ q pj . (4.11)
Spherical system. As an example, in a spherical coordinates/components system, it holds 1 i, j ∈ {φ, θ}, κ R = 1 r , γ φ = 1 tan φ , γ θ = 0, and hence

Γ r ij =   0 0 0 0 -1 r 0 0 0 -1 r   , Γ φ ij =   0 1 r 0 0 0 0 0 0 -1 r tan φ   , Γ θ ij =   0 0 1 r 0 0 1 r tan φ 0 0 0   . (4.12)
Hence, the torsion reads

T r ij = 0, T φ ij =   0 1 r 0 -1 r 0 0 0 0 0   , T θ ij =   0 0 1 r 0 0 1 r tan φ -1 r -1 r tan φ 0   . (4.13)
Accordingly the covariant derivative reads

v i j ij =   ∂ r v r 1 r ∂ φ v r - v φ r 1 r sin φ ∂ θ v r -v θ r ∂ r v φ 1 r ∂ φ v φ + vr r 1 r sin φ ∂ θ v φ -v θ r tan φ ∂ r v θ 1 r ∂ φ v θ 1 r sin φ ∂ θ v θ + vr r + v φ r tan φ   ij . (4.14)
Hence the curl of a vector, ( Curl v) i := ikj v k j , writes by (4.5) and (4.14) as

( Curl v) i =   1 r tan φ v θ + 1 r ∂ φ v θ -1 r sin φ ∂ θ v φ 1 r sin φ ∂ θ v r -1 r v θ -∂ r v θ 1 r v φ + ∂ r v φ -1 r ∂ φ v r   i . ( 4
T ij k = ∂ k T ij -Γ l ik T lj -Γ l jk T il . (4.16) Let T ij = v i j .
Then by (4.16) one has v i jk := (u i j ) k and hence

v i jk g i = ∂ k v i j g i -(Γ l ik v l j + Γ l jk v i l )g i = ∂ k v i j g i + v i j Γ i lk g l -(Γ l ik v l j + Γ l jk v i l )g i = ∂ k v i j g i -Γ l jk v i l g i
, where (4.3) and a change of dumb indices have been used. Therefore,

∂ k (∂ j v) = v i jk + Γ l jk v i l g i . (4.17)
In particular,

v i jk -v i kj = (∂ k ∂ j -∂ l ∂ k )v • g i -T l jk v i l , (4.18)
1 Here, φ denotes the polar, and θ the azimuthal coordinate, respectively.

that following [START_REF] Dubrovin | Modern geometry -methods and applications, Part 1[END_REF] can be rewritten by means of the Riemann curvature as

v i jk -v i kj = Riem iqkj v q -T l jk v i l . (4.19)
Remark that in spherical coordinates and by (4.13) and (4.9), (4.18) yields (with a slight abuse of notations, since one writes q R instead of R as covariant differentiation indice)

v i rφ -v i φr = v i rθ -v i θr = v i θφ -v i φθ = 0, (4.20)
that is, the second covariant derivatives commute in spherical coordinates/components. In particular, one has ljk v i jk = ljk (v i j ) k = 0 in spherical coordinates/components. Note that the identity Curl ∇u = 0 holds simply by Schwarz lemma if u is a scalar, whereas its vector counterpart Curl ∇v reads

Curl ∇v = Curl (∇v) CART ij e i ⊗ e j = Curl (∇v) CURV ij g i ⊗ g j ,
that is, in the Cartesian system, Curl (∇v) CART ij e i ⊗ e j =ljk (v i,j ) ,k e i ⊗ e l = 0, by Schwarz lemma. However, in a general curvilinear system, one has

Curl (∇v) CURV ij g i ⊗ g j = ∇ CURV (∇v) CURV ij × g i ⊗ g j + (∇v) CURV ij Curl CURV (g i ⊗ g j ) = 0,
and hence

∇ CURV (∇v) CURV ij × (g i ⊗ g j ) = (v i j ) k g i ⊗ lkj g l = -(∇v) CURV ij Curl CURV (g i ⊗ g j ). (4.21)
Remark that in the spherical system, the first term on the RHS of (4.21) vanishes by (4.20).

Summarizing, Eq. ( 4.21) shows that the non-commutation operator in the RHS of (4.19) is related in the curvilinear system to the curl of the basis diads. This fact will appear crucial in the calculations of the following sections.

4.3.

Expression of the incompatibility in the curvinormal basis. Now, the incompatibility operator on a second-rank tensor T is defined as 2 inc T := Curl Curl t T, (4.22) which in a Cartesian system is equivalent to writing componentwise inc T = ikm jln ∂ k ∂ l T mn e i ⊗ e j .

In a general curvilinear system, Eq. (4.22) shows that it suffices to express the curl of a tensor and apply twice the curl operator. In fact, (4.22) rewrites as

inc T = Curl Curl t T CART ij e i ⊗ e j = Curl Curl t T CURV ij g i ⊗ g j , with Curl T CURV ij g i ⊗ g j = ∇ CURV T CURV ij × g i ⊗ g j + T CURV ij Curl CURV (g i ⊗ g j ). (4.23)
Remark that, as compared with the Cartesian system case, the second term in the RHS is nonvanishing and requires to compute the curl of the basis diads. Summarizing, one has

Curl t T CURV ij g i ⊗ g j = S CURV ij g i ⊗ g j , (4.24)
for some components S CURV ij obtained by rearranging (4.23). Hence, the incompatibility in the curvilinear system writes as

inc T = ∇ CURV S CURV ij × g i ⊗ g j + S CURV ij Curl CURV (g i ⊗ g j ) = ∇ CURV S CURV ij × g i ⊗ g j + S CURV ij ∇ CURV g i × g j + g i ⊗ Curl CURV g j = η ij g i ⊗ g j .(4.25)
Obviously, η ij is symmetric as soon as T is, by the symmetry property of its Cartesian counterpart ikm jln ∂ k ∂ l T mn . Moreover, its explicit expression only requires to determine the gradient of scalar S CURV ij in the curvilinear system, which is expressed by means of the tangent vectors as

∇ CURV S CURV ij = DS CURV ij (x)[g l ]g l = g l • ∇S CURV ij g l
2 For some authors, e.g. L. E. Malvern [START_REF] Malvern | Introduction to the mechanics of a continuous medium[END_REF] the incompatibility is defined with a minus sign.

(see also Eq. (4.7)), together with the curvilinear differentials of the basis tensors to be found in Section 3.3 (by means of Eq. (4.1), (4.3) and (4.5)), and expressed by means of κ R , γ R and ξ, which are intrinsic numbers of the boundary as related to the choice of the basis. Note that use has been made of identity

Curl a ⊗ b = ∇a × b + a ⊗ Curl b.
(4.26)

Incompatibility in the spherical system

We recall that we adopt the convention that φ be the polar and θ be the azimuthal angles. Moreover x = rg r , with r the radius. The spherical system consists of the triad {g r , g φ , g θ }, with according to our conventions, N = g r , τ A = g φ and τ B = g θ , the latter two being tangent to the sphere of radius r and normal vector g r .

We consider a symmetric tensor T whose general form reads

T = T rr g r ⊗ g r + T φφ g φ ⊗ g φ + T θθ g θ ⊗ g θ + 2T rφ g r g φ + 2T rθ g r g θ + 2T φθ g φ g θ . ( 5
Curl (g r ⊗ g r ) = - 1 r (g φ ⊗ g θ -g θ ⊗ g φ ), Curl g r ⊗ g φ = 1 r (g r ⊗ g θ -g θ ⊗ g r ), Curl g r ⊗ g θ = - 1 r (g r ⊗ g φ -g φ ⊗ g r ) + 1 r tan φ g r ⊗ g r , Curl g φ ⊗ g r = 1 r g r ⊗ g θ + 1 r tan φ g θ ⊗ g φ , Curl g φ ⊗ g φ = - 1 r tan φ g θ ⊗ g r + 1 r g φ ⊗ g θ , Curl g φ ⊗ g θ = - 1 r (g r ⊗ g r + g φ ⊗ g φ ) + 1 r tan φ g φ ⊗ g r , Curl g θ ⊗ g r = - 1 r g r ⊗ g φ - 1 r tan φ g φ ⊗ g φ , Curl g θ ⊗ g φ = 1 r (g r ⊗ g r + g θ ⊗ g θ ) + 1 r tan φ g φ ⊗ g r , Curl g θ ⊗ g θ = 1 r tan φ g θ ⊗ g r - 1 r g θ ⊗ g φ . (5.2)
5.2. inc T rr g r ⊗ g r . Let us first compute the curl of T = T rr g r ⊗ g r , by using the formulation

Curl (T ij g i ⊗ g j ) = ∇T ij × (g i ⊗ g j ) + T ij Curl (g i ⊗ g j ).
One has

Curl T = ∇T rr × (g r ⊗ g r ) + T rr Curl (g r ⊗ g r ) = - 1 r ∂ φ T rr g r ⊗ g θ + 1 r sin φ ∂ θ T rr g r ⊗ g φ - T rr r (g φ ⊗ g θ -g θ ⊗ g φ ). (5.3) Hence Curl Curl t T = -∇( ∂ φ T rr r ) × (g θ ⊗ g r ) + ∇( ∂ φ T rr r sin φ ) × (g φ ⊗ g r ) + ∇( T rr r ) × (g φ ⊗ g θ -g θ ⊗ g φ ) - ∂ φ T rr r Curl (g θ ⊗ g r ) + ∂ φ T rr r sin φ Curl (g φ ⊗ g r ) + T rr r Curl (g φ ⊗ g θ -g θ ⊗ g φ ) inc T = - 2T rr r 2 g r ⊗ g r -( ∂ r T rr r - ∂ φ T rr r 2 tan φ - ∂ 2 θ T rr r 2 sin 2 φ )g φ ⊗ g φ + ∂ 2 φ T rr r 2 - ∂ r T rr r g θ ⊗ g θ + 2∂ φ T rr r 2 g r g φ + 2∂ θ T rr r 2 sin φ g r g θ + 2 r 2 sin φ ( ∂ θ T rr tan φ -∂ θ ∂ φ T rr )g φ g θ .
(5.4) 5.3. Complete expression of the incompatibility. Collecting (A.1) and all the computations of the Appendix, one arrives at the general formula.

( inc T ) rr = - 2T rr r 2 + ∂ 2 θ T φφ + T φφ r 2 sin φ - ∂ φ T φφ r 2 tan φ + ∂ r T φφ r + 2T φφ r 2 + ∂ 2 φ T θθ r 2 + 2∂ φ T θθ r 2 tan φ + ∂ r T θθ r - 2∂ φ (T rφ sin φ) r 2 sin φ - 2∂ θ T rθ r 2 sin φ - 2 sin φ ∂ θ ∂ φ T φθ r 2 + ∂ θ T φθ r 2 tan φ , (5.5)
( inc T ) φφ = -( ∂ r T rr r - ∂ φ T rr r 2 tan φ - ∂ 2 θ T rr r 2 sin 2 φ ) + (∂ 2 r T θθ + 2∂ r T θθ r ) - 4T rφ r 2 tan φ -2 ∂ θ T rφ r 2 sin φ + ∂ θ ∂ r T rθ r sin φ , (5.6)
( inc T ) θθ = ∂ 2 φ T rr r 2 - ∂ r T rr r + ∂ 2 r T φφ + 2∂ r T φφ r -2 ∂ r ∂ φ T rφ r + ∂ φ T rφ r 2 , (5.7)
( inc T ) rφ = ∂ φ T rr r 2 + ∂ r T φφ r tan φ - ∂ r ∂ φ T θθ r + 2∂ r T θθ r tan φ - ∂ 2 θ T rφ r 2 sin 2 φ + 2T rφ r 2 + ∂ φ (∂ θ T rθ sin φ) r 2 sin 2 φ + ∂ φ ∂ r T φθ r sin φ , (5.8)
( inc T ) rθ = ∂ θ T rr r 2 sin φ - ∂ r ∂ θ T φφ r sin φ + 1 sin φ ∂ φ ∂ θ T rφ r 2 - ∂ θ T rφ r 2 tan φ - ∂ φ (∂ φ T rθ sin φ) r 2 sin 2 φ + T rθ r 2 - T rθ r 2 tan 2 φ + 2∂ r T φθ r tan φ + ∂ r ∂ φ T φθ r , (5.9)
( inc T ) φθ = 1 r 2 sin φ ( ∂ θ T rr tan φ -∂ θ ∂ φ T rr ) + ∂ r ∂ θ T rφ r sin φ + ∂ θ T rφ r 2 sin φ + ∂ r ∂ φ T rθ r - ∂ r T rθ r tan φ - T rθ r 2 tan φ + ∂ φ T rθ r 2 -∂ 2 r T θφ + 2∂ r T φθ r .

Concluding remarks

In this note a method to compute the incompatibility operator in a system of curvilinear components/coordinates is proposed. Moreover an in-extenso expression of the incompatibility is given in the spherical system. It has been shown that the incompatibility of the elastic strain is directly linked to the dislocation density of a solid. Therefore our method will allow us, in a first step, to compute the energy related to dislocations in spherical inclusions. In a second step, to optimize the location of these inclusions in the body with a view to minimizing or maximizing certain cost functionals. This will be the purpose of future work. To conclude, a simple application of our full expression will be given.

Application: determining the dislocation-induced force in linearized elasticity. Recall the general form of second-order free energy Ψ := W e ( ) + Ŵe ( Curl t , div ) + Ŵdefect (κ, Curl κ, div κ) + Wdefect ( 0 ). Let ϕ e (u) = Ω W e (∇u)dx and ϕ defect (F ) = Ω Wdefect ( inc F )dx. The Fréchet derivative of ϕ e at u in the direction v reads Dϕ e (u) [v] = Ω A∇u∇vdx = -div (A∇u) • vdx, that is, the differential ϕ e (u) := Dϕ e (u) = f , by Riesz theorem. Also, Dϕ defect (F )[V ] = Ω inc W defect ( inc F ) • V dx and we set G := ϕ defect (F ) = inc W defect ( inc F ) which is symmetric and divergence free. Assume also that G is independent of 0 . Now we would like to solve (2.8) with Ω the unit sphere, in the simplified case where A = αI 4 , and taking G = g θ ⊗ g θ , that is, we seek F such that inc (α inc F ) = g θ ⊗ g θ .

The solution of inc T = g θ ⊗ g θ is found by (5.7) as T (r)g r ⊗ g r with T (r) = -1 2 r 2 + c. Moreover, by (5.5) F = F (r)g r ⊗ g r is a solution with F (r) = 1 4α r 4 -c 2α r 2 . Uniqueness is obtained by imposing the natural homogeneous Dirichlet conditions F = Curl t F × N = 0 on ∂Ω, i.e., at r = 1, that is, for c = 1/2. Thus 0 = -1 2α (r 2 -1)g r ⊗ g r . Now the displacement u is obtained by solving -div (A∇u) = f + div (A 0 ), where div (A 0 ) = -2r -1 r g r is a radial force due to the presence of dislocations.

Appendix A. Other terms of the incompatibility A.1. inc T φφ g φ ⊗ g φ . Let us compute the curl of T = T φφ g φ ⊗ g φ . One has

Curl T = ∇T φφ × (g φ ⊗ g φ ) + T φφ Curl (g φ ⊗ g φ ) = ∂ r T φφ g φ ⊗ g θ - 1 r sin φ ∂ θ T φφ g φ ⊗ g r + T φφ r (g φ ⊗ g θ - 1 tan φ g θ ⊗ g r ).
Hence

Curl Curl t T = ∇(∂ r T φφ ) × (g θ ⊗ g φ ) -∇( ∂ θ T φφ r sin φ ) × (g r ⊗ g φ ) + ∇ T φφ r × g θ ⊗ g φ -∇ T φφ r tan φ × g r ⊗ g θ + ∂ r T φφ Curl (g θ ⊗ g φ ) - ∂ θ T φφ r sin φ Curl (g r ⊗ g φ ) + T φφ r Curl (g θ ⊗ g φ - 1 tan φ g r ⊗ g θ ) inc T = ∂ 2 θ T φφ + T φφ r 2 sin φ - ∂ φ T φφ r 2 tan φ + ∂ r T φφ r + 2T φφ r 2 g r ⊗ g r + (∂ 2 r T φφ + 2∂ r T φφ r )g θ ⊗ g θ + 2∂ r T φφ r tan φ g r g φ - 2∂ r ∂ θ T φφ r sin φ g r g θ . (A.1) A.2. inc T θθ g θ ⊗ g θ . Let us compute the curl of T = T θθ g θ ⊗ g θ . One has Curl T = ∇T θθ × (g θ ⊗ g θ ) + T θθ Curl (g θ ⊗ g θ ) = -∂ r T θθ g θ ⊗ g φ + ∂ θ T θθ r g θ ⊗ g r + T θθ r ( 1 tan φ g θ ⊗ g r -g θ ⊗ g φ ).
Hence

Curl Curl t T = -∇(∂ r T θθ ) × (g φ ⊗ g θ ) + ∇( ∂ φ T θθ r ) × (g r ⊗ g θ ) + ∇( T θθ r ) × ( 1 tan φ g r ⊗ g θ -g φ ⊗ g θ ) -∂ r T θθ Curl (g φ ⊗ g θ ) + ∂ φ T θθ r Curl (g r ⊗ g θ ) + T θθ r Curl ( 1 tan φ g r ⊗ g θ -g φ ⊗ g θ ) inc T = ∂ 2 φ T θθ r 2 + 2∂ φ T θθ r 2 tan φ + ∂ r T θθ r g r ⊗ g r + (∂ 2 r T θθ + 2∂ r T θθ r )g φ ⊗ g φ -2 ∂ r ∂ φ T θθ r + 2∂ r T θθ r tan φ g r g φ .
A.3. inc 2T rφ g r g φ . Let us compute the curl of T = 2T rφ g r g φ . One has

Curl T = 2∇T rφ × (g r g φ ) + 2T rφ Curl (g r g φ ) = ∂ r T rφ g r ⊗ g θ + ∂ φ T rφ r g φ ⊗ g θ - ∂ θ T rφ r sin φ g r ⊗ g r + ∂ θ T rφ r sin φ g φ ⊗ g φ ) + T rφ r (2g r ⊗ g θ -g θ ⊗ g r + 1 tan φ g θ ⊗ g φ ).
Hence Curl Curl t T = ∇(∂ r T rφ ) × (g θ ⊗ g r ) -∇( ∂ φ T rφ r ) × (g θ ⊗ g φ ) -∇( ∂ θ T rφ r sin φ ) × (g r ⊗ g r ) + ∇( ∂ θ T rφ r sin φ ) × (g φ ⊗ g φ ) + 2∇( T rφ r ) × (g θ ⊗ g r ) -∇( T rφ r ) × (g r ⊗ g θ ) + ∇( T rφ r tan φ ) × (g φ ⊗ g θ ) + ∂ r T rφ Curl (g θ ⊗ g r ) -∂ φ T rφ r Curl (g θ ⊗ g φ ) -∂ θ T rφ r sin φ Curl (g r ⊗ g r )

+ ∂ θ T rφ r sin φ Curl (g φ ⊗ g φ ) + T rφ r Curl (2g θ ⊗ g r -g r ⊗ g θ + 1 tan φ g φ ⊗ g θ ) inc T = -2 ∂ φ T rφ r 2 + T rφ r 2 tan φ g r ⊗ g r - 4T rφ r 2 tan φ g φ ⊗ g φ -2 ∂ r ∂ φ T rφ r + ∂ φ T rφ r 2 g θ ⊗ g θ -2 ∂ 2 θ T rφ r 2 sin 2 φ + 2T rφ r 2 g r g φ + 2 ∂ φ ∂ θ T rφ r 2 sin φ - ∂ θ T rφ r 2 sin φ tan φ g r g θ + 2 ∂ r ∂ θ T rφ r sin φ + ∂ θ T rφ r 2 sin φ g φ g θ .
A.4. inc 2T rθ g r g θ . Let us compute the curl of T = 2T rθ g r g θ . One has Curl T = 2∇T rθ × (g r g θ ) + 2T rθ Curl (g r g θ ) = -∂ r T rθ g r ⊗ g φ + ∂ φ T rθ r (g r ⊗ g r -g θ ⊗ g θ ) + ∂ θ T rθ r sin φ g θ ⊗ g φ + T rθ r g φ ⊗ g r -2g r ⊗ g φ + 1 tan φ (g r ⊗ g r -g φ ⊗ g φ ) .

Hence Curl Curl t T = -∇(∂ r T rθ ) × (g φ ⊗ g r ) + ∇( ∂ φ T rθ r ) × (g r ⊗ g r -g θ ⊗ g θ ) + ∇( ∂ θ T rθ r sin φ ) × (g φ ⊗ g θ ) + ∇( T rθ r ) × g r ⊗ g φ -2g φ ⊗ g r + 1 tan φ (g r ⊗ g r -g φ ⊗ g φ ) -∂ r T rθ Curl (g φ ⊗ g r ) + ∂ φ T rθ r Curl (g r ⊗ g r -g θ ⊗ g θ ) + ∂ θ T rθ r sin φ Curl (g φ ⊗ g θ ) + T rθ r Curl g r ⊗ g φ -2g φ ⊗ g r + 1 tan φ (g r ⊗ g r -g φ ⊗ g φ ) inc T = -2∂ θ T rθ r 2 sin φ g r ⊗ g r -2

∂ θ T rφ r 2 sin φ + ∂ θ ∂ r T rθ r sin φ g φ ⊗ g φ + 2 sin φ ∂ φ ∂ θ T rθ r 2 + ∂ θ T rθ r 2 tan φ g r g φ -2 ∂ 2 φ T rθ r 2 + T rθ r 2 - T rθ r 2 tan 2 φ + ∂ φ T rφ r 2 tan φ g r g θ + 2 ∂ r ∂ φ T rθ r - ∂ r T rθ r tan φ - T rθ r 2 tan φ + ∂ φ T rθ r 2 g φ g θ .
A.5. inc 2T φθ g φ g θ . Let us compute the curl of T = 2T φθ g φ g θ . One has Curl T = 2∇T φθ × (g φ g θ ) + 2T φθ Curl (g φ g θ ) = ∂ r T φθ (g θ ⊗ g θ -g φ ⊗ g φ ) + ∂ φ T φθ r g φ ⊗ g r -∂ θ T θφ r sin φ g θ ⊗ g r + T φθ r g θ ⊗ g θ -g φ ⊗ g φ + 2 tan φ g φ ⊗ g r . (A.2) Hence Curl Curl t T = ∇(∂ r T φθ ) × (g θ ⊗ g θ -g φ ⊗ g φ ) + ∇( ∂ φ T φθ r ) × (g r ⊗ g φ ) -∇( ∂ θ T φθ r sin φ ) × (g r ⊗ g φ ) + ∇( T φθ r ) × g θ ⊗ g θ -g φ ⊗ g φ + 2 tan φ g r ⊗ g φ + ∂ r T φθ Curl (g θ ⊗ g θ -g φ ⊗ g φ ) + ∂ φ T φθ r Curl (g r ⊗ g φ ) -∂ θ T φθ r sin φ Curl (g r ⊗ g θ )

+ T φθ r Curl g θ ⊗ g θ -g φ ⊗ g φ + 2 tan φ g r ⊗ g φ inc T = - 2 sin φ ∂ θ ∂ φ T φθ r 2 + ∂ θ T φθ r 2 tan φ g r ⊗ g r + 2 ∂ φ ∂ r T φθ r sin φ g r g φ + 2 2∂ r T φθ r tan φ + ∂ r ∂ φ T φθ r g r g θ -2 ∂ 2 r T θφ + 2∂ r T φθ r g φ g θ .

3. 3 .

 3 Extended tangent vectors in Ω and their curvilinear differentials. Let d be defined in W by

. 1 ) 5 . 1 .

 151 Curl of the diads. By virtue of (4.14) and (4.15), and recalling (4.26), let us first compute the curl of the base diads.


Chris Budd

Felix Dux

Chattering and related behaviour in impact oscillators

One of the most interesting properties of an impacting system is the possibility of an infinite number of impacts occurring in a finite time (such as a ball bouncing to rest on a table). Such behaviour is usually called chatter. In this paper we make a systematic study of chattering behaviour for a periodically forced, single-degree-offreedom impact oscillator with a restitution law for each impact. We show that chatter can occur for such systems and we compute the sets of initial data which always lead to chatter. We then show how these sets determine the intricate form of the domains of attraction for various types of asymptotic periodic motion. Finally, we deduce the existence of periodic motion which includes repeated chattering behaviour and show how this motion is related to certain types of chaotic behaviour.

Introduction

The single degree of freedom impact oscillator comprising a periodically driven oscillating particle, impacting and rebounding against a rigid obstacle has a remarkably rich dynamical structure. Some of the dynamics, such as periodic and certain forms of chaotic behaviour, can be explained by the usual theory of smooth dynamical systems [START_REF] Thompson | Nonlinear dynamics and chaos[END_REF]) but much of it is quite different and due to the essentially discontinuous nature of the impacting process.

Previous papers studying impact oscillators have, in general, either considered solutions which are periodic and repeat after a fixed number of impacts (often one) or are chaotic with an infinite number of essentially random impacts in an infinite time. By restricting the analysis to solutions which repeat after one impact, it is possible to study the bifurcations that these undergo as the parameters in the system are varied. Some examples of these are given in Shaw & Holmes (1983) and [START_REF] Thompson | Nonlinear dynamics and chaos[END_REF] and have subsequently been developed by [START_REF] Whiston | Singularities in vibro-impact dynamics[END_REF], [START_REF] Nordmark | Non-periodic motion caused by grazing incidence in an impact oscillator[END_REF] and [START_REF] Budd | The effect of frequency and clearance variations in single degree of freedom impact oscillators[END_REF].

In this paper we shall take an opposite approach and examine the dynamical behaviour of those solutions of an impact oscillator which have a very large number or indeed infinite number of impacts in a finite time. One example of this form of behaviour is very familiar: an elastic ball released under gravity from above an obstacle will bounce on the obstacle, acting as a simple impact oscillator. If there is only a small energy loss at each impact then it will bounce a very large number (in an idealized case an infinite number) of times before coming to rest in a finite time. Similar systems arise in many applications in engineering where driven components may collide very frequently against each other or against rigid obstacles, resulting in a high wear rate [START_REF] Goyda | A study of the impact dynamics of loosely supported heat exchanger tubes[END_REF]. Engineers often refer to such behaviour as chatter and throughout this paper we refer to the idealized mathematical description of this as chatter as well. In this paper we shall examine the dynamical behaviour related to chatter in a mathematical idealization of an impact oscillator. However, we expect that our conclusions are also relevant to more general systems.

A vibrating mechanical component undergoing chatter can be very simply modelled as an oscillating particle which is on a periodically forced spring and is repeatedly impacting with the obstacle. In this case, chatter is likely to occur if the forcing is toward the obstacle but the particle is rebounding away from it. In such circumstances, the particle may still come to rest on the obstacle after a large number of impacts, but due to the forcing it will move away from the obstacle at a later time and the chattering behaviour will then form part of a more complex motion. It is this form of behaviour that we study in this paper. We shall examine other forms of behaviour closely related to chatter, such as cases where there are a large number of impacts but the particle does not come to rest (incomplete chatter), examples where chatter is closely related to chaotic behaviour and the effect that chatter has upon the domains of attraction of periodic solutions. Surprisingly the study of systems with a large number of impacts closely spaced in time is analytically rather simpler than those with fewer, more widely spaced impacts. Indeed we shall see that studying the rather special dynamics of those solutions of impact oscillators with such a large number of impacts can give very valuable insights into the qualitative form of the more general solutions.

To set the scene we now describe the precise form of the idealized impact oscillator we study in this paper. The general such oscillator is illustrated in figure 1 and comprises a particle on a linear spring at position x(t) ~ owhich is periodically forced at a frequency w. When x = othe mass impacts with an obstacle and we model this by a simple restitution law so that dxj dt I after collision =r dxj dt I before collision' where r < 1 is the coefficient of restitution.

The complete system is then

d 2 xj~t 2 ~X~ j(t),_ X< 0",} X-7 rx, x-o-.
(1)

Here various constants such as the mass of the particle, the stiffness of the spring and the magnitude of the forcing have, without loss of generality, been scaled to one (see [START_REF] Budd | The effect of frequency and clearance variations in single degree of freedom impact oscillators[END_REF] for the details of this scaling). The forcing functionj(t) will be taken to be periodic with frequency w and to have the property that j(t)-ochanges sign.

An appropriate such function is

j(t) = cos (wt),
and for most of this paper we shall assume that j(t) has this form.

The system (1) departs from a physical model in several respects. In particular, we have not considered damping between impacts and have taken very simple models both for the linear motion on the spring and also for the linear form of the impact. These simplifications make the analysis of the behaviour much easier. Certain forms of behaviour will be rather different in the idealized model: for example, during chatter the particle will impact an infinite number of times in the idealized model of an impact but only a large but finite number of times in a more realistic model. Although this does affect the detailed dynamics, the most significant consequence of chattering-that the particle eventually comes to rest-occurs in both models and it is this which has the dominant effect upon the subsequent motion.

The system (1) has three parameters w, cr and rand its behaviour depends subtly upon the values these take and a review of some of this is given in Shaw & Holmes (1983) and [START_REF] Budd | The effect of frequency and clearance variations in single degree of freedom impact oscillators[END_REF].

If the particle is released from the obstacle at a time t and a velocity x =rv < 0 then it will in general describe a trajectory x(t) which comprises a sequence of smooth motions between impacts. In figure 2 we illustrate three such motions for the case cr = 0 and r = 0.8 with w = 2.6, 2.8 and 2.7 respectively. In the first figure, following an initial transient, the motion becomes and stays periodic. In the second the motion is chaotic and following the initial transient is governed by a strange attractor. In the third the particle repeatedly collides with the obstacle and then comes to rest, sticking to it while its acceleration, j(t)-cr is positive, and being released from it when its acceleration is first negative.

It is the motion illustrated in figure 2 (c) which is of most interest to us here and this figure demonstrates the two phenomena of the particle coming to rest after a chattering sequence and then sticking to the obstacle. In figure 3 we give a close up of this part of the motion. A trajectory of the form illustrated in figure 2 (c) can either be part of a transient or it can in fact repeat and form part of a periodic motion.

There is a closely related phenomenon illustrated in figure 2 (b) where several low velocity impacts occur close together in an incompletely chattering sequence. These sequences frequently form part of a chaotic motion.

To study the phenomena of chatter, incomplete chatter and sticking it is convenient to describe the system (1) by a map P, and to do this we use the impact map introduced in Whiston (1992) and Shaw & Holmes (1983). To define this map we make two observations of system (1). First the trajectory x(t) of the particle is uniquely defined by the velocityrv 0 at which it leaves the obstacle and the phase ¢ 0 of the time t 0 at which it leaves, where ¢ 0 = (t 0 ) mod (2rt/w).

(2)

Second, as has been shown in [START_REF] Whiston | Singularities in vibro-impact dynamics[END_REF], if v 0 is strictly positive then the trajectory necessarily impacts the obstacle again with velocity v 1 ~ 0 at a phase ¢ 1 .

Following the impact the motion continues with initial velocityrv 1 • Thus we may define a map P acting on the phase space [0,2rt/w] x W by

If v 1 > 0, then by simply running time backwards we obtain a unique inverse for P.

As we are interested in examples where the particle comes to rest on the obstacle, we must extend this definition of P to the cases when v 0 = 0. If v 0 = 0 and x > 0 so x(t)

(b) \ 1"'\ M n f"\ M r N x(t) (c) x(t) ~ \ N N \ N v Figure 2
. Time series of (a) periodic motion arising when w = 2.6 and 0' = 0; (b) chaotic motion arising when w = 2.8 and 0' = 0; (c) chattering motion followed by sticking and subsequent release of the particle arising when w = 2.7 and 0' = 0.

Figure 3. A close up of the chattering motion, followed by sticking and release at the time ¢a• that f( ¢ 0 ) -O" > 0 then the particle remains stuck to the obstacle until the first point <Pa > 0 at which where we assume that j(t) is such that this point occurs. The set of values of¢ for which i is positive when x = (}is termed by [START_REF] Whiston | Singularities in vibro-impact dynamics[END_REF] the sticking region of the map P and we shall assume that this is given by

I= [¢p, 2njw] U [0, tPal
If(} = 0 and f(t) = cos (wt) then ¢a= nj2w and tPp = 3nj2w; moreover if(}> 1 then I is the empty set.

When combined with chatter, the sticking region plays a major role in the dynamics of the impact oscillator. This is because it has a large domain of attraction :?i' comprising those trajectories which come to rest following a chattering sequence. Part of the purpose of this paper is to construct the set :?i'. This set occupies a significant proportion of the domain of definition of the map P and consequently is likely to be robust to changes in the model given in (1). It is highly significant that for all cpEJ, P(¢, 0) = P(¢'", 0), and hence P is not invertible over the set I.

By far the most significant aspect of P, realized by Shaw & Holmes (1983) and developed by [START_REF] Whiston | Singularities in vibro-impact dynamics[END_REF], D. Chillingworth (unpublished notes), [START_REF] Nordmark | Non-periodic motion caused by grazing incidence in an impact oscillator[END_REF] and Foale & Bishop (1992) is that it is discontinuous on a one dimensional set S 1 in the phase space and, more importantly, that it introduces considerable stretching onto the phase space in a neighbourhood of S 1 . We shall review this behaviour in the next section. Away from the setS\ the map P is smooth and it behaves very like any other nonlinear map of two variables, with fixed points that undergo various standard forms of bifurcation. However, the discontinuity and stretching result in the many new forms of behaviour remarked upon in the literature.

A very good illustration of this richness of the behaviour of impact oscillators and the role played by chatter is provided by studying the domain of attraction of the fixed points of P for certain parameter values. If, for example, w = 2.6, (} = 0 and r = 0.8 then P has both the fixed point A = (!w¢jn, v) = (0.48214, 0.50544) and the period-6 orbit B, where B = {(!w¢/n, v)} = { (0.285, 0.554), (0.486, 0.907), (0.632, 0.443),} (0.272, 0.563), (0.481, 0.934), (0.635, 0.472).

Both A andB correspond to periodic motions of the impact oscillator of periods 2njw and 12njw respectively. A large number of numerical simulations indicate strongly that these are the only two forms of stable asymptotic motion of the impact oscillator for these ranges of the parameter values. In figure 4 we plot the domains of attraction for these two orbits for a range of initial conditions, [ ¢, v] E [0, 2nj w] x [0, 4]. Also indicated on this figure are the orbits A, B and the sticking region I.

The figure has a beautiful and complex form demonstrating the interesting fact that even if a system has only two periodic asymptotic states, the asymptotic motion of the particle can be very sensitive to its initial condition. The sensitivity is caused by the large and intertwined domains of attraction of A and B resulting from the stretching and discontinuity of the phase space introduced by P. The complexity of the domain of attraction can be understood by examining the structure of the set S 1 , its preiterates and the intersection of these sets with the domain of attraction, :?i' of the sticking region!. We can see part of this structure by observing from figure 4 that there is a neighbourhood f!2J of I all of which is attracted to B as is each of its preimages under P. This neighbourhood is indeed a subset of :F which is composed of f!2J and its pre-images. Moreover, the domains of attraction of A and B close to the boundary, 'f?, of f!2J have a very fine structure with close by points being attracted to either attractor. The images and pre-images of both f!2J and its boundary form large loops in the phase space and these can also be seen in figure 4. A major purpose of this paper is an explanation of this structure in detail. Similar structures are observed for the domains of attraction of different orbits for other parameter values and are due to the same underlying mechanism.

The outline of the rest of this paper is as follows. In §2 we review some of the theory of the stretching behaviour of the map P and of the discontinuity set 8 1 . In §3 we demonstrate the existence of both chattering and incomplete chattering behaviour, and construct an invariant set f!2J for P such that all orbits starting in f!2J lead to chattering motion followed by sticking. By examining the behaviour of P close to the boundary of f!2J we show that it can be simplified to a one dimensional map with a consequent simplification in the analysis of the dynamics. In §4 we compute the pre-images of ::0 to calculate the domain of attraction ff of the sticking region, demonstrating the existence of loops in the boundary of this set. Using this information we deduce the form of the domains of attraction for the attractors of P.

In §5 we show how chatter can form part of a periodic motion and show how this is related to certain types of chaotic behaviour. Finally, in §6 we draw some conclusions from this work.

The stretching associated with P

The map P which was defined in § 1 has many interesting properties and some of these are reviewed in [START_REF] Whiston | Singularities in vibro-impact dynamics[END_REF], [START_REF] Nordmark | Non-periodic motion caused by grazing incidence in an impact oscillator[END_REF] andFoale &Bishop (1992). Of most interest to us here are the two facts that it can have periodic points x such that pm(x) = x for some integer m, and that it is discontinuous on a set Sl, and smooth away from 8 1 .

The periodic points also correspond to (subharmonic) trajectories which repeat after a time 2nnjw and we label them by the indices (m, n). As parameters such as w and a" are varied, the periodic points display similar behaviour to the periodic points of smooth maps (including period-doubling and saddle-node bifurcations) until one of the iterates of x intersects the line 8 1 .

For the remainder of this section we discuss some of the basic properties of 8 1 and the behaviour of P in a neighbourhood of 8 1 . Although it is rather technical we include it because of the insight it gives into the general dynamics of the particle. (This discussion will review the results presented in [START_REF] Nordmark | Non-periodic motion caused by grazing incidence in an impact oscillator[END_REF], [START_REF] Whiston | Singularities in vibro-impact dynamics[END_REF] and Budd & Dux (1994) and will necessarily be very incomplete.)

We suppose that we have three impacts at points A, B, 0 which have similar phases and velocities. We suppose further thatP(A) = (¢ 1 , v 1 ) where v 1 is close to zero and that P(B) = (¢ 2 , 0) where ¢ 1 is close to ¢ 2 , but that the trajectory with initial condition 0 has a maximum with phase ¢m close to ¢ 2 so that the next impact occurs at a rather different phase ¢ 3 , so that P(O) = (¢ 3 , v 3 ).

The map Pis then discontinuous at the point B. In a sense this discontinuity is artificial in that P 2 (A) ~ P 2 (B) ~ P(O) and, indeed, the trajectories starting from A, B and 0 are similar. However, the discontinuity introduces considerable sensitivity to the initial conditions (or stretching in the phase space) and this plays a major role in the overall behaviour of the impact oscillator.

(a) Local behavio'Ur of P close to S

To make this behaviour precise we introduce the following notation first used by [START_REF] Whiston | Singularities in vibro-impact dynamics[END_REF] for the discontinuity set and its iterates and pre-iterates: and so that

8 1 = {(¢, v): P(¢, v) = (¢ 2 , 0)}, sn = {(¢, v): pn(¢, v) = (¢ 2 , 0)} wn = {(¢, v): (¢, v) = pn(¢2, 0)}, WI= P 2 (S 1 ).
The sets wn and sn are dual; indeed by exploiting the time reversal symmetry of the impact oscillator we have

8 1 = {(¢, v): [(2rcjw)-¢, rv] E WI}.
Close to 8 1 , the map P causes considerable stretching of the phase space and the effect of this is illustrated in figure 5 and discussed in detail in [START_REF] Whiston | Singularities in vibro-impact dynamics[END_REF], The subset Rcomprises the initial data which is mapped to low velocity impacts and the effect of these is to stretch Rto a set intersecting the set v = 0 which is parallel to the vector ( = ~), whereN = j{¢ 2 ) -(J' is the local acceleration of the particle at the point P(B). Indeed the stretching factor approaches infinity as we approach B. On the further action of the map P, the set P(R-) is contracted in the direction ( _: N) and is mapped to a set which is locally tangent to Wl, intersecting W in a set which includes the point P 2 {B). The map from Rto P 2 {R-) introduces a pronounced stretching between two data points connected by a line orthogonal to 8 1 , but not between two points connected by a line parallel to 8 1 . In contrast, the subset R+, comprises those points which must miss an impact at the point P(B). This set is not stretched at all by P, and the set P(R+) approaches W transversally, with P(R+) and P 2 {R-) lying on either side of W. (We note that the stretching behaviour applies to much more general systems than the (linear) impact oscillators discussed here.)

(b) The global behaviour of the sets sn and wn

It was first recognized in [START_REF] Whiston | Singularities in vibro-impact dynamics[END_REF] that the global dynamics of the map P is largely governed by the global structure of the sets sn and wn. To gain some insight into this structure we display, in figure 6, the set

8 1 U 8 2 U 8 3 U W U ff 12 U WS 2 v 0 0.5 ~wcp!'lt Figure 6. The set S 1 U S 2 U S 3 U W 1 U W2 U wa.
1.0 when r = 0.8, 0' = 0 and w = 2.6. (It is interesting to compare this figure with the domain of attraction for the same parameter values presented in figure 4.) Some of the structure of these sets, particularly as 0' is varied, can be complex and [START_REF] Whiston | Singularities in vibro-impact dynamics[END_REF] and D. Chillingworth (unpublished notes) devote much space to describing it. However, the general outline of figure 6 can be described fairly simply. To do this we first consider the set Wl and its iterates wn. The nature of the sets sn is essentially equivalent to that of wn and may be deduced in a similar manner by reversing the direction of time in the impact oscillator and replacing r by 1/r.

The sets Wn are given by Wn = pn(¢, 0) for 0 < ¢ < 2njw. If ifJEl then the particle remains stuck to the obstacle and we may consistently define wn = {¢, 0) for all values of n. At the point¢ = ¢a the particle becomes unstuck and the point pn( ¢a, 0) is than an end-point for wn. Equivalently the point p-n(¢p, 0) is an end point for sn. As ¢ increases away from ¢a, the curves wn describe smooth loci until (for example) Wl intersects 8 1 at the point C when (say)¢= ¢ 8 . At this point we apply the results of the last section to calculate the locus of W2. Indeed wn divides into two disconnected curves for ¢ < ¢ 8 and ¢;?: ¢ 8 one of which (given by ¢ < ¢ 8 ) approaches (but does not intersect) Wl at the point A indicated on figure 6 and the other intersects the line {v = 0} when¢= ¢ 8 . Similarly, wa intersects W1 tangentially at A when¢;?: ¢ 8 . By reversing time a very similar behaviour is observed when 8 1 intersects fV1 in that S 2 is divided into two subsets, one of which approaches 8 1 transversally at the indicated point B and the other intersects the line {v = 0}. Moreover, the set S 3 intersects 8 1 tangentially at B.

As ¢ increases to ¢p the sets wn tend toward the point (¢p, 0) and have an approximately parabolic form close to this point. Similarly the sets sn tend toward parabolae passing through the point (¢", 0). We describe this behaviour in detail in §3.

The simple observations account for much of the structure of figure 6 and this basic structure is, in general, preserved for other values of w, (]'and also when damping between impacts is included in the model. In the following sections we see how, by including chatter into the dynamics, we can predict much of the global behaviour of P (and the forms of the domains of attraction) from these observations of the forms of wn and sn.

Chatter and the invariant chattering region

We now turn our attention to the main theme of this paper: namely the existence of chattering orbits and the role played by chatter on the dynamics of the impact oscillator. Loosely speaking, a chattering motion is one that involves an infinite number of low velocity impacts with the obstacle before the particle comes to rest against it with some phase ¢ 00 in the sticking region. Similar phenomena have been observed in the (chaotic) motion of billiard balls [START_REF] Katok | Invariant manifolds, entropy and billiards: smooth maps with singularities[END_REF] where they are referred to as 'dead end orbits'.

To be precise, we define an orbit to include chatter if it contains a sequence ( ¢n, vn) such that as n-+oo,

where ¢ 00 EI. We can also have various forms of incomplete chatter in which there are a large number of low velocity impacts in a short time interval, but at some critical value of N, ¢ N > ¢" and, as the acceleration of the particle is now away from the obstacle, the phase of the next impact ¢N+l is then rather different from ¢N• Any motion which includes chatter must include the point(¢", 0). Hence if chatter repeats then the orbit is necessarily periodic. However, the motion can still be very complex, especially if the image of the point (¢", 0) is close to a saddle-point of the mapP.

In this section we show that chattering orbits define a basin of attraction !!J for the sticking region I, formed by the pre-images of the orbits which include chatter.

Moreover, !!J is bounded by a set ((5 which is an accumulation set for the sets sn as n-+oo.

For simplicity, our discussion will be mainly concerned with the forcing function j(t) =cos (wt). However, many of the conclusions apply as well to more general differentiable, periodic forcing functionsj(t) such thatj(t)-(]'changes sign at the two points ¢" and ¢p•

(a) The dynamics of a chattering motion

To demonstrate the existence of a chattering motion when j(t) =cos (wt) we presume that the particle has an initially low velocity for some value of ¢ 0 EI, where ¢ is not too close to the end points of I. In this case we have a positive acceleration toward the obstacle and we may presume that the time to the next impact is small. If we take dxjdt = -rv 0 ~ 1 at phase ¢ 0 then d 2 xjdt 2 = cos(w¢ 0 )-(J', d 3 xjdt 3 = -wsin(w¢ 0 )+rv 0 , etc. Thus, if A = ¢-rp 0 , we may express the local behaviour of the particle in terms of a Taylor series so that x = O"-rv 0 A+!(cos(wrp 0 )-0")A 2 -!(wsin(wrp 0 )-rv 0 )A 3 +0(A 4 ),

(3)

x =-rv 0 +(cos (w¢ 0 ) -O") A -!(w sin (w¢ 0 ) -rv 0 ) A 2 + O(A 3 ).

(4)

The time to the next impact is given by the first value of A for which x(A) = 0". If v 0 is small then A is also small and to leading order in v 0 the value of A is given by A = 2rv 0 j[ cos ( w¢ 0 ) -0"].

(

) 5
This approximation is reasonable provided that (w sin (w¢ 0 ) + rv 0 ) A ~ (cos (w¢ 0 ) -O"), so that v 0 ~ (cos ( w¢ 0 ) -0") 2 / ( w sin ( w¢ 0 ) + rv 0 ).

(6)

We assume this to be the case for the present and return to it in the next section.

It follows from (5) that the next impact occurs at the phase ¢ 1 , where ¢ 1 = rp 0 +2rv 0 /[cos(wrp 0 )-0"].

(7)

Combining (4) and (5) gives the velocity at the next impact as

(8)
Thus, provided that condition (6) is met, (7) and (8) give a simplification of the impact map P for small velocities v 0 and for rpEl. A similar simplification can be made for a wide class of forcing functions j(t).

For any value of rp 00 El this map has a fixed point given by (rp,v)=(¢ 00 ,0) corresponding to a chattering trajectory which has come to rest at this point. If rfln is close to rp 00 and vn close to 0 then we can linearize the map about this point to give rfln+l ~ rfln + 2rvnf[ cos (wrp 00 ) -0"], vn+l ~ rvn-

Thus the local behaviour of the iterates of Pis given by

Vn ~ rnv 0 , rfln ~ rp 00 -2rn+lv 0 /[(1-r)(cos(wrp 00 )-0")];

(9) (10) an infinite number of collisions of ever decreasing velocity occurring in a finite time.

The formulae ( 9) and ( 10) give the approximate dynamics of P close to the fixed point and it is evident from this that the points (r/Jn, vn) lie on the straight line passing through (¢ 00 , 0) given by vn = !(1-r) (cos (w¢ 00 )-0") (¢ 00 -rfln)fr.

(11)

The gradient of these lines is -!(1-r) (cos (w¢ 00 ) -O")jr which takes its maximum value when rp 00 = 0 and its minimum of zero when rp 00 = rfla or rp 00 = ¢p• When studying the trajectory of a chattering motion we should, in general, expect to see a set of points which approaches the line v = 0 along such a straight line. In figure 7 we illustrate such a trajectory for which the first few iterations do not have a clear structure but eventually do become asymptotic to a straight line. question by calculating a neighbourhood[!) of I which is mapped into itself by P and has the property that any trajectory starting in it inevitably experiences chatter before coming to rest. Surprisingly, we can calculate this region by examining the behaviour of the map P close to the end points(¢'", 0) and (¢p, 0) of I. We do this by considering the action of P on a series of parabolae ~parametrized by A passing through the point (¢'", 0). The structure of the map P close to ¢rx is then completely characterized by the effect it has on each such parabola and we express this in the following lemma. In other words, a parabolic set of initial conditions is mapped to another such set. (This result remains true for any sufficiently smooth forcing function.)

Proof. We again use the Taylor series decomposition

x = tr-rv 0 A +!(cos (w¢ 0 ) -tr) A 2 -!(w sin (w¢ 0 ) -rv 0 ) A 3 + O(A 4 ).
We now presume that (¢ 0 ,v 0 )E~ and that ¢ 0 = ¢rx-£, where cos (w¢J-tr = 0 and c is small. Then To solve this to find A such that x = tr, we put A= cT, where Tis of order 1, giving 0 = -rwA. +!wsin (w¢'") T-~wsin (w¢J T 2 + O(c 2 ).

If we then ignore terms of O(c 2 ) we find that T satisfies the quadratic equation ! sin (w¢J T 2 -!sin (w¢J T + rA = 0.

The corresponding equation for x is x = -rE 2 wA + Ewsin (w¢J A -!wsin (w¢J A 2 + O(c 2 A).

Substituting A = cT gives the velocity v 1 at the next impact as (12)

We take the negative root in this expression to ensure that T gives the phase of the first impact after ¢ 0 . Combining our expressions gives the formula for j in the Lemma.

D

The significance of this result is that close to <Pa we can simplify the map P to a one dimensional map from the set ~ to the set ~ characterized by the single function f.

After some manipulation, we can show that j has the following properties:

j(O) = 0; df/dA(O) = r < 1; d 2 f/dA 2 (A)
> 0 for all A< sin(w¢J/3r; d 2 J/dA 2 (0) = 20 r 2 /[3 sin (w¢J]; j(A) -+00, T-+ 1 as A-+ sin (w¢J/3r.

The general form ofj(A) is illustrated in figure 8, together with the line g(A) = A, and the intersection of these demonstrates the existence of a fixed point of f.

Combining these simple results gives the following important lemma.

Lemma 2. (i) The Junction j(A) has a unique non-zero fixed point A 00 where

! sin (w<fia)/r > A 00 > 0 such that j(A 00 ) = A 00 and a further fixed point A= 0.

(ii) If sin (w¢a)f3r > A > A 00 then j(A) > A.

(iii) If A< A 00 thenj(A) <A andr(A)-+0 as n-+oo.

By using the local information on j close to the origin we have that for small A Lemma 2 is a most important result because if we are sufficiently close to the point (¢'", 0), so that the approximations we have made are negligible, we can deduce from it the existence of a curve which is invariant under the action of P.

Corollary. If¢ is close to ¢rx then the curve is invariant under the action of P.

Any chattering orbit for which ¢ 00 = ¢rx must ultimately lie on the parabola T 00 and we illustrate this in figure 9. Indeed, a point on T 00 will move monotonically to the right, converging to the point (¢'", 0).

In fact, the arc of the parabola T 00 with end-point (¢'",0) is a limit for an invariant curve of the mapP as ¢--'J>¢rx• Further away from ¢rx thetermsof0(£ 2 ), which we have neglected in the expressions for x and v, become important and while P continues to have an invariant curve, this set departs from being a simple parabola. (We note that this is not a rigorous proof of the existence of such a curve which would need more careful arguments than those given here such as a fixed point construction in an appropriate function space. However, the arguments here do demonstrate how such a curve can be constructed.) To compute this invariant curve we take a small arc of But T must also satisfy the quadratic equation Hence to leading order, 8 1 has a component given by the parabola T;,_ where 1 A1 =~sin (w¢Jjr >!sin (w¢J/r > A00 • Proceeding as before, sn has a component given by T;,_ where n An = r 1 (An-1) and by the previous calculation, An> A 00 ,An < An_ 1 and An-l>-A 00 as n-?-00.

D

The above lemma gives a precise description of the sets sn close to the point (¢a, 0) which is in accord with the observations recorded in figure 6.

We conclude that under the assumptions about the fixed points of P, the curve CtJ and the line v = 0 bound a region ~ with the following properties.

Theorem 1. (i)

The region ~ is an invariant region for P and any orbit starting within it experiences chatter before sticking.

(ii) The boundary of~ is an invariant curve CIJ.

(iii) CtJ is an accumulation set for the curves sn as n _,.. oo.

We can say a little more about the shape of ~ by considering the points of intersection between CtJ and the sets wn considered in the last section. These sets have a structure which, in a neighbourhood of the point (¢p, 0) is very similar to that of sn and is given by the following result. (iii) As n _,..oo, wn _,.. {v = 0}.

Proof. The proof of this result is similar to that for 8 1 and follows by taking a Taylor series expansion close to (¢p, 0). For simplicity we shall omit it.

D

The map g(A) is very similar to f(A); however it has the significant difference that dgjdA (0) = r, d 2 g/dA 2 (A)< 0, VA> 0.

Thus it has no fixed points other than A = 0 and hence An_,.. 0 so that wn_,..{v = 0} as n-?-00. The resulting form of !(!, and of the sets sn and wn in a neighbourhood of I is illustrated in figure 11.

To complete our description of the invariant set!'» we note from figure 11 that the sets S 1 and W1 intersect at a point A 1 which is close to!(/. Similarly, W1 intersects the sets sn at a sequence of such points An accumulating at the point A 00 where W1 intersects 1(/.

We define Bn to be the pre-image of An so that Bn = p-l (An) and

Bn ESn+l n {v = 0}.
It follows thatBn-+B 00 , whereB 00 = (¢>y,O) is the point where!(/ intersects the line v = 0, so that if Bn = (¢>n, 0) then ¢>n -+¢>1' <¢>pas n-+oo. We can also deduce the structure of!(! close to the point B 00 , by studying the pre-image of the set sn close to its intersection with Wl. It follows from the discussion in §2 (b) that this is a portion of the set sn+l passing through Bn, and that close to Bn this set approximates a straight line of gradient 1/[r(o--cos (<f>n))]. Taking the limit as n-+oo we deduce that close to Boo the curve!(/ approximates a straight line of gradient 1/[r(o--cos (¢>y))].

(c) Incomplete chatter We conclude this section by looking at those trajectories which have a large number of low velocity impacts but do not ultimately stick to the obstacle. Such a trajectory must include a point which lies close to the curve !(/ but not inside the region !'». The iterates of this point will then initially lie close to !(/ and will have low velocity impacts, but will eventually leave a neighbourhood of !(/ after a finite number of iterations. The closer the initial point is to !(/, the more low impacts will be included in the motion. We can estimate the curve dividing those low velocity impacts that lead to further such impacts at a nearby time from those that lead to high velocity impacts at a rather later time, by the set S 1 . Indeed, those points lying to the 'left' of S 1 are mapped to a neighbourhood of!(/, whereas those lying to the 'right' are mapped to that part of fVl which is close toP( if> a, 0). We illustrate this in figure 12 by plotting the iterates of P corresponding to a trajectory which includes a sequence of impacts with incomplete chatter.

Domains of attraction

For various values of w and o-, the map P has distinct periodic attractors (for example w = 2.6, o-= 0) and, as we saw in the introduction, the domains of attraction of each such attractor can have a complex form. Much of this complexity can, however, be understood by examining the images and pre-images of the set 5» constructed in the last section. This follows partly from the following simple result.

Lemma 5. All points in 5» and its pre-images have the same asymptotic behaviour as the single point (if>a, 0). Thus, if a trajectory starting from (if> a, 0) is attracted to a periodic orbit then so will all trajectories with initial conditions in 5» and its pre-iterates.

We now calculate the form taken by the pre-iterates of 5». The point A1 ES 1 n fVl constructed in §3 has the pre-iterateB 1 ES 2 n (v = 0). By applying the results of §2 (b) we may deduce that B 1 has, in turn, the pre-iterate 0 1 lying on the intersection of sa with S 1 and which is also a limit point of the set S 2 . If we now compute the pre-image of that portion of S 2 lying between B 1 and A 2 we see that it forms a part of sa which intersects S 1 at the point 0 1 and approaches S 1 again at the point 0 2 which is the preimage of B 2 such that 0 2 ES 1 n S 4 . Indeed, sa forms a loop around the end-point of S 1 which is the pre-image of the point (r/>p, 0). Following the discussion in §2 it is evident that sa intersects S 1 tangentially at 0 1 and approaches it transversally at 0 2 .

A similar argument applied to the portion of sn lying between Bn_ 1 and An, shows that the pre-iterate of this set forms a loop around the end of S 1 intersecting it tangentially at Cn_ 1 and approaching it at Cn. We illustrate this behaviour in figure 13. components of the set $' which, by the arguments above, will form a series of loops around the end-points of the sets sn. This short discussion explains the loops visible in figure 4 which are, in turn, the boundaries of the components of the set$'.

To see more precisely why this behaviour occurs, we consider the images of the sets Vn, defined to be the subsets of WI lying between the points An and An+l but not including the point An+1• A simple calculation shows that the set pm(Vn) is a subset of wm+l that lies between sn-m and sn-m+l if m < n, and between S 1 and the line {v = 0} if m = n, intersecting {v = 0} at the point Dn = pn(Anl• Furthermore, the set pn+l(Vn) is a subset of wn+ 2 that intersects WI tangentially at the point En = pn+l(An), loops round the end point, E 00 , of WI, given by P((if>a,O)), and approaches WI transversely at the point En+1• We illustrate this behaviour in figures 15 and 16. These figures are closely related to the strange attractors of various forms of chaotic behaviour and we shall return to this in the next section.

D 1 =P(A 1 )
Figure 16. The approximately triangular set L1 1 and two of its iterates.

The behaviour of these curves is very similar to that of the curves sn close to the endpoint of S 1 . However, there is a significant difference in that as n--+00,

Dn-+(¢rx,O), En-+P(¢a,O).

The approximately triangular region Lin, bounded by the sets Vn, sn, sn+l and the point (¢a, 0), is mapped by pn onto the region bounded by pn(Vn), S 1 and the line v = 0 and by pn+I onto the region Gn, on the end of Wl, that is bounded by the loop made by wn+ 2 and that part of JV1 lying between the points En and En+1• Furthermore, as n-+oo the sets Gn accumulate onto the single point E 00 •

The precise form of the domain of attraction of P close to C(/, and for the preiterates of these points, is then determined by the form of this domain in a neighbourhood of the single pointE oo-If we consider the special case in which P has precisely two attracting states A and B, then in this case two possibilities may arise.

(i) Eoo is interior to the domain of attraction of A (or B).

(ii) E 00 is on the boundary of the domain of attraction of A (or B) (or, perhaps, is very close to the boundary).

In case (i) the sets G n (which decrease in size as n --+00) will ultimately all lie inside the domain of attraction of A (or B). Consequently, close to C(/, all trajectories will be attracted to A (or B). In contrast, in case (ii) the boundary of the domain of attraction of A (or B) will intersect G n for all n and hence, some points of Vn will be attracted to A and some will be attracted to B. As a consequence each region Vn will contain 'stripes' which are attracted to A or B. This behaviour is illustrated in figure 17 and accounts for the structure observed in figure 14 and also in figure 4.

Periodic chatter and chaos

Finally we study the occurrence of trajectories in which chatter is repeated. Such trajectories are necessarily periodic. Closely related to these are chaotic trajectories with recurrent examples of interrupted chatter. A trajectory with chatter necessarily evolves in subsequent time as the iterates Fn = pn( efta, 0) of the single point (efta, 0). These iterates are also the 'end-points' of the sets wn. Such a trajectory will be periodic if FN lies inside £» for some value of N < oo (or alternatively if FN-lEtff etc.). As the sets£» (and Iff) have non-zero measure it is likely that this will occur for a range of values of the parameters r, wand rr. For example, if rr and rare fixed then we would expect to see an interval of values [w 1 , w 2 ] at which periodic chatter is observed, such that at w 1 and w 2 then FN E C(f is on the boundary of £», and for w 1 < w < w 2 then FN lies inside £». Indeed if r = 0.8 and rr = 0 then such an interval occurs for 2.692 < w < 2.701 and N = 5.

The resulting iterates Fn for the particular case w = 2.7 are illustrated in figure 18.

It follows, by construction, that if w 1 < w < w 2 , then the end-point E 00 ( = F 1 ) of WI must lie interior to the domain of attraction of the periodically chattering orbit. First, this implies that any such an orbit must be stable. Second the domain of attraction must then be of the form (i) described in §4. Indeed the domain of attraction for the periodically chattering orbit is rather large and is illustrated in figure 19.

We now suppose that w is slightly greater than w 2 so that FN lies just outside!!).

Instead of experiencing chatter the further iterates of FN will form an incomplete chattering sequence. Indeed there will be a further point pN+M such that pN+M ~ (¢"',0) and pN+M+l ~F1 . In this case pzN+M will be close to FN but, owing to the effects of stretching will be further away from !!) than FN is. There will be a sequence of such points F 3 N+ 2 M, F 4 N+aM etc. each giving an incomplete chattering sequence until the points are too far from !!) for incomplete chatter to be observed. This behaviour is illustrated in figure 20. The subsequent behaviour of the points Fk is then going to be very different from the chattering sequence and will probably be attracted to a very different asymptotic state. We conclude that the periodic chatter observed for wE [ w 1 , w 2 ] does not deform continuously into another form of periodic motion when w leaves this interval but is likely to be replaced by a very different form of asymptotic state. Indeed, if we continue to consider the above example and let w increase from 2.701, we find that the iterates of P have a chaotic behaviour. The strange attractor of this motion when w = 2.8 is illustrated in figure 21 (the corresponding trajectory is given in figure 2 (b)). Some of the structure of the strange attractor can be understood in terms of trajectories which include incomplete chatter. In particular it is shown in [START_REF] Whiston | Singularities in vibro-impact dynamics[END_REF] that the strange attractor of a chaotic orbit is, in general, close to the set Woo= U~~1 wn. Now, the analysis in the previous section, summarized in figure 15, showed that as a result of the intersection of WI with Cfl, the set Woo includes a series of loops connected to WI and accumulating on the point P(¢a, 0). Orbits mapped onto these loops necessarily include sequences of interrupted chatter. Inspection of figure 21 shows that the strange attractor certainly inherits a great deal of the structure of woo with the loops close to the end point of JVl clearly visible. The related sequences of incomplete chatter are visible in figure 2 (c). Very similar behaviour can be observed in the strange attractors of the chaotic trajectories of orbits with different parameter values and will occur in general whenever one of the iterates of the point (¢a, 0) lies close to the curve Cf/. For the case o-= 0 this occurs when w is close to an odd integer and partly explains the fact that chaos is frequently observed for these values (see the discussions in Shaw & Holmes 1983; [START_REF] Budd | The effect of frequency and clearance variations in single degree of freedom impact oscillators[END_REF].

Conclusions

We have shown that chatter can arise quite naturally in the idealized impact oscillator and that it helps explain both the forms of the domain of attraction of periodic attractors and the strange attractors of chaotic motions. The discussions in this paper have been of a rather idealized model, but it is significant that the set :?2 calculated is large and will therefore be stable to perturbations of the model which make it more realistic. (Although we would have to modify the precise definition of chatter in the case of different impact laws.) As the dynamical significance of chatter is largely related to the interaction of :?2 and its iterates with the set 8 1 , it is probable that some of the predictions of this paper are also true for more realistic models. We leave this as a subject for future investigation.

Figure 1 .

 1 Figure 1. The general form of the impact oscillator.

Figure 4 .

 4 Figure4. The domains of attraction for the (1, 1) (light shade) and the (6, 6) (dark shade) periodic orbits when w = 2.6 and 0" = 0. The (1, 1) orbit is indicated by a cross and the (6, 6) orbit by six circles.

  Figure 5. The stretching of the phase space produced by the map P close to the discontinuity set 8 1 . Nordmark (1991), Foale & Bishop (1992) and Budd & Dux(1994). In this figure we consider a rectangle R = R+ U R-of initial data which intersects the set 8 1 and is bounded by line segments through the points A, BES and 0 which were considered earlier. Thus the subsetsR+ andR-lie on either side of 8 1 andR-includes part of 8 1 .

Figure 7 .

 7 Figure 7. The iterates of P for an orbit which includes chatter.

Lemma 1 .

 1 Let~= {(¢,v): v = wA.(¢'"-¢) 2 ,¢ ~ ¢J; then if (¢-¢rxl is sufficiently small,where tt = j(A.) and j(A.) = ( -rA. +sin (w¢'") T-! sin (w¢J T 2 )/(1-T) 2 , where = ~-h/[9-24rA.jsin (w¢'")].

  cos (w¢0 ) -tr = cos (w(¢'"-c))-tr = we sin (w¢J + O(c 2 ), wsin(w¢ 0 ) = wsin(w¢rx)+O(c). Now, by definition, v 0 = WA£ 2 • Thus x = tr-r£ 2 wA.A + !cw sin (w¢J A 2 -!(w sin (w¢J-rwA.c 2 ) A 3 + O(c 2 A 2 ) + O(cA 3 ) + O(A 4 ).

  v 1 = E 2 w[ -rA+sin(w¢a)T-!sin(w¢JT 2 ]. Now, the phase of the next impact is ¢ 1 = ¢ 0 +A, so that ¢1 = <Pa-E+Tc = <fia-E(1-T), or <Pa-¢1 = c(1-T). Combining this with (12) gives v1 = Wf1(¢a-¢1)2, (¢1,v1)E~, where 11 = [ -rA+sin(w¢a)T-!sin(w¢JT 2 ]/[1-TJ 2 . Now, T satisfies the quadratic equation T 2 -3T+ [6ritjsin (w¢J] = 0. So T = ~±h/[9-24ri\.jsin(w¢J].

Figure 8 .

 8 Figure 8. The function J(it).

Figure 9 .

 9 Figure 9. The parabola Too and the locus of a chattering orbit for which ¢oo = ¢a•

T 2 -

 2 3T + [6rAjsin (w¢J] = 0. Combining these two equations gives T = ~ and A= ~sin (w¢Jjr.

Lemma 4 .

 4 Let~= {(v, ¢): v = wA(¢-¢p) 2 , ¢;?! ¢p}. (i) If¢ is sufficiently close to ¢p then P(~) = ~' where and T satisfies the quadratic equation T 2 +3T-[6rAjsin(w¢a)J = 0.(ii) In a neighbourhood of (¢p,O) v) = ~n' where An= g(An_1 ).

Figure 11 .

 11 Figure 11. The form taken by the sets sn and wn in the neighbourhood of the set~.

Figure 12 .

 12 Figure12. The iterates of a trajectory which includes a sequence of incomplete chatter.

Figure 14 .Figure 15 .

 1415 Figure 14. The 'striped' nature of the domain of attraction close to the region !?fl.

vFigure 18 .

 18 Figure 17. The intersection of the sets wn with the domain of attraction of point A shown shaded and the resulting pre-iterates of this region.

Figure 19 .

 19 Figure19. The domain of attraction of the (2, 2) (shaded), (8, 6) (white) and periodically chattering (black) orbits that occur when w = 2.7 and tJ = 0.

Figure 21 .

 21 Figure21. The strange attractor for the chaotic motion that arises when w = 2.8 and 17 = 0.

the parabola T 00 with (rpa, 0) as an end-point and iterate this set backwards under the action of p-1 . This defines a complex set of curves which are not all connected and some of these are illustrated in figure 10.

From this figure we identify an arc, also with end-point (rpa, 0), which further intersects the line v = 0 at the point (¢1', 0) where

The image of (rpy, 0) under Pis the first point where the backward iterates of the arc of T 00 intersects the set JVI. We denote this curve by ((5, and the region enclosed between itself and the line v = 0 by !?fl. We make the assumption at this point that the map P has no fixed points lying either on ((5 or in !?fl and the closely related assumption that the set 8 1 does not intersect ((5. It follows from these assumptions that Pis smooth and invertible within !?fl. Now, as both ((5 and the line I are invariant under the action of P, it then follows that !?fl will be an invariant set (under the forward action of P) and, by the assumption on the fixed point of P, any trajectory with an initial point in !?fl must inevitably lead to chatter.

We can show further that, under these assumptions, the curve ((5 is an accumulation set for the curves sn. We establish this in the following lemma. Lemma 3. (i) If¢ is close to rfla, then 8 1 has a component which coincides with the parabola 'qsin(wg\.)fr• (ii) As n -+00' sn has a component which coincides with the set of parabolae TA where i\.n is monotonically decreasing, f(i\.n) = i\.n_ 1 , i\. 1 =~sin (wrpaJjr and i\.n-+ i\. 00 as nn -+00.

Proof. The discontinuity set 8 1 is the set of points (¢ 0 , v 0 ) for which P(r/Jo,Vo) = (¢1,0).

If we assume that ¢ 0 ~ rfla and we set v 1 = 0 in ( 12 The resulting curves accumulate onto the pre-image of the portion of C(J which lies between Boo and Aw-This curve loops around the endpoint of S 1 and intersects it twice at the point 0 00 which is the limit of the points 0 n as n -+00. The loop encloses a non-empty region Iff containing the point p-1 (r/Jp,O) such that P(cff) is that subset of 5» bounded by C(J, fVl and the line v = 0. Any point in Iff will form an initial point of a trajectory with chatter. The pre-images of Iff are similarly non-empty regions containing the points p-n(r/Jp, 0). Any point in such a region is ultimately mapped into 5» and hence is attracted onto the same trajectory. We denote the union of these regions, together with 5», by ff, which is then the domain of attraction of the sticking region I. The general form of ff is clearly visible in figure 10 and also in the domains of attraction of the two periodic attractors for P given when w = 2.6, o-= 0 which is presented in figure 4. In the latter example the set ff is visible as a union of sets, each centred on the end points of a set sn and all attracted to the (6, 6) orbit.

We can extend our understanding of the domains of attraction of the fixed points of P by considering the asymptotic behaviour of the trajectories with incomplete chatter which have an initial point close to C(J. As we have shown, C(J is an accumulation of the sets sn and these sets lead to considerable stretching of the phase space. Consequently, we would expect that the domain of attraction will have a complex form close to C(J. Our discussions in §2 showed that a set of data 'parallel' to a set sn is not stretched by the action of P whereas a set of data which intersects sn transversely is considerably stretched by P. This simple observation indicates that two points joined by a chord parallel to sn are likely to have the same asymptotic behaviour and be mapped to the same attractor, whereas two points joined by a chord which intersects sn may be mapped to different attractors. The resulting form of the domain of attraction of the two at tractors given when w = 2.6, o-= 0 in a neighbourhood of C(J is illustrated schematically in figure 14 and takes the form of a series of 'stripes' roughly parallel to C(J. The pre-iterates of this neighbourhood of C(J are in turn a series of stripes surrounding the various