Numerical solution to a free boundary problem for the Stokes equation using the coupled complex boundary method in shape optimization settings

TL;DR

This paper introduces a novel boundary reformulation for free boundary Stokes problems using complex Robin boundary conditions, enabling shape optimization via a domain variation method with numerical validation in 2D and 3D.

Contribution

It proposes a new complex boundary reformulation and cost functional for free boundary Stokes problems, with explicit shape derivatives and a shape optimization algorithm.

Findings

Effective in 2D and 3D numerical tests

Outperforms classical Dirichlet tracking methods

Provides explicit shape derivatives for free boundary identification

Abstract

A new reformulation of a free boundary problem for the Stokes equations governing a viscous flow with overdetermined condition on the free boundary is proposed. The idea of the method is to transform the governing equations to a boundary value problem with a complex Robin boundary condition coupling the two boundary conditions on the free boundary. The proposed formulation give rise to a new cost functional that apparently has not been exploited yet in the literature, specifically, and at least, in the context of free surface problems. The shape derivatives of the cost function constructed by the imaginary part of the solution in the whole domain in order to identify the free boundary is explicitly determined. Using the computed shape gradient information, a domain variation method from a preconditioned steepest descent algorithm is applied to solve the shape optimization problem. Numerical results illustrating the applicability of the method is then provided both in two and three spatial dimensions. For validation and evaluation of the method, the numerical results are compared with the ones obtained via the classical tracking Dirichlet data.

Equations

$$\left\{\begin{array}[]{rcll}-\alpha\Delta\bm{\mathbf{u}}+\nabla p&=&\bm{\mathbf{f}}&\quad\text{in $\Omega$},\\ \nabla\cdot\bm{\mathbf{u}}&=&\bm{\mathbf{0}}&\quad\text{in $\Omega$},\\ \bm{\mathbf{u}}&=&\bm{\mathbf{g}}&\quad\text{on $\Gamma$},\\ -p\bm{\mathbf{n}}+\alpha\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}&=&\bm{\mathbf{0}}&\quad\text{on $\Sigma$},\\ \bm{\mathbf{u}}\cdot\bm{\mathbf{n}}&=&0&\quad\text{on $\Sigma$},\end{array}\right.$$

$$J_{KV}(\Omega ):=\frac{1}{2}{\int }_{\Omega }{\left|\nabla \left({\mathbf{u}}_{\text{N}}-{\mathbf{u}}_{\text{D}}\right)\right|}^{2}dx\to \inf ,$$

$$\left\{\begin{array}[]{rcll}-\alpha\Delta\bm{\mathbf{u}}_{D}+\nabla p_{D}&=&\bm{\mathbf{f}}&\quad\text{in $\Omega$},\\ \nabla\cdot\bm{\mathbf{u}}_{D}&=&\bm{\mathbf{0}}&\quad\text{in $\Omega$},\\ \bm{\mathbf{u}}_{D}&=&\bm{\mathbf{g}}&\quad\text{on $\Gamma$},\\ \bm{\mathbf{u}}_{D}\cdot\bm{\mathbf{n}}&=&\bm{\mathbf{0}}&\quad\text{on $\Sigma$},\\ \alpha\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}_{D}}\cdot\bm{\mathbf{\tau}}&=&\bm{\mathbf{0}}&\quad\text{on $\Sigma$},\end{array}\right.$$

$$\left\{\begin{array}[]{rcll}-\alpha\Delta\bm{\mathbf{u}}_{N}+\nabla p_{N}&=&\bm{\mathbf{f}}&\quad\text{in $\Omega$},\\ \nabla\cdot\bm{\mathbf{u}}_{N}&=&\bm{\mathbf{0}}&\quad\text{in $\Omega$},\\ \bm{\mathbf{u}}_{N}&=&\bm{\mathbf{g}}&\quad\text{on $\Gamma$},\\ -p_{N}\bm{\mathbf{n}}+\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}_{N}}&=&\bm{\mathbf{0}}&\quad\text{on $\Sigma$}.\end{array}\right.$$

$$J_D(\Omega )$$

$$J_N(\Omega )$$

$$J_{L^2}(\Omega ):=\frac{1}{2}{\int }_{\Omega }|{\mathbf{u}}_{\text{N}}-{\mathbf{u}}_{\text{D}}{|}^{2}dx.$$

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

$$\left\{\begin{array}[]{rcll}-\alpha\Delta\bm{\mathbf{u}}_{r}+\nabla p_{r}&=&\bm{\mathbf{f}}&\quad\text{in $\Omega$},\\ \nabla\cdot\bm{\mathbf{u}}_{r}&=&\bm{\mathbf{0}}&\quad\text{in $\Omega$},\\ \bm{\mathbf{u}}_{r}&=&\bm{\mathbf{g}}&\quad\text{on $\Gamma$},\\ -p_{r}\bm{\mathbf{n}}+\alpha\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}_{r}}&=&(\bm{\mathbf{u}}_{i}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}&\quad\text{on $\Sigma$},\end{array}\right.$$

$$\left\{\begin{array}[]{rcll}-\alpha\Delta\bm{\mathbf{u}}_{i}+\nabla p_{i}&=&\bm{\mathbf{0}}&\quad\text{in $\Omega$},\\ \nabla\cdot\bm{\mathbf{u}}_{i}&=&\bm{\mathbf{0}}&\quad\text{in $\Omega$},\\ \bm{\mathbf{u}}_{i}&=&\bm{\mathbf{0}}&\quad\text{on $\Gamma$},\\ -p_{i}\bm{\mathbf{n}}+\alpha\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}_{i}}&=&-(\bm{\mathbf{u}}_{r}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}&\quad\text{on $\Sigma$},\end{array}\right.$$

$$\nabla :=\left(\begin{array}{c}\frac{\partial }{\partial x_1}\\ ⋮\\ \frac{\partial }{\partial x_d}\end{array}\right)\text{and}\frac{\partial }{\partial \mathbf{n}}=\mathbf{n}\cdot \nabla =\left(\begin{array}{c}{n}_1\\ ⋮\\ n_d\end{array}\right)\cdot \left(\begin{array}{c}\frac{\partial }{\partial x_1}\\ ⋮\\ \frac{\partial }{\partial x_d}\end{array}\right)={\mathbf{n}}^{\top }\nabla =\sum _{i=1}^d{n}_i\frac{\partial }{\partial x_i},$$

$$W^{m,2}(\Omega )=\{u\in L^2(\Omega )\mid D^{\beta }u\in L^2(\Omega ),\text{for 0⩽∣β∣⩽m}\},$$

$$\Vert \cdot {\Vert }_{H^m(\Omega )}=\sum _{|\beta |⩽m}{\int }_{\Omega }|D^{\beta }u{|}^{2}dx.$$

$$H^m(\Omega {)}^d:=\{\mathbf{u}=(u_1,u_2,\dots ,u_d{)}^{\top }\mid u_i\in H^m(\Omega ),\text{for i=1,2,…,d}\}.$$

$$\Vert \mathbf{u}{\Vert }_{H^m(\Omega {)}^d}={\left(\sum _{i=1}^d\Vert u_i{\Vert }_{H^m(\Omega )}^2\right)}^{1/2}.$$

$$\text{for all u,v∈Hm(Ω)d},((\mathbf{u},\mathbf{v}){)}_{m,\Omega ,d}=\sum _{j=1}^d(\mathbf{u},\overline{\mathbf{v}}{)}_{m,\Omega }\text{and}{\left|\left|\left|\mathbf{v}\right|\right|\right|}_{m,\Omega ,d}=\sqrt{((\mathbf{v},\mathbf{v}){)}_{m,\Omega ,d}}.$$

$$X$$

$$Q$$

$$\text{ensurestackMath}{\text{stackon}[0.5pt]V\circ }_0$$

$$=\left\{\mathbf{v}\in \text{ensurestackMath}\text{stackon}[0.5pt]V\circ \mid (\nabla \cdot \mathbf{v},q)=0,\forall q\in Q\right\}$$

$$=\left\{\mathbf{v}\in \text{ensurestackMath}\text{stackon}[0.5pt]V\circ \mid \nabla \cdot \mathbf{v}=0\right\},$$

$$\text{ensurestackMath}{\text{stackon}[0.5pt]V\circ }_{\perp }$$

$$\nabla \mathbf{v}:\nabla \overline{\mathbf{w}}=\left(\sum _{i,j=1}^d\frac{\partial {\mathbf{v}}_j}{\partial x_i}\frac{\partial {\overline{\mathbf{w}}}_j}{\partial x_i}\right)\in \mathbb{R}.$$

$${\left|\left|\left|\mathbf{v}\right|\right|\right|}_X={\left|\left|\left|\mathbf{v}\right|\right|\right|}_{1,\Omega ,d}\text{and}{\left|\left|\left|q\right|\right|\right|}_Q={\left|\left|\left|q\right|\right|\right|}_{0,\Omega },$$

$$\mathcal{T}$$

$${\mathcal{T}}_{\perp }$$

$$\{\begin{array}{rl}a(\phi ,\psi )& ={\int }_{\Omega }\alpha \nabla \phi :\nabla \overline{\psi }dx+i{\int }_{\Sigma }(\phi \cdot \mathbf{n})(\overline{\psi }\cdot \mathbf{n})d\sigma ,\forall \phi ,\psi \in V_{\Gamma },\\ b(\phi ,\lambda )& =-{\int }_{\Omega }\bar{\lambda }(\nabla \cdot \phi )dx,\forall \phi \in V_{\Gamma },\forall \lambda \in Q,\\ F(\psi )& ={\int }_{\Omega }\mathbf{f}\cdot \overline{\psi }dx,\forall \psi \in V_{\Gamma }.\end{array}$$

$$\left\{\begin{array}[]{rcll}a({\bm{\mathbf{u}}},{\bm{\mathbf{\varphi}}})+b(\bm{\mathbf{\varphi}},p)&=&F(\bm{\mathbf{\varphi}}),&\quad\forall{\bm{\mathbf{\varphi}}}\in{V_{\Gamma}},\\ b(\bm{\mathbf{u}},\lambda)&=&0,&\quad\forall\lambda\in{Q}.\\ \end{array}\right.$$

$${\left|\left|\left|\mathbf{u}\right|\right|\right|}_X,{\left|\left|\left|p\right|\right|\right|}_Q\lesssim {\left|\left|\left|\mathbf{f}\right|\right|\right|}_{0,\Omega }.$$

$$\Re (a(\phi ,\phi ))⩾c_a{\left|\left|\left|\phi \right|\right|\right|}_X,\forall \phi \in V_{\Gamma },$$

$$\underset{\begin{array}{c}\lambda \in Q\\ \lambda \ne 0\end{array}}{\inf }\underset{\begin{array}{c}\phi \in V_{\Gamma }\\ \phi \ne 0\end{array}}{\sup }\frac{b(\phi ,\lambda )}{{\left|\left|\left|\phi \right|\right|\right|}_X{\left|\left|\left|\lambda \right|\right|\right|}_Q}⩾{\beta }_0.$$

$$|a(\phi ,\psi )|$$

$$|F(\phi )|$$

$$|b(\phi ,\lambda )|$$

$$\Re (a(\phi ,\phi ))=\alpha {\left|\left|\left|\nabla \phi \right|\right|\right|}_{0,\Omega }^2⩾c_a{\left|\left|\left|\phi \right|\right|\right|}_{1,\Omega }^2,$$

Full text

Keywords: coupled complex boundary method, free surface flow, shape optimization, shape derivatives, rearrangement method, and adjoint method

1. 1 Introduction

2. 2 CCBM in shape optimization settings

3. 2.1 The coupled complex boundary method formulation

4. 2.2 Well-posedness of the state problem

5. 2.3 The proposed shape optimization formulation

6. 3 Computation of the shape derivative

7. 3.1 Some concepts from shape calculus

8. 3.2 Some identities from trangential shape calculus

9. 3.3 Computation of the shape gradient

10. 4 Numerical approximation

11. 4.1 Numerical algorithm

12. 4.2 Numerical examples

13. 4.2.1 Example in two dimensions

14. 4.2.2 Example in three dimensions

15. 5 Conclusions and Future Work

16. A Appendices

17. A Lemmata proofs

18. A.1 Proof of Lemma 3.1

19. A.2 Proof of Lemma 3.2

20. B Computations of some identities

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

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

23. C Computation of the shape gradient of the cost via chain rule

1 Introduction

We are interested in a free boundary problem for fluid flows governed by the Stokes equations which arises in many applications such as in magnetic shaping processes where the shape of the fluid is determined by the Lorentz force. In such setting, the model is described by the Stokes flow equations and a pressure balance equation on the free boundary neglecting surface tension effects [BPST17]. Depending on the application, two model problems can be considered in this context. First, an interior case where the fluid is confined in a mould and has an internal unknown boundary and, second, an exterior case where a part of the fluid boundary adheres to a solid and the remaining part is free and is in contact with the ambient air. Here, we pay attention to the second case for $d$-dimensional geometries, $d\in\{2,3\}$.

Main Problem. In this work, we are particularly interested in the free surface problem analog of the Bernoulli free boundary problem or Alt-Caffarelli problem (see [FR97, AC81]) where the Laplace operator is replaced by the Stokes operator. More precisely, we are interested in the problem formulated as follows (see [Kas14] and also [BPST17]). Given a simply connected bounded $A\subset\mathbb{R}^{d}$ with boundary $\Gamma:=\partial A$ domain $\Omega$. The fluid is considered levitating around $A$, which is influenced by a gravity-like force $\bm{\mathbf{f}}=(f_{1},f_{2},\ldots,f_{d})^{\top}$, and occupies then the domain $\Omega=B\setminus\overline{A}$, where $B$ is a larger bounded (simply connected) domain with boundary $\Sigma:=\partial B$ containing $\overline{A}$. The incompressible viscous flow occupying $\Omega$, the velocity field $\bm{\mathbf{u}}=(u_{1},u_{2},\ldots,u_{d})^{\top}$, and the pressure $p$ is then supposed to satisfy the following overdetermined system of Stokes equations in non-dimensional form:

[TABLE]

where $\alpha:={Re}^{-1}$, and $Re>0$ is the Reynolds number, and $\partial_{\bm{\mathbf{n}}}{}:=\bm{\mathbf{n}}\cdot\nabla$ denotes the derivative in the (outward unit) normal vector $\bm{\mathbf{n}}$. Here, the boundary data $\bm{\mathbf{g}}$ is a prescribed velocity which in some sense expresses the motion of $\Gamma$ while the boundary condition imposed on the free surface $\Sigma$ indicate zero ambient pressure and negligible surface tension effects222The zero-surface tension assumption is a typical setup in the literature which not only simplifies the discussion, but also allows one to ignore technical difficulties resulting from higher derivative terms.. Additional assumptions concerning on the density of external forces $\bm{\mathbf{f}}$ and on $\bm{\mathbf{g}}$ will be given in next section. Meanwhile, the slip boundary condition $\bm{\mathbf{u}}\cdot\bm{\mathbf{n}}=0$ on $\Sigma$ means that inflow or outflow of the fluid cannot happen on the free boundary while tangential velocities on the said boundary can be non-zero. Such condition is appropriate for problems that involve free boundaries and situations where the usual no-slip condition $\bm{\mathbf{u}}=\bm{\mathbf{0}}$ is not valid (e.g., flows past chemically reacting walls see [BJ67, Lia01]). Other realistic situations where no-slip and slip boundary conditions appear naturally are as follows. The no-slip condition has long been firmly established for moderate pressures and velocities, not only by direct observations (see, e.g., [GR86] or [Gal94]), but also through comparisons between numerical simulations and results from actual experiments of a large class of complex flow problems. Early experimental investigations, however, showed that low-temperature slip past solid surfaces. This happens, in particular, for sufficiently large Knudsen numbers where velocity slip occurs at the wall surface. Such a physical phenomenon has also been observed in hydraulic fracturing and biological fluids [Lia01] which are examples of nonlinear fluid flows. A slip boundary condition also applies to free surfaces in free boundary problems. A well-known example where the condition arises is the so-called coating problem; see, e.g., [Bab63, SS81].

Clearly, motivated by numerous applications such as ship hydrodynamics [AD94] and thin film manufacturing processes [VPLG09], numerical solutions of flows that are partially bounded by a freely moving boundary are of great practical importance. In such problems, not only the flow variables, but also the flow domain has to be determined – hence the word ‘free boundary’ is used for the unknown boundary. Due to the complexity of simultaneously resolving these unknowns, a numerical solution which can be obtained through different methods [WSUR96] has to be carried out iteratively [KT99].

Known approaches. In this investigation, we shall resolve the free boundary problem via a novel shape optimization method. Because free surface problems (FSPs) such as (1) consist of overdetermined boundary conditions on the unknown part of the boundary, they can be reformulated into an equivalent shape optimization problem; see, e.g., [Kas14, BPST17].333The same technique – but in the context of optimal shape design problems – of finding the boundary that minimizes a norm of the residual of one of the free surface conditions, subject to the boundary value problem with the remaining free surface conditions imposed, has also been used for potential free surface flows in [VRK01, BS03]. Such approach can be carried out in several ways. A typical strategy is to choose one of the boundary conditions on the unknown boundary to obtain a well-posed state equation and then track the remaining boundary data in an appropriate norm – typically in $L^{2}$ (see, e.g., [EH09, EH10, HIK*+*09, HKKP03, IKP06] for the exterior case of the Bernoulli problem and [Kas14] for free surface problems). Another one is to utilize the so-called Kohn-Vogelius cost functional which consists of two auxiliary problems each posed with one of the boundary condition on the free boundary (see [ABP*+*13, Bac13, EH12] for the Bernoulli problem and for free surface problems [Kas14, BPST17]). For the latter formulation, one considers the minimization problem

[TABLE]

where the state variables $\bm{\mathbf{u}}_{N}:=\bm{\mathbf{u}}_{N}(\Omega)$ and $\bm{\mathbf{u}}_{D}=\bm{\mathbf{u}}_{D}(\Omega)$ respectively satisfy the following well-posed systems of partial differential equations (PDEs):

[TABLE]

[TABLE]

In (3), $\bm{\mathbf{\tau}}$ denotes the tangential vector to a domain $\Omega$. Meanwhile, for tracking-boundary-data cost functional approach, the following minimization problems can be considered:

[TABLE]

where, of course, $\bm{\mathbf{u}}_{\text{D}}$ and $\bm{\mathbf{u}}_{\text{N}}$ satisfy problems (3) and (4), respectively.

Remark 1.1**.**

If $(\Omega,(\bm{\mathbf{u}},p))$ is a solution of (1), then $\bm{\mathbf{u}}_{D}=\bm{\mathbf{u}}_{N}=\bm{\mathbf{u}}$. Hence, $J_{KV}(\Omega)=J_{D}(\Omega)=J_{N}(\Omega)=0$. Conversely, if $J_{i}(\Omega)=0$, $i\in\{KV,D,N\}$, are defined by (2), (5), and (6), respectively, then for $i=KV$ ($J_{KV}(\Omega)=0$), then we have that $\bm{\mathbf{u}}_{D}=\bm{\mathbf{u}}_{N}=\bm{\mathbf{u}}$ and $\bm{\mathbf{u}}_{D}=\bm{\mathbf{u}}_{N}=\bm{\mathbf{u}}$ is a solution of problem (1). For $i=D$ ($J_{D}(\Omega)=0$), we have that $\bm{\mathbf{u}}_{N}\cdot\bm{\mathbf{n}}=0$ almost everywhere on $\Sigma$ and $\bm{\mathbf{u}}=\bm{\mathbf{u}}_{N}$ is a solution of problem (1). Similarly, the same holds for $i=N$ ($J_{N}(\Omega)=0$). Therefore, $(\Omega,(\bm{\mathbf{u}},p))$ is a solution of (1) if and only if $J_{i}(\Omega)=0$, $i\in\{KV,D,N\}$.

Remark 1.2**.**

The identity $J_{KV}(\Omega)=0$ is equivalent to the existence of a constant $\lambda_{0}\in\mathbb{R}$ such that $(\bm{\mathbf{u}}_{D},p_{D})=(\bm{\mathbf{u}}_{D},p_{N}+\lambda_{0})$ [BPST17].

Remark 1.3**.**

In addition to the abovementioned classical approaches, there is another domain-integral-type penalization that may be considered. More exactly, one may opt to penalize, instead of $J_{KV}$, by the cost functional

[TABLE]

Compared with $J_{KV}$, however, the shape gradient of the above cost function is more complicated in structure and would be numerically expensive to evaluate. In fact, it can be shown that in this case one would need to solve five real BVPs: two state equations, two adjoint equations (see [LP12] for the case of the Bernoulli problem), and one to approximate the mean curvature of the free boundary implicitly.

In all of the minimization problems presented above, the infimum has always to be taken over all sufficiently smooth domains. On the other hand, we emphasize that the cost function $J_{N}$ requires more regularity for $p_{D}$ and $\bm{\mathbf{u}}_{\text{D}}$ to be well defined. As a consequence, it may be impractical to utilize this functional in numerical experiments where high regularity of the state variables is not guaranteed. For the feasibility of these approaches, with $\bm{\mathbf{g}}\equiv\bm{\mathbf{0}}$, we refer the readers to [Kas14]. In that paper, Kasumba computed explicitly the shape derivatives of the functionals $J_{KV}$, $J_{D}$, and $J_{N}$ via the chain rule approach. Then, he utilizes the gradient information in a boundary variation method in a preconditioned steepest descent algorithm to solve the given shape optimization problems.

New strategy and novelty. In this investigation, we want to offer yet another shape optimization technique to solve (1). More exactly, we want to propose here a novel application of the so-called coupled complex boundary method or CCBM for solving stationary free surface problems. The point of departure of the method is somewhat similar to [Tii97], but applies the concept of complex PDEs. The idea is straightforward: we couple the Dirichlet and Neumann data in a Robin boundary condition in such a way that the Dirichlet data and the Neumann data are the respective real and imaginary parts of the Robin boundary condition. In effect, the boundary conditions that have to be satisfied on the free boundary are transformed into one condition that needs to be satisfied on the domain. Consequently, the formulation requires the introduction of a new cost function that has not been studied yet in the literature (see (20)). Moreover, the new reformulation seems to be advantageous in the sense that it leads to a volume integral and produces more regular adjoint state than those of boundary-data tracking-type cost functionals.

The coupled complex boundary method was first introduced by Cheng et al. in [CGHZ14] as a method to solve an inverse source problem (see also [CGH16]) and was then used to deal with a Cauchy problem stably in [CGH16]. It is then applied to an inverse conductivity problem with one measurement in [GCH17] and also to parameter identification in elliptic problems in [ZCG20]. In a much recent paper, CCBM was also applied in solving inverse obstacle problems by Afraites in [Afr22]. Moreover, in the context of free boundary problems, the method was first applied by Rabago in [Rab22] to solve the exterior case of the Bernoulli problem. As far as we are concern, there is no work yet dealing with FSPs via CCBM, except in [OCNN22] where the authors applied CCBM as a numerical approach to solve an inverse Cauchy Stokes problem. Therefore, the novelty of this work comes from the application of CCBM in treating free surface problems, especially as a numerical resolution to the free boundary problem (1) to which we found some merit over the classical shape optimization methods of tracking the boundary data using a least-squares misfit cost function.

Contributions to the literature. The main contributions and highlights of this study are listed as follows.

• •

As opposed to classical shape optimization approaches for free surface problems [Kas14, BPST17], our state problem is a complex PDE system (see equation (7) in the next section), and to the best our knowledge, the Stokes equation with a Robin condition having complex Robin coefficient has not been considered yet in any of the previous investigations – at least in the context of shape optimization. The problem therefore is new and hence warrants some initial examination. Although the approach to establish the well-posedness of the problem follows the same argumentations as in the real case, the discussion of the topic (which we provide here) is not yet available in the existing literature. See subsection 2.2 for the discussion.

• •

As a consequence of the CCBM formulation, a new cost functional for the free surface problem – which apparently has not been examined yet in the literature – is first introduced in this work; see equation (20).

• •

The study focuses on the rigorous computation of the first-order shape derivative of the cost functional associated with CCBM (refer to Proposition 3.6) using only the Hölder continuity of the state variables (see Lemma 3.12) – differing from the classical approach which uses either the material or shape derivatives of the states (see, e.g., [Kas14, BPST17, Rab22]). In fact, our approach bypasses the computation of the said derivatives.

• •

The proof of the Hölder continuity of the state variables proceeds in a slightly different manner from the usual strategy found, for instance, in [BP13, IKP06, IKP08]; see the proof of Lemma 3.12. By applying the same technique used in the aforementioned papers, it seems difficult to eliminate the pressure variables to get a uniform estimate for the difference between the velocity solution of the transformed and the steady case of the Stokes problem. To get around this difficulty, we instead use the corresponding expansions of some specific expressions appearing in the computation. See the proof of Lemma 3.12.

• •

As far as we are aware of, previous numerical investigations on the free boundary problem with the Stokes equation (such as that of [Kas14]) only dealt with two dimensional cases. In this study, the proposed shape optimization method is also tested in solving a test problem in three dimensions. In fact, it will be shown in the numerical part of the paper that the proposed formulation has a smoothing effect – specially in the case of three dimensions – in approximating a solution to the overdetermined system of complex partial differential equations. See the discussion on subsection 4.2.2.

The remainder of the paper is organized as follows. In Section 2, we will demonstrate how the free surface problem (1) is formulated into a shape optimization problem via the coupled complex boundary method. The well-posedness of the CCBM formulation is also discussed in this section. Meanwhile, we devote Section 3 in computing the boundary integral expression of the shape gradient of the cost associated with the proposed shape optimization problem. The section of the paper starts with a brief discussion of some ideas from shape calculus needed in the study, followed by a short list of identities from tangential shape calculus. For the main result, which pertains to the shape gradient of the cost, the expression will be characterized rigorously via the rearrangement method in the spirit of [IKP06, IKP08]. In Section 4, the continuous formulation is discretized, and a numerical algorithm based on Sobolev gradient method for solving the discrete shape optimization problem is developed and implemented. The section therefore is divided into two parts: the discussion of the iterative scheme (subsection 4.1) and the presentation of the numerical experiments carried out in two and three dimensions (subsection 4.2). The main content of the paper ends with Section 5 where we issue a short conclusion about the study and a statement of future work.

Additionally, the paper also contains three appendices wherein we provide the proofs of some lemmas (Appendix A), show the computation of some identities used in the investigation (Appendix B), and demonstrate the derivation of the shape gradient expression via the chain rule approach (Appendix C).

2 CCBM in shape optimization settings

We present here the proposed CCBM formulation of (1) and discuss the well-posedness of the state problem.

2.1 The coupled complex boundary method formulation

The coupled complex boundary method suggests to write the boundary conditions on the free boundary as one condition. This means to consider the complex boundary value problem

[TABLE]

where $i=\sqrt{-1}$ stands for the unit imaginary number. Letting $(\bm{\mathbf{u}},p):=(\bm{\mathbf{u}}_{r}+i\bm{\mathbf{u}}_{i},p_{r}+ip_{i})$ denote the solution of (7), then it can be verified that the real vector-valued functions $\bm{\mathbf{u}}_{r}$ and $\bm{\mathbf{u}}_{i}$, and real scalar valued functions $p_{r}$ and $p_{i}$, respectively satisfy the real PDE systems:

[TABLE]

[TABLE]

From now on, if there is no confusion, we represent the real and imaginary parts of any function (vector-valued or scalar valued) by attaching the function with the subscript $\cdot\,_{r}$ and $\cdot\,_{i}$, respectively.

Remark 2.1**.**

Observe from the previous systems of real PDEs that if $\bm{\mathbf{u}}_{i}=\bm{\mathbf{0}}$ and $p_{i}=0$ in $\Omega$, then we have $\bm{\mathbf{u}}_{i}=\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}_{i}}=0$ and $p_{i}=0$ on $\Sigma$, and thus $\bm{\mathbf{u}}_{r}\cdot\bm{\mathbf{n}}=0$ on $\Sigma$. From (8) and (9), we see that the pair $(\Omega,\bm{\mathbf{u}}_{r})$ solves the original free boundary problem (1). Conversely, if $(\Omega,\bm{\mathbf{u}}_{r})$ is the solution to the free boundary problem (1), then clearly $\bm{\mathbf{u}}_{r}$ and $\bm{\mathbf{u}}_{i}$ satisfy (8) and (9).

We infer from the previous remark that the original free boundary problem (1) can be recast into an equivalent shape optimization problem given as follows.

Problem 2.2**.**

Given a fixed interior boundary $\Gamma$ and a function $\bm{\mathbf{f}}$, find an annular domain $\Omega$, with the exterior boundary denoted by $\Sigma:=\partial\Omega\setminus\Gamma$, and a function $\bm{\mathbf{u}}:=\bm{\mathbf{u}}(\Omega)$ such that $\bm{\mathbf{u}}_{i}=0$ in $\Omega$ and $\bm{\mathbf{u}}=\bm{\mathbf{u}}_{r}+i\bm{\mathbf{u}}_{i}$ solves the PDE system (7).

Unless otherwise specified, we assume throughout the paper that $\Omega$ is a bounded, non-empty, connected subset of $\mathbb{R}^{d}$, and of class (at least) $\mathcal{C}^{1,1}$.

2.2 Well-posedness of the state problem

We discuss here the well-posedness of the PDE system (7). To this end, for simplicity, we carry out the analysis on the basis of homogenous Dirichlet boundary conditions on fixed boundary $\Gamma$. Extension to non-homogeneous Dirichlet boundary conditions can be accomplished by standard techniques.

To proceed, we first introduce some notations.

Notations. Let us start with the notations for some operators and operations. We recall, for clarity, the following standard operators/operations

[TABLE]

where $\bm{\mathbf{n}}=(n_{1},\ldots,n_{d})^{\top}$ is the outward unit normal vector. Also, we denote the inner product444Another notation for the inner product with different usage will also be introduced shortly in the succeeding discussion. of two (column) vectors $\bm{\mathbf{a}},\bm{\mathbf{b}}\in\mathbb{R}^{d}$ in $\mathbb{R}^{d}$ by $\bm{\mathbf{a}}\cdot\bm{\mathbf{b}}:=\bm{\mathbf{a}}^{\top}\bm{\mathbf{b}}\equiv\langle\bm{\mathbf{a}},\bm{\mathbf{b}}\rangle_{\mathbb{R}^{d}}=\langle\bm{\mathbf{a}},\bm{\mathbf{b}}\rangle$, respectively (where the latter notation is used when there is no confusion).

For a vector-valued function $\bm{\mathbf{u}}:=(u_{1},u_{2},\ldots,u_{d})^{\top}:\Omega\to\mathbb{R}^{d}$, the gradient of $\bm{\mathbf{u}}$, denoted by $\nabla\bm{\mathbf{u}}$, a second-order tensor defined as $\nabla\bm{\mathbf{u}}=(\nabla\bm{\mathbf{u}})_{ij}:=(\partial u_{j}/\partial x_{i})_{i,j=1,\ldots,d}$, where $(\nabla\bm{\mathbf{u}})_{ij}$ is the entry at the $i$th row and $j$th column. Meanwhile, the Jacobian of $\bm{\mathbf{u}}$, denoted by $D\bm{\mathbf{u}}$ (the total derivative of $\bm{\mathbf{u}}$), is the transpose of the gradient (i.e., $D\bm{\mathbf{u}}=(D\bm{\mathbf{u}})_{ij}=(\partial u_{i}/\partial x_{j})_{i,j=1,\ldots,d}=\nabla^{\top}\bm{\mathbf{u}}$). Accordingly, we write the $\bm{\mathbf{n}}$-directional derivative of $\bm{\mathbf{u}}$ as $\partial_{\bm{\mathbf{n}}}\bm{\mathbf{u}}:=(D\bm{\mathbf{u}})\bm{\mathbf{n}}$. Other notations are standard and will only be emphasized for clarity.

The notations for the function spaces used in the paper are as follows. We denote by $W^{m,p}(\Omega)$ the standard real Sobolev space with the norm $\|\cdot\|_{W^{m,p}(\Omega)}$. Let $W^{0,p}(\Omega)=L^{p}(\Omega)$, and particularly, we write $H^{m}(\Omega)$ to represent the space

[TABLE]

where $D^{\beta}$ is the weak (or distributional) partial derivative operator and $\beta$ is a multi-index. Its corresponding inner product is given by $(\cdot,\cdot)_{m,\Omega}$ and norm

[TABLE]

For vector-valued functions, we define the Sobolev space $H^{m}(\Omega)^{d}$ as follows:

[TABLE]

Its associated norm is given by

[TABLE]

Similar definition is given when $\Omega$ is replaced by the boundary $\partial\Omega$.

We let $\bm{\mathbf{\mathsf{H}}}^{m}(\Omega)^{d}$ be the complex version of $H^{m}(\Omega)^{d}$ with the inner product $(\!(\cdot,\cdot)\!)_{m,\Omega,d}$ and norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\,\cdot\,\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{m,\Omega,d}$ defined respectively as follows:

[TABLE]

Also, for ease of writing, we will use the following notations in the paper:

[TABLE]

In above, the operation ‘$:$’ stands for the tensor scalar product and is defined as

[TABLE]

Because $\ensurestackMath{\stackon[0.5pt]{V}{\mkern 1.0mu\circ}}_{0}$ is a closed subspace of $\ensurestackMath{\stackon[0.5pt]{V}{\mkern 1.0mu\circ}}$, we have the decomposition $\ensurestackMath{\stackon[0.5pt]{V}{\mkern 1.0mu\circ}}=\ensurestackMath{\stackon[0.5pt]{V}{\mkern 1.0mu\circ}}_{0}\oplus\ensurestackMath{\stackon[0.5pt]{V}{\mkern 1.0mu\circ}}_{\perp}$. Also, with the above definitions, we sometimes write

[TABLE]

and drop $d$ when there is no confusion.

Lastly, throughout the paper, $c$ will denote a generic positive constant which may have a different value at different places. Also, we occasionaly use the symbol “$\lesssim$” which means that if $x\lesssim y$, then we can find some constant $c>0$ such that $x\leqslant cy$. Of course, $y\gtrsim x$ is defined as $x\lesssim y$. Other notations are standard and will only be emphasized for clarity.

Next we quote the following lemma which is essential in our argumentation.

Lemma 2.3**.**

Let $\Omega$ be a connected set and consider the following maps

[TABLE]

*Then, the following results holds: $\mathcal{T}$ is surjective, $\mathcal{T}_{\perp}\in\mathscr{L}(\ensurestackMath{\stackon[0.5pt]{V}{\mkern 1.0mu\circ}},\ensurestackMath{\stackon[0.5pt]{Q}{\mkern 1.0mu\circ}})$, $\mathcal{T}_{\perp}$ is bijective, and $\mathcal{T}_{\perp}\in\operatorname{Isom}(\ensurestackMath{\stackon[0.5pt]{V}{\mkern 1.0mu\circ}}_{\perp},\ensurestackMath{\stackon[0.5pt]{Q}{\mkern 1.0mu\circ}})$, i.e., there exists $\mathcal{T}_{\perp}^{-1}:\ensurestackMath{\stackon[0.5pt]{Q}{\mkern 1.0mu\circ}}\to\ensurestackMath{\stackon[0.5pt]{V}{\mkern 1.0mu\circ}}_{\perp}$, $\mathcal{T}_{\perp}^{-1}\in\mathscr{B}:=\mathscr{L}(\ensurestackMath{\stackon[0.5pt]{Q}{\mkern 1.0mu\circ}},\ensurestackMath{\stackon[0.5pt]{V}{\mkern 1.0mu\circ}}_{\perp})$. *

Proof.

The proof of the lemma follows similar arguments as in the real case (see, e.g.,[GR86, Chap. 1.2]), so we omit it. ∎

Now we exhibit the well-posedness of the complex PDE system (7). On this purpose, we introduce the following forms:

[TABLE]

Hence, we can state the weak formulation of (7) as follows: find $({\bm{\mathbf{u}}},p)\in V_{\Gamma}\times Q$ such that

[TABLE]

Our first proposition is given next.

Proposition 2.4**.**

For a given $\bm{\mathbf{f}}\in L^{2}(\Omega)^{d}$, the weak formulation of (7) given above admits a unique solution $(\bm{\mathbf{u}},p)\in{V_{\Gamma}}\times Q$ which depends continuously on the data. Moreover, we have

[TABLE]

The proof of the proposition is based on standard arguments as in the real case. Indeed, the validity of the statement follows from the continuity of the sesquilinear form $a(\cdot,\cdot)$ on $X\times X$, its coercivity in $X$, i.e., to show that there exists a constant $c_{a}>0$ such that

[TABLE]

the continuity of $F$ in $X$, and on an inf-sup condition. For the last condition, it needs to be shown that the bilinear form $b$ satisfies the condition that there is a constant $\beta_{0}>0$ such that

[TABLE]

Let us confirm the above results for completeness. Let $\bm{\mathbf{f}}\in L^{2}(\Omega)^{d}$ and $\bm{\mathbf{\varphi}},\bm{\mathbf{\psi}}\in V_{\Gamma}$, ${\lambda}\in{Q}$.

Continuity We have the following computations

[TABLE]

Coercivity We also have the inequality condition

[TABLE]

for some constant $c_{a}>0$, which confirms (14).

Inf-Sup condition Next, we want to show that we can find $\beta_{0}>0$ such that

[TABLE]

Let ${\lambda}\in Q$ be fixed arbitrarily. We define ${\lambda}=\ensurestackMath{\stackon[0.5pt]{\lambda}{\mkern 1.0mu\circ}}+{\lambda}_{\ast}$ with $\ensurestackMath{\stackon[0.5pt]{\lambda}{\mkern 1.0mu\circ}}\in\ensurestackMath{\stackon[0.5pt]{Q}{\mkern 1.0mu\circ}}$ and ${\lambda}_{\ast}:=\dfrac{1}{|\Omega|}\displaystyle\int_{\Omega}{{\lambda}}{\,{d}x}\in\mathbb{C}$. By Lemma 2.3, there exists $\ensurestackMath{\stackon[0.5pt]{\bm{\mathbf{\varphi}}}{\mkern 1.0mu\circ}}\in\ensurestackMath{\stackon[0.5pt]{V}{\mkern 1.0mu\circ}}_{\perp}$ such that $-\nabla\cdot\ensurestackMath{\stackon[0.5pt]{\bm{\mathbf{\varphi}}}{\mkern 1.0mu\circ}}=\ensurestackMath{\stackon[0.5pt]{\lambda}{\mkern 1.0mu\circ}}$ ($=\mathcal{T}_{\perp}\ensurestackMath{\stackon[0.5pt]{\bm{\mathbf{\varphi}}}{\mkern 1.0mu\circ}}$). We let $\widetilde{\bm{\mathbf{\varphi}}}_{\ast}$ be a fixed function in $V_{\Gamma}$ such that $\int_{\Sigma}{\widetilde{\bm{\mathbf{\varphi}}}_{\ast}\cdot\bm{\mathbf{n}}}{\,{d}\sigma}\neq 0$ and define

[TABLE]

With the above function, we consider $\bm{\mathbf{\varphi}}=\ensurestackMath{\stackon[0.5pt]{\bm{\mathbf{\varphi}}}{\mkern 1.0mu\circ}}+t_{0}{\lambda}_{\ast}\widetilde{\bm{\mathbf{\varphi}}}_{\ast}\in V_{\Gamma}$ where $t_{0}>0$. This leads to the following inequality

[TABLE]

For the numerator $N_{1}$, we have the following computations

[TABLE]

Here, $\langle 1\rangle$ is obtained from the fact that $\ensurestackMath{\stackon[0.5pt]{\bm{\mathbf{\varphi}}}{\mkern 1.0mu\circ}}\in\ensurestackMath{\stackon[0.5pt]{V}{\mkern 1.0mu\circ}}_{0}$ and the assumption that $-\nabla\cdot\ensurestackMath{\stackon[0.5pt]{\bm{\mathbf{\varphi}}}{\mkern 1.0mu\circ}}=\ensurestackMath{\stackon[0.5pt]{\lambda}{\mkern 1.0mu\circ}}$, $\langle 2\rangle$ follows from Green’s theorem, together with (17), $\langle 3\rangle$ is due to the assumption that $\widetilde{\bm{\mathbf{\varphi}}}_{\ast}$ is fixed, so ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widetilde{\bm{\mathbf{\varphi}}}_{\ast}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{1,\Omega}$ equates to some constant $c_{0}>0$, while inequality $\langle 4\rangle$ is obtained because of Peter-Paul inequality555Here, of course, $\varepsilon_{0}>0$ is chosen such that $c_{1}>0$. For example, this condition holds if we choose $\varepsilon_{0}=c_{0}^{-1}>0$ and $t_{0}=\varepsilon c_{0}^{-1}=c_{0}^{-2}>0$. applied to the product $|{\lambda}_{\ast}|{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\ensurestackMath{\stackon[0.5pt]{\lambda}{\mkern 1.0mu\circ}}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{Q}$.

Let us further estimate below the sum $h({\lambda})$. Note that we have the following calculations

[TABLE]

Thus, there actually exist constants $c_{2},c_{3}>0$ such that the following inequalities hold

[TABLE]

By this estimate, we get

[TABLE]

for the numerator in (18).

For the denominator $N_{2}$ in (18), we have the following estimations

[TABLE]

In above, we have used the fact that $\mathcal{T}_{\perp}^{-1}$ is a bounded linear operator. Meanwhile, $\left\|\cdot\right\|_{\mathscr{B}}$ denotes the operator norm for linear operators in $\mathscr{B}=\mathscr{L}(\ensurestackMath{\stackon[0.5pt]{Q}{\mkern 1.0mu\circ}},\ensurestackMath{\stackon[0.5pt]{V}{\mkern 1.0mu\circ}}_{\perp})$. Combining the above estimate with (19), we finally get

[TABLE]

This proves (16), and thus (15).

2.3 The proposed shape optimization formulation

To solve Problem 2.2, we introduce the cost functional

[TABLE]

where $(\bm{\mathbf{u}}_{i},p_{i})$ is subject to the state problem (7). Notice that compared to $J_{KV}$, the cost function $J$ only requires the solution of a single complex PDE problem to be solved. The optimization problem we consider here is the problem of minimizing $J(\Omega)$ over a set of admissible domains $\mathcal{O}_{\text{ad}}$, where $\mathcal{O}_{\text{ad}}$ is essentially the set of $\mathcal{C}^{k,1}$, $k\geqslant 1$, $k\in\mathbb{N}$, (non-empty) doubly connected domains with (fixed) interior boundary $\Gamma$ and (free) exterior boundary $\Sigma$. In other words, we consider the shape optimization problem that reads as follows: find $\Omega$ such that

[TABLE]

We note that it is actually enough to consider $\Gamma$ to be only Lipschitz regular to derive the shape derivative of the cost function, but for simplicity we also assume it to be $\mathcal{C}^{k,1}$ regular. We emphasize that, in this paper, we will not tackle the question of existence of optimal shape solution to (21). Instead, we will tacitly assume the existence of solution to the free surface problem (1) and adopt the shape optimization formulation (21) to resolve (1) numerically.

To numerically solve the optimization problem $J(\Omega)\to\inf$, we will apply a shape-gradient-based descent method based on finite element method (FEM). We will not, however, employ any kind of adaptive mesh refinement in our numerical scheme as opposed to [RA19a, RA19b, RA20]. In this way, we can further assess the stability of the new method in comparison with the classical least-squares method of tracking the Dirichlet data. The expression for the shape derivative of the cost will be exhibited in the next section using shape calculus [DZ11, HP18, MS76, Sim80, SZ92]. However, as opposed to [Kas14, BPST17], our strategy to obtain the shape gradient expression does not require the knowledge of the shape derivative of the state, nor its material derivative.

Throughout the rest of the paper, we shall refer to our proposed shape optimization method simply by CCBM.

3 Computation of the shape derivative

The main purpose of this section is to characterize the first-order shape derivative of the cost function $J$ given in Proposition 3.6.

3.1 Some concepts from shape calculus

Before we derive the shape gradient of the cost function $J$, we briefly review in this section some important concepts and results from shape calculus.

Consider a (convex) bounded hold-all domain $U$ of class $\mathcal{C}^{k,1}$, $k\geqslant 1$, strictly containing $\overline{\Omega}$ (the closure of $\Omega$) and define $T_{t}$ as the perturbation of the identity $id$666Here, $id$ is also used to denote the identity matrix in $d$-dimension. If there is no confusion, this abuse of notation is used throughout the paper. given by the map

[TABLE]

where $\bm{\mathbf{\theta}}:=(\theta_{1},\theta_{2},\ldots,\theta_{d})^{\top}\in\mathbb{R}^{d}$ is a $t$-independent deformation field belonging to the admissible space

[TABLE]

and $k\in\mathbb{N}$, later on specified as what is needed. Hereinafter, $t$ is assumed sufficiently small such that $T_{t}$ is a diffeomorphism from $\Omega\in\mathcal{C}^{k,1}$ onto its image.888In fact, $T_{t}:U\longmapsto U$ is a $\mathcal{C}^{k,1}$ diffeomorphism for sufficiently small $t$. In other words, we suppose that the reference domain $\Omega$ and its perturbation $\Omega_{t}$ have the same topological structure and regularity under the transformation $T_{t}$. On this purpose, particularly for $k=1$ (thus, $\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}^{1}$), we let $\varepsilon\in\mathbb{R}^{+}$ be sufficiently small such that $[t\mapsto T_{t}]\in\mathcal{C}^{1}(\mathcal{I},\mathcal{C}^{1,1}(\overline{U})^{d})$, where $\mathcal{I}:=[0,\varepsilon]$. Throughout the paper, when not specified, it is always assumed that (at least) $k=1$ so that $T_{t}$ becomes a $\mathcal{C}^{1,1}$ diffeomorphism.

By definition, we also have the perturbation $\Sigma_{t}:=T_{t}(\Sigma)$ and $\Gamma_{t}:=T_{t}(\Gamma)\equiv\Gamma$ where the latter identity is due to the fact that $\bm{\mathbf{\theta}}=\bm{\mathbf{0}}$ on $\Gamma$. In addition, of course, $\Omega_{0}=\Omega$ and $\Sigma_{0}=\Sigma$. Accordingly, the set of all admissible domains $\mathcal{O}_{ad}$ is given as follows

[TABLE]

The functional $J:\mathcal{O}_{ad}\to\mathbb{R}$ has a directional first-order Eulerian derivative at $\Omega$ in the direction of the field $\bm{\mathbf{\theta}}$ if the limit

[TABLE]

exists (see, e.g., [DZ11, Sec. 4.3.2, Eq. (3.6), p. 172]). The functional $J$ is said to be shape differentiable at $\Omega$ in the direction of $\bm{\mathbf{\theta}}$ if the map $\bm{\mathbf{\theta}}\mapsto{d}J(\Omega)[\bm{\mathbf{\theta}}]$ is linear and continuous. In this case, we refer to the map as the shape gradient of $J$.

We introduce a few more notations for ease of writing. We denote by $DT_{t}$ the Jacobian matrix of $T_{t}$ (i.e., $(DT_{t})_{ij}=\partial(T_{t})_{i}/\partial x_{j}$) and write the inverse and inverse transpose of this matrix by $(DT_{t})^{-1}$ and $(DT_{t})^{-\top}:=((DT_{t})^{\top})^{-1}$, respectively. For convenience, we define

[TABLE]

For $t\in\mathcal{I}$, $I_{t}$ is positive. At $t=0$, it is evident that $I_{0}=1$, $A_{0}=id$ (because $DT_{0}^{-1}=id$ and $(DT_{0})^{-\top}=id$), and $B_{0}=1$. Moreover, on $\mathcal{I}$, the maps $t\mapsto I_{t}$, $t\mapsto A_{t}$, and $t\mapsto B_{t}$ are continuously differentiable. That is, for $t\in\mathcal{I}$, we have (see, e.g., [IKP06, IKP08])

[TABLE]

Furthermore, $[t\mapsto I_{t}]\in\mathcal{C}^{1}(\mathcal{I},\mathcal{C}^{0,1}(\overline{U}))$ and $[t\mapsto A_{t}]\in\mathcal{C}(\mathcal{I},\mathcal{C}(\overline{U})^{d\times d})$.

The derivatives of the maps in (26) are respectively given as follows:

[TABLE]

where ${\operatorname{div}}_{\Sigma}\bm{\mathbf{\theta}}$ denotes the tangential divergence of the vector $\bm{\mathbf{\theta}}$ on $\Sigma$. Meanwhile, for simplicity we write $\theta_{n}$ instead of $\bm{\mathbf{\theta}}\cdot\bm{\mathbf{n}}$.

Additionally to the above results, we specifically point out that (d/dt)DT_{t}\big{\rvert}_{t=0}=D\bm{\mathbf{\theta}}\in\mathcal{C}^{0,1}(\overline{U})^{d\times d} which implies that $[t\mapsto DT_{t}]\in\mathcal{C}^{1}(\mathcal{I},\mathcal{C}^{0,1}(\overline{U})^{d\times d})$. On the other hand, we have $[t\mapsto T_{t}^{-1}]\in\mathcal{C}(\mathcal{I},\mathcal{C}^{1}(\overline{U})^{d})$ from which we see that $[t\mapsto(DT_{t})^{-\top}]\in\mathcal{C}^{1}(\mathcal{I},\mathcal{C}(\overline{U})^{d\times d})$.

We also assume that for $t\in\mathcal{I}$, we have

[TABLE]

for all $\xi\in\mathbb{R}^{d}$, for some constants $\Lambda_{1}$, $\Lambda_{2}$, $\Lambda_{3}$, and $\Lambda_{4}$ ($\Lambda_{1}<\Lambda_{2}$, $\Lambda_{3}<\Lambda_{4}$).

Furthermore, as we intend to refer any function $\bm{\mathbf{\varphi}}_{t}:\Omega_{t}\to\mathbb{R}^{d}$ to the reference domain using the transformation $T_{t}$, we shall make use of the notation $\bm{\mathbf{\varphi}}^{t}$ which is essentially the composite function $\bm{\mathbf{\varphi}}_{t}\circ T_{t}:\Omega\to\mathbb{R}^{d}$.

Lastly, to close the section we list some properties of the Jacobian of the map $T_{t}$ and state an auxiliary result that will be helpful with our investigation. Their proofs are provided in Appendix A.

Lemma 3.1**.**

Let $\Omega$ and $U$ be nonempty bounded open connected subsets of $\mathbb{R}^{d}$ (defined as before). Also, let $T_{t}$ be given by (22) and $\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}^{k}$, $k\geqslant 1$. Then, we have the following results:

• •

$DT_{t}$* and $DT_{t}^{-1}$ are Lipschitz continuous in $U$.*

• •

there exists some constant $C>0$ such that $\lvert DT_{t}^{-1}(x)\rvert_{\infty}<C<\infty$, for all $x\in\overline{U}$.

• •

$\lim_{t\to 0}\dfrac{1}{t}\left(DT_{t}-I\right)=D\bm{\mathbf{\theta}}$* and $\lim_{t\to 0}\dfrac{1}{t}\left[(DT_{t})^{-1}-I\right]=-D\bm{\mathbf{\theta}}$.*

Lemma 3.2**.**

Let $\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}^{k}$, $k\geqslant 1$. Then, we have the following

1. (i)

for any $f\in L^{p}(U)$, $p\geqslant 2$, we have

[TABLE] 2. (ii)

for any $f\in W^{1,p}(U)$, $p\geqslant 2$, we have

[TABLE] 3. (iii)

for $f\in W^{2,p}(U)$, $p\geqslant 1$, the mapping $t\mapsto f\circ T_{t}$ from $\mathcal{I}\to W^{1,p}(\Omega)$ is differentiable at $t=0$ and the derivative is given by

[TABLE] 4. (iv)

for $f\in H^{1}(U)$. Then, the map $t\to I_{t}f\circ T_{t}$ from $\mathcal{I}$ to $L^{2}(\Omega)$ is differentiable at $t=0$ and the derivative is given by

[TABLE]

For $f\in L^{p}(U)$, $p\geqslant 1$, the limit $\lim_{t\to 0}f\circ T_{t}=f$ in $L^{p}(\Omega)$, where $\overline{\Omega}\subset U$, also holds; see [IKP06, Lem. 3.4].

Here we note that similar results to the previous lemma for the case of vector-valued functions also hold. For instance, for $\bm{\mathbf{\varphi}}=(\varphi_{1},\varphi_{2},\ldots,\varphi_{d})^{\top}\in H^{2}(U)^{d}$, we have

[TABLE]

Moreover, the map $t\mapsto\bm{\mathbf{\varphi}}\circ T_{t}$ from $\mathcal{I}\to H^{1}(\Omega)^{d}$ is differentiable at $t=0$ and we have

[TABLE]

Furthermore, the mapping $t\to I_{t}\bm{\mathbf{\varphi}}\circ T_{t}$ from $\mathcal{I}$ to $L^{2}(\Omega)^{d}$ is differentiable at $t=0$ and the following limit holds

[TABLE]

where

[TABLE]

and $\otimes$ denotes the outer product (i.e., $\bm{\mathbf{\varphi}}\otimes\bm{\mathbf{\theta}}$ gives us the matrix $\bm{\mathbf{\varphi}}\bm{\mathbf{\theta}}^{\top}=(\varphi_{j}\theta_{k})_{jk}$ where $j,k=1,2,\ldots,d$).

The proof of the above results follows the same argumentation in validating Lemma 3.2 which is provided in Appendix A.2.

3.2 Some identities from trangential shape calculus

As we will see in the next section, the expression for the shape gradient of $J$ contains expressions from tangential shape calculus (see, e.g., [DZ11, Chap. 9, Sec. 5] and see also [HP18, Sec. 5.4.3, pp. 216–221] for a detailed discussion). So, to facilitate our investigation, we summarize below some notations that will appear in our subsequent discussions.

Definition 3.3** ([DZ11, Chap. 9, Sec. 5.2, eq. (5.17) – (5.19), p. 497]).**

For any domain $\Omega\subset\mathbb{R}^{d}$ with $\mathcal{C}^{1,1}$ smooth boundary $\Gamma:=\partial\Omega$, the following expressions are well-defined:

• •

the tangential gradient of $\psi\in\mathcal{C}^{1}(\Gamma)$ is given by

[TABLE]

• •

the tangential Jacobian matrix of a vector function $\bm{\mathbf{\varphi}}\in\mathcal{C}^{1}(\Gamma;\mathbb{R}^{d})$ is given by

[TABLE]

• •

the tangential divergence of a vector function $\bm{\mathbf{\varphi}}\in\mathcal{C}^{1}(\Gamma;\mathbb{R}^{d})$ is given by

[TABLE]

In above definitions, $\tilde{\psi}$ and $\tilde{\bm{\mathbf{\varphi}}}$ are any corresponding $\mathcal{C}^{1}$ extensions of $\psi$ and $\bm{\mathbf{\varphi}}$ into a neighborhood of $\Gamma$.

With the above definitions, the formulas given in the next lemma can easily be shown.

Lemma 3.4** ([DZ11, Chap. 9, Sec. 5.5, eq. (5.26) – (5.27), p. 498]).**

Consider a $\mathcal{C}^{1,1}$ domain $\Omega$ with boundary $\Gamma:=\partial\Omega$. Then, for $\psi\in H^{1}(\Gamma)$ and $\bm{\mathbf{\varphi}}\in\mathcal{C}^{1}(\Gamma;\mathbb{R}^{d})$ the following identities hold:

• •

tangential divergence formula (see also **[SZ92, Lem. 2.63, eq. (2.140), p. 91]**):

[TABLE]

• •

tangential Stokes’ formula:

[TABLE]

• •

tangential Green’s formula:

[TABLE]

As a consequence of the previous identity, we have

[TABLE]

Another version of the tangential Green’s formula (whose proof can be found, for instance, in [MS76]) is given in the next lemma (see, e.g., [SZ92, Chap. 2, eq. (2.144), p. 92]).

Lemma 3.5**.**

Let $U$ be a bounded domain of class $\mathcal{C}^{1,1}$ and $\Omega\subset U$ with boundary $\Gamma:=\partial\Omega$. Consider the functions $\bm{\mathbf{\varphi}}\in\mathcal{C}^{1,1}(\overline{U};\mathbb{R}^{d})$ and $\psi\in W^{2,1}(U)$. Then, the following equation holds

[TABLE]

For a smooth enough function $\bm{\mathbf{\varphi}}$, we denote its second normal derivative as $\partial_{\bm{\mathbf{n}}\bm{\mathbf{n}}}^{2}{\bm{\mathbf{\varphi}}}$. Let us denote by $\mathcal{M}_{d\times d}$ the space of the matrix of size $d\times d$. The tangential differential operators will be represented with the subscript $\cdot\,_{\Sigma}$. Particularly, for $\bm{\mathbf{\varphi}}\in\bm{\mathbf{W}}^{1,1}(\partial\Omega)$ and $\bm{\mathbf{M}}\in W^{1,1}(\partial\Omega;\mathcal{M}_{d\times d})$, the following operators are defined on $\Gamma:=\partial\Omega$:

• •

tangential vectorial gradient operator:

[TABLE]

• •

tangential vectorial divergence operator [DZ11, Chap. 9, Sec. 5.2, eq. (5.11), p. 496]:

[TABLE]

• •

Laplace-Beltrami operator [DZ11, Chap. 9, Sec. 5.3, eq. (5.12), p. 496]:

[TABLE]

In relation to the last definition, the Laplace-Beltrami operator can be decomposed in the following way (see, e.g., [DZ11, p. 28] or [HP18, Prop. 5.4.12, eq. (5.59), p. 220])

[TABLE]

In the next section, we will exhibit the computation of the shape gradient of $J$.

3.3 Computation of the shape gradient

We compute the shape derivative of $J$ rigorously via rearrangement method in the spirit of [IKP06] given that $\bm{\mathbf{f}}\in H^{1}(U)^{d}$, the $\Omega$ of class $\mathcal{C}^{1,1}$, and $\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}^{1}$. This approach has the advantage of bypassing the need to calculate the material derivative of the state. Moreover, instead of the $\mathcal{C}^{1}$ continuity of the map $t\mapsto u^{t}$ around [math], we only need the Hölder continuity of the state $u^{t}$ – the composition $u_{t}\circ T_{t}$ which is defined on $\Omega$. Our main result is given in the following proposition.

Proposition 3.6**.**

Let $\bm{\mathbf{f}}\in H^{1}(U)^{d}$, $\Omega$ be of class $\mathcal{C}^{1,1}$, and $\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}^{1}$. Then, $J$ is shape differentiable and its shape derivative is given by ${d}J(\Omega)[\bm{\mathbf{\theta}}]=\int_{\Sigma}{{g_{\Sigma}}\bm{\mathbf{n}}\cdot\bm{\mathbf{\theta}}}{\,{d}\sigma}$ where101010Here $\nabla_{\Gamma}$ is the tangential gradient operator on $\Gamma$. The intrinsic definition of the operator is given, for instance, in [DZ11, Chap. 5., Sec. 5.1, p. 492].

[TABLE]

In above expression of the shape gradient of $J$, $\bm{\mathbf{B}}[\theta_{n}]$ is given by

[TABLE]

and $\kappa$ denotes the mean curvature of the free boundary $\Sigma$, $\bm{\mathbf{u}}=\bm{\mathbf{u}}_{r}+i\bm{\mathbf{u}}_{i}$ and $p=p_{r}+ip_{i}$ is the unique pair of solution to (7), and the pair of adjoints $(\bm{\mathbf{v}},q)$, where $\bm{\mathbf{v}}=\bm{\mathbf{v}}_{r}+i\bm{\mathbf{v}}_{i}$ and $q=q_{r}+iq_{i}$, uniquely solves the adjoint system

[TABLE]

The proof of the above result relies on several lemma that we first prove below. To start, we present the weak formulation of (39). On this purpose, we introduce the following forms:

[TABLE]

Hence, we may write the weak formulation of (39) as follows: find $({\bm{\mathbf{v}}},q)\in V_{\Gamma}\times Q$ such that

[TABLE]

Remark 3.7**.**

The well-posedness of the above problem can be verified using similar arguments issued in subsection 2.2 for the state problem (7) (provided, of course, (13) admits a unique solution $({\bm{\mathbf{u}}},p)\in V_{\Gamma}\times Q$). Particularly, one can show that an inf-sup condition (cf. (15)) holds for the above variational problem.

Remark 3.8**.**

For $\Omega$ of class $\mathcal{C}^{k+1,1}$, $k$ a non-negative integer, it can be shown that the weak solution $(\bm{\mathbf{u}},p)\in X\times Q$ to the variational problem (13) is also $\bm{\mathbf{\mathsf{H}}}^{k+2}(\Omega)^{d}\times H^{k+1}(\Omega)$. In particular, $\bm{\mathbf{u}}_{i}\in H^{k+2}(\Omega)^{d}$ and $p_{i}\in H^{k+1}(\Omega)$. Consequently, we find that the weak solution $(\bm{\mathbf{v}},q)$ of problem (43) is not only in $V_{\Gamma}\times Q$, but is also an element of $\bm{\mathbf{\mathsf{H}}}^{k+4}(\Omega)^{d}\times H^{k+3}(\Omega)$.

At this juncture, we introduce the following sesquilinear forms $a^{t},\tilde{a}^{t}\in V_{\Gamma}\times V_{\Gamma}\to\mathbb{R}$ and linear forms, $b^{t}:V_{\Gamma}\times Q\to\mathbb{R}$ and $F^{t},\tilde{F}^{t}:V_{\Gamma}\to\mathbb{R}$ (which are essentially the transformed versions of the forms listed in (10) and (40)) defined as follows:

[TABLE]

The first of the several lemmas that we need is given as follows.

Lemma 3.9**.**

For $t\in\mathcal{I}$, the sesquilinear forms $a^{t}$ and $\tilde{a}^{t}$ defined on $X\times X$ are bounded and coercive on $X\times X$. In addition, the linear forms $F^{t}(\bm{\mathbf{\varphi}})$ and $\tilde{F}^{t}(\bm{\mathbf{\varphi}})$ are also bounded. Moreover, the bilinear form $b$ satisfies the condition that there is a constant $\beta_{1}>0$ such that

[TABLE]

Proof.

The proof is similar to the arguments used in the discussion issued in subsection 2.2. This time, however, one has to take into account the properties of $A_{t}$ and $B_{t}$ given in (26) and the bounds given in (28) to prove the given results. ∎

The next lemma is concerned about the well-posedness of the transported perturbed version of (13).

Lemma 3.10**.**

The pair $(\bm{\mathbf{u}}^{t},p^{t})=(\bm{\mathbf{u}}_{r}^{t}+i\bm{\mathbf{u}}_{i}^{t},p_{r}^{t}+ip_{i}^{t})$ uniquely solves in $V_{\Gamma}\times Q$ the variational problem

[TABLE]

Proof.

The functions ${\bm{\mathbf{u}}_{t}}\in\bm{\mathbf{\mathsf{H}}}^{1}(\Omega_{t})^{d}$ and ${p_{t}}\in\bm{\mathbf{\mathsf{L}}}^{2}(\Omega_{t})$ solve the variational equation

[TABLE]

By change of variables (cf. [DZ11, pp. 482–484]), together with the identities

[TABLE]

and the transformation [DZ11, Thm. 4.4, p. 488]

[TABLE]

we get

[TABLE]

In light of the notations listed in (44), we recover (48).

The rest of the proof is similar to the arguments used in subsection 2.2 combined with the properties of $A_{t}$ and $B_{t}$ issued in (26) and (28), and together with Lemma 3.9. Concerning uniqueness of solution, the proof is also standard, so we omit it. This proves the lemma. ∎

Lemma 3.11**.**

For $t\in\mathcal{I}$, the solutions $(\bm{\mathbf{u}}^{t},p^{t})$ of (48) are uniformly bounded in $X\times Q$. More precisely, for all $t\in\mathcal{I}$, there exists a constant $c>0$ independent of $t$ such that

[TABLE]

Proof.

For the (uniform) boundedness of $\bm{\mathbf{u}}^{t}$ in $X$ for $t\in\mathcal{I}$, the result is established by taking $\bm{\mathbf{\varphi}}=\bm{\mathbf{u}}^{t}\in V_{\Gamma}$ in (48), and then applying the properties of $A_{t}$ and $B_{t}$ given in (26) and (28), and the coercivity of $a^{t}$ on $X\times X$ (Lemma 3.9). Indeed, from the first equation in (48) together with the second equation $b^{t}(\bm{\mathbf{u}}^{t},\lambda)=0$ with $\lambda=p^{t}\in Q$, we have

[TABLE]

This yields the following inequality

[TABLE]

where the constant $c>0$ appearing on the leftmost side of the inequality is the coercivity constant of $a^{t}$ whose existence is ensured by Lemma 3.9. To further get an estimate for the rightmost side of the above inequality, we note of the following calculation

[TABLE]

where in the last inequality we used the bound for $I_{t}$ given in (28). Hence, we have

[TABLE]

This gives us a priori estimate

[TABLE]

for some constant $c>0$, which shows that $\bm{\mathbf{u}}^{t}$ is bounded in $X$ for $t\in\mathcal{I}$.

For the boundedness of $p^{t}$ in $Q$ for $t\in\mathcal{I}$, we make use of (45), which is equivalent to the following

[TABLE]

We let $\lambda=p^{t}\in Q$. Then, from (3.10), the bounds in (28), and equation (49), we have the following calculation

[TABLE]

Therefore, using the estimate for ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{\mathbf{u}}^{t}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{X}$ given in (50), we arrive at

[TABLE]

for some constant $c>0$. The desired inequalities then follow from the fact that $\bm{\mathbf{f}}\in H^{1}(U)^{d}$. This proves the lemma. ∎

To complete our preparations, we next prove the following lemma concerning the Hölder continuity of the state variables $\bm{\mathbf{u}}^{t}$ and $p^{t}$ with respect to $t$.

Lemma 3.12**.**

Let $(\bm{\mathbf{u}},p)\in V_{\Gamma}\times Q$ be the solution of (13). Then, the following limit holds

[TABLE]

where $({\bm{\mathbf{u}}^{t}},p^{t})\in V_{\Gamma}\times Q$ solves (48), for $t\in\mathcal{I}$.

Proof.

Let us first note that the assumption that $\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}^{1}$ implies that

[TABLE]

which implies further that $|\bm{\mathbf{\theta}}|_{\infty}$ and $|D\bm{\mathbf{\theta}}|_{\infty}$ are both finite.

Hereafter, we proceed in two steps:

Step 1.

We first show that $\lim_{t\to 0}\bm{\mathbf{u}}^{t}=\bm{\mathbf{u}}$ in $X$ and $\lim_{t\to 0}p^{t}=p$ in $Q$.

Step 2.

Then, we validate our claim that $\lim_{t\to 0}\dfrac{1}{\sqrt{t}}\left({\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{X}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|p^{t}-p\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{Q}\right)=0$.

Step 1. We consider the difference between the variational equations (48) and (13) and define $\bm{\mathbf{y}}^{t}:=\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}}$ and $r^{t}:=p^{t}-p$. By making $\varepsilon>0$ smaller if necessary, it can be shown that the following expansions hold for sufficiently small $t\in\mathcal{I}$,

[TABLE]

where $\lambda_{max}$ is the dominant eigenvalue (i.e., the eigenvalue with the largest modulus) of $D\bm{\mathbf{\theta}}$. Below we use the above expansions and specifically write $M_{t}^{\top}$ as follows

[TABLE]

Henceforth, we denote

[TABLE]

Now, let us consider the variational equation

[TABLE]

Applying (51) to the above equation, we get

[TABLE]

Taking $\lambda=r^{t}=p^{t}-p\in Q$, we get

[TABLE]

On the other hand, we also have the equation

[TABLE]

which is equivalent to

[TABLE]

where

[TABLE]

and $a:V_{\Gamma}\times V_{\Gamma}\to\mathbb{R}$ and $b:V_{\Gamma}\times Q\to\mathbb{R}$ are respectively the sesquilinear and linear forms given in (10).

By choosing $\bm{\mathbf{\varphi}}=\bm{\mathbf{y}}^{t}\in V_{\Gamma}$ in (54) and utilizing identity (53), it follows that

[TABLE]

Therefore, we have the estimate

[TABLE]

where

[TABLE]

Now, let us take $\lambda=r^{t}=p^{t}-p\in Q$ in the inf-sup condition (16) and consider equation (54). We have

[TABLE]

Therefore, we have the inequality

[TABLE]

Going back (56) and utilizing the above estimate, we get

[TABLE]

We apply Peter-Paul inequality to $\left({m}_{t}+t{\Xi}^{t}\beta_{0}^{-1}\max\{|\alpha|,1\}\right){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{\mathbf{y}}^{t}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{X}$ to obtain the following estimate

[TABLE]

for some constant $\varepsilon_{1}>0$. We choose (and fixed) $\varepsilon_{1}$ such that $\bar{c}:=c(c_{a},\varepsilon_{1}):=2c_{a}-\varepsilon_{1}>0$ so that, from our first estimate (60), we have

[TABLE]

We note that, at $t=0$, $I_{t}|M_{t}\bm{\mathbf{n}}|^{-1}=(1)(|id\bm{\mathbf{n}}|^{-1})=1$. Thus, in view of (57) together with Lemma 3.11, (29) and the vector-valued version of Lemma 3.2 (refer to (31)), we see that ${m}_{t}\to 0$ as $t\to 0$. Moreover, it is not hard to see that ${\Xi}^{t}$ is (uniformly) bounded for all $t\in\mathcal{I}$ because of Lemma 3.11. Therefore, we can conclude – by Lebesgue’s dominated convergence theorem – that

[TABLE]

Similarly, based from the above discussion and from (58), we know that the terms on the right-hand side of the said inequality vanish as $t\to 0$. Thus, we also have the limit

[TABLE]

Step 2. Before we proceed to the last part of the proof, we note that, for sufficiently small $t>0$, $\frac{1}{t}\bm{\mathbf{y}}^{t}\in V_{\Gamma}$ and $\frac{1}{t}r^{t}\in Q$. In addition, we recall that the derivatives of $I_{t}$ and $M_{t}$ with respect to $t$ exists in $L^{\infty}(\Omega)$ and $L^{\infty}(\Omega)^{d\times d}$, respectively.

Now, to finish the proof, we go back to the computations in the previous step (referring particularly to (61)), to obtain, after dividing by $t>0$,

[TABLE]

Observe that, we have

[TABLE]

Thus, using (26), (27), (31), Lemma 3.11, (62), and (63), we deduce that the following limit holds

[TABLE]

Similarly, from estimate (59), we know that (after dividing by $\sqrt{t}>0$) the following inequality holds

[TABLE]

Again, applying (26), (27), (31), and Lemma 3.11, but now combined with (64), we infer that

[TABLE]

In above argumentations of taking the limit (equations (64) and (65)), it is understood that we are taking a sequence $\{t_{n}\}$ such that $\lim_{n\to\infty}t_{n}=0$ and we have exploited the fact that the following sequences

[TABLE]

are (uniformly) bounded in $L^{\infty}$ (because of (26), (27), and (28)). In addition, Lemma 3.11 was employed above in the manner that, for the sequence $\{t_{n}\}$ (passing to a subsequence if necessary), the solutions $(\bm{\mathbf{u}}^{t_{n}},p^{t_{n}})$ of (48) are uniformly bounded in $X\times Q$. Note that these results allow us to interchange the limit and the supremum in obtaining, particularly, equation (65). Finally, combining (64) and (65) concludes the lemma. ∎

We are now in the position to prove Proposition 3.6 without using the shape derivative of $\bm{\mathbf{u}}$ and $p$. Before we dealt with the proof, we make a few preparations by proving the following four lemmas which primarily contain several identities. These results will serve useful in deriving the boundary expression for the shape gradient in accordance with Hadamard-Zolésio structure theorem [DZ11, Chap. 9, Sec. 3.4, Thm. 3.6, p. 479].

Lemma 3.13**.**

The following equation is satisfied by the state solution $\bm{\mathbf{u}}\in X$ of (13):

[TABLE]

Proof.

Let $(\bm{\mathbf{u}}^{t},p^{t})\in V_{\Gamma}\times Q$ be as in Lemma 3.10. We recall equation (52) and divide both sides by $t>0$, to obtain the following equation

[TABLE]

Using the properties of $\rho(t)$ and $R(t)$, we get – by letting $t$ goes to zero in (66) and in view of the last paragraph in the proof of the previous lemma – the following limit

[TABLE]

Because $\operatorname{div}\bm{\mathbf{u}}=0$ in $\Omega$, we arrive at the desired equation, ending the proof of the lemma. ∎

We next prove three identities which prove to be useful in simplifying the expression for the shape gradient.

Lemma 3.14**.**

Let $\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}^{1}$ and $\bm{\mathbf{\varphi}}$ be sufficiently smooth vector such that $\operatorname{div}\bm{\mathbf{\varphi}}=0$ in $\Omega\in\mathcal{C}^{1,1}$. Then, we have the following identity

[TABLE]

Proof.

To obtain the identity, we recall the definition of some of the tangential operators in Definition 3.3 and consider some identities (on $\Sigma$) given as follows:

[TABLE]

We have the following computations:

[TABLE]

This proves the lemma. ∎

Note that, by using tangential divergence formula, we can further write (67) into

[TABLE]

Lemma 3.15**.**

Let $\bm{\mathbf{\theta}}=(\theta_{1},\ldots,\theta_{d})^{\top}\in\bm{\mathbf{\Theta}}^{1}$ and $\bm{\mathbf{\varphi}}=(\varphi_{1},\ldots,\varphi_{d})^{\top}$ be a sufficiently smooth vector in $\Omega$. Then, the following identity holds

[TABLE]

Proof.

Let us first note that

[TABLE]

Therefore, direct computation of the divergence of $\nabla\bm{\mathbf{\varphi}}^{\top}\bm{\mathbf{\theta}}$ gives us

[TABLE]

proving the identity. ∎

Lemma 3.16**.**

Let $\bm{\mathbf{\theta}}=(\theta_{1},\ldots,\theta_{d})^{\top}\in\bm{\mathbf{\Theta}}^{1}$ and $\bm{\mathbf{\varphi}}=(\varphi_{1},\ldots,\varphi_{d})^{\top}$ be a given vector. Then, the following equation holds

[TABLE]

Proof.

We first note of the following identities:

[TABLE]

Therefore, we have the following computations

[TABLE]

which proves the equation. ∎

We now provide the proof of Proposition 3.6.

Proof of Proposition 3.6.

The proof essentially proceeds in two parts. First, we evaluate the limit $\lim_{t\to 0}\frac{1}{t}\left(J(\Omega_{t})-J(\Omega)\right)$. Then, using the regularity of the domain as well as the state and adjoint variables (to be introduced below), we characterized the boundary integral expression for the computed limit.

We begin by applying the following domain transformation formula (see, e.g., [SZ92, p. 77] or [DZ11, Chap. 9, Sec. 4.1, eq. (4.2), p. 482]):

[TABLE]

which holds for $\varphi_{t}\in L^{1}(\Omega_{t})$ combined with the identity $\eta^{2}-\zeta^{2}=(\eta-\zeta)^{2}+\zeta(\eta-\zeta)$ to obtain the following sequence of calculations:

[TABLE]

In view of (27)1 and using Lemma 3.12, we infer that

[TABLE]

Moreover, (27)1 also reveals that

[TABLE]

At this point, we put into use the expansion of the divergence of the product of a scalar function $\psi$ and a vector field $\bm{\mathbf{\varphi}}$ given by

[TABLE]

and then apply Green’s theorem to obtain the following equivalent form of the previous computed expression:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

and

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

• •

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

[TABLE]

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

[TABLE]

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:

[TABLE]

• •

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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)):

[TABLE]

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

[TABLE]

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

[TABLE]

or equivalently,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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)

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 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 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)

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.

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:

[TABLE]

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

[TABLE]

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

[TABLE]

By the triangle inequality, we know that

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

[TABLE]

Now, considering the expansion (compare (51))

[TABLE]

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:

[TABLE]

This implies that

[TABLE]

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]):

[TABLE]

(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

[TABLE]

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,

[TABLE]

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

[TABLE]

or equivalently,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

from which it can be deduced that

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

Moreover, we have the following computations

[TABLE]

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

[TABLE]

Thus, on $\Sigma$, we get

[TABLE]

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]):

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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)

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

to obtain

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

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$:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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))

[TABLE]

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))

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

as desired. ∎