Günter Leugering

Jan Sokolowski

Keywords: Signorini problem, variational inequality, shape functional, shape sensitivity, topological sensitivity, domain decomposition, Steklov-Poincaré operator, contact problems

may come

Passive control of singularities by topological optimization. The second order mixed shape derivatives of energy functionals for variational inequalities

Passive control of singularities by topological optimization. The second order mixed shape derivatives of energy functionals for variational inequalities Günter Leugering, Jan Soko lowski and Antoni Zochowski Abstract. The class of nonsmooth shape optimization problems for variational inequalities is considered. The variational inequalities model elliptic boundary value problems with the Signorini type unilateral boundary conditions. The shape functionals are given by the first order shape derivatives of the elastic energy. In such a way the singularities of weak solutions to elliptic boundary value problems can be characterized. An example in solid mechanics is given by the Griffith's functional, which is defined in the plane elasticity to measure SIF, the so-called stress intensity factor, at the crack tips. Thus, the topological optimization can be used for passive control of singularities of weak solutions to variational inequalities.

The Hadamard directional differentiability of metric projection onto the positive cone in fractional Sobolev spaces is employed to the topological sensitivity analysis of weak solutions of nonlinear elliptic boundary value problems. The first order shape derivatives of energy functionals in the direction of specific velocity fields depend on the solutions to variational inequalities in a subdomain. The domain decomposition technique is used in order to separate the unilateral boundary conditions and the energy asymptotic analysis.

The topological derivatives of nonsmooth integral shape functionals for variational inequalities are derived. The singular geometrical doamin perturbations in an elastic body Ω are approximated by the regular perturbations of bilinear forms in variational inequality, without any loss of precision for the purposes of the second order shape-topological sensitivity analysis. The second-order shapetopological directional derivatives are obtained for the Laplacian and for the linear elasticity in two and three spatial dimensions. In the proposed method of sensitivity analysis, the singular geometrical perturbations ǫ → ωǫ ⊂ Ω centered at x ∈ Ω are replaced by regular perturbations of bilinear forms supported on the manifold ΓR = {|x -x| = R} in an elastic body, with R > ǫ > 0. The obtained expressions for topological derivatives are easy to compute and therefore useful in numerical methods of topological optimization for contact problems.

Introduction

Topological derivatives of shape functionals Ω → J(Ω) are introduced in [START_REF] Soko Lowski | On the topological derivative in shape optimization[END_REF] for linear elliptic boundary value problems defined in singularly perturbed domains ǫ → Ω(ǫ), where ǫ → 0 is a small parameter which governs the size of small hole or inclusion in the bounded domain Ω ⊂ R d , d = 2, 3. The topological derivatives are given by expressions depending on pointwise values of solutions as well as of its gradients [START_REF] Novotny | Topological derivative in shape optimization[END_REF]. Therefore, the obtained expressions for topological derivatives are not well defined on the energy spaces associated with the boundary value problems under considerations.

In this paper the topological sensitivity analysis of solutions to variational inequalities is performed by the domain decomposition technique. The regular perturbations defined on the energy space ε → ε d b(Γ R ; •, •) for bilinear forms ε → a(Ω(ǫ); •, •) of boundary value problems, with respect to small parameter ǫ → 0, are introduced. Such perturbations are given by line integrals in two spatial dimensions, or by surface integrals in three spatial dimensions. As a result, the topological derivatives of shape functionals can be derived for solutions of variational inequalities posed in the intact domain Ω.

In order to derive the topological derivatives by an application of the domain decomposition technique the artificial interface Σ ⊂ Ω is introduced and Ω := Ω 1 ∪ Σ ∪ Ω 2 is decomposed into two subdomains.

For the boundary value problem under considerations such a decomposition is indeed useful. In some applied problems we are interested in the influence of singular perturbations in subdomain Ω 1 on the behaviour of solutions in subdomain Ω 2 . The functional under consideration is the elastic energy E(Ω) of whole domain Ω. The mixed second order derivatives of shape-topological or topologicalshape types for the elastic energy are evaluated. The shape sensitivity analysis is performed e.g., in Ω 2 , then the asymptotic analysis is performed in the second subdomain Ω 1 . In the framework of shape-topological sensitivity analysis the velocity method is applied in order to determine the shape functional J(Ω) := dE(Ω; V ), where V is the specific vector field in derivation of V → dE(Ω; V ). Then the asymptotic expansion of ǫ → J(Ω ǫ ) is evaluated. In the framework of topological-shape sensitivity analysis, first the asymptotic expansion of ǫ → E(Ω ǫ ) is performed, and the first order term of such an expansion is called the topological derivative. It turns out (25; 22) that the topological derivative of energy functional is not well defined for arbitrary elements from the energy space of the elasticity boundary value problems under considerations. Therefore, we introduce the equivalent representations of topological derivatives which are well defined in the energy space. These representations can be used as well to modify the state equations by replacing the singular domain perturbations by the regular perturbations of bilinear forms in variational setting.

The asymptotic expansion of the energy functional performed in one subdomain, e.g., Ω 1 , can be used in the second subdomain Ω 2 to evaluate the asymptotic expansion of the Steklov-Poincaré operator on the interface between subdomains. The method is justified by the fact that the first order expansion of the energy functional in the subdomain leads to the first order asymptotic expansion of the Dirichlet-to-Neumann mapping on the interface between subdomains. Thus, the first order expansion of the Steklov-Poincaré operator on the interface for the second subdomain is obtained. In this way the first order expansion of the energy functional in the truncated domain Ω 2 is derived. The precision of the obtained expansion is sufficient (27; 28) to replace the original energy functional by its first order expansion, provided the obtained expression is well defined on the energy space. Furthermore, the first order approximation of the energy functional in Ω is established. We point out that another method of approximation of the state equation by using the so-called self-adjoint extensions of the elliptic operators can be considered (20; 21).

Asymptotic approximation for variational inequalities

The proposed domain decomposition method is important for variational inequalities. The asymptotic analysis of solutions to variational inequalities is more involved (2) compared to the analysis of solutions to linear elliptic boundary value problems.

The variational inequality under consideration results from the minimization problem of quadratic functional

v → I(v) = 1 2 a(v, v) -L(v) (1.1)
over a convex, closed subset K ⊂ H of the Hilbert space H called the energy space. The function space H := H(Ω) is a Sobolev space which contains the functions defined over a domain Ω ⊂ R d , d = 2, 3. The singular geometrical perturbation ω ǫ centered at x ∈ Ω of the domain Ω is denoted by Ω ǫ , the size of perturbation is governed by a small parameter ǫ → 0. As an example can serve the hole or inclusion at the origin B ǫ := {|x| < ǫ}.

The quadratic functional defined on

H := H(Ω ǫ ) becomes v → I ǫ (v) = 1 2 a ǫ (v, v) -L ǫ (v) (1.2)
with the minimizers denoted by u ǫ ∈ K := K(Ω ǫ ). The expansion of associated energy functional

ǫ → E(Ω ǫ ) := I ǫ (u ǫ ) = 1 2 a ǫ (u ǫ , u ǫ ) -L ǫ (u ǫ ) (1.3)
is considered at ǫ = 0. Namely, we are looking for its asymptotic expansion

E(Ω ǫ ) = E(Ω) + ǫ d T ( x) + o(ǫ d ), (1.4)
where x → T ( x) is the topological derivative (25; 22). We show that there are regular perturbations of bilinear form defined on the energy space

H(Ω), v → b(v, v)
such that the perturbed quadratic functional defined on the unperturbed function space H(Ω)

v → I ǫ (v) = 1 2 a(v, v) + ǫ d b(v, v) -L(v) (1.5)
furnishes the first order expansion (1.4). In our applications to contact problems in linear elasticity it turns out that the bilinear form

v → b(v, v) is supported on Γ R := {|x -x| = R} ⊂ Ω with R > ǫ > 0.
Remark 1. The contact problems in elastic bodies are modeled by variational inequalities

u ∈ K : a(u, v -u) ≥ L(v -u) ∀v ∈ K. (1.6)
For the sensitivity analysis in singularly perturbed geometrical domains, the weak solutions of contact problems ǫ → u ǫ are given by perturbed variational inequalities

u ∈ K : a(u, v -u) + ǫ d b(u, v -u) ≥ L(v -u) ∀v ∈ K, (1.7)
where ǫ → 0 measures the size of singular perturbation. This is the main contribution of the paper. Therefore, we need the form of ǫ d b(u, vu) in order to apply our method of sensitivity analysis to numerical methods of topological optimization.

Variational inequalities are used to model contact problems in elasticity. It is known that the solutions to variational inequalities are Lipschitz continuous with respect to the shape [START_REF] Soko | Introduction to shape optimization, Shape sensitivity analysis[END_REF]. In general, the state governed by a variational inequality is not Fréchet differentiable with respect to the shape. For a class of variational inequalities described by the unilateral constraints in Sobolev spaces of Dirichlet type the metric projection onto the constraints turns out to be Hadamard differentiable [START_REF] Frémiot | On the analysis of boundary value problems in nonsmooth domains[END_REF]. This property is used in order to obtain the first order directional differentiability of the associated shape functionals.

In order to show the second order shape differentiability for variational inequalities, we have to restrict ourselves to energy-type shape functionals. The energy functional is the so-called marginal function and it is Fréchet differentiable with respect to the shape [START_REF] Frémiot | On the analysis of boundary value problems in nonsmooth domains[END_REF]. The first order shape derivative of the energy functional in the direction of a specific velocity vector field is considered as the shape functional for topological optimization. Thus, its topological derivative is evaluated.

The possible applications of shape-topological derivatives include the control of singularities of solutions to variational inequalities by insertion of elastic inclusions far from the singularities.

We describe the shape-topological differentiability of the energy shape functional for the Signorini problem in two spatial dimensions. The same idea can be used for the frictionless contact problems in linear elasticity.

Let us consider the Signorini problem posed in Ω ⊂ R 2 , with boundary ∂Ω = Γ ∪ Γ 0 , and

Γ c ⊂ Γ. Denote H 1 Γ0 (Ω) = {v ∈ H 1 (Ω) | v = 0 on Γ 0 ⊂ ∂Ω}.
The solution u ∈ K minimizes the quadratic functional

I(v) = 1 2 a(Ω; v, v) -(f, v) Ω over the cone K = {v ∈ H 1 Γ0 (Ω) | v ≥ 0 on Γ c ⊂ Γ ⊂ ∂Ω}.
The shape functional is the energy

E(Ω) = 1 2 a(Ω; u, u) -(f, u) Ω , where
a(Ω; u, u) = Ω ∇u • ∇udx, (f, u) Ω = Ω f udx.
We assume that Γ ∩ Γ 0 = ∅. Let Γ t 0 := T t (V )(Γ 0 ) be the boundary variations [START_REF] Soko | Introduction to shape optimization, Shape sensitivity analysis[END_REF] of the Dirichlet boundary Γ 0 .

Let us consider the decomposition of

Ω = Ω 1 ∪ Σ ∪ Ω 2 , Ω 1 ∩ Ω 2 = ∅, such that Γ 0 ⊂ ∂Ω 1 and Γ c ⊂ ∂Ω 2 .
It means that the boundary variations as well as the topological asymptotic analysis are performed in Ω 1 , and the unilateral conditions are prescribed in the second subdomain Ω 2 .

The shape derivative of the energy functional with respect to the boundary variations of Γ 0 can be written in distributed form [START_REF] Soko | Introduction to shape optimization, Shape sensitivity analysis[END_REF]

dE(Ω; V ) = Ω1 A ′ (0) • ∇u, ∇u dx
where A ′ (0) = div V I -DV -DV ⊤ , under the assumption that the velocity field V is supported in a small neighborhood of Γ 0 and that supp V ∩ supp f = ∅.

The second shape functional for the purposes of topological optimization is simply defined by

J(Ω) := Ω1 A ′ (0) • ∇u, ∇u dx. (1.8)
We are going to determine the topological derivatives of Ω → J(Ω) for insertion of small inclusions in Ω 1 far from Γ 0 . In this way we could control the possible singularities on Γ 0 by topology optimization in Ω. We consider the domain decomposition method for purposes of the shape-topological differentiability of energy shape functionals. First, the domain Ω is split into two subdomains Ω 1 , Ω 2 and the interface Σ. The differentiability with respect to small parameter of the Dirichlet-to-Neumann map which lives on the boundary Σ ⊂ ∂Ω 1 is established. This map is called the Steklov-Poincaré operator for subdomain Ω 2 .

Once, the derivative of the energy functional is given, we can proceed with the subsequent topological optimization problem. For topological optimization another decomposition Ω := Ω R ∪ Γ R ∪ Ω c is introduced. The small inclusion ω ε centered at the origin x := O is located in subdomain Ω R ⊂ Ω with the interface Γ R ⊂ ∂Ω R .

Applications of Steklov-Poincaré operators in asymptotic analysis

We analyse the precision of proposed method of approximation for variational inequalities in singularly perturbed geometrical domains. We assume for simplicity that the singular perturbation is a disc B ǫ = {|x| < ǫ}.

The Signorini variational inequality in Ω

ǫ := Ω \ B ǫ , u ǫ ∈ K(Ω ǫ ) : a(Ω ǫ ; u ǫ , v -u ǫ ) -L(Ω ǫ ; v -u ǫ ) ≥ 0 ∀v ∈ u ǫ ∈ K(Ω ǫ ), (2.1)
can be considered in the truncated domain Ω c := Ω \ B R for R > ǫ > 0, R small enough. It is assumed that the source or linear form v → L(Ω; v) := (f, v) Ωc is supported in Ω c . Hence the restriction u ǫ ∈ K(Ω c ) of u ǫ ∈ K(Ω ǫ ) to the truncated domain is given by the solution to variational inequality

u ǫ ∈ K(Ω c ) : a(Ω c ; u ǫ , v -u ǫ ) + A ǫ (u ǫ ), v -u ǫ -L(Ω c ; v -u ǫ ) ≥ 0 ∀v ∈ K(Ω ǫ ), (2.2)
where A ǫ stands for the Steklov-Poincaré operator which replaces the portion of bilinear form defined over the ring C(R, ǫ) := {R > |x| > ǫ}.

Proposition 2. Assume that the Steklov-Poincaré operator admits the one-term expansion

A ǫ (v), v = A(v), v + ǫ 2 B(v), v + o(ǫ 2 ; v, v) (2.3)
with the compact remainder o(ǫ 2 ; v, v), then we can replace in (2.2) the Steklov-Poincaré operator by its one term approximation

u ǫ ∈ K(Ω c ) : a(Ω c ; u ǫ , v -u ǫ ) + A( u ǫ ), v -u ǫ + ǫ 2 B( u ǫ ), v -u ǫ -L(Ω c ; v -u ǫ ) ≥ 0 ∀v ∈ K(Ω ǫ ), (2.4) with the estimate u ǫ -u ǫ = o(ǫ 2 ). (2.5)
Remark 3. From Proposition 2 it follows that for the shape-topological differentiability of the energy functional we can consider the variational inequality

u ǫ ∈ K(Ω) : a(Ω; u ǫ , v -u ǫ ) + ǫ 2 B( u ǫ ), v -u ǫ -L(Ω; v -u ǫ ) ≥ 0 ∀v ∈ K(Ω), (2.6)
since u ǫ -u ǫ = o(ǫ 2 ) in Ω c .
In this way, the aproximation (1.5) of quadratic functional (1.2) is justified for the first order topological derivatives of variational inequalities in truncated domains.

For the quadratic functional (1.1) and the associated boundary value problem, the bilinear form

v → b(Γ R ; v, v) := B(v), v
is determined. The linear operator B is obtained from the one term expansion of the Steklov-Poincaré operator A ǫ , the expansion results from the energy expansion in the subdomain Ω R . Therefore, the perturbed quadratic functional (1.3) can be replaced by its approximation given by (1.5). For the Signorini problem in two spatial dimensions it means that the variational inequality is obtained for minimization of perturbed functional (1.3) over the energy space in unperturbed domain Ω, and the associated energy functional

E ǫ (Ω) = 1 2 a(Ω; u ǫ , u ǫ ) + ǫ 2 2 b(Γ R ; u ǫ , u ǫ ) -(f, u ǫ ) Ω ,
is evaluated for the solution of variational inequality

u ǫ ∈ K(Ω) : a(Ω; u ǫ , v -u ǫ ) + ǫ 2 b(Γ R ; u ǫ , v -u ǫ ) -(f, v -u ǫ ) Ω ≥ 0 ∀v ∈ K(Ω).

Asymptotic analysis by domain decomposition method

In order to apply the domain decomposition technique to topological differentiability ω ǫ → J ǫ (Ω) in topologically perturbed domains Ω := Ω ǫ for the shape functionals Ω → J(Ω) we need the appropriate results on topological differentiability ǫ → B ǫ of the Steklov-Poincaré pseudodifferential boundary operators defined on the artificial interface Σ. In the particular case of holes ǫ → ω ǫ the notation is straighforward, with the singularly perturbed domain Ω ǫ := Ω \ ω ǫ and with the shape functional to be analysed with respect to small parameter ǫ → J ǫ (Ω) := J(Ω \ ω ǫ ). In the case of inclusions ǫ → ω ǫ the shape functional depends on the characteristic functions ǫ → χ ǫ of the domain perturbation ω ǫ .

For inclusions the state solution ǫ → u ǫ ∈ H(Ω) is obtained by solving boundary value problems with operator coefficients depending on small parameter ǫ → 0. In both cases the asymptotics of Steklov-Poincaré operators are obtained by asymptotic analysis of the energy functional for linear elliptic boundary value problems in subdomains Ω 2 which contains the perturbations ǫ → ω ǫ .

Let us consider the direct method of sensitivity analysis in subdomain Ω 1 which contains the contact subset Γ c ⊂ ∂Ω. This is possible due to the conical differentiability of metric projection onto the convex set K which is valid under some assumptions (eg., the convex, closed cone K is polyhedric in the Dirichlet space In the case of the Signorini problem in two spatial dimensions the direct method of asymptotic analysis for the shape functional (1.8)

H(Ω) (7)). Ω Γ R ε Γ 0 Γ c
J ǫ (Ω ǫ ) := Ω1 A ′ (0) • u ǫ , u ǫ dx
can be described as follows for the disc ω ǫ := B(ǫ) = {|x| < ǫ} located at the origin.

1. We solve the variational inequality in Ω 1 : determine u ∈ K and its coincidence set Ξ := {x ∈ Γ c : u(x) = 0}. Thus, the convex cone

S = {v ∈ H 1 Γ0 (Ω) : v ≥ 0 on Ξ a(Ω; u, v) = (f, v
) Ω } used in conical differentiability of the element u with respect to the shape can be determined. 2. The asymptotic analysis of solutions to variational inequality in singularly perturbed domain Ω(ǫ) : Ω \ B(ǫ) with respect to small parameter ǫ → 0 which governs the size of the hole B(ǫ) leads to the expansion

u ǫ = u + ǫ 2 q + o(ǫ 2 )
obtained by the domain decomposition method with the Steklov-Poincaré boundary operators, where q ∈ S : a(Ω; q, vq) + ǫ 2 Bq, vq R ≥ ∀v ∈ S.

The shape functional

J ǫ (Ω ǫ ) := Ω1 A ′ (0) • u ǫ , u ǫ dx
can be expanded in Ω 1 , the expansion is valid in the whole domain Ω,

J ǫ (Ω ǫ ) = Ω A ′ (0) • u, u dx + 2ǫ 2 Ω A ′ (0) • q, u dx + o(ǫ 2 ),
however the obtained expression for the topological derivative may not be constructive in numerical methods. We want to obtain an equivalent expression, when possible, which replaces the topological derivative

T (O) = 2 Ω A ′ (0) • q, u dx
in the first order expansion of the energy functional for Signorini problem. In the linear boundary value problems such an expression can always be obtained by the introduction of an appropriate adjoint state. We point out that for variational inequalities the existence of an adjoint state in general cannot be expected.

Asymptotic analysis of boundary value problems in rings or spherical shells 4.1. Elasticity boundary value problems

In this section we shall consider asymptotic corrections to the energy function corresponding to the elasticity system or Laplace equation in IR d , where d = 2, 3. The change of the energy is caused by creating a small ball-like void of variable radius ǫ in the interior of the domain Ω, with homogeneous Neumann boundary condition on its surface. We assume that this void has its centre at the origin O. In order to eliminate the variability of the domain, we take as

Ω R the open ball B(O, R) = B(R)
with fixed R. In this way the void B(ǫ) is surrounded by B(R) ⊂ int Ω. We denote also the ring or spherical shell as

C(R, ǫ) = B(R) \ B(ǫ), Ω(R) = Ω \ B(R) and Γ R = ∂B(R).
Using these notations we define our main tool, namely the Dirichlet-to-Neumann mapping for linear elasticity or the Steklov-Poincaré operator

A ǫ : H 1/2 (Γ R ) -→ H -1/2 (Γ R )
by means of the boundary value problem:

(1 -2ν)∆w + grad div w = 0, in C(R, ǫ), (4.1)

w = v on Γ R , σ(w).n = 0 on Γ ǫ so that A ǫ v = σ(w).n on Γ R . (4.
2) Domain decomposition -Steklov-Poincaré operator. Let u R be the restriction of u to Ω(R) and γ R ϕ the projection of ϕ on Γ R . We may then define the functional

I R ǫ (ϕ ǫ ) = 1 2 Ω(R) σ(ϕ ǫ ) : ε(ϕ ǫ ) dx - ΓN h.ϕ ǫ ds+ (4.3) + 1 2 ΓR (A ǫ γ R ϕ ǫ ).γ R ϕ ǫ ds
and the solution u R ǫ as a minimal argument for

I R ǫ (u R ǫ ) = inf ϕ ǫ ∈K⊂Vǫ I R ǫ (ϕ ǫ ), (4.4)
Here lies the essence of the domain decomposition concept: we have replaced the the variable domain by a fixed one, at the price of introducing variable boundary operator A ǫ . The above expressions have even simpler form in case of a single Laplace equation. It is enough to replace the displacement by the scalar function u, elasticity operator by -∆, and σ(u) := grad u, ε(u) := grad u, σ(u).n := ∂u/∂n.

The goal is to find the expansion

A ǫ = A + ǫ d B + R ǫ , (4.5)
where the remainder R ǫ is of order o(ǫ d ) in the operator norm in the space L(H 1/2 (Γ R ), H -1/2 (Γ R )), and the operator B is regular enough, namely it is bounded and linear:

B ∈ L L 2 (Γ R ), L 2 (Γ R ) .
Under this assumption the following propositions hold.

Proposition 4. Assume that (4.5) holds in the operator norm. Then strong convergence takes place

u R ǫ → u R (4.6)
in the norm of H 1 (Ω(R)).

Proposition 5. The energy functional has the representation

I R ǫ (u R ǫ ) = I R (u R ) + ǫ d B(u R ), u R R + o(ǫ 3 ) , (4.7)
where o(ǫ d )/ǫ d → 0 with ǫ → 0 in the same energy norm.

Here I R (u R ) denotes the functional I R ǫ on the intact domain, i.e. ǫ := 0 and A ǫ := A, applied to truncation of u.

Generally, the energy correction for both elasticity system and Laplace operator has the form 22), [START_REF] Soko Lowski | Optimality conditions for similtaneous topology and shape optimization[END_REF]. Here II is the fourth order identity tensor, and I is the second order identity tensor. This approach is important for variational inequalities since it allows us to derive the formulas for topological derivatives which are similar to the expressions obtained for the corresponding linear boundary value problems.

B(u R ), u R R = -c d e u (O
IP = 1 -ν 7 -5ν (10II - 1 -5ν 1 -2ν I ⊗ I) see (

Explicit form of the operator B for the Laplacian in two spatial dimensions.

If the function u is harmonic in a ball B(R) ⊂ IR 2 , of radius R > 0 and centre at x 0 = O, then the exact expressions for the first order derivatives of u take on the following form ( 27)

u /1 (O) = 1 πR 3 ΓR u • x 1 ds, u /2 (O) = 1 πR 3 ΓR u • x 2 ds.
Since the line integrals on Γ R are well defined for functions in L 2 (Γ R ), it follows that the operator B can be extended to the bounded operator on L 2 (Γ R ),

B ∈ L(L 2 (Γ R ) → L 2 (Γ R )).
The symmetric bilinear form for this operator, given by

Bu, v R = - 1 2πR 6 ΓR ux 1 ds ΓR vx 1 ds + ΓR ux 2 ds ΓR vx 2 ds
is continuous for all u, v ∈ L 2 (Γ R ). In fact, the bilinear form

L 2 (Γ R ) × L 2 (Γ R ) ∋ (u, v) → b(Γ R ; u, v) ∈ R
is continuous with respect to the weak convergence because of the simple structure

b(Γ R ; u, v) = l 1 (u)l 1 (v) + l 2 (u)l 2 (v) u, v ∈ L 1 (Γ R ) with two linear forms v → l i (v), i = 1, 2, l i (u) = 1 √ 2π R -3
ΓR ux i ds defined as line integrals on Γ R . This gives an additional regularity for the regular non-local perturbation B of the pseudo-differential Steklov-Poincaré boundary operator A ǫ .

4.3. Explicit form of the operator B for the Laplacian in three spatial dimensions.

Similarly as in two spatial dimensions, for harmonic functions in IR 3 it may be proved [START_REF] Soko | Modelling of topological derivatives for contact problems[END_REF] that

u /1 (O) = 3 4πR 4 S(R) ux 1 ds, u /2 (O) = 3 4πR 4 S(R) ux 2 ds, u /3 (O) = 3 4πR 4 S(R) ux 3 ds.
Using this one can easily write down the bilinear form

b(Γ R ; u, v) = Bu, v R = l 1 (u)l 1 (v) + l 2 (u)l 2 (v) + l 3 (u)l 3 (v)
where

l i (u, v) = 3 8π R -4 S(R) ux i ds.
From the computational point of view, the effort in comparison to the two dimensional case grows similarly as the difficulty of computing integrals over circle versus integrals over sphere.

Explicit form of the operator B for elasticity in two spatial dimensions.

Let us denote for the plain stress case

k = λ + µ λ + 3µ .
It has been proved in [START_REF] Soko | Modelling of topological derivatives for contact problems[END_REF] that the following exact formulae hold

ε 11 (O) + ε 22 (O) = 1 πR 3 ΓR (u 1 x 1 + u 2 x 2 ) ds, ε 11 (O) -ε 22 (O) = 1 πR 3 ΓR (1 -9k)(u 1 x 1 -u 2 x 2 ) + 12k R 2 (u 1 x 3 1 -u 2 x 3 2 ) ds, 2ε 12 (O) = 1 πR 3 ΓR (1 + 9k)(u 1 x 2 + u 2 x 1 ) - 12k R 2 (u 1 x 3 2 + u 2 x 3 1 ) ds.
These expressions are easy to compute numerically, but contain additional integrals of third powers of x i . Therefore, strains ε ij (O) may be expressed as linear combinations of integrals over circle which have the form

ΓR u i x j ds, ΓR u i x 3 j ds.
The same is true, due to Hooke's law, for stresses σ ij (O). They may then be substituted into expression for the operator B, yielding

B(u R ), v R R = - 1 2 c 2 IPσ(u) : ε(v).
These formulas are quite similar to the ones obtained for Laplace operator and easy to compute numerically.

4.5. Explicit form of the operator B for elasticity in three spatial dimensions It turns out that similar situation holds in three spatial dimensions, but obtaining the formulas is more difficult. Assuming given values of u on Γ R , the solution of elasticity system in B(R) may be expressed, following partially the derivation from ( 16) (pages 285 and later), as

u = ∞ n=0 [U n + (R 2 -r 2 )k n (ν)grad div U n ]. (4.8)
where

k n (ν) = 1/2[(3 -2ν)n -2(1 -ν)
] and r = x . In addition

U n = 1 R n [a n0 d n (x) + n m=1 (a nm c m n (x) + b nm s m n (x))]. (4.9)
The vectors

a n0 = (a 1 n0 , a 2 n0 , a 3 n0 ) ⊤ , a nm = (a 1 nm , a 2 nm , a 3 nm ) ⊤ , b nm = (b 1 nm , b 2 nm , b 3
nm ) ⊤ are constant and the set of functions

{d 0 ; d 1 , c 1 1 , s 1 1 ; d 2 , c 1 2 , s 1 2 , c 2 2 , s 2 2 ; d 3 , c 1 3 , s 1 3 , c 2 3 , s 2 3 , c 3
3 , s 3 3 ; . . .} constitutes the complete system of orthonormal harmonic polynomials on Γ R , related to Laplace spherical functions, see next paragraph. Specifically,

c l k (x) = P l,c k (x) P l,c k R , s l k (x) = P l,s k (x) P l,s k R , d k = P k ((x) Pk R .
For example,

c 2 3 (x) = 1 R 4 7 240π (15x 2 1 x 3 -15x 2 2 x 3 ), If the value of u on Γ R is assumed as given, then, denoting φ, ψ R = ΓR φψ ds,
we have for n ≥ 0, m = 1..n, i = 1, 2, 3:

a i n0 = R n u i , d n (x) R , (4.10)
a i nm = R n u i , c m n (x) R , b i nm = R n u i , s m n (x) R .
Since we are looking for ε ij (O), only the part of u which is linear in x is relevant. It contains two terms:

û = U 1 + R 2 k 3 (ν)grad div U 3 . (4.11)
For any f (x), grad div (af ) = H(f ) • a, where a -constant vector and H(f ) is the Hessian matrix of f . Therefore

û = 1 R [a 10 d 1 (x) + a 11 c 1 1 (x) + b 11 s 1 1 (x))] (4.12) + R 2 k 3 (ν) 1 R 3 H(d 3 )(x)a 30 + 3 m=1 H(c m 3 )(x)a 3m + H(s m 3 )(x)b 3m
From the above we may single out the coefficients standing at x 1 , x 2 , x 3 in u 1 , u 2 , u 3 . For example, Observe that

ε 11 (O) = 1 R 3 3 4π a 1 11 + 1 R 5 k 3 (ν) -3 7 4π a 3
ε 11 (O) + ε 22 (O) + ε 33 (O) = 1 R 3 3 4π R u 1 , c 1 1 R + R u 2 , s 1 1 R + R u 3 , d 1 R and c 1 1 = 1 R 2 3 4π x 1 , s 1 1 = 1 R 2 3 4π x 2 , d 1 = 1 R 2 3 4π
x 3 , exactly the same as for the case of Laplace equation. This should be expected, since tr ε is a harmonic function.

As a result, the operator B may be defined by the formula

Bu, u R = -c 3 IPσ(u((O))) : ε(u((O)))
but the right-hand side consists of integrals of u multiplied by first and third order polynomials in x i over Γ R resulting from (4.10). This is a very similar situation as in two spatial dimensions. Thus, the new expressions for strains make possible to rewrite B in the form possessing the desired regularity.

Laplace spherical polynomials.

For n = 1:

P1 (x) = x 3 , P 1,c 1 (x) = x 1 , P 1,s 1 (x) = x 2 , P1 R = P 1,c 1 R = P 1,s 1 R = R 2 4π 3 ,
and for n = 3:

P3 (x) = x 3 3 - 3 2 x 2 2 x 3 - 3 2 x 2 1 x 3 , P3 R = R 4 4π 7 , P 1,c 3 (x) = 6x 1 x 2 3 - 3 2 x 3 1 - 3 2 x 1 x 2 2 , P 1,c 3 R = R 4 24π 7 , P 1,s 3 (x) = 6x 2 x 2 3 - 3 2 x 3 2 - 3 2 x 2 1 x 2 , P 1,s 3 R = R 4 24π 7 , P 2,c 3 (x) = 15x 2 1 x 3 -15x 2 2 x 3 , P 2,c 3 R = R 4 240π 7 , P 2,s 3 (x) = 15x 1 x 2 x 3 , P 2,s 3 R = R 4 60π 7 , P 3,c 3 (x) = 15x 3 1 -45x 1 x 2 2 , P 3,c 3 R = R 4 1440π 7 , P 3,s 3 (x) = 45x 2 1 x 2 -15x 3 2 , P 3,s 3 R = R 4 1440π 7 ,

Asymptotic analysis of Steklov-Poincaré operators in reinforced rings in two spatial dimensions

In this section the similar asymptotic analysis of elliptic boundary value problems in subdomain Ω R ∈ IR 2 is performed, but we modify the situation, assuming that the hole is filled only partially, different material constituting a fixed part of it. In this way, we may consider double asymptotic transition, where both the size of the hole, as well as the proportion of the different material contained in it can vary. Mechanically this situation corresponds e.g. to the hole with hardened walls.

The analysis is based again on exact representation of solutions and allows to obtain the perturbation of solutions, using the fact that these solutions may be considered as minimizers of energy functional. The method is also suitable for double asymptotic expansions of solutions as well as energy form. The ultimate goal is to use obtained formulas in the evaluation of topological derivatives for elliptic boundary value problems.

Model problem

Let us consider the the domain Ω containing the hole with boundary made of modified material as depicted in Fig. 1. For simplicity the hole is located at the origin of the coordinate system. In order to write down the model problem, we introduce some notations.

B s = { x ∈ IR 2 | x < s } C s,t = { x ∈ IR 2 | s < x < t } Γ s = { x ∈ IR 2 | x = s } Ω s = Ω \ B s
Then the problem in the intact domain Ω has the form

k 1 ∆w 0 = 0 in Ω w 0 = g 0 on ∂Ω (5.1)
The model problem in the modified domain reads:

k 1 ∆w ρ = 0 in Ω ρ w ρ = g 0 on ∂Ω w ρ = v ρ on Γ ρ k 2 ∆v ρ = 0 in C λρ,ρ (5.2) k 2 ∂v ρ ∂n 2 = 0 on Γ λρ k 1 ∂w ρ ∂n 1 + k 2 ∂v ρ ∂n 2 = 0 on Γ ρ ,
where n 1 -exterior normal vector to Ω ρ , n 2 -exterior normal vector to C λρ,ρ , and 0 < λ < 1.

We want to investigate the influence of the small ring-like inclusion made of another material on the difference w ρw 0 in Ω R , where Γ R surrounds C λρ,ρ and R is fixed. We assume that ρ → 0+ and λ is considered temporarily constant.

If we define

u ρ = w ρ in Ω ρ v ρ in C λρ,ρ
then the problem (5.2) reduces to finding minimum of the energy functional

E 1 (u ρ ) = 1 2 Ωρ k 1 ∇u ρ • ∇u ρ dx + 1 2 C λρ,ρ k 2 ∇u ρ • ∇u ρ dx (5.3) for u ρ ∈ H 1 (Ω ρ ), u ρ = g 0 on ∂Ω.
This expression may be rewritten as

E 1 (u ρ ) = 1 2 ΩR k 1 ∇w ρ • ∇w ρ dx+ + 1 2 Cρ,R k 1 ∇w ρ • ∇w ρ dx+ + 1 2 C λρ,ρ k 2 ∇v ρ • ∇v ρ dx.
Using integration by parts we obtain

E 1 (u ρ ) = 1 2 ΩR k 1 ∇w ρ • ∇w ρ dx+ + 1 2 Γρ w ρ k 1 ∂w ρ ∂n 1 + v ρ k 2 ∂v ρ ∂n 2 ds+ + 1 2 ΓR k 1 w ρ ∂w ρ ∂n 3 ds,
where n 3 -exterior normal to Ω R . Hence, due to boundary and transmission condition,

E 1 (u ρ ) = 1 2 ΩR k 1 ∇w ρ • ∇w ρ dx + 1 2 ΓR k 1 w ρ ∂w ρ ∂n 3 ds (5.4)

Steklov-Poincaré operator

Observe that E 1 (w 0 ) corresponds to the problem (5.1). Therefore the main goal is to find the Steklov-Poincaré operator

A λ,ρ : w ∈ H 1/2 (Γ R ) -→ ∂w ρ ∂n 3 ∈ H -1/2 (Γ R ) (5.5)
where the normal derivative is computed from auxiliary problem

k 1 ∆w ρ = 0 in C ρ,R w ρ = w on Γ R w ρ = v ρ on Γ ρ k 2 ∆v ρ = 0 in C λρ,ρ k 2 ∂v ρ ∂n 2 = 0 on Γ λρ k 1 ∂w ρ ∂n 1 + k 2 ∂v ρ ∂n 2 = 0 on Γ ρ (5.6)
The geometry of domains of definition for functions is shown in Fig. 2. Now let us adopt the polar coordinate system around origin and assume the Fourier series form for w on Γ R .

Γ ρ Γ R Γ λρ w ρ v ρ
w = C 0 + ∞ k=1 (A k cos kϕ + B k sin kϕ) (5.7)
The general form of the solution w ρ is

w ρ = A w + B w log r + ∞ k=1 (w c k (r) cos kϕ + w s k (r) sin kϕ) , (5.8)
where

w c k (r) = A c k r k + B c k 1 r k , w s k (r) = A s k r k + B s k 1 r k . Similarly for v ρ : v ρ = A v + B v log r + ∞ k=1 (v c k (r) cos kϕ + v s k (r) sin kϕ) , (5.9)
where

v c k (r) = a c k r k + b c k 1 r k , v s k (r) = a s k r k + b s k 1 r k .
Additionally, we denote the Fourier expansion of v ρ on Γ ρ by

v ρ = c 0 + ∞ k=1 (a k cos kϕ + b k sin kϕ) (5.10)
From boundary conditions on Γ λρ it follows easily B v = 0, A v = c 0 , and then

B w = 0, A w = A v = c 0 = C 0 . There remains to find a k , b k , a c k , b c k , a s k , b s k , A c k , B c k , A s k , B s k assuming A k , B k as given.

Asymptotic expansion

In order to eliminate the above mentioned coefficients we consider first the terms at cos kϕ. From boundary and transmission conditions we have for k = 1, 2, . . .

A c k R k + B c k 1 R k = A k A c k ρ k + B c k 1 ρ k -a k = 0 a c k ρ k + b c k 1 ρ k -a k = 0 a c k (λρ) k-1 -b c k 1 (λρ) k+1 = 0 k 1 A c k ρ k-1 -k 1 B c k 1 ρ k+1 -k 2 a c k ρ k-1 + k 2 b c k 1 ρ k+1 = 0 (5.11)
This may be rewritten in the matrix form: grouping unknown parameters into a vector

p k = [A c k , B c k , a c k , b c k , a k ] ⊤ we obtain T (k 1 , k 2 , R, λ, ρ)p k = R k A k e 1
where

T =       R 2k 1 0 0 0 ρ 2k 1 0 0 -ρ k 0 0 (λρ) 2k 1 -ρ k 0 0 (λρ) 2k -1 0 k 1 ρ 2k -k 1 -k 2 ρ 2k k 2 0       (5.12)
where e 1 = [1, 0, 0, 0, 0] ⊤ . It is easy to see that

p k = p 0 k A k + ρ 2k p 1 k A k + o(ρ 2k ) (5.13)
where

p 0 k = lim ρ→0+ lim λ→0+ p k (k 1 , k 2 , R, λ, ρ) A k
and p 0 k = [1/R k , 0, 0, 0, 0] ⊤ , which corresponds to the ball B R filled completely with material k 1 . Similar reasoning may be conducted for terms containing sin kϕ.

As a result,

A λ,ρ = A 0,0 + ρ 2 A 1 λ,ρ (k 1 , k 2 , R, λ, ρ, A 1 , B 1 ) + o(ρ 2 ). (5.14)
The exact form of

A 1 λ,ρ (k 1 , k 2 , R, λ, ρ, A 1 , B 1
) is obtained from inversion of matrix T , but, what is crucial, it is linear in both A 1 and B 1 . They in turn are computed as line integrals

A 1 (w) = 1 πR 2 ΓR wx 1 ds, B 1 (w) = 1 πR 2 ΓR wx 2 ds.
As a result, for computing u ρ we may use the following energy form

E(u ρ ) = 1 2 Ω k 1 ∇u ρ • ∇u ρ dx+ + ρ 2 Q(k 1 , k 2 , R, λ, ρ, A 1 , B 1 ) + o(ρ 2 ), (5.15)
where

A 1 = A 1 (u ρ ), B 1 = B 1 (u ρ ) and Q is a quadratic function of A 1 , B 1 .
This constitutes a regular perturbation of the energy functional which allows computing perturbations of any functional depending on this solution and caused by small inclusion of the described above form.

Extension to linear elasticity

Let us consider the plane elasticity problem in the ring C R,ρ . We use polar coordinates (r, θ) with e r pointing outwards and e θ perpendicularly in the counter-clockwise direction. Then there exists an exact representation of both solutions, using the complex variable series. It has the form ( 12), ( 16),( 19)

σ rr -iσ rθ = 2ℜφ ′ -e 2iθ (zφ ′′ + ψ ′ ) σ rr + iσ θθ = 4ℜφ ′ 2µ(u r + iu θ ) = e -iθ (κφ -z φ′ -ψ).
(5.16)

The functions φ, ψ are given by complex series

φ = A log(z) + k=+∞ k=-∞ a k z k ψ = -κ Ā log(z) + k=+∞ k=-∞ b k z k .
(5.17)

Here µ -the Lame constant, ν -the Poisson ratio, κ = 3 -4ν in the plain strain case, and κ = (3ν)/(1 + ν) for plane stress.

Similarly as in the simple case described in former sections, the displacement data may be given in the form of Fourier series,

2µ(u r + iu θ ) = k=+∞ k=-∞ A k e ikθ
(5.18)

The traction-free condition on some circle means σ rr = σ rθ = 0. From (5.16),(5.17) we get for displacements the formula

2µ(u r + iu θ ) = 2κAr log(r) 1 z - Ā 1 r z+ + p=+∞ p=-∞ [κra p+1 -(1 -p)ā 1-p r -2p+1 -b-(p+1) r -2p-1 ]z p .
(5.19)

Similarly we obtain representation of tractions on some circle

σ rr -iσ rθ = 2A 1 z + (κ + 1) 1 r 2 Āz+ + p=+∞ p=-∞ (1 -p)[(1 + p)a p+1 + ā1-p r -2p + 1 r 2 b p-1 ]z p .
(5.20)

As we see, in principle it is possible to repeat the same procedure again, glueing solutions in two rings together and eliminating the intermediary Dirichlet data on the interface. The only difference lies in considerably more complicated calculations, see e.g. ( 9). This could be applied for making double asymptotic expansion, in term of both ρ and λ. However, in our case λ does not need to be small in comparison to ρ.

Asymptotic expansions of the Steklov-Poincaré operators and perturbations of bilinear forms in particular cases

The explicit form of solutions in B R allows us to conclude that for

w ρ H 1/2 (ΓR) ≤ Λ 0
the correction to the energy functional contains part proportional to ρ d and the remainder of order o(ρ d ). This in turn ( 27), ( 28) implies the possibility of representation

w ρ = w 0 + ρ 2 q + o(ρ 2 ) in H 1 (Ω R )
for both standard and contact problems, justifying computations of topological derivatives.

It is well known that the singularities of solutions to Partial Differential Equations due to the singularities of geometrical domains can be characterized by specific shape derivatives of the associated energy shape functionals [START_REF] Frémiot | On the analysis of boundary value problems in nonsmooth domains[END_REF]. Therefore, the influence of topological changes in domains on the singularities can be measured by the appropriate second-order topological derivatives of the energy functionals. It means that we evaluate the shape derivatives of the energy functional by using the velocity field method, and subsequently the second order topological derivatives of the functionals by an application of the domain decomposition method,

• the portion Γ 0 of the boundary with the homogeneous Dirichlet boundary conditions is deformed to obtain t → T t (V )(Γ 0 ) as well as t → E(Ω t ) for the energy shape functional; as a result the first order shape derivative J(Ω) := dE(Ω; V ) is obtained in the distributed form as a volume integral.

• the second order derivative of the energy functional is evaluated with respect to small paramater ε → 0, the parameter governs the size of small inclusion with the material defined by a contrast parameter γ ∈ [0, ∞).

We consider the energy shape functional Ω → E(Ω) for Signorini problems for the Laplacian as well for the frictionless contact. The shape derivative J(Ω) := dE(Ω; V ) of this functional is evaluated with respect to the boundary variations of the portion Γ 0 ⊂ ∂Ω. In another words the velocity vector field V is supported in a small neighbourhood of Γ 0 . The topological derivatives of J(Ω) are evaluated with respect to nucleation of small inclusions far from Γ 0 . domain decomposition method is applied in order to obtain the robust expressions for topological derivatives.

Directional differentiability of metric projection onto positive cone in fractional Sobolev spaces

Let us consider the subdomain Ω c := Ω \ Ω R with the contact zone Γ c in the scalar case as well as in an elastic body, see Figure 3. We recalll that the convex cone for the contact problem in elasticity with linearized nonpenetration conditions takes the form

K := {v ∈ H 1 (Ω c ) : v n ∈ K(Γ c ) ⊂ H 1/2 (Γ c )},
where K(Γ c ) is the positive cone in the fractional Sobolev space H 1/2 (Γ c ). The particular case is the space H 1/2 00 (Γ c ) for Γ c ⊂ Σ and the homogeneous Dirichlet conditions on the complement Σ \ Γ c , for the cracks.

Therefore, we establish the Hadamard differentiability (11; 18) of the metric projection in the Dirichlet space H 1/2 (Γ c ) onto its positive cone [START_REF] Frémiot | On the analysis of boundary value problems in nonsmooth domains[END_REF].

Let us consider the directional differentiability of the metric projection onto the positive cone in the fractional Sobolev spaces H 1/2 (Γ c ).

Sensitivity analysis of variational inequalities. The conical differentiability of solutions to variational inequalities for the contact problem follows from the abstract result given by Theorem 14. The general result [START_REF] Soko | Introduction to shape optimization, Shape sensitivity analysis[END_REF] is adapted here to our setting within the domain decomposition framework. Thus, the bilinear form a(•, •)+b t (•, •) defined in the subdomain Ω c is introduced, where b t (•, •) is the contribution from the Steklov-Poincaré operator on Γ R = ∂Ω R . The real parameter t > 0 governs the shape perturbations of the inclusion t → ω t in Ω R , where t → 0 governs the topological changes of Ω R in the framework of asymptotic analysis.

Two boundary value problems in two subdomains are coupled by the transmission conditions on the interface Γ R . The linear boundary value problem in Ω R furnishes the expansions of the Steklov-Poincaré operators resulting from perturbations of the inclusion in the interior of the subdomain. The sensitivity analysis of solutions to variational inequality in Ω c is performed for compact perturbations of nonlocal boundary conditions on the interface. As a result, the weak solution to the unilateral elasticity boundary value problem under considerations is directionally differentiable with respect to the parameter t → 0 which governs the perturbations of the inclusion far from the contact set. Now, we provide the precise result on the conical differentiability of solutions to variational inequalities (11; 18; 29) (see also [START_REF] Frémiot | On the analysis of boundary value problems in nonsmooth domains[END_REF]) which is given here without the proof.

Let K ⊂ H be a convex and closed subset of a Hilbert space H, and let •, • denote the duality pairing between H ′ and H, where H ′ denotes the dual of H. Let us assume that there are given symmetric bilinear forms a(•, •) + b t (•, •) : H × H → R parametrized by t ≥ 0, and the linear form f ∈ H ′ , such that Condition 13. 1. There are 0 < α ≤ M such that

|a(u, v) + b t (u, v)| ≤ M u v , α u 2 ≤ a(v, v) + b t (v, v) ∀u, v ∈ H (7.1)
uniformly with respect to t ∈ [0, t 0 ). Furthermore, there exists Q ′ ∈ L(H; H ′ ) such that

Q t = Q + tQ ′ + o(t) in L(H; H ′ ) , (7.2)
where

Q t ∈ L(H; H ′ ) a(φ, ϕ) + b t (φ, ϕ) = Q t (φ), ϕ ∀φ, ϕ ∈ H.
2. The set K ⊂ H is convex and closed, and the solution operator H ′ ∋ f → P(f ) ∈ H for (7.6)

P(f ) ∈ K : a(P(f ), ϕ -P(f )) ≥ f, ϕ -P(f ) ∀ϕ ∈ K (7.3)
is differentiable in the sense that ∀h ∈ H ′ : P(f + sh) = P(f ) + sP ′ (h) + o(s) in H (7.4) for s > 0, s small enough, where the mapping P ′ : H ′ → H is continuous and positively homogeneous, in addition, the remainder o(s) is uniform with respect to the direction h ∈ H ′ on compact subsets of H ′ .

Let us consider the unique solutions u t = P t (f ) to variational inequalities depending on a parameter t ∈ [0, t 0 ), t 0 > 0,

u t ∈ K : a(u t , ϕ -u t ) + b t (u t , ϕ -u t ) ≥ f, ϕ -u t ∀ϕ ∈ K . (7.5)
In particular, for t = 0

u ∈ K : a(u, ϕ -u) + b(u, ϕ -u) ≥ f, ϕ -u ∀ϕ ∈ K , (7.6)
with u = P(f ) a unique solution to (7.6). The mapping t → u t is strongly differentiable in the sense of Hadamard at 0 + , and its derivative is given by a unique solution of the auxiliary variational inequality [START_REF] Soko | Introduction to shape optimization, Shape sensitivity analysis[END_REF].

Theorem 14. Assume that Condition 13 is satisfied. Then the solutions to the variational inequality (7.5) are right-differentiable with respect to t at t = 0, i.e. for t > 0, t small enough,

u t = u + tu ′ + o(t) in H , (7.7)
where u ′ = P ′ (-Q ′ u) . (7.8)

Metric projection onto positive cone in H

1/2 00 (Γ c ) For boundary value problems in domains with contact conditions, the unilateral conditions are prescribed on the contact set for the normal component of the displacement field. Hence, the normal component of the displacement field belongs to the positive cone in the fractional Sobolev space H 1/2 (Γ c ). The sensitivity analysis of variational inequalities for Signorini problems was reduced in [START_REF] Soko | Introduction to shape optimization, Shape sensitivity analysis[END_REF] to the directional differentiability of the metric projection onto the positive cone in a fractional space which is the Dirichlet space. This result is further extended in [START_REF] Frémiot | On the analysis of boundary value problems in nonsmooth domains[END_REF] to some crack problem. The method is also used in the present paper, however for the purposes of sensitivity analysis of contact problems. Sensitivity analysis of the crack problem. We are going to explain how the results obtained in [START_REF] Soko | Introduction to shape optimization, Shape sensitivity analysis[END_REF] for the Signorini problem in linear elasticity can be extended to the crack problems with unilateral constraints. To this end, the abstract analysis performed in [START_REF] Frémiot | On the analysis of boundary value problems in nonsmooth domains[END_REF] for the differentiability of the metric projection onto the cone of nonnegative elements in the Dirichlet space is employed.

The framework for analysis is established in function spaces over Ω := Ω + ∪ Σ ∪ Ω -, where Σ is a C 1,1 regular curve without intersections. The regularity assumption can be weakened, if necessary.

Let Γ c ⊂ Σ be the segment {(x 1 , 0) : 0 < x 1 < 1} included in the curve Σ. We denote by n the unit normal vector field on Σ which points out of Ω + , and by n the unit normal vector field on ∂Γ c orthogonal to n . We consider deformations of the crack in the direction of the vector field V colinear with n in the neighbourhood of the crack tip

A = (1, 0) ∈ Ω c ⊂ R 2 .
In the Sobolev space defined on the cracked domain Ω c , the elements enjoy jumps over the crack which are denoted by v := v +v -, and we have the regularity property of traces v ∈ H 1/2 00 (Γ c ). In our geometry of Ω c , the Sobolev space H 1/2 00 (Γ c ) coincides with the linear subspace of H 1/2 (Σ) H 1/2 00 (Γ c ) = {ϕ ∈ H 1/2 (Σ) : ϕ = 0 q.e. on Σ \ Γ c }, where q.e. means quasi-everywhere with respect to the capacity, see e.g. [START_REF] Plotnikov | Compressible Navier-Stokes Equations. Theory and Shape Optimization[END_REF] for the definition and elementary properties of the capacity useful for the existence of optimal shapes in shape optimization problems with nonlinear PDE's constraints.

In order to investigate the properties of the metric projection in the space of admissible displacement fields onto the convex cone

K := {v ∈ H 1 (Ω c ) : v n ≥ 0},
where H 1 (Ω c ) := H 1 (Ω c ; R 2 ), we need to show that the positive convex cone

K = {ϕ ∈ H 1/2 00 (Γ c ) : ϕ ≥ 0 on Γ c }.
is polyhedric in the sense of (7; 11; 18).

We consider here the rectilinear crack Γ c in two spatial dimensions. The scalar product in

H 1/2 00 (Γ c ) := H 1/2 00 (0, 1) is defined ϕ, ψ c = Γc Γc (ϕ(x) -ϕ(y))(ψ(x) -ψ(y))
|x -y| 2 dxdy (7.9)

+ Γc ϕ(x)ψ(x) + ϕ(x)ψ(x) dist(x, ∂Γ c ) dx
Polyhedricity of the positive cone in H 1/2 00 (Γ c ). In order to show the polyhedricity of the nonnegative cone K in H := H 1/2 00 (0, 1), it is enough to check the property

ϕ + , ϕ - c ≤ 0 ∀v ∈ H 1/2 00 (0, 1)
which is straightforward, here ϕ + (x) = max {v(x), 0}. The full proof of polyhedricity in such a case is provided in [START_REF] Frémiot | On the analysis of boundary value problems in nonsmooth domains[END_REF]. It is easy to check that the polyhedricity with respect to the scalar product implies the polyhedricity with respect to a bilinear form which is equivalent to the scalar product.

Theorem 15. Let us consider the variational inequality for the metric projection of f + th ∈ H onto K

u t ∈ K : u t -f -th, v -u t ≥ 0 ∀v ∈ K,
where f, h ∈ H are given, denote by Ξ{u} = {x ∈ Γ c : u(x) = 0}. Then

u t = u + tq(h) + o(t; h) in H,
where the remainder o(t; h) is uniform on compact subsets of H, and the conical diffferential of the metric projection q := q(h) is given by the unique solution to the variational inequality q ∈ S(u) : qh, vq ≥ 0 ∀v ∈ S(u)

and the closed convex cone S(u) = {ϕ ∈ H : ϕ ≥ 0 q.e. on Ξ{u}, uf, ϕ = 0}.

Rectilinear crack in two spatial dimensions

In this section the general method of shape-topological sensitivity analysis is presented in the domain Ω := Ω c ∪ Γ R ∪ Ω R , where the first subdomain Ω c contains the rectilinear cracks Γ c and the second subdomain Ω R contains the inclusion ω.

We denote by Ω in := Ω c ∪ Γ c , the first subdomain in the elastic body without the crack. We assume that there is a regular C 1,1 -curve Σ ⊂ Ω in , without intersections, which contains the rectilinear crack Γ c := {(x 1 , 0) : 0 ≤ x 1 ≤ 1}.

To simplify the presentation, let us consider a torus Ω := T := T 2 with 2π-periodic coordinates x = (x 1 , x 2 ).

The deformations of the subdomain Ω c are defined by the vector field (x, t) → V (x, t) = (v(x, t), 0), where the

C ∞ 0 (Ω + ) function (x, t) → v(x, t) is supported in [1 -δ, 1 + δ] 2 × [-t 0 , t 0 ] ⊂ Ω + ⊂ R 2 × R and v(x, t) ≡ 1 on [1 -δ/2, 1 + δ/2] 2 × [-t 0 /2, t 0 /2].
In our notation, the real variable t ∈ R is a parameter.

It means that the vector field V deforms the reference domain Ω + c to t → T t (V )(Ω + c ) just by moving the tip of the crack X = (1, 0) → x(t) = (x 1 (t), 0) in the direction of the x 1 -axis. The mapping T t : X → x(t) is given by the system of equations dx dt (t) = V (x(t), t), x(0) = X.

The boundary value problem of linear isotropic elasticity in Ω c is defined by the variational inequality

u ∈ K : a(u, v -u) ≥ (f, v -u) ∀u ∈ K , (8.1)
where

K = {v ∈ H 1 (Ω c ) : v • n := (v + -v -) • n ≥ 0 on Γ c }, (8.2)
here v = v +v -is the jump of the displacement field over the crack Γ c . The bilinear form

a(u, v) = Ωc   µ 2 2 j,k=1 (∂ j u k + ∂ k u j )(∂ j v k + ∂ k v j ) + λdiv udiv v   dx
is associated with the operator

Lu := -µ∆u -(λ + µ)grad div u. (8.3)
The deformation tensor 2ε(u) = ∂ j u k + ∂ k u j as well as the stress tensor σ(u) = associated with the displacement field u are useful in the description of the boundary value problems in linear elasticity.

The energy functional E(Ω c ) = 1/2a(u, u) -(f, u) Ωc is twice differentiable [START_REF] Frémiot | On the analysis of boundary value problems in nonsmooth domains[END_REF] in the direction of a vector field V , for the specific choice of the field V = (v, 0). The first order shape derivative

V → dE(Ω c ; V ) = 1 t lim t→0 (E(T t (Ω c )) -(E(Ω c ))
can be interpreted as the derivative of the elastic energy with respect to the crack length, we refer the reader to [START_REF] Khludnev | Griffith's formulae for elasticity systems with unilateral conditions in domains with cracks[END_REF] for the proof, the same result for the Laplacian is given in (14; 15).

Theorem 16. We have

dE(Ω c ; V ) = 1 2 Ωc {divV • ε ij (u) -2E ij (V ; u)} σ ij (u) - Ωc div(V f i )u i . (8.4)
Now we restrict our consideration to the perturbation of the crack tip only in the direction which coincides with the crack direction. The derivative is evaluated in the framework of the velocity method (29) for a specific velocity vector field V selected in such a way that the result dE(Ω c ; V ) is independent of the field V and it depends only on the perturbation of the crack tip. That is why, this derivative is called the Griffith's functional J(Ω c ) := dE(Ω c ; V ) defined for the elastic energy in a domain with crack. We are interested in the dependence of this functional on domain perturbations far from the crack. As a result, shape and topological derivatives of the nonsmooth Griffith's shape functional are obtained with respect to the boundary variations of an inclusion.

Green formulae and Steklov-Poincaré operators

The Steklov-Poincaré operator on the interface for the domain Ω c ∪ Γ R ∪ Ω R is defined by the Green formula, first as the Dirichlet-to-Neumann map in Ω R , then it is used on the interface as nonlocal boundary operator. Therefore, we recall here the Green formula for linear elasticity operators in two and three spatial dimensions.

We start with analysis in two spatial dimensions. To simplify the presentation let us consider the reference domain without a crack in the form of the torus T := T 2 with 2π-periodic coordinates x = (x 1 , x 2 ).

For the purpose of shape-topological sensitivity analysis we assume that the elastic body without the crack is decomposed into two subdomains, Ω in and Ω R , separated from each other by the interface Γ R .

Thus, the elastic body with the crack Γ c is written as

Ω := Ω c ∪ Γ R ∪ Ω R .
The rectilinear crack Γ c ⊂ Σ ⊂ Ω in is an open set, where the ficticious interface Σ ⊂ Ω in is a closed C 1,1 -curve without intersections. In our notation Ω c = Ω in \ Γ c . The bilinear form of the linear isotropic elasticity is associated with the operator Lu := -µ∆u -(λ + µ)grad div u for given Lame coefficients µ > 0, λ ≥ 0. The displacement field u in the elastic body Ω is given by the unique solution of the variational inequality u ∈ K : a(u, vu) ≥ (f, vu) ∀u ∈ K , (8.5) where

K = {v ∈ H 1 (Ω c ) : v • n := (v + -v -) • n ≥ 0 on Γ c }, (8.6)
here v = v +v -is the jump of the displacement field over the crack Γ c . Given the unique solution u ∈ K of the variational inequality and the admissible vector field V compactly supported in Ω c , we consider the associated shape functional (8.4) is differentiable in the direction of the velocity field W compactly supported in a neighbourhood of the inclusion ω ⊂ Ω R , then the Griffith's functional (8.7) is directionally differentiable in the direction of the velocity field W . Therefore, the second order directional shape derivative dE(Ω; V, W ) of the energy functional in Ω in the direction of fields V, W is obtained.

This result can be proved by the domain decomposition technique :

• the shape differentiability of the energy functional in the subdomain Ω R implies the differentiability of the associated Steklov-Poincaré operator defined on the Lipschitz curve given by the interface Ω R ∩ Ω c with respect to the scalar parameter t → 0 which governs the boundary variations of the inclusion ω; • the expansion of the Steklov-Poincaré nonlocal boundary pseudodifferential operator obtained in the subdomain Ω R is used in the boundary conditions for the variational inequality defined in the cracked subdomain Ω c and leads to the conical differential of the solution to the unilateral problem in the subdomain; • the one term expansion of the solution to the unilateral problem is used in the Griffith's functional in order to obtain the directional derivative with respect to the boundary variations of the inclusion.

Remark 18. For the cicular inclusion ω := {x ∈ Ω R : |x -y| < r 0 }, r 0 > 0, the scalar parameter t → 0 which governs the shape perturbations of ∂ω in the direction of a field W (29) can be replaced by the parameter r → r 0 . Thus, the moving domain t → ω t is replaced by the moving domain r → {x ∈ Ω R : |x -y| < r}. In this way the shape sensitivity analysis (29) for r 0 > 0 and the topological sensitivity analysis [START_REF] Novotny | Topological derivative in shape optimization[END_REF] for r 0 = 0 + are performed in the same framework for the simple case of circular inclusion.

Shape and topological derivatives of elastic energy in two spatial dimensions for an inclusion

In the subdomain Ω c the unique weak solutions ε → u := u ε of the elasticity boundary value subproblem are given by the variational inequality u ∈ K : a(Ω c ; u, vu) + b ε (Γ R ; u, vu) ≥ (f, vu) Ωc ∀v ∈ K.

In order to differentiate the solution mapping for this variational inequality, it is required to differentiate the bilinear form ε → b ε (Γ R ; u, v), which is performed in this section.

Figure 1 .

 1 Figure 1. The domain with the hole and the surrounding circle.

  ), where c d = vol(B(1)) with B(1) being the unit ball in IR d . The energy-like density function e u (O) has the form: • In case of the Laplace operator e u (O) = 1 2 ∇u R (O) 2 for both d = 2 and d = 3, see (27). • In case of the elasticity system e u (O) = 1 2 IPσ(u R (O) : ε(u R (O), where for d = 2 and plain stress IP = 1 1ν (4II -I ⊗ I) and for d = 3

Figure 2 .

 2 Figure 2. Domains of definition for w ρ and v ρ .

  evaluated in Ω c , which is called the Griffith's functionalJ(Ω c ) := dE(Ω c ; V ) .(8.7)Let ω ⊂ Ω R be an elastic inclusion. Introduce the family of inclusions t → ω t ⊂ Ω R governed by the velocity field W compactly supported in Ω R . The elastic energy in Ω R with the inclusion ω t is denoted byω t → E t (Ω R ) := 1 2 a t (Ω R ; u, u) -(f, u) ΩR .Its shape derivative dE(Ω R ; W ) in the direction W is obtained by differentiation at t = 0 of the functiont → E t (Ω R ) := 1 2 a t (Ω R ; u, u) -(f, u) ΩR .Proposition 17. Assume that the energy shape functional in the subdomain Ω R ,

	ω → E(Ω R ) :=	1 2	a(Ω

R ; u, u) -(f, u) ΩR

Acknowledgement. This work has been supported by the DFG EC315 'Engineering of Advanced Materials' and by the ANR-12-BS0-0007 Optiform.

In order to present the results, we are going to consider a simple geometry of the contact zone Γ c . In the general setting the results can be obtained in the similar way.

Therefore, we consider the subset B = {|x| < R}, x = (x 1 , • • • , x d ) ⊂ Ω, of an elastic body Ω, with the contact set Γ c := {x = (x ′ , x d ) ∈ R d : x d = 0, |x ′ | < R/2} and Σ defined by an extension of the subset Σ := {x = (x ′ , x d ) ∈ B : x d = 0}. In such a case, the unit normal vector to the contact set n := (0, . . . , 0, 1) is constant on Γ c , and the unit tangent vector to n on the boundary ∂Γ c is n := (n 1 , . . . , n d-1 , 0). For the displacement field u = (u 1 , • • • , u d ) it follows that u n = u d , hence, the unilateral constraints for the normal component over the contact set H

Thus, the convex cone of admissible displacements for the contact problem takes the form

and our analysis of the metric projection is reduced to the positive cone in H

Remark 6. We recall that in general for a domain Ω with the boundary Γ, the Sobolev spaces H 1 (Ω) and H 1/2 (Γ) are (1; 10) the so-called Dirichlet spaces. It means that for the scalar product a(•, •), with v + := sup {v, 0} and v -:= sup {-v, 0}, the property a(v + , v -) ≤ 0 holds for all elements of the Sobolev spaces.

Remark 7. The metric projection in the Dirichlet space onto the cone of nonnegative elements is considered for the purpose of sensitivity analysis of solutions to frictionless contact problems in [START_REF] Soko | Introduction to shape optimization, Shape sensitivity analysis[END_REF]. This result is extended to the crack problem. In order to avoid unnecessary technicalities, we restrict ourselves to a model problem. Now, we consider the Hadamard differentiability of metric projection in Dirichlet space onto the cone of positive elements, and recall the result on its conical differentiability.

Consider the convex, closed cone

With the element u = P K (f ) we associate the convex cone

and denote by T K (u) the closure of C K (u) in H 1/2 (Σ). On the other hand [START_REF] Frémiot | On the analysis of boundary value problems in nonsmooth domains[END_REF] there is a nonnegative Radon measure m such that for all v ∈ H 1/2 (Σ) we have the equality

We recall the result on polyhedricity of the positive cone in a Dirichlet space [START_REF] Frémiot | On the analysis of boundary value problems in nonsmooth domains[END_REF].

Lemma 9. The convex cone

is dense in the closed, convex cone

Proof. Using the property of the Dirichlet space

Remark 10. In [START_REF] Frémiot | On the analysis of boundary value problems in nonsmooth domains[END_REF] the tangent cone T K (u) is derived for u ∈ K, in the case of the positive cone K = {v ∈ H : v ≥ 0} in the Dirichlet space H equipped with the scalar product (u, v) H . We have

The convex cone S := T K (u) ∩ m ⊥ is important for our applications. It is obtained in ( 7)

The following result on the directional differentiability of metric projection holds for polyhedric convex sets (11; 18). Lemma 11. Let K be a polyhedric cone. For t > 0, t small enough,

where

and the remainder o(t; h) is uniform on compact subsets of H 1/2 (Σ). Hence, the directional derivative of the metric projection is uniquely determined by the variational inequality q := P S (h) ∈ S : (qh, vq) 1/2,Σ ≥ 0 ∀v ∈ S.

For a contact set Γ c ⊂ Σ we introduce the following convex cones

For the variational problems with unilateral conditions for the normal component of the displacement vector field over the contact set, the convex cones K(Γ c ) and K(Σ) are employed in order to show the polyhedricity of the cone of admissible displacements.

Remark 12. The proof of Lemma 9 applies as well to the convex cone

00 (Γ c ), hence, a nonnegative distribution is a Radon measure. In addition, contraction operates (4) for the scalar product (7.9) in H 1/2 00 (Γ c ). Let us note that the scalar products in H 1/2 (Σ) and in H 1/2 00 (Γ c ) are not the same, the latter is a weighted space.

We recall an abstract result on shape sensitivity analysis of variational inequalities. 9.1. Shape and topological derivatives of the energy functional in Ω R with respect to the inclusion ω In order to evaluate the topological derivative of energy functional in isotropic elasticity, the shape sensitivity analysis is combined with the asymptotic analysis [START_REF] Novotny | Topological derivative in shape optimization[END_REF]. In this section the small parameter is denoted by ε → 0, and the circular inclusion ε → ω ε := B ε is considered.

The general shape of inclusion ε → ω ε can be considered in the same way for shape sensitivity analysis [START_REF] Soko | Introduction to shape optimization, Shape sensitivity analysis[END_REF] and the asymptotic analysis [START_REF] Novotny | Topological derivative in shape optimization[END_REF].

For the sake of simplicity, the subscript R is omitted, thus, we denote Ω := Ω R , since the inclusion is located in the subdomain Ω R . We also allow for the Neumann Γ N and Dirichlet Γ D pieces of the boundary ∂Ω := ∂Ω R , thus, ∂Ω R := Γ N ∪ Γ D ∪ Γ.

Thus, we evaluate the shape and topological derivative [START_REF] Novotny | Topological derivative in shape optimization[END_REF] of the total potential energy associated to the plane stress linear elasticity problem, considering the nucleation of a small inclusion, represented by B ε ⊂ Ω, as the topological perturbation. In this way the expansion of the Steklov-Poincaré operator on the interface Γ := Γ R is obtained. Then by Green's formula E(Ω; u) = T (u), u Γ . (9.2) In the case of an inclusion ω ε ⊂ Ω, the formula becomes E ε (Ω; u) = T ε (u), u Γ .

(9.3)

Hence, the expansion of the energy functional in Ω, on the left hand side of (9.3) with respect to the parameter ε → 0 can be used in order to determine the associated expansion of the Steklov-Poincaré operator u → T (u) on the right hand side of (9.3). Therefore, let us consider the smooth domain Ω with the boundary ∂Ω := Γ N ∪ Γ D ∪ Γ, here Γ is the interface on which the Steklov-Poincaré operator introduced in our domain decomposition method is defined.