Substructure-Based Topology Optimization for Symmetric Hierarchical Lattice Structures

Abstract: This work presents a topology optimization method for symmetric hierarchical lattice structures with substructuring. In this method, we define two types of symmetric lattice substructures, each of which contains many finite elements. By controlling the materials distribution of these elements, the configuration of substructure can be changed. And then each substructure is condensed into a super-element. A surrogate model based on a series of super-elements can be built using the cubic B-spline interpolation. Here, the relative density of substructure is set as the design variable. The optimality criteria method is used for the updating of design variables on two scales. In the process of topology optimization, the symmetry of microstructure is determined by self-defined microstructure configuration, while the symmetry of macro structure is determined by boundary conditions. In this proposed method, because of the educing number of degree of freedoms on macrostructure, the proposed method has high efficiency in optimization. Numerical examples show that both the size and the number of substructures have essential influences on macro structure, indicating the effectiveness of the presented method.


Introduction
Symmetrical structure is a kind of widely used structure. In topology optimization [1,2], the structure is designed according to symmetry, so as to realize the multi-functional design of macro structure [3][4][5]. The symmetrical structure is more suitable for manufacturing. It is helpful to avoid the problem of low manufacturability caused by the homogenization method [6].
Because of its scale separation assumption [7,8], the material is often unconnected in the multiscale structure [9,10], and the designed structures are not manufacturable in practices using this method. Moreover, due to the large overhang hand angle and length [11][12][13], the optimized configuration only exists in theory. To overcome the barrier between topology optimization and manufacturing, many scholars try to find the topology optimization method of integrated manufacturing. Osanov et al. [14] discussed the manufacturability of multiscale designs and presented the fundamental findings for designing the periodic structure. Wang et al. [15] presented an optimization method for the cellular lattice structure based on the shape metamorphosis, and this optimized structure can be fabricated in additive manufacturing. Based on the thermal conductivity of lattice structure, Cheng et al. [4] proposed a concurrent optimization method based on the boundary conditions to optimize the lattice infill and design-dependent movable features. To verify the accuracy and efficiency of the optimization method, Wang et al. [16] presented a periodic lattice optimization method using Isogeometric analysis. Based on the symmetry of microstructure, most of works tackled the design method to optimize the two scales structure. In order to ensure the material connectivity among microstructures, the shape geometries or physical responses interpolation method has been adopted to design the macro/micro structures [17,18]. To interpolate a prototype microstructure and obtain a family of graded connectable microstructures, Wang et al. [19] proposed a shape metamorphosis technology to design the multiscale structure. Based on the homogenization method, Groen and Sigmund [20] presented an advanced projection scheme to achieve high-resolution and connectable microstructures. To validate the accuracy of the asymptotic homogenization method, Arabnejad et al. [21] determined the effective elastic moduli and yield strength of six lattice structures to optimize the multiscale structure. In order to further extend the design space, Li et al. [22] presented a special composite structure with microstructures of identical topology with different length/width ratios. Considering the adjacent cellular materials, Radman et al. [23] proposed a progressive design strategy to design the cellular materials using the globally filtered sensitivities.
This paper presents a topology optimization design method for symmetric structures based on substructure. According to the static condensation [24,25] or domain decomposition [26], the substructure-based method is suitable for the large-scale computing of structural topology optimization. In the calculation of the condensed stiffness matrix, the calculation efficiency can be doubled due to the symmetry of this stiffness matrix [27,28]. In this paper, we also design two types of microstructures with symmetry. It should be noted that the symmetry of the condensed stiffness matrix is not related to the symmetry of microstructure.
The structure and content of this paper is organized as follows. The static condensation of the substructure is described in detail in Section 2, the construction method for the substructure based symmetric microstructure is elaborated in Section 3, and the topology optimization model and the sensitivity analysis for this proposed topology optimization method is expressed in Section 4. To demonstrate the capability of the proposed method, three numerical examples are presented in Section 5. Finally, a summary of some open problems related to this paper is discussed in Section 6.

The Static Condensation of Substructure
In the finite element method, the finite mesh in any region can be condensed into a super-element shown in Figure 1. And the strain and stress fields of its internal nodes can be represented by their boundary nodes [27,28]. Based on the symmetry of microstructure, most of works tackled the design method to optimize the two scales structure. In order to ensure the material connectivity among microstructures, the shape geometries or physical responses interpolation method has been adopted to design the macro/micro structures [17,18]. To interpolate a prototype microstructure and obtain a family of graded connectable microstructures, Wang et al. [19] proposed a shape metamorphosis technology to design the multiscale structure. Based on the homogenization method, Groen and Sigmund [20] presented an advanced projection scheme to achieve high-resolution and connectable microstructures. To validate the accuracy of the asymptotic homogenization method, Arabnejad et al. [21] determined the effective elastic moduli and yield strength of six lattice structures to optimize the multiscale structure. In order to further extend the design space, Li et al. [22] presented a special composite structure with microstructures of identical topology with different length/width ratios. Considering the adjacent cellular materials, Radman et al. [23] proposed a progressive design strategy to design the cellular materials using the globally filtered sensitivities.
This paper presents a topology optimization design method for symmetric structures based on substructure. According to the static condensation [24,25] or domain decomposition [26], the substructure-based method is suitable for the large-scale computing of structural topology optimization. In the calculation of the condensed stiffness matrix, the calculation efficiency can be doubled due to the symmetry of this stiffness matrix [27,28]. In this paper, we also design two types of microstructures with symmetry. It should be noted that the symmetry of the conden ed stiffness matrix is not related to the symmetry of microstructure.
The structure and content of this paper is organized as follows. The static condensation of the substructure is described in detail in Section 2 he construction method for the substructure based symmetric microstructure is elaborated in Section 3, and the topology optimization model and the sensitivity analysis for this proposed topology optimization method is expressed in Section 4. To demonstrate the capability of the proposed method, three numerical examples are presented in Section 5. Finally, a summary of some open problems related to this paper is discussed in Section 6.

he static condensation of substructure
In the finite element method, the finite mesh in any region can be condensed into a superelement shown in Figure 1. And the strain and stress fields of its internal nodes can be represented by their boundary nodes [27,28]. Here, based on the internal and boundary modes, the equilibrium equation of a substructure can be established as follows: where and are the stiffness matrix of the boundary and internal nodes, respectively.
are the coupling sub-stiffness matrixes. and are the response sub-vector of the nodes.
Here, applying the Guyan method [24] the displacement vector can be expressed as follows: Here, based on the internal and boundary modes, the equilibrium equation of a substructure can be established as follows: where K mm and K ss are the stiffness matrix of the boundary and internal nodes, respectively. K sm , K ms are the coupling sub-stiffness matrixes. U m and U s are the response sub-vector of the nodes.
Here, applying the Guyan method [24] the displacement vector can be expressed as follows: where U s = TU m . Substituting Equation (2) into Equation (1), the first row of Equation (1) yields the following: Substituting Equation (5) into Equations (2) and (3), the condensed stiffness matrix can be obtained, Because only boundary nodes are reserved, the size of the condensed matrix is much smaller than the original matrix.

Substructure-Based Symmetric Microstructure Model
In order to construct the stiffness with different densities, a series of condensed substructures with a similar configuration for sampling densities using the finite element method shown in Figure 2 are considered. The density of the associated substructure is controlled by the parameter t ∈ R. Here, the self-defined substructure only contains two types of material: full of materials or void materials. In the process of the proposed optimization method, the more elements of the substructure, the larger the stiffness matrix of the corresponding surrogate model, and the longer the optimization iteration time. Considering the efficiency of iterative calculation, we fix the number of elements of substructure here [29,30]. Considering the computation efficiency, all substructures are discretized into 100 × 100 square bilinear elements.
Because only boundary nodes are reserved, the size of the condensed matrix is much smaller than the original matrix.

Substructure-based symmetric microstructure model
In order to construct the stiffness with different densities, a series of condensed substructures with a similar configuration for sampling densities using the finite element method shown in Figure   2 are considered. The density of the associated substructure is controlled by the parameter t ∈ .
Here, the self-defined substructure only contains two types of material: full of materials or void materials. In the process of the proposed optimization method, the more elements of the substructure, the larger the stiffness matrix of the corresponding surrogate model, and the longer the optimization iteration time. Considering the efficiency of iterative calculation, we fix the number of elements of substructure here [29,30]. Considering the computation efficiency, all substructures are discretized into 100×100 square bilinear elements. Here, each substructure has the same boundary nodes. According to Eq. ( ), we consider the constructed a series of condensed stiffness matrices for sampling densities For the convenience of calculation, each condensed stiffness matrix is restored in a column vector of length attributed. Here, we assume that a condensed stiffness matrix at any density can be expressed as follows: where coefficients i i α ρ are scale functions of density .
i ρ can be expressed in matrix form: Here, each substructure has the same boundary nodes. According to Equation (4), we consider the constructed a series of condensed stiffness matrices K * for R sampling densities k * For the convenience of calculation, each condensed stiffness matrix is restored in a column vector of length N 2 attributed. Here, we assume that a condensed stiffness matrix at any density can be expressed as follows: where coefficients α i (ρ i ) are scale functions of density ρ. K * (ρ i ) can be expressed in matrix form: According to symmetry, each element in a matrix can be written as follows: The coefficient can be obtained using use cubic spline interpolation method: where . The condensed stiffness matrix K * (ρ i ) can be reconstructed explicitly, Here, we consider a generalized model for substructure stiffness augmented with a penalized density: In the proposed model for minimizing the compliance of structure, the penalty factor p≥1 means that the associated matrix with intermediate densities is penalized.

The Topology Optimization Model
According to Equation (10), the global stiffness matrix can be re-expressed as follows:

Sensitivity Analysis
The Derivative of the Objective and Constraint with Respect to the Densities ρ i can be obtained, Based on Equation (12), ∂K * /∂ρ j can be obtained as follows: Here, to update the density design variables [30], we choose the optimality criteria (OC) as the optimality criteria: And the parameter D i can be found as follows: where the κ is Lagrange multiplier. ∂V/∂ρ i = v i is the volume of i-th substructure, which is an constant value in the SIMP method. ∂c/∂ρ i is the sensitivity number given in the Equation (13).
In this paper, we assume that the design result is the optimal structure when the difference objective function value between two iterations is less than 10 −3 .

Numerical Example
Here, we set the Young's modulus to E=1, and the Poisson's coefficient to 0.3 in our given examples.

Clamped Ends Beam
The design domain of the problem is a rectangle of size 5m × 1m that fixed at both sides, as shown in Figure 3. A load F=1 is placed at the center. The volume fraction of solid material to be preserved is set to 30%. symmetry 2020, x, x; doi: FOR PEER REVIEW www.mdpi.com/journal/symmetry And the parameter i D can be found as follows: where the κ is Lagrange multiplier.
i i V v ρ ∂ ∂ = is the volume of i-th substructure, which is an constant value in the SIMP method. i c ρ ∂ ∂ is the sensitivity number given in the Section 3.2.
In this paper, we assume that the design result is the optimal structure when the difference objective function value between two iterations is less than 10 -3 .

Numerical example
Here, we set the Young's modulus to E=1, and the Poisson's coefficient to 0.3 in our given examples.

Clamped ends beam
The design domain of the problem is a rectangle of size × that fixed at both sides, as shown in Figure 3. A load F=1 is placed at the center. The volume fraction of solid material to be preserved is set to 30%.  In this design example, three discretization resolutions are considered to divide the whole design domain: N x × N y = 20 × 4, N x × N y = 30 × 6 and N x × N y = 60 × 12. in the process of topology optimization, the filer radiuses are set as 1.1 times of the substructure length. To avoid the checkerboard pattern issue, we consider the optimal configuration under different penalty factors 1, 2, and 3, respectively. The optimized results for the three discretization resolutions are shown in Figures 4-9. From the optimal result, the following can be observed: (1) The optimization results are all symmetrical structures, independent of the discretization resolutions. (2) The symmetry of microstructure does not affect the symmetry of macrostructure. It can be observed that the compliance value of the design with case B lattice is lower than the design with Case A lattice. And this compliance value decrease with the increase in the number of substructures.

The Cantilever Design
The cantilever design problem is given in Figure 10. The ratio of the horizontal and vertical dimension of this design domain is 2. The left-end of the cantilever is clamped and a unit vertical load is applied on the middle point of the right edge. The volume fraction of solid material to be preserved is set as 0.4. symmetry 2020, 17, x FOR PEER REVIEW 8 of 12

The cantilever design
The cantilever design problem is given in Figure . The ratio of the horizontal and vertical dimension of this design domain is 2. The left-end of the cantilever is clamped and a unit vertical load is applied on the middle point of the right edge. The volume fraction of solid material to be preserved is set as 0.4. The penalty is set to p=1,2,3. The design results of the two resolutions for the two cases are shown in Figures 11 and 12. It can be concluded that the compliance value increases when the penalty increasing. For the same discretization resolution, the objective function value with case A is larger than that of case B. Here, the whole design domain is divided into N x × N y = 8 × 4 and N x × N y = 16 × 8, respectively. The penalty is set to p = 1, 2, 3. The design results of the two resolutions for the two cases are shown in  The penalty is set to p=1,2,3. The design results of the two resolutions for the two cases are shown in Figures 11 and 12. It can be concluded that the compliance value increases when the penalty increasing. For the same discretization resolution, the objective function value with case A is larger than that of case B.
x y N N ,p ,c .

5.3.D Simply supported cube design
In this example, a cube is clamped at four bottom corners. And the design domain is D . = × × , as shown in Figure 13. A load is applied to the center of the top surface. Here, we set the volume faction as 30%. To improve the calculation efficiency, the substructure is discretized into × × eight nodes cubic elements Figure 13. Definition of the 3D simply supported cube design problem.
To simplify, the given design domain is only discretized into substructures. Due to the large number of intermediate density substructures, it is unnecessary to obtain the clarify structural pattern design. We design the structure with the penalty p=2 and 3. Similar to the example, the filter radius is set as times of the substructure length. And the corresponding design results are summarized in Figures 14 and 15. Form these results, the compliance value with case B is more likely to converge to compact topology pattern, and this value

D Simply Supported Cube Design
In this example, a cube is clamped at four bottom corners. And the design domain is D = 1 m × 0.5 m × 1 m, as shown in Figure 13. A load is applied to the center of the top surface. Here, we set the volume faction as 30%. To improve the calculation efficiency, the substructure is discretized into 30 × 30 × 30 eight nodes cubic elements.

5.3.D Simply supported cube design
In this example, a cube is clamped at four bottom corners. And the design domain is D . = × × , as shown in Figure 13. A load is applied to the center of the top surface. Here, we set the volume faction as 30%. To improve the calculation efficiency, the substructure is discretized into × × eight nodes cubic elements Figure 13. Definition of the 3D simply supported cube design problem.
To simplify, the given design domain is only discretized into substructures. Due to the large number of intermediate density substructures, it is unnecessary to obtain the clarify structural pattern design. We design the structure with the penalty p=2 and 3. Similar to the example, the filter radius is set as times of the substructure length. And the corresponding design results are summarized in Figures 14 and 15. Form these results, the compliance value with case B is more likely to converge to compact topology pattern, and this value for the case A lattice substructure is a little smaller than case B with the same penalty factor. With an Figure 13. Definition of the 3D simply supported cube design problem.
To simplify, the given design domain is only discretized into N x × N y × N z = 40 × 20 × 40 substructures. Due to the large number of intermediate density substructures, it is unnecessary to obtain the clarify structural pattern design. We design the structure with the penalty p = 2 and 3.
Symmetry 2020, 12, 678 9 of 11 Similar to the example, the filter radius is set as times of the substructure length. And the corresponding design results are summarized in Figures 14 and 15. Form these results, the compliance value with case B is more likely to converge to compact topology pattern, and this value for the case A lattice substructure is a little smaller than case B with the same penalty factor. With an increasing penalty factor, the design becomes more compact with a clear structural pattern. Here, we can conclude that the optimization results are symmetric as long as the boundary constraints are symmetric. And the layout of microstructure also has symmetry, which is consistent with that of the macrostructure. Here, we can conclude that the optimization results are symmetric as long as the boundary constraints are symmetric. And the layout of microstructure also has symmetry, which is consistent with that of the macrostructure. Here, we can conclude that the optimization results are symmetric as long as the boundary constraints are symmetric. And the layout of microstructure also has symmetry, which is consistent with that of the macrostructure.

Conclusions
In this paper, we present a substructure-based topology optimization for symmetric lattice structures with substructuring to minimize compliance of the multi-scale structure. Form the design result, all optimization design are symmetrical structures, which is unrelated to the discretization resolutions. Moreover, the symmetry of microstructure does not affect the symmetry of macrostructure. In our proposed method, the reconstruct stiffness matrix only calculates the element in one upper or lower triangle area using the cubic spline interpolation based on the symmetry of the stiffness matrix. Three numerical examples show that this proposed method for the reconstructing matrix can be applied to two-and three-dimensional topology optimization problems.