{'fraction_non_alphanumeric': 0.05410830201571575, 'fraction_numerical': 0.021502391527160915, 'mean_word_length': 3.868295218295218, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A phase field approach to shape optimization in Navier--Stokes flow with integral state constraint

We consider the shape optimization of an object in Navier--Stokes flow by employing a combined phase field and porous medium approach, along with additional perimeter regularization. By considering integral control and state constraints, we extend the results of earlier works concerning the existence of optimal shapes and the derivation of first order optimality conditions. The control variable is a phase field function that prescribes the shape and topology of the object, while the state variables are the velocity and the pressure of the fluid. In our analysis, we cover a multitude of constraints which include constraints on the center of mass, the volume of the fluid region, and the drag of the object. Finally, we present numerical results of the optimization problem that is solved using the variable metric projection type (VMPT) method proposed by Blank and Rupprecht, where we consider one example of topology optimization without constraints and one example of maximizing the lift of the object with a state constraint, as well as a comparison with earlier results for the drag minimization.


Introduction
Fundamental to the design of aircraft and cars, as well as any technologies that would involve an object traveling within a fluid, such as wind turbines and drug delivery in biomedical applications, is the consideration of hydrodynamic forces acting on the object, for example the drag and lift forces. The desire to construct an object with minimal drag or with maximal lift-to-drag ratio naturally leads to the notion of shape optimization in fluids, in which the problem can often be formulated in terms of an optimal control problem with PDE constraints.
Let us assume that Ω ⊂ R d , d = 2, 3, is a bounded domain with Lipschitz boundary, and contains a non-permeable object B. We will denote the boundary of B by Γ ∶= ∂B ∩Ω with the outer unit normal ν, and assume that Γ ∩ ∂Ω = ∅, i.e., the object B never touches the external boundary. A fluid is present in the complement region E ∶= Ω ∖ B, and we assume that the velocity u and the pressure p of the fluid in the region E obey the stationary Navier-Stokes equations with no-slip conditions on Γ, namely, −µ∆u + (u ⋅ ∇)u + ∇p = f in E, (1.1a) div u = 0 in E, (1.1b) Here f denotes the external body force, µ denotes the (constant) viscosity, and g models the inflow and outflow on the boundary ∂Ω such that ∫ ∂Ω g ⋅ ν ∂Ω dH d−1 = 0, where ν ∂Ω denotes the outer unit normal on ∂Ω.
Our present contribution is motivated from a previous numerical study [14] for the shape optimization problem of maximizing the lift-to-drag ratio subject to the PDE constraint (1.1). In two spatial dimensions, the classical formulation of the lift-to-drag ratio is defined as where u ∞ is the flow direction, u ⊥ ∞ is the perpendicular vector and H d−1 is the Hausdorff measure on the set Γ. In [14], using a phase field approximation which we will detail below, the authors obtain an optimal shape similar to a non-symmetric airfoil with a small angle of attack. However, a chief obstacle to a rigorous mathematical treatment of the problem is that it is unknown if the lift-to-drag ratio (1.2) is bounded from above (as we want to maximize the ratio). Furthermore, due to the fractional form of the lift-to-drag ratio (1.2), we also observe fractions entering in the associated adjoint system and optimality conditions computed by the formal Lagrangian method, leading to severe complications in the numerical implementation.
One idea is to study a related problem involving maximizing the lift while the drag is constrained to be below a certain threshold, i.e., max Γ u ⊥ ∞ ⋅ µ ∇u + (∇u) ⊺ − p I ν dH d−1 subject to (1.1) and where D > 0 is a threshold for the drag. In this case, the problematic fractional form is replaced and analysis can be performed on the optimization problem. In exchange we now have to deal with (integral) state constraints, and the difficulty lies in establishing the existence of the associated Lagrange multipliers.
To fix the setting for this paper, we now introduce a design function ϕ ∶ Ω → {±1}, where {ϕ = 1} = E describes the fluid region and {ϕ = −1} = B is its complement. The natural function space for the design functions is the space of bounded variations that take values ±1, i.e., ϕ ∈ BV (Ω, {±1}), which implies that the fluid region E has finite perimeter P Ω (E). If ϕ is a function of bounded variation, its distributional derivative Dϕ is a finite Radon measure which can be decomposed into a positive measure Dϕ and a S d−1 -valued function ν ϕ ∈ L 1 (Ω, Dϕ ) d , where S d−1 denotes the (d − 1)-dimensional sphere. The total variation for ϕ ∈ BV (Ω, {±1}), denoted by Dϕ (Ω), satisfies Dϕ (Ω) = 2P Ω ({ϕ = 1}), and thus we can express the Hausdorff measure H d−1 on the set Γ as 1 2 Dϕ (Ω). Furthermore, the S d−1 -valued function ν ϕ can be considered as a generalized normal on the set ∂{ϕ = 1} (see [1,10,16] for a more detailed introduction to the theory of sets of finite perimeter and functions of bounded variation). For we consider the following general shape optimization problem with perimeter regularization: min (ϕ,u,p) J 0 (ϕ, u, p) ∶= Ω b(x, u, ∇u, p, ϕ) dx In addition, for fixed m 1 , m 2 ∈ N ∪ {0}, we impose the m 1 integral equality constraints and m 2 integral inequality constraints: where for each 1 ≤ i ≤ m 1 + m 2 , for functions K i ∶ Ω×R d ×R d×d ×R×R → R and L i ∶ Ω×R d×d ×R → R d . The parameter γ > 0 in (1.3) is the weighting factor for the perimeter regularization, which is given by the term 1 2 Dϕ , and in light of the above discussion regarding the measure 1 2 Dϕ representing the Hausdorff measure on Γ, we see that the functions b and {K i } m 1 +m 2 i=1 model objectives and constraints in the bulk phases E and B, while h and {L i ⋅ ν ϕ } m 1 +m 2 i=1 model constraints on the interface Γ. Examples of b, h, K i and L i are given below, and it is noteworthy to point out that there is no dependence on u in L i as the no-slip condition (1.4d) ensures that u = 0 on Γ. However, the gradient ∇u may not vanish on Γ, which leads to its appearance in the surface constraints.
The appearance of the perimeter regularization γ 2 Dϕ (Ω) in (1.3) is motivated from the well-known difficulties regarding the mathematical treatment of shape optimizationin particular the existence of minimizers/optimal shapes are not guaranteed [23,26,36]. However, if the shape optimization problem is additionally supplemented with a perimeter regularization, then positive results concerning existence of optimal shapes have been obtained (see for instance [34]).
Let us now give some examples of functions b, h, K and L (where we neglect the index i for convenience) that are of relevance. For a subset A ⊂ Ω, we use the notation χ A (x) to denote the characteristic function of A, i.e., χ A (x) = 1 if x ∈ A and χ A (x) = 0 if x ∈ Ω ∖ A. In particular, one can think of the design function ϕ as ϕ(x) = −1 + 2χ E (x) which satisfies ϕ(x) = 1 for x ∈ E and ϕ(x) = −1 for x ∈ Ω ∖ E = B. Hence, in the following examples for the function b, we can use 1 2 (1 + ϕ) as a restriction to the region E and similarly, 1 2 (1 − ϕ) as a restriction to the region B: • the total potential power of the fluid • the construction cost of the object 1−ϕ 2 w(x), where w denotes a cost function per unit volume, • the least square approximation to a target velocity profile u tar and a target pressure profile p tar in an observation region Q ⊂ E. Here δ 1 and δ 2 denote nonnegative constants.
An example for the surface cost h which has practical applications is the hydrodynamic force component in the direction of the unit vector a, which is given as where I denotes the identity tensor. The drag of the object is given when a is parallel to the flow direction u ∞ , while the lift of the object is given when a = u ⊥ ∞ , the unit vector perpendicular to the flow direction.

Examples of integral constraints that are of interests include
• volume constraints on the amount of fluid -setting K 1 = ϕ − β 1 , L 1 = 0, K 2 = β 2 − ϕ and L 2 = 0 for fixed constants −1 < β 1 ≤ β 2 < 1 leads to inequality constraints: or equivalently • the prescribed mass of the object -setting is a mass density and M > 0 is a target/maximal mass leads to the inequality constraint (1.10) • the prescribed center of mass of the object (with uniform mass density) at a point y in the interior of Ω, i.e., y ∉ ∂Ω -setting K i = 1−ϕ 2 (x i − y i ) and L i = 0 for 1 ≤ i ≤ d leads to the equality constraints • the prescribed drag of the object -setting L 1 = −µ(∇u+(∇u) ⊺ )a+pa, where a is the unit vector parallel to the flow direction u ∞ and D > 0 is a maximal drag value leads to the inequality constraint (1.12) In the examples of the cost functional described above, the problem involving minimizing the drag of the object has received much attention and is well-studied in the literature, see [6,3,28,29,32] and the references therein. For the formal derivation of shape derivatives with general volume and boundary objective functionals in Navier-Stokes flow, we refer the reader to [31], but to the authors' best knowledge, the shape optimization problem with integral state constraint has not received much attention.
In this work, we study a phase field approximation of the problem (1.3)-(1.6), which is derived in Sec. 2. Under appropriate assumptions we prove the existence of minimizers to the phase field optimization problem and derive the first order optimality conditions. The main difficulty we encounter is establishing the existence of Lagrange multipliers, which is achieved via constraint qualifications such as the Zowe-Kurcyusz constraint qualification [39], for some of the integral state constraints mentioned above. We give two examples: one involves constraints that depend only on the design function ϕ which are the volume (1.9), the prescribed mass (1.10) and the center of mass (1.11), while the second example involves the total potential power (1.7). We encounter some technical difficulties regarding the drag constraint (1.12), and can only show the existence of Lagrange multipliers if the threshold D is sufficiently large. Via numerical simulations we give a proof of concept showing that with the help of the phase field approach shape and topology optimization for fluid flow taking state constraints can be solved. For large Reynolds number problems more efficient numerical solution methods have to be devised in the future.
The rest of the paper is organized as follows: In Sec. 2, we present the phase field approximation of (1.3)-(1.6) that utilizes the porous-medium approach of Borrvall and Petersson [7], and state several preliminary results on the state equations. Then, in Sec. 3, we state the assumptions on b, h, K i and L i that allows us to establish the existence of minimizers to the phase field shape optimization problem. Sufficient conditions on the differentiability of b, h, K i and L i are outlined in Sec. 4 which lead to the existence of Lagrange multipliers, the solvability of the adjoint system, and the derivation of the necessary optimality conditions. We verify the aforementioned conditions in Sec. 5 for two specific examples of integral constraints; the first example involves constraints on the mass, center of mass and volume, while the second example involves a constraint on the total potential power. Lastly, in Sec. 6 we briefly outline our numerical approach to solving the optimality conditions, and present several numerical simulations.

Phase field formulation
One approach to tackle shape optimization problems that can yield rigorous mathematical results is to employ a phase field approximation, similar in spirit to Bourdin and Chambolle [8] that was applied to topology optimization (see also [4,27,35,38] and the reference cited therein), and has been recently used for drag minimization in stationary Stokes flow [12] and in stationary Navier-Stokes flow [11,13,14,24].
The approach we take in this paper is similar to the previous works [11,12,14], in which we relax the condition that the design function ϕ takes only values in {±1} (i.e., ϕ ∈ BV (Ω, {±1})) and now allow ϕ to be a function with values in R and inherits H 1 (Ω) regularity. In particular, we change the admissible space of design functions from subsets of BV (Ω, {±1}) to subsets of H 1 (Ω). This leads to the development of interfacial layers {−1 < ϕ < 1} in between the fluid region E = {ϕ = 1} and the object region B = {ϕ = −1}. This interfacial layer replaces the boundary Γ of B and a parameter ε > 0 is associated to the thickness of the interfacial layer. The idea is to reformulate the original shape optimization problem (1.3)-(1.6) by taking into account the above modification of the design functions. For the perimeter regularization, we can use the scaled Ginzburg-Landau energy functional where Ψ is a potential with equal minima at ϕ = ±1, to approximate the perimeter functional P Ω . The positive constant c 0 is dependent only on the potential Ψ via the relation and it is well-known that 1 2c 0 E ε approximates ϕ ↦ 1 2 Dϕ (Ω) = P Ω ({ϕ = 1}) in the sense of Γ-convergence [25].
By introducing an interfacial region between the fluid and the object, we have relaxed the non-permeability assumption of the object in the vicinity of its boundary. Therefore, we use the so-called porous medium approach of Borrvall and Petersson [7] and replace the object B with a porous medium of small permeability (α ε ) −1 ≪ 1. A function α ε (ϕ) is introduced to interpolate between the inverse permeabilities of the fluid region α ε (1) = 0 and the porous medium α ε (−1) = α ε , which satisfies α ε → ∞ as ε → 0.
With this, we extend the state equations from E to the whole domain Ω by the addition of the porous-medium term α ε (ϕ)u: (2.2c) We note that this additional term vanishes in the fluid region, and in the limit ε → 0, one expects the velocity u in the object region to vanish. We point out that in the modified state equations ( h(x, ∇u, p, ∇ϕ) dx to the surface objective involving the function h and from this point onward we will assume that h is one-homogeneous with respect to its last variable. By a similar argument, we see that is a phase field approximation of the surface integral constraint involving L i . Recall that in the modified state equations (2.2) the porous-medium term α ε (ϕ)u serves to enforce the condition that the velocity u in the object region should vanish in the limit ε → 0. In earlier works motivated by the paper of Borrvall and Petersson [7] the authors of [12,18] also added to the phase field objective functional whereα ε is a function with similar properties as α ε , i.e.,α ε (1) = 0 andα ε (−1) → ∞ as ε → 0. In fact, in the rigorous analysis of the phase field approximation in Stokes flow, the addition of penalization term (2.3) to the objective functional does indeed lead to the velocity field vanishing in the object region as ε → 0 (see [12, §3] and [18, §6.3] for more details). In this paper we consider including both elements in the analytical treatment of the optimization problem. It is also possible to considerα ε = 0, however in this case no rigorous results on the sharp interface limit ε → 0 are known. Before we state the phase field optimization problem, let us mention that for the analysis we assume that the function α ε in the porous-medium term satisfies the properties: (A0) α ε ∈ C 1,1 (R) is non-negative, and there exist constants s a , s b ∈ R with s a ≤ −1 and In particular, for arbitrary φ and its truncationφ ∶= max(s a , min(s b , φ)) we see that α ε (φ) = α ε (φ), and so the state equations (2.2) for φ andφ are equivalent. Hence, without loss of generality, we now search for optimal design functions ϕ exhibiting H 1 (Ω)-regularity and satisfies the pointwise bounds s a ≤ ϕ ≤ s b a.e. in Ω.
Taking into account the above discussions, we arrive at the following phase field approximation to the optimal control problem (1.3)-(1.6): satisfying the weak formulation of (2.2): for all v ∈ H 1 0 (Ω) ∶= {h ∈ H 1 (Ω) h = 0 on ∂Ω}, along with m 1 equality and m 2 inequality integral constraints

Preliminaries on the state equations
Since the porous medium Navier-Stokes equation (2.2) have been analyzed in detail in previous works [11,18], we recall some useful results in this section.
is non-negative and satisfies (2.4), for every ϕ ∈ L 1 (Ω) there exists at least one pair (u, p) ∈ H 1 g,σ (Ω) × L 2 0 (Ω) such that (2.6) is satisfied. Furthermore, there exists a positive constant C = C(µ, α ε , f , g, Ω) independent of ϕ such that u H 1 (Ω) + p L 2 (Ω) ≤ C. (2.8) The estimate (2.8) can be established by testing the weak form (2.6) with u−G, where G ∈ H 1 g,σ (Ω) is a vector field depending on the inflow/outflow g and the domain Ω, and satisfying certain properties. Furthermore, this computation also shows that the constant C depends on ε only through the function α ε .
By the above existence result, we can define a set-valued solution operator for any ϕ ∈ L 1 (Ω). Next, we state a continuity property of the solution operator. .
Then, there exists a subsequence, denoted by the same index, and functions u ∈ H 1 (Ω), p ∈ L 2 (Ω) such that with the property that (u, p) ∈ S ε (ϕ). Furthermore, it holds that In general, we do not have uniqueness of solutions to (2.6), but there is a conditional uniqueness result.
That is, there is exactly one solution of (2.6) corresponding to ϕ ∈ L 1 (Ω).
The additional assumption (2.10) on the solution u ∈ S ε (ϕ) to ensure uniqueness of the state equations can be achieved for small data f and g or with high viscosity µ. However, there are also settings in which (2.10) can be justified a posteriorly [11]. For the subsequent analysis, more precisely in showing the differentiability of the solution operator S ε and the derivation of the optimality conditions, we require that S ε is a oneto-one mapping. Hence, throughout the rest of the paper we assume that (2.10) holds. Alternatively, instead of assuming (2.10), we can work with an isolated local solution to (2.6), for which the subsequent analysis is valid in a neighbourhood of this isolated local solution.
We now state the Fréchet differentiability of the solution operator S ε .
Then, there exists a neighborhood N of ϕ ε in H 1 (Ω) ∩ L ∞ (Ω) such that for every δ ∈ N , the solution operator consists of exactly one pair, and we may write which is the weak formulation of the following the linearized state system

Existence of a minimizer
We make the following assumptions for the potential Ψ and the functionsα ε , b, h, K i , and L i .
such that for a.e. x ∈ Ω it holds for any r ≥ 0, p ≥ 2 in two-dimensions and 2 ≤ p ≤ 6 in three-dimensions, Carathéodory function that is onehomogeneous with respect to its last variable and there exist non-negative functions for functions k i ∈ L 1 (Ω) and Carathéodory functions such that for a.e. x ∈ Ω it holds for any r ≥ 0, p ∈ [2, ∞) in two-dimensions and p ∈ [2, 6) in three-dimensions, and there exist non-negative functions z 6 ∈ L 1 (Ω), z 7 , z 8 ∈ L ∞ (Ω) such that for a.e.
x ∈ Ω it holds, for all s ∈ R and A ∈ R d×d .
(A7) We assume that the set is non-empty, where the phase field integral state constraints G k , for 1 ≤ k ≤ m 1 +m 2 , are of the form (2.7) and we recall the set Φ is defined as (Ω) are bounded from below, B is weakly lower semicontinuous, and for all ϕ n ⇀ ϕ in The particular forms of b and K i are motivated from the discussions in Sec. 1, where z and y would typically be functions of the form 1+ϕ 2 , and the function k i would be of the form D Ω −1 . Furthermore, the set K ad is the set of admissible design functions whose elements satisfy the m 1 equality integral constraints and m 2 inequality integral constraints. While we assume the non-emptiness of K ad for the general setting here, later in Sec. 5 we show for two examples that the corresponding set K ad is indeed non-empty.
The following weakly closed property is useful for showing the existence of minimizers to the optimal control problem (2.5)-(2.7).
For each 1 ≤ i ≤ m 1 + m 2 , by the continuity of L i with respect to its second and third variables, it holds that L i (x, ∇u n j , p n j ) → L i (x, ∇u, p) a.e. in Ω. Using the growth conditions in (A6), the strong convergences for {u n j , p n j } j∈N and the generalized Lebesgue dominated convergence theorem leads to Together with the weak convergence ∇ϕ n j to ∇ϕ in L 2 (Ω), we have Note that s a ≤ ϕ n j , ϕ ≤ s b a.e. in Ω for all j ∈ N, and thus there exists a constant M > 0 such that sup x∈Ω y i (x, ϕ n j ) , y i (x, ϕ) ≤ M for all n ∈ N. Using the splitting we can show that lim n→∞ G i (ϕ n j , u n j , p n j ) = G i (ϕ, u, p) once we demonstrate that I 1 , I 2 → 0 as n → ∞. This would then imply that ϕ ∈ K ad . Using the growth conditions in (A5) for K i , the strong convergences for {(u n j , p n j )} j∈N and the generalized Lebesgue dominated convergence theorem yields that Then, the assertion that I 1 → 0 as j → ∞ follows from the above strong convergence in L 1 (Ω) and the boundedness of , and the application of the usual Lebesgue dominating convergence theorem yields and hence I 2 → 0 as n → ∞.
We state the existence result for a minimizer of the problem (2.5)-(2.7). Proof. By (A8), (B + H) K ad ×H 1 g,σ (Ω)×L 2 0 (Ω) is bounded from below by a constant C 0 ∈ R. Then, by the non-negativity ofα ε and Ψ, we find that there exists a constant C 1 ∈ R such that J ε ∶ K ad × H 1 g,σ (Ω) × L 2 0 (Ω) → R is bounded from below by C 1 . Thus, we can choose a minimizing sequence (ϕ n , u n , p n ) n∈N ⊂ K ad × H 1 g,σ (Ω) × L 2 0 (Ω) such that (u n , p n ) ∈ S ε (ϕ n ) for all n ∈ N and Then, for arbitrary η > 0, there exists N ∈ N such that for n > N , The above estimate implies that Thus, we may choose a subsequence (ϕ n k ) k∈N such that ϕ n k → ϕ strongly in L p (Ω) and almost everywhere in Ω for 2 ≤ p < ∞ in two-dimensions and 2 ≤ p < 6 in three-dimensions. Furthermore, by Lem. 3.1 we also have that ϕ ∈ K ad , and by Lem.
for some (u, p) ∈ S ε (ϕ), and The continuity of Ψ together with the fact that (ϕ n k ) k∈N ⊂ L ∞ (Ω) implies (Ψ(ϕ n k )) k∈N is a bounded sequence in L ∞ (Ω). The application of the dominated convergence theorem yields that Ψ(ϕ n k ) converges strongly to Ψ(ϕ) in L 1 (Ω) as k → ∞. Furthermore, by the weak lower semicontinuity assumptions of B and H, and the weak lower semicontinuity of the mapping ϕ ↦ ∇ϕ 2 L 2 (Ω) , we find that From this point onwards, for fixed ε > 0, we denote a minimizer to the optimal control problem (2.5)-(2.7) as ϕ ε with corresponding unique solution (u ε , p ε ) to the state equation (2.6).

Optimality conditions
We use the notation D j f to denote the partial derivative of f with respect to its jth variable. Furthermore, the notation D (i,j) f ≤ P means that the partial derivatives D i f and D j f satisfy D i f ≤ P and D j f ≤ P . To obtain optimality conditions, we make the following assumptions on the differentiability of B, z, h, K i , y i , and L i .
, and the partial derivatives exist for all w ∈ R d , s ∈ R, A ∈ R d×d , and a.e. x ∈ Ω as Carathéodory functions with where p ≥ 2 in two dimensions and p ∈ [2,6] in three dimensions.

Fréchet differentiability of the objective functional
We now show that, as a mapping from N ⊂ H 1 (Ω)∩L ∞ (Ω) → R, j ε is Fréchet differentiable at ϕ ε . As Lem. 2.4 guarantees the Fréchet differentiability of the solution operator S ε (ϕ ε ) as a mapping from N to H 1 (Ω) × L 2 (Ω), we focus on the dependence of J ε on the first variable.

Necessary optimality conditions
Now we can formulate the first order necessary optimality conditions for our optimal control problem.

Verification of constraint qualification
In this section, we consider a model problem of minimizing the drag subject to constraints on the mass, center of mass and volume of the object. More precisely, in a bounded domain Ω ⊂ R 2 with Lipschitz boundary, we study the following optimal control problem min (ϕ,u,p) Ω subject to (ϕ, u, p) solving the porous-medium Navier-Stokes equations (2.6) and the following integral constraints: where a is a constant unit vector parallel to the flow direction u ∞ , M > 0 is a given positive constant representing an upper bound on the mass of the object, ρ(x) ∈ L ∞ (Ω) is a non-negative mass density, and β ∈ (−1, 1) so that the object is constraint to occupy a maximal volume of 1−β 2 Ω . The constraints G 1 (ϕ) = 0 and G 2 (ϕ) = 0 imply that the centre of mass for the object is located at the origin in R 2 (which we can assume to hold without loss of generality by translating the domain Ω). We point out that one can also consider more general surface objective functionals h that are one-homogeneous with respect to the last variable, as well as volume objective functionals b, however we consider this particular example of drag minimization as a practical application of our present approach. Furthermore, in this example we have chosen to neglect the penalization termα ε u 2 in the objective functional.
It is straightforward to check that the function is bounded from below for ϕ ∈ H 1 (Ω) ∩ L ∞ (Ω) with s a ≤ ϕ ≤ s b a.e. in Ω, u ∈ H 1 g,σ (Ω) and p ∈ L 2 0 (Ω), and satisfies due to the product of weak-strong convergence: lim n→∞ H(ϕ n , u n , p n ) = H(ϕ, u, p) for sequences ϕ n ⇀ ϕ in H 1 (Ω), u n → u in H 1 (Ω) and p n → p in L 2 (Ω). Hence, (A8) is also fulfilled. A short computation shows that and as a is a constant vector, one can infer that (B1) (specifically (4.1)) is also fulfilled. Then, it remains to verify (A5), (A6), (B2) and (C1) for the existence of Lagrange multipliers for the integral constraints G 1 , . . . , G 4 , and show that the admissible set K ad is non-empty. For the latter, note that we have the trivial example φ ≡ 1 ∈ K ad which corresponds to the case where there is no object in the domain Ω. In the following we will construct a non-trivial example in order to rule out the possibility where K ad = {1}, which would imply the solution to the shape optimization problem is to have no object at all. We can always choose a function φ ∈ H 1 (Ω), −1 ≤ φ ≤ 1 a.e. in Ω such that which is equivalent to choosing an object {φ = −1} with its centre of mass at the origin with volume bounded above by 1−β 2 Ω . Note that the mapping φ ↦ ∫ Ω 1 2 ρ(x)(1 − φ) dx is continuous, and thus we can always decrease the volume of the object region {φ = −1} to ensure the mass is bounded from above by the constant M . This ensures that φ ∈ K ad and hence (A7) is satisfied.
As Ω is a bounded domain, the functions x 1 , x 2 are bounded. Then, upon setting we observe that (A5), (A6) and (B2) are fulfilled by the above choices. Then, by Thm. 3.2 we are guaranteed the existence of a minimizer ϕ ε to the optimal control problem. To verify the main assumption (C1) and derive the optimality conditions, we have to show that for an arbitrary z = (z 1 , z 2 , z 3 , z 4 ) ⊺ ∈ Y = R 4 , there exists one function ψ * ∈ Φ, along with non-negative constants τ 1 , . . . , τ 4 , ξ 1 , ξ 2 , η 1 , η 2 such that the following four conditions are fulfilled simultaneously: Due to their nature as equality constraints, we can use the fact that G 1 (ϕ ε ) = G 2 (ϕ ε ) = 0 to simplify (5.1a) into We first argue for (5.2). As the origin 0 ∉ ∂Ω, this implies that Ω has non-empty intersections with the four quadrants of R 2 , which we denote by Q 1 = {x 1 , x 2 > 0}, Q 2 = {x 1 < 0, x 2 > 0}, Q 3 = {x 1 , x 2 < 0} and Q 4 = {x 1 > 0, x 2 < 0}. If z 1 (resp. z 2 ) is zero, we choose τ 1 (resp. τ 2 ) to be zero. Thus, it is sufficient to focus on the case where z 1 and z 2 are non-zero, and in this case we consider a function ψ * ∈ Φ not identically equal to 1 with β < ψ * ≤ 1 a.e. in Ω such that the non-empty set A ∶= supp(1 − ψ * ) has Lebesgue measure and satisfies Then, we set and thanks to the fact that ψ * ≤ 1 a.e. in Ω, the function 1 − ψ * is non-negative in Ω and only positive in A. The location of A implies that the integrand (1 − ψ * )x i has the same sign as z i for i = 1, 2, and so τ i is positive for i = 1, 2. The condition on the Lebesgue measure of A is used to satisfy the mass constraint. For the inequality constraint, we have to show that the same function ψ * considered above simultaneously satisfies (5.1b) and (5.1c). We argue for the mass constraint, and the volume constraint follows along a similar argument. There are two cases to consider: suppose the inequality constraint G 3 (ϕ ε ) is not active for the minimizer ϕ ε , i.e., ϕ ε satisfies Then, we can choose τ 3 = 0 and it holds that Hence, we have fulfilled (5.1b) without making use of the function ψ * . On the other hand, if G 3 (ϕ ε ) is active, i.e., ∫ Ω 1 2 ρ(x)(1 − ϕ ε ) dx = M , the condition (5.1b) simplifies to A short calculation using (5.3) shows that the quantity in the bracket is positive, and so For the volume constraint (5.1c) we again divide the argument into two cases: if G 4 (ϕ ε ) is inactive, then (5.1c) holds automatically without the use of the function ψ * , and if G 4 (ϕ ε ) is active, then using ψ * > β yields the desired result.
As a consequence, (C1) is fulfilled and we obtain the existence of Lagrange multipliers λ 1 , λ 2 ∈ R, λ 3 , λ 4 ∈ R ≥0 . By Thm. 4.2 the first order optimality condition is together with the complementary slackness conditions Remark 5.1. We point out that the mass constraint G 3 (ϕ) = M − ∫ Ω 1 2 ρ(x)(1 − ϕ) dx ≥ 0 can also be thought of as a constraint on a construction cost, where the value ρ(x) > 0 represents the cost of building the object at the point x ∈ Ω, and M denotes a maximal cost.
Let us now consider a similar model problem but with the single integral constraint on the total potential power (1.7). More precisely, in a bounded domain Ω ⊂ R 2 with Lipschitz boundary, we study the following optimal control problem min (ϕ,u,p) Ω subject to (ϕ, u, p) solving the porous-medium Navier-Stokes equations (2.6) (with zero body force f = 0) and the following integral constraint: where a = −u ⊥ ∞ is the negative unit vector perpendicular to the flow direction u ∞ and D > 0 is a given positive constant representing an upper bound on the total potential power. Then, upon setting we see that (A5), (A6) and (B2) are fulfilled. In this setting we observe that the trivial example ϕ ≡ −1 belongs to the admissible set of design functions K ad . A non-trivial example can be found if the domain Ω is sufficiently large or the viscosity µ is sufficiently small. Indeed, let ϕ be a function in H 1 (Ω) with −1 ≤ ϕ ≤ 1 a.e. in Ω but not identically equal to 1 or −1. Denote by u the unique velocity field associated to the state equation (2.2), then by (2.10) it holds that where from (2.11) the constant K Ω is K Ω = 1 2 Ω 1 2 in two dimensions. Note that the above upper bound is independent of ϕ. For any fixed positive constant D, we can take a sufficiently large domain Ω or sufficiently small viscosity µ, so that µ 3 2K 2 Ω ≤ D. Then, this implies that K ad is non-empty and thus (A7) is fulfilled. Furthermore, following the arguments above, we deduce by Thm. 3.2 that there exists at least one minimizer ϕ ε to the optimal control problem. Writing G(ϕ ε ) = G(ϕ ε , S ε (ϕ ε )), to verify the assumption (C1), we have to show for an arbitrary z ∈ R, there exists one function ψ * ∈ Φ along with non-negative constants τ, ξ, η such that Observe that for this particular setting (with a large domain Ω or small viscosity µ), the inequality D ≥ µ 3 2K Ω 2 holds independently of the minimizer ϕ ε , and thus by (5.4) G(ϕ ε ) > 0 holds. In particular, the constraint is always inactive, and we do not need to find the function ψ * as Furthermore, as the constraint is always inactive the complementary slackness condition (4.8) implies that the associated Lagrange multiplier λ is zero. Hence, by Thm. 4.2 the first order optimality condition is

Numerical implementation and simulations
Let us now describe how we can use the above results to compute optimal shapes and topologies in given flow settings. Since our optimization variable is a phase field, and thus has the natural regularity ϕ ∈ H 1 (Ω) ∩ L ∞ (Ω), we use the variable metric projection type (VMPT) method proposed in [5] to solve the resulting minimization problems. A standard projected gradient method can not be used for the constraint minimization problem due to the fact that H 1 (Ω)∩L ∞ (Ω) is not a Hilbert space. The VMPT method uses derivative information which can be represented with the help of the adjoint variables as specified in (4.9). For the potential function Ψ we use the double-obstacle free energy, namely From this we obtain the constraint ϕ ε ≤ 1, and c 0 = π 2 , where c 0 is the constant defined in (2.1). Although the double-obstacle potential (6.1) does not satisfy (A1), the analysis is not affected once we choose s a = −1 and s b = 1, so that ϕ ε ≤ 1 and the potential becomes Ψ(ϕ ε ) = 1 2 (1 − ϕ 2 ε ). We refer the reader to [11,13] which also uses the double-obstacle potential (6.1). For the porous-medium term α ε (ϕ ε ) in the state equations (2.2a) we choose

Spatial discretization
We use finite elements for the numerical discretization of the minimization problem. We use piecewise linear and globally continuous finite elements for the representation of ϕ ε , p ε and π ε and piecewise quadratic and globally continuous finite elements for u ε and q ε on a conforming triangulation of the domain Ω. It is well-known that in phase field applications the variable ϕ ε changes rapidly across the interfacial layers, and an adaptive concept for its spatial resolution is indispensable. Hence, for the mesh generation we use the Dual Weighted Residual (DWR) method [2] where our implementation is guided by [20]. This generates adaptive meshes which well resolve the interfacial regions, and also well reflect the underlying flow physics, compare also [19]. The DWR approach is only applicable if for a given triangulation an optimal solution is already found and uses this information to calculate error indicators.
For fast calculations, it is desirable to use coarse meshes. In the core of the VMPT method we solve projection-type problems using a primal-dual-active-set strategy (PDAS).
Here the active set corresponds to degrees of freedom with ϕ ε = 1. Thus in every step of the PDAS we solve the problems on the inactive set ϕ ε < 1 only. Note that the integral constraints have to be fulfilled by changing the phase field on the inactive set only. If this set contains too few degrees of freedom, the PDAS is not successful in solving the projection-type problem and thus the algorithm breaks down.
To overcome this numerical issue on coarse meshes, we additionally require that a given amount, say 2%, of the phase field's degrees of freedom are inactive. If this is not the case, we use mesh adaptation that is based on ϕ ε only, namely we use the jumps of the normal derivatives of ϕ ε across edges as proposed in [13] to generate new degrees of freedom inside the interface to be able to proceed with the PDAS.
We stop the adaptation loop as soon as a given maximum number of degrees of freedom is reached.

Topology optimization -a tube through heavy ground
Although we have mainly focused on shape optimization with the phase field approach in this paper, we point out that using a phase field variable for the representation of the unknown shape also allows us to deal with situations where no a priori toplogical information is available. In particular, the phase field approach is capable of topology optimization, as done in [4,27,35,38]. Here we consider the situation where the domain Ω = (0, 1) 2 contains several impermeable rocks and we would like to search for a tube that connects the inflow at the bottom to the outflow at the top, see Fig. 1.
Constructing a tube through the rocks is expensive and therefore a tube that avoids these regions is desired. So this is a setting where we want to minimize the cost of an object. The inflow and the outflow regions as well as the location of the rocks are a priori known. We define the inflow and outflow conditions as For the objective functional we define a 'rock' centered at m with radius σ and associated cost c as We consider the functions The optimization problem (2.5) then becomes subject to ϕ ∈ Φ, u ∈ H 1 g,σ (Ω), p ∈ L 2 0 (Ω) satisfying (2.6), and α ε was defined earlier in (6.2). For this example we do not apply any integral constraints, as it serves to demonstrate the strength of the phase field approach in being able to deal with situations where no a priori topological information is known. Having the solution to the unconstrained problem at hand, one might reduce for example the size of the tube by imposing additional volume constraints. However, specifying such constraints beforehand might lead to inadmissible situations.
We start the optimization procedure with no prior information, i.e., ϕ 0 ε ≡ 0, on a homogeneous coarse grid with mesh size h = 1 20 yielding 685 degrees of unknowns for ϕ ε . We stop the solution and adaptation procedures as soon as an optimal solution with more than 100000 degrees of freedom is found. The numerical parameters are: σ = 0.15, c = 50, ε = 0.01, α = 5, µ = 0.02 and γ = 0.001. To stress the advantages of our approach in Fig. 2 we show ϕ ε at various stages of the optimization procedure.

Reproduction of results on drag minimization from earlier works
We now reproduce the numerical results for the surface formulation of drag minimization presented by the authors in [14]. The key distinction is that in [14] a gradient flow approach is used to solve the optimality conditions, leading to a non-linear time-dependent equation of Cahn-Hilliard type for ϕ ε . However, here we employ the VMPT method to solve the optimization problem, which reads min (ϕ,u,p) subject to ϕ ∈ Φ, u ∈ H 1 g,σ (Ω), p ∈ L 2 0 (Ω) satisfying (2.6) and the volume constraint (see (1.9)) Ω ϕ dx ≤ β 2 Ω for β 2 ∈ (−1, 1). (6.3) Figure 2: The iterations 20,30,40,80,120,190 of the VMPT to minimize the cost of a tube through heavy ground. We see that after 40 steps already an optimal structure is found and that in subsequent steps mostly the resolution of the structure is improved. Let us also note that at iteration 20 we have only 923 degrees of freedom for ϕ ε and in iteration 40 still only 1398 degrees of freedom. The final iteration has 125069 degrees of freedom.
We use the parameters from [14], namely Ω = (0, 1. To be able to use only a small number of unknown as long as possible, we start the optimization with ε = 0.008 and a maximum number of allowed degrees of freedom of 10000. We halve the value of ε as soon as an optimal solution is found with current maximum allowed number of degrees of freedom and increase this value by 20%, resulting in 45000 unknowns for the final result. In Fig. 3 we show the optimal shape for different values of ε, namely ε ∈ {0.008, 0.004, 0.002, 0.001, 0.0005, 0.00025}. In Table 1 we show the diffuse interface drag and the sharp interface drag by evaluation with a = (1.0, 0.0) ⊺ over Γ = {ϕ ε = 0}. We reproduce the results found in [14] where a gradient flow approach is applied that is based on an artificial time evolution. We stress that, in using a gradient flow approach, the interface has to be resolved in each time step of the temporal evolution, which leads to a large numerical effort. To be precise, while the results shown here are found in a few hours using the VMPT method, the results in [14] required several days of calculation using the gradient flow.   [14] and different values of ε. Note that α ε (−1) → ∞ for ε → 0, i.e., the object becomes less permeable and thus the drag increases with ε → 0. In [14] for ε = 0.00025 we observed F D ε = 3.9117 × 10 −2 and F D = 3.9499 × 10 −2 .

Comparison of volume and surface formulations for drag
Let us point out that the hydrodynamic force component (6.5) in its classical representation as a surface integral over Γ can be expressed in terms of a volume integral over the fluid region E, and this reformulation has been used extensively in numerical simulations, see [9, §5.1], [15, §2.2], and [22, §9]. Given the unit vector a ≠ 0, let η be a smooth vector field such that η = a on Γ and η = 0 on ∂Ω. (6.6) This can be done since it is assumed in the introduction that Γ does not intersect with ∂Ω. Then, by taking the scalar product of (1.4a) with η, we obtain Integrating by parts and noting that the boundary integrals over ∂Ω vanish owning to Here we have also used that u has no tangential component on Γ due to the no-slip condition u = 0 on Γ, and together with the divergence-free condition, we obtain that (∇u) ⊺ ν = 0 on Γ (see [14, §2] for more details). This implies that we can also consider the following function as the volume formulation of the drag (if a is parallel to the flow direction) Alternatively, using integration by parts and the boundary conditions η = 0 on ∂Ω and u = 0 on Γ, we see that and so we may also use the function as a volume formulation for the drag. The corresponding phase field approximations of (6.8) and (6.9) have the exact same form. However, for our numerical investigations, we use the formation (6.9) instead of (6.8).
those obtained with the surface formulation (6.4) for the drag. For the above optimization problem, we consider the same setup as in Sec. 6.3 and set η ≡ a on (0.15, 1.0)×(0.13, 0.27).
Using the surface formulation (6.4) we observe that for larger values of α (the constant in (6.2)) an interfacial region { ϕ ε < 1} that is neither fluid nor object appearing in front of the object. A similar behaviour was observed in the previous work [14] with another minimization algorithm. In any case, a sufficiently impermeable object can be obtained by using smaller values of ε. We stress that, in Sec. 6.3 for ε = 0.00025 the the velocity u ε inside the object is five orders of magnitude smaller than outside the object (see [14, Fig. 1]).
On the other hand, using the volume formulation (6.9) we have to define the extension of the unit vector field a, namely the vector field η which has to vanish at ∂Ω. We define η as the solution of a Poisson problem on Ω with η = a on a square around the object and η = 0 on ∂Ω. That is, let S denote a square such that {ϕ ε = −1} ⊂ S and ∂S ∩ ∂Ω = ∅, then we solve −∆η = 0 in Ω ∖ S, η = 0 on ∂Ω, η = a in S. (6.10) For small values of α, we observe that the object splits and the solid is collected close to the inflow outflow boundaries. We believe this behavior is due to the following: On the one hand, due to the boundary condition η = 0 on ∂Ω, the magnitude η is small close to the inflow and outflow boundaries, which results in small drag forces. On the other hand, for α small, the porous-medium penalization term ∫ Ω α ε (ϕ) u 2 dx is small, and thus the value of the objective functional can be reduced by placing material in regions where η is small. Therefore, in contrast to the surface formulation (6.4), large values of α are needed for the volume formulation to obtain reasonable optimal shapes, which additionally allows us to construct sufficiently impermeable objects when we use larger values of ε. We use the setup from Sec. 6.3 with only one modification, that we set µ = 0.01. In Fig. 4 the optimal shapes of the objects using the surface and the volume formulations of the drag are shown. We observe that the front of the object with both formulations is rather similar, while the surface formulation leads to a less pronounced rear. The corresponding drag values in sharp interface evaluation (6.5) as defined in Sec. 6.3 are F D = 0.106467052 (volume formulation) and F D = 0.106470276 (surface formulation).
As described above, using the volume formulation we can use larger values for α to model objects with smaller permeability. To show the influence of α in Fig. 5 we show the optimal shape for the above parameters, but using a larger value α = 1 and µ = 0.01 (left) and µ = 0.001 (right). For µ = 0.01 we observe, that we get a sharper rear of the Figure 5: The optimized shape of the object using the volume formulation and α = 1 with µ = 0.01 (left) and µ = 0.001 (right). Compared to Fig. 4 we observe a sharper rear and for µ = 0.001 a symmetric airfoil shape emerges. For µ = 0.01 the drag is F D = 0.205542595 and the velocity inside the object is of order 10 −4 , which is two orders of magnitude smaller than in the case α = 0.03. For µ = 0.001 the drag is F D = 0.041090517 and the velocity inside the object is of order 10 −6 . object, while the magnitude of the velocity inside the object is of order 10 −4 , which is two orders of magnitudes smaller than in the case α = 0.03. We also mention that the shapes obtained here bear similarities to the optimized shape for the minimization of the dissipative energy, as presented in [13,Figs. 4 and 5]. For α = 1, and µ = 0.001 we observe a symmetric airfoil shape. Figure 6: The optimal shape for the maximization of the lift of an object, under a constraint on the dissipative power. We observe an inclined shape.
In Fig. 6 we show the resulting optimal shape of the object. As expected we observe an inclined structure in order to maximize lift, but due to the constraint on the potential power, the angle of attack is restricted. This is consistent with previous results in [14, Fig. 2].

Conclusion
In this paper, we formulate and analyze a phase field approximation for an abstract shape optimization problem subject to stationary Navier-Stokes flow with general objective functionals and integral state constraints. We provide examples for the objective functionals and integral constraints that are of practical relevance, and we establish the existence of minimizers, and derive the first order optimality conditions. A crucial point in the analysis is to show the existence of Lagrange multipliers corresponding to the integral constraints. In the general setting we assume that the Zowe-Kurcyusz constraint qualification holds, and verify these assumptions for two specific examples. The first involves integral constraints only in the variable ϕ, while the second involves the state variable u. The optimality conditions are solved using the VMPT method and several simulations are performed. We demonstrate that the proposed phase field approach can handle topology optimization, and compare the results of drag minimization with previous works. Lastly, we consider an example with an integral constraint involving the state variables, namely maximization of lift with constraint on the potential power. The optimal shapes obtained are consistent with previous works on the lift-to-drag ratio for fluid flow with small Reynolds number.