Periodic cellular structures with excellent mechanical properties widely exist in nature. A generative design and optimization method for triply periodic level surface (TPLS)-based functionally graded cellular structures is developed in this work. In the proposed method, by controlling the density distribution, the designed TPLS-based cellular structures can achieve better structural or thermal performances without increasing its weight. The proposed technique can be divided into four steps. First, the modified 3D implicit functions of the triply periodic minimal surfaces are developed to design different types of cellular structures parametrically and generate spatially graded cellular structures. Second, the numerical homogenization method is employed to calculate the elastic tensor and the thermal conductivity tensor of the cellular structures with different densities. Third, the optimal relative density distribution of the object is computed by the scaling laws of the TPLS-based cellular structures added optimization algorithm. Finally, the relative density of the numerical results of structure optimization is mapped into the modified parametric 3D implicit functions, which generates an optimum lightweight cellular structure. The optimized results are validated subjected to different design specifications. The effectiveness and robustness of the obtained structures is analyzed through finite element analysis and experiments. The results show that the functional gradient cellular structure is much stiffer and has better heat conductivity than the uniform cellular structure.

## Introduction

Additive manufacturing (AM) allows designers to create and fabricate more complex shapes and topologies. However, current design for manufacturing techniques have been developed largely for conventional manufacturing processes [1]. Therefore, to take advantage of the design freedom enabled by AM processes, design for additive manufacturing (DfAM) methods need to be developed [1,2]. Representative DfAM methodologies include cellular structure design, topology optimization (TO), manufacturability analysis for AM, and other design approaches that can make use of the AM-enabled features [1]. Among them, an artificially designed cellular structure possesses many superior properties that make it a promising solution for various applications, such as improved stiffness [3], heat transfer [4,5], energy absorption [6], and so on. In addition, cellular structures exist widely in nature. For example, the trabecular structures of human bones, the skulls of birds, and the microstructures of bamboo all consist of gradient cellular structures [7]. These naturally evolved gradient structures distribute the material reasonably to the areas where performance needs to be optimized. These excellent performances include high specific stiffness, impact resistance, and high thermal conductivity. Inspired by gradient cellular structures in nature, we will focus on the design and optimization of functionally graded cellular structures (FGCS) in this work.

The traditional cellular structures are arranged by the regular cell structures or composed of random foam structures [3]. For instance, the designers choose the regular cellular structure in commercial software with a given volume fraction or generate the random cellular structure by the Voronoi algorithm. However, there are still some issues in designing cellular structures for AM. Usually, the lightweight design of the components use only uniform cellular structures for infilling, which cannot satisfy a reasonable stress distribution. Besides, although some optimization methods [8,9] control the cellular infill density through stress distribution based on finite element analysis (FEA), such optimizations cannot obtain the optimal solution. In contrast, the TO method can effectively optimize the distribution of materials in space. To further optimize the cellular infilling, the gradient cellular structures based on TO calculation design methods [10] have been rapidly developed. However, most methods simply match the stress size to the relative density of the cellular structures and the structural properties of the cellular itself are not considered when numerical calculating is done. In addition, these methods generally employ truss-based lattice structures, making it difficult to obtain smooth geometries that can be used for direct 3D printing. Therefore, the above problems need to be optimized and solved.

In this work, a generative design and optimization approach for the generation of functionally graded triply periodic level surface (TPLS)-based cellular structures is developed. It optimizes the mechanical stiffness or thermal conduction by redistributing the relative density of cellular structures in a design domain and under given boundary conditions. Notably, to generate graded cellular structures, a TPLS-based function representation (F-Rep) modeling method that can obtain geometrically continuous cellular structures is proposed. One of the distinctive benefits of adopting this modeling method is that the assigned parameters can control the shape and volume fraction distribution of a TPLS-based cellular structure easily. The reason why the TPLS-based cellular structures are chosen for this work is that they offer the following potential advantages over the strut-based lattice structures: (1) the open cell structure allows easy removal of powders or liquid resins; (2) the self-supporting helical structures that do not require any supporting material during the AM process [9,11]; (3) a high specific surface area and low relative density while maintaining a higher specific strength than regular diamond cellulars in certain circumstances; and (4) a functional gradient structure can be generated by the TPLS-based cellular structure easily by controlling the level-set parameter. Furthermore, the effective mechanical and thermal conductive properties of the TPLS-based cellular structures can be calculated via the computational homogenization method [12,13]. The obtained mechanical or thermal conductive properties are added to the topology optimization process to generate gradient cellular structures with effective physical meaning. Finally, for validation purposes, a variety of design specifications are given for optimization, and the effectiveness and robustness of the obtained structures are analyzed through FEA and further experiments.

## Related Works

In the field of DfAM, many researchers focused on the optimization of geometric shapes and topology structures with specific mechanical properties. In this section, a brief review related to the design and optimization of cellular structures for AM is made.

*TPLS-based cellular structures* derived from the triply periodic minimal surface [14] are of special interest to many scientists. A typical characteristic is that the mathematical expression represents an equation composed of trigonometric functions [14], which have no intersecting or folded surfaces. “Triply periodic” means that the structure can be patterned in the 3D space, and “level surface” means that it can adjust the size of the separated area to control the density of the cellular structures. The main advantage of TPLS-based cellular structures is that its opening cell structure has high specific surface area, high porosity, and low relative density while maintaining excellent mechanical properties [15,16]. Furthermore, TPLS has been discovered to exist in the nature stable, such as the gyroid structure in butterfly wing scales [17] and the Schwartz Diamond structure in diamonds [18]. Previously, the stiffness and strength [15,16] and manufacturability by AM [9,11] of TPLS-based cellular structures were investigated numerically and experimentally. Besides, several researchers have employed TPLS-based cellular structures in a range of applications to enhance functionalities through their topology properties relationship. For example, they have used the gyroid structure as a supporting structure for AM [11], TPLS-based cellular structures as a bone tissue engineering scaffold [19], TPLS-based cellular structures with varying transport of heat and electricity [20], and TPLS-based composite structures with high conductivities [21], and so on.

*Functionally graded cellular structure* is a type of functionally graded material (FGM). FGM is mainly designed by having a variable proportion of composite materials with random or regular microstructures. FGCS typically refers to the variable-density cellular or porous structures design according to the physical and mechanical properties [22–24]. Here, we concentrate on discussing FGCS design and optimization based on topology optimization. Brackett et al. [10] first proposed a solid isotropic material with penalization (SIMP)-based density mapping approach to optimize the intermediate densities to 2D lattice structures of varying volume fraction. Panesar et al. [25] used the density of topology optimization directly mapping into the density of the infilled lattice structures. Liu and Shapiro [26] proposed a method to generate a spatially varying structure of FGM with target properties from computed tomography image or topology optimization. Wu et al. [27] presented a self-supporting “rhombic cell” topology gradient infill design system for 3D printing object optimization. Tang et al. [28] proposed a bidirectional evolutionary structural optimization (BESO)-based method for truss-like lattice structure optimization. Alzahrani et al. [29] designed truss-like cellular structures using the relative density mapping (RDM) method. Daynes et al. [30] presented a combined methodology of topology optimization and size optimization. In addition, Cheng et al. [31] took the anisotropic property of cellular materials into account by the homogenization method. In addition, Wang et al. [32] used an efficient homogenization-based topology optimization method to generate a variable-density honeycombs cellular structure for natural frequency optimization. Cheng et al. [4,5] generated the graded lattice structures to optimize the heat conduction of the object based on level-set methods.

*Topology optimization* is a mathematical technique that optimizes the distribution of material in a specific design domain, with the goal of maximizing system performance for a given set of loads and constraints. The popular topology optimization approaches include the homogenization-based method [33], SIMP [34], BESO [35], the level-set method [36], the moving morphable component/moving morphable void (MMC/MMV) methods [37–39], and various algorithms derived from them. The MMC method is used in conjunction with the homogenization method to generate infilling gradient microstructure faster [40]. The BESO algorithm is a finite element-based topology optimization approach, where inefficient material is iteratively removed from a structure while efficient material is simultaneously added to the structure [35,41]. The basic principle of the SIMP algorithm is to calculate the evolution of the relative density distribution of the elements over the design space [34,42]. For a comprehensive review of SIMP and BESO topology optimization techniques, see the recent survey articles [41,42].

Our work is inspired by these works, yet it avoids the use of truss-based lattice structures to generate structural compositions and avoids the density calculated from original topology optimization directly. On the one hand, truss-like cellular structures are limited in the functionally graded modeling method, and some topologies are limited in manufacture by AM. In comparison, it is easy to generate a continuous gradient cellular structure using the TPLS-based implicit modeling method. On the other hand, most of the traditional topology optimization methods are calculated for isotropic materials. However, few cellular structures are completely isotropic. Most of them are anisotropic. In our work, the elastic tensor and the thermal conductivity tensor scaling laws of TPLS-based cellular structures are added to the optimization algorithm. This allows the optimized relative density distributions of the cellular structures distributed in a component to possess an accurate physical meaning.

## Overview

The proposed design methodology for FGCS optimization mainly consists of four steps, which is illustrated by an example of the torque transfer optimization problem in the flow chart in Fig. 1. In step 1, the modified 3D implicit functions of several triply periodic minimal surfaces are developed to parametrically design different types of cellular structures and generate spatially graded cellular structures (see Sec. 4). In step 2, the numerical homogenization method is used to calculate the elastic tensor and the thermal conductivity tensor of TPLS-based cellular structures with different densities. According to the results, the equation of density is established to demonstrate the equivalent elastic tensor and the thermal conductivity tensor of the specific cellular structures (see Sec. 5). In step 3, the penalization function in the original SIMP-based optimization method is replaced by the calculated homogenized model. The accurate volume fraction of each unit is obtained by figuring out the tensor scaling laws added to the topology optimization problem (see Sec. 6.1). In step 4, the relative density of the numerical results of structure optimization is mapped into the parametric TPLS-based implicit functions, which realize an optimum lightweight cellular structure (see Sec. 6.2). Finally, FEA simulation is conducted to verify the structural robustness of our optimized results in Sec. 7, before conclusions are drawn in Sec. 8.

## TPLS-Based Cellular Structures Design

### Mathematical Description.

*F*:

*R*

^{3}→

*R*at point (

*x*,

*y*,

*z*) ∈

*R*

^{3}, which satisfies the equation

*t*is the level parametric. The form of the surface is controlled by

*F*(

*x*,

*y*,

*z*), while the parameter

*t*determines the volume fraction of the region. A periodic surface

*F*(

*x*,

*y*,

*z*) can generally be defined as [43]

*ϕ*(

*hkl*) denotes the structure magnitude factor reflecting the symmetry of the structure; the set of allowed

*hkl*values express the reciprocal cellular vectors of a given lattice [44];

*α*

_{hkl}represents the phase shift, and

*L*is the edge length of periods. Therefore, specific periodic structures and phases can be constructed based on the implicit form.

*x*,

*y*, and

*z*coordinates of a point in three-dimensional Euclidean space. The close approximation to the G, D, P, and W level surfaces [44] are composed of trigonometric functions

*F*:

*R*

^{3}→

*R*at point (

*x*,

*y*,

*z*) ∈

*R*

^{3}, and they are defined as follows:

*t*is a variable and determines the volume fractions pertaining to the regions separated by the surface. The cellular structures are subjected to translation and scaling. Moreover, they are pinched off into unconnected blobs as the value of the parameter

*t*becomes sufficiently negative or positive. The relative density, or the volume fraction of the solid phase, is an important parameter for cellular materials. It is usually used as a variable to express other properties, such as elastic modulus and yield stress. It is essential to acquire the relationship between parameter

*t*and relative density

*ρ*

_{TPLS}for TPLS-based cellular structures generated from Eqs. (3)–(6), respectively. The data plotted in Fig. 3 reflect how the relative density of G, D, P, and W structures change with the parameter

*t*. The best fitting of data gives the scaling laws of relative density for G, D, P, and W structures, and they are given as follows:

The most obvious trend in Fig. 3 is that the relative density increases gradually with an increase in parameter *t*. It indicates the linear relationships for both G- and D-based cellular structures and the quadratic polynomial relations for both P- and W-based cellular structures. According to the error analysis of the fitting, *R*^{2} values are very close to 1.0, so the fitting accuracy of Eqs. (7)–(10) is very high.

### Generating Graded Cellular Structures.

*t*is modified as a function:

*t*(

*x*,

*y*,

*z*) in

*R*

^{3}. Therefore, the TPLS-based density graded cellular can be expressed as follows:

*t*(

*x*,

*y*,

*z*) determines the density of a certain point (

*x*,

*y*,

*z*) ∈

*R*

^{3}, and it can be represented with an expression of high degree to achieve more complex density distribution. In this work, the calculated stress vector or temperature vector is transferred as a 3D matrix

*ρ*_{i,j,k}to generate a complex density distribution. According to Eqs. (7)–(10), the value of each parameter

*t*corresponding to the density three-dimensional matrix

*ρ*_{i,j,k}can be obtained, and then the parameter matrix

*t*_{i,j,k}can be established in 3D space. For example, Fig. 4 shows an example of a linear grading for G, D, P, and W unit cells where density varies from 0.1 at the right to 0.7 at the left of the cellular structures. In addition, material properties such as stress, Young’s modulus, and heat distribution are functions of their relative density in

*R*

^{3}[45]. Accordingly, if the mechanical property of a point in space is determined, the density value of the point can be obtained. The density values for all points can form a three-dimensional matrix, which corresponds to the values in

*t*(

*x*,

*y*,

*z*). Therefore, functional gradient structures that meet the mechanical or thermal properties can be generated by adding this 3D matrix into Eqs. (11)–(14). In the next step, the developed parametric functionally graded TPLS-based cellular generation algorithm will be combined with structural optimization techniques.

## Homogenized Model of TPLS-Based Cellular Structures

In this work, the effective properties of lattice material can be calculated by the asymptotic homogenization (AH) approach [46]. The AH method with a sound mathematical foundation has been widely used to calculate the effective properties of the periodic microlattice structure, the random porous structure, composite material, and also density-based topology optimization [46]. The AH approach is based on the analysis of a representative volume element of the cellular material at mesoscale. As shown in Fig. 5, it is the basic principle of the AH method. Suppose a body Ω* ^{ε}* is infilled with a periodic cellular structure subjected to the traction

**at the traction boundary**

*t**Γ*, a displacement

_{t}**at the displacement boundary**

*d**Γ*, and a body force

_{d}**is distributed in the domain. The cellular structure is considered to be a continuum solid domain Ω with equivalent properties when the AH method is used. Therefore, the computational cost is significantly reduced compared with the full-scale simulation of the explicitly represented cellular structure.**

*f*### Effective Mechanical Elasticity Tensor.

*u*is the superposition of the macroscale quantity

*u*

_{0}(

*x*,

*y*) and a small periodically fluctuated mesoscale quantity

*u*

_{1}(

*x*,

*y*), which can be represented by the first-order asymptotic expansion [46]:

*u*

_{0}is the average displacement on the macroscopic scale, and

*u*

_{1}and

*u*

_{2}represent the perturbation displacement on the microscale [46]. Homogenized elasticity tensor $CTPLSH$ of a periodic unit cell can be rewritten in an equivalent discrete form in terms of elemental energy:

*n*finite elements. The

*V*is the unit cell volume,

*V*is the volume of element

_{e}*e*,

**I**is a six times identity matrix,

**B**

*is the element strain-displacement matrix,*

_{e}**χ**

*is the matrix containing the element displacement vectors, and*

_{e}**D**

*is the constitutive matrix for the element, which for the basic solid material of the object. The global stiffness matrix*

_{e}**K**of the unit cellular microstructure can be described as follows:

**f**

*that corresponds to macroscopic volumetric straining can be expressed as follows:*

^{ij}**ɛ**

^{ij}are chosen to be

**χ**

^{ij}is the global displacement vector of the TPLS-based cellular structures. When the global displacements have been obtained, the homogenized constitutive matrix $CijklH$ can be found as follows:

### Effective Thermal Conductivity Tensor.

*n*is the number of finite elements of the discretized unit cellular structure,

*V*is the volume of the finite element,

_{e}**I**is a six times identity matrix, $Bet$ is the element temperature gradient matrix,

**T**

*is the matrix containing the element nodal temperature vectors, and*

_{e}**k**

*is the element thermal conductivity tensor. For the thermal finite element analysis method, the thermal stiffness matrix*

_{e}**K**

*of the unit structure can be expressed as follows:*

^{t}**t**

*are chosen to be*

^{i}**T**

^{i}is the global temperature vector of the TPLS-based cellular structures. When the temperature vectors have been obtained, the entries in the homogenized thermal constitutive $kijH$ can be found as

*V*is the unit cell volume. $ket$ is the thermal conductivity tensor for each element

_{e}*e*. The tensor is a symmetric matrix of the form

### Scaling Laws of the Tensors.

Relative density *ρ*_{r} is an important characteristic of cellular structures. The “Gibson–Ashby” model [45] suggests that the relative elastic tensor of a cellular structure is a function of its relative density. This function is known as “scaling law.” This law can equally be applied to TPLS-based cellular structures [16,45]. Based on this theory, the relationship between the elastic or thermal tensors, and their relative density of the TPLS-based cellular structures is fitted, and this relationship is referred to as “tensor scaling laws” in this work.

*ρ*

_{r}(

*ρ*

_{r}=

*ρ*

_{c}/

*ρ**, where

*ρ*

_{c}is the density of the cellular structure and

*ρ** is the density of the solid material) is the relative density of a cellular structure and

*a*_{1},

*a*_{2}, and

*a*_{3}are the constant symmetric matrices to be determined. For the anisotropic case, Eq. (30) gives the relationship of the mechanical properties and relative density. In the structure optimization procedure, the elastic tensor and the thermal conductivity tensor scaling laws become one when

*ρ*

_{r}is fully solid and becomes zero when

*ρ*

_{r}is void. Therefore, the fitted function passes two points: (1,1) and (0,0). Hence, the tensor scaling laws of cellular structures can be modified as

The constant symmetric matrices *a*_{1} and *a*_{2} are obtained after fitting them to the homogenization calculation results. The scaling laws of the G, D, P, and W-based cellular structures are plotted in Fig. 6. The results show that the G-, D-, and W-based cellular structures have a relatively strong shear resistance, and the P-based cellular structure has a relatively strong pressure resistance. Therefore, the G, D, and W structures are suitable for optimized structural design with shear resistance requirements, while the P structure is suitable for optimized structural design with compressive requirements. In addition, the heat conduction performance of the G structure is better than that of other structures.

According to the calculated results, the elastic tensor and the thermal conductivity tensor scaling laws for the two constants of TPLS-based cellular structures are presented in Table 1, where $C11111H\u2217$, $C1122H\u2217$, $C1212H\u2217$, and $k11H\u2217$ are the corresponding elastic constants and thermal conductivity constants of the composing material.

## Structural Optimization

### Objective Functions.

General topology optimization methods provide an optimal solution for retaining useful material within a prescribed design domain to achieve the required structural or thermal performances. The elements are considered as isotropic materials. However, almost all types of cellular structures belong to cubic symmetric geometry. The cubic symmetric geometric structure is a highly symmetrical geometrical structure. Furthermore, their mechanical properties are usually anisotropic and are symmetrical only in the direction of the three axes. Therefore, when the anisotropic cellular structure has a design material, Liu et al. [48] proposed a porous anisotropic material with penalization topology optimization method. Wang et al. [49] proposed a topological optimization with connectable graded microstructures, and the anisotropic properties of the microstructure are taken into account. Besides, Liu et al. [50], based on the MMC/MMV topology optimization framework, designed an orientation controllable lattice structure with an anisotropic mechanical performance. In this section, the minimum compliance cellular structure is obtained by adding the calculated elastic tensor or thermal conductivity tensor scaling laws into the optimization algorithm. Therefore, the structure optimization problem can be treated as the optimization of density distribution of the cellular structures of a mechanical part.

*c*(

**) is the objective function,**

*ρ***is the vector of the element relative densities in the design domain and the element densities that constitute**

*ρ***are the design variables,**

*ρ***U**is the displacement vector, $KTPLS(\rho )$ is the global stiffness matrix, and

**f**is the external load vector.

*ρ*is the relative density of element

_{e}*e*, while

*v*is its corresponding element volume. The design variable ranges from

_{i}*ρ*

_{min}to

*ρ*

_{max}. The first constraint is the effective elastic tensor scaling law, which is established in Sec. 5. The second constraint is the Hooke’s law equilibrium equation. The third constraint limits the total design volume to

*V*. The final constraint is the minimum and maximum design variables for the relative density. The optimization problem is solved by the optimality criteria algorithm in this work.

*c*(

^{t}**) is the objective function,**

*ρ***is the vector of the element relative densities in the design domain, $KTPLSt(\rho )$ is the global conductivity matrix,**

*ρ***T**is the global nodal temperature vector, and

**f**

^{t}is the applied heat load vector.

*ρ*is the relative density of element

_{e}*e*, while

*v*is its corresponding element volume. The design variable ranges from

_{i}*ρ*

_{min}to

*ρ*

_{max}. The first constraint is the thermal conductivity tensor scaling law, which is established in Sec. 5.2. The second constraint is the Fourier’s law equilibrium equation. The third constraint is the limit design volume. The final constraint is the minimum and maximum design variables for the relative density.

### Relative Density Mapping Method.

*X*∈

*R*

^{3}is given with target volume constraint, boundary conditions, and the cellular structure

*x*∈

*R*

^{3}used for infill. Second, the effective cellular structure mechanical properties P(

*X*) can be calculated via map

*h*. In this work, map

*h*represents the numerical homogenization method. Third, the established scaling law of cellular structures is added in the optimization procedure, and this process is called mapping

*g*. Therefore, the target density distribution of optimized material property fields P(

*Y*) is achieved, and it is converted into a three-dimensional matrix. In the $hmap\u22121$ mapping process, the density of each element in the optimized output is assigned to all nodes corresponding to that element. Then, the resulting density on each node is calculated by using Eq. (34) to get the average density for all elements of the shared pattern.

*ρ*

_{ijk}is the synthetic density of the node point (

*i*,

*j*,

*k*) and $Xijk\u2217$ is the density of the elements that share the node point.

*N*

_{node}is the number of elements sharing the node (

*N*

_{node}can be 1, 2, 4, 8). Thus, Eq. (34) provides the average mapping density at each element node. After the density at the element node is obtained, the number of cells filled up in each direction can be controlled by linear interpolation across the entire element volume along three dimensions. For example, interpolation is performed between two points (

*x*,

_{a}*y*,

*z*) and (

*x*,

_{b}*y*,

*z*) adjacent in the

*x*-axis direction, and the corresponding densities are

*ρ*

_{x}_{a}and

*ρ*

_{x}_{b}. Then, the distance between the two points is divided into

*n*points {

*x*

_{1},

*x*

_{2}, …,

*x*

_{n}}, each point is (

*x*,

_{i}*y*,

*z*), and the density of each interpolation point is

In the same way, arbitrary density interpolation can be performed throughout the design space to control the number of cells in a certain direction. Lastly, the matrix can be exported to the FGCS modeling method in Sec. 4.1, and the resulting structure *Y* can be generated. By combining the mapping *h* with the mapping *g*, a comprehensive mapping method *f* = *hgh*^{−1} is established. Therefore, the finally optimized results can be directly obtained.

## Results and Discussion

In this work, the proposed gradient cellular structures design and optimization technique is programmed by the matlab code. To illustrate the effectiveness of the optimized gradient cellular structures herein, different design parameters are analyzed, comparison with existing typical gradient lattice design methods is made, and structural robustness is discussed, separately. The gyroid cellular structure is selected as a representative for the design and optimization in this section. For mechanical performance, it is assumed that the Young’s modulus for the solid material that constructs the cellular materials is 2.15 Gpa and the Poisson’s ratio is 0.3. For thermal conductive performance, for the solid material, the Young’s modulus of 60 GPa, the density of 2.5 × 10^{−6} kg/mm^{3}, and the Poisson ratio of 0.35, the thermal conductivity k = 0.235 W/(mm K) are assumed.

### Size-Effect Estimation.

The homogenization theory is based on a periodic distribution of cell structures that are much smaller than the macroscopic scale. The ideal condition is that the microscale unit size is much smaller than the macroscopic scale. However, when considering the results of actual optimization to meet the minimum size of existing AM, we need to find a minimum cell size that best satisfies the homogenization theory. In this paper, the influence of the size and boundary of the cellular structure on homogenization accuracy is discussed through numerical studies.

As shown in Fig. 8(a), a cubic block of 30 × 30 × 30 mm is employed as a benchmark with a nonperiodic boundary to study the influence of element size. The fixed boundary conditions are subjected to the bottom of the cube, while a load of 1000 N is applied to the upper surface while the other surfaces are unstressed, as shown in Fig. 8(b). For thermal conductivity, the temperature is 20° C at the top surface, a heat flux of 0.1 W/mm^{2} is applied to the bottom surface, and all other surfaces are perfectly insulated in Fig. 8(c). The cube is infilled with a gyroid-based cellular structure and the volume fraction is 0.3, and the constitutive model in Eq. (31) is used in the simulation. The analysis process is performed using the OptiStruct 2017 software. Accordingly, the maximum total deformation of the homogenized model is 0.092 mm and the maximum temperature of the homogenized model is 121.5° C.

To test the influence of the unit size on the nonperiodic boundaries in this benchmark, the number of unit cells in one direction is set as 1, 2, 3, 4, 5, 6, 8, 10, and 12 (Figs. 8(d)–8(l)). Then, a full-scale simulation based on these models is performed to compare with the homogenized model. Figure 9(a) shows the maximum deformation of the upper surface relative to the number of unit cells in one direction. It can be seen from the results that the number of unit cells along one direction is larger than six; then the maximum total deformation gradually converges to the homogenized result. Besides, as can be seen in Fig. 9(b), as the unit cell size decreases, the temperature distribution of the bottom surface gradually approaches the homogenized result. In particular, when the number of unit cells in the heat transfer direction is larger than six, the difference between the homogenized model and the full-scale simulation is less than 2%. Therefore, six-unit cells along one direction of the design component are considered to be the minimum number in this work. In summary, the cell size and the relative density ranges have been appropriately determined to maximize the accuracy of the AH model while addressing manufacturability issues.

### Structural Mechanical Performance

#### Example.

In this section, the classic 3D cantilever beam design problem is solved to demonstrate the effectiveness of the proposed methodology. Figure 10(a) shows that the dimensions of the beam model are 60 mm × 20 mm × 20 mm. The loading and constraints are also indicated in Fig. 10(a). The design domain is discretized using a 300 × 100 × 100 uniform grid. A mesh of 15,125 eight-node hexahedron elements is used to discretize the design domain and solve the FEA model. A load of 100 N is initially distributed over the object, as shown in the figure. The selected material parameters are the modulus of elasticity *E* = 2.15 GPa and Poisson’s ratio *v* = 0.3. The elastic tensor scaling law of the gyroid-based cellular structure applied to the optimization procedure is expressed in Eq. (30), and the correlation coefficients can be queried in Table 1. The target volume fractions are set to 0.3, 0.4, and 0.5, respectively. In addition, the minimum and maximum relative densities are set as uniform in Fig. 10(c), set as *ρ*_{min} = 0.1 and *ρ*_{max} = 0.9 in Fig. 10(e), and set as *ρ*_{min} = 0.1 and *ρ*_{max} = 0.8 in Fig. 10(f).

Figure 10(b) shows the optimized cellular structure by using gyroid to infill the design domain. The optimized structure has several notable characteristics. First, compared with the traditional topological optimization, the cellular structure is filled in the whole design area, which makes the overall structure robust, which will be discussed in the following sections. Second, the cellular structure is not only the endurance structure but also can be used as the support structure of the 3D printing process to improve manufacturability. Finally, it can be seen from the cross section area that the cellular structures with different densities can be smoothly connected (Fig. 10(c)), which means the smooth continuous structure can effectively transfer mechanical properties. In particular, the effects of different filling densities and different density intervals on the design results are also illustrated in Figs. 10(d)–10(f). The main conclusions are as follows: First, according to Figs. 10(e) and 10(f), the stiffness of the optimized structure with a density of 0.1–0.9 intervals is better than that of the density of 0.1–0.8 intervals. Second, by comparing Fig. 10(d) with Figs. 10(e) and 10(f), it can be seen that the optimized filling method makes the overall structure to have a better stiffness than that of uniform structures. For example, Fig. 11 illustrates the stiffness obtained by homogenization result and the full-scale simulation of Case #1 and Case #2. in Fig. 10. Table 2 lists the corresponding values of the stiffness of the linear elastic region for the nine cases in Fig. 10. From Fig. 11 and Table 2, it can be seen that the predicted stiffness value from the homogenization model is closer to the value of the full-scale simulation, and the stiffness of our optimized results is significantly better than that of the uniform form. This means that the proposed homogenization model can be effectively used to obtain a density variable lattice structure for optimizing the stiffness of an object.

#### Robustness.

The cellular structures possess an advantage of being able to tolerate localized damage. Parts in the manufacturing or working process are prone to local damage caused by stress concentration or an accidental event. These forms of damage mainly include accidental collision, fatigue failure, corrosion, and manufacturing defects. With traditional topology optimization, the structure retains the stressed material, and the rest is deleted. Therefore, when the optimized results are partially destroyed, the mechanical properties of the global structure will be seriously affected. The advantage of the density optimized cellular structures used in this work is that even if parts are locally broken, the overall stiffness can persist.

First, the half Messerschmitt-Bölkow-Blohm (MBB) beam (60 mm × 20 mm × 20 mm) is optimized by three different optimization strategies with the same materials and boundary conditions. They are the SIMP-based topology optimization method, the uniform gyroid-based cellular lightweight design method, and our proposed method, respectively (as shown in Fig. 12(a)). From the results, we can see that the stiffness of the topology optimization method is the best, and our method is slightly smaller than it. Moreover, compared with the uniform design method, the above two approaches have apparent advantages in stiffness. However, the SIMP-optimized results need much more support structures (vertical structures in Fig. 12(a)) for the AM process.

Second, to demonstrate the effectiveness of damage tolerance of distributed cellular structures, several simplified local deficiencies for the design results in the same place are added, as shown in Fig. 12(b). Subsequently, FEA simulations are completed separately for the original optimized models and the models with the inflicted flaws in different areas. The analysis process is performed using the OptiStruct 2017 software. The material properties used in the simulation are as follows: *E* = 2.15 GPa, *v* = 0.3, and *F* = 100 N. The maximum displacement of these models is counted in Fig. 12(c). For example, the displacement of the damaged structures (deficiency in A) and the nondamaged ones is compared; the displacement changes by 192.3% for the SIMP method, while it changes only by 33.6% for our optimized infill method. Furthermore, from the deformation results caused by different damaged areas, it can be seen that the displacement of the SIMP method fluctuates considerably. However, our proposed FGCS-based optimization method has less influence on its displacement under different failure modes. Although the displacement of a uniform structure under different failure modes is less affected, the overall deformation is always larger than that of other design strategies. The above result indicates that the proposed infill optimization technique produces a less damage-sensitive design than that of the standard SIMP method.

#### Comparison.

In this section, the proposed elastic tensor scaling law added optimization method will be compared with the cellular relative density mapping based on topology optimization with isotropic material [10,25]. In addition, it will be compared with the truss structure-based relative density optimization method [28,29]. The above methods correspond to “Method-1”, “Method-2”, and “Method-3” in Table 3, respectively. The first column of Table 3 is the design domain with the boundary conditions, as well as the items that need to be compared. The three models of Method-1, Method-2 and Method-3 are optimized with the same material properties, volume constraints and the compliance values for the two iterations are less than 0.01 as the same convergence condition. Besides, the maximum stress is obtained by FEA of the different optimized results.

According to the compliance values, it indicates that our proposed method is 1.35 times stiffer than Method-2 and 1.12 times stiffer than Method-3. Moreover, from the results of the finite element analysis, it can be seen that our proposed method has a maximum von Mises stress that is 52.6% of the maximum stress in Method-2 and 32.9% of the maximum stress in Method-3. This occurs because there is no direct material property relationship between the optimization calculation process in Method-2 and the filled cellular structure. In addition, the geometrical shape in Method-1 is smoother and more continuous than Method-3. Therefore, Method-1 is not an easy-to-produce stress concentration phenomenon. The above results show that the proposed method has better structural performance.

### Thermal Conductive Performance.

Heat sinks are passive heat exchangers that are used widely in electronic devices to cool key components such as central processing units or graphic cards in computers. In this section, several heat-sink designs are presented. Their thermal conductive behaviors are compared. The initial geometry and boundary conditions are given in Fig. 13(a). It consists of a design domain 180 × 180 × 30 mm and has a sink on the middle of the top surface. A uniform surface heat flux of magnitude 0.1 W/mm^{2} is applied on the top surface and a forced convection exists over the structural boundaries. The optimization goal of the heat-sink is to minimize the thermal compliance satisfy a volume constraint of 30% of the design domain. For the solid material, the Young’s modulus of 60 GPa, the density of 2.5 × 10^{−6} kg/mm^{3}, and the Poisson ratio of 0.35 are assumed. Other parameters are set as follows: thermal conductivity k = 0.235 W/(mm K) and prescribed temperature T_{0} = 20° C on the bottom, left, and right boundaries. Figure 13(b) illustrates a reference heat-sink design with uniform cellular structures and the corresponding thermal conductive simulation result. It has a maximum temperature of 66.1 °C. Figs. 13(c) and 13(d) show the cellular relative density distribution optimized results by our proposed method. The minimum and maximum relative densities are set as *ρ*_{min} = 0.1 and *ρ*_{max} = 0.8 in Fig. 13(c) and as *ρ*_{min} = 0.1 and *ρ*_{max} = 0.9 in Fig. 13(d). It can be seen from the simulation results that the optimized results are superior to the thermal conductivity of the homogeneous structure. Besides, in the optimal design results of the same volume fraction, the density range of 0.1 to 0.9 has better heat conduction performance than the density range of 0.1–0.8. Furthermore, the results of the homogenization method are very close to the full-scale simulation results in Table 4. It implies that the proposed homogenization model can be effectively used to obtain a density variable lattice structure for optimizing the temperature distribution.

## Conclusions

In this work, a generative design algorithm combining the TPLS-based graded cellular structure modeling method with structure optimization is developed to optimize the parts for additive manufacturing. In contrast to the traditional uniform cellular design and standard topology optimization method, our technique not only optimizes the stiffness of the structure, but also makes the optimized result structural robustness concerning damage tolerant. Besides, the thermal conductivity of the cellular structures can also be effectively optimized by our method. The conclusions of this paper are as follows:

We presented an FGCS modeling methodology based on the TPLS-based cellular structures. It can be used to map material properties into 3D space.

A structure optimization strategy was established to achieve the optimum density distribution of cellular structures. The scaling laws of the elastic tensor (thermal conductivity tensor) of the TPLS-based cellular structures were calculated by using the numerical homogenization method, and it was added to the optimization process. This allows the optimized relative density distributions of the cellular structures distributed in a component to possess an accurate physical meaning.

In future work, we will expand our approach to the study of multiple physical coupling optimizations by using gradient cellular structures.

## Acknowledgment

This work is supported by financial support from the Program of China Scholarships Council (No. 201706830039); National Science Foundation of China (Grant No. 51775273); Natural Science Foundation of Jiangsu Province, China (No. BK20161487); and Six talent peaks project in Jiangsu Province, China (Grant No. GDZB-034).