Hole Seeding in Level Set Topology Optimization via Density Fields

09/24/2019 ∙ by Jorge L. Barrera, et al. ∙ University of Colorado Boulder 0

Two approaches that use a density field for seeding holes in level set topology optimization are proposed. In these approaches, the level set field describes the material-void interface while the density field describes the material distribution within the material phase. Both fields are optimized simultaneously by coupling them through either a single abstract design variable field or a penalty term introduced into the objective function. These approaches eliminate drawbacks of level set topology optimization methods that rely on seeding the initial design domain with a large number of holes. Instead, the proposed approaches insert holes during the optimization process where beneficial. The dependency of the optimization results on the initial hole pattern is reduced, and the computational costs are lowered by keeping the number of elements intersected by the material interface at a minimum. In comparison to level set methods that use topological derivatives to seed small holes at distinct steps in the optimization process, the proposed approaches introduce holes continuously during the optimization process, with the hole size and shape being optimized for the particular design problem. The proposed approaches are studied using the extended finite element method for spatial discretization, and the solid isotropic material with penalization for material interpolation using fictitious densities. Their robustness with respect to algorithmic parameters, dependency on the density penalization, and performance are examined through 2D and 3D benchmark linear elastic numerical examples, and a geometrically complex mass minimization with stress constraint design problem.



There are no comments yet.


page 17

page 18

page 19

page 22

page 23

page 25

page 28

page 34

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Topology optimization (TO) via level set (LS) methods typically relies on shape sensitivities along material interfaces to evolve the design. Considering material-void problems, holes cannot nucleate, only merge, split or disappear. To facilitate topological changes in the optimization process, either the initial design domain is seeded with an array of holes or holes are introduced at distinct stages of the optimization process. In this paper, two alternative hole seeding approaches that nucleate holes during the design process are proposed.

Figure 1: Stages of the optimization process for a 2D Beam LS TO problem that minimizes the strain energy () using different initial holes seeding patterns.

In absence of any additional mechanism that introduces holes in the course of the optimization process, appropriate initial hole seeding is critical to avoid suboptimal designs (van Dijk et al. (2013)). In cases where the geometry and the associated mechanical properties of a particular hole pattern violate design constraints (e.g., mass, stress, stiffness), often the majority of holes merge in the early stages of the evolution of the design. Premature hole merging can restrict the final design to simplistic geometries and, in certain cases, trigger the appearance of large disconnected subdomains. The latter can result in ill-posed problems and may compromise the stability of the TO process. Even small variations in seeding parameters (i.e., number, size, shape or arrangement of holes) can degenerate the evolution of the design and severely affect its performance (Wang et al. (2007a)). An example of this scenario is shown in Fig. 1 for a strain energy minimization problem subject to a mass constraint (Bendsøe and Sigmund (2004)). Stages of the design process with different initial seeding patterns are shown to demonstrate that early merging of holes (bottom row) can result in final designs with poor performance ( difference).

In some cases, the shortcomings mentioned above can be mitigated by using a “large enough” number of holes as an initial guess, see Villanueva and Maute (2014). Excessive hole seeding, however, requires a fine mesh and a considerable number of design iterations to substantially change the shape of the initial hole pattern and alter conceptually the geometry. Furthermore, overseeding may produce more refined features without necessarily converging to one specific final design nor impacting performance significantly. This case is displayed in Fig. 2 where, despite noticeable differences in the final designs, similar performances are achieved for initial designs using 12, 24, 48, and 96 holes.

Figure 2: Initial designs with different number of holes seeded (left) and their corresponding final designs (right) for a 2D beam strain energy () minimization problem.

The aforementioned issues have motivated LS approaches capable of nucleating holes. These methods not only circumvent the computational burden present when using a large number of initial holes, but they also eliminate the difficulty of finding a suitable initial seeding pattern.

Originated from the bubble method (Eschenauer et al. (1994)), topological derivatives can be used to seed new holes during the optimization process (Allaire et al. (2005); Burger et al. (2004); Wang et al. (2004); Norato et al. (2007); Andreasen et al. (2019)). These derivatives evaluate the influence of introducing infinitesimal holes at any point in the design domain during the optimization problem. Typically, finite size holes are inserted at locations of minimal topological derivative at distinct steps in the optimization process. Strategies based on topological sensitivity formulations that nucleate holes in regions of low strain energy and stress, among others, have been proposed (Sethian and Wiegmann (2000); Belytschko et al. (2003); Sokolowski and Zochowski (1999); Park and Youn (2008). However, regardless of the underlying formulation, using topological derivatives as a guide for hole nucleation introduces user-defined parameters that control the number, frequency, and shape of the seeded holes. This can lead to adding arbirtrarily-shaped holes either excessively or scarcely and, as a result, compromise efficiency, robustness and performance (Allaire et al. (2004)). A systematic procedure for deciding if and when nucleating a hole is currently lacking.

Alternative hole seeding mechanisms that use topological derivatives in a non-classical sense, as well as schemes that do not rely on these derivatives, have also been developed. For example, a topological derivative field can be used to directly extract a geometry if interpreted as a level set field (LSF). In Suresh (2010, 2013); Suresh and Takalloozadeh (2013); Deng and Suresh (2015, 2016), holes are nucleated by redefining the isocontour that separates the material from the void subdomains. If an implicit LS TO framework is used, the solution of the Hamilton-Jacobi equation can be approximated to allow for the nucleation of holes, as demonstrated by Wang et al. (2007b). Alternatively, as shown by Dunning and Kim (2013), a second LSF can be used to construct a pseudo third dimension to nucleate holes in 2D problems.

The class of approaches explored in this paper uses density methods (Bendsøe and Sigmund (2004)) to nucleate holes in an informed manner. Since these methods can start from an unbiased design (i.e., homogeneous density distribution) and quickly generate the conceptual layout of an optimal design (Sigmund and Maute (2013)), shortcomings associated with initial hole seeding are avoided. In these approaches, regions of low density indicate that material is no longer needed and can be used to nucleate holes. Furthermore, sensitivities in the entire design domain characteristic of density-based TO are computed in addition to the shape sensitivities of the LS TO problem, achieving an improved efficiency compared to classical LS TO approaches.

Early attempts to use density methods to advance a well-defined material interface can be found in the context of (concurrent) shape and topology optimization methods. In Kumar and Gossard (1996), a “shape density” field was used to remove parts of the design domain for which densities were under a prescribed threshold. Maute and Ramm (1995); Bletzinger and Maute (1997) extensively explored sequentially alternating between shape and topology optimization using separate design and analysis meshes. In their work, post-processed densities define the material interface, which is described by splines (see also Maute and Ramm (1997); Maute et al. (1998)). Smoothing algorithms to realize a crisp material interface from intermediate densities, however, introduces inaccuracies in the description of the geometry.

More recently, density-based and LS-based TO methods have been combined for purposes other than hole nucleation. In Kang and Wang (2013), a classical solid isotropic material with penalization (SIMP) approach is used for material interpolation while a LS approach excludes void domains from the TO process. However, the LSF remains fixed during the design process. In an explicit LS TO setting developed by Jansen (2019), feature size control is achieved via geometric constraints on a density field. The combined LS-density approach in Geiss and Maute (2018); Geiss et al. (2019b) couples the LS and density fields to optimize the material distribution within the solid phase. However, the density field is not used for hole nucleation.

In this paper, two LS TO approaches that nucleate holes informed by a density field are proposed. The goal of these two approaches is to seed holes continuously in the optimization process. The location, shape and size of the holes evolve for the specific optimization problem at hand. In the first approach, a single abstract design variable field governs both the LS and density variables. In the second approach, two independent design variable fields are introduced to define the LS and density fields, respectively. The LSF is one-way coupled with the density field. Although feature size control based on the density field could be achieved, as demonstrated in Geiss (2019); Barrera et al. (2019), this paper focuses on hole nucleation only. Note that, similar to the topological derivatives field approach developed by Suresh (2010, 2013); Suresh and Takalloozadeh (2013); Deng and Suresh (2015, 2016), here an isocontour of a field indicates where a hole needs to be created. However, among other differences, in this paper this indicator density field evolves during the optimization process.

In the material-void problems considered in this study, a LSF is used to distinguish between phases through crisp, well-defined interfaces. This LSF is parametrized by an explicit LS method (van Dijk et al. (2013); Sigmund and Maute (2013)). In addition, the density field interpolates the material properties within the material phase using the SIMP scheme (Bendsøe and Sigmund (2004)). The weak form of the governing equations are discretized by the extended finite element method (XFEM, Belytschko et al. (2009)). Both approaches are studied for linear elastic problems in 2D and 3D, and considering different optimization formulations.

The remainder of the paper is organized as follows: Section 2 presents the basic components of the LS TO framework employed. Section 3 explains the two coupling schemes adopted. Section 4 summarizes the discretization/analysis method and describes the general optimization problem formulation. Numerical examples are provided in Section 5; and Section 6 concludes this paper with directions for future work.

2 Level Set-based Topology Optimization

2.1 Geometry description

The material layout composed of two phases in a design domain, , is described by the LSF:


where and are the material domains of phases and , respectively, such that . The interface, denoted by , corresponds to the zero LS isocontour, . The problems studied in this work consider solid-void domains where solid is assigned to phase I and void to phase II.

The LSF, bounded between lower, , and upper, , bounds, is a function of a filtered LS design variables field, , introduced below.

2.2 LS design variables

In the LS TO framework employed here, the geometry of the solid-void interface is defined by a vector of LS optimization variables,

. Here, an LS optimization variable is assigned to each node of a structured mesh, and represents the number of nodes in such mesh.

To increase numerical stability and enhance convergence of the optimization problem, the linear filter presented in Kreissl and Maute (2012) is applied to to obtain a vector of filtered LS coefficients, . In this distance-based filter, a filtered LS coefficient, , at node , is defined as a function of its neighboring LS optimization variables, , at nodes , using the following expression:


The number of nodes within the filter radius, is denoted by ; and is the Euclidean distance between nodes i and j.

The filtered LS coefficients, , are used to parametrize the LS design variable field, , following the expression:


where , with the Sobolev space denoted by . The field is interpolated by bi-linear and tri-linear shape functions, , on quadrilateral and hexahedral meshes in 2D and 3D, respectively. The LSF, , is defined as a function of the design variable field, , as formulated in Section 3 separately for each of the hole seeding approaches considered in this paper.

Unlike implicit approaches where a Hamilton-Jacobi-type equation needs to be solved, the LSF is defined as an explicit function of the LS design variables, . This provides the advantage of solving the optimization problem via a mathematical programming technique. The reader is referred to Sigmund and Maute (2013) and van Dijk et al. (2013) for more details.

2.3 LS regularization

To avoid spurious oscillations in the LSF, the regularization scheme of Geiss et al. (2019a) is adopted. This approach promotes a uniform spatial gradient of the LSF near the solid-void interface while converging to either a positive or negative target value away from the interface. A truncated signed distance field used as target field, , is constructed via the heat method (Crane et al. (2017)). This entails solving two linear finite element problems (i.e., a heat conduction problem and a Poisson’s problem) on the entire design domain. The target field is enforced through a penalty in the objective, , of the following form:


where is the difference between the upper, , and lower, , bounds in the target LSF. In Geiss et al. (2019a), the weights and were kept constant in the entire design domain. In the current work, to balance the influence of the regularization components in the vicinity and away from the material interface, these weights are defined as:


The weights control the mismatch between design and target LSF near and far from the material interface, respectively. The weights do the same for the LS spatial gradient component in Eq. 4. The parameter , defined as:


controls the region of influence of the regularization through , as seen in Fig. 3. Larger values decrease and thus, the weight of the regularization in the vicinity of the interface.

Figure 3: LS regularization penalty as function of the parameter in the vicinity of the interface for a mesh with element length .

3 Hole Nucleation via Density Methods

Density methods start typically from an unbiased, i.e., homogeneous, density distribution and quickly converge to the conceptual layout of the optimized design. This attractive feature is used here to instruct the LS optimization problem to remove material where it is not needed. To this end, we introduce a nodal fictitious density field, bounded between 0 and 1 in the LS TO problem to nucleate holes in a mechanically-informed manner.

Following the same scheme presented in the previous section for the LSF, is defined as an explicit function of a filtered density design variables field, ; i.e., . The field is discretized using Eq. 3 and a vector of filtered density coefficients, . Similarly, the linear filter in Eq. 2 is employed to define in terms of a vector of nodal density optimization variables, .

In this section we present two approaches to couple the LS and density design variables. In both cases, the couplings of the LS and density fields are constructed to achieve that (i) a positive LSF is associated with high densities to represent material, i.e., ; and (ii) a negative LSF is linked to densities close to zero to identify void, i.e., . In the first approach, both the LS and density fields, and , are both functions of a single set of optimization variables; i.e., . We identify this approach as Single-Field Coupling (SFC). The second approach interpolates and using two independent sets of optimization variables. This approach is termed Two-Field Coupling (TFC).

3.1 Single-field coupling strategy

In the first proposed, we assume a single vector of abstract design variables, , to compute the vectors of filtered coefficients, , using Eq. 2, and interpolate using Eq. 3. Here, the LSF is defined by the following linear relation:


Similarly, the density field is obtained by using:


In both Eqs. 7 and 8, is a threshold that alters the LS value used to differentiate between solid and void. In this work, is kept constant, i.e., . The scaling parameter in Eq. 7 takes into account the mesh size, and is set to , where is the element length. The density field described by Eq. 8 is only defined in the material subdomain, .

3.2 Two-field coupling strategy

The second coupling strategy uses two sets of independent LS and density optimization variables, and . Consequently, the optimization process in this case operates on optimization variables; i.e., . Here, the LS and density fields are assumed to be identical to the LS and density design variable fields, respectively; i.e., and .

In this scheme, holes are seeded in the course of the optimization process by promoting a negative LSF in regions of low density through a penalty term, , defined as:


This coupling penalty is active in a region bounded by the and thresholds, described below. High LS values in regions of low densities are penalized in a decreasing manner from 1 at to 0 for . To serve as a hole nucleation mechanism, it is intentionally formulated by a discontinuous function that has zero gradients with respect to the density optimization variables and exhibits non-zero gradients with respect to the LS optimization variables only below .

3.2.1 LS and density thresholds

The LS threshold, , deactivates the penalization for . This threshold is set to a value below zero, i.e., . Setting would result in issues associated with the robustness of the hole nucleation process. An explanation of these issues is provided in Section 3.2.3.

The threshold is a parameter that decreases from an initial value, , to zero during the optimization process. A decreasing is desired since intermediate densities at a later stage of the design process could promote the creation of spurious holes, see Section 3.4.1. Note that should be lower than the prescribed initial density, , in the first optimization iteration. Setting would activate the penalty in the entire domain and could result in moving the entire LSF from material to void.

Through a continuation scheme active in the first design iterations, the density threshold, , is computed at every continuation step using:


where is the design iteration index. As seen in Fig. 4, decreases to gradually inhibit the nucleation of holes. The coupling is inactive for . The continuation step size, , is assumed to be constant. The exponent controls the rate at which the shift increases in each continuation step, and is set to to attenuate changes in this threshold at early stages of the continuation scheme.

Figure 4: Continuation scheme in Eq. 10 for the density threshold of the TFC approach.

3.2.2 Smooth coupling formulation

To obtain a differentiable penalty formulation with respect to the LSF, the expression in Eq. 9 is approximated by a smooth function as follows:


The nondimensional parameter smoothes the transition of the penalty, and is typically between . An illustration of this smooth two-field density-LS coupling scheme is presented in Fig. 5.

Figure 5: Landscape of the smooth coupling penalization in Eq. 11 for the TFC approach.

The TFC strategy introduced here represents an improvement over the definition of the coupling penalty in Geiss (2019). Numerical studies showed that the approach in Geiss (2019) is prone to oscillations of the LSF in low density regions. The spurious interaction between density and LS fields is avoided by the discontinuous formulation of the penalty term in Eq. 11. Since the gradients of the penalty term with respect to the density optimization variables are zero, this formulation does not provide information that would promote the removal of a hole through increasing the local densities. However, a hole can still be removed through shape changes. Thus, the coupling penalty formulation in Eq. 11 should be considered a hole nucleation scheme.

Note that the choice of variables in this coupling strategy is not unique. Either the sets of LS and density optimization variables, or any choice of filtered and/or projected quantities derived from them could be discretized and coupled. Our choice of using and is motivated by an increased in resolution and design freedom (compared to elemental variables), as well as smoothness in the optimization problem.

3.2.3 Hole nucleation process

Figure 6 depicts the effect of the TFC penalty in a simple 1D example. The density, LS, and coupling penalty fields, i.e., , , and , are plotted over at different stages of the nucleation process of a hole. Initially, in Fig. 6(a), the coupling penalty is inactive in the material domain, i.e., and everywhere. As the local density value decreases, the coupling penalty is activated along defined by the region of the density field smaller than . As a result, the LSF in is lowered as seen in Fig. 6(b). The nonlinear distribution of the penalty (see Eq. 11) is maximum at , and its effect diminishes as the LSF approaches . Eventually, in Fig. 6(c), a hole is nucleated by the coupling penalty, which continues decreasing the LSF (to a lesser extent) while it is above . Note that, as mentioned in Section 3.2.1, using can compromise the hole nucleation process. Setting may introduce a robustness issue since the LSF could oscillate between positive and negative values, nucleating and removing small holes intermittently. If , the coupling terminates prematurely and the LSF is pushed back into the material domain before crossing the zero isocontour, preventing a hole from being nucleated. This is a consequence of the LS regularization scheme detailed in Section 2.3. Hence, sufficiently away from the zero isocontour is used in this work. In Fig. 6(d), the LSF below along is no longer affected by the penalty, but continues decreasing until it reaches due to the effect of the LS regularization. At this stage, the density field is converged and the penalty is acting only in a small portion of the domain contained under and above (i.e., in ). Finally, Fig. 6(e) shows the regularized LSF with a hole nucleated in a region of low density. The penalization is no longer active since the process has been completed. A study of the influence of the two thresholds, and is presented in Section 5.1.3.

Figure 6: Evolution of the density, LS, and penalty fields at all stages of the nucleation process via the TFC approach.

3.3 Single versus two-field coupling

The two coupling strategies presented above enable hole seeding during the optimization process. In both cases, the design process consists of a handover from a pure density problem that does not need hole seeding but can produce intermediate density at convergence to a pure LS problem that relies on hole nucleation but features a crisply defined material interface. In terms of the gradients, only design sensitivities with respect to the density variables exist everywhere within the solid domain until holes are nucleated. Once that happens, non-zero LS shape sensitivities are also present in the vicinity of the material interface. Non-zero density sensitivities in an LS TO approach not only accelerate convergence, but also facilitate convergence to (local) minima with satisfactory performance and mitigate issues associated with initial seeding (see Fig. 1). Note that non-zero sensitivities in the entire design domain can also exist as a result of the LS regularization method used. The method adopted here is one of such cases (see Section 2.3).

In both hole nucleation schemes, the computational costs associated with physical analyses are reduced. The mesh resolution is no longer dominated by the need to resolve an initial hole seeding. Hence, the proposed method might allow using coarser meshes depending on the feature sizes in the optimized design. In the early stages of the optimization process, complex intersection configurations due to user-defined initial hole seeding are avoided. This reduces the number of intersected elements and thus, the computational cost when using immersed boundary techniques for the physical analyses. Furthermore, the number of total design iterations might decrease due to faster convergence. However, in contrast to evaluating LS shape sensitivities by only considering intersected elements, the contributions of all elements to the sensitivity equations within the solid domain need to be computed. The increased cost of the TFC approach over the SFC formulation associated with the increased number of optimization variables is insignificant because of the class of optimization algorithms and the type of sensitivity analysis used in this manuscript; see Section 5.

3.4 Material interpolation via augmented SIMP

Transition regions between high and low density material are commonly found in density-based methods. For example, SIMP-like approaches typically require a suitable formulation of the optimization problem and large penalties to converge to [0-1] material distributions. However, a penalization effect on intermediate densities, and thus their complete removal from the design domain, cannot be guaranteed (Sigmund and Maute (2013)).

3.4.1 Shifted density field

To avoid intermediate densities, we introduce a density shift, , to compute a shifted density field, , using the following expression:


The parameter increases from an initial value, , to 1 during the optimization process. It has the effect of shifting the minimum density in the material domain, as shown in Fig. 7 for the SFC approach. The larger , the more restricted the density distribution becomes. If , a uniform density field is obtained.

Figure 7: Effect of density shift, , in the SFC approach.

By shifting the densities, not only a constant density field is achieved in the final design, but ill-conditioning of the discretized governing equations due to large differences in material properties is alleviated. The latter is achieved because areas/volumes of low density and low material stiffness are removed by controlling the minimum density in the solid domain. To attain well-conditioned systems, unless a smaller value is required to satisfy a mass or volume constraint. Moreover, the density shift boosts the LS sensitivities of new holes since, by advancing the shift, the problem transitions from a pure SIMP to a pure LS formulation. Thus, the density sensitivities decrease, and the LS shape sensitivities increase.

A continuation scheme similar to the one used for the density threshold (see Eq. 10), is employed to update the density shift:


Eq. 13 ensures the elimination of intermediate densities for , see Fig. 8. Every iterations, the density shift is updated from an initial value, , to 1. Numerical studies suggest that setting the exponent , together with an appropriate number of continuations steps (i.e., ) sufficiently mitigate the effect of jumps in the densities and promote a smooth evolution of the optimization problem.

Figure 8: Continuation scheme for density shift, .

3.4.2 Material interpolation

The classical SIMP approach, as presented by Bendsøe and Sigmund (2004), is employed to relate the shifted fictitious density field, , to physical material properties. However, alternatives such as the Rational Approximation of Material Properties (RAMP) method (Stolpe and Svanberg (2001)) are equally suitable. The material density, , is interpolated using:


and the Young’s modulus, , is computed using the following power law:


The properties of the bulk material density and Young’s modulus are denoted by and . The Young’s modulus is penalized by the SIMP exponent, denoted by . A constant density within an element is interpolated by averaging the shifted nodal densities in such element.

While material properties are interpolated using a shifted filtered density field, , the hole seeding strategies presented earlier in this section couple the unshifted filtered density field, , to the LSF, . In the TFC approach, using the shifted densities in the formulation of the penalty term would terminate the coupling too early since hole nucleation capabilities are lost once the density shift exceeds the density threshold, i.e., . Also, later in the optimization process, does not imply . Hence, intermediate densities might still be present in , without affecting the interpolation of the material properties and thus, the optimization problem. For this reason, the density threshold is decreased gradually approaching zero (see Section 3.2.1). Otherwise, spurious holes might be nucleated in regions of high shifted densities in designs that are almost converged. A motivation for using the augmented material interpolation scheme presented in this section, together with an analysis of its influence, is discussed in Section 5.1.1.

4 Optimization Framework

The physical responses of the systems are predicted by the XFEM in this paper. Since LS TO is not restricted to a particular immersed boundary technique, alternative approaches such as CutFEM (Burman et al. (2015, 2019)) or HIFEM (Soghrati and Barrera (2016)) could also be used for physical and sensitivity analyses.

A generalized Heaviside enrichment strategy (Makhija and Maute (2014); Terada et al. (2003); Tran et al. (2011)) is employed to avoid artificial coupling of disconnected material. The response is consistently interpolated in elemental subdomains with the same phase. The discretized state variable field, , at node of a two material problem (, ) is approximated as:


where the Heaviside function, , is determined by the LSF as:


The maximum number of enrichment levels is denoted by , is the nodal shape function and is the Kronecker delta which selects the active enrichment level for node . ensures that displacements at node

are only interpolated by a single set of degrees of freedom (DOFs),

, such that the partition of unity principle is satisfied (Babuška and Melenk (1997)). The number of nodes per element is denoted by . For more details about the generalized Heaviside enrichment strategy used in this work, the reader is referred to Makhija and Maute (2014). The described XFEM framework has been successfully applied to various TO multiphysics problems (e.g., in Maute et al. (2015); Coffin and Maute (2016); Villanueva and Maute (2017); Behrou et al. (2017); Pizzolato et al. (2017)).

4.1 Governing equations

The following augmented residual equation with stabilization terms is considered in this work:


in which the weak form of the linear elastic governing equation, , is defined as:


The displacement field and the test function are denoted by and , respectively. Traction forces, , are applied along the boundary

. The material tensor for isotropic linear elasticity,

, together with the infinitesimal strain tensor, , define the Cauchy stress .

The , , and terms correspond to the weakly enforced essential boundary and interface conditions, face-oriented ghost stabilization, and selective structural springs, respectively.

4.1.1 Weak enforcement of boundary and interface conditions

The unsymmetric version of Nitsche’s method (Nitsche (1971)) is employed to weakly enforce Dirichlet boundary and interface conditions (Burman and Hansbo (2012)). These conditions are applied using:


The outward normal vector on a domain or interface boundary is denoted by . The jump operator with respect to a target state variable is defined as:


The integrals in (20) correspond to the standard consistency, adjoint consistency, and the Nitsche penalty terms, respectively. The third term is scaled by the Young’s modulus, , the element size, , and the penalty factor . The latter provides additional control over the accuracy at which a boundary condition (BC) is enforced.

4.1.2 Face-oriented ghost stabilization

Face-oriented ghost penalization, as proposed by Schott and Wall (2014) and Burman et al. (2015), is used in the vicinity of the interface. This stabilization technique prevents numerical instabilities due to vanishing zones of influence of certain DOFs occurring when the material interface moves too close to a node. Independent of the intersection configuration, ill-conditioning is mitigated by applying the following virtual work-based formulation (Geiss and Maute (2018); Geiss (2019)):


Element faces in the vicinity of the material interface for which at least one of the two adjacent elements is intersected are included in (Villanueva and Maute (2017)). Outward facing normals of these shared element faces are denoted by . The influence of the ghost penalty term presented above is controlled by the penalty factor . The virtual work-based formulation is adopted instead of the one proposed by Burman and Hansbo (2014) since it allows for different material properties in adjacent intersected elements. In this work, the ghost penalty term is computed based on elementally constant material properties that are evaluated at element centroids (un-intersected elements) or sub-phase centroids (intersected elements).

4.1.3 Selective structural springs

To suppress rigid body motion associated with disconnected material domains that may emerge and develop, selective structural springs as proposed in Villanueva and Maute (2017) are added to the governing equations. An additional stiffness is applied only to solid disconnected subdomains via the following residual component:


The parameter denotes the spring stiffness scaling and is non-zero only for the free-floating pieces of material. An auxiliary indicator field obtained by solving a diffusion problem is employed to identify such disconnected subdomains. More details of the application of this approach to structural problems can be found in Geiss and Maute (2018); Geiss et al. (2019b).

4.1.4 Evaluation of stresses

In the final example in Section 5, a gradient stabilized scalar stress field, , is post-processed via the XFEM informed smoothing procedure described in Sharma and Maute (2018). The additional set of state variables is computed by solving the residual equation:


The field and the parameter represent the von Mises stress and the ghost penalty weight, respectively. Overestimation of stresses is prevented by penalizing the jump in spatial stress gradients across elemental faces (second term in Eq. 24).

4.2 General optimization problem formulation

In this paper the following formulation of the optimization problem is considered:


The objective is minimized over the vector of admissible optimization variables, , defined in Sections 3.1 and 3.2 for the SFC and TFC, respectively; and the vector of state variables, , with , and being the number of state variables.

The first component of the objective represents the quantity of interest to be minimized, (e.g., strain energy, mass). The second term is the normalized perimeter control penalty,


which prevents the emergence of irregular geometric features. corresponds to the perimeter of the design domain . The LS regularization penalty, , is included in the objective to avoid spurious oscillations in the LSF (see Section 2.3). Note that promotes a positive LSF while the coupling penalty attempts to lower it to a negative value to nucleate a hole. By reducing the weights and in Eq. 5, these competing effects are mitigated. This flexibility is exploited in the third example in Section 5. The final component denotes the normalized two-field coupling penalty,


with , as formulated in Eq. 11. This coupling penalty can also be defined as a constraint. We choose to add it to the objective as it favors a looser coupling of the density and LS fields. In our experience, this approach is sufficient to nucleate holes and prevents over-constraining the design.

The weights , , , and are chosen such that all the penalty contributions in the objective are significantly lower () than . Since a strong perimeter control penalty might prevent the nucleation of small holes (Wang et al. (2007a)), its contribution is kept below by manipulating the weight. Both, a constant low weight and gradual increase of the perimeter contribution within the continuation scheme (i.e., ), are considered in the numerical examples. In addition, for the remainder of the optimization process (i.e., ), the perimeter penalty weight is increased to promote a smooth geometry in the final design. For the SFC approach, is set to 0.0 since a coupling penalty is not required.

The design needs to satisfy a set of problem dependent inequality constraints,

(e.g., target mass, maximum allowable stress, maximum eigenvalue). The constraints are defined for each problem studied in Section


5 Numerical Examples

The proposed SFC and TFC approaches are studied in this section using single material, solid-void linear elastic problems in 2D and 3D. Algorithmic parameter dependencies and performance are investigated with a 2D structural plane stress problem under uniform pressure. A beam problem in 3D is used to assess the influence of the SIMP penalization on the robustness of the coupling strategies. And finally, the overall behavior of these approaches is examined with a geometrically complex engineering problem considering stress constraints.

The optimization problems are solved using the Globally Convergent Method of Moving Asymptotes (GCMMA, Svanberg (2002)) with no inner iterations. The adjoint method, as detailed in Sharma et al. (2017), is used for the sensitivity analysis. The relative change in objective between two consecutive design iterations is less or equal to at the end of all the problems shown in this section. Relevant optimization parameters common to all problems presented in this work are summarized in Table 1. The weights for the LS regularization penalty in Examples 1 and 2 are , and for the third problem, , .

Parameter Value
LSF upper bound
LSF lower bound
Initial constant LSF
LSF regularization control
Filter radius in 2D
Filter radius in 3D
Nitsche penalty factor
Ghost penalty factor
Spring stiffness factor
Stress ghost penalty factor (Ex. 3)
Table 1: Parameters common to all numerical examples function of the element size .

At the beginning of the optimization process, for both the SFC and TFC approaches, a uniform LSF of is prescribed such that the entire design domain is filled with homogeneous porous material. In all examples, the continuation scheme described in Section 3.4 is used to update the shift in the densities unless specified otherwise; see Fig. 8. Furthermore, when using the TFC approach, the density threshold is decreased using the continuation scheme detailed in Section 3.2; see Fig. 4.

The governing equations are discretized using the XFEM approach outlined in Section 4. The parameters used for interpolation of material properties in Examples 1 and 2 are shown in Table 2 in self-consistent units, and for the last example in Table 5 in SI units. The example problems consist of a one-way coupled set of governing equations; i.e., the diffusion problem describing the auxiliary indicator field, the stabilized linear elasticity equations (Eq. 18), the stabilized stress projection (Eq. 24), and the heat method equations. Taking advantage of this simple coupling scheme, a Block Gauss-Seidel approach (Elfving (1980)) is employed for both the physical and the adjoint sensitivity analyses. The systems of equations of the first and second examples were solved using the Multifrontal Massively Parallel Solver (MUMPS, Amestoy et al. (2006)). The third example was solved iteratively via Trilinos algebraic multigrid solver (Heroux et al. (2005)) for the selective springs heat problem, the structural response, and the stresses; and ILU preconditioner for the LS regularization heat method problem.

Property Value
Initial fictitious density ()
SIMP exponent = 2.0
Young’s Modulus () x
Young’s Modulus () x
Poisson Ratio ( and ) = 0.4
Material Density ()
Material Density ()
Table 2: Material properties interpolation parameters for numerical examples 1 and 2.

5.1 Example 1: structure subject to uniform pressure

A 2D structure subject to uniform pressure load is studied to demonstrate the hole nucleation capabilities of the proposed approaches. The problem setup with loads and BCs is shown in Fig. 9. A vertical traction load of is applied over the top boundary of the design domain. The bottom left corner is clamped by prescribing weakly zero displacement in the and directions. Design symmetry about the axis is assumed in the middle of the design domain (i.e., zero displacement in the direction along symmetry plane). Solid material is prescribed and excluded from the design domain in the vicinity of the Dirichlet BCs and where the traction load is applied, as highlighted in dark grey color in Fig. 9. The domain of size is discretized by a structured mesh with uniform element size .

Figure 9: Problem setup of the 2D structural design problem. Due to symmetry, only half of the domain is simulated.

The compliance minimization with mass constraint optimization problem reads:


The strain energy of the initial design is denoted by , and the objective weights for the TFC setup are =[ 0.93, 0.01, 0.05, 0.01]. In the SFC setup, . A mass constraint of is enforced.

The default continuation parameters for updating the density shift (in both approaches) and the density threshold (in the TFC approach) are summarized in Table 3, and used for all example configurations unless specified otherwise. Both and are updated every optimization iterations in a total of eight continuation steps. No density shift is applied at the beginning of the optimization process. All simulations start with a uniform density of , except the portions excluded from the design domain, which are prescribed a fixed density of 1.0. Note that the SIMP exponent is set to 2.0 to mitigate mesh dependency, avoid premature convergence, and allow for nucleating many holes. An analysis of the influence of this parameter is provided in Example 2.

Parameter Value
Continuation step size
Number of design iterations in continuation
Maximum number of design iterations
Initial density shift
Initial density threshold
Continuation density shift exponent
Continuation density threshold exponent
Table 3: Default continuation parameters for the 2D structural problem.

5.1.1 Influence of density shift

In this first study, the effect of the density shift described in Section 3.4 is investigated. The final designs for both the single and two-field coupling approaches with and without density shift are shown in Fig. 10. Only the material phase, , is visible; and it is colored by the shifted filtered densities, . The void regions are identified by dotted patterns with white background. The SFC and TFC results are shown in rows one and two, respectively, without and with a density shift (left and right columns). The insets highlight regions of the final designs where significant changes in the density distributions are observed.

On the left side of Fig. 10 it can be seen that both methods are unable to fully remove intermediate densities in absence of a density shift. In the SFC approach, since both the LSFs and fictitious densities are function of the same abstract design variables (Eqs. 7 and 8), intermediate densities exist in the vicinity of the material interface by construction. In the TFC, despite the density field can evolve independently of the LSF, grey areas are also observed.

The results in the right column of Fig. 10 show that the gradual shift scheme adopted guarantees a uniform density distribution of 1 in the material subdomain. Moreover, even though the topologies being noticeably different, the performance of the final design is not compromised by shifting the densities. On the contrary, a considerable improvement is observed in the SFC results. Note that as an alternative to improve the convergence to [0-1] designs, the SIMP penalization could be increased and combined with a projection scheme (e.g., see Wang et al. (2011)). However, those approaches do not guarantee that the solid phase consists only of bulk material in the optimized design. See Example 2 for a more detailed discussion on this matter.

Figure 10: Final designs of the 2D structural problem for both SFC and TFC approaches with and without the augmented SIMP material interpolation.

Initial density shift

To better understand the influence of an initial density shift in the hole nucleation process, a comparison of the design evolutions with is presented. All parameters but the initial density shift are the same as before. Given that the numerical examples in this study show similar findings for both the SFC and TFC approaches, we only show results for the SFC approach.

Figure 11: Snapshots of the evolution of the 2D structural design using the SFC approach from left to right are shown with different initial density shift per row.

Snapshots of the optimization process at = [50, 100, 150, 500] for = [0.1, 0.2, 0.3] are provided in Fig. 11. A fast and clean hole nucleation process is observed in all cases. Note that densities in the void phase do not contribute to neither the stiffness nor the mass of the design. Therefore, sensitivities associated with these variables only exist in the material domain (). Material is removed from the design domain by either nucleating a hole due to low densities or shrinking/expanding the material interface front. A larger initial density shift promotes the generation of more holes at an early stage, as seen in the first column of Fig. 11. Nucleating less holes initially (guided by the density sensitivities) is, however, balanced by the moving material interface (result of shape sensitivities of the LS problem). Eventually, similar topologies at a later stage of the optimization process are observed. Although most holes are nucleated within the first third of the optimization process, the hole nucleation process is active, and holes are nucleated until the density shift enforces the maximum density at .

The performance of the final design is not compromised by introducing a density shift. As is reported on the last column of Fig. 11, the variation of the final strain energy is less than 1%. Considering the results of the previous subsection where (see Fig. 10), there is a small improvement on performance with an increasing shift.

5.1.2 Influence of continuation strategy

The final designs of different continuation settings for the density shift using the SFC approach are shown in Fig. 12. A total of six (left) and eight (right) continuation steps are used in optimization problems with (a) 400, (b) 300, and (c) 200 design iterations. The continuation scheme parameters for each case are provided on top of the final designs. These converged designs demonstrate that the proposed approaches provide the flexibility of choosing different continuation step sizes, , and total number of continuation iterations, , to accelerate the design process. Similar performances are achieved for all configurations. However, reducing the continuation step size decreases the number of holes nucleated. An excessively small continuation step size may cause the optimization process to prematurely converge to a suboptimal design. Fig.13 shows the evolution of the objective and constraints of the designs in Fig. 12, separated by total number of design iterations. Here, represent the maximum number of design iterations allowed. It can be seen that a larger number of continuation steps smoothes the transition of the density shift at the cost of more frequent modifications to the optimization problem.

Numerical experiments with the problems studied in this paper suggest that using five or more continuation steps of provides an adequate seeding capabilities and circumvents large jumps in the densities updates.

Figure 12: Final designs of 2D structural problem for multiple continuation step sizes and number of design iterations with a continuation strategy.
Figure 13: Evolution of objective and constraint of the 2D structural problem for two continuation step sizes and multiple total number of design iterations.

5.1.3 TFC optimization problem characteristics

Fig. 14 shows (a) the evolution of the objective and mass constraint, (b) the objective components contributions, and (c) the LS and density sensitivity contributions to the optimization problem for the TFC results with density shift of the 2D structural problem in Fig. 10. Small oscillations are observed early on in the design process due to holes being nucleated; see Fig. 14(a). Since this occurs in low density regions, the nucleation process has a reduced effect on the objective and thus, the observed oscillations are small and quickly vanish. A smooth design evolution is achieved in a later stage. This is in contrast to using alternatives like topological derivatives for which, in our experience, the size and number of holes nucleated may produce large jumps in the optimization process. A slight increase in the objective at originates from increasing the weight of the perimeter control penalty. However, its impact on the strain energy convergence is negligible. Small jumps in the mass constraint, , at every continuation step update are consequence of sudden increments in mass due to the density shift.

Figure 14: Evolution of (a) objective and mass constraint, (b) logarithm of contributions of the components of the objectives; and (b) implicit sensitivities contributions (logarithms of norms) of the density and LS design variables for the 2D structural problem using the TFC approach.

Objective components contributions

The logarithm of the contributions of each objective component (i.e., , , , and ; see Eq. 28), are provided in Fig. 14(b). The optimization problem is at all times completely dominated by the strain energy. The coupling, LS regularization and perimeter control penalties are on average two orders of magnitude lower. As a result, the evolution of the objective is predominantly smooth despite considerable oscillations and periods of inactivity of ; see Fig. 14(b). The intermittent non-zero TFC penalty contribution is a direct consequence of its definition in Eq. 11. Intermediate densities larger than may shift up or down, when a hole is nucleated or due to the moving interface front. Sudden increments in the coupling penalty represent an intermediate density dropping below .

LS and Density sensitivities

The interplay of the LS and density variables in the TFC approach is examined with Fig. 14(c). The logarithm of the norms of the sensitivities with respect to the LS and density variables show that initially the density field drives the optimization process. Then, gradually with the introduction of holes, the shape sensitivities take over. The influence of the density sensitivities vanishes once the shifted densities have reached the maximum density at . A pure LS problem with shape sensitivities only is recovered for .

5.1.4 TFC thresholds

As explained in Section 3.2, the TFC coupling penalty is bounded by the LS () and the density () thresholds. To avoid robustness issues in the hole nucleation process, is set to a number below the zero isocontour ( in this example). As long as this condition is satisfied, no significant changes are observed in the hole nucleation process. Similarly, the density threshold is constrained by the uniform density field used at the beginning of the optimization process. However, unlike the LS threshold, different settings of the density threshold can have a considerable impact in the initial stages of the hole nucleation process, as is demonstrated below.

Snapshots of the evolution of the 2D structural design problem are presented in Fig. 15 to assess the influence of setting the initial density threshold closer to the initial constant density field. The designs at =[50, 100, 150, 500] for =[, , ] separated by rows, show significant variations in the hole nucleation process. As expected, setting closer to promotes the nucleation of holes earlier given that the penalty is active earlier in the design process. The final designs for different values exhibit slightly different topologies but similar performances. Also, as seen for the SFC approach previously, there is an interplay between the densities and the LS interface front until the density shift removes all intermediate densities. As long as intermediate densities exist, a hole can be nucleated if it is advantageous.

Figure 15: Snapshots of the evolution of the 2D structural problem using the TFC approach from left to right are shown using different initial density threshold per row.

Overall, similar performances are achieved using a wide range of choices for parameters in both the single- and two-field strategies, demonstrating the robustness of these approaches with respect to algorithmic parameters.

5.2 Example 2: beam

A 3D beam problem is used to investigate the effect of the SIMP penalizations in the density field, as well as the behavior of both approaches compared to LS-based with initial hole seeding and density-based TO. The design domain of size is simply supported weakly at all four corners at the bottom face. A traction load is applied over the center of the top face. Due to symmetry, only one quarter of the full domain is analyzed and optimized. The problem setup is shown in Fig. 16.

Figure 16: Problem setup for the classical Beam problem in 3D. Due to symmetry, only one quarter of the domain is simulated.

The optimization problem formulation remains the same as in the previous example (see Eq. 28). The initial objective weights are = [0.92, 0.001, 0.01, 0.05]. The perimeter weight, , is updated at every continuation step from 0.001 to 0.01, using Eq. 13 with an exponent of 3.0. A mass constraint of is enforced.

The density field is shifted using the augmented SIMP continuation scheme discussed in Section 3.4; and the density threshold for the TFC approach is gradually reduced to zero using Eq. 10. The continuation parameters are listed in Table 4. The design domain is initialized with a uniform density of everywhere but along the prescribed BCs, where the density is set to 1.0 for the duration of the optimization process.

Parameter Value
Continuation step size
Number of design iterations in continuation
Maximum number of design iterations
Continuation density threshold exponent
Continuation density shift exponent
Initial density threshold
Initial density shift
Table 4: Continuation parameters for the beam problem.

5.2.1 SIMP density penalization

Large exponents (i.e., ) reduce intermediate densities at the cost of more nonlinear optimization problems and poorer conditioned finite element problems. Contrary to these methods, to converge to a [0-1] density field we employ the density shift scheme presented in Section 3.4. To examine the reliance of the proposed approaches on the penalization of the densities for the material description (see Eq. 15), different exponents in the power law are considered.

Figure 17: Final designs for both the SFC (left), and TFC (right) approaches using SIMP exponents of .
Figure 18: Evolution of the strain energy in the 3D beam problem using the (a) SFC and (b) TFC approaches. Snapshots of designs before and after updating continuation parameters are colored by the Young’s modulus.

The final designs of the beam problem using the single and two-field strategies with penalizations of =[1.0, 2.0, 3.0] are shown in Fig. 17. Their strain energy is also reported to evaluate structural performance. All designs converge to a [0-1] density distribution as a consequence of the density shift. Note that, although not shown, intermediate densities were also observed without the shifting scheme, in accordance with the findings in Example 1 for the 2D problem. For all values, similar truss-like structures with partial shear-webs are obtained. When using the TFC approach, the performances for the designs does not differ significantly with changing . For the SFC approach, however, the strain energy of the optimized design increases noticeably with increasing value. Due to the tight coupling for density and LS fields, and the penalizations of intermediate densities, a more truss-like structure with inferior performance is developed as the density penalization is increased.

Fig. 18 contains the objective evolution using both approaches together with snapshots of the quarter design domain colored by the Young’s modulus to visualize the effect of the penalization. Designs before and after updating the shifted density field with a continuation scheme are presented for all . A considerable initial increase in the strain energy, , is observed in all cases at early stages of the optimization process. The magnitude of the peaks depends on the strength of the density penalization since more material with intermediate densities is removed as increases. A small penalization slows the hole nucleation process. In the TFC approach, a low initial density threshold further delays hole nucleation; see also Section 5.1.3.

We found that increasing the density shift parameter, , in the continuation strategy can compromise robustness if a small is employed, as seen in Fig. 18(a) for the SFC approach. The evolution of the strain energy is mostly smooth with small influences of the parameter for = [2.0, 3.0]. However, due to a more diffuse density field with intermediate densities, applying the density shift can result in abrupt design changes, as seen for . Since the shift increments become larger with every continuation step if (see Fig. 8), this effect is more pronounced later in the optimization process. This can be especially detrimental in designs that favor thin features and for the SFC approach where intermediate densities are found in the vicinity of the interface by construction. The snapshots of the designs at and for in Fig. 18(a) illustrate this scenario. In contrast, a less noticeable impact of updating the density shift on the design evolution is observed in the TFC approach, as seen in Fig.18(b). In this case the densities transition to 1.0 faster, especially near the interface.

Overall, the stability and robustness of both coupling methods rely to a certain extent on penalizing the intermediate densities. For the problems studied in this paper, values for SIMP penalization in the range of =[2.0-3.0] allow for the formation of fine geometric features and a smoother optimization process while avoiding ill-conditioning. For other problems where a SIMP penalization is insufficient, alternative techniques for removing intermediate densities might need to be considered.

5.2.2 Comparing TO approaches

The mesh dependency and performance of the SFC and TFC in comparison to LS TO with initial hole seeding, LSO (LS only), and standard SIMP are investigated. Final designs obtained by all four TO approaches are shown Fig. 19 for three levels of refinement, i.e., element lengths of . In the LSO setups, 12 holes were initially seeded and a total of 300 design iterations without continuation were employed. The final SIMP designs are post-processed by applying an isovolume filter on the density field such that a total volume of 1.2 is preserved in all designs. Filtered densities are penalized by an exponent of 2.0 in the standard SIMP, SFC and TFC results. The same filter radius is used in all cases; see Table 1.

Figure 19: Final designs using meshes with element sizes =[4.0, 2.0,1.0] of classical SIMP, LS-XFEM only (LSO) with initial hole seeding, and using the SFC and TFC approaches.
Figure 20: Snapshots of void domain for LSO with initial hole seeding, and using the SFC and TFC approaches for .

The SIMP results transition from a truss-like design to a shear-web structure as the mesh is refined (see first row of Fig. 19). Suboptimal performances are attained in these designs due to intermediate densities not completely removed with a constant SIMP exponent of 2.0. The coarser the mesh, the more noticeable is this effect. A SIMP approach that projects the density field (e.g., see Wang et al. (2011)) could have been used instead to mitigate the intermediate densities issue. However, to keep consistency in the comparison and assess the resemblance of the density problem and the coupling approaches in terms of their topologies, the same density penalization scheme is used in all cases.

The LSO final designs are continuous walls connecting the two flanges at the top and bottom of the beam with a thickness of about the element edge length. Wrinkles are observed in all three levels of refinement. However, they are reduced as the mesh is refined. These spatial oscillations in the topology are typically mitigated by using either a larger filter radius or higher order interpolation functions. Nevertheless, in order to keep consistency in the comparison with the proposed approaches, the same filter radius and linearly interpolated fields are used in all runs.

As reported in the last two rows of Fig. 19, the SFC and TFC approaches experience smaller variations in final topologies and performance with mesh refinement. In both cases, the density field promotes final designs with shear-webs that are typically seen for LS-XFEM results with initial hole seeding. Despite using a density field, they retain the desired LS-XFEM feature of having an increased geometrical resolution in design with coarser meshes. Furthermore, the TFC approach seems to provide more consistent results with mesh refinement. Note also that the SFC and TFC strategies have an inherent tendency to form bulkier features due to the SIMP penalization. In settings where the initial hole pattern avoids triggering suboptimal final designs, like the one shown in this example, pure LS results are marginally better than using the SFC and TFC strategies. However, this is might not the case in general.

Fig. 20 shows snapshots of the void domain to highlight the hole nucleation process of the LSO, SFC, and TFC beam designs in Fig. 19 using a mesh of size . The evolution of the material removed by each approach is shown at =[5, 25, 50, 100]. The initial hole pattern in the LSO approach favors the removal of material along the external boundaries of the domain, which eventually results in a structure with internal shear-webs. In contrast, the holes nucleated in both the SFC and TFC are concentrated inside the design domain. Thus, in this particular case, the hole seeding promotes the formation of the shear-webs along the external boundaries.

5.3 Example 3: bracket

Finally, the proposed hole nucleation approaches are applied to a complex engineering design problem. The complexity of this design problem stems from the 3D geometry of the design domain and the multiple objective components and constraints.

The goal of this design problem is to find a structure that supports a payload given a set of supports and bolts for attaching the structure to the payload. Fig. 21 shows the design domain () colored in light grey, while part of the non-design domains for the payload box (), the supports (), and bolts () colored in dark grey. Bolts and supports are modeled via hollow cylinders. Side and bottom views are provided in addition to a 3D view to better visualize the constrained subdomains. A uniform pressure load acts on the top surface of the payload box. The entire structure is subject to a body force in direction, representing an equivalent shock loading. The material properties and load conditions for this static problem can be found in Table 5.

Property Value
Young’s Modulus (, solid) x [N/]
Young’s Modulus (, void) x [N/]
Material Density (, solid) x [kg/]
Material Density (, void) [kg/]
Poisson Ratio ( and ) = 0.342 [-]
Pressure load () 1.2x [N/]
Maximum stress () 398.7 [N/]
Table 5: Material properties and load conditions employed in bracket problem.
Parameter Value
Continuation step size
Number of design iterations in continuation
Maximum number of design iterations
Continuation density threshold exponent
Continuation density shift exponent
Initial density threshold
Initial density shift
Table 6: Continuation parameters for the bracket problem.
Figure 21: Bracket problem design and non-design domains in multiple views.
Figure 22: 3D and plane views of the locally refined mesh used in the bracket problem.

The objective of this optimization problem is to minimize the mass, , of the structure that supports the payload box. The optimization problem formulation is the following:


The first term in the objective is the mass of the design domain to be minimized. The second term is a strain energy component added to prevent an overly aggressive removal of mass early in the optimization process and promote a sufficiently stiff structure. The penalty term, , is formulated as follows:


with the field defined as:


based on the previously defined subdomains highlighted in Fig. 21. Through the field defined in Eq. 31, the LS penalization in Eq. 30 (i) forces the support structure to stay within the design domain, (ii) prevents contact between the payload box and the support structure except for the bolts, and (iii) avoids altering the LSF in the payload box, attachment bolts and the supports. The remaining penalty terms in the objective (, , and ) are defined in Section 4.2. Similar to the beam example, the perimeter weight, , is updated at every continuation step from 0.0001 to 0.01, using Eq. 13 with an exponent of 3.0.

In addition to considering the mass in the objective, a mass constraint (, with ) is imposed to represent an upper limit on the mass. The stress constraint, , enforces the volume of solid in which the stress exceeds the maximum von Mises stress, , to zero. The stress penalty,


is an integral over the volume of the scalar stress field, , defined as: