level set formulation for topology optimization of functionally graded lattice structures

Computers and Structures 231 (2020) 106205

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

A hybrid density/level set formulation for topology optimization of functionally graded lattice structures M. Jansen ⇑, O. Pierard Cenaero, Rue des Frères Wright 29, 6041 Gosselies, Belgium

a r t i c l e

i n f o

Article history: Received 26 September 2019 Accepted 13 January 2020

Keywords: Topology optimization Additive manufacturing Lattice structures Homogenization Level set XFEM

a b s t r a c t Additive manufacturing is capable of producing lightweight structures by introducing mesoscale lattice structures in the design without significant additional manufacturing costs. Nevertheless, designing efficient structures including complex lattices remains a challenging task. This paper presents a strategy for the design of additively manufactured structures that is able to simultaneously optimize the shape of the macroscale structure and the functionally graded lattices inside this design. Assuming a separation of scales between the lattice cell size and the macroscale structure, an effective material model based on numerical homogenization is adopted to model the lattice behavior in the analysis of the macroscale structure. While this material model is parametrized by volume densities to mimic functionally graded lattices, the shape of the macroscale structure is represented by a level set function and modeled by the extended finite element method inside the design domain. By adopting an explicit approach to the level set optimization the design problem can be formulated as a single nonlinear programming problem which can be solved with standard optimization algorithms for topology optimization. The proposed formulation avoids degenerate lattice members by excluding small densities in the optimization. Furthermore, the distinct representation of the macroscale structure enables to easily include geometric constraints on the macroscale level such as a minimum feature size or other constraints specific to the additive manufacturing process. The effectiveness of the method is demonstrated in a number of examples. Ó 2020 Elsevier Ltd. All rights reserved.

1. Introduction The technological progress in recent years has made additive manufacturing (AM) a competitive alternative to traditional subtractive manufacturing techniques. This has led to a widespread adaptation of AM in manufacturing industries such as the automotive and aerospace sectors. AM is capable of producing highlycomplex parts without significant additional costs as long as the properties of the manufacturing process are taken into account in the design process. Lattice structures are a good example of complex structures that can be easily produced by metal AM processes in order to generate lightweight and efficients parts. These microor mesoscale structures have a number of interesting features such as a high strength/mass ratio and decent energy absorption and thermal insulation properties, which make them useful for different physical applications [1]. Furthermore, careful design of functionally graded lattice structures enables to generate efficient designs while also improving the manufacturability of the part. ⇑ Corresponding author. E-mail address: [email protected] (M. Jansen). https://doi.org/10.1016/j.compstruc.2020.106205 0045-7949/Ó 2020 Elsevier Ltd. All rights reserved.

Traditionally, software packages for the design of AM allow the user to manually replace bulk parts of the design by lattice structures in order to reduce the weight of the part. This process is typically accompanied by a numerical simulation, e.g. using the finite element method (FEM), to verify the structural integrity of the altered design. Topology optimization has emerged as a powerful tool to automate and improve the design of structures produced by additive manufacturing. This synergy between AM and topology optimization has therefore received considerable attention in the scientific community; several review papers [2,3] provide overviews of these recent developments in the field. Regarding lattice structures, topology optimization methods can automate (parts of) the manual trial-and-error process described above and help the designer in finding efficient solutions. This work also focuses on topology optimization for the automatic design of efficient parts containing functionally graded lattice structures. Only a selective overview of the existing methods is provided here due to the fast succession of new developments being published on this topic in recent years. Incorporating lattice structures in a model of the macroscale structure is considered as a multiscale problem since the

2

M. Jansen, O. Pierard / Computers and Structures 231 (2020) 106205

dimensions of the lattices are typically relatively small compared to the size of the macroscale part. Based on the assumption of a separation of scales, numerical homogenization is often used to model the lattice structures by means of equivalent material properties. There has always been a strong link between topology optimization and homogenization with many works as proof (e.g. [4,5]). The hierarchical method developed by Rodrigues et al. [6], which relies on homogenization to simultaneously optimize the material distribution and the material properties, is particularly relevant in the context of the current work. Several researchers [7–9] have used numerical homogenization to optimize the functional grading of lattice structures. These works use a single density variable to parametrize the functional grading of isotropic lattice structures. Several methods [10,11] have been developed that allow for the dimensions of the members in the different principal directions of the cubic cell to vary independently. By also permitting the unit cells to rotate, these methods are able to generate lattices with orthotropic properties that are truly adapted to the macroscale loads. White et al. [12] used neural networks in order to model more complex orthotropic unit cells with a larger number of design variables. Note that there are also alternative strategies to deal with the multiscale problem that do not use numerical homogenization. For example, Alexandersen and Lazarov [13] proposed to use specialized preconditioners in order to perform multiscale topology optimization without any separation of scale. Wu et al. [14] developed a method based on finite element sub-structuring in order to incorporate lattice structures in topology optimization. Infill material is another technique that is often used in AM to reduce weight. Although these methods are typically used for Fused Deposition Modeling (FDM) with plastic materials, the infill strategies in topology optimization are strongly related to the problem of functionally graded lattice structures. The methods of Wu et al. [15,16] directly optimize the infill at the macroscale by incorporating geometric constraints similar to maximum length constraints in density-based topology optimization. Dapogny et al. [17] proposed a method to optimize the anisotropic infill in FDM parts. Their level set-based formulation optimizes the macroscale structure and its infill in two sequential steps. Wang et al. [18] also investigated a two step sequential approach to optimize both the macroscale shape and the microscale structures in a level set (LS) formulation. This work presents a method to optimize both the shape of the macroscale structure and the functionally graded lattices inside the part. The mechanical behavior of the lattice structures inside the part are approximated by an equivalent material model based on numerical homogenization. The current work limits itself to a single type of unit cell with fixed orientation where the geometric properties of the cells are parametrized by a single volume density. On the other hand, the macroscale shape is described by a level set approach which is augmented with the extended finite element method (XFEM or LS XFEM) [19] to track the shape on a nonconformal mesh of the design domain. In topology optimization, the XFEM is able to maintain a sharp material/void interface represented by a level set on a non-conformal mesh of the design domain [20–23]. The level set is parametrized explicitly by a set of independent optimization variables. This explicit approach enables to formulate a single nonlinear optimization problem in both the density and shape variables. Solving this problem simultaneously optimizes the macroscale shape and its functional grading. Furthermore, this hybrid formulation combining both level set and density-based variables is able to benefit from the advantages of both methods. In particular, possessing a clearly defined macroscale shape simplifies the inclusion of different manufacturing constraints.

The paper is organized as follows. The following section briefly reviews the parametric level set method that is used to describe and optimize a macroscopic structure consisting of a solid linear elastic material. Next, the optimization of functionally graded lattice structures is considered. This section consists of two parts. First, a computational model for lattice structures on the macroscale is constructed using numerical homogenization. Second, the functional grading of the lattice structures is parametrized using volume densities and a density-based optimization problem is formulated to optimize the spatially varying lattice structures. Section 4 presents the main contribution of this work where the level set parametrization of the macroscopic structure is combined with the density-based lattice model. This hybrid formulation is used to simultaneously optimize the macroscopic structure and lattice distribution inside this structure. Finally, the proposed hybrid strategy is illustrated in a number of examples. 2. Level set XFEM topology optimization of the macroscale structure This section summarizes the parametric level set method that is used in this work to optimize macroscale structures consisting of solid linear elastic material. 2.1. Level set shape parametrization Topology optimization typically confines the structure that is optimized to a certain design domain D
Rd (d 2 f2; 3g). The body of the structure is represented by the material domain X
D,
is considered as while the remainder of the design domain D n X the void phase. Level set methods use a scalar function / : D ! R to distinguish the material domain X inside the design domain D as follows:

8 ðmaterialÞ > < /ðx; pÞ < 0 8 x 2 X ðinterfaceÞ /ðx; pÞ ¼ 0 8 x 2 C > :
ðvoidÞ /ðx; pÞ > 0 8 x 2 D n X

ð1Þ

The definition (1) of the level set function differentiates the material and void phase inside the design domain while the interface C between both domains is described by the implicit iso-zero of the level set function. Fig. 1 illustrates these different domains and level set regions in the well-known example of the MBB beam. Here a socalled parametric or explicit approach to the level set method is adopted [24]. This means the level set function is parametrized by

Fig. 1. Representation of the different domains in the design domain D by means of the level set function /: (a) the material domain X and the interface C inside the design domain D; (b) the level set function /.

M. Jansen, O. Pierard / Computers and Structures 231 (2020) 106205

a vector of independent optimization variables p. The relation between the level set function and the variables p is discussed further on in this section. 2.2. Extended finite element method The applications in this work involve the optimization of structures in an elastostatic setting. The static equilibrium of a continuum structure is described in terms of the displacement field u : X ! Rd by the Navier-Cauchy equations:

div ðC : ðuÞÞ þ b ¼ 0; in X ¼
t; on Cn

ðC : ðuÞÞ  n

ð2Þ

¼ 0; on Cd

u

where  is the linear strain operator, b are the body forces per unit volume and
t are the externally applied tractions on the Neumann boundary Cn . The structure is fixed on the Dirichlet boundary Cd . Finally, C is the symmetric fourth-order constitutive tensor for the linear elastic material of the structure. Since the material domain X is described implicitly by the level set function (1), the extended finite element method for implicit material/void interfaces [25,26] is used to solve the elasticity problem (2). The domain of the displacement field is extended to the fixed design domain: u : D ! Rd . The weak form of the problem (2) is likewise expressed on the complete design domain D while limiting the integration over the material domain by means of the Heaviside function H:

Z D

Hð/Þððv Þ : C : ðuÞ  v  bÞ dx 

Z Cn

v 
t ds ¼ 0; 8v 2 U

ð3Þ

n X Ni ðxÞui ;

8x 2 D

ð4Þ

i

Kð/Þud ¼ f

ð5Þ

where f is the external load vector and ud is the global displacement vector containing all the nodal displacements ui . Given a mesh of the design domain D with ne elements De , the assembly of the stiffness matrix K is expressed as a sum over the element contributions: ne Z X e

In the framework of the XFEM, the level set function / is also discretized on the mesh of the design domain D using first-order shape functions:

/ðx; pÞ ¼

Hð/ÞBT CB dx

n X

Ni ðxÞ/i ðpÞ

ð7Þ

i¼1

where the nodal values /i are a function of the independent optimization variables p. Different strategies are used in the literature to parametrize the level set in terms of the optimization variables [24]. In this work, a design variable pi is assigned to each node of the mesh and collected in the vector p 2 Rn of optimization variables. The level set values are computed by a linear filtering operation which performs a local weighted averaging of the optimization variables p [23]:

Pn j¼1

/i ¼ Pn

jij v j pj jij v j

ð8Þ

j¼1

where the discretized nodal weights v j and the filter weights jij are based on a conic kernel function with the user-defined filter radius R/ :

vj ¼ jij

Z

Nj ðxÞ dx     ¼ max R/  xi  xj 2 ; 0

ð9Þ

D

ð10Þ

where xi are the coordinates of node i. In the following, we will often refer to the linear filter (8)–(10) using the following shorthand notation:

where / represents the vector of all nodal level set values. Besides smoothing the level set function /, the filter operation also diffuses inaccuracies in the numerical shape sensitivities in order to improve the convergence of the optimization algorithm [29]. As usual in topology optimization, the structure is optimized for minimum compliance with respect to external loads while limiting the total volume of the structure. Assuming there are no body loads, the compliance of the structure for the external loads is defined as:

Z

F¼

where ui are the nodal displacements and n is the total number of nodes in the mesh. The linear finite element system is obtained by introducing the interpolation (4) for the (virtual) displacements u and v in the weak form (3):

K¼

2.3. Optimization problem

ð11Þ

where v are virtual displacements in the set of admissible displacements U satisfying the homogeneous Dirichlet conditions. The extension of the domain in Eq. (3) enables to work with a fixed mesh of the design domain in the finite element discretization. The displacement field is discretized on this mesh using standard first-order shape functions N i [27]:

uðxÞ ¼

3

ð6Þ

De

where B is the discretized linear strain matrix. It can be seen that the element stiffness (6) reverts to the standard finite element expression for elements that are entirely in the material domain and regular Gauss integration can be used to evaluate these contributions. In the elements cut by the iso-zero level set / ¼ 0, the Heaviside function is taken into account by a specialized integration scheme based on a sub-triangulation strategy. Fig. 2 illustrates this modified integration rule for a quadrilateral element.

Cn

u 
t ds

ð12Þ

where the displacement u are the solution of the weak form (3). The total volume of the structure is expressed as follows:

Z

V¼

X

Z

dx ¼

Hð/Þ dx

ð13Þ

D

Note that the numerical evaluation of this integral again relies on the XFEM integration strategy described above. Using the definitions (12) and (13) and the parametrization (11) of the level set function, the minimum compliance problem is formulated as the following nonlinear programming problem with optimization variables p:

minn F ðuðpÞ; pÞ p2R

s:t: Gvf ðpÞ ¼ V ðpÞ  vf V D 6 0  1 6 pi 6 1 i ¼ 1; . . . ; n

ð14Þ

where the volume V of the structure is limited to the fraction vf of the volume V D of the design domain D. The lower and upper bounds on the design variables p are chosen arbitrarily. The numerical solution of the optimization problem is discussed in Section 4, but first the modeling and optimization of functionally graded lattice structures is considered.

4

M. Jansen, O. Pierard / Computers and Structures 231 (2020) 106205

Fig. 2. XFEM integration scheme for a quadrilateral element cut by the iso-zero level set based on sub-triangulation. Integration points for a second order Gaussian quadrature rule are indicated by a cross symbol  (source: [28]).

3. Optimization of functionally graded lattice structures

q ¼ 1  ð1  2rl Þ2 ¼ 4rl  4r2l

This section discusses the optimization of the functional grading of lattices in a macroscale structure. Assuming a separation of scales between the meso- or microscale lattices and the macro structure, the lattice behavior in the macroscale analysis is modeled by an equivalent material law based on numerical homogenization. Note that apart from the design regularization using nonlinear projection filters, the presented strategy is in fact quite similar to existing methods in the literature [7–9].

Note however that in the present implementation the r l ðqÞ relation is determined numerically in both 2D and 3D by computing jXl j from the FE mesh used during the numerical homogenization routine explained in the following. Numerical homogenization methods for linear elastic media are well described in the literature; see, for example, References [30,31]. The effective properties are determined by analyzing the behavior of the RVE when a set of constant macroscopic strain
ð0ijÞ are applied. Here the components of the constant strain fields 

3.1. Numerical homogenization of lattice structures


ð0ijÞ are prescribed as: tensors 

The goal of numerical homogenization is to determine an equivH

alent constitutive model C that can be used to model the behavior of the lattice structures in an elastostatic analysis of the macroscopic structure. The homogenization analysis is performed by considering a representative volume element (RVE) which corresponds to the cube (or square in 2D) enclosing a unit cell. This work focuses on the distribution of lattices limited to a single base cell with a fixed topology in order to avoid the complications that arise when neighboring lattice cells are geometrically incompatible. Fig. 3 shows the RVEs for a 2D and 3D cubic unit cell, where the domain Xl indicates the material domain of the unit cell inside the domain XRVE of the RVE. The geometry of the unit cell is described by its length Ll and the radius Rl of the bars (see Fig. 3). Assuming the length Ll of the unit cell is fixed, the functional grading is varied by changing the relative radii rl ¼ Rl =Ll of the members of the cell. Moreover, the radius rl can be related to the relative volume density q ¼ jXl j=jXRVE j of the unit cell in order to express the functional grading problem in a form similar to a standard topology optimization problem. For the 2D cell this relation is found in analytical form:

ijÞ 
ð0;kl ¼



ð15Þ

1

if kl ¼ ij or lk ¼ ij

0

else

ð16Þ

A set of elastostatic problems are solved on the RVE where Dirichlet boundary conditions are imposed on the exterior boundary corresponding to the displacements induced by the constant strain
ð0ijÞ : fields 

   div C :  uðijÞ ¼ 0; u

¼


ð0ijÞ

ðijÞ

in Xl  x; on CRVE \ Cl

ð17Þ

where CRVE and Cl are the external boundaries of respectively the RVE domain XRVE and the lattice material domain Xl . Given the   microscopic strain fields ðklÞ ¼  uðklÞ corresponding to the solutions uðklÞ of the different problems (17), the components of the homogenized material tensor are determined as follows:

R ðklÞ 1 C Hijkl ¼ jXRVE E  dx j Xl ijpq pq R 1 ¼ jXRVE rðklÞ dx j Xl ij

ð18Þ

It is clear from its geometry that the effective properties of the cubic unit cell corresponds to a cubic material law, i.e. an orthotropic law with additional 90 rotational symmetry around its three principal axes. A cubic material law has only three independent components, i.e. C H1111 ; C H1122 , and C H1212 . Furthermore, it suffices to consider two
ð011Þ load cases in problem (17), corresponding to the unit normal 
ð012Þ strain, to completely determine the constitutive and shear  tensor:

Z C H1111 ¼

Xl

Z C H1122 ¼

Xl

Z C H1212 ¼

Xl

rð1111Þ dx

ð19Þ

rð2211Þ dx

ð20Þ

rð1212Þ dx

ð21Þ

A routine has been implemented in order to determine the relation Fig. 3. RVEs for periodic cubic lattice structures: (a) periodic lattice structure in 2D, (b) RVE for a 2D cubic cell and (c) RVE for a 3D cubic cell.

C Hijkl ðqÞ numerically. The following steps are performed for different values of the relative bar radius r l :

5

M. Jansen, O. Pierard / Computers and Structures 231 (2020) 106205

Fig. 4. Finite element meshes of the cubic unit cell for different relative volume fractions.

1. Perform the finite element analyses of the problems (17) for
ð012Þ ;
ð011Þ and  unit strain fields  H

2. Determine C by evaluating the integrals (19)–(21); 3. Determine r l ðqÞ by computing jXl j numerically; The meshes for the different unit cell configurations are generated by the Gmsh library [32]. Fig. 4 shows some discretizations of the 3D cell for different density values. The finite element analyses are done using our in-house finite element software. Fig. 5 illustrates some post-processing results of these analyses. The routine described above generates for each material prop  ðkÞ  erty a sequence of data points such as qðkÞ ; C H . In the fol1111 q lowing, a piecewise linear interpolation based on these data sets is used to evaluate the functions C H ijkl ðqÞ. Finally, Fig. 6 illustrates the

density dependence of the constitutive stiffness components by H showing the normalized values C H ijkl ¼ C ijkl =C ijkl as a function of

the volume density q. The dots in these graphs correspond to the data points calculated by means of homogenization. Only the first and last point, corresponding to pure void (q ¼ 0) and bulk material (q ¼ 1) have added afterwards. A small ersatz stiffness

(¼ 109 C ijkl ) is attributed to the void phase (q ¼ 0) in order to ensure a positive-definite stiffness matrix. 3.2. Optimization of functional grading

The parametrized material model of the previous section is used for the optimization of the functional grading of the lattices inside a macroscale structure. A field of volume densities q : D ! ½0; 1 is defined in order to parametrize the functional grading of the lattice structures inside the design domain D. Assuming the design domain D is entirely filled with lattices, the elastic behavior of the structure is found by solving the following weak form of static equilibrium:

Z D

Fig. 5. Illustration of the displacements uðijÞ of the 3D cubic unit cell for boundary
ð011Þ and (b) shear strain conditions corresponding to a uniform (a) normal strain  
ð012Þ .

ðvÞ : CH ðqÞ : ðuÞ dx 

Z Cn

v 
t ds ¼ 0; 8v 2 U

ð22Þ

where the homogenized constitutive tensor CH ðqÞ as determined in Section 3.1 is used. The elasticity problem (22) is solved by means of the standard finite element method. A finite element mesh of the domain D is generated and the displacements are discretized as in Eq. (4). The density field q is also interpolated using the same piecewise linear shape functions in order to ensure a continuous variation of the functional grading:

qðxÞ ¼

n X Ni ðxÞqi ;

8x 2 D

ð23Þ

i¼1

The displacement (4) and density (23) interpolations are used together with standard Gaussian quadrature to assemble the following global linear finite element system:

KðqÞud ¼ f

ð24Þ

where the vector q 2 Rn collects all the nodal density variables qi . In the field of topology optimization it is well known that some design regularization on the physical densities q has to be incorporated in order to obtain a well-posed optimization problem with mesh independent solutions [33]. Moreover, it is preferable to limit the spatial variation of the functional grading (i.e. densities) due to manufacturing limitations. In this work, this is achieved by including a density filter in the design parametrization [34,35]. A new set
2 Rn is introduced and the volume denof nodal design variables q sities q are determined by applying the linear filter function (11):

ð25Þ Fig. 6. Normalized homogenized material properties for the cubic unit cell as a function of the relative volume density q in (a) 2D and (b) 3D.

The spatial variation of the density field is then limited by choosing an appropriate filter radius Rq .

6

M. Jansen, O. Pierard / Computers and Structures 231 (2020) 106205


Þ; q
Þ minn F ðuðq

The functional grading is optimized for minimum compliance similar to the macroscale design problem (14). Given the density field q the total volume of material is defined as:

Z

V¼ D

q dx

ð26Þ

Using the parametrization of the density field (25), the solution of Eq. (22) in the compliance (12) and the volume functional (26), the minimum compliance problem for the functional grading is expressed as follows:


Þ
Þ; q minn F ðuðq q
2R

ð27Þ


Þ  vf V D 6 0
Þ ¼ V ðq s:t: Gvf ðq

ql 6 q
i 6 qu i ¼ 1; . . . ; n where the lower and upper bounds ql and qu on the volume densities are chosen such that 0 < ql < qu < 1. The choice of the bounds ql and qu should be based on the geometrical limitations of the lattice unit cells, i.e., in order to avoid degenerate lattice members and very small holes in the optimized design. For example, for the 2D cubic cell, the lower bound ql corresponding to a given minimum member radius Rl can be determined using Eq. (15). Since the density field is limited to the interval ½ql ; qu , problem (27) can only be used to optimize the functional grading without changing the topology of the macroscale structure; i.e. introduce holes (q ¼ 0) or solid material (q ¼ 1). In order to perform topology optimization including lattice structures, Groen and Sigmund [10] proposed a nonlinear projection filter that is able to generate material/void density values (q 2 f0; 1g) in combination with intermediate densities q 2 ½ql ; qu  corresponding to lattice struc
is introduced for this tures. A third field of intermediate densities q purpose, and the volume densities q are determined in two steps.
are filtered using the linear filter (25) First, the design variables q to generate the intermediate densities . Second, the filter projects densities below the threshold ql to 0 and densities above the threshold qu to 1 by means of the following nonlinear function:

q ¼ q~ ð1  Hc ðq~ ; qu ; bÞÞHc ðq~ ; ql ; bÞ þ

~ b1þq ~ ; qu ; bÞ Hc ðq b

ð28Þ

where b is the steepness parameter of the following differentiable approximation Hc of the Heaviside function [36]:

~ ; qu ; bÞ ¼ H c ðq

~  q u ÞÞ tanhðbqu Þ þ tanhðbðq tanhðbqu Þ þ tanhðbð1  qu ÞÞ

ð29Þ

The projection function (28) is illustrated in Fig. 7 for different values of the parameter b. As explained in Ref. [10], a continuation on the parameter b is typically used during the optimization due to the highly nonlinear behavior of the projection (28). With the inclusion of the projection filter (28), it is possible to extend the bounds in the optimization problem (27) and allow
2 ½0; 1: the design variables to vary over the complete interval q

Fig. 7. Nonlinear function (28) that projects densities close to the bounds.

q
2R

ð30Þ


Þ  vf V D 6 0
Þ ¼ V ðq s:t: Gvf ðq
i 6 1 i ¼ 1; . . . ; n 06q

Note however that a limitation of this single field solution is that the filter (28) is not able to represent a sharp transition between pure solid phase to void phase: in combination with a finite filter radius Rq in the density filter (25) the transition from solid to void will always contain a region of intermediate lattice densities. This phenomenon will be illustrated in the example of Section 5.1 where the present density formulation including the projection filter is compared to the hybrid formulation presented in the following section. As a result of this limitation, the exterior surface of the macroscale structure is not well defined, which complicates the definition of geometric or manufacturing constraints for the design. 4. Simultaneous optimization of the macroscale structure and its functional grading 4.1. Formulation of the optimization problem In order to alleviate some of the disadvantages of the single density field formulation described in the previous section, an alternative formulation is presented in this section. In this method, the description of the shape of the macroscale structure and its functional grading are separated. The density formulation of Section 3 is combined with the parametric LS XFEM formulation of Section 2 in order to leverage the advantages of both methods. More specifically, the level set parametrization (1) is used to describe the shape of the macroscale structure represented by the material domain X, while the distribution of lattice structures inside this domain X is described by a density field q. It is important to note however that in the same spirit of the XFEM (cf. Section 2) the density field q is still defined on the complete design domain D. In this way the simple linear interpolation (23) can be retained. Furthermore, this also greatly simplifies the design sensitivity analysis as will be discussed in Section 4.2. In this hybrid formulation, the total volume of the design can be expressed by combining Eqs. (13) and (26) as follows:

Z

V¼ D

qHð/Þ dx

ð31Þ

Similarly, the static equilibrium of a structure described by the hybrid formulation is found by introducing the homogenized material model of Section 3.1 in the XFEM weak form (3):

Z

D

Z   Hð/Þ ðv Þ : CH ðqÞ : ðuÞ dx 

Cn

v 
t ds ¼ 0; 8v 2 U

ð32Þ

Integral expressions in the hybrid formulation, such as Eqs. (31) and (32), are evaluated numerically using the same XFEM strategy as outlined in Section 2 in order to deal with the presence of the Heaviside function H. In accordance to Sections 2 and 3, the physical design fields in the hybrid formulation, i.e. the level set function / and the physical volume densities q, are parametrized by a set of independent design variables. In the hybrid formulation, the domain X of the macroscale structure is unambiguously distinguished from the void domain by the level set function /. As a result, degenerate lattice members can be eliminated simply by setting an appropriate lower bound ql on the volume densities q. However, this does not solve the problem of very small holes related to densities that are close the upper bound q ¼ 1. In order to resolve this issue, the physical densities q are still parametrized by means of a two-step ~ ðq
ÞÞ, as explained at the end of Section 3 but with a filter, i.e. qðq modified projection function. First, the intermediate densities

M. Jansen, O. Pierard / Computers and Structures 231 (2020) 106205


. are obtained by filtering the design variables q Second, a modified version of the projection filter (28) is used which only projects densities above a threshold value qu to solid material (q ¼ 1):

q ¼ q~ ð1  Hc ðq~ ; qu ; bÞÞ þ

~ b1þq ~ ; qu ; bÞ H c ðq b

ð33Þ

This modified projection function is illustrated in Fig. 8 for different values of the parameter b. An alternative to this approach would be to consider the functional grading (ql 6 q 6 qu ) and the solid material (q ¼ 1) as distinct phases. The resulting three-phase problem could be modeled for example by a colored level set approach [37]. This alternative formulation is not investigated in this work. The level set parametrization /ðpÞ using the linear filter (8) remains unchanged. This explicit level set approach is the main factor that enables to formulate the hybrid density/level set formulation as a single nonlinear optimization problem. As a result, the macroscale structure and the functional grading inside this structure can be optimized simultaneously by solving this optimization problem. The hybrid density/level set optimization problem is formulated as follows in terms of the independent optimization vari
; pÞ: ables ðq

min n

q
2R ;p2Rn


; pÞ; q
; pÞ F ðuðq


; pÞ  vf V D 6 0
; pÞ ¼ V ðq s:t: Gvf ðq

ð34Þ

ql 6 q
i 6 1 i ¼ 1; . . . ; n 1 6 pi 6 1 i ¼ 1; . . . ; n
; pÞ that are used in the compliance F where the displacements uðq are the solution of the weak form (32). The nonlinear optimization problem (34) can be solved by means of standard nonlinear programming algorithms [38]. In particular, the examples in this work (see Section 5) use the Method of Moving Asymptotes (MMA) [39] which is a gradient-based algorithm that is well adapted for topology optimization. 4.2. Sensitivity analysis Since the optimization problem (34) is solved by means of a gradient-based optimization algorithm, an efficient method is required for computing the design sensitivities of the objective and constraint functions with respect to the optimization variables
; pÞ. ðq First the design sensitivities with respect to the physical design variables are considered, i.e. the volume densities q and the level set function / parametrizing the material domain X. Since the density field q is defined on the complete design domain D, it is independent of the material domain X. As a result, the sensitivity analysis for the density variables q and the shape parameters / can be performed independently. This means that the partial derivatives with respect to the density variables q can be derived

Fig. 8. Modified projection filter (33) that projects densities close to the upper bound 1.

7

following the same steps as for regular density-based topology optimization [33]. Likewise, the partial derivatives with respect to the nodal level set values / are found following standard techniques for parametric shape derivatives, see e.g. [24]. Since these elements are well known in the topology optimization community, this section only presents the final formulas for the sensitivity analysis of the volume (31) and compliance (12) without going into detail. Using the interpolations (23) and (7) for the density field q and level set function /, the partial derivatives of the total volume (31) with respect to density and level set nodal values are expressed as:

Z @V ¼ Ni Hð/Þ dx @ qi D Z @V N ¼ q i ds @/i kr/k C/

ð35Þ ð36Þ

where the parametric shape derivative (36) is only integrated on the movable part C/
C of the interface C. The sensitivity analysis of the compliance (12) is based on the adjoint method, see e.g. [33,40]. Assuming that the displacements u solve the weak form (32), the partial derivatives of the compliance (12) are given as:

! ! Z @F @CH ¼  Hð/Þ ðuÞ : Ni : ðuÞ dx @ qi @q D Z   N @F i dx ¼ ðuÞ : CH ðqÞ : ðuÞ @/i kr/k C/

ð37Þ ð38Þ

Once the partial derivatives with respect to the physical design variables (q; /) have been computed, it is a matter of applying the chain rule of differentiation to obtain the sensitivities with respect to the
; p). For example, the sensitivities of the optimization variables (q volume V with respect to the level set parameters p are determined as: n @V X @V @/k ¼ @pi k¼1 @/k @pi

ð39Þ

where @/k =@pi are found by differentiation of the linear filter in Eq. (8). The design sensitivities with respect to the
are expressed similarly by differentiation optimization variables q and the projection function (33): of the linear filter n X ~k @V @V @ qk @ q ¼
i ~
i @q @ q @ q @ q k k k¼1

ð40Þ

 
; Rq =@ q ~ k =@ q
i ¼ @F k q ~ k are found by analyt
i and @ qk =@ q where @ q ical differentiation of the projection function (33). 4.3. Manufacturability of the macroscale structure An important advantage of the proposed formulation is that the shape of the macroscale structure is unambiguously defined by the level set parametrization of the material domain X. This simplifies controlling the structural complexity of the optimized design and imposing manufacturing constraints in the optimization problem. This section briefly discusses some geometric constraints that can be used to control the structural complexity in the hybrid density/level set formulation. More complex manufacturing constraints such as limitations on overhanging surfaces or internal cavities are beyond the scope of this article and will be the subject of future work. Note that this section is only concerned with the features of the macroscale structure. It is assumed that the two-step filter strategy as outlined in Section 4.1 suffices to ensure the manufacturability of the functional grading.

8

M. Jansen, O. Pierard / Computers and Structures 231 (2020) 106205

With respect to the level set parametrization, it is well known that the filter (8) is not able to limit the complexity of the macroscale structure represented by the material domain X [29]. Moreover, some form of design regularization has to be included in the problem in order to obtain a well-posed optimization problem [40]. Two geometric constraints that can resolve these issues are briefly considered here. The classical solution in level set methods is to include the perimeter, which measures the size of the interface C between the material and void phase, as either a constraint or a penalty term in the optimization. The overall complexity as measured by the perimeter is typically not a limiting factor in additive manufacturing. However, the perimeter is revelant in the post-machining phase which often involves a surface finishing step. For example, in metal additive manufacturing polishing of the exterior surface is often required and the cost of this process is directly related to the size of outer surface. A constraint on the perimeter can be expressed as follows:

Z

GP ¼

C

ds  Pmax 6 0

ð41Þ

where P max is the maximum allowed perimeter. The perimeter provides a simple measure of the global complexity of a design. However, due to this global nature it is often difficult to choose an appropriate allowable value or penalty factor based on the physical limitations of the manufacturing process. The geometrical resolution of a manufacturing process is often expressed in terms of the minimal feature size that can be produced. It is therefore often useful to directly limit the minimum local feature size and hole size in the optimization. Here the minimum length scale constraints [28] that were recently developed for the parametric level set approach are adopted as constraints on the macroscale structure. These constraints are in fact an adaptation of similar constraints that were originally developed for three-field density-based topology optimization [41]. Their main advantage is that they provide a simple means to impose a minimum length scale in both the material and void phase of the design. The application of these constraints is only briefly summarized here, and the interested reader is referred to these references for a detailed discussion. The main idea is that a minimum length scale is achieved if the design can withstand a certain amount of material erosion and dilation. Using the properties of the filter (11), the erosion and dilation operations are understood as variations of the threshold level set value representing the material/ void interface; i.e. the nominal zero threshold value is allowed to h i vary in the interval /e ; /d where /e and /d represent respectively the erosion and dilation threshold. The level set function should therefore attain at least these threshold values in the respective phases in order to preserve the topology of the design during the boundary variations. These requirements are expressed by the following set of geometric constraints:

Z 1 2 Ið/Þ½/  /e þ dx  min 6 0 VD X Z h i2 1 GV ¼ Ið/Þ /d  / dx  min 6 0 þ V D DnX


GS ¼

ð42Þ ð43Þ

where ½þ ¼ max ½; 0 and min is a tolerance value. The function I serves as an approximate indicator for the inflection zones of the level set function:

  Ið/Þ ¼ exp ckr/k2

ð44Þ

where the constraint parameter c is determined based on the filter radius R/ and the mesh size [41]. As shown in References [28,41], the minimum length scale imposed by these constraints for a

h i particular filter radius R/ and interval /e ; /d can be estimated from the equivalent one-dimensional case. Note that the definition of both the perimeter constraint (41) and the minimum length scale constraints (42) and (43) only depend on the level set function / and are entirely independent of the functional grading q. 4.4. Post-processing of the optimized design This section describes a simple strategy to post-process the optimized hybrid design, consisting of the density field q and the level set /, in order to generate a detailed geometric model of the design which includes the individual lattice truss members. The method creates a new level set function /l : D ! R that represents the structure including the individual lattice cells. Since the homogenized material properties and the relative volume densities are independent of the absolute cell size Ll , the designer has to select an appropriate cell size. In theory, the cell size should be significantly smaller than the macroscale feature size in order to ensure the validity of the separation of scales used in the homogenization. However, in practice, homogenized material properties often provide adequate results even for relatively large cell sizes. For example, Hollister and Kikuchi [42] reported acceptable results for a ratio of 0:2 between the cell and macroscale size. The procedure starts by generating a new mesh Dhf of the design domain D that is sufficiently fine to capture the geometry of the lattice cells; i.e. the element size should be a fraction of the lattice cell size Ll . For example, in the 2D problems in the following section an element size Ll =10 is typically used. The optimized fields ðq; /Þ are interpolated on this fine mesh Dhf . A field Rl ðxÞ : D ! R of approximate lattice radii is created using the density field q and the function relating density to lattice geometry (e.g. Eq. (15) in 2D). Since the orientation and the size Ll of the lattice cells is fixed a priori, the lattice truss forms a regular grid inside the material domain X and the location of the axes of the bars is known in closed form. For example, the bars in the global x-direction are represented by the following set:

X l1 ¼ fse1 þ iLl e2 þ jLl e3 js 2 R; i 2 N; j 2 Ng

ð45Þ

where ei are the unit-vectors of the Cartesian coordinate system. For each point x 2 D, the distance to the nearest lattice axis in each la

direction is represented by the vector-valued function d : D ! Rd : la

d ðxÞ ¼ jx  bx=Ll eLl j

ð46Þ

where be is the nearest integer function. The unsigned distance la

function dm ðxÞ to the lattice axes is then defined as: la dm ðxÞ

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2  2 la la ¼ minkx  yk ¼ di1 ðxÞ þ di2 ðxÞ

ð47Þ

y2X l

where X l ¼

Sd i

X li , and the indices i1 and i2 refer to the two smallest la

components of d . Assuming for now that the entire design domain is filled with lattices, a level set function that delineates the domain la

of lattice truss members is found by offsetting dm ðxÞ with the local radius Rl ðxÞ:

~ l ðxÞ ¼ dla ðxÞ  Rl ðxÞ / m

ð48Þ

~ l ðxÞ assumes the entire design domain D is filled with The function / lattices. The solution still has to be limited to the material domain X. Here a simple strategy is followed to clip the lattice structures to the material domain while also providing adequate closure of

M. Jansen, O. Pierard / Computers and Structures 231 (2020) 106205

the lattice structures that are cut on the interface C by projecting them onto this interface. First, a signed-distance function /d : D ! R with respect to the material domain X is computed, for example, by applying the fast marching method [43] to the optimized level set function /. The distance function /d is used to evall/

uate the level set function dm which approximates the projection of the regular lattice grid members in the vicinity of the interface (/d ¼ 0) onto this interface. For this purpose, the level set /d is lin
assuming that the length earized locally at the current point x ¼ x scale of variation of the interface is much larger than the lattice cell size. This linearization is used to transform the global coordinates y 2 Rd to the following local coordinates s 2 Rd :


Þðy  x
Þ
Þ ¼ T / ðx sðy; x

ð49Þ



T where T/ ¼ e/1 ; e/2 ; e/3 is the transformation matrix of the local Cartesian coordinate system whose third unit vector corresponds to the normal on the level set surface e/3 ¼ r/d and the first two axes
Þ. e/1 and e/2 are in the tangent plane on the surface /d ¼ /d ðx Since the first two components of s represent the coordinates in the tangent plane of the level set /d , the minimal in-plane distance sm to the lattice bars projected onto this surface is found by solving the following optimization problem:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sm ðxÞ ¼ min s21 ðy; xÞ þ s22 ðy; xÞ

ð50Þ

y2X l

pffiffiffi s:t:js3 ðy; xÞ  /d ðxÞj 6 3Lc =2 where the constraint on s3 ensures that only the lattice bars that are sufficiently close to the interface /d ¼ 0 are considered in the prol/

jection. The distance function dm is then defined by adding the out-of-plane distance /d ðxÞ as follows: l/

dm ðxÞ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðsm ðxÞÞ2 þ ð/d ðxÞÞ2

ð51Þ

Note that the problem (50) represents a set of linear least squares problems (one for each lattice bar enclosing the point x) which can be solved analytically. Furthermore, in two-dimensional probl/

lems this simplifies to dm ðxÞ ¼ j/d ðxÞj since the projected distance sm ¼ 0. The different functions are combined into the final level set function /l representing the lattice structures limited to the material domain X:



/l ðxÞ ¼

/d ðxÞ

if/d ðxÞ > 0

/ls ðxÞ

if/d ðxÞ 6 0

ð52Þ

where

n h i o la l/ /ls ðxÞ ¼ max /d ðxÞ; min dm ðxÞ; dm ðxÞ  Rl ðxÞ

ð53Þ

Fig. 9 illustrates the effect of this projection step on a functionally graded cube for different distance functions /d . It can be seen that the lattice structures are properly projected onto the iso-zero of this level set function. Note that the signed distance functions in these examples have been chosen solely for illustrative purposes. Some of these functions create macroscale features with a size similar to the lattice cell size, while in practice, the lattice cells should be small compared to the macroscale structure to ensure the validity of the homogenization approach. This post-processing procedure has been scripted using Gmsh [32] and the mesh file library Meshio [44]. The function /l obtained following these steps can be used, for example, in an XFEM analysis on the fine mesh Dhf of the design domain. In the examples of Section 5, such an XFEM analysis of the post-processed design will be performed as a sort of verification of the optimized design. Due to the large number of degrees of freedom associated with the fine

9

mesh Dhf , the computational cost and memory consumption of this analysis is typically significantly larger than for the finite element analysis used in the optimization (based on the equivalent material material). In particular, whereas the linear system of the latter can be solved efficiently by a direct linear solver on a modern laptop, the fine scale analysis is done using an iterative GMRES linear solver in the examples. 5. Examples In this section the proposed formulation is illustrated in a number of examples. Two academic optimization problems are considered first. The properties of the base material in these examples are a Young’s modulus E ¼ 1 and a Poisson’s ratio m ¼ 0:3 (where we conveniently omit units assuming consistent dimensions are used). Furthermore, plane stress conditions are assumed in the twodimensional problems. The final example considers an industrial design problem of an aircraft bracket. Note that all the finite element meshes and post-processing views were created using Gmsh [32]. 5.1. Double-clamped beam The first example considers a two-dimensional design problem of a beam to compare the proposed hybrid formulation (34) of Section 4 with the pure density-based formulation (30) of Section 3. The design domain and boundary conditions for the problem are shown in Fig. 10. The beam is clamped at the left and right edge of the design domain. The external load consists of a distributed load p ¼ 0:1 applied on the top side of the domain. The design domain with length L ¼ 60 is meshed with 240  80 bilinear square finite elements. The volume fraction of the structure is limited to vf ¼ 0:5 in the minimum compliance problem. The design is optimized with functionally graded cubic lattice cells modeled by the density-dependent material model of Section 3.1. The parametrization of the density field q uses a filter radius Rq ¼ 1:6 and the lower and upper bounds ql ¼ 0:2 and qu ¼ 0:9. Note that in the density-based formulation (30) both bounds ql and qu are used in the projection filter (28). On the other hand, in the hybrid formulation (34) the lower bound ql is used in the bound constraints of the optimization problem, while the upper bound qu is used in the modified projection filter (33). In the hybrid formulation, a radius R/ ¼ 1:6 is used in the level set filter (8). In addition, the complexity of the macroscale shape is controlled by including the minimum length scale constraints (42) and (43) with the thresholds /e ¼ 0:2 and /d ¼ 0:2. Fig. 11 shows the results for both formulations. The density formulation is initialized from a uniform density field q ¼ 0:5 (Fig. 11 (b)). On the other hand, the hybrid formulation is initialized by a level set (Fig. 11(a)) that generates a material domain with the circular hole pattern shown in Fig. 11(c). This punctured material domain is filled entirely with solid material (q ¼ 1). Fig. 11(d) and (e) show the optimized designs obtained by both methods. It can be seen that both designs contain lattice structures to redistribute the load to the support structure. Although the topology and functional grading of both designs are similar, the density design has a higher compliance (F ¼ 89:35) than the hybrid design (F ¼ 80:40). As explained in Section 3, the projection filter (28) in the density formulation is not able to generate a sharp transition from solid material to void material. This transition zone of intermediate lattice structures is the main reason for the reduced performance. As a final step, the optimized hybrid design of Fig. 11(e) is postprocessed using the method of Section 4.4 to generate a new numerical model of the design including lattices. The lattice cell

10

M. Jansen, O. Pierard / Computers and Structures 231 (2020) 106205

Fig. 9. Projection of the lattice structures corresponding to a constant volume density q ¼ 0:3 onto the solid domain represented by different signed distance functions /d shown in the top row. The bottom row shows the solid domain corresponding to the final level set /l .

size is set to Ll ¼ 0:5. In this two-dimensional problem, a suffi-

Fig. 10. The double-clamped beam: problem setting.

ciently fine mesh Dhf is therefore generated by using square elements with a uniform element size of 0:05. Fig. 12 shows the different intermediate results to generate the final design on this mesh. Fig. 12(a) is simply the projection of the hybrid density/level set design onto this fine mesh. Fig. 12(b) shows the signed-distance level set /d with respect to the material domain X, which is used to generate the level set function /l including the individual lattice members in Fig. 12(c). The corresponding solid domain modeled by the LS XFEM approach is shown in Fig. 12(d). Finally, it can be noted that the compliance of this complex LS XFEM model 12(d)

Fig. 11. Results for the double-clamped beam. The left column shows the results for the density-based formulation with: (b) the initial density field and (d) the optimized density field. The right column shows the results for the hybrid density/level set formulation with: (a) the initial level set, (c) the initial design and (e) the optimized design. In Figures (c) and (e) the density field is shown in the LS XFEM material domain X.

M. Jansen, O. Pierard / Computers and Structures 231 (2020) 106205

11

Fig. 12. Post-processing and verification of the optimized hybrid design for the double-clamped beam. The following results are shown on the fine mesh Dhf used in the verification calculation: (a) the optimized design of the hybrid density/level set formulation; (b) the signed-distance level set /d w.r.t. the material domain of the optimized design in (a); (c) the level set function /l including the individual lattice members; (d) the LS XFEM material domain of the post-processed design.

is almost equal to that of the hybrid density/level set model (Fig. 12 (a)) with values of respectively F ¼ 80:10 and F ¼ 80:37.

R/ ¼ 1:6 in filter (8). In addition, two cases of geometric constraints are considered: [1] the minimum length scale constraints (42) and (43) with

5.2. MBB beam This example considers the minimum compliance design of the well-known MBB beam in order to illustrate the application of the different manufacturing constraints on the macroscale structure of Section 4.3. Furthermore, the efficiency of cubic lattice structures in such a minimal compliance design with geometric constraints is investigated. For this purpose, the hybrid formulation (34) including isotropic cubic lattices is compared to the optimized design consisting only of solid material, which is obtained by solving the LS XFEM problem (14). Fig. 13 shows the design domain with the boundary conditions for the elastic problem: only half of the beam is considered based on the symmetry of the problem. The displacement boundary conditions consist of a rolling support in the bottom right corner and symmetry conditions on the left side of the domain. A unit point load P ¼ 1 is applied in the top left corner of the domain. A regular mesh with 240  80 bilinear square elements is used to discretize the design domain with length L ¼ 60. The lattice structures that are used in the hybrid formulation are expected to mainly contribute in bearing the shear forces in the mechanical problem. Therefore, the cubic lattice cells used in the functional grading are rotated by 45 with respect to the global coordinate system. This rotation is taken into account in the optimization by applying the same rotation to the homogenized constitutive tensor CH . The corresponding density field is parametrized with a filter radius Rq ¼ 1:6 and the lower and upper bounds ql ¼ 0:2 and qu ¼ 0:9. The maximum volume fraction is limited to vf ¼ 0:5 in the optimization. The level set / is parametrized with a filter radius

Fig. 13. The MBB beam: problem setting.

thresholds /e ¼ 0:4 and /d ¼ 0:4 are included; [2] the constraints of Case 1 and a constraint (41) on the total perimeter with Pmax ¼ 160 are included; Fig. 14 shows the results for the different optimization problems. As mentioned earlier, the performance of the lattice structures is assessed by solving the same optimization with only solid material (Eq. (14)). The optimization is initialized from the level set function shown in Fig. 14(a), which is complemented by a uniform density q ¼ 0:5 in the hybrid formulation (Fig. 14(b)). Fig. 14(c)–(d) show the optimized designs without and with lattices for Case 1 of the constraints. The designs have a very similar macroscopic geometry, and their performance is almost equal: F ¼ 190:16 for the solid design versus F ¼ 189:83 for the functionally graded design. Furthermore, the design with functional grading (Fig. 14(d)) only contains some small zones with lattice structures. This result is not very surprising since it is wellknown that isotropic cubic lattices are not highly efficient from a stiffness point of view [45]. Including functional grading therefore does not significantly improve the performance of the design. However, when considering the more strict geometric constraints of Case 2, the optimized designs in Fig. 14(e)-(f) differ considerably with larger zones of lattices appearing in the functionally graded design (Fig. 14(f)). In addition, the performance of the functionally graded design (F ¼ 196:33) is significantly better than that of the solid design (F ¼ 209:67). This simple example illustrates how the functional grading can help improve the manufacturability of the optimized design. As in the previous example, the functionally graded design of Fig. 14(f) is again post-processed using the method of Section 4.4. Since the lattice cell size is set to Ll ¼ 0:5, a uniform mesh Dhf with 14400  4800 square elements is used for this purpose. Fig. 15 again shows the different intermediate steps of the postprocessing method. The projection of the hybrid density/level set design onto the fine mesh is shown in Fig. 15(a). Fig. 15(b) and (c) show respectively the signed-distance level set /d with respect to the material domain X and the level set function /l including the individual lattice members. The final geometry as modeled by the LS XFEM approach based on the level set /l is shown in Fig. 15(d). It can be seen that the level set cutting procedure using the interme-

12

M. Jansen, O. Pierard / Computers and Structures 231 (2020) 106205

Fig. 14. Results for the MBB beam. The left column shows the results for the level set formulation with solid material. The right column shows the results for the hybrid density/level set formulation including functional grading. Figures (a) and (b) show the initial design for the optimization: (a) the level set / which is augmented by (b) the density field q in the functionally graded case. Figures (c) and (d) show the optimized designs for Case 1. Figures (e) and (f) show the optimized designs for Case 2. In Figures (b), (d) and (f) the density field is shown in the LS XFEM material domain X.

Fig. 15. Post-processing and verification of the optimized hybrid design for the MBB beam. The following results are shown on the fine mesh Dhf used in the verification calculation: (a) the optimized design of the hybrid density/level set formulation; (b) the signed-distance level set /d w.r.t. the material domain of the optimized design in (a); (c) the level set function /l including the individual lattice members; (d) the LS XFEM material domain of the post-processed design.

diate level set /d avoids degenerate lattice members and adequately closes the lattice structures in the vicinity of the interface of the material domain. Furthermore, the performance of the homogenized model (F ¼ 196:28) corresponds well to that of the post-processed LS XFEM design (F ¼ 196:15). 5.3. Industrial bracket In order to illustrate the proposed method in a threedimensional setting, the final example considers the optimization of an engine bracket. The optimization problem is based on a public design challenge promoted by GradCAD, Inc. [46]. The problem setting is depicted in Fig. 16(a). The design domain is based on the STEP file that is publicly available on the GrabCAD website [46] to

which the interested reader is referred for a detailed description of the geometry. The bounding box size of the domain is approximately 178:5  108:3  62:5 mm measured in the global coordinate system in Fig. 16. Load condition 3 of the original design problem is used in the minimum compliance problem. In this load case, the displacements are fixed in the four support points, while a uniform traction p ¼ ½6:626; 0; 7:490 N=mm2 is applied on the inner surfaces of the two load points. The elastic properties of the base material are a Young’s modulus 210000 N=mm2 and a Poisson’s ratio m ¼ 0:3. The finite element model uses the mesh shown in Fig. 16(b) consisting of 272458 tetrahedral elements with a maximum element size of 2 mm. This example immediately considers the hybrid density/level set formulation to simultaneously optimize the macroscale shape and functional grading in the

M. Jansen, O. Pierard / Computers and Structures 231 (2020) 106205

13

Fig. 16. Industrial bracket: (a) problem setting and (b) finite element mesh of the design domain.

bracket. The allowed volume fraction in the minimum compliance optimization (34) is set to vf ¼ 0:35. In addition, the design is fixed to solid material in the colored zones in Fig. 16(b) around the support and loads points. The density field for the functional grading uses a filter radius Rq ¼ 5:0 mm with lower and upper bounds ql ¼ 0:1 and qu ¼ 0:8. The macroscale structure is parametrized by a level set with filter radius R/ ¼ 14:0 mm. In addition, the design complexity is limited globally and locally by including both the perimeter constraint (41) with Pmax ¼ 18000 mm2 and the minimum length scale constraint (42) and (43) with thresholds /e ¼ 0:4 and /d ¼ 0:4. As a side note, to the author’s knowledge, these are also the first published results where the current type of minimum length scale constraints [41,28] are applied successfully on an unstructured 3D mesh. Fig. 17 summarizes the results of the optimization. The initial design for the optimization is shown in Fig. 17(a): the level set is initialized in order to generate a number of (arbitrary) holes in the material domain X while the density field is set uniformly to

q ¼ vf . The optimized design is shown in Fig. 17(b). It can be seen that the tension structure connecting the load points to the supports is mainly composed of solid material. Below the load points there is a part consisting of lattice structures that connects the two sides of these solid struts. A lattice cell size Ll ¼ 3:0 mm is chosen for the post-processing of the optimized design. An adequately fine mesh Dhf has to be generated first for the post-processing phase. Adopting a uniform mesh size based on the lattice cell size for the complete design domain quickly leads to an enormous number of elements in this three-dimensional problem. Since the element size should only be sufficiently small in those places where lattices (q < 1) are located, a mesh with a non-uniform element size based on the optimized density field q is used. Gmsh [32] has a very convenient feature for this purpose that allows to control the mesh size by a user-imported field. The mesh generated in this way (Fig. 18(a)) contains 2446935 tetrahedral elements. The optimized design projected onto this mesh is shown in Fig. 18(b). Fig. 18(c) shows the level set function /l as determined by the procedure described in

Fig. 17. Results for the industrial bracket: (a) the initial design and (b) the optimized design. The density field is shown in the LS XFEM material domain X.

14

M. Jansen, O. Pierard / Computers and Structures 231 (2020) 106205

Fig. 18. Post-processing of the optimized hybrid design for the industrial bracket: (a) the fine mesh Dhf used in the post-processing step; (b) the optimized design of the hybrid density/level set formulation projected on the mesh Dhf ; (c) the level set function /l including the individual lattice members; (d) the LS XFEM material domain of the post-processed design.

Fig. 19. Post-processing of the optimized hybrid design for the industrial bracket: close-ups of the lattice zones in the post-processed design.

Section 4.4. This level set function is used in a subsequent LS XFEM analysis which includes the individual lattice members. The material domain for this XFEM analysis is shown in Fig. 18(d). In addition, Fig. 19 contains some close-ups of the lattice sections. It can be seen that the simple post-processing strategy generates acceptable lattice structures which match the original interface of the optimized design, even on the unstructured tetrahedral mesh. Furthermore, the compliances based on the homogenized material model (F ¼ 84:82 Nmm) and the detailed XFEM model (F ¼ 84:34 Nmm) correspond well.

6. Conclusions This work presents a hybrid density/level set formulation for the topology optimization of mechanical parts that contain functionally graded lattice structures. The shape of the macroscale

structure is modeled by a LS XFEM approach, while the functional grading inside the structure is modeled by a homogenized material law that is parametrized by volume densities. By adopting an explicit LS formulation, a single nonlinear optimization problem can be formulated. Solving this problem simultaneously optimizes the macroscale shape and the functional grading. This hybrid formulation roughly doubles the number of optimization variables in the optimization problem, but the associated increase in computational cost is insignificant compared to the cost of the finite element analysis. Furthermore, the formulation simplifies the inclusion of different manufacturing constraints imposed on the macroscale structure by maintaining a clear distinction between the material and void domain. The proposed method was applied successfully in a number of two-dimensional design problems, and a three-dimensional case of an industrial bracket. A simple post-processing method was applied to the optimized density/level set design to generate a

M. Jansen, O. Pierard / Computers and Structures 231 (2020) 106205

detailed level set representation including all the individual lattice members. This level set model was used to validate the optimized design by means of a high-fidelity XFEM analysis. Our future work will focus on using the current hybrid density/ level set formulation in different applications where functionally graded lattice structures might be useful. We intend to investigate how functionally graded lattices can improve the manufacturability for different types of manufacturing constraints such as overhang limitations, internal cavities and thermal build-up. The physical applications could also be extended to cases such as buckling-limited problems and thermal cooling problems where lattice structures are known to be effective solutions. Finally, the density parametrization of the lattice properties should also be generalized to enable the optimization of orthotropic lattice structures that are adapted to the macroscale loads. Declaration of Competing Interest None. Acknowledgments This research was supported by the SIMATHA project which is part of the IAWATHA convention sponsored by the European Fund for Regional Development and the Walloon Region. References [1] Chu C, Graf G, Rosen D. Design for additive manufacturing of cellular structures. Comput-Aided Des Appl 2008;5(5):686–96. https://doi.org/ 10.3722/cadaps.2008.686-696. arXiv: https://www.tandfonline.com/doi/pdf/ 10.3722/cadaps.2008.686-696. [2] Liu J, Gaynor A, Chen S, Kang Z, Suresh K, Takezawa A, et al. Current and future trends in topology optimization for additive manufacturing. Struct Multidiscip Optim 2018;57(6):2457–83. https://doi.org/10.1007/s00158-018-1994-3. [3] Meng L, Zhang W, Quan D, Shi G, Tang L, Hou Y, et al. From topology optimization design to additive manufacturing: today’s success and tomorrow’s roadmap. Arch Comput Methods Eng 2019;112:841–54. https:// doi.org/10.1007/s11831-019-09331-1. [4] Bendsøe M, Kikuchi N. Generating optimal topologies in structural design using a homogenization method. Comput Methods Appl Mech Eng 1988;71 (2):197–224. https://doi.org/10.1016/0045-7825(88)90086-2. [5] Allaire G. Shape optimization by the homogenization method. New York, USA: Springer; 2002. [6] Rodrigues H, Guedes J, Bendsœ M. Hierarchical optimization of material and structure. Struct Multidiscip Optim 2002;24(1):1–10. https://doi.org/10.1007/ s00158-002-0209-z. [7] Zhang P, Toman J, Yu Y, Biyikli E, Kirca M, Chmielus M, et al. Efficient designoptimization of variable-density hexagonal cellular structure by additive manufacturing: Theory and validation. J Manuf Sci Eng 2015;137(2). https:// doi.org/10.1115/1.4028724. pp. 021004–021004–8. [8] Cheng L, Zhang P, Biyikli E, Bai J, Robbins J, To A. Efficient design optimization of variable-density cellular structures for additive manufacturing: theory and experimental validation. Rapid Prototyp J 2017;23(4):660–77. https://doi.org/ 10.1108/RPJ-04-2016-0069. [9] Cheng L, Bai J, To A. Functionally graded lattice structure topology optimization for the design of additive manufactured components with stress constraints. Comput Methods Appl Mech Eng 2019;344:334–59. https://doi.org/10.1016/j. cma.2018.10.010. [10] Groen J, Sigmund O. Homogenization-based topology optimization for highresolution manufacturable microstructures. Int J Numer Meth Eng 2017;113 (8):1148–63. https://doi.org/10.1002/nme.5575. arXiv: https://onlinelibrary. wiley.com/doi/pdf/10.1002/nme.5575. [11] Allaire G, Geoffroy-Donders P, Pantz O. Topology optimization of modulated and oriented periodic microstructures by the homogenization method. Comput Math Appl https://doi.org/10.1016/j.camwa.2018.08.007. [12] White D, Arrighi W, Kudo J, Watts S. Multiscale topology optimization using neural network surrogate models. Comput Methods Appl Mech Eng 2019;346:1118–35. https://doi.org/10.1016/j.cma.2018.09.007. [13] Alexandersen J, Lazarov B. Topology optimisation of manufacturable microstructural details without length scale separation using a spectral coarse basis preconditioner. Comput Methods Appl Mech Eng 2015;290:156–82. https://doi.org/10.1016/j.cma.2015.02.028. [14] Wu Z, Xia L, Wang S, Shi T. Topology optimization of hierarchical lattice structures with substructuring. Comput Methods Appl Mech Eng 2019;345:602–17. https://doi.org/10.1016/j.cma.2018.11.003.

15

[15] Wu J, Clausen A, Sigmund O. Minimum compliance topology optimization of shell – infill composites for additive manufacturing. Comput Methods Appl Mech Eng 2017;326(Supplement C):358–75. https://doi.org/10.1016/j. cma.2017.08.018. [16] Wu J, Aage N, Westermann R, Sigmund O. Infill optimization for additive manufacturing – approaching bone-like porous structures. IEEE Trans Visual Comput Graph 2018;24(2):1127–40. https://doi.org/10.1109/ TVCG.2017.2655523. [17] Dapogny C, Estevez R, Faure A, Michailidis G. Shape and topology optimization considering anisotropic features induced by additive manufacturing processes. Comput Methods Appl Mech Eng 2019;344:626–65. https://doi.org/10.1016/j. cma.2018.09.036. [18] Wang Y, Zhang L, Daynes S, Zhang H, Feih S, Wang M. Design of graded lattice structure with optimized mesostructures for additive manufacturing. Mater Des 2018;142:114–23. https://doi.org/10.1016/j.matdes.2018.01.011. [19] Fries T-P, Belytschko T. The extended/generalized finite element method: an overview of the method and its applications. Int J Numer Meth Eng 2010;84 (3):253–304. https://doi.org/10.1002/nme.2914. [20] Wei P, Wang M, Xing X. A study on X-FEM in continuum structural optimization using a level set model. Comput Aided Des 2010;42(8):708–19. https://doi.org/10.1016/j.cad.2009.12.001. [21] Van Miegroet L. Generalized shape optimization using XFEM and level set description Ph.D. thesis. University of Liege, Aerospace and Mechanical Engineering Department; 2012. [22] Abdi M, Ashcroft I, Wildman R. High resolution topology design with IsoXFEM. In: Proceedings of the 25th annual international solid freeform fabrication symposium. p. 1288–303. [23] Kreissl S, Maute K. Levelset based fluid topology optimization using the extended finite element method. Struct Multidiscip Optim 2012;46 (3):311–26. https://doi.org/10.1007/s00158-012-0782-8. [24] van Dijk N, Maute K, Langelaar M, van Keulen F. Level-set methods for structural topology optimization: a review. Struct Multidiscip Optim 2013;48 (3):437–72. https://doi.org/10.1007/s00158-013-0912-y. [25] Daux C, Moës N, Dolbow J, Sukumar N, Belytschko T. Arbitrary branched and intersecting cracks with the extended finite element method. Int J Numer Meth Eng 2000;48(12):1741–60. https://doi.org/10.1002/1097-0207 (20000830)48:12<1741::AID-NME956>3.0.CO;2-L. [26] Belytschko T, Parimi C, Moës N, Sukumar N, Usui S. Structured extended finite element methods for solids defined by implicit surfaces. Int J Numer Meth Eng 2003;56(4):609–35. https://doi.org/10.1002/nme.686. [27] Bathe K. Finite element procedures. 2nd ed. Englewood Cliffs, NJ: PrenticeHall; 1996. [28] Jansen M. Explicit level set and density methods for topology optimization with equivalent minimum length scale constraints. Struct Multidiscip Optim 2019;59(5):1775–88. https://doi.org/10.1007/s00158-018-2162-5. [29] Makhija D, Maute K. Numerical instabilities in level set topology optimization with the extended finite element method. Struct Multidiscip Optim 2014;49 (2):185–97. https://doi.org/10.1007/s00158-013-0982-x. [30] Zohdi T, Wriggers P. An introduction to computational micromechanics. Lecture notes in applied and computational mechanics, 1st ed., vol. 20. Verlag Berlin Heidelberg: Springer; 2005. https://doi.org/10.1007/978-3-54032360-0. [31] Hassani B, Hinton E. A review of homogenization and topology optimization I—homogenization theory for media with periodic structure. Comput Struct 1998;69(6):707–17. https://doi.org/10.1016/S0045-7949(98) 00131-X. [32] Geuzaine C, Remacle J-F. Gmsh: a 3-D finite element mesh generator with built-in pre- and post-processing facilities. Int J Numer Meth Eng 2009;79 (11):1309–31. https://doi.org/10.1002/nme.2579. [33] Bendsøe M, Sigmund O. Topology optimization: Theory, Methods and Applications. 2nd ed. Berlin: Springer; 2004. [34] Bourdin B. Filters in topology optimization. Int J Numer Meth Eng 2001;50 (9):2143–58. https://doi.org/10.1002/nme.116. [35] Bruns T, Tortorelli D. Topology optimization of non-linear elastic structures and compliant mechanisms. Comput Methods Appl Mech Eng 2001;190(26– 27):3443–59. https://doi.org/10.1016/S0045-7825(00)00278-4. [36] Wang F, Lazarov B, Sigmund O. On projection methods, convergence and robust formulations in topology optimization. Struct Multidiscip Optim 2011;43:767–84. https://doi.org/10.1007/s00158-010-0602-y. [37] Wang M, Wang X. Color level sets: a multi-phase method for structural topology optimization with multiple materials. Comput Methods Appl Mech Eng 2004;193(6):469–96. https://doi.org/10.1016/j.cma.2003.10.008. [38] Nocedal J, Wright S. Numerical optimization. New York, USA: Springer; 1999. [39] Svanberg K. The method of moving asymptotes – a new method for structural optimization. Int J Numer Meth Eng 1987;24:359–73. https://doi.org/10.1002/ nme.1620240207. [40] Allaire G, Jouve F, Toader A. Structural optimization using sensitivity analysis and a level-set method. J Comput Phys 2004;194(1):363–93. https://doi.org/ 10.1016/j.jcp.2003.09.032. [41] Zhou M, Lazarov B, Wang F, Sigmund O. Minimum length scale in topology optimization by geometric constraints. Comput Methods Appl Mech Eng 2015;293:266–82. https://doi.org/10.1016/j.cma.2015.05.003. [42] Hollister S, Kikuchi N. A comparison of homogenization and standard mechanics analyses for periodic porous composites. Comput Mech 1992;10 (2):73–95. https://doi.org/10.1007/BF00369853.

16

M. Jansen, O. Pierard / Computers and Structures 231 (2020) 106205

[43] Sethian J. Fast marching methods and level set methods for propagating interfaces. Tech rep, Department of Mathematics, University of California, Berkeley, California; 1998. [44] Schlömer N et al. nschloe/meshio: 2.3.7; May 2019. https://doi.org/10.5281/ zenodo.2702446.

[45] Clausen A. Topology optimization for additive manufacturing, Ph.D. thesis, DTU Department of Mechanical Engineering; 2016. [46] GrabCAD, GE jet engine bracket challenge, [ONLINE] Available at: https://grabcad.com/challenges/ge-jet-engine-bracket-challenge; 2012. .