C Toulemonde

C Gontier

Sticking motions of impact oscillators

Sticking motions of impact oscillators

ABSTRACT. -The present paper is devoted to the study of impact oscillators subjected to harmonic excitation. The interest is essentially centred on a family of behaviours which has never been systematically investigated before: the class of sticking periodic responses of impact oscillators.

A formulation of this kind of motion is presented and a Poincare application is built. A method is defined in order to produce its characteristics and an analytic differentiation of the responses is used to evaluate local stability.

A methodology of analysis, based on a Predictor-Corrector method, is presented and applied to single and multiple degree of freedom systems: bifurcation diagrams and parameter space partitionings are developed. 0 Elsevier, Paris 1. Introduction

Among the various aspects of nonlinear vibrations, the problem of impacting structures deserves special attention in view of its scientific interest and its technical usefulness.

Indeed, the specific challenge of this problem comes from the strong nonlinearity resulting from the impact. Many researchers have notice the tremendous variety of dynamic responses of these systems, even under a simple harmonic excitation [START_REF] Thompson | Nonlinrur Dynamics und Chuos[END_REF][START_REF] Moon | Experimental Study of Chaotic Vibrations in a Pin-jointed Space Truss Structure[END_REF]. As to its technical applications, the impact oscillator provides a realistic model for a very common situation in vibrations. The vibration of a structure with clearance, the response of an untightened bolted joint, the motion of a body travelling between guides are examples of such situations. In general, the dynamics of an impacting structure is so complex that most authors have restricted their analysis to the simplest model, i.e. the single degree of fredom (d.o.f.) oscillator. In this framework they were able to confirm the classic scenario of period doubling bifurcation (Holmes, 1983;[START_REF] Chow | Bifurcation of Subharmonics[END_REF]). Some specific features such as grazing bifurcations have also been identified [START_REF] Nordmark | Non-Periodic Motions Caused by Grazing Incidence in an Impact Oscillator[END_REF][START_REF] Ivanov | Impact Oscillations: Linear Theory of Stability and Bifurcations[END_REF]. A few studies were carried out using analytical approach but most studies proceeded from time integration of the dynamics [START_REF] Paoli | Analyse nume'rique cle r~ihrutions ur'bc contruinres nnik&rale[END_REF].

In spite of all this research there is a very common phenomenon which strangely enough has drawn very little interest. We will refer below to this phenomenon in terms of "sticking motion". To describe roughly this phenomenon, let us assume that the impact is approximately governed by a restitution rule. Depending on the mechanical conditions and on the dynamic state of the body, a structure may undergo an infinite sequence of decreasing impacts, although finite in time, before it remains for a time in contact with the obstacle. In the case of multi-d.o.f. oscillators, this does not mean that the structure is motionless, but during a finite time contact is maintained with the stop, as if the body was stuck to it. Moreover, this kind of motion can converge asymptotically on a periodic solution including one or several sticking phases per period.

Apparently, this phenomenon has been the object of very few studies in the past. A particular case [START_REF] Shaw | Periodically forced linear oscillator with impacts: chaos and long-period motions[END_REF] with a small or zero restitution ratio was analysed, and the domains of attraction of infinite impact pariodic motions (Budd, 1994) have been numerically analysed. However, these studies were also restricted to the one d.o.f. oscillator. To our knowledge, a general investigation of this phenomenon is still missing, and is the subject of the present paper, whose aim, in general, is to define a suitable procedure for the analysis of these systems. First, a general methodology for handling the dynamics of multi-d.o.f. oscillators will be proposed. Then, by applying the procedure to the case of sticking motions, a method to determine the periodic solutions will be defined, and their stability will be calculated. Finally, the whole procedure will be applied to a few examples, that will serve as validation tests, and as a means of understanding some specific features of the motions.

2. Dynamic equations of an impacting structure

I. THE MECHANICAL SYSTEM

In this paper, the system under consideration is a multiple d.o.f. impact oscillator. This oscillator has a linear structure that can impact on a single rigid stop and it is subjected to harmonic excitation (Fig. 1). The type of nonlinear displacment constraint is determined by the position of the stop and by its orthogonal distance 20 to the equilibrium position of the structure.

When the displacement constraint :cg is positive the static system is always tensioned. For small excitation magnitudes this preload constrains one degree of freedom of the f sin(ot) rigid stop A/SOLIDS. VOL. 17, No 2. 1998 structure onto the stop. In the case of a multi d.o.f. structure, linear harmonic motion continues while contact is maintained with the stop. When the displacement contraint ,~a is negative clearance exists in the system. In both cases, for small excitation magnitudes the structure behaves linearly and nonlinear behaviour does not exist. The dynamic behaviour changes when the stop is reached by the oscillator.

For such a dynamic system, a very large variety of asymptotic solutions can be encountered: finite impact periodic motions, infinite impact periodic motions, quasiperiodic motions, chaotic motions (Gontier, 1997). Infinite impact periodic motions and sticking motions in general are the subject of the present study.

Referring to Figure 1, let us note two major points of the system, A and B. Harmonic excitation is applied at point A in a constant direction. The stop is impacted by the structure at point B. As zn is small, it is assumed that the contact point B on the surface is constant and that the reaction of the stop has a constant direction orthogonal to the structure at B.

DESCRIPTION OF THE STICKING MOTIONS

Let us define more accurately the so-called sticking motions. For such dynamic responses, contact at point B is maintained during a finite time per period. This dynamic phenomenon may result from different situations. In the case of a purely plastic shock, i.e., a shock involving no rebound, the post-impact velocity at point B is set to zero. In that case, if the incident acceleration is such that the structure impacts with the stop, i.e. if the acceleration is negative, the oscillator will be stuck on the stop. This phenomenon can be observed when the impact rule is a restitution rule whose coefficient is zero.

However, in most impacting systems, the post-impact velocity follows an elasto-plastic impact rule. Generally, this rule can be approached by a restitution rule with a non zero restitution coefficient. Even in this case, sticking periodic motions can be encountered, but now they are observed after an infinite sequence of decreasing impacts, although finite in time. Once the limit of the sequence is reached, it can be said that the structure sticks onto the stop. Such dynamic behaviour frequently occurs in impact oscillators, espacially in the preload case, where it is encountered for rather small values of excitation magnitude.

THE FREE FLIGHT PHASE

During infinite impact periodic motions, the dynamics of the structure undergoes an infinite sequence of free flight phases. For the sake of simplicity, let us write the dynamics directly in the eigenbase of the linear structure. The equation of motion of a free flight phase is therefore The vectors and the matrices of Eq. ( 1) have an arbitrary dimension n. Therefore, the theory is applicable to any multi-d.o.f. structure.

A part of the kinetic energy is lost when the structure impacts the stop, which happens when 1:~ (t) is equal to ~(1. In the present study a restitution rule is applied in this case (2) where T is the restitution coefficient, whose value lies between zero and unity. This formulation is equivalent to assuming that the impact duration is very small and that the reaction of the stop is an impulse. During sticking motions, and when the restitution coefficient is not zero, the sticking phenomenon will occur after an infinite sequence of impacts. theoretically, an infinite number of free flight phases will be encountered until the sticking conditions are reached.

THE STICKING PHASE

Sticking occurs when the position xB(t) of point B is equal to 20, its velocity :iu (t) is equal to 0 and the contact reaction of the stop is positive. This situation may be the consequence of a purely plastic impact, or of an infinite sequence of decreasing impacts. Stucking ends when the reaction of the stop is null and when it changes sign.

We denote the reaction of the stop as R(t). The sticking phase motion equation is . .

X + RX = T;lf sin(d) + TBf2(t):

In Eq. (3) TB is the projection of a unit force applied at B onto the eigenbase.

During the sticking phase, the contact reaction R(t) can be easily calculated. The position and the velocity at point B remain constant during this period, therefore acceleration at point B is zero and the calculation of reaction R(t) results directly from the expression Zg(t) = 0. Let us first assume relation between scalar displacement n:~(t) and displacement field X(t):

(4) XB((t) = X;X(t)

In Eq. (4) TB is a vector and (') is the transpose operator. In the same manner, the equation of motion (3) can be projected and the reaction of the stop can be extracted: The result (5) will be used to modify the equation of motion (3) and to determine the end of the sticking phase when the reaction is zero.

The eigenbase of the structure being linear, vector Xg can be chosen such that II x:, = [l.. . I], i.e. such that zo(t) = xxi(t) EUROPEAN JOURNALOFMECHANICS.A/SOLIDS,VOL. 17, No 2. 1998

RESTRICTED EQUATION OF MOTION OF THE STICKING PHASE

During the sticking phase, one of the degrees of freedom of the structure is constrained by the stop. The dynamics of the constrained system is ruled by a (n -1)-sized dynamic equation. First, by using Eq. (3) and Eq. ( 5) it is possible to give a new n-sized formulation of the dynamics:

(6) Eq. ( 6) only holds when the reaction of the stop is positive. Once the dynamics described by Eq. ( 6) are integrated, the reaction can be calculated with the help of expression (5). As shown by this relation, the reaction of the stop is a linear function of the displacement coordinates X. Therefore, in this example the stiffness matrix is modified by the expression of the contact reaction and its eigenvalues are the square of the eigenpulsations of the linear structure constrained at point B.

In order to achieve the complete transformation of the dynamic equation a number of other conditions must be added. During the sticking phase, only n -1 components of the n-sized vector X are free. The last coordinate of X can be expressed as a sum of the first n -1 components: I, -1 (7) x:o(t) = 20 e X(n) =z(jcx(4 i=l

We denote the vector of the first n -1 coordinates of vector X as Y. By using this new vector, the n-sized equation of motion (6) can be reduced to the following (n -l)-sized equation The last two subsections were devoted respectively to the two equations of motion, Eq. ( 1) and Eq. (3). Now the impact rule (2) will also be expressed in detail. In the case of a multidimensional structure the restitution rule has to be applied carefully and relation (2) cannot be translated directly using the velocity field 2.

The stop is rigid, therefore the shock stiffness is very high and the impact duration may be very short, as described in Eq. ( 2). The reaction can be approximated by an impulse, and the relation between the pre-impact and post-impact velocity fields can be deduced from the time integration, in the sense of distributions, of the acceleration field during the shock:

(10) A@+) = k(t-) + p . TB In Eq. ( 10) p is the magnitude of the impulse to which the structure is subjected, which can be calculated by taking the impact rule (2) Eq. ( 10) and relation (4) into account: Then, the relation between the two velocity fields results from Eq. ( 10) and Eq. ( 11).

The multidimensional restitution rule can now be rewritten as This paper is mainly focused on the study of sticking periodic motions. Depending on the value of the restitution coefficient, two different situations can be encountered. First, when the restitution coefficient is zero, i.e. when the impact rule is a plastic impact rule, each free flight phase is followed by a sticking phase (Fig. 2). The motion is an alternate sequence of free flight phases and sticking phases.

xB Cm) -~. When the restitution coefficient is not null, Le. when the impact rule is elasto-plastic, the motion encounters an infinite number of impacts before any sticking phase (Fig. 3). From a pratical point of view, it is imposible to calculate this infinite amount of impacts, but it is possible to approach it as the limit of a geometrical series.

I 3.6E-4 3.6E-4 -~~~ I 3.4E-4 3.4E-4 -_ / 3.2E-4 3.2E-4 - I 3.OE-4 3.OE-4 7 I 4.OE-2 4.OE-2 6.OE-2 6.OE-2 S.OE-2 S.OE-2 1 1 .OE-1 .OE-1 1.2E- 1.2E-1 1 1.4E-1.4E-1 1 1.6E- 1.6E-1 1 time (s) time (s)
XB Cm)

DISCRETISATION OF THE FREE FLIGHT DYNAMICS

The equation of motion (1) of free flight motion is deterministic, thus its continuous dynamics can be completely defined by the initial conditions. Expressed in the structure ElJROPEAN JOURNAL OF MECHANICS.

A/SOLIDS. VOL. 17. No 2, 1998 eigenbase, the uihole displacement and velocity fields between two successive impacts (impact number j and impact number j + 1) can be written as the sum X(t) = U(t) + fJCJ2j.j COS(LIJj(tt,j)) + b;,,, sin(wj(t -t,))]X; i=l In Eq. ( 13) U(t) is the forced harmonic displacement field and tj is the time of impact numbered .j. Set (Xj)jE[l.rr] is the set of eigenvectors of the free structure.

The motion defined in Eq. ( 13) depends only on t,j and on the set of 2n variables (ai,,j, bj,j)jclr.,l.

But tj is the time of impact j and it must satisfy the relation "B(t)) = XLX(tj) = :cg and thus depends on the variables (ai,j: bi,,j)jcl1.,,]. Therefore, only 2n variables of the 2n + 1 variables describing the dynamics in Eq. ( 13) are free. The dimension of the free flight dynamics is then 2r1. Let us define Vj as the 2r1 dimensional degree of freedom vector describing the free flight phase numbered j: The free flight continuous motion dynamics can then be reduced to a discrete dynamics relationship, equivalent to a recurrence relation between vector Vj and vector ?,+I. Vj+l can then be expressed explicitly using vj and the two successive impact times t,j and t,)+l : EUROPEAN JOURNAL OF MECHANICS, A/SOLIDS, VOL. 17, No 2, 1998

DISCRETISATION OF THE STICKING PHASE

In the same manner, the dynamics of the sticking phase can be discretised. The equation of motion (8) of the sticking system is also deterministic and its dynamics can also be completely defined by the discrete initial conditions. We denote the set of the 7~ -1 eigenvectors of the linear homogeneous dynamic system of Eq. ( 8) as (k;:)iEll.,,--l~; we also denote the set of eigenpulsations of this system as (~j);~lr.,~-rl. During the sticking phase, the restricted displacement and velocity fields in the structure can be written as the sum as an implicit function of the set (e;. j';)i~[l,,l-lI. Consequently, the state of the system at the end of the sticking phase can be completely described by the following vector: In the same fashion, taking expression ( 18) and relation R(ti) = 0 into account, it can be proved that WI is an implicit function of WO. In the case of a one dimensional impact oscillator, Wo is reduced to t{, and WI cannot be defined, but the problem can be solved very easily. In that case, the time, displacement and velocity fields at the end of the constrained phase are defined by an implicit and constant function of ti.

Y(t) = YO + U'(t) + C[c; cos(wj(t -tb)) + di sin(wj(t -&))]X;

FIXED POINT FORMALISM

In the last two subsections, the two continuous dynamics of the system, in its free state and in its constrained state, have been discretised. Now, the complete sequence of motion has be discretised. Let us denote F as the function which describes the free flight motion The explicit expression of F, as given by Eq. (151, is in fact an implicit function. In the same manner, the sticking phase as described by Eq. ( 20) can be denoted as function G:

(22) Some other functions must be defined in order to complete the transition the degrees of freedom of the two different types of dynamics. We first describe the relation between EUROPEAN JOURNAL OF MECHANICS, A/SOLIDS, VOL. 17, No 2, 1998 the degrees of fredom at the end of sticking phase WI and the degrees of fredom of the first free flight phase VI: ( 23We also define the relation between the degrees of freedom at the end of last free flight V, and the degrees of fredom at the beginning of the next sticking phase ( 241The complete sequence can { Hz : @" --) jJih-1 Kc ++ wo 1 be defined by linking the last four described functions. Let Wt,j be the vector describing the state of the structure at the end of the sticking phase numbered j. The dynamics between the beginning and the end of two successive sticking phases, i.e. between W1.j and Wl.,j+l, can be described by a implicit function 25) application P appears as a composed function, and defines a Poincare map associated with sticking motions. Every continuous orbit of this kind can be made equivalent to the series (Wl.j),jc-v of restricted degree of fredom vectors. As defined above, a Poincare map P is a very convenient representation of the dynamics, and is well adapted to the analysis of the asymptotic behaviour of the system. For instance, an infinite impact periodic motion with only one sticking phase per period can be defined by the following fixed point equation: where WI defines the restricted degree of freedom vector at the end of each sticking phase. In the case of a one dimensional impact oscillator Eq. ( 25) and Eq. ( 26) cannot be defined. As mentioned in subsection 3.3, it is possible to discretise the sticking phase, i.e. to calculate HI o G. The Poincare application P can then be replaced by Hz o FX o HI o G.

RESOLUTION METHOD

The simplest method to determine stable periodic motions is based on iterative integration of the dynamics until stabilisation is found. It is well known that the efficiency of these methods depends heavily on the stability of the motions and the choice of the initial conditions. Although sticking periodic motions are generally stable, the method is in any case expensive and often not very accurate. A more direct method based on analytic resolution is proposed here. Using the fixed point formalism introduced in the last subsection, the problem of periodic motions reduces to finding the roots of a Boundary Value Problem

ANALYTIC APPROACH TO INFINITE SERIES OF IMPACTS

In Eq. ( 25), the Poincare map appears as a chain of several functions. In the case of an infinite impact periodic motion the theory implies that application F will be iterated an infinite number of times. Pratically, since this calculation cannot be carried out exactly, the problem must be simplified in order for it to be solved.

Qualitatively, in the case of infinite impact periodic motions, it can be shown that the sequence of impacts converges towards a geometric series. The ratio of this series is equal to the restitution coefficient T. Let (tj)je[l.+K,[ denote the series of impact times in this infinite sequence of impacts. The series of impact times and the series of impact velocities satisfy the following relations:

Series (Vj)jc[l.+x[ of the degrees of freedom satisfies relation (15). Using this recurrence relation and the asymptotic geometric behaviour described in Eq. ( 28) the limit V,, of The new application defined by equation ( 30) will be used in place of the original function.

EUROPEAN JOURNAL OF MECHANICS.

A/SOLIDS, VOL. 17. No 2. 1998

Another way to solve this problem is to suppose a purely plastic impact rule as soon as the impact velocity is smaller than a given value.

FIRST ORDER DIFFERENTIATION TERMS

As mentioned in section 3.4, equation ( 25) can be reduced to a root finding problem, (26). The roots of this system are found by applying a Newton algorithm. Although the whole nonlinear application is quite heavy to handle it is possible to calculate its first order derivatives, i.e. its gradient, in an analytic manner. The first order derivative of 4 is composed of five different linear applications, DHr, DHz, DF, DG and DL.

The five implicit functions have been analytically differentiated with respect to their inputs. In the case of DHI and DHz the calculation in very simple. The expression of function F is given by Eq. ( 15) and its differentiation requires the calculation of a number of first order terms, such as dtj and dt,;+l, that can be performed by using the relation ( 31) "B(t,j) = XB'U(t,j) + [I 0 ". 1 ()]I/; = X(J By differentiating Eq. ( 31) the first order term dtj appears as a linear function of dyi. DG and DL are calculated in the same manner.

ANALYTIC AND EXPLICIT RESOLUTION IN A SIMPLE CASE

In the simple case of a one degree of freedom impact oscillator, when the restitution coefficient T is null, the characteristics of the periodic sticking motions can be easily calculated. For the sake of simplicity, we consider a non dimensional impact oscillator, i.e. in the case where zu and the first eigenpulsation are equal to 1.

By definition time ti is the time at the end of the sticking phase, but it is also the time at the beginning of the next free flight phase. ti thus satisfies the following conditions:

With the help of the previous relation it is easy to predict that no periodic response with sticking phases exist for values of amplitude f less than 1 because in that case the trigonometric equations would not have any real solutions.

Taking these initial conditions into account, the free flight phase motion can be explicitly integrated. Time tb then satisfies the condition

(33) z&) = -f$ cos(tb -t'l> -wm sin(f/o -() + ' 1-J ' ' 1 sin(&) = 1
EUROPEAN JOURNAL OF MECHANICS.

A/SOLIDS, VOL. 17, N" 2, 1998

13 A new sticking phase at time ti,. As a consequence, tb must also satisfy the relation .Y(t[,) < 0. In the two following limit cases, by developing Eq. ( 33), the root of Eq. ( 33) can be calculated explicitly: In general Eq. ( 33) can be used to determine the zone of existence in parameter space of a large class of sticking periodic motions. In the case of a nondimensional single degree of freedom impact oscillator with zero restitution coefficient, subharmonic periodic responses can be encountered. These responses occur when solution $ of Eq. ( 33) is such that t[, -"i is larger than the period z (Fig. 4). In that case, a subharmonic coefficient g of the subharmonic sticking perioWdic responses is introduced. In the following equation we denote the smaller amplitude for which these subharmonic orbits exist as &,icking. The responses corresponding to this limit are such that (35) Taking condition (33) and equation ( 35) into account, the critical amplitude isticking can be implicitly calculated:

(36) J f&king -I When the subharmonic coefficient q is equal to 1 the value of fsti,-king is constant and equal to 1. The limit determined by this expression has asymptotes in the parametric plane defined by the two parameters (w, f):

(37) w------+zq-1 fsticking -t+x .&ticking __i d-+-X f such that 4@-Z = (2q -1)~ -2 arcsin

In the same manner, the larger amplitude for which these subharmonic responses exist can be calculated analytically. We assign a value j'free-sis+t to this amplitude, which is determined by introducing relation tb -ti = F in equation ( 33):

(38)

This expression also has two limits :

(39)

When the restitution coefficient 1' is not null but small, the conclusions about the existence of the subharmonic responses are always valid, as the duration of the infinite series of impacts is very small. When the restitution coefficient is larger, and when the subharmonic coefficients is greater than one, the conditions required by stricking periodic motions are rarely fulfilled.

Stability of infinite impact periodic motions

PERTURBATION OF THE DYNAMICS

The motion of the impact oscillators defined by Eq. ( 1) and Eq. ( 3) is highly nonlinear, and so is their dynamics. Therefore, systems may have more than a single response to a given periodic excitation. This phenomenon is called coexistence [START_REF] Guckenheimer | Non-lineur O.scillation.s, Dyamical Systems, und B(furcutions of Vector Fielris[END_REF]. For a given set of external parameters, infinite impact periodic motions have been proved exist together with other finite impact motions (Fig. 5) (T: w, f, Q) = (0.75, 10: 1.1, 1). Depending on the value of the control parameters the existence, as well as the stability of the response, may change. This section is devoted to an evaluation of the stability of infinite impact periodic responses. evaluation can be carried out by using the first order perturbation of impacting dynamics. The relation between the end of sticking phase j and the end of sticking phase j + 1, i.e. between Wl*,j and Wl.,j+r is defined by a Poincare application P.

Local stability can be calculated by using a perturbation method. The first order perturbation term dIVl,j+l , which is the image of first order perturbation term &VI,;, is calculated by differentiation of the dynamic, i.e. through differentiation of P. This first order perturbation is equivalent to calculating the Jacobian matrix DP of the Poincare application. This Jacobian matrix can be determined with the help of the first order terms described in subsection 3. When the impact oscillator is one dimensional, the Jacobian DP cannot be defined. In that case ti is uniquely defined, and it is easy to prove that these motions are always stable.

EUROPEAN JOURNAL OF MECHANICS, A/SOLIDS, VOL. 17, No 2, 1998

LOSS OF LOCAL STABILITY, LOCAL BIFURCATION

Once the Jacobian is calculated, the discussion regarding local stability depends on the eigenvalues of matrix DP. This (2n -2) square matrix is real but not symmetric, so some of its 2r1 -2 eigenvalues may be complex and conjugate. When all of them are inscribed in the unit circle of the complex plane the motion is stable. Conversely if only one of the eigenvalues lies outside this circle, it is unstable. If at least one eigenvalue lies on the unit circle, the dynamics of the system may bifurcate (Gumovski, 1980).

In the latter case, depending on the phase of eigenvalue A, three kinds of bifurcation may be encountered. When X = -1 the dynamics crosses over a flip or subharmonic bifurcation. In that case the previous motion still exists but becomes unstable and the dynamics of the system stabilises on a response whose period is twice the previous one [START_REF] Guckenheimer | Non-lineur O.scillation.s, Dyamical Systems, und B(furcutions of Vector Fielris[END_REF]. When X + 1 the periodic solution disappears into a saddle node bifurcation [START_REF] Guckenheimer | Non-lineur O.scillation.s, Dyamical Systems, und B(furcutions of Vector Fielris[END_REF]. When X = trig 6 8 and cp # Z!Z$ (27r) (Gumovski, 1980), the notion bifurcates into a Hopf bifurcation and the periodic motion turns into a quasiperiodic one.

In the case of a one dimensional impact oscillator, sticking periodic motions are always stable, and it can be proved that such local bifurcations cannot be encountered in their dynamics.

"RISING-BIFURCATION"

Infinite impact periodic responses may lose their stability through a new kind of dynamic phenomenon in which the sticking phase can be interrupted. The sticking phase of the notion ends when the contact reaction is zero and when it changes sign. This reaction is given in Eq. ( 9). The contact reaction is a function of time but it is also a function of the restricted displacement field Y. Therefore, in the case of a multidimensional impact oscillator, i.e. when Y can be defined, the end of the sticking phase is not uniquely defined as in the one dimensional case.

When varying a control parameter, the characteristics of a sticking response can change suddenly if a new zero contact reaction appears between tb and fi. We refer to this phenomenon as a "rising bifurcation" (Fig. 6) (r, w, f, .7:0) = (0.75, 148, 28.5, X2 lo-'). This bifurcation is qualitatively similar to a 'grazing bifurcation' [START_REF] Nordmark | Non-Periodic Motions Caused by Grazing Incidence in an Impact Oscillator[END_REF] which can be encountered in finite impact periodic dynamics.

When a rising bifurcation is encountered, the original sticking periodic dynamics does not exist any more. This kind of bifurcation may lead to other kinds of dynamics, e.g. other kinds of sticking periodic dynamics (Fig. 7), or sticking chaotic dynamics. For example periodic motions with two or three sticking phases per period may result from these bifurcations. Similar to grazing bifurcations, rising bifurcations would merit specific investigation, but will be no longer analysed in the present paper.

Methods of investigation

ANALYTIC PREDICTOR-CORRECTOR METHOD

The method an results detailed in the two sections provide information regarding the existence and stability of sticking responses of multi-d.o.f. impact oscillators, but these results correspond to a given value (T, w, f, XJ:O) of the control parameters. The purpose of constructing diagrams is to extend these results to a wider area inside the parameter space. Many numerical studies have been carried out, including bifurcation diagrams (Thompson, 1982). These calculations were limited to single d.o.f. impact oscillators. At the same time, analytic studies [START_REF] Shaw S | The Transition to Chaos in a Simple Mechanical System[END_REF] were made, but they were restricted to single impact periodic motions of one dimensional oscillators. The authors (Gontier, 1997), (Toulemonde, 1997) have carried out an analytic study of finite impact periodic motions of these kinds of oscillators.

It is useful to remark that the displacement contraint 20 and magnitude f are redundant parameters, where /L = (T, W, f) denotes the reduced parameter vector of the system.

The analytic method used to build bifurcation diagrams is based on the continuation method and results from the differentiation of the Poincare application defined by Eq. ( 25) (Allgower, 1990). The aim of the present investigation is to explore the phase parameter space from two points of view: by determining bifurcation diagrams; and constructing a partitionoing of the parameter space. The first order differentiation of the periodic motions in Eq. ( 26) is the main tool for this investigation, and the procedure defined is a continuation method.

Essentially, there are two different types of continuation methods (Allgower, 1990). The Predictor Corrector (P.C.) methods and the Piecewise Linear (P.L.) methods. The idea of P.C. methods is to trace implicitly defined curves numerically by generating a sequence of points along the curve while satisfying a chosen tolerance criterion. In P.L. methods, the curve is linearlised, piece by piece, then fitted to the system equations. Essentially, the continuation method can be used as a convex homotopic method applied to the solution of the nonlinear Eq. ( 27). Equation ( 27) can be rewritten showing the relation between the parameter vector 11, defined above and the unknown vector A/SOLIDS, Vol.. 17, No 2, 1998

When vector /I, follows a curve in the parameter space, Eq. ( 41) implicitly defines a curve Wr (lb) in the space of unknowns. We assume the existence of a given starting point (/rn. Wr,n) satisfying Eq. ( 41). The aim of this method is to generate a sequence of points (p,]: Wr.,r),,=E~ satisfying Eq. ( 41) with a chosen tolerance criterion.

Assume that (I] : WI.]) is a root of the multidimensional function 4. A new point (/j,,,+i . Wi.,l+l) along the implicitly defined curve is obtained by first making a Predictor

Step, generally an Euler Predictor. This step is based on the first order differentiation of Eq. ( 41):

(42)

In Eq. ( 42) IL represents the "step-size". The manner in which h is to be chosen depends on the application of the P.C. method; for example, it may be controlled by an arclength method. However, solution (l~>+r j IV; j+l) is only a first order prediction and it must be corrected (bj+r , Wr,j+r ) is then obtained by applying a Corrector step to (/~i+r > W,'.j+l). An obvious way to obtain this correction is to use a Newton-like method. Defining the correction by (/l,j+r, Wr,j+r) = (/L\+~: Wi,,j+l) + (dp,, &Vi), the first Corrector step is such that The solution of Eq. ( 43) is not unique, so that another relation between dp and dW1 must be looked for. This relation defines a correction strategy, which depends on the application of the continuation method. In view of constructing a bifurcation diagram, choosing dp = 0 was found to be sufficient. However, when partitioning parameter space, a strategy at 1~1 or 1~2 constant can be deficient because of the shape of the curves, especially if they have asymptotic branches. So another equation is added to Eq. ( 43) which defines the corrector step. This equation determines d,u and then dW1 is derived from Eq. (43).

ADVANTAGES OF THE P.C. METHOD

Many authors have obtained bifurcation diagrams by time integration of the dynamics. Actually, numerical approaches usually produce rather noisy diagrams, whose inaccuracies are inherent in the weak stability of the dynamics and the integration process. There are two reasons why more accurate results are obtained here. First, in some cases, as in the present situation, the dynamics of infinite impact periodic motions is totally described in an analytic manner. This represents a dramatic improvement in accuracy over a step-by-step numerical integration.

The second reason derives from the considerable loss of accuracy of numerical methods when approaching bifurcation or when motions are not stable. In these zones, periodic solutions are considerably less attractive (Thompson, 1982) and numerical methods EUROPEAN JOURNAL OF MECHANICS, A/SOLIDS. VOL.. 17, No 2. 1998 hardly converge. Conversely, both the determination of infinite impact periodic motions described in section 3 and the continuation method described in this section are not sensitive to local instability, as the points in the bifurcation diagram remain defined by an analytic method.

Another advantage of the continuation method is that it requires less CPU time than direct integration methods. The computation of a neighbouring point requires only the calculation of the two differential elements of Eq. ( 41). A very good approximation of the neighbouring point is obtained with only one Predictor step and less than two Corrector steps, which take very little time to perform on a standard 486 PC.

Moreover the continuation method is characteristically efficient and accurate, especially when applied to analytic functions. For a given set of control parameters the method investigates the finite impact periodic responses and it also able to compute the Floquet multipliers analytically. These results are very useful in bifurcation diagrams, as they are the key to an accurate prediction of local bufircations and subsequently of new branches.

BIFURCATION DIAGRAMS

The method described in the last two subsections makes it possible to slide along a branch of a bifurcation diagram. Thus the diagrams presented here will stop at the point at which the motions bifurcate in a non-periodic way-quasiperiodicity occurs after Hopf bifurcation and chaos may appear after rising bifurcation. In the bifurcation diagrams presented in this article the sticking proportion is projected onto the vertical axis and the magnitude f of the harmonic excitation is used as a control parameter.

PARTITIONING THE PARAMETER SPACE

The P.C. method can also be applied to the partitioning of a parameter space, for instance the (w, f) plane. In the parameter plane, curves are drawn representing the motions with a constant sticking proportion. The partitioning lines correspond to constant t'ltb sticking proportion motions, defined as a ratio of ---.

Thus starting from a given point (wa, fu) in the plane associated with a given pro;ortion. a curve can be drawn by a continuation process at constant proportion, representing a set of equivalent sticking orbits.

The algorithm used to calculate these curves is described in subsection 5.1 where a new constant proportion condition is added to Eq. ( 42) during the Predictor step. In the same fashion the Corrector Step must take the error in the calculation of the proportion into account. In this subsection the above methods are applied to the one d.o.f. oscillator. As stated above, if a sticking periodic response exists in the dynamics of a one d.o.f. system, it can EUROPEAN JOURNAL OF MECHANICS.

A/SOLIDS. VOL. 17, No 2, 1998 be easily proved that this motion is locally stable. Local bifurcation such as subharmonic bifurcation or Hopf bifurcation cannot be encountered in this kind of dynamics.

As a consequence, the cascades of flip bifurcations observed in the dynamics of finite impact motions and referenced in the literature [START_REF] Shaw S | The Transition to Chaos in a Simple Mechanical System[END_REF]Gontier, 1997) do not exist in the sticking dynamics of single degree of freedom impact oscillators.

Depending on the sign of displacement constraint :I:(), two classes of impact oscillators and therefore two classes of dynamic behaviour can be encountered. Sticking periodic motions have thus been proved to exist in both cases, depending on whether the systems have clearance or are preloaded.

(i) In the case of a preloaded system, i.e. when z:() > 0, sticking periodic motions exist for every value of the the excitation pulsation. They exist for excitation magnitude values larger than the value of the perload force induced by .r:o; below this force, sticking motions cannot be encountered. These responses exist either with a plastic impact or elastic restitution rule (Fig. 5b).

By applying the method described in section 5, bifurcation diagrams can be constructed. For a given value of excitation pulsation, the excitation magnitude f can be regarded as a control parameter of the dynamic system. It can be proved that Sticking periodic motions can be proved to be non existent for values of .f smaller than the preload induced by JO. But they begin to appear when f is slightly larger. These first responses give rise to a time ratio of sticking dynamics slightly less than 1. It is then observed that the sticking proportion decreases as a function of amplitude (Fig. 8). For values of .f larger than a given level, sticking motions disappear and finite impact responses occur (cf: Eq. ( 38)). This result shows that for a given value of excitation pulsation, sticking period one responses exist only for a finite segment of excitation magnitude. Infinite impact responses of impact oscillators with preload can also be represented in a two dimensional (w, f) parameter space, where sets of sticking motions with the same sticking ratio can be easily represented (Fig. 9). It must be noted that these curves are very smooth and that they are not interrupted by any kind of bifurcation.

EUROPEAN JOURNAL OF MECHANICS.

A/SOLIDS. VOL. 17, N" 2. 1998 f 00 For large values of excitation pulsation, the constant sticking ratio curves are parallel to the pulsation axis. This means the qualitative behaviour of these motions does not depend on the pulsation when it is sufficiently large. For small values of excitation pulsation and for large values of excitation magnitude, the constant sticking ratio curves are parallel to the magnitude axis. In case the qualitative aspect of the impact response does not depend on the excitation magnitude. Partitioning (Fig. 9) clearly shows that sticking motions exist in a restricted area of the parameter space and that they will not be observed in the upper right quarter of the (w, f) plane. It should be noted that the limit provided by these numerical calculations have been predicted by the purely analytic investigation of section 3.8.

The conclusions of section 3.8 prove that subharmonic sticking periodic responses exist and that the responses predicted by these analytic calculations have only one free-flight phase per period. When the restitution coefficient is null, some other kinds of orbits exist, some of which have two sticking phases per period (Fig. 10). EUROPEAN JOURNAL OF MECHANICS. A/SOLIDS. VOL. 17, No 2. I998 (ii) In the case of a mechanical system with clearance, i.e. when .r:u < 0, sticking periodic responses have been proved to exist either with a purely plastic impact rule or with an elasto-plastic impact rule (Fig. 1 Bifurcation diagrams have also been constructed, and for a given value of the excitation pulsation, the field of excitation magnitudes have been explored. The least value of the amplitude for which sticking motions exist here corresponds to a very low sticking proportion. Although local bifurcation cannot be encountered in these dynamics, grazing bifurcations have been found. Here, the single-periodic sticking responses can disappear, and lead to two-periodic sticking responses (Fig. 12). This diagram shows that such motions can be observed for large values of excitation magnitude. A parameter space partitioning has also been constructed for such systems (Fig. 13), where the curves are much more complex than in the preload case, and different classes EUROPEAN JOURNAL OF MECHANICS. A/SOLIDS. VOL. 17, N" 2. 1998 of curves can be observed, which seem to be asymptotically parallel to the magnitude axis. It is easy to prove here that the reaction cannot be maintained at a positive value for half the period, and consequently that the sticking ratio is always smaller than 50%.

6.2. SYSTEMS WITH SEVERAL D.O.F.

The methods described above have been extended to some structures with several d.o.f. There are no major differences in shape between multidimensional responses and motions in the one d.o.f. case. Sticking periodic motions either exist with a purely plastic or with an elasto-plastic impact rule (Fig. 23).

In the sticking dynamics of multi-d.o.f. systems, periodic motions can lose their stability in different ways. For instance, local bifurcation can be encountered, as in the case of a three d.o.f. system a Hopf bifurcation is encountered and a sticking periodic motion is turned into a sticking quasiperiodic response (Fig. 14). Such dynamic behaviour could not have been observed in the single d.o.f. sticking dynamics. Subharmonic bifurcations have not been encountered in these sticking dynamics. This is a major difference with finite impact periodic dynamics where they have been very EUROPEAN JOURNAL OF MECHANICS. A/SOLIDS. VOL. 17, No 2. 1998 25 frequently observed (Gontier, 1997;[START_REF] Shaw S | The Transition to Chaos in a Simple Mechanical System[END_REF] either in one or two dimensional systems (Toulemonde, 1997).

The new kind of bifurcation, called "rising bifurcation" and predicted in subsection 4.3, which cannot be predicted with the help of eigenvalues, can be found in a system with at least two d.o.f. Although it has no major influence on the shape of the motions (Fig. 6) it has a very strong influence on the shape of bifurcation diagrams whose smoothness can be broken by this kind of bifurcation. The sticking proportion decreases very sharply following rising bifurcations (Fig. 15) and some non-periodic motions may result from them. This result can be compared to the existence of chaotic motions after grazing bifurcations (Nordmark,199 1 ), as immediately afterwards a succession of chaotic (Fig. 16) and periodic windows can be observed. The phenomena cited above are typical of the sticking behaviours of multiple degree of freedom impact oscillators, which cannot be encountered in the sticking dynamics of EUROPEAN JOURNAL OF MECHANICS.

A/SOLIDS, VOL. 17, NC 2, 1998 single degree of freedom impacting systems. And it has been proved in section 4.3 that rising bifurcations (Fig. 15) cannot be observed in the dynamics of a one dimensional system. Also, the chaotic behaviour (Fig. 16) can be observed only in multidimensional cases because it has been proved that the sticking motions of the single degree of freedom impact oscillators are always periodic and stable.

A parameter space partitioning has also been constructed for a multiple d.o.f. case, where constant sticking proportion curves are interrupted by the existence of rising bifurcations (Fig. 17). These bifurcations can turn periodic responses into chaotic ones so that the periodic branch cannot be continued. Thus it can be noted that the lower the sticking ratio, the more often the curve is interrupted.

Conclusions

The phenomenon of sticking periodic motions of impact oscillators has been analysed and a general methodology of approach suitable for multi-d.o.f. oscillators as well as one d.o.f. oscillators has been defined; bifurcation diagrams and parameter space partitionings have been determined. In the case of multi-d.o.f. oscillators some specific features such as quasiperiodic sticking motions and a new kind of bifurcation called "rising bifurcation" have been identified. All these examples have indicated the efficiency of the method used to investigate these sorts of problems. At the present time, the tool developed here should be used more extensively to provide a better understanding of the sticking solutions of impact oscillators. Rising bifurcations, transitions from unstable solutions to stable solutions, the limit of chaotic responses are some questions among many which would deserve a more thorough investigation. From a general point of view, multiple d.o.f. impact oscillator behaviour has to a large extent been unexplored up to now.

Fig

  Fig. I. -Multidimensional impacting structure. The equilibrium state of the structure is outlined by the dashed curve

  the time derivative and X is the displacement dt field in the structure. R is the diagonal matrix of the squares of the eigenpulsations (w,", . . *, wl', ). zg (t) is the resulting displacement at point B in the direction Bn:. f and w are respectively the magnitude and the pulsion of the harmonic excitation applied EUROPEAN JOURNAL OF MECHANICS.A/SOLIDS, VOL. 17, No 2, 1998 3 at point A; T-4 is the projection of a unit force applied at A onto the eigenbase of the structure.

  MECHANICS, A/SOLIDS, VOL. 17, No 2, 1998 5The reaction of the stop can also be calculated using the restricted displacement field Y: of a one dimensional impact oscillator, the restricted displacement field Y cannot be defined and Eq. (7) does not exist. In this case, the entire structure is stationary during the sticking phases, and the expression of the reaction reduced to the simple relation

Fig. 2 .

 2 Fig. 2. -Plastic impact rule, sticking periodic motion (two d.o.f. impact oscillator with preload).

  Fig. 3. -Elastic impact rule, infinite impact periodic motion (two d.o.f. impact oscillator with preload)

  and tj+l are determined by the implicit relations 2~ (tj) = zu and z~(t,+r ) = zu. They are therefore two implicit functions of the degree of freedom vector Vj. Consequently, vj+r can be defined as an implicit function of vi.

  Ii(t) = U'(t) + Cwi[--c; sin(w:(t -tb)) + d; cos(w:(t -t{I))]X; i=l Me [tb; t'l[ In Eq. (16) U'(t) is the forced displacement field of the blocked structure induced by the harmonic excitation; Ya is the constant response associated with the constant excitation term of Eq. (8); t; and ti are respectively the initial and final time of the sticking phase. Thus the motion defined in Eq. (16) depends only on tb and on the set of 2n -2 variables h di)iE[l,M]* w e now define IV0 as the 2n -1 dimensional vector describing the begining of the sticking phase: (17) wo = At time ti the restricted displacement and velocity fields in the structure can also be expressed with the help of the eigenbase of the constrained structure (18) L II-1 Y(t:) = Yo + U'(t',) + C e;K i=l II -1 I+:) = iqt'l) + cLdp;r; i=l As shown by relation (9) the reaction of the stop depends only on the time and on the restricted displacement field Y. Time ti and restricted degree of fredom vector Y (ti) are then linked by the relation R(ti) = 0. By taking Eq. (18) into account, it can be expressed EUROPEAN JOURNAL OF MECHANICS, A/SOLIDS, VOL. 17, No 2, 1998

  of the sticking phase can then be reduced to a discrete relation between Wo and WI: (w;(t'l -tb)) sin(wi(t\ -ti,)) -sin(w:(ti -t[))) cos(~i(ti -tb)) 0 0 . . . i

  /-2 --f $jp-2 W.j ++ W,~+I = (G 0 HL 0 FX' 0 HI)(WJ} In Eq. (

  (27) qqW1) = P(Wl) -Wl = 0 EUROPEAN JOURNAL OF MECHANICS. A/SOLIDS. VOL. 17, No 2, 1998 11 This equation is easily solved by using an algorithm based on the Newton method requiring local determination of the gradient D,,,$ of 4.

  the application defined by Eq. (29). A good approximation of the complete sequence of impacts Fm can be obtained with the help of the previous expression(30) Lo Fx'-l -F" k ++x

  Fig. 4. -Subharmonic periodic response with sticking phases; the subharmonic coefficient is equal to 2. ('1.. w. ,f) = (0. 10. 6).

  Fig. 5. -a) Coexistence: single impact period one motion, single d.o.f. impact oscillator, (1.. w. f) = (0.75 , 10. 1.1). b) Coexistence: infinite impact period one motion. single d.o.f. impact oscillator, (7.. w. f) = (0.75, 10. 1.1).

  t-z ---) $&&t-2 dW1.j H dWl,j+l = (DG o DHa o DL 0 DFX' 0 DHl)(dWl,;)

Fig

  Fig. 6. -a) Infinite impact periodic reaction before a rising bifurcation. b) Infinite impact reaction after a rising bifurcation. c) Detail of a sticking motion during a rising bifurcation.

  . SYSTEMS WITH ONLY ONE D.O.F.

  Fig. 8. ~ Bifurcation diagram, one dimensional system with prebad. (1.. ti) = (O.i5. 2.X)

Fig

  Fig 9. -Parameter space partioning, one dimensional system with preioad, I' = 0.75. .A, \. "Y.Y. L,

  Fig. 10. -Periodic response with two free-flight phases per period. purely plastic impact rule. (r. w. ,f) = (0. 10. 4.8).

  Fig. I I -Sticking periodic motion, plastic impact rule. one dimensional system with clearance, (v. w. j') = (0. 1.4. 10)

  Fig. 12. -Bifurcation diagram. one dimensional system with clearance, (I,, u/) = (0, 1.,4)

  Fig. 13. -Parameter space partioning, one dimensional system with clearance, I' = 0.

  Fig. 14. -Quasiperiodic sticking response, three d.o.f. impact oscillator.

  Fig. IS. -Bifurcation diagram, two dimensional system with preload.A rising bifurcation is encountered in this diagram.

  Fig. 17. -Parameter space partitioning, two dimensional system. The curves are interrupted by rising bifurcation.

Acknowledgements.

-The authors would like to thank Michel Panet and Alain Leger at E.D.F., for their assistance.

EUROPEAN JOURNAL OF MECHANICS, A/SOLIDS. Vol.. 17, No 2. 1998


C S Chang

P.-Y Hicher

Model for Granular Materials with Surface Energy Forces

Keywords: Granular material, Stress-strain relations, Anisotropy, Elastoplasticity, Moon

In light of environmental differences ͑ such as gravitational fields, surface temperatures, atmospheric pressures, etc.͒, the mechanical behavior of the subsurface soil on the Moon is expected to be different from that on the Earth. Before any construction on the Moon can be envisaged, a proper understanding of soil properties and its mechanical behavior in these different environmental conditions is essential. This paper investigates the possible effect of surface-energy forces on the shear strength of lunar soil. All materials, with or without a net surface charge, exhibit surface-energy forces, which act at a very short range. Although, these forces are negligible for usual sand or silty sand on Earth, they may be important for surface activated particles under extremely low lunar atmospheric pressure. This paper describes a constitutive modeling method for granular material considering particle level interactions. Comparisons of numerical simulations and experimental results on Hostun sand show that the model can accurately reproduce the overall mechanical behavior of soils under terrestrial conditions. The model is then extended to include surface-energy forces between particles in order to describe the possible behavior of lunar soil under extremely low atmospheric pressure conditions. Under these conditions, the model shows that soil has an increase of shear strength due to the effect of surface-energy forces. The magnitude of increased shear strength is in reasonable agreement with the observations of lunar soil made on the Moon's surface.

Introduction

The soil on the Moon formed by space weathering processes and dynamic impacts of micrometeorites is composed of aggregates of minerals, rocklets, and glasses, welded together into agglutinates. The glass, which is abundant in lunar soil, is present in the form of sharp, abrasive, interlocking, fragile shards and fragments. The lunar subsurface soil can be classified as dry silty sand ͑Perko et al. 2001͒. However, observations and measurements conducted during Surveyor ͑1966-1968͒, Apollo ͑1968-1972͒, and Luna ͑1959-1976͒ missions indicate that lunar soil has an unusually high cohesion in comparison to the same type of soils tested under the Earth's conditions. According to Buzz Aldrin, astronaut on Apollo 11, the soil on the Moon ''has the cohesive property that wet sand would have [START_REF] Costes | Apollo 11: Soil mechanics results[END_REF]. A stable vertical trench excavated during the Surveyor 7 mission is shown in Fig. 1͑a͒ and a footprint left by Surveyor 3 demonstrates the fineness and cohesiveness of the lunar soil ͓Fig. 1͑b͔͒.

The fine particles of lunar soil are the products of the continual impact on the surface by meteoroids, which smash and grind rocks into soil and weld soil into new rocks. The writers we are interested in the behavior of soils in a typical regolith of the Moon-a surface layer that consists of loose sand and rock fragments, which overlies the solid rock. The range of grain size distribution is summarized in a soil gradation curve given in Fig.

2 [START_REF] Costes | Apollo 11: Soil mechanics results[END_REF][START_REF] Mitchell | Soil mechanical properties at the Apollo 14 site[END_REF][START_REF] Scott | Lunar soil mechanics[END_REF] based on the core samples obtained from the surface layer of the Moon ͑Apollo 11 and 12͒. According to the particle size distribution, the lunar soil can be classified as silty sand following the Unified Soil Classification System, ASTM standards.

The writers' knowledge on engineering properties of lunar soil are derived mainly from operation performed during the Surveyor and Appolo programs, including bulk and core samples returned to the Earth, trench excavation, ground penetration, and some simple physical experiments. In laboratory studies, it was found that soils tested under a chamber with ultrahigh vacuum and high temperature had an increase up to 13°in angle of internal friction and 1.1 kPa in cohesion ͑Bromwell 1966; [START_REF] Nelson | Environmental effects on engineering properties of simulated lunar soils[END_REF][START_REF] Perko | Surface cleanliness effect on lunar soil shear strength[END_REF]. This increase of soil strength is consistent with the observed stable vertical trend excavation in dry silty sand on the Moon.

For soil on the Earth, cohesion is originated from electrostatic forces and surface-energy forces, which includes van der Waals forces ͑intermolecular potential energy͒. The surface-energy forces act at very short range. Comparing to the gravitational forces, the surface-energy forces are negligible on the Earth for soils with particle sizes larger than around 0.06 mm. Thus, terrestrially, the forces are negligible for sand or silty sand. However, as particle size decreases to the size of clay, the surface-energy forces have a greater influence on its strength. The effect of surface-energy forces on friction and adhesion properties of small particles can be found in the literature of powder technology, particulate technology, and interface sciences ͑Derjaguin 1934; [START_REF] Derjaguin | Effect of contact deformations on adhesion of particles[END_REF][START_REF] Johnson | Surface energy and contact of elastic solids[END_REF][START_REF] Molerus | Effect of interparticle cohesive forces on flow Behavior of Powders[END_REF]Thorton and Ning 1998;[START_REF] Tomas | Assessment of mechanical properties of cohesive particulate solids. Part 1: Particle contact constitutive model[END_REF] Due to various environmental differences on the Moon ͑such as the gravitational fields, temperature extremes, low to practically nonexistent atmospheric pressure, plus the space weatheringm, which includes surface activating effects like UV radiation, the solar wind and galactic particle fluxes͒, it seems plausible to consider surface-energy forces as a dominant cause for altering the mechanical behavior of lunar soil with the size of silty sand. An issue that may be addressed is to what order of magnitudes the surface-energy forces can influence the strength and deformation of lunar soils.

In this paper, a model is developed that accounts for the surface-energy forces between particles and is capable of modeling its effect on the shear strength of the soil assembly. For this purpose, a microstructural modeling approach is adopted ͑Chang and Hicher 2005; Hicher and Chang 2005͒. Models using this approach can also be found in the work by Jenkins ͑1988͒, Walton ͑1987͒, Rothenburg and Selvadurai ͑1981͒, Chang ͑1988͒, Emeriault and Cambou ͑1996͒, Liao et al. ͑2000͒, Kruyt and Rothernburg ͑2002͒, among others. Elastoplastic models using this approach can be found in Chang et al. ͑1992a,b͒, Suiker and[START_REF] Suiker | Modelling failure and deformation of an assembly of spheres with frictional contacts[END_REF].

It is noted that the above-mentioned microstructural continuum approach ͑Chang and Hicher 2005͒ is not a "complete" micromechanical approach. As the complete microstructural details for a particle assembly are not known, the complex deformation behavior of the particle assembly cannot be derived directly from the first-principles theory based solely on particle-level properties. Two major simplified assumptions have been applied. 1. Besides the parameters of particle-level properties, the present model adopts some macroscale parameters, which are related to the degree of interlocking and porosity for the soil assembly. These macroscale parameters are adopted because interparticle behavior is not exclusively a local phenomenon; it is considerably influenced by the degree of interlocking and porosity for a soil assembly. In the present model, these parameters are based on the critical void ratio concept, and are empirically determined from testing on soil specimens. 2. Particle rotation is neglected in the model derivation. Particle rotation is a significant factor that influences the stress-strain behavior of assembly. A model considering particle rotation would result in the use of higher-order stresses and strain such as micropolar strain and couple/Cosserat stress ͑Chang and Liao 1990; [START_REF] Chang | Elastic material constants for isotropic granular solids with particle rotation[END_REF][START_REF] Chang | Micromechanical modeling of constructive relations for granular material[END_REF]Gao 1996, Luding et al. 2001;[START_REF] Chang | On virtual work and stress in granular media[END_REF][START_REF] Tejchman | Modeling of a cyclic plane strain compression-extension test in granular bodies within a polar hypoplasticity[END_REF]. Currently, there is no stress-strain model with explicit consideration of particle rotation that can successfully simulate complex behavior of soil. The difficulties of establishing such a continuum model mainly arise from the lack of experimental results on measured particle rotation field and related couple stress for an assembly of soil particles. With the previous discussion on restrictions and simplified assumptions, the present model should be regarded as a model of semiempirical nature. It is necessary to calibrate parameters from soil tests in order to predict complex behavior with good quantitative agreement. The main advantage of the present model over conventional continuum constitutive models is its microscale consideration, which allows one to conveniently extend the model for new phenomenon at the particle level, such as the surface-energy forces.

However, adding the surface-energy term from first-principle physics into an otherwise semiempirical model may elicit the following concern: how can one part of the physics be treated from first principles if the other parts of the physics are treated only empirically? The answer is that even though the surface-energy forces are added from first principle physics, it does not alter the empirical nature of the model. The model is aimed to predict the "order of magnitude" of the surface-energy effect.

In what follows, we first describe the formulation of this model, which takes account of the surface-energy forces between particles. The model performance is then demonstrated through the results of a simulation of fine Hostun sand tested under terrestrial conditions ͑i.e., without surface-energy forces͒. The fine sand has a mean particle diameter of 0.5 mm, loaded to failure in triaxial testing cells under different confining stresses. Finally, the effect of surface-energy forces between particles is introduced to predict the stress-strain-strength behavior of lunar soil.

Stress-Strain Model

In this model, we envision a granular material as a collection of particles. The deformation of a representative volume of the material is generated by the mobilization of contact particles in all contacts. Thus, the stress-strain relationship can be derived as an average of the mobilization behavior of local contact planes in all contacts. For a contact plane in the ␣th orientation, the local forces f i ␣ and the local movements ␦ i ␣ can be denoted as follows:

f i ␣ = ͕f n ␣ , f s ␣ , f t ␣ ͖ and ␦ i ␣ = ͕␦ n ␣ , ␦ s ␣ , ␦ t ␣ ͖
, where the subscripts n, s, and t represent the components in the three directions of the local coordinate system. The direction normal to the plane is denoted as n; the other two orthogonal directions, s and t, are tangential to the plane. Rotation of particles is not considered here.

The forces and movements at contact planes of all contacts are suitably superimposed to obtain the macroscopic stress strain tensors. The macroscopic stiffness tensor is obtained on the condition that the rate of energy dissipation expressed in terms of the macrostress and strain must be equivalent to that expressed in terms of microforces and movements. In such formulation, it has usually been assumed that the microstructure is statically constrained, which means that the forces on each contact plane are assumed equal to the resolve components of the macroscopic stress tensor. Another equally simple possibility is to assume a kinematically constrained microstructure, in which the movements, rather than the forces on a contact plane, are the resolved components of the macroscopic strain tensor.

The kinematically constrained models are more popularly used in concrete models. The overriding reason for using this approach is that, in case of strain softening, it is easier to construct a stable model using a kinematic rather than a static constraint. However, the kinematic constraint gives a more restrictive deformation pattern, thus providing a stiffer result, especially when granular materials are subjected to high deviatoric stresses ͑Chang and Misra 1990͒. Further discussion on this issue can be found in the work by Chang and Gao ͑1996͒, Kruyt and Rothenburg ͑2002͒, and Kruyt ͑2003͒.

Interparticle Behavior van der Waals Forces

The surface-energy forces pull soil particles together, thus increase the shear strength of soil. The physics of the surfaceenergy force between two particles is reviewed in this section. The adhesive forces between two solid grains mainly result from electrostatic force and van der Waals forces. The electrostatic component in lunar soil is assumed to be negligible ͑Perko et al. 2001͒. Therefore, in this paper, the surface-energy forces are calculated from van der Waals energy fields. There are several situations that can alter the surface properties and surface energy, such as monolayers of water or passivation layers. For simplicity, lunar soil grains are represented by spheres of equal radii that are separated by thin layers of adsorbed air molecules.

van der Waals forces between two bodies are derived from the dispersion interaction energy between two identical atoms or molecules ͑Israelachvili 1992͒:

W͑r͒ = - C r 6 ͑1͒
where r = distance between the two atoms and C = London dispersion coefficient.

Then, with the assumption of additivity of these attractive forces, the interaction energy of a molecule located at a distance D of the planar surface of a solid made up of like molecules is the sum of its interaction with all the molecules in the body:

W͑D͒ = - C 6D 3 ͑2͒
Where = number of atoms per unit volume. The corresponding van der Waals force is

F = ␦W͑D͒ ␦D = C 2D 4 ͑3͒
It is then possible to calculate the interaction energy between two solids as, e.g., two spheres: and1 and2 = number of atoms per unit volume for the two bodies. For two flat surfaces, one obtains the interaction energy per unit area:

W = - AR 1 R 2 6D͑R 1 + R 2 ͒ ͑4͒ D = distance between two spheres; R 1 and R 2 = radii of the two spheres; A = 2 C 1 2 = Hamaker constant;
W = - A 12D 2 ͑5͒
The van der Waals force between two solids in contact can then be computed from the interaction energy. One obtains the following two cases; Between two spheres:

f = 2W͑D͒ R 1 R 2 R 1 + R 2 ͑6͒
Eq. ͑6͒ is called the Derjaguin approximation ͑Israelachvili 1992͒.

In the case of two identical spheres, Eq. ͑6͒ becomes

f = AR 12D 2 ͑7͒
Between two flat surfaces, per unit area:

f = A 6D 3 ͑8͒
Let us now consider two soil grains as two identical elastic spheres. If subjected to an external force f, the two spheres will deform and create a flat circular contact area with a radius a

a = ͩ 3͑1 -p 2 ͒R 4E p ͪ 1/3 f 1/3 ͑9͒
where E p and p = Young's modulus and Poisson's ratio of the particles, respectively. The derivation of Eq. ͑8͒ can be found in Under these simplified conditions, van der Waals forces acting on the two particles can be considered as the sum of two terms, one due to the interaction between two flat surfaces of area S = a 2 and one along the remaining surface of the two spheres. Using the Derjaguin approximation for the second term, we obtain the expression of the van der Waals force between two particles:

f = A 6D 3 a 2 + AR 12D 2 ͑10͒
where A = Hamaker coefficient; D = thickness of molecules layer between two particles; and R = radius of particles. Hamaker constant A was estimated to be 4.3ϫ 10 -20 J for lunar soil and 1.5 ϫ 10 -20 J for terrestrial quartz sand ͑Perko et al. 2001͒. The thickness of molecules layer between two particles D are highly dependent on the atmospheric pressure and composition. On the Moon, the atmospheric pressure is nearly zero, which can lead to a very thin layer of molecules between two particles compared to that under a terrestrial environment. Therefore, according to Eq. ͑10͒, it is reasonable to expect that the surface-energy forces between particles are much higher than those between particles under a terrestrial environment.

It is noted that if the radius a of the contact area in Eq. ͑10͒ increases with confining stress of a specimen. Thus, the surfaceenergy force also increseases with confining stress, which indicates that the surface-energy force will contribute to the shear strength not only on the cohesive component but also on the frictional component.

The orientation of a contact plane between two particles is defined by the vector perpendicular to this plane. On each contact plane, an auxiliary local coordinate can be established as shown in Fig. 3. The surface-energy force f given in Eq. ͑10͒ represents only the force magnitude. Its direction is normal to the interparticle contact plane given by n i ␣ . Thus in a vector form, the surface-energy force f i ␣͑SE͒ = fn i ␣ . Between two particles, the interparticle force f i ␣ can be decomposed into two components: ͑1͒ due to applied load on the boundary of the soil assembly f i ␣͑A͒ , and ͑2͒ due to surface-energy forces between particles f i ␣͑SE͒ . Thus,

f i ␣ = f i ␣͑A͒ + f i ␣͑SE͒ ͑11͒
The surface-energy force is a function of contact area as shown in Eq. ͑10͒. However, the contact area is, in turn, a function of interparticle force f i ␣ ͓see Eq. ͑9͔͒, thus the expression of Eq. ͑11͒ is implicit and nonlinear in nature.

Interparticle Force-Displacement Relationship

Elastic Part

The contact stiffness of a contact plane includes normal stiffness, k n ␣ , and shear stiffness,

k r ␣ ͑assuming k r ␣ = k t ␣ = k s ␣ ͒.
The elastic stiffness tensor is defined by

f i ␣ = k ij ␣e ␦ j ␣e ͑12͒
which can be related to the contact normal and shear stiffness

k ij ␣e = k n ␣ n i ␣ n j ␣ + k r ␣ ͑s i ␣ s j ␣ + t i ␣ t j ␣ ͒ ͑ 13͒
where n, s, and t = three orthogonal unit vectors that form the local coordinate system ͑see Fig. 3͒. The vector n is outward normal to the contact plane. Vectors s and t are on the contact plane. The value of the stiffness for two elastic spheres can be estimated from Hertz-Mindlin's formulation ͑Mindlin 1969͒. For sand grains, a revised form was adopted ͑Chang et al., 1989͒, given by

k n = k n0 ͩ f n G g l 2 ͪ n , k r = k r0 ͩ f n G g l 2 ͪ n

͑14͒

where G g = elastic modulus for the grains; f n = contact force in normal direction; l = branch length between the two particles; k n0 , k r0 , and n = material constants. For two spherical particles, the branch length is the same as the particle size l = d. If the Hertz-Mindlin's contact formulation is used, the value of k n0 in Eq. ͑14͒ can be expressed as

k n0 = G g d 2 ͩ ͱ 12 1 -g ͪ 2/3

͑15͒

Plastic Part

The elastic sliding behavior between two particles does not have a coupling effect ͑i.e., the sliding direction is on the contact plane; no movement component is in the normal direction.͒. However, the plastic sliding often occurs with an upward or downward movement, thus the shear dilation/contraction takes place. The stress dilatancy is a well-recognized phenomenon in sand ͓see discussions in the work by Taylor ͑1948͒, Rowe ͑1962͒, and Goddard and Bashir ͑1990͔͒, and should be correctly modeled. The dilatancy effect has been described by ͑Chang and Hicher 2005͒

␦ ˙n p ⌬ ˙p = f r f n -tan 0 ͑16͒
where 0 = material constant which, in most cases, can be considered equal to the interparticle friction angle . Note that the shear force T and the rate of plastic sliding ⌬ ˙p are defined as

f r = ͱ f s 2 + f t 2 , ⌬ ˙p = ͱ ͑␦ ˙s p ͒ 2 + ͑␦ ˙t p ͒ 2

͑17͒

The yield function is assumed to be of Mohr-Coulomb type,

F͑f i ,͒ = f r -f n ͑⌬ p ͒ = 0 ͑18͒
where ͑⌬ P ͒ = isotropic hardening/softening parameter. When F Ͼ 0, it indicates loading, otherwise unloading. The hardening function, defined as a hyperbolic relationship between and ⌬ p , involves two material constants: p and k p0 ͑Chang and Hicher 2005͒.

= k p0 tan p ⌬ p ͉f n ͉tan p + k p0 ⌬ p ͑19͒
The value of asymptotically approaches tan p . The initial slope of the hyperbolic curve is k p0 . On the contact plane, under a yield condition, the direction of plastic shear sliding ⌬ ˙p follows the associated flow rule, thus, it is perpendicular to the yield surface. However, the plastic movement in the direction normal to

Interlocking Influence

One of the important elements to be adopted in granular modeling is the critical state concept. Under critical state, the granular material remains at a constant volume whereas it is subjected to a continuous distortion. The void ratio corresponding to this state is denoted as e c .

The critical void ratio e c is a function of the mean stress. The relationship has traditionally been written as follows: The interparticle friction angle is a constant for the material. However, the peak friction angle, p , on a contact plane between two particles is dependent on the degree of interlocking of neighboring particles, which can be related to the state of packing void ratio e by tan p = ͩ e c e ͪ m tan ͑21͒

where m = material constant ͑Biarez and Hicher 1994͒.

For dense packing, the peak frictional angle p is greater than . When the packing structure dilates, the degree of interlocking and the peak frictional angle are reduced, which results in a strain-softening phenomenon.

Elastoplastic Force-Displacement Relationship

With the elements discussed previously, the relations between the rate of force displacement for two particles can be derived that includes both elastic and plastic behavior, given by

f ˙i ␣ = k ij ␣p ␦ ˙j ␣ ͑22͒
The detailed expression of k ij ␣p can be found in Chang and Hicher ͑2005͒.

Stress-Strain Relationship

Macro-and Microrelationships

The stress-strain relationship for an assembly can be determined from integrating the behavior of all interparticle contacts. During the integration process, a relationship is required to link the macro-and microvariables. Using the static hypotheses proposed by Liao et al. ͑2000͒, we obtain the relation between the macrostrain rate and interparticle displacement rate ͑here, we do not consider the finite-strain condition͒

u ˙j,i = A ik -1 ͚ ␣=1 N ␦ ˙j ␣ l k ␣ ͑23͒
where ␦ ˙j = relative displacement rate between two contact particles and the branch vector l k = vector joining the centers of two contact particles. It is noted that contact particles include both direct contact and indirect contact of neighboring particles associate with a Voronoi polyhedron as discussed by Cambou et al. ͑2000͒. For convenience, we let N be the total number of contacts.

The variables ␦ ˙j ␣ and l k ␣ are defined, respectively, as the values of ␦ ˙j and l k associated with the ␣th contact. The fabric tensor A ik in Eq. ͑23͒ is defined as

A ik = ͚ ␣=1 N l i ␣ l k ␣ ͑24͒
Using the principle of energy balance and using Eq. ͑23͒, the mean force rate on the ␣th contact is

f ˙j ␣ = ˙ij A ik -1 l k ␣ V ͑25͒
In Eq. ͑25͒, the stress increment ˙ij can be obtained by the contact forces and branch vectors for all contacts ͑Christofferson et al. 1981; Rothenburg and Selvadurai 1981͒:

˙ij = 1 V ͚ ␣=1 N f ˙j ␣ l i ␣ ͑26͒
Applying the defined contact force in Eq. ͑25͒, Eq. ͑26͒ is unconditionally satisfied. Because of its approximation nature, Eq. ͑25͒ can be viewed as an averaged solution, in which the interparticle force can be regarded as the mean value for forces on all contact planes of the same orientation. For convenience, we also regard the branch length as the mean value for all contact planes of the same orientation.

The interparticle force in of Eq. ͑26͒ consists of two parts: f i ␣͑A͒ due to applied load on the boundary of the soil assembly, and f i ␣͑SE͒ due to surface-energy forces between particles ͓see Eq. ͑11͔͒. Thus the stress can be separated into two parts:

˙ij = 1 V ͚ ␣=1 N f ˙j ␣͑A͒ l i ␣ + 1 V ͚ ␣=1 N f ˙j ␣͑SE͒ l i ␣ ͑27͒
The second part of Eq. ͑27͒ represents the stress induced by surface-energy forces denoted as

͑ ˙ij ͒ SE = 1 V ͚ ␣=1 N f ˙j ␣͑SE͒ l i ␣ ͑28͒
It is noted that this term is not analogous to the usual concept of cohesion for bulk material. In the present model, the surfaceenergy stress depends on the packing structure and is a tensor rather than a scalar. Only for an isotropic distribution of the branch lengths l ␣ , the surface-energy stress can be reduced to an isotropic tensor.

As mentioned earlier, the surface-energy force is a function of contact area, which is in turn a function of interparticle force. Thus all equations in the second section are nonlinear in nature. The set of nonlinear equations for the stress-strain relationship is discussed in next section.

Computation Scheme

The problem is defined by the following. Initially, we know the global variables ͑ ij and ij ͒ for the assembly and the local variables ͑f j ␣ and ␦ j ␣ ͒ for each contact orientation. For a given loading increment, which can be stress control, strain control, or mixed mode, out of the 12 variables ͑ ˙ij and ˙ij ͒, 6 of them are unknown. The objective is to determine all global variables ͑ ij and ij ͒ and local variables ͑f j ␣ and ␦ j ␣ ͒ at the end of load increment. For a system with N interparticle orientations, the number of unknown is 3N for f ˙j ␣ and 3N for ␦ ˙j ␣ . The total number of unknown is 3N +3N +6.

The following constraints must be satisfied: 1. The local constitutive equation, i.e., Eq. ͑22͒. As there are three equations for each contact plane orientation, the total number of equations is 3N, N being the total number of interparticle orientations. 2. Static hypothesis between global stress and local forces, i.e., Eq. ͑25͒: the number of equations is 3N. 3. Strain definition between global strain and local displacement, i.e., Eq. ͑23͒. The number of equations is 6 ͑strain is symmetric͒. The total number of unknowns is the same as the total number of equations. Therefore, the solution can be determined.

Using Eqs. ͑22͒, ͑23͒, and ͑25͒, the following relationship between stress and strain can be obtained:

u ˙i,j = C ijmp ˙mp ͑29͒
where

C ijmp = A ik -1 A mn -1 V ͚ ␣=1 N ͑k jp ep ͒ -1 l k ␣ l n ␣
The summation in Eq. ͑29͒ can be replaced by an integral over orientations. The integral can lead to a closed-form solution for the elastic modulus of randomly packed equal-size particles ͑Chang et al. 1995͒. However, in the elastic plastic behavior, due to the nonlinearity nature of the local constitutive equations, a numerical calculation with iterative process is necessary to carry out the summation in Eq. ͑29͒. In order to facilitate the numerical calculation, the orientations are selected to coincide with the locations of Gauss integration points in a spherical coordinate. Summation over these orientations with the Gauss weighting factor for each orientation is equivalent to determining the integral over orientations. Results were found to be more accurate by using a set of fully symmetric integration points. From a study of the performance of using different numbers of orientations, we found N ജ 74 to be adequate ͑Chang and Hicher 2005͒.

For a strain-controlled test, Eq. ͑29͒ is not useful especially at the after-peak range of strain softening. In this case, a method of "elastic predictor-plastic corrector" was adopted to obtain the solution. For a mix-mode loading condition, an additional process of distributing the unbalance stresses was also needed ͑Chang and Hicher 2005͒.

Summary of Parameters

One can summarize the material parameters as: 1. Normalized contact number per unit volume: Nl 3 / V.

Mean particle size, d:

• Interparticle elastic constants: k n0 , k r0 , and n;

• Interparticle friction angle: and m;

• Interparticle hardening rule: k p0 and 0 ; and • Critical state for the packing: and ⌫ or e ref and p ref .

Besides critical state parameters, all other parameters are for interparticles. Standard values for k p0 and 0 are the following: k p0 = k n0 and 0 = and a typical ratio k r0 / k n0 = 0.4 can be generally assumed. Therefore, only six parameters have to be derived from experimental results and they can all be determined from the stress-strain curves obtained from drained compression triaxial tests.

Regarding the effect of the particle size, Eq. ͑15͒ shows that the interparticle stiffness k n0 is influenced by the mean particle size. As in this model, the tangential stiffness k r0 and the plastic stiffness k p0 are directly related to k n0 , they are also particle size dependent. Considering an assembly of particles under an externally applied stress, according to Eqs. ͑14͒ and ͑15͒ and using the expression of contact force in terms of external stress, the interparticle stiffness is proportional to the mean particle size ͑see Liao et al. 2000͒:

k n ϰ G g 1-n dͩ Nl 3 V ͪ -n n ͑30͒
Although the magnitude of the interparticle stiffness is dependent on the particle size, the elastic shear modulus of a packing is not particle-size dependent. This can be derived by the fact that the shear modulus of a packing is proportional to

k n / d, thus it leads to ͑Liao et al. 2000͒ G ϰ G g 1-n ͑͒ n ͩ Nl 3 V ͪ 1-n

͑31͒

Note that the dimensionless density factor Nl 3 / V is a function of the void ratio only. Thus, the shear modulus of the packing is not particle-size dependent.

Results of Numerical Simulation

Triaxial Tests on Terrestrial Soil without Surface Energy Forces

A series of drained triaxial tests on fine Hostun Sand were performed by Al Mahmoud ͑1999͒. The gradation curve tested sand is shown in Fig. 4. The mean size of the particles for fine Hostun sand is d = 0.5 mm. It is classified as uniform fine sand. The minimum and maximum void ratio for this sand is e min = 0.575 and e max = 0.943, respectively.

It is assumed here that surface-energy forces can be neglected under terrestrial conditions. This assumption will be discussed in the next section. The model needs a small number of input parameters, such as mean particle size, particle stiffness, interparticle friction, initial porosity tensor, and an initial degree of interlocking. The interparticle elastic constant k n0 is assumed as 61, 000 N / mm. Considering the grain size distribution curve, the contact numbers per volume Nl 3 / V are 2 for dense sand and 0.9 for loose sand ͑Hicher and Chang 2005͒. The evolution of Nl 3 / V due to new contact generation with neighboring particles was not considered.

The value of k r0 / k n0 is commonly about 0.4, corresponding to a Poisson's ratio for Hostun Sand = 0.2 and the exponent n = 0.5 ͑Biarez and Hicher 1994͒. From test results, we were able to = 34°. The equation governing the dilatancy rate requires the determination of the parameter 0 . This parameter represents the concept of "phase transformation" as defined by Ishihara and Towhata ͑1983͒ or "characteristic state" as defined by Luong ͑1980͒. A value of 0 = was retained.

The peak friction angle is not an intrinsic parameter, but varies with the void ratio according to Eq. ͑21͒. A value of m = 1.5 was determined from the test results. The values of k p0 are directly connected to the elastic properties through the relation:

k p0 = k n .
The set of parameters for fine Hostun sand is presented in Table 1.

The tests were performed at different confining pressures on samples prepared at different initial void ratios. Typical results are presented in Figs. 5 and6, which show the triaxial testing results for both dense and loose specimens made of Hostun sand. The stress-strain curves are plotted for dense sand in Fig. 5, and loose sand in Fig. 6. The model performance can be demonstrated by comparing the predicted and measured macrobehavior.

One can see the combined influence of the initial void ratio and of the mean effective stress on the stress-strain curves and the volumetric change. The stress-strain curve has a peak, which increases with confining stress. For dilatant materials, the deviatoric stress reduces after peak and converges toward a constant state of stress corresponding to the critical state. In practice, it is difficult to reach this state because of strain localization, especially in dense materials. The critical state was estimated from the results on loose specimens.

Triaxial Tests with Consideration of Surface Energy Forces

In order to evaluate the possible effect of lunar environmental conditions on the stress-strain response of soils, the influence of the surface-energy forces for different values of interparticle distance D was examined. The stress-strain curves were predicted using the same soil parameters given in Table 1. Given the introduction of van der Waals forces, shear strength is expected to be higher. Fig. 7 shows the influence of the distance D between the particles on the shear strength for soil specimens under 20 kPa confining stress, where q 0 represents the shear strength for soil without surface-energy forces and q represents the shear strength with the effect of these forces.

It is noted that when the distance D is greater than 2 nm, the effect of the van der Waals forces is practically negligible. This distance D is a function of the amount of molecules, which can be adsorbed on the solid surface. The adsorbate thickness may be estimated by means of potential theory ͑Adamson 1990; [START_REF] Perko | Surface cleanliness effect on lunar soil shear strength[END_REF]. Molecule adsorption on a solid surface is conditioned by temperature, gas pressure, and atmospheric composition. On Earth, adsorption conditions are easily fulfilled due to the presence of a high atmospheric pressure. Therefore, the distance D can easily exceed several nanometers such that the effect of the van der Waals force is negligible. Under lunar atmospheric condition, the thickness of the adsorbed molecule layer is likely to be very thin. For the case of 0.3 nm thickness ͑D = 0.6 nm͒, the comparisons of the stress-strain curves in the Earth's and the Moon is atmosphere environments are given in Fig. 8 for dense soil and Fig. 9 for loose soil under a 20 kPa confining stress. Assuming the thickness D is 0.3-0.6 nm, the computed shear strength for lunar soils is 12-15% higher than that of the Earth's soil. The stiffness of the lunar soil is also higher. By plotting the maximum shear strength envelope in the Mohr plane for specimens under different confining stresses, the effect of van der Waals forces for the dense sand would result in an increase of 2 -3 kPa in cohesion, whereas an increase of 0.5 kPa in cohesion is sufficient to hold a 1 m trench cut in lunar gravity condition.

Summary and Conclusion

In the present study we introduce surface-energy forces between particles and estimate its effect on the shear strength of soil. The microstructural approach is adopted because it allows us to model interactions at particle level.

In the model, a simple elastic-plastic behavior was assumed on each contact plane. The elastic part is based on the Hertz-Mindlin's contact formulation, whereas the plastic part is based on a Mohr-Coulomb friction law with an isotropic hardening assumption and a nonassociated flow rule. The interaction among interparticle planes is assumed to be dependent on the degree of interlocking and porosity of the assembly of soil. The effect of porosity is modeled by a phenomenological approach using the concepts of critical state such that strain softening behavior can be modeled for dense materials. On the whole, the model requires a limited number of parameters, which can easily be determined from conventional triaxial tests. Ability of the model to reproduce the main features of sand behavior has been demonstrated. Model simulations were compared with drained triaxial tests results at different initial void ratios and different confining stresses leading to contractant or dilatant behavior of the sand specimens. The comparison demonstrated that the model is capable of reproducing the general trend for both loose and dense sands.

The model was applied to predict the shear strength of lunar soil under extremely low atmospheric pressure. Model predictions indicate that soil under extremely low atmospheric pressure has an increase of shear strength by 12-15% higher than the one which would be present under Earth's atmospheric pressure. The magnitude is in the same order as the measured increase of shear strength for lunar soil simulants tested under the usual atmospheric pressure and under a chamber with ultrahigh vacuum. However, the predicted shear strength for lunar soil due to surface-energy forces is only a hypothesis, which must await further testing on the Moon for validation/calibration.

Fig. 2 .

 2 Fig. 2. Approximated gradation range for Apollo 11 and 12 soil core samples

Fig. 1 .

 1 Fig. 1. ͑a͒ Trench excavated by Surveyor 7 Spacecraft in Lunar subsurface soil ͑Perko et al. 2001, ASCE͒; ͑b͒ the footprint of Spacecraft ͑Surveyor 3-Apollo 12͒ ͑NASA Image ID number AS12-H-48-7110͒

Fig.

  Fig. Local coordinate at interparticle contact

  e c = ⌫ -log͑pЈ͒ or e c = e ref -logͩ pЈ p ref ͪ ͑20͒ where ⌫ and = two material constants; pЈ = mean stress of the packing, and ͑e ref , p ref ͒ = reference point on the critical state line.

Fig. 4 .

 4 Fig. 4. Gradation curve for fine Hostun sand

Fig. 6 .

 6 Fig. 5. Stress-strain curves for dense Hostun sand

Fig. 8 .Fig. 9 .

 89 Fig. 7. Effect of distance between particles on shear strength

  Derjaguin et al. ͑1975͒, Johnson ͑1971͒, Dahneke ͑1972͒, and Valverde et al. ͑2001͒. It is noted that for simplicity, we consider only the case of elastic contact flattening ͓Eq. ͑9͔͒. For the case of elastic-plastic flattening, refer to Molerus ͑1978͒, Thornton and Ning ͑1998͒, and Tomas ͑2001͒.

Table 1 .

 1 Model Parameters for Fine Hostun Sand

		p ref
	e ref	͑MPa͒		͑deg͒	0 ͑deg͒	m	k p0 / k n
	1.05	0.01	0.14	34	34	1.5	0.02

Dept. of Civil and Environmental Engineering, Univ. of Massachusetts, Amherst, MA 01002. E-mail: chang@ecs.umass.edu

2 Research Institute in Civil and Mechanical Engineering, UMR CNRS 6183, Ecole Centrale Nantes-Univ. of Nantes, France.