Stefano Bosi
email: stefano.bosi@univ-evry.fr.

Thomas Seegmuller

Carine Nourry

Alain Venditti

Optimal Cycles and Social Inequality: What Do We Learn from the Gini Index? *

Keywords: JEL classification: D31, E32, O41 Endogenous cycles, two-sector models, heterogeneous agents, Gini index

One of the plausible explanations for macroeconomic fluctuations relies on the occurrence of endogenous deterministic cycles. In the last three decades, most of the relevant literature has rested on the assumption of a representative agent but, recently, a few papers have investigated the role of consumers' heterogeneity on endogenous fluctuations. Our article aims at taking a step forward in order to give a more suitable interpretation. To keep things as simple as possible, we introduce heterogeneous households in a two-sector optimal growth model and we study how wealth heterogeneity affects the occurrence of endogenous cycles. In contrast to previous results, we relate the existence of such cycles to the most commonly used inequality measure, the Gini index, and analyze the impact of consumers' heterogeneity on this index.

Introduction

One of the puzzling questions the economists are addressed is whether income inequality can be pointed out as a catalyzer of macroeconomic instability, or, more explicitly, whether less developed countries with higher income inequality are less sheltered from economic fluctuations and suffer a wider macroeconomic volatility. In this paper, we aim at providing a theoretical answer, by considering 1 a usual neoclassical dynamic framework, where volatility is interpreted in terms of endogenous cycles.

Well-known features of the one-sector optimal growth models are the uniqueness of the steady state and the saddle-path stability. The monotonic convergence to the steady state is no longer ensured in multi-sector models. [START_REF] Benhabib | Competitive Equilibrium Cycles[END_REF] have shown that oscillations and two-period cycles can arise in a two-sector framework. In their model, the consumption and the investment good are produced through different technologies. Such a productive heterogeneity plays an essential role in the emergence of cycles, since the consumption good is required to be more capital intensive than the investment good. Nevertheless, most of this literature, which is concerned with the existence of competitive cycles and endogenous fluctuations, fails to introduce heterogeneity on the consumption side of the economy and takes the shortcut of a representative agent.

Recently, [START_REF] Herrendorf | Ruling Out Multiplicity and Indeterminacy: The Role of Heterogeneity[END_REF] have addressed some criticism to the occurrence of expectations-driven fluctuations under the assumption of a representative agent. Surprisingly, they prove that a degree of agents' heterogeneity can rule out indeterminacy. However, they deal with a quite peculiar model, where nominal rigidities matter. Thus, the robustness of their findings has been the subject of further studies since. A few authors have introduced heterogenous households in dynamic general equilibrium models. [START_REF] Ghiglino | Inequalities and Fluctuations in a Dynamic General Equilibrium Model[END_REF], [START_REF] Ghiglino | Wealth Inequality and Dynamic Stability[END_REF] and [START_REF] Ghiglino | Wealth Inequality and Macroeconomic Volatility in Two-Sector Economies[END_REF] have focused on the influence of heterogeneity on the occurrence of optimal cycles and non-monotonic trajectories. In these papers, households are heterogeneous, because they are represented by different utility functions and are unequally endowed with capital and labor. In order to show the effects of heterogeneity on the equilibrium dynamic properties, they exploit the market perfection and characterize directly the Pareto-optimal allocations without caring about the decentralized solution, by using the existence and continuity of the welfare weights with respect to initial conditions (see [START_REF] Kehoe | Determinacy of Equilibria in Dynamic Models with Finitely Many Consumers[END_REF] and [START_REF] Santos | Differentiability and Comparative Analysis in Discrete-Time Infinite-Horizon Optimization[END_REF]). [START_REF] Ghiglino | Wealth Inequality and Dynamic Stability[END_REF] and [START_REF] Ghiglino | Wealth Inequality and Macroeconomic Volatility in Two-Sector Economies[END_REF] find that the influence of inequality and heterogeneity on macroeconomic (in)stability depends on the concavity (convexity) of the absolute risk tolerance, that is the inverse of the absolute risk aversion. Despite the findings are rich and precise, they remain difficult to understand and to interpret in terms of heterogeneity at an individual level: these authors don't consider popular aggregate inequality measures. In addition, several original results obtained in these papers no longer apply to more consensual utility functions characterized by constant elasticities of intertemporal substitution in consumption.

Our paper takes into account these criticisms, by targeting and reexamining the connection between consumers' heterogeneity and stability properties in an optimal growth two-sector model, where consumers' preferences are characterized by a constant elasticity of intertemporal substitution. We will discuss not only the role of wealth inequality on dynamics, but also that of preferences heterogeneity. We connect the emergence of optimal cycles to the well-known in-equality measure, the Gini index, and analyze how it is affected by heterogeneity in preferences.

We focus directly on the decentralized equilibrium, instead of considering the Pareto allocations. The stationary solution is shown to depend on the sharing of initial endowments and, therefore, is formally characterized by the presence of a unit root. Such a feature is not a curiosity, but rather a general property of an heterogeneous world.

In line with [START_REF] Benhabib | Competitive Equilibrium Cycles[END_REF], we prove that endogenous cycles cannot arise when the investment good is relatively more capital intensive, while non-monotonic trajectories can occur if the pure consumption good is more capital intensive. In this case we are able to determine the critical distribution of capital stock beyond which two-period cycles and period-doubling flip bifurcations occur. Moreover, in contrast to the existing papers, we introduce a commonly used measure of inequality, the Gini index, the equilibrium value of which mainly depends on the wealth sharing. We find the critical Gini degree corresponding to the critical distribution of capital and we explain how it varies in response to changes in preferences heterogeneity, more precisely when heterogeneity in consumers' preferences raises. More heterogeneity in consumers' preferences raises or lowers the Gini index, according to the interplay between the substitution and income effects, the relative importance of which is determined by the consumers' average elasticity of intertemporal substitution.

In sum, optimal cycles occur under well-known technological conditions in our two-sector economy. However, in contrast to the existing literature, we provide clear-cut insights about the role of consumers' heterogeneity on the emergence of persistent cycles. Indeed, heterogeneity in preferences interferes with the technological source of cycles, through the unequal capital distribution. Eventually, we connect the appearance of cycles to the Gini index, one of the most commonly used measures of inequality. Hence, the role of heterogeneity on optimal cycles is revisited from a point of view perhaps much simpler than the existing approaches, but probably more understandable.

The rest of the paper is organized as follows. Section 2 and 3 are devoted to present, respectively, the behavior of heterogeneous consumers and producers. In section 4, we determine the equilibrium and analyze the steady state. In section 5, we define the measure of inequality we will use later: the Gini index. Section 6 presents the stability properties of the steady state and the conditions under which cycles occur. We also relate the flip bifurcation to the Gini index. In section 7, the impact of heterogeneous preferences on the critical Gini index is discussed. Concluding remarks are provided in the last section, whereas technicalities are gathered in the Appendix.

Heterogeneous consumers

We consider an economy with a constant population l of infinite-lived agents over a discrete time t = 0, 1, . . . , ∞. There are two types of consumers (i = 1, 2), who differ by their capital endowments and their preferences. More precisely, there is a continuum [0, l 1 ] of agents of type 1 (say, the poor), each one endowed with an amount k1 of capital, and a continuum (l 1 , l] of agents of type 2 (the rich), endowed with an amount k2 > k1 . Whatever the class he belongs, each household supplies inelastically one unit of labor at each period.

Two goods are produced and exchanged in this economy: a consumption good, which plays also the role of numéraire, and an investment good, which is sold at the price p t . The real interest rate r t and the real wage w t are measured in consumption good units, while ∆ ≡ 1δ ∈ (0, 1), where δ ∈ (0, 1) is the rate of capital depreciation. At each period, a consumer of type i faces a budget constraint:

cit + p t ³ kit+1 -∆ kit ´≤ r t kit + w t (1)
where cit represents his consumption at time t. His intertemporal preferences are summarized by an utility function:

∞ X t=0 β t c1-1/σi it 1 -1/σ i (2)
where β ∈ (0, 1) is a common discount factor to both agents' types. σ i > 0 (σ i 6 = 1) is the elasticity of intertemporal substitution in consumption and captures heterogeneity in preferences as soon as σ 1 6 = σ 2 . The consumer maximizes his utility function (2), taking the budget constraint (1) as given. We derive the Euler equation:

cit+1 cit = µ β r t+1 + ∆p t+1 p t ¶ σi (3)
and the budget constraint (1), now binding. The transversality condition is also required to be satisfied by the sequence of individual demands: lim t→∞ λ t kit+1 = 0.

Aggregating equations (3) across the individuals of the same type and setting

c i ≡ ci l i , k i ≡ ki l i with i = 1, 2, we find two equations: 1 c it+1 c it = µ β r t+1 + ∆p t+1 p t ¶ σ i (4) r t k it + w t l i = c it + p t (k it+1 -∆k it ) (5)

Heterogeneous producers

The amounts of the pure consumption good and the capital good are given, respectively, by c and y and are produced with two different technologies:

c = F 1 ¡ k 1 , l 1 ¢ and y = F 2 ¡ k 2 , l 2 ¢
, where k j and l j denote the inputs processed in each sector. The aggregate factors demands are given by k = k 1 + k2 and l = l 1 + l 2 . The derivatives of F j with respect to k j and l j are respectively denoted F j 1 and F j 2 . We further assume: Assumption 1 Each production function F j : R 2 + → R + , j = 1, 2, is C 2 , increasing in each argument, concave, homogeneous of degree one and such that, for any x > 0,

F j (0, x) = F j (x, 0) = 0, F j 1 (0, x) = F j 2 (x, 0) = +∞, F j 1 (+∞, x) = F j 2 (x, +∞) = 0.
Under this assumption, we can apply the implicit function theorem to y =

F 2 ¡ k 2 , l 2 ¢ , in order to obtain a locally explicit function k 2 ¡ y, l 2 ¢ . Replacing k 1 = k -k 2 ¡ y, l 2 ¢
and l 1 = ll 2 in the production function of the consumption sector and maximizing the output c with respect to l 2 , given k, y and l, we get eventually the value c = max

l 2 F 1 ¡ k -k 2 ¡ y, l 2 ¢ , l -l 2 ¢ ≡ T (k, y, l)
, which defines a technological frontier. Applying the envelope theorem to the value, we compute the slope of the transformation frontier:

T 1 (k, y, l) = r ≡ r (k, y, l) (6) T 2 (k, y, l) = -p ≡ -p (k, y, l) (7) T 3 (k, y, l) = w ≡ w (k, y, l) (8)
where T i is the partial derivative of the function T with respect to the ith argument. 2 The concavity of the frontier is also characterized, but the reader is referred, for the sake of brevity, to the seminal [START_REF] Benhabib | Competitive Equilibrium Cycles[END_REF] or to [START_REF] Bosi | Competitive Equilibrium Cycles with Endogenous Labor[END_REF] for more computational details:

T 11 = ∂r (k, y, l) /∂k < 0 T 22 = -∂p (k, y, l) /∂y = b 2 T 11 < 0 (9) T 12 = ∂r (k, y, l) /∂y = T 21 = -∂p(k, y, l)/∂k = -bT 11 where b ≡ ¡ k 2 /l 2 -k 1 /l 1 ¢ l 2
/y is the relative capital intensity difference across sectors: a simple, but indispensable parameter to understand the technological heterogeneity. We notice that b determines the sign of T 12 and that the investment good is more capital intensive when b is positive, while, on the opposite, the consumption good is more capital intensive when b is negative.

Equilibrium

An intertemporal equilibrium requires, at each period, that the inputs markets on the one side and the outputs markets on the other side clear. Equilibrium on the capital market requires

k 1 t + k 2 t = k 1t + k 2t = k t , while on the labor market requires l 1 t + l 2 t = l 1 + l 2 = l. Since T (k t , y t , l
) represents the supply of consumption good and y t the produced amount of investment good, both the markets clear, if and only if:

c 1t + c 2t = T (k t , y t , l) (10) y t = k t+1 -∆k t (11)
The consumers' and the producers' optimal plans (4)-( 5) and ( 6)-( 8), respectively, jointly with all these equilibrium conditions, define a competitive equilibrium, the properties of which will be characterized in Section 6.

Dynamic system

Combining consumers' systems (4)-( 5), producers' solutions ( 6)-( 7) and equilibrium conditions ( 10)-( 11), we obtain, eventually, a three-dimensional dynamic system:

y t = k t+1 -∆k t (12) c 1t+1 c 1t = µ β • ∆T 2 (k t+1 , y t+1 , l) -T 1 (k t+1 , y t+1 , l) T 2 (k t , y t , l) ¸ ¶σ 1 (13) T (k t+1 , y t+1 , l) -c 1t+1 T (k t , y t , l) -c 1t = µ β • ∆T 2 (k t+1 , y t+1 , l) -T 1 (k t+1 , y t+1 , l) T 2 (k t , y t , l) ¸ ¶σ 2 (14)
with k 0 > 0 given. These three equations determine the sequence (k t , y t , c 1t ) ∈ R 3 ++ for all t, which governs entirely the dynamics of the economy.

Steady state

A steady state is defined as a solution k t = k t+1 = k, y t = y t+1 = y and c 1t = c 1t+1 = c 1 , that satisfies equations ( 12)-( 14). From equation ( 12), we deduce y = δk. Putting this value in equation ( 13), we obtain:

- T 1 (k, δk, l) T 2 (k, δk, l) = ρ (15)
where ρ ≡ 1/β -∆. As in the two-sector model with a representative household, we get a standard result:

Proposition 1 Under Assumption 1, there exists a unique steady state k > 0 solution of equation ( 15).

Proof. See [START_REF] Becker | Ramsey Equilibrium in a Two-Sector Model with Heterogeneous Households[END_REF].

A solution k corresponds to the Modified Golden Rule: there are no imperfection and dynamics are optimal. The stationary aggregate consumption is defined by c 1 + c 2 = T (k, δk, l).

We observe that neither the aggregate capital k, nor r, p, w are affected by the initial conditions, while c 1 is, as well as the equilibrium transition path. Indeed, using ( 5)-( 8), c 1 can be explicitly written:

c 1 = [(T 1 + δT 2 ) k 1 + T 3 l 1 ].
The stationary sharing (k 1 , k 2 ) and, then, (c 1 , c 2 ) will depend on the initial condition (k 10 , k 20 ).

Before defining the inequality measure, let's introduce a useful notations to simplify formulas. We make use of the share of consumption good consumed by the first agent's type:

c 1 / (c 1 + c 2 ) or, equivalently, κ ≡ c 1 /T ∈ (0, 1) (16)
Market clearing equations allow us to explicitly compute κ. Substituting equations ( 5), ( 6), ( 7) and the equilibrium condition [START_REF] Becker | Ramsey Equilibrium in a Two-Sector Model with Heterogeneous Households[END_REF] in ( 16), we get:

κ = λ -(λ -κ) (1 -δ/ρ) ε 1 (17)
where

ε 1 ≡ kT 1 /T ∈ (0, 1) (18)
while

κ ≡ k 1 / (k 1 + k 2 ) (19) λ ≡ l 1 / (l 1 + l 2 ) (20)
are the poor's shares of capital and labor supply, respectively. We further observe that since k1 ≤ k2 and k i = ki l i , we obtain k 1 /k 2 ≤ l 1 /l 2 or, equivalently,

0 ≤ κ ≤ λ (21) which implies that κ ∈ [0, λ].
5 Inequality measure: the Gini index

In this section, we apply the most widespread of social inequality: the Gini index. The existing papers about the role of heterogeneity on the occurrence of optimal cycles [START_REF] Ghiglino | Inequalities and Fluctuations in a Dynamic General Equilibrium Model[END_REF], [START_REF] Ghiglino | Wealth Inequality and Dynamic Stability[END_REF] and [START_REF] Ghiglino | Wealth Inequality and Macroeconomic Volatility in Two-Sector Economies[END_REF]), instead of taking in account this simple measure, provide more complicate indicators, such as the absolute risk tolerance of the social utility function. 3 Ghiglino (2005) and [START_REF] Ghiglino | Wealth Inequality and Macroeconomic Volatility in Two-Sector Economies[END_REF] show that inequality favors macroeconomic volatility, when the absolute risk tolerance is a strictly convex function. In the opposite case, equality favors macroeconomic volatility. However, in order to measure income inequality in a country and compare this phenomenon over time and among countries, economists usually make use of Lorentz curves and Gini indexes. A Lorentz curve plots the cumulative percentages of total income received against the cumulative percentages of recipients, starting with the poorest agent. The Gini index is calculated as the area between a Lorentz curve and the line of absolute equality, expressed as a percentage of the triangle under the line, and is a pure number, independent on the measure unit. Thus a Gini index of 0 represents perfect equality and the Lorentz curve coincides with the straight line of absolute equality.

To have a rough idea about income inequalities, we provide the Gini index for some countries in accordance with the OECD statistics. In 1989 this measure was about 36% for the OECD countries. The Slovak Republic and Brazil experienced a Gini rate of, respectively, 20% in 1992 and 63% in 1989. France was at 33% in 1989, while United States and China were at, respectively, 40% in 1994 and 41% in 1995.

Formal simplicity and economic pertinence of the index justify why we build it up as a reduced key parameter in order to study the dynamic properties of the model and local bifurcations.

In general, after ranking the agents in [0, l] according to their increasing income and denoting the cumulative income by R (h) we define the Gini index as

g = 1 -2 R l 0 R (h) dh R (l) l (22)
Let's recall that, in our economy, there are l 1 (poor) agents endowed with e k 1 and l 2 = ll 1 (rich) agents endowed with e k 2 > e k 1 . The aggregate income R i of the ith class is given by:

R i ≡ rk i + wl i (23) with R 2 /l 2 > R 1 /l 1 . Therefore, in our linear case, R (h) is given by hR 1 /l 1 , if 0 ≤ x ≤ l 1 , and by R 1 + (h -l 1 ) R 2 /l 2 , if l 1 < h ≤ l.
Replacing this function in ( 22), we obtain a simple formula:

g = λ -R 1 / (R 1 + R 2 ) (24)
Using notations ( 18), ( 19), ( 20) and ( 23), the ratio R 1 / (R 1 + R 2 ) can be computed from the definitions of prices ( 6) and ( 7) on one hand, and the equilibrium conditions ( 10) and ( 11) on the other hand:

R 1 R 1 + R 2 = ρκ + δε 1 κ ρ + δε 1 ( 25
)
where κ is given by (17). Combining ( 24) and ( 25), we find the equilibrium expression of the Gini index:

g = ρε 1 ρ + δε 1 (λ -κ) = ρ ρ -δ ρ ρ + δε 1 (λ -κ) (26)
6 Saddle-path stability and optimal cycles

In this section, we first analyze the stability properties of the steady state and establish the conditions for the occurrence of a flip bifurcation, corresponding to a critical value of capital distribution between the consumers' types. Afterwards, we will connect this result with the Gini index. We linearize the dynamic system ( 12)-( 14) in a neighborhood of the steady state. We obtain the linear form (dk t+1 /k, dy t+1 /y, dc 1t+1 /c 1 ) T = J (dk t /k, dy t /y, dc 1t /c 1 ) T where the Jacobian matrix is explicitly given by:

J =   -β (1 + ∆b) σ 1 δβb (1 + ∆b) σ 1 (ρε 2 ) -1 β (1 + ∆b) A -ε 1 δ [ε 1 /ρ -βb (1 + ∆b) A] κ 1 0 0   -1   -bσ 1 δb 2 σ (ρε 2 ) -1 bA -ε 1 δ ¡ ε 1 /ρ -b 2 A ¢ κ ∆ δ 0   with A ≡ (1 -κ)
σ 2 ρε 2 and ε 2 ≡ kT 11 /T 1 < 0, the elasticity of the interest rate with respect to the capital. The characteristic polynomial gets the form P (α) = α 3 -T α 2 + M α -D, where T , M and D denote, respectively, the trace, the sum of the second order principal minors and the determinant. In order to locate the characteristic roots, we evaluate the polynomial at α = -1, 0, 1.

P (1) = 0 (27) P (0) = -det J = -1/β (28) P (-1) = -2 σ (1 + (1 + ∆) b) [1 + (1/β + ∆) b] -2 (1 + β) E σb (1 + ∆b) -βE (29)
where

σ ≡ κσ 1 + (1 -κ) σ 2 (30)
is the average elasticity of intertemporal substitution and

E ≡ (ε 1 /ε 2 ) / (βρ) 2 < 0 is a technological feature.
The polynomial feature deserves some comments. First, we observe that D = 1/β as in the Cass-Koopmans benchmark and in [START_REF] Benhabib | Competitive Equilibrium Cycles[END_REF]: consumers' heterogeneity and social inequality are not imperfections and the equilibrium remains optimal even when trajectories display cyclical properties. Second, P (1) = 0 entails that α = 1 is the unit root we were finding, deeply related to the existence of a continuum of steady states (stationary consumption sharing), which are determined, as seen above, by the initial capital sharing.

When there are no market imperfections and the investment sector is more intensive in capital, there is no space for optimal cycles. Cycles appear only when the consumption sector becomes more intensive in capital. The robustness of such a conjecture is now confirmed.

In order to present the main proposition, we set the critical values for the average elasticity of intertemporal substitution:

σ * ≡ 2 (1 + β) E [1 + (1 + ∆) b] [1 + (1/β + ∆) b] (31) σ * * ≡ βE b (1 + ∆b)
We observe that, if -1/ (1 + ∆) < b < -1/ (1/β + ∆), then 0 < σ * * < σ * . We are now able to characterize the stability properties of the steady state and the occurrence of endogenous fluctuations:

Proposition 2 In order to ensure prices positivity, we assume b < 1/ρ. Then, the stability properties of the steady state are the following.

(1) If -∞ < b < -1/∆, the steady state is a saddle and there is monotonic convergence.

(2) If -1/∆ < b < -1/ (1 + ∆), the steady state is a saddle. There is monotonic convergence for σ < σ * * and oscillating convergence for σ * * < σ.

(

) If -1/ (1 + ∆) < b < -1/ (1/β + ∆) 3
, the steady state is a saddle and there is monotonic convergence for σ < σ * * , the steady state is a saddle and there is oscillating convergence for σ * * < σ < σ * and the steady state is a source for σ * < σ. When σ goes through σ * , the system undergoes a flip bifurcation.

(4) If -1/ (1/β + ∆) < b < 0, the steady state is a saddle. There is monotonic convergence for σ < σ * * and oscillating convergence for σ * * < σ.

(5) If 0 ≤ b < 1/ρ, the steady state is a saddle and there is monotonic convergence.

Proof. See the Appendix.

We observe that Proposition (2) focuses on the saddle-path stability (monotonic or oscillating), instability and occurrence of persistent cycles. We provide in the Appendix a detailed presentation of the relative position of the eigenvalues with respect to the unit circle in all these cases.

We notice that when the capital good is more capital-intensive (b > 0) or we deal with the one-sector specification (b = 0), case 5 holds and convergence turns out to be monotonic: there is no room for endogenous cycles. On the contrary, when b < 0, endogenous fluctuations can occur. Indeed, optimal cycles arise through a flip bifurcation when b ∈ (-1/ (1 + ∆) , -1/ (1/β + ∆)). As in the two-sector model with a representative agent, the existence of endogenous fluctuations thus requires a more capital intensive consumption good.

The intuition for this result, initially provided by [START_REF] Benhabib | Competitive Equilibrium Cycles[END_REF] in a two-sector model with homogeneous agents, comes from the Rybcszynski and Stolper-Samuelson effects (see also [START_REF] Drugeon | On the Production Possibility Frontier in Multi-Sectorial Economies[END_REF] or [START_REF] Bosi | Competitive Equilibrium Cycles with Endogenous Labor[END_REF]). 4 The occurrence of deterministic fluctuations follows a straightforward rationale. For simplicity, let's consider an instantaneous increase in the capital stock k t . This results in two opposing forces. (i) Since the consumption good is more capital intensive than the investment good, the trade-off in production becomes more favorable to the consumption good. The Rybcszynski theorem implies also a decrease of the output of the capital good y t . This tends to lower the investment and the capital stock in the next period k t+1 . (ii) In the next period the decrease of k t+1 implies again through the Rybcszynski effect an increase of the output of the capital good y t+1 . This mechanism is explained by the fact that the decrease of k t+1 improves the trade-off in production in favor of the investment good which is relatively less intensive in capital. Therefore this tends to increase the investment and the capital stock in period t + 2, k t+2 . Notice also that the rise of y t+1 implies a decrease of the rental rate r t+1 and an increase in w t+1 , correlated through the Stolper-Samuelson effect to an increase of the price p t+1 .

In our analysis, the existence of a cycle of period two is associated to the critical value σ * of the elasticity of intertemporal substitution in consumption. Using equation ( 30), it corresponds to this flip bifurcation a critical value χ * = (σ 2σ * )/(σ 2σ 1 ) of the share of the consumption good consumed by the first type of consumers. Finally, considering equation ( 17), the flip bifurcation can be associated to a level κ * = λ -(λχ * )/[(1δ/ρ)² 1 ] of the wealth or capital inequality between both types of consumers. Now, we want to relate the Gini index with the occurrence of a flip bifurcation. In order to do that, we define a critical Gini index g * associated to the critical bifurcation value σ * . Using equations ( 26), we obtain:

g * = ρ ρ -δ ρ ρ + δε 1 µ λ - σ 2 -σ * σ 2 -σ 1 ¶ (32)
Clearly, for the feasibility of the flip bifurcation κ * , κ * and g * are required to belong to (0, λ). 5

Heterogeneity, inequalities and fluctuations

In this section, we begin by analyzing the role of consumers' heterogeneity on the Gini index evaluated at the steady state. Then, we will show how an increase of heterogeneity in preferences can affect the critical Gini index, i.e. the degree of inequality compatible with endogenous fluctuations.

Heterogeneity and inequality

On the one hand, we observe that the aggregate steady state doesn't depend on σ 1 and σ 2 , because these parameters only affect the two Euler equations ( 13) and ( 14), while the aggregate variable k is given by ( 15) and does not depend on the parameters of heterogeneity. More precisely, the aggregate production (c, y) = (T (k, δk, l) , δk), and derivatives T 1 (k, δk, l) and T 11 (k, δk, l) are not affected by heterogeneity. This entails also that the reduced parameters b, ε 1 , ε 2 and σ * are invariant with respect to the heterogeneity measures.

5 It is easy to show that κ * , κ * and g * are together in (0, λ), if and only if

λ • 1 -ε 1 µ 1 - δ ρ ¶¸< σ 2 -σ * σ 2 -σ 1 < λ
On the other hand, at equilibrium, according to expression (26) and given the poor's share λ in total population, the Gini index can be viewed as a function of the share κ of capital in the poor's hands. In other words, social inequalities are equivalently captured by g or κ. From (26) we derive a straightforward relation between these indices.

Corollary 3 The impact of the equilibrium capital sharing on the Gini index is negative:

∂g ∂κ = - ρε 1 ρ + δε 1 < 0
In order to understand this simple mechanism, we need to keep in mind that capital is the only source of wealth inequality in our model. If we raise the relative part of capital owned by the poor, we lower the social inequality.

In order to stress the role of heterogeneity in preferences, let's now introduce symmetric tastes. For the sake of simplicity, to describe the role of heterogeneity on local dynamics and to capture the impact of heterogeneity on inequality through simple derivatives, we fix the midpoint

σ ≡ (σ 1 + σ 2 ) /2 (33)
and we define an heterogeneity measure

h ≡ |σ 1 -σ 2 | /2 (34)
When σ 1 < σ 2 , current and future consumption are less substitutable for the first type of agent, while they are more substitutable for the second one. The converse holds, when σ 1 > σ 2 . A simple result holds:

Corollary 4 ∂g/∂h = 0.

Not surprisingly, heterogeneity in preferences has no influence on the steady state value of the Gini index, since, as seen above, the elasticities of intertemporal substitution in consumption σ i don't affect the steady state. In what follows, we show why, conversely, heterogeneity in preferences matters for the critical value of the Gini index through which two-period cycles occur.

Heterogeneity and critical inequality

This section is devoted to the role of heterogeneity in consumers' preferences on the occurrence of endogenous cycles. In the main proposition, we provide precise conditions that rely an increase of heterogeneity in preferences to a variation of the inequality level compatible with the emergence of endogenous fluctuations.

Proposition 5 Taking the midpoint σ as given, we raise the heterogeneity measure h and we obtain the following findings.

(1) If σ 1 < σ 2 , then ∂g * /∂h < 0 for σ < σ * , while ∂g * /∂h > 0 for σ > σ * .

(2) If σ 1 > σ 2 , then ∂g * /∂h > 0 for σ < σ * , while ∂g * /∂h < 0 for σ > σ * .

Proof. See the Appendix.

Results in Proposition 5 deserve a clear-cut interpretation of the role of heterogeneity on the critical degree of inequality, giving endogenous cycles. Two main questions arise: (i) what is the effect of the relative behavior of poor and rich (σ 1 ≶ σ 2 )?; (ii) what is the role of the average behavior (σ ≶ σ * )?

First, assume the rich to smooth more easily consumption over time (σ 1 < σ 2 ). A raise in σ-heterogeneity, urges the rich to save more and, thus, to become richer, and the poor to save less. Social inequality increases and a compensating adjustment in terms of the capital sharing κ is required, in order to restore the flip bifurcation, that is to bring back to zero P (-1), the characteristic polynomial evaluated at -1. If σ < σ * (large income effects), the positive adjustment of κ is large enough to entails a decrease in the degree of social inequality, according to Corollary 3, which dominates the initial raise: ∂g * /∂h < 0. Conversely when σ > σ * (small income effects), the increase in κ, restoring the flip, turns out to be smaller and the related decrease in the Gini index doesn't thwart the initial raise: ∂g * /∂h > 0.

Second, when σ 1 > σ 2 , the mechanism is reversed, since, now, more heterogeneity means the poor save more and the rich less, and, in consequence, we observe an initial drop of the inequality index. If σ < σ * , we have a large and negative compensation in κ and a rise in the inequality index larger than the initial fall: ∂g * /∂h > 0. Conversely, small income effects (σ > σ * ) will entail a decrease of the critical inequality degree, at the very end: ∂g * /∂h < 0.

Concluding remarks

Since the seminal work of Benahbib and Nishimura (1985), a well-known condition necessary to the appearance of endogenous cycles in a two-sector economy is a more capital intensive consumption sector. The taking into consideration of heterogeneity on the producers' side is now paralleled by the introduction of heterogeneity on the consumers' side. Our concern is understanding the interplay between these forces. More precisely, after recognizing the origin of fluctuations in productive heterogeneity, the main issue we address, is to see how the second type of diversity affects the bifurcation value through which cycles appear.

When consumers are heterogeneous there is room for thinking and studying social inequality. In our framework we are able to reformulate the bifurcation value as a critical degree of income inequality and, to keep things more understandable than in previous works, we write down the bifurcation as a Gini index.

More heterogeneity in consumers' preferences raises or lowers the critical Gini index according to the relative consumers' ability to smooth consumption over time and to the average elasticity of intertemporal substitution. The economic intuition behind is also provided.

We think that this methodology could apply also to other dynamic models, where the first welfare theorem no longer holds. A further step could be a study on the way heterogeneity amplifies the effects of market imperfections and changes conditions for equilibrium determinacy.

Appendix

Proof of equations ( 6) and ( 7) We maximize F 1 ¡ kk 2 ¡ y, l 2 ¢ , ll 2 ¢ with respect to l 2 , taking k, y and l as parameters. The optimal function l 2 (k, y, l) we obtain, is replaced in k 2 ¡ y, l 2 ¢ to get k 2 ¡ y, l 2 (k, y, l) ¢ , and the value T (k, y, l) is computed, replacing these optimal factor demands k 2 and l 2 into F 1 ¡ kk 2 , ll 2 ¢ . The envelope theorem, applied to the value, gives the slopes of the transformation frontier:

T 1 (k, y, l) = F 1 1 (35) T 2 (k, y, l) = -F 1 1 ∂k 2 /∂y (36) T 3 (k, y, l) = F 1 2 (37)
Profit maximization and perfect factors mobility entail F 1 1 = r = pF 2 1 and F 1 2 = w = pF 2 2 . Substituting these equalities in ( 35)-( 37) and using ∂y/∂k 2 = F 2 1 , eventually we obtain equations ( 6)-( 8).

Proof of Proposition 2

In order to characterize the local stability properties, we need to study the characteristic roots, that is the intersections of the characteristic polynomial P (α) with the axis of abscissas. From equations ( 27) and ( 28), we know that P (1) = 0, P (0) < 0, while the form of the characteristic polynomial implies lim α→+∞ P (α) = +∞ and lim α→-∞ P (α) = -∞. Let, without loss of generality, α 1 = 1, be the first eigenvalue. Since α 1 α 2 α 3 = 1/β, then the product of the other eigenvalues is α 2 α 3 > 1 and at least one characteristic root, say α 3 , lies outside the unit circle. Since the polynomial P is continuous, P (-1) > 0 implies that α 2 , lies in (-1, 0) and α 3 < -1. The steady state is a saddle point. If P (-1) < 0, then either (i) α 2 , α 3 < -1 or (ii) 0 < α 2 < 1 < α 3 or (iii) α 2 , α 3 > 1 or (iv) α 2 , α 3 are complex and conjugated with |α 2 | = |α 3 | > 1. Case (ii) is characterized by saddle point stability and monotonic convergence, while the other cases correspond to source configurations. Case (ii) is equivalent to a negative slope of the characteristic polynomial at 1, that is P 0 (1) < 0. When P (-1) = 0, α 2 goes through -1 and the system generically undergoes a flip bifurcation.

First, assume σb (1 + ∆b) -βE > 0 (38) Inequality (38) includes the case of a more capital intensive investment sector (b > 0) and the one-sector case (b = 0). Since b < 1/ρ, in order to ensure prices positivity (see [START_REF] Benhabib | Competitive Equilibrium Cycles[END_REF]) and always 1/ρ < 1/δ, (38) implies, in general,

P 0 (1) = - σ (1 -ρb) (1 -δb) σb (1 + ∆b) -βE < 0
Proposition 2 sums up all these findings, without detailing the relative position of the eigenvalues, but discussing the crucial points of saddle-path stability and monotonic convergence.

Proof of Proposition 5 Let σ be fixed according to (33) and g * be defined by (32). If σ 1 < σ 2 , then h = (σ 2σ 1 ) /2 and (σ 2σ * ) / (σ 2σ 1 ) = (h + σσ * ) / (2h). We have:

∂g * ∂h = 1 2 ρ ρ -δ ρ ρ + δε 1 σ -σ * h 2
We deduce that ∂g * /∂h < 0, if and only if σ < σ * . Similarly, if σ 1 > σ 2 , h = (σ 1σ 2 ) /2 and (σ 2σ * ) / (σ 2σ 1 ) = (h + σ * -σ) / (2h), and therefore:

∂g * ∂h = 1 2 ρ ρ -δ ρ ρ + δε 1 σ * - σ h 2
We conclude that ∂g * /∂h > 0, if and only if σ < σ * .

Notice that l

is defined by l 2 ≡ l -l 1 .

The derivation of equations (6), (7) and (8) is provided in the Appendix.

The social absolute risk tolerance is defined as the inverse of the social absolute risk aversion.

The Rybcszynski effect focuses on quantities to explain the emergence of cycles: if the endowment of a factor increases, then the production of the good that uses this factor relatively more intensively raises and the production of the other good decreases. The Stolper-Samuelson effect is a dual theory referring to price mechanics: when the price of the good j increases, the equilibrium price of the factor more intensively used in the production of good j increases too, while the price of the other factor decreases.

Therefore, we obtain 0 < α 2 < α 1 = 1 < α 3 , that is a saddle-path stability, jointly with a monotonic convergence.

Second, suppose, to the converse:

Then P (-1) > 0, if and only if

In order to locate the eigenvalues with respect to the unit circle, we require to know at first the position of the relative capital intensity difference b with respect to its ranked critical points: -∞, -1/∆, -1/ (1 + ∆) , -1/ (1/β + ∆) , 0, 1/ρ, and, second, the position of the aggregate elasticity of intertemporal substitution σ with respect to its own critical points σ * , σ * * . Crossing the intervals, we obtain an exhaustive list of cases and subcases.

(1) -∞ < b < -1/∆. Then (38) holds, that is 0 < α 2 < α 1 = 1 < α 3 (saddle-path stability and monotonic convergence).

(

Two subcases hold according to the values of the elasticity of intertemporal substitution.

(2.1) If σ < σ * * , then (38) holds, that is 0 < α 2 < α 1 = 1 < α 3 (saddle-path stability and monotonic convergence).

(2.2) If σ * * < σ, then ( 39) and ( 40) are satisfied: P (-1) > 0 and α 3 < -1 < α 2 < 0 < α 1 = 1 (saddle-path stability and oscillating convergence).

(

In this case, under b < 1/ρ (prices positivity) always σ * * < σ * . Three subcases matter according to the values of the elasticity of intertemporal substitution.

(3.1) If σ < σ * * , then (38) is satisfied. Therefore 0 < α 2 < α 1 = 1 < α 3 (saddle-path stability and monotonic convergence).

(3.2) If σ * * < σ < σ * , then (39) and ( 40) are satisfied, that is P (-1) > 0 and α 3 < -1 < α 2 < 0 < α 1 = 1 (saddle-path stability and oscillating convergence).

(3.3) If σ * * < σ * < σ, then (39) holds, while (40) doesn't. Therefore P (-1) < 0 and, since P 0 (1) > 0, then, as above, either (i) α 2 , α 3 < -1 or (iii) α 2 , α 3 > 1 or (iv) α 2 , α 3 are complex and conjugated with |α 2 | = |α 3 | > 1. All these cases correspond to source configurations. When σ goes through σ * , the system undergoes a flip bifurcation. By continuity with respect to the case 3.2, we have that, in a neighborhood of σ * , as soon as σ increases through σ * , then α 2 decreases through -1, while α 3 < -1 and α 1 = 1.

(4) Let -1/ (1/β + ∆) < b < 0. Two subcases matter according to the values of the elasticity of intertemporal substitution.

(4.1) If σ < σ * * , then (38) holds, that is 0 < α 2 < α 1 = 1 < α 3 (saddle-path stability and monotonic convergence).

(4.2) If σ * * < σ, then (39) and ( 40) are satisfied, that is P (-1) > 0 and α 3 < -1 < α 2 < 0 < α 1 = 1 (saddle-path stability and oscillating convergence).

(5) Finally, set 0 ≤ b < 1/ρ. Then (38) holds as in point (1). Therefore: 0 < α 2 < α 1 = 1 < α 3 : saddle-path stability and monotonic convergence.


H Hamidane
email: hamidanehafsa@gmail.com

S El Faiz

M Guerbaoui

A Ed -Dahak

A Lachhab

B Bouchikhi

IMAT3E'18 International Meeting on Advanced Technologies in Energy and Electrical Engineering Pole placement problem for linear systems with constrained control. (An Overview)

An overview of some of the methods used for modern systems and control theory related to the pole placement problems is given. These methods are used in order to develop and synthesis stabilizing controller which studies the stability of the system. As a particular case; linear time invariant systems in continuous time domain

( ) ( ) ( ) x t Ax t Bu t   (1)
Where To investigate the stability of such linear systems, we reconsider the pole placement problems methods used in order to choose a state feedback control law:

; ( )

m n u Fx F    (2)
such that the closed loop system

( ) ( ) ( ) x t A BF x t   (3)
has a specified poles, regarding the constraints notion. These constraints are usually represented as inequalities that can affect the systems performance if they are not taken into consideration when synthesising and developing the controller. For that aim it is preferable to develop the exact mathematical models representing these saturations that can affect the systems state or control law

In the present framework we review the case where the control law u of the system is constrained to evolve in the domain m R  Ω defined as:

IMAT3E'18

International Meeting on Advanced Technologies in Energy and Electrical Engineering

min max min max { / - , ,- 0} m u R u u u u u      Ω
Where Ω represents a non symmetric polyhedral set.

In this regard a necessary concept to take into account when involving constraints concept, is the positive invariance approach [1, 2, 4, 8], thus the main objective of such an approach is to obtain a stabilizing regulator F as a unique solution to the nonlinear algebraic equation:

FA FBF HF   (4)
Such that the model of equation ( 3) remains valid. Based on the cited equation ( 4) an algorithm developed in [START_REF] Elfaiz | Commande des systèmes linéaires sous contraintes:stabilisation et placement de poles robuste via LMI[END_REF] will be presented as well to describe the purpose of such an equation in synthesizing the control law of linear constrained systems.

In addition to the mentioned non linear equation, another type of equations that contributes positively in the development of the pole placement solutions is known as

Sylvester equation [START_REF] Benzaouia | Stabilization of linear systems with saturation: a Sylvester equation approach[END_REF][START_REF] Benzaouia | The resolution of equation XA+ XBX= HX and the pole assignment problem[END_REF][START_REF] Wang | Constraint generalized Sylvester matrix equations[END_REF][START_REF] Elfaiz | Commande des systèmes linéaires sous contraintes:stabilisation et placement de poles robuste via LMI[END_REF] , this equation depends mainly on equation ( 4) to find a control law F for the linear system and it has the form:

AX XH B    (5)
X is the unique solution of equation ( 5), For both non linear and Sylvester equations, the reduced order technique [START_REF] Castelan | A reduced-order framework applied to linear systems with constrained controls[END_REF] is also engaged as a method that gives the possibility to solve some classes of constrained linear systems by the projection of the trajectories of the original system onto a subspace associated

with the undesirable open-loop eigenvalues which enable to decreases the complexity of some computations.

We also introduce linear matrix inequality (LMI) technique [START_REF] Ammar | Automatic tuning of paper machines cross-direction controllers via linear matrix inequalities[END_REF][START_REF] Boyd | Linear matrix inequalities in system and control theory[END_REF] due to its efficiency and simplicity in resolving the pole placement and concomitant problems; in an efficient numerical manner.

  , x and u represent respectively the state and control vectors of the system; are engaged in this framework to show dynamic systems pole placement development in the state space form.

  this order another another algorithm presented in will be considered and clarified.

Engineering

Finally and for the same aim some of the numerical examples treated in [START_REF] Elfaiz | Commande des systèmes linéaires sous contraintes:stabilisation et placement de poles robuste via LMI[END_REF] corroborating and illustrating the usefulness of the proposed synthesis techniques are also suggested. In this order we involve some of the most recently used softwares such Matlab, including latest and various LMI Toolboxes, one to be cite Yalmip Toolbox,to reexecute these examples.


Félix Darve

Cédric Lambert

Felix Darve

Cedric Lambert

Continuous and discrete modelling of failure in geomechanics

come

Continuous and discrete modeling of failure in geomaterials

Introduction

One of the most basic question in geomechanics is how to define, analyse and simulate failure. The general way to introduce the phenomenon of failure is to note that in ele mentary laboratory experiments, some loading paths are leading to limit stress states. Along these paths, these stress states can not be passed : experimentally they constitute asymptotical limit states. By restricting our study to elasto-plastic behaviors, the general expression of rate independent constitutive equations is given by :

(1.1)

where the elasto-plastic tangent operator � depends on the previous stress-strain history (characterized by state variables and memory parameters x ) and on the direction of dE (characterized by unit matrix ;:, = �/11�11).

In the six-dimensional associated spaces it comes :

(cx,(J = 1, ... , 6) (1.2) v , = de,/11 � 11 and 11 � 11 = )deo:deo:.

From a mathematical point of view, limit stress points are characterized by : du = 0 with 11 � 11 =1-0

Condition ( 1.3) implies :

det N ( v ) = 0 and N (v)de = 0 �x � "' x ,... .. ., ,... .. ..,

(1.3)

( 1.4)

If the constitutive relation is incrementally linear, f! is independent on �. The first equa tion of (1.4) represents precisely the plastic limit criterion because this is the equation of an hyper-surface in the six-dimensional stress space. The second equation represents the material flow rule, since this equation gives the direction of de when the plastic limit � surface has been reached. For incrementally piecewise linear constitutive relations, where N depends in a discontinuous manner from � (inside "tensorial zones" , see the chapter by Darve and Servant), det N ( v ) = 0 is still the equation of an hyper-surface in the stress space, possibly in-

� x �
volving several plastic mechanisms. Besides, N ( v )de = 0 is a singular generalized flow

� X � �
rule with vertices which are locally constituted by pyramids [START_REF] Hill | Eigenmodal deformations in elastic-plastic continua[END_REF]. This classical elasto-plastic view of failure has given rise to the wide domain of the "limit analysis" , where the whole body at a failure state is assumed to have reached the plastic limit condition.

However, in practice, various modes of failure are encountered : localized modes in shear bands [START_REF] Rice | The localization of plastic deformation[END_REF][START_REF] Vardoulakis | An incrementally non-linear constitutive law of second order and its ap plication to localization[END_REF]Darve F. , 1984), in compaction bands or in dilation bands, diffuse modes due to geometric instabilities (buckling, ... ) or with a chaotic displacement field, ... For non-associated materials (essentially because of the non-symmetry of the elasto-plastic matrix f!J, the elasto-plastic theory shows effectively that, by considering proper bifurcation criteria, these localized or diffuse modes of failure can appear before the plastic limit condition [START_REF] Bigoni | Uniqueness and Localization-I. Associative and Nonasso ciative Elastoplasticity[END_REF].

In relation with these various modes of failure, different bifurcation criteria are available.

With respect to shear band formation by plastic strain localizations, Rice's criterion (1976) based on the description of such an incipient shear band of normal n corresponds to vanishing values of the so-called "acoustic tensor"

(1 .5)

For non-associated materials, equation (1.5) can be satisfied before plasticity criterion

(1.4). This has been verified experimentally for dense sand [START_REF] Desrues | Shear band initiation in granular materials : experimentation and theory[END_REF]. Strain localization corresponds to a bifurcation of the strain mode from a diffuse one to a strictly discontinuous one. This kind of bifurcation can be called "discontinuous bifurcation" . "Continuous bifurcations" [START_REF] Darve | Homogeneous bifurcation in soils[END_REF] opposed to the previous one correspond also to a failure mode, but without strain localization. Such failure mode is called "diffuse failure" and this is the response path which is subjected to a bifurcation with a loss of constitutive uniqueness at the bifurcation point. According to the control mode (stress, strain or mixed control) of the loading path, different response paths are possible from the bifurcation point. For certain control modes, the loading is no more "controllable" in Nova's sense [START_REF] Nova | Controllability of the incremental response of soils specimens subjected to arbitrary loading programmes[END_REF]. These continuous bifurcations and the related diffuse modes of failure can be detected by Hill's sufficient condition of stability [START_REF] Hill | A general theory of uniqueness and stability in elastic-plastic solids[END_REF] which corresponds, for unstable states, to vanishing values of second order work, that is to say vanishing values of the determinant of the symmetric part of matrix N for incrementally linear constitutive relations. Section 2 of this chapter is devoted to recall briefly Lyapunov's definition of stability [START_REF] Lyapunov | Probleme general de la stabilite des mouvements[END_REF]) and Hill's condition, then to apply them to the well known case of the failure of undrained loose sands. The generalization to axisymmetric paths and to plane strain conditions allows to determine numerically and analytically the "unstable" domain (i.e. the failure domain) and, inside this domain, the "unstable" stress directions which form cones for both these cases.

In section 3, this analysis is applied to axisymmetric q-constant loading path, since some experiments (Chu .J. & Leong W.K. ( 2003)) have shown failure modes without localization along these paths in the case of loose sands.

Section 2 and 3 have illustrated the application of Hill's condition in the fr amework of continuum mechanics. It is also possible to consider a discrete form of second order work [START_REF] Mandel | Conditions de stabilite et postulat de Drucker[END_REF] and to apply it to boundary value problems described by discrete mechanics. In this perspective the most known examples of typical failures are the so called "grain avalanches" . This phenomenon has been extensively studied specially by physicists of granular media (see Hermann H . .J. et al. (1998) for example). It has been also considered as a paradigm for the concept of "self-organized criticality" [START_REF] Bak | How nature works[END_REF]. In section 4 of this paper, grain avalanches are analyzed in the strict framework of discrete mechanics. The, more or less local, failures will be characterized by bursts of kinetic energy (computed by a discrete element method). Strong spatial and temporal correlations between these bursts and the peaks of negative values of second order work are exhibited. These results represent another illustration of the deep link which relates certain failure modes and second order work, which could constitute a novel indicator of failure.

Material instabilities

The most basic definition of stability was certainly proposed by [START_REF] Lyapunov | Probleme general de la stabilite des mouvements[END_REF]. Applied to the field of continuum mechanics it states that :

"A stress-strain state, for a given material after a given loading history, is called stable, if for every positive scalar E, a positive number 7](E) exists such that for all incremental loading bounded by 7], the associated responses remain bounded by E"

According to this definition, all limit stress states (as defined in the introduction) are unstable. Indeed if a "small" incremental reverse stress from a limit stress state is considered, a "small" (more or less elastic) response will be induced, while a "small" incremental additional stress beyond the limit state produces a "large" strain response. Lyapunov's definition of stability is interesting because it states clearly that material instabilities can be conjectured in elasto-plastic media. Since some limit stress states are strictly inside Mohr-Coulomb plastic limit surface (because of non-associativity) -an example is given in section 2.1 -some instabilities can be expected before reaching Mohr-Coulomb plastic limit condition. But, even if this first conclusion is essential, Lyapunov's definition is not a proper tool to investigate instabilities in geomaterials. For that, we are now considering Hill's sufficient condition of stability [START_REF] Hill | A general theory of uniqueness and stability in elastic-plastic solids[END_REF]. According to Hill, a stress-strain state is unstable if it exists at least one loading direction which can be pursued in an infinitesimal manner without any input of energy from outside. In practice it means that the deformation process will continue "on one's own" under dead loads . Indeed, in some situations, failure needs some energy from outside of soil boby (for example, the weight of a structure on a foundation), while in other situations energy is not needed (landslides, rockfalls, ... ). Finally Hill's condition of stability states that a stress-strain state is stable if, for all If the equivalence between Lyapunov's definition and "non-controllable" paths in Nova's sense [START_REF] Nova | Controllability of the incremental response of soils specimens subjected to arbitrary loading programmes[END_REF] is accepted, then non-controllability and Hill's condition are equivalent in classical elasto-plasticity [START_REF] Nova | Controllability of the incremental response of soils specimens subjected to arbitrary loading programmes[END_REF]. For hypoplastic relations (as defined by [START_REF] Chambon | CLoE, a new rate-type consti tutive model for geomaterials. Theoretical basis and implementation[END_REF]) the positiveness of second order work assures the uniqueness [START_REF] Chambon | Existence and uniqueness theorems for boundary value problems involving incrementally non-linear models[END_REF]. In the general case of incrementally non-linear constitutive relations, counter-example might exist. In the next sub-section, all that is illustrated by the well known case of the undrained triaxial behavior of loose sands.

U ndrained paths on loose sands

The undrained triaxial behavior of loose sands is typically represented in figure 1, where q is passing through a maximum which is situated strictly inside Mohr-Coulomb plastic limit surface. The experiments show that, if a "small" additional axial force is applied at q peak, there is a sudden failure of the whole sample without any localization pattern. Because of this sudden failure for a small perturbation at q peak, q peak is an unstable state according to Lyapunov's definition. For axisymmetric conditions and coinciding principal axes, the second order work is equal to :

(2.2)

,.

� �---.. ----�--� � � --. �00 --� .� ----------� p' (kPa) (2.3)

Thus according to Hill's condition , q peak has to be considered as an unstable state ( [START_REF] Darve | Constitutive instabilities in incrementally non-linear mod elling[END_REF]).

Let us see now that we can define at q peak a generalized plasticity criterion and a generalized flow rule (as defined in equation ( 1.4)).

Because q-E1 and Ev-E3 (with : Ev = E1 + 2c3) are conjugate variables with respect to energy : a1E1 + 2a3E3 = qE1 + Eva3, the constitutive relation can be written under the following form :

[::J = � [�:�] (2.4)
For an undrained loading, the isochoric condition is fulfilled, thus :

dEv = 0 (2.5)
As dq is vanishing at q peak, it comes : detP = 0 and

(2.6)

These relations generalize equations (1.4). The failure mode can be viewed as the eigen mode associated to the vanishing eigenvalue of matrix P.

In conclusion it is emphasized that q peak is a proper failure state (satisfying all the cri teria associated to the notion of failure) which is situated strictly inside Mohr-Coulomb plastic limit surface and is related to a diffuse mode of failure.

In several previous papers [START_REF] Darve | Instabilities in Granular Materials and Application to Land slides[END_REF], 2001;[START_REF] Laouafa | Modelling of slope failure by a material instability mecha nism[END_REF] we have exhibited numerically the boundaries of the unstable domain in ax isymmetric conditions, in plane strain conditions and more generally for 3D conditions. Besides, second order work is typically a directional quantity in the stress space. Indeed, if the constitutive relation is given by : with (2.7)

(2.8)

Thus there are cones of unstable stress directions, where the second order work takes negative values.

Next sub-section 2.2 is devoted to the analytical determination of the equations of the unstable domain and of the cones of unstable stress directions in the case of an incre mentally piecewise linear constitutive relation for axisymmetric conditions. in 8 "tensorial zones" . Let us restrict the presentation to the axisymmetric case with 4 tensorial zones. In such a case, the expression of the model is given by : with :

fO V g l± l -v £--;:;r E3 1-V3 ± � 3 (2 .9) (2.10)
where index (+) means "compression" (da; > 0) and (-) "extension" (dai < 0) in the axial direction "i" of "generalized" triaxial paths (two lateral constant stresses aj and ak)• Ei are "generalized" tangent Yo ung moduli and V/ "generalized" tangent Poisson's ratios defined by : According to the signs of da1 and da3 (4 possibilities), there are 4 linear relations between � and d;:, (see details in chapter 1).

In order to investigate the sign of second order work at a given stress-strain state in axisymmetric conditions, it is necessary to consider all the directions for example in the axisymmetric strain space (i.e. Rendulic Plane). For that it is introduced axisymmetric proportional strain loading paths defined incrementally by : positive constant dc3 ( axisymmetric condition) 0 (R constant for a given path) R = 1 corresponds to the undrained path which is characterized by the isochoric condition (no volume variation) .

For 0 < R < 1, the imposed path is dilatant, while for 1 < R it is contractant . For this class of paths, the proper conjugate variables are (0"1 -0"3j R) versus E1 and (El+ 2RE3) versus (0"3/R) since: 0"1f1 + 20"3f3 = (0"1-0"3jR)E1 + 0"3(f1 + 2RE3)/R. We have verified [START_REF] Laouafa | Modelling of slope failure by a material instability mecha nism[END_REF] numerically that , for a range of small R values, curve (0"1 -0"3j R) versus E1 is passing through a maximum. This range of R values depends, of course, on the initial density of the sand (R must be smaller for a higher density). These curves can be considered as generalizations of q-E1 diagrams for undrained triaxial paths.

It is easy to verify that, at ( 0"1 -0"3j R) peak, Lyapunov's definition is verified.

As regards Hill's condition, the second order work is equal to :

(2.11)

With the constraint : dc1 + 2Rdc3 = 0, it is clear that the second order work is vanishing at (0"1 -0"3/ R) peak and, after it, takes negative values. The constitutive equation can be written like that :

(2.12) if we assume that d0"1 and d0"3 are negative (as given by the numerical results). At ( 0"1 -0"3/ R) peak, the bifurcation criterion is given by :

detS = 0 (2.13)
if � is the constitutive matrix of relation (2.12). (2.13) gives the following equation in R

E-

( E -) 2 -1 (1 -v; 3 -) R 2 -2 11, 3 -+ _ 1 y: l -R + 1 = 0 E- 3 1 E- 3 3 3 (2.14)
which is an equation of second degree in R. It has two real solutions in R, if discriminant � is strictly positive :

(2.15)

or equivalently :

( 1-v; 3 --2V?-v; 1 -) 2 .6. = -2 3 _ 1 3 det Q 8 > 0 E 3 �
where Q 8 is the symmetric part of the constitutive matrix Q given by : which relates [ �� c3 ] to [ ��� 3 ] for da1 < 0 and da3 < 0 (see equation 2. 10).

The unstable domain is thus characterized by : det Q 8 < 0 and its boundary by : detQ 8 = 0 Indeed, for incrementally linear constitutive relations, The discussion is thus more intricate, because the first unstable stress direction (which is given by the unique solution of equation (2. 14) when discriminant .6. is vanishing) must belong to the tensorial zone defined by the signs of da1 and da3. Anyway the equation of the boundary of the unstable domain (if it exists) is given by : detQ 8 = 0 where Q 8 is the symmetric part of constitutive matrix Q associated to the proper tensorial zone. At this boundary, the unique solution of (2. 14) is equal to :

"-' rv rv !"V I"V rv I"V I"V rv rv f"'V I"V In such a simple case, d 2 W is vanishing only if det Q 8 is vanishing (if the eigenvalues
R --&1 - E! Vl-+ E 3V1 3 - c -2&3 -2E 1 (1-V3 3

-)

In Rendulic plane, the first unstable stress direction is then given by :

(2. 17)

(2. 18) with 8a1 < 0 and 8a3 < 0. These stress directions are plotted in figure 5.

The associated mode of rupture is characterized by the following rupture rule :

8 [ dE1 ] _ [0] � da3/R -0 with detS = 0
where S is the constitutive matrix defined by equation (2.12). It comes :

� E!dE1 + (2�[ vl--�) da3 = o (2.19)
At the bifurcation point (or equivalently at the limit stress-strain state) characterized by

{ da1 -da3jR = 0 dE1 + 2RdE3 = 0
there are an infinite number of solutions (loss of uniqueness), but all these solutions must fulfilled constraint (2.19).

Note that these conclusions are perfectly similar to the plastic limit conditions discussed in the introduction (see general equations The strain proportional loading paths have been simulated for various values of R. Some maxima for (a1-a3/R) have been found for a certain range of R values (see figure 2). Finally the cones of unstable stress directions have been plotted for a dense sand (figure 3) and a loose sand (figure 4). Note that for high stress levels there are two cones of unstable stress directions. Indeed the previous analytical results have just shown that there are possibly two solutions in R (what means one cone) inside one tensorial zone. Thus several cones can exist, when there are several tensorial zones.

In the next section, the case of plain strain conditions is considered.

Plane strain conditions

For plane strain conditions (E2 = 0) and coinciding stress-strain principal axes, the same reasoning can be developed. The octo-linear constitutive model takes the following form :

(2.20) ,200

'""' .,. ������--�--�--�--�--� � . � ""' � � � � � -� p (kPa) •200
•-••• , . .... . .. .. , .......... , . ......... , . .. . .. ... ;.:" ' ' ' ' ' ' ' ' : !/' . , .

•�•• � ••• ••• •• • , m•;7$ ••f •• • j •:• •m T � z:J��-:T • . . • •= 'o .. '"' '"" '""
... ... "t',

•: Following indices ( +,-), 2.20 corresponds to 8 different linear relations between (dc1, dc3) and ( d0"1, d0"3) by taking into account the sign of d0" 2 .

• •• : ... ...-r-r-=--t . ' • • • • •t>;•; •• -�•r• "� . . ' -� .. . �•�• � • - T• • • - --�-�-;/'
Let us consider proportional strain loading paths defined by : positive constant 0 0 (R constant for a given path)

The conjugate variables are (0"1 --']t) versus c1 and (c1 + Rc3) versus ( -]t ).

The second order work is equal to :

(2.21)

and will vanish at (0"1 -�) peak. If the tensorial zone d0"1 < 0, d0" 2 < 0, d0"3 < 0 is considered (this assumption has to be verified a posteriori), the octo-linear model can be written under the following form :

(2.22)

with d = 1-v l 2 �v 2 1 � -v l 3 �vl� -v 2 3 �v 3 2� -v l 2�v 2 3 �v 3 1 � -v 2 l �v l 3 �v3 2 �0
Taking into account the loading constraint dc1 + Rdc3 = 0, and, at (0"1 --]t) peak, it comes the bifurcation criterion given by the vanishing determinant of constitutive matrix of relation 2.22 :

E1(1-V 2 3�V 3 2�)R 2-[E;3(V 1 3� + V 1 2�V 2 3�) + E!(V}� + �' 3 2�V 2 1�)] R +E3(1-V 1 2 �V 2 1 �) = o (2.

23)

As for axisymetric conditions an equation of second degree in R is found, which implies the existence of cones of unstable directions in a given tensorial zone (if the unstable directions belong to the tensorial zone which has been assumed a priori), when the discriminant is positive. This discriminant is equal to :

� = (E3)2(v l 3 � + v l 2 �v 2 3 �) 2 + (E 1)2(v 3 1 � + v 3 2 �vl�) 2 -2E1 E3[d + (1-V 2 3 �V 3 2 �)( 1-V 1 2 �V 2 1 �) ]
which can be put under the following form :

(2.24)

(2.25)

where pS is the symmetric part of the constitutive matrix of relation 2.20 (with indices � "-" ).

Once more, the boundary of the unstable domain is given by : detP5 = 0 (2.26)

(if the unstable directions correspond to da1 < 0 and da3 < 0)

The first unstable direction corresponds to :

or

R --( _ dE _ l ) -_ E -"-3 ---'-( V, � l "- 3 -_ + _ V, _, ? _-_ v; ....: :2 3 _ -_:_ ) _ + _ E ..,..=-1....:. . ( V -<2 3 '-- l -_ + _V .. ..:e l _-_ v; -"'-2 1 _ ---'--) c - dE3 c - 2E! (l-V 2 3 -V 3 2 -)
(2 .27)

Figures 5,6 and 7 present the "first" unstable stress states (according to the level of stress deviator) respectively for the axisymmetric conditions (figure 5), for plane strain conditions (figure 6) and for 3D conditions in a deviatoric stress plane (figure 7). The stress curve or surface, which is obtained by linking these stress states, corresponds to the boundary of the "unstable" domain. Inside this stress unstable domain one can find various kinds of bifurcations and losses of uniqueness, which means various failure modes by material instabilities (localized failure, diffuse failure ... ) and geometric instabilities.

These failure modes depend on the previous stress-strain history, on the current incre mental stress direction and on the loading mode (stress controlled, strain controlled ... ). (2.29)

It follows that there are cones of unstable stress directions, as plotted in figures 3 and 4.

The first unstable direction is progressively opening into a cone when the stress deviator ,. Is mcreasing. The limit directions of these cones are given by th solutions of equation 2.14 for axisymetric conditions and the solutions of 2.23 in plane strain conditions. It is worth to note in figures 3 and 4 that the vertical stress direction (which corresponds to drained triaxial tests) is never inside the cones. Along drained triaxial compressions or extensions it is impossible to obtain diffuse failure, which is a well known experimental fact.

3 q-constant loading paths an example

Analysis

At a given stress-strain state, q-constant axisymmetric loading paths consist in main taining q constant by an incrementally isotropic unloading defined by : dCY1 = dCY 2 = With our octo-linear model, the bifurcation criterion is given by :

E- 2 -1 (1-v 3 --v 1 -) + 1 -2v 3 -= o E- 3 3 1 3
The rupture rule corresponds to :

E!dcl + (2!f vl--1 ) da3 = 0 (3.6) (3.7)
The unstable stress direction is obviously given by : da1 = da3 and it is parallel to the hydrostatic line.

Finally one question is remaining open : are the volume variations possessing a maximum for a certain range of initial densities and a certain range of stress levels ? To answer to this question in the next subsection, we have simulated q-constant loading paths with our constitutive relations.

3.2

Numerical modelling q-constant paths have been simulated with our octo-linear model and with the non linear one. Figures 8 and9 present the results for loose Hostun sand and for initial isotropic pressures respectively equal to 100 kPa and 500 kPa. When the second or der work is negative, the computed points are replaced by different symbols (triangles, squares, circles, depending on the value of q, see figure 8). The diagrams are composed of two parts. A first part is constituted by a triaxial compression until reaching various q values. Then q is maintained constant. These results allow to validate our previous assumptions. For example, lateral pressure is decreasing monotonously. On the contrary the volume variations are indeed possessing a maximum. Of course it is verified that from this maximum, the second order work is negative.

Figure 10 shows the results obtained with dense Hostun sand. The same comments as previously are pertinent, except the fact that the volume variations do not reach any maximum. Note that the computations are stopped when Mohr-Coulomb plastic limit condition is reached. The second order work is always strictly positive. These results could have been drawn from figures 3 and 4. q-constant loading paths are parallel to the hydrostatic axis. We see that this direction is included inside the cone of unstable directions for low qjp values for loose Hostun sand and not at all for any qjp value for the dense one. Modelling has been also performed with the incrementally non-linear model with similar conclusions (see figure 11 for such a comparison).

Several conclusions can be drawn. The most important one is certainly the fact that a diffuse mode of failure, which was repeatedly noted in experiments (Chu J . & Leong W.K., 2003), can be properly analyzed by Hill's condition of stability. It has been also shown that the bifurcation criterion and the failure rule are the same for undrained (i.e. isochoric) paths and q-constant paths. Finally a generalization of the notion of limit state has been brought out. For q-constant paths the stresses are monotonously decreasing (dCJ1 = dCJ 2 = dCJ3 < 0) while the dilatancy is reaching a maxi mum value. If one injects at this maximum a small amount of water inside the sample, a sudden failure is occurring without any localization pattern. Thus we can expect the existence of a limit surface in the strain space (dual of the classical limit surface in the stress space). Indeed the strains are limited on the hydrostatic axis or for one-dimensional compressions. Next section is now devoted to investigate a possible application of Hill's condition in discrete mechanics. •,!%) The initial isotropic pressure is equal to 100 kPa.

•� �0 �•, •0• <I ••� . z •U .. . .. "" � . .., .. ••• 0 ,, �6 : G A •o ..

Grain avalanches

Granular media are characterized by a double nature : continuous and discrete.

For some phenomena (experiments on homogeneous samples, in situ behavior at large scale, ... ) they behave like continuous media , the assumption of matter continuity repre sents a proper description, stresses, strains and constitutive relations are pertinent tools and the framework of continuum mechanics is adequate.

For other phenomena (granular segregation, mixing of two granular populations, ... ) they behave like discrete media, each grain has to be considered as a specific individual, forces, displacements, rotations and local granular interactions are pertinent tools, while the framework of discrete mechanics is convenient. Grain avalanches along a slope are certainly a consequence of discrete nature of granular media. Thus, in this section, we have to introduce discrete notions, discrete tools and finally to perform discrete compu tations.

1 Discrete second order work and kinetic energy

In this section we consider only two-dimensional granular materials which are consti tuted by piled cylinders (the so-called Schneebeli material) . All the variables are thus defined in a two-dimensional space.

A given grain is subjected to several forces due to the other grains in contact and due

----+ �
to gravity. The resultant force is called F. Under the force increment dF the grain is

----+
moving of dl. All these forces are also developing torques along the cylinder axes. The resultant torque is noted C, whose consequence is a rotation fl around the axis.

Thus the discrete second order work is defined by (Mandel J. , 1966)

d 2 W = dF.df + dCdfl ( 4.1)
The objective is to investigate a possible link between discrete second order work and failure in discrete media, as an extension of the previous results which have shown such a link between second order work and diffuse failure.

In discrete media, failure at free surface is generally characterized by granular avalanches. During these avalanches the kinetic energy is varying a lot from zero values to maxima, then again to zero values. Thus kinetic energy has been chosen as a proper variable to characterize a grain avalanche from a mechanical point of view. For a given grain its expression is :

1

1 e = -mv 2 + -Iw 2 c 2 2 ( 4.2)
where m is the grain mass, v its speed, I the inertia moment around the cylinder axis and w the rotation speed.

Second order work and kinetic energy have been defined for a given grain. It is also interesting to consider their values for the whole assembly of grains as follows :

(4.3)

where N is the total number of grains.

Discrete material and discrete numerical method

The granular medium is a two-dimensional material constituted by wood piled cylin ders. The diameters are equal to 13, 18 and 28 mm. There are around 350 particles for a given experiment simulated by a discrete numerical method. The large dimensions of the grains allow an automatic image acquisition and data processing. From this photo graphic technique it is possible to extract incremental displacement and rotation fields (or the corresponding velocity fields) for the whole granular assembly.

The experiments have been performed in a plane shear box [START_REF] Joer | 11'2E : A new shear apparatus to study the behavior of granular media[END_REF] which allows to apply loading programs at the boundary of the slope. In this paper two cases are presented in a more detailed manner, even if the conclusions are also valid for the other cases not considered here. The first case is constituted by an active failure of the slope by pulling the retaining vertical wall outside of the granular assembly (see figure 12( a)), the second case is a slope "loaded at his top" , where avalanches are induced by a vertical loading on the flat top of the slope (see figure 12(b)). The configurations of the grain assemblies, as they are appearing on figure 12, are the real experimental ones. These initial configurations are exactly reproduced numerically in order to have a realistic granular assembly. Indeed the numerical experiments have shown the great influence of the initial configuration on the following deformations and failures. With this procedure (i.e. the simulation of experimental initial configurations), a good agreement has been found between the experimental displacement/rotation fields and the numerical ones [START_REF] Darve | Discrete instabilities in granular materials[END_REF]. This question is not discussed here, where it is chosen to focus on the mechanical analysis of the numerical results. The numerical discrete (1979), because the numerical algorithm for time integration is implicit and also because the interactions between grains are not described by springs, slides and dash pots. In the Contact Dynamics method, the grains are assumed perfectly rigid. Thus there is a kinematic constraint constituted by the non-penetrability condition (the so called Signorini condition, see figure 13(a)). The granular interactions are described by a simple coulombian friction, as represented on figure 13(b). If J.L is the friction coefficient between grains, a tangential displacement U t is allowed only if the tangential component R t of the interaction force is equal to : ± J.L Rn, where Rn is the normal component (figure 13(b)). For these simulations the friction coefficient between grains has been estimated to 0.5317 and between grains and walls to 0.25. Note that both these coefficients are the only mechanical parameters introduced in the computation. The other data have a geometric nature and are constituted by the initial configuration of the grain assembly, as discussed before. , where Rn and R t are respectively the normal and tangential components of the interaction force and J.L the friction coefficient. J is the "distance" (roughly) between two grains and U t the tangential inter-granular displacement.

Grain avalanches by active failure of the slope

In this case the number of grains is equal to : N = 335, the displacement speed (see figure 12(a)) to w-2 m/s (the physical time has however not any influence on the com putation). There are 1080 increments of computation (the time increment is: 5.610-3s).

Figure 14 presents a comparison between the measured displacement fields (on the left) and the computed values (on the right) for 4 values of the plate displacement. A first conclusion is obvious : the utilized discrete element method is able to capture the main characteristics of the successive grain avalanches thanks to the taking into account of the initial configuration of the grain assembly in the experiments. A second aspect is the fact that various failure modes are visible on figure 14. Some modes are deep and seem to be localized because of parallel displacement vectors (as for rigid body motions), while others are surficial and correspond to classical grain avalanches. These different failure modes (local or more global) will be characterized later on by bursts of kinetic energy of various amplitudes. Figure 15 shows a spatial comparison between the negative values of discrete second or der work and the kinetic energy of the grains. The dimensions of the square symbols are proportional to the intensity of the scalar variables d 2 W and e c . The main failure mode is surficial with erratic fields of kinetic energy confirming the essential features of grain avalanches. However the most important point to be noticed is the spatial correlation between both the area where grains have negative second order work and where grains have bursts of kinetic energy. It is worthy to note also that the correlation is not attached to grains (a grain with or without negative second order work may or not move) but to area which can be characterized as "unstable" . This correlation is also a temporal one. Let us consider now the global quantities D 2 W and E c determined by adding the local values of each grain. Figure 16 presents such a comparison. Each point on both graphs corresponds to one computation, which means that there are exactly 1080 points. All the values of D 2 W have been plotted : the positive as the negative ones. Some conclusions are obvious from these results like the fact that the global second order work is varying suddenly from positive values to negative ones. Moreover, large positive val ues are associated to strong negative ones (the diagram is more or less symmetric with respect to the axis D 2 W = 0). Indeed, the second order work is linked to energy, which can not be extracted in a finite manner from the medium. As regards the global kinetic energy variations, figure 16 shows successive 5 large failures for numbers of increments equal to 0, 90, 310, 490 and 1040. Smaller events (i.e. more surficial grain avalanches) occur for increment numbers of 200, 520, 710 and 760. Now let us go to the global second order work variations on same figure. We find 5 large minima for the negative values of D 2 W and also 4 smaller relative minima for respectively the same 9 increment numbers. For example the small grain avalanches at increment 710 and 760 are clearly noticeable on D 2 W graph by relative minima. This strong correlation is the manifestation of the deep link between failure and second order work. Another example is presented in the next sub section.

4.4

Grain avalanches by loading at the top Here the number of grains is equal to : N = 377. The vertical displacement rate of the plate is : 10-3 m/s (see figure 12(b)). The loading programme is still divided into 1000 increments. The friction coefficients are the same as before.

Figure 17 shows the incremental displacements fields for 8 values of the vertical displace ment of the plate. Two different failure modes appear here. For plate displacements equal to 2.1, 25.2, 51.1 mm it seems that there is a deep localized failure mechanism exhibiting a rotation of a rather solid body. On the contrary, for plate displacements equal to 5.6, 23. 1, 31.5, 35 mm the failure mechanism is rather a grain avalanche with a chaotic displacement field for the involved grains. From this rough analysis it appears difficult to conclude for the last case (plate displacement equal to 57.4 mm). On figure 18 the spatial comparison between the area where the local second order work takes negative values and the one where the local kinetic energy is finite is presented. As before (see figure 15), a correlation between both areas seems to exist.

Figure 19 shows the temporal correlation between the bursts of global kinetic energy and the minima of negative values for the global second order work. From the kinetic energy variations there are 5 visible events for plate displacements of about 7, 14, 32, 46 and 65 mm. 5 events are also clearly distinguishable on the second order work diagram for 5, 14, 31, 46 and 65 mm.

Conclusion

From an experimental point of view first, it has been shown that the question of fail me in geomaterials has to be revisited. The peak of q in undrained conditions for a loose sand does satisfy neither the Mohr-Coulomb's plastic limit condition nor a strain localization criterion, but the failure is predicted by Hill's condition. Secondly this kind of failure states (located strictly inside the plastic limit surface and a localization condition) can be expected from the theory of non-associated hardening elasto-plasticity. Thirdly the nu merical analysis carried out with an incrementally piece-wise linear and an incrementally non-linear constitutive relation have shown the existence of a whole stress domain where Hill's condition of stability is lost. In this domain, failure is characterized by bifurcation phenomena, losses of uniqueness and instabilities in Lyapunov's sense. Various modes d2WT max =0.001 9282 J. r:fwT min =-0.0032095 J.

,: _ � -: _I�-�--: �_c: --1�-�� -: t_ � -: _ I� : . -: : t . . � �l• : _ � -: r•: _ � -: -I: : : : : : : : : : : : : :

- : : : : : :

: : : : :

- 4 o �-----1 � 00 -----2 � 00 -----3� 00 -----4 � 00 -----5 �0 -0 ----6� 0-0 ----7�0 -0----8�0-0----9�0-0--�1 000
Increment EC max =31 .2689 J. EC min =0 J. : : : : : : : 2 0 ............... : ................ ?""'"'"""""""'"''?''''''"'""""""''!"'''''""'""""""'!""""'''''''''""'!""'""""""""""""'!'""'''''"'''''''!''""""""''''""""'!""""''''''""""' of failure can develop (localized, diffuse, due to geometric instabilities ... ). It seems that diffuse modes of failure can be analyzed by Hill's condition of stability. The equation of the boundaries of the failure domain has been given and discussed, as the equations of the cones of unstable stress directions. Diffuse modes of failure can appear only if the stress state lies inside the failure domain, if the loading direction is included inside the cone of unstable directions and if proper loading modes (stress controlled) are considered. For boundary value problems (cases not considered in this paper) , the global second order work (i.e. integrated on the whole body) and the symmetric part of the matrix of rigidity have to be analyzed.

The specific case of q-constant loading paths has been considered. It has been shown that the unexplained diffuse modes, exhibited in triaxial experiments on a loose sans along these paths, can be predicted by Hill's condition of stability.

In a second part, the phenomenon (of discrete nature) of grain avalanches has been tack led. The grain motions have been characterized by local and global kinetic energies. Surprisingly strong correlations between the bursts of kinetic energy and the negative values of second order work have been derived. As a general conclusion, it appears that the sign of second order work (under both its continuous and its discrete forms) could constitute a proper criterion to detect certain modes of failure, which appear before the localized failure but for specific loading condi tions given in the paper in the framework of continuum mechanics.

The application to natural hazards driven by gravity (landslides, debris flows, rockfalls, snow avalanches) is of great interest since some of them can not be explained in a classical plasticity framework. relations (symbols "o" and "+" for the incrementally non-linear and the octolinear models respectively). Figure 6 Lowest unstable stress levels in plane strain conditions for the loose and dense sand, for 2 different constitutive relations. Figure 9 Simulation of q-constant loading paths for loose Hostun sand by the octo linear constitutive model. Points are replaced by other symbols when the second order work takes negative values from the maximum of volume variations. The initial isotropic pressure is equal to 500 kPa. Figure 10 Simulation of q-constant loading paths for dense Hostun sand by the octo linear constitutive model. The second order work is never vanishing. Figure 11 Simulation of q-constant loading paths for loose Hostun sand by the octo linear and the non-linear constitutive model. Points are placed by other symbols when the second order work takes negative values from the maximum of volume variations. The initial isotropic pressure is equal to 100 kPa. Figure 12 Configurations of the slope and loading programmes to induce granular avalanches in a hi-dimensional medium. Figure 13 Illustrations of the non-penetrability condition 13(a) and of the coulombian friction 13(b), where Rn and R t are respectively the normal and tangential com ponents of the interaction force and f-l the friction coefficient. 8 is the "distance" (roughly) between two grains and U t the tangential inter-granular displacement .

  (d ;; , �) linked by the constitutive relation, the second order work is strictly positive : (2 .1) For hyperelastic-plastic constitutive relations, Drucker's postulate ( d ;; � P > 0, see Drucker D.C. (1959)) implies Hill's condition, since : d ;; � e 2 0 in hyperelasticity, and : dCJ de = dCJ dce + dCJ dE:P, with the usual decomposition of the incremental strain into t"Vrv rvrv an elastic part and a plastic part in classical elasto-plasticity. The link between Lyapunov's definition and Hill's condition has not been established in a general manner. The equivalence has been proved by Koi:ter W.T. (1969) in the case of elasticity.

Figure 1 .

 1 Figure 1. Example of an unstable state related to a diffuse mode of failure strictly inside Mohr-Coulomb plastic limit condition : typical undrained triaxial behavior of a loose sand.

  Analytical explicit computations are feasible with our incrementally piecewise linear model,called "the octo-linear model" (see chapter 1 by Darve and Servant and[START_REF] Darve | Incremental law for sands and clays. Simulations of mono tonic and cyclic tests[END_REF]) because it is constituted by 8 linear relations between � and d ,::

  = dadc: = t daMda = t daM 8 da, with : de= Mda and M 8 symmetric part of M.

  rv

�

  of �s are varying in a continuous manner from strictly positive initial values, what is always verified for realistic constitutive matrices). For an incrementally piecewise linear relation (as considered here with the octo-linear model), we have 4 different constitutive matrices Q following the signs of da1 and da3.

  (1.3)and (1.4)), where : 0 0 defines the plastic limit stress state, det f! = 0 is the plasticity criterion, and : Vl+ dE1 + dE3 = 0 is the flow rule. What is remarkable in equation (2.19) is the fact that it is a mixed condition, even if essentially this is the generalization of the notion of flow rule. All these results have been verified numerically with the octo-linear model and with the thoroughly incrementally non-linear relation (for a presentation of this relation see chapter 1 byDarve and Servant and Darve F. et al. (1995)).

Figure 2 .

 2 Figure 2. Proportional strain paths simulated by the incrementally non-linear model in the case of Hostun dense sand. ( a1 -1f) possesses a maximum for low R values, RE {0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0}.

  .. 4 • ••• •:::. •• ••. .,...-:. .. ...... . . ' .... / .. :..cC.';-

Figure 3 .

 3 Figure 3. Cones of unstable stress directions and boundary of the unstable domain in Rendulic plane for dense Hostun sand in the cases of the octo-linear model (on the left) and of the non-linear one (on the right).

  ... � --•• • ---�- ... �� • :: • � :c!\ "�� -, .;"" ' •:•;..;. �-�" • • • : 7�1 :__, . ! ... . .. . � ... ,.)/. ��--r �

Figure 4 .

 4 Figure 4. Cones of unstable stress directions and boundary of the unstable domain in Rendulic plane for loose Hostun sand in the cases of the octo-linear model (on the left) and of the non-linear one (on the right).

Figure 5 .

 5 Figure 5. Domains of instability in Rendulic Plane and axisymmetric conditions for the loose sand on the left , for the dense sand on the right and for 2 different constitutive relations (symbols "o" and "+" for the incrementally non-linear and the octolinear models respectively).

  2=t 7};�(7};)7};

Figure 6 .Figure 7 .

 67 Figure 6. Lowest unstable stress levels in plane strain conditions for the loose and dense sand, for 2 different constitutive relations.

  Note that equations (3.4) and (3.5) are exactly equivalent to equations (2.6), what means that bifurcation criterion and rupture rule are the same for q-constant loading path and for undrained loading. Thus the results of the subsection 2.2 apply with constant R equal to 1.

Figure 8 .

 8 Figure 8. Simulation of q-constant loading paths for loose Hostun sand by the octo linear constitutive model. Points are replaced by other symbols when the second order work takes negative values from the maximum of volume variations. The initial isotropic pressure is equal to 100 kPa.

Figure 9 .Figure 10 .Figure 11 .

 91011 Figure 9. Simulation of q-constant loading paths for loose Hostun sand by the octo linear constitutive modeL Points are replaced by other symbols when the second order work takes negative values from the maximum of volume variations. The initial isotropic pressure is equal to 500 kPa.

Figure 12 .

 12 Figure 12. Configurations of the slope and loading programmes to induce granular avalanches in a two-dimensional medium.

  essentially different from the one proposed by Cundall P.A. & Strack O.D.L.

Figure 13 .

 13 Figure 13. Illustrations of the non-penetrability condition 13(a) and of the coulombian friction 13(b), where Rn and R t are respectively the normal and tangential components of

Figure 14 .Figure 16 .

 1416 Figure 14. Experimental (on the left) and numerical (on the right) displacements fields for various plate displacements (see figure 12(a))

Figure 18 .

 18 Figure 18. Spatial comparison between the negative values of second order work d2W (on the left) and the values of kinetic energy ec (on the right).

  1'-.............. L . ............. L ............... L ............... � ................ l ................ L .............. L. ............. J . ................ L. . . . . . . . . . . .. ••oooooooooooooloooooooooooooooo oooooooooonooo'1oooooooooooooooo-:-oooooooooooooooot"••ooooooooooooool•ooooooooooooooool••ooooooooooooo

  ---r-----r-----r-----r-----r-----r-----r-----r-----r-----, ............. � ................ � ................ � . . . . . . .......... � . . .............. � ................ � ................ : . ................ : ................. : . ..............

Figure 19 .

 19 Figure 19. Comparison between global second order work and global kinetic energy as functions of the vertical plate displacement (or equivalently the number of time incre ments).

Figure 1

 1 Figure1Example of an unstable state related to a diffuse mode of failure strictly inside Mohr-Coulomb plastic limit condition : typical undrained triaxial behavior of a loose sand. Figure2Proportional strain paths simulated by the incrementally non-linear model in

Figure 2

 2 Figure1Example of an unstable state related to a diffuse mode of failure strictly inside Mohr-Coulomb plastic limit condition : typical undrained triaxial behavior of a loose sand. Figure2Proportional strain paths simulated by the incrementally non-linear model in the case of Hostun dense sand. ( 0"1 -�) possesses a maximum for low R values, RE {0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0} .

Figure 3

 3 Figure 3 Cones of unstable stress directions and boundary of the unstable domain in Rendulic plane for dense Hostun sand in the cases of the octo-linear model (on the left) and of the non-linear one (on the right).

Figure 4

 4 Figure 4 Cones of unstable stress directions and boundary of the unstable domain in Rendulic plane for loose Hostun sand in the cases of the octo-linear model (on the left) and of the non-linear one (on the right).

Figure 5

 5 Figure 5 Domains of instability in Rendulic Plane and axisymmetric conditions for the loose sand on the left , for the dense sand on the right and for 2 different constitutive

Figure 7

 7 Domain of instability in deviatoric plane for the loose and dense sand, and for 2 different constitutive relations.

Figure 8

 8 Simulation of q-constant loading paths for loose Hostun sand by the octo linear constitutive model. Points are replaced by other symbols when the second order work takes negative values from the maximum of volume variations. The initial isotropic pressure is equal to 100 kPa.

Figure 14

 14 Figure 14 Field of numerical incremental displacements at different values of rotation angle a.

Figure 15

 15 Figure 15 Spatial comparison between negative values of second order work d 2 W (on the left) and the values of kinetic energy ec (on the right).

Figure 16

 16 Figure 16 Comparison between global second order work and global kinetic energy as functions of rotation angle a. Figure 17 Field of incremental displacements at different values of vertical plate dis placement u.

Figure 18

 18 Figure 18 Spatial comparison between the negative values of second order work d 2 W (on the left) and the values of kinetic energy ec (on the right) .

Figure 19

 19 Figure19Comparison between global second order work and global kinetic energy as functions of the vertical plate displacement .

  ............... i ................ t . . . ............. j ................. j ................ i . . .............. + ................ t . ............... j ................. j ............. ..

	5: st' 1 �	: �	: �	: �	: i	: i	: i	: i	: i	: i

Acknowledgements

The supports of European projects DIGA (5th PCRD), LESSLOSS (6th PCRD) and french national projects PIR (RGCU) and ACI CatNat (2002) are gratefully acknowl edged.

da3 = negative constant .

The proper conjugate variables to analyse this path are : cv -a3 and q -c1, where : cv = c1 + 2c3, because we have :

The second order work can thus be written as follows : d 2 W = dcvda3 + dqdc1, and with the constraint : dq = 0, it comes :

(3.2)

If cv passes through a minimum value, this minimum will be an unstable state according to Hill's condition. On the contrary a3 is monotonously decreasing.

As the relative volume variations �� is equal to ( -cv) : �� = -cv, the minimum for cv is a maximum for ��.

Thus a "small" additional positive volume variation (for a volume variation controlled loading) induces a sudden failure of the sample. It can be concluded that this maximum of volume variations is unstable following Lyapunov's definition of stability. Now if the constitutive relation is written under the following form (for axisymmetric conditions) : the analysis can be carried out as previously.

Along the path the constraint : dq = 0 is fulfilled and, at cv minimum, dcv is zero.

Thus, the bifurcation criterion is given by :