p_{i}\rvert^{2})\theta_{n}}{\,{d}\sigma}$
$\displaystyle=-\int_{\Omega}{(\bm{\mathbf{u}}_{i}\cdot\nabla\bm{\mathbf{u}}_{i}^{\top}\bm{\mathbf{\theta}}+p_{i}\bm{\mathbf{\theta}}\cdot\nabla
p_{i})}{\,{d}x}+\frac{1}{2}\int_{\Sigma}{(\lvert\bm{\mathbf{u}}_{i}\rvert^{2}+\lvert
p_{i}\rvert^{2})\theta_{n}}{\,{d}\sigma}.$

The computations of the remaining two expressions $\dot{J}_{4}(0)$ and
$\dot{J}_{8}(0)$ need much more work. This require the utilization of the
adjoint system (43). Indeed, since
$\bm{\mathbf{y}}^{t}=\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}}\in V_{\Gamma}$ and
$r^{t}=p^{t}-p\in Q$, we can write

$\displaystyle{J}_{4}(t)+{J}_{8}(t)$
$\displaystyle\qquad=\Im\left\\{\int_{\Omega}{\alpha\nabla\overline{\bm{\mathbf{v}}}\cdot\nabla\bm{\mathbf{y}}^{t}}{\,{d}x}+i\int_{\Sigma}{(\overline{\bm{\mathbf{v}}}\cdot\bm{\mathbf{n}})(\bm{\mathbf{y}}^{t}\cdot\bm{\mathbf{n}})}{\,{d}\sigma}-\int_{\Omega}{\overline{q}\operatorname{div}\bm{\mathbf{y}}^{t}}{\,{d}x}-\int_{\Omega}{r^{t}\operatorname{div}\bm{\mathbf{v}}}{\,{d}x}\right\\}$
$\displaystyle\qquad\equiv\Im\left\\{a(\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}},\bm{\mathbf{v}})+b(\bm{\mathbf{v}},p_{t}-p)+b(\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}},q)\right\\}.$

Therefore, in light of (54) with $\bm{\mathbf{\varphi}}=\bm{\mathbf{v}}\in
V_{\Gamma}$, we get

${J}_{4}(t)+{J}_{8}(t)=\Im\\{\Phi^{t}(\bm{\mathbf{v}})+b(\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}},q)\\},$

where $\Phi^{t}(\bm{\mathbf{v}})$ is given by (55) with
$\bm{\mathbf{\varphi}}=\bm{\mathbf{v}}$. Using (27), (31), and Lemma 3.13, we
obtain the following limit

$\displaystyle\dot{J}_{4}(0)+\dot{J}_{8}(0)$ $\displaystyle=\lim_{t\to
0}\frac{1}{t}\left[\Im\\{\Phi^{t}(\bm{\mathbf{v}})+b(\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}},q)\\}\right]$
(73) $\displaystyle=\Im\Big{\\{}-\int_{\Omega}{\alpha
A\nabla{\bm{\mathbf{u}}}:\nabla{\overline{\bm{\mathbf{v}}}}}{\,{d}x}$
$\displaystyle\quad\qquad-i\int_{\Sigma}{(\operatorname{div}\bm{\mathbf{\theta}}+D\bm{\mathbf{\theta}}\bm{\mathbf{n}}\cdot\bm{\mathbf{n}})({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})({\overline{\bm{\mathbf{v}}}}\cdot\bm{\mathbf{n}})}{\,{d}\sigma}$
$\displaystyle\quad\qquad-i\int_{\Sigma}{[(-D\bm{\mathbf{\theta}}){\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}}]({\overline{\bm{\mathbf{v}}}}\cdot\bm{\mathbf{n}})}{\,{d}\sigma}-i\int_{\Sigma}{({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})[(-D\bm{\mathbf{\theta}}){\overline{\bm{\mathbf{v}}}}\cdot\bm{\mathbf{n}}]}{\,{d}\sigma}$
$\displaystyle\quad\qquad+\int_{\Omega}{(\nabla\cdot\bm{\mathbf{\theta}})p(\nabla\cdot\overline{\bm{\mathbf{v}}})}{\,{d}x}+\int_{\Omega}{p[(-D\bm{\mathbf{\theta}}):\nabla\overline{\bm{\mathbf{v}}}]}{\,{d}x}$
$\displaystyle\quad\qquad+\int_{\Omega}{[\bm{\mathbf{\nabla}}\cdot(\bm{\mathbf{f}}\otimes\bm{\mathbf{\theta}})]\cdot\overline{\bm{\mathbf{v}}}}{\,{d}x}-\int_{\Omega}{\overline{q}(D\bm{\mathbf{\theta}}:\nabla\bm{\mathbf{u}})}{\,{d}x}\Big{\\}}$
$\displaystyle=:\Im\left\\{\sum_{i=1}^{8}K_{i}\right\\},$

where $\bm{\mathbf{\nabla}}\cdot(\bm{\mathbf{f}}\otimes\bm{\mathbf{\theta}})$
is evaluated as in (32). Before we go further with our computation, we gather
in the next few lines some identities that will be useful in our calculations.
We put into use the weak formulation of the state and the adjoint state
problem given in (13) and (43), respectively. In these variational equations –
since $\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}$ and
$\bm{\mathbf{u}},\bm{\mathbf{v}}\in\bm{\mathbf{\mathsf{H}}}^{2}(\Omega)^{d}$ –
we can respectively take
$\bm{\mathbf{\varphi}}=\nabla\bm{\mathbf{v}}^{\top}\bm{\mathbf{\theta}}\in
V_{\Gamma}$ and
$\bm{\mathbf{\psi}}=\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}}\in
V_{\Gamma}$, apply integration by parts (twice), and then use the boundary
conditions on $\Sigma$ to obtain the following equations:

$\displaystyle-\int_{\Omega}{\alpha\Delta\bm{\mathbf{u}}\cdot\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}}}{\,{d}x}+\int_{\Sigma}{\alpha\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\cdot\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}}}{\,{d}\sigma}+\int_{\Omega}{\nabla
p\cdot\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}}}{\,{d}x}-\int_{\Sigma}{p\bm{\mathbf{n}}\cdot\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}}}{\,{d}\sigma}$
(74)
$\displaystyle\quad=-\int_{\Omega}{\alpha\Delta\bm{\mathbf{u}}\cdot\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}}}{\,{d}x}+\int_{\Sigma}{\alpha\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\cdot\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}}}{\,{d}\sigma}-\int_{\Omega}{p\operatorname{div}(\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}})}{\,{d}x}$
$\displaystyle\quad=\int_{\Omega}{\bm{\mathbf{f}}\cdot\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}}}{\,{d}x}-i\int_{\Sigma}{(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\cdot\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}}}{\,{d}\sigma}$

and

$\displaystyle-\int_{\Omega}{\alpha\Delta\overline{\bm{\mathbf{v}}}\cdot\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}}}{\,{d}x}+\int_{\Sigma}{\alpha\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{v}}}}\cdot\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}}}{\,{d}\sigma}+\int_{\Omega}{\nabla\overline{q}\cdot\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}}}{\,{d}x}-\int_{\Sigma}{\overline{q}\bm{\mathbf{n}}\cdot\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}}}{\,{d}\sigma}$
(75)
$\displaystyle\quad=-\int_{\Omega}{\alpha\Delta\overline{\bm{\mathbf{v}}}\cdot\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}}}{\,{d}x}+\int_{\Sigma}{\alpha\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{v}}}}\cdot\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}}}{\,{d}\sigma}-\int_{\Omega}{\overline{q}\operatorname{div}(\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}})}{\,{d}x}$
$\displaystyle\quad=\int_{\Omega}{\bm{\mathbf{u}}_{i}\cdot\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}}}{\,{d}x}-i\int_{\Sigma}{(\overline{\bm{\mathbf{v}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\cdot\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}}}{\,{d}\sigma}.$

Now, let us simplify or expand some of the integrals $K_{i}$ above.

  * •

Rewriting $K_{7}$. For the last integral, we note of the following expansion
of the integrand:

$\displaystyle[\bm{\mathbf{\nabla}}\cdot(\bm{\mathbf{f}}\otimes\bm{\mathbf{\theta}})]\cdot\overline{\bm{\mathbf{v}}}$
$\displaystyle=[\bm{\mathbf{\nabla}}\cdot(\bm{\mathbf{f}}\bm{\mathbf{\theta}}^{\top})]\cdot\overline{\bm{\mathbf{v}}}$
$\displaystyle=\sum_{j=1}^{d}\left(\bm{\mathbf{\theta}}\cdot\nabla
f_{j}+f_{j}\operatorname{div}\bm{\mathbf{\theta}}\right)\bar{v}_{j}$
$\displaystyle=\sum_{j=1}^{d}\bm{\mathbf{\theta}}\cdot\nabla
f_{j}\bar{v}_{j}+\left(\sum_{j=1}^{d}f_{j}\bar{v}_{j}\right)\operatorname{div}\bm{\mathbf{\theta}}$
$\displaystyle=\sum_{j=1}^{d}\sum_{k=1}^{d}\frac{\partial f_{j}}{\partial
x_{k}}\theta_{k}\bar{v}_{j}+\left(\bm{\mathbf{f}}\cdot\overline{\bm{\mathbf{v}}}\right)\operatorname{div}\bm{\mathbf{\theta}}$
$\displaystyle=\sum_{i=1}^{d}\sum_{k=1}^{d}\frac{\partial f_{j}}{\partial
x_{k}}\theta_{k}\bm{\mathbf{e}}_{i}\cdot\sum_{l=1}^{d}\bar{v}_{l}\bm{\mathbf{e}}_{l}+\left(\bm{\mathbf{f}}\cdot\overline{\bm{\mathbf{v}}}\right)\operatorname{div}\bm{\mathbf{\theta}}$
$\displaystyle=\nabla\bm{\mathbf{f}}^{\top}\bm{\mathbf{\theta}}\cdot\overline{\bm{\mathbf{v}}}+\left(\bm{\mathbf{f}}\cdot\overline{\bm{\mathbf{v}}}\right)\operatorname{div}\bm{\mathbf{\theta}}.$

At this juncture, we utilize the identity (71) to get the following equation

$\displaystyle\int_{\Omega}{\left[(\bm{\mathbf{f}}\cdot\overline{\bm{\mathbf{v}}})\operatorname{div}\bm{\mathbf{\theta}}+\bm{\mathbf{f}}\cdot\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}}+\nabla\bm{\mathbf{f}}^{\top}\bm{\mathbf{\theta}}\cdot\overline{\bm{\mathbf{v}}}\right]}{\,{d}x}$
$\displaystyle=\int_{\Omega}{\left(\bm{\mathbf{f}}\cdot\overline{\bm{\mathbf{v}}}\operatorname{div}\bm{\mathbf{\theta}}+\bm{\mathbf{\theta}}\cdot\nabla(\bm{\mathbf{f}}\cdot\overline{\bm{\mathbf{v}}})\right)}{\,{d}x}$
$\displaystyle=\int_{\Omega}{\operatorname{div}\left[(\bm{\mathbf{f}}\cdot\overline{\bm{\mathbf{v}}})\bm{\mathbf{\theta}}\right]}{\,{d}x}.$

Thus, by the previous computation together with the divergence theorem, we may
write $K_{7}$ as follows

$K_{7}=\int_{\Omega}{[\bm{\mathbf{\nabla}}\cdot(\bm{\mathbf{f}}\otimes\bm{\mathbf{\theta}})]\cdot\overline{\bm{\mathbf{v}}}}{\,{d}x}=\int_{\Sigma}{(\bm{\mathbf{f}}\cdot\overline{\bm{\mathbf{v}}})\theta_{n}}{\,{d}\sigma}-\int_{\Omega}{\bm{\mathbf{f}}\cdot\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}}}{\,{d}x}=:K_{71}+K_{72}.$
(76)

  * •

Rewriting $K_{5}$. Using the definition of the tangential divergence of a
vector function (see Definition 3.3) and based from the discussion given in
Step 1 of the proof of Proposition C.1 issued in Appendix C (refer to equation
(C.92)), we have

$\displaystyle\int_{\Omega}{(\nabla\cdot\bm{\mathbf{\theta}})p(\nabla\cdot\overline{\bm{\mathbf{v}}})}{\,{d}x}$
$\displaystyle\qquad=\int_{\Omega}{\\{\operatorname{div}[p(\nabla\cdot\overline{\bm{\mathbf{v}}})\bm{\mathbf{\theta}}]-(\bm{\mathbf{\theta}}\cdot\nabla)[p(\nabla\cdot\overline{\bm{\mathbf{v}}})]\\}}{\,{d}x}$
$\displaystyle\qquad=\int_{\Sigma}{p(\nabla\cdot\overline{\bm{\mathbf{v}}})\theta_{n}}{\,{d}\sigma}-\int_{\Omega}{(\bm{\mathbf{\theta}}\cdot\nabla)[p(\nabla\cdot\overline{\bm{\mathbf{v}}})]}{\,{d}x}$
$\displaystyle\qquad=\int_{\Sigma}{[p\operatorname{div}_{\Sigma}\overline{\bm{\mathbf{v}}}+\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{v}}}}\cdot(p\bm{\mathbf{n}})]\theta_{n}}{\,{d}\sigma}-\int_{\Omega}{(\bm{\mathbf{\theta}}\cdot\nabla)[p(\nabla\cdot\overline{\bm{\mathbf{v}}})]}{\,{d}x}$
$\displaystyle\qquad=\int_{\Sigma}{[\kappa
p\bm{\mathbf{n}}\theta_{n}-\nabla_{\Sigma}(p\theta_{n})]\cdot\overline{\bm{\mathbf{v}}}}{\,{d}\sigma}+\int_{\Sigma}{\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{v}}}\cdot
p\bm{\mathbf{n}}\theta_{n}}}{\,{d}\sigma}-\int_{\Omega}{(\bm{\mathbf{\theta}}\cdot\nabla)[p(\nabla\cdot\overline{\bm{\mathbf{v}}})]}{\,{d}x}$
$\displaystyle\qquad=\int_{\Sigma}{[\kappa
p\bm{\mathbf{n}}\theta_{n}-\nabla_{\Sigma}(p\theta_{n})]\cdot\overline{\bm{\mathbf{v}}}}{\,{d}\sigma}+\int_{\Sigma}{\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{v}}}\cdot
p\bm{\mathbf{n}}\theta_{n}}}{\,{d}\sigma}$
$\displaystyle\qquad\qquad\qquad-\int_{\Omega}{(\nabla\cdot\overline{\bm{\mathbf{v}}})\bm{\mathbf{\theta}}\cdot\nabla
p}{\,{d}x}-\int_{\Omega}{p(\bm{\mathbf{\theta}}\cdot\nabla)(\nabla\cdot\overline{\bm{\mathbf{v}}})}{\,{d}x}$
$\displaystyle\qquad=\int_{\Sigma}{[\kappa
p\bm{\mathbf{n}}\theta_{n}-\nabla_{\Sigma}(p\theta_{n})]\cdot\overline{\bm{\mathbf{v}}}}{\,{d}\sigma}+\int_{\Sigma}{\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{v}}}\cdot
p\bm{\mathbf{n}}\theta_{n}}}{\,{d}\sigma}$
$\displaystyle\qquad\qquad\qquad+\int_{\Omega}{(\bm{\mathbf{\theta}}\cdot\nabla
p)p_{i}}{\,{d}x}-\int_{\Omega}{p(\bm{\mathbf{\theta}}\cdot\nabla)(\nabla\cdot\overline{\bm{\mathbf{v}}})}{\,{d}x}$
$\displaystyle\qquad=K_{51}+K_{52}+K_{53}+K_{54}.$

In the last equality above (i.e., for the expression $K_{53}$), we applied the
identity

$-\int_{\Omega}{(\bm{\mathbf{\theta}}\cdot\nabla
p)(\nabla\cdot\overline{\bm{\mathbf{v}}})}{\,{d}x}=\int_{\Omega}{(\bm{\mathbf{\theta}}\cdot\nabla
p)p_{i}}{\,{d}x},$

which was obtained by taking $\mu=\bm{\mathbf{\theta}}\cdot\nabla p\in Q$ in
the second equation of (43).

  * •

Rewriting $K_{2}$. Using the identity
$\operatorname{div}_{\Sigma}\bm{\mathbf{\theta}}=\operatorname{div}\bm{\mathbf{\theta}}-D\bm{\mathbf{\theta}}\bm{\mathbf{n}}\cdot\bm{\mathbf{n}}$
on $\Sigma$ and the tangential formula, we can expand $K_{2}$ as follows:

$\displaystyle K_{2}$
$\displaystyle=-i\int_{\Sigma}{(\operatorname{div}\bm{\mathbf{\theta}}+D\bm{\mathbf{\theta}}\bm{\mathbf{n}}\cdot\bm{\mathbf{n}})({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})({\overline{\bm{\mathbf{v}}}}\cdot\bm{\mathbf{n}})}{\,{d}\sigma}$
$\displaystyle=-i\int_{\Sigma}{\operatorname{div}_{\Sigma}\bm{\mathbf{\theta}}({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})({\overline{\bm{\mathbf{v}}}}\cdot\bm{\mathbf{n}})}{\,{d}\sigma}-2i\int_{\Sigma}{D\bm{\mathbf{\theta}}\bm{\mathbf{n}}\cdot\bm{\mathbf{n}}({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})({\overline{\bm{\mathbf{v}}}}\cdot\bm{\mathbf{n}})}{\,{d}\sigma}$
$\displaystyle=-i\int_{\Sigma}{\theta_{n}\kappa({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\cdot{\overline{\bm{\mathbf{v}}}}}{\,{d}\sigma}+\int_{\Sigma}{\nabla_{\Sigma}[({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})({\overline{\bm{\mathbf{v}}}}\cdot\bm{\mathbf{n}})]\cdot\bm{\mathbf{\theta}}}{\,{d}\sigma}$
$\displaystyle\qquad-2i\int_{\Sigma}{D\bm{\mathbf{\theta}}\bm{\mathbf{n}}\cdot\bm{\mathbf{n}}({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})({\overline{\bm{\mathbf{v}}}}\cdot\bm{\mathbf{n}})}{\,{d}\sigma}$
$\displaystyle=:K_{21}+K_{22}+K_{23}.$

  * •

Rewriting $K_{1}$. For the first integral, we make use of the following
formula:

$\displaystyle-\int_{\Omega}{A\nabla\varphi\cdot\nabla\overline{\psi}}{\,{d}x}$
(77)
$\displaystyle\qquad=-\int_{\Omega}{(\Delta\varphi)\bm{\mathbf{\theta}}\cdot\nabla\overline{\psi}}{\,{d}x}-\int_{\Omega}{(\Delta\overline{\psi})\bm{\mathbf{\theta}}\cdot\nabla\varphi}{\,{d}x}+\int_{\Sigma}{\partial_{\bm{\mathbf{n}}}{\varphi}(\bm{\mathbf{\theta}}\cdot\nabla\overline{\psi})}{\,{d}\sigma}$
$\displaystyle\qquad\quad\qquad+\int_{\Sigma}{\partial_{\bm{\mathbf{n}}}{\overline{\psi}}(\bm{\mathbf{\theta}}\cdot\nabla\varphi)}{\,{d}\sigma}-\int_{\Sigma}{(\nabla\overline{\psi}\cdot\nabla\varphi)\theta_{n}}{\,{d}\sigma},$

which holds for all functions $\varphi$,
$\psi\in\bm{\mathbf{\mathsf{V}}}(\Omega)\cap\bm{\mathbf{\mathsf{H}}}^{2}(\Omega)$
and $\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}^{1}$ (see, e.g., [DZ11] or
the proof of Lemma 32 in [BP13]). Hence, with reference to equation (74) and
(LABEL:eq:adjoint_tested), we get the following computations

$\displaystyle-\int_{\Omega}{\alpha
A\nabla{\bm{\mathbf{u}}}:\nabla{\overline{\bm{\mathbf{v}}}}}{\,{d}x}=\sum_{j=1}^{d}\int_{\Omega}{\alpha
A\nabla u_{j}\cdot\nabla\overline{v}_{j}}{\,{d}x}$
$\displaystyle\qquad=\sum_{j=1}^{d}\Bigg{\\{}\int_{\Omega}{(-\alpha\Delta{u}_{j})\bm{\mathbf{\theta}}\cdot\nabla\overline{{v}}_{j}}{\,{d}x}-\int_{\Omega}{(\alpha\Delta\overline{{v}}_{j})\bm{\mathbf{\theta}}\cdot\nabla{u}_{j}}{\,{d}x}+\int_{\Sigma}{\alpha\partial_{\bm{\mathbf{n}}}{{u}_{j}}(\bm{\mathbf{\theta}}\cdot\nabla\overline{{v}}_{j})}{\,{d}\sigma}$
$\displaystyle\qquad\quad\qquad+\int_{\Sigma}{\alpha\partial_{\bm{\mathbf{n}}}{\overline{{v}}_{j}}(\bm{\mathbf{\theta}}\cdot\nabla{u}_{j})}{\,{d}\sigma}-\int_{\Sigma}{\alpha(\nabla\overline{{v}}_{j}\cdot\nabla{u}_{j})\theta_{n}}{\,{d}\sigma}\Bigg{\\}}$
$\displaystyle\qquad=\int_{\Omega}{\bm{\mathbf{f}}\cdot\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}}}{\,{d}x}-i\int_{\Sigma}{(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\cdot\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}}}{\,{d}\sigma}$
$\displaystyle\qquad\quad\qquad+\int_{\Omega}{\bm{\mathbf{u}}_{i}\cdot\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}}}{\,{d}x}-i\int_{\Sigma}{(\overline{\bm{\mathbf{v}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\cdot\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}}}{\,{d}\sigma}$
$\displaystyle\qquad\quad\qquad-\int_{\Sigma}{\alpha(\nabla_{\Sigma}\bm{\mathbf{u}}:\nabla_{\Sigma}\overline{\bm{\mathbf{v}}})\theta_{n}}{\,{d}\sigma}-\int_{\Sigma}{\alpha\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\cdot\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{v}}}}\theta_{n}}{\,{d}\sigma}$
$\displaystyle\qquad\quad\qquad+\int_{\Omega}{p\operatorname{div}(\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}})}{\,{d}x}+\int_{\Omega}{\overline{q}\operatorname{div}(\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}})}{\,{d}x}$
$\displaystyle\qquad=:K_{11}+K_{12}+K_{13}+K_{14}+K_{15}+K_{16}+K_{17}+K_{18}.$

Notice here that $K_{11}$ will cancel out with $K_{72}$. Moreover, in light of
Lemma 3.15, the sum $K_{6}+K_{54}+K_{17}$ actually equates to zero. Similarly,
using equation (69) from the said lemma, we know that

$K_{18}=\int_{\Omega}{\overline{q}\operatorname{div}(\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}})}{\,{d}x}=\int_{\Omega}{\overline{q}\left[(\bm{\mathbf{\theta}}\cdot\nabla)(\nabla\cdot\bm{\mathbf{u}})+D\bm{\mathbf{\theta}}:\nabla\bm{\mathbf{u}}\right]}{\,{d}x}.$

However, since $\operatorname{div}\bm{\mathbf{u}}=0$ in $\Omega$, then we
actually have
$K_{18}=\int_{\Omega}{\overline{q}D\bm{\mathbf{\theta}}:\nabla\bm{\mathbf{u}}}{\,{d}x}$.
This integral, in addition, also vanishes when combined with $K_{8}$.

In above, we have obtained some expression where in we can utilize Lemma 3.14.
Looking at the sum $K_{3}+K_{14}$, we can apply the identity (68) to get the
following equation

$\displaystyle K_{3}+K_{14}$
$\displaystyle=i\int_{\Sigma}{(\overline{\bm{\mathbf{v}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\cdot
D\bm{\mathbf{\theta}}\bm{\mathbf{u}}}{\,{d}\sigma}-i\int_{\Sigma}{(\overline{\bm{\mathbf{v}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\cdot\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}}}{\,{d}\sigma}$
$\displaystyle=i\int_{\Sigma}{(\overline{\bm{\mathbf{v}}}\cdot\bm{\mathbf{n}})\operatorname{div}_{\Sigma}(\theta_{n}\bm{\mathbf{u}})}{\,{d}\sigma}+i\int_{\Sigma}{(\overline{\bm{\mathbf{v}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\cdot(D\bm{\mathbf{\theta}}\bm{\mathbf{n}}\otimes\bm{\mathbf{n}})\bm{\mathbf{u}}}{\,{d}\sigma}$
$\displaystyle\qquad-i\int_{\Sigma}{(\overline{\bm{\mathbf{v}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{\theta}}\cdot\nabla_{\Sigma}(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})}{\,{d}\sigma}$
$\displaystyle=:J_{1}+J_{2}+J_{3}.$

Meanwhile, for the sum $K_{4}+K_{12}$, we have

$\displaystyle K_{4}+K_{12}$
$\displaystyle=i\int_{\Sigma}{(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\cdot
D\bm{\mathbf{\theta}}\overline{\bm{\mathbf{v}}}}{\,{d}\sigma}-i\int_{\Sigma}{(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\cdot\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}}}{\,{d}\sigma}$
$\displaystyle=i\int_{\Sigma}{(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\nabla_{\Sigma}\theta_{n}\cdot\overline{\bm{\mathbf{v}}}}{\,{d}\sigma}+i\int_{\Sigma}{(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\cdot(D\bm{\mathbf{\theta}}\bm{\mathbf{n}}\otimes\bm{\mathbf{n}})\overline{\bm{\mathbf{v}}}}{\,{d}\sigma}$
$\displaystyle\qquad\qquad\qquad-i\int_{\Sigma}{(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\bm{\mathbf{\theta}}\cdot\nabla_{\Sigma}(\overline{\bm{\mathbf{v}}}\cdot\bm{\mathbf{n}})}{\,{d}\sigma}-i\int_{\Sigma}{(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\cdot\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{v}}}}\theta_{n}}{\,{d}\sigma}$
$\displaystyle=:H_{1}+H_{2}+H_{3}+H_{4}.$

Combining this with the sum $K_{16}+K_{52}$ and noting that
$-p\bm{\mathbf{n}}+\alpha\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}+i(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}=0$
on $\Sigma$, we further get

$\displaystyle K_{4}+K_{12}+K_{16}+K_{52}$
$\displaystyle=H_{1}+H_{2}+H_{3}-\int_{\Sigma}{[-p\bm{\mathbf{n}}+\alpha\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}+i(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}]\cdot\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{v}}}}\theta_{n}}{\,{d}\sigma}$
$\displaystyle=H_{1}+H_{2}+H_{3}.$

At this point, we underline the observation that the sum $K_{22}+J_{3}+H_{3}$
vanishes. Moreover, employing Lemma 3.16, we also observe that

$\displaystyle
i\int_{\Sigma}{(\overline{\bm{\mathbf{v}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\cdot(D\bm{\mathbf{\theta}}\bm{\mathbf{n}}\otimes\bm{\mathbf{n}})\bm{\mathbf{u}}}{\,{d}\sigma}+i\int_{\Sigma}{(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\cdot(D\bm{\mathbf{\theta}}\bm{\mathbf{n}}\otimes\bm{\mathbf{n}})\overline{\bm{\mathbf{v}}}}{\,{d}\sigma}$
$\displaystyle\qquad\qquad-2i\int_{\Sigma}{D\bm{\mathbf{\theta}}\bm{\mathbf{n}}\cdot\bm{\mathbf{n}}({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})({\overline{\bm{\mathbf{v}}}}\cdot\bm{\mathbf{n}})}{\,{d}\sigma}=0,$

i.e., the sum $J_{2}+H_{2}+K_{23}$ also disappears.

Finally, summarizing our computations – and after applying the tangential
Green’s formula (see Lemma 3.4) and using the fact that
$\nabla_{\Sigma}\bm{\mathbf{u}}\bm{\mathbf{n}}=\bm{\mathbf{0}}$ (see Lemma
C.7) – we get the following equivalent expression for the sum
$\dot{J}_{2}(0)+\dot{J}_{6}(0)+\dot{J}_{4}(0)+\dot{J}_{8}(0)$ (recall
equations (73) and (73)):

$\displaystyle\dot{J}_{2}(0)+\dot{J}_{6}(0)+\dot{J}_{4}(0)+\dot{J}_{8}(0)$
$\displaystyle\quad=-\int_{\Omega}{[\bm{\mathbf{u}}_{i}\cdot\nabla\bm{\mathbf{u}}_{i}^{\top}\bm{\mathbf{\theta}}+p_{i}(\bm{\mathbf{\theta}}\cdot\nabla
p_{i})]}{\,{d}x}+\frac{1}{2}\int_{\Sigma}{(\lvert\bm{\mathbf{u}}_{i}\rvert^{2}+\lvert
p_{i}\rvert^{2})\theta_{n}}{\,{d}\sigma}$
$\displaystyle\qquad+\Im\Bigg{\\{}\int_{\Omega}{\bm{\mathbf{u}}_{i}\cdot\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}}}{\,{d}x}+\int_{\Omega}{p_{i}(\bm{\mathbf{\theta}}\cdot\nabla
p)}{\,{d}x}$
$\displaystyle\qquad\quad+\int_{\Omega}{\Big{[}\bm{\mathbf{f}}\theta_{n}-\nabla_{\Sigma}(p\theta_{n})+\operatorname{div}_{\Sigma}{[\alpha(\nabla_{\Sigma}\bm{\mathbf{u}})\theta_{n}]}+i\operatorname{div}_{\Sigma}(\theta_{n}\bm{\mathbf{u}})\bm{\mathbf{n}}+i({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})\nabla_{\Sigma}\theta_{n}\Big{]}\cdot\overline{\bm{\mathbf{v}}}}{\,{d}x}\Bigg{\\}}$
$\displaystyle\quad=\int_{\Sigma}{\left[\Im\left\\{\bm{\mathbf{B}}[\theta_{n}]\cdot\overline{\bm{\mathbf{v}}}\right\\}+\frac{1}{2}\left(|\bm{\mathbf{u}}_{i}|^{2}+|p_{i}|^{2}\right)\theta_{n}\right]}{\,{d}\sigma},$

as desired. This ends the proof of Proposition 3.6. ∎

######  Remark 3.17.

The computed expression (34) for the shape gradient of $J$ obtained via
rearrangement method clearly agrees with the structure of the same derivative
derived through the classical chain rule approach issued in Appendix C. In
fact, after some manipulations, the expression can equivalently be expressed
as

$\bm{\mathbf{B}}[\theta_{n}]=\left[\alpha\nabla\bm{\mathbf{u}}+({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}}-p)id\right]\nabla_{\Sigma}\theta_{n}-\left[-\partial_{\bm{\mathbf{n}}}{p}\bm{\mathbf{n}}+\partial_{\bm{\mathbf{n}}\bm{\mathbf{n}}}^{2}{\bm{\mathbf{u}}}+i{(\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}}\right]\theta_{n}-i\left({\bm{\mathbf{u}}}\cdot\nabla_{\Sigma}\theta_{n}\right)\bm{\mathbf{n}};$

see (C.93) in Appendix C.

The following results can be drawn easily from (34), (39), and Remark 2.1.

######  Corollary 3.18 (Necessary optimality condition).

Let the domain $\Omega^{\ast}$ be such that the state
$(\bm{\mathbf{u}},p)=(\bm{\mathbf{u}}(\Omega^{\ast}),p(\Omega^{\ast}))$, where
$\bm{\mathbf{u}}=\bm{\mathbf{u}}_{r}+i\bm{\mathbf{u}}_{i}$ and
$p=p_{r}+ip_{i}$, satisfies the free surface problem (1), i.e., we have

$-p\bm{\mathbf{n}}+\alpha\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}=\bm{\mathbf{0}}\quad\text{and}\quad\bm{\mathbf{u}}\cdot\bm{\mathbf{n}}=0\quad\text{on
$\Sigma^{\ast}$},$

or equivalently,

$\bm{\mathbf{u}}_{i}=\bm{\mathbf{0}}\quad\text{and}\quad p_{i}=0\quad\text{on
$\Omega^{\ast}$},$ (78)

with $(\bm{\mathbf{u}},p)$ satisfying (7). Then, the domain $\Omega^{\ast}$ is
stationary for the shape problem (21) (the minimization problem
$J(\Omega)=\frac{1}{2}\int_{\Omega}{\left(|\bm{\mathbf{u}}_{i}|^{2}+|p_{i}|^{2}\right)}{\,{d}x}\to\inf$,
where $(\bm{\mathbf{u}}_{i},p_{i})$ is subject to the state problem (7)). That
is, it fulfills the necessary optimality condition

${d}J(\Omega^{\ast})[\bm{\mathbf{\theta}}]=0,\quad\text{for all
$\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}^{1}$}.$ (79)

###### Proof.

By the assumption that $\bm{\mathbf{u}}_{i}=\bm{\mathbf{0}}$ and $p_{i}=0$ on
$\Omega^{\ast}$, one obtains – in view of (39) – that
$\bm{\mathbf{v}}_{i}=\bm{\mathbf{0}}$ and $q_{i}=0$ on $\Omega^{\ast}$. Thus,
it follows that $g_{\Sigma}=0$ on $\Sigma^{\ast}$ which therefore gives us the
conclusion ${d}J(\Omega^{\ast})[\bm{\mathbf{\theta}}]=0$, for any
$\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}^{1}$. ∎

In connection with the previous conclusion, we point out that the solutions of
the necessary condition (79) might exist such that the states does not satisfy
equation (78) (refer, for example to Figure 2 and Figure 3 in terms of
numerics). Nonetheless, only in the case of exact matching of boundary data a
stationary domain $\Omega^{\ast}$ is a global minimum because
$J(\Omega^{\ast})=0$.

##  4 Numerical approximation

The numerical resolution to our proposed shape optimization method to (1) is
realized through a Sobolev gradient-based method. We implement the idea of the
said approach via finite element method, along the lines of the author’s
previous work, see [RA20]. Accordingly, this section is divided into two main
subsections. In the first part, we provide some details of the algorithm and
then in the second part we give a concrete test example of the problem.

###  4.1 Numerical algorithm

For completeness, we give below the important details of our algorithm.

Choice of descent direction. As alluded above, the problem will be solved via
finite element method. To this end, we make use of the Riesz representation of
the shape gradient (34). The reason for this is that, the gradient given by
$g_{\Sigma}$ is only supported on the free boundary $\Sigma$, and if we use
directly $g_{\Sigma}\bm{\mathbf{n}}$ as the descent vector
$\bm{\mathbf{\theta}}$, the perturbation may produce some unwanted
oscillations on the free boundary which, of course, we do not want to happen.
So, we apply an extension-regularization technique by taking the descent
direction $\bm{\mathbf{\theta}}$ as the solution in
$\bm{\mathbf{H}}^{1}_{\Gamma,\bm{\mathbf{0}}}(\Omega)$ to the variational
problem

$a(\bm{\mathbf{\theta}},\bm{\mathbf{\varphi}})=-\int_{\Sigma}{g_{\Sigma}\bm{\mathbf{n}}\cdot\bm{\mathbf{\varphi}}}{\,{d}\sigma},\quad\text{for
all
$\bm{\mathbf{\varphi}}\in\bm{\mathbf{H}}^{1}_{\Gamma,\bm{\mathbf{0}}}(\Omega)$},$
(80)

where $a$ is the
${\bm{\mathbf{H}}^{1}(\Omega)}(:=H^{1}(\Omega;\mathbb{R}^{d})$)-inner product
in $d$-dimension, $d\in\\{2,3\\}$. In this manner, the Sobolev gradient
$\bm{\mathbf{\theta}}$ [Neu97] becomes a smoothed preconditioned extension of
$-g_{\Sigma}\bm{\mathbf{n}}$ over the entire domain $\Omega$. Further
discussion on discrete gradient flows for various shape optimization problems
can be found in [DMNV07].

The main algorithm. To compute the $k$th domain $\Omega^{k}$, a typical
procedure is given as follows:

1\. Initilization

Choose an initial guess for $\Sigma^{0}$ (this gives us $\Omega^{0}$).

2\. Iteration

For $k=0,1,2,\ldots$, do the following

  1. 2.1

solve the state and adjoint state systems on the current domain $\Omega^{k}$;

  2. 2.2

choose $t^{k}>0$, and compute the Sobolev gradient $\bm{\mathbf{\theta}}^{k}$
on $\Omega^{k}$;

  3. 2.3

update the current domain by setting
$\Omega_{k+1}=(\operatorname{id}+t^{k}\bm{\mathbf{\theta}}^{k})\Omega^{k}$.

3\. Stop Test

Repeat the Iteration part until convergence.

######  Remark 4.1 (Step-size computation).

In Step 2.2, the step size $t^{k}$ is computed via a backtracking line search
procedure using the formula $t^{k}=\mu
J(\Omega^{k})/|\bm{\mathbf{\theta}}^{k}|^{2}_{\mathbf{H}^{1}(\Omega^{k})}$ at
each iteration, where $\mu>0$ is a given real number. This procedure is based
on an Armijo-Goldstein-like condition for shape optimization method proposed
in [RA20, p. 281]. Here, however, the step size parameter $\mu$ is not limited
to the interval $(0,1)$ in contrast to [RA20]. We do this to maximize the
initial choice of the length of the descent vector and simply rely on the
backtracking procedure to get a suitable step size.

######  Remark 4.2 (Evaluating the mean curvature).

In Step 2.2, the mean curvature $\kappa$ of $\Sigma$ is evaluated as the
divergence of some vector $\bm{\mathbf{N}}$, where
$\bm{\mathbf{N}}\in{\bm{\mathbf{H}}^{1}(\Omega)}$ is an extension of
$\bm{\mathbf{n}}$ satisfying the equation
$c_{\text{N}}\int_{\Omega}{\nabla\bm{\mathbf{N}}:\nabla\bm{\mathbf{\varphi}}}{\,{d}x}+\int_{\Sigma}{\bm{\mathbf{N}}\cdot\bm{\mathbf{\varphi}}}{\,{d}\sigma}=\int_{\Sigma}{\bm{\mathbf{n}}\cdot\bm{\mathbf{\varphi}}}{\,{d}\sigma}$,
for all $\bm{\mathbf{\varphi}}\in{\bm{\mathbf{H}}^{1}(\Omega)}$, for some
small real number $c_{\text{N}}>0$131313This follows the same idea in
computing the Riesz representation of the shape gradient.. We point that it is
a requirement that $\bm{\mathbf{N}}$ is a unitary extension of
$\bm{\mathbf{n}}$ for the formula
$\kappa=\operatorname{div}_{\Sigma}\bm{\mathbf{n}}=\operatorname{div}\bm{\mathbf{N}}-(\nabla\bm{\mathbf{N}}\bm{\mathbf{n}})\bm{\mathbf{n}}=\operatorname{div}\bm{\mathbf{N}}$
(on $\Sigma$) to hold; see, e.g., [SZ92, Lem. 2.14, p. 92]. So, the said
extension of the (outward) unit normal vector to $\Sigma$ may not be accurate.
However, the approximation $\nabla\cdot\bm{\mathbf{N}}$ of $\kappa$ is enough
to get a good descent direction for the optimization procedure. Moreover, such
strategy easily applies to three dimensional cases. A more accurate numerical
computation of the mean curvature is of course possible, but here we are
satisfied with our results using this approximation technique for $\kappa$. In
our numerical experiments, we take $c_{\text{N}}=10^{-8}$.

######  Remark 4.3 (Stopping condition).

We terminate the algorithm after it reached a maximum number of iterations.

######  Remark 4.4 (Incorporating the shape Hessian information in the
numerical procedure).

The convergence behavior of our scheme can of course be improved by
incorporating the shape Hessian information in the numerical procedure. The
drawback, however, of a second-order method is that, typically, it demands
additional computational burden and time to carry out the calculation,
especially when the Hessian is complicated [NR00, Sim89]. Here, we will not
employ a second-order method to numerically solve the optimization problem as
the first-order method already provides an excellent approximation of the
optimal solution.

###  4.2 Numerical examples

Let us now illustrate the feasibility and applicability of the coupled complex
boundary method in the context of shape optimization for solving concrete
examples of the free surface problem (1). We first test our method in two
dimensions, and carry out a comparison with the least-squares method of
tracking the Dirichlet data (see Problem (5)). Also, since most of the
previous studies in free surface problems (see [Kas14]) only dealt with
problems in two dimensions, we also put our attention on testing the proposed
method to three dimensional cases.

######  Remark 4.5 (Details of the computational setup and environment).

The numerical simulations conducted here are all implemented in the
programming software FreeFem++, see [Hec12]. Every variational problems
involved in the procedure (except those that correspond to (7) and (39)) are
solved using $\mathbb{P}_{1}$ finite element discretization. Moreover, all
mesh deformations are carried out without any kind of adaptive mesh refinement
as opposed to what has been usually done in earlier works, see [RA20]. In this
way, we can assess the stability of CCBM in comparison with the classical
approach of tracking the Dirichlet data in least-squares sense. The
computations are all performed on a MacBook Pro with Apple M1 chip computer
having 16GB RAM processors.

######  Remark 4.6.

We note that the advantage of employing a Sobolev-gradient method is that we
can avoid the generation of a new triangulation of the domain at every
iterative step. This is achieve by moving not only the boundary nodes, but
also the internal nodes of the mesh triangulation at every iteration. Doing so
only requires one to generate an initial profile for the computational mesh.
To carry out such approach, one simply need to solve the discretized version
of (80) and then move the domain in the direction of the resulting vector
field scaled with the pseudo-time step size $t^{k}$. That is, we find
${{\bm{\mathbf{\theta}}_{h}^{k}}}\in\mathbb{P}_{1}(\Omega_{h}^{k})^{d}$ which
satisfies the following system of equations

$-\Delta{\bm{\mathbf{\theta}}_{h}^{k}}+{\bm{\mathbf{\theta}}_{h}^{k}}=\bm{\mathbf{0}}\
\ \text{in
$\Omega_{h}^{k}$},\qquad{\bm{\mathbf{\theta}}_{h}^{k}}=\bm{\mathbf{0}}\ \
\text{on
$\Gamma^{h}$},\qquad{\nabla{\bm{\mathbf{\theta}}_{h}^{k}}}\cdot{\bm{\mathbf{n}}_{h}^{k}}=-g_{\Sigma}^{k}{\bm{\mathbf{n}}_{h}^{k}}\
\ \text{on $\Sigma_{h}^{k}$},$

where we suppose a polygonal domain $\overline{\Omega_{h}^{k}}$ and its
triangulation
$\mathcal{T}_{h}(\overline{\Omega_{h}^{k}})=\\{K^{k}_{l}\\}^{N_{e}}_{l=1}$
($K^{k}_{l}$ is a closed triangle for $d=2$, or a closed tetrahedron for
$d=3$) are given, and $\mathbb{P}_{1}(\Omega_{h}^{k})^{d}$ denotes the
$\mathbb{R}^{d}$-valued piecewise linear function space on
$\mathcal{T}_{h}(\overline{\Omega_{h}^{k}})$. Accordingly, we update the
domain or, equivalently, move the nodes of the mesh by defining
$\Omega_{k+1}^{h}$ and
$\mathcal{T}_{h}(\overline{\Omega_{k+1}^{h}})=\\{K^{k+1}_{l}\\}^{N_{e}}_{l=1}$
respectively as
$\overline{\Omega^{k+1}_{h}}:=\left\\{x+t^{k}{\bm{\mathbf{\theta}}_{h}^{k}}(x)\
\middle|\ x\in\overline{\Omega_{h}^{k}}\right\\}$ and
$K^{k+1}_{l}:=\left\\{x+t^{k}{\bm{\mathbf{\theta}}_{h}^{k}}(x)\ \middle|\ x\in
K^{k}_{l}\right\\}$, for all $k=0,1,\ldots$.

We are now ready to give our first numerical example.

####  4.2.1 Example in two dimensions

In this example, we consider a gravity-like force
$\bm{\mathbf{f}}=(-10x,-10y)$ that keeps the fluid on top of the circular
domain of radius $0.4$. The fluid flowing in the domain is triggered by an
initial velocity. As the fluid comes to rest, a free surface position that is
concentric with the circular domain, is attained. On the fixed boundary
$\Gamma$, a homogenous Dirichlet boundary conditions is imposed (i.e.,
$\bm{\mathbf{g}}=\bm{\mathbf{0}}$), and we set the value of $\alpha$ to
$0.01$. Meanwhile, we set the initial location (or initial guess) of the fluid
free-boundary $\Sigma^{0}$ as follows

$\Sigma^{0}:=\left\\{(x,y)\in\mathbb{R}^{2}\mid\frac{x^{2}}{1}+\frac{y^{2}}{1.1^{2}}=1\right\\}.$

We discretize the computational domain with triangular elements generated by
the (built-in) bi-dimensional anisotropic mesh generator of FreeFEM++ [Hec12].
Here, the exterior (free) boundary is discretized with $70$ nodal points while
the interior (fixed) boundary is discretized with $30$ discretization points.
The Stokes equations (7) (as well as the adjoint system (39)) are then
discretized using the Galerkin finite-element method. We use Taylor-Hood
elements ($\mathbb{P}_{2}$-$\mathbb{P}_{1}$ finite elements) for the
approximation of the velocity and pressure. As a result, we obtain a set of
linear algebraic equations that may be written in matrix form as
$\bm{\mathbf{K}}\bar{\bm{\mathbf{u}}}=\bm{\mathbf{F}}$ where $\bm{\mathbf{K}}$
is the global system matrix, $\bar{\bm{\mathbf{u}}}$ is the global vector of
unknowns (velocities and pressures), and $\bm{\mathbf{F}}$ is a vector that
includes the effects of body forces and boundary conditions. This linear
system is solved using the default LU-solver of FreeFEM++.

The results of the experiment are shown in Figure 4.2.1 which depict a cross
comparison of the computed shapes (top figure), histories of cost values (left
upper plot), and the histories of $H^{1}$ gradient norms (left lower plot).
Observe that the solution attained by the proposed approach CCBM coincides
with the solution of the classical least-squares approach TD (tracking-
Dirichlet). Moreover, the convergence behavior of the two methods are actually
comparable as evident in the plots of the histories of cost values and $H^{1}$
gradient norms. At this point, we do not see much advantage of using the
proposed method. Nonetheless, the advantage we claim will be more apparent
when it is applied in the case of three dimension.

Meanwhile, we also plotted the initial and final imaginary parts of the Stokes
velocity and pressure profiles in Figure 2 as well as the corresponding values
for the adjoint solutions in Figure 3. Clearly, the imaginary parts reduced to
very small magnitudes as it obtains an approximation shape solution. For the
computed Stokes and adjoint pressure at the approximate optimal solution, the
(maximum) magnitudes are or order $10^{-5}$ and $10^{-7}$, respectively. This
numerically confirms Remark 2.1 and the optimality condition issued in
Corollary 3.18.

Figure 1: Cross comparison of computed shapes (left figure), histories of cost
values (left upper plot), and $H^{1}$ gradient norms (left lower plot)

Figure 2: Imaginary parts of the initial and final flow field and pressure of
the Stokes system

Figure 3: Imaginary parts of the initial and final flow field and pressure of
the adjoint system

####  4.2.2 Example in three dimensions

Now we consider a test case in three spatial dimensions. The assumptions are
similar to the previous example. That is, we consider a gravity-like force
$\bm{\mathbf{f}}=(-10x,-10y,-10z)\in\mathbb{R}^{3}$. This force is assumed to
keep the fluid to surround an object that is spherical in shape having radius
equal to $0.4$. Again, the fluid flowing in the domain is triggered by an
initial velocity. The value of $\alpha$ is again set to $0.01$ and this time
we take the initial guess to be the boundary of a three-petal flower like
shape that is given exactly in the figure shown in Figure 4.

We again look at the situations where we have coarse and fine mesh for the
computational domain. For the latter (initial) setup, the tetrahedrons have
maximum and minimum mesh width of $h_{\max}\approx 0.665$ and $h_{\min}\approx
0.047$, respectively. Meanwhile, for the test experiment with fine mesh we
take $h_{\max}\approx 0.487$ and $h_{\min}\approx 0.047$. The rest of the
computational setup is similar to the case of two dimensions.

For coarse mesh, the computational results are shown in Figure 4.2.2 – Figure
LABEL:fig:figure6d. It is evident from Figure 4.2.2 that the sequence of shape
approximations due to CCBM is different from the classical approach TD.
Nevertheless, as we observed in Figure 6, the computed optimal shapes obtained
from the two methods almost coincide. Meanwhile, Figure 7 illustrates the mesh
profile of the computed optimal shape for each method which also show a close
similarity between the two results. The next two figures, Figure 4.2.2 –
Figure 4.2.2,141414For the final pressure profile of the Stokes and the
adjoint solutions, the maximum magnitude is found to be of order $10^{-3}$ and
$10^{-4}$, respectively. depict the initial and final imaginary parts of the
Stokes’ and adjoint’s flow field and pressure profiles (magnitude) for the
coarse-mesh experiment. These numerical results – similar to the observation
made for two dimensional case – corroborate the statement issued in Remark 2.1
and the conclusion drawn in Corollary 3.18.

On the other hand, the results for a finer computational mesh are shown in
Figure 4.2.2 – Figure 12. In Figure 4.2.2 we notice a well-behaved (in the
sense that the shape approximations are smooth) evolution of the initial
domain to the optimal domain due to CCBM as opposed to the case of TD. Observe
from the latter (left plots in Figure 4.2.2), that there are obvious dents in
the computed optimal shape which we do not expect to appear. Moreover, it
becomes more evident as we made a cross comparison between the computed
optimal shapes obtained from the two methods in Figure 11. For the (initial
and) final imaginary parts of the Stokes’ and adjoint flow field and pressure
profiles (magnitude) for the case of finer mesh, the maximum magnitude is
found to be of smaller values than that of the coarse mesh, as expected.

Additionally, it is obvious from Figure 12 that the optimal shape obtained via
TD is less rounded that the one achieve using CCBM. These results clearly show
the merit of employing CCBM which we actually expect since it utilizes
(naturally as a consequence of the formulation) a volume-integral.

Finally, in Figure 13, we plot the graph of the histories of costs and
gradient norms for the two methods. It seems that, for coarse mesh (see left
plot in Figure 13), CCBM converges faster to a stationary point compared to
TD. Meanwhile, for finer mesh, it appears that graph corresponding to the
histories of cost and gradient norm values for TD stops abruptly (as seen in
the right plot of Figure 13). This is primarily due to the instability
occurring in the computational mesh which was observed in previous figures.
Based from these results, it seems that CCBM is more robust compared to TD (as
expected).

Figure 4: Mesh profile (coarse mesh) of the initial guess viewed on different
positions

Figure 5: Shape histories of the free boundary computed using TD (upper plots)
and CCBM (lower plots) with finer mesh viewed on the plane $xz$ (leftmost
column), $yx$ (middle column), and $yz$ (rightmost column)

Figure 6: Cross comparison of computed shapes (case of coarser meshes) viewed
on the plane $xz$ (leftmost plot), $yx$ (middle plot), and $yz$ (rightmost
plot)

Figure 7: Mesh profile of computed shapes (TD: left plot, CCBM: right plot)
with coarse mesh

Figure 8: Stokes’ flow field and pressure profiles (magnitude) at initial
configuration (upper two top rows) and at the computed shape (lower two bottom
rows) under coarse mesh viewed on $xz$ (leftmost column), $xy$ (middle
column), and $yz$ (rightmost column) plane.

Figure 9: Adjoint’s flow field and pressure profiles (magnitude) at initial
configuration (upper two top rows) and at the computed shape (lower two bottom
rows) under coarse mesh viewed on $xz$ (leftmost column), $xy$ (middle
column), and $yz$ (rightmost column) plane.

Figure 10: Shape histories of the free boundary computed using TD (upper
plots) and CCBM (lower plots) with finer mesh viewed on the plane $xz$
(leftmost column), $yx$ (middle column), and $yz$ (rightmost column)

Figure 11: Cross comparison of computed shapes (case of finer meshes) viewed
on the plane $xz$ (leftmost plot), $yx$ (middle plot), and $yz$ (rightmost
plot)

Figure 12: Mesh profile of computed shapes (TD: left plot, CCBM: right plot)
with finer mesh viewed on different planes

Figure 13: Histories of cost and gradient norm values for coarse (left) and
finer mesh (right)

##  5 Conclusions and Future Work

In this investigation, we have developed a coupled complex boundary method in
shape optimization framework to resolve a free surface problem involving the
Stokes equation. The shape gradient of the cost which was obtained naturally
from the proposed formulation was delicately computed under a mild regularity
assumption on the domain and without using the shape derivative of the states.
Using the shape gradient information, a Sobolev gradient-based descent scheme
was formulated in order to solve the problem numerically via finite element
method. The realization and implementation of the method and scheme put
forward in this article in carried out in two and three dimensions. Numerical
results showed that the new method has some advantages over the classical
approach of tracking the Dirichlet data in least-squares sense. Moreover, it
seems that the method is more accurate (in the sense that it achieves the
expected optimal shape solution) compared to the classical boundary tracking
method, as expected.

For future work, we propose to calculate and examine the shape Hessian of the
cost functional to investigate the ill-posedness of the proposed shape
optimization problem. The said expression could then be used in a shape Newton
method to numerically solve the minimization problem. On the other hand, an
application of CCBM in solving inverse obstacle problems under shape
optimization settings will also be a subject of our next investigation.

## References

  * [ABP+13] A. Ben Abda, F. Bouchon, G. H. Peichl, M. Sayeh, and R. Touzani.  A Dirichlet-Neumann cost functional approach for the Bernoulli problem.  J. Eng. Math., 81:157–176, 2013.
  * [AC81] A. Alt and L. A. Caffarelli.  Existence and regularity for a minimum problem with free boundary.  J. Reine. Angew. Math., 325:105–144, 1981.
  * [AD94] B. Alessandrini and G. Delhommeau.  Simulation of three-dimensional unsteady viscous free surface flow around a ship model.  Int. J. Numer. Methods Fluids, 19(4):321–342, 1994.
  * [AF03] R. A. Adams and J. J. F. Fournier.  Sobolev Spaces, volume 140 of Pure and Applied Mathematics.  Academic Press, Amsterdam, 2003.
  * [Afr22] L. Afraites.  A new coupled complex boundary method (CCBM) for an inverse obstacle problem.  Discrete Contin. Dyn. Syst. Ser. S, 15(1):23 – 40, 2022.
  * [Bab63] I. Babuška.  The theory of small changes in the domain of existence in the theory of partial differential equations and its applications.  Differential Equations and their Applications. Academic Press, New York, 1963.
  * [Bac13] J. B. Bacani.  Methods of shape optimization in free boundary problems.  PhD thesis, Karl-Franzens-Universität-Graz, Graz, Austria, 2013.
  * [BJ67] G. J. Beavers and D. D. Joseph.  Boundary conditions of a naturally permeable wall.  J. Fluid Mech., 30:197–207, 1967.
  * [BP13] J. B. Bacani and G. H. Peichl.  On the first-order shape derivative of the Kohn-Vogelius cost functional of the Bernoulli problem.  Abstr. Appl. Anal., 2013:19 pp. Article ID 384320, 2013.
  * [BPST17] F. Bouchon, G. H. Peichl, M. Sayeh, and R. Touzani.  A free boundary problem for the Stokes equation.  ESAIM Control Optim. Calc. Var., 23:195–215, 2017.
  * [BS03] E. H. Van Brummelen and A. Segal.  Numerical solution of steady free-surface flows by the adjoint optimal shape design method.  Int. J. Numer. Methods Fluids, 41(1):3–27, 2003.
  * [CGH16] X. L. Cheng, R. F. Gong, and W. Han.  A coupled complex boundary method for the cauchy problem.  Inverse Probl. Sci. Eng., 24(9):1510–1527, 2016.
  * [CGHZ14] X. L. Cheng, R. F. Gong, W. Han, and X. Zheng.  A novel coupled complex boundary method for solving inverse source problems.  Inverse Problems, page Article 055002, 2014.
  * [DL98] R. Dautray and J.-L. Lions.  Mathematical Analysis and Numerical Methods for Science and Technology, volume 2.  Springer, 1998.
  * [DMNV07] G. Doǧan, P. Morin, R.H. Nochetto, and M. Verani.  Discrete gradient flows for shape optimization and applications.  Comput. Methods Appl. Mech. Engrg., 196:3898–3914, 2007.
  * [DZ11] M. C. Delfour and J.-P. Zolésio.  Shapes and Geometries: Metrics, Analysis, Differential Calculus, and Optimization, volume 22 of Adv. Des. Control.  SIAM, Philadelphia, 2nd edition, 2011.
  * [EH09] K. Eppler and H. Harbrecht.  Tracking Neumann data for stationary free boundary problems.  SIAM J. Control Optim., 48:2901–2916, 2009.
  * [EH10] K. Eppler and H. Harbrecht.  Tracking the Dirichlet data in ${L}^{2}$ is an ill-posed problem.  J. Optim. Theory Appl., 145:17–35, 2010.
  * [EH12] K. Eppler and H. Harbrecht.  On a Kohn-Vogelius like formulation of free boundary problems.  Comput. Optim. App., 52:69–85, 2012.
  * [FR97] M. Flucher and M. Rumpf.  Bernoulli’s free-boundary problem, qualitative theory and numerical approximation.  J. Reine. Angew. Math., 486:165–204, 1997.
  * [Gal94] G. P. Galdi.  An introduction to the mathematical theory of the Navier-Stokes equations.  Springer Tracts in Natural Philosophy. Springer-Verlag, New York, 1994\.
  * [GCH17] R. Gong, X. Cheng, and W. Han.  A coupled complex boundary method for an inverse conductivity problem with one measurement.  Appl. Anal., 96(5):869–885, 2017.
  * [GR86] V. Girault and P. A. Raviart.  Finite element methods for Navier–Stokes equations.  Springer, Berlin, 1986.
  * [GT88] D. Gilbarg and N. S. Trudinger.  Elliptic Partial Differential Equations of Second Order.  Springer-Verlag, Berlin, Heidelberg, 1988.
  * [Hec12] F. Hecht.  New development in FreeFem++.  J. Numer. Math., 20:251–265, 2012.
  * [HIK+09] J. Haslinger, K. Ito, T. Kozubek, K. Kunisch, and G. H. Peichl.  On the shape derivative for problems of Bernoulli type.  Interfaces Free Bound., 11:317–330, 2009.
  * [HKKP03] J. Haslinger, T. Kozubek, K. Kunisch, and G. H. Peichl.  Shape optimization and fictitious domain approach for solving free-boundary value problems of Bernoulli type.  Comput. Optim. Appl., 26(3):231–251, 2003.
  * [HP18] A. Henrot and M. Pierre.  Shape Variation and Optimization: A Geometrical Analysis, volume 28 of Tracts in Mathematics.  European Mathematical Society, Zürich, 2018.
  * [IKP06] K. Ito, K. Kunisch, and G. H. Peichl.  Variational approach to shape derivative for a class of Bernoulli problem.  J. Math. Anal. Appl., 314(2):126–149, 2006.
  * [IKP08] K. Ito, K. Kunisch, and G. H. Peichl.  Variational approach to shape derivatives.  ESAIM Control Optim. Calc. Var., 14:517–539, 2008.
  * [JH13] C. R. Johnson and R. A. Horn.  Matrix Analysis.  Cambridge University Press, 2013.
  * [Kas14] H. Kasumba.  Shape optimization approaches to free-surface problems.  Int. J. Numer. Meth. Fluids, 74:818–845, 2014.
  * [KT99] K. T. Kärkkäinen and T. Tiihonen.  Free surfaces: shape sensitivity analysis and numerical methods.  Int. J. Numer .Methods Eng., 44(8):1079–1098., 1999.
  * [Lia01] A. Liakos.  Discretization of the Navier-Stokes equations with slip boundary condition.  Num. Meth. for Partial Diff. Eq., 1(1–18), 2001.
  * [LP12] A. Laurain and Y. Privat.  On a Bernoulli problem with geometric constraints.  ESAIM Control Optim. Calc. Var., 18:157–180, 2012.
  * [MS76] F. Murat and J. Simon.  Sur le contrôle par un domaine géométrique.  Research report 76015, Univ. Pierre et Marie Curie, Paris, 1976.
  * [Neu97] J. W. Neuberger.  Sobolev Gradients and Differential Equations.  Springer-Verlag, Berlin, 1997.
  * [NR00] A. Novruzi and J.-R. Roche.  Newton’s method in shape optimisation: a three-dimensional case.  BIT Numer. Math., 40:102–120, 2000.
  * [OCNN22] H. Ouaissa, A. Chakib, A. Nachaoui, and M. Nachaoui.  On numerical approaches for solving an inverse Cauchy Stokes problem.  Appl. Math. Optim., 85(Art. 3):37 pp., 2022.
  * [RA19a] J. F. T. Rabago and H. Azegami.  An improved shape optimization formulation of the Bernoulli problem by tracking the Neumann data.  J. Eng. Math., 117:1–29, 2019.
  * [RA19b] J. F. T. Rabago and H. Azegami.  A new energy-gap cost functional cost functional approach for the exterior Bernoulli free boundary problem.  Evol. Equ. Control Theory, 8(4):785–824, 2019.
  * [RA20] J. F. T. Rabago and H. Azegami.  A second-order shape optimization algorithm for solving the exterior Bernoulli free boundary problem using a new boundary cost functional.  Comput. Optim. Appl., 77(1):251–305, 2020.
  * [Rab22] J. F. T. Rabago.  On the new coupled complex boundary method in shape optimization framework for solving stationary free boundary problems.  Math. Control Relat. Fields, 2022.
  * [Sim80] J. Simon.  Differentiation with respect to the domain in boundary value.  Numer. Funct. Anal. Optim., 2:649–687, 1980.
  * [Sim89] J. Simon.  Second variation for domain optimization problems.  In F. Kappel, K. Kunisch, and W. Schappacher, editors, Control and Estimation of Distributed Parameter Systems, volume 91 of International Series of Numerical Mathematics, pages 361–378, Basel, 1989. Birkhäuser.
  * [SS81] H. Saito and L. E. Scriven.  Study of coating flow by the finite element method.  J. Comput. Phys., 42:53–76, 1981.
  * [SS11] S. A. Sauter and C. Schwab.  Boundary Element Methods.  Springer-Verlag, Berlin, Heidelberg, 2011.
  * [SZ92] J. Sokołowski and J.-P. Zolésio.  Introduction to Shape Optimization: Shape Sensitivity Analysis.  Springer Series in Computational Mathematics. Springer-Verlag, Berlin, Heidelberg, 1992.
  * [Tii97] T. Tiihonen.  Shape optimization and trial methods for free boundary problems.  RAIRO Modél. Math. Anal. Numér., 31:805–825, 1997.
  * [VPLG09] O. Volkov, B. Protas, W. Liao, and D. W. Glander.  Adjoint-based optimization of thermo-fluid phenomena in welding processes.  J. Eng. Math., 65(3):201–220, 2009.
  * [VRK01] E. H. VanBrummelen, H. C. Raven, and B. Koren.  Efficient numerical solution of steady free-surface Navier-Stokes flow.  J. Comput. Phys., 174(1):120–137, 2001.
  * [WSUR96] S. Wei, R. W. Smith, H. S. Udaykumar, and M. M. Rao.  Computational Fluid Dynamics with Moving Boundaries.  Taylor & Francis, Inc., Bristol, PA, USA, 1996.
  * [ZCG20] X. Zheng, X. Cheng, and R. Gong.  A coupled complex boundary method for parameter identification in elliptic problems.  Int. J. Comput. Math., 97(5):998–1015, 2020.

##  Appendix A Appendices

###  A Lemmata proofs

####  A.1 Proof of Lemma 3.1

###### Proof.

The desired results basically follows from the fact that $T_{t}$ is (at least)
a $\mathcal{C}^{1,1}$ diffeomorphism. Meanwhile, we note in particular that
the second statement of the lemma is already guaranteed by our assumption that
$A_{t}$ is bounded as underlined in (28). Nevertheless, we briefly verify our
claim as follows. Let $\varepsilon>0$ be sufficiently small and
$t\in\mathcal{I}=[0,\varepsilon]$. By reducing $\varepsilon$ if necessary, we
can assume without lost of generality that
$|tD\bm{\mathbf{\theta}}|_{\infty}<1$ for $t\in\mathcal{I}$. This permits us
to write $(DT_{t})^{-1}$ as a Neumann series:

$DT_{t}^{-1}(x)=(DT_{t})^{-1}\circ
T_{t}^{-1}(x)=\sum_{j=0}^{\infty}(-t)^{j}(D\bm{\mathbf{\theta}})^{j}(T_{t}^{-1}(x))$

for each $x\in\overline{U}$. We can then estimate its norm as follows:

$\lvert
DT_{t}^{-1}(x)\rvert_{\infty}\leqslant\sum_{j=0}^{\infty}\lvert(-t)^{j}(D\bm{\mathbf{\theta}})^{j}(T_{t}^{-1}(x))\rvert\leqslant\frac{1}{1-\varepsilon\lvert
D\bm{\mathbf{\theta}}\rvert_{\infty}}=:C,$ (A.81)

which gives us a choice for $C>0$. ∎

####  A.2 Proof of Lemma 3.2

###### Proof.

The proof of (i) for the case $p=2$ can be found in [DZ11, p. 529] (see also
the proof of [DZ11, Thm. 6.1, p. 567]). For the case $p>2$, the argumentation
is similar. On the other hand, (ii) can be obtained from triangle inequality
and the application of (i). Indeed, it is enough to show that

$\lim_{t\to 0}\|\nabla(f\circ T_{t}-f)\|_{L^{p}(U)}=\lim_{t\to
0}\|(DT_{t})^{\top}((\nabla f)\circ T_{t})-\nabla f\|_{L^{p}(U)}.$

By the triangle inequality, we know that

$\displaystyle\|(DT_{t})^{\top}((\nabla f)\circ T_{t})-\nabla f\|_{L^{p}(U)}$
$\displaystyle\qquad\qquad\leqslant\|[(DT_{t})^{\top}-id]((\nabla f)\circ
T_{t})\|_{L^{p}(U)}+\|\nabla f\circ T_{t}-\nabla f\|_{L^{p}(U)}.$

For the latter summand, we have $\|\nabla f\circ T_{t}-\nabla
f\|_{L^{p}(U)}\to 0$ as $t\to 0$ because of (i). Meanwhile, we recall that
$DT_{0}=id$ and that we have $[t\mapsto
A_{t}]\in\mathcal{C}(\mathcal{I},\mathcal{C}(\overline{U})^{d\times d})$. So,
in particular, the map $(DT_{t})^{\top}\to id$ holds in
$\mathcal{C}(\overline{U})^{d\times d}$. Thus, the first summand also vanishes
as $t\to 0$. For the proof of (iii), we may refer to [IKP06, Proof of Lemma
3.5]. Finally, (iv) is a consequence of (iii) which is easily shown as follows

$\displaystyle\frac{1}{t}\left(I_{t}f\circ
T_{t}-f\right)=\frac{1}{t}\left(I_{t}-1\right)f^{t}+\frac{1}{t}\left(f^{t}-f\right)\
\ \longrightarrow\ \
f\operatorname{div}\bm{\mathbf{\theta}}+Df\bm{\mathbf{\theta}}=\operatorname{div}(f\bm{\mathbf{\theta}})\quad\text{as
$t\to 0$}.$

This proves Lemma 3.2. ∎

###  B Computations of some identities

####  B.1 Expansion of $I_{t}$

We look at the expansion of the determinant
$I_{t}:=\operatorname{det}(DT_{t})$. Note that the Jacobian of
$T_{t}:=id+t\bm{\mathbf{\theta}}$,
$\bm{\mathbf{\theta}}:=(\theta_{1},\theta_{2},\ldots,\theta_{d})^{\top}\in\bm{\mathbf{\Theta}}^{1}$
($T_{t}:\Omega\to\mathbb{R}^{d}$, $d\in\mathbb{N}\setminus\\{1\\}$) given by
$DT_{t}=(DT_{t})_{ij}=:(M_{ij})$ has entries of the form

$M_{ij}=\delta_{ij}+t\frac{\partial V_{i}}{\partial
x_{j}}=:\delta_{ij}+tm_{ij},$

where $\delta_{ij}$ denotes the Kronecker delta function.

Let $\mathcal{S}_{d}$ be the set of all permutations of
$N_{d}:=\\{1,\ldots,d\\}$ and $\operatorname{sgn}$ be the signum of the
permutation $\sigma$ of $N_{d}$ (i.e, it is equal to $+1$ or $-1$ according to
whether the minimum number of transpositions (pairwise interchanges) necessary
to achieve it starting from $N_{d}$ is even or odd. Moreover, let
$\mathcal{I}_{d}:=\\{\sigma\in\mathcal{I}_{d}:\sigma(j)=j,j\in{N}^{\prime}_{d}\subseteq
N_{d}\\}$ and $\iota$ be the identity permutation. Then, by definition of
determinant (see, e.g., [JH13, Eq. (0.3.2.1), p. 29]), we have the following
computations

$\displaystyle I_{t}$
$\displaystyle=\sum_{\sigma\in\mathcal{S}_{d}}\left(\operatorname{sgn}\sigma\prod_{i=1}^{d}M_{i\sigma(i)}\right)$
$\displaystyle=\sum_{\sigma=\iota}\prod_{i=1}^{d}\left(1+t\frac{\partial
V_{i}}{\partial
x_{i}}\right)+\sum_{\sigma\in\mathcal{I}_{d}\setminus\\{\iota\\}}\left(\operatorname{sgn}\sigma\prod_{i=1}^{d}M_{i\sigma(i)}\right)+\sum_{\sigma\in\mathcal{S}_{d}\setminus\mathcal{I}_{d}}\left(\operatorname{sgn}\sigma\prod_{i=1}^{d}M_{i\sigma(i)}\right)$
$\displaystyle=:S_{1}+S_{2}+S_{3}.$

Observe that we may write, for some function
$\rho_{1}\in\mathcal{C}^{0}:=\mathcal{C}(\mathbb{R},\mathcal{C}^{0,1}(U))$,
the first summand as
$S_{1}=1+t\operatorname{div}\bm{\mathbf{\theta}}+t^{2}\rho_{1}(t,\bm{\mathbf{\theta}})$.
In addition, we can write the second sum as
$t^{2}\rho_{2}(t,\bm{\mathbf{\theta}})$, for some function
$\rho_{2}\in\mathcal{C}^{0}$, since each term $S_{2}$ consists of at least two
factors of $tm_{ij}$, $i\neq j$, $i,j\in N_{d}$,
$d\in\mathbb{N}\setminus\\{1\\}$. Meanwhile, all terms of $S_{3}$ have factors
of the form $tm_{ij}$, $i\neq j$. So, the sum $S_{3}$ can be expressed as
$t^{d}\rho_{\ast}(\bm{\mathbf{\theta}})$ which can be written as
$t^{2}\rho_{3}(t,\bm{\mathbf{\theta}})$, for some function
$\rho_{3}\in\mathcal{C}^{0}$. All together, we observe that
$I_{t}=1+t\operatorname{div}\bm{\mathbf{\theta}}+t^{2}\tilde{\rho}(t,\bm{\mathbf{\theta}})$,
for some function $\tilde{\rho}\in\mathcal{C}^{0}$.

####  B.2 Derivative of $B_{t}|M_{t}\mathbf{n}|^{-2}$

Let us compute the derivative of $B_{t}|M_{t}\bm{\mathbf{n}}|^{-2}$. We recall
that $B_{t}=I_{t}|M_{t}\bm{\mathbf{n}}|$ and $M_{t}=(DT_{t})^{-\top}$, and
$B_{0}=1$ and
$|M_{0}\bm{\mathbf{n}}|=|id\bm{\mathbf{n}}|=|\bm{\mathbf{n}}|=1$. We also
remember that
$\frac{d}{dt}B_{t}\big{\rvert}_{t=0}=\operatorname{div}\bm{\mathbf{\theta}}\big{\rvert}_{\Sigma}-({D}\bm{\mathbf{\theta}}\bm{\mathbf{n}})\cdot\bm{\mathbf{n}}$.

Let us consider two (column) vectors
$\bm{\mathbf{a}}:=(a_{1},a_{2},\ldots,a_{d})^{\top},\bm{\mathbf{b}}:=(b_{1},b_{2},\ldots,b_{d})^{\top}\in\mathbb{R}^{d}$.
Then, we have the following computations151515Here, the expression $O(t)$
represents a generic remainder term.

$\displaystyle\frac{d}{dt}\left(\lvert\bm{\mathbf{a}}-t\bm{\mathbf{b}}+O(t^{2})\rvert^{2}\right)$
$\displaystyle=\frac{d}{dt}\left(\sum_{i=1}^{d}(a_{i}-tb_{i}+O(t^{2}))^{2}\right)$
$\displaystyle=\sum_{i=1}^{d}2(a_{i}-tb_{i}+O(t^{2}))\cdot(-b_{i}+O(t))$
$\displaystyle=2\sum_{i=1}^{d}(-a_{i}b_{i}+O(t))=2(-\bm{\mathbf{a}}\cdot\bm{\mathbf{b}}+O(t)).$

Now, considering the expansion (compare (51))

$M_{t}\bm{\mathbf{n}}=(DT_{t})^{-\top}\bm{\mathbf{n}}=(id+tD\bm{\mathbf{\theta}})^{-\top}\bm{\mathbf{n}}=[id-t(D\bm{\mathbf{\theta}})^{\top}+\tilde{R}(t)]\bm{\mathbf{n}},$

where $\tilde{R}(t)$ is a square matrix of size $d\times d$ and is of order
$O(t^{2})$, and in view of (26) – knowing that $M_{t}\bm{\mathbf{n}}$ is
differentiable with respect to $t$ – we obtain the following calculations
based on the previous computation:

$\frac{d}{dt}\lvert
M_{t}\bm{\mathbf{n}}\rvert^{2}\Big{|}_{t=0}=\frac{d}{dt}\left(\lvert\bm{\mathbf{n}}-t(D\bm{\mathbf{\theta}})^{\top}+\tilde{R}(t)\rvert^{2}\right)\Big{|}_{t=0}=-2\bm{\mathbf{n}}\cdot[({D}\bm{\mathbf{\theta}})^{\top}\bm{\mathbf{n}}]=-2({D}\bm{\mathbf{\theta}}\bm{\mathbf{n}})\cdot\bm{\mathbf{n}}.$

This implies that

$\frac{d}{dt}\left(\frac{B_{t}}{\lvert
M_{t}\bm{\mathbf{n}}\rvert^{2}}\right)\Big{|}_{t=0}=\frac{\left(\dfrac{d}{dt}B_{t}\right)\lvert
M_{t}\bm{\mathbf{n}}\rvert^{2}-B_{t}\left(\dfrac{d}{dt}\lvert
M_{t}\bm{\mathbf{n}}\rvert^{2}\right)}{\lvert
M_{t}\bm{\mathbf{n}}\rvert^{4}}\Big{|}_{t=0}=\operatorname{div}\bm{\mathbf{\theta}}\big{\rvert}_{\Sigma}+({D}\bm{\mathbf{\theta}}\bm{\mathbf{n}})\cdot\bm{\mathbf{n}}.$

###  C Computation of the shape gradient of the cost via chain rule

To validate the expression for the shape gradient, we give below the
computation of the expression $g_{\Sigma}$ under a $\mathcal{C}^{2,1}$
regularity assumption on the domain, supposing in addition that
$\bm{\mathbf{f}}\in H_{loc}^{1}(\mathbb{R}^{d})^{d}$, specifically, we assume
$\bm{\mathbf{f}}\in H^{1}(U)^{d}$, where $U$ is a fixed convex bounded open
set in $\mathbb{R}^{d}$ such that $U\supset\overline{\Omega}$. We note that
the given regularity guarantees the existence of the material and the shape
derivative of the state and because of this, the shape gradient of the cost
can easily be established using Hadamard’s domain differentiation formula:
(see, e.g., [DZ11, Thm. 4.2, p. 483]), [HP18, eq. (5.12), Thm. 5.2.2, p. 194]
or [SZ92, eq. (2.168), p. 113]):

$\displaystyle\left.\left\\{\frac{{d}}{{d}t}\int_{\Omega_{t}}{f(t,x)}{\,{d}x_{t}}\right\\}\right|_{t=0}=\int_{\Omega}{\frac{\partial}{\partial
t}f(0,x)}{\,{d}x}+\int_{\partial\Omega}{f(0,\sigma)\theta_{n}}{\,{d}\sigma}$
(C.82)

(and, of course, with the assumption that the perturbation of $\Omega$
preserves its regularity).

The next result is a restatement of Proposition 3.6, but with higher
regularity assumption.

######  Proposition C.1.

Let $\Omega\in\mathcal{C}^{2,1}$ and
$\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}^{2}$. Then, the shape derivative
of $J$ at $\Omega$ along $\bm{\mathbf{\theta}}$ is given by
$dJ(\Omega)[\bm{\mathbf{\theta}}]=\int_{\Sigma}{g_{\Sigma}{\theta_{n}}}{\,{d}\sigma}$,
where $g_{\Sigma}$ is the expression in (34).

Before we prove the above proposition, we briefly prepare the following
lemmata which will be useful in the derivation of the shape derivative.

######  Lemma C.2.

Let $\bm{\mathbf{n}}$ be the outward unit normal to $\Sigma$ and
$\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}^{2}$. Then, the shape derivative
of $\bm{\mathbf{n}}$ denoted by $\bm{\mathbf{n}}^{\prime}$ is given by

$\bm{\mathbf{n}}^{\prime}=-\nabla_{\Sigma}\theta_{n}.$ (C.83)

###### Proof.

See part of the proof of Proposition 5.4.14 in [HP18, Chap. 5., Sec. 4.4, p.
222]. ∎

######  Lemma C.3.

Let $\Omega\in\bm{\mathbf{\Theta}}^{2}$ and $\bm{\mathbf{n}}$ be the outward
unit normal to $\Sigma$. Then,

$D\tilde{\bm{\mathbf{n}}}\bm{\mathbf{n}}=(D\tilde{\bm{\mathbf{n}}})^{\top}\bm{\mathbf{n}}=(\nabla\tilde{\bm{\mathbf{n}}})\bm{\mathbf{n}}=0\quad\text{on}\
\Sigma,$

where $\tilde{\bm{\mathbf{n}}}$ is a unitary $\mathcal{C}^{1}$ extension of
the vector field $\bm{\mathbf{n}}$ on $\Sigma$.

###### Proof.

Because $\Omega$ is $\mathcal{C}^{2,1}$ regular, then by Proposition 5.4.8 of
[HP18, p. 218] (see also [GT88, Lem. 16.1, p. 390]), there exists a
$\mathcal{C}^{1}$ unitary extension
$\tilde{\bm{\mathbf{n}}}:=(\tilde{n}_{1},\ldots,\tilde{n}_{d})^{\top}$161616The
same is used in Lemma C.5 and Lemma C.6. of $\bm{\mathbf{n}}$. So, in an open
neighborhood of $\Sigma$, we have
$|\tilde{\bm{\mathbf{n}}}|^{2}=\langle\tilde{\bm{\mathbf{n}}},\tilde{\bm{\mathbf{n}}}\rangle=1$.
Thus, for each $j=1,\ldots,d$, we have the following computation

$\dfrac{\partial}{\partial
x_{j}}(|\tilde{\bm{\mathbf{n}}}|^{2})=\dfrac{\partial}{\partial
x_{j}}\left(\sum_{i=1}^{d}\tilde{n}_{i}^{2}\right)=2\sum_{i=1}^{d}\tilde{n}_{i}\dfrac{\partial\tilde{n}_{i}}{\partial
x_{j}}=0,$

or equivalently,

$\nabla(|\tilde{\bm{\mathbf{n}}}|^{2})=\nabla\left(\sum_{i=1}^{d}\tilde{n}_{i}^{2}\right)=2\sum_{i=1}^{d}\tilde{n}_{i}\nabla\tilde{n}_{i}=2\sum_{i=1}^{d}\tilde{n}_{i}\left(\sum_{j=1}^{d}\frac{\partial\tilde{n}_{i}}{\partial
x_{j}}\bm{\mathbf{e}}_{j}\right)=\bm{\mathbf{0}},$ (C.84)

where
$\bm{\mathbf{e}}_{j}:=(\overbrace{0,\ldots,0,1}^{j},0,\ldots,0)^{\top}\in\mathbb{R}^{d}$
is the $j$th vector of the canonical basis in $\mathbb{R}^{d}$. Now, on the
other hand, we have

$\displaystyle(D\tilde{\bm{\mathbf{n}}})^{\top}\tilde{\bm{\mathbf{n}}}=\begin{pmatrix}\dfrac{\partial\tilde{n}_{1}}{\partial
x_{1}}&\ldots&\dfrac{\partial\tilde{n}_{d}}{\partial x_{1}}\\\
\vdots&\ddots&\vdots\\\ \dfrac{\partial\tilde{n}_{1}}{\partial
x_{d}}&\ldots&\dfrac{\partial\tilde{n}_{d}}{\partial
x_{d}}\end{pmatrix}\begin{pmatrix}\tilde{n}_{1}\\\ \vdots\\\
\tilde{n}_{d}\end{pmatrix}=\begin{pmatrix}\displaystyle\sum_{j=1}^{d}\dfrac{\partial\tilde{n}_{j}}{\partial
x_{1}}\tilde{n}_{j}\\\ \vdots\\\
\displaystyle\sum_{j=1}^{d}\dfrac{\partial\tilde{n}_{j}}{\partial
x_{d}}\tilde{n}_{j}\end{pmatrix}=\sum_{k=1}^{d}\left(\sum_{j=1}^{d}\dfrac{\partial\tilde{n}_{j}}{\partial
x_{k}}\tilde{n}_{j}\right)\bm{\mathbf{e}}_{k}.$ (C.85)

Combining (C.84) and (C.85), we deduce that

$\langle(D\tilde{\bm{\mathbf{n}}})\bm{\mathbf{n}},\bm{\mathbf{n}}\rangle=\langle\bm{\mathbf{n}},(D\tilde{\bm{\mathbf{n}}})^{\top}\bm{\mathbf{n}}\rangle=0\quad\text{on}\
\Gamma,$

from which we infer the conclusion. ∎

######  Remark C.4.

In light of the previous lemma, and recalling the definition of the tangential
Jacobian matrix from Definition 3.3, we see that for the $\mathcal{C}^{1}$
unitary extension $\tilde{\bm{\mathbf{n}}}$ of $\bm{\mathbf{n}}$, we clearly
have the identity
$D\tilde{\bm{\mathbf{n}}}\big{|}_{\Sigma}=D\bm{\mathbf{n}}=D_{\Sigma}\bm{\mathbf{n}}=(D_{\Sigma}\bm{\mathbf{n}})^{\top}$
(refer to [DZ11, Chap. 9, Sec. 5.2, eq. (5.17) – (5.19), p. 497], for the last
equation).

######  Lemma C.5.

Let $\Omega\in\bm{\mathbf{\Theta}}^{2}$ and $\bm{\mathbf{n}}$ be the outward
unit normal to $\Sigma$. Then, for the solution $\bm{\mathbf{u}}$ of (7), we
have

$\frac{\partial(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})}{\partial\bm{\mathbf{n}}}=\frac{\partial\bm{\mathbf{u}}}{\partial\bm{\mathbf{n}}}\cdot\bm{\mathbf{n}}\qquad\text{on
$\Sigma$}.$

###### Proof.

First, we note that $\dfrac{\partial
n_{j}}{\partial\bm{\mathbf{n}}}=\bm{\mathbf{n}}^{\top}\nabla
n_{j}=\displaystyle\sum_{k=1}^{d}n_{k}\dfrac{\partial n_{j}}{\partial x_{k}}.$
Moreover, from the proof of the previous lemma, we know that

$\displaystyle(D\tilde{\bm{\mathbf{n}}})\tilde{\bm{\mathbf{n}}}=\sum_{k=1}^{d}\left(\sum_{j=1}^{d}\dfrac{\partial\tilde{n}_{k}}{\partial
x_{j}}\tilde{n}_{j}\right)\bm{\mathbf{e}}_{k}=\bm{\mathbf{0}}$ (C.86)

from which it can be deduced that

$\sum_{j=1}^{d}{n}_{j}\dfrac{\partial{n}_{k}}{\partial x_{j}}=0\quad\text{on
$\Sigma$},$

for all $k=1,\ldots,d$. In addition, we have the following identity

$\frac{\partial\bm{\mathbf{u}}}{\partial\bm{\mathbf{n}}}\cdot\bm{\mathbf{n}}=(D\bm{\mathbf{u}}){\bm{\mathbf{n}}}\cdot\bm{\mathbf{n}}=\sum_{k=1}^{d}\left(\sum_{j=1}^{d}\dfrac{\partial
u_{k}}{\partial
x_{j}}{n}_{j}\right)\bm{\mathbf{e}}_{k}\cdot\sum_{i=1}^{d}n_{i}\bm{\mathbf{e}}_{i}=\sum_{k=1}^{d}\left(\sum_{j=1}^{d}\dfrac{\partial
u_{k}}{\partial x_{j}}{n}_{j}\right)n_{k},$ (C.87)

where the last equation follows from the fact that
$\bm{\mathbf{e}}_{i}\cdot\bm{\mathbf{e}}_{j}=\delta_{ij}$ and $\delta_{ij}$ is
the Kronecker delta function. Thus, we have the following computations:

$\displaystyle\frac{\partial(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})}{\partial\bm{\mathbf{n}}}$
$\displaystyle=\frac{\partial}{\partial\bm{\mathbf{n}}}\left(\sum_{k=1}^{d}u_{k}n_{k}\right)=\sum_{k=1}^{d}\left(\frac{\partial
u_{k}}{\partial\bm{\mathbf{n}}}n_{k}+u_{k}\frac{\partial
n_{k}}{\partial\bm{\mathbf{n}}}\right)=\sum_{k=1}^{d}\left(\sum_{j=1}^{d}n_{j}\dfrac{\partial
u_{k}}{\partial x_{j}}\right)n_{k}.$

The desired identity then follows by comparing the previous equation with
(C.87). ∎

######  Lemma C.6.

Let $\Omega\in\bm{\mathbf{\Theta}}^{2}$ and $\bm{\mathbf{n}}$ be the outward
unit normal to $\Sigma$. Then, for the solution $\bm{\mathbf{u}}$ of (7), we
have

$\frac{\partial}{\partial\bm{\mathbf{n}}}\left[(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\right]=\left(\frac{\partial\bm{\mathbf{u}}}{\partial\bm{\mathbf{n}}}\cdot\bm{\mathbf{n}}\right)\bm{\mathbf{n}}\qquad\text{on
$\Sigma$}.$

###### Proof.

Let $\tilde{\bm{\mathbf{n}}}$ be the $\mathcal{C}^{1}$ extension of
$\bm{\mathbf{n}}$ as before and denote
$\bm{\mathbf{b}}=a\tilde{\bm{\mathbf{n}}}:=\left(\sum_{i=1}^{d}u_{i}\tilde{n}_{i}\right)\tilde{\bm{\mathbf{n}}}$
with
$\bm{\mathbf{b}}:=(b_{1},\ldots,b_{d})=(a\tilde{n}_{1},\ldots,a\tilde{n}_{d})$.
Note that

$(D\bm{\mathbf{b}})\tilde{\bm{\mathbf{n}}}=\sum_{k=1}^{d}\left(\sum_{j=1}^{d}\dfrac{\partial
b_{k}}{\partial x_{j}}\tilde{n}_{j}\right)\bm{\mathbf{e}}_{k}.$ (C.88)

Moreover, we have the following computations

$\displaystyle\frac{\partial b_{k}}{\partial
x_{j}}=\frac{\partial(a\tilde{n}_{k})}{\partial
x_{j}}=\left(\sum_{i=1}^{d}\left(\frac{\partial u_{i}}{\partial
x_{j}}\tilde{n}_{i}+u_{i}\frac{\partial\tilde{n}_{i}}{\partial
x_{j}}\right)\right)\tilde{n}_{k}+\left(\sum_{i=1}^{d}u_{i}\tilde{n}_{i}\right)\frac{\partial\tilde{n}_{k}}{\partial
x_{j}}.$

Inserting the above expression to (C.88) and using Lemma C.5 (see also
(C.86)), we get

$\displaystyle(D\bm{\mathbf{b}})\tilde{\bm{\mathbf{n}}}$
$\displaystyle=\sum_{k=1}^{d}\left(\sum_{j=1}^{d}\left[\left(\sum_{i=1}^{d}\frac{\partial
u_{i}}{\partial
x_{j}}\tilde{n}_{i}\right)\tilde{n}_{k}\right]\tilde{n}_{j}\right)\bm{\mathbf{e}}_{k}+\sum_{k=1}^{d}\left(\sum_{j=1}^{d}\left[\left(\sum_{i=1}^{d}u_{i}\tilde{n}_{i}\right)\frac{\partial\tilde{n}_{k}}{\partial
x_{j}}\right]\tilde{n}_{j}\right)\bm{\mathbf{e}}_{k}$
$\displaystyle=\sum_{k=1}^{d}\left(D\bm{\mathbf{u}}\tilde{\bm{\mathbf{n}}}\cdot\bm{\mathbf{n}}\right)\tilde{n}_{k}\bm{\mathbf{e}}_{k}+(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\sum_{k=1}^{d}\left(\sum_{j=1}^{d}\frac{\partial\tilde{n}_{k}}{\partial
x_{j}}\tilde{n}_{j}\right)\bm{\mathbf{e}}_{k}$
$\displaystyle=\left(D\bm{\mathbf{u}}\tilde{\bm{\mathbf{n}}}\cdot\tilde{\bm{\mathbf{n}}}\right)\tilde{\bm{\mathbf{n}}}.$

Thus, on $\Sigma$, we get

$\frac{\partial}{\partial\bm{\mathbf{n}}}\left[(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\right]=\left(\frac{\partial\bm{\mathbf{u}}}{\partial\bm{\mathbf{n}}}\cdot\bm{\mathbf{n}}\right)\bm{\mathbf{n}},$

as desired. ∎

######  Lemma C.7.

Let $\Omega$ be a sufficiently smooth domain with boundary
$\Gamma:=\partial\Omega$. Assume that
$\bm{\mathbf{\varphi}}\in\bm{\mathbf{W}}^{1,1}(\partial\Omega)$. Then, the
following identity holds (see, e.g., [DZ11, Chap. 9, Sec. 5.4, eq. (5.20), p.
497]):

$\nabla_{\Gamma}\bm{\mathbf{\varphi}}\bm{\mathbf{n}}=\bm{\mathbf{0}}.$

###### Proof.

Assume that $\bm{\mathbf{\varphi}}\in\bm{\mathbf{W}}^{1,1}(\Gamma)$. By
definition
$\nabla_{\Gamma}\bm{\mathbf{\varphi}}:=\nabla\bm{\mathbf{\varphi}}-(\nabla\bm{\mathbf{\varphi}}\bm{\mathbf{n}})\otimes\bm{\mathbf{n}}$,
so we have

$\nabla_{\Gamma}\bm{\mathbf{\varphi}}\bm{\mathbf{n}}=\left[\nabla\bm{\mathbf{\varphi}}-(\nabla\bm{\mathbf{\varphi}}\bm{\mathbf{n}})\otimes\bm{\mathbf{n}}\right]\bm{\mathbf{n}}=\nabla\bm{\mathbf{\varphi}}\bm{\mathbf{n}}-\left[(\nabla\bm{\mathbf{\varphi}}\bm{\mathbf{n}})\otimes\bm{\mathbf{n}}\right]\bm{\mathbf{n}}.$

Computing
$\left[(\nabla\bm{\mathbf{\varphi}}\bm{\mathbf{n}})\otimes\bm{\mathbf{n}}\right]\bm{\mathbf{n}}$,
we get

$\displaystyle\left[(\nabla\bm{\mathbf{\varphi}}\bm{\mathbf{n}})\otimes\bm{\mathbf{n}}\right]\bm{\mathbf{n}}$
$\displaystyle=\left[\left(\sum_{j=1}^{d}\left(\sum_{i=1}^{d}\dfrac{\partial\varphi_{i}}{\partial
x_{j}}n_{i}\right)\bm{\mathbf{e}}_{j}\right)\otimes\sum_{k=1}^{d}n_{k}\bm{\mathbf{e}}_{k}\right]\sum_{l=1}^{d}n_{l}\bm{\mathbf{e}}_{l}$
$\displaystyle=\sum_{j=1}^{d}\left(\sum_{i=1}^{d}\dfrac{\partial\varphi_{i}}{\partial
x_{j}}n_{i}\right)\bm{\mathbf{e}}_{j}\sum_{k=1}^{d}\sum_{l=1}^{d}n_{l}n_{k}\bm{\mathbf{e}}_{k}^{\top}\bm{\mathbf{e}}_{l}$
$\displaystyle=\sum_{j=1}^{d}\left(\sum_{i=1}^{d}\dfrac{\partial\varphi_{i}}{\partial
x_{j}}n_{i}\right)\bm{\mathbf{e}}_{j}\sum_{k=1}^{d}\sum_{l=1}^{d}n_{l}n_{k}\delta_{kl}$
$\displaystyle=\sum_{j=1}^{d}\left(\sum_{i=1}^{d}\dfrac{\partial\varphi_{i}}{\partial
x_{j}}n_{i}\right)\bm{\mathbf{e}}_{j}$
$\displaystyle=\nabla\bm{\mathbf{\varphi}}\bm{\mathbf{n}},$

from which the desired result clearly follows. ∎

###### Proof of Proposition C.1.

Let us assume that $\Omega$ is of class $\mathcal{C}^{2,1}$ and
$\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}^{2}$. Using classical regularity
theory, we have $\bm{\mathbf{u}}_{r}$, $\bm{\mathbf{u}}_{i}\in
H^{3}(\Omega)^{d}$ and $p_{r}$, $p_{i}\in H^{2}(\Omega)$. Because we have
sufficient regularity for $\bm{\mathbf{u}}$, $\Omega$, and
$\bm{\mathbf{\theta}}$, then we can apply formula (C.82) to obtain – noting
that $\bm{\mathbf{\theta}}=\bm{\mathbf{0}}$ on $\Gamma$ – the derivative

${{d}}J(\Omega)[\bm{\mathbf{\theta}}]=\int_{\Omega}{\left(\bm{\mathbf{u}}_{i}\bm{\mathbf{u}}_{i}^{\prime}+p_{i}p_{i}^{\prime}\right)}{\,{d}x}+\frac{1}{2}\int_{\Sigma}{\left(|\bm{\mathbf{u}}_{i}|^{2}+|p_{i}|^{2}\right)\theta_{n}}{\,{d}\sigma}=:\mathbb{I}_{1}+\mathbb{I}_{2}.$
(C.89)

Hereafter, we proceed in four steps:

Step 1.

We establish the strong form of the shape derivatives
$\bm{\mathbf{u}}^{\prime}$ and $p^{\prime}$ which is characterized by the
complex PDE system (C.98).

Step 2.

We prove the differentiability of $J(\Omega)$ in the direction of
$(\delta\tilde{\bm{\mathbf{u}}},\delta\tilde{p})\in X\times Q$.

Step 3.

We justify the structure of the adjoint system (39).

Step 4.

We obtain the expression for the shape gradient via the adjoint method.

Step 1. We recall the variational equation (13). Because $\Omega$,
$\bm{\mathbf{u}}$, $p$, and $\bm{\mathbf{\theta}}$ are regular enough, then
the solution $(\bm{\mathbf{u}},p)\in X\times Q$ is shape differentiable and we
can differentiate (13) (formally) to get

$\displaystyle\int_{\Omega}{\alpha\nabla{\bm{\mathbf{u}}^{\prime}}:\nabla{\overline{\bm{\mathbf{\psi}}}}}{\,{d}x}+i\int_{\Sigma}{({\bm{\mathbf{u}}^{\prime}}\cdot\bm{\mathbf{n}})({\overline{\bm{\mathbf{\psi}}}}\cdot\bm{\mathbf{n}})}{\,{d}\sigma}-\int_{\Omega}{p^{\prime}(\nabla\cdot\overline{\bm{\mathbf{\psi}}})}{\,{d}x}-\int_{\Omega}{\bar{\lambda}(\nabla\cdot{\bm{\mathbf{u}}^{\prime}})}{\,{d}x}$
(C.90)
$\displaystyle\qquad=-\int_{\Sigma}{\alpha(\nabla{\bm{\mathbf{u}}}:\nabla{\overline{\bm{\mathbf{\psi}}}}-p(\nabla\cdot\overline{\bm{\mathbf{\psi}}}))\theta_{n}}{\,{d}\sigma}$
$\displaystyle\qquad\qquad-i\int_{\Sigma}{({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}}^{\prime})({\overline{\bm{\mathbf{\psi}}}}\cdot\bm{\mathbf{n}})}{\,{d}\sigma}-i\int_{\Sigma}{({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})({\overline{\bm{\mathbf{\psi}}}}\cdot\bm{\mathbf{n}}^{\prime})}{\,{d}\sigma}$
$\displaystyle\qquad\qquad-i\int_{\Sigma}{\left[\frac{\partial}{\partial\bm{\mathbf{n}}}\left(({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\right)+\kappa\left({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\right)\right]\cdot\overline{\bm{\mathbf{\psi}}}\theta_{n}}{\,{d}\sigma}+\int_{\Sigma}{\bm{\mathbf{f}}\cdot\overline{\bm{\mathbf{\psi}}}\theta_{n}}{\,{d}\sigma},$

where the shape derivative $\bm{\mathbf{n}}^{\prime}$ of the normal vector
$\bm{\mathbf{n}}$ is given by (C.83).

From the previous equation we can derive a boundary value problem for
$(\bm{\mathbf{u}}^{\prime},p^{\prime})$. Namely, choosing
$\bm{\mathbf{\psi}}\in\bm{\mathbf{\mathcal{C}}}^{\infty}_{0}(\Omega)^{d}$ and
$\lambda\in\bm{\mathbf{\mathcal{C}}}^{\infty}_{0}(\Omega)$171717For clarity,
we mention here that there is slight abuse of notations. Specifically,
$\bm{\mathbf{\psi}}\in\bm{\mathbf{\mathcal{C}}}^{\infty}_{0}(\Omega)^{d}$ and
$\lambda\in\bm{\mathbf{\mathcal{C}}}^{\infty}_{0}(\Omega)$ means that
$\bm{\mathbf{\psi}}_{r},\bm{\mathbf{\psi}}_{i}\in\mathcal{C}^{\infty}_{0}(\Omega)^{d}$
and $\lambda_{r},\lambda_{i}\in\mathcal{C}^{\infty}_{0}(\Omega)$ where
$\mathcal{C}^{\infty}_{0}(\Omega)^{d}$ and $\mathcal{C}^{\infty}_{0}(\Omega)$
denotes the usual space of infinitely differentiable (vector-valued and
scalar-valued) functions, respectively. reveals that (via integration by
parts)

$-\alpha\Delta\bm{\mathbf{u}}^{\prime}+\nabla
p^{\prime}=\bm{\mathbf{0}}\quad\text{in
$\Omega$}\qquad\text{and}\qquad\nabla\cdot\bm{\mathbf{u}}^{\prime}=\bm{\mathbf{0}}\quad\text{in
$\Omega$},$ (C.91)

which hold in the distributional sense. That is, we have obtained the first
equation above from

$\langle-\alpha\Delta\bm{\mathbf{u}}^{\prime}+\nabla
p^{\prime},\overline{\bm{\mathbf{\psi}}}\rangle_{[\bm{\mathbf{\mathcal{C}}}^{\infty}_{0}(\Omega)^{d}]^{\prime},\bm{\mathbf{\mathcal{C}}}^{\infty}_{0}(\Omega)^{d}}=\int_{\Omega}{\left(\alpha\nabla{\bm{\mathbf{u}}^{\prime}}:\nabla{\overline{\bm{\mathbf{\psi}}}}-p^{\prime}(\nabla\cdot\overline{\bm{\mathbf{\psi}}})\right)}{\,{d}x}=0,$

and the second one from the argument, i.e., we check that
$\nabla\cdot\bm{\mathbf{u}}^{\prime}=\bm{\mathbf{0}}$ in
$[\bm{\mathbf{\mathcal{C}}}^{\infty}_{0}(\Omega)^{d}]^{\prime}$.

Meanwhile, because $\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}^{2}$, i.e.,
$\bm{\mathbf{\theta}}$ vanishes on $\Gamma$, the boundary condition on
$\Gamma$ satisfied by $\bm{\mathbf{u}}^{\prime}$ easily follows. That is, we
have $\bm{\mathbf{u}}^{\prime}=\bm{\mathbf{0}}$ on $\Gamma$.

Now, for the succeeding arguments, we underline here the fact that
$\Omega\in\mathcal{C}^{2,1}$, and since we have the
$\bm{\mathbf{\mathsf{H}}}^{3}(\Omega)^{d}\times H^{2}(\Omega)$ regularity for
$(\bm{\mathbf{u}},p)$, we know that
$\nabla\bm{\mathbf{u}}^{\prime}\in\bm{\mathbf{L}}^{2}(\Omega)^{d}$ and $\nabla
p^{\prime}\in\bm{\mathbf{L}}^{2}(\Omega)^{d}$.

We next exhibit the boundary condition on $\Sigma$. We choose181818In fact,
here, we choose a test function
$\bm{\mathbf{\psi}}\in\bm{\mathbf{\mathsf{H}}}^{2}(\Omega)^{d}$, and because
$\Omega\in\mathcal{C}^{2,1}$, it follows that – by Stein’s extension theorem
(see, e.g., [AF03, Thm. 5.24, p. 154]) – we can construct an extension of
$\bm{\mathbf{\psi}}$ in $\bm{\mathbf{\mathsf{H}}}^{2}(\mathbb{R}^{d})^{d}$
(still denoted by $\bm{\mathbf{\psi}}$).
$\bm{\mathbf{\varphi}}\in\bm{\mathbf{\mathcal{C}}}^{\infty}(\Sigma)^{d}$ and
$\mu\in\bm{\mathbf{\mathcal{C}}}^{\infty}(\Sigma)$. Accordingly, we can find
an extension
$\bm{\mathbf{\psi}}\in\bm{\mathbf{\mathcal{C}}}^{\infty}(\Omega)^{d}$ and
$\lambda\in\bm{\mathbf{\mathcal{C}}}^{\infty}(\Omega)$ such that
$\bm{\mathbf{\psi}}\big{|}_{\Sigma}=\bm{\mathbf{\varphi}}$ and
$\lambda\big{|}_{\Sigma}=\mu$, and
$\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{\psi}}}\big{|}_{\Sigma}=\bm{\mathbf{0}}$
and $\partial_{\bm{\mathbf{n}}}{\lambda}\big{|}_{\Sigma}=0$. Applying
integration by parts in the left hand side of (C.90) and using (C.91), we get

$\displaystyle\int_{\Sigma}{\left[\alpha\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}^{\prime}}+i({\bm{\mathbf{u}}^{\prime}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}-p^{\prime}\bm{\mathbf{n}}\right]\cdot\overline{\bm{\mathbf{\psi}}}}{\,{d}\sigma}$
$\displaystyle\qquad\qquad\qquad=-\int_{\Sigma}{\left[\alpha\nabla{\bm{\mathbf{u}}}:\nabla{\overline{\bm{\mathbf{\psi}}}}-p(\nabla\cdot\overline{\bm{\mathbf{\psi}}})\right]\theta_{n}}{\,{d}\sigma}+\int_{\Sigma}{\bm{\mathbf{B}}_{1}[\theta_{n}]\cdot{\overline{\bm{\mathbf{\psi}}}}}{\,{d}\sigma},$

where

$\bm{\mathbf{B}}_{1}[\theta_{n}]:=-i[({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}}^{\prime})\bm{\mathbf{n}}+({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}^{\prime}]-i\left[\partial_{\bm{\mathbf{n}}}{\left(({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\right)}+\kappa({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\right]\theta_{n}+\bm{\mathbf{f}}\theta_{n}.$

At this juncture, let us recall the definition of tangential gradient and
divergence from subsection 3.2, to deduce the identity

$\int_{\Sigma}{p\nabla\cdot\overline{\bm{\mathbf{\psi}}}\theta_{n}}{\,{d}\sigma}=\int_{\Sigma}{\left[p\operatorname{div}_{\Sigma}\overline{\bm{\mathbf{\psi}}}+\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{\psi}}}}\cdot(p\bm{\mathbf{n}})\right]\theta_{n}}{\,{d}\sigma},$

and make use of the formulas (cf. equation in remark found in [SZ92, p. 87])

$\displaystyle(\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\otimes\bm{\mathbf{n}}):(\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{\psi}}}}\otimes\bm{\mathbf{n}})$
$\displaystyle=\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\cdot\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{\psi}}}};$
$\displaystyle\nabla\bm{\mathbf{u}}:(\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{\psi}}}}\otimes\bm{\mathbf{n}})$
$\displaystyle=\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\cdot\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{\psi}}}}=\nabla\overline{\bm{\mathbf{\psi}}}:(\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\otimes\bm{\mathbf{n}}),$

to obtain

$\nabla\bm{\mathbf{u}}:\nabla\overline{\bm{\mathbf{\psi}}}=\nabla_{\Sigma}\bm{\mathbf{u}}:\nabla_{\Sigma}\overline{\bm{\mathbf{\psi}}}+\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\cdot\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{\psi}}}}.$

Moreover, using the tangential Green’s formula (see Lemma 3.4) and applying
integration by parts on the boundary $\Sigma$, we can write

$\displaystyle\int_{\Sigma}{p\theta_{n}\operatorname{div}_{\Sigma}\overline{\bm{\mathbf{\psi}}}}{\,{d}\sigma}$
$\displaystyle=\int_{\Sigma}{\left[\kappa
p\theta_{n}\bm{\mathbf{n}}-\nabla_{\Sigma}(p\theta_{n})\right]\cdot\overline{\bm{\mathbf{\psi}}}}{\,{d}\sigma}=:\int_{\Sigma}{\bm{\mathbf{B}}_{2}[\theta_{n}]\cdot\overline{\bm{\mathbf{\psi}}}}{\,{d}\sigma},$
(C.92)
$\displaystyle-\int_{\Sigma}{\alpha(\nabla_{\Sigma}{\bm{\mathbf{u}}}:\nabla_{\Sigma}{\overline{\bm{\mathbf{\psi}}}})\theta_{n}}{\,{d}\sigma}$
$\displaystyle=\int_{\Sigma}{\overline{\bm{\mathbf{\psi}}}\cdot\left\\{\operatorname{div}_{\Sigma}{[\alpha(\nabla_{\Sigma}\bm{\mathbf{u}})\theta_{n}]}-[\alpha\kappa(\nabla_{\Sigma}\bm{\mathbf{u}})\theta_{n}\bm{\mathbf{n}}]\right\\}}{\,{d}\sigma}$
$\displaystyle=\int_{\Sigma}{\overline{\bm{\mathbf{\psi}}}\cdot\left\\{\operatorname{div}_{\Sigma}{[\alpha(\nabla_{\Sigma}\bm{\mathbf{u}})\theta_{n}]}\right\\}}{\,{d}\sigma}=:\int_{\Sigma}{\overline{\bm{\mathbf{\psi}}}\cdot\bm{\mathbf{B}}_{3}[\theta_{n}]}{\,{d}\sigma},$

where the last equality is due to the fact that
$\nabla_{\Sigma}\bm{\mathbf{u}}\bm{\mathbf{n}}=\bm{\mathbf{0}}$, see Lemma
C.7. Using these identities, together with the fact that
$-\alpha\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}+p\bm{\mathbf{n}}=i(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}$
and $\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{\psi}}}=\bm{\mathbf{0}}$ on
$\Sigma$, we get

$\displaystyle\int_{\Sigma}{\left(\alpha\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}^{\prime}}+i({\bm{\mathbf{u}}^{\prime}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}-p^{\prime}\bm{\mathbf{n}}\right)\cdot\overline{\bm{\mathbf{\psi}}}}{\,{d}\sigma}$
$\displaystyle\qquad=-\int_{\Sigma}{\alpha(\nabla_{\Sigma}{\bm{\mathbf{u}}}:\nabla_{\Sigma}{\overline{\bm{\mathbf{\psi}}}}-p(\nabla_{\Sigma}\cdot\overline{\bm{\mathbf{\psi}}}))\theta_{n}}{\,{d}\sigma}+i\int_{\Sigma}{(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\theta_{n}\cdot\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{\psi}}}}}{\,{d}\sigma}+\int_{\Sigma}{\bm{\mathbf{B}}_{1}\cdot{\overline{\bm{\mathbf{\psi}}}}}{\,{d}\sigma}$
$\displaystyle\qquad=\int_{\Sigma}{\bm{\mathbf{B}}[\theta_{n}]\cdot{\overline{\bm{\mathbf{\psi}}}}}{\,{d}\sigma},$

where

$\displaystyle\bm{\mathbf{B}}[\theta_{n}]$
$\displaystyle:=\bm{\mathbf{B}}_{1}[\theta_{n}]+\bm{\mathbf{B}}_{2}[\theta_{n}]+\bm{\mathbf{B}}_{3}[\theta_{n}]$
$\displaystyle\
=\left\\{-i[({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}}^{\prime})\bm{\mathbf{n}}+({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}^{\prime}]-i\left[\partial_{\bm{\mathbf{n}}}{\left(({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\right)}+\kappa({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\right]\theta_{n}+\bm{\mathbf{f}}\theta_{n}\right\\}$
$\displaystyle\qquad+\left\\{\kappa
p\theta_{n}\bm{\mathbf{n}}-\nabla_{\Sigma}(p\theta_{n})\right\\}+\left\\{\operatorname{div}_{\Sigma}{[\alpha(\nabla_{\Sigma}\bm{\mathbf{u}})\theta_{n}]}\right\\}$
$\displaystyle\
=\bm{\mathbf{f}}\theta_{n}-\nabla_{\Sigma}(p\theta_{n})+\operatorname{div}_{\Sigma}{[\alpha(\nabla_{\Sigma}\bm{\mathbf{u}})\theta_{n}]}+i\left[\left({\bm{\mathbf{u}}}\cdot\nabla_{\Sigma}\theta_{n}\right)\bm{\mathbf{n}}+({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})\nabla_{\Sigma}\theta_{n}\right]$
$\displaystyle\qquad-\left[i\left(\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}}\right)\bm{\mathbf{n}}+i\kappa({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}-\kappa
p\bm{\mathbf{n}}\right]\theta_{n},$

where the latter equation follows from Lemma C.6 and Lemma C.2.

At this juncture, we comment that
$\bm{\mathbf{B}}[\theta_{n}]\in\bm{\mathbf{H}}^{1/2}(\Sigma)^{d}$. Indeed,
since $\Omega$ is of class $\mathcal{C}^{2,1}$, the (outward) unit normal
vector $\bm{\mathbf{n}}$ is $\mathcal{C}^{1,1}(N^{\varepsilon})$ regular and
$\kappa\in\mathcal{C}^{0,1}(N^{\varepsilon})\subset
W^{1,\infty}(N^{\varepsilon})\subset H^{1}(N^{\varepsilon})\subset
H^{1/2}(\partial\Omega)$, where $N^{\varepsilon}$ is a small neighborhood of
$\partial\Omega$ (see, e.g., [DZ11, Sec. 7.8]). Therefore, (using the density
in $\bm{\mathbf{L}}^{2}(\Omega)^{d}$ of the traces on $\Sigma$ of functions in
$\bm{\mathbf{\mathsf{H}}}^{2}(\Omega)^{d}$) we get

$\displaystyle\alpha\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}^{\prime}}+i({\bm{\mathbf{u}}^{\prime}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}-p^{\prime}\bm{\mathbf{n}}=\bm{\mathbf{B}}[\theta_{n}],\qquad\text{on
$\Sigma$}.$

In the next several lines, we simplify the expression
$\bm{\mathbf{B}}[\theta_{n}]$. To do this, we recall some identities from
Lemma 3.4 and the definitions of some tangential operators, specially the
Laplace-Beltrami operator and its decomposition (see (33)). Also, we note
that, in $\Omega$, we have $-\alpha\Delta\bm{\mathbf{u}}+\nabla
p=\bm{\mathbf{f}}$. So, we get the following equations on $\Sigma$:

$\displaystyle\nabla_{\Sigma}(p\theta_{n})$
$\displaystyle=\theta_{n}\nabla_{\Sigma}p+p\nabla_{\Sigma}(\theta_{n}),$
$\displaystyle\operatorname{div}_{\Sigma}(\alpha\nabla_{\Sigma}\bm{\mathbf{u}}\theta_{n})$
$\displaystyle=\alpha\theta_{n}\operatorname{div}_{\Sigma}(\nabla_{\Sigma}\bm{\mathbf{u}})+\alpha\nabla_{\Sigma}\bm{\mathbf{u}}\nabla_{\Sigma}\theta_{n}$
$\displaystyle=\theta_{n}(\alpha\Delta_{\Sigma}\bm{\mathbf{u}})+\alpha\nabla\bm{\mathbf{u}}\nabla_{\Sigma}\theta_{n}$
$\displaystyle=-\theta_{n}(\bm{\mathbf{f}}+\kappa\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}+\partial_{\bm{\mathbf{n}}\bm{\mathbf{n}}}^{2}{\bm{\mathbf{u}}}-\nabla
p)+\alpha\nabla\bm{\mathbf{u}}\nabla_{\Sigma}\theta_{n}.$

Because
$p\bm{\mathbf{n}}-\alpha\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}=i(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}$
on $\Sigma$, the latter equation further implies that

$\displaystyle\kappa
p\bm{\mathbf{n}}\theta_{n}+\left[\bm{\mathbf{f}}\theta_{n}+\operatorname{div}_{\Sigma}(\alpha\nabla_{\Sigma}\bm{\mathbf{u}}\theta_{n})\right]$
$\displaystyle=\kappa
p\bm{\mathbf{n}}\theta_{n}+\left[-\theta_{n}(\kappa\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}+\partial_{\bm{\mathbf{n}}\bm{\mathbf{n}}}^{2}{\bm{\mathbf{u}}}-\nabla
p)+\alpha\nabla\bm{\mathbf{u}}\nabla_{\Sigma}\theta_{n}\right]$
$\displaystyle=i\kappa(\bm{\mathbf{u}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\theta_{n}-\theta_{n}(\partial_{\bm{\mathbf{n}}\bm{\mathbf{n}}}^{2}{\bm{\mathbf{u}}}-\nabla
p)+\alpha\nabla\bm{\mathbf{u}}\nabla_{\Sigma}\theta_{n}.$

These computations, together with Lemma C.6 and the identities

$\nabla
p-\nabla_{\Sigma}p=\partial_{\bm{\mathbf{n}}}{p}\bm{\mathbf{n}}\quad\text{and}\quad\nabla_{\Sigma}\theta_{n}\cdot\bm{\mathbf{u}}-\theta_{n}(D\bm{\mathbf{u}}\bm{\mathbf{n}}\cdot\bm{\mathbf{n}})=\\!\nabla_{\Sigma}\theta_{n}\cdot\bm{\mathbf{u}}+\theta_{n}\operatorname{div}_{\Sigma}\bm{\mathbf{u}}=\operatorname{div}_{\Sigma}(\theta_{n}\bm{\mathbf{u}}),$

on $\Sigma$, lead us to conclude that the expression
$\bm{\mathbf{B}}[\theta_{n}]$ can be written equivalently in the following
ways:

$\displaystyle\bm{\mathbf{B}}[\theta_{n}]$
$\displaystyle=\bm{\mathbf{f}}\theta_{n}-\nabla_{\Sigma}(p\theta_{n})+\operatorname{div}_{\Sigma}{[\alpha(\nabla_{\Sigma}\bm{\mathbf{u}})\theta_{n}]}+i\operatorname{div}_{\Sigma}(\theta_{n}\bm{\mathbf{u}})\bm{\mathbf{n}}+i({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})\nabla_{\Sigma}\theta_{n}$
(C.93)
$\displaystyle\qquad-\kappa\left[-p\bm{\mathbf{n}}+i({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}\right]\theta_{n}$
$\displaystyle=\left[\alpha\nabla\bm{\mathbf{u}}+({\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}}-p)id\right]\nabla_{\Sigma}\theta_{n}-\left[-\partial_{\bm{\mathbf{n}}}{p}\bm{\mathbf{n}}+\partial_{\bm{\mathbf{n}}\bm{\mathbf{n}}}^{2}{\bm{\mathbf{u}}}+i{(\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}}\right]\theta_{n}$
$\displaystyle\qquad-i\left({\bm{\mathbf{u}}}\cdot\nabla_{\Sigma}\theta_{n}\right)\bm{\mathbf{n}}.$

In summary, the shape derivative $(\bm{\mathbf{u}}^{\prime},p^{\prime})$ of
$(\bm{\mathbf{u}},p)$ which solves equation (7) is given by

$\left\\{\begin{array}[]{rcll}-\alpha\Delta\bm{\mathbf{u}}^{\prime}+\nabla
p^{\prime}&=&\bm{\mathbf{0}}&\quad\text{in $\Omega$},\\\
\nabla\cdot\bm{\mathbf{u}}^{\prime}&=&\bm{\mathbf{0}}&\quad\text{in
$\Omega$},\\\ \bm{\mathbf{u}}^{\prime}&=&\bm{\mathbf{0}}&\quad\text{on
$\Gamma$},\\\
-p^{\prime}\bm{\mathbf{n}}+\alpha\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}^{\prime}}+i(\bm{\mathbf{u}}^{\prime}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}&=&\bm{\mathbf{B}}[\theta_{n}]&\quad\text{on
$\Sigma$},\end{array}\right.$ (C.98)

where $\bm{\mathbf{B}}[\theta_{n}]$ is given by (C.93).

Step 2. Because we have sufficient regularity on the unknown variables, the
directional derivatives (see, e.g., [DZ11, Chap. 2, Sec. 2])
$J^{\prime}(\Omega)\delta\tilde{\bm{\mathbf{u}}}$ and
$J^{\prime}(\Omega)\delta\tilde{p}$ exist and are easily computed as

$\displaystyle J^{\prime}(\Omega)\delta\tilde{\bm{\mathbf{u}}}$
$\displaystyle=\frac{d}{d\epsilon}\left.\left(\frac{1}{2}\int_{\Omega}{\left(|\bm{\mathbf{u}}_{i}+\epsilon\delta\tilde{\bm{\mathbf{u}}}|^{2}+|p_{i}|^{2}\right)}{\,{d}x}\right)\right|_{\epsilon=0}=\int_{\Omega}{\bm{\mathbf{u}}_{i}\delta\tilde{\bm{\mathbf{u}}}}{\,{d}x},$
$\displaystyle J^{\prime}(\Omega)\delta\tilde{p}$
$\displaystyle=\frac{d}{d\epsilon}\left.\left(\frac{1}{2}\int_{\Omega}{\left(|\bm{\mathbf{u}}_{i}|^{2}+|p_{i}+\epsilon\delta\tilde{p}|^{2}\right)}{\,{d}x}\right)\right|_{\epsilon=0}=\int_{\Omega}{p_{i}\delta\tilde{p}}{\,{d}x}.$

Step 3. Now, to derive and justify the structure of the adjoint system (39),
we let $\bm{\mathbf{\bar{x}}}:=(\bm{\mathbf{u}},p)\in\mathcal{X}:=X\times Q$,
$\bm{\mathbf{\xi}}:=(\bm{\mathbf{\psi}},\mu)\in\mathcal{X}$,
$\bm{\mathbf{\bar{y}}}:=(\bm{\mathbf{v}},q)\in\mathcal{X}$, and define the
operator
$\mathcal{E}_{\bm{\mathbf{\bar{x}}}}({\bm{\mathbf{\bar{x}}}},\Omega)\bm{\mathbf{\xi}}\in\mathscr{L}(X,X^{\ast})$
where (cf. (40))

$\langle\mathcal{E}_{\bm{\mathbf{\bar{x}}}}({\bm{\mathbf{\bar{x}}}},\Omega)\bm{\mathbf{\xi}},\bm{\mathbf{\bar{y}}}\rangle_{\mathcal{X}^{\ast},\mathcal{X}}:=\tilde{a}({\bm{\mathbf{\psi}}},{\bm{\mathbf{v}}})+b(\bm{\mathbf{v}},\mu)+b(\bm{\mathbf{\psi}},q).$

The above operator is bijective if and only if for every
$\bm{\mathbf{\varphi}}\in X^{\ast}$ and $\lambda\in Q^{\ast}$, there exists a
unique solution $\bm{\mathbf{\bar{y}}}:=(\bm{\mathbf{v}},q)\in\mathcal{X}$ to
the equation (cf. equation (43))

$\langle\mathcal{E}_{\bm{\mathbf{\bar{x}}}}({\bm{\mathbf{\bar{x}}}},\Omega)\bm{\mathbf{\xi}},\bm{\mathbf{\bar{y}}}\rangle_{\mathcal{X}^{\ast},\mathcal{X}}=(\bm{\mathbf{\varphi}},\bm{\mathbf{\psi}})+(\mu,\lambda),$

for all $\bm{\mathbf{\xi}}\in\mathcal{X}$. Similar to what we have noted in
Remark 3.7, the existence of a unique solution to the above equation – with
$\bm{\mathbf{\varphi}}=\bm{\mathbf{u}}_{i}\in V_{\Gamma}$ and
$\lambda=p_{i}\in Q$ – can be established using similar arguments issued for
the well-posedness of the state problem (7) (see subsection 2.2). So, we omit
the details of the proof.

Since $(\bm{\mathbf{u}},p)\in(V_{\Gamma}\times
Q)\cap(\bm{\mathbf{\mathsf{H}}}^{3}(\Omega)^{d}\times H^{2}(\Omega))$, then by
the previous step together with Remark 3.7, there exists a unique adjoint
state $(\bm{\mathbf{v}},q)\in(V_{\Gamma}\times
Q)\cap(\bm{\mathbf{\mathsf{H}}}^{3}(\Omega)^{d}\times H^{2}(\Omega))$ such
that the variational equation (43) holds. By standard arguments as in Step 1 –
applying integration by parts and/or Green’s formula – we recover (39).

Step 4. Now, we are ready to express the shape derivative of $J$ computed in
(C.89) via the adjoint method – eliminating the shape derivative of the states
$\bm{\mathbf{u}}^{\prime}$ and $p^{\prime}$ appearing in $\mathbb{I}_{1}$. To
do this, let us first consider the following weak formulation of (C.98): find
$({\bm{\mathbf{u}}^{\prime}},p^{\prime})\in V_{\Gamma}\times Q$ such that

$\left\\{\begin{array}[]{rcll}a({\bm{\mathbf{u}}^{\prime}},{\bm{\mathbf{\varphi}}})+b(\bm{\mathbf{\varphi}},p^{\prime})&=&\displaystyle\int_{\Sigma}{\bm{\mathbf{B}}[\theta_{n}]\cdot\overline{\bm{\mathbf{\varphi}}}}{\,{d}\sigma},&\quad\forall{\bm{\mathbf{\varphi}}}\in{V_{\Gamma}},\\\\[5.0pt]
b(\bm{\mathbf{u}}^{\prime},\lambda)&=&0,&\quad\forall\lambda\in{Q}.\\\
\end{array}\right.$ (C.101)

Following similar arguments carried out in subsection 2.2, the existence of
(unique) weak solution to the above problem is a consequence of the complex
version of the Lax-Milgram lemma [DL98, Thm. 1, p. 376] (see also [SS11, Lem.
2.1.51, p. 40]).

Now, by taking $(\bm{\mathbf{\varphi}},\lambda)=(\bm{\mathbf{v}},q)\in
V_{\Gamma}\times Q$, we get

$\left\\{\begin{array}[]{rcl}a({\bm{\mathbf{u}}^{\prime}},\bm{\mathbf{v}})+b(\bm{\mathbf{v}},p^{\prime})&=&\displaystyle\int_{\Sigma}{\bm{\mathbf{B}}[\theta_{n}]\cdot\overline{\bm{\mathbf{v}}}}{\,{d}\sigma},\\\\[5.0pt]
b(\bm{\mathbf{u}}^{\prime},q)&=&0.\\\ \end{array}\right.$

On the other hand, let us take
$(\bm{\mathbf{\psi}},\mu)=(\bm{\mathbf{u}}^{\prime},p^{\prime})\in
V_{\Gamma}\times Q$ in the variational equation (43) of the adjoint system
(39). This lead us to

$\left\\{\begin{array}[]{rcll}\tilde{a}({\bm{\mathbf{v}}},\bm{\mathbf{u}}^{\prime})+b(\bm{\mathbf{u}}^{\prime},q)&=&\tilde{F}(\bm{\mathbf{u}}^{\prime}),\\\
b(\bm{\mathbf{v}},p^{\prime})&=&(p^{\prime},p_{i}).\end{array}\right.$

Taking the complex conjugate of both sides of the equations above, and then
combining it with the previous two, yields the following identity

$\int_{\Omega}{\left(\bm{\mathbf{u}}_{i}\bm{\mathbf{u}}^{\prime}+p_{i}p^{\prime}\right)}{\,{d}x}=\int_{\Sigma}{\bm{\mathbf{B}}[\theta_{n}]\cdot\overline{\bm{\mathbf{v}}}}{\,{d}\sigma}.$

Comparing the imaginary parts of both sides of the above equation gives us

$\mathbb{I}_{1}=\int_{\Omega}{\left(\bm{\mathbf{u}}_{i}\bm{\mathbf{u}}_{i}^{\prime}+p_{i}p_{i}^{\prime}\right)}{\,{d}x}=\Im\left\\{\int_{\Omega}{\left(\bm{\mathbf{u}}_{i}\bm{\mathbf{u}}^{\prime}+p_{i}p^{\prime}\right)}{\,{d}x}\right\\}=\Im\left\\{\int_{\Sigma}{\bm{\mathbf{B}}[\theta_{n}]\cdot\overline{\bm{\mathbf{v}}}}{\,{d}\sigma}\right\\}.$

In conclusion, the shape derivative of $J$ at $\Omega$ in the direction of the
vector filed $\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}^{2}$ is given by

${{d}}J(\Omega)[\bm{\mathbf{\theta}}]=\Im\left\\{\int_{\Sigma}{\bm{\mathbf{B}}[\theta_{n}]\cdot\overline{\bm{\mathbf{v}}}}{\,{d}\sigma}\right\\}+\frac{1}{2}\int_{\Sigma}{\left(|\bm{\mathbf{u}}_{i}|^{2}+|p_{i}|^{2}\right)\theta_{n}}{\,{d}\sigma},$

as desired. ∎