∎

11institutetext: V.A. Kovtunenko 22institutetext: Institute for Mathematics
and Scientific Computing, Karl-Franzens University of Graz, NAWI Graz,
Heinrichstraße 36, 8010 Graz, Austria 22email<EMAIL_ADDRESS>
_and_ Lavrentyev Institute of Hydrodynamics, Siberian Division of the Russian
Academy of Sciences, 630090 Novosibirsk, Russia 33institutetext: K. Kunisch
44institutetext: Institute for Mathematics and Scientific Computing, Karl-
Franzens University of Graz, NAWI Graz, Heinrichstraße 36, 8010 Graz, Austria
44email<EMAIL_ADDRESS>
_and_ Radon Institute, Austrian Academy of Sciences, RICAM Linz,
Altenbergerstraße 69, 4040 Linz, Austria

# Shape Derivative for Penalty-Constrained Nonsmooth–Nonconvex Optimization:
Cohesive Crack Problem

Victor A. Kovtunenko  Karl Kunisch

###### Abstract

A class of non-smooth and non-convex optimization problems with penalty
constraints linked to variational inequalities (VI) is studied with respect to
its shape differentiability. The specific problem stemming from quasi-brittle
fracture describes an elastic body with a Barenblatt cohesive crack under the
inequality condition of non-penetration at the crack faces. Based on the
Lagrange approach and using smooth penalization with the Lavrentiev
regularization, a formula for the shape derivative is derived. The explicit
formula contains both primal and adjoint states and is useful for finding
descent directions for a gradient algorithm to identify an optimal crack shape
from a boundary measurement. Numerical examples of destructive testing are
presented in 2D.

###### Keywords:

Shape optimization Optimal control Variational inequality Penalization
Lagrange method Lavrentiev regularization Free discontinuity problem Non-
penetrating crack Quasi-brittle fracture Destructive physical analysis

###### MSC:

35R37 49J40 49Q10 74RXX

††journal: JOTA

##  1 Introduction

We develop a shape derivative of geometry-dependent least-squares functions
for a class of non-smooth and non-convex optimization problems. The shape
optimization problem is constrained by a penalty equation linked to a
variational inequality (VI). The specific problem describes non-penetrating
cracks with cohesion in the framework of quasi-brittle fracture and
destructive physical analysis (DPA).

Within the general theory for optimal control of VI Bar/84 ; MP/84 , the main
challenge consists in the derivation of optimality conditions. It can be
studied by proper approximation of VI by regularized equations and taking the
limit as the regularization parameter tends to zero. The corresponding methods
for optimal control of obstacle problems can be found in Ber/97 ; HHL/14
using augmented Lagrangians, in e.g. HK/09 ; IK/03  for a Moreau–Yosida
regularization, and in NT/09  based on a Lavrentiev regularization, for the
latter see HKR/16 ; Lav/67 . Furthermore we cite BLR/15 ; CCK/13  for control
of non-smooth and non-convex functionals, KV/07  for boundary control, and
GJKS/18 ; ZMK/21  for control of quasi- and hemi-VI. Shape optimization for
free-interface identification with obstacle-type VI using adjoints was
developed recently by FSW/18 ; LSW/20 . The common difficulty is a lack of
regularity that needs assumptions on a solution in order to take the limit
Sch/22 .

Relying on linearized relations, a crack identification problem was treated
e.g. in BB/13 . We can refer also to ATS/21 ; HIIS/19  for relevant shape
optimization problems in acoustics, to GGK/21  in nonlinear flows subject to
the divergence-free constraint, to KO/20  for over-determined and to HKKP/03
for Bernoulli-type free boundary problems. In the case of non-penetrating
cracks (which are inequality-constrained), the shape differentiability of the
bulk energy was proved in FHLRS/09 ; KK/00  for rectilinear cracks and used
for optimal shape design in KS/18 ; KMZ/08 ; LI/19 ; LSZ/15 . For curvilinear
cracks, adopting the theorem of Correa–Seeger CS/85  on directional
differentiability of Lagrangians the shape derivative was derived in Kov/06 ;
KK/07 , and in Rud/04  using $\Gamma$-convergence.

For the non-penetrating Barenblatt crack that we investigate here, the study
of the objective function and its optimal control with respect to the crack
shape has a number of challenging tasks that we address below. The subsequent
Sections 3–7 follow Tasks (i)–(v), which for convenience are summarized and
explained in the following Sections 2.

##  2 Modeling tasks

Let $t\mapsto\Omega_{t}$ by a parameter (time)-dependent geometry with a crack
$\Gamma_{t}$ along an interface (the breaking line) $\Sigma_{t}$. Denote by
$\nu_{t}$ a normal vector to the surface $\Sigma_{t}$. Motivated by
applications in fracture mechanics (see e.g. BMP/09 ), we consider a total
energy functional $u\mapsto\mathcal{E}:V(\Omega_{t})\mapsto\mathbb{R}$, which
is given in a Hilbert space $V(\Omega_{t})$ by the sum

$\mathcal{E}(u;\Omega_{t})=\mathcal{B}(u;\Omega_{t})+\mathcal{S}([\\![u]\\!];\Sigma_{t}),$
(2.1)

where the bulk term $\mathcal{B}$ is convex, typically, quadratic. The term
$\mathcal{S}$ describes a surface energy according to the Barenblatt idea of a
cohesion zone and depends on the jump $[\\![u]\\!]$ expressing a possible
discontinuity across the interface $\Sigma_{t}$ field $u$. The latter term is
non-convex. The condition of non-penetration (see KK/00 ; KS/97 ) for the
normal opening $\nu_{t}\cdot[\\![u]\\!]\geq 0$ describes the feasible set
$K(\Omega_{t})\subset V(\Omega_{t})$ which is a convex cone. For
differentiable maps $u\mapsto\mathcal{E}$, the first order optimality
condition for the minimization of $\mathcal{E}(u;\Omega_{t})$ over $u\in
K(\Omega_{t})$ results in a VI

$u_{t}\in
K(\Omega_{t}),\quad\langle\partial_{u}\mathcal{E}(u_{t};\Omega_{t}),u-u_{t}\rangle\geq
0\quad\text{for all }u\in K(\Omega_{t}).$ (2.2)

It constitutes a non-convex problem for a solid with a non-penetrating crack
(see Kov/05 ).

For comparison, the classic Griffith model of brittle fracture simplifies
$\mathcal{S}$ to be constant, and a crack $\Gamma_{t}$ to be predefined at the
interface $\Sigma_{t}$. This simplification results in a square-root
singularity of the displacement $u_{t}$ and infinite stress at the crack tip
(front) $\partial\Gamma_{t}$. This is the main disadvantage of the Griffith
model, we refer to CFMT/00  for a discussion. A model, consistent with the
physics of quasi-brittle fracture for non-constant $\mathcal{S}$, was
suggested by Barenblatt Bar/62 . It takes into account the surface cohesion
from the meso-level such that the interface surfaces close in a smooth way,
and thus allow healing of the crack. Indeed, after solving problem (2.2)
according to Barenblatt, the set of points where an opening
$[\\![u_{t}]\\!]\not=0$ occurs, determines the a-priori unknown crack
$\Gamma_{t}$ along the interface $\Sigma_{t}$. This is the complement to the
closed part of the interface where $[\\![u_{t}]\\!]=0$.

The main challenge of the direct problem (2.2) concerns the term $\mathcal{S}$
in (2.1). From an optimization point of view, minimization over feasible $u\in
K(\Omega_{t})$ of $\mathcal{E}$ with a non-smooth surface density
$[\\![u]\\!]\mapsto\mathcal{S}$ (when not a $C^{1}$-function) leads to a hemi-
VI (2.2). The hemi-VI approach was analyzed theoretically and numerically in
HKK/11 ; OG/14  and used in Kov/11 ; KS/06 ; LPSS/13  to describe a quasi-
static crack propagation. A quadratic function $\mathcal{S}$ describing
adhesive cracks was studied in FIR/20 . In the present paper, we study
$C^{2}$-smooth surface energies $\mathcal{S}$ that are small compared to the
bulk term $\mathcal{B}$ in (2.1), see assumption (4.14) below, which is
consistent with meso-level modeling.

Our ultimate aim is to identify the free-interface $\Sigma_{t}$ by a shape
optimization approach as described in GKK/20 . For this task, we introduce the
VI-constrained least-squares misfit from a given measurement $z$ at an
observation boundary $\Gamma^{\rm O}_{t}$:

$\mathcal{J}(u_{t};\Omega_{t})=\frac{1}{2}\int_{\Gamma^{\rm
O}_{t}}|u_{t}-z|^{2}\,dS_{x}+\rho|\Sigma_{t}|\quad\text{such that $u_{t}$
solves \eqref{1.2}},$ (2.3)

where the regularization uses parameter $\rho>0$. This constitutes a
nonsmooth–nonconvex optimization problem.

Our current work focuses on the following tasks.

### Task (i): $C^{2}$-approximation of $\mathcal{E}$.

To provide a shape derivative of $\mathcal{J}$ defined in (2.3) a continuously
differentiable approximation of VI (2.2) is needed. The standard penalization
of non-penetration $\nu_{t}\cdot[\\![u]\\!]\geq 0$ by
$-[\nu_{t}\cdot[\\![u]\\!]]^{-}/\varepsilon$ has only $C^{0}$-regularity. Here
the regularization parameter $\varepsilon>0$ is small, and $u=[u]^{+}-[u]^{-}$
implies the decomposition into positive $[u]^{+}=\max(0,u)$ and negative
$[u]^{-}=-\min(0,u)$ parts. Therefore, we suggest a $C^{1}$-penalization by
the normal compliance $\beta_{\epsilon}(\nu_{t}\cdot[\\![u]\\!])$ based on the
Lavrentiev regularization (see Theorem 4.1). This results in a
$C^{2}$-approximation of $\mathcal{E}$ for the $\varepsilon$-approximation of
(2.1)–(2.3) by

$\mathcal{J}(u^{\varepsilon}_{t};\Omega_{t})=\frac{1}{2}\int_{\Gamma^{\rm
O}_{t}}|u^{\varepsilon}_{t}-z|^{2}\,dS_{x}+\rho|\Sigma_{t}|,\quad\text{where
}\partial^{\varepsilon}_{u}\mathcal{E}(u^{\varepsilon}_{t};\Omega_{t})=0,$
(2.4)

and the penalty equation involves the operator
$\partial^{\varepsilon}_{u}\mathcal{E}$ introduced as follows

$\langle\partial^{\varepsilon}_{u}\mathcal{E}(u;\Omega_{t}),v\rangle:=\langle\partial_{u}\mathcal{E}(u;\Omega_{t}),v\rangle+\int_{\Sigma_{t}}\beta_{\epsilon}(\nu_{t}\cdot[\\![u]\\!])\,(\nu_{t}\cdot[\\![v]\\!])\,dS_{x}.$
(2.5)

### Task (ii): adjoint-based optimality conditions.

Applying to the penalty-constrained least-square misfit (2.4) a Lagrange
multiplier approach (see IK/08 ), we can define an $\varepsilon$-dependent
Lagrangian
$(u,v)\mapsto\mathcal{L}^{\varepsilon}:V(\Omega_{t})^{2}\mapsto\mathbb{R}$ as

$\mathcal{L}^{\varepsilon}(u,v;\Omega_{t})=\mathcal{J}(u;\Omega_{t})-\langle\partial^{\varepsilon}_{u}\mathcal{E}(u;\Omega_{t}),v\rangle.$
(2.6)

The primal (inf-sup) problem: for fixed $v^{\varepsilon}_{t}\in
V(\Omega_{t})$, find $u^{\varepsilon}_{t}\in V(\Omega_{t})$ such that

$\mathcal{L}^{\varepsilon}(u^{\varepsilon}_{t},v;\Omega_{t})\leq\mathcal{L}^{\varepsilon}(u^{\varepsilon}_{t},v^{\varepsilon}_{t};\Omega_{t})\quad\text{for
all }v\in V(\Omega_{t}).$ (2.7)

Since $\mathcal{L}^{\varepsilon}$ is affine in $v$, the first order optimality
condition is given by

$u^{\varepsilon}_{t}\in
V(\Omega_{t}),\quad\langle\partial^{\varepsilon}_{u}\mathcal{E}(u^{\varepsilon}_{t};\Omega_{t}),u\rangle=0\quad\text{for
all }u\in V(\Omega_{t}).$

The dual (sup-inf) problem (see (ET/76, , Chapter 6)) reads: for fixed
$u^{\varepsilon}_{t}\in V(\Omega_{t})$, find $v^{\varepsilon}_{t}\in
V(\Omega_{t})$ such that

$\mathcal{L}^{\varepsilon}(u^{\varepsilon}_{t},v^{\varepsilon}_{t};\Omega_{t})\leq\mathcal{L}^{\varepsilon}(u,v^{\varepsilon}_{t};\Omega_{t})\quad\text{for
all }u\in V(\Omega_{t}).$

Note that $\mathcal{L}^{\varepsilon}$ with respect to $u$ is not a linear
continuous functional on the dual space $V(\Omega_{t})^{\star}$.

The corresponding nonlinear optimization theory was developed in e.g. IK/08 ;
MZ/79 ; ZK/79  as follows. If the variation
$\partial_{u}(\partial^{\varepsilon}_{u}\mathcal{E})\in\mathscr{L}(V(\Omega_{t}),V(\Omega_{t})^{\star})$
with respect to $u$ in (2.5) exists, and the associated adjoint operator
$[\partial_{u}(\partial^{\varepsilon}_{u}\mathcal{E})(u^{\varepsilon}_{t};\Omega_{t})]^{\star}\in\mathscr{L}(V(\Omega_{t}),V(\Omega_{t})^{\star})$
satisfying

$\langle[\partial_{u}(\partial^{\varepsilon}_{u}\mathcal{E})(u^{\varepsilon}_{t};\Omega_{t})]^{\star}v,u\rangle=\langle\partial_{u}(\partial^{\varepsilon}_{u}\mathcal{E})(u^{\varepsilon}_{t};\Omega_{t})u,v\rangle\quad\text{for
all }u,v\in V(\Omega_{t})$

is surjective with respect to $u^{\varepsilon}_{t}$, then the optimality
condition is given by

$v^{\varepsilon}_{t}\in
V(\Omega_{t}),\quad\langle[\partial_{u}(\partial^{\varepsilon}_{u}\mathcal{E})(u^{\varepsilon}_{t};\Omega_{t})]^{\star}v^{\varepsilon}_{t},v\rangle=0\quad\text{for
all }v\in V(\Omega_{t}).$ (2.8)

For the abstract theory associated to adjoint operators we cite EL/13 ; Ibr/06
; MAS/96 . To justify (2.8), we shall linearize $\partial_{u}\mathcal{E}$
around the primal solution $u^{\varepsilon}_{t}$ to (2.7) (see Theorem 4.2)
and suggest a suitable linearized functional
$(u,v)\mapsto\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u,v):V(\Omega_{t})^{2}\mapsto\mathbb{R}$
such that

$\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u^{\varepsilon}_{t},v;\Omega_{t})=\mathcal{L}^{\varepsilon}(u^{\varepsilon}_{t},v;\Omega_{t})\quad\text{for
}v\in V(\Omega_{t}).$ (2.9)

### Task (iii): shape derivative.

Our purpose is to calculate a shape derivative of the mapping
$t\mapsto\mathcal{J}(u^{\varepsilon}_{t};\Omega_{t})$ that is expressed by the
one-sided limit (see DZ/11 ; SZ/92 ):

$\partial_{t}\mathcal{J}(u^{\varepsilon}_{t};\Omega_{t})=\lim_{s\to
0^{+}}\frac{1}{s}\bigl{(}\mathcal{J}(u^{\varepsilon}_{t+s};\Omega_{t+s})-\mathcal{J}(u^{\varepsilon}_{t};\Omega_{t})\bigr{)}.$
(2.10)

If a saddle-point $(u^{\varepsilon}_{t},v^{\varepsilon}_{t})\in
V(\Omega_{t})^{2}$ based on (2.7) and (2.9) exists, then the optimal value
misfit function defined in (2.4) is evidently equal to the optimal value
Lagrange function

$\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u^{\varepsilon}_{t},v^{\varepsilon}_{t};\Omega_{t})=\mathcal{L}^{\varepsilon}(u^{\varepsilon}_{t},v^{\varepsilon}_{t};\Omega_{t})=\mathcal{J}(u^{\varepsilon}_{t};\Omega_{t})-\langle\partial^{\varepsilon}_{u}\mathcal{E}(u^{\varepsilon}_{t};\Omega_{t}),v^{\varepsilon}_{t}\rangle\quad\text{subject
to }\\\
\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u^{\varepsilon}_{t},v;\Omega_{t})\leq\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u^{\varepsilon}_{t},v^{\varepsilon}_{t};\Omega_{t})\leq\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u,v^{\varepsilon}_{t};\Omega_{t})\\\
\text{for all }(u,v)\in V(\Omega_{t})^{2}.$ (2.11)

Henceforth, we have the following identity for the shape derivative according
to (2.10):

$\partial_{t}\mathcal{J}(u^{\varepsilon}_{t};\Omega_{t})=\partial_{t}\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u^{\varepsilon}_{t},v^{\varepsilon}_{t};\Omega_{t})\\\
=\lim_{s\to
0^{+}}\frac{1}{s}\bigl{(}\tilde{\mathcal{L}}^{\varepsilon}(s,u^{\varepsilon}_{t},u^{\varepsilon}_{t+s},v^{\varepsilon}_{t+s};\Omega_{t})-\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u^{\varepsilon}_{t},v^{\varepsilon}_{t};\Omega_{t})\bigr{)}.$
(2.12)

In order to construct a proper $\tilde{\mathcal{L}}^{\varepsilon}$, using a
diffeomorphic coordinate transformation $y=\phi_{s}(x)$ such that
$\phi_{s}:\Omega_{t}\mapsto\Omega_{t+s}$ (see (SZ/92, , Chapter 2)), the
bijection $V(\Omega_{t+s})\mapsto V(\Omega_{t})$, $u\mapsto u\circ\phi_{s}$
provides the perturbed Lagrangian as

$\widetilde{\mathcal{L}}^{\varepsilon}(s,u\circ\phi_{s},u\circ\phi_{s},v\circ\phi_{s};\Omega_{t})=\mathcal{L}^{\varepsilon}(u,v;\Omega_{t+s})\quad\text{for
}(u,v)\in V(\Omega_{t+s})^{2}.$ (2.13)

Then the results of Delfour–Zolesio DZ/11  on shape differentiabiliy can be
applied to justify the limit in (2.12), see respective Theorem 5.1 and its
Corollary 1.

### Task (iv): limit as $\varepsilon\to 0^{+}$.

Taking the limit as $\varepsilon\to 0^{+}$ in relations (2.11) we shall prove
the optimality conditions (see Theorem 6.1 and its Corollary 3). However, we
cannot pass to the limit in (2.12) due to the presence of the unbounded term
$\beta_{\varepsilon}^{\prime}$. We conjecture that the limit problem (2.2) is
not differentiable. This agrees with the assertion that VIs are not Fréchet
differentiable with respect to shape (see LSZ/15 ). Therefore, in the
numerical treatment we rely on the approximation (2.12) with small
$\varepsilon>0$ for the shape derivative $\partial_{t}\mathcal{J}$.

### Task (v): shape optimization.

Commonly adopted in shape optimization, the gradient method needs a descent
direction minimizing the objective map $\Omega_{t}\mapsto\mathcal{J}$ such
that $\partial_{t}\mathcal{J}<0$. This can be attained by a proper choice of
the transformation $\phi_{s}$ entering implicitly in formula (2.12) (see
Corollary 2). Realizing the optimization algorithm for crack shape
identification, from our numerical tests we report the following feature.
Those parts of the crack faces which are in contact (where the non-penetration
constraint is active) are hidden from identification. To identify a crack
needs its faces to be open (that is, VI turns into unconstrained equation) in
accordance with the concept of destructive physical analysis (DPA).

##  3 Cohesive crack problem

We start with a detailed description of the geometry. Let
$\Omega\subset\mathbb{R}^{d}$, $d=2,3$, be a fixed hold-all domain with
Lipschitz boundary $\partial\Omega$. For the time-parameter
$t\in(t_{0},t_{1})$, $t_{0}<t_{1}$, we consider a parameter-dependent geometry
$\Omega_{t}=(\Gamma^{\rm D}_{t},\Gamma^{\rm N}_{t},\Gamma^{\rm
O}_{t},\Sigma_{t})$ defined as follows. For brevity we use a single notation
$\Omega_{t}$ for the collection of geometric objects describing a broken
domain $\Omega\setminus\Sigma_{t}$ by means of the Dirichlet, Neumann,
observation boundaries, and the breaking line, respectively.

The outer boundary is split into two variable parts such that
$\partial\Omega=\overline{\Gamma^{\rm D}_{t}}\cup\overline{\Gamma^{\rm
N}_{t}}$ and $\Gamma^{\rm D}_{t}\cap\Gamma^{\rm N}_{t}=\emptyset$ with normal
vector $n_{t}$ outward to $\Omega$. The observation boundary is $\Gamma^{\rm
O}_{t}\subset\Gamma^{\rm N}_{t}$. The domain is split into two variable sub-
domains $\Omega^{\pm}_{t}$ with Lipschitz boundaries
$\partial\Omega^{\pm}_{t}$ and outward normal vectors $n^{\pm}_{t}$ such that
$n^{\pm}_{t}=n_{t}$ at $\partial\Omega$. The conditions $\Gamma^{\rm
D}_{t}\cap\partial\Omega^{+}_{t}\not=\emptyset$ and $\Gamma^{\rm
D}_{t}\cap\partial\Omega^{-}_{t}\not=\emptyset$ are needed to guarantee the
Korn–Poincare inequality. These two domains are separated by a breaking
manifold (the free-interface)
$\Sigma_{t}=\partial\Omega^{+}_{t}\cap\partial\Omega^{-}_{t}$ with normal
direction $\nu_{t}=n^{-}_{t}=-n^{+}_{t}$ such that
$\Omega=\Omega^{+}_{t}\cup\Omega^{-}_{t}\cup\Sigma_{t}$. An example geometry
of $\Omega_{t}$ is sketched in 2D in Figure 1.

$\;\Gamma^{\rm D}_{t}\;$$\Omega^{+}_{t}$$\Omega^{-}_{t}$$n_{t}$$\Gamma^{\rm
D}_{t}\;$$n_{t}$$\Gamma^{\rm N}_{t}\;$$n_{t}$$\Gamma^{\rm
N}_{t}\;$$n_{t}$$\Sigma_{t}$$\nu_{t}$$\tau_{t}$ Figure 1: An example
configuration of variable geometry $\Omega_{t}$ in 2D.

We assume that these geometric properties are preserved for all
$t\in(t_{0},t_{1})$ under suitable shape perturbations, which we specify below
in Section 5.

For fixed $t$, we consider a linear elastic body that occupies the
disconnected domain
$\Omega\setminus\Sigma_{t}=\Omega^{+}_{t}\cup\Omega^{-}_{t}$. By this,
$d$-dimensional vectors of displacement $u(x)$ at points
$x\in\Omega\setminus\Sigma_{t}$ admit discontinuity across $\Sigma_{t}$
resulting in the jump
$[\\![u]\\!]=u|_{\Sigma_{t}\cap\partial\Omega^{+}_{t}}-u|_{\Sigma_{t}\cap\partial\Omega^{-}_{t}}$.
For further use we employ an orthogonal decomposition of admissible
$[\\![u]\\!]$ into the normal component with factor $\nu_{t}\cdot[\\![u]\\!]$
and the tangential vector $[\\![u]\\!]_{\tau_{t}}$ at the interface such that

$[\\![u]\\!]=(\nu_{t}\cdot[\\![u]\\!])\nu_{t}+[\\![u]\\!]_{\tau_{t}},\quad\nu_{t}\cdot[\\![u]\\!]\geq
0\quad\text{on }\Sigma_{t}.$ (3.1)

The latter inequality in (3.1) describes the non-penetration, see KK/00 .

The essential issue of modeling is to introduce a density at $\Sigma_{t}$ for
the surface energy $\mathcal{S}$ in (2.1) that is consistent with physics.
Based on the decomposition (3.1), we set

$\mathcal{S}([\\![u]\\!];\Sigma_{t})=\int_{\Sigma_{t}}\bigl{\\{}\alpha_{\rm
f}([\\![u]\\!]_{\tau_{t}})+\alpha_{\rm
c}(\nu_{t}\cdot[\\![u]\\!])\bigr{\\}}\,dS_{x}.$ (3.2)

The former, shear-induced term in (3.2), is associated with friction between
the crack surfaces. Let the mapping $\xi\mapsto\alpha_{\rm
f}(\xi):\mathbb{R}^{d}\mapsto\mathbb{R}$, and its first and second derivatives
be uniformly continuous functions, satisfying for constants $K_{\rm f},K_{\rm
f1},K_{\rm f2}\geq 0$ and all $\xi$,

$-K_{\rm f}|\xi|\leq\alpha_{\rm f}(\xi),\quad|\nabla\alpha_{\rm f}(\xi)|\leq
K_{\rm f1},\quad|\nabla^{2}\alpha_{\rm f}(\xi)|\leq K_{\rm f2}.$ (3.3)

For example, we have in mind a standard regularization of the Coulomb law (see
e.g. (SM/91, , Section 4.3.3)) with the positive, convex function

$\alpha_{\rm f}(\xi)=F_{\bf b}\sqrt{\delta^{2}+|\xi|^{2}},$ (3.4)

where $\delta>0$ is small, and $F_{\bf b}>0$ is the friction bound. In this
case, $K_{\rm f}=0$, $K_{\rm f1}=F_{\bf b}$, and $K_{\rm f2}=F_{\bf
b}/\delta$. For convenience, the function $\alpha_{\rm f}(s)$ in one variable
$s\in\mathbb{R}$ together with its first two derivatives are depicted in
Figure 2.

$\alpha_{\rm f}(s)$$s$$0$$\delta$

$\alpha^{\prime}_{\rm f}(s)$$s$$0$$-K_{\rm f1}$$K_{\rm f1}$

$\alpha^{\prime\prime}_{\rm f}(s)$$s$$0$$K_{\rm f2}$

Figure 2: Example graphics of $\alpha_{\rm f},\alpha^{\prime}_{\rm
f},\alpha^{\prime\prime}_{\rm f}$ in 1d.

The latter term in (3.2) associates cohesion between the crack surfaces. Let
$s\mapsto\alpha_{\rm c}(s):\mathbb{R}\mapsto\mathbb{R}$ and its second
derivative be uniformly continuous functions, and let there exist constant
$K_{\rm c},K_{\rm c1},K_{\rm c2}\geq 0$ such that

$-K_{\rm c}|s|\leq\alpha_{\rm c}(s),\quad|\alpha_{\rm c}^{\prime}(s)|\leq
K_{\rm c1},\quad|\alpha_{\rm c}^{\prime\prime}(s)|\leq K_{\rm c2}.$ (3.5)

From the physics literature (e.g. Kit/08 ) we suggest the following generic
function

$\alpha_{\rm c}(s)=K_{\rm c}\,\frac{s}{\kappa+|s|^{m}},$ (3.6)

where $K_{\rm c}>0$ is related to the fracture toughness, and $\kappa>0$,
$m\geq 1$ are parameters. In this case, $K_{\rm c1}$ and $K_{\rm c2}$ are
proportional to $K_{\rm c}$. The example of $\alpha_{\rm c},\alpha_{\rm
c}^{\prime},\alpha_{\rm c}^{\prime\prime}$ for $m=4$ is depicted in Figure 3.

$\alpha_{\rm c}(s)$$s$$0$

$\alpha_{\rm c}^{\prime}(s)$$s$$0$$-K_{\rm c1}$$K_{\rm c1}$

$\alpha_{\rm c}^{\prime\prime}(s)$$s$$0$$-K_{\rm c2}$$K_{\rm c2}$

Figure 3: Example graphics of $\alpha_{\rm c},\alpha_{\rm
c}^{\prime},\alpha_{\rm c}^{\prime\prime}$ as $\kappa=1$ and $m=4$.

In particular, the left plot in Figure 3 depicts the typical softening
phenomenon for growing $s$. It is worth noting that the left branch of
$\alpha_{\rm c}(\nu_{t}\cdot[\\![u]\\!])$ for $\nu_{t}\cdot[\\![u]\\!]<0$
implies a normal compliance and it is avoided when the non-penetration
$\nu_{t}\cdot[\\![u]\\!]\geq 0$ in (3.1) holds.

The symmetric $d$-by-$d$ tensors of linearized strain $\epsilon$ and the
Cauchy stress $\sigma$ are given by

$\epsilon(u)=\frac{1}{2}(\nabla u+\nabla
u^{\top}),\quad\sigma(u)=C\epsilon(u),$ (3.7)

where $(\nabla u)=(\partial u_{i}/\partial x_{j})$ for $i,j=1,\ldots,d$, the
transposition $(\,\cdot\,)^{\top}$ swaps columns for rows. A symmetric fourth
order tensor of elastic coefficients $C(x)\in W^{1,\infty}(\Omega)^{d\times
d\times d\times d}$, such that $C_{ijkl}=C_{jikl}=C_{klij}$ for
$i,j,k,l=1,\ldots,d$, is positive definite and fulfills the Korn–Poincare
inequality: there exists $K_{\rm KP}>0$ such that

$\int_{\Omega\setminus\Sigma_{t}}\sigma(u)\cdot\epsilon(u)\,dx\geq K_{\rm
KP}\|u\|^{2}_{H^{1}(\Omega\setminus\Sigma_{t})^{d}}\quad\text{ for }u\in
V(\Omega_{t}).$ (3.8)

over the Sobolev space

$V(\Omega_{t})=\\{u\in H^{1}(\Omega^{+}_{t})^{d}\cap
H^{1}(\Omega^{-}_{t})^{d}|\quad u=0\text{ on }\Gamma^{\rm D}_{t}\\}.$ (3.9)

For a boundary traction vector $g\in H^{1}(\partial\Omega)^{d}$, we consider
the following bulk energy

$\mathcal{B}(u;\Omega_{t})=\frac{1}{2}\int_{\Omega\setminus\Sigma_{t}}\sigma(u)\cdot\epsilon(u)\,dx-\int_{\Gamma^{\rm
N}_{t}}g\cdot u\,dS_{x}.$ (3.10)

The feasible set corresponding to the non-penetration condition in (3.1) reads

$K(\Omega_{t})=\\{u\in V(\Omega_{t})|\quad\nu_{t}\cdot[\\![u]\\!]\geq 0\text{
on }\Sigma_{t}\\},$ (3.11)

which is a convex, closed cone.

######  Theorem 3.1 (Well-posedness of cohesive crack problem)

There exists a solution to the non-convex, constrained minimization problem:
find $u_{t}\in K(\Omega_{t})$ such that

$\mathcal{E}(u_{t};\Omega_{t})=\min_{u\in
K(\Omega_{t})}\mathcal{E}(u;\Omega_{t}),$ (3.12)

where the total energy $\mathcal{E}$ according to (3.2) and (3.10) is given by

$\mathcal{E}(u;\Omega_{t})=\frac{1}{2}\int_{\Omega\setminus\Sigma_{t}}\sigma(u)\cdot\epsilon(u)\,dx-\int_{\Gamma^{\rm
N}_{t}}g\cdot u\,dS_{x}+\int_{\Sigma_{t}}\bigl{\\{}\alpha_{\rm
f}([\\![u]\\!]_{\tau_{t}})+\alpha_{\rm
c}(\nu_{t}\cdot[\\![u]\\!])\bigr{\\}}\,dS_{x}.$ (3.13)

The solution satisfies the first-order optimality condition (2.2) in the form
of VI:

$\int_{\Omega\setminus\Sigma_{t}}\sigma(u_{t})\cdot\epsilon(u-u_{t})\,dx+\int_{\Sigma_{t}}\bigl{\\{}\nabla\alpha_{\rm
f}([\\![u_{t}]\\!]_{\tau_{t}})\cdot[\\![u-u_{t}]\\!]_{\tau_{t}}\\\
+\alpha_{\rm
c}^{\prime}(\nu_{t}\cdot[\\![u_{t}]\\!])\,(\nu_{t}\cdot[\\![u-u_{t}]\\!])\bigr{\\}}\,dS_{x}\geq\int_{\Gamma^{\rm
N}_{t}}g\cdot(u-u_{t})\,dS_{x}$ (3.14)

for all test functions $u\in K(\Omega_{t})$. For smooth solutions the boundary
value relations hold:

$\displaystyle{\rm div}\,\sigma(u_{t})=0$ $\displaystyle\text{ in
}\Omega\setminus\Sigma_{t},$ $\displaystyle u_{t}=0\text{ on }\Gamma^{\rm
D}_{t},\quad\sigma(u_{t})n=g$ $\displaystyle\text{ on }\Gamma^{\rm N}_{t},$
$\displaystyle[\\![\sigma(u_{t})\nu_{t}]\\!]=0,\quad(\sigma(u_{t})\nu_{t})_{\tau_{t}}=\nabla\alpha_{\rm
f}([\\![u_{t}]\\!]_{\tau_{t}}),$ $\displaystyle\nu_{t}\cdot[\\![u_{t}]\\!]\geq
0,\quad\nu_{t}\cdot(\sigma(u_{t})\nu_{t})\leq\alpha_{\rm
c}^{\prime}(\nu_{t}\cdot[\\![u_{t}]\\!]),$
$\displaystyle(\nu_{t}\cdot[\\![u_{t}]\\!])\,\bigl{\\{}\nu_{t}\cdot(\sigma(u_{t})\nu_{t})-\alpha_{\rm
c}^{\prime}(\nu_{t}\cdot[\\![u_{t}]\\!])\bigr{\\}}=0$ $\displaystyle\text{ on
}\Sigma_{t},$ (3.15)

for the decomposition of vector
$\sigma(u_{t})\nu_{t}=\bigl{(}\nu_{t}\cdot(\sigma(u_{t})\nu_{t})\bigr{)}\nu_{t}+(\sigma(u_{t})\nu_{t})_{\tau_{t}}$
according to (3.1). The last two lines in (3.1) are the complementarity
conditions. If both $\alpha_{\rm f}$ and $\alpha_{\rm c}$ were convex (that is
not $\alpha_{\rm c}$ in (3.6)), then the solution $u_{t}$ to (3.12) and (3.14)
would be unique.

###### Proof

On the right-hand side of (3.13), the first, quadratic in $u$ integral term
over $\Omega\setminus\Sigma_{t}$, is strongly positive by the Korn–Poincare
inequality (3.8). Using the Cauchy–Schwarz inequality, the other boundary
integral terms over $\Sigma_{t}$ and $\Gamma^{\rm N}_{t}$ are bounded from
below by a sub-linear in $u$ function

$\int_{\Sigma_{t}}\bigl{\\{}\alpha_{\rm f}([\\![u]\\!]_{\tau_{t}})+\alpha_{\rm
c}(\nu_{t}\cdot[\\![u]\\!])\bigr{\\}}\,dS_{x}-\int_{\Gamma^{\rm N}_{t}}g\cdot
u\,dS_{x}\\\ \geq-\bigl{(}K_{\rm
f}\|[\\![u]\\!]_{\tau_{t}}\|_{L^{2}(\Sigma_{t})^{d}}+K_{\rm
c}\|\nu_{t}\cdot[\\![u]\\!]\|_{L^{2}(\Sigma_{t})}\bigr{)}\sqrt{|\Sigma_{t}|}-\|g\|_{L^{2}(\Gamma^{\rm
N}_{t})^{d}}\|u\|_{L^{2}(\Gamma^{\rm N}_{t})^{d}}$ (3.16)

by virtue of the properties for $\alpha_{\rm f}$, $\alpha_{\rm c}$ in (3.3),
(3.5). Therefore, estimating the jump by
$\|[\\![u]\\!]\|^{2}_{L^{2}(\Sigma_{t})}\leq
2\|u\|^{2}_{L^{2}(\Sigma_{t}\cap\partial\Omega^{+}_{t})}+2\|u\|^{2}_{L^{2}(\Sigma_{t}\cap\partial\Omega^{-}_{t})}$
and applying the trace inequality we have

$\|u\|_{L^{2}(\partial\Omega^{\pm}_{t})^{d}}\leq\|u\|_{H^{1/2}(\partial\Omega^{\pm}_{t})^{d}}\leq
K_{\rm tr}\|u\|_{H^{1}(\Omega^{\pm}_{t})^{d}},\quad u\in
H^{1}(\Omega^{\pm}_{t})^{d}.$ (3.17)

Then we get that $\mathcal{E}$ is radially unbounded, and thus coercive. The
functions $\alpha_{\rm f}$ and $\alpha_{\rm c}$ are uniformly continuous,
hence preserving $L^{2}$-convergence (see BJ/61 ). Using the compactness of
the embedding of the traces of $u$ at
$\Sigma_{t}\cap\partial\Omega^{\pm}_{t}$, from $H^{1}(\Omega^{\pm}_{t})$ into
$L^{2}(\partial\Omega^{\pm}_{t})$, it follows that the mapping
$u\mapsto\mathcal{E}(u)$ from $V(\Omega_{t})\mapsto\mathbb{R}$ is weakly lower
semi-continuous.

Let $\\{u^{n}\\}$, $n\in\mathbb{N}$, be an infimal sequence in
$K(\Omega_{t})$. The coercivity of $\mathcal{E}$ implies the boundedness of
$\\{u^{n}\\}$ in $V(\Omega_{t})$. Then, on a subsequence $\\{u^{n_{k}}\\}$,
there exists an accumulation point $u_{t}$ such that $u^{n_{k}}\rightharpoonup
u_{t}$ weakly in $H^{1}(\Omega\setminus\Sigma_{t})^{d}$ as $n_{k}\to\infty$.
By weak closedness of $K(\Omega_{t})$ we have $u_{t}\in K(\Omega_{t})$. Taking
the limit inferior of $\mathcal{E}(u^{n_{k}})$, the weak lower semi-continuity
of $\mathcal{E}$ implies that $u_{t}$ attains the minimum in (3.12). Applying
standard variational arguments implies the optimality condition (3.14) and
(3.1), see details in (KK/00, , Section 1.4). Moreover, if $\alpha_{\rm f}$,
$\alpha_{\rm c}$ were convex, then the integral over $\Sigma_{t}$ in (3.13) is
monotone. This would lead to uniqueness of $u_{t}$ as solution to (3.14),
which is then necessarily the unique solution for (3.12). ∎

Next we approximate the VI (3.14) by a penalty method. By itself penalization
is a self-contained physical model allowing compliance, see And/99  for the
discussion.

##  4 Lavrentiev based regularization and saddle-point problem

Let $\varepsilon_{0}>0$. For $\varepsilon\in(0,\varepsilon_{0})$, the standard
penalization of the inequality constraint $s\geq 0$ by $-[s]^{-}/\varepsilon$
has only a generalized derivative $\mathcal{H}(-s)/\varepsilon$, where
$\mathcal{H}$ is the Heaviside step function such that $\mathcal{H}(s)=1$ for
$s>0$, otherwise $\mathcal{H}(s)=0$ for $s\leq 0$. We suggest a Lavrentiev
based $C^{1}$-regularization by the normal compliance $\beta_{\varepsilon}$ as
follows. Let the function
$s\mapsto\beta_{\varepsilon}(s):\mathbb{R}\mapsto\mathbb{R}$ be concave and
differentiable, with $\beta$ and $\beta^{\prime}$ uniformly continuous, and
let there exist $K_{\beta},K_{\beta 1}\geq 0$ such that

$\bigl{|}\beta_{\epsilon}(s)+\frac{[s]^{-}}{\varepsilon}\bigr{|}\leq
K_{\beta},\quad 0\leq\beta^{\prime}_{\epsilon}(s)\leq\frac{K_{\beta
1}}{\varepsilon}.$ (4.1)

We assume that the following conditions hold, which describe relaxed
complementarity and compliance, respectively:

$\beta_{\epsilon}(s)[s]^{+}\geq-\varepsilon
K_{\beta},\quad\beta_{\epsilon}(s)[s]^{-}\leq-\frac{([s]^{-})^{2}}{\varepsilon}+\varepsilon
K_{\beta}.$ (4.2)

$s$$-\varepsilon$$\varepsilon$$0$$-1$$\beta_{\epsilon}(s)$

$s$$-\varepsilon$$\varepsilon$$0$$\frac{1}{\varepsilon}$$\beta^{\prime}_{\epsilon}(s)$

$s$$-\varepsilon$$\varepsilon$$0$$\beta^{\prime\prime}_{\epsilon}(s)$

Figure 4: Example graphics of
$\beta_{\varepsilon},\beta^{\prime}_{\varepsilon},\beta^{\prime\prime}_{\varepsilon}$
for fixed $\varepsilon$.

For example, we construct the following mollification of minimum function

$\beta_{\epsilon}(s)=\begin{cases}s/\varepsilon&\text{ for }s<-\varepsilon\\\
-\exp\bigl{(}2(s+\varepsilon)/(s-\varepsilon)\bigr{)}&\text{ for
}-\varepsilon\leq s<\varepsilon\\\ 0&\text{ for }s\geq\varepsilon\end{cases}$
(4.3)

which is depicted in Figure 4 together with its two derivatives.

###### Lemma 1

For $\beta_{\varepsilon}$ from (4.3), the properties (4.1) and (4.2) hold true
with $K_{\beta}=K_{\beta 1}=1$. Moreover,
$\beta^{\prime\prime}_{\varepsilon}\leq 0$ implies that
$\beta^{\prime}_{\varepsilon}\geq 0$ decreases monotonically, and
$\beta_{\varepsilon}\leq 0$ is concave and increases monotonically.

###### Proof

The properties (4.1) can be easily checked. To verify the first inequality in
(4.2), from (4.3) we deduce that $\beta_{\epsilon}(s)[s]^{+}=0$ for
$s\geq\varepsilon$. Here we use the complementary condition $[s]^{-}[s]^{+}=0$
and $[s]^{+}=0$ for $s<0$. We further have
$\beta_{\epsilon}(s)\geq-[s]^{-}/\varepsilon-K_{\beta}$ according to the first
estimate in (4.1). Henceforth, after multiplication with
$[s]^{+}\in[0,\varepsilon)$, the lower bound
$\beta_{\epsilon}(s)[s]^{+}\geq-\varepsilon K_{\beta}$ holds for $0\leq
s<\varepsilon$.

Similarly, $\beta_{\epsilon}(s)[s]^{-}=-([s]^{-})^{2}/\varepsilon$ for
$s<-\varepsilon$ in (4.3), and $\beta_{\epsilon}(s)[s]^{-}=0$ due to
$[s]^{-}=0$ for $s\geq 0$. The first estimate in (4.1), that is
$\beta_{\epsilon}(s)\leq-[s]^{-}/\varepsilon+K_{\beta}$, after multiplication
with $[s]^{-}\in(0,\varepsilon]$ leads to the upper bound
$\beta_{\epsilon}(s)[s]^{-}\leq-([s]^{-})^{2}/\varepsilon+\varepsilon
K_{\beta}$ for $-\varepsilon\leq s<0$. This proves the second inequality in
(4.2). ∎

Using Lemma 1 we obtain the existence result for the penalized cohesive crack
problem.

######  Theorem 4.1 (Well-posedness of $\varepsilon$-regularized cohesive
crack problem)

There exists a solution to the penalty problem: find $u^{\varepsilon}_{t}\in
V(\Omega_{t})$ such that

$\int_{\Omega\setminus\Sigma_{t}}\sigma(u^{\varepsilon}_{t})\cdot\epsilon(u)\,dx+\int_{\Sigma_{t}}\bigl{\\{}\nabla\alpha_{\rm
f}([\\![u^{\varepsilon}_{t}]\\!]_{\tau_{t}})\cdot[\\![u]\\!]_{\tau_{t}}\\\
+[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!])\,(\nu_{t}\cdot[\\![u]\\!])\bigr{\\}}\,dS_{x}=\int_{\Gamma^{\rm
N}_{t}}g\cdot u\,dS_{x}$ (4.4)

for all test functions $u\in V(\Omega_{t})$. For smooth solutions the boundary
value relations hold:

$\displaystyle{\rm div}\,\sigma(u^{\varepsilon}_{t})=0$ $\displaystyle\text{
in }\Omega\setminus\Sigma_{t},$ $\displaystyle u^{\varepsilon}_{t}=0\text{ on
}\Gamma^{\rm D}_{t},\quad\sigma(u^{\varepsilon}_{t})n=g$ $\displaystyle\text{
on }\Gamma^{\rm N}_{t},$
$\displaystyle[\\![\sigma(u^{\varepsilon}_{t})\nu_{t}]\\!]=0,\;(\sigma(u^{\varepsilon}_{t})\nu_{t})_{\tau_{t}}=\nabla\alpha_{\rm
f}([\\![u^{\varepsilon}_{t}]\\!]_{\tau_{t}}),\;\nu_{t}\cdot(\sigma(u^{\varepsilon}_{t})\nu_{t})=[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!])$
$\displaystyle\text{ on }\Sigma_{t}.$ (4.5)

If both $\nabla\alpha_{\rm f}$ and $\alpha_{\rm c}^{\prime}$ were monotone,
then the solution $u^{\varepsilon}_{t}$ to (4.4) would be unique.

###### Proof

We apply arguments similar to those in the proof of Theorem 3.1. From the
properties of $\nabla\alpha_{\rm f}$ in (3.3), and $\alpha_{\rm c}^{\prime}$
from (3.5), the fact that
$\beta_{\epsilon}(s)s\geq([s]^{-})^{2}/\varepsilon-2\varepsilon K_{\beta}$ by
(4.2), and using the Cauchy–Schwarz, Korn–Poincare (3.8), and trace
inequalities (3.17), similarly to (3.16) we deduce the uniform lower bound

$\int_{\Omega\setminus\Sigma_{t}}\sigma(u)\cdot\epsilon(u)\,dx+\int_{\Sigma_{t}}\bigl{\\{}\nabla\alpha_{\rm
f}([\\![u]\\!]_{\tau_{t}})\cdot[\\![u]\\!]_{\tau_{t}}+[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![u]\\!])\,(\nu_{t}\cdot[\\![u]\\!])\bigr{\\}}\,dS_{x}\\\
-\int_{\Gamma^{\rm N}_{t}}g\cdot u\,dS_{x}\geq K_{\rm
KP}\|u\|^{2}_{H^{1}(\Omega\setminus\Sigma_{t})^{d}}-K_{t{\rm
fc}1}\|u\|_{H^{1}(\Omega\setminus\Sigma_{t})^{d}}-2\varepsilon
K_{\beta}|\Sigma_{t}|,$ (4.6)

where

$K_{t{\rm fc}1}:=\bigl{(}\|g\|_{L^{2}(\Gamma^{\rm N}_{t})^{d}}+(K_{\rm
f1}+K_{\rm c1})\sqrt{2|\Sigma_{t}|}\bigr{)}K_{\rm tr}.$ (4.7)

Therefore, the operator associated to (4.4), denoted following (2.5) by
$\partial^{\varepsilon}_{u}\mathcal{E}:V(\Omega_{t})\mapsto
V(\Omega_{t})^{\star}$, is coercive. We have $\nabla\alpha_{\rm f}$ and
$[\alpha_{\rm c}^{\prime}+\beta_{\varepsilon}]$ are uniformly continuous, and
thus preserve $L^{2}$-convergence, the operator
$\partial^{\varepsilon}_{u}\mathcal{E}$ is weakly continuous in the following
sense. If $u^{n}\rightharpoonup u_{t}$ weakly in
$H^{1}(\Omega\setminus\Sigma_{t})^{d}$ as $n\to\infty$ (hence $u^{n}\to u_{t}$
strongly in $L^{2}(\partial\Omega\cup\Sigma^{\pm}_{t})^{d}$ by compactness),
then for each $u\in V(\Omega_{t})$ the following convergence holds

$\int_{\Omega\setminus\Sigma_{t}}\sigma(u^{n})\cdot\epsilon(u)\,dx+\int_{\Sigma_{t}}\bigl{\\{}\nabla\alpha_{\rm
f}([\\![u^{n}]\\!]_{\tau_{t}})\cdot[\\![u]\\!]_{\tau_{t}}+[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![u^{n}]\\!])\,(\nu_{t}\cdot[\\![u]\\!])\bigr{\\}}\,dS_{x}\\\
\to\int_{\Omega\setminus\Sigma_{t}}\sigma(u_{t})\cdot\epsilon(u)\,dx+\int_{\Sigma_{t}}\bigl{\\{}\nabla\alpha_{\rm
f}([\\![u_{t}]\\!]_{\tau_{t}})\cdot[\\![u]\\!]_{\tau_{t}}+[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![u_{t}]\\!])\,(\nu_{t}\cdot[\\![u]\\!])\bigr{\\}}\,dS_{x}.$

Therefore, applying a Galerkin approximation and the Brouwer fixed point
theorem (see Fra/94 ), a solution to the variational problem (4.4) can be
argued. Its uniqueness under the monotony assumption (that is not $\alpha_{\rm
c}^{\prime}$ in (3.6)), and the boundary value formulation (4.1) can be
derived in a standard way. ∎

Next, for a given observation $z\in H^{1}(\partial\Omega)^{d}$, we consider
the $\varepsilon$-dependent least-squares misfit function from (2.4), where
$u^{\varepsilon}_{t}$ satisfies (4.4):

$\mathcal{J}(u^{\varepsilon}_{t};\Omega_{t})=\frac{1}{2}\int_{\Gamma^{\rm
O}_{t}}|u^{\varepsilon}_{t}-z|^{2}\,dS_{x}+\rho|\Sigma_{t}|.$ (4.8)

From the fundamental theorem of calculus, we have the following
representations

$\displaystyle\nabla\alpha_{\rm
f}([\\![u]\\!]_{\tau_{t}})=\int_{0}^{1}\nabla^{2}\alpha_{\rm
f}([\\![ru]\\!]_{\tau_{t}})\,[\\![u]\\!]_{\tau_{t}}\,dr+\nabla\alpha_{\rm
f}(0),$ $\displaystyle[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![u]\\!])=\int_{0}^{1}[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru]\\!])\,(\nu_{t}\cdot[\\![u]\\!])\,dr+[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](0)$ (4.9)

for differentiable $\nabla\alpha_{\rm f},\alpha_{\rm
c}^{\prime},\beta_{\varepsilon}$. Let us fix a solution $u^{\varepsilon}_{t}$
to the variational equation (4.4). Based on (4) we introduce a quadratic
Lagrangian (compare to $\mathcal{L}^{\varepsilon}$ in (2.6)) linearized around
$u^{\varepsilon}_{t}$

$\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u,v;\Omega_{t})=\frac{1}{2}\int_{\Gamma^{\rm
O}_{t}}|u-z|^{2}\,dS_{x}+\rho|\Sigma_{t}|-\int_{\Omega\setminus\Sigma_{t}}\sigma(u)\cdot\epsilon(v)\,dx+\int_{\Gamma^{\rm
N}_{t}}g\cdot v\,dS_{x}\\\
-\int_{\Sigma_{t}}\Bigl{\\{}\Bigl{(}\int_{0}^{1}\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})\,[\\![u]\\!]_{\tau_{t}}\,dr+\nabla\alpha_{\rm
f}(0)\Bigr{)}\cdot[\\![v]\\!]_{\tau_{t}}\\\ +\Bigl{(}\int_{0}^{1}[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])\,(\nu_{t}\cdot[\\![u]\\!])\,dr+[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](0)\Bigr{)}(\nu_{t}\cdot[\\![v]\\!])\Bigr{\\}}\,dS_{x},$
(4.10)

and a saddle point problem corresponding to (2.11): for all $(u,v)\in
V(\Omega_{t})^{2}$,

$\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u^{\varepsilon}_{t},v;\Omega_{t})\leq\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u^{\varepsilon}_{t},v^{\varepsilon}_{t};\Omega_{t})\leq\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u,v^{\varepsilon}_{t};\Omega_{t}).$
(4.11)

Then (4.8) can be expressed equivalently in the primal-dual form (2.11) as

$\mathcal{J}(u^{\varepsilon}_{t};\Omega_{t})=\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u^{\varepsilon}_{t},v^{\varepsilon}_{t};\Omega_{t}),$
(4.12)

where according to (4) the optimal value of the Lagrangian at the solution is

$\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u^{\varepsilon}_{t},v^{\varepsilon}_{t};\Omega_{t})=\frac{1}{2}\int_{\Gamma^{\rm
O}_{t}}|u^{\varepsilon}_{t}-z|^{2}\,dS_{x}+\rho|\Sigma_{t}|-\int_{\Omega\setminus\Sigma_{t}}\sigma(u^{\varepsilon}_{t})\cdot\epsilon(v^{\varepsilon}_{t})\,dx+\int_{\Gamma^{\rm
N}_{t}}g\cdot v^{\varepsilon}_{t}\,dS_{x}\\\
-\int_{\Sigma_{t}}\bigl{\\{}\nabla\alpha_{\rm
f}([\\![u^{\varepsilon}_{t}]\\!]_{\tau_{t}})\cdot[\\![v^{\varepsilon}_{t}]\\!]_{\tau_{t}}+[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!])\,(\nu_{t}\cdot[\\![v^{\varepsilon}_{t}]\\!])\bigr{\\}}\,dS_{x}.$
(4.13)

######  Theorem 4.2 (Well-posedness of $\varepsilon$-regularized saddle-point
problem)

Assume that the cohesion is small in the sense that constant $K_{\rm
f2},K_{\rm c2}$ in (3.3), (3.5) are sufficiently small so that

$K_{\rm fc2}:=K_{\rm KP}-(K_{\rm f2}+K_{\rm c2})2K_{\rm tr}^{2}>0,$ (4.14)

where $K_{\rm KP}$, $K_{\rm tr}$ are from (3.8), (3.17). Then there exists a
unique saddle-point $(u^{\varepsilon}_{t},v^{\varepsilon}_{t})\in
V(\Omega_{t})^{2}$ in (4.11). Its primal component $u^{\varepsilon}_{t}$
solves (4.4). The dual component $v^{\varepsilon}_{t}$ is a solution to the
adjoint equation corresponding to fixed $u^{\varepsilon}_{t}$:

$\langle
A_{\varepsilon}(u^{\varepsilon}_{t})v,v^{\varepsilon}_{t}\rangle:=\int_{\Omega\setminus\Sigma_{t}}\sigma(v)\cdot\epsilon(v^{\varepsilon}_{t})\,dx+\int_{\Sigma_{t}}\int_{0}^{1}\Bigl{\\{}\Bigl{(}\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})\,[\\![v]\\!]_{\tau_{t}}\Bigr{)}\cdot[\\![v^{\varepsilon}_{t}]\\!]_{\tau_{t}}\\\
+[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])\,(\nu_{t}\cdot[\\![v]\\!])(\nu_{t}\cdot[\\![v^{\varepsilon}_{t}]\\!])\Bigr{\\}}\,dr\,dS_{x}=\int_{\Gamma^{\rm
O}_{t}}(u^{\varepsilon}_{t}-z)\cdot v\,dS_{x}$ (4.15)

for all test functions $v\in V(\Omega_{t})$. For smooth solutions the boundary
value relations hold:

$\displaystyle{\rm div}\,\sigma(v^{\varepsilon}_{t})=0$ $\displaystyle\text{
in }\Omega\setminus\Sigma_{t},$ $\displaystyle v^{\varepsilon}_{t}=0\text{ on
}\Gamma^{\rm
D}_{t},\quad\sigma(v^{\varepsilon}_{t})n=u^{\varepsilon}_{t}-z\text{ on
}\Gamma^{\rm O}_{t},\quad\sigma(v^{\varepsilon}_{t})n=0$ $\displaystyle\text{
on }\Gamma^{\rm N}_{t}\setminus\Gamma^{\rm O}_{t},$
$\displaystyle[\\![\sigma(v^{\varepsilon}_{t})\nu_{t}]\\!]=0,\quad(\sigma(v^{\varepsilon}_{t})\nu_{t})_{\tau_{t}}=\int_{0}^{1}\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})\,[\\![v^{\varepsilon}_{t}]\\!]_{\tau_{t}}\,dr,$
$\displaystyle\nu_{t}\cdot(\sigma(v^{\varepsilon}_{t})\nu_{t})=\int_{0}^{1}[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])\,(\nu_{t}\cdot[\\![v^{\varepsilon}_{t}]\\!])\,dr$
$\displaystyle\text{ on }\Sigma_{t}$ (4.16)

implying linear, Robin-type boundary conditions at the interface.

###### Proof

The saddle-point problem consists of two sub-problems: the former and the
latter inequalities in (4.11). Since the Lagrangian
$\tilde{\mathcal{L}}^{\varepsilon}$ from (4.10) is linear in $v$, the primal
maximization problem (the former inequality in (4.11)) is equivalent to the
first order optimality condition (4.4). Its solvability is proven in Theorem
4.1. Since $\tilde{\mathcal{L}}^{\varepsilon}$ from (4.10) is quadratic and
convex in $u$, the dual minimization problem (the latter inequality in (4.11))
is the optimality condition expressed by the adjoint equation (4.15).

Now we prove the solution existence for (4.15). For fixed
$u^{\varepsilon}_{t}$, the left-hand side of (4.15) forms a linear continuous
operator $A_{\varepsilon}(u^{\varepsilon}_{t}):V(\Omega_{t})\mapsto
V^{\star}(\Omega_{t})$. Indeed, using the Cauchy–Schwarz inequality and the
upper bounds for $\nabla^{2}\alpha_{\rm f}$, $\alpha_{\rm c}^{\prime\prime}$,
$\beta^{\prime}_{\varepsilon}$ in (3.3), (3.5), (4.1), the operator is bounded
from above, hence continuous. Recalling the symmetry of the elasticity
coefficients $C$ and the Hessian matrix $\nabla^{2}\alpha_{\rm f}$, the
operator is self-adjoint. Applying the Cauchy–Schwarz, Korn–Poincare (3.8) and
trace inequalities (3.17), due to the boundedness of $\nabla^{2}\alpha_{\rm
f}$, $\alpha_{\rm c}^{\prime\prime}$, $\beta^{\prime}_{\varepsilon}\geq 0$ in
(3.3), (3.5), (4.1), similarly to (4.6), we estimate uniformly from below

$\langle A_{\varepsilon}(u^{\varepsilon}_{t})u,u\rangle\geq K_{\rm
KP}\|u\|^{2}_{H^{1}(\Omega\setminus\Sigma_{t})^{d}}-\int_{\Sigma_{t}}\bigl{\\{}K_{\rm
f2}\bigl{|}[\\![u]\\!]_{\tau_{t}}\bigr{|}^{2}+K_{\rm
c2}\bigl{|}\nu_{t}\cdot[\\![u]\\!]\bigr{|}^{2}\bigr{\\}}\,dS_{x}\\\ \geq
K_{\rm fc2}\|u\|^{2}_{H^{1}(\Omega\setminus\Sigma_{t})^{d}}.$ (4.17)

Here $K_{\rm fc2}>0$ due to assumption (4.14). In this case,
$A_{\varepsilon}(u^{\varepsilon}_{t})$ is uniformly positive. Because
$\nabla^{2}\alpha_{\rm f}$ and $[\alpha_{\rm
c}^{\prime\prime}+\beta_{\varepsilon}^{\prime}]$ are assumed uniformly
continuous, they preserve $L^{2}$-convergence, and the operator
$A_{\varepsilon}(u^{\varepsilon}_{t})$ is weakly lower semi-continuous by the
compactness similar to arguments presented in the proof of Theorem 4.1.
According to the Lax–Milgram theorem, the variational equation (4.15) has a
unique solution. We derive straightforwardly its boundary value formulation
(4.2).

Since the variational equation (4.4) can be rewritten in the equivalent form

$\langle
A_{\varepsilon}(u^{\varepsilon}_{t})u^{\varepsilon}_{t},u\rangle+\int_{\Sigma_{t}}\bigl{(}\nabla\alpha_{\rm
f}(0)\cdot[\\![u]\\!]_{\tau_{t}}+[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](0)(\nu_{t}\cdot[\\![u]\\!])\bigr{)}\,dS_{x}=\int_{\Gamma^{\rm
N}_{t}}g\cdot u\,dS_{x}$

for all $u\in V(\Omega_{t})$, by assumption (4.14) its solution
$u^{\varepsilon}_{t}$ is unique, too. ∎

##  5 Shape derivative

Let us fix a flow and its inverse

$s\mapsto(\phi_{s},\phi^{-1}_{s})\in
C^{1}([t_{0}-t_{1},t_{1}-t_{0}];W^{1,\infty}(\overline{\Omega})^{d})^{2}.$
(5.1)

This defines an associated coordinate transformation $y=\phi_{s}(x)$ and its
inverse $x=\phi^{-1}_{s}(y)$. For every fixed $t\in(t_{0},t_{1})$, we suppose
that for $s\in[t_{0},t_{1}]-t$ it forms a diffeomorphism

$\phi_{s}:\Omega_{t}\mapsto\Omega_{t+s},\;x\mapsto
y,\quad\phi^{-1}_{s}:\Omega_{t+s}\mapsto\Omega_{t},\;y\mapsto x,$ (5.2)

where the perturbed geometry $\Omega_{t+s}=(\Gamma^{\rm D}_{t+s},\Gamma^{\rm
N}_{t+s},\Gamma^{\rm O}_{t+s},\Sigma_{t+s})$ describes the broken domain
$\Omega\setminus\Sigma_{t+s}$. From (5.1), a time-dependent kinematic velocity
$\Lambda(t,x)\in C([t_{0},t_{1}];W^{1,\infty}(\overline{\Omega})^{d})$ is
assumed defined by the formula

$\Lambda(t+s,y):={\textstyle\frac{d}{ds}}\phi_{s}(\phi^{-1}_{s}(y)).$ (5.3)

If a stationary velocity is given explicitly by $\Lambda(x)\in
W^{1,\infty}(\overline{\Omega})^{d}$ with $n\cdot\Lambda=0$ at
$\partial\Omega$, thus preserving the hold-all domain, then $\Lambda$
determines the flow (5.1) by unique solutions to the autonomous ODE systems:

$\left\\{\begin{array}[]{rl}{\textstyle\frac{d}{ds}}\phi_{s}=\Lambda(\phi_{s})&\;\text{for
$s\not=0$},\\\ \phi_{s}=x&\;\text{for
$s=0$},\end{array}\right.\quad\left\\{\begin{array}[]{rl}{\textstyle\frac{d}{ds}}\phi^{-1}_{s}=-\Lambda(\phi^{-1}_{s})&\;\text{for
$s\not=0$},\\\ \phi^{-1}_{s}=y&\;\text{for $s=0$},\end{array}\right.$ (5.4)

which build a semi-group of transformations.

The following properties (T1)–(T4) are needed to prove shape
differentiability.

(T1)

We assume that the map $u\mapsto u\circ\phi_{s}$ is bijective between the
function spaces

$V(\Omega_{t+s})\mapsto V(\Omega_{t}).$ (5.5)

Based on assumption (5.5), the perturbed objective $(t_{0}-t,t_{1}-t)\times
V(\Omega_{t})$, $(s,\tilde{u})\mapsto\tilde{\mathcal{J}}$ and Lagrangian
$(t_{0}-t,t_{1}-t)\times V(\Omega_{t})^{2}$,
$(s,\tilde{u},\tilde{v})\mapsto\tilde{\mathcal{L}}^{\varepsilon}$, are well-
defined for $(u,v)\in V(\Omega_{t+s})^{2}$ when transformed to the reference
geometry $\Omega_{t}$ by setting

$\tilde{\mathcal{J}}(s,u\circ\phi_{s};\Omega_{t})=\mathcal{J}(u;\Omega_{t+s}),\quad\tilde{\mathcal{L}}^{\varepsilon}(s,u\circ\phi_{s},u\circ\phi_{s},v\circ\phi_{s};\Omega_{t})=\mathcal{L}^{\varepsilon}(u,v;\Omega_{t+s}).$
(5.6)

At $s=0$ relations (5.6) imply that for $(\tilde{u},\tilde{v})\in
V(\Omega_{t})^{2}$

$\tilde{\mathcal{J}}(0,\tilde{u};\Omega_{t})=\mathcal{J}(\tilde{u};\Omega_{t}),\quad\tilde{\mathcal{L}}^{\varepsilon}(0,\tilde{u},\tilde{u},\tilde{v};\Omega_{t})=\mathcal{L}^{\varepsilon}(\tilde{u},\tilde{v};\Omega_{t}).$
(5.7)

According to (2.11) we look for a saddle-point
$(\tilde{u}^{\varepsilon}_{t+s},\tilde{v}^{\varepsilon}_{t+s})\in
V(\Omega_{t})^{2}$ satisfying the inequalities

$\tilde{\mathcal{L}}^{\varepsilon}(s,u^{\varepsilon}_{t},\tilde{u}^{\varepsilon}_{t+s},\tilde{v};\Omega_{t})\leq\tilde{\mathcal{L}}^{\varepsilon}(s,u^{\varepsilon}_{t},\tilde{u}^{\varepsilon}_{t+s},\tilde{v}^{\varepsilon}_{t+s};\Omega_{t})\leq\tilde{\mathcal{L}}^{\varepsilon}(s,u^{\varepsilon}_{t},\tilde{u},\tilde{v}^{\varepsilon}_{t+s};\Omega_{t})$
(5.8)

for all $(\tilde{u},\tilde{v})\in V(\Omega_{t})^{2}$. In the case of
$\mathcal{J}$ from (4.8), applying the coordinate transformation (5.2) we
derive explicitly the objective function

$\tilde{\mathcal{J}}(s,\tilde{u};\Omega_{t})=\frac{1}{2}\int_{\Gamma^{\rm
O}_{t}}|\tilde{u}-z\circ\phi_{s}|^{2}\,\omega^{\rm
b}_{s}dS_{x}+\rho\int_{\Sigma_{t}}\omega^{\rm b}_{s}\,dS_{x},$ (5.9)

where

$\omega^{\rm d}_{s}:=\det(\nabla\phi_{s})\text{ in
}\Omega\setminus\Sigma_{t},\quad\omega^{\rm
b}_{s}:=|(\nabla\phi_{s}^{-\top}\circ\phi_{s})n^{\pm}_{t}|\omega^{\rm
d}_{s}\text{ at }\partial\Omega^{\pm}_{t}$ (5.10)

denote the Jacobians, and set the perturbed Lagrangian according to (4.10) as

$\tilde{\mathcal{L}}^{\varepsilon}(s,u^{\varepsilon}_{t},\tilde{u},\tilde{v};\Omega_{t})=\tilde{\mathcal{J}}(s,\tilde{u};\Omega_{t})\\\
-\int_{\Omega\setminus\Sigma_{t}}\bigl{(}(C\circ\phi_{s})E(\nabla\phi_{s}^{-1}\circ\phi_{s},\tilde{u})\cdot
E(\nabla\phi_{s}^{-1}\circ\phi_{s},\tilde{v})\bigr{)}\,\omega^{\rm
d}_{s}dx+\int_{\Gamma^{\rm N}_{t}}(g\circ\phi_{s})\cdot\tilde{v}\,\omega^{\rm
b}_{s}dS_{x}\\\
-\int_{\Sigma_{t}}\Bigl{\\{}\Bigl{(}\int_{0}^{1}\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})[\\![\tilde{u}]\\!]_{\tilde{\tau}_{t+s}}\,dr+\nabla\alpha_{\rm
f}(0)\Bigr{)}\cdot[\\![\tilde{v}]\\!]_{\tilde{\tau}_{t+s}}\\\
+\Bigl{(}\int_{0}^{1}[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])(\tilde{\nu}_{t+s}\cdot[\\![\tilde{u}]\\!])\,dr+[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](0)\Bigr{)}(\tilde{\nu}_{t+s}\cdot[\\![\tilde{v}]\\!])\Bigr{\\}}\,\omega^{\rm
b}_{s}dS_{x}.$ (5.11)

In (5.11), the following decomposition at $\Sigma_{t}$ was used in accordance
with (3.1):

$[\\![\tilde{u}]\\!]_{\tilde{\tau}_{t+s}}:=[\\![\tilde{u}]\\!]-(\tilde{\nu}_{t+s}\cdot[\\![\tilde{u}]\\!])\,\tilde{\nu}_{t+s},\quad\tilde{\nu}_{t+s}:=\nu_{t+s}\circ\phi_{s}.$
(5.12)

Further in view of the chain rule
$\nabla_{y}u=(\nabla\phi_{s}^{-T}\circ\phi_{s})\nabla(u\circ\phi_{s})$, there
appears the expression

$E(M,\tilde{u}):=\frac{1}{2}(M^{\top}\nabla\tilde{u}+\nabla\tilde{u}^{\top}M),\quad
M\in\mathbb{R}^{d\times d},$ (5.13)

for which $E(I,\tilde{u})=\epsilon(\tilde{u})$ according to (3.7). For more
details of the derivation, see Kov/06 ; KK/07 ; KO/20 .

######  Lemma 2 (T2)

The asymptotic expansion in the first argument of $\tilde{\mathcal{J}}$ from
(5.9) is given by

$\tilde{\mathcal{J}}(s,\tilde{u};\Omega_{t})=\mathcal{J}(\tilde{u};\Omega_{t})+{\rm
O}(|s|),$ (5.14)

and the expansion of $\tilde{\mathcal{L}}^{\varepsilon}$ from (5.11) by:

$\tilde{\mathcal{L}}^{\varepsilon}(s,u^{\varepsilon}_{t},\tilde{u},\tilde{v};\Omega_{t})=\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},\tilde{u},\tilde{v};\Omega_{t})+s{\textstyle\frac{\partial}{\partial
s}}\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},\tilde{u},\tilde{v};\Omega_{t})+{\rm
o}(|s|)$ (5.15)

holds as $s\to 0$. The partial derivative
$(t_{0},t_{1})-t\mapsto\mathbb{R},\tau\mapsto{\textstyle\frac{\partial}{\partial
s}}\tilde{\mathcal{L}}^{\varepsilon}$ in (5.15) is a continuous function and
exhibits the explicit representation

${\textstyle\frac{\partial}{\partial
s}}\tilde{\mathcal{L}}^{\varepsilon}(\tau,u^{\varepsilon}_{t},\tilde{u},\tilde{v};\Omega_{t})=\int_{\Gamma^{\rm
O}_{t}}\Bigl{(}\frac{1}{2}{\rm
div}_{\tau_{t}}\Lambda|_{t+\tau}|\tilde{u}-z|^{2}-\nabla
z\Lambda|_{t+\tau}\cdot(\tilde{u}-z)\Bigr{)}dS_{x}+\rho\int_{\Sigma_{t}}{\rm
div}_{\tau_{t}}\Lambda|_{t+\tau}dS_{x}\\\
-\int_{\Omega\setminus\Sigma_{t}}\Bigl{(}\bigl{(}{\rm
div}\Lambda|_{t+\tau}C+\nabla
C\Lambda|_{t+\tau}\bigr{)}\epsilon(\tilde{u})\cdot\epsilon(\tilde{v})-\sigma(\tilde{u})\cdot
E(\nabla\Lambda|_{t+\tau},\tilde{v})-\sigma(\tilde{v})\cdot
E(\nabla\Lambda|_{t+\tau},\tilde{u})\Bigr{)}dx\\\ +\int_{\Gamma^{\rm
N}_{t}}\bigl{(}{\rm div}_{\tau_{t}}\Lambda|_{t+\tau}g+\nabla
g\,\Lambda|_{t+\tau}\bigr{)}\cdot\tilde{v}\,dS_{x}-\int_{\Sigma_{t}}\Bigl{\\{}\int_{0}^{1}\bigl{(}\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})\,[\\![\tilde{u}]\\!]_{\nabla\tau_{t}\Lambda|_{t+\tau}}\bigr{)}\cdot[\\![\tilde{v}]\\!]_{\tau_{t}}dr\\\
+\Bigl{(}\int_{0}^{1}\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})\,[\\![\tilde{u}]\\!]_{\tau_{t}}dr+\nabla\alpha_{\rm
f}(0)\Bigr{)}\cdot\bigl{(}{\rm
div}_{\tau_{t}}\Lambda|_{t+\tau}[\\![\tilde{v}]\\!]_{\tau_{t}}+[\\![\tilde{v}]\\!]_{\nabla\tau_{t}\Lambda|_{t+\tau}}\bigr{)}\\\
+\Bigl{(}\int_{0}^{1}[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])\,(\nu_{t}\cdot[\\![\tilde{u}]\\!])\,dr+[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](0)\Bigr{)}\bigl{(}({\rm
div}_{\tau_{t}}\Lambda|_{t+\tau}\nu_{t}+\nabla\nu_{t}\Lambda|_{t+\tau})\cdot[\\![\tilde{v}]\\!]\bigr{)}\\\
+\int_{0}^{1}[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])\,\bigl{(}\nabla\nu_{t}\Lambda|_{t+\tau}\cdot[\\![\tilde{u}]\\!]\bigr{)}(\nu_{t}\cdot[\\![\tilde{v}]\\!])\,dr\Bigr{\\}}\,dS_{x}.$
(5.16)

In (5.16) the notation $\nabla\tau_{t}\Lambda$ and $\nabla\nu_{t}\Lambda$ at
$\Sigma_{t}$ stands for

$[\\![\tilde{u}]\\!]_{\nabla\tau_{t}\Lambda}:=-(\nu_{t}\cdot[\\![\tilde{u}]\\!])\nabla\nu_{t}\Lambda-(\nabla\nu_{t}\Lambda\cdot[\\![\tilde{u}]\\!])\nu_{t},\quad\nabla\nu_{t}\Lambda:=\bigl{(}(\nabla\Lambda\,\nu_{t})\cdot\nu_{t}\bigr{)}\nu_{t}-\nabla\Lambda^{\top}\nu_{t},$
(5.17)

and the tangential divergence is defined as

${\rm div}_{\tau_{t}}\Lambda={\rm
div}\Lambda-(\nabla\Lambda\,n^{\pm}_{t})\cdot n^{\pm}_{t}\text{ at
}\partial\Omega^{\pm}_{t}.$ (5.18)

The proof of Lemma 2 is presented in Appendix A.

######  Lemma 3 (T3)

The set of saddle points
$(\tilde{u}^{\varepsilon}_{t+s},\tilde{v}^{\varepsilon}_{t+s})$ for (5.8) is a
singleton for all $s\in[t_{0},t_{1}]-t$, and
$(\tilde{u}^{\varepsilon}_{t},\tilde{v}^{\varepsilon}_{t})=(u^{\varepsilon}_{t},v^{\varepsilon}_{t})$
as $s=0$.

The proof is given in in Appendix B and follows the arguments in the proof of
Theorem 4.2, which treats a particular case of the saddle-point problem (5.8)
as $s=0$.

######  Lemma 4 (T4)

There exists a subsequence $s_{k}\to 0$ as $k\to\infty$, such that

$(\tilde{u}^{\varepsilon}_{t+s_{k}},\tilde{v}^{\varepsilon}_{t+s_{k}})\to(u^{\varepsilon}_{t},v^{\varepsilon}_{t})\quad\text{strongly
in $V(\Omega_{t})^{2}$ as $s_{k}\to 0$}.$ (5.19)

The proof of Lemma 4 is technical. It is presented in Appendix C.

Based on the properties (T1)–(T4) we establish the main result of this
section.

######  Theorem 5.1 (Shape differentiability of $\varepsilon$-regularized
optimization problem)

Under assumption (4.14), the shape derivative (see its definition (2.10) and
existence criterion (2.12)) can be expressed by the partial derivative from
(5.16) as

$\partial_{t}\mathcal{J}(u^{\varepsilon}_{t};\Omega_{t})={\textstyle\frac{\partial}{\partial
s}}\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u^{\varepsilon}_{t},v^{\varepsilon}_{t};\Omega_{t})=\int_{\Gamma^{\rm
O}_{t}}\Bigl{(}\frac{1}{2}{\rm
div}_{\tau_{t}}\Lambda\,|u^{\varepsilon}_{t}-z|^{2}-\nabla
z\Lambda\cdot(u^{\varepsilon}_{t}-z)\Bigr{)}dS_{x}\\\
-\int_{\Omega\setminus\Sigma_{t}}\Bigl{(}\bigl{(}{\rm div}\Lambda\,C+\nabla
C\Lambda\bigr{)}\epsilon(u^{\varepsilon}_{t})\cdot\epsilon(v^{\varepsilon}_{t})-\sigma(u^{\varepsilon}_{t})\cdot
E(\nabla\Lambda,v^{\varepsilon}_{t})-\sigma(v^{\varepsilon}_{t})\cdot
E(\nabla\Lambda,u^{\varepsilon}_{t})\Bigr{)}dx\\\ +\int_{\Gamma^{\rm
N}_{t}}\bigl{(}{\rm div}_{\tau_{t}}\Lambda\,g+\nabla g\,\Lambda\bigr{)}\cdot
v^{\varepsilon}_{t}\,dS_{x}-\int_{\Sigma_{t}}\Bigl{\\{}\nabla\alpha_{\rm
f}([\\![u^{\varepsilon}_{t}]\\!]_{\tau_{t}})\cdot\bigl{(}{\rm
div}_{\tau_{t}}\Lambda\,[\\![v^{\varepsilon}_{t}]\\!]_{\tau_{t}}+[\\![v^{\varepsilon}_{t}]\\!]_{\nabla\tau_{t}\Lambda}\bigr{)}\\\
+\int_{0}^{1}\bigl{(}\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})\,[\\![u^{\varepsilon}_{t}]\\!]_{\nabla\tau_{t}\Lambda}\bigr{)}\cdot[\\![v^{\varepsilon}_{t}]\\!]_{\tau_{t}}dr+[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!])\bigl{(}({\rm
div}_{\tau_{t}}\Lambda\,\nu_{t}+\nabla\nu_{t}\Lambda)\cdot[\\![v^{\varepsilon}_{t}]\\!]\bigr{)}\\\
+\int_{0}^{1}[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])\,\bigl{(}\nabla\nu_{t}\Lambda\cdot[\\![u^{\varepsilon}_{t}]\\!]\bigr{)}(\nu_{t}\cdot[\\![v^{\varepsilon}_{t}]\\!])\,dr\Bigr{\\}}\,dS_{x}+\rho\int_{\Sigma_{t}}{\rm
div}_{\tau_{t}}\Lambda\,dS_{x},$ (5.20)

where $(u^{\varepsilon}_{t},v^{\varepsilon}_{t})\in V(\Omega_{t})^{2}$ is a
saddle-point to (4.11).

###### Proof

Indeed, due to (T1)–(T4) all assumptions in Delfour–Zolesio (DZ/11, , Chapter
10, Theorem 5.1) are satisfied. Details of the proof can be found in KO/20 . ∎

######  Corollary 1 (Hadamard representation of the $\varepsilon$-dependent
shape derivative)

Assume that the solution of (4.4) and (4.15) satisfies
$(u^{\varepsilon}_{t},v^{\varepsilon}_{t})\in H^{2}(\Omega^{+}_{t})^{2d}\cap
H^{2}(\Omega^{-}_{t})^{2d}$. Introducing the decomposition into normal and
tangential components according to

$\Lambda=(n_{t}\cdot\Lambda)n_{t}+\Lambda_{\tau_{t}},\quad\nabla=(n_{t}\cdot\nabla)n_{t}+\nabla_{\tau_{t}},\quad\mathcal{D}=(n_{t}\cdot\mathcal{D})n_{t}+\mathcal{D}_{\tau_{t}},$
(5.21)

the following equivalent representation of the shape derivative (5.20) holds
in terms of boundary integrals in 2D:

${\textstyle\frac{\partial}{\partial
s}}\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u^{\varepsilon}_{t},v^{\varepsilon}_{t};\Omega_{t})\\\
=\int_{\Gamma^{\rm
D}_{t}}(\tau_{t}\cdot\Lambda)\tau_{t}\cdot\mathcal{D}_{1}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})\,dS_{x}+\int_{\Sigma_{t}}\bigl{(}(\tau_{t}\cdot\Lambda)\tau_{t}\cdot\mathcal{D}^{\varepsilon}_{2}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})+(\nu_{t}\cdot\Lambda)\mathcal{D}^{\varepsilon}_{3}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})\bigr{)}\,dS_{x}\\\
+(\tau_{t}\cdot\Lambda)[\\![\mathcal{D}^{\varepsilon}_{4}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})]\\!]_{\partial\Sigma_{t}}+(\tau_{t}\cdot\Lambda)\mathcal{D}_{5}(u^{\varepsilon}_{t})|_{\partial\Gamma^{\rm
O}_{t}}+(\tau_{t}\cdot\Lambda)[\\![\mathcal{D}_{6}(v^{\varepsilon}_{t})]\\!]_{\partial\Gamma^{\rm
N}_{t}\cap\Sigma_{t}},$ (5.22)

where $\tau_{t}$ is a tangential vector at the boundary, and in 3D:

$=\int_{\Gamma^{\rm
D}_{t}}\Lambda_{\tau_{t}}\cdot\mathcal{D}_{1}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})_{\tau_{t}}dS_{x}+\int_{\Sigma_{t}}\bigl{(}\Lambda_{\tau_{t}}\cdot\mathcal{D}^{\varepsilon}_{2}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})_{\tau_{t}}+(\nu_{t}\cdot\Lambda)\mathcal{D}^{\varepsilon}_{3}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})\bigr{)}\,dS_{x}\\\
+{\displaystyle\int_{\partial\Sigma_{t}}}\\!\\!\\!\\!(b_{t}\cdot\Lambda)[\\![\mathcal{D}^{\varepsilon}_{4}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})]\\!]dL_{x}+{\displaystyle\int_{\partial\Gamma^{\rm
O}_{t}}}\\!\\!\\!\\!(b_{t}\cdot\Lambda)\mathcal{D}_{5}(u^{\varepsilon}_{t})dL_{x}+{\displaystyle\int_{\partial\Gamma^{\rm
N}_{t}\cap\Sigma_{t}}}\\!\\!\\!\\!(b_{t}\cdot\Lambda)[\\![\mathcal{D}_{6}(v^{\varepsilon}_{t})]\\!]dL_{x},$
(5.23)

where $b_{t}=\tau_{t}\times n_{t}$ is a binomial vector within the moving
frame at the respective boundary. The terms in (5.22) and (5.23) are

$\mathcal{D}_{1}(\tilde{u},\tilde{v}):=\nabla\tilde{u}^{\top}\sigma(\tilde{v})n_{t}+\nabla\tilde{v}^{\top}\sigma(\tilde{u})n_{t},\quad\mathcal{D}^{\varepsilon}_{2}(\tilde{u},\tilde{v}):=-[q_{\rm
f}+q^{\varepsilon}_{\rm c}]_{\tau_{t}}(\tilde{u},\tilde{v}),\\\
\mathcal{D}^{\varepsilon}_{3}(\tilde{u},\tilde{v}):=[\\![\sigma(\tilde{u})\cdot\epsilon(\tilde{v})]\\!]+\rho\varkappa_{t}-\varkappa_{t}[p_{\rm
f}+p^{\varepsilon}_{\rm c}](\tilde{u},\tilde{v})-\nu_{t}\cdot[\nabla(p_{\rm
f}+p^{\varepsilon}_{\rm c})+q_{\rm f}+q^{\varepsilon}_{\rm
c}](\tilde{u},\tilde{v}),\\\
\mathcal{D}^{\varepsilon}_{4}(\tilde{u},\tilde{v}):=\rho-[p_{\rm
f}+p^{\varepsilon}_{\rm
c}](\tilde{u},\tilde{v}),\quad\mathcal{D}_{5}(\tilde{u}):=\frac{1}{2}|\tilde{u}-z|^{2},\quad\mathcal{D}_{6}(\tilde{v}):=g\cdot\tilde{v},$
(5.24)

with the curvature $\varkappa_{t}={\rm div}_{\tau_{t}}\nu_{t}$ at
$\Sigma_{t}$. The expressions along $\Sigma_{t}$ are defined by

$p_{\rm f}(\tilde{u},\tilde{v}):=\nabla\alpha_{\rm
f}([\\![\tilde{u}]\\!]_{\tau_{t}})\cdot[\\![\tilde{v}]\\!]_{\tau_{t}},\quad
p^{\varepsilon}_{\rm c}(\tilde{u},\tilde{v}):=[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![\tilde{u}]\\!])\,(\nu_{t}\cdot[\\![\tilde{v}]\\!]),$
(5.25)

and next

$q_{\rm
f}(\tilde{u},\tilde{v}):=[\\![\nabla\tilde{v}]\\!]^{\top}\nu_{t}\bigl{(}\nu_{t}\cdot\nabla\alpha_{\rm
f}([\\![\tilde{u}]\\!]_{\tau_{t}})\bigr{)}+[\\![\nabla\tilde{u}]\\!]^{\top}\nu_{t}\Bigl{(}\nu_{t}\cdot\int_{0}^{1}\nabla^{2}\alpha_{\rm
f}([\\![ru^{t}_{\varepsilon}]\\!]_{\tau_{t}})[\\![\tilde{v}]\\!]_{\tau_{t}}dr\Bigr{)}\\\
+\nabla([\\![\tilde{u}]\\!]_{\tau_{t}})^{\top}\int_{0}^{1}\bigl{(}\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})-\nabla^{2}\alpha_{\rm
f}([\\![u^{\varepsilon}_{t}]\\!]_{\tau_{t}})\bigr{)}[\\![\tilde{v}]\\!]_{\tau_{t}}\,dr,\\\
q^{\varepsilon}_{\rm
c}(\tilde{u},\tilde{v}):=\nabla(\nu_{t}\cdot[\\![\tilde{u}]\\!])^{\top}\int_{0}^{1}\bigl{(}[\alpha^{\prime\prime}_{\rm
c}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])-[\alpha^{\prime\prime}_{\rm
c}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!])\bigr{)}(\nu_{t}\cdot[\\![\tilde{v}]\\!])\,dr.$
(5.26)

The proof of Corollary 1 is given in Appendix D.

We remark that the additional $H^{2}$-regularity is available when a piecewise
$C^{2,0}$-boundaries $\partial\Omega^{\pm}_{t}$ exclude singular points (e.g.
in 2D when the boundary parts meet each other with an $\pi/2$-angle as in
Figure 1).

######  Corollary 2 (Descent direction for the $\varepsilon$-dependent
optimization)

A descent direction for the perturbed $\tilde{\mathcal{L}}^{\varepsilon}$ in
(5.15) is provided by the following choice of the velocity

$\tau_{t}\cdot\Lambda=-k_{1}\tau_{t}\cdot\mathcal{D}_{1}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})\text{
at $\Gamma^{\rm
D}_{t}$},\quad\tau_{t}\cdot\Lambda=-k_{2}\tau_{t}\cdot\mathcal{D}^{\varepsilon}_{2}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})\text{
and
}\nu_{t}\cdot\Lambda=-k_{3}\mathcal{D}^{\varepsilon}_{3}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})\text{
at $\Sigma_{t}$},\\\
\tau_{t}\cdot\Lambda=-k_{4}[\\![\mathcal{D}^{\varepsilon}_{4}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})]\\!]\text{
at
$\partial\Sigma_{t}$},\;\tau_{t}\cdot\Lambda=-k_{5}\mathcal{D}_{5}(u^{\varepsilon}_{t})\text{
at $\partial\Gamma^{\rm
O}_{t}$},\;\tau_{t}\cdot\Lambda=-k_{6}[\\![\mathcal{D}_{6}(v^{\varepsilon}_{t})]\\!]\text{at
$\partial\Gamma^{\rm N}_{t}\cap\Sigma_{t}$},\\\ n_{t}\cdot\Lambda=0\text{ at
$\partial\Omega$}$ (5.27)

in 2D, and in 3D respectively

$\Lambda_{\tau_{t}}=-k_{1}\mathcal{D}_{1}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})_{\tau_{t}}\text{
at $\Gamma^{\rm
D}_{t}$},\quad\Lambda_{\tau_{t}}=-k_{2}\mathcal{D}^{\varepsilon}_{2}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})_{\tau_{t}}\text{and
}\nu_{t}\cdot\Lambda=-k_{3}\mathcal{D}^{\varepsilon}_{3}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})\text{
at $\Sigma_{t}$},\\\
b_{t}\cdot\Lambda=-k_{4}[\\![\mathcal{D}^{\varepsilon}_{4}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})]\\!]\text{
at
$\partial\Sigma_{t}$},\;b_{t}\cdot\Lambda=-k_{5}\mathcal{D}_{5}(u^{\varepsilon}_{t})\text{
at $\partial\Gamma^{\rm
O}_{t}$},\;b_{t}\cdot\Lambda=-k_{6}[\\![\mathcal{D}_{6}(v^{\varepsilon}_{t})]\\!]\text{at
$\partial\Gamma^{\rm N}_{t}\cap\Sigma_{t}$},\\\ n_{t}\cdot\Lambda=0\text{ at
$\partial\Omega$},$ (5.28)

with $k_{i}\geq 0$, $i=1,\ldots,6$, and not all simultaneously equal to zero.

###### Proof

Direct substitution of (5.27) into (5.22) in 2D, respectively (5.28) into
(5.23) in 3D, provides that ${\textstyle\frac{\partial}{\partial
s}}\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u^{\varepsilon}_{t},v^{\varepsilon}_{t};\Omega_{t})<0$.
∎

Corollary 2 is of practical importance since it provides well-posedness of
gradient schemes (see Algorithm 1) based on the descent direction from (5.27)
and (5.28).

##  6 The limit as $\varepsilon\to 0^{+}$

In the following we derive the limit relations as $\varepsilon\to 0^{+}$. We
recall that all results involving the dual variable $v^{\varepsilon}_{t}$
assume that (4.14) holds true.

######  Lemma 5 (Uniform estimate)

The following a-priori estimate holds uniformly in
$\varepsilon\in(0,\varepsilon_{0})$:

$\|u^{\varepsilon}_{t}\|_{H^{1}(\Omega\setminus\Sigma_{t})^{d}}+\frac{1}{\sqrt{\varepsilon}}\|\bigl{[}\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!]\bigr{]}^{-}\|_{L^{2}(\Sigma_{t})}+\|v^{\varepsilon}_{t}\|_{H^{1}(\Omega\setminus\Sigma_{t})^{d}}\leq
K,\quad K\geq 0.$ (6.1)

Consequently, there exists a subsequence $\varepsilon_{k}\to 0$ as
$k\to\infty$ and an accumulation point $(u_{t},v_{t})\in K(\Omega_{t})\times
V(\Omega_{t})$ such that

$(u^{\varepsilon_{k}}_{t},v^{\varepsilon_{k}}_{t})\to(u_{t},v_{t})\text{
weakly in $H^{1}(\Omega\setminus\Sigma_{t})^{2d}$,
$H^{1/2}(\partial\Omega^{\pm}_{t})^{2d}$, strongly in
$L^{2}(\partial\Omega^{\pm}_{t})^{2d}$}.$ (6.2)

###### Proof

Passing $s\to 0$ due to the convergences (C.6) and (C.7) and using the lower
bound
$\beta_{\epsilon}(\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!])(\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!])\geq([\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!]]^{-})^{2}/\varepsilon-2\varepsilon
K_{\beta}$ due to (4.2), in the limit we improve the uniform a-priori estimate
(C.5) and get (6.1). Consequently (6.2) follows by a standard compactness
argument. Moreover, $[\nu_{t}\cdot[\\![u^{\varepsilon_{k}}_{t}]\\!]]^{-}\to 0$
ensures $\nu_{t}\cdot[\\![u_{t}]\\!]\geq 0$ at $\Sigma_{t}$, hence $u_{t}\in
K(\Omega_{t})$. ∎

Let $u_{t}\in K(\Omega_{t})$ be a solution to the VI (3.14) in Theorem 3.1.
According to (4.10) we introduce the $\varepsilon$-independent Lagrangian
$(u,v)\mapsto\mathcal{L}:V(\Omega_{t})^{2}\mapsto\mathbb{R}$ as

$\mathcal{L}(u_{t},u,v;\Omega_{t}):=\frac{1}{2}\int_{\Gamma^{\rm
O}_{t}}|u-z|^{2}\,dS_{x}+\rho|\Sigma_{t}|-\int_{\Omega\setminus\Sigma_{t}}\sigma(u)\cdot\epsilon(v)\,dx+\int_{\Gamma^{\rm
N}_{t}}g\cdot v\,dS_{x}\\\
-\int_{\Sigma_{t}}\Bigl{\\{}\Bigl{(}\int_{0}^{1}\nabla^{2}\alpha_{\rm
f}([\\![ru_{t}]\\!]_{\tau_{t}})\,[\\![u]\\!]_{\tau_{t}}dr+\nabla\alpha_{\rm
f}(0)\Bigr{)}\cdot[\\![v]\\!]_{\tau_{t}}\\\ +\Bigl{(}\int_{0}^{1}\alpha_{\rm
c}^{\prime\prime}(\nu_{t}\cdot[\\![ru_{t}]\\!])\,(\nu_{t}\cdot[\\![u]\\!])\,dr+\alpha_{\rm
c}^{\prime}(0)\Bigr{)}(\nu_{t}\cdot[\\![v]\\!])\Bigr{\\}}\,dS_{x}.$ (6.3)

Based on Lemma 5 we prove the following.

######  Theorem 6.1 (Limit optimality conditions)

(i) There exists a pair $(u_{t},\lambda_{t})\in V(\Omega_{t})\times
H^{1/2}(\Sigma_{t})^{\star}$ which satisfies the variational equation

$\int_{\Omega\setminus\Sigma_{t}}\sigma(u_{t})\cdot\epsilon(u)\,dx+\int_{\Sigma_{t}}\bigl{\\{}\nabla\alpha_{\rm
f}([\\![u_{t}]\\!]_{\tau_{t}})\cdot[\\![u]\\!]_{\tau_{t}}+\alpha_{\rm
c}^{\prime}(\nu_{t}\cdot[\\![u_{t}]\\!])\,(\nu_{t}\cdot[\\![u]\\!])\bigr{\\}}\,dS_{x}\\\
+\langle\lambda_{t},\nu_{t}\cdot[\\![u]\\!]\rangle_{\Sigma_{t}}=\int_{\Gamma^{\rm
N}_{t}}g\cdot u\,dS_{x}$ (6.4)

for all test functions $u\in V(\Omega_{t})$, simultaneously with the
complementary relations

$\nu_{t}\cdot[\\![u_{t}]\\!]\geq 0,\quad\lambda_{t}\leq
0,\quad\langle\lambda_{t},\nu_{t}\cdot[\\![u_{t}]\\!]\rangle_{\Sigma_{t}}=0,$
(6.5)

where $\langle\,\cdot\,,\,\cdot\,\rangle_{\Sigma_{t}}$ stands for the duality
pairing between $H^{1/2}(\Sigma_{t})$ and its dual space
$H^{1/2}(\Sigma_{t})^{\star}$. The first component $u_{t}\in K(\Omega_{t})$
solves the VI (3.14), and according to (3.1) the second, $\lambda_{t}$,
satisfies

$\lambda_{t}=\nu_{t}\cdot(\sigma(u_{t})\nu_{t})-\alpha_{\rm
c}^{\prime}(\nu_{t}\cdot[\\![u_{t}]\\!])\quad\text{at $\Sigma_{t}$}.$ (6.6)

(ii) Under the assumption (4.14), an adjoint pair $(v_{t},\mu_{t})\in
V(\Omega_{t})\times H^{1/2}(\Sigma_{t})^{\star}$ exists and satisfies the
adjoint equation

$\int_{\Omega\setminus\Sigma_{t}}\sigma(v)\cdot\epsilon(v_{t})\,dx+\int_{\Sigma_{t}}\int_{0}^{1}\Bigl{\\{}\Bigl{(}\nabla^{2}\alpha_{\rm
f}([\\![ru_{t}]\\!]_{\tau_{t}})\,[\\![v]\\!]_{\tau_{t}}\Bigr{)}\cdot[\\![v_{t}]\\!]_{\tau_{t}}\\\
+\alpha_{\rm
c}^{\prime\prime}(\nu_{t}\cdot[\\![ru_{t}]\\!])\,(\nu_{t}\cdot[\\![v]\\!])(\nu_{t}\cdot[\\![v_{t}]\\!])\Bigr{\\}}\,dr\,dS_{x}+\langle\mu_{t},\nu_{t}\cdot[\\![v]\\!]\rangle_{\Sigma_{t}}=\int_{\Gamma^{\rm
O}_{t}}(u_{t}-z)\cdot v\,dS_{x}$ (6.7)

for all test functions $v\in V(\Omega_{t})$, such that the compatibility
relation holds:

$\langle\lambda_{t}-\beta_{\varepsilon}(0),\nu_{t}\cdot[\\![v_{t}]\\!]\rangle_{\Sigma_{t}}=\langle\mu_{t},\nu_{t}\cdot[\\![u_{t}]\\!]\rangle_{\Sigma_{t}},$
(6.8)

where $\beta_{\varepsilon}(0)=-exp(-2)$ in (4.3) does not depend on
$\varepsilon$. In case $v_{t}$ is smooth, the following boundary value
relations hold:

$\displaystyle{\rm div}\,\sigma(v_{t})=0$ $\displaystyle\text{ in
}\Omega\setminus\Sigma_{t},$ $\displaystyle v_{t}=0\text{ on }\Gamma^{\rm
D}_{t},\quad\sigma(v_{t})n=u_{t}-z\text{ on }\Gamma^{\rm
O}_{t},\quad\sigma(v_{t})n=0$ $\displaystyle\text{ on }\Gamma^{\rm
N}_{t}\setminus\Gamma^{\rm O}_{t},$
$\displaystyle[\\![\sigma(v_{t})\nu_{t}]\\!]=0,\quad(\sigma(v_{t})\nu_{t})_{\tau_{t}}=\int_{0}^{1}\nabla^{2}\alpha_{\rm
f}([\\![ru_{t}]\\!]_{\tau_{t}})\,[\\![v_{t}]\\!]_{\tau_{t}}\,dr,$
$\displaystyle\nu_{t}\cdot(\sigma(v_{t})\nu_{t})=\int_{0}^{1}\alpha_{\rm
c}^{\prime\prime}(\nu_{t}\cdot[\\![ru_{t}]\\!])\,(\nu_{t}\cdot[\\![v_{t}]\\!])\,dr+\mu_{t}$
$\displaystyle\text{ on }\Sigma_{t}.$ (6.9)

(iii) The quadruple $(u_{t},v_{t},\lambda_{t},\mu_{t})$ constitutes an
accumulation point as $\varepsilon_{k}\to 0$:

$u^{\varepsilon_{k}}_{t}\to u_{t}\text{ strongly in
$H^{1}(\Omega\setminus\Sigma_{t})^{d}$},\quad
v^{\varepsilon_{k}}_{t}\rightharpoonup v_{t}\text{ weakly in
$H^{1}(\Omega\setminus\Sigma_{t})^{d}$},$ (6.10)
$\beta_{\varepsilon_{k}}(\nu_{t}\cdot[\\![u^{\varepsilon_{k}}_{t}]\\!])\to\lambda_{t}\text{
strongly in $H^{1/2}(\Sigma_{t})^{\star}$},$ (6.11)
$\int_{0}^{1}\beta^{\prime}_{\varepsilon_{k}}(\nu_{t}\cdot[\\![ru^{\varepsilon_{k}}_{t}]\\!])\,(\nu_{t}\cdot[\\![v^{\varepsilon_{k}}_{t}]\\!])\,dr\rightharpoonup\mu_{t}\text{
$\star$-weakly in $H^{1/2}(\Sigma_{t})^{\star}$}.$ (6.12)

###### Proof

(i) Taking the limit in (4.4) with the help of the weak convergence
$u^{\varepsilon_{k}}_{t}\rightharpoonup u_{t}$ in (6.2) we get

$\lim_{\varepsilon_{k}\to
0}\int_{\Sigma_{t}}\beta_{\varepsilon_{k}}(\nu_{t}\cdot[\\![u^{\varepsilon_{k}}_{t}]\\!])\,(\nu_{t}\cdot[\\![u]\\!])\,dS_{x}=\int_{\Gamma^{\rm
N}_{t}}g\cdot
u\,dS_{x}-\int_{\Omega\setminus\Sigma_{t}}\sigma(u_{t})\cdot\epsilon(u)\,dx\\\
-\int_{\Sigma_{t}}\bigl{\\{}\nabla\alpha_{\rm
f}([\\![u_{t}]\\!]_{\tau_{t}})\cdot[\\![u]\\!]_{\tau_{t}}+\alpha_{\rm
c}^{\prime}(\nu_{t}\cdot[\\![u_{t}]\\!])\,(\nu_{t}\cdot[\\![u]\\!])\bigr{\\}}\,dS_{x}=:\langle\lambda_{t},\nu_{t}\cdot[\\![u]\\!]\rangle_{\Sigma_{t}}.$
(6.13)

This implies the $\star$-weak convergence
$\beta_{\varepsilon_{k}}(\nu_{t}\cdot[\\![u^{\varepsilon_{k}}_{t}]\\!])\rightharpoonup\lambda_{t}$
in $H^{1/2}(\Sigma_{t})^{\star}$, equation (6.4), and $\lambda_{t}\leq 0$ in
(6.5) due to $\beta_{\varepsilon_{k}}\leq 0$ in Lemma 1. Testing (6.13) with
$u=u^{\varepsilon_{k}}_{t}$ and using (4.2) such that

$\int_{\Sigma_{t}}\beta_{\varepsilon_{k}}(\nu_{t}\cdot[\\![u^{\varepsilon_{k}}_{t}]\\!])\,(\nu_{t}\cdot[\\![u^{\varepsilon_{k}}_{t}]\\!])\,dS_{x}\geq\frac{1}{\varepsilon_{k}}\int_{\Sigma_{t}}(\bigl{[}\nu_{t}\cdot[\\![u^{\varepsilon_{k}}_{t}]\\!]\bigr{]}^{-})^{2}\,dS_{x}-2\varepsilon_{k}K_{\beta}\geq-2\varepsilon_{k}K_{\beta}\to
0,$

after passage $\varepsilon_{k}\to 0$, we get in the limit
$\langle\lambda_{t},\nu_{t}\cdot[\\![u_{t}]\\!]\rangle_{\Sigma_{t}}\geq 0$. On
the other hand we have
$\langle\lambda_{t},\nu_{t}\cdot[\\![u_{t}]\\!]\rangle_{\Sigma_{t}}\leq 0$
because $\lambda_{t}\leq 0$ and the non-penetration
$\nu_{t}\cdot[\\![u_{t}]\\!]\geq 0$, which together lead to the equality in
(6.5). Substituting $\lambda_{t}$ with the expression (6.6) at $\Sigma_{t}$ we
derive the VI (3.14) and its boundary value formulation (3.1). Thus, $u_{t}\in
K(\Omega_{t})$ yields a solution of the cohesive crack problem.

(ii) The limit of the adjoint equation (4.15) using the convergences in (6.2)
is

$\lim_{\varepsilon_{k}\to
0}\int_{\Sigma_{t}}\int_{0}^{1}\beta^{\prime}_{\varepsilon_{k}}(\nu_{t}\cdot[\\![ru^{\varepsilon_{k}}_{t}]\\!])\,(\nu_{t}\cdot[\\![v]\\!])\,(\nu_{t}\cdot[\\![v^{\varepsilon_{k}}_{t}]\\!])\,dr\,dS_{x}\\\
=\int_{\Gamma^{\rm O}_{t}}(u_{t}-z)\cdot
v\,dS_{x}-\int_{\Omega\setminus\Sigma_{t}}\sigma(v)\cdot\epsilon(v_{t})\,dx-\int_{\Sigma_{t}}\int_{0}^{1}\Bigl{\\{}\Bigl{(}\nabla^{2}\alpha_{\rm
f}([\\![ru_{t}]\\!]_{\tau_{t}})\,[\\![v]\\!]_{\tau_{t}}\Bigr{)}\cdot[\\![v_{t}]\\!]_{\tau_{t}}\\\
+\alpha_{\rm
c}^{\prime\prime}(\nu_{t}\cdot[\\![ru_{t}]\\!])\,(\nu_{t}\cdot[\\![v]\\!])(\nu_{t}\cdot[\\![v_{t}]\\!])\Bigr{\\}}\,dS_{x}=:\langle\mu_{t},\nu_{t}\cdot[\\![v]\\!]\rangle_{\Sigma_{t}}.$
(6.14)

The convergence in (6.14) implies (6.12) and the adjoint equation (6.7).
Derivation of the boundary value relations (6.1) is standard. According to (4)
we have

$\langle\beta_{\varepsilon_{k}}(\nu_{t}\cdot[\\![u^{\varepsilon_{k}}_{t}]\\!]),\nu_{t}\cdot[\\![v^{\varepsilon_{k}}_{t}]\\!]\rangle_{\Sigma_{t}}=\langle\int_{0}^{1}\beta^{\prime}_{\varepsilon_{k}}(\nu_{t}\cdot[\\![ru^{\varepsilon_{k}}_{t}]\\!])\,(\nu_{t}\cdot[\\![u^{\varepsilon_{k}}_{t}]\\!])\,dr+\beta_{\epsilon}(0),\nu_{t}\cdot[\\![v^{\varepsilon_{k}}_{t}]\\!]\rangle_{\Sigma_{t}},$

hence based on (6.11) and (6.12) we derive in the limit the compatibility
equation (6.8).

(iii) The weak convergences in (6.10) are proved in Lemma 5. To justify the
strong convergence $u^{\varepsilon_{k}}_{t}-u_{t}\to 0$, we subtract (6.4)
from (4.4), test the difference with $u=u^{\varepsilon_{k}}_{t}-u_{t}$ and
rearrange the terms as follows

$\int_{\Omega\setminus\Sigma_{t}}\sigma(u^{\varepsilon_{k}}_{t}-u_{t})\cdot\epsilon(u^{\varepsilon_{k}}_{t}-u_{t})\,dx=-\int_{\Sigma_{t}}\bigl{\\{}\bigl{(}\nabla\alpha_{\rm
f}([\\![u^{\varepsilon_{k}}_{t}]\\!]_{\tau_{t}})-\nabla\alpha_{\rm
f}([\\![u_{t}]\\!]_{\tau_{t}})\bigr{)}\cdot[\\![u^{\varepsilon_{k}}_{t}-u_{t}]\\!]_{\tau_{t}}\\\
+\bigl{(}[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon_{k}}](\nu_{t}\cdot[\\![u^{\varepsilon_{k}}_{t}]\\!])-[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon_{k}}](\nu_{t}\cdot[\\![u_{t}]\\!])\bigr{)}(\nu_{t}\cdot[\\![u^{\varepsilon_{k}}_{t}-u_{t}]\\!])\bigr{\\}}\,dS_{x}\\\
-\langle\beta_{\varepsilon_{k}}(\nu_{t}\cdot[\\![u_{t}]\\!])-\lambda_{t},\nu_{t}\cdot[\\![u^{\varepsilon_{k}}_{t}-u_{t}]\\!]\rangle_{\Sigma_{t}}.$
(6.15)

Using the monotony of $\beta_{\varepsilon_{k}}$ and the uniform boundedness
$-1<\beta_{\varepsilon_{k}}(0)\leq\beta_{\varepsilon_{k}}(\nu_{t}\cdot[\\![u_{t}]\\!])\leq
0$ for $\nu_{t}\cdot[\\![u_{t}]\\!]\geq 0$, the strong convergence in (6.10)
follows upon taking the limit in (6.15) as $\varepsilon_{k}\to 0$, see (6.2).
Consequently, from (4.4) and (6.13) we conclude the strong convergence in
(6.11). This finishes the proof. ∎

Based on assertion (iii) of Theorem 6.1 we get the following.

######  Corollary 3 (Limit optimization problems)

For the fixed $(\lambda_{t},\mu_{t})\in(H^{1/2}(\Sigma_{t})^{\star})^{2}$ from
Theorem 6.1 and Lagrangian $\mathcal{L}$ from (6.3), the pair
$(u_{t},v_{t})\in V(\Omega_{t})^{2}$ solving optimality conditions (6.4),
(6.5) and (6.7) satisfies the primal problem:

$\mathcal{L}(u_{t},u_{t},v;\Omega_{t})-\langle\lambda_{t},\nu_{t}\cdot[\\![v]\\!]\rangle_{\Sigma_{t}}\leq\mathcal{L}(u_{t},u_{t},v_{t};\Omega_{t})-\langle\lambda_{t},\nu_{t}\cdot[\\![v_{t}]\\!]\rangle_{\Sigma_{t}}$
(6.16)

for all $v\in V(\Omega_{t})$, and the dual problem:

$\mathcal{L}(u_{t},u_{t},v_{t};\Omega_{t})-\langle\mu_{t},\nu_{t}\cdot[\\![u_{t}]\\!]\rangle_{\Sigma_{t}}-\langle\beta_{\epsilon}(0),\nu_{t}\cdot[\\![v_{t}]\\!]\rangle_{\Sigma_{t}}\\\
\leq\mathcal{L}(u_{t},u,v_{t};\Omega_{t})-\langle\mu_{t},\nu_{t}\cdot[\\![u]\\!]\rangle_{\Sigma_{t}}-\langle\beta_{\epsilon}(0),\nu_{t}\cdot[\\![v_{t}]\\!]\rangle_{\Sigma_{t}}\quad\text{for
all }u\in V(\Omega_{t}).$ (6.17)

By the virtue of compatibility (6.8), the corresponding optimal value function
for the objective $\mathcal{J}$ in (2.3) has the equivalent representations
using the adjoint equation as follows:

$\mathcal{J}(u_{t};\Omega_{t})=\mathcal{L}(u_{t},u_{t},v_{t};\Omega_{t})-\langle\lambda_{t},\nu_{t}\cdot[\\![v_{t}]\\!]\rangle_{\Sigma_{t}}\\\
=\mathcal{L}(u_{t},u_{t},v_{t};\Omega_{t})-\langle\mu_{t},\nu_{t}\cdot[\\![u_{t}]\\!]\rangle_{\Sigma_{t}}-\langle\beta_{\epsilon}(0),\nu_{t}\cdot[\\![v_{t}]\\!]\rangle_{\Sigma_{t}}.$
(6.18)

###### Proof

Indeed, taking the limit $\varepsilon_{k}\to 0$ in the saddle-point problem
(4.11) with the Lagrangian $\mathcal{\tilde{L}}^{\varepsilon_{k}}$ from
(4.13), and observing (6.10)–(6.12), the inequalities (6.16), (6.17) follow.
From the $\varepsilon$-dependent representation (4.12) of the optimal value
function $\mathcal{J}$ and by using the compatibility (6.8) we derive the
limit formula (6.18). ∎

We finish by noting the difficulty that, in general, we can pass to the limit
as $\varepsilon\to 0^{+}$ neither in the term
$\int_{0}^{1}\beta^{\prime}_{\varepsilon}(\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])\,dr$
in the Lagrangian $\tilde{\mathcal{L}}^{\varepsilon}$ in (4.10), nor in the
term $\beta^{\prime}_{\varepsilon}(\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!])$
in the shape derivative ${\textstyle\frac{\partial}{\partial
s}}\tilde{\mathcal{L}}^{\varepsilon}$ in (5.20) and (5.26). Otherwise, if

$\eta_{t}=\lim_{\varepsilon\to
0^{+}}\int_{0}^{1}\beta^{\prime}_{\varepsilon}(\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])\,dr$

exists, then the compatibility properties
$\lambda_{t}=(\nu_{t}\cdot[\\![u_{t}]\\!])\eta_{t}+\beta_{\varepsilon}(0)$ and
$\mu_{t}=(\nu_{t}\cdot[\\![v_{t}]\\!])\eta_{t}$ which are stronger than (6.8)
hold. For a factorization of $\lambda_{t}$ and $\mu_{t}$, additional solution
regularity, as in the particular case of obstacle problems, could be helpful,
see Bar/84 ; HK/09 ; MP/84 .

##  7 Shape optimization of breaking line

We apply the theoretical results to a numerical example in 2D.

As a true shape to be identified within an admissibility set
$\\{\Sigma_{t}\\}$ we take the piecewise-linear line

$\Sigma=\\{x_{1}\in(0,1),\,x_{2}=\psi(x_{1})\\},\quad\psi(x_{1})=\min(0.3,x_{1}/3+0.1),$

which breaks the rectangle $\Omega=(0,1)\times(0,0.5)$ into two parts
$\Omega^{\pm}$. Let the boundary $\partial\Omega$ be split symmetrically into
the fixed Dirichlet part $\Gamma^{\rm
D}=\\{x_{1}\in\\{0,1\\},\,x_{2}\in(0,0.5)\\}$ and the Neumann part
$\Gamma^{\rm N}=\\{x_{1}\in(0,1),\,x_{2}\in\\{0,0.5\\}\\}$. For an isotropic
elastic body occupying $\Omega$ we set the material parameters: Young modulus
$E_{\rm Y}=73000$ (mPa), Poisson ratio $\nu_{\rm P}=0.34$, and the
corresponding Lamé parameters $\mu_{\rm L}=E_{\rm Y}/(2(1+\nu_{\rm P}))$,
$\lambda_{\rm L}=2\mu_{\rm L}\nu_{\rm P}/(1-2\nu_{\rm P})$. For the matrix $C$
of isotropic elastic coefficients the stress-strain relations are

$\sigma_{ij}=2\mu_{\rm L}\epsilon_{ij}+\lambda_{\rm
L}(\epsilon_{11}+\epsilon_{22})\delta_{ij},\quad i,j=1,2.$

We rely on the approximation of $\nu_{t}\cdot[\\![u]\\!]$ by
$[\\![u]\\!]_{2}:=[\\![u_{2}]\\!]$, and
$[\\![u]\\!]_{\tau_{t}}=[\\![u]\\!]_{1}\tau_{t}$ with
$[\\![u]\\!]_{1}:=[\\![u_{1}]\\!]$ at $\Sigma_{t}$, which is reasonable for
flat shapes. For a friction function in one variable $\alpha_{\rm f}(s)=F_{\bf
b}\sqrt{\delta^{2}+s^{2}}$ such that $\nabla\alpha_{\rm f}=\tau_{t}\alpha_{\rm
f}^{\prime}$, and $\alpha_{\rm c}(s)=K_{\rm c}s/(\kappa+|s|)$, applying to the
body the traction force

$g_{1}=0,\quad g_{2}(x)=(1-7x_{1}/4)(4x_{2}-1)\mu_{\rm L},$

according to Theorem 3.1 there exists a solution $z\in
H^{1}(\Omega\setminus\Sigma)^{2}$ such that $z=0$ on $\Gamma^{\rm D}$,
$[\\![z]\\!]_{2}\geq 0\text{ on }\Sigma$, and satisfying the VI (3.14):

$\int_{\Omega\setminus\Sigma}\sigma(z)\cdot\epsilon(u-z)\,dx+\int_{\Sigma}\bigl{\\{}\alpha_{\rm
f}^{\prime}([\\![z]\\!]_{1})\,[\\![u-z]\\!]_{1}\\\ +\alpha_{\rm
c}^{\prime}([\\![z]\\!]_{2})\,[\\![u-z]\\!]_{2}\bigr{\\}}\,dS_{x}\geq\int_{\Gamma^{\rm
N}}g\cdot(u-z)\,dS_{x}$ (7.1)

for all test functions $u\in H^{1}(\Omega\setminus\Sigma)^{2}$ such that $u=0$
on $\Gamma^{\rm D}$ and $[\\![u]\\!]_{2}\geq 0\text{ on }\Sigma$. Let the
observation boundary be $\Gamma^{\rm O}=\Gamma^{\rm N}$. We insert the
solution $z$ of (7.1) as a measurement into the objective function
$\mathcal{J}$ in (2.3) and consider the shape optimization problem: find
$\Sigma_{t}$ from the feasible set
$\mathfrak{S}=\\{x\in\Omega:\,x_{1}\in(0,1),x_{2}=\psi(x_{1})\in(0,0.5),\,\psi\in
C^{0,1}(0,1)\\}$ such that

$\min_{\Sigma_{t}\in\mathfrak{S}}\mathcal{J}(u_{t};\Omega_{t})=\frac{1}{2}\int_{\Gamma^{\rm
O}_{t}}|u_{t}-z|^{2}\,dS_{x}+\rho|\Sigma_{t}|,\quad\text{where $u_{t}$ solves
\eqref{2.13}}.$ (7.2)

Evidently, the trivial minimum in (7.2) is attained as $\Sigma_{t}=\Sigma$ and
$u_{t}=z$. To avoid the inverse crime, we use two different meshes for $z$,
and for $u_{t}$ when solving the inverse problem.

Now we discretize the problem. For fixed $t$, let $\Omega^{1}_{t,h}$,
$\Omega^{2}_{t,h}$ be triangular meshes with grid size $h>0$ in
$\Omega^{1}_{t}$, $\Omega^{2}_{t}$, which are compatible at the interface such
that
$\Sigma_{t,h}:=\Sigma_{t}\cap\partial\Omega^{1}_{t,h}=\Sigma_{t}\cap\partial\Omega^{2}_{t,h}$.
At the interface $\Sigma_{t,h}$ the nonlinear functions are set: friction
$\alpha_{\rm f}$ from (3.4) with $F_{\rm b}=10^{-5}$ (mPa); cohesion
$\alpha_{\rm c}$ from (3.6) with $m=1$, $K_{\rm c}=10^{-3}$ (mPa$\cdot$m),
$\kappa=10^{-2}$ (m). The parameters $\delta$, $h$ are assumed sufficiently
small such that we rely on the discretization:

$(\alpha_{\rm f})_{h}(s)=F_{\rm b}|s|,\quad(\alpha_{\rm
f}^{\prime})_{h}(s)=F_{\rm b}\,{\rm sgn}(s),\quad(\alpha_{\rm
f}^{\prime\prime})_{h}(s)=0;\\\ (\alpha_{\rm c})_{h}(s)=\frac{K_{\rm
c}}{\kappa}\min(\kappa,|s|),\quad(\alpha_{\rm c}^{\prime})_{h}(s)=\frac{K_{\rm
c}}{\kappa}\,{\rm ind}\\{|s|<\kappa\\},\quad(\alpha_{\rm
c}^{\prime\prime})_{h}(s)=0.$ (7.3)

Figure 5: Computed true solution $z_{h}$ to (7.1) within current configuration
(a); componentwise in (b), (c).

After FE-discretization of problem (7.1) according to (7.3) on a grid of size
$h=10^{-2}$, we solve it by a primal-dual active set (PDAS) iterative
algorithm developed in HKK/11 . The reference numerical solution $z_{h}$
obtained after 4 iterations is plotted in Figure 5. In plot (a) we present the
grid in the so-called current or deformed configuration $x+z(x)$ for
$x\in\Omega\setminus\Sigma$ under the traction force $g$ prescribed at
$\Gamma^{\rm N}$. Here we observe an open part of $\Sigma$ which is the
complement to the cohesion part (where $[\\![z]\\!]_{2}<\kappa$) with contact
(where $[\\![z]\\!]_{2}=0$) marked by colors in finite elements adjacent to
the interface. In plots (b), (c) of Figure 5 the solution components
$(z_{h})_{1}$, $(z_{h})_{2}$ in the reference configuration
$\Omega\setminus\Sigma$ are depicted.

According to Theorem 4.1 we approximate the VI (3.14) by the
$\varepsilon$-regularized cohesive crack problem (4.4). For sufficiently small
$\varepsilon$ fixed, the compliance $\beta_{\varepsilon}$ from (4.3) is
discretized as

$(\beta_{\varepsilon})_{h}(s)=\frac{1}{\varepsilon}\min(0,s),\quad(\beta_{\varepsilon}^{\prime})_{h}(s)=\frac{1}{\varepsilon}\,{\rm
ind}\\{s<0\\}.$ (7.4)

Let $V_{t,h}$ be the finite element (FE) space of piecewise-linear functions
such that

$V_{t,h}\subset V(\Omega_{t,h})=\\{u\in H^{1}(\Omega^{+}_{t,h})^{2}\cap
H^{1}(\Omega^{-}_{t,h})^{2}|\quad u=0\text{ on }\Gamma^{\rm D}\\}.$

Then the discretization of the penalty equation (4.4) becomes: find
$u^{\varepsilon}_{t,h}\in V_{t,h}$ such that

$\int_{\Omega\setminus\Sigma_{t,h}}\sigma(u^{\varepsilon}_{t,h})\cdot\epsilon(u_{h})\,dx+\int_{\Sigma_{t,h}}\bigl{\\{}(\alpha_{\rm
f}^{\prime})_{h}\bigl{(}[\\![u^{\varepsilon}_{t,h}]\\!]_{1}\bigr{)}\cdot[\\![u_{h}]\\!]_{1}\\\
+[(\alpha_{\rm
c}^{\prime})_{h}+(\beta_{\varepsilon})_{h}]\bigl{(}[\\![u^{\varepsilon}_{t,h}]\\!]_{2}\bigr{)}\,[\\![u_{h}]\\!]_{2}\bigr{\\}}\,dS_{x}=\int_{\Gamma^{\rm
N}}g\cdot u_{h}\,dS_{x},$ (7.5)

and due to (7.3) the discrete adjoint equation (4.15) reads: find
$v^{\varepsilon}_{t,h}\in V_{t,h}$ such that

$\int_{\Omega\setminus\Sigma_{t,h}}\sigma(v_{h})\cdot\epsilon(v^{\varepsilon}_{t,h})\,dx+\int_{\Sigma_{t,h}}\int_{0}^{1}(\beta^{\prime}_{\varepsilon})_{h}([\\![ru^{\varepsilon}_{t,h}]\\!]_{2})\,[\\![v_{h}]\\!]_{2}[\\![v^{\varepsilon}_{t,h}]\\!]_{2}\,dr\,dS_{x}\\\
=\int_{\Gamma^{\rm N}}(u^{\varepsilon}_{t,h}-z_{h})\cdot v_{h}\,dS_{x}$ (7.6)

for all test functions $u_{h},v_{h}\in V_{t,h}$.

After solving problems (7.5) and (7.6), since $\Gamma^{\rm D}$ and
$\Gamma^{\rm N}=\Gamma^{\rm O}$ are fixed in this example, according to
Corollary 2 we calculate $\mathcal{D}^{\varepsilon}_{3}$ at the moving
boundary $\Sigma_{t,h}$, and $\mathcal{D}_{1}$ at $\Sigma_{t,h}\cap\Gamma^{\rm
D}$:

$(\mathcal{D}_{1})_{t,h}=[\\![\nabla(u^{\varepsilon}_{t,h})^{\top}\sigma(v^{\varepsilon}_{t,h})+\nabla(v^{\varepsilon}_{t,h})^{\top}\sigma(u^{\varepsilon}_{t,h})]\\!]\tau_{t}(2x_{1}-1),\\\
(\mathcal{D}^{\varepsilon}_{3})_{t,h}=[\\![\sigma(u^{\varepsilon}_{t,h})\cdot\epsilon(v^{\varepsilon}_{t,h})]\\!]+\varkappa_{t}\bigl{(}\rho-(p_{\rm
f})_{t,h}-(p^{\varepsilon}_{\rm c})_{t,h}\bigr{)}-\nu_{t}\cdot\bigl{(}(\nabla
p_{\rm f})_{t,h}+(\nabla p^{\varepsilon}_{\rm c})_{t,h}\bigr{)},$ (7.7)

where $\rho=1/\mu_{\rm L}$ is set, $(q_{\rm f})_{t,h}=(q^{\varepsilon}_{\rm
c})_{t,h}=0$ by the virtue of (7.3), (7.4). Relying on a flat shape
approximation we take $\nabla\nu_{t}=\nabla\tau_{t}=0$ and

$(p_{\rm f})_{t,h}=(\alpha_{\rm
f}^{\prime})_{h}([\\![u^{\varepsilon}_{t,h}]\\!]_{1})\,[\\![v^{\varepsilon}_{t,h}]\\!]_{1},\quad(p^{\varepsilon}_{\rm
c})_{t,h}=[(\alpha_{\rm
c}^{\prime})_{h}+(\beta_{\varepsilon})_{h}]([\\![u^{\varepsilon}_{t,h}]\\!]_{2})\,[\\![v^{\varepsilon}_{t,h}]\\!]_{2},\\\
(\nabla p_{\rm f})_{t,h}=[\\![\nabla
v^{\varepsilon}_{t,h}]\\!]^{\top}\tau_{t}\,(\alpha_{\rm
f}^{\prime})_{h}([\\![u^{\varepsilon}_{t,h}]\\!]_{1}),\\\ (\nabla
p^{\varepsilon}_{\rm c})_{t,h}=[\\![\nabla
v^{\varepsilon}_{t,h}]\\!]^{\top}\nu_{t}\,[(\alpha_{\rm
c}^{\prime})_{h}+(\beta_{\varepsilon})_{h}]([\\![u^{\varepsilon}_{t,h}]\\!]_{2})+[\\![\nabla
u^{\varepsilon}_{t,h}]\\!]^{\top}\nu_{t}\,(\beta_{\varepsilon}^{\prime})_{h}([\\![u^{\varepsilon}_{t,h}]\\!]_{2})\,[\\![v^{\varepsilon}_{t,h}]\\!]_{2}.$
(7.8)

The discrete velocity $\Lambda_{H}$ at interface $\Sigma_{t}$ is defined on a
coarse grid of size $H>0$. According to Corollary 3 we get a descent direction
by setting $(\Lambda_{H})_{1}=0$ and

$(\Lambda_{H})_{2}=\frac{k}{\sqrt{h}}(2x_{1}-1)\nu_{t}\cdot(\mathcal{D}_{1})_{t,h}\text{
at $\Sigma_{t,h}\cap\Gamma^{\rm
D}$},\quad(\Lambda_{H})_{2}=-k(\mathcal{D}^{\varepsilon}_{3})_{t,h}\text{ at
$\Sigma_{t,h}\setminus\Gamma^{\rm D}$},$ (7.9)

where the scaling $k=0.1h/\|(\Lambda_{H})_{2}\|_{C(\overline{\Sigma_{t,h}})}$
is chosen, and the weight $1/\sqrt{h}$ at $\Gamma^{\rm D}$ was found
empirically as in GKK/20 . Based on formulas (7.7)–(7.9) we formulate the
shape optimization algorithm of breaking line identification for the
discretized version of (4.8) .

######  Algorithm 1 (breaking line identification)

  * (0)

Initialize the constant grid function $\psi^{(0)}_{H}=0.25$ at points
$s_{H}\in[0,1]$. Determine the line segment
$\Sigma^{(0)}=\\{x_{1}\in(0,1),\,x_{2}=\psi^{(0)}(x_{1})\\}$, where
$\psi^{(0)}$ is the linear interpolate of $\psi^{(0)}_{H}$; set $n=0$.

  * (1)

Set the interface $\Sigma_{t,h}=\Sigma^{(n)}$ and construct triangulations
$\Omega^{1}_{t,h}$, $\Omega^{2}_{t,h}$; find solutions
$u^{\varepsilon}_{t,h}$, $v^{\varepsilon}_{t,h}$ of the discrete penalty and
adjoint equations (7.5), (7.6).

  * (2)

Calculate a velocity $(\Lambda_{H})_{2}$ by formula (7.9); update the grid
function

$\psi^{(n+1)}_{H}=\psi^{(n)}_{H}+(\Lambda_{H})_{2}\quad\text{at points
$s_{H}\in[0,1]$}.$ (7.10)

From linear interpolation $\psi^{(n+1)}$ of $\psi^{(n+1)}_{H}$ determine the
piecewise-linear segment

$\Sigma^{(n+1)}=\\{x_{1}\in(0,1),\,x_{2}=\psi^{(n+1)}(x_{1})\\}.$ (7.11)

  * (3)

If stopping criterion holds, then STOP; else set $n=n+1$ and go to Step (1).

For 11 equidistant points $s_{H}$ as $H=0.1$, the numerical result of
Algorithm 1 after $\\#n=200$ iterations (the stopping criterion) is depicted
in Figure 6.

Figure 6: Iterations of $\Sigma^{(n)}$ (a); objective function ratio
$\mathcal{J}^{(n)}/\mathcal{J}^{(0)}$ (b); shape error ratio (c).

In plot (a) the selected iterations $n=0,10,20,40,100,200$ of $\Sigma^{(n)}$
from (7.11) are drawn in $\Omega$ in comparison with the true interface
$\Sigma$ (the thick solid line). In plot (b) of Figure 6 we plot the ratio
$\mathcal{J}^{(n)}/\mathcal{J}^{(0)}$ of the objective function during
iterations of $\Sigma_{t,h}=\Sigma^{(n)}$, where we recall

$\mathcal{J}^{(n)}(u^{\varepsilon}_{t,h};\Omega\setminus\Sigma^{(n)})=\frac{1}{2}\int_{\Gamma^{\rm
O}}|u^{\varepsilon}_{t,h}-z_{h}|^{2}\,dS_{x}+\rho|\Sigma^{(n)}|\quad\text{subject
to }\eqref{6.5}.$ (7.12)

The computed ratio attains as minimum $0,6\%$. In plot (c) of Figure 6 the
ratio of shape error $\|\Sigma^{(n)}-\Sigma\|/\|\Sigma^{(0)}-\Sigma\|$ is
plotted versus $n\in[0,200]$, where according to (7.10)

$\|\Sigma^{(n)}-\Sigma\|:=\|\psi^{(n)}-\psi\|_{C([0,1])}.$ (7.13)

Here the accuracy of shape identification attains only $46\%$. It is worth
noting that the computation is presented for small penalty parameter
$\varepsilon=10^{-8}$, while insufficiently small value $\varepsilon=10^{-5}$
causes some increase of the ratio curves after reaching the minimum; see
Figure 6 (b), (c).

From the simulation we conclude the following. In Figure 6 (a) it can be
observed that the left part of curve $\Sigma$, where the constraints are
inactive (see Figure 5 (a)), is recovered well by the identification Algorithm
1, whereas the right part of interface, where either contact or cohesion
occurs, the initialized $\Sigma^{(0)}$ is almost not modified during the
iteration.

To remedy the hidden part, we apply to the same physical and geometrical
configuration the traction force $g_{2}(x)=(1-5x_{1}/4)(4x_{2}-1)\mu_{\rm L}$,
which is more stretching than the one from Figure 5 (a). Because of that, the
whole $\Sigma$ is open, neither contact nor cohesion occur at the interface
(see Figure 7 (a)). The corresponding result of Algorithm 1 is depicted in
Figure 7.

Figure 7: The true solution $z_{h}$ (a); iterations of $\Sigma^{(n)}$ (b);
objective function ratio and shape error ratio (c).

Here plot (b) presents the selected iterations of $\Sigma^{(n)}$, and plot (c)
shows the objective function ratio $\mathcal{J}^{(n)}/\mathcal{J}^{(0)}$
together with the shape error ratio
$\|\Sigma^{(n)}-\Sigma\|/\|\Sigma^{(0)}-\Sigma\|$. The former ratio attains
the minimum $0,25\%$, and the latter one $23\%$ of accuracy. Now we see in
Figure 7 (b) that the whole curve $\Sigma$ is recovered well compared to that
from Figure 6 (a).

##  8 Conclusions

The Barenblatt’s crack model assuming cohesion at a breaking line is stated as
the variational inequality due to the non-penetration condition and penalized
using smooth Lavrentiev’s approximation. For the geometry-dependent least-
square function describing misfit of the solution from a boundary measurement,
the expression of shape derivative is derived in an analytical form. On its
basis, from our numerical simulation we make a conclusion that the suggested
breaking line identification algorithm is consistent within the setup of
destructive physical analysis (DPA).

### Data availability statement

Data sharing not applicable to this article as no datasets were generated or
analysed during the current study.

###### Acknowledgements.

The research was supported by the ERC advanced grant 668998 (OCLOC) under the
EU’s H2020 research program.

## References

  * (1) Alekseev, G.V., Tereshko, D.A., Shestopalov, Yu.V.: Optimization approach for axisymmetric electric field cloaking and shielding. Inverse Probl. Sci. Eng. 49, 40–55 (2021)
  * (2) Andersson, L.-E.: A quasistatic frictional problem with a normal compliance penalization term. Nonlinear Anal. 37, 689–705 (1999)
  * (3) Barbu, V.: Optimal Control of Variational Inequalities, Pitman, London (1984)
  * (4) Barenblatt, G.I.: The mathematical theory of equilibrium cracks in brittle fracture. Adv. Appl. Mech. 7, 55–129 (1962)
  * (5) Bartle, R.G., Joichi, J.T.: The preservation of convergence of measurable functions under composition. Proc. Amer. Math. Soc. 12 122–126 (1961)
  * (6) Bellis, C., Bonnet, M.: Qualitative identification of cracks using 3D transient elastodynamic topological derivative: formulation and FE implementation. Comput. Meth. Appl. Mech. Eng. 253, 89–105 (2013)
  * (7) Bergounioux, M.: Use of augmented Lagrangian methods for the optimal control of obstacle problems. J. Optim. Theory Appl. 95, 101–126 (1997)
  * (8) Bratov, V.A., Morozov, N.F., Petrov ,Yu.V.: Dynamic Strength of Continuum, St.Petersburg University (2009)
  * (9) Bredies, K., Lorenz, D.A., Reiterer, S.: Minimization of non-smooth, non-convex functionals by iterative thresholding. J. Optim. Theory Appl. 165, 78–112 (2015)
  * (10) Casas, E., Clason, C., Kunisch, K.: Parabolic control problems in measure spaces with sparse solutions. SIAM J. Control Optim. 51, 28–63 (2013)
  * (11) Charlotte, M., Francfort, G., Marigo, J.-J., Truskinovsky, L.: Revisiting brittle fracture as an energy minimization problem: comparison of Griffith and Barenblatt surface energy models. In: Benallal, A. (ed.): Continuous Damage and Fracture, pp. 1–12. Elsevier, Paris (2000)
  * (12) Correa, R., Seeger, A.: Directional derivative of a minimax function. Nonlinear Anal. Theory Methods Appl. 9, 834–862 (1985)
  * (13) Delfour, M.C., Zolésio, J.-P.: Shape and Geometries: Metrics, Analysis, Differential Calculus, and Optimization, SIAM, Philadelphia (2011)
  * (14) Ekeland, I., Temam, R.: Convex Analysis and Variational Problems, North-Holland, Amsterdam (1976)
  * (15) Estep, D., Lee, S.: Adaptive error control during gradient search for an elliptic optimization problem. Appl. Anal. 92, 1434–1448 (2013)
  * (16) Franců, J.: Weakly continuous operators. Applications to differential equations. Appl. Math. 39, 45–56 (1994)
  * (17) Fremiot, G., Horn, W., Laurain, A., Rao, M., Sokolowski, J.: On the Analysis of Boundary Value Problems in Nonsmooth Domains, Dissertationes Mathematicae 462, Inst. Math. Polish Acad. Sci., Warsaw (2009)
  * (18) Führ, B., Schulz, V., Welker, K.: Shape optimization for interface identification with obstacle problems. Vietnam J. Math. 46, 967–985 (2018)
  * (19) Furtsev, A.I., Itou, H., Rudoy, E.M.: Modeling of bonded elastic structures by a variational method: Theoretical analysis and numerical simulation. Int. J. Solids Struct. 182-183, 100–111 (2020)
  * (20) Ghilli, D., Kunisch, K., Kovtunenko, V.A.: Inverse problem of breaking line identification by shape optimization. J. Inverse Ill-posed Probl. 28, 119–135 (2020)
  * (21) González Granada, J.R., Kovtunenko, V.A., A shape derivative for optimal control of the nonlinear Brinkman–Forchheimer equation. J. Appl. Numer. Optim. 3, 243–261 (2021)
  * (22) Gwinner, J., Jadamba, B., Khan, A.A., Sama, M.: Identification in variational and quasi-variational inequalities. J. Convex Anal. 25, 545–569 (2018)
  * (23) Haslinger, J., Kozubek, T., Kunisch, K., G. Peichl, G.: Shape optimization and fictitious domain approach for solving free boundary problems of Bernoulli type. Comput. Optim. Appl. 26, 231–251 (2003)
  * (24) Hauptmann, A., Ikehata, M., Itou, H., Siltanen, S.: Revealing cracks inside conductive bodies by electric surface measurements. Inverse Probl. 35, 025004 (2019)
  * (25) Hintermüller, M., Hoppe, R.H.W., Löbhard, C.: Use of augmented Lagrangian methods for the optimal control of obstacle problems. ESAIM: COCV 20, 524–546 (2014)
  * (26) Hintermüller, M., Kovtunenko, V.A., Kunisch, K.: Obstacle problems with cohesion: A hemi-variational inequality approach and its efficient numerical solution. SIAM J. Optim. 21, 491–516 (2011)
  * (27) Hintermüller, M., Kopacka, J.: Mathematical programs with complementarity constraints in function space: C- and strong stationarity and a path-following algorithm. SIAM J. Control Optim. 20, 868–902 (2009)
  * (28) Hofmann, B., Kaltenbacher, B., Resmerita, E.: Lavrentiev’s regularization method in Hilbert spaces revisited. Inverse Probl. Imag. 10, 741–764 (2016)
  * (29) Ibragimov, N.H.: Integrating factors, adjoint equations and Lagrangians. J. Math. Anal. Appl. 318, 742–757 (2006)
  * (30) Ito, K., Kunisch, K.: Semi-smooth Newton methods for state-constrained optimal control problems. Systems Control Lett. 50, 221–228 (2003)
  * (31) Ito, K., Kunisch, K.: Lagrange Multiplier Approach to Variational Problems and Applications, SIAM, Philadelphia (2008)
  * (32) Khludnev, A.M., Kovtunenko, V.A.: Analysis of Cracks in Solids, WIT-Press, Southampton, Boston (2000)
  * (33) Khludnev, A.M., Shcherbakov, V.V.: A note on crack propagation paths inside elastic bodies. Appl. Math. Lett. 79, 80–84 (2018)
  * (34) Khludnev, A.M., Sokolowski, J.: Modelling and Control in Solid Mechanics, Birkhäuser, Basel (1997)
  * (35) Kitamura, K.: Crack surface energy: temperature and force dependence. Materials Trans. 49, 643–649 (2008)
  * (36) Knees, D., Mielke, A., Zanini, C.: On the inviscid limit of a model for crack propagation. Math. Models Meth. Appl. Sci. 18, 1529–1569 (2008)
  * (37) Kovtunenko, V.A.: Nonconvex problem for crack with nonpenetration. Z. Angew. Math. Mech. 85, 242–251 (2005)
  * (38) Kovtunenko, V.A.: Primal-dual methods of shape sensitivity analysis for curvilinear cracks with non-penetration. IMA J. Appl. Math. 71, 635–657 (2006)
  * (39) Kovtunenko, V.A.: A hemivariational inequality in crack problems. Optimization 60, 1071–1089 (2011)
  * (40) Kovtunenko, V.A., Kunisch, K.: Problem of crack perturbation based on level sets and velocities. Z. Angew. Math. Mech. 87, 809–830 (2007)
  * (41) Kovtunenko, V.A., Ohtsuka, K.: Shape differentiability of Lagrangians and application to overdetermined problems. In: Itou, H., Hirano, S., Kimura, M., Kovtunenko, V.A., Khludnev, A.M. (eds.): Mathematical Analysis of Continuum Mechanics and Industrial Applications III (Proc. CoMFoS18), Ser. Mathematics for Industry 34, pp. 97–110, Springer, Singapur (2020)
  * (42) Kovtunenko, V.A., Sukhorukov, I.V. : Optimization formulation of the evolutionary problem of crack propagation under quasibrittle fracture. Appl. Mech. Tech. Phys. 47, 704–713 (2006)
  * (43) Kunisch, K., Vexler, B.: Constrained Dirichlet boundary control in $L^{2}$ for a class of evolution equations. SIAM J. Control Optim. 46, 1726–1753 (2007)
  * (44) Lavrentiev, M.M.: Some Improperly Posed Problems of Mathematical Physics, Springer, Berlin, Heidelberg (1967)
  * (45) Lazarev, N., Itou, H.: Optimal location of a rigid inclusion in equilibrium problems for inhomogeneous Kirchhoff–Love plates with a crack. Math. Mech. Solids 24, 3743–3752 (2019)
  * (46) Leugering, G., Prechtel, M., Steinmann, P., Stingl, M.: A cohesive crack propagation model: mathematical theory and numerical solution. Commun. Pure Appl. Anal. 12, 1705–1729 (2013)
  * (47) Leugering, G., Sokolowski, J., Zockowski, A.: Shape- and topology optimization for passive control of crack propagation. In: Pratelli, A., Leugering, G. (eds): New Trends in Shape Optimization. Int. Ser. Numer. Math. 166, pp. 141–197, Birkhäluser, Cham (2015)
  * (48) Luft, D., Schulz, V., Welker, K.: Efficient techniques for shape optimization with variational inequalities using adjoints. SIAM J. Optim. 30, 1922–1953 (2020)
  * (49) Marchuk, G.I., Agoshkov, V.I., Shutyaev, V.P.: Adjoint Equations and Perturbation Algorithms in Nonlinear Problems, CRC Press, Boca Raton (1996)
  * (50) Maurer, H., Zowe, J.: First and second-order necessary and sufficient optimality conditions for infinite-dimensional programming problems. Math. Program. 16, 98–110 (1979)
  * (51) Mignot, F., Puel, J.P.: Optimal control in some variational inequalities. SIAM J. Control Optim. 22, 466–476 (1984)
  * (52) Neitzel, I., Tröltzsch, F.: On regularization methods for the numerical solution of parabolic control problems with pointwise state constraints. ESAIM: COCV 15, 426–453 (2009)
  * (53) Ovcharova, N., Gwinner, J.: A study of regularization techniques of nondifferentiable optimization in view of application to hemivariational inequalities. J. Optim. Theory Appl. 162, 754–778 (2014)
  * (54) Rudoy, E.M.: Differentiation of energy functionals in two-dimensional elasticity theory for solids with curvilinear cracks. J. Appl. Mech. Techn. Phys. 54, 843–852 (2004)
  * (55) Shcherbakov, V.V.: Shape derivatives of energy and regularity of minimizers for shallow elastic shells with cohesive cracks. Nonlinear Anal. Real World Appl. 65, 103505 (2022)
  * (56) Sofonea, M., Matei, A.: Mathematical Models in Contact Mechanics, Cambridge Univ. Press (1991)
  * (57) Sokolowski, J., Zolesio, J.-P.: Introduction to Shape Optimization: Shape Sensitivity Analysis, Springer, Berlin, Heidelberg (1992)
  * (58) Zeng, S.D., Migórski, S., Khan, A.A.: Nonlinear quasi-hemivariational inequalities: existence and optimal control. SIAM J. Control Optim. 59, 1246–1274 (2021)
  * (59) Zowe, J., Kurcyusz, S.: Regularity and stability for the mathematical programming problem in Banach spaces. Appl. Math. Optim. 5, 49–62 (1979)

##  Appendix A Proof of Lemma 2

As $s\to 0$, the following asymptotic expansion of terms in (5.12)–(5.10)
holds (see e.g. (SZ/92, , Chapter 2)):

$z\circ\phi_{s}=z+s\nabla z\Lambda+{\rm o}(s),\quad g\circ\phi_{s}=g+s\nabla
g\Lambda+{\rm o}(s),\quad C\circ\phi_{s}=g+s\nabla C\Lambda+{\rm o}(s),\\\
\nabla\phi_{s}^{-1}\circ\phi_{s}=I-s\nabla\Lambda+{\rm o}(s),\quad
E(\nabla\phi_{s}^{-1}\circ\phi_{s},\tilde{u})=\epsilon(\tilde{u})-sE(\nabla\Lambda,\tilde{u})+{\rm
o}(s),\\\ \omega^{\rm d}_{s}=1+s\,{\rm div}\Lambda+{\rm o}(s),\quad\omega^{\rm
b}_{s}=1+s\,{\rm div}_{\tau_{t}}\Lambda+{\rm o}(s)\\\
\nu_{t+s}\circ\phi_{s}=\nu_{t}+s\nabla\nu_{t}\Lambda+{\rm
o}(s),\quad[\\![\tilde{u}]\\!]_{\tilde{\tau}_{t+s}}=[\\![\tilde{u}]\\!]_{\tau_{t}}+s[\\![\tilde{u}]\\!]_{\nabla\tau_{t}\Lambda}+{\rm
o}(s)$ (A.1)

for $\tilde{u}\in V(\Omega_{t})$. It is worth noting that
$\nabla\nu_{t}\Lambda$ and $\nabla\tau_{t}\Lambda$ from (5.17) are just a
notation used for short, which does not require existence of the gradients
here. The tangential divergence ${\rm div}_{\tau_{t}}\Lambda$ is defined in
(5.18).

Inserting representations (A.1) into the objective
$\tilde{\mathcal{J}}(s,\tilde{u};\Omega_{t})$ and the perturbed Lagrangian
$\tilde{\mathcal{L}}^{\varepsilon}(s,u^{\varepsilon}_{t},\tilde{u},\tilde{v};\Omega_{t})$
given by (5.9), (5.11), we derive their expansions (5.14), (5.15) with respect
to $s$. The asymptotic term ${\textstyle\frac{\partial}{\partial
s}}\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},\tilde{u},\tilde{v};\Omega_{t})$
is from (5.16) at $\tau=0$ (implying that $\Lambda|_{t}=\Lambda$). Since
$\Lambda|_{t+\tau}$ and $\nabla\Lambda|_{t+\tau}$ are continuous functions of
the argument $t+\tau$, the partial derivative
$\tau\mapsto{\textstyle\frac{\partial}{\partial
s}}\tilde{\mathcal{L}}^{\varepsilon}(\tau,\,\cdot\,)$ in (5.16) is continuous.
This finishes the proof.

##  Appendix B Proof of Lemma 3

The first inequality in (5.8) implies the optimality condition
$\partial_{v}\tilde{\mathcal{L}}^{\varepsilon}(s,u^{\varepsilon}_{t},\tilde{u}^{\varepsilon}_{t+s},\tilde{v}^{\varepsilon}_{t+s};\Omega_{t})=0$,
that is

$\int_{\Omega\setminus\Sigma_{t}}\bigl{(}(C\circ\phi_{s})E(\nabla\phi_{s}^{-1}\circ\phi_{s},\tilde{u}^{\varepsilon}_{t+s})\cdot
E(\nabla\phi_{s}^{-1}\circ\phi_{s},\tilde{v})\bigr{)}\,\omega^{\rm d}_{s}dx\\\
+\int_{\Sigma_{t}}\Bigl{\\{}\Bigl{(}\int_{0}^{1}\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})[\\![\tilde{u}^{\varepsilon}_{t+s}]\\!]_{\tilde{\tau}_{t+s}}\,dr+\nabla\alpha_{\rm
f}(0)\Bigr{)}\cdot[\\![\tilde{v}]\\!]_{\tilde{\tau}_{t+s}}+\Bigl{(}\int_{0}^{1}[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])(\tilde{\nu}_{t+s}\cdot[\\![\tilde{u}^{\varepsilon}_{t+s}]\\!])\,dr\\\
+[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](0)\Bigr{)}(\tilde{\nu}_{t+s}\cdot[\\![\tilde{v}]\\!])\Bigr{\\}}\,\omega^{\rm
b}_{s}dS_{x}=\int_{\Gamma^{\rm
N}_{t}}(g\circ\phi_{s})\cdot\tilde{v}\,\omega^{\rm b}_{s}dS_{x}\quad\text{for
all }\tilde{v}\in V(\Omega_{t}).$ (B.1)

According to the asymptotic representation (5.15) and the mean value theorem,
using the operator $A_{\varepsilon}$ from (4.15) it is possible to express the
equation (B.1) in the form

$\langle
A_{\varepsilon}(u^{\varepsilon}_{t})\tilde{u}^{\varepsilon}_{t+s},\tilde{v}\rangle+\int_{\Sigma_{t}}\bigl{(}\nabla\alpha_{\rm
f}(0)\cdot[\\![\tilde{v}]\\!]_{\tau_{t}}+[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](0)(\nu_{t}\cdot[\\![\tilde{v}]\\!])\bigr{)}\,dS_{x}\\\
=\int_{\Gamma^{\rm
N}_{t}}g\cdot\tilde{v}\,dS_{x}+sR_{v}(\alpha^{v}_{s},\tilde{u}^{\varepsilon}_{t+s},\tilde{v})\quad\text{for
all $\tilde{v}\in V(\Omega_{t})$},\quad\alpha^{v}_{s}\in(0,s),$ (B.2)

with a bounded, bilinear residual $R_{v}:V(\Omega_{t})^{2}\mapsto\mathbb{R}$.
Under assumption (4.14) the operator $A_{\varepsilon}(u^{\varepsilon}_{t})$ is
coercive (see (4.17)) and weakly continuous. Thus by the Brouwer fixed point
theorem, for small $s$ the variational equation (B.2) has a unique solution
$\tilde{u}^{\varepsilon}_{t+s}\in V(\Omega_{t})$.

Similarly, the optimality condition
$\partial_{u}\tilde{\mathcal{L}}^{\varepsilon}(s,u^{\varepsilon}_{t},\tilde{u}^{\varepsilon}_{t+s},\tilde{v}^{\varepsilon}_{t+s};\Omega_{t})=0$
reads as

$\int_{\Omega\setminus\Sigma_{t}}\bigl{(}(C\circ\phi_{s})E(\nabla\phi_{s}^{-1}\circ\phi_{s},\tilde{u})\cdot
E(\nabla\phi_{s}^{-1}\circ\phi_{s},\tilde{v}^{\varepsilon}_{t+s})\bigr{)}\,\omega^{\rm
d}_{s}dx\\\
+\int_{\Sigma_{t}}\int_{0}^{1}\Bigl{\\{}\Bigl{(}\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})[\\![\tilde{u}]\\!]_{\tilde{\tau}_{t+s}}\Bigr{)}\cdot[\\![\tilde{v}^{\varepsilon}_{t+s}]\\!]_{\tilde{\tau}_{t+s}}+[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])\,(\tilde{\nu}_{t+s}\cdot[\\![\tilde{u}]\\!])(\tilde{\nu}_{t+s}\cdot[\\![\tilde{v}^{\varepsilon}_{t+s}]\\!])\Bigr{\\}}\,\omega^{\rm
b}_{s}dr\,dS_{x}\\\ =\int_{\Gamma^{\rm
O}_{t}}(\tilde{u}^{\varepsilon}_{t+s}-z\circ\phi_{s})\cdot\tilde{u}\,\omega^{\rm
b}_{s}dS_{x}\quad\text{for all $\tilde{u}\in V(\Omega_{t})$}.$ (B.3)

The second inequality in (5.8) admits the decomposition for a weight
$\alpha^{u}_{s}\in(0,s)$:

$\langle
A_{\varepsilon}(u^{\varepsilon}_{t})\tilde{u},\tilde{v}^{\varepsilon}_{t+s}\rangle=\int_{\Gamma^{\rm
O}_{t}}(\tilde{u}^{\varepsilon}_{t+s}-z)\cdot\tilde{u}\,dS_{x}+sR_{u}(\alpha^{u}_{s},\tilde{v}^{\varepsilon}_{t+s},\tilde{u})\quad\text{for
all $\tilde{u}\in V(\Omega_{t})$},$ (B.4)

with bounded bilinear $R_{u}:V(\Omega_{t})^{2}\mapsto\mathbb{R}$, thus
possesses a unique solution $\tilde{v}^{\varepsilon}_{t+s}\in V(\Omega_{t})$,
for $s$ small enough.

##  Appendix C Proof of Lemma 4

### Uniform estimate of $\tilde{u}^{\varepsilon}_{t+s}$.

Testing the variational equation (B.1) with
$\tilde{v}=\tilde{u}^{\varepsilon}_{t+s}$ and applying the asymptotic
expansion (B.2) it follows

$\int_{\Omega\setminus\Sigma_{t}}\sigma(\tilde{u}^{\varepsilon}_{t+s})\cdot\epsilon(\tilde{u}^{\varepsilon}_{t+s})\,dx+\int_{\Sigma_{t}}\Bigl{\\{}\Bigl{(}\int_{0}^{1}\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})[\\![\tilde{u}^{\varepsilon}_{t+s}]\\!]_{\tau_{t}}\,dr+\nabla\alpha_{\rm
f}(0)\Bigr{)}\cdot[\\![\tilde{u}_{t+s}]\\!]_{\tau_{t}}\\\
+\Bigl{(}\int_{0}^{1}[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])(\nu_{t}\cdot[\\![\tilde{u}^{\varepsilon}_{t+s}]\\!])\,dr+[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](0)\Bigr{)}(\nu_{t}\cdot[\\![\tilde{u}_{t+s}]\\!])\Bigr{\\}}\,dS_{x}\\\
=\int_{\Gamma^{\rm
N}_{t}}g\cdot\tilde{u}^{\varepsilon}_{t+s}\,dS_{x}+sR_{v}(\alpha^{v}_{s},\tilde{u}^{\varepsilon}_{t+s},\tilde{u}^{\varepsilon}_{t+s}).$
(C.1)

We apply to (C.1) the Cauchy–Schwarz, Korn–Poincare (3.8) and trace
inequalities (3.17). By the virtue of boundedness of $\nabla\alpha_{\rm f}$,
$\nabla^{2}\alpha_{\rm f}$, $\alpha_{\rm c}^{\prime}$, $\alpha_{\rm
c}^{\prime\prime}$, $\beta_{\varepsilon}$ and
$\beta^{\prime}_{\varepsilon}\geq 0$ in (3.3), (3.5), (4.1), we derive the
estimate:

$(K_{\rm
fc2}-C_{1}|s|)\|\tilde{u}^{\varepsilon}_{t+s}\|_{H^{1}(\Omega\setminus\Sigma_{t})^{d}}\leq\sqrt{2}K_{\rm
tr}\bigl{(}\|g\|_{L^{2}(\Gamma^{\rm N}_{t})^{d}}+(K_{\rm f1}+K_{\rm
c1}-\beta_{\epsilon}(0))\sqrt{|\Sigma_{t}|}\bigr{)}+C_{1}|s|,\quad C_{1}>0,$
(C.2)

uniform in $\varepsilon$ and $s\leq s_{0}$ for sufficiently small $s_{0}>0$,
where $K_{\rm fc2}:=K_{\rm KP}-(K_{\rm f2}+K_{\rm c2})2K_{\rm tr}^{2}>0$ due
to the assumption (4.14).

### Uniform estimate of $\tilde{v}^{\varepsilon}_{t+s}$.

We test the variational equation (B.3) with
$\tilde{u}=\tilde{v}^{\varepsilon}_{t+s}$. and apply (B.4):

$\int_{\Omega\setminus\Sigma_{t}}\sigma(\tilde{v}^{\varepsilon}_{t+s})\cdot\epsilon(\tilde{v}^{\varepsilon}_{t+s})\,dx+\int_{\Sigma_{t}}\int_{0}^{1}\Bigl{\\{}\Bigl{(}\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})[\\![\tilde{v}^{\varepsilon}_{t+s}]\\!]_{\tau_{t}}\Bigr{)}\cdot[\\![\tilde{v}^{\varepsilon}_{t+s}]\\!]_{\tau_{t}}\\\
+[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])(\nu_{t}\cdot[\\![\tilde{v}^{\varepsilon}_{t+s}]\\!])^{2}\Bigr{\\}}\,dr\,dS_{x}=\int_{\Gamma^{\rm
O}_{t}}(\tilde{u}^{\varepsilon}_{t+s}-z)\cdot\tilde{v}^{\varepsilon}_{t+s}\,dS_{x}+sR_{u}(\alpha^{u}_{s},\tilde{v}^{\varepsilon}_{t+s},\tilde{v}^{\varepsilon}_{t+s}).$
(C.3)

With the help of Cauchy–Schwarz, Korn–Poincare and trace inequalities (3.8),
(3.17), due to the bondedness of $\nabla^{2}\alpha_{\rm f}$, $\alpha_{\rm
c}^{\prime\prime}$, $\beta^{\prime}_{\varepsilon}\geq 0$ in (3.3), (3.5),
(4.1), from (C.3) we derive the uniform estimate: there exists $C_{2}>0$ such
that

$(K_{\rm
fc2}-C_{2}|s|)\|\tilde{v}^{\varepsilon}_{t+s}\|_{H^{1}(\Omega\setminus\Sigma_{t})^{d}}\leq\sqrt{2}K_{\rm
tr}\|\tilde{u}^{\varepsilon}_{t+s}-z\|_{L^{2}(\Gamma^{\rm
O}_{t})^{d}}+C_{2}|s|.$ (C.4)

Thus, for small $|s|<K_{\rm fc2}/\min(C_{1},C_{2})$ relations (C.2) and (C.4)
together give

$\|\tilde{u}^{\varepsilon}_{t+s}\|_{H^{1}(\Omega\setminus\Sigma_{t})^{d}}+\|\tilde{v}^{\varepsilon}_{t+s}\|_{H^{1}(\Omega\setminus\Sigma_{t})^{d}}\leq
K,\quad K\geq 0.$ (C.5)

### Weak convergence of
$(\tilde{u}^{\varepsilon}_{t+s},\tilde{v}^{\varepsilon}_{t+s})$.

By the virtue of the uniform estimate (C.5), there exists a subsequence
$s_{k}\to 0$ as $k\to\infty$, and a weak accumulation point
$(\tilde{u}^{\varepsilon}_{t},\tilde{v}^{\varepsilon}_{t})\in
V(\Omega_{t})^{2}$ such that

$(\tilde{u}^{\varepsilon}_{t+s_{k}},\tilde{v}^{\varepsilon}_{t+s_{k}})\rightharpoonup(\tilde{u}^{\varepsilon}_{t},\tilde{v}^{\varepsilon}_{t})\quad\text{weakly
in
$H^{1}(\Omega\setminus\Sigma_{t})^{2d},H^{1/2}(\partial\Omega^{\pm}_{t})^{2d}$
as $s_{k}\to 0$}.$ (C.6)

By the compactness of embedding of the boundary traces it follows that

$(\tilde{u}^{\varepsilon}_{t+s_{k}},\tilde{v}^{\varepsilon}_{t+s_{k}})\to(\tilde{u}^{\varepsilon}_{t},\tilde{v}^{\varepsilon}_{t})\quad\text{strongly
in $L^{2}(\partial\Omega^{\pm}_{t})^{2d}$ as $s_{k}\to 0$}.$ (C.7)

Next we take the limit in (B.1) and (B.3) with $s=s_{k}$ as $k\to\infty$. Due
to the uniform continuity of $\nabla\alpha_{\rm f}$, $\alpha_{\rm
c}^{\prime},\beta_{\varepsilon}$ and $\nabla^{2}\alpha_{\rm f}$, $\alpha_{\rm
c}^{\prime\prime},\beta^{\prime}_{\varepsilon}$, and using (4) we arrive at
the variational equations (4.4) and (4.15), respectively. Therefore,
$(\tilde{u}^{\varepsilon}_{t},\tilde{v}^{\varepsilon}_{t})=(u^{\varepsilon}_{t},v^{\varepsilon}_{t})$.

### Strong convergence of $\tilde{u}^{\varepsilon}_{t+s}$.

With the help of asymptotic relation (C.1) and equation (4.4) with
$u=u^{\varepsilon}_{t}$, using the Korn–Poincare inequality (3.8), we
rearrange the terms as follows

$K_{\rm
KP}\|\tilde{u}^{\varepsilon}_{t+s}-u^{\varepsilon}_{t}\|^{2}_{H^{1}(\Omega\setminus\Sigma_{t})^{d}}\leq\int_{\Omega\setminus\Sigma_{t}}\sigma(\tilde{u}^{\varepsilon}_{t+s}-u^{\varepsilon}_{t})\cdot\varepsilon(\tilde{u}^{\varepsilon}_{t+s}-u^{\varepsilon}_{t})\,dx\\\
=\int_{\Omega\setminus\Sigma_{t}}\bigl{\\{}\sigma(\tilde{u}^{\varepsilon}_{t+s})\cdot\varepsilon(\tilde{u}^{\varepsilon}_{t+s})-\sigma(u^{\varepsilon}_{t})\cdot\varepsilon(u^{\varepsilon}_{t})-2\sigma(\tilde{u}^{\varepsilon}_{t+s}-u^{\varepsilon}_{t})\cdot\varepsilon(u^{\varepsilon}_{t})\bigr{\\}}\,dx=\int_{\Gamma^{\rm
N}_{t}}g\cdot(\tilde{u}^{\varepsilon}_{t+s}-u^{\varepsilon}_{t})\,dS_{x}\\\
-2\int_{\Omega\setminus\Sigma_{t}}\sigma(\tilde{u}^{\varepsilon}_{t+s}-u^{\varepsilon}_{t})\cdot\varepsilon(u^{\varepsilon}_{t})\,dx-\int_{\Sigma_{t}}\bigl{\\{}\bigl{(}\nabla\alpha_{\rm
f}([\\![\tilde{u}^{\varepsilon}_{t+s}]\\!]_{\tau_{t}})\cdot[\\![\tilde{u}^{\varepsilon}_{t+s}]\\!]_{\tau_{t}}-\nabla\alpha_{\rm
f}([\\![u^{\varepsilon}_{t}]\\!]_{\tau_{t}})\cdot[\\![u^{\varepsilon}_{t}]\\!]_{\tau_{t}}\bigr{)}\\\
+\bigl{(}[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![\tilde{u}^{\varepsilon}_{t+s}]\\!])(\nu_{t}\cdot[\\![\tilde{u}^{\varepsilon}_{t+s}]\\!]-[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!])(\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!]\bigr{)}\bigr{\\}}\,dS_{x}+{\rm
O}(|s|).$ (C.8)

Taking the limit in (C.8) as $s_{k}\to 0$, due to the convergence established
in (C.6) and (C.7), we conclude that

$\|\tilde{u}^{\varepsilon}_{t+s_{k}}-u^{\varepsilon}_{t}\|_{H^{1}(\Omega\setminus\Sigma_{t})^{d}}\to
0\quad\text{as $s_{k}\to 0$}.$ (C.9)

### Strong convergence of $\tilde{v}^{\varepsilon}_{t+s}$.

We subtract equation (4.15) from (B.3) and use asymptotic expansions (A.1)
such that

$\int_{\Omega\setminus\Sigma_{t}}\varepsilon(\tilde{v})\cdot\sigma(\tilde{v}^{\varepsilon}_{t+s}-v^{\varepsilon}_{t})\,dx=\int_{\Sigma_{t}}\int_{0}^{1}\Bigl{\\{}\Bigl{(}\nabla^{2}\alpha_{\rm
f}([\\![r\tilde{u}^{\varepsilon}_{t+s}]\\!]_{\tau_{t}})[\\![\tilde{v}^{\varepsilon}_{t+s}]\\!]_{\tau_{t}}-\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})[\\![v^{\varepsilon}_{t}]\\!]_{\tau_{t}}\Bigr{)}\cdot[\\![\tilde{v}]\\!]_{\tau_{t}}\\\
+\bigl{(}[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![r\tilde{u}^{\varepsilon}_{t+s}]\\!])(\nu_{t}\cdot[\\![\tilde{v}^{\varepsilon}_{t+s}]\\!])-[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])(\nu_{t}\cdot[\\![v^{\varepsilon}_{t}]\\!])\bigr{)}(\nu_{t}\cdot[\\![\tilde{v}]\\!])\Bigr{\\}}\,dr\,dS_{x}+{\rm
O}(|s|).$ (C.10)

Applying to (C.10) the Cauchy–Schwarz inequality, due to the properties of
$\nabla^{2}\alpha_{\rm f}$, $\alpha_{\rm c}^{\prime\prime}$,
$\beta^{\prime}_{\varepsilon}$ in (3.3), (3.5), (4.1), we obtain the upper
bound

$\int_{\Omega\setminus\Sigma_{t}}\varepsilon(\tilde{v})\cdot\sigma(\tilde{v}^{\varepsilon}_{t+s}-v^{\varepsilon}_{t})\,dx\leq
K_{\rm
f2}\|[\\![\tilde{v}^{\varepsilon}_{t+s}-v^{\varepsilon}_{t}]\\!]_{\tau_{t}}\|_{L^{2}(\Sigma_{t})^{d}}\|[\\![\tilde{v}]\\!]_{\tau_{t}}\|_{L^{2}(\Sigma_{t})^{d}}\\\
+\int_{0}^{1}\Bigl{\\{}\bigl{\|}\nabla^{2}\alpha_{\rm
f}([\\![r\tilde{u}^{\varepsilon}_{t+s})]\\!]_{\tau_{t}}-\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})\bigr{\|}_{L^{2}(\Sigma_{t})^{d\times
d}}\|[\\![v^{\varepsilon}_{t}]\\!]_{\tau_{t}}\|_{L^{4}(\Sigma_{t})^{d}}\|[\\![\tilde{v}]\\!]_{\tau_{t}}\|_{L^{4}(\Sigma_{t})^{d}}\\\
+\bigl{\|}[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![r\tilde{u}^{\varepsilon}_{t+s})]\\!]-[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])\bigr{\|}_{L^{2}(\Sigma_{t})}\|\nu_{t}\cdot[\\![v^{\varepsilon}_{t}]\\!]\|_{L^{4}(\Sigma_{t})}\|\nu_{t}\cdot[\\![\tilde{v}]\\!]\|_{L^{4}(\Sigma_{t})}\Bigr{\\}}\,dr\\\
+\bigl{(}K_{\rm c2}+\frac{K_{\beta
1}}{\varepsilon}\bigr{)}\|\nu_{t}\cdot[\\![\tilde{v}^{\varepsilon}_{t+s}-v^{\varepsilon}_{t}]\\!]\|_{L^{2}(\Sigma_{t})}\|\nu_{t}\cdot[\\![\tilde{v}]\\!]\|_{L^{2}(\Sigma_{t})}+C|s|,\quad
C>0.$ (C.11)

By the Sobolev embedding theorem the continuity property holds:

$\|u\|_{L^{4}(\partial\Omega^{\pm}_{t})^{d}}\leq K_{\rm
emb}\|u\|_{H^{1/2}(\partial\Omega^{\pm}_{t})^{d}},\quad u\in
H^{1}(\Omega^{\pm}_{t})^{d},\quad d=2,3.$ (C.12)

Then (C.12), Korn–Poincare and trace inequalities (3.8), (3.17), together with
convergences (C.6), (C.7) guarantee that for fixed $\varepsilon$:

$K_{\rm
KP}\|\tilde{v}^{\varepsilon}_{t+s_{k}}-v^{\varepsilon}_{t}\|_{H^{1}(\Omega\setminus\Sigma_{t})^{d}}\leq\sup_{\tilde{v}\in
V(\Omega_{t})}\frac{1}{\|\tilde{v}\|_{H^{1}(\Omega\setminus\Sigma_{t})^{d}}}\int_{\Omega\setminus\Sigma_{t}}\varepsilon(\tilde{v})\cdot\sigma(\tilde{v}^{\varepsilon}_{t+s_{k}}-v^{\varepsilon}_{t})\,dx\to
0\text{ as $s_{k}\to 0$}.$ (C.13)

The proof of Lemma 4 is complete.

##  Appendix D Proof of Corollary 1

Let $(u^{\varepsilon}_{t},v^{\varepsilon}_{t})\in
H^{2}(\Omega^{+}_{t})^{2d}\cap H^{2}(\Omega^{-}_{t})^{2d}$ be a solution to
(4.4) and (4.15). We integrate by parts the domain integral over
$\Omega\setminus\Sigma_{t}$ from (5.20) at $\tau=0$ so that

$I(\Omega\setminus\Sigma_{t}):=-\int_{\Omega^{\pm}_{t}}\bigl{(}({\rm
div}\Lambda C+\nabla
C\Lambda)\epsilon(u^{\varepsilon}_{t})\cdot\epsilon(v^{\varepsilon}_{t})-\sigma(u^{\varepsilon}_{t})\cdot
E(\nabla\Lambda,v^{\varepsilon}_{t})-\sigma(v^{\varepsilon}_{t})\cdot
E(\nabla\Lambda,u^{\varepsilon}_{t})\bigr{)}dx\\\
=-\int_{\partial\Omega^{\pm}_{t}}\Lambda\cdot\bigl{(}n^{\pm}_{t}\sigma(u^{\varepsilon}_{t})\cdot\epsilon(v^{\varepsilon}_{t})-\nabla(u^{\varepsilon}_{t})^{\top}\sigma(v^{\varepsilon}_{t})n^{\pm}_{t}-\nabla(v^{\varepsilon}_{t})^{\top}\sigma(u^{\varepsilon}_{t})n^{\pm}_{t}\bigr{)}\,dS_{x}=\int_{\Sigma_{t}}\Lambda\cdot\bigl{(}\nu_{t}[\\![\sigma(u^{\varepsilon}_{t})\cdot\epsilon(v^{\varepsilon}_{t})]\\!]\\\
-[\\![\nabla(u^{\varepsilon}_{t})^{\top}\sigma(v^{\varepsilon}_{t})]\\!]\nu_{t}-[\\![\nabla(v^{\varepsilon}_{t})^{\top}\sigma(u^{\varepsilon}_{t})]\\!]\nu_{t}\bigr{)}\,dS_{x}+\int_{\Gamma^{\rm
D}_{t}\cup\Gamma^{\rm
N}_{t}}\\!\\!\Lambda\cdot\bigl{(}\nabla(u^{\varepsilon}_{t})^{\top}\sigma(v^{\varepsilon}_{t})n_{t}+\nabla(v^{\varepsilon}_{t})^{\top}\sigma(u^{\varepsilon}_{t})n_{t}\bigr{)}\,dS_{x},$

where we use the assumption $n_{t}\cdot\Lambda=0$ at $\partial\Omega$. Using
boundary conditions from (4.1), (4.2) and the notation $\mathcal{D}_{1}$ from
(5.24) it follows that

$I(\Omega\setminus\Sigma_{t})=\int_{\Sigma_{t}}\Lambda\cdot\Bigl{(}\nu_{t}[\\![\sigma(u^{\varepsilon}_{t})\cdot\epsilon(v^{\varepsilon}_{t})]\\!]-[\\![\nabla
v^{\varepsilon}_{t}]\\!]^{\top}\bigl{(}\nabla\alpha_{\rm
f}([\\![u^{\varepsilon}_{t}]\\!]_{\tau_{t}})+[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!])\,\nu_{t}\bigr{)}\\\
-[\\![\nabla
u^{\varepsilon}_{t}]\\!]^{\top}\int_{0}^{1}\bigl{(}\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})\,[\\![v^{\varepsilon}_{t}]\\!]_{\tau_{t}}+[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])\,(\nu_{t}\cdot[\\![v^{\varepsilon}_{t}]\\!])\,\nu_{t}\bigr{)}\,dr\Bigr{)}\,dS_{x}\\\
+\int_{\Gamma^{\rm
O}_{t}}\Lambda\cdot\bigl{(}\nabla(u^{\varepsilon}_{t})^{\top}(u^{\varepsilon}_{t}-z)\bigr{)}\,dS_{x}+\int_{\Gamma^{\rm
N}_{t}}\Lambda\cdot(\nabla(v^{\varepsilon}_{t})^{\top}g)\,dS_{x}+\int_{\Gamma^{\rm
D}_{t}}\Lambda\cdot\mathcal{D}_{1}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})\,dS_{x}.$
(D.1)

After substitution of (D.1) into (5.20), the integrand at $\Sigma_{t}$ is
gathered in the expression:

$I_{\Sigma_{t}}:=-{\rm div}_{\tau_{t}}\Lambda\,\bigl{\\{}\nabla\alpha_{\rm
f}([\\![u^{\varepsilon}_{t}]\\!]_{\tau_{t}})\cdot[\\![v^{\varepsilon}_{t}]\\!]_{\tau_{t}}+[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!])\,(\nu_{t}\cdot[\\![v^{\varepsilon}_{t}]\\!])\bigr{\\}}\\\
+\Lambda\cdot\bigl{\\{}\nu_{t}[\\![\sigma(u^{\varepsilon}_{t})\cdot\epsilon(v^{\varepsilon}_{t})]\\!]-\bigl{(}[\\![\nabla
v^{\varepsilon}_{t}]\\!]^{\top}-(\nu_{t}\cdot[\\![v^{\varepsilon}_{t}]\\!])\nabla\nu_{t}^{\top}-\nabla\nu_{t}^{\top}[\\![v^{\varepsilon}_{t}]\\!]\nu_{t}^{\top}\bigr{)}\nabla\alpha_{\rm
f}([\\![u^{\varepsilon}_{t}]\\!]_{\tau_{t}})\\\ -\bigl{(}[\\![\nabla
v^{\varepsilon}_{t}]\\!]^{\top}\nu_{t}+\nabla\nu_{t}^{\top}[\\![v^{\varepsilon}_{t}]\\!]\bigr{)}\,[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!])-\bigl{(}[\\![\nabla
u^{\varepsilon}_{t}]\\!]^{\top}\nu_{t}+\nabla\nu_{t}^{\top}[\\![u^{\varepsilon}_{t}]\\!]\bigr{)}\int_{0}^{1}[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])(\nu_{t}\cdot[\\![v^{\varepsilon}_{t}]\\!])\,dr\\\
-\bigl{(}[\\![\nabla
u^{\varepsilon}_{t}]\\!]^{\top}-(\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!])\nabla\nu_{t}^{\top}-\nabla\nu_{t}^{\top}[\\![u^{\varepsilon}_{t}]\\!]\nu_{t}^{\top}\bigr{)}\int_{0}^{1}\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})[\\![v^{\varepsilon}_{t}]\\!]_{\tau_{t}}\,dr\bigr{\\}}.$
(D.2)

In order to combine like terms, we exploit the calculus

$\Lambda\cdot\nabla(\xi\cdot\eta)=\Lambda\cdot(\nabla\xi^{\top}\eta+\nabla\eta^{\top}\xi)=\eta\cdot\nabla\xi\Lambda+\xi\cdot\nabla\eta\Lambda\quad\text{for
$\xi,\eta\in\mathbb{R}^{d}$}.$ (D.3)

With the help of (D.3), the gradient of the product due to friction term is
calculated:

$p_{\rm f}(\tilde{u},\tilde{v}):=\nabla\alpha_{\rm
f}([\\![\tilde{u}]\\!]_{\tau_{t}})\cdot[\\![\tilde{v}]\\!]_{\tau_{t}},\quad\nabla
p_{\rm
f}(\tilde{u},\tilde{v})=\nabla([\\![\tilde{v}]\\!]_{\tau_{t}})^{\top}\nabla\alpha_{\rm
f}([\\![\tilde{u}]\\!]_{\tau_{t}})+\nabla([\\![\tilde{u}]\\!]_{\tau_{t}})^{\top}\nabla^{2}\alpha_{\rm
f}([\\![\tilde{u}]\\!]_{\tau_{t}})\,[\\![\tilde{v}]\\!]_{\tau_{t}},$ (D.4)

where
$\nabla([\\![\tilde{u}]\\!]_{\tau_{t}})^{\top}=[\\![\nabla\tilde{u}]\\!]^{\top}-(\nu_{t}\cdot[\\![\tilde{u}]\\!])\nabla\nu_{t}^{\top}-\nabla(\nu_{t}\cdot[\\![\tilde{u}]\\!])\nu_{t}^{\top}$
at $\Sigma_{t}$ according to (3.1). Similarly, we compute the gradient for the
cohesive term

$p^{\varepsilon}_{\rm c}(\tilde{u},\tilde{v}):=[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![\tilde{u}]\\!])\,(\nu_{t}\cdot[\\![\tilde{v}]\\!]),\quad\nabla
p^{\varepsilon}_{\rm
c}(\tilde{u},\tilde{v})=([\\![\nabla\tilde{v}]\\!]^{\top}\nu_{t}+\nabla\nu_{t}^{\top}[\\![\tilde{v}]\\!])\,[\alpha_{\rm
c}^{\prime}+\beta_{\varepsilon}](\nu_{t}\cdot[\\![\tilde{u}]\\!])\\\
+([\\![\nabla\tilde{u}]\\!]^{\top}\nu_{t}+\nabla\nu_{t}^{\top}[\\![\tilde{u}]\\!])\,[\alpha_{\rm
c}^{\prime\prime}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![\tilde{u}]\\!])\,(\nu_{t}\cdot[\\![\tilde{v}]\\!]).$
(D.5)

By (D.4) and (D.5), the integrand (D.2) is expressed as

$I_{\Sigma_{t}}=-{\rm div}_{\tau_{t}}\Lambda\,[p_{\rm f}+p^{\varepsilon}_{\rm
c}](u^{\varepsilon}_{t},v^{\varepsilon}_{t})+\Lambda\cdot\bigl{\\{}\nu_{t}[\\![\sigma(u^{\varepsilon}_{t})\cdot\epsilon(v^{\varepsilon}_{t})]\\!]-\nabla[p_{\rm
f}+p^{\varepsilon}_{\rm
c}](u^{\varepsilon}_{t},v^{\varepsilon}_{t})-[\\![\nabla
v^{\varepsilon}_{t}]\\!]^{\top}\nu_{t}\bigl{(}\nu_{t}\cdot\nabla\alpha_{\rm
f}([\\![u^{\varepsilon}_{t}]\\!]_{\tau_{t}})\bigr{)}\\\ -[\\![\nabla
u^{\varepsilon}_{t}]\\!]^{\top}\nu_{t}\Bigl{(}\nu_{t}\cdot\int_{0}^{1}\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})[\\![v^{\varepsilon}_{t}]\\!]_{\tau_{t}}dr\Bigr{)}-\nabla(\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!])^{\top}\int_{0}^{1}\Bigl{(}[\alpha^{\prime\prime}_{\rm
c}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![ru^{\varepsilon}_{t}]\\!])\\\
-[\alpha^{\prime\prime}_{\rm
c}+\beta^{\prime}_{\varepsilon}](\nu_{t}\cdot[\\![u^{\varepsilon}_{t}]\\!])\Bigr{)}(\nu_{t}\cdot[\\![v^{\varepsilon}_{t}]\\!])\,dr\bigr{\\}}-\nabla([\\![u^{\varepsilon}_{t}]\\!]_{\tau_{t}})^{\top}\int_{0}^{1}\bigl{(}\nabla^{2}\alpha_{\rm
f}([\\![ru^{\varepsilon}_{t}]\\!]_{\tau_{t}})-\nabla^{2}\alpha_{\rm
f}([\\![u^{\varepsilon}_{t}]\\!]_{\tau_{t}})\bigr{)}[\\![v^{\varepsilon}_{t}]\\!]_{\tau_{t}}\,dr.$
(D.6)

Introducing for short the notation of $q_{\rm f},q^{\varepsilon}_{\rm c}$ in
(5.26) which is based on (D.6), we rearrange the terms in the shape derivative
in the form

${\textstyle\frac{\partial}{\partial
s}}\tilde{\mathcal{L}}^{\varepsilon}(0,u^{\varepsilon}_{t},u^{\varepsilon}_{t},v^{\varepsilon}_{t};\Omega_{t})=\frac{1}{2}\int_{\Gamma^{\rm
O}_{t}}\bigl{(}{\rm
div}_{\tau_{t}}\Lambda\,|u^{\varepsilon}_{t}-z|^{2}+\Lambda\cdot\nabla(|u^{\varepsilon}_{t}-z|^{2})\bigr{)}dS_{x}+\rho\int_{\Sigma_{t}}{\rm
div}_{\tau_{t}}\Lambda\,dS_{x}\\\ +\int_{\Sigma_{t}}\bigl{\\{}-{\rm
div}_{\tau_{t}}\Lambda\,[p_{\rm f}+p^{\varepsilon}_{\rm
c}](u^{\varepsilon}_{t},v^{\varepsilon}_{t})+\Lambda\cdot\bigl{(}\nu_{t}[\\![\sigma(u^{\varepsilon}_{t})\cdot\epsilon(v^{\varepsilon}_{t})]\\!]-[\nabla(p_{\rm
f}+p^{\varepsilon}_{\rm c})+q_{\rm f}+q^{\varepsilon}_{\rm
c}](u^{\varepsilon}_{t},v^{\varepsilon}_{t})\bigr{)}\bigr{\\}}\,dS_{x}\\\
+\int_{\Gamma^{\rm N}_{t}}\bigl{(}{\rm div}_{\tau_{t}}\Lambda(g\cdot
v^{\varepsilon}_{t})+\Lambda\cdot\nabla(g\cdot
v^{\varepsilon}_{t})\bigr{)}\,dS_{x}+\int_{\Gamma^{\rm
D}_{t}}\Lambda\cdot\mathcal{D}_{1}(u^{\varepsilon}_{t},v^{\varepsilon}_{t})\,dS_{x}.$
(D.7)

Since the tangential velocity, its tangential divergence, and the curvature
are equal to

$\Lambda_{\tau_{t}}=\Lambda-(n^{\pm}_{t}\cdot\Lambda)n^{\pm}_{t},\quad{\rm
div}_{\tau_{t}}\Lambda_{\tau_{t}}={\rm
div}_{\tau_{t}}\Lambda-(n^{\pm}_{t}\cdot\Lambda)\varkappa^{\pm}_{t},\quad\varkappa^{\pm}_{t}={\rm
div}_{\tau_{t}}n^{\pm}_{t}\text{ at }\partial\Omega^{\pm}_{t},$ (D.8)

for smooth $p$ the integration along a boundary
$\Gamma_{t}\subset\partial\Omega^{\pm}_{t}$ is given by the formula (see e.g.
(SZ/92, , (2.125))):

$\int_{\Gamma_{t}}({\rm div}_{\tau_{t}}\Lambda\,p+\Lambda\cdot\nabla
p)\,dS_{x}=\int_{\Gamma_{t}}(n_{t}\cdot\Lambda)(\varkappa_{t}p+n_{t}\cdot\nabla
p)\,dS_{x}+\begin{cases}(\tau_{t}\cdot\Lambda)p|_{\partial\Gamma_{t}}&\text{in
2D,}\\\
{\displaystyle\int_{\partial\Gamma_{t}}}(b_{t}\cdot\Lambda)p\,dL_{x}&\text{in
3D.}\end{cases}$ (D.9)

In (D.9) $\tau_{t}$ is a tangential vector at $\partial\Gamma_{t}$ positively
oriented to $n_{t}$ in 2D, and $b_{t}=\tau_{t}\times n_{t}$ is a binomial
vector within the moving frame at $\partial\Gamma_{t}$ in 3D. Applying (D.9)
to (D.7), decomposing the vectors in (5.21) into the normal and tangential