Phase field topology optimisation for 4D printing

Harald Garcke ¹¹1Fakultät für Mathematik, Universität Regensburg, 93040 Regensburg, Germany (H[email protected]). Kei Fong Lam ²²2Department of Mathematics, Hong Kong Baptist University, Kowloon Tong, Hong Kong (a[email protected]). Robert Nürnberg ³³3Department of Mathematics, University of Trento, Trento, Italy (r[email protected]). Andrea Signori ⁴⁴4Department of Mathematics, Politecnico di Milano, 20133 Milano, Italy (a[email protected]).

Abstract

This work concerns a structural topology optimisation problem for 4D printing based on the phase field approach. The concept of 4D printing as a targeted evolution of 3D printed structures can be realised in a two-step process. One first fabricates a 3D object with multi-material active composites and apply external loads in the programming stage. Then, a change in an environmental stimulus and the removal of loads cause the object deform in the programmed stage. The dynamic transition between the original and deformed shapes is achieved with appropriate applications of the stimulus. The mathematical interest is to find an optimal distribution for the materials such that the 3D printed object achieves a targeted configuration in the programmed stage as best as possible.

Casting the problem as a PDE-constrained minimisation problem, we consider a vector-valued order parameter representing the volume fractions of the different materials in the composite as a control variable. We prove the existence of optimal designs and formulate first order necessary conditions for minimisers. Moreover, by suitable asymptotic techniques, we relate our approach to a sharp interface description. Finally, the theoretical results are validated by several numerical simulations both in two and three space dimensions.

Keywords: 4D printing, Printed active composites, Topology optimisation, Phase field, Linear elasticity, Optimal control AMS (MOS) Subject Classification: 49J20, 49K40, 49J50

1 Introduction

Four-dimensional (4D) printing [2, 52, 61] entails the combination of additive manufacturing (3D printing) and active material technologies to create printed composites capable of morphing into different configurations in response to various environmental stimuli. First designs of such composites consist of active material components, such as piezoelectric ceramics, hydrogels or shape memory polymers [24], in the form of fibers integrated within a passive elastomeric matrix [37]. These multi-material active composites were originally difficult to manufacture, owing to the fragility of the materials involved [44]. However, with 3D printing techniques it is nowadays feasible to fabricate these active composites to a high degree of precision, resulting in so-called printed active composites (PACs) [50]. For an overview of other 4D printing strategies besides PACs in the construction of smart materials allowing direct stimuli-responsive transformations, we refer to [67].

The shape shifting functionality of the active components enables the self-actuating and self-assembling potentials of PACs, allowing them to fold, bend, twist, expand and contract when a stimulus is applied, and return to their original configurations after the stimulus is removed. This property has led to the fabrication of intelligent active hinges and origami-like objects [35, 37], mesh structures [24, 64] and self-actuated deformable solids [59] in the form of in the form of graspers and smart key-lock systems. We refer to the review article [52] and the references cited therein for more applications of 4D printing. The shape memory behaviour of the PACs can be programmed in a two-step cycle: The first (programming) step involves deforming the structure from its permanent shape to a metastable temporary shape, and the second (recovery) step involves applying an appropriate stimulus so that the structure regains its original shape. A typical stimulus is heat (in combination with light [36] or water [7]), in which programmed PACs alter their shapes when the temperature rises above or drops below a critical value.

With the advances in the state-of-the art 3D printing technologies, the designs of PACs need not be limited to the conventional fibre-matrix architectures first considered in [37]. In particular, the distribution of active and passive materials in the designs can take on more complicated geometries to better fulfil the intended functionalities of the PACs. This opens up the possibility of a computational design approach guided by a structural topology optimisation framework. In the context of 3D printing, see, e.g., the review [26], this framework has been applied to explore optimising support structures to overhang regions [3, 46, 51], as well as self-support designs respecting the overhang angle constraints [20, 32, 43, 45]. For active materials and active composites, [40, 54] studied how to pattern thin-film layers within a multi-layer structure with the aim of generating large shape changes via spatially varying eigenstrains within the microstructures, while [50] aimed to optimise the microstructures of PACs matching various target shapes after a thermomechanical training and activation cycle. Later works incorporated nonlinear thermoelasticity [39, 59], thermo-mechanical cycles of shape memory polymers [12], reversible deformations [49], as well as multi-material designs [65] within the topology optimisation framework.

In many of the aforementioned contributions, the topology optimisation is implemented numerically with the level-set method or the solid isotropic material with penalisation (SIMP) approach. In this work we employ an alternative approach based on the phase field methodology [17], which allows a straightforward extension to the multiphase setting [13, 62] involving multiple (possibly distinct) types of active materials within the design. In particular, this opens up the design to multiphase PACs that can memorise more than two shapes [38, 47, 58, 63, 66]. The phase field-based structural topology optimisation approach has been popularised in recent years by many authors, with applications in nonlinear elasticity [55], stress constraints [19], compliance optimisation [15, 60], elastoplasticity [4], eigenfrequency maximisation [31, 60], graded-material design [21], shape optimisation in fluid flow [28, 29, 30] and more recently for 3D printing with overhang angle constraints [32].

Taking inspiration from the setting of Maute et al. [50], we formulate a structural topology optimisation problem for a multiphase PAC with the objective of finding optimal distributions of active and passive materials so that the resulting composite matches targeted shapes as close as possible. An additional perimeter regularisation term, in the form of a multiphase Ginzburg–Landau functional, is added, and our contribution involves a mathematical analysis of the resulting multiphase structural topology optimisation problem with emphasis on the rigorous derivation of minimisers and optimality conditions. A sharp interface asymptotic analysis is performed to obtain a set of optimality conditions applicable in a level set-based shape optimisation framework. We perform numerical simulations in two and three spatial dimensions to show the optimal distributions of active and passive components in order to match with various target shapes for the PACs.

The rest of this paper is organised as follows: in Section 2 we formulate the phase field structural optimisation problem to be studied, and present several preliminary mathematical results. In Sections 3 and 4 we analyse the design optimisation problem and establish analytical results concerning minimisers and optimality conditions. The sharp interface limit is explored in Section 5 and, finally, in Section 6 we present the numerical discretisation and several simulations of our approach.

2 Problem formulation

Within a bounded domain $\Omega\subset{\mathbb{R}}^{d}$ , $d\in\{2,3\}$ , with Lipschitz boundary $\Gamma:=\partial\Omega$ , we assume there are $L$ types of linearly elastic materials, whose volume fractions are encoded with the help of a vectorial phase field variable ${\bm{\varphi}}=(\varphi_{1},\dots,\varphi_{L}):\Omega\to\Delta^{L}$ , where $\Delta^{L}$ denotes the Gibbs simplex in ${\mathbb{R}}^{L}$ :

\Delta^{L}:=\Big{\{}\bm{x}=(x_{1},...,x_{L})\in{\mathbb{R}}^{L}\,:\,\sum_{i=1}^{L}x_{i}=1,\,x_{i}\geq 0\text{ for all }i\in\{1,\dots,L\}\Big{\}}.

For our application to PACs, we take $\varphi_{L}$ as the volume fraction of the passive elastic material, and $\varphi_{1},\dots,\varphi_{L-1}$ as the volume fractions of (possibly different) active elastic materials. Note that in the two-phase case $L=2$ , we simply have ${\bm{\varphi}}=(\varphi_{1},\varphi_{2})$ , and due to the relation $\varphi_{1}+\varphi_{2}=1$ we may instead use the scalar difference function $\varphi:=\varphi_{1}-\varphi_{2}$ to encode ${\bm{\varphi}}$ via the relation ${\bm{\varphi}}=(\frac{1}{2}(1+\varphi),\frac{1}{2}(1-\varphi))$ . This particular scenario will be employed later on, when dealing with the connection between the problem we are going to analyse and the corresponding sharp interface limit in Section 5, as well as for the numerical simulations presented in Section 6.

The shape shifting mechanism considered in [50, 68] involves two levels of temperature and one set of external loads, with one temperature $T_{H}$ higher than a critical transition temperature $T_{g}$ of the active materials (e.g., the glass transition temperature for shape memory polymers), and the other temperature $T_{L}$ lower than the critical temperature. The printed composite is first heated to $T_{H}$ , and the shape memory cycle starts at $T_{H}$ and proceeds as follows: First, external loads are applied to deform the printed composite while the temperature remains at $T_{H}$ , with the new configuration being known as the programming stage (or Stage 1). Then, the temperature is decreased while the loads are maintained on the printed composite, which are then removed once the temperature reached $T_{L}$ . The resulting shape at $T_{L}$ is the desired shape and we denote it as the programmed stage (or Stage 2). Increasing the temperature to $T_{H}$ enables the printed composite to recover its original shape, and this ends the shape memory cycle, see Figure 1 for the thermo-mechanical processing steps involving the two stages.

Refer to caption — Figure 1: Schematics of the the shape memory cycle from [50] involving a programming stage (Stage 1) and a programmed stage (Stage 2).

To capture the above behaviour, following [50] we consider a model for each stage. In the programming stage (Stage 1), we consider an elastic displacement ${\overline{{\bf u}}}:\Omega\to{\mathbb{R}}^{d}$ and decompose the domain boundary $\Gamma$ into a partition $\Gamma={\rm cl}(\overline{\Gamma}_{D})\cup{\rm cl}(\overline{\Gamma}_{N})$ with relative open subsets $\overline{\Gamma}_{D}$ and $\overline{\Gamma}_{N}$ such that $\overline{\Gamma}_{D}\cap\overline{\Gamma}_{N}=\emptyset$ and $\overline{\Gamma}_{D}\neq\emptyset$ , where ${\rm cl}(A)$ denotes the closure of a set $A$ , and we assign a prescribed displacement $\overline{\bm{U}}$ on $\overline{\Gamma}_{D}$ and surface loads $\overline{{\bf g}}$ on $\overline{\Gamma}_{N}$ . Under a linearised elasticity setting, the balance of momentum yields the following system of equations for the displacement ${\overline{{\bf u}}}$ :


$\displaystyle-\mathop{\rm div}\nolimits\big{(}\overline{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}})\big{)}$	$\displaystyle=\overline{{\bf F}}$	$\displaystyle\quad\text{ in }\Omega,$	(2.1a)
$\displaystyle{\overline{{\bf u}}}$	$\displaystyle=\overline{\bm{U}}$	$\displaystyle\quad\text{ on }\overline{\Gamma}_{D},$	(2.1b)
$\displaystyle\big{(}\overline{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}})\big{)}{\bf n}$	$\displaystyle=\overline{{\bf g}}$	$\displaystyle\quad\text{ on }\overline{\Gamma}_{N},$	(2.1c)

with a phase-dependent elasticity tensor $\overline{{\mathbb{C}}}$ , body force $\overline{{\bf F}}$ , outer unit normal ${\bf n}$ , and symmetrised gradient $\mathcal{E}({\overline{{\bf u}}})$ . One example of $\overline{{\mathbb{C}}}({\bm{\varphi}})$ is

\overline{{\mathbb{C}}}({\bm{\varphi}}(\bm{x}))=\sum_{i=1}^{L}\overline{{\mathbb{C}}}_{i}\varphi_{i}(\bm{x})\quad\text{ for }{\bm{\varphi}}(\bm{x})=(\varphi_{1}(\bm{x}),\dots,\varphi_{L}(\bm{x}))\in\Delta^{L},\quad\bm{x}\in\Omega,

with constant tensors $\overline{{\mathbb{C}}}_{i}$ , $1\leq i\leq L$ .

After the change in temperature from $T_{H}$ to $T_{L}$ and after the programming loads in Stage 1 have been removed, the PAC experiences deformations due to residual stresses generated during the thermomechanical processing steps. When the temperature falls below $T_{g}$ , the active elastic materials undergo a phase transition from a soft rubbery state to a glassy state that has a higher Young’s modulus. We introduce a new variable ${\widehat{{\bf u}}}:\Omega\to{\mathbb{R}}^{d}$ to denote the displacement in the programmed stage (Stage 2), and as in [50], model the strains from the programming stage (Stage 1) as eigenstrains for ${\widehat{{\bf u}}}$ . These eigenstrains are present only in the regions of active elastic materials, which we model with a fixity function $\chi:{\mathbb{R}}^{L}\to[0,\infty)$ . The shape fixity for a shape memory material is the ratio (expressed as a percentage) between the strain in the stress-free state after the programming step and the maximum strain [1]. For example, if the deformation elongates the material, the fixity quantifies the ability of the material to hold the temporary elongated length when the stress is removed. It is clear from the definition that for a passive elastic material, the fixity is zero, and so we set that $\chi=0$ in the region $\{\varphi_{L}=1\}$ of the passive elastic material. Decomposing the domain boundary $\Gamma$ into a possibly different partition $\Gamma={\rm cl}(\widehat{\Gamma}_{D})\cup{\rm cl}(\widehat{\Gamma}_{N})$ with relative open subsets $\widehat{\Gamma}_{D}$ and $\widehat{\Gamma}_{N}$ such that $\widehat{\Gamma}_{D}\cap\widehat{\Gamma}_{N}=\emptyset$ and $\widehat{\Gamma}_{D}\neq\emptyset$ , where we assign a prescribed displacement $\widehat{\bm{U}}$ on $\widehat{\Gamma}_{D}$ and surface loads $\widehat{{\bf g}}$ on $\widehat{\Gamma}_{N}$ , the equations for the programmed stage (Stage 2) read as


$\displaystyle-\mathop{\rm div}\nolimits\big{(}\widehat{{\mathbb{C}}}({\bm{\varphi}})(\mathcal{E}({\widehat{{\bf u}}})-\chi({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}}))\big{)}$	$\displaystyle=\widehat{{\bf F}}$	$\displaystyle\text{ in }\Omega,$	(2.2a)
$\displaystyle{\widehat{{\bf u}}}$	$\displaystyle=\widehat{\bm{U}}$	$\displaystyle\text{ on }\widehat{\Gamma}_{D},$	(2.2b)
$\displaystyle\big{(}\widehat{{\mathbb{C}}}({\bm{\varphi}})(\mathcal{E}({\widehat{{\bf u}}})-\chi({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}}))\big{)}{\bf n}$	$\displaystyle=\widehat{{\bf g}}$	$\displaystyle\text{ on }\widehat{\Gamma}_{N},$	(2.2c)

with a phase-dependent elasticity tensor $\widehat{{\mathbb{C}}}$ and body force $\widehat{{\bf F}}$ . In the above, the change in the elasticity tensor from $\overline{{\mathbb{C}}}$ in Stage 1 to $\widehat{{\mathbb{C}}}$ in Stage 2 encodes the change in the elastic properties of the PAC when the temperature changes from $T_{H}$ to $T_{L}$ . Similarly to [50], here we have neglected the strains arising from thermal expansion in (2.1) and (2.2).

In the next section, under a suitable functional framework, we demonstrate that (2.1) and (2.2) are uniquely solvable, with the solution depending continuously on ${\bm{\varphi}}$ . Since ${\bm{\varphi}}$ controls the distribution of the passive and active elastic materials, it is natural to ask for specific material distributions that optimise certain cost functionals related to the design of PACs. Motivated from [50], we primarily focus on the following cost functional

\displaystyle J({\bm{\varphi}},{\widehat{{\bf u}}}):=\frac{1}{2}\int_{\Gamma^{\rm tar}}\Big{(}W({\widehat{{\bf u}}}-\bm{u}^{\rm tar})\Big{)}\cdot({\widehat{{\bf u}}}-\bm{u}^{\rm tar})\,\mathrm{d}\mathcal{H}^{d-1}+\gamma\int_{\Omega}\varepsilon|\nabla{\bm{\varphi}}|^{2}+\frac{1}{\varepsilon}\Psi({\bm{\varphi}})\,\mathrm{dx},

(2.3)

where $\gamma>0$ is a weighting factor, ${\widehat{{\bf u}}}$ is a solution to (2.2) depending on ${\bm{\varphi}}$ (and also on ${\overline{{\bf u}}}$ , a solution to (2.1)), $W\in{\mathbb{R}}^{d\times d}$ is a fixed weighting matrix, $\Gamma^{\rm tar}$ is a subset of the boundary $\widehat{\Gamma}_{N}$ , $\varepsilon>0$ is a fixed constant related to the thickness of the interfacial regions $\{0<\varphi_{i}<1\}$ , $i\in\{1,...,L\}$ , $\mathcal{H}^{d-1}$ indicates the standard $(d-1)$ -dimensional Hausdorff measure, and $\Psi:{\mathbb{R}}^{L}\to{\mathbb{R}}$ is a non-negative multi-well potential that attains its minimum at the corners $\{\bm{e}_{1},\dots,\bm{e}_{L}\}$ (the unit vectors in ${\mathbb{R}}^{L}$ ) of the Gibbs simplex $\Delta^{L}$ . The first term in (2.3) consists of a target shape matching term, where we like to match the displacement ${\widehat{{\bf u}}}$ in Stage 2 with a prescribed deformation ${\bf u}^{\rm tar}$ over the surface $\Gamma^{\rm tar}\subset\widehat{\Gamma}_{N}$ by minimising the squared difference weighted by a matrix $W$ . The second term is the well-known Ginzburg–Landau functional in the multiphase setting that serves as a form of perimeter regularisation. It provides some form of regularity to our design solutions and penalises designs that have large interfaces between the different phases of elastic materials.

For our problem we introduce the design space

\mathcal{U}_{\rm ad}=\Big{\{}{\bm{\varphi}}\in H^{1}(\Omega,{\mathbb{R}}^{L})\,:\,{\bm{\varphi}}(\bm{x})\in\Delta^{L}\text{ for a.e.~{}}\bm{x}\in\Omega\Big{\}},

and our design problem can be formulated as the following

\displaystyle{\bf(P)}\quad

\displaystyle\min_{{\bm{\varphi}}\in\mathcal{U}_{\rm ad}}J({\bm{\varphi}},{\widehat{{\bf u}}})\text{ subject to }({\bm{\varphi}},{\overline{{\bf u}}},{\widehat{{\bf u}}})\text{ is a solution to }\eqref{bu:sys}-\eqref{hu:sys}.

Remark 2.1.

For the existence theory for optimal designs to $(\bf P)$ , it is also possible to consider a more general form of the cost functional:

	$\displaystyle J({\bm{\varphi}},{\overline{{\bf u}}},{\widehat{{\bf u}}})$	$\displaystyle=\int_{\Omega}h_{\Omega}(\bm{x},{\overline{{\bf u}}},{\widehat{{\bf u}}})\,\mathrm{dx}+\int_{\overline{\Gamma}_{N}}\overline{h}({\bm{s}},{\overline{{\bf u}}})\,\mathrm{d}\mathcal{H}^{d-1}+\int_{\widehat{\Gamma}_{N}}\widehat{h}(\bm{s},{\widehat{{\bf u}}})\,\mathrm{d}\mathcal{H}^{d-1}$
		$\displaystyle\quad+\gamma\int_{\Omega}\varepsilon\|\nabla{\bm{\varphi}}\|^{2}+\frac{1}{\varepsilon}{\Psi({\bm{\varphi}})}\,\mathrm{dx},$

with Carathéodory functions $h_{\Omega}$ , $\overline{h}$ and $\widehat{h}$ satisfying (see [14, Remark 5])

$\displaystyle\|h_{\Omega}(\bm{x},\bm{v},\bm{w})\|$	$\displaystyle\leq a_{1}(\bm{x})+b_{1}(\bm{x})\|\bm{v}\|^{p}+b_{2}(\bm{x})\|\bm{w}\|^{p}\quad$	$\displaystyle\text{ for all }\bm{v},\bm{w}\in{\mathbb{R}}^{d},\text{ a.e.~{}}\bm{x}\in\Omega,$
$\displaystyle\|\overline{h}(\bm{s},\bm{v})\|$	$\displaystyle\leq a_{2}(\bm{s})+c_{1}(\bm{s})\|\bm{v}\|^{2}\quad$	$\displaystyle\text{ for all }\bm{v}\in{\mathbb{R}}^{d},\text{ a.e.~{}}\bm{s}\in\overline{\Gamma}_{N},$
$\displaystyle\|\overline{h}(\bm{s},\bm{w})\|$	$\displaystyle\leq a_{3}(\bm{s})+c_{2}(\bm{s})\|\bm{w}\|^{2}\quad$	$\displaystyle\text{ for all }\bm{w}\in{\mathbb{R}}^{d},\text{ a.e.~{}}\bm{s}\in\widehat{\Gamma}_{N},$

for any $2\leq p<\infty$ if $d=2$ and $2\leq p<6$ if $d=3$ , with functions $a_{1}\in L^{1}(\Omega)$ , $b_{1},b_{2}\in L^{\infty}(\Omega)$ , $a_{2}\in L^{1}(\overline{\Gamma}_{N})$ , $a_{3}\in L^{1}(\widehat{\Gamma}_{N})$ , $c_{1}\in L^{\infty}(\overline{\Gamma}_{N})$ and $c_{2}\in L^{\infty}(\widehat{\Gamma}_{N})$ . In our current setting we have $h_{\Omega}=\overline{h}=0$ and $\widehat{h}(\bm{s},{\widehat{{\bf u}}})=\frac{1}{2}\mathrm{\mathcal{X}}_{\Gamma^{\rm tar}}(\bm{s})W({\widehat{{\bf u}}}-\bm{u}^{\rm tar})\cdot({\widehat{{\bf u}}}-\bm{u}^{\rm tar})$ .

Remark 2.2.

It is also possible to consider mass constraints for ${\bm{\varphi}}$ of the form

\displaystyle\frac{1}{|\Omega|}\int_{\Omega}{\bm{\varphi}}\,\mathrm{dx}\leq\bm{\alpha}\quad\text{ or }\quad\frac{1}{|\Omega|}\int_{\Omega}{\bm{\varphi}}\,\mathrm{dx}\geq\bm{\beta}\quad\text{ or }\quad\bm{\beta}\leq\frac{1}{|\Omega|}\int_{\Omega}{\bm{\varphi}}\,\mathrm{dx}\leq\bm{\alpha},

for fixed vectors $\bm{\alpha},\bm{\beta}\in\Delta^{L}$ (possibly also $\bm{\alpha}=\bm{\beta}$ ), where in the above the inequalities are taken component-wise. These are convex constraints and thus when included into the definition of $\mathcal{U}_{\rm ad}$ , the design space remains a closed and convex set. Then, in the corresponding necessary optimality condition, associated Lagrange multipliers will appear, see [13, 14] for more details.

Notation.

For a Banach space $X$ , we denote its topological dual by $X^{*}$ , and the corresponding duality pairing by $\langle\cdot,\cdot\rangle_{X}$ . For any $p\in[1,\infty]$ and $k>0$ , the standard Lebesgue and Sobolev spaces over $\Omega$ are denoted by $L^{p}:=L^{p}(\Omega)$ and $W^{k,p}:=W^{k,p}(\Omega)$ with the corresponding norms $\|\cdot\|_{L^{p}(\Omega)}$ and $\|\cdot\|_{W^{k,p}(\Omega)}$ . In the special case $p=2$ , these become Hilbert spaces and we employ the notation $H^{k}:=H^{k}(\Omega)=W^{k,2}(\Omega)$ with the corresponding norm $\|\cdot\|_{H^{k}(\Omega)}$ . For our subsequent analysis, we introduce the spaces

\displaystyle\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}):=\left\{{\bf v}\in H^{1}(\Omega,{\mathbb{R}}^{d}):{\bf v}=\bm{0}\quad\text{a.e. on $\overline{\Gamma}_{D}$}\right\},

and

\displaystyle\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}):=\left\{{\bf v}\in H^{1}(\Omega,{\mathbb{R}}^{d}):{\bf v}=\bm{0}\quad\text{a.e. on $\widehat{\Gamma}_{D}$}\right\}.

For brevity, the corresponding norms are denoted by the same symbol $\mathopen{\|}\cdot\mathclose{\|}_{H^{1}(\Omega)}$ if no confusion may arise. Vectors, matrices, and vector- or matrix-valued functions will be denoted by bold symbols. Furthermore, for a subset $\Gamma_{N}\subset\Gamma$ , we consider the function space

H^{1/2}_{00}(\Gamma_{N},{\mathbb{R}}^{d}):=\left\{{\bf v}\in H^{1/2}(\Gamma_{N},{\mathbb{R}}^{d}):\widetilde{\bf v}\in H^{1/2}(\Gamma,{\mathbb{R}}^{d})\right\},

where $\widetilde{\bf v}$ denotes the trivial extension of ${\bf v}$ to $\Gamma$ , and we endow it with the norm

\displaystyle\mathopen{\|}{\bf v}\mathclose{\|}_{H^{1/2}_{00}(\Gamma_{N})}:=\mathopen{\|}\widetilde{\bf v}\mathclose{\|}_{H^{1/2}(\Gamma)},\quad{\bf v}\in H^{1/2}_{00}(\Gamma_{N},{\mathbb{R}}^{d}).

We highlight that the above definition is not redundant as in general the trivial extension of a $H^{1/2}(\Gamma_{N},{\mathbb{R}}^{d})$ function does not belong to $H^{1/2}(\Gamma,{\mathbb{R}}^{d})$ . Besides, we remark that $H^{1/2}_{00}(\Gamma_{N},{\mathbb{R}}^{d})$ is a Hilbert space and

\displaystyle(H^{1/2}_{00}(\Gamma_{N},{\mathbb{R}}^{d}),L^{2}(\Gamma_{N},{\mathbb{R}}^{d}),H^{1/2}_{00}(\Gamma_{N},{\mathbb{R}}^{d})^{*})

forms a Hilbert triple (see, e.g., [48]).

For the forthcoming analysis we make the following structural assumptions.

(A1)

The domain $\Omega\subset{\mathbb{R}}^{d}$ , $d\in\{2,3\}$ , is a bounded domain with $C^{1,1}$ or convex boundary $\Gamma=\partial\Omega$ that admits a partition $\Gamma={\rm cl}(\overline{\Gamma}_{D})\cup{\rm cl}(\overline{\Gamma}_{N})$ with relative open subsets $\overline{\Gamma}_{D}$ and $\overline{\Gamma}_{N}$ such that $\overline{\Gamma}_{D}\cap\overline{\Gamma}_{N}=\emptyset$ and $\overline{\Gamma}_{D}\neq\emptyset$ . Here, ${\rm cl}(A)$ denotes the closure of the set $A$ . The same assumptions are made for $\widehat{\Gamma}_{D}$ and $\widehat{\Gamma}_{N}$ .

(A2)

The elasticity tensor ${\overline{{\mathbb{C}}}}$ is assumed to be a tensor-valued function

\displaystyle{\overline{{\mathbb{C}}}}:{\mathbb{R}}^{L}\to{\mathbb{R}}^{d\times d\times d\times d},

with ${\overline{{\mathbb{C}}}}_{ijkl}\in{C^{1,1}({\mathbb{R}}^{L},{\mathbb{R}})}$ , $i,j,k,l\in\{1,\dots,d\}$ . Moreover, it fulfils the symmetry conditions

\displaystyle{\overline{{\mathbb{C}}}}_{ijkl}={\overline{{\mathbb{C}}}}_{klij}={\overline{{\mathbb{C}}}}_{ijlk}={\overline{{\mathbb{C}}}}_{jikl}\quad\text{for all $i,j,k,l\in\{1,\dots,d\}$,}

and there exist positive constants $C_{0},C_{1}$ , and $C_{2}$ such that, for all ${\bm{\varphi}},{\bf h}\in{\mathbb{R}}^{L}$ ,

		$\displaystyle C_{0}\|{{\bf A}}\|^{2}\leq{\overline{{\mathbb{C}}}}({\bm{\varphi}}){{\bf A}}:{{\bf A}}\leq C_{1}\|{{\bf A}}\|^{2}\qquad$	$\displaystyle\forall{{\bf A}}\in{\mathbb{R}}^{d\times d}_{\rm sym}\setminus\{\bm{0}\},$		(2.4)
		$\displaystyle\|{\overline{{\mathbb{C}}}}^{\prime}({\bm{\varphi}}){\bf h}\,{\bf A}:{\bf B}\|\leq C_{2}\|{\bf h}\|\|{\bf A}\|\|{\bf B}\|\qquad$	$\displaystyle\forall{{\bf A},{\bf B}}\in{\mathbb{R}}^{d\times d}_{\rm sym}\setminus\{\bm{0}\},$		(2.5)

where ${\bf A}:{\bf B}=\sum_{i,j=1}^{n}A_{ij}B_{ij}$ , $|{\bf A}|=\sqrt{{\bf A}:{\bf A}}$ , and for every ${\bf h}=(h_{1},\dots,h_{L})\in{\mathbb{R}}^{L}$ ,

\displaystyle[{\overline{{\mathbb{C}}}}^{\prime}({\bm{\varphi}}){\bf h}]_{ijkl}:=\sum_{m=1}^{L}\partial_{m}{\overline{{\mathbb{C}}}}_{ijkl}({\bm{\varphi}})h_{m}\quad\text{ for all }i,j,k,l\in\{1,\dots,d\}.

The set ${\mathbb{R}}^{d\times d}_{\rm sym}$ consists of the symmetric $({d\times d})$ -matrices. The same assumptions are made for $\widehat{{\mathbb{C}}}$ .

(A3)

The multiwell potential $\Psi$ possesses the form

$\displaystyle\Psi:{\mathbb{R}}^{L}\to{\mathbb{R}}\cup\{+\infty\},\quad\Psi=\widetilde{\Psi}+I_{\Delta},$

where $\widetilde{\Psi}\in{C^{1,1}({\mathbb{R}}^{L})}$ and the indicator function $I_{\Delta}$ of the simplex $\Delta^{L}$ is defined as

$\displaystyle I_{\Delta}({\bm{\varphi}})=\begin{cases}0\quad&\text{if ${\bm{\varphi}}\in\Delta^{L}$},\\ +\infty\quad&\text{otherwise}.\end{cases}$
(A4)

The function $\chi:{\mathbb{R}}^{L}\to{\mathbb{R}}$ is ${C^{1,1}({\mathbb{R}}^{L})}$ and there exist a positive constant $\chi_{0}$ such that

$\displaystyle 0\leq\chi({\bm{\varphi}})\leq\chi_{0}\qquad\text{ for all }{\bm{\varphi}}\in{\mathbb{R}}^{L}.$

(A5)

The data of the problems satisfy

	$\displaystyle\overline{{\bf F}},\widehat{{\bf F}}$	$\displaystyle\in L^{2}(\Omega,{\mathbb{R}}^{d}),\quad\overline{{\bf U}}\in H^{1/2}(\overline{\Gamma}_{D},{\mathbb{R}}^{d}),\quad\widehat{{\bf U}}\in H^{1/2}(\widehat{\Gamma}_{D},{\mathbb{R}}^{d}),$
	$\displaystyle\overline{{\bf g}}$	$\displaystyle\in H^{1/2}_{00}(\overline{\Gamma}_{N},{\mathbb{R}}^{d})^{},\quad\widehat{{\bf g}}\in H^{1/2}_{00}(\widehat{\Gamma}_{N},{\mathbb{R}}^{d})^{}.$

(A6)

The target displacement ${\bf u}^{\rm tar}\in L^{2}(\Gamma^{\rm tar},{\mathbb{R}}^{d})$ , where $\Gamma^{\rm tar}\subset\widehat{\Gamma}_{N}$ .

It is worth pointing out that condition (A3) entails that $\Psi({\bm{\varphi}})=\widetilde{\Psi}({\bm{\varphi}})$ for every ${\bm{\varphi}}\in\mathcal{U}_{\rm ad}$ .

3 Analysis of the design optimisation problem

3.1 Linear elasticity system with mixed boundary conditions

In this section we provide a preliminary well-posedness result for the following linear elasticity system with mixed boundary conditions


$\displaystyle-\mathop{\rm div}\nolimits\big{(}{\mathbb{C}}\mathcal{E}({\bf u})+{\mathbb{F}}\big{)}$	$\displaystyle={\bf f}\qquad$	$\displaystyle\text{in $\Omega$},$	(3.1a)
$\displaystyle{\bf u}$	$\displaystyle=\bm{0}\qquad$	$\displaystyle\text{on $\Gamma_{D}$},$	(3.1b)
$\displaystyle({\mathbb{C}}\mathcal{E}({\bf u})+{\mathbb{F}}){\bf n}$	$\displaystyle={\bf g}\qquad$	$\displaystyle\text{on $\Gamma_{N}$}.$	(3.1c)

The well-posedness of the system (3.1) is formulated as follows.

Proposition 3.1.

In addition to (A1), suppose that ${\mathbb{C}}\in L^{\infty}(\Omega,{\mathbb{R}}^{d\times d\times d\times d})$ fulfils the symmetry conditions

\displaystyle{\mathbb{C}}_{ijkl}

\displaystyle={\mathbb{C}}_{klij}={\mathbb{C}}_{ijlk}={\mathbb{C}}_{jikl}\quad

for all

i,j,k,l\in\{1,\dots,d\}

and satisfies

\displaystyle\lambda|{{\bf A}}|^{2}\leq({\mathbb{C}}{{\bf A}}):{{\bf A}}\leq\Lambda|{{\bf A}}|^{2}\quad\text{ for all }{\bf A}\in{\mathbb{R}}^{d\times d}_{\rm sym}\setminus\{\bm{0}\}\,\,\text{and a.e. in $\Omega$},

(3.2)

for some positive constants $\lambda$ and $\Lambda$ . Then, for every ${\mathbb{F}}\in L^{2}(\Omega,{\mathbb{R}}^{d\times d})$ , ${\bf f}\in L^{2}(\Omega,{\mathbb{R}}^{d})$ and ${\bf g}\in H^{1/2}_{00}(\Gamma_{N},{\mathbb{R}}^{d})^{*}$ , there exists a unique weak solution ${\bf u}\in H^{1}_{D}(\Omega,{\mathbb{R}}^{d}):=\{\bm{v}\in H^{1}(\Omega,{\mathbb{R}}^{d})\,:\,\bm{v}=\bm{0}\text{ a.e.~{}in }\Gamma_{D}\}$ to the elasticity system (3.1) in the sense that

\displaystyle\int_{\Omega}{\mathbb{C}}\mathcal{E}({\bf u}):\mathcal{E}({\bf v})\,\mathrm{dx}=-\int_{\Omega}({\mathbb{F}}:\mathcal{E}({\bf v})-{\bf f}\cdot{\bf v})\,\mathrm{dx}+\mathopen{\langle}{\bf g},{\bf v}\mathclose{\rangle}_{H^{1/2}_{00}(\Gamma_{N})}\quad\forall{\bf v}\in H^{1}_{D}(\Omega,{\mathbb{R}}^{d}),

(3.3)

and there exists a positive constant $C$ independent of ${\bf u}$ such that

\displaystyle\mathopen{\|}{\bf u}\mathclose{\|}_{H^{1}(\Omega)}\leq C\left(\mathopen{\|}{\mathbb{F}}\mathclose{\|}_{L^{2}(\Omega)}+\mathopen{\|}{\bf f}\mathclose{\|}_{L^{2}(\Omega)}+\mathopen{\|}{\bf g}\mathclose{\|}_{H^{1/2}_{00}(\Gamma_{N})^{*}}\right).

(3.4)

Note that whenever ${\bf g}\in L^{2}(\Gamma_{N},{\mathbb{R}}^{d})\subset H^{1/2}_{00}(\Gamma_{N},{\mathbb{R}}^{d})^{*}$ , we can identify the duality product in (3.3) as the standard boundary integral, that is,

\displaystyle\mathopen{\langle}{\bf g},{\bf v}\mathclose{\rangle}_{H^{1/2}_{00}(\Gamma_{N})}=\int_{\Gamma_{N}}{\bf g}\cdot{\bf v}\,\mathrm{d}\mathcal{H}^{d-1}.

Proof of Proposition 3.1.

The variational equality (3.3) admits a unique solution by a direct application of the Lax–Milgram theorem (cf., e.g., [5]). In this direction, we set

	$\displaystyle{\bf V}=H^{1}_{D}(\Omega,{\mathbb{R}}^{d})$
	$\displaystyle a(\cdot,\cdot):{\bf V}\times{\bf V}\to{\mathbb{R}},\quad a({\bf u},{\bf v}):=\int_{\Omega}{\mathbb{C}}\mathcal{E}({\bf u}):\mathcal{E}({\bf v})\,\mathrm{dx},$
	$\displaystyle\langle{\bf F},{\bf v}\rangle:=-\int_{\Omega}({\mathbb{F}}:\mathcal{E}({\bf v})-{\bf f}\cdot{\bf v})\,\mathrm{dx}+\mathopen{\langle}{\bf g},{\bf v}\mathclose{\rangle}_{H^{1/2}_{00}(\Gamma_{N})}.$

It is worth noticing that it readily follows from the assumptions on ${\mathbb{F}}$ , ${\bf f}$ , and ${\bf g}$ that ${\bf F}\in{\bf V}^{*}$ . With the above notation, (3.3) can then be rewritten as the variational problem

\displaystyle a({\bf u},{\bf v})=\mathopen{\langle}{\bf F},{\bf v}\mathclose{\rangle}\quad\forall{\bf v}\in{\bf V}.

Thus, to apply the Lax–Milgram theorem, it is sufficient to show the bilinear form $a(\cdot,\cdot)$ is continuous and coercive in ${\bf V}$ . By (3.2) we have

\displaystyle|a({\bf u},{\bf v})|\leq C\mathopen{\|}{\bf u}\mathclose{\|}_{H^{1}(\Omega)}\mathopen{\|}{\bf v}\mathclose{\|}_{H^{1}(\Omega)}\leq C\mathopen{\|}{\bf u}\mathclose{\|}_{{\bf V}}\mathopen{\|}{\bf v}\mathclose{\|}_{{\bf V}}\quad\forall{\bf u},{\bf v}\in{{\bf V}},

while (3.2) and Korn’s inequality yield the ${\bf V}$ -coercivity:

\displaystyle|a({\bf v},{\bf v})|\geq\lambda\mathopen{\|}\mathcal{E}({\bf v})\mathclose{\|}_{L^{2}(\Omega)}^{2}\geq C(\lambda,{C_{K}})\mathopen{\|}{\bf v}\mathclose{\|}_{{\bf V}}^{2}\quad\forall{\bf v}\in{{\bf V}},

with a constant ${C_{K}}$ arising from Korn’s inequality. Thus, the existence and uniqueness of ${\bf u}\in H^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ solving (3.3), as well as the estimate (3.4), readily follow from the Lax–Milgram theorem. ∎

We end this section with another abstract result that will be useful for the subsequent analysis. Consider the following problem with inhomogeneous data on the Dirichlet boundary:

$\displaystyle-\mathop{\rm div}\nolimits\big{(}{\mathbb{C}}\mathcal{E}({\bf u})\big{)}$	$\displaystyle=\bm{0}\qquad$	$\displaystyle\text{ in }\Omega,$	(3.5)
$\displaystyle{\bf u}$	$\displaystyle={\cal U}\qquad$	$\displaystyle\text{ on }\Gamma_{D},$
$\displaystyle\big{(}{\mathbb{C}}\mathcal{E}({\bf u})\big{)}{\bf n}$	$\displaystyle=\bm{0}\qquad$	$\displaystyle\text{ on }\Gamma_{N}.$

Well-known theory yields that, for every ${\cal U}\in H^{1/2}(\Gamma_{D},{\mathbb{R}}^{d})$ , there exists a unique weak solution ${\bf u}\in H^{1}(\Omega,{\mathbb{R}}^{d})$ . The proof follows similarly to the above as a direct consequence of the Lax–Milgram theorem. This allows us to introduce the associated solution operator, that we call the extension operator

\displaystyle{\cal H}:H^{1/2}(\Gamma_{D},{\mathbb{R}}^{d})\to H^{1}(\Omega,{\mathbb{R}}^{d}),\quad{\cal H}:{\cal U}\mapsto{\bf u},\quad

(3.6)

where ${\bf u}$ is the unique weak solution to system (3.5).

3.2 Well-posedness of the state systems

Similarly to (3.5) and (3.6), we can introduce the extension operators $\overline{\cal H}$ and $\widehat{\cal H}$ related to $\overline{{\mathbb{C}}}$ , $\overline{\Gamma}_{D}$ and $\widehat{{\mathbb{C}}}$ , $\widehat{\Gamma}_{D}$ , respectively. Then defining the functions $\overline{H}:=\overline{\cal H}(\overline{\bm{U}}),\widehat{H}:=\widehat{\cal H}(\widehat{\bm{U}})\in H^{1}(\Omega,{\mathbb{R}}^{d})$ allows us to transform (2.1) and (2.2) into the equivalent problems

$\displaystyle-\mathop{\rm div}\nolimits\big{(}\overline{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}}^{\rm new}+{\overline{H}})\big{)}$	$\displaystyle=\overline{{\bf F}}\qquad$	$\displaystyle\text{ in }\Omega,$	(3.7)
$\displaystyle{\overline{{\bf u}}}^{\rm new}$	$\displaystyle=\bm{0}\qquad$	$\displaystyle\text{ on }\overline{\Gamma}_{D},$
$\displaystyle\big{(}\overline{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}}^{\rm new}+{\overline{H}})\big{)}{\bf n}$	$\displaystyle=\overline{{\bf g}}\qquad$	$\displaystyle\text{ on }\overline{\Gamma}_{N},$

and

$\displaystyle-\mathop{\rm div}\nolimits\big{(}\widehat{{\mathbb{C}}}({\bm{\varphi}})(\mathcal{E}({\widehat{{\bf u}}}^{\rm new}+{\widehat{H}})-\chi({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}}^{\rm new}+{\overline{H}})\big{)}$	$\displaystyle=\widehat{{\bf F}}\qquad$	$\displaystyle\text{ in }\Omega,$	(3.8)
$\displaystyle{\widehat{{\bf u}}}^{\rm new}$	$\displaystyle=\bm{0}\qquad$	$\displaystyle\text{ on }\widehat{\Gamma}_{D},$
$\displaystyle\big{(}\widehat{{\mathbb{C}}}({\bm{\varphi}})(\mathcal{E}({\widehat{{\bf u}}}^{\rm new}+{\widehat{H}})-\chi({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}}^{\rm new}+{\overline{H}})\big{)}{\bf n}$	$\displaystyle=\widehat{{\bf g}}\qquad$	$\displaystyle\text{ on }\widehat{\Gamma}_{N},$

where we set

{\overline{{\bf u}}}={\overline{{\bf u}}}^{\rm new}+{\overline{H}},\quad{\widehat{{\bf u}}}={\widehat{{\bf u}}}^{\rm new}+{\widehat{H}}.

For a cleaner presentation, we abuse notation and use the same variables ${\overline{{\bf u}}}$ and ${\widehat{{\bf u}}}$ to denote ${\overline{{\bf u}}}^{\rm new}$ and ${\widehat{{\bf u}}}^{\rm new}$ . Then, the well-posedness of (2.1) and (2.2) (equivalently (3.7) and (3.8)) are formulated as follows.

Theorem 3.1.

Under (A1)–(A5), for every ${\bm{\varphi}}\in L^{\infty}(\Omega,{\mathbb{R}}^{L})$ , there exists a unique solution pair $({\overline{{\bf u}}},{\widehat{{\bf u}}})\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})\times\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ satisfying

		$\displaystyle\int_{\Omega}\overline{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}}):\mathcal{E}(\overline{{\bm{\zeta}}})\,\mathrm{dx}=\int_{\Omega}\overline{{\bf F}}\cdot\overline{{\bm{\zeta}}}\,\mathrm{dx}-\int_{\Omega}\overline{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}({\overline{H}}):\mathcal{E}(\overline{{\bm{\zeta}}})\,\mathrm{dx}+\mathopen{\langle}\overline{{\bf g}},\overline{{\bm{\zeta}}}\mathclose{\rangle}_{H^{1/2}_{00}(\overline{\Gamma}_{N})},$		(3.9)
		$\displaystyle\int_{\Omega}\widehat{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}({\widehat{{\bf u}}}):\mathcal{E}(\widehat{{\bm{\zeta}}})\,\mathrm{dx}-\int_{\Omega}\widehat{{\mathbb{C}}}({\bm{\varphi}})\chi({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}}):\mathcal{E}(\widehat{{\bm{\zeta}}})\,\mathrm{dx}=\int_{\Omega}\widehat{{\bf F}}\cdot\widehat{{\bm{\zeta}}}\,\mathrm{dx}$
		$\displaystyle\quad+\int_{\Omega}\widehat{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}({\widehat{H}}):\mathcal{E}(\widehat{{\bm{\zeta}}})\,\mathrm{dx}-\int_{\Omega}\widehat{{\mathbb{C}}}({\bm{\varphi}})\chi({\bm{\varphi}})\mathcal{E}({\overline{H}}):\mathcal{E}(\widehat{{\bm{\zeta}}})\,\mathrm{dx}+\mathopen{\langle}\widehat{{\bf g}},\widehat{{\bm{\zeta}}}\mathclose{\rangle}_{H^{1/2}_{00}(\widehat{\Gamma}_{N})}$		(3.10)

for all $\overline{{\bm{\zeta}}}\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ and $\widehat{{\bm{\zeta}}}\in\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ . Moreover, there exists a positive constant $C$ , independent of ${\bm{\varphi}}$ , such that

\displaystyle\mathopen{\|}{\overline{{\bf u}}}\mathclose{\|}_{H^{1}(\Omega)}+\mathopen{\|}{\widehat{{\bf u}}}\mathclose{\|}_{H^{1}(\Omega)}\leq C.

(3.11)

Proof.

For (3.9) we invoke Proposition 3.1 with the specifications

\displaystyle{\bf u}={\overline{{\bf u}}},\quad{\mathbb{C}}=\overline{{\mathbb{C}}}({\bm{\varphi}}),\quad{\mathbb{F}}=\overline{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}({\overline{H}}),\quad{\bf f}=\overline{{\bf F}},\quad{\bf g}=\overline{{\bf g}},\quad\Gamma_{D}=\overline{\Gamma}_{D},\quad\Gamma_{N}=\overline{\Gamma}_{N},

and for (3.10) we consider

	$\displaystyle{\bf u}={\widehat{{\bf u}}},\quad{\mathbb{C}}=\widehat{{\mathbb{C}}}({\bm{\varphi}}),\quad{\mathbb{F}}=\widehat{{\mathbb{C}}}({\bm{\varphi}})(\mathcal{E}({\widehat{H}})-\chi({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}}+{\overline{H}})),\quad{\bf f}=\widehat{{\bf F}},\quad{\bf g}=\widehat{{\bf g}},$
	$\displaystyle\Gamma_{D}=\widehat{\Gamma}_{D},\quad\Gamma_{N}=\widehat{\Gamma}_{N},$

to obtain the existence and uniqueness of solutions ${\overline{{\bf u}}}$ and ${\widehat{{\bf u}}}$ . Lastly the estimate (3.11) can be obtained from (3.4) and the uniform boundedness of the tensors $\overline{{\mathbb{C}}}$ and $\widehat{{\mathbb{C}}}$ in (A2). ∎

Theorem 3.2.

Under (A1)–(A5), for $i=1,2$ , let ${\bm{\varphi}}_{i}\in L^{\infty}(\Omega,{\mathbb{R}}^{L})$ with $\mathopen{\|}{\bm{\varphi}}_{i}\mathclose{\|}_{L^{\infty}(\Omega)}\leq R$ , where $R>0$ is fixed, and let $({\overline{{\bf u}}}_{i},{\widehat{{\bf u}}}_{i})$ denote the unique solutions to systems (3.7) and (3.8) corresponding to ${\bm{\varphi}}_{i}$ , but with the same data ${\overline{{\bf F}}}$ , $\widehat{{\bf F}}$ , $\overline{{\bf g}}$ , $\widehat{{\bf g}}$ , ${\overline{H}}$ and ${\widehat{H}}$ . Then, there exists a positive constant $C$ , independent of the differences, such that

	$\displaystyle\mathopen{\\|}{\overline{{\bf u}}}_{1}-{\overline{{\bf u}}}_{2}\mathclose{\\|}_{H^{1}(\Omega)}$	$\displaystyle\leq C\mathopen{\\|}{\bm{\varphi}}_{1}-{\bm{\varphi}}_{2}\mathclose{\\|}_{L^{\infty}(\Omega)},$		(3.12)
	$\displaystyle\mathopen{\\|}{\widehat{{\bf u}}}_{1}-{\widehat{{\bf u}}}_{2}\mathclose{\\|}_{H^{1}(\Omega)}$	$\displaystyle\leq C\mathopen{\\|}\mathcal{E}({\overline{{\bf u}}}_{1}-{\overline{{\bf u}}}_{2})\mathclose{\\|}_{L^{2}(\Omega)}+C\mathopen{\\|}{\bm{\varphi}}_{1}-{\bm{\varphi}}_{2}\mathclose{\\|}_{L^{\infty}(\Omega)}.$		(3.12)

Proof.

To start, let us set

\displaystyle{\bm{\varphi}}:={\bm{\varphi}}_{1}-{\bm{\varphi}}_{2},\quad{\overline{{\bf u}}}:={\overline{{\bf u}}}_{1}-{\overline{{\bf u}}}_{2},\quad{\widehat{{\bf u}}}:={\widehat{{\bf u}}}_{1}-{\widehat{{\bf u}}}_{2}.

Then, we consider the difference between the variational equalities (3.9)–(3.10) written for $({\bm{\varphi}}_{1},{\overline{{\bf u}}}_{1},{\widehat{{\bf u}}}_{1})$ and for $({\bm{\varphi}}_{2},{\overline{{\bf u}}}_{2},{\widehat{{\bf u}}}_{2})$ to infer that

		$\displaystyle\int_{\Omega}\big{(}\overline{{\mathbb{C}}}({\bm{\varphi}}_{1})\mathcal{E}({\overline{{\bf u}}}_{1})-\overline{{\mathbb{C}}}({\bm{\varphi}}_{2})\mathcal{E}({\overline{{\bf u}}}_{2})\big{)}:\mathcal{E}(\overline{{\bm{\zeta}}})\,\mathrm{dx}=0,$		(3.13)
		$\displaystyle\int_{\Omega}\big{(}\widehat{{\mathbb{C}}}({\bm{\varphi}}_{1})\mathcal{E}({\widehat{{\bf u}}}_{1})-\widehat{{\mathbb{C}}}({\bm{\varphi}}_{2})\mathcal{E}({\widehat{{\bf u}}}_{2})\big{)}:\mathcal{E}(\widehat{{\bm{\zeta}}})\,\mathrm{dx}$
		$\displaystyle\quad-\int_{\Omega}\big{(}\widehat{{\mathbb{C}}}({\bm{\varphi}}_{1})\chi({\bm{\varphi}}_{1})\mathcal{E}({\overline{{\bf u}}}_{1})-\widehat{{\mathbb{C}}}({\bm{\varphi}}_{2})\chi({\bm{\varphi}}_{2})\mathcal{E}({\overline{{\bf u}}}_{2})\big{)}:\mathcal{E}(\widehat{{\bm{\zeta}}})\,\mathrm{dx}=0,$		(3.14)

for all $\overline{{\bm{\zeta}}}\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ and $\widehat{{\bm{\zeta}}}\in\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}).$ Choosing $\overline{{\bm{\zeta}}}={\overline{{\bf u}}}$ and invoking condition (2.4) yields

	$\displaystyle C_{0}\mathopen{\\|}\mathcal{E}({\overline{{\bf u}}})\mathclose{\\|}_{L^{2}(\Omega)}^{2}$	$\displaystyle\leq-\int_{\Omega}(\overline{{\mathbb{C}}}({\bm{\varphi}}_{1})-\overline{{\mathbb{C}}}({\bm{\varphi}}_{2}))\mathcal{E}({\overline{{\bf u}}}_{2}):\mathcal{E}({\overline{{\bf u}}})\,\mathrm{dx}$		(3.15)
		$\displaystyle\leq\frac{C_{0}}{2}\mathopen{\\|}\mathcal{E}({\overline{{\bf u}}})\mathclose{\\|}_{L^{2}(\Omega)}^{2}+C\mathopen{\\|}{\bm{\varphi}}\mathclose{\\|}_{L^{\infty}(\Omega)}^{2}.$		(3.15)

By Korn’s inequality we infer

\displaystyle\mathopen{\|}{\overline{{\bf u}}}_{1}-{\overline{{\bf u}}}_{2}\mathclose{\|}_{H^{1}(\Omega)}\leq C\mathopen{\|}{\bm{\varphi}}_{1}-{\bm{\varphi}}_{2}\mathclose{\|}_{L^{\infty}(\Omega)}.

(3.16)

Then, inserting $\widehat{{\bm{\zeta}}}={\widehat{{\bf u}}}$ in (3.14), and invoking the Lipschitz continuity of $\widehat{{\mathbb{C}}}$ and $\chi$ from (A2) and (A4), as well as (3.16), yields

		$\displaystyle C_{0}\mathopen{\\|}\mathcal{E}({\widehat{{\bf u}}})\mathclose{\\|}_{L^{2}(\Omega)}^{2}$		(3.17)
		$\displaystyle\quad\leq-\int_{\Omega}(\widehat{{\mathbb{C}}}({\bm{\varphi}}_{1})-\widehat{{\mathbb{C}}}({\bm{\varphi}}_{2}))\mathcal{E}({\widehat{{\bf u}}}_{2}):\mathcal{E}({\widehat{{\bf u}}})+(\widehat{{\mathbb{C}}}({\bm{\varphi}}_{1})-\widehat{{\mathbb{C}}}({\bm{\varphi}}_{2}))\chi({\bm{\varphi}}_{1})\mathcal{E}({\overline{{\bf u}}}_{1}):\mathcal{E}({\widehat{{\bf u}}})\,\mathrm{dx}$
		$\displaystyle\qquad-\int_{\Omega}\widehat{{\mathbb{C}}}({\bm{\varphi}}_{2})(\chi({\bm{\varphi}}_{1})-\chi({\bm{\varphi}}_{2}))\mathcal{E}({\overline{{\bf u}}}_{1}):\mathcal{E}({\widehat{{\bf u}}})\,\mathrm{dx}-\int_{\Omega}\widehat{{\mathbb{C}}}({\bm{\varphi}}_{2})\chi({\bm{\varphi}}_{2})\mathcal{E}({\overline{{\bf u}}}):\mathcal{E}({\widehat{{\bf u}}})\,\mathrm{dx}$
		$\displaystyle\quad\leq\frac{C_{0}}{2}\mathopen{\\|}\mathcal{E}({\widehat{{\bf u}}})\mathclose{\\|}_{L^{2}(\Omega)}^{2}+C\Big{(}\mathopen{\\|}\mathcal{E}({\overline{{\bf u}}})\mathclose{\\|}_{L^{2}(\Omega)}^{2}+\mathopen{\\|}{\bm{\varphi}}\mathclose{\\|}_{L^{\infty}(\Omega)}^{2}\Big{)}.$

Applying Korn’s inequality leads to (3.12). ∎

Corollary 3.1.

Suppose that (A1)–(A5) hold. Let $R>0$ and let $\{{\bm{\varphi}}_{n}\}_{n\in{\mathbb{N}}}$ be a sequence of functions in $L^{\infty}(\Omega,{\mathbb{R}}^{L})$ such that $\|{\bm{\varphi}}_{n}\|_{L^{\infty}(\Omega)}\leq R$ for all $n\in{\mathbb{N}}$ and ${\bm{\varphi}}_{n}\to{\bm{\varphi}}_{*}$ strongly in $L^{1}(\Omega,{\mathbb{R}}^{L})$ as $n\to\infty$ . Let $({\overline{{\bf u}}}_{n},{\widehat{{\bf u}}}_{n})$ denote the solutions to (3.9) and (3.10) corresponding to data ${\bm{\varphi}}_{n}$ , $\overline{{\bf F}}$ , $\widehat{{\bf F}}$ , $\overline{{\bf g}}$ , $\widehat{{\bf g}}$ , ${\overline{H}}$ and ${\widehat{H}}$ . Then, it holds that, as $n\to\infty$ ,

{\overline{{\bf u}}}_{n}\to{\overline{{\bf u}}}_{*}\text{ strongly in }\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}),\quad{\widehat{{\bf u}}}_{n}\to{\widehat{{\bf u}}}_{*}\text{ strongly in }\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}),

where $({\overline{{\bf u}}}_{*},{\widehat{{\bf u}}}_{*})\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})\times\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ are the unique solutions to (3.9) and (3.10) corresponding to data ${\bm{\varphi}}_{*}$ , $\overline{{\bf F}}$ , $\widehat{{\bf F}}$ , $\overline{{\bf g}}$ , $\widehat{{\bf g}}$ , ${\overline{H}}$ and ${\widehat{H}}$ .

Proof.

From (3.11) we infer that ${\overline{{\bf u}}}_{n}$ and ${\widehat{{\bf u}}}_{n}$ are bounded in $\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ and $\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ , respectively, and thus there exist limit functions $({\overline{{\bf u}}}_{*},{\widehat{{\bf u}}}_{*})\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})\times\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ such that, along a non-relabelled subsequence, ${\overline{{\bf u}}}_{n}\rightharpoonup{\overline{{\bf u}}}_{*}$ in $\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ and ${\widehat{{\bf u}}}_{n}\rightharpoonup{\widehat{{\bf u}}}_{*}$ in $\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ . To obtain strong convergence, in (3.13) and (3.14) we substitute ${\bm{\varphi}}_{1}={\bm{\varphi}}_{n}$ , ${\bm{\varphi}}_{2}={\bm{\varphi}}_{*}$ , ${\overline{{\bf u}}}_{1}={\overline{{\bf u}}}_{n}$ , ${\overline{{\bf u}}}_{2}={\overline{{\bf u}}}_{*}$ , ${\overline{{\bf u}}}={\overline{{\bf u}}}_{n}-{\overline{{\bf u}}}_{*}$ and ${\widehat{{\bf u}}}={\widehat{{\bf u}}}_{n}-{\widehat{{\bf u}}}_{*}$ . Then, by virtue of the dominated convergence theorem, we infer the strong convergence, as $n\to\infty$ ,

\displaystyle({\overline{{\mathbb{C}}}}({\bm{\varphi}}_{n})-{\overline{{\mathbb{C}}}}({\bm{\varphi}}_{*}))\mathcal{E}({\overline{{\bf u}}}_{*})\to\bm{0}\quad\text{ in }L^{2}(\Omega,{\mathbb{R}}^{d\times d}).

Hence, in the analogue of the first inequality in (3.15) we see that the integral on the right-hand side converges to zero, which implies via Korn’s inequality that

\mathopen{\|}{\overline{{\bf u}}}_{n}-{\overline{{\bf u}}}_{*}\mathclose{\|}_{H^{1}(\Omega)}\to 0.

By the generalised dominated convergence theorem, we have the strong convergences

\displaystyle\begin{cases}({\widehat{{\mathbb{C}}}}({\bm{\varphi}}_{n})-{\widehat{{\mathbb{C}}}}({\bm{\varphi}}_{*}))\chi({\bm{\varphi}}_{n})\mathcal{E}({\overline{{\bf u}}}_{n})&\to\bm{0}\\ {\widehat{{\mathbb{C}}}}({\bm{\varphi}}_{*})(\chi({\bm{\varphi}}_{n})-\chi({\bm{\varphi}}_{*}))\mathcal{E}({\widehat{{\bf u}}}_{n})&\to\bm{0}\\ \widehat{{\mathbb{C}}}({\bm{\varphi}}_{*})\chi({\bm{\varphi}}_{*})\mathcal{E}({\overline{{\bf u}}}_{n}-{\overline{{\bf u}}}_{*})&\to\bm{0}\end{cases}\text{ in }L^{2}(\Omega,{\mathbb{R}}^{d\times d}),

and thus, in the analogue of the first inequality in (3.17) we see that the integral on the right-hand side converges to zero, leading to the assertion

\mathopen{\|}{\widehat{{\bf u}}}_{n}-{\widehat{{\bf u}}}_{*}\mathclose{\|}_{H^{1}(\Omega)}\to 0.

Thus, by combining the weak convergences with the above norms convergence the claim follows. ∎

The above analysis for systems (3.7) and (3.8) allows us to define some solution operators. Namely, we introduce

\displaystyle\mathcal{S}:L^{\infty}(\Omega,{\mathbb{R}}^{L})\to\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}),\quad\mathcal{S}:{\bm{\varphi}}\mapsto{\widehat{{\bf u}}}={\widehat{{\bf u}}}({\bm{\varphi}}),

(3.18)

as well as the intermediate operators:

	$\displaystyle\mathcal{S}_{1}:L^{\infty}(\Omega,{\mathbb{R}}^{L})\to L^{\infty}(\Omega,{\mathbb{R}}^{L})\times\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}),\quad$	$\displaystyle\mathcal{S}_{1}:{\bm{\varphi}}\mapsto(\mathcal{S}_{1}^{1}({\bm{\varphi}}),\mathcal{S}_{1}^{2}({\bm{\varphi}}))=({\bm{\varphi}},{{\overline{{\bf u}}}}),$
	$\displaystyle\mathcal{S}_{2}:L^{\infty}(\Omega,{\mathbb{R}}^{L})\times\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})\to\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}),\quad$	$\displaystyle\mathcal{S}_{2}:({\bm{\varphi}},{\overline{{\bf u}}})\mapsto{\widehat{{\bf u}}},$

where ${\overline{{\bf u}}}={\overline{{\bf u}}}({\bm{\varphi}})$ and ${\widehat{{\bf u}}}={\widehat{{\bf u}}}({\bm{\varphi}},{\overline{{\bf u}}})$ are the unique solutions obtained from Theorem 3.1. Then, the overall solution operator $\mathcal{S}$ in (3.18) is simply the composition of the intermediate mappings $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ , i.e., $\mathcal{S}=\mathcal{S}_{2}\circ\mathcal{S}_{1}$ . In particular, we can define the reduced cost functional

J_{\rm red}:\mathcal{U}_{\rm ad}\to{\mathbb{R}},\quad J_{\rm red}:{\bm{\varphi}}\mapsto J({\bm{\varphi}},\mathcal{S}({\bm{\varphi}})).

3.3 Existence of optimal designs

Theorem 3.3.

Under (A1)–(A6), the optimisation problem (P) admits at least one solution.

Proof.

As the proof is nowadays standard with the direct method of the calculus of variations, let us briefly outline the main points. Consider a minimising sequence $\{{\bm{\varphi}}_{n}\}_{n\in{\mathbb{N}}}\subset\mathcal{U}_{\rm ad}$ for the reduced cost functional $J_{\rm red}$ , which satisfies

\displaystyle\lim_{n\to\infty}J_{\rm red}({\bm{\varphi}}_{n})=\inf_{{\bm{\varphi}}\in\mathcal{U}_{\rm ad}}J_{\rm red}({\bm{\varphi}})\geq 0.

This yields that $\{{\bm{\varphi}}_{n}\}_{n\in{\mathbb{N}}}$ is bounded in $H^{1}(\Omega,{\mathbb{R}}^{L})\cap L^{\infty}(\Omega,{\mathbb{R}}^{L})$ . By standard compactness arguments, since $\mathcal{U}_{\rm ad}$ is closed and convex, we obtain a limit function ${\bm{\varphi}}^{*}\in\mathcal{U}_{\rm ad}$ such that ${\bm{\varphi}}_{n}\to{\bm{\varphi}}^{*}$ weakly* in $H^{1}(\Omega,{\mathbb{R}}^{L})\cap L^{\infty}(\Omega,{\mathbb{R}}^{L})$ along a non-relabelled subsequence. Consequently, by (3.11) the sequence $\{{\widehat{{\bf u}}}_{n}=\mathcal{S}({\bm{\varphi}}_{n})\}_{n\in{\mathbb{N}}}$ is bounded in $\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ , and on invoking Corollary 3.1 there exists a limit function ${\widehat{{\bf u}}}^{*}\in\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ such that, along a non-relabelled subsequence, ${\widehat{{\bf u}}}_{n}\to{\widehat{{\bf u}}}^{*}$ strongly in $\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ as $n\to\infty$ . Continuity of the boundary trace operator gives ${\widehat{{\bf u}}}_{n}\to{\widehat{{\bf u}}}^{*}$ strongly in $L^{2}(\widehat{\Gamma}_{N},{\mathbb{R}}^{d})$ , and thus

\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}_{n}-\bm{u}^{\rm tar})\cdot({\widehat{{\bf u}}}_{n}-\bm{u}^{\rm tar})\,\mathrm{d}\mathcal{H}^{d-1}\to\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}^{*}-\bm{u}^{\rm tar})\cdot({\widehat{{\bf u}}}^{*}-\bm{u}^{\rm tar})\,\mathrm{d}\mathcal{H}^{d-1}.

By Fatou’s lemma and the a.e. convergence of ${\bm{\varphi}}_{n}$ to ${\bm{\varphi}}^{*}$ , we have

\displaystyle\liminf_{n\to\infty}\|\Psi({\bm{\varphi}}_{n})\|_{L^{1}(\Omega)}\geq\|\Psi({\bm{\varphi}}^{*})\|_{L^{1}(\Omega)},

and using also the weak lower semicontinuity of the $L^{2}$ -norm, we infer that

J_{\rm red}({\bm{\varphi}}^{*})\leq\lim_{n\to\infty}J_{\rm red}({\bm{\varphi}}_{n})=\inf_{{\bm{\varphi}}\in\mathcal{U}_{\rm ad}}J_{\rm red}({\bm{\varphi}}).

This shows that ${\bm{\varphi}}^{*}$ is a solution to $(\bf P)$ . ∎

4 Optimality conditions

To derive the first order necessary optimality conditions for ${\bm{\varphi}}^{*}$ , we first study the linearised system for the linearised variables introduced below, and use adjoint variables to provide a simplification of the optimality condition.

4.1 Linearised systems and Fréchet differentiability

Here, we analyse some differentiability properties of the solutions operators $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ introduced above. This will help us to formulate the first order optimality conditions of (P).

Theorem 4.1.

The solution operator $\mathcal{S}_{1}$ is Fréchet differentiable at ${\bm{\varphi}}$ as a mapping from $L^{\infty}(\Omega,{\mathbb{R}}^{L})$ into $L^{\infty}(\Omega,{\mathbb{R}}^{L})\times\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ . Moreover, it holds that

\displaystyle D\mathcal{S}_{1}\in{\cal L}(L^{\infty}(\Omega,{\mathbb{R}}^{L}),L^{\infty}(\Omega,{\mathbb{R}}^{L})\times\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})),

and its directional derivative at ${\bm{\varphi}}\in L^{\infty}(\Omega,{\mathbb{R}}^{L})$ along a direction ${\bf h}\in L^{\infty}(\Omega,{\mathbb{R}}^{L})$ is given by

\displaystyle D\mathcal{S}_{1}({\bm{\varphi}})({\bf h})=({\bf h},{\overline{{\bf v}}}),

(4.1)

where ${\overline{{\bf v}}}\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ is the unique weak solution to the following system:


$\displaystyle-\mathop{\rm div}\nolimits\big{(}\overline{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}\mathcal{E}({\overline{{\bf u}}})+\overline{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}({\overline{{\bf v}}})\big{)}$	$\displaystyle=\bm{0}\qquad$	$\displaystyle\text{ in }\Omega,$	(4.2a)
$\displaystyle{\overline{{\bf v}}}$	$\displaystyle=\bm{0}\qquad$	$\displaystyle\text{ on }\overline{\Gamma}_{D},$	(4.2b)
$\displaystyle(\overline{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}\mathcal{E}({\overline{{\bf u}}})+\overline{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}({\overline{{\bf v}}})){\bf n}$	$\displaystyle=\bm{0}\qquad$	$\displaystyle\text{ on }\overline{\Gamma}_{N},$	(4.2c)

in the sense

\displaystyle\int_{\Omega}\overline{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}({\overline{{\bf v}}}):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}+\int_{\Omega}\overline{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}\mathcal{E}({\overline{{\bf u}}}):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}=0

(4.3)

for all ${\bm{\zeta}}\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ , where ${\overline{{\bf u}}}$ is the unique solution to (3.9) associated to ${\bm{\varphi}}$ obtained from Theorem 3.1.

Proof.

Firstly, the unique solvability of the linearised system (4.2) follows directly from the application of Proposition 3.1 upon choosing

\displaystyle{\bf u}={\overline{{\bf v}}},\quad{\mathbb{C}}=\overline{{\mathbb{C}}}({\bm{\varphi}}),\quad{\mathbb{F}}=\overline{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}\mathcal{E}({\overline{{\bf u}}}),\quad{\bf f}={\bf g}=\bm{0},\quad\Gamma_{D}=\overline{\Gamma}_{D},\quad\Gamma_{N}=\overline{\Gamma}_{N}.

Next, we take ${\bm{\varphi}}\in\mathcal{U}_{\rm ad}$ and ${\bf h}\in L^{\infty}(\Omega,{\mathbb{R}}^{L})$ such that ${\bm{\varphi}}^{{\bf h}}:={\bm{\varphi}}+{\bf h}\in\mathcal{U}_{\rm ad}$ , and set ${\overline{{\bf u}}}^{\bf h}=\mathcal{S}_{1}^{2}({\bm{\varphi}}^{\bf h})$ , i.e., $({\bm{\varphi}}^{\bf h},{\overline{{\bf u}}}^{\bf h})=\mathcal{S}_{1}({\bm{\varphi}}^{\bf h})$ . Since the first component of $\mathcal{S}_{1}$ is just the identity in $L^{\infty}(\Omega,{\mathbb{R}}^{L})$ , we only need to investigate the Fréchet differentiability of the second component $\mathcal{S}_{1}^{2}$ . In this direction, we denote

\displaystyle{\bf w}:={\overline{{\bf u}}}^{\bf h}-{\overline{{\bf u}}}-{\overline{{\bf v}}}\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})

with ${\overline{{\bf v}}}$ being the unique solution to the linearised system (4.2) associated to ${\bm{\varphi}}$ and ${\bf h}$ . Our aim is to show the existence of a positive constant $C$ such that

\displaystyle\|{\bf w}\|_{H^{1}(\Omega)}=\|\mathcal{S}_{1}^{2}({\bm{\varphi}}^{\bf h})-\mathcal{S}_{1}^{2}({\bm{\varphi}})-D\mathcal{S}_{1}^{2}({\bm{\varphi}}){\bf h}\|_{H^{1}(\Omega)}\leq C\|{\bf h}\|_{L^{\infty}(\Omega)}^{2},

(4.4)

which would then imply the Fréchet differentiability of the operator $\mathcal{S}_{1}^{2}$ . To this end, we subtract from (3.9) for ${\bm{\varphi}}^{\bf h}$ the sum of (3.9) for ${\bm{\varphi}}$ and (4.3) for ${\overline{{\bf v}}}$ to obtain

	$\displaystyle\int_{\Omega}\overline{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}({\bf w}):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}+\int_{\Omega}(\overline{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\overline{{\mathbb{C}}}({\bm{\varphi}}))(\mathcal{E}({\overline{{\bf u}}}^{\bf h})-\mathcal{E}({\overline{{\bf u}}})):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}$
	$\displaystyle\quad+\int_{\Omega}[\overline{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\overline{{\mathbb{C}}}({\bm{\varphi}})-\overline{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}]\mathcal{E}({\overline{{\bf u}}}):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}=0\quad\forall{\bm{\zeta}}\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}).$

Choosing ${\bm{\zeta}}={\bf w}$ and using (2.4) we infer that

	$\displaystyle C_{0}\mathopen{\\|}\mathcal{E}({\bf w})\mathclose{\\|}^{2}_{L^{2}(\Omega)}$	$\displaystyle\leq-\int_{\Omega}(\overline{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\overline{{\mathbb{C}}}({\bm{\varphi}}))(\mathcal{E}({\overline{{\bf u}}}^{\bf h})-\mathcal{E}({\overline{{\bf u}}})):\mathcal{E}({\bf w})\,\mathrm{dx}$
		$\displaystyle\quad-\int_{\Omega}[\overline{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\overline{{\mathbb{C}}}({\bm{\varphi}})-\overline{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}]\mathcal{E}({\overline{{\bf u}}}):\mathcal{E}({\bf w})\,\mathrm{dx}.$

Lipschitz continuity of $\overline{{\mathbb{C}}}$ yields a positive constant $C_{\overline{{\mathbb{C}}}^{\prime}}$ such that

\displaystyle|\overline{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\overline{{\mathbb{C}}}({\bm{\varphi}})-\overline{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}|\leq|{\bf h}|\int_{0}^{1}|\overline{{\mathbb{C}}}^{\prime}({\bm{\varphi}}+s{\bf h})-\overline{{\mathbb{C}}}^{\prime}({\bm{\varphi}})|\,\mathrm{d}s\leq C_{\overline{{\mathbb{C}}}^{\prime}}|{\bf h}|^{2},

(4.5)

and keeping in mind the estimate obtained from Theorem 3.2:

\displaystyle\mathopen{\|}{\overline{{\bf u}}}^{\bf h}-{\overline{{\bf u}}}\mathclose{\|}_{H^{1}(\Omega)}\leq C\mathopen{\|}{\bf h}\mathclose{\|}_{L^{\infty}(\Omega)},

(4.6)

we find that

	$\displaystyle\int_{\Omega}(\overline{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\overline{{\mathbb{C}}}({\bm{\varphi}}))(\mathcal{E}({\overline{{\bf u}}}^{\bf h})-\mathcal{E}({\overline{{\bf u}}})):\mathcal{E}({\bf w})\,\mathrm{dx}$
	$\displaystyle\quad\leq C\mathopen{\\|}{\bf h}\mathclose{\\|}_{L^{\infty}(\Omega)}\mathopen{\\|}\mathcal{E}({\overline{{\bf u}}}^{\bf h})-\mathcal{E}({\overline{{\bf u}}})\mathclose{\\|}_{L^{2}(\Omega)}\mathopen{\\|}\mathcal{E}({\bf w})\mathclose{\\|}_{L^{2}(\Omega)}\leq\frac{C_{0}}{4}\mathopen{\\|}\mathcal{E}({\bf w})\mathclose{\\|}_{L^{2}(\Omega)}^{2}+C\mathopen{\\|}{\bf h}\mathclose{\\|}_{L^{\infty}(\Omega)}^{4},$
	$\displaystyle\int_{\Omega}[\overline{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\overline{{\mathbb{C}}}({\bm{\varphi}})-\overline{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}]\mathcal{E}({\overline{{\bf u}}}):\mathcal{E}({\bf w})\,\mathrm{dx}$
	$\displaystyle\quad\leq C\mathopen{\\|}{\bf h}\mathclose{\\|}_{L^{\infty}(\Omega)}^{2}\mathopen{\\|}\mathcal{E}({\overline{{\bf u}}})\mathclose{\\|}_{L^{2}(\Omega)}\mathopen{\\|}\mathcal{E}({\bf w})\mathclose{\\|}_{L^{2}(\Omega)}\leq\frac{C_{0}}{4}\mathopen{\\|}\mathcal{E}({\bf w})\mathclose{\\|}_{L^{2}(\Omega)}^{2}+C\mathopen{\\|}{\bf h}\mathclose{\\|}_{L^{\infty}(\Omega)}^{4}.$

Then, by Korn’s inequality we infer (4.4), and whence the claimed Frechét differentiability of $\mathcal{S}_{1}$ . ∎

Before presenting the Fréchet differentiability of $\mathcal{S}_{2}$ , let us provide a formal discussion. Recall that ${\widehat{{\bf u}}}=\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}})$ and thus the directional derivative $D\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}})({\bf h},{\bf k})$ of $\mathcal{S}_{2}$ at $({\bm{\varphi}},{\overline{{\bf u}}})$ along a direction $({\bf h},{\bf k})$ will be given by

D\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}})({\bf h},{\bf k})=D_{{\bm{\varphi}}}\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}})({\bf h})+D_{{\overline{{\bf u}}}}\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}})({\bf k})

where $D_{{\bm{\varphi}}}$ and $D_{{\overline{{\bf u}}}}$ represent the partial derivatives with respect to ${\bm{\varphi}}$ and ${\overline{{\bf u}}}$ , respectively. Hence, we expect $\widehat{\bm{\vartheta}}:=D\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}})({\bf h},{\bf k})$ to be a sum of two functions ${\widehat{{\bf v}}}:=D_{{\bm{\varphi}}}\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}})({\bf h})$ and ${\widehat{{\bf w}}}:=D_{{\overline{{\bf u}}}}\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}})({\bf k})$ , and the Fréchet differentiability of $\mathcal{S}_{2}$ can be established by demonstrating

\displaystyle\|\mathcal{S}_{2}({\bm{\varphi}}+{\bf h},{\overline{{\bf u}}}+{\bf k})-\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}})-\widehat{\bm{\vartheta}}\|_{H^{1}(\Omega)}\leq C\|({\bf h},{\bf k})\|_{L^{\infty}(\Omega)\times H^{1}(\Omega)}^{2}.

(4.7)

The result is formulated as follows.

Theorem 4.2.

The solution operator $\mathcal{S}_{2}$ is Fréchet differentiable at $({\bm{\varphi}},{\overline{{\bf u}}})$ as a mapping from $L^{\infty}(\Omega,{\mathbb{R}}^{L})\times\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ into $\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ . Furthermore,

\displaystyle D\mathcal{S}_{2}\in{\cal L}(L^{\infty}(\Omega,{\mathbb{R}}^{L})\times\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}),\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})),

and its directional derivative at $({\bm{\varphi}},{\overline{{\bf u}}})\in L^{\infty}(\Omega,{\mathbb{R}}^{L})\times\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ along a direction $({\bf h},{\bf k})\in L^{\infty}(\Omega,{\mathbb{R}}^{L})\times\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ is given by

\displaystyle\widehat{\bm{\vartheta}}:=D\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}})({\bf h},{\bf k})=D_{\bm{\varphi}}\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}})({\bf h})+D_{\overline{{\bf u}}}\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}})({\bf k})=:{\widehat{{\bf v}}}+{\widehat{{\bf w}}},

where $\widehat{\bm{\vartheta}}\in\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ is the unique weak solution to the following system:


$\displaystyle-\mathop{\rm div}\nolimits\big{(}\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}(\mathcal{E}({\widehat{{\bf u}}})-\chi({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}}))\big{)}$
$\displaystyle\quad-\mathop{\rm div}\nolimits\big{(}\widehat{{\mathbb{C}}}({\bm{\varphi}})(\mathcal{E}(\widehat{\bm{\vartheta}})-\chi^{\prime}({\bm{\varphi}}){\bf h}\mathcal{E}({\overline{{\bf u}}})-\chi({\bm{\varphi}})\mathcal{E}({\bf k}))\big{)}$	$\displaystyle=\bm{0}\qquad$	$\displaystyle\text{ in }\Omega,$	(4.8a)
$\displaystyle\widehat{\bm{\vartheta}}$	$\displaystyle=\bm{0}\qquad$	$\displaystyle\text{ on }\widehat{\Gamma}_{D},$	(4.8b)
$\displaystyle\big{(}\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}(\mathcal{E}({\widehat{{\bf u}}})-\chi({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}}))\big{)}{\bf n}$
$\displaystyle\quad+\big{(}\widehat{{\mathbb{C}}}({\bm{\varphi}})(\mathcal{E}(\widehat{\bm{\vartheta}})-\chi^{\prime}({\bm{\varphi}}){\bf h}\mathcal{E}({\overline{{\bf u}}})-\chi({\bm{\varphi}})\mathcal{E}({\bf k}))\big{)}{\bf n}$	$\displaystyle=\bm{0}\qquad$	$\displaystyle\text{ on }\widehat{\Gamma}_{N},$	(4.8c)

in the sense that

		$\displaystyle\int_{\Omega}\widehat{{\mathbb{C}}}({\bm{\varphi}})(\mathcal{E}(\widehat{\bm{\vartheta}})-\chi^{\prime}({\bm{\varphi}}){\bf h}\mathcal{E}({\overline{{\bf u}}})-\chi({\bm{\varphi}})\mathcal{E}({\bf k})):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}$		(4.9)
		$\displaystyle\quad+\int_{\Omega}\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}(\mathcal{E}({\widehat{{\bf u}}})-\chi({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}})):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}=0$		(4.9)

for all ${\bm{\zeta}}\in\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ , where ${\overline{{\bf u}}}$ is the unique solution to (3.9) associated to ${\bm{\varphi}}$ and ${\widehat{{\bf u}}}$ is the unique solution to (3.10) associated to $({\bm{\varphi}},{\overline{{\bf u}}})$ . Moreover, ${\widehat{{\bf v}}}\in\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ and ${\widehat{{\bf w}}}\in\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ are the unique solutions to the following equations

	$\displaystyle\int_{\Omega}\widehat{{\mathbb{C}}}({\bm{\varphi}})(\mathcal{E}({\widehat{{\bf v}}})-\chi^{\prime}({\bm{\varphi}}){\bf h}\mathcal{E}({\overline{{\bf u}}})):\mathcal{E}({\bm{\zeta}})+\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}(\mathcal{E}({\widehat{{\bf u}}})-\chi({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}})):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}=0,$		(4.10)
	$\displaystyle\int_{\Omega}\widehat{{\mathbb{C}}}({\bm{\varphi}})(\mathcal{E}({\widehat{{\bf w}}})-\chi({\bm{\varphi}})\mathcal{E}({\bf k})):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}=0,$		(4.11)

for all ${\bm{\zeta}}\in\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ .

Proof.

Unique solvability of (4.9), (4.10) and (4.11) are obtained by application of Theorem 3.1, and thus we focus on demonstrating the crucial estimate (4.7). Let ${\bm{\varphi}}\in\mathcal{U}_{\rm ad}$ , ${\bf h}\in L^{\infty}(\Omega,{\mathbb{R}}^{L})$ such that ${\bm{\varphi}}^{\bf h}={\bm{\varphi}}+{\bf h}\in\mathcal{U}_{\rm ad}$ and ${\bf k}\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ . Setting

{\widehat{{\bf u}}}=\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}}),\quad{\widehat{{\bf u}}}^{{\bf k}}:=\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}}+{\bf k}),\quad{\widehat{{\bf u}}}^{{\bf h},{\bf k}}:=\mathcal{S}_{2}({\bm{\varphi}}+{\bf h},{\overline{{\bf u}}}+{\bf k}),

then, setting $\bm{\xi}:={\widehat{{\bf u}}}^{{\bf h},{\bf k}}-{\widehat{{\bf u}}}-\widehat{\bm{\vartheta}}$ , (4.7) is equivalent to

\|\bm{\xi}\|_{H^{1}(\Omega)}\leq C\|({\bf h},{\bf k})\|_{L^{\infty}(\Omega)\times H^{1}(\Omega)}^{2}.

We observe that by subtracting from (3.10) for $({\bm{\varphi}}+{\bf h},{\overline{{\bf u}}}+{\bf k},{\widehat{{\bf u}}}^{{\bf h},{\bf k}})$ the sum of (3.10) for $({\bm{\varphi}},{\overline{{\bf u}}},{\widehat{{\bf u}}})$ and (4.9) for $\widehat{\bm{\vartheta}}$ , we obtain

		$\displaystyle\int_{\Omega}\big{(}\widehat{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})(\mathcal{E}({\widehat{{\bf u}}}^{{\bf h},{\bf k}})-\chi({\bm{\varphi}}+{\bf h})\mathcal{E}({\overline{{\bf u}}}+{\bf k})\big{)}:\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}$		(4.12)
		$\displaystyle\quad-\int_{\Omega}\widehat{{\mathbb{C}}}({\bm{\varphi}})(\mathcal{E}({\widehat{{\bf u}}})-\chi({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}})):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}$
		$\displaystyle\quad-\int_{\Omega}\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}(\mathcal{E}({\widehat{{\bf u}}})-\chi({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}})):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}$
		$\displaystyle\quad-\int_{\Omega}\widehat{{\mathbb{C}}}({\bm{\varphi}})(\mathcal{E}(\widehat{\bm{\vartheta}})-\chi^{\prime}({\bm{\varphi}}){\bf h}\mathcal{E}({\overline{{\bf u}}})-\chi({\bm{\varphi}})\mathcal{E}({\bf k})):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}=0\quad\forall{\bm{\zeta}}\in\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}).$

Before proceeding with some computations, let us point out the following identities:

		$\displaystyle\widehat{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})\mathcal{E}({\widehat{{\bf u}}}^{{\bf h},{\bf k}})-\widehat{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}({\widehat{{\bf u}}})-\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}\mathcal{E}({\widehat{{\bf u}}})-\widehat{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}(\widehat{\bm{\vartheta}})$		(4.13)
		$\displaystyle\quad=\widehat{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}(\bm{\xi})+(\widehat{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\widehat{{\mathbb{C}}}({\bm{\varphi}}))(\mathcal{E}({\widehat{{\bf u}}}^{{\bf h},{\bf k}})-\mathcal{E}({\widehat{{\bf u}}}^{\bf h}))$
		$\displaystyle\qquad+(\widehat{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\widehat{{\mathbb{C}}}({\bm{\varphi}}))(\mathcal{E}({\widehat{{\bf u}}}^{{\bf h}})-\mathcal{E}({\widehat{{\bf u}}}))$
		$\displaystyle\qquad+[\widehat{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\widehat{{\mathbb{C}}}({\bm{\varphi}})-\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}]\mathcal{E}({\widehat{{\bf u}}})=:\widehat{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}(\bm{\xi})+{\cal Y}_{1},$

and

		$\displaystyle-\widehat{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})\chi({\bm{\varphi}}+{\bf h})\mathcal{E}({\overline{{\bf u}}}+{\bf k})+\widehat{{\mathbb{C}}}({\bm{\varphi}})\chi({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}})+\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}\chi({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}})$		(4.14)
		$\displaystyle\qquad+\widehat{{\mathbb{C}}}({\bm{\varphi}})\chi^{\prime}({\bm{\varphi}}){\bf h}\mathcal{E}({\overline{{\bf u}}})+\widehat{{\mathbb{C}}}({\bm{\varphi}})\chi({\bm{\varphi}})\mathcal{E}({\bf k})$
		$\displaystyle\quad=-[\widehat{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\widehat{{\mathbb{C}}}({\bm{\varphi}})-\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}]\chi({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}})$
		$\displaystyle\qquad-\widehat{{\mathbb{C}}}({\bm{\varphi}})[\chi({\bm{\varphi}}+{\bf h})-\chi({\bm{\varphi}})-\chi^{\prime}({\bm{\varphi}}){\bf h}]\mathcal{E}({\overline{{\bf u}}})$
		$\displaystyle\qquad-(\widehat{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\widehat{{\mathbb{C}}}({\bm{\varphi}}))(\chi({\bm{\varphi}}+{\bf h})-\chi({\bm{\varphi}}))\mathcal{E}({\overline{{\bf u}}})$
		$\displaystyle\qquad-(\widehat{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\widehat{{\mathbb{C}}}({\bm{\varphi}}))(\chi({\bm{\varphi}}+{\bf h})-\chi({\bm{\varphi}}))\mathcal{E}({\bf k})$
		$\displaystyle\qquad-\widehat{{\mathbb{C}}}({\bm{\varphi}})(\chi({\bm{\varphi}}+{\bf h})-\chi({\bm{\varphi}}))\mathcal{E}({\bf k})$
		$\displaystyle\qquad-(\widehat{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\widehat{{\mathbb{C}}}({\bm{\varphi}}))\chi({\bm{\varphi}})\mathcal{E}({\bf k})=:{\cal Y}_{2}.$

Similar to (4.5), we have, for positive constants $C_{\chi^{\prime}}$ and $C_{\widehat{{\mathbb{C}}}^{\prime}}$ , that

	$\displaystyle\|\chi({\bm{\varphi}}+{\bf h})-\chi({\bm{\varphi}})-\chi^{\prime}({\bm{\varphi}}){\bf h}\|$	$\displaystyle\leq\|{\bf h}\|\int_{0}^{1}\|\chi^{\prime}({\bm{\varphi}}+s{\bf h})-\chi^{\prime}({\bm{\varphi}})\|\,\mathrm{d}s\leq C_{\chi^{\prime}}\|{\bf h}\|^{2},$		(4.15)
	$\displaystyle\|\widehat{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\widehat{{\mathbb{C}}}({\bm{\varphi}})-\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}\|$	$\displaystyle\leq\|{\bf h}\|\int_{0}^{1}\|\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}+s{\bf h})-\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}})\|\,\mathrm{d}s\leq C_{\widehat{{\mathbb{C}}}^{\prime}}\|{\bf h}\|^{2},$		(4.15)

and upon choosing ${\bm{\zeta}}=\bm{\xi}$ in (4.12), we infer that

\displaystyle\mathopen{\|}\mathcal{E}(\bm{\xi})\mathclose{\|}_{L^{2}(\Omega)}\leq C\|{\cal Y}_{1}\|_{L^{2}(\Omega)}+C\|{\cal Y}_{2}\|_{L^{2}(\Omega)}.

(4.16)

Then, employing the Young and Hölder inequalities, (4.6), as well as the stability estimates

\mathopen{\|}{\widehat{{\bf u}}}^{{\bf h},{\bf k}}-{\widehat{{\bf u}}}^{{\bf h}}\mathclose{\|}_{H^{1}(\Omega)}\leq C\mathopen{\|}{\bf k}\mathclose{\|}_{H^{1}(\Omega)},\quad\mathopen{\|}{\widehat{{\bf u}}}^{{\bf h}}-{\widehat{{\bf u}}}\mathclose{\|}_{H^{1}(\Omega)}\leq C\mathopen{\|}{\bf h}\mathclose{\|}_{L^{\infty}(\Omega)}

deduced from (3.12), we find that

	$\displaystyle\\|{\cal Y}_{1}\\|_{L^{2}(\Omega)}$	$\displaystyle\leq C\\|{\bf h}\\|_{L^{\infty}(\Omega)}\Big{(}\mathopen{\\|}{\widehat{{\bf u}}}^{{\bf h},{\bf k}}-{\widehat{{\bf u}}}^{\bf h}\mathclose{\\|}_{H^{1}(\Omega)}+\mathopen{\\|}{\widehat{{\bf u}}}^{{\bf h}}-{\widehat{{\bf u}}}\mathclose{\\|}_{H^{1}(\Omega)}+\mathopen{\\|}{\bf h}\mathclose{\\|}_{L^{\infty}(\Omega)}\mathopen{\\|}\mathcal{E}({\widehat{{\bf u}}})\mathclose{\\|}_{L^{2}(\Omega)}\Big{)}$
		$\displaystyle\leq C\mathopen{\\|}({\bf h},{\bf k})\mathclose{\\|}_{L^{\infty}(\Omega)\times H^{1}(\Omega)}^{2},$
	$\displaystyle\\|{\cal Y}_{2}\\|_{L^{2}(\Omega)}$	$\displaystyle\leq C\\|{\bf h}\\|_{L^{\infty}(\Omega)}\Big{(}\\|{\bf h}\\|_{L^{\infty}(\Omega)}\\|\mathcal{E}({\overline{{\bf u}}})\\|_{L^{2}(\Omega)}+\mathopen{\\|}{\bf h}\mathclose{\\|}_{L^{\infty}(\Omega)}\\|{\bf k}\\|_{H^{1}(\Omega)}+\\|{\bf k}\\|_{H^{1}(\Omega)}\Big{)}$
		$\displaystyle\leq C\mathopen{\\|}({\bf h},{\bf k})\mathclose{\\|}_{L^{\infty}(\Omega)\times H^{1}(\Omega)}^{2},$

and thus we obtain by Korn’s inequality

\mathopen{\|}\bm{\xi}\mathclose{\|}_{H^{1}(\Omega)}\leq C\mathopen{\|}({\bf h},{\bf k})\mathclose{\|}_{L^{\infty}(\Omega)\times H^{1}(\Omega)}^{2},

which verifies the Fréchet differentiability of $\mathcal{S}_{2}$ . Furthermore, it is clear that by the uniqueness of the solutions, the sum ${\widehat{{\bf v}}}+{\widehat{{\bf w}}}$ is equal to $\widehat{\bm{\vartheta}}$ . To establish the identification of the partial derivatives $D_{{\bm{\varphi}}}\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}})({\bf h})={\widehat{{\bf v}}}$ and $D_{{\overline{{\bf u}}}}\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}})({\bf k})={\widehat{{\bf w}}}$ , we argue as follows: Consider ${\bf k}=\bm{0}$ , so that from (4.11) we obtain that ${\widehat{{\bf w}}}=\bm{0}$ and $\widehat{\bm{\vartheta}}={\widehat{{\bf v}}}$ . Then, in (4.12) with ${\bf k}=\bm{0}$ , we now have for $\bm{\xi}={\widehat{{\bf u}}}^{{\bf h}}-{\widehat{{\bf u}}}-{\widehat{{\bf v}}}$ the estimate (4.16), where

	$\displaystyle{\cal Y}_{1}$	$\displaystyle=(\widehat{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\widehat{{\mathbb{C}}}({\bm{\varphi}}))(\mathcal{E}({\widehat{{\bf u}}}^{{\bf h}})-\mathcal{E}({\widehat{{\bf u}}}))+[\widehat{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\widehat{{\mathbb{C}}}({\bm{\varphi}})-\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}]\mathcal{E}({\widehat{{\bf u}}}),$
	$\displaystyle{\cal Y}_{2}$	$\displaystyle=-[\widehat{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\widehat{{\mathbb{C}}}({\bm{\varphi}})-\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}]\chi({\bm{\varphi}})\mathcal{E}({\overline{{\bf u}}})$
		$\displaystyle\quad-\widehat{{\mathbb{C}}}({\bm{\varphi}})[\chi({\bm{\varphi}}+{\bf h})-\chi({\bm{\varphi}})-\chi^{\prime}({\bm{\varphi}}){\bf h}]\mathcal{E}({\overline{{\bf u}}})$
		$\displaystyle\quad-(\widehat{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\widehat{{\mathbb{C}}}({\bm{\varphi}}))(\chi({\bm{\varphi}}+{\bf h})-\chi({\bm{\varphi}}))\mathcal{E}({\overline{{\bf u}}}),$

where we made use of ${\widehat{{\bf u}}}^{{\bf h},\bm{0}}={\widehat{{\bf u}}}^{{\bf h}}$ . A short calculation shows that

\|{\widehat{{\bf u}}}^{{\bf h}}-{\widehat{{\bf u}}}-{\widehat{{\bf v}}}\|_{H^{1}(\Omega)}\leq C\|{\bf h}\|_{L^{\infty}(\Omega)}^{2},

which entails that $D_{{\bm{\varphi}}}\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}})({\bf h})={\widehat{{\bf v}}}$ . The other identification $D_{{\overline{{\bf u}}}}\mathcal{S}_{2}({\bm{\varphi}},{\overline{{\bf u}}})({\bf k})={\widehat{{\bf w}}}$ is in fact more straightforward as $\mathcal{S}_{2}({\bm{\varphi}},\cdot)$ is a linear operator. This completes the proof. ∎

4.2 Adjoint systems

We now move to the investigation of some adjoint systems which, as typically happens in constrained minimisation problems, will allow us to simplify the formulation of the optimality conditions of $(\bf P)$ .

Theorem 4.3.

Under (A1)–(A6), for every ${\bm{\varphi}}\in L^{\infty}(\Omega,{\mathbb{R}}^{L})$ , there exists a unique solution $\widehat{{\bf q}}\in\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ to


$\displaystyle-\mathop{\rm div}\nolimits\big{(}\widehat{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}(\widehat{{\bf q}})\big{)}$	$\displaystyle=\bm{0}\qquad$	$\displaystyle\text{ in }\Omega,$	(4.17a)
$\displaystyle\widehat{{\bf q}}$	$\displaystyle=\bm{0}\qquad$	$\displaystyle\text{ on }\widehat{\Gamma}_{D},$	(4.17b)
$\displaystyle(\widehat{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}(\widehat{{\bf q}})){\bf n}$	$\displaystyle=W({\widehat{{\bf u}}}-{\bf u}^{\rm tar})\mathrm{\mathcal{X}}_{\Gamma^{\rm tar}}\qquad$	$\displaystyle\text{ on }\widehat{\Gamma}_{N},$	(4.17c)

in the sense

\displaystyle\int_{\Omega}\widehat{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}(\widehat{\bm{q}}):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}=\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}-{\bf u}^{\rm tar})\cdot{\bm{\zeta}}\,\mathrm{d}\mathcal{H}^{d-1}\quad\forall{\bm{\zeta}}\in\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}).

(4.18)

Moreover, there exists a unique solution $\overline{{\bf p}}\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ to


$\displaystyle-\mathop{\rm div}\nolimits\big{(}\overline{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}(\overline{{\bf p}})-\widehat{{\mathbb{C}}}({\bm{\varphi}})\chi({\bm{\varphi}})\mathcal{E}(\widehat{{\bf q}}))\big{)}$	$\displaystyle=\bm{0}\qquad$	$\displaystyle\text{ in }\Omega,$	(4.19a)
$\displaystyle\overline{{\bf p}}$	$\displaystyle=\bm{0}\qquad$	$\displaystyle\text{ on }\overline{\Gamma}_{D},$	(4.19b)
$\displaystyle(\overline{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}(\overline{{\bf p}})-\widehat{{\mathbb{C}}}({\bm{\varphi}})\chi({\bm{\varphi}})\mathcal{E}(\widehat{{\bf q}})){\bf n}$	$\displaystyle=\bm{0}\qquad$	$\displaystyle\text{ on }\overline{\Gamma}_{N},$	(4.19c)

in the sense

\displaystyle\int_{\Omega}\overline{{\mathbb{C}}}({\bm{\varphi}})\mathcal{E}(\overline{\bm{p}}):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}=\int_{\Omega}\widehat{{\mathbb{C}}}({\bm{\varphi}})\chi({\bm{\varphi}})\mathcal{E}(\widehat{\bm{q}}):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}\quad\forall{\bm{\zeta}}\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}),

(4.20)

where $\widehat{{\bf q}}\in\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ is the unique solution to (4.18).

As the proof of existence and uniqueness is simply an application of Theorem 3.1, we omit the details.

Remark 4.1.

Notice that the adjoint variable $\overline{\bm{p}}$ to the Stage 1 deformation ${\overline{{\bf u}}}$ is dependent on the adjoint variable $\widehat{\bm{q}}$ to the Stage 2 deformation ${\widehat{{\bf u}}}$ . This backwards dependence has some parallels with the adjoint systems associated to time-dependent state equations, which have to be solved backwards in time.

4.3 First order optimality conditions

Theorem 4.4.

Under (A1)–(A6), let ${\bm{\varphi}}^{*}\in\mathcal{U}_{\rm ad}$ be a minimiser to $J_{\rm red}$ with corresponding states ${\overline{{\bf u}}}^{*}=\mathcal{S}_{1}^{2}({\bm{\varphi}}^{*})$ , ${\widehat{{\bf u}}}^{*}=\mathcal{S}({\bm{\varphi}}^{*})$ , and adjoint variables $(\overline{\bm{p}},\widehat{\bm{q}})$ as unique solutions to (4.18) and (4.20) corresponding to $({\bm{\varphi}}^{*},{\overline{{\bf u}}}^{*},{\widehat{{\bf u}}}^{*})$ . Then, it necessarily holds that

	$\displaystyle-\int_{\Omega}\overline{{\mathbb{C}}}^{\prime}({\bm{\varphi}}^{})(\bm{\phi}-{\bm{\varphi}}^{})\mathcal{E}({\overline{{\bf u}}}^{*}):\mathcal{E}(\overline{{\bf p}})\,\mathrm{dx}$
	$\displaystyle\quad-\int_{\Omega}\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}^{})(\bm{\phi}-{\bm{\varphi}}^{})(\mathcal{E}({\widehat{{\bf u}}}^{})-\chi({\bm{\varphi}}^{})\mathcal{E}({\overline{{\bf u}}}^{*})):\mathcal{E}(\widehat{{\bf q}})\,\mathrm{dx}$
	$\displaystyle\quad+\int_{\Omega}\widehat{{\mathbb{C}}}({\bm{\varphi}}^{})\chi^{\prime}({\bm{\varphi}}^{})(\bm{\phi}-{\bm{\varphi}}^{})\mathcal{E}({\overline{{\bf u}}}^{}):\mathcal{E}(\widehat{{\bf q}})\,\mathrm{dx}$
	$\displaystyle\quad+2\gamma\varepsilon\int_{\Omega}\nabla{\bm{\varphi}}^{}\cdot\nabla(\bm{\phi}-{\bm{\varphi}}^{})\,\mathrm{dx}+\frac{\gamma}{\varepsilon}\int_{\Omega}\widetilde{\Psi}_{,{\bm{\varphi}}}({\bm{\varphi}}^{})\cdot(\bm{\phi}-{\bm{\varphi}}^{})\,\mathrm{dx}\geq 0\quad\forall\bm{\phi}\in\mathcal{U}_{\rm ad},$		(4.21)

where we set $\widetilde{\Psi}_{,{\bm{\varphi}}}$ as the vector of partial derivatives of $\widetilde{\Psi}$ .

Proof.

As $\mathcal{U}_{\rm ad}$ is a non-empty, closed and convex set, standard results in optimal control and convex analysis yield that the first order necessary optimality condition for ${\bm{\varphi}}^{*}$ is

\displaystyle\langle DJ_{\rm red}({\bm{\varphi}}^{*}),\bm{\phi}-{\bm{\varphi}}^{*}\rangle\geq 0\quad\forall\bm{\phi}\in\mathcal{U}_{\rm ad},

which reads, in view of Theorem 4.1 and Theorem 4.2, as

		$\displaystyle\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}^{*}-{\bf u}^{\rm tar})\cdot\widehat{\bm{\vartheta}}\,\mathrm{d}\mathcal{H}^{d-1}$		(4.22)
		$\displaystyle\quad+2\gamma\varepsilon\int_{\Omega}\nabla{\bm{\varphi}}^{}\cdot\nabla(\bm{\phi}-{\bm{\varphi}}^{})\,\mathrm{dx}+\frac{\gamma}{\varepsilon}\int_{\Omega}\widetilde{\Psi}_{,{\bm{\varphi}}}({\bm{\varphi}}^{})\cdot(\bm{\phi}-{\bm{\varphi}}^{})\,\mathrm{dx}\geq 0\quad\forall{\bm{\varphi}}\in\mathcal{U}_{\rm ad},$		(4.22)

where $\widehat{\bm{\vartheta}}$ is the unique solution to (4.9) with ${\bm{\varphi}}={\bm{\varphi}}^{*}$ , ${\overline{{\bf u}}}={\overline{{\bf u}}}^{*}$ , ${\widehat{{\bf u}}}={\widehat{{\bf u}}}^{*}$ , ${\bf h}=\bm{\phi}-{\bm{\varphi}}^{*}$ and ${\bf k}={\overline{{\bf v}}}$ from (4.3) (also with ${\bf h}=\bm{\phi}-{\bm{\varphi}}^{*}$ ). To simplify the above relation, the procedure is to compare the equalities obtained from (4.3) with ${\bm{\zeta}}=\overline{\bm{p}}$ , (4.9) with ${\bf k}={\overline{{\bf v}}}$ and ${\bm{\zeta}}=\widehat{\bm{q}}$ , (4.18) with ${\bm{\zeta}}=\widehat{\bm{\vartheta}}$ and (4.20) with ${\bm{\zeta}}={\overline{{\bf v}}}$ . A short calculation then shows

	$\displaystyle\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}^{*}-{\bf u}^{\rm tar})\cdot\widehat{\bm{\vartheta}}\,\mathrm{d}\mathcal{H}^{d-1}$	$\displaystyle=-\int_{\Omega}\overline{{\mathbb{C}}}^{\prime}({\bm{\varphi}}^{})({\bm{\varphi}}-{\bm{\varphi}}^{})\mathcal{E}({\overline{{\bf u}}}^{*}):\mathcal{E}(\overline{{\bf p}})\,\mathrm{dx}$
		$\displaystyle\quad-\int_{\Omega}\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}^{})({\bm{\varphi}}-{\bm{\varphi}}^{})(\mathcal{E}({\widehat{{\bf u}}}^{})-\chi({\bm{\varphi}}^{})\mathcal{E}({\overline{{\bf u}}}^{*})):\mathcal{E}(\widehat{{\bf q}})\,\mathrm{dx}$
		$\displaystyle\quad+\int_{\Omega}\widehat{{\mathbb{C}}}({\bm{\varphi}}^{})\chi^{\prime}({\bm{\varphi}}^{})({\bm{\varphi}}-{\bm{\varphi}}^{})\mathcal{E}({\overline{{\bf u}}}^{}):\mathcal{E}(\widehat{{\bf q}})\,\mathrm{dx},$

which allows us to remove the dependence of $\widehat{\bm{\vartheta}}$ from (4.22) and leads to (4.21). ∎

5 Sharp interface asymptotics

In this section we study the behaviour of solutions in the sharp interface limit $\varepsilon\to 0$ . For $\varepsilon>0$ , we denote

	$\displaystyle E_{\varepsilon}({\bm{\varphi}})$	$\displaystyle=\varepsilon\int_{\Omega}\|\nabla{\bm{\varphi}}\|^{2}\,\mathrm{dx}+\frac{1}{\varepsilon}\int_{\Omega}\Psi({\bm{\varphi}})\,\mathrm{dx},$
	$\displaystyle G({\bm{\varphi}})$	$\displaystyle=\frac{1}{2}\int_{\Gamma^{\rm tar}}W(\mathcal{S}({\bm{\varphi}})-{\bf u}^{\rm tar})\cdot(\mathcal{S}({\bm{\varphi}})-{\bf u}^{\rm tar})\,\mathrm{d}\mathcal{H}^{d-1},$

where $\mathcal{S}$ is the solution operator defined in (3.18), so that the corresponding reduced functional can be expressed as the sum $J_{\rm red}^{\varepsilon}({\bm{\varphi}})=G({\bm{\varphi}})+\gamma E_{\varepsilon}({\bm{\varphi}})$ . The asymptotic behaviour of solutions can be studied under the framework of $\Gamma$ -convergence. In order to state the result some preparation is needed. A function ${\bm{\varphi}}\in L^{1}(\Omega,{\mathbb{R}}^{L})$ is termed a function of bounded variation in $\Omega$ , written as ${\bm{\varphi}}\in\mathrm{BV}(\Omega,{\mathbb{R}}^{L})$ if there exists a matrix-valued measure $D{\bm{\varphi}}$ of dimension $L\times d$ on $\Omega$ such that

\sum_{j=1}^{L}\int_{\Omega}\varphi_{j}(\mathop{\rm div}\nolimits{\bm{\psi}})_{j}\,\mathrm{dx}=-\sum_{j=1}^{L}\sum_{i=1}^{d}\int_{\Omega}\psi_{i}^{j}dD_{i}\varphi_{j},

for all $\bm{\psi}=(\psi_{i}^{j})_{1\leq i\leq d,1\leq j\leq L}$ where $\psi_{i}^{j}\in C^{1}_{c}(\Omega)$ . Let $T=\{\bm{e}_{1},\dots,\bm{e}_{L}\}$ where $\Psi^{-1}(0)=T$ . For ${\bm{\varphi}}\in\mathrm{BV}(\Omega,T)$ we set, for $i\in\{1,...,L\}$ ,

E_{\bm{\varphi}}^{i}=\{\bm{x}\in\Omega\,:\,{\bm{\varphi}}(\bm{x})=\bm{e}_{i}\}

and define the essential boundary $\partial^{*}E_{{\bm{\varphi}}}^{i}$ as

\partial^{*}E_{\bm{\varphi}}^{i}=\Big{\{}\bm{x}\in{\mathbb{R}}^{d}\,:\,\lim_{\rho\to 0^{+}}\frac{|E_{\bm{\varphi}}^{i}\cap B_{\rho}(\bm{x})|}{|B_{\rho}(\bm{x})|}\notin\{0,1\}\Big{\}}

where, for any $\rho>0$ , $B_{\rho}(\bm{x})$ is the $\rho$ -ball in ${\mathbb{R}}^{d}$ centered in $\bm{x}$ , i.e., $B_{\rho}(\bm{x})=\{\bm{y}\in{\mathbb{R}}^{d}\,:\,|\bm{y}-\bm{x}|<\rho\}$ . Consider the extended functionals

	$\displaystyle{\mathbb{E}}_{\varepsilon}({\bm{\varphi}})$	$\displaystyle:=\begin{cases}E_{\varepsilon}({\bm{\varphi}})&\text{ if }{\bm{\varphi}}\in H^{1}(\Omega,{\mathbb{R}}^{L}),\\ +\infty&\text{ elsewhere in }L^{1}(\Omega,{\mathbb{R}}^{L}),\end{cases}$
	$\displaystyle{\mathbb{E}}_{0}({\bm{\varphi}})$	$\displaystyle:=\begin{cases}\sum_{i,j=1,\,i<j}^{L}b_{ij}\mathcal{H}^{d-1}(\Omega\cap\partial^{}E_{{\bm{\varphi}}}^{i}\cap\partial^{}E_{{\bm{\varphi}}}^{j})&\text{ if }{\bm{\varphi}}\in\mathrm{BV}(\Omega,T),\\ +\infty&\text{ elsewhere in }L^{1}(\Omega,{\mathbb{R}}^{L}),\end{cases}$

with constants $b_{ij}$ defined as

b_{ij}=\inf\Big{\{}\int_{0}^{1}\Psi^{1/2}(\bm{\gamma}(t))|\bm{\gamma}^{\prime}(t)|\,\mathrm{d}t\,:\,\bm{\gamma}\in C^{1}([0,1];\Delta^{L}),\,\bm{\gamma}(0)=\bm{e}_{i},\,\bm{\gamma}(1)=\bm{e}_{j}\Big{\}}.

Then, the $\Gamma$ -convergence of ${\mathbb{E}}_{\varepsilon}$ to ${\mathbb{E}}_{0}$ as $\varepsilon\to 0$ is expressed as the following assertions:

•

Liminf property. If $\{{\bm{\varphi}}^{\varepsilon}\}_{\varepsilon>0}$ is a sequence such that $\liminf_{\varepsilon\to 0}\mathbb{E}_{\varepsilon}({\bm{\varphi}}^{\varepsilon})<\infty$ and ${\bm{\varphi}}^{\varepsilon}\to{\bm{\varphi}}^{0}$ in $L^{1}(\Omega,{\mathbb{R}}^{L})$ , then ${\bm{\varphi}}^{0}\in\mathrm{BV}(\Omega,T)$ with $\mathbb{E}_{0}({\bm{\varphi}}^{0})\leq\liminf_{\varepsilon\to 0}\mathbb{E}_{\varepsilon}({\bm{\varphi}}^{\varepsilon})$ .
•

Limsup property. For any ${\bm{\varphi}}^{0}\in L^{1}(\Omega,T)$ , there exists a sequence $\{{\bm{\varphi}}^{\varepsilon}\}_{\varepsilon>0}\subset H^{1}(\Omega,{\mathbb{R}}^{L})$ , such that ${\bm{\varphi}}^{\varepsilon}\to{\bm{\varphi}}^{0}$ in $L^{1}(\Omega,{\mathbb{R}}^{L})$ and $\limsup_{\varepsilon\to 0}{\mathbb{E}}_{\varepsilon}({\bm{\varphi}}^{\varepsilon})\leq\mathbb{E}_{0}({\bm{\varphi}}^{0})$ .
•

Compactness property. Let $\{{\bm{\varphi}}^{\varepsilon}\}_{\varepsilon>0}$ be a sequence such that $\sup_{\varepsilon}\mathbb{E}_{\varepsilon}({\bm{\varphi}}^{\varepsilon})<\infty$ . Then, there exists a non-relabelled subsequence and a function ${\bm{\varphi}}^{0}\in\mathrm{BV}(\Omega,T)$ such that ${\bm{\varphi}}^{\varepsilon}\to{\bm{\varphi}}^{0}$ in $L^{1}(\Omega,{\mathbb{R}}^{L})$ .

For a proof we refer to [8, Thm. 2.5 and Prop. 4.1], see also [11, Thm. 3.1 and Rmk. 3.3] with the choice $f(z,X)=|X|^{2}$ .

5.1 Convergences of minimisers

Lemma 5.1.

For each $\varepsilon\in(0,1]$ , let ${\bm{\varphi}}^{\varepsilon}\in\mathcal{U}_{\rm ad}$ denote a minimiser to the extended reduced cost functional $J_{\rm red}^{\varepsilon}({\bm{\varphi}})=G({\bm{\varphi}})+\gamma\mathbb{E}_{\varepsilon}({\bm{\varphi}})$ . Then, there exists a non-relabelled subsequence $\varepsilon\to 0$ and a limit function ${\bm{\varphi}}_{0}\in\mathrm{BV}(\Omega,T)$ such that ${\bm{\varphi}}^{\varepsilon}\to{\bm{\varphi}}_{0}$ strongly in $L^{1}(\Omega,{\mathbb{R}}^{L})$ , $\lim_{\varepsilon\to 0}J_{\rm red}^{\varepsilon}({\bm{\varphi}}^{\varepsilon})=J_{\rm red}^{0}({\bm{\varphi}}_{0})$ , where

\displaystyle J_{\rm red}^{0}({\bm{\varphi}})=G({\bm{\varphi}})+\gamma\mathbb{E}_{0}({\bm{\varphi}})\text{ for }{\bm{\varphi}}\in\mathrm{BV}(\Omega,T),

and ${\bm{\varphi}}_{0}$ is a minimiser to $J_{\rm red}^{0}$ .

Proof.

The proof relies on the $\Gamma$ -convergence of the Ginzburg–Landau functional and the stability of $\Gamma$ -convergence under continuous perturbations. By Corollary 3.1, and the continuity of the trace operator, we see that $G$ is continuous. For arbitrary $\bm{\psi}\in\mathrm{BV}(\Omega,T)$ , we invoke the limsup property to find a sequence $\{\bm{\psi}^{\varepsilon}\}_{\varepsilon>0}$ such that $\bm{\psi}^{\varepsilon}\to\bm{\psi}$ strongly in $L^{1}(\Omega,{\mathbb{R}}^{L})$ and $\limsup_{\varepsilon\to 0}{\mathbb{E}}_{\varepsilon}(\bm{\psi}^{\varepsilon})\leq{\mathbb{E}}_{0}(\bm{\psi})<\infty$ . Continuity of $G$ implies $G(\bm{\psi}^{\varepsilon})\to G(\bm{\psi})$ as $\varepsilon\to 0$ , leading to

\limsup_{\varepsilon\to 0}J_{\rm red}^{\varepsilon}(\bm{\psi}^{\varepsilon})\leq J_{\rm red}^{0}(\bm{\psi})<\infty.

As ${\bm{\varphi}}^{\varepsilon}$ minimises $J_{\rm red}^{\varepsilon}$ , we see that

\limsup_{\varepsilon\to 0}J_{\rm red}^{\varepsilon}({\bm{\varphi}}^{\varepsilon})\leq\limsup_{\varepsilon\to 0}J_{\rm red}^{\varepsilon}(\bm{\psi}^{\varepsilon})\leq J_{\rm red}^{0}(\bm{\psi})<\infty.

By the non-negativity of $G$ , the above estimate implies $\sup_{\varepsilon\in(0,1]}{\mathbb{E}}_{\varepsilon}({\bm{\varphi}}^{\varepsilon})<\infty$ , and by the compactness property we deduce that there exists a limit function ${\bm{\varphi}}_{0}\in\mathrm{BV}(\Omega,T)$ such that ${\bm{\varphi}}^{\varepsilon}\to{\bm{\varphi}}_{0}$ strongly in $L^{1}(\Omega,{\mathbb{R}}^{L})$ along a non-relabelled subsequence. Continuity of $G$ then gives $G({\bm{\varphi}}^{\varepsilon})\to G({\bm{\varphi}}_{0})$ and invoking the liminf property we subsequently infer that

J_{\rm red}^{0}({\bm{\varphi}}_{0})\leq\liminf_{\varepsilon\to 0}J_{\rm red}^{\varepsilon}({\bm{\varphi}}^{\varepsilon})\leq J_{\rm red}^{0}(\bm{\psi}).

As $\bm{\psi}$ is arbitrary, this shows that ${\bm{\varphi}}_{0}$ is a minimiser of $J_{\rm red}^{0}$ . We now return to the beginning of the proof and consider using the limsup property to construct a sequence $\{\bm{\varphi}^{\varepsilon}\}_{\varepsilon>0}$ that converges strongly to ${\bm{\varphi}}_{0}$ in $L^{1}(\Omega,{\mathbb{R}}^{d})$ . Then, following a similar argument we arrive at

J_{\rm red}^{0}({\bm{\varphi}}_{0})\leq\liminf_{\varepsilon\to 0}J_{\rm red}^{\varepsilon}({\bm{\varphi}}^{\varepsilon})\leq\limsup_{\varepsilon\to 0}J_{\rm red}^{\varepsilon}({\bm{\varphi}}^{\varepsilon})\leq J_{\rm red}^{0}({\bm{\varphi}}_{0}),

which provides the claimed assertion $\lim_{\varepsilon\to 0}J_{\rm red}^{\varepsilon}({\bm{\varphi}}^{\varepsilon})=J_{\rm red}^{0}({\bm{\varphi}}_{0})$ . ∎

5.2 Formally matched asymptotic expansions

We turn our attention towards the optimality condition (4.21) and study its sharp interface limit $\varepsilon\to 0$ with the method of formally matched asymptotic expansions, where we assume the functions ${\bm{\varphi}}^{\varepsilon}$ , ${\overline{{\bf u}}}^{\varepsilon}$ , ${\widehat{{\bf u}}}^{\varepsilon}$ , $\overline{{\bf p}}^{\varepsilon}$ , and $\widehat{{\bf q}}^{\varepsilon}$ admit asymptotic expansions in powers of $\varepsilon$ . From Lemma 5.1 we saw that ${\bm{\varphi}}^{\varepsilon}$ converges to a function ${\bm{\varphi}}_{0}\in\mathrm{BV}(\Omega,T)$ as $\varepsilon\to 0$ , and thus for $0<\varepsilon<1$ , we expect ${\bm{\varphi}}^{\varepsilon}$ to change its values rapidly on a length scale proportional to $\varepsilon$ . This inspires us to consider two asymptotic expansions of $({\bm{\varphi}}^{\varepsilon},{\overline{{\bf u}}}^{\varepsilon},{\widehat{{\bf u}}}^{\varepsilon},\overline{{\bf p}}^{\varepsilon},\widehat{{\bf q}}^{\varepsilon})$ in the bulk and interfacial regions (to be defined below), and the procedure is to match these expansions in an intermediate region to deduce the expected equations in the sharp interface limit. We follow the ideas in [13] that treats a similar system of equations, and refer the reader to, e.g., [18, 25, 32, 33] for more details regarding the methodology.

Recalling $T=\Psi^{-1}(0)=\{\bm{e}_{1},\dots,\bm{e}_{L}\}$ as the set of corners of the Gibbs simplex $\Delta^{L}$ , we partition the domain $\Omega$ into regions $\Omega^{i}$ , $i=1,\dots,L$ , where $\Omega^{i}=\{\bm{x}\in\Omega\,:\,{\bm{\varphi}}_{0}(\bm{x})=\bm{e}_{i}\}$ . Then, we assume the functions ${\bm{\varphi}}^{\varepsilon}$ , ${\overline{{\bf u}}}^{\varepsilon}$ , ${\widehat{{\bf u}}}^{\varepsilon}$ , $\overline{{\bf p}}^{\varepsilon}$ , and $\widehat{{\bf q}}^{\varepsilon}$ are sufficiently smooth and admit the following asymptotic expansion in $\varepsilon$ :

	$\displaystyle{\bm{\varphi}}^{\varepsilon}(\bm{x})=\sum_{k=0}^{\infty}\varepsilon^{k}{\bm{\varphi}}_{k}(\bm{x}),\quad{\overline{{\bf u}}}^{\varepsilon}(\bm{x})=\sum_{k=0}^{\infty}\varepsilon^{k}{\overline{{\bf u}}}_{k}(\bm{x}),\quad{\widehat{{\bf u}}}^{\varepsilon}(\bm{x})=\sum_{k=0}^{\infty}\varepsilon^{k}{\widehat{{\bf u}}}_{k}(\bm{x}),$
	$\displaystyle\overline{{\bf p}}^{\varepsilon}(\bm{x})=\sum_{k=0}^{\infty}\varepsilon^{k}\overline{{\bf p}}_{k}(\bm{x}),\quad\widehat{{\bf q}}^{\varepsilon}(\bm{x})=\sum_{k=0}^{\infty}\varepsilon^{k}\widehat{{\bf q}}_{k}(\bm{x}),$

for all points $\bm{x}\in\Omega$ away from the interfaces $\Gamma_{ij}=\partial\Omega_{i}\cap\partial\Omega_{j}$ for $i,j\in\{1,\dots,L\}$ , $i\neq j$ . This is known as the outer expansion. Furthermore, we assume that

{\bm{\varphi}}_{k}(\bm{x})\in T\Sigma^{L}:=\Big{\{}\bm{v}=(v_{1},\dots,v_{L})\in{\mathbb{R}}^{L}\,:\,\sum_{i=1}^{L}v_{i}=0\Big{\}},\quad k\geq 1,

where $T\Sigma^{L}$ is the tangent space of the affine hyperplane $\Sigma^{L}=\{\bm{v}\in{\mathbb{R}}^{L}\,:\,\sum_{i=1}^{L}v_{i}=1\}$ , so that by the above construction ${\bm{\varphi}}^{\varepsilon}(\bm{x})\in\Delta^{L}$ for $\varepsilon$ sufficiently small. We assume that there are constant elasticity tensors $\overline{{\mathbb{C}}}_{i}$ and $\widehat{{\mathbb{C}}}_{i}$ for $i=1,\dots,L$ , satisfying the standard symmetric conditions and are positive definite. Then, for ${\bm{\varphi}}=(\varphi_{1},\dots,\varphi_{L})$ such that ${\bm{\varphi}}(\bm{x})\in\Delta^{L}$ , we consider

\overline{{\mathbb{C}}}({\bm{\varphi}})=\sum_{i=1}^{L}\overline{{\mathbb{C}}}_{i}\varphi_{i},\quad\widehat{{\mathbb{C}}}({\bm{\varphi}})=\sum_{i=1}^{L}\widehat{{\mathbb{C}}}_{i}\varphi_{i}.

Then, substituting the outer expansions into the state systems (2.1), (2.2) and the adjoint systems (4.17) and (4.19), to leading order we obtain the following equations for $i=1,\dots,L$ :

		$\displaystyle\begin{cases}-\mathop{\rm div}\nolimits\big{(}\overline{{\mathbb{C}}}_{i}\mathcal{E}({\overline{{\bf u}}}_{0})\big{)}=\overline{{\bf F}}&\text{ in }\Omega^{i},\\ {\overline{{\bf u}}}_{0}=\overline{\bm{U}}&\text{ on }\overline{\Gamma}_{D}\cap\partial\Omega^{i},\\ (\overline{{\mathbb{C}}}_{i}\mathcal{E}({\overline{{\bf u}}}_{0})){\bf n}=\overline{{\bf g}}&\text{ on }\overline{\Gamma}_{N}\cap\partial\Omega^{i},\end{cases}$		(5.1)
		$\displaystyle\begin{cases}-\mathop{\rm div}\nolimits\big{(}\widehat{{\mathbb{C}}}_{i}(\mathcal{E}({\widehat{{\bf u}}}_{0})-\chi_{i}\mathcal{E}({\overline{{\bf u}}}_{0}))\big{)}=\widehat{{\bf F}}&\text{ in }\Omega^{i},\\ {\widehat{{\bf u}}}_{0}=\widehat{\bm{U}}&\text{ on }\widehat{\Gamma}_{D}\cap\partial\Omega^{i},\\ (\widehat{{\mathbb{C}}}_{i}(\mathcal{E}({\widehat{{\bf u}}}_{0})-\chi_{i}\mathcal{E}({\overline{{\bf u}}}_{0}))){\bf n}=\widehat{{\bf g}}&\text{ on }\widehat{\Gamma}_{N}\cap\partial\Omega^{i},\end{cases}$		(5.1)

where $\chi_{i}:=\chi(\bm{e}_{i})$ , and

		$\displaystyle\begin{cases}-\mathop{\rm div}\nolimits\big{(}\widehat{{\mathbb{C}}}_{i}\mathcal{E}(\widehat{{\bf q}}_{0})\big{)}=\bm{0}&\text{ in }\Omega^{i},\\ \widehat{{\bf q}}_{0}=\bm{0}&\text{ on }\widehat{\Gamma}_{D}\cap\partial\Omega^{i},\\ (\widehat{{\mathbb{C}}}_{i}\mathcal{E}(\widehat{{\bf q}}_{0})){\bf n}=W({\widehat{{\bf u}}}_{0}-{\bf u}^{\rm tar})\mathrm{\mathcal{X}}_{\Gamma^{\rm tar}}&\text{ on }\widehat{\Gamma}_{N}\cap\partial\Omega^{i},\end{cases}$		(5.2)
		$\displaystyle\begin{cases}-\mathop{\rm div}\nolimits\big{(}\overline{{\mathbb{C}}}_{i}\mathcal{E}(\overline{{\bf p}}_{0})-\widehat{{\mathbb{C}}}_{i}\chi_{i}\mathcal{E}(\widehat{{\bf q}}_{0}))\big{)}=\bm{0}&\text{ in }\Omega^{i},\\ \overline{{\bf p}}_{0}=\bm{0}&\text{ on }\overline{\Gamma}_{D}\cap\partial\Omega^{i},\\ (\overline{{\mathbb{C}}}_{i}\mathcal{E}(\overline{{\bf p}}_{0})-\widehat{{\mathbb{C}}}_{i}\chi_{i}\mathcal{E}(\widehat{{\bf q}}_{0})){\bf n}=\bm{0}&\text{ on }\overline{\Gamma}_{N}\cap\partial\Omega^{i}.\end{cases}$		(5.2)

It then remains to furnish the above with boundary conditions for $({\overline{{\bf u}}}_{0},{\widehat{{\bf u}}}_{0},\overline{{\bf p}}_{0},\widehat{{\bf q}}_{0})$ on the interfaces $\Gamma_{ij}$ for $i,j\in\{1,\dots,L\}$ , $i<j$ , which we assume are smooth hypersurfaces that can be obtained in the limit $\varepsilon\to 0$ . These boundary conditions can be inferred with the help of a corresponding inner expansion for $({\bm{\varphi}}^{\varepsilon},{\overline{{\bf u}}}^{\varepsilon},{\widehat{{\bf u}}}^{\varepsilon},\overline{{\bf p}}^{\varepsilon},\widehat{{\bf q}}^{\varepsilon})$ in the interfacial regions bordering two bulk regions $\Omega^{i}$ and $\Omega^{j}$ . To this end, we focus on one particular interface $\Gamma_{ij}$ and introduce a change of coordinates with the help of a local parameterisation $\bm{\gamma}:U\subset{\mathbb{R}}^{d-1}\to{\mathbb{R}}^{d}$ of $\Gamma_{ij}$ , where $U$ is a spatial parameter domain.

Close to $\bm{\gamma}(U)$ , consider the signed distance function $d$ such that $d(\bm{x})>0$ if $\bm{x}\in\Omega^{j}$ and $d(\bm{x})<0$ if $\bm{x}\in\Omega^{i}$ , so that the unit normal $\bm{\nu}$ to $\Gamma_{ij}$ points from $\Omega^{i}$ to $\Omega^{j}$ . Introducing the rescaled signed distance $z=\frac{d}{\varepsilon}$ , a local parameterization of $\bm{x}\in{\mathbb{R}}^{d}$ close to ${\bm{\gamma}}(U)$ can be given as

\bm{x}={\cal G}^{\varepsilon}(\bm{s},z)=\bm{\gamma}(\bm{s})+\varepsilon z\bm{\nu}(\bm{s}),\quad\bm{s}\in U\subset{\mathbb{R}}^{d-1},\,z\in{\mathbb{R}}.

This representation allows us to infer the following expansions for gradients, divergences and Laplacians [34]:

	$\displaystyle\nabla_{x}b$	$\displaystyle=\frac{1}{\varepsilon}\partial_{z}\widehat{b}\bm{\nu}+\nabla_{\Gamma_{ij}}\widehat{b}+\mathcal{O}(\varepsilon),\quad\mathop{\rm div}\nolimits_{x}\bm{j}=\frac{1}{\varepsilon}\partial_{z}\widehat{\bm{j}}\cdot\bm{\nu}+\mathop{\rm div}\nolimits_{\Gamma_{ij}}\widehat{\bm{j}}+\mathcal{O}(\varepsilon),$
	$\displaystyle\nabla_{x}\bm{j}$	$\displaystyle=\frac{1}{\varepsilon}\partial_{z}\widehat{\bm{j}}\otimes\bm{\nu}+\nabla_{\Gamma_{ij}}\widehat{\bm{j}}+\mathcal{O}(\varepsilon),\quad\Delta_{x}b=\frac{1}{\varepsilon^{2}}\partial_{zz}\widehat{b}-\frac{1}{\varepsilon}\kappa_{\Gamma_{ij}}\partial_{z}\widehat{b}+\mathcal{O}(1),$

for scalar functions $b(\bm{x})=\widehat{b}(\bm{s}(\bm{x}),z(\bm{x}))$ and vector functions $\bm{j}(\bm{x})=\widehat{\bm{j}}(\bm{s}(\bm{x}),z(\bm{x}))$ , along with the curvature $\kappa_{\Gamma_{ij}}$ of $\Gamma_{ij}$ , the surface gradient operator $\nabla_{\Gamma_{ij}}$ , and the surface divergence operator $\mathop{\rm div}\nolimits_{\Gamma_{ij}}$ on $\Gamma_{ij}$ . Then, for points close by $\Gamma_{ij}$ , we assume an inner expansion of the form

	$\displaystyle{\bm{\varphi}}^{\varepsilon}(\bm{x})=\sum_{k=0}^{\infty}\varepsilon^{k}\bm{\Phi}_{k}(\bm{s},z),\quad{\overline{{\bf u}}}^{\varepsilon}(\bm{x})=\sum_{k=0}^{\infty}\varepsilon^{k}\overline{\bm{U}}_{k}(\bm{s},z),\quad{\widehat{{\bf u}}}^{\varepsilon}(\bm{x})=\sum_{k=0}^{\infty}\varepsilon^{k}\widehat{\bm{U}}_{k}(\bm{s},z),$
	$\displaystyle\overline{{\bf p}}^{\varepsilon}(\bm{x})=\sum_{k=0}^{\infty}\varepsilon^{k}\overline{\bm{P}}_{k}(\bm{s},z),\quad\widehat{{\bf q}}^{\varepsilon}(\bm{x})=\sum_{k=0}^{\infty}\varepsilon^{k}\widehat{\bm{Q}}_{k}(\bm{s},z).$

Lastly, we assume in a tubular neighborhood of $\Gamma_{ij}$ the outer expansions and the inner expansions hold simultaneously. Within this intermediate region certain matching conditions relating the outer expansions to the inner expansions must hold. For a scalar function $b(\bm{x})$ admitting an outer expansion $\sum_{k=0}^{\infty}\varepsilon^{k}b_{k}(\bm{x})$ and an inner expansion $\sum_{k=0}^{\infty}\varepsilon^{k}B_{k}(\bm{s},z)$ , it holds that (see [34, Appendix D])

	$\displaystyle B_{0}(\bm{s},z)$	$\displaystyle\to\begin{cases}\lim_{\delta\searrow 0}b_{0}(\bm{x}+\delta\bm{\nu}(\bm{x}))=:b_{0}^{+}(\bm{x})&\text{ for }z\to+\infty,\\ \lim_{\delta\searrow 0}b_{0}(\bm{x}-\delta\bm{\nu}(\bm{x}))=:b_{0}^{-}(\bm{x})&\text{ for }z\to-\infty,\end{cases}$
	$\displaystyle\partial_{z}B_{0}(\bm{s},z)$	$\displaystyle\to 0\text{ as }z\to\pm\infty,$
	$\displaystyle\partial_{z}B_{1}(\bm{s},z)$	$\displaystyle\to\begin{cases}\lim_{\delta\searrow 0}(\nabla b_{0})(\bm{x}+\delta\bm{\nu}(\bm{x}))\cdot\bm{\nu}(\bm{x})=:\nabla b_{0}^{+}\cdot\bm{\nu}&\text{ for }z\to+\infty,\\ \lim_{\delta\searrow 0}(\nabla b_{0})(\bm{x}-\delta\bm{\nu}(\bm{x}))\cdot\bm{\nu}(\bm{x})=:\nabla b_{0}^{-}\cdot\bm{\nu}&\text{ for }z\to-\infty,\end{cases}$

for $\bm{x}\in\Gamma_{ij}$ . Consequently, we denote the jump of a quantity $b$ across $\Gamma_{ij}$ as

[b]_{-}^{+}:=\lim_{\delta\searrow 0}b(\bm{x}+\delta\bm{\nu}(\bm{x}))-\lim_{\delta\searrow 0}b(\bm{x}-\delta\bm{\nu}(\bm{x}))\quad\text{ for }\bm{x}\in\Gamma_{ij}.

Note that the above matching conditions also apply to vectorial functions. We introduce the orthogonal projection

\bm{P}_{T\Sigma}:{\mathbb{R}}^{L}\to T\Sigma^{L},\quad\bm{P}_{T\Sigma}{\bm{\varphi}}={\bm{\varphi}}-\Big{(}\frac{1}{L}\sum_{i=1}^{L}\varphi_{i}\Big{)}\bm{1}

where $\bm{1}:=(1,\dots,1)^{\top}$ , so that the optimality condition (4.21) can be expressed in the following strong form

		$\displaystyle-2\gamma\varepsilon\Delta{\bm{\varphi}}^{\varepsilon}+\bm{P}_{T\Sigma}\Big{(}\frac{\gamma}{\varepsilon}\widetilde{\Psi}_{,{\bm{\varphi}}}({\bm{\varphi}}^{\varepsilon})+\widehat{{\mathbb{C}}}({\bm{\varphi}}^{\varepsilon})\chi^{\prime}({\bm{\varphi}}^{\varepsilon})\mathcal{E}({\overline{{\bf u}}}^{\varepsilon}):\mathcal{E}(\widehat{{\bf q}}^{\varepsilon})\Big{)}$		(5.3)
		$\displaystyle\quad-\bm{P}_{T\Sigma}\Big{(}\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}^{\varepsilon})(\mathcal{E}({\widehat{{\bf u}}}^{\varepsilon})-\chi({\bm{\varphi}}^{\varepsilon})\mathcal{E}({\overline{{\bf u}}}^{\varepsilon})):\mathcal{E}(\widehat{{\bf q}}^{\varepsilon})+\overline{{\mathbb{C}}}^{\prime}({\bm{\varphi}}^{\varepsilon})\mathcal{E}({\overline{{\bf u}}}^{\varepsilon}):\mathcal{E}(\overline{\bm{p}}^{\varepsilon})\Big{)}=\bm{0}.$		(5.3)

We substitute the inner expansions into the equations (2.1a), (2.2a), (4.17a), (4.19a), and (5.3) and collect terms of the same order. Then, we perform some computations to deduce the boundary conditions posed on $\Gamma_{ij}$ . As the subsequent analysis is similar to that performed in [13, 32], we omit most of the straightforward details. In the sequel we use the notation

(\bm{B})^{\mathrm{sym}}=\frac{1}{2}(\bm{B}+\bm{B}^{\top}),\quad\mathcal{E}_{\bf X}:=(\partial_{z}{\bf X}_{1}\otimes\bm{\nu}+\nabla_{\Gamma_{ij}}{\bf X}_{0})^{\mathrm{sym}},

for second order tensors $\bm{B}$ and for ${\bf X}\in\{\overline{\bm{U}},\widehat{\bm{U}},\overline{\bm{P}},\widehat{\bm{Q}}\}$ .

To leading order $\mathcal{O}(\frac{1}{\varepsilon^{2}})$ , equations (2.1a) and (4.17a) yield

\partial_{z}\Big{(}\overline{{\mathbb{C}}}(\bm{\Phi}_{0})(\partial_{z}\overline{\bm{U}}_{0}\otimes\bm{\nu})^{\mathrm{sym}}\bm{\nu}\Big{)}=\bm{0},\quad\partial_{z}\Big{(}\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})(\partial_{z}\widehat{\bm{Q}}_{0}\otimes\bm{\nu})^{\mathrm{sym}}\bm{\nu}\Big{)}=\bm{0}.

Multiplying by $\overline{\bm{U}}_{0}$ and $\widehat{\bm{Q}}_{0}$ , respectively, integrating over $z\in{\mathbb{R}}$ , integrating by parts and applying the matching conditions allow us to deduce that $\partial_{z}\overline{\bm{U}}_{0}=\partial_{z}\widehat{\bm{Q}}_{0}=\bm{0}$ , and hence both $\overline{\bm{U}}_{0}$ and $\widehat{\bm{Q}}_{0}$ are constant in $z$ . Applying matching conditions we infer that

[{\overline{{\bf u}}}_{0}]_{i}^{j}=[\widehat{\bm{q}}_{0}]_{i}^{j}=\bm{0}\quad\text{ on }\Gamma_{ij}.

Then, to leading order $\mathcal{O}(\frac{1}{\varepsilon^{2}})$ , equations (2.2a) and (4.19a) yield

\partial_{z}\Big{(}\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})(\partial_{z}\widehat{\bm{U}}_{0}\otimes\bm{\nu})^{\mathrm{sym}}\bm{\nu}\Big{)}=\bm{0},\quad\partial_{z}\Big{(}\overline{{\mathbb{C}}}(\bm{\Phi}_{0})(\partial_{z}\overline{\bm{P}}_{0}\otimes\bm{\nu})^{\mathrm{sym}}\bm{\nu}\Big{)}=\bm{0},

on account of the fact that $\partial_{z}\overline{\bm{U}}_{0}=\partial_{z}\widehat{\bm{Q}}_{0}=\bm{0}$ . Hence, we also obtain

[{\widehat{{\bf u}}}_{0}]_{i}^{j}=[\overline{\bm{p}}_{0}]_{i}^{j}=\bm{0}\quad\text{ on }\Gamma_{ij}.

To first order $\mathcal{O}(\frac{1}{\varepsilon})$ , we get from (2.1a) and (4.17a) that

\displaystyle\partial_{z}\Big{(}\overline{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\overline{\bm{U}}}\bm{\nu}\Big{)}=\bm{0},\quad\partial_{z}\Big{(}\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\widehat{\bm{Q}}}\bm{\nu}\Big{)}=\bm{0}.

(5.4)

Integrating over $z\in{\mathbb{R}}$ and using the matching conditions yields

[\overline{{\mathbb{C}}}\mathcal{E}({\overline{{\bf u}}}_{0})\bm{\nu}]_{i}^{j}=\bm{0},\quad[\widehat{{\mathbb{C}}}\mathcal{E}(\widehat{\bm{q}}_{0})\bm{\nu}]_{i}^{j}=\bm{0}\quad\text{ on }\Gamma_{ij}.

Similarly, from equations (2.2a) and (4.19a), we obtain to first order $\mathcal{O}(\frac{1}{\varepsilon})$ that

\displaystyle\partial_{z}\Big{(}\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})(\mathcal{E}_{\widehat{\bm{U}}}-\chi(\bm{\Phi_{0}})\mathcal{E}_{\overline{\bm{U}}})\bm{\nu}\Big{)}=\bm{0},\quad\partial_{z}\Big{(}\overline{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\overline{\bm{P}}}\bm{\nu}-\chi(\bm{\Phi}_{0})\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\widehat{\bm{Q}}}\bm{\nu}\Big{)}=\bm{0}.

(5.5)

Integrating over $z\in{\mathbb{R}}$ and applying the matching conditions, we obtain

[\widehat{{\mathbb{C}}}(\mathcal{E}({\widehat{{\bf u}}}_{0})-\chi({\bm{\varphi}}_{0})\mathcal{E}({\overline{{\bf u}}}_{0}))]_{i}^{j}\bm{\nu}=\bm{0},\quad[\overline{{\mathbb{C}}}\mathcal{E}(\overline{\bm{p}}_{0})-\chi({\bm{\varphi}}_{0})\widehat{{\mathbb{C}}}\mathcal{E}(\widehat{\bm{q}}_{0})]_{i}^{j}\bm{\nu}=\bm{0}\quad\text{ on }\Gamma_{ij}.

Turning now to the optimality condition (5.3), we use the fact that $\partial_{z}\overline{\bm{U}}_{0}=\partial_{z}\widehat{\bm{U}}_{0}=\partial_{z}\overline{\bm{P}}_{0}=\partial_{z}\widehat{\bm{Q}}_{0}=\bm{0}$ to see that the elasticity terms do not contribute to leading order $\mathcal{O}(\frac{1}{\varepsilon^{2}})$ and first order $\mathcal{O}(\frac{1}{\varepsilon})$ . Hence, to first order $\mathcal{O}(\frac{1}{\varepsilon})$ we obtain from (5.3) that

\displaystyle 2\partial_{zz}\bm{\Phi}_{0}-\bm{P}_{T\Sigma}\Big{(}\widetilde{\Psi}_{,{\bm{\varphi}}}(\bm{\Phi}_{0})\Big{)}=\bm{0}.

Following [18], $\bm{\Phi}_{0}$ can be chosen independent of $s$ and as the solution to the above ordinary differential system such that $\lim_{z\to-\infty}\bm{\Phi}_{0}(z)=\bm{e}_{i}$ and $\lim_{z\to+\infty}\bm{\Phi}_{0}(z)=\bm{e}_{j}$ . To the next order $\mathcal{O}(1)$ , we obtain

		$\displaystyle-2\gamma\partial_{zz}\bm{\Phi}_{1}+\gamma\bm{P}_{T\Sigma}\Big{(}\widetilde{\Psi}_{,{\bm{\varphi}}{\bm{\varphi}}}(\bm{\Phi}_{0})\bm{\Phi}_{1}\Big{)}+2\gamma\kappa_{\Gamma_{ij}}\partial_{z}\bm{\Phi}_{0}+\bm{P}_{T\Sigma}\Big{(}\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\chi^{\prime}(\bm{\Phi}_{0})\mathcal{E}_{\overline{\bm{U}}}:\mathcal{E}_{\widehat{\bm{Q}}}\Big{)}$		(5.6)
		$\displaystyle\quad-\bm{P}_{T\Sigma}\Big{(}\widehat{{\mathbb{C}}}^{\prime}(\bm{\Phi}_{0})(\mathcal{E}_{\widehat{\bm{U}}}-\chi(\bm{\Phi}_{0})\mathcal{E}_{\overline{\bm{U}}}):\mathcal{E}_{\widehat{\bm{Q}}}+\overline{{\mathbb{C}}}^{\prime}(\bm{\Phi}_{0})\mathcal{E}_{\overline{\bm{U}}}:\mathcal{E}_{\overline{\bm{P}}}\Big{)}=\bm{0},$		(5.6)

where $\widetilde{\Psi}_{,{\bm{\varphi}}{\bm{\varphi}}}$ denotes the Hessian matrix of $\widetilde{\Psi}$ , and we have used that $\partial_{z}\bm{\Phi}_{0}\in T\Sigma^{L}$ . Note that by the fact that $\partial_{z}{\bf X}_{0}=\bm{0}$ for ${\bf X}\in\{\overline{\bm{U}},\widehat{\bm{U}},\overline{\bm{P}},\widehat{\bm{Q}}\}$ , and by the symmetry of the elasticity tensors ${\mathbb{C}}_{ijkl}={\mathbb{C}}_{jikl}$ , we have the relations

\displaystyle\partial_{z}\mathcal{E}_{\bf X}=(\partial_{zz}{\bf X}_{1}\otimes\bm{\nu})^{\mathrm{sym}},\quad{\mathbb{C}}\mathcal{E}_{\bf Y}:\partial_{z}\mathcal{E}_{\bf X}=(({\mathbb{C}}\mathcal{E}_{\bf Y})\bm{\nu})\cdot\partial_{zz}{\bf X}_{1},

(5.7)

for any ${\bf X},{\bf Y}\in\{\overline{\bm{U}},\widehat{\bm{U}},\overline{\bm{P}},\widehat{\bm{Q}}\}$ . To obtain a solution $\bm{\Phi}_{1}$ , a solvability condition has to hold, which can be derived by multiplying (5.6) with $\partial_{z}\bm{\Phi}_{0}$ and integrating over $z$ . Using the relations $\bm{P}_{T\Sigma}(\partial_{z}\bm{\Phi}_{0})=\partial_{z}\bm{\Phi}_{0}$ , $\bm{P}_{T\Sigma}(\partial_{z}\bm{\Phi}_{1})=\partial_{z}\bm{\Phi}_{1}$ , after integrating by parts and applying the matching conditions, we obtain

	$\displaystyle 0$	$\displaystyle=\int_{-\infty}^{\infty}\gamma\underbrace{(-2\partial_{zz}\bm{\Phi}_{0}+\bm{P}_{T\Sigma}\Big{(}\widetilde{\Psi}_{,{\bm{\varphi}}}(\bm{\Phi}_{0}))\Big{)}}_{=\bm{0}}\cdot\,\partial_{z}\bm{\Phi}_{1}dz+2\gamma\kappa_{\Gamma_{ij}}\int_{-\infty}^{\infty}\|\partial_{z}\bm{\Phi}_{0}\|^{2}dz$		(5.8)
		$\displaystyle\quad+\int_{-\infty}^{\infty}\partial_{z}[\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\chi(\bm{\Phi}_{0})]\mathcal{E}_{\overline{\bm{U}}}:\mathcal{E}_{\widehat{\bm{Q}}}-\partial_{z}\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\widehat{\bm{U}}}:\mathcal{E}_{\widehat{\bm{Q}}}-\partial_{z}\overline{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\overline{\bm{U}}}:\mathcal{E}_{\overline{\bm{P}}}\;dz.$		(5.8)

We employ the identities obtained from (5.4), (5.5) and (5.7) to obtain that

\displaystyle-\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\widehat{\bm{U}}}:\partial_{z}\mathcal{E}_{\widehat{\bm{Q}}}+\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\widehat{\bm{U}}}\bm{\nu}\cdot\partial_{zz}\widehat{\bm{Q}}_{1}=\bm{0},\quad\partial_{z}(\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\widehat{\bm{Q}}}\bm{\nu})=\bm{0},

as well as $\partial_{z}\bm{\nu}=\bm{0}$ to see that

	$\displaystyle-\partial_{z}\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\widehat{\bm{U}}}:\mathcal{E}_{\widehat{\bm{Q}}}$
	$\displaystyle\quad=-\partial_{z}\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\widehat{\bm{U}}}:\mathcal{E}_{\widehat{\bm{Q}}}-\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\widehat{\bm{U}}}:\partial_{z}\mathcal{E}_{\widehat{\bm{Q}}}-\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\widehat{\bm{Q}}}:\partial_{z}\mathcal{E}_{\widehat{\bm{U}}}$
	$\displaystyle\qquad+\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\widehat{\bm{U}}}\bm{\nu}\cdot\partial_{zz}\widehat{\bm{Q}}_{1}+\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\widehat{\bm{Q}}}\bm{\nu}\cdot\partial_{zz}\widehat{\bm{U}}_{1}$
	$\displaystyle\qquad+\partial_{z}\Big{(}\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\widehat{\bm{Q}}}\Big{)}\bm{\nu}\cdot\partial_{z}\widehat{\bm{U}}_{1}+\partial_{z}\Big{(}\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})(\mathcal{E}_{\widehat{\bm{U}}}-\chi(\bm{\Phi}_{0})\mathcal{E}_{\overline{\bm{U}}})\Big{)}\bm{\nu}\cdot\partial_{z}\widehat{\bm{Q}}_{1}$
	$\displaystyle\quad=-\partial_{z}\Big{(}\widehat{{\mathbb{C}}}(\bm{\Phi_{0}})\mathcal{E}_{\widehat{\bm{U}}}:\mathcal{E}_{\widehat{\bm{Q}}}-\widehat{{\mathbb{C}}}(\bm{\Phi_{0}})\mathcal{E}_{\widehat{\bm{Q}}}\bm{\nu}\cdot\partial_{z}\widehat{\bm{U}}_{1}-\widehat{{\mathbb{C}}}(\bm{\Phi_{0}})\mathcal{E}_{\widehat{\bm{U}}}\bm{\nu}\cdot\partial_{z}\widehat{\bm{Q}}_{1}\big{)}\Big{)}$
	$\displaystyle\qquad-\partial_{z}\Big{(}\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\chi(\bm{\Phi}_{0})\mathcal{E}_{\overline{\bm{U}}}\Big{)}\bm{\nu}\cdot\partial_{z}\widehat{\bm{Q}}_{1}.$

Via a similar calculation we infer

	$\displaystyle-\partial_{z}\overline{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\overline{\bm{U}}}:\mathcal{E}_{\overline{\bm{P}}}$
	$\displaystyle\quad=-\partial_{z}\Big{(}\overline{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\overline{\bm{U}}}:\mathcal{E}_{\overline{\bm{P}}}-\overline{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\overline{\bm{P}}}\bm{\nu}\cdot\partial_{z}\overline{\bm{U}}_{1}-\overline{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\overline{\bm{U}}}\bm{\nu}\cdot\partial_{z}\overline{\bm{P}}_{1}\Big{)}$
	$\displaystyle\qquad-\partial_{z}\Big{(}\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\chi(\bm{\Phi}_{0})\mathcal{E}_{\widehat{\bm{Q}}}\Big{)}\bm{\nu}\cdot\partial_{z}\overline{\bm{U}}_{1},$
	$\displaystyle\partial_{z}\big{[}\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\chi(\bm{\Phi}_{0})\big{]}\mathcal{E}_{\overline{\bm{U}}}:\mathcal{E}_{\widehat{\bm{Q}}}$
	$\displaystyle\quad=\partial_{z}\Big{(}\chi(\bm{\Phi}_{0})\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\overline{\bm{U}}}:\mathcal{E}_{\widehat{\bm{Q}}}\Big{)}$
	$\displaystyle\qquad-\chi(\bm{\Phi}_{0})\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\overline{\bm{U}}}\bm{\nu}\cdot\partial_{zz}\widehat{\bm{Q}}_{1}-\chi(\bm{\Phi}_{0})\widehat{{\mathbb{C}}}(\bm{\Phi}_{0})\mathcal{E}_{\widehat{\bm{Q}}}\bm{\nu}\cdot\partial_{zz}\overline{\bm{U}}_{1},$

and thus, setting $b_{ij}:=\int_{-\infty}^{\infty}2|\partial_{z}\bm{\Phi}_{0}|^{2}dz$ , and applying matching conditions for $\mathcal{E}_{\bf X}$ and $\partial_{z}{\bf X}_{1}$ with ${\bf X}\in\{\overline{\bm{U}},\widehat{\bm{U}},\overline{\bm{P}},\widehat{\bm{Q}}\}$ , we obtain from (5.8) the solvability condition

$\displaystyle 0$	$\displaystyle=\gamma\kappa_{\Gamma_{ij}}b_{ij}+[\chi({\bm{\varphi}}_{0})\widehat{{\mathbb{C}}}\mathcal{E}({\overline{{\bf u}}}_{0}):\mathcal{E}(\widehat{\bm{q}}_{0})]_{i}^{j}-[\widehat{{\mathbb{C}}}\mathcal{E}({\widehat{{\bf u}}}_{0}):\mathcal{E}(\widehat{\bm{q}}_{0})]_{i}^{j}-[\overline{{\mathbb{C}}}\mathcal{E}({\overline{{\bf u}}}_{0}):\mathcal{E}(\overline{\bm{p}}_{0})]_{i}^{j}$	(5.9)
	$\displaystyle\quad+[\widehat{{\mathbb{C}}}\mathcal{E}(\widehat{\bm{q}}_{0})\bm{\nu}\cdot(\nabla{\widehat{{\bf u}}}_{0})\bm{\nu}+\widehat{{\mathbb{C}}}\mathcal{E}({\widehat{{\bf u}}}_{0})\bm{\nu}\cdot(\nabla\widehat{\bm{q}}_{0})\bm{\nu}]_{i}^{j}-[\chi({\bm{\varphi}}_{0})\widehat{{\mathbb{C}}}\mathcal{E}({\overline{{\bf u}}}_{0})\bm{\nu}\cdot(\nabla\widehat{\bm{q}}_{0})\bm{\nu}]_{i}^{j}$
	$\displaystyle\quad+[\overline{{\mathbb{C}}}\mathcal{E}(\overline{\bm{p}}_{0})\bm{\nu}\cdot(\nabla{\overline{{\bf u}}}_{0})\bm{\nu}+\overline{{\mathbb{C}}}\mathcal{E}({\overline{{\bf u}}}_{0})\bm{\nu}\cdot(\nabla\overline{\bm{p}}_{0})\bm{\nu}]_{i}^{j}-[\chi({\bm{\varphi}}_{0})\widehat{{\mathbb{C}}}\mathcal{E}(\widehat{\bm{q}}_{0})\bm{\nu}\cdot(\nabla{\overline{{\bf u}}}_{0})\bm{\nu}]_{i}^{j}$

that has to hold on $\Gamma_{ij}$ . Thus, the sharp interface limit consists of the equations (5.1) and (5.2) posed in $\Omega^{i}$ , $1\leq i\leq L$ , furnished by the boundary conditions (5.9) and

	$\displaystyle[{\widehat{{\bf u}}}_{0}]_{i}^{j}=\bm{0},\quad[\overline{\bm{p}}_{0}]_{i}^{j}=\bm{0},\quad[{\overline{{\bf u}}}_{0}]_{i}^{j}=\bm{0},\quad[\widehat{\bm{q}}_{0}]_{i}^{j}=\bm{0},$
	$\displaystyle[\overline{{\mathbb{C}}}\mathcal{E}({\overline{{\bf u}}}_{0})\bm{\nu}]_{i}^{j}=\bm{0},\quad[\widehat{{\mathbb{C}}}\mathcal{E}(\widehat{\bm{q}}_{0}))\bm{\nu}]_{i}^{j}=\bm{0},\quad[\overline{{\mathbb{C}}}\mathcal{E}({\overline{{\bf u}}}_{0})\bm{\nu}]_{i}^{j}=\bm{0},\quad[\widehat{{\mathbb{C}}}\mathcal{E}(\widehat{\bm{q}}_{0}))\bm{\nu}]_{i}^{j}=\bm{0}$

on $\Gamma_{ij}$ , $1\leq i<j\leq L$ .

Remark 5.1.

It is possible to consider the sharp interface limit near a triple junction where three regions meet. We refer to [13, 18, 53] for more details regarding the asymptotic analysis around a triple junction.

5.3 Rigorous convergence in the two-phase setting

In the two phase case $L=2$ , since ${\bm{\varphi}}=(\varphi_{1},\varphi_{2})\in\Delta^{2}$ , we may use the difference $\varphi:=\varphi_{2}-\varphi_{1}\in[-1,1]$ to encode the vector ${\bm{\varphi}}=(\frac{1}{2}(1-\varphi),\frac{1}{2}(1+\varphi))$ , so that $|\nabla{\bm{\varphi}}|^{2}=\frac{1}{2}|\nabla\varphi|^{2}$ . Hence, the problem (P) can be re-phrased in terms of the scalar function $\varphi$ ranging between $-1$ and $1$ , and it suffices to consider the following

	$\displaystyle(\bm{P}^{\varepsilon})\quad\min_{\varphi\in\mathcal{U}_{\rm ad}}J_{\rm red}^{\varepsilon}(\varphi)$	$\displaystyle=\frac{1}{2}\int_{\Gamma^{\mathrm{tar}}}W(\mathcal{S}(\varphi)-\bm{u}^{\mathrm{tar}})\cdot(\mathcal{S}(\varphi)-{\bf u}^{\rm tar})\,\mathrm{d}\mathcal{H}^{d-1}$
		$\displaystyle\quad+\gamma\int_{\Omega}\frac{\varepsilon}{2}\|\nabla\varphi\|^{2}+\frac{1}{\varepsilon}\Psi(\varphi)\,\mathrm{dx},$

where $\mathcal{U}_{\rm ad}=\{f\in H^{1}(\Omega)\,:\,f\in[-1,1]\text{ a.e.~{}in }\Omega\}$ , and, as no confusion can arise, we use the short-hand notations $\mathcal{S}(\varphi)$ and $\Psi(\varphi)$ for the functions $\mathcal{S}({\bm{\varphi}})$ and $\Psi({\bm{\varphi}})$ evaluated at ${\bm{\varphi}}=(\frac{1}{2}(1-\varphi),\frac{1}{2}(1+\varphi))$ . On recalling (A3), we hence assume that

\Psi(s)=\widetilde{\Psi}(s)+I_{[-1,1]}(s)\text{ for }s\in{\mathbb{R}},

(5.10)

for a $\widetilde{\Psi}\in C^{1,1}({\mathbb{R}})$ . To simplify the calculations, we consider $\overline{\bm{U}}=\widehat{\bm{U}}=\bm{0}$ (homogeneous Dirichlet data), $\overline{{\bf F}},\widehat{{\bf F}}\in H^{1}(\Omega,{\mathbb{R}}^{d})$ , $\overline{{\bf g}}\in H^{2}(\overline{\Gamma}_{N},{\mathbb{R}}^{d})$ and $\widehat{{\bf g}}\in H^{2}(\widehat{\Gamma}_{N},{\mathbb{R}}^{d})$ . We now consider deriving an alternative set of optimal conditions for a minimiser $\varphi^{*}\in\mathcal{U}_{\rm ad}$ based on geometric variations. To this end, we consider the following admissible transformations and their corresponding velocity fields.

Definition 5.1.

The space $\mathcal{V}_{\mathrm{ad}}$ of admissible velocity fields is defined as the set of all ${\bf V}\in C^{0}([-\tau,\tau]\times\overline{\Omega},{\mathbb{R}}^{d})$ , where $\tau>0$ is a fixed small constant, such that it holds:

•

${\bf V}(t,\cdot)\in C^{2}(\overline{\Omega},{\mathbb{R}}^{d})$ , and $\exists C>0$ such that $\|{\bf V}(\cdot,\bm{y})-{\bf V}(\cdot,\bm{x})\|_{C^{0}([-\tau,\tau],{\mathbb{R}}^{d})}\leq C\|\bm{x}-\bm{y}\|$ for all $\bm{x},\bm{y}\in\overline{\Omega}$ ;
•

${\bf V}(t,\bm{x})\cdot{\bf n}(\bm{x})=0$ for all $\bm{x}\in\partial\Omega$ ;
•

${\bf V}(t,\bm{x})=\bm{0}$ for all $\bm{x}\in\overline{\Gamma}_{D}\cup\widehat{\Gamma}_{D}$ .

Then, the space $\mathcal{T}_{\mathrm{ad}}$ of admissible transformations is defined as the set of solutions to the ordinary differential equations

\partial_{t}T_{t}(\bm{x})={\bf V}(t,T_{t}(\bm{x})),\quad T_{0}(\bm{x})=\bm{x}

with ${\bf V}\in\mathcal{V}_{\mathrm{ad}}$ .

Notice that by the second property it holds $\overline{T_{t}(\Omega)}=\overline{\Omega}$ for all $t\in[-\tau,\tau]$ . Let ${\bf V}\in\mathcal{V}_{\mathrm{ad}}$ be an admissible velocity field with corresponding transformation $T\in\mathcal{T}_{\mathrm{ad}}$ . For $\varphi\in\mathcal{U}_{\rm ad}$ we define $\varphi^{t}:=\varphi\circ T_{t}^{-1}$ , along with the unique solutions $({\overline{{\bf u}}}^{t},{\widehat{{\bf u}}}^{t})\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})\times\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ , where ${\overline{{\bf u}}}^{t}=\mathcal{S}_{1}^{2}(\varphi^{t})$ and ${\widehat{{\bf u}}}^{t}=\mathcal{S}(\varphi^{t})$ .

Setting $(\varphi_{0},{\overline{{\bf u}}}_{0},{\widehat{{\bf u}}}_{0})=(\varphi,{\overline{{\bf u}}},{\widehat{{\bf u}}})$ , by following a similar proof to [14, Lem. 25], we define for $\tau_{0}>0$ sufficiently small the function ${\bm{F}}_{1}:(-\tau_{0},\tau_{0})\times\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})\to(\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}))^{*}$ by

	$\displaystyle{\bm{F}}_{1}(t,{\bf u})[{\bf v}]$	$\displaystyle=\int_{\Omega}\big{(}\overline{{\mathbb{C}}}(\varphi)(\nabla T_{t}^{-1}\nabla{\bf u})^{\rm sym}:(\nabla T_{t}^{-1}\nabla{\bf v})^{\rm sym}-((\overline{{\bf F}}\circ T_{t})\cdot{\bf v})\big{)}\det(\nabla T_{t})\,\mathrm{dx}$
		$\displaystyle\quad-\int_{\overline{\Gamma}_{N}}((\overline{{\bf g}}\circ T_{t})\cdot{\bf v}\det(\nabla T_{t})\\|\nabla T_{t}^{-\top}{\bf n}\\|\,\mathrm{d}\mathcal{H}^{d-1}.$

Using a change of variables $y=T_{t}(\bm{x})$ , the relation

\nabla T_{t}^{-1}\nabla(\widetilde{{\bf u}}\circ T_{t})=(\nabla\widetilde{{\bf u}})\circ T_{t}\quad\text{ for }\widetilde{{\bf u}}:T_{t}(\Omega)\to{\mathbb{R}}^{d},

and also [57, Prop. 2.47] for the boundary transformation, we obtain

	$\displaystyle{\bm{F}}_{1}(t,{\bf u})[{\bf v}]$	$\displaystyle=\int_{T_{t}(\Omega)}\overline{{\mathbb{C}}}(\varphi\circ T_{t}^{-1})\mathcal{E}_{y}({\bf u}\circ T_{t}^{-1}):\mathcal{E}_{y}({\bf v}\circ T_{t}^{-1})\,\mathrm{d}y$
		$\displaystyle\quad-\int_{T_{t}(\Omega)}\overline{{\bf F}}\cdot({\bf v}\circ T_{t}^{-1})\,\mathrm{d}y-\int_{T_{t}(\overline{\Gamma}_{N})}\overline{{\bf g}}\cdot({\bf v}\circ T_{t}^{-1})\,\mathrm{d}\mathcal{H}^{d-1}_{y},$

where $\mathcal{E}_{y}({\bf u})=\frac{1}{2}(\nabla_{y}{\bf u}+(\nabla_{y}{\bf u})^{\top})$ for ${\bf u}:T_{t}(\Omega)\to{\mathbb{R}}^{d}$ and $\,\mathrm{d}\mathcal{H}^{d-1}_{y}$ denotes the $(d-1)$ -dimensional Hausdorff measure related to $\bm{y}$ . From the properties of the mapping $T_{t}$ , we find that $T_{t}(\overline{\Gamma}_{N})=\overline{\Gamma}_{N}$ and $\widetilde{{\bf v}}:={\bf v}\circ T_{t}^{-1}\in\overline{H}^{1}_{D}(T_{t}(\Omega),{\mathbb{R}}^{d})$ . Hence, from the above identity we observe that

{\bm{F}}_{1}(t,{\overline{{\bf u}}}^{t}\circ T_{t})[{\bf v}]=\int_{\Omega}\overline{{\mathbb{C}}}(\varphi^{t})\mathcal{E}_{y}({\overline{{\bf u}}}^{t}):\mathcal{E}_{y}(\widetilde{{\bf v}})-\overline{{\bf F}}\cdot\widetilde{{\bf v}}\,\mathrm{d}y-\int_{\overline{\Gamma}_{N}}\overline{{\bf g}}\cdot\widetilde{{\bf v}}\,\mathrm{d}\mathcal{H}^{d-1}_{y}=0

for all ${\bf v}\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ . Denoting by $D_{\bf u}{\bm{F}}_{1}$ as the partial derivative of ${\bm{F}}_{1}$ with respect to its second argument, we find that $D_{\bf u}{\bm{F}}_{1}(0,{\bf u}):\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})\to(\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}))^{*}$ is given by

\Big{(}D_{\bf u}{\bm{F}}_{1}(0,{\bf u})[{\bf v}]\Big{)}[{\bf w}]=\int_{\Omega}\overline{{\mathbb{C}}}(\varphi)\mathcal{E}({\bf v}):\mathcal{E}({\bf w})\,\mathrm{dx}

for all ${\bf v},{\bf w}\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ , where we have used the relation $\nabla T_{t}^{-1}|_{t=0}=\mathbb{I}$ the identity matrix. As $D_{\bf u}{\bm{F}}_{1}(0,{\bf u})$ is an isomorphism by the Lax–Milgram theorem, the application of the implicit function theorem allows us to infer that the mapping

t\mapsto({\overline{{\bf u}}}^{t}\circ T_{t})

is differentiable at $t=0$ with derivative $\dot{{\overline{{\bf u}}}}[{\bf V}]:=\partial_{t}|_{t=0}({\overline{{\bf u}}}^{t}\circ T_{t})\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ being the unique solution to the distributional equation

D_{\bf u}{\bm{F}}_{1}(0,{\overline{{\bf u}}})\big{[}\dot{{\overline{{\bf u}}}}[{\bf V}]\big{]}=-\partial_{t}{\bm{F}}_{1}(0,{\overline{{\bf u}}})\text{ in }(\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}))^{*},

which reads as

		$\displaystyle\int_{\Omega}\overline{{\mathbb{C}}}(\varphi)\mathcal{E}(\dot{{\overline{{\bf u}}}}[{\bf V}]):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}=\int_{\Omega}\overline{{\mathbb{C}}}(\varphi)(\nabla{\bf V}(0)\nabla{\overline{{\bf u}}})^{\mathrm{sym}}:\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}$		(5.11)
		$\displaystyle\quad+\int_{\Omega}\overline{{\mathbb{C}}}(\varphi)\mathcal{E}({\overline{{\bf u}}}):(\nabla{\bf V}(0)\nabla{\bm{\zeta}})^{\mathrm{sym}}-\overline{{\mathbb{C}}}(\varphi)\mathcal{E}({\overline{{\bf u}}}):\mathcal{E}({\bm{\zeta}})\mathop{\rm div}\nolimits{\bf V}(0)\,\mathrm{dx}$
		$\displaystyle\quad+\int_{\Omega}\Big{(}\nabla\overline{{\bf F}}{\bf V}(0)+\overline{{\bf F}}\mathop{\rm div}\nolimits({\bf V}(0))\Big{)}\cdot{\bm{\zeta}}\,\mathrm{dx}$
		$\displaystyle\quad+\int_{\overline{\Gamma}_{N}}\Big{(}\nabla\overline{{\bf g}}{\bf V}(0)+\overline{{\bf g}}\Big{(}\mathop{\rm div}\nolimits({\bf V}(0))-{\bf n}\cdot\nabla{\bf V}(0){\bf n}\Big{)}\Big{)}\cdot{\bm{\zeta}}\,\mathrm{d}\mathcal{H}^{d-1}$

for all ${\bm{\zeta}}\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ . In the above, we have made use of the following relations (see [57, Lem. 2.31, Prop. 2.36, Lem. 2.49])

	$\displaystyle\partial_{t}\nabla T_{t}\|_{t=0}=\nabla{\bf V}(0),\quad\partial_{t}\nabla T_{t}^{-1}\|_{t=0}=-\nabla{\bf V}(0),$
	$\displaystyle\partial_{t}\det\nabla T_{t}\|_{t=0}=\mathop{\rm div}\nolimits{\bf V}(0),\quad\partial_{t}(f\circ T_{t})\|_{t=0}=\nabla f\cdot{\bf V}(0),$
	$\displaystyle\partial_{t}(\det(\nabla T_{t})\\|(\nabla T_{t})^{-1}\cdot{\bf n}\\|)\|_{t=0}=\mathop{\rm div}\nolimits{\bf V}(0)-{\bf n}\cdot\nabla{\bf V}(0){\bf n}.$

Furthermore, substituting ${\bm{\zeta}}=\dot{{\overline{{\bf u}}}}[{\bf V}]$ into (5.11), by means of Korn’s inequality and the smoothness of ${\bf V}(0)$ , we obtain the estimate

\displaystyle\mathopen{\|}\dot{{\overline{{\bf u}}}}[{\bf V}]\mathclose{\|}_{H^{1}(\Omega)}\leq C\Big{(}\|{\overline{{\bf u}}}\|_{H^{1}(\Omega)}+\|\overline{{\bf F}}\|_{H^{1}(\Omega)}+\|\overline{{\bf g}}\|_{H^{2}(\overline{\Gamma}_{N})}\Big{)}.

(5.12)

Via a similar procedure, for a small $\tau_{0}>0$ , we consider ${\bm{F}}_{2}:(-\tau_{0},\tau_{0})\times\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})\to(\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}))^{*}$ defined as

	$\displaystyle{\bm{F}}_{2}(t,{\bf u})[{\bf v}]$	$\displaystyle=\int_{\Omega}\big{(}\widehat{{\mathbb{C}}}(\varphi)(\nabla T_{t}^{-1}\nabla{\bf u})^{\rm sym}:(\nabla T_{t}^{-1}\nabla{\bf v})^{\rm sym}-((\widehat{{\bf F}}\circ T_{t})\cdot{\bf v})\big{)}\det(\nabla T_{t})\,\mathrm{dx}$
		$\displaystyle\quad-\int_{\widehat{\Gamma}_{N}}((\widehat{{\bf g}}\circ T_{t})\cdot{\bf v}\det(\nabla T_{t})\\|\nabla T_{t}^{-\top}{\bf n}\\|\,\mathrm{d}\mathcal{H}^{d-1}$
		$\displaystyle\quad-\int_{\Omega}\chi(\varphi)\widehat{{\mathbb{C}}}(\varphi)(\nabla T_{t}^{-1}\nabla({\overline{{\bf u}}}^{t}\circ T_{t}))^{\rm sym}:(\nabla T_{t}^{-1}\nabla{\bf v})^{\rm sym}\det(\nabla T_{t})\,\mathrm{dx}.$

Then, by a change of variables $\bm{y}=T_{t}(\bm{x})$ , we find that

	$\displaystyle{\bm{F}}_{2}(t,{\widehat{{\bf u}}}^{t}\circ T_{t})[{\bf v}]$	$\displaystyle=\int_{\Omega}\widehat{{\mathbb{C}}}(\varphi^{t})\mathcal{E}_{y}({\widehat{{\bf u}}}^{t}):\mathcal{E}_{y}(\widetilde{{\bf v}})-\widehat{{\bf F}}\cdot\widetilde{{\bf v}}\,\mathrm{d}y-\int_{\widehat{\Gamma}_{N}}\widehat{{\bf g}}\cdot\widetilde{{\bf v}}\,\mathrm{d}\mathcal{H}^{d-1}_{y}$
		$\displaystyle\quad-\int_{\Omega}\chi(\varphi^{t})\widehat{{\mathbb{C}}}(\varphi^{t})\mathcal{E}_{y}({\overline{{\bf u}}}^{t}):\mathcal{E}_{y}(\widetilde{{\bf v}})\,\mathrm{d}y=0$

where $\widetilde{{\bf v}}={\bf v}\circ T_{t}^{-1}\in\widehat{H}^{1}_{D}(T_{t}(\Omega),{\mathbb{R}}^{d})$ . Since $D_{\bf u}{\bm{F}}_{2}(0,{\bf u}):\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})\to(\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}))^{*}$ given by

\Big{(}D_{\bf u}{\bm{F}}_{2}(0,{\bf u})[{\bf v}]\Big{)}[{\bf w}]=\int_{\Omega}\widehat{{\mathbb{C}}}(\varphi)\mathcal{E}({\bf v}):\mathcal{E}({\bf w})\,\mathrm{dx}\quad\text{ for all }{\bf v},{\bf w}\in\widehat{H}_{D}^{1}(\Omega,{\mathbb{R}}^{d})

is an isomorphism, by the implicit function theorem we infer that the mapping

t\mapsto({\widehat{{\bf u}}}^{t}\circ T_{t})

is differentiable at $t=0$ with derivative $\dot{{\widehat{{\bf u}}}}[{\bf V}]:=\partial_{t}|_{t=0}({\widehat{{\bf u}}}^{t}\circ T_{t})\in\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ being the unique solution to the distributional equation

D_{\bf u}{\bm{F}}_{2}(0,{\widehat{{\bf u}}})\big{[}\dot{{\widehat{{\bf u}}}}[{\bf V}]\big{]}=-\partial_{t}{\bm{F}}_{2}(0,{\widehat{{\bf u}}})\text{ in }(\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}))^{*},

which reads as

		$\displaystyle\int_{\Omega}\widehat{{\mathbb{C}}}(\varphi)\mathcal{E}(\dot{{\widehat{{\bf u}}}}[{\bf V}]):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}=\int_{\Omega}\widehat{{\mathbb{C}}}(\varphi)(\nabla{\bf V}(0)\nabla{\widehat{{\bf u}}}-\chi(\varphi)\nabla{\bf V}(0)\nabla{\overline{{\bf u}}})^{\mathrm{sym}}:\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}$		(5.13)
		$\displaystyle\quad+\int_{\Omega}\widehat{{\mathbb{C}}}(\varphi)(\mathcal{E}({\widehat{{\bf u}}})-\chi(\varphi)\mathcal{E}({\overline{{\bf u}}})):(\nabla{\bf V}(0)\nabla{\bm{\zeta}})^{\mathrm{sym}}+\chi(\varphi)\widehat{{\mathbb{C}}}(\varphi)\mathcal{E}(\dot{{\overline{{\bf u}}}}[{\bf V}]):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}$
		$\displaystyle\quad-\int_{\Omega}\widehat{{\mathbb{C}}}(\varphi)(\mathcal{E}({\widehat{{\bf u}}})-\chi(\varphi)\mathcal{E}({\overline{{\bf u}}})):\mathcal{E}({\bm{\zeta}})\mathop{\rm div}\nolimits{\bf V}(0)\,\mathrm{dx}$
		$\displaystyle\quad+\int_{\Omega}\Big{(}\nabla\widehat{{\bf F}}{\bf V}(0)+\widehat{{\bf F}}\mathop{\rm div}\nolimits({\bf V}(0))\Big{)}\cdot{\bm{\zeta}}\,\mathrm{dx}$
		$\displaystyle\quad+\int_{\widehat{\Gamma}_{N}}\Big{(}\nabla\widehat{{\bf g}}{\bf V}(0)+\widehat{{\bf g}}\Big{(}\mathop{\rm div}\nolimits({\bf V}(0))-{\bf n}\cdot\nabla{\bf V}(0){\bf n}\Big{)}\Big{)}\cdot{\bm{\zeta}}\,\mathrm{d}\mathcal{H}^{d-1}$

for all ${\bm{\zeta}}\in\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ . Substituting ${\bm{\zeta}}=\dot{{\widehat{{\bf u}}}}[{\bf V}]$ into (5.13), then using (5.12) and Korn’s inequality, we obtain the estimate

	$\displaystyle\mathopen{\\|}\dot{{\widehat{{\bf u}}}}[{\bf V}]\mathclose{\\|}_{H^{1}(\Omega)}$	$\displaystyle\leq C\Big{(}\\|{\widehat{{\bf u}}}\\|_{H^{1}(\Omega)}+\\|{\overline{{\bf u}}}\\|_{H^{1}(\Omega)}\Big{)}$		(5.14)
		$\displaystyle\quad+C\Big{(}\\|\overline{{\bf F}}\\|_{H^{1}(\Omega)}+\\|\overline{{\bf g}}\\|_{H^{2}(\overline{\Gamma}_{N})}+\\|\widehat{{\bf F}}\\|_{H^{1}(\Omega)}+\\|\widehat{{\bf g}}\\|_{H^{2}(\widehat{\Gamma}_{N})}\Big{)}.$		(5.14)

The next result details an optimality condition for minimisers $\varphi^{\varepsilon}$ to $(\bm{P}^{\varepsilon})$ obtained via geometric variations.

Theorem 5.1.

Assume (A1)–(A6), and additionally suppose that $\overline{{\bf F}},\widehat{{\bf F}}\in H^{1}(\Omega,{\mathbb{R}}^{d})$ , $\overline{\bm{U}}=\widehat{\bm{U}}=\bm{0},\overline{{\bf g}}\in H^{2}(\overline{\Gamma}_{N},{\mathbb{R}}^{d})$ , $\widehat{{\bf g}}\in H^{2}(\widehat{\Gamma}_{N},{\mathbb{R}}^{d})$ and ${\bf u}^{\rm tar}\in H^{2}(\Gamma^{\rm tar},{\mathbb{R}}^{d})$ hold. Let $\varphi^{\varepsilon}\in\mathcal{U}_{\rm ad}$ be a minimiser to $(\bm{P}^{\varepsilon})$ , with corresponding solutions ${\overline{{\bf u}}}^{\varepsilon}=\mathcal{S}_{1}^{2}(\varphi^{\varepsilon})$ and ${\widehat{{\bf u}}}^{\varepsilon}=\mathcal{S}(\varphi^{\varepsilon})$ . For every admissible velocity ${\bf V}\in\mathcal{V}_{\mathrm{ad}}$ , let $(\dot{{\overline{{\bf u}}}^{\varepsilon}}[{\bf V}],\dot{{\widehat{{\bf u}}}^{\varepsilon}}[{\bf V}])\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})\times\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ denote the unique solutions to (5.11) and (5.13) with $(\varphi,{\overline{{\bf u}}},{\widehat{{\bf u}}})=(\varphi^{\varepsilon},{\overline{{\bf u}}}^{\varepsilon},{\widehat{{\bf u}}}^{\varepsilon})$ . Then, for all ${\bf V}\in\mathcal{V}_{\mathrm{ad}}$ , it holds that

$\displaystyle 0$	$\displaystyle=\frac{1}{2}\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}^{\varepsilon}-{\bf u}^{\rm tar})\cdot({\widehat{{\bf u}}}^{\varepsilon}-{\bf u}^{\rm tar})\Big{(}\mathop{\rm div}\nolimits{\bf V}(0)-{\bf n}\cdot\nabla{\bf V}(0){\bf n}\Big{)}\,\mathrm{d}\mathcal{H}^{d-1}$	(5.15)
	$\displaystyle\quad+\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}^{\varepsilon}-{\bf u}^{\rm tar})\cdot(\dot{{\widehat{{\bf u}}}^{\varepsilon}}[{\bf V}]-(\nabla{\bf u}^{\rm tar}){\bf V}(0))\,\mathrm{d}\mathcal{H}^{d-1}$
	$\displaystyle\quad+\gamma\int_{\Omega}\Big{(}\frac{\varepsilon}{2}\|\nabla\varphi^{\varepsilon}\|^{2}+\frac{1}{\varepsilon}\Psi(\varphi^{\varepsilon})\Big{)}\mathop{\rm div}\nolimits{\bf V}(0)-\varepsilon\nabla\varphi^{\varepsilon}\cdot\nabla{\bf V}(0)\nabla\varphi^{\varepsilon}\,\mathrm{dx}.$

Remark 5.2.

With more regularity, it is possible to relate (5.15) to the optimality condition (4.21), see [14, 41] for more details.

Proof.

For any ${\bf V}\in\mathcal{V}_{\mathrm{ad}}$ , let $T\in\mathcal{T}_{\mathrm{ad}}$ be the associated transformation and consider the scalar function $g(t):=J_{\rm red}^{\varepsilon}(\varphi^{\varepsilon}\circ T_{t}^{-1})$ for $t\in(-\tau_{0},\tau_{0})$ , where $\tau_{0}$ is sufficiently small. As $\varphi^{\varepsilon}$ is a minimiser of $J_{\rm red}^{\varepsilon}$ , we have

g^{\prime}(0)=\frac{d}{dt}J_{\rm red}^{\varepsilon}(\varphi^{\varepsilon}\circ T_{t}^{-1})|_{t=0}=0.

The directional derivative

\frac{d}{dt}E_{\varepsilon}(\varphi^{\varepsilon}\circ T_{t}^{-1})|_{t=0}=\int_{\Omega}\Big{(}\frac{\varepsilon}{2}|\nabla\varphi^{\varepsilon}|^{2}+\frac{1}{\varepsilon}\Psi(\varphi^{\varepsilon})\Big{)}\mathop{\rm div}\nolimits{\bf V}(0)-\varepsilon\nabla\varphi^{\varepsilon}\cdot\nabla{\bf V}(0)\nabla\varphi^{\varepsilon}\,\mathrm{dx}

can be obtained as in [41, Lem. 7.5]. Denoting by ${\widehat{{\bf u}}}^{\varepsilon}(t)=\mathcal{S}(\varphi^{\varepsilon}\circ T_{t}^{-1})$ , then the derivative of $G(\varphi^{\varepsilon}\circ T_{t}^{-1})$ can be obtained by a standard change of variables:

	$\displaystyle\frac{d}{dt}G(\varphi^{\varepsilon}\circ T_{t}^{-1})$
	$\displaystyle\quad=\frac{1}{2}\frac{d}{dt}\int_{\Gamma^{\rm tar}}W(({\widehat{{\bf u}}}^{\varepsilon}(t)-{\bf u}^{\rm tar})\circ T_{t})\cdot(({\widehat{{\bf u}}}^{\varepsilon}(t)-{\bf u}^{\rm tar})\circ T_{t})\det(\nabla T_{t})\\|\nabla T_{t}^{-\top}{\bf n}\\|\,\mathrm{d}\mathcal{H}^{d-1}$
	$\displaystyle\quad=\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}^{\varepsilon}-{\bf u}^{\rm tar})\cdot\frac{d}{dt}({\widehat{{\bf u}}}^{\varepsilon}(t)\circ T_{t}-{\bf u}^{\rm tar}\circ T_{t})\|_{t=0}\,\mathrm{d}\mathcal{H}^{d-1}$
	$\displaystyle\qquad+\frac{1}{2}\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}^{\varepsilon}-{\bf u}^{\rm tar})\cdot({\widehat{{\bf u}}}^{\varepsilon}-{\bf u}^{\rm tar})\frac{d}{dt}(\det(\nabla T_{t})\\|\nabla T_{t}^{-\top}{\bf n}\\|)\,\mathrm{d}\mathcal{H}^{d-1}$
	$\displaystyle\quad=\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}^{\varepsilon}-{\bf u}^{\rm tar})\cdot(\dot{{\widehat{{\bf u}}}^{\varepsilon}}[{\bf V}]-(\nabla{\bf u}^{\rm tar}){\bf V}(0))\,\mathrm{d}\mathcal{H}^{d-1}$
	$\displaystyle\qquad+\frac{1}{2}\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}^{\varepsilon}-{\bf u}^{\rm tar})\cdot({\widehat{{\bf u}}}^{\varepsilon}-{\bf u}^{\rm tar})\Big{(}\mathop{\rm div}\nolimits{\bf V}(0)-{\bf n}\cdot\nabla{\bf V}(0){\bf n}\Big{)}\,\mathrm{d}\mathcal{H}^{d-1}$

leading to (5.15). ∎

The convergence of (5.15) to the sharp interface limit $\varepsilon\to 0$ is formulated as follows.

Theorem 5.2.

Suppose the hypotheses of Theorem 5.1 hold, and let $\varphi^{\varepsilon}\in\mathcal{U}_{\rm ad}$ be a minimiser to $(\bf P^{\varepsilon})$ . For any ${\bf V}\in\mathcal{V}_{\mathrm{ad}}$ with corresponding transformation $T\in\mathcal{T}_{\mathrm{ad}}$ , there exists a non-relabelled subsequence $\varepsilon\to 0$ such that

	$\displaystyle\varphi^{\varepsilon}$	$\displaystyle\to\varphi_{0}\text{ in }L^{1}(\Omega),\quad J_{\rm red}^{\varepsilon}(\varphi^{\varepsilon})\to J_{\rm red}^{0}(\varphi_{0})\text{ in }{\mathbb{R}},$
	$\displaystyle\dot{{\overline{{\bf u}}}^{\varepsilon}}[{\bf V}]$	$\displaystyle\rightharpoonup\dot{{\overline{{\bf u}}}}_{0}[{\bf V}]\text{ in }\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}),\quad\dot{{\widehat{{\bf u}}}^{\varepsilon}}[{\bf V}]\rightharpoonup\dot{{\widehat{{\bf u}}}}_{0}[{\bf V}]\text{ in }\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d}),$

where $\varphi_{0}\in\mathrm{BV}(\Omega,\{-1,1\})$ is a minimiser to the reduced functional $J_{\rm red}^{0}(\varphi)=G(\varphi)+\gamma TV(\varphi)$ , where the total variation $TV(\varphi)$ for $\varphi\in\mathrm{BV}(\Omega)$ is defined as

\displaystyle TV(\varphi)=\sup\Big{\{}\int_{\Omega}\varphi\mathop{\rm div}\nolimits\bm{\phi}\,\mathrm{dx}\text{ s.t. }\bm{\phi}\in C^{1}_{0}(\Omega,{\mathbb{R}}^{d}),\,\|\bm{\phi}\|_{L^{\infty}(\Omega)}\leq 1\Big{\}}.

Furthermore, $\dot{{\overline{{\bf u}}}}_{0}[{\bf V}]$ and $\dot{{\widehat{{\bf u}}}}_{0}[{\bf V}]$ satisfy (5.11) and (5.13), respectively, with $(\varphi,{\overline{{\bf u}}},{\widehat{{\bf u}}})$ replaced by $(\varphi_{0},{\overline{{\bf u}}}_{0}=\mathcal{S}_{1}^{2}(\varphi_{0}),{\widehat{{\bf u}}}_{0}=\mathcal{S}_{2}(\varphi_{0}))$ . Lastly, it holds that

\displaystyle\frac{d}{dt}J_{\rm red}^{\varepsilon}(\varphi^{\varepsilon}\circ T_{t}^{-1})|_{t=0}\to\frac{d}{dt}J_{\rm red}^{0}(\varphi_{0}\circ T_{t}^{-1})|_{t=0}\text{ in }{\mathbb{R}},

(5.16)

where

		$\displaystyle\frac{d}{dt}J_{\rm red}^{0}(\varphi_{0}\circ T_{t}^{-1})\|_{t=0}$		(5.17)
		$\displaystyle\quad=\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}_{0}-{\bf u}^{\rm tar})\cdot(\dot{{\widehat{{\bf u}}}}_{0}[{\bf V}]-(\nabla{\bf u}^{\rm tar}){\bf V}(0))\,\mathrm{d}\mathcal{H}^{d-1}$
		$\displaystyle\qquad+\frac{1}{2}\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}_{0}-{\bf u}^{\rm tar})\cdot({\widehat{{\bf u}}}_{0}-{\bf u}^{\rm tar})\Big{(}\mathop{\rm div}\nolimits{\bf V}(0)-{\bf n}\cdot\nabla{\bf V}(0){\bf n}\Big{)}\,\mathrm{d}\mathcal{H}^{d-1}$
		$\displaystyle\qquad+\gamma\int_{\Omega}\Big{(}\mathop{\rm div}\nolimits{\bf V}(0)-\mu\cdot\nabla{\bf V}(0)\mu\Big{)}d\|\mathrm{D}\mathrm{\mathcal{X}}_{\{\varphi_{0}=1\}}\|,$

with $\mu=\frac{\mathrm{D}\mathrm{\mathcal{X}}_{\{\varphi_{0}=1\}}}{|\mathrm{D}\mathrm{\mathcal{X}}_{\{\varphi_{0}=1\}}|}$ as the generalised unit normal on the set $\{\varphi_{0}=1\}$ .

Remark 5.3.

With more regularity, it is possible to relate (5.17) with the solvability condition (5.9) in the two-phase setting, see [14, 41] for more details.

Proof.

The first two assertions on the convergence of $\varphi^{\varepsilon}$ and $J_{\rm red}^{\varepsilon}(\varphi^{\varepsilon})$ come from Lemma 5.1. Consequently, by the calculations in the proof of [27, Thm. 4.2] we infer the convergence, as $\varepsilon\to 0$ ,

	$\displaystyle\int_{\Omega}\Big{(}\frac{\varepsilon}{2}\|\nabla\varphi^{\varepsilon}\|^{2}+\frac{1}{\varepsilon}\Psi(\varphi^{\varepsilon})\Big{)}\mathop{\rm div}\nolimits{\bf V}(0)-\varepsilon\nabla\varphi^{\varepsilon}\cdot\nabla{\bf V}(0)\nabla\varphi^{\varepsilon}\,\mathrm{dx}$
	$\displaystyle\quad\to\int_{\Omega}\Big{(}\mathop{\rm div}\nolimits{\bf V}(0)-\mu\cdot\nabla{\bf V}(0)\mu\Big{)}d\|\mathrm{D}\mathrm{\mathcal{X}}_{\{\varphi_{0}=1\}}\|.$

Next, from (5.11) and (5.13), we see that $\dot{{\overline{{\bf u}}}^{\varepsilon}}[{\bf V}]$ and $\dot{{\widehat{{\bf u}}}^{\varepsilon}}[{\bf V}]$ satisfy

\displaystyle\int_{\Omega}\overline{{\mathbb{C}}}(\varphi^{\varepsilon})\mathcal{E}(\dot{{\overline{{\bf u}}}^{\varepsilon}}[{\bf V}]):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}=\mathcal{F}_{\varphi^{\varepsilon},{\overline{{\bf u}}}^{\varepsilon}}({\bm{\zeta}}),\quad\int_{\Omega}\widehat{{\mathbb{C}}}(\varphi^{\varepsilon})\mathcal{E}(\dot{{\widehat{{\bf u}}}^{\varepsilon}}[{\bf V}]):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}=\mathcal{F}_{\varphi^{\varepsilon},{\overline{{\bf u}}}^{\varepsilon},{\widehat{{\bf u}}}^{\varepsilon}}({\bm{\zeta}}),

where $\mathcal{F}_{\varphi^{\varepsilon},{\overline{{\bf u}}}^{\varepsilon}}({\bm{\zeta}})$ and $\mathcal{F}_{\varphi^{\varepsilon},{\overline{{\bf u}}}^{\varepsilon},{\widehat{{\bf u}}}^{\varepsilon}}({\bm{\zeta}})$ denotes the right-hand sides of (5.11) and (5.13), respectively. Thanks to Corollary 3.1, along a non-relabelled subsequence, ${\overline{{\bf u}}}^{\varepsilon}\to{\overline{{\bf u}}}_{0}\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ and ${\widehat{{\bf u}}}^{\varepsilon}\to{\widehat{{\bf u}}}_{0}\in\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ strongly as $\varepsilon\to 0$ . Hence, together with the dominated convergence theorem, it is clear that, as $\varepsilon\to 0$ ,

\mathcal{F}_{\varphi^{\varepsilon},{\overline{{\bf u}}}^{\varepsilon}}({\bm{\zeta}})\to\mathcal{F}_{\varphi_{0},{\overline{{\bf u}}}_{0}}({\bm{\zeta}}),\quad\mathcal{F}_{\varphi^{\varepsilon},{\overline{{\bf u}}}^{\varepsilon},{\widehat{{\bf u}}}^{\varepsilon}}({\bm{\zeta}})\to\mathcal{F}_{\varphi_{0},{\overline{{\bf u}}}_{0},{\widehat{{\bf u}}}_{0}}({\bm{\zeta}}).

On the other hand, we infer from (5.12) and (5.14) that $\dot{{\overline{{\bf u}}}^{\varepsilon}}[{\bf V}]$ and $\dot{{\widehat{{\bf u}}}^{\varepsilon}}[{\bf V}]$ are uniformly bounded in $\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ and $\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ , which then implies the existence of limit functions $\dot{{\overline{{\bf u}}}}_{0}[{\bf V}]\in\overline{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ and $\dot{{\widehat{{\bf u}}}}_{0}[{\bf V}]\in\widehat{H}^{1}_{D}(\Omega,{\mathbb{R}}^{d})$ , where

	$\displaystyle\int_{\Omega}\overline{{\mathbb{C}}}(\varphi^{\varepsilon})\mathcal{E}(\dot{{\overline{{\bf u}}}^{\varepsilon}}[{\bf V}]):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}$	$\displaystyle\to\int_{\Omega}\overline{{\mathbb{C}}}(\varphi_{0})\mathcal{E}(\dot{{\overline{{\bf u}}}}_{0}[{\bf V}]):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx},$
	$\displaystyle\int_{\Omega}\widehat{{\mathbb{C}}}(\varphi^{\varepsilon})\mathcal{E}(\dot{{\widehat{{\bf u}}}^{\varepsilon}}[{\bf V}]):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}$	$\displaystyle\to\int_{\Omega}\widehat{{\mathbb{C}}}(\varphi_{0})\mathcal{E}(\dot{{\widehat{{\bf u}}}}_{0}[{\bf V}]):\mathcal{E}({\bm{\zeta}})\,\mathrm{dx}.$

Lastly, using the compactness of the boundary-trace operator, we obtain, as $\varepsilon\to 0$ ,

	$\displaystyle\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}^{\varepsilon}-{\bf u}^{\rm tar})\cdot(\dot{{\widehat{{\bf u}}}^{\varepsilon}}[{\bf V}]-(\nabla{\bf u}^{\rm tar}){\bf V}(0))\,\mathrm{d}\mathcal{H}^{d-1}$
	$\displaystyle\qquad+\frac{1}{2}\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}^{\varepsilon}-{\bf u}^{\rm tar})\cdot({\widehat{{\bf u}}}^{\varepsilon}-{\bf u}^{\rm tar})\Big{(}\mathop{\rm div}\nolimits{\bf V}(0)-{\bf n}\cdot\nabla{\bf V}(0){\bf n}\Big{)}\,\mathrm{d}\mathcal{H}^{d-1}$
	$\displaystyle\quad\to\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}_{0}-{\bf u}^{\rm tar})\cdot(\dot{{\widehat{{\bf u}}}}_{0}[{\bf V}]-(\nabla{\bf u}^{\rm tar}){\bf V}(0))\,\mathrm{d}\mathcal{H}^{d-1}$
	$\displaystyle\qquad+\frac{1}{2}\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}_{0}-{\bf u}^{\rm tar})\cdot({\widehat{{\bf u}}}_{0}-{\bf u}^{\rm tar})\Big{(}\mathop{\rm div}\nolimits{\bf V}(0)-{\bf n}\cdot\nabla{\bf V}(0){\bf n}\Big{)}\,\mathrm{d}\mathcal{H}^{d-1}.$

This shows (5.16) and completes the proof. ∎

6 Numerical simulations

In this section we present the finite element discretisation and showcase several numerical simulations in two and three dimensions for the two-phase case. Namely, we have $L=2$ and we consider the optimal distribution of a single type of active material within a passive material.

6.1 Finite element discretisation

We assume that $\Omega$ is a polyhedral domain and let $\mathcal{T}_{h}$ be a regular triangulation of $\Omega$ into disjoint open simplices. Associated with $\mathcal{T}_{h}$ are the piecewise linear finite element spaces

\displaystyle S^{h}=\left\{\zeta\in C^{0}(\overline{\Omega}):\,\zeta_{|_{o}}\in P_{1}(o)\,\forall o\in\mathcal{T}_{h}\right\}\quad\text{and}\quad\SS^{h}=S^{h}\times\cdots\times S^{h}=[S^{h}]^{d},

where we denote by $P_{1}(o)$ the set of all affine linear functions on $o$ , cf. [22]. In addition we define

\mathcal{V}^{h}=\left\{\zeta\in S^{h}:|\zeta|\leq 1\text{ in }\overline{\Omega}\right\},

(6.1)

as well as

\overline{\SS}^{h}_{D}=\left\{\bm{\eta}\in\SS^{h}:\bm{\eta}={\bf 0}\text{ on }\overline{\Gamma}_{D}\right\},\quad\widehat{\SS}^{h}_{D}=\left\{\bm{\eta}\in\SS^{h}:\bm{\eta}={\bf 0}\text{ on }\widehat{\Gamma}_{D}\right\}.

We also let $(\cdot,\cdot)$ denote the $L^{2}$ –inner product on $\Omega$ , and let $(\cdot,\cdot)^{h}$ be the usual mass lumped $L^{2}$ –inner product on $\Omega$ associated with $\mathcal{T}_{h}$ . In addition, $\langle{\bf A},{\bf B}\rangle_{{\mathbb{C}}}=({\mathbb{C}}{\bf A},{\bf B})$ for any fourth order tensor ${\mathbb{C}}$ and any matrices ${\bf A}$ and ${\bf B}$ . Finally, $\tau$ denotes a chosen uniform time step size.

We now introduce finite element approximations of the state equations (3.9) and (3.10), adjoint systems (4.18) and (4.20), as well as the optimality conditions (4.21) with a pseudo-time evolution based on an $L^{2}$ -gradient flow approach. In particular, we consider the obstacle potential (5.10) with $\widetilde{\Psi}(s)=\frac{1}{2}(1-s^{2})$ , which leads to a variational inequality.

The fully discrete numerical scheme is formulated as follows: Given $\varphi_{h}^{n-1}\in\mathcal{V}^{h}$ , find $({\overline{{\bf u}}}_{h}^{n},{\widehat{{\bf u}}}_{h}^{n},\widehat{{\bf q}}_{h}^{n},\overline{\bm{p}}_{h}^{n},\varphi_{h}^{n})\in\overline{\SS}^{h}_{D}\times\widehat{\SS}^{h}_{D}\times\widehat{\SS}^{h}_{D}\times\overline{\SS}^{h}_{D}\times\mathcal{V}^{h}$ such that


	$\displaystyle\langle\mathcal{E}({\overline{{\bf u}}}_{h}^{n}),\mathcal{E}({\bm{\zeta}})\rangle_{\overline{{\mathbb{C}}}(\varphi_{h}^{n-1})}=\left(\overline{{\bf F}},{\bm{\zeta}}\right)^{h}+\int_{\overline{\Gamma}_{N}}\overline{{\bf g}}\cdot{\bm{\zeta}}\,\mathrm{d}\mathcal{H}^{d-1}\qquad\forall{\bm{\zeta}}\in\overline{\SS}^{h}_{D},$		(6.2a)
	$\displaystyle\langle\mathcal{E}({\widehat{{\bf u}}}_{h}^{n})-\chi(\varphi_{h}^{n-1})\mathcal{E}({\overline{{\bf u}}}_{h}^{n}),\mathcal{E}({\bm{\zeta}})\rangle_{\widehat{{\mathbb{C}}}(\varphi_{h}^{n-1})}=\big{(}\widehat{{\bf F}},{\bm{\zeta}}\big{)}^{h}+\int_{\widehat{\Gamma}_{N}}\widehat{{\bf g}}\cdot{\bm{\zeta}}\,\mathrm{d}\mathcal{H}^{d-1}\qquad\forall{\bm{\zeta}}\in\widehat{\SS}^{h}_{D},$		(6.2b)
	$\displaystyle\langle\mathcal{E}(\widehat{{\bf q}}_{h}^{n}),\mathcal{E}({\bm{\zeta}})\rangle_{\widehat{{\mathbb{C}}}(\varphi_{h}^{n-1})}=\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}_{h}^{n}-{\bf u}^{\rm tar})\cdot{\bm{\zeta}}\,\mathrm{d}\mathcal{H}^{d-1}\qquad\forall{\bm{\zeta}}\in\widehat{\SS}^{h}_{D},$		(6.2c)
	$\displaystyle\langle\mathcal{E}(\overline{\bm{p}}_{h}^{n}),\mathcal{E}({\bm{\zeta}})\rangle_{\overline{{\mathbb{C}}}(\varphi_{h}^{n-1})}=\langle\chi(\varphi_{h}^{n-1})\mathcal{E}(\widehat{{\bf q}}_{h}^{n}),\mathcal{E}({\bm{\zeta}})\rangle_{\widehat{{\mathbb{C}}}(\varphi_{h}^{n-1})}\qquad\forall{\bm{\zeta}}\in\overline{\SS}^{h}_{D},$		(6.2d)
	$\displaystyle\left(\tfrac{\varepsilon}{\tau}(\varphi_{h}^{n}-\varphi_{h}^{n-1})-\tfrac{\gamma}{\varepsilon}\varphi_{h}^{n},\zeta-\varphi_{h}^{n}\right)^{h}+\gamma\varepsilon(\nabla\varphi_{h}^{n},\nabla(\zeta-\varphi_{h}^{n}))$
	$\displaystyle\qquad\quad+\langle\chi^{\prime}(\varphi_{h}^{n-1})\mathcal{E}({\overline{{\bf u}}}_{h}^{n}),(\zeta-\varphi_{h}^{n})\mathcal{E}(\widehat{{\bf q}}_{h}^{n})\rangle_{\widehat{{\mathbb{C}}}(\varphi_{h}^{n-1})}$
	$\displaystyle\qquad\quad-\langle\mathcal{E}({\widehat{{\bf u}}}_{h}^{n})-\chi(\varphi_{h}^{n-1})\mathcal{E}({\overline{{\bf u}}}_{h}^{n}),(\zeta-\varphi_{h}^{n})\mathcal{E}(\widehat{{\bf q}}_{h}^{n})\rangle_{\widehat{{\mathbb{C}}}^{\prime}(\varphi_{h}^{n-1})}$
	$\displaystyle\quad\qquad-\langle\mathcal{E}({\overline{{\bf u}}}_{h}^{n}),(\zeta-\varphi_{h}^{n})\mathcal{E}(\overline{\bm{p}}_{h}^{n})\rangle_{\overline{{\mathbb{C}}}^{\prime}(\varphi_{h}^{n-1})}\geq 0\qquad\forall\zeta\in\mathcal{V}^{h}.$		(6.2e)

We implemented the scheme (6.2) with the help of the finite element toolbox ALBERTA, see [56]. To increase computational efficiency, we employ adaptive meshes, which have a finer mesh size $h_{f}$ within the diffuse interfacial regions and a coarser mesh size $h_{c}$ away from them, see [9, 10] for a more detailed description.

Clearly, we first solve the linear systems (6.2a) in order to obtain ${\overline{{\bf u}}}_{h}^{n}$ , then (6.2b) for ${\widehat{{\bf u}}}_{h}^{n}$ , then (6.2c) for $\widehat{{\bf q}}_{h}^{n}$ , then (6.2d) for $\overline{\bm{p}}_{h}^{n}$ , and finally the variational inequality (6.2e) for $\varphi_{h}^{n}$ . In two space dimensions we employ the package LDL, see [23], together with the sparse matrix ordering AMD, see [6], in order to solve the linear systems (6.2a)–(6.2d). In three space dimensions, on the other hand, we use a preconditioned conjugate gradient algorithm, with a $W$ -cycle multigrid step as preconditioner. In order to solve the variational inequality (6.2e) we employ a nonlinear multigrid method similar to [42].

For the computational domain we will choose $\overline{\Omega}=[0,L_{1}]\times[-\frac{1}{2},\frac{1}{2}]$ in two dimensions, and $\overline{\Omega}=[0,L_{1}]\times[-\frac{1}{2}L_{2},\frac{1}{2}L_{2}]\times[-\frac{1}{2},\frac{1}{2}]$ in three dimensions with positive lengths $L_{i}$ , $i\in\{1,2\}$ , given below. For the physical parameters we loosely follow the settings in [50]. In particular, for the forcings we choose $\overline{{\bf F}}=\widehat{{\bf F}}={\bf 0}$ throughout, as well as $\widehat{\bm{g}}={\bf 0}$ and

\overline{\bm{g}}(\bm{x})=\begin{cases}g{\bm{e}}_{1}&x_{1}=L_{1},\\ {\bf 0}&x_{1}<L_{1},\end{cases}\qquad\text{with}\quad g=0.1.

(6.3)

For the interpolated elasticity tensors we choose $\overline{{\mathbb{C}}}(s)=\frac{1}{2}[(1+s)\overline{{\mathbb{C}}}_{+}+(1-s)\overline{{\mathbb{C}}}_{-}]$ , where the two tensors $\overline{{\mathbb{C}}}_{\pm}$ are defined through the Young’s moduli $\overline{E}_{\pm}$ and Poisson ratios $\overline{\nu}_{\pm}$ via

\overline{E}_{+}=3,\ \overline{E}_{-}=0.7,\quad\overline{\nu}_{+}=\overline{\nu}_{-}=0.45.

(6.4)

Similarly, $\widehat{{\mathbb{C}}}(s)=\frac{1}{2}[(1+s)\widehat{{\mathbb{C}}}_{+}+(1-s)\widehat{{\mathbb{C}}}_{-}]$ , with the Young’s moduli and Poisson ratios

\widehat{E}_{+}=13,\ \widehat{E}_{-}=0.6,\quad\widehat{\nu}_{+}=\widehat{\nu}_{-}=0.45.

(6.5)

For the fixity function we choose

\chi(s)=\tfrac{2}{5}(1+s).

(6.6)

For the visualisation of the progress of the discrete gradient flow computations, we define the discrete cost functional

	$\displaystyle J^{h}(\varphi_{h}^{h},{\widehat{{\bf u}}}_{h}^{n})$	$\displaystyle=\gamma\left(\frac{\varepsilon}{2}\|\nabla\varphi_{h}^{n}\|^{2}+\frac{1}{\varepsilon}\psi(\varphi_{h}^{n}),1\right)+{\frac{1}{2}}\int_{\Gamma^{\rm tar}}W({\widehat{{\bf u}}}_{h}^{n}-\bm{u}^{\rm tar})\cdot({\widehat{{\bf u}}}_{h}^{n}-\bm{u}^{\rm tar})\,\mathrm{d}\mathcal{H}^{d-1}$
		$\displaystyle=:\gamma\mathcal{E}^{h}(\varphi_{h}^{n})+\mathcal{E}^{\rm h,tar}({\widehat{{\bf u}}}_{h}^{n}).$		(6.7)

As the initial data $\varphi_{h}^{0}$ we choose a random mixture with mean zero. Choosing other initial data, including random mixtures with positive or negative mean had no visible influence on the obtained optimal shapes.

6.2 Numerical simulations in two dimensions

For the target shapes we consider the parabolic profile

	$\bm{u}^{\rm tar}(x_{1},x_{2})=u^{\rm tar}(x_{1}){\bm{e}}_{2},\quad u^{\rm tar}(x_{1})=c^{\rm tar}(x_{1})^{2},$	(6.8a)
with $c^{\rm tar}>0$ , and the cosine profile
	$\bm{u}^{\rm tar}(x_{1},x_{2})=u^{\rm tar}(x_{1}){\bm{e}}_{2},\quad u^{\rm tar}(x_{1})=c^{\rm tar}\Big{(}1-\cos\frac{k^{\rm tar}\pi x_{1}}{L_{1}}\Big{)}$	(6.8b)

with $c^{\rm tar}>0$ and $k^{\rm tar}>0$ . In Figure 2 we plot two examples for the profiles in (6.8) for the domain length $L_{1}=12$ .

In each of the Figures 3, 4, 5 and 6, we provide visualisations of the numerical solution $\varphi_{h}^{n}$ at various pseudo-times (black denotes the passive material $\{\varphi_{h}^{n}=-1\}$ and grey denotes the active material $\{\varphi_{h}^{n}=1\}$ ), and the corresponding displacements ${\overline{{\bf u}}}_{h}^{n}$ (in red) and ${\widehat{{\bf u}}}_{h}^{n}$ (in green). As mentioned before, each time the gradient flow is started from a random mixture $\varphi_{h}^{0}$ with zero mean. We also provide pseudo-time plots of the cost functional $J^{h}(\varphi_{h}^{n},{\widehat{{\bf u}}}_{h}^{n})$ , the proportion of the elastic $\mathcal{E}^{\rm h,tar}({\widehat{{\bf u}}}_{h}^{n})$ and interfacial $\gamma\mathcal{E}^{h}(\varphi_{h}^{n})$ energies, as well as log-plots of the elastic energies. The parameter details are summarised in the Table 1.

Figure	Profile	Domain	$W$	$\Gamma^{\rm tar}$	$c^{\rm tar}$	$k^{\rm tar}$
3	Parabolic (6.8a)	$[0,6]\times[-\frac{1}{2},\frac{1}{2}]$	$\mathrm{Id}$	$\partial\Omega$	$0.075$	-
4	Cosine (6.8b)	$[0,6]\times[-\frac{1}{2},\frac{1}{2}]$	$\bm{e}_{2}\otimes\bm{e}_{2}$	$\partial\Omega$	$0.25$	$2$
5	Parabolic (6.8a)	$[0,12]\times[-\frac{1}{2},\frac{1}{2}]$	$\mathrm{Id}$	$\partial\Omega$	$0.02$	-
6	Cosine (6.8b)	$[0,12]\times[-\frac{1}{2},\frac{1}{2}]$	$\mathrm{Id}$	$\partial_{\rm top}\Omega$	$1$	$1.5$

Table 1: Parameter details for numerical simulations in 2D.

For all the presented simulations we choose the parameters $\varepsilon=\frac{1}{8\pi}$ and $\gamma=0.01$ . In each case the cost functional decreases monotonically, but the proportions of the two energies (elastic vs interfacial) differ from case to case.

The first experiment in Figure 3 is for the parabolic profile on the domain $[0,6]\times[-\frac{1}{2},\frac{1}{2}]$ . We observe that in the optimal distribution of material, the active phase occupies most of the lower domain. This ensures that in the programmed stage, the printed active composite is able to attain the desired target shape.

Not surprisingly, a very different distribution of material is obtained when changing the target shape functional to enforce a cosine profile at the programmed stage. As can be seen from Figure 3, the optimal distribution of the active material is now given by an elongated region that connects the lower left of the domain with the upper right.

It turns out that on longer (or thinner) domains, far less active material is needed to achieve significant deformations at the programmed state. For example, in Figure 4 a miniscule amount of active material, spread in several connected components arranged at the bottom of the domain, is sufficient to result in the desired parabolic target shape.

Similarly, we observe in Figure 5 that three strategically placed small amounts of active material guarantee a cosine profile at the programmed stage for the printed active composite.

6.3 Numerical simulations in three dimensions

For the target shapes we consider a function of the form

\bm{u}^{\rm tar}(x_{1},x_{2},x_{3})=u^{\rm tar}(x_{1},x_{2}){\bm{e}}_{3}

and choose $u^{\rm tar}$ from one of the following:

u^{\rm tar}(x_{1},x_{2})=c^{\rm tar}(x_{1})^{2}.

(6.9)

u^{\rm tar}(x_{1},x_{2})=c^{\rm tar}\Big{(}1-\cos\frac{k^{\rm tar}\pi x_{1}}{L_{1}}\Big{)}.

(6.10)

u^{\rm tar}(x_{1},x_{2})=c^{\rm tar}x_{1}x_{2}.

(6.11)

Notice that (6.9) and (6.10) are simply the three dimensional analogues of the parabolic profile (6.8a) and the cosine profile (6.8b), respectively. On the other hand, (6.11) yields a linear profile in $x_{1}$ with a twisting in the $x_{2}$ -direction.

In each of the Figures 7, 8, and 9, 10 we provide visualisations of the numerical solution $\varphi_{h}^{n}$ at various pseudo-times (black denotes the passive material $\{\varphi_{h}^{n}=-1\}$ and grey denotes the active material $\{\varphi_{h}^{n}=1\}$ ), and the corresponding displacements ${\widehat{{\bf u}}}_{h}^{n}$ (darker colours indicating lower values of ${\widehat{{\bf u}}}_{h}^{n}\cdot\bm{e}_{3}$ and lighter colours indicating higher values of ${\widehat{{\bf u}}}_{h}^{n}\cdot\bm{e}_{3}$ ). We also provide pseudo-time plots of the energy functionals, similarly to the 2D simulations in the previous subsection. The parameter details are summarised in the Table 2.

Figure	Profile	Domain	$W$	$\Gamma^{\rm tar}$	$c^{\rm tar}$	$k^{\rm tar}$
7	(6.9)	$[0,6]\times[-\frac{3}{2},\frac{3}{2}]\times[0,1]$	$\mathrm{Id}$	$\partial_{\rm top}\Omega$	$0.075$	-
8	(6.10)	$[0,12]\times[-\frac{3}{2},\frac{3}{2}]\times[0,1]$	$\bm{e}_{3}\otimes\bm{e}_{3}$	$\partial_{\rm top}\Omega$	1	2
9	(6.10)	$[0,12]\times[-\frac{3}{2},\frac{3}{2}]\times[0,1]$	$\mathrm{Id}$	$\partial_{\rm top}\Omega$	$1$	$2$
10	(6.11)	$[0,6]\times[-3,3]\times[0,1]$	$\bm{e}_{3}\otimes\bm{e}_{3}$	$\partial_{\rm right}\Omega$	$0.1$	-

Table 2: Parameter details for numerical simulations in 3D.

In the first three figures we take $\varepsilon=\frac{1}{4\pi}$ and for $\gamma$ choose either $0.1$ or $0.01$ . In each simulation the cost functional decreases monotonically, but the proportions of the two energies (elastic vs interfacial) differ from case to case.

The first simulation is a direct 3D analogue for the computation previously shown in Figure 3, the only difference being that here we restrict the set $\Gamma^{\rm tar}$ to the upper part of the boundary $\partial\Omega$ . As expected, the observed results are very close to the ones seen previously in the 2D setting. In particular, the active phase occupies the lower half of the domain, with a dip towards the right end of the domain.

On more elongated domains we once again observe that relatively little active material can result in large deformations at the programmed stage. For example, the strategic placement of the active component seen in Figure 8 yields a large cosine profile deformation.

It is interesting to note that by simply changing the weighting matrix $W$ in the target energy functional, we obtain a completely different optimal distribution of active material. This can be seen in Figure 9, where the only change to the previous simulation is $W=\mathrm{Id}$ , rather than $W=\bm{e}_{3}\otimes\bm{e}_{3}$ . Now there are just two connected components for the active phase, one at the lower left part of the domain, and one at the middle of the top of the domain. Yet the obtained deformations at the programmed stage are very similar.

Our final numerical simulation is for a twisted target shape. In particular, we use the target function (6.11) with $c^{\rm tar}=0.1$ on the domain $\overline{\Omega}=[0,6]\times[-3,3]\times[0,1]$ , with $W=\bm{e}_{3}\otimes\bm{e}_{3}$ and $\Gamma^{\rm tar}=\partial_{\rm right}\Omega$ . In Figure 10 we provide visualisations of the numerical solution $\varphi_{h}^{n}$ at various pseudo-times, and the corresponding displacements ${\widehat{{\bf u}}}_{h}^{n}$ (here darker colours indicate lower values of $|{\widehat{{\bf u}}}_{h}^{n}|$ and lighter colours indicate higher values of $|{\widehat{{\bf u}}}_{h}^{n}|$ ). For this computation we take $\varepsilon=\frac{1}{2\pi}$ and $\gamma=0.02$ . At the optimal configuration we observe an elaborate distribution of the active material, which yields a twisted shape of the component at the programmed stage.

Acknowledgments

The authors HG and AS gratefully acknowledge the support by the Graduiertenkolleg 2339 IntComSin of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 321821685. The work of KFL is supported by the Research Grants Council of the Hong Kong Special Administrative Region, China [Project No.: HKBU 14302218 and HKBU 12300321].

References

[1] S.A. Abdullah, A. Jumahat, N.R. Abdullah and L. Frormann. Setermination of Shape Fixity and Shape Recovery Rate of Carbon Nanotube-flled Shape Memory Polymer Nanocomposites. Procedia Eng. 41 (2012) 1641–1646
[2] M.H. Ali, A. Abilgaziyev and D. Adair. 4D printing: a critical review of current developments, and future prospects. Int. J. Adv. Manuf. Technol. 104 (2019) 701–717
[3] G. Allaire and B. Bogosel. Optimizing supports for additive manufacturing. Struct. Multidiscip. Optim. 58 (2018) 2493–2515
[4] S. Almi and U. Stefanelli. Topology optimization for incremental elastoplasticity: A phase-field approach. SIAM J. Control Optim. 59 (2021) 339–364
[5] H.W. Alt: Linear Functional Analysis, an Application Oriented Introduction. Springer, London, (2016)
[6] P.R. Amestoy, T.A. Davis and I.S. Duff. Algorithm 837: AMD, an approximate minimum degree ordering algorithm. ACM Trans. Math. Software 30 (2004) 381–388
[7] S.E. Bakarich, E. Gorkin and G.M. Spinks. 4D printing with mechanically robust, thermally actuating hydrogels. Macromol. Rapid Commun. 36 (2015) 1211–1217
[8] S. Baldo. Minimal interface criterion for phase transitions in mixtures of Cahn–Hilliard fluids. Ann. Inst. H. Poincaré Anal. Non Linéaire 7 (1990) 67–90
[9] L. Baňas and R. Nürnberg. Finite element approximation of a three dimensional phase field model for void electromigration. J. Sci. Comp. 37 (2008) 202–232
[10] J.W. Barrett, R. Nürnberg and V. Styles. Finite element approximation of a phase field model for void electromigration. SIAM J. Numer. Anal. 42 (2004) 738–772
[11] G. Bellettini, A. Braides and G. Riey. Variational approximation of anisotropic functionals on partitions. Annali di Matematica 184 (2005) 75–93
[12] A. Bhattacharyya and K.A. James. Topology optimization of shape memory polymers structures with programmable morphology. Struct. Multidiscip. Optim. 63 (2021) 1863–1887
[13] L. Blank, H. Garcke, M. Hassan Farshbaf-Shaker and V. Styles. Relating phase field and sharp interface approaches to structural topology optimization. ESAIM: Control Optim. Calc. Var. 20 (2014) 1025–1058
[14] L. Blank, H. Garcke, C. Hecht and C. Rupprecht. Sharp interface limit for a phase field model in structural optimization. SIAM J. Control Optim. 54 (2016) 1558–1584
[15] L. Blank, H. Garcke, L. Sarbu, T. Srisupattarawanit, V. Styles and A. Voigt. Phase-field approaches to structural topology optimization, in Constrained optimization and optimal control for partial differential equations, vol. 160 of Internat. Ser. Numer. Math., Birkhäuser/Springer Basel AG, Basel, 2012, 245–256
[16] J.F. Blowey and C.M. Elliott. Curvature dependent phase boundary motion and parabolic double obstacle problems, in Degenerate diffusions (Minneapolis, MN, 1991), vol. 47 of IMA Vol. Math. Appl., Springer, New York, 1993, 19–60
[17] B. Bourdin and A. Chambolle. Design-dependent loads in topology optimization. ESAIM: Control Optim. Cal. Var. 9 (2003) 19–48
[18] L. Bronsard, H. Garcke and B. Stoth. A multi-phase Mullins-Sekerka system: matched asymptotic expansions and an implicit time discretization for the geometric evolution problem. SIAM J. Appl. Math. 60 (1999) 295–315
[19] M. Burger and R. Stainko. Phase-field relaxation of topology optimization with local stress constraints. SIAM J. Control Optim. 45 (2006) 1447–1466
[20] S. Cacace, E. Cristiani and L. Rocchi. A Level Set Based Method for Fixing Overhangs in 3D Printing. Appl. Math. Model. 44 (2017) 446–455
[21] M. Carraturo, E. Rocca, E. Bonetti, D. Hömberg, A. Reali and F. Auricchio. Graded-material design based phase-field and topology optimization. Comput. Mech. 64 (2019) 1589–1600
[22] P.G. Ciarlet: The Finite Element Method for Elliptic Problems. North-Holland Publishing Co., Amsterdam, (1978)
[23] T.A. Davis, Algorithm 849: a concise sparse Cholesky factorization package. ACM Trans. Math. Software 31 (2005) 587–591
[24] Z. Ding, C. Yuan, X. Peng, T. Wang, JH.J. Qi and M.L. Dunn. Direct 4D printing via active composite materials. Sci. Adv. 3 (2017) e1602890
[25] P.C. Fife and O. Penrose. Interfacial dynamics for thermodynamically consistent phase-field models with nonconserved order parameter. J. Differential Equations 16 (1995) 1–49
[26] J. Jiang, X. Xu and J. Stringer. Support Structures for Additive Manufacturing: A Review. J. Mauf. Mater. Process. 2 (2008) 64 (23 pages)
[27] H. Garcke. The $\Gamma$ -limit of the Ginzburg–Landau energy in an elastic medium. AMSA 18 (2008) 345–379
[28] H. Garcke and C. Hecht. Apply a phase field approach for shape optimization of a stationary Navier-Stokes flow. ESAIM: Control Optim. Cal. Var. 22 (2016) 309–337
[29] H. Garcke and C. Hecht. Shape and topology optimization in Stokes flow with a phase field approach. Appl. Math. Optim. 73 (2016) 23–70
[30] H. Garcke, C. Hecht, M. Hinze, C. Kahle and K.F. Lam. Shape optimization for surface functionals in Navier–Stokes flow using a phase field approach. Interfaces Free Bound. 18 (2016) 219–261
[31] H. Garcke, P. Hüttl and P. Knopf. Shape and topology optimization involving the eigenvalues of an elastic structure: A multi-phase-field approach. Adv. Nonlinear Anal. 11 (2022) 159–197
[32] H. Garcke, K.F. Lam, R. Nürnberg and A. Signori. Overhang penalization in additive manufacturing via phase field structural topology optimization with anisotropic energies. Preprint ArXiv:2111.14070
[33] H. Garcke, B. Nestler and B. Stoth. On anisotropic order parameter models for multi0phase systems and their sharp interface limits. Phys. D 115 (1998) 87–108
[34] H. Garcke and B. Stinner. Second order phase field asymptotics for multi-component systems. Interfaces Free Bound. 8 (2006) 131–157
[35] Q. Ge, C.K. Dunn, H.J. Qi and M.L. Dunn. Active origami by 4D printing. Smart Mater. Struct. 23 (2014) 094007
[36] O. Kuksenok and A.C. Balazs. Stimuli-responsive behavior of composites integating thermo-responsive gels with photo-responsive fibers. Mater. Horizons 3 (2016) 53–62
[37] Q. Ge, H.J. Qi and M.L. Dunn. Active materials by four-dimension printing. Appl. Phys. Lett. 103 (2013) 131901
[38] Q. Ge, A.H. Sakhaei, H. Lee, C.K. Dunn, N.X. Fang and M.L. Dunn. Multimaterial 4D Printing with Tailorable Shape Memory Polymers. Sci. Rep. 6 (2016) 31110
[39] M.J. Geiss, N. Boddeti, O. Weeger, K. Maute and M.L. Dunn. Combined level-set-XFEM-density topology optimization of four-dimensional printed structures undergoing large deformation. J. Mech. Des. 141 (2019) 051405
[40] M. Howard, J. Pajot, K. Maute and M.L. Dunn. A Computational Design Methodology for Assembly and Actuation of Thin-Film Structures via Patterning of Eigenstrains. J. Microelectromech. Syst. 18 (2009) 1137–1148
[41] C. Hecht. Shape and topology optimization in fluids using a phase field approach an application in structural optimization. PhD Thesis University of Regensburg, Regensburg, Germany 2014 https://epub.uni-regensburg.de/29869/1/Dissertation_ClaudiaHecht.pdf
[42] R. Kornhuber. Monotone multigrid methods for elliptic variational inequalities II. Numer. Math. 72 (1996) 481–499
[43] M. Leray, L. Merli, F. Torti, M. Mazur and M. Bandt. Optimal topology for additive manufacturing: A method for enabling additive manufacture of support-free optimal structures. Mater. Des. 63 (2014) 678–690
[44] J. Liu, A.T. Gaynor, S. Chen, Z. Kang, K. Suresh, A. Takezawa, L. Li, J. Kato, J. Tang, C.C.L. Wang, L. Cheng, X. Liang and A.C. To. Current and future trends in topology optimization for additive manufacturing. Struct. Multi. Optim. 57 (2018) 2457–2483
[45] J. Liu and H. Yu. Self-Supporting Topology Optimization With Horizontal Overhangs for Additive Manufacturing. J. Manuf. Sci. Eng. 142 (2020) 091003 (14 pages)
[46] M. Langelaar. Combined optimization of part topology, support structure layout and build orientation for additive manufacturing. Struct. Multidiscip. Optim. 57 (2018) 1985–2004
[47] H. Li, X. Gao and Y. Luo. Multi-shape memory polymers achieved by the spatio-assembly of 3D printable thermoplastic building blocks. Soft Matter 12 (2016) 3226–3233
[48] J.L. Lions and E. Magenes. Non-homogeneous boundary value problems and applications. Springer, 1972.
[49] T.S. Lumpe and K. Shea. Computational design of 3D-printed active lattice structures for reversible shape morphing. J. Mater. Res. 36 (2021) 3642–3655
[50] K. Maute, A. Tkachuk, J. Wu, H.J. Qi, Z. Ding and M.L. Dunn. Level Set Topology Optimization of Printed Active Composites. J. Mech. Des. 137 (2015) 111402 (13 pages)
[51] A.M. Mirzendehdel and K. Suresh. Support Structure Constrained Topology Optimization for Additive Manufacturing. Comput. Aided Des. 1 (2016) 1–13
[52] F. Momeni, S.M. Hassani, N.X. Liu and J. Ni. A review of 4D printing. Mater. Des. 122 (2017) 42–79
[53] N.C. Owen, J. Rubinstein, P. Sternberg. Minimizers and gradient flows for singularly perturbed bi-stable potentials with a Dirichlet condition. Proc. R. Soc. Lond. A 429 (1990) 505–532
[54] J.M. Pajot, K. Maute, Y. Zhang and M.L. Dunn. Design of patterned multilayer films with eigenstrains by topology optimization. Int. J. Solids Struct. 43 (2006) 1832–1853
[55] P. Penzler, M. Rumpf and B. Wirth. A phase-field model for compliance shape optimization in nonlinear elasticity. ESAIM Control Optim. Cal. Var. 18 (2012) 229–258
[56] A. Schmidt and K.G. Siebert. Design of Adaptive Finite Element Software: The Finite Element Toolbox ALBERTA. vol. 42 of Lecture Notes in Computational Science and Engineering, Springer-Verlag, Berlin, 2005
[57] J. Sokolowski and J.-P. Zolesio. Introduction to Shape Optimization. Springer series in computational mathematics, Vol. 16, Springer-Verlag, Berlin Heidelberg, 1991
[58] L. Sun and W.M. Huang. Mechanisms of the multi-shape memory effect and temperature memory effect in shape memory polymers. Soft Matter 6 (2010) 4403–4406
[59] Y. Sun, W. Ouyang, Z. Liu, N. Ni, Y. Savoye, P. Song and L. Liu. Computational Design of Self-actuated Deformable Solids via Shape Memory Material. IEEE Trans. Vis. Comput. Graph. —bf 28 (2022) 2577–2588
[60] A. Takezawa, S. Nishiwaki and M. Kitamura. Shape and topology optimization based on the phase field method and sensitivity analysis. J. Comput. Phys. 229 (2010) 2697–2718
[61] S. Tibbits. 4D printing: Multi-material shape change. Archit. Des. 84 (2014) 116–121
[62] M.Y. Wang and S. Zhou. Phase field: a variational method for structural topology optimization. CMES Comput. Model. Eng. Sci. 6 (2004) 547–566
[63] X. Wan, Y. He, Y. Liu and J. Leng. 4D printing of multiple shape memory polymer and nanocomposites with biocompatible programmable and selectively actuated properties. Addit. Manuf. 53 (2022) 102689
[64] G. Wang, H. Yang, Z. Yan, N.G. Ulu, Y. Tao, J. Gu, L.B. Kara and L. Yao. 4DMesh: 4d printing morphing non-developable mesh surfaces. in the 31st Annual ACM Symposium on User Interface Software and Technology (2018) 623–635
[65] Y. Wei, P. Huang, Z. Li, P. Wang and X. Feng. Design of active materials distributions for four-dimensional printing based on multi-material topology optimization. Smart Mater. Struct. 30 (2021) 095002
[66] J. Wu, C. Yuan, Z. Ding, M. Isakov, Y. Mao, T. Wang, M.L. Dunn and H.J. Qi. Multi-shape active composites by 3D printing of digital shape memory polymers. Sci. Rep. 6 (2016) 24224
[67] C. Yuan, T. Lu and T.J. Wang. Mechanics-based design strategies for 4D printing: A review. Forces Mech. 7 (2022) 100081.
[68] Q. Zhang, K. Zhang and G. Hu. Smart three-dimensional lightweight structure triggered from a thin composite sheet via 3D printing technique. Sci. Rep. 6 (2016) 22431

$\displaystyle\|h_{\Omega}(\bm{x},\bm{v},\bm{w})\|$	$\displaystyle\leq a_{1}(\bm{x})+b_{1}(\bm{x})\|\bm{v}\|^{p}+b_{2}(\bm{x})\|\bm{w}\|^{p}\quad$	$\displaystyle\text{ for all }\bm{v},\bm{w}\in{\mathbb{R}}^{d},\text{ a.e.~{}}\bm{x}\in\Omega,$
$\displaystyle\|\overline{h}(\bm{s},\bm{v})\|$	$\displaystyle\leq a_{2}(\bm{s})+c_{1}(\bm{s})\|\bm{v}\|^{2}\quad$	$\displaystyle\text{ for all }\bm{v}\in{\mathbb{R}}^{d},\text{ a.e.~{}}\bm{s}\in\overline{\Gamma}_{N},$
$\displaystyle\|\overline{h}(\bm{s},\bm{w})\|$	$\displaystyle\leq a_{3}(\bm{s})+c_{2}(\bm{s})\|\bm{w}\|^{2}\quad$	$\displaystyle\text{ for all }\bm{w}\in{\mathbb{R}}^{d},\text{ a.e.~{}}\bm{s}\in\widehat{\Gamma}_{N},$

	$\displaystyle\mathopen{\\|}{\overline{{\bf u}}}_{1}-{\overline{{\bf u}}}_{2}\mathclose{\\|}_{H^{1}(\Omega)}$	$\displaystyle\leq C\mathopen{\\|}{\bm{\varphi}}_{1}-{\bm{\varphi}}_{2}\mathclose{\\|}_{L^{\infty}(\Omega)},$		(3.12)
	$\displaystyle\mathopen{\\|}{\widehat{{\bf u}}}_{1}-{\widehat{{\bf u}}}_{2}\mathclose{\\|}_{H^{1}(\Omega)}$	$\displaystyle\leq C\mathopen{\\|}\mathcal{E}({\overline{{\bf u}}}_{1}-{\overline{{\bf u}}}_{2})\mathclose{\\|}_{L^{2}(\Omega)}+C\mathopen{\\|}{\bm{\varphi}}_{1}-{\bm{\varphi}}_{2}\mathclose{\\|}_{L^{\infty}(\Omega)}.$		(3.12)

	$\displaystyle\|\chi({\bm{\varphi}}+{\bf h})-\chi({\bm{\varphi}})-\chi^{\prime}({\bm{\varphi}}){\bf h}\|$	$\displaystyle\leq\|{\bf h}\|\int_{0}^{1}\|\chi^{\prime}({\bm{\varphi}}+s{\bf h})-\chi^{\prime}({\bm{\varphi}})\|\,\mathrm{d}s\leq C_{\chi^{\prime}}\|{\bf h}\|^{2},$		(4.15)
	$\displaystyle\|\widehat{{\mathbb{C}}}({\bm{\varphi}}+{\bf h})-\widehat{{\mathbb{C}}}({\bm{\varphi}})-\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}){\bf h}\|$	$\displaystyle\leq\|{\bf h}\|\int_{0}^{1}\|\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}+s{\bf h})-\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}})\|\,\mathrm{d}s\leq C_{\widehat{{\mathbb{C}}}^{\prime}}\|{\bf h}\|^{2},$		(4.15)

	$\displaystyle\\|{\cal Y}_{1}\\|_{L^{2}(\Omega)}$	$\displaystyle\leq C\\|{\bf h}\\|_{L^{\infty}(\Omega)}\Big{(}\mathopen{\\|}{\widehat{{\bf u}}}^{{\bf h},{\bf k}}-{\widehat{{\bf u}}}^{\bf h}\mathclose{\\|}_{H^{1}(\Omega)}+\mathopen{\\|}{\widehat{{\bf u}}}^{{\bf h}}-{\widehat{{\bf u}}}\mathclose{\\|}_{H^{1}(\Omega)}+\mathopen{\\|}{\bf h}\mathclose{\\|}_{L^{\infty}(\Omega)}\mathopen{\\|}\mathcal{E}({\widehat{{\bf u}}})\mathclose{\\|}_{L^{2}(\Omega)}\Big{)}$
		$\displaystyle\leq C\mathopen{\\|}({\bf h},{\bf k})\mathclose{\\|}_{L^{\infty}(\Omega)\times H^{1}(\Omega)}^{2},$
	$\displaystyle\\|{\cal Y}_{2}\\|_{L^{2}(\Omega)}$	$\displaystyle\leq C\\|{\bf h}\\|_{L^{\infty}(\Omega)}\Big{(}\\|{\bf h}\\|_{L^{\infty}(\Omega)}\\|\mathcal{E}({\overline{{\bf u}}})\\|_{L^{2}(\Omega)}+\mathopen{\\|}{\bf h}\mathclose{\\|}_{L^{\infty}(\Omega)}\\|{\bf k}\\|_{H^{1}(\Omega)}+\\|{\bf k}\\|_{H^{1}(\Omega)}\Big{)}$
		$\displaystyle\leq C\mathopen{\\|}({\bf h},{\bf k})\mathclose{\\|}_{L^{\infty}(\Omega)\times H^{1}(\Omega)}^{2},$

	$\displaystyle-\int_{\Omega}\overline{{\mathbb{C}}}^{\prime}({\bm{\varphi}}^{})(\bm{\phi}-{\bm{\varphi}}^{})\mathcal{E}({\overline{{\bf u}}}^{*}):\mathcal{E}(\overline{{\bf p}})\,\mathrm{dx}$
	$\displaystyle\quad-\int_{\Omega}\widehat{{\mathbb{C}}}^{\prime}({\bm{\varphi}}^{})(\bm{\phi}-{\bm{\varphi}}^{})(\mathcal{E}({\widehat{{\bf u}}}^{})-\chi({\bm{\varphi}}^{})\mathcal{E}({\overline{{\bf u}}}^{*})):\mathcal{E}(\widehat{{\bf q}})\,\mathrm{dx}$
	$\displaystyle\quad+\int_{\Omega}\widehat{{\mathbb{C}}}({\bm{\varphi}}^{})\chi^{\prime}({\bm{\varphi}}^{})(\bm{\phi}-{\bm{\varphi}}^{})\mathcal{E}({\overline{{\bf u}}}^{}):\mathcal{E}(\widehat{{\bf q}})\,\mathrm{dx}$
	$\displaystyle\quad+2\gamma\varepsilon\int_{\Omega}\nabla{\bm{\varphi}}^{}\cdot\nabla(\bm{\phi}-{\bm{\varphi}}^{})\,\mathrm{dx}+\frac{\gamma}{\varepsilon}\int_{\Omega}\widetilde{\Psi}_{,{\bm{\varphi}}}({\bm{\varphi}}^{})\cdot(\bm{\phi}-{\bm{\varphi}}^{})\,\mathrm{dx}\geq 0\quad\forall\bm{\phi}\in\mathcal{U}_{\rm ad},$		(4.21)