W P M H Heemels

M K Camlibel

B Brogliato

J M Schumacher

Observer-Based Control of Linear Complementarity Systems

In this paper, we will present observer and output-based controller design methods for linear complementarity systems (LCS) employing a passivity approach. Given various inherent properties of LCS, such as the presence of state jumps, mode dynamics described by DAEs, and regions ("invariants") for certain modes being lower dimensional, several proposed observers and controllers for other classes of hybrid dynamical systems do not apply. We will provide sufficient conditions for the observer design for a LCS, which is effective also in the presence of state jumps. Using the certainty equivalence approach we obtain output-based controllers for which we will derive a separation principle.

Introduction

Complementarity systems form a class of hybrid dynamical systems that received considerable attention in recent years [START_REF] Van Der Schaft | Complementarity modeling of hybrid systems[END_REF][START_REF] Heemels | Linear complementarity systems[END_REF][START_REF] Camlibel | On linear passive complementarity systems[END_REF][START_REF] Shen | Semicopositive lineaer complementarity systems[END_REF][START_REF] Camlibel | Complementarity methods in the analysis of piecewise linear dynamical systems[END_REF][START_REF] Brogliato | Some perspectives on the analysis and control of complementarity systems[END_REF][START_REF] Brogliato | On the equivalence between complementarity systems, projected systems and differential inclusions[END_REF][START_REF] Heemels | Projected dynamical systems in a complementarity formalism[END_REF][START_REF] Camlibel | Switched networks and complementarity[END_REF][START_REF] Camlibel | Consistency of a time-stepping method for a class of piecewise-linear networks[END_REF]. The linear complementarity system LCS(A, B, C, D, E, F ) is given by ẋ(t) = Ax(t) + Bw(t) + Eu(t) (1a)

z(t) = Cx(t) + Dw(t) + F u(t) ( 1b) 0 ≤ z(t) ⊥ w(t) ≥ 0, (1c)
where the inequalities are interpreted componentwise and ⊥ indicates the orthogonality between the vectors z(t) and w(t), i.e, z (t)w(t) = 0. The complementarity conditions (1c) constitute a particular set of equalities and inequalities, which are related to the well-known relations between the constraint variables and Lagrange multipliers in the Karush-Kuhn-Tucker conditions for optimality, the voltage-current relationship of ideal diodes, the conditions between unilateral constraints and reaction forces in constrained mechanics, etc. As such, the complementarity framework includes mechanical systems with unilateral constraints, constrained optimal control problems, switched electrical circuits, etc.

Although LCS has its own peculiarities, it has connections to other classes of hybrid systems. Indeed, observe that (1c) implies that w i (t) = 0 or z i (t) = 0 for all i ∈ m := {1, . . . , m}. As a consequence, the system (1) has 2 m modes. Each mode can be characterized by the active index set J ⊆ m, which indicates 1 z i = 0, i ∈ J, and w i = 0, i ∈ J c , where J c := {i ∈ m | i ∈ J}. For mode J the dynamics is given by the linear differential and algebraic equations (DAEs) ẋ(t) = Ax(t) + Bw(t) + Eu(t), (2a) z(t) = Cx(t) + Dw(t) + F u(t), (2b) zi(t) = 0, i ∈ J, and wi(t) = 0, i ∈ J c , (2c)

The evolution of system (1) will be governed by [START_REF] Heemels | Linear complementarity systems[END_REF] for mode J as long as the remaining inequalities ("the invariant" in the terminology of hybrid automata [START_REF] Lynch | Hybrid I/O automata[END_REF][START_REF] Branicky | A unified framework for hybrid control: model and optimal control theory[END_REF][START_REF] Lygeros | Dynamical properties of hybrid automata[END_REF]) in (1c) zi(t) ≥ 0, i ∈ J c and wi(t) ≥ 0, i ∈ J [START_REF] Camlibel | On linear passive complementarity systems[END_REF] are satisfied. Impending violation of (3) will trigger a mode change. As a consequence, during the evolution in time of the system several mode dynamics will be active successively. This indicates that LCS might be recast within the hybrid automaton framework [START_REF] Lynch | Hybrid I/O automata[END_REF][START_REF] Branicky | A unified framework for hybrid control: model and optimal control theory[END_REF][START_REF] Lygeros | Dynamical properties of hybrid automata[END_REF]. However, with exception only of the very simplest cases, the reformulation of LCS dynamics into the hybrid automaton framework leads to voluminous and opaque system descriptions. This effect is already evident in the example worked out in [START_REF] Heemels | Modelling, wellposedness, and stability of switched electrical networks[END_REF], which concerns an electrical circuit with two diodes. Alternatively one can rewrite LCS in the formulation ẋ ∈ F (x) for x ∈ C and x + ∈ G(x) for x ∈ D, as used for instance in [START_REF] Cai | Smooth Lyapunov functions for hybrid systems. Part I: Existence is equivalent to robustness[END_REF]. Again, the transformation is in general cumbersome, and the resulting system data F , G, C and D do generally not satisfy the assumptions adopted in [START_REF] Cai | Smooth Lyapunov functions for hybrid systems. Part I: Existence is equivalent to robustness[END_REF] (cf. the example below). Structural properties become harder to study and compactness of descriptions is lost when LCS are translated to such generic frameworks.

To show further links between LCS and other (sub)classes of hybrid dynamical systems, let us consider an LCS with one complementarity pair and F = 0:

ẋ = Ax + bw + eu; z = c x + dw; 0 ≤ z ⊥ w ≥ 0 ( 4
)
where

A ∈ R n×n , b ∈ R n×1 , c ∈ R n×1 , d ∈ R
, and 0 = e ∈ R n×1 . As either z = 0 or w = 0, this system has two modes. When d > 0, one can rewrite (4) as

ẋ = Ax + eu if c x ≥ 0, (A -bd -1 c )x + eu if c x ≤ 0, (5)
which is a piecewise linear (PWL) system [START_REF] Sontag | Nonlinear regulation: The piecewise linear approach[END_REF]. When d = 0 and c b > 0, we have

ẋ = Ax + eu if (c x > 0) or (c x = 0 and c Ax + c eu ≥ 0), P (Ax + eu) if c x = 0 and c Ax + c eu < 0. ( 6
)
where P = Ib(c b) -1 c . In this case one has also a bimodal PWL system, but the second subsystem 'lives' on a lower dimensional subspace given by c x = 0, which is a situation hardly studied within the realm of PWL systems. Note also that for c x < 0, there is no smooth evolution possible and state jumps will occur. In this situation, LCS can also be considered as differential inclusions (DIs) with normal cones as their set values (see e.g. [START_REF] Brogliato | On the equivalence between complementarity systems, projected systems and differential inclusions[END_REF]). These DIs do in general not satisfy the boundedness conditions of its values nor the upper semicontinuity properties as often used within the context of DIs. In the case when d = c b = 0, c ab > 0 and e = 0 (no external inputs) the flow set, i.e. the set of states from which the system can continue with a smooth solution temporarily [START_REF] Cai | Smooth Lyapunov functions for hybrid systems. Part I: Existence is equivalent to robustness[END_REF], is given by all x 0 such that (c x 0 , c Ax 0 ) 0 (see [START_REF] Heemels | Linear complementarity systems[END_REF]Thm. 6.8]), where denotes the lexicographic ordering. This indicates that the flow set has no simple closedness properties. This indicates that although LCS have connections to PWL systems and other classes of hybrid systems, they also have their own peculiarities. For instance, the presence of state jumps (think of impacts in constrained mechanical systems) in continuous-time LCS, differentiates them from much of the work done for continuous-time PWL systems. Although for discrete-time LCS strong equivalence links have been established in [START_REF] Heemels | Equivalence of hybrid dynamical models[END_REF] with piecewise affine systems [START_REF] Sontag | Nonlinear regulation: The piecewise linear approach[END_REF] and other classes of hybrid models such as min-max-plus-scaling systems [START_REF] De Schutter | On model predictive control for max-min-plusscaling discrete event systems[END_REF] and mixed logic dynamic systems [START_REF] Bemporad | Control of systems integrating logic, dynamics, and constraints[END_REF], in the continuous-time framework, which is the natural habitat for most of the LCS applications, such broad equivalence relations are out of the question. There are relations though of LCS to other specific classes of nonsmooth systems such as the mentioned "normal cone DIs" and projected dynamical systems [START_REF] Brogliato | On the equivalence between complementarity systems, projected systems and differential inclusions[END_REF][START_REF] Heemels | Projected dynamical systems in a complementarity formalism[END_REF].

The attention that LCS received recently is not surprising given the broad range of interesting applications. The research [START_REF] Van Der Schaft | Complementarity modeling of hybrid systems[END_REF][START_REF] Heemels | Linear complementarity systems[END_REF][START_REF] Camlibel | On linear passive complementarity systems[END_REF][START_REF] Shen | Semicopositive lineaer complementarity systems[END_REF][START_REF] Camlibel | Complementarity methods in the analysis of piecewise linear dynamical systems[END_REF][START_REF] Brogliato | Some perspectives on the analysis and control of complementarity systems[END_REF][START_REF] Brogliato | On the equivalence between complementarity systems, projected systems and differential inclusions[END_REF][START_REF] Heemels | Projected dynamical systems in a complementarity formalism[END_REF][START_REF] Camlibel | Switched networks and complementarity[END_REF][START_REF] Camlibel | Consistency of a time-stepping method for a class of piecewise-linear networks[END_REF] focussed on several fundamental system-theoretic issues like well-posedness, discretization (simulation), controllability, observability, stabilizability and stability. In this paper the emphasis will be on observer-based controller design, a topic that is hardly touched upon for LCS. We will adopt a "certainty equivalence control" approach, where one designs output feedback controllers that generate the control input via a state feedback law using an estimate of the state, which is obtained from an observer. For linear systems, the separation principle gives a formal justification of this approach. Due to the absence of a general controllerobserver separation principle for nonlinear systems and certainly hybrid systems, the observer-based controller may not result in a stable closed-loop system.

Several interesting papers are available on observer design for hybrid systems, especially in the context of switched and piecewise linear systems, e.g. [START_REF] Alessandri | Design of Luenberger observers for a class of hybrid linear systems[END_REF][START_REF] Iulia Bara | State estimation for a class of hybrid systems[END_REF][START_REF] Petterson | Switched state jump observers for switched systems[END_REF][START_REF] Juloski | Observer design for a class of piecewise linear systems[END_REF][START_REF] Pavlov | Convergent piecewise affine systems: analysis and design[END_REF]. Unfortunately, these results do not apply to LCS as LCS typically exhibit lower dimensional regions and state jumps. Observer and observer-based controller design methods for Lur'e type systems as studied in [START_REF] Arcak | Observer based control of systems with slope-restricted nonlinearities[END_REF][START_REF] Arcak | Certainty-equivalence output-feedback design with circle-criterion observers[END_REF][START_REF] Fan | Observer design for systems with multivariable monotone nonlinearities[END_REF][START_REF] Rajamani | Observers for Lipschitz nonlinear systems[END_REF] are also related to LCS. Indeed, one can consider LCS as a kind of Lur'e type systems in which the linear system (1a)-(1b) is interconnected with the nonsmooth and unbounded complementarity relations in (1c). Typically, the results in [START_REF] Arcak | Observer based control of systems with slope-restricted nonlinearities[END_REF][START_REF] Arcak | Certainty-equivalence output-feedback design with circle-criterion observers[END_REF][START_REF] Fan | Observer design for systems with multivariable monotone nonlinearities[END_REF][START_REF] Rajamani | Observers for Lipschitz nonlinear systems[END_REF] study locally Lipschitz slope restricted nonlinearities in the feedback path. As such their conditions do not allow for the non-smoothness and setvalued nonlinearities (and even state jumps) as induced by the complementarity relations. Observer designs for differential inclusions with bounded set-values are treated in [START_REF] Osorio | Dissipative design of observers for multivalued nonlinear systems[END_REF]. Since complementarity conditions are unbounded, [START_REF] Osorio | Dissipative design of observers for multivalued nonlinear systems[END_REF] does not cover LCS.In summary, although there are various interesting approaches for hybrid observer design, none of them includes (all) the peculiarities of LCS.

Before focussing on controller and observer design for LCS, we will explain the solution concepts for LCS and also extend available well-posedness (existence and uniqueness of solutions) results for LCS to include external inputs. This will result in global existence results thereby excluding the Zeno phenomenon of livelock (an infinite number of discrete actions on one time instant) and providing continuations beyond accumulation points of mode switching times. Next, we will present passivity conditions for state feedback design and observer design for LCS. Interestingly, this means that we will present methods for observer design for systems without knowing the mode and allowing for state resets. Next we will present a separation principle for this class of hybrid dynamical systems.

Preliminaries

R denotes the real numbers, R + := [0, ∞) the nonnegative real numbers, L 2 (T ) the square integrable functions on a time-interval T ⊆ R, and B the Bohl functions (i.e. functions having strictly proper rational Laplace transforms) defined on R + . Note that sines, cosines, exponentials, polynomials and their sums and products are all Bohl functions. The distribution δ (i) t stands for the i-th distributional derivative of the Dirac impulse supported at t. The dual cone of a set Q ⊆ R n is defined by Q * = {x ∈ R n |x y ≥ 0 for all y ∈ Q}. For a positive integer m, the set m is defined as {1, 2, . . . , m} and 2 m denotes the collection of all subsets of m. A vector u ∈ R k is called nonnegative, denoted by u ≥ 0, if u i ≥ 0 for all i ∈ k. This means that inequalities for vectors are interpreted componentwise. The orthogonality u y = 0 between two vectors u ∈ R k and y ∈ R k is denoted by u ⊥ y. As usual, we say that (A, B, C)

(or sometimes (A, B, C, D)) is minimal, when the matrices [B AB . . . A n-1 B] and [C A C . . . (A ) n-1 C ] have full rank. A matrix M ∈ R k×k is called positive definite (not necessarily symmetric), if x M x > 0 for all x = 0. It is called nonnegative definite, if x M x ≥ 0 for all x ∈ R k . For a matrix M ∈ R k×l we denote its kernel by kerM := {x ∈ R l | M x = 0} and its im- age by imM := {M x | x ∈ R l }. Finally, for two linear subspaces V 1 and V 2 we write V 1 ⊕ V 2 = V, if V = V 1 + V 2 = {v 1 + v 2 | v 1 ∈ V 1 , v 2 ∈ V 2 } and V 1 ∩ V 2 = {0}.

Linear Complementarity Problem

We define the linear complementarity problem LCP(q, M ) (see [START_REF] Cottle | The Linear Complementarity Problem[END_REF] for a survey) with data q ∈ R k and M ∈ R k×k by the problem of finding z ∈ R k such that 0 ≤ z ⊥ q + M z ≥ 0. The solution set of LCP(q, M ) will be denoted by SOL(q, M ). The notation K(M ) will denote the set {q | LCP(q, M ) is solvable}.

Let a matrix M of size k×k and two subsets I and J of k of the same cardinality be given. The (I, J)-submatrix of M is the submatrix M IJ := (M ij ) i∈I,j∈J . The (I, J)-minor is defined as the determinant of M IJ . The (I, I)-submatrices and -minors are also known as the principal submatrices and the principals minors. M is called a P-matrix if all its principal minors are positive. P-matrices play an important role in linear complementarity problems, as the following result is well known (cf. [30, thm. 3.3.7]).

Theorem 1. For a given matrix M ∈ R k×k , the problem LCP(q, M ) has a unique solution for all vectors q ∈ R k if and only if M is a P-matrix.

Passivity of a Linear System

We will recall the notion of passivity as it is defined in [START_REF] Willems | Dissipative dynamical systems[END_REF] for a linear system

Σ(A, B, C, D) : ẋ(t) = Ax(t) + Bu(t); y(t) = Cx(t) + Du(t). ( 7
)
Definition 1. [START_REF] Willems | Dissipative dynamical systems[END_REF] The system Σ(A, B, C, D) given by ( 7) is said to be passive if there exists a function V : R n → R + (called a storage function), such that

V (x(t0)) + t 1 t 0 u (t)y(t)dt ≥ V (x(t1))
holds for all t 1 ≥ t 0 , and for all solutions (u, x, y) ∈ L m+n+m 2 (t 0 , t 1 ) of [START_REF] Brogliato | On the equivalence between complementarity systems, projected systems and differential inclusions[END_REF].

Next, we quote a very well-known characterization of passivity. Moreover, V (x) = 1 2 x P x defines a quadratic storage function if and only if P satisfies the linear matrix inequalities [START_REF] Heemels | Projected dynamical systems in a complementarity formalism[END_REF].

Next we will define strict passivity of (7). A,B,C,D) is called strictly passive, if there are P and ε > 0 with P = P > 0 and

Definition 2. Σ(

A P + P A + εP P B -C B P -C -(D + D ) ≤ 0. ( 9
)
The following technical assumption will be used below. Its former part is standard in the literature on dissipative systems, see e.g. [START_REF] Willems | Dissipative dynamical systems[END_REF]. Both conditions are related to removing specific kinds of redundancy in the system description.

Linear Complementarity Systems and Initial Solution

We now are going to study the behaviour of the LCS given in [START_REF] Van Der Schaft | Complementarity modeling of hybrid systems[END_REF]. Before presenting a global solution concept that incorporates the switching of modes, we will first concentrate on what we call "initial solutions," which are trajectories that satisfy the dynamics of one mode only and that satisfy the inequality conditions possibly only in their beginning. We will employ the theory of distributions in formalizing the solution concept, since the abrupt changes in the trajectories (e.g. impacts in mechanics) can be modeled adequately by Dirac impulses (see also Example 1 below). To do so, we recall the concepts of Bohl distribution and initial solution [START_REF] Heemels | Linear complementarity systems[END_REF].

Definition 3. We call w a Bohl distribution, if w = w imp + w reg with w imp = l i=0 w -i δ (i)
0 for w -i ∈ R and w reg ∈ B. We call w imp the impulsive part of w and w reg the regular part of w. The space of all Bohl distributions is denoted by B imp .

Note that Bohl distributions have rational Laplace transforms. It seems natural to call a (smooth) Bohl function w ∈ B initially nonnegative if there exists an ε > 0 such that w(t) ≥ 0 for all t ∈ [0, ε). Note that a Bohl function w is initially nonnegative if and only if there exists a σ 0 ∈ R such that its Laplace transform satisfies ŵ(σ) ≥ 0 for all σ ≥ σ 0 . Hence, there is a connection between small time values for time functions and large values for the indeterminate s in the Laplace transform. This fact is closely related to the well-known initial value theorem. The definition of initial nonnegativity for Bohl distributions will be based on this observation (see also [START_REF] Heemels | Linear complementarity systems[END_REF]). Definition 4. We call a Bohl distribution w initially nonnegative, if its Laplace transform ŵ(s) satisfies ŵ(σ) ≥ 0 for all sufficiently large real σ.

To relate the definition to the time domain, note that a scalar-valued Bohl distribution w without derivatives of the Dirac impulse (i.e. w imp = w 0 δ for some w 0 ∈ R) is initially nonnegative if and only if either w 0 > 0, or w 0 = 0 and there exists an ε > 0 such that w reg (t) ≥ 0 for all t ∈ [0, ε). With these notions we can recall the concept of an initial solution [START_REF] Heemels | Linear complementarity systems[END_REF]. Loosely speaking, an initial solution to (1) with initial state x 0 and Bohl input u ∈ B k is a triple (w, x, z) ∈ B m+n+m imp satisfying (2) for some mode I and satisfying (3) either on a time interval of positive length or on a time instant at which Dirac distributions are active.

At this point we only allow Bohl functions (combinations of sines, cosines, exponentials and polynomials) as inputs. In the global solution concept we will allow the inputs to be concatenations of Bohl functions (i.e., piecewise Bohl), which may be discontinuous. Definition 5. The distribution (w, x, z) ∈ B m+n+m imp is said to be an initial solution to [START_REF] Van Der Schaft | Complementarity modeling of hybrid systems[END_REF] with initial state x 0 and input u ∈ B k if 1. ẋ = Ax+Bw+Eu+x 0 δ 0 and z = Cx+Dw+F u as equalities of distributions.

2. there exists an J ⊆ m such that w i = 0, i ∈ I c and z i = 0, i ∈ J as equalities of distributions. 3. w and z are initially nonnegative.

The statements 1 and 2 in the definition above express that an initial solution satisfies the dynamics (2) for mode J on the time-interval R + and the initial condition x(0) = x 0 .

Example 1. Consider the system ẋ(t) = w(t), z(t) = x(t) together with (1c). This represents an electrical network consisting of a capacitor connected to a diode. The current is equal to w and the voltage across the capacitor is equal to z = x. For initial state x(0) = 1, (w, x, z) with w = 0 and z(t) = x(t) = 1 for all t ∈ R is an initial solution. This corresponds to the case that the diode is blocking and there is no (nonzero) current in the network. To show that the distributional framework is convenient, consider the initial state x(0) = -1, for which (w, x, z) with w = δ, x(t) = z(t) = 0, t > 0 is the unique initial solution. This corresponds to an instantaneous discharge of the capacitor at time instant 0 resulting in a state jump from x(0) = -1 to 0 at time 0 induced by the impulse.

This circuit example indicates that there is a clear physical interpretation of the impulses in initial solutions to model the abrupt changes in trajectories of electrical circuits (see also [START_REF] Camlibel | On linear passive complementarity systems[END_REF][START_REF] Heemels | Modelling, wellposedness, and stability of switched electrical networks[END_REF]). Also for mechanical systems, these initial solutions induce state jumps with a clear physical meaning (related to inelastic restitution laws) as shown in [START_REF] Heemels | Linear complementarity systems[END_REF].

Initial and Local Well-Posedness

In this section, we are interested in existence and uniqueness of initial solutions, and we will provide an extension to a local well-posedness result. Definition 6. A rational matrix H(s) ∈ R l×l (s) is said to be of index r, if it is invertible as a rational matrix and s -r H -1 (s) is proper. It is said to be totally of index r, if all its principal submatrices H JJ (s) for J ⊂ l are of index r.

The theorem below is an extension of a result proven in [START_REF] Camlibel | Complementarity methods in the analysis of piecewise linear dynamical systems[END_REF] that used F = 0. [START_REF] Van Der Schaft | Complementarity modeling of hybrid systems[END_REF]

Theorem 4. Consider an LCS with external inputs given by

such that G(s) = C(sI -A) -1 B + D is totally of index 1 and G(σ) is a P-matrix for sufficiently large σ. Define Q D := SOL(0, D) = {v ∈ R m | 0 ≤ v ⊥ Dv ≥ 0} and let K(D) be the set {q ∈ R m | LCP (q, D) solvable}.
1. For arbitrary initial state x 0 ∈ R n and any input u ∈ B k , there exists exactly one initial solution, which will be denoted by (w x0,u , x x0,u , z x0,u ).

No initial solution contains derivatives of the Dirac distribution. Moreover,

w x0,u imp = w 0 δ 0 ; x x0,u imp = 0; z x0,u imp = Dw 0 δ 0 for some w 0 ∈ Q D .
3. For all x 0 ∈ R n and u ∈ B k it holds that Cx 0 + F u(0) + CBw 0 ∈ K(D).

4. The initial solution (w x0,u , x x0,u , z x0,u ) is smooth (i.e., has a zero impulsive part) if and only if Cx 0 + F u(0) ∈ K(D).

The impulsive part w x0,u imp = w 0 δ 0 of the initial solution induces a state jump from x 0 to x 0 + Bw 0 , which in circuits might correspond to infinitely large currents related to instantaneous discharges of capacitors (see Example 1) and in mechanical systems to infinite reaction forces that cause resets in the velocities of impacting bodies. For brevity we skip the proof (see e.g. [START_REF] Camlibel | On linear passive complementarity systems[END_REF] for the case F = 0.) Theorem 4 gives explicit conditions for existence and uniqueness of solutions.The second statement indicates that derivatives of Dirac distributions are absent in the behaviour of LCS of index 1. The fourth statement gives necessary and sufficient condition for an initial solution to be smooth. In particular, a LCS satisfying the conditions of Theorem 5 is "impulse-free" (no state jumps), if Q D =SOL(0, D) = {0} (or, in terms of [START_REF] Cottle | The Linear Complementarity Problem[END_REF], if D is an R 0 -matrix). Note that in this case Q * D = R m . In case the matrix [C F] has full row rank, this condition is also necessary. Other sufficient conditions, that are more easy to verify, are D being a positive definite matrix, or Ker(D + D ) ∩ R m + = {0}. For the general case of (not necessarily passive) systems with transfer functions totally of index 1, Theorem 4 implies that K(D) = R m (in terms of [START_REF] Cottle | The Linear Complementarity Problem[END_REF] this means that the matrix D should be a so-called Q-matrix) is sufficient for the system being impulse-free. One condition that guarantees this is that w T Dw > 0 for all w ∈ R m + and w = 0. Note that the first statement in Theorem 4 by itself does not immediately guarantee the existence of a solution on a time interval with positive length. The reason is that an initial solution with a non-zero impulsive part may only be valid at the time instant on which the Dirac distribution is active. If the impulsive part of the (unique) initial solution is equal to w 0 δ 0 , the state after re-initialization is equal to x 0 +Bw 0 . From this "next" initial state again an initial solution has to be determined, which might in principle also have a non-zero impulsive part, which results in another state jump. As a consequence, the occurrence of infinitely many jumps at t = 0 without any smooth continuation on a positive length time interval is not immediately excluded (sometimes called "livelock" in hybrid systems theory). However, Theorem 4 excludes this kind of Zeno phenomenon: if smooth continuation is not directly possible from x 0 , it is possible after one re-initialization. Indeed, since C(x 0 + Bw 0 ) + F u(0) = Cx 0 + F u(0) + CBw 0 ∈ K(D), it follows from the fourth claim that the initial solution corresponding to x 0 + Bw 0 and input u is smooth. This initial solution satisfies the (in)equalities in [START_REF] Van Der Schaft | Complementarity modeling of hybrid systems[END_REF] on an interval of the form (0, ε) with ε > 0 by definition and hence, we proved a local existence and uniqueness result. However, we still have to show global existence of solutions as other kinds of Zeno behaviour (accumulation of mode switching times) might prevent this.

Global Well-Posedness

Before we can formulate a global well-posedness theorem, we need to define a class of allowable input functions and the global solution concept.

Definition 7. A function u : R + → R is called a piecewise Bohl function, if -u is right-continuous, i.e. lim t↓τ u(t) = u(τ ) for all τ ∈ R + -for all τ ∈ R + there are ε > 0 and v ∈ B such that u(t) = v(t) for all t ∈ (τ, τ + ε) -u is locally bounded in the sense that for any interval [0, T ] there is a constant M such that |u(t)| ≤ M for all t ∈ [0, T ].
We denote this function space as PB. A distribution u : R + → R is called piecewise Bohl with first order impulses, if u is the sum of a piecewise Bohl function u reg and the distribution u imp = θ∈Γ w θ δ θ , where Γ = {τ i } i is a finite or countable subset of R + , which is isolated1 . We denote this distribution space by PB 0 . Definition 8. Let (w, x, z) ∈ PB m+n+m 0 be given with (absolutely) continuous and satisfies [START_REF] Van Der Schaft | Complementarity modeling of hybrid systems[END_REF] for almost all t ∈ (a, b) 2. For each θ ∈ Γ the corresponding impulse (w θ δ θ , x θ δ θ , z θ δ θ ) is equal to the impulsive part of the unique initial solution2 to (1) with initial state x reg (θ-) := lim t↑θ x reg (t) (taken equal to x 0 for θ = 0) and input t → u(tθ). 3. For times θ ∈ Γ it holds that x reg (θ+) = x reg (θ-) + Bw θ with w θ the multiplier of Dirac pulse supported at θ. Theorem 6. Consider an LCS with external inputs given by [START_REF] Van Der Schaft | Complementarity modeling of hybrid systems[END_REF] Proof. The proof follows along similar lines as the proof of the case F = 0 as given in the thesis [START_REF] Camlibel | Complementarity methods in the analysis of piecewise linear dynamical systems[END_REF] by carefully incorporating the presence of impulses (see also [START_REF] Camlibel | On linear passive complementarity systems[END_REF] for the passive case).

w imp = θ∈Γ w θ δ θ , x imp = θ∈Γ x θ δ θ , z imp = θ∈Γ z θ δ θ for w θ ∈ R m , x θ ∈ R n and w θ ∈ R m for θ ∈ Γ
such that G(s) = C(sI -A) -1 B + D is totally of index 1 and G(σ)
This theorem implies that if F u is continuous, jumps of the state can only occur at the initial time instant 0. Global existence: for each x 0 there exists a global solution to (1) and moreover, for each T ≥ 0 all solutions (w, x, z) ∈ PB m+n+m 0 to (1) for initial state x 0 defined on [0, T ) can be continued to a global solution on [0, ∞); Lyapunov stability: for each ε > 0 there exists a δ > 0 such that x reg (t) < ε, for all t ≥ 0 when x 0 < δ, where (w, x, z) ∈ PB m+n+m 0 is a global solution in the sense of Definition 8 to (1) for initial state x 0 ; Attractiveness: lim t→∞ x reg (t) = 0 for any global solution (w, x, z) ∈ PB 0 .

In this section we aim at designing a state feedback controller u(t) = Kx(t) [START_REF] Camlibel | Consistency of a time-stepping method for a class of piecewise-linear networks[END_REF] that renders the system (1) GAS. For necessary and sufficient conditions of "passifiability" for the case D = 0 (i.e. finding K such that (A + EK, B, C + F K, 0) is strictly passive), see [START_REF] Arcak | Observer based control of systems with slope-restricted nonlinearities[END_REF].

Theorem 8. Consider the LCS (1) and the state feedback [START_REF] Camlibel | Consistency of a time-stepping method for a class of piecewise-linear networks[END_REF] and suppose that Assumptions 3 and 7 hold. Then the closed-loop system (1)-( 10) is GAS.

In [START_REF] Camlibel | On linear passive complementarity systems[END_REF] it was shown that a sufficient condition for GAS is the strict passivity of the underlying system. The above theorem on state feedback design is a consequence of that result. After these preparatory steps, we continue with the main results of this paper related to observer design and output-based controller design.

Observer Design

Consider the LCS (1) and assume that only

y(t) = Gx(t) ∈ R p (11)
is measured instead of the complete state being available for feedback. Based on this output measurement, we aim at estimating the continuous state x(t) of ( 1) using an observer. For the observer design to be meaningful, we have to assume some conditions of the existence of solutions to the observed system (1).

Assumption 9. For any initial state x 0 and for any input function u ∈ PB k there exists a global solution (w, x, z) ∈ PB m+n+m 0 to system (1) in the sense of Definition 8.

The well-posedness theory derived before can be used to guarantee this property. For instance, Theorem 6 shows that under the assumption that G(s) = C(sI -A) -1 B + D is totally of index 1 and G(σ) is a P-matrix for sufficiently large σ, Assumption 9 is indeed satisfied.

We propose the following observer for LCS [START_REF] Van Der Schaft | Complementarity modeling of hybrid systems[END_REF] with measured output [START_REF] Lynch | Hybrid I/O automata[END_REF]:

ẋ(t) = Ax(t) + B ŵ(t) + Eu(t) + L(y(t) -ŷ(t)) (12a) ẑ(t) = C x(t) + D ŵ(t) + F u(t) + M (y(t) -ŷ(t)) (12b) 0 ≤ ẑ(t) ⊥ ŵ(t) ≥ 0. (12c) ŷ(t) = Gx(t), ( 12d
)
where we have two observer gains (L and M ).

Assumption 10. (A -LG, B, C -M G, D) is strictly passive and minimal.

First of all, one has to show that the observer structure [START_REF] Branicky | A unified framework for hybrid control: model and optimal control theory[END_REF] Since we proved global existence of x, we can consider the observation error e := xx, which evolves according to the following dynamics

ė(t) = (A -LG)e(t) + Bw(t) -B ŵ(t) (13a) z(t) = Cx(t) + Dw(t) + F u(t) (13b) ẑ(t) = C x(t) + D ŵ(t) + F u(t) + M (y(t) -ŷ(t)) (13c) 0 ≤z(t) ⊥ w(t) ≥ 0, and 0 ≤ ẑ(t) ⊥ ŵ(t) ≥ 0 (13d)
Theorem 12. Consider the error dynamics [START_REF] Lygeros | Dynamical properties of hybrid automata[END_REF] such that Assumption 3, Assumption 9 and Assumption 10 hold. Then the error dynamics is GAS 3 .

Proof. See the report [START_REF] Heemels | Observer-based control of linear complementarity systems[END_REF].

The theorem shows (under the given hypothesis) that the observer (12) recovers asymptotically the state of the LCS (1) based on the output [START_REF] Lynch | Hybrid I/O automata[END_REF], even when the state of the system or observer exhibits state jumps. Since the jumps in both the observer and the observed plant are triggered by discontinuities in the external signal (due to the low index of the underlying linear system), the time instants of the jumps coincide for the observer and controller, which is exploited in the proof. For higher index systems (e.g. mechanical systems with unilateral constraints) this property is lost, which complicates observer design significantly.

Separation Principle: Observer-Based Controller

To design an output-based controller for (1) with ( 11), we will employ a "certainty equivalence" approach by using the estimate x obtained from the designed observer in the state feedback controller, i.e. u(t) = K x(t) = Kx(t)-Ke(t). The closed-loop system consisting of the system (1), the observer ( 12) and the controller u(t) = K x(t) becomes

ẋ ė = (A + EK) -EK 0 A -LG =:A cl x e + B 0 B -B =:B cl w ŵ (14a) z ẑ = C + F K -F K C + F K -(C + F K -MG) =:C cl x e + D 0 0 D =:D cl w ŵ (14b) 0 ≤ z ẑ ⊥ w ŵ ≥ 0 (14c)
We focus on the so-called basic observer, which is the observer [START_REF] Branicky | A unified framework for hybrid control: model and optimal control theory[END_REF] with M = 0 implying that there is only an innovation (output injection) term in the differential equation (12a) and not in the complementarity relation (12b). We will return to the extended observer with M = 0 in Remark 1. We will start by proving closed-loop well-posedness.

Theorem 13. Consider the system (14) such that Assumptions 3, 7 and 10 with M = 0 are satisfied. The system [START_REF] Heemels | Modelling, wellposedness, and stability of switched electrical networks[END_REF] has for each initial condition x 0 and e 0 a unique global solution in the sense of Definition 8. Moreover, only on time 0, there can be a discontinuity in the state trajectory x.

Proof. See the report [START_REF] Heemels | Observer-based control of linear complementarity systems[END_REF].

Now we can state a separation principle for LCS.

Theorem 14. [Separation principle]

Consider the closed-loop LCS [START_REF] Heemels | Modelling, wellposedness, and stability of switched electrical networks[END_REF]. If Assumptions 3, 7 and 10 with M = 0 are satisfied, then the LCS (14) is GAS.

Proof. See the report [START_REF] Heemels | Observer-based control of linear complementarity systems[END_REF].

Once the observer [START_REF] Branicky | A unified framework for hybrid control: model and optimal control theory[END_REF] is included in an observer-based control configuration for the LCS (1), jumps in the state variable of ( 14) can only take place at the initial time (under the given hypothesis of the theorem above). The reason is that for the 'open' LCS (1) the state jumps are triggered by the external signal u, while the closed-loop system ( 14) is a 'closed' system without external inputs. Note, however, that when the observer is applied to an 'open' LCS [START_REF] Van Der Schaft | Complementarity modeling of hybrid systems[END_REF], the state of (1) is still recovered asymptotically (under the hypothesis of Theorem 12) even when state jumps (triggered by discontinuities in u) remain to be persistently present. The fact that jumps only occur at the initial time is related to the low index of the underlying linear system of the LCS [START_REF] Heemels | Modelling, wellposedness, and stability of switched electrical networks[END_REF]. In general for systems of higher index (like constrained mechanical systems) discontinuities are not only externally triggered by exogenous signals, but also by internal events (impacts).

Remark 1. The extended observer case (i.e. M = 0) can be covered in a similar manner as above under the assumption that F = 0.

Conclusion

We presented observer and output-based controller design methods for linear complementarity systems (LCS) employing a passivity approach. We provided sufficient conditions for the observer design for a LCS, which is effective also in the presence of state jumps. Using the certainty equivalence approach we obtained output-based controllers for which we provided a separation principle in case the basic observer ("M = 0") is used or there is no direct feedthrough of the input in the complementarity conditions ("F = 0"). Future work will involve the study of the full separation principe for both F and M nonzero. Another important line of future research is the observer and observer-based control design of LCS for which the underlying linear system is of higher index (such as constrained mechanical systems.)

Theorem 2 .

 2 [START_REF] Willems | Dissipative dynamical systems[END_REF] Assume that (A, B, C) is minimal. Then the following statements are equivalent: 1. Σ(A, B, C, D) is passive. 2. The following matrix inequalities have a solution P = P > 0 and A P + P A P B -C B P -C -(D + D ) ≤ 0. (8)

Assumption 3 .

 3 (A, B, C) is minimal and B D + D has full column rank.

Theorem 5 .

 5 Consider an LCS with (A, B, C, D) is passive and Assumption 3 holds. Then G(s) = C(sI -A) -1 B + D is totally of index 1 and both G(σ) and D + s -1 CB are P-matrices for sufficiently large σ. In this case the results of Theorem 4 hold for (1) with K(D) = Q * D .

  and some Γ . Then we call (w, x, z) a (global) solution to LCS[START_REF] Van Der Schaft | Complementarity modeling of hybrid systems[END_REF] with input function u ∈ PB and initial state x 0 , if the following properties hold. 1. For any interval (a, b) such that (a, b) ∩ Γ = ∅ the restriction x reg | (a,b) is

  is a P-matrix for sufficiently large σ. The LCS (1) has a unique (global) solution (w, x, z) ∈ PB k+n+k 0 for any initial state x 0 and input u ∈ PB k . Moreover, x imp = 0 and impulses in (w, z) can only show up at the initial time and times for which F u is discontinuous (i.e. Γ in Definition 8 can be taken as a subset of {0} ∪ Γ d F u ).

Assumption 7 .

 7 (A + EK, B, C + F K, D) is strictly passive and minimal.

6

  Stability and State Feedback Design for LCS Let us first define formally what we mean by stability of LCS.

	Definition 9. The LCS (1) without inputs (i.e. E = F = 0) is called globally
	asymptotically stable (GAS), if

The set Γ ⊆ R is called isolated, if for all τ ∈ Γ there is an ε > 0 such that Γ ∩ (τε, τ + ε) = {τ }.

Note that we shift time over θ to be able to use the definition of an initial solution, which is given for an initial condition at t = 0.

Note that the definition of GAS has to be slightly generalized to allow for exogenous signals.


Energy dissipation and dispersion effects in granular media

Zhen Zhao and Caishan Liu * State Key Laboratory for Turbulence and Complex Systems, College of Engineering, Peking University, Beijing, China 100871 Bernard Brogliato † INRIA, Bipop Team-Project, Inovallée Montbonnot, 655 avenue de l 'Europe,France The strong interactions between particles will make the energy within the granular materials propagate through the network of contacts and be partly dissipated. Establishing a model that can clearly classify the dissipation and dispersion effects is crucial for the understanding of the global behaviors in the granular materials. For particles with rate-independent material, the dissipation effects come from the local plastic deformation and can be constrained at the energy level by using energetic restitution coefficients. On the other hand, the dispersion effects should depend on the intrinsic nature of the interaction law between two particles. In terms of a bistiffness compliant contact model that obeys the energetical constraint defined by the energetic coefficients, our recent work related to the issue of multiple impacts indicates that the propagation of energy during collisions can be represented by a distributing law. In particular, this law shows that the dispersion effects are dominated by the relative contact stiffness and the relative potential energy stored at the contact points. In this paper, we will apply our theory to the investigation of the wave behavior in granular chain systems. The comparisons between our numerical results and the experimental ones by Falcon et al., ͓Eur. Phys. J. B 5, 111 ͑1998͔͒ for a column of beads colliding against a wall show very good agreement and confirm some conclusions proposed by Falcon et al. Other numerical results associated with the case of several particles impacting a chain, and the collisions between two so-called solitary waves in a Hertzian type chain are also presented.

I. INTRODUCTION

Granular materials suffering shocks exhibit complex behaviors ͓1-5͔. In particular, some ordered phenomena such as so-called "solitary wave" in one-dimensional ͑1D͒ chain systems ͓6-30͔, different pattern formation in vibrated layers ͓2,31-35͔, as well as the collective behaviors in elongated particles ͓20,36-40͔, could be observed if the microscopic properties of particles are appropriate. Yet understanding the relationship between microscopic sizes and macroscopic behaviors in granular matters still remains a challenge in mechanics and physics due to the dissipative and dispersive nature of forces acting on the interacting grains. Part of the energy within the whole system will be dissipated due to the local interplays at the contacts, such as the plastic deformation and friction as well as the adhesive forces induced by the complex interaction with an interstitial fluid ͑water or air͒ ͓2,12,13,31͔. Meanwhile most of it will propagate and disperse to make more particles participate into the impacts during a shock process ͓6,8,41͔. A better and deeper understanding of the mechanisms related to the dissipation and dispersion effects of the energy is fundamental for the studies of granular matters.

It is obvious that the local contact model plays a significant role for discovering the properties of the dissipation and dispersion of the energy in granular matters. Except for friction at tangential contacts, the dissipation effects during im-pacts come from the local plastic deformation or the viscoelasticity of the rate-dependent materials for the normal contacts between particles ͓31͔. In order to encompass the dissipation effects in molecular-dynamics simulations, the simplest way of addressing the specifics in the course of impacts is Newton's kinematic coefficient of restitution, which originates from the concept of a binary collision between two identical spherical particles. This method can dramatically simplify the numerical simulation for the hugedimension granular systems because the impact process is ignored by setting the ratio between velocities after and before the collision with a fixed value ͓39,41͔. However, the approximation often yields some nonphysical behaviors such as the inelastic collapse or the divergence of the number of collisions in a finite time due to the well-known knowledge that the value of the Newton's coefficient is not a constant. In order to consider the influence of the impact velocity on the local plastic deformation between particles, a modification could be found in recent papers ͓34,35͔ by setting the coefficient of restitution with a velocity-dependent value. Such velocity-dependent restitution coefficient models have recently shown to be important in numerical and experimental studies. However, this improvement cannot essentially reflect the dissipation mechanism during the impacts because the value of the kinematic coefficient is not only related to the local plastic deformation, but also to the couplings among various contacts.

A more direct way of modeling the effects of energy dissipation is the so-called soft-particle algorithm in moleculardynamics simulation, in which all the forces acting on a particle from either walls or other particles or external forces, are calculated based on the positions of the particles. Once the forces are found, time is advanced by the explicit integration of the corresponding Newton equations of motion. There exist various models for calculating normal contact forces: some of them based on the viscoelasticity ͓14,42͔, others on the elastoplasticity of the material properties of the contact particles ͓31͔. The shortcomings of the soft-particle method are on one hand that many parameters are required, consequently, the simulation is impractical, and that the ordinary differential equations ͑ODE͒-stiff problem may arise in time-scale integration. This may result in the divergence of the numerical simulation away from the realistic phenomena, and then mislead the numerical observation.

The energy dispersion is another characteristic of granular systems during impacts, and depends on the intrinsic nature of the interaction law between two particles ͓8͔. In a uniform 1D chain system with Hertz-type interaction law, the socalled solitary wave phenomena are widely discovered in plenty of experiments ͓13,16,20͔, in which the energy within the system will be confined in an approximate fixed space with about five-particle length. For the chain with a linear elastic law, however, the dispersive effects will be enlarged in comparison with the case for nonlinear Hertz-type contact, where the initial impulsive wave will move as a spreading wave with wavelength increase and its amplitude decrease during propagation ͓23,29͔. Moreover, the profile of the pulse generated in a granular chain will be significantly affected by the mass and the material properties of the impacting particle. For instance, the so-called solitary wave may be disintegrated into a train of waves or oscillations when ones change the impacting particle ͓11,12,14,19,20͔. The reason for that can be attributed to the phenomena of multiple compression and/or expansion cycles or the repeating impacts appearing at the contact between the impacting particle and the chain. In fact, the relative properties among various contacts ͑such as the mass ratios and stiffness ratios between particles͒ should be responsible for the potential structures of the impulsive behavior in granular materials ͓6͔. Therefore, an underlying physical law depending on the interaction law should be established for the energy dispersion in the granular systems, such that the reason for the formation of some ordered phenomena in granular media could be discovered.

In ͓6͔, we proposed a method for solving frictionless multiple impacts in multibody systems, in which the interaction law between particles takes the power-type form for the relationship between the acting force and the indentation. The Stronge energetic restitution coefficient ͓43,44͔ is used to describe the extent of energy dissipation for dry grains with rate-independent material ͑such as the stainless steel͒. In comparison with Newton's coefficient, which usually mixes the dissipation effects and the dispersion effects together, the energetic restitution coefficient defines the dissipation in an energy level by using a ratio between input and output energy for a single compression and/or expansion cycle. This coefficient assumes that the extent of dissipation is independent on the dynamical process of multiple impacts, thus separates the energy dissipation from the energy dispersion. The dissipated energy at a contact is obtained according to the energetical coefficient and the input energy accumulated in the compression phase, which depends on the interaction law between the particles. The dissipated energy for the con-tact with multiple compression and/or expansion cycles at the same contact point can also be obtained by using the energetic coefficient, in which the dissipated energy in each loading and/or unloading cycle is separately calculated according to the input energy of the cycle.

The dispersion effects for particles with power-type interaction laws are also well discovered in ͓6͔ by transferring the time scale into the normal impulsive scale to describe a shock process, such as the way Darboux-Keller's method is used for the frictional collision with a single contact ͓5͔. In particular, the theoretical analysis indicates that the evolution of the energy is dominated by a distributing law that is closely related to the relative stiffnesses and the relative potential energies stored at various contacts. This law presents a clear scenario for the dispersion effects in a shock process, and may benefit us with the understanding of the complex behaviors in granular systems. Moreover, advantages in the numerical integration could also be brought in using this model because the description of the impacts is established at the energy level. Unlike the method of so-called soft-particle algorithm where the calculation for the contact forces is needed prior to the integration of the second-order ODEs, the method proposed in ͓6͔ uses the normal impulse as an independent variable, thus avoids the problems of stiffness encountered with second-order nonlinear ODEs by various numerical schemes ͓45͔.

In this paper the theory proposed in ͓6͔ is applied to investigate the shock behavior in granular chain systems. The friction effects between particles are excluded from this paper and will be studied in the future. The effects of energy dispersion and dissipation in granular media will be emphasized. The organization is as follows: Sec. II will introduce a brief overview of the theoretical aspects. In Sec. III, we will provide a comparison between the obtained numerical results and the experiments presented in ͓8͔, in which the influence of the property of the wall on the collision, the distribution of the postimpact velocities of the beads after the collision, as well as the contact force felt at the wall will be investigated. The numerical results for the case of a Hertz-type granular chain with the collisions of two solitary waves and the one of a granular chain impacted by n particles will be presented in Sec. IV. Summaries and conclusions are given in the final section.

II. BRIEF INTRODUCTION OF THE MULTI-IMPACT DYNAMICS

A. The impulsive dynamics For a granular system composed of N spherical particles with s frictionless contacts, the equations of motion take the following form if no rotational degrees of freedom is considered:

M͑q͒q ¨-h͑q,q ˙,t͒ -͚ i=1 s w i ͑q͒ i = 0, ͑1͒
where q is the coordinate set of the mass centers of the particles, M͑•͒ is a symmetric mass matrix, h͑•, • , •͒ contains the inertial and applied forces. For each contact point i, i is the scalar value of the normal contact forces along the common normal to the surfaces of the contacting balls. The connection between normal contact forces and the generalized forces is defined by w i ͑q͒, which is related to the Jacobian matrices of the contact points ͓46͔.

The kinematic state of a contact is determined by the distance ␦ i ͑q͒ between the contacting bodies. Clearly ␦ i ͑q͒ represents the normal elastic deformation when contact is established. The relative velocity of the ith contact point is expressed as

␦ ˙i͑q͒ = w i T ͑q͒q ˙. ͑2͒
For a 1D granular chain consisting of n aligned balls with respective radius r i , ␦ i ͑q͒ can be expressed as

␦ i = ͑x i + r i ͒ -͑x i+1 -r i+1 ͒ = q i -q i+1 + ͑r i + r i+1 ͒, ͑3͒
where x i is the displacement of the mass center for the ball B i , ␦ ˙i Ͼ 0 corresponds to compression phases, ␦ ˙i Ͻ 0 corresponds to expansion phases. Let ͓t 0 , t f ͔ denote the time interval of the impact, which can be further divided into much smaller intervals ͓t i , t i+1 ͔. According to the Darboux-Keller model ͓5͔ ͑Sec. 4.2.5͒, an integration over ͓t i , t i+1 ͔ has to be done in order to achieve a representation of the equations of motion at the impulse level. Thus a set of differential equations with respect to the normal impulses can be obtained,

͵ t i t i+1 ͓Mq ¨-h -W͔dt = M͓q ˙͑t i+1 ͒ -q ˙͑t i ͔͒ -W͓P͑t i+1 ͒ -P͑t i+1 ͔͒ = M • dq ˙-W • dP = 0. ͑4͒
The terms M, W remain unchanged during the integration thanks to the assumption of constant configuration q on ͓t 0 , t f ͔. The vector h consists of finite, nonimpulsive terms and therefore vanishes by the integration ͑see ͓5͔, Chap. 1͒. The quantities dq ˙and dP are the changes of generalized velocities and normal impulses during ͓t i , t i+1 ͔, respectively. It is obvious that the normal impulses dP are not independent since they are changed during the same interval of time. In order to make the impact differential equations of multiple impacts solvable, a connection among the normal impulses that contains the wave effects should be discovered. Therefore, adding one way or the other some compliance effects at the contact points is necessary when dealing with multiple impacts.

B. The bistiffness compliant contact model

Establishing a realistic model for the particles in contact with each other is still a tough task. Usually, one often uses a simple relationship to relate the interaction force and the overlap between two particles, in which a contact process is characterized by compression and expansion phases ͓43,47͔. The compression ͑or loading͒ phase corresponds to a process of transferring the external kinetic energy into the potential energy stored at the contact point, while the expansion ͑un-loading͒ phase will release the potential energy to make particles generate new movements. It is obvious that there may exist a certain dissipation of energy during a compression and/or expansion cycle. In order to consider the dissipation effects due to the plastic deformation, Lankarani and Nikravesh ͓44͔ proposed a bistiffness compliant contact model, in which different force and/or indentation relationships are adopted for the compression and expansion phases.

Let us consider the case of the interaction between two particles with a power-type law, i.e., the force and/or indentation relationship for the compression phase at the contact point j can be expressed as

c,j = K j ͑␦ c,j ͒ , ͑5͒
where K j is the contact stiffness, and c,j is the contact force during the compression phase. The subscripts c and e represent quantities related to the compression and expansion phases, respectively. The exponent determine the kind of contacts between bodies ͑ = 3 2 is for Hertz contact, =1 is linear elasticity͒.

Let us assume the form of the interaction law is not changed even though some plastic deformation is generated in the compression phase. The force and/or indentation relationship for the expansion phase is

e,j = m,j ͩ ␦ e,j -␦ r,j ␦ m,j -␦ r,j ͪ , ͑6͒
where ␦ r,j is the plastic deformation, and m,j and ␦ m,j correspond to the maxima of the normal contact force and normal deformation at the end of the compression phase, which corresponds to the values when ␦ ˙j = 0. For simplicity we omit the subscript j in the following expressions.

In order to constrain the local energy loss due to the permanent plastic deformation ͑the undeterminate scalar ␦ r ͒ generated in the compression phase, the energetic restitution coefficient ͓43͔ is applied on the compliant contact model. The integration of expressions ͑5͒ and ͑6͒ leads to

W c = ͵ 0 ␦ m c ͑␦ c ͒d␦ c = 1 + 1 K j ͑␦ m ͒ +1 , W r = ͵ ␦ m ␦ r
e ͑␦ e ͒d␦ e = -1

+ 1 K j ͑␦ m ͒ ͑␦ m -␦ r ͒.

͑7͒

Based on the definition of the energetic coefficient, W c and W r can be connected by e s 2 =-W r / W c , where W c ജ 0 and W r ഛ 0. In the case of a constant e s , we have from Eq. ͑7͒

␦ r = ␦ m ͑1 -e s 2 ͒. ͑8͒
The indentation ␦ m corresponds to the first instant t c such that ␦ ˙͑t c ͒ = 0, so it is not a parameter of the impact dynamics. From Eq. ͑8͒ neither is ␦ r . When the energetic coefficient takes a constant value, the energetic definition defines a linear relationship between the local plastic deformation ␦ r and the maximum indentation ␦ m , such that the dissipation of energy is reflected in an average level.

C. The potential energy at a contact point

Since the potential energy at a contact point will be accumulated during the compression phase and released in the expansion phase, we must calculate the potential energy in the two different processes. Let us pick a point p in the curve of the compression phase. The potential energy E p at p is equal to the work done by the contact force along the path from start of the contact to the position p. At any instant during the compression process, the contact force c ͑͒ can always be expressed as c ͑͒ = dP͑͒ d , in which P͑͒ is the normal impulse accumulated in the time interval ͓0,͔ corresponding to a deformation ␦͑͒. Thus, we have

E p ͑␦ p ͒ = ͵ 0 ␦ p dP͑͒ d␦͑͒ d . ͑9͒
Since P͑͒ and ␦͑͒ can be connected by a one-to-one mapping during the compression phases, we can use the variable P͑͒ to replace ␦͑͒ as the integral variable. Thus,

E p "P͑t͒… = ͵ 0 P͑t͒ ␦ j "P͑͒…dP͑͒ ͑ 10͒
where P͑t͒ is the normal impulse that is needed to make the indentation change from zero to ␦ p by obeying the relationship ͑5͒.

When the compression process finishes, the accumulation of energy at the contact point j will be ended, and the total accumulated potential energy can be expressed as

W c ͑␦ m ͒ = ͵ 0 ␦ m c ͑␦͒d␦. ͑11͒
After that the potential energy will be released through an expansion phase, in which the total energy that can be recovered through the expansion phase will take the value of e s 2 W c based on the energetical constraint.

Similarly to the compression phase, the residual potential energy for a dynamical process reaching a position R located at the expansion curve can be expressed as

E R = ͵ 0 P c ␦ ˙"P͑͒…dP͑͒ + 1 e s 2 ͵ P c P t ␦ ˙"P͑͒…dP͑͒, ͑12͒
where P c is the normal impulse when the compression phase finishes ͑so ␦ = ␦ m ͒, P͑͒ is an integral variable that is related to the normal impulse of the contact force experiencing a time interval ͓0,͔ by obeying the energetic relationship defined by the compliant contact model. The proof for establishing the expression ͑12͒ can be found in ͓6͔.

The first term in Eq. ͑12͒ corresponds to the total potential accumulated in the compression phase, while the second one means that the potential energy, 1 e s 2 ͐ P c P t ␦ ˙(P͑͒)dP͑͒, has been released at position R. For a full compression and/or expansion cycle, there is e s 2 W c potential energy to be recovered in the expansion phase. If a contact experiences a multiple loading and/or unloading cycles ͑a new compressional phase may restart before the expansion phase finishes͒, we also show in ͓6͔ that the transformation between the kinetic and potential will also obey the similar relationship expressed in Eq. ͑12͒. Each cycle will dissipate a part of the energy that is related to the restitution coefficient e s and the potential energy at the end of each compressional phase.

D. The distributing rule for the bistiffness compliant contact model

In terms of the compliant model expressed in Eq. ͑6͒ and the scale transformation of converting the time scale into the impulse scale, it is shown in ͓6͔ that the normal contact force at the contact point j during the compression and expansion phases can always be expressed as c,j "P j ͑t͒… = ͑1 + ͒ /͑+1͒ K j 1/͑+1͒ ͓E"P j ͑t͒…͔ /͑+1͒ . ͑13͒

Since the contact forces can always be expressed as the differential of the normal impulse with respect to time, from Eq. ͑13͒ the ratios between the increments of the normal impulses among various contact points j can be expressed as

dP i dP j = ␥ ij 1/͑+1͒ "E ij ͑P i , P j ͒… ͑+1͒/ , i = 1,2, ... ,s,i j,

͑14͒

when all the contact points have the same coefficient, and where ␥ ij = K i / K j and E ij = E i / E j are related to the ratios of the contact stiffnesses and the potential energy between contact points i and j. Equation ͑14͒ reflects the distribution of the increments of normal impulses in space since dP i and dP j vary within the same time interval. It is clear that the mode of the potential energy accumulation ͑the exponent in the force and/or indentation relationship͒ and the relative quantities of the contact stiffness among various contacts will dominate the transformation of energy within the network of contacts.

The outcomes of the multiple impacts can be obtained by combining Eq. ͑14͒ with the impulsive differential equations ͑4͒ and selecting a normal impulse at a contact point as an independent "timelike" variable. A guideline is also proposed in ͓6͔ for the selection of the independent impulse, in which the corresponding contact point is termed as the "primary contact point." This point should correspond to the point that takes the maximum potential energy among various contacts.

Since the impulse scale is used as the independent variable, the absolute values of the contact stiffnesses are not needed for the calculation of the postimpact velocities of each particle. However, if the value of the contact stiffness K j can be estimated effectively, the value of the contact force j at the instant of P j ͑t͒ can be obtained according to the potential energy E j (P j ͑t͒), as shown in Eq. ͑13͒. Furthermore, the time t related to the normal impulse P j ͑t͒ can also be extracted by using t = P j ͑t͒ / j . This is a postprocessing associated with the impulsive differential equations, i.e., we store the information of the potential energy at each numerical step for the ODEs, then carry out an algebraic operation to obtain the values of the contact force and the corresponding time.

The multi-impact Darboux-Keller dynamics is summarized as follows:

͑1͒ Contact parameters: ␥ ij , e j , 1ഛ i ഛ s, 1ഛ j ഛ s, ͑ = 1-linear stiffness-or= 3 2 -Hertz contact-, or other suitable values͒, E 0,j for the initial precompression cases.

͑2͒ Dynamical equation:

M dq ḋP i = W dP dP i if E ji ͑P j , P i ͒ ഛ 1 for all j i ͑15͒ with dP j dP i = ␥ ji 1/͑+1͒ "E ji ͑P j , P i ͒… ͑+1͒/ , ͑16͒ E ji = E j ͑P j ͒ E i ͑P i ͒ , 1 ഛ i ഛ s, 1 ഛ j ഛ s, ͑17͒ E j ͑P j ͒ = E 0,j + 1 Tra ͵ 0 P j ͑t͒ w j T q ˙dP j . ͑18͒
Tra is a parameter to transfer the work done by the normal impulse into the potential energy, in which Tra = 1 if ␦ ˙Ͼ 0, and Tra = e j 2 if ␦ ˙Ͻ 0. When E j (P j ͑t f ͒) = 0 and ␦ ˙Ͻ 0, the con- tact j will be open at time P j ͑t f ͒. The details of the numerical method and its implementation are in ͓6͔.

III. THE COMPARISON BETWEEN NUMERICAL AND EXPERIMENTAL RESULTS

In this section the experimental results reported in ͓8͔ are compared to the numerical results obtained by the approach outlined in Sec. II. In ͓6͔ numerous comparisons between compliant models ͑Hertz or linear springs͒ and the Darboux-Keller model of Sec. III are presented. They already prove that the Darboux-Keller model possesses very interesting properties. The results of this section confirm these positive conclusions.

A. A column of beads colliding against a wall

Figure 1 shows the granular system studied in ͓8͔, in which N identical stainless steel beads, each one 8 mm in diameter, are put together to form a column colliding against a fixed wall. The column of beads leaves the wall with a height h and starts its free fall to collide against the wall. The number of beads may vary from N = 1 up to N = 40 for experiments, and the contact forces felt at the wall are recorded by a force sensor. In order to reproduce the experimental results reported in ͓8͔, we use the bistiffness relationships expressed in Eqs. ͑5͒ and ͑6͒ to model the interactions at contacts, in which the local energy loss is taken into account by the restitution coefficient e s,s for particle to particle and e s,p for particle to wall, respectively. Let us number the contact bodies in a sequence from the wall assigned with 0 to the last bead with N. The contact points are denoted in a sequence with number 1 for the contact between the wall ͑0th contact body͒ and the 1th bead, and number N for the contact between ͑N -1͒th and Nth beads. According to Hertz theory, the exponent equals 3 / 2 for the nonconforming contact mode between bodies, and the contact stiffness K i at the contact point i is expressed as

K i = 4 3 R i * 1/2 E i * , i = 1, ... ,N, ͑19͒
where

R i * = R i R i-1 R i +R i-1 and 1 E i * = 1-i-1 2 E i-1 + 1-i 2
E i are the equivalent contact radius and the equivalent Young modulus between adjacent contact bodies, respectively. The mass of a stainless steel bead used in experiment is m = 2.05ϫ 10 -3 kg. The Young modulus and Poisson ratio for stainless steel are E s = 21.6ϫ 10 10 N / m 2 and s = 0.276, respectively. Thus, the value of the contact stiffness, K i ͑i =2, ... ,N͒ for spheresphere contact is 6.9716ϫ 10 9 N / m 3/2 . For the contact between the bead and the wall made of stainless steel, the value of the contact stiffness K 1 for the sphere-plane contact is 9.858ϫ 10 9 N / m 3/2 .

B. The contact force felt at the wall during the collision

Based on the description in ͓8͔ for the experiments, two cases are investigated in the following numerical simulations: a column of N beads free fall to collide against the wall with two different height h = 3.1 mm and h = 5.1 mm, corresponding to the impacting velocities of the column ͑v imp = ͱ 2gh͒ 0.246 m / s and 0.316 m / s, respectively.

For the columns of beads with N changing from 1 to 8, colliding against the wall with h = 3.1 mm, Fig. 2 shows the outcomes of the contact forces at the wall. It is clear that our numerical results can well reproduce the experimental phenomena discovered in ͓8͔, Fig. 4a and4b, in which the maximum force F max is approximately equal to 52 N and very little influenced by N, while the duration of impact increases linearly with N. Moreover, the contact force oscillates with a period of P Ϸ 32 s ͑the experimental value in ͓8͔ is P = 32.4Ϯ 1 s͒, and is damped similar to the response of the motion describing the damped free vibrations of a system with viscous damping.

Let us set the fall height h = 5.1 mm for the column of beads with various number of N =5,6, ... ,12, colliding against the wall. The numerical results obtained from our method are shown in Fig. 3, in which the maximum force F max is approximately equal to 71 N and the period P for the oscillation of the contact force is around 31 s. This also well agrees with the experimental results in ͓8͔, Fig. 4c and4d. Since the value of the restitution coefficient is often estimated from an independent experiment, it may be interesting to investigate its influence on the contact force at the wall by slightly modifying e s . Let us set both the values for particleparticle coefficient of restitution and the particle-wall one as e s = 0.92. Figures 4 and5 present the numerical results for the contact forces felt at the wall that is collided by the column of beads with two different fall heights h = 3.1 and 5.1 mm, respectively. Comparing them with Figs. 2 and3, we find that the little change of the coefficient of restitution has no identifiable influence on the value and the shape of the curves of the contact forces at the wall. The proposed model and numerical algorithm has good robustness ͑insensitivity͒ properties.

C. Distribution of the postimpact velocities

As illustrated in ͓8͔, at the end of the collision the beads in the column will detach from each other due to the energy redistribution within the system. In this section we will apply our method to investigate the phenomena of the detachment effect found in ͓8,48͔ by presenting the postimpact velocities of the beads after the collision.

The definition of the end of the column-wall collision given in ͓8͔, in which the collision is assumed to be finished at the time at which the bead at the bottom leaves the wall, is different from ours: we use the energy calibration to define the end of the column-wall collision, corresponding to the instant when all the potential energy among various contacts completely released. Figure 6 shows the postimpact velocities of the beads obtained from the simulation for the column free falling from a height h = 5.1 mm. It is clear from this figure that the upper beads of the column go away with velocities greater than v imp = 0.316 m / s, while the lower beads have velocities smaller than v imp . Moreover, we can also find that the postimpact velocity of the last ball will increase with N, and reach the maximum value at N = 8, then decrease when N is further enlarged. Unlike a uniform chain impacted by an identical particle, where a regular wave could be formed and confined within five particles, the wave behavior in the column is complex and unregular similar to the phenomena discovered in the next section, in which n-ball chains are impacted by m-ball chains ͑hence N = m + n͒. However, we can still observe the following interesting phenomena: If N Ͻ 8, the beads will separate one after the other from the top of the column. For N Ͼ 8, some beads separate, while others are leaving in a cluster since the discrepancies between the postimpact velocities of the beads are very small. If N Ͼ 8, the local energy loss will make the postimpact velocity of the top-bead decrease. Due to the local energy loss, it is imaginable that the top bead will not bounce anymore if the column is long enough.

The effective coefficient of restitution for the whole chain, proposed in ͓48͔ and then applied in ͓8͔, is defined as

e ef f = ͚ i=1 N v i + ͚ i=1 N v i -= - 1 Nv imp ͚ i=1 N v i + , ͑20͒
where v i + corresponds to the postimpact velocity of ith beads. It is noteworthy that the value of v i + is different from the one of v i f in ͓8͔, Eq. ͑7͒, where all the values for each v i f correspond to the instant when the force felt at the wall equals zero. Figure 7 presents the value of e ef f varying with N for a column free fall with a height h = 5.1 mm. The tendency of the variation of e ef f agrees well with ͓8͔, Fig. 11. However, when N Ͼ 4 there is an obvious discrepancy between the value of e ef f obtained from our numerical results and the one shown in ͓8͔, Fig. 11. can be attributed to the different for v i + in the two figures for the calculation of e ef f . Comparing Figs. 2 and4, with Figs. 3 5, respectively, we can find that the variation of the value of the restitution coefficient e s,s has little influence on the value of the contact force at the wall. However, it will change the of the postimpact velocities of beads ͑the dispersion effect͒ that seems not to be ignorable. Let only change the value of e s,s from 0.96 to 0.92 while keeping other parameters with the same values as used for Fig. 7. We can find Fig. 8 that the curve of e ef f has an identifiable variation due to the little change of e s,s . This means that the local energy loss will influence the distribution of the postimpact velocities of the beads. This can confirmed Fig. 9, e s,s is assigned as 0.7, and the effective coefficient of the whole column, e ef f for N = 20 is near to corresponding to the column with a

D. Influence of the mechanical properties of the wall

The mechanical property of the wall will significantly influence the collision processes as exhibited in ͓8͔. Further investigations can also be found in recent studies in the physics community ͓12,15,16͔. For instance, Job et al. ͓15͔ performed a precise experiment to reflect the dependence of the impulsive wave on the wall mechanical properties. Similar phenomena are also found in the experiments by Daraio and co-workers ͓12,13͔.

In ͓8͔, Falcon et al. carried out different experiments by setting the column of beads colliding against the wall with different materials. According to the output signal from the force sensor they found that the rigidity of the wall will influence the collision behavior. For example, experimental results show that the maximum contact force felt at the wall will depend on the number N when a soft material sheet is stuck on the sensor. Since the mechanical parameters for the materials covered in the wall are not provided in ͓8͔, we use the duration for one bead colliding with a plane to estimate the contact stiffness between the bead and the plane covered with various material. According to Hertz theory, the contact time 1 between a bead dropped from a height h and a plane assumed to be of infinite mass can be approximately expressed as

1 = 2.94ͩ 5m 4K 1 ͪ 2/5 v imp -1/5 . ͑21͒
For a bead with mass m = 2.05ϫ 10 -3 kg falling from a height h = 2.9 mm ͑v imp = ͱ 2gh = 0.2385 m / s͒, Table I pre- sents the values of K 1 extracted from the max ͑N =1͒ in ͓8͔, Table 2 and expression ͑21͒.

We can find from Table I that the value of K 1 extracted from the experiment of the stainless steel bead colliding with a stainless steel plane is larger than the theoretical value ͑K 1 = 9.858ϫ 10 9 N / m 3/2 ͒, so it is not possible to provide a precise comparison for the contact force between our numerical results and the experimental ones in ͓8͔, Fig. 13. At best one can compare some qualitative features of the collision due to the variation of the material in the sheet.

Let us first investigate the contact force at the wall covered with polyvinyl chloride ͑PVC͒ ͑index 2 in Table I͒ by setting a column of beads with various N free falling from a height h = 2.9 mm. From Fig. 10, the experimental phenom-

N = 1 N = 2 N = 3 N = 4 N = 5 N = 6 N = 7 N = 8 N = 9 N = 10 N = 14 N = 11 N = 12 N = 13 N = 15 N = 16
FIG. 10. The contact force at the wall with PVC during the collision of a column with N =1,2, . . . ,16, 2.9 mm, e s,s = 0.92, e s,p = 0.7.

N = 1 N = 2 N = 3 N = 4 N = 5 N = 6 N = 7 N = 8 N = 9 N = 10 N = 11 = 12 N = 13 N = 14 N = 15 N = FIG.
11. The contact force at the wall with beech during the collision of a column of N =1,2, . . . ,16, h = 2.9 mm, e s,s = 0.92, e s,p = 0.92. ena discovered in ͓8͔ are well reproduced in our numerical results, in which the maximum force increases with N at low values of N, until it becomes independent of N, finally the contact force will slightly decrease with N similar to a damping response.

If the wall is covered with softer material, it is imaginable that the number N related to the maximum independent force will increase, and the amplitude among all the maximum forces will decrease due to the small value of the contact stiffness between the bead and the wall. Figures 11 and12 show the results obtained from our numerical method for the column of beads with various N-chains colliding against a wall covered respectively with beech and rubber materials. The qualitative characteristics discovered in ͓8͔ are well reproduced even though the values of the contact stiffness K 1 are roughly estimated with possible large errors compared with practical situations.

In summary, the local energy loss can be easily taken into account by using energetical constraints and the information related to the contact force and the duration of impact could also be extracted from the solutions of the differential impulsive equations through simple algebraic operations. The comparisons between our numerical results and the experimental ones in ͓8͔ show good agreement and confirm some conclusions proposed therein.

IV. OTHER NUMERICAL RESULTS RELATED TO GRANULAR CHAINS

A. n-ball chains impacted by m-ball chains

As proved theoretically and experimentally by many authors, an impulsive wave can be generated when a chain endures an impulse ͓7,15,16,23,25,29͔. In particular, a com- pressive wave with regular profile involving five particles can formed when a uniform ball-chain system is impacted by an identical particle. However, this regularity of the wave behavior will be destroyed when the chain is collided with by more than one impacting particles.

Let us denote a system with 10 balls that are identical and satisfy the Hertz contact relationship. The parameter ␤ is defined as the ratio between the numbers of the m impacting balls and the n impacted balls. For instance, ␤ =1:9 means that the chain with nine stationary balls is impacted by one particle ͑i.e., n = 9 and m =1͒, ␤ =2:8 is related to situation that two particles collide against eight stationary balls ͑i.e., n = 8 and m =2͒, etc. We will investigate the impulsive behaviors by changing ␤ from 1:9 to 4:6. The initial velocities of the impacting balls are set as v 0 =1 m/ s and no energy is dissipated during impacts: e s,i = 1 for all contacts i.

Figure 13 presents the evolutions of the velocities and the normal impulses for the case of the chain with parameter ␤ =1:9. Clearly a regular wave can be formed and only three balls move forward when the impact finishes, that is near to the half length of the wavelength. The values related to the final velocities for each ball are presented in Table II, in which the last row shows the percentages of the forward kinetic energy ͑the sum of the kinetic energies of the balls with positive post-impact velocity͒ and of the energy of the final ball with respect to the total energy of the system, respectively. In this case, the percentage of the forward kinetic energy is about 99%, and the one of the final ball is near to 97.4%.

When the chain is impacted by two identical balls, there are four balls finally moving forward, and the kinetic energy will be centralized in the last two balls. It is obvious that the wave behavior is not regular due to the complex coupling among various contacts. Figures 14 and15 show the evolutions of the postvelocities of balls and the normal impulses for balls 1 to 5 and for balls 5 to 10, respectively. Although most of the kinetic energy of the system is still kept in the balls moving forward, the distribution of the postimpact velocities of the balls is much different from the one in the chain impacted by one identical ball. The forward kinetic energy is not concentrated in the last two balls, but will disperse in a length with about four balls.

As m increases, the interactions between the contacts becomes more complex, and the number of balls rebounding backwards decreases. Figures 16 and 17 present the evolutions of the velocities and of the normal impulses for the balls from Nos. 1 to 5 and from Nos. 6 to 10 in the chain with ␤ =3:7, respectively. Observing the values presented in Table II we can find that the forward kinetic energy changes little, while the number of balls moving forward is five and the kinetic energy is dispersed in a wider scope.

When m increases to four, the postimpact velocity of the last ball changes very little in comparison with the cases of a chain impacted by two balls or three balls, while the forward kinetic energy will disperse in the chain with six balls as shown in Table II. The evolutions of the velocities of each ball and the impulses among various contacts are shown in Figs. 18 and19.

In conclusion, as m increases, the ratio ͑forward postimpact kinetic energy͒ ͑total kinetic energy͒ Ϸ 99 %, and the number of balls with a positive postimpact velocity increases as m + 2. These results confirm the proposed scheme is able to reproduce the dispersive effects observed for such impacts ͓8,48͔.

B. Two solitary waves collision

Collisions between two identical solitary waves traveling in opposite directions are also investigated in the physics community ͓17,49͔. A question to be concerned with herein is whether at the point of intersection, the opposing solitary waves could be annihilated completely, i.e., whether the center of the central grain suffers any motion at any time. Let us set a chain with 15 identical balls that are stationary and keep contact with adjacent balls. The contact relationship between the balls is assumed to be nonlinear elastic and satisfies Hertz assumption ͑i.e., =3/ 2͒. Suppose that two balls with the same velocity, ͉v 0 ͉ =1 m/ s collide against the ends of this chain from two opposite directions. Two impulsive waves will be generated at the ends of this symmetric structure, and are propagated through the chain ͑shown in Fig. 20͒. The five curves correspond to the velocities of the balls at five different instants ͑different P j ͒, where the solid line represents the initial condition of the chain, while the other four lines describe the impulsive wave behavior at different times.

Initially the imparting balls at the ends of this chain only influence the adjacent balls with limited length, and then two impulsive waves are formed due to the coupling between adjacent balls ͑self-organized behaviors among the balls͒, then travel through the chain with an approximately fixed length and a constant amplitude ͑similar to a solitary wave͒. Meanwhile, the balls behind the wave will rebound backwards with small velocities, confirming the scattering phenomenon or bead detachment effect ͓8͔. The length of the impulsive wave is confined in five balls, that well coincides with conclusions obtained from numerical results, theoretical analysis as well as the experiments reported in ͓9,15,29͔. The two impulsive waves originated from the ends of this chain will meet at the chain center, in which the located particle always keeps motionless due to the symmetric structure of the chain. After that, the waves with unchanged shape will reflect from the center of the chain to make particles separate from each other.

If the number of balls in the chain is even, the two impulsive waves will interfere with each other when the waves reach the center of the chain. This phenomenon is well reflected in Fig. 21 obtained by setting the chain with 16 identical balls. When the two waves met at the center of the chain, the velocities of the particles involved in the impulsive wave will abruptly decrease, and thus the momenta resided in particles will be redistributed to form a different shape of wave profiles. However, the variation of the wave profile has little influence on the final outcomes of the velocities after impacts according to our numerical simulation. It is imaginable that the waves will be split or oscillate when the center particle is different from the other ones, and socalled secondary solitary waves will be generated after the two solitary waves meet with one another. These phenomena have been pointed out in ͓17,49͔.

V. SUMMARY AND CONCLUSION

The energy of colliding bodies during granular systems will be transferred and redistributed through the network of contact points. Moreover, part of the energy will be dissipated at the vicinity of contacts because of the plastic deformation or other factors. Establishing a model that can clearly separate the dissipation effects and the dispersion of energy among various contacts is crucial for the understanding of the rich and complex behaviors exhibited in granular media.

According to the method proposed in ͓6͔, the energetic coefficient seems to be a good candidate in the energy sense to confine the local energy loss, while the dispersion effects should depend on a distributing law that is related to the relative stiffness and the mode of the potential energy accumulation between two particles. This method clearly separates dissipation and dispersion effects during multiple impacts and gives a legible for the dynamics of multiple impacts. The comparisons between our numerical predictions and the experimental results in ͓8͔ show good agreement, thus validating the proposed algorithm. This method may highlight a way of better understanding the complex phenomena exhibited in granular systems and paves the way to controlling the global behavior through adjusting the microproperties of the particles.

  FIG. 1. A column of beads colliding against a wall.

FIG. 2 .FIG. 3 .

 23 FIG. 2.The contact force at the wall during the collision of a column of N =1,2, . . . ,8 beads with fall height h = 3.1 mm, e s,s = 0.96, e s,p = 0.92.

FIG. 5 .

 5 FIG.4. The contact force at the wall during the collision of column of N =1,2, . . . ,8 beads with fall height h = 3.1 mm, e s,s = 0.92, e s,p = 0.92.

  FIG.6. ͑Color online͒ Velocity of the ith bead, at the end of the collision with fall height h = 5.1 mm, as a function of the bead's number i, for different values of N =5,6, . . . ,12 beads, and e s,s = 0.96, e s,p = 0.92.

  FIG. 9. The coefficient e ef f for the column with different N dropped from a height h = 5.1 mm, and setting e s,s = 0.7 and e s,p = 0.92.

  the bottom ball and the wall (N)

-

  FIG.12. The contact force at the wall with rubber during the collision of a column of N =1,2, . . . ,16, h = 2.9 mm, s,s = 0.92, e s,p = 0.92

  FIG.14. The evolutions of the velocities and the normal impulses for chain with eight particles impacted by two balls: ␤ =2:8 ͑balls 1 to 5͒.

  FIG.15. The evolutions of the velocities and the normal impulses a chain with eight particles impacted by two balls: ␤ =2:8 ͑balls 6 to 10͒.

  FIG. 17. The evolutions of the velocities and the normal impulses for chain with seven particles impacted by 3 ␤ =3:7 ͑balls 6 to 10͒.

  FIG. 19. The evolutions of the velocities and the normal impulses for a chain with six particles impacted by 4 balls: ␤ =4:6 ͑ball 6 to 10͒.

ACKNOWLEDGMENT

  This work was supported by the National Science Foundation of China ͑Grant No. 10772002͒. ͓1͔ C. Goldenberg, and I. Goldhirsch, Nature ͑London͒ 435, 188 ͑2005͒. ͓2͔ I. S.

TABLE I .

 I The value of the contact stiffness K 1 estimated from experimental results for the contact time between a bead dropped from a height h and a plane assumed of infinite mass in ͓8͔ Table2.

	Index No.	Material	max ͑s͒	K 1 ͑N / m 3/2 ͒
	1	steel		1.15ϫ 10 10
	2	PVC	67.5	3.67ϫ 10
	3	beech	127.5	7.5ϫ 10 7
	4	rubber	221	1.89ϫ 10 7

TABLE II .

 II The final velocities for each ball with different ␤. Sign "ϩ" means "moves forward," while "Ϫ" means "moves backwards." The first figure in the last row represents the ratio between the forward postimpact kinetic energy and the total kinetic energy; the second figure is the percentage of kinetic energy in the last ball.

	␤	͑1:9͒	͑ 2:8͒	͑ 3:7͒	͑ 4:6͒
	q ˙1 +	-0.0710	-0.1126	-0.1441	-0.1706
	q ˙2 +	-0.0303	-0.0481	-0.0612	-0.0729
	q ˙3 +	-0.0159	-0.0248	-0.0312	-0.0373
	q ˙4 +	-0.0089	-0.0133	-0.0169	-0.0146
	q ˙5 +	-0.0052	-0.0068	-0.0054	0.1274
	q ˙6 +	-0.0030	-0.0022	0.0996	0.4648
	q ˙7 +	-0.0018	0.0497	0.4043	0.4847
	q ˙8 +	0.0025	0.2893	0.5108	0.9111
	q ˙9 +	0.1467	0.6570	1.0118	1.0928
	q ˙10 +	0.9869	1.2118	1.2323	1.2145
	Percentage ͑99.5, 97.4͒ ͑99.3, 73.4͒ ͑99.2, 50.6͒ ͑99.2, 36.9͒

  Aranson and Lev S. Tsimring, Rev. Mod. Phys. 78, 641 ͑2006͒. ͓3͔ M. V. Hecke, Science 317, 49 ͑2007͒. ͓4͔ V. Narayan, S. Ramaswamy, and N. Menon, Science 317, 105 ͑2007͒. ͓5͔ B. Brogliato, Nonsmooth Mechanics, 2nd ed. ͑Springer, London, 1999͒. ͓6͔ C. Liu, Z. Zhao, and B. Brogliato, INRIA Research report, 2008 ͑unpublished͒, http://hal.inria.fr ͓7͔ C. Coste, E. Falcon, and S. Five, Phys. Rev. E 56, 6104 ͑1997͒. ͓8͔ E. Falcon, C. Laroche, S. Five, and C. Coste, Eur. Phys. J. B 5, FIG. 21. ͑Color online͒ The impulsive wave shape in the 16-ball system at different P ti that represents the different instants scaled by impulse.

		0.6 0.8 1		Initial	P t1	P t2		P t3	P t4
	Velocities of the ball	-0.4 -0.2 0 0.2 0.4								P t5	Two waves merged
																		Initial
		-0.6
		-0.8									P t4	P t3		P t2	P t1
		-1
		0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17
									Numbers of the balls


Ling Rong Kong

Wei Hong

Huang

Pierre-Yves Hicher

Dong Mei Zhang

Microstructural analysis of Shanghai remolded clay during one-dimensional consolidation test

Keywords: microstructure, clay, consolidation test, SEM

The microstructure of soil has a significant influence on the mechanical behavior. In this paper, photos of Scanning Electron Microscope (SEM) are taken on Shanghai remolded clay samples under different consolidated pressures during onedimensional consolidation test. The evolutions of microstructure of remolded clay (including aggregate shape, aggregate mean surface and orientation distribution) from both horizontal cut plane and vertical cut plane of the specimen are analyzed and compared.

INTRODUCTION

Shanghai clay is characterized with high water content, high void ratio, low shear strength, high compressibility and low permeability. Because of the economical boom in Shanghai, the clay deposits in Shanghai become more and more solicited by different constructions. In general, soil properties are decided by its macroscopic behavior which can be represented by its stress-strain relationship, which is lacking of considering the soil microstructure. The microstructure is an important aspect of soil. The microstructure during one-dimensional consolidation has been studied by Martin and Ladd [1975], Mcconnachie [1974], Matsuo [1977], Hicher et al. [2000] and etc.

For better understanding the shanghai clay behavior, in this paper remolded Shanghai soft clay specimens are studied by standard oedometer tests. The SEM observation is made on samples loaded at different stress levels. The quantitative analysis of the microstructure of Shanghai remolded clay are studied.

SHANGHAI CLAY MICROSTRUCTURE PROPERTY

In this paper, the soil samples are from Nanxiang-Baili Unistate of Metro 11. The depth is 8∼8.5 m. Physical properties of the selected soil are listed in Table 1. In order to investigate the deformation mechanism, three conventional oedometer tests were carried out for samples under 132, 264 and 800 kPa.

Observations for both horizontal planes and vertical planes of shanghai remolded clay specimens under different consolidation pressures are shown in Fig. 2. The particles of remolded specimen are dominant on aggregates. The void between particles exists. The particle size with the diameter is distributed between 0∼10 µm, but the diameter of large formed aggregates is more than 20 µm. The form of aggregates is associated by face to face. Increasing the consolidation pressure, the particles and aggregates become dense for both horizontal and vertical planes. The aggregates are reduced into small-size aggregates when remolded. As the same time, there are some particles and small aggregates move to adjacent particles and aggregates to form new big aggregates. In addition, the particles and aggregates rearrange to be privilege on the direction under remolded force, where the orientation in vertical plane is more obvious than the horizontal plane.

QUANTITATIVE ANALYSIS OF SHANGHAI REMOLDED CLAY STRUCTURE EVOLUTION

The SEM photos can be used to quantitatively measure the evolution of clay structure elements, mainly of aggregates and single particles. The analysis is based on the following parameters:

1. Particle orientation F(a). F m M i F F ( ) % ) ) = × 100 (1)
where F i (α) means one direction distribution frequency, in each SEM photo, the plane was divided into a given number quadrants, each of them equal to 15°; m means the particle numbers in the ith quadrant; M means the total numbers in one SEM photo. In order to study particle orientation, the rose diagrams from SEM photos was carried out. 2. I s shape index, which is the value of the ratio between the long axis A and the short axis B of the individual particle or individual aggregate (the particle or aggregate is considered as an ellipse). If I s equal to 1, it means that the particle or aggregate is round or square. 3. S mean the mean surface:

S mean = πAB/2
(2)

Particle orientation analysis

For Shanghai remolded clay, shown in Fig. 3(a), we could see that the particles or aggregates distributed in every quadrant for both horizontal and vertical planes, but the particles or aggregates have two peak values: 0°∼15° and 75°∼90°. The particles or aggregates are composed of mainly frames for both in horizontal and vertical planes. The particle orientation is obviously changing with the increasing of consolidation pressure. For example, in Fig. 3(b) when the consolidation pressure is 264 kPa, there are more particles in 0°∼15° and 120°∼135°, but less particles in 75°∼90°. This means that the particles rotate when increasing the consolidation pressure, which results in particles existing in some privileged directions. The phenomenon is more remarkable for 800 kPa. It can be concluded that more anisotropy due to reorientation of particle under higher consolidation pressure.

Particle shape index analysis

The shape index histograms of aggregates for Shanghai remolded clay were shown in Fig. 4. The shape index dominant on 1.5∼2 for horizontal plane shows that the particle shape is mainly ellipse; the shape index dominant on 1.5∼2 and 2.5∼3 for vertical plane shows that the particle shape is more ellipse compared with horizontal plane. The particle shape is obviously ellipse for both horizontal and vertical planes with the increasing of consolidation pressures. The value of shape index is 2.5∼3. When the consolidation pressure is 800 kPa, the phenomenon is more remarkable. There are more ellipse shaped aggregates and less other shaped aggregates. During the one-dimensional compression, aggregates slap and change their shapes, which is in accordance with the anisotropy of the aggregates.

Particle mean surface analysis

The particle mean surface of aggregate for Shanghai remolded clay is shown in Fig. 5. The particle mean surface or equal diameter decreases with the increasing of consolidation pressure for both horizontal and vertical planes. We observed that big aggregates gradually break into two or more aggregates.

CONCLUSION

SEM photos show that the aggregate of Shanghai remolded clay is likely ellipse-sphere shape for both horizontal plane and vertical plane. Under the one-dimensional compression, the shape of aggregate changes more obviously for vertical plane rather than for horizontal plane. The one-dimensional compression results in a more orientated microstructure in the vertical plane rather than in the horizontal plane. The increase in direction frequency due to loading is mainly the results of aggregates breakage, with smaller aggregates easier orientated to the loading direction. Comparisons between SEM photos demonstrate that the Shanghai remolded clay has the characters of gradual damage, inherent and induced anisotropy during one-dimensional compression.

Figure 1 .

 1 Figure 1. Compressive curves up to different consolidation pressures for remolded samples.

Figure 2 .

 2 Figure 2. Remolded specimen SEM photos under different consolidation pressure (a-d for vertical planes and e-h for horizontal planes).

Figure 3 .

 3 Figure 3. Rose diagrams of aggregate orientation of Shanghai remolded clayunder different consolidation pressure.

Figure 4 .

 4 Figure 4. Shape index histograms of aggregate of Shanghai remolded clay under different consolidation pressure.

Figure 5 .

 5 Figure 5. Mean area of particles under different consolidation pressure.

Table 1 .

 1 Physical properties of soil.5

	depth/m	e 0			W/%	W L /%	W P /%	γ kg/m -3
	8∼8.5	1.12∼1.16		39.4∼43.7	51.6	25.8	1.71∼1.74
			1.2
							132 kPa
			1.1				264 kPa 800 kPa
		Void ratio e	0.8 0.9 1.0
			0.7
			0.6	1	10	100	1000
						Log σ' V