A Le Van

J Royer

Boundary formulation for three-dimensional anisotropie crack problems

Keywords: anisotropie elasticity, crack, integral equation

Boundary integra-differentiai equations for three-dimensional anisotropie cracked bodies are derived. Both the cases of the infinite body (with an embedded crack) and a finite body with an embedded or surface crack are considered. Detailed mathematical conditions for the results to be va/id are specified throughout.

Introduction

In the field of numerical techniques, the most widely known is undoubtedly the domain finite element method (FEM). In recent years, the structural components engineers have to cape with become more and more complex, and the number of degrees of freedom entailed by a sufficiently refined mesh is continually increasing to obtain more and more accurate results. lt was precisely for this reason that the advent of the boundary integral equation method (BIEM) is of great significance since it allows the structure to be meshed over the boundary only, at least when the body forces are absent. In any case, the BIEM allows the solving of problems with smaller matrices because the unknowns are exclusively boundary quantities.

The pioneering works in the BIEM were those of Rizzo, 1 Cruse and Rizzo, 2 and Cruse. 3 • 4 The regularized expression of the BlE was first given by Rizzo and Shippy. 5 Since then the boundary method bas been extensively developed showing its ability to deal with various types of mechanical problems. In the field of fracture mechanics, the BlE, or more precisely the boundary integra-differentiai equation (BIDE), bas been formulated for elastostatic, thermoelastic as weil as elastodynamic problems. 6 -10 Regularized BIDEs have also been given elsewhere. 11 -14 A complete review of different works in the BlE field can be found in the paper of Tanaka et al., 15 though related to regularization techniques.

As a matter of fact, whereas most of the materials are more or less anisotropie, the BIEM was mainly developed for isotropie materials. This can be accounted for by two reasons. First, equations for the isotropie case are simpler to solve from the theoretical and numerical standpoints alike. Second, the elastie material properties are much more difficult to be determined experimentally in anisotropy. Nevertheless, increasing use is being made of anisotropie structural components, and this requires more efficient studies for this class of materials. Using the decomposition of the Dirac function into plane waves, Vogel and Rizzo 16 derived the integral representation of the fundamental solution for a general anisotropie elastic three-dimensional (3-D) continuum. Later on, an efficient numerical implementation for anisotropie problems was proposed by Wilson and Cruse. 17 In fracture mechanics, Sladek and Sladek 18 and Balas et al. 19 (pp. 50-52) have discussed the boundary formulation for anisotropie cracked bodies, and the corresponding BIDE bas been proposed, conjecturing that sorne results in isotropy can be extended to anisotropy.

In the first part of this paper, the integral representation of the fundamental solution for an anisotropie elastie medium derived by Vogel and Rizzo 16 is briefly reviewed and its basic properties are investigated. In particular, relations describing the limit behavior of the fundamental solution near a closed or open surface are presented. Ali the results obtained are generalizations of the well-known ones in the isotropie case. As an application, in the second part of the paper, the BlE and BIDE are derived for the prob\em of ani"&otropic cracked bOO\'è'i>. Th<ë oou'i\-dâr'.f formulation includes both the cases of the infinite body (with an embedded crack) and a finite body with an embedded or surface crack. Throughout the paper, emphasis is made on the mathematical conditions fo r the results to be obtained.

The fondamental solution and its basic properties

In this section, the fundamental solution for a homogeneous anisotropie linear elastic medium is recalled and its basic properties are given. The components of the fourthorder elastic tensor C written in a fixed base (e 1 e 2 e 3 ) of the 3-D space g are material constants verifying the usual symmetries:

(1)

Let us introduce the following definition where the summation convention is implied over repeated subscripts, which ali have the range (1, 2, 3).

Definition

Given a unit vector em of the base (e 1 e 2 e 3 ), the fundamental solution related to em and denoted by U(em, x, y) is the solution of the partial differentiai equation in the infinite space e:

div[C: grad U(em, x ,y )] + 8(y -x)em = 0 (2)
where My -x) is the Dirac function.

The corresponding stress tensor is defined as:

= C: grad U(em, x, y)
where the tensor product a ® b is defined by (a ® b ); 1 = a,.bi. The differentiation in equation ( 9) is performed with respect to variable y, and tu(em. x,y)(y, n y) denotes the stress vector at point y with respect to normal n Y and corresponding to the displacement field U( e m, x, y). Now the fundamental displacement tensor is defined by:

E(x,y) = U(em,x,y) ®em, i.e., E,./x, y)= lf;(ei, x, y) (4)
This, in turn, gives rise to the third-order tensor of the fundamental stress:

(5) Eventually, the Kupradze tensor is defined by 20 (p. 99):

The relationships between the above-defined tensors are straightforward:

T;k(x , y, n y ) = 'I;/ek, x, y )n/y) = Dijk(x , y )n/y) [START_REF] Weaver | Three-dimensional crack analysis. !nt[END_REF] Conversely:

Because tensors I and T are functions of D , any equation in the sequel cao be expressed in terms of D alone. In practice, however, the simultaneous use of notations T and D proves to be more convenient.

The equilibrium equation in terms of the fundamental solution reads

'tf k' \1 i' div I ( e k, x, y) + 8 (y -x) • e k = 0, aD k k,-'-1 (x, y)+ô(y -x)•ô;k =O a yi i.e.,
which gives rise to the so-called rigid body identity:

where n is any bounded region and S its boundary. [START_REF] Sladek | Dynamic stress intensity factors studied by boundary integra-differentiai equations[END_REF] We now come to the expression of the fundamental solution and its derivatives. For this purpose, let us introduce the following notations. Given two points x and y of space g , let (a 1 , a 2 , a) be an orthogonal base with the third vector a equal to e, = (y -x)jlly -xli, a 1 and a 2 being arbitrary. This based allows us to define the spherical coordinates (x , lj!) E [0, 7T[ X [0, 27T[ such that the line x= 0 coïncides with e" whereas the origin for 1/1 ts arbitrary.

Let us define the tensor Q as Q,.k(() = C ijkt (j( 1 , Ç i is the j-component of vector t in the global ftxed base (ep e 2 , e 3 ) . The inverse of Q is denoted by P:

P(() = Q -1 ({), i.e ., 1
P;k( () = 2 det Q 8 kpq eirsQp,Qqs
where ekpq is the permutation symbol. Both P and Q are symmetric tensors. Since a unit vector ~ can be entirely determined by its spherical coordinates (x, 1/1) in the base (op Oz, a), tensors P( n and Q( n are also written as P(a , x , 1/1) and Q(a, x, 1/1 ).

The following theorem was proved by Vogel and Rizzo 16 by decomposing the Dirac function into plane waves.

Theorem

The fondamental solution for an anisotropie elastic medium is given by: [START_REF] Bui | Régularisation des équations intégrales de l'élastostatique et de l'élastodynamique[END_REF] in which r = li y -xli and P(a, 1/J) stands for P(a, x= 7T /2, "').

1 i21T Vm, U(em, x, y)= -2 - P(a, 1/J) dr/rem 87T r o 1 f21T =E(x,y)=- 8 2 P(a,!/f)d!/1 7T r o
The integral is taken along the unit circle in the plane normal to a =er and passing through x.

The derivatives of tensor E(x, y) can be computed by [START_REF] Bui | Regular BlE for three-dimensional cracks in elastodynamics[END_REF] where r,; =dr j dy; = er •e; and Q(a, 1/J) stands for Q(a, x= 7T /2, "').

The following properties result from the above theorem.

• The symmetry of E:

E(x,y) =ET(x,y) (13)
• By applying the same reasoning made for E(x, y ) to E(y, x), we obtain the variable interchange properties: • The singularity of the fundamental solution:

E(x,y) =E(y,x) aE aE -( x, y)=--(x, y)
E(x,y) = o( ~) D(x,y) = o( r\ ) T(x ,y,ny) = o(, 1 2 )
asr=lly -xii~O

(16) (17) (18)
• The asymptotic behavior at infinity:

E(x,y) = o( ~) D(x,y) = o( ,\) T (x, y, n) = 0 (, 1 2 )
as r = Il y -x Il ~ oo [START_REF] Halas | Stress Analysis by Boundary Equation Methods[END_REF] Of course, in isotropie or transversely isotropie cases, relations [START_REF] Bui | Régularisation des équations intégrales de l'élastostatique et de l'élastodynamique[END_REF] and [START_REF] Bui | Regular BlE for three-dimensional cracks in elastodynamics[END_REF] simplify appreciably, yielding wellknown expressions. 21 The fundamental solution for a general anisotropie medium is known only by the integral representation [START_REF] Bui | Régularisation des équations intégrales de l'élastostatique et de l'élastodynamique[END_REF], but the expression in closed form is not available in general since the actual integration is not possible for any time of anisotropy. Nevertheless, the above properties prove to be sufficient to establish the so-called limit theorems in what follows, without the knowledge of any closed form expressions for E and T whatever.

Limit theorems

The purpose of this section is to present sorne results about the limit behavior of the fondamental solution when the load point x approaches a point y 0 belonging to a given surface S. These results hold in the case of general anisotropy, thus constituting the generalization of the well-known ones in the isotropie case. They ali are based on the !emma below, which alone exploits the integral representation [START_REF] Bui | Regular BlE for three-dimensional cracks in elastodynamics[END_REF] of the fondamental solution.

First let us introduce sorne notations. Given a point y 0 and a unit vector Dyo• let n be the plane passing through y 0 and perpendicular to Dyo• The normal Dyo defines two si des of plane n, which will be recognized by the signs + and -, the + side being the half-space containing the point y 0 +nyo• Let 8 be an arbitrary positive number, B(y 0 , 8) the bail centered at y 0 , and with radius 8. The plane n divides the boundary a B(y 0 , 8) of the bail into two parts denoted by S 1 (y 0 , 8) and Sz(y 0 , 8 ), situated in the sides + and -, respectively.

L em ma V y 0 , Vnyo• "18, Vk , l,p , q , [START_REF] Kupradze | Three-Dimensiona/ Prob/ems of the Mathematical Theory of Elasticity and Thermoelasticity[END_REF] where the normal n Y is outward to B(y 0 , e ).

Proof:

Re garding the integral over S 2 (y 0 , 8 ), relation ( 12) yields:

1 f7r /27T { f27T = -87T2 ll=O cp=o'.t r,q 1/f=Oppk(a,l/f)d!/J ( / 27T aQ,m + Snq-r.nr.q) PP 1 (a, 1/J )--(a, 1/J) 1/J=O
oan XPmk(a, 1/J) dl/f} sin ()d0dcp [START_REF] Willis | The elastic interaction energy of dislocation loops in anisotropie media[END_REF] in which r =il y-y 0 ll, r.q = (yq-Yoq)jr, e, =(y-y 0 )jr,

the spherical coordinates ( 0, cp) E [0, 7T[x[O, 27T
[ are so defined that the Iine (} = 0 coïncides with n whereas the origin for cp is arbitrary. As for the integral Y~ver S (y 8)

it is identical to that of ( 21), only the bounds fo~ o' ar~ different, since now (} varies from 7Tj2 to 7T. Using the variable change ()' = 7T-0, cp' = 7Tcp and introducing n~w notations: y' denoting the symmetrical point of y wtth respect to y 0 , r' = lly'y 0 il, r' = ar'jay' =(y' -

) / 1 d 1 ,q q q
Yoq r an a = -a, we have:

r aEpk ], -a-(y 0 ,y)n 1 (y)dYS SICYo• e) Yq 1 Jo -7T = --- J (-r') 8 7T 2 0' = 7T 12 cp' = 7T ,/ x {c-r,' q) f 2 7T ppk( -a''"') dl/f .P=O + ( 8 1 1 )f21T p ( 1 ) aQ/m nq-'.n',q pl -a'"' a( 1) 1/f=O -an X (-a', 1/1 )Pmk( -a', 1/1) dl/f} sin(}' d(}' dq/ (22)
Eventually, by comparing relations [START_REF] Willis | The elastic interaction energy of dislocation loops in anisotropie media[END_REF] and ( 22) and by noting that: '

we obtain the proof of the lemma. Applying relations (5) and ( 6) to relation [START_REF] Kupradze | Three-Dimensiona/ Prob/ems of the Mathematical Theory of Elasticity and Thermoelasticity[END_REF], we directly deduce the following propositions.

Proposition (23) Proposition Vy 0 , Vnyo• 'V8, Vi,j, k, l, (24)
1t should be noted that the above propositions involve the plane II passing through ;>' o an~ normal to vector "ro• Moreover, the value of 8 JS arb1trary, not necessarily small. !n t~e subsequent applications where y 0 belongs to an arbitranly shaped surface S, the neighborhood of y in S is not plane, this is why we shall have to take the li':nit 8~0.

Relations ( 23) and ( 24) allow us to prove the so-called limit theorems below. Let us first agree about notations for the orientation of a surface. Of course, any surface S

(closed or open) considered here is assumed to be orientable. This implies that, for any point y 0 ES '\. as, we can locally define two sides of S which we label side + and si de -, ail the normals to S being directed from si de -to side +.

Theorem

Let S be a surface (closed or open) and u a vector field defined on S. If: • Let us denote r =li y -xli, r 0 = lly 0 -xli and p =li yy 0 ll. From (18) and hypothesis (ii), we have:

T(x, y, n y) = O(ljr 2 ) = 0[1/( p 2 + rJ)], u(y)-u(y 0 ) = 0( p 13 ),
and dYS = 0( p). Hence:

lim J T(x, y, ny)[u(y)-u(y 0 )] dYS = 0(8f3)
x--+y 0 ± S(y 0 , e)

• It remains to investigate the last integral in equation (26). For brevity, let us denote:

where the normal ny for y E S 1 (y 0 ,8)USiy 0 ,8) is outward to the bali B(y 0 , 8) while we recall that n Y for y E S(y 0 , 8) is directed from side -to side +.

The rigid body identity (10) yields:

(

) 27
where 1 is the unit tensor. Similarly, we have:

(28)

On the other band, from hypothesis (i), S 1 (y 0 , 8) and S 2 (y 0 , 8) tend to two symmetrical hemispheres as e --+ 0 e---> 0

Solving equations ( 27) to (29) gives

1 limA±=+-. e -->0 - 2
The theorem is proved.

Theo rem

Let S be a surface (closed or open) and t a vector field defined on S. If:

(i) SE C 1 • a, Û <a~ 1 (ii) tE C 0 • f3(S), 0 < f3 ~ 1 th en (30)
where n y o is the normal vector at y 0 to S and the product of the third-order tensor D and vector t is the second-arder tensor defined by

(D • t)ij = Dijmtm.

Proof'

The proof is similar to that of equation ( 25), bearing in mind that, in view of relations [START_REF] Weaver | Three-dimensional crack analysis. !nt[END_REF] and ( 16):

The following theorem requires a somewhat stronger condition for the function u.

Theorem

Assuming that where the symbol .9f represents the differentiai operator defined as: =.9fijm( a y, x' y, n y)um(y) = CiJk!Dmnk(x, y) .gn,( a y, n y)um(Y) (32) 9Jn 1 (ay, ny) is the tangential differentiai operator defined as (33) (there is no possible confusion of the normal ny with the subscript nE {1,2,3}), the symbol ay recalls that the differentiation is carried out with respect to variable y.

(i) SEC 1 •a,O<a~l (ii) u E C 1 •f3 (S), 0 < f3 ~ 1, i.

Proof"

Invoking arguments similar to those in the proof of (25), we can write:

The main work is the study of Let us also denote where the normal n Y on different surfaces is defined as in the proof of (25).

• First let S 1 (y 0 , e) be involved by writing X [nnum,l(y 0 ) -n 1 um,n(y 0 )] dyS • Djo (34) where n 1 = n 1 (y ). For the sake of brevity, the right-hand side of the previous equation will be written in the form of a difference (35) In virtue of the rigid body identity [START_REF] Le Van | Integral equations for three-dimensional problems[END_REF], we have -Cijmnum,n<Yo)njo Then, substituting equations (36) and (37) into (35) gives:

x Dij/x, y)( -njonq + njnq 0 ) dYS

• 0 m,n<Yo)
where But, from equation (18) and hypothesis (i):

Th us,

• The same treatment can be applied to S 2 (y 0 , e ), with special attention to the normal orientation. The corresponding results for the ± cases are interchanged with respect to those obtained with S 1 (y 0 , e), and we are led to

-A ±+ C =E + O(ea) (39)
where • Furthermore, according to hypothesis (i), S 1 (y 0 , e) and S 2 (y 0 , e) tend to two symmetrical hemispheres as e ~ 0, so that equation ( 24 where as outward with respect to the body considered, in accordance with the usual convention. In the sequel, ali the equations will be written using s;;, so that normal n;

is taken as the reference one. This section gives the displacement at any interior point of the body in terms of boundary quantities. Two cases are considered: the infinite body (with an embedded crack) and a finite body with an embedded or surface crack.

1 nfinite body

Consider the infinite body 0 containing the crack Scr = S~ Us:;. The following hypotheses referred to as the regularity conditions are assumed. The radiation condition for the body forces is assumed as: f(y)=O(r - 2 -8 ) when r~oo, where ô is a positive constant smaller than 1.

Regularity conditions

Note that condition (ii) is identically verified if the body forces f is confined to a finite region. Given any point x, let B(x, R) denote the bali centered at x with radius R large enough for B(x, R) to include the crack. Let SR be the exterior surface of B(x, R) and O.R, be the bounded region within s;;u s;; and SR, with boundary sc--;.u s;;u Sw The second condition (ii) ensures that limR_,cofnRE(x,y)f(y)dYV is bounded.

The following theorem holds.

Theorem

Provided the regularity conditions are fulfilled, we can write the integral representation of the displacement as:

u(x)= j_[E(x , y)It(y)
Sc.

+T 7 (x,y,ny)Àu(y)] d YS

+ * j 0 E(x,y)f(y)dyV (42)
in which It(y) is the sum of the stress vectors on the crack faces: It(y) = t(y +,n;) + t(y-,n ; ) (if the crack is loaded symmetrically, i.e., t(y+,n ; ) = -t(y-,n;), then It(y) = 0) Àu is the displacement jump through the crack: Âu(y) = u(y+)u(y-), and the asterisk denotes an improper integral.

The proof of the theorem does not formally differ from that in the isotropie case.

2 Finite body

Consider a finite body 0 with outer boundary SB, containing a crack Scr• The following theorem can be proved in a similar way as in the infinite body case.

Theorem

The integral representation of the displacement reads:

V xE n " (SB u scr), u(x) = J:_[ E ( x,y)It(y) sn +T 7 (x,y,ny)Âu(y)] dYS + J [ E(x,

y)t(y , n y)

So -T 7

(x ,y, n y)u( y )] dyS

+ * lr/(x, y )f(y)dYV (43) Note that the regularity conditions are useless in the finite body case.

Integral representation of the stress

In practical purposes, it is often necessary to compute the complete stress tensor at any point inside the body. This section is thus devoted to the integral representation of the internai stress.

Infinite body Theo rem

Assuming that:

(i) the regularity conditions are fulfilled.

(ii) S;;E C 1 •a, 0 <a~ 1 (iii) Àu E C 1 (5;;)
we have the integral representation of the stress:

Vx~Scr • u( x ) = J: {D(y,x)It(y ) Sc~ + * j D( y, x)f( y ) drV (44) n Proof"
The integral representation of the displacement (42) holds because of hypothesis (i). Given a point x ~ Scr> the derivative auk(x)j ax 1 will be investigated to obtain the stress u(x). Ali integrais in equation ( 42) can be differentiated behind the integral sign. Indeed, since x is interior to fl the surface integral is regular. Moreover, relation [START_REF] Sladek | Three-dimensional crack analysis for an anisotropie body[END_REF] implies that the kernel of the volume integral behaves as 1/ r as r ~O. Thus:

(45)
The first and third integrais containing aEmd ax 1 do not require any further transformations. As regards the second integral, exploiting hypotheses (ii) and (iii) it can be transformed by means of the so-called regularization the-orem (see Appendix): (46) In writing equation (46), use bas been made of the boundary condition along the crack front: ~u(y) = 0 for y E aSer• The theorem is proved.

Note the order of variables x and y in the D-terms of equation (44): D(y, x) instead of D(x, y). On the other band, the volume integral in equation ( 44) is improper convergent on account of the regularity condition fE O(r-2-o).

2 Finite body Theo rem

Assuming that:

(i) (ii) s-ec 1 •()' s EC 1 •"'0<aa'<l cr ' B ' ' - ~u E C 1 (S;), u E C 1 (S 8 )
we have the integral representation of the stress for a finite cracked body:

'V xE 0" (Sa U Scr), u(x) = j {D(y,x)It(y)
s~ + Bf( a y, x, y, ny) ~u(y )} dYS + j {D(y,x)t(y,n y)

Sa

-.9f(ay,x,y,nY)u(y)}dYS

+ * f D(y, x)f(y) dy V n (47)
In the case of a surface crack, Sa must be replaced by Sa " L throughout, where L = as;n Sais referred to as the surface line.

Pro of'

Only the case of the surface crack, which needs special attention, is investigated; the proof for the embedded crack can be deduced in an obvions way. The reasoning is essentially the same as in the infinite body case, the difference being in that here there are two surfaces: the closed surface Sa upon which the displacement u and the stress vector tare defined and the open surface s; upon which the displacement jump ~u and the stress sum It are defined.

The surface tine L being positively oriented with respect to the normal n; of sc -;, let us introduce the closed contour L -U L + made up to two arcs L-and L + such that both of them coïncide with L whereas the orientation of L-(respectively L + ) is the opposite to (respectively the same as) that of L (Figure 1).

Assume for definiteness that the contour L-U L + is oriented negatively with respect to the outward normal to the exterior boundary Sa, as shown in Figure 1. In fact, it can be easily verified that the final result does not actually depend on this assumption. Hypotheses (i) and (ii) make it possible to apply the regularization theorem (see Appendix) successively to surfaces S; and S 8 " L the boundary of which is L -u L +: =surface integral over sc-; 48) and (49), it turns out that, although line integrais appear in the case of a surface crack, their very sum is zero. The theorem is proved.

Notes

(a) In the case of a surface crack, it has been assumed

that u E C 1 (Sa "L), not u E C 1 (Sa). Indeed, the surface line L corresponds to an incision in the boundary Sa, giving rise to a displacement discontinuity on S 3 along L and making the hypothesis u E C 1 (S

Boundary integro-differential equation

It is well-known in isotropy that the BlE obtained in usual way by taking the limit 0 '-Scr 3 x---+ y 0 E Scr in equation ( 42) is nonunique and insufficient for the determination of the unknowns on the crack. Indeed, the loading on the crack faces is involved in the BlE only through the sum '!.t so that, if the stress vector is continuous, i.e., !.t = 0, there is no information about the loading at ali. The purpose of this section is to overcome this deficiency by obtaining an adequate BIDE, possibly combined with a BlE.

Infinite body Theo rem

Under the following assumptions: where by x ---+ y 0 ± are meant the upper and lower limits with respect to normal n; 0 , respectively. According to the hypotheses, the integral representation (44) of the stress holds. Then multiplying it by n; 0 and taking the limit as x---+ y6' gives (50). The theorem is proved.

If the loads t+ and r are specified, the solution of the BIDE (50) gives the displacement jump du , and then u+ and uusing the integral representation of the displacement (42). Conversely, if du is prescribed, the BIDE (50)

does not allow to obtain r+ and r separately, unless an additional information is supplied, e.g., t+= -t-. This is better accounted for by the particular case of the plane crack considered below.

Plane crack in the isotropie medium

From equation (50) immediately results the following corollary.

Corollary

Consider a plane crack imbedded in the infinite medium and lying in the plane y 3 =O. Assuming that:

(i) the medium is isotropie (ii) the regularity conditions are fulfilled (iii) !luE C 1 (S,;) we have the BIDE for a plane crack:

'fly o E Scr 'Va E {1, 2}, 1 2[ta(yr), e3) -ta(yt,-e3 )] 1-2v 1 = 8 (
(51a) 1 z-ltiy 0 ,e3) -tiyti, -eJ)] 1-2v 1 = - ( )pvj 2Ita(y)r,adyS 87T 1-V Sc~ r IL 1 + 4 ( )pvj-zra~u 3 ady S 7T 1-v s~ r ' ' + W/y 0 ) (51b)
where r= lly-Yoll, r,; =(y;-Yo)fr,

W(yo) = * fo_D(y,yo)f(y)dYV.e3

There is no coupling between (ta, ,lua)a E 11 , 21 and (t 3 , ~u 3 ) if and only if the loading on the crack is symmetrical, i.e., It =O. The case of the plane crack allows us to go deeper into the question of the aforesaid indetermination of t+ and t-. For simplicity, let us consider a circular crack free of body forces (W = 0). The crack is symmetrical with respect to its plane y 3 =O. The prescription of ,lu 3 can be realized by inserting a rigid wedge, axisymmetrical and dissymetrical with respect to the plane y 3 = 0, between the crack faces. It is assumed that there is no friction between the wedge and the crack faces, then Ita = 0, a E {1, 2} and equation (51b) cao be written in the abbreviated form (lj2)(tjt;) = function of ,lu 3 •

The knowledge of ,lu 3 implies that the difference tj-t:; is known, but not tj and 13 separately. This is confirmed when the wedge is turned upside down: the same quantity ~u 3 is prescribed, the difference tj-t3 remains the same, as can be easily verified, but tj and t3 are modified individually. This simple example shows that knowing Au in equation ( 50) is not enough to determine the loads t+ and t -.

2 Finite body Theo rem

Consider a finite body n containing an embedded or surface crack Scr• Assuming that:

(i) s-EC 1 • " s ECI,a' O<a a'<1 cr " B ' ' - (ii) ,luE C 1 •'\S;), u E C 1 (SB), 0 < {3 ~ 1
we have the following system of BlE and BIDE:

+ * j {E(y 0 ,y)t(y,ny) Vyo E scr t(tCy 0 ,n; 0 ) -t(yti,n; 0 )] =pu j {D(y , y 0 )It(y)

s.;

-~(ay , y 0 ,y,nY)u(y)}dY S • n; 0

+ * jnD(y,yo)f(y)dYV.n ;o (53)
where in the case of a surface crack, SB must be replaced by SB '-. L throughout.

Proof' Equation ( 53) is obtained in a similar way as equation (50). To prove equation ( 52), let us transform the integral representation of the displacement (43) by means of the rigid body identity (10):

+ j {E(x,y)t(y,n y )

Ss

-TT(x, y ,ny)[u(y) -u(x)]}dYS

+ * /,E(x,y) f (y)dyV=O !1 (54)
Now, we proceed to the limit in equation (54) as n'-. (SB u s cr) 3 x~ Yo E SB. From hypothesis (ii) and relation [START_REF] Sladek | Three-dimensional crack analysis for an anisotropie body[END_REF], the integral over SB is weakly singular white that over s.; is regular. Thus, the limit is performed by replacing x by y 0 in equation (54). The theorem is proved.

The solution of the system (52) and ( 53) give (u, t) on SB and (Au, It) on scr• The coupling between these equations accounts for the interaction between the outer boundary and the crack. In the case of a surface crack, the system provides no equations along the surface tine, since y 0 $ L. However, the lacking equations are compensated for by expressing the displacement compatibility at the surface line. With notations introduce in equation (49), we cao write at every point on L : Also, it should be noted that the principal value integral in the BIDE (50) or (53) cao be numerically transformed in the manner indicated by Balas et al. (Ref. 19, p. 166) into a regular one which in turn can be easily computed.

Conclusions

The specifie feature of the formulation has been the general anisotropy of the medium. The results have consisted of:

(1) the limit theorems (25), (30), and (31), which describe the limit behavior of the fundamental solution when the Joad point x approaches a point y 0 belonging to a given surface S, closed or open. These theorems are generalization of the well-known ones in the isotropie case. It is noteworthy that they have been obtained without even knowing the closed form expression of the fundamental solution. A minimum amount of basic properties, relations [START_REF] Bui | Régularisation des équations intégrales de l'élastostatique et de l'élastodynamique[END_REF] to [START_REF] Kupradze | Three-Dimensiona/ Prob/ems of the Mathematical Theory of Elasticity and Thermoelasticity[END_REF], has been enough to complete the proofs. (2) the BIDEs for crack problems, equation (50) for the infinite body with an embedded crack, and the system of coupled equations, (52) and ( 53), for a finite body with an embedded or surface crack, which clearly shows the interaction between the outer boundary and the crack.

From the numerical points of view, the success of the boundary formulation in anisotropie problems depends on an efficient computation of the integral representations [START_REF] Bui | Régularisation des équations intégrales de l'élastostatique et de l'élastodynamique[END_REF] and ( 12) of the fundamental solution. Wilson and Cruse 17 have shawn that the anisotropie solution cao be numerically evaluated with essentially arbitrary accuracy. Further improvements in the numerical scheme should allow to more precisely calculate the fundamental solution as weil as to significantly save the CPU time.

Nomenclature

c fourth-order elasticity tensor em a unit vector

•

  Using relation (5) and the foregoing properties, we obtain similar relations for I and D: I(ek,x,y) = -I(ek,y,x) and D(x,y) = -D(y,x) ai ai -(ek,x ,y) = --(ek ,x,y) and ax[

1 =

 1 (i) Sis a Lyapunov surface: SEC 1 • ", O<a~l, which means that 3C > 0, \/y, y' ES, liny'-nyll ~ Clly' -y ll" (ii) u satisfies the Holder condition on S: u e C 0 • f3(S), 0 < {3 ~ 1, i.e., 3C > 0, \lx,y ES, llu(y)-u(x)ll ~ Clly -xii 13 th en lim Jrcx,y , nY)u(y)d s x--->y<f S Y ± 2u(yo) +pu fs T(y 0 , y, n y)u(y) dYS (25) where by x~ y 0 ± are meant the limits as x approaches y 0 , x belonging to the side + and the side -: r~spectively .. The symbol pu denotes a Cauchy pnnc1pal value mtegral. Relation (25) also holds if tensor T is replaced by its transpose Tr. Prooj-By denoting S(y 0 , 8) =Sn B(y 0 , 8 ), the left-hand side of (25) can be recast as • The first integral is continuous at y 0 when y-Yo• It tends to the Cauchy principal value integral appearing in the right-hand side of (25) as 8 -O.

  e., ali the derivatives of u belong to the class C 0 • 13 (S) we have the property of continuity across the boundary:

  ± can be transformed successively as follows. Using equation (5), Now, since S(y 0 , e) U S 1 (y 0 , e) is a closed surface, we have Then, owing to the symmetry of E, relation (14), and relation (5): Using equation (10) and the symmetry CmniJ = CiJmn relation (1), we get and a-= lim cmnpqf x--->ycr S(y 0 , e)US 1 (y 0 , e) X D;j/njonq-njnq 0 ) dYS •um,n(y 0 )

  (i) The radiation condition for unknown displacement and stress is assumed as: u(y)=o(l), O'(y) = o(ljr), i.e., t(y, ny) = o(ljr ), when r = /ly-xli~ oo, x being a fixed point. (ii)

Figure 1 .

 1 Figure 1 . Definition of the closed contour L-u L + in the case of a surface crack.

1 (

 1 (i) the regularity conditions are fulfilled (ii) S,;E C 1 • ", 0 < a :5: 1 (iii) du E C 1 • i3(S,; ), 0 < (3 :5: from hypotheses (ii) and (iii), it follows that 2-t = -[C: (graddu)]•n ; EC 0 •/3 '(sc~), 0<(3':5:1) the BIDE writes Vyo E scr -}[ t(y 0 , n; 0 ) -t(yt, n; 0 )] =pu J {D(y,y 0 )2.t(y) sc-; +~(ay,y 0 ,y,n y ) du(y)} dyS•n;o + * J D(y, y 0 )f(y)dYV• n ; 0 (50) n Proof' Equations (30) and (31) are valid in view of hypotheses (ii) and (iii), yielding and

Ss-

  T T(y 0 , y, ny)[u(y)-u(y 0 )]}dYS + * jE(y 0 ,y)f(y)dyV=O fl (52)

  8 ) impossible. (b) The theorem cannat be proved by considering Sa U Scr as a single closed surface. The crack geometry clearly indicates that, contrary to s,;, scr = s:;u s,; cannat be assumed to belong to the class C 1 • ".

Appendix

Regularization theorem

Let S be a surface (open or closed) with edge aS. If:

where fgnt is the tangent differentiai operator defined by equation (33). If S is a closed surface, the line integral along as is zero.

Proof"

Relations ( 7) and ( 17) give Now, it is easy to verify the following formula of compound derivatives: for any differentiable functions and using the equilibrium equation ( 9) with x =1= y ES, the theorem is proved.


domain: math spi spi.meca.vibr

text:
Wai Chin

Edward Ott

Helena E Nusse

Celso Grebogi

Grazing bifurcations in impact oscillators

Impact oscillators demonstrate interesting dynamical features. In particular, new types of bifurcations take place as such systems evolve from a nonimpacting to an impacting state (or vice versa), as a system parameter varies smoothly. These bifurcations are called grazing bifurcations. In this paper we analyze the different types of grazing bifurcations that can occur in a simple sinusoidally forced oscillator system in the presence of friction and a hard wall with which the impacts take place. The general picture we obtain exemplifies universal features that are predicted to occur in a wide variety of impact oscillator systems.

I. INTRODUCfiON

We say that a system is an impact oscillator if it has an oscillating object that impacts frequently with some other object [START_REF] Moon | Chaotic vibrations[END_REF][START_REF] Thompson | [END_REF]. Impact oscillators occur in many technological situations. For example, mechanical devices are often engineered with loose fitting joints to accommodate thermal expansion, and the dynamical behavior of such systems often leads to impacts in the joint. In addition, many machines inevitably suffer from effects of vibroimpacts. A common feature shared by models of these systems is the smoothness of the systems between the impacts. Shaw and Holmes studied a piecewise linear, sinusoidally forced impact oscillator for various values of the forcing frequency [3][4][5][6]. Whiston originally showed the importance of grazing impacts (i.e., zero velocity impacts) of the global dynamics [7,8). Recently, Nordmark expanded (to first order) solutions in the neighborhood of a grazing orbit for a simple physical system (described below) and obtained a two dimensional map representing the dynamics of lm orbit in the neighborhood of the grazing state [9]. Nordmark also studied the dynamics of this map, obtaining several important results [9,10]; we will give details later. Nusse, Ott, and Yorke [11) obtained results for the dynamics of the one dimensional limit of a two dimensional map equivalent to the map derived by N ordmark. Budd, Dux, and Lamba considered sinusoidally forced impact oscillators, studying such features as chattering, intermittency, the effect of frequency and clearance variations, and the scaling of Lyapunov exponents at nonsmooth bifurations [12].

In this paper, we use the simple physical system shown in Fig. 1 as a prototype impact oscillator. This is the system considered in [9,10). A mass m is attached to a linear spring with spring constant k that is fixed to the wall on the right hand side. There is a sinusoidal external force F 0 sinwt acting on the mass. The friction force is proportional to the velocity of the mass with coefficient Jl• Here s represents the position of the mass m, and t is the time derivative of 5, which is the velocity of the mass m. A hard wall stands at the position Sc• When the amplitude of oscillation is sufficiently small, there are no impacts between the mass m and the wall at Sc• and the dynamics of the system is the same as that of a forced damped harmonic oscillator without the wall at Sc• As the amplitude of oscillation is increased, the mass m begins to have impacts with the wall, first with very low velocity. The bordering state between the impacting and nonimpacting is called a grazing impact, i.e., when the mass contacts with the wall at Sc with zero velocity. Interesting new bifurcations are observed at grazing, and they are called grazing bifurcations [9,10) considered in [11,13] and discussed in Sec. VI). The purpose of this paper is to present an analysis of grazing bifurcations for the system in Fig. 1. It is anticipated that these results are universal in that they apply to many systems in which impacts occur.

The two dimensional map derived for the system in Fig. 1 by Nordmark in [9] is equivalent to the following map, which will henceforth be referred to as the Nordmark map: Here xn and Yn are transformed coordinates in the position-velocity space (s,t> evaluat.ed at times tn, where wt n = 2n 1T, and w is the frequency of the external forcing (see Fig. 1). The quantity 71 is the restitution coefficient of the impacts. The relation of r and a to the intrinsic properties of the oscillator such as the quantities k,m,w,ll in Fig. 1 is given in Sec. II. The parameter pis related to F 0 • Equations (1) govern the system if there is no impact between time tn and tn +I• Otherwise, if an impact takes place between tn and tn + I• then Eqs. (2) govern the system. Note that the Nordmark map is continuous at xn =0, but that its Jacobian matrix of partial derivatives is singular _at xn =0 [in particular, axn + ,Jaxn = -1 /(2y Xn) for Xn > 0]. This singularity at xn =0 is responsible for the new bifurcations studied in this paper. The map is normalized so that for fixed r and a, the long-time behavior is such that the orbit does not impact with the wall at Sc for p < 0, is in the grazing state for p=O, and may impact with the wall at Sc for p>O.

Thus if we vary p through zero with fixed r and a, the Nordmark map describes the dynamics of an orbit in the neighborhood of the grazing state if lpl << 1. Since the map is obtained by expansion of solutions in the neighborhood of the grazing state, its dynamics is related to the physical system only for lp I << 1. However, since we are interested in the bifurcations at p=O (i.e., the grazing bifurcations), the map is expected to capture the universal properties of impact oscillators near grazing. That is, other, physically different systems, when suitably normalized and expanded about the grazing state, should also yield Eqs. [START_REF] Moon | Chaotic vibrations[END_REF] and [START_REF] Thompson | [END_REF].

In what follows we shall be concerned with the bifurcation phenomena for the Nordmark map that occur as the bifurcation parameter p is increased through p=O (grazing incidence) with r and a held fixed. Depending on the values of r and a ( 0 < r < l' a < l + r for physically admissible systems), we observe three basic bifurcation scenarios listed as cases l-3 below. One of our goals will be to give an analysis to delineate the ( r, a) parameter space into regions in which the bifurcations in each case take place. solely by Eqs. [START_REF] Moon | Chaotic vibrations[END_REF]. In terms of the system in Fig. 1, this period-1 attractor of the map corresponds to a forced periodic orbit where the mass never impacts the wall. Referring to Fig. 2(a) we see that for subcase (a), as p is decreased from positive values, we encounter windows of stable periodic behavior, and each such window is followed by a band of chaos and then a window of stable periodic behavior whose period is one higher than the period in the previous window. Asp decreases, there is an infinite cascade of such windows of ever decreasing width in p and ever increasing period, accumulating on p=O+. To make this phenomenology clearer we plot again in Fig. 3(a) the bifurcation diagram for the same values of ( r, a) as in Fig. 2(a), but now using the variables

x /p vs lnp. We clearly see in this figure that there are six successive period addings with period 3 occurring on the right of the figure and period 9 occurring on the left. Numerically, we find no evidence of any stable periodic orbits other than those in the reversed period adding cascade. Currently we believe that the p intervals between a period m window and a period m + 1 window are occupied entirely by a chaotic attractor. Now refer to Fig. We see three successive period addings in this figure, with period 2 occurring at the right of the figure and period 5 occurring at the left. We have derived a scaling rule for the widths of the periodic windows in terms of r and a, applicable to both subcases (a) and (b). The stable periodic orbits in our period adding cascades are numerically observed to be of a very special type. In particular, if the period of the orbit is m, then the orbit spends one iterate in x > 0 and the other m -1 iterates in x < 0. We call such a periodic orbit maximal. In terms of the system in Fig. 1, a maximal periodic orbit of the map corresponds to a forced periodic orbit where the mass impacts with the wall exactly once per period.

Case 2: Bifurcation from a stable period-I orbit in p < 0 to a chaotic attractor as p increases through zero. An example of a bifurcation diagram for this case is shown in Fig. 4(a). We see that as soon asp is increased through zero (corresponding to the occurrence of impacts in Fig. 1), chaos appears. Numerically, we find for Fig. 4(a) that there is no evidence of any window of stable periodic behavior throughout the entire range between p=O and the value of p at which the stable period-2 orbit first ap- pears. In general, case 2 is defined as follows: as p increases from zero there is an interval of p values occupied entirely by a chaotic attractor, and this interval terminates at the appearance of a periodic orbit of some period M 0 • In Fig. 4(a), M 0 =2, but other values of the period M 0 occur depending on the values of r and a.

Figure 4(b) shows a case where M 0 = 3. Indeed, we observed numerically that M 0 ~ oo as the boundary in ( r, a) space between case 1 and case 2 is approached from the case 2 side. Case 3: Collision of an unstable period-M maximal orbit (which is a regular saddle, and is created, together with a stable period-M maximal orbit, in a saddle-node bifurcation in p < 0) and the period-] orbit at p = 0. When plotting a bifurcation diagram, the regular saddle, of course, does not show up. One observes that the attractor is a stable period-1 orbit for p < 0 and it is a stable period-M maximal orbit (which is created in the saddle-node bifurcation) for p > 0. Loosely speaking, we will say that there is a (discontinuous) "bifurcation" from a stable period-1 orbit to a stable period-M maximal orbit as p increases through zero. To explain the basic phenomenology of this case, imagine that the orbit is initialized on the period-1 attractor for some negative value of p, and p is then increased very slowly with time. While p remains negative, the orbit tracks the location of the period-1 orbit since the period-1 orbit is attracting for p < 0. However, when p increases through zero, the period-1 orbit becomes unstable and the orbit goes to some other attractor away from the period-1 orbit. We find that this other attractor is always a stable period-M maximal orbit, which is created in a saddle-node bifurcation in p < 0. The unstable period~M maximal orbit created in the same saddle-node bifurcation collides with the period-1 orbit at p=O. Furthermore, we :find that at p=O, depending on the parameters ( r, a), there exists either only one stable maximal periodic orbit or two stable maximal periodic orbits. When two stable maximal periodic orbits coexist, their periods differ by l. In the cases where two stable maximal orbits coexist, it is always the maximal orbit of lower period that the orbit goes to from the period-1 orbit as p increases slowly from negative to positive values. We call this the "observed" maximal orbit and we say that the period-1 orbit "bifurcates" to this observed maximal periodic orbit as p increases through zero. Figure 5(a) shows a bifurcation diagram for ( r, a) in the region where only a single period-3 stable maximal orbit exists at p=O (this is typical of what happens for other periods). We see that the period-3 stable maximal orbit is born in a saddle-node bifurcation at some negative p value, p = p 3 < 0. (The location of the period-3 saddle is indicated by the dashed lines in the figure.) For p 3 < p < 0, the stable period-1 orbit coexists with the pair of stable and unstable maximal period-3 orbits created at p=p 3 • As p~o-, the unstable period-3 maximal orbit collapses onto the period-1 orbit. The stable period-3 maximal orbit continues to exist in p > 0 and the period-1 orbit becomes a flip saddle in p > 0. In addition, we want to point out that the period-3 maximal saddle and the period-1 orbit are involved in the local bifurcation that occurs at p=O, while the stable period-3 maximal orbit is not (since it is bounded away from the origin). On the other hand, for p > 0, the solutions will converge to the stable period-3 maximal orbit that is created at p 3 • Therefore, we call the bifurcation a "bifurcation" from a period-1 attractor to a period-3 attractor. Figure 5(b) shows a bifurcation diagram for ( r, a) in the region where period-3 and period-4 stable maximal orbits coexist at p=O. Now two stable maximal orbits are created in p < 0, the period 3 in a saddle-node bifurcation at p = p 3 < 0, and the period 4 in a saddle-node bifurcation at p=p 4 <0, where p 4 <p 3 • Both stable maximal orbits continue to exist asp becomes positive, but, as already discussed, only the period 3 will be observed to bifurcate from the period-one orbit as p increases through zero. Later on in Sec. VA, it will be explained why this bifurcation to the lower period orbit is observed.

It should be noted that in all three cases above, the stable period-1 orbit that exists in p < 0 becomes a flip saddle in p > 0. That is, suppose 11 and K are eigenvalues A stable period-3 maximal orbit and an unstable period-3 maximal orbit are simultaneously created in a saddle-mode bifurcation at p = p 3 < 0. The dashed curves indicate the locations of the unstable period-3 orbit. The stable period-1 orbit and the stable period-3 orbit are on the solid curves. One can see that the unstable period-3 orbit collapses onto the stable period-1 orbit atp=O. (b) Bifurcation diagram for (y,a)=(0.80,0.55) and -? = 1. A pair of stable and unstable period-3 maximal orbits are created in a saddle-node bifurcation at p = p3 < 0, and a pair of stable and unstable period-4 maximal orbits are created in a saddle-node bifurcation at p=p4 <0. The unstable periodic orbits are not shown.

The region of ( r, a) space corresponding to systems with non-negative friction [JL ?:: 0 in Fig. 1 and v?:: 0 in Eq. ( 3)] is shown in Fig. 6, where the parameter values corresponding to the various cases in Figs. 2-5 are labeled as points. (The region shown shaded is unphysical and corresponds to negative friction.) As shown subsequently, the requirement of positive friction leads to the restrictions 0 < r < 1 and a < I + r. This region is divided into two parts by the parabolic curve K given by r =a 2 /4. The part above curve K (i.e., regions I and II) corresponds to overdamping (i.e., the linear harmonic oscillator that results from Fig. 1 with the wall removed is overdamped). This leads to real eigenvalues for the Jacobian matrix of the linear map in Eqs. [START_REF] Moon | Chaotic vibrations[END_REF]. The part below curve K (i.e., region Ill) corresponds to underdamped systems [or systems that have complex conjugate eigenvalues for the Jacobian matrix of the linear map in Eqs.

(1)]. Systems with ( r ,a) in region I have grazing bifurcations from a period-1 attractor to a reversed infinite period adding cascade (case I). Systems with ( r, a) in region II experience grazing bifurcations from a period-1 attractor to a chaotic attractor (case 2). Regions I and II are separated by the straight line segment a= ( t )y + f extending downward and leftward from its tangency point with the curve Kat (y,a)=(~,f) (see Fig. 6). stable period-M maximal orbits exist at p=O for M = 3, 4, 5, 6 with -? = 1. As already noted, when two such orbits coexist, only the one of lower period will be observed to bifurcate from the period-1 orbit with slowly increasing p. Figure 7(b) [obtained by assigning the overlap regions of the (y,a) space in Fig. 7(a) to the lower period] shows regions for which the observed bifurcating orbit has period M. Regions corresponding to higher M appear in a similar way and accumulate on the curve K as Moo. It will become clear in Sees. IV and V that the delineation of the regions in Fig. 6 is valid for all 0 ~ ?-~ 1, while the results presented in Figs. 7 is obtained with?-= 1.

Nordmark [9,10] has previously discussed scaling for case 1 and obtained case 2. The existence of the two subcases within case 1, like our treatment of the existence and stability of maximal orbits for cases 1 and 2 (see Sec. IV), is new. All the results reported for case 3, and the delineation of the ( r, a ) parameter space corresponding to each case, are also new.

In Sec. II, explicit relations between y, a, and the physical parameters of the model (Fig. 1) are obtained. Section III derives expressions for maximal periodic orbits. Sections IV and V contain an analysis of the Nordmark map for lp I << 1. This includes the existence and stability conditions of maximal periodic orbits for all ( r, a) located in the physically allowed regions in Fig. 6. The analysis results in the division of ( r, a) space into regions corresponding to the different types of grazing bifurcations in the system as explained above, as well as a scaling law for widths of the windows with high periods. Special attention is devoted to the limiting behavior as one approaches the boundaries in Fig. 6. This will, for example, show how an infinite period adding cascade results as one approaches the boundary to region I from one of the other regions. In Sec. VI, we discuss the results of this paper on grazing bifurcations in view of some general results on border-collision bifurcations obtained by Nusse and Yorke in [13]. Conclusions are presented in Sec. VII.

II. RELATION BETWEEN PHYSICAL QUANTITIES AND PARAMETERS OF THE NORD MARK MAP

In this section we study the relation between the parameters r and a and the physical parameters of Fig. 1, namely the mass m, the spring constant k, the frequency of the external forcing {t), and the friction coefficient J.L• With these expressions we will be able to understand the physical meaning of results obtained from our analysis of the Nordmark map, which is in terms of y and a. For the physical system in Fig. 1, the equation of motion without impacts with the wall at 5, is

d 2 ~ +v!!l_ +0 2 5=F 0 sin21Tf, (3)
dT dT

where we have introduced the quantities v=21TJ.L/m{t), From now on, we assume that ~-40 2 *0. Then we also have the time derivative of 5 (i.e., velocity of the mass) from Eq. ( 4):

0 2 =41T 2 k/m{t)
. • slr s21 5=P+C1s1e +C 2 s 2 e .

Hence for T=n, the state vector in the (5,~) space is and for T=n + 1, we have

I s< n + o j I 1 1 ] ~(n+1) = St Sz
where B= s e sl

I [ s<n l ] [P] =B ~(n) +(1-B) p , (5)
The matrix B has the same set of eigenvalues as the Jacobian matrix (6) of the linear map in Eqs. [START_REF] Moon | Chaotic vibrations[END_REF]. We denote the eigenvalues of matrices A and B by A. 1 and A. 2 with I A. 1 1 2: I A. 2 1. From

(5) we have

s s At= e 1' Az = e 2 . ( 7
)
From (6) we have

A. =a+v'~ I 2 (8) _a-v'~ Az- 2 (9)
Combining these relations, we obtain explicit expressions of the parameters r and a in terms of the physical parameters (10) a=A.t+A.z=e'l+/2=2e-vl2cosh [ v'v2~402]. (11) For positive friction v > 0, we have from (10) and (11) that O<r<1, O<a<1+r. (12) This also yields I A. 1 1 < 1 by ( 8) and ( 9) and corresponds to the unshaded region of the ( y, a) space in Fig. 6. Points on the curve K in Fig. 6 satisfy the relation a 2 -4y =0 [or, equivalently, v-40 2 =0 by Eqs. (7)- (11)] and correspond to systems with critical damping. Points above the curve K correspond to overdamped systems (i.e;, systems with real eigenvalues A 1 and A 2 ) and points below the curve K correspond to underdamped systems (i.e., systems with complex conjugate eigenvalues A 1 and A 2 ).

Also notice that y is related to the friction coefficient by Eq. (10). In the limit of large friction coefficient, v~ 00' we have r ~o. and the two dimensional map in Eqs. ( 1) and ( 2) reduces to the one dimensional map studied in [11], { axn +p for xn ::SO,

Xn+I= - -v'xn+P forxn>O.
Case 1 [along with subcases (a) and (b)] and case 2 were found [ 11] to occur for this one dimensional map for the a value ranges evident by examining the a axis (i.e., r =0) in Fig. 6 [i.e., case l(a) occurs for t <a< f; case l(b) occurs for 0 <a< f; and case 2 occurs for t <a< 1 ].

In the opposite limit of zero dissipation (i.e., v=O and -?= 1 ), the map given by Eqs.

(1) and ( 2) becomes area preserving. This case has been studied in [12].

III. MAXIMAL PERIODIC ORBITS

We study the grazing bifurcations at p=O for Eqs. [START_REF] Moon | Chaotic vibrations[END_REF] and (2) in the physically admissible region of the ( r, a)

parameter space as characterized in (12). For all these values of ( r, a), the system has a stable period-1 orbit for small negative p values, which becomes a flip saddle for small positive p values.

Our numerical experiments indicate that only one type of stable periodic orbit is involved in the bifurcations at p=O. We call such orbits the maximal periodic orbits. Here a maximal periodic orbit is a periodic orbit for which exactly one point per period is in the region x > 0. Our strategy is to find the range of p values in which a period-m maximal orbit exists and the range of p values in which the same orbit is stable.

Let ( x n ,y n ) represent a point on a trajectory of the N ordmark map. We use the notation (13) Then Eqs. ( 1) can be written as For a maximal orbit of period m with trajectory points ~1 .~2, ••• • ~m• we assume x 1 >0 so that x 2 ,x 3 , ••• ,xm are negative and ~m +I =~1 • Thus ~k+I= Ak-1~2+(1+ A+ A2+ ... + Ak-2)p, [START_REF] Nusse | Dynamics: Numerical Explorations[END_REF] where k = 2, 3, ... , m. The eigenvalues of A are given in Eqs. ( 8) and ( 9). Select, for example, the vectors and as eigenvectors corresponding to the eigenvalues A 1 and A 2 , respectively. Write S= [ -~2 - 1 Al] .

-r-?xi

Eq. ( 14) can now be written as a (2) k +1 is related to points on the period-m maximal orbit via Eq. ( 19).

[ Af-1 0 I -~ +yi +p ~k+I=S o A~-1 s-I -r-?xi +S'I'<k -2)g-Ip or [ Af-l 0 ] -~ +yi s-~~k+I= o A~-~ s-I -r-?xi +'l'(k-l)g-Ip.

IV. ANALYSIS IN THE CASE OF REAL EIGENVALUES OF EQS. (1)

In this section we assume that ( y, a) lie in the region where the eigenvalues At and A 2 are real (i.e., in the unshaded region above curve Kin Fig. 6).

By dividing both sides of (20) by Ar(k-u and using the notation &'kl=a/Ar(k-ll, with a standing for any variable, Eq. (20) takes the form

a _(l)(k) k+! a -(2)(k) k +I 0 0 r ~: r-1 (21)
From Eqs. (8) and (9), we have 0 < A 2 < A 1 < 1 for points that are not on the curve K. Thus (A 2 /A 1 )k and (A 1 )k both approach zero ask goes to infinity. Also in the expression of\P'k -I l [cf. Eq. ( 17)],

1-Ak For the period-m maximal orbit we have x~m~t =ximl, then Eq. ( 25) can be solved, yielding (for p 2: 0)

1 tP(k-1)=1+A•+A2+ ... +Ak-1= _ _ 1 --+ --k--+ 00 ' I I I I 1-Ai 1-Ai'
Vx\ml = At Vl+4Rm-l At-A2 2 (26)
where

R = _._if_m_l_ [ At ~tA2 ] 2 m 1-a+r A
Notice that Rm and p always have the same signs; see also (12). Equation ( 25) also indicates that for large m, the period-m maximal orbit is stable if

I d-(m) I ' Xm+t At 1 - < 1 dx\ml -At-'A2 2-v'iF (27)
Substituting Eq. ( 26) for V ximl into (27), we find that the stability condition ( 27) is equivalent to the inequality

P >pst ~l_ [ At ]2(1-A )(1-A )}._2(m-ll (
(l-A.l)2(1-A.2) 2(m-1) p<pm~ 2 At • (AI-'}..2) ( 29
)
The expressions (28) and ( 29) are key results for our subsequent discussions. We see that a period-m maximal orbit exists for p < p";, and is stable for p > p";,. Let I m denote the interval [p";, ,p"; ]. The stable period-m window appears only if the interval I m exists, i.e., if p"; > p";,. From the expressions of p";, and p"; in (28) and (29), we see that the interval Im shifts to the left (but never reaches zero) for larger values of m. Thus windows of high period may appear as p---+0+. This is in agreement with the phenomena seen in Figs. 234.

There are two distinct situations: case 1 in which Im exists for all large m, and case 2 in which for every integer m > M 0 (for some threshold M 0 ) the interval I m does not exist. The first case implies bifurcations from the period-1 attractor in p < 0 to a reversed infinite period adding cascade in p > 0, and the second case implies bifurcations from the period-1 attractor in p < 0 to a chaotic attractor in p > 0. We discuss these two cases separately as follows.

A. Case 1: bifurcation from a period-1 attractor to a reversed infinite period adding cascade

The interval I m exists for all large m if p"; > p";, for all large m. Then for any period m, there is an interval Im for which the period-m window appears if p E Im. Hence there is a reversed infinite cascade of period adding windows as p---+0+. Using the expressions of p";, and p"; in (28) and (29), we find that p"; > p";, for all m sufficiently large if At < }. By Eq. ( 8) this gives the upper border of region I as shown in Fig. 6, i.e., the straight line segment, a=tr+}'

(30)
extending leftward and downward from its point of tangency with the curve K, ( r, a)= ( f, t ).

Furthermore, from (28) and (29), we can deduce a scaling law for the window widths as p---+0+. In particular,

IIm+tl =A.2 ( 31 ) lim I 1 '
where I I m I = p";-p";, (assuming p"; > p";, for all large m) is the width of the period-m window. This scaling agrees accurately with our numerical results, and it applies (for large m) to all systems with ( y, a) in region I in Fig. 6. This scaling law was also obtained in [9].

As indicated in Figs. 2 and3, there are two different types of reversed cascades of period adding windows. For the first type, the system is chaotic between successive periodic windows in the bifurcation diagram, as in Figs. 2(a) and 3(a). Numerical experiments show no evidence of stable periodic orbits for the p values between the successive maximal periodic windows. For the second type, successive periodic windows overlap, and the system presumably does not have chaotic attractors, as in Figs. 2(b) and 3(b). The first type corresponds to the case in which the neighboring intervals I m and I m + 1 have no intersection, as schematically shown in Fig. 8(a). The system is presumably chaotic for p"; + 1 < p < p~. The second type corresponds to the case in which the neighboring intervals I m and I m + 1 overlap for large m, as schematically illustrated in Fig. 8(b). The period m and m + 1 orbits coexist for p~ <p <p"; +I• and we call this hysteresis. The border between these two types of cascades is p";, =p"; +I• which by (28) and (29) reduces to A. 1 = t• By Eq. ( 8) we find that this border is given by the segment of the line a=4r+t,

extending leftward and downward from its point of tangency with the curve K, ( y, a)= ( i6, t ). This is the dashed line in region I in Fig. u k+! II as shown in Fig. 6, i.e., the border between regions I and II.

V. ANALYSIS FOR CASE 3

The matrix A has complex conjugate eigenvalues if ( y, a) falls below the curve K in Fig. 6. Let A. 1 =rei 0 =r(cosO+i sinO), A. 2 =re -;o=r(cosO-i sinO) , where r=Vy, O=sin-1 [ 4y4~a2 r/2 (33) Note that 8=0 on the curve K (given by y =a 2 /4).

In the region below the curve K in Fig. 6, we observe grazing bifurcations from a stable period-1 orbit to a stable period-M maximal orbit as p increases through zero (as described in case 3 in Sec. 1). Recall that the actual bifurcation is a collision of an unstable period-M maximal orbit and the period-1 orbit. For ( r, a) values below but very close to the curve K, the grazing bifurcations involve orbits with high periods. In particular, M--+ oo as 8--+0 (i.e., as the curve K is approached from below). When e is not small (i.e., when M is not large), we concentrate only on the local bifurcation that occurs at x =y =0 as p--+0-; while for small e (i.e., near the curve K) we are able to do more. Asp increases from zero, the grazing bifurcation to the period-M maximal orbit is either followed by chaos or by a reversed period adding cascade starting with a period M -I window. For large M, we thus also investigate the occurrence and scaling properties of stable maximal periodic orbits in p > 0. This allows us to obtain an understanding of how phenomena below the curve K match on to those above the curve K (in particular, how the bifurcation from a stable period-1 orbit to a stable period-M maximal orbit of case 3 goes over to the bifurcation to the reversed infinite period adding cascade of case 1 and the interval of chaos extending from p=O of case 2 as the curve K is crossed from below.) A. When 6 is not necessarily small Our goal in this section is to find the regions under the curve Kin Fig. 6 corresponding to different values of the integer M. In this case of complex conjugate eigenvalues, we divide both sides of Eq. (20) by r 2 <k -I l instead of by A.i<k-u. This time we use the notation a<kl=a/r 2 <k-n, where a stands for any variable. Thus Eq. (20) now takes the form

(34)
Since the integer M need not be large, we cannot make the approximation rM-1 -o. Thus the terms of order rk-l may not be dropped, and the quantity f/!~k-ll=1+A.;+A.r+ • • • +A.~-1 (where i=1,2) in the expression w<k-1) [cf. Eq. ( 17)] may not be approximated by 1 /( 1 -A;). Then, substituting Eq. (34) into Eq. ( 19), we find that for a period-m maximal orbit with orbit points (x 1 ,y 1 ),(x 2 ,y 2 ), ••• , (xm,Yml For the part of' ( r ,a) parameter space in Fig. 6, numerical computations show that the product c:c; is nonpositive. With this in mind, we discuss the two kinds of period-m maximal orbits depending on the sign of the quantity C;' as follows:

Case (i): Cx"' < 0
In this case, both solutions x\';'l and x~2l can exist.

The expressions for ~ and ~ in (40) and ( 41)

indicate that a pair of period-m orbits, corresponding to x~';'l and x\2l, respectively, are created in a saddle-node bifurcation at some negative p value satisfying l-4C;'C;'tfm>=o.

The orbit corresponding to x\';'l only exists for p <0

[since the right hand side of Eq. ( 40) is negative for p > 0] and is numerically observed to always be unstable. In particular, it collapses onto the origin as p-0-. On the other hand, the orbit corresponding to x\2l continues to exist up to some positive p value and is observed to remain stable. Figure 5(a) is an example of this scenario. We denote the period of the pair of maximal orbits created in p < 0 by M, and the negative p value satisfying Eq. ( 42) form =Mby PM• Case {ii): Cx"' > 0

In this case the orbits corresponding to xi';'> can exist only if p > 0 (so that ~ is real and positive) and are observed to be stable only for positive p values bounded away from zero. The solution x~2l does not exist in this case. Later, in Sec. V B, we consider the periodic orbits in p > 0 that correspond to V x\';'l in the limit that the curve K is approached from below.

From cases (i) and (ii) we see that a pair of stable and unstable period-M maximal orbits are created in p < 0 in systems with ( r, a) satisfying C~ < 0 and xl!' 1 < 0 for 1 ~ k ~M -1 in Eq. ( 35) for p=O. To delineate the regions of ( r, a) space satisfying these conditions for fixed M with -r2 = 1, we take a grid in the region below the curve K in Fig. 6 We find that the regions corresponding to M and M + 1 overlap. For example, in Fig. 7(a) we see that between the regions where only M =3 occurs and only M =4 occurs, there is a region where both M =3 and M =4 occur. For systems with (y,a) located in the overlap of the regions corresponding to M and M + 1, a pair of period-M maximal orbits are created in a saddle-node bifurcation at p=pM <0 and a pair of period M + 1 maximal orbits are created in a saddle-node bifurcation at Substituting the expression of () in terms of r and a given in (33) into Eq. ( 43), we have a2 r= .

4[ l-sin 2 ( 1T /M)] (44)

Thus as M ___,. oo, the regions of successively larger M accumulate on the curve K (given by y =a 2 /4).

For any fixed pair ( y ,a) in region III [so ( y ,a) is below the curve K], we make the following three claims.

Claim 1. At most two pairs of maximal orbits, with periods differing by 1, can be created at negative p values. That is, at most two regions, corresponding to M and M + 1, respectively, can have overlap.

Claim 2. When the parameters (y,a) are in the region where only one maximal stable periodic orbit is created in p < 0, that orbit is the one that will be observed to bifurcate from the stable period-1 orbit asp increases slowly through zero.

To explain what we mean by "observed" in Claim 2, assume that p is initially negative and that the orbit is initially on the period-1 orbit. Now imagine that p is allowed to drift slowly upward with time. For p < 0 the orbit will track the location of the period-1 orbit because the period-1 orbit is stable. However, when p becomes positive, the period-1 orbit becomes unstable, and the orbit will go to some other attractor. Claim 2 is that the other attractor to which the orbit goes is always the stable maximal periodic orbit.

Claim 3. When the parameters (y,a) are in an overlap region such that maximal stable orbits of period M and M + 1 are both created in p < 0, the lower period stable maximal orbit (i.e., period M) is the one that will be observed to "bifurcate" from the stable period-1 orbit asp increases slowly through zero. Thus, Claims 2 and 3 are that the regions of the ( r, a) parameter space corresponding to observed bifurcations to a stable period-M maximal orbit are as shown in Fig. 7(b), where Fig. 7(b) is obtained from Fig. 7(a) by incorporating overlap regions into the regions of lower period.

Derivation supporting Claim 1. Let f/>1,f/>2, ••• ,f/>M be points on a period-M maximal orbit, assuming x 1 > 0 so that x 2 ,x 3 , ... ,xM<O and 4>M+J=tf> 1• Then for k=2,3, ... ,M we have 4>k+I=Af/>k+p by Eqs. [START_REF] Moon | Chaotic vibrations[END_REF], where A is defined in Eq. ( 6). For p=O, 4>k+ 1=Atf>k• We have seen that the matrix A can be written as where Sand s-1 are as in Eqs. ( 15) and ( 16). The matrix is similar to the rotation matrix

. , tf>M

(where tf>/c=P-1 S-1 tf>k, k=l,2, ... ,M) that are on a circular arc with an angle () between two successive points (except the one between 4>; and 4>; ), as schematically drawn for M = 5 in Fig. 9. For a maximal periodic orbit the points f/>;,f/>3, ... , 4>M are in the left half plane, while the point tf>i=tf>M+t is the only point in the right half plane. Write (3 1 for the angle between the half line through 0 and 4>2 and the positive y axis; and write (3 2 for the angle between the half line through 0 and 4>M and the negative y axis (see Fig. 9 for M=5). Now, let 0<(3 1 <()and 0<(3 2 <fJ be given. This implies the following. If (3 1 +(3 2 ~ (), then in addition to the stable maximal period-M orbit the Nordmark map may have a stable maximal periodic orbit of period M -1 but no other stable maximal periodic orbit. Furthermore, if (3 1 + (3 2 2:: (), then in addition to the stable maximal period-M orbit the Nordmark map may have a stable maximal periodic orbit of period M + 1 but no other stable maximal periodic orbit. Therefore, the Nordmark map has at most two maximal periodic orbits, with periods differing by 1 and coexisting at p = 0. Hence y i FIG. 9. Schematic positions of (transformed) points or a period-5 maximal orbit.

at most two pairs of maximal orbits can be created at negative p values.

Support for Claims 2 and 3. From now on, U 0 denotes a suitable region that includes the origin in its interior. We observed numerically that after its birth at PM <0, the stable period-M maximal orbit and its basin of attraction (.13M) are embedded in a region U 0 • This region U 0 was originally occupied by the basin 13 1 of the period-1 orbit. In particular, for PM <p <0, 13M and 131 share the region U 0 that was occupied entirely by 13 1 before the birth of 13M. (Here, if two maximal orbits coexist, M denotes the lower of the two periods.) As p increases toward zero, the area occupied by 13M increases and the area occupied by 13 1 diminishes. Roughly speaking, as p--.0-the region that was originally occupied by 13 1 is gradually taken over by 13M. In particular, 13 1 shrinks to a finite number of curves emanating from the origin as p--.0-. Meanwhile, the region U 0 occupied by 13M and 13 1 combined, as well as the basins of attraction of other stable periodic orbits, are not significantly altered. A region U that includes U 0 may also have points that belong to basins of stable nonmaximal periodic orbits. As an example illustrating Claim 3, the point (y,a)=(0.9,0.5) with ,-2= 1 falls in the overlap of regions M =3 and M =4 in Fig. 7(a). A pair of period-3 maximal orbits are created in a saddle-node bifurcation at p = p 3 and a pair of period-4 maximal orbits are created in a saddle-node bifurcation at p=p 4 , where p 4 <p 3 <0. A stable period-7 orbit (which is not a maximal orbit) also exists. Figures 10 plot basins of different attractors for the system with (y,a)=(0.9,0.5) and ,-2=1: the grey regions are occupied by 13 4 and the black regions are occupied by 13 (b) and lO(c), the origin (i.e., the position of the period-1 orbit at p=O) is embedded in 13 3 and is bounded away from 13 7 and 2 4 • Therefore, if p increases slightly, trajectories that start near the origin will converge to the stable period-3 maximal orbit. This explains why one observes a bifurcation from a fixed point attractor to a period-3 attractor when p traverses zero. This pattern of distribution of 13 1 , 13M, and 2M+! (M =3 in Figs. 10) is typical of systems with (y,a) in the regions where both period M and period M + 1 maximal orbits are born at negative p values. Therefore only the stable period-M maximal orbit is observed to bifurcate from the period-one orbit as p increases through zero for systems with ( y, a) located in the overlap of regions M and M+l. observed to bifurcate from the period-1 orbit as p increases through zero. Also notice that .23 7 no longer exists in Figs. ll(b) and ll(c). [Increasing a to 1.0 we find that the point (y,a)=(0.9,1.0) is in the M=4 region of Fig. 7(a) (i.e., it falls above the overlap of the M = 3 and M = 4 regions). Then only a pair of period-4 maximal orbits are created in p < 0, and almost every point in the phase space belongs to .13 4 at p=O.] Therefore we conclude that in the complex eigenvalue case, the maximal orbit that is observed to bifurcate from the period-1 orbit as p increases through zero has period M if (y,a) falls between the curves Me=rr and (M -1 )e=rr [cf. Eq. ( 43)], which includes the region where periods M and M + 1 coexist, as shown in Figs. 7.

B. Limit of small (}

From the analysis in Sec. VA, we see that the maximal orbit observed to bifurcate from the period-1 orbit at p=O has period

M = Int [ ; l + 1 , (45)
where Int means integer part of. Thus M---+ oo as e-o [where e is defined in (33)]. Therefore, in the limit of small e, terms of order ,M-I approach zero and the quantity t/!~M -o= 1 +A;+ A.f + • • • +A.:' -I (where i = 1, 2) can now be approximated by 1 /(1-A;). Then Eqs. ( 35 Basins of attraction for the Nordmark map with ?= 1, p=O, and (a) (y,a)=(0.9,0.6). The grey region is occupied by the basin of the stable period-4 maximal orbit, the black region by the basin of the stable period-7 nonmaximal orbit, and the white region by the basin of the stable period-3 maximal orbit. (b) (y,a)=(0.9,0.7). The gray region is occupied by the basin of the stable period-4 maximal orbit, and the white region by the basin of the stable period-3 maximal orbit. The stable period-7 nonmaximal orbit no longer exists at this parameter value. (c) (y,a)=(0.9,0.8). For clarity, now the white region is the basin of the table period-4 maximal orbit, and the black region is the basin of the stable period-3 maximal orbit. The stable period-7 orbit does not exist for this parameter value. [Notice that Qm and p always have the same signs; see also (12).]

Using the approximated expressions above, we reexamine the two kinds of period-m maximal orbits depending on the sign of C;' (as in Sec. VA).

Case (i): C""' < 0 By Eq. ( 48), in this case sinm () < 0 and m = M, where M is defined in Eq. ( 45). Then we have from Eqs. ( 49) and (50) that

_ m _ -sinMO-Vsin 2 M0+4QM Y x (i-• - 2 sinO (51) (52)
According to Eq. (44), regions of successively higher Min Figs. 7 approach the curve K (i.e., the parabola r =a 2 /4). The approximated equation ( 46) indicates that for large M only the period-M maximal orbit can exist at p=O, so that xjf~1 >0 and Xk~\ <0 for 1 ~k ~M-1.

(The phenomenon of the coexistence of two stable maximal orbits of periods M and M + 1 at p=O is lost due to the approximations made for large M.) Moreover, with the approximated expressions ( 46), (51), and (52), we can easily show that for l~M the maximal period-M orbit corresponding to V xW> is always unstable and the maximal period-M orbit corresponding to V x\fl is always stable (which is observed numerically in Sec. VA). In particular, since sinM () < 0, -sinMO-V sin 2 M0+4QM for all negative p values (since x \~l exists only for p < 0) satisfying sin 2 M0+4QM > 0, and

[ dx if-l. 1 } -sinM () 0 < dx\Ml x\~> = 2 sinO~ -sinMO < 1 -sinMO+V sin 2 M0+4QM
for all p value." satisfying sin 2 M()+ 4QM > 0.

Notice, however, that for o-o, i.e., when (y,a) approaches the curve K from below in Fig. 6, we have c:'-o and c: -o. Then Eq. (37) implies that the period-M maximal orbits shrink to the origin, and thus do not exist. Virfi = _ _ _ _ __.:.... _ _ _ _ ..=.::::..

sinO (53)

A /{m) _ -sinmO-V sin where fi! is the stability threshold in the complex case.

where If:. is the existence threshold in the complex case.

Thus a stable period-m maximal orbit exists in the pin-

terval J m = [
Notice that for o-o, (57) is satisfied for all large m if r < 4/9. Thus for r < f, as p increases from zero, the grazing bifurcation from the stable period-1 orbit to the stable period-M maximal orbit is followed by a reversed period adding cascade starting with a period M -1 window. The stable periodic orbits in this cascade are the maximal periodic orbits corresponding to x\'fl. Hence, for ()=0, the maximal orbits corresponding to xi'f> form a reversed infinite period adding cascade in p > 0. This corresponds to points on the curve K in Fig. 6 with y < f, which is the right hand side border of region I. Moreover, the neighboring intervals J m and J m + 1 do not overlap if {J',: + 1 < p~, which is similar to the schematic illustration in Fig. 8(a), while the intervals J m and J m + 1 overlap if {J'; + 1 > p~, which is similar to the schematic illustration in Fig. 8(b). By the expressions of p~ and if:. in (55) and ( 56), the condition if:.+ 1 > p~ yields

1 sinm () 1 8 sin() l 1 2 l 12 y < 16 sin(m + l)(J = 16 cos -tan(m + l)(J (58)
For 8-o, (58) is satisfied for all large m if y < -i6• Then there is an integer N (dependent on y) such that all windows of periods m > N overlap with each other. This corresponds to points on the curve K in Fig. 6 with y < -i6• which is the right hand side border of region l(b).

Similarly, in the limit of 8-o, (57) and the opposite of (58) combined correspond to the right hand side border of region l(a), where the system is chaotic between successive windows with high periods; while the opposite of (57) corresponds to the right hand side border of region II, where the system is chaotic in an interval in p > 0 extending from zero (recall that the period-M maximal orbits collapse to the origin as e-o). Thus we see how the phenomena below the curve K approach those above K as 8-0.

VI. DISCUSSION OF THE GRAZING BIFURCATION IN VIEW OF ORBIT INDEX

In the bifurcation theory for maps, attention has been focused on differentiable maps when one or more eigenvalues of a fixed point (or periodic point) cross the unit circle. When this occurs, the nature of the fixed point changes. For example, a fixed point attractor becomes a saddle (possibly a flip saddle) or a repeller. The Nordmark map, however, is piecewise smooth and is not differentiable at x =0. In particular, the Jacobian matrix of the Nordmark map changes discontinuously at x =0 and becomes singular as

X -0+ [ 3x,. + J3x,. = 1 /( 2~)
]. This singularity is responsible for the new bifurcation phenomena studied in this paper. According to [13], the fixed point of the Nordmark map, which is an attractor located in x < 0 for p < 0 and a flip saddle located in x > 0 for p > 0, is a border crossing fixed point; and the grazing bifurcations that take place at p = 0 are examples of border-collision bifurcations. In the rest of this section, we present the precise definitions of the terms used above, introduce the border-collision bifurcation theorem obtained by Nusse and Yorke in [13], and discuss the grazing bifurcations analyzed in this paper in the context of the more general results on border-collision bifurcations.

A map is smooth if it has a continuous derivative. Here we examine maps that are piecewise smooth, and restrict our attention to those that are smooth in two regions of the plane with the border between these regions being a smooth curve. Let r be a smooth curve that divides the plane into two regions denoted by R 1 and R 2 •

We say that a map F from the phase space R 2 to itself is piecewise smooth if (i) F is continuous and (ii) F is smooth in both the regions R 1 and R 2 • Let F(•,ll)=FI-' be a one-parameter family of piecewise smooth maps from the phase space R 2 to itself, which depends smoothly on the parameter ll• where 1l varies in a certain interval on the real line. Let El-' denote a fixed point of Fl-' defined on -E < 1l < E for some E > 0. The position of E 1-' depends continuously on ll• We say Ell is a border crossing fixed point if it crosses the border r between the two regions R 1 and R 2 as 1l is varied. Assume that the crossing occurs at ll=O. A periodic orbit P is a border crossing orbit if it includes a point that is a border crossing fixed point under some iterate of the map. If, furthermore, there exists a neighborhood U of the orbit P such that P is the only periodic orbit in U at 1l = 0, then P is an isolated border crossing orbit.

For a general approach we need the concept of the "orbit index" of a periodic orbit [START_REF] Mallet-Paret | [END_REF]. The orbit index is a number associated with a periodic orbit, and this number is useful in understanding the allowable patterns of bifurcations the orbit undergoes. We say an orbit of period p is typical if its Jacobian matrix (i.e., the Jacobian matrix of the pth iterate of the map at a point on the orbit) exists and neither + 1 or -1 is an eigenvalue (of the Jacobian matrix). For typical orbits, the orbit index is -1, 0, or + 1. The orbit index is a bifurcation invariant with respect to, as in our cases, the periodic orbits that collapse onto the fixed point E ~-' as ll-o. That is, the sum of the orbit indices of the periodic orbits that collapse onto the fixed point Ell as ll-ois equal to the corresponding sum as ll-o+. Suppose a typical periodic orbit P of a map F has (minimum) period p. The orbit index of P depends on the eigenvalues of the Jacobian matrix AP of the map FP at one ofthe points on P. Let m be the number of real eigenvalues of AP smaller than -1, and let n be the number of real eigenvalues of AP greater than + 1. The orbit index I p of P is defined by lp=O if m is odd, I P = -l if m is even and n is odd , I P = + 1 if both m and n are even .

If the orbit itidex = -1, then the orbit is a regular saddle.

If the orbit index =0, then the orbit is a flip saddle. The typical orbits with orbit index + 1 are repellers and attractors and orbits with nonreal eigenvalues. The definition of the drbit index is technical when a point of the orbit lies on the boundary r since the Jacobian matrix of the map does not exist for points on the boundary.

It is unnecessary to define the orbit index on r since we consider orbits for J.L=foO.

For a moment, assume that E J. L is in the interior of the region R 1 (or the region R 2 ), and denote 71 and K for the eigenvalues of the Jacobian matrix DFJ.L(EJ.L). If neither of the two eigenvalues 71 and K is on the unit circle, then the fixed point E J. L is a flip saddle (and has index 0) if 71 < -1 <K < 1; E~' is a regular saddle (and has index -1) if -1 < 71 < 1 < K; E is a repeller (and has index + 1 ) if both l11l > 1 and IKI> 1; and El' is an attractor (and has index + 1) if both l11l < 1 and IKI < 1. (Note that E~' has orbit index + 1 if the eigenvalues are not real.) Hence, a typical fixed point is a flip saddle, a regular saddle, a repeller, or an attractor. Now we introduce the border-collision bifurcation theorem. Let the regions R 1 and R 2 , the map F J.L, and the fixed point (or periodic point) E I' be as above. Suppose there exists a number E > 0 such that (i) E 0 is on the border of the two regions R 1 and R 2 , (ii) for -E < J.L < 0 the fixed point E J. L is in the region R 1 , and its index is I 1 , and (iii) for 0 < J.L < E the fixed point E J. L is in the region R 2 , and its index is I 2 • If I 1 and I 2 are different, then (as to be stated next) some bifurcation must occur at E 0 so that the sum of the orbit indices is invariant as J.L crosses zero. The following "border-collision bifurcation theorem" is obtained in [13]:

Border-collision bifurcation theorem. For each two dimensional piecewise smootp map depending smoothly on a parameter J.L, if the index of an isolated border crossing orbit changes as J.L crosses zero, then at J.L =0 a bifurcation occurs at this point, a bifurcation involving at least one additional periodic orbit.

This result says that additional fixed points or periodic points must bifurcate from E 0 at J.L =0 if the index of E J. L changes as J.L crosses zero. Since this bifurcation occurs while the fixed point (or periodic point) collides with the border of the regions R 1 and R 2 , we call it a bordercollision bifurcation. In other words, a border-collision bifurcation is a bifurcation at a fixed point (or periodic point) on the border of the two regions, when the orbit index of the fixed point (or periodic point) before the collision with the border is different from the orbit index of the fixed point (or periodic point) after the collision.

Therefore, if the orbit index of a fixed point (or periodic orbit) is different before and after it crosses the border r, the following two things can possibly occur.

(i) There are additional periodic orbits which collapse onto E 0 as J.L~O+ and/or J.L~O-, whose individual indices are such as to make the total orbit index conserved.

(ii) One or more chaotic sets collapse onto E 0 as J.L~O+ and/or J.L~O-. (Since the orbit index for a chaotic set is not defined, conservation of the index is no longer an issue.)

The border separating the two regions in which the Nordmark map is smooth is the line x =0. The fixed point in the Nordmark map crosses the border x =0 at p=O; thus it is a border crossing fixed point. For p<O, the fixed point is in the region x < 0 and is stable (thus it has orbit index I 1 = + 1 ); for p > 0, the fixed point is in the region x > 0 and is a flip saddle (thus it has orbit index I 2 =0). We find numerically that it is an isolated saddle at p=O. Since I 1 =foi 2 , by the border-collision bifurcation theorem, there must be other orbits (periodic or chaotic) that collapse onto the fixed point at p=O and are involved in the bifurcation there.

For systems with ( y, a) in the region I (Case 1) in Fig. 6, a reversed infinite period adding cascade collapses onto the fixed point [located at (0,0)] as p~o+. Conservation of the orbit index before and after the bifurcation at p = 0 is not violated in this case, because the orbit index of the stable periodic orbits in the infinite cascade is + 1 and the orbit index of the chaotic sets in the infinite cascade is not defined. For systems with ( y, a) in the region II (Case 2) in Fig. 6, a chaotic set collapses onto the fixed point as p~o +, whose orbit index is not defined. As for systems with (y,a) in the region III (Case 3), which is below the curve K in Fig. 6, an unstable period-M maximal orbit (which is a regular saddle and thus has orbit index -1) collapses onto the fixed point as p~o-. The value of M is determined by Eq. (45). Therefore the sum of the orbit indices is invariant (and equals zero) as the fixed point crosses the border x =0. Hence the grazing bifurcations studied in this paper are border-collision bifurcations involving different types of orbits that collapse onto the fixed point as p~O+ or p~o-.

VIII. CONCLUSION

We have observed three major types of grazing bifurcations: (i) bifurcations from a stable period-1 orbit to a reversed infinite period adding cascade; (ii) bifurcation from a stable period-! orbit to attracting chaos occupying a full interval of the bifurcation parameter; and (iii) collision of an unstable maximal periodic orbit and a period-1 orbit, which is observed to be a local bifurcation from a stable period-1 orbit to a stable maximal periodic orbit. These bifurcations are "unconventional" in that they do not occur in smooth systems. Since the Nordmark map represents the dynamics of typical systems that have low-velocity impacts and that are smooth between the impacts, the bifurcations studied in this paper are expected to be universal for such systems.

  FIG. 1. Our model physical system.

Case 1 : 3 FIG. 2 .

 132 FIG. 2. (a) Bifurcation diagram for (y,a)=(0.05,0.65) and ?= l. (b) Bifurcation diagram for (y,a)=(O.Ol,0.25) and ?= l. We use the following steps to produce a bifurcation dia• gram. (i) Set p to the leftmost value Pmin in the figure. (In (a), we start with p=pmin = -0.05.) (ii) Set initial point (x 0 ,yo) to an arbitrary point. (iii) Iterate the map 10000 times without plotting anything, to eliminate transient behavior. (iv) Iterate the map another 300 times and plot the resulting 300 values of x. This is the x position of the points on the attractor. (v) In• crement p by a small amount (in this figure, p-+p+ ~),and set the new initial point (x 0 ,y0 ) to the last point produced in the last step, and return to step (iii). Continue until p reaches the rightmost value Pmar. in the figure (in (a), Pmax =0.10). (vi) If Pmaxis reached, go to step (v), except now decrease p by a small amount every time (here p-+p-~) until Pmin is reached again. Step (vi) enables us to plot the x positions of coexisting attractors. The same steps are used to produce Figs. 3-5. The numbers of iterations in steps (iii) and (iv), and the amount of increment in steps (v) and (vi) are varied for each figure.

  FIG. 3. (a) Bifurcation diagram for (y,a)=(0.05,0.65) and r = 1 for small positive p values. We plot X I p vs lnp to take a close look at the dymamics for small positive values of p. (b) Bifurcation diagram for (y,a)=(0.01,0.25) and r= 1, X /p VS lnp.

FIG. 4 .

 4 FIG. 4. (a) Bifurcation diagram for (y,a)=(O.l5,1.0) and r = 1. The highest stable periodic orbit in p > 0 has period M 0 =2. (b) Bifurcation diagram for (y,a)=(0.15,0.95) and r = l. The highest stable periodic orbit in p > 0 has period M 0 =3.

o p3 p 2 -FIG. 5 .

 25 FIG. 5. (a) Bifurcation diagram for ( y, a)= (0. 50, 0. 20) and ~ = 1. A stable period-3 maximal orbit and an unstable period-3 maximal orbit are simultaneously created in a saddle-mode bifurcation at p = p 3 < 0. The dashed curves indicate the locations of the unstable period-3 orbit. The stable period-1 orbit and the stable period-3 orbit are on the solid curves. One can see that the unstable period-3 orbit collapses onto the stable period-1 orbit atp=O. (b) Bifurcation diagram for (y,a)=(0.80,0.55) and

I 8 y 1 . 0 FIG. 6 .

 8106 Fig. 5bI

4 FIG. 7 .

 47 Fig.11c Fig. 11b Fig. 11a Figs. 10 { Maximal M=3 orbits exist 0.0 ~~_.__~~L...L.~_.__......_._~L...L.~-' in p<O 0.0 0.2 0.4 0.6 0.8 1.0 y

  every integer 1 ~ 0 we can write where (17) Since by Eqs. (2), -1-~+yl+p] ~z-

  Now let then .k =So-k, and we have [ xk l [ a~ll+a~2l l Yk = -A2a~u-A 1 a~2 l 09 ) Using these notations, we rewrite Eq. (18) as follows: (20) This is the mapping from to [ a~1~1 I a~2~1 fork= 1,2, ... , m, where [ an~]

where i = 1 , 2 .

 12 Hence for large k Eq. (21) small terms fork>> 1, -(2)(k)_ -A2 p ak+l-(At-A2)(1-A2) AI(k-Il (24) By Eq. (19), Eqs. (23) and (24) give A --~kl x(k) = -1 • / x(k) + _.£:.p _ _ k+l At-A2 v t 1-a+y (23) (25)

6 .FIG. 8 .

 68 FIG. 8. (a) Schematic illustration of the situation in which intervals Im+ 1 and Jm have no intersection. The system is presumably chaotic for p"; + 1 < p < p~. (b) Schematic illustration of the situation in which intervals Im + 1 and Im overlap.Period m + 1 and period m orbits coexist for p~ < p < p"; + 1• This is called hysteresis.

+ 1 sinO 1 sinO y 1 _ 1 _~

 111 where x 1 >0 and x 2 ,x 3 , ••• ,xm <0, we have x<kl =-sinkO -Vx\kl + sinkO ( <k-UI2)y<kl k sin( k -1 )0 ( kl2-r2)x<kl sinO r I + [ A.lf/!~k-I)_A.2f/!~k-1) }t~k), A. 1 -A. 2 ylkl = vr sin(k -1 )0 -Vx'P k +! sinO sin(k -1)0 ( k/2)rr<kl sinO r J' I + sin(k -2)0 ( <k+llll-r2>x<kl sinO r I I .t.<k -I)_ .t.<k -I) I k = 1,2, ... , m. For the period-m maximal orbit we set x~ml =x~m~1 and y\ml =y~m~1 in Eqs. (35) and (36); then we have cmx<ml+~ /x<ml +Cmp.,.i.ml=O X sin(m -1 )0 ml2-r2+ sin(m -1 )0 mt2 X sinmO sinmO r sinmO r + sin 2 ( m -1 )0 m-r2 sin( m -2 )0 m-r2 /.,.<ml _ -l-V1-4C;'C;'tfm> v ""'12 -2C;' (41) Notice that the right hand sides of the solutions [Eqs. (40) and (41)] are required to be real and positive.

  and numerically determine from Eqs. (35)-(41) the regions yielding C~ <0 and xl!l 1 <0 for 1 ~ k ~ M -1 with p =0 and -r2 = 1. Since all points in the region corresponding to M are required to yield c: < 0, and since Eq. (38) indicates that c:-+ oo for sinMO-o+ while c~--oo for sinMO-o-, the upper boundary of the region in the ( r, a) parameter space corresponding to M is given by M0=1r (43) for all 0 ::=:: -r2 ::=:: 1. Results for M = 3, 4, 5, 6 are shown in Fig. 7(a). (The conditions are never satisfied forM= 1,2.)

  p=pM+I <0, where PM+! <PM• Figure 5(b) is an example of this forM= 3.

  [ cos() -sin() J sin() cos() • That is, there exists an invertible matrix P such that [ e;e 0 ]= [cos() -sin()) _ 1 0 e -,e p sin() cosO p • Hence at p=O, ( cos() tf>k +I= Atf>k =rSP sin() -sin()] p-Is-If/> cos() k Multiplying both sides by p-Is -I, we have -1 -1 -[cos() -sin() l -I -1 (P S 4>k+I)-r sin() cos() (P S 4>k> • Thus the points f/> 1 ,f/> 2 , ..• ,tf>M can be linearly transformed to points f/>;,f/>;, ..

  7 • In Fig. lO(a), p= -0.02 <p 3 , the period-3 maximal orbits are not yet born, and the white region is occupied by 13 1 • The point marked with a + sign in Fig. lO(a) is the position of the stable period-1 orbit. In Figs. 10(b) and lO(c) [where Fig. lO(c) is an enlargement of the small box in Fig. lO(b)], p=O, the period-1 orbit has just lost its stability, and the white region is occupied by 2 3 • The three points marked with + signs in Figs. lO(b) and lO(c) are the positions of the stable period-3 maximal orbit. Notice that in Figs. lO

Figures 11 (

 11 Figures 11(a)-11 (c) are similar plots of basins of attraction as Figs. 10(a)-10(c). (Both Figs. 10 and 11 are produced with the DYNAMICS software [14].) In Figs. 11(a)-ll(c), ,-l=t and p=O with (y,a)=(0.9,0.6) in

FIG. 10 .

 10 FIG. 10. Basins of attraction for the Nordmark map with (y,a)=(0.9,0.5), r=l, and (a) p=-0.02, (b) p=O, (c) p=O. Figure lO(c) is an enlargement of the small box in Fig. lO(b). The grey regions in (a)-(c) are occupied by the basin of the stable period-4 maximal orbit; the black regions in (a)-(c) are occupied by the basin of the stable period-7 nonmaximal orbit. The white regions in (a) is the basin of the stable period-1 orbit. The white regions in (b) and (c) are the basin of the stable period-3 maximal orbit.

Fig. 11 (

 11 Fig. 11(a), (y,a)=(0.9,0.7) in Fig. 1l(b), and (y,a)=(0.9,0.8) in Fig. ll(c). The points (y,a) in Figs. 1l(a)-1l(c) all fall in the overlap of regions M = 3 and M =4 in Fig. 7(a), so both stable period-3 and period-4 maximal orbits exist at p=O. Comparing Fig. lO(b), Fig. ll(a), Fig. 11(b), and Fig. ll(c), we see that at p=O the territories occupied by .23 3 and .23 7 combined are reduced and the area of .23 3 decreases to zero as ( r, a) moves upwards approaching the upper boundary of the overlap of the regions M = 3 and M =4 in Fig. 7(a). Thus we see how M = 4 supplants M = 3 as the stable maximal orbit

1 -

 1 FIG. 11. Basins of attraction for the Nordmark map with ?= 1, p=O, and (a) (y,a)=(0.9,0.6). The grey region is occupied by the basin of the stable period-4 maximal orbit, the black region by the basin of the stable period-7 nonmaximal orbit, and the white region by the basin of the stable period-3 maximal orbit. (b) (y,a)=(0.9,0.7). The gray region is occupied by the basin of the stable period-4 maximal orbit, and the white region by the basin of the stable period-3 maximal orbit. The stable period-7 nonmaximal orbit no longer exists at this parameter value. (c) (y,a)=(0.9,0.8). For clarity, now the white region is the basin of the table period-4 maximal orbit, and the black region is the basin of the stable period-3 maximal orbit. The stable period-7 orbit does not exist for this parameter value.

  ) x!M) =2 , ()-/;;fMi 1 11 sm Y XiT = -sinM() ___ ___;;;:;:;==:===== > 1 '

Case

  (ii): C:' > 0 By Eq. (48), in this case sinmO>O; hence m <M, where M is defined by Eq. (45). Then we have from Eqs. (49) and (50) that -sinmO+V sin 2 m0+4Qm

1 Vr

 1 On the other hand, a maximal orbit with period m can exist if x 2 ,x 3 , ••• ,xm <0. Equation (46) indicates that Xm =max{xk: 2~k ~m} for large m. Hence a maximal orbit with period m can exist if Xm < 0, which means that Ym + 1 = -rxm > 0 by Eqs. (1). Using Eqs. (47) and (53), this reduces to <~x-[ 1-a+r sin 2 (m -1)() p Pm -r sin20 1-a+r sin(m -1)()sinm0 l m-

  {f~ ,{f;.] if {f;. > {f~. Using the expressions of

	ps~ and {f;. in (55) and (56), the condition {f;. >If~ yields
	4 [ sin(m -1)() ] 2 _ 4 [ () r<-=-oos -9 sinm () 9	sinO ] 2 tanm ()

ACKNOWLEDGMENTS

This work was supported by the Department of Energy (Office of Scientific Computing, Office of Basic Energy Sciences) and by the Office of Naval Research (Physics).


Marion Schroedter-Homscheidt
email: marion.schroedter-homscheidt@dlr.de.

Antti Arola

Niels Killius

Mireille Lefèvre

Laurent Saboret

William Wandji

Lucien Wald

Etienne Wey

The Copernicus Atmosphere Monitoring Service (CAMS) Radiation Service in a nutshell

The Copernicus Programme and its Atmosphere Monitoring Service (CAMS)

Copernicus, previously known as GMES (Global Monitoring for Environment and Security), is the European Programme for the establishment of a European capacity for Earth Observation [START_REF]Copernicus website[END_REF] with respect to land, marine, and atmosphere monitoring, emergency management, security, and climate change. The atmosphere service of Copernicus combines state-of-the-art atmospheric modeling on aerosols with Earth observation data to provide information services covering European air quality, global atmospheric composition, climate, and UV and solar energy [START_REF]Copernicus Atmosphere Monitoring Service (CAMS) website[END_REF]. Besides the radiation service, it provides information on -among othersultra-violet radiation and aerosol concentration on a global scale.

The CAMS Radiation Service

Within the radiation service, existing historical databases HelioClim-3 and SOLEMI for monitoring incoming surface solar irradiance have been further developed. The new service is jointly provided by DLR, Armines, and Transvalor. The Monitoring Atmospheric Composition and Climate (MACC) project series has been preparing for the service provision, which is now operational as part of the Copernicus programme. Data are made available both via the Copernicus portal and the SODA portal [3]. A User's Guide [START_REF]User's Guide[END_REF] has been created during the MACC project and will be updated on a yearly basis. The scientific algorithm is described in [START_REF] Lefèvre | McClear: a new model estimating downwelling solar radiation at ground level in clear-sky conditions[END_REF] and [START_REF] Qu | Fast radiative transfer parameterization for assessing the surface solar irradiance: The Heliosat-4 method[END_REF].

The radiation service consists of an all-sky radiation time series service taking satellite-based cloud parameters into account and a clear-sky radiation time series service for cloud-free skies. Quality of the service is ensured by regular input quality control, regular quarterly benchmarking against ground stations, and regular monitoring of the consistency in order to detect possible trends.

Following the Copernicus data policy, all data is provided free for any use after a registration giving a name and email address. It is not allowed to sell the data directly without modification, but data may be used for any purpose and value-added data may be part of any commercial usage.

Details of the all-sky radiation service are:

• Period of record: Feb 2004-present, updates are made continuously, data is provided with up to 2 days delay

• Temporal resolution: 1-minute, 15-minute, hour, day, month

• Spatial coverage: Europe/Africa/Middle East/Eastern part of South America/Atlantic Ocean.

• Spatial resolution: Spatial resolution is the original pixel of the Meteosat Second Generation image (approx. 3 km at satellite nadir and 5 km at mid-latitude).

• Data elements and sources: Global, direct, diffuse, and direct at normal incidence irradiances; global, direct, diffuse and direct normal irradiances in cloud free conditions; verbose mode with all atmospheric input parameters used for clouds, aerosols, ozone, water vapor and the surface reflective properties. The paper will summarize the new service capabilities and illustrate quality control and validation results.

Fig. 1 :

 1 Fig. 1: Web interface of the CAMS radiation service