A phase field approach to shape optimization in Navier--Stokes flow with integral state constraint

TL;DR

This paper develops a phase field method for shape optimization in Navier--Stokes flow, incorporating integral state constraints and regularization, with numerical demonstrations on topology optimization and lift maximization.

Contribution

It extends existing shape optimization frameworks by integrating phase field, porous medium models, and multiple constraints, including center of mass, volume, and drag, with numerical validation.

Findings

Successful implementation of the phase field approach for shape optimization.

Numerical results demonstrate effectiveness in topology optimization and lift maximization.

Comparison shows advantages over previous drag minimization methods.

Abstract

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.

Figures

Tables

Table 1. Table 1: The diffuse ( F ε D superscript subscript 𝐹 𝜀 𝐷 F_{\varepsilon}^{D} ) and sharp ( F D superscript 𝐹 𝐷 F^{D} ) drag for the parameters from [ 14 ] and different values of ε 𝜀 \varepsilon . Note that α ε ​ ( − 1 ) → ∞ → subscript 𝛼 𝜀 1 \alpha_{\varepsilon}(-1)\to\infty for ε → 0 → 𝜀 0 \varepsilon\to 0 , i.e., the object becomes less permeable and thus the drag increases with ε → 0 → 𝜀 0 \varepsilon\to 0 . In [ 14 ] for ε = 0.00025 𝜀 0.00025 \varepsilon=0.00025 we observed F ε D = 3.9117 × 10 − 2 superscript subscript 𝐹 𝜀 𝐷 3.9117 superscript 10 2 F_{\varepsilon}^{D}=3.9117\times 10^{-2} and F D = 3.9499 × 10 − 2 superscript 𝐹 𝐷 3.9499 superscript 10 2 F^{D}=3.9499\times 10^{-2} .

$\epsilon$	0.008	0.004	0.002
$F_{\epsilon }^D$	$1.0570\times {10}^{-2}$	$1.9806\times {10}^{-2}$	$2.8370\times {10}^{-2}$
$F^D$	$1.1103\times {10}^{-2}$	$2.0519\times {10}^{-2}$	$2.9025\times {10}^{-2}$
$\epsilon$	0.001	0.0005	0.00025
$F_{\epsilon }^D$	$3.4255\times {10}^{-2}$	$3.8184\times {10}^{-2}$	$4.0739\times {10}^{-2}$
$F^D$	$3.4777\times {10}^{-2}$	$3.8572\times {10}^{-2}$	$4.1012\times {10}^{-2}$

Equations

$$-\mu \Delta \mathbf{u}+(\mathbf{u}\cdot \nabla )\mathbf{u}+\nabla p$$

$$\mathrm{div}\mathbf{u}$$

$$\mathbf{u}$$

$$\frac{{\int }_{\Gamma }{\mathbf{u}}_{\infty }^{\perp }\cdot \left(\mu \left(\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^{\top }\right)-p\mathrm{I}\right)\mathbf{\nu }\mathrm{d}{\mathcal{H}}^{d-1}}{{\int }_{\Gamma }{\mathbf{u}}_{\infty }\cdot \left(\mu \left(\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^{\top }\right)-p\mathrm{I}\right)\mathbf{\nu }\mathrm{d}{\mathcal{H}}^{d-1}},$$

$$\max {\int }_{\Gamma }{\mathbf{u}}_{\infty }^{\perp }\cdot \left(\mu \left(\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^{\top }\right)-p\mathrm{I}\right)\mathbf{\nu }\mathrm{d}{\mathcal{H}}^{d-1}$$

$$\text{ subject to eqrefSI:NS and }{\int }_{\Gamma }{\mathbf{u}}_{\infty }\cdot \left(\mu \left(\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^{\top }\right)-p\mathrm{I}\right)\mathbf{\nu }\mathrm{d}{\mathcal{H}}^{d-1}\le D,$$

$$\left|\mathrm{D}\phi \right|(\Omega )=2P_{\Omega }(\{\phi =1\}),$$

$$\underset{(\phi ,\mathbf{u},p)}{\min }{\mathcal{J}}_0(\phi ,\mathbf{u},p)$$

$$+{\int }_{\Omega }\frac{1}{2}h(x,\nabla \mathbf{u},p,{\mathbf{\nu }}_{\phi })\mathrm{d}\left|\mathrm{D}\phi \right|+\frac{\gamma }{2}\left|\mathrm{D}\phi \right|(\Omega ),$$

$$-\mu \Delta \mathbf{u}+(\mathbf{u}\cdot \nabla )\mathbf{u}+\nabla p$$

$$\mathrm{div}\mathbf{u}$$

$$\mathbf{u}$$

$$G_i(\phi ,\mathbf{u},p)=0\text{ for }1\le i\le m_1,G_i(\phi ,\mathbf{u},p)\ge 0\text{ for }{m}_1+1\le i\le m_1+m_2,$$

$$G_i(\phi ,\mathbf{u},p)$$

$$+{\int }_{\Omega }\frac{1}{2}{\mathbf{L}}_i(x,\nabla \mathbf{u},p)\cdot {\mathbf{\nu }}_{\phi }\mathrm{d}\left|\mathrm{D}\phi \right|,$$

$$\frac{1+\phi }{2}\left(\frac{\mu }{2}{\left|\nabla \mathbf{u}\right|}^2-\mathbf{f}\cdot \mathbf{u}\right),$$

$$\frac{1+\phi }{2}{\chi }_{\mathcal{Q}}(x)({\delta }_1{\left|p-p_{\mathrm{tar}}\right|}^2+{\delta }_2{\left|\mathbf{u}-{\mathbf{u}}_{\mathrm{tar}}\right|}^2)$$

$$\mathbf{a}\cdot \left(\mu \left(\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^{\top }\right)-p\mathrm{I}\right){\mathbf{\nu }}_{\phi },$$

$$\displaystyle\color[rgb]{0,0,0}G_{1}(\varphi)=\int_{\Omega}\varphi-\beta_{1}\,\mathrm{dx}\,\geq 0,\quad G_{2}(\varphi)=\int_{\Omega}\beta_{2}-\varphi\,\mathrm{dx}\,\geq 0,\color[rgb]{0,0,0}$$

$$\frac{{\beta }_1+1}{2}\left|\Omega \right|\le {\int }_{\Omega }\frac{1+\phi }{2}\mathrm{dx}=\left|E\right|\le \frac{{\beta }_2+1}{2}\left|\Omega \right|$$

$$\iff$$

$$\displaystyle\color[rgb]{0,0,0}G_{1}(\varphi)=M-\int_{\Omega}\tfrac{1}{2}\rho(x)(1-\varphi)\,\mathrm{dx}\,\geq 0\Leftrightarrow\color[rgb]{0,0,0}\int_{\Omega}\tfrac{1}{2}\rho(x)(1-\varphi)\,\mathrm{dx}\,\leq M,$$

$$G_i(\phi )={\int }_{\Omega }\frac{1}{2}(1-\phi )(x_i-y_i)\mathrm{dx}=0\text{ for }i=1,2,\dots ,d,$$

$$G_1(\phi ,\mathbf{u},p)=D-\mathbf{a}\cdot {\int }_{\Omega }\left(\mu (\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^{\top }){\mathbf{\nu }}_{\phi }-p{\mathbf{\nu }}_{\phi }\right)\frac{1}{2}\mathrm{d}\left|\mathrm{D}\phi \right|\ge 0$$

$$\iff \mathbf{a}\cdot {\int }_{\Omega }\left(\mu (\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^{\top }){\mathbf{\nu }}_{\phi }-p{\mathbf{\nu }}_{\phi }\right)\frac{1}{2}\mathrm{d}\left|\mathrm{D}\phi \right|\le D.$$

$$\frac{1}{2c_0}{\mathcal{E}}_{\epsilon }(\phi )=\frac{1}{2c_0}{\int }_{\Omega }\frac{\epsilon }{2}{\left|\nabla \phi \right|}^2+\frac{1}{\epsilon }\Psi (\phi )\mathrm{dx},$$

$$c_0:=\frac{1}{2}{\int }_{-1}^1\sqrt{2\Psi (s)}\mathrm{ds},$$

$$\overline{{\alpha }_{\epsilon }}\to \infty \text{ as }\epsilon \to 0.$$

$${\alpha }_{\epsilon }(\phi )\mathbf{u}-\mu \Delta \mathbf{u}+(\mathbf{u}\cdot \nabla )\mathbf{u}+\nabla p$$

$$\mathrm{div}\mathbf{u}$$

$$\mathbf{u}$$

$${\int }_{\Omega }b(x,\mathbf{u},\nabla \mathbf{u},p,\phi )\mathrm{dx}\text{ and }{\int }_{\Omega }{K}_i(x,\mathbf{u},\nabla \mathbf{u},p,\phi )\mathrm{dx}$$

$${\int }_{\Omega }\frac{1}{2}h(x,\nabla \mathbf{u},p,{\mathbf{\nu }}_{\phi })\mathrm{d}\left|\mathrm{D}\phi \right|\text{ and }{\int }_{\Omega }\frac{1}{2}{\mathbf{L}}_i(x,\nabla \mathbf{u},p)\cdot {\mathbf{\nu }}_{\phi }\mathrm{d}\left|\mathrm{D}\phi \right|$$

$${\int }_{\Omega }\frac{1}{2}h(x,\nabla \mathbf{u},p,\nabla \phi )\mathrm{dx}$$

$${\int }_{\Omega }\frac{1}{2}{\mathbf{L}}_i(x,\nabla \mathbf{u},p)\cdot \nabla \phi \mathrm{dx}$$

Full text

Key words. Topology optimization, shape optimization, phase field approach, Navier–Stokes flow, integral state constraints

AMS subject classification. 35Q35, 35Q56, 35R35, 49Q10, 49Q12, 65M22, 65M60, 76S05

1 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 $\Omega\subset\mathbb{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 $\Gamma:=\partial B\cap\Omega$ with the outer unit normal $\bm{\nu}$, and assume that $\Gamma\cap\partial\Omega=\emptyset$, i.e., the object $B$ never touches the external boundary. A fluid is present in the complement region $E:=\Omega\setminus B$, and we assume that the velocity $\bm{u}$ and the pressure $p$ of the fluid in the region $E$ obey the stationary Navier–Stokes equations with no-slip conditions on $\Gamma$, namely,

[TABLE]

Here $\bm{f}$ denotes the external body force, $\mu$ denotes the (constant) viscosity, and $\bm{g}$ models the inflow and outflow on the boundary $\partial\Omega$ such that $\int_{\partial\Omega}\bm{g}\cdot\bm{\nu}_{\partial\Omega}\,\mathrm{d}\mathcal{H}^{d-1}\,=0$, where $\bm{\nu}_{\partial\Omega}$ denotes the outer unit normal on $\partial\Omega$.

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

[TABLE]

where $\bm{u}_{\infty}$ is the flow direction, $\bm{u}_{\infty}^{\perp}$ is the perpendicular vector and $\mathcal{H}^{d-1}$ is the Hausdorff measure on the set $\Gamma$. 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.,

[TABLE]

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 $\varphi:\Omega\to\{\pm 1\}$, where $\{\varphi=1\}=E$ describes the fluid region and $\{\varphi=-1\}=B$ is its complement. The natural function space for the design functions is the space of bounded variations that take values $\pm 1$, i.e., $\varphi\in BV(\Omega,\{\pm 1\})$, which implies that the fluid region $E$ has finite perimeter $P_{\Omega}(E)$. If $\varphi$ is a function of bounded variation, its distributional derivative $\mathrm{D}\varphi$ is a finite Radon measure which can be decomposed into a positive measure $\left|\,\mathrm{D}\varphi\,\right|$ and a $S^{d-1}$-valued function $\bm{\nu}_{\varphi}\in L^{1}(\Omega,\left|\,\mathrm{D}\varphi\,\right|)^{d}$, where $S^{d-1}$ denotes the $(d-1)$-dimensional sphere. The total variation for $\varphi\in BV(\Omega,\{\pm 1\})$, denoted by $\left|\,\mathrm{D}\varphi\,\right|(\Omega)$, satisfies

[TABLE]

and thus we can express the Hausdorff measure $\mathcal{H}^{d-1}$ on the set $\Gamma$ as $\frac{1}{2}\left|\,\mathrm{D}\varphi\,\right|(\Omega)$. Furthermore, the $S^{d-1}$-valued function $\bm{\nu}_{\varphi}$ can be considered as a generalized normal on the set $\partial\{\varphi=1\}$ (see [1, 10, 16] for a more detailed introduction to the theory of sets of finite perimeter and functions of bounded variation).

For functions $b:\Omega\times\mathbb{R}^{d}\times\mathbb{R}^{d\times d}\times\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ and $h:\Omega\times\mathbb{R}^{d\times d}\times\mathbb{R}\times\mathbb{R}^{d}\to\mathbb{R}$, we consider the following general shape optimization problem with perimeter regularization:

[TABLE]

subject to $\varphi\in BV(\Omega,\{\pm 1\})$ and $(\bm{u},p)\in\bm{H}^{1}(E)\times L^{2}(E)$ fulfilling

[TABLE]

In addition, for fixed $m_{1},m_{2}\in\mathbb{N}\cup\{0\}$, we impose the $m_{1}$ integral equality constraints and $m_{2}$ integral inequality constraints:

[TABLE]

where for each $1\leq i\leq m_{1}+m_{2}$,

[TABLE]

for functions $K_{i}:\Omega\times\mathbb{R}^{d}\times\mathbb{R}^{d\times d}\times\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ and $\bm{L}_{i}:\Omega\times\mathbb{R}^{d\times d}\times\mathbb{R}\to\mathbb{R}^{d}$. The parameter $\gamma>0$ in (1.3) is the weighting factor for the perimeter regularization, which is given by the term $\frac{1}{2}\left|\,\mathrm{D}\varphi\,\right|$, and in light of the above discussion regarding the measure $\frac{1}{2}\left|\,\mathrm{D}\varphi\,\right|$ representing the Hausdorff measure on $\Gamma$, we see that the functions $b$ and $\{K_{i}\}_{i=1}^{m_{1}+m_{2}}$ model objectives and constraints in the bulk phases $E$ and $B$, while $h$ and $\{\bm{L}_{i}\cdot\bm{\nu}_{\varphi}\}_{i=1}^{m_{1}+m_{2}}$ model constraints on the interface $\Gamma$. Examples of $b$, $h$, $K_{i}$ and $\bm{L}_{i}$ are given below, and it is noteworthy to point out that there is no dependence on $\bm{u}$ in $\bm{L}_{i}$ as the no-slip condition (1.4d) ensures that $\bm{u}=\bm{0}$ on $\Gamma$. However, the gradient $\nabla\bm{u}$ may not vanish on $\Gamma$, which leads to its appearance in the surface constraints.

The appearance of the perimeter regularization $\frac{\gamma}{2}\left|\,\mathrm{D}\varphi\,\right|(\Omega)$ in (1.3) is motivated from the well-known difficulties regarding the mathematical treatment of shape optimization - in 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 $\bm{L}$ (where we neglect the index $i$ for convenience) that are of relevance. For a subset $A\subset\Omega$, we use the notation $\chi_{A}(x)$ to denote the characteristic function of $A$, i.e., $\chi_{A}(x)=1$ if $x\in A$ and $\chi_{A}(x)=0$ if $x\in\Omega\setminus A$. In particular, one can think of the design function $\varphi$ as $\varphi(x)=-1+2\chi_{E}(x)$ which satisfies $\varphi(x)=1$ for $x\in E$ and $\varphi(x)=-1$ for $x\in\Omega\setminus E=B$. Hence, in the following examples for the function $b$, we can use $\frac{1}{2}(1+\varphi)$ as a restriction to the region $E$ and similarly, $\frac{1}{2}(1-\varphi)$ as a restriction to the region $B$:

• •

the total potential power of the fluid

[TABLE]

• •

the construction cost of the object $\frac{1-\varphi}{2}w(x)$, where $w$ denotes a cost function per unit volume,

• •

the least square approximation

[TABLE]

to a target velocity profile $\bm{u}_{\mathrm{tar}}$ and a target pressure profile $p_{\mathrm{tar}}$ in an observation region $\mathcal{Q}\subset E$. Here $\delta_{1}$ and $\delta_{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 $\bm{a}$, which is given as

[TABLE]

where $\,\bm{\mathrm{I}}\,$ denotes the identity tensor. The drag of the object is given when $\bm{a}$ is parallel to the flow direction $\bm{u}_{\infty}$, while the lift of the object is given when $\bm{a}=\bm{u}_{\infty}^{\perp}$, 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}=\varphi-\beta_{1}$, $\bm{L}_{1}=\bm{0}$, $K_{2}=\beta_{2}-\varphi$ and $\bm{L}_{2}=\bm{0}$ for fixed constants $-1<\beta_{1}\leq\beta_{2}<1$ leads to inequality constraints:

[TABLE]

or equivalently

[TABLE]

• •

the prescribed mass of the object - setting $K_{1}=M\left|\,\Omega\,\right|^{-1}-\frac{1-\varphi}{2}\rho(x)$, $\bm{L}_{1}=\bm{0}$, where $\rho(x)$ is a mass density and $M>0$ is a target/maximal mass leads to the inequality constraint

[TABLE]

• •

the prescribed center of mass of the object (with uniform mass density) at a point $y$ in the interior of $\Omega$, i.e., $y\notin\partial\Omega$ - setting $K_{i}=\frac{1-\varphi}{2}(x_{i}-y_{i})$ and $\bm{L}_{i}=\bm{0}$ for $1\leq i\leq d$ leads to the equality constraints

[TABLE]

• •

the prescribed drag of the object - setting $\bm{L}_{1}=-\mu(\nabla\bm{u}+(\nabla\bm{u})^{\top})\bm{a}+p\bm{a}$, $K_{1}=D\left|\,\Omega\,\right|^{-1}$, where $\bm{a}$ is the unit vector parallel to the flow direction $\bm{u}_{\infty}$ and $D>0$ is a maximal drag value leads to the inequality constraint

[TABLE]

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 $\varphi$ 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 $\bm{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 $\bm{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.

2 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 $\varphi$ takes only values in $\{\pm 1\}$ (i.e., $\varphi\in BV(\Omega,\{\pm 1\})$) and now allow $\varphi$ to be a function with values in $\mathbb{R}$ and inherits $H^{1}(\Omega)$ regularity. In particular, we change the admissible space of design functions from subsets of $BV(\Omega,\{\pm 1\})$ to subsets of $H^{1}(\Omega)$. This leads to the development of interfacial layers $\{-1<\varphi<1\}$ in between the fluid region $E=\{\varphi=1\}$ and the object region $B=\{\varphi=-1\}$. This interfacial layer replaces the boundary $\Gamma$ of $B$ and a parameter $\varepsilon>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

[TABLE]

where $\Psi$ is a potential with equal minima at $\varphi=\pm 1$, to approximate the perimeter functional $P_{\Omega}$. The positive constant $c_{0}$ is dependent only on the potential $\Psi$ via the relation

[TABLE]

and it is well-known that $\frac{1}{2c_{0}}\mathcal{E}_{\varepsilon}$ approximates $\varphi\mapsto\frac{1}{2}\left|\,\mathrm{D}\varphi\,\right|(\Omega)=P_{\Omega}(\{\varphi=1\})$ in the sense of $\Gamma$-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 $(\overline{\alpha_{\varepsilon}})^{-1}\ll 1$. A function $\alpha_{\varepsilon}(\varphi)$ is introduced to interpolate between the inverse permeabilities of the fluid region $\alpha_{\varepsilon}(1)=0$ and the porous medium $\alpha_{\varepsilon}(-1)=\overline{\alpha_{\varepsilon}}$, which satisfies

[TABLE]

With this, we extend the state equations from $E$ to the whole domain $\Omega$ by the addition of the porous-medium term $\alpha_{\varepsilon}(\varphi)\bm{u}$:

[TABLE]

We note that this additional term vanishes in the fluid region, and in the limit $\varepsilon\to 0$, one expects the velocity $\bm{u}$ in the object region to vanish. We point out that in the modified state equations (2.2) we solve for a velocity field $\bm{u}$ and pressure field $p$ that are defined on the fixed domain $\Omega$. Furthermore, in the objective functional (1.3) and in the integral state constraints (1.6), the terms

[TABLE]

require no modification when we consider the phase field setting. For the surface terms such as

[TABLE]

arising from the objective functional (1.3) and the integral state constraints (1.6), we employ the idea in [14] to reformulate them. Assuming the function $h$ is one-homogeneous with respect to its last variable, which is true for the case of the hydrodynamic force (1.8), we use the vector-valued measure with density $\frac{1}{2}\nabla\varphi$ as an approximation to $\frac{1}{2}\bm{\nu}_{\varphi}\,\mathrm{d}\,\left|\,\mathrm{D}\varphi\,\right|$, see [14, §3.2] for more details. This in turn gives us the phase field approximation

[TABLE]

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

[TABLE]

is a phase field approximation of the surface integral constraint involving $\bm{L}_{i}$.

Recall that in the modified state equations (2.2) the porous-medium term $\alpha_{\varepsilon}(\varphi)\bm{u}$ serves to enforce the condition that the velocity $\bm{u}$ in the object region should vanish in the limit $\varepsilon\to 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

[TABLE]

a penalization term

[TABLE]

where $\hat{\alpha}_{\varepsilon}$ is a function with similar properties as $\alpha_{\varepsilon}$, i.e., $\hat{\alpha}_{\varepsilon}(1)=0$ and $\hat{\alpha}_{\varepsilon}(-1)\to\infty$ as $\varepsilon\to 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 $\varepsilon\to 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 $\hat{\alpha}_{\varepsilon}=0$, however in this case no rigorous results on the sharp interface limit $\varepsilon\to 0$ are known.

Before we state the phase field optimization problem, let us mention that for the analysis we assume that the function $\alpha_{\varepsilon}$ in the porous-medium term satisfies the properties:

1. $(\mathrm{A0})$

$\alpha_{\varepsilon}\in C^{1,1}(\mathbb{R})$ is non-negative, and there exist constants $s_{a},s_{b}\in\mathbb{R}$ with $s_{a}\leq-1$ and $s_{b}\geq 1$ such that

[TABLE]

In particular, for arbitrary $\phi$ and its truncation $\tilde{\phi}:=\max(s_{a},\min(s_{b},\phi))$ we see that $\alpha_{\varepsilon}(\phi)=\alpha_{\varepsilon}(\tilde{\phi})$, and so the state equations (2.2) for $\phi$ and $\tilde{\phi}$ are equivalent. Hence, without loss of generality, we now search for optimal design functions $\varphi$ exhibiting $H^{1}(\Omega)$-regularity and satisfies the pointwise bounds $s_{a}\leq\varphi\leq s_{b}$ a.e. in $\Omega$.

Taking into account the above discussions, we arrive at the following phase field approximation to the optimal control problem (1.3)-(1.6):

[TABLE]

subject to

[TABLE]

satisfying the weak formulation of (2.2):

[TABLE]

for all $\bm{v}\in\bm{H}^{1}_{0}(\Omega):=\{\bm{h}\in\bm{H}^{1}(\Omega)\,|\,\bm{h}=\bm{0}\text{ on }\partial\Omega\}$, along with $m_{1}$ equality and $m_{2}$ inequality integral constraints

[TABLE]

of the form

[TABLE]

for $1\leq i\leq m_{1}+m_{2}$.

2.1 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.

Lemma 2.1** ([14, Lem. 4.3]).**

Suppose $\alpha_{\varepsilon}\in C^{1,1}(\mathbb{R})$ is non-negative and satisfies (2.4), for every $\varphi\in L^{1}(\Omega)$ there exists at least one pair $(\bm{u},p)\in\bm{H}^{1}_{\bm{g},\sigma}(\Omega)\times L^{2}_{0}(\Omega)$ such that (2.6) is satisfied. Furthermore, there exists a positive constant $C=C(\mu,\alpha_{\varepsilon},\bm{f},\bm{g},\Omega)$ independent of $\varphi$ such that

[TABLE]

The estimate (2.8) can be established by testing the weak form (2.6) with $\bm{u}-\bm{G}$, where $\bm{G}\in\bm{H}^{1}_{\bm{g},\sigma}(\Omega)$ is a vector field depending on the inflow/outflow $\bm{g}$ and the domain $\Omega$, and satisfying certain properties. Furthermore, this computation also shows that the constant $C$ depends on $\varepsilon$ only through the function $\alpha_{\varepsilon}$.

By the above existence result, we can define a set-valued solution operator

[TABLE]

for any $\varphi\in L^{1}(\Omega)$. Next, we state a continuity property of the solution operator.

Lemma 2.2** ([14, Lem. 4.4 and 4.5]).**

Let $(\varphi_{k})_{k\in\mathbb{N}}\subset L^{1}(\Omega)$ be a sequence with corresponding solution $(\bm{u}_{k},p_{k})\in\bm{S}_{\varepsilon}(\varphi_{k})\subset\bm{H}^{1}(\Omega)\times L^{2}(\Omega)$ for each $k\in\mathbb{N}$. Suppose there exists $\varphi\in L^{1}(\Omega)$ such that

[TABLE]

Then, there exists a subsequence, denoted by the same index, and functions $\bm{u}\in\bm{H}^{1}(\Omega)$, $p\in L^{2}(\Omega)$ such that

[TABLE]

with the property that $(\bm{u},p)\in\bm{S}_{\varepsilon}(\varphi)$. Furthermore, it holds that

[TABLE]

In general, we do not have uniqueness of solutions to (2.6), but there is a conditional uniqueness result.

Lemma 2.3** ([11, Lem. 5], [18, Lem. 12.2]).**

If there exists $\bm{u}\in\bm{S}_{\varepsilon}(\varphi)$ with

[TABLE]

where

[TABLE]

Then, $\bm{S}_{\varepsilon}(\varphi)=\{(\bm{u},p)\}$. That is, there is exactly one solution of (2.6) corresponding to $\varphi\in L^{1}(\Omega)$.

The additional assumption (2.10) on the solution $\bm{u}\in\bm{S}_{\varepsilon}(\varphi)$ to ensure uniqueness of the state equations can be achieved for small data $\bm{f}$ and $\bm{g}$ or with high viscosity $\mu$. 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 $\bm{S}_{\varepsilon}$ and the derivation of the optimality conditions, we require that $\bm{S}_{\varepsilon}$ is a one-to-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 $\bm{S}_{\varepsilon}$.

Lemma 2.4** ([14, Lem. 4.8]).**

Let $\varphi_{\varepsilon}\in H^{1}(\Omega)\,\cap L^{\infty}(\Omega)$ be given such that $\bm{S}_{\varepsilon}(\varphi_{\varepsilon})=\{(\bm{u}_{\varepsilon},p_{\varepsilon})\}$. Then, there exists a neighborhood $N$ of $\varphi_{\varepsilon}$ in $H^{1}(\Omega)\cap L^{\infty}(\Omega)$ such that for every $\delta\in N$, the solution operator consists of exactly one pair, and we may write $\bm{S}_{\varepsilon}:N\to\bm{H}^{1}(\Omega)\times L^{2}(\Omega)$. This mapping is differentiable at $\varphi_{\varepsilon}$ with

[TABLE]

where $(\bm{w}_{\varepsilon},r_{\varepsilon})$ is the unique solution to

[TABLE]

which is the weak formulation of the following the linearized state system

[TABLE]

3 Existence of a minimizer

We make the following assumptions for the potential $\Psi$ and the functions $\hat{\alpha}_{\varepsilon}$, $b$, $h$, $K_{i}$, and $\bm{L}_{i}$.

1. $(\mathrm{A1})$

Let $\Psi\in C^{1,1}(\mathbb{R})$ be a non-negative function such that $\Psi(s)=0$ if and only if $s=\pm 1$, and that there exist positive constants $c_{1}$, $c_{2}$, $t_{0}$ such that

[TABLE] 2. $(\mathrm{A2})$

The function $\hat{\alpha}\in C^{1,1}(\mathbb{R})$ satisfies the same assumptions as $\alpha_{\varepsilon}$, i.e., $\hat{\alpha}_{\varepsilon}$ is non-negative, with $\hat{\alpha}_{\varepsilon}(1)=0$, $\hat{\alpha}_{\varepsilon}(-1)\to\infty$ as $\varepsilon\to 0$, and fulfills (2.4). 3. $(\mathrm{A3})$

The function $b:\Omega\times\mathbb{R}^{d}\times\mathbb{R}^{d\times d}\times\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ is a Carathéodory function of the form

[TABLE]

for some Carathéodory functions $B:\Omega\times\mathbb{R}^{d}\times\mathbb{R}^{d\times d}\times\mathbb{R}\to\mathbb{R}$, $z:\Omega\times\mathbb{R}\to\mathbb{R}$ and there exist non-negative functions $f_{b}\in L^{1}(\Omega)$, $\{f_{b,i}\}_{i=1}^{4}\subset L^{\infty}(\Omega)$ such that for a.e. $x\in\Omega$ it holds for any $r\geq 0$, $p\geq 2$ in two-dimensions and $2\leq p\leq 6$ in three-dimensions,

[TABLE]

for all $s,t\in\mathbb{R}$, $\bm{w}\in\mathbb{R}^{d}$ and $\bm{A}\in\mathbb{R}^{d\times d}$. 4. $(\mathrm{A4})$

The function $h:\Omega\times\mathbb{R}^{d\times d}\times\mathbb{R}\times\mathbb{R}^{d}\to\mathbb{R}$ is a Carathéodory function that is one-homogeneous with respect to its last variable and there exist non-negative functions $f_{h}\in L^{1}(\Omega)$, $\{f_{h,k}\}_{k=1}^{3}\subset L^{\infty}(\Omega)$ such that for a.e. $x\in\Omega$ it holds,

[TABLE]

for all $s,t\in\mathbb{R}$, $\bm{w}\in\mathbb{R}^{d}$ and $\bm{A}\in\mathbb{R}^{d\times d}$. 5. $(\mathrm{A5})$

For each $1\leq i\leq m_{1}$, the function $K_{i}:\Omega\times\mathbb{R}^{d}\times\mathbb{R}^{d\times d}\times\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ is a Carathéodory function of the form

[TABLE]

for functions $k_{i}\in L^{1}(\Omega)$ and Carathéodory functions $\mathcal{K}_{i}:\Omega\times\mathbb{R}^{d}\times\mathbb{R}^{d\times d}\times\mathbb{R}\to\mathbb{R}$, $y_{i}:\Omega\times\mathbb{R}\to\mathbb{R}$ and there exist non-negative functions $z_{1}\in L^{1}(\Omega)$, $z_{2},z_{3},z_{4},z_{5}\in L^{\infty}(\Omega)$ such that for a.e. $x\in\Omega$ it holds for any $r\geq 0$, $p\in[2,\infty)$ in two-dimensions and $p\in[2,6)$ in three-dimensions,

[TABLE]

for all $s,t\in\mathbb{R}$, $\bm{w}\in\mathbb{R}^{d}$ and $\bm{A}\in\mathbb{R}^{d\times d}$. 6. $(\mathrm{A6})$

For each $1\leq j\leq m_{2}$, the function $\bm{L}_{j}:\Omega\times\mathbb{R}^{d\times d}\times\mathbb{R}\to\mathbb{R}^{d}$ is a Carathéodory function and there exist non-negative functions $z_{6}\in L^{1}(\Omega)$, $z_{7},z_{8}\in L^{\infty}(\Omega)$ such that for a.e. $x\in\Omega$ it holds,

[TABLE]

for all $s\in\mathbb{R}$ and $\bm{A}\in\mathbb{R}^{d\times d}$. 7. $(\mathrm{A7})$

We assume that the set

[TABLE]

is non-empty, where the phase field integral state constraints $G_{k}$, for $1\leq k\leq m_{1}+m_{2}$, are of the form (2.7) and we recall the set $\Phi$ is defined as $\{f\in H^{1}(\Omega)\,|\,s_{a}\leq f\leq s_{b}\text{ a.e. in }\Omega\}$. 8. $(\mathrm{A8})$

The functionals $\mathcal{B}:H^{1}(\Omega)\times\bm{H}^{1}(\Omega)\times L^{2}(\Omega)\to\mathbb{R}$ and $\mathcal{H}:H^{1}(\Omega)\times\bm{H}^{1}(\Omega)\times L^{2}(\Omega)\to\mathbb{R}$ defined as

[TABLE]

satisfy $\mathcal{B}|_{\mathbb{K}_{ad}\times\bm{H}^{1}_{\bm{g},\sigma}(\Omega)\times L^{2}_{0}(\Omega)}$ and $\mathcal{H}|_{\mathbb{K}_{ad}\times\bm{H}^{1}_{\bm{g},\sigma}(\Omega)\times L^{2}_{0}(\Omega)}$ are bounded from below, $\mathcal{B}$ is weakly lower semicontinuous, and for all $\varphi_{n}\rightharpoonup\varphi$ in $H^{1}(\Omega)$, $\bm{u}_{n}\to\bm{u}$ in $\bm{H}^{1}(\Omega)$, $p_{n}\to p$ in $L^{2}(\Omega)$, it holds that

[TABLE]

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 $\frac{1+\varphi}{2}$, and the function $k_{i}$ would be of the form $D\left|\,\Omega\,\right|^{-1}$. Furthermore, the set $\mathbb{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 $\mathbb{K}_{ad}$ for the general setting here, later in Sec. 5 we show for two examples that the corresponding set $\mathbb{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).

Lemma 3.1**.**

Under ($\mathrm{A5}$) and ($\mathrm{A6}$), let $\{\varphi_{n}\}_{n\in\mathbb{N}}$ be a sequence in $\mathbb{K}_{ad}$ such that $\varphi_{n}\rightharpoonup\varphi\in H^{1}(\Omega)$ for some $\varphi\in H^{1}(\Omega)$, then $\varphi\in\mathbb{K}_{ad}$.

Proof.

Let $\{\varphi_{n}\}_{n\in\mathbb{N}}$ be a sequence in $\mathbb{K}_{ad}$ with weak limit $\varphi\in H^{1}(\Omega)$. It suffices to show that if $(\bm{u},p)\in\bm{S}_{\varepsilon}(\varphi)$, then $G_{i}(\varphi,\bm{u},p)=0$ for $1\leq i\leq m_{1}$ and $G_{m_{1}+j}(\varphi,\bm{u},p)\geq 0$ for $1\leq j\leq m_{2}$, which then implies that $\varphi\in\mathbb{K}_{ad}$.

Let $\{(\bm{u}_{n},p_{n})\}_{n\in\mathbb{N}}\subset\bm{H}^{1}_{\bm{g},\sigma}(\Omega)\times L^{2}_{0}(\Omega)$ be the corresponding solutions to (2.2) for $\varphi_{n}$, i.e., for each $n\in\mathbb{N}$, $(\bm{u}_{n},p_{n})\in\bm{S}_{\varepsilon}(\varphi_{n})$. Since $\varphi_{n}\rightharpoonup\varphi\in H^{1}(\Omega)$, by compactness we have strong convergence along subsequences $\varphi_{n_{j}}\to\varphi$ in $L^{p}(\Omega)$ for $p\in[1,\infty)$ in two dimensions and $p\in[1,6)$ in three dimensions. Consequently, we also have $\varphi_{n_{j}}\to\varphi$ a.e. in $\Omega$ and hence $s_{a}\leq\varphi\leq s_{b}$ a.e. in $\Omega$. Furthermore, by the assertions of Lem. 2.2, the corresponding solutions $\{(\bm{u}_{n_{j}},p_{n_{j}})\}_{j\in\mathbb{N}}$ satisfy $\bm{u}_{n_{j}}\to\bm{u}$ in $\bm{H}^{1}(\Omega)$ and $p_{n_{j}}\to p$ in $L^{2}(\Omega)$ where $(\bm{u},p)\in\bm{S}_{\varepsilon}(\varphi)$.

For each $1\leq i\leq m_{1}+m_{2}$, by the continuity of $\bm{L}_{i}$ with respect to its second and third variables, it holds that $\bm{L}_{i}(x,\nabla\bm{u}_{n_{j}},p_{n_{j}})\to\bm{L}_{i}(x,\nabla\bm{u},p)$ a.e. in $\Omega$. Using the growth conditions in ($\mathrm{A6}$), the strong convergences for $\{\bm{u}_{n_{j}},p_{n_{j}}\}_{j\in\mathbb{N}}$ and the generalized Lebesgue dominated convergence theorem leads to

[TABLE]

Together with the weak convergence $\nabla\varphi_{n_{j}}$ to $\nabla\varphi$ in $\bm{L}^{2}(\Omega)$, we have

[TABLE]

Note that $s_{a}\leq\varphi_{n_{j}},\varphi\leq s_{b}$ a.e. in $\Omega$ for all $j\in\mathbb{N}$, and thus there exists a constant $M>0$ such that $\sup_{x\in\Omega}\left(\left|\,y_{i}(x,\varphi_{n_{j}})\,\right|,\left|\,y_{i}(x,\varphi)\,\right|\right)\leq M$ for all $n\in\mathbb{N}$. Using the splitting

[TABLE]

we can show that $\lim_{n\to\infty}G_{i}(\varphi_{n_{j}},\bm{u}_{n_{j}},p_{n_{j}})=G_{i}(\varphi,\bm{u},p)$ once we demonstrate that $I_{1},I_{2}\to 0$ as $n\to\infty$. This would then imply that $\varphi\in\mathbb{K}_{ad}$. Using the growth conditions in ($\mathrm{A5}$) for $\mathcal{K}_{i}$, the strong convergences for $\{(\bm{u}_{n_{j}},p_{n_{j}})\}_{j\in\mathbb{N}}$ and the generalized Lebesgue dominated convergence theorem yields that

[TABLE]

Then, the assertion that $I_{1}\to 0$ as $j\to\infty$ follows from the above strong convergence in $L^{1}(\Omega)$ and the boundedness of $y_{i}(x,\varphi_{n_{j}})$ in $L^{\infty}(\Omega)$. Meanwhile, dominating the sequence $\{\mathcal{K}_{i}(x,\bm{u},\nabla\bm{u},p)\;y_{i}(x,\varphi_{n_{j}})\}_{j\in\mathbb{N}}$ by the function $\|z_{5}\|_{L^{\infty}(\Omega)}M\left|\,\mathcal{K}_{i}(x,\bm{u},\nabla\bm{u},p)\,\right|\in L^{1}(\Omega)$, and the application of the usual Lebesgue dominating convergence theorem yields

[TABLE]

and hence $I_{2}\to 0$ as $n\to\infty$. ∎∎

We state the existence result for a minimizer of the problem (2.5)-(2.7).

Theorem 3.2**.**

Under Assumptions ($\mathrm{A1}$)-($\mathrm{A8}$), there exists at least one minimizer to the problem (2.5)-(2.7).

Proof.

By ($\mathrm{A8}$), $(\mathcal{B}+\mathcal{H})|_{\mathbb{K}_{ad}\times\bm{H}^{1}_{\bm{g},\sigma}(\Omega)\times L^{2}_{0}(\Omega)}$ is bounded from below by a constant $C_{0}\in\mathbb{R}$. Then, by the non-negativity of $\hat{\alpha}_{\varepsilon}$ and $\Psi$, we find that there exists a constant $C_{1}\in\mathbb{R}$ such that $\mathcal{J}_{\varepsilon}:\mathbb{K}_{ad}\times\bm{H}^{1}_{\bm{g},\sigma}(\Omega)\times L^{2}_{0}(\Omega)\to\mathbb{R}$ is bounded from below by $C_{1}$. Thus, we can choose a minimizing sequence $(\varphi_{n},\bm{u}_{n},p_{n})_{n\in\mathbb{N}}\subset\mathbb{K}_{ad}\times\bm{H}^{1}_{\bm{g},\sigma}(\Omega)\times L^{2}_{0}(\Omega)$ such that $(\bm{u}_{n},p_{n})\in\bm{S}_{\varepsilon}(\varphi_{n})$ for all $n\in\mathbb{N}$ and

[TABLE]

Then, for arbitrary $\eta>0$, there exists $N\in\mathbb{N}$ such that for $n>N$,

[TABLE]

The above estimate implies that $\{\varphi_{n}\}_{n\in\mathbb{N}}\subset\mathbb{K}_{ad}$ is bounded uniformly in $H^{1}(\Omega)\cap L^{\infty}(\Omega)$. Thus, we may choose a subsequence $(\varphi_{n_{k}})_{k\in\mathbb{N}}$ such that $\varphi_{n_{k}}\to\varphi$ strongly in $L^{p}(\Omega)$ and almost everywhere in $\Omega$ for $2\leq p<\infty$ in two-dimensions and $2\leq p<6$ in three-dimensions. Furthermore, by Lem. 3.1 we also have that $\varphi\in\mathbb{K}_{ad}$, and by Lem. 2.2, there is a subsequence $(\bm{u}_{n_{k}},p_{n_{k}})_{n\in\mathbb{N}}\subset\bm{H}^{1}_{\bm{g},\sigma}(\Omega)\times L^{2}_{0}(\Omega)$ such that

[TABLE]

for some $(\bm{u},p)\in\bm{S}_{\varepsilon}(\varphi)$, and

[TABLE]

The continuity of $\Psi$ together with the fact that $(\varphi_{n_{k}})_{k\in\mathbb{N}}\subset L^{\infty}(\Omega)$ implies $(\Psi(\varphi_{n_{k}}))_{k\in\mathbb{N}}$ is a bounded sequence in $L^{\infty}(\Omega)$. The application of the dominated convergence theorem yields that $\Psi(\varphi_{n_{k}})$ converges strongly to $\Psi(\varphi)$ in $L^{1}(\Omega)$ as $k\to\infty$. Furthermore, by the weak lower semicontinuity assumptions of $\mathcal{B}$ and $\mathcal{H}$, and the weak lower semicontinuity of the mapping $\varphi\mapsto\|\nabla\varphi\|_{\bm{L}^{2}(\Omega)}^{2}$, we find that

[TABLE]

and so $(\varphi,\bm{u},p)\in\mathbb{K}_{ad}\times\bm{H}^{1}_{\bm{g},\sigma}(\Omega)\times L^{2}_{0}(\Omega)$ is a minimizer of (2.5)-(2.7). ∎∎

From this point onwards, for fixed $\varepsilon>0$, we denote a minimizer to the optimal control problem (2.5)-(2.7) as $\varphi_{\varepsilon}$ with corresponding unique solution $(\bm{u}_{\varepsilon},p_{\varepsilon})$ to the state equation (2.6).

4 Optimality conditions

We use the notation $\mathrm{D}_{j}f$ to denote the partial derivative of $f$ with respect to its $j$th variable. Furthermore, the notation $\left|\,\mathrm{D}_{(i,j)}f\,\right|\leq P$ means that the partial derivatives $\mathrm{D}_{i}f$ and $\mathrm{D}_{j}f$ satisfy $\left|\,\mathrm{D}_{i}f\,\right|\leq P$ and $\left|\,\mathrm{D}_{j}f\,\right|\leq P$. To obtain optimality conditions, we make the following assumptions on the differentiability of $B$, $z$, $h$, $\mathcal{K}_{i}$, $y_{i}$, and $\bm{L}_{i}$.

1. $(\mathrm{B1})$

In addition to ($\mathrm{A3}$) assume further that $x\mapsto B(x,\bm{w},\bm{A},s)$, $x\mapsto z(x,t)$ and $x\mapsto h(x,\bm{A},s,\bm{w})$ belong to $W^{1,1}(\Omega)$ for all $\bm{w}\in\mathbb{R}^{d}$, $\bm{A}\in\mathbb{R}^{d\times d}$, $s,t\in\mathbb{R}$, and the partial derivatives

[TABLE]

exist for all $\bm{w}\in\mathbb{R}^{d}$, $s\in\mathbb{R}$, $\bm{A}\in\mathbb{R}^{d\times d}$, and a.e. $x\in\Omega$ as Carathéodory functions with

[TABLE]

for some non-negative functions $\tilde{a}\in L^{2}(\Omega)$, $\tilde{c}\in L^{\frac{p}{p-1}}(\Omega)$, $\tilde{b}_{1},\tilde{b}_{2},\tilde{b}_{3}\in L^{\infty}(\Omega)$, where $p\geq 2$ in two dimensions and $p\in[2,6]$ in three dimensions. 2. $(\mathrm{B2})$

For each $1\leq i\leq m_{1}+m_{2}$, in addition to ($\mathrm{A5}$) assume further that $x\mapsto\mathcal{K}_{i}(x,\bm{w},\bm{A},s)$, $x\mapsto k_{i}(x)$, $x\mapsto y_{i}(x,t)$, and $x\mapsto\bm{L}_{i}(x,\bm{A},s)$ belong to $W^{1,1}(\Omega)$ for all $\bm{w}\in\mathbb{R}^{d}$, $\bm{A}\in\mathbb{R}^{d\times d}$, $s,t\in\mathbb{R}$ and the partial derivatives

[TABLE]

exist for all $\bm{w}\in\mathbb{R}^{d}$, $s\in\mathbb{R}$, $\bm{A}\in\mathbb{R}^{d\times d}$, and a.e. $x\in\Omega$ as Carathéodory functions. Moreover, we assume that

[TABLE]

for some non-negative functions $\tilde{a}\in L^{2}(\Omega)$, $\tilde{c}\in L^{\frac{p}{p-1}}(\Omega)$ and $\tilde{b}_{1},\tilde{b}_{2},\tilde{b}_{3}\in L^{\infty}(\Omega)$, where $p\geq 2$ in two dimensions and $p\in[2,6]$ in three dimensions.

Under ($\mathrm{B1}$) and using [37, §4.3.3] or [17, Thm. 1 and 3], the Nemytskii operators

[TABLE]

are well-defined and the operator

[TABLE]

is continuously Fréchet differentiable (see [17, Thm. 7] or [37, §4.3.3] with $p=r=2$ and $q=1$). Hence, we find that

[TABLE]

is continuously Fréchet differentiable with derivative at $(\varphi_{\varepsilon},\bm{u}_{\varepsilon},p_{\varepsilon})$ in the direction $(\eta,\bm{v},s)$ given as

[TABLE]

Here we use the notation

[TABLE]

where the partial derivatives are evaluated at $(x,\nabla\bm{u}_{\varepsilon},p_{\varepsilon},\nabla\varphi_{\varepsilon})$. With a similar argument, the mappings

[TABLE]

for $1\leq i\leq m_{1}+m_{2}$, are continuously Fréchet differentiable, with derivatives at $(\varphi_{\varepsilon},\bm{u}_{\varepsilon},p_{\varepsilon})$ in the direction $(\eta,\bm{v},s)$ given as

[TABLE]

4.1 Fréchet differentiability of the objective functional

Due to the well-posedness of the state equations, we may now write the problem (2.5)-(2.7) as a minimizing problem for a reduced objective functional defined on an open set in $H^{1}(\Omega)\cap L^{\infty}(\Omega)$ with the help of Lem. 2.4. Let $(\varphi_{\varepsilon},\bm{u}_{\varepsilon},p_{\varepsilon})\in\mathbb{K}_{ad}\times\bm{H}^{1}_{\bm{g},\sigma}(\Omega)\times L^{2}_{0}(\Omega)$ denote a minimizer of (2.5)-(2.7), obtained from Thm. 3.2. By Lem. 2.4, there exists a neighborhood $N\subset H^{1}(\Omega)\cap L^{\infty}(\Omega)$ of $\varphi_{\varepsilon}$ such that for every $\psi\in N$, (2.6) is uniquely solvable. We define the reduced functional $j_{\varepsilon}:N\to\mathbb{R}$ by

[TABLE]

We now show that, as a mapping from \color[rgb]{0,0,0}N\subset\color[rgb]{0,0,0}H^{1}(\Omega)\cap L^{\infty}(\Omega)\to\mathbb{R}, $j_{\varepsilon}$ is Fréchet differentiable at $\varphi_{\varepsilon}$. As Lem. 2.4 guarantees the Fréchet differentiability of the solution operator $\bm{S}_{\varepsilon}(\varphi_{\varepsilon})$ as a mapping from \color[rgb]{0,0,0}N\color[rgb]{0,0,0} to $\bm{H}^{1}(\Omega)\times L^{2}(\Omega)$, we focus on the dependence of $\mathcal{J}_{\varepsilon}$ on the first variable.

Fix $\varphi\in H^{1}(\Omega)$, then by ($\mathrm{A2}$), $\hat{\alpha}_{\varepsilon}$ and $\hat{\alpha}^{\prime}_{\varepsilon}$ are uniformly bounded and so

[TABLE]

is a well-defined mapping from $H^{1}(\Omega)\subset L^{6}(\Omega)$ to $L^{6}(\Omega)$. By [37, §4.3.3], we see that $\hat{\alpha}_{\varepsilon}$ defines a Fréchet differentiable Nemytskii operator as a mapping from $L^{6}(\Omega)$ to $L^{3}(\Omega)$. Meanwhile, the assumption $\Psi\in C^{1,1}(\mathbb{R})$ and [37, Lem. 4.12] imply that $\Psi(\varphi)$ is continuously Fréchet differentiable Nemytskii operator as a mapping from $L^{\infty}(\Omega)$ to $L^{\infty}(\Omega)$. Combined with the Fréchet differentiability of the mapping $H^{1}(\Omega)\ni\varphi\mapsto\int_{\Omega}\left|\,\nabla\varphi\,\right|^{2}\,\mathrm{dx}\,$, $\mathcal{B}$ and $\mathcal{H}$, we obtain that $j_{\varepsilon}:N\to\mathbb{R}$ is Fréchet differentiable.

4.2 Existence of Lagrange multipliers

To show the existence of Lagrange multipliers for the integral constraints, we make use of the Zowe–Kurcyusz constraint qualification (ZKCQ), see [39] and [37, §6.1.2] for more details. For this purpose, we introduce the notation

[TABLE]

and recall the set

[TABLE]

Then, $\Phi$ is a closed convex subset of $H^{1}(\Omega)$ and $\mathbb{K}$ is a closed convex cone in $\mathbb{Y}$ with vertex at the origin, i.e., $\delta_{1}\mathbb{K}+\delta_{2}\mathbb{K}\subset\mathbb{K}$ for $\delta_{1},\delta_{2}>0$. In the notation of [39], we introduce the sets

[TABLE]

Fix $1\leq i\leq m_{1}+m_{2}$ and an arbitrary function $\zeta\in\Phi$. Convexity of $\Phi$ implies that $\varphi_{\varepsilon}+t(\zeta-\varphi_{\varepsilon})\in\Phi$ for sufficiently small values of $t$. Then, denoting the linearized state variables associated to $\delta=\zeta-\varphi_{\varepsilon}$ as $(\bm{w}_{\varepsilon},r_{\varepsilon})$ (see Lem. 2.4), the mapping $\phi\mapsto\mathcal{G}_{i}(\phi)=G_{i}(\phi,\bm{S}_{\varepsilon}(\phi))$ is continuously Fréchet differentiable at $\varphi_{\varepsilon}$ with derivative in the direction $\zeta-\varphi_{\varepsilon}$ given as (see also (4.4))

[TABLE]

where $\mathrm{D}_{2}K_{i}$, $\mathrm{D}_{3}K_{i}$, $\mathrm{D}_{4}K_{i}$ and $\mathrm{D}_{5}K_{i}$ are evaluated at $(x,\bm{u}_{\varepsilon},\nabla\bm{u}_{\varepsilon},p_{\varepsilon},\varphi_{\varepsilon})$. Then, it holds that

[TABLE]

The existence of bounded Lagrange multipliers $\bm{\lambda}:=(\lambda_{1},\dots,\lambda_{m_{1}+m_{2}})\in\mathbb{K}^{+}:=\{\bm{y}\in\mathbb{Y}\,|\,\bm{y}\cdot\bm{\eta}=0\;\forall\bm{\eta}\in\mathbb{K}\}$ satisfying

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the duality pairing between $H^{1}(\Omega)$ and its dual, follows if $\varphi_{\varepsilon}$ is a regular point in the sense of [39], or equivalently the so-called Zowe–Kurcyusz constraint qualification

[TABLE]

has to hold. We now make the following assumption:

1. $(\mathrm{C1})$

For any $\bm{z}\in\mathbb{Y}=\mathbb{R}^{m_{1}+m_{2}}$, there exists a function $\psi_{*}\in\Phi$, vectors $\bm{\tau}\in\mathbb{Y}$, $\bm{\xi},\bm{\eta}\in\mathbb{R}^{m_{2}}$ such that $\tau_{i}\geq 0$, $\xi_{j}\geq 0$, $\eta_{j}\geq 0$ for $1\leq i\leq m_{1}+m_{2}$ and $1\leq j\leq m_{2}$, and

[TABLE]

Then, under ($\mathrm{C1}$) and using [39, Thm. 3.1 and 4.1] there exist $\lambda_{1},\dots,\lambda_{m_{1}}\in\mathbb{R}$ and $\lambda_{m_{1}+1},\dots,\lambda_{m_{1}+m_{2}}\in\mathbb{R}_{\geq 0}$ such that

[TABLE]

holds with the following complementary slackness conditions for the inequality constraints

[TABLE]

We mention that (4.6) is equivalent (see [39, §3] and [21, Thm. 1.56]) to the following interior point/linearized Slater condition (which is also commonly known as the Robinson regularity condition [30]):

[TABLE]

4.3 Adjoint system

We now introduce the Lagrangian $\mathbb{L}:(H^{1}(\Omega)\cap L^{\infty}(\Omega))\times\bm{H}^{1}(\Omega)\times L^{2}(\Omega)\times\bm{H}^{1}_{0}(\Omega)\times L^{2}(\Omega)\to\mathbb{R}$ as

[TABLE]

where $\lambda_{i}$ is the Lagrange multiplier for the integral constraint $\mathcal{G}_{i}(\varphi)$ and $\theta$ is a Lagrange multiplier for the constraint $\int_{\Omega}p\,\mathrm{dx}\,=0$ for the pressure. A formal computation of $\mathrm{D}_{\bm{u}}\mathbb{L}$ and $\mathrm{D}_{p}\mathbb{L}$ yields the following adjoint system for the minimizer $\varphi_{\varepsilon}$:

[TABLE]

where $\mathrm{D}_{(2,3,4)}b$ are evaluated at $(x,\bm{u}_{\varepsilon},\nabla\bm{u}_{\varepsilon},p_{\varepsilon},\varphi_{\varepsilon})$, $\mathrm{D}_{(2,3)}h$ are evaluated at $(x,\nabla\bm{u}_{\varepsilon},p_{\varepsilon},\nabla\varphi_{\varepsilon})$, $\mathrm{D}_{(2,3,4)}K_{i}$ are evaluated at $(x,\bm{u}_{\varepsilon},\nabla\bm{u}_{\varepsilon},p_{\varepsilon},\varphi_{\varepsilon})$, and $\mathrm{D}_{(2,3)}\bm{L}_{i}$ are evaluated at $(x,\nabla\bm{u}_{\varepsilon},p_{\varepsilon})$, and upon integrating the divergence equation for $\bm{q}_{\varepsilon}$, we obtain

[TABLE]

Let us also recall from (3.1) and (3.2) that

[TABLE]

We now show that the adjoint system is well-posed.

Lemma 4.1**.**

Let $\varphi_{\varepsilon}\in H^{1}(\Omega)\cap L^{\infty}(\Omega)$ be the minimizer obtained from Thm. 3.2 and $(\bm{u}_{\varepsilon},p_{\varepsilon})=\bm{S}_{\varepsilon}(\varphi_{\varepsilon})$. Furthermore, let $\{\lambda_{i}\}_{i=1}^{m_{1}+m_{2}}$ be the Lagrange multipliers associated to the integral state constraints $\{\mathcal{G}_{i}(\varphi_{\varepsilon})\}_{i=1}^{m_{1}+m_{2}}$. Then, under ($\mathrm{B1}$), ($\mathrm{B2}$) and ($\mathrm{C1}$), there exists a unique weak solution pair $(\bm{q}_{\varepsilon},\pi_{\varepsilon})\in\bm{H}^{1}_{0}(\Omega)\times L^{2}(\Omega)$ to the adjoint system (4.9) in the following sense

[TABLE]

for all $\bm{v}\in\bm{H}^{1}_{0,\sigma}(\Omega)$.

Proof.

For convenience, we use the notation

[TABLE]

so that (4.9b) reads as $\,\mathrm{div}\,\bm{q}_{\varepsilon}=g$ in $\Omega$. Then, by ($\mathrm{B1}$), ($\mathrm{B2}$) and (4.10), we see that $g$ belongs to the function space $L^{2}_{0}(\Omega)$.

Applying [33, Lem. II.2.1.1], we find a vector field $\bm{G}\in\bm{H}^{1}_{0}(\Omega)$ such that

[TABLE]

for some constant $C>0$ depending only on $\Omega$. We define the bilinear form $a:\bm{H}^{1}_{0,\sigma}(\Omega)\times\bm{H}^{1}_{0,\sigma}(\Omega)\to(\bm{H}^{1}_{0,\sigma}(\Omega))^{\prime}$ by

[TABLE]

Using (2.10), Poincaré’s inequality, Hölder’s inequality, the boundedness of $\alpha_{\varepsilon}$ and properties of the trilinear form $b(\bm{u},\bm{v},\bm{w}):=\int_{\Omega}(\bm{u}\cdot\nabla)\bm{v}\cdot\bm{w}\,\mathrm{dx}\,$ (see [14, Lem. 4.1]), it can be shown similar to [14, Proof of Lem. 4.9] that $a(\cdot,\cdot)$ is a bounded and coercive bilinear form. Furthermore, defining

[TABLE]

and applying ($\mathrm{B1}$), ($\mathrm{B2}$), the fact that $\bm{G},\bm{u}_{\varepsilon}\in\bm{H}^{1}(\Omega)$ and Sobolev embeddings leads to the deduction that $\bm{F}(\bm{v})$ is a bounded linear form on $\bm{H}^{1}_{0,\sigma}(\Omega)$. Thus, by the Lax–Milgram theorem, we obtain a unique $\hat{\bm{q}}\in\bm{H}^{1}_{0,\sigma}(\Omega)$ such that

[TABLE]

This implies that the solution $\bm{q}_{\varepsilon}:=\hat{\bm{q}}+\bm{G}\in\bm{H}^{1}_{0}(\Omega)$ satisfies the weak formulation (4.11) with

[TABLE]

The existence of a unique adjoint pressure $\pi_{\varepsilon}\in L^{2}(\Omega)$ follows from standard results, see for instance [33, Lem. II.2.2.1]. Thus $(\bm{q}_{\varepsilon},\pi_{\varepsilon})$ is the unique weak solution to the adjoint system (4.9). ∎

4.4 Necessary optimality conditions

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

Theorem 4.2**.**

Let $\varphi_{\varepsilon}\in\mathbb{K}_{ad}$ be a minimizer of (2.5)-(2.7) with corresponding (unique) state variables $(\bm{u}_{\varepsilon},p_{\varepsilon})=\bm{S}_{\varepsilon}(\varphi_{\varepsilon})$, $\bm{u}_{\varepsilon}\in\bm{H}^{1}_{\bm{\sigma},\sigma}(\Omega)$, $p_{\varepsilon}\in L^{2}_{0}(\Omega)$. Furthermore, let $\{\lambda_{i}\}_{i=1}^{m_{1}+m_{2}}$ be the Lagrange multipliers associated to the integral state constraints $\{\mathcal{G}_{i}(\varphi_{\varepsilon})\}_{i=1}^{m_{1}+m_{2}}$, and $(q_{\varepsilon},\pi_{\varepsilon})$ be the unique solution to the adjoint system (4.9). Then, under ($\mathrm{B1}$), ($\mathrm{B2}$) and ($\mathrm{C1}$), the following optimality system is fulfilled:

[TABLE]

where $\mathrm{D}_{4}h$ is evaluated at $(x,\nabla\bm{u}_{\varepsilon},p_{\varepsilon},\nabla\varphi_{\varepsilon})$, $\bm{L}_{i}$ is evaluated at $(x,\nabla\bm{u}_{\varepsilon},p_{\varepsilon})$, and $\mathrm{D}_{5}b=B(x,\bm{u}_{\varepsilon},\nabla\bm{u}_{\varepsilon},p_{\varepsilon})\;\mathrm{D}_{2}z(x,\varphi_{\varepsilon})$, $\mathrm{D}_{5}K_{i}=\mathcal{K}_{i}(x,\bm{u}_{\varepsilon},\nabla\bm{u}_{\varepsilon},p_{\varepsilon})\;\mathrm{D}_{2}y_{i}(x,\varphi_{\varepsilon})$.

Proof.

In Sec. 4.1 we have shown that the reduced functional

[TABLE]

is Fréchet differentiable with respect to $\varphi_{\varepsilon}$, and in Sec. 4.2 we derived the gradient equation (4.7). We now want to rewrite (4.7) into a more convenient form using the adjoint system. For $\zeta\in\Phi$, let $(\bm{w}_{\varepsilon},r_{\varepsilon})$ denote the unique solution to the linearized state equations (2.12) corresponding to $\delta=\zeta-\varphi_{\varepsilon}$. Then, computing the derivative of $j_{\varepsilon}$ at $\varphi_{\varepsilon}$ in the direction $\delta$ leads to

[TABLE]

where in the above and for the rest of the proof $\{\mathrm{D}_{i}b\}_{i=2}^{5}$ are evaluated at $(x,\bm{u}_{\varepsilon},\nabla\bm{u}_{\varepsilon},p_{\varepsilon},\varphi_{\varepsilon})$ and $\{\mathrm{D}_{i}h\}_{i=2}^{4}$ are evaluated at $(x,\nabla\bm{u}_{\varepsilon},p_{\varepsilon},\nabla\varphi_{\varepsilon})$. Using the adjoint state $\bm{q}_{\varepsilon}$ as a test function in (2.12) (with $\delta=\zeta-\varphi_{\varepsilon}$) leads to

[TABLE]

where $g=\,\mathrm{div}\,q_{\varepsilon}$ as in (4.12). Using the linearized state $\bm{w}_{\varepsilon}$ as a test function in the adjoint system (4.11) leads to

[TABLE]

where we have used $\bm{w}_{\varepsilon}\in\bm{H}^{1}_{0,\sigma}(\Omega)$ to deduce that $\int_{\Omega}(\nabla q_{\varepsilon})^{\top}\cdot\nabla\bm{w}_{\varepsilon}\,\mathrm{dx}\,=0$ (see for instance [14, (4.29)]). Upon comparing terms in (4.17) and (4.18) we find that

[TABLE]

Using that $r_{\varepsilon}\in L^{2}_{0}(\Omega)$, $\,\mathrm{div}\,\bm{u}_{\varepsilon}=0$ in $\Omega$, $\bm{q}_{\varepsilon}=\bm{w}_{\varepsilon}=\bm{0}$ on $\partial\Omega$, and thus

[TABLE]

we can simplify (4.19) into

[TABLE]

and upon rearranging we obtain

[TABLE]

Substituting (4.20) into (4.16), we obtain

[TABLE]

Together with the gradient equation (4.7) and the distributional derivatives (4.5), we then obtain (4.15). ∎

Remark 4.1**.**

In the case where there is only a volume constraint, i.e., $m_{1}+m_{2}=1$ with $\mathcal{G}(\varphi):=\int_{\Omega}\varphi-\beta\,\mathrm{dx}\,$ for a fixed constant $\beta\in(-1,1)$, the existence of Lagrange multipliers using the Zowe–Kurcyusz constraint qualification has been shown in [18, Proof of Thm. 7.1] (for the case of inequality constraint), see also [12, Proof of Thm. 3] for another argument using geometric variations. For the case of equality constraint, we refer to [14, Proof of Thm. 4.10] which is based on a different argument.

5 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 $\Omega\subset\mathbb{R}^{2}$ with Lipschitz boundary, we study the following optimal control problem

[TABLE]

subject to $(\varphi,\bm{u},p)$ solving the porous-medium Navier–Stokes equations (2.6) and the following integral constraints:

[TABLE]

where $\bm{a}$ is a constant unit vector parallel to the flow direction $\bm{u}_{\infty}$, $M>0$ is a given positive constant representing an upper bound on the mass of the object, $\rho(x)\in L^{\infty}(\Omega)$ is a non-negative mass density, and $\beta\in(-1,1)$ so that the object is constraint to occupy a maximal volume of $\frac{1-\beta}{2}\left|\,\Omega\,\right|$.

The constraints $\mathcal{G}_{1}(\varphi)=0$ and $\mathcal{G}_{2}(\varphi)=0$ imply that the centre of mass for the object is located at the origin in $\mathbb{R}^{2}$ (which we can assume to hold without loss of generality by translating the domain $\Omega$). 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 $\hat{\alpha}_{\varepsilon}\left|\,\bm{u}\,\right|^{2}$ in the objective functional.

It is straightforward to check that the function

[TABLE]

fulfils ($\mathrm{A4}$) by the application of the Young’s inequality. Furthermore, it is shown in [14, Proof of Thm. 4.1 and Rmk. 4.2] that the functional

[TABLE]

is bounded from below for $\varphi\in H^{1}(\Omega)\cap L^{\infty}(\Omega)$ with $s_{a}\leq\varphi\leq s_{b}$ a.e. in $\Omega$, $\bm{u}\in\bm{H}^{1}_{\bm{g},\sigma}(\Omega)$ and $p\in L^{2}_{0}(\Omega)$, and satisfies due to the product of weak-strong convergence:

[TABLE]

for sequences $\varphi_{n}\rightharpoonup\varphi$ in $H^{1}(\Omega)$, $\bm{u}_{n}\to\bm{u}$ in $\bm{H}^{1}(\Omega)$ and $p_{n}\to p$ in $L^{2}(\Omega)$. Hence, ($\mathrm{A8}$) is also fulfilled. A short computation shows that

[TABLE]

and as $\bm{a}$ is a constant vector, one can infer that ($\mathrm{B1}$) (specifically (4.1)) is also fulfilled. Then, it remains to verify ($\mathrm{A5}$), ($\mathrm{A6}$), ($\mathrm{B2}$) and ($\mathrm{C1}$) for the existence of Lagrange multipliers for the integral constraints $\mathcal{G}_{1},\dots,\mathcal{G}_{4}$, and show that the admissible set $\mathbb{K}_{ad}$ is non-empty.

For the latter, note that we have the trivial example $\phi\equiv 1\in\mathbb{K}_{ad}$ which corresponds to the case where there is no object in the domain $\Omega$. In the following we will construct a non-trivial example in order to rule out the possibility where $\mathbb{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 $\phi\in H^{1}(\Omega)$, $-1\leq\phi\leq 1$ a.e. in $\Omega$ such that

[TABLE]

which is equivalent to choosing an object $\{\phi=-1\}$ with its centre of mass at the origin with volume bounded above by $\tfrac{1-\beta}{2}\left|\,\Omega\,\right|$. Note that the mapping $\phi\mapsto\int_{\Omega}\frac{1}{2}\rho(x)(1-\phi)\,\mathrm{dx}\,$ is continuous, and thus we can always decrease the volume of the object region $\{\phi=-1\}$ to ensure the mass is bounded from above by the constant $M$. This ensures that $\phi\in\mathbb{K}_{ad}$ and hence ($\mathrm{A7}$) is satisfied.

As $\Omega$ is a bounded domain, the functions $x_{1}$, $x_{2}$ are bounded. Then, upon setting

[TABLE]

we observe that ($\mathrm{A5}$), ($\mathrm{A6}$) and ($\mathrm{B2}$) are fulfilled by the above choices. Then, by Thm. 3.2 we are guaranteed the existence of a minimizer $\varphi_{\varepsilon}$ to the optimal control problem. To verify the main assumption ($\mathrm{C1}$) and derive the optimality conditions, we have to show that for an arbitrary $z=(z_{1},z_{2},z_{3},z_{4})^{\top}\in\mathbb{Y}=\mathbb{R}^{4}$, there exists one function $\psi_{*}\in\Phi$, along with non-negative constants $\tau_{1},\dots,\tau_{4},\xi_{1},\xi_{2},\eta_{1},\eta_{2}$ such that the following four conditions are fulfilled simultaneously:

[TABLE]

Due to their nature as equality constraints, we can use the fact that $\mathcal{G}_{1}(\varphi_{\varepsilon})=\mathcal{G}_{2}(\varphi_{\varepsilon})=0$ to simplify (5.1a) into

[TABLE]

We first argue for (5.2). As the origin $\bm{0}\notin\partial\Omega$, this implies that $\Omega$ has non-empty intersections with the four quadrants of $\mathbb{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 $\tau_{1}$ (resp. $\tau_{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 $\psi_{*}\in\Phi$ not identically equal to $1$ with $\beta<\psi_{*}\leq 1$ a.e. in $\Omega$ such that the non-empty set $A:=\mathrm{supp}(1-\psi_{*})$ has Lebesgue measure

[TABLE]

and satisfies

[TABLE]

Then, we set

[TABLE]

and thanks to the fact that $\psi_{*}\leq 1$ a.e. in $\Omega$, the function $1-\psi_{*}$ is non-negative in $\Omega$ and only positive in $A$. The location of $A$ implies that the integrand $(1-\psi_{*})x_{i}$ has the same sign as $z_{i}$ for $i=1,2$, and so $\tau_{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 $\psi_{*}$ 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 $\mathcal{G}_{3}(\varphi_{\varepsilon})$ is not active for the minimizer $\varphi_{\varepsilon}$, i.e., $\varphi_{\varepsilon}$ satisfies $\int_{\Omega}\frac{1}{2}\rho(x)(1-\varphi_{\varepsilon})\,\mathrm{dx}\,<M$. Then, we can choose $\tau_{3}=0$ and it holds that

[TABLE]

Hence, we have fulfilled (5.1b) without making use of the function $\psi_{*}$. On the other hand, if $\mathcal{G}_{3}(\varphi_{\varepsilon})$ is active, i.e., $\int_{\Omega}\frac{1}{2}\rho(x)(1-\varphi_{\varepsilon})\,\mathrm{dx}\,=M$, the condition (5.1b) simplifies to

[TABLE]

A short calculation using (5.3) shows that the quantity in the bracket is positive, and so

[TABLE]

which implies that (5.1b) is fulfilled. Indeed, we see that

[TABLE]

For the volume constraint (5.1c) we again divide the argument into two cases: if $\mathcal{G}_{4}(\varphi_{\varepsilon})$ is inactive, then (5.1c) holds automatically without the use of the function $\psi_{*}$, and if $\mathcal{G}_{4}(\varphi_{\varepsilon})$ is active, then using $\psi_{*}>\beta$ yields the desired result.

As a consequence, ($\mathrm{C1}$) is fulfilled and we obtain the existence of Lagrange multipliers $\lambda_{1},\lambda_{2}\in\mathbb{R}$, $\lambda_{3},\lambda_{4}\in\mathbb{R}_{\geq 0}$. By Thm. 4.2 the first order optimality condition is

[TABLE]

together with the complementary slackness conditions

[TABLE]

Remark 5.1**.**

We point out that the mass constraint $\mathcal{G}_{3}(\varphi)=M-\int_{\Omega}\tfrac{1}{2}\rho(x)(1-\varphi)\,\mathrm{dx}\,\geq 0$ can also be thought of as a constraint on a construction cost, where the value $\rho(x)>0$ represents the cost of building the object at the point $x\in\Omega$, 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 $\Omega\subset\mathbb{R}^{2}$ with Lipschitz boundary, we study the following optimal control problem

[TABLE]

subject to $(\varphi,\bm{u},p)$ solving the porous-medium Navier–Stokes equations (2.6) (with zero body force $\bm{f}=\bm{0}$) and the following integral constraint:

[TABLE]

where $\bm{a}=-\bm{u}_{\infty}^{\perp}$ is the negative unit vector perpendicular to the flow direction $\bm{u}_{\infty}$ and $D>0$ is a given positive constant representing an upper bound on the total potential power. Then, upon setting

[TABLE]

we see that ($\mathrm{A5}$), ($\mathrm{A6}$) and ($\mathrm{B2}$) are fulfilled. In this setting we observe that the trivial example $\varphi\equiv-1$ belongs to the admissible set of design functions $\mathbb{K}_{ad}$. A non-trivial example can be found if the domain $\Omega$ is sufficiently large or the viscosity $\mu$ is sufficiently small. Indeed, let $\varphi$ be a function in $H^{1}(\Omega)$ with $-1\leq\varphi\leq 1$ a.e. in $\Omega$ but not identically equal to $1$ or $-1$. Denote by $\bm{u}$ the unique velocity field associated to the state equation (2.2), then by (2.10) it holds that

[TABLE]

where from (2.11) the constant $K_{\Omega}$ is $K_{\Omega}=\frac{1}{2}\left|\,\Omega\,\right|^{\frac{1}{2}}$ in two dimensions. Note that the above upper bound is independent of $\varphi$. For any fixed positive constant $D$, we can take a sufficiently large domain $\Omega$ or sufficiently small viscosity $\mu$, so that $\frac{\mu^{3}}{2K_{\Omega}^{2}}\leq D$. Then, this implies that $\mathbb{K}_{ad}$ is non-empty and thus ($\mathrm{A7}$) is fulfilled. Furthermore, following the arguments above, we deduce by Thm. 3.2 that there exists at least one minimizer $\varphi_{\varepsilon}$ to the optimal control problem. Writing $\mathcal{G}(\varphi_{\varepsilon})=G(\varphi_{\varepsilon},\bm{S}_{\varepsilon}(\varphi_{\varepsilon}))$, to verify the assumption ($\mathrm{C1}$), we have to show for an arbitrary $z\in\mathbb{R}$, there exists one function $\psi_{*}\in\Phi$ along with non-negative constants $\tau,\xi,\eta$ such that

[TABLE]

Observe that for this particular setting (with a large domain $\Omega$ or small viscosity $\mu$), the inequality $D\geq\frac{\mu^{3}}{2K_{\Omega^{2}}}$ holds independently of the minimizer $\varphi_{\varepsilon}$, and thus by (5.4) $\mathcal{G}(\varphi_{\varepsilon})>0$ holds. In particular, the constraint is always inactive, and we do not need to find the function $\psi_{*}$ as

[TABLE]

Furthermore, as the constraint is always inactive the complementary slackness condition (4.8) implies that the associated Lagrange multiplier $\lambda$ is zero. Hence, by Thm. 4.2 the first order optimality condition is

[TABLE]

6 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 $\varphi\in H^{1}(\Omega)\cap L^{\infty}(\Omega)$, 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}(\Omega)\cap L^{\infty}(\Omega)$ 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 $\Psi$ we use the double-obstacle free energy, namely

[TABLE]

From this we obtain the constraint $\left|\,\varphi_{\varepsilon}\,\right|\leq 1$, and $c_{0}=\frac{\pi}{2}$, where $c_{0}$ is the constant defined in (2.1). Although the double-obstacle potential (6.1) does not satisfy ($\mathrm{A1}$), the analysis is not affected once we choose $s_{a}=-1$ and $s_{b}=1$, so that $\left|\,\varphi_{\varepsilon}\,\right|\leq 1$ and the potential becomes $\Psi(\varphi_{\varepsilon})=\frac{1}{2}(1-\varphi_{\varepsilon}^{2})$. We refer the reader to [11, 13] which also uses the double-obstacle potential (6.1). For the porous-medium term $\alpha_{\varepsilon}(\varphi_{\varepsilon})$ in the state equations (2.2a) we choose

[TABLE]

with a fixed positive constant $\overline{\alpha}$, and ($\mathrm{A0}$) is fulfilled with $s_{a}=-1$ and $s_{b}=1$. We choose

[TABLE]

so that ($\mathrm{A2}$) is also satisfied. For the remaining part of this section, we denote both variables by $\alpha_{\varepsilon}$, set $\bm{f}=\bm{0}$ in (2.2), and define

[TABLE]

6.1 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 $\varphi_{\varepsilon}$, $p_{\varepsilon}$ and $\pi_{\varepsilon}$ and piecewise quadratic and globally continuous finite elements for $\bm{u}_{\varepsilon}$ and $\bm{q}_{\varepsilon}$ on a conforming triangulation of the domain $\Omega$.

It is well-known that in phase field applications the variable $\varphi_{\varepsilon}$ 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 $\left|\,\varphi_{\varepsilon}\,\right|=1$. Thus in every step of the PDAS we solve the problems on the inactive set $\left|\,\varphi_{\varepsilon}\,\right|<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 $\varphi_{\varepsilon}$ only, namely we use the jumps of the normal derivatives of $\varphi_{\varepsilon}$ 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.

6.2 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 $\Omega=(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

[TABLE]

For the objective functional we define a ‘rock’ centered at $\bm{m}$ with radius $\sigma$ and associated cost $c$ as

[TABLE]

We consider the functions

[TABLE]

where

[TABLE]

The optimization problem (2.5) then becomes

[TABLE]

subject to $\varphi\in\Phi$, $\bm{u}\in\bm{H}^{1}_{\bm{g},\sigma}(\Omega)$, $p\in L^{2}_{0}(\Omega)$ satisfying (2.6), and $\alpha_{\varepsilon}$ 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., $\varphi_{\varepsilon}^{0}\equiv 0$, on a homogeneous coarse grid with mesh size $h=1/20$ yielding 685 degrees of unknowns for $\varphi_{\varepsilon}$. 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: $\sigma=0.15$, $c=50$, $\varepsilon=0.01$, $\overline{\alpha}=5$, $\mu=0.02$ and $\gamma=0.001$. To stress the advantages of our approach in Fig. 2 we show $\varphi_{\varepsilon}$ at various stages of the optimization procedure.

6.3 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 $\varphi_{\varepsilon}$. However, here we employ the VMPT method to solve the optimization problem, which reads

[TABLE]

subject to \varphi\in\color[rgb]{0,0,0}\Phi\color[rgb]{0,0,0}, $\bm{u}\in\bm{H}^{1}_{\bm{g},\sigma}(\Omega)$, $p\in L^{2}_{0}(\Omega)$ satisfying (2.6) and the volume constraint (see (1.9))

[TABLE]

We use the parameters from [14], namely $\Omega=(0,1.7)\times(0,0.4)$, $\varepsilon=0.00025$, $\overline{\alpha}=0.03$, $\mu=0.001$ and $\gamma=0.01$. The boundary velocity is set to $\bm{g}=(1,0)^{\top}$ to stay close to the analysis and we initialize the optimization with $\varphi_{\varepsilon}^{0}(x):=-R[(0.5,0.2)^{\top},0.25,-1](x)$, i.e., a ball around $m=(0.5,0.2)^{\top}$ with radius $r=0.25$. For the volume constraint, we choose $\beta_{2}=\beta=0.975$, i.e., $\int_{\Omega}\varphi\,\mathrm{dx}\,\leq 0.663$.

To be able to use only a small number of unknown as long as possible, we start the optimization with $\varepsilon=0.008$ and a maximum number of allowed degrees of freedom of 10000. We halve the value of $\varepsilon$ 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 $\varepsilon$, namely $\varepsilon\in\{0.008,0.004,0.002,0.001,0.0005,0.00025\}$. In Table 1 we show the diffuse interface drag

[TABLE]

and the sharp interface drag

[TABLE]

by evaluation with $\bm{a}=(1.0,0.0)^{\top}$ over $\Gamma=\{\varphi_{\varepsilon}=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.

6.4 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 $\Gamma$ 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 \bm{a}\color[rgb]{0,0,0}\neq\bm{0}\color[rgb]{0,0,0}, let $\bm{\eta}$ be a smooth vector field such that

[TABLE]

This can be done since it is assumed in the introduction that $\Gamma$ does not intersect with $\partial\Omega$. Then, by taking the scalar product of (1.4a) with $\bm{\eta}$, we obtain

[TABLE]

Integrating by parts and noting that the boundary integrals over $\partial\Omega$ vanish owning to $\bm{\eta}=\bm{0}$ on $\partial\Omega$ yields

[TABLE]

Here we have also used that $\bm{u}$ has no tangential component on $\Gamma$ due to the no-slip condition $\bm{u}=\bm{0}$ on $\Gamma$, and together with the divergence-free condition, we obtain that $(\nabla\bm{u})^{\top}\bm{\nu}=\bm{0}$ on $\Gamma$ (see [14, §2] for more details). This implies that we can also consider the following function as the volume formulation of the drag (if $\bm{a}$ is parallel to the flow direction)

[TABLE]

Alternatively, using integration by parts and the boundary conditions $\bm{\eta}=\bm{0}$ on $\partial\Omega$ and $\bm{u}=\bm{0}$ on $\Gamma$, we see that

[TABLE]

and so we may also use the function

[TABLE]

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

The aim of this section is to compare results of Sec. 6.3 with the following drag minimization problem

[TABLE]

subject to $\varphi\in\Phi$, $\bm{u}\in\bm{H}^{1}_{\bm{g},\sigma}(\Omega)$, $p\in L^{2}_{0}(\Omega)$ satisfying (2.6) and the volume constraint (6.3). In particular, we compare the optimal shapes obtained with volume formulation (6.9) and 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 $\bm{\eta}\equiv\bm{a}$ on $(0.15,1.0)\times(0.13,0.27)$.

Using the surface formulation (6.4) we observe that for larger values of $\overline{\alpha}$ (the constant in (6.2)) an interfacial region $\{\left|\,\varphi_{\varepsilon}\,\right|<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 $\varepsilon$. We stress that, in Sec. 6.3 for $\varepsilon=0.00025$ the the velocity $\left|\,\bm{u}_{\varepsilon}\,\right|$ 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 $\bm{a}$, namely the vector field $\bm{\eta}$ which has to vanish at $\partial\Omega$. We define $\bm{\eta}$ as the solution of a Poisson problem on $\Omega$ with $\bm{\eta}=\bm{a}$ on a square around the object and $\bm{\eta}=0$ on $\partial\Omega$. That is, let $S$ denote a square such that $\{\varphi_{\varepsilon}=-1\}\subset S$ and $\partial S\cap\partial\Omega=\emptyset$, then we solve

[TABLE]

For small values of $\overline{\alpha}$, 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 $\bm{\eta}=\bm{0}$ on $\partial\Omega$, the magnitude $\left|\,\bm{\eta}\,\right|$ is small close to the inflow and outflow boundaries, which results in small drag forces. On the other hand, for $\overline{\alpha}$ small, the porous-medium penalization term $\int_{\Omega}\alpha_{\varepsilon}(\varphi)\left|\,\bm{u}\,\right|^{2}\,\mathrm{dx}\,$ is small, and thus the value of the objective functional can be reduced by placing material in regions where $\left|\,\bm{\eta}\,\right|$ is small. Therefore, in contrast to the surface formulation (6.4), large values of $\overline{\alpha}$ 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 $\varepsilon$.

We use the setup from Sec. 6.3 with only one modification, that we set $\mu=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 $\overline{\alpha}$ to model objects with smaller permeability. To show the influence of $\overline{\alpha}$ in Fig. 5 we show the optimal shape for the above parameters, but using a larger value $\overline{\alpha}=1$ and $\mu=0.01$ (left) and $\mu=0.001$ (right). For $\mu=0.01$ we observe, that we get a sharper rear of the 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 $\overline{\alpha}=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 $\overline{\alpha}=1$, and $\mu=0.001$ we observe a symmetric airfoil shape.

6.5 Maximizing the lift with constraints on the total potential power

We give an example of dealing with a state constraint, namely we consider the maximization of the lift of an object under the constraint that the total potential power is bounded by some given value. This is a non-linear constraint on the state variables of the constraint optimization problem, namely the velocity field. To treat the highly non-linear potential power constraint we use Moreau–Yosida relaxation.

The optimization problem we solve is

[TABLE]

subject to $\varphi\in\Phi$, $\bm{u}\in\bm{H}^{1}_{\bm{g},\sigma}(\Omega)$, $p\in L^{2}_{0}(\Omega)$ satisfying (2.6), where to realize the maximization of lift, we set $\bm{a}=(-1,0)^{\top}$ as the negative unit vector perpendicular to the flow direction, and the parameter $s>0$ penalizes violation of the constraint that the total potential power of the fluid region must be less than or equal to a prescribed value $D$:

[TABLE]

The integral constraint for this optimization problem are volume constraints of the form $0.663\leq\int_{\Omega}\varphi_{\varepsilon}\,\mathrm{dx}\,\leq 0.665$ and the center of mass is fixed at $(0.5,0.2)^{\top}$.

The set up is similar to that in Sec. 6.3, where we set $\Omega=(0.0,1.7)\times(0.0,0.4)$, $\bm{g}=(1,0)^{\top}$, $\varphi_{\varepsilon}^{0}(x):=-R[(0.5,0.2)^{\top},0.25,-1](x)$, i.e., a circle around $\bm{m}=(0.5,0.2)^{\top}$ with radius $r=0.25$. For the penalization parameter we choose $s=100$ and further numerical parameters are $\varepsilon=0.02$, $\overline{\alpha}=2$, $\mu=0.01$, $\gamma=0.001$, and $D=0.06$.

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].

7 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 $\varphi$, while the second involves the state variable $\bm{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.