Prof François Bouchut

Christiane Lafargue

Anh Đức

Ơn Việt

Sao

Subbie Bích

Ngô

Chú Thuận

Hân

Em Như Sơn

Thúy Giang

Lời Cuối

for his continuous support of my PhD study and research, particularly for his patience, motivation, enthusiasm, and encouragement

Constitutive rheological laws for granular materials

Classical Physics contain several very well-known equations called "constitutive equations": Hooke's law, Fick's law, Ohm's law... These laws do not illustrate directly fundamental laws of nature, but rather a relation between two physical quantities, that should hold for any steady or dynamical evolution. In fluid dynamics, a crucial constitutive relation characterizing a material is called "Rheology". It means merely relations between stress and strain rate, and eventually other quantities. The history of this issue traces back to the 17th century. Sir Isaac Newton observed that each fluid has its own resistance to deformation at a given rate. It usually corresponds to the informal conception of the "thickness" of a fluid. For example, when pouring out a can of honey and another one containing water, it can be observed that flow of water can easily deform to a greater extent than the flow of honey does; thus one has tendency to conclude that Honey is "thicker" than water. To debunk this misconception, the quantity "viscosity" was created. Newton claimed that under certain ideal conditions each fluid has its own constant viscosity, for example, the statement "honey has a higher viscosity than water" is considered to be more precise. Now, if the fluid is examined under the assumption that its viscosity remains constant then it can be said intuitively that the flow with the higher velocity is "stronger". Scientifically one concludes that with higher strain rate (deformation over time), the fluids possesses higher stress tensor (internal forces that neighbouring particles of a continuous material exert on each other). There are fluids such as water, air, alcohol, glycerol, and thin motor oil, etc.. that have a linear stress-strain rate relation. Such single-phase fluids made up of simple, small-weight molecules are called "Newtonian fluids". Mathematically, the rheological relation can be formulated for Newtonian fluids as σ = ηDu, where σ denotes the stress tensor and Du is the strain rate. However most fluids in the real world are Non-Newtonian fluids. Roughly speaking this means that the fluid has a shear stress nonlinearly proportional to the deformation of the media. A quite general constitutive equation can then be formulated as σ = F (Du) where F is a nonlinear function.

Bingham fluid

Even in the world of Non-Newtonian fluids, there are certain classifications based on the nonlinearity of the stress -strain rate relation, such as: Fluids that exhibit a logarithmic increase in shear stress with shear strain rate are called shear-thinning fluids or pseudo-plastic fluids (i.e. blood, paint, ketchup). Meanwhile for other types of fluids which are called shear-thickening fluids In 1913, [START_REF] Bingham | Fluidity and plasticity[END_REF][START_REF] Bingham | An investigation of the laws of plastic flow[END_REF] a professor of the department of Chemistry at Lafayette College, named Eugene Cook Bingham introduced the first mathematical form describing the motion of a fluid that behaves as a rigid body at low stress but flows as a viscous fluid at high stress. This kind of Non-Newtonian fluid, also called "fluid by threshold", is named after him Bingham fluid, and constitutes a very important class of Non-Newtonian fluids. The study of Bingham materials is a crucial branch exhibiting the fundamental character of viscoplasticity. That is the reason why the simulation of Bingham fluids is always necessary and plays a key role in fully understanding the viscoplastic material; however, the work of modeling all the properties of viscoplastic media can not rely solely on the study of the Bingham model. Besides, modeling and simulating all the properties of viscoplastic media seems to be an impossible task. In general, besides the Non-Newtonian behaviour, this thesis will focus on the fact that viscoplastic fluids are also characterized by irreversible deformations inducing dissipation.

Why studying viscoplastic materials is important? Viscoplastic media are present in our real life much more frequently than it is thought. An impressive example which proves that people cope with viscoplastic fluids everyday without even being aware of: Every morning, say you want to get some toothpaste out of the tube; when placed merely upside down, the paste will not come out even under the sole influence of the gravity or when squeezed with insufficient force. At this point, the toothpaste is considered to have an infinite viscosity and acts as a solid body. On the other hand, when the force is exerted beyond the certain threshold, the toothpaste will become runnier and will behave more "liquidy". There are numerous examples of viscoplastic media, namely butter, foams, pastes, slurries, oils, ceramics, etc. The study upon viscoplastic fluid has a vast number of applications in: the prediction of the plastic Figure 1.2: Snow avalanche, pyroclast are also considered as viscoplastic media (Source: Internet) collapse of structures, dynamic systems exposed to high strain rates, crash simulations, etc. We also can find their application in the petroleum industry: mitigation of paraffin wax, cements or drilling mud from crude oil. Moreover, a large variety of products in food industry (milk products, jam, chocolate confections, etc.) exhibit the fundamental characteristics of viscoplasticity. Their applications are also found in the study of environmental disasters (as in Figure 1.2): avalanche, ejection of volcanic magma, mudslides, and so forth.

This thesis is indeed motivated by the application to granular materials that are involved in the modeling of landslides. Their description by appropriate viscoplastic models is the subject of a research program between the Laboratoire d'Analyse et de Mathématiques Appliquées (Université Gustave Eiffel) and the Institut de Physique du Globe de Paris. A first relevant rheology for granular materials has been proposed by [START_REF] Jop | A constitutive law for dense granular flows[END_REF], the so called µ(I) rheology. This model does however not well represent real granular flows, even in experimental laboratory context. Additionally it has been proved by [START_REF] Barker | Well-posed and ill-posed behaviour of the µ(I)-rheology for granular flow[END_REF] that this model is strongly ill-posed. In numerical computations it can show shear band instabilities that are stronger and stronger when refining the mesh. The original µ(I) law being an incompressible model, it has been then recently shown [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF] that considering compressible models could be a way to bypass the ill-posedness. Moreover volume variations, that do not exceed 10% in dry materials, are suspected to play nevertheless a key role in the dynamics. Moreover in the presence of an interstitial fluid like water (wet flows), volume and pore pressure variations are ubiquitous and their mathematical description are absolutely necessary. This is what is called "dilatancy", and a relevant law has been proposed in [START_REF] Roux | Texture-dependent rigid plastic behaviour[END_REF], that describes how the rheology should involve the divergence of the velocity field. Hence it is necessary to build models and develop simulations tools for viscoplastic compressible models with variable volume fraction. These models must include a yield stress to simulate the static and flowing parts of the material.

Progress in the mathematical study and numerical simulation of viscoplastic media

The rigorous analysis of viscoplastic models is at the moment mostly restricted to Bingham's law or direct generalizations. Over the past 40 years most studies assumed rheological models with a threshold of the stress tensor of Bingham or Herchel-Buckley type [START_REF] Herschel | Measurement of consistency as applied to rubber-benzene solutions[END_REF], where the viscosity follows a power law, and depends only on the norm of the strain rate. Mathematical studies of compressible viscoplastic models can be found in [START_REF] Málek | Compressible generalized Newtonian fluids[END_REF][START_REF] Sadovskaya | Mathematical modeling in mechanics of granular materials[END_REF][START_REF] Han | Numerical analysis of the primal problem[END_REF][START_REF] Ae Mamontov | Existence of global solutions to multidimensional equations for Bingham fluids[END_REF]. But most viscoplastic fluids are not like this, a dependency in the trace of the strain rate is important to include so called dilatancy effects, which are not present in most studies that consider incompressible models. The subject of this thesis is to discuss more complex models with nonlinearity that can include dilatancy effects. In particular the pure plastic models are part of it: Cam-Clay model, Drucker dilatant model or degenerate Bingham model. There is a theoretical issue in particular concerning the existence of solutions to these models and their characterization.

Numerically simulating viscoplastic flows is not straightforward. The difficulty is due to the presence of unknown interfaces separating the yielded and the unyielded regions, which are difficult to track. The problem can be written as a set of nonlinear variational inequalities. The case where it is well studied case is the Bingham case with viscosity [START_REF] Duvaut | Les inéquations en mécanique et en physique. dunod, paris[END_REF][START_REF] Glowinski | Handbook of numerical analysis[END_REF]. The most well-known method is the regularization method [START_REF] Temam | Problèmes mathématiques en plasticité[END_REF], but the augmented lagrangian method is also very much used [START_REF] Roquet | An adaptive finite element method for Bingham fluid flows around a cylinder[END_REF]. These methods have not been extended to general nonlinear laws, and moreover the case without viscosity has not been much considered except [START_REF] Bouchut | Convergence of conforming approximations for inviscid incompressible Bingham fluid flows and related problems[END_REF].

Mathematical description of viscoplastic materials 1.2.1 Conservation of mass

Consider Ω ⊂ R N where N = 1, 2, 3 an open, bounded domain occupied by the medium. The conservation of mass expresses that the rate of increase of mass equals to the mass influx through the boundary of a control volume. Writing this law upon the control volume V ⊂ Ω which is an arbitrary open, regular subdomain, this gives

d dt V ρ(t, x) dx = - ∂V ρ(t, x)u(t, x) • n(x) ds,
where n is the outward unit normal vector. The sign of the normal component of the velocity determines whether the fluid flows goes in or out of the control volume. Applying the Stokes formula this yields the concise form V ∂ρ ∂t + div(ρu) dx = 0.

Since this relation can be applied for all arbitrary neighborhood of a point x ∈ Ω and for any t ∈ (0, T ), we obtain the local form of the conservation of mass in the differential form (known as the continuity equation)

∂ρ ∂t + div(ρu) = 0 in (0, T ) × Ω , (1.2.1)
where ρ = ρ(t, x) is the density of the flow, and u = u(t, x) = (u i (t, x)) 1≤i≤N is the velocity vector field at point x and time t.

Conservation of momentum

The second relation is based on one of the most well-known law in physics: Newton's third law: the rate of increase of momentum equals to the net influx of momentum plus the exerted force. By recalling the unit momentum per volume as ρu, we have the momentum in any control volume V would be V ρu dx.

The net inflow momentum can be written according to Stokes' formula

- ∂V ρu(u • n) ds = V div(ρu ⊗ u) dx .
The forces exerted upon the control volume V can be classified in two types of forces

• External force, in particular the total body force due to gravity. It can be formulated as

V ρ g,
where g is the gravity vector that has constant magnitude.

• Internal force, caused by the deformation of the fluid, that can be written as

∂V σ tot n(x) ds,
where σ tot (t, x) is the symmetric tensor of total stress.

Finally we obtain the following formula which expresses the conservation of momentum

V ∂(ρu) ∂t + div(ρu ⊗ u) -div σ tot dx = V ρ g dx .
The local form of the conservation of momentum is therefore ∂(ρu) ∂t + div(ρu ⊗ u) -div σ tot = ρ g in (0, T ) × Ω.

(1.2.2)

We shall assume from now on that the total stress can be decomposed as the sum of a compressible pressure law and a material dependent stress:

σ tot = -p th (ρ) Id +σ.

(1.2.3)

Constitutive rheological law

The aim of this subsection is to discuss the heart of our subject: the law for σ. Let us start by some familiar notions in matrix theory. The trace of a matrix is defined as

Tr(σ) := N i=1 σ ii .
Let us denote by Id := (δ ij ) 1≤i,j≤N the identity matrix, where δ ij is the Kronecker symbol. Then we can always decompose a (symmetric) tensor σ as

σ = 1 N Tr(σ) Id +σ , (1.2.4)
where σ (the "deviatoric" part of σ) satisfies Tr(σ ) = 0. We then introduce the pressure p by Let us remark that there is a general difficulty in defining the notion of pressure in fluids mechanics. It is surprising that such a classical notion remains not fully clear. Yet there is no concrete argument for formulating the pressure p as (1.2.5). One may wonder: Why do we have to decompose the diagonal part of the stress tensor? At least we can remark that if we define the inner product in the space of the tensors as

σ : τ = N i,j=1 σ ij τ ij ,
then it follows that the decomposition (1.2.4) is orthogonal. The associated Fröbenius norm |σ| = (σ : σ) 1/2 has the property to be a physical invariant: it can be computed with the same formula in any orthonormal basis.

The constitutive law expresses the relationship between the stress tensor σ and the strain rate tensor

Du = ∇u + (∇u) t 2 .
(1.2.7)

Let us denote by M s N ×N (R) the space of real symmetric matrices of size N × N . Suppose that F : M s N ×N (R) → R is convex, lower semi-continuous, and proper, where R := R ∪ {+∞}. Let us recall that F is called proper if it is not identically +∞ i.e. there exists at least one D such that F (D) < +∞. We are going to consider rheologies of viscoplastic type defined by a relation σ ∈ ∂F (Du), (1.2.8) where ∂F denotes the subdifferential, which is defined below. Knowing the monotonicity of the subdifferential of a convex function, (1.2.8) characterizes the monotonicity of the stress-strain rate relation. As we shall see this is an important property that ensures the existence of a variational formulation for the momentum evolution (1.2.2) (at least without inertial terms).

Since the above assumptions on F will be repeated many times, let us name them as Hypothesis 1. F : M s N ×N (R) → R is a convex, proper and lower semi-continuous function.

Due to the lack of regularity for F , the notion sub-gradient/subdifferential ∂F which generalizes the classical gradient/differential, plays a crucial role. We denote by dom(F ) the set of Ds that satisfy F (D) < ∞. Definition 1.2.1 (Subgradient). Suppose that F satisfies Hypothesis 1. A symmetric tensor σ is called a subgradient of F at D if

F (D) ≥ F (D) + σ : (D -D), ∀D ∈ M s N ×N (R).
Definition 1.2.2 (Subdifferential). Suppose that F satisfies Hypothesis 1. The subdifferential ∂F (D) of F at D is the set of all subgradients:

∂F (D) = σ ∈ M s N ×N (R) | F (D) ≥ F (D) + σ : (D -D) ∀D ∈ M s N ×N (R) .
Taking into account this definition in the constitutive relation (1.2.8), we see that we obtain a set of inequalities, a so called variational formulation (here in a local setting, i.e. not integrated in space).

Let us introduce some well-known viscoplastic models that enter the considered class, i.e. they can be written as (1.2.8). Examples

• Bingham model

It takes the form (1.2.8) with

F (D) = η |D| 2 2 + σ 0 |D| ∀D ∈ M s N ×N (R), (1.2.9)
where η := η(ρ) > 0 represents the viscosity of the fluid, and σ 0 := σ 0 (ρ) > 0 is the yield stress. In other words the relation between σ and D is

   σ = η + σ 0 |D| D if D = 0, |σ| ≤ σ 0 if D = 0.
(1.2.10)

Thus the material is solid (D = 0) for |σ| ≤ σ 0 , and liquid (D = 0) for |σ| > σ 0 . The case without viscosity η = 0 is special: σ depends only on the direction of D, not on its magnitude. In this case it is a pure plastic model.

• Herschel-Bulkley model

In this model F is taken as

F (D) = η 1 + d |D| 1+d + σ 0 |D|,
where η > 0, d > 0. Taking the subdifferential of F gives

σ = (η|D| d-1 + σ 0 |D| -1 )D if D = 0, |σ| ≤ σ 0 if D = 0.
Again σ 0 is the yield stress, η is the consistency parameter, d > 0 is the flow index. Roughly speaking, for d < 1 the fluid exhibits shear-thinning properties whereas for d > 1 it is shearthickening. The case d = 1 corresponds to the Bingham model.

• Newtonian model

The simplest case in the Herschel-Bulkley model is to take n = 1 and σ 0 = 0, i.e. F (D) = η|D| 2 /2. The model then reduces to the well-known Newtonian model σ = ηD (Navier-Stokes equations).

The previous models are rather simple: F depends only on |D|. There exist more elaborate models which are relevant for granular media: Cam-Clay model, Drucker-dilatant model, degenerate Bingham model. These will be discussed in Chapter

p = f 1 (Tr(D), |D |), |σ | = f 2 (Tr(D), |D |), (1.2.12)
with f 1 , f 2 two scalar functions. The mechanical description of rheology is more often written as

|σ | = g 1 (p, |D |), Tr(D) = g 2 (p, |D |), (1.2.13)
with g 1 , g 2 two scalar functions that represent respectively the internal friction law and the dilatancy law. It is not so easy to pass from (1.2.13) to (1.2.12).

Initial and boundary conditions

The problem (1.2.1), (1.2.2) is completed by initial conditions on the density and the velocity field,

ρ(t = 0, x) = ρ 0 (x), x ∈ Ω, u(t = 0, x) = u 0 (x), x ∈ Ω. (1.2.14)
One has also to set conditions on the boundary of Ω. There are several boundary conditions that are physically relevant: Dirichlet, Neumann, periodic, slip, friction conditions. To give an example, let us consider that the boundary is decomposed as several parts and impose a different condition on each of them, see Figure 1.3. It can be

Dirichlet u = 0 in (0, T ) × Γ D , Neumann σn = 0 in (0, T ) × Γ N , (1.2.15) slip u • n = 0 and (σn) × n = 0 in (0, T ) × Γ S ,
where ∂Ω = Γ D ∪ Γ N ∪ Γ S . For applications to granular flows it is useful to have a moving free boundary. In this thesis we shall not consider this case.

Incompressible and compressible models

One states that "The compressibility (change in volume due to change in pressure) of a liquid is inversely proportional to its volumic modulus of elasticity, also known as the bulk modulus". In simple terms, the compressibility of the fluid is related to the variation of its density ρ. When the density is constant, the continuity equation (

                                   ∂ρ ∂t + div(ρu) = 0 in (0, T ) × Ω, ∂(ρu) ∂t + div(ρu ⊗ u) + ∇p th (ρ) -div σ = ρg in (0, T ) × Ω, σ ∈ ∂F (Du) in (0, T ) × Ω, u = 0 in (0, T ) × Γ D , σn = 0 in (0, T ) × Γ N , u • n = 0 and (σn) × n = 0 in (0, T ) × Γ S , ρ(t = 0, x) = ρ 0 (x), u(t = 0, x) = u 0 (x), x ∈ Ω.
(1.2.17)

where F satisfies Hypothesis 1. To simplify we can consider that p th (ρ) = κρ γ where γ > 1, a classical behaviour in gas dynamics. Physically relevant values for gases are γ = 5 3 or 7 5 . The previous system of equations is classical in fluid mechanics, except the rheological behaviour. In the case of a Newtonian material σ = ηDu, this is the compressible Euler system of gas dynamics. Existence of solutions to this system in the Newtonian case is now known in arbitrary dimension for γ not too small. However this is still mostly an open problem for non Newtonian rheologies, see however [START_REF] Feireisl | Global weak solutions to a class of non-Newtonian compressible fluids[END_REF][START_REF] Ae Mamontov | Existence of global solutions to multidimensional equations for Bingham fluids[END_REF], and [START_REF] Basov | Generalized solutions to the equations of compressible Bingham flows[END_REF] for the 1d case. Notice that the existence for the 2d incompressible inviscid Bingham including inertial terms has been established in [START_REF] Lions | Remarks on some nonlinear evolution problems arising in Bingham flows[END_REF], see also [START_REF] Chupin | Existence theorem for homogeneous incompressible Navier-Stokes equation with variable rheology[END_REF].

Numerical scheme and time algorithm

In order to solve numerically the problem (1.2.17), we use the standard time splitting algorithm. Let N t be a positive integer. We consider a constant timestep ∆t > 0 and define the discrete times by t n = n∆t, n ∈ {0, 1, . . . , N t }.

We denote by u n the approximation of the velocity field u at time t n , i.e. u n ≈ u(t n ), and similarly for the density ρ n ≈ ρ(t n ). We use a similar notation for the stress tensor σ n ≈ σ(t n ) or other useful quantities. We shall also need quantities u n+ 1 2 , ρ n+ 1 2 that are intermediate between n and n + 1. In order to update the values ρ n and u n we proceed as follows. Suppose that ρ n , u n are known. As a first step we apply the finite volume method (FVM) to solve the compressible Euler system, and get u n+ 1 2 , ρ n+ 1 2 . In other words we have a finite volume discretization of

           ρ n+ 1 2 -ρ n ∆t + div(ρ n u n ) = 0, ρ n+ 1 2 u n+ 1 2 -ρ n u n ∆t + div(ρ n u n ⊗ u n + p n th Id) = 0. (1.3.1)
In the second step we obtain u n+1 , ρ n+1 by using the finite element method (FEM) to solve

           ρ n+1 -ρ n+ 1 2 ∆t = 0, ρ n+1 u n+1 -ρ n+ 1 2 u n+ 1 2 ∆t -div σ n+1 = f n+ 1 2 , σ n+1 ∈ ∂F (Du n+1 ), (1.3.2)
with f n+1/2 = ρ n+1/2 g. Since the first equation gives trivially ρ n+1 = ρ n+1/2 , only the second equation comes into play, it is a viscoplastic problem with a space dependent weight (the density). We also need boundary conditions for each of the systems (1.3.2) and (1.3.1). We will discuss this in more detail in Chapter 5.

The viscoplastic problem

In the previous subsection we have seen that the time splitting algorithm leads to the viscoplastic problem (1.3.2). When the weight ρ n+1/2 is constant (= 1) this can be written as a steady problem

αu -div σ = f, σ ∈ ∂F (Du), (1.4.1)
with α = 1/∆t and f = f n+1/2 + u n+1/2 ∆t . This viscoplastic problem can be formally obtained by minimizing the functional

J(v) = α Ω |v| 2 2 + Ω F (Dv) - Ω f • v. (1.4.2)
This minimization can be proved to be equivalent to the variational formulation

α Ω u • (v -u) + Ω F (Dv) ≥ Ω F (Du) + Ω f • (v -u), for all functions v. (1.4.3)
This problem has been studied a lot [START_REF] Duvaut | Les inéquations en mécanique et en physique. dunod, paris[END_REF][START_REF] Temam | Problèmes mathématiques en plasticité[END_REF][START_REF] Trémolières | Numerical analysis of variational inequalities[END_REF], but usually only the Bingham model with positive viscosity is considered, or eventually the Herschel-Bulkley model. The case of inviscid Bingham has been considered in [START_REF] Bouchut | Convergence of conforming approximations for inviscid incompressible Bingham fluid flows and related problems[END_REF], where the time dependent case is also treated. In this thesis we consider a quite general F , with or without viscosity. Numerical methods to solve (1.4.1) are numerous, and we refer to Chapter 2 for the description of classical ones.

Space discretisation

As mentioned above the general scheme uses a splitting algorithm in order to solve two simpler problems, one by the FVM and the other by the FEM. We use the approach of [START_REF] Bouchut | Convection and total variation flow[END_REF][START_REF] Bouchut | Convection and total variation flowerratum and improvement[END_REF] that deals with a nonlinear hyperbolic scalar conservation law regularised by the total variation flow operator.

Here we consider a compressible fluid model together with viscoplastic rheology. First we have to recall the space discretisation used in those articles, that allows us to create two dual meshes, one for the FVM and the other for the FEM. The finite element mesh, denoted by T h , is a conforming mesh of the open, bounded domain Ω ⊂ R N of size h: T h is a finite set of disjoint open simplices such that ∪ K∈T h K = Ω, and h is the maximum value of the diameter of all K ∈ T h . The mesh is conforming in the sense that for two distinct elements K, L of T h , K ∩ L is either empty or a simplex included in an affine subset with dimension strictly lower than N , whose vertices are simultaneously vertices of K and L. The finite set of the vertices of the mesh is denoted by {x i , i ∈ V }. Each element of T h is assumed to be nonobtuse; which means that the angles between any two facets are less than or equal to π/2. For any K ∈ T h , we denote by V K ⊂ V the set of the N + 1 indices of the vertices of K, and by E K the set of the N + 1 faces of K.

The finite volume mesh, denoted by D h , is a polyhedral mesh of Ω such that the interface between two cells is a finite union of faces. The mesh D h is a dual mesh of T h in the sense that each cell of D h contains one and only one node of T h . For any i ∈ V , the cell of D h containing the node x i is denoted by

Q i . We assume that ∀i ∈ V, Q i ⊂ K∈T h ,i∈V K K.
(

Besides, N i is the set containing the indices of the neighbouring cells of Q i , E h is the set of couples (i, j) such that Q i and Q j are neighbours and i < j, Γ i,j is the interface between two neighbour cells Q i and Q j , n i,j is the unit normal vector to Γ i,j pointing toward Q j . The unknown function u(t, x) is reconstructed simultaneously from the values u h = (u i ) i∈V at the points x i for i ∈ V , using a continuous piecewise affine reconstruction denoted by ûh , and using a piecewise constant reconstruction denoted by ūh . More explicitly, for all

u h ∈ R V we have ûh ∈ C(Ω), ûh | K is affine for each cell K ∈ T h , ûh (x i ) = u i ∀i ∈ V, ūh ∈ L 1 loc (Ω), ūh (x) = u i ∀x ∈ Q i , i ∈ V.
Similarly, we use the same notation for ρh , ρh . The numerical scheme approximating (1.2.17) is defined by

1 -Initialisation of u 0 h ∈ R V , ρ 0 h ∈ R V : u 0 i = u init (x i ) ∀i ∈ V, ρ 0 i = ρ init (x i ) ∀i ∈ V.
2 -Finite volume step: Suppose that (u n h , ρ n h ) are known. We compute the numerical flux by a first-order explicit formula

F n i,j = F (U n i , U n j , n i,j ) (1.5.2)
corresponding to each pair of neighbor cells Q i , Q j , where U i = (ρ i , ρ i u i ). We assume that the scheme is conservative, i.e. F n j,i = -F n i,j . We impose the CFL condition (for Courant, Friedrichs, Levy [START_REF] Courant | On the partial difference equations of mathematical physics[END_REF]) on the timestep to prevent the blow up of the numerical values, under the form ∆ta n ij ≤ h, (1.5.3) at each interface between two cells Q i , Q j , where a ij is an approximation of the propagation speed, that should be deduced from the choice of the numerical flux. We compute ρ

n+ 1 2 h , u n+ 1 2 h by U n+ 1 2 i = U n i - ∆t |Q i | j∈N i |Γ i,j |F n i,j . (1.5.4) 3 -Finite element step: We set ρ n+1 h = ρ n+ 1 2 h
. Then, define

V h := {v h ∈ C(Ω) such that vh | K is affine for each cell K ∈ T h }, Λ h := {σ h ∈ L ∞ (Ω) such that σ h | K is constant for each K ∈ T h }.
The values ρ

n+ 1 2 h , u n+ 1 2 h
being known, we look for (û n+1 h , σ n+1 h

) ∈ V h × Λ h such that Ω ρ n+1/2 h ūn+1 h - ūn+ 1 2 h ∆t • vh dx + Ω σ n+1 h : Dv h dx = Ω fh • vh dx ∀v h ∈ V h , (1.5.5) σ n+1 h ∈ ∂F (Dû n+1 h ) a.e. in Ω. (1.5.6)

Sketch of the results

Chapter 2 describes mathematical tools that are useful to study (1.4.1), and classical numerical methods to solve this problem. A particular attention is given to the description of pure plastic models. Then Chapter 3 concerns the rigorous analysis of the theoretical problem (1.4.1). Under technical assumptions on F such as subquadratic growth and superlinearity we prove the existence of solutions to the associated variational problem (1.4.3). This is obtained in the viscous as well as in the inviscid cases. We establish Euler-Lagrange characterizations of these solutions, i.e. the equivalence with the local formulation (1.4.1). No regularity is assumed on F , thus yield stress rheologies are included. Then Chapter 4 introduces an iterative but explicit method in the sense that there is no linear system to solve, inherited from the minimizing of total variation functionals used in imaging [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. It is applicable to any type of nonlinearity, and includes a kind of projection on some convex sets, that already appears in the augmented Lagrangian method or in the time regularized approach of [START_REF] Chupin | A bi-projection method for Bingham type flows[END_REF]. We prove the convergence of the method discretized in space with finite elements. Numerical tests performed in Chapter 5 confirm the theoretical results. Finally Chapter 6 is a work on a lubrication equation for a simplified model of shear-thinning fluid, that has been performed at the CEMRACS 2019 in collaboration with Khawla Msheik, Meissa M 2019Baye and François James.

Chapter 2

Convex analysis and approximation methods for nonlinear viscoplastic models

In this section we consider the steady viscoplastic model

       σ ∈ ∂F (Du), αu -div σ = f in Ω, u = 0 on Γ, (2.0.1)
where α is a positive constant, F satisfies Hypothesis 1 and f is a given source. We are going to introduce briefly the state-of-the-art concerning the numerical resolution of (2.0.1), and introduce useful tools of convex analysis.

There are various methods concerning the numerical simulation of viscoplastic flows, in particular for the Bingham fluid when F (D) = |D|. These results have been the subject of intense investigations from the early 1970 up to now. There are mainly two historical methods to approximate viscoplastic models. The first is based on the regularization procedure, and the second is based on the use of duality formulations, as the augmented Lagrangian algorithm. Due to the impossibility to present all the methods, we shall introduce here two methods within the two above classes. Another class of methods have been formulated recently. These methods are based on the idea to formulate the stress condition in (2.0.1) via a fixed point involving a projection. This is done in [START_REF] Saramito | A damped Newton algorithm for computing viscoplastic fluid flows[END_REF][START_REF] Chupin | A bi-projection method for Bingham type flows[END_REF]. The history of the development of this kind of method is depicted in [START_REF] Saramito | Progress in numerical simulation of yield stress fluid flows[END_REF]. In Chapter 4 we shall introduce a more particular projection method that is inherited from image processing [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], that is able to deal with general nonlinearities F .

Preliminaries

Throughout the next chapters we shall often use some basic results of convex analysis and optimization. This section is dedicated to remind and list several classical results. One can easily find them in [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF][START_REF] Tyrell | Convex analysis[END_REF].

We consider a Hilbert space H, that will be taken later on as either H = M s N ×N (R) the space of symmetric square matrices of size N , or an infinite dimensional functional space. The scalar product of two elements D, σ ∈ H will be denoted by D : σ. Definition 2.1.1 (Legendre -Fenchel transform). For F : H → R proper, we can define F * : H → R, the conjugate function of F , by

F * (σ) = sup D∈H (σ : D -F (D)), σ ∈ H, (2.1.1)
where we denote R = R ∪ {+∞}.

Here the property of F to be proper means that the domain of F , i.e. the set where F has a finite value, is not empty. In other words, F is proper if and only if F is not identically +∞.

Proposition 2.1.2. Suppose that F is convex, proper and lower semi-continuous (l.s.c.). Then F * is convex, proper and l.s.c. Moreover F * * = F .

We recall the definition of the subdifferential ∂F of the function F , for any D ∈ H, (e) For σ 1 ∈ ∂F (D 1 ), σ 2 ∈ ∂F (D 2 ), we have

∂F (D) = σ ∈ H | F (D) ≥ F (D) + σ : (D -D), ∀D ∈ H . (2.1.2) Lemma 2.1.3. Let be F a proper, convex, l.s.c. function. Then (a) If ∂F (D) is nonempty then F (D) < ∞. (b) For D, σ ∈ H, σ ∈ ∂F (D) ⇐⇒ D ∈ ∂F * (σ), ( 2
(σ 1 -σ 2 ) : (D 1 -D 2 ) ≥ 0. (2.1.5) (f ) inf F = -F * (0) and inf F * = -F (0).
F * (σ) ≥ D : σ -F (D) + σ : (D -D), ∀σ ∈ H, ∀D ∈ H. (2.1.8)
Then using the definition of F * (σ) we obtain

F * (σ) ≥ F * (σ) + D : (σ -σ), ∀σ ∈ H, (2.1.9)
which means that D ∈ ∂F * (σ). The converse follows from the fact that F * * = F . Finally the properties F (D) < ∞ and F * (σ) < ∞ follow from (a). (c) The definition of F * * (D) is

F * * (D) = sup σ∈H (σ : D -F * (σ)).
To say that the sup is attained at σ means that (2.1.12)

σ : D -F * (σ) ≤ σ : D -F * (σ),
F (D) ≥ F (D 1 ) + σ 1 : (D -D 1 ), ∀D ∈ H, F (D) ≥ F (D 2 ) + σ 2 : (D -D 2 ), ∀D ∈ H. ( 2
F (D 2 ) ≥ F (D 1 ) + σ 1 : (D 2 -D 1 ), F (D 1 ) ≥ F (D 2 ) + σ 2 : (D 1 -D 2 ). ( 2
Then F ε is finite everywhere, and (a) The infimum is attained at the unique point

D = prox ε F (D), prox ε F (D) = argmin D∈H F (D) + |D -D| 2 2ε . ( 2
F ε (D) = D -prox ε F (D) ε ∈ ∂F (prox ε F (D)).
(2.1.17)

(f ) The Moreau envelope F ε is convex, and F ε (D) ↑ F (D) pointwise as ε → 0.

Proof. (a) The proof is classical and can be found in [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF][START_REF] Tyrell | Convex analysis[END_REF]. (b) Denote D = (Id +ε∂F ) -1 (D). We have

D -D ε ∈ ∂F ( D) ⇐⇒ D ∈ ∂F * D -D ε ⇐⇒ D ε ∈ 1 ε ∂F * D -D ε + D -D ε .
It follows that

D ε ∈ Id + ∂F * ε D -D ε ⇐⇒ D -D ε = Id + ∂F * ε -1 D ε .
We deduce that

D = ε Id + ∂F * ε -1 D ε + D = (Id +ε∂F ) -1 (D) + ε Id + ∂F * ε -1 D ε . Setting r = 1/ε, D = εσ, we obtain (2.1.16). (c) Denote D 1 = prox ε F (D 1 ) and D 2 = prox ε F (D 2 )
. Then from (2.1.14) one has

D 1 -D 1 ε ∈ ∂F ( D 1 ), D 2 -D 2 ε ∈ ∂F ( D 2 ).
Using the monotonicity of the subdifferential (2.1.5), one has

(D 1 -D 2 -D 1 + D 2 ) : ( D 1 -D 2 ) ≥ 0, (D 1 -D 2 ) : ( D 1 -D 2 ) ≥ | D 1 -D 2 | 2 , |D 1 -D 2 || D 1 -D 2 | ≥ | D 1 -D 2 | 2 , |D 1 -D 2 | ≥ | D 1 -D 2 |. (2.1.18) Thus |D 1 -D 2 | ≥ | prox ε F ( D 1 ) -prox ε F ( D 2 )|, (2.1.19) which proves the claim. (d) Inserting D 2 = 0, D 1 = D into (2.1.19) gives |D| ≥ | prox ε F (D) -prox ε F (0)| ≥ | prox ε F (D)| -| prox ε F (0)|, |D| + | prox ε F (0)| ≥ | prox ε F (D)|. Since prox ε F (0) = argmin D F (D) + |D| 2 2ε = C ε , we get | prox ε F (D)| ≤ C R,ε where |D| ≤ R. Thus prox ε F (D) is bounded when D is bounded. Next, considering D 1 = prox ε F (D 1 ) and D 2 = prox ε F (D 2 ), one has F ε (D 1 ) ≤ F ( D 2 ) + 1 2ε | D 2 -D 1 | 2 = F ε (D 2 ) + 1 2ε (| D 2 -D 1 | 2 -| D 2 -D 2 | 2 ) = F ε (D 2 ) + 1 2ε (D 1 + D 2 -2 D 2 ) : (D 1 -D 2 ).
Therefore

F ε (D 1 ) -F ε (D 2 ) ≤ 1 2ε (D 1 + D 2 -2 D 2 ) : (D 1 -D 2 ) ≤ CR,ε |D 1 -D 2 |,
when D 1 , D 2 are in the ball B R of center 0 and radius R. Similarly, one has

F ε (D 2 ) -F ε (D 1 ) ≤ CR,ε |D 1 -D 2 |, thus for D 1 , D 2 ∈ B R , one has |F ε (D 1 ) -F ε (D 2 )| ≤ CR,ε |D 1 -D 2 |. This proves that F ε is locally Lipschitz continuous.
(e) We want to prove that F ε is differentiable. We consider as above D 1 , D 2 ∈ H. Setting

M ε = D 1 + D 2 -2 D 1 2ε , one has F ε (D 2 ) -F ε (D 1 ) ≤ 1 2ε (D 1 + D 2 -2 D 1 ) : (D 2 -D 1 ) = (D 2 -D 1 ) : M ε . Defining L(D) = D -prox ε F (D) ε = D -D ε , one has F ε (D 2 ) -F ε (D 1 ) -L(D 1 ) : (D 2 -D 1 ) ≤ (M ε -L(D 1 )) : (D 2 -D 1 ) = |D 2 -D 1 | 2 2ε . ( 2
L(D 2 ) -L(D 1 ) : (D 1 -D 2 ) = 1 ε D 2 -D 2 -D 1 + D 1 : (D 1 -D 2 ) = - 1 ε |D 1 -D 2 | 2 + 1 ε ( D 1 -D 2 ) : (D 1 -D 2 )
≤0.

(2.1.21)

Reversing the role of D 1 and D 2 in (2.1.20) and using the previous inequality yields

F ε (D 1 ) -F ε (D 2 ) -L(D 1 ) : (D 1 -D 2 ) ≤ |D 2 -D 1 | 2 2ε .
Hence

F ε (D 2 ) -F ε (D 1 ) -L(D 1 ) : (D 2 -D 1 ) ≤ |D 2 -D 1 | 2 2ε . (2.1.22)
We deduce that F ε is differentiable and

F ε (D) = L(D) = D -prox ε F (D) ε . Using (2.1.18) we write (D 1 -D 2 -D 1 + D 2 ) : ( D 1 -D 2 -D 1 + D 2 ) ≥ (D 1 -D 2 -D 1 + D 2 ) : (D 2 -D 1 ), thus |D 1 -D 2 -D 1 + D 2 | 2 ≤ (D 1 -D 2 -D 1 + D 2 ) : (D 2 -D 2 ), which implies that |D 1 -D 2 -D 1 + D 2 | ≤ |D 1 -D 2 |. ( 2
F (D) + |D -D| 2 /2ε ≥|D -D| 2 /2ε + A + σ 0 : D -|σ 0 ||D -D| ≥|D -D| 2 /4ε + A + σ 0 : D -ε|σ 0 | 2 ≥r 2 /4ε + A + σ 0 : D -ε|σ 0 | 2 .
Thus for ε small enough this will be larger than λ, and it follows that F ε (D) ≥ λ, which finishes the proof that F ε (D) tends to F (D) as ε tends to 0. Remark: Since F is convex, proper and lower semi-continuous, it follows that F is continuous on the interior of its domain [START_REF] Tyrell | Convex analysis[END_REF].

The regularization approach

One of the intriguing obstacles in the modelling of viscoplastic media is the presence of unknown interfaces separating the yielded and the un-yielded regions, that are difficult to track. In (2.0.1) they correspond to locations where the stress σ ∈ ∂F (Du) switches from a place where F is differentiable to one that is not. To avoid the associated numerical difficulties, the regularization method replaces the nonlinearity F by a smooth one, with a small approximation parameter ε. Then the material behaves as a fluid in the entire computational domain, and we have a differentiable, convex, proper, l.s.c .function F ε . Note that here F ε can be any function that approximates F . It can be the Moreau envelope of F defined by (2.1.12), or any other approximation. The approximate system is then

αu -div F ε (Du) = f in Ω, u = 0 on Γ.
This system can be solved by an iterative method with linearization, which is possible since F ε is smooth. The classical theory of monotone operators is developed in [START_REF] Brezis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF][START_REF] Trémolières | Numerical analysis of variational inequalities[END_REF][START_REF] Duvaut | Les inéquations en mécanique et en physique. dunod, paris[END_REF]. The solution to the former problem can be formulated via a minimization problem as

inf u∈H Ω α |u| 2 2 + F (Du) -f • u dx,
where H is an appropriate Hilbert space of functions u(x) such that u → F (Du) is convex, proper and l.s.c. There is here a difficulty to define such space for a general nonlinearity F , in particular concerning the lower semi-continuity. This will be studied in Chapter 3 for quite general nonlinearities F . In the regularization approach we have the approximate differentiable nonlinear optimization problem

inf uε∈H Ω α |u ε | 2 2 + F ε (Du ε ) -f • u ε dx .
In fact the regularization approach has until now only been considerd for particular nonlinearities F , that is for the Herschel

or inf u∈H Ω α |u| 2 2 + η |Du| 2 2 + σ 0 |Du| -f • u dx .
In the approximate problem we look for

u ε ∈ H realizing inf uε∈H Ω α |u ε | 2 2 + η |Du ε | 2 2 + σ 0 ( |Du ε | 2 + ε 2 -ε) -f • u ε dx,
where ε > 0. It corresponds to solve the well-posed nonlinear parabolic problem

αu ε -η div Du ε -σ 0 div Du ε |Du ε | 2 + ε 2 = f in Ω. (2.2.2)
It can be solved by an iterative algorithm,

αu k+1 ε -η div Du k+1 ε -σ 0 div Du k+1 ε |Du k ε | 2 + ε 2 = f in Ω, (2.2.3)
each iteration giving rise to the resolution of a linear problem. An overall space discretization has to be performed, usually by finite elements. The regularization method has certain advantages, namely it can use classical numerical schemes and it can be implemented easily in many existing codes. One of the key points in this method is finding the optimal value for ε. Indeed the smaller is ε the smaller is the error done in the law (one can prove that u -u ε L 2 √ ε), but larger is the number of iterations needed. Thus for a given tolerance one has to choose the corresponding ε, as well as the corresponding space size ∆x, as usual. Therefore some optimal relation between ε and ∆x has to be found. In [START_REF] Bouchut | Convergence of conforming approximations for inviscid incompressible Bingham fluid flows and related problems[END_REF] the authors propose an optimal relation of the form ε ∼ (∆x) 2 for P 1 finite elements (or ε ∼ (∆x) 4 in some cases when the solution less regular). A good formula for the choice of ε is

ε ∼ 10 -2 ∆x 2 |Du| 2 α σ 0 , (2.2.4)
in the case when the viscosity is small. Most of the time the regularization runs quite fast. There are many studies using this classical approach, as [START_REF] Bercovier | A finite-element method for incompressible non-Newtonian flows[END_REF][START_REF] Chatzimina | Cessation of Couette and Poiseuille flows of a Bingham plastic and finite stopping times[END_REF][START_REF] Mitsoulis | Flow of Bingham plastics in a lid-driven square cavity[END_REF]. Some studies [START_REF] Lusso | Explicit solutions to a free interface model for the static/flowing transition in thin granular flows[END_REF][START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a µ(I)-rheology[END_REF] depicted certain drawbacks of this method: the yield/plug zones (where Du = 0) can be not well resolved, or the asymptotic behaviour as t → ∞ can be wrong. This is the motivation for developing alternative approaches.

Augmented Lagrangian method

The second approach, especially designed to overcome the difficulty of non-differentiability of the constitutive law, strongly uses the theory of variational inequalities that one can find in Duvaut and Lions [START_REF] Duvaut | Les inéquations en mécanique et en physique. dunod, paris[END_REF]. The viscoplastic problem can be formulated as a saddle point problem with Lagrange multipliers. Then the well-known Augmented Lagrangian method for solving constrained optimization problems, first introduced by Fortin and Glowinski [START_REF] Fortin | Augmented Lagrangian methods in quadratic programming[END_REF], can be applied. It is also is called ADMM algorithm (Alternate directions method of multipliers), see [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF].

Recall first that the viscoplastic problem (2.0.1) can be formulated via the convex-concave Lagrangian

L(σ, u) = Ω α |u| 2 2 + σ : Du -F * (σ) -f • u dx, (2.3.1)
as the two equivalent problems

inf u sup σ L(σ, u) = sup σ inf u L(σ, u). (2.3.2)
The principle of the method is to introduce an additional variable γ that is an approximation of Du, and to set the Lagrangian multiplier corresponding to the constrain γ -Du = 0. Recalling that F * (σ) = sup γ (σ : γ -F (γ)), the problem (2.3.2) can be rewritten as

sup σ inf u,γ L(σ, u, γ), (2.3.3)
with

L(σ, u, γ) = Ω α |u| 2 2 + σ : (Du -γ) + F (γ) -f • u dx . (2.3.4)
The augmented formulation is to rewrite the problem as

sup σ inf u,γ L aug (σ, u, γ), (2.3.5)
with the augmented Lagrangian

L aug (σ, u, γ) = L(σ, u, γ) + r 2 Du -γ 2 L 2 , (2.3.6)
where r > 0 is a free parameter. The augmented term does not change the minimality condition, nor the value of the minimum, but enforces the stability of the numerical methods. The algorithm proposed by Fortin and Glowinski [START_REF] Fortin | Augmented Lagrangian methods in quadratic programming[END_REF] contains the following steps:

• Initialize with given σ 0 , γ 0 .

• For k ≥ 0, assuming that σ k , γ k are known, calculate u k by u k := argmin u L aug (σ k , u, γ k ), i.e.

u k solves αu k -div σ k -div r(Du k -γ k ) = f. (2.3.7) • Calculate γ k+1 by γ k+1 := argmin γ L aug (σ k , u k , γ), i.e. ∂F (γ k+1 ) σ k -r(γ k+1 -Du k ). (2.3.8)
• Calculate σ k+1 by going into the direction of the gradient of L aug with respect to σ, i.e.

σ k+1 = σ k + r(Du k -γ k+1 ). (2.3.9)
The first step consists in solving an elliptic problem, which can be done by the finite elements method (FEM). The two latter steps are nonlinear but local, they can be calculated pointwise efficiently. With the definition (2.4.1) of P r and the Moreau identity (2.1.16), the update of γ and σ can be written equivalently as

γ k+1 = (Id +∂F/r) -1 σ k + rDu k r = σ k + rDu k -σ k+1 r , σ k+1 = P r (σ k + rDu k ). (2.3.10)
The Augmented Lagrangian method was introduced in the 70s but is still used, one can find it in the works of Saramito et al. [START_REF] Bresch | Augmented Lagrangian method and compressible visco-plastic flows: applications to shallow dense avalanches[END_REF]. However, the rate of convergence is still limited in comparison with the regularization method when accurate determination of solid zones is not required. Another delicate issue in this method is the determination of an appropriate value of r in order to have a good convergence rate. An improvement with this respect is the Bermudez-Moreno scheme. Besides, the preconditioned version of the ADMM method is equivalent to the primal-dual algorithm proposed in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], which is studied in Chapter 4.

Projection approach

We would like here to introduce a class of methods that have a formulation which is close to the Augmented Lagrangian method, retaining the essential projection step (2.3.10). In particular, the algorithms used in imaging such as [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] fall into this class. The most challenging task when solving the viscoplastic model (2.0.1) is to solve the constraint σ ∈ ∂F (Du), since a single value of Du leads to several possible values of σ. In the projection approach, the main idea is to formulate the constraint via a fixed point on σ involving a projection. This kind of formulation is used in [START_REF] Chupin | A bi-projection method for Bingham type flows[END_REF][START_REF] Chalayer | A bi-projection method for incompressible Bingham flows with variable density, viscosity, and yield stress[END_REF] in the context of incompressible Navier-Stokes equations with a viscous Bingham viscoplastic law.

In the incompressible setting, this approach is well suited in conjunction with the projection on free divergence vector fields, that is described in [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] and that is useful for the simulation of turbulent flows. The approach of [START_REF] Chupin | A bi-projection method for Bingham type flows[END_REF] relies strongly on the presence of viscosity and can be interpreted as a time regularization. In this subsection, a projection formulation is introduced for a quite general nonlinearity F , and viscosity is not necessary.

As in Section 2.1 we consider a convex, proper and l.s.c. function F defined on a Hilbert space H, that will be taken later on as either the space M s N ×N (R) of symmetric square matrices of size N or eventually a space of functions. Then the proximal operator of F * , as defined in Proposition 2.1.4, is well defined. We consider thus for r > 0 the operator P r defined as

P r (σ) = (Id +r∂F * ) -1 (σ), σ ∈ H. ( 2
P r (σ + rD) = σ ⇐⇒ σ + rD ∈ (Id +r∂F * )(σ) ⇐⇒ D ∈ ∂F * (σ) ⇐⇒ σ ∈ ∂F (D).
(d) It follows from (c) by taking D = 0.

The previous lemma enables to formulate the constraint σ ∈ ∂F (D) of (2.0.1) as σ being a fixed point of the map σ → P r (σ + rD), which is a 1-Lipschitz mapping. Then one can think of an iteration procedure in order to get a solution σ. However the mapping is not a contraction (which is related to the fact that there could be several solutions), and thus the convergence is slow, if it holds. This is studied in Chapter 4 when the iteration is coupled with the momentum equation. Indeed according to Lemma 2.4.1, the viscoplastic model (2.0.1) can be rewritten as

αu -div σ = f, σ = P r (σ + rDu), (2.4.5)
where r is any positive constant. Note that the set ∂F (0) plays a particular role here since σ ∈ ∂F (0) corresponds to the stresses σ that are admissible in solid zones i.e. where Du = 0.

An issue is how to compute P r if F is known (but F * is not known explicitly). One can think of using the Moreau identity (2.1.16). When F has some regularity outside the origin we can use the following formulas. Proposition 2.4.2. Assume that F is a convex, proper, l.s.c. function on H, and that F is finite everywhere, and F is differentiable outside the origin.

(2.4.6)

Then for σ / ∈ ∂F (0), there exists a unique D = 0 such that F (D) + rD = σ. For r > 0 we then have

P r (σ) = σ if σ ∈ ∂F (0), F ((F + r Id) -1 (σ)) if σ / ∈ ∂F (0).
(2.4.7)

Proof. According to Moreau's identity (2.1.16), we have Hence we get

P r (σ) = (Id +r∂F * ) -1 (σ) = σ -r Id + ∂F r -1 σ r . ( 2
F ((F + r Id) -1 (σ)) + r Id + ∂F r -1 σ r = F (D) + rD = σ.
Therefore inserting it into (2.4.8), for σ / ∈ ∂F (0) we conclude that P r (σ) = F ((F + r Id) -1 (σ)).

Pure plastic models

We provide here particular properties for pure plastic models, which means that F is homogeneous of degree 1. The mechanical interpretation is that in this case the law does not include any viscous effect. These models are also called "rate independent" because when σ ∈ ∂F (D), σ depends only on the direction of D, but not on its magnitude.

Proposition 2.5.1. Suppose that F is a convex, l.s.c, proper function on H. Then F is homogeneous of degree 1 i.e. F (λD) = λF (D), for all λ > 0 and D ∈ H,

if and only if there exists A ⊂ H convex, closed and nonempty such that

F * (σ) = 0 if σ ∈ A, ∞ if σ / ∈ A.
(2.5.2)

Moreover, we have then

F (D) = sup σ∈A σ : D. (2.5.3)
Proof. Let us start from recalling the Fenchel conjugate function

F * (σ) = sup D (σ : D -F (D)). If F
is homogeneous, then for all σ ∈ H and λ > 0, we have

F * (λσ) = sup D (λσ : D -F (D)) = sup D λσ : D λ -F D λ = sup D σ : D - F (D) λ = 1 λ sup D (λσ : D -F (D)) = 1 λ F * (λσ).
Taking λ = 2 we deduce that F * (2σ) = 0 or ∞, for all σ. It means that F * (σ) = 0 or ∞ for all σ.

Hence F * takes the form

F * (σ) = 0 if σ ∈ A, ∞ if σ / ∈ A,
for some set A ⊂ H. We know that F * is convex, l.s.c, and proper. Since F * is proper it implies that A = ∅. Since F * is convex we get that A is necessarily convex. Lastly, to prove that A is closed, take an arbitrary sequence

σ n ∈ A converging to σ ∈ H. Because F * is lower semi-continuous, F * (σ) ≤ limF * (σ n ) = 0. Thus F * (σ) = 0 and σ ∈ A. It proves that A is closed.
Conversely, if we define F * by (2.5.2) for some A convex, closed and non-empty, as previously we have that F * is convex, lower semi-continuous and proper. Then

F (D) = F * * (D) = sup σ (D : σ -F * (σ)) = sup σ∈A (σ : D). For all λ > 0, D ∈ H, we have F (λD) = sup σ∈A (λσ : D) = λ sup σ∈A (σ : D) = λF (D).
Hence F is homogeneous of degree 1. Moreover the formula (2.5.3) is proved.

Lemma 2.5.2. Suppose that F is convex, l.s.c, proper, and homogeneous of degree 1. Then we have the following properties: (a)

F (0) = 0. (b) The set A in (2.5.2) is given by A = ∂F (0). (c) A is bounded if and only if dom F = H. Proof. (a) Taking D = 0, λ = 2 in (2.5.1), we deduce F (0) = 2F (0). Thus F (0) = 0 or F (0) = +∞. If F (0) < ∞ then we are done. Otherwise, since F is proper, there exists D = 0 such that F (D) < ∞. Since F is l.s.c, we have F (0) ≤ lim λ→0 + F (λD) = lim λ→0 + λF (D) = 0.
This implies that F (0) = 0. (b) We have the equivalences

σ ∈ A ⇐⇒ F * (σ) ≤ 0 ⇐⇒ sup D (σ : D -F (D)) ≤ 0 ⇐⇒ F (D) ≥ σ : D ∀D ⇐⇒F (D) -F (0) ≥ σ : (D -0) ∀D ⇐⇒ σ ∈ ∂F (0).
Hence A = ∂F (0). (c) Because of (2.5.3), the property that F (D) < ∞ for all D ∈ H means that A is weakly bounded, which is equivalent to A bounded by the Banach-Steinhaus theorem.

Proposition 2.5.3. Suppose that F * is defined as (2.5.2) for some A convex, closed and nonempty. Then we have:

(a) If σ ∈ int(A), then ∂F * (σ) = {0}. (b) If σ / ∈ A, then ∂F * (σ) = ∅. (c) If σ ∈ ∂A, then D ∈ ∂F * (σ) ⇐⇒ D : (τ -σ) ≤ 0 ∀τ ∈ A. (2.5.4)
When D satisfies (2.5.4), we say that D is an outer direction of A at σ ∈ ∂A.

(d) In the case of a finite dimensional space H, if A has a C 1 boundary and intA = ∅, then for

σ ∈ ∂A, ∂F * (σ) = R + n(σ),
where n(σ) is the exterior normal of ∂A at the point σ.

Proof. (a) Assume that σ ∈ int(A). Since F * ≥ 0 and F * (σ) = 0, it follows that 0 ∈ ∂F * (σ).

Conversely suppose that D ∈ ∂F * (σ). Then there exists ε > 0 such that σ + εD ∈ A. It follows that F * (σ + εD) = 0. Since D ∈ ∂F * (σ),

F * (σ + εD) ≥ F * (σ) + D : (σ + εD -σ), thus 0 ≥ ε|D| 2 and D = 0. Therefore ∂F * (σ) ⊂ {0}. (b) If σ / ∈ A, we have F * (σ) = ∞,
F * (τ ) ≥ F * (σ) + D : (τ -σ) ∀τ, 0 ≥ D : (τ -σ) ∀τ ∈ A.
This proves the claim. (d) The case of C 1 boundary is obvious and not detailed here.

Remark: Some references call a function of the form (2.5.2) the indicator of the closed convex set A, denoted by Id A . And the subdifferential ∂F * (σ) of an indicator of a closed convex set A, for σ ∈ ∂A, is called a "normal cone". Proposition 2.5.4. Suppose that F is convex, l.s.c, proper and homogeneous of degree 1.

(a) If D = 0, then ∂F (D) = A. (b) If D = 0, then ∂F (D) = {σ ∈ ∂A | D is an outer direction of A at σ}.
Proof. It follows from Proposition 2.5.3 by using the fact that D ∈ ∂F * (σ) ⇐⇒ σ ∈ ∂F (D).

Remark: In the case of C 1 boundary as stated in Proposition 2.5.3(d), when we have the relation σ ∈ ∂F (D) and D varies (D = 0), σ remains in ∂A. Then a small variation δD of D induces a small variation δσ of σ, that verifies δσ : D = 0. In such situation, mechanics call the law F "associated". Slightly different presentations of pure plastic models are possible, see [START_REF] Carstensen | Some properties of the dissipative model of strain-gradient plasticity[END_REF].

Theorem 2.5.5. Suppose that F is convex, l.s.c, proper and homogeneous of degree 1. Let A = ∂F (0). If we define P r (σ) = (Id +r∂F * ) -1 (σ) for r > 0, then P r (σ) is independent of r and is the orthogonal projection of σ on A. To find the projection P r (σ) when σ / ∈ A, we have to find σ b ∈ ∂A and n b an outer direction of A at σ b , such that

σ = σ b + n b .
(2.5.5)

Then P r (σ) = σ b .
Proof. We recall that for a nonempty closed convex set A we can define uniquely an orthogonal projection P A on A such that for any σ ∈ H, P A (σ) is characterized by

(σ -P A (σ)) : (σ -P A (σ)) ≤ 0 ∀σ ∈ A. (2.5.6)
By definition of the projection P r , there exists some D ∈ H such that P r (σ) + rD = σ with P r (σ) ∈ ∂F (D).

(2.5.7)

Since D ∈ ∂F * (P r (σ)), we have that ∂F * (P r (σ)) = ∅ and that by Proposition 2.5.3 we obtain P r (σ) ∈ A. Thus either P r (σ) ∈ int(A) or P r (σ) ∈ ∂A.

If P r (σ) ∈ int(A) then D = 0 and P r (σ) = σ = P A (σ). If P r (σ) ∈ ∂A, since D ∈ ∂F * (P r (σ)) we have according to Proposition 2.5.3(c) 0 ≥ D : (σ -P r (σ)) ∀σ ∈ A, 0 ≥ σ -P r (σ) r : (σ -P r (σ)) ∀σ ∈ A.
This means that P r (σ) is the orthogonal projection of σ on the convex set A. This proves that

P r = P A .
For the formula (2.5.5), consider σ ∈ A. From the definition one has P r (σ) = (Id +r∂F * ) -1 (σ), which means that

P r (σ) + r∂F * (P r (σ)) σ. (2.5.8)
It follows that ∂F * (P r (σ)) = ∅. By Proposition 2.5.3 and since σ ∈ A we get that P r (σ) ∈ ∂A. Then (2.5.8) means that σ -P r (σ) is an outer direction of A at P r (σ), which proves (2.5.5).

Examples

In this subsection we apply the results of the previous section to various pure plastic models, corresponding to different choices of F (D) defined on the space H = M s N ×N (R) of symmetric square matrices of size N .

Bingham model

In this model we take F (D) = |D|. This is the most simple choice. Using that A = ∂F (0), we obtain that σ ∈ A is characterized by

σ : D ≤ |D| ∀D.
It follows that A = {σ : |σ| ≤ 1}. We can represent A geometrically as a half circle in the variables p, |σ | with σ = -p Id +σ , Tr σ = 0, as shown on Figure 2.1a. In [START_REF] Chupin | A bi-projection method for Bingham type flows[END_REF], L. Chupin and T. Dubois used the associated projection P A on A. It is given simply by

P A (σ) =    σ if |σ| ≤ 1, σ |σ| if |σ| > 1.
(2.6.1)

Degenerate Bingham model

In this model we take (2.6.5)

F (D) = |D |,
A O p |σ | -σ -P r (σ) (a) Bingham model A O p |σ | -σ -P r (σ) -σ -P r (σ) (b) Degenerate Bingham model
     σ = rD + D |D | if |D | = 0, Tr(σ -rD) = 0, |σ -rD| ≤ 1 if |D | = 0.
(2.6.6)

If σ = rD + D |D | with |D | = 0, then Tr(σ) Id N + σ = r Tr(D) Id N + D + D |D | , (2.6.7)
which is equivalent to

     r Tr(D) = Tr(σ), rD + D |D | = σ . (2.6.8)
By the second equation |σ | > 1, A and σ are co-linear, and taking the absolute value we obtain

     Tr(D) = Tr(σ) r , |D | = |σ | -1 r , D = |σ | -1 r σ |σ | .
(2.6.9)

Then D = Tr(σ) rN Id + |σ | -1 r σ |σ | . ( 2

Therefore we obtain in any case

Id + ∂F r

-1 σ r = Tr(σ) rN Id + |σ | -1 r + σ |σ | . (2.6.11)
With (2.6.5) we deduce that

P A (σ) =      σ if |σ | ≤ 1, σ |σ | if |σ | > 1.

Cam-Clay model

We consider the well-known Cam-Clay viscoplastic model. Defining the set A or the function F is equivalent, and we choose here to first define A. Afterwards we shall deduce F , as (2.6.21). Writing σ = -p Id +σ with Tr(σ ) = 0, we define A in the variables σ and p as a half-ellipse,

A := σ | 0 ≤ p ≤ p 0 , |σ | 2 sin 2 δ ≤ p(p 0 -p) = σ | 0 ≤ p ≤ p 0 , |σ | 2 sin 2 δ + p - p 0 2 2 ≤ p 2 0 4 ,
(2.6.12) where p 0 > 0 and sin δ > 0 are given, see Figure 2.2. Then

F (D) = sup 0≤p≤p 0 sup |σ |≤sin δ √ p 0 p-p 2 (-p Tr(D) + σ : D ) = sup 0≤p≤p 0 -p Tr(D) + |D | sin δ p 0 p -p 2 .
(2.6.13)

We set g(p) = -p Tr(D) + |D | sin δ p 0 p -p 2 . Writing that the derivative of g vanishes at the maximum of the concave function g, we obtain

-Tr(D) + |D | sin δ p 0 -2p 2 p 0 p -p 2 = 0, |D | sin δ(p 0 -2p) = Tr(D)2 p 0 p -p 2 .
(2.6.14)

We remark that it implies that Tr(D) has the same sign as p 0 -2p, thus

Tr(D) > 0 ⇒ p < p 0 2 , Tr(D) < 0 ⇒ p > p 0 2 .
(2.6.15)

Then we take the square in (2.6.14), giving

|D | 2 sin 2 δ(p 2 0 -4p 0 p + 4p 2 ) = 4 Tr(D) 2 (p 0 p -p 2 ), (2.6.16
)

p 2 -pp 0 + |D | 2 sin 2 δ 4(|D | 2 sin 2 δ + Tr(D) 2 )
p 2 0 = 0, (2.6.17)

p = p 0 2 ± p 0 2 Tr(D) |D | 2 sin 2 δ + Tr(D) 2 .
(2.6.18)

Recalling from (2.6.15) that Tr(D) has the same sign as p 0 -2p, we get that the sign is minus, thus The relation between p and Y can also be inverted as

p = p 0 2 1 - Y √ 1 + Y 2 , with Y = Tr(D) |D | sin δ . (2.6.19) A -σ -P r (σ) p |σ | O p 0
Y = p 0 -2p 2 p 0 p -p 2 .
(2.6.20)

Hence (2.6.13) yields

F (D) = - p 0 2 1 - Y √ 1 + Y 2 Tr(D) + |D | sin δ p 0 2 √ 1 + Y 2 .
Using that Tr(D) = Y |D | sin δ, it can be recast as

F (D) = p 0 2 sin δ|D | √ 1 + Y 2 + Y .
Thus

F (D) = p 0 2 sin δ|D | Y + √ 1 + Y 2 , with Y = Tr(D) |D | sin δ .
(2.6.21)

We have to mention that according to Lemma 2.1.3(c), the point σ where the supremum (2.6.13) is attained gives ∂F (D).

For |D | = 0 one has F (D) = 0 if Tr(D) ≥ 0, F (D) = -p 0 Tr(D) if Tr(D) ≤ 0.

Finding the orthogonal projection on the half ellipse

Since the orthogonal projection on A is involved in Theorem 2.5.5, it is useful to compute it. For a stress σ b on the boundary of A, one has an outer normal n b , such that for a small variation δσ b one has n b : δσ b = 0. In order to find the orthogonal projection of σ = -p Id +σ on A (σ ∈ A), according to Theorem 2.5.5 one has to find σ b ∈ ∂A and λ ≥ 0 such that σ = σ b + λn b . Then and it follows that we can take .6.26) This formula could also be deduced from the computation of the supremum above, by Proposition 2.5.3(d), with the help of Lemma 2.1.3(b,c), and using the inversion formula (2.6.20). Thus in order to find the orthogonal projection of σ = -p Id +σ on A (σ ∈ A), one has to find σ b ∈ ∂A and λ ≥ 0 such that

P A (σ) = σ b . One has σ b = -p b Id +σ b sin δ p b (p 0 -p b ), ( 2
n b = σ b + sin δ N p 0 -2p b 2 p b (p 0 -p b ) Id . ( 2
σ ≡ -p Id +σ = -p b Id +σ b sin δ p b (p 0 -p b ) + λ σ b + sin δ N p 0 -2p b 2 p b (p 0 -p b ) Id .
(2.6.27)

Taking the trace and the deviator parts we get that σ and σ b are aligned (thus σ b = σ /|σ |) and the equations

-p = -p b + λ sin δ N p 0 -2p b 2 p b (p 0 -p b ) , (2.6.28
)

|σ | = λ + sin δ p b (p 0 -p b ). (2.6.29)
These two equations determine 0 ≤ p b ≤ p 0 and λ ≥ 0, that should be unique under the assumption that σ ∈ A i.e. either p < 0, or p > p 0 , or 0 ≤ p ≤ p 0 and |σ | > sin δ p(p 0 -p). Inserting the expression of λ from (2.6.29) into (2.6.28) yields

-p = -p b + sin δ N p 0 -2p b 2 p b (p 0 -p b ) |σ | -sin δ p b (p 0 -p b ) , (2.6.30)
which is a quartic equation in p b . This can be rewritten equivalently in terms of Y b related to p b by (2.6.19), as

sin δ N 2 p 0 |σ |Y b + 1 - sin 2 δ N Y b 1 + Y 2 b + 2p p 0 -1 = 0. ( 2
       p = 0 if sin ψ|D | -Tr(D) < 0, 0 ≤ p ≤ p 0 if sin ψ|D | -Tr(D) = 0, p = p 0 if sin ψ|D | -Tr(D) > 0.
(2.6.34)

In particular, having a pressure p < p 0 implies that Tr(D) ≥ 0. This justifies the name "Druckerdilatant": as long as the pressure does not take the maximal value, the material dilates. Moreover, as long as 0 < p < p 0 , we must have the dilatancy law Tr(D) = sin ψ|D |. (2.6.36)

If follows that

σ b = σ |σ | , |σ | = λ + p b sin ψ, -p = -p b + λ sin ψ N .
(2.6.37)

Plugging the value of λ into the equation on p we obtain

-p = -p b + (|σ | -p b sin ψ) sin ψ N , (2.6.38)
thus

p b = p + |σ | sin ψ N 1 + sin 2 ψ N , σ b = p + |σ | sin ψ N 1 + sin 2 ψ N -Id + sin ψ σ |σ | . (2.6.39)
This formula is valid if 0 ≤ p b ≤ p 0 , and λ ≥ 0, which gives the conditions (2.6.41)

0 ≤ p + |σ | sin ψ N ≤ 1 + sin 2 ψ N p 0 ,
It follows that

σ b = σ |σ | , |σ | = µ.
(2.6.42)

Finally σ b = -p 0 Id +σ . This is valid when p ≥ p 0 and |σ | ≤ p 0 sin ψ.

For the other values of σ, σ b is one of the corners. Thus finally, for σ ∈ A,

P A (σ) =                            0 if p + sin ψ N |σ | ≤ 0, p + sin ψ N |σ | 1 + sin 2 ψ N -Id + sin ψ σ |σ | if 0 ≤ p + sin ψ N |σ | ≤ 1 + sin 2 ψ N p 0 , p sin ψ ≤ |σ |, -p 0 Id +p 0 sin ψ σ |σ | if p + sin ψ N |σ | ≥ (1 + sin 2 ψ N )p 0 , |σ | ≥ p 0 sin ψ, -p 0 Id +σ if |σ | ≤ p 0 sin ψ, p ≥ p 0 .
(2.6.43)

A logarithmic model

In this model we define A as

A := σ | 0 ≤ p ≤ ep 0 and |σ | ≤ λp 1 -ln p p 0 , (2.6.44)
for some given p 0 > 0, λ > 0, and where e = exp(1), see Knowing that we have the constrain that 0 ≤ p ≤ ep 0 , we deduce that

• If p max ≤ ep 0 then the supremum is attained at p = p max .

• If p max > ep 0 then the supremum is attained at p = ep 0 .

We check that

p max ≤ ep 0 ⇐⇒ p 0 exp -Tr(D) λ|D | ≤ ep 0 ⇐⇒ -Tr(D) λ|D | ≤ 1 ⇐⇒ Tr(D) ≥ -λ|D |.

It follows that

• If Tr(D) ≥ -λ|D | then the supremum is attained at p = p max .

• If Tr(D) < -λ|D | then the supremum is attained at p = ep 0 .

Then for Tr(D) ≥ -λ|D |, taking p = p max gives

F (D) = λ p max |D | - Tr(D) λ -|D |p max ln p max p 0 = λ p 0 exp - Tr(D) λ|D | |D | - Tr(D) λ + |D |p 0 exp - Tr(D) λ|D | Tr(D) λ|D | = λp 0 |D | exp - Tr(D) λ|D | . Otherwise if Tr(D) < -λ|D |, taking p = ep 0 gives F (D) = λ ep 0 |D | - Tr(D) λ -|D |ep 0 = -ep 0 Tr(D).
(2.6.47)

Finally we obtain the general formula

F (D) =      λp 0 |D | exp - Tr(D) λ|D | if Tr(D) ≥ -λ|D |, -ep 0 Tr(D) if Tr(D) ≤ -λ|D |.
-p = -p b -γ λ N ln p b p 0 , σ = γ + λp b 1 -ln p b p 0 σ b . (2.6.55) If follows that σ b = σ |σ | , |σ | = γ + λp b 1 -ln p b p 0 , ( 2
.6.56)

p = p b + γ λ N ln p b p 0 . (2.6.57)
Inserting the expression of γ from (2.6.56) into (2.6.57) yields

p = p b + λ N ln p b p 0 |σ | -λp b 1 -ln p b p 0 . (2.6.58)
We observe graphically on Remark: In the three last models, the assumption that F is convex, lower semi-continuous, proper is not obvious from its expression. Nevertheless it follows by construction from Proposition 2.5.1 that depicts the relation between the function F and the set A, knowing that A is convex, closed and non-empty, which is true when it is defined by an inequality of the form |σ | ≤ f (p) with f concave. We remark that in the Bingham model, when σ ∈ A the pressure p = -Tr(σ)/N can take negative values. The other models presented above verify p ≥ 0 when σ ∈ A, which is more physical. One could also ask the stronger condition σ ≤ 0 in the sense of symmetric matrices. A sufficient condition for that is |σ | ≤ Existence of the solution for the steady viscoplastic model with nonlinear rheology Before going to the numerical approximation of (2.0.1), an existence theory is necessary. Moreover the analysis of Chapter 4 shows that we need a formulation in the local sense (not only in integral form). This is sometimes called the Euler-Lagrange equation, and it does not follow directly from the variational theory.

In the case of the (possibly viscous) incompressible Bingham model, i.e. F (D) = η 2 |D| 2 + σ 0 |D|, one can find in the literature some discussions on the existence of solutions for the inviscid case η = 0, as in [START_REF] Lions | Remarks on some nonlinear evolution problems arising in Bingham flows[END_REF] or [START_REF] Bouchut | Convergence of conforming approximations for inviscid incompressible Bingham fluid flows and related problems[END_REF]. In all the classical results of [START_REF] Duvaut | Inequalities in mechanics and physics[END_REF] the analysis relies on the viscous term, and strongly depends on the assumption η > 0. Meanwhile in [START_REF] Bouchut | Convergence of conforming approximations for inviscid incompressible Bingham fluid flows and related problems[END_REF] the results are established for the case η = 0. These results show that the two cases (viscous or inviscid) have to be treated by a slighly different analysis, relying on the theory of variational inequalities and monotone operators [START_REF] Brezis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF], but with a different choice of Hilbert space.

In this chapter we consider the viscoplastic model (2.0.1) with σ ∈ ∂F (Du) where F is convex, proper, l.s.c. on M s N ×N (R) the space of symmetric square matrices of size N , and satisfies some growth conditions. We recall that the problem can be formulated as

inf u∈H Ω α |u| 2 2 + F (Du) -f • u dx, (3.0.1)
where H is an appropriate Hilbert space of functions u(x) with values in R N . The key issue is to define the space H for a general nonlinearity F , in such a way that u → F (Du) is convex, proper and l.s.c. on H. Then one has to prove also the Euler-Lagrange equations, i.e. the local formulation (2.0.1).

With viscosity

We first consider the case with additional viscosity η > 0, i.e.

inf u∈H Ω α |u| 2 2 + η 2 |Du| 2 + F (Du) -f • u dx . (3.1.1)
Another way to formulate it would be to consider the problem (3.0.1) with a coercive nonlinearity F (D) + η 2 |D| 2 . We recall that Hypothesis 1 means that F is convex, l.s.c. and proper. This will always be assumed. Then from Section 2.1 we have that F is necessarily bounded below by an affine function. With the notation of Proposition 2.1.4 we have the following result.

Lemma 3.1.1. Suppose that F satisfies Hypothesis 1. Then | prox ε F (0)| is bounded independently of ε as long as 0 < ε ≤ ε 0 . Proof. Since F is proper, there exists D 0 such that F (D 0 ) < ∞. Then since prox ε F is 1-Lipschitz, one has | prox ε F (0)-prox ε F (D 0 )| ≤ |D 0 |. Thus proving that prox ε F (0) is bounded independently of ε is equivalent to proving that prox ε F (D 0 ) is bounded independently of ε. The definition of the Moreau envelope gives F ε (D 0 ) = inf D F (D) + |D -D 0 | 2 2ε ≤ F (D 0 ). (3.1.2)
We have that F is lower bounded by an affine function, thus F (D) ≥ A : D -B for some matrix A and some real number B. It follows that

F (D) + |D -D 0 | 2 2ε ≥ A : (D -D 0 ) + A : D 0 -B + |D -D 0 | 2 2ε ≥ - |D -D 0 | 2 4ε -ε|A| 2 + A : D 0 -B + |D -D 0 | 2 2ε ≥ |D -D 0 | 2 4ε + A : D 0 -B -ε 0 |A| 2 . (3.1.3) Therefore F ε (D 0 ) ≥ A : D 0 -B -ε 0 |A| 2 .
Moreover there is some M > 0 such that the above quantity is larger than

F (D 0 ) for |D -D 0 | √ ε ≥ M . It follows that F (D) + |D -D 0 | 2 2ε has to attain its infimum inside the ball |D -D 0 | ≤ M √ ε. Since ε ≤ ε 0 this set is bounded by |D 0 | + M √ ε 0 . Consequently, | prox ε F (D 0 )| ≤ |D 0 | + M √ ε 0 for all 0 < ε ≤ ε 0 .
In the following we shall consider several growth assumptions.

Hypothesis 2. For any D ∈ M s N ×N (R) and σ ∈ ∂F (D) one has

|σ| ≤ C(1 + |D|),
where C is a positive constant independent of D and σ.

Lemma 3.1.2. Under Hypothesis 1, 2, we have (a) the derivative of

F ε verifies |F ε (D)| ≤ C (1 + |D|), (3.1.4)
where C is a positive constant independent of D and ε ≤ ε 0 , (b) F is finite everywhere and for all D

|F ε (D)| ≤ C (1 + |D| 2 ), |F (D)| ≤ C (1 + |D| 2 ), (3.1.5)
where C is a constant independent of D and ε ≤ ε 0 .

Proof. (a) According to Proposition 2.1.4(e) we have

F ε (D) ∈ ∂F (prox ε F (D)), thus |F ε (D)| ≤ C(1 + | prox ε F (D)|).
Besides, Proposition 2.1.4(c) ensures that prox ε F is 1-Lipschitz, giving

|F ε (D)| ≤ C(1 + | prox ε F (0)| + |D|). ( 3
F * (σ) ≥ sup D σ : D -C (1 + |D| 2 ) . (3.1.11)
The supremum is attained at the value of D given by D = σ/(2C ). We deduce that

F * (σ) ≥ |σ| 2 2C -C 1 + |σ| 2 4C 2 , ( 3
| ≤ R, |D 2 | ≤ R. Let us then set D 3 = D 2 + R D 2 -D 1 |D 2 -D 1 | . (3.1.13)
We have then |D 3 | ≤ 2R and

D 2 = R |D 2 -D 1 | + R D 1 + |D 2 -D 1 | |D 2 -D 1 | + R D 3 , (3.1.14)
so that by convexity of F ,

F (D 2 ) ≤ R |D 2 -D 1 | + R F (D 1 ) + |D 2 -D 1 | |D 2 -D 1 | + R F (D 3 ). (3.1.15)
We deduce with (3.1.8) that

F (D 2 ) -F (D 1 ) ≤ |D 2 -D 1 | |D 2 -D 1 | + R (F (D 3 ) -F (D 1 )) ≤ |D 2 -D 1 | R 2M, (3.1.16) With M = C (1 + 4R 2 ). Therefore F (D 2 ) -F (D 1 ) ≤ |D 2 -D 1 |2C 1 R + 4R . (3.1.17)
The case D 2 = D 1 being trivial, we have now that this inequality holds for all D 1 , D 2 such that

|D 1 | ≤ R and |D 2 | ≤ R Now, given D 1 and D 2 , if max(|D 1 |, |D 2 |) ≤ 1/2 then we can take R = 1/2, so that 1/R + 4R = 4. Otherwise we must have R ≥ 1/2, which implies that 1/R + 4R ≤ 2 + 4R.
Thus in any case we get for all D 1 , D 2 We consider now an open bounded subset Ω of R N . We shall need a variant of the Sobolev space H 1 (Ω), adapted to the symmetric derivative Du = (∇u + (∇u) t )/2 where u takes values in R N . We define Proof. (a) ψ is convex. If u, v ∈ H 1s 0 (Ω) and θ ∈ (0, 1), then since F is convex

F (D 2 ) -F (D 1 ) ≤ |D 2 -D 1 |8C 1 + max(|D 1 |, |D 2 |) . ( 3
H 1s (Ω) = u ∈ L 2 (Ω), with values in R N , such that Du ∈ L 2 (Ω) , ( 3
ψ(θu + (1 -θ)v) = Ω F (θDu + (1 -θ)Dv) ≤ Ω θF (Du) + (1 -θ)F (Dv) = θψ(u) + (1 -θ)ψ(v).
F (Du n ) - ϕ(Du n ) ≥ 0, one has Ω lim(F (Du n ) -ϕ(Du n )) ≤ lim Ω (F (Du n ) -ϕ(Du n )).
Since ϕ(Du n ) → ϕ(Du) in L 2 and a.e., we deduce that Reciprocally, assume that (3.1.25) hold. Take v = u + θ(v 0 -u) where θ ∈ [0, 1] and v 0 ∈ H 1s 0 (Ω). Then Dv = Du + θ(Dv 0 -Du). Taking the square of both sides and the integral over Ω we get

Ω limF (Du n ) ≤ lim Ω F (Du n ). ( 3
(Ω) 1 ∆t Ω u • (v -u) + η 2 Ω |Dv| 2 - η 2 Ω |Du| 2 + ψ(v) ≥ ψ(u) + G(v -u), (3.1.25) ∀v ∈ H 1s 0 (Ω) 1 ∆t Ω u • (v -u) + η Ω Du : (Dv -Du) + ψ(v) ≥ ψ(u) + G(v -u). ( 3
Ω |Dv| 2 - Ω |Du| 2 = 2θ Ω Du : (Dv 0 -Du) + θ 2 Ω |Dv 0 -Du| 2 . (3.1.27)
Multiplying by η 2 , substituting in (3.1.25) and using that by convexity of ψ one has ψ u + θ(v 0 -u) ≤

(1 -θ)ψ(u) + θψ(v 0 ) (and observing that (3.1.25) implies that ψ(u) < ∞), we obtain

θ ∆t Ω u • (v 0 -u) + θη Ω Du : (Dv 0 -Du) + η 2 θ 2 Ω |Dv 0 -Du| 2 + θψ(v 0 ) -θψ(u) ≥ θG(v 0 -u).
Dividing by θ > 0 and letting θ → 0, we recover (3.1.26).

Remark: (a) Since

Ω Du : (Dv -Du) = 1 2 Ω |Dv| 2 - Ω |Du| 2 - Ω |Dv -Du| 2
, the variational inequality (3.1.26) can be written as

∀v ∈ H 1s 0 (Ω) 1 ∆t Ω u•(v -u)+ η 2 Ω |Dv| 2 - η 2 Ω |Du| 2 - η 2 Ω |Dv -Du| 2 +ψ(v) ≥ ψ(u)+G(v -u).
This inequality looks stronger than (3.1.25). However, Lemma 3.1.6 says that there is indeed equivalence with (3.1.25).

(b) The variational formulation (3.1.26) is powerful and plays a crucial role to establish estimates, see [START_REF] Duvaut | Inequalities in mechanics and physics[END_REF] and references therein. However it can be difficult to establish it when numerical approximations are involved. Indeed when using test functions in finite elements spaces, a large number of inequalities arise. Then it is useful to use the Euler-Lagrange local formulation, see Theorem 3.1.9 below.

Proposition 3.1.7. Assume that F satisfies Hypothesis 1, 2. Then there is a unique solution u ∈ H 1s 0 (Ω) to ∀v ∈ H 1s 0 (Ω)

1 ∆t Ω u • (v -u) + Ω F (Dv) - Ω F (Du) + η Ω Du : (Dv -Du) ≥ G(v -u), (3.1.28)
where G is a bounded linear form on H 1s 0 (Ω), and η, ∆t are positive constants (the correspondence with (3.1.1) is α = 1/∆t).

Proof. With ψ coming from Definition 3.1.4, the set of inequalities (3.1.28) identifies to (3.1.26). Thus according to Lemma 3.1.6 it is equivalent to (3.1.25). Using the argument of the proof of Lemma 3.1.6 applied to the first term u • (v -u) we obtain that it is also equivalent to ∀v ∈ H 1s 0 (Ω)

1 ∆t Ω 1 2 (|v| 2 -|u| 2 ) + η 2 Ω |Dv| 2 - η 2 Ω |Du| 2 + ψ(v) -ψ(u) ≥ G(v -u),
which means that u is a minimum over H 1s 0 (Ω) of the functional

J(v) = 1 ∆t Ω 1 2 |v| 2 + η 2 Ω |Dv| 2 + ψ(v) -G(v). (3.1.29)
This functional is the sum of the square of the H 1s norm and a l.s.c. function. Thus by Proposition 2.1.4(a), there is a unique minimum in H 1s 0 (Ω).

The Moreau envelope F ε of F enables to build uniformly bounded approximations to our variational problem (3.1.28). Lemma 3.1.8. Assume that F satisfies Hypothesis 1, 2. Then for any ε > 0 there is a unique solution u ε ∈ H 1s 0 (Ω) to ∀v ∈ H 1s 0 (Ω)

1 ∆t Ω u ε • (v -u ε ) + Ω F ε (Dv) - Ω F ε (Du ε ) + η Ω Du ε : (Dv -Du ε ) ≥ G(v -u ε ), ( 3
Ω |u ε | 2 + η Ω |Du ε | 2 + Ω F ε (Du ε ) ≤ G(u ε ) + Ω F ε (0). ( 3
(v) = f n + u n ∆t , v + η Du n , Dv for some f n ∈ L 2 (Ω), u n ∈ H 1s (Ω), is equivalent to finding (u, σ) ∈ H 1s 0 (Ω) × L 2 (Ω) such that    u -u n ∆t -div σ -η div Du = -η div Du n + f n in Ω,
σ ∈ ∂F (Du) a.e. in Ω.

(3.1.32)

Moreover one has F * (σ) ∈ L 1 (Ω).

Proof. Assume that (3.1.32) holds, the first equation being understood in the sense of distributions.

Then by density we can apply it to test functions in H 1s 0 (Ω). Taking v -u as test function, where v ∈ H 1s 0 (Ω), we obtain denoting by ., . the L 2 (Ω) scalar product, Reciprocally, assume that (3.1.28) holds. Following Proposition 2.1.4, we can consider the sequence of Moreau envelopes F ε of F, which are convex and continuously differentiable and satisfy F ε ↑ F pointwise as ε → 0. We define u ε as the unique solution to the problem in which we replace F by F ε , i.e. u ε solves (3.1.30). At fixed ε we take then v = u ε + sw in (3.1.30), where w ∈ H 1s 0 (Ω) and s is a real number. Since F ε is differentiable, we have

1 ∆t u -u n , v -u + σ, Dv -Du + η Du, Dv -Du = η Du n , Dv -Du + f n , v -u . ( 3
F ε (Du ε + sDw) -F ε (Du ε ) s -→ F ε (Du ε ) : Dw a.e. as s -→ 0.
Moreover, since F ε is 1/ε-Lipschitz, we also have for some 0 ≤ θ s (x) ≤ 1

F ε (Du ε + sDw) -F ε (Du ε ) s = F ε (Du ε + θ s sDw) : Dw, F ε (Du ε + sDw) -F ε (Du ε ) s ≤ |F ε (0)| + 1 ε |Du ε + θ s sDw| |Dw| ≤ C + |Du ε | + |s||Dw| ε |Dw|.
Therefore by Lebesgue's dominated convergence theorem, [START_REF] Goldstein | Adaptive primaldual hybrid gradient methods for saddle-point problems[END_REF] where s > 0 and w ∈ H 1s 0 (Ω) and dividing by s, one gets

F ε (Du ε + sDw) -F ε (Du ε ) s -→F ε (Du ε ) : Dw in L 1 (Ω) as s → 0. Taking v = u ε + sw in (3.1.
1 ∆t u ε , w + Ω F ε (Du ε + sDw) -F ε (Du ε ) s + η Ω Du ε : Dw ≥ G(w), thus letting s → 0 1 ∆t u ε , w + Ω F ε (Du ε ) : Dw + η Ω Du ε : Dw ≥ G(w). (3.1.36)
Similarly by taking s < 0 we get

1 ∆t u ε , w + Ω F ε (Du ε ) : Dw + η Ω Du ε : Dw ≤ G(w), (3.1.37)
and it follows that

1 ∆t u ε , w + Ω F ε (Du ε ) : Dw + η Ω Du ε : Dw = G(w) ∀w ∈ H 1s 0 (Ω). (3.1.38) Thus (u ε , σ ε ) is a solution to u ε -u n ∆t -div σ ε -η div Du ε = -η div Du n + f n , (3.1.39) σ ε = F ε (Du ε ). (3.1.40)
In order to get (3.1.32) it remains to pass to the limit as ε → 0. According to Lemma 3.1.8, u ε and Du ε are bounded in L 2 (Ω) uniformly in ε ≤ ε 0 . Therefore with Lemma 3.1.2(a), σ ε is bounded in L 2 . Then according to Lemma 2.1.3(d) we have

Du ε : σ ε = F ε (Du ε ) + (F ε ) * (σ ε ). ( 3
(F ε ) * (σ ε ) is bounded in L 1 (Ω) independently of ε.
Hence after extraction of a subsequence, there exists u ∈ H 1s 0 (Ω) and σ ∈ L 2 (Ω) such that

u ε u, Du ε Du, σ ε σ, (3.1.43)
where denotes the weak convergence in L 2 , and additionally u ε → u locally strongly. Letting ε → 0 in (3.1.39), we obtain

u -u n ∆t -div σ -η div Du = -η div Du n + f n . (3.1.44)
Now we would like to prove that σ ∈ ∂F (Du) a.e. Since σ ε = F ε (Du ε ) a.e., we have

∀W F ε (W ) ≥ F ε (Du ε ) + σ ε : (W -Du ε ) a.e.
For ε ≤ ε 0 , we have

F ε (Du ε ) ≥ F ε 0 (Du ε ), thus ∀W F ε (W ) ≥ F ε 0 (Du ε ) + σ ε : (W -Du ε ) a.e.
Multiplying by ϕ ∈ C ∞ c (Ω), ϕ ≥ 0, and taking the integral over Ω we obtain

Ω ϕF ε (W ) ≥ Ω ϕF ε 0 (Du ε ) + Ω ϕσ ε : W - Ω ϕσ ε : Du ε . (3.1.45)
But multiplying (3.1.39) by ϕu ε , we get

Ω u ε ∆t • u ε ϕ + Ω σ ε : D(ϕu ε ) + η Ω Du ε : D(ϕu ε ) = G(ϕu ε ), Ω ϕ |u ε | 2 ∆t + Ω ϕσ ε : Du ε + Ω σ ε : (u ε ⊗∇ϕ)+η Ω ϕ|Du ε | 2 +η Ω Du ε : (u ε ⊗∇ϕ) = G(ϕu ε ). (3.1.46)
Similarly, taking ϕu as test function in (3.1.44), we get

Ω ϕ |u| 2 ∆t + Ω ϕσ : Du + Ω σ : (u ⊗ ∇ϕ) + η Ω ϕ|Du| 2 + η Ω Du : (u ⊗ ∇ϕ) = G(ϕu). (3.1.47)
Since u ε → u locally strongly, we can pass to the limit in the terms involving ∇ϕ. Thus taking the limit ε → 0 in (3.1.46) and comparing to (3.1.47), we obtain

Ω ϕσ ε : Du ε - Ω ϕσ : Du + η Ω ϕ|Du ε | 2 -η Ω ϕ|Du| 2 -→ 0.
Since Du ε Du, we have This is true for all ϕ ∈ C ∞ c (Ω) such that ϕ ≥ 0, thus

Ω ϕ|Du| 2 ≤ lim Ω ϕ|Du ε | 2 .
F (W ) ≥ F (Du) + σ : W -σ : Du a.e.
This is true for all W . But since F is finite everywhere and convex on a finite dimensional space, it is continuous, and consequently for a.e. x ∈ Ω the inequality holds for all W . We conclude that σ ∈ ∂F (Du) a.e. in Ω. Therefore (u, σ) satisfies (3.1.32). We have already proved that this implies that u is the solution to (3.1.28), which is unique according to Proposition 3.1.7.

Finally the last assertion F * (σ) ∈ L 1 (Ω) follows from Lemma 2.1.3(d) that ensures that σ : Du = F (Du) + F * (σ) a.e.

Remark: In the particular case of a Bingham fluid, i.e. F (D) = |D|, the proof of Euler-Lagrange equations was provided in [START_REF] Duvaut | Inequalities in mechanics and physics[END_REF][START_REF] Trémolières | Numerical analysis of variational inequalities[END_REF]. In this case, additionally to the regularization method that we have used here, the use of the Hahn-Banach theorem is possible. A stronger version of the Euler-Lagrange equations is also available for the Bingham problem, see [START_REF] Carstensen | A natural nonconforming FEM for the Bingham flow problem is quasi-optimal[END_REF].

Without viscosity

We consider now the case without viscosity (3.0.1), that corresponds to the unmodified problem (2.0.1). We shall need several growth assumptions on F .

Hypothesis 3. F is lower bounded, which means that there is a constant E such that for all

D ∈ M s N ×N (R) E + F (D) ≥ 0. (3.2.1)
According to Lemma 2.1.3(f ), this assumption is equivalent to F * (0) < ∞.

Hypothesis 4. F is finite everywhere and there is a constant C > 0 such that for all D

F (D) ≤ C(1 + |D| 2 ). (3.2.2)
Hypothesis 5. F is superlinear, which means that

F (D) |D| -→ ∞ as |D| -→ ∞. (3.2.3)
We consider again the Moreau envelope F ε of F , as defined in Proposition 2.1.4. We have the following result. According to (2.1.12) one has

F ε (D) = inf D H ε (D, D), (3.2.7)
with

H ε (D, D) = F (D) + |D -D| 2 2ε . (3.2.8)
Let λ > 0 be given. Then according to (3.2.6) there exists R > 0 such that for all |D| ≥ R, one has F (D) ≥ λ|D|. Assuming that ε ≤ ε 0 for some fixed ε 0 > 0, one can assume that R > λε 0 . Then let us take D such that |D| ≥ R. We have for all D satisfying |D| ≥ R,

H ε (D, D) ≥ λ|D| + |D -D| 2 2ε . (3.2.9)
The right-hand side is larger than its minimum over all D, which is attained at D = (|D|λε) + D/|D|. Since |D| ≥ R > λε, plugging this value in the right hand side of (3.2.9) yields

∀|D| ≥ R, H ε (D, D) ≥ λ(|D| -λε) + λ 2 ε 2 . (3.2.10)
For |D| ≤ R, one has

H ε (D, D) ≥ -E + (|D| -R) 2 2ε . (3.2.11)
We deduce that inf We consider now an open bounded subset Ω of R N , and functions u ∈ L 2 (Ω) that take values in R N . We use as previously the symmetric derivative Du = (∇u + (∇u) t )/2 which is in general a distribution. Definition 3.2.2. Suppose that F satisfies Hypothesis 1, 3, 4, and define for Ω an open bounded subset of R N ,

D H ε (D, D) ≥ min λ|D| -λ 2 ε 0 , -E + (|D| -R) 2 /2ε 0 . ( 3
K F (Ω) = u ∈ L 2 (Ω) such that Du ∈ L 1 (Ω) and F (Du) ∈ L 1 (Ω) , (3.2

.13)

and

K F 0 (Ω) ⊂ K F (Ω) by u ∈ K F 0 (Ω) if and only if Du = Du in R N , (3.2.14)
where for a function w defined in Ω, w denotes the extension by 0 of w, i.e. the function defined over the whole space R N as w = w in Ω, w = 0 outside Ω.

Notice that since by Hypothesis 3 F is lower bounded by a constant -E, for any u ∈ L 2 (Ω) such that Du ∈ L 1 (Ω) one can consider Ω F (Du) ∈ R. Then this integral is finite if and only if Moreover we have then

F (Du) ∈ L 1 (Ω). Note also that C ∞ c (Ω) ⊂ K F 0 (Ω) ⊂ K F (Ω).
Ω F (Du) ≤ lim Ω F (Du n ). (3.2.16
)

Additionally, if u n ∈ K F 0 (Ω) then u ∈ K F 0 (Ω). Proof. (a) If u, v ∈ K F (Ω) and 0 ≤ θ ≤ 1, we have that w = (1 -θ)u + θv ∈ L 2 , Dw = (1 -θ)Du + θDv ∈ L 1
, and by convexity of

F Ω F (Dw) ≤ Ω (1 -θ)F (Du) + θF (Dv) = (1 -θ) Ω F (Du) + θ Ω F (Dv) < ∞. (3.2.17) Hence w ∈ K F (Ω). This proves that K F (Ω) is convex. Next if u, v ∈ K F 0 (Ω)
, by linearity of (3.2.14) one has also w ∈ K F 0 (Ω), proving the convexity of

K F 0 (Ω). Finally 0 ∈ C ∞ c (Ω) ⊂ K F 0 (Ω) ⊂ K F (Ω). (b) Since F is superlinear (Hypothesis 5), there is some R 0 > 0 such that for any R ≥ R 0 , for all |D| ≥ R one has |D| ≤ K R F (D), where K R > 0 verifies K R → 0 as R → ∞. Indeed one has just to take 1 K R = inf |D|≥R F (D)
|D| . We deduce that Du n is bounded in L 1 . Indeed,

Ω |Du n | = |Dun|≥R |Du n | + |Dun|<R |Du n | ≤ K R |Dun|≥R F (Du n ) + R|Ω| = K R Ω F (Du n ) -K R |Dun|<R F (Du n ) + R|Ω| ≤ K R Ω F (Du n ) + K R E|{|Du n | < R}| + R|Ω|, (3.2.18)
and it is enough to take a single R to get the boundedness in L 1 . Then the same estimate (3.2.18) gives that

|Dun|≥R |Du n | → 0, as R → ∞, uniformly in n, (3.2.19)
which gives that Du n is uniformly equi-integrable since Ω is bounded. Therefore we can extract a subsequence u n such that u n u in weak L 2 , for some u ∈ L 2 , and such that Du n w weakly in L 1 , where w ∈ L 1 . This implies that Du n , ϕ → w, ϕ for all ϕ ∈ C ∞ c (Ω). But since u n u in weak L 2 we have Du n , ϕ = -u n , div ϕ → -u, div ϕ . It follows that Du = w ∈ L 1 and Du n

Du weakly in L 1 . One can check that the map which to w(x) ∈ L 1 (Ω) associates F (w) is convex, l.s.c. on L 1 strong. It follows that it is also l.s.c. on L 1 -weak. Therefore, since Du n Du in L 1 -weak, we deduce that

Ω F (Du) ≤ lim Ω F (Du n ) < ∞.
This gives that u ∈ K F (Ω) and concludes the proof in the case of K F (Ω).

Finally, for the case u n ∈ K F 0 (Ω), the convergences u n u in weak L 2 (Ω) and Du n Du in weak L 1 (Ω) imply that u n u in weak L 2 (R N ) and Du n Du in weak L 1 (R N ). Thus we can pass to the limit in (3.2.14) and obtain u ∈ K F 0 (Ω).

Lemma 3.2.4. Assume that F satisfies Hypothesis 1, 3, 4, and that Ω is strictly star-shaped. Then for u ∈ K F 0 (Ω) there exists a sequence

u k ∈ C ∞ c (Ω) such that u k → u in L 2 (Ω) strong, Du k → Du in L 1 (Ω) strong, Ω F (Du k ) → Ω F (Du). (3.2.20)
Proof. The property for Ω to be star-shaped means that there exists a point x 0 ∈ Ω such that for any x ∈ Ω, one has [x 0 , x] ⊂ Ω. It is strictly star-shaped if additionally for any x ∈ ∂Ω,

[x 0 , x) ⊂ Ω. Consider u ∈ K F 0 (Ω). Then (3.2.14) ensures that u ∈ L 2 (R N ), Du ∈ L 1 (R N ). Moreover R N F (Du) -F (0) = Ω F (Du) -F (0) = Ω F (Du) -|Ω|F (0). (3.2.21) Consider then for λ > 0 v λ (x) = 1 λ u(x 0 + λ(x -x 0 )). (3.2.22) Then v λ ∈ L 2 (R N ), Dv λ ∈ L 1 (R N ) and as λ → 1 one has v λ → u in L 2 (R N ), Dv λ → Du in L 1 (R N ), and
R N F (Dv λ ) -F (0) = R N F Du(x 0 + λ(x -x 0 )) -F (0) = λ -N R N F Du -F (0) . (3.2.23)
Consider then a smoothing sequence ρ δ (x) on R N , and define w δ = ρ δ * v λ , that indeed depends on δ and λ. Then we have that w δ ∈ C ∞ c (R N ), and as δ → 0 at fixed λ, w δ → v λ in L 2 (R N ), Dw δ → Dv λ in L 1 (R N ). Moreover by Jensen's inequality

F (Dw δ (x)) = F Dv λ (x-y)ρ δ (y)dy ≤ F Dv λ (x-y) ρ δ (y)dy = (ρ δ * F Dv λ ) (x). (3.2.24) one has also (1 -θ)u + θv ∈ K F 0 since K F 0 is convex, and D((1 -θ)u + θv) = (1 -θ)Du + θDv. Then since F is convex ψ((1 -θ)u + θv) = Ω F ((1 -θ)Du + θDv) ≤ Ω (1 -θ)F (Du) + θF (Dv) = (1 -θ)ψ(u) + θψ(v).
(b) ψ is proper. Since F is finite everywhere one has F (0) < ∞, thus using that 0 ∈ K F 0 we deduce that ψ(0) = |Ω|F (0) < ∞. (c) ψ is l.s.c. Suppose that u n → u in L 2 (Ω). We have to prove that ψ(u) ≤ lim ψ(u n ). One can assume that lim ψ(u n ) < ∞, and extracting a subsequence if necessary, one can assume that ψ(u n ) tends to some limit, that is necessarily finite since ψ is lower bounded. Then for n large enough ψ(u n ) is finite, thus u n ∈ K F 0 (Ω). Applying Lemma 3.2.3(b) we deduce that u ∈ K F 0 (Ω) and (3.2.16) yields that ψ(u) ≤ lim ψ(u n ). Proposition 3.2.7. Assume Hypothesis 1, 3, 4, 5 and consider ∆t > 0, u n , f n ∈ L 2 (Ω). Then there is one and only one solution u ∈ L 2 (Ω) to the problem

∀v ∈ L 2 (Ω) 1 ∆t u, v -u + ψ(v) ≥ ψ(u) + G(v -u), (3.2

.34)

which moreover satisfies u = argmin

v∈L 2 (Ω) J(v), (3.2.35)
where J : L 2 (Ω) → R is defined by (Ω) satisfying div σ ∈ L 2 (Ω) and F * (σ) ∈ L 1 (Ω). Then we have

J(v) = 1 ∆t Ω |v| 2 2 + ψ(v) -G(v), (3.2
∀w ∈ K F 0 (Ω) Ω σ : Dw ≥ - Ω (div σ) • w. (3.2.37)
The meaning of the integral σ : Dw in (3.2.37) needs to be explained. Under the assumptions on F , both F and F * are lower bounded by a constant. Thus integrals like F (Dw) or F * (σ) are well defined, as finite or +∞. The properties F (Dw) ∈ L 1 (which holds since w ∈ K F 0 (Ω)) or F * (σ) ∈ L 1 (which is assumed here) mean exactly that these respective integrals are finite. Then according to the inequality (2.1.4) one has

σ : Dw ≤ F (Dw) + F * (σ) a.e. ( 3
σ : Dw k = - Ω (div σ) • w k . (3.2.39)
Then according to the Fenchel-Young inequality we have Remark: If F is even, i.e. F (-D) = F (D), then in the previous lemma we can apply the same result to -w instead of w, thus we obtain that there is indeed equality in (3.2.37). In the general case, we do not know if there could be a strict inequality in (3.2.37).

g k ≡ F (Dw k ) + F * (σ) -σ : Dw k ≥ 0 a.
Ω g k = Ω F (Dw k ) + Ω F * (σ) + Ω (div σ) • w k , (3.2.43) thus lim
Ω g k = Ω F * (σ) + Ω (div σ) • w + Ω F (
Theorem 3.2.9 (Euler-Lagrange equations). With the same assumptions as in Proposition 3.2.7 and if Ω is strictly star-shaped, the solution u ∈ L 2 (Ω) to the problem (3.2.34) is characterized by the existence of σ ∈ L 2 (Ω) such that Knowing that there is existence and uniqueness of a solution to (3.2.34), to conclude the theorem it remains only to prove that there is a solution (u, σ) ∈ L 2 × L 2 to (3.2.46), (3.2.47), (3.2.48), (3.2.49).

u ∈ K F 0 (Ω), F * (σ) ∈ L 1 (Ω), (3.2.46) u -u n ∆t -div σ = f n in Ω, ( 3
Ω u -u n ∆t • u + Ω F (Du) + Ω F * (σ) ≤ Ω f n • u. ( 3
∈ K F 0 (Ω) Ω u -u n ∆t • v + Ω σ : Dv ≥ Ω f n • v. ( 3
1 ∆t u -u n , v -u + σ, Dv -Du ≥ f n , v -u . ( 3
According to Lemma 3.1.3(b), Hypothesis 2 is satisfied. Therefore for η > 0 we can apply Proposition 3.1.7 and Theorem 3.1.9. It gives

u η ∈ H 1s 0 (Ω), σ η ∈ L 2 (Ω) satisfying u η -u n ∆t -div σ η -η div Du η = f n in Ω, ( 3
(Ω), 1 ∆t Ω u η • w + Ω σ η : Dw + η Ω Du η : Dw = Ω f n + u n ∆t • w. (3.2.58)
Taking w = u η we get

1 ∆t Ω |u η | 2 + η Ω |Du η | 2 + Ω σ η : Du η = Ω f n + u n ∆t • u η . ( 3
Ω |Du η | 2 + Ω E+F (Du η ) + Ω F (0)+F * (σ η ) = Ω f n + u n ∆t •u η +|Ω|(E+F (0)). (3.2.61) Writing that Ω f n + u n ∆t • u η ≤ 1 2∆t Ω |u η | 2 + ∆t 2 Ω f n + u n ∆t 2 ,
we obtain the estimate

1 2∆t Ω |u η | 2 +η Ω |Du η | 2 + Ω E+F (Du η ) + Ω F (0)+F * (σ η ) ≤ ∆t 2 Ω f n + u n ∆t 2 +|Ω|(E+F (0)).
(3.2.62) Since F ≥ -E and F * ≥ -F (0) all the terms on the left-hand side are nonnegative, and we therefore get bounds on these quantities uniformly with respect to η. In particular we have Remark: The condition (3.2.48), i.e. σ ∈ ∂F (Du) a.e. could be removed in the statement of the previous theorem, since as the end of the proof shows it is a consequence of the other conditions. Notice that the solution u is unique, but σ associated to u may be nonunique.

u η 2 L 2 ≤ ∆t 2 f n + u n ∆t 2 L 2 +2∆t|Ω|(E +F (0)), Ω F (Du η ) ≤ ∆t 2 f n + u n ∆t 2 L 2 +|Ω|F (0). (3.2.63) We claim that u η ∈ K F 0 (Ω). Indeed since u η ∈ H 1s 0 (Ω), there exists a sequence u k ∈ C ∞ c (Ω) such that u k → u η in L 2 (Ω) and Du k → Du η in L 2 (Ω). Then u k ∈ K F 0 (Ω),
Remark: For applications to granular flows, it is important not to assume any special dependency of F (D) with respect to D, such as dependency only on |D|. In particular, the above proofs would simplify somehow for an even nonlinearity F , i.e. verifying F (-D) = F (D), but we cannot assume that, as the examples of Section 2.6 show. The assumptions we have made on F are only on its asymptotic growth. Only the assumption of superlinearity Hypothesis 5 is a bit restrictive, since it excludes a linear behavior at infinity. For a 1-homogeneous function such as described in Section 2.5 one should be able to get well-posedness results with the approach of [START_REF] Bouchut | Convergence of conforming approximations for inviscid incompressible Bingham fluid flows and related problems[END_REF]. In this case Du could be a measure instead of an L 1 function, which modifies significantly the proofs. However it would be desirable to consider a nonlinearity F that can be superlinear in certain directions, and asymptotically linear in other directions. At the time being we do not know how to treat this case.

The superlinearity assumption can indeed be related to the finiteness of F * , as the following lemma shows. Conversely suppose that F is not superlinear. Then

∃C > 0, ∀M > 0, ∃|D| > M such that F (D) ≤ C|D|. (3.2.75)
Chapter 4

Explicit primal-dual algorithm

In this chapter we introduce a numerical scheme based on an algorithm used in image processing [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. We give convergence estimates on both the space continuous formulation and the approximation by finite elements (FEM). The scheme is intended to solve the viscoplastic problem with or without viscosity, the well-posedness of which has been established in the previous chapter. Two algorithms are introduced, one is the first-order primal dual algorithm as in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], the other is the acceleration scheme with the theoretical higher order convergence rate O(1/n 2 ) whereas the former has the rate O(1/n) in terms of the primal-dual gap, where n is the number of iterations. According to Chambolle et al. [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], our inviscid minimization problem

inf u∈H Ω α |u| 2 2 + F (Du) + η |Du| 2 2 -f • u dx, (4.0.1)
with η ≥ 0, can be considered formally as a particular case of the problem inf u∈H h(Ku) + g(u) , (4.0.2)

where g, h are convex, proper, l.s.c. and K : H -→ L 2 (Ω) is a linear map. For us

Ku = Du, h(Ku) = Ω F (Du), g(u) = Ω α |u| 2 2 + η |Du| 2 2 -f • u . (4.0.3)
As before denote by ., . the L 2 duality. The principle of the algorithm is to rewrite (4.0.2) as a primal problem (P) and a dual problem (D) as follows

(P) = inf u∈H h(Ku) + g(u) = inf u∈H sup σ∈L 2 σ, Ku -h * (σ) + g(u) (4.0.4) = inf u∈H sup σ∈L 2 u, K * σ -h * (σ) + g(u) ≥ sup σ∈L 2 inf u∈H u, K * σ -h * (σ) + g(u) = sup σ∈L 2 -g * (-K * σ) -h * (σ) = (D). (4.0.5)
Here K * denotes the adjoint operator of K. In our case of Ku = Du, we have K * σ = -div σ. This formulation corresponds to (2.3.1), (2.3.2), and here the Lagrangian is

L(σ, u) = σ, Ku -h * (σ) + g(u). (4.0.6)
The solution (u, σ) is characterized by ∂ u L(σ, u) 0 and ∂ σ L(σ, u) 0, which can be written ∂g(u) -K * σ, ∂h(Ku) σ. (4.0.7)

Then given u and σ, a way to measure the error from the exact solution is to compute the primaldual gap

h(Ku) + g(u) + g * (-K * σ) + h * (σ) = g(u) + g * (-K * σ) -u, -K * σ + h(Ku) + h * (σ) -σ, Ku ≥ 0, (4.0.8)
that vanishes only when (u, σ) is a solution. We can notice that when the first relation of (4.0.7) is satisfied (which corresponds to the momentum conservation (3.2.47)), the primal-dual gap reduces to (3.2.49).

The idea of the algorithm is to characterize a solution (u, σ) by

Ku ∈ ∂h * (σ), -K * σ ∈ ∂g(u), ⇐⇒ σ ∈ ∂h(Ku), u ∈ ∂g * (-K * σ), ⇐⇒ σ = (Id +r∂h * ) -1 (σ + rKu), u = (Id +τ ∂g) -1 (u -τ K * σ), (4.0.9)
where r > 0, τ > 0 are two parameters. It leads us to the iterative algorithm of [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF],

       σ k+1 = (Id +r∂h * ) -1 (σ k + rK ūk ), u k+1 = (Id +τ ∂g) -1 (u k -τ K * σ k+1 ), ūk+1 = u k+1 + θ(u k+1 -u k ), (4.0.10)
for k ≥ 0, where θ ∈ [0, 1]. According to [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] the stability condition is

rτ K 2 ≤ 1. (4.0.11)
The iteration is completed with initial values u 0 , σ 0 , and we set u -1 = u 0 so that applying the ū formula of (4.0.10) also to k = -1 gives ū0 = u 0 .

The main algorithm is for the choice θ = 1, then ūk = 2u k -u k-1 .

Continuous formulation

Taking into account the definitions (4.0.3), the algorithm (4.0.10) with θ = 1 can be written

   σ k+1 = P r σ k + r(2Du k -Du k-1 ) , u k+1 -u n ∆t -div σ k+1 + u k+1 -u k τ -div ηDu k+1 = -div ηDu n + f n , (4.1.1)
where P r is defined by (2.4.1). One can see the formal consistency of (4.1.1) with

   û -u n ∆t -div σ -div ηDû = -div ηDu n + f n , σ ∈ ∂F (Dû), (4.1.2)
since according to Lemma 2.4.1(c), σ ∈ ∂F (Dû) if and only if σ = P r (σ + rDû).

When the viscosity is positive η > 0, at each iteration (4.1.1) we have to solve an elliptic problem. On the contrary in the inviscid case η = 0 each iteration is fully explicit. The term (u k+1 -u k )/τ can be thought as a relaxation term that enforces the convergence. One can consider also a more regularizing algorithm, replacing it by -div(Du k+1 -Du k )/τ .

We define H = L 2 (Ω) if η = 0, and H = H 1s 0 (Ω) if η > 0, with the norm v 2

H 1s 0 = v 2 L 2 + Dv 2 L 2 .
According to Theorems 3.1.9 and 3.2.9, under some hypotheses on F and Ω there is a solution (û, σ) ∈ H × L 2 (Ω) to (4.1.2), that satisfies F (Du) ∈ L 1 (Ω) and F * (σ) ∈ L 1 (Ω).

In order to get error estimates we consider that v → Dv is a bounded operator, i.e. there exists a constant L ≥ 0 such that

Dv L 2 ≤ L v L 2 . (4.1.3)
Of course this is not true for the inviscid case since then H = L 2 , but when using discrete approximations this becomes true with a constant L depending on the approximation space. In the case of taking the relaxation term as -div(Du k+1 -Du k )/τ the assumption (4.1.3) can be simply replaced by the definition L ≡ 1. In order to understand the mechanism of the estimates, we state a formal result with the assumption (4.1.3).

Proposition 4.1.1 (Formal). If rτ L 2 ≤ 1, then the sequence (u k , σ k ) defined by (4.1.1) verifies u k → û in L 2 as k → ∞, where (û, σ) is the solution to (4.1.2). Moreover, if η > 0 then Du k → Dû in L 2 as k → ∞. If rτ L 2 < 1 then σ k is bounded in L 2 , and σ k+1 -σ k → 0 in L 2 .
Remark: The convergence of σ k to σ is a difficult issue related to the non-uniqueness of σ (as shown in [START_REF] He | Steady Bingham fluid flow in cylindrical pipes: a time dependent approach to the iterative solution[END_REF] in the case of a Bingham fluid). However div σ is unique.

Proof. We use the shorthand notation • for the norm in L 2 (Ω). In the proof we temporary ignore the regularity of u k+1 , σ k+1 and their integrability. By the definition (2.4.1) of P r we have

σ k+1 = (Id +r∂F * ) -1 (σ k + r(2Du k -Du k-1
)), which is equivalent to

∂F * (σ k+1 ) σ k -σ k+1 r + 2Du k -Du k-1 .
It follows that

∀σ ∈ L 2 F * (σ) ≥ F * (σ k+1 ) + σ k -σ k+1 r + 2Du k -Du k-1 , σ -σ k+1 .
Using the identity 2a

• b = |a| 2 + |b| 2 -|a -b| 2 , one gets F * (σ) ≥ F * (σ k+1 ) + σ k -σ k+1 2 2r + σ -σ k+1 2 2r - σ k -σ 2 2r + 2Du k -Du k-1 , σ -σ k+1 . (4.1.4)
Multiplying both sides of the momentum equation of (4.1.1) by u -u k+1 and taking the integral over Ω, we get for any u ∈ H

u k+1 -u n ∆t , u -u k+1 + σ k+1 , Du -Du k+1 + u k+1 -u k τ , u -u k+1 -f n , u -u k+1 -η Du n -Du k+1 , Du -Du k+1 = 0.
Applying again the quadratic identity above, we deduce that for any σ ∈ L 2

u 2 -u k+1 2 -u -u k+1 2 2∆t -f n -div ηDu n + u n ∆t , u -u k+1 + σ k+1 , Du -σ, Du k+1 + σ -σ k+1 , Du k+1 - u k+1 -u k 2 2τ - u -u k+1 2 2τ + u -u k 2 2τ + η 2 Du 2 - η 2 Du k+1 2 - η 2 Du -Du k+1 2 = 0.
Now we have according to the Young inequality, for any λ > 0,

σ k+1 -σ k , Du k -Du k-1 ≤ λ Du k -Du k-1 2 2 + σ k+1 -σ k 2 2λ . (4.1.11) Therefore σ k -σ 2 2r + u k -û 2 2τ + λ Du k -Du k-1 2 2 + σ k -σ, Du k -Du k-1 ≥ σ k+1 -σ 2 2r + u k+1 -û 2 2τ + λ Du k+1 -Du k 2 2 + σ k+1 -σ, Du k+1 -Du k + 1 2r - 1 2λ σ k+1 -σ k 2 + u k+1 -u k 2 2τ -λ Du k+1 -Du k 2 2 + u k+1 -û 2 ∆t + η Du k+1 -Dû 2 .
(4.1.12)

Supposing λ ≥ r, we set

a k = σ k -σ 2 2r + u k -û 2 2τ + λ Du k -Du k-1 2 2 + σ k -σ, Du k -Du k-1 ≥ 0. ( 4
a k ≥ a k+1 + 1 2r - 1 2λ σ k+1 -σ k 2 + 1 2τ L 2 - λ 2 Du k+1 -Du k 2 + u k+1 -û 2 ∆t + η Du k+1 -Dû 2 . (4.1.14)
With the condition τ λL 2 ≤ 1, or in other words

r ≤ λ ≤ 1 τ L 2 , (4.1.15)
we deduce that the sequence {a k } is nonincreasing, and it follows that

∞ k=0 u k+1 -û 2 < ∞, η ∞ k=0 Du k+1 -Dû 2 < ∞. (4.1.16)
Consequently u k → û, and

Du k → Dû if η > 0.
If rτ L 2 < 1 then one can take λ > r satisfying (4.1.15), and according to the Young inequality one has

a k ≥ 1 2r - 1 2λ σ k -σ 2 . (4.1.17)
Since a k is nonincreasing it is bounded, and it follows that σ k -σ is bounded, thus σ k is bounded in L 2 . We also have from (4.1.14) that

1 2r - 1 2λ ∞ k=0 σ k+1 -σ k 2 < ∞, (4.1.18)
which implies that σ k+1 -σ k → 0.

Remark: In the case of taking the relaxation term as -div(Du k+1 -Du k )/τ , all the terms in 1/τ in (4.1.12) are modified, the norms v of a quantity v are replaced by Dv , thus the assumption (4.1.3) is not necessary, we have simply to take L = 1.

Remark:

The result says nothing about the 2018best 2019 choice of the parameters r, τ . Meanwhile the convergence rate depends a lot on this choice, this is discussed in Chapter 5. Another heuristic efficient approach is proposed in [START_REF] Goldstein | Adaptive primaldual hybrid gradient methods for saddle-point problems[END_REF]. A particular feature of this primal 2013dual algorithm is that it can be accelerated when the function g is uniformly convex. In our case

g(u) = Ω α |u| 2 2 + η |Du| 2 2
-f • u is α-uniformly convex (recall that α = 1/∆t). Thus according to [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] one can get convergence at rate O(1/k 2 ) for the primal-dual gap, which means O(1/k) for the velocity error.

Proposition 4.1.2 (Formal). Define the accelerated algorithm as

• Initialization: Choose u 0 , σ 0 , r 0 , τ 0 such that r 0 τ 0 L 2 ≤ 1, and set u -1 = u 0 .

•

Iteration (k ≥ 0): Update r k , τ k , u k , σ k as                        σ k+1 = P r k σ k + r k (Du k + θ k (Du k -Du k-1 )) , u k+1 -u n ∆t -div σ k+1 + u k+1 -u k τ k -div ηDu k+1 = -div ηDu n + f n , τ k+1 = θ k+1 τ k , r k+1 = r k θ k+1 , with θ k+1 = 1 1+2 τ k ∆t . (4.1.19) Then u k → û in L 2 as k → ∞, with u k -û L 2 = O(1/k), where (û, σ) is the solution to (4.1.2). Moreover, if r 0 τ 0 L 2 < 1, then σ k is bounded in L 2 and σ k+1 -σ k L 2 → 0.
Proof. Note that the value of θ 0 is not defined, but does not matter since it appears as a factor of Du 0 -Du -1 = 0. We proceed as in the previous proposition. The inequality (4.1.10) becomes

σ k -σ 2 2r k + u k -û 2 2τ k ≥ σ k+1 -σ 2 2r k + u k+1 -û 2 2τ k + σ k+1 -σ k 2 2r k + u k+1 -u k 2 2τ k + u k+1 -û 2 ∆t + η Du k+1 -Dû 2 + σ k+1 -σ, Du k+1 -(Du k + θ k (Du k -Du k-1 )) . (4.1.20)
The last term can be decomposed as

σ k+1 -σ, Du k+1 -(Du k + θ k (Du k -Du k-1 )) = σ k+1 -σ, Du k+1 -Du k -θ k σ k -σ, Du k -Du k-1 -θ k σ k+1 -σ k , Du k -Du k-1 . (4.1.21)
According to the Young inequality, for any λ k > 0 we have

θ k σ k+1 -σ k , Du k -Du k-1 ≤ λ k θ 2 k Du k -Du k-1 2 2 + σ k+1 -σ k 2 2λ k . (4.1.22) Hence σ k+1 -σ 2 2r k + u k+1 -û 2 2τ k + σ k+1 -σ, Du k+1 -Du k + 1 r k - 1 λ k σ k+1 -σ k 2 2 + u k+1 -u k 2 2τ k + u k+1 -û 2 ∆t + η Du k+1 -Dû 2 ≤ σ k -σ 2 2r k + u k -û 2 2τ k + θ k σ k -σ, Du k -Du k-1 + λ k θ 2 k Du k -Du k-1 2 2 . (4.1.23)
Dividing by τ k and taking λ k ≥ r k , we get

σ k+1 -σ 2 2r k τ k + 1 2τ 2 k + 1 ∆tτ k u k+1 -û 2 + 1 τ k σ k+1 -σ, Du k+1 -Du k + η τ k Du k+1 -Dû 2 + 1 r k - 1 λ k σ k+1 -σ k 2 2τ k + 1 2τ 2 k u k+1 -u k 2 - λ k 2τ k Du k+1 -Du k 2 + λ k τ k Du k+1 -Du k 2 2 ≤ σ k -σ 2 2r k τ k + u k -û 2 2τ 2 k + θ k τ k σ k -σ, Du k -Du k-1 + λ k θ 2 k τ k Du k -Du k-1 2 2 . (4.1.24)
We notice that

r k+1 τ k+1 = r k τ k , r k τ k = r k+1 τ k+1 θ 2 k+1 , 1 τ k = θ k+1 τ k+1 , (4.1.25) 1 2τ 2 k + 1 ∆tτ k = 1 2τ 2 k 1 + 2 τ k ∆t = 1 2τ 2 k θ 2 k+1 = 1 2τ 2 k+1
, and we define {λ k } k≥0 by r 0 ≤ λ 0 ≤ 1/(τ 0 L 2 ) and the update formula

λ k+1 = λ k τ k+1 τ k θ 2 k+1 . (4.1.26)
It follows that λ k+1 /r k+1 = λ k /r k , and λ k ≥ r k for all k. By setting

a k = σ k -σ 2 2r k τ k + u k -û 2 2τ 2 k + λ k θ 2 k τ k Du k -Du k-1 2 2 + θ k τ k σ k -σ, Du k -Du k-1 ≥ 0, (4.1.27) using that Du k+1 -Du k ≤ L u k+1 -u k we obtain a k+1 + η τ k Du k+1 -Dû 2 + 1 r k - 1 λ k σ k+1 -σ k 2 2τ k + 1 τ 2 k L 2 - λ k τ k Du k+1 -Du k 2 2 ≤ a k . (4.1.28) Since λ k τ k L 2 = (λ 0 /r 0 )r k τ k L 2 = λ 0 τ 0 L 2 ≤ 1, we have that {a k } is a nonincreasing nonnegative sequence. As a consequence, a 0 ≥ a k ≥ u k -û 2 2τ 2 k . (4.1.29) Therefore u k → û in L 2 at rate O(τ k ).
According to [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] one has

τ k ∼ ∆t k for large k, thus u k → û at rate O(1/k).
If r 0 τ 0 L 2 < 1, then we can take λ 0 such that r 0 < λ 0 < 1/(τ 0 L 2 ). Since λ k /r k and λ k τ k are constant we have then r k < λ k < 1/(τ k L 2 ) and

a k ≥ 1 r k - 1 λ k σ k -σ 2 2τ k + u k -û 2 2τ 2 k . (4.1.30) Since (1/r k -1/λ k )/τ k is a positive constant, it follows that σ k -σ is bounded, so that σ k is bounded in L 2 . Moreover from (4.1.28) we have ∞ k=0 1 r k - 1 λ k σ k+1 -σ k 2 2τ k < ∞, (4.1.31)
which proves that σ k+1 -σ k → 0.

Remark: Since τ k → 0 as k → ∞, we have r k → ∞ and θ k → 1. Remark: It is not clear how to choose τ 0 , r 0 (satisfying r 0 τ 0 ≤ 1/L 2 ) and k in order to reach a given accuracy at the lowest cost, especially if L is large, which is the case when a space discretisation comes into play as in the next subsection.

Finite element approximation

We suppose now that Ω is a polyhedral domain, and that we have a mesh T h as described in Section 1.5. We define

V h := {v h ∈ C(Ω) such that v h | K is affine for each cell K ∈ T h , v h | ∂Ω = 0}, (4.2.1) Λ h := {σ h ∈ L ∞ (Ω) such that σ h | K is constant for each K ∈ T h }. (4.2.2)
Following Theorem 3.1.9, we consider the problem of finding (û h , σh )

∈ V h × Λ h such that Ω ûh -u n ∆t • v h + Ω σh : Dv h + η Ω Dû h : Dv h = η Ω Du n : Dv h + Ω f n • v h ∀v h ∈ V h , (4.2.3) σh ∈ ∂F (Dû h ) a.e. in Ω. (4.2.4)
As in Chapter 3 we have to consider the integrated version, which is: find ûh ∈ V h such that

∀v h ∈ V h 1 ∆t ûh -u n , v h -û h + Ω F (Dv h )- Ω F (Dû h )+η Dû h -Du n , Dv h -Dû h ≥ f n , v h -û h , (4.2
.5) where ., . still denotes the L 2 duality. Proposition 4.2.1. We assume that η ≥ 0, and that F is convex and finite everywhere. Then there exists one and only one solution ûh to the problem (4.2.5), and moreover it satisfies

ûh = argmin v h ∈V h J(v h ), (4.2.6)
where J is defined by

J(v) = 1 ∆t Ω 1 2 |v| 2 + η 2 Ω |Dv| 2 + Ω F (Dv) -G(v), (4.2.7) with G(v) = f n + u n ∆t , v + η Du n , Dv .
Proof. We note that V h ⊂ H 1s 0 (Ω). Moreover for v h ∈ V h , Dv h ∈ L ∞ . Since F is convex and finite everywhere on a finite-dimensional space, it is continuous and bounded on bounded sets. It follows that v h → Ω F (Dv h ) is a finite valued convex function on V h . Since V h has finite dimension, this functional is thus continuous. Since G is a linear form on V h , the result follows from the classical Proposition 2.1.4(a). Proposition 4.2.2. We assume that η ≥ 0, and that in the viscous case (η > 0) F satisfies Hypothesis 1, 2, 3, whereas in the inviscid case (η = 0), F satisfies Hypothesis 1, 3, 4, 5. We denote

H = H 1s 0 (Ω) if η > 0, H = L 2 (Ω) if η = 0.
Let û ∈ H be the solution to the problem (3.1.28) if η > 0, respectively (3.2.34) if η = 0, and ûh be the solution to the problem (4.2.5). Then

ûh -û 2 L 2 ∆t + η Dû h -Dû 2 L 2 ≤ 2 R û, (4.2.8)
and

Ω F (Dû h ) - Ω F (Dû) ≤ 6 M R û, ( 4

.2.9)

where

R û = inf v h ∈V h √ M v h -û 2 L 2 ∆t + η Dv h -Dû 2 L 2 1/2 + max 0, Ω F (Dv h ) - Ω F (Dû) , (4.2.10) M = 2|Ω|(E + F (0)) + ∆t f n + u n ∆t 2 L 2 + η Du n 2 L 2 . (4.2.11)
Proof. According to Propositions 3.1.7 and 3.2.7, û is the unique function in H such that

∀v ∈ H 1 ∆t û, v -û + η Dû, Dv -Dû + ψ(v) ≥ ψ(û) + G(v -û), (4.2.12)
where ψ is defined on H either by (3.1.22) or by (3.2.33). Thus

ψ(v) is either Ω F (Dv) if η > 0, or if η = 0, Ω F (Dv) for v ∈ K F 0 (Ω), +∞ if v ∈ K F 0 (Ω). We have V h ⊂ H 1s 0 (Ω) ⊂ K F 0 (Ω), and ψ(v) is equal to Ω F (Dv) for v ∈ V h . The problem (4.2.12) is also equivalent to minimizing J over H, with ∀v ∈ H, J(v) = 1 ∆t Ω 1 2 |v| 2 + η 2 Ω |Dv| 2 + ψ(v) -G(v). (4.2.13)
This value of J reduces to (4.2.7) when v ∈ V h . The problem (4.2.5) with solution ûh ∈ V h can be written also

∀v h ∈ V h 1 ∆t ûh , v h -ûh + η Dû h , Dv h -Dû h + ψ(v h ) ≥ ψ(û h ) + G(v h -ûh ). ( 4

.2.14)

In order to prove the estimates, we follow [START_REF] Bouchut | Convergence of conforming approximations for inviscid incompressible Bingham fluid flows and related problems[END_REF] and we denote

f = f n + u n ∆t . Setting v h = 0
v h ∈ V h 1 ∆t ûh , û-û h +η Dû h , Dû-Dû h +ψ(û)-ψ(û h )+R(v h ) ≥ f, û-û h +η Du n , Dû-Dû h , (4.2.18) with R(v h ) = 1 ∆t ûh , v h -û +η Dû h , Dv h -Dû +ψ(v h )-ψ(û)-f, v h -û -η Du n , Dv h -Dû . (4.2.19)
Hence, taking the infimum gives

1 ∆t ûh , û -ûh + η Dû h , Dû -Dû h + ψ(û) -ψ(û h ) + inf v h ∈V h R(v h ) ≥ f, û -ûh + η Du n , Dû -Dû h . ( 4
R(v h ) ≤ 1 ∆t ûh v h -û + η Dû h Dv h -Dû + ψ(v h ) -ψ(û) + + f v h -û + η Du n Dv h -Dû ≤ ûh 2 ∆t + η Dû h 2 1/2 v h -û 2 ∆t + η Dv h -Dû 2 1/2 + ∆t f 2 + η Du n 2 1/2 v h -û 2 ∆t + η Dv h -Dû 2 1/2 + ψ(v h ) -ψ(û) + ≤ 2 √ M v h -û 2 ∆t + η Dv h -Dû 2 1/2 + ψ(v h ) -ψ(û) + .
(4.2.24) Taking the infimum over v h ∈ V h and taking into account that û ∈ K F 0 (Ω) so that ψ(û) = F (Dû), it yields inf

v h ∈V h R(v h ) ≤ 2 R û. ( 4
ψ(û h ) -ψ(û) ≤ f, ûh -û + η Du n , Dû h -Dû + 1 ∆t ûh , û -ûh + η Dû h , Dû -Dû h + inf v h ∈V h R(v h ) ≤ ∆t f 2 + η Du n 2 1/2 + ûh 2 ∆t + η Dû h 2 1/2 × ûh -û 2 ∆t + η Dû h -Dû 2 1/2 + inf v h ∈V h R(v h ) ≤ 2 √ M ûh -û 2 ∆t + η Dû h -Dû 2 1/2 + inf v h ∈V h R(v h ). ( 4

.2.26)

Similarly using (4.2.21),

ψ(û) -ψ(û h ) ≤ 1 ∆t û, ûh -û + η Dû, Dû h -Dû -f, ûh -û -η Du n , Dû h -Dû ≤ ∆t f 2 + η Du n 2 1/2 + û 2 ∆t + η Dû 2 1/2 × ûh -û 2 ∆t + η Dû h -Dû 2 1/2 . ( 4
ψ(û h ) -ψ(û) ≤ 2 √ M ûh -û 2 ∆t + η Dû h -Dû 2 1/2 + inf v h ∈V h R(v h ). (4.2.28)
Taking into account (4.2.22) we deduce that

ψ(û h ) -ψ(û) ≤ 2 M inf vh ∈V h R(v h ) + inf v h ∈V h R(v h ). (4.2.29) But taking v h = 0 in (4.2.19) we have inf v h ∈V h R(v h ) ≤ R(0) = - 1 ∆t ûh , û -η Dû h , Dû + ψ(0) -ψ(û) + f, û + η Du n , Dû ≤ M + |Ω|(E + F (0)) + M ≤ 5 2 M. ( 4
ψ(û h ) -ψ(û) ≤ 2 √ M + 5 2 M inf v h ∈V h R(v h ) ≤ √ 2 2 + 5 2 M R û ≤ 6 M R û. ( 4
(Ω) if η > 0, or in L 2 (Ω) if η = 0. Moreover Dû h
Dû in weak L 1 (Ω) and Ω F (Dû h ) → Ω F (Dû).

Proof. Since û ∈ H and û ∈ K F 0 (Ω) with ψ(û) < ∞ if η = 0, for any ε > 0 there exists w ∈ C ∞ c (Ω) such that w -û H ≤ ε and

F (Dw) -F (Dû) ≤ ε. (4.2.32)
Indeed if η > 0, by definition of H 1s 0 we can find a sequence w k ∈ C ∞ c such that w k -û H 1s 0 → 0. Then since F is subquadratic one has that F (Dw k ) → F (Dû), and the result follows. In the case η = 0, the result follows from Lemma 3.2.4. Then, w being given in C ∞ c , for h small enough there is a function Remark: In Corollary 4.2.3 the convergence holds as h → 0 at ∆t fixed. The analysis in the case of small ∆t for a time-dependent problem is slighly different and is done in [START_REF] Bouchut | Convergence of conforming approximations for inviscid incompressible Bingham fluid flows and related problems[END_REF], at least in the case when F attains its minumum at 0, i.e. E + F (0) = 0.

v h ∈ V h such that v h -w H 1 (Ω) ≤ ε and F (Dv h ) -F (Dw) ≤ ε, ( 4
Proposition 4.2.4 (Discrete Euler-Lagrange equations). We assume that η ≥ 0, and that F is convex, finite everywhere and lower bounded. Then the solution ûh ∈ V h to the problem (4.2. Now in the general case when F is not differentiable, for any ε > 0 we can consider the Moreau envelope F ε of F of Proposition 2.1.4. Then F ε is continuously differentiable, and converges monotonically to F as ε → 0. We can apply the existence result to F ε , thus there exist

(û ε h , σε h ) ∈ V h × Λ h satisfying Ω ûε h -u n ∆t • v h + Ω σε h : Dv h + η Ω Dû ε h : Dv h = η Ω Du n : Dv h + Ω f n • v h ∀v h ∈ V h , (4.2.39)
σε h ∈ ∂F ε (Dû ε h ) a.e. in Ω. (4.2.40)
We have that ûε h also satisfies (4.2.5) with F replaced by F ε . Since F is lower bounded by a constant -E one has that F ε is lower bounded by the same constant. Thus taking v h = 0 as test function, we can do the same estimate (4.2.17) as in Proposition 4.2.2, and get

ûε h 2 2∆t + η 2 Dû ε h 2 ≤ |Ω|(E + F ε (0)) + ∆t 2 f 2 + η 2 Du n 2 . (4.2.41)
We have F ε (0) ≤ F (0), thus ûε h is bounded in H independently of ε. But since ûε h ∈ V h and V h is finite dimensional, we deduce that ûε h is bounded in V h , and thus that Dû ε h is bounded in L ∞ . Since ∂F (D) remains bounded when D lies in a bounded set, using Lemma 3.1.1 and Proposition 2.1.4(c)(e) we deduce that F ε (D) remains bounded when D lies in a bounded set, independently of ε ≤ ε 0 . It follows that σε h is bounded in L ∞ independently of ε. Extracting a subsequence if necessary, (V h and Λ h are finite dimensional), we get ûε

h → ûh ∈ V h , σε h → σh ∈ Λ h as ε → 0. The spaces V h and Λ h being finite dimensional, we have thus Dû ε h → Dû h in L ∞ , σε h → σh in L ∞ .
Thus we can pass to the limit in (4.2.39) and get (4.2.3). Then according to (4.2.40) one has for all D

F ε (D) ≥ F ε (Dû ε h ) + σε h : (D -Dû ε h ) a.e. in Ω. (4.2.42)
The 'almost everywhere' means indeed that it holds on all cells K ∈ T h , since the quantities are constant on each cell. For ε ≤ ε 0 we have

F ε (Dû ε h ) ≥ F ε 0 (Dû ε h ). Thus letting ε → 0 we get F (D) ≥ F ε 0 (Dû h ) + σh : (D -Dû h ) a.e. in Ω. (4.2.43)
Finally we let ε 0 → 0 and get that σh ∈ ∂F (Dû h ) a.e., i.e. (

We can now give a complement to Corollary (

F * (σ h ) → Ω F * (σ).
Proof. Since σh satisfies (4.2.4), we have σh :

Dû h = F (Dû h ) + F * (σ h ) a.e. in Ω. (4.2.44)
The functions involved in this identity are indeed constant in each cell

K ∈ T h . If follows that F * (σ h ) ∈ L ∞ (Ω). Then taking v h = ûh in (4.2.
3) and using (4.2.44) we obtain ûh

2 ∆t + η Dû h 2 + Ω F (Dû h ) + Ω F * (σ h ) = η Du n , Dû h + f n + u n ∆t , ûh . (4.2.45)
Since F and F * are lower bounded, this gives bounds on ûh H , F (Dû h ), F * (σ h ) independent of h. By Lemma 3.1.3(a) we deduce that σh is bounded in L 2 (Ω) independently of h. Thus after extraction of a subsequence, there is some σ ∈ L 2 (Ω) such that σh σ in weak L 2 . Since σ → F * (σ) is convex and l.s.c. on L 2 , thus it is also l.s.c. on weak L 2 . Therefore passing to the limit in (4.2.45) and using Corollary (4.2.3) we obtain

û 2 ∆t + η Dû 2 + Ω F (Dû) + Ω F * (σ) ≤ η Du n , Dû + f n + u n ∆t , û . (4.2.46)
Thus if η = 0 we get (3.2.49) and (3.2.46). Next, for ϕ ∈ C ∞ c (Ω) and ε > 0, as in Corollary (4.2.3), for h small enough there exists v h ∈ V h such that v h -ϕ H 1 ≤ ε. The formulation (4.2.3) ensures then that for h small enough

Ω ûh -u n ∆t • ϕ + Ω σh : Dϕ + η Ω Dû h : Dϕ -η Ω Du n : Dϕ - Ω f n • ϕ ≤ Cε, (4.2.47)
where C is a constant independent of h and ε. Letting h → 0 we get

Ω û -u n ∆t • ϕ + Ω σ : Dϕ + η Ω Dû : Dϕ -η Ω Du n : Dϕ - Ω f n • ϕ ≤ Cε. (4.2.48)
Since this holds for any ε > 0 we conclude that

û -u n ∆t -div σ -η div Dû = -η div Du n + f n , (4.2.49)
in the sense of distributions in Ω. This proves the first equation of (3.1.32) in the case η > 0 or (3.2.47) in the case η = 0. Then we have to prove the second equation of (3.1.32) in the case η > 0, or (3.2.48) in the case η = 0, which is any case writes σ ∈ ∂F (Dû) a.e. in Ω. We have proved in Theorem 3.2.9 that for η = 0 this condition was consequence of the others, thus it is not necessary to prove it. In the case η > 0 we can use the same argument. By density of C ∞ c in H 1s 0 we can take û as test function in (4.2.49), giving

û 2 ∆t + η Dû 2 + σ, Dû = η Du n , Dû + f n + u n ∆t , û . (4.2.50)
Comparing to (4.2.46) we deduce that

Ω F (Dû) + Ω F * (σ) ≤ σ, Dû , (4.2.51)
which proves that σ ∈ ∂F (Dû) a.e. in Ω. Finally we have F (Dû) + F * (σ) = σ : Dû a.e. in Ω, and we deduce that there is equality in (4.2.46) (in the case η = 0 this was already proved in Theorem 3.2.9). Therefore comparing to (4.2.45) we conclude that Ω F * (σ h ) → Ω F * (σ).

Next we prove the convergence of the algorithm (4.1.1), corresponding to θ = 1, applied at the discrete level. This means that we define a sequence (

u k h , σ k h ) ∈ V h × Λ h for all k ≥ 0 by the iteration formula: for k ≥ 0    σ k+1 h = P r σ k h + r(2Du k h -Du k-1 h ) a.e. in Ω, u k+1 h -u n ∆t , v h + σ k+1 h , Dv h + u k+1 h -u k h τ , v h + η Du k+1 h -Du n , Dv h = f n , v h , ∀v h ∈ V h , (4.2.52
) where P r is defined by (2.4.1) and ., . still denotes the L 2 scalar product. Here r > 0 and τ > 0 are parameters. We observe that the first equation in (4.2.52) deals with data that are constant in each cell K ∈ T h . Thus this formula is applied independently in each cell. By the Riesz theorem the second equation determines a unique solution

u k+1 h ∈ V h . Given u 0 h ∈ V h and σ 0 h ∈ Λ h we set u -1 h = u 0 h . Then the above iteration formulas define (u k h , σ k h ) ∈ V h × Λ h
for all k ≥ 0. One can see the formal consistency with (4.2.3), (4.2.4).

The operator v h → Dv h being linear, it is bounded on V h since V h is finite dimensional,

Dv h L 2 ≤ L h v h L 2 ∀v h ∈ V h . ( 4
∂F * (σ k+1 h ) σ k h -σ k+1 h r + 2Du k h -Du k-1 h a.e. in Ω. (4.2.55)
We deduce that F * (σ k+1 h ) < ∞ and thus since it is piecewise constant, F * (σ k+1 h ) ∈ L ∞ and

∀σ h ∈ Λ h F * (σ h ) ≥ F * (σ k+1 h ) + σ k h -σ k+1 h r + 2Du k h -Du k-1 h , σ h -σ k+1 h . (4.2.56)
Using the quadratic identity, it follows that for any

σ h ∈ Λ h F * (σ h ) ≥ F * (σ k+1 h ) + σ k h -σ k+1 h 2 2r + σ h -σ k+1 h 2 2r - σ k h -σ h 2 2r + 2Du k h -Du k-1 h , σ h -σ k+1 h . (4.2.57)
Taking

v h = u h -u k+1 h
in the momentum equation of (4.2.52) we get for any

u h ∈ V h u k+1 h -u n ∆t , u h -u k+1 h + σ k+1 h , Du h -Du k+1 h + u k+1 h -u k h τ , u h -u k+1 h +η Du k+1 h -Du n , Du h -Du k+1 h -f n , u h -u k+1 h = 0. ( 4
∀σ h ∈ Λ h F * (σ h ) ≥ F * (σ h ) + Dû h , σ h -σh . (4.2.59)
Then we follow the proof of Proposition 4.1.1. We take some λ such that

r ≤ λ ≤ 1 τ L 2 h , (4.2.60)
which is possible according to (4.2.54). Defining a k as (4.1.13), we obtain the inequality (4.1.14), and the estimates (4.1.16). We deduce that u k h → ûh in L 2 . In the case of strict inequality rτ L 2 h < 1, one can choose λ with strict inequalities in (4.2.60). Then (4.1.17), (4.1.18) hold, and it follows that

σ k h is bounded in L 2 and σ k+1 h -σ k h → 0 as k → ∞.
Since Λ h is finite dimensional, after extraction of a subsequence one has σ k h → σh ∈ Λ h (that may be different from the one previously considered). Passing to the limit in (4.2.52) and using that u k h → ûh and σ k+1 h -σ k h → 0 we obtain (4.2.3), (4.2.4).

One can also consider the accelerated algorithm (4. 1.19), that can be written as defining

(u k h , σ k h ) ∈ V h × Λ h by • Initialization: Choose u 0 h ∈ V h , σ 0 h ∈ Λ h , r 0 , τ 0 such that r 0 τ 0 L 2 h ≤ 1, and set u -1 h = u 0 h . • Iteration (k ≥ 0): Update r k , τ k , u k h , σ k h as                        σ k+1 h = P r k σ k h + r k (Du k h + θ k (Du k h -Du k-1 h )) a.e. in Ω, u k+1 h -u n ∆t , v h + σ k+1 h , Dv h + u k+1 h -u k h τ k , v h + η Du k+1 h -Du n , Dv h = f n , v h , ∀v h ∈ V h , τ k+1 = θ k+1 τ k , r k+1 = r k θ k+1 , with θ k+1 = 1 1+2 τ k ∆t . (4.2.61)
Theorem 4.2.7. We assume that η ≥ 0, and that F is convex, finite everywhere and lower bounded. We denote by ûh ∈ V h the solution to (4.2.5) obtained in Proposition 4.2.1, and we assume that

r 0 τ 0 L 2 h ≤ 1, (4.2

.62)

where L h is such that (4.2.53) holds. Then the sequence defined by the algorithm (4.2.61) verifies

u k h → ûh as k → ∞ with u k h -ûh L 2 = O(1/k). Moreover if there is strict inequality in (4.2.62) then σ k h is bounded in Λ h ,

and up to extraction of a subsequence one has

σ k h → σh in Λ h , where (û h , σh ) ∈ V h × Λ h solves (4.2.3), (4.2.4).
Proof. It is identical to the one of the previous theorem, taking into account the arguments in the proof of Proposition 4.1.2.

Remark: For the inviscid case η = 0, we would like to use a really explicit algorithm, thus not having to invert the mass matrix in order to get u k+1 h from the second line of (4.2.52). With the notation of Section 1.5 we thus introduce the scalar product over

V h × V h u h , v h h = Ω u h • v h . ( 4

.2.63)

We then use this scalar product instead of the L 2 scalar product on V h × V h . In particular, the second line of (4.2.52) is changed into

u k+1 h -u n ∆t , v h h + σ k+1 h , Dv h + u k+1 h -u k h τ , v h h + η Du k+1 h -Du n , Dv h = f n , v h h , ∀v h ∈ V h , (4.2
.64) where u n and f n are approximations in V h of u n and f n respectively. Note that we do not modify the scalar product of matrices: the terms involving D are unmodified. Then with this formulation there is no linear system to invert at all, each iteration is the multiplication of the vector of unknowns by a (sparse) matrix. The definition of L h has to be modified however, as

Dv h L 2 ≤ L h v h , v h 1/2 h , ∀v h ∈ V h . ( 4

.2.65)

Remark: It is not straightforward to use a higher order approximation in space, for example by taking V h = P 2 for u and Λ h = P 1 for σ, since the projection step (first line of (4.2.52)) would not operate in the discrete space Λ h . Thus one would have to consider a further projection to Λ h as

σ k+1 h , σ h = P r σ k h + r(2Du k h -Du k-1 h ) , σ h , ∀σ h ∈ Λ h . (4.2.66)
Then we cannot ensure the condition σh ∈ ∂F (Dû h ) a.e. for the limit (û h , σh ) obtained as k → ∞.

x 0 x 1

x 2 x Nx x Nx+1 h/2 h h h h/2
= x 0 < h 2 = x 1 < 3h 2 = x 2 < • • • < N x - 1 2 h = x Nx < L = x Nx+1 , (5.1.11)
where h = L N x , see Figure 5.1. The two boundary cells have size half of the size of the internal cells. We denote

u n i ≈ u(t n , x i ), ρ n i ≈ ρ(t n , x i ).
In our scheme the unknown function u(t, x) is reconstructed simultaneously from the values u h = (u i ) i∈{0,••• ,Nx+1} using a continuous piecewise affine reconstruction denoted by ûh , and a piecewise constant reconstruction denoted by ūh . Thus for all u h ∈ R Nx+2 we have

ûh ∈ C([0, L]), ûh is affine for each interval [x i , x i+1 ], ûh (x i ) = u i ∀i ∈ {0, • • • , N x + 1}, (5.1.12) ūh ∈ L 1 loc ((0, L)), ūh (x) = u i ∀x ∈ ((i -1)h, ih), i ∈ {1, • • • , N x }. (5.1.13)
We use the same notation for ρ h , ρh . Numerical scheme. The numerical scheme approximating (5.1.1) is given by 1

-Initialisation of u 0 h ∈ R Nx+2 , ρ 0 h ∈ R Nx+2 : u 0 i = u init (x i ) ∀i ∈ {0, • • • , N x + 1}, ρ 0 i = ρ init (x i ) ∀i ∈ {0, • • • , N x + 1}.
2-Finite volume step: Suppose that (u n h , ρ n h ) are known. We define u

n+ 1 2 h , ρ n+ 1 2 h which approximate (u n+ 1 2 , ρ n+ 1 2
), the solution to 1.3.1:

U n+ 1 2 i = U n i - ∆t h F i+ 1 2 -F i-1 2 , (5.1.14)
where U n i := (ρ n i , ρ n i u n i ) and similarly with U

n+ 1 2 i
. We consider here only first-order explicit three points schemes where

F i+ 1 2 = F (U n i , U n i+1 ). (5.1.15)
The function F (U l , U r ) is called the numerical flux. There are lots of different choices for F , some of being well-known and widely used, as the upwind, Lax-Friedrichs or Siliciu schemes. We will not discuss further this here, and just apply a known numerical flux such as the Suliciu one, see [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF].

There is however an associated CFL condition (for Courant, Friedrichs, Levy [START_REF] Courant | On the partial difference equations of mathematical physics[END_REF]) on the timestep to prevent the blow up of the numerical values, under the form ∆ta ≤ h,

where a is an approximation of the propagation speed.

We have to use compatible boundary conditions between the finite volume and finite element steps. We use the Dirichlet boundary condition in the finite element step and the "wall" condition for the finite volume step. This means that, the basic unknowns being for i = 1, . . . , N x , we set before the finite volume step

ρ n 0 = ρ n 1 , ρ n Nx+1 = ρ n Nx , u n 0 = -u n 1 , u n Nx+1 = -u n Nx .
(5.1.17)

3-Finite element step: for illustrative purpose we use the regularisation method. Suppose that (u

n+ 1 2 h , ρ n+ 1 2 h
) are known. We define u n+1 h , ρ n+1 h which approximate u n+1 , ρ n+1 the solution to (1.3.2), as

ρ n+1 h = ρ n+1/2 h and L 0 ρ n+1 h ūn+1 h - ūn+ 1 2 h ∆t vh dx +σ 0 L 0 ∂ x ûn+1 h |∂ x ûn+1 h | ∂ x vh dx = L 0 fh vh dx ∀v h ∈ V h . (5.1.18)
Using the regularisation method, we replace this by the approximation, for a small positive ε,

L 0 ρ n+1 h ūn+1 h - ūn+ 1 2 h ∆t vh dx +σ 0 L 0 ∂ x ûn+1 h |∂ x ûn+1 h | 2 + ε 2 ∂ x vh dx = L 0 fh vh dx ∀v h ∈ V h . (5.1.19)
Note that ε is chosen "optimally" as

ε ∼ 10 -2 h 2 ||∂ x u|| 2 ρ σ 0 1 T , (5.1.20)
see [START_REF] Lusso | Two-dimensional simulation by regularization of free surface viscoplastic flows with Drucker-Prager yield stress and application to granular collapse[END_REF]. We use the fixed point method to get u n+1 h . We initialize

u n+1,0 h = u n+1/2 h
, and when

u n+1,k h is known we compute u n+1,k+1 h by L 0 ρ n+1 h ūn+1,k+1 h - ūn+ 1 2 h ∆t vh dx +σ 0 L 0 ∂ x ûn+1,k+1 h |∂ x ûn+1,k h | 2 + ε 2 ∂ x vh dx = L 0 fh vh dx ∀v h ∈ V h .
(5.1.21) The unknown ûn+1,k+1 h has to vanish at 0 and L, thus it remains only the unknown values u n+1,k+1 i for i = 1, . . . , N x . By choosing vh = vj h the function such that vj

h (x k ) = 1 if k = j and 0 if k = j, it yields L 0 ρ n+1 h ūn+1,k+1 h - ūn+ 1 2 h ∆t vj h = Nx i=1 hρ n+1/2 i u n+1,k+1 i -u n+ 1 2 i ∆t v j i = hρ n+1/2 j u n+1,k+1 j -u n+ 1 2 j ∆t , j = 1, . . . , N x , (5.1.22)
L 0 ∂ x ûn+1,k+1 h |∂ x ûn+1,k h | 2 + ε 2 ∂ x vj h =                                u n+1,k+1 1 -u n+1,k+1 2
h 2 ε 2 + u n+1,k 1 -u n+1,k 2 2 + u n+1,k+1 1 h 2 4 ε 2 +(u n+1,k 1 ) 2 j = 1, u n+1,k+1 j -u n+1,k+1 j+1 h 2 ε 2 + u n+1,k j+1 -u n+1,k j 2 + u n+1,k+1 j -u n+1,k+1 j-1 h 2 ε 2 + u n+1,k j -u n+1,k j-1 2 2 ≤ j ≤ N x -1, -u n+1,k+1 Nx-1 +u n+1,k+1 Nx h 2 ε 2 + u n+1,k Nx -u n+1,k Nx-1 2 + u n+1,k+1 Nx h 2 4 ε 2 +(u n+1,k Nx ) 2 j = N x , (5.1.23) L 0 f vj h =                          h f 0 12 + f 1 2 + f 2 6 j = 1, h f j-1 6 + 2f j 3 + f j+1 6 2 ≤ j ≤ N x -1, h f Nx-1 6 + f N 2 + f Nx+1 12 j = N x .
(5.1.24)

To find u n+1,k+1 h , one has to solve the system of linear equations AX = F as

          a 1 -b 1 -b 1 a 2 -b 2 -b 2 a 3 -b 3 . . . -b Nx-2 a Nx-1 -b Nx-1 -b Nx-1 a Nx                       u n+1,k+1 1 . . . . . . . . . . . . u n+1,k+1 Nx             =             F 1 . . . . . . . . . . . . F Nx             , (5.1.25)
where

b i = σ 0 ∆t h 2 ε 2 + u n+1,k i+1 -u n+1,k i 2 , i = 1, . . . , N x -1, a 1 = σ 0 ∆t h 2 ε 2 + u n+1,k 2 -u n+1,k 1 2 + σ 0 ∆t h 2 4 ε 2 + u n+1,k 1 2 + hρ n+1/2 1 = b 1 + σ 0 ∆t h 2 4 ε 2 + u n+1,k 1 2 + hρ n+1/2 1 , a Nx = σ 0 ∆t h 2 ε 2 + u n+1,k Nx-1 -u n+1,k Nx 2 + σ 0 ∆t h 2 4 ε 2 + u n+1,k Nx 2 + hρ n+1/2 Nx = b Nx-1 + σ 0 ∆t h 2 4 ε 2 + u n+1,k Nx 2 + hρ n+1/2 Nx , a i = σ 0 ∆t h 2 ε 2 + u n+1,k i-1 -u n+1,k i 2 + σ 0 ∆t h 2 ε 2 + u n+1,k i -u n+1,k i+1 2 + hρ n+1/2 i = b i-1 + b i + hρ n+1/2 i , i = 2, . . . , N x -1, F 1 = h∆t f 0 12 + f 1 2 + f 2 6 + hρ n+1/2 1 u n 1 , F Nx = h∆t f Nx-1 6 + f N 2 + f Nx+1 12 + hρ n+1/2 Nx u n Nx , F i = h∆t f i-1 6 + 2f i 3 + f i+1 6 + hρ n+1/2 i u n i , i = 2, . . . , N x -1.
This tridiagonal system is resolved classically writing

u n+1,k+1 i = λ i u n+1,k+1 i+1 + r i , for i = 1, . . . , N x -1.
(5.1.26)

The system (5.1.25) indeed gives the recursive formulas

λ 1 = b 1 a 1 , λ i = b i a i -b i-1 λ i-1
for i = 2, . . . , N x -1, (5.1.27)

r 1 = F 1 a 1 , r i = F i + b i-1 r i-1 a i -b i-1 λ i-1 for i = 2, . . . , N x .
(5.1.28)

We find u n+1,k+1 by setting u n+1,k+1 Nx = r Nx and applying (5.1.26). The previous procedure enables to find u n+1,k+1 from the knowledge of u n+1,k . We stop the iteration procedure when a stopping criterion is satisfied, u n+1,k+1 -u n+1,k ≤ ε tol .

The numerical results are shown on seen in the left subfigure of Figure 5.2 there is no significant difference between the computed ρ (red line) and the exact ρ (green line) at the final time T = 1. We have the same conclusion for u. In fact, the L 1 error in this case is 3.872e-3 i.e. (ρ

N T h -ρ(T ) + ûN T h -u(T )) L 1 (0,L) =3.872e-3.

2D compressible Bingham model

Steady case

In this subsection the primal-dual algorithm is considered to solve the two-dimensional steady Bingham equation

αu -div √ 2 Du |Du| = f, for (x, y) ∈ (-1, 1) × (-1, 1). (5.2.1)
An analytical solution is built under the form u(x, y) = Φ(r) -y x where r = x 2 + y 2 . Then |Du|/ √ 2 = r|∂ r Φ|/2 and we have the equation

αΦ(r) - ∂ r (sgn ∂ r Φ) r - 2sgn ∂ r Φ r 2 -y x = f. (5.2.2)
We consider For 0 ≤ r ≤ 1/6,

sgn ∂ r Φ = (12r -36r 2 ) 2 , Φ = 1, f = α -2 × 12 2 (1 -3r)(2 -9r) -y x . For 1/6 ≤ r ≤ 1/3, sgn ∂ r Φ = 1, Φ = 6r, f = 6αr -2/r 2 -y x . For 1/3 ≤ r ≤ 1/2, sgn ∂ r Φ = cos(π(6r -2)), Φ = 2, f = 2α + 6π sin(π(6r -2)) r - 2 cos(π(6r -2)) r 2 -y x .
For 1/2 ≤ r ≤ 5/6, sgn ∂ r Φ = -1, Φ = 5 -6r, f = α(5 -6r) + 2/r 2 -y x .

For 5/6 ≤ r ≤ 1,

sgn ∂ r Φ = - 1 + cos(π(6r -5)) 2 , Φ = 0, f = -3π sin(π(6r -5)) r + 1 + cos(π(6r -5)) r 2 -y x .
For 1 ≤ r, sgn ∂ r Φ = 0, Φ = 0, f = 0.

Space discretisation: We consider a rectangular domain Ω = (-L x , L x ) × (-L y , L y ), and we discretize it in both directions similarly as we did in the 1d case, see Figure 5.1. We denote by N x + 2 and N y + 2 the numbers of points in the horizontal and vertical directions. We define the spatial steps h x = 2L x N x , h y = 2L y N y and the grid points as

x 0 = -L x , x Nx+1 = L x , y 0 = -L y , y Ny+1 = L y , (5.2.3)
x i = -L x + i - 1 2 h x , i = 1, • • • , N x , (5.2.4)
y i = -L y + j - 1 2 h y , j = 1, • • • , N y , (5.2.5)
x i+ 1 2 = 1 2 (x i + x i+1 ), i = 0, • • • , N x , (5.2.6
)

y i+ 1 2 = 1 2 (y i + y i+1 ), j = 0, • • • , N y .
(5.2.7)

For each couple of integers i,j such that 0 ≤ i ≤ N x + 1, 0 ≤ j ≤ N y + 1, we define the rectangular cells as R i,j = (x i , x i+1 ) × (y j , y j+1 ).

(5.2.8)

In order to get a mesh of triangles, these cells are cut in two along one diagonal, see Figure 5.3(left) giving the cells K ∈ T h .

We use the primal-dual algorithm, under the form (4.2.52) modified as (4.2.64) (mass lumped scheme). We have here F (D) = √ 2|D|, thus u k h and σ k h being known we look for u k+1 h and σ k+1 h solution to

σ k+1 h = P r (σ k h + r(2Dû k h -Dû k-1 h )), (5.2.9) α Ω ūk+1 h • vh dx + Ω σ k+1 h : Dv h dx + Ω ūk+1 h -ūk h τ • vh dx = Ω fh • vh dx ∀v h ∈ V h , (5.2.10)
where

P r (σ) =    σ if |σ| ≤ √ 2, √ 2 σ |σ| otherwise.
(5.2.11)

(x i-1 , y i ) (x i+1 , y i ) (x i+1 , y i-1 ) (x i , y i-1 ) (x i , y i+1 ) (x i-1 , y i+1 ) Qn K 5 K 6 K 1 K 2 K 3 K 4 n Figure
Dv h L 2 ≤ L h Ω |v h | 2 1/2
, ∀v h ∈ V h .

(5.2.17)

Since obviously L h is of the order of 1/h, the stability condition means that rτ ≤ Ch 2 for some appropriate constant C. We test various values of the product rτ in the case h x = h y = 1/20 with the stopping criteria ε tol = 10 -6 and r = τ = √ Ch, α = 1. As can be seen in Table 5.1, the error increases as C decreases, even though it decreases significantly in the case C = 1/8 in comparison with C = 1/6, and then decreases slowly with smaller values of C. On the other hand, the number of iteration increases proportionally to the increase of 1/C. Besides, the algorithm does not converge for values higher than 1/6. Thus we choose C = 1/8 as the optimal value for h = 1/20. Testing with h = 1/20, 1/40, 1/80 the optimal C is found to be 1/6 with respect to both the number of iterations and the error. Optimal choice for the stopping criterion ε tol for the primal-dual algorithm In the first primal-dual algorithm we have set a stopping criterion ûk+1 h -ûk h L 2 ≤ ε tol for the iteration loop. A smaller ε tol means a smaller error, but more iterations. According to Table 5.2 showing the error, we can find the optimal stopping value for each value of h. In this test we fix

α = 1, r = τ = h √ 8 .
Optimal choice of the parameters r, τ for the primal-dual algorithm Starting from the term σ+rDu, we expect that r ∼ [σ] [Du] ∼ h [σ] [u] , where the square brackets represent Table 5.2: Computed error for different values of h and ε tol . Boldface numbers correspond to the optimal choice of ε tol , the largest for which we obtain the floor value of the error.

the order of magnitude of a quantity. Hence the expected order of magnitude of r is

r ∼ h [σ] [u] .
(5.2.18)

In the iteration scheme one can expect that the term αu has the same order of magnitude as the penalty term u k+1 -u k τ . Thus τ is comparable to 1 α . On the other hand, from the equation

αu -div σ = f , one has α ∼ [div σ] [u] ∼ [σ] h[u]
. Hence an expected order of magnitude of τ is which the error is the smallest is also the choice for which the number of iterations is the smallest.

τ ∼ h[u] [σ] . ( 5
Considering that here [u] = 2, the previous considerations lead to the "optimal" choice

r = 1 8 h [σ] [u] , τ = h[u] [σ] . ( 5

.2.20)

If we replace [u] by [f ] α and considering that here [f ] α = 32, this leads rather to the formula

r = 2h α[σ] [f ] , τ = 1 16 h[f ] α[σ] . (5.2.21)
Primal-dual algorithm without mass lumping Consider now using the original primal-dual iteration method (4.2.52) without the mass lumping modification (4.2.64). In other words we use the FEM scalar product instead of the FVM scalar product. The main difference between both methods is that now one has to solve a linear system of equations. Its cost is known to be generically O(n 2 log n) for a general matrix, and O(n log n) in case of a "good" sparse matrix, with n the size of the matrix.

For dealing with the method without mass lumping we use the software FreeFEM++, that automatically generates a mesh. For it we choose the number of intervals n = 40, we take u init = σ init = 0, τ = r = h/6. We make the test for a scalar function (i.e. for the so called TV minimization problem), with various values of ε tol . The results are as follows.

• Without mass lumping ε = 10 In conclusion, the primal-dual algorithm without or with mass lumping give the same quality of approximation for the same number of iterations. However the mass-lumped method runs faster because there is no linear system to solve. The gain is not big here because the mass matrix is almost diagonal and it is not much costly to invert it.

Accelerated scheme

We evaluate here the accelerated primal-dual algorithm (4.2.61) (with mass lumping). For this method we have to choose r 0 and τ 0 .

Optimal choice for r 0 , τ 0 Similarly as in the previous non-accelerated scheme we try several values. The error and number of iterations are reported in Tables 5.6, 5.7, 5.8. We take here α = 20. We can see that the most efficient choice is (τ 0 , r 0 ) = ( h 2 , h 4 ).

Comparison between the two methods

We choose α = 100 and consider the normal and accelerated primal-dual algorithms. Following Tables 5.9, 5.10, 5.11, 5.12, for any value of h the accelerated scheme with the best parameters (τ 0 , r 0 ) = ( h 2 , h 4 ) gives only a slightly smaller error than the normal scheme with the best parameters (τ, r) = ( h 16 , 2h). Meanwhile, the accelerated scheme uses more iterations to attain this value of the (τ 0 , r 0 ) ) for both methods, the accelerated scheme performs better. We conclude that finding optimal values for r, τ (or r 0 , τ 0 ) is more important that choosing the normal or accelerated scheme. For the normal scheme the formula (5.2.21) performs quite well since [f ]/α = 2 here.

h 2 8 ,

Implementation of various types of boundary conditions

In the previous tests we always implemented the Dirichlet boundary condition when solving the viscoplastic part (1.3.2). It is however possible to formulate and implement the finite element formulation in the case of other boundary conditions, by adding degrees of freedom on the boundary of the domain. One can treat in particular Neumann or slip boundary conditions. The implementation of a friction condition is a bit more difficult.

We consider the viscoplastic model

αu -div σ = f, (5.2.22) (τ, r) h 2 8 , 1 h 16 , 2h h 8 , h h 4 , h 2 h √ 8 , h √ 8 h 2 , h 4 h, h 8 2h, h 16

Comparison with the regularization and augmented Lagrangian methods

We implement the primal-dual algorithm, its accelerated version, the regularization method, and the augmented Lagrangian method within the same FreeFEM++ software, so as to compare them in particular concerning their execution time. The test case is the radial one of Subsection 5.2.1, with α = 100. The results are shown in Tables 5.13, 5.14, 5.15, 5.16. We observe that for a given space resolution, the error that we obtain with each method (taking its best parameters) is more or less the same, the accelerated primal-dual method having a slightly smaller error. On the contrary the number of iterations varies a lot. More iterations are necessary with the primaldual algorithm, and in particular with the accelerated algorithm (which is unexpected since it is supposed to converge faster). We observe that comparing the unstructured mesh used here with the Cartesian mesh used with the Fortran code, the primal-dual algorithm takes more or less the same number of iterations, whereas the accelerated version takes more iterations in the unstructured mesh configuration than in the Cartesian one (compare the last columns of Tables 5.13, 5.14 to Table 5.11). Then we have to compare execution time, which is the real issue since each step of the iteration does not have the same cost for all the methods (regularisation and augmented Lagrangian methods have linear systems to solve). We can see that the primal-dual algorithm takes more time

Time dependent case

We now consider the Bingham problem in the time dependent case, i.e. (5.2.1) with αu replaced by ∂ t u. We build an exact solution as previously with u(t, x, y) = Φ(t, r) -y x with r = x 2 + y 2 .

Then the problem can be written

∂ t Φ(t, r) - ∂ r sgn ∂ r Φ(t, r) r - 2sgn ∂ r Φ(t, r) r 2 -y x = f. (5.2.38)
An exact solution is built under the form

• For 0 ≤ r ≤ 1/6, sgn ∂ r Φ = (12r -36r 2 ) 2 , Φ = t, f = 1 -2 × 12 2 (1 -3r)(2 -9r) -y x . • For 1/6 ≤ r ≤ 1/3, sgn ∂ r Φ = 1, Φ = 6rt, f = 6r -2/r 2 -y x . • For 1/3 ≤ r ≤ 1/2, sgn ∂ r Φ = cos(π(6r -2)), Φ = 2t, f = 2 + 6π sin(π(6r -2)) r - 2 cos(π(6r -2)) r 2 -y x .
• For 1/2 ≤ r ≤ 5/6, sgn ∂ r Φ = -1, Φ = (5 -6r)t, f = 5 -6r + 2/r 2 -y x .

• For 5/6 ≤ r ≤ 1,

sgn ∂ r Φ = - 1 + cos(π(6r -5)) 2 , Φ = 0, f = -3π sin(π(6r -5)) r + 1 + cos(π(6r -5)) r 2 -y x .
• For 1 ≤ r, sgn ∂ r Φ = 0, Φ = 0, f = 0.

Comparison between various values of the timestep

Considering that the problem is of parabolic type we have a natural value of the timestep

∆t ∼ h 2 [Du] [σ] . ( 5
F ij = -F ji . (5.3.6)
This ensures that i |Q i |U n i is time independent. We take the numerical flux of the form

F ij = F (U i , U j , n ij ), (5.3.7)
which means that given the direction n ij we have to solve the system (5.3.3) in the direction n ij . The consistency of the numerical scheme is that F (U, U, n) = n x F x (U ) + n y F y (U ). More details can be found in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF].

Formulation with dual meshes

The location of the degrees of freedom are the (x i , y j ) as defined by (5.2.3)-(5.2.7). The primal mesh is the FEM mesh made of triangles with (x i , y j ) as nodes. The dual mesh is the FV mesh made of rectangles around each (x i , y j ).

Introduction

The lubrication equation is quite a classical simplification of the incompressible Navier-Stokes system. It is obtained for thin films of fluid, when viscous effects balance the pressure force. This occurs for instance for thin films of oil, hence the name of the equation. The study of this approximation goes back to Reynolds in 1886 [START_REF] Reynolds | On the theory of lubrication and its application to Mr. Beauchamp Tower's experiments, including an experimental determination of the viscosity of olive oil[END_REF]. Several scalings are involved to obtain this model. First the aspect ratio between the thickness of the film and the characteristic length of the substrate must be small, say δ. Simultaneously, the time scale has to be of order 1/δ. This is the so-called long wave regime, and is classically used in the shallow-water approximation. The lubrication equation requires another assumption of balance between the viscous effects and the pressure effects, which amounts to neglect all kinematic effects. This simplified flow is known as the Stokes flow. The lubrication equation itself is then obtained by integration over the fluid thickness. We are interested here in the lubrication model for a class of non Newtonian fluids. Several fluids are known to depart from the usual Newtonian rheology, where the deviatoric stress tensor is a linear function of the strain rate tensor, thus defining the dynamical viscosity of the fluid. The lubrication equation for Newtonian fluids has been studied for instance by Huppert [START_REF] Herbert | The propagation of two-dimensional and axisymmetric viscous gravity currents over a rigid horizontal surface[END_REF]. Non Newtonian fluids arise in several applications in engineering, biology, geophysics... In particular, viscoplastic or pseudoplastic fluids are involved in various geological problems, for instance lava flows, mudslides and avalanches. We refer to [START_REF] Ancey | Plasticity and geophysical flows: a review[END_REF] for a review on the subject. A model which is widely used is the so-called Bingham-plastic model. This model involves a yield stress, namely a threshold on strain rate: for values of the strain rate above this threshold the fluid behaves like a viscous fluid, for values below, it looks like a solid. This can be thought of as an infinite viscosity fluid. We refer to the papers by Liu and Mei [START_REF] Ko | Slow spreading of a sheet of Bingham fluid on an inclined plane[END_REF] and Balmforth et al. [START_REF] Neil | Viscoplastic flow over an inclined surface[END_REF] for the study of such fluids in the lubrication approximation. Both papers contain also a complete bibliography. Liu and Mei also introduced in [START_REF] Liu | Approximate equations for the slow spreading of a thin sheet of Bingham plastic fluid[END_REF] a perturbed Bingham model, which is actually a two viscosities model, with a high viscosity for small deformations. When this viscosity goes to ∞ the Bingham model is recovered, thus giving a fluid mechanics interpretation of this solid behaviour. This is precisely the two viscosities model we investigate here. First we describe the mathematical model we use, namely the incompressible Navier-Sokes equations in a time-dependent domain, since we consider a free-boundary problem. In particular we explain in some details all the scalings involved. Next, we turn to the lubrication equation itself, which is a one-dimensional equation, obtained by averaging the previous ones along the thickness. Finally we provide a few numerical illustrations based on a finite volume scheme.

Mathematical model

In this section we set up the model. The starting point is the incompressible Navier-Stokes system. We limit ourselves in this paper to the two-dimensional case, thus aiming at a one-dimensional lubrication equation. Similar computations can be performed in three space dimensions. The domain we consider is Ω t defined by f b (x) < z < ϕ(t, x), for t > 0 and x ∈ (-∞, +∞), where f b is given topography, and ϕ is a free surface. The notation we use is gathered in Figure 6.1.

z * ϕ(t, x) f b (x) substrate x z 0 h * h * h(t, x)
Figure 6.1: Notation for the two viscosities fluid: ϕ is the free surface; f b is the topography of the substrate; z * is the ordinate which separates "small deformations" (white zone) from "large deformations" (green zone), see Section 6.3 below. We introduce the thicknesses h

= ϕ -f b , h * = ϕ -z * , h * = h -h * .
The incompressible Navier-Stokes equations are

∂ x u + ∂ z v = 0, (6.2.1)
∂ t u + u∂ x u + v∂ z u = - 1 ρ ∂ x p + ∂ x τ xx + ∂ z τ xz , (6.2.2)
∂ t v + u∂ x v + v∂ z v = -g - 1 ρ ∂ z p + ∂ x τ zx + ∂ z τ zz , (6.2.3)
where ρ is the density of the fluid, U = (u, v) is the velocity field, and the stress tensor σ is written as the sum of a volumetric stress tensor, involving the pressure p, and a deviatoric stress tensor τ :

σ = -pId + τ, τ = τ xx τ xz τ zx τ zz .
The density ρ is assumed to be constant here, and the tensor σ will be defined in Section 6.2.1 below. fluid, that allows to handle explicit computations. The main feature of this kind of fluids is a nonlinear viscosity, decreasing with the strain rate. Mimicking the Bingham model, which is based on a threshold on the shear stress, we consider a model with a threshold on the strain rate: the viscosity is equal to some large µ B for small deformations, that is γ < γ c , where γ c > 0 is a given constant, and to another value µ for large deformations, γ > γ c . Such models were introduced by Liu and Mei [START_REF] Liu | Approximate equations for the slow spreading of a thin sheet of Bingham plastic fluid[END_REF], and the limit case ν B → ∞, which leads to a Bingham fluid, is studied in [START_REF] Ko | Slow spreading of a sheet of Bingham fluid on an inclined plane[END_REF] and [START_REF] Neil | Viscoplastic flow over an inclined surface[END_REF]. Notice that using (6.2.8) the threshold γ c on γ can be replaced by a threshold γ c = γ c / √ 2 on ε .

A multidimensional formulation for these simplified pseudo-plastic fluids is therefore

τ P P =    2ρν B ε if ε γ c , 2ρν ε + 2ρ(ν B -ν)γ c ε ε if ε > γ c , (6.2.9)
where we have introduced the kinematic viscosities ν and ν B . A particular limit case is ν B → ∞, which leads to a Bingham type fluid. To view this, it is convenient to define the following quantities (see Figure 6.3 below for an illustration in 1 dimension)

τ c = ν B γ c , τ * = (ν B -ν)γ c = (1 -ν/ν B )τ c , (6.2.10)
so that definition (6.2.9) can be rewritten

τ P P =    2ρν B ε if ε τ c /ν B , 2ρν ε + 2ρ(1 -ν/ν B )τ c ε ε if ε > τ c /ν B . (6.2.11)
It is clear on this formulation that the relevant limit is ν B → +∞, together with γ c → 0, keeping ν B γ c = τ c . In doing so, we recover the classical Bingham stress tensor, with threshold τ c :

τ Bing =    any τ s.t. τ τ c if ε = 0, 2ρν ε + 2ρτ c ε ε if ε > 0.
Finally, notice that in the pseudo-plastic (or shear thinning) context, we consider 0 < ν < ν B , but similar computations can be performed in any case. It is convenient for the scalings below to rewrite expression (6.2.11) using an equivalent kinematic viscosity ν eq , which satisfies ν ν eq ν B :

τ P P = ρν eq ε,

where

ν eq =    2ν B if ε γ c , 2ν + 2(1 -ν/ν B ) τ c ε if ε > γ c . ( 6
.2.12)

Scalings

We introduce now the scaling laws, namely thin layer, or more precisely long wave approximation, and slow motion, in order to finally obtain the lubrication model. This kind of scalings is already present e.g. in [START_REF] Neil | Viscoplastic flow over an inclined surface[END_REF] in the context of a visco-plastic fluid. Hence we propose the following family of scalings: we introduce a first set of characteristic scales, namely dimensions 0 and h 0 , characteristic velocities u 0 and v 0 , and a characteristic time t 0 . The quantities 0 and u 0 correspond to the horizontal direction, h 0 and v 0 to the vertical one. The aspect ratio δ = h 0 / 0 will be an important parameter, assumed to be small in the thin layer case. Dimensionless variables are then defined by

x = 0 x, z = h 0 z, t = t 0 t u = u 0 ū, v = v 0 v.
First, we rewrite the incompressibility equation (6.2.1) in the rescaled variables. We obtain u 0 0

∂ x ū + v 0 h 0 ∂ z v = 0,
and following the least degeneracy principle [START_REF] Van Dyke | Perturbation methods in fluid mechanics/annotated edition[END_REF], this implies u 0 / 0 = v 0 /h 0 , or equivalently 0 /h 0 = u 0 /v 0 . Thus v 0 /u 0 = δ, so that in the thin layer approximation v 0 is also small compared to u 0 .

We turn now to the kinematic part of the equation. Using u 0 / 0 = v 0 /h 0 , we readily obtain

∂ t u + u∂ x u + v∂ z u = u 0 t 0 ∂t ū + u 0 v 0 h 0 ū∂ x ū + u 0 v 0 h 0 v∂ z ū.
Once again we apply the least degeneracy principle and obtain t 0 = 0 /u 0 = h 0 /v 0 , or, as expected, u 0 = 0 /t 0 and v 0 = h 0 /t 0 . We proceed in the same way for the momentum equation in v and finally obtain

∂ t u + u∂ x u + v∂ z u = u 0 v 0 h 0 (∂t ū + ū∂ x ū + v∂ z ū) = δ u 2 0 h 0 (∂t ū + ū∂ x ū + v∂ z ū) , (6.2.13)
∂ t v + u∂ x v + v∂ z v = u 0 v 0 0 (∂tv + ū∂ x v + v∂ z v) = δ 2 u 2 0 h 0 (∂tv + ū∂ x v + v∂ z v) , (6.2.14)
where we have emphasized the aspect factor δ = h 0 / 0 = v 0 /u 0 . Following Balmforth [START_REF] Neil | Viscoplastic flow over an inclined surface[END_REF], we rescale the pressure and the stress tensor by

p = ρgh 0 p, τ = ρν u 0 h 0 τ . ( 6
δ 2 u 2 0 h 0 (∂tv + ū∂ x v + v∂ z v) = -g∂ z p -g + ν u 0 h 2 0 (δ∂ x τxz + ∂ z τzz ) . ( 6

.2.17)

Hence the condition ε > γ c leads us to define a dimensionless threshold γc = ( √ 2h 0 /u 0 )γ c = (h 0 /u 0 )γ c , so that the condition ε > γ c becomes |∂ z ū| > γc . Thus we get νeq

= ν eq ν =      2 ν B ν if |∂ z ū| γc , 2 + 2(1 -ν/ν B ) τ c ν h 0 u 0 √ 2 |∂ z ū| if |∂ z ū| > γc . (6.2.29)
We introduce a dimensionless viscosity νB and a dimensionless yield stress B by setting

νB = ν B ν 1, B = √ 2τ c h 0 νu 0 , (6.2.30)
so that the dimensionless deviatoric stress tensor becomes (see Figure 6.3)

τxz = νB ∂ z ū if |∂ z ū| γc , ∂ z ū + (1 -1/ν B )B sgn(∂ z ū) if |∂ z ū| > γc . (6.2.31)
This is the model proposed by Liu and Mei in [START_REF] Liu | Approximate equations for the slow spreading of a thin sheet of Bingham plastic fluid[END_REF]. The blue curve is the Bingham limit: ν B → +∞, γ c → 0 with γ c ν B = τ c . The green dashed line is the pure Newtonian limit ν B → ν.

As concerns the boundary conditions, we notice that the no-slip and kinematic boundary conditions remain unchanged by the scaling. In contrast, the continuity of the stress tensor across the free surface ϕ is greatly simplified. Recalling that ϕ = h 0 φ, we indeed obtain

σ • n = -p∂ x η + τ xx ∂ x ϕ -τ xz -p + τ xz ∂ x ϕ -τ zz = ρ   -δ pgh 0 ∂ x φ + δν u 0 h 0 τxx ∂ x φ -ν u 0 h 0 τxz -pgh 0 + δν u 0 h 0 τxz ∂ x φ -ν u 0 h 0 τzz   .
Now making use of (6.2.26), we obtain

σ • n = ρ     δ (δν eq u 0 h 0 ∂ x ū -gh 0 p)∂ x φ -δν eq u 0 2h 0 ∂ z v -ν eq u 0 2h 0 ∂ z ū δν eq u 0 2h 0 δ 2 ∂ x v + (∂ z ū ∂ x φ -∂ z v -gh 0 p     .
Letting δ go to zero gives therefore

[p]| φ = 0, [∂ z ū]| φ = 0, (6.2.32)

where [g] φ is the jump of some function g across φ. In other words, we recover separately the continuity of the pressure and the continuity of ∂ z u.

Lubrication equation

The so-called lubrication equation is obtained by integrating equations (6.2.1) along the vertical direction. The long wave and slow motion assumptions imply that we obtain a single nonlinear equation on the depth ϕ. Similar computations were performed by Liu and Mei [START_REF] Liu | Approximate equations for the slow spreading of a thin sheet of Bingham plastic fluid[END_REF], for a twoviscosity model, in order to justify the Bingham case, which corresponds to ν B → ∞ in our context. For Bingham fluids, we refer to Liu and Mei [START_REF] Ko | Slow spreading of a sheet of Bingham fluid on an inclined plane[END_REF], and more recently to Balmforth [START_REF] Neil | Viscoplastic flow over an inclined surface[END_REF]. The final equation is obtained through three steps we present in detail now.

We recall the equations we obtained in the preceding section, dropping the bars for clarity. First we have the hydrostatic relation

∂ z p = -1, f b z ϕ. (6.3.1)
Next, the dimensionless Stokes equation (6.2.23)

∂ z (τ xz ) = -∂ x p, f b z ϕ, (6.3.2)
where τ xz is the dimensionless deviatoric stress tensor defined by (6.2.31). These equations are coupled with the following boundary conditions (in these relations, t and x are hidden parameters):

• on the free surface z = ϕ p(ϕ) = 0, ∂ z u(ϕ) = 0, (6.3.3)

• on z = f b u(f b ) = u b , v(f b ) = v b . (6.3.4)
Concerning first the pressure, using the boundary condition on the free surface we obtain the usual hydrostatic approximation p(z) = ϕ -z, f b z ϕ. (6.3.5)

The averaged equation we look for is obtained by integrating in z the incompressibility equation, or mass conservation,

∂ x u + ∂ z v = 0.
This is quite classical, see e.g. [START_REF] Neil | Viscoplastic flow over an inclined surface[END_REF] in the same slow motion context, or [START_REF] James | Towards a new friction model for shallow water equations through an interactive viscous layer[END_REF] for shallow water approximation. We obtain v(t, x, ϕ) = v(t, x, f b ) - The kinematic boundary condition on z = ϕ leads to v(t, x, ϕ) -u(t, x, ϕ)∂ x ϕ = ∂ t ϕ = ∂ t h in the first equation. For z = f b , we make use of the no-slip boundary condition (6.3.4), to obtain the following averaged equation The first step towards the computation of the flux is to obtain the vertical velocity profile. The general structure of this profile is as follows. We have τ xz = F (∂ z u), where F is a continuous, oneto-one, increasing function, with F (0) = 0, see (6.2.31) and Figure 6.3. From (6.3.2) and (6.3.5) we are led to solve ∂ z F (∂ z u) = ∂ x ϕ. Since F (∂ z u) = 0 for z = ϕ (or equivalently ∂ z u = 0) we get F (∂ z u) = ∂ x ϕ(z -ϕ), so that F (∂ z u) is monotone (increasing if ∂ x ϕ 0, decreasing if not). Because F is increasing, ∂ z u is monotone as well, in particular, since ∂ z u = 0 for z = ϕ, its sign remains constant. Therefore |∂ z u| is decreasing in z (increasing with depth).

∂ t h + ∂ x
The threshold in formula (6.2.31) eventually splits the fluid in two layers. Let z * be defined by |∂ z u(z * )| = γ c . Provided z * ∈]f b , ϕ[ (see below for precise formulas), we have a "small deformation", that is |∂ z u(z)| < γ c , region for z ∈]z * , ϕ[, because |∂ z u| is decreasing from 0 for increasing depth. Similarly for z ∈]f b , z * [ we have |∂ z u(z)| > γ c , so that finally, according to (6.2.31), the velocity is ruled by the system of equations Notice that the curve z = z * is not a physical interface, yet we have continuity of the stress tensor, or equivalently here continuity of ∂ z u. Now the computations are quite easy. We integrate once the first equation between ϕ and z * , to obtain These definitions ensure that z * f b and h * h, and are valid for ∂ x ϕ = 0 with the convention B/0 = ∞. Notice that z * can be equal to f b for weak slopes (small ∂ x ϕ), or small depths (small h). Conversely, z * → ϕ when |∂ x ϕ| goes to ∞.

∂ z u = 1 ν B ∂ x ϕ(z -ϕ).
Integrating once again between z * and ϕ, we obtain the velocity profile for ϕ z z * :

u(z) = 1 2ν B ∂ x ϕ(ϕ -z) 2 + K,
where the constant K will be determined later. Notice for further use that by construction (6.3.10) We turn now to the lower layer, z * z f b . The fluid here has dimensionless viscosity 1, and we use the boundary conditions (6.3.10) for z = z * , and no slip (6.3.4) at z = f b . First we get, using (6.3.10),

∂ z u(z * ) = 1 ν B ∂ x ϕ(z * -ϕ) = γ c .
∂ z u = ∂ x ϕ(z -z * ) - 1 ν B ∂ x ϕh * ,
next, integrating once again between f b and z * ,

u = 1 2 ∂ x ϕ(z * -z) 2 - 1 ν B ∂ x ϕh * z + L,
where L is computed using (6.3.4), leading to

L = u b - 1 2 ∂ x ϕ(z * -f b ) 2 + 1 ν B ∂ x ϕh * f b , so that u = 1 2 ∂ x ϕ (z * -z) 2 -(z * -f b ) 2 - 1 ν B ∂ x ϕh * (z -f b ) + u b , (6.3.11)
Finally, we use the continuity of the velocity at z = z * to obtain the constant K:

1 2ν B ∂ x ϕ(ϕ -z * ) 2 + K = - 1 2 ∂ x ϕ(z * -f b ) 2 - 1 ν B ∂ x ϕh * (z * -f b ) + u b .
The velocity profile is therefore given by

u(z) =            ∂ x ϕ 2 (z * -z) 2 -(z * -f b ) 2 - ∂ x ϕ ν B h * (z -f b ) + u b , f b z z * , ∂ x ϕ 2 (ϕ -z) 2 -(ϕ -z * ) 2 - ∂ x ϕ 2 (z * -f b ) 2 - ∂ x ϕ ν B h * (z * -f b ) + u b , z * z ϕ.
(6.3.12) Notice that for ν B = 1, easy computations show that the profile is the same in the two layers, namely u = ∂ x ϕ 2 (z -f b )(z -f b -2h) + u b , which is as expected the usual parabolic profile for a Newtonian fluid.

On the other hand, letting ν B → +∞, and γ c → 0, keeping ν B γ c = τ c , we recover formally the Bingham fluid velocity, as in Balmforth [START_REF] Neil | Viscoplastic flow over an inclined surface[END_REF]: Therefore the flux we are looking for is given by It is easy once again to check on this formula that we recover the usual cubic flux -∂ x ϕ 3 h 3 for the Newtonian fluid ν B = 1. On the other hand the limit case ν B → ∞ gives back Balmforth's formula Inserting (6.3.13) in the conservation equation (6.3.8) leads to the following advection-diffusion equation: Notice that 0 < D(h, ∂ x h) h 3 /(3ν B ).

u(z) =          - ∂ x ϕ 2
∂ t h + ∂ x (u b h) = v b + ∂ x (D(h, ∂ x h)∂ x (h + f b )) , (6.3

Numerical illustrations

We turn now to numerical examples to illustrate the behaviour of the two-viscosity fluid. The point here is not to give an accurate specific scheme, which is an interesting perspective since the diffusion term may degenerate, but is beyond the scope of this work. We merely apply here a simple finite volume strategy. The infinite space domain is replaced by some finite computational domain [a, b].

Since we do not want to cope with boundary conditions here, we merely impose a free flux on the boundaries, which is compatible with the examples we choose. Positive time and space steps ∆t and ∆x being given, we introduce the usual notation t n = ∆t, n 0, and x j = j∆x, 0 j J, where J = (b -a)/∆x. An approximation of the depth h is saught for in the form

h n+1 j = h n k - ∆t ∆x (F n j+1/2 -F n j-1/2 ) + ∆t ∆x (G n j+1/2 -G n j-1/2 ),
where F n j+1/2 is the numerical advection flux, and G n j+1/2 the numerical diffusion flux, both computed at interface x j+1/2 . In the following we denote u n j the discretized basal velocity, and f j the discrete topography, which are both given functions.

Figure 1 . 1 :

 11 Figure 1.1: Time independent stress -strain rate relation for different types of fluids

p = - 1 N

 1 Tr(σ). (1.2.5) Combining (1.2.4) and (1.2.5) we get σ = -p Id +σ , Tr(σ ) = 0. (1.2.6)

2 . 2 N+Figure 1 . 3 :

 2213 Figure 1.3: A domain with several boundary parts

.1. 3 )

 3 and when this holds one has F (D) < ∞ and F * (σ) < ∞. (c) σ ∈ ∂F (D) if and only if the sup in the definition of F * * (D) = F (D) is attained at σ. (d) For D, σ ∈ H, one has the Fenchel-Young inequality [49] σ : D ≤ F (D) + F * (σ).(2.1.4)Moreover, equality holds if and only if σ ∈ ∂F (D).

  Proof. (a) If ∂F (D) is nonempty, there exists some σ ∈ ∂F (D). Then by definition of the subdifferential, F (D) ≥ F (D) + σ : (D -D), ∀D ∈ H. (2.1.6) Since F is proper there is at least one D such that F (D) is finite. Hence the previous inequality yields that F (D) < ∞. (b) Assume that σ ∈ ∂F (D). Then the inequality (2.1.6) holds. According to the definition of F * , we have F * (σ) ≥ D : σ -F (D), ∀σ ∈ H.

  .1.10) Taking D = D 2 in the first inequality and D = D 1 in the second one, we get

.1. 11 )Proposition 2 . 1 . 4 (

 11214 Since ∂F (D 1 ) and ∂F (D 2 ) are nonempty, we haveF (D 1 ) < ∞, F (D 2 ) < ∞.Comparing the two inequalities of (2.1.11) we conclude that 0 ≥ σ 1 : (D 2 -D 1 ) + σ 2 : (D 1 -D 2 ), which is the claim. (f) The definition (2.1.1) of F * gives F * (0) = sup(-F ) = -inf F , and reversing the role of F and F * we get the result. Moreau envelope and proximal operator). Let F : H → R a convex, proper, l.s.c. function. Then for any ε > 0 and D ∈ H we can define the Moreau envelope of F as F ε (D) := inf D∈H F (D) + |D -D| 2 2ε .

Remark:

  Any convex, proper and lower semi-continuous F is lower bounded by an affine function. It follows from Proposition 2.1.4(a): by taking arbitrary D ∈ H and ε > 0 we get some D ∈ H satisfying (2.1.14), which implies that ∂F ( D) = ∅. Thus there is some σ ∈ ∂F ( D), and the definition (2.1.2) of ∂F ( D) ensures that F is lower bounded by the affine function D → F ( D) + σ : (D -D).

.4. 8 ) 8 )

 88 For the case σ ∈ ∂F (0), we have Id + ∂F r we obtain P r (σ) = σ. For σ / ∈ ∂F (0) there exists a unique D such that ∂F (D) + rD σ. Then since σ / ∈ ∂F (0), D cannot be zero, and thus F (D) + rD = σ.

  and by Lemma 2.1.3(a) it follows that ∂F * (σ) = ∅. (c) If σ ∈ ∂A, a vector D ∈ H verifies D ∈ ∂F * (σ) if and only if

4 )

 4 with D the deviator of D, i.e. D = D -Tr(D) Id /N . Then since A = ∂F (0), the stresses σ ∈ A are characterized by σ : D ≤ |D | ∀D, It implies that A = {σ = -p Id +σ | p = 0, |σ | ≤ 1}. Next we compute the projection P r (σ) from (2.4.1), (2.4.8), P r (σ) = σ -r Id + ∂F r -1 σ r .

Figure 2 . 1 :

 21 Figure 2.1: Set A and orthogonal projection

.6. 10 )

 10 Besides, if Tr(σ -rD) = 0, |σ -rD| ≤ 1 with |D | = 0, then D = Tr(σ) rN Id, and |σ | ≤ 1.

Figure 2 . 2 :

 22 Figure 2.2: Set A in the Cam-Clay model, and orthogonal projection.

.6. 22 )

 22 with Tr(σ b ) = 0, |σ b | = 1, and 0 ≤ p b ≤ p 0 . Assume that 0 < p b < p 0 . Then for a small variation δσ b of σ b one has δσ b = -Id +σ b sin δ p 0 -2p b 2 p b (p 0 -p b ) δp b + δσ b sin δ p b (p 0 -p b ). (2.6.23) Then since Tr(δσ b ) = 0, σ b : δσ b = 0, we compute δσ b : Id = -N δp b , δσ b : σ b = sin δ p 0 -2p b Thus we get δσ b : σ b + sin δ N p 0 -2p b 2 p b (p 0 -p b ) Id = 0, (2.6.25)

.6. 31 )Figure 2 . 3 :

 3123 Figure 2.3: Set A in the Drucker-dilatant model with 4 possible cases of projections.

  Finding the orthogonal projection on the triangle We take σ = -p Id +σ ∈ A, and we have to find P A (σ) = σ b ∈ ∂A, such that σ = σ b + λn b , with n b the external normal to ∂A at σ b , and λ ≥ 0. Let us first consider the case of the inclined boundary. Then σ b = -p b Id +p b sin ψσ b , (2.6.35) with Tr(σ b ) = 0, |σ b | = 1. Then δσ b = -Id δp b + sin ψσ b δp b + p b sin ψδσ b . We have then Tr(δσ b ) = -N δp b , δσ b : σ b = sin ψδp b . Thus δσ b : (σ b + Id sin ψ/N ) = 0. Therefore we can take n b = σ b + Id sin ψ/N . We write that σ = σ b + λn b , which gives -p = -p b + λ sin ψ N , σ = (λ + p b sin ψ)σ b .

  and p sin ψ ≤ |σ |. (2.6.40) Next let us consider the case of the right boundary. In this case σ b = -p 0 Id +µσ b , with Tr(σ b ) = 0, |σ b | = 1, 0 ≤ µ ≤ p 0 sin ψ. Then δσ b = δµσ b + µδσ b , Tr(δσ b ) = 0. It follows that we can take n b = -Id. Then writing that σ = -p Id +σ = σ b + λn b we obtain -p = -p 0 -λ, σ = µσ b .

Figure 2 . 4 .Figure 2 . 4 : 0 = λ sup 0≤p≤ep 0 p

 242400 Figure 2.4: Set A in the logarithmic model, and orthogonal projection.

( 2 .

 2 6.48) Moreover, using Lemma 2.1.3 we obtain that for D = 0, ∂F (D) is a single point σ obtained by σ = λp(1 -ln p p 0 ) D |D | and -Tr(D) λ|D | if Tr(D) ≥ -λ|D |, ep 0 if Tr(D) ≤ -λ|D |. (2.6.49) In the case D = 0 and Tr(D) = 0, ∂F (D) is a single point σ defined as σ = 0 if Tr(D) > 0, σ = -ep 0 Id if Tr(D) < 0. We observe in particular that for D = 0 and Tr(D) ≥ -λ|D |, we have σ = p λ + Tr(D) |D | D |D | , (2.6.50) which is a friction law of the type expected in granular materials: proportional to the pressure with a coefficient which is the sum of a constant and the dilatancy ratio Tr(D)/|D |. The dilatancy law is Tr(D) = -λ|D | ln p p 0 , (2.6.51) with a critical pressure p 0 . Finding the orthogonal projection on A for the logarithmic model We take σ = -p Id +σ ∈ A, and we have to find P A (σ) = σ b ∈ ∂A, such that σ = σ b + γn b , with n b the external normal to ∂A at σ b and γ ≥ 0. One has σ b = -p b Id +λp b 1 -ln p b p 0 σ b , (2.6.52) with Tr(σ b ) = 0, |σ b | = 1. Then δσ b = -Id δp b -λ ln p b p 0 σ b δp b + λp b 1 -ln p b p 0 δσ b . (2.6.53) Then since Tr(δσ b ) = 0, σ b : δσ b = 0, we compute δσ b : Id = -N δp b , δσ b : σ b = -λ ln p b p 0 δp b . (2.6.54) Thus we get δσ b : σ b -λ N ln p b p 0 Id = 0, and it follows that n b = σ b -λ N ln p b p 0 Id. Writing that σ = σ b + γn b we obtain

Figure 2 .

 2 4 that for σ = -p Id +σ ∈ A, the projection P A (σ) = σ b is given by • If |σ | = 0 and p ≤ 0 then σ b = 0, • If |σ | ≤ N λ (p -ep 0 ) then σ b = -ep 0 Id, • If |σ | > 0 and |σ | > N λ (p -ep 0 ) then σ b is given by (2.6.52) with σ b = σ |σ | and p b ∈ (0, ep 0 ) solution to (2.6.58) such that p b 1 -ln p b p 0 ≤ |σ | λ .

NN - 1 p

 1 (it is also necessary if N = 2). It indicates that the set A should be cut below this line, similarly as in the Drucker-dilatant model. In this case, and if 0 ∈ A, we see geometrically that |D | < N -1 N Tr(D) ⇒ ∂F (D) = {0}. A related property is that | Tr(D)| < N N -1 |D | implies that the eigenvalues of D are not all of the same sign.

  .1.6) With Lemma 3.1.1 we obtain (3.1.4). (b) The inequality (3.1.4) implies that |F ε (D) -F ε (0)| ≤ C |D|(1 + |D|). (3.1.7) Following (3.1.3), F ε is lower bounded by an affine function independent of ε ≤ ε 0 . Thus F ε (0) is bounded below by a constant independent of ε. Moreover according to Proposition 2.1.4(f) we have F ε (0) ↑ F (0) as ε → 0. If F ε (0) were unbounded, (3.1.7) would imply that F ε (D) → ∞ for all D, which is not the case since F is proper. Thus F (0) < ∞, |F ε (0)| remains bounded, and (3.1.7) gives |F ε (D)| ≤ C (1 + |D| 2 ). Since F ε (D) ↑ F (D) we deduce that F (D) < ∞ with the same bound. Lemma 3.1.3. Let F satisfy Hypothesis 1 and assume that F is finite everywhere and there exists a constant C > 0 such that for all D |F (D)| ≤ C (1 + |D| 2 ). (3.1.8) Then (a) For all σ one has F * (σ) ≥ 1 4C |σ| 2 -C . (3.1.9) (b) For any D and σ such that σ ∈ ∂F (D) one has |σ| ≤ 8C (1 + |D|), (3.1.10) so that F satisfies Hypothesis 2. Proof. (a) Using the definition of F * and (3.1.8), one has

  .1.12) which yields (3.1.9). (b) Consider D 1 , D 2 , and a number R > 0 such that D 1 = D 2 and |D 1

  .1.18) Consider now D 1 and σ given, such that σ ∈ ∂F (D 1 ). It follows that for all D 2 , F (D 2 ) ≥ F (D 1 ) + σ : (D 2 -D 1 ). (3.1.19) With (3.1.18) we deduce σ : (D 2 -D 1 ) ≤ |D 2 -D 1 |8C 1 + max(|D 1 |, |D 2 |) . (3.1.20) Taking D 2 in a neighbourhood of D 1 but such that D 2 -D 1 takes all the possible directions, we conclude that |σ| ≤ 8C (1 + |D 1 |), which concludes (3.1.10). Remark: Lemmas 3.1.2 and 3.1.3 prove that for a nonlinearity F satisfying Hypothesis 1, one has equivalence between Hypothesis 2 and (3.1.8).

  (b) ψ is lower semi-continuous. Suppose u n → u in H 1s 0 . Then one has Du n → Du in L 2 . Thus after extraction of a subsequence, Du n → Du almost everywhere. Since F is l.s.c, one has F (Du) ≤ limF (Du n ) almost everywhere. According to Lemma 3.1.2(b), F (Du n ) is bounded in L 1 . Taking the integral over Ω, we get Ω F (Du) ≤ Ω limF (Du n ). (3.1.23) Next, F (D) is lower bounded by an affine function ϕ(D). Appying Fatou's Lemma to

.1. 26 ) 2 ≥ |Du| 2 2 +

 2622 Proof. Since g(D) = |D| 2 /2 for all D ∈ M s N ×N (R) is a convex function and g (D) = D, we have |Dv| 2 Du : (Dv -Du). This gives directly that (3.1.26) implies (3.1.25).

  .1.30) where G is a bounded linear form on H 1s 0 (Ω), and η, ∆t are positive constants. Moreover u ε , Du ε are bounded in L 2 (Ω) uniformly in ε ≤ ε 0 . Proof. According to Proposition 2.1.4 and Lemma 3.1.2(a), F ε satisfies the same hypotheses 1, 2 as F . Thus applying Proposition 3.1.7 to F ε we get the existence and uniqueness of u ε . Then inserting v = 0 in (3.1.30) we get 1 ∆t

.1. 33 )FF

 33 Since σ ∈ ∂F (Du) a.e. we have F (Dv) ≥ F (Du) + σ : (Dv -Du) a.e. Thus Ω (Dv) ≥ Ω (Du) + σ, Dv -Du . (3.1.34) Using this in (3.1.33) it follows that ∀v ∈ H 1s 0 (Ω) 1 ∆t u, v -u + Ω F (Dv) -Ω F (Du) + η Du, Dv -Du ≥ G(v -u). (3.1.35) Hence (3.1.28) holds.

  Consequently we obtain lim Ω ϕσ ε : Du ε ≤ Ω ϕσ : Du. Hence, taking the lim as ε → 0 in (3.1.45), we get Ω

Lemma 3 . 2 . 1 .

 321 Under Hypothesis 1, 3, 4, we have (a) F ε is uniformly lower bounded, E + F ε (D) ≥ 0, (3.2.4) (b) For all ε > 0 and all D, F ε (D) ≤ C(1 + |D| 2 ), (3.2.5) (c) The estimates obtained in Lemma 3.1.3 are valid. (d) If hypothesis 5 is satisfied, then it is also satisfied for F ε , uniformly for ε ≤ ε 0 . Proof. (a) By (3.2.1) one has F (D) ≥ -E, and the definition (2.1.12) of F ε gives immediately F ε (D) ≥ -E. (b) One has F ε (D) ≤ F (D), so that (3.2.2) gives immediately (3.2.5). (c) Because of (3.2.1) and (3.2.2), F satisfies |F (D)| ≤ max(C, E)(1 + |D| 2 ). Therefore the assumptions of Lemma 3.1.3 are satisfied. (d) Hypothesis 5 means that ∀λ > 0, ∃R > 0 inf |D|≥R F (D) |D| ≥ λ. (3.2.6)

.2. 12 )

 12 It follows with (3.2.7) that F ε (D)/|D| ≥ λ/2 as soon as |D| ≥ R, for some R large enough and depending only on λ, ε 0 , E, R. This concludes that F ε (D)/|D| → ∞ as |D| → ∞, uniformly in ε as long as ε ≤ ε 0 .

Lemma 3 . 2 . 3 .

 323 When F satisfies Hypothesis 1, 3, 4, we have (a) K F (Ω) and K F 0 (Ω) are convex and contain 0. (b) Assuming additionally Hypothesis 5, if we have a sequence u n ∈ K F (Ω) such that u n L 2 and Ω F (Du n ) are bounded, then Du n L 1 is bounded and there exist u ∈ K F (Ω) and a subsequence u n such that u n u in L 2 (Ω) weak, Du n Du in L 1 (Ω) weak. (3.2.15)

FF

  .2.53) Since σ ∈ ∂F (Du) a.e. we have F (Dv) ≥ F (Du) + σ : (Dv -Du) a.e. Thus Ω (Du) ≥ G(v -u). (3.2.55) Hence (3.2.34) holds for all v ∈ K F 0 (Ω). But if v ∈ L 2 (Ω) and v ∈ K F 0 (Ω), then ψ(v) = ∞, and (3.2.34) holds also trivially. Finally (3.2.34) holds for all v ∈ L 2 (Ω).

Lemma 3 . 2 . 10 .

 3210 Let F be a convex, proper, l.s.c. function on a Hilbert space. If F is superlinear i.e.F (D) |D| → ∞ as |D| → ∞, then F * is finite everywhere. In finite dimension the converse is true. Proof. Suppose that F is superlinear. One hasF * (σ) = sup D (σ : D -F (D)). (3.2.71)For a given σ, since F is superlinear there exist M > 0 (depending on |σ|) such that|D| ≥ M =⇒ F σ : D -F (D) ≤ σ : D -|D||σ| ≤ 0. (3.2.73) But since F is lower bounded by an affine function, there is some β such that F (D) ≥ β for all D such that |D| ≤ M . It follows that ∀|D| ≤ M, σ : D -F (D) ≤ M |σ| -β. (3.2.74) Therefore with (3.2.73) we deduce that F * (σ) < ∞.

Figure 5 . 1 :

 51 Figure 5.1: Locations x i of the degrees of freedom and size of the cells

Figure 5 .Figure 5 . 2 :

 552 Figure 5.2: One-dimensional compressible Bingham model: comparison between exact and approximate solution computed by the regularisation method. Left: ρ, right: u.

  Error between u x h exact and u x h computed. It is everywhere less than 4%.

Figure 5 . 4 :

 54 Figure 5.4: Primal-dual algorithm with structured mesh for h = 1/80

- 4 .•

 4 With mass lumping ε = 10 -2 err = 0.770504 n.o.iter = 92 execution time = 12.4226s ε = 10 -3 err = 0.112977 n.o.iter = 419 execution time = 52.3369s ε = 10 -4 err = 0.111199 n.o.iter = 752 execution time = 93.9013s ε = 10 -5 err = 0.111853 n.o.iter = 1077 execution time = 134.989s

Figure 5 . 6 :

 56 Figure 5.6: Numerical results for slip boundary condition

Figure 6 . 3 :

 63 Figure 6.3: Simplified shear-thinning model. We consider a piecewise linear approximation (in black) of the "theoretical" pseudoplastic law (in red). Parameters τ * and τ c are defined by (6.2.10). The blue curve is the Bingham limit: ν B → +∞, γ c → 0 with γ c ν B = τ c . The green dashed line is the pure Newtonian limit ν B → ν.

∂

  x u(z) dz (6.3.6) = v(t, x, f b ) -∂ x ϕ f b u(z) dz + u(t, x, ϕ)∂ x ϕ -u(t, x, f b )∂ x f b .(6.3.7)

  dz can be computed explicitly as a function of ϕ, by integrating twice equation (6.3.2).

ν

  B ∂ zz u = ∂ x ϕ, z * z ϕ, ∂ zz u = ∂ x ϕ, f b z z * ,where for the second equation we have used that ∂ z u has a constant sign. These equations are complemented with the boundary conditions∂ z u = 0, z = ϕ ; u = 0, z = f b .

  This leads toγ = |∂ z u S | = 1 ν B |∂ x ϕ|(ϕ -z),so that the value of z * and the thickness h * of this layer are given byz * = max ϕ -B |∂ x ϕ| , f b , h * = ϕ -z * = min B |∂ x ϕ|, h . (6.3.9)

(y * ) 2 -- ∂ x ϕ 2 (y * ) 2 +

 222 (y * -(z -f b )) 2 + u b , f b z z * , u b , z * z ϕ,where we have sety * = z * -f b = h -h * .It is now straightforward to obtain the flux in (6.3.8), sinceϕ f b u(z) dz = z * f b u(z) dz + ϕ z * u(z) dz.We have on the one handz * f b u(z) dz = -∂ x ϕ 3 (y * ) 3 -∂ x ϕ 2ν B h * (y * ) 2 + u b y * , on the other hand ϕ z * u T (z) dz = -∂ x ϕ 3ν B (h * ) 3 -∂ x ϕ 2 (y * ) 2 h * -∂ x ϕ ν B y * (h * ) 2 + u b h * .

  y * ) 2 h * + 3 ν B y * (h * ) 2 + 1 ν B (h * ) 3 + u b h. (6.3.13)

  ) dz = -∂ x ϕ 6 (y * ) 2 (y * -3h).

Figure 6 . 4 :

 64 Figure 6.4: Comparison between the three models, time evolution. Yield stress B = 0.1, ν B = 100

Figure 6 . 5 :Figure 6 . 6 :

 6566 Figure 6.5: Comparison between the three models, time evolution. Yield stress B = 0.5, ν B = 100

Figure 6 . 7 :

 67 Figure 6.7: Convergence towards the Bingham model. Yield stress B = 1.25. From top to bottom: ν B = 10, 100, 1000.

Figure 6 . 8 :

 68 Figure 6.8: Lift of a magma bubble -Yield stress B = 1.25 -ν B = 100 -Top left: basal velocity -Top right: Timelapse oh thickness h and low viscosity layer y * -Next 6 pictures: time evolution of the three models

Figure 6 . 9 :

 69 Figure 6.9: Descent of a magma bubble -Yield stress B = 1.25 -ν B = 100 -Top left: basal velocity -Top right: Time-lapse of thickness h and low viscosity layer y * -Next 6 pictures: time evolution of the three models

  Imposing this condition, the strain rate D = Du is no longer an arbitrary symmetric matrix, but it has vanishing trace: Tr(D) = 0. Therefore in this case the constitutive relation (1.2.8) needs to be understood with the subdifferential holding in the vector space of symmetric trace free matrices. It follows that the pressure p is not determined by the constitutive relation. Instead it appears as a Lagrange multiplier for the constraint (1.2.16). Thus the incompressible formulation is slightly different from the compressible formulation described above. In this thesis we consider only the compressible formulation.

	1.2.1) reduces to
	div u = 0.	(1.2.16)

1.2.6 Problem statement and state of the art Combining (1.2.1), (1.2.2), (1.2.3), (1.2.8), (1.2.14), (1.2.15), our viscoplastic problem writes: Find (ρ, u) satisfying

  ∀σ ∈ H. When σ 1 ∈ ∂F (D 1 ) and σ 2 ∈ ∂F (D 2 ), we have

	Comparing to (2.1.2) this means exactly that D ∈ ∂F * (σ), or by (b) that σ ∈ ∂F (D).
	(d) According to the definition (2.1.1) of F * , for a given σ one has F * (σ) ≥ σ : D -F (D) for all D,
	and the inequality (2.1.4) follows. Equality means that D → σ : D -F (D) attains its maximum
	at D, which by (c) (reverting the role of F and F * ) is equivalent to D ∈ ∂F * (σ), which gives the
	claim with the property (b).
	(e)

  -Bulkley model or the Bingham model, with or without viscosity. Let us here consider the Bingham fluid with, which corresponds to F (D) = σ 0 |D| + η|D| 2 /2, for some yield stress σ 0 > 0 and viscosity η ≥ 0. This particular problem can be considered as a model for quasi-Newtonian viscous fluids [53, Chapter 2]. Then we have to find u ∈ H such that

	αu -η div Du -σ 0 div	Du |Du|	= f in Ω,	(2.2.1)

  .1.21) with the scalar product u, v L 2 + Du, Dv L 2 . Then H 1s (Ω) is a Hilbert space. Lemma 3.1.5. Assume that F satisfies Hypothesis 1, 2. Then ψ is finite everywhere, convex and lower semi-continuous on H 1s 0 (Ω).

		We shall denote
	by H 1s 0 (Ω) the closure of C ∞ c (Ω) in H 1s (Ω).
	Definition 3.1.4. Suppose that F satisfies Hypothesis 1, 2 and define for Ω an open bounded subset
	of R N , ψ : H 1s 0 (Ω) → R as
	∀u ∈ H 1s 0 (Ω),	ψ(u) =

Ω F (Du). (3.1.22) Remark: The shorthand notation Ω F (Du) is used, meaning Ω F (Du(x)) dx. According to Lemma 3.1.2(b) and since Ω is bounded, the integral is always finite.

  .1.41) According to Lemma 2.1.3(f), (F ε ) * is lower bounded by a constant independent of ε. Taking w = u ε in (3.1.38) we get that : σ ε ≤ G(u ε ).(3.1.42)The right-hand side is bounded because u ε and Du ε are bounded in L 2 (Ω), and with (3.1.41) we deduce that Ω (F ε ) * (σ ε ) is upper bounded. Since (F ε ) * is lower bounded by a constant we conclude that

Ω

Du ε

  Proof of Lemma 3.2.8. Let w ∈ K F 0 (Ω). Then according to Lemma 3.2.4, there exists a sequencew k ∈ C ∞ c (Ω) such that w k → w in L 2 , Dw k → Dw in L 1 ,and F (Dw k ) → F (Dw). Extracting if necessary a subsequence we have Dw k → Dw a.e. in Ω. Since w k ∈ C ∞ c one has

	Ω

.

2.38)

It follows that σ : Dw is upper bounded by an L 1 function, and thus that σ : Dw is well-defined, as finite or -∞. The equation (3.2.37) says in particular that this integral must be finite, i.e. σ : Dw ∈ L 1 (Ω).

  (σ) ∈ L 1 thus g k ∈ L 1 . Since F is continuous we have that g k → g a.e, with

		e. in Ω.	(3.2.40)
	By assumption F g = F (Dw) + F * (σ) -σ : Dw ≥ 0.	(3.2.41)
	By Fatou's Lemma we have
	g ≤ lim	g k .	(3.2.42)
	Ω	Ω
	Using (3.2.39) we have

*

  Assume first that u ∈ L 2 and σ ∈ L 2 satisfy (3.2.46), (3.2.47), (3.2.48), (3.2.49). Then (3.2.47) gives that div σ ∈ L 2 , and according to Lemma 3.2.8 one has ∀v

.2.49)

Proof.

  Dû, Dû h -Dû + ψ(û h ) -ψ(û) ≥ f, ûh -û + η Du n , Dû h -Dû . (4.2.21)

	Taking v = ûh in (4.2.12), we get
	1 ∆t û, ûh -û + η Adding the inequalities (4.2.20) and (4.2.21) yields
	û -ûh ∆t	2	+ η Dû -Dû h	2 ≤ inf v h ∈V h	R(v h ).	(4.2.22)
	We notice that M in (4.2.11) is twice the right-hand side of (4.2.17), thus
			ûh ∆t	2	+ η Dû h	2 ≤ M.	(4.2.23)
	Therefore we can estimate (4.2.19) as
							.2.20)

  .2.27) But similarly as (4.2.17), taking v = 0 in (4.2.12) yields that û satisfies the same bound (4.2.23) as ûh . Thus together with (4.2.26) we obtain, noticing that by (4.2.22) the infimum is nonnegative,

  .2.33) because again v → F (Dv) is continuous on H 1s (Ω). From (4.2.32), (4.2.33) we get that for h small enough we have a function v h ∈ V h such that

	v h -û H ≤ 2ε and	F (Dv h ) -F (Dû) ≤ 2ε.	(4.2.34)
	This proves that R û → 0 as h → 0. With the estimates (4.2.8), (4.2.9) we conclude the convergence

of ûh to û in H and the convergence of F (Dû h ). About the weak L 1 convergence of Dû h , if η > 0 this is obvious since we have convergence in L 2 . If η = 0 we apply Lemma 3.2.3 and get that after extraction of a subsequence, Dû h converges in weak L 1 . Since the limit is necessarily Dû, it proves that the convergence holds without extracting any subsequence.

  5) obtained in Proposition 4.2.1 is characterized by the existence of σh ∈ Λ h such that the local equations (4.2.3), (4.2.4) hold. Proof. Assume first that (û h , σh ) ∈ V h × Λ h satisfy (4.2.3), (4.2.4). Consider w h ∈ V h . Then by (4.2.4) one hasF (Dw h ) ≥ F (Dû h ) + σh : (Dw h -Dû h ) a.e. in Ω.Using the same argument for t < 0 we obtain the converse inequality, we we deduce that (4.2.38) is indeed an equality. Since ûh ∈ V h is continuous and piecewise affine, Dû h is piecewise constant, i.e. Dû h ∈ Λ h . Setting σh = F (Dû h ) ∈ Λ h , we obtain (4.2.3), (4.2.4).

	(4.2.35)
	Integrating over Ω we obtain

Ω F (Dw h ) ≥ Ω F (Dû h ) + σh , Dw h -Dû h . (

4

.2.36) Taking in (4.2.3) v h = w h -ûh and using the previous inequality we get (4.2.5) with the test function w h . Conversely, knowing that there is existence and uniqueness for (4.2.5), we only have to prove that there exists a solution (û h , σh ) ∈ V h × Λ h to (4.2.3), (4.2.4). Let us consider first the case when F is continuously differentiable. Then we have a solution ûh to (4.2.5). Consider a test function w h ∈ V h , and take v h = ûh + tw h for t = 0, in (4.2.5). We obtain 1 ∆t ûh -u n , tw h + η Dû h -Du n , tDw h + Ω F (Dû h + tDw h ) -Ω F (Dû h ) ≥ f n , tw h . (4.2.37) Dividing by t > 0 and letting t → 0 we obtain using Lebesgue's theorem since Dû h , Dw h belong to L ∞ , 1 ∆t ûh -u n , w h + η Dû h -Du n , Dw h + Ω F (Dû h ) : Dw h ≥ f n , w h . (4.2.38)

  (3.2.48), (3.2.49) if η = 0. Moreover one has Ω

. Proposition 4.2.5. With the same assumptions as in Proposition 4.2.2 and if Ω is strictly starshaped, when h → 0 one has ûh → û in H. According to Proposition 4.2.4 there exists some σh ∈ Λ h such that the local equations (4.2.3), (4.2.4) hold. Then after extraction of a subsequence one has σh σ in weak L 2 (Ω), where (û, σ) is a solution to (3.1.32) if η > 0, or to (3.2.46), (3.2.47),

  .2.53) Theorem 4.2.6. We assume that η ≥ 0, and that F is convex, finite everywhere and lower bounded. We denote by ûh ∈ V h the solution to (4.2.5) obtained in Proposition 4.2.1, and we assume that Then the sequence defined by the algorithm (4.2.52)verifies u k h → ûh in V h as k → ∞.Moreover if there is strict inequality in (4.2.54) then σ k h is bounded in Λ h , and up to extraction of a subsequence one has σ k h → σh in Λ h , where (û h , σh ) ∈ V h × Λ h solves (4.2.3), (4.2.4), as obtained in Proposition 4.2.4. Proof. We proceed as exposed in the formal proof of Proposition 4.1.1. Due to the definition (2.4.1) of P r , one has

	rτ L 2 h ≤ 1.	(4.2.54)

Table 5 .

 5 1: Test with various values of rτ Optimal constant C in the formula rτ = Ch 2 In the convergence of the primal-dual algorithm, Theorem 4.2.6 requires the stability condition rτ L 2 h ≤ 1 where L h is defined according to (4.2.65) by

		5.3: Space discretisation
	rτ	Number of iterations Error
	> h 2 /6	NaN	NaN
	h 2 /6	4103	1.0683
	h 2 /8	4407	0.1068
	h 2 /16	5076	0.1067
	h 2 /32	6270	0.1067

Table 5 .

 5 αu -div σ = f one has merely [f ] ∼ α[u]. For a time-dependent problem one can use instead [u] = [u 0 ] + t[f ], where u 0 is the initial datum. In Tables 5.3, 5.4, 5.5, we give the error and the cost for several choices of r and τ and different values of h. The optimal value in the three tables is r = h 16 and τ = 2h. Notice that the choice for 3: Number of iterations and error for different choices of τ , r for h = 1/20, ε tol = 10 -3 .

	.2.19)

Table 5 .

 5 4: Number of iterations and error for different choices of τ , r for h = 1/40, ε tol = 10 -4 .

Table 5 .

 5 5: Number of iterations and error for different choices of τ , r for h = 1/80, ε tol = 10 -5 .

		-2 err = 0.332847 n.o.iter = 153 execution time = 26.4396s
	ε = 10 -3 err = 0.113881 n.o.iter = 418 execution time = 69.0468s
	ε = 10 -4 err = 0.10936 n.o.iter = 746 execution time = 112.164s
	ε = 10 -5 err = 0.109946 n.o.iter = 1055 execution time = 172.84s
	(τ, r)	h 2 8 , 1	h 16 , 2h	h 8 , h	h 4 , h 2	h √ 8 , h √ 8	h 2 , h 4	h, h 8	2h, h 16	1, h 2 8
	N o of iterations 31915	5174	3218	2466	2162	1858	1283	1153	18882
	Error	0.5222 4.9833e-2 4.6645e-2 4.6131e-2 4.5556e-2 4.4453e-2 4.0351e-2 3.2291e-2 3.3465e-2

Table 5 .

 5 6: Number of iterations and error for the Accelerated scheme with different choices of τ 0 , r 0 for h = 1/20, ε tol = 10 -3 .

	(τ 0 , r 0 )	h 2 8 , 1	h 16 , 2h	h 8 , h	h 4 , h 2	h √ 8 , h √ 8	h 2 , h 4	h, h 8	2h, h 16	1, h 2 8
	N o of iterations	1843	517	484	365	365	365	367	368	370
	Error	0.2497 5.1105e-2 2.9266e-2 1.8845e-2 1.6860e-2 1.6701e-2 1.7926e-2 1.8984e-2 1.9894e-2

Table 5 .

 5 7: Number of iterations and error for the Accelerated scheme with different choices of τ 0 , r 0 for h = 1/40, ε tol = 10

Table 5 .

 5 1 8: Number of iterations and error for the Accelerated scheme with different choices of τ 0 , r 0 for h = 1/80, ε tol = 10 -5 .

			h 16 , 2h	h 8 , h	h 4 , h 2	h √ 8 , h √ 8	h 2 , h 4	h, h 8	2h, h 16	1, h 2 8
	N o of iterations 13172	2284	1812	1660	1659	1660	1661	1663	1666
	Error	0.1581 2.4531e-2 1.5232e-2 1.0248e-2 9.5125e-3 9.4371e-3 9.8587e-3 1.0252e-2 1.0655e-2

Table 5 .

 5 9: Number of iterations and error for normal and accelerated schemes with different choices of τ , r for h = 1/20, ε tol = 10 -3 , α = 100

	(τ, r)	h 2 8 , 1	h 16 , 2h	h 8 , h	h 4 , h 2	h √ 8 , h √ 8	h 2 , h 4	h, h 8	2h, h 16	1, h 2 8
	N o of iterations	932	99	120	161	195	243	406	693	3
	Error	2.7232e-2 1.0971e-2 1.3108e-2 1.6009e-2 1.7693e-2 1.9275e-2 2.2429e-2 2.6756e-2 0.1324
	(τ 0 , r 0 )	h 2 8 , 1	h 16 , 2h	h 8 , h	h 4 , h 2	h √ 8 , h √ 8	h 2 , h 4	h, h 8	2h, h 16	1, h 2 8
	N o of iterations	985	244	192	166	164	164	165	166	167
	Error	0.1116 2.3446e-2 1.4898e-2 1.0437e-2 9.5771e-3 9.3094e-3 9.3612e-3 9.4723e-3 9.5632e-3

Table 5 .

 5 10: Number of iterations and error for the normal and accelerated schemes with different choices of τ , r for h = 1/40, ε tol = 10 -4 , α = 100 error, and this is worse for smaller h. However if we choose (τ, r) = ( h

	√	8 , h √ 8

Table 5 .

 5

	1, h 2 8

11: Number of iterations and error for the normal and accelerated schemes with different choices of τ , r for h = 1/80, ε tol = 10 -5 , α = 100

  The case with or without degrees of freedom on the boundary (depending on the type of boundary conditions) are shown on Figures 5.7 and 5.8.

	5.3.2 Numerical results for the 2D Euler/Bingham model We now consider the 2D Euler transport together with the viscoplastic model Bingham rheology ∂ t ρ+ div(ρu) = 0, (5.3.8) error number of iterations ε tol 20 0.763 31 10 -2 40 0.2219 60 10 -4 80 0.11 120 10 -5 160 6.706e -2 241 10 -6 ∂ N 320 3.817e -2 482 10 -7

t (ρu)+ div(ρu ⊗ u + p(ρ)I) -div Du |Du| = f, (5.3.9) Figure 5.10: Euler/Bingham model: approximate and exact velocity u x Figure 5.11: Euler/Bingham model: approximate and exact velocity u y

Table 5 .

 5 19: Error and number of iterations with the 2d code

	N	error estimate
	32	7.07e -2
	75	4.541e -2
	150	2.695e -2
	300	3.87e -3
	600	3.03e -3
	1200	1.849e -3

Table 5 .

 5 20: Error obtained with the 1d code

  Figure 6.2: Qualitative behaviour of various types of fluids. Left: stress vs shear stress -Right: apparent viscosity vs shear stress. The Bingham type fluids can be viewed as enjoying infinite apparent viscosity below the threshold τ c .

	Shear stress τ [Pa]	B i n g h a m B i n g h a m	p s e u d o p l a s t i c	N e w t o n i a n d i l a t a n t	s -1 ] Viscosity ν [m 2	d il a ta n t
						Newtonian
	τ c		ps eu do pl as ti c
						p s e u d o p l a s t i c
			Shear rate γ [s -1 ]		Shear rate γ [s -1 ]

p b (p 0 -p b ) δp b .(2.6.24)

Acknowledgements

First and foremost, I would like to express my deepest and most sincere gratitude to my supervisor

gia đình, vừa là tấm gương, của em, anh chị, bản thân. Cảm trong những năm qua. tới gia đình cô Thúy, chú Tiến, cô Nội, một gia đình thứ hai, có thêm nhiều Nội, Minh Hà, ra cho em đón em mỗi khi em trở về. thầy, không bắt đầu.

iii Abstract This thesis is motivated by a research program between the LAMA (Mathematics, Univ. Gustave Eiffel) and the Institut de

It follows that

Now, for a given integer n one can find λ n > 1 such that

Since Ω is strictly star-shaped, one has ∂Ω n ∩ ∂Ω = ∅. Furthermore, v λn has support in Ω n , and dist(Ω n , ∂Ω) > 0. Therefore for δ small enough, w δ ∈ C ∞ c (Ω). Thus we can find such small δ = δ n so that

Taking into account that both w δn and u have support in Ω, (3.

(3.2.33) Proposition 3.2.6. Under Hypothesis 1, 3, 4, 5, we have that ψ is convex, proper and lower semi-continuous on L 2 (Ω).

Proof. (a) ψ is convex. Assume that u, v ∈ L 2 (Ω) and θ ∈ (0, 1). If u ∈ K F 0 or v ∈ K F 0 , then ψ(u) = ∞ or ψ(v) = ∞, and it follows obviously that ψ((1 -θ)u + θv) ≤ ∞ = (1 -θ)ψ(u) + θψ (v) We can find a subsequence such that D i

V , for some unit vector V (here we use that we are in finite dimension, otherwise V could be null). Then for σ such that |σ| > C and σ |σ| = V , one has F * (σ) = ∞, which finishes the proof.

Taking the opposite we obtain

Hence, adding (4.1.5) to (4.1.4) one gets for all (u, σ)

Besides, (û, σ) being the exact solution to (4.1.2), we have

and applying the same estimate as in (4.1.5),

From (4.1.7), (4.1.8) one has for all (u, σ)

Substituting σ = σ, u = û in (4.1.6), we can then subtract to (4.1.9) with σ = σ k+1 and u = u k+1 , to obtain Numerical experiments

1D Bingham model with Euler transport

In the numerical experiments, additionally to the viscoplastic rheology we include transport (inertial) terms. With such terms, incompressible viscoplastic models were considered years ago, and lots of explicit or reference solutions are available, as well as laboratory experiments. One of the obstacles in the numerical evaluation of compressible models with Euler transport terms is the lack of explicit solutions. In this first test we propose an analytical solution for the one-dimensional compressible Euler equations with Bingham rheology. This will be used as a particular 1d solution to our 2d model. In the 1d model the unknown (ρ(t, x), u(t, x)) has to satisfy

where p(ρ) = 1 2 ρ 2 , Ω = (0, L) and f (t, x) is a given force term. We take L = 4, σ 0 = 1, T = 1 and we build an exact solution such that u has the form

(5.1.2) Indeed ρ and f will be derived from this choice. In order to find the density we look for the flow X(s, t, x) that satisfies dX ds = u(s, X), X(t, t, x) = x.

(5.1.3)

Classically, a function V (t, x) satisfies the linear transport equation

if and only if V (t, x) = V 0 (X(s = 0, t, x)). In our case of u given by (5.1.2) a (not so simple) computation gives for s, t ≥ 0 and 0 ≤ x ≤ 4

+ if ≤ x ≤ 4. Then by taking s = 0 and assuming t ≤ 2 we get

(5.1.5)

Taking the spatial derivative of (5.1.4), we get ∂ t ρ + ∂ x (ρu) = 0 with ρ = ∂ x V . Thus ρ(t, x) is given by ρ(t, x) = ρ 0 (X(s = 0, t, x))∂ x X(s = 0, t, x). Taking ρ 0 (x) = 1 we obtain

(

Then we compute

and we take sgn

which implies that

(5.1.9)

Therefore f is finally deduced as

We denote u h = (u i ) i∈I = (u x i , u y i ) i∈I where I is a set of indices corresponding to the nodes in the domain. For n ∈ I, the finite volume cell around the node n is denoted by Q n . Then one has

Then, writing vh in terms of basis function as vh = n∈I v n ϕ n where v n is a constant vector and where ϕ n is an affine function such that ϕ(x n ) = 1 and ϕ(x m ) = 0 for all m = n, one has

By choosing vh = 1 0 ϕ n in (5.2.10) we get for n ∈ I

Similarly by choosing vh = 0 1 ϕ n it gives for n ∈ I

Therefore we can compute explicitly

(5.2.16)

An illustration is on Figure 5.

3(right).

Primal-dual algorithm on a structured mesh using Fortran. To avoid any geometric complexity we use the above structured mesh together with Fortran coding for the primal-dual algorithm (without acceleration). The results are shown on Figure 5.4.

(τ, r) We have the corresponding variational formulation

For Dirichlet boundary conditions we take

For Neumann boundary conditions σn = 0, we take

For slip boundary conditions u • n = 0, (σn) × n = 0, we take

Dirichlet condition. We use the radial solution

where r = x 2 + y 2 , Ω = (-1, 1) × (-1, 1), that has been used in subsection 5.2.1.

Neumann condition. We use the radial solution

which is a solution to (5.2.22) with α = 1 and with Neumann condition if we add 3 5 to the previous f . We use 161 points in each direction. The stopping criterion is 10 -6 . The numerical results give the error 4.63.10 -2 , and the number of iterations 3740. This last condition holds trivially. Similarly we get u x | x=-1 = 0. We take

(5.2.36)

(5.2.37) than the regularisation method, but is faster than the augmented Lagrangian method. In fact the augmented Lagrangian method can take more iterations. This maybe due to the fact that the penalty parameter r is not chosen optimally, indeed this is difficult to know the best value. We conclude that for this test the regularisation method is the best, knowing that in FreeFEM++ linear systems are solved very efficiently. Moreover here α = 100 thus the problem is not stiff. We have to remember also that in some cases the regularisation method has difficulties to well localise the solid zones, what the primal-dual algorithm is supposed to do well. We consider two choices τ = r = h 3 or τ = h∆t 3 , r = h 3∆t . We take t max = 0.5, h = 1/20, ε tol = 10 -3 . According to Tables 5.17 and 5.18 we can see that the first choice is better in terms of error, but the second one is less costly in terms of iterations. In both cases the value of the error depends very little on the timestep, this maybe due to the fact that the solution u is linear in time. More experiments would be necessary in order to better achieve the balance between τ and r, and decide whether the first or the second choice is the best.

Compressible Euler equations with Bingham viscoplasticity

Numerical scheme for the hyperbolic part

As described in the introduction, the compressible Euler equations with Bingham viscoplasticity can be solved by the splitting method. For the Euler hyperbolic part we have to solve

with p = p(ρ). We shall take p(ρ) = ρ 2 /2. This can be written as

where U = (ρ, ρu) T and F (U ) = (ρu, ρu ⊗ u -pI) T . In order to solve the system (5.3.3) we use the finite volume method, which at first order consists in updating the value of U n i corresponding to each cell Q i by the formula

where ∆t is the timestep, As we can observe on Table 5.19, by choosing the optimal value of ε tol the scheme is first-order accurate with respect to the number of points in each direction, as with the 1d code (Table 5.20). This is true even though u y exact = 0 but u y approx is not zero.

Chapter 6

A lubrication equation for a simplified model of shear-thinning fluid Boundary conditions are: z = ϕ: fluid-atmosphere interface. We have continuity of the stress tensor at the free surface, together with a kinematic boundary condition. Since the atmosphere can be viewed as an ideal fluid, the stress tensor can be taken equal to zero above ϕ. Hence we get

z = f b : interface between the fluid and the substrate, which is fixed. This is a material interface, on which we have the no-slip boundary condition

Here (u b , v v ) is the so-called basal velocity. Often in fluid mechanics the basal velocity is zero, but for geophysical applications it can actually be the driving force, and thus depend on (t, x).

Rheology

For a fluid, the deviatoric stress tensor τ is usually a function of the strain rate tensor

A Newtonian fluid is characterized by a linear relation, defining the viscosity of the fluid, which is assumed here to be isotropic and constant. Therefore we introduce the dynamical viscosity coefficient µ, and define the Newtonian stress tensor by

where ν = µ/ρ is the kinematic viscosity. Fluids that do not follow this kind of constitutive law are non-Newtonian. In the general case, the material invariance principle implies that the stress tensor depends only on the similarity invariants of the strain rate tensor, in particular the coefficients of its characteristic polynomial. In dimension 2 there are only two such coefficients εI and εII . Namely εI is the trace of the matrix and εII its determinant. For an incompressible fluid, the trace is zero, and moreover we have

This allows to define the strain rate γ as

In a similar way we can check that the Frobenius norm of ε, that is ε

A very sketchy illustration of the possible behaviours of non-Newtonian fluids is given in Figure 6.2. We will be mostly interested in this work in the so-called pseudoplastic case, that is the red curve in Figure 6.2, for which experimental evidence can be given, see [START_REF] David | Demonstration of upper and lower Newtonian fluid behaviour in a pseudoplastic fluid[END_REF]. This kind of models are also used in geophysics, see [START_REF] Bercovici | 7.07-the generation of plate tectonics from mantle dynamics[END_REF][START_REF] Paul | Self-consistent generation of tectonic plates in time-dependent, threedimensional mantle convection simulations 2. Strain weakening and asthenosphere[END_REF][START_REF] Zhong | Role of faults, nonlinear rheology, and viscosity structure in generating plates from instantaneous mantle flow models[END_REF] We wish to give a simplified model for this pseudo-plastic At this stage, we introduce two classical dimensionless quantities, namely the Froude and Reynolds numbers, defined from the characteristic horizontal velocity u 0 , the vertical extension h 0 , and the viscosity for large deformations ν:

We divide the previous two equations by u 2 0 /h 0 , and noticing that τ = ρReu 2 0 τ , we obtain

The idea now is to send δ to zero, thus implementing the thin layer assumption, but in a regime where the Reynolds number Re is kept of order 1, together with a balance between viscosity and gravity forces. Therefore we set F r 2 = δRe. (6.2.21)

This readily gives

which is the scaling proposed in [START_REF] Neil | Viscoplastic flow over an inclined surface[END_REF]. It introduces another characteristic velocity, namely u 0 = (gh 2 0 )/ν. The latter equality shows that this is indeed a slow motion scaling, thus we meet the initial requirement.

Inserting (6.2.21) in equations (6.2. [START_REF] Chalayer | A bi-projection method for incompressible Bingham flows with variable density, viscosity, and yield stress[END_REF]) and (6. Now we compute τ from (6.2.12). We start by rewriting ε in rescaled variables

From this we easily deduce τ = 1 2

We define a dimensionless equivalent viscosity by νeq = ν eq /ν, and rewrite equation (6.2.23)

.2.27)

We turn now to the expression of νeq . We first notice that, using (6.2.25)

The advection flux is merely an upwind flux

For the diffusive flux, we write G n j+1/2 = D n j+1/2 K n j+1/2 , where K n j+1/2 is the approximate value of the slope ∂

and D n j+1/2 is a discretization of (6.3.15). To obtain it we need to compute h * and y * at the interface. Accordingly to (6.3.16), we put

The time step ∆t is actually updated at each time step using the CFL condition

, where σ < 1.

The following simulations have been performed with J = 200 cells in the interval [-1, 1], together with σ = 0.9. All figures are gathered at the end of the paper. The first set of simulations concerns the collapse of a square-shaped stack on a horizontal flat bottom: h 0 (x) = 1 for x ∈] -1/3, 1/3[, 0 elsewhere, with zero basal velocity (u b = v b = 0). We first propose a comparison between the two viscosities model and the high viscosity and low viscosity models. The small deformation viscosity is ν B = 100 (recall that ν = 1), and the yield stress is 0.1 in Figure 6.4, and 0.5 in Figure 6.5. These figures are complemented by Figure 6.6 where we display for four values of the yield stress B a time-lapse of the evolution of both the total thickness of the fluid h (plain lines) and the thickness of the low velocity layer (dashed lines).

For B = 0.1, the fluid clearly behaves similarly as the low viscosity fluid in the early stages, then eventually it slows down, when the low viscosity layer tends to disappear, see Figure 6.6, top left. With a yield stress B = 0.5, the two viscosities model stays in-between the other two, as expected, faster that the high viscosity model, slower than the low viscosity one, see Figure 6.5. However, one can check that the front hardly moves between t = 10 and t = 50, indicating that the fluid tends to behave as the high velocity one. This is made more explicit in Figure 6.6, top right, where for t = 10 and t = 50 the low viscosity layer is very small. In general, the thickness y * decreases with time, faster when B is larger. It is hardly observable fot t = 50 when B = 2.5, indicating that the fluid is almost completely driven by the high viscosity.

Using the same initial data, we check the convergence of the two viscosities model towards the Bingham fluid when ν B goes to ∞. We take a yield stress B = 1.25, and ν B = 10, 100, 1000. As expected, the behaviour becomes close to the Bingham fluid, yet it departs from it for larger times, see Figure 6.7.

We turn now to a different context, closer to the situation in geophysics. The flow here is no longer purely gravity driven, it is actually dragged along by a non zero basal velocity. The idea here is that our pseudo-plastic fluid is a very crude model of some planetary lithosphere, below which lies the mantle. The basal velocity is the upper trace of convection currents in the mantle, which are supposed to be the main drivers of plate tectonics. The initial thickness is constant equal to 1, and we use two basal velocities U b (x) = (u b (x), 0), where These velocities crudely correspond respectively to the vertical motion of a magma bubble, which generates local perturbations of the velocity. The first one corresponds to some bubble lift, with negative velocity on the left and positive on the right. It generates some kind of a valley surrounded by mountains, see Figure 6.8. Conversely, the descent of a bubble reverses the velocities, and produces a mountain surrounded with valleys, Figure 6.9. We notice in both cases that the small velocity model has very little influence on the time evolution, and that the two viscosity model leads to rather sharp angles in the thickness.