A Shape Optimization Problem Constrained with the Stokes Equations to Address Maximization of Vortices

john.simon@stu.kanazawa-u.ac.jp notsu@se.kanazawa-u.se.ac

Abstract.

We study an optimization problem that aims to determine the shape of an obstacle that is submerged in a fluid governed by the Stokes equations. The mentioned flow takes place in a channel, which motivated the imposition of a Poiseuille-like input function on one end and a do-nothing boundary condition on the other. The maximization of the vorticity is addressed by the $L^{2}$ -norm of the curl and the det-grad measure of the fluid. We impose a Tikhonov regularization in the form of a perimeter functional and a volume constraint to address the possibility of topological change. Having been able to establish the existence of an optimal shape, the first order necessary condition was formulated by utilizing the so-called rearrangement method. Finally, numerical examples are presented by utilizing a finite element method on the governing states, and a gradient descent method for the deformation of the domain. On the said gradient descent method, we use two approaches to address the volume constraint: one is by utilizing the augmented Lagrangian method; and the other one is by utilizing a class of divergence-free deformation fields.

Key words and phrases:

Shape optimization, Stokes equations, augmented Lagrangian, Rearrangment method.

1991 Mathematics Subject Classification:

Primary: 49Q10,76D55, 35Q93; Secondary: 76U99

This work is supported by JSPS KAKENHI Grant Numbers JP18H01135, JP20H01823, JP20KK0058 and JP21H04431, and JST CREST Grant Number JPMJCR2014 for HN; and by the Japanese Government (MEXT) Scholarship for JSS

^∗ Corresponding author: john.simon@stu.kanazawa-u.ac.jp

John Sebastian H. Simon^∗

Graduate School of Natural Science and Technology

Kanazawa University

Kanazawa 920-1192, Japan

Hirofumi Notsu

Faculty of Mathematics and Physics

Kanazawa University

Kanazawa 920-1192, Japan

Introduction

The study of fluid dynamics is among the most active areas in mathematics, physics, engineering, and most recently in information theory [13]. An interesting area in this field is the study of turbulent flows and on controlling the emergence of such phenomena [7, 30, 10].

Even though turbulence represents chaos, harnessing it in a controlled manner can be useful, for instance, in optimal mixing problems. Among these studies is [8] where an optimal stirring strategy to maximize mixing is studied, while in [24] the authors studied the best type of input functions that will provide a better mixing in the fluid.

Recently, Goto, K. et.al.[13] utilized the generation of vortices in the fluid for information processing tasks. In the said literature, the authors found out that the length of the twin-vortex affects the memory capacity in the context of physical reservoir computing.

For these reasons, the current paper is dedicated to studying a shape optimization problem intended to maximize the turbulence of a flow governed by the two-dimensional Stokes equations. Here, the shape of an obstacle submerged in a fluid – that flows through a channel – is determined to maximize vorticity. Furthermore, the vorticity is quantified in two ways: one is by the curl of the velocity of the fluid; and the other is by rewarding the unfolding of complex eigenvalues of the gradient tensor of the velocity field.

To be precise, our goal is to study the following shape optimization problem

\min\{\mathcal{G}({\operatorname{\boldsymbol{u}}},\Omega);|\Omega|=m,\Omega\subset D\},

(1)

where $D$ is a hold-all domain which is assumed to be a fixed bounded connected domain in $\mathbb{R}^{2}$ , $\Omega\subset D$ is an open bounded domain, $\mathcal{G}({\operatorname{\boldsymbol{u}}},\Omega):=J({\operatorname{\boldsymbol{u}}},\Omega)+\alpha P(\Omega),$ $J$ and $P$ are vortex and perimeter functionals given by

J({\operatorname{\boldsymbol{u}}},\Omega)=-\int_{\Omega}\frac{\gamma_{1}}{2}|\nabla\times{\operatorname{\boldsymbol{u}}}|^{2}+\gamma_{2}h(\mathrm{det}(\nabla{\operatorname{\boldsymbol{u}}}))\operatorname*{d\!}x,\ P(\Omega)=\int_{\Gamma_{\rm f}}\operatorname*{d\!}s,

$h$ is such that

h(t)=\left\{\begin{matrix}0&\text{if }t\leq 0,\\ t^{3}/(t^{2}+1)&\text{if }t>0,\end{matrix}\right.

$|\Omega|:=\int_{\Omega}\operatorname*{d\!}x$ , and $\alpha,\gamma_{1},\gamma_{2}$ and $m$ are positive constants. Here, ${\operatorname{\boldsymbol{u}}}$ is the velocity field governed by a fluid flowing through a channel with an obstacle (see Figure 1 for reference). The flow of the fluid is reinforced by a divergence-free input function $\operatorname{\boldsymbol{g}}$ on the left end of the channel, which is an inflow boundary denoted by $\Gamma_{\rm{in}}$ , whilst an outflow boundary condition is imposed to the fluid on the right end of the channel denoted by $\Gamma_{\rm{out}}$ . The boundary $\Gamma_{\rm f}$ is the free surface and is the boundary of the submerged obstacle in the fluid. The remaining boundaries of the channel are wall boundaries denoted by $\Gamma_{\rm{w}}$ , upon which – together with the obstacle boundary $\Gamma_{\rm f}$ – a no-slip boundary condition is imposed on the fluid. We refer the reader to [11, 15, 21, 25, 33] for shape optimization problems involving fluids among others.

Figure 1. Set up of the domain.

We also point out that usually, the curl $\nabla\times{\operatorname{\boldsymbol{u}}}$ is enough to quantify the vorticity of the fluid. However, as pointed out by Kasumba, K. and Kunisch, K. in [23], the magnitude of such quantity may still be high on laminar flows. As such, on the same literature the authors proposed the second term in the vortex functional for the quantification of rotation of fluids. The impetus is that the rotational cores of fluids mostly occur near regions where the eigenvalues of $\nabla{\operatorname{\boldsymbol{u}}}$ are complex [4, 19, 22], and that the eigenvalues are complex when $\mathrm{det}\nabla{\operatorname{\boldsymbol{u}}}>0$ .

Convention: When we talk about a domain $\Omega$ , we always consider it having the boundary $\partial\Omega=\overline{\Gamma}_{\rm in}\cup\overline{\Gamma}_{\rm w}\cup\overline{\Gamma}_{\rm out}\cup\overline{\Gamma}_{\rm f}$ with $\mathrm{meas}(\Gamma_{i}\cap\Gamma_{j})=0$ for $i,j\in\{{\rm in,w,out,f}\},i\neq j$ , and that $\mathrm{dist}(\partial\Omega\backslash\Gamma_{\rm f},\Gamma_{\rm f})>0$ .

One of the challenges in this problem is the non-convexity of the functional $J$ , which leads into the possible non-existence of solutions for the shape optimization problem. Fortunately, we will be able to show later on that the functional $J$ is continuous with respect to domain variations. Furthermore, the Tikhonov regularization in the form of the perimeter functional $P$ is taken into consideration to circumvent the issue of possible topological changes of the domain.

Meanwhile, for the equality constraint $|\Omega|=m$ , we note that the use of an augmented Lagrangian generates smoother solutions as compared to the common Lagrangian. This is due to the regularizing effect of the quadratic augmentation. As a further consequence, this quadratic term also acts as a penalizing term for the violation of the equality constraint. Aside from the augmented Lagrangian, we also propose an analogue of the $H^{1}$ -gradient method that satisfies the incompressibility condition to preserve the volume on each iterate of the gradient-descent method.

This paper is organized as follows: in the next section, the system that governs the fluid and the functional spaces to be used for the analysis will be introduced. The variational formulation of the governing system and the existence of its solution will also be shown in the next section. Section 2 is dedicated to the analysis of the existence of an optimal shape. Here, we employ the $L^{\infty}$ -topology of characteristic functions for convergence of deformed domains, to ensure the volume preservation.

We shall derive the necessary condition for the optimization problem in Section 3. This is done by investigating the sensitivity of the objective functional $\mathcal{G}$ with respect to shape variations. In particular, the shape variations are induced by the velocity method [34], and the shape sensitivity is analyzed using the rearrangement method which was formalized by Ito, K., Kunisch, K., and Peichl, G. in [20]. In Section 4, we provide numerical algorithms through the help of the necessary conditions formulated in Section 3, and provide some illustrations of how these algorithms are implemented. To objectively observe the effects of the final shapes to fluid vortices, we look at the generation of twin-vortices in a dynamic nonlinear fluid flow. Finally, we draw conclusions and possible problems on the last section.

1. Preliminaries

1.1. Governing equations and necessary functional spaces

Let $\Omega\subset D$ be open, bounded and of class $C^{1,1}$ , the motion of the fluid is described by its velocity ${\operatorname{\boldsymbol{u}}}$ and pressure $p$ which satisfies the stationary Stokes equations given by:

\displaystyle\left\{\begin{aligned} \,-\nu\Delta{\operatorname{\boldsymbol{u}}}+\nabla p&=\operatorname{\boldsymbol{f}}&&\text{ in }\Omega,\\ \,\operatorname*{div}{\operatorname{\boldsymbol{u}}}&=0&&\text{ in }\Omega,\\ \,{\operatorname{\boldsymbol{u}}}&=\operatorname{\boldsymbol{g}}&&\text{ on }\Gamma_{\rm in},\\ \,{\operatorname{\boldsymbol{u}}}&=0&&\text{ on }\Gamma_{\rm{f}}\cup\Gamma_{\rm w},\\ \,-p\operatorname{\boldsymbol{n}}+\nu\partial_{\operatorname{\boldsymbol{n}}\!}{\operatorname{\boldsymbol{u}}}&=0&&\text{ on }\Gamma_{\rm out},\end{aligned}\right.

(2)

where $\nu>0$ is a constant, $\operatorname{\boldsymbol{f}}$ is an external force, $\operatorname{\boldsymbol{g}}$ is a divergence-free input function acting on the boundary $\Gamma_{\rm in}$ , $\operatorname{\boldsymbol{n}}$ is the outward unit normal vector on $\partial\Omega$ , and $\partial_{\operatorname{\boldsymbol{n}}}:=\operatorname{\boldsymbol{n}}\cdot\nabla$ is the $\operatorname{\boldsymbol{n}}$ -directional derivative on $\partial\Omega$ . The condition ${\operatorname{\boldsymbol{u}}}=0$ on $\Gamma_{\rm{f}}\cup\Gamma_{\rm w}$ corresponds to no-slip boundary condition, while the Neumann condition on $\Gamma_{\rm out}$ is the so-called do-nothing boundary condition.

The analysis will take place in the Sobolev spaces $W^{k,p}(D)$ for $D\subset\mathbb{R}^{2}$ , $k\geq 0$ and $p\geq 1$ . Note that if $p=2$ , then $H^{k}(D)=W^{k,2}(D)$ and $H^{0}(D)=L^{2}(D)$ . For any domain $\Omega\subset\mathbb{R}^{2}$ , we define

\mathcal{W}(\Omega)=\{\boldsymbol{\varphi}\in C^{\infty}(\Omega)^{2}:\boldsymbol{\varphi}=0\text{ on a neighborhood of }\Gamma_{0}:=\partial\Omega\backslash\Gamma_{\rm out}\}.

We denote by $\operatorname{\mathbf{H}}_{\Gamma_{0}}^{r}(\Omega)$ the closure of $\mathcal{W}(\Omega)$ with respect to the space $H^{r}(\Omega)^{2}$ .

We also consider the following bilinear forms $a(\cdot,\cdot)_{\Omega}:H^{2}(\Omega)^{2}\times H^{2}(\Omega)^{2}\to\mathbb{R}$ and $b(\cdot,\cdot)_{\Omega}:H^{2}(\Omega)^{2}\times L^{2}(\Omega)\to\mathbb{R}$ defined by

\displaystyle a({\operatorname{\boldsymbol{u}}},\operatorname{\boldsymbol{v}})_{\Omega}=\int_{\Omega}\nabla{\operatorname{\boldsymbol{u}}}:\nabla\operatorname{\boldsymbol{v}}\operatorname*{d\!}x,\ b(\operatorname{\boldsymbol{v}},q)_{\Omega}=-\int_{\Omega}q\operatorname{\mathrm{div}}\operatorname{\boldsymbol{v}}\operatorname*{d\!}x.

Furthermore, for any measurable set $\mathcal{D}\subset\mathbb{R}^{d}$ ( $d=1,2$ ) we shall denote $(\cdot,\cdot)_{\mathcal{D}}$ as the $L^{2}(\mathcal{D})$ , $L^{2}(\mathcal{D})^{2}$ or $L^{2}(\mathcal{D})^{2\times 2}$ inner product.

1.2. Weak formulation and existence of solutions

Let $\operatorname{\boldsymbol{g}}\in H^{2}(D)^{2}$ satisfy the following properties:

\displaystyle\left\{\begin{aligned} &\operatorname*{div}\operatorname{\boldsymbol{g}}=0\text{ in }\Omega;\quad\operatorname{\boldsymbol{g}}=0\text{ on }\Gamma_{\rm f}\cup\Gamma_{\rm wall};\\ &\int_{\Gamma_{\rm out}}\operatorname{\boldsymbol{g}}\cdot\operatorname{\boldsymbol{n}}\operatorname*{d\!}s=-\int_{\Gamma_{\rm in}}\operatorname{\boldsymbol{g}}\cdot\operatorname{\boldsymbol{n}}\operatorname*{d\!}s\geq 0,\end{aligned}\right.

(3)

and $\operatorname{\boldsymbol{f}}\in L^{2}(D)^{2}$ . The existence of the function $\operatorname{\boldsymbol{g}}\in H^{2}(D)^{2}$ that satisfies (3), can be easily established by utilizing for instance [12, Lemma 2.2].

By letting $\tilde{\operatorname{\boldsymbol{u}}}={\operatorname{\boldsymbol{u}}}-\operatorname{\boldsymbol{g}}$ , we consider the variational form of the state equation (2) given by: For a given domain $\Omega$ , find $(\tilde{\operatorname{\boldsymbol{u}}},p)\in\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(\Omega)\times L^{2}(\Omega)$ that satisfies, for any $(\operatorname{\boldsymbol{\varphi}},q)\in\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(\Omega)\times L^{2}(\Omega)$ , the following equations

\left\{\begin{aligned} \nu a(\tilde{\operatorname{\boldsymbol{u}}},\operatorname{\boldsymbol{\varphi}})_{\Omega}+b(\operatorname{\boldsymbol{\varphi}},p)_{\Omega}&=(\operatorname{\boldsymbol{f}},\operatorname{\boldsymbol{\varphi}})_{\Omega}-\nu(\nabla\operatorname{\boldsymbol{g}},\nabla\operatorname{\boldsymbol{\varphi}})_{\Omega}\\ b(\tilde{\operatorname{\boldsymbol{u}}},q)_{\Omega}&=0.\end{aligned}\right.

(4)

Any pair $(\tilde{\operatorname{\boldsymbol{u}}},p)\in\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(\Omega)\times L^{2}(\Omega)$ that solves the variational equation (4) is said to be a weak solution to the Stokes equation (2). The existence of the solution $\tilde{\operatorname{\boldsymbol{u}}}$ , is summarized below.

Theorem 1.1.

Let $\Omega$ be of class $C^{1,1},$ $\operatorname{\boldsymbol{f}}\in L^{2}(D)^{2}$ , and $\operatorname{\boldsymbol{g}}\in H^{2}(D)^{2}$ satisfy (3). The solution $(\tilde{\operatorname{\boldsymbol{u}}},p)\in\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(\Omega)\times L^{2}(\Omega)$ to the variational problem (4) exists such that

\|\tilde{\operatorname{\boldsymbol{u}}}\|_{\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(\Omega)}+\|p\|_{L^{2}(\Omega)}\leq c(\|\operatorname{\boldsymbol{f}}\|_{L^{2}(\Omega)^{2}}+\|\operatorname{\boldsymbol{g}}\|_{H^{1}(\Omega)^{2}}).

(5)

for some constant $c>0$ .

The proof of the theorem above can be done by utilizing the fact that the operator $a(\cdot,\cdot)_{\Omega}$ is coercive, $b(\cdot,\cdot)_{\Omega}$ satisfies the inf-sup condition [12], and that the right hand side of the first equation of (4) can be seen as an action of an element of $\operatorname{\mathbf{H}}^{-1}(\Omega)$ .

Remark 1.

(i) The energy estimate can be extended to the hold-all domain $D$ , i.e.,

\|\tilde{\operatorname{\boldsymbol{u}}}\|_{\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(\Omega)}+\|p\|_{L^{2}(\Omega)}\leq c(\|\operatorname{\boldsymbol{f}}\|_{L^{2}(D)^{2}}+\|\operatorname{\boldsymbol{g}}\|_{H^{1}(D)^{2}}),

where $c>0$ is dependent on $D$ but not on $\Omega$ . (ii) As expected, a more regular domain yields a more regular solution. In particular, if $\Omega$ is of class $C^{k,1}$ for $k\geq 0$ then the weak solution to (4) satisfies ${\operatorname{\boldsymbol{u}}}\in\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(\Omega)\cap\,H^{k+1}(\Omega)^{2}$ , and as a consequence of the Rellich-Kondrachov embedding theorem [9, Chapter 5, Theorem 6] the solution is in $C(\overline{\Omega})^{2}$ . Furthermore, if we instead assume that the outer boundary covers a convex polygonal domain while the inner boundary $\Gamma_{\rm f}$ is of class $C^{1,1}$ , the same regularity of the solution is obtained when $k=1$ .

2. Existence of optimal shapes

2.1. Cone property and the set of admissible domains

We begin by assuming, without loss of generality, that $\Gamma_{\rm out}\subset\partial D$ . We define the set of admissible domains as a subset of the collection of domains inside the hold-all domain $D$ that possesses the cone property. Furthermore, we define the topology on the said admissible set by the convergence of the corresponding characteristic function of each domain in the $L^{\infty}$ -topology. As highlighted by Henrot and Privat in [18] and is rigorously discussed in [6] and [17], this approach helps in the preservation of volume which is among the goals in our exposition. Several authors utilized parametrizing the free-boundary (c.f. [31, 32] and the references therein) to define the set of admissible domains, however free-boundary parametrization might lead to generating domains with varying volumes.

We shall adapt the definition of the cone property as in [3]. In what follows, $(\cdot,\cdot)$ and $\|\cdot\|$ denote the inner product and the norm in $\mathbb{R}^{2}$ , respectively.

Definition 2.1.

Let $h>0$ , $2h>r>0$ , $\theta\in[0,2\pi]$ , and $\xi\in\mathbb{R}^{2}$ such that $\|\xi\|=1$ .

(i) The cone of angle $\theta$ , height $h$ and axis $\xi$ is defined by

C(\xi,\theta,h)=\{x\in\mathbb{R}^{2};(x,\xi)>\|x\|\cos\theta,\|x\|<h\}.

(ii) A set $\Omega\subset\mathbb{R}^{2}$ is said to satisfy the cone property if and only if for all $x\in\partial\Omega$ , there exists $C_{x}=C(\xi_{x},\theta,h)$ , such that for all $y\in B(x,r)\cap\Omega$ we have $y+C_{x}\subset\Omega$ .

From this definition, the set of admissible domains as

\operatorname{\mathcal{O}_{ad}}:=\{\Omega\subset D;\Omega\text{ satisfies the }cone\ property\text{ and }|\Omega|=m\}.

This set of admissible domains has been established to be non-empty (see the proof of Proposition 4.1.1 in [17]), which exempts us from the futility of the analyses we will be discussing.

A sequence $\{\Omega_{n}\}_{n}\subset\operatorname{\mathcal{O}_{ad}}$ is said to converge to $\Omega\in\operatorname{\mathcal{O}_{ad}}$ if

\chi_{\Omega_{n}}\operatorname{\mathrel{\ensurestackMath{\stackon[1pt]{\rightharpoonup}{\scriptstyle\ast}}}}\chi_{\Omega}\text{ in }L^{\infty}(D),

where the function $\chi_{A}$ for a set $A\subset\mathbb{R}^{2}$ refers to the characteristic function defined by

\chi_{A}(x)=\left\{\begin{matrix}1&\text{if }x\in A,\\ 0&\text{if }x\not\in A.\end{matrix}\right.

Remark 2.

As explained by Henrot, A., et.al., [17, Proposition 2.2.1], the weak^∗ convergence in $L^{\infty}$ implies that the convergence also holds in the space $L^{p}_{loc}(\mathbb{R}^{2})$ and thus $\chi_{\Omega}$ is almost everywhere a characteristic function.

We shall denote the collection of characteristic functions of elements of $\operatorname{\mathcal{O}_{ad}}$ as $\mathcal{U}_{ad}$ , i.e., $\mathcal{U}_{ad}=\{\chi_{\Omega};\Omega\in\operatorname{\mathcal{O}_{ad}}\}$ , and whenever we mention $\mathcal{U}_{ad}$ we take into account the weak^∗ topology in $L^{\infty}(D)$ . We refer the reader to [6, Chapter 5] and [17, Chapter 2 Section 3] for a more detailed discussion on the topology of characteristic functions of finite measurable domains.

The compactness of $\operatorname{\mathcal{O}_{ad}}$ follows from the fact that it is closed and relatively compact – as defined by Chenais, D. [3] – with respect to the topology on $\mathcal{U}_{ad}$ . One can also read upon the proof in [17, Proposition 2.4.10]. We shall not discuss the proof of such properties, nevertheless they are summarized on the lemma below.

Lemma 2.2.

The set $\operatorname{\mathcal{O}_{ad}}$ is compact with respect to the topology on $\mathcal{U}_{ad}$ .

Another important implication of the cone property is the existence of a uniform extension operator.

Lemma 2.3 (cf. [3]).

Let $d=\{1,2\}$ and $m\in\mathbb{N}\cup\{0\}$ , there exists $K>0$ such that for all $\Omega\in\operatorname{\mathcal{O}_{ad}}$ , there exists

\mathcal{E}^{d}_{\Omega}:H^{m}(\Omega)^{d}\to H^{m}(D)^{d}

which is linear and continuous such that $\displaystyle\max\{\|\mathcal{E}^{d}_{\Omega}\|\}_{d=1,2}\leq K$ .

These result will be instrumental for proving some vital properties, for example when we show that the domain-to-state map is continuous.

Remark 3.

The set of admissible domains can be identified as a collection of Lipschitzian domains. Since we only consider the hold-all domain $D$ to possess a bounded boundary we refer the reader to [17, Theorem 2.4.7] for the proof of such property.

2.2. Well-posedness of the Optimization Problem

Before showing that the optimal shape indeed exists, we take note that from Theorem 1.1 the following map is well-posed:

\Omega\in\operatorname{\mathcal{O}_{ad}}\mapsto(\tilde{\operatorname{\boldsymbol{u}}}(\Omega),p(\Omega))\in\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(\Omega)\times L^{2}(\Omega).

This implies that we can write the objective functional in the manner that it solely depends on the domain $\Omega$ , i.e.,

\mathcal{J}(\Omega):=\mathcal{G}({\operatorname{\boldsymbol{u}}}(\Omega),\Omega),

where ${\operatorname{\boldsymbol{u}}}(\Omega)=\tilde{{\operatorname{\boldsymbol{u}}}}+\operatorname{\boldsymbol{g}}.$ However, the mentioned well-posedness of the domain-to-state map is insufficient to prove the existence of the minimizing domain. In fact, we need the continuity of the map $\Omega\mapsto(\tilde{\operatorname{\boldsymbol{u}}}(\Omega),p(\Omega))$ .

Proposition 1.

Let $\{\Omega_{n}\}_{n}\subset\operatorname{\mathcal{O}_{ad}}$ be a sequence that converges to $\Omega\in\operatorname{\mathcal{O}_{ad}}$ . Suppose that for each $\Omega_{n}$ , $(\tilde{\operatorname{\boldsymbol{u}}}_{n},p_{n})\in\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(\Omega_{n})\times L^{2}(\Omega_{n})$ is the weak solution of the Stokes equations on the respective domain; then the extensions $(\overline{{\operatorname{\boldsymbol{u}}}}_{n},\overline{p}_{n}):=(\mathcal{E}^{2}_{\Omega_{n}}\tilde{\operatorname{\boldsymbol{u}}}_{n},\mathcal{E}^{1}_{\Omega_{n}}p_{n})\in\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(D)\times L^{2}(D)$ coverges to a state $(\overline{\operatorname{\boldsymbol{u}}},\overline{p})\in\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(D)\times L^{2}(D)$ , such that $(\tilde{\operatorname{\boldsymbol{u}}},p)=(\overline{\operatorname{\boldsymbol{u}}},\overline{p})\big{|}_{\Omega}\in\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(\Omega)\times L^{2}(\Omega)$ is a solution to (4).

Proof.

From the uniform extension property, there exist $K>0$ such that

\|\overline{\operatorname{\boldsymbol{u}}}_{n}\|_{\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(D)}+\|\overline{p}\|_{L^{2}(D)}\leq K(\|\tilde{\operatorname{\boldsymbol{u}}}_{n}\|_{\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(\Omega_{n})}+\|p_{n}\|_{L^{2}(\Omega_{n})})\text{ for all }\Omega_{n}.

Furthermore, from Theorem 1.1

\|\tilde{\operatorname{\boldsymbol{u}}}_{n}\|_{\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(\Omega_{n})}+\|p_{n}\|_{L^{2}(\Omega_{n})}\leq c(\|\operatorname{\boldsymbol{f}}\|_{L^{2}(D)^{2}}+\|\operatorname{\boldsymbol{g}}\|_{H^{1}(D)^{2}}).

This implies uniform boundedness of $\{(\overline{\operatorname{\boldsymbol{u}}}_{n},\overline{p}_{n})\}_{n}$ in $\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(D)\times L^{2}(D)$ . Hence, by Rellich–Kondrachov’s theorem, there exists a subsequence of $\{(\overline{\operatorname{\boldsymbol{u}}}_{n},\overline{p}_{n})\}_{n}$ , which shall also be denoted as $\{(\overline{\operatorname{\boldsymbol{u}}}_{n},\overline{p}_{n})\}_{n}$ , and an element $(\overline{\operatorname{\boldsymbol{u}}},\overline{p})\in\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(D)\times L^{2}(D)$ such that the following properties hold:

\displaystyle\left\{\begin{aligned} &\overline{\operatorname{\boldsymbol{u}}}_{n}\rightharpoonup\overline{\operatorname{\boldsymbol{u}}}&&\text{in }H^{1}(D)^{2},\\ &\overline{\operatorname{\boldsymbol{u}}}_{n}\to\overline{\operatorname{\boldsymbol{u}}}&&\text{in }L^{2}(D)^{2},\\ &\overline{p}_{n}\rightharpoonup\overline{p}&&\text{in }L^{2}(D).\end{aligned}\right.

(6)

Passing through the limit. We shall now show that the limit $(\overline{\operatorname{\boldsymbol{u}}},\overline{p})$ , when restricted to the domain $\Omega$ , solves (4). Indeed, since $(\tilde{\operatorname{\boldsymbol{u}}}_{n},{p}_{n})=(\overline{\operatorname{\boldsymbol{u}}}_{n},\overline{p}_{n})\big{|}_{\Omega_{n}}$ solves (4) in $\Omega_{n}$ , then for any $(\operatorname{\boldsymbol{\varphi}},q)\in\operatorname{\mathbf{H}}^{1}_{\Gamma_{0}}(D)\times L^{2}(D)$

\left\{\begin{aligned} \nu(\chi_{n}\nabla\overline{\operatorname{\boldsymbol{u}}}_{n},\nabla\operatorname{\boldsymbol{\varphi}})_{D}+b(\operatorname{\boldsymbol{\varphi}},\chi_{n}\overline{p}_{n})_{D}&=(\chi_{n}\operatorname{\boldsymbol{f}},\operatorname{\boldsymbol{\varphi}})_{D}-\nu(\chi_{n}\nabla\operatorname{\boldsymbol{g}},\nabla\operatorname{\boldsymbol{\varphi}})_{D},\\ b(\overline{\operatorname{\boldsymbol{u}}}_{n},\chi_{n}q)_{D}&=0.\end{aligned}\right.

(7)

Furthermore, since $\chi_{n}:=\chi_{\Omega_{n}}\operatorname{\mathrel{\ensurestackMath{\stackon[1pt]{\rightharpoonup}{\scriptstyle\ast}}}}\chi:=\chi_{\Omega}$ in $\mathcal{U}_{ad}$ and from (6), we can easily show that

\left\{\begin{aligned} \nu(\chi\nabla\overline{\operatorname{\boldsymbol{u}}},\nabla\operatorname{\boldsymbol{\varphi}})_{D}+b(\operatorname{\boldsymbol{\varphi}},\chi\overline{p})_{D}&=(\chi\operatorname{\boldsymbol{f}},\operatorname{\boldsymbol{\varphi}})_{D}-\nu(\chi\nabla\operatorname{\boldsymbol{g}},\nabla\operatorname{\boldsymbol{\varphi}})_{D},\\ b(\overline{\operatorname{\boldsymbol{u}}},\chi q)_{D}&=0.\end{aligned}\right.

(8)

Thus, $(\tilde{\operatorname{\boldsymbol{u}}},p)=(\overline{\operatorname{\boldsymbol{u}}},\overline{p})\big{|}_{\Omega}$ is a solution to (4). ∎

As usual in proving the existence of solution for minimization problems, the arguments will be based on the lower semicontinuity, and the existence of a minimizing sequence based on the boundedness below of the objective functional.

Note that the objective functional $\mathcal{G}$ can be dissected into three different integrals, namely $J_{1}(\operatorname{\boldsymbol{u}},\Omega)=\frac{\gamma_{1}}{2}\|\nabla\times\operatorname{\boldsymbol{u}}\|_{L^{2}(\Omega)}^{2}$ , $J_{2}({\operatorname{\boldsymbol{u}}},\Omega)={\gamma_{2}}\|h(\mathrm{det}(\nabla{\operatorname{\boldsymbol{u}}}))\|_{L^{1}(\Omega)}$ , and $P(\Omega)=\int_{\Gamma_{\rm f}}\operatorname*{d\!}s$ . Thus, to establish that $\mathcal{G}$ is bounded below, we can show that $J_{1}$ and $J_{2}$ are uniformly bounded. This implies that, since the boundary $\Gamma_{\rm f}$ is strictly inside the bounded hold-all domain $D$ , $C\leq\mathcal{G}$ for some constant $C$ .

Note that $J_{1}$ can be estimated by the $H^{1}(\Omega)^{2}$ norm of the state ${\operatorname{\boldsymbol{u}}}$ , i.e.,

\displaystyle J_{1}({\operatorname{\boldsymbol{u}}},\Omega)

\displaystyle=\frac{\gamma_{1}}{2}\|\nabla\times{\operatorname{\boldsymbol{u}}}\|_{L^{2}(\Omega)}^{2}\leq\frac{\gamma_{1}}{2}\int_{\Omega}|\nabla{\operatorname{\boldsymbol{u}}}|^{2}\operatorname*{d\!}x\leq\frac{\gamma_{1}}{2}\|{\operatorname{\boldsymbol{u}}}\|^{2}_{H^{1}(\Omega)^{2}}.

Meanwhile, since $t^{3}/(t^{2}+1)\leq t$ for any $t$ , we get the following estimate,

J_{2}({\operatorname{\boldsymbol{u}}},\Omega)\leq\gamma_{2}\|\mathrm{det}(\nabla{\operatorname{\boldsymbol{u}}})\|_{L^{1}(\Omega)}\leq\frac{1}{2}\|{\operatorname{\boldsymbol{u}}}\|_{H^{1}(\Omega)^{2}}^{2},

where Young’s inequality has been employed for the last inequality.

From Proposition 1, we have established the continuity of $\Omega\mapsto\tilde{\operatorname{\boldsymbol{u}}}$ , which also implies the continuity of the map $\Omega\mapsto{\operatorname{\boldsymbol{u}}}=\tilde{\operatorname{\boldsymbol{u}}}+\operatorname{\boldsymbol{g}}\in H^{1}(\Omega)^{2}$ . Thus, from the recently established estimates for $J_{1}$ and $J_{2}$ with respect to ${\operatorname{\boldsymbol{u}}}$ , we can infer that the respective functionals are continuous (and hence upper semicontinuous) with respect to the elements of $\operatorname{\mathcal{O}_{ad}}\,$ . Furthermore, the functionals are uniformly bounded, since for each $i=1,2$ ,

\displaystyle J_{i}({\operatorname{\boldsymbol{u}}},\Omega)

\displaystyle\leq c\|{\operatorname{\boldsymbol{u}}}\|_{H^{1}(\Omega)^{2}}^{2}=c\|\tilde{\operatorname{\boldsymbol{u}}}+\operatorname{\boldsymbol{g}}\|_{H^{1}(\Omega)^{2}}^{2}\leq\tilde{c}(\|\operatorname{\boldsymbol{f}}\|_{L^{2}(D)^{2}}+\|\operatorname{\boldsymbol{g}}\|_{H^{1}(D)^{2}})^{2}.

Given all these facts, we can now show the existence of an optimal shape rendering our problem well-posed.

Theorem 2.4.

Let $\alpha,\gamma_{1},\gamma_{2}>0$ , ${\operatorname{\boldsymbol{f}}}\in L^{2}(D)^{2}$ , and $\operatorname{\boldsymbol{g}}\in H^{2}(D)^{2}$ satisfy (3). Then, there exists $\Omega^{*}\in\operatorname{\mathcal{O}_{ad}}$ such that

\mathcal{J}(\Omega^{*})=\min_{\Omega\in\operatorname{\mathcal{O}_{ad}}}\mathcal{J}(\Omega).

Proof.

From the fact that $\mathcal{J}$ is bounded below, there exists a sequence $\{\Omega_{n}\}\subset\operatorname{\mathcal{O}_{ad}}$ such that

\mathcal{J}^{*}:=\inf_{\Omega\in\operatorname{\mathcal{O}_{ad}}}\mathcal{J}(\Omega)=\lim_{n\to\infty}\mathcal{J}(\Omega_{n}).

As a consequence of Lemma 2.2, there exists a subsequence of $\{\Omega_{n}\}$ , which shall be denoted in the same manner, and a domain $\Omega^{*}$ such that $\Omega_{n}\to\Omega^{*}$ in $\operatorname{\mathcal{O}_{ad}}\,$ . By definition, $\mathcal{J}^{*}\leq\mathcal{J}(\Omega^{*})$ . However, since $-J_{1}$ , $-J_{2}$ , and $P$ (for $P$ , see Proposition 2.3.7 in [17]) are lower semicontinuous,

	$\displaystyle\mathcal{J}(\Omega^{*})$	$\displaystyle=P(\Omega^{})-(J_{1}(\Omega^{})+J_{2}(\Omega^{*}))\leq\liminf_{n\to\infty}[P(\Omega_{n})-(J_{1}(\Omega_{n})+J_{2}(\Omega_{n}))]$
		$\displaystyle=\liminf_{n\to\infty}\mathcal{J}(\Omega_{n})=\mathcal{J}^{*}.$

Therefore, $\Omega^{*}$ is our desired minimizer. ∎

3. Shape Sensitivity Analysis

In this section, given the previously established existence of the optimal shape, we shall discuss the necessary condition for the optimization problem. This is done by investigating the sensitivity of the objective functional $\mathcal{G}$ with respect to shape variations. We start this section by introducing the velocity method, where we consider domain variations generated by a given velocity.

3.1. Identity perturbation

In what follows, we consider a family of autonomous deformation fields ${\operatorname{\boldsymbol{\theta}}}$ belonging to ${\operatorname{\mathbf{\Theta}}}:=\{{\operatorname{\boldsymbol{\theta}}}\in C^{1,1}(\overline{D};\mathbb{R}^{2});{\operatorname{\boldsymbol{\theta}}}=0\text{ on }\partial D\cup\Gamma_{\rm in}\cup\Gamma_{\rm w}\}$ . An element ${\operatorname{\boldsymbol{\theta}}}\in\operatorname{\mathbf{\Theta}}$ generates an identity perturbation operator $T_{t}:\overline{D}\to\mathbb{R}^{2}$ defined by

\displaystyle\begin{aligned} T_{t}(x)=x+t{\operatorname{\boldsymbol{\theta}}}(x),\quad\forall x\in\overline{D}.\end{aligned}

(9)

With this operator, a domain $\Omega\subset D$ is perturbed so that for some $\tau:=\tau({\operatorname{\boldsymbol{\theta}}})>0$ we have a family of perturbed domains $\{\Omega_{t}=T_{t}(\Omega);t<\tau\}$ . Here, $\tau$ is chosen so that $\mathrm{det}\nabla T_{t}>0$ , where the following statement shows that $\tau$ is indeed dependent to the velocity field ${\operatorname{\boldsymbol{\theta}}}$ .

Lemma 3.1 (cf. [2, 14]).

Let ${\operatorname{\boldsymbol{\theta}}}\in\operatorname{\mathbf{\Theta}}$ and $T_{t}$ be the generated transformation by means of (9) and denote $J_{t}:=\mathrm{det}\nabla T_{t}$ . Then, we have the following:

(i)

$J_{t}=1+t\operatorname{\mathrm{div}}{\operatorname{\boldsymbol{\theta}}}+t^{2}\mathrm{det}\nabla{\operatorname{\boldsymbol{\theta}}}$ ;
(ii)

there exists $\tau_{0}=\tau_{0}({\operatorname{\boldsymbol{\theta}}}),\,\alpha_{1},\alpha_{2}>0$ independent of $t$ and $x$ such that

$0<\alpha_{1}\leq J_{t}(x)\leq\alpha_{2},\quad\forall t\in[0,\tau_{0}],\forall x\in D.$

We further mention that by the definition of ${\operatorname{\boldsymbol{\theta}}}$ , i.e., ${\operatorname{\boldsymbol{\theta}}}\equiv 0$ on $\Gamma_{\rm out}$ , $\Gamma_{\rm in}$ , and $\Gamma_{\rm w}$ , then these boundaries are part of the perturbed domains $\Omega_{t}$ . To be precise, we have

\partial\Omega_{t}=\Gamma_{\rm out}\cup\Gamma_{\rm in}\cup\Gamma_{\rm w}\cup T_{t}(\Gamma_{\rm f}).

Additionally, a domain that has at most $C^{1,1}$ regularity preserves its said regularity with this transformation, this means that for $0\leq t\leq\tau$ , $\Omega_{t}$ has $C^{1,1}$ regularity given that the initial domain $\Omega$ is a $C^{1,1}$ domain. Lastly, we note that $\{\Omega_{t}=T_{t}(\Omega);t<\tau\}\subset\operatorname{\mathcal{O}_{ad}}$ due to Remark 3.

Before we move further in this exposition, let us look at some vital properties of $T_{t}$ .

Lemma 3.2 (cf [6, 34]).

Let ${\operatorname{\boldsymbol{\theta}}}\in\operatorname{\mathbf{\Theta}}$ , then for sufficiently small $\tau>0$ , the map $T_{t}$ defined in (9) satisfies the following properties:

$\bullet$

$[t\mapsto T_{t}]\in C^{1}([0,t_{0}];C^{2,1}(\overline{D},\mathbb{R}^{2}));\ \quad\bullet\ [t\mapsto T_{t}^{-1}]\in C([0,t_{0}];C^{2,1}(\overline{D},\mathbb{R}^{2}));$
$\bullet$

$[t\mapsto J_{t}]\in C^{1}([0,t_{0}];C^{1,1}(\overline{D}));\quad\qquad\!\bullet\ M_{t},M_{t}^{\top}\in C^{1,1}(\overline{D},\mathbb{R}^{2\times 2});$
$\bullet$

$\frac{d}{dt}J_{t}\big{|}_{t=0}=\operatorname{\mathrm{div}}\operatorname{\boldsymbol{\theta}};\quad\qquad\qquad\qquad\quad\ \bullet\ \frac{d}{dt}M_{t}\big{|}_{t=0}=-\nabla\operatorname{\boldsymbol{\theta}},$

where $M_{t}(x)=(\nabla T_{t}(x))^{-1}$ .

Let us also recall Hadamard’s identity which will be indespensible for the discussion of the necessary conditions.

Lemma 3.3.

Let $f\in C([0,\tau];W^{1,1}(D))$ and suppose that $\frac{\partial f}{\partial t}(0)\in L^{1}(D)$ , then

\displaystyle\frac{d}{dt}\int_{\Omega_{t}}f(t,x)\operatorname*{d\!}x\Big{|}_{t=0}=\int_{\Omega}\frac{\partial f}{\partial t}(0,x)\operatorname*{d\!}x+\int_{\Gamma_{\rm f}}f(0,x){\operatorname{\boldsymbol{\theta}}}\cdot\operatorname{\boldsymbol{n}}\operatorname*{d\!}s.

Proof.

See Theorem 5.2.2 and Proposition 5.4.4 of [17]. ∎

3.2. Rearrangement Method

In this part, we determine the shape derivative of the objective functional with respect to the variations generated by the transformation $T_{t}$ . This approach gets rid of the tedious process of solving first the shape derivative of the state solutions, then solving the shape derivative of the objective functional.

To start with, we consider a Hilbert space $Y(\Omega)$ and an operator

F:Y(\Omega)\times\operatorname{\mathcal{O}_{ad}}\to Y(\Omega)^{\prime},

where the equation $\langle F(y,\Omega),\phi\rangle_{Y(\Omega)^{\prime}\times Y(\Omega)}=0$ corresponds to a variational problem in $\Omega.$

Suppose that the free-boundary is denoted by $\Gamma_{\rm f}\subset\partial\Omega$ , the said method deals with the shape optimization

\min_{\Omega\in\operatorname{\mathcal{O}_{ad}}}J(y,\Omega):=\int_{\Omega}j(y)dx+\alpha\int_{\Gamma_{\rm f}}\operatorname*{d\!}s.

subject to

\displaystyle F(y,\Omega)=0\text{ in }Y(\Omega)^{\prime}.

(10)

We define the Eulerian derivative of $J$ at $\Omega$ in the direction ${\operatorname{\boldsymbol{\theta}}}\in\operatorname{\mathbf{\Theta}}$ by

dJ(y,\Omega){\operatorname{\boldsymbol{\theta}}}=\lim_{t\searrow 0}\frac{J(y_{t},\Omega_{t})-J(y,\Omega)}{t},

where $y_{t}$ solves the equation $F(y_{t},\Omega_{t})=0$ in $Y(\Omega_{t})^{\prime}$ . If $dJ(y,\Omega){\operatorname{\boldsymbol{\theta}}}$ exists for all ${\operatorname{\boldsymbol{\theta}}}\in\operatorname{\mathbf{\Theta}}$ and that $dJ(y,\Omega)$ defines a bounded linear functional on $\operatorname{\mathbf{\Theta}}$ then we say that $J$ is shape differentiable at $\Omega$ .

The so-called rearrangement method is given as below:

Lemma 3.4 (cf [20]).

Suppose $j\in C^{1,1}(\mathbb{R}^{2},\mathbb{R})$ and that the following assumptions hold:

(A1)

There exists an operator $\tilde{F}:Y(\Omega)\times[0,\tau]\to Y(\Omega)^{\prime}$ such that $F(y_{t},\Omega_{t})=0$ in $Y(\Omega_{t})^{\prime}$ is equivalent to

$\displaystyle\tilde{F}(y^{t},t)=0\text{ in }Y(\Omega)^{\prime},$ (11)

with $\tilde{F}(y,0)=F(y,\Omega)$ for all $y\in Y(\Omega).$

(A2)

Let $y,v\in Y(\Omega)$ . Then $F_{y}(y,\Omega)\in\mathcal{L}(Y(\Omega),Y(\Omega)^{\prime})$ satisfies

\langle F(v,\Omega)-F(y,\Omega)-F_{y}(y,\Omega)(v-y),z\rangle_{Y(\Omega)^{\prime}\times Y(\Omega)}=\mathcal{O}(\|v-y\|_{Y(\Omega)}^{2}),

for all $z\in Y(\Omega).$

(A3)

Let $y\in Y(\Omega)$ be the unique solution of (10). Then for any $f\in Y(\Omega)^{\prime}$ the solution of the following linearized equation exists:

\displaystyle\langle F_{y}(y,\Omega)\delta y,z\rangle_{Y(\Omega)\times Y(\Omega)^{\prime}}=\langle f,z\rangle_{Y(\Omega)\times Y(\Omega)^{\prime}}\text{ for all }z\in Y(\Omega).

(A4)

Let $y^{t},y\in Y(\Omega)$ be the solutions of (11) and (10), respectively. Then $\tilde{F}$ and $F$ satisfy

\lim_{t\searrow 0}\frac{1}{t}\langle\tilde{F}(y^{t},t)-\tilde{F}(y,t)-F(y^{t},\Omega)+F(y,\Omega),z\rangle_{Y(\Omega)^{\prime}\times Y(\Omega)}=0

for all $z\in Y(\Omega).$

Let $y\in Y(\Omega)$ be the solution of (10), and suppose that the adjoint equation, for all $z\in Y(\Omega)$

\displaystyle\langle F_{y}(y,\Omega)z,p\rangle_{Y(\Omega)^{\prime}\times Y(\Omega)}=(j^{\prime}(y),z)_{L^{2}(\Omega)}

(12)

has a unique solution $p\in Y(\Omega)$ . Then the Eulerian derivative of $J$ at $\Omega$ in the direction ${\operatorname{\boldsymbol{\theta}}}\in\operatorname{\mathbf{\Theta}}$ exists and is given by

\displaystyle\begin{aligned} dJ(y,\Omega){\operatorname{\boldsymbol{\theta}}}=&\,-\frac{d}{dt}\langle\tilde{F}(y,t),p\rangle_{Y(\Omega)^{\prime}\times Y(\Omega)}\Big{|}_{t=0}\\ &\,+\int_{\Omega}j(y)\operatorname{\mathrm{div}}{\operatorname{\boldsymbol{\theta}}}\operatorname*{d\!}x+\alpha\int_{\Gamma_{\rm f}}\kappa{\operatorname{\boldsymbol{\theta}}}\cdot\operatorname{\boldsymbol{n}}\operatorname*{d\!}s,\end{aligned}

(13)

where $\kappa$ is the mean curvature of $\Gamma_{\rm f}$ .

Before we go further, we note that in usual practice the element $y^{t}\in Y(\Omega)$ is solved using the composition $y^{t}:=y_{t}\circ T_{t}$ , and utilizes the function space parametrization (see [6, Chapter 10 Section 2.2] for example, for further details). Furthermore, assumption (A3) is originally written as the Hölder continuity of solutions $y^{t}\in Y(\Omega)$ of (11) with respect to the time parameter $t\in[0,\tau]$ . Fortunately, Ito, K., et.al.[20], have shown that assumption (A3) implies the continuity. We cite the said result in the following lemma.

Lemma 3.5.

Suppose that $y\in Y(\Omega)$ solves (10) and $y^{t}\in Y(\Omega)$ is the solution to (11). Assume furthermore that (A3) holds. Then, $\|y^{t}-y\|_{X(\Omega)}=o(t^{\frac{1}{2}})$ as $t\searrow 0$ .

We applying the rearrangement method to the velocity-pressure operator $E(\cdot)_{\Omega}:X(\Omega)\to X(\Omega)^{\prime}$ defined by

	$\displaystyle\langle E({\operatorname{\boldsymbol{u}}},p)_{\Omega},(\operatorname{\boldsymbol{\varphi}},\psi)\rangle_{X(\Omega)^{\prime}\times X(\Omega)}=$	$\displaystyle\,\nu a(\operatorname{\boldsymbol{u}},\operatorname{\boldsymbol{\varphi}})_{\Omega}+b(\operatorname{\boldsymbol{\varphi}},p)_{\Omega}+b({\operatorname{\boldsymbol{u}}},q)_{\Omega}$
		$\displaystyle-(\operatorname{\boldsymbol{f}},\operatorname{\boldsymbol{\varphi}})_{\Omega}+a(\operatorname{\boldsymbol{g}},\operatorname{\boldsymbol{\varphi}})_{\Omega}$

where $X(\Omega):=\operatorname{\mathbf{H}}^{1}_{\Gamma_{0}}(\Omega)\times L^{2}(\Omega)$ .

According to Theorem 1.1, there exists a unique $(\tilde{\operatorname{\boldsymbol{u}}},p)\in X(\Omega)$ such that for any $(\operatorname{\boldsymbol{\varphi}},\psi)\in X(\Omega)$

\displaystyle\langle E(\tilde{\operatorname{\boldsymbol{u}}},p)_{\Omega},(\operatorname{\boldsymbol{\varphi}},\psi)\rangle_{X(\Omega)^{\prime}\times X(\Omega)}=0.

(14)

Notation.

Moving forward we shall use the following notations $X=X(\Omega)$ , $X_{t}=X(\Omega_{t})$ .

Our goal is to characterize, and of course show the existence of the Eulerian derivative of the objective functional

\displaystyle\mathcal{G}({\operatorname{\boldsymbol{u}}},\Omega)=\alpha\int_{\Gamma_{\rm f}}\operatorname*{d\!}s-\int_{\Omega}\frac{\gamma_{1}}{2}|\nabla\times{\operatorname{\boldsymbol{u}}}|^{2}+\gamma_{2}g(\mathrm{det}(\nabla{\operatorname{\boldsymbol{u}}}))\operatorname*{d\!}x,

(15)

where ${\operatorname{\boldsymbol{u}}}=\tilde{\operatorname{\boldsymbol{u}}}+\operatorname{\boldsymbol{g}}$ , and $\tilde{\operatorname{\boldsymbol{u}}}\in\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(\Omega)$ is the first component of the solution of (14).

From the deformation field $T_{t}$ , we let $(\tilde{\operatorname{\boldsymbol{u}}}_{t},p_{t})\in X_{t}$ be the solution of the equation

\displaystyle\langle E(\tilde{\operatorname{\boldsymbol{u}}}_{t},p_{t})_{\Omega_{t}},(\operatorname{\boldsymbol{\varphi}}_{t},\psi_{t})\rangle_{X_{t}^{\prime}\times X_{t}}=0,\quad\forall(\operatorname{\boldsymbol{\varphi}}_{t},\psi_{t})\in X_{t}.

(16)

By perturbing the equation (16) back to the reference domain $\Omega$ we get the operator $\tilde{E}:X\times[0,\tau]\to X^{\prime}$ defined by

\displaystyle\begin{aligned} \langle\tilde{E}((\operatorname{\boldsymbol{u}},p),t),(\operatorname{\boldsymbol{\varphi}},\psi)\rangle_{X^{\prime}\times X}&=\nu a_{t}({\operatorname{\boldsymbol{u}}},\operatorname{\boldsymbol{\varphi}})+b_{t}(\operatorname{\boldsymbol{\varphi}},p)+b_{t}({\operatorname{\boldsymbol{u}}},\psi)\\ &-(J_{t}\operatorname{\boldsymbol{f}}^{t},\operatorname{\boldsymbol{\varphi}})_{\Omega}+\nu a_{t}({\operatorname{\boldsymbol{g}}},\operatorname{\boldsymbol{\varphi}})\end{aligned}

(17)

where $\operatorname{\boldsymbol{f}}^{t}=\operatorname{\boldsymbol{f}}\circ T_{t}\in L^{2}(\Omega)^{2}$ and – by denoting $(M_{t}^{\top})_{k}$ the $k^{th}$ row of $M_{t}^{\top}$ , $\operatorname{\boldsymbol{v}}_{k}$ the $k^{th}$ component of a vector $\operatorname{\boldsymbol{v}}$ , and by using Einstein convention over $k=1,2$ – the bilinear forms are defined as follows

	$\displaystyle a_{t}({\operatorname{\boldsymbol{u}}},\operatorname{\boldsymbol{\varphi}})=$	$\displaystyle\,(J_{t}M_{t}\nabla{\operatorname{\boldsymbol{u}}},M_{t}\nabla\operatorname{\boldsymbol{\varphi}})_{\Omega},$
	$\displaystyle b_{t}(\operatorname{\boldsymbol{v}},\psi)=$	$\displaystyle\,-(J_{t}\psi,(M_{t}^{\top})_{k}\nabla\operatorname{\boldsymbol{v}}_{k})_{\Omega}.$

Using similar arguments as in [20], it can be easily shown that $E(\cdot,\cdot)_{\Omega}$ and $\tilde{E}$ satisfy (A1) and (A4).

Furthermore, the Fréchet derivative $E^{\prime}\in\mathcal{L}(X,\mathcal{L}(X,X^{\prime}))$ of the velocity-pressure operator at a point $({\operatorname{\boldsymbol{u}}},p)\in X$ can be easily obtained as

\langle E^{\prime}({\operatorname{\boldsymbol{u}}},p)(\delta{\operatorname{\boldsymbol{u}}},\delta p),(\operatorname{\boldsymbol{\varphi}},\psi)\rangle_{X^{\prime}\times X}=\nu a(\delta{\operatorname{\boldsymbol{u}}},\operatorname{\boldsymbol{\varphi}})_{\Omega}+b(\operatorname{\boldsymbol{\varphi}},\delta p)_{\Omega}+b(\delta{\operatorname{\boldsymbol{u}}},q)_{\Omega},

for all $(\operatorname{\boldsymbol{\varphi}},\psi)\in X$ . Due to the linearity, we further infer that

\langle E({\operatorname{\boldsymbol{v}}},q)_{\Omega}-E({\operatorname{\boldsymbol{u}}},p)_{\Omega}-E^{\prime}({\operatorname{\boldsymbol{u}}},p)({\operatorname{\boldsymbol{v}}}-{\operatorname{\boldsymbol{u}}},q-p),(\operatorname{\boldsymbol{\varphi}},\psi)\rangle_{X^{\prime}\times X}=0

for any $({\operatorname{\boldsymbol{v}}},q),({\operatorname{\boldsymbol{u}}},p),(\operatorname{\boldsymbol{\varphi}},\psi)\in X$ . From the coercivity of $a(\cdot,\cdot)_{\Omega}$ , and because $b(\cdot,\cdot)_{\Omega}$ satisfies the inf-sup condition, then for any $\Phi\in X^{\prime}$ the variational problem

\langle E^{\prime}(\tilde{\operatorname{\boldsymbol{u}}},p)(\delta{\operatorname{\boldsymbol{u}}},\delta p),(\operatorname{\boldsymbol{\varphi}},\psi)\rangle_{X^{\prime}\times X}=\langle\Phi,(\operatorname{\boldsymbol{\varphi}},\psi)\rangle_{X^{\prime}\times X},\quad\forall(\operatorname{\boldsymbol{\varphi}},\psi)\in X,

is well-posed. These imply that the operators $E(\cdot,\cdot)_{\Omega}$ and $E^{\prime}({\operatorname{\boldsymbol{u}}},p)$ satisfy assumptions (A2) and (A3).

We now characterize the Eulerian derivative of $\mathcal{G}$ by utilizing Lemma 3.4.

Theorem 3.6.

Suppose that $\Omega\subset D$ is a domain that is of class $C^{1,1}$ , and let $(\tilde{\operatorname{\boldsymbol{u}}},p)\in(\operatorname{\mathbf{H}}_{\Gamma_{0}}^{1}(\Omega)\cap H^{2}(\Omega)^{2})\times L^{2}(\Omega)$ be the solution to (14). Then, the Eulerian derivative of $\mathcal{G}$ exists and is given by

	$\displaystyle d\mathcal{G}((\operatorname{\boldsymbol{u}},p),\Omega){\operatorname{\boldsymbol{\theta}}}=$	$\displaystyle\,\int_{\Gamma_{\rm f}}\bigg{[}\alpha\kappa-\frac{\gamma_{1}}{2}\|\nabla\times{\operatorname{\boldsymbol{u}}}\|^{2}-\gamma_{2}h(\mathrm{det}(\nabla{\operatorname{\boldsymbol{u}}}))$
		$\displaystyle+\frac{\partial{\operatorname{\boldsymbol{u}}}}{\partial{\operatorname{\boldsymbol{n}}}}\left(\frac{\partial{\operatorname{\boldsymbol{v}}}}{\partial{\operatorname{\boldsymbol{n}}}}+\gamma_{1}(\nabla\times{\operatorname{\boldsymbol{u}}}){\operatorname{\boldsymbol{\tau}}}+\gamma_{2}P({\operatorname{\boldsymbol{u}}})\right)\bigg{]}{\operatorname{\boldsymbol{\theta}}}\cdot{\operatorname{\boldsymbol{n}}}\operatorname*{d\!}s,$

where $\operatorname{\boldsymbol{\tau}}$ is the unit tangential vector on $\Gamma_{\rm f}$ , ${\operatorname{\boldsymbol{u}}}=\tilde{\operatorname{\boldsymbol{u}}}+{\operatorname{\boldsymbol{g}}}\in H^{2}(\Omega)^{2}$ , $(\operatorname{\boldsymbol{v}},\pi)\in X$ solves the adjoint equation

\displaystyle\langle E^{\prime}(\tilde{\operatorname{\boldsymbol{u}}},p)({\operatorname{\boldsymbol{\varphi}}},\psi),({\operatorname{\boldsymbol{v}}},\pi)\rangle_{X^{\prime}\times X}=-(\gamma_{1}\vec{\nabla}\times(\nabla\times{\operatorname{\boldsymbol{u}}})+\gamma_{2}R(\operatorname{\boldsymbol{u}}),{\operatorname{\boldsymbol{\varphi}}})_{\Omega},

(18)

for any $(\operatorname{\boldsymbol{\varphi}},\psi)\in X$ ,

R(\operatorname{\boldsymbol{u}})=\begin{bmatrix}-\nabla\times(h^{\prime}(\mathrm{det}(\nabla{\operatorname{\boldsymbol{u}}}))\nabla u_{2})\\ \nabla\times(h^{\prime}(\mathrm{det}(\nabla{\operatorname{\boldsymbol{u}}}))\nabla u_{1})\end{bmatrix}\!,P(\operatorname{\boldsymbol{u}})=\begin{bmatrix}h^{\prime}(\mathrm{det}(\nabla{\operatorname{\boldsymbol{u}}}))\left(\frac{\partial u_{2}}{\partial x_{2}}n_{x_{1}}-\frac{\partial u_{2}}{\partial x_{1}}n_{x_{2}}\right)\\ h^{\prime}(\mathrm{det}(\nabla{\operatorname{\boldsymbol{u}}}))\left(\frac{\partial u_{1}}{\partial x_{1}}n_{x_{2}}-\frac{\partial u_{1}}{\partial x_{2}}n_{x_{1}}\right)\end{bmatrix}

and the curl of a scalar valued function is defined by $\vec{\nabla}\times\psi=\left(\frac{\partial\psi}{\partial x_{2}},-\frac{\partial\psi}{\partial x_{1}}\right)^{\top}.$

Remark 4.

In strong form, we can write the adjoint equation (18) as follows:

\displaystyle\left\{\begin{aligned} \,-\nu\Delta{\operatorname{\boldsymbol{v}}}+\nabla\pi&=-[\gamma_{1}\vec{\nabla}\times(\nabla\times{\operatorname{\boldsymbol{u}}})+\gamma_{2}R({\operatorname{\boldsymbol{u}}})]&&\text{ in }\Omega,\\ \,\operatorname*{div}{\operatorname{\boldsymbol{v}}}&=0&&\text{ in }\Omega,\\ \,{\operatorname{\boldsymbol{v}}}&=0&&\text{ on }\partial\Omega\backslash\Gamma_{\rm out},\\ \,-\pi{\operatorname{\boldsymbol{n}}}+\nu\partial_{\operatorname{\boldsymbol{n}}}{\operatorname{\boldsymbol{v}}}&=0&&\text{ on }\Gamma_{\rm out}.\end{aligned}\right.

(19)

Proof of Theorem 3.6 From the symmetry of the operator $E^{\prime}(\tilde{\operatorname{\boldsymbol{u}}},p)\in\mathcal{L}(X^{\prime},X)$ , the unique existence of the solution to the adjoint equation (18) can be established easily using the same arguments as in (4) – this of course uses the fact that $\operatorname{\boldsymbol{u}}\in H^{2}(\Omega)^{2}$ , and hence the well-definedness of the right-hand side of (18). Thus, the existence of the Eulerian derivative of $\mathcal{G}$ is assured.

We now characterize the derivative $d\mathcal{G}(({\operatorname{\boldsymbol{u}}},p),\Omega){\operatorname{\boldsymbol{\theta}}}$ . We begin by solving for the expression $\frac{d}{dt}\langle\tilde{E}((\tilde{\operatorname{\boldsymbol{u}}},p),t),({\operatorname{\boldsymbol{v}}},\pi)\rangle_{X^{\prime}\times X}\big{|}_{t=0}$ . By pushing $\tilde{E}(({\operatorname{\boldsymbol{u}}},p),t)$ toward the perturbed domain $\Omega_{t}$ and by using Lemma 3.3 yield the following

	$\displaystyle\frac{d}{dt}\langle\tilde{E}((\tilde{\operatorname{\boldsymbol{u}}},p),t),$	$\displaystyle({\operatorname{\boldsymbol{v}}},\pi)\rangle_{X^{\prime}\times X}\Big{\|}_{t=0}$
	$\displaystyle=$	$\displaystyle\,[\nu(\nabla{\operatorname{\boldsymbol{u}}}:\nabla{\operatorname{\boldsymbol{v}}},{\operatorname{\boldsymbol{\theta}}}\cdot{\operatorname{\boldsymbol{n}}})_{\Gamma_{\rm f}}-(p\operatorname{\mathrm{div}}{\operatorname{\boldsymbol{v}}},{\operatorname{\boldsymbol{\theta}}}\cdot{\operatorname{\boldsymbol{n}}})_{\Gamma_{\rm f}}-(\psi\operatorname{\mathrm{div}}{\operatorname{\boldsymbol{u}}},{\operatorname{\boldsymbol{\theta}}}\cdot{\operatorname{\boldsymbol{n}}})_{\Gamma_{\rm f}}]$
		$\displaystyle+[\nu a({\operatorname{\boldsymbol{u}}},\phi_{v})_{\Omega}+b(\phi_{v},p)_{\Omega}+b({\operatorname{\boldsymbol{u}}},\phi_{\pi})_{\Omega}-(\operatorname{\boldsymbol{f}}\phi_{v})_{\Omega}]$
		$\displaystyle+[\nu a(\phi_{u},{\operatorname{\boldsymbol{v}}})_{\Omega}+b(\phi_{u},\pi)_{\Omega}+b({\operatorname{\boldsymbol{v}}},\phi_{q})_{\Omega}]$

where $\phi_{u}=-\nabla{\operatorname{\boldsymbol{u}}}^{\top}{\operatorname{\boldsymbol{\theta}}}$ , $\phi_{v}=-\nabla{\operatorname{\boldsymbol{v}}}^{\top}{\operatorname{\boldsymbol{\theta}}}$ , $\phi_{p}=-\nabla p\cdot{\operatorname{\boldsymbol{\theta}}}$ , and $\phi_{\pi}=-\nabla\pi\cdot{\operatorname{\boldsymbol{\theta}}}$ , which were generated by using the identity $\frac{d}{dt}(\varphi\circ T_{t}^{-1})=-\nabla\varphi\cdot{\operatorname{\boldsymbol{\theta}}},$ for $\phi\in H^{1}(D)$ [34].

Using Green’s identities, the divergence-free property of the variables ${\operatorname{\boldsymbol{u}}}$ and ${\operatorname{\boldsymbol{v}}}$ and Green’s curl identity [26, Theorem 3.29], we obtain the following derivative

	$\displaystyle\frac{d}{dt}\langle\tilde{E}$	$\displaystyle((\tilde{\operatorname{\boldsymbol{u}}},p),t),({\operatorname{\boldsymbol{v}}},\pi)\rangle_{X^{\prime}\times X}\Big{\|}_{t=0}$
	$\displaystyle=$	$\displaystyle\,-\left(\nu\partial_{\operatorname{\boldsymbol{n}}}{\operatorname{\boldsymbol{u}}}\cdot\partial_{\operatorname{\boldsymbol{n}}}{\operatorname{\boldsymbol{v}}},{\operatorname{\boldsymbol{\theta}}}\cdot{\operatorname{\boldsymbol{n}}}\right)_{\Gamma_{\rm f}}+\gamma_{1}\Big{[}(j_{1}^{\prime}({\operatorname{\boldsymbol{u}}}),\nabla{\operatorname{\boldsymbol{u}}}^{\top}{\operatorname{\boldsymbol{\theta}}})_{\Omega}-((\nabla\times{\operatorname{\boldsymbol{u}}}){\operatorname{\boldsymbol{\tau}}}\cdot\partial_{\operatorname{\boldsymbol{n}}}{\operatorname{\boldsymbol{u}}},{\operatorname{\boldsymbol{\theta}}}\cdot{\operatorname{\boldsymbol{n}}})_{\Gamma_{\rm f}}\Big{]}$
		$\displaystyle+\gamma_{2}\Big{[}(j_{2}^{\prime}({\operatorname{\boldsymbol{u}}}),\nabla{\operatorname{\boldsymbol{u}}}^{\top}{\operatorname{\boldsymbol{\theta}}})_{\Omega}-(P({\operatorname{\boldsymbol{u}}})\cdot\partial_{\operatorname{\boldsymbol{n}}}{\operatorname{\boldsymbol{u}}},{\operatorname{\boldsymbol{\theta}}}\cdot{\operatorname{\boldsymbol{n}}})_{\Gamma_{\rm f}}\Big{]}.$

Finally, we now use (13) to determine the derivative $d\mathcal{G}((\operatorname{\boldsymbol{u}},p),\Omega){\operatorname{\boldsymbol{\theta}}}$ to get

	$\displaystyle d\mathcal{G}((\operatorname{\boldsymbol{u}},p)$	$\displaystyle,\Omega){\operatorname{\boldsymbol{\theta}}}=\big{(}\partial_{\operatorname{\boldsymbol{n}}}{\operatorname{\boldsymbol{u}}}\cdot(\nu\partial_{\operatorname{\boldsymbol{n}}}{\operatorname{\boldsymbol{v}}}+(\nabla\times{\operatorname{\boldsymbol{u}}}){\operatorname{\boldsymbol{\tau}}}+P({\operatorname{\boldsymbol{u}}})),{\operatorname{\boldsymbol{\theta}}}\cdot{\operatorname{\boldsymbol{n}}}\big{)}_{\Gamma_{\rm f}}$
		$\displaystyle+\alpha(\kappa,{\operatorname{\boldsymbol{\theta}}}\cdot\operatorname{\boldsymbol{n}})_{\Gamma_{\rm f}}-\frac{\gamma_{1}}{2}\big{[}2(j_{1}^{\prime}({\operatorname{\boldsymbol{u}}}),\nabla{\operatorname{\boldsymbol{u}}}^{\top}{\operatorname{\boldsymbol{\theta}}})_{\Omega}+(j_{1}({\operatorname{\boldsymbol{u}}}),\operatorname{\mathrm{div}}{\operatorname{\boldsymbol{\theta}}})_{\Omega}\big{]}$
		$\displaystyle-\gamma_{2}\big{[}(j_{2}^{\prime}({\operatorname{\boldsymbol{u}}}),\nabla{\operatorname{\boldsymbol{u}}}^{\top}{\operatorname{\boldsymbol{\theta}}})_{\Omega}+(j_{2}({\operatorname{\boldsymbol{u}}}),\operatorname{\mathrm{div}}{\operatorname{\boldsymbol{\theta}}})_{\Omega}\big{]}.$

Before we write the final form of the shape derivative, we note that

\displaystyle\int_{\Omega}\operatorname{\mathrm{div}}(j_{1}({\operatorname{\boldsymbol{u}}}){\operatorname{\boldsymbol{\theta}}})\operatorname*{d\!}x=2(j_{1}^{\prime}({\operatorname{\boldsymbol{u}}}),\nabla{\operatorname{\boldsymbol{u}}}^{\top}{\operatorname{\boldsymbol{\theta}}})_{\Omega}+(j_{1}({\operatorname{\boldsymbol{u}}}),\operatorname{\mathrm{div}}{\operatorname{\boldsymbol{\theta}}})_{\Omega},

and

\displaystyle\int_{\Omega}\operatorname{\mathrm{div}}(j_{2}({\operatorname{\boldsymbol{u}}}){\operatorname{\boldsymbol{\theta}}})\operatorname*{d\!}x=(j_{2}^{\prime}({\operatorname{\boldsymbol{u}}}),\nabla{\operatorname{\boldsymbol{u}}}^{\top}{\operatorname{\boldsymbol{\theta}}})_{\Omega}+(j_{2}({\operatorname{\boldsymbol{u}}}),\operatorname{\mathrm{div}}{\operatorname{\boldsymbol{\theta}}})_{\Omega}.

Hence, we have the following form

	$\displaystyle d\mathcal{G}((\operatorname{\boldsymbol{u}},p),\Omega){\operatorname{\boldsymbol{\theta}}}=$	$\displaystyle\,\int_{\Gamma_{\rm f}}\left[\alpha\kappa+\frac{\partial{\operatorname{\boldsymbol{u}}}}{\partial{\operatorname{\boldsymbol{n}}}}\cdot\left(\frac{\partial{\operatorname{\boldsymbol{v}}}}{\partial{\operatorname{\boldsymbol{n}}}}+\gamma_{1}(\nabla\times{\operatorname{\boldsymbol{u}}}){\operatorname{\boldsymbol{\tau}}}+\gamma_{2}P({\operatorname{\boldsymbol{u}}})\right)\right]{\operatorname{\boldsymbol{\theta}}}\cdot\operatorname{\boldsymbol{n}}\operatorname*{d\!}s$
		$\displaystyle-\frac{\gamma_{1}}{2}\int_{\Omega}\operatorname{\mathrm{div}}(j_{1}({\operatorname{\boldsymbol{u}}}){\operatorname{\boldsymbol{\theta}}})\operatorname{d\!}x-\gamma_{2}\int_{\Omega}\operatorname{\mathrm{div}}(j_{2}({\operatorname{\boldsymbol{u}}}){\operatorname{\boldsymbol{\theta}}})\operatorname{d\!}x.$

Therefore, from the assumed regularity of the domain, and from the divergence theorem

	$\displaystyle d\mathcal{G}((\operatorname{\boldsymbol{u}},p),\Omega){\operatorname{\boldsymbol{\theta}}}=$	$\displaystyle\,\int_{\Gamma_{\rm f}}\left[\alpha\kappa+\frac{\partial{\operatorname{\boldsymbol{u}}}}{\partial{\operatorname{\boldsymbol{n}}}}\cdot\left(\frac{\partial{\operatorname{\boldsymbol{v}}}}{\partial{\operatorname{\boldsymbol{n}}}}+\gamma_{1}(\nabla\times{\operatorname{\boldsymbol{u}}}){\operatorname{\boldsymbol{\tau}}}+\gamma_{2}P({\operatorname{\boldsymbol{u}}})\right)\right]{\operatorname{\boldsymbol{\theta}}}\cdot\operatorname{\boldsymbol{n}}\operatorname*{d\!}s$
		$\displaystyle-\frac{\gamma_{1}}{2}\int_{\Gamma_{\rm f}}j_{1}({\operatorname{\boldsymbol{u}}}){\operatorname{\boldsymbol{\theta}}}\cdot{\operatorname{\boldsymbol{n}}}\operatorname{d\!}s-\gamma_{2}\int_{\Gamma_{\rm f}}j_{2}({\operatorname{\boldsymbol{u}}}){\operatorname{\boldsymbol{\theta}}}\cdot{\operatorname{\boldsymbol{n}}}\operatorname{d\!}s.$

Remark 5.

The resulting form of the shape derivative of $\mathcal{G}$ coincides with the Zolesio–Hadamard Structure Theorem [6, Corollary 9.3.1], i.e., we were able to write

d\mathcal{G}((\operatorname{\boldsymbol{u}},p),\Omega){\operatorname{\boldsymbol{\theta}}}=\int_{\Gamma_{\rm f}}\nabla{G}({\operatorname{\boldsymbol{\theta}}}\cdot{\operatorname{\boldsymbol{n}}})\operatorname*{d\!}s.

In this case, we shall call $\nabla G$ the shape gradient of $\mathcal{G}$ which will be useful in the numerical implementations that we shall illustrate in the proceeding sections.

4. Numerical Realization

In this section we shall use the recently formulated necessary condition to solve the minimization problem numerically. However, the said condition does not take into account the volume constraint. With this reason, we shall resort to two approaches, one is the augmented Lagrangian method, which will be based on the implementation of Dapogny, C., et.al.[5], and the other one is by generating divergence-free velocity fields.

4.1. Augmented Lagrangian Method

By writing the shape optimization problem (1) as the equality constrained optimization

\displaystyle\min_{\Omega\in\operatorname{\mathcal{O}_{ad}}}\mathcal{J}(\Omega)\text{ subject to }\mathcal{F}(\Omega):=|\Omega|-m=0,

we consider minimizing the augmented Lagrangian

\mathcal{L}(\Omega,\ell,b):=\mathcal{J}(\Omega)-\ell\mathcal{F}(\Omega)+\frac{b}{2}(\mathcal{F}(\Omega))^{2}.

In this definition, $\ell$ corresponds to the Lagrange multiplier with respect to the equality constraint, while $b>0$ is a penalty parameter. For a more detailed discussion on augmented Lagrangian methods one can refer to [28, Section 17.3], and to [5] for implementation to shape optimization.

We can solve the shape gradient of $\mathcal{L}$ at a given domain $\Omega\in\operatorname{\mathcal{O}_{ad}}$ - which we shall denote as $\nabla L$ - by utilizing the shape gradient of $\mathcal{J}$ and Lemma 3.3. This yields

\nabla L=\nabla G-\ell+b(|\Omega|-m),

where the term $|\Omega|-m$ is evaluated as $\int_{\Omega}\operatorname*{d\!}x-m.$

4.2. Gradient Descent and Normal Vector Regularization

Similarly with Remark 5, one can write the shape derivative of $\mathcal{L}$ as

d\mathcal{L}(\Omega,\ell,b){\operatorname{\boldsymbol{\theta}}}=\int_{\Gamma_{\rm f}}\nabla L({\operatorname{\boldsymbol{\theta}}}\cdot{\operatorname{\boldsymbol{n}}})\operatorname*{d\!}s.

This gives us an intuitive form of the velocity field as either $\operatorname{\boldsymbol{\theta}}=-\nabla G{\operatorname{\boldsymbol{n}}}$ (or $\operatorname{\boldsymbol{\theta}}=-\nabla L{\operatorname{\boldsymbol{n}}}$ for the augmented Lagrangian), which implies that

d\mathcal{G}(\Omega,\ell,b){\operatorname{\boldsymbol{\theta}}}=-\|\nabla G\|^{2}_{L^{2}(\Gamma_{\rm f})}\ (\text{or }d\mathcal{L}(\Omega,\ell,b){\operatorname{\boldsymbol{\theta}}}=-\|\nabla L\|^{2}_{L^{2}(\Gamma_{\rm f})})

However, such choice of $\operatorname{\boldsymbol{\theta}}$ may lead to unwanted behavior of the deformed domain especially for iterative schemes. To circumvent this issue, several authors utilized extensions of the velocity field in the domain $\Omega$ . In [5], the authors considered an elasticity extension with the tangential gradient for solving the field. Rabago et.al. [31] on the other hand cosidered a Neumann-type extension for solving the mentioned velocity field. In our context we shall consider two approaches based on $H^{1}$ gradient method [1]. In particular, for sufficiently small parameter $\gamma>0$ , and $\lambda\in\{0,1\}$ we shall use the following $H^{1}$ gradient method: Solve for $(\operatorname{\boldsymbol{\theta}},\vartheta)\in H^{1}(\Omega)^{2}\times L^{2}(\Omega)$ ( $\operatorname{\boldsymbol{\theta}}\in H^{1}(\Omega)^{2}$ if $\lambda=0$ ) that satisfies the variational equation

\displaystyle\left\{\begin{aligned} \gamma a({\operatorname{\boldsymbol{\theta}}},{\operatorname{\boldsymbol{\varphi}}})_{\Omega}+({\operatorname{\boldsymbol{\theta}}},{\operatorname{\boldsymbol{\varphi}}})_{\Omega}+\lambda b(\vartheta,{\operatorname{\boldsymbol{\varphi}}})_{\Omega}&=-(\nabla K{\operatorname{\boldsymbol{n}}},{\operatorname{\boldsymbol{\varphi}}})_{\Gamma_{\rm f}}\\ \lambda b(\psi,{\operatorname{\boldsymbol{\theta}}})_{\Omega}&=0\end{aligned}\right.

(20)

for all $(\operatorname{\boldsymbol{\varphi}},\psi)\in H^{1}(\Omega)^{2}\times L^{2}(\Omega)$ ( $\operatorname{\boldsymbol{\varphi}}\in H^{1}(\Omega)^{2}$ if $\lambda=0$ ). Here $\nabla K$ is $\nabla G$ if $\lambda=1$ , and $\nabla L$ if $\lambda=0$ .

Note that when $\lambda=0$ , (20) reduces to the usual $H^{1}$ method, while the method with $\lambda=1$ is motivated by the fact that the solution $\operatorname{\boldsymbol{\theta}}$ satisfies the divergence-free property, and thus the preservation of volume.

Since we consider domains which are $C^{1,1}$ -regular, we can solve the mean curvature as $\kappa=\operatorname{\mathrm{div}}_{\Gamma_{\rm f}}{\operatorname{\boldsymbol{n}}}=\operatorname{\mathrm{div}}{\operatorname{\boldsymbol{N}}}$ [17, Proposition 5.4.8], where $\operatorname{\boldsymbol{N}}$ is any extension of $\operatorname{\boldsymbol{n}}$ in $\Omega$ . Having this in mind, we solve for $\operatorname{\boldsymbol{N}}\in H^{1}(\Omega)^{2}$ by means of the equation

\displaystyle\varepsilon a({\operatorname{\boldsymbol{N}}},{\operatorname{\boldsymbol{\varphi}}})_{\Omega}+({\operatorname{\boldsymbol{N}}},{\operatorname{\boldsymbol{\varphi}}})_{\Gamma_{\rm f}}=({\operatorname{\boldsymbol{n}}},{\operatorname{\boldsymbol{\varphi}}})_{\Gamma_{\rm f}}\ for\ all\ \operatorname{\boldsymbol{\varphi}}\in H^{1}(\Omega)^{2}.

(21)

for sufficiently small $\varepsilon>0.$

4.3. Deformation Algorithm

Given an initial domain $\Omega_{0}$ , we shall solve the minimization problem iteratively. In particular, we shall update an iterate $\Omega_{k}$ by the identity perturbation operator, i.e., $\Omega_{k+1}=(I+t_{k}{\operatorname{\boldsymbol{\theta}}}_{k})\Omega_{k}$ . Here, ${\operatorname{\boldsymbol{\theta}}}_{k}$ is solved using (20) but on $\Omega_{k}$ instead of $\Omega$ . On the other hand, we shall solve for $t_{k}$ using a line-search method starting at an initial guess

t_{k}^{0}=\beta\frac{\mathcal{J}(\Omega_{k})}{\|\operatorname{\boldsymbol{\theta}}_{k}\|_{L^{2}(\Gamma_{\rm f})^{2}}^{2}}\text{ or }\beta\frac{\mathcal{L}(\Omega_{k})}{\|\operatorname{\boldsymbol{\theta}}_{k}\|_{L^{2}(\Gamma_{\rm f})^{2}}^{2}}

(22)

where $\beta>0$ is a fixed, sufficiently small parameter. This step size is updated whenever $\mathcal{J}(\Omega_{k+1})>\mathcal{J}(\Omega_{k})$ (or $\mathcal{L}(\Omega_{k+1})>\mathcal{L}(\Omega_{k})$ ) or if the deformed domain $\Omega_{k+1}$ exhibits mesh degeneracy (see [5, Section 3.4] for more details).

As for the augmented Lagrangian $\mathcal{L}$ , the multipliers $\ell$ and $b$ are iteratively updated as well. Inspired by Framework 17.3 in [28] for the $\ell$ -update, and by [5, Section 3.3], we use the following rule

\displaystyle\ell_{n+1}=\ell_{n}-b_{n}\mathcal{F}(\Omega_{k}),\text{ and }b_{n+1}=\left\{\begin{aligned} \tau b_{n}&\text{ if }b_{n}<\overline{b}\\ b_{n}&\text{ otherwise,}\end{aligned}\right.

(23)

where $\overline{b}>0$ is a given target parameter and $\tau>1$ is an increase-iteration parameter. One may refer to [5, Section 3.3] for a more detailed rationale on the mentioned updates.

In summary, we have the following algorithm:

Initialize:

\Omega_{0}

\beta

\tau

\gamma

\ell_{0}

, and

b_{0}

;

Determine the mean curvature

\kappa=\operatorname{\mathrm{div}}\operatorname{\boldsymbol{N}}

using (21);

Solve the respective states

({\operatorname{\boldsymbol{u}}}_{0},p_{0})

and

({\operatorname{\boldsymbol{v}}}_{0},\pi_{0})

from (14) and (18) on

\Omega_{0}

;

for k=1,…,max iteration do

Solve the value of the objective function

\mathcal{J}(\Omega_{k})

( or the Lagrangian

\mathcal{L}(\Omega_{k})

);

Solve the the deformation field

{\operatorname{\boldsymbol{\theta}}}_{k}\,

from (20) on

\Omega_{k}

;

Set

\Omega_{k+1}=(Id+t_{k}{\operatorname{\boldsymbol{\theta}}}_{k}\,)\Omega_{k}

, given that

t_{k}

is obtained from (22) and if

\mathcal{J}(\Omega_{k+1})<\mathcal{J}(\Omega_{k})

( or

\mathcal{L}(\Omega_{k+1})<\mathcal{L}(\Omega_{k})

);

Determine new mean curvature

\kappa

;

Solve new solutions

(\operatorname{\boldsymbol{u}}_{k+1},p_{k+1})

and

(\operatorname{\boldsymbol{v}}_{k+1},q_{k+1})

;

if $|\mathcal{J}(\Omega_{k+1})-\mathcal{J}(\Omega_{k+1})\text{( or }|\mathcal{L}(\Omega_{k+1})-\mathcal{L}(\Omega_{k+1})|\text{)}<$ tolerance then

break;

else

continue;

Update

\ell_{k+1}

and

b_{k+1}

using (23);

Algorithm 1 aL-algorithm if

\lambda=0

, dF-algorithm if

\lambda=1

4.4. Numerical Examples

Using the algorithms discussed previously, we will now show some numerically implemented examples, which are all run using the programming software FreeFem++ [16] on a Intel Core i7 CPU $@$ 3.80 GHz with 64GB RAM.

In all these examples, we shall consider a Poiseuille-like input function given by $\operatorname{\boldsymbol{g}}=(1.2(0.25-x_{2}^{2}),0)^{\top}$ . Furthermore, we consider the rectangular channel whose corners are $(0,-0.5)$ , $(0,0.5)$ , $(2,0.5)$ , and $(2,-0.5)$ . On the other hand we consider a circular initial shape for the free-boundary $\Gamma_{\rm f}$ that is centered at $(0.325,0)$ and with radius $r=0.13$ (see Figure 2). For simplicity, we shall also consider a fairly viscous fluid, i.e., $\nu=1/100.$

Refer to caption — Figure 2. Initial geometry of the domain with the refined mesh.

The volume constant $m$ is set as the volume of the initial domain $\Omega_{0}$ , which is refined according to the magnitude of the state solution with the minimum mesh size $hmin=1/50$ and a maximum mesh size $hmax=1/30$ . All the governing equations – the state equation, adjoint equation, state responsible for the deformation fields, and the smooth approximation of the normal vector – are solved using the UMFPACK option. The tolerance, furthermore, is selected to be $tol=1\times 10^{-6}.$

4.4.1. Augmented Lagrangian Method

The implementation of aL-algorithm is divided into two parts. Namely, we consider when $\gamma_{1}=1$ and $\gamma_{2}=0$ (which we call as the curl_aL-problem), $\gamma_{1}=0$ and $\gamma_{2}=1$ (which we call as the detgrad_aL-problem).

The values of the constants $\alpha$ , $\ell_{0}$ , $b_{0}$ , $\tau$ , and $\overline{b}$ are chosen so as to maintain stability of the method.

Implementation of curl_aL-problem ( $\gamma_{1}=1,\gamma_{2}=0$ ). For this implementation, we have chosen the numerical parameters shown in Table 1.

Table 1. Parameter values for curl_aL-problem

Parameter	Value	Parameter	Value
$\alpha$	6.0	$\ell_{0}$	20
$b_{0}$	$1\times 10^{-4}$	$\tau$	1.05
$\overline{b}$	10

Using these parameters, the evolution of the free-boundary is shown in Figure 3(left) below. We see that the free-boundary evolves into a bean-shaped surface.

As expected, we see that the value of the objective functional decreases on each iteration (see Figure 3(upper right)). In fact, the value of the objective functional decreases by $10.47\%$ . Lastly, we see from Figure 3(lower right) that the volume increases by $0.58\%$ .

Implementation of detgrad_aL-problem ( $\gamma_{1}=0,\gamma_{2}=1$ ). In this particular scenario, we consider the parameter values shown in Table 2.

Table 2. Parameter values for detgrad_aL-problem

Parameter	Value	Parameter	Value
$\alpha$	1.3	$\ell_{0}$	.5
$b_{0}$	$1\times 10^{-2}$	$\tau$	1.05
$\overline{b}$	10

Unlike the curl_aL-problem, the shape does not evolve into a bean-shaped obstacle. The evolution acts in a non-symmetrical manner (see Figure 4(left)) due to the fact that the right hand-side of the adjoint equation with the current choice of $\gamma_{1}$ and $\gamma_{2}$ is also non-symmetric. Even though this problem runs seven iterations, it can be observed that after two iterations, the shapes are almost identical as shown in Figure 4(left). This phenomena can also be observed on the trend of the decrease of the value of the objective functional, i.e., the values of the functional from the second iteration up to the last iteration differ minimally.

The percentage change from the initial shape to the second iterate-shape is $3.67\%$ , while from the second iterate to the final shape it is $0.03\%$ . In general, from the initial shape to the final shape the percentage difference is $3.7\%.$ We also point out that the minute difference of iterate 2 and the final shape may be caused by the fact that the changes are caused by mesh adaptation used in our program.

4.4.2. Divergence-free deformation fields

Unlike the implementation of aL-algorithm, the minimization problem is solved by dF-algorithm in three parts. In particular, we shall illustrate the following modificatoins:

$\bullet$

curl_dF-problem with parameters $\gamma_{1}=1$ , and $\gamma_{2}=0$ ;
$\bullet$

detgrad_dF-problem with parameters $\gamma_{1}=0$ , and $\gamma_{2}=1$ ;
$\bullet$

mixed_dF-problem with parameters $\gamma_{1},\gamma_{2}\geq 1$ .

Implementation of curl_dF-problem ( $\gamma_{1}=1,\gamma_{2}=0$ ). Note that in using dF-algorithm, we are more relaxed in the choice of the parameters. In particular, we are just concerned with the parameter $\alpha$ which we choose as $\alpha=5$ .

As we can see in Figure 5(left) the method converges to a bean shaped boundary, just like the generated shape in curl_aL-problem. However, in the current case, the method steers the shape to have the inner domain bounded by $\Gamma_{\rm f}$ lose its convexity on the left side.

One can also infer from Figure 5(upper right) that the relative difference - with respect to the initial shape - of the value of the objective function at the optimal shape is $1.85\%$ , while of the volume is $0.12\%$ .

If we compare this algorithm with aL-algorithm with $\alpha=5$ and $\ell_{0}=13,13.5,14$ , we deduce from Figure 6(lower right) that our current method preserves the volume better. We also see from Figure 6(upper right)) that the objective function values using dF-algorithm are lower than the ones from aL-algorithm except from the case $\ell_{0}=14$ , which fails horribly with volume preservation.

We also see in Figure 6(left) what we meant we we said that the inner domain bounded by $\Gamma_{\rm f}$ loses its convexity on the right side, i.e., the curvature of the shapes generated by aL-algorithm are more pronounced compared to our current algorithm.

Remark 6.

(i.) The reason why we started the simulation of aL-algorithm with $\alpha=5$ at $\ell_{0}=13$ is that if we try a smaller value of $\ell_{0}$ , the scheme becomes too stringent and will not generate any domain perturbation.

(ii.) One interesting feature in Figure 6 is that even after ten iterations, the objective function and volume trends for the case $\ell_{0}=14$ seem to have not converged yet. In fact, even after 20 iterations, the scheme still finds a descent direction. This highlights an issue with the balance between regularization and stability, which is an interesting venture to delve into if one is keen on studying stability of numerical methods. One may also refer to [27] for some insights on the issue of regularity and stability in general.

Implementation of detgrad_dF-problem ( $\gamma_{1}=0,\gamma_{2}=1$ ). Similar with our previous implementation, we only needed to choose a value for the parameter $\alpha$ , which in this case is $\alpha=1.0.$

As illustrated in Figure 7(left), the final perturbed domain differs minimally from the initial domain. Furthermore, we can see from Figure 8(left) that the deformation is less extensive as compared to the final shapes generated by the aL-algorithm. One can also infer from Figure 7(upper and lower right) that the value of the objective function and of the volume at the final shape differs relatively (w.r.t. the initial shape) by $1.30\%$ and $0.03\%$ , respectively. This implies that the volume preservation is significantly satisfied by the dF-algorithm, while maintaining a significant decrease on the objective function.

If we compare the current algorithm to aL-algorithm with the same value of $\alpha$ , we see from Figure 8(lower right) that the only value of $\ell_{0}$ that can preserve the initial volume better is $\ell_{0}=0.1$ . We also see from Figure 8(upper right) that since the difference between the final value of the objective function at the final shape using dF-algorithm and aL-algorithm with $\ell_{0}=0.1$ is very minimal, the aL-algorithm is more preferable. The only drawback though is that aL-algorithm depends on more parameters, hence it is more subject to instabilities.

The implementations of the aL-algorithm with $\ell_{0}=0.2,0.4$ exhibit a significant decrease on the objective functional as compared to the dF-algorithm, however they both violate the volume constraint significantly.

To illustrate another case where the aL-algorithm is more preferable for the minimization of the detgrad functional, Figure 9 shows an implementation where $\alpha=1.2$ , and $\ell_{0}=0.1$ for the aL-algorithm.

As can be easily observed in Figure 9(right), aL-algorithm yields better result when it comes to the minimization of the objective functional and preservation of the volume. In particular, the objective functional values are much lower for aL-algorithm, and the volume difference between the initial and final shape are much closer compared to the values generated by dF-algorithm.

Implementation of mixed_dF-problem ( $\gamma_{1},\gamma_{2}\geq 1$ ). In this implementation, we consider several configurations for reasons we shall explain later. In particular, we shall consider the value of parameters as shown in Table 3.

Table 3. Parameter Values for mixed_dF-problem

Configuration	$\alpha$	$\gamma_{1}$	$\gamma_{2}$
1	6.0	1.0	1.0
2	7.0	1.0	2.0
3	8.0	1.0	3.0
4	9.0	1.0	4.0
5	10.0	1.0	5.0
6	11.0	1.0	6.0
7	12.0	1.0	7.0
8	13.0	1.0	8.0
9	14.0	1.0	9.0
10	15.0	1.0	10.0

For the first configuration, Figure 10 shows how the initial shape evolves into a bean shaped boundary similar with that of curl_dF-problem. Furthermore, the value of the objective functional at the final shape differs relatively from the initial shape by $0.87\%$ while the volume is preserved with a relative difference of $0.30\%$ .

We note that in this configuration the final shape is almost the same with the final shape generated by curl_dF-problem, as shown in Figure 11(left).

To understand why this phenomena occurs, we note that the objective function in our configurations can be written as

\mathcal{J}(\Omega)=\left[5\int_{\Gamma_{\rm f}}\operatorname*{d\!}s-\frac{1}{2}\int_{\Omega}|\nabla\times{\operatorname{\boldsymbol{u}}}|\operatorname*{d\!}x\right]+k\left[\int_{\Gamma_{\rm f}}\operatorname*{d\!}s-\int_{\Omega}h(\mathrm{det}(\nabla{\operatorname{\boldsymbol{u}}}))\operatorname*{d\!}x\right],

where $k$ is the configuration number. We shall call the first bracket as the curl part of the mixed_df-problem while the other one the detgrad part. Evaluating the two parts numerically at the initial shape, we can see in Figure 11(right) that the value of the curl part is higher than the detgrad part, with the respective values $2.45$ and $0.65$ . This results to having the optimization problem prioritize minimizing the curl part of the mixed_dF-problem, and this results into having the similar shape as that of the curl_dF-problem. Intriguingly though, when compared to the values of the objective functional for the curl_dF-problem and the detgrad_dF-problem, the values of the curl part on each iterate is much higher than the objective function for curl_dF-problem (see Figure 11(upper right)) while the values of the detgrad part is lower than the objective function for detgrad_dF-problem (see Figure 11(lower right)).

To investigate the influence of the curl_dF-problem and the detgrad_dF-problem to the solution obtained by mixed_dF-problem even further, we implement configuration 2. In this example, the value of the detgrad part is higher than the one considered on the previous configuration.

As can be observed on Figure 12(left) the final shape for the mixed_dF-problem is now more distinct from that of the curl_dF-problem. In particular, we now see a bump emerging on the right side of the final shape for mixed_dF-problem. The emergence of such bump may be attributed to its presence on the final shape of detgrad_dF-problem. In fact, the value of the detgrad part is now higher compared to the previous configuration (see Figure 12(upper right) and (lower right)).

Skipping to the fourth configuration, we see that value of the detgrad part is higher than the curl part. From Figure 13(left), the final shape of the mixed_df-problem is now greatly influenced by the shape of detgrad_df-problem. In particular, it can be observed that the right side of the shape generated by mixed_df-problem converges to the shape induced by the detgrad_df-problem. However, we can still observe the effect of the shape by the curl_df-problem on the left side.

To see the convergence more clearly we utilize Hausdorff distance. Figure 14(upper right) illustrates the increase of the Hausdorff distance between the solutions of the mixed_dF-problem and the curl_dF-problem. Meanwhile, Figure 14(lower right) shows a decreasing trend in the Hausdorff distance between the solutions of the mixed_dF-problem and the detgrad_dF-problem.

4.5. Effects of the Shape Solutions to Generation of Twin-Vortex

We conclude this section by briefly illustrating the effects of the final shapes to twin-vortices induced by a time-dependent Navier–Stokes equations right before the shedding of Karman vortex, and is solved using a stabilized Lagrange–Galerkin method [29]. Since similar behaviors can be observed if we use aL-algorithm or dF-algorithm, we only illustrate the cases using dF-algorithm.

Figure 15 shows the comparison of the twin-vortices using the initial domain and the domain obtained from curl_dF-problem. The figure shows that the solution to the curl_dF-problem generates a longer vortex compared to that of the initial domain.

Similarly, the solution of the detgrad_dF-problem generates a longer vortex as compared to the initial domain as shown in Figure 15.

Lastly, Figure 17 shows that the vortex generated using the final shape of curl_df-problem is longer than that of the detgrad_df-problem.

Remark 7.

As mentioned before, the generation of vortices gives a possible application of our shape optimization problem to a branch of machine learning known as physical reservoir computing. According to Goto, K., et.al.[13], the length or the long diameter of the twin-vortices has a correlation with the memory capacity - as defined in the aforementioned reference - of the physical reservoir computer. With our previous observations, we can naively conjecture that the shapes generated by the implementations above will cause a better memory capacity for the physical reservoir computer. However, since these are just visual observations, we need more experiments to verify that these shapes indeed generate a good computational ability.

5. Conclusion

We presented a minimization problem that is intended to maximize the vorticity of the fluid governed by the Stokes equations. Furthermore, we considered a Tikhonov-type regularization - in the form of the perimeter functional - and a volume constraint. The shape sensitivity analysis is carried out by utilizing the rearrangement method which gave us the necessary conditions in the form of the Zolesio–Hadamard Structure. We then implemented the necessary conditions in two ways, one is by utilizing the augmented Lagrangian, and by utilizing a new method for solving the deformation field that possesses the divergence-free condition. From these methods, we showed some numerical examples. Lastly, we illustrated how the shape solutions affect the flow, and we briefly mentioned a possible application of this problem in the field of physical reservoir computers.

We end this note by pointing out that the investigation in the stability of the numerical approaches with respect to the parameters were only scratched on the surface and can be explored with more rigor. Furthermore, the visual observations on the emergence of the twin-vortex are not enough to give a good conclusion with regards to the memory capacity of the physical reservoir computer. With this, one can investigate even further and try some more tests as carried out in [13].

References

[1] H. Azegami. Solution of shape optimization problem and its application to product design. In Hiromichi Itou, Masato Kimura, Vladimír Chalupecký, Kohji Ohtsuka, Daisuke Tagami, and Akira Takada, editors, Mathematical Analysis of Continuum Mechanics and Industrial Applications, pages 83–98, Singapore, 2017. Springer Singapore.
[2] J. Bacani and G. Peichl. On the first-order shape derivative of the Kohn–Vogelius cost functional of the Bernoulli problem. Abstract and Applied Analysis, 2013:1–19, 2013.
[3] D. Chenais. On the existence of a solution in a domain identification problem. Journal of Mathematical Analysis and Applications, 52(2):189 – 219, 1975.
[4] M. S. Chong, A. E. Perry, and B. J. Cantwell. A general classification of three‐dimensional flow fields. Physics of Fluids A: Fluid Dynamics, 2(5):765–777, 1990.
[5] C. Dapogny, P. Frey, F. Omnès, and Y. Privat. Geometrical shape optimization in fluid mechanics using freefem++. Structural and Multidisciplinary Optimization, 58(6):2761–2788, 2018.
[6] M. Delfour and J.-P. Zolesio. Shapes and Geometries: Metrics, Analysis, Differential Calculus, and Optimization. Society for Industrial and Applied Mathematics, 3600 Market Street, 6th Floor, Philadelphia, PA 19104-2688 USA, 2 edition, 2011.
[7] M. Desai and K. Ito. Optimal controls of Navier–Stokes equations. SIAM Journal on Control and Optimization, 32(5):1428–1446, 1994.
[8] M. F. Eggl and Peter J. Schmid. Mixing enhancement in binary fluids using optimised stirring strategies. Journal of Fluid Mechanics, 899:A24, 2020.
[9] L.C. Evans. Partial Differential Equations. Graduate studies in mathematics. American Mathematical Society, 1998.
[10] T. L. B. Flinois and T. Colonius. Optimal control of circular cylinder wakes using long control horizons. Physics of Fluids, 27(8):087105, 2015.
[11] Z. M. Gao, Y. C. Ma, and H. W. Zhuang. Shape optimization for navier–stokes flow. Inverse Problems in Science and Engineering, 16(5):583–616, 2008.
[12] V. Girault and P.-A. Raviart. Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms, volume 5 of Springer Series in Computational Mathematics. Springer, Berlin, Heidelberg, 1986.
[13] K. Goto, K. Nakajima, and H. Notsu. Twin vortex computer in fluid flow. New Journal of Physics, 23(6):063051, jun 2021.
[14] J. Haslinger, K. Ito, T. Kozubek, K. Kunisch, and G. Peichl. On the shape derivative for problems of Bernoulli type. Interfaces and Free Boundaries, 11(2):231–251, 2006.
[15] J. Haslinger, J. Málek, and J. Stebel. Shape optimization in problems governed by generalised navier-stokes equations: Existence analysis. Control and cybernetics, 34:283–303, 03 2005.
[16] F. Hecht. New development in freefem++. J. Numer. Math., 20(3-4):251–265, 2012.
[17] A. Henrot and M. Pierre. Shape Variation and Optimization: A Geometrical Analysis. EMS Tracts in Mathematics, 28 edition, 2014.
[18] A. Henrot and Y. Privat. What is the optimal shape of a pipe? Archive for Rational Mechanics and Analysis, 196(1):281–302, 2010.
[19] J. C. R. Hunt, A. A. Wray, and P. Moin. Eddies, streams, and convergence zones in turbulent flows. In Proc. 1988 Summer Program of Center for Turbulence Research Program, 1998.
[20] Ito, K., Kunisch, K., and Peichl, G. H. Variational approach to shape derivatives. ESAIM: COCV, 14(3):517–539, 2008.
[21] Y. Iwata, H. Azegami, T. Aoyama, and E. Katamine. Numerical solution to shape optimization problems for non-stationary navier-stokes problems. JSIAM Letters, 2:37–40, 2010.
[22] Jinhee Jeong and Fazle Hussain. On the identification of a vortex. Journal of Fluid Mechanics, 285:69–94, 1995.
[23] H. Kasumba and K. Kunisch. Vortex control in channel flows using translational invariant cost functionals. Computational Optimization and Applications, 52(3):691–717, 2012.
[24] G. Mather, I. Mezić, S. Grivopoulos, U. Vaidya, and L. Petzold. Optimal control of mixing in Stokes fluid flows. Journal of Fluid Mechanics, 580:261–281, 2007.
[25] B. Mohammadi and O. Pironneau. Applied Shape Optimization for Fluids. NUMERICAL MATHEMATICS AND SCIENTIFIC COMPUTATION. Oxford University Press, 2010.
[26] P. Monk. Finite Element Methods for Maxwell’s Equations. Oxford University Press, 2003.
[27] M. T. Nair, M. Hegland, and R. S. Anderssen. The trade-off between regularity and stability in tikhonov regularization. Mathematics of Computation, 66(217):193–206, 1997.
[28] J. Nocedal and S. Wright. Numerical Optimization, volume 2 of Springer Series in Operations Research and Financial Engineering. Springer-Verlag New York, 2006.
[29] Notsu, H. and Tabata, M. Error estimates of a stabilized lagrange-galerkin scheme for the navier-stokes equations. ESAIM: M2AN, 50(2):361–380, 2016.
[30] M. Pošta and T. Roubíček. Optimal control of Navier–Stokes equations by Oseen approximation. Computers & Mathematics with Applications, 53(3):569–581, 2007. Recent Advances in the Mathematical Analysis of Non-Linear Phenomena.
[31] J. F. T. Rabago and H. Azegami. An improved shape optimization formulation of the bernoulli problem by tracking the neumann data. Journal of Engineering Mathematics, 117(1):1–29, 2019.
[32] J. F. T. Rabago and H. Azegami. A second-order shape optimization algorithm for solving the exterior Bernoulli free boundary problem using a new boundary cost functional. Computational Optimization and Applications, 2020.
[33] S. Schmidt and V. Schulz. Shape derivatives for general objective functions and the incompressible navier-stokes equations. Control and Cybernetics, 39(3):677–713, 2010.
[34] J. Sokolowski and J.-P. Zolesio. Introduction to Shape Optimization: Shape Sensitivity Analysis. Springer Series in Computational Mathematics. Springer-Verlag Berlin Heidelberg, 1 edition, 1992.