Numerical analysis of a topology optimization problem for Stokes flow
NNumerical analysis of a topology optimization problemfor Stokes flow
I. P. A. Papadopoulos a,1, ∗ , E. S¨uli a a Mathematical Institute, University of Oxford, Oxford, OX2 6GG, UK
Abstract
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 andshowed that a suitably chosen finite element method will converge in a weak(-*) sense to an unspecified solution. In this work, we prove novel regularityresults and extend their numerical analysis. In particular, given an isolatedlocal minimizer to the analytical problem, we show that there exists a sequenceof finite element solutions, satisfying necessary first-order optimality conditions,that strongly converges to it. We also provide the first numerical investigationinto convergence rates.
Keywords: topology optimization, Stokes flow, regularity, finite elementmethod, nonconvex variational problem, multiple solutions
1. Introduction
Topology optimization has become an effective technique in structural andadditive manufacturing and has found multiple uses in medicine, architecture,and engineering [3, 20, 24]. The objective is to find the optimal distribution ofa fluid or solid within a given domain that minimizes a problem-specific costfunctional [5, 10]. In this paper we consider a model for topology optimizationfor fluids proposed by Borrvall and Petersson [11]. Their seminal work hasbecome the foundation for a number of developments in recent years [4, 6,7, 26, 17, 15, 14, 12, 19, 2, 23, 30]. Their goal was to minimize the powerdissipation of a fluid that satisfies both the Stokes equations and a volumeconstraint restricting the proportion of the domain that the fluid can occupy. ∗ Corresponding author
Email addresses: [email protected] (I. P. A. Papadopoulos), [email protected] (E. S¨uli) I. P. is supported by the EPSRC Centre for Doctoral Training in Partial DifferentialEquations: Analysis and Applications [grant number EP/L015811/1] and The MathWorks,Inc.
Preprint submitted to arXiv February 23, 2021 a r X i v : . [ m a t h . NA ] F e b n their paper, they derived generalized Stokes equations , which incorporatethe classical velocity and pressure terms but also introduce a variable, ρ , thatrepresents the material distribution of the fluid over the given domain. Thepresence of fluid is indicated by a value of one in the material distributionwhereas absence of fluid is represented by a value of zero. It would be ideal for ρ : Ω → { , } , in order to remove any ambiguity in the solutions, however, ingeneral, this is a numerically intractable objective. In the Borrvall–Peterssonmodel ρ : Ω → [0 , α , which favors solutions where ρ is close to zero or one. From thegeneralized Stokes equations, Borrvall and Petersson formulated an infinite-dimensional nonconvex optimization problem with inequality, PDE and boxconstraints.In the original paper [11], it is shown that a minimizing velocity and mate-rial distribution to the optimization problem exist [11, Th. 3.1]; however, theminimizer is not necessarily unique [11, Sec. 4.5]. It is also shown that thereexist finite element solutions that converge to a minimizer of the problem [11,Th. 3.2]. The proven convergence is weak in the approximation of the velocityand weak-* in the material distribution, with no results for the pressure. Inaddition, Borrvall and Petersson show that the approximation of material dis-tribution strongly converges to a solution in L s (Ω b ), s ∈ [1 , ∞ ), where Ω b is anymeasurable subset of Ω in which the analytical material distribution is equalto zero or one a.e. [11, Sec. 3.3]. Weak-* convergence permits large oscillationsin the material distribution, called checkerboarding, which could occur in areaswhere the material distribution is not zero or one under the current results.However, in practice, checkerboarding is not observed in these regions. Sincethere can be multiple solutions, the nature of the convergence is ambiguous.In particular, it is not clear if there are sequences of finite element solutionsconverging to every analytical solution.Our goal is to extend and refine the analysis of Borrvall and Petersson. Weshow that, given an isolated local minimizer to the analytical problem, thereexists a sequence of finite element solutions, satisfying the necessary first-orderoptimality conditions, that strongly converges. In particular, we strengthen theconvergence from weak convergence in H (Ω) d to strong covergence in H (Ω) d for the velocity, and from weak-* convergence in L ∞ (Ω) to strong convergencein L s (Ω), s ∈ [1 , ∞ ) for the material distribution. Moreover, in the case of ahomogeneous Dirichlet boundary condition, we show that the material distri-bution is weakly differentiable inside any compact subset of the support of thevelocity; more specifically ρ ∈ H ( U θ ), for any θ >
0, where U θ is any measur-able subset of Ω in which | u | ≥ θ > U θ . This analysis confirms thatcheckerboarding cannot occur under mild assumptions on the model. Hence,isolated local minimizers of the problem can be well approximated using thefinite element method. We conclude with a numerical investigation into theconvergence 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 Stokesequations are satisfied but we also show that the material distribution satisfies a2ariational inequality. In Section 3, we show that under moderate assumptions,the material distribution is weakly differentiable in the case of a homogeneousDirichlet boundary condition for the velocity. We tackle the issue of multiplelocal minima in Section 4, by considering closed balls around isolated localminimizers. In that section, we also prove that for each isolated local minimizerthere exists a sequence of finite element solutions to the discretized first-orderoptimality conditions, which strongly converges to the analytical solution. InSection 5, we computationally investigate the convergence of sequences of finiteelement approximations to their respective analytical solutions.
2. Existence and necessary first-order optimality conditions
The topology optimization problem of Borrvall and Petersson [11] is as fol-lows: find the velocity, u , and the material distribution, ρ , that solve the mini-mization problem min ( v ,η ) ∈ H g , div (Ω) d × C γ J ( v , η ) , (BP)where, J ( v , η ) := 12 (cid:90) Ω (cid:0) α ( η ) | v | + ν |∇ v | − f · v (cid:1) d x,H g (Ω) d := { v ∈ H (Ω) d : v | ∂ Ω = g on ∂ Ω } ,H g , div (Ω) d := { v ∈ H g (Ω) d : div( v ) = 0 a.e. in Ω } ,C γ := (cid:26) η ∈ L ∞ (Ω) : 0 ≤ η ≤ , (cid:90) Ω η d x ≤ γ | Ω | (cid:27) , where Ω ⊂ R d is a Lipschitz domain with dimension d = 2 or d = 3, f ∈ L (Ω) d , ν is the (constant) viscosity, and γ ∈ (0 ,
1) is the volume fraction. Here, α is theinverse permeability, modeling the influence of the material distribution on theflow. For values of ρ close to one, α ( ρ ) is small permitting fluid flow; for smallvalues of ρ , α ( ρ ) is very large, restricting fluid flow. The function α is assumedto have the following properties:(A1) α : [0 , → [ α, α ] with 0 ≤ α and α < ∞ ;(A2) α is convex and monotonically decreasing;(A3) α (0) = α and α (1) = α ,generating a superposition operator also denoted α : C γ → L ∞ (Ω; [ α, α ]). Typ-ically, in the literature α takes the form [11] α ( ρ ) = ¯ α (cid:18) − ρ ( q + 1) ρ + q (cid:19) , (2.1)where q > q →∞ α ( ρ ) = ¯ α (1 − ρ ). Fur-thermore, | ∂ Ω is to be understood in the boundary trace sense [13, Ch. 5.5],3 ∈ H / ( ∂ Ω) d , g = on Γ ⊂ ∂ Ω, with H d − (Γ) >
0, i.e. Γ has nonzero( d − c p such that (cid:107) v (cid:107) L (Ω) ≤ c p (cid:107)∇ v (cid:107) L (Ω) for all v ∈ H (Ω) d with v | Γ = . Remark 1.
The integral in (BP) is well defined. Indeed, since α is assumed tobe convex it is Borel measurable; also since ρ ∈ C γ is Lebesgue measurable, thecomposition α ( ρ ) : Ω → [ α, α ] is Lebesgue measurable. The next theorem is due to Borrvall and Petersson [11, Th. 3.1].
Theorem 1.
Suppose that Ω ⊂ R d is a Lipschitz domain, with d = 2 or d = 3 and α satisfies properties (A1)–(A3). Suppose in addition that α ∈ C ([0 , α, ¯ α ]) . Then, there exists a pair ( u , ρ ) ∈ H g , div (Ω) d × C γ that locallyminimizes J , as defined in (BP). Proposition 1.
Suppose α also satisfies(A4) α is twice continuously differentiable.Then J : H (Ω) d × L s (Ω) → R is Fr´echet differentiable with respect to u and ρ , where < s ≤ ∞ in two dimensions and ≤ s ≤ ∞ in three dimensions.Moreover, for all v ∈ H (Ω) d and η ∈ C γ we have that (cid:104) J (cid:48) u ( u , ρ ) , v (cid:105) = (cid:90) Ω α ( ρ ) u · v + ν ∇ u : ∇ v − f · v d x, (2.2) (cid:104) J (cid:48) ρ ( u , ρ ) , η − ρ (cid:105) = 12 (cid:90) Ω α (cid:48) ( ρ ) | u | ( η − ρ ) d x, (2.3) where J (cid:48) u ( u , ρ ) denotes the Fr´echet deriviative of J with respect to u and J (cid:48) ρ ( u , ρ ) denotes the Fr´echet deriviative of J with respect to ρ .Proof. Proposition 1 follows from the definition of Fr´echet differentiability, themean value theorem, and the Sobolev embedding theorem.
Remark 2.
It can be checked that if α is ( n +1) -times continuously differentiablethen J is n -times Fr´echet differentiable with respect to u and ρ . In the result that follows, we are required to distinguish between differenttypes of local minimizers.
Definition 1 (Isolated local minimizer) . Let Z be a Banach space. The function z ∈ Z is an isolated local minimizer of the functional J : Z → R if there existsan open neighborhood E of z such that there are no other minimizers containedin E . The following proposition is the main result of this section. We show thatif ( u , ρ ) is an isolated local minimizer of the optimization problem (BP), thenthe minimizer also satisfies first-order optimality conditions consisting of twoequations and a variational inequality. 4 roposition 2. Suppose that Ω ⊂ R d is a Lipschitz domain, with d = 2 or d = 3 and α satisfies properties (A1)–(A4). By Theorem 1 there exists a localminimizer ( u , ρ ) ∈ H g , div (Ω) d × C γ . If the local minimizer is also an isolatedlocal minimizer, then there exists a unique Lagrange multiplier p ∈ L (Ω) suchthat the following necessary first-order optimality conditions hold: a ρ ( u , v ) + b ( v , p ) = l f ( v ) for all v ∈ H (Ω) d , (FOC1) b ( u , q ) = 0 for all q ∈ L (Ω) , (FOC2) c u ( ρ, η − ρ ) ≥ for all η ∈ C γ , (FOC3) where L (Ω) := (cid:26) q ∈ L (Ω) : (cid:90) Ω q d x = 0 (cid:27) , and a ρ ( u , v ) := (cid:90) Ω [ α ( ρ ) u · v + ν ∇ u : ∇ v ] d x, l f ( v ) := (cid:90) Ω f · v d x,b ( v , q ) := − (cid:90) Ω q div( v ) d x, c u ( ρ, η ) := 12 (cid:90) Ω α (cid:48) ( ρ ) η | u | d x. Proof.
We will first show that (FOC1)–(FOC2) are satisfied by generalizingarguments, used for the Stokes system with a homogeneous Dirichlet boundarycondition, found in [27]. For ease of notation we define X g := H g (Ω) d , X := H (Ω) d , V g := H g , div (Ω) d , V := H , div (Ω) d and M := L (Ω). The respectivedual spaces of X , V and M are denoted with ∗ . We also define the associatedoperators, A ∈ L ( X g , X ∗ ), B ∈ L ( X g , M ) and B ∈ L ( X , M ) by (cid:104) A u , v (cid:105) := a ρ ( u , v ) , (cid:104) B w , q (cid:105) := b ( w , q ) , and (cid:104) B v , q (cid:105) := b ( v , q ) . (2.4)We note that ker( B ) = V . 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 , wesee that u + t v ∈ V g , t ∈ R . If ( u , ρ ) ∈ V g × C γ is an isolated local minimum,then, by definition, there exists an r > w , η ) ∈ V g × C γ ,( w , η ) (cid:54) = ( u , ρ ) that satisfies (cid:107) u − w (cid:107) H (Ω) + (cid:107) ρ − η (cid:107) L ∞ (Ω) ≤ r (2.5)we have that J ( u , ρ ) < J ( w , η ). Hence, for any given v ∈ V , if 0 < t ≤ r/ (cid:107) v (cid:107) H (Ω) , the following inequality holds1 t ( J ( u + t v , ρ ) − J ( u , ρ )) ≥ . (2.6)By Proposition 1, J is Fr´echet differentiable, and therefore also Gateaux differ-entiable, with respect to u . Hence as t → + , we see that (cid:104) J (cid:48) u ( u , ρ ) , v (cid:105) ≥ v ∈ V . (2.7)5y considering the same reasoning with t <
0, we deduce that (cid:104) J (cid:48) u ( u , ρ ) , v (cid:105) = 0 for all v ∈ V . (2.8)From Proposition 1, we know that J (cid:48) u ( u , ρ ) = A u − f and hence A u − f ∈ V ◦ where V ◦ := (ker( B )) ◦ = { h ∈ X ∗ : (cid:104) h, v (cid:105) = 0 for all v ∈ V } . (2.9)We know that the operator B satisfies the following equivalent version of theinf-sup condition [18, Ch. 1, Sec. 4.1, Lem. 4.1]:there exists a β > q ∈ M, (cid:107) B ∗ q (cid:107) X ∗ ≥ β (cid:107) q (cid:107) M , (2.10)where B ∗ is the dual operator of B , defined by (cid:104) v , B ∗ q (cid:105) = (cid:104) B v , q (cid:105) . Thisimplies that B ∗ is injective (and therefore bijective) from M into Im B ∗ . Fur-thermore, it also implies that ( B ∗ ) − is continuous. Consider f ∈ Im B ∗ ; then,there exists a q ∈ M such that f = B ∗ q and (cid:107) ( B ∗ ) − f (cid:107) M ≤ β (cid:107) f (cid:107) X ∗ . (2.11)Therefore, Im B ∗ is closed.Since Im B ∗ is closed, by Banach’s closed range theorem, we know thatIm B ∗ = (ker( B )) ◦ = V ◦ . Hence since, A u − f ∈ V ◦ , there exists a p ∈ M suchthat A u + B ∗ p = f . (2.12)Since B ∗ is injective, p is also unique. Since u ∈ V g , we have that B u = 0.Hence (FOC1) and (FOC2) hold.We will now show that (FOC3) holds. We note that C γ is a convex subsetof a linear space. For any given ζ, η ∈ C γ and t ∈ [0 , ζ + t ( η − ζ ) ∈ C γ . Since ( u , ρ ) is a local minimizer, it follows that for each η ∈ C γ , if 0 < t ≤ r/ (cid:107) η − ρ (cid:107) L ∞ (Ω) , with r as in (2.5), then1 t ( J ( u , ρ + t ( η − ρ )) − J ( u , ρ )) ≥ . (2.13)From Proposition 1, we know that J is Fr´echet differentiable, and therefore alsoGateaux differentiable, with respect to ρ . Hence, by taking the limit as t → c u ( ρ, η − ρ ) = (cid:104) J (cid:48) ρ ( u , ρ ) , η − ρ (cid:105) ≥ η ∈ C γ . (2.14)Therefore (FOC3) holds.The following lemma will be used in the proof of the next proposition.6 emma 2. Consider a nontrivial function η ∈ C γ and the measurable non-empty set E ⊂⊂ supp( η ) , where supp denotes the support of a function, i.e. η > a.e. in E . Then, there exists an (cid:15) (cid:48) > such that, for all (cid:15) ∈ (0 , (cid:15) (cid:48) ] , there existsa set E (cid:15) ⊆ E , | E (cid:15) | > where η > (cid:15) a.e. in E (cid:15) .Proof. For a contradiction, suppose that there exists no such (cid:15) (cid:48) such that E (cid:15) (cid:48) exists. This implies that for all n ≥ , | E \ ˆ E n | = 0 , (2.15)where ˆ E n := { ≤ η ≤ /n a.e. in E } . We see that ∅ = E \ ˆ E ⊆ E \ ˆ E ⊆· · · ⊆ E \ ˆ E n ⊆ · · · , i.e. E \ ˆ E n is ascending. By (2.15) we note thatlim n →∞ | E \ ˆ E n | = 0 . (2.16)Moreover, ∪ ∞ n =1 E \ ˆ E n = lim n →∞ E \{ ≤ η ≤ /n a.e. in E } = E \{ η = 0 a.e. in E } = E \ ∅ = E. (2.17)Now we see that0 < | E | = | ∪ ∞ n =1 E \ ˆ E n | = | lim n →∞ E \ ˆ E n | = lim n →∞ | E \ ˆ E n | = 0 , (2.18)where the first equality follows from (2.17), the third equality follows from thecontinuity of the Lebesgue measure, and the fourth equality follows from (2.16).(2.18) is a contradiction and, therefore, such an (cid:15) (cid:48) > E (cid:15) = E (cid:15) (cid:48) for all 0 < (cid:15) ≤ (cid:15) (cid:48) , we conclude that the statement holds for all (cid:15) ∈ (0 , (cid:15) (cid:48) ].The following proposition is a property of isolated local minimizers that willbe useful for the numerical analysis of the finite element method. Proposition 3.
Suppose that Ω ⊂ R d is a Lipschitz domain, with d = 2 or d = 3 , and α satisfies properties (A1)–(A4). Further assume that the localminimizer ( u , ρ ) ∈ H g , div (Ω) d × C γ of (BP) is isolated . Then, supp( ρ ) ⊆ U ,where U := supp( u ) .Proof. By definition of an isolated local minimizer, there exists an r > w , η ) ∈ H g , div (Ω) d × C γ , ( w , η ) (cid:54) = ( u , ρ ) that satisfies (cid:107) u − w (cid:107) H (Ω) + (cid:107) ρ − η (cid:107) L ∞ (Ω) ≤ r, we have that J ( u , ρ ) < J ( w , η ). For a contradiction, suppose that there existsa set E ⊂ Ω, E ∩ U = ∅ , of positive measure, where ρ > E . By Lemma2, there exists an (cid:15) ∈ (0 , r ) such that there exists a set E (cid:15) ⊆ E , | E (cid:15) | > ρ > (cid:15) a.e. in E (cid:15) . Define ˜ ρ as˜ ρ := (cid:40) ρ a.e. in Ω \ E (cid:15) ,ρ − (cid:15) a.e. in E (cid:15) . (2.19)7s ρ ∈ C γ , also ˜ ρ ∈ C γ . We note that (cid:107) ρ − ˜ ρ (cid:107) L ∞ (Ω) = (cid:107) (cid:15) (cid:107) L ∞ ( E (cid:15) ) < r and,therefore, ( u , ˜ ρ ) lies inside the minimizing neighborhood of the ( u , ρ ). However, J ( u , ˜ ρ ) = J ( u , ρ ) as ρ and ˜ ρ only differ on the set E (cid:15) , but u = a.e. in E (cid:15) ⊆ E by assumption. This contradicts the assertion that ( u , ρ ) is an isolated localminimizer.
3. Regularity of ρ In the this section we show that ρ ∈ C γ possesses higher regularity in thecase of a homogeneous Dirichlet boundary condition on u and if α satisfies astronger (but not restrictive) convexity assumption. Theorem 3 (Weak differentiability of ρ ) . Suppose that the domain Ω ⊂ R d isbounded, the boundary is Lipschitz, and that the data g = on ∂ Ω . Considera local minimizer, ( u , ρ ) ∈ H , div (Ω) d × C γ , of (BP) such that u is not thetrivial solution and there exists a closed subset U θ ⊂ Ω with non-empty interioron which | u | is bounded below by a positive constant, θ > . Suppose that(A1)–(A4) hold and that α ∈ C ([0 , is strongly convex, i.e.,(A5) There exists a constant α (cid:48)(cid:48) min > such that α (cid:48)(cid:48) ( y ) ≥ α (cid:48)(cid:48) min > for all y ∈ [0 , .Then, ∇ ρ exists in U θ and ρ ∈ C γ ∩ H ( U θ ) . Remark 3.
The assumption (A5) excludes the case where α is linear. This isconsistent with previous theory, as Borrvall and Petersson [11] showed that if α is linear then ρ is a linear combination of Heaviside functions and thus ρ / ∈ H (Ω) . However, the assumptions (A1)–(A5) do include α ( ρ ) = ¯ α (1 − ρ ( q + 1) / ( ρ + q )) ,where the lower bound in (A5) is α (cid:48)(cid:48) min = 2¯ αq/ ( q + 1) .Proof of Theorem 3. If we can bound the L -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 ( U θ )for 1 ≤ k ≤ d .The variational inequality on ρ states that12 (cid:90) Ω α (cid:48) ( ρ ) | u | ( η − ρ ) d x ≥ . (3.1)We define U ⊂ Ω as U := supp( u ) and fix an open, bounded and connecteddomain ˆΩ such that ˆΩ = Ω if U ⊂⊂ Ω and Ω ⊂⊂ ˆΩ otherwise. In the casewhere U is not a compact subset of Ω, we extend u and ρ by zero to the wholeof R d . Since the trace of u is zero on the boundary, the extension of u by zerolives in H ( ˆΩ) d . Let 0 < | h | < (1 / U, ∂ ˆΩ) and choose k ∈ { , . . . , d } . Wedefine ρ h as ρ h ( x ) = (cid:40) ρ ( x + he k ) for x ∈ Ω − he k , x ∈ R d \ (Ω − he k ) .
8e define the difference quotient, D hk , in the k -th coordinate direction, as D hk ρ ( x ) = ρ ( x + he k ) − ρ ( x ) h , h ∈ R \{ } , x ∈ ˆΩ . Let η = ( ρ h + ρ − h ) /
2. We note that η ∈ C γ , since0 ≤ ρ h ≤
12 a . e . in Ω , and 0 ≤ ρ − h ≤
12 a . e . in Ω , which implies that 0 ≤ η ≤ (cid:90) Ω η d x = 12 (cid:90) Ω ρ h + ρ − h d x = 12 (cid:90) Ω − he k ∩ Ω ρ d x + 12 (cid:90) Ω+ he k ∩ Ω ρ d x ≤ (cid:90) Ω ρ d x ≤ γ | Ω | . If we multiply (3.1) through by 4 and divide by h we see that1 h (cid:90) Ω α (cid:48) ( ρ ) | u | ( ρ h + ρ − h − ρ ) d x ≥ . (3.2)We note that, D − hk ( D hk ρ ) = ρ − ρ − h h − ρ h − ρh − h = ρ h + ρ − h − ρh . Hence, because u is zero outside of Ω, (3.2) is equivalent to (cid:90) ˆΩ α (cid:48) ( ρ ) | u | ( D − hk ( D hk ρ )) d x ≥ . (3.3)In order to obtain a first-order difference quotient, we will perform the finitedifference analogue of integration by parts to shift the D − hk operator from D hk ρ to α (cid:48) ( ρ ) | u | . We note that, by definition, the left-hand side of (3.3) is equal to − h (cid:90) ˆΩ ( α (cid:48) ( ρ ) | u | )( x ) (cid:0) ( D hk ρ )( x − he k ) − ( D hk ρ )( x ) (cid:1) d x, (3.4)which by a change of variables is equal to − h (cid:18)(cid:90) ˆΩ − he k ( α (cid:48) ( ρ ) | u | )( x + he k )( D hk ρ )( x )d x − (cid:90) ˆΩ ( α (cid:48) ( ρ ) | u | )( x )( D hk ρ )( x )d x (cid:19) . We note that U ⊂⊂ ˆΩ and | h | < (1 / U, ∂ ˆΩ), which implies that U ⊂⊂ ˆΩ − he k . Therefore, (cid:90) ˆΩ − he k ( α (cid:48) ( ρ ) | u | )( x + he k )( D hk ρ )( x )d x = (cid:90) U − he k ( α (cid:48) ( ρ ) | u | )( x + he k )( D hk ρ )( x )d x = (cid:90) ˆΩ ( α (cid:48) ( ρ ) | u | )( x + he k )( D hk ρ )( x )d x. (3.5)9herefore, from (3.3)–(3.5) we see that (cid:90) ˆΩ D hk ( α (cid:48) ( ρ ) | u | )( D hk ρ ) d x ≤ . (3.6)Now we wish to rewrite D hk ( α (cid:48) ( ρ ) | u | ) in a form that we can decouple from D hk ρ in order to be able to bound (3.6) above and below. Now, D hk ( α (cid:48) ( ρ ) | u | )( x ) = 1 h (cid:0) α (cid:48) ( ρ ( x + he k )) | u ( x + he k ) | − α (cid:48) ( ρ ( x )) | u ( x ) | (cid:1) = 12 h (cid:0) α (cid:48) ( ρ ( x + he k )) (cid:0) | u ( x + he k ) | − | u ( x ) | (cid:1)(cid:1) + 12 h (cid:0) α (cid:48) ( ρ ( x )) (cid:0) | u ( x + he k ) | − | u ( x ) | (cid:1)(cid:1) + 12 h (cid:0) | u ( x + he k ) | ( α (cid:48) ( ρ ( x + he k )) − α (cid:48) ( ρ ( x ))) (cid:1) + 12 h (cid:0) | u ( x ) | ( α (cid:48) ( ρ ( x + he k )) − α (cid:48) ( ρ ( x ))) (cid:1) = 12 (cid:0) α (cid:48) ( ρ h ) + α (cid:48) ( ρ ) (cid:1) D hk ( | u | ) + 12 (cid:0) | u h | + | u | (cid:1) D hk ( α (cid:48) ( ρ )) . Therefore, from (3.6) we see that (cid:90) ˆΩ (cid:20) (cid:0) | u h | + | u | (cid:1) D hk ( α (cid:48) ( ρ ))+ 12 (cid:0) α (cid:48) ( ρ h ) + α (cid:48) ( ρ ) (cid:1) D hk | u | (cid:21) D hk ( ρ ) d x ≤ . (3.7)Now, 12 (cid:0) | u h | + | u | (cid:1) D hk ( α (cid:48) ( ρ )) + 12 (cid:0) α (cid:48) ( ρ h ) + α (cid:48) ( ρ ) (cid:1) D hk | u | = 1 h (cid:90) dd s (cid:20) α (cid:48) (cid:0) sρ h + (1 − s ) ρ (cid:1) (cid:0) | u h | + | u | (cid:1) + 12 (cid:0) α (cid:48) ( ρ h ) + α (cid:48) ( ρ ) (cid:1) (cid:12)(cid:12) s u h + (1 − s ) u (cid:12)(cid:12) (cid:21) d s = 1 h (cid:90) (cid:2) α (cid:48)(cid:48) (cid:0) sρ h + (1 − s ) ρ (cid:1)(cid:3) d s (cid:124) (cid:123)(cid:122) (cid:125) =: A (cid:0) | u h | + | u | (cid:1) ( ρ h − ρ )+ 12 h (cid:0) α (cid:48) ( ρ h ) + α (cid:48) ( ρ ) (cid:1) (cid:90) (cid:2) (cid:0) s u h + (1 − s ) u (cid:1)(cid:3) d s (cid:124) (cid:123)(cid:122) (cid:125) =: B · ( u h − u ) . Hence from (3.7) we find that12 (cid:90) ˆΩ A ( | u h | + | u | ) | D hk ρ | + ( α (cid:48) ( ρ h ) + α (cid:48) ( ρ )) B · ( D hk u ) D hk ρ d x ≤ . (3.8)10ubtracting the second term on the left-hand side in (3.8) from both sides, takingabsolute values on the right-hand side, using the Cauchy–Schwarz inequality andmultiplying by 2, we see that (cid:90) ˆΩ A ( | u h | + | u | ) | D hk ρ | d x ≤ (cid:90) ˆΩ | B || α (cid:48) ( ρ h ) + α (cid:48) ( ρ ) || D hk u || D hk ρ | d x. (3.9)Furthermore we note that A ≥ α (cid:48)(cid:48) min and B = (cid:90) (cid:2) (cid:0) s u h + (1 − s ) u (cid:1)(cid:3) d s = 2 (cid:20) s u h + (cid:18) s − s (cid:19) u (cid:21) = u h + u . Hence, using Cauchy’s inequality and Young’s inequality, we see that α (cid:48)(cid:48) min (cid:90) U − he k | u h | | D hk ρ | d x + α (cid:48)(cid:48) min (cid:90) U | u | | D hk ρ | d x ≤ (cid:90) ˆΩ | u + u h || α (cid:48) ( ρ h ) + α (cid:48) ( ρ ) || D hk u || D hk ρ | d x ≤ (cid:90) U − he k | u h || α (cid:48) ( ρ h ) + α (cid:48) ( ρ ) || D hk u || D hk ρ | d x + (cid:90) U | u || α (cid:48) ( ρ h ) + α (cid:48) ( ρ ) || D hk u || D hk ρ | d x ≤ (cid:15) (cid:90) U − he k | u h | | D hk ρ | d x + (cid:15) (cid:90) U | u | | D hk ρ | d x + 12 (cid:15) (cid:90) U − he k | α (cid:48) ( ρ h ) + α (cid:48) ( ρ ) | | D hk u | d x + 12 (cid:15) (cid:90) U | α (cid:48) ( ρ h ) + α (cid:48) ( ρ ) | | D hk u | d x. (3.10)By fixing (cid:15) = α (cid:48)(cid:48) min , from (3.10) we see that, α (cid:48)(cid:48) min (cid:90) U | u | | D hk ρ | d x ≤ α (cid:48)(cid:48) min (cid:90) U | u | | D hk ρ | d x + α (cid:48)(cid:48) min (cid:90) U − he k | u h | | D hk ρ | d x ≤ α (cid:48)(cid:48) min (cid:90) ˆΩ | α (cid:48) ( ρ h ) + α (cid:48) ( ρ ) | | D hk u | d x. (3.11)Now | α (cid:48) ( ρ h ) + α (cid:48) ( ρ ) | is bounded above by 4 sup ζ ∈ C γ | α (cid:48) ( ζ ) | which is in-dependent of h . Consider a set ˜Ω ⊂ R d such that ˆΩ ⊂⊂ ˜Ω. We note that u ∈ H ( ˜Ω) d . By applying Theorem 3 in [13, pg. 294], we see that (cid:90) U | u | | D hk ρ | d x ≤ ˜ C (Ω) sup ζ ∈ C γ | α (cid:48) ( ζ ) | ( α (cid:48)(cid:48) min ) (cid:107)∇ u (cid:107) L (˜Ω) ≤ ˆ C (Ω) sup ζ ∈ C γ | α (cid:48) ( ζ ) | ( α (cid:48)(cid:48) min ) (cid:107)∇ u (cid:107) L (Ω) ≤ C < ∞ , (3.12)11here ˜ C , ˆ C and C are constants. The bound is independent of h and k . Because,by hypothesis, there exists a subset U θ ⊂ Ω such that, | u | ≥ θ > U θ ,we see from (3.12) that, because U θ ⊂ U = supp( u ), also θ (cid:90) U θ | D hk ρ | d x ≤ (cid:90) U θ | u | | D hk ρ | d x ≤ (cid:90) U | u | | D hk ρ | d x ≤ C. (3.13)Estimate (3.13) implies thatsup h (cid:107) D hk ρ (cid:107) L ( U θ ) < ∞ . (3.14)From (3.14) we see that there exists a function η k ∈ L ( U θ ) and a subsequence h i → D h i k ρ (cid:42) η k weakly in L ( U θ ) . Finally, we wish to identify η k with ∂ x k ρ . First choose any smooth and com-pactly supported function, φ ∈ C ∞ c ( U θ ). We note that (cid:90) U θ ρ ∂ x k φ d x ≤ C (cid:107) ρ (cid:107) L ∞ ( U θ ) (cid:107) φ (cid:107) W , ∞ ( U θ ) < ∞ . Since ∂ x k φ is compactly supported in U θ , it follows that (cid:90) U θ ρ ∂ x k φ d x = (cid:90) ˆΩ ρ ∂ x k φ d x. Hence (cid:90) U θ ρ ∂ x k φ d x = lim h i → (cid:90) ˆΩ ρ D − h i k φ d x = − lim h i → (cid:90) ˆΩ ( D h i k ρ ) φ d x = − lim h i → (cid:90) U θ ( D h i k ρ ) φ d x = − (cid:90) U θ η k φ d x. Hence η k = ∂ x k ρ a.e. in U θ for k = 1 , . . . , d . Therefore, from (3.12) we see byweak lower semicontinuity that (cid:90) U θ | ∂ x k ρ | d x ≤ C (Ω , sup ζ ∈ C γ | α (cid:48) ( ζ ) | , α (cid:48)(cid:48) min ) , (3.15)for some constant C . We conclude that ρ ∈ H ( U θ ) ∩ C γ , θ >
4. Finite element approximation
We will be approximately solving (FOC1)–(FOC3) by approximating theanalytical solution with finite element functions. Borrvall and Petersson [11,Sec. 3.3] considered a piecewise constant finite element approximation of thematerial distribution coupled with an inf-sup stable quadrilateral finite element12pproximation of the velocity and the pressure. They showed that such approx-imations of the velocity and material distribution converge to an unspecifiedsolution ( u , ρ ) of (BP) in the following sense [11]: u h (cid:42) u weakly in H (Ω) d ,ρ h ∗ (cid:42) ρ weakly-* in L ∞ (Ω) ,ρ h → ρ strongly in L s (Ω b ) , s ∈ [1 , ∞ ) , 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 approximationof the pressure, p . Since there can be multiple local minimizers, it is unclearwhich solution the sequence of finite element solutions is converging to. In thissection, we consider any suitable conforming mixed finite element space suchthat the velocity and pressure spaces are inf-sup stable. We prove that, for every isolated minimizer of (BP), there exist sequences of finite element solu-tions to the discretized first-order optimality conditions that strongly convergeto the minimizer as the mesh size tends to zero. More specifically, we showthat, for each isolated analytical local minimizer, there exist (possibly different)sequences of finite element solutions ( u h , ρ h , p h ) that converge to it strongly in H (Ω) d × L s (Ω) × L (Ω) as h →
0, where s ∈ [1 , ∞ ). We emphasize that theresults hold in the case where the local minima are isolated.Consider the conforming finite element spaces X h ⊂ H (Ω) d , C γ,h ⊂ C γ ,and M h ⊂ L (Ω). Let X ,h := { v h ∈ X h : v h | ∂ Ω = } .In general, it will not be possible to represent the boundary data g exactlyin the velocity finite element space. Hence, for each h , we instead considerboundary data g h (which can be represented) and assume that(F1) g h → g strongly in H / ( ∂ Ω) d .We now define the space X g h ,h := { v h ∈ X h : v h | ∂ Ω = g h } . We will also assumethat:(F2) X ,h and M h satisfy the following inf-sup condition for some c b > c b ≤ inf q h ∈ M h \{ } sup v h ∈ X ,h \{ } b ( v h , q h ) (cid:107) v h (cid:107) H (Ω) (cid:107) q h (cid:107) L (Ω) . (4.1)(F3) The finite element spaces are dense in their respective function spaces, i.e.,for any ( v , η, q ) ∈ H (Ω) d × C γ × L (Ω),lim h → inf w h ∈ X h (cid:107) v − w h (cid:107) H (Ω) = lim h → inf ζ h ∈ C γ,h (cid:107) η − ζ h (cid:107) L (Ω) = lim h → inf r h ∈ M h (cid:107) q − r h (cid:107) L (Ω) = 0 . Theorem 4 (Convergence of the finite element method) . Let Ω ⊂ R d be apolygonal domain in two dimensions or a polyhedral Lipschitz domain in three imensions. Suppose that (A1)–(A5) hold and there exists an isolated localminimizer ( u , ρ ) ∈ H g , div (Ω) d × C γ of (BP). Moreover, assume that, for θ > , U θ is the subset of Ω where | u | ≥ θ a.e. in U θ and suppose that there exists a θ (cid:48) > such that U θ is closed and has non-empty interior for all θ ≤ θ (cid:48) . 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 (Ω) d , C γ,h ⊂ C γ ,and M h ⊂ L (Ω) and suppose that the assumptions (F1)–(F3) hold.Then, there exists an ¯ h > such that, for ¯ h ≥ h → , there is a sequence ofsolutions ( u h , ρ h , p h ) ∈ X g h ,h × C γ,h × M h to the following discretized first-orderoptimality conditions a ρ h ( u h , v h ) + b ( v h , p h ) = l f ( v h ) for all v h ∈ X ,h , (FOC1 h ) b ( u h , q h ) = 0 for all q h ∈ M h , (FOC2 h ) c u h ( ρ h , η h − ρ h ) ≥ for all η h ∈ C γ,h , (FOC3 h ) such that u h → u strongly in H (Ω) d , ρ h → ρ strongly in L s (Ω) , s ∈ [1 , ∞ ) ,and p h → p strongly in L (Ω) as h → . In Proposition 4, by fixing a ball around an isolated local minimizer, weshow that finite element minimizers of a modified optimization problem convergeweakly in H (Ω) d × L (Ω) to the analytical isolated minimizer. From this wededuce that there exists a subsequence of finite element solutions ( u h ) thatconverges strongly to the analytical isolated minimizer in L (Ω) d . We thenstrengthen the convergence of ρ h to strong convergence in L s (Ω), s ∈ [1 , ∞ ),in Proposition 5 and strengthen the convergence of u h to strong convergencein H (Ω) d in Proposition 6. In Proposition 7, we prove that there exists an¯ h > h > h →
0, of strongly converging finiteelement solutions that also satisfy discretized first-order optimality conditions.Finally, in Proposition 8, we show that the Lagrange multiplier, p h ∈ M h , thatsatisfies the discretized first-order optimality conditions, converges strongly in L (Ω) to the analytical Lagrange multiplier.We now fix an isolated local minimizer ( u , ρ ) of (BP). We define the radiusof the basin of attraction as the largest value r such that ( u , ρ ) is the uniquelocal minimizer in B r,H (Ω) × L (Ω) ( u , ρ ) ∩ ( H g , div (Ω) d × C γ ), where B r,H (Ω) × L (Ω) ( u , ρ ):= { v ∈ H (Ω) d , η ∈ C γ : (cid:107) u − v (cid:107) H (Ω) + (cid:107) ρ − η (cid:107) L (Ω) ≤ r } . (4.2)We also define B r,H (Ω) ( u ) and B r,L (Ω) ( ρ ) by B r,H (Ω) ( u ) := { v ∈ H (Ω) d : (cid:107) u − v (cid:107) H (Ω) ≤ r } , (4.3) B r,L (Ω) ( ρ ) := { η ∈ C γ : (cid:107) ρ − η (cid:107) L (Ω) ≤ r } . (4.4)We note that( H g , div (Ω) d ∩ B r/ ,H (Ω) ( u )) × ( C γ ∩ B r/ ,L (Ω) ( ρ )) ⊂ B r,H (Ω) × L (Ω) ( u , ρ ) ∩ ( H g , div (Ω) d × C γ )14nd hence ( u , ρ ) is also the unique minimizer in ( H g , div (Ω) d ∩ B r/ ,H (Ω) ( u )) × ( C γ ∩ B r/ ,L (Ω) ( ρ )).Moreover, we define the spaces V g h ,h and V ,h by V g h ,h := { v h ∈ X g h ,h : b ( v h , q h ) = 0 for all q h ∈ M h } ,V ,h := { v h ∈ X ,h : b ( v h , q h ) = 0 for all q h ∈ M h } . Proposition 4.
Suppose that the conditions of Theorem 4 hold. Fix an isolatedminimizer ( u , ρ ) of (BP). Consider the finite-dimensional optimization problem:find ( u h , ρ h ) that minimizes min ( v h ,η h ) ∈ ( V g h,h ∩ B r/ ,H ( u )) × ( C γ,h ∩ B r/ ,L ( ρ )) J ( v h , η h ) . (BP h ) Then, a minimizer ( u h , ρ h ) of (BP h ) exists and there exist subsequences (up torelabeling) such that u h (cid:42) u weakly in H (Ω) d , (4.5) u h → u strongly in L (Ω) d , (4.6) ρ h ∗ (cid:42) ρ weakly-* in L ∞ (Ω) , (4.7) ρ h (cid:42) ρ weakly in L s (Ω) , s ∈ [1 , ∞ ) . (4.8) Proof.
The functional J is continuous and ( V g h ,h ∩ B r/ ,H (Ω) ( u )) × ( C γ,h ∩ B r/ ,L (Ω) ( ρ )) is a finite-dimensional, closed and bounded set and, therefore,sequentially compact by Weierstrass’ theorem. Hence J obtains its infimumin ( V g h ,h ∩ B r/ ,H (Ω) ( u )) × ( C γ,h ∩ B r/ ,L (Ω) ( ρ )) and therefore, a minimizer( u h , ρ h ) exists.By a corollary of Kakutani’s Theorem, if a Banach space is reflexive thenevery norm-closed, bounded and convex subset of the Banach space is weaklycompact and thus, by the Eberlein–ˇSmulian theorem, sequentially weakly com-pact. It can be checked that H (Ω) d ∩ B r/ ,H (Ω) ( u ) and C γ ∩ B r/ ,L (Ω) ( ρ )are norm-closed, bounded and convex subsets of the reflexive Banach spaces H (Ω) d and L (Ω), respectively. Therefore, H (Ω) d ∩ B r/ ,H (Ω) ( u ) is weaklysequentially compact in H (Ω) d and C γ ∩ B r/ ,L (Ω) ( ρ ) is weakly sequentiallycompact in L (Ω).Hence we can extract subsequences, ( u h ) and ( ρ h ) of the sequence generatedby the global minimizers of (BP h ) such that u h (cid:42) ˆ u ∈ H (Ω) d ∩ B r/ ,H (Ω) ( u ) weakly in H (Ω) d , (4.9) ρ h (cid:42) ˆ ρ ∈ C γ ∩ B r/ ,L (Ω) ( ρ ) weakly in L (Ω) . (4.10)By assumption (F3), there exists a sequence of finite element functions ˜ ρ h ∈ C γ,h that strongly converges to ρ in L (Ω). Moreover let ˜ u h ∈ V g h ,h be a finiteelement function taken from the sequence of finite element functions that satisfy ˜ u h → u strongly in H (Ω) d . Such a sequence is shown to exist in [11, Lemma3.1]. 15e now wish to identify the limits ˆ u and ˆ ρ . Consider the following bound: | J ( ˜ u h , ˜ ρ h ) − J ( u , ρ ) |≤ (cid:90) Ω | ( α ( ρ ) − α (˜ ρ h )) | u | | + | α (˜ ρ h )( | u | − | ˜ u h | ) | d x + (cid:90) Ω ν (cid:12)(cid:12) |∇ u | − |∇ ˜ u h | | (cid:12)(cid:12) + 2 | f · ( u − ˜ u h ) | d x ≤ L α (cid:107) u (cid:107) L (Ω) (cid:107) ˜ ρ h − ρ (cid:107) L (Ω) + ¯ α (cid:107) ˜ u h − u (cid:107) L (Ω) ( (cid:107) ˜ u h − u (cid:107) L (Ω) + 2 (cid:107) u (cid:107) L (Ω) )+ ν (cid:107) ˜ u h − u (cid:107) H (Ω) ( (cid:107) ˜ u h − u (cid:107) H (Ω) + 2 (cid:107) u (cid:107) H (Ω) )+ 2 (cid:107) f (cid:107) L (Ω) (cid:107) ˜ u h − u (cid:107) L (Ω) , (4.11)where L α denotes the Lipschitz constant for α . From (4.11) we see that J ( ˜ u h , ˜ ρ h ) → J ( u , ρ ) as h → . Furthermore, for sufficiently small h > ˜ u h , ˜ ρ h ) ∈ ( V g h ,h ∩ B r/ ,H (Ω) ( u )) × ( C γ,h ∩ B r/ ,L (Ω) ( ρ )) . Therefore, J ( u h , ρ h ) ≤ J ( ˜ u h , ˜ ρ h ) . (4.12)By taking the limit as h → ˜ u h and ˜ ρ h to u and ρ , respectively, we see thatlim h → J ( u h , ρ h ) ≤ J ( u , ρ ) . (4.13)(F1) implies that u h | ∂ Ω = g h → g strongly in H / ( ∂ Ω) d . (4.14)By assumption (F3), for every q ∈ L (Ω), there exists a sequence of ˜ q h ∈ M h such that ˜ q h → q strongly in L (Ω). Since u h (cid:42) ˆ u weakly in H (Ω) d and u h ∈ V g h ,h , we see that b ( ˆ u , q ) = lim h → b ( u h , ˜ q h ) + lim h → b ( u h , q − ˜ q h ) = 0 for all q ∈ L (Ω) . (4.15)Hence ˆ u is pointwise divergence-free and together with (4.14), we deduce thatˆ u ∈ H g , div (Ω) d ∩ B r/ ,H (Ω) ( u ). By construction ˆ ρ ∈ C γ ∩ B r/ ,L (Ω) ( ρ ).Since ( u , ρ ) is the unique local minimizer of ( H g , div (Ω) d ∩ B r/ ,H (Ω) ( u )) × ( C γ ∩ B r/ ,L (Ω) ( ρ )), from (4.13), we see thatlim h → J ( u h , ρ h ) = J ( u , ρ ) . (4.16)16ince ( u , ρ ) is the unique minimizer in the spaces we consider, we can identifythe limits ˆ u and ˆ ρ as u and ρ , respectively, and state that u h (cid:42) u weakly in H (Ω) d and ρ h (cid:42) ρ weakly in L (Ω). By the Rellich–Kondrachov theorem, wecan extract a further subsequence such that u h → u strongly in L (Ω) d .We note that by the Banach–Alaoglu theorem, the closed unit ball of the dualspace of a normed vector space, (for example L (Ω)), is compact in the weak-*topology. Hence we can also find a subsequence such that ρ h ∗ (cid:42) ˆ ρ ∈ C γ ∩ { η : (cid:107) ρ − η (cid:107) L ∞ (Ω) ≤ r/ } weakly-* in L ∞ (Ω). By the uniqueness of the weak limit,we can identify ˆ ρ = ρ a.e. in Ω and, thus, we deduce that ρ h ∗ (cid:42) ρ weakly-* in L ∞ (Ω). Consequently, ρ h (cid:42) ρ weakly in L s (Ω) for all s ∈ [1 , ∞ ). Corollary 1.
Suppose that the conditions of Theorem 4 hold. Fix any isolatedlocal minimizer of (BP) and let Ω b be any measurable subset of Ω of positivemeasure 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 convergestrongly in L s (Ω b ) to the isolated local minimizer of (BP), where s ∈ [1 , ∞ ) .Proof. We have shown that there exists a sequence of finite element minimizers( u h , ρ h ) of (BP h ) that converge to the isolated local minimizer ( u , ρ ). In partic-ular ρ h ∗ (cid:42) ρ weakly-* in L ∞ (Ω) and ρ h (cid:42) ρ weakly in L s (Ω) for all s ∈ [1 , ∞ ).The result is then deduced by following the proof of Corollary 3.2 in [29]. Proposition 5.
Suppose that the conditions stated in Theorem 4 hold. Fix anisolated minimizer ( u , ρ ) of (BP). Then there exists a subsequence of minimiz-ers, ( ρ h ) , of (BP h ) such that ρ h → ρ strongly in L s (Ω) , s ∈ [1 , ∞ ) . (4.17) Proof.
We note that C γ,h ∩ B r/ ,L (Ω) ( ρ ) is a convex set, and hence for any η h ∈ C γ,h ∩ B r/ ,L (Ω) ( ρ ), t ∈ [0 , ρ h + t ( η h − ρ h ) ∈ C γ,h ∩ B r/ ,L (Ω) ( ρ ).Since ( u h , ρ h ) is a minimizer of (BP h ), by the arguments used in Proposition 2we deduce that (cid:90) Ω α (cid:48) ( ρ h ) | u h | ( η h − ρ h )d x ≥ η h ∈ C γ,h ∩ B r/ ,L (Ω) ( ρ ) . (4.18)Hence, (FOC3) and (4.18) imply that for all η ∈ C γ and η h ∈ C γ,h ∩ B r/ ,L (Ω) ( ρ )we have that (cid:90) Ω α (cid:48) ( ρ ) | u | ρ d x ≤ (cid:90) Ω α (cid:48) ( ρ ) | u | η d x, (4.19) (cid:90) Ω α (cid:48) ( ρ h ) | u h | ρ h d x ≤ (cid:90) Ω α (cid:48) ( ρ h ) | u h | η h d x. (4.20)By subtracting (cid:82) Ω α (cid:48) ( ρ ) | u | ρ h d x from (4.19) and (cid:82) Ω α (cid:48) ( ρ h ) | u h | ρ d x from (4.20),we see that (cid:90) Ω α (cid:48) ( ρ ) | u | ( ρ − ρ h ) d x ≤ (cid:90) Ω α (cid:48) ( ρ ) | u | ( η − ρ h )d x, (4.21) (cid:90) Ω α (cid:48) ( ρ h ) | u h | ( ρ h − ρ ) d x ≤ (cid:90) Ω α (cid:48) ( ρ h ) | u h | ( η h − ρ ) d x. (4.22)17umming (4.21) and (4.22) and rearranging the left-hand side, we see that (cid:90) Ω ( α (cid:48) ( ρ ) − α (cid:48) ( ρ h )) | u | ( ρ − ρ h )d x + (cid:90) Ω α (cid:48) ( ρ h )( | u | − | u h | )( ρ − ρ h )d x ≤ (cid:90) Ω α (cid:48) ( ρ ) | u | ( η − ρ h )d x + (cid:90) Ω α (cid:48) ( ρ h ) | u h | ( η h − ρ )d x. (4.23)By fixing η = ρ h ∈ C γ and subtracting the second term on the left-hand side of(4.23) from both sides we deduce that (cid:90) Ω ( α (cid:48) ( ρ ) − α (cid:48) ( ρ h )) | u | ( ρ − ρ h )d x ≤ (cid:90) Ω α (cid:48) ( ρ h ) | u h | ( η h − ρ )d x + (cid:90) Ω α (cid:48) ( ρ h )( | u h | − | u | )( ρ − ρ h )d x. (4.24)By an application of the mean value theorem, we note that there exists a c ∈ (0 ,
1) such that (cid:90) Ω ( α (cid:48) ( ρ ) − α (cid:48) ( ρ h )) | u | ( ρ − ρ h )d x = (cid:90) Ω α (cid:48)(cid:48) ( ρ h + c ( ρ − ρ h )) | u | ( ρ − ρ h ) d x. (4.25)By (A5) and the definition of U θ we bound (4.25) from below: (cid:90) Ω α (cid:48)(cid:48) ( ρ h + c ( ρ − ρ h )) | u | ( ρ − ρ h ) d x ≥ (cid:90) U θ α (cid:48)(cid:48) ( ρ h + c ( ρ − ρ h )) | u | ( ρ − ρ h ) d x ≥ α (cid:48)(cid:48) min θ (cid:107) ρ − ρ h (cid:107) L ( U θ ) . (4.26)Now we bound the right-hand side of (4.24) as follows, (cid:90) Ω α (cid:48) ( ρ h ) | u h | ( η h − ρ )d x + (cid:90) Ω α (cid:48) ( ρ h )( | u h | − | u | )( ρ − ρ h )d x ≤ α (cid:48) max ( (cid:107) u (cid:107) L (Ω) + (cid:107) u − u h (cid:107) L (Ω) ) (cid:107) ρ − η h (cid:107) L (Ω) + α (cid:48) max (cid:107) ρ − ρ h (cid:107) L q (Ω) (cid:107) u + u h (cid:107) L q (cid:48) (Ω) (cid:107) u − u h (cid:107) L (Ω) , (4.27)where 2 < q (cid:48) < ∞ in two dimensions, 2 < q (cid:48) ≤ q = 2 q (cid:48) / ( q (cid:48) − (cid:107) u + u h (cid:107) L q (cid:48) (Ω) ≤ (cid:107) u (cid:107) L q (cid:48) (Ω) + (cid:107) u h (cid:107) L q (cid:48) (Ω) ≤ (cid:107) u (cid:107) H (Ω) + (cid:107) u h (cid:107) H (Ω) ≤ ˆ C < ∞ , (4.28)where the second inequality holds thanks to the Sobolev embedding theorem.Combining (4.24)–(4.28) we see that (cid:107) ρ − ρ h (cid:107) L ( U θ ) ≤ C (cid:0) (cid:107) ρ − η h (cid:107) L (Ω) + (cid:107) ρ − ρ h (cid:107) L q (Ω) (cid:107) u − u h (cid:107) L (Ω) (cid:1) , (4.29)18here C = C ( α (cid:48) max , α (cid:48)(cid:48) min , θ, (cid:107) u (cid:107) L (Ω) , ˆ C ). By assumption (F3), there exists asequence of finite element functions ˜ ρ h ∈ C γ,h such that ˜ ρ h → ρ strongly in L (Ω). Thanks to the strong convergence, we note that for sufficiently small h , ˜ ρ h ∈ C γ,h ∩ B r/ ,L (Ω) ( ρ ). Hence we can fix η h = ˜ ρ h . By Proposition 4, weknow that u h → u strongly in L (Ω) d and since ρ ∈ C γ , ρ h ∈ C γ,h ⊂ C γ , then (cid:107) ρ − ρ h (cid:107) L q (Ω) ≤ | Ω | /q (cid:107) ρ − ρ h (cid:107) L ∞ (Ω) ≤ | Ω | /q . Therefore, the right-hand sideof (4.29) tends to zero as h →
0. Hence, we deduce that ρ h → ρ strongly in L ( U θ ) , θ > . (4.30)Now we note that (cid:107) ρ − ρ h (cid:107) L (Ω) = (cid:107) ρ − ρ h (cid:107) L ( U θ ) + (cid:107) ρ − ρ h (cid:107) L ( U \ U θ ) + (cid:107) ρ − ρ h (cid:107) L (Ω \ U ) . (4.31)If U \ U θ or Ω \ U are empty, we neglect the corresponding term in (4.31) withno loss of generality. Suppose Ω \ U is non-empty. By definition of U , u = a.e. in Ω \ U . By Proposition 3, this implies that ρ = 0 a.e. in Ω \ U . Therefore,Corollary 1 implies that ρ h → ρ strongly in L (Ω \ U ) . (4.32)Suppose U \ U θ is non-empty. Since, ρ, ρ h ∈ C γ we see that (cid:107) ρ − ρ h (cid:107) L ( U \ U θ ) ≤ | U \ U θ | / → θ → . (4.33)Therefore, by first taking the limit as h → θ →
0, (4.30)–(4.33) imply that ρ h → ρ strongly in L (Ω).Since (cid:107) ρ − ρ h (cid:107) L (Ω) ≤ | Ω | / (cid:107) ρ − ρ h (cid:107) L (Ω) , we see that ρ h → ρ strongly in L (Ω). Hence, for any s ∈ [1 , ∞ ), (cid:90) Ω | ρ − ρ h | s d x = (cid:90) Ω | ρ − ρ h | s − | ρ − ρ h | d x ≤ s − (cid:107) ρ − ρ h (cid:107) L (Ω) , (4.34)which implies that ρ h → ρ strongly in L s (Ω). Proposition 6.
Suppose that the conditions stated in Theorem 4 hold. Fix anisolated minimizer ( u , ρ ) of (BP). Then, there exists a subsequence of minimiz-ers, ( u h ) , of (BP h ) such that u h → u strongly in H (Ω) d . (4.35) Proof.
We note that the set W h := V g h ,h ∩ B r/ ,H (Ω) ( u ) is convex. Hence itcan be shown that minimizers of (BP h ) satisfy the variational inequality a ρ h ( u h , w h − u h ) − l f ( w h − u h ) ≥ w h ∈ W h . (4.36)From Proposition 2 we deduce that, for all w h ∈ W h , a ρ ( u , w h − u h ) + b ( w h − u h , p ) = l f ( w h − u h ) . a ρ h ( u h , u h − w h ) ≤ a ρ ( u , u h − w h ) + b ( u h − w h , p ) . (4.37)By subtracting a ρ h ( w h , u h − w h ) from both sides we see that a ρ h ( u h − w h , u h − w h ) ≤ a ρ ( u , u h − w h ) − a ρ h ( w h , u h − w h ) + b ( u h − w h , p − q h ) . We note that a ρ h is coercive and bounded with constants, c a = ν/ ( c p + 1) , C a = max { ν, α } , and b is bounded with constant C b . Hence, (cid:107) u h − w h (cid:107) H (Ω) ≤ c a a ρ h ( u h − w h , u h − w h ) ≤ c a ( a ρ ( u , u h − w h ) − a ρ h ( w h , u h − w h ) + b ( u h − w h , p − q h ))= 1 c a (cid:18)(cid:90) Ω α ( ρ h )( u − w h ) · ( u h − w h ) + ( α ( ρ ) − α ( ρ h )) u · ( u h − w h )d x + (cid:90) Ω ν ∇ ( u − w h ) : ∇ ( u h − w h )d x + b ( u h − w h , p − q h ) (cid:19) ≤ c a ¯ α (cid:107) u − w h (cid:107) L (Ω) (cid:107) u h − w h (cid:107) L (Ω) + 1 c a (cid:107) ( α ( ρ ) − α ( ρ h )) u (cid:107) L (Ω) (cid:107) u h − w h (cid:107) L (Ω) + νc a | u − w h | H (Ω) | u h − w h | H (Ω) + C b c a (cid:107) u h − w h (cid:107) H (Ω) (cid:107) p − q h (cid:107) L (Ω) . Hence, (cid:107) u h − w h (cid:107) H (Ω) ≤ C (cid:0) (cid:107) u − w h (cid:107) H (Ω) + (cid:107) ( α ( ρ ) − α ( ρ h )) u (cid:107) L (Ω) + (cid:107) p − q h (cid:107) L (Ω) (cid:1) , where C = C (¯ α, ν, C a , c a , C b ) is a constant. This implies that, for all w h ∈ W h , (cid:107) u − u h (cid:107) H (Ω) ≤ C (cid:0) (cid:107) u − w h (cid:107) H (Ω) + (cid:107) ( α ( ρ ) − α ( ρ h )) u (cid:107) L (Ω) + (cid:107) p − q h (cid:107) L (Ω) (cid:1) . (4.38)where C (cid:48) = C (cid:48) (¯ α, ν, C a , c a , C b , L α ). For sufficiently small h , we note that ˜ u h ∈ W h (where ˜ u h is defined in the proof of Proposition 4) and ˜ u h → u strongly in H (Ω) d . Moreover by assumption (F3), there exists a sequence of finite elementfunctions ˜ p h ∈ M h that converges to p strongly in L (Ω). Suppose w h = ˜ u h and q h = ˜ p h . From Proposition 5, we know that there exists a subsequence (notindicated) such that ρ h → ρ strongly in L (Ω). We now observe that (cid:107) ( α ( ρ ) − α ( ρ h )) u (cid:107) L (Ω) ≤ L α (cid:107) ρ − ρ h (cid:107) L (Ω) (cid:107) u (cid:107) L (Ω) (4.39)where L α is the Lipschitz constant for α . Hence by taking the limit as h → u h → u strongly in H (Ω) d .20n the following proposition, we show that (up to a subsequence) mini-mizers of (BP h ) also satisfy the first-order optimality conditions that are thefinite-dimensional analogue of the first-order optimality conditions associatedwith (BP). This allows us to consider the finite-dimensional optimization prob-lem over the whole set V g h ,h × C γ,h , rather than the restricted set ( V g h ,h ∩ B r/ ,H (Ω) ( u )) × ( C γ,h ∩ B r/ ,L (Ω) ( ρ )). Proposition 7.
Suppose that the conditions stated in Theorem 4 hold. Then,there exists an ¯ h > such that for all h < ¯ h , there exists a unique Lagrangemultiplier p h ∈ M h such that the functions ( u h , ρ h ) that locally minimize (BP h )satisfy the first-order optimality conditions (FOC1 h )–(FOC3 h )Proof. From Proposition 6, we know that u h → u strongly in H (Ω) d . Hence bydefinition of strong convergence, there exists an ¯ h > h ≤ ¯ h , (cid:107) u − u h (cid:107) H (Ω) ≤ r/
4. Therefore, for each v h ∈ V ,h , if | t | < r/ (4 (cid:107) v h (cid:107) H (Ω) )then u h + t v h ∈ V g h ,h ∩ B r/ ,H (Ω) ( u ). Now we can follow the reasoning of theproof of Proposition 2 (adding the subscript h where necessary) to deduce theexistence of a unique p h ∈ M h such that (FOC1 h )–(FOC3 h ) hold. Proposition 8.
Suppose that the conditions stated in Theorem 4 hold. Then,there is a subsequence of the unique p h ∈ M h defined in Proposition 7 thatconverges strongly in L (Ω) to the p ∈ L (Ω) that solves (FOC1)–(FOC3) forthe given isolated local minimizer ( u , ρ ) .Proof. The inf-sup condition (F2) for M h and X ,h implies that, for any q h ∈ M h , c b (cid:107) q h − p h (cid:107) L (Ω) ≤ sup w h ∈ X ,h \{ } b ( w h , q h − p h ) (cid:107) w h (cid:107) H (Ω) = sup w h ∈ X ,h \{ } b ( w h , p − p h ) + b ( w h , q h − p ) (cid:107) w h (cid:107) H (Ω) ≤ sup w h ∈ X ,h \{ } | b ( w h , p − p h ) | + | b ( w h , q h − p ) |(cid:107) w h (cid:107) H (Ω) = sup w h ∈ X ,h \{ } | a ρ ( u , w h ) − a ρ h ( u h , w h ) | + | b ( w h , q h − p ) |(cid:107) w h (cid:107) H (Ω) ≤ (cid:107) ( α ( ρ ) − α ( ρ h )) u (cid:107) L (Ω) + (¯ α + 1) (cid:107) u − u h (cid:107) H (Ω) + C b (cid:107) p − q h (cid:107) L (Ω) . Hence, (cid:107) p − p h (cid:107) L (Ω) ≤ C (cid:0) (cid:107) ( α ( ρ ) − α ( ρ h )) u (cid:107) L (Ω) + (cid:107) u − u h (cid:107) H (Ω) + (cid:107) p − q h (cid:107) L (Ω) (cid:1) , where C = C ( c b , C b , ¯ α, L α ). By assumption (F3), there exists a sequence offinite element functions, ˜ p h ∈ M h that satisfies ˜ p h → p strongly in L (Ω).Let q h = ˜ p h . We have already shown that u h → u strongly in H (Ω) d inProposition 6. Similarly, in the proof of Proposition 6 we also showed that21 ( α ( ρ ) − α ( ρ h )) u (cid:107) L (Ω) →
0. Hence we conclude that p h → p strongly in L (Ω). Proof of Theorem 4.
Fix an isolated minimizer ( u , ρ ) of (BP) and its uniqueassociated Lagrange multiplier p . By the results of Propositions 4, 5, 6, and 7,there exists a mesh size ¯ h such that for, h < ¯ h , there exists a sequence of finiteelement solutions ( u h , ρ h , p h ) ∈ V g h ,h × C γ,h × M h satisfying (FOC1 h )–(FOC3 h )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 6implies that u h → u strongly in H (Ω) d , and Proposition 8 implies that p h → p strongly in L (Ω).
5. Numerical results
The main goal of this section is to experimentally verify the existence ofstrongly converging sequences that were proven to exist in Section 4. In allexamples the systems were discretized with the finite element method usingFEniCS [25]. The computational domains are triangulated with simplices andwe define the mesh size, h , as the maximum diameter of all the simplices in thetriangulation. The solutions were computed using the deflated barrier method[28], an algorithm that can systematically discover multiple solutions of topologyoptimization problems, via deflation [16], by solving the discretized first-orderoptimality conditions, (FOC1 h )–(FOC3 h ). The resulting linear systems arisingin the deflated barrier method were solved by a sparse LU factorization withMUMPS [8] and PETSc [9].There are no known analytical solutions for choices of the inverse permeabil-ity, α , used in practice. Hence, errors are measured with respect to a heavily-refined finite element solution, which is constructed as follows; first the finiteelement solutions are computed on a mesh with mesh size h = 0 . h ), that attains the smallest objective func-tional value for J , within the basin of attraction of the isolated local minimizer.For sufficiently small h , a minimizer satisfying this selection mechanism mustexist. However, it is not necessarily unique and numerically enforcing such acondition can be difficult. In order to promote convergence to the minimizerof (BP h ) with the smallest value J , we interpolate the heavily-refined finite el-ement solutions onto coarser meshes as initial guesses for the deflated barriermethod. This strategy was effective in practice. The effects of the second obser-vation (P2) are harder to test. However, in Section 5.2, we attempt to minimizemesh bias by measuring errors on unstructured meshes. Code availability:
For reproducibility, an implementation of the deflatedbarrier method as well as scripts to generate the convergence plots and solutionscan be found at https://bitbucket.org/papadopoulos/deflatedbarrier/ .The version of the software used in this paper is archived on Zenodo [1].
Consider the optimization problem (BP), with a homogeneous Dirichletboundary condition on u , Ω = (0 , , volume fraction γ = 1 /
3, viscosity ν = 1and a forcing term given by f ( x, y ) = (cid:40) (10 , (cid:62) if 3 / < x < /
10 and 3 / < y < / , (0 , (cid:62) otherwise . (5.1)The inverse permeability, α , is as given in (2.1), with α = 2 . × and q = 1 / q is a penalty parameter which controls the levelof intermediate values (between zero or one) in the optimal design. Figure 1depicts the material distribution of three minimizers. One local minimizer is inthe shape of a figure eight and the two Z symmetric global minimizers are inthe shape of annuli. Since the domain is convex, g = , and f ∈ L (Ω) d , then,by the regularity results proven in the Appendix, u ∈ H (Ω) and p ∈ H (Ω).The conditions of Theorem 3 hold and, therefore, ρ ∈ H ( U θ ) for every θ >
0. Inthis particular example, the support of ρ is compactly contained in the supportof the velocity in all three solutions. Therefore, we conclude that ρ ∈ H (Ω).Consider a ( P ) × P Taylor–Hood finite element discretization for the ve-locity and the pressure, and a P piecewise constant finite element discretizationfor the material distribution. Since all three solutions are isolated local mini-mizers, by Theorem 4, there exists a sequence of finite element solutions to thediscretized first-order optimality conditions that strongly converges to the figureeight solution, and different sequences of different finite element solutions that23 igure 1: The material distribution of a local (left) and the global (middle and right) min-imizers of the discontinuous-forcing optimization problem. Black corresponds to a value of ρ = 0 and white corresponds to a value of ρ = 1, with the gray regions indicating intermediatevalues. The arrows indicate the velocity profile of the solutions. strongly converge to the two annulus solutions. Their existence is confirmed inFigure 2.Since ρ ∈ H (Ω) and we are using a P finite element discretization, ana¨ıve prediction for the convergence rate of the L -norm error of the materialdistribution 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 (Ω) d , and weare using a ( P ) finite element discretization, a prediction for the convergencerate of the velocity is O ( h r ) and O ( h t +1 ), r, t ∈ [1 , H -norm and L -norm errors of the velocity, respectively. The rates observed in the toppanels of Figure 2 match this prediction. The H -norm error is decreasingat a rate slightly faster than O ( h / ) for all three solutions and the L -normerror convergence rate is O ( h ) for the figure eight solution and O ( h / ) for theannuli solutions. We hypothesize that the upper limit convergence rate with r = t = 2 is not achieved because of the lack of regularity in the forcing term, f ∈ H s (Ω) d , s < /
2, and the relatively slower rate of the convergence of thematerial distribution. Finally, since the minimum regularity of the pressure is p ∈ H (Ω) and the discretization is P , a prediction for the convergence rate ofthe L -norm error is O ( h r ), 1 ≤ r ≤
2. Initially, the convergence rate is O ( h / )which matches our na¨ıve prediction. However, on finer meshes, the convergencerate increases. We hypothesize that this speedup is artificial and is caused bythe lack of resolution of the refined finite element solutions that are being usedas proxies for the analytical solutions in the error norm estimate. Qualitatively,it can be checked that mesh refinement in areas where the discretized materialdistribution lies between 1/10 and 9/10 is an ineffective strategy for improvingthe approximation of the analytical pressure over the whole domain. Consider the optimization problem (BP) [11, Sec. 4.5], with two prescribedflow inputs and two prescribed outputs, where Ω = (0 , / × (0 , γ = 1 / f = (0 , (cid:62) and ν = 1 and the boundary conditions on u are given by the24 − − × − × − × − × − h − − k u − u h k L ( Ω ) L (Ω)-norm error of the velocity Upper annulusFigure eightLower annulus O ( h ) O ( h ) − − × − × − × − × − h − − k u − u h k H ( Ω ) H (Ω)-norm error of the velocity Upper annulusFigure eightLower annulus O ( h ) O ( h ) − − × − × − × − × − h − × − × − × − × − k ρ − ρ h k L ( Ω ) L (Ω)-norm error of the material distribution Upper annulusFigure eightLower annulus O ( h ) − − × − × − × − × − h − − k p − p h k L ( Ω ) L (Ω)-norm error of the pressure Upper annulusFigure eightLower annulus O ( h ) O ( h ) Figure 2: The convergence of u h , ρ h , and p h in the discontinuous-forcing problem for thefigure eight and annulus solutions on structured meshes, with a P × ( P ) × P discretizationfor ( ρ h , u h , p h ). boundary data g ( x, y ) = (cid:0) − y − / , (cid:1) (cid:62) if 2 / ≤ y ≤ / , x = 0 or 3 / , (cid:0) − y − / , (cid:1) (cid:62) if 1 / ≤ y ≤ / , x = 0 or 3 / , (0 , (cid:62) elsewhere on ∂ Ω . (5.2)The function α is as given in (2.1), with α = 2 . × and q = 1 /
10. Figure3 visualizes the setup of the problem and depicts the material distribution ofthe two minimizers. One local minimizer is a straight channel solution and theglobal minimizer is in the form of a double-ended wrench.We employ a ( P ) × P Taylor–Hood finite element discretization for thevelocity and the pressure, and a P continuous piecewise linear finite elementdiscretization for the material distribution. This example satisfies all the condi-tions of Theorem 4 and, for both minimizers, we numerically verify that thereexists a sequence of finite element solutions that strongly converges to it inFigure 4.As mentioned earlier in this section, it may be possible to find a sequence ofmesh sizes, ( h i ), such that there exist two different sequences of finite elementsolutions that strongly converge to the same isolated minimizer. In Figure 5, we25 igure 3: The setup of the double-pipe problem is shown on the left. The middle and rightfigures depict the material distribution of the straight channel and double-ended wrench so-lutions, respectively. Black corresponds to a value of ρ = 0 and white corresponds to a valueof ρ = 1, with the gray regions indicating intermediate values. − h − − − k u − u h k L ( Ω ) L (Ω)-norm error of the velocity Straight channelsDouble-ended wrench − h − k u − u h k H ( Ω ) H (Ω)-norm error of the velocity Straight channelsDouble-ended wrench − h − − k ρ − ρ h k L ( Ω ) L (Ω)-norm error of the material distribution Straight channelsDouble-ended wrench − h − k p − p h k L ( Ω ) L (Ω)-norm error of the pressure Straight channelsDouble-ended wrench
Figure 4: The convergence of u h , ρ h , and p h for the double-pipe problem for both the straightchannel and double-ended wrench solutions on an unstructured mesh with a P × ( P ) × P discretization for ( ρ h , u h , p h ). depict two different straight channel finite element solutions that exist on thesame unstructured mesh where h = 0 .
04. Both solutions satisfy the discretizedfirst-order optimality conditions (FOC1 h )–(FOC3 h ) and both locally minimize J ( v h , η h ). Choosing one over the other would change the convergence patternof the strongly converging sequence. This may cause difficulty in practice, asoptimization strategies are unlikely to discover the discretized global minimumwithout additional selection mechanisms.26 igure 5: Two different straight channel finite element solutions of the double-pipe optimiza-tion problem that exist on the same unstructured mesh where h = 0 .
04. The differences canbe spotted at the midway point of the top channel.
6. Conclusions
In this work we have studied the fluid topology optimization model of Bor-rvall and Petersson [11]. In the case of a homogeneous Dirichlet boundarycondition and under a mild convexity assumption on the inverse permeabilityterm, α , we have shown that the material distribution necessarily lives in theSobolev space H inside any compact subset of the support of the velocity.Moreover, we have formally treated the nonconvexity of the optimization prob-lem (including the case of inhomogeneous Dirichlet boundary conditions) andhave shown that, given an analytical isolated local minimizer, there exists asequence of discretized solutions, satisfying the associated first-order optimalityconditions, that strongly converges to the minimizer in the appropriate spaces.We have numerically verified that these sequences exist and have discussed theobserved convergence rates. Appendix. Elliptic regularity of the generalized Stokes system withDirichlet boundary dataLemma 5.
Let the domain Ω be either a convex polygon in two dimensionsor a convex polyhedron in three dimensions and consider the triple ( u , ρ, p ) ∈ H g (Ω) d × C γ × L (Ω) that satisfies (FOC1)–(FOC3). Suppose that the forcingterm f ∈ L (Ω) d and the boundary datum g is the boundary trace of a function ˆ g ∈ H (Ω) d on the boundary ∂ Ω , where div(ˆ g ) = 0 a.e. in Ω . Then u ∈ H (Ω) d and p ∈ H (Ω) .Proof. The idea of the proof is to reduce the system (FOC1)–(FOC2) to thestandard Stokes system with a homogeneous Dirichlet boundary condition andthen invoke the regularity results of Kellogg and Osborn [21] in two dimensionsand Kozlov et al. [22] in three dimensions.Let w := u − ˆ g . Since the trace operator is a linear operator, we see that w | ∂ Ω = ( u − ˆ g ) | ∂ Ω = g − g = . Since u ∈ H (Ω) d and ˆ g ∈ H (Ω) d , then w ∈ H (Ω) d . 27y substituting w into (FOC1)–(FOC2), we see that (FOC1)–(FOC2) isequivalent to finding ( w , p ) ∈ H (Ω) d × L (Ω) that satisfies for all ( v , q ) ∈ H (Ω) d × L (Ω): (cid:90) Ω ∇ w : ∇ v − p div( v ) d x = (cid:90) Ω ( f − α ( ρ )( w + ˆ g )) · v − ∇ ˆ g : ∇ v d x, (A.1) (cid:90) Ω q div( w + ˆ g ) d x = 0 . (A.2)Define ˆ f as ˆ f := f − α ( ρ )( w + ˆ g ) + ∆ ˆ g . Since f ∈ L (Ω) d , α ( ρ ) ∈ L ∞ (Ω), w ∈ H (Ω) d , and ˆ g ∈ H (Ω) d , then ˆ f ∈ L (Ω) d . By assumption div( ˆ g ) = 0a.e. in Ω and, hence, by an application of integration by parts on the final termon the right-hand side of (A.1), we see that (A.1)–(A.2) is equivalent to finding( w , p ) ∈ H (Ω) d × L (Ω) that satisfies for all ( v , q ) ∈ H (Ω) d × L (Ω): (cid:90) Ω ∇ w : ∇ v − p div( v ) d x = (cid:90) Ω ˆ f · v d x, (A.3) (cid:90) Ω q div( w ) d x = 0 . (A.4)We note that (A.3)–(A.4) is the standard Stokes system with a homogeneousDirichlet boundary condition and forcing term ˆ f ∈ L (Ω) d . Therefore, by theelliptic regularity of the Stokes system [21, 22], w ∈ H (Ω) d and p ∈ H (Ω).Since u = w + ˆ g and ˆ g ∈ H (Ω) d , we conclude that u ∈ H (Ω) d . References [1]
Software used in ‘Numerical analysis of a topology optimization problem forStokes flow’ , 2021, https://doi.org/10.5281/zenodo.4514054 .[2]
N. Aage, T. H. Poulsen, A. Gersborg-Hansen, and O. Sig-mund , Topology optimization of large scale Stokes flow problems , Struc-tural and Multidisciplinary Optimization, 35 (2008), pp. 175–180, https://doi.org/10.1007/s00158-007-0128-0 .[3]
L. Adam, M. Hinterm¨uller, D. Peschka, and T. M. Surowiec , Optimization of a multiphysics problem in semiconductor laser design ,SIAM Journal on Applied Mathematics, 79 (2019), pp. 257–283, https://doi.org/10.1137/18M1179183 .[4]
J. Alexandersen and C. S. Andreasen , A review of topology optimi-sation for fluid-based problems , Fluids, 5 (2020), p. 29, https://doi.org/10.3390/fluids5010029 .[5]
G. Allaire , Shape optimization by the homogenization method , vol. 146,Springer Science & Business Media, 2012, https://doi.org/10.1007/978-1-4684-9286-6 . 286]
D. H. Alonso, L. F. N. de S´a, J. S. R. Saenz, and E. C. N.Silva , Topology optimization applied to the design of 2D swirl flow devices ,Structural and Multidisciplinary Optimization, 58 (2018), pp. 2341–2364, https://doi.org/10.1007/s00158-018-2078-0 .[7]
D. H. Alonso, J. S. R. Saenz, and E. C. N. Silva , Non-newtonianlaminar 2D swirl flow design by the topology optimization method , Struc-tural and Multidisciplinary Optimization, (2020), pp. 1–23, https://doi.org/10.1007/s00158-020-02499-2 .[8]
P. R. Amestoy, I. S. Duff, J.-Y. L’Excellent, and J. Koster , Afully asynchronous multifrontal solver using distributed dynamic scheduling ,SIAM Journal on Matrix Analysis and Applications, (2001).[9]
S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune,K. Buschelman, L. Dalcin, A. Dener, V. Eijkout, W. Gropp,R. Tran Mills, T. Munson, K. Rupp, P. Sana, B. Smith,S. Zampini, H. Zhang, and H. Zhang , PETSc Users Manual , Tech.Report ANL-95/11 - Revision 3.11, Argonne National Laboratory, 2019, .[10]
M. P. Bendsøe and O. Sigmund , Topology Optimization , SpringerBerlin Heidelberg, Berlin, Heidelberg, 2004, https://doi.org/10.1007/978-3-662-05086-6 .[11]
T. Borrvall and J. Petersson , Topology optimization of fluids inStokes flow , International Journal for Numerical Methods in Fluids, 41(2003), pp. 77–107, https://doi.org/10.1002/fld.426 .[12]
Y. Deng, Y. Wu, and Z. Liu , Topology Optimization Theory for LaminarFlow: Applications in Inverse Design of Microfluidics , Springer Singapore,Singapore, 2018, https://doi.org/10.1007/978-981-10-4687-2 .[13]
L. C. Evans , Partial Differential Equations , American Mathematical So-ciety, 2 ed., 2010.[14]
A. Evgrafov , Topology optimization of slightly compressible fluids ,ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift f¨urAngewandte Mathematik und Mechanik: Applied Mathematics andMechanics, 86 (2006), pp. 46–62, https://doi.org/10.1002/zamm.200410223 .[15]
A. Evgrafov , State space Newton’s method for topology optimization ,Computer Methods in Applied Mechanics and Engineering, 278 (2014),pp. 272–290, https://doi.org/10.1016/j.cma.2014.06.005 .[16]
P. E. Farrell, ´A. Birkisson, and S. W. Funke , Deflation tech-niques for finding distinct solutions of nonlinear partial differential equa-tions , SIAM Journal on Scientific Computing, 37 (2015), pp. A2026–A2045, https://doi.org/10.1137/140984798 .2917]
A. Gersborg-Hansen, O. Sigmund, and R. B. Haber , Topol-ogy optimization of channel flow problems , Structural and Multidisci-plinary Optimization, 30 (2005), pp. 181–192, https://doi.org/10.1007/s00158-004-0508-7 .[18]
V. Girault and P.-A. Raviart , Finite element methods for Navier–Stokes equations: theory and algorithms , vol. 5, Springer-Verlag Berlin Hei-delberg, 1986, https://doi.org/10.1007/978-3-642-61623-5 .[19]
J. K. Guest and J. H. Pr´evost , Topology optimization of creep-ing fluid flows using a Darcy–Stokes finite element , International Journalfor Numerical Methods in Engineering, 66 (2006), pp. 461–484, https://doi.org/10.1002/nme.1560 .[20]
I. G. Jang and I. Y. Kim , Computational study of Wolff ’s law withtrabecular architecture in the human proximal femur using topology op-timization , Journal of Biomechanics, 41 (2008), pp. 2353–2361, https://doi.org/10.1016/j.jbiomech.2008.05.037 .[21]
R. B. Kellogg and J. E. Osborn , A regularity result for the Stokesproblem in a convex polygon , Journal of Functional Analysis, 21 (1976),pp. 397–431, https://doi.org/10.1016/0022-1236(76)90035-5 .[22]
V. A. Kozlov, V. G. Maz’ya, and C. Schwab , On singularities of so-lutions to the Dirichlet problem of hydrodynamics near the vertex of a cone ,Journal f¨ur die reine und angewandte Mathematik, 456 (1994), pp. 65–97, https://doi.org/10.1515/crll.1994.456.65 .[23]
S. Kreissl, G. Pingen, and K. Maute , Topology optimization for un-steady flow , International Journal for Numerical Methods in Engineering,87 (2011), pp. 1229–1253, https://doi.org/10.1002/nme.3151 .[24]
J. Liu, A. T. Gaynor, S. Chen, Z. Kang, K. Suresh, A. Takezawa,L. Li, J. Kato, J. Tang, C. C. L. Wang, L. Cheng, X. Liang, andA. C. To , Current and future trends in topology optimization for additivemanufacturing , Structural and Multidisciplinary Optimization, 57 (2018),pp. 2457–2483, https://doi.org/10.1007/s00158-018-1994-3 .[25]
A. Logg, K.-A. Mardal, and G. Wells , Automated solution of differ-ential equations by the finite element method: The FEniCS book , SpringerScience & Business Media, 84 (2012).[26]
L. H. Olesen, F. Okkels, and H. Bruus , A high-level programming-language implementation of topology optimization applied to steady-stateNavier–Stokes flow , International Journal for Numerical Methods in Engi-neering, 65 (2006), pp. 975–1001, https://doi.org/10.1002/nme.1468 .[27]
W. S. O˙za´nski , The Lagrange multiplier and the stationary Stokes equa-tions , Journal of Applied Analysis, 23 (2017), pp. 137–140, https://doi.org/10.1515/jaa-2017-0017 . 3028]
I. P. A. Papadopoulos, P. E. Farrell, and T. M. Surowiec , Com-puting multiple solutions of topology optimization problems , SIAM Journalon Scientific Computing, (Accepted for publication), (2021).[29]
J. Petersson , A finite element analysis of optimal variable thicknesssheets , SIAM Journal on Numerical Analysis, 36 (1999), pp. 1759–1778, https://doi.org/10.1137/S0036142996313968 .[30]
L. F. N. S´a, J. S. Romero, O. Horikawa, and E. C. N. Silva , Topology optimization applied to the development of small scale pump ,Structural and Multidisciplinary Optimization, 57 (2018), pp. 2045–2059, https://doi.org/10.1007/s00158-018-1966-7https://doi.org/10.1007/s00158-018-1966-7