{'fraction_non_alphanumeric': 0.03054772484750028, 'fraction_numerical': 0.01278927166451654, 'mean_word_length': 4.709557151950532, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
CutFEM topology optimization of 3D laminar incompressible flow problems

This paper studies the characteristics and applicability of the CutFEM approach as the core of a robust topology optimization framework for 3D laminar incompressible flow and species transport problems at low Reynolds number (Re<200). CutFEM is a methodology for discretizing partial differential equations on complex geometries by immersed boundary techniques. In this study, the geometry of the fluid domain is described by an explicit level set method, where the parameters of a level set function are defined as functions of the optimization variables. The fluid behavior is modeled by the incompressible Navier-Stokes equations. Species transport is modeled by an advection-diffusion equation. The governing equations are discretized in space by a generalized extended finite element method. Face-oriented ghost-penalty terms are added for stability reasons and to improve the conditioning of the system. The boundary conditions are enforced weakly via Nit\-sc\-he's method. The emergence of isolated volumes of fluid surrounded by solid during the optimization process leads to a singular analysis problem. An auxiliary indicator field is modeled to identify these volumes and to impose a constraint on the average pressure. Numerical results for 3D, steady-state and transient problems demonstrate that the CutFEM analyses are sufficiently accurate, and the optimized designs agree well with results from prior studies solved in 2D or by density approaches.


Introduction
The performance of a broad range of engineering applications is dependent on the characteristics of internal and external flows. Such applications include flows for momentum, energy and species transport. The performance of these systems can be improved by finding the optimal geometry of the fluid-solid interface. The increase of computational power and the development of improved numerical schemes allow for the systematic design of flow problems via mathematical optimization methods. In this paper, we study the influence of conceptual changes to the design geometry of fluid flow problems within the context of topology optimization.
Initial ventures into fluid optimization focused on shape optimization [2,3]. Topology optimization of fluid problems was pioneered by [4], who adopted the concept of density methods to Stokes flows. Density methods were originally developed for problems in solid mechanics [5,6], and describe the geometry of a body by the density distribution of a fictitious two-phase material. The density, or volume fraction, varies continuously to describe a smooth transition between two distinct phases, e.g. fluid and solid. Fluid topology optimization was extended to a Darcy-Stokes flow model by [7], and to the Navier-Stokes equations by [8] and [9], among others.
Common to most density methods in fluid topology optimization is the use of a Brinkman model [10] to weakly enforce no-slip conditions at the fluid-solid interface. This approach may suffer from errors if the optimization process converges to a "0-1" material distribution and the fluid-solid interface is not aligned with the mesh. In this case, the interface is approximated by a stair-step pattern that hinders the precise determination of the optimum geometry and may affect the accuracy of the flow solution. This issue is more pronounced for coarse meshes and high Reynolds numbers [11], and can be mitigated by mesh refinement [12,13,14]. However, the required level of mesh refinement may lead to high computational costs. Adaptive mesh refinement may affect the convergence of the optimization process if a gradient-based optimization algorithm is applied [15].
The shortcomings of density methods have promoted the development of Level Set Methods (LSMs) for topology optimization, c.f. [16,17]. External and internal phase boundaries are described implicitly by the zero level set isosurfaces of a Level Set Function (LSF), φ (x), where x is the position vector [18,19,20]. LSMs are well-suited for topology optimization because smooth differentiable changes in the LSFs lead to changes in the topology of the body, such as holes merging or splitting [21,22,23]; however, these changes may lead to discontinuities in the physical response [24].
The LSF is typically discretized on a fixed background mesh and updated in the optimization process by solving the Hamilton-Jacobi equations [19,20]. An alternate approach that is utilized in this work is to define the parameters of the discretized LSF as explicit functions of the optimization variables. The resulting parameter optimization problem is solved by standard nonlinear programming (NLP) methods, which provide flexibility with respect to the class of optimization problems one can model [17]. In the context of flow topology optimization, NLP schemes allow, for example, including additional optimization variables in addition to the parameters of the discretized LSF, such as the position and shape of the inlets and outlets. For a detailed discussion of the LSM, the reader is referred to the comprehensive review by [17] and [25].
Several methods exist to represent the geometry and the material distribution described by the LSF in the mechanical model. The Ersatz material method [19,26] interpolates the physical properties of a fictitious material by using either an element-wise constant volume fraction or the LSF value at a point [27]. In a fluid-solid problem, the Ersatz material approach defines typically the local permeability of the Brinkman model. While the Ersatz material approach eases the computational complexity, the method faces the same issues as density methods in regards to enforcing boundary conditions.
In this work, we utilize the eXtended Finite Element Method (XFEM) to describe the material distribution in the mechanical model. The XFEM is an immersed boundary technique that does not require a mesh that conforms to the phase boundaries. The method was built upon the concept of partition of unity developed by [28], and it was originally used to model crack propagation [29]. The XFEM augments the standard finite element interpolation space with additional enrichment functions to capture discontinuities in either the state variables or their spatial gradients within an element. The XFEM decomposes the elements cut by the zero level set isosurfaces into subdomains and interfaces that it uses to integrate the weak form of the governing equations. This approach avoids the need for material interpolation schemes used in density methods because each subdomain has a distinct phase. Boundary conditions on the interface are imposed weakly via penalty methods [30], stabilized Lagrange multipliers [31], or via Nitsche's method [32]. For a general overview of the XFEM, the reader is referred to [33]. In the context of topology optimization of fluid flow problems, the enforcement of no-slip boundary conditions along the phase boundaries via the XFEM and a stabilized Lagrange multiplier method was adopted by [34] for a Navier-Stokes fluid model; a hydrodynamic Boltzmann model was employed in combination with a level set-based interface representation for generalized topology optimization of fluids by [35].
One challenge of the XFEM is that an ill-conditioned system of equations results when the ratio of the phases volumes in an intersected element is very small or very large. Such interface configurations are often unavoidable when using fixed meshes in topology optimization. In general, the ill-conditioning impedes the convergence of solvers for nonlinear problems and reduces the performance of iterative linear solvers. Several approaches have been proposed to avoid this ill-conditioning issue, such as the geometric preconditioner of [36], the Jacobi preconditioner of [37], the preconditioners of [38] and [39] based on a Cholesky decomposition, and face-oriented ghost-penalty methods [40]. Face-oriented ghost-penalty methods have been studied in the context of fluid flow problems, where discontinuities in the spatial gradients of the velocities and the pressure across the common facets of intersected elements are penalized [41,42]. Here, we explore the characteristics and the performance of the ghost-penalty methods in the context of flow topology optimization.
A second challenge of the XFEM in the context of flow topology optimization is that isolated volumes of fluid surrounded by solid may emerge during the optimization process, as shown in Figure 1. These "puddles" produce a singular analysis problem because the absolute value of the pressure field is not governed. A similar issue was encountered by [43] in the context of linear elasticity, where isolated regions of solid material could undergo rigid-body motion. The issue was avoided by placing the entire domain on a system of soft springs. A similar approach is adopted and studied here for incompressible flow, where we augment our fluid model with a penalty formulation to enforce a constraint on the average pressure. However, unlike [43], where the penalty was applied to the entire domain, we model an auxiliary indicator field to detect these isolated volumes of fluid and only apply the penalty there.
A third challenge is the computational cost of performing topology optimization on 3D flows. Recent developments into large-scale studies using high-resolution density methods can be found in [44], [45], and [46]. As proof of concept, here we focus on laminar flows to study the characteristics of our proposed approach, avoiding thus the additional costs associated with studying large Reynolds number flows on highly refined meshes.
The combination of the LSM, the XFEM, and face-oriented ghost-penalty methods is referred to as CutFEM [1]. CutFEM can handle problems with a complex geometry description through the LSM. The boundary and interface conditions are part of the discrete formulation through the XFEM. Face-oriented ghost-penalty methods are added for stability. In this paper, we introduce and study the characteristics of the CutFEM approach as the core of a robust and broadly applicable topology optimization framework. We focus on 3D laminar flow problems modeled by the incompressible Navier-Stokes equations, as well as species transport problems. We consider two-phase problems, where one phase represents fluid, and the other represents solid. We expand the CutFEM by using a generalized Heaviside enrichment strategy for the XFEM, and adding terms to control the pressure in isolated volumes of fluid surrounded by solid. This work builds upon previous studies by [11] and [35], but it enhances the flexibility, accuracy, and robustness of the XFEM flow analysis by weakly enforcing boundary conditions via a Nitsche formulation that accounts for viscousand convection-dominated flow regimes, and applying face-oriented ghost-penalty methods to the velocity and pressure fields. We will show through numerical examples that these enhancements allow solving 3D steady-state and transient fluid topology optimization problems at low Reynolds number (Re < 200), such that relevant physical phenomena are sufficiently resolved and the boundary conditions are enforced well on the boundary geometries generated in the course of the optimization process.
The paper is structured as follows: Section 2 describes the formulation of our optimization problems. Section 3 details the description of the geometry using the design variables and the LSM. In Section 4, we present the governing equations and the discretization by the XFEM. Numerical examples are studied in Section 5.
The conclusions drawn from this study are presented in Section 6.

Optimization Model
The optimization problems considered here are formulated with respect to an objective and one or more inequality constraints for some desired functionality. The objective and constraints are defined in terms of design criteria, such as drag, power dissipation, fluid volume, etc. These design criteria can depend explicitly on the state and optimization variables, e.g. drag, or only on the optimization variables, e.g. fluid volume.
As the state variables may depend on time, the design criteria are either integrated over a specific period or at a given instance in time. The optimization problems of interest are formulated as follows: where s is the vector of optimization variables, of size N s , and u (t) is the vector of time-dependent discrete state variables, of size N u . The objective function Z is the integral of the time-dependent function z over the interval [t 1 , t 2 ]. The function g i is the i-th optimization constraint, and N g is the number of inequality constraints. In general, the constraints g i can be formulated analogously to the objective. However, in this study we consider state-dependent constraints that depend on the flow solution only at time t 2 . The optimization variables s i are bounded by lower and upper limits, s L i and s U i , respectively. The state variables satisfy the residual of the governing equations, R (s, u (t)) = 0. Note that we consider the state variables dependent on the optimization variables and, thus, solve the optimization problem in a reduced space over just s.

Geometry
In our CutFEM approach, the parameters of a discretized LSF are defined as explicit functions of the optimization variables. In the following subsections, we describe the concepts behind the LSM, and the parametrization of the LSF with respect to the optimization variables.

Level Set Method
The LSM describes the geometry of a body immersed in a domain by the zero isosurfaces of the LSF, φ (x). Considering a two-phase, fluid-solid problem, such as the one displayed in Figure 1, the material distribution is defined as follows: where Ω s is occupied by the solid phase, Ω f is occupied by the fluid phase, and Γ fs is the fluid-solid interface. The external boundaries are defined by the outer surface of the design domain, Ω d = Ω s ∪Ω f , and are denoted by Γ m , where the superscript m denotes the material phase; see Figure 1.

Parametrization of the Level Set Function
We parametrize the LSF to describe a combination of geometric primitives and free-form shapes. In this study, geometric primitives, such as cylinders, are used to represent outflow ports; see Subsection 3.2.2. The geometry of the fluid-solid interface within the design domain is not restricted to a particular shape and topology. Free-form shapes are parametrized by a finite element approach; see Subsection 3.2.1.
Using either parametrization scheme, the LSF is mapped onto the XFEM mesh. The level set values at the nodes are evaluated using the parametrized LSF. Within an XFEM element, the LSF is interpolated by standard trilinear finite element shape functions. Note that this approach restricts the number an element edge can be intersected by the fluid-solid interface to at most one. Therefore, the geometry resolution is limited by the size of an XFEM element, and convergence issues have been observed if sub-element-size features tend to emerge in the optimization process [24].
To prevent the formation of sub-element-size features, [47] proposes to add a constraint on the spatial gradient of the LSF in combination with restricting the nodal level set values to ±h/2, where h is the element size. As we did not observe issues with the formation of sub-element-size features in the numerical results presented in Section 5, we did not include the feature size control of [47] in the formulations of our optimization problems, although our approach does allow adding additional constraints.

Topology Optimization
To allow for the emergence of free-form geometries within the design domain, the LSF is parametrized by local shape functions defined on a finite element mesh. In general, this mesh can differ from the XFEM mesh. To ease the computational complexity, here we use the XFEM mesh to discretize the LSF to describe a free-form geometry.
One optimization variable, s i , for i = 1 . . . N n , is assigned to each node of the XFEM mesh, where N n is the number of nodes. The LSF value of the i-th node, φ i , is defined via the following linear filter: where r φ is the filter radius. The linear filter was used previously in the studies of [34,35] to widen the zone of influence of the optimization variables, and to improve the convergence rate. Furthermore, the filter may promote (but does not guarantee) smooth shapes of the phase boundaries; however, in contrast to density or sensitivity filters used in density methods, the filter above does not guarantee control of the minimum feature size [43].

Ports Geometry
In contrast to Hamilton-Jacobi-based updated schemes, formulating and solving the topology optimization problems via NLPs provide the flexibility to introduce additional optimization variables that do not stem from a finite element discretization of the LSF. Here we take advantage of this feature and define the geometry of cylindrical outflow ports in terms of optimization variables. The LSF of the j-th cylinder, φ c,j (x), is defined in its local coordinate system,x, as follows: wherex 1 andx 2 are the in-plane coordinates of the cylinder,x 3 points along the cylinder axis, and r c,j is the radius of the cylinder. The transformation from the global to the local coordinate system is defined by the rotation matrix T and the location of the cylinder center x c . The radius, the orientation of the local coordinate system, and the position of the cylinder center can be defined as functions of the optimization variables.
The nodal level set value, φ i , in the XFEM mesh is defined by the minimum value among all ports, ignoring the LSF associated with the free-form geometry, defined by (3). The minimum is approximated by the Kreisselmeier-Steinhauser function [48] to ensure differentiability of the formulation with respect to the cylinder parameters.

Analysis
The main challenge in optimizing the topology of fluid flow and convective species transport problems is the modeling and numerical prediction of the flow and species fields. This section introduces the weak form of the governing equations and outlines the spatial discretization schemes. A schematic of the problem setup is shown in Figure 1. The external surfaces and the fluid-solid interface, upon which Dirichlet boundary conditions are applied, are denoted by Γ m D and Γ fs D , respectively. The surfaces Γ m N and Γ fs N are defined analogously for the Neumann boundary conditions.

Governing Equations
In this study, we model the flow by the incompressible Navier-Stokes equations, which describe the transport of momentum and the conservation of mass. Species transport is described by coupling an advectiondiffusion equation to our flow model. An indicator field is introduced to identify isolated volumes of fluid where a constraint on the average fluid pressure is enforced. The governing equations in the fluid phase are summarized subsequently. Note that the solid phase is considered "void" in this study; that is, we do not model any physical phenomena in the solid phase.

Incompressible Navier-Stokes Equations
The residual of the weak form of the incompressible Navier-Stokes equations, denoted as r u,p , is decomposed into volumetric and surface contributions: where r Ω u,p and rΩ u,p are the residuals of the volumetric contributions, non-stabilized and stabilized, respectively. The term r Ω p,ψ is used to enforce a constraint on the average pressure in isolated volumes of fluid surrounded by solid. The terms r D u,p and r fs u,p enforce the Dirichlet boundary conditions on the external surfaces, and at the fluid-solid interface, respectively. The residual of the Neumann conditions on the external boundaries is r N u,p . The stabilization term, rΩ u,p , depends on the discretization scheme, and the ghost-penalty term, r GP u,p , depends on the face-oriented ghost-penalty formulation; both are defined in Sections 4.4 and 4.5, respectively.
The non-stabilized volumetric contribution is formulated as: The first integral describes the momentum equations, with admissible test functions v i ; the second integral models the incompressibility condition, with admissible test function q. The fluid velocity is denoted by u i , p is the pressure, ρ is the density, ij (u) is the strain rate tensor given by: and σ ij (u, p) is the Cauchy stress tensor for Newtonian fluids: where µ is the dynamic viscosity.
The boundary conditions are defined as: whereû fs i andû i are the prescribed velocities,t i is the traction, and n f j is the normal on the surface pointing outwards. Dirichlet boundary conditions are enforced weakly on the fluid-solid interface, and on the external surfaces via Nitsche's method [49]. The formulation adopted here is the one described in [42]. The surface residual of the external Dirichlet boundaries is defined as: where γ N,u is a penalty parameter. The first integral of (12) results from the integration by parts of the momentum equations; its two terms are denoted as the pressure and viscous standard consistency terms, respectively. The second integral is the addition of pressure and viscous adjoint consistency terms. The third integral introduces an additional penalty term that ensures coercivity of the viscous part of the formulation, and balances the lack of coercivity that is introduced by the viscous standard and adjoint consistency terms [42]. The terms β p and β µ determine whether the adjoint consistency terms use a symmetric formulation (β p = +1, β µ = +1), or a skew-symmetric formulation (β p = −1, β µ = −1). In this work, we use the symmetric variant for the viscous adjoint consistency term because it leads to smaller errors compared to the skew-symmetric variant, as reported by [50]. For the pressure adjoint consistency term, we use the skewsymmetric variation because it consistently controls the mass conservation, u i n f i = 0 [51,42]. A similar treatment is applied on the fluid-solid interface for the term r fs u,p by usingû fs i instead ofû i in (12). The formulation of the penalty parameter γ N,u is taken from [42] and defined as: The terms in (13) account for viscous-dominated and convective-dominated flows, respectively. The term α N,u is a constant problem-dependent penalty term, and the term u ∞ is the infinity norm evaluated at each integration point. The influence of this penalty term will be studied later in this paper.
The residual contribution from external Neumann boundary conditions is defined as: An auxiliary indicator field is introduced to identify isolated volumes of fluid surrounded by the solid domain. These "puddles" lead to an ill-conditioned system of equations because the absolute value of the pressure is not governed. To stabilize the system, we add the following penalty formulation to enforce an average pressure: where k p is a scaling factor. The fieldψ serves as a binary indicator with a value of 1 in isolated fluid regions and zero everywhere else, as shown in Figure 2; the equations to model the field will be defined in Section 4.2. The influence of the term (15) on the conservation of mass will be studied in Section 5.3.

Advection-Diffusion Equation
Species transport is modeled by an advection-diffusion equation. Similar to the Navier-Stokes equations (6), we denote the weak form as r u,c , and decompose it into volumetric and surface contributions: where r Ω u,c and rΩ u,c are the residuals of the volumetric contributions, non-stabilized and stabilized, respectively; the term r D u,c represents the Dirichlet boundary conditions at the external boundaries; the terms r N u,c and r N,fs u,c describe the Neumann boundary conditions at the external boundaries, and at the fluid-solid interface, respectively; and the term r GP u,c models the ghost-penalty formulation. Similar to (5), the residual rΩ u,c is defined in Section 4.4, and the residual r GP u,c is defined in Section 4.5. The non-stabilized volumetric residual contribution, r Ω u,c , is defined as: where c is the species concentration, w is an admissible test function, u i is the vector of fluid velocities,q Ω is the volumetric source, and J i (c) is the diffusive flux defined as: where k is the isotropic diffusivity.
Dirichlet and Neumann boundary conditions are imposed on the external boundaries, Γ f , and Neumann conditions on the fluid-solid interface, Γ fs , as: where n fs i is the normal vector on the fluid-solid interface pointing towards the solid phase,ĉ is a prescribed concentration value, andq Γ andq fs Γ are prescribed species flux values. The weak enforcement of Dirichlet boundary conditions is modeled using Nitsche's method [49]. The residual contributions of the Dirichlet conditions at the external surfaces are defined as: where the penalty parameter α N,c is a problem-dependent constant. The first, second, and third terms correspond to the standard consistency term, the adjoint consistency term, and the penalty term of the Nitsche formulation, respectively.
The Neumann contributions at the fluid-solid interface and the external boundaries are defined as:

Auxiliary Indicator Field
The auxiliary field is modeled as a linear diffusion problem in the fluid domain. The residual of the weak form is defined as: where the residual of the Dirichlet boundary conditions, r D ψ , is formulated in the same way as in (22). Dirichlet boundary conditions are imposed on all inlets and outlets, by settingψ to 0. Adiabatic boundary conditions are imposed on the fluid-solid interface.
The volumetric residual contribution, r Ω ψ is defined as: where ξ is an admissible test function, h ψ is the convection coefficient, and ψ ∞ is the reference indicator value. The diffusion coefficient in (18) is set to 1. The parameters h ψ and ψ ∞ are set to 0.01 and 1, respectively, so that fluid channels that are connected to the inlets and outlets will have an indicator field value close to 0, while the fluid "puddles" will have a value close to 1.
A smooth-Heaviside projection scheme is applied to the indicator field to project the values of the solution either to 0 or to 1, and is defined as:ψ where the sharpness of the projection increases with the parameter k w , and k t ψ ∞ is the threshold. We adopt the values of k w = 1000 and k t = 0.99 to effectively turn theψ term into a binary switch, where a value of 0 corresponds to a point connected to the inlet and/or outlet ports, and a value of 1 corresponds to an isolated volume of fluid.

Spatial Discretization
The governing equations in the fluid phase are discretized in space by the XFEM. This study adopts a generalized enrichment strategy based on the Heaviside-step enrichment of [52], which interpolates consistently the solution fields in the presence of small features, and does not suffer from artificial coupling of disconnected phases. This particular approach has been used by [53,43], [34], [36], and [35], who considered linear elasticity, incompressible Navier-Stokes, linear diffusion, and advection-diffusion problems, respectively.
Here we discretize the fluid state variables, the species concentration, and the indicator field in the fluid domain by the XFEM. Using the symbol ω to represent any of these state variables, the approximation of ω within an element,ω, can be written as follows: where l is the enrichment level, N l is the maximum number of enrichment levels, N i (x) are the nodal basis functions, ω i l is the degrees-of-freedom of enrichment level l at node i in the fluid phase, respectively. The Kronecker delta, δ i lk , selects the degree-of-freedom that is active at node i. At any given point, only one degree-of-freedom per node is used to interpolate the solution, ensuring that the partition of unity is satisfied.
Multiple enrichment levels, i.e. sets of shape functions, may be necessary to interpolate the state variables in multiple, physically disconnected regions of the same phase, c.f. [54], [55], and [53]. When interpolating the level set field by element-wise linear functions in a structured grid, a maximum of 14 enrichment levels is needed in 3D problems [43]. In order to accurately integrate the weak form of the governing equations by Gauss quadrature, intersected elements are decomposed into tetrahedrons using Delaunay triangulation. The reader is referred to [53] for more details on the particular XFEM implementation used in this study.
The Heaviside-step enrichment formulation (27) is ill-defined for cases in which the fluid-solid interface lies exactly on a node, i.e. the level set value φ i at node i equals 0. To avoid this issue, we adopt the level set perturbation approach outlined in [56] and [36]. If the magnitude of the level set value at a node is smaller than some critical value, φ c , the level set value is modified to a shift value, φ s . This perturbation results in the fluid-solid interface moving away from the node, solving the singularity issue. In this study, we adopt the values of φ c = φ s = 10 −6 × h. Numerical studies have shown that the influence of this perturbation is negligible for the problems considered here [57].

Subgrid Stabilization
The convective terms in the incompressible Navier-Stokes and advection-diffusion equations may cause spurious node-to-node velocity oscillations. Furthermore, the equal-order approximations used for u i and p may give rise to spurious pressure oscillations. To prevent these numerical instabilities, we augment the incompressible Navier-Stokes equations with the Streamline Upwind Petrov-Galerkin (SUPG) and the Pressure Stabilized Petrov-Galerkin (PSPG) stabilization formulations introduced by [58]. The stabilized volumetric residual contribution of (5), rΩ u,p , is defined as: where Ω e denotes the set of all elements in the domain Ω, and the stabilization terms τ SUPG,u and τ PSPG are defined in [58].
The stabilized volumetric residual contribution of the advection-diffusion model (16), rΩ u,c , uses the SUPG method, and is defined as: where the stabilization terms τ SUPG,c is defined in [59].

Face-oriented Ghost-penalty Methods
As the geometry of the design evolves during the optimization process, the fluid-solid interface geometry may lead to intersection configurations where certain degrees-of-freedom interpolate in very small subdomains. This produces an ill-conditioning of the system, which manifests itself through an increase in the condition number of the linearized system, and may slow down or prevent the convergence of the nonlinear problem.
To guarantee stability, as well as to improve the conditioning of the system, face-oriented ghost-penalty stabilization terms are used in the vicinity of the fluid-solid interface, c.f. [40]. The ghost-penalty terms for the residual contribution of the incompressible Navier-Stokes equations are defined as: where r GP,µ u,p , r GP,p u,p , and r GP,u u,p are the viscous, pressure, and convective ghost-penalty formulations, respectively.
We adopt the viscous face-oriented ghost-penalty formulation as proposed by [41]: where γ GP,µ is a penalty parameter defined as: and α GP,µ is a constant scaling factor. The jump operator is defined as: and is evaluated at the facet between two adjacent elements, Ω 1 e and Ω 2 e . This formulation overcomes the issue of having a small ratio of volumes on elements bisected by the interface because the domain of integration is the entire edge, regardless of the intersection configuration. As illustrated in Figure 3, the set Ξ belonging to the domain Ω contains all facets F in the immediate vicinity of the fluid-solid interface, for which at least one of the two adjacent elements is cut by the interface.
To control pressure instabilities due to a violated inf-sup condition for equal-order approximations used for u i and p [42], a pressure ghost-penalty stabilization term is applied: where γ GP,p is a penalty parameter defined as: and accounts for the viscous and convective flow regimes, c.f. [40]. The term α GP,p is a constant scaling parameter.
For high Reynolds number flows, a convective ghost-penalty formulation was proposed by [60] to have sufficient control over the convective derivative, u i ∇u i , of the incompressible Navier-Stokes equations. This formulation is defined as: where the parameter γ GP,u is a penalty factor defined as: and α GP,u is a constant scaling parameter. Additional ghost-penalty measures have been proposed in the literature, for example, to control instabilities arising from the incompressibility constraint. However, these additional formulations are not considered in this paper because previous studies have not revealed any further improvement for the laminar flow situations analyzed here, as stated by [42]. To stabilize the concentration field (17), we use the formulation from [61]: where γ GP,c is a penalty parameter defined as: and α GP,c is a scaling constant.
The ghost-penalty formulation for the auxiliary indicator field is identical to (38), except that it operates on a different admissible test function, ξ, a different set of degrees-of-freedom, ψ, and a different scaling factor, α GP,ψ .
The values for α GP,µ , α GP,p , α GP,u , and α GP,c are set on a per-problem basis. The value of α GP,ψ is set to 0.05 for all numerical examples, in accordance to the parameter used by [61] for a linear diffusion field.

Numerical Examples
In the following, we study the characteristics of the proposed CutFEM topology optimization framework for steady-state and transient laminar flow and species transport problems in 3D. We first compare the accuracy and convergence of the analysis of the proposed CutFEM framework against a body-fitted problem taken from the literature. Then, we study the effects on the mass conservation caused by the penalty formulation introduced in (15). Finally, we apply topology optimization to several problems in order to study the characteristics of the framework with respect to different physical phenomena. Unless otherwise stated, geometric and material parameters are given in non-dimensional and self-consistent units.
In this study, time integration is performed by a two-step backward differentiation scheme. A sufficiently large time step is chosen to simulate steady-state conditions. For all time steps, equilibrium is satisfied by solving the associated nonlinear system of equations via Newton's method. Linear problems are solved by the Generalized Minimal RESidual (GMRES) iterative method [62], with an Incomplete LU factorization with dual Threshold (ILUT) preconditioner [63]. For all examples studied in the following, the nonlinear and linear problems are considered converged if the relative residuals are less then 10 −6 .
The optimization problems (1) are solved by a gradient-based algorithm in the reduced space. The gradients of the objective and constraint functions with respect to the optimization variables, s i , are computed via the adjoint method. In this work, we adopt the discrete adjoint formulation for nonlinear fluid and coupled systems of [11] and [64]. The problems are solved via the Globally Convergent Method of Moving Asymptotes (GCMMA) of [65]. The GCMMA parameters are given in Table 1. The optimization problem is considered converged if the change of the objective function relative to the previous objective value is less than 10 −6 , and all constraints are satisfied.

Optimization criteria
Here we summarize the design criteria used to formulate the objectives and the constraints in the optimization problems studied in this section. These include: Drag Coefficient: The drag coefficient is used to qualify the forces of laminar flow on a surface, and is defined as: where e i is a unit vector pointing in the direction of the flow velocity, u c is the characteristic velocity, and L c is the characteristic length.
Mass Flow Rate: The mass flow rate criterion computes the mass of the fluid that passes through a surface per unit of time, and is defined as:ṁ Pressure difference: The total pressure criterion measures the sum of the static and dynamic pressures over a surface: Volume: The volumes of the fluid and solid domains are computed as: Surface Area: The surface area criterion is computed at the fluid-solid interface, and is defined as: Reducing the value of this measure, either with a penalty in the objective functional or with a constraint, discourages the emergence of small geometric features and oscillatory shapes in the optimization problem. While a surface area penalty does not allow explicitly controlling the local shape and the feature size [43], it has been reported effective in regularizing flow optimization problems [17,35].
Target Scalar Value: To measure the maximum difference between a current species concentration and a target species concentration, c ref , over the fluid phase, we use the Kreisselmeier-Steinhauser function [48]: A large value of β KS increases the accuracy of approximating the maximum, but may result in large design sensitivities that affect the convergence of the optimization problem.

Verification of the CutFEM analysis
In this first example, we seek to verify the analysis capabilities of our CutFEM optimization framework. The validity of our approach is measured in terms of the accuracy of the flow solutions with respect to solutions from the literature, and in terms of the convergence of the flow solutions with respect to different levels of mesh refinement and different values for the Nitsche penalty parameter, α N,u . We study the 3D-1Z problem from [66], which considers a 3D laminar steady-state flow around a cylinder. No indicator fields are modeled in this example because this configuration does not include isolated volumes of fluid and the geometry does not change during the analysis.
The problem setup is shown in Figure 4. The inflow condition is:   where the maximum inflow velocity is set to u m = 0.45. A traction-free boundary condition is imposed on the outlet. No-slip boundary conditions are imposed on the surface of the cylinder and on all other planes. The characteristic velocity is defined by the average inflow velocity, that is u c = 4/9 u m = 0.2. The characteristic length, L c , is defined as the diameter of the cylinder, 0.1, which yields a Reynolds number of 20. Note that in our CutFEM approach the cylinder is represented by the level set function (4).
The following quantities are used to monitor mesh convergence: the drag coefficient (40) around the cylinder, c D , and the total pressure drop (42) between the inlet and the outlet planes, T in − T out . The numerical solutions provided in the study by [66] give the values of 6.05 and 6.25 as the lower and upper bounds for the drag, respectively. These bounds were computed from the numerical results provided by several research groups through different numerical schemes, such as Finite Difference, Finite Volume, and Finite Element Methods, among others, and mesh convergence studies. No reference solution is provided for the total pressure difference; however, given that this measure is used as the objective of the optimization problems in several numerical examples ahead, we include it in the convergence study.
The mesh for the body-fitted problem utilizes a boundary layer around the cylinder, with 128 elements on the surface and 64 layers. The width of the first layer is 1.0832 × 10 −9 , and the exponential growth factor of each subsequent layer is 1.2. The number of elements on the inlet and outlet surfaces is 32 × 32 and 16 × 16, respectively. We do not consider symmetry, and solve the problem on the entire flow domain. For the body-fitted mesh, all boundary conditions are imposed in the strong form, and no face-oriented ghost-penalty formulation is applied. The total number of degrees-of-freedom is 1, 556, 066. The remaining parameters used for this problem are shown in Table 2.
The meshes for the CutFEM convergence study are constructed using a local hierarchical refinement strategy along the fluid-solid interface. We deliberately do not perform adaptive mesh refinement nor generate boundary layer meshes in the CutFEM approach, as we do not consider these techniques in the topology optimization problems ahead. The first level of mesh refinement has an element size of h = 1.14 × 10 −2 Value Element size h = 1.14 × 10 −2 , 3.79 × 10 −3 , 1.27 × 10 −3 , 4.22 × 10 −4 Characteristic velocity u c = 0.2 Characteristic length L c = 0.1 Dynamic viscosity µ = 0.001 Density ρ = 1 Nitsche velocity penalty α N,u = 10, 10 2 , 10 3 , 10 4 Viscous ghost-penalty α GP,µ = 0.05 Pressure ghost-penalty α GP,p = 0.005 Convective ghost-penalty α GP,u = 0.05 Pressure constraint parameter k p = 0  for elements located at x < 1 and h = 3.41667 × 10 −2 for all other elements. Subsequently, we only refine elements that are intersected by the surface of the cylinder; we perform this process thrice. In consequence, we generate the four refinement levels shown in Figure 5. Note that the hierarchical mesh refinement utilized here does not provide the same resolution as the boundary layer meshing scheme above.
The hierarchical mesh refinement scheme leads to a larger number of elements compared to the boundary layer approach. To reduce the computational cost, we model half of the flow domain and impose symmetry boundary conditions by settingû z = 0 along the plane z = 0.205. The inflow and outflow conditions are the same as in the body-fitted problem setup. All boundary conditions are enforced weakly. The total number of degrees-of-freedom for the finest mesh is 8, 019, 736. The face-oriented ghost-penalty methods increase the bandwidth of the sparse matrix of the system, which may increase the linear solve time; see also discussion on the computational costs below. The remaining parameters used for this problem are shown in Table 3.
The results for the body-fitted and the CutFEM problems are shown in Figure 6. The drag coefficient for the body-fitted problem is 6.169, well within the lower and upper bounds established in the study by [66]. The total pressure difference is 0.0213. The results of the CutFEM analysis are shown in Figure 6a.   6b and display a higher convergence rate than the drag coefficient. The total pressure drop characterizes the global flow solution, while the drag is a local measure along the cylinder surface that depends on the spatial gradients of the flow field. The total pressure difference reaches the same solution as its body-fitted counterpart. Spurious mass flow rates through the surface of the cylinder are shown in Figure 6c. We observe that as the mesh is refined, and the value of the Nitsche velocity penalty is increased, the spurious mass flow converges to zero within numerical precision.
As the mesh is refined, the influence of the Nitsche velocity parameter in (13) vanishes. The relative maximum difference between the body-fitted problem and the CutFEM solutions at the finest mesh is 0.4% for the drag coefficient and 0.03% for the total pressure drop. These results suggest that our CutFEM framework is sufficiently accurate for topology optimization purposes. However, our discretization scheme may suffer from inaccuracies in local quantities that strongly depend on the resolution of boundary layer phenomena. If such quantities are used in the formulation of the optimization problem, we recommend that the performance of the optimized design is verified using a body-fitted mesh with a resolved boundary layer.
To provide insight into the computational costs, wall clock timing results for individual steps in the forward analysis with the body-fitted mesh and the CutFEM mesh with h = 1.27 × 10 −3 are given in Table 4, using 60 and 120 cores. These particular meshes are chosen as they have comparable numbers of elements and unconstrained DOFs. For the CutFEM analysis, the case with α N,u = 10 is considered; other values  for α N,u lead to similar results. The wall clock times for assembling and solving the linearized systems are averaged over all Newton iterations. For both, the body-fitted and the CutFEM analysis the same convergence criteria for the nonlinear and linear problems are applied, as defined above; the ILUT fill is set to 7.0. The computations were performed on a cluster equipped with Quad-Core AMD Opteron 2378 Processors. These results show that the computational effort for building and updating the CutFEM model is insignificant when compared to the time spent for solving the nonlinear system. Assembly of the CutFEM model is more costly than the body-fitted mesh due to the integration of the ghost-penalty terms and the increased number of integration points of intersected elements. However, the assembly operation scales well with the number of cores. The timing results presented in Table 4

Verification of the Average Pressure Constraint
In this second example, we seek to verify the penalty formulation in (15) with respect to the accuracy to which mass conservation is satisfied. We model a steady-state flow through a bent pipe and measure the relative mass flow rate difference between the inlet and the outlet. The problem setup is shown in Figure  7. The presence of an isolated spherical inclusion of fluid causes a singular analysis problem because the absolute value of the pressure is not governed within it. We compare flow solutions with and without the spherical inclusion and study the influence of the pressure penalty parameter, k p .
The design domain is discretized by a uniform structured mesh, and the geometries of the pipe and the spherical inclusion are described by the LSFs presented in [1]. We only model half of the domain, and apply symmetry boundary conditions along the plane z = 0.5. We study the pressure penalty formulation in (15) and, for illustration purposes, a formulation in which the pressure penalty is applied over the entire fluid domain, rather than on the isolated "puddles" exclusively.
The problem parameters are given in Table 5. The inflow condition is:    No fluid "puddle", no penalty Fluid "puddle", penalty over Ω f Fluid "puddle", penalty by indicator field Figure 8: Comparison of the influence of the average pressure constraint on the relative mass flow rate difference between the inlet and the outlet for the bent pipe example.
where y c and z c are the coordinates at the center of the inflow, as defined in Figure 7. The characteristic velocity is u c = 200, and the characteristic length is defined as the diameter of the pipe, L c = 0.2, for a Reynolds number of 40. A traction-free boundary condition is imposed on the outlet. No-slip boundary conditions are imposed on the surfaces of the pipe and of the sphere.
The relative error between the inlet and outlet mass flow rates is shown in Figure 8. The flow without the sphere results in a relative mass flux error of 0.03%. We obtain the same error when the sphere is included and if we apply the pressure penalty formulation (15) exclusively to the isolated volume of fluid through the use of the indicator field in (24). The error is insensitive to the penalty parameter value, which is set to k p = 1 for all numerical examples below. Applying the penalty formulation over the entire domain, similar to the approach used by [43] for linear elasticity problems, can cause a significant error in the mass conservation if a large value of k p is chosen; see Figure 8.

Design of a Manifold with Multiple Outlets
In this example, we apply the CutFEM framework to the design of a steady-state flow bend, with multiple inlets and outlets. We seek to minimize the total pressure drop between the inlets and the outlets at steadystate while controlling the mass flow rates through the outlet ports. The example is the 3D analog to the 2D problem found in [67]. The problem setup is shown in Figure 9. The design domain has two inlets (on the left and right), and four outlets (on all other planes). The inflow condition is formulated in the same way as (47), the characteristic velocity is u c = 200, and the characteristic length is L c = 1. Traction-free boundary conditions are imposed on the outlets. No-slip boundary conditions are imposed on the fluid-solid interface. We only model an eighth of the domain, and symmetry boundary conditions are imposed on the planes with x = 3.5, y = 3.5, and z = 3.5.
In addition to minimizing the total pressure drop, we penalize the surface area of the fluid-solid interface. The use of the surface area as a contribution to the objective function has been applied previously to species transport topology optimization by [35] to improve the smoothness of the final design, and to regularize the optimization problem. The objective is defined as: where the superscript "0" denotes the values of the initial design, the subscript i denotes the i-th inlet or outlet, and w S is a constant scaling factor. The design is subject to a 5% volume constraint of the fluid domain to suppress trivial solutions, and to promote the formation of distinct fluid channels: Further, we wish to impose a constraint such that the amount of mass flow exiting through each outlet is the same. Given that the GCMMA algorithm does not allow equality constraints, we impose inequality constraints with lower and upper limits on the mass flow rates. The upper and lower bounds are set to 25% ± 1.25%, respectively, where the tolerance value of ±1.25% was chosen in order to not overconstrain the optimization problem. The constraints are defined as follows: whereṁ in,i andṁ out,i are the mass flow rates at the i-th inlet and i-th outlet, respectively, and N out is the number of outlets.
The remaining parameters are given in Table 6. The design domain is initialized with 5 × 5 × 5 spherical solid inclusions of radii 0.5, as shown in Figure 10. In our experience, the flow topology optimization problems studied here are rather insensitive to the initial design as long as the number of inclusions is sufficiently large. Additional mechanisms for seeding solid inclusions could be added to the proposed CutFEM framework, such as topological derivatives [68]; however, they are outside the scope of this study.
The converged design after 99 iterations is shown in Figure 11, which resembles conceptually the 2D results from [67]. Figure 12a shows the convergence plots of the objective and the volume constraint. The initialization process described above leads to an initial design that violates the volume constraint. Initially, the objective increases while the volume constraint value decreases. Once the constraint is satisfied, the objective decreases until a feasible minimum is found. The objective changes from a normalized value of 1.1, where the mass constraint was violated, to 1.56. Figures 12b and 12c show the convergence plots for the lower and upper bounds of the mass inequality constraints. We can observe that the mass flow rate constraints are also satisfied, and the amount of fluid flow exiting through each outlet is virtually the same. The mass flow rate at the inlets is 6.24, while the rates at each of the outlets is 1.56.

Design of a Manifold with Variable Outlets
In this example, we study a fluid flow problem with one inlet and multiple outlets, similar to the one from Example 5.4, with the caveat that the outlets are allowed to vary in position and shape. The problem setup is shown in Figure 13. Fluid flows into the domain through the top inlet, and exits through the 2 outlets on the left and right, and the 2 outlets on the bottom. The outlets are described by cylinders using the level set parametrization of (4). We introduce three additional design variables per outlet that define the in-plane coordinates at the center of the outlets and their radii; the lower and upper bounds are given in Table 8. The interior of the domain is parametrized by (3), and the lower and upper bounds of the design variables are the same as in Table 6. The inflow, outflow, and the interface boundary conditions are the same as in Example 5.4. The problem parameters are given in Table 7. We model a quarter of the domain, and impose symmetry boundary conditions along the x = 1.2 and the z = 0.25 planes.
We seek to minimize the total pressure drop between the inlet and outlet ports and the surface area of the fluid-solid interface, and impose a 20% fluid volume constraint (49). Further, we formulate constraints on the mass flow rates such that the left and right outlets have 33.333% each, and the bottom outlets have 16.667% each, of the fluid flow entering through the inlet. Similar to (50), we formulate these limits as inequality constraints and set the tolerances to ∆ṁ 1,2 = ∓0.333% for the left and right outlets and to ∆ṁ 3,4 = ∓0.167% for the bottom outlets. As we assume symmetry and model only one quarter of the design domain, constraints for only one outlet on a vertical face, i.e. Outlet 1 or 2, and one outlet at the bottom face, i.e. Outlet 3 or 4, need to be imposed; see Figure 13. This leads to the following four additional constraints: We impose these constraints via a continuation method and gradually contract the tolerances from initially ∆ṁ 1,2 = ∓1.667% and ∆ṁ 3,4 = ∓0.833% to their target values in the course of the optimization process.
The design domain is initialized with 3 × 6 × 1 spherical solid inclusions of radii 0.125. The optimized design after 350 iterations is shown in Figure 14. The mass flow through each outlet, and the radii and coordinates of the outlets are given in Table 9. During the optimization process, the left and right outlets increase their  Table 9: Mass flow rate, fraction of the mass flow entering through the inlet, radii, and coordinates of the outlets in the variable outlets example. radii and move up in the y-direction. The bottom outlets merge and form a single outlet. The numbers in Table 9 reveal that these outlets do not fully overlap, and that they have moved away from the x = 1.2 and z = 0.25 planes of symmetry. The relative difference between their final positions and the axes of symmetry is small nonetheless, in the order of 0.1%. This value got increasingly closer to 0 as the mass flow rate tolerances of (52) were lowered. The mass flow exiting through each of the left and right outlets is 33.2% of the mass flow entering through the inlet, and the mass flow exiting through each of the bottom outlets is 16.8%; these values are within the lower and upper bounds of the mass flow limits.

Design of a Micromixer
In this example, we apply our CutFEM framework to the modeling and optimizing of a micromixer at a low Reynolds number and steady-state conditions. The example is the 3D analog to the 2D problem found in [35], and is similar to the micromixer studies with flow topology optimization found in [69], [70] and [71]. The problem setup is shown in Figure 15. A "red" fluid and a "blue" fluid enter the design domain through the left inlet and exit it through the lower right side. We assume that the fluids are ideally miscible and have identical flow properties. The "red" fluid is represented by a species concentration value ofĉ = 1, and the "blue" fluid byĉ = 0. We do not consider diffusion of the species field through the solid phase. The inflow condition is formulated using the same approach as in Example 5.3. No-slip boundary conditions are applied at the fluid-solid interface, and a traction-free boundary condition is imposed on the outlet. We only model one half of the domain, and symmetry boundary conditions are imposed on the z = 1.5 plane. An adiabatic condition is imposed on the fluid-solid interface for the species field.  The objective uses the target scalar value formulation from (45) at steady-state, and is defined as: with β KS = 400, and c ref = 0.5. A small surface area penalty with w S = 0.001 is applied to regularize the optimization problem. The design is subject to a volume constraint (49) of 35% to suppress trivial solutions, and to promote the formation of distinct fluid channels. Similar to [35], a constraint is imposed on the maximum pressure drop to prevent the formation of small geometric features: The problem is initialized with 15 × 5 × 5 spherical solid inclusions of radii 0.1, similar to the previous examples. The remaining parameters are given in Table 10.
The optimized material layout, along with the species field and the velocity magnitude, is shown in Figure  16. The average species concentration at the outlet is 0.46. Analogous to the 2D results from [35], the length of the channel layout increases by creating an intricate wavy design, thereby increasing the path traveled by the fluids. This lengthening mechanism is the key to enhance the mixing of the fluids in 3D laminar flows. The channel does not contain internal isolated solid particles. The number of iterations is rather large: 1, 500 iterations are required to form the channel layout and fine-tune its shape. This behavior has also been seen by [35], who attributes this to the interplay of localized sensitivities along the fluid-solid interface and the volume constraint on the fluid phase. In contrast to the density approach of [72], we do not obtain numerical artifacts in the optimized design.

Design of an Oscillating Pump
In the final example, we study the applicability of the CutFEM framework to problems with transient behavior. We consider the design of a simplified fluid pump. The example is the 3D analog to the 2D       problem found in [73]. The problem setup is shown in Figure 17. The basic idea is to prescribe a harmonically oscillating velocity at the left port, and optimize the design domain to maximize the amount of fluid that is transported through the port at the top. The port at the bottom represents a reservoir from which additional fluid can enter the domain.
where N t is the total number of time iterations. The design domain is subject to a volume constraint (49) of 15% to suppress trivial solutions and to promote the formation of smooth fluid channels.
We apply traction-free boundary conditions at the top and bottom ports. The flow velocity at the left port is prescribed as: where the inflow frequency is set to v = 100π. No-slip boundary conditions are imposed on the fluid-solid interface. We only model one half of the domain, and symmetry boundary conditions are imposed on the plane z = 0.5. We model a single pumping cycle with 40 time steps per cycle. Numerical experiments showed that this number of time steps is sufficient to model the harmonic behavior of the flow. The design is initialized with the same layout as in Example 5.4. The remaining problem parameters are given in Table  11.
The design converged in 225 iterations and the optimized layout is shown in Figure 18. Similar to the results from [73], the design exhibits a narrowing of the inflow channel, and the formation of a central reservoir from which fluid flows towards the pumping outlet during the outflow cycle. The pipe joining the bottom port to the top outlet forms naturally during the optimization process. During the inflow cycle, fluid flows  from the left port towards the outlet port, while additional fluid is pulled from the bottom port towards the central reservoir. An area of recirculation then forms around the central reservoir. During the outflow phase, fluid flows from the central reservoir towards the outlet port, while fluid from the bottom port flows in a vortex like path around the left inlet; this fluid is then transported towards the outlet port during the next inflow phase. This process is shown by streamlines in Figure 19. To verify the harmonic behavior of the pump for the optimized design, we modeled 5 pumping cycles of 40 time steps each; the corresponding mass flow rates for all ports are shown in Figure 20. We can observe that the optimized design maintains its oscillatory behavior. The average mass flow rate through the top port of the optimized design is −0.328, with the negative sign indicating that the flow leaves the domain. In contrast to the results shown by [73], we did not observe numerical artifacts in the optimized material layout of our pump design, such as isolated regions of fluid flow.

Conclusions
We presented and studied a design optimization framework based on the CutFEM method [1], which combines the LSM, the XFEM, and face-oriented ghost-penalty methods. This framework builds upon previous optimization studies that utilized the LSM and the XFEM. The approaches presented in this paper add flexibility, accuracy, and robustness through weakly enforcing boundary conditions via a Nitsche formulation and by applying face-oriented ghost-penalty methods to the velocity and pressure fields, discretized by a generalized Heaviside-step enrichment strategy. An auxiliary indicator field was introduced to identify isolated volumes of fluid, and to enforce a constraint on the average pressure on these "puddles" to prevent a singular analysis problem.
The geometry of the fluid-solid interface was described by an explicit LSM where the parameters of the parametrized LSF are defined as functions of the optimization variables. The resulting optimization problem was solved in the reduced space by a nonlinear programming method. The flexibility of this scheme allowed us to define optimization variables in addition to the parameters resulting from a finite element discretization of the LSF. We demonstrated this feature by defining the location and shape of the outlet ports as optimization variables.
The accuracy of our CutFEM approach was verified quantitatively through comparison against benchmark studies, and qualitatively through application to topology optimization problems. The design optimization problems converged to intuitive designs and resembled well the results from previous 2D studies. A rigorous comparison of the proposed approach against topology optimization methods using Brinkman penalization was beyond the scope of this paper. However, such a comparison should be performed in future studies to    As demonstrated by the numerical studies presented in this paper, our CutFEM approach is accurate, robust, and applicable to a broad range of design problems for low Reynolds number flows. However, the computational costs can be significant, in particular if accurate flow solutions are required. This is in parts due to the enlarged bandwidth of the linear systems caused by the face-oriented ghost-penalty formulation, but more importantly due to the lack of an adaptive mesh refinement strategy. Future research may focus on integrating mesh refinement strategies into our CutFEM framework such the boundary layer phenomena are captured accurately and efficiently.