- Email: [email protected]
VCUT level set method for topology optimization of functionally graded cellular structures
Available online at www.sciencedirect.com
ScienceDirect Comput. Methods Appl. Mech. Engrg. 354 (2019) 487–505 www.elsevier.com/locate/cma
VCUT level set method for topology optimization of functionally graded cellular structures Hongming Zonga , Hui Liua,b , Qingping Maa , Ye Tiana , Mingdong Zhouc , Michael Yu Wanga ,∗ a
Department of Mechanical and Aerospace Engineering, Hong Kong University of Science and Technology, Kowloon, Hong Kong Department of Engineering Mechanics, School of Civil Engineering, Wuhan University, Wuhan, Hubei, People’s Republic of China c Shanghai Key Laboratory of Digital Manufacture for Thin-walled Structures, School of Mechanical Engineering, Shanghai Jiao Tong University, Shanghai, People’s Republic of China b
Received 15 November 2018; received in revised form 3 May 2019; accepted 20 May 2019 Available online 31 May 2019
Abstract This paper presents a variable cutting (VCUT) level set based topology optimization method to design functionally graded cellular structures (FGCS). A variable and continuous cutting function by interpolating with a set of height variables is proposed to generate functionally graded cellular structures, which offers a novel tool to optimize the macroscopic graded pattern. Due to the continuity of the cutting function, perfect geometric connections between adjacent cells are guaranteed without imposing extra constraints in the optimization. In addition, a solid covering skin can be easily attached to the FGCS by means of the variable cutting function. Three FGCS design problems, including graded density control, compliance minimization and layered cellular structure design, are presented to demonstrate the effectiveness and applicability of the proposed method. The optimized 2D and 3D macroscopic FGCS designs exhibit well-performing structural layout with fully connected micro-scale geometries. c 2019 Elsevier B.V. All rights reserved. ⃝ Keywords: Topology optimization; Functionally graded cellular structures; Level set method; Microstructure connectivity
1. Introduction The cellular structures composed of downscale microstructures are known for superior performance with relatively low densities [1]. By tailoring the configurations in the micro scale, excellent multifunctional responses can be achieved, such as stress resistance, heat exchange and shock absorption [2]. Besides, recent advances in additive manufacturing (AM) provide additional design freedom in fabricating these structures with complex and multiscale geometric features [3,4], which also stimulates the development of cellular structures. In particular, the functionally graded cellular structures [5] with spatially varying properties have gained increased popularity in engineering over the past decades. Compared with the periodically repeated cellular structures, progressively changing the geometric feature shapes or sizes of the embedded microstructures [6] allows FGCS ∗
Corresponding author. E-mail address: [email protected] (M.Y. Wang).
https://doi.org/10.1016/j.cma.2019.05.029 c 2019 Elsevier B.V. All rights reserved. 0045-7825/⃝
488
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
to flexibly adjust and customize the layout according to the loading and boundary conditions in practice. For example, the gradual variation in composition can reduce the stress concentration near the interface between different phases [7] and thus can be used to relieve the “stress shielding” phenomenon for orthopedic implants. A decreased density distribution of honeycomb cellular structure in the crushing direction can also be adopted to enhance the energy absorption at early crushing stages [8]. Other applications can also be found in for instance thermal insulation [9] and acoustical wave attenuation [10]. As a result, FGCS have attracted considerable attentions from the aerospace, automobile and medical industries. However, the FGCS common in biological systems ranging from the bird beaks, bones to plant stems cannot be easily designed due to the difficulties in modeling the configurations that typically span over several length scales and in determining a suitable graded pattern at the meantime. A widely used method for realizing this purpose is to first characterize the FGCS via explicit parameters and then run parametric optimization to obtain a graded distribution [11,12]. However, it is not easy to use an explicit parametric technique to represent microstructures involving general types of internal features and geometric variations. Therefore, many newly designed artificial materials [13,14] cannot be adopted by the FGCS, which limits the structural performance for multiple functionalities. Alternatively, topology optimization is a suitable design tool for FGCS because of its capability in generating innovative structural and material layouts with optimized performance. Interestingly, the very early homogenization method [15] is developed based on a similar concept of FGCS, where a hole is used in each cell and parametric optimization is conducted to determine the sizes of the holes. After three decades of development, various methods, e.g. the solid isotropic material with penalization method (SIMP) [16,17], the evolutionary structural optimization [18] and the level set method [19,20] have been proposed. Particularly, Alzahrani et al. [21] proposed a mapping technology to design FGCS, where the intermediate densities in SIMP are mapped to lattice unit cells with variable volume fractions. Schury et al. [22] used free material optimization method to produce manufacturable graded structures. Li et al. [23] developed a material-structure integrated topology optimization scheme for producing layer-wise FGCS. Daynes et al. [24] presented a novel biomimetic methodology for designing spatially graded lattice structures. Recently, Li et al. [25] proposed a concurrent topology optimization method that can generate cellular structures with multiple patches of microstructures rationally distributed in the structure, and Li et al. [26] presented a gyroid-based FGCS design method for optimizing lightweight lattice structures with uniformly varying densities in the design space. When designing FGCS, it is crucial to guarantee the connectivity between graded cells. If homogenization method [27] is applied directly without any geometric constraints, the scale separation assumption will result in unconnected regions, which cannot bear any load and will cause difficulties in additive manufacturing. Regarding this problem, Zhou and Li [28] made an in-depth discussion on this issue and developed three methods to guarantee the connectivity. Radman et al. [29] applied the filtering technique in every single microstructure design and preserved the connectivity in a sequential way. Though these implementations can relieve the disconnection problem to some extent, they can be hardly extended for complex graded patterns determined by the macroscopic loading and boundary conditions. Wang et al. [30] and Carmer et al. [31] used metamorphosis technology and interpolation scheme to generate connectable graded microstructures for macroscopic design. However, mismatch is observed in the obtained FGCS, which may cause stress concentrations and degrade the final structural performance. This paper develops a novel variable cutting level set based topology optimization method for FGCS design. A basic periodic level set function in which the microstructures are embedded is first constructed and distributed in a replicated way. Then a variable cutting function is introduced for generating the cellular structure with various microstructures and optimizing the graded pattern according to the requirement. Due to the periodicity of the basic level set function and the connectivity of the variable cutting function, the obtained FGCS are guaranteed to have at least C0 continuity. The proposed approach is demonstrated by several 2D and 3D FGCS design examples of graded density control, compliance minimization and layered structure design problems, respectively. 2. Modeling on VCUT level set method As shown in Fig. 1(a), the considered design domain with length L and height W for FGCS design is discretized into M unit cells with m nodes in the macro scale. For each unit cell c, it is further discretized in the micro scale as shown in Fig. 1(b). Based on such a discretization, the VCUT level set model is proposed to represent the graded cellular structures with two components: a basic level set function for the geometric details in the micro scale and a variable cutting function for the adjustable graded pattern in the macro scale.
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
489
Fig. 1. Design domain discretization for FGCS: (a) macro scale and (b) micro scale.
Fig. 2. General level set method for representing dynamic structural interfaces. (a) The level set function φ(x, t) and zero cutting plane, and (b) structural boundaries represented by the level set function.
2.1. Basic level set function In the general level set method, the dynamic structural interfaces are implicitly expressed by the zero level set of a higher-dimensional level set function φ(x, t) [19,20], as illustrated in Fig. 2. The geometry model can be represented by: ⎧ ⎨φ(x, t) > 0 ∀x ∈ Ω \∂Ω φ(x, t) = 0 ∀x ∈ ∂Ω (1) ⎩ φ(x, t) < 0 ∀x ∈ D\Ω , where x ∈ D is any point in the fixed design domain D and ∂Ω are the boundaries of the solid domain Ω . Considering the unit cell in a cellular structure, the level set method offers a flexible and convenient way for representing complex internal micro features. Theoretically, any type of microstructures can be digitally built based on the method and some 3D unit cells are illustrated in Fig. 3. One can use the signed distance function (Fig. 3(a– d)) [32,33] or the minimal surface function (Fig. 3(e)) [32] as the basic level set function φp (x) for a unit cell. In general, the signed distance function within the unit cell Dc can be calculated by: ⎧ ⎨+ min ∥x − x b ∥ ∀x ∈ Ωp x b ∈∂Ωp φp (x) = (2) ⎩− min ∥x − x b ∥ ∀x ∈ Dc \Ωp , x b ∈∂Ωp
where ∂Ωp denotes the boundaries of the microstructure. In addition, the unit cell with corresponding level set function φp (x) can also be obtained by using the inverse homogenization method within the level set framework. Some of the optimized unit cells with desired properties like negative Poisson’s ratio are given in Fig. 3(c, d) [33]. With a set of basic level set functions, extra operations such as the basic Boolean operations including intersection, union and difference between any two cells can be easily realized by a series of min/max operations
490
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
Fig. 3. Unit cells represented by level set method [32,33], including (a) octet truss, (b) carbon truss, (c, d) optimized microstructures with negative Poisson’s ratio, (e) Schwarz P surface unit cell and (f) artificial unit cell by union of a Schwarz P surface unit cell with a body-centered cubic unit cell.
on the level set functions. Fig. 3(f) is a unit cell generated by union of a Schwarz P surface unit cell with a body-centered cubic unit cell. These operations provide great design freedom for engineers to customize the unit cells according to practical design requirements. When the basic level set function φp (x) is determined, it could be mapped into macro elements in a replicated way, which will form a global level set function φg (x). Note that the φp (x) used in the VCUT level set method should possess geometric periodicity, which means that φp (x) on the opposite boundaries of the unit cell should be identical. This geometric periodicity may restrict the design space of its properties to some extent. However, it is still flexible enough for the unit cell to achieve various functionalities [23,29,30], including some extreme [14] or extraordinary properties [13]. 2.2. Variable cutting function As illustrated in Fig. 2, it is easy to notice that the boundaries are determined simultaneously by the level set function and the cutting function. Typically, the zero cutting plane is selected and fixed during the optimization. By evolving the level set function, the boundaries keep moving till an optimized structure is obtained. Alternatively, with a predefined level set function, optimization can also be performed by changing the cutting function globally [34]. A similar concept was used for graded cellular structure design in [30], where a series of unit cells were generated by adjusting the heights of the cutting planes as plotted in Fig. 4. The graded cellular structures were optimized by determining the heights in each macro element. However, as indicated in their work, mismatch on the common boundaries of neighboring cells cannot be avoided based on this modeling scheme. A simple illustration of this problem is given in Fig. 5, where the cellular structure is made up of 2 by 2 macro elements. Because the cutting planes are element-wisely distributed, there will be an offset between two neighboring cutting planes if the heights are different. As a result, mismatches may occur at the boundaries of adjacent cells. In order to achieve a natural geometric connection without the above mismatch issue, a variable and continuous cutting function is introduced in this paper. As shown in Fig. 1, a set of height variables are defined on the nodes of the entire macro mesh like the nodal density in [35], and they are represented by a time-dependent global vector: [ ]T H(t) = h 1 (t) h 2 (t) h 3 (t) . . . h m (t) , (3) where m denotes the total number of the macro mesh nodes. The height vector of a macro mesh element, for example c, is a subset of the global height vector: [ ]T hc (t) = h c1 (t) h c2 (t) h c3 (t) h c4 (t) = Sc H(t), (4) with Sc denoting its selection matrix in the common notation. The cutting function within the unit cell can be formulated by interpolating the hc (t) using the common shape functions. For convenience, a shape function vector N c (x) is defined as: [ ]T (5) N c (x) = N1c (x) N2c (x) N3c (x) N4c (x) , and the cutting function controlled by the height variables is given as: Cc (x, t) = N Tc (x)hc (t) = N Tc (x)Sc H(t).
(6)
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
491
Fig. 4. Signed distance function cut by planes at different heights. (a) The basic level set function generated on ∂Ωp in (b), and (b–e) generated microstructures.
Fig. 5. Illustration on the generation of mismatch in [30]. (a) The global level set function φg (x) and element-wisely distributed cutting planes, and (b) the generated cellular structure and mismatches at the boundaries of adjacent cells.
With the variable cutting function, the level set function for the geometry model of unit cell c becomes: ⎧ ⎨φc (x, t) = φp (x) − Cc (x, t) > 0 ∀x ∈ Ωc \∂Ωc φc (x, t) = φp (x) − Cc (x, t) = 0 ∀x ∈ ∂Ωc ⎩ φc (x, t) = φp (x) − Cc (x, t) < 0 ∀x ∈ Dc \Ωc .
(7)
It can be found that the geometric features in the macro elements are determined by the cutting function, or essentially by the cutting heights H(t). The cutting effect for a single cell is illustrated in Fig. 6. For any macro element, if the corresponding hc (t) is high, the element becomes a void element; if it is low, the element becomes a solid element; if it lays in between, the microstructure will exhibit geometric variations from its basic geometry, as shown in Fig. 6(b). In addition, since the heights do not have to be the same, the obtained unit cell can exhibit directional stiffness according to the external loading and boundary conditions and thus possesses a larger design space compared with the unit cells in Fig. 4. It should be noted that, in the current definition of the variable cutting function by Eq. (6), Cc (x, t) is mainly used to control the cell gradient in the macro scale and the potential topological variation of a cell in the micro scale
492
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
Fig. 6. Illustration on the variable cutting process. (a) The basic level set function φp (x) and the variable cutting function Cc (x, t) controlled by heights defined on four nodes, (b) the generated microstructure.
is determined by φp (x). However, if a more flexible function form is chosen, then the cutting function may be able to change the topology in the micro scale with the topological derivatives [36] provided. In this way, the structural topology in the micro scale becomes more flexible, which can help to enlarge the design space and improve the structural performance further. By assembling the cutting functions of unit cells together, a global cutting function Cg (x, t) with at least C0 continuity is formed. Combined with previous global level set function φg (x), the cellular structure can be modeled by: ⎧ ⎨φ(x, t) = φg (x) − Cg (x, t) > 0 ∀x ∈ Ω \∂Ω φ(x, t) = φg (x) − Cg (x, t) = 0 ∀x ∈ ∂Ω (8) ⎩ φ(x, t) = φg (x) − Cg (x, t) < 0 ∀x ∈ D\Ω , where the global boundaries ∂Ω : = {x|x ∈ ∂Ωc , c = 1, 2, . . . , M} are the assembly of cellular boundaries ∂Ωc . Since Cg (x, t) is variable and can be adjusted by changing H(t), cellular structures can yield variations in geometric features and thus graded patterns can be expected. Meanwhile, since neighboring cells share the same pair of height variables on common boundaries, the global cutting function Cg (x, t) has at least C0 continuity across the cells. Considering that the global level set function φg (x) also has C0 continuity because of the periodicity of φp (x), the resulted level set function φ(x, t) is C0 continuous. Therefore, the VCUT method can naturally guarantee the full geometric continuity for microstructures in macro elements without imposing any extra constraints. As an example, a 2 by 2 cellular structure is illustrated in Fig. 7, where the heights are set to be different and thus a graded distribution is observed. No geometric mismatches are observed at the connections between any two neighboring cells. Further, if higher-order interpolating functions are applied, e.g., using p-type shape functions, the cutting function Cg (x, t) will exhibit a geometric continuity higher than C0 . As a result, the cellular structure represented by the level set function can have higher-order continuity, e.g. C1 continuity, at the geometric boundaries between the neighboring cells. 3. Topology optimization using the VCUT level set method With the proposed VCUT model, the topology optimization problem is formulated as: min. J (u, Ω ) Ω
s.t.
⎧ ⎪ ⎪ ⎪ ⎨
a(u, w, ∫ Ω ) = l(w, Ω ) g(Ω ) =
∀w ∈ U (9)
dx − VD f v ≤ 0 Ω
⎪ ⎪ ⎪ ⎩
u = u0 Eε(u) · n = τ
in Γu in Γτ ,
where J(u, Ω ) is the objective functional, which can be defined in the following form without loss of generality: ∫ J (u, Ω ) = F(u)dx, (10) Ω
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
493
Fig. 7. Illustration on generation of graded cellular structure based on the VCUT level set method. (a) The global level set function φg (x) and the global cutting function Cg (x, t), (b) graded and connected cellular structure.
where g(Ω ) is the volume constraint, VD is the volume of the design domain and fv is the allowed maximum volume fraction; U is the kinematically admissible displacement space; a(u, w, Ω ) and l(w, Ω ) are the terms of the state equation of the linear elastic system: ∫ εT (u)Eε(w)dx, (11) a(u, w, Ω ) = Ω ∫ ∫ l(w, Ω ) = f T wdx + τ T wdx. (12) Ω
Γτ
The derivative of the objective function J(u, Ω ) is related to the variation of the solid domain Ω . In the VCUT method, since the solid domain is controlled by the height variable H(t), the relationship between the derivative of J(u, Ω ) and the variation of H(t) is derived in this section, which is used to optimize the graded pattern of the cellular structure. 3.1. Sensitivity analysis The Lagrange function of the optimization problem (9) is formulated as: L = J + a − l + λg,
(13)
where λ is the Lagrange multiplier. Following the two lemmas listed in [37], the shape derivative of Eq. (13) with respect to the pseudo time t is calculated by: L˙ = J˙ + a˙ − l˙ + λg, ˙
(14)
where: ∫ F ′ (u)dx + F(u)Vn dx, Ω ∂Ω ∫ ∫ ∫ a(u, ˙ w, Ω ) = εT (u′ )Eε(w)dx + εT (u)Eε(w′ )dx + εT (u)Eε(w)Vn dx, Ω Ω ∂Ω ∫ ∫ ∫ ∫ [ ] ( T )T T ′ T T ′ ˙l(w, Ω ) = f w dx + f wVn dx + τ w dx + ∇τ w n + τ T wk Vn dx, ∂Ω Γτ Γτ ∫ Ω g(Ω ˙ )= Vn dx, J˙(u, Ω ) =
∫
∂Ω
(15) (16) (17) (18)
494
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
where Vn is the normal velocity of the boundaries ∂Ω ; k is twice the mean curvature of the boundary defined as k = div n. Since the body force is assumed zero on the boundary ∂Ω and Γτ is fixed, the derivative of the virtual work can be simplified as: ∫ ∫ ˙ l(w, Ω) = f T w ′ dx + τ T w ′ dx. (19) Ω
Γτ
By Substituting all these terms into Eq. (14), the derivative of the Lagrange function can be reformulated as: {∫ } ∫ ∫ T ′ T ′ T ′ ˙L = ε (u)Eε(w )dx − f w dx − τ w dx Ω Ω {∫ } Γτ ∫ (20) + F ′ (u)dx + εT (u′ )Eε(w)dx Ω∫ {∫Ω } ∫ + F(u)Vn dx + εT (u)Eε(w)Vn dx + λ Vn dx . ∂Ω
∂Ω
∂Ω
Since a(u, w, Ω ) = l(w, Ω ) ∀w ∈ U and w ′ ∈ U, the first part of the derivative is canceled by: ∫ ∫ ∫ T ′ T ′ τ T w′ dx, f w dx + ε (u)Eε(w )dx = By solving the following adjoint equation: ∫ ∫ F ′ (u)dx = − εT (u′ )Eε(w)dx, Ω
(21)
Γτ
Ω
Ω
∀u′ ∈ U,
(22)
Ω
the derivative of the Lagrange equation can be obtained as: ∫ ∫ [ ] T ˙L = F(u) + ε (u)Eε(w) + λ Vn dx = G(u, w)Vn dx. ∂Ω
Moreover, by discretizing the global boundaries ∂Ω into cellular boundaries ∂Ωc , the derivative becomes: ∫ M ∫ ∑ c ˙L = G(u, w)Vn dx = G(u, w)Vnc dx, ∑ M c=1 ∂Ωc
(23)
∂Ω
c=1
(24)
∂Ωc
where Vnc is the normal velocity that drives ∂Ωc to move in the cell. 3.2. Update scheme for height variables When evolving the structural frontiers ∂Ωc in unit cells, the following Hamilton–Jacobi equation must be satisfied: dH(t) ∂φc (x, t) dφc (x, t) = − ∥∇φc (x, t)∥ Vnc = −N Tc (x)Sc − ∥∇φc (x, t)∥ Vnc = 0, ∀x ∈ ∂Ωc . (25) dt ∂t dt The relationship between the normal velocity Vnc and the design variables H(t) can be described as: 1 dH(t) N T (x)Sc , ∀x ∈ ∂Ωc . ∥∇φc (x, t)∥ c dt By substituting the normal velocity of Eq. (26) into Eq. (24), the derivative becomes: M ∫ ∑ 1 dH(t) G(u, w) N Tc (x)Sc L˙ = − dx ∥∇φ (x, t)∥ dt c c=1 ∂Ωc ] ∫ M [ dH(t) T ∑ T G(u, w) =− Sc N c (x)dx . dt ∂Ωc ∥∇φc (x, t)∥ c=1 Vnc = −
(26)
(27)
To guarantee the steepest descent searching direction, the derivative of the design variable H(t) is defined as: ] ∫ M [ dH(t) ∑ T G(u, w) = Sc N c (x)dx . (28) dt ∂Ωc ∥∇φc (x, t)∥ c=1
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
495
Fig. 8. Workflow of the VCUT level set method.
Fig. 9. Basic level set functions used in this paper, including (a) φp1 (x) and its zero level set (b), (c) φp2 (x) and its zero level set (d), and (e) zero level set of φp3 (x).
and the height variables are updated by: dH(ti ) ∆t, (29) dt where ti is the ith time step and ∆t is the time step size. The workflow of the proposed VCUT level set method is presented in Fig. 8. H(ti+1 ) = H(ti ) +
4. Numerical implementation In this paper, three basic level set functions are utilized for demonstration, including two of φp1 (x) in Fig. 9(a) and φp2 (x) in Fig. 9(c) for generating 2D unit cells and one of φp3 (x) in Fig. 9(e) for generating 3D unit cells. The φp2 (x) is constructed by the Schwarz P minimal surface formulation [32] as: φp2 (x) = cos(2π x) + cos(2π y),
∀x = (x, y) ∈ Dc .
(30)
The φp1 (x) and φp3 (x) are signed distance functions built based on Eq. (2), where truss central axes are chosen as ∂Ωp . For example, two diagonals of the square in Fig. 4(b) are taken as the boundaries ∂Ωp for constructing φp1 (x). Similarly, four diagonals of the cube are used for building φp3 (x). Then, All the basic level set functions are normalized to −1 and 1 to guarantee the generality of the numerical program. During the optimization, the state equation in (9) needs to be solved to obtain the displacement field u. Since the geometric features in the cellular structure typically span over several length orders, it is difficult to solve the state equation directly by the ordinary finite element method. In this study, an efficient multiscale finite element method (MsFEM) [38] is applied for mechanical analysis. In the MsFEM, the multiscale shape function is constructed to capture the micro-scale features in unit cells and thus the state equation can be solved on the macro scale. In addition, the construction process is parallelized to ensure a good computational efficiency.
496
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
In Eq. (28), the boundary integral is calculated to determine the derivatives. To avoid the time-consuming remeshing process when considering the exact boundaries during the analysis, the problem (9) is reformulated with the Heaviside function and Dirac function δ(φ) as in [19]. Under this formulation, the derivative of Eq. (28) becomes: ] ∫ M [ dH(t) ∑ T = Sc G(u, w)N c (x)δ(φc )dx (31) dt Dc c=1 where δ(φ) is chosen to be the same as that in [19]. In addition, the Lagrange multiplier λ is selected and updated using the scheme in [39]. To illustrate the generality of the proposed VCUT level set method, three design problems, including graded density distribution, compliance minimization and layered structure design, are studied. For different problems, the optimization setup (9) or the graded pattern may be different. Hence, the formulation of the derivative in Eq. (31) should be modified accordingly. Graded density control The first problem is to obtain FGCS with a given graded density distribution. It is formulated as an unconstrained optimization problem and the objective function can be defined as: M
J (Ω ) =
1∑ (Vc − f c∗ VDc )2 , 2 c=1
(32)
where Vc is the volume of the microstructure in unit cell c, f c∗ indicates the predefined target volume ratio of cell c and VDc is the volume of the unit cell domain. The height derivative under this objective function can be easily derived as: ] ∫ M [ ) dH(t) ∑ ( = Vc − f c∗ VDc STc N c (x)δ(φc )dx . (33) dt Dc c=1 Compliance minimization problem For the compliance minimization problem, the objective function is: ∫ 1 εT (u)Eε(u)dx. J (u, Ω ) = 2 Ω The corresponding adjoint equation is: ∫ ∫ εT (u′ )Eε(u)dx, ∀u′ ∈ U, εT (u′ )Eε(w)dx = −
(34)
(35)
Ω
Ω
and thus the adjoint displacement is simply w = -u. The derivative can then be calculated by: ] } ∫ [ M { dH(t) ∑ T 1 − εT (u)Eε(u) + λ N c (x)δ(φc )dx . = Sc dt 2 Dc c=1
(36)
Layered cellular structure design Layered cellular structure is one typical type of FGCS, where the geometric features in microstructures vary in a multilayered manner. Here, the layered structures are designed for a minimized structural compliance. Different from the previous example, extra operations must be applied in order to obtain layered distribution. Similar to the pattern repetition constraints [40], the obtained derivatives and height variables in the same layer are averaged in each iteration step. Suppose there are n nodes in each layer, the averaged node height in the l-th layer can be calculated by: hl =
1 n
ln ∑
hi .
(37)
i=(l−1)n+1
This process is illustrated in Fig. 10. Since the structure is still modeled based on the VCUT level set method, this operation will not influence the connection of the generated cellular structure.
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
497
Fig. 10. Process for generating layered cellular structure within the VCUT method. (a) The height variables before the operation, and (b) the height variables after the operation.
Fig. 11. (a) The predefined density distribution, and (b) FGCS obtained by the VCUT method.
5. Numerical examples In this section, several 2D and 3D numerical examples on the three design problems indicated in Section 4 are used to demonstrate the proposed VCUT level set method for FGCS design. For all examples, the Young’s moduli for the solid material is Es = 1 and for the artificial weak material is Ev = 10−6 . The latter is used to avoid numerical singularity when solving the state equation. The Poisson’s ratios for both materials are 0.3. Unless otherwise specified, the optimization terminates when the relative change of 10 successive objective functions is below 0.1%, or the maximum 300 iteration steps is reached. For 2D problems, the same design domain D with W = 200 and L = 600 is used, which is discretized into 60 × 20 unit cells and each unit cell is further discretized by 10 × 10 4-node elements. For 3D problems, the same design domain D with the dimension of 128 × 32 × 64 is employed, which is discretized into 16 × 4 × 8 unit cells and each unit cell is further discretized into 8 × 8 × 8 8-node elements. 5.1. Graded density control In this problem, since the objective function is evaluated based on the absolute differences between the obtained design and the target design, a different terminate condition is used that the objective function is less than 0.01. The target density distribution ranging from 0.2 to 0.8 is given in Fig. 11(a) and the φp1 (x) is used for illustration. After 77 iteration steps, the objective function reduces to 9.56 × 10−3 and the produced design is shown in Fig. 11(b). Note that, although an iterative optimization process is carried out here, the computational cost is very low because no mechanical analysis is involved. The resulted design exhibits an obvious graded variation in both geometric features and unit cell densities. In addition, the magnified plot indicates that each unit cell is connected well as expected. With this simple example, the VCUT level set method is proved to be capable of controlling the density distribution of cellular structure, which provides designers with a convenient and efficient way to adjust the graded pattern in practical usage.
498
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
Fig. 12. Loading and boundary conditions for 2D cantilever beam design.
Fig. 13. Iteration history of 2D cantilever beam design based on φp1 (x) using the VCUT level set method, and (a–d) the intermediated designs during the optimization.
5.2. Compliance minimization Case 1 2D cantilever beam A benchmark cantilever beam design problem is studied in this example, where the loading and boundary conditions are illustrated in Fig. 12. The beam is fixed at left-hand side and loaded vertically at the right-bottom point with F = 1. The volume constraint is set to fv = 0.5. Fig. 13 shows the iteration history of the design process based on φp1 (x). After 157 iteration steps, the optimization converges to the solution as shown in Fig. 14(a) with the objective function J = 1.86 × 102 . With the same problem setup, the optimization with the other basic level set function φp2 (x) converges after 109 iteration steps with the objective function J = 1.95 × 102 . The obtained FGCS is plotted in Fig. 14(b). Controlled by the variable cutting functions, for instance in Fig. 15 for the case of φp1 (x), two solutions present a similar and rational graded pattern to bear the vertical load at the corner. For areas that contribute little to bearing the load, the cutting function moves beyond the maximum value of φp (x) and generates the void regions. For the load path areas, the cutting function moves below the minimum value of φp (x) and generates the solid region. For areas in between, microstructures with low densities are generated to support the two main solid frames, similar to the solutions obtained by the classical SIMP method with penalty factor of p = 1. In the classical SIMP method, the effective properties of intermediate densities with penalty factor of p = 1 are better than those with penalty factor of p = 3. As a result, the optimized results with p = 1 typically have a great percentage of gray regions. In our case, the cells with porous microstructures play a role of the elements with intermediate densities. As indicated in [30], this kind of graded cellular structures can provide a better mechanical performance. In addition, the geometric features of the microstructures are perfectly connected when the cells change gradually, as shown in the magnified plots of Fig. 14. Note that, in the transition areas where the cells fade to void, the cutting function will intersect with the maximum value of the basic level set function. As a result, small burrs appear on the boundaries. To achieve smooth boundaries, a solid covering skin can be attached to the FGCS by means of the global cutting function, which can also act as a post-processing to facilitate the real manufacturing [4].
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
499
Fig. 14. Obtained FGCS based on (a) φp1 (x) and (b) φp2 (x) using the VCUT level set method, and (c) and (d) the FGCS with the covering skins based on (a) and (b).
Fig. 15. Variable cutting function for the case of φp1 (x).
Since the cutting function is a continuous function and the burrs only occur when Cg (x) cuts across the maximum value of φp (x), the boundaries of the FGCS in the macroscopic scale ∂ΩFGCS can be determined by: Cg (x) − m p = 0,
∀x ∈ ∂ΩFGCS ,
(38)
where mp denotes the maximum value of φp (x). Based on the ∂ΩFGCS , a new signed distance function φFGCS can be built according to Eq. (2). Considering the boundaries of the design domain, another signed distance function φD can also be built. Finally, the function for generating the whole covering skin can be calculated by a min operation as: φskin (x) = min(φFGCS , φD ),
(39)
and the new level set function φnew (x) to represent FGCS with covering skins can be represented by: φnew (x) = max(min(φskin , −(φskin − θ)), φ),
(40)
where θ is the parameter to control the thickness of the skin. By choosing θ = 3 for illustration, φnew (x) based on the solutions in Fig. 14(a, b) are calculated and then substituted into Eq. (1) to determine the solid areas. Two final FGCS with covering skins are illustrated in Fig. 14(c, d). Case 2 2D simply supported beam Fig. 16 defines the loading and boundary conditions for the 2D simply supported beam design. The left-bottom point of the beam is fixed and the vertical degree of freedom at the right-bottom point is constrained. A vertical load
500
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
Fig. 16. Loading and boundary conditions for 2D simply supported beam design.
Fig. 17. Obtained FGCS based on (a) φp1 (x) and (b) φp2 (x) using the VCUT level set method, and (c) and (d) the FGCS with skins based on (a) and (b).
with F = 1 is applied at the center of the upper edge. The volume constraint is set to fv = 0.5. The optimization based on φp1 (x) converges to J = 26.4 after 114 iteration steps and that based on φp2 (x) ends up with J = 27.9 after 120 iteration steps. The obtained FGCS as well as those with covering skins are presented in Fig. 17. The two designs exhibit different topologies in the central area, where the cutting functions of the corresponding designs are compared in Fig. 18. Such a difference can be mainly attributed to the different intrinsic geometric features of the basic level set functions. Simple models consisting of two unit cells as shown in Fig. 19 are used for explanation. Since the height variables are defined on the nodes, when the height on node A is set to 1, the two cells lose connection around this node immediately. As a comparison, when the height on node B equals 1, the two cells are still connected at the common boundary and thus it is not easy to generate holes around the center line. In addition, two unit cells located around the load paths are extracted as shown in Fig. 17(a). Since the FGCS are symmetric with respect to the vertical center line, configurations of the two cells are mirrored with each other. Even with the same density value, they actually exhibit directional stiffness to follow the macroscopic loadings. Compared with the methods that characterize the unit cells with only density variables, the VCUT level set method has a larger design space and potential to further improve the mechanical performance of FGCS. Case 3 3D cantilever beam A 3D cantilever beam as shown in Fig. 20 is optimized based on the basic level set function of φp3 (x). The beam is fixed at the left face and loaded at the center of right-bottom edge with F = 100. The volume fraction is set to fv = 0.4. After 182 iteration steps, the optimization converges to J = 3.16 × 104 and the evolving history is plotted in Fig. 21. The final design is presented in Fig. 22 and that with a covering skin is shown in Fig. 23. The 3D microstructures are distributed in a graded way for bearing the load and those which do not contribute are wiped out. Since the thickness is relatively small, the result shows a similar configuration with that in Fig. 14. Nevertheless, since the
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
501
Fig. 18. Variable cutting function for the case of (a) φp1 (x) and (b) φp2 (x).
Fig. 19. Simple cellular structures consisting of two unit cells based on (a) φp1 (x) and (b) φp2 (x).
Fig. 20. Loading and boundary conditions for 3D cantilever beam design.
concentrated load is applied at the center, the variations in the thickness direction can still be observed. In addition, the connections between graded cells are verified through the magnified plots, where no mismatch is observed. 5.3. Layered cellular structure design Case 1 2D layered cellular structure Using the same loading and boundary conditions in Fig. 12, layered cellular structures are optimized based on the basic level set function φp1 (x). Fig. 24 gives the iteration history of the optimization, which converges after 151 steps with the objective function J = 2.57 × 102 . The final solution is presented in Fig. 25(a), where the upper and bottom layers end up with solid microstructures and the middle regions are filled with gradually varied microstructures. This graded distribution pattern is determined by the cutting function as illustrated in Fig. 25(b),
502
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
Fig. 21. Iteration history of 3D cantilever beam design.
Fig. 22. Optimized 3D cantilever beam using the VCUT level set method.
Fig. 23. Optimized 3D cantilever beam with covering skin.
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
503
Fig. 24. Iteration history of 2D layered cellular structure design.
Fig. 25. (a) Obtained 2D layered cellular structure based on φp1 (x) and (b) corresponding cutting function.
where the values at the two ends are below −1, which is the minimal value of φp1 (x), and then gradually increase to a value close to 1 in the middle. Case 2 3D layered cellular structure The same 3D cantilever beam problem in Fig. 20 is solved within the layered design scheme. The optimization terminates after 277 steps with J = 4.71 × 104 , and its iteration history is illustrated in Fig. 26. The obtained solution is shown in Fig. 27, which is similar to that in Fig. 25(a). Two solid plates are generated and graded cells are distributed in between to support them. It should be noted that, the compliance in the layered cellular structures is higher than that of Section 5.2 Case 3, because the operation in Eq. (37) essentially restricts the design space. Even though, according to the magnified plot in Fig. 27, it is proved that the operation has no influence on the microstructure connection. Geometric features in graded cells connect well on the common boundaries since the layered cellular structure is still modeled by the VCUT method and thus the connection is guaranteed. 6. Conclusions In this paper, a variable cutting level set topology optimization method is developed for designing functionally graded cellular structures. Controlled by the heights defined on macro element nodes, the continuous cutting function can be optimized to generate the graded pattern in the macro scale for the cellular structures. In addition, due to the continuity of the cutting function and the periodicity of the basic level set function, the obtained FGCS are ensured to have at least C0 continuity. Based on the method, three design problems, including graded density control, compliance minimization and layered cellular structure design, are studied. Several 2D and 3D numerical examples are given, and rational graded patterns as well as perfect connections among cells are observed, which demonstrates the effectiveness of the proposed design method. In this study, three predefined basic level set functions are used for illustrating the ability of the VCUT method in determining optimized graded patterns. In the future, a more general concurrent design method can be developed by combining the VCUT method with the inverse homogenization method. Given specific loading and boundary conditions, the inverse homogenization method can be used to generate the optimized basic level set function in the micro scale, and the VCUT method can then be adopted simultaneously to obtain the optimized graded pattern in the macro scale. By such a concurrent design method, further performance improvement can be expected. In
504
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
Fig. 26. Iteration history of 3D layered cellular structure design.
Fig. 27. Obtained 3D layered cellular structure using the VCUT level set method.
addition, only the bilinear or trilinear shape function is adopted as the interpolating function to form the cutting function in this paper. In the future, a higher-order interpolating function can be used to build a more flexible cutting function and to achieve higher-order continuity. Acknowledgment The authors acknowledge the support from National Natural Science Foundation of China (Grant No. 51705311). References [1] L.J. Gibson, M.F. Ashby, Cellular Solids: Structure and Properties, Cambridge University Press, 1999. [2] H.N.G. Wadley, Multifunctional periodic cellular metals, Phil. Trans. R. Soc. A 364 (1838) (2006) 31–68. [3] M.G. Rashed, M. Ashraf, R.A.W. Mines, et al., Metallic microlattice materials: a current state of the art on manufacturing, mechanical properties and applications, Mater. Des. 95 (2016) 518–533. [4] Y. Wang, L. Zhang, S. Daynes, et al., Design of graded lattice structure with optimized mesostructures for additive manufacturing, Mater. Des. 142 (2018) 114–123. [5] D. Mahmoud, M.A. Elbestawi, Lattice structures and functionally graded materials applications in additive manufacturing of orthopedic implants: a review, J. Manuf. Mater. Process. 1 (2) (2017) 13. [6] M. Niino, S. Maeda, Recent development status of functionally gradient materials, ISIJ Int. 30 (9) (1990) 699–703.
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
505
[7] K. Shah, I. ul Haq, A. Khan, et al., Parametric study of development of inconel-steel functionally graded materials by laser direct metal deposition, Mater. Des. 54 (2014) 531–538. [8] A. Ajdari, H. Nayeb-Hashemi, A. Vaziri, Dynamic crushing and energy absorption of regular, irregular and functionally graded cellular structures, Int. J. Solids Struct. 48 (3–4) (2011) 506–516. [9] B.L. Wang, Y.W. Mai, X.H. Zhang, Thermal shock resistance of functionally graded materials, Acta Mater. 52 (17) (2004) 4961–4972. [10] R. Samadhiya, A. Mukherjee, S. Schmauder, Characterization of discretely graded materials using acoustic wave propagation, Comput. Mater. Sci. 37 (1–2) (2006) 20–28. [11] D.W. Rosen, Computer-aided design for additive manufacturing of cellular structures, Comput.-Aided Des. Appl. 4 (5) (2007) 585–594. [12] F. Tamburrino, S. Graziosi, M. Bordegoni, The design process of additive manufactured Meso-Scale Lattice Structures: a review, J. Comput. Inf. Sci. Eng. (2018). [13] P. Vogiatzis, S. Chen, X. Wang, et al., Topology optimization of multi-material negative Poisson’s ratio metamaterials using a reconciled level set method, Comput. Aided Des. 83 (2017) 15–32. [14] E. Andreassen, B.S. Lazarov, O. Sigmund, Design of manufacturable 3D extremal elastic microstructure, Mech. Mater. 69 (1) (2014) 1–10. [15] M.P. Bendsœ, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Comput. Methods Appl. Mech. Engrg. 71 (2) (1988) 197–224. [16] M. Zhou, G.I.N. Rozvany, The COC algorithm, Part II: topological, geometrical and generalized shape optimization, Comput. Methods Appl. Mech. Engrg. 89 (1–3) (1991) 309–336. [17] M.P. Bendsøe, O. Sigmund, Material interpolation schemes in topology optimization, Arch. Appl. Mech. 69 (9–10) (1999) 635–654. [18] Y.M. Xie, G.P. Steven, A simple evolutionary procedure for structural optimization, Comput. Struct. 49 (5) (1993) 885–896. [19] M.Y. Wang, X. Wang, D. Guo, A level set method for structural topology optimization, Comput. Methods Appl. Mech. Engrg. 192 (1–2) (2003) 227–246. [20] G. Allaire, F. Jouve, A.M. Toader, Structural optimization using sensitivity analysis and a level-set method, J. Comput. Phys. 194 (1) (2004) 363–393. [21] M. Alzahrani, S.K. Choi, D.W. Rosen, Design of truss-like cellular structures using relative density mapping method, Mater. Des. 85 (2015) 349–360. [22] F. Schury, M. Stingl, F. Wein, Efficient two-scale optimization of manufacturable graded structures, SIAM J. Sci. Comput. 34 (6) (2012) B711–B733. [23] H. Li, Z. Luo, N. Zhang, et al., Integrated design of cellular composites using a level-set topology optimization method, Comput. Methods Appl. Mech. Engrg. 309 (2016) 453–475. [24] S. Daynes, S. Feih, W.F. Lu, et al., Optimisation of functionally graded lattice structures using isostatic lines, Mater. Des. 127 (2017) 215–223. [25] H. Li, Z. Luo, L. Gao, et al., Topology optimization for concurrent design of structures with multi-patch microstructures by level sets, Comput. Methods Appl. Mech. Engrg. 331 (2018) 536–561. [26] D. Li, W. Liao, N. Dai, et al., Optimal design and modeling of gyroid-based functionally graded cellular structures for additive manufacturing, Comput. Aided Des. 104 (2018) 87–99. [27] E. Andreassen, C.S. Andreasen, How to determine composite material properties using numerical homogenization, Comput. Mater. Sci. 83 (2014) 488–495. [28] S. Zhou, Q. Li, Design of graded two-phase microstructures for tailored elasticity gradients, J. Mater. Sci. 43 (15) (2008) 5157. [29] A. Radman, X. Huang, Y.M. Xie, Topology optimization of functionally graded cellular materials, J. Mater. Sci. 48 (4) (2013) 1503–1510. [30] Y. Wang, F. Chen, M.Y. Wang, Concurrent design with connectable graded microstructures, Comput. Methods Appl. Mech. Engrg. 317 (2017) 84–101. [31] A.D. Cramer, V.J. Challis, A.P. Roberts, Microstructure interpolation for macroscopic design, Struct. Multidiscip. Optim. 53 (3) (2016) 489–500. [32] A. Panesar, M. Abdi, D. Hickman, et al., Strategies for functionally graded lattice structures derived using topology optimisation for additive manufacturing, Additive Manuf. 19 (2018) 81–94. [33] H. Zong, H. Zhang, Y. Wang, et al., On two-step design of microstructure with desired Poisson’s ratio for AM, Mater. Des. 159 (2018) 90–102. [34] P.D. Dunning, Design parameterization for topology optimization by intersection of an implicit function, Comput. Methods Appl. Mech. Engrg. 317 (2017) 993–1011. [35] Y. Wang, S. Arabnejad, M. Tanzer, et al., Hip implant design with three-dimensional porous architecture of optimized graded density, J. Mech. Des. 140 (11) (2018) 111406. [36] Q. Xia, T. Shi, L. Xia, Stable hole nucleation in level set based topology optimization by using the material removal scheme of BESO, Comput. Methods Appl. Mech. Engrg. 343 (2019) 438–452. [37] Q. Xia, M.Y. Wang, Simultaneous optimization of the material properties and the topology of functionally graded structures, Comput. Aided Des. 40 (6) (2008) 660–675. [38] H. Liu, Y. Wang, H. Zong, et al., Efficient structure topology optimization by using the multiscale finite element method, Struct. Multidiscip. Optim. 58 (4) (2018) 1411–1430. [39] P. Wei, Z. Li, X. Li, et al., An 88-line MATLAB code for the parameterized level set method based topology optimization using radial basis functions, Struct. Multidiscip. Optim. 58 (2) (2018) 831–849. [40] Y. Liu, Z. Li, P. Wei, et al., Parameterized level-set based topology optimization method considering symmetry and pattern repetition constraints, Comput. Methods Appl. Mech. Engrg. 340 (2018) 1079–1101.
ScienceDirect Comput. Methods Appl. Mech. Engrg. 354 (2019) 487–505 www.elsevier.com/locate/cma
VCUT level set method for topology optimization of functionally graded cellular structures Hongming Zonga , Hui Liua,b , Qingping Maa , Ye Tiana , Mingdong Zhouc , Michael Yu Wanga ,∗ a
Department of Mechanical and Aerospace Engineering, Hong Kong University of Science and Technology, Kowloon, Hong Kong Department of Engineering Mechanics, School of Civil Engineering, Wuhan University, Wuhan, Hubei, People’s Republic of China c Shanghai Key Laboratory of Digital Manufacture for Thin-walled Structures, School of Mechanical Engineering, Shanghai Jiao Tong University, Shanghai, People’s Republic of China b
Received 15 November 2018; received in revised form 3 May 2019; accepted 20 May 2019 Available online 31 May 2019
Abstract This paper presents a variable cutting (VCUT) level set based topology optimization method to design functionally graded cellular structures (FGCS). A variable and continuous cutting function by interpolating with a set of height variables is proposed to generate functionally graded cellular structures, which offers a novel tool to optimize the macroscopic graded pattern. Due to the continuity of the cutting function, perfect geometric connections between adjacent cells are guaranteed without imposing extra constraints in the optimization. In addition, a solid covering skin can be easily attached to the FGCS by means of the variable cutting function. Three FGCS design problems, including graded density control, compliance minimization and layered cellular structure design, are presented to demonstrate the effectiveness and applicability of the proposed method. The optimized 2D and 3D macroscopic FGCS designs exhibit well-performing structural layout with fully connected micro-scale geometries. c 2019 Elsevier B.V. All rights reserved. ⃝ Keywords: Topology optimization; Functionally graded cellular structures; Level set method; Microstructure connectivity
1. Introduction The cellular structures composed of downscale microstructures are known for superior performance with relatively low densities [1]. By tailoring the configurations in the micro scale, excellent multifunctional responses can be achieved, such as stress resistance, heat exchange and shock absorption [2]. Besides, recent advances in additive manufacturing (AM) provide additional design freedom in fabricating these structures with complex and multiscale geometric features [3,4], which also stimulates the development of cellular structures. In particular, the functionally graded cellular structures [5] with spatially varying properties have gained increased popularity in engineering over the past decades. Compared with the periodically repeated cellular structures, progressively changing the geometric feature shapes or sizes of the embedded microstructures [6] allows FGCS ∗
Corresponding author. E-mail address: [email protected] (M.Y. Wang).
https://doi.org/10.1016/j.cma.2019.05.029 c 2019 Elsevier B.V. All rights reserved. 0045-7825/⃝
488
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
to flexibly adjust and customize the layout according to the loading and boundary conditions in practice. For example, the gradual variation in composition can reduce the stress concentration near the interface between different phases [7] and thus can be used to relieve the “stress shielding” phenomenon for orthopedic implants. A decreased density distribution of honeycomb cellular structure in the crushing direction can also be adopted to enhance the energy absorption at early crushing stages [8]. Other applications can also be found in for instance thermal insulation [9] and acoustical wave attenuation [10]. As a result, FGCS have attracted considerable attentions from the aerospace, automobile and medical industries. However, the FGCS common in biological systems ranging from the bird beaks, bones to plant stems cannot be easily designed due to the difficulties in modeling the configurations that typically span over several length scales and in determining a suitable graded pattern at the meantime. A widely used method for realizing this purpose is to first characterize the FGCS via explicit parameters and then run parametric optimization to obtain a graded distribution [11,12]. However, it is not easy to use an explicit parametric technique to represent microstructures involving general types of internal features and geometric variations. Therefore, many newly designed artificial materials [13,14] cannot be adopted by the FGCS, which limits the structural performance for multiple functionalities. Alternatively, topology optimization is a suitable design tool for FGCS because of its capability in generating innovative structural and material layouts with optimized performance. Interestingly, the very early homogenization method [15] is developed based on a similar concept of FGCS, where a hole is used in each cell and parametric optimization is conducted to determine the sizes of the holes. After three decades of development, various methods, e.g. the solid isotropic material with penalization method (SIMP) [16,17], the evolutionary structural optimization [18] and the level set method [19,20] have been proposed. Particularly, Alzahrani et al. [21] proposed a mapping technology to design FGCS, where the intermediate densities in SIMP are mapped to lattice unit cells with variable volume fractions. Schury et al. [22] used free material optimization method to produce manufacturable graded structures. Li et al. [23] developed a material-structure integrated topology optimization scheme for producing layer-wise FGCS. Daynes et al. [24] presented a novel biomimetic methodology for designing spatially graded lattice structures. Recently, Li et al. [25] proposed a concurrent topology optimization method that can generate cellular structures with multiple patches of microstructures rationally distributed in the structure, and Li et al. [26] presented a gyroid-based FGCS design method for optimizing lightweight lattice structures with uniformly varying densities in the design space. When designing FGCS, it is crucial to guarantee the connectivity between graded cells. If homogenization method [27] is applied directly without any geometric constraints, the scale separation assumption will result in unconnected regions, which cannot bear any load and will cause difficulties in additive manufacturing. Regarding this problem, Zhou and Li [28] made an in-depth discussion on this issue and developed three methods to guarantee the connectivity. Radman et al. [29] applied the filtering technique in every single microstructure design and preserved the connectivity in a sequential way. Though these implementations can relieve the disconnection problem to some extent, they can be hardly extended for complex graded patterns determined by the macroscopic loading and boundary conditions. Wang et al. [30] and Carmer et al. [31] used metamorphosis technology and interpolation scheme to generate connectable graded microstructures for macroscopic design. However, mismatch is observed in the obtained FGCS, which may cause stress concentrations and degrade the final structural performance. This paper develops a novel variable cutting level set based topology optimization method for FGCS design. A basic periodic level set function in which the microstructures are embedded is first constructed and distributed in a replicated way. Then a variable cutting function is introduced for generating the cellular structure with various microstructures and optimizing the graded pattern according to the requirement. Due to the periodicity of the basic level set function and the connectivity of the variable cutting function, the obtained FGCS are guaranteed to have at least C0 continuity. The proposed approach is demonstrated by several 2D and 3D FGCS design examples of graded density control, compliance minimization and layered structure design problems, respectively. 2. Modeling on VCUT level set method As shown in Fig. 1(a), the considered design domain with length L and height W for FGCS design is discretized into M unit cells with m nodes in the macro scale. For each unit cell c, it is further discretized in the micro scale as shown in Fig. 1(b). Based on such a discretization, the VCUT level set model is proposed to represent the graded cellular structures with two components: a basic level set function for the geometric details in the micro scale and a variable cutting function for the adjustable graded pattern in the macro scale.
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
489
Fig. 1. Design domain discretization for FGCS: (a) macro scale and (b) micro scale.
Fig. 2. General level set method for representing dynamic structural interfaces. (a) The level set function φ(x, t) and zero cutting plane, and (b) structural boundaries represented by the level set function.
2.1. Basic level set function In the general level set method, the dynamic structural interfaces are implicitly expressed by the zero level set of a higher-dimensional level set function φ(x, t) [19,20], as illustrated in Fig. 2. The geometry model can be represented by: ⎧ ⎨φ(x, t) > 0 ∀x ∈ Ω \∂Ω φ(x, t) = 0 ∀x ∈ ∂Ω (1) ⎩ φ(x, t) < 0 ∀x ∈ D\Ω , where x ∈ D is any point in the fixed design domain D and ∂Ω are the boundaries of the solid domain Ω . Considering the unit cell in a cellular structure, the level set method offers a flexible and convenient way for representing complex internal micro features. Theoretically, any type of microstructures can be digitally built based on the method and some 3D unit cells are illustrated in Fig. 3. One can use the signed distance function (Fig. 3(a– d)) [32,33] or the minimal surface function (Fig. 3(e)) [32] as the basic level set function φp (x) for a unit cell. In general, the signed distance function within the unit cell Dc can be calculated by: ⎧ ⎨+ min ∥x − x b ∥ ∀x ∈ Ωp x b ∈∂Ωp φp (x) = (2) ⎩− min ∥x − x b ∥ ∀x ∈ Dc \Ωp , x b ∈∂Ωp
where ∂Ωp denotes the boundaries of the microstructure. In addition, the unit cell with corresponding level set function φp (x) can also be obtained by using the inverse homogenization method within the level set framework. Some of the optimized unit cells with desired properties like negative Poisson’s ratio are given in Fig. 3(c, d) [33]. With a set of basic level set functions, extra operations such as the basic Boolean operations including intersection, union and difference between any two cells can be easily realized by a series of min/max operations
490
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
Fig. 3. Unit cells represented by level set method [32,33], including (a) octet truss, (b) carbon truss, (c, d) optimized microstructures with negative Poisson’s ratio, (e) Schwarz P surface unit cell and (f) artificial unit cell by union of a Schwarz P surface unit cell with a body-centered cubic unit cell.
on the level set functions. Fig. 3(f) is a unit cell generated by union of a Schwarz P surface unit cell with a body-centered cubic unit cell. These operations provide great design freedom for engineers to customize the unit cells according to practical design requirements. When the basic level set function φp (x) is determined, it could be mapped into macro elements in a replicated way, which will form a global level set function φg (x). Note that the φp (x) used in the VCUT level set method should possess geometric periodicity, which means that φp (x) on the opposite boundaries of the unit cell should be identical. This geometric periodicity may restrict the design space of its properties to some extent. However, it is still flexible enough for the unit cell to achieve various functionalities [23,29,30], including some extreme [14] or extraordinary properties [13]. 2.2. Variable cutting function As illustrated in Fig. 2, it is easy to notice that the boundaries are determined simultaneously by the level set function and the cutting function. Typically, the zero cutting plane is selected and fixed during the optimization. By evolving the level set function, the boundaries keep moving till an optimized structure is obtained. Alternatively, with a predefined level set function, optimization can also be performed by changing the cutting function globally [34]. A similar concept was used for graded cellular structure design in [30], where a series of unit cells were generated by adjusting the heights of the cutting planes as plotted in Fig. 4. The graded cellular structures were optimized by determining the heights in each macro element. However, as indicated in their work, mismatch on the common boundaries of neighboring cells cannot be avoided based on this modeling scheme. A simple illustration of this problem is given in Fig. 5, where the cellular structure is made up of 2 by 2 macro elements. Because the cutting planes are element-wisely distributed, there will be an offset between two neighboring cutting planes if the heights are different. As a result, mismatches may occur at the boundaries of adjacent cells. In order to achieve a natural geometric connection without the above mismatch issue, a variable and continuous cutting function is introduced in this paper. As shown in Fig. 1, a set of height variables are defined on the nodes of the entire macro mesh like the nodal density in [35], and they are represented by a time-dependent global vector: [ ]T H(t) = h 1 (t) h 2 (t) h 3 (t) . . . h m (t) , (3) where m denotes the total number of the macro mesh nodes. The height vector of a macro mesh element, for example c, is a subset of the global height vector: [ ]T hc (t) = h c1 (t) h c2 (t) h c3 (t) h c4 (t) = Sc H(t), (4) with Sc denoting its selection matrix in the common notation. The cutting function within the unit cell can be formulated by interpolating the hc (t) using the common shape functions. For convenience, a shape function vector N c (x) is defined as: [ ]T (5) N c (x) = N1c (x) N2c (x) N3c (x) N4c (x) , and the cutting function controlled by the height variables is given as: Cc (x, t) = N Tc (x)hc (t) = N Tc (x)Sc H(t).
(6)
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
491
Fig. 4. Signed distance function cut by planes at different heights. (a) The basic level set function generated on ∂Ωp in (b), and (b–e) generated microstructures.
Fig. 5. Illustration on the generation of mismatch in [30]. (a) The global level set function φg (x) and element-wisely distributed cutting planes, and (b) the generated cellular structure and mismatches at the boundaries of adjacent cells.
With the variable cutting function, the level set function for the geometry model of unit cell c becomes: ⎧ ⎨φc (x, t) = φp (x) − Cc (x, t) > 0 ∀x ∈ Ωc \∂Ωc φc (x, t) = φp (x) − Cc (x, t) = 0 ∀x ∈ ∂Ωc ⎩ φc (x, t) = φp (x) − Cc (x, t) < 0 ∀x ∈ Dc \Ωc .
(7)
It can be found that the geometric features in the macro elements are determined by the cutting function, or essentially by the cutting heights H(t). The cutting effect for a single cell is illustrated in Fig. 6. For any macro element, if the corresponding hc (t) is high, the element becomes a void element; if it is low, the element becomes a solid element; if it lays in between, the microstructure will exhibit geometric variations from its basic geometry, as shown in Fig. 6(b). In addition, since the heights do not have to be the same, the obtained unit cell can exhibit directional stiffness according to the external loading and boundary conditions and thus possesses a larger design space compared with the unit cells in Fig. 4. It should be noted that, in the current definition of the variable cutting function by Eq. (6), Cc (x, t) is mainly used to control the cell gradient in the macro scale and the potential topological variation of a cell in the micro scale
492
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
Fig. 6. Illustration on the variable cutting process. (a) The basic level set function φp (x) and the variable cutting function Cc (x, t) controlled by heights defined on four nodes, (b) the generated microstructure.
is determined by φp (x). However, if a more flexible function form is chosen, then the cutting function may be able to change the topology in the micro scale with the topological derivatives [36] provided. In this way, the structural topology in the micro scale becomes more flexible, which can help to enlarge the design space and improve the structural performance further. By assembling the cutting functions of unit cells together, a global cutting function Cg (x, t) with at least C0 continuity is formed. Combined with previous global level set function φg (x), the cellular structure can be modeled by: ⎧ ⎨φ(x, t) = φg (x) − Cg (x, t) > 0 ∀x ∈ Ω \∂Ω φ(x, t) = φg (x) − Cg (x, t) = 0 ∀x ∈ ∂Ω (8) ⎩ φ(x, t) = φg (x) − Cg (x, t) < 0 ∀x ∈ D\Ω , where the global boundaries ∂Ω : = {x|x ∈ ∂Ωc , c = 1, 2, . . . , M} are the assembly of cellular boundaries ∂Ωc . Since Cg (x, t) is variable and can be adjusted by changing H(t), cellular structures can yield variations in geometric features and thus graded patterns can be expected. Meanwhile, since neighboring cells share the same pair of height variables on common boundaries, the global cutting function Cg (x, t) has at least C0 continuity across the cells. Considering that the global level set function φg (x) also has C0 continuity because of the periodicity of φp (x), the resulted level set function φ(x, t) is C0 continuous. Therefore, the VCUT method can naturally guarantee the full geometric continuity for microstructures in macro elements without imposing any extra constraints. As an example, a 2 by 2 cellular structure is illustrated in Fig. 7, where the heights are set to be different and thus a graded distribution is observed. No geometric mismatches are observed at the connections between any two neighboring cells. Further, if higher-order interpolating functions are applied, e.g., using p-type shape functions, the cutting function Cg (x, t) will exhibit a geometric continuity higher than C0 . As a result, the cellular structure represented by the level set function can have higher-order continuity, e.g. C1 continuity, at the geometric boundaries between the neighboring cells. 3. Topology optimization using the VCUT level set method With the proposed VCUT model, the topology optimization problem is formulated as: min. J (u, Ω ) Ω
s.t.
⎧ ⎪ ⎪ ⎪ ⎨
a(u, w, ∫ Ω ) = l(w, Ω ) g(Ω ) =
∀w ∈ U (9)
dx − VD f v ≤ 0 Ω
⎪ ⎪ ⎪ ⎩
u = u0 Eε(u) · n = τ
in Γu in Γτ ,
where J(u, Ω ) is the objective functional, which can be defined in the following form without loss of generality: ∫ J (u, Ω ) = F(u)dx, (10) Ω
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
493
Fig. 7. Illustration on generation of graded cellular structure based on the VCUT level set method. (a) The global level set function φg (x) and the global cutting function Cg (x, t), (b) graded and connected cellular structure.
where g(Ω ) is the volume constraint, VD is the volume of the design domain and fv is the allowed maximum volume fraction; U is the kinematically admissible displacement space; a(u, w, Ω ) and l(w, Ω ) are the terms of the state equation of the linear elastic system: ∫ εT (u)Eε(w)dx, (11) a(u, w, Ω ) = Ω ∫ ∫ l(w, Ω ) = f T wdx + τ T wdx. (12) Ω
Γτ
The derivative of the objective function J(u, Ω ) is related to the variation of the solid domain Ω . In the VCUT method, since the solid domain is controlled by the height variable H(t), the relationship between the derivative of J(u, Ω ) and the variation of H(t) is derived in this section, which is used to optimize the graded pattern of the cellular structure. 3.1. Sensitivity analysis The Lagrange function of the optimization problem (9) is formulated as: L = J + a − l + λg,
(13)
where λ is the Lagrange multiplier. Following the two lemmas listed in [37], the shape derivative of Eq. (13) with respect to the pseudo time t is calculated by: L˙ = J˙ + a˙ − l˙ + λg, ˙
(14)
where: ∫ F ′ (u)dx + F(u)Vn dx, Ω ∂Ω ∫ ∫ ∫ a(u, ˙ w, Ω ) = εT (u′ )Eε(w)dx + εT (u)Eε(w′ )dx + εT (u)Eε(w)Vn dx, Ω Ω ∂Ω ∫ ∫ ∫ ∫ [ ] ( T )T T ′ T T ′ ˙l(w, Ω ) = f w dx + f wVn dx + τ w dx + ∇τ w n + τ T wk Vn dx, ∂Ω Γτ Γτ ∫ Ω g(Ω ˙ )= Vn dx, J˙(u, Ω ) =
∫
∂Ω
(15) (16) (17) (18)
494
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
where Vn is the normal velocity of the boundaries ∂Ω ; k is twice the mean curvature of the boundary defined as k = div n. Since the body force is assumed zero on the boundary ∂Ω and Γτ is fixed, the derivative of the virtual work can be simplified as: ∫ ∫ ˙ l(w, Ω) = f T w ′ dx + τ T w ′ dx. (19) Ω
Γτ
By Substituting all these terms into Eq. (14), the derivative of the Lagrange function can be reformulated as: {∫ } ∫ ∫ T ′ T ′ T ′ ˙L = ε (u)Eε(w )dx − f w dx − τ w dx Ω Ω {∫ } Γτ ∫ (20) + F ′ (u)dx + εT (u′ )Eε(w)dx Ω∫ {∫Ω } ∫ + F(u)Vn dx + εT (u)Eε(w)Vn dx + λ Vn dx . ∂Ω
∂Ω
∂Ω
Since a(u, w, Ω ) = l(w, Ω ) ∀w ∈ U and w ′ ∈ U, the first part of the derivative is canceled by: ∫ ∫ ∫ T ′ T ′ τ T w′ dx, f w dx + ε (u)Eε(w )dx = By solving the following adjoint equation: ∫ ∫ F ′ (u)dx = − εT (u′ )Eε(w)dx, Ω
(21)
Γτ
Ω
Ω
∀u′ ∈ U,
(22)
Ω
the derivative of the Lagrange equation can be obtained as: ∫ ∫ [ ] T ˙L = F(u) + ε (u)Eε(w) + λ Vn dx = G(u, w)Vn dx. ∂Ω
Moreover, by discretizing the global boundaries ∂Ω into cellular boundaries ∂Ωc , the derivative becomes: ∫ M ∫ ∑ c ˙L = G(u, w)Vn dx = G(u, w)Vnc dx, ∑ M c=1 ∂Ωc
(23)
∂Ω
c=1
(24)
∂Ωc
where Vnc is the normal velocity that drives ∂Ωc to move in the cell. 3.2. Update scheme for height variables When evolving the structural frontiers ∂Ωc in unit cells, the following Hamilton–Jacobi equation must be satisfied: dH(t) ∂φc (x, t) dφc (x, t) = − ∥∇φc (x, t)∥ Vnc = −N Tc (x)Sc − ∥∇φc (x, t)∥ Vnc = 0, ∀x ∈ ∂Ωc . (25) dt ∂t dt The relationship between the normal velocity Vnc and the design variables H(t) can be described as: 1 dH(t) N T (x)Sc , ∀x ∈ ∂Ωc . ∥∇φc (x, t)∥ c dt By substituting the normal velocity of Eq. (26) into Eq. (24), the derivative becomes: M ∫ ∑ 1 dH(t) G(u, w) N Tc (x)Sc L˙ = − dx ∥∇φ (x, t)∥ dt c c=1 ∂Ωc ] ∫ M [ dH(t) T ∑ T G(u, w) =− Sc N c (x)dx . dt ∂Ωc ∥∇φc (x, t)∥ c=1 Vnc = −
(26)
(27)
To guarantee the steepest descent searching direction, the derivative of the design variable H(t) is defined as: ] ∫ M [ dH(t) ∑ T G(u, w) = Sc N c (x)dx . (28) dt ∂Ωc ∥∇φc (x, t)∥ c=1
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
495
Fig. 8. Workflow of the VCUT level set method.
Fig. 9. Basic level set functions used in this paper, including (a) φp1 (x) and its zero level set (b), (c) φp2 (x) and its zero level set (d), and (e) zero level set of φp3 (x).
and the height variables are updated by: dH(ti ) ∆t, (29) dt where ti is the ith time step and ∆t is the time step size. The workflow of the proposed VCUT level set method is presented in Fig. 8. H(ti+1 ) = H(ti ) +
4. Numerical implementation In this paper, three basic level set functions are utilized for demonstration, including two of φp1 (x) in Fig. 9(a) and φp2 (x) in Fig. 9(c) for generating 2D unit cells and one of φp3 (x) in Fig. 9(e) for generating 3D unit cells. The φp2 (x) is constructed by the Schwarz P minimal surface formulation [32] as: φp2 (x) = cos(2π x) + cos(2π y),
∀x = (x, y) ∈ Dc .
(30)
The φp1 (x) and φp3 (x) are signed distance functions built based on Eq. (2), where truss central axes are chosen as ∂Ωp . For example, two diagonals of the square in Fig. 4(b) are taken as the boundaries ∂Ωp for constructing φp1 (x). Similarly, four diagonals of the cube are used for building φp3 (x). Then, All the basic level set functions are normalized to −1 and 1 to guarantee the generality of the numerical program. During the optimization, the state equation in (9) needs to be solved to obtain the displacement field u. Since the geometric features in the cellular structure typically span over several length orders, it is difficult to solve the state equation directly by the ordinary finite element method. In this study, an efficient multiscale finite element method (MsFEM) [38] is applied for mechanical analysis. In the MsFEM, the multiscale shape function is constructed to capture the micro-scale features in unit cells and thus the state equation can be solved on the macro scale. In addition, the construction process is parallelized to ensure a good computational efficiency.
496
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
In Eq. (28), the boundary integral is calculated to determine the derivatives. To avoid the time-consuming remeshing process when considering the exact boundaries during the analysis, the problem (9) is reformulated with the Heaviside function and Dirac function δ(φ) as in [19]. Under this formulation, the derivative of Eq. (28) becomes: ] ∫ M [ dH(t) ∑ T = Sc G(u, w)N c (x)δ(φc )dx (31) dt Dc c=1 where δ(φ) is chosen to be the same as that in [19]. In addition, the Lagrange multiplier λ is selected and updated using the scheme in [39]. To illustrate the generality of the proposed VCUT level set method, three design problems, including graded density distribution, compliance minimization and layered structure design, are studied. For different problems, the optimization setup (9) or the graded pattern may be different. Hence, the formulation of the derivative in Eq. (31) should be modified accordingly. Graded density control The first problem is to obtain FGCS with a given graded density distribution. It is formulated as an unconstrained optimization problem and the objective function can be defined as: M
J (Ω ) =
1∑ (Vc − f c∗ VDc )2 , 2 c=1
(32)
where Vc is the volume of the microstructure in unit cell c, f c∗ indicates the predefined target volume ratio of cell c and VDc is the volume of the unit cell domain. The height derivative under this objective function can be easily derived as: ] ∫ M [ ) dH(t) ∑ ( = Vc − f c∗ VDc STc N c (x)δ(φc )dx . (33) dt Dc c=1 Compliance minimization problem For the compliance minimization problem, the objective function is: ∫ 1 εT (u)Eε(u)dx. J (u, Ω ) = 2 Ω The corresponding adjoint equation is: ∫ ∫ εT (u′ )Eε(u)dx, ∀u′ ∈ U, εT (u′ )Eε(w)dx = −
(34)
(35)
Ω
Ω
and thus the adjoint displacement is simply w = -u. The derivative can then be calculated by: ] } ∫ [ M { dH(t) ∑ T 1 − εT (u)Eε(u) + λ N c (x)δ(φc )dx . = Sc dt 2 Dc c=1
(36)
Layered cellular structure design Layered cellular structure is one typical type of FGCS, where the geometric features in microstructures vary in a multilayered manner. Here, the layered structures are designed for a minimized structural compliance. Different from the previous example, extra operations must be applied in order to obtain layered distribution. Similar to the pattern repetition constraints [40], the obtained derivatives and height variables in the same layer are averaged in each iteration step. Suppose there are n nodes in each layer, the averaged node height in the l-th layer can be calculated by: hl =
1 n
ln ∑
hi .
(37)
i=(l−1)n+1
This process is illustrated in Fig. 10. Since the structure is still modeled based on the VCUT level set method, this operation will not influence the connection of the generated cellular structure.
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
497
Fig. 10. Process for generating layered cellular structure within the VCUT method. (a) The height variables before the operation, and (b) the height variables after the operation.
Fig. 11. (a) The predefined density distribution, and (b) FGCS obtained by the VCUT method.
5. Numerical examples In this section, several 2D and 3D numerical examples on the three design problems indicated in Section 4 are used to demonstrate the proposed VCUT level set method for FGCS design. For all examples, the Young’s moduli for the solid material is Es = 1 and for the artificial weak material is Ev = 10−6 . The latter is used to avoid numerical singularity when solving the state equation. The Poisson’s ratios for both materials are 0.3. Unless otherwise specified, the optimization terminates when the relative change of 10 successive objective functions is below 0.1%, or the maximum 300 iteration steps is reached. For 2D problems, the same design domain D with W = 200 and L = 600 is used, which is discretized into 60 × 20 unit cells and each unit cell is further discretized by 10 × 10 4-node elements. For 3D problems, the same design domain D with the dimension of 128 × 32 × 64 is employed, which is discretized into 16 × 4 × 8 unit cells and each unit cell is further discretized into 8 × 8 × 8 8-node elements. 5.1. Graded density control In this problem, since the objective function is evaluated based on the absolute differences between the obtained design and the target design, a different terminate condition is used that the objective function is less than 0.01. The target density distribution ranging from 0.2 to 0.8 is given in Fig. 11(a) and the φp1 (x) is used for illustration. After 77 iteration steps, the objective function reduces to 9.56 × 10−3 and the produced design is shown in Fig. 11(b). Note that, although an iterative optimization process is carried out here, the computational cost is very low because no mechanical analysis is involved. The resulted design exhibits an obvious graded variation in both geometric features and unit cell densities. In addition, the magnified plot indicates that each unit cell is connected well as expected. With this simple example, the VCUT level set method is proved to be capable of controlling the density distribution of cellular structure, which provides designers with a convenient and efficient way to adjust the graded pattern in practical usage.
498
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
Fig. 12. Loading and boundary conditions for 2D cantilever beam design.
Fig. 13. Iteration history of 2D cantilever beam design based on φp1 (x) using the VCUT level set method, and (a–d) the intermediated designs during the optimization.
5.2. Compliance minimization Case 1 2D cantilever beam A benchmark cantilever beam design problem is studied in this example, where the loading and boundary conditions are illustrated in Fig. 12. The beam is fixed at left-hand side and loaded vertically at the right-bottom point with F = 1. The volume constraint is set to fv = 0.5. Fig. 13 shows the iteration history of the design process based on φp1 (x). After 157 iteration steps, the optimization converges to the solution as shown in Fig. 14(a) with the objective function J = 1.86 × 102 . With the same problem setup, the optimization with the other basic level set function φp2 (x) converges after 109 iteration steps with the objective function J = 1.95 × 102 . The obtained FGCS is plotted in Fig. 14(b). Controlled by the variable cutting functions, for instance in Fig. 15 for the case of φp1 (x), two solutions present a similar and rational graded pattern to bear the vertical load at the corner. For areas that contribute little to bearing the load, the cutting function moves beyond the maximum value of φp (x) and generates the void regions. For the load path areas, the cutting function moves below the minimum value of φp (x) and generates the solid region. For areas in between, microstructures with low densities are generated to support the two main solid frames, similar to the solutions obtained by the classical SIMP method with penalty factor of p = 1. In the classical SIMP method, the effective properties of intermediate densities with penalty factor of p = 1 are better than those with penalty factor of p = 3. As a result, the optimized results with p = 1 typically have a great percentage of gray regions. In our case, the cells with porous microstructures play a role of the elements with intermediate densities. As indicated in [30], this kind of graded cellular structures can provide a better mechanical performance. In addition, the geometric features of the microstructures are perfectly connected when the cells change gradually, as shown in the magnified plots of Fig. 14. Note that, in the transition areas where the cells fade to void, the cutting function will intersect with the maximum value of the basic level set function. As a result, small burrs appear on the boundaries. To achieve smooth boundaries, a solid covering skin can be attached to the FGCS by means of the global cutting function, which can also act as a post-processing to facilitate the real manufacturing [4].
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
499
Fig. 14. Obtained FGCS based on (a) φp1 (x) and (b) φp2 (x) using the VCUT level set method, and (c) and (d) the FGCS with the covering skins based on (a) and (b).
Fig. 15. Variable cutting function for the case of φp1 (x).
Since the cutting function is a continuous function and the burrs only occur when Cg (x) cuts across the maximum value of φp (x), the boundaries of the FGCS in the macroscopic scale ∂ΩFGCS can be determined by: Cg (x) − m p = 0,
∀x ∈ ∂ΩFGCS ,
(38)
where mp denotes the maximum value of φp (x). Based on the ∂ΩFGCS , a new signed distance function φFGCS can be built according to Eq. (2). Considering the boundaries of the design domain, another signed distance function φD can also be built. Finally, the function for generating the whole covering skin can be calculated by a min operation as: φskin (x) = min(φFGCS , φD ),
(39)
and the new level set function φnew (x) to represent FGCS with covering skins can be represented by: φnew (x) = max(min(φskin , −(φskin − θ)), φ),
(40)
where θ is the parameter to control the thickness of the skin. By choosing θ = 3 for illustration, φnew (x) based on the solutions in Fig. 14(a, b) are calculated and then substituted into Eq. (1) to determine the solid areas. Two final FGCS with covering skins are illustrated in Fig. 14(c, d). Case 2 2D simply supported beam Fig. 16 defines the loading and boundary conditions for the 2D simply supported beam design. The left-bottom point of the beam is fixed and the vertical degree of freedom at the right-bottom point is constrained. A vertical load
500
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
Fig. 16. Loading and boundary conditions for 2D simply supported beam design.
Fig. 17. Obtained FGCS based on (a) φp1 (x) and (b) φp2 (x) using the VCUT level set method, and (c) and (d) the FGCS with skins based on (a) and (b).
with F = 1 is applied at the center of the upper edge. The volume constraint is set to fv = 0.5. The optimization based on φp1 (x) converges to J = 26.4 after 114 iteration steps and that based on φp2 (x) ends up with J = 27.9 after 120 iteration steps. The obtained FGCS as well as those with covering skins are presented in Fig. 17. The two designs exhibit different topologies in the central area, where the cutting functions of the corresponding designs are compared in Fig. 18. Such a difference can be mainly attributed to the different intrinsic geometric features of the basic level set functions. Simple models consisting of two unit cells as shown in Fig. 19 are used for explanation. Since the height variables are defined on the nodes, when the height on node A is set to 1, the two cells lose connection around this node immediately. As a comparison, when the height on node B equals 1, the two cells are still connected at the common boundary and thus it is not easy to generate holes around the center line. In addition, two unit cells located around the load paths are extracted as shown in Fig. 17(a). Since the FGCS are symmetric with respect to the vertical center line, configurations of the two cells are mirrored with each other. Even with the same density value, they actually exhibit directional stiffness to follow the macroscopic loadings. Compared with the methods that characterize the unit cells with only density variables, the VCUT level set method has a larger design space and potential to further improve the mechanical performance of FGCS. Case 3 3D cantilever beam A 3D cantilever beam as shown in Fig. 20 is optimized based on the basic level set function of φp3 (x). The beam is fixed at the left face and loaded at the center of right-bottom edge with F = 100. The volume fraction is set to fv = 0.4. After 182 iteration steps, the optimization converges to J = 3.16 × 104 and the evolving history is plotted in Fig. 21. The final design is presented in Fig. 22 and that with a covering skin is shown in Fig. 23. The 3D microstructures are distributed in a graded way for bearing the load and those which do not contribute are wiped out. Since the thickness is relatively small, the result shows a similar configuration with that in Fig. 14. Nevertheless, since the
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
501
Fig. 18. Variable cutting function for the case of (a) φp1 (x) and (b) φp2 (x).
Fig. 19. Simple cellular structures consisting of two unit cells based on (a) φp1 (x) and (b) φp2 (x).
Fig. 20. Loading and boundary conditions for 3D cantilever beam design.
concentrated load is applied at the center, the variations in the thickness direction can still be observed. In addition, the connections between graded cells are verified through the magnified plots, where no mismatch is observed. 5.3. Layered cellular structure design Case 1 2D layered cellular structure Using the same loading and boundary conditions in Fig. 12, layered cellular structures are optimized based on the basic level set function φp1 (x). Fig. 24 gives the iteration history of the optimization, which converges after 151 steps with the objective function J = 2.57 × 102 . The final solution is presented in Fig. 25(a), where the upper and bottom layers end up with solid microstructures and the middle regions are filled with gradually varied microstructures. This graded distribution pattern is determined by the cutting function as illustrated in Fig. 25(b),
502
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
Fig. 21. Iteration history of 3D cantilever beam design.
Fig. 22. Optimized 3D cantilever beam using the VCUT level set method.
Fig. 23. Optimized 3D cantilever beam with covering skin.
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
503
Fig. 24. Iteration history of 2D layered cellular structure design.
Fig. 25. (a) Obtained 2D layered cellular structure based on φp1 (x) and (b) corresponding cutting function.
where the values at the two ends are below −1, which is the minimal value of φp1 (x), and then gradually increase to a value close to 1 in the middle. Case 2 3D layered cellular structure The same 3D cantilever beam problem in Fig. 20 is solved within the layered design scheme. The optimization terminates after 277 steps with J = 4.71 × 104 , and its iteration history is illustrated in Fig. 26. The obtained solution is shown in Fig. 27, which is similar to that in Fig. 25(a). Two solid plates are generated and graded cells are distributed in between to support them. It should be noted that, the compliance in the layered cellular structures is higher than that of Section 5.2 Case 3, because the operation in Eq. (37) essentially restricts the design space. Even though, according to the magnified plot in Fig. 27, it is proved that the operation has no influence on the microstructure connection. Geometric features in graded cells connect well on the common boundaries since the layered cellular structure is still modeled by the VCUT method and thus the connection is guaranteed. 6. Conclusions In this paper, a variable cutting level set topology optimization method is developed for designing functionally graded cellular structures. Controlled by the heights defined on macro element nodes, the continuous cutting function can be optimized to generate the graded pattern in the macro scale for the cellular structures. In addition, due to the continuity of the cutting function and the periodicity of the basic level set function, the obtained FGCS are ensured to have at least C0 continuity. Based on the method, three design problems, including graded density control, compliance minimization and layered cellular structure design, are studied. Several 2D and 3D numerical examples are given, and rational graded patterns as well as perfect connections among cells are observed, which demonstrates the effectiveness of the proposed design method. In this study, three predefined basic level set functions are used for illustrating the ability of the VCUT method in determining optimized graded patterns. In the future, a more general concurrent design method can be developed by combining the VCUT method with the inverse homogenization method. Given specific loading and boundary conditions, the inverse homogenization method can be used to generate the optimized basic level set function in the micro scale, and the VCUT method can then be adopted simultaneously to obtain the optimized graded pattern in the macro scale. By such a concurrent design method, further performance improvement can be expected. In
504
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
Fig. 26. Iteration history of 3D layered cellular structure design.
Fig. 27. Obtained 3D layered cellular structure using the VCUT level set method.
addition, only the bilinear or trilinear shape function is adopted as the interpolating function to form the cutting function in this paper. In the future, a higher-order interpolating function can be used to build a more flexible cutting function and to achieve higher-order continuity. Acknowledgment The authors acknowledge the support from National Natural Science Foundation of China (Grant No. 51705311). References [1] L.J. Gibson, M.F. Ashby, Cellular Solids: Structure and Properties, Cambridge University Press, 1999. [2] H.N.G. Wadley, Multifunctional periodic cellular metals, Phil. Trans. R. Soc. A 364 (1838) (2006) 31–68. [3] M.G. Rashed, M. Ashraf, R.A.W. Mines, et al., Metallic microlattice materials: a current state of the art on manufacturing, mechanical properties and applications, Mater. Des. 95 (2016) 518–533. [4] Y. Wang, L. Zhang, S. Daynes, et al., Design of graded lattice structure with optimized mesostructures for additive manufacturing, Mater. Des. 142 (2018) 114–123. [5] D. Mahmoud, M.A. Elbestawi, Lattice structures and functionally graded materials applications in additive manufacturing of orthopedic implants: a review, J. Manuf. Mater. Process. 1 (2) (2017) 13. [6] M. Niino, S. Maeda, Recent development status of functionally gradient materials, ISIJ Int. 30 (9) (1990) 699–703.
H. Zong, H. Liu, Q. Ma et al. / Computer Methods in Applied Mechanics and Engineering 354 (2019) 487–505
505
[7] K. Shah, I. ul Haq, A. Khan, et al., Parametric study of development of inconel-steel functionally graded materials by laser direct metal deposition, Mater. Des. 54 (2014) 531–538. [8] A. Ajdari, H. Nayeb-Hashemi, A. Vaziri, Dynamic crushing and energy absorption of regular, irregular and functionally graded cellular structures, Int. J. Solids Struct. 48 (3–4) (2011) 506–516. [9] B.L. Wang, Y.W. Mai, X.H. Zhang, Thermal shock resistance of functionally graded materials, Acta Mater. 52 (17) (2004) 4961–4972. [10] R. Samadhiya, A. Mukherjee, S. Schmauder, Characterization of discretely graded materials using acoustic wave propagation, Comput. Mater. Sci. 37 (1–2) (2006) 20–28. [11] D.W. Rosen, Computer-aided design for additive manufacturing of cellular structures, Comput.-Aided Des. Appl. 4 (5) (2007) 585–594. [12] F. Tamburrino, S. Graziosi, M. Bordegoni, The design process of additive manufactured Meso-Scale Lattice Structures: a review, J. Comput. Inf. Sci. Eng. (2018). [13] P. Vogiatzis, S. Chen, X. Wang, et al., Topology optimization of multi-material negative Poisson’s ratio metamaterials using a reconciled level set method, Comput. Aided Des. 83 (2017) 15–32. [14] E. Andreassen, B.S. Lazarov, O. Sigmund, Design of manufacturable 3D extremal elastic microstructure, Mech. Mater. 69 (1) (2014) 1–10. [15] M.P. Bendsœ, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Comput. Methods Appl. Mech. Engrg. 71 (2) (1988) 197–224. [16] M. Zhou, G.I.N. Rozvany, The COC algorithm, Part II: topological, geometrical and generalized shape optimization, Comput. Methods Appl. Mech. Engrg. 89 (1–3) (1991) 309–336. [17] M.P. Bendsøe, O. Sigmund, Material interpolation schemes in topology optimization, Arch. Appl. Mech. 69 (9–10) (1999) 635–654. [18] Y.M. Xie, G.P. Steven, A simple evolutionary procedure for structural optimization, Comput. Struct. 49 (5) (1993) 885–896. [19] M.Y. Wang, X. Wang, D. Guo, A level set method for structural topology optimization, Comput. Methods Appl. Mech. Engrg. 192 (1–2) (2003) 227–246. [20] G. Allaire, F. Jouve, A.M. Toader, Structural optimization using sensitivity analysis and a level-set method, J. Comput. Phys. 194 (1) (2004) 363–393. [21] M. Alzahrani, S.K. Choi, D.W. Rosen, Design of truss-like cellular structures using relative density mapping method, Mater. Des. 85 (2015) 349–360. [22] F. Schury, M. Stingl, F. Wein, Efficient two-scale optimization of manufacturable graded structures, SIAM J. Sci. Comput. 34 (6) (2012) B711–B733. [23] H. Li, Z. Luo, N. Zhang, et al., Integrated design of cellular composites using a level-set topology optimization method, Comput. Methods Appl. Mech. Engrg. 309 (2016) 453–475. [24] S. Daynes, S. Feih, W.F. Lu, et al., Optimisation of functionally graded lattice structures using isostatic lines, Mater. Des. 127 (2017) 215–223. [25] H. Li, Z. Luo, L. Gao, et al., Topology optimization for concurrent design of structures with multi-patch microstructures by level sets, Comput. Methods Appl. Mech. Engrg. 331 (2018) 536–561. [26] D. Li, W. Liao, N. Dai, et al., Optimal design and modeling of gyroid-based functionally graded cellular structures for additive manufacturing, Comput. Aided Des. 104 (2018) 87–99. [27] E. Andreassen, C.S. Andreasen, How to determine composite material properties using numerical homogenization, Comput. Mater. Sci. 83 (2014) 488–495. [28] S. Zhou, Q. Li, Design of graded two-phase microstructures for tailored elasticity gradients, J. Mater. Sci. 43 (15) (2008) 5157. [29] A. Radman, X. Huang, Y.M. Xie, Topology optimization of functionally graded cellular materials, J. Mater. Sci. 48 (4) (2013) 1503–1510. [30] Y. Wang, F. Chen, M.Y. Wang, Concurrent design with connectable graded microstructures, Comput. Methods Appl. Mech. Engrg. 317 (2017) 84–101. [31] A.D. Cramer, V.J. Challis, A.P. Roberts, Microstructure interpolation for macroscopic design, Struct. Multidiscip. Optim. 53 (3) (2016) 489–500. [32] A. Panesar, M. Abdi, D. Hickman, et al., Strategies for functionally graded lattice structures derived using topology optimisation for additive manufacturing, Additive Manuf. 19 (2018) 81–94. [33] H. Zong, H. Zhang, Y. Wang, et al., On two-step design of microstructure with desired Poisson’s ratio for AM, Mater. Des. 159 (2018) 90–102. [34] P.D. Dunning, Design parameterization for topology optimization by intersection of an implicit function, Comput. Methods Appl. Mech. Engrg. 317 (2017) 993–1011. [35] Y. Wang, S. Arabnejad, M. Tanzer, et al., Hip implant design with three-dimensional porous architecture of optimized graded density, J. Mech. Des. 140 (11) (2018) 111406. [36] Q. Xia, T. Shi, L. Xia, Stable hole nucleation in level set based topology optimization by using the material removal scheme of BESO, Comput. Methods Appl. Mech. Engrg. 343 (2019) 438–452. [37] Q. Xia, M.Y. Wang, Simultaneous optimization of the material properties and the topology of functionally graded structures, Comput. Aided Des. 40 (6) (2008) 660–675. [38] H. Liu, Y. Wang, H. Zong, et al., Efficient structure topology optimization by using the multiscale finite element method, Struct. Multidiscip. Optim. 58 (4) (2018) 1411–1430. [39] P. Wei, Z. Li, X. Li, et al., An 88-line MATLAB code for the parameterized level set method based topology optimization using radial basis functions, Struct. Multidiscip. Optim. 58 (2) (2018) 831–849. [40] Y. Liu, Z. Li, P. Wei, et al., Parameterized level-set based topology optimization method considering symmetry and pattern repetition constraints, Comput. Methods Appl. Mech. Engrg. 340 (2018) 1079–1101.