A shape design problem for the Navier–Stokes flow with a convective boundary condition

In this study, a shape optimization problem for the two-dimensional stationary Navier–Stokes equations with an artificial boundary condition is considered. The fluid is assumed to be flowing through a rectangular channel, and the artificial boundary condition is formulated so as to take into account the possibility of ill-posedness caused by the usual do-nothing boundary condition. The goal of the optimization problem is to maximize the vorticity of the said fluid by determining the shape of an obstacle inside the channel. Meanwhile, the shape variation is limited by a perimeter functional and a volume constraint. The perimeter functional was considered to act as a Tikhonov regularizer and the volume constraint is added to exempt us from possible topological changes. The shape derivative of the objective functional was formulated using the rearrangement method, and this derivative was later on used for gradient descent methods. Additionally, an augmented Lagrangian method and a class of solenoidal deformation fields were considered to take into account the goal of volume preservation. Lastly, numerical examples based on gradient descent and the volume preservation methods are presented.


Introduction
Partial differential equation (PDE) constrained optimization is among the most active fields in mathematics, mainly due to its applicability in engineering and physics. In particular, the system of the Navier-Stokes equations is the theme of a lot of papers because of its ability to closely mimic physical fluid flow. From this reason, several authors including Colin and Fabrie (2010), Desai andIto (1994), andFürsikov et al. (2000) have studied optimal control problems where the controls are designed to steer the fluid dynamics according to a given objective.
In this paper, we are interested in an optimization problem where the control is the shape of the domain instead of functionals. Several authors have considered such problems due to its applicability for example in aeronautics (Mohammadi and Pironneau 2004). Moubachir and Zolesio (2006) wrote an extensive introductory literature for shape optimization problems which are governed by fluid flows and are mostly described through the Navier-Stokes equations, and the objectives are mostly formulated as tracking functionals. For a more applied approach to fluid shape design problems, we refer the reader to Mohammadi and Pironneau (2010). Meanwhile, Kasumba and Kunisch (2012) analyzed a vorticity minimization problem that steers the fluid flowing through a channel to exhibit laminar flow. In the said literature, the authors considered a bounded domain, which compels the imposition of an input function on one end of the channel, and a do-nothing boundary condition on the other to capture outflow. However, due to the nonlinear nature of the Navier-Stokes equations, the outflow condition the authors considered can cause nonexistence of weak solutions. For this reason, we propose an artificial boundary condition that can capture the outflow while making sure of the well-posedness of the state equations.
We also mention that although most shape design problems involving fluid flow are formulated as to minimize vorticity or drag, there are literatures where vorticity maximization is the goal. This type of maximization problems are mostly considered applicable, for example, in optimal mixing problems (Eggl and Schmid 2020;Mather et al. 2007). Recently, Goto et al. (2021) also showed that emergence of vortices has an interesting result in the field of information theory.
With all these, we consider the following shape optimization problem min{G(u, ) : | | = m, ⊂ D}, where D is a hold-all domain which is assumed to be a fixed bounded connected domain in R 2 , ⊂ D is an open bounded domain, G(u, ) := J (u, ) + α P( ), J and P are vortex and perimeter functionals, respectively, given by J (u, ) = − γ 2 |∇ × u| 2 dx, and P( ) = f ds, m > 0 is a given constant, and | | := dx. Here, u is the velocity field governed by a fluid flowing through a channel with an obstacle. The flow of the fluid is reinforced by an input function g on the left end of the channel, which is denoted by in , whilst an outflow boundary condition is imposed to the fluid on the right end of the channel denoted by out . The boundary f is the boundary of the submerged obstacle in the fluid. The remaining boundaries of the channel are denoted by w , upon which-together with the obstacle boundary f -a no-slip boundary condition is imposed on the fluid. For the condition on out , let us point out that the usual candidate for capturing outflow is the so-called do-nothing condition, however such condition is not sufficient to ensure the existence of solutions to the Navier-Stokes equations. Aside from the do-nothing condition, one formulation caught our attention. An artificial boundary condition called the directional do-nothing shows to be a good candidate for imposing outflow while ensuring well-posedness of the Navier-Stokes system. However, this condition generates complications for optimization problems, especially when an adjoint approach is considered. In this paper, we decided to consider a boundary condition, which we call a convective boundary condition. The said conditions will be discussed even further in the coming section.
From here, when we talk about a domain , we always consider it having the boundary ∂ = in ∪ w ∪ out ∪ f , and that dist(∂ \ f , f ) > 0.
Due to the non-convexity of the vorticity functional J , possible non-existence of an optimal solution is a dilemma. Nevertheless, we shall show that J is sequentially continuous with respect to the state variables. Combine this with the continuity of the map → u with respect to some topology in the set of admissible domains, we circumvent the issue of existence of optimal shapes. We also propose a Tikhonov regularization in the form of the perimeter functional and a volume constraint to handle possible topological changes 1 .
To numerically solve the optimization problem, we shall utilize a gradient descent method based on the first order Eulerian derivative of the shape functional G. Although there are several methods to determine the derivative with respect to shape variations-for example, Gao et al. (2008) utilized minimax formulation and chain rule based on the Piola transform (see also the work of Schmidt and Schulz (2011) for computation of shape derivatives to general functionals)-the formulation of the shape derivative will be aided by the so-called rearrangement method which was formalized by Ito et al. (2008). However, this computation will not take into account the volume constraint. This compels us to utilize two methods: first is the augmented Lagrangian method based on (Nocedal and Wright 2006, Framework17.3) which was applied to shape optimization problems by Dapogny et al. (2018); and the other one is by using solenoidal deformation fields. Comparing the solutions of the two methods, we recognize a common profile-although the shapes are virtually different-which is a bean-shaped obstacle.
This paper is organized as follows: in the next section, we introduce the governing state equations, and its corresponding variational formulation. In Sect. 3, we establish the existence of optimal shapes, where we shall use the L ∞ -topology on the set of characteristic functions. The sensitivity of the objective functional will be analyzed in Sect. 4. Section 5 is dedicated to the discussion of the numerical treatment by which we shall make use of the derivative formulated in Sect. 4. In this section, we shall also show the convergence of numerical solutions to a manufactured exact solution with respect to the Hausdorff measure, and the H 1 and L 2 norms. Concluding remarks and possible future works will be discussed in Sect. 6.

Let
⊂ D be an open and bounded domain, the motion of the fluid is described by the velocity u and the pressure p which satisfy the stationary Navier-Stokes equations given by: Here, g is a divergence-free input function acting on the boundary in . The condition u = 0 on f ∪ w corresponds to no-slip boundary condition. For boundary out , we choose an appropriate condition that will correspond to an outflow condition and that will give us a good energy estimate that is crucial in showing the existence of the velocity u and for showing continuity with respect to domain variations. Usually, the condition imposed to capture fluid outflow is what we call the do-nothing boundary condition. Mathematically, this condition is written by letting the product of the stress tensor and the vector normal to the boundary out equal to zero, i.e., − pn + ν∂ n u = 0, where ∂ n u := (∇u)n, and n is the outward unit vector normal to out . This condition has been considered as the outflow character in a lot of inquiries including that of Gresho (1991), among others. This profile, however, does not ensure the possibility of obtaining a good energy estimate, let alone the existence of a weak solution. As such, several alternative weak formulations has been formulated to address this issue (Bruneau and Fabrie 1996;Heywood et al. 1992;Zhou and Saito 2016). To illustrate the severity of the issue, we note that to solve the existence of a weak solution for the Navier-Stokes equations coupled with the do-nothing boundary condition the computations will give us the following expression out (u · n)|v| 2 ds. Now, if u and v are identical, which will be a case when trying to establish the existence, specifically on the coercivity of the left-hand side of the weak formulation, it is imperative to have a knowledge on the sign of the above quantity. This will be challenging due to its cubic form. This problem is addressed and circumvented by Bruneau and Fabrie (1996) by introducing several versions of artificial boundary conditions for the Neumann condition. In this paper, we shall utilize one of the conditions introduced in the said literature, as well as establish its well-posedness.
We denote by H r 0 ( ) the closure of W( ) with respect to the norm in H r ( ) 2 . By utilizing Meyers-Serrin type arguments it can be easily shown that To deal with the incompressibility condition, we shall take into account the following solenoidal spaces The space W( ) is dense in V( ), and the spaces V( ) and H( ) satisfy the Galfand triple property, i.e., V( ) → H( ) → V( ) * , where the first imbedding is dense, continuous and compact. Lastly, for simplicity of notations, we shall denote by (·, ·) D the L 2 (D), L 2 (D) 2 , or L 2 (D) 2×2 inner products, where D is any measurable set.

Weak formulation and existence of solutions
Before we begin, let us first settle the issue on the regularity of the domain. Although the assumption that is of class C 1 is sufficient to ensure the existence of weak solution to the Navier-Stokes equations, which we shall show later, this resulting solution will be established to be at most first differentiable. With this reason, if we want a higher regularity for the solution then additional regularity on the domain is needed to be imposed. So from here on out, we shall say that the domain satisfies (H ) if either of the following assumptions is satisfied: • the outer surface is a boundary of a convex polygon and the inner surface is of class C 1,1 .

(H )
The reason for mentioning the above regularity assumptions on is twofold: one-as referred to before-is for the regularity of the weak solution; and the other reason is to ensure that the extension of the input function g that satisfies the properties below exists ⎧ ⎪ ⎨ ⎪ ⎩ div g = 0 in , g = 0 on f ∪ w , and out g · n ds = − in g · n ds ≥ 0.
Of course, if we assume minimum regularity on the domain, e.g., is of class C 1 , then we can easily show, by the virtue of (Girault and Raviart 1986, Lemma 2.2), that the extension of the input function g ∈ H 1/2 ( in ) 2 that satisfies (2) exists and is in H 1 ( ) 2 . However, this is not enough if we would wish to have higher regularity on the extension. To be precise, if we wish to extend g ∈ H 3/2 ( in ) 2 in that satisfies (2) and that the extension is in H 2 ( ) 2 , then the minimum regularity is for the domain to satisfy (H ). To simplify things later, we shall use the same notation for the input function and its extension. That is, when we refer to the assumption that g ∈ H m ( ) 2 and satisfies (2), we are implicitly saying that the input function is in H m−1/2 ( in ) 2 . By lettingũ = u − g, we consider the following variational problem: For a given domain , findũ ∈ V( ) such that forms, respectively, defined as follows: where the members of A are defined below: for all ϕ ∈ V( ). Note that the variational equation (3) is achieved from assuming that u satisfies the boundary condition We further mention that this boundary equation is among the artificial conditions proposed by Bruneau and Fabrie (1996) for (a) = a. As far as we are aware, this condition has never been analyzed, let alone be considered in an optimization problem. The problem arises in-as we have mentioned before-showing the coercivity of the operators of the left hand side of the variational equation and due to the non-homogeneous Dirichlet data on in . More obvious formulations have been considered by several authors (Zhou and Saito 2016;Braack and Mucha 2014) by considering the nonlinear part to be (u · n) − u, i.e., where the notation (·) − is defined as This will automatically give a coercive left hand side, and hence the well-posedness of the state equations. However, this formulation will not be easy to handle when a variational approach for the optimization problem is utilized, and in the numerical treatment of the adjoint problem, to be specific. Fortunately, we shall prove an analogous result to that of (Girault and Raviart 1986, Lemma IV.2.3) but takes into account the boundary integral on out that will help to render the variational problem (3) well-posed. The mentioned analogy is shown in Lemma 4 and is proven in the Appendix.
We start the analysis by showing-under appropriate conditions-that ∈ V( ) * .
Proof The proof will utilize the compact embeddings H 1 ( ) 2 → L q (∂ ) 2 and H 1 ( ) 2 → L q ( ) 2 for q ≥ 2 (see for example (Adams and Fournier 2003, Part I, Theorem 6 where c = max{1, ν,c 1 /2} andc 1 is dependent on the L ∞ -norm of the outward normal vector n which is bounded due to the regularity assumptions on the domain. Lastly, the linearity follows from the linearity of the action itself. Any functionũ ∈ V( ) that solves the variational equation (3) is said to be a weak solution to the Navier-Stokes equations (SE) with the boundary condition (BC1). The existence of the solutionũ is summarized below.
The proof of the theorem will make use of the following lemmata whose proofs can be easily redone and can be based for example on (Temam 2001, Lemma II.1.8) and (Girault and Raviart 1986, Lemma IV.2.3).
Lemma 3 Let u, v, w ∈ V( ), the trilinear form b satisfies the following properties.
Lemma 4 For any γ > 0 there exists w 0 = w 0 (γ ) ∈ H 1 ( ) 2 such that Remark 1 We note that the assumption that the input function g is in H 1 ( ) 2 is valid due to Lemma 4. Meaning to say, w 0 = g not only on 0 but also inside . In particular, we also get the estimate Furthermore, let us also point out that even though the input function is a fixed Dirichlet profile on the boundary in , its extension is dependent on the domain . This implies that the extension-which we chose to denote similarly as g-in will also be sensitive with domain deformations. Aside from the lemmata, it is noteworthy to point out the well definedness of the trilinear form b, i.e., [(u · ∇)v] · w ∈ L 1 ( ) whenever u, v, w ∈ V( ).

Proof of Theorem 2
The proof utilizes Galerkin method by considering the finite dimensional subspaces V n ⊂ V( ). We consider the projected problem on V n , i.e., we find solutions u n ∈ V n that satisfies for all ϕ ∈ V n the equation To show the existence of solutionsũ n ∈ V n , we shall employ (Girault and Raviart 1986, Corollary IV.1.1), and assume for now the coercivity of B and its sequential weak continuity, which we shall prove later. Let us begin by introducing n : V n → V n defined by ( n (u), ϕ) = B(u, ϕ) − , ϕ V( ) * ×V( ) for any ϕ ∈ V n , where the notation (·, ·) is the inner-product in V n . From this operator, we notice thatũ n ∈ V n is a solution to (8) if and only if n (ũ n ) = 0.
By the sequential weak continuity of B we get that n is continuous. Since V n is finite dimensional, we can use (Girault and Raviart 1986, Corollary IV.1.1) if we can find a constant α > 0 such that ( n (u), u) ≥ 0 whenever u ∈ V n satisfies u V ( ) = α.
We focus on the first convergence, since the others can be done similarly. Before we begin, we note that H 1 ( ) is compactly embedded to L q (∂ ), for q ≥ 2 (see for example (Adams and Fournier 2003, Part I, Theorem 6.3)). This implies thatũ n →ũ andũ n · n →ũ · n in L q ( out ) 2 and L q ( out ), respectively. Hence, we have the following computation: The right side approaches zero from the mentioned convergences in L q ( out ) 2 and L q ( out ), which proves the claim. Therefore,ũ solves (3) with the energy estimate Remark 2 (i) Even though we just assumed f ∈ L 2 ( ) 2 and g ∈ H 1 ( ) 2 , as we shall see later, we can have these functions extended to the hold all domain. Furthermore, since ⊂ D, the energy estimate can be extended to the hold-all domain D, i.e., and c > 0 is dependent on D, thanks to Faber-Krahn inequality. (ii) As previously mentioned, it can be shown that a more regular domain yields a more regular solution. In particular, if satisfies (H ), then the weak solution to (3) satisfiesũ ∈ V( ) ∩ H 2 ( ) 2 , and as a consequence of the Rellich-Kondrachov embedding theorem (Evans 1998, Chapter 5, Theorem 6) the solution is in C( ) 2 .
The existence of the distribution p is a direct consequence of De Rham's theorem (Temam 2001, Proposition I.1.1) in such a way that p ∈ L 2 ( ). Furthermore, if we assume that then one can show that the solution is unique.
In strong form, if we assume appropriate regularity on , additionally if g ∈ H 2 (D) 2 (a property which we shall implicitly assume from now on), andũ ∈ V( ) ∩ H 2 ( ) 2 solves the variational equation (3), then u =ũ + g ∈ H 2 ( ) 2 is regarded as the solution to the equations, Before closing this section, we look at some operators which will be handy when we go to the investigation of sensitivities. First, let us introduce the Stokes operator A : D(A ) ⊂ H( ) → H( ) defined by A u = −P u for u ∈ D(A ). Here P : L 2 ( ) 2 → H( ) is the Leray projection that is associated with the decomposition L 2 ( ) 2 = H( ) ⊕ ∇ L 2 ( ). As a direct consequence, if the domain satisfies (H ) then D(A ) = V( ) ∩ H 2 ( ) 2 , hence item (ii) of the remark above. We mention, furthermore, that the Stokes operator is a self-adjoint positive operator with dense domain and compact inverse, this in turn, gives us the orthonormal basis of H( ). This orthonormal basis makes the Galerkin method in the proof of Theorem 2 possible to utilize. With this operator, we can write the first member of A as follows, We also introduce the operator B : For any u ∈ V( ), we shall also use the notation B u := B (u, u). Lastly, we consider the boundary operator C : The well-definedness and continuity of the operator C can be established by utilizing the Rellich-Kondrachov embedding (Adams and Fournier 2003, Part I, Theorem 6.3).

Existence of optimal shapes
The purpose of this section is to establish the existence of a solution to the shape optimization problem. To be able to do this, a set of admissible domains O ad will be considered to ensure the existence of state solutions as well as ensure that-assuming that it exists-the solution is embedded to the hold-all domain D. In this set of admissible domains, a topology will be endowed upon which O ad itself is compact. Utilizing this compactness property, standard sequential arguments will then be followed to establish the promised existence of the optimal shape.

Cone property and the set of admissible domains
Before we begin, we assume that out ⊂ ∂ D, which is not going to be an issue on the generality of the problem since out is not a part of the free-boundary. Now, to define the topology on the set of admissible domains, we shall resort to the collection of domains that satisfy the cone property. Note that one way to define the topology on the set of admissible domains is by parametrizing the free-boundary. Examples of such parametrizations are found in Azegami (2019, 2020) and the references therein. However defining the domains in such way will be problematic for the current problem since the free-boundary parametrization might lead to generating domains with varying volumes.
To pass through the challenge of preserving the volume, a classical way to define the topology on the set of admissible domains is by considering the collection of characteristic functions, which is a fact highlighted by Henrot and Privat (2010) and is rigorously discussed by Delfour and Zolesio (2011) and Henrot and Pierre (2014). But before we delve into this topology, let us first look at the manner by which the domains are defined. In particular, we shall consider domains which satisfy the cone property as defined by Chenais (1975).
(i.) The cone of angle θ , height h and axis ξ is defined as where (·, ·) and · denote the inner product and Euclidean norm in R 2 , respectively. (ii.) A set ⊂ R 2 is said to satisfy the cone property if and only if for all x ∈ ∂ , there From this definition, we now define the set of admissible domains as O ad := { ⊂ D : satisfies the cone propert y and | | = m}.
where the function χ A for a set A ⊂ R 2 refers to the characteristic function defined by Remark 3 (i.) This set of admissible domains has been established to be non-empty (see the proof of Proposition 4.1.1 by Henrot and Pierre (2014)), which exempts us from the futility of the analyses we will be going through in the succeeding sections. (ii.) Note that this convergence may seem too weak in the sense that the weak * convergence only assures us that the limit χ only satisfies χ (x) ∈ [0, 1], however as pointed out by Henrot and Pierre (2014) in Proposition 2.2.1, this convergence holds in the space L p loc (R 2 ) and thus χ is an almost everywhere characteristic function. We refer the reader to (Delfour and Zolesio 2011, Chapter 5) and (Henrot and Pierre 2014, Chapter 2 Sect. 3) for a more detailed discussion on the topology of characteristic functions of finite measurable domains.
Before we mention the compactness of the set O ad , let us first look at an important implication of the cone property, i.e., the existence of a uniform extension operator. This is given by the lemma below.
Lemma 5 (Chenais (1975)) There exists K > 0 such that for all ∈ O ad , there exists operators This lemma will be utilized in several occasions, for example when we show that the domain-to-state map is continuous. Meanwhile, the compactness of O ad follows from the fact that it is closed and relatively compact-as defined by Chenais (1975)-with respect to the weak * -L ∞ topology on U ad := {χ : ∈ O ad }. One can also read upon the proof in (Henrot and Pierre 2014, Proposition 2.4.10). We shall not discuss the proof of such properties, nevertheless they are summarized on the lemma below.

Lemma 6 The set O ad is compact with respect to the topology on U ad .
Aside from these properties, it is noteworthy to mention the set of admissible domains can be identified as Lipschitzian domains as well. Since we only consider the hold-all domain D to be possessing a bounded boundary, a proof of this property can be found in (Henrot and Pierre 2014, Theorem 2.4.7).

Well-posedness of the optimization problem
Now that we have defined a good topology on the set of admissible domains, this subsection is dedicated to establishing the existence of the solution to the shape optimization problem.
Recall that the functional J is written as a function of the state solution u and of the domain . Fortunately, we point out that for each ∈ O ad there exists a solution u ∈ V( ) to the weak formulation (3), hence the map → u. This implies that the well-posedness of the optimization will depend on the continuity of the domain-to-state map, which we shall briefly establish shortly. For now, let us consider the velocity-pressure formulation of (3) given by finding (ũ, (11) follows from the existence of solution to (3) and because d(·, ·) satisfies the inf-sup condition (Bertoluzza et al. 2017). Furthermore, the following energy estimate holds The impetus for introducing the variational equation (11) is to make sure that the divergence-free property of the states will be preserved on any domain in O ad and so that the said property will still be reflected when the uniform extension property in Lemma 5 is utilized.
Proof From the uniform extension property, there exists K > 0 such that Furthermore, from the energy estimate (12) From the uniform boundedness of {(u n , p n )} in H 1 (D) 2 × L 2 (D) and by the virtue of the Rellich-Kondrachov and Banach-Alaoglu theorems, there exists a subsequence of {(u n , p n )}, which we denote in the same manner, and an element (u, Passing through the limit. The next step is to show that (ũ, p) solves (11) in . By the assumed domain convergence, we have χ n := χ n * χ := χ in L ∞ (D).
Furthermore, since (ũ n , p n ) solves (11) in n , then for any Here, for any function ϑ : D → R, the trilinear form A ϑ (·, ·, ·) D can be dissected into several components, namely Our goal is to show that (ũ, p) = (u, p) is a solution to (11) by establishing that the following system holds by utilizing (13) and the weak- * limit of the characteristic functions on (14).
Using the same arguments as in Theorem 2, it can be easily shown that Furthermore, from the assumed convergence of the characteristic functions, we infer that Lastly, since u n u in H 1 (D) 2 and p n p in L 2 (D), then By letting (ϕ, q) ∈ H 1 0 ( ) 2 × L 2 ( ), and defining (ψ, φ) ∈ H 1 (D) 2 × L 2 (D) by (ψ, φ) = (ϕ, q) in and (ψ, φ) = (0, 0) in D\ , we get the variational equation (11) in . Before we end this proof, we note that on (16), a quite challenging and questionable inference is the convergence a χ n (u n , ψ) → a χ (u, ψ). The problem with this convergence is that as of now, we were only able to mention the weak convergence of the velocity fields in H 1 (D) 2 and the weak * convergence of the characteristic functions in L ∞ (D). Fortunately, we can actually show that χ n ∇u n converges strongly to χ∇u in L 2 (D) 2×2 , which establishes the convergence in question.
Indeed, by substituting (ψ, φ) = (u n , p) into equation (14), by recalling the strong convergence of u n to u in L 2 (D) 2 and the convergence χ n * χ in L ∞ (D), and since Here, we also used the implicit assumption that g ∈ H 2 (D) 2 ; hence, the convergence on the second line.
Just to reiterate, Proposition 7 proves that the map → (ũ, p) is continuous. This property will be instrumental to prove that the optimal shape exists. In fact, we shall use the fact that we can write the objective functional as a function that depends solely on the elements of O ad , and show that it is continuous with respect to this collection as well. To be precise, we prove the existence of a solution to the shape optimization problem on the theorem below. Before we begin, let us introduce the following notations which were made possible from the well definedness of the map → u: Theorem 8 Suppose that the assumptions for Theorem 2 together with the uniqueness assumption (9) hold; then, there exists * ∈ O ad such that Proof First, we show that the objective functional G is lower semicontinuous. This is done by showing that the component J is continuous (hence upper-semicontinuous) with respect to the state variable u =ũ + g, and using the lower-semicontinuity of the perimeter functional (Henrot and Pierre 2014, Proposition 2.3.7).
We shall then show that G is bounded from below, which will imply the existence of a minimizing sequence of domains whose evaluations will converge to an infimum value of the objective functional. Using the compactness of O ad , we shall then show that this sequence converges to a domain such that its evaluation coincides with the infimum.
Note that we can estimate J by the H 1 norm of the state u by virtue of the following computation: Furthermore, since Proposition 7 implies that the map →ũ + g is continuous, the map → J ( ) is also continuous, and hence upper-semicontinuous. Thus, the map → G( ) is lower-semicontinuous, i.e., for any sequence { n } n ⊂ O ad that converges to an element * ∈ O ad , then Step 2: Existence of a minimizing sequence. First, using estimate (17), we can show that J is uniformly bounded. Indeed, from (17) and the energy estimate (5), we get Furthermore, since any domain ∈ O ad is bounded, then the inner boundary f is also bounded and that its perimeter can be bounded below uniformly, say α P( ) ≥ P * for any ∈ O ad . Therefore, G is bounded from below, i.e., for any ∈ O ad Hence, there exists a sequence { n } n ⊂ O ad such that Step 3: Existence of a minimizer for G.
Since O ad is compact, the sequence { n } n from Step 2 has a subsequence-which we shall denote similarly-that converges to an element * ∈ O ad . Hence, from the lowersemicontinuity of G, we establish that * is the minimizer of G:

Shape sensitivity analysis
This section is dedicated to investigate the sensitivity of the objective functional G with respect to domain variations. We start this section by introducing the identity perturbation method, where we consider domain variations generated by a given autonomous velocity.

Identity perturbation
We consider a family of autonomous velocity fields θ belonging to : . We note that this operator is generated as a direct consequence of the velocity method discussed by Sokolowski and Zolesio (1992) and Henrot and Pierre (2014).
A given domain ⊂ D is perturbed by means of the identity perturbation operator so that for some t 0 := t 0 (θ ) > 0 we get the family of perturbed domains By the definition of θ , θ ≡ 0 on out , in , and w , this implies that these boundaries are part of the perturbed domains t . To be precise, we have Additionally, a domain that has at most C 1,1 regularity preserves its said regularity with this transformation, this means that for 0 ≤ t ≤ t 0 , t has C 1,1 regularity given that the initial domain is a C 1,1 domain.
Before we move further in this exposition, let us look at some vital properties of T t .
Lemma 9 (Sokolowski and Zolesio (1992); Delfour and Zolesio (2011)) Let θ ∈ , then for sufficiently small t 0 > 0, the identity perturbation operator T t satisfies the following properties: Let us recall Hadamard's identity which will be integral for solving the necessary conditions.
Proof See Theorem 5.2.2 and Proposition 5.4.4 of Henrot and Pierre (2014).

Rearrangement method
To investigate the sensitivity of the objective functional with respect to shape variations generated by the transformation T t , we shall resort to a variational approach formalized by Ito et al. (2008), which is known by many as the rearrangement method. This approach gets rid of the tedious process of solving first the sensitivity of the state solutions, then solving the shape derivative of the objective functional. Aside from the convenience the rearrangement method poses, we also mention that using the usual methods-such as the chain rule, minmax formulation, etc.-will not take into account the linearization of the state on the fixed boundary out . This in turn will render the linearized state and the adjoint equation ill-posed. This problem, thankfully, is resolved by the rearrangement method which is focused on the Frechét derivative of the state operator.
To start with, we consider a Hilbert space Y ( ) and an operator Suppose that the free boundary is denoted by f ⊂ ∂ , and g : Y ( ) → R, the said method deals with the shape optimization We define the Eulerian derivative of J at in the direction θ ∈ by where y t solves the equation E(y t , t ) = 0 in Y ( t ) * . If dJ (y, )θ exists for all θ ∈ and that dJ (y, ) defines a bounded linear functional on then we say that J is shape differentiable at .
(A3) Let y ∈ Y ( ) be the unique solution of (19). Then for any f ∈ Y ( ) * the solution of the following linearized equation exists: (A4) Let y t , y ∈ Y ( ) be the solutions of (20) and (19), respectively. ThenẼ and E satisfy Let y ∈ Y ( ) be the solution of (19), and suppose that the adjoint equation, for all z ∈ Y ( ) has a unique solution p ∈ Y ( ). Then, the Eulerian derivative of J at in the direction θ ∈ exists and is given by where κ is the mean curvature of the surface f .
Let us take note that one should be cautious when dealing with the variables y t , y t in (A1), in particular, one should remember the membership of these variables, i.e., y t ∈ Y ( t ) and y t ∈ Y ( ). Furthermore, equation (20), usually, is derived by pulling the equation back into the domain by the change of variables induced by T −1 t . We also mention that originally, assumption (A3) is written as the Hölder continuity of solutions y t ∈ Y ( ) to (20) with respect to the time parameter t ∈ [0, τ ] (Ito et al. 2008). Fortunately, from the same paper, Ito et al. (2008) have shown that assumption (A3) implies the aforementioned continuity. We cite the said result in the following lemma.
Our goal is to characterize, and of course show the existence of the Eulerian derivative of the objective functional Now, from the deformation field T t , we let (ũ t , p t ) ∈ X t be the solution of the equation We perturb equation (25) back to the reference domain which gives us the operatorẼ : where-by denoting (M t ) k the k th row of M t , and v k the k th component of a vector v-the components are defined as follows and the element t ∈ H −1 ( ) is defined by By construction, this operator satisfies (A1), in particular if (ũ t , p t ) ∈ X t solves (25), then the translated element (ũ t , Indeed, for sufficiently small values of t > 0 the coercivity of A t (·) − 1 2 C(·) can be easily verified, as well as the inf-sup condition for the bilinear form d t .
For property (A2), we introduce the linearization of the Stokes operator A : V ∩ H 2 ( ) 2 → V * , and of the bilinear operators B : V × V → V * and C : V × V → H − 1 2 ( out ) 2 which were briefly discussed at the end of Sect. 2.

Proposition 13
The Fréchet derivative of the operators A , B and C at the point u ∈ V in the direction δu ∈ V are given as follows: where we used the notations Cu = C(u, u), Proof We only expose the parts where the nonlinearity occurs, i.e., we show that and Indeed, we have the following computations On the last inequality, we used Rellich-Kondrachov embedding H 1 ( ) 2 → L q (∂ ) 2 for q ≥ 2.
To verify (A4), we note that Thus, by dividing the previous computation by t and by utilizing Lemmata 9 and 12, we infer that E(·) andẼ satisfy (A4).
Having been able to show that the operator E(·) andẼ satisfy (A1)-(A4), and from the fact that the objective functional mostly consists of squared-norms, i.e., P, J ∈ C ∞ (X , R), the last remaining task to assure the existence of the Eulerian derivative of G is the unique existence of the solution (v, π) ∈ X to the adjoint problem where G corresponds to the derivative of the integrand of G with respect to u whose action is given as here the difference between the two curls are as follows: given a vector valued function ϕ = (ϕ 1 , ϕ 2 ) we have ∇ × ϕ = ∂ϕ 2 ∂ x 1 − ∂ϕ 1 ∂ x 2 , while if ψ is a scalar valued function then ∇ × ψ = ∂ψ ∂ x 2 , − ∂ψ ∂ x 1 . The unique existence of the adjoint variables (v, π) ∈ X can quite easily be established following the arguments of Theorem 2. Furthermore, from the regularity of the domain, the adjoint solution v satisfies v ∈ V ∩ H 2 ( ) 2 , and that the solution can be looked at as the solution to the system To finally characterize the shape derivative of G, we start by evaluating where (ũ, p), (v, π) ∈ X solve (23) and (28), respectively. To do this, we write the operator E with the pushed-forward form on t and utilize Lemma 10, that is, we determine the derivative of We note that in the expression ofẼ above, we added the term d(g Sinceũ = v = 0 on f , and divũ = divv = 0, the integrals on the boundary f vanish except for (∇ũ : ∇v, θ · n) f and (∇ g : ∇v, θ · n) f . Furthermore, since (u, p) = (ũ + g, p), we get To further simplify the expression above, we take into account the facts that (u, p) satisfies (10), and (v, π) on the other hand solves (29). We begin the simplification on the terms with ψ v and ψ π , which we shall denote by I 1 .
The computation above utilized the divergence-free property ofũ and the identity (∇v) θ = ∂ n v(θ · n) on f . Furthermore, the boundary integral on ∂ is simplified into just the integral on f since θ = 0 on ∂ \ f . Similarly, if we denote by I 2 the terms containing ψ u and ψ p , we have the following simplification, We note that the last equality is achieved using Green's curl identity (Monk 2003, Theorem 3.29). Combining the expressions obtained from simplifying I 1 and I 2 , and since u = v = 0 on f , then we can further simplify D tẼ as follows: Since divu = divv = 0 in , and u = v = 0 on f , by utilizing the definition of tangential divergence (Sokolowski and Zolesio 1992, p. 82) the integrals (π n · ∂ n u, θ · n) f and ( pn · ∂ n v, θ · n) f both equate to zero. Hence, we get the following derivative, With the computation of the derivative above, we finally characterize the Eulerian derivative of G which is shown in the following theorem.
From the two previous identities, we rewrite the second line in (31) as follows: Therefore, from the assumed regularity of the domain, and by employing divergence theorem, we obtain dG ((u, p), Remark 4 (i) The shape derivative of G as we have formulated it agrees with the Zolesio-Hadamard Structure Theorem (Delfour and Zolesio 2011, Corollary 9.3.1), in it we were able to write the derivative in the form dG ((u, p), )θ = f ∇G(θ · n) ds.
In this case, we shall call ∇G the shape gradient of the objective functional. Furthermore, this form gives us an intuitive gradient descent direction given by θ = −∇G. This fact-together with the challenges with this chosen direction-will further be explored in the subsequent parts of the paper.
(ii) If one observes the adjoint equation (29), we can easily see why (BC2) will not be easily handled. In fact, if such condition is imposed instead of (BC1), it would be impossible to write the adjoint equation in its strong form. Furthermore, the weak form would include the term (ϕ · n) − , where ϕ is a test function. This expression is quite hard to treat numerically due to its discontinuous nature.

Numerical realization
We shall discuss the numerical implementation of the shape optimization problem in this section. We start by discussing the resolution on solving the nonlinearity on the state equations, then we proceed by introducing gradient descent methods based on the Eulerian derivative of the objective functional. The said gradient descent methods will include the rectification of the volume preservation issue, as required in the formulation of the problem. Lastly, we shall show the convergence-with respect to domain discretization-of the final shapes to a manufactured solution in terms of the final deformation fields and in terms of the Hausdorff distance.

Newton implementation of the Navier-Stokes equations
The nonlinearity is one of the challenges not only in the analysis but also in the numerical implementation of the Navier-Stokes equations. In this subsection, we shall discuss how this is resolved by means of a Newton's method.
We begin by reiterating the fact that if (ũ, p) ∈ X is the unique solution of (10), then for any F ∈ X * , there exists a unique solution (δu, δ p) ∈ X to This implies that E (ũ, p) ∈ L(X , X * ) is an isomorphism. Furthermore, by the inverse function theorem there exists a closed ball X ((ũ, p); ε) centered at (ũ, p) ∈ X with radius ε > 0 such that (ũ, p) is an isolated nonsingular solution.
With all the facts presented above, we propose the following Newton's algorithm: we start with an initial element (ũ 0 , p 0 ), and we generate the following sequence {(ũ k , p k )} k using the difference equation for all (ϕ, ψ) ∈ X . In strong form, by letting u k =ũ k + g, we can write (32) as For a given f ∈ L 2 ( ) 2 , and g ∈ H 2 ( ) 2 , one can show-by following the arguments of Girault and Raviart (1986)-that if u ∈ V ∩ H 2 ( ) 2 is the solution to (10) there exists c < 1 such that the following convergence estimate holds Solving numerically, we shall approximate the solution to (10) using (33), with the stopping criterion u k+1 − u k H 1 ( ) 2 < ε for sufficiently small ε > 0. It is noteworthy to mention that Newton's method for the Navier-Stokes equations with (BC1) is naturally constructed thanks to the absence of (·) − , on the other hand the linearization will not be easily handled if (BC2) is used.

Gradient descent methods for deformation fields
As mentioned before, in this part we discuss gradient methods based on the Eulerian derivative of the objective functional. Furthermore, due to the volume constraint we shall utilize two methods, namely the use of an augmented Lagrangian method based on Nocedal and Wright (2006), and of divergence-free deformation fields.

Augmented Lagrangian method
Note that the optimization problem can be written as the equality constrained optimization by With this reason, we formulate the augmented Lagrangian given as where > 0 is a Lagrangian multiplier, and b > 0 is a regularizing parameter. We note that the quadratic term-aside from its regularizing effect-acts as a more strict penalizing term as compared to the usual Lagrangian methods. This method was formalized in the context of shape optimization by Dapogny et al. (2018). So to minimize the objective functional while not neglecting the constraint F = 0, we instead minimize the augmented Lagrangian L. By following the same arguments as for solving the Eulerian derivative of G, we can formulate the derivative of L and is solved as

Smooth extensions of the deformation fields and the unit normal vector
As mentioned in Remark 5, an intuitive gradient descent direction on the boundary f is either θ = −∇G, or θ = −∇L if one chooses to minimize the augmented Lagrangian. However, this choice of θ may cause irregularities to the deformed boundary f . For this reason, methods of approximating the deformation field θ and extending it to the domain have been proposed, for example Azegami and Takeuchi (2006) developed a seminal approach on such smoothing method. In our current problem, we shall adapt such smoothing extension for minimizing the augmented Lagrangian. In particular, we shall be tasked to solve for θ ∈ H 1 ( ) 2 := {ϕ ∈ H 1 ( ) 2 ; ϕ = 0 on ∂ \ f } that solves the Robin equation Meanwhile, we shall also utilize what was proposed by Simon and Notsu (2022b) by minimizing the objective functional itself but uses divergence-free deformation fields, i.e., we shall solve for (θ , ϑ) ∈ H 1 ( ) 2 × L 2 ( ) 2 that solves the following Stokes equation for all (ϕ, ψ) ∈ H 1 ( ) 2 × L 2 ( ) 2 . We note that on both methods ε, ε 1 , ε 2 > 0 are chosen sufficiently small, so that −ε∂ n θ + θ = −∇L and −ε 1 ∂ n θ + θ + ε 2 ϑ n = −∇G on f , and hence are respective approximations for θ = −∇L and θ = −∇G on f . Using the same narrative as above, we shall determine a smooth extension of the unit normal vector n on f , by finding the vector N that solves The impetus for solving for such extension is the appearance of the mean curvature κ in the shape gradients. According to (Henrot and Pierre 2014, Proposition 5.4.8), the mean curvature of the surface f may be determined using the identity κ = div f n = divN for any extension N, given that the surface is sufficiently smooth. Since our domain is C 1,1 -regular, we have the freedom to apply such identity.

Gradient descent iterations
Having been able to lay out the necessary ingredients, we finally write here the iterative scheme of approximating the shape solution. Recall that we introduced the domain variations that utilizes the identity perturbation operator. One of the reason why we take advantage of such perturbation is its usability and simplicity for iterative numerical methods. In particular, for a given initial domain 0 we shall generate a sequence { k } k by means of the difference equation k+1 = k + t k θ k . In this iteration, we note that the only changing boundary is the free-boundary f , so we denote the k th iteration of the mentioned boundary by k f and the k th iteration of ∇G (or ∇ L) by ∇G k (resp. ∇ L k ), i.e., the states and the mean curvature are all evaluated in k instead of , and the deformation fields θ k are determined by solving either (35) or (36), but with k , k f and ∇G k (resp. ∇ L k ) instead of , f , and ∇G (resp. ∇ L), respectively.
The time step on the other hand, just like most gradient descent methods, is determined using a line search method, i.e., for fixed parameters 0 < σ 1 , σ 2 < 1, j ∈ {0, 1, 2, . . .}, and by denotingt (or L( k ) for the case of the augmented Lagrangian ), we take t k =t k j , where j is the smallest integer such that G( k +t k j θ k ) < G( k ) ( resp. L( k +t k j θ k ) < L( k ) ), and such that the perturbed domain k +t k j θ k does not exhibit mesh degeneracy. Lastly, the parameters , b > 0 in the augmented Lagrangian will also be updated iteratively based on (Nocedal and Wright 2006, Framework 17.3). Given initial parameters 0 , b 0 > 0, we generate the iterates k , b k using the following rules where τ > 1 and b > 0 are given parameters. See also the work of Allaire and Pantz (2006) for a detailed discussion of volume preserving Lagrange multiplier schemes. Summarizing these methods, the algorithm 2 below provides the steps for solving an approximate shape solution: Initialization: Choose the parameters ν, α, γ , ε, ε 1 , ε 2 , ε 3 , σ 1 , σ 2 , τ and b; initialize the domain 0 , the solution u 0 (using Newton's method (33)), and the parameters 0 and b 0 (only if the augmented Lagrangian is being implemented).
Step 1: Evaluate G( k ) (or L( k )); solve for the adjoint variable v k from (29), the mean curvature κ k from (37), and the deformation field θ k using (36), or (35) in the case of the augmented Lagrangian method.
Step 2: Set k+1 = k +t k θ k , where t k is chosen appropriately as discussed above.
Step 3: Solve the new solution u k+1 and evaluate G( k+1 is accepted as the approximate solution. Step 4: (Only for augmented Lagrangian method) Update the parameters k+1 , b k+1 using (39).

Numerical examples
This part shows examples of implementations of the algorithm previously discussed. We start with implementing the augmented Lagrangian method, and show its volume preserving limitations by simulating different values of the Lagrangian multiplier 0 . We then simulate examples using the divergence-free deformation fields, and compare the solutions with the augmented Lagrangian method based on the convergence rate (number of iterations) and the volume preserving property. Lastly, we end by showing convergence of shape solutions to a manufactured solution based on the Hausdorff measure on the boundaries and on the H 1 -convergence of the deformation fields. The simulations were all ran using FreeFem++ (Hecht 2012) on an Intel Core i7 CPU @ 3.80 GHz with 64GB RAM. The state, adjoint, and the resolution of the deformation fields (for the divergence-free deformation fields) are solved using triangular Taylor-Hood (P2-P1) finite elements, while the resolution of the mean curvature and of the deformation fields using the augmented Lagrangian method are respectively solved using P1 and P2 finite elements, all of which are solved with the UMFPACK solver. The input function g is a Poiseuille-like function given by g = (1.2(0.25 − x 2 2 ), 0) on in , for simplicity the external force is taken as f = 0, the viscosity parameter is chosen as ν = 1/100, and the Tikhonov parameter is chosen as α = 7 to ensure that G( * ) ≈ 0 at the solution domain * . As for the domain, we consider a rectangular outer boundary and a circular initial free(inner)-boundary, as shown in Fig. 2.
Meanwhile, the domain variation based on the deformation fields will be dealt with by using the movemesh command in FreeFem++, and the possible mesh degeneracy will be circumvented (aside from the choice of the step size) by utilizing the combination of the checkmovemesh and adaptmeshmesh commands.
Lastly, we mention that for the first two subparts of this subsection the mesh size will be taken uniformly and has size h = 1/60, i.e., the diameter of all the triangles in the domain triangulation will be taken as 1/60. Furthermore, the stopping criterion of shape approximation is decided to be tol = 10 −4 .

Augmented Lagrangian method
We start the implementation of the augmented Lagrangian with 0 = 54 since this value is relatively effective in our goal of preserving the volume of the domain compared to other values of 0 that we simulated. We note also that in some figures, we shall utilize scaling of the computed values, such scaling is done in a way that the maximum value is at y = 1, and the minimum value is at y = 0, i.e.,for an iteration k the value in the trend is computed as where F is either the objective function or the volume.
As shown in Fig. 3A, the circular initial shape evolves into a bean-shaped surface but exhibits more convexity than the usual bean shape. The method converges after thirty iterations, where-as seen in Fig. 3B-the value of the objective functional starts with the value G( 0 ) = 0.949 and converges to the approximate solution with value G( 30 ) = 3.49×10 −4 . With regards to the volume preservation, we observe a sudden increase on the first iteration. Nevertheless, the sudden increase on the first iteration is corrected starting from the second iteration, this is caused by the manner by which we update the Lagrange multiplier and the regularizing parameter b. However, the decrease continues until the volume, which started with the value | 0 | = 1.94699, is decreased up to the volume of | 30 | = 1.94649. This decrease on the volume would still be considered a fair volume preservation, with the relative percentage difference 3 equal to 2.568 × 10 −2 %.
For the next illustration, we show the effects of varying the initial value of the Lagrange multiplier 0 . From Fig. 4A, the domain bounded by the surface f gets larger as the value of 0 decreases. This observation implies that the stringency of the equality constraint F ( ) = 0 3 The relative percentage difference is computed with respect to the initial volume, i.e., is stronger as the value of the Lagrange multiplier 0 increases, but the satisfaction of the said constraint is only up to a maximum value. This can be observed in Fig. 4C, where, starting from 0 = 16, the variation 4 increases from around -0.02 to the value zero. However, at 0 = 56 the variation becomes positive and thus the constraint becomes too constricting which-when higher values of 0 are simulated-causes the region inside the surface f to be smaller, and fails the volume constraint. We also plotted with a differently colored hollow point the variation from the initial volume to the volume of the final shape for 0 = 54 which was the parameter value chosen in our first illustration. From the trend of the variational line we can make a conjecture that, given all other parameters are taken constant, the best value of the Lagrange multiplier that will satisfy F = 0 is around 0 = 55. We also observe in Fig. 4B that even though all objective function trends converge to zero, the decrease for 0 = 56 is the least steep. This is because when the value of 0 is smaller, the shape is more allowed to change and thus the algorithm is more relaxed to immediately decrease the value of the objective functional.

Divergence-free method
For the divergence-free method, we get the freedom from tediously choosing an appropriate parameter 0 that will give a volume preserving deformation fields.
We can also observe a bean-shaped surface develop as we go further in each iteration (see Fig. 5A). Contrary to the shape generated by the augmented Lagrangian method, the surface obtained in this current method bounds a domain whose convexity on the left side is lost. Nevertheless, as we can see from Fig. 5B, the objective functional tends to zero as the method reaches its tolerance. Although, a drawback in this method is the apparent slow convergence, as the method converges only after 94 iterations. Fortunately, the volume preservation in this To see even further the difference between the augmented Lagrangian and the divergencefree method with regards to the final shapes, objective value and volume trends, we refer to Fig. 7.
Before we move on, let us first verify if we are fulfilling one of the main goals of the problem which is maximizing the vorticity of the fluid. Even though we are observing a decrease in the objective functional in Figs. 3B and 5B, it is important that we confirm that the descent does not only affect the perimeter of the free-boundary but also minimizes the negative of the L 2 -norm of the curl of the velocity field.
As can be observed in Fig. 8, the perimeter functional is increasing in each iteration, while the quantity ∇ × u 2 L 2 ( ) is also increasing. This verifies and satisfies our goal of maximizing the vorticity of the fluid. Nevertheless, it would be foolish to assume that it is okay to neglect the perimeter part of the objective functional, or even remove it in the definition of the objective functional. As a matter of fact, we recall that this portion of the functional helped us in regularizing the domain variation, and the negligence on this part of the objective functional may cause topological changes in the domain, i.e., additional holes may emerge due to lack of the regularization in the perimeter. Furthermore, the choice of α > 0 in the regularization helps to make sure that G( k ) ≥ 0 which implies that our choice of step size given in (38) is always nonnegative.
To visually observe the effect of the shape solutions on the fluid vortex, we simulate a dynamic version 5,6 of the Navier-Stokes equations and observe the length and width of the twin-vortex right before Karman vortex shedding. Figure 9 shows the comparison of the twin-vortex induced by the initial geometry and the final shape using the divergence-free method. Note that the simulations were done using k (number of iteration)  . 9 The figure shows the generation of twin-vortex before the shedding of Karman vortex. The upper part of the figure corresponds to the upper half of the twin-vortex using the initial shape, while the lower part is done using the final shape from the divergence-free method similar contour levels for the comparison to be reliable. We observe that the length and the width of the vortex is improved when the final shape from the optimization problem is used. The improved length and width of the vortex using the final shape from the augmented Lagrangian method behaves similarly as when the divergence-free method is used, hence we skip such illustration.

Convergence with respect to domain triangulations
For our last set of illustrations, we shall show the convergence of the approximate shape solutions to a reference solution according to the domain triangulation. The reference solution is obtained by using higher order polynomial basis functions and finer domain triangulation. Here, we chose (P4-P3) Taylor-Hood finite elements for the state, adjoint, and the divergencefree deformation fields, while we utilized P3 finite elements for the mean curvature and the augmented Lagrangian deformation fields, we also considered the mesh size h = 1/160. The impetus for generating such solution is due to the notion that higher order polynomials and finer mesh generates a solution that is very near with the exact solution. Hence, we can investigate how approximations due to coarser meshes and lower order polynomials converge to an exact solution by looking at how the approximations converge to the reference solution.
By denoting f,e and f,h the free-boundaries of the reference and the approximate domain solutions, respectively, we shall first show the convergence in terms of the Hausdorff measure 7 between f,e and f,h . We consider the mesh sizes h = 1/10, 1/20, . . . , 1/100 and utilize the same finite elements considered in the previous illustrations. Figure 10 shows the results using the divergence-free method, which we implemented with the different values of mesh size as discussed above. In Fig. 10 A we see the comparison of the final shapes f,h , and the manufactured final shape of the boundary f,e . It can be observed that indeed the shape generated with mesh size h = 1/10 is much coarse as compared to the other cases, and becomes smoother as h becomes smaller. Furthermore, the approximated boundaries get closer to the reference boundary as h → 0.
Aside from the qualitative observation of the domain convergence, Fig. 10B shows how the Hausdorff measure between the approximate and exact boundaries becomes smaller as we decrease the mesh size. As for the preservation of the volume, we can see that the divergencefree method preserves the domain volume efficiently. In fact, we can see from Fig. 10C that the domain variations do not exceed 9 × 10 −4 .
As for the augmented Lagrangian method, Fig. 11 shows the implementation with different values of h as discussed before, but with a fixed Lagrange multiplier 0 = 54.
As we can observe, the approximate solutions converge to the reference solution. Starting with a coarse boundary, the solutions smoothen while converging to the said reference shape (see Fig. 11A).
If we compare it with the approximate solutions generated by the divergence-free method, we can see that the solutions to the divergence-free method can be imagined as results of applying pressure to an inflated balloon. Meanwhile, the approximations using the augmented Lagrangian method are obtained by slicing a bread into portions. This situation is also reflected on the volume preserving capacity of the two methods, as slicing does not conserve the volume of the bread while squeezing/applying pressure retains the volume of the fluid/air in a balloon. This observation can also be witnessed in Fig. 11C, where the volume preservation for the augmented Lagrangian method is much less efficient for higher values of h, while the volume for the shapes generated by the divergence-free method seems to be consistently close to the initial volume regardless of the mesh size.  As for the quantitative convergence of the shapes-by virtue of Hausdorff measure-we can see in Fig. 11B that the distance between the reference and approximate solutions tends to go smaller as the mesh size approaches zero, hence the convergence witnessed in Fig. 11A is quantitatively verified.    For our last set of illustrations, we shall show the experimental order of convergence of the solutions of the state and adjoint equations, and of the deformation fields, all of which are evaluated at the final shapes.
For k = 0, 1, 2, 3 we shall consider the following mesh sizes h k = 10 −1 × 2 −k . Furthermore, we shall respectively denote by u k , v k , and θ k the solutions to the state and adjoint equations, and the deformation fields at the approximate shape solution k when the triangulation has a mesh size h k . Meanwhile, we shall denote by u * , v * , and θ * the counterparts of the quantities previously mentioned but solved on the reference solution * .
For an approximate state ϕ k and an exact state ϕ * , the experimental order of convergence is solved by computing the following quantity where D k is the discretization of the hold-all domain D with mesh size h k and contains the nodes of the approximate solution k . We note that it is imperative to measure the difference ϕ k − ϕ * in D k to make sure that the approximation ϕ k is well defined, while the reference solution is interpolated according to the nodes of the discretized hold-all domain. We also mention that the impetus for solving the eoc for the deformation fields θ k , is that these fields can be looked at as the controls for an optimal control problem. In this way, we verify that the error estimates for our type of controls does not agree with the usual error estimates for optimal control problems 8 , which is due to the evolving domain. From Table 1, we see that using the divergence free method, the error reduction for the state, adjoint and the deformation fields using the L 2 -norm is at most of order 3/2, while for the H 1 -norm is at most linear. Unlike the error estimates for elliptic problems (or even the Stokes/Navier-Stokes equations) with fixed domain where the convergence rate can be easily obtained as quadratic and linear for the L 2 and H 1 norms, respectively.
Using the augmented Lagrangian method on the other hand yields an almost quadratic and linear error reduction for the deformation fields on the L 2 and H 1 norms, respectively.    This convergence may be attributed to the slicing phenomena for approximating the freeboundary that we talked about earlier. However, the irregularity and slow convergence that was apparent in the divergence-free method may still be witnessed in the convergence rates for the state and adjoint solutions. With these observations, we conclude this section by giving some possible reason for the anomaly with the above convergence rates. We note that the approximate shape solution k is not necessarily-or even close to-the discretization of the exact shape * 9 . If we look at Fig. 10A for instance, we can see that the approximate boundary f,h when h = 1/10 is very much distinct from the free-boundary f,e of the exact shape solution. For this reason, the usual error estimates with respect to domain discretizations may not hold. In fact, k is the solution to the fully discretized problem min k G k ( k ) subject to system (33), where all the integrals including that of the objective functional are evaluated using Gaussian quadrature. Now, if the exact shape * is discretized with same mesh size h = 10 −1 × 2 −k , which we denote by * k , and solve the Navier-Stokes equations numerically and denote its solution as u * k , then we can achieve u * − u k H m ( * ) = O(h 2−m ) as an order of convergence. 10 Furthermore, we infer from using triangle inequality that Since the first term on the right-hand side of (41) can be estimated by O(h 2−m ), we conjecture that the second term is what impedes us to fully realize the rate of convergence that is enjoyed by solutions for usual optimal control problems.  Kasumba and Kunisch (2012) for such formulation. Note also that theoretical existence of solution is not guaranteed. d See Braack and Mucha (2014) and Bruneau and Fabrie (1996) e See Braack and Mucha (2014) and Bruneau and Fabrie (1996) f see Remark 4 g see Remark 4

Conclusion
To summarize, we started the exposition by introducing an artificial boundary condition to make sure of the existence of the solutions to the Navier-Stokes equations. We later on showed the existence of states and shape solutions. The improved regularity of the state solution based on an additional regularity assumption on the domain was also briefly discussed. During the course of the discussions, we also unraveled some advantages of the outflow boundary condition we considered as opposed to the artificial conditions considered in the literature. Such properties are summarized in the Table 3. For more information about the difference between the outflows induced by the artificial boundary conditions, we refer the reader to Simon and Notsu (2022). We proceeded by formulating the Eulerian derivative of the objective functional. In the shape derivative, the volume constraint seems to have been neglected which compelled us to adapt an augmented Lagrangian method and to use a class of deformation fields that satisfy the divergence-free property. The final shapes yielded by the two methods were then compared. We observed that the volume constraint is more satisfied when the latter method is utilized. However, we also saw that this method needed much more iterations as compared to the former.
We ended by showing the convergence of numerical solutions to a manufactured solution based on the domain discretization. The convergence was measured in terms of the Hausdorff distance, and of the H 1 and L 2 norms of the state and adjoint solutions, and of the deformation fields. We found out that the experimental convergence rates do not reflect the same convergence rates of typical finite element error estimates for elliptic problems defined in fixed domains. With this in mind, we close this exposition by recommending to future authors the formulation of the theoretical error estimates. In particular, and if possible, we hope to see the development of estimating the gap between the states u k and u * .
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
To establish (A3), we shall refer to the following Hardy inequality (H. Hardy and Pólya 1934, Theorem 330): Lemma 15 Let p > 1, and δ ∈ (0, ∞]. For any u ∈ W 1, p (0, δ) such that u(0) = 0, the following inequality holds: Indeed, from the properties of θ ε , by denoting ε out := {s ∈ out : d(s, 0 ) ≤ 2e −1/ε }, we get w 0,ε L 2 ( out ) 2 = ε out |w 0,ε | 2 ds Let us next obtain an estimate for the first term on the last line of (A5). First let us note that ε out consists of finite number of separate segments, but such segments would always include an endpoint that is attached to w . In our particular set-up, ε out consists of two segments. Nevertheless, we shall only focus on one segment which we shall denote as ε out,1 . This is because we only need to take advantage of the property that φ = 0 on w to utilize Lemma 15 which is satisfied by φ in all segments. The estimate for the first term on the last line of (A5) will then be generated by adding all the estimates computed from each segment which are solved using the same arguments as for ε out,1 . We start by taking a ball B ε centered at the intersection w ∩ ε out,1 with radius ε. From this, we introduce a coordinate system {(x 1 , x 2 )} where the origin is at the center of the ball B ε , the axis of x 1 is tangential to ε out,1 and the axis of x 2 has the direction inward normal to ε out,1 . Since out is at least of class C 1 , there exists h ∈ C 1 (R; R), such that ε out,1 = {(x 1 , x 2 ) : x 2 = h(x 1 ); x 1 ∈ [0, ε]}. From this, we write the integral as follows: ε out |φ| 2 d(s, 0 ) 2 ds = ε 0 |φ(x 1 , h(x 1 ))| 2 d((x 1 , h(x 1 )), 0 ) 2 d(x 1 , h(x 1 )).

Proposition 17
Suppose that the assumptions in Proposition 16 hold. Then, the unique solution δu ∈ V to (A11) exists.
From Lemma 4, by using the quantity N , and since ũ V ≤ 2 ν V we get the following absolute estimate Using the above estimate to (A12) yields Lastly, the uniqueness assumption (9) implies that ν 2 − 4N V > 0 and therefore yields the coercivity.