<figure><img src="image_1.png" /></figure>

HAL
open science

# Shape optimization involving the Tresca friction law in a 2D linear elastic model

Loïc Bourdin, Fabien Caubet, Aymeric Jacob de Cordemoy

► To cite this version:

Loïc Bourdin, Fabien Caubet, Aymeric Jacob de Cordemoy. Shape optimization involving the Tresca friction law in a 2D linear elastic model. 2024. hal-04785047

HAL Id: hal-04785047

https://hal.science/hal-04785047v1

Preprint submitted on 15 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Shape optimization involving the Tresca friction law in a 2D linear elastic model

Loïc Bourdin,\* Fabien Caubet,† Aymeric Jacob de Cordemoy‡

November 14, 2024

## Abstract

The aim of this work is to analyse a shape optimization problem in a mechanical friction context. Precisely we perform a shape sensitivity analysis of a Tresca friction problem, that is, a boundary value problem involving the usual linear elasticity equations together with the (nonsmooth) Tresca friction law on a part of the boundary. We prove that the solution to the Tresca friction problem admits a directional shape derivative which moreover coincides with the solution to a boundary value problem involving tangential Signorini's unilateral conditions. Then an explicit expression of the shape gradient of the Tresca energy functional is provided (which allows us to provide numerical simulations illustrating our theoretical results). Our methodology is not based on any regularization procedure, but rather on the twice epi-differentiability of the (nonsmooth) Tresca friction functional which is analyzed thanks to a change of variables which is well-suited in the two-dimensional case. The obstruction in the higher-dimensional case is discussed.

**Keywords:** Shape optimization, shape sensitivity analysis, contact mechanics, Tresca's friction law, Signorini's unilateral conditions, variational inequalities, twice epi-differentiability.

**AMS Classification:** 49Q10, 49Q12, 49J40, 74M10, 74M15, 74P10.

## 1 Introduction

Shape optimization problems involving (nonsmooth) mathematical models from contact mechanics (including for instance Signorini's unilateral conditions, Tresca's friction law, etc.) have already been investigated in the literature (see, e.g., [0, 16, 22, 21, 23] and references therein). They can be treated by using, for example, regularization procedures (see [8, 13, 14]) or dualization procedures (see [34, Chapter 4] and [35]). In order to avoid distorting the physical meaning of the contact models or obtaining abstract results involving dual elements, we have introduced in a recent series of papers [3, 4, 10, 11, 24] a new methodology based on the notion of *twice epi-differentiability* from the nonsmooth analysis literature (see, e.g., [29, 30]). In particular, this methodology has been successfully applied in [4] in order to analyze a shape optimization problem involving the Tresca friction law. Precisely, thanks to the twice epi-differentiability of the (nonsmooth) Tresca friction functional, we proved that the solution to the corresponding Tresca friction problem admits a directional shape derivative, which moreover coincides with the solution to a

<footnote>\*Institut de recherche XLIM. UMR CNRS 7252. Université de Limoges, France. loic.bourdin@unilim.fr</footnote>

<footnote>†Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LMAP, UMR 5142, 64000 Pau, France. fabien.caubet@univ-pau.fr</footnote>

<footnote>‡Sorbonne Université, Université Paris Cité, CNRS, INRIA, Laboratoire Jacques-Louis Lions, LJLL, F-75005 Paris, France. aymeric.jacob_de_cordemoy@sorbonne-universite.fr</footnote>
boundary value problem involving Signorini's unilateral conditions, and we provided an explicit expression of the shape gradient of the associated Tresca energy functional (which allowed us to provide numerical simulations illustrating our theoretical results).

However, as indicated in its title, the above paper [4] deals only with the scalar case (which has no physical sense from the point of view of contact mechanics). Therefore the objective of the present paper is to discuss the applicability of our methodology to the elastic case (which is the natural framework in contact mechanics). Here we would like to insist on the fact that this extension to the elastic case is not a simple replica of our previous paper [4]. Indeed, in addition to the obvious and significant difficulties inherent in calculations, our methodology leads in the elastic case to a major technical obstruction (that does not appear in the scalar case). The main contribution of the present paper is to show that a well-suited change of variables allows to overcome this obstruction in the two-dimensional elastic case. However, as discussed later in this introduction, the higher-dimensional elastic case remains an open challenge.

This long introduction is divided into several paragraphs in order to highlight (as concisely as possible) the major technical obstruction that appears in the application of our methodology in the general elastic case and how it can be overcome in the two-dimensional case.

**Description of the shape optimization problem.** In the sequel we will use standard notations, terminologies and assumptions that are precised in Section 2.1. Let $d \ge 2$, $f \in H^1(\mathbb{R}^d, \mathbb{R}^d)$, $g \in H^2(\mathbb{R}^d, \mathbb{R})$ such that $g > 0$ a.e. on $\mathbb{R}^d$, and $\Omega_{\text{ref}}$ be a nonempty connected bounded open subset of $\mathbb{R}^d$ with a $C^1$-boundary $\Gamma_{\text{ref}} := \text{bd}(\Omega_{\text{ref}})$ (see Remark 2.3 for comments on this $C^1$-regularity assumption) such that $\Gamma_{\text{ref}} = \Gamma_D \cup \Gamma_{T_{\text{ref}}}$, where $\Gamma_D$ and $\Gamma_{T_{\text{ref}}}$ are two measurable (with positive measure) disjoint subsets of $\Gamma_{\text{ref}}$. In this paper we consider the shape optimization problem with volume constraint given by

$$ \underset{\substack{\Omega \in \mathcal{U}_{\text{ref}} \\ |\Omega| = |\Omega_{\text{ref}}|}}{\text{minimize}} \quad \mathcal{J}(\Omega), \qquad (1.1) $$

where the set of admissible shapes is defined by

$$ \mathcal{U}_{\text{ref}} := \left\{ \Omega \subset \mathbb{R}^d \mid \begin{array}{l} \Omega \text{ nonempty connected bounded open subset of } \mathbb{R}^d \\ \text{with a } C^1\text{-boundary } \Gamma := \text{bd}(\Omega) \text{ such that } \Gamma_D \subset \Gamma \end{array} \right\}, $$

where $\mathcal{J}: \mathcal{U}_{\text{ref}} \to \mathbb{R}$ is the Tresca energy functional defined by

$$ \mathcal{J}(\Omega) := \frac{1}{2} \int_{\Omega} \text{Ae}(u_{\Omega}) : \text{e}(u_{\Omega}) + \int_{\Gamma_{T}} g \|u_{\Omega_{\tau}}\| - \int_{\Omega} f \cdot u_{\Omega}, $$

for all $\Omega \in \mathcal{U}_{\text{ref}}$, where $u_{\Omega} \in H_D^1(\Omega, \mathbb{R}^d)$ stands for the unique weak solution to the Tresca friction problem (see, e.g., [16, Chapter 3 Section 5.2] or [11, Section 2.1.3]) given by

$$ \left\{ \begin{array}{ll} -\text{div}(\text{Ae}(u)) - f = 0 & \text{in } \Omega, \\ u = 0 & \text{on } \Gamma_D, \\ \sigma_n(u) = 0 & \text{on } \Gamma_T, \\ \|\sigma_{\tau}(u)\| \le g & \text{and } u_{\tau} \cdot \sigma_{\tau}(u) + g \|u_{\tau}\| = 0 & \text{on } \Gamma_T, \end{array} \right. \qquad (\text{TP}_{\Omega}) $$

where $\Gamma := \text{bd}(\Omega)$, $\Gamma_T := \Gamma \setminus \Gamma_D$ and

$$ H_D^1(\Omega, \mathbb{R}^d) := \{v \in H^1(\Omega, \mathbb{R}^d) \mid v = 0 \text{ a.e. on } \Gamma_D\}. $$

The tangential boundary conditions on $\Gamma_T$ in $(\text{TP}_{\Omega})$ are known as the Tresca friction law. Finally recall that the unique weak solution $u_{\Omega} \in H_D^1(\Omega, \mathbb{R}^d)$ to $(\text{TP}_{\Omega})$ is characterized by the variational
inequality

$$ \int_{\Omega} \mathrm{Ae}(u_{\Omega}) : e(v - u_{\Omega}) + \int_{\Gamma_{\mathrm{T}}} g \|v_{\tau}\| - \int_{\Gamma_{\mathrm{T}}} g \|u_{\Omega_{\tau}}\| \geq \int_{\Omega} f \cdot (v - u_{\Omega}), \quad \forall v \in \mathrm{H}_{\mathrm{D}}^{1}(\Omega, \mathbb{R}^{d}), $$

and can be expressed as $u_{\Omega} = \mathrm{prox}_{\phi_{\Omega}}(F_{\Omega})$, where $F_{\Omega} \in \mathrm{H}_{\mathrm{D}}^{1}(\Omega, \mathbb{R}^{d})$ is the unique solution to the Dirichlet-Neumann problem (see, e.g., [11, Section 2.1.1]) given by

$$ \left\{ \begin{array}{ll} -\mathrm{div}(\mathrm{Ae}(F)) - f & = \quad 0 \quad \text{in } \Omega, \\ F & = \quad 0 \quad \text{on } \Gamma_{\mathrm{D}}, \\ \mathrm{Ae}(F)\mathrm{n} & = \quad 0 \quad \text{on } \Gamma_{\mathrm{T}}, \end{array} \right. \qquad (1.2) $$

and $\mathrm{prox}_{\phi_{\Omega}} : \mathrm{H}_{\mathrm{D}}^{1}(\Omega, \mathbb{R}^{d}) \to \mathrm{H}_{\mathrm{D}}^{1}(\Omega, \mathbb{R}^{d})$ stands for the proximal operator (see Definition A.1) associated with the (convex) *Tresca friction functional* $\phi_{\Omega} : \mathrm{H}_{\mathrm{D}}^{1}(\Omega, \mathbb{R}^{d}) \to \mathbb{R}$ defined by

$$ \phi_{\Omega}: \mathrm{H}_{\mathrm{D}}^{1}(\Omega, \mathbb{R}^{d}) \to \mathbb{R} \\ v \mapsto \int_{\Gamma_{\mathrm{T}}} g \|v_{\tau}\|. $$

**Application of the classical strategy from (smooth) shape optimization literature.** To deal with the numerical treatment of the above shape optimization problem, a suitable expression of the shape gradient of $\mathcal{J}$ is required. For this purpose, we follow the classical strategy developed in (smooth) shape optimization literature (see, e.g., [7, 22]). Consider $\Omega_0 \in \mathcal{U}_{\mathrm{ref}}$ and a direction $\theta \in \mathcal{C}_{\mathrm{D}}^{2,\infty}(\mathbb{R}^d, \mathbb{R}^d)$ where

$$ \mathcal{C}_{\mathrm{D}}^{2,\infty}(\mathbb{R}^d, \mathbb{R}^d) := \{ \theta \in \mathcal{C}^2(\mathbb{R}^d, \mathbb{R}^d) \cap W^{2,\infty}(\mathbb{R}^d, \mathbb{R}^d) \mid \theta = 0 \text{ on } \Gamma_{\mathrm{D}} \}. $$

For any $t \ge 0$ sufficiently small such that $\mathrm{id} + t\theta$ is a $\mathcal{C}^2$-diffeomorphism of $\mathbb{R}^d$, where $\mathrm{id} : \mathbb{R}^d \to \mathbb{R}^d$ stands for the identity map, we denote by $\Omega_t := (\mathrm{id} + t\theta)(\Omega_0) \in \mathcal{U}_{\mathrm{ref}}$ and by $u_t := u_{\Omega_t} \in \mathrm{H}_{\mathrm{D}}^{1}(\Omega_t, \mathbb{R}^d)$ (note that $u_t$ is defined on the moving domain $\Omega_t$). To get an expression of the shape gradient of $\mathcal{J}$ at $\Omega_0$ in the direction $\theta$, defined by $\mathcal{J}'(\Omega_0)(\theta) := \lim_{t \to 0^+} \frac{\mathcal{J}(\Omega_t) - \mathcal{J}(\Omega_0)}{t}$ (if it exists), the usual first step consists in introducing $\bar{u}_t := u_t \circ (\mathrm{id} + t\theta) \in \mathrm{H}_{\mathrm{D}}^{1}(\Omega_0, \mathbb{R}^d)^t$ (note that $\bar{u}_t$ is defined on the fixed domain $\Omega_0$) and obtaining an expression of the derivative (if it exists) of the map $t \in \mathbb{R}_+ \mapsto \bar{u}_t \in \mathrm{H}_{\mathrm{D}}^{1}(\Omega_0, \mathbb{R}^d)$ at $t = 0$, denoted by $\bar{u}'_0 \in \mathrm{H}_{\mathrm{D}}^{1}(\Omega_0, \mathbb{R}^d)$ and called *directional material derivative*. Then the *directional shape derivative* is defined by $u'_0 := \bar{u}'_0 - \nabla u_0 \theta$ which corresponds (roughly speaking) to the derivative of the map $t \in \mathbb{R}_+ \mapsto u_t \in \mathrm{H}_{\mathrm{D}}^{1}(\Omega_t, \mathbb{R}^d)$ at $t = 0$. We refer to Remark 2.12 for a short discussion on the terminology *directional* that has been added with respect to the classical literature on shape optimization.

To get an expression of the directional material derivative, we use the change of variables $\mathrm{id} + t\theta$ and the equality

$$ n_t \circ (\mathrm{id} + t\theta) = \frac{(\mathrm{I} + t\nabla\theta^\top)^{-1}n_0}{\|(\mathrm{I} + t\nabla\theta^\top)^{-1}n_0\|}, $$

(see, e.g., [34, Chapter 2 Proposition 2.48]) in order to prove that $\bar{u}_t \in \mathrm{H}_{\mathrm{D}}^{1}(\Omega_0, \mathbb{R}^d)$ is the unique solution to the parameterized variational inequality

$$ \int_{\Omega_0} \mathrm{J}_t \mathrm{A} \left[ \nabla \bar{u}_t (\mathrm{I} + t \nabla \theta)^{-1} \right] : \nabla (v - \bar{u}_t) (\mathrm{I} + t \nabla \theta)^{-1} \\ + \int_{\Gamma_{\mathrm{T}_0}} g_t \mathrm{J}_{\mathrm{T}_t} \left\| v - \left( v \cdot \frac{(\mathrm{I} + t \nabla \theta^\top)^{-1} n_0}{\|(\mathrm{I} + t \nabla \theta^\top)^{-1} n_0\|^2} \right) (\mathrm{I} + t \nabla \theta^\top)^{-1} n_0 \right\| $$
$$ -\int_{\Gamma_{T_0}} g_t J_{T_t} \left\| \bar{u}_t - \left( \bar{u}_t \cdot \frac{(I + t \nabla \theta^\top)^{-1} n_0}{\|(I + t \nabla \theta^\top)^{-1} n_0\|^2} \right) (I + t \nabla \theta^\top)^{-1} n_0 \right\| \\ \geq \int_{\Omega_0} f_t J_t \cdot (v - \bar{u}_t), \quad \forall v \in H_D^1(\Omega_0, \mathbb{R}^d), \quad (1.3) $$

where $f_t := f \circ (\mathrm{id} + t\theta) \in H^1(\mathbb{R}^d, \mathbb{R}^d)$, $g_t := g \circ (\mathrm{id} + t\theta) \in H^2(\mathbb{R}^d, \mathbb{R})$, $J_t := \det(I + t\nabla\theta) \in L^\infty(\mathbb{R}^d, \mathbb{R})$ is the Jacobian determinant, $J_{T_t} := \det(I + t\nabla\theta) \|(I + t\nabla\theta^\top)^{-1} n_0\| \in C^0(\Gamma_0, \mathbb{R})$ is the tangential Jacobian and I is the identity matrix of $\mathbb{R}^{d \times d}$. Thus we get that

$$ \bar{u}_t = \operatorname{prox}_{\bar{\phi}_t}(\bar{F}_t), \quad (1.4) $$

where $\bar{F}_t \in H_D^1(\Omega_0, \mathbb{R}^d)$ is the unique solution to the parameterized variational equality

$$ \int_{\Omega_0} J_t A \left[ \nabla \bar{F}_t (I + t \nabla \theta)^{-1} \right] : \nabla v (I + t \nabla \theta)^{-1} = \int_{\Omega_0} f_t J_t \cdot v, \quad \forall v \in H_D^1(\Omega_0, \mathbb{R}^d), $$

and $\operatorname{prox}_{\bar{\phi}_t}: H_D^1(\Omega_0, \mathbb{R}^d) \to H_D^1(\Omega_0, \mathbb{R}^d)$ is the proximal operator associated with the parameterized convex functional $\bar{\phi}_t: H_D^1(\Omega_0, \mathbb{R}^d) \to \mathbb{R}$ defined by

$$ \bar{\phi}_t: H_D^1(\Omega_0, \mathbb{R}^d) \longrightarrow \mathbb{R} \qquad (1.5) \\ v \longmapsto \int_{\Gamma_{T_0}} g_t J_{T_t} \left\| v - \left( v \cdot \frac{(I + t \nabla \theta^\top)^{-1} n_0}{\|(I + t \nabla \theta^\top)^{-1} n_0\|^2} \right) (I + t \nabla \theta^\top)^{-1} n_0 \right\|. $$

**Application of our methodology and facing a major obstruction.** Now the difficulty (that does not appear in standard smooth shape optimization problems) is that, from (1.4), the differentiability of the map $t \in \mathbb{R}_+ \mapsto \bar{u}_t \in H_D^1(\Omega_0, \mathbb{R}^d)$ at $t = 0$ is related to the differentiability (in a generalized sense) of the parameterized proximal operator $\operatorname{prox}_{\bar{\phi}_t}$. For this purpose, the methodology that we have developed in [3, 4, 10, 11, 24] invokes the notion of twice epi-differentiability for convex functions (introduced by Rockafellar in [29]) which ensures the *proto-differentiability* of the corresponding proximal operators. Actually, since the work by Rockafellar deals only with nonparameterized convex functions, we used instead the recent work [2] in which the notion of twice epi-differentiability has been extended to parameterized convex functions. The content of Proposition A.5 in Appendix A (extracted from [2, Theorem 4.15]) invites us to analyze the twice epi-differentiability of the parameterized convex function $\bar{\phi}_t$ in order to obtain from (1.4) the differentiability of the map $t \in \mathbb{R}_+ \mapsto \bar{u}_t \in H_D^1(\Omega_0, \mathbb{R}^d)$ at $t = 0$ and a characterization of the directional material derivative $\bar{u}'_0$.

However, at this step, a major technical obstruction appears in the elastic case (that does not appear in the scalar case in our previous paper [4]). Indeed the twice epi-differentiability of the parameterized convex functional $\bar{\phi}_t$ is naturally related to the twice epi-differentiability of the parameterized convex integrand that appears in (1.5). As a reminder, in the scalar case [4], the parameterized convex functional $\bar{\phi}_t$ is given by the (simpler) expression

$$ \bar{\phi}_t: H_D^1(\Omega_0, \mathbb{R}) \longrightarrow \mathbb{R} \qquad [\text{scalar case}] \\ v \longmapsto \int_{\Gamma_{T_0}} g_t J_{T_t} |v|, $$

and the twice epi-differentiability of the (simple) parameterized convex integrand can be analyzed. On the contrary, in the elastic case, the heavy expression of the parameterized convex integrand in (1.5) does not allow, to our best knowledge, a tractable analysis of its twice epi-differentiability.
At this step of our researches, we arrive to the conclusion that, even if our methodology based on the notion of twice epi-differentiability allows to analyze a shape optimization problem involving the Tresca friction law in the scalar case (see our previous paper [4]), we are not able, at least for now, to pursue our methodology in the general elastic case.

**Overcoming the major obstruction in the two-dimensional case d = 2.** In the two-dimensional case d = 2, one can fix $\tau_0 \in C^0(\Gamma_0, \mathbb{R}^2)$ an oriented (with an orientation arbitrarily fixed) orthonormal vector to $n_0 \in C^0(\Gamma_0, \mathbb{R}^2)$ and get that $(I + t\nabla\theta)\tau_0 \cdot (I + t\nabla\theta^\top)^{-1}n_0 = 0$ on $\Gamma_0$. Therefore Inequality (1.3) can be rewritten as

$$ \int_{\Omega_0} J_t A \left[ \nabla \bar{u}_t (I + t \nabla \theta)^{-1} \right] : \nabla (v - \bar{u}_t) (I + t \nabla \theta)^{-1} + \int_{\Gamma_{T_0}} \frac{g_t J_{T_t}}{\|(I + t \nabla \theta) \tau_0\|} |v \cdot (I + t \nabla \theta) \tau_0| \\ - \int_{\Gamma_{T_0}} \frac{g_t J_{T_t}}{\|(I + t \nabla \theta) \tau_0\|} |\bar{u}_t \cdot (I + t \nabla \theta) \tau_0| \geq \int_{\Omega_0} f_t J_t \cdot (v - \bar{u}_t), \quad \forall v \in H_D^1(\Omega_0, \mathbb{R}^2). $$

Therefore we introduce $\bar{u}_t := (I + t\nabla\theta^\top)\bar{u}_t \in H_D^1(\Omega_0, \mathbb{R}^2)$ which satisfies

$$ \int_{\Omega_0} J_t A \left[ \nabla \left( (I + t \nabla \theta^\top)^{-1} \bar{u}_t \right) (I + t \nabla \theta)^{-1} \right] : \nabla \left( (I + t \nabla \theta^\top)^{-1} (v - \bar{u}_t) \right) (I + t \nabla \theta)^{-1} \\ + \int_{\Gamma_{T_0}} \frac{g_t J_{T_t}}{\|(I + t \nabla \theta) \tau_0\|} |v \cdot \tau_0| - \int_{\Gamma_{T_0}} \frac{g_t J_{T_t}}{\|(I + t \nabla \theta) \tau_0\|} |\bar{u}_t \cdot \tau_0| \\ \geq \int_{\Omega_0} (I + t \nabla \theta)^{-1} f_t J_t \cdot (v - \bar{u}_t), \quad \forall v \in H_D^1(\Omega_0, \mathbb{R}^2), \quad (1.6) $$

and thus can be expressed as $\bar{u}_t = \text{prox}_{\bar{\phi}_t}(\bar{F}_t)$ where $\bar{F}_t \in H_D^1(\Omega_0, \mathbb{R}^2)$ is the unique solution to the parameterized variational equality

$$ \int_{\Omega_0} J_t A \left[ \nabla \left( (I + t \nabla \theta^\top)^{-1} \bar{F}_t \right) (I + t \nabla \theta)^{-1} \right] : \nabla \left( (I + t \nabla \theta^\top)^{-1} v \right) (I + t \nabla \theta)^{-1} = \\ \int_{\Omega_0} (I + t \nabla \theta)^{-1} f_t J_t, \quad \forall v \in H_D^1(\Omega_0, \mathbb{R}^2), \quad (1.7) $$

and $\text{prox}_{\bar{\phi}_t}: H_D^1(\Omega_0, \mathbb{R}^2) \to H_D^1(\Omega_0, \mathbb{R}^2)$ is the proximal operator associated with the parameterized convex functional $\bar{\phi}_t: H_D^1(\Omega_0, \mathbb{R}^2) \to \mathbb{R}$ defined by

$$ \bar{\phi}_t: H_D^1(\Omega_0, \mathbb{R}^2) \to \mathbb{R} \\ v \mapsto \int_{\Gamma_{T_0}} \frac{g_t J_{T_t}}{\|(I + t \nabla \theta) \tau_0\|} |v \cdot \tau_0|. \quad (1.8) $$

As we will see in the present paper (see Section 2.2), the parameterized convex integrand in the expression (1.8) (simpler than the one in (1.5)) allows a tractable analysis of its twice epi-differentiability, and therefore allows to continue our methodology (but only in the two-dimensional case d = 2). Precisely, thanks to the twice epi-differentiability of the parameterized convex functional $\bar{\phi}_t$, we are able to obtain from Proposition A.5 a characterization of the derivative of the map $t \in \mathbb{R}_+ \mapsto \bar{u}_t \in H_D^1(\Omega_0, \mathbb{R}^d)$ at $t = 0$, denoted by $\bar{u}'_0 \in H_D^1(\Omega_0, \mathbb{R}^d)$, and then to deduce successively a characterization of the directional material derivative given by $\bar{u}'_0 = \bar{u}'_0 - \nabla\theta^\top u_0 \in H_D^1(\Omega_0, \mathbb{R}^2)$, then a characterization of the directional shape derivative given by $u'_0 := u'_0 - \nabla u_0 \theta$, and finally an expression of the shape gradient $J'(\Omega_0)(\theta)$. These results are summarized in the next paragraph.
**Remark 1.1.** As described above, the changes of variables used in this paper lead to $\bar{u}_t = \text{prox}_{\phi_t}(\bar{F}_t)$ where the proximal operator is defined on the Hilbert space $H^1(\Omega_0, \mathbb{R}^2)$ endowed with the parameterized scalar product given by

$$ (v_1, v_2) \in (H_D^1(\Omega_0, \mathbb{R}^2))^2 \mapsto \\ \int_{\Omega_0} J_t A \left[ \nabla \left( (I + t \nabla \theta^\top)^{-1} v_1 \right) (I + t \nabla \theta)^{-1} \right] : \nabla \left( (I + t \nabla \theta^\top)^{-1} v_2 \right) (I + t \nabla \theta)^{-1} \in \mathbb{R}, $$

thus Proposition A.5 cannot be applied. This difficulty can be overcome by adding the t-independent scalar product $\langle \cdot, \cdot \rangle_{H^1(\Omega_0, \mathbb{R}^2)}$ (see (2.1)) to both members of Inequality (1.6) which leads us to replace $\bar{F}_t$ by another solution that satisfies a more complex variational equality than Equality (1.7). Actually, this difficulty also appears in the three-dimensional case $d = 3$ which can be overcome in the same manner, and also (in an easier way) in the scalar case in our previous paper [4].

**Main results in the two-dimensional case $d = 2$.** We summarize here our main theoretical results (given in Theorems 2.8 and 2.13). However we present the directional material and shape derivatives, and the shape gradient of $J$, under some additional regularity assumptions, precisely in the framework of Corollaries 2.9, 2.11 and 2.14, because their expressions are more elegant in that case. Furthermore, to ease the notations, we will use the notations $n := n_0$ and $\tau := \tau_0$.

(i) Under some appropriate assumptions described in Corollary 2.9, the map $t \in \mathbb{R}_+ \mapsto \bar{u}_t \in H_D^1(\Omega_0, \mathbb{R}^2)$ is differentiable at $t = 0$, and the directional material derivative $\bar{u}'_0 \in H_D^1(\Omega_0, \mathbb{R}^2)$ is the unique weak solution to the tangential Signorini problem (see, e.g., [11, Section 2.1.2]) given by

$$ \left\{ \begin{array}{l} -\text{div}(\text{Ae}(\bar{u}'_0)) + \text{div}(\text{Ae}(\nabla u_0 \theta)) = 0 \quad \text{in } \Omega_0, \\ \bar{u}'_0 = 0 \quad \text{on } \Gamma_D, \\ \sigma_n(\bar{u}'_0) - \xi^m(\theta)_n = 0 \quad \text{on } \Gamma_{T_0}, \\ \sigma_\tau(\bar{u}'_0) + p(\theta) \frac{u_{0\tau}}{\|u_{0\tau}\|} - \xi^m(\theta)_\tau = 0 \quad \text{on } \Gamma_{T_0^{u_{0,g}}}, \\ \bar{u}'_{0\tau} + (\nabla \theta^\top u_0)_\tau = 0 \quad \text{on } \Gamma_{T_0^{u_{0,g}}}, \\ (\bar{u}'_{0\tau} + (\nabla \theta^\top u_0)_\tau) \in \mathbb{R}_- \frac{\sigma_\tau(u_0)}{g} \\ \left(\sigma_\tau(\bar{u}'_0) - p(\theta) \frac{\sigma_\tau(u_0)}{g} - \xi^m(\theta)_\tau\right) \cdot \frac{\sigma_\tau(u_0)}{g} \le 0 \\ (\bar{u}'_{0\tau} + (\nabla \theta^\top u_0)_\tau) \cdot \left(\sigma_\tau(\bar{u}'_0) - p(\theta) \frac{\sigma_\tau(u_0)}{g} - \xi^m(\theta)_\tau\right) = 0 \quad \text{on } \Gamma_{T_0^{u_{0,g}}}. \end{array} \right. $$

where $\xi^m(\theta) := ((\text{Ae}(u_0))\nabla\theta^\top + \text{A}(\nabla u_0 \nabla\theta) + (\nabla\theta - \text{div}(\theta)\text{I})\text{Ae}(u_0)) \in L^2(\Gamma_{T_0}, \mathbb{R}^2)$ and $p(\theta) := \nabla g \cdot \theta + g(\text{div}_\tau(\theta) - \nabla\theta\tau \cdot \tau) \in L^2(\Gamma_{T_0})$, and where $\Gamma_{T_0}$ is decomposed, up to a null set, as $\Gamma_{T_0^{u_{0,g}}} \cup \Gamma_{T_0^{u_{0,g}}} \cup \Gamma_{T_0^{u_{0,g}}} (\text{see details in Theorem 2.8}).$ We emphasize the notable fact that the boundary conditions which appear on $\Gamma_{T_0^{u_{0,g}}}$ are called tangential Signorini's unilateral conditions because they are very close to the classical Signorini unilateral conditions (see, e.g., [32, 33]) except that, here, they are concerned with the tangential components (instead of the normal components in the classical case).

(ii) We deduce in Corollary 2.11 that, under appropriate assumptions, the directional shape derivative, defined by $u'_0 := \bar{u}'_0 - \nabla u_0 \theta \in H_D^1(\Omega_0, \mathbb{R}^2)$, is the unique weak solution to the
tangential Signorini problem given by

$$
\left\{
\begin{array}{l}
-\operatorname{div}(\operatorname{Ae}(u'_0)) = 0 \quad \text{in } \Omega_0, \\
u'_0 = 0 \quad \text{on } \Gamma_D, \\
\sigma_n(u'_0) - \xi^s(\theta)_n = 0 \quad \text{on } \Gamma_{T_0}, \\
\sigma_\tau(u'_0) + p(\theta) \frac{u_{0\tau}}{\|u_{0\tau}\|} - \xi^s(\theta)_\tau = 0 \quad \text{on } \Gamma_{T_{0N}^{u_0,g}}, \\
u'_{0\tau} - W(\theta)_\tau = 0 \quad \text{on } \Gamma_{T_{0D}^{u_0,g}}, \\
(u'_{0\tau} - W(\theta)_\tau) \in \mathbb{R}_- \frac{\sigma_\tau(u_0)}{g} \\
\text{and } \left(\sigma_\tau(\bar{u}'_0) - p(\theta) \frac{\sigma_\tau(u_0)}{g} - \xi^s(\theta)_\tau\right) \cdot \frac{\sigma_\tau(u_0)}{g} \le 0 \\
\text{and } (u'_{0\tau} - W(\theta)_\tau) \cdot \left(\sigma_\tau(\bar{u}'_0) - p(\theta) \frac{\sigma_\tau(u_0)}{g} - \xi^s(\theta)_\tau\right) = 0 \quad \text{on } \Gamma_{T_{0S}^{u_0,g}}.
\end{array}
\right.
$$

where $W(\theta) := -\nabla\theta^\top u_0 - \nabla u_0\theta \in H^1(\Omega_0, \mathbb{R}^2)$ and

$$
\begin{align*}
\xi^s(\theta) := \theta \cdot n (\partial_n (\mathrm{Ae}(u_0)n) - \partial_n (\mathrm{Ae}(u_0))n) &+ \mathrm{Ae}(u_0) \nabla_{\tau}(\theta \cdot n) - \nabla(\mathrm{Ae}(u_0)n)\theta \\
&+ (\nabla\theta - \mathrm{div}_{\tau}(\theta)\mathrm{I}) \mathrm{Ae}(u_0)n \in L^2(\Gamma_{T_0}, \mathbb{R}^2),
\end{align*}
$$

(iii) Finally the two previous items are used to obtain Corollary 2.14 asserting that, under appropriate assumptions, the shape gradient of $\mathcal{J}$ at $\Omega_0$ in the direction $\theta$ is given by

$$
\mathcal{J}'(\Omega_0)(\theta) = \int_{\Gamma_{T_0}} \theta \cdot n \left( \frac{\mathrm{Ae}(u_0) : \mathrm{e}(u_0)}{2} - f \cdot u_0 - \sigma_{\tau}(u_0) \cdot \partial_n(u_0) + \|u_{0\tau}\| (Hg + \partial_n g) \right) \\
+ \int_{\Gamma_{T_0}} u_{0n} \sigma_{\tau}(u_0) \cdot \tau (\nabla_{\tau}\theta_{\tau} - \nabla\theta_{\tau}) \cdot n,
$$

where $H$ stands for the mean curvature of $\Gamma_0$. One can notice that $\mathcal{J}'(\Omega_0)$ depends only on $u_0$ (and not on $u'_0$). Hence its expression is explicit and linear with respect to the direction $\theta$ and allows us to exhibit a descent direction for $\mathcal{J}$ at $\Omega_0$ (see Section 3 for details). Then, using this descent direction together with a basic Uzawa algorithm to take into account the volume constraint, we perform in Section 3 numerical simulations to solve the shape optimization problem (1.1) on a toy example.

**Obstruction in the higher-dimensional case $d \ge 3$ and additional comments.** In the higher-dimensional case $d \ge 3$, the parameterized convex functional is given by (1.5) and, to the best of our knowledge, we did not find any change of variables in order to simplify the expression of its integrand. Therefore, for almost all $s \in \Gamma_{T_0}$, we would have to investigate the twice epi-differentiability of the parameterized convex map

$$
x \in \mathbb{R}^d \mapsto g_t J_{T_t} \left\| x - \left( x \cdot \frac{(I + t \nabla \theta^\top)^{-1} n}{\|(I + t \nabla \theta^\top)^{-1} n\|^2} \right) (I + t \nabla \theta^\top)^{-1} n \right\| \in \mathbb{R}_+,
$$

that we did not succeed to prove and to compute. This is an highly nontrivial work and an interesting topic for further researches.

**Remark 1.2.** During our bibliographical researches, we discovered a surprising hypothesis made in the paper [35] that is concerned, as the present work, with a shape optimization problem in a two dimensional case, involving the Tresca friction law in a linear elastic model, but with a different
methodology based on dualization. With our notations and framework, this hypothesis consists in assuming that, for sufficiently small $t \ge 0$, it holds that

$$ (I + t\nabla\theta)n_0 = \frac{(I + t\nabla\theta^\top)^{-1}n_0}{\|(I + t\nabla\theta^\top)^{-1}n_0\|^2}. $$

With this hypothesis, we have observed that our methodology can be continued (even in the elastic case) by applying a well-suited change of variables in the expression of $\bar{\phi}_t$ in order to get the simpler expression $v \in H_D^1(\Omega_0, \mathbb{R}^2) \to \int_{\Gamma_{T_0}} g_t J_{T_t} \|(I + t\nabla\theta)\tau_0\| |v \cdot \tau_0|$, that allows a tractable analysis of the twice epi-differentiability of its integrand but which differs from (1.8). In this regard, we refer to our paper [24] (concerned with a shape optimization problem involving Signorini's unilateral conditions in the elastic case) in which the same change of variables has been applied (but without assuming the above hypothesis which is useless in the context of [24]). However we emphasize that the above hypothesis is not satisfactory. Indeed, one can easily construct numerous counterexamples for which the above equality is not satisfied. For example, in the two-dimensional case $d = 2$, one can easily construct a situation where

$$ \nabla\theta = \begin{pmatrix} 1 & 1 \\ -1 & 0 \end{pmatrix} \quad \text{and} \quad n_0 = \begin{pmatrix} 1 \\ 0 \end{pmatrix}, $$

for which the above equality is not true. We refer to [25, Remark A.2.2. p.198] for additional comments on this hypothesis and what its implies on the direction $\theta$.

**Remark 1.3.** We mention that our methodology has already been successfully applied in the elastic case (in any dimension) in our previous paper [11] (with also the emergence of tangential Signorini's unilateral conditions), but in order to solve an optimal control problem. To be clear, we underline that the obstruction encountered in the present paper does not appear in the (simpler) context of [11].

**Organization of the paper.** The paper is organized as follows. Section 2 is the core of the present work, where the main results are stated and proved. In Section 3, numerical simulations are performed to solve the shape optimization problem (1.1) on a toy example. Finally, Appendices A and B are dedicated to recalls on the notion of twice epi-differentiability and on differential geometry respectively.

## 2 Main results in the two-dimensional case $d = 2$

The whole paper is now dedicated to the two-dimensional case $d = 2$. Furthermore, all along the present section, variational equalities and variational inequalities will be involved, as well as boundary value problems (Dirichlet-Neumann problem, tangential Signorini problem, Tresca friction problem). We refer to [11, Section 2.1] for notions of strong/weak solutions, existence/uniqueness results and proximal expressions of the solutions.

This section is organized as follows. In Section 2.1, we precise the notations, terminologies and assumptions used in our linear elastic model involving the Tresca friction law. Section 2.2 is the technical part of the present paper. There, the notion of twice epi-differentiability and Proposition A.5 are used in order to characterize the derivative $\bar{u}'_0$ as the unique solution to a variational inequality (see Proposition 2.5). From that result, we deduce in Section 2.3 a characterization of the directional material derivative $\bar{u}'_0$ as the unique solution to a variational inequality (see Theorem 2.8). Then, under additional regularity assumptions, we characterize the directional material derivative $\bar{u}'_0$ and the directional shape derivative $u'_0$ as the unique weak solutions to tangential
Signorini problems (see Corollaries 2.9 and 2.11). Finally, in Section 2.4, we provide an expression of the shape gradient $J'(\Omega)(\theta)$ of the Tresca energy functional $J$ (see Theorem 2.13 and Corollary 2.14).

## 2.1 Precisions on the notations, terminologies and assumptions used in our model

Consider the shape optimization problem (1.1) in the two-dimensional case $d = 2$ and let us give some precisions on the notations, terminologies and assumptions used in our model.

First, the notation $\cdot$ stands for the standard inner product on $\mathbb{R}^2$ and $||\cdot||$ for the corresponding Euclidean norm. We denote by $B(0, 1)$ the unit open ball of $\mathbb{R}^2$ centered at 0, with its boundary denoted by $bd(B(0, 1))$. Finally we denote by $: \text{ the scalar product on } \mathbb{R}^{2 \times 2}$ defined by } B : C = \sum_{i=1}^{2} B_i \cdot C_i \text{ for all } B, C \in \mathbb{R}^{2 \times 2}, \text{ where } B_i \in \mathbb{R}^2 \text{ (resp. } C_i \in \mathbb{R}^2) \text{ stands for the transpose of the } i\text{-th line of } B \text{ (resp. } C) \text{ for all } i \in \{1, 2\}$.

Second, in the Tresca friction problem ($TP_\Omega$) for some $\Omega \in \mathcal{U}_{\text{ref}}$, recall that $A \in L^\infty(\Omega, \mathbb{R}^{2^4})$ stands for the stiffness tensor, assumed to be linear with constant coefficients (denoted by $a_{ijkl}$ for all $(i, j, k, l) \in \{1, 2\}^4$), and $e$ is the infinitesimal strain tensor defined by $e: v \in H^1(\Omega, \mathbb{R}^2) \to (\nabla v + \nabla v^\top)/2 \in L^2(\Omega, \mathbb{R}^{2 \times 2})$. In this paper we assume that there exists a constant $\alpha > 0$ such that all coefficients of A and $e$ (denoted by $\epsilon_{ij}$ for all $(i, j) \in \{1, 2\}^2$) satisfy

$$a_{ijkl} = a_{jikl} = a_{lkij} \quad \text{and} \quad \sum_{i=1}^{2} \sum_{j=1}^{2} \sum_{k=1}^{2} \sum_{l=1}^{2} a_{ijkl} \epsilon_{ij}(v_1)(x) \epsilon_{kl}(v_2)(x) \geq \alpha \sum_{i=1}^{2} \sum_{j=1}^{2} \epsilon_{ij}(v_1)(x) \epsilon_{ij}(v_2)(x),$$

for all $v_1, v_2 \in H^1(\Omega, \mathbb{R}^2)$ and for a.e. $x \in \Omega$. From the symmetry assumption on A, note that $Ae(v) = A\nabla v$ for all $v \in H^1(\Omega, \mathbb{R}^2)$. Moreover, since $\Gamma_D$ has a positive measure, it follows that

$$\langle \cdot, \cdot \rangle_{H_D^1(\Omega, \mathbb{R}^2)} : (H_D^1(\Omega, \mathbb{R}^2))^2 \to \mathbb{R} \quad (2.1) \\ (v_1, v_2) \mapsto \int_{\Omega} Ae(v_1) : e(v_2),$$

is a scalar product on $H_D^1(\Omega, \mathbb{R}^2)$ (see, e.g., [16, Chapter 3]) and we denote by $||\cdot||_{H_D^1(\Omega, \mathbb{R}^2)}$ the corresponding norm.

Finally, for any $\Omega \in \mathcal{U}_{\text{ref}}$, the notation $n \in C^0(\Gamma, \mathbb{R}^2)$ stands for the outward-pointing unit normal vector to $\Gamma$ and, since we deal with the two-dimensional case $d = 2$, we can fix $\tau \in C^0(\Gamma, \mathbb{R}^2)$ an oriented (with an orientation arbitrarily fixed) orthonormal vector to $n$. For any $v \in L^2(\Gamma, \mathbb{R}^2)$, one has the decomposition $v = v_n n + v_\tau$, where $v_n := v \cdot n \in L^2(\Gamma, \mathbb{R})$ and $v_\tau := v - v_n n = (v \cdot \tau)\tau \in L^2(\Gamma, \mathbb{R}^2)$. In particular, if the stress vector $Ae(v)n$ belongs to $L^2(\Gamma, \mathbb{R}^2)$ for some $v \in H^1(\Omega, \mathbb{R}^2)$, then $Ae(v)n = \sigma_n(v)n + \sigma_\tau(v)$, where $\sigma_n(v) := Ae(v)n \cdot n \in L^2(\Gamma, \mathbb{R})$ is the normal stress and $\sigma_\tau(v) := Ae(v)n - \sigma_n(v)n = (Ae(v)n \cdot \tau)\tau \in L^2(\Gamma, \mathbb{R}^2)$ is the shear stress.

## 2.2 Twice epi-differentiability and the derivative $\bar{u}'_0$

As in Introduction, let $\Omega_0 \in \mathcal{U}_{\text{ref}}$ and $\theta \in C_D^{2,\infty}(\mathbb{R}^2, \mathbb{R}^2)$ be fixed for the whole section. For any $t \ge 0$ sufficiently small, we denote by $\Omega_t := (\text{id} + t\theta)(\Omega_0) \in \mathcal{U}_{\text{ref}}$ and by $u_t := u_{\Omega_t} \in H_D^1(\Omega_t, \mathbb{R}^2)$. Then we introduce $\bar{u}_t := u_t \circ (\text{id} + t\theta) \in H_D^1(\Omega_0, \mathbb{R}^2)$ and $\bar{u}_t := (\text{I} + t\nabla\theta^\top)\bar{u}_t \in H_D^1(\Omega_0, \mathbb{R}^2)$. In the sequel, to ease the notations, we denote by $n := n_0$ and $\tau := \tau_0$.

In this section our objective is to get a characterization of the derivative $\bar{u}'_0 \in H_D^1(\Omega_0, \mathbb{R}^2)$ (if it exists). For this purpose, our methodology relies on the equality $\bar{u}_t = \text{prox}_{\phi_t}(\bar{F}_t)$ (established in Introduction) and on the application of Proposition A.5. However, as mentioned in Remark 1.1,
the equality $\bar{u}_t = \text{prox}_{\bar{\phi}_t}(\bar{F}_t)$ holds true, but considered on the Hilbert space $H_D^1(\Omega_0, \mathbb{R}^2)$ endowed with the parameterized scalar product given by

$$ (v_1, v_2) \in (H_D^1(\Omega_0, \mathbb{R}^2))^2 \mapsto \\ \int_{\Omega_0} J_t A \left[ \nabla \left( (I + t \nabla \theta^\top)^{-1} v_1 \right) (I + t \nabla \theta)^{-1} \right] : \nabla \left( (I + t \nabla \theta^\top)^{-1} v_2 \right) (I + t \nabla \theta)^{-1} \in \mathbb{R}. $$

Therefore one cannot apply Proposition A.5 directly. To overcome this difficulty, as in [4], let us add $\langle \bar{u}_t, v - \bar{u}_t \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)}$ to both members of Inequality (1.6). Then, by using the equality B : CD = BD$^\top$: C which is true for all B, C, D $\in \mathbb{R}^{2 \times 2}$, we obtain that

$$ \langle \bar{u}_t, v - \bar{u}_t \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)} + \int_{\Gamma_{T_0}} \frac{g_t J_{T_t}}{\|(I + t \nabla \theta) \tau\|} |v \cdot \tau| - \int_{\Gamma_{T_0}} \frac{g_t J_{T_t}}{\|(I + t \nabla \theta) \tau\|} |\bar{u}_t \cdot \tau| \\ \geq - \int_{\Omega_0} J_t A \left[ \nabla \left( (I + t \nabla \theta^\top)^{-1} \bar{u}_t \right) (I + t \nabla \theta)^{-1} \right] (I + t \nabla \theta^\top)^{-1} : \nabla \left( (I + t \nabla \theta^\top)^{-1} (v - \bar{u}_t) \right) \\ + \int_{\Omega_0} (I + t \nabla \theta)^{-1} f_t J_t \cdot (v - \bar{u}_t) + \int_{\Omega_0} \text{Ae}(\bar{u}_t) : e(v - \bar{u}_t), \quad \forall v \in H_D^1(\Omega_0, \mathbb{R}^2), $$

and thus we get the equality $\bar{u}_t = \text{prox}_{\bar{\phi}_t}(E_t)$ where $E_t \in H_D^1(\Omega_0, \mathbb{R}^2)$ is the unique solution to the parameterized variational equality

$$ \langle E_t, v \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)} = \int_{\Omega_0} (I + t \nabla \theta)^{-1} f_t J_t \\ - \int_{\Omega_0} J_t A \left[ \nabla \left( (I + t \nabla \theta^\top)^{-1} \bar{u}_t \right) (I + t \nabla \theta)^{-1} \right] (I + t \nabla \theta^\top)^{-1} : \nabla \left( (I + t \nabla \theta^\top)^{-1} v \right) \\ + \int_{\Omega_0} \text{Ae}(\bar{u}_t) : e(v), \quad \forall v \in H_D^1(\Omega_0, \mathbb{R}^2), $$

considered on the Hilbert space $H_D^1(\Omega_0, \mathbb{R}^2)$ endowed with the nonparameterized scalar product $\langle \cdot, \cdot \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)}$.

Furthermore, to be in accordance with the notations of Proposition A.5, we introduce the parameterized convex functional $\Phi: \mathbb{R}_+ \times H_D^1(\Omega_0, \mathbb{R}^2) \to \mathbb{R}$ defined by

$$ \Phi: \mathbb{R}_+ \times H_D^1(\Omega_0, \mathbb{R}^2) \to \mathbb{R} \qquad (2.2) \\ (t, v) \mapsto \Phi(t, v) := \bar{\phi}_t(v) = \int_{\Gamma_{T_0}} \frac{g_t J_{T_t}}{\|(I + t \nabla \theta) \tau\|} \|v_\tau\|. $$

Hence, from the equality $\bar{u}_t = \text{prox}_{\bar{\phi}(t, \cdot)}(E_t)$ satisfied on the Hilbert space $H_D^1(\Omega_0, \mathbb{R}^2)$ endowed with the nonparameterized scalar product $\langle \cdot, \cdot \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)}$, we are now in a satisfactory setting in order to apply Proposition A.5 (if its assumptions are satisfied of course). The first step is to analyze the differentiability of the map $t \in \mathbb{R}_+ \mapsto E_t \in H_D^1(\Omega_0, \mathbb{R}^2)$ at $t = 0$. For this purpose, let us recall from [22] that:

(i) the map $t \in \mathbb{R}_+ \mapsto J_t \in L^\infty(\mathbb{R}^2)$ is differentiable at $t = 0$ with derivative given by $\text{div}(\theta)$;

(ii) the map $t \in \mathbb{R}_+ \mapsto (I + t \nabla \theta)^{-1} \in L^\infty(\mathbb{R}^2, \mathbb{R}^{2 \times 2})$ is differentiable at $t = 0$ with derivative given by $-\nabla \theta$;

(iii) the map $t \in \mathbb{R}_+ \mapsto (I + t \nabla \theta^\top)^{-1} \in L^\infty(\mathbb{R}^2, \mathbb{R}^{2 \times 2})$ is differentiable at $t = 0$ with derivative given by $-\nabla \theta^\top$;
(iv) the map $t \in \mathbb{R}_+ \mapsto (I + t\nabla\theta)^{-1} f_t J_t \in L^2(\mathbb{R}^2, \mathbb{R}^2)$ is differentiable at $t = 0$ with derivative given by $\text{div}(f\theta^\top) - \nabla\theta f$;

(v) the map $t \in \mathbb{R}_+ \mapsto \frac{g_t J_{T_t}}{\|(I+t\nabla\theta)\tau\|} \in L^2(\Gamma_{T_0})$ is differentiable at $t = 0$ with derivative given by $p(\theta) := \nabla g \cdot \theta + g(\text{div}_\tau(\theta) - \nabla\theta\tau \cdot \tau)$.

**Lemma 2.1.** The map $t \in \mathbb{R}_+ \mapsto E_t \in H_D^1(\Omega_0, \mathbb{R}^2)$ is differentiable at $t = 0$ and its derivative, denoted by $E'_0 \in H_D^1(\Omega_0, \mathbb{R}^2)$, is the unique solution to the variational equality given by

$$ \langle E'_0, v \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)} = \int_{\Omega_0} (\text{div}(f\theta^\top) - \nabla\theta f) \cdot v \\ + \int_{\Omega_0} \left( (\text{Ae}(u_0))\nabla\theta^\top + \text{A}(\nabla u_0 \nabla\theta) + \text{Ae}(\nabla\theta^\top u_0) - \text{div}(\theta)\text{Ae}(u_0) \right) : \nabla v \\ + \int_{\Omega_0} \text{Ae}(u_0) : e(\nabla\theta^\top v), \quad \forall v \in H_D^1(\Omega_0, \mathbb{R}^2). \quad (2.3) $$

*Proof.* Using the Riesz representation theorem, we denote by $Z \in H_D^1(\Omega_0, \mathbb{R}^2)$ the unique solution to the above variational inequality (2.3). From linearity and using differentiability results (i), (ii), (iii), (iv), one gets

$$ \left\| \frac{E_t - E_0}{t} - Z \right\|_{H_D^1(\Omega_0, \mathbb{R}^2)} \le C(\Omega_0, A, \theta) \left( \left\| \frac{(I + t \nabla \theta)^{-1} f_t J_t - f}{t} - (\text{div}(f\theta^\top) - \nabla \theta f) \right\|_{L^2(\mathbb{R}^2, \mathbb{R}^2)} \right. \\ \left. + \| \bar{u}_t - u_0 \|_{H_D^1(\Omega_0, \mathbb{R}^2)} + \frac{o(t)}{t} \| \bar{u}_t \|_{H_D^1(\Omega_0, \mathbb{R}^2)} \right), $$

for all $t > 0$ sufficiently small, where $C(\Omega_0, A, \theta) > 0$ is a constant which depends on $\Omega_0$, A and $\theta$, and where $o$ stands for the standard Bachmann-Landau notation, with $\frac{|o(t)|}{t} \to 0$ when $t \to 0^+$. Therefore, to conclude the proof, we only need to prove the continuity of the map $t \in \mathbb{R}_+ \mapsto \bar{u}_t \in H_D^1(\Omega_0, \mathbb{R}^2)$ at $t = 0$. For this purpose, take $v = u_0$ in the variational formulation of $\bar{u}_t$ and $v = \bar{u}_t$ in the variational formulation of $u_0$ to get that

$$ \| \bar{u}_t - u_0 \|_{H_D^1(\Omega_0, \mathbb{R}^2)} \le C(\Omega_0, A, \theta) \left( \| (I + t \nabla \theta^\top) f_t J_t - f \|_{L^2(\mathbb{R}^2, \mathbb{R}^2)} + \left\| \frac{g_t J_{T_t}}{\|(I + t \nabla \theta)\tau\|} - g \right\|_{L^2(\Gamma_{T_0}, \mathbb{R})} \right. \\ \left. + \| \bar{u}_t \|_{H_D^1(\Omega_0, \mathbb{R}^2)} (t + o(t)) \right), $$

for all $t \ge 0$ sufficiently small. Hence, to conclude the proof, we only need to prove that the map $t \in \mathbb{R}_+ \mapsto \|\bar{u}_t\|_{H_D^1(\Omega_0, \mathbb{R}^2)} \in \mathbb{R}$ is bounded for $t \ge 0$ sufficiently small. For this purpose, take $v = 0$ in the variational formulation of $\bar{u}_t$ to get that

$$ \| \bar{u}_t \|_{H_D^1(\Omega_0, \mathbb{R}^2)}^2 \le C(\Omega_0, A, \theta) \left( \| (I + t \nabla \theta^\top) f_t J_t \|_{L^2(\mathbb{R}^2, \mathbb{R}^2)} + \left\| \frac{g_t J_{T_t}}{\|(I + t \nabla \theta)\tau\|} \right\|_{L^2(\Gamma_{T_0}, \mathbb{R})} \right) \| \bar{u}_t \|_{H_D^1(\Omega_0, \mathbb{R}^2)} $$
$$+ C(\Omega_0, A, \theta) \|\bar{u}_t\|_{H_D^1(\Omega_0, \mathbb{R}^2)}^2 (t+o(t)),$$

for all $t \ge 0$ sufficiently small. Thus one deduces

$$\|\bar{u}_t\|_{H_D^1(\Omega_0, \mathbb{R}^2)} \le \frac{C(\Omega_0, A, \theta) \left( \|(I+t\nabla\theta^\top) f_t J_t\|_{L^2(\mathbb{R}^2, \mathbb{R}^2)} + \left\| \frac{g_t J_{T_t}}{\|(I+t\nabla\theta)^\top\tau\|} \right\|_{L^2(\Gamma_{T_0})} \right)}{1 - C(\Omega_0, A, \theta) (t+o(t))},$$

for all $t \ge 0$ sufficiently small, and using the continuity of the map $t \in \mathbb{R}_+ \mapsto (I+t\nabla\theta^\top)f_t J_t \in L^2(\mathbb{R}^2, \mathbb{R}^2)$ (see (iv)) and of the map $t \in \mathbb{R}_+ \mapsto \frac{g_t J_{T_t}}{\|(I+t\nabla\theta)^\top\tau\|} \in L^2(\Gamma_{T_0})$ (see (v)), the proof is complete. $\square$

Now the second step is to investigate the twice epi-differentiability of the parameterized convex functional $\Phi$ defined in (2.2), as we did in our previous paper [11] from which the next two lemmas are extracted. Precisely, to derive the next two lemmas, one has to apply [11, Propositions 2.18 and 2.23] on the particular case given by the expression (2.2). For the needs of these lemmas, and to avoid any confusion, we recall that the notation $\partial$ stands for the notion of subdifferential (see Appendix A). We also introduce the notation $x_{\tau(s)} := (x \cdot \tau(s))\tau(s)$ for all $x \in \mathbb{R}^2$ and all $s \in \Gamma_0$. Similarly we will use, for all $s \in \Gamma_0$, the tangential norm map given by $\|\cdot|_{\tau(s)}\| : x \in \mathbb{R}^2 \mapsto \|x_{\tau(s)}\| = |x \cdot \tau(s)| \in \mathbb{R}^+$.

**Lemma 2.2** (Second-order difference quotient functions of $\Phi$). For all $t > 0$, $u \in H_D^1(\Omega_0, \mathbb{R}^2)$ and $v \in \partial\Phi(0, \cdot)(u)$, it holds that

$$\Delta_t^2 \Phi(u | v)(\varphi) = \int_{\Gamma_{T_0}} \Delta_t^2 G(s)(u(s) | \sigma_{\tau}(v)(s))(\varphi(s)) \, ds, \quad (2.4)$$

for all $\varphi \in H_D^1(\Omega_0, \mathbb{R}^2)$, where, for almost all $s \in \Gamma_{T_0}$, $\Delta_t^2 G(s)(u(s)|\sigma_{\tau}(v)(s))$ stands for the second-order difference quotient functions of $G(s)$ at $u(s) \in \mathbb{R}^2$ for $\sigma_{\tau}(v)(s) \in \partial G(s)(0, \cdot)(u(s)) = g(s)\partial\|\cdot|_{\tau(s)}\|(u(s))$, with $G(s)$ defined by

$$G(s) : \mathbb{R}_+ \times \mathbb{R}^2 \to \mathbb{R} \\ (t, x) \mapsto G(s)(t, x) := \frac{g_t(s) J_{T_t}(s)}{\|(I + t \nabla \theta(s)) \tau(s)\|} \|x_{\tau(s)}\|.$$

**Remark 2.3.** The assumption that $\Gamma_0$ is of class $C^1$ is made to ensure that $n \in C^0(\Gamma, \mathbb{R})$ which gives us, for all $s \in \Gamma_{T_0}$ and all $x \in \mathbb{R}^2$, the continuity of the map $s \in \Gamma_{T_0} \mapsto \|x_{\tau(s)}\| \in \mathbb{R}_+$. This property is used in the proof of [11, Proposition 2.18] (precisely, in the proof of [11, Lemma 2.16]).

**Lemma 2.4** (Second-order epi-derivative of $G(s)$). Assume that, for almost all $s \in \Gamma_{T_0}$, $g$ has a directional derivative at $s$ in any direction. Then, for almost all $s \in \Gamma_{T_0}$, the map $G(s)$ is twice epi-differentiable at any $x \in \mathbb{R}^2$ for all $y \in \partial G(s)(0, \cdot)(u(s)) = g(s)\partial\|\cdot|_{\tau(s)}\|(x)$ with

$$D_e^2 G(s)(x | y)(z) := \begin{cases} p(\theta)(s) \frac{x_{\tau(s)}}{\|x_{\tau(s)}\|} \cdot z & \text{if } x_{\tau(s)} \neq 0, \\ t \mathrm{N}_{\overline{B(0,1) \cap (\mathbb{R}n(s))^{\perp}}} \left( \frac{-y}{g(s)} \right)(z) + p(\theta)(s) \frac{y}{g(s)} \cdot z & \text{if } x_{\tau(s)} = 0, \end{cases}$$

for all $z \in \mathbb{R}^2$, where $p(\theta) \in L^2(\Gamma_{T_0})$ is defined as the derivative at $t = 0$ of the map $t \in \mathbb{R}_+ \mapsto \frac{g_t J_{T_t}}{\|(I+t\nabla\theta)^\top\tau\|} \in L^2(\Gamma_{T_0})$ (see (v)), $N_{\overline{B(0,1) \cap (\mathbb{R}n(s))^{\perp}}} (\frac{y}{g(s)})$ is the normal cone to $\overline{B(0,1)} \cap (\mathbb{R}n(s))^{\perp}$ at $\frac{y}{g(s)}$ given by

$$N_{\overline{B(0,1) \cap (\mathbb{R}n(s))^{\perp}}} \left( \frac{y}{g(s)} \right) = \begin{cases} \mathbb{R}n(s) & \text{if } \frac{y}{g(s)} \in \mathbb{B}(0,1) \cap (\mathbb{R}n(s))^{\perp}, \\ \mathbb{R}n(s) + \mathbb{R}_{+} \frac{y}{g(s)} & \text{if } \frac{y}{g(s)} \in \text{bd}(\mathbb{B}(0,1)) \cap (\mathbb{R}n(s))^{\perp}, \end{cases}$$
and $\iota_{N_{\overline{B(0,1)} \cap (\mathbb{R}n(s))^{\perp}}(\frac{y}{g(s)})}^{}$ stands for the indicator function of $N_{\overline{B(0,1)} \cap (\mathbb{R}n(s))^{\perp}}(\frac{y}{g(s)})$ which is defined by $\iota_{N_{\overline{B(0,1)} \cap (\mathbb{R}n(s))^{\perp}}(\frac{y}{g(s)})}(z) := 0$ if $z \in N_{\overline{B(0,1)} \cap (\mathbb{R}n(s))^{\perp}}(\frac{y}{g(s)})$, and $\iota_{N_{\overline{B(0,1)} \cap (\mathbb{R}n(s))^{\perp}}(\frac{y}{g(s)})}^{}(z) := +\infty$ otherwise.

We are now in a position to prove, under appropriate assumptions, that the map $t \in \mathbb{R}_+ \mapsto \bar{u}_t \in H_D^1(\Omega_0, \mathbb{R}^2)$ is differentiable at $t = 0$.

**Proposition 2.5** (The derivative $\bar{u}'_0$). Assume that:

(H1) for almost all $s \in \Gamma_{T_0}$, $g$ has a directional derivative at $s$ in any direction.

(H2) the parameterized convex functional $\Phi$ defined in (2.2) is twice epi-differentiable (see Definition A.4) at $u_0$ for $E_0 - u_0 \in \partial\Phi(0, \cdot)(u_0)$, with

$$D_e^2 \Phi(u_0 | E_0 - u_0)(\varphi) = \int_{\Gamma_{T_0}} D_e^2 G(s)(u_0(s) | \sigma_{\tau}(E_0 - u_0)(s))(\varphi(s)) \, ds, \quad (2.5)$$

for all $\varphi \in H_D^1(\Omega_0, \mathbb{R}^2)$.

Then the map $t \in \mathbb{R}_+ \mapsto \bar{u}_t \in H_D^1(\Omega_0, \mathbb{R}^2)$ is differentiable at $t = 0$ and its derivative $\bar{u}'_0 \in K_0$ is the unique solution to the variational inequality

$$\begin{aligned} \langle \bar{u}'_0, \varphi - \bar{u}'_0 \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)} &\geq - \int_{\Omega_0} \operatorname{div} (\operatorname{div}(\mathrm{Ae}(u_0)) \theta^{\top}) \cdot (\varphi - \bar{u}'_0) \\ &\quad + \int_{\Omega_0} ((\mathrm{Ae}(u_0)) \nabla \theta^{\top} + \mathrm{A}(\nabla u_0 \nabla \theta) + \mathrm{Ae}(\nabla \theta^{\top} u_0) - \operatorname{div}(\theta) \mathrm{Ae}(u_0)) : \nabla (\varphi - \bar{u}'_0) \\ &\quad + \langle \mathrm{Ae}(u_0) \mathrm{n}, \nabla \theta^{\top} (\varphi - \bar{u}'_0) \rangle_{H^{-1/2}(\Gamma_0, \mathbb{R}^2) \times H^{1/2}(\Gamma_0, \mathbb{R}^2)} - \int_{\Gamma_{T_{0_N}^{u_0,g}}} p(\theta) \frac{u_{0_{\tau}}}{\|u_{0_{\tau}}\|} \cdot (\varphi_{\tau} - \bar{u}'_{0_{\tau}}) \\ &\quad + \int_{\Gamma_{T_{0_S}^{u_0,g}}} p(\theta) \frac{\sigma_{\tau}(u_0)}{g} \cdot (\varphi_{\tau} - \bar{u}'_{0_{\tau}}), \quad \forall \varphi \in K_0, \end{aligned}$$

where $K_0$ is the nonempty closed convex subset of $H_D^1(\Omega_0, \mathbb{R}^2)$ given by

$$K_0 = \left\{ \varphi \in H_D^1(\Omega_0, \mathbb{R}^2) \mid \varphi_{\tau} = 0 \text{ a.e. on } \Gamma_{T_{0_N}^{u_0,g}} \text{ and } \varphi_{\tau} \in \mathbb{R}_{-} \frac{\sigma_{\tau}(u_0)}{g} \text{ a.e. on } \Gamma_{T_{0_S}^{u_0,g}} \right\}, \quad (2.6)$$

and where $\Gamma_{T_0}$ is decomposed, up to a null set, as $\Gamma_{T_0_N^{u_0,g}} \cup \Gamma_{T_0_D^{u_0,g}} \cup \Gamma_{T_0_S^{u_0,g}}$ with

$$\begin{aligned} \Gamma_{T_0_N^{u_0,g}} &:= \{s \in \Gamma_{T_0} \mid u_{0_{\tau}}(s) \neq 0\}, \\ \Gamma_{T_0_D^{u_0,g}} &:= \{s \in \Gamma_{T_0} \mid u_{0_{\tau}}(s) = 0 \text{ and } \frac{\sigma_{\tau}(u_0)(s)}{g(s)} \in \mathrm{B}(0,1) \cap (\mathbb{R}n(s))^{\perp}\}, \\ \Gamma_{T_0_S^{u_0,g}} &:= \{s \in \Gamma_{T_0} \mid u_{0_{\tau}}(s) = 0 \text{ and } \frac{\sigma_{\tau}(u_0)(s)}{g(s)} \in \mathrm{bd}(\mathrm{B}(0,1)) \cap (\mathbb{R}n(s))^{\perp}\}. \end{aligned}$$

*Proof.* From Hypotheses (H1), (H2) and Lemma 2.4, it follows that

$$\begin{aligned} D_e^2 \Phi(u_0 | E_0 - u_0)(\varphi) = & \int_{\Gamma_{T_0_N^{u_0,g}}} p(\theta) \frac{u_{0_{\tau}}}{\|u_{0_{\tau}}\|} \cdot \varphi_{\tau} + \int_{\Gamma_{T_0} \setminus \Gamma_{T_0_N^{u_0,g}}} p(\theta) \frac{\sigma_{\tau}(E_0 - u_0)}{g} \cdot \varphi_{\tau} \\ & + \int_{\Gamma_{T_0} \setminus \Gamma_{T_0_N^{u_0,g}}} \iota_{N_{\overline{B(0,1)} \cap (\mathbb{R}n(s))^{\perp}}(\frac{\sigma_{\tau}(E_0 - u_0)(s)}{g(s)})}^{}(\varphi(s)) \, ds, \end{aligned}$$
which can be rewritten as

$$D_e^2 \Phi(u_0 | E_0 - u_0)(\varphi) = \int_{\Gamma_{T_0 N}^{u_0, g}} p(\theta) \frac{u_{0_\tau}}{\|u_{0_\tau}\|} \cdot \varphi_\tau + \int_{\Gamma_{T_0} \setminus \Gamma_{T_0 N}^{u_0, g}} p(\theta) \frac{\sigma_\tau (E_0 - u_0)}{g} \cdot \varphi_\tau + \iota_{K_0}(\varphi),$$

for all $\varphi \in H_D^1(\Omega_0, \mathbb{R}^2)$, where $K_0$ is the nonempty closed convex subset of $H_D^1(\Omega_0, \mathbb{R}^2)$ defined by

$$K_0 := \left\{ \varphi \in H_D^1(\Omega_0, \mathbb{R}^2) \mid \varphi(s) \in N_{\overline{B(0,1) \cap (\mathbb{R}n(s))}^\perp} \left( \frac{\sigma_\tau (E_0 - u_0)(s)}{g(s)} \right) \quad \text{for almost all } s \in \Gamma_{T_0} \setminus \Gamma_{T_0 N}^{u_0, g} \right\}.$$

Moreover, since $E_0 = F_{\Omega_0} \in H_D^1(\Omega_0, \mathbb{R}^2)$ is solution to the Dirichlet-Neumann problem (1.2) with $\Omega = \Omega_0$, then $\sigma_\tau(E_0) = 0$ a.e. on $\Gamma_{T_0}$. Thus it follows that $K_0$ is given by (2.6). Now, since $D_e^2 \Phi(u_0|E_0-u_0)$ is a proper function on $H_D^1(\Omega_0, \mathbb{R}^2)$ and the map $t \in \mathbb{R}_+ \mapsto E_t \in H_D^1(\Omega_0, \mathbb{R}^2)$ is differentiable at $t=0$ with its derivative $E'_0 \in H_D^1(\Omega_0, \mathbb{R}^2)$ being solution to the variational inequality (2.3), we apply Proposition A.5 to deduce that the map $t \in \mathbb{R}_+ \mapsto \bar{u}_t \in H_D^1(\Omega_0, \mathbb{R}^2)$ is differentiable at $t=0$, and its derivative $\bar{u}'_0 \in H_D^1(\Omega_0, \mathbb{R}^2)$ satisfies

$$\bar{u}'_0 = \text{prox}_{D_e^2 \Phi(u_0|E_0-u_0)}(E'_0),$$

which, from the definition of the proximal operator (see Definition A.1), leads to

$$\langle E'_0 - \bar{u}'_0, \varphi - \bar{u}'_0 \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)} \le D_e^2 \Phi(u_0 | E_0 - u_0)(\varphi) - D_e^2 \Phi(u_0 | E_0 - u_0)(\bar{u}'_0),$$

for all $\varphi \in H_D^1(\Omega, \mathbb{R}^2)$. Hence we get that

$$\langle E'_0 - \bar{u}'_0, \varphi - \bar{u}'_0 \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)} \le +\iota_{K_0}(\varphi) - \iota_{K_0}(\bar{u}'_0) \\ + \int_{\Gamma_{T_0 N}^{u_0, g}} p(\theta) \frac{u_{0_\tau}}{\|u_{0_\tau}\|} \cdot (\varphi_\tau - \bar{u}'_{0_\tau}) + \int_{\Gamma_{T_0} \setminus \Gamma_{T_0 N}^{u_0, g}} p(\theta) \frac{\sigma_\tau (E_0 - u_0)}{g_0} \cdot (\varphi_\tau - \bar{u}'_{0_\tau}),$$

for all $\varphi \in H_D^1(\Omega, \mathbb{R}^2)$. Since $\varphi_\tau = 0$ a.e. on $\Gamma_{T_0 N}^{u_0, g}$ for all $\varphi \in K_0$, one deduces that $\bar{u}'_0 \in K_0$ satisfies

$$\langle \bar{u}'_0, \varphi - \bar{u}'_0 \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)} \ge \int_{\Omega_0} (\text{div}(f\theta^\top) - \nabla\theta f) \cdot (\varphi - \bar{u}'_0) \\ + \int_{\Omega_0} ((\text{Ae}(u_0))\nabla\theta^\top + \text{A}(\nabla u_0 \nabla\theta) + \text{Ae}(\nabla\theta^\top u_0) - \text{div}(\theta)\text{Ae}(u_0)) : \nabla(\varphi - \bar{u}'_0) \\ + \int_{\Omega_0} \text{Ae}(u_0) : e(\nabla\theta^\top(\varphi - \bar{u}'_0)) - \int_{\Gamma_{T_0 N}^{u_0, g}} p(\theta) \frac{u_{0_\tau}}{\|u_{0_\tau}\|} \cdot (\varphi_\tau - \bar{u}'_{0_\tau}) \\ + \int_{\Gamma_{T_0 S}^{u_0, g}} p(\theta) \frac{\sigma_\tau(u_0)}{g} \cdot (\varphi_\tau - \bar{u}'_{0_\tau}),$$

for all $\varphi \in K_0$. Using the equality $-\text{div}(\text{Ae}(u_0)) = f$ in $H^1(\Omega_0, \mathbb{R}^2)$ and the divergence formula (see Proposition B.1), the proof is complete. $\square$

**Remark 2.6.** Note that Hypothesis (H2) corresponds to the inversion of the symbols ME-lim and $\int_{\Gamma_{T_0}}$ in Equality (2.4), which is an open question in general. We refer to [10, Appendix A] and [11, Remark 2.26] for additional comments and some sufficient conditions for this inversion.
**Remark 2.7.** In Proposition 2.5 (and in its proof), note that the set $\mathcal{K}_0$ corresponds to the set $\mathcal{K}_{u_0, \frac{\sigma_{\tau}(E_0-u_0)}{g}}$ with the notations introduced in our previous paper [11] (see [11, proof of Theorem 2.25]).

## 2.3 Directional material derivative $\bar{u}'_0$ and directional shape derivative $u'_0$

From Proposition 2.5 and since $\bar{u}_t = (I + t\nabla\theta^\top)^{-1}\bar{u}_t$ for all $t \ge 0$, it is possible now to state and prove the first main result of this paper that characterizes the directional material derivative $\bar{u}'_0$.

**Theorem 2.8 (Directional material derivative $\bar{u}'_0$).** Consider the framework of Proposition 2.5. Then the map $t \in \mathbb{R}_+ \mapsto \bar{u}_t \in H_D^1(\Omega_0, \mathbb{R}^2)$ is differentiable at $t = 0$ and its derivative $\bar{u}'_0 \in \mathcal{K}_0 - \nabla\theta^\top u_0$ (that is, the directional material derivative) is the unique solution to the variational inequality

$$ \langle \bar{u}'_0, v - \bar{u}'_0 \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)} \geq - \int_{\Omega_0} \mathrm{div} \left( \mathrm{div} \left( \mathrm{Ae}(u_0) \right) \theta^\top \right) \cdot (v - \bar{u}'_0) \\ + \int_{\Omega_0} \left( \left( \mathrm{Ae}(u_0) \right) \nabla\theta^\top + \mathrm{A} \left( \nabla u_0 \nabla\theta \right) - \mathrm{div}(\theta) \mathrm{Ae}(u_0) \right) : \nabla(v - \bar{u}'_0) \\ + \langle \mathrm{Ae}(u_0) \mathrm{n}, \nabla\theta^\top (v - \bar{u}'_0) \rangle_{H^{-1/2}(\Gamma_0, \mathbb{R}^2) \times H^{1/2}(\Gamma_0, \mathbb{R}^2)} - \int_{\Gamma_{T_{0N}^{u_0,g}}} p(\theta) \frac{u_{0\tau}}{\|u_{0\tau}\|} \cdot (v_\tau - \bar{u}'_{0\tau}) \\ + \int_{\Gamma_{T_{0S}^{u_0,g}}} p(\theta) \frac{\sigma_\tau(u_0)}{g} \cdot (v_\tau - \bar{u}'_{0\tau}), \quad \forall v \in \mathcal{K}_0 - \nabla\theta^\top u_0, \quad (2.7) $$

where

$$ \mathcal{K}_0 - \nabla\theta^\top u_0 = \left\{ v \in H_D^1(\Omega_0, \mathbb{R}^2) \mid \begin{array}{l} v_\tau = -(\nabla\theta^\top u_0)_\tau \quad \text{a.e. on } \Gamma_{T_{0D}^{u_0,g}} \\ \text{and } (v_\tau + (\nabla\theta^\top u_0)_\tau) \in \mathbb{R}_- \frac{\sigma_\tau(u_0)}{g} \quad \text{a.e. on } \Gamma_{T_{0S}^{u_0,g}} \end{array} \right\}. $$

*Proof.* Since $\bar{u}_t = (I + t\nabla\theta^\top)^{-1}\bar{u}_t$ for all $t \ge 0$, one deduces from Proposition 2.5 that the map $t \in \mathbb{R}_+ \mapsto \bar{u}_t \in H_D^1(\Omega_0, \mathbb{R}^2)$ is differentiable at $t = 0$ with $\bar{u}'_0 = \bar{u}'_0 - \nabla\theta^\top u_0 \in H_D^1(\Omega_0, \mathbb{R}^2)$. Moreover, from the variational inequality satisfied by $\bar{u}'_0$, one deduces that

$$ \langle \bar{u}'_0 + \nabla\theta^\top u_0, \varphi - \nabla\theta^\top u_0 - \bar{u}'_0 \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)} \geq - \int_{\Omega_0} \mathrm{div} \left( \mathrm{div} \left( \mathrm{Ae}(u_0) \right) \theta^\top \right) \cdot \left( \varphi - \nabla\theta^\top u_0 - \bar{u}'_0 \right) \\ + \int_{\Omega_0} \left( \left( \mathrm{Ae}(u_0) \right) \nabla\theta^\top + \mathrm{A} \left( \nabla u_0 \nabla\theta \right) + \mathrm{Ae} \left( \nabla\theta^\top u_0 \right) - \mathrm{div}(\theta) \mathrm{Ae}(u_0) \right) : \nabla \left( \varphi - \nabla\theta^\top u_0 - \bar{u}'_0 \right) \\ + \langle \mathrm{Ae}(u_0) \mathrm{n}, \nabla\theta^\top \left( \varphi - \nabla\theta^\top u_0 - \bar{u}'_0 \right) \rangle_{H^{-1/2}(\Gamma_0, \mathbb{R}^2) \times H^{1/2}(\Gamma_0, \mathbb{R}^2)} \\ - \int_{\Gamma_{T_{0N}^{u_0,g}}} p(\theta) \frac{u_{0\tau}}{\|u_{0\tau}\|} \cdot \left( \varphi_\tau - (\nabla\theta^\top u_0)_\tau - \bar{u}'_{0\tau} \right) + \int_{\Gamma_{T_{0S}^{u_0,g}}} p(\theta) \frac{\sigma_\tau(u_0)}{g} \cdot \left( \varphi_\tau - (\nabla\theta^\top u_0)_\tau - \bar{u}'_{0\tau} \right), $$

for all $\varphi \in \mathcal{K}_0$, and this is also

$$ \langle \bar{u}'_0 + \nabla\theta^\top u_0, v - \bar{u}'_0 \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)} \geq - \int_{\Omega_0} \mathrm{div} \left( \mathrm{div} \left( \mathrm{Ae}(u_0) \right) \theta^\top \right) \cdot (v - \bar{u}'_0) $$
$$ + \int_{\Omega_0} ((\mathrm{Ae}(u_0)) \nabla \theta^\top + \mathrm{A}(\nabla u_0 \nabla \theta) + \mathrm{Ae}(\nabla \theta^\top u_0) - \mathrm{div}(\theta) \mathrm{Ae}(u_0)) : \nabla (v - \bar{u}'_0) \\ + \langle \mathrm{Ae}(u_0) \mathrm{n}, \nabla \theta^\top (v - \bar{u}'_0) \rangle_{\mathrm{H}^{-1/2}(\Gamma_0, \mathbb{R}^2) \times \mathrm{H}^{1/2}(\Gamma_0, \mathbb{R}^2)} \\ - \int_{\Gamma_{T_{0N}^{u_0,g}}} p(\theta) \frac{u_{0\tau}}{\|u_{0\tau}\|} \cdot (v_\tau - \bar{u}'_{0\tau}) + \int_{\Gamma_{T_{0S}^{u_0,g}}} p(\theta) \frac{\sigma_\tau(u_0)}{g} \cdot (v_\tau - \bar{u}'_{0\tau}), $$

for all $v \in \mathcal{K}_0 - \nabla \theta^\top u_0$, which concludes the proof. $\square$

<figure><img src="image_2.png" /></figure>

The presentation of Theorem 2.8 can be improved under additional regularity assumptions.

**Corollary 2.9.** Consider the framework of Proposition 2.5 with the additional assumption that $u_0 \in H^3(\Omega_0, \mathbb{R}^2)$. Then the directional material derivative $\bar{u}'_0 \in \mathcal{K}_0 - \nabla \theta^\top u_0$ is the unique weak solution to the tangential Signorini problem given by

$$ \left\{ \begin{array}{l} -\mathrm{div}(\mathrm{Ae}(\bar{u}'_0)) + \mathrm{div}(\mathrm{Ae}(\nabla u_0 \theta)) = 0 \quad \text{in } \Omega_0, \\ \bar{u}'_0 = 0 \quad \text{on } \Gamma_D, \\ \sigma_n(\bar{u}'_0) - \xi^m(\theta)_n = 0 \quad \text{on } \Gamma_{T_0}, \\ \sigma_\tau(\bar{u}'_0) + p(\theta) \frac{u_{0\tau}}{\|u_{0\tau}\|} - \xi^m(\theta)_\tau = 0 \quad \text{on } \Gamma_{T_{0N}^{u_0,g}}, \\ \bar{u}'_{0\tau} + (\nabla \theta^\top u_0)_\tau = 0 \quad \text{on } \Gamma_{T_{0D}^{u_0,g}}, \\ (\bar{u}'_{0\tau} + (\nabla \theta^\top u_0)_\tau) \in \mathbb{R} - \frac{\sigma_\tau(u_0)}{g} \\ \text{and } (\sigma_\tau(\bar{u}'_0) - p(\theta) \frac{\sigma_\tau(u_0)}{g} - \xi^m(\theta)_\tau) \cdot \frac{\sigma_\tau(u_0)}{g} \le 0 \\ \text{and } (\bar{u}'_{0\tau} + (\nabla \theta^\top u_0)_\tau) \cdot (\sigma_\tau(\bar{u}'_0) - p(\theta) \frac{\sigma_\tau(u_0)}{g} - \xi^m(\theta)_\tau) = 0 \quad \text{on } \Gamma_{T_{0S}^{u_0,g}}. \end{array} \right. $$

where $\xi^m(\theta) := ((\mathrm{Ae}(u_0))\nabla\theta^\top + \mathrm{A}(\nabla u_0 \nabla\theta) + (\nabla\theta - \mathrm{div}(\theta)\mathrm{I})\mathrm{Ae}(u_0)) \in L^2(\Gamma_{T_0}, \mathbb{R}^2)$.

*Proof.* Since $u_0 \in H^2(\Omega_0, \mathbb{R}^2)$ and $\theta \in C_D^{2,\infty}(\mathbb{R}^2, \mathbb{R}^2)$, it holds that

$$ \mathrm{div}((\mathrm{Ae}(u_0))\nabla\theta^\top + \mathrm{A}(\nabla u_0 \nabla\theta) - \mathrm{div}(\theta)\mathrm{Ae}(u_0)) \in L^2(\Omega_0, \mathbb{R}^2). $$

Thus, using the divergence formula (see Proposition B.1) in Inequality (2.7), we get that

$$ \langle \bar{u}'_0, v - \bar{u}'_0 \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)} \ge \int_{\Gamma_{T_0}} \xi^m(\theta) \cdot (v - \bar{u}'_0) \\ - \int_{\Omega_0} \mathrm{div}((\mathrm{div}(\mathrm{Ae}(u_0))\theta^\top + (\mathrm{Ae}(u_0))\nabla\theta^\top + \mathrm{A}(\nabla u_0 \nabla\theta) - \mathrm{div}(\theta)\mathrm{Ae}(u_0)) \cdot (v - \bar{u}'_0) \\ - \int_{\Gamma_{T_{0N}^{u_0,g}}} p(\theta) \frac{u_{0\tau}}{\|u_{0\tau}\|} \cdot (v_\tau - \bar{u}'_{0\tau}) + \int_{\Gamma_{T_{0S}^{u_0,g}}} p(\theta) \frac{\sigma_\tau(u_0)}{g} \cdot (v_\tau - \bar{u}'_{0\tau}), \quad (2.8) $$

for all $v \in \mathcal{K}_0 - \nabla \theta^\top u_0$. Furthermore, one has $\mathrm{div}(\mathrm{Ae}(\nabla u_0 \theta)) \in L^2(\Omega_0, \mathbb{R}^2)$ from the fact that $u_0 \in H^3(\Omega_0, \mathbb{R}^2)$. Thus, using the equality

$$ \mathrm{div}(\mathrm{Ae}(\nabla u_0 \theta)) = \mathrm{div}((\mathrm{div}(\mathrm{Ae}(u_0))\theta^\top + (\mathrm{Ae}(u_0))\nabla\theta^\top + \mathrm{A}(\nabla u_0 \nabla\theta) - \mathrm{div}(\theta)\mathrm{Ae}(u_0)), \text{ in } L^2(\Omega_0, \mathbb{R}^2), \text{ it follows that} $$
$$ \langle \bar{u}'_0, v - \bar{u}'_0 \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)} \geq - \int_{\Omega_0} \mathrm{div}(\mathrm{Ae}(\nabla u_0 \theta)) \cdot (v - \bar{u}'_0) + \int_{\Gamma_{T_0}} \xi^m(\theta) \cdot (v - \bar{u}'_0) \\ - \int_{\Gamma_{T_{0,N}^{u_0,g}}} p(\theta) \frac{u_{0,\tau}}{\|u_{0,\tau}\|} \cdot (v_\tau - \bar{u}'_{0,\tau}) + \int_{\Gamma_{T_{0,S}^{u_0,g}}} p(\theta) \frac{\sigma_\tau(u_0)}{g} \cdot (v_\tau - \bar{u}'_{0,\tau}), $$

for all $v \in \mathcal{K}_0 - \nabla\theta^\top u_0$, which corresponds to the weak formulation of the expected tangential Signorini problem (see [11, Section 2.1.2] for details). □

**Remark 2.10.** Note that, from the proof of Corollary 2.9, one can get, under the weaker assumption $u_0 \in H^2(\Omega_0, \mathbb{R}^2)$, that the directional material derivative $\bar{u}'_0$ is the solution to the variational inequality (2.8) which is, from [11, Section 2.1.2], the weak formulation of a tangential Signorini problem, with the source term given by $-\mathrm{div}(\mathrm{div}(\mathrm{Ae}(u_0))\theta^\top + (\mathrm{Ae}(u_0))\nabla\theta^\top + \mathrm{A}(\nabla u_0 \nabla\theta) - \mathrm{div}(\theta)\mathrm{Ae}(u_0)) \in L^2(\Omega_0, \mathbb{R}^2)$.

Thanks to Corollary 2.9, we are now in a position to characterize the directional shape derivative $u'_0$.

**Corollary 2.11 (Directional shape derivative $u'_0$).** Consider the framework of Proposition 2.5 with the additional assumptions that $u_0 \in H^3(\Omega_0, \mathbb{R}^2)$ and that $\Gamma_0$ is of class $C^3$. Then the directional shape derivative, defined by $u'_0 := \bar{u}'_0 - \nabla u_0 \theta \in \mathcal{K}_0 - \nabla\theta^\top u_0 - \nabla u_0 \theta$, is the unique weak solution to the tangential Signorini problem

$$ \left\{ \begin{array}{ll} -\mathrm{div}(\mathrm{Ae}(u'_0)) = 0 & \mathrm{in} \ \Omega_0, \\ u'_0 = 0 & \mathrm{on} \ \Gamma_D, \\ \sigma_n(u'_0) - \xi^s(\theta)_n = 0 & \mathrm{on} \ \Gamma_{T_0}, \\ \sigma_\tau(u'_0) + p(\theta) \frac{u_{0,\tau}}{\|u_{0,\tau}\|} - \xi^s(\theta)_\tau = 0 & \mathrm{on} \ \Gamma_{T_{0,N}^{u_0,g}}, \\ u'_{0,\tau} - W(\theta)_\tau = 0 & \mathrm{on} \ \Gamma_{T_{0,D}^{u_0,g}}, \\ (u'_{0,\tau} - W(\theta)_\tau) \in \mathbb{R}_- \frac{\sigma_\tau(u_0)}{g} & \\ \text{and } (\sigma_\tau(\bar{u}'_0) - p(\theta) \frac{\sigma_\tau(u_0)}{g} - \xi^s(\theta)_\tau) \cdot \frac{\sigma_\tau(u_0)}{g} \le 0 & \\ \text{and } (u'_{0,\tau} - W(\theta)_\tau) \cdot (\sigma_\tau(\bar{u}'_0) - p(\theta) \frac{\sigma_\tau(u_0)}{g} - \xi^s(\theta)_\tau) = 0 & \mathrm{on} \ \Gamma_{T_{0,S}^{u_0,g}}. \end{array} \right. $$

$$ \text{where } W(\theta) := -\nabla\theta^\top u_0 - \nabla u_0 \theta \in H^1(\Omega_0, \mathbb{R}^2), $$

$$ \xi^s(\theta) := \theta \cdot n (\partial_n (\mathrm{Ae}(u_0)n) - \partial_n (\mathrm{Ae}(u_0)n)) + \mathrm{Ae}(u_0) \nabla_\tau (\theta \cdot n) - \nabla (\mathrm{Ae}(u_0)n)\theta \\ + (\nabla\theta - \mathrm{div}_\tau(\theta)I) \mathrm{Ae}(u_0)n \in L^2(\Gamma_{T_0}, \mathbb{R}^2), $$

where $\partial_n (\mathrm{Ae}(u_0)n) := \nabla (\mathrm{Ae}(u_0)n)n$ stands for the normal derivative of $\mathrm{Ae}(u_0)n$, and $\partial_n (\mathrm{Ae}(u_0))$ is the matrix whose the i-th line is the transpose of the vector $\partial_n (\mathrm{Ae}(u_0)_i) := \nabla (\mathrm{Ae}(u_0)_i)n$, where $\mathrm{Ae}(u_0)_i$ is the transpose of the i-th line of the matrix $\mathrm{Ae}(u_0)$, for all $i \in \{1, 2\}$.

**Proof.** Since $u'_0 := \bar{u}'_0 - \nabla u_0 \theta$, one deduces from the weak formulation of $\bar{u}'_0$ and the divergence formula (see Proposition B.1) that

$$ \langle u'_0, v - u'_0 \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)} \geq \\ \int_{\Omega_0} (\mathrm{div}(\mathrm{Ae}(u_0))\theta^\top + (\mathrm{Ae}(u_0))\nabla\theta^\top + \mathrm{A}(\nabla u_0 \nabla\theta) - \mathrm{Ae}(\nabla u_0 \theta)) : \nabla(v - u'_0) $$
$$ - \int_{\Omega_0} \mathrm{div}(\theta) \mathrm{Ae}(u_0) : e(v - u'_0) + \int_{\Gamma_{T_0}} \mathrm{Ae}(u_0) \mathrm{n} \cdot \nabla \theta^\top (v - u'_0) - \int_{\Gamma_0} (\theta \cdot \mathrm{n}) \mathrm{div}(\mathrm{Ae}(u_0)) \cdot (v - u'_0) \\ - \int_{\Gamma_{T_{0_N}^{u_0,g}}} p(\theta) \frac{u_{0_\tau}}{\|u_{0_\tau}\|} \cdot (v_\tau - u'_{0_\tau}) + \int_{\Gamma_{T_{0_S}^{u_0,g}}} p(\theta) \frac{\sigma_\tau(u_0)}{g} \cdot (v_\tau - u'_{0_\tau}), $$

for all $v \in \mathcal{K}_0 - \nabla\theta^\top u_0 - \nabla u_0 \theta$. Moreover, one has

$$ \int_{\Omega_0} \mathrm{div}(\mathrm{Ae}(u_0)) \theta^\top : \nabla \varphi = \int_{\Omega_0} \mathrm{div}(\mathrm{Ae}(u_0)) \cdot \nabla \varphi \theta = - \int_{\Omega_0} \mathrm{Ae}(u_0) : \nabla(\nabla \varphi \theta) + \int_{\Gamma_0} \mathrm{Ae}(u_0) \mathrm{n} \cdot \nabla \varphi \theta, $$

and also

$$ - \int_{\Omega_0} \mathrm{div}(\theta) \mathrm{Ae}(u_0) : e(\varphi) = \int_{\Omega_0} \theta \cdot \nabla(\mathrm{Ae}(u_0) : e(\varphi)) - \int_{\Gamma_0} \theta \cdot \mathrm{n}(\mathrm{Ae}(u_0) : e(\varphi)), $$

for all $\varphi \in C^\infty(\bar{\Omega}_0, \mathbb{R}^2)$. Therefore, using the equality

$$ ((\mathrm{Ae}(u_0)) \nabla\theta^\top + A(\nabla u_0 \nabla\theta) - \mathrm{Ae}(\nabla u_0 \theta)) : \nabla\varphi + \theta \cdot \nabla(\mathrm{Ae}(u_0) : e(\varphi)) - \mathrm{Ae}(u_0) : \nabla(\nabla\varphi\theta) = 0, $$

which holds true a.e. on $\Omega_0$, one deduces from the divergence formula that

$$ \int_{\Omega_0} (\mathrm{div}(\mathrm{Ae}(u_0)) \theta^\top + (\mathrm{Ae}(u_0)) \nabla\theta^\top + A(\nabla u_0 \nabla\theta) - \mathrm{Ae}(\nabla u_0 \theta)) : \nabla\varphi \\ - \int_{\Omega_0} \mathrm{div}(\theta) \mathrm{Ae}(u_0) : e(\varphi) + \int_{\Gamma_0} \nabla\theta (\mathrm{Ae}(u_0)\mathrm{n}) \cdot \varphi - \int_{\Gamma_0} (\theta \cdot \mathrm{n}) \mathrm{div}(\mathrm{Ae}(u_0)) \cdot \varphi - \int_{\Gamma_{T_{0_N}^{u_0,g}}} p(\theta) \frac{u_{0_\tau}}{\|u_{0_\tau}\|} \cdot \varphi_\tau \\ + \int_{\Gamma_{T_{0_S}^{u_0,g}}} p(\theta) \frac{\sigma_\tau(u_0)}{g} \cdot \varphi_\tau = \int_{\Gamma_0} \theta \cdot \mathrm{n}(-\mathrm{Ae}(u_0) : e(\varphi) - \mathrm{div}(\mathrm{Ae}(u_0)) \cdot \varphi) \\ + \int_{\Gamma_0} \nabla\varphi^\top (\mathrm{Ae}(u_0)\mathrm{n}) \cdot \theta + \nabla\theta(\mathrm{Ae}(u_0)\mathrm{n}) \cdot \varphi - \int_{\Gamma_{T_{0_N}^{u_0,g}}} p(\theta) \frac{u_{0_\tau}}{\|u_{0_\tau}\|} \cdot \varphi_\tau + \int_{\Gamma_{T_{0_S}^{u_0,g}}} p(\theta) \frac{\sigma_\tau(u_0)}{g} \cdot \varphi_\tau, $$

for all $\varphi \in C^\infty(\bar{\Omega}_0, \mathbb{R}^2)$. Furthermore, since $\Gamma_0$ is of class $C^3$ and $u_0 \in H^3(\Omega_0, \mathbb{R}^2)$, $\mathrm{Ae}(u_0)\mathrm{n}$ can be extended into a function defined in $\Omega_0$ such that $\mathrm{Ae}(u_0)\mathrm{n} \in H^2(\Omega_0, \mathbb{R}^2)$. Thus, it holds that $\mathrm{Ae}(u_0)\mathrm{n} \cdot \varphi \in W^{2,1}(\Omega_0, \mathbb{R}^2)$, for all $\varphi \in C^\infty(\bar{\Omega}_0, \mathbb{R}^2)$, and one can use Proposition B.2 to get that

$$ \int_{\Gamma_0} \theta \cdot \mathrm{n}(-\mathrm{Ae}(u_0) : e(\varphi) - \mathrm{div}(\mathrm{Ae}(u_0)) \cdot \varphi) + \int_{\Gamma_0} \nabla\varphi^\top (\mathrm{Ae}(u_0)\mathrm{n}) \cdot \theta + \nabla\theta (\mathrm{Ae}(u_0)\mathrm{n}) \cdot \varphi \\ - \int_{\Gamma_{T_{0_N}^{u_0,g}}} p(\theta) \frac{u_{0_\tau}}{\|u_{0_\tau}\|} \cdot \varphi_\tau + \int_{\Gamma_{T_{0_S}^{u_0,g}}} p(\theta) \frac{\sigma_\tau(u_0)}{g} \cdot \varphi_\tau \\ = \int_{\Gamma_0} \theta \cdot \mathrm{n}(-\mathrm{Ae}(u_0) : e(\varphi) - \mathrm{div}(\mathrm{Ae}(u_0)) \cdot \varphi + \partial_n (\mathrm{Ae}(u_0)\mathrm{n} \cdot \varphi) + H\mathrm{Ae}(u_0)\mathrm{n} \cdot \varphi) \\ - \int_{\Gamma_0} (\nabla(\mathrm{Ae}(u_0)\mathrm{n})\theta - \nabla\theta(\mathrm{Ae}(u_0)\mathrm{n}) + \mathrm{div}_\tau(\theta)\mathrm{Ae}(u_0)\mathrm{n}) \cdot \varphi - \int_{\Gamma_{T_{0_N}^{u_0,g}}} p(\theta) \frac{u_{0_\tau}}{\|u_{0_\tau}\|} \cdot \varphi_\tau \\ + \int_{\Gamma_{T_{0_S}^{u_0,g}}} p(\theta) \frac{\sigma_\tau(u_0)}{g} \cdot \varphi_\tau, $$
where $H$ is the mean curvature of $\Gamma_0$. By Proposition B.3 it follows that

$$ \int_{\Gamma_0} \theta \cdot n (-\operatorname{div}(\operatorname{Ae}(u_0)) + H \operatorname{Ae}(u_0)n) \cdot \varphi = \int_{\Gamma_0} \operatorname{Ae}(u_0) : \nabla_{\tau}(\varphi(\theta \cdot n)) - (\theta \cdot n) \partial_n (\operatorname{Ae}(u_0)) n \cdot \varphi, $$

for all $\varphi \in C^\infty(\bar{\Omega}_0, \mathbb{R}^2)$. Therefore, using the following two equalities

$$ \operatorname{Ae}(u_0) : \nabla_{\tau}(\varphi(\theta \cdot n)) = \theta \cdot n (\operatorname{Ae}(u_0) : \nabla_{\tau}\varphi) + \operatorname{Ae}(u_0) \nabla_{\tau}(\theta \cdot n) \cdot \varphi, \text{ a.e. on } \Gamma_0, $$

and

$$ \operatorname{Ae}(u_0) : \nabla_{\tau}\varphi = \operatorname{Ae}(u_0) : e(\varphi) - \nabla\varphi^{\top}(\operatorname{Ae}(u_0)n) \cdot n \text{ a.e on } \Gamma_0, $$

one gets

$$ \begin{aligned} & \int_{\Gamma_0} \theta \cdot n (-\operatorname{Ae}(u_0) : e(\varphi) - \operatorname{div}(\operatorname{Ae}(u_0)) \cdot \varphi + \partial_n (\operatorname{Ae}(u_0)n \cdot \varphi) + H \operatorname{Ae}(u_0)n \cdot \varphi) \\ & \quad - \int_{\Gamma_0} (\nabla(\operatorname{Ae}(u_0)n)\theta - \nabla\theta(\operatorname{Ae}(u_0)n) + \operatorname{div}_{\tau}(\theta)\operatorname{Ae}(u_0)n) \cdot \varphi - \int_{\Gamma_{T_{0_N}^{u_0,g}}} p(\theta) \frac{u_{0_{\tau}}}{\|u_{0_{\tau}}\|} \cdot \varphi_{\tau} \\ & \quad + \int_{\Gamma_{T_{0_S}^{u_0,g}}} p(\theta) \frac{\sigma_{\tau}(u_0)}{g} \cdot \varphi_{\tau} = \int_{\Gamma_0} (\theta \cdot n (\partial_n (\operatorname{Ae}(u_0)n) - \partial_n (\operatorname{Ae}(u_0))n) + \operatorname{Ae}(u_0) \nabla_{\tau}(\theta \cdot n)) \cdot \varphi \\ & \quad + \int_{\Gamma_0} (-\nabla(\operatorname{Ae}(u_0)n)\theta + (\nabla\theta - \operatorname{div}_{\tau}(\theta)I) \operatorname{Ae}(u_0)n) \cdot \varphi - \int_{\Gamma_{T_{0_N}^{u_0,g}}} p(\theta) \frac{u_{0_{\tau}}}{\|u_{0_{\tau}}\|} \cdot \varphi_{\tau} + \int_{\Gamma_{T_{0_S}^{u_0,g}}} p(\theta) \frac{\sigma_{\tau}(u_0)}{g} \cdot \varphi_{\tau}, \end{aligned} $$

and thus

$$ \begin{aligned} & \int_{\Omega_0} (\operatorname{div}(\operatorname{Ae}(u_0)) \theta^{\top} + (\operatorname{Ae}(u_0)) \nabla\theta^{\top} + A(\nabla u_0 \nabla\theta) - \operatorname{Ae}(\nabla u_0 \theta)) : \nabla\varphi - \int_{\Omega_0} \operatorname{div}(\theta) \operatorname{Ae}(u_0) : e(\varphi) \\ & + \int_{\Gamma_{T_0}} \operatorname{Ae}(u_0)n \cdot \nabla\theta^{\top}\varphi - \int_{\Gamma_0} (\theta \cdot n) \operatorname{div}(\operatorname{Ae}(u_0)) \cdot \varphi - \int_{\Gamma_{T_{0_N}^{u_0,g}}} p(\theta) \frac{u_{0_{\tau}}}{\|u_{0_{\tau}}\|} \cdot \varphi_{\tau} + \int_{\Gamma_{T_{0_S}^{u_0,g}}} p(\theta) \frac{\sigma_{\tau}(u_0)}{g} \cdot \varphi_{\tau} \\ & = \int_{\Gamma_0} (\theta \cdot n (\partial_n (\operatorname{Ae}(u_0)n) - \partial_n (\operatorname{Ae}(u_0))n) + \operatorname{Ae}(u_0) \nabla_{\tau}(\theta \cdot n)) \cdot \varphi \\ & + \int_{\Gamma_0} (-\nabla(\operatorname{Ae}(u_0)n)\theta + (\nabla\theta - \operatorname{div}_{\tau}(\theta)I) \operatorname{Ae}(u_0)n) \cdot \varphi - \int_{\Gamma_{T_{0_N}^{u_0,g}}} p(\theta) \frac{u_{0_{\tau}}}{\|u_{0_{\tau}}\|} \cdot \varphi_{\tau} + \int_{\Gamma_{T_{0_S}^{u_0,g}}} p(\theta) \frac{\sigma_{\tau}(u_0)}{g} \cdot \varphi_{\tau}, \end{aligned} $$

for all $\varphi \in C^\infty(\bar{\Omega}_0, \mathbb{R}^2)$. Finally, one deduces from the density of $C^\infty(\bar{\Omega}_0, \mathbb{R}^2)$ in $H^1(\Omega_0, \mathbb{R}^2)$ that

$$ \begin{aligned} \langle u'_0, v - u'_0 \rangle_{H_D^1(\Omega_0, \mathbb{R}^2)} \ge & \int_{\Gamma_{T_0}} (\theta \cdot n (\partial_n (\operatorname{Ae}(u_0)n) - \partial_n (\operatorname{Ae}(u_0))n) + \operatorname{Ae}(u_0) \nabla_{\tau}(\theta \cdot n)) \cdot (v - u'_0) \\ & + \int_{\Gamma_{T_0}} (-\nabla(\operatorname{Ae}(u_0)n)\theta + (\nabla\theta - \operatorname{div}_{\tau}(\theta)I) \operatorname{Ae}(u_0)n) \cdot (v - u'_0) - \int_{\Gamma_{T_{0_N}^{u_0,g}}} p(\theta) \frac{u_{0_{\tau}}}{\|u_{0_{\tau}}\|} \cdot (v_{\tau} - u'_{0\tau}) \\ & + \int_{\Gamma_{T_{0_S}^{u_0,g}}} p(\theta) \frac{\sigma_{\tau}(u_0)}{g} \cdot (v_{\tau} - u'_{0\tau}), \end{aligned} $$

for all $v \in K_0 - \nabla\theta^\top u_0 - \nabla u_0 \theta$, which corresponds to the weak formulation of the expected tangential Signorini problem (see [11, Section 2.1.2] for details). $\square$

**Remark 2.12.** Note that $\bar{u}'_0$ and $u'_0$ are not linear with respect to the direction $\theta$. This nonlinearity is standard in shape optimization for variational inequalities (see, e.g., [4, 23] or [34, Section 4]), and justifies the names of *directional* material and shape derivatives.
## 2.4 Shape gradient of the Tresca energy functional

Thanks to the characterization of the directional material and shape derivatives obtained in the previous section, we are now in a position to derive an expression of the shape gradient of the Tresca energy functional $\mathcal{J}$ at $\Omega_0$ in the direction $\theta$.

**Theorem 2.13.** Consider the framework of Proposition 2.5. Then the Tresca energy functional $\mathcal{J}$ admits a shape gradient at $\Omega_0$ in the direction $\theta$ given by

$$ \mathcal{J}'(\Omega_0)(\theta) = \int_{\Omega_0} \mathrm{div}(\theta) \frac{\mathrm{Ae}(u_0) : \mathrm{e}(u_0)}{2} - \int_{\Omega_0} \mathrm{div}(\mathrm{Ae}(u_0)) \cdot \nabla u_0 \theta - \int_{\Omega_0} \mathrm{Ae}(u_0) : \nabla u_0 \nabla \theta \\ - \int_{\Gamma_{T_0}} \theta \cdot \mathrm{n}(f \cdot u_0) - \langle \mathrm{Ae}(u_0)\mathrm{n}, \nabla \theta^\top u_0 \rangle_{\mathrm{H}^{-1/2}(\Gamma_0, \mathbb{R}^2) \times \mathrm{H}^{1/2}(\Gamma_0, \mathbb{R}^2)} \\ + \int_{\Gamma_{T_0}^{u_0,g}} (\nabla g \cdot \theta + g(\mathrm{div}_\tau(\theta) - \nabla \theta \tau \cdot \tau)) \|u_{0\tau}\| $$

*Proof.* By taking $v = u_t$ in the variational inequality satisfied by $u_t$, one can obtain that

$$ \mathcal{J}(\Omega_t) = -\frac{1}{2} \int_{\Omega_t} \mathrm{Ae}(u_t) : \mathrm{e}(u_t). $$

Following the usual strategy developed in (smooth) shape optimization literature (see, e.g., [7, 22]) to compute the shape gradient of $\mathcal{J}$ at $\Omega_0$ in the direction $\theta$, one gets

$$ \mathcal{J}'(\Omega_0)(\theta) = -\frac{1}{2} \int_{\Omega_0} \mathrm{div}(\theta) \mathrm{Ae}(u_0) : \mathrm{e}(u_0) + \int_{\Omega_0} \mathrm{Ae}(u_0) : \nabla u_0 \nabla \theta - \langle \bar{u}'_0, u_0 \rangle_{\mathrm{H}^1_{\mathrm{D}}(\Omega_0, \mathbb{R}^2)}. $$

Moreover, from the variational inequality satisfied by $\bar{u}'_0$ (see (2.7)), the divergence formula (see Proposition B.1) and since $\bar{u}'_0 \pm u_0 \in \mathcal{K}_0 - \nabla \theta^\top u_0$, it follows that

$$ \langle \bar{u}'_0, u_0 \rangle_{\mathrm{H}^1_{\mathrm{D}}(\Omega_0, \mathbb{R}^2)} = \int_{\Omega_0} (\mathrm{div}(\mathrm{Ae}(u_0)) \theta^\top + (\mathrm{Ae}(u_0)) \nabla \theta^\top + \mathrm{A}(\nabla u_0 \nabla \theta) - \mathrm{div}(\theta) \mathrm{Ae}(u_0)) : \nabla u_0 \\ + \int_{\Gamma_{T_0}} \theta \cdot \mathrm{n}(f \cdot u_0) + \langle \mathrm{Ae}(u_0)\mathrm{n}, \nabla \theta^\top u_0 \rangle_{\mathrm{H}^{-1/2}(\Gamma_0, \mathbb{R}^2) \times \mathrm{H}^{1/2}(\Gamma_0, \mathbb{R}^2)} - \int_{\Gamma_{T_0}^{u_0,g}} p(\theta) \|u_{0\tau}\|. $$

Then, since $u_{0\tau} = 0$ a.e. on $\Gamma_{T_0}^{u_0,g}$ and using the equality $\mathrm{div}(\mathrm{Ae}(u_0)) \theta^\top : \nabla u_0 = \mathrm{div}(\mathrm{Ae}(u_0)) \cdot \nabla u_0 \theta$ which holds true a.e. on $\Omega_0$, we conclude the proof. $\square$

As we did for the directional material derivative, the presentation of Theorem 2.13 can be improved under additional assumptions.

**Corollary 2.14.** Consider the framework of Proposition 2.5 with the additional assumptions that $u_0 \in \mathrm{H}^3(\Omega_0, \mathbb{R}^2)$, $\Gamma_0$ is of class $\mathcal{C}^3$ and almost every point of $\Gamma_{T_0}$ belongs to the relative interior $\mathrm{int}\Gamma_0(\Gamma_{T_0})$. Then the Tresca energy functional $\mathcal{J}$ admits a shape gradient at $\Omega_0$ in the direction $\theta$ given by

$$ \mathcal{J}'(\Omega_0)(\theta) = \int_{\Gamma_{T_0}} \theta \cdot \mathrm{n} \left( \frac{\mathrm{Ae}(u_0) : \mathrm{e}(u_0)}{2} - f \cdot u_0 - \sigma_\tau(u_0) \cdot \partial_n(u_0) + \|u_{0\tau}\| (Hg + \partial_n g) \right) \\ + \int_{\Gamma_{T_0}} u_{0n} \sigma_\tau(u_0) \cdot \tau (\nabla \tau \theta_\tau - \nabla \theta \tau) \cdot \mathrm{n}, $$

where $H$ is the mean curvature of $\Gamma_0$.
*Proof.* Since $u_0 \in H^2(\Omega_0, \mathbb{R}^2)$, it follows from Theorem 2.13 that

$$
\begin{aligned}
\mathcal{J}'(\Omega_0)(\theta) = & -\frac{1}{2} \int_{\Omega_0} \theta \cdot \nabla (\mathrm{Ae}(u_0) : \mathrm{e}(u_0)) + \int_{\Gamma_0} \theta \cdot n \frac{\mathrm{Ae}(u_0) : \mathrm{e}(u_0)}{2} + \int_{\Omega_0} \mathrm{Ae}(u_0) : \mathrm{e}(\nabla u_0 \theta) \\
& - \int_{\Gamma_0} \mathrm{Ae}(u_0) n \cdot \nabla u_0 \theta - \int_{\Omega_0} \mathrm{Ae}(u_0) : \nabla u_0 \nabla \theta - \int_{\Gamma_{T_0}} \theta \cdot n(f \cdot u_0) - \int_{\Gamma_{T_0}} \mathrm{Ae}(u_0) n \cdot \nabla \theta^\top u_0 \\
& + \int_{\Gamma_{T_{0,N}^{u_0,g}}} (\nabla g \cdot \theta + g(\mathrm{div}_\tau(\theta) - \nabla \theta \tau \cdot \tau)) \|u_{0\tau}\|.
\end{aligned}
$$

Moreover, since

$$ \mathrm{Ae}(u_0) : \mathrm{e}(\nabla u_0 \theta) = \mathrm{Ae}(u_0) : \nabla u_0 \nabla \theta + \frac{1}{2} \theta \cdot \nabla (\mathrm{Ae}(u_0) : \mathrm{e}(u_0)) \quad \text{a.e. on } \Omega_0, $$

one deduces that

$$
\begin{aligned}
\mathcal{J}'(\Omega_0)(\theta) = & \int_{\Gamma_0} \theta \cdot n \left( \frac{\mathrm{Ae}(u_0) : \mathrm{e}(u_0)}{2} \right) - \int_{\Gamma_0} \mathrm{Ae}(u_0) n \cdot \nabla u_0 \theta - \int_{\Gamma_{T_0}} \theta \cdot n(f \cdot u_0) - \int_{\Gamma_0} \mathrm{Ae}(u_0) n \cdot \nabla \theta^\top u_0 \\
& + \int_{\Gamma_{T_{0,N}^{u_0,g}}} (\nabla g \cdot \theta + g(\mathrm{div}_\tau(\theta) - \nabla \theta \tau \cdot \tau)) \|u_{0\tau}\|.
\end{aligned}
$$

Furthermore, since almost every point of $\Gamma_{T_0}$ belongs to $\mathrm{int}_{\Gamma_0}(\Gamma_{T_0})$, it follows that $u_0$ is a strong solution to the Tresca friction problem (see [11, Definition 2.11 and Proposition 2.13]). Thus, from the Tresca friction law, one has $\|u_{0\tau}\| = -\frac{\sigma_{\tau}(u_0) \cdot u_{0\tau}}{g}$ a.e. on $\Gamma_{T_0}$ and, since $u_{0\tau} = 0$ on $\Gamma_{T_{0D}^{u_0,g}} \cup \Gamma_{T_{0S}^{u_0,g}}$ and $\theta = 0$ on $\Gamma_D$, one gets that

$$
\begin{aligned}
\mathcal{J}'(\Omega_0)(\theta) = & \int_{\Gamma_0} \theta \cdot n \left( \frac{\mathrm{Ae}(u_0) : \mathrm{e}(u_0)}{2} - f \cdot u_0 \right) - \int_{\Gamma_0} \mathrm{Ae}(u_0) n \cdot \nabla u_0 \theta - \int_{\Gamma_0} \mathrm{Ae}(u_0) n \cdot \nabla \theta^\top u_0 \\
& - \int_{\Gamma_0} \left( g \sigma_\tau(u_0) \cdot u_{0\tau} \frac{\nabla g}{g^2} \cdot \theta + \sigma_\tau(u_0) \cdot u_{0\tau} \mathrm{div}_\tau(\theta) \right) + \int_{\Gamma_0} \sigma_\tau(u_0) \cdot u_{0\tau} \nabla \theta \tau \cdot \tau. \quad (2.9)
\end{aligned}
$$

Moreover, since $\Gamma_0$ is of class $C^3$, $n \in C^2(\Gamma_0, \mathbb{R}^2)$ can be extended over $\mathbb{R}^2$ such that $n \in C^2(\mathbb{R}^2, \mathbb{R}^2)$ and $\|n\| = 1$ over $\mathbb{R}^2$ (see, e.g., [22, Chapter 5 Section 5.4]). It follows that $\tau \in C^2(\Gamma_0, \mathbb{R}^2)$ can also be extended over $\mathbb{R}^2$ such that $\tau \in C^2(\mathbb{R}^2, \mathbb{R}^2)$ and $\|\tau\| = 1$ over $\mathbb{R}^2$. Moreover, since $u_0 \in H^3(\Omega_0, \mathbb{R}^2)$, the shear stress $\sigma_\tau(u_0) = \mathrm{Ae}(u_0)n - (\mathrm{Ae}(u_0)n \cdot n)n$ can be extended into a function defined over $\Omega_0$ such that $\sigma_\tau(u_0) \in H^2(\Omega_0, \mathbb{R}^2)$. Thus $\sigma_\tau(u_0) \cdot \tau \in H^2(\Omega_0, \mathbb{R}^2)$, $u_0 \cdot \tau \in H^2(\Omega_0, \mathbb{R}^2)$ and one deduces that

$$
\begin{aligned}
\mathcal{J}'(\Omega_0)(\theta) = & \int_{\Gamma_0} \theta \cdot n \left( \frac{\mathrm{Ae}(u_0) : \mathrm{e}(u_0)}{2} - f \cdot u_0 \right) - \int_{\Gamma_0} \mathrm{Ae}(u_0) n \cdot \nabla u_0 \theta - \int_{\Gamma_0} \mathrm{Ae}(u_0) n \cdot \nabla \theta^\top u_0 \\
& - \int_{\Gamma_0} \left( \nabla((\sigma_\tau(u_0) \cdot \tau) u_0 \cdot \tau) \cdot \theta + \sigma_\tau(u_0) \cdot u_{0\tau} \mathrm{div}_\tau(\theta) \right) + \int_{\Gamma_0} \sigma_\tau(u_0) \cdot \tau \nabla(u_0 \cdot \tau) \cdot \theta \\
& + \int_{\Gamma_0} g u_0 \cdot \tau \nabla \left( \frac{\sigma_\tau(u_0) \cdot \tau}{g} \right) \cdot \theta + \int_{\Gamma_0} \sigma_\tau(u_0) \cdot u_{0\tau} \nabla \theta \tau \cdot \tau,
\end{aligned}
$$

and, since $(\sigma_\tau(u_0) \cdot \tau) u_0 \cdot \tau \in W^{1,2}(\Omega_0, \mathbb{R})$, one can apply Proposition B.2 to get that

$$
\mathcal{J}'(\Omega_0)(\theta) = \int_{\Gamma_0} \theta \cdot n \left( \frac{\mathrm{Ae}(u_0) : \mathrm{e}(u_0)}{2} - f \cdot u_0 - H \sigma_{\tau}(u_0) \cdot u_{0\tau} - \partial_n ((\sigma_{\tau}(u_0) \cdot \tau) u_0 \cdot \tau) \right)
$$
$$ - \int_{\Gamma_0} \mathrm{Ae}(u_0) \mathrm{n} \cdot \nabla u_0 \theta - \int_{\Gamma_0} \mathrm{Ae}(u_0) \mathrm{n} \cdot \nabla \theta^\top u_0 + \int_{\Gamma_0} \sigma_\tau(u_0) \cdot \tau \nabla(u_0 \cdot \tau) \cdot \theta \\ + \int_{\Gamma_0} g u_0 \cdot \tau \nabla \left( \frac{\sigma_\tau(u_0) \cdot \tau}{g} \right) \cdot \theta + \int_{\Gamma_0} \sigma_\tau(u_0) \cdot u_0 \tau \nabla \theta \tau \cdot \tau. $$

Moreover, since $\sigma_n(u_0) = 0$ a.e. on $\Gamma_{T_0}$, note that

$$ - \int_{\Gamma_0} \mathrm{Ae}(u_0) \mathrm{n} \cdot \nabla \theta^\top u_0 + \int_{\Gamma_0} \sigma_\tau(u_0) \cdot u_{0\tau} \nabla \theta \tau \cdot \tau = - \int_{\Gamma_{T_0}} (\sigma_\tau(u_0) \cdot \tau) \nabla \theta \tau \cdot u_0 \\ + \int_{\Gamma_{T_0}} \sigma_\tau(u_0) \cdot u_{0\tau} \nabla \theta \tau \cdot \tau = - \int_{\Gamma_{T_0}} \sigma_\tau(u_0) \cdot u_{0\tau} \nabla \theta \tau \cdot \tau - \int_{\Gamma_{T_0}} (\sigma_\tau(u_0) \cdot \tau) u_{0n} \nabla \theta \tau \cdot \mathrm{n} \\ + \int_{\Gamma_{T_0}} \sigma_\tau(u_0) \cdot u_{0\tau} \nabla \theta \tau \cdot \tau = - \int_{\Gamma_{T_0}} (\sigma_\tau(u_0) \cdot \tau) u_{0n} \nabla \theta \tau \cdot \mathrm{n}, $$

and also that

$$ \int_{\Gamma_0} \theta \cdot \mathrm{n} (-\partial_n ((\sigma_\tau(u_0) \cdot \tau) u_0 \cdot \tau)) - \int_{\Gamma_0} \mathrm{Ae}(u_0) \mathrm{n} \cdot \nabla u_0 \theta + \int_{\Gamma_0} \sigma_\tau(u_0) \cdot \tau \nabla(u_0 \cdot \tau) \cdot \theta \\ = \int_{\Gamma_{T_0}} \theta \cdot \mathrm{n} (-\partial_n(u_0 \cdot \tau) \sigma_\tau(u_0) \cdot \tau - \partial_n(\sigma_\tau(u_0) \cdot \tau) u_0 \cdot \tau) - \int_{\Gamma_{T_0}} (\sigma_\tau(u_0) \cdot \tau) \tau \cdot \nabla u_0 \theta \\ + \int_{\Gamma_{T_0}} \theta \cdot \mathrm{n} (\sigma_\tau(u_0) \cdot \tau) \nabla(u_0 \cdot \tau) \cdot \mathrm{n} + \int_{\Gamma_{T_0}} (\sigma_\tau(u_0) \cdot \tau) \nabla(u_0 \cdot \tau) \cdot \theta_\tau \\ = \int_{\Gamma_{T_0}} \theta \cdot \mathrm{n} (-\partial_n(\sigma_\tau(u_0) \cdot \tau) u_0 \cdot \tau) - \int_{\Gamma_{T_0}} \theta \cdot \mathrm{n} (\sigma_\tau(u_0) \cdot \tau) \tau \cdot \nabla u_0 \mathrm{n} - \int_{\Gamma_{T_0}} (\sigma_\tau(u_0) \cdot \tau) \theta_\tau \cdot \nabla u_0 \tau \\ + \int_{\Gamma_{T_0}} (\sigma_\tau(u_0) \cdot \tau) \nabla u_0^\top \tau \cdot \theta_\tau + \int_{\Gamma_{T_0}} (\sigma_\tau(u_0) \cdot \tau) \nabla \tau^\top u_0 \cdot \theta_\tau \\ = - \int_{\Gamma_{T_0}} \theta \cdot \mathrm{n} (\partial_n(\sigma_\tau(u_0) \cdot \tau) u_0 \cdot \tau + \sigma_\tau(u_0) \cdot \partial_n(u_0)) + \int_{\Gamma_{T_0}} (\sigma_\tau(u_0) \cdot \tau) u_{0n} \nabla \tau^\top \mathrm{n} \cdot \theta_\tau, $$

since $||\tau||=1$ on $\mathbb{R}^2$, thus $(\nabla\tau)\tau \cdot \tau = 0$ on $\Gamma_{T_0}$. Hence one has

$$ \mathcal{J}'(\Omega_0)(\theta) = \int_{\Gamma_{T_0}} g u_0 \cdot \tau \nabla \left( \frac{\sigma_\tau(u_0) \cdot \tau}{g} \right) \cdot \theta + \int_{\Gamma_{T_0}} u_{0n} \sigma_\tau(u_0) \cdot \tau (\nabla \tau \theta_\tau - \nabla \theta \tau) \cdot \mathrm{n} \\ + \int_{\Gamma_{T_0}} \theta \cdot \mathrm{n} \left( \frac{\mathrm{Ae}(u_0):\mathrm{e}(u_0)}{2} - f \cdot u_0 - H \sigma_\tau(u_0) \cdot u_{0\tau} - u_0 \cdot \tau \partial_n (\sigma_\tau(u_0) \cdot \tau) - \sigma_\tau(u_0) \cdot \partial_n(u_0) \right). $$

Now let us focus on the first term. Since $u_{0\tau} = 0$ on $\Gamma_{T_0^{u_{0,g}}} \cup \Gamma_{T_0^{u_{0,g}}}$, we have

$$ \int_{\Gamma_{T_0}} g u_0 \cdot \tau \nabla \left( \frac{\sigma_\tau(u_0) \cdot \tau}{g} \right) \cdot \theta = \int_{\Gamma_{T_0^{u_{0,g}}}} g u_0 \cdot \tau \nabla \left( \frac{\sigma_\tau(u_0) \cdot \tau}{g} \right) \cdot \theta. $$

Let us introduce two disjoint subsets of $\Gamma_{T_0^{u_{0,g}}}$ given by

$$ \Gamma_{T_0^{u_{0,g}}} := \{s \in \Gamma_{T_0} \mid u_0(s) \cdot \tau(s) > 0\} \quad \text{and} \quad \Gamma_{T_0^{u_{0,g}}} := \{s \in \Gamma_{T_0} \mid u_0(s) \cdot \tau(s) < 0\}. $$

It follows that $\Gamma_{T_0^{u_{0,g}}} = \Gamma_{T_0^{u_{0,g}}} \cup \Gamma_{T_0^{u_{0,g}}}$, with $\sigma_\tau(u_0) \cdot \tau = -g$ a.e. on $\Gamma_{T_0^{u_{0,g}}}$, and $\sigma_\tau(u_0) \cdot \tau = g$ a.e. on $\Gamma_{T_0^{u_{0,g}}}$. Moreover, since $u_0 \in \mathbb{H}^3(\Omega, \mathbb{R}^2)$, we get from Sobolev embeddings (see, e.g., [1,
Chapter 4]) that $u_0$ is continuous over $\Gamma_{T_0}$, thus $\Gamma_{T_0^{u_0,g}}$ and $\Gamma_{T_0^{u_0,g}}$ are open subsets of $\Gamma_{T_0}$. Hence $\nabla_{\tau}(\frac{\sigma_{\tau}(u_0)\cdot\tau}{g}) = 0$ a.e. on $\Gamma_{T_0^{u_0,g}} \cup \Gamma_{T_0^{u_0,g}}$, and one deduces that

$$ \int_{\Gamma_{T_0}} gu_0 \cdot \tau \nabla \left( \frac{\sigma_{\tau}(u_0) \cdot \tau}{g} \right) \cdot \theta = \int_{\Gamma_{T_0}} \theta \cdot n \left( gu_0 \cdot \tau \nabla \left( \frac{\sigma_{\tau}(u_0) \cdot \tau}{g} \right) \cdot n \right), $$

that is

$$ \int_{\Gamma_{T_0}} gu_0 \cdot \tau \nabla \left( \frac{\sigma_{\tau}(u_0) \cdot \tau}{g} \right) \cdot \theta = \int_{\Gamma_{T_0}} \theta \cdot n \left( u_0 \cdot \tau \partial_n (\sigma_{\tau}(u_0) \cdot \tau) - \frac{\sigma_{\tau}(u_0) \cdot u_{0\tau}}{g} \partial_n g \right). $$

Then, using the Tresca friction law, one has $\sigma_{\tau}(u_0) \cdot u_{0\tau} = -g||u_{0\tau}||$ a.e. on $\Gamma_{T_0}$, which concludes the proof. $\square$

**Remark 2.15.** Under the weaker condition $u_0 \in H^2(\Omega_0, \mathbb{R}^2)$, one can follow the proof of Corollary 2.14 and obtain that the shape gradient of $\mathcal{J}$ is given by Equality (2.9).

# 3 Numerical illustration

In this section our objective is to numerically solve a toy example of the shape optimization problem (1.1), by making use of the theoretical results established in this work. The numerical simulations have been performed using Freefem++ software [20] with P1-finite elements and standard affine mesh. We could use the expression of the shape gradient of $\mathcal{J}$ obtained in Theorem 2.13 but, in order to simplify the computations, we chose to use the expression provided in Corollary 2.14 under additional assumptions that we assumed to be true at each iteration.

## 3.1 Numerical methodology

Consider an initial shape $\Omega_0 \in \mathcal{U}_{\text{ref}}$. Note that Corollary 2.14 allows to exhibit a descent direction $\theta_0$ of the Tresca energy functional $\mathcal{J}$ at $\Omega_0$, by finding the unique solution $\theta_0 \in H_D^1(\Omega_0, \mathbb{R}^2)$ to the variational equality

$$ \int_{\Omega_0} (\nabla \theta_0 : \nabla \theta + \theta_0 \cdot \theta) = -\mathcal{J}'(\Omega_0)(\theta), \quad \forall \theta \in H_D^1(\Omega_0, \mathbb{R}^2), $$

since it satisfies $\mathcal{J}'(\Omega_0)(\theta_0) = -\int_{\Omega_0} (||\nabla \theta_0||^2 + ||\theta_0||^2) \le 0$.

In order to numerically solve the shape optimization problem (1.1) on a given example, we have to deal with the volume constraint $|\Omega| = |\Omega_{\text{ref}}| > 0$. For this purpose, the Uzawa algorithm (see, e.g., [7, Chapter 3]) is used, and one refers to [4, Section 4] for methodological details.

Let us mention that the Tresca friction problem is numerically solved using an adaptation of iterative switching algorithms (see [5]). This algorithm operates by checking at each iteration if the Tresca boundary conditions are satisfied and, if they are not, by imposing them and restarting the computation (see [3, Appendix C p.25] for detailed explanations). We also precise that, for all $j \in \mathbb{N}^*$, the difference between the Tresca energy functional $\mathcal{J}$ at the iteration $20 \times j$ and at the iteration $20 \times (j - 1)$ is computed. The smallness of this difference is used as a stopping criterion for the algorithm. Finally the curvature term $H$ is numerically computed by extending the normal $n$ into a function $\tilde{n}$ which is defined on the whole domain $\Omega_0$. Then the curvature is given by $H = \text{div}(\tilde{n}) - \nabla(\tilde{n}) \cdot n$ (see, e.g., [22, Proposition 5.4.8]).
## 3.2 A toy example and numerical results

In this section, let $f \in H^1(\mathbb{R}^2, \mathbb{R}^2)$ be defined by

$$f: \begin{array}{rcl} \mathbb{R}^2 & \to & \mathbb{R}^2 \\ (x, y) & \mapsto & f(x, y) := (-5x \exp(x) \quad 0.6 \exp(x^2)) \eta(x, y), \end{array}$$

and $g \in H^2(\mathbb{R}^2, \mathbb{R})$ be defined by

$$g: \begin{array}{rcl} \mathbb{R}^2 & \to & \mathbb{R} \\ (x, y) & \mapsto & g(x, y) := (1 + \sin(-y\frac{\pi}{2}) + 10^{-3}) \eta(x, y), \end{array}$$

where $\eta \in C^\infty(\mathbb{R}^2, \mathbb{R})$ is a cut-off function chosen appropriately so that $f$ belongs to $H^1(\mathbb{R}^2, \mathbb{R}^2)$, $g \in H^2(\mathbb{R}^2, \mathbb{R})$ and $g > 0$ on $\mathbb{R}^2$. The reference shape $\Omega_{\text{ref}} \subset \mathbb{R}^2$ is an ellipse centered at $(0, 0) \in \mathbb{R}^2$, with semi-major axis $a = 1.1$ and semi-minor axis $b = 1/a$, and the fixed part $\Gamma_D$ is given by

$$\Gamma_D := \left\{ (a \cos \gamma, b \sin \gamma) \in \Gamma_{\text{ref}} \mid \gamma \in \left[ \frac{2\pi}{3}, \frac{4\pi}{3} \right] \cup \left[ \frac{5\pi}{3}, \frac{7\pi}{3} \right] \right\}.$$

We refer to Figure 1. The volume constraint is $|\Omega_{\text{ref}}| = \pi$ and the initial shape is $\Omega_0 := \Omega_{\text{ref}}$.

In the sequel we consider that, for all $\Omega \in \mathcal{U}_{\text{ref}}$, the Cauchy stress tensor $\sigma$, defined by $\sigma(v) := \text{Ae}(v)$ for all $v \in H_D^1(\Omega, \mathbb{R}^2)$, satisfies

$$\sigma(v) = 2\mu e(v) + \lambda \text{tr}(e(v)) I,$$

for all $v \in H_D^1(\Omega, \mathbb{R}^2)$, where $\text{tr}(e(v))$ is the trace of the matrix $e(v)$, and $\mu \ge 0, \lambda \ge 0$ are Lamé parameters (see, e.g., [31]). From a physical point of view, this assumption corresponds to *isotropic* elastic solids. In the sequel we consider the arbitrary data $\mu = 0.5$ and $\lambda = 0$.

We present now the numerical results obtained for this toy example using the numerical methodology described in Section 3.1.

<figure><img src="image_3.png" /><figcaption>Figure 1: $\Omega_{\text{ref}}$ and its boundary $\Gamma_{\text{ref}} = \Gamma_D \cup \Gamma_{T_{\text{ref}}}$.</figcaption></figure>

In Figure 2 is represented the initial shape (left) and the shape which numerically solves Problem (1.1) (right). On top are the vector values of the solution $u$ to the Tresca friction problem ($\text{TP}_{\Omega}$). On the initial shape, note that, on the bottom blue boundary, the norm of the shear stress is strictly inferior at the friction threshold $g$, thus $u_{\tau} = 0$, while the top black boundary shows some points where the norm of the shear stress reaches the friction threshold.

Figure 3 shows the values of $J$ (left) and the volume $|\Omega|$ of the shape (right) with respect to the iterations. We observe that $J$ is lower at the final shape, than at the initial shape, with some oscillations due to the Lagrange multiplier in order to satisfy the volume constraint.
<figure><img src="image_4.png" /><figcaption>Figure 2: Initial shape (left) and the shape minimizing $\mathcal{J}$ under the volume constraint $|\Omega| = \pi$ (right).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Values of the energy functional (left) and of the volume (right) with respect to the iterations.</figcaption></figure>

## A Reminders on twice epi-differentiability

For notions and results recalled in this appendix, we refer to standard references from nons-mooth analysis literature such as [12, 26, 28] and [30, Chapter 12]. In what follows, $(\mathcal{H}, \langle \cdot, \cdot \rangle_{\mathcal{H}})$ stands for a general real Hilbert space. The domain and the epigraph of an extended real-valued function $\psi: \mathcal{H} \to \mathbb{R} \cup \{\pm\infty\}$ are respectively defined by

$$ \mathrm{dom}(\psi) := \{x \in \mathcal{H} \mid \psi(x) < +\infty\} \quad \text{and} \quad \mathrm{epi}(\psi) := \{(x, \beta) \in \mathcal{H} \times \mathbb{R} \mid \psi(x) \le \beta\}. $$

Recall that $\psi$ is said to be proper if $\mathrm{dom}(\psi) \ne \emptyset$ and $\psi(x) > -\infty$ for all $x \in \mathcal{H}$, and that $\psi$ is convex (resp. lower semi-continuous) if and only if $\mathrm{epi}(\psi)$ is a convex (resp. closed) subset of $\mathcal{H} \times \mathbb{R}$. When $\psi$ is proper, we denote by $\partial\psi: \mathcal{H} \rightrightarrows \mathcal{H}$ its convex subdifferential operator, defined by

$$ \partial\psi(x) := \{y \in \mathcal{H} \mid \forall z \in \mathcal{H}, \langle y, z - x \rangle_{\mathcal{H}} \le \psi(z) - \psi(x)\}, $$

when $x \in \mathrm{dom}(\psi)$, and by $\partial\psi(x) := \emptyset$ when $x \notin \mathrm{dom}(\psi)$. The notion of proximal operator has been introduced by J.J. Moreau in 1965 (see [27]) as follows.

**Definition A.1 (Proximal operator).** The proximal operator associated with a proper, lower semi-continuous and convex function $\psi: \mathcal{H} \to \mathbb{R} \cup \{+\infty\}$ is the map $\mathrm{prox}_{\psi}: \mathcal{H} \to \mathcal{H}$ defined by

$$ \mathrm{prox}_{\psi}(x) := \underset{y \in \mathcal{H}}{\operatorname{argmin}} \left[ \psi(y) + \frac{1}{2} \|y - x\|_{\mathcal{H}}^2 \right] = (\mathrm{id} + \partial\psi)^{-1}(x), $$

for all $x \in \mathcal{H}$, where $\mathrm{id}: \mathcal{H} \to \mathcal{H}$ stands for the identity operator.
Recall that, if $\psi: \mathcal{H} \to \mathbb{R} \cup \{+\infty\}$ is a proper, lower semi-continuous and convex function, then its subdifferential $\partial\psi$ is maximal monotone (see, e.g., [28]), and thus its proximal operator $\text{prox}_\psi: \mathcal{H} \to \mathcal{H}$ is well-defined, single-valued and nonexpansive, i.e. Lipschitz continuous with modulus 1 (see, e.g., [12, Chapter II]).

The notion of twice epi-differentiability introduced by R.T. Rockafellar in 1985 (see [29]) is defined as the Mosco epi-convergence of second-order difference quotient functions. In what follows, we provide reminders and backgrounds on these notions for the reader's convenience. For more details, we refer to [30, Chapter 7, Section B] for the finite-dimensional case and to [15] for the infinite-dimensional case. In the sequel, the strong (resp. weak) convergence of a sequence in $\mathcal{H}$ will be denoted by $\to$ (resp. $\dasharrow$) and all limits with respect to $t$ will be considered for $t \to 0^+$.

**Definition A.2 (Mosco convergence).** The outer, weak-outer, inner and weak-inner limits of a parameterized family $(S_t)_{t>0}$ of subsets of $\mathcal{H}$ are respectively defined by

$$ \limsup S_t := \{x \in \mathcal{H} \mid \exists (t_n)_{n \in \mathbb{N}} \to 0^+, \exists (x_n)_{n \in \mathbb{N}} \to x, \forall n \in \mathbb{N}, x_n \in S_{t_n}\}, $$

$$ \text{w-lim sup } S_t := \{x \in \mathcal{H} \mid \exists (t_n)_{n \in \mathbb{N}} \to 0^+, \exists (x_n)_{n \in \mathbb{N}} \to x, \forall n \in \mathbb{N}, x_n \in S_{t_n}\}, $$

$$ \liminf S_t := \{x \in \mathcal{H} \mid \forall (t_n)_{n \in \mathbb{N}} \to 0^+, \exists (x_n)_{n \in \mathbb{N}} \to x, \exists N \in \mathbb{N}, \forall n \ge N, x_n \in S_{t_n}\}, $$

$$ \text{w-lim inf } S_t := \{x \in \mathcal{H} \mid \forall (t_n)_{n \in \mathbb{N}} \to 0^+, \exists (x_n)_{n \in \mathbb{N}} \to x, \exists N \in \mathbb{N}, \forall n \ge N, x_n \in S_{t_n}\}. $$

The family $(S_t)_{t>0}$ is said to be Mosco convergent if $\text{w-lim sup } S_t \subset \liminf S_t$. In that case, all the previous limits are equal and we write

$$ \text{M-lim } S_t := \liminf S_t = \limsup S_t = \text{w-lim inf } S_t = \text{w-lim sup } S_t. $$

**Definition A.3 (Mosco epi-convergence).** Let $(\psi_t)_{t>0}$ be a parameterized family of functions $\psi_t: \mathcal{H} \to \mathbb{R} \cup \{\pm\infty\}$ for all $t > 0$. We say that $(\psi_t)_{t>0}$ is Mosco epi-convergent if $(\text{epi}(\psi_t))_{t>0}$ is Mosco convergent in $\mathcal{H} \times \mathbb{R}$. Then we denote by ME-lim $\psi_t: \mathcal{H} \to \mathbb{R} \cup \{\pm\infty\}$ the function characterized by its epigraph epi (ME-lim $\psi_t$) := M-lim epi $(\psi_t)$ and we say that $(\psi_t)_{t>0}$ Mosco epi-converges to ME-lim $\psi_t$.

The notion of twice epi-differentiability was originally introduced in [29] for nonparameterized convex functions. However, the framework of the present paper requires an extended version to parameterized convex functions which has been developed in [2]. To provide reminders on this extended notion, when considering a function $\Psi: \mathbb{R}_+ \times \mathcal{H} \to \mathbb{R} \cup \{+\infty\}$ such that, for all $t \ge 0$, $\Psi(t, \cdot): \mathcal{H} \to \mathbb{R} \cup \{+\infty\}$ is a proper function, we will make use of the following two notations: $\partial\Psi(0, \cdot)(x)$ stands for the convex subdifferential operator at $x \in \mathcal{H}$ of the function $\Psi(0, \cdot)$, and, for each $t \ge 0$, $\Psi^{-1}(t, \mathbb{R}) := \{x \in \mathcal{H} \mid \Psi(t, x) \in \mathbb{R}\}$ and $\Psi^{-1}(\cdot, \mathbb{R}) := \cap_{t \ge 0} \Psi^{-1}(t, \mathbb{R})$.

**Definition A.4 (Twice epi-differentiability depending on a parameter).** Let $\Psi: \mathbb{R}_+ \times \mathcal{H} \to \mathbb{R} \cup \{+\infty\}$ be a function such that, for all $t \ge 0$, $\Psi(t, \cdot): \mathcal{H} \to \mathbb{R} \cup \{+\infty\}$ is a proper lower semi-continuous convex function. Then $\Psi$ is said to be twice epi-differentiable at $x \in \Psi^{-1}(\cdot, \mathbb{R})$ for $y \in \partial\Psi(0, \cdot)(x)$ if the family of second-order difference quotient functions $(\Delta_t^2 \Psi(x|y))_{t>0}$ defined by

$$ \begin{array}{rcl} \Delta_t^2 \Psi(x|y): & \mathcal{H} & \to \mathbb{R} \cup \{+\infty\} \\ & z & \mapsto \Delta_t^2 \Psi(x|y)(z) := \frac{\Psi(t, x + tz) - \Psi(t, x) - t \langle y, z \rangle_{\mathcal{H}}}{t^2}, \end{array} $$

for all $t > 0$, is Mosco epi-convergent. In that case, we denote by

$$ D_e^2 \Psi(x|y) := \text{ME-lim} \Delta_t^2 \Psi(x|y), $$

which is called the second-order epi-derivative of $\Psi$ at $x$ for $y$.
The next proposition (which can be found in [2, Theorem 4.15]) is the key point to derive our main results in the present work.

**Proposition A.5.** Let $\Psi: \mathbb{R}_+ \times \mathcal{H} \to \mathbb{R} \cup \{+\infty\}$ be a function such that, for all $t \ge 0$, $\Psi(t, \cdot): \mathcal{H} \to \mathbb{R} \cup \{+\infty\}$ is a proper, lower semi-continuous and convex function. Let $F: \mathbb{R}_+ \to \mathcal{H}$ and $u: \mathbb{R}_+ \to \mathcal{H}$ be defined by

$$u(t) := \operatorname{prox}_{\Psi(t, \cdot)}(F(t)),$$

for all $t \ge 0$. If the conditions

(i) $F$ is differentiable at $t=0$;

(ii) $\Psi$ is twice epi-differentiable at $u(0)$ for $F(0) - u(0) \in \partial\Psi(0, \cdot)(u(0))$;

(iii) $D_e^2\Psi(u(0)|F(0) - u(0))$ is a proper function on $\mathcal{H}$;

are satisfied, then $u$ is differentiable at $t=0$ with

$$u'(0) = \operatorname{prox}_{D_e^2\Psi(u(0)|F(0)-u(0))}(F'(0)).$$

# B Reminders on differential geometry

In this appendix, let $\Omega$ be a nonempty bounded connected open subset of $\mathbb{R}^2$ with a Lipschitz boundary $\Gamma := \partial\Omega$ and $n$ be the outward-pointing unit normal vector to $\Gamma$. The next proposition, known as divergence formula, can be found in [6, Theorem 4.4.7 p.104]. The following two propositions are also useful in the present paper and their proofs can be found in [22].

**Proposition B.1 (Divergence formula).** Consider the space

$$H_{\operatorname{div}}(\Omega, \mathbb{R}^{2 \times 2}) := \{w \in L^2(\Omega, \mathbb{R}^{2 \times 2}) \mid \operatorname{div}(w) \in L^2(\Omega, \mathbb{R}^2)\},$$

where $\operatorname{div}(w)$ is the vector whose the $i$-th component is defined by $\operatorname{div}(w)_i := \operatorname{div}(w_i) \in L^2(\Omega, \mathbb{R})$, and where $w_i \in L^2(\Omega, \mathbb{R}^2)$ is the transpose of the $i$-th line of $w$, for all $i \in \{1, 2\}$. If $w \in H_{\operatorname{div}}(\Omega, \mathbb{R}^{2 \times 2})$, then $w$ admits a normal trace, denoted by $w_n \in H^{-1/2}(\Gamma, \mathbb{R}^2)$, satisfying

$$\int_{\Omega} \operatorname{div}(w) \cdot v + \int_{\Omega} w : \nabla v = \langle w_n, v \rangle_{H^{-1/2}(\Gamma, \mathbb{R}^2) \times H^{1/2}(\Gamma, \mathbb{R}^2)}, \quad \forall v \in H^1(\Omega, \mathbb{R}^2).$$

**Proposition B.2.** Assume that $\Gamma$ is of class $C^2$ and let $\theta \in C^1(\mathbb{R}^2, \mathbb{R}^2)$. It holds that

$$\int_{\Gamma} (\theta \cdot \nabla v + v \operatorname{div}_{\tau}(\theta)) = \int_{\Gamma} \theta \cdot n(\partial_n v + H v), \quad \forall v \in W^{2,1}(\Omega, \mathbb{R}),$$

where $\operatorname{div}_{\tau}(\theta) := \operatorname{div}(\theta) - (\nabla\theta n \cdot n) \in L^{\infty}(\Gamma)$ is the tangential divergence of $\theta$, $\partial_n v := \nabla v \cdot n \in L^1(\Gamma, \mathbb{R})$ stands for the normal derivative of $v$, and $H$ stands for the mean curvature of $\Gamma$.

**Proposition B.3.** Assume that $\Gamma$ is of class $C^2$ and let $w \in H^2(\Omega, \mathbb{R}^{2 \times 2})$. It holds that

$$\operatorname{div}(w) = \operatorname{div}_{\tau}(w_{\tau}) + H w_n + (\partial_n w) n \quad \text{a.e. on } \Gamma,$$

where $\operatorname{div}_{\tau}(w_{\tau}) \in L^2(\Gamma, \mathbb{R}^2)$ is the vector whose the $i$-th component is defined by $\operatorname{div}_{\tau}(w_{\tau})_i := \operatorname{div}_{\tau}((w_i)_{\tau}) \in L^2(\Gamma, \mathbb{R})$, where $(w_i)_{\tau} := w_i - (w_i \cdot n)n \in L^2(\Gamma, \mathbb{R}^2)$, and where $\partial_n w \in L^2(\Gamma, \mathbb{R}^{2 \times 2})$ is the matrix whose the $i$-th line is the transpose of the vector $\partial_n w_i := (\nabla w_i)n \in L^2(\Gamma, \mathbb{R}^2)$, for all $i \in \{1, 2\}$. Moreover, it holds that

$$\int_{\Gamma} v \cdot \operatorname{div}_{\tau}(w_{\tau}) = - \int_{\Gamma} w : \nabla_{\tau} v, \quad \forall v \in H^2(\Omega, \mathbb{R}^2),$$

where $\nabla_{\tau} v$ is the matrix whose the $i$-th line is the transpose of the tangential gradient $\nabla_{\tau} v_i := \nabla v_i - (\partial_n v_i)n \in H^{1/2}(\Gamma, \mathbb{R}^2)$, for all $i \in \{1, 2\}$.
References

[1] R. Adams and J. Fournier. Sobolev Spaces, volume 140 of Pure and Applied Mathematics. Elsevier, 2003.

[2] S. Adly and L. Bourdin. Sensitivity analysis of variational inequalities via twice epi-differentiability and proto-differentiability of the proximity operator. *SIAM Journal on Optimization*, 28(2):1699–1725, 2018.

[3] S. Adly, L. Bourdin, and F. Caubet. Sensitivity analysis of a Tresca-type problem leads to Signorini's conditions. *ESAIM: COCV*, 28(29), 2022.

[4] S. Adly, L. Bourdin, F. Caubet, and A. Jacob de Cordemoy. Shape optimization for variational inequalities: the scalar Tresca friction problem. *SIAM Journal on Optimization*, 33(4):2512–2541, 2023.

[5] J. M. Aitchison and M. W. Poole. A numerical algorithm for the solution of Signorini problems. *J. Comput. Appl. Math.*, 94(1):55–67, 1998.

[6] G. Allaire. Analyse numerique et optimisation. Mathématiques et Applications. Éditions de l'École Polytechnique, 2007.

[7] G. Allaire. Conception optimale de structures. Mathématiques et Applications. Springer-Verlag Berlin Heidelberg, 2007.

[8] G. Allaire, F. Jouve, and A. Maury. Shape optimisation with the level set method for contact problems in linearised elasticity. *The SMAI journal of computational mathematics*, 3:249–292, 2017.

[9] P. Beremlijski, J. Haslinger, M. Kočvara, R. Kučera, and J. V. Outrata. Shape optimization in three-dimensional contact problems with Coulomb friction. *SIAM J. Optim.*, 20(1):416–444, 2009.

[10] L. Bourdin, F. Caubet, and A. Jacob de Cordemoy. Sensitivity analysis of a scalar mechanical contact problem with perturbation of the Tresca's friction law. *J Optim Theory Appl*, 192:856–890, 2022.

[11] L. Bourdin, F. Caubet, and A. Jacob de Cordemoy. Sensitivity analysis and optimal control for a friction problem in the linear elastic model. *Appl Math Optim*, 90(29), 2024.

[12] H. Brezis. *Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert*, volume 5 of *North-Holland Mathematics Studies*, No. 5. Notas de Matemática (50). North-Holland Publishing Co., Amsterdam, 1973.

[13] B. Chaudet-Dumas and J. Deteix. Shape derivatives for the penalty formulation of elastic contact problems with tresca friction. *SIAM Journal on Control and Optimization*, 58(6):3237–3261, 2020.

[14] B. Chaudet-Dumas and J. Deteix. Shape derivatives for an augmented lagrangian formulation of elastic contact problems. *ESAIM: Control, Optimisation and Calculus of Variations*, 27(14):23, 2021.

[15] C. N. Do. Generalized second-order derivatives of convex functions in reflexive banach spaces. *Transactions of the American Mathematical Society*, 334(1):281–301, 1992.
[16] G. Duvaut and J.-L. Lions. *Inequalities in mechanics and physics*, volume 219 of *Grundlehren der Mathematischen Wissenschaften*. Springer-Verlag, Berlin-New York, 1976. Translated from the French by C. W. John.

[17] P. Fulmański, A. Laurain, J.-F. Scheid, and J. Sokołowski. A level set method in shape and topology optimization for variational inequalities. *Int. J. Appl. Math. Comput. Sci.*, 17(3):413–430, 2007.

[18] J. Haslinger and A. Klarbring. Shape optimization in unilateral contact problems using generalized reciprocal energy as objective functional. *Nonlinear Anal.*, 21(11):815–834, 1993.

[19] J. Haslinger and R. A. E. Mäkinen. *Introduction to shape optimization*, volume 7 of *Advances in Design and Control*. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2003.

[20] F. Hecht. New development in Freefem++. *J. Numer. Math.*, 20(3-4):251–265, 2012.

[21] C. Heinemann and K. Sturm. Shape optimization for a class of semilinear variational inequalities with applications to damage models. *SIAM Journal on Mathematical Analysis*, 48(5):3579–3617, 2016.

[22] A. Henrot and M. Pierre. *Shape Variation and Optimization : a Geometrical Analysis*. Tracts in Mathematics Vol. 28. European Mathematical Society, 2018.

[23] M. Hintermüller and A. Laurain. Optimal shape design subject to elliptic variational inequalities. *SIAM Journal on Control and Optimization*, 49(3):1015–1047, 2011.

[24] A. Jacob de Cordemoy. Shape optimization for contact problem involving Signorini unilateral conditions. Submitted, 2024.

[25] A. Maury. *Shape optimization for contact and plasticity problems thanks to the level set method*. Thèse, Université Pierre et Marie Curie - Paris VI, 2016.

[26] G. J. Minty. Monotone (nonlinear) operators in Hilbert space. *Duke Mathematical Journal*, 29(3):341–346, 1962.

[27] J. J. Moreau. Proximité et dualité dans un espace hilbertien. *Bulletin de la Société Mathématique de France*, 93:273–299, 1965.

[28] R. T. Rockafellar. On the maximal monotonicity of subdifferential mappings. *Pacific Journal of Mathematics*, 33(1):209 – 216, 1970.

[29] R. T. Rockafellar. Maximal monotone relations and the second derivatives of nonsmooth functions. *Ann. Inst. H. Poincaré Anal. Non Linéaire*, 2(3):167–184, 1985.

[30] R. T. Rockafellar and R. J.-B. Wets. *Variational Analysis*, volume 317 of *Grundlehren der mathematischen Wissenschaften*. Springer-Verlag Berlin Heidelberg, 1998.

[31] J. Salençon. *Handbook of continuum mechanics (general concepts, thermoelasticity)*. Springer-Verlag, Berlin, 2001.

[32] A. Signorini. Sopra alcune questioni di statica dei sistemi continui. *Annali della Scuola Normale Superiore di Pisa - Classe di Scienze*, Ser. 2, 2(2):231–251, 1933.

[33] A. Signorini. Questioni di elasticità non linearizzata e semilinearizzata. *Rend. Mat. Appl., V. Ser.*, 18:95–139, 1959.
[34] J. Sokołowski and J. Zolésio. *Introduction to shape optimization*, volume 16 of *Springer Series in Computational Mathematics*. Springer-Verlag, Berlin, 1992.

[35] J. Sokołowski and J. Zolésio. Shape sensitivity analysis of contact problem with prescribed friction. *Nonlinear Analysis: Theory, Methods & Applications*, 12(12):1399–1411, 1988.