M Jean

Frictional contact in collections of rigid or deformable bodies: numerical simulation of geomaterial motions

The purpose of this paper is to present some general numerical methods for treating dynamical problems involving unilateral contact and dry friction. Some examples of applications related to the structural response of rigid or deformable geomaterials such as, rocks, soils, collections of blocks, granular materials, are given. Emphasis is put on Coulomb's dry friction law. This law is relevant for a large class of applications to geomaterials. It accounts for the main features of dry friction. It may be easily improved without drastic changes in the proposed methods. The frictional problems appear to be strongly non linear, and call for the techniques of nonsmooth mechanics. Convex Analysis is widely used to formulate friction equations and numerical algorithms.

INTRODUCTION

When modelling a mechanical problem, some mechanical variables are selected and the equations governing the system are formulated. Some equations derive from mechanical principles such as the equations of motion or the equations of continuum mechanics which are universally accepted. Other equations such as constitutive laws for a material, or frictional contact laws, are often complex and difficult to define and may be unreliable. Finally some numerical method is chosen to fmd approximate solutions to the system of equations. When the main interest is devoted to the description of physical phenomena, numerical results are of assistance. There is some inclinacy to take numerical results for granted, especially when graphical outputs are produced showing some similarity with observed phenomena. The numerical results should not be blindly accepted when sensitive mechanical phenomena are studied, and the influence of the modelling process and in particular the choice of frictional 1 contact laws should be examined. The influence of computational parameters, such as the time step and convergence criteria are also important Many complicated phenomena are involved when contact and friction occur between rigid or deformable bodies. Reproducible behaviours may be experimentally difficult to obtain. Moreover in some cases reproducibility should not be expected due to either insufficient control of experimental circumstances or due to lack of information. Frictional contact laws are usually written as relations between the local stress at some point of the interface and the relative velocity. Formulating such laws and selecting particular expressions depend on the time and length scales of the investigation, the scope of applications, the expected results, and the methods used to solve the equations of the problem. For instance, time and length scales decide whether the materials should be considered as either rigid or deformable.

Generally, the equations of motion govern the evolution, but if inertia effects are negligible with respect to external and internal forces, one may end with a quasi-static problem. When a collection of rigid bodies is under consideration, for instance a granular flow, or a wall made of blocks under seismic excitation, these problems have to be treated in a dynamic sense. In applications such that strain-stress experiments in granular material, or progressive loading of a wall, interest is focussed on the quasi-static behaviour, while supposedly negligible dynamical effects certainly occur. In such cases it is difficult to propose consistent equations governing the system, other than the equations of motion.

In this paper, finite dimensional systems which may be collections of rigid bodies or models of continuous media will be considered. Finite dimensional models of continuous media are obtained through such methods as finite elements methods. The question of the choice of the fmite elements is not within the scope of this paper, which is devoted to numerical treatment of contact and friction.

Frictional contact laws are presented in part l. The equations governing the problems are written in part 2, and the discretization processes are presented in part 3. The method used to compute the frictional contact forces is presented in part 4. Some examples are given in part 5.

FRICTIONAL CONTACT LAWS

Complicated phenomena, occuring at the so-called microscopic scale, take place between contacting bodies or edges of discontinuities. Surfaces in contact are found to be rough, and contact is established through asperities which suffer damage during compression and sliding, so fragments of material are generated. Molecular interactions may occur as well. When investigating the structural response of the system, such processes are omitted. A frictional contact law should appear as a relation between tangential and normal stresses and the sliding velocity at the point of contact. This law may be constructed from experimental data. It may also emerge from theoretical analysis based on multiple length-scale such as the homogeneization technique. In such techniques a neighborhood of the contacting zone is considered, including parts of the two bodies, and sometimes also a third body interacting with the two bodies. When the thickness of this neighborhood vanishes the homogeneization technique yields a constitutive law at the interface, or a contact element in finite element methods, (Licht,[START_REF] Licht | Comportement asymptotique d'une bande dissipative mince de faible rigidite[END_REF]). This law or this contact element accounts for interfaces features and also for the response of the materials.

When large sliding displacements or deformations are investigated, other parameters such as the sliding displacement have to be introduced to take into account wear or structural phenomena such as dilatancy in granular materials. When such a degree of complexity has to be attained, and when micro-scale variables are purposely ignored for the developpement of the model, a phenomenological model might prove valuable, (Selvadurai, this volum, Boulon, [START_REF] Boulon | Basic features of soil-structure interface behaviour[END_REF]).

When developping a frictional contact law, it is first generally assumed that motions are smooth, i.e. the velocities are continuous functions of time and forces may be described as piecewise continuous functions of time. Nevertheless when frictional contact occurs, the velocities may have jumps and reactions appear as impulses. The proper mathematical tool to describe jumps of the velocity is the concept of a function U with bounded variation on the considered time interval. This secures the existence of the left limit u-, i.e. the value before the impact, and the right limit u+, the value after the impact. Derivatives of such functions are measures, such as the Lebesgue measure and the Dirac measure at each point of discontinuities. The densities with respect to the Lebesgue measure describe the usual forces, and the densities with respect to the Dirac measure describe the impulses. Lebesgue forces are usually applied. Impulses are exerted when rigid bodies collide. In some applications, it seems relevant to consider the bodies to be rigid. When deformable bodies are subjected to impact, they are kept into contact during a certain interval of time, the stiffer the bodies, the shorter the interval. When performing a double time-scale analysis, the stiffness coefficients of the bodies are introduced as increasing parameters, so that the bodies tend to become rigid, the duration of the contact might vanish, and in most cases the contact locus reduces to isolated points. It may happen that kinematic variables, such as the relative velocities before the impact, and after the impact, and the time density of impulse at the contact point, are found to satisfy a relation, which is referred to as a shock law. For instance some elementary models, like contacting particles or rigid spheres, interacting through a frictionless thin elastic layer, yield the well known restitution shock law when the stiffness tends to infinity: the normal relative velocity 'UN+ after the impact and the normal relative velocity 'UN-before the impact, satisfy the requirement 'UN +=-e'UN -. where eE [0,1] is Newton's classical coefficient of restitution. Furthermore, if the contact episode obeys Coulomb's law, the same restitution law still emerges while the sliding velocity 'Ur + after the impact and the time density of impulse 9{., are found to satisfy Coulomb's law, (Jean,[START_REF] Jean | Dynamics of rigid bodies with dry friction and partially elastic collisions[END_REF]). If a single contact is involved this frictional shock law proves to be dissipative. If other shocks occur the dissipativity cannot be proven, except for the inelastic case e=O. Moreau [START_REF] Moreau | Unilateral Contact and Dry Friction in Finite Freedom Dynamics[END_REF] has introduced the dissipation coefficient o: the linear combination (1-o)/2'lr+ (l+o)/2'ZJ+, and the time density of impulse 9{. satisfy a relation similar to Coulomb's law. The dissipativity of such a law is proven. Generally, there is no evidence, neither experimental, nor theoretical, arising from a double time-scale analysis for instance, that a shock law may be exhibited or even more be relevant A number of authors have mentioned cases where such laws fail to be dissipative (Stronge, (5)). These remarks apply to collection of rocks or stone blocks.

In this paper, only the inelastic shock law will be used as outlined below. In paragraphs 1.1, ... ,1.5, motions are assumed to be smooth. The extension to motions with discontinuities will be discussed at the paragraph 1.6.

Local variables

For the sake of simplicity we consider the 2-dimensional case of a body 0 , candidate for contact with some opposing body 0', at some timet These bodies may be rigid or deformable ones. Some special particles P of CJO are defined as particles candidates for contact. For instance if 0 is described by a mesh, the nodes of the boundary may be defined as candidates for contact Another criterion which may be used in general situations is that the particle P be a proximal point to 0'. An opposing particle P' to a candidate P is defined as a particle of the boundary oO' of 0' which is an onhogonal projection of P. A local frame is defined by its origin P', by the unit normal vector N directed from P' toward P, and by a tangent vector T to CJO' at P', (figure 1). Such a frame might be easily constructed when the boundaries of 0 and 0 ' are smooth. However, a candidate for contact P might have several projection P'. Bodies 0 and 0' should be partitioned so that a candidate P for contact be associated with a unique pair of objects, namely a candidate for contact object, and an opposing object, allowing one to construct a unique opposing particle P' and a unique local frame. So it may happen, that a particle P be several times candidate for contact, defmed by a single candidate object, but various opposing objects. The following local variables are defined: N coordinate of P in the local frame, or gap, components of the reaction from 0' exerted on P, components of the relative velocity of P with respect to 0', Vp-Vp• ; when contact occurs, 'Uris the sliding velocity.

When two continuous deformable bodies 0 , 0', are contacting, particles of some continuous subset of ao are possibly candidates for contact The reaction from 0' exerted on each particle of ao candidate for contact, is described by a line (2-dimensional case), surface (3-dimensional case) density of force, i.e. a stress vector. Since numerical methods are the subject of this paper, only a finite number of degrees of freedom is considered. When dealing with such finite models, a fmite number of candidates for contact is selected, for instance nodes or Gauss points in a finite mesh, and reaction forces are exerted on these candidates, possibly affected with some weight coefficients. To derive the relations governing these forces from those frictional contact laws adopted for the contacting continuous media is a mathematical and numerical analysis problem, (Jean,[START_REF] Jean | Unilateral conctact with dry friction: time and space discrete variables formulation[END_REF]). Here, for the sake of simplicity, it will be assumed that frictional contact laws are applied to nodes of meshes, which proves to be a correct approximation for linear elements, and is considered to be correct when appropriate weight coefficients are used.

Unilateral contact

The main features of unilateral contact are:

Impenetrability: 'IN ~ 0 .
ii Bodies are assumed not to attract themselves when contacting: 'IN = 0 :::

) 'l{.N ~ 0 .
iii lfP is is not in contact with 0, the reaction force is still mentioned, but with a null value:

'IN > 0 :::) '/{.N = 0 .

The relations i, ii, iii, may be summarized all together in one of the three equivalent relations,

'IN~ 0 'l{.N ~ 0 and 'IN ~ 0 and
Another equivalent form is

'IN '/{.N = 0 ' 'IN <SN -'l{.N ) ~ 0 • (JN-qN)'l{.N ~ O •
there exists p > 0 such that 'l{.N = proj R+ ('l{.N -P 'IN) ,

(if this relation is true for some p > 0, it is actually true for any p > 0).

(

) (2) 1
(4)

The symbol projE x denotes the orthogonal projection of x on a convex set E. The relation (1) is known as a complementarity condition, or Signorini's condition.

(x) = 0 if xe E, \!IE( X) = -toO if X~ E;
the variables ~ and -~N are conjugate with respect to the pair of conjugate functions

\!'*R+ = \!'R-, \!'R+ •
In this paper the mathematical details are not considered. This property is only mentioned to emphasize the fact that some useful mathematical properties, such as monotonicity, underlie Signorini's condition. For brevity, Signorini's condition will be referred to as follows:

S (~, ~N) is true. Figure 2.
Signorini's condition graph.

Friction law

Coulomb's law is first presented since it accounts for the main features of dry friction. In

the 2-dimensional case, ~T e [-Jl~N, ~N] , Jl friction coefficient, ' l.Lr > O => ~T = -Jl~N ' ' l.Lr < O => ~T = Jl~N ( 7
)
'UT is the sliding velocity. The graph of this relation is displayed on figure 3. Another equivalent form is the so-called principle of maximal dissipation,

\:f .5-re C (8)
where C denotes the interval [-Jl~N, Jl~N] . This relation is in tum equivalent to the following:

there exists p > 0 such that ~T = proj C (~T-P ' l.Lr) •

(if this relation is true for some p > 0, then it is true for any p > 0).

Here again the above relations are equivalent ways of expressing the convex analysis property:

the variables -'Ur and ~T are conjugate with respect to the pair of functions 'l'*c, 'l'c.

The function 'l'*c conjugate of'l'c actually equals 'I'* c (-'Ur) = I'll Tl .

'l'*c may be viewed as a dissipation "pseudo-potential". The same holds in the 3-dimensional case. The convex set C then equals the disk, C = { ~ : 11~11:::;; g } , g = Jl ~N . For brevity, Coulomb's law will be referred to as follows:

Regular forms of frictional contact laws

The graphs of figures 2 and 3 are not the graphs of mappings, since lfN is neither a function of ~N nor ~N a function of lfN. Similarly, 'Uris neither a function of ~T nor ~T a function of 'llr Convex analysis allows one to deal with such graphs. Usual techniques of regular nonlinear analysis may be applied only to graphs of mappings. A classical example of such graphs of the latter sort is displayed on figures 4 and 5. The graph on figure 4, shows the normal reaction force opposing interpenetration as a linear function of the negative gap. The slope of the graph, namely a stiffness coefficient, is supposed to be large enough to restrict penetration at an acceptable level. Such a mechanical behaviour appears realistic if one figures out that the boundaries of the contacting bodies are coated with a thin elastic layer, or if the possible asperities are elastic. On the graph on figure 5, when the sliding velocity ~ is vanishing, the friction force is proportional and opposite to this velocity. This is viscous damping, with a viscosity coefficient large enough to ensure a reasonably small sliding velocity. For numerical purposes, the sliding velocity is approximated as the ratio, .1.y.1.t, where .1-t is the time step and ~is the increment of tangential displacement. Thus, the friction force appears as proportional to the displacement from a reference position, the end position of the previous time step. This is intepreted as the action from elastic layers or elastic asperities. Nevertheless, especially when contact involves high pressures and large sliding distances, more complicated phenomena should be expected, like plasticity and wear. This suggests that the graphs on figures 4 and 5 should be smoothed. Besides, using smooth graphs allows one to apply smooth nonlinear analysis, (Oden, Martins, [START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF]). Nevertheless, reproducible experiments which could produce reliable values of tangential or normal elastic stiffness coefficients, tangential viscosity coefficients, or any physical value related to frictional contact, are still unavailable. With the elastic shear behaviour, it is more advisable to introduce one more variable, the shear elastic displacement ~LT, together with Coulomb's law, to avoid any error when interpretating the graph on figure 5. Suppose a tangential loading causes an elastic shear ~~ and a relative velocity 'liT with the same sign. When the loading is reversed, the relative velocity changes its sign while the elastic shear is still ~~. The graph on figure 5 does not make any difference between ~ and ' ZLr , and shows a friction force with the incorrect sign.

Complicated phenomena such as wear, solid lubrication, existence of a joint lying between the contacting bodies, may produce a frictional behaviour rather discrepant from Coulomb's law. For instance the graph of the friction law displayed on figure 6, distinguishes a static friction coefficient and a dynamical one. In this example, the friction force 9{_T may be considered as the sum of two terms: 9{_T = ~ + ~ ; ~ is a friction force obeying Coulomb's law while ~ = fi 'Ur) is a smooth function of the sliding velocity ' ZLr.

Relaxed contact and friction Jaws, thick graphs

As it has been mentioned in paragraphs 1.2 and 1.3, Signorini's condition and Coulomb's law have interesting properties in the context of Convex Analysis. Many techniques such as quasi-variational inequalities, differential inclusion, piecewise continuous mapping fixed point theories, may be used to deal with such laws. This is a reason to favour them when constructing numerical algorithms. It has been noticed that these laws are not adequate to describe complex phenomena. A way to overcome the inherent uncertainty of the situation is to define relaxed frictional contact laws. A pair 'lN• ~• is said to satisfy a relaxed Signorini's unilateral condition, up to some given gap margin ~'lN• and some given reaction force margin ~~. if there exists a pair o'lN• o~. with lo'lNI~'lN• lo~~~~. such that the pair 'lN+o'lN, ~+5~. satisfies Signorini's condition. The set of pairs qN, ~• satisfying the relaxed Signorini's condition is displayed on figure 7, as what may be called a thick graph. In the same way one may define a relaxed Coulomb's law, up to some given sliding velocity margin A~ and some given friction force margin A~. The corresponding thick graph is displayed on figure 8. These relaxed frictional contact laws account for the lack of experimental information. They also account for numerical computation errors, since the computed values obtained from an algorithm satisfy the asserted laws up to certain margins only. For the applications of the paragraph 5, Signorini's condition and Coulomb's law up to certain margins are adopted.

Non smooth motions

Only smooth motions have been considered in the foregoing. The variables used to describe a frictional contact law are, the gap ~• the sliding velocity ' ZLr , and the reaction force 1(. The frictional contact law is referred to under the form,

S (~, !l(N) is true, C!l(N ('Ur, !l(T) is true.

Signorini's condition

Coulomb's law When non smooth motions occur, i.e. discontinuous velocities, the relation to be obtained should involve only those elements significant for bounded variation functions, i.e. the right and left limits, 'U+, 'U-, and the densities of impulse. Signorini's condition is positively homogeneous with respect to ~ , i.e. if a pair, ~• ~• verifies S, the pair ~• A.~ , verifies also S, where A. is any positive number. The relation S may thus be readily extended to the case where ~ is a density of impulse. The relation C is also positively homogeneous with respect to 1(, i.e. if a pair, 'Ur• 1(, verifies C, the pair 'Ur• A.!!(, verifies also C. The relation C may thus also be extended to the case where 1( is a density of impulse. But what interpretation should then be given to ' ZLr? When one of the two contacting bodies is deformable, Coulomb's law is extended as When two rigid bodies are contacting, a double time-scale analysis should be performed, to understand or model the short physical process producing the steep changes of velocity, and see if a shock law may finally be exhibited. This is not always the case. For the applications of the present paper, inelastic shocks are assumed namely: and Coulomb's law is extended as

C~N ( 'lLy +, ~T) is true.
It may be proved that Signorini's condition together with the inelastic shock law is a complete law in the sense of Moreau [4]. This law is found equivalent to

~JN>O=> ~=0, 1JN = 0 => 'UN+ ~ 0 ~ ~ 0 'UN+ ~ = 0 . (Ibis)
The proof is partly based on the property that when a contact occurs at some time t, a strictly negative value 'UN+ would yield strictly negative values of the gap ' IN in some interval ]t, t+'t].

EQUATIONS OF THE FRICTIONAL CONTACT PROBLEM

Kinematic equations

A mechanical system with n degrees of freedom is to be described by a variable qe Rn. When the system is a discrete model of a continuous medium obtained through such a method as the fmite element method, the variable q is for instance the vector of node displacements. If t ~ q(t) is a smooth motion, the first and second time derivative are denoted respectively, q, q.

Let pa be a candidate for contact (upper Greek indices are numbering candidates for contact, a=l, ... , x. where xis the number of candidates). The relative velocity ~of the candidate pa with respect to the antagonistic body may be written,

~ = H*(q)a q' (9a)
where H*(q)a is the transpose of a linear mapping H(q) from Rn into R2 (2-dimensional case) or R3 (3-dimensional case). This relation is obtained from mere kinematic considerations. It involves the matrices giving the components of vectors in the local frame from those in the general frame. The dual relation is an expression of the representative Ra of the local reaction force ~.for the system of parameter q, (lOa)

The following property allows one to relate the gap to the nonnal relative velocity, (lla) Considering all candidates for contact, introducing the 2x (or 3x) vectors 'l1 = ('lli, ... ,'Ua, .. ,'UX), '/{_= ('l{_l, ... ,'I<?•, .. ,~). the above relations may be summarized as

'l1 = H*(q) q , R = H(q) '1(, 'UN = qN • (9) (10) (11)
The linear mapping H(q) from R2X into Rn is injective if and only if H*(q) from Rn into R2X is surjective. The mapping H(q) is injective, if there exists no system of local reaction forces '1( with a null representative R; the mapping H*(q) is smjective, if any system of relative velocity 'l1 is the image of some q. In other words H(q) is an injective mapping or H*(q) a surjective mapping if constraints are kinematically independent This is not the general case, for instance when two rigid bodies have several contact points. Another example is when several nodes of the same element are contacting the same boundary element.

The equation of Dynamics

Lagrange's equation governing the smooth motion of the mechanical system may be written as, M(q) q = F(t, q, q ) , [START_REF] Kruyt | Towards Micro-Mechanical Constitutive Relations for Granular Materials[END_REF] where M(q) is the inertia matrix, F represents external forces, and quadratic inertia tenns with respect to q. When the system is a discrete model of a continuous medium obtained through such a method as the finite element method, F stands for external or internal forces as well.

When frictional contact occurs, the representative R of frictional contact forces '1( has to be written at the righthand side of the dynamical equations together with the other forces, M(q) q = F(t, q, q) + R. [START_REF] Sinopoli | Dynamic analysis of a stone column excited by a sine wave ground motion[END_REF] This equation has to be written together with the kinematic relations, ( 9), ( 10), [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], and the frictional contact law. So far smooth motions have been considered. When contact and friction occurs, the velocities are not generally smooth, but must be expected to be bounded variation functions, while the reaction forces'/{_, R, are densities with respect to positive real measures as already discussed in paragraph 1.6. The equation ( 13) must be understood in the sense of distributions, or written as a differential equation in the sense of measures, (dt is the Lebesgue measure on the real line R, d9 is a positive measure, standing for the Lebesgue measure as well as for the Dirac measure), M(q) dq = F(t, q, q )dt + Rd8.

(13bis)

The kinematic relations may be readily extended to the case where q and 'll are bounded variation functions, and 9{. and R are densities of measures. The question to establish whether the frictional contact laws may be extended to the case of bounded variation functions and density of measures is not merely a mathematical matter and has already been discussed in 1.6. The equations of the frictional contact problem are as follows: the equation of Dynamics, M(q) dq = F(t, q, q )dt + Rd8, the initial conditions, q(to)=qo. the kinematic relations, 'll = H*(q) q ' R = H(q) 9{., 'UN = qN ' the frictional contact Jaws, for a=l, ... , x.

S (iJNa, 9{_N<l

) is true, Signorini's condition (or relaxed Signorini's condition) C9{.N ('Urn+, 9{.T) is true. The question of the existence and uniqueness of the solutions of these equations is not within the scope of this paper. The literature provides theorems of existence dealing with the frictionless case or with regularized frictional contact laws, for the dynamical or the quasi-static case. If the friction coefficient is "small enough", there exists a unique solution. With Coulomb's law, even when the constraints are kinematically independent, (H(q) is an injective mapping), the existence and uniqueness of a solution is not ensured, and a locking phenomenon may occur. This occurs a fortiori when the constraints are kinematically dependent, for example when collections of rigid bodies are concerned.

DISCRETIZED EQUATIONS

The equation of Dynamics

To begin with, the non contact case is considered. M(q) dq = F(t, q, q )dt .

One considers a time step, ] tr, tr+J], h = tr+J-tr . The symbol q(I) denotes an approximate value of q(tr) and q(l+ 1) an approximate value of q(tr+ J); the meaning is similar using q(n and q(tr). q(l+ 1) and q(tr+t). Integrating both sides of the equation of Dynamics on the time interval and using the approximations, q(I+ 1)-q(I) JM(q) dq JF(t, q, q )dt = hq(I+l) , = M(q(I+1) )(q(I+1)-q(I)) , = hF(i,q(I+1),q(I+l)), t=tr+I , allows one to derive a discretized form of the equation of Dynamics, M(q(l+ 1))(q(l+ 1)-q(l)) = hF(t,q(l+ 1 )),{I(I+ 1))

(1)

The terms q(l+ 1), q(l+ 1), are obtained as limits of sequences qk, qk, defined by the following algorithm. The matrix M(q) is assumed to be a slowly varying function of q. The term F may be predicted using a first order expansion, as usually done when performing a Newton Raphson method. One sets Mk= M(qk) ,

~ * 0 0
Vk = -aF( t , qk , qk ) I ()q , ""\ • 0 Kk = -aF( t , qk , qk) I ()q ,

• 0 Fk = F( t , qk , qk ) ,

~ = Mk + h Vk + h2Kk , Wk = 1(\'k-1 .
The matrices Mk , Vk , Kk , are supposed symetric positive definite so that fu has a symetric positive definite inverse matrix. One chooses, qk -q(I) = h qk , qk+l -q(I) = h qk+l .

An approximate value ofF is * 0 * 0 0 0 F(t, qk+1, qk+1) = F(t, qk + (qk+1 -qk ), qk + (qk+1-qk)). Fk+1:::Fk-hKk(qk+1-qk)-Vk(qk+1-qk)

The algorithm can be written as qk+1= qk + Wk(-Mk ( qk-q(l)) +hFk) , qk+1-q(l) = h qk+1 .

A more general statement is,

0 0 # 0 0 qk+1= qk + Wk(-Mk ( qk-q(I)) +hFk) , qk+1-q(l) = h qk+1 ' #
where Wk is a syrnetric positive defmite matrix. The usual initial values are, qO = q(l) , qO = q(D .

If the sequences qk, qk, converge to some values q(l+ 1 ), q(I+ 1), they are found to satisfy the equation (1). If Mk = M is a constant matrix, and F(t, q, q) = Kq + Uq + P(t), # where, K, U, are symetric positive definite matrices, when the choice Wk = (M+hU+h2K)-1 is made, the algorithm converges within a single iteration.

Remark 1: A first order time expansion has been used to approximate q, q. Higher order approximations are most often used in dynamical computations. It is certainly very valuable as far as smooth motions are concerned. But when contact and friction occurs, the velocity q is not differentiable anymore in the usual sense. It is a function with discontinuities. At a point of discontinuity, the derivative q in the sense of distributions is a Dirac measure. So, a second order approximation scheme has to manage with approximate values of q, which are quite perturbating large numbers. Approximations such as, q(I+1)-q(l) = h(l-8) q(I) +h8 q(I+l), may also be used giving numerically stable algorithms provided Se]l/2,1]. # Remark 2: When the choice Wk = !Clk-1 is made, the corresponding algorithm is implicit Notice that if the mass matrix Mk is vanishing the algorithm is found to be the one generated by the Newton-Raphson method for quasi-static formulations. When Mk = M is a constant diagonal # matrix, choosing Wk = M-1 , and performing only a single iteration on the index k, one gets a fast explicit scheme. When large elasto-plastic deformations are involved, implicit methods may be favoured, since the equation of Dynamics is then satisfied with external and internal forces computed at the end of the time step, while when using explicit methods these forces are computed at the end of the previous time step, with the risk of error accumulation. The implicit method is most costly since derivatives ofF have to be computed and a large linear system has to be solved for each iteration k, while a straightforward computation is performed when using the explicit scheme. Nevertheless, the time step should be small enough to ensure numerical stability and sufficient accuracy, while the time step is usually larger when using implicit methods which prove to be unconditionally stable. Of course many other choices of correcting # matrices Wk are possible.

The discretized equations of the frictional contact problem

When contact and friction occur, the representative of the reaction forces must be writen together with other forces. Approximate values of the relative velocities and gaps, at the end of the time step, are denoted, 'U(I+l), qN(I+l). The right limit of the relative velocity 'lJ+ is involved in the writing of Coulomb's law and of the inelastic shock law, paragraph 1.6. The choice is made to consider 'U(I+ 1) as an approximation of 'll'". The local reaction impulses and the representatives of these impulses, during the time step, are denoted, h~l+ 1), hR(I+ 1). The approximated value R(l+l) is sought as the limit of a sequence Rk satisfying

0 0 # 0 0 qk+l = qk + Wk(-Mk ( qk-q(l)) +hFk +hRk+l), qk+l-q(l) = h qk+l .
Notice that the unknown is hR(I+l) and not hR(I). Notice also that Rk+l has been written at the second hand rather than Rk. The scheme is thus implicit with respect to the frictional contact variables. Approximate values ~I+ 1), 'U(I+ 1), ~(I+ 1), are sought as the limits of sequences, ~. 'llk, ~k . The kinematic relatrions are written using the discretized form,

'Uk+l = H*(qk+l) qk+I, Rk+l = H(qk+l) ~+1 , ~k+l = ~(I)+ h'UNk+l .
The last relation comes from 2.1 [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. The term ~a(l) is the gap for the contact a at the beginning of the time step. When the radius of curvature of contacting regions is large enough with respect to the distance run within a time step, the local frames do not change much, and an explicit computation of values of H and H* may be done once for all, for instance H(qk), H(qk), or merely H(q(l)), H*(q(l)). Such values are noted Hk, Hk•. The discretized contact and friction relations are

S (~ak+l, 1LNak+l) is true, C1LNk+l ('Urak+l, 1LTak+l) is true, 1LNak+l => 'UNak+l = 0 .
(16) (15)
A number of methods may be used to solve the system of equations of paragraph 3, such as variational inequalities techniques or fixed point techniques, (see a review, Jean, [START_REF] Jean | Simulation numerique des problemes de contact avec frottement, Materiaux et Techniques[END_REF]). The method presented here is the one used to compute the examples of paragraph 5. Using the kinematic relations, (15), ( 16), (17), and the equation of Dynamics [START_REF] Raous | Experimental analysis of the rocking of a rigid block[END_REF] When the inequality -1 < Jla < 1 is not satisfied, the solution is not unique. In the 3- dimensional case, the solution is not straightforward and it is numerically computed using a generalized Newton method which yields the solution within a few iterations, (Jean,[START_REF] Jean | Numerical methods for three dimensional dynamical problems[END_REF]).

4.4. The solution when several contacts are involved The unknowns ua., ~a. I a.=1, ... , x. (actually UO.k+l, ~a.k+l 'but the indices k+l, k' are omitted here) are sought as limits of sequences ua.p, ~a.p (actually subsequences Wk+l,p, ~k+l,p ), (Jean,Moreau,[START_REF] Jean | Unilaterality and dry friction in the Dynamics of rigid bodies collections[END_REF]). Suppose that provisional values of the local reactions ~Pare adopted for the candidates~. The equation (18) writes, ua.p+l = U lib a.+ L wa.~ h~~P + wa.a. h~p+l ~~ ' together with Signorini's condition and Coulomb's law written for the candidate a.. A value of ua.p+l , ~a.p+l, is readily obtained from the formulas of paragraph 4.2. All candidates for contact are successively examined, repeatedly, until some convergence criterion is satisfied. Namely the approximate values are prescribed to satisfy relaxed Signorini's condition and Coulomb's law up to certain given margins. The convergence of this relaxation algorithm has not been mathematically established. Nevertheless, the numerical convergence is satisfactory.

EXAMPLES OF GEOMA TERIAL MOTIONS

Compression of geologic layers

Two viscoelastic layers, 10 kms deep, 20 kms long, are subjected on the left edge to a constant velocity equal to v0=0.33 J0-10 mls toward the right edge, while only vertical displacements are allowed on the right and left edges, and only horizontal displacements occur along the bottom edge. The layers are subjected to gravity. The viscoelatic layers obey a simplified behaviour law, Maxwell's law cr*=2J.1D+A.tr(D)I-2)'1ldev(cr), A.=vEJ(l +V)(l-2v), Jl=El2(1 +V), 11=1/2y, cr is the Cauchy stress tensor and cr* is the Green-Naghdi derivative. The stiffness coefficient E is equal to IQIO Pa, the Poisson coefficient v=0.25. The upper layer has a viscosity coefficient 11 equal to IQ23 Paxs, and the lower layer, quite less viscous, has a viscosity coefficient 11 equal to IQ21 Paxs. The experiment duration equals 106 years. A fault (crack) crosses the two layers. The layers are discretized with 1890 linear triangular elements TI, and 2062 degrees of freedom; 25 nodes are candidate for contact on each edge of the fault. Signorini's condition and Coulomb's law are assumed on the fault, with a friction coefficient equal to 0.2. An explicit scheme with respect to q, q, is used to integrate the equation of Dynamics, and an implicit scheme is used to compute the reaction forces as described in paragraph 4. One is interested in the quasi-static evolution of the layers. Nodes are equipped with fictitious adaptive masses in order to allow one the use of large time steps, (about 3000 steps for the duration of the experience), the ratio of inertia forces to internal forces being controled during the computation. The final state is displayed on figures 10, 11, 12. The figure 10 displays the sliding velocity between the edges of the fault (divided by vO) versus the depth. It shows that no sliding occurs on the lower part of the fault, and the sliding velocity rapidly increases, from the lower layer, up to the surface of the upper layer. Large deviatoric stresses are developed in the neighborhood of the fault when crossing from a layer to the other one, figure 11. Large deviatoric strains are developed in the neighborhood of the fault in the lowe r layer, figure 12. This example has been computed by J. Chery, Laboratoire de G~opbysique, Montpellier.

Bi-directional strain-stress experiment of a Schneebeli material

The Schneebeli material consists of a collection of rigid disks or rolls. It is widely used as an approximate model of soil, either in experiments, or in numerical simulations. A number of authors have used the numerical program TRUBAL, or improved versions of this program, originally developped by Cundall, [II], to perform numerical simulations of classical mechanical tests. In the TRUBAL program, the contact is established through a system of springs and dampers, and the tangential stress follows Coulomb's law. An explicit integration scheme is used for the equation of Dynamics. Numerical strain-stress experiments presented in this paper have been performed with the program LMGC developed by the author, using relaxed Signorini's condition, relaxed Coulomb's law and an implicit Euler integration scheme for the equation of Dynamics, as described in the previous section. The results obtained by Yemmas, [START_REF] Yemmas | Simulation numerique de materiaux granulaires. These Mecanique[END_REF], are compared with those obtained by Cambou, Mahboubi, Ecole Centrale de Lyon, France, using TRUBAL. This comparison was the object of a research program within the GRECO Geomateriaux (a report is to be published). Using TRUBAL, the question of the influence of the regularizing coefficients, namely the tangential and normal stiffness and viscosity coefficients has to be investigated. This question has been widely surveyed by Kruyt, [START_REF] Kruyt | Towards Micro-Mechanical Constitutive Relations for Granular Materials[END_REF]. It seems that for the strain-stress experiment, the results do not depend too much on these coefficients. When using LMGC, the question of the permitted margins is also raised. To provide a comparison, an interpenetration less than l/100 of the minimal disk radius was allowed. The facilities offered by the two programs are different and it is so far impossible to perform exactly the same numerical experiments. For instance, the samples of material are not generated in the same way, boundaries conditions are not applied in the same way, and the stress tensor is not computed in the same way. The sample used for the purpose by LMGC contains 256 disks (rolls with a unit lenght of I m), 48 disks with radius 1.6 mm, 80 disks with radius 1.05 mm, 128 disks with radius 0.65 mm, figure 13. Samples with 1024 disks have also been used. The upper frictionless wall, is submitted to a constant vertical force P=3.300 N, while the right hand side frictionless wall, moves with a constant velocity V=lO crn/s and induces a resulting reaction force R. The two other walls are fixed and frictionless. The ratio of P to the box width L, is referred to as cr2, and the ratio of R to the height box as crl. The initial width box is Lo. The response crl/cr2 versus the strain (L-Lo)!Lo is rather erratic. The disks are moving discontinuously. A system of rigid disks may produce locking, and the deformation is possible only if some ease is allowed between disks, for instance allowing elastic deformation in TRUBAL, or allowing margins of interpenetration in LMGC. The figure 13 displays the internal or global angle of friction cp, defined by sin$= (crl-cr2)/(crl+cr2) versus the local angle of friction <p between disks. The experimental results were obtained by Abriak, Ecole des Mines de Douai, (to be published). These angles are computed using "the maximum value" of crl/cr2. The figure shows that when the local angle <p becomes large, the global angle <1> seems to reach a steady value. One may also note that when <p is large enough, the void ratio increases when the strain (L-Lo)/Lo increases. Disks instead of sliding on each other, rather roll, the sample expands, and more ease is allowed between disks. Consequently <1> does not increase as much as <p.

A temple-like building made of rigid blocks

This example has been computed by Moreau (private communication), with a program us ing an algorithm as described in this paper. Blocks are rigid and inelastic shocks are assumed to occur between comers and edges. A horizontal deviation of the ground (1 meter) is applied for 1 second, figure 13. This is an example of numerical simulation performed in order to understand the behaviour of ancient buildings made of blocks without any joint when subjected to earthquakes. More complicated examples may be found in [START_REF] Jean | Unilaterality and dry friction in the Dynamics of rigid bodies collections[END_REF]. Examples with a single block are also interesting, since they allow one to make comparisons with rigid body models or finite e lement viscoelastic m odels, wi th analytical solutions obtained for special motions. (Sinopoli,[START_REF] Sinopoli | Dynamic analysis of a stone column excited by a sine wave ground motion[END_REF]) or e xperiments, (Raous,[START_REF] Raous | Experimental analysis of the rocking of a rigid block[END_REF]). The underlying question, already raised paragraph I, is the existence of a relevant shock law for blocks. The computer program UDEC, originated by Cundall, [START_REF] Cundall | A computer model for simulating progressive large scale movements of blocky rock systems[END_REF], deals with such collections of blocks. As in TRUBAL, the contact is established through a system of springs and dampers, a nd the tangential stress follows Coulomb's law. Similarly one has to enquire which normal and tangentia l stiffness and viscosity coefficients are to be adopted in an analysis to obtain a realistic frictional contact model.

DETAIL at time Is

Figure 1 .

 1 Figure 1. Local frame.

  qN '/{.= ('l{.T, 'l{.N) 'l.l= ('UT ''UN)

Figure 3 .

 3 Figure 3. Coulomb's law graph.

Figure 4

 4 Figure 4. A regularized unilateral contact law.

Figure 5

 5 Figure 5. A regularized friction law.Figure 6. Static and dynamic friction coefficients.

Figure 6 .

 6 Figure 5. A regularized friction law.Figure 6. Static and dynamic friction coefficients.

Figure 7 .Figure 8 .

 78 Figure 7. Thick graph of Signorini's condition.

Figure 9 .

 9 Figure 9. Equations of the frictional contact problem.

Signorini's condition

  Coulomb's law inelastic shock when rigid bodies are contacting Due to the positive homogeneity of Signorini's condition and Coulomb's law, the impulse h~+l may be choosen as a variable as well as the density of impulse ~+I. To summarize, + Wk(-Mk ( qk-q(I)) +hFk +hRk+l), qk+l-q(I) = h qk+l , the initial conditions, the kinematic relations 'Uk+l = Hk* qk+l, Rk+I = Hk 9ac+I , ~k+l = ~(I)+ h'UNk+l , the frictional contact laws, for a=l, ... , x.S (~<Xk+l , 9{.Nak+l) is true, C9{.Nk+l ('Ur<Xk+l, 9{.Tak+l) is true, 9{.Nk+l => 'UNak+l = 0 .Signorini's conditionCoulomb's law inelastic shock when rigid bodies are contacting 4. SOLVING THE DISCRETIZED EQUATIONS 4.1. The system of equations with local variables as primary unknowns[START_REF] Raous | Experimental analysis of the rocking of a rigid block[END_REF]

'

  VIib k = Hk*(qk + Wk( -Mk ( qk-q(I)) +hFk)).(18)The term 'l1ib k is the free relative velocity, i.e. the relative velocity of candidates when no reactions are applied, 9Gc = 0. The relation(18) shows the linear dependance between the unknowns, 'Uk+l, h9ac+I, obtained through the linearized form of the equation of Dynamics and the kinematic equations. The equation (18) is written together with Signorini's condition and Coulomb's law. The unknowns of the system are thus 'Uk+l and 9ac+l. Sub iterations are needed to solve this system. It is understood that quantities with indices k have provisional if Ulii>N :s; 0 and fT-llfN :s; 0 then h~T = Jl h~N , h~N = -UnbN I (1-Jla) UJm , i.e backward sliding; if UnbN :s; 0 and fT + Jl!N :s; 0 et fT-Jl!N ~ 0 then h~T = -h , h~N = -!N , i.e. contact without sliding.

Figure 11 .

 11 Second invariant of the deviatoric stress; scale from 20 to 220 MPa.

Figure 12 .

 12 Figure 12. Second invariant of the deviatoric strain; scale from 0.01 to 0.09.

Figure 14 .

 14 Figure 14. Temple horizontally shakcd.

values which are known. In the following, the indices k, k+l, shall be omitted to simplify. The system to be solved is then

Rewriting Signorini's condition

Signorini's condition may be written UNa 2:: 0 ~a 2:: 0 UNa~a= 0 where UNa = ~a(l)/h +'UNa .

When the candidate a defines a contact between two rigid bodies, the form (Ibis), paragraph 1.6, equivalent to Signorini's condition together with the inelastic shock law is used: if a contact is not expected, ~a = 0 if a contact is expected, UNa 2:: 0 ~a 2:: 0 UNa ~a= 0 where UNa = 'UNa .

A criterion has to be chosen, in order to decide if a contact is expected. Two possible criteria are: i) A value of the free gap, i.e. the predicted gap when reactions are not exerted is If this free gap is negative a contact is assumed to be expected. A prediction at half the time step may also be used. This criterion decides explicitly the status of contact ii) A value of the free gap, i.e. the predicted gap when the reaction on the candidate a is null, and when the reactions on other candidates ~;ta have provisional values, is If this free gap is negative a contact is assumed to be expected. This criterion, of the implicit kind, is examined each time values of '1( a are being computed.

To summarize, the following form of Signorini's condition is written According to the definition of uNa' in the case where one of the contacting body is deformable this relation accounts for the usual Signorini's condition S (qNa.' '!(NO.), and when the two contacting bodies are rigid, it accounts for the usual Signorini's condition together with the inelastic shock law .

The solution for the case of a single contact

The case of a single contact is first discussed. The index a is omitted here. One definesU = (UT'UN ), UT = 'l.JT, UlibT = 'l.iiibr and UN is defined as in 4.2: UN= ~ = llN(I)!h + 'UN, UiibN = llN(I)Ih + 'l.iiibN when one of the contacting body is deformable; UN= 'UN, UiihN = 'lllibN when two rigid bodies are contacting. Using the forms ( 4) and ( 9) of Signorini's condition and Coulomb's law, the system may be written, u -Unb -'Wh'l( = o . 'l(N -proj R+ (~ -pUN) = 0, 'l(T -proj C ('l(T-p UT) = 0, C = [jl'l(N'-jl'l(N] in the 2-dirnensional case, C is the disk with center 0 and radius jl'l(N in the 3-dirnensional case.

This system has the form <l>(X)=O, where X=( U ,1(), and <l> is a piecewise continuous linear function. In the 2-dimensional case, the solution may be exhibited in a straightforward manner: one sets It is assumed -I < jla < I . Then the system has a unique solution: if UlibN > 0 then h'l(T = 0 , h'l(N = 0 , i.e. no contact if UlibN ~ 0 and fT + llfN 2: 0 then h'l(T = -jl h'l(N , h'l(N = -UlibN I (I +jla) ~ , i.e. forward sliding; (degree) crl/cr2 2.5 0.3 void rate