A regularization scheme for explicit level-set XFEM topology optimization

Abstract

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.

Key words： level-set regularization; explicit level-sets; XFEM; CutFEM; topology optimization; heat method; signed distance field; nonlinear structural mechanics; fluid mechanics

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 [¹–³]. 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) [⁴], the method has gained great popularity in the areas of image processing, computer graphics, computational geometry and computational physics [⁵–⁷]. LSMs describe geometry changes by evolving an implicit boundary, conventionally defined as the zero-level iso-contour of a level-set function (LSF), \varphi \left(\mathit{X}\right). When applied to TO, LSMs enable a clear and unambiguous definition of the material interface [²,⁸]. van Dijk et al. [⁸] 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 [⁹–¹²], and ii) explicit methods where a parametrized LSF is updated using mathematical programming techniques [⁸,¹³–¹⁷].

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 [⁸]. Perimeter regularization is used to obtain a well-posed optimization problem [¹⁸] whereas Tikhonov regularization is used to control the smoothness of the LSF gradient [¹⁹,²⁰]. These methods, however, do not guarantee a unique LSF and may lead to a flat LSF [⁸]. 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 [⁸].

An LSF with a uniform gradient along the interface, e.g., \left|\nabla \varphi \left(\mathit{X}\right)\right|=\mathrm{1}, is desired for both explicit and implicit approaches [⁵,⁸]. 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 [²¹]. To alleviate this issue, several SDF regularization techniques have been proposed [²²–²⁹]. For the implicit LSMs, typically an auxiliary HJ equation [²⁵] is solved, or a fast marching method [²⁶] 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. [⁵,³⁰]. Even though LSF reinitialization is widely used, it slightly moves the zero LS iso-contour during the reinitialization process [⁵,³¹] 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 [¹³,²²–²⁴,³²]. However, these measures do not allow for sufficient control of the LSF gradient away from the interface and may lead to oscillations [²⁷,²⁸]. Double-well energy functionals were introduced in Refs. [²⁷–²⁹] 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+\mathrm{1} (Fig. 1(b)), without these intersections being necessarily beneficial for improving the performance of the design.

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

Full size|PPT slide

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 {\varphi }_{\mathrm{Bnd}} (i.e., -{\varphi }_{\mathrm{Bnd}}\le \varphi \le {\varphi }_{\mathrm{Bnd}}) at which an interface exists and areas at which no interface is present. The local regularization scheme enforces \left|\nabla \varphi \right|=\mathrm{1} anywhere, where -{\varphi }_{\mathrm{Bnd}}\le \varphi \le {\varphi }_{\mathrm{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 \varphi \left(\mathit{X}\right) 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) [³³,³⁴]. 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 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 \varphi (\mathit{X}). The material layout within a design domain {\mathrm{\Omega }}_{\mathrm{D}} composed of two distinct phases is given by

(1)

where {\mathrm{\Omega }}_{\mathrm{I}} and {\mathrm{\Omega }}_{\mathrm{II}} are the material domains of Phases I and II, respectively such that {\mathrm{\Omega }}_{\mathrm{D}}={\mathrm{\Omega }}_{\mathrm{I}}\cup {\mathrm{\Omega }}_{\mathrm{II}}. The interface between them is denoted by {\mathrm{\Gamma }}_{\mathrm{I},\mathrm{II}} which corresponds to the zero LS iso-contour \varphi (\mathit{X})=\mathrm{0.0}. In explicit LS-based TO, the nodal values of the discretized LS field {\varphi }_i(\mathit{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:

(2){\varphi }_i={\displaystyle \frac{{\displaystyle \sum _{j=\mathrm{1}}^{N_n}{w}_{ij}{s}_j}}{{\displaystyle \sum _{j=\mathrm{1}}^{N_n}{w}_{ij}}}},w_{ij}=\max (\mathrm{0},r_{\mathrm{f}}-\left|{\mathit{X}}_i-{\mathit{X}}_j\right|),

where N_{\mathrm{n}} is the number of nodes within a filter radius {\mathit{r}}_{\mathrm{f}} and \left|{\mathit{X}}_i-{\mathit{X}}_j\right| is the Euclidean distance between nodes i and j. This linear filtering scheme Eq. (2) initially proposed by Ref. [³⁵] 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:

(3)\begin{array}{l}\underset{s}{\min }z(\mathit{s},\mathit{u})=w_{\mathrm{1}}{\displaystyle \frac{z_{\mathrm{1}}(\mathit{s},\mathit{u})}{z_{\mathrm{1}}({\mathit{s}}^{\mathrm{0}},{\mathit{u}}^{\mathrm{0}})}}+w_{\mathrm{2}}{\displaystyle \frac{{\displaystyle {\int }_{{\mathrm{\Gamma }}_{\mathrm{I,II}}}\mathrm{d}}A}{{\displaystyle {\int }_{{\mathrm{\Gamma }}_{{}_{\mathrm{I,II}}}^{\mathrm{0}}}\mathrm{d}}A}}+w_{\mathrm{3}}{\displaystyle \frac{p(\mathit{s})}{p({\mathit{s}}^{\mathrm{0}})}},\\ \mathrm{s}\mathrm{.t}.g_{\mathrm{1}}(\mathit{s})={\displaystyle \frac{{\mathrm{\Omega }}_{\mathrm{I}}}{{\mathrm{\Omega }}_{\mathrm{I}}+{\mathrm{\Omega }}_{\mathrm{II}}}}-{\gamma }_V\le \mathrm{0},\\ \mathit{s}\in \mathrm{\Pi }=\{{\mathrm{ℝ}}^{N_s}|{\mathit{s}}_{\mathrm{L}}\le \mathit{s}\le {\mathrm{s}}_{\mathrm{U}}\},\\ \mathit{u}\in {\mathrm{ℝ}}^{N_u}\mathrm{.}\end{array}

where \mathit{s} denotes the vector of design variables and \mathit{u} is the vector of state variables. The objective function consists of a weighted contribution z_{\mathrm{1}} that characterizes the physical performance (e.g., strain energy, total fluid pressure drop) and two weighted penalty contributions, such that w_{\mathrm{1}}+w_{\mathrm{2}}+w_{\mathrm{3}}=\mathrm{1.0}. 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_{\mathrm{3}}=\mathrm{0.0}. All optimization problems considered in this work are subject to a volume constraint {\gamma }_V on {\mathrm{\Omega }}_{\mathrm{I}} which will prevent trivial solutions. The lower and upper bounds of the design variables are denoted by {\mathit{s}}_{\mathrm{L}}and {\mathit{s}}_{\mathrm{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 [³⁶] and fast sweeping approaches [³⁷]. These methods however require non-trivial implementations within a FE-based software platform and present issues with parallelization [³³]. 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 \theta (\mathit{X}) is:

(4)\dot{\theta }=\Delta \theta

, where \mathrm{\Delta } is the Laplace operator, and the temperature time derivative is denoted by \dot{\theta }. The initial and Dirichlet boundary conditions at the XFEM interface are {\theta }^{\mathrm{0}}(\mathit{X})=\mathrm{0.0} and \theta ({\mathrm{\Gamma }}_{\mathrm{I,II}})=\mathrm{1.0}, respectively. Adiabatic boundary conditions are applied to the temperature field at the domain boundary.

The distance field {\varphi }_D(\mathit{X}) is obtained by solving a Poisson’s equation with a volumetric flux that depends on the normalized gradient of the temperature field \theta (\mathit{X}). The governing equation of the distance field {\varphi }_D(\mathit{X}) is

(5)\mathrm{\Delta }{\varphi }_{\mathrm{D}}=\nabla \cdot \left(-{\displaystyle \frac{\nabla \theta }{|\nabla \theta |}}\right),

with the Dirichlet boundary condition {\varphi }_{\mathrm{D}}({\mathrm{\Gamma }}_{\mathrm{I,II}})=\mathrm{0.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 {\varphi }_{\mathrm{SD}}(\mathit{X}) is then computed by multiplying the distance field with the sign of the original design LSF \varphi (\mathit{X}) (Fig. 2(c)). Mathematically, this is stated as

(6){\varphi }_{\mathrm{SD}}=\mathrm{sign}(\varphi )\hspace{0.17em}{\varphi }_{\mathrm{D}}.

Note that due to enforcing {\varphi }_{\mathrm{D}}=\mathrm{0.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 \varphi (\mathit{X}) which is an explicit function of the design variables s (Eq. (2)) and the reconstructed SDF {\varphi }_{\mathrm{SD}}(\mathit{X}). The design LSF \varphi (\mathit{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.

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

Full size|PPT slide

The weak form of the residual equation of the temperature field \theta (\mathit{X}) is:

(7)R_{\theta }={\displaystyle {\int }_{\mathrm{\Omega }}\delta \theta \hspace{0.17em}\dot{\theta }\hspace{0.17em}\mathrm{d}V}+{\displaystyle {\int }_{\mathrm{\Omega }}\delta \nabla \theta \hspace{0.17em}\nabla \theta \hspace{0.17em}\mathrm{d}V=\mathrm{0}},

where the admissible test functions are denoted by \delta \theta. The time derivative at the current time step m+1 is approximated using an implicit Euler backward scheme as

(8){\dot{\theta }}^{m+\mathrm{1}}\approx {\displaystyle \frac{{\theta }^{m+\mathrm{1}}-{\theta }^m}{\mathrm{\Delta }t}},

where {\theta }^mis the temperature field at the previous time step m and \mathrm{\Delta }t is the time step size. To obtain an accurate distance field, Crane et al. [³³] recommended a time step size in the order of \mathrm{\Delta }t=h^{\mathrm{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 \theta (\mathit{X}) of Eq. (7). The time step size\mathrm{\Delta }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:

(9)\mathrm{\Delta }t\ge {\left({\displaystyle \frac{L}{\ln {\theta }_L}}\right)}^{\mathrm{2}},

where L is defined as L=L_{\max }\sqrt{d} and {\theta }_L represents a small temperature value at the far end of the design domain, e.g., {\theta }_L=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{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=\mathrm{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 {\varphi }_{\mathrm{D}}(\mathit{X}) and the signed distance field {\varphi }_{\mathrm{SD}}(\mathit{X}) in two sequential steps. The residual equation for the signed distance field {\varphi }_{\mathrm{SD}}(\mathit{X}) is obtained by integration by parts of Eq. (5) and stated as

(10)R_{{\varphi }_{\mathrm{SD}}}={\displaystyle {\int }_{\mathrm{\Omega }}\delta \nabla {\varphi }_{\mathrm{SD}}}(\nabla {\varphi }_{\mathrm{SD}}-G)\hspace{0.17em}\mathrm{d}V=\mathrm{0},

where the coupling term G is computed as the normalized temperature gradient times the negative of the sign of the design LS field:

(11)G=-\mathrm{sign}(\varphi )\left({\displaystyle \frac{\nabla \theta }{|\nabla \theta |}}\right).

It should be noted that Eq. (10) does not contain any boundary contribution, as \nabla {\varphi }_{\mathrm{SD}}=\mathit{G} is assumed over the outer domain boundary.

The extended finite element method

The XFEM [³⁸] 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 [³⁹]. 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 [⁴⁰] is employed where the degrees of freedom within each unique sub-domain are approximated by standard FE shape functions. The HM state variables u=\{\theta ,{\varphi }_{\mathrm{SD}}\} at Node i are therefore approximated as

(12)u_i(\mathit{X})={\displaystyle \sum _{m=\mathrm{1}}^M\left(H(-\varphi (\mathit{X})){\displaystyle \sum _{k=\mathrm{1}}^{N_{\mathrm{n}}^e}{N}_k}(\mathit{X}){\delta }_{mq}^k{u}_{im}^k+H(\varphi (\mathit{X})){\displaystyle \sum _{k=\mathrm{1}}^{N_{\mathrm{n}}^e}{N}_k}(\mathit{X}){\delta }_{mq}^k{u}_{im}^k\right)},

where the Heaviside function H is a function of the LS value and is defined as

(13)

The maximum number of enrichment levels is denoted by M, N_{\mathrm{n}}^e is the number of nodes per element and N_k(\mathit{X}) is the elemental shape function. The Kronecker delta {\delta }_{mq}^k 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. [⁴¹,⁴²].

Face-oriented ghost penalization, as proposed by Refs. [⁴³,⁴⁴], is used to stabilize the XFEM discretization. Numerical instabilities arise in the XFEM when the material interface {\Gamma }_{\mathrm{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 ill-conditioning 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

(14)R_G=h{\gamma }_G{\displaystyle \sum _{F\in F_{\mathrm{cut}}}{\displaystyle {\int }_F〚\delta \nabla u\hspace{0.17em}{\mathit{N}}_{\mathrm{e}}〛〚\nabla u\hspace{0.17em}{\mathit{N}}_{\mathrm{e}}〛\mathrm{d}A}}=\mathrm{0},

where {\gamma }_G is the ghost penalization parameter and F_{\mathrm{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 [⁴⁵]. The jump operator is defined as 〚•〛={•|}_{{\mathrm{\Omega }}_{{\mathrm{e}}_{\mathrm{1}}}}-{•|}_{{\mathrm{\Omega }}_{{\mathrm{e}}_{\mathrm{2}}}}. The ghost penalty is evaluated along all faces between two adjacent elements, {\mathrm{\Omega }}_{{\mathrm{e}}_{\mathrm{1}}} and {\mathrm{\Omega }}_{{\mathrm{e}}_2}. The outward facing normal vector between {\mathrm{\Omega }}_{{\mathrm{e}}_{\mathrm{1}}} and {\mathrm{\Omega }}_{{\mathrm{e}}_2} is denoted by {\mathit{N}}_{\mathrm{e}}. This form of stabilized XFEM is also referred to as CutFEM in Ref. [⁴⁶].

Boundary conditions and interface conditions are applied weakly in this work using the unsymmetrical version of Nitsche’s method [⁴⁷]. Weakly enforced boundary and interface conditions are essential in LS-based 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

(15)R_N^p=-{\displaystyle {\int }_{\mathrm{\Gamma }}\delta u^p\hspace{0.17em}\nabla u^p\hspace{0.17em}\mathit{N}\hspace{0.17em}dA}+{\displaystyle {\int }_{\mathrm{\Gamma }}\nabla u^p\hspace{0.17em}\mathit{N}\hspace{0.17em}〚u^p〛\mathrm{d}A}+{\gamma }_N{\displaystyle {\int }_{\mathrm{\Gamma }}〚\delta u^p〛〚u^p〛\mathrm{d}A}=\mathrm{0},

where the phase index is denoted by p=\{\mathrm{I,II}\} 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 {\gamma }_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 \varphi (\mathit{X}) with the SDF {\varphi }_{\mathrm{SD}}(\mathit{X}), the following penalty formulation is proposed to achieve a continuous LS regularization:

(16)p={\displaystyle \frac{{\displaystyle {\int }_{{\mathrm{\Omega }}_{\mathrm{D}}}{(\varphi -\tilde{\varphi })}^{\mathrm{2}}}\mathrm{d}V}{{\displaystyle {\int }_{{\mathrm{\Omega }}_{\mathrm{D}}}{\varphi }_{\mathrm{Bnd}}^{\mathrm{2}}\mathrm{d}}V}}+{\displaystyle \frac{{\displaystyle {\int }_{{\mathrm{\Omega }}_{\mathrm{D}}}{\left|\nabla \varphi -\nabla \tilde{\varphi }\right|}^{\mathrm{2}}}\mathrm{d}V}{{\displaystyle {\int }_{\mathrm{\Omega }}\mathrm{d}}V}},

where {\varphi }_{\mathrm{Bnd}} denotes an upper (lower) bound for the target LSF. The penalty measures the difference between the design LSF \varphi (\mathit{X}) and a target LSF \tilde{\varphi }(\mathit{X}), as well as the difference in the spatial gradients. The target LSF \tilde{\varphi }(\mathit{X}) is constructed from the SDF {\varphi }_{\mathrm{SD}}(\mathit{X}). As the design LSF \varphi (\mathit{X}) converges to the target LSF \tilde{\varphi }(\mathit{X}), both contributions in Eq. (16) vanish. Penalizing the difference in the spatial gradients of the design and target LSFs, in addition to the difference in their values, increases the accuracy and avoids spatial oscillations. Both contributions are integrated over the entire design domain {\mathrm{\Omega }}_{\mathrm{D}}={\mathrm{\Omega }}_{\mathrm{I}}\cup {\mathrm{\Omega }}_{\mathrm{II}}. In the authors’ experience, equal weighting of the two contributions in Eq. (16) provides a good balance between matching the target LS value and avoiding oscillations. 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:

(17)\widehat{\varphi }(\mathit{X})=\{\begin{array}{ll}\begin{array}{l}{\varphi }_{\mathrm{Bnd}}\\ -{\varphi }_{\mathrm{Bnd}}\end{array}\hfill & \begin{array}{l}\mathrm{if}{\varphi }_{\mathrm{SD}}(\mathit{X})\ge {\varphi }_{\mathrm{Bnd}}\\ \mathrm{if}{\varphi }_{\mathrm{SD}}(\mathit{X})\le -{\varphi }_{\mathrm{Bnd}}\end{array}\hfill \\ {\varphi }_{\mathrm{SD}}(\mathit{X})\hfill & \mathrm{else}\hfill \end{array}.

The truncated target LSF \tilde{\varphi }(\mathit{X}) matches the SDF in the vicinity of the interface. Away from the interface, the lower bound -{\varphi }_{\mathrm{Bnd}} is matched in {\mathrm{\Omega }}_{\mathrm{I}} and the upper bound {\varphi }_{\mathrm{Bnd}} is matched in {\mathrm{\Omega }}_{\mathrm{II}}. To avoid the non-differentiability of the piecewise linear target LSF \tilde{\varphi }(\mathit{X}), a smooth target LSF \tilde{\varphi }(\mathit{X}) is used to approximate Eq. (17). This is achieved by the following sigmoid function:

(18)\tilde{\varphi }=\left({\displaystyle \frac{\mathrm{2}}{\mathrm{1}+{\mathrm{e}}^{\frac{-\mathrm{2}{\varphi }_{\mathrm{SD}}}{{\varphi }_{\mathrm{Bnd}}}}}}-\mathrm{1}\right){\varphi }_{\mathrm{Bnd}}.

A comparison between the piecewise linear target LSF and the smooth target LSF is illustrated in Fig. 3 for a one-dimensional interface configuration.

Fig.3 Piecewise and smooth approximation of the design LSF

Full size|PPT slide

The target LSF \tilde{\varphi }(\mathit{X}) depends implicitly on the geometry of the interface, defined by the zero iso-contour of the design LSF \varphi (\mathit{X}) that depends explicitly on the optimization variables. The implicit dependence of \tilde{\varphi }(\mathit{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_{\mathrm{3}} for the LS regularization term in the formulation of the objective function Eq. (3) is large and the implicit dependency of the target LSF \tilde{\varphi }(\mathit{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., ; 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 LSF \tilde{\varphi }(\mathit{X}) as a prescribed field which is reconstructed at every design iteration of the optimization process. Using this interpretation of \tilde{\varphi }(\mathit{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,\pm {\varphi }_{\mathrm{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 [⁸].

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\le \epsilon with \epsilon \ll \mathrm{1.0}. 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⁶ 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) [⁴⁸]. Bilinear four node quadrilateral elements are used in 2D, and trilinear eight node hexahedral 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) [⁴⁹], with no inner iterations.

Tab.1 Parameters used for all numerical examples with h denoting the element size

Parameter	Value
Weak boundary condition penalty Eq. (15)	{\gamma }_N=\mathrm{100}/h
Ghost penalty Eq. (13)	{\gamma }_G=\mathrm{0.001}
Perimeter penalty weight Eq. (3)	w_{\mathrm{2}}=\mathrm{0.01}
Lower bound of s	s_L=-\mathrm{3}h
Upper bound of s	s_U=+\mathrm{3}h
Target bound of LSF	{\varphi }_{\mathrm{Bnd}}=\mathrm{2}h
Filter radius used in 2D	r_f=\mathrm{1.6}h
Filter radius used in 3D	r_f=\mathrm{2.4}h

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 \mathrm{1}\times {\mathrm{10}}^{-\mathrm{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. [⁵⁰]. Selective structural springs [⁵¹] 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 {\varphi }_{\mathrm{Bnd}} for the target LSF \tilde{\varphi }(\mathit{X}) is within the range of values of the discretized design LSF \varphi (\mathit{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

(19)\mathit{R}={\displaystyle {\int }_{{\mathrm{\Omega }}_I}\delta \mathit{\epsilon }\hspace{0.17em}\mathit{\sigma }\mathrm{d}V}-{\displaystyle {\int }_{{\mathrm{\Gamma }}_T}\delta \mathit{u}\hspace{0.17em}\overline{\mathit{T}}\mathrm{d}A}=\mathbf{0},

with \mathit{u}=\overline{\mathit{u}} on {\Gamma }_u, where the displacement vector is denoted by \mathit{u}. The surface traction applied on {\Gamma }_T is \overline{\mathit{T}}. The infinitesimal strain tensor is defined by \mathit{\epsilon }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\nabla {\mathit{u}}^{\mathrm{T}}+\nabla \mathit{u}), and the Cauchy stress is \mathit{\sigma }=\mathit{D}\mathit{\epsilon }. The fourth-order constitutive 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 \lambda and \mu as follows:

(20)D_{ijkl}=\lambda {\delta }_{ij}{\delta }_{kl}+\mu ({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}),

where {\delta }_{ij} is the Kronecker delta. The Lamé constants can be expressed in terms of the Young’s modulus E and the Poisson’s ratio \nu as

(21)\lambda ={\displaystyle \frac{E\nu }{(\mathrm{1}+\nu )(\mathrm{1}-\mathrm{2}\nu )}},\mu ={\displaystyle \frac{E}{\mathrm{2}(\mathrm{1}+\nu )}}.

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. [⁵²].

Tab.2 Parameters used for the linear elastic design problems

Parameter	Value
Young’s modulus	E=\mathrm{2}\times {\mathrm{10}}^{\mathrm{3}}
Poisson’s ratio	\nu =\mathrm{0}
LS regularization weight	w_{\mathrm{3}}=\mathrm{0.01}
Element edge length	h=\mathrm{1.0}

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. [²] with the density method. The initial design of size \mathrm{80}\times \mathrm{40} with boundary conditions is shown in Fig. 4(a). Only one half of the design is analyzed and optimized, with weakly enforcing symmetry boundary conditions along the vertical center line. The top edge of the domain is clamped while a traction load of T_{X_{\mathrm{2}}}=-\mathrm{30} is applied in X_{\mathrm{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 a LS regularization penalty subject to a volume constraint of {\gamma }_V=\mathrm{0.16}. A relative GCMMA optimization step size of 0.1 and a time step size of \mathrm{\Delta }t=\mathrm{4.0} 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.

Fig.4 (a) Initial design of the 2D hanging bar problem and (b) final design with boundary conditions and dimensions

Full size|PPT slide

Fig.5 Evolution of (a) objective and (b) constraint with and without regularization for the 2D hanging bar

Full size|PPT slide

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₁ at the top of the design domain, i.e., X_{\mathrm{2}}=\mathrm{40.0}. When no LS regularization is used, slight 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.

Fig.6 (a) Comparison of the design LSFs with and without regularization at different snapshots during the optimization process (see Fig. 4), final design LSF (b) without regularization and (c) with regularization

Full size|PPT slide

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 be still observed even at convergence. In contrast, LS regularization achieves a design LSF with LS values at the target boundary values {\varphi }_{\mathrm{Bnd}} and a unit gradient in the vicinity of the interface. Due to the exclusion of the bottom center from the design domain, slight irregularities in the final iso-contour 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 modelled 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 domain. The dimensions of the design domain are \mathrm{80}\times \mathrm{40}\times \mathrm{80}. A volume constraint of {\gamma }_V=\mathrm{0.035} is enforced, and a relative GCMMA step size of 0.2 is used. A time step size of \mathrm{\Delta }t=\mathrm{6.0} is used in the HM.

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

Fig.7 (a) Initial design of the 3D hanging bar problem and (b) final design with boundary conditions and dimensions. The diagonal area highlighted in red represents the plane in which the design LSFs with and without regularization are being compared in Fig. 8

Full size|PPT slide

Fig.8 Comparison of the final design LSFs (a) without LS regularization and (b) with regularization

Full size|PPT slide

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

(22)\mathit{R}={\displaystyle {\int }_{{\mathrm{\Omega }}_I}\delta \mathit{F}\mathit{P}\mathrm{d}V}-{\displaystyle {\int }_{{\mathrm{\Gamma }}_T}\delta \mathit{u}\overline{\mathit{T}}\mathrm{d}A}=\mathrm{0},

where \mathit{F}=\partial \mathit{x}/\partial \mathit{X} is the deformation gradient tensor, and \mathit{x}=\mathit{u}+\mathit{X} describes the relationship between reference \mathit{X} and current \mathit{x} coordinates. The first Piola-Kirchhoff stress is denoted by \mathit{P}. A hyperelastic Saint-Venant Kirchhoff constitutive model for homogeneous, isotropic compressible materials is used, which is formulated as:

(23)\mathit{S}=\mathrm{2}\mu \mathit{E}+\lambda \mathrm{tr}(\mathit{E})\mathit{I},

</m:math>where \mathit{S} is the second Piola-Kirchhoff stress tensor and \mathit{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.

Tab.3 Parameters used for the hyperelastic design problems

Parameter	Value
Young’s modulus	E=\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{3}}
Poisson’s ratio	\nu =\mathrm{0.4}
LS regularization weight	w_{\mathrm{3}}=\mathrm{0.01}
Element edge length	h=\mathrm{1.0}

The Green-Lagrange strain tensor is defined as:

(24)\mathit{E}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathit{C}-\mathit{I}),

where the right Cauchy-Green deformation tensor Cis computed as \mathit{C}={\mathit{F}}^{\mathrm{T}}\mathit{F}. Finally, the first Piola-Kirchhoff stress is obtained from

(25)\mathit{P}=\mathit{F}\mathit{S}.

For more details on the formulation and verification of the nonlinear XFEM formulation used in this work, the interested reader is referred to Ref. [⁵³]. 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 \mathrm{240.0}\times \mathrm{40.0}. A traction load of T_{X_{\mathrm{2}}}=-\mathrm{10.0} is applied at the top center of the domain over a length of 3.0 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 modelled 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. [⁵⁴], a plane strain model is used in 2D. The objective of the optimization problem is to minimize the strain energy with a 10% penalty weight on the perimeter and a 1.0% penalty weight on the LS regularization. The optimization problem is subject to a volume constraint of {\gamma }_V=\mathrm{0.6}, and an optimization step size of 0.025 is used. A single time step of size \mathrm{\Delta }t=\mathrm{35.0} is used in the HM.

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

Full size|PPT slide

Fig.10 Comparison of the warped final design LSFs (a) without LS regularization, (b) with regularization, and (c) SDF of the 2D beam

Full size|PPT slide

Figure 9(b) shows a comparison of the zero LS iso-contour 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, \varphi (\mathit{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 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.

The influence of different weights w₃ 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_{\mathrm{3}}=\{\mathrm{0.01},\mathrm{0.05},\mathrm{0.1},\mathrm{0.5}\}. 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.

Fig.11 Evolution of (a) strain energy and (b) regularization penalty for different penalization weights

Full size|PPT slide

For a weighting parameter in the range of \mathrm{1.0}\%\le w_{\mathrm{3}}\le \mathrm{10.0}\%, 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.0%), the optimization problem changes noticeably and the physical performance of the optimized design is affected. The LS regularization term 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, a LS regularization weight up to 10.0% 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 LSF \tilde{\varphi }(\mathit{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 a LS regularization weight of w_{\mathrm{3}}=\mathrm{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_{\mathrm{3}}=\mathrm{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.

Fig.12 Spurious void inclusions within the structure as a result of including implicit design sensitivities of the target LSF with a regularization penalty weight of 10.0%

Full size|PPT slide

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 \mathrm{240.0}\times \mathrm{40.0}\times \mathrm{40.0}. An element edge length of h=\mathrm{2.0} is used for this example, along with a time step size of \mathrm{\Delta }t=\mathrm{51.0} in the HM. Analogous to the 2D configurations, a traction load of T_{X_{\mathrm{2}}}=-\mathrm{2.0} is applied within a circular region of radius 2.0 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 {\gamma }_V=\mathrm{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.

Fig.13 (a) Initial design of the 3D beam problem with boundary conditions and dimensions; (b) final design. The central area highlighted in red represents the plane in which the design LSFs are being compared in Fig. 14

Full size|PPT slide

Fig.14 Comparison of the design LSFs at the mid-plane of the beam (a) without and (b) with LS regularization

Full size|PPT slide

Tab.4 Parameters used for the fluid design problem

Parameter	Value
Reynolds number	\mathrm{Re}=\mathrm{66.0}
Fluid density	\rho =\mathrm{1.0}
LS regularization weight	w_{\mathrm{3}}=\mathrm{0.05}
Element edge length	h=\mathrm{0.25}

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 -{\varphi }_{\mathrm{Bnd}}=-\mathrm{2.0} 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. [⁵⁵].) or local mesh refinement would be required to properly regularize the design LSF.

Example for incompressible Navier-Stokes flow

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

(26)\mathit{R}={\displaystyle {\int }_{{\mathrm{\Omega }}_I}[\delta \mathit{v}\rho (\mathit{v}\nabla \mathit{v})+\mathit{\epsilon }(\nabla v)\mathit{\sigma }(\mathit{v},p)]\mathrm{d}V}+{\displaystyle {\int }_{{\mathrm{\Omega }}_I}\delta p\nabla \mathit{v}\mathrm{d}V},

where v is the velocity vector, p is the pressure, \rho is the density. The admissible test functions for velocity and pressure are denoted by \delta \mathit{v} and \delta p, respectively. The infinitesimal strain rate tensor \delta \mathit{v} is defined as:

(27)\epsilon (\mathit{v})={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\nabla {\mathit{v}}^{\mathrm{T}}+\nabla \mathit{v}).

The Cauchy stress tensor for Newtonian fluids is denoted by \mathit{\sigma }(\mathit{v},p) and is defined as:

(28)\mathit{\sigma }(v,p)=-p\mathit{I}+\mathrm{2}{\mu }_{\mathrm{D}}\mathit{\epsilon }(\mathit{v}),

where {\mu }_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. [⁴⁵]. The fluid domain boundary is decomposed into the fluid-void interface {\mathrm{\Gamma }}_{\mathrm{I},\mathrm{I}\mathrm{I}}, and Dirichlet and Neumann external boundaries, {\mathrm{\Gamma }}_u and {\mathrm{\Gamma }}_T, respectively. No-slip conditions are applied weakly at the fluid-void interface; the other boundary conditions are problem dependent and are specified below.

Fluid nozzle in 3D

An extension of the fluid nozzle problem studied by Refs. [⁵⁶–⁵⁸], 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 \mathrm{5.0}\times \mathrm{5.0}\times \mathrm{5.0} cube with a \mathrm{0.75}\times \mathrm{5.0}\times \mathrm{5.0} non-design domain downstream from the inlet face. A parabolic velocity profile with a maximum inlet velocity of 30.0 in X_{\mathrm{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_{\mathrm{1}}-X_{\mathrm{2}} and X_{\mathrm{1}}-X_{\mathrm{3}} symmetry planes.

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

Full size|PPT slide

Assuming a constant density of \rho =\mathrm{1.0} and a dynamic viscosity of {\mu }_D=\mathrm{1.0}, the flow conditions correspond to a Reynolds number of \mathrm{Re}=\mathrm{66.0} with the reference velocity being the average inlet velocity and the reference length being the edge length of the design domain, i.e., L_{\mathrm{Ref}}=\mathrm{5.0}. 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. [⁴⁵].

The objective of this nozzle design problem is the minimization of the total pressure drop between inlet and outlet along with a 1.0% perimeter penalty and a 5.0% LS regularization penalty. The optimization problem is subject to a {\gamma }_V=\mathrm{0.3} volume constraint on the fluid phase. A relative GCMMA step size of 0.01 and a single time step of \mathrm{\Delta }t=\mathrm{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 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.

Fig.16 Comparison of the design LSFs across the diagonal of the fluid nozzle final design (a) without and (b) with LS regularization

Full size|PPT slide

Fig.17 Evolution of (a) normalized objective and (b) LS regularization penalty for different number of time steps

Full size|PPT slide

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 }=\mathrm{1.0}. 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. [³³], 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.0%, 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 eliminates the need for reinitialization. Furthermore, the proposed method mitigates robustness issues observed with regularization schemes that solely operate on the value or the spatial gradients of the design LSF. The capabilities of the proposed method were demonstrated through numerical examples in 2D and 3D, including problems in linear and nonlinear elasticity and fluid mechanics.

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 [⁵⁹] or topological derivatives.

Acknowledgements

The first, second and fourth authors acknowledge the support of the United States National Science Foundation (CMMI-1463287). The third author acknowledge the support of the SUTD Digital Manufacturing and Design (DManD) Centre supported by the National Research Foundation of Singapore. The fourth author acknowledges the support of the Air Force Office of Scientific Research (Grant No. FA9550-16-1-0169) and from the Defense Advanced Research Projects Agency (DARPA) under the TRADES program (agreement HR0011-17-2-0022). The opinions and conclusions presented in this paper are those of the authors and do not necessarily reflect the views of the sponsoring organizations.

Open Access

This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

References

1 Bendsøe M P, Sigmund O. Topology Optimization. Berlin: Springer, 2004

2 Sigmund O, Maute K. Topology optimization approaches: A comparative review. Structural and Multidisciplinary Optimization, 2013, 48(6): 1031–1055

3 Deaton J D, Grandhi R V. A survey of structural and multidisciplinary continuum topology optimization: Post 2000. Structural and Multidisciplinary Optimization, 2014, 49(1): 1–38

4 Osher S, Sethian J A. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 1988, 79(1): 12–49

5 Osher S, Fedkiw R. Level Set Methods and Dynamic Implicit Surfaces. Vol. 153. New York: Springer, 2003

6 Osher S, Paragios N. Geometric Level Set Methods in Imaging, Vision, and Graphics. New York: Springer, 2003

7 Gibou F, Fedkiw R, Osher S. A review of level-set methods and some recent applications. Journal of Computational Physics, 2018, 353: 82–109

8 van Dijk N P, Maute K, Langelaar M, Level-set methods for structural topology optimization: A review. Structural and Multidisciplinary Optimization, 2013, 48(3): 437–472

9 Osher S J, Santosa F. Level set methods for optimization problems involving geometry and constraints I. Frequencies of a two-density inhomogeneous drum. Journal of Computational Physics, 2001, 171(1): 272–288

10 Allaire G, Jouve F, Toader A M. A level-set method for shape optimization. Comptes Rendus Mathématique, 2002, 334(12): 1125–1130 (in French)

11 Wang M Y, Wang X, Guo D. A level set method for structural topology optimization. Computer Methods in Applied Mechanics and Engineering, 2003, 192(1–2): 227–246

12 Allaire G, Jouve F, Toader A M. Structural optimization using sensitivity analysis and a level-set method. Journal of Computational Physics, 2004, 194(1): 363–393

13 Wang S, Wang M Y. Radial basis functions and level set method for structural topology optimization. International Journal for Numerical Methods in Engineering, 2006, 65(12): 2060–2090

14 Luo Z, Wang M Y, Wang S, A level set-based parameterization method for structural shape and topology optimization. International Journal for Numerical Methods in Engineering, 2008, 76(1): 1–26

15 Kreissl S, Pingen G, Maute K. An explicit level set approach for generalized shape optimization of fluids with the lattice Boltzmann method. International Journal for Numerical Methods in Fluids, 2011, 65(5): 496–519

16 van Dijk N P, Langelaar M, van Keulen F. Explicit level-set-based topology optimization using an exact Heaviside function and consistent sensitivity analysis. International Journal for Numerical Methods in Engineering, 2012, 91(1): 67–97

17 de Ruiter M J, van Keulen F. Topology optimization using a topology description function. Structural and Multidisciplinary Optimization, 2004, 26(6): 406–416

18 Haber R B, Jog C S, Bendsøe M P. A new approach to variable-topology shape design using a constraint on perimeter. Structural Optimization, 1996, 11(1): 1–12

19 Yamada T, Izui K, Nishiwaki S, A topology optimization method based on the level set method incorporating a fictitious interface energy. Computer Methods in Applied Mechanics and Engineering, 2010, 199(45–48): 2876–2891

20 Otomori M, Yamada T, Izui K, Level set-based topology optimisation of a compliant mechanism design using mathematical programming. Mechanical Science, 2011, 2(1): 91–98

21 Gomes J, Faugeras O. Reconciling distance functions and level sets. Journal of Visual Communication and Image Representation, 2000, 11(2): 209–223

22 Zhu B, Zhang X. A new level set method for topology optimization of distributed compliant mechanisms. International Journal for Numerical Methods in Engineering, 2012, 91(8): 843–871

23 Li C, Xu C, Gui C, Level set evolution without re-initialization: A new variational formulation. In: Proceedings of 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, (CVPR’05). San Diego: IEEE, 2005, Vol. I, 430–436

24 Coffin P, Maute K. Level set topology optimization of cooling and heating devices using a simplified convection model. Structural and Multidisciplinary Optimization, 2016, 53(5): 985–1003

25 Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow. Journal of Computational Physics, 1994, 114(1): 146–159

26 Osher S. Book review: Level set methods and fast marching methods: Evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science. Mathematics of Computation, 2001, 70(233): 449–451

27 Li C, Xu C, Gui C, Distance regularized level set evolution and its application to image segmentation. IEEE Transactions on Image Processing, 2010, 19(12): 3243–3254

28 Zhu B, Zhang X, Fatikow S. Structural topology and shape optimization using a level set method with distance-suppression scheme. Computer Methods in Applied Mechanics and Engineering, 2015, 283: 1214–1239

29 Jiang L, Chen S. Parametric structural shape & topology optimization with a variational distance-regularized level set method. Computer Methods in Applied Mechanics and Engineering, 2017, 321: 316–336

30 Burger M, Osher S J. A survey in mathematics for industry a survey on level set methods for inverse problems and optimal design. European Journal of Applied Mathematics, 2005, 16(2): 263–301

31 Hartmann D, Meinke M, Schröder W. The constrained reinitialization equation for level set methods. Journal of Computational Physics, 2010, 229(5): 1514–1535

32 Fu J, Li H, Xiao M, Topology optimization of shell-infill structures using a distance regularized parametric level-set method. Structural and Multidisciplinary Optimization, 2018, 1–14 (in press)

33 Crane K, Weischedel C, Wardetzky M. Geodesics in heat. ACM Transactions on Graphics, 2013, 32(5): 1–11

34 Crane B K, Weischedel C, Wardetzky M. The heat method for distance computation. Communications of the ACM, 2017, 60(11): 90–99

35 Kreissl S, Maute K. Levelset based fluid topology optimization using the extended finite element method. Structural and Multidisciplinary Optimization, 2012, 46(3): 311–326

36 Sethian J A. Fast marching methods. SIAM Review, 1999, 41(2): 199–235

37 Wong T, Leung S. A fast sweeping method for eikonal equations on implicit surfaces. Journal of Scientific Computing, 2016, 67(3): 837–859

38 Daux C, Moës N, Dolbow J, Arbitrary branched and intersecting cracks with the extended finite element method. International Journal for Numerical Methods in Engineering, 2000, 48(12): 1741–1760

39 Fries T P, Belytschko T. The extended/generalized finite element method: An overview of the method and its applications. International Journal for Numerical Methods in Engineering, 2010, 84(3): 253–304

40 Hansbo A, Hansbo P. A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Computer Methods in Applied Mechanics and Engineering, 2004, 193(33–35): 3523–3540

41 Makhija D, Maute K. Numerical instabilities in level set topology optimization with the extended finite element method. Structural and Multidisciplinary Optimization, 2014, 49(2): 185–197

42 Tran A B, Yvonnet J, He Q C, A multiple level set approach to prevent numerical artefacts in complex microstructures with nearby inclusions within XFEM. International Journal for Numerical Methods in Engineering, 2011, 85(11): 1436–1459

43 Burman E, Hansbo P. Fictitious domain methods using cut elements: III. A stabilized Nitsche method for Stokes’ problem. Mathematical Modelling and Numerical Analysis, 2014, 48(3): 859–874

44 Schott B, Rasthofer U, Gravemeier V, A face-oriented stabilized Nitsche-type extended variational multiscale method for incompressible two-phase flow. International Journal for Numerical Methods in Engineering, 2015, 104(7): 721–748

45 Villanueva C H, Maute K. CutFEM topology optimization of 3D laminar incompressible flow problems. Computer Methods in Applied Mechanics and Engineering, 2017, 320(Suppl C): 444–473

46 Burman E, Claus S, Hansbo P, CutFEM: Discretizing geometry and partial differential equations. International Journal for Numerical Methods in Engineering, 2015, 104(7): 472–501

47 Nitsche J A. On a variational principle for the solution of Dirichlet problems under the use of subspaces which are subject to no boundary conditions. Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg, 1971, 36(1): 9–15 (in German)

48 Amestoy P R, Guermouche A, L’Excellent J Y, Hybrid scheduling for the parallel solution of linear systems. Parallel Computing, 2006, 32(2): 136–156

49 Svanberg K. The method of moving asymptote—A new method for structural optimization. International Journal for Numerical Methods in Engineering, 1987, 24(2): 359–373

50 Sharma A, Villanueva H, Maute K. On shape sensitivities with Heaviside-enriched XFEM. Structural and Multidisciplinary Optimization, 2017, 55(2): 385–408

51 Geiss M J, Maute K. Topology optimization of active structures using a higher-order level-set-XFEM-density approach. In: Proceedings of 2018 Multidisciplinary Analysis and Optimization Conference, AIAA AVIATION Forum, (AIAA 2018-4053). Atlanta, 2018

52 Villanueva C H, Maute K. Density and level set-XFEM schemes for topology optimization of 3-D structures. Computational Mechanics, 2014, 54(1): 133–150

53 Geiss M J, Bodeti N, Weeger O, Combined level-set-XFEM-density topology optimization of 4D printed structures undergoing large deformation. Journal of Mechanical Design, 2018 (in press)

54 Sharma A, Maute K. Stress-based topology optimization using spatial gradient stabilized XFEM. Structural and Multidisciplinary Optimization, 2018, 57(1): 17–38

55 Guo X, Zhang W, Zhong W. Explicit feature control in structural topology optimization via level set method. Computer Methods in Applied Mechanics and Engineering, 2014, 272: 354–378

56 Kreissl S, Pingen G, Maute K. Topology optimization for unsteady flow. International Journal for Numerical Methods in Engineering, 2011, 87(13): 1229–1253

57 Makhija D, Maute K. Level set topology optimization of scalar transport problems. Structural and Multidisciplinary Optimization, 2015, 51(2): 267–285

58 Borrvall T, Petersson J. Topology optimization of fluids in Stokes flow. International Journal for Numerical Methods in Fluids, 2003, 41(1): 77–107

59 Dunning P D, Alicia Kim H. A new hole insertion method for level set based structural topology optimization. International Journal for Numerical Methods in Engineering, 2013, 93(1): 118–134