Numerical analysis of a topology optimization problem for Stokes flow

T. Borrvall and J. Petersson [Topology optimization of fluids in Stokes flow, International Journal for Numerical Methods in Fluids 41 (1) (2003) 77--107] developed the first model for topology optimization of fluids in Stokes flow. They proved the existence of minimizers in the infinite-dimensional setting and showed that a suitably chosen finite element method will converge in a weak(-*) sense to an unspecified solution. In this work, we prove novel regularity results and extend their numerical analysis. In particular, given an isolated local minimizer to the infinite-dimensional problem, we show that there exists a sequence of finite element solutions, satisfying necessary first-order optimality conditions, that strongly converges to it. We also provide the first numerical investigation into convergence rates.


Introduction
Topology optimization has become an effective technique in structural and additive manufacturing and has found multiple uses in medicine, architecture, and engineering [1,2,3]. The objective is to find the optimal distribution of a fluid or solid within a given domain that minimizes a problem-specific cost functional [4,5]. In this paper we consider a model for topology optimization for fluids proposed by Borrvall and Petersson [6]. Their seminal work has become the foundation for a number of developments in recent years [7,8,9,10,11,12,13,14,15,16,17,18]. Their goal was to minimize the power dissipation of a fluid that satisfies both the Stokes equations and a volume constraint restricting the proportion of the domain that the fluid can occupy. In their paper, they derived generalized Stokes equations, which incorporate the classical velocity and pressure terms but also introduce a variable, ρ, that represents the material distribution of the fluid over the given domain. The presence of fluid is indicated by a value of one in the material distribution whereas absence of fluid is represented by a value of zero. It would be ideal for ρ : Ω → {0, 1}, in order to remove any ambiguity in the solutions, however, in general, this is a numerically intractable objective. In the Borrvall-Petersson model ρ : Ω → [0, 1], but the model is regularized with an inverse permeability term, α, which favors solutions where ρ is close to zero or one. From the generalized Stokes equations, Borrvall and Petersson formulated an infinite-dimensional nonconvex optimization problem with inequality, PDE and box constraints. There have been studies on the numerics of a Ginzburg-Landau regularization of the objective functional that can be shown to recover solutions with sharp transitions in the material distribution [19,20]. Notably, Garcke et al. [21] derived a posteriori error estimators designed to resolve the interfaces in the material distribution for the Navier-Stokes extension to the Borrvall-Petersson problem. As far as we are aware, there exist only a couple of results dealing with weak(-*) convergence of discretized solutions, as the mesh size tends to zero, to solutions of the Borrvall-Petersson problem on the whole domain [6,22]. Moreover, there have been no results concerning strong convergence nor the convergence to all the isolated minimizers of the problem.
In the original paper [6], it is shown that a minimizing velocity and material distribution to the optimization problem exist [6,Th. 3.1]; however, the minimizer is not necessarily unique [6,Sec. 4.5]. It is also shown that there exist finite element solutions that converge to a minimizer of the problem [6,Th. 3.2]. The proven convergence is weak in the approximation of the velocity and weak-* in the material distribution, with no results for the pressure. In addition, Borrvall and Petersson show that the approximation of material distribution strongly converges to a solution in L s (Ω b ), s ∈ [1, ∞), where Ω b is any measurable subset of Ω in which the material distribution that solves the infinite-dimensional problem is equal to zero or one a.e. [6,Sec. 3.3]. Weak-* convergence permits large oscillations in the material distribution, called checkerboarding, which could occur in areas where the material distribution is not zero or one under the current results. However, in practice, checkerboarding is not observed in these regions. Since there can be multiple solutions, the nature of the convergence is ambiguous. In particular, it is not clear if there are sequences of finite element solutions converging to every solution of the infinite-dimensional problem.
Our goal is to extend and refine the analysis of Borrvall and Petersson. We show that, given an isolated minimizer to the infinite-dimensional problem, there exists a sequence of finite element solutions, satisfying the necessary first-order optimality conditions, that strongly converges. In particular, we strengthen the convergence from weak convergence in H 1 (Ω) d to strong convergence in H 1 (Ω) d for the velocity, and from weak-* convergence in L ∞ (Ω) to strong convergence in L s (Ω), s ∈ [1, ∞) for the material distribution. Moreover, in the case of a ho-mogeneous Dirichlet boundary condition, we show that the material distribution is weakly differentiable inside any compact subset of the support of the velocity; more specifically ρ ∈ H 1 (U θ ), for any θ > 0, where U θ is any measurable subset of Ω in which |u| 2 ≥ θ > 0 a.e. in U θ . This analysis confirms that checkerboarding cannot occur under mild assumptions on the model. Hence, isolated minimizers of the problem can be well approximated using the finite element method. We conclude with a numerical investigation into the convergence of the finite element solutions.
By first considering the optimization problem, we derive necessary first-order optimality conditions in Section 2. By construction, the generalized Stokes equations are satisfied but we also show that the material distribution satisfies a variational inequality. In Section 3, we show that under moderate assumptions, the material distribution is weakly differentiable in the case of a homogeneous Dirichlet boundary condition for the velocity. We tackle the issue of multiple local minima in Section 4, by considering closed balls around isolated local minimizers. In that section, we also prove that for each isolated minimizer there exists a sequence of finite element solutions to the discretized first-order optimality conditions, which strongly converges to the solution of the infinite-dimensional problem. In Section 5, we computationally investigate the convergence of sequences of finite element approximations to the respective solutions of the infinite-dimensional problem.

Existence and necessary first-order optimality conditions
The topology optimization problem of Borrvall and Petersson [6] is as follows: find the velocity, u, and the material distribution, ρ, that solve the minimization problem where, in Ω}, is the (constant) viscosity, and γ ∈ (0, 1) is the volume fraction. Here, α is the inverse permeability, modeling the influence of the material distribution on the flow. For values of ρ close to one, α(ρ) is small permitting fluid flow; for small values of ρ, α(ρ) is very large, restricting fluid flow. The function α is assumed to have the following properties: (A2) α is convex and monotonically decreasing; (A3) α(0) = α and α(1) = α, generating an operator also denoted α : C γ → L ∞ (Ω; [α, α]). Typically, in the literature α takes the form [6] α where q > 0 is a penalty parameter, so that lim q→∞ α(ρ) =ᾱ(1 − ρ). Furthermore, | ∂Ω is to be understood in the boundary trace sense [23,Ch. 5.5], The next theorem is due to Borrvall  Then, there exists a pair (u, ρ) ∈ H 1 g,div (Ω) d × C γ that locally minimizes J, as defined in (BP).

Remark 2.
Although a solution exists, it is not necessarily unique, as observed in the numerical examples in Section 5, since the optimization problem is nonconvex. The nonconvexity is caused by the term α(ρ)|u| 2 in (BP). A rough argument examines the second partial Fréchet derivative of J with respect to u and ρ (assuming that it exists), i.e., for suitable variations v, η, in two dimensions and 3/2 ≤ s ≤ ∞ in three dimensions. Moreover, for all v ∈ H 1 0 (Ω) d and η ∈ C γ we have that where J u (u, ρ) denotes the Fréchet derivative of J with respect to u and J ρ (u, ρ) denotes the Fréchet derivative of J with respect to ρ. Remark 3. It can be checked that if α is (n+1)-times continuously differentiable then J is n-times Fréchet differentiable with respect to u and ρ.
The following proposition is the main result of this section. We show that if (u, ρ) is a minimizer of the optimization problem (BP), then the minimizer also satisfies first-order optimality conditions consisting of two equations and a variational inequality.
, 3}, and α satisfies properties (A1)-(A4). Consider a local or global minimizer (u, ρ) ∈ H 1 g,div (Ω) d × C γ of (BP). Then, there exists a unique Lagrange multiplier p ∈ L 2 0 (Ω) such that the following necessary first-order optimality conditions hold: where L 2 0 (Ω) := q ∈ L 2 (Ω) : Proof. We will first show that (FOC1)-(FOC2) are satisfied by generalizing arguments, used for the Stokes system with a homogeneous Dirichlet boundary condition, found in [25]. For ease of notation we define X g := H 1 g (Ω) d , X 0 := The respective dual spaces of X 0 , V 0 and M are denoted with * . We also define the associated operators, A ∈ L(X g , X * 0 ), B ∈ L(X g , M ) and B 0 ∈ L(X 0 , M ) by We note that ker(B 0 ) = V 0 . From Theorem 1, we know that there exists a pair (u, ρ) ∈ V g × C γ that is a local minimizer for (BP). For any given v ∈ V 0 , we see that u + tv ∈ V g , t ∈ R. If (u, ρ) ∈ V g × C γ is a minimizer, then, by definition, there exists an r > 0 such that, for any (w, η) ∈ V g × C γ , (w, η) = (u, ρ) that satisfies we have that J(u, ρ) ≤ J(w, η). Hence, for any given v ∈ V 0 , if 0 < t < r/ v H 1 (Ω) , the following inequality holds By Proposition 1, J is Fréchet differentiable, and therefore also Gateaux differentiable, with respect to u. Hence as t → 0 + , we see that By considering the same reasoning with t < 0, we deduce that From Proposition 1, we know that J u (u, ρ) = Au − f and hence there exists a β > 0 such that, for all q ∈ M, where B * 0 is the dual operator of B 0 , defined by v, B * 0 q = B 0 v, q . This implies that B * 0 is injective (and therefore bijective) from M into ImB * 0 . Furthermore, it also implies that (B * 0 ) −1 is continuous. Consider f ∈ ImB * 0 ; then, there exists a q ∈ M such that f = B * 0 q and Therefore, ImB * 0 is closed.
Since ImB * 0 is closed, by Banach's closed range theorem, we know that Since B * 0 is injective, p is also unique. Since u ∈ V g , we have that Bu = 0. Hence (FOC1) and (FOC2) hold.
We will now show that (FOC3) holds. We note that C γ is a convex subset of a linear space. For any given ζ, η ∈ C γ and t ∈ [0, 1], we therefore have that ζ + t(η − ζ) ∈ C γ . Since (u, ρ) is a local minimizer, it follows that for each η ∈ C γ , if 0 < t < r/ η − ρ L ∞ (Ω) , with r as in (2.6), then (2.14) From Proposition 1, we know that J is Fréchet differentiable, and therefore also Gateaux differentiable, with respect to ρ. Hence, by taking the limit as t → 0, we see that Therefore (FOC3) holds.
The following lemma will be used in the proof of the next proposition.
Lemma 2. Consider a nonzero function η ∈ C γ and the measurable non-empty set E ⊂⊂ supp(η), where ⊂⊂ denotes that the containment is compact and supp denotes the support of a function, i.e. η > 0 a.e. in E. Then, there exists an > 0 such that, for all ∈ (0, ], there exists a set E ⊆ E, |E | > 0 where η > a.e. in E . Proof. For a contradiction, suppose that there exists no such such that E exists. This implies that for all n ≥ 0, |E\Ê n | = 0, (2.16) whereÊ n := {0 ≤ η ≤ 1/n a.e. in E}. We see that ∅ = E\Ê 1 ⊆ E\Ê 2 ⊆ · · · ⊆ E\Ê n ⊆ · · · , i.e. E\Ê n is nondecreasing. Note that the limit of a nondecreasing sequence of sets (A n ) can be defined as lim n→∞ A n := ∪ n≥1 A n . By (2.16) we note that lim n→∞ |E\Ê n | = 0.  Now we see that where the first equality follows from (2.18), the second equality follows from the definition of the limit of a nondecreasing sequence of sets, the third equality follows from the continuity of the Lebesgue measure, and the fourth equality follows from (2.17). (2.19) is a contradiction and, therefore, such an > 0 must exist. By choosing E = E for all 0 < ≤ , we conclude that the statement holds for all ∈ (0, ].
In the result that follows, we are required to distinguish between different types of global and local minimizers.
Definition 1 (Strict minimizer). Let Z be a Banach space and suppose that z 0 ∈ Z is a local or global minimizer of the functional J : Z → R. We say that z 0 is a strict minimizer if there exists an open neighborhood E ⊂ Z of z 0 such that J(z 0 ) < J(z) for all z = z 0 , z ∈ E.
Definition 2 (Isolated minimizer). Let Z be a Banach space and suppose that z 0 ∈ Z is a local or global minimizer of the functional J : Z → R. We say that z 0 is isolated if there exists an open neighborhood E ⊂ Z of z 0 such that there are no other minimizers contained in E.

Remark 4.
If z is an isolated minimizer, then it is also a strict minimizer.
The following proposition is a property of strict minimizers that will be useful for the numerical analysis of the finite element method.

Regularity of ρ
In this section we show that ρ ∈ C γ possesses higher regularity in the case of a homogeneous Dirichlet boundary condition on u and if α satisfies a stronger (but not restrictive) convexity assumption.
Proof of Theorem 3. Let ∂ x k denote the partial derivative with respect to x k . If we can bound the L 2 -norm of the difference quotients of ρ, in all coordinate directions in U θ , above by constants independent of h, then, by taking the weak limit, we deduce that ∂ x k ρ exists as an element of L 2 (U θ ) for 1 ≤ k ≤ d.
The variational inequality on ρ states that We define U ⊂ Ω as U := supp(u) and fix an open, bounded and connected domainΩ such thatΩ = Ω if U ⊂⊂ Ω and Ω ⊂⊂Ω otherwise. In the case where U is not a compact subset of Ω, we extend u and ρ by zero to the whole of R d . Since the trace of u is zero on the boundary, the extension of u by zero Analogously, we define u h as We define the difference quotient, D h k , in the k-th coordinate direction, as which implies that 0 ≤ η ≤ 1 a.e. in Ω and If we multiply (3.1) through by 4 and divide by h 2 we see that We note that, Hence, because u is zero outside of Ω, (3.2) is equivalent to In order to obtain a first-order difference quotient, we will perform the finite difference analogue of integration by parts to shift the D −h k operator from D h k ρ to α (ρ)|u| 2 . We note that, by definition, the left-hand side of (3.3) is equal to which by a change of variables is equal to We note that U ⊂⊂Ω and |h| < (1/2)dist(U, ∂Ω), which implies that U ⊂⊂ Ω − he k . Therefore, Now we wish to rewrite D h k (α (ρ)|u| 2 ) in a form that we can decouple from D h k ρ in order to be able to bound (3.6) above and below. Now, Therefore, from (3.6) we see that Now, Hence from (3.7) we find that Subtracting the second term on the left-hand side in (3.8) from both sides, taking absolute values on the right-hand side, using the Cauchy-Schwarz inequality and multiplying by 2, we see that Furthermore we note that A ≥ α min and Hence, using Cauchy's inequality and Young's inequality, we see that (3.10) By fixing = α min , from (3.10) we see that, Now |α (ρ h ) + α (ρ)| 2 is bounded above by 4 sup ζ∈Cγ |α (ζ)| 2 which is independent of h. Consider a setΩ ⊂ R d such thatΩ ⊂⊂Ω. We note that u ∈ H 1 (Ω) d . By applying Theorem 3 in [23, pg. 294], we see that whereC,Ĉ and C are constants. The bound is independent of h and k. Because, by hypothesis, there exists a subset U θ ⊂ Ω such that, (3.14) From (3.14) we see that there exists a function η k ∈ L 2 (U θ ) and a subsequence h i → 0 such that, Finally, we wish to identify η k with ∂ x k ρ. First choose any smooth and compactly supported function, φ ∈ C ∞ c (U θ ). We note that Hence η k = ∂ x k ρ a.e. in U θ for k = 1, . . . , d. Therefore, from (3.13) we see by weak lower semicontinuity that for some constant C. We conclude that ρ ∈ H 1 (U θ ) ∩ C γ , θ > 0.

Remark 6.
The assumption that the boundary datum g = 0 on ∂Ω is required so as to expand the domain of integration from Ω in (3.2) toΩ in (3.3) in order to perform the finite difference analogue of integration by parts in (3.4).
A homogeneous Dirichlet boundary condition on u is rarely imposed in practice. However, we observe during numerical experiments that ρ possesses additional regularity in the case of inhomogeneous Dirichlet boundary conditions, and we hypothesize that the results can be generalized to that case.

Finite element approximation
We will be approximately solving (FOC1)-(FOC3) by approximating the solutions of the infinite-dimensional problem with finite element functions. Borrvall and Petersson [6, Sec. 3.3] considered a piecewise constant finite element approximation of the material distribution coupled with an inf-sup stable quadrilateral finite element approximation of the velocity and the pressure. They showed that such approximations of the velocity and material distribution converge to an unspecified solution (u, ρ) of (BP) in the following sense [6]: where Ω b is any measurable subset of Ω where ρ is equal to zero or one a.e. There are no proven convergence results for the finite element approximation of the pressure, p. Since there can be multiple local minimizers, it is unclear which solution the sequence of finite element solutions is converging to. In this section, we consider any suitable conforming mixed finite element space such that the velocity and pressure spaces are inf-sup stable. We prove that, for every isolated minimizer of (BP), there exist sequences of finite element solutions to the discretized first-order optimality conditions that strongly converge to the minimizer as the mesh size tends to zero. More specifically, we show that, for each isolated minimizer of the infinite-dimensional problem, there exist (possibly different) sequences of finite element We emphasize that the results hold in the case where the local minima are isolated. Consider the conforming finite element spaces . In general, it will not be possible to represent the boundary data g exactly in the velocity finite element space. Hence, for each h, we instead consider boundary data g h (which can be represented) and assume that (F1) g h → g strongly in H 1/2 (∂Ω) d .
We now define the space X g h ,h := {v h ∈ X h : v h | ∂Ω = g h }. We will also assume that: (F2) X 0,h and M h satisfy the following inf-sup condition for some c b > 0, independent of h, .
Theorem 4 (Convergence of the finite element method).
Let Ω ⊂ R d be a polygonal domain in two dimensions or a polyhedral Lipschitz domain in three dimensions. Suppose that (A1)-(A5) hold and there exists an isolated local minimizer (u, ρ) ∈ H 1 g,div (Ω) d × C γ of (BP). Moreover, assume that, for θ > 0, U θ is the subset of Ω where |u| 2 ≥ θ a.e. inŪ θ and suppose that there exists a θ > 0 such thatŪ θ is closed and has non-empty interior U θ for all θ ≤ θ . Let p denote the unique Lagrange multiplier associated with (u, ρ) such that (u, ρ, p) satisfy the first-order optimality conditions (FOC1)-(FOC3).
Consider the conforming finite element spaces X h ⊂ H 1 (Ω) d , C γ,h ⊂ C γ , and M h ⊂ L 2 0 (Ω) and suppose that the assumptions (F1)-(F3) hold. Then, there exists anh > 0 such that, forh ≥ h → 0, there is a sequence of solutions (u h , ρ h , p h ) ∈ X g h ,h ×C γ,h ×M h to the following discretized first-order optimality conditions In Proposition 4, by fixing a ball around an isolated local minimizer, we show that finite element minimizers of a modified optimization problem converge weakly in H 1 (Ω) d × L 2 (Ω) to the isolated minimizer of the infinite-dimensional problem. From this we deduce that there exists a subsequence of finite element solutions (u h ) that converges strongly to the isolated minimizer of the infinitedimensional problem in L 2 (Ω) d . We then strengthen the convergence of ρ h to strong convergence in L s (Ω), s ∈ [1, ∞), in Proposition 5 and strengthen the convergence of u h to strong convergence in H 1 (Ω) d in Proposition 6. In Proposition 7, we prove that there exists anh > 0 such that there is a subsequence, h > h → 0, of strongly converging finite element solutions that also satisfy discretized first-order optimality conditions. Finally, in Proposition 8, we show that the Lagrange multiplier, p h ∈ M h , that satisfies the discretized first-order optimality conditions, converges strongly in L 2 (Ω) to the Lagrange multiplier for the infinite-dimensional problem.
We now fix an isolated minimizer (u, ρ) of (BP). We define the radius of the basin of attraction as the largest value r such that (u, ρ) is the unique local minimizer in B r, We also define B r,H 1 (Ω) (u) and B r,L 2 (Ω) (ρ) by We note that and hence (u, ρ) is also the unique minimizer in (H 1 g,div (Ω) d ∩ B r/2,H 1 (Ω) (u)) × (C γ ∩ B r/2,L 2 (Ω) (ρ)).
Moreover, we define the spaces V g h ,h and V 0,h by Remark 7. In the context of the material distribution, it might be a more natural choice to assume that ρ is isolated with respect to the L ∞ -norm. Assuming ρ is isolated with respect to the L 2 -norm is a stronger isolation assumption, as it cannot be guaranteed that if η ∈ C γ lives in an isolated neighborhood with respect to the L 2 -norm, η ∈ B r,L 2 (Ω) (ρ), then, there exists an r * > 0, such that η ∈ B r * ,L ∞ (Ω) (ρ) where B r * ,L ∞ (Ω) (ρ) is an isolated neighborhood with respect to the L ∞ -norm. We make this stronger isolation assumption as simple and continuous functions are not dense in L ∞ (Ω), but are dense in L 2 (Ω). This has implications in the assumption (F3) and, subsequently, in the remaining results. However, as far as we are aware, the L 2 -isolation assumption is valid for all practical problems found in the literature, in particular it holds for both examples found in Section 5.
Remark 8. We note that balls centred at ρ ∈ C γ are equivalent if measured against any L s -norm for s ∈ [1, ∞). More precisely, for any s, q ∈ [1, ∞), if η ∈ C γ ∩ B r,L q (Ω) (ρ), there exists an r * > 0, depending on s ∈ [1, ∞), such that η ∈ C γ ∩ B r * ,L s (Ω) (ρ). Hence, the assumption that the material distribution is isolated with respect to L 2 -norm is not a stronger assumption than being isolated with respect to the L s -norm provided s ∈ [1, ∞).
In Propositions 4-7 and Corollary 1, we fix an isolated minimizer (u, ρ) of (BP) and suppose that the conditions of Theorem 4 hold.
Hence we can extract subsequences, (u h ) and (ρ h ) of the sequence generated by the global minimizers of (BP h ) such that By assumption (F3), there exists a sequence of finite element functionsρ h ∈ C γ,h that strongly converges to ρ in L 2 (Ω). Moreover letũ h ∈ V g h ,h be a finite element function taken from the sequence of finite element functions that satisfỹ u h → u strongly in H 1 (Ω) d . Such a sequence is shown to exist in [6, Lemma 3.1].
Corollary 1 (Strong convergence of ρ h in L s (Ω b ) ). Let Ω b be any measurable subset of Ω of positive measure on which ρ is equal to zero or one a.e. (if such a set exists). Then, there exists a sequence of finite element minimizers, ρ h , of (BP h ) that converge strongly in L s (Ω b ) to the isolated local minimizer of (BP), where s ∈ [1, ∞).
which implies that ρ h → ρ strongly in L s (Ω).
There exists a subsequence of minimizers, (u h ), of (BP h ) such that Proof. We note that the set W h := V g h ,h ∩ B r/2,H 1 (Ω) (u) is convex. Hence, by following the same arguments for deriving the variational inequality on ρ as in Proposition 2, it can be shown that minimizers of (BP h ) satisfy the variational inequality We note that in the next proposition (Proposition 7), we will show that (4.37) can be strengthened to an equality. However, this result is not currently available at this point and an equality does not follow from the arguments in Proposition 2. We note though that an inequality is sufficient for the subsequent arguments. From Proposition 2 we deduce that, for all w h ∈ W h , Hence By subtracting a ρ h (w h , u h − w h ) from both sides we see that We note that a ρ h is coercive with constant c a = ν/(c 2 p + 1), and b is bounded with constant C b . Hence, Hence, where C = C (ᾱ, ν, c a , C b , L α ). For sufficiently small h, we note thatũ h ∈ W h (whereũ h is defined in the proof of Proposition 4) andũ h → u strongly in H 1 (Ω) d . Moreover by assumption (F3), there exists a sequence of finite element functionsp h ∈ M h that converges to p strongly in L 2 (Ω). Suppose w h =ũ h and q h =p h . From Proposition 5, we know that there exists a subsequence (not indicated) such that ρ h → ρ strongly in L 4 (Ω). We now observe that where L α is the Lipschitz constant for α. Hence by taking the limit as h → 0, we deduce that u h → u strongly in H 1 (Ω) d .
In the following proposition, we show that (up to a subsequence) minimizers of (BP h ) also satisfy the first-order optimality conditions that are the finite-dimensional analogue of the first-order optimality conditions associated with (BP). This allows us to consider the finite-dimensional optimization problem over the whole set V g h ,h × C γ,h , rather than the restricted set (V g h ,h ∩ B r/2,H 1 (Ω) (u)) × (C γ,h ∩ B r/2,L 2 (Ω) (ρ)).
Proposition 7 (Discretized first-order optimality conditions). There exists an h > 0 such that for all h <h, there exists a unique Lagrange multiplier p h ∈ M h such that the functions (u h , ρ h ) that locally minimize (BP h ) satisfy the firstorder optimality conditions (FOC1 h )-(FOC3 h ) Proof. From Proposition 6, we know that u h → u strongly in H 1 (Ω) d . Hence by definition of strong convergence, there exists anh > 0 such that, for all h ≤h, . Now we can follow the reasoning of the proof of Proposition 2 (adding the subscript h where necessary) to deduce the existence of a unique p h ∈ M h such that (FOC1 h )-(FOC3 h ) hold.
Proposition 8 (Strong convergence of p h in L 2 (Ω)). There is a subsequence of the unique p h ∈ M h defined in Proposition 7 that converges strongly in L 2 (Ω) to the p ∈ L 2 0 (Ω) that solves (FOC1)-(FOC3) for the given isolated local minimizer (u, ρ).
Proof. The inf-sup condition (F2) for M h and X 0,h implies that, for any q h ∈ M h , Hence, where C = C(c b , C b ,ᾱ, L α , ν). By assumption (F3), there exists a sequence of finite element functions,p h ∈ M h that satisfiesp h → p strongly in L 2 (Ω). Let q h =p h . We have already shown that u h → u strongly in H 1 (Ω) d in Proposition 6. Similarly, in the proof of Proposition 6 we also showed that (α(ρ) − α(ρ h ))u L 2 (Ω) → 0. Hence we conclude that p h → p strongly in L 2 (Ω).
Proof of Theorem 4. Fix an isolated minimizer (u, ρ) of (BP) and its unique associated Lagrange multiplier p. By the results of Propositions 4, 5, 6, and 7, there exists a mesh sizeh such that for, h <h, there exists a sequence of finite that converges to (u, ρ, p). By taking a subsequence if necessary (not indicated), Proposition 5, implies that ρ h → ρ strongly in L s (Ω), s ∈ [1, ∞), Proposition 6 implies that u h → u strongly in H 1 (Ω) d , and Proposition 8 implies that p h → p strongly in L 2 (Ω).

Numerical results
The main goal of this section is to experimentally verify the existence of strongly converging sequences that were proven to exist in Section 4. In all examples the systems are discretized with the finite element method using FEn-iCS [28]. The computational domains are triangulated with simplices and we define the mesh size, h, as the maximum diameter of all the simplices in the triangulation. The solutions are computed using the deflated barrier method [29]. The deflated barrier method reformulates the discretized first-order optimality conditions (FOC1 h )-(FOC3 h ) as a mixed complementarity problem. The mixed complementarity problem is then solved with a primal-dual active set strategy [30], a Newton-like solver that incorporates the box constraints on ρ h . Typically, nonlinear convergence is difficult to achieve from a naïve initial guess. Hence, the functional J in (BP) is augmented with barrier-like terms. Continuation of the barrier parameter to zero then recovers the solution to the original first-order optimality conditions. A key feature of the deflated barrier method, as required in this work, is the ability to systematically discover multiple solutions of topology optimization problems from the same initial guess. This is achieved via the deflation technique [31,32], [29,Sec. 3.2]. Deflation prevents a Newton-like solver from converging to an already discovered solution by modifying the discretized first-order optimality conditions with a deflation operator. Deflation is extremely cheap to implement (effectively at the same cost as two inner products) and does not affect the conditioning of the linear systems that are solved during the run of the Newton-like solver. The resulting linear systems arising in the deflated barrier method are solved by a sparse LU factorization with MUMPS [33] and PETSc [34].
There are no known solutions for choices of the inverse permeability, α, used in practice. Hence, errors are measured with respect to a heavily-refined finite element solution, which is constructed as follows; first the finite element solutions are computed on a mesh with mesh size h = 0.028, using the deflated barrier method. Next, the mesh is adaptively refined three times in areas where the material distribution is between 1/10 and 9/10. Each time the mesh is refined, the coarse-mesh solution is interpolated onto the finer mesh as an initial guess and the first-order optimality conditions are re-solved using the deflated barrier method.
In principle, there can be infinitely many different subsequences of finite element solutions that strongly converge to the same minimizer of the infinitedimensional problem at different convergence rates. Separate subsequences cause difficulties in the interpretation of the convergence plots as they present themselves as oscillations in the error. This is observed in practice and appears to be caused by at least the following two reasons: (P1) Multiple finite element solutions can exist on the same mesh that represent the same solution of the infinite-dimensional problem, e.g. Figure 5; (P2) A fine mesh can align worse than a coarser mesh with the jumps in the material distribution that solves the infinite-dimensional problem.
Observation (P1) is not surprising in the context of nonlinear PDEs and nonconvex variational problems. In such cases, an additional selection mechanism is required in order to favor one particular solution over others coexisting on the same mesh. Selection mechanisms are problem-dependent. In the case of nonlinear hyperbolic conservation laws the entropy condition plays this role. In the present context, one might propose choosing the solution, minimizing the modified optimization problem (BP h ), that attains the smallest objective functional value for J, within the basin of attraction of the isolated local minimizer. For sufficiently small h, a minimizer satisfying this selection mechanism must exist. However, it is not necessarily unique and numerically enforcing such a condition can be difficult. In order to promote convergence to the minimizer of (BP h ) with the smallest value J, we interpolate the heavily-refined finite element solutions onto coarser meshes as initial guesses for the deflated barrier method. This strategy was effective in practice. The effects of the second observation (P2) are harder to test. However, in Section 5.2, we attempt to minimize mesh bias by measuring errors on unstructured meshes.
Code availability: For reproducibility, an implementation of the deflated barrier method as well as scripts to generate the convergence plots and solutions can be found at https://bitbucket.org/papadopoulos/deflatedbarrier/. The version of the software used in this paper is archived on Zenodo [35].
The inverse permeability, α, is as given in (2.1), with α = 2.5×10 4 and q = 1/10, which satisfies (A1)-(A5). Here q is a penalty parameter which controls the level of intermediate values (between zero or one) in the optimal design. Figure 1 depicts the material distribution of three minimizers. One local minimizer is in the shape of a figure eight and the two Z 2 symmetric global minimizers are in the shape of annuli. Since the domain is convex, g = 0, and f ∈ L 2 (Ω) d , then, by the regularity results proven in the Appendix, u ∈ H 2 (Ω) 2 and p ∈ H 1 (Ω). The conditions of Theorem 3 hold and, therefore, ρ ∈ H 1 (U θ ) for every θ > 0. In this particular example, the support of ρ is compactly contained in the support of the velocity in all three solutions. Therefore, we conclude that ρ ∈ H 1 (Ω). Consider a (P 2 ) 2 × P 1 Taylor-Hood finite element discretization for the velocity and the pressure, and a P 0 piecewise constant finite element discretization for the material distribution. Since all three solutions are isolated local minimizers, by Theorem 4, there exists a sequence of finite element solutions to the discretized first-order optimality conditions that strongly converges to the figure eight solution, and different sequences of different finite element solutions that strongly converge to the two annulus solutions. Their existence is confirmed in Figure 2.
Since ρ ∈ H 1 (Ω) and we are using a P 0 finite element discretization, a naïve prediction for the convergence rate of the L 2 -norm error of the material distribution is O(h). This rate is observed in the bottom left panel of Figure  2. Moreover, since the minimum regularity of the velocity is u ∈ H 2 (Ω) d , and we are using a (P 2 ) 2 finite element discretization, a prediction for the expected convergence rates of the velocity are O(h) and O(h 2 ) for the H 1 -norm and L 2 -norm errors of the velocity, respectively.
In the standard Stokes system, the regularity of u is related to the regularity of the forcing term f ∈ H s (Ω) d , such that u ∈ H s+2 (Ω) d (assuming the domain and boundary data are also suitably regular). Here, the regularity of the forcing term satisfies s < 1/2. If we assume that the velocity has the additional regularity u ∈ H s+2 (Ω) d , s ∈ (0, 1/2), in this context, a prediction for the upper limit of the convergence rate is O(h r ) and O(h t+1 ), for some r, t ∈ [1, s + 1], for the H 1 -norm and L 2 -norm errors of the velocity, respectively. The rates observed in the top panels of Figure 2 match this prediction. The H 1 -norm error is decreasing at a rate slightly faster than O(h 3/2 ) for all three solutions and the L 2 -norm error convergence rate is O(h 2 ) for the figure eight solution and O(h 5/2 ) for the annuli solutions. We hypothesize that the upper limit of the convergence rate of the L 2 -norm error of the velocity is bounded by the relatively slower rate of the convergence of the material distribution.
Finally, since the minimum regularity of the pressure is p ∈ H 1 (Ω) and the discretization is P 1 , a prediction for the convergence rate of the L 2 -norm error is O(h r ), for some r ∈ [1, s + 1]. Initially, the convergence rate is O(h 3/2 ) which matches our naïve prediction. However, on finer meshes, the convergence rate increases. We hypothesize that this speedup is artificial and is caused by the lack of resolution of the refined finite element solutions that are being used as proxies for the solutions of the infinite-dimensional problem in the error norm estimate. Qualitatively, it can be checked that mesh refinement in areas where the discretized material distribution lies between 1/10 and 9/10 is an ineffective strategy for improving the approximation of the pressure that solves the infinitedimensional problem over the whole domain.  Figure 2: The convergence of u h , ρ h , and p h in the discontinuous-forcing problem for the figure eight and annulus solutions on structured meshes, with a P 0 × (P 2 ) 2 × P 1 discretization for (ρ h , u h , p h ).  The function α is as given in (2.1), with α = 2.5 × 10 4 and q = 1/10. Figure  3 visualizes the setup of the problem and depicts the material distribution of the two minimizers. One local minimizer is a straight channel solution and the global minimizer is in the form of a double-ended wrench. We employ a (P 2 ) 2 × P 1 Taylor-Hood finite element discretization for the velocity and the pressure, and a P 1 continuous piecewise linear finite element discretization for the material distribution. This example satisfies all the conditions of Theorem 4 and, for both minimizers, we numerically verify that there exists a sequence of finite element solutions that strongly converges to it in Figure 4.

Double-pipe
As mentioned earlier in this section, it may be possible to find a sequence of mesh sizes, (h i ), such that there exist two different sequences of finite element solutions that strongly converge to the same isolated minimizer. In Figure 5, we depict two different straight channel finite element solutions that exist on the same unstructured mesh where h = 0.04. Both solutions satisfy the discretized first-order optimality conditions (FOC1 h )-(FOC3 h ) and both locally minimize J(v h , η h ). Choosing one over the other would change the convergence pattern of the strongly converging sequence. This may cause difficulty in practice, as optimization strategies are unlikely to discover the discretized global minimum without additional selection mechanisms.

Conclusions
In this work we have studied the fluid topology optimization model of Borrvall and Petersson [6]. In the case of a homogeneous Dirichlet boundary condition and under a mild convexity assumption on the inverse permeability term, α, we have shown that the material distribution necessarily lives in the Sobolev space H 1 inside any compact subset of the support of the velocity. Moreover, we have formally treated the nonconvexity of the optimization problem (including the case of inhomogeneous Dirichlet boundary conditions) and have shown that, given an isolated minimizer of the infinite-dimensional problem, there exists a sequence of discretized solutions, satisfying the associated first-order   optimality conditions, that strongly converges to the minimizer in the appropriate spaces. We have numerically verified that these sequences exist and have discussed the observed convergence rates.
Appendix. Elliptic regularity of the generalized Stokes system with Dirichlet boundary data Lemma 5. Let the domain Ω be either a convex polygon in two dimensions or a convex polyhedron in three dimensions and consider the triple (u, ρ, p) ∈ H 1 g (Ω) d × C γ × L 2 0 (Ω) that satisfies (FOC1)-(FOC3). Suppose that the forcing term f ∈ L 2 (Ω) d and the boundary datum g is the boundary trace of a function g ∈ H 2 (Ω) d on the boundary ∂Ω and satisfies ∂Ω g·n ds = 0. Then u ∈ H 2 (Ω) d and p ∈ H 1 (Ω).
Proof. The idea of the proof is to reduce the system (FOC1)-(FOC2) to the standard Stokes system with a homogeneous Dirichlet boundary condition and invoke the regularity results of Kellogg and Osborn [36] and in three dimensions the results found in Kozlov et al. [37] and Maz'ya and Shaposhnikova [38].