List of authors: T. Guignon 1 M. Hacene 2 J.M. Gratien 3 Porous media flow simulation in petroleum industry is used in reservoir simulation or basin modeling to plan fields development or characterize potential new fields. It's also used in CO2 geological sequestration to prove feasibility and stock safety for long time period. One of the main computation done in theses simulations is solving one or more large sparse unstructured linear systems at each time step. These systems come from finite volume discretization of Darcy flow equation. In reservoir simulation time spent in linear solver can be up to 80% of the total simulation time so solver efficiency is a key point in simulation performance.

The classical method for solving theses systems is an iterative Krylov method like BiCGstab. Since most systems are not well conditioned we use some preconditioner like ILU0, ILU(k) or CPR-AMG for difficult problems which are numerically efficient but not really scalable on distributed memory parallel computers.

Using hardware accelerators like GPU is another way to improve linear solver performance because theses accelerators provide a very high peak floating point and memory performances compared to CPU. But theses improvements comes with some constraints on parallel algorithms we want to deploy on theses accelerators. Mainly datas and computations must have a very regular structure and algorithms must exhibit a very fine grain parallelism to achieve high performance computation. These constraints are very similar to those we have with data parallel programing models and data parallel architectures that where used in the 90's.

Iterative Krylov method like BiCGStab rely mainly on two main software components: a sparse matrix vector product and preconditioner application ,for example solving L.U.x = y for ILU. So the focus is on deploying these components to GPU. Other solver parts, simple vector algebra, can be easily done with a dense BLAS GPU library like CUBLAS.

We will first discuss sparse matrix vector product for unstructured sparse matrix which is now well understood. We present on figure some spmv 4 results with different GPU systems compared to CPU. Linear systems come from reservoir simulation study cases and for most of them there at least is a 10x increase in performance over the fastest CPU tested.

In second we will discuss preconditioning. We will focus on preconditioner application (building is still done on CPU). We will first study L.U.x = y on GPU where L and U are linear system incomplete factors. Unstructured sparse triangular solve required in L.U.x = y does not exhibit much parallelism with the linear system "natural ordering". Thus for an efficient triangular solve on GPU we must reorder system to make it with large diagonal blocks and extra diagonal sparse blocks A u and A l :

P.A.P -1 =      D 1 A u D 2 . . . A l D n     
This reordering is obtained by coloring vertices of the adjacency graph such as a vertex (an equation) has no neighbour of his own color and group vertices of the same color. After ILU solving U.x = y (or L.x = y, with unit diagonal) expose a high level of fine grain parallelism since we have to chain diagonal inversion and spmv:

     D 1 . . . U 1 . . . D 2 U 2 . . . . . . D n      .      x 1 x 2 . . . x n      =     

1 1 Figure 1 :

 11 Figure 1: Comparing CPU spmv, IFP spmv and NVIDIA spmv on various GPU configurations.

Table 1 :

 1 Comparing various CPU and GPU preconditioner on full CO2 sequestration simulation, 114000 grid cells

	y 1 y 2 . . . y n	    


Jacques Giacomoni -Habib

Mâagli -Paul Sauvy

EXISTENCE OF COMPACT SUPPORT SOLUTIONS FOR A QUASILINEAR AND SINGULAR PROBLEM

Let Ω be a C 2 bounded domain of R N , N ≥ 2. We consider the following quasilinear elliptic problem:

where p > 1 and ∆ p u def = div `|∇u| p-2 ∇u ´denotes the p-Laplacian operator. In this paper, λ > 0 is a real parameter, the exponents q and r satisfy -1 < r < q < p -1 and K : Ω -→ R is a positive function having a singular behaviour near the boundary ∂Ω. Precisely, K(x) = δ(x) -k L(δ(x)) in Ω, with 0 < k < p, L a positive perturbation function and δ(x) the distance of x ∈ Ω to ∂Ω.

By using a sub-and super-solution technique, we discuss the existence of positive solutions or compact support solutions of (P λ ) in respect to the blow-up rate k. Precisely, we prove that if k < 1 + r, (P λ ) has at least one positive solution for λ > 0 large enough, whereas it has only compact support solutions if k ≥ 1 + r.

Introduction

Let Ω be a C 2 bounded domain of R N , N ≥ 2. We discuss the existence of weak solutions in W 1,p 0 (Ω) ∩ L ∞ (Ω) to (P λ ) -∆ p u = K(x)(λu q -u r ) in Ω, u = 0 on ∂Ω, u ≥ 0 in Ω.

u ∈ W 1,p 0 (Ω) ∩ L ∞ (Ω) is a weak solution to (P λ ) if for all test functions ϕ ∈ D(Ω), Ω |∇u| p-2 ∇u.∇ϕ dx = Ω K(x) (λu q -u r ) ϕ dx.

(1. [START_REF] Allegretto | A Picone's identity for the p-Laplacian and applications[END_REF] In the equation in (P λ ), λ > 0 is a positive parameter, -1 < r < q < p -1 and K ∈ C(Ω) is a positive function having a singular behaviour near the boundary ∂Ω. Precisely, K(x) = δ(x) -k L(δ(x)) in Ω, with 0 < k < p and L ∈ C 2 ((0, 2d]) a positive function, with d def = diam(Ω), defined as follows:

L(t) = exp 2d t z(s) s ds , (1.2)
1 with z ∈ C([0, 2d]) ∩ C 1 ((0, 2d]) and z(0) = 0. Let us note that (1.2) implies that lim t→0 + tL (t) L(t) = 0 (1.3) and for all ε > 0, lim t→0 + t ε L(t) = 0 (1.4) and lim t→0 + t -ε L(t) = +∞.

(1.5)

The above asymptotics of L force ∀ξ > 0 , lim

t→0 + L(ξt) L(t) = 1.
Then L belongs to the Karamata class [START_REF] Karamata | Über die Hardy-Littlewoodsche Umkehrung des Abelschen Stätigkeitssatzes[END_REF]. Setting K the class of functions satisfying (1.2), we get the following properties: if L 1 , L 2 ∈ K and if p ∈ R, then L 1 .L 2 ∈ K and L 1 p ∈ K. In the present paper, we investigate first the following issues for the problem (P λ ): existence of non-trivial weak solutions according to λ > 0, Hölder regularity of weak solutions. Next, we study further the properties of non-trivial solutions. Since the non-linearity in the right-hand side is a singular absorption term near the boundary, a non-trivial weak solution may not be positive everywhere in Ω and compact support (non-trivial) weak solutions or compactons (solutions with zero normal derivative at the boundary) may exist for stronger singularities, that is for large k > 0 whereas for small k > 0 any non-trivial weak solution is positive. Then, the natural question is to determine the borderline condition on the parameter k, which gives the strength of the singular potential K, between existence of positive weak solutions and existence of free boundary weak solutions. The existence of compact support solutions is important in the applications, in particular in biology models (population dynamics and epidemiology models for instance) and was investigated quite intensely for nonlinear reaction diffusion equations with absorption in the last decades. In particular, concerning the case where the equation involves a quasilinear and degenerate operator, we can refer to the result in Vázquez [START_REF] Vásquez | A strong maximum principle for some quasilinear elliptic equations[END_REF] where under a suitable condition about the behaviour of the non-linearity near the origin, a strong maximum principle is proved and consequently the positivity of solutions. The given condition is sharp in the sense that for different situations where this condition is not satisfied, the existence of free boundary solutions is shown. In Alvarez-Díaz [START_REF] Alvarez | On the behaviour of the free boundary of some nonhomogeneous elliptic problems[END_REF] (see also Díaz [4] for related results on the subject), the authors consider a class of nonhomogeneous reaction-diffusion equations with strong absorption and study the behaviour of the solution near the free boundary. In particular, a non degeneracy property (the solution grows faster than some function of the distance to the free boundary) is obtained when the growth of the reaction term near the boundary satisfies some estimate by below. In Il'yasov-Egorov [START_REF] Il | Hopf boundary maximum principle violation for semilinear elliptic equations[END_REF], the authors consider a semilinear equation with a similar (and non singular) conflicting nonlinearity as in the equation in (P λ ) and the existence of compactons is proved using the fibering method. An interesting feature of this result is that the Hopf lemma is violated for such kind of equations. In the present work, we consider the extension case where the equation involves a p-Laplace operator and a singular potential in the right-hand side and show that a more complex situation occurs in respect to the non singular case.

In the next section, we give the main results proved in this paper. These results extend a previous work due to Haitao [7] in the semilinear case (p = 2) and which involves a smaller class of nonlinearities.

Main results

The main results of our paper concerning the problem (P λ ) are stated below:

Theorem 2.1. When k < 1 + r, there exists a constant Λ 1 > 0 such that:

(1) For λ > Λ 1 , (P λ ) admits a positive weak solution.

(2) Any weak solution of (P λ ) is C 1,β Ω for some β ∈ (0, 1).

(3) For λ < Λ 1 , (P λ ) has no positive solution.

Theorem 2.2. Let r > 0 and one of the two following conditions be satisfied:

1 + r > q and k ∈ » 1 + r, 1 + (p -1)(r + 1) p -q + r « , (2.1) 1 + r ≥ q and k ∈ [1 + r, 2 + r). (2.2)
Then, there exists Λ 2 > 0 such that:

(1) For λ > Λ 2 , (P λ ) has a compact support weak solution u λ .

(2) Any weak solution of (P λ ) is C 1,β Ω for some β ∈ (0, 1).

(3) For λ < Λ 2 , (P λ ) has no non-trivial solution.

The outline of the paper is as follows. Before giving the proofs of those theorems, we establish some useful preliminary results in the next section. The proof of Theorem 2.1 is given in section 4 and the proof of Theorem 2.2 is given in section 5. The technical results stated in section 3 are finally proved in appendix A and B. The related regularity results are a consequense of the general regularity results stated in appendix C.

Throughout this paper, we will use the following notations:

(1) To p ∈ (1, +∞) we associate

p def = p p-1 . (2) For x ∈ Ω, δ(x) def = dist(x, Ω) = inf y∈Ω d(x, y). (3) d def = diam(Ω) = sup x,y∈Ω d(x, y).
(4) Let ω be a non-empty set of Ω and f, g : ω -→ [0, +∞]. We write

f (x) ∼ g(x) in ω if there exist two positive constants C 1 and C 2 such that ∀x ∈ ω, C 1 f (x) ≤ g(x) ≤ C 2 f (x). (5) Let ω ⊂ R N , L N (ω) denotes the N -dimentionnal Lebegue's mesure of ω. (6) Let ε > 0, we define Ω ε def = {x ∈ Ω, δ(x) < ε}. (7) ν : ∂Ω -→ R N denotes the outward normal associated to Ω. (8) For v ∈ W 1,p 0 (Ω), we write v def = ∇v L p (Ω) = Ω |∇v| p dx 1 p . ( 9
)
The function ϕ 1 ∈ W 1,p 0 (Ω) denotes the positive and renormalized (i.e. ϕ 1 L p (Ω) = 1) eigenfunction corresponding to the first eigenvalue of -∆ p ,

λ 1 def = inf Ω |∇v| p dx, v ∈ W 1,p 0 (Ω), Ω |v| p dx = 1 .
It is a weak solution of the following eigenvalue problem:

-∆ p u = λ 1 u p-1 in Ω, u = 0 on ∂Ω, u ≥ 0 in Ω.
Using Moser iterations and the regularity result in Lieberman [START_REF] Lieberman | Boundary regularity for solutions of degenerate elliptic equations[END_REF], ϕ 1 ∈ C 1,α Ω for some α ∈ (0, 1). Moreover the strong maximum and boundary principles from Vásquez [START_REF] Vásquez | A strong maximum principle for some quasilinear elliptic equations[END_REF], guarantee that ϕ 1 satisfies those two properties:

(a) There exist two positive constants K 1 and K 2 only depending on p and Ω such that:

∀x ∈ Ω, K 1 δ(x) ≤ ϕ 1 (x) ≤ K 2 δ(x). (2.

3)

(b) There exist ε * > 0 and δ * > 0 only depending on p and Ω such that:

∀x ∈ Ω δ * , |∇ϕ 1 (x)| > ε * .
(2.4)

Preliminary results

3.1.

A non-existence lemma.

Lemma 3.1. When k < 1 + r, there exists λ * > 0 such that (P λ ) has no non-trivial solution for λ ≤ λ * .

Proof: Let us define

λ1,K def = inf v∈W 1,p 0 (Ω) v =0 Z Ω |∇v| p dx Z Ω K(x)|v| p dx .
From the Hardy's inequality, there exists a constant C > 0 only depending on Ω and p such that for all v ∈ W 1,p 0 (Ω),

Ω |v| p δ(x) p dx ≤ C Ω |∇v| p dx. Since k < p, λ 1,K > 0. Let u ∈ W 1,p 0 (Ω) ∩ L ∞ (Ω)
be a non-trivial solution of (P λ ), then from (1.1) taking u ∈ W 1,p 0 (Ω) as a test function we get,

0 < λ 1,K Ω K(x)u p dx ≤ Ω |∇u| p dx = Ω K(x) λu q+1 -u r+1 dx. (3.1)
This inequality is impossible for

λ ≤ λ * def = min{1, λ 1,K }. Indeed, -if u(x) ≤ 1, λu q+1 -u r+1 ≤ 0 as soon as λ ≤ 1, -if u(x) > 1, K(x) λu q+1 -u r+1 < λ 1,K K(x)u q+1 as soon as λ ≤ λ 1,K . Then, either L N ({x ∈ Ω, u(x) > 1}) = 0 and we get 0 < λ 1,K Ω K(x)u p dx ≤ 0, or λ 1,K Ω 1 {u>1} K(x)u p dx ≤ λ 1,K Ω 1 {u>1} K(x)u q+1 dx,
which contradicts q < p -1 .

3.2.

Construction of a sub-solution for (P λ ). Lemma 3.2. When k < 1 + r, there exist M > 0, λ * > 0 and τ > 1 such that u λ def = M ϕ 1 τ is a sub-solution of (P λ ) in Ω, provided that λ ≥ λ * .

Proof: Let M > 0 and τ > 1, then we define

u λ = M ϕ 1 τ in Ω. A straightforward computation yields -∆ p u λ = -(M τ ) p-1 (p -1)(τ -1)|∇ϕ 1 | p ϕ 1 (τ -1)(p-1)-1 -λ 1 ϕ 1 τ (p-1)
and

K(x) (λu λ q -u λ r ) = -K(x) (M r ϕ 1 τ r -λM q ϕ 1 τ q ) .
By properties (2.3) and (2.4) of the function ϕ 1 , if we let

δ 0 def = min ( δ * , ε * K2 " (τ -1)(p -1) 2λ1 « 1 p , 1 K2 " 1 2λM q-r « 1 τ (q-r)
) , both of the above expressions are negative on Ω δ 0 . Moreover,

∆ p u λ (x) ∼ (M τ ) p-1 (τ -1)δ(x) (τ -1)(p-1)-1 in Ω δ 0 and K(x) (u λ r -λu λ q ) ∼ M r L(δ(x))δ(x) τ r-k in Ω δ 0 .
Since k < 1 + r, we can choose a constant τ > 1 satisfying

(τ -1)(p -1) -1 < τ r -k.
Hence, for M > 0 large enough we get -∆ p u λ ≤ K(x) λu q λ -u r λ in Ω δ 0 . In Ω \ Ω δ 0 , K, ϕ 1 and |∇ϕ 1 | are bounded, therefore there exists λ * > 0 such that for λ ≥ λ * , -∆ p u λ ≤ K(x) λu q λ -u r λ in Ω \ Ω δ 0 . Thus, u λ is a sub-solution of (P λ ) in Ω for M large enough and λ ≥ λ * .

Construction of a super-solution for (P λ ).

We consider the following problem:

(Q) -∆ p v = K(x)v q in Ω, v = 0 on ∂Ω, v > 0 in Ω,
with q, p and K satisfying the above assumptions.

Lemma 3.3.

(

) If k ∈ (0, 1 + q), (Q) has a unique solution v ∈ W 1,p 0 (Ω) ∩ C Ω satisfying v(x) ∼ δ(x) in Ω. (2) If k = 1+q, (Q) has a unique solution v ∈ W 1,p 0 (Ω)∩C Ω satisfying v(x) ∼ δ(x) Z 2d δ(x) L(t) t dt ! 1 p-k 1
in Ω.

(

) If k ∈ 1 + q, 1 + q + p-(1+q) p , (Q) has a unique solution v ∈ W 1,p 0 (Ω)∩ C Ω satisfying v(x) ∼ δ(x) p-k p-(1+q) L (δ(x)) 1 p-(1+q) in Ω. (4) If k ∈ 1 + q + p-(1+q) p , p , (Q) has a unique solution v ∈ W 1,p loc (Ω) ∩ C 0 Ω satisfying v(x) ∼ δ(x) p-k p-(1+q) L (δ(x)) 1 p-(1+q) in Ω. 3
(5) If k = p and if L satisfies the following condition:

2d 0 t -1 L(t) 1 p-1 dt < +∞, (3.2) (Q) has a unique solution v ∈ W 1,p loc (Ω) ∩ C 0 Ω satisfying v(x) ∼ Z δ(x) 0 t -1 L(t) 1 p-1 dt ! p-1 p-(1+q)
in Ω.

Proof: See section A in appendix.

From a solution of (Q), we can easily construct a super-solution of (P λ ). Indeed, let us consider v ∈ W 1,p loc (Ω) ∩ C 0 Ω the solution of (Q) given by lemma 3.3. Then, u λ def = M v is a super-solution of (P λ ) in Ω as soon as M ≥ λ 1 p- (1+q) . Particularly when k < 1 + r and λ ≥ λ * , chosing M sufficiently large u λ ∈ W 1,p 0 (Ω) ∩ C Ω and is a super-solution of (P λ ) in Ω satisfying

u λ ≤ u λ and u λ (x) ∼ δ(x) in Ω.
Now let us state a non-existence result for the problem (Q).

Proposition 3.1. Let v ∈ W 1,p 0 (Ω) ∩ C Ω be a positive sub-solution of (Q) in Ω and let us suppose there exists ε > 0 such that Ω K(x)ϕ 1 p-1+ε dx = +∞. (3.3) Therefore, for any η > 0, (Q) has no weak solution v ∈ W 1,p loc (Ω) ∩ C 0 Ω such that v ≥ ηv in Ω. Proof: See section B in appendix. Corollary 3.1. If k > p, there is no non-trivial weak solution of (Q). 4. Proof of Theorem 2.1 4.1. Existence of a C 1,β positive solution when λ ≥ λ * . Proposition 4.1. When k < 1 + r, provided λ ≥ λ * , (P λ ) has a weak solution u λ ∈ W 1,p 0 (Ω) ∩ L ∞ (Ω). Furthermore, there exists β ∈ (0, 1) such that u λ ∈ C 1,β Ω .
Proof: In the equation of (P λ ), the expression h λ (x, v) def = K(x)(λv q -v r ) involves a singular term K(x) which blows up as δ(x) → 0, so we can not directly apply the sub-and super-solution method on Ω. To overcome this difficulty, we apply a sub-and super-solution method in a sequense of

subdomains of Ω. Let us introduce (Ω k ) k∈N * ⊂ Ω an increasing sequence of smooth subdomains of Ω such that Ω k -→ k→∞ Ω in the Hausdorff topology with ∀k ∈ N * , 1 k + 1 < dist(∂Ω, ∂Ω k ) < 1 k .
Then, for all k ∈ N * we consider the following problem:

(P k ) -∆ p u k = K(x)(λu k q -u k r ) in Ω k , u k = u λ on ∂Ω k , u k ≥ 0 in Ω k . By definition of Ω k , there exists C k > 0 such that ∀v ∈ I k def = min Ω k u λ , max Ω k u λ , sup x∈Ω k ∂h λ ∂v (x, v) ≤ C k .
As a consequense, there exists µ k > 0 such that for all x ∈ Ω k , the func-

tion v -→ h λ (x, v) + µ k v p-1 is increasing on I k .
Therefore by sub-and super-solution method, (P k ) has a solution u k ∈ W 1,p (Ω k ). Indeed, we can construct the following iterative monotone scheme: for all n ∈ N * , let u k,n ∈ W 1,p (Ω k ) be the weak solution of

(P k,n ) -∆ p u k,n + µ k (u k,n ) p-1 = h λ (x, u k,n-1 ) + µ k (u k,n-1 ) p-1 in Ω k , u k,n = u λ on ∂Ω k , u k,n ≥ 0 in Ω k with u k,0 = u λ in Ω k . By induction on n ∈ N, (P k,n ) has a unique solution u k,n ∈ W 1,p (Ω k ). By using the weak comparison principle (u k,n ) n∈N satisfies u λ ≤ u k,n ≤ u k,n+1 ≤ u λ in Ω k . Consequently, for all n ∈ N * , h λ (x, u k,n-1 ) + µ k (u k,n-1 ) p-1 -(u k,n ) p-1 ∈ L ∞ (Ω k )
and since u λ is smooth in Ω, we can state by a regularity result due to Lieberman [START_REF] Lieberman | Boundary regularity for solutions of degenerate elliptic equations[END_REF]

(see Theorem 1) that (u k,n ) n∈N ⊂ C 1,γ Ω k for some γ ∈ (0, 1). Moreover there exists a constant C > 0 only depending on γ, Ω k , u λ L ∞ (Ω k ) and u λ L ∞ (Ω k ) such that u k,n C 1,γ (Ωk) ≤ C. From Ascoli- Arzelà theorem, there exist u k ∈ C 1 (Ω k ) and a subsequence (u k,m ) m∈N such that u k,m -→ m→∞ u k in C 1 Ω k . Passing to the limit when n → +∞ in (P k,n ),
u k is a weak solution of (P k ).

For all k ∈ N, we define ũk def = 1 Ω k .u k in order to extend u k on Ω by zero. We prove that (ũ k ) k∈N is an increasing sequence in Ω. Indeed, since Ω k ⊂ Ω k+1 , if we compare u k+1 with every term of (u k,n ) n∈N in Ω k , using the weak comparison principle we get

∀n ∈ N, u k,n ≤ ũk+1 in Ω k .
Passing to the limit in the above inequality, ũk (x)

k∈N is non-decreasing for all x ∈ Ω. Therefore there exists u λ ∈ L ∞ (Ω) such that ũk -→ k→∞ u λ a.e.
in Ω and

u λ ≤ u λ ≤ u λ in Ω. ( 4
(x) [λ(u λ ) q -(u λ ) r ] ∈ W -1,p (Ω), which im- plies that u λ ∈ W 1,p 0 (Ω).
* ≤ Λ 1 ≤ λ * < +∞. By definition of Λ 1 , for any λ > Λ 1 there exists µ ∈ (Λ 1 , λ) such that (P µ ) has a positive solution u µ ∈ W 1,p 0 (Ω) ∩ L ∞ (Ω). Moreover using proposition C.1, u µ ∈ C 1,β Ω . Since u µ is a sub-solution to (P λ ), we prove that u µ ≤ u λ in Ω. Indeed, K(x) > 0 in Ω, so there exists δ 0 > 0 such that -∆ p u µ ≤ 0 ≤ -∆ p (C 0 ϕ 1 ) in Ω δ 0 , with C 0 > 0 large enough to satisfy u µ ≤ C 0 ϕ 1 on ∂Ω δ 0 .
By the weak comparison principle, u µ ≤ C 0 ϕ 1 in Ω δ 0 . Morever u µ and ϕ 1 are bounded in Ω \ Ω δ 0 , thus u µ ≤ Cϕ 1 in Ω for some C > 0. Therefore chosing M sufficiently large in the definition of u λ , we get u µ ≤ u λ in Ω. Finally, applying again sub-and super-solution technique as in step 1, we get a solution u λ ∈ C 1,β Ω of (P λ ) .

Let k ∈ 1 + r, 1 + q + p-(1+q) p . Then, under condition Ω K(x) (u λ ) r+1 dx < +∞, (5.1)
there exists λ * * > 0 such that the problem (P λ ) has a non-trivial weak solution

u λ ∈ W 1,p 0 (Ω) ∩ L ∞ (Ω) as soon as λ > λ * * . Remark 5.1. Since u λ ∈ W 1,p 0 (Ω), by Hardy's inequality u λ δ(x) ∈ L p (Ω)
. So using Hölder's inequality, assumption (5.1) in Theorem 5.1 is satisfied if

L(δ(x))δ(x) r+1-k ∈ L α (Ω), where α = p r+1 > 1.
And this last condition is satisfied if

k < 1 + r + p -(1 + r) p . (5.2)
So (5.2) implies (5.1); but this condition is not sharp and can be weakened by using the precise behaviour of u λ given in lemma 3.3. Indeed,

(1) if k ∈ [1 + r, 1 + q), u λ (x) ∼ δ(x) in Ω. Therefore condition (5.1) is satisfied if k < 2 + r. (5.3) (2) if k = 1 + q, condition (5.1) is also satisfied if k < 2 + r.
(

) if k ∈ " 1 + q, 1 + q + p-(1+q) p " , u λ ∼ δ(x) p-k p-(1+q) " L (δ(x)) " 1 p-(1+q) in Ω. Therefore, condition (5.1) is satisfied if k < 1 + (p -1)(r + 1) p -q + r . (5.4) 3
Remark that if 1 + r > q, (5.3) is always true for k ∈ [1 + r, 1 + q] and since

1 + q < 1 + (p -1)(r + 1) p -q + r ⇐⇒ r + 1 > q, (5.5)
condition (2.1) implies (5.1). Similarly if 1 + r ≤ q, by equivalence (5.5), (5.4) is never satisfied for k ∈

" 1 + q, 1 + q + p-(1+q) p "
and condition (2.2) implies (5.1). We can easily check that both conditions (2.1) and (2.2) are weaker than (5.2). Moreover, let us suppose one of the following conditions be satisfied:

1 + r > q and k ∈ " 1 + (p -1)(r + 1) p -q + r , 1 + q + p -(1 + q) p « , 1 + r ≥ q and k ∈ " 2 + r, 1 + q + p -(1 + q) p « .
Then, using lemma 3.3 again, condition (5.1) is not satisfied, which guarantees the "sharpness" of conditions (2.1) and (2.2).

In the proof of Proposition 5.1, we will need the following well known lemma.

Lemma 5.1. Let x, y ∈ R N and • , • the standard scalar product in R N . Then there exists a constant C p > 0 such that

|x| p-2 x -|y| p-2 y, x -y ≥    C p |x -y| p if p ≥ 2, C p |x -y| 2 (|x| + |y|) 2-p if 1 < p < 2.
Proof: See Lemma 4.2 in Lindqvist [START_REF] Lindqvist | On the equation div (|∇u| p-2 ∇u) + λ|u| p-2 u = 0[END_REF] or Lemma A.0.5 in Peral [START_REF]Multiplicity of solutions for the p-laplacian[END_REF].

Proof of proposition 5.1: Let us introduce the functionnal

I λ (v) = 1 p Ω |∇v| p dx + 1 r + 1 Ω K(x)|v| r+1 dx - λ q + 1 Ω K(x)|v| q+1 dx,
with v ∈ W 1,p 0 (Ω). Let ϕ 0 = 0 ∈ D(Ω) be a non-negative function. Therefore there exists λ * * > 0 such that I λ (ϕ 0 ) < 0 for λ > λ * * . Let us fix a constant M > 1 such that M u λ ≥ ϕ 0 in Ω and introduce the cut-off function f λ defined in Ω × R by:

f λ (x, v) =    K(x) [λ(M u λ ) q -(M u λ ) r ] if v > M u λ (x), K(x) [λ|v| q -|v| r ] if v ∈ [0, M u λ (x)] , 0 if v < 0. The function v -→ f λ (x, v) is a Carathéodory function. For (x, v) ∈ Ω × R, let us set F λ (x, v) = v 0 f λ (x, t
) dt and consider the functional E λ defined as follows:

∀v ∈ W 1,p 0 (Ω), E λ (v) = 1 p Ω |∇v| p dx - Ω F λ (x, v(x)) dx.
A straightforward computation gives us

E λ (v) = 1 p Ω |∇v| p dx - λ q + 1 A(v, q) + 1 r + 1 A(v, r) -λB(v, q) + B(v, r) - r r + 1 C(r) + λ q q + 1 C(q), (5.6)
with

A(v, s) def = Ω 1 {0≤v≤M u λ } K(x)|v| s+1 dx, B(v, s) def = Ω 1 {v≥M u λ } K(x) (M u λ ) s vdx and C(s) def = Ω 1 {v≥M u λ } K(x)(M u λ ) s+1 dx.
Let ε > 0 and v ∈ W 1,p 0 (Ω), then we split the integral A(v, q) in Ω \ Ω ε and Ω ε :

A(v, q) = Ω\Ω ε 1 {0≤v≤M u λ } K(x)|v| q+1 dx + Ω ε 1 {0≤v≤M u λ } K(x)|v| q+1 dx def = A Ω\Ωε (v, q) + A Ω ε (v, q).
Since in Ω \ Ω ε , K is bounded, from the embedding W 1,p 0 (Ω) → L q+1 (Ω), there exists a constant C 1 such that

A Ω\Ω ε (v, q) ≤ C 1 v q+1 .
(5.7)

In Ω ε , by Hölder's inequality we have,

A Ωε (v, q) ≤ A Ωε (v, r) τ Ωε 1 {0≤v≤M u λ } K(x)|v| p dx 1-τ , with τ = p-(1+q) p-(1+r) < 1.
Using inequality (1.4) and Hardy's Inequality, we finally obtain, for ε small enough

A Ωε (v, q) ≤ C 2 ε 1 2 (p-k)(1-τ ) A Ωε (v, r) τ Ωε |v| p δ(x) p dx 1-τ ≤ C 2 ε 1 2 (p-k)(1-τ ) (τ A(v, r) + C 3 (1 -τ ) v p ) .
(5.8)

From the above arguments and since

B(v, q) = Z Ω\Ω ε 1 {v≥M u λ } K(x) (M u λ ) q vdx + Z Ω ε 1 {v≥M u λ } K(x) (M u λ ) q vdx def = B Ω\Ωε (v, q) + B Ω ε (v, q),
we also get

B Ω\Ωε (v, q) ≤ C 4 v (5.9)
and

B Ω ε (v, q) ≤ C 5 ε 1 2 (p-k)(1-τ ) (τ B(v, r) + C 6 (1 -τ ) v p ) .
(5.10)

Using inequalities (5.7) to (5.10),

E λ (v) ≥ 1 2p v p -λ C 1 q + 1 v q+1 -λC 4 v + 1 2 B(v, r) + 1 2(r + 1) A(v, r) - r r + 1 C(r) + λ q q + 1 C(q), (5.11)
for ε > 0 sufficiently small. With condition (5.1), this inequality proves that E λ is coercive and bounded from below on W 1,p 0 (Ω). So let us define

c λ def = inf v∈W 1,p 0 (Ω) E λ (v) and let (v n ) n∈N ⊂ W 1,p 0 (Ω) be a minimizing sequence of E λ , that is to say E λ (v n ) -→ n→+∞ c λ . E λ (v n ) n∈N is bounded, therefore by inequality (5.11) (v n ) n∈N is bounded in W 1,p
0 (Ω). Thus, there exist u λ ∈ W 1,p 0 (Ω) and a subsequence (v m ) m∈N such that v m -→ m→+∞ u λ weakly in W 1,p 0 (Ω), strongly in L q+1 (Ω) and in L 1 (Ω) and a.e. in Ω. Then we get

u λ ≤ lim inf m→+∞ v m .
(5.12) Using Fatou's Lemma and inequality (5.1), it follows that

1 r A(u λ , r) + B(u λ , r) ≤ lim inf m→+∞ 1 r A(v m , r) + B(v m , r) < +∞. ( 5.13)
Again from Fatou's lemma and inequalities (5.8),(5.10) and (5.12),

λ q + 1 A Ω ε (u λ , q) + λB Ω ε (u λ , q) ≤ lim inf m→+∞ " λ q + 1 A Ω ε (v m , q) + λB Ω ε (v m , q) « ≤ C 7 ε 1 2 (p-k)(1-τ ) . (5.14) Since v m -→ m→+∞ u λ in L q+1 (Ω) and in L 1 (Ω), A Ω\Ω ε (v m , q) -→ m→+∞ A Ω\Ω ε (u λ , q) and B Ω\Ω ε (v m , q) -→ m→+∞ B Ω\Ω ε (u λ , q).
(5.15)

Gathering the estimates (5.12) to (5.15) and using (5.6), we obtain:

c λ = lim inf m→+∞ E(v m ) ≥ E λ (u λ ) -C 7 ε 1 2 (p-k)(1-τ ) ≥ c λ -C 7 ε 1 2 (p-k)(1-τ ) .
Passing to the limit as ε → 0, we finally get

E λ (u λ ) = c λ . By definition of c λ , u λ satisfies E λ (u λ ) = min v∈W 1,p 0 (Ω) E λ (v)
and since E λ is Gâteaux differentiable, u λ satisfies the Euler-Lagrange equation associated to E λ :

∀v ∈ W 1,p 0 (Ω), Ω |∇u λ | p-2 ∇u λ .∇v dx = Ω f λ (x, u λ )v dx.
In particular, setting v = (u λ ) -∈ W 1,p 0 (Ω), by weak maximum principle it follows that u λ ≥ 0 a.e. in Ω. Moreover, since M u λ is a super-solution of (P λ ), for all non-negative v ∈ W 1,p 0 (Ω),

Ω |∇(M u λ )| p-2 ∇(M u λ ).∇v dx ≥ Ω K(x) [λ(M u λ ) q -(M u λ ) r ] v dx. Setting v = (u λ -M u λ ) + ∈ W 1,p 0 (Ω), we obtain 0 = Ω f λ (x, u λ ) -K(x) [λ(M u λ ) q -(M u λ ) r ] (u λ -M u λ ) + dx ≥ Ω |∇u λ | p-2 ∇u λ -|∇(M u λ )| p-2 ∇(M u λ ) .∇ (u λ -M u λ ) + dx.
Using lemma 5.1, ∇ (u λ -M u λ ) + = 0 a.e. in Ω and by Poincaré's inequality u λ ≤ M u λ a.e. in Ω. Finally

I λ (u λ ) = E λ (u λ ) = min v∈W 1,p 0 (Ω) E λ (v) ≤ E λ (ϕ 0 ) = I λ (ϕ 0 ) < 0,
therefore u λ is a non-trivial weak solution of (P λ ) .

Compacted support of the solution.

In this section we define (5.16) in such a way that g λ is positive and increasing on the interval (0, a * ). Let us start by stating a result which guarantees the existence of an appropriate super-solution of (P λ ) near the boundary. Lemma 5.2. Let u λ ∈ W 1,p 0 (Ω) ∩ L ∞ (Ω) be a weak solution of (P λ ). Then u λ ∈ C Ω and there exist δ * > 0, M > 0 and α ∈ (1, p ) such that

g λ (t) def = t r -λt q , t ∈ [0, +∞) and a * def = r λq 1 q-r
u λ (x) ≤ M ϕ 1 (x) α in Ω δ * .
In the proof of this lemma, we will use the following weak comparison principle: Proposition 5.2. Let us consider the Dirichlet problems:

-∆ p u -b(x, u) = f in Ω, u = f on ∂Ω (5.17) and -∆ p v -b(x, v) = g in Ω, v = g on ∂Ω.
(5.18)

Assume that f ≤ g in L p (Ω), f ≤ g in W 1 p ,p (∂Ω), u, v ∈ W 1,p (Ω)
are any weak solutions of the Dirichlet problems (5.17) and (5.18), respectively and b(x, •) : R → R is non-increasing for a.e. x ∈ Ω. Then, u ≤ v in Ω.

Proof: See proposition 2.3 in Cuesta-Takáč [START_REF] Cuesta | A strong comparison principle for positive solutions of degenerate elliptic equations[END_REF].

Proof of Lemma 5.2: With the previous notations, the set ω * def = {x ∈ Ω, u λ (x) ≤ a * } contains a neigbourhood of ∂Ω and there exists δ 0 > 0 such that Ω δ 0 ⊂ ω * . Since u λ is bounded, there exists C * > 0 large enough such that u λ ≤ C * ϕ 1 on ∂Ω δ 0 . Hence, u λ and C * ϕ 1 satisfy

-∆ p u λ ≤ -∆ p (C * ϕ 1 ) in Ω δ 0 , u λ ≤ C * ϕ 1 on ∂Ω δ 0 . (5.19)
Therefore, by the weak comparison principle u λ ≤ C * ϕ 1 in Ω δ 0 . From this estimate and the interior regularity result of Serrin [START_REF] Serrin | Local behavior of solutions of quasi-linear equations[END_REF], u λ ∈ C Ω .

Let M > 0 and α ∈ (1, p ), we want to construct a super-solution v to (P λ ) near the boundary such that v def = M ϕ 1 α . Similarly to the proof of Lemma 3.2, there exists a δ 1 > 0 only depending on Ω, p, M and α such that:

∆ p v ∼ (M α) p-1 (α -1)(p -1)δ(x) (α-1)(p-1)-1 in Ω δ 1 (5.20) and K(x)(v r -λv q ) ∼ M r L(δ(x))δ(x) αr-k in Ω δ 1 .
(5.21)

Precisely,

δ 1 def = min ( δ * , ε * K 2 " (α -1)(p -1) 2λ 1 « 1 p , 1 K 2 " 1 2λ « 1 α(q-r) " 1 M « 1 α ) ,
where ε * and δ * are defined in (2.3). By definition of δ 1 , choosing α > 1 small enough,

δ 1 = ε * K 2 (α-1)(p-1) 2λ 1
1 p and we can impose

M ≤ 2 4 inf δ≤δ 1 L(δ)δ -(α(p-r-1)-(p-k)) α p-1 (α -1) 3 5 1 p-(1+r)
.

(5.22)

Then, by estimates (5.20) and (5.21), v is a super-solution of (P λ ) in Ω δ 1 . Moreover, if we set

δ2 def = min 8 > < > : δ0, a * C * K2 , 0 @ a * α p-1 (α -1) C1 inf δ<δ 1 L(δ) 1 A p-(1+r) p-k 9 > = > ; , u λ ≤ a * and v ≤ a * in Ω δ 2 .
Finally, setting δ * def = min{δ 1 , δ 2 } and choosing α close enough to 1, u λ and v satisfy:

   -∆ p v -K(x)g λ (v) ≥ 0 in Ω δ * , -∆ p u λ -K(x)g λ (u λ ) = 0 in Ω δ * , v ≥ u λ on ∂Ω δ * .
Note that the third assertion is a consequense of (5.19) and (5.22). Since v → -K(x)g λ (v) is non-increasing in (0, a * ) for all x ∈ Ω δ * , applying the weak comparison principle of proposition 5.2, it follows that

u λ ≤ v in Ω δ * . Proposition 5.3. Let k ∈ [1 + r, p) and let u λ ∈ W 1,p 0 (Ω) ∩ L ∞ (Ω)
be a weak solution of (P λ ), therefore u λ has a compact support.

proof: For s ∈ [0, a * ], we define G λ (s) def = s 0 g λ (t) dt. Since r < p -1, we have a * 0 G λ (s) -1 p ds < +∞.
(5.23) Note that this above equation is close to condition (2) in Vásquez [START_REF] Vásquez | A strong maximum principle for some quasilinear elliptic equations[END_REF] and implies that u λ may be not positive everywhere in Ω. Precisely, let us fix ε < a * (small) and

δ ε def = ε M 1 α with M and α defined in Lemma 5.2 in such a way that u λ < ε in Ω δ ε . Let us define for t ∈ [0, a * ], h(t) def = ε t G λ (s) -1 p ds. h is a C 2 bijection from [0, a * ] to [h(a * ), h(0)] and h (t) = -G λ (t) -1 p < 0, for t ∈ (0, a * ).
Then h -1 is also twice differentiable on (h(a * ), h(0)) and we get,

(h -1 ) (y) = -G λ (h -1 (y)) 1 p and (h -1 ) (y) = 1 p g λ h -1 (y) G λ h -1 (y) 2-p p for y ∈ (h(a * ), h(0)).

Now let us define

j(x) def = inf ∂Ω δ ε ϕ 1 ϕ 1 (x) g λ s 2 s 2 -1 p ds, for x ∈ Ω δε and J(x) def = min {j(x), h(0)} , for x ∈ Ω δ ε . Remark that j(x) > h(a * ) for x ∈ Ω δ ε provided ε is sufficiently small. Indeed, h(a * ) ≤ C 1 ε p-(r+1) p -(a * ) p-(r+1) p < -C 2 ε p-(r+1) αp ≤ j(x), for x ∈ Ω δε ,
with C 1 and C 2 two positive constants independant of ε.

With all this notations, we finally define the function w in Ω δε by

w(x) def = h -1 (J(x)) , for x ∈ Ω δ ε .
In other words,

ε w(x) G λ (s) -1 p ds = J(x), for x ∈ Ω δε .
Using this last caracterisation, w is non-negative in Ω δ ε and w ≤ a * in Ω δ ε . Moreover, w vanishes when δ(x) is small. Indeed, for s ∈ (0, a * ),

g λ (s)s > G λ (s) > s s 2 g λ (t) dt > s 2 g λ s 2 .
Then for ε ∈ (0, a * ),

ε 0 G λ (s) -1 p ds < ε 0 s 2 g λ s 2 -1 p ds. So, j(x) -→ δ(x)→0 inf ∂Ω δ ε ϕ 1 0 g λ s 2 s 2 -1 p ds ≥ K 1 δ ε 0 g λ s 2 s 2 -1 p ds > ε 0 G λ (s) -1 p ds,
for ε > 0 small, from which, together with the definitions of J and w, w has a compact support in Ω.

Then, to complete the proof, it's enough to show that

u λ ≤ w in Ω δ ε . Because of the compacted support of w, J ∈ W 1,p (Ω δ ε ), thus w = h -1 • J ∈ W 1,p (Ω δ ε ) and satisfies ∇w = -G λ (w) 1 p ∇j in D (Ω δ ε ) .
Then,

∆ p w + G λ (w) 1 p ∆ p j = 1 p g λ (w)|∇j| p in D (Ω δ ε ) .
In this equation, provided ε is sufficiently small

1 p |∇j| p = 1 p |∇ϕ 1 | p ϕ 1 2 g λ ϕ 1 2 -1 K(x) in Ω δ ε and ∆ p j = 1 p |∇ϕ 1 | p ϕ 1 2 g λ ϕ 1 2 -1 p 1 2 g λ ϕ 1 2 + ϕ 1 4 (g λ ) ϕ 1 2 +λ 1 ϕ 1 p-1 ϕ 1 2 g λ ϕ 1 2 1 p ≥ 0 in Ω δ ε .
Hence, ∆ p w ≤ K(x)g λ (w) in Ω δ ε . Moreover, since g λ ≥ 0 on ∂Ω δ ε , we have u λ (x) ≤ ε ≤ w(x) on ∂Ω δε . Therefore, by the weak comparison principle of proposition 5.2,

u λ (x) ≤ w(x) in Ω δ ε .
Proof of Theorem 2.2: Since u λ is compactly supported in Ω, inequality (3.1) is also satisfied when k ≥ 1 + r, which implies the existence of a critical parameter λ * * > 0 such that (P λ ) has no non-trivial solution for λ < λ * * . Thanks to above propositions 5.1 and 5.3 and remark 5.1, using the C 1,β regularity result of Lieberman [START_REF] Lieberman | Boundary regularity for solutions of degenerate elliptic equations[END_REF] and the same arguments as in paragraph 4.2, Theorem 2.2 follows.

Appendix A. Proof of Lemma 3.3

A.1. When 0 ≤ k < 1 + q. By (1.5), v def = mϕ 1 is a sub-solution of (Q) in Ω for m > 0 small enough. Now let us define f (x) def = M δ(x) -(k-q) L (δ(x)) in Ω, with M > 0. Let (k -q) + < ε < 1, therefore 0 < f (x) ≤ C 1 δ(x) -ε in Ω. Thus if we consider the problem Q -∆ p v = f in Ω, v = 0 on ∂Ω, v > 0 in Ω, by a result of Giacomoni, Schindler and Takáč [6] Q has a unique solution v ∈ C 1,α Ω , with α ∈ (0, 1) and v ∼ δ(x) in Ω. Therefore, -∆ p v ≥ K(x)v q
in Ω for M > 0 sufficiently large. Hence we get both sub-and super-solution of the problem Q which behave like the distance function δ(x) in Ω. Using the same sub-and super-solution method as section 4, we get a solution

v ∈ W 1,p 0 (Ω) ∩ C Ω satisfying v(x) ∼ δ(x) in Ω. Now let w ∈ W 1,p 0 (Ω) ∩ C Ω be a solution to Q satisfying w(x) ∼ δ(x) in Ω. Then we can define C * def = sup {C > 0, Cw ≤ v in Ω} ∈ R.
It's easy to see that C * > 0 and C * w ≤ v in Ω, so for all x ∈ Ω we get

-∆ p (C * ) q p-1 w(x) = K(x) (C * w(x)) q ≤ K(x)v(x) q = -∆ p (v(x)).
If we suppose C * < 1, the weak maximum principle implies (C * )

q p-1 w ≤ v in Ω. But (C * ) q p-1 > C * because C * < 1 and q < p -1, therefore C * w < (C * ) q p-1 w v in Ω, wich contradicts the definition of C * . So C * ≥ 1 and w ≤ C * w ≤ v in Ω.
Interchanging the role of w and v, we finally get that w = v and this proves the uniqueness of the solution of Q in the convex set

Λ 1 def = v ∈ W 1,p 0 (Ω) ∩ C Ω , v(x) ∼ δ(x) in Ω . A.2. When 1 + q ≤ k ≤ p. For t ∈ (0, d] we define Θ(t) def = exp 2d t y(s) s ds ,
with y ∈ C ([0, 2d]) ∩ C 1 ((0, 2d]) such that y(0) = 0 and lim t→0 + ty (t)

y(t) = 0 in order to satisfy lim t→0 + tΘ (t) Θ(t)
= 0 and lim

t→0 + tΘ (t) Θ (t) = -1. (A.1) Let β ∈ [0, 1], for x ∈ Ω we also define w(x) def = ϕ 1 (x) β Θ (ϕ 1 (x)) in Ω.
A direct computation of -∆ p w in Ω gives us

-∆ p w = " Θ(ϕ 1 ) " p-1 ϕ 1 (β-1)(p-1)-1 " β + ϕ1Θ (ϕ1) Θ(ϕ 1 ) « p-2 × h" β + ϕ 1 Θ (ϕ 1 ) Θ(ϕ 1 ) " λ1ϕ1 p + (p -1)|∇ϕ1| p " β(1 -β) -2β ϕ 1 Θ (ϕ 1 ) Θ(ϕ 1 ) - ϕ 2 1 Θ (ϕ 1 ) Θ(ϕ 1 )
"i .

From now, we will distinguish the following cases.

A.2.1. Case 1: 0 < β < 1.

There exists ε > 0 sufficiently small such that for x ∈ Ω ε ,

β 2 ≤ β + ϕ 1 (x)Θ (ϕ 1 (x)) Θ(ϕ 1 (x)) ≤ 3β 2 and β(1 -β) 2 ≤ β(1 -β) -2β ϕ 1 (x)Θ (ϕ 1 (x)) Θ(ϕ 1 (x)) - ϕ 1 (x) 2 Θ (ϕ 1 (x)) Θ(ϕ 1 (x)) ≤ 3 2 β(1 -β).
Therefore we get

-∆ p w(x) ∼ Θ(ϕ 1 (x)) p-1 ϕ 1 (x) (β-1)(p-1)-1 in Ω, which implies -∆ p w(x) w(x) -q ∼ Θ(ϕ 1 (x)) p-(1+q) ϕ 1 (x) (β-1)(p-1)-1-qβ in Ω. When 1 + q < k < p, if we choose β = p-k p-(1+q) ∈ (0, 1) and y(t) = z(t) p-(1+q) for t ∈ [0, 2d], w satisfies -∆ p w(x) w(x) -q ∼ K(x) in Ω, w = 0 on ∂Ω, w > 0 in Ω.
Therefore there exist C 1 , C 2 > 0 such that C 1 w and C 2 w are respectively sub-and super-solutions of the problem (Q). Thus, (Q) has a solution

v ∈ W 1,p loc (Ω) ∩ C 0 Ω satisfying v(x) ∼ δ(x) p-k p-(1+q) L (δ(x)) 1 p-(1+q) in Ω. (A.2)
Using the same arguments as section A.1, we get the uniqueness of the solution in the set 1+q) in Ω .

Λ 2 def = v ∈ W 1,p loc (Ω) ∩ C 0 Ω , v(x) ∼ δ(x) p-k p-(1+q) L (δ(x)) 1 p-(
Moreover, u λ ∈ W 1,p 0 (Ω) if and only if the right hand term in the equation of problem (Q) is W -1,p (Ω), i.e. if and only if there exists a constant C > 0 such that

∀v ∈ W 1,p 0 (Ω), Ω K(x)u λ (x) q v(x)dx ≤ C v .
Using estimate (A.2), Hardy's and Hölder's inequalities and property (1.4) of the perturbation L, this condition is satisfied if k < 1+q + p-(1+q) p . Moreover, by (A.2), we have

Ω K(x) (u λ (x)) q+1 dx < +∞ only if k ≤ 1+q + p-(1+q) p . Then, as soon as k > 1+q + p-(1+q) p , u λ / ∈ W 1,p 0 (Ω). A.2.2. Case 2: β = 1.
The computation of -∆ p w becomes

-∆ p w = Θ (ϕ 1 ) (Θ(ϕ 1 )) p-2 1 + ϕ 1 Θ (ϕ 1 ) Θ(ϕ 1 ) p-2 × 1 + Θ(ϕ 1 ) ϕ 1 Θ (ϕ 1 ) λ 1 ϕ 1 p + (p -1)|∇ϕ 1 | p -2 - ϕ 1 Θ (ϕ 1 ) Θ (ϕ 1
) .

We choose Θ such that

C 1 ϕ 1 p ≤ - Θ(ϕ 1 )ϕ 1 p-1 Θ (ϕ 1 ) ≤ C 2 ϕ 1 p-1
near the boundary, that is equivalent to require

C 1 t ≤ - Θ(t) Θ (t) ≤ C 2 , for t > 0 small enough. (A.3) Hence, -∆ p w(x) w(x) -q ∼ -Θ (ϕ 1 (x)) Θ (ϕ 1 (x)) p-2-q ϕ 1 (x) -q in Ω.
To get

-∆ p w(x) w(x) -q ∼ ϕ 1 (x) -k L (ϕ 1 (x)) in Ω, we require t -(1+q) y(t)Θ(t) p-(q+1) ∼ t -k 2d z(s) s ds in (0, d].
This condition can be satisfied only if k = 1 + q. Then taking

Θ(t) = 2d t s -1 L(s)ds 1 p-(1+q)
, 0 < t ≤ d, Θ satisfies conditions (A.1) and (A.3). Thus, if k = 1 + q and

w(x) = ϕ 1 (x) 2d ϕ 1 (x) t -1 L(t)dt 1 p-(1+q)
in Ω, there exist C 1 , C 2 > 0 such that C 1 w and C 2 w are respectively sub-and super-solutions of (Q). Thus, there exists a solution

v ∈ W 1,p 0 (Ω) ∩ C Ω of (Q) satisfying v(x) ∼ δ(x) 2d δ(x) t -1 L(t)dt 1 p-(1+q) in Ω.
Using the same argument as section A.1 we get the uniqueness of the solution in the set

Λ 3 def =    v ∈ W 1,p 0 (Ω) ∩ C Ω , v(x) ∼ δ(x) 2d δ(x) t -1 L(t)dt 1 p-(1+q) in Ω    .
A.2.3. Case 3: β = 0.

In this case, we get

-∆ p w = ϕ(x) -1 Θ (ϕ 1 ) p-1 λ 1 ϕ 1 p -(p -1)|∇ϕ 1 | p ϕ 1 Θ (ϕ 1 ) Θ(ϕ 1 ) .
Hence,

-∆ p w(x) w(x) -q ∼ ϕ 1 (x) -1 Θ (ϕ 1 (x)) p-1 Θ(ϕ 1 (x)) -q in Ω.
Similarly as the previous case, to get

ϕ 1 (x) -1 Θ (ϕ 1 (x)) p-1 Θ(ϕ 1 (x)) -q ∼ ϕ 1 (x) -k L(ϕ(x)) in Ω,
we require

t -p Θ(t) p-(1+q) -y(t) p-1 ∼ t -k exp 2d t z(s) s dt in (0, t].
This condition can be satisfied only if k = p. Then if condition (3.2) holds and if we choose

Θ(t) = exp 2d t y(s) s ds = C t 0 s -1 L(s) 1 p-1 ds p-1 p-(1+q) , 0 < t ≤ d,
we get that Θ satisfies conditions (A.1) and (A.3). Thus if k = p and

w(x) = C ϕ 1 (x) 0 t -1 L(t) 1 p-1 dt p-1 p-(1+q)
, there exists C 1 , C 2 > 0 such that C 1 w and C 2 w are respectively sub-and super-solutions of (Q) and there exists a solution

v ∈ W 1,p loc (Ω) ∩ C 0 Ω of (Q) satisfying v(x) ∼ δ(x) 0 s -1 L(s) 1 p-1 ds p-1 p-(1+q) in Ω.
Using the same argument as section A.1, we get the uniqueness of the solution in the set

Λ 4 def = 8 < : v ∈ W 1,p loc (Ω) ∩ C 0 `Ω´, v(x) ∼ Z δ(x) 0 s -1 L(s) 1 p-1 ds ! p-1 p-(1+q) in Ω 9
R(u, v) def = |∇u| p -|∇v| p-2 ∇v.∇ u p v p-1 .
Moreover, L(u, v) = 0 in Ω if and only if there exists C > 0 such that u = Cv in Ω.

Proof: See Theorem 1.1 in Allegretto-Huang [START_REF] Allegretto | A Picone's identity for the p-Laplacian and applications[END_REF].

Lemma B.2. (Díaz-Saa inequality) For i = 1, 2 let w i ∈ L ∞ (Ω) such that w i ≥ 0 a.e. in Ω, w i 1 p ∈ W 1,p (Ω), ∆ p w i 1 p ∈ L ∞ (Ω) and w 1 = w 2 on ∂Ω. Moreover if w 1 w 2 , w 2 w 1 ∈ L ∞ (Ω), we have the inequality Ω   -∆ p w 1 1 p w 1 p-1 p + ∆ p w 2 1 p w 2 p-1 p   (w 1 -w 2 ) dx ≥ 0. (B.1)
Futhermore, (B.1) becomes an equality if and only if there exists C > 0 such that w 1 = Cw 2 a.e. in Ω.

Proof: See Lemme 2 in Díaz-Saa [START_REF] Díaz | Existence et unicité de solutions positives pour certaines équations elliptiques quasilinéaires[END_REF].

Proof of proposition 3.1: We argue by contradiction. If proposition 3.1 does not hold, there exist v ∈ W 1,p loc (Ω) ∩ C 0 Ω weak solution of (Q), η > 0 and ε > 0 satisfying v ≥ ηu a.e. in Ω and (3.3).

B.1.

Step 1: when q ≥ 0. We consider the following perturbated problem:

(Q n ) -∆ p v = K n (x)v q , v > 0 in Ω, v = 0 on ∂Ω, where (K n ) n∈N ⊂ L ∞ (Ω) is increasing sequense satisfying K n -→ n→+∞ K a.e.
in Ω. We will prove there exists a unique solution of (Q n ) in W 1,p 0 (Ω) and show that this solution is C 1,α Ω for some α ∈ (0, 1). Let us consider the functional,

I n (u) def = Ω |∇u| p dx, u ∈ V def = w ∈ W 1,p 0 (Ω), Ω K n (x)w q+1 dx = 1 .
Since V is a compact subset of W 1,p 0 (Ω), there exists a non-negative and non-trivial ṽn ∈ W 1,p 0 (Ω) satisfying

I n (ṽ n ) = min u∈V I n (u).
Therefore, from the Lagrange multiplier rule, there exists a Lagrange multiplier λ n > 0 such that

-∆ p ṽn = λ n K n (x) (ṽ n ) q
in Ω, ṽn = 0 on ∂Ω.

By homogeneity of the p-Laplacian operator, if we define

v n def = (λ n ) 1 p-(1+q) ṽn ∈ W 1,p 0 (Ω), v n satisfies -∆ p v n = K n (x)v n q , in Ω, v n = 0 on ∂Ω,
Since q < p -1 and K n ∈ L ∞ (Ω), using Moser iterations we prove that v n ∈ L ∞ (Ω) and due to the well known regularity result in Liebermann [START_REF] Lieberman | Boundary regularity for solutions of degenerate elliptic equations[END_REF], v n ∈ C 1,α Ω for a certain α ∈ (0, 1). The positivity of v n comes from the strong maximum principle in Vásquez [START_REF] Vásquez | A strong maximum principle for some quasilinear elliptic equations[END_REF] and v n is a solution of (Q n ). Now, let us prove the uniqueness of a such solution. Therefore, for that, we use the Díaz-Saa inequality (B.1). So let u n ∈ C 1,α Ω be an other solution of (Q n ), then

Ω -∆ p u n u n p-1 + ∆ p v n v n p-1 (u n p -v n p ) dx ≥ 0, (B.2) which implies Ω K n (x) 1 u n p-(1+q) - 1 v n p-(1+q) (u n p -v n p ) dx = 0.
Then inequality (B.2) becomes an equality, therefore by lemma B.1 there exits C > 0 such that u n = Cv n in Ω. Futhermore, by homogeneity arguments,

-∆ p (Cv n ) = K n (x)(Cv n ) q in Ω if C = 1, so u n = v n
in Ω and we get the uniqueness. Now, we will prove that for all n ∈ N v n ≤ v. For that, we apply a suband super-solution method in a compact subset of Ω. So let us fix n ∈ N and define (Ω m ) m∈N * an increasing sequense of smooth subdomains of Ω such that Ω m -→ m→∞ Ω in the Hausdorff topology with

∀m ∈ N * , 1 m + 1 < dist(∂Ω, ∂Ω m ) < 1 m .
Then we consider the following sequense of problems:

(Q n,m ) -∆ p u = K n (x)v q in Ω m , v = ηv on ∂Ω m , v > 0 in Ω m ,
with v ∈ W 1,p 0 (Ω)∩C Ω the sub-solution of (P). Since v ∈ W 1,p loc (Ω)∩L ∞ (Ω) and v ≥ ηv in Ω by hypothesis, using the same arguments as in the proof of Proposition 4.1, for all m ∈ N there exists v n,m ∈ W 1,p

(Ω m ) ∩ C Ω m unique solution of (Q n,m ). Movereover, v n,m satisfies ηv ≤ v n,m ≤ v in Ω. Now if we define ṽn,m def = 1 Ωm .v n,m
in Ω in order to extend v n in Ω by zero, the sequense (ṽ n,m ) m∈N * is an increasing sequense which converges pointwise to an element u n ∈ W 1,p loc (Ω) ∩ C 0 Ω solution of (Q n ), by similar arguments as in the proof of Proposition 4.1. Then, the uniqueness argument implies u n = v n in Ω and then ∀n ∈ N, v n ≤ v in Ω.

B.2.

Step 2: when q < 0. Let us define the following problem:

(Q n ) -∆ p v = K n (x) v + 1 n q , v > 0 in Ω, v = 0 on ∂Ω.
Using a similar method as step 1, we get the existence and the uniqueness of a sequense of weak solutions of (Q n ) in W 1,p 0 (Ω). (2) When q < 0 arguing similarly as in the first case, we get with values in W -1,p (Ω). Moreover, the components a i of the vector field a : Ω × R N → R N , a = (a 1 , • • • , a N ), are functions of x and η ∈ R N , such that for i, j ∈ {1, • • • , N }, a i ∈ C(Ω × R N ) and ∂a i ∂η j ∈ C(Ω × (R N \ {0})). We assume that a satisfies the following ellipticity and growth conditions:

β p Ω |∇ϕ 1 | p ϕ 1 (β-1)p dx ≥ Ω K n (x)
(H1) There exist some constants κ ∈ [0, 1], γ, Γ ∈ (0, +∞) and α ∈ (0, 1), such that for all x, y ∈ Ω, all η ∈ R N \ {0} and ξ ∈ R N , proof: The proof is very simlar to the proof of Therorem B.1 in Giacomoni, Schindler and Takáč [START_REF] Giacomoni | Sobolev versus Hölder local minimizers and existence of multiple solutions for a singular quasilinear equation[END_REF]. In particular, condition (C.6) replacing the growth condition (B.9) imply the estimate (B.17) in [START_REF] Giacomoni | Sobolev versus Hölder local minimizers and existence of multiple solutions for a singular quasilinear equation[END_REF].

The C 1,β regularity of u λ is directy proved applying this proposition with a(x, η) def = |η| p-2 η and f (x) def = K(x) [λu λ (x) q -u λ (x) r ] for x ∈ Ω and η ∈ R N .

Example 1 . 1 .

 11 Let m ∈ N * and A > 0 large enough. Let us defineL(t) = m n=1 log n A t µn , t ∈ (0, 2d] with log n def = log • • • • • log (n times) and µ n > 0. Then L ∈ K.

. 1 )

 1 It follows that ũk -→ k→∞ u λ in D (Ω) and u λ satisfies (1.1). Using inequality (4.1) and Hardy's inequality, K

  Finally applying proposition C.1 of the appendix, we get the C 1,β Ω regularity of u λ . 4.2. Existence of Λ 1 . Let us define Λ 1 def = inf {λ > 0, (P λ ) has a positive solution} . By Lemma 3.1 and the first step of this proof, λ

Proof of Theorem 2 . 1 : 5 . 2 5. 1 .

 21521 The proof follows from proposition 4.1 and subsection 4.2. Proof of Theorem 2.Existence of a solution under condition (2.1) or (2.2). Proposition 5.1.

1 |∇v| p- 2

 12 Proof of proposition 3.To prove this proposition, we need the two following lemmas: Lemma B.1. (Picone's Identity) Let u, v ∈ C 1 Ω two positive functions satisfying the Hopf 's lemma. Then, L(u, v) def = |∇u| p + (p -1) ∇v.∇u satisfies L(u, v) ≥ 0 in Ω and L(u, v) = R(u, v) where

B. 3 .( 1 )Ω|∇ϕ 1 | p ϕ 1 (ϕ 1 βpv

 31111 Step 3: Applying Picone's Identity withu = ϕ 1 β ∈ C 1 Ω , where β = p-1+ε p and v = u n ∈ C 1 Ω , we get 0 ≤ Ω |∇u| p -|∇v| p-2 ∇v.∇ u p v p-1 dx. (B.3) When q ≥ 0 this inequality becomesβ p n p-(q+1) dx. Therefore, passing to the limit when n → +∞, there exists C > 0 such thatinf y∈ω 1 v(y) p-(q+1) ω K(x)ϕ 1 p-1+ε dx ≤ C, ∀ω ⊂⊂ Ω.This inequality does not hold for ω close enough to Ω, i.e when dist(Ω, ω) is sufficiently small, because Ω K(x)ϕ 1 p-1+ε dx = +∞, by hypothesis.

ϕ 1 βp v n + 1 n 1 (

 111 p-(q+1) dx.Therefore passing to the limit when n → +∞,inf y∈ω v(y) + 1) p-(q+1) ω K(x)ϕ 1 p-1+ε dx ≤ C, ∀ω ⊂⊂ Ω.And we conclude as above.Appendix C. C 1,β regularityWe consider the following quasilinear elliptic boundary value problem,(P) -div(a(x, ∇u)) = f (x) in Ω, u = 0 on Ω.In this equation, f ∈ L ∞ loc (Ω) and div(a(x, ∇u)) i (x, ∇u(x)), for x ∈ Ω and u ∈ W 1,p 0 (Ω) (C.1)

  a i (x, 0) = 0, for i = 1, • • • , N, (C.2) N i,j=1 ∂a i ∂η j (x, η)ξ i ξ j ≥ γ(κ + |η|) p-2 |ξ| 2 , (C.3) N i,j=1 ∂a i ∂η j (x, η) ≤ Γ(κ + |η|) p-2 , (C.4) N i=1 |a i (x, η) -a i (y, η)| ≤ Γ(1 + |η|) p |x -y| α . (C.5)We remark that condition (C.2) through (C.5) are motivated by the elliptic boundary value problem,(P) -∆ p u = f (x) in Ω, u = 0 on Ω.Finally, we impose the following growth condition on the function f :(H2) f ∈ L ∞ loc (Ω)and there exist some constants c > 0 and ε ∈ (0, 1) such that, for almost all x ∈ Ω,|f (x)| ≤ cδ(x) -ε . (C.6)Proposition C.1. Assume that a(x, η) satisfies the structural hypotheses (C.2) through (C.5) and f (x) satisfies the growth hypothesis (C.6). Let u ∈ W 1,p 0 (Ω) be a weak solution of the problem (P). In addition assume0 ≤ u(x) ≤ Cδ(x) for almost all x ∈ Ω, (C.7)where C > 0 and u ∈ W 1,p 0 (Ω) is a weak solution ofP -div(a(x, ∇u)) = f (x) in Ω, u = 0 on Ω.Then there exist constants β ∈ (0, α) and M ≥ 0 only depending on Ω, N, p, on γ, Γ, α in (C.2) through (C.5), on the constants c, ε in (C.6) and on the constant C in (C.7), but not on κ ∈ [0, 1], such that u satisfies u ∈ C 1,β Ω and u C 1,β (Ω) ≤ M.


P O Bouchard

DAMAGE AND DISCRETE CRACK PROPAGATION MODELLING -SOME RESULTS AND CHALLENGES FOR 2D AND 3D CONFIGURATIONS

In this paper we present a technique to model damage and fracture mechanics using remeshing in 2D and 3D configurations. We use the finite element software FORGE2® and FORGE3® which can deal with elastic, elastoplastic and elastic-viscoplastic materials in large deformation. Coupled damage models (Lemaître and Gurson) have been implemented to model the progressive mechanical degradation of the material during its deformation. Once damage reaches a critical value, fracture has to be modelled. 2D remeshing technique and crack propagation criteria are presented to model automatic discrete crack propagation for different configurations. Extension to 3D modelling of fracture is also discussed.

INTRODUCTION

Damage and Fracture Mechanics have been studied for many years now. During the last 50 years, numerous different damage models and crack propagation techniques have been introduced and implemented in different finite element software. On a macro-scale, damage is characterized by a progressive loss of rigidity in the mechanical behaviour of a material. From a micro-scale point of view, damage can be described by the three well-known stages of nucleation of micro-voids, growth of these voids and finally coalescence leading to fracture. Among the numerous models presented in the literature, the evolution of damage can be coupled to the mechanical properties of the material or not. Once damage reaches a critical value, a crack becomes initiated and has to be modelled. Numerical modelling of crack propagation has been a real challenge for many years now. The initial and more natural approach was based on a discrete crack propagation modelling approach using remeshing techniques. To avoid numerical difficulties due to geometric topological changes and remeshing during crack propagation, numerous different techniques have been introduced: embedded crack models, mesh-free techniques, XFEM and so on. Each approach has its own advantages and drawbacks, and I shall not in this paper discuss or compare each technique to the others. Depending on the application, one approach can give better results than the others. Sometimes, for instance in high shearing processes such as blanking, it is even sufficient to use a simple "kill-element" technique. After a brief description of the finite element software used for this study ( §2), we introduce the notion of damage and its coupling with the material's mechanical behaviour ( §3). Finally, discrete crack propagation modelling is presented and applied to 2D and 3D configurations examples ( §4).

FORGE2® AND FORGE3®

FORGE2® and FORGE3® have been developed to model large deformation of elastic, elastoplastic and elastic-viscoplastic materials. They are based on a mixed velocity-pressure formulation. The so-called Mini-element is used in FORGE3® (Arnold [START_REF] Arnold | A stable finite element for Stokes equations[END_REF]). It is based on linear isoparametric tetrahedra and a bubble function is added at element level in order to satisfy the Brezzi/Babuska condition. The space discretization based on this element associated to the incremental formulation of the virtual work principle lead to a set of discretized non-linear equations. The well-known iterative Newton-Raphson linearization method is used. A one step Euler scheme enables to compute the solution at time t+δt when the solution at time t is known. FORGE2® and FORGE3® can both handle multimaterial structures using a nodal incremental form of the penalty technique (Pichelin [START_REF] Pichelin | An application of a master-slave algorithm for solving 3D contact problems between deformable bodies in forming processes[END_REF]). Besides this, an automatic adaptive remesher enables to deal with large deformation without losing accuracy. The well-known Delaunay Triangulation is used in 2D and a topological remeshing technique enables to deal with automatic 3D remeshing (Gruau [3]).

DAMAGE MODELLING

Continuum damage mechanics is a constitutive theory that describes the progressive loss of material integrity due to the propagation and coalescence of microcracks, microvoids, and similar defects. Beyond a certain value of strain, void nucleation and void growth appear in the material: this phenomenon, called damage, allows to model the ductile fracture of materials. When these voids reach a critical size, they coalesce and give raise to instabilities or cracks propagation. Damage models are generally based on the study of void growth, using different parameters such as triaxiality, maximal principal stress, plastic strain and so on. The first models to be used were uncoupled damage models, which means that the damage law does not influence the mechanical properties of the material. The damage parameter is computed using an integral of a strain and stress function, and its distribution can be computed in a postprocessing step. This approach is easy to introduce in a numerical software, but is quite unrealistic because the damage evolution does not influence the material properties. In order to better represent the evolution of damage in materials, coupled damage models have been proposed. In this approach, damage and mechanical properties are directly linked and the material fracture is modeled by a progressive decrease of the global response of the structure. Contrary to the uncoupled approach, coupled damage models are quite difficult to introduce in numerical software, but are closer to the physical phenomenon of micromechanical fracture of ductile materials. Some of these models use the notion of effective stress which represents the actual stress transmitted by the bulk material between the microdefects. Another frequently used approach consists in introducing a damage variable f v which represents the volumetric fraction of voids in the material. The parameter f v is then used in the constitutive laws of the material and interacts with the others state variables. The damage model of Tvergaard-Needleman, based on the model introduced by Gurson 3 belongs to these approaches. The well-known Lemaître model [START_REF] Lemaître | A continuum damage mechanics model for ductile fracture[END_REF] and Gurson-Tvergaard-Needleman model [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth : Part I-Yield criteria and flow rules of porous ductile media[END_REF][START_REF] Tvergaard | Material Failure by void growth to coalescence[END_REF] have been implemented in FORGE3®. Details on these models can be found in [START_REF] Bouchard | Damage and Adaptive Remeshing applied to 3D modeling of blanking and Milling[END_REF]. Once damage reaches a critical value, fracture has to be taken into account. An easy way to represent fracture is the so-called "kill element" technique. When the damage parameter reaches a critical value inside an element, the element mechanical contribution to the stiffness matrix is set to zero. Coupled with adaptive remeshing, this technique enables to model fracture easily in 3D configurations. In Figure 1, a Lemaître damage model with adaptive remeshing and kill element enables to model a blanking process with accuracy. Figure 1. Blanking simulation: isovalues of damage during fracture However this technique involves a loss of volume during the simulation and the stress singularity at the crack tip is not properly represented. In some cases, when the crack path has to be correctly modelled, a more accurate technique is preferable.

CRACK PROPAGATION MODELLING

Our approach is based on discrete crack propagation, and we use automatic remeshing. Crack initiation is based both on critical stress or critical damage. Once a crack is initiated the propagation direction is computed using one of the following criteria: maximum circumferential stress criterion, minimum strain energy density criterion or maximum strain energy release rate criterion (Bouchard [START_REF] Bouchard | Crack propagation using an advanced remeshing technique[END_REF][START_REF] Bouchard | Numerical modelling of crack propagation: automatic remeshing and comparison of different criteria[END_REF]). In this last criterion, the strain energy release rate is computed using the Gθ method which appears to be both efficient and accurate (Bouchard [START_REF] Bouchard | Gtheta Method Applied to Crack Propagation Modelling[END_REF]). However, in high shearing configurations, when mode II becomes predominant, such criteria are no more available. In such cases, a maximum shear stress criterion has to be used (Bouchard [11]). Once the crack propagation direction has been predicted, a new outline representing the crack advance is added to the previous outline of the part. Then, a remeshing stage -based on the Delaunay triangulation -is performed and the crack opens naturally due to tensile stresses applied on its edges (Figure 2.): Figure 2. propagation of a crack in a 2D mesh The example in Figure 3, performed experimentally by Sumi, shows the growth of a crack starting from a fillet in a structural member. The crack propagation depends on the welding residual stresses and the bending stiffness of the structure. For the sake of simplicity, residual stresses are not taken into account. The bending stiffness of the structure is modified by varying the size of the bottom I-beam presented in figure 2 from 15mm to 315mm with an intermediate value of 115mm. A linear elastic plane strain simulation is performed with a Young modulus and a Poisson ratio of E=200GPa, and ν=0.3 respectively. The initial crack length is a 0 =5mm, and the maximum circumferential stress criterion is used to compute the crack path. Numerical simulations in figure 3 show the important influence of the bottom I-beam rigidity on the crack path: Figure 3. Crack growth from a fillet and influence of the support beam rigidity on the crack path.

The following example concerns a pre-cracked part with an inclusion. In [START_REF] Gruau | 3D tetrahedral, unstructured and anisotropic mesh generation with adaptation to natural and multidomain metric to appear[END_REF] we have studied the propagation of a crack in a planar part with an off-center hole. It has been shown that the crack is attracted by the hole since it creates a stress drop in this region. In the present example, we replace the hole by an inclusion, and we study the influence of this inclusion on the crack path. A rectangular part is pre-cracked and submitted to a tensile test. This part contains an inclusion which may be more rigid or less rigid than the matrix. If E matrix is the Young modulus associated with the matrix, and E incl the one associated with the inclusion, we define R as the ratio : R=E matrix /E incl . The numerical simulation is performed in plane strain, and the maximum strain energy release rate is used. This example shows the ability of the software to deal with multimaterial applications. Figure 4.a shows that for a soft inclusion -the inclusion is less rigid than the matrix -the crack is still attracted by the inclusion. The crack reorientation is however less pronounced than the one obtained with a hole. Conversely, if the inclusion is more rigid (see figure 4.b), the crack is repulsed. In 3D configurations, discrete crack propagation modelling is far more complex from a topological point of view since we have to deal with complex non planar surfaces representing crack faces. Finite element software that manages to model 3D discrete crack propagation are extremely rare (Carter [12], Schöllmann [START_REF] Schöllmann | Development of a new software for adaptive crack growth simulations in 3D structures[END_REF]). Our approach is still based on a modification of the part topology coupled with surface and volume remeshing. Once the new crack front has been localized, the surface mesh of the part is modified in order to insert new surface elements representing the new faces of the crack advance. If necessary, a topological improvement of the surface mesh can be performed at this stage. A volume remeshing of the new geometry is then performed. The 3D topological remeshing technique used here is both efficient and well adapted for such complex geometries [START_REF] Gruau | 3D tetrahedral, unstructured and anisotropic mesh generation with adaptation to natural and multidomain metric to appear[END_REF]. Figure 5 shows a crack propagating in the same 3D T-shape structure presented in figure 3. This example is purely intended for numerical feasibility purposes, since the crack propagation direction is determined by the user in this example. However, it demonstrates the ability of the 3D topological remesher to deal with complex 3D modifications such as crack propagation.

Figure 5. 3D crack growth from a fillet 5 CONCLUSION Recent advances in finite element software enable to model 3D complex forming processes or structural analysis applications. In such simulations the use of accurate damage models and of fracture modelling is sometimes necessary. Coupled damage models enable to couple the evolution of damage with the mechanical behaviour of materials. Once a critical damage value is reached, fracture has to be modelled explicitly. Our approach is based on discrete crack propagation and automatic remeshing. In 2D, different crack propagation criteria have been implemented and enable to model automatically crack propagation with accuracy. In 3D

Soft inclusion a)

Hard inclusion b)

Figure 4 .

 4 Figure 4. Crack propagation in a part with a) a soft inclusion R=10, b) a hard inclusion R=0.1

configurations, discrete crack propagation is more complex. In some cases, a simple kill element technique is accurate enough to model fracture. When the crack path becomes important, 3D discrete crack propagation is needed. We have extended our 2D developments to 3D configurations to model crack propagation. The first results are promising and have to be generalized to more complex geometries.


Irinel-Constantin Morȃrescu
email: constantin.morarescu@ensem.inpl-nancy.fr.

Carlos Canudas-De-Wit

Highway traffic model-based density estimation

Keywords: Highway traffic analysis, observer, switching mode model I.-C

The travel time spent in traffic networks is one of the main concerns of the societies in developed countries. A major requirement for providing traffic control and services is the continuous prediction, for several minutes into the future. This paper focuses on an important ingredient necessary for the traffic forecasting which is the real-time traffic state estimation using only a limited amount of data. Simulation results illustrate the performances of the proposed state-estimation technique.

I. INTRODUCTION

The increased travel time in congested sections has a dramatic social and economic impact. This led to an increasing research on freeway traffic control and development of intelligent transportation systems which are able to provide continuous forecasting of the traffic. A nice survey on the existing techniques for short-term traffic flow prediction can be found in [3]. One of the prerequisites for continuous prediction is an efficient real-time traffic conditions estimation using only a limited amount of data [START_REF] Muñoz | Traffic density estimation with the cell transmission model[END_REF], [START_REF] Muralidharan | Imputation of ramp flow data for freeway traffic simulation[END_REF]. The real-time density estimation is a challenging problem since the traffic is described by a system which is observable only when a segment situated between two vehicle detector stations (sensors) is entirely congested or entirely free.

Simulation modeling is a popular tool for analyzing transportation problems. Several studies focus on the validation of different microscopic or macroscopic models. Once the models validated they are used for open-loop stateestimation, forecasting and control ( [START_REF] Jia | The PeMS Algorithms for Accurate, Real-Time Estimates of g-factors and Speeds from Single-Loop Detectors[END_REF], [START_REF] Lin | Validating the Basic Cell Transmission Model on a Single Freeway Link[END_REF]). Other papers present imputation techniques to determine the missing onramp and off-ramp flows [START_REF] Muralidharan | Imputation of ramp flow data for freeway traffic simulation[END_REF].

In [START_REF] Gomes | Behavior of the Cell Transmission Model and Effectiveness of Ramp Metering[END_REF] it was shown that the traffic dynamics has two globally asymptotically stable equilibrium points which correspond to strictly feasible and infeasible demand (the link is either entirely free or entirely congested). This is the reason why approximating the dynamics using random switches between the mentioned two modes, the density is fairly well estimated. The drawback is, as observed in [START_REF] Kurzhanskiy | Set-valued estimation of freeway traffic density[END_REF], the system is no longer a conservation law (vehicles may appear/disappear in the network). In order to overcome this inconvenient we use in this paper a deterministic constrained model that reduce the number of possible affine dynamics of the system and preserve the number of vehicles in the network. Moreover, we use the vehicles conservation law to guaranty that the estimation error does not increase during the unobservable modes. This model is used to recover the state of the traffic network and precisely localize the eventual congestion front. The state of the network is recovered using what we call forward/backward observers.

The highway network is designed as a sequence of nodes relied by links. Since the sensors are located to the node level, one can easily determine if the node is over saturated or under saturated. Thus, the main concern for the traffic conditions is the estimation of the density inside the links. In order to better locate the congestions appearing into the network, each link is partitioned in several cells. It is worth noting here that the estimation problem is decentralized to the link level. In other words the density of the cells belonging to a link is estimated using only the data provided by the sensors located on the link boundary.

The structure of the paper is as follows: in Section II we introduce the deterministic constrained model that describe the density dynamics. Section III is devoted to observability of the system under consideration and in Section IV we propose an observer design. Section V focuses on the density estimation for each cell of a highway segment. In this section we also study the global observability of the traffic state. Simulation results are provided in Section VI before some concluding remarks.

II. CONSTRAINED SWITCHING MODEL FOR TRAFFIC

ESTIMATION

The traffic dynamics models are based on the car conservation principle. The simplest continuous macroscopic traffic model, involving only the density ρ, is the LWR model introduced in [START_REF] Lighthill | On kinematic waves II: A theory of traffic flow on long crowded roads[END_REF], [START_REF] Richards | Shock Waves on the Highway[END_REF]. The constitutive assumption of this model, motivated by experimental data, is that the vehicles tend to travel at an equilibrium speed v = v(ρ) where ρ represents the density of a specific section at a specific time. Thus the equilibrium speed depends implicitly on the location and on the time. Since the flow is defined as ϕ(ρ) = ρv(ρ), one can depict an equilibrium flow function ϕ = ϕ(ρ) called the fundamental diagram in traffic engineering.

In the sequel, we use the macroscopic traffic flow model called the switching mode model (SMM) derived from the cell transmission model (CTM) proposed by Daganzo [START_REF] Daganzo | The Cell Transmission Model: A Dynamic Representation of Highway Traffic Consistent with the Hydrodynamic Theory[END_REF]. The SMM is a piecewise linear state-dependent model in which the flow on each interface is a trade-off between the supply and the demand

ϕ i = min{D i-1 , S i } (1)
with

D i-1 = min{v i-1 ρ i-1 , ϕ m,i-1 }, S i = min{ϕ m,i , w i (ρ m,i -ρ i )}
where ϕ m,i is the maximum flow allowed by the capacity of cell i, ρ m,i is the jam density (i.e. the maximum density that can be reached), v i corresponds to the free flow speed and w i is the congestion wave speed in cell i. All these parameters can be the same for all cells or allowed to vary for each cell. It is noteworthy that D i-1 is the flow that can be delivered by the cell i -1 while S i is the flow that can be received by the cell i. Definition 1: A cell i ∈ {2, . . . , N } is considered free if it is able to accept the flow delivered by its upstream neighboring cell (i.e. ϕ i = D i-1 ) and is considered congested if it is not free (i.e. ϕ i = S i ). The first cell is free if ϕ 1 ≤ S 1 and is congested otherwise.

The state of the system is given by the vector ρ = (ρ 1 , . . . , ρ N ), the measured data used by the system are the upstream and downstream flows (ϕ u , ϕ d ). In order to simplify the analysis we consider that only one congestion wave may exist in the highway segment. Thus one can have only N + 1 modes since the congestions always appear at cell N and propagate upstream. Furthermore, the front wave moves downstream when the congestion disappears. Denoting by F the free state of a cell and by C the congested one, two adjacent cells can be in one of the following situations: FF, FC or CC. Let us introduce the index s(k) ∈ {0, 1, . . . , N } in order to precise the mode of the entire highway segment. This index roughly locates the congestion front. Precisely, s(k) = i ∈ {0, 1, . . . , N } if and only if the first i cells are free while the last Ni are congested (see Table I for illustration).

s(k) Cell 1 Cell 2 . . . Cell N 0 C C . . . C 1 F C . . . C 2 F F . . . C . . . . . . . . . . . . . . . N F F . .

. F TABLE I OPERATING MODE TABLE

With this notation the system dynamics writes as:

     ρ(k + 1) = A s(k) ρ(k) + Bϕ(k) + B m,s(k) ρ m s(k + 1) = s(k) + f (ρ(k), ϕ(k)) y(k) = h(ρ(k)) (2) where ϕ = (ϕ u , ϕ d ) is the input, ρ m = (ρ m,1 , . . . , ρ m,N ), h(ρ(k)) =      w 1 (ρ m,1 -ρ 1 (k)), if s(k) = 0 v N ρ N (k), if s(k) = N 0, otherwise (3) and f (ρ(k), u(k)) =      -1 if C -(ρ(k), s(k)) 0 if C 0 (ρ(k), s(k), ϕ(k)) 1 if C + (ρ(k), s(k)) (4) with C -(ρ(k), s(k)) = (s(k) > 0)∧ v s(k)-1 ρ s(k)-1 (k) > w s(k) (ρ m,s(k) -ρ s(k) (k)) C 0 (ρ(k), s(k), ϕ(k)) = (s(k) = 0) ∧ ϕ u (k) = w 1 (ρ m,1 -ρ 1 (k)) ∨ (s(k) = N ) ∧ v N -1 ρ N -1 (k) ≤ w N (ρ m,N -ρ N (k)) ∨ (0 < s(k) < N )∧ v s(k)-1 ρ s(k)-1 (k) ≤ w s(k) (ρ m,s(k) -ρ s(k) (k)) ∧ v s(k) ρ s(k) (k) ≥ w s(k)+1 (ρ m,s(k)+1 -ρ s(k)+1 (k)) C + (ρ(k), s(k)) = (s(k) < N )∧ v s(k) ρ s(k) (k) < w s(k)+1 (ρ m,s(k)+1 -ρ s(k)+1 (k))
It is worth noting here that the function f (ρ(k), ϕ(k)) formalize the conditions characterizing the forward/backward motion of the congestion front. Precisely, the cell s(k) becomes congested in the moment when the conditions C -(ρ(k), s(k)) hold true. The cell s(k) + 1 becomes free when the conditions C + (ρ(k), s(k)) hold true. When the conditions C 0 (ρ(k), s(k)) are verified the front of congestion is kept inside the cell s(k) sufficiently far from the interface between cells s(k) -1 and s(k). All these conditions are based on the interface flows adjoint to the cell s(k). In order to define the matrices

A i ∈ R N ×N , B i ∈ R N ×2 , B m,i ∈ R N ×N , ∀i ∈ {0, 1, . . . , N } used in (2)
we introduce the following notation:

Γ i :=      1 -T l1 v 1 0 . . . 0 0 T l2 v 1 1 -T l2 v 2 . . . 0 0 . . . . . . . . . . . . . . . 0 0 . . . T li v i-1 1      ∈ R i×i , ∀i = 1, . . . , N
For all i = 1, . . . , N -1 one defines the matrix

∆ i ∈ R (N -i)×(N -i) by         1 -T l i+1 wi+1 T l i+1 wi+2 0 . . . 0 0 1 -T l i+2 wi+2 T l i+2 wi+3 . . . 0 . . . . . . . . . . . . . . . 0 0 0 . . . T l N -1 wN 0 0 0 . . . 1 -T l N wN                            A N = Γ N , A 0 = 1 T l1 w 2 , 0 N -2 0 N -1 ∆ 1 A i =   Γ i 0 i-1,N -i ( T li w i+1 , 0 N -i-1 ) 0 N -i,i ∆ i   , ∀ 1 ≤ i ≤ N -1 B =      T l1 0 0 0 . . . . . . 0 -T l N      ,                B m,N = 0 N,N , B m,0 = I N -A 0 B m,i =   0 i-1,i 0 i-1,N -i 0 i (-T li w i+1 , 0 N -i-1 ) 0 N -i,i I N -i -∆ i   , ∀ 1 ≤ i ≤ N -1
where 0 represents either a column vector (one index) or a matrix (two indices) having all the components equal zero.

III. OBSERVABILITY

The observability for different SMM modes (see [START_REF] Muñoz | Traffic density estimation with the cell transmission model[END_REF] and the references therein) is summarized in Table II on this table we define the notions of forward/backward observability. In order to better explain these concepts let us consider a simple case study where the highway segment has only one cell.

ρ(k + 1) =              ρ(k) + T l ϕ u (k) -vρ(k) , if FF ρ(k) + T l w(ρ m -ρ(k)) -ϕ d (k) , if CC ρ(k) + T l ϕ u (k) -ϕ d (k) , if FC y(k) = h(ρ(k))
where FF means the upstream and the downstream flows are free, CC means the upstream and the downstream flows are congested while FC means the upstream flow is free and the downstream flow is congested.

In FF case ϕ u is not restricted by the density of the cell while ϕ d = vρ which is state dependent. In other words ϕ u is the input and ϕ d the output of the system. The observability matrix is equal with v. Since v = 0, the observability condition is satisfied and we say the system is backward observable (i.e. using downstream measurements).

In CC case ϕ d is measured and ϕ u = w(ρ m -ρ). Thus we have a reversed situation in which ϕ d is the input and ϕ u is the output of the system. In this case the observability matrix is w. Since w = 0, the observability condition is satisfied and we say the system is forward observable (i.e. using upstream measurements). Finally in FC case neither ϕ u , nor ϕ d depend on the density of the cell. The observability matrix is 0 and the density is unobservable but as we shall see the system is open-loop stable.

IV. OBSERVER DESIGN

An open-loop estimator for system (2) should be described by

ρ(k + 1) = A ŝ(k) ρ(k) + Bϕ(k) + B m,ŝ(k) ρ m ŝ(k + 1) = ŝ(k) + f (ρ(k), ϕ(k)) (5)
But, computations similar with those provided in the proof of Proposition 3 below, show that using this scheme the estimation error will remain always constant. Therefore, if the initial estimation is bad the algorithm is worthless.

In the next section we propose an observer that use the measures (v u , v d ) to detect the moment when the congestion front pases over the sensor. Doing so, we design the input and the output of the system according to the operating mode in order to assure the estimation error decreasing. The observer design is illustrated in figure 2 (see also [START_REF] Kurzhanskiy | Set-valued estimation of freeway traffic density[END_REF] below).

Input/Output selector

Mode selector Observer

ϕ d (k) ϕ u (k) u(k) y(k) ρ(k) ŝ(k)
Fig. 2. The proposed observer architecture.

V. TRAFFIC STATE ESTIMATION A. Density estimation for a highway segment

In the sequel we consider a highway segment partitioned in N cells. We suppose that all the cells have the same length (i.e. l i = l). We also consider the estimation dynamics

         ρ(k + 1) = A ŝ(k) ρ(k) + B φ(k) + B m,ŝ(k) ρ m + L ŝ(k) (y(k) -ŷ(k)) ŝ(k + 1) = ŝ(k) + f (ρ(k), φ(k)) ŷ(k) = h(ρ(k)) (6)
where φ will be defined according to different possible situations. As in the single cell case we study three situations. Situation 1: backward observer for the FF mode. The upstream and the downstream flows are free. In this case ϕ u is measured and not restricted by the density of the first cell while ϕ d = v N ρ N . Thus, the density dynamics is given by:

ρ(k + 1) = A N ρ(k) + Bϕ(k) + B m,N ρ m (7)
Let us consider φ(k) = (ϕ u , v N ρN ) and the initial estimation ρ(0) = (0, 0, . . . , 0) .

Proposition 1: If the upstream and the downstream flows are free the open-loop estimation error is exponentially decreasing and the decreasing rate is given by

max i=1,n 1 - T l i v i .
Proof: As before the estimation error at step k will be denoted by ρ(k) := ρ(k) -ρ(k). Straightforward computation shows that

ρ(k + 1) = E 1 ρ(k), E 1      1 -T l1 v 1 0 . . . 0 0 T l2 v 1 1 -T l2 v 2 . . . 0 0 . . . . . . . . . . . . . . . 0 0 . . . T l N v N -1 1 -T l N v N     
The eigenvalues of E 1 , and of E 1 as well, are 1 -T li v i , i ∈ {1, . . . , N }. Therefore, denoting by ||E 1 || the spectral norm of E 1 and by |x| the Euclidean norm of the vector x, one obtains:

|ρ(k + 1)| ≤ ||E 1 |||ρ(k)| = λ max (E 1 E 1 )|ρ(k)| ≤ max i=1,n 1 - T l i v i |ρ(k)|
where λ max (E 1 E 1 ) stands for the largest eigenvalue of the symmetric matrix

E 1 E 1 . Remark 1: Using the definition of h(•) (see (3)), when s(k) = N one gets y(k) -ŷ(k) = v N ρ(k).
Thus, in the FF case the closed-loop error dynamics is given by

ρ(k + 1) = (E 1 -L F • C F )ρ(k)
where L F = L N v N := ( F,1 , F,2 , . . . , F,N ) and C F = (0, . . . , 0, 1). The eigenvalues and the spectral norm of E 1 -

L F • C F =      1 -T l1 v 1 0 . . . 0 -F,1 T l2 v 1 1 -T l2 v 2 . . . 0 -F,2 . . . . . . . . . . . . . . . 0 0 . . . T l N v N -1 1 -T l N v N -F,N     
can be arbitrarily decreased by choosing an appropriate observer gain L F .

It is worth noting that E 1 is a non-negative matrix and the choice of ρ(0) = (0, 0, . . . , 0) will assure the nonnegativity of the vector ρ. Therefore, we always underestimate the traffic state in the FF mode. Situation 2: forward observer for the CC mode. The upstream and the downstream flows are congested. In this case ϕ d is measured and does not depend on ρ N . On the other hand ϕ u = w 1 (ρ m,1 -ρ 1 ). The density dynamics is given by:

ρ(k + 1) = A 0 ρ(k) + Bϕ(k) + B m,0 ρ m (8)
In this situation we consider φ(k) = (w 1 (ρ m,1 -ρ1 ), ϕ d ) and the initial estimation ρ(0) = ρ m .

Proposition 2: If the upstream and the downstream flows are congested the open-loop estimation error is exponentially decreasing and the decreasing rate is given by

max i=1,n 1 - T l i w i .
Proof: The dynamics of the estimation error is given in this case by

ρ(k + 1) = E 2 ρ(k), E 2         1 -T l1 w 1 T l1 w 2 0 . . . 0 0 1 -T l2 w 2 T l2 w 3 . . . 0 . . . . . . . . . . . . . . . 0 0 0 . . . T l N -1 w N 0 0 0 . . . 1 -T l N w N        
The eigenvalues of E 2 , and of E 2 as well, are 1 -T li w i , i ∈ {1, . . . , N }. Denoting by ||E 2 || the spectral norm of E 2 one obtains:

|ρ(k + 1)| ≤ ||E 2 |||ρ(k)| = λ max (E 2 E 2 )|ρ(k)| ≤ max i=1,n 1 - T l i w i |ρ(k)|
where as before λ max (E 2 E 2 ) stands for the largest eigenvalue of the symmetric matrix E 2 E 2 .

Remark 2: Using the definition of h(•) (see (3)), when s(k) = 0 one obtains y(k) -ŷ(k) = -w 1 ρ(k). Thus, in the CC case the closed-loop error dynamics is given by

ρ(k + 1) = (E 2 -L C • C C ) ρ(k)
where L C = -L 0 w 1 := ( C,1 , C,2 , . . . , C,N ) and C C = (1, 0, . . . , 0). The eigenvalues and the spectral norm of

E 2 - L C • C C =         1 -T l1 w 1 -C,1 T l1 w 2 0 . . . 0 -C,2 1 -T l2 w 2 T l2 w 3 . . . 0 . . . . . . . . . . . . . . . -C,N -1 0 0 . . . T l N -1 w N -C,N 0 0 . . . 1 -T l N w N        
can be arbitrarily decreased by choosing an appropriate observer gain L C .

It is worth noting that E 2 is a non-negative matrix and the choice of ρ(0) = ρ m will assure the non-positivity of the vector ρ. Therefore, we always overestimate the traffic state in congested flow mode. Situation 3: Open-loop estimator for the coupled FC mode. The upstream flow is free and the downstream flow is congested. In this case neither ϕ u , nor ϕ d depend on the density of the first cell or last cell. Furthermore, in this situation we cannot detect how many cells are congested.

As even in a single cell case we are not able to provide an algorithm for the density estimation, we show that (6) keeps the estimation error constant during this transition phase between Situation 1 and Situation 2.

Remark 3: If we initialize the system when both the upstream and downstream flows are either free or congested, the density of each cell can be very well approximated before encountering the situation 3. On the other hand, if the length of the highway segment is not very large we can propagate the estimation error for a short period before switching to Situation 1 or Situation 2 and starting to decrease it. In this situation we consider φ(k) = ϕ(k) and we initialize with 0 the densities of the cells assumed free and with the corresponding component of the vector ρ m the densities of the cells assumed congested. As we shall see our assumption on the number of free/congested has no impact on the state estimation algorithm.

Proposition 3: If the upstream flow is free and the downstream flow is congested the estimation error provided by the equations ( 2) and ( 6) does not increase.

Proof: For all i belonging to {1, . . . , N } let us introduce the following quantity:

K i (k) = T l v i ρ i (k) -w i+1 (ρ m,i+1 -ρ i+1 (k)) (9)
Without any loss of generality let us suppose that s(k) = i ∈ {1, . . . , N } and ŝ(k) = i + j ∈ {1, . . . , N }, j ≥ 0. Therefore, (2) and ( 6) lead to

ρ(k + 1) = A i ρ(k) + Bϕ(k) + B m,i ρ m ρ(k + 1) = A i+j ρ(k) + B φ(k) + B m,i+j ρ m Therefore, since B i = B i+j and ϕ(k) = φ(k) one obtains ρ(k + 1) = A i ρ(k) -A i+j ρ(k) + B m,i ρ m -B m,i+j ρ m = A i+j ρ(k) + j =1 (A i+ -1 -A i+ )ρ(k)+ + j =1 (B m,i+ -1 -B m,i+ )ρ m = A i+j ρ(k) +              0 i-2 K i-1 (k) K i (k) -K i-1 (k) K i+1 (k) -K i (k) . . . K i+j (k) -K i+j-1 (k) -K i+j (k) 0 N -i-j             
Thus, Sum(ρ(k+1)) = Sum(A i+j ρ(k)). On the other hand

Sum(A i x) = Sum(Γ i (x 1 , . . . , x i ) ) + T l w i+1 x i+1 + Sum(∆ i (x i+1 , . . . , x N ) ) = Sum((x 1 , . . . , x i ) ) + Sum((x i+1 , . . . , x N ) ) = Sum(x), ∀x ∈ R n , ∀i ∈ {1, . . . , N }
We conclude that Sum(ρ(k + 1)) = Sum(ρ(k)) which means that the estimation error is constant in average inside the highway segment but it may be distributed in different way at each time-step. Precisely, the estimation error will decrease inside the cells where both the inflow and outflow are either free or congested for the real and the estimation model in the same time and it will accumulate in the other cells.

B. Global observability for traffic state

From the previous sections it is clear that the hybrid system (2) is not observable since some of its mode are unobservable. Nevertheless, during the unobservable mode the state estimation error remains constant. This means that using only partial data we are able to asymptotically reconstruct the state of the system if the following Assumption is satisfied.

Assumption 1: The coupled FC mode periods are shorter than δ. It is noteworthy that Assumption 1 is satisfied in practice and δ can be fixed by studying the behavior of the network during several weeks.

Definition 2: We say that a system is globally observable if there exists a non-increasing function V : R → R + characterizing the error estimation and some fixed strictly positive constants δ ∈ Z and α < 1 such that V (k + pδ) < α p V (k), ∀p ∈ Z.

Proposition 4: If Assumption 1 holds, the system (2) is globally observable using the observer [START_REF] Kurzhanskiy | Set-valued estimation of freeway traffic density[END_REF].

Proof:

Let us consider V (k) = |ρ(k)| and α = max {||E 1 -LC N ||, ||E 2 -LC 0 ||} < 1.
Since we are not able to assure the decreasing of V during the coupled FC mode periods we shall consider the observer gain L = 0 during these periods. Doing so we get V (k + 1) ≤ V (k) for all k ≥ 0. On the other hand Assumption 1 assures as that during δ steps at least one time the forward or the backward observability situation is encountered. Taking into account the definition of α one obtains V (k + δ) < αV (k) which leads straightforwardly to the global observability of the system (2).

VI. SIMULATION RESULTS

In the sequel, the theoretical results are illustrated by some simulations. The fundamental diagram of the network has been done using the technique described in [2]. Precisely, we consider a highway segment with five identical cells. The jam density is ρ m,i = 200 V eh/Km, the free flow speed is v i = 90 Km/h and the front of congestion speed is w i = 16 Km/h.

The macroscopic simulation has been done during 140 minutes. The upstream flow was sequentially increased while the downstream flow was set constant to a half of the maximal capacity of the road. This has been done in order to create a congestion which expand backwards. When all the cells have been congested we have drastically decreased the upstream flow inducing the congestion vanishing. It is worth noting here that the macroscopic simulation accurately reproduce the measured densities(see Figure 3).

Figure 4 emphasize the behavior of the network and the congestion front motion. Figure 5 shows the estimation results when the estimation starts after 50 minutes when the fifth cell is already congested. Therefore, we have initialized the densities of the first four cells to zero and the density of the last cell to the jam density ρ m,5 = 200 V eh/Km. The estimated densities approach the simulated ones during the CC mode.

VII. CONCLUSIONS

In this paper we proposed a strategy for real-time density estimation for traffic networks. To this aim, we introduced a deterministic constrained macroscopic model which reduce the number of possible affine dynamics of the system and preserve the number of vehicles in the network. This model is used to recover the state of the traffic network and precisely localize the eventual congestion front. The state of the network is recovered using what we call forward/backward observers. We pointed out that during unobservable modes the estimation error is preserved due to vehicle conservation law. Numerical simulations show the efficiency of the proposed strategy.

Fig. 1 .

 1 Fig. 1. The switching mode model based on cell transition model.

4 Fig. 3 .

 43 Fig.3. The density evolution given by the macroscopic simulation is smoother but it accurately reproduces the behavior of the measured density.

5 Fig. 4 .

 54 Fig. 4. The congestion front appears downstream and propagates backward from the fifth cell to the first one.

  . Based

	Upstream Cells Downstream Cells	Observable with
	Free-flow	Free-flow	Downstream measurement
	Congested	Congested	Upstream measurement
	Congested	Free-flow	Up. and Down. measurement
	Free-flow	Congested	Unobservable
		TABLE II
	OBSERVABILITY FOR DIFFERENT SMM MODES

This work has been done during the post-doctoral fellowship of I.-C.