A novel lattice structure topology optimization method with extreme anisotropic lattice properties

This research presents a lattice structure topology optimization (LSTO) method that significantly expands the design space by creating a novel candidate lattice that assesses an extremely large range of effective material properties. About the details, topology optimization is employed to design lattices with extreme directional tensile or shear properties subject to different volume fraction limits and the optimized lattices are categorized into groups according to their dominating properties. The novel candidate lattice is developed by combining the optimized elementary lattices, by picking up one from each group, and then parametrized with the elementary lattice relative densities. In this way, the LSTO design space is greatly expanded for the ever increased accessible material property range. Moreover, the effective material constitutive model of the candidate lattice subject to different elementary lattice combinations is pre-established so as to eliminate the tedious in-process repetitive homogenization. Finally, a few numerical examples and experiments are explored to validate the effectiveness of the proposed method. The superiority of the proposed method is proved through comparing with a few existing LSTO methods. The options of concurrent structural topology and lattice optimization are also explored for further enhancement of the mechanical performance.


Introduction
With the rapid development of advanced material design and manufacturing technology, lattices have been increasingly employed in structural design for their excellent properties, such as high specific stiffness and energy absorbability (Zhang et al., 2019a(Zhang et al., , 2020a, negative Poisson's ratio (Alomarah et al., 2020), negative thermal expansion coefficient, and many others (Han et al., 2015;Takezawa et al., 2015;Takezawa & Kobashi, 2017). The complex micro-or mesoscale geometry of the lattices con-tributes to the above-mentioned outstanding physical properties. Hence, design of the lattice geometry, especially the varying lattice shape and topology infilling a macroscale part, deserves careful investigations (Xia & Breitkopf, 2017;Chen et al., 2019;Xing et al., 2020).
Topology optimization, as a generative structural design method, features in searching the optimal material distribution within a given design domain to maximize the structural physical performance (Bendsøe & Sigmund, 2004). The design target is not restricted to the macroscale material distribution , and a more attractive feature lies in the support for twoscale structure design, i.e. concurrently optimize the macroscale part and the micro-or mesoscale representative material unit (Li et al., , 2021. Hence, lattice structure topology optimization (LSTO) has been extensively studied in the past few decades Groen et al., 2020;Wu et al., 2021;Xiao et al., 2021).
At the early stage, the research on LSTO focused on monoscale structure design. An inverse homogenization method (IHM) is proposed to design the representative material units to achieve specific constitutive parameters (Sigmund, 1994). Based on IHM, extensive researches have been carried out on material-level topology optimization to realize specific or extreme constitutive properties that are normally nonachievable in nature. Sigmund and Torquato (1997) introduced a three-phase topology optimization method to design composites with extremal thermal expansion coefficients. Guest and Prévost (2006) optimized multifunctional materials for maximized stiffness and fluid permeability. Metamaterials with negative Poisson's ratio (Wang et al., 2014;Alomarah et al., 2020), functionally graded properties (Radman et al., 2013), maximum bulk or shear modulus (Huang et al., 2011), and so forth have been thoroughly investigated.
To enhance the macroscale part performance, other than the monoscale metamaterial design, research attention has been drawn to the field of multiscale topology optimization. Generally speaking, under the assumption of scale separation, multiscale topology optimization targets at optimizing the geometric structure in macro-and microscales simultaneously. One typical hypothesis of the multiscale design methods is that the entire macroscale part is composed of homogeneous microstructures; i.e. the same representative material unit repetitively infills the part. Liu et al. (2008) presented a concurrent topology optimization method to simultaneously optimize both the macroscopic material distribution and microscopic representative material unit. The optimized material units consisted of anisotropic truss-like structures. Subject to similar concurrent topology optimization methods, the fundamental frequency optimization (Niu et al., 2009;Vicente et al., 2016), maximum stiffness optimization Gao et al., 2019a), robust optimization (Guo et al., 2015;Zheng et al., 2019), nonlinear solid behavior optimization , and multi-objective optimization (Yan et al., , 2016 have been widely studied. Due to the assumption of homogeneous microstructures, the computational time is not significantly increased compared with the traditional macroscale-only topology optimization, but the great design space of multiscale topology optimization enabled by the varying microstructures cannot be explored yet. Then, new assumptions arise about the macrostructure consisting of heterogeneous microstructures. Apparently, the varying microstructures provide more flexible material options to address the local loading conditions, and thus result in better mechanical performance enhancement Nguyen et al., 2021). The hierarchical model for concurrent material and topology optimization of 2D structures and threedimensional (3D) structures were proposed by Rodrigues et al. (2002) and Coelho et al. (2008), respectively. Then, Xia and Breitkopf (2014) constructed an FE2 nonlinear multiscale analysis framework for concurrent design of the materials and structures. Sivapuram et al. (2016) presented a trial-and-error criterion to divide the geometrical domain into subdomains and have each subdomain composed of a unique microstructure. Xu and Cheng (2018) proposed a two-scale concurrent topology opti-mization method with multiple microstructures selected based on the principal stress orientation. Wang and Kang (2019) developed a novel multiscale topology optimization framework for composite structure design by combining the velocity field level set method and the Solid Isotropic Material with Penalization (SIMP) method. Xu et al. (2021) contribute a multiscale topology optimization method that includes both the shell and latticelattice interface layers to ensure the connectivity between different microstructures. More researches on the optimization of macroscale part composed of heterogeneous microstructures can be found in Li et al. (2016Li et al. ( , 2018b, Gao et al. (2019b), andZhang et al. (2020b, c).
On the other hand, multiscale topology optimization with heterogeneous microstructures has several conspicuous limitations. First is the huge computational cost, since the elementwise homogenization is needed for each iteration for every macroscale element; second is the manufacturing difficulty due to the disconnectivity issue between adjacent heterogeneous microstructures (Zhang et al., 2019b;Liu et al., 2020); the last issue lies in the intensive post-processing cost to generate an effective CAD model and the huge size of the generated CAD model (in Gigabytes) that altogether hinder the manufacturing attempts.
In recent years, the design of concurrent multiscale structures with parametrized lattice microstructures is also popular (Ypsilantis et al., 2021). Li et al. (2018a) presented an LSTO method for designing functionally graded gyroid lattice structures. Cheng et al. (2019) proposed a functional gradient LSTO method for additive manufacturing (AM) components with stress constraints. Wang et al. (2018Wang et al. ( , 2020 optimized 2D and 3D gradient lattice structures composed of cubic and X-shaped lattices, in which the two lattices are combined inside each material unit to expand the design space. Wu and Li (2020) employed the extended multiscale finite element method to evaluate the effective properties of lattice microstructures and then developed the conformal LSTO method. Yu et al. (2020) developed a strength-constrained multiscale topology optimization method by simultaneously incorporating the von Mises and Tsai-Hill failure criteria. More importantly, the numerous homogenization processes are performed beforehand to construct the parametrized material constitutive model that significantly reduces the in-process computational burden of multiscale topology optimization. Meanwhile, connectivity between the heterogeneous lattices is guaranteed due to consistent lattice topology because in most cases, only the lattice densities are varying. The above advantages have led to many successful industrial application examples of LSTO. However, the majority of existing LSTO studies focus on parametrized optimization attempts with the classic lattice unit types, e.g. cubic, Xshape, cubic center, gyroid, and some others. Given most of the selectable lattice types, only varying the lattice densities or geometric parameters is insufficient to provide enough property variations to address the complex loading conditions spanning across a macroscale part, which severely restricts the structural performance level achievable through LSTO.
To address the above issues, this paper proposes a novel strategy to generate three groups of sample lattices through topology optimization that magnify the tensile modulus in the x direction, the tensile modulus in the y direction, and the shear modulus in the x-y plane, respectively. Each group includes topologically optimized lattices with different material volume fractions. Then, the candidate lattice for LSTO is constructed through combination of the topologically optimized sample lattices; i.e. one lattice sample is selected from each group and then, the three selections are combined through union Boolean operations to form the candidate lattice for LSTO. This novel candidate lattice is controlled by the density ratios of the three composing lattices, and hence, its material constitutive model can be empirically represented by these density ratio parameters through sampling and interpolation. In this manner, the accessible properties of the candidate lattice are greatly expanded since ratios of the three composing lattices can be flexibly arranged. Therefore, parametrized LSTO can be performed with greatly expanded anisotropic effective property ranges while not requiring the repetitive homogenization calculations, leading to ever enhanced structural performance improvement for LSTO.
The remainder of this paper is organized in the following manner: Section 2 introduces the energy-based homogenization method. Section 3 presents the idea of combining topologically optimized lattices to form the novel parametrized candidate lattice with extremely large property ranges. Meanwhile, the LSTO problem is formulated and the related sensitivities are derived. Sections 4 and 5 demonstrate the numerical examples and mechanical testing results, respectively, to show the effectiveness of the proposed method. LSTO of two practical parts and their AM are showcased in Section 6. Section 7 explores the concurrent structural topology and lattice optimization to further improve the structural stiffness. Discussions and conclusions are provided in Section 8.

Homogenization Theory
Homogenization refers to a method that can replace the composite with an equivalent material model in order to resolve the difficulty in analysing the boundary value problem with high heterogeneities (Hassani & Hinton, 1998;Dong et al., 2019). The effective elasticity tensor D H i jkl of the unit cell can be calculated through homogenization (Cheng et al., 2013) as where D pqrs is the locally varying elastic properties. Y is the volume of the representative unit cell, ε 0(i j) pq is the linearly independent unit test strain field, and ε * (i j) pq denotes the unknown microscale strain field. In finite element analysis, the homogenized effective elasticity tensor D H can be reformulated through summing up the element-wise calculations through the following expression (Gao et al., 2018): where u * (i j) e is the elemental displacements under periodic boundary conditions, and k e is the stiffness matrix of the finite element discretizing the representative unit cell.
With the energy-based homogenization method, the unit test strains are directly imposed on the boundaries of the periodic unit cell, and the effective elasticity properties are interpreted as the summation of element-wise strain energies. Then, the superimposed displacement field (u 0 e − u * e ) can be transformed into the induced displacement field u n e (Gibiansky & Sigmund, 2000): Hence, the effective elasticity tensor D H is calculated as Topology optimization can be applied to design the periodic unit cells for extreme effective properties (Gao et al., 2018). The objective function is stated as where J denotes the objective function and λ i jkl represents the weight factor. For instance, in the 3D case, the maximization of the material bulk modulus corresponds to and the maximization of the material shear modulus corresponds to Therefore, the periodic unit cells with extreme effective properties can be designed by topology optimization; see Fig. 1 for two examples.

Design of lattice structures with extreme mechanical properties
As illustrated in the last section, the optimal lattice structures with extreme properties (such as maximum tensile modulus or shear modulus) can be designed through topology optimization.
In this research, we focus on the lattice properties within the x-y plane and therefore, three groups of lattices subject to different volume fraction constraints are designed through topology optimization to have the maximum tensile modulus in the x direction, the maximum tensile modulus in the y direction, and the maximum shear modulus inside the x-y plane, respectively, as demonstrated in Fig. 2. Meanwhile, the initial designs for elementary lattice optimization are identical and only the volume upper limits are varied for each group of lattices. This setup ensures the connectivity of the varying-density elementary lattices, especially for the shear-stiff elementary lattices. Then, the connectivity of the candidate lattices formed through union Boolean operations of the elementary lattices is also guaranteed.
The relative density ρ λ , which denotes the ratio of the optimized cellular structure volume to the solid cellular structure volume, can be expressed as where V domain is the total volume of the unit cell and V opt is the volume of the optimized lattice structure. The relative densities ρ λ of the three groups of optimized lattices are described as: ρ x for the optimized lattices with maximum tensile modulus in the x direction, ρ y for the optimized lattices with maximum tensile The effective elasticity matrix modulus in the y direction, and ρ xy for the optimized lattices with maximum shear modulus in the x-y plane. Within a macroscale part, the stress state varies across the material domain and thus, it is significant to diversify the lattices' effective properties to address the varying local loading conditions. Therefore, a novel candidate lattice is developed through a selective combination of lattices from the abovedeveloped three lattice groups. Specifically, one optimized lattice is selected from each group and then, the three selected lattices are combined through union Boolean operations to create a new lattice (see Fig. 3). Apparently, the new lattice can have excellent x-direction tensile modulus, y-direction tensile modulus, x-y plane shear modulus, or even concurrently all of the above. When performing superimposition, the extreme properties of the elementary lattices, such as the tensile or shear stiffness, are not conflicting or diminishing each other because of only adding up material through the union Boolean operations. The excellent superimposition effect is also verified by checking the effective elasticity matrices. Hence, the complex loading conditions inside a macroscale part can be addressed by varying the compositions of the proposed candidate lattice. Note that the proposed candidate lattice may have discontinuous material distribution (refer to the first row in Fig. 3), especially for the shear-dominated lattices, because in the current research, loadings inside the x-y plane are majorly focused while the zdirection tensile modulus is not included when optimizing the composing lattices. Therefore, we propose to add four corner stripes along the z direction to guarantee the connectivity. As shown in Fig. 3, the strips consume a small portion of the total material volume and therefore, do not obviously reduce the property-to-density ratio.

Parametrized property modeling for the proposed candidate lattice
The last subsection illustrates the idea of creating a new candidate lattice through a combination of the topologically optimized lattices. Then, it is necessary to parametrize the new candidate lattice and build its constitutive model, i.e. build the mathematical relationship between the elastic properties of the new lattice and the design variables describing its composition. The relative densities (ρ λ , λ = x, y, or xy) of the composing lattices are employed as the design variables. Without loss of generality, it is assumed that the solid material that constructs the lattices is isotropic. The Young's modulus and the Poisson's ratio of the material are chosen as E 0 = 100 MPa and υ = 0.3, respectively. The effective elastic matrix of the proposed candidate lattice can be obtained by numerical homogenization method as shown in Table 1.
As shown in Table 1, the candidate lattice remains symmetric with the change of density parameters (ρ x , ρ y , and ρ xy ) and there exist nine independent constants in the effective elasticity matrix. For the convenience of expression, the nine independent constants are denoted as follows: The configuration of the new lattice is controlled by the density ratio parameters, and therefore, the elastic constants are developed as functions of the three density ratio parameters. Specifically, second-order polynomials are employed to model the relationship between the effective elastic matrix components and the three relative density parameters. So, D i j can be expressed as D i j = a 1 + a 2 ρ x + a 3 ρ y + a 4 ρ xy + a 5 ρ 2 x + a 6 ρ 2 y + a 7 ρ 2 xy + a 8 ρ x ρ y + a 9 ρ x ρ xy + a 10 ρ y ρ xy , where a i (i = 1-10) denotes the coefficients obtained by fitting the homogenization data. In order to accurately describe the fitting relationship, 1000 sample points are used to fit the polynomial function in (10), i.e. 10 relative density options for each lattice group. The fitting process involves three parameter variables of ρ x , ρ y , and ρ xy (10 sample densities for each) and is accomplished by the least squares-based regression analysis. To describe the volume fraction of the newly developed candidate lattice, a polynomial interpolation function for the density ratio (ρ v ) is also constructed, which could be written as Similarly, b i (i = 1-10) are the polynomial coefficients. In order to show the fitting results, the homogenized values and the fitting values of D i j and ρ v are plotted in Fig. 4, and it is observed that the fitting results coincide well with the theoretical calculations. In addition, the normalized maximum fitting errors are calculated and shown in Table 2. For the majority, the maximum error is smaller than 0.025 while the only exception is 0.0613 for D 33 .
So far, we have obtained the interpolation results of the elastic components D i j and the total density ρ v with design variables (ρ x , ρ y , and ρ xy ). These results would help avoid the repeated numerical homogenization efforts during optimization and thus greatly save computational cost.

Formulation of the LSTO problem
In this paper, the minimum compliance problem of 3D structures is considered. The objective of the minimum compliance problem is to find the material density distributionx that minimizes the structural deformation under the prescribed support and loading condition (Liu & Tovar, 2014). The objective function      could be written as where C denotes the global compliance of the macrostructure. U and F are the nodal displacement and force vectors of the macrostructure, respectively. The global stiffness matrix K(x) of macrostructure can be obtained by assembling the element stiffness matrix k i as follows: where B is the strain-displacement matrix. Then, the topology optimization problem is formulated as below: wherex is the collection of design variables. v i is the volume of the i-th element. The total volume V * has the upper limitV. As the lower and upper values of the design variables, ρ min and ρ max can be set according to the requirement.

Sensitivity analysis
The derivatives of the objective function can be derived by Taking derivatives on both sides of the equilibrium equation: Then, we can derive the following expression: For the design-independent optimization problems as considered, it follows ∂ F /∂x = 0. By substituting (17) into (15), we obtain Combining (13), the sensitivity of the objective function with respect to the design variablex can be expressed as Then, the derivative of the effective elastic matrix D H can be calculated from the quantitative relationship between the elastic matrix constants and the design variables as follows:

Numerical implementation
In this paper, the topology optimization problem formulated in (13) is solved by the method of moving asymptotes (Svanberg, 1987). A sensitivity filtering technique (Sigmund & Petersson, 1998) is imposed in the optimization. The flowchart of the proposed method is illustrated in Fig. 5. The whole procedure involves three stages: construction of the interpolation models, concurrent topology optimization, and establishment of the geometric database of candidate lattice structures for CAD model generation. In the first stage, the new candidate lattice with accessible extreme anisotropic properties is constructed by combing pre-optimized elementary lattices. Then, the equivalent properties of the new candidate lattice are parametrized and estimated through data fitting. In the second stage, concurrent topology optimization is performed to obtain the optimal relative density distributions of the elementary lattices, thus obtaining the composition distribution of the candidate lattices. In the last stage, the optimized candidate lattices are reconstructed by picking up the pre-established geometric models of the preoptimized elementary lattices from the database. The geometric database of elementary lattices uses 0.01 as the density interval, and theoretically more than 90 000 candidate lattices with extreme anisotropic properties can be formed. The part CAD model is built by assembling the 3D binary matrices of the candidate lattices and then transforming it into an STL model with programs available in Vogiatzis et al. (2017).

A 3D three-point bending beam
In the first example, a 3D three-point bending beam is optimized. The design domain and boundary condition are shown in Fig. 6. The dimension of the beam is 48 mm × 12 mm × 8 mm.
A 60 N load is uniformly distributed along the centerline of the top surface. The two foot corners are supported to restrict the vertical displacement. Because of the symmetry condition, only one half of the structure is optimized. A mesh with 2304 (24 × 12 × 8) 8-node hexahedron elements is employed to discretize the design domain. The objective function is to minimize the macrostructure compliance under a total volume constraint of 0.5. Considering manufacturability constraints, the upper and lower limits of the lattice densities are set to be 0.18 ≤ ρ v ≤ 0.93, corresponding to the upper and lower limits of the design variables of 0.06 ≤ ρ λ ≤ 0.5. The initial setup of the optimization problem is the uniform infill with the proposed candidate lattices with ρ v = 0.5 and ρ λ = 0.22. The convergence curves of the structural compliance and total volume fraction are plotted in Fig. 7. It is clear that the optimization process converges rapidly with less than 100 iterations. Slight fluctuations appear since both the microstructure's equivalent constitutive model and relative density are modeled with a highly nonlinear relationship on the design variables (density ratios of the composing lattices), which increases the converging difficulty. Traditionally, for LSTO, the microstructure densities are directly applied as the design variables so that the nonlinearity is not so strong.
The optimized structure has a compliance of 217 mJ by the proposed method. Given the fact that the finite element calculation is based on a fitting mathematical model, the compliance of this optimized design is recalculated by using strictly homogenized elasticity models to access the precision of the fitting function, and the bending beam has a recalculated compliance of 215 mJ, showing a very small deviation from the fitting function-based result. Figure 8 shows the evolution histories of the total density ratio and the design variables. In general, the optimized density distribution is similar to what have been obtained in literature (Liu & Tovar, 2014;Cheng et al., 2019;Wang et al., 2020). As illustrated in the figure, most of the materials distribute around the areas of high strain energy densities. The lattices through optimization are majorly composed of the x-direction tensile-stiff elementary lattices and the shear-stiff elementary lattices, which makes sense given the overall loading condition of the three-point bending beam. To better show the optimized compositions of the lattices, the full-scale beam structure is reconstructed as shown in Fig. 9. In Figs 8 and 9, only the right half of the bending beam is shown because of symmetry.
For the purpose of comparison, three other designs obtained by different methods under the same boundary conditions and volume fraction are given. Table 3 shows the geometric models of these specimens: Model in Case A is uniformly filled with cubic lattice; model in Case B is obtained by a concurrent design method with heterogeneous lattice microstructures composed of varying cubic and X-shaped lattice densities ; the macrostructure infilled with topologically optimized homogeneous microstructures is given in Case C. It is obvious that the multiscale design achieved with the proposed method as demonstrated in Case D has shown the outstanding compliance performance.

A 3D L-bracket structure
In this example, a 3D L-bracket structure is investigated. The dimension and the boundary conditions of the L-bracket structure are shown in Fig. 10. The whole design domain is meshed with 2560 8-node hexahedron elements. Similarly, the optimization problem minimizes the structure compliance with the volume fraction limit of 0.5. The model is clamped at its top surface, and a 30 N load is applied as shown in Fig. 10. Other parameters remain the same as the first example. Figure 11 plots the evolutionary histories of the structural compliance and the total volume fraction obtained by the proposed LSTO method. The compliance curve eventually converges to 399 mJ. Due to the complexity of the involved multivariate design space, slight fluctuations are observed during the convergence process, similar to the first example. Then, the evolution histories of the total lattice densities and the composing lattice densities are shown in Fig. 12. It can be ob- served that inside the bottom half of the final optimized Lbracket structure, the high-density lattices are majorly composed of the x-direction tensile-stiff elementary lattices and the shear-stiff elementary lattices for the bending-dominated deformation behavior. The y-direction tensile-stiff elementary lattices play an important role forming the lattices for the top half of the optimized L-bracket structure. The above observation is consistent with static structural analysis results. The full-scale L-bracket model is reconstructed as shown in Fig. 13.
Meanwhile, for the comparison purpose, numerical designs from the proposed and other optimization methods are presented in Table 4. The structural compliance of Case D designed with the proposed LSTO method is the smallest among these four cases, which verifies the effectiveness of the proposed LSTO method.

Mechanical Testing
In this section, experimental testing is performed on the optimized three-point bending beam model to further validate the effectiveness of the proposed LSTO method. As the shape is highly intricate, the optimized model is manufactured through AM Huang et al., 2020;Zhang & Moon, 2021). The digital light processing (DLP) is one of the AM technologies that uses liquid photopolymer resin as the material and ultraviolet (UV) light to convert liquid resin into solid structures (Gardan, 2019). It shows higher printing accuracy than the extrusionbased polymer AM. Hence, the models are fabricated using a DLP-type 3D printer (Rayshape Shape 1) with its standard white resin to physically realize the optimized structural details. Figure 14 shows the reconstructed CAD models and the 3D printed structures. The size of the unit cells is uniformly configured to be 5 mm × 5 mm × 5 mm and the overall structure size is 120 mm × 30 mm × 20 mm. In order to prevent contingency of the experiments, two samples of each specimen are printed for testing. A universal test machine was used for the bending test, and the punch speed was set as 2 mm/min (see Fig. 14). All data of the experiments from bending to crushing are shown in Fig. 15a. Figure 15b captures the linear stage of the deformation process and moves the starting points of the linear segments to the origin of the coordinate system, which demonstrates the stiffness of the designs from different optimization strategies. It is obvious that the optimized beam in Case D is much stiffer than other beams, which is consistent with the numerical calculation results.
At the end of this section, a special discussion is made on the Case C from multiscale topology optimization wherein constant-thickness plate-based unit cells are derived. In fact, for 2.5D loading conditions, the plate-based representative unit cells (Fig. 16a) perform excellently given both the tensile and shear properties, wherein the effective elasticity model is almost linearly dependent on the relative density. Hence, performing variable-thickness optimization based on the plate-based unit cells can achieve even better structural performance. Referring to Fig. 16a, the derived structural compliance is 182 mJ that is smaller than that of Case D in Table 3. Apparently, the proposed

Case
Compliance Macrostructure Microstructure candidate lattice, made of topologically optimized elementary lattices, can get close to but can hardly reach the same performance level of the plate-based unit cell, owing to the local optimum issue. On the other hand, the plate unit cell-based optimization result falls short in buckling strength. As shown in Fig. 16, the optimization result is transformed into an equivalent CAD model for testing, and to guarantee the structural integrity, four bars are added to connect the two separated plates inside the unit cell. Note that the bar structures are not counted into the total volume when performing optimization and two diameter options are considered of 0.5 mm (noted as Case S1) and 1 mm (Case S2), respectively. Then, the mechanical testing results are shown in Fig. 16c and d. It is observed that the structural stiffness values of both Case S1 and Case S2 are better than those of Case D; however, the structural strengths of Case S1 and Case S2 are only 1/4 and 1/2 of Case D because of the buckling failures. The strength could be even worse if not involving the connection bars. Hence, in conclusion, LSTO results with the proposed candidate lattice do significantly enhance the structural stiffness performance compared with conventional LSTO results, and the resulted structural strength is strong as well, especially compared with the plate unit cell-based structures.

Applications to Part Design
In order to illustrate the applicability of the proposed LSTO method for practical cases, this section gives a design case of an automobile suspension control arm. The design requirements are as follows: (1) The total weight of the part should be reduced by at least 20%.
The complete design process can be divided into five steps as shown in Fig. 17. The first step is to identify the design domain of the part. Then, the design domain is voxelized for finite element analysis and to incorporate the element-wise design variables. Next, conduct the proposed LSTO on the voxelized structure to reduce the weight by 40%, and according to the optimized lattice density distribution, generate the CAD model for the optimized lattice structure as shown in Fig. 17d. Then, postprocess the CAD model to add a constant-thickness solid layer to accurately remodel the part shape while eliminating the zig-zag lattice boundary. The final design result after CAD model regeneration is shown in Fig. 17e. It measures 48 750 mm 3 in volume, with a reduction of 24.92% from the initial design. The last step fabricates the optimized lattice structure through AM.  The stress distribution of the original suspension control arm part is simulated in ANSYS through static structure analysis. Figure 18 illustrates the relationship between the analysis results and the optimization results. As shown in Fig. 18a, highdensity lattices are observed at the positions with large equivalent stresses. A large number of x-direction tensile-stiff elementary lattices distribute inside the upper left beam of the control arm structure, which matches the x-direction normal stress distribution in the static simulation result as shown in Fig. 18b. Distributions of other composing lattices also match with the related stress component distributions from the static simulation, referring to Fig. 18c and d.
Another part design example of the quadcopter arm is shown in Fig. 19. The original design is shown in Fig. 19a and the  boundary condition for its working state is demonstrated in Fig. 19b. The optimized lattice structure with a constantthickness solid shell is shown in Fig. 19c, which measures 62.3% of the original design. Then, the DLP-based manufacturing re-sult is shown in Fig. 19e. The optimized and manufactured models from an alternative perspective can be seen in Fig. 19d and f, respectively, wherein a quarter of the model is cut away to show the internal lattice distribution.

Extension to Incorporate Macroscopic Topological Changes
This section explores further structural performance improvement by incorporating macroscopic topological changes with the proposed LSTO method. The MBB case is concentrated and a finer mesh is applied wherein the design domain is discretized into 121 500 cube elements. The boundary condition and meshing results are shown in Fig. 20. In addition, the maximum allowable volume fraction is 50%. First, the proposed LSTO method is applied and the optimization result of equivalent density distribution is demonstrated in Fig. 21. The optimized structural compliance is 65 mJ. Obviously, the structural performance can be further enhanced if the elements reaching the lower density limit are removed, since they contribute little to the overall structural stiffness while it is more effective to allocate the materials to the more stressed areas. Thus, one more density variable ρ a (0 ≤ ρ a ≤ 1) is applied to rebuild the constitutive model interpolation, as presented in equation (21): Here, C 0 = 1 and C min = 10 −9 are defined to prevent the stiffness matrix singularity. p ( p = 3) is the penalization factor to ensure the clear-cut lattice-void solution (Sigmund, 2001). With the above modification, ρ a = 0 removes the inefficient lattice elements while ρ a = 1 defines the lattice domain for variabledensity lattice infill. The optimized structural compliance is improved to 56 mJ with a 13.85% reduction from the original LSTO result. Distributions of the total lattice density, the topological variable, and the composing lattice densities are shown in Fig. 22, and a detailed illustration of the optimized lattice structure is shown in Fig. 23.
Furthermore, considering that the lattice density cannot reach the full solid state, the elements reaching the upper density limit can perform better if switching to solids is enabled. Hence, an alternative modification is made to the LSTO method to enable design of hybrid solid-lattice-void structures. Specifications of the modified design method are presented as follows. In the first stage, the SIMP-based topology optimization (Liu & Tovar, 2014) (macroscale topology optimization) subject to a 50% volume fraction constraint is performed, wherein the stiffness penalization parameter is set to p = 1. Hence, an optimized density field with a vast number of intermediate densities is derived as shown in Fig. 24a. Then, the design domain is divided into three subdomains according to the density distribution, where the areas with larger than 0.95 element densities are defined as the solid subdomain, the areas with smaller than 0.1 element densities are defined as the void subdomain, and the rest is defined as the lattice subdomain; refer to Fig. 24b. In the second stage, the solid and void subdomains remain unchanged while only the lattice subdomain is further optimized with the LSTO formulation of equation (14). Consequently, the final design is a mix of solids, lattices, and voids; i.e. the highly stressed elements can reach the solid state, the lowly stressed elements can reach the void state, and the rests are made of lattices accessing extreme anisotropic properties. Lattice structure optimization, conventionally subject to the manufacturability restriction, only employs intermediate lattice densities, e.g. ranging from 0.18 to 0.93, which, however, cannot access the fully solid or void state of the elements. Hence, classifying the design domain into solid, lattice, and void subdomains addresses the above issue by expanding the design scenario from lattice-only to solid-lattice-void hybrid structures. As shown in Fig. 25, the solid-lattice-void structure realizes a further  reduced compliance of 53 mJ. Finally, a comparison is made on the proposed method and the optimal solid-material structure using the same amount of materials. The comparative numerical results are shown in Table 5. It can be seen that the structural stiffness of the optimal solid-material structure is slightly better than the proposed lattice infill design. However, the solid materials are not porous and they do not possess the multifunctional properties of lattice materials. It is also possible that the solidmaterial structure is weaker in load-bearing capability because of the lower buckling strength of the thin structure members.

Discussion and Conclusion
This paper presents an LSTO method to design the 3D hierarchical structures. This proposed LSTO method features in the development of a novel candidate lattice for optimization through a combination of topologically optimized elementary lattices. The elementary lattices are optimized to be either x-direction tensile-stiff, y-direction tensile-stiff, or shear-stiff within the x-y plane. Therefore, the developed candidate lattice can have excellent x-direction tensile modulus, y-direction tensile modulus, x-y plane shear modulus, or even concurrently all of the above by flexibly arranging the relative density ratios of the Finally, a two-stage concurrent structural topology and lattice optimization method is developed that successfully reduces the structural compliance from 65 mJ of the lattice-only optimization result to 53 mJ of the hybrid solid-lattice-void structure, leading to a 18.46% reduction by incorporating the topological changes into the MBB example.  In particular, this paper addressed 3D structure optimization cases under 2.5D loading conditions, while the current candidate lattices may not perform equally well for more complex 3D loading conditions. We would admit this is-sue as a limitation of the current method, while further study will be performed in our forthcoming research to include more types of elementary lattices to address the real 3D loadings.  On the other hand, the proposed lattice structure may include several internal enclosed cavities when the density ratios of the composing lattices increase beyond certain thresholds. These interior enclosed cavities increase the difficulty of manufacturing since the unsolidified liquid resins attach to the cavity surface that cannot be washed away during post-processing. Therefore, improvement of the proposed LSTO method to concurrently design the liquid resin elimination channels deserves careful investigations.