Frédéric Golay
email: golay@isitv.univ-tln.fr

Pierre Seppecher

Locking materials and the topology of optimal shapes

We study a problem of structural optimisation using the fictitious material approach. This is connected with the equilibrium of locking materials, which can be approximated by strongly non linear materials. A finite element simulation allows us to experiment with some conjectures about the topology of the optimal solutions.

Introduction

Structural optimisation is a major objective in the conception of industrial systems. In particular, for a given desired performance, engineers may need to minimise the mass of a structure by using an adapted geometry. Early methods consist of progressively tuning some geometric parameters of a postulated initial shape. The optimal design obtained in this way depends strongly on the initial one: the optimisation process cannot change the fundamental topological characteristics of the shape. Methods which allow such modifications [START_REF] Allaire | Shape optimisation by the homogenisation method[END_REF][START_REF] Bendsoe | Optimisation of Structural Topology, Shape and Material[END_REF] are called topological optimisation. Here we consider a structure submitted to a given load and we try to get a structure with the smallest volume for a desired compliance (i.e. global stiffness, or stored elastic energy for the given load). It is now well known that this problem may admit no solution in the classical sense: the optimal 'shape' may consist of an intricate mixture of material and holes. The optimal solution has to be understood in the framework of homogenisation theory and the optimisation problem has to be set at the very beginning in a relaxed form. However, as the set of all possible effective stiffness tensors for a given volume fraction of the material is rather intricate [START_REF] Allaire | Optimal design for minimum weight and compliance in plane stress using extremal microstructure[END_REF][START_REF]Etude théorique et numérique de quelques problèmes d'optimisation de formes à l'aide de méthodes d'homogénéisation[END_REF][START_REF] Bendsoe | Optimisation of Structural Topology, Shape and Material[END_REF], it seems difficult to find analytic solutions or mathematical properties for optimal designs. That is why we consider a very similar but much simpler problem: we assume that the stiffness tensor depends linearly on the volume fraction (or the mass density) of the material. In 3-D this assumption is not realistic as the stiffness of a composite material with a given volume fraction is weaker than the stiffness obtained assuming a linear dependence: the problem we consider is then called 'fictitious material optimisation'. In 2-D, the problem corresponds to the optimisation of the thickness of a plate submitted to a plane load (plane stress). These linear dependence or different 'artificial mixture' rules are discussed in [START_REF] Duysinx | Optimisation topologique: du milieu continu à la structure élastique[END_REF].

The optimal designs obtained using homogenisation theory or the fictitious approach are qualitatively similar. Hopefully our results could then be extended to the optimal structures obtained by the homogenisation method.

We show that the problem of fictitious material optimisation is equivalent to the equilibrium problem of a perfect locking material (i.e. a material whose strain tensor must belong to a given bounded set, in this case its internal energy density vanishes). This convex minimisation enables us to give non-trivial analytical solutions (section 4.1). Its regularisation leads to a simple non-linear elasticity problem (section 3). Some 2-D numerical solutions are presented (section 5) and compared to analytical solutions. The numerical solutions enable us to test some conjectures about the topology of optimal solutions. In the case when the support of the optimal design is a connected set, it seems that this support is a simply-connected bounded set. In particular, we can prove that circular holes cannot be present in an optimal design. This conjecture, if confirmed, is remarkable: the topological complexity is not introduced by the optimisation process, but only by a later penalisation process (such a penalisation is generally used to get out of the relaxed formulation and to obtain a classical solution [START_REF] Allaire | Shape optimisation by the homogenisation method[END_REF]. Only this penalisation process justifies the name of 'topological optimisation' for the whole process.

Optimal design of membranes

Mathematical formulation

Let us consider a thin plane membrane of variable thickness h(x), x belonging to some bounded domain ⊂ R 2 . Its volume V is then given by:

V = h(x) ds.
(1)

Assuming that the membrane is made of an homogeneous isotropic elastic material with Young's modulus E and Poisson coefficient ν, its elastic energy under plane stress conditions is:

1 2 ε v (x) : D : ε v (x)h(x) ds, (2)
where v ∈ H 1 ( , R 2 ) denotes the displacement field of the membrane, ε v the plane strain tensor ε v = 1 2 (∇v + ∇ t v) and D is the two-dimensional fourth-order stiffness tensor:

D ij kl = Eν 1 -ν 2 δ ij δ kl + E 1 + ν δ ik δ jl .
At equilibrium the displacement field u minimises the potential energy:

J (v, h) = 1 2 ε v : D : ε v h ds - L F • v d , ( 3
)
where F denotes a line density of applied forces along a line L. Other types of measures could be considered for F , but, for the sake of simplicity, this will not be done here. For equilibrium to be possible, we assume that the resulting force and moment vanish ( L F d = L x ∧ F d = 0). Such conditions could be replaced by Dirichlet conditions for the displacement on some part of .

The elastic energy at equilibrium can be rewritten, in terms of the equilibrium displacement field u,

E(h) = -inf v J (v, h) = -J (u, h) = 1 2 ε u : D : ε u h ds = 1 2 L F • u d . ( 4
)
Assume that we need to maximise the global stiffness of the structure (or which is equivalent, to minimise the volume of the structure) for the given single loading F . An optimal design of the membrane is a function h(x) which, for a given volume V, minimises the equilibrium elastic energy. The optimisation problem reads:

E = inf h E(h); h(x) 0; h(x) ds = V , (5)
or, using (3) and ( 4):

E = -sup h inf v 1 2 ε v : D : ε v h ds - L F • v d ; h 0; h ds = V . ( 6
)
Let us remark that, even in this relaxed formulation, the existence of an optimal design (a function h minimising (6)) is generally not assured. For some special loading F the optimal solution may concentrate along a line: indeed, as proved in [START_REF] Bouchitté | Characterisation of optimal shapes and masses through Monge-Kantorovich equation[END_REF], the optimal solution h has to be sought in the space of measures. Examples of this phenomenon are given in (Bouchitté et al., 1997-a) for conductivity optimisation and can easily be extended to our problem. In this paper we assume that forces are such that h and the associated displacement field u exist and are sufficiently regular.

Note that the line L, where external forces are applied, may be different from the boundary ∂ of and may not even coincide with a part of the boundary. This allows us to consider the case = R 2 or, in numerical studies, to consider a domain much larger than the support of applied forces. Then we can obtain optimal structures which are not subjected to any geometrical constraint.

Perfect locking materials

We show now that the problem (6) is equivalent to the resolution of the equilibrium of a perfect locking material. We search, a priori, the supremum in (6) in the convex compact set M of all non-negative measures h on ¯ which verify h = V. By a density argument, the infimum can be sought in C ∞ ( ¯ ), so (6) can be rewritten:

E = inf h∈M sup v∈C ∞ ( ¯ ) -J (v, h) . (7)
As, on C ∞ ( ¯ ) × M, -J is concave with respect to the first variable and convex (linear) lower semicontinuous with respect to the second one, the 'lop sided minimax theorem' (cf. chapter 6 of [START_REF] Aubin | Applied Nonlinear Analysis[END_REF]) allows us to interchange inf and sup in ( 7):

E = sup v∈C ∞ ( ¯ ) inf h∈M -J (v, h) = -inf v∈C ∞ ( ¯ ) sup h∈M J (v, h). ( 8
)
This last supremum in h is easy to compute (it is enough to concentrate h where ε v : D : ε v ∞ reaches its maximum), and we get:

E = -inf v∈C ∞ ( ¯ ) - L F • v d + ε v : D : ε v ∞ V 2 . ( 9
)
In order to get rid of the influence of the total volume V, we can consider, for each v ∈ C ∞ ( ¯ ), the quantities

s = ε v : D : ε v 1/2
∞ and w = s -1 v and write (9) as:

E = - inf w∈C ∞ ( ¯ ),s 0 -s L F • w d + Vs 2 2 , ε v : D : ε v ∞ = 1 . ( 10
)
As the minimisation in s is straightforward, we obtain:

E = - 1 2V inf w L -F • w d ; ε w : D : ε w ∞ = 1 2 , ( 11
)
or, in an equivalent way,

E = - 1 2V inf w L -F • w d ; ε w : D : ε w ∞ 1 2 . ( 12
)
In this last formulation, the minimisation problem corresponds to the equilibrium of a perfect locking material [START_REF] Demengel | On locking materials[END_REF]: the strain tensor of such a material has to satisfy the constraint ε w : D : ε w ∞ 1, then its volume energy vanishes. These results also hold in the more general case when the elastic material is a non-linear one. This has been proved by G. [START_REF] Bouchitté | Characterisation of optimal shapes and masses through Monge-Kantorovich equation[END_REF].

Note that formulation ( 12) can be obtained by using the stress based formulation of [START_REF] Allaire | Optimal design for minimum weight and compliance in plane stress using extremal microstructure[END_REF]. In our simpler case the set of all stiffness tensors for a given thickness h, is reduced to the singleton hD. Then, the variational problem (2.16) of [START_REF] Allaire | Optimal design for minimum weight and compliance in plane stress using extremal microstructure[END_REF] is easily explicated: in our case this problem becomes convex and corresponds to the dual formulation of (12).

A similar formulation has already been used to study optimal truss structures [START_REF] Michell | The limits of economy of material in frame structures[END_REF][START_REF] Lagache | Treillis de volume minimal dans une région donnée[END_REF][START_REF] Rozvany | Structural Design via Optimality Criteria[END_REF]. Indeed, problems (P2) or (P3) in [START_REF] Lagache | Treillis de volume minimal dans une région donnée[END_REF] clearly correspond to the equilibrium of a perfect locking material. In the case of optimal trusses, the strain tensor has to verify the constraint |||ε w ||| 1 (where |||ε||| denotes the highest singular value of ε) while in our case the constraint is ε w : D : ε w ∞ 1. The locking material formulation for optimal trusses has been particularly fruitful and optimal solutions can be described in terms of Hencky nets [START_REF] Strang | Hencky-Prandtl nets and constrained Michell trusses[END_REF]). Here we investigate how the results obtained for optimal trusses can be extended to the problem (12). We cannot obtain a simple structure like Hencky nets, nevertheless the formulation (12) enables us to get some analytical solutions and some topological properties of optimal designs. Formulation (9) has already been used in finite dimension by [START_REF] Ben-Tal | A new method for optimal truss topology design[END_REF] for the numerical analysis of optimal trusses and by [START_REF] Zowe | Free material optimisation via mathematical programming[END_REF] for the finite element study of fictitious material. We will first regularise our convex minimisation problems before discretising them.

Numerical formulation

Problem regularisation

In order to solve numerically (9) or ( 12), let us regularise the term containing the norm • ∞ . In the first formulation, we approximate the infinity norm by a p-norm where p is a large parameter. Indeed:

ε v : D : ε v ∞ = lim p→∞ ε v : D : ε v p = lim p→∞ (ε v : D : ε v ) p ds 1/p . ( 13
)
Then we deal with the following optimisation problem:

lim p→∞ inf v J p (v) , ( 14
)
where

J p (v) := V (ε v : D : ε v ) p ds 1/p -2 L F • v d . ( 15
)
As the derivative of J p is given by:

J p (u) • v = V (ε u : D : ε u ) p ds 1/p-1 2(ε u : D : ε u ) p-1 (ε u : D : ε v ) ds -2 L F • v d ,
the variational formulation of the optimisation problem becomes:

H (u)ε u : D : ε v ds = L F • vd , ( 16
)
where

H (u) := V(ε u : D : ε u ) p-1 (ε u : D : ε u ) p ds 1/p-1 . ( 17
)
This variational formulation corresponds to a non-local elastic problem. Indeed the stiffness tensor H D, through H , depends on the solution u on the whole domain.

For a given H , this formulation corresponds again to the equilibrium of an elastic plate of thickness H under plane stress. The comparison of ( 17) and ( 4) shows that H corresponds to the optimal design h we are looking for. Here H depends on the global strain energy.

In the second formulation (12), introducing the indicator function defined on R + by:

I (t) := ∞ if t > 1, 0 otherwise, ( 18
)
the optimisation problem becomes:

inf w - L F • w d + I (ε w : D : ε w ) ds . ( 19
)
Let us now approximate the indicator function I by the function I p , defined by I p (t) := (2p) -1 t p , p being a large parameter. The minimisation problem becomes:

lim p→∞ inf w Jp (w) , ( 20
)
where

Jp (w) := 1 2p (ε w : D : ε w ) p ds - L F • w d . ( 21
)
As the derivative of Jp is given by:

J p (u) • v = (ε u : D : ε u ) p-1 (ε u : D : ε v ) ds - L F • v d , ( 22
)
the variational formulation of the optimisation problem is written again:

H (u)ε u : D : ε v ds = L F • v d , ( 23
)
where

H (u) := (ε u : D : ε u ) p-1 . ( 24
)
This variational formulation corresponds to a simple local elastic problem. Again the comparison of ( 23) and ( 4) shows that H corresponds to the optimal design h we are looking for. Here H depends simply on the local strain energy density. Note that the problems ( 16)-( 17) and ( 23)-( 24) are equivalent up to a multiplicative constant. They are both non linear, but convex, elastic equilibrium problems. They lead to the same optimal design. Similar formulations have been established by Ben-Tal and [START_REF] Ben-Tal | A new method for optimal truss topology design[END_REF] for the study of optimal trusses or by [START_REF] Petersson | A finite element analysis of optimal variable thickness sheets[END_REF] for the study (very closely connected to ours) of optimal variable thickness sheets.

Finite element formulation

We use a finite element formulation [START_REF] Batoz | Modélisation des Structures par Éléments Finis 1[END_REF] for our problems. Let us denote {u} the vector of the nodal unknowns and, for each element e, N e the interpolation functions. Then the displacement field is approximated on e by: u(x) = N e (x) {u}, which, by differentiation, defines B e such that:

ε u (x) = B e (x) {u}.
Then we write on each element e:

ε u : D : ε v = {u} T B e T [D] B e {v} = {u} T K e {v}.
Here [K e ] := [B e ] T [D][B e ] is the elementary elastic stiffness matrix. Decomposing the domain with the supports V e of the finite element and denoting:

E p := (ε u : D : ε u ) p ds = e V e
{u} T K e {u} p ds, the finite element formulation of ( 16) is written as:

VE 1/p-1 p e V e
{u} T K e {u} p-1 {v} T K e {u} ds = e ∂V e {v} T N e {F } d .

As this equality holds for any {v}, the residue:

{R p } {u} := - e ∂V e N e {F } d + V e V e {u} T K e {u} p ds 1/p-1 × e V e
{u} T K e {u} p-1 K e {u} ds, (25) must vanish.

In the same way the residue of the finite element formulation of (23) becomes:

Rp {u} := e V e
{u} T K e {u} p-1 K e {u} dse ∂V e N e {F } d .

(26)

Computational method

As, for a given value of p = 1, the problems:

{R p } {u} = {0}, ( 27
) Rp {u} = {0}, (28)
are non linear, we use the Newton-Raphson method. At each iteration j we solve the following linear system:

∂{R} ∂{u} u j -1 {δu} = -{R} u j -1 , u j = {δu} + u j -1 . ( 29
)
For p = 1 the problem is linear. Increasing p progressively, we initialise the Newton-Raphson algorithm by choosing for {u 0 } the displacement field computed with a smaller value of p.

The differentiation of formula (28) leads to the following tangent matrix:

∂{ Rp } ∂{u} = e V e
{u} T K e {u} p-1 K e ds

+ 2(p -1) e V e
{u} T K e {u} p-2 K e {u} K e {u} T ds,

which we can use directly while the differentiation of formula ( 27) leads to a full tangent matrix:

∂{R p } ∂{u} = V e V e {u} T K e {u} p ds 1/p-1 e V e {u} T K e {u} p-1 K e ds + 2(p -1)V e V e {u} T K e {u} p ds 1/p-1 × e V e
{u} T K e {u} p-2 K e {u} K e {u} T ds (31)

+ 2(1 -p)V e V e {u} T K e {u} p ds 1/p-2 × e V e
{u} T K e {u} p-1 K e {u} ds e V e {u} T K e {u} p-1 K e {u} ds T .

As using such a full tangent matrix is numerically expensive, it is better to replace it by a truncated one. This can be done by neglecting the last term in (31). We have tested both methods and their performances are similar.

Comparison with an analytical solution

An analytical optimal design

In order to test our numerical procedure, we construct in this section a non-trivial analytical solution: we consider the case of an elastic material with a vanishing Poisson's coefficient. Then D ij k = 2µδ ik δ j and the constraint for the displacement field u takes the simple form:

ε u : D : ε u = 2µε u : ε u 1. ( 32
)
Using adapted units, we consider 2µ = 1. Then, in support of the optimal design {x, h(x) > 0} where the constraint is active (the inequality is actually an equality), the strain tensor has to verify:

ε 2 11 + ε 2 22 + 2ε 2 12 = 1, ( 33
)
where ε ij denote the components of ε u in an orthogonal basis (x 1 , x 2 ).

On the other hand, as ε u is a symmetric part of a gradient, it has to verify the compatibility condition:

∂ 2 ε 11 ∂x 2 2 + ∂ 2 ε 22 ∂x 2 1 -2 ∂ 2 ε 12 ∂x 1 ∂x 2 = 0. (34)
Finally, u must correspond to the equilibrium displacement field of the optimal design. In a domain free from external forces, this condition states that there exists a scalar field h(x) such that div(hε u ) = 0, i.e.:

hdiv(ε u ) + grad(h) • ε u = 0, or div(ε u ) + grad ln(h) • ε u = 0. ( 35
)
The existence of such a field h is assured, when ε u is invertible, if ε u satisfies:

rot ε -1 u • div(ε u ) = 0. ( 36
)
It is difficult to describe the general properties of the solutions of equations ( 33), ( 34) and ( 36), but it is easy to find particular non-trivial solutions.

Example 1: Let us determine, for instance, a solution ε u diagonal in the Cartesian coordinate system (x 1 , x 2 ), i.e. such that ε 12 (x) = 0 everywhere. Then equations ( 33), ( 34) and (36) become:

ε 2 11 + ε 2 22 = 1, ∂ 2 ε 11 ∂x 2 2 + ∂ 2 ε 22 ∂x 2 1 = 0, ∂ 2 ∂x 1 ∂x 2 ln ε 22 ε 11 = 0 (37)
a solution of which is:

ε 11 = -1 + x 2 x 1 4 -1/2 , ε 22 = -1 + x 1 x 2 4 -1/2 . ( 38
)
The thickness h of the optimal plate is then given, up to a multiplicative constant, by equation ( 35). We have:

h(x) = γ x -4 1 + x -4 2 . ( 39
)
Let us assume now that the plate is contained in the square := [a, b] 2 with 0 < a < b. It is easy to determine the external density of forces F which are applied on the boundary of the domain. Indeed the equilibrium condition reads F = hε u • n, where n denotes the external normal to the boundary ∂ of . For the considered domain, F is parallel to n: F = -pn. The plate is therefore subjected to the non-constant pressure field p:

p(x) = x -2 2 along the vertical edges {a} × [a, b] and {b} × [a, b], p(x) = x -2 1 along the horizontal edges [a, b] × {a} and [a, b] × {b}. ( 40
)
The design given by ( 39) is the optimal design of a plate contained in the square and subjected to the pressure field given by ( 40), the multiplicative constant γ being determined by the volume constraint. The associated strain field is given by ( 38).

Example 2: It is also possible to obtain solutions with cylindrical symmetry by assuming that the strain tensor is diagonal in the natural basis (e r , e θ ) associated with the polar coordinates (r, θ). Such solutions (for which ε rθ = 0) correspond to more realistic loads than (40). Their computation is more intricate than in the previous example and we only give the results. The parametric description of the optimal designs has the form:

r = ae t/2 cos(t) -sin(t) 1/2 , h = be -t cos(t) -1 , ( 41
)
where the constants a and b are determined by the loading. Indeed, if we consider the annulus {r 1 < r < r 2 } subjected to an internal constant pressure p 1 and an external constant pressure p 2 , the quantity:

p := hε rr = be -t , (42)
coincides with p 1 and p 2 for r = r 1 and r = r 2 respectively. This enables us to determine the constants a and b and the values t 1 and t 2 of the parameter t which correspond to r = r 1 and r = r 2 . It is remarkable that this determination is only possible if the ratio p 2 /p 1 does not exceed critical values (p 2 /p 1 ∈ [k 1 , k 2 ]) which depend on the geometrical parameter r 2 /r 1 (for instance, when

r 2 /r 1 is very large, [k 1 , k 2 ] [0.1, 2.2]
). Outside of this interval, the optimal design is not entirely made of a regular 2-D plate: part of the mass concentrates along the circle {r = r 1 }.

Numerical confrontation

In order to validate the numerical procedure presented in section 3, we compare its results with the analytical solutions. Let us consider the square domain := [a, b] 2 submitted to the boundary forces described by (40) (see figure 1). Then the optimal design is a plate the thickness of which is given by (39).

We compute the finite element solution for increasing values of the parameter p. Generally we use the increment 2 for p. We stop at the maximum value before the Newton's method diverges. In order to check the numerical scheme, we study three different regular finite element grids with quadratic interpolation (20×20, 40×40, 80×80).

Analysing the relative difference between the numerical and analytical solutions at every point of the domain (cf. figure 2a), we find that it is lower than 1%. We can conclude that the numerical scheme presented here is efficient for solving the problem (12). The scheme's convergence is the same in both formulations ( 27) or (28). For a given value of p, we initialise the Newton's scheme by the solution obtained previously for a smaller value of p. Beginning with the value p = 1 (which corresponds to a linear elastic problem) assures the convergence of the scheme. Many iterations (about 20) are needed to cross the first non-linearity (p = 3). Afterward, for p > 3, the convergence is obtained within five iterations.

Analysing the L 2 -norm of the difference between the numerical and analytical solutions, we observe a good agreement, which improves as p increases (cf. figure 2b). Even when p is small (about 10), the design obtained is close to the exact solution: it is pertinent to compute first a solution on a large mesh for a rough design, then to refine the mesh in order to obtain the best precision. Indeed, the mesh quality governs the precision and the maximal value of the parameter p we can reach.

We have also considered the annulus {r 1 < r < r 2 } subjected to the pressures given by ( 42) (when their ratio does not exceed the critical values). Again we get a very good agreement with the analytical solution (41).

Topology of optimal solutions

The goal of topological optimisation is to obtain structures, the topology of which has not been postulated a priori.

What can we say about the topology of the optimal design obtained using a fictitious material approach? Essentially, we are limited to conjectures: first, it seems reasonable, and our numerical results show, that the support of the structure (the location of non vanishing thickness) is bounded. We explored every equilibrium system of three point forces (cf. figure 3). These numerical experiments have been computed with the same relatively rough mesh and a moderate value of p (p = 10): the design is not quite precise but qualitatively pertinent. The gray scale used in the figure corresponds to the thickness of the plate. The domain in which the optimal design has been sought is much larger than the support of the obtained designs. Thus, these designs are not limited by any geometrical constraint. As a matter of fact, the optimal design seems to be included in a disk containing the support of the forces. Here circles seem to play a central role for a still unknown reason. See figure 4 in which we consider a system of vertical forces applied at three points of an horizontal line: the optimal design is a disk with variable thickness. This figure can be compared to figure 2.7, p. 106 of [START_REF] Bendsoe | Optimisation of Structural Topology, Shape and Material[END_REF] in which, in a different context and under geometrical constraints, circles also appear.

Secondly, if connected, the structure seems to be simply connected. The restriction to connected structures is essential as many examples can be found of non connected optimal structures. On the other hand, this conjecture does not hold when the design domain is subjected to some geometrical constraint. For instance, we consider in figure 5 the case of three symmetrical forces acting on the wedges of an equilateral triangle. The thickness of the optimal design is small at the center of the structure but does not vanish (figure 6 represents the thickness of the design along one of its axis of symmetry obtained for different values of the parameter p). This fact is numerically clear: we verified that this minimal thickness converges to a positive value when refining the mesh and increasing the parameter p. Moreover this phenomenon is not an artefact of our regularisation procedure: indeed, the previous comparison between numerical and analytical solutions has showed that the numerical solutions are sufficiently accurate. Finally, we notice that, outside of the optimal design, the numerical solutions tend to zero without any ambiguity.

We explored different values for the Poisson coefficient: the optimal design depends slightly on this coefficient, but there is not any qualitative difference. Note that the triangle made by three elastic 2-D bars (or any other equivalent truss structure) which is an optimal design for the Michell criterion cannot be optimal for our minimisation problem. Indeed, the optimal design we obtain has a smaller compliance. This fact has been rigorously proved in [START_REF] Bouchitté | Characterisation of optimal shapes and masses through Monge-Kantorovich equation[END_REF].

We conjecture that the structure is simply connected. Up till now, we can only prove that there is no circular hole in the optimal structure: if (h, u) is a regular solution of problem (12) in a domain containing the disk B(0, R), then it is not possible for h to vanish inside a disk B of radius R 0 and to be positive in its vicinity (for all x in {x, R 0 < x < R 1 R}). Indeed, as the field ε u is assumed to be regular, it verifies the equilibrium condition div(hε u ) = 0 in the whole disk B(0, R 1 ), and at each point of the boundary ∂B of B we have:

ε u • n = 0,
n denoting the external normal of ∂B.

The optimality condition, ε u 2 = 1 where h > 0, yields to:

τ • ε u • τ = 1 on ∂B ( 43
)
(τ denoting a unit vector tangent to ∂B) and

ε u 2 1 in B. ( 44
)
To prove that there is no field u in B verifying these two conditions, let us consider the following minimisation problem:

inf B ε u (x) 2 dx, u ∈ M , ( 45
)
where M denotes the space of those functions of the usual Sobolev space H 1 (B) such that u| ∂B ∈ H 1 (∂B) and τ • ε u • τ = 1 on ∂B. This is a linear elastic equilibrium problem. The boundary condition is not usual but it has been proved in (Bouchitté et al., 1997-b) that this minimisation problem is well-posed. Therefore, there exists a unique solution u 0 , which due to the symmetry of B, is radial. Then, it is easy to verify that u 0 (x) = x and the infimum of ( 45) is B 2 dx. The mean value of ε u 2 is larger or equal to 2 for every u in M: no element of M can verify condition (44).

Final remarks

In this study we show how shape optimisation and locking materials are closely linked. First, this enables us to have mechanic intuition of optimal designs, but overall this leads to a classical numerical approach. We only have to search for the equilibrium configuration of a strongly non-linear elastic medium. The optimisation problem we considered is very particular. Many practical applications need a more sophisticated model: one should replace the fictitious energy (2) by an energy resulting from homogenisation theory; one should optimise the stiffness under multiple loads; one should consider materials with different behaviors; one should discuss the case of mass concentration along domains of lower dimension. All these variants lead to significant difficulties. Our problem is pertinent for optimisation of 2-D structures, it raises some central questions and enables us to test some conjectures. Our numerical experiments show that, as expected, the optimal structures are bounded. To our knowledge, this quite intuitive fact is not proved even for the Michell optimisation criterion. Our experiments show also, less expectedly, that the optimal structure is simply connected. This makes a fundamental difference between Michell trusses and fictitious optimal designs. Before a penalisation step,

Figure 1 .

 1 Figure 1. Boundary conditions for the numerical test.

Figure 2 .

 2 Figure 2. Difference between analytical and numerical solutions.

Figure 3 .

 3 Figure 3. Optimal designs for some systems of three forces.

Figure 4 .

 4 Figure 4. Optimal design for a particular system of three forces.

Figure 5 .

 5 Figure 5. Optimal design for the system of three symmetrical forces.

Figure 6 .

 6 Figure 6. Thickness along the axis of symmetry.


Frédéric Golay
email: frederic.golay@univ-tln.fr

Damien Lachouette

Stéphane Bonelli
email: stephane.bonelli@aix.cemagref.fr

Pierre Seppecher
email: seppecher@imath.fr

Interfacial erosion: A three-dimensional numerical model

Keywords:

The aim of this study was to develop a numerical model to simulate the surface erosion occurring at a fluid/soil interface subject to a flow process. We used a penalization procedure to compute flow around obstacles. A fictitious domain method allowed the use of fast and efficient finite volumes approximations on Cartesian meshes and avoided unstructured body-fitted meshes. The water/soil interface evolution was described with a Level Set function. Several numerical simulations confirmed the model's ability to predict the interfacial erosion of soils.

Introduction

The present study deals with the internal erosion process subject to tangential water flow. The present approach differs from previous theories in that our description deals with a sharp fluid/soil interface. The basic balance equations and erosion constitutive law were developed by Brivois et al. [START_REF] Brivois | Soil erosion in the boundary layer flow along a slope: A theoretical study[END_REF] and Bonelli and Brivois [START_REF] Bonelli | The scaling law in the hole erosion test with a constant pressure drop[END_REF]. For the sake of simplification, sedimentation and deposition processes are neglected here, and the flow is assumed to be dilute, Lachouette et al. [START_REF] Lachouette | One-dimensional modeling of piping flow erosion[END_REF]. A penalization procedure, incorporating a fictitious domains method, Angot et al. [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF], is used to compute the flow around obstacles. The evolution of the water/soil interface is described with a Level Set function. Spatial derivatives are computed with an order 5 WENO (weighted essentially non-oscillatory) scheme, while time derivatives are computed with a fourth order (RK4) Runge-Kutta scheme. A finite volume formulation is developed on a Cartesian MAC mesh.

Governing equations

Here, we focus on the surface erosion under tangential flow of cohesive soils (like clay), considered as impervious materials. Our description deals with singular (or discontinuous) fluid/soil interfaces. As the erosion process is much slower than the fluid flow and as the eroded particles are small (10-100 μm), the fluid flow is assumed to be a dilute flow and we neglect the variation of the fluid density and redeposition. These hypotheses were validated in our previous experimental tests (such as the "Hole Erosion Test") (Bonelli and Brivois [START_REF] Bonelli | The scaling law in the hole erosion test with a constant pressure drop[END_REF]). The fluid flow is described by an incompressible Navier-Stokes system of equations:

∇.u = 0 (1) ρ ∂u ∂t + ρ(u.∇)u -∇. μ ∇u + ∇ T u + ∇ p = 0 ( 2 )
Here, u denotes the vector velocity field, ρ the fluid density, μ the viscosity and p the pressure. The mass and momentum conservation (Eqs. ( 1) and ( 2)) are valid in the fluid domain Ω f . A part of its boundary, the interface between the fluid and the soil Γ , is in movement because of the erosion process. Let us denote c Γ the celerity of the interface Γ and n its unit normal vector oriented from the soil to the fluid. At the interface the mass conservation implies a jump condition:

Jρ(c Γ -u) • nK = 0 (3)
As we consider surface erosion between an impervious soil (u = 0) and a tangential fluid flow (u • n = 0), the mass flux crossing the interface and evacuated by the dilute flow (of density ρ) is defined as follows:

ṁ = ρc Γ • n (4)
It is then necessary to describe this erosion process. For our experiments and in line with Brivois et al. [START_REF] Brivois | Soil erosion in the boundary layer flow along a slope: A theoretical study[END_REF] and others (Foster and Fell [5], Wan and Fell [START_REF] Wan | Investigation of rate of erosion of soils in embankment dams[END_REF]), we use a constitutive erosion law driven by tangential shear stress, written in the form of a threshold law as follows:

ṁ = k er (|τ | -τ c ) if |τ | > τ c 0 o t h e r w i s e ( 5
)
where τ is the tangential (shear) stress exerted on the soil, while k er and τ c are positive parameters of the erosion law.

Mathematical model

We use a spatial discretization which is not boundary-fitted to the fluid domain. The main idea consists in incorporating the original domain of study in a geometrically bigger and simpler one called "fictitious domain". This "penalized method", is popular because it allows the use of fairly structured meshes. The fictitious domain method is now widely used, for example, in fluid-porous-solid flows and has been mathematically justified (see e.g. Angot [START_REF] Angot | Analysis of singular perturbations on the Brinkman problem for fictitious domain models of viscous flows[END_REF], Angot et al. [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF], Diaz-Goano et al. [START_REF] Diaz-Goano | A fictitious domain/finite element method for particulate flows[END_REF], Khadra et al. [START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF]). Let us take a computation domain Ω (fictitious domain) with the previous fluid domain Ω f and the soil domain Ω s . The penalized method allows describing the behavior of the two sub-domains with a set of equations defined in the whole domain introducing a new term in the balance equation ( 2):

ρ ∂u ∂t + ρ(u.∇)u -∇. μ ∇u + ∇ T u = -∇ p -H μ K s u (6)
where K s is a penalization coefficient (a very small number having the dimension of a permeability) and H is the characteristic function of the soil domain, which is unity within Ω s and zero elsewhere. As the last term of the right-hand side of ( 6) is zero in the fluid domain, it corresponds to the Navier-Stokes equation. In the soil domain, penalization leads to very low values of u, and the left-hand side of ( 6) is negligible, as well as the viscous part of the stress tensor.

As we consider a sharp interface, the Level Set method is an appropriate formulation for describing the interface. It was introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature dependent speed: Algorithm for tracking material interface[END_REF] and is widely used in multi-fluid flows (e.g. Chang et al. [START_REF] Chang | A level set formulation of Eulerian interface capturing methods for incompressible fluid flows[END_REF], Olson and Kreiss [START_REF] Olson | A conservative level set method for two phase flow[END_REF], Sussman et al. [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flows[END_REF]) and the fictitious domain approach (e.g. Chantalat et al. [START_REF] Chantalat | Level-set, penalization and Cartesian meshes: A paradigm for inverse problems and optimal design[END_REF]). We introduce a function ψ , that is negative in the fluid domain Ω f and positive in the soil domain Ω s . Interface Γ is represented by the zero level set of ψ , and the value of H depends only on ψ . A relevant approach is to choose this function as the signed distance to the interface. The displacement of the interface is then driven by a transport equation:

∂ψ ∂t + c.∇ψ = 0 ( 7
)
where c is a suitable extension on the whole domain Ω of the celerity of the interface c Γ defined in (4), which is defined only on the interface. Note that the normal to the interface is given by the gradient of the Level Set function.

Numerical modeling

As the erosion process is much slower than the main fluid flow, these two phenomena are split, by using different time scales. For the fluid flow (1), ( 6), we use a finite volume method on a staggered grid for spatial discretization. This method, known as the Marker and Cell (MAC) method formulated by Harlow and Welch [START_REF] Harlow | Numerical calculation of time dependent viscous incompressible flow with free surface[END_REF], is now widely used in many applications (see, for example, Patil and Tiwari [START_REF] Patil | Computation of flow past complex geometries using MAC algorithm on body-fitted coordinates[END_REF]). The time derivative is approached with an explicit Euler scheme. In order to deal with the constraint of divergence-free velocity, we use the augmented Lagrangian method (see e.g. Vincent and Caltagirone [START_REF] Vincent | A one-cell local multigrid method for solving unsteady incompressible multiphase flows[END_REF], Galusinski and Vigneaux [START_REF] Galusinski | On stability condition for bifluid flows with surface tension: Application to microfluidics[END_REF]). This iterative predictor/corrector method consists in solving an optimization problem by solving a velocity-pressure saddle point with an Uzawa algorithm, in order to converge towards a solution which satisfies the incompressibility constraint. The linear algebra system is solved numerically with an iterative Bi-Conjugate Gradient Stabilized algorithm (BiCGSTAB). This algorithm is particularly effective for managing computer memory, because it does away with the need to build and store the linear matrix. We simply need to express the product of the matrix and the vector of the unknowns. In [START_REF] Angot | Analysis of singular perturbations on the Brinkman problem for fictitious domain models of viscous flows[END_REF], which describes the interface motion, the spatial gradients are approached with a high order Weighted Essentially Non-Oscillating (WENO5) scheme introduced by Jiang and Peng [START_REF] Jiang | Weighted ENO schemes for Hamilton-Jacobi equations[END_REF], while the time discretization uses a fourth order Runge-Kutta scheme (RK4).

The computation of the shear stress at the interface is crucial because it determines the rate of the erosion [START_REF] Foster | Assessing embankment dam filters that do not satisfy design criteria[END_REF]. The numerical estimation of the velocity gradient can be difficult because it is discontinuous across the interface. As the velocity almost vanishes in the soil, we obtain an accurate estimation of the shear stress by computing it with the fluid velocity only (on the fluid side on the interface), i.e. by extending this velocity in a few cells in the soil near the interface, in a way similar to the "ghost-fluid" method (Fedkiw [20]). To improve accuracy we propose a better approximation of the shear stress at the interface. Indeed the velocity gradient is computed at the center x of the cell and not at the interface. That is why we extrapolate the gradient (or the stress) to the interface: τ (xψn) ≈ τ (x)ψn.∇τ [START_REF] Diaz-Goano | A fictitious domain/finite element method for particulate flows[END_REF] Recall that ψ denotes the distance to the interface in [START_REF] Diaz-Goano | A fictitious domain/finite element method for particulate flows[END_REF].

The different parts of the numerical implementation have been validated. Among other tests, the flow solver has been tested on a Couette flow; the transport equation has been validated by the Zalesak's disk test; the erosion process has been successfully compared in the case of the "Hole Erosion Test" to the analytical solution developed by Bonelli and Brivois [START_REF] Bonelli | The scaling law in the hole erosion test with a constant pressure drop[END_REF].

Results

Erosion of a soil cylinder

We considered a test case: the erosion of a fixed cylinder of soil (radius 0.3 m, length 2 m, height 1 m, mesh 100 × 50) in a channel.

The fluid flow was approached by a stationary Stokes model. We imposed a pressure gradient from the left to the right and a wall condition at the top and the bottom of the computational domain (ρ f = 1000 kg/m 3 , μ f = 0.001 Pa/s, k er = 0.001 s/m, τ c = 10-20 Pa, r 1 = r 2 = 1, cfl = 0.8, δ P = 0.1 Pa).

Fig. 1 shows the evolution of the cylinder due to interfacial erosion. We then verified that the shape and flow were symmetrical, as was expected, and that the cylinder was mostly eroded at the upper and lower part of the interface. A stationary Stokes model was used because the erosion process is much slower than the flow process, as highlighted by the erosion characteristic time (Bonelli and Brivois [START_REF] Bonelli | The scaling law in the hole erosion test with a constant pressure drop[END_REF]). A parametric study showed that the flow can be assumed to be stationary until the erosion coefficient reaches a low value (k er 10 -3 ). A parametric study of the mesh underlined the consistency of the formulation.

Erosion around a bridge

In order to test several three-dimensional computations, we consider the erosion of soil around a bridge pier (Figs. 2 and3). A wall condition is applied at the bottom of the computational box, a symmetric condition at the top, a periodic condition on the sides and a pressure gradient is imposed between the inlet and the outlet. We use a coarse mesh 150 × 50 × 50 and a Stokes solver. In a realistic case, a free surface should be considered (e.g. Golay and Helluy [START_REF] Golay | Numerical schemes for low Mach wave breaking[END_REF]) and the full Navier-Stokes equation should be used. These developments are out of the scope of this study: here we only intend to check the possibility of simulating an erosion process in a 3D case.

Two Level Set functions were used: one for the interface with the bridge pier, which was a non-erodible material, and another for the soil-fluid interface. The erosion process stopped as the shear stress at the interface reached the threshold shear stress. This numerical simulation required about 2 days on a personal computer (Xeon 3 GHz), so it proved the feasibility of such computations and thus more realistic ones.

Conclusion

The aim of this study was to draw up a numerical model for simulating interfacial erosion under a tangential flow. We presented an erosion model with a unified model using a penalization procedure (fictitious domain method), thereby avoiding prohibitive remeshing methods. The water/soil interface evolution was described with a Level Set function. A finite volume formulation was developed classically on a Cartesian MAC mesh to approximate the flow, in order to implement appropriate fast solvers. The ability of the model to predict the interfacial erosion of soils was confirmed by presenting several three-dimensional simulations. In further works, we intend to improve the fluid solver, for example with turbulent models or parallel formulations, in order to achieve more realistic simulations. For turbulent flows, the challenge will be to give a sufficiently accurate description of the boundary layers in order to get correct values of the shear stress and consequently of the erosion rate.

Fig. 1 .

 1 Fig. 1. Evolution of the shape of the cylinder through time.

Fig. 2 .

 2 Fig. 2. Erosion around a bridge pier.

Fig. 3 .

 3 Fig. 3. Erosion around a bridge pier at t = 0, 1.9 × 10 7 , 4.2 × 10 7 , 5.5 × 10 7 , 6.6 × 10 7 , 9.5 × 10 7 .

Acknowledgements

This work was funded by the French National Research Agency (ANR) through the COSINUS program (CARPEINTER project n • ANR-08-COSI-002), and by the Région Provence Alpes Côte d'Azur.