A regularization scheme for explicit level-set XFEM topology optimization

Regularization of the level-set (LS) field is a critical part of LS-based topology optimization (TO) approaches. Traditionally this is achieved by advancing the LS field through the solution of a Hamilton-Jacobi equation combined with a reinitialization scheme. This approach, however, may limit the maximum step size and introduces discontinuities in the design process. Alternatively, energy functionals and intermediate LS value penalizations have been proposed. This paper introduces a novel LS regularization approach based on a signed distance field (SDF) which is applicable to explicit LS-based TO. The SDF is obtained using the heat method (HM) and is reconstructed for every design in the optimization process. The governing equations of the HM, as well as the ones describing the physical response of the system of interest, are discretized by the extended finite element method (XFEM). Numerical examples for problems modeled by linear elasticity, nonlinear hyperelasticity and the incompressible Navier-Stokes equations in two and three dimensions are presented to show the applicability of the proposed scheme to a broad range of design optimization problems.


Introduction
Topology optimization (TO), with its large design freedom, has emerged as a powerful design tool for a variety of applications including structural mechanics, fluid flow and heat transfer [1][2][3]. The two most commonly used TO approaches are density-based methods and level-set (LS)based implicit boundary methods. Since the introduction of the level-set method (LSM) [4], the method has gained great popularity in the areas of image processing, computer graphics, computational geometry and computational physics [5][6][7]. LSMs describe geometry changes by evolving an implicit boundary, conventionally defined as the zero-level iso-contour of a level-set function (LSF), fðX Þ. When applied to TO, LSMs enable a clear and unambiguous definition of the material interface [2,8]. van Dijk et al. [8] classified LSMs into two broad categories based on the LSF update procedure: i) Implicit methods where some form of the Hamilton-Jacobi (HJ) equation is used to evolve the LSF based on a velocity field defined by shape sensitivities that are in turn governed by the physics [9][10][11][12], and ii) explicit methods where a parametrized LSF is updated using mathematical programming techniques [8,[13][14][15][16][17].
By construction, shape sensitivities only exist in the vicinity of the domain boundary, i.e., zero iso-contour of the LSF, and depend on the spatial LSF gradient. Furthermore, locally too flat or too steep LSF gradients affect the stability and the rate of convergence of LSMs, while a uniform and uniquely defined LSF improves those features. To this end, several regularization schemes for LSMs have been proposed [8]. Perimeter regularization is used to obtain a well-posed optimization problem [18] whereas Tikhonov regularization is used to control the smoothness of the LSF gradient [19,20]. These methods, however, do not guarantee a unique LSF and may lead to a flat LSF [8]. For both, explicit LS descriptions and implicit LSMs using the HJ equation, there exists a strong need for better regularization schemes to improve convergence of the optimization process and avoid convergence to local minima [8].
An LSF with a uniform gradient along the interface, e.g., jrfðX Þj ¼ 1, is desired for both explicit and implicit approaches [5,8]. A uniform, unit norm gradient is a unique property of the signed distance function (SDF), and hence the SDF is commonly used to reinitialize the LSF. However, as it evolves, the LSF quickly loses its SDF characteristics [21]. To alleviate this issue, several SDF regularization techniques have been proposed [22][23][24][25][26][27][28][29]. For the implicit LSMs, typically an auxiliary HJ equation [25] is solved, or a fast marching method [26] is used to reinitialize the LSF intermittently during the optimization process. For a more detailed discussion of HJ methods, the interested reader is referred to Refs. [5,30]. Even though LSF reinitialization is widely used, it slightly moves the zero LS iso-contour during the reinitialization process [5,31] and therefore affects the convergence of the design optimization process. Moreover, if the HJ equations are solved by an explicit time integration scheme, the design step size is limited by the Courant-Friedrichs-Lewy stability criterion.
To overcome the LSF reinitialization issues discussed above, an energy functional that penalizes the deviation from a unit norm LSF gradient is most commonly added to the objective function [13,[22][23][24]32]. However, these measures do not allow for sufficient control of the LSF gradient away from the interface and may lead to oscillations [27,28]. Double-well energy functionals were introduced in Refs. [27][28][29] to enforce a unit norm gradient near the interface and a zero LSF gradient away from the interface, thus leading to an LSF that has an SDF characteristic near the interface and is constant away from it. The fundamental limitation of these local LS regularization approaches is that they operate directly on the local design LS value, or its gradient, and lack information about the minimum distance of a point to the interface. Thus, they cannot distinguish between points that have the same LS value but differ in their distance to the interface. This inability may cause undesirable material inclusions away from the original interface as illustrated in Fig. 1. Enforcing a unit norm gradient of the LSF at points with LS values close to zero but away from the interface at iteration n ( Fig. 1(a)) may create new intersections at iteration n þ 1 ( Fig. 1(b)), without these intersections being necessarily beneficial for improving the performance of the design.
Using the local design LS value, either a zero gradient is achieved away from the interface or a unit norm gradient is enforced in the vicinity of it (see Fig. 1(a) for a simplified one-dimensional LSF and the corresponding interface). Due to the local enforcement of the different targets, the local measure lacks the ability to distinguish between areas where the LSF is within the LS bound f Bnd (i.e., -f Bnd £f£f Bnd ) at which an interface exists and areas at which no interface is present. The local regularization scheme enforces jrfj ¼ 1 anywhere, where -f Bnd £f£ f Bnd and thus has the tendency to create spurious material inclusions ( Fig. 1(b)). Limiting the enforcement of the unit norm gradient to only intersected elements does not reliably alleviate this issue. Based on the authors' experience, adding energy functionals to the objective while formulating a well-posed optimization problem is challenging and requires problem dependent fine-tuning of parameters defining the regularization.
To overcome the shortcomings of the previously discussed regularization schemes, this paper introduces a novel LS regularization approach for explicit LS-based TO. This approach penalizes the difference between the design LSF fðX Þ and a target LSF. The target LSF is constructed at every design iteration from the SDF, which is obtained using an extension of the heat method (HM) [33,34]. The SDF is computed at every design iteration for the current interface geometry and treated as a prescribed target field for the design LSF. The governing equations of the physics models and the HM are discretized in space by the extended finite element method (XFEM). The advantage of this novel approach is that a smooth and unique target field is used as reference for a locally enforced LS regularization. Using a differentiable penalty formulation to match the design LSF with the computed target LSF alleviates the need for reinitialization and therefore does not introduce discontinuities in the optimization process. The implementation of the HM is straightforward and can easily be coupled with explicit LSMs. The effectiveness of the proposed scheme is demonstrated via LS-based TO numerical examples in two (2D) and three (3D) spatial dimensions. Linear elasticity, nonlinear hyperelasticity and fluid flow examples are studied to demonstrate the general applicability of the novel LS regularization scheme.
The remainder of the paper is organized as follows: Section 2 gives a brief summary of the explicit LS-based TO framework; Section 3 introduces the HM; Section 4 Fig. 1 Effect of a local LS regularization scheme on a simplified one-dimensional design problem for two distinct design iterations: (a) n and (b) n + 1 discusses the XFEM discretization; Section 5 provides details of the explicit LS regularization scheme; numerical examples are presented in Section 6; and Section 7 concludes and summarizes the paper.

Explicit level-set topology optimization
The LSM implicitly describes the geometry of a body and its evolution using a scalar LSF fðX Þ. The material layout within a design domain Ω D composed of two distinct phases is given by where Ω I and Ω II are the material domains of Phases I and II, respectively, such that between them is denoted by Γ I,II which corresponds to the zero LS iso-contour fðX Þ ¼ 0. In explicit LS-based TO, the nodal values of the discretized LS field f i ðX Þ are defined as an explicit function of the design variables.
Here, the design variable field is discretized using linear finite element (FE) shape functions, and each node j is assigned one design variable s j . The LSF function is defined by filtering the discretized design variable field as follows: where N n is the number of nodes within a filter radius r f and jX i -X j j is the Euclidean distance between Nodes i and j. This linear filtering scheme Eq. (2) initially proposed by Ref. [35] for LSMs increases the area of influence of every design variable and therefore enhances convergence of the optimization problem. The LSF is then used to discretize distinct physical sub-domains using the XFEM (see Section 4). The optimization problems considered in this work are formulated as where s denotes the vector of design variables and u is the vector of state variables. The objective function consists of a weighted contribution z 1 that characterizes the physical performance (e.g., strain energy, total fluid pressure drop) and two weighted penalty contributions, such that w 1 þ w 2 þ w 3 ¼ 1. All objective function contributions are normalized by their values of the initial design, denoted by the superscript 0. The first penalty contribution minimizes the perimeter to avoid the emergence of irregular geometric features, and the second penalty is the newly introduced LS regularization that is discussed in detail in Section 5. When no LS regularization is used, w 3 ¼ 0. All optimization problems considered in this work are subject to a volume constraint γ V on Ω I which prevents trivial solutions. The lower and upper bounds of the design variables are denoted by s L and s U , respectively. The number of design variables is denoted by N s and the number of state variables is N u . The state variables u are governed by a set of discretized partial differential equations. These equations are satisfied for each design in the optimization process.

The heat method
The basis of the proposed regularization scheme for explicit LSMs is the construction of the SDF at every design iteration. Most commonly used approaches to obtain the SDF are fast marching methods [36] and fast sweeping approaches [37]. These methods however require non-trivial implementations within an FE-based software platform and present issues with parallelization [33]. In the current work, the SDF is obtained using an extension of the HM. First, a transient heat conduction equation is solved on the entire design domain, with a heat source at the material interface ( Fig. 2(a)). The strong form of the governing equation of the temperature field ðX Þ is where Δ is the Laplace operator, and the temperature time derivative is denoted by _ . The initial and Dirichlet boundary conditions at the XFEM interface are 0 ðX Þ ¼ 0 and ðΓ I,II Þ ¼ 1, respectively. Adiabatic boundary conditions are applied to the temperature field at the domain boundary.
The distance field f D ðX Þ is obtained by solving a Poisson's equation with a volumetric flux that depends on the normalized gradient of the temperature field ðX Þ. The governing equation of the distance field f D ðX Þ is with the Dirichlet boundary condition f D ðΓ I,II Þ ¼ 0 at the material interface and homogeneous Neumann boundary conditions at the domain boundary. The solution of Eq. (5) is illustrated in Fig. 2(b). The SDF f SD ðX Þ is then computed by multiplying the distance field with the sign of the original design LSF fðX Þ (Fig. 2(c)).
Mathematically, this is stated as Note that due to enforcing f D ¼ 0 at the interface, constructing the SDF from the distance field via Eq. (6) does not introduce spurious fluctuations.
From a design optimization point of view, it is important to distinguish between the design LSF fðX Þ which is an explicit function of the design variables s (Eq. (2)) and the reconstructed SDF f SD ðX Þ. The design LSF fðX Þ determines the decomposition of the design domain into distinct phases and the material interface, which is the starting point of the HM and the SDF computation.
The weak form of the residual equation of the temperature field ðX Þ is where the admissible test functions are denoted by δ. The time derivative at the current time step m + 1 is approximated using an implicit Euler backward scheme as where m is the temperature field at the previous time step m and Δt is the time step size. To obtain an accurate distance field, Crane et al. [33] recommended a time step size in the order of Δt ¼ h 2 where h is the element edge length. To increase computational efficiency of the HM, only a single time step is used for solving the temperature field ðX Þ of Eq. (7). The time step size Δt is set sufficiently large to obtain a non-zero temperature gradient in the entire design domain and a meaningful SDF. This greatly reduces the computational overhead compared to a fully transient problem while only slightly effecting the accuracy of the obtained SDF away from the interface. For selecting a sufficiently large, single time step size, the following guideline can be used: where L is defined as L ¼ L max ffiffiffi d p and L represents a small temperature value at the far end of the design domain, e.g., L ¼ 1Â10 -4 . The maximum side length of the design domain bounding box is denoted by L max and d is the spatial dimensionality, e.g., d ¼ 2 in 2D. When employing the HM within a TO process, computational efficiency is more important than a high accuracy of the SDF. As discussed in Section 6.3.1, numerical studies have shown that solving Eq. (7) only for a single time step does not impede the functionality of the LS regularization but significantly simplifies the application of the HM for TO.
It is not necessary to compute the distance field f D ðX Þ and the signed distance field f SD ðX Þ in two sequential steps. The residual equation for the signed distance field f SD ðX Þ is obtained by integration by parts of Eq. (5) and stated as where the coupling term G is computed as the normalized temperature gradient times the negative of the sign of the design LS field: It should be noted that Eq. (10) does not contain any boundary contribution, as rf SD ¼ G is assumed over the outer domain boundary.

The extended finite element method
The XFEM [38] is used for discretization of the physics and the HM governing equations on a non-conforming background mesh. Being an immersed boundary method, it alleviates the need for re-meshing, which can be challenging and computationally costly during an optimization process.
Enrichment of the classical FE approximation spaces with additional shape functions is used for interpolation into disconnected sub-domains [39]. Multiple levels of enrichment are used to avoid spurious coupling or load transfer between disconnected material sub-domains. In this work, a generalized Heaviside enrichment strategy [40] is employed where the degrees of freedom within each Fig. 2 Construction of the SDF using the HM. (a) A heat distribution is obtained from heat sources at the material interface; (b) the normalized temperature gradient is utilized to compute a distance field; from that, (c) the SDF is obtained unique sub-domain are approximated by standard FE shape functions. The HM state variables u ¼ f, f SD g at Node i are therefore approximated as where the Heaviside function H is a function of the LS value and is defined as The maximum number of enrichment levels is denoted by M, N e n is the number of nodes per element and N k ðX Þ is the elemental shape function. The Kronecker delta δ k mq selects the active enrichment Level q for Node k such that the partition of unity principle is satisfied. More details regarding Heaviside enriched XFEM can be found in Refs. [41,42].
Face-oriented ghost penalization, as proposed by Refs. [43,44], is used to stabilize the XFEM discretization. Numerical instabilities arise in the XFEM when the material interface Γ I,II moves too close to a FE node, leading to a vanishing zone of influence of certain degrees of freedom. Face-oriented ghost stabilization cures this illconditioning independent of the intersection configuration. For stabilization of the solution fields, face-oriented ghost penalization is applied in the vicinity of the interface. It is formulated as where γ G is the ghost penalization parameter and F cut contains all element faces in the immediate vicinity of the material interface for which at least one of the two adjacent elements is intersected [45]. The jump operator is defined as〚〛¼ j Ω e 1 -j Ω e 2 . The ghost penalty is evaluated along all faces between two adjacent elements, Ω e 1 and Ω e 2 . The outward facing normal vector between Ω e 1 and Ω e 2 is denoted by N e . This form of stabilized XFEM is also referred to as CutFEM in Ref. [46]. Boundary conditions and interface conditions are applied weakly in this work using the unsymmetrical version of Nitsche's method [47]. Weakly enforced boundary and interface conditions are essential in LSbased XFEM TO where the material phase of a domain boundary at which Dirichlet boundary conditions are applied may change. The weakly enforced conditions are applied using where the phase index is denoted by p ¼ fI,IIg and N denotes the normal vector on the domain or interface boundary. The first term in Eq. (15) is the standard consistency, the second term is the adjoint consistency, and the last term is a penalty term on the jump of the state variables. The Nitsche penalty parameter is denoted by γ N . The same XFEM approach, stabilization and application of Dirichlet boundary conditions via Nitsche's method as outlined in this section is also used for discretization of all physics governing equations discussed in Section 6.

Explicit level-set regularization
The SDF obtained by the HM is used for regularization of the design LSF during the optimization process. Instead of reinitializing the design LSF fðX Þ with the SDF f SD ðX Þ, the following penalty formulation is proposed to achieve a continuous LS regularization: where f Bnd denotes an upper (lower) bound for the target LSF. The penalty measures the difference between the design LSF fðX Þ and a target LSFfðX Þ, as well as the difference in the spatial gradients. The target LSFfðX Þ is constructed from the SDF f SD ðX Þ. As the design LSF fðX Þ converges to the target LSFfðX Þ, both contributions in Eq. (16)  To avoid length scale dependence, the penalty terms are normalized.
In this work, upper and lower bounds are imposed on the design LSF. To obtain a bounded design LSF away from the material interface, the following piecewise linear LSF is proposed as a target: The truncated target LSFfðX Þ matches the SDF in the vicinity of the interface. Away from the interface, the lower bound -f Bnd is matched in Ω I and the upper bound f Bnd is matched in Ω II . To avoid the non-differentiability of the piecewise linear target LSFfðX Þ, a smooth target LSF fðX Þ is used to approximate Eq. (17). This is achieved by the following sigmoid function: A comparison between the piecewise linear target LSF and the smooth target LSF is illustrated in Fig. 3 for a onedimensional interface configuration.
The target LSFfðX Þ depends implicitly on the geometry of the interface, defined by the zero iso-contour of the design LSF fðX Þ that depends explicitly on the optimization variables. The implicit dependence offðX Þ is described by the governing equations of the HM. In general, these implicit and explicit dependencies need to be considered for computing consistent design sensitivities of the LS regularization penalty Eq. (16). However, if the weight w 3 for the LS regularization term in the formulation of the objective function Eq. (3) is large and the implicit dependency of the target LSFfðX Þ on the optimization variables is accounted for, the evolution of the design may be dominated by the LS regularization and the optimization process may converge to a design with poor physical performance. To overcome this issue, the LS regularization weight should be chosen small, e.g., w 3 <0:1; a motivation for this recommendation will be presented in Section 6.2.1.1.
In addition, the authors found it advantageous to consider the target LSFfðX Þ as a prescribed field which is reconstructed at every design iteration of the optimization process. Using this interpretation offðX Þ, the penalty term Eq. (16) depends only explicitly on the optimization variables and the implicit sensitivities are ignored. As it will be shown in Section 6.2.1.2, this approach reduces the influence of the LS regularization term on the evolution of the design LSF in the vicinity of the zero iso-contour, i.e., the interface geometry. In addition, ignoring the implicit sensitivity contributions enhances the convergence of the optimization process to a design with optimized physical performance. Furthermore, the computational cost is noticeably reduced by omitting the computation of the adjoint solution of the HM. The LS regularization mainly controls the slope of the LSF along the interface and ensures that the LSF converges to either upper or lower bounds, AEf Bnd , away from the interface. The reader may note that the LS regularization penalty Eq. (16) provides non-zero sensitivities in the entire design domain. This is usually not the case in LS-based TO using the XFEM, where shape sensitivities only exist in the vicinity of the interface [8].
Instead of introducing the LS regularization by the penalty term Eq. (16) into the formulation of the optimization problem, one could also add a constraint in the form of p£ε with ε< <1. As the authors obtained satisfactory results with the penalty formulation for a wide range of penalty weights (see Section 6), the constraint formulation has not been considered in this work.

Numerical examples
Numerical examples considering different physical phenomena in 2D and 3D are presented in this section to study the characteristics of the proposed regularization approach. The examples include structural design problems modeled by linear and nonlinear elasticity and a flow problem described by the incompressible Navier-Stokes equations at steady-state.
Common to all examples is that the governing equations, including the ones of the HM, are discretized by the XFEM. An iterative Newton-Raphson scheme is used to solve the nonlinear problems that are considered as converged when a relative nonlinear residual norm drop greater than 10 6 is achieved. A single load increment is used in all numerical examples. The linearized sub-systems are solved using the Multifrontal Massively Parallel Solver (MUMPS) [48]. Bilinear four node quadrilateral elements are used in 2D, and trilinear eight node hexahedral Fig. 3 Piecewise and smooth approximation of the design LSF elements are used in 3D. The same discretization is used for the design LSF. The parameter optimization problem Eq. (3) is solved using a gradient-based nonlinear programming scheme, namely the Globally Convergent Method of Moving Asymptotes (GCMMA) [49], with no inner iterations.
The optimization problem is considered converged if the constraint is satisfied and the relative change in objective between two consecutive design iterations is less than 1Â10 -5 . The design sensitivities are obtained using the adjoint method. For more details on design sensitivities using XFEM, the interested reader is referred to Ref. [50]. Selective structural springs [51] are applied for all structural problems to suppress rigid body motion of isolated material domains that may emerge in solid-void LS-based TO. The initial seeding of the design domain with holes is performed using primitives in the form of squares in 2D and cubes in 3D.
The parameters used for the following numerical examples are listed in Table 1. Self-consistent units are used for all numerical examples and therefore not stated explicitly. The bounds for the design and target LSF are set as a function of the element edge length h. Note that the bound f Bnd for the target LSFfðX Þ is within the range of values of the discretized design LSF fðX Þ. As discussed in Section 3, the temperature field in the HM is advanced in time by a single time step, unless otherwise stated. A staggered solution approach is used to separately solve the two partial differential equations of the HM in a one-way coupled fashion. This improves computational efficiency as each sub-problem is linear and of smaller size.

Examples for linear elasticity
The physical response in the first set of examples is described by linear elasticity. The weak form of the governing equation is with u ¼ u on Γ u , where the displacement vector is denoted by u. The surface traction applied on Γ T is T . The infinitesimal strain tensor is defined by ε ¼ 1 2 ðru T þ ruÞ, and the Cauchy stress is σ ¼ Dε. The fourth-order material tensor is denoted by D, and for isotropic, linear elastic homogeneous materials considered in this work it is defined in terms of the Lamé constants l and as follows: where δ ij is the Kronecker delta. The Lamé constants can be expressed in terms of the Young's modulus E and the Poisson's ratio as The problem parameters used for linear elastic problems are listed in Table 2. For more details about the linear elastic XFEM formulation used in this section, the reader is referred to Ref. [52]. 6

.1.1 Hanging bar in 2D
As a first design example, a 2D linear elastic plane stress "hanging bar" design optimization problem is solved. This example problem is a modified version of the two-bar truss solid-void problem that was studied in Ref. [2] with the density method. The initial design of size 80 Â 40 with boundary conditions is shown in Fig. 4(a). Only one half of the design is analyzed and optimized, with weakly enforced symmetry boundary conditions along the vertical center line. The top edge of the domain is clamped while a traction load of T X 2 ¼ -30 is applied in X 2 direction at the center of the bottom edge over a length of 12. The region at the center of the bottom edge, at which the load is applied, is excluded from the design domain. The optimization problem Eq. (3) is to minimize the strain energy with a perimeter penalty and an LS regularization penalty subject to a volume constraint of γ V ¼ 0:16. A relative GCMMA optimization step size of 0.1 and a time step size of Δt ¼ 4 is used in the HM.
The final design is shown in Fig. 4(b) and consists of only a single vertical bar, to transfer the applied traction load at the bottom to the support at the top of the domain.   Figure 5 shows a comparison of the evolution of objective and constraint with and without LS regularization. Early on in the design process, oscillations are observed without LS regularization, while with LS regularization a smooth design evolution is obtained. Moreover, the design problem converges significantly faster when LS regularization is applied: About 300 design iterations versus about 500 design iterations. Since the LS regularization contribution vanishes at the optimum, the regularized design problem converges to the same objective and constraint values as without regularization. Figure 6(a) shows snapshots of the design LSF at discrete steps during the optimization process, comparing the evolution of the design LSFs with and without LS regularization. The LSFs are plotted over X 1 at the top of the design domain, i.e., X 2 ¼ 40. When no LS regularization is used, irregularities of the design LSF are observed as the design is changed and a quick degradation of the slope of the LSF at the interface is seen. With regularization, a non-oscillatory LSF is obtained in the entire design domain. Even though the design problem without regularization also eventually converges, the LSF is noisy and at some locations close to zero. Due to numerical noise, LS values close to zero often create unintended isolated material islands in the vicinity of the interface. These may cause ill-conditioning of the analysis and potential divergence of the design problem. Thus, when using the proposed LS regularization, a much larger optimization step size can be used, owing to the enhanced stability of the optimization problem. Figure 6 on the right shows the design LSFs at design iteration 200 without regularization applied ( Fig. 6(b)) and with LS regularization (Fig. 6(c)). Even though the same zero iso-contour is obtained, the non-regularized LSF is noisy, and the initial design can still be observed even at convergence. In contrast, LS regularization achieves a design LSF with LS values at the target boundary values f Bnd and a unit norm gradient in the vicinity of the interface. Due to the exclusion of the bottom center from the design domain, slight irregularities in the final isocontour and in the regularized LSF are seen in this region.

Hanging bar in 3D
A 3D configuration of the 2D problem discussed in Section 6.1.1 is considered here. The initial design with loads and boundary conditions is shown in Fig. 7(a). Due to the symmetry of the design problem, only one quarter of the domain is modeled and optimized. Again, the area at the bottom of the domain at which the traction load is applied (circular area of radius 8.5) is excluded from the design As in the 2D problem, the optimization process converges to a single vertical bar (Fig. 7(b)). Figure 8 shows the design LSFs of the final design obtained without and with LS regularization. The 2D plane shown here is taken along the diagonal of the design domain as indicated in red in Fig. 7(b).
It can clearly be seen that when LS regularization is used the LS values are at the target bounds away from the interface, while a smoothed LSF with a unit norm gradient is obtained near the solid-void interface. Without regularization, the LSF of the initial design is clearly preserved at convergence, and large spatial oscillations exist throughout the entire design domain. As before, slight oscillations in the design LSF are observed in the vicinity of the loaded edge since this domain is excluded from the design domain. Overall, increased stability and higher convergence rates are observed for this initial set of design problems when the LS regularization is applied.

Examples for nonlinear hyperelasticity
To demonstrate the applicability of the proposed LS regularization scheme to design problems with increased complexity, examples are considered next where the structural response is described by a finite strain hyperelastic model. The weak form of the governing equation is stated as  Fig. 8 where F ¼ ∂x=∂X is the deformation gradient tensor, and x ¼ u þ X describes the relationship between reference X and current x coordinates. The first Piola-Kirchhoff stress is denoted by P. A hyperelastic Saint Venant-Kirchhoff constitutive model for homogeneous, isotropic compressible materials is used, which is formulated as where S is the second Piola-Kirchhoff stress tensor and E is the Green-Lagrange strain tensor. The second order identity tensor is denoted by I. The Lamé constants defined in Eq. (21) are used. The Green-Lagrange strain tensor is defined as: where the right Cauchy-Green deformation tensor C is computed as C ¼ F T F. Finally, the first Piola-Kirchhoff stress is obtained from P ¼ FS: (25) For more details on the formulation and verification of the nonlinear XFEM formulation used in this work, the interested reader is referred to Ref. [53]. The parameters listed in Table 3 are used for all hyperelastic design optimization problems, unless stated otherwise.

Beam in 2D
First, we consider the design of a beam-type structure in 2D. The initial design with loads and boundary conditions is shown in Fig. 9(a). The design domain is of size 240Â40. A traction load of T X 2 ¼ -10 is applied at the top center of the domain over a length of 3 while the structure is simply supported on either ends on the bottom of the domain. Due to the symmetry of the design problem, only half of the domain is modeled and optimized. Symmetry boundary conditions are applied weakly. The structural response is described by the hyperelastic model outlined above and discretized by the XFEM. Following the work of Ref. [54], a plane strain model is used in 2D. The objective of the optimization problem is to minimize the strain energy with a 1% penalty weight on the perimeter and a 1% penalty weight on the LS regularization. The optimization problem is subject to a volume constraint of   Fig. 9 (a) Initial design of the 2D beam problem with boundary conditions and dimensions; (b) comparison of the zero LS iso-contours of the final designs without and with LS regularization γ V ¼ 0:6, and an optimization step size of 0.025 is used. A single time step of size Δt ¼ 35 is used in the HM. Figure 9(b) shows a comparison of the zero LS isocontour of the final beam designs obtained without and with LS regularization. The typical truss-like structure is obtained for both methods with only slight differences in the final geometries. These differences can be attributed to different evolutions of the LSFs during the optimization process. The proposed LS regularization scheme leads to an increased convergence behavior due to a uniform LS gradient in the vicinity of the interface. This is also reflected in a slightly (0.1%) lower strain energy of the regularized design versus the non-regularized design.
The design LSFs, fðX Þ, at the final optimization iteration are shown in Figs. 10(a) and 10(b) without LS regularization and with regularization, respectively. Both LSFs are warped for illustration purposes. An oscillatory LSF is obtained without regularization, while with LS regularization the optimization process converges to a smoothly truncated design LSF. The regularized LSF shows a unit norm gradient in the vicinity of the material interface while having a zero gradient away from the interface. Figure 10(c) shows the SDF obtained by the HM for the final design of the 2D beam. It can be seen that overall the SDF is well resolved. Only in areas with small geometric features, with a size of h, the XFEM discretization insufficiently resolves the SDF. Consequently, the LS regularization suffers in these areas from a degraded target LSF due to the limited resolution of spatial discretization.

Influence of the LS regularization penalty weight
The influence of different weights w 3 for the LS regularization penalty is studied in Fig. 11. Figure 11(a) shows the evolution of strain energy and Fig. 11(b) shows the LS regularization penalty for regularization weights of w 3 ¼ f0:01,0:05,0:1,0:5g. With an increased LS regularization penalty weight, the minimization of the LS regularization term is favored early in the design process, while the minimization of strain energy is given less importance. The reader may note small jumps in the evolution of strain energy and the LS regularization penalty, for example, at iteration 350 and iteration 400. The jumps are caused by thin structural members disconnecting. The design iteration at which this happens depends on the weight of the LS regularization term.
For a weighting parameter in the range of 1%£w 3 £10%, both the strain energy and the regularization penalty assume similar values after about 200 design iterations. If the LS regularization weight is too large (e.g., 50%), the optimization problem changes noticeably and the physical performance of the optimized design is affected. The LS regularization term then dominates the overall objective and the physics contribution is of lesser importance ( Fig. 11 (a)). In the authors' experience with the current problem and other design problems, an LS regularization weight up to 10% provides a good balance between sufficient regularization while not impairing the performance of the optimized design. As discussed in Section 5, the proposed regularization scheme considers the target LSFfðX Þ as a prescribed field and ignores the implicit contributions of the penalty term Eq. (16) to the design sensitivities. Only the explicit dependency of the design LSF on the optimization variables is accounted for in the sensitivity analysis. To illustrate the benefits of this approach, the influence of including the implicit design sensitivities is investigated. The implicit contributions are computed by the adjoint approach. Figure 12 shows the optimized beam design obtained with an LS regularization weight of w 3 ¼0:1 and including implicit sensitivities of the target LSF. Due to a fairly large weight of the regularization on the objective, the implicit design sensitivities influence significantly the evolution of the zero LS iso-contour. The design evolution is predominantly influenced by the regularization scheme and insufficiently driven by the physics performance. This leads to spurious void inclusions, premature convergence, and poor physical performance of the optimized structure (Fig. 12). For a sufficiently low regularization penalty (e.g., w 3 ¼0:01) these issues are not observed, and the design convergence is indistinguishable from the one where the implicit design sensitivities are omitted. Thus, to prevent an undesired influence of the regularization on the design evolution and to gain computational efficiency, it is recommended to ignore design sensitivities of the target LSF on the design variables and to use a low penalty weight for the regularization term.

Beam in 3D
The hyperelastic beam example is studied next in 3D to demonstrate the applicability of the proposed LS regularization scheme to 3D problems where the geometry undergoes significant changes during the optimization process. Figure 13(a) shows the initial design with boundary conditions for a design domain of size 240Â40Â40. An element edge length of h¼2 is used for this example, along with a time step size of Δt ¼ 51 in the HM. Analogous to the 2D configurations, a traction load of T X 2 ¼ -2 is applied within a circular region of radius 2 at the center of the top face of the domain. The structure is simply supported at all four corners at the bottom face of the design domain. The loading and the support domains are excluded from the design domain, respectively.  Two-fold symmetry is exploited for the mechanical and the design problem. A relative optimization step size of 0.0125 is used and a volume constraint of γ V ¼0:3 is enforced through a continuation approach. Figure 13(b) depicts the final design obtained after 400 optimization iterations, using the proposed LS regularization. As before, a smooth evolution of objective and constraint is achieved when employing the LS regularization scheme. Figure 14 shows the design LSFs at convergence obtained without the LS regularization and with LS regularization, extracted in the center of the design domain (indicated by the red plane in Fig. 13(b)). As before, a clear difference can be observed with respect to the smoothness of the LSF and the uniformity of the spatial gradient along the zero LS iso-contour. While the non-regularized LSF is shallow, the regularized LSF quickly approaches the LS bounds away from the interface. Figure 14 shows thin vertical members in the optimal design, which represent shear webs in between the top and bottom flanges of the beam. Because their thickness is in the order of the mesh size h, the XFEM discretization provides insufficient resolution of both the stress and strain fields, as well as the SDF. Due to the inability of the discretization to capture the SDF in areas of small features, the target LSF is not well developed and, therefore, the LS regularization suffers. The result of this can be seen in Fig. 14(b) where the lower LS bound of -f Bnd ¼ -2 is not reached by the design LSF within the thin vertical members of the structure. While this effect has already been observed in the 2D beam example (Section 6.2.1) this issue is more pronounced here. This is not an inherent drawback of the proposed LS regularization scheme, but rather stems for the underlying XFEM discretization and its limitations to resolve features with a size in the order of h. Minimum feature size control (e.g., Ref. [55]) or local mesh refinement would be required to properly regularize the design LSF.  The goal of the final example is to demonstrate the applicability of the proposed LS regularization scheme to a flow problem, modeled by the incompressible Navier-Stokes equations at steady state. Furthermore, the influence of the accuracy with which the HM is time integrated is studied.
The non-stabilized weak form of the volumetric contribution (R) of the governing equation is stated as where v is the velocity vector, p is the pressure, and is the density. The admissible test functions for velocity and pressure are denoted by δv and δp, respectively. The infinitesimal strain rate tensor εðvÞ is defined as The Cauchy stress tensor for Newtonian fluids is denoted by σðv,pÞ and is defined as σðv,pÞ ¼ -pI þ 2 D εðvÞ, (28) where D is the dynamic viscosity. The governing Eq. (26) is augmented by sub-grid stabilization to suppress spurious velocity and pressure oscillations, as well as by ghost penalization. For more details on the XFEM discretization and the corresponding stabilization techniques, the reader is referred to Ref. [45]. The fluid domain boundary is decomposed into the fluid-void interface Γ I,II , and Dirichlet and Neumann external boundaries, Γ u and Γ T , respectively. No-slip conditions are applied weakly at the fluid-void interface; the other boundary conditions are problem dependent and are specified below.
An extension of the fluid nozzle problem studied by Refs. [56][57][58], to 3D is studied here. The design domain with boundary conditions and the initial design are shown in Fig. 15(a). The computational domain is a 5Â5Â5 cube with a 0:75Â5Â5 non-design domain downstream from the inlet face. A parabolic velocity profile with a maximum inlet velocity of 30 in X 1 direction is applied to the inlet face, and zero pressure is enforced weakly at the circular outlet of radius 1.25. Both, the inlet domain and the circular outlet face are excluded from the design domain. Only one quarter of the design domain is modeled. Slip conditions are imposed on the X 1 -X 2 and X 1 -X 3 symmetry planes.
Assuming a constant density of ¼1 and a dynamic viscosity of D ¼1, the flow conditions correspond to a Reynolds number of Re ¼ 66 with the reference velocity being the average inlet velocity and the reference length being the edge length of the design domain, i.e., L Ref ¼5. A ghost stabilization parameter of 0.005 is used for stabilization of the pressure and a ghost penalization parameter of 0.05 is used for stabilization of velocities; for details see Ref. [45].
The objective of this nozzle design problem is the minimization of the total pressure drop between inlet and outlet along with a 1% perimeter penalty and a 5% LS regularization penalty. The optimization problem is subject to a γ V ¼0:3 volume constraint on the fluid phase. A relative GCMMA step size of 0.01 and a single time step of Δt ¼0:1 is used for the HM. The main problem parameters are listed in Table 4.
The final design obtained after 90 design iterations is shown in Fig. 15(b). These results agree with the ones presented in the literature. As before, a smooth evolution of objective and constraint is obtained when LS Fig. 15 (a) Initial design of the 3D fluid nozzle with boundary conditions and dimensions; (b) final nozzle design. The diagonal highlighted in red represents the plane in which the design LSFs are being compared regularization is used. The design LSFs at convergence without and with the LS regularization are compared in Fig. 16. Again, a rather oscillatory LSF with flat gradients near the interface is obtained when no LS regularization is employed. With LS regularization, the design LSF assumes the target bounds away from the interface while having a unit norm gradient in the vicinity of the fluid-void interface. The final designs obtained without and with the HM do not differ significantly. However, improved numerical stability and robustness during the optimization process was observed due to the regularization of the LSF.
To provide insight into the dependence of the optimization results on the accuracy with which the temperature field in the HM is time integrated (Eq. (7)), the number of time steps in the Euler backward scheme Eq. (8) is varied. The total time is kept constant at t max ¼1. Comparisons of objective and LS regularization penalty evolution for different number of time steps of the HM Eq. (7) are shown in Fig. 17. No significant differences are observed when solving the HM with multiple time steps and reduced time step sizes. This confirms the observations by Ref. [33], and shows that a single time step is sufficient for LS regularization using the HM.

Conclusions
A regularization scheme for explicit LS XFEM design optimization in 2D and 3D was presented. The regularization scheme augments the objective function by a penalty term that measures the difference between the design LSF and a target LSF, both in regard to the field value and its spatial derivatives. The target LSF has a unit norm gradient in the vicinity of the interface and assumes either an upper or lower bound away from the interface, depending on the material phase. The target LSF is constructed from the SDF that is computed by an XFEM discretization of the HM at every design iteration for the current interface geometry. Numerical experiments on 2D and 3D problems in solid and fluid mechanics showed that the proposed regularization scheme is largely insensitive to the penalty weight for the regularization term. As long as the weights are less than 10%, the LS regularization does not influence noticeably the final design. A small influence on the design evolution has been observed for larger penalty weights. Furthermore, it was observed that it is beneficial to ignore the dependence of the target LSF on the interface geometry for computing the design sensitivities.
Omitting the sensitivities of the target LSF on the design variables leads to an improved convergence to a design with improved physical performance and reduces the computational cost. The numerical results further suggest that the temperature field of the HM can be computed by a single time step without significantly affecting the accuracy of the SDF. The time step size is a function of the domain length. Good results were obtained with a time step size determined by Eq. (9).
Comparing the results obtained with and without the proposed regularization scheme suggests that the proposed scheme significantly improves the convergence of the optimization process and mitigates issues in the XFEM analysis due to the emergence of small inclusions of one phase within a domain occupied by another phase. The scheme mitigates irregular interface evolution and promotes a uniform LS gradient at the zero LS iso-contour. It  The numerical studies presented in this paper have revealed a few shortcomings of the proposed method that need to be addressed in future work. These include overcoming the limited resolution of a fixed XFEM background mesh with linear interpolation. When features of dimensions in the order of the element edge length of the background mesh emerge in the optimization process, the resolution of a linear interpolation is insufficient to accurately compute the target LSF. Therefore, the performance of the regularization is reduced. This could be addressed by adding feature size control to the design problem, locally refining the background mesh or by using higher order spatial discretizations. In addition, an increase in computational cost was observed due to the need for solving two additional partial differential equations in the HM. Depending on the complexity of the physics model, this additional cost may become significant when compared to runs without the LS regularization scheme. However, the additional cost is partially offset by an increased convergence rate and by the ability to use larger optimization step sizes. Future work needs to address ways to improve the computational efficiency of the scheme. For example, the temperature and SDF fields in the HM could be solved only approximately, using a few iterations of an iterative linear solver. Finally, due to the regularization of the LSF, the nucleation of new holes is impaired. In order to mitigate the dependency of the final design on the initial seeding, systematic approaches for the creation of new holes during the optimization process need to be explored in combination with the proposed LS regularization. This includes but is not limited to using additional LSFs [59] or topological derivatives.