Topology optimization design of three-dimensional multi-material and multi-body structure based on irregular cellular hybrid cellular automata method

In recent years, Hybrid Cellular Automata Method (HCAM) has been successfully applied to solve structural topology optimization problems. However, there was no report on HCAM research of three-dimensional composite structure composed of multiple materials and multiple bodies. Therefore, in this paper, three-dimensional non-cube cells of irregular size (such as tetrahedral cells with adaptive changes inside length) and Finite Element Method (FEM) are introduced to extend HCAM, which is better and more flexibly to fit complex geometric shapes. Furthermore, a better structure configuration of multi-material and multi-body structure is obtained. The typical example study showed that the proposed topology optimization method could effectively remove the redundant materials of multi-material and multi-body structure, and the optimized structure configuration could still meet the requirements of the original condition after geometric reconstructed. Thus it provided a reference for the intelligent design of other products.

In recent years, HCAM has attracted a large number of scholars. It is simple to understand, has good convergence, high computational efficiency and reduces numerical instability. It can be used to solve multi-objective problems and provides a new analytical idea and technical means for the effective topology optimization design of the structure 10,11 . Tovar et al. 12 have applied HCAM in the optimization design of vehicle collision mitigation, energy collector shell and blast mitigating armor structures based on a novel method for designing crashworthy structures with controlled energy absorption and harmonic vibration loading method 13 . Aiming at the problem of topology optimization of nonlinear structures with large deformation, Guo et al. 14 proposed a HCAM based multi-domain continuum structure topology optimization method, and then introduced the strain-based and multi-domain topology optimization algorithm into HCAM to study vehicle collision avoidance under large deformation 15 . In order to solve the static load requirements under normal operating conditions and energy absorption requirements targeting passive safety in crash events existing in the process of vehicle design, Aulig et al. 16 applied the proportional energy weighting method based on priority approach to HCAM, and applied it to the optimization design of the conceptual structure of industrial vehicle body. Afrousheh et al. 17 improved the search efficiency of HCAM by introducing the variable neighborhood radius concept, so as to achieve the purpose of improving the energy absorption of vehicle structure under high collision. This method realized an intelligent search strategy of improved HCAM. Da et al. 18 used HCAM to carry out topology optimization research on materials with extreme properties, and obtained the optimal topology configuration of the design object without obtaining sensitivity information through energy homogenization. Deng et al. 19 applied the HCAM in the thermal topology optimization design of the spindle structure under temperature-structural field coupling condition, after topology optimization, the spindle configuration not only reduced the material, but also improved its thermal characteristics. However, this method is only used for solving multi-field coupling topology optimization of 2D plane structures, and has not been carried out for the topology optimization of 3D structures. Jia et al. 20 proposed a topology optimization method for two-scale non-uniform microstructures based on HCAM, and applied the method to the optimization design of 2D and three-dimensional two-scale non-uniform microstructures.
In general, the research on 3D multi-material and multi-structure topology optimization using HCAM has not been reported. This paper proposes a 3D topology optimization method based on irregular cells of HCAM to study the 3D continuum structure with a variety of materials, which was aim to obtain the optimal structure configuration. Finally, the effectiveness of the method is verified by comparing with other algorithms. This paper discusses the new idea and gives the numerical calculation process based on the method to provide a reference for product intelligent design.

3D irregular cell HCAM model
3D HCAM topology optimization model based on irregular cell. Now 3D HCAM is commonly based on regular cube cells, which is shown in Fig. 1. Figure 1a is a 3D von Neumann type neighbor cell ( N = 6 ). Only six immediate neighbors are considered, and these neighboring cells have common surfaces with the central cell. Figure 1b shows a 3D Radial type neighbor cell ( N = 18 ). The neighbor cell with common edge or surface with the central cell is considered to be the neighbor cell of the central cell. Figure 1c is a 3D Moore type neighbor cell ( N = 26 ), and any cell with a common surface or edge or vertex with the central cell can be considered as the neighbor cell of the central cell.
In many cases, cube cells cannot cover the design area uniformly. In practical engineering analysis and design, irregular cells can be used to deal with complex geometric figures. In key areas of concern, cells need to be locally refinement. In this paper, the traditional cube cells of regular size are expanded to three-dimensional non-cube cells of irregular size (such as tetrahedral cells with adaptive side length), which could fit complex geometric shapes better and more flexibly. Thus, better solutions were obtained when solving topological configurations, as shown in Fig. 2. In some places, the shape of non-cube cells appears, even though the cube cells is set, as shown in Fig. 2a. At the same time, introducing the idea of local grid refinement in complex model fitting, a sharp change in the stress concentration or strain area for local cell refinement, which do not need to use the whole structure of fine cell. And it can be used to avoid excessive number of cellular automata and cause cell with the difficulty of the mapping between Finite Element mesh, speed up the calculation, increase the authenticity of the results. As shown in Fig. 3, only four adjacent neighboring cells for the tetrahedral von Neumann type neighboring cells (the number of neighboring cells is 4) are considered. And these adjacent cells had common surfaces with the  where, J is the number of cellular states; a J i (t) is the Jth state defined for discrete time t and discrete position i. The state at a time (t+1) under each cell is jointly determined by its state at time t and the state of its neighboring cell. The Jth cell's status update rule can be expressed as where, a i+�1 (t),…,a i+N (t) is the neighborhood cell of cellular i; N is the number of cells in the neighborhood of a cell; R J i is the local update rules of cell i, which is the same for all positions. The neighbor cell query rules of 3D tetrahedral irregular cell HCAM are as follows: all the cells and their nodes are taken out at first, and then the nodes of each cell are compared. If two cells have three identical nodes, the two cells are defined as neighbors.
In the iterative process of topology optimization using HCAM, the design variable x i (t) of each cell can be defined as cell density, geometric size, elastic modulus, etc., and the state field variable S i (t) of each cell can be defined as stress, strain, strain energy density or a function of them, Then, the state a i (t) of the cell at time t can be expressed as The design variable x i (t) is determined according to the material model. In this paper, the relative density of the element is taken as the design variable, and the definition of the field variable S i (t) depends on the optimization problem to be solved. Then, in a discrete region, the strain energy E of the whole structure is  where, ε x , ε y , ε z , γ xy , γ yz , γ xz are the x-direction strain, y-direction strain, z-direction strain and shearing strain of the node respectively; σ x1 , σ y1 , σ z1 are the normal stresses in x-direction, y-direction and z-direction respectively; τ xy , τ yz , τ xz are the shear stress on surfaces xy, yz and xz. The strain energy density is the local rigidity index of the structure. Thus, the strain energy density is used as the field variable in this paper, and the optimization criterion refers to the uniform distribution of the strain energy density. From this standard, the CA state can be defined as The topology optimization design of the structure is carried out by using HCAM, and the definition of strain energy density U i (t) can be determined by the specific optimization objective function. The objective of optimization is to minimize the difference between the mean value of cell strain energy density U i and the set value of state field variable, the optimization problem can be expressed by Eq. (7) where, e i is the error signal; U * is the set value of strain energy density. In the equation, there is a requirement of lower limit x min , so that the singular matrix can be avoided in the process of FEA.
At this point, the state of the cell, the strain energy of the whole structure and the strain energy density U i (t) are the same as the regular cell model, but the U i (t) of the strain energy density of cell i is changed, which is expressed by Eq. (8): where, U j is the strain energy density of adjacent cell j which is adjacent to cell i; S is the average area of each cell in CA; S i is the area of cell i; S j is the area of neighboring cell j which is adjacent to cell i.
In order to expand the application range of HCAM optimization model and accelerate the convergence of topology optimization design results, so as to adapt to structural topology optimization with fixed and variable loads, the target value of strain energy density U * in the process of topology optimization is calculated by using Eq. (9) where, aif is weight coefficient; V 0 is the volume of the original structure; E t is the total strain energy of the structure when iteration to t times.

Convergence criterion.
In the iterative process of topology optimization, whether the structure meets the criteria of convergence is determined by the type of design rules used in the constantly updated design variables. The change of structure mass, �M(t) is taken as the convergence criterion, then the structure mass is expressed as where, M(t) is the total mass of the structure from iteration to t times; ρ i (t) is the density of cell i when iteration to t time; ρ is solid material density. When the cell relative density x i (t) is 1, the cell is composed of solid material; When the cell relative density x i (t) is 0, the cell can be considered to have no material.
When �M(t) in mass goes to zero, that means there's no change in mass. Then it is determined that the optimization process has reached the convergence standard. At the same time, in order to avoid the partial phenomenon that small changes in quality may occur in the application process of HCAM and cause large changes in structure, and to reduce the accidental uncertainty, the convergence criterion uses the average change of two consecutive iterations to determine where, M 0 is the total mass of the initial structure. Mass ratio F is proposed to reflect the specific optimization degree of the current iteration number structure, then the mass ratio F(t) at time t can be expressed as follows Local control rule. The material distribution in the design domain is determined by the local control rule, which seeks to minimize the difference e i (t) between the state field variable U i (t) and the set value of the state field variable U * when iteration to t times, so that the calculation result tends to the target value. Linear control rule is adopted in this paper to express the change value of the relative cell density of the material as follows where, f (e i (t)) is a local control rule, C p is a linear weight coefficient, usually a normal number.

Topology optimization design of 3D monomer structure
In order to verify the effectiveness and universality of HCAM based on tetrahedral irregular cell topology optimization algorithm, topology optimization design is carried out by using different three-dimensional monomer structures and working cases.
Topology optimization design of L-shaped bracket. The L-shaped bracket structure is shown in  Fig. 4. A fixed constraint is imposed on the Surface A of the L-shaped bracket, and the concentrated load P = 6.5 N along the negative direction of the Z-axis is set on the right edge of the transverse arm. In this case study, tetrahedral elements are used as the cell shape of HCAM, and tetrahedral von Neumann type is selected as the neighboring cells, when calculating by FEM, selecting the mesh element SOLID186, a total of 66,005 tetrahedral cells (grids) are generated by dividing the structure, and the cell size is consistent with the grid size. The optimization goal is to remove 30% of the material, The initial design variable x i (0) = 1 , the convergence error ε = 10 −4 , the weight coefficient aif = 1.75 , the local control rules are linear control rules, the linear weight coefficient C p = 1.
Then, the number of iterations to solve the result is 38 based on HCAM calculation results, and the final topology optimization configuration obtained is shown in Fig. 5a. It is very similar to the topology optimization configuration obtained by ANSYS Topological Optimization Module in Fig. 5b.
The topology optimization configuration obtained above is output to CAD software for geometric reconstructed of the model. The reconstructed model basically kept the original topology structure of the optimization result, and the results are shown in Fig. 6.
In order to verify whether the geometrically reconstructed model could meet the requirements of working conditions, static simulation calculation is carried out under the same boundary conditions. At the same time, the stress and displacement comparison results are obtained as shown in Figs. 7 and 8. As can be seen from the stress contour figures in Fig. 7, the maximum stress position of the original model and the reconstructed model Figure 4. L-shaped bracket structure model.  Fig. 8, the maximum displacement of the original model appeared below the right-most end of the L-shaped bracket cross arm, and the maximum displacement of the reconstructed model appears at the upper right end of the L-shaped bracket cross arm, which are 1.29 × 10 −2 mm and 1.92 × 10 −2 mm respectively.
For the examples of L-shaped support structure under the same working conditions, the mass ratio F is set to 30%, 50% and 70% respectively, and the results are compared with those obtained by FEM. The results are shown in Table 1.
The performance indexes of each generation are calculated, and the fitting curve of HCAM performance indicators is obtained, as shown in Fig. 9. It could be seen that the topology optimization results of the three mass ratios tend to converge. As reported in the literature 11 , hundreds of thousands of iterations for same conventional topology optimization methods are required to achieve convergence. When using the HCAM presented in this work, the number of iterations can be reduced to about 60 iterations.
By comparing the optimization results of HCAM and FEM under different mass ratios, it could be seen that the topology optimization configuration of 3D structures with different mass ratios could be obtained effectively using HCAM. And it is very similar to that obtained by FEM.
Topology optimization design of the regular triangular yoke plate structure. The regular triangular yoke plate structure is shown in Fig. 10, where its side length l = 60 mm, the thickness h = 2 mm, fillet radius R = 10 mm, round hole radius r = 5 mm. The modulus of elasticity of the material E 0 = 200 GPa, Poisson's ratio μ = 0.3, density ρ = 7850 kg/m 3 , and yield limit σ s = 250 MPa . The constraints and load settings are   Fig. 10. Fixed constraints are imposed on the Surface A and Surface B of the circular hole Surface. And concentrated load P = 800 N in the negative direction along the Z axis is set on the Surface C of the circular hole Surface. Meanwhile, tetrahedral cell are used as the cell shape of HCAM, and tetrahedral von Neumann type is selected as the neighboring cells. During the Finite Element calculation, the mesh element is chosen as SOLID186 (high-order 3D 20-node solid structure element), and a total of 108,297 tetrahedral grids (cells) are generated after dividing the structure, and the cell size is consistent with the grid size so that map to each other. The optimization goal is to remove 50% of the material. The initial design variable x i (0) = 1 , the convergence error ε= 10 −4 , the weight coefficient aif = 0.5 , the local control rule adopts the linear control rule, and the linear weight coefficient C p = 0.5. . After HCAM calculation, the number of iterations is 26, and the final topology optimization configuration obtained is shown in Fig. 11a, which is very similar to the topology optimization configuration obtained by using the ANSYS topology optimization module in Fig. 11b.
The topology optimization configuration obtained above is output to CAD software for geometric reconstructed of the model. The reconstructed model basically kept the original topology structure of the optimization result, and the results are shown in Fig. 12.  12 MPa, and that of the reconstructed model is about 119.01 MPa, when 50% of the material is removed, the maximum stress increased by about 46.7%, but the yield limit of the material is 250 MPa, and the safety factor n is 2, which could still meet the strength requirements. Meanwhile, as can be seen from the displacement contour figures in Fig. 14, the maximum displacement position of the original model and the reconstructed model is basically the same, both appearing on the rounded corner surface outside the circular hole Surface C. The maximum displacement of the initial scheme is 9.69 × 10 −3 mm , and the maximum displacement of the reconstructed structure is 12.40 × 10 −3 mm.
The material removal ratio F of the triangular yoke plate structure is set to 30%, 70% and 80% of the original structure respectively, and HCAM and FEM are used to carry out the comparative analysis of its topology optimization. The results are shown in Table 2.
The performance indexes of each generation are calculated, and the fitting curve of HCAM performance indicators is obtained, as shown in Fig. 15. It could be seen that the topology optimization results of the three mass ratios tend to converge less than 25 iterations. By comparing the optimization results, it could be seen that the topology optimization configuration are very similar obtained by HCAM and FEM.

3D multi-material and multi-body HCAM model
Finite element analysis model for contact problem. In 3D multi-material and multi-body analysis, the contact relationship between objects is unavoidable, and FEM is the most effective method to deal with the contact problem. As shown in Fig. 16a, object A and object B are in contact. Assuming that S a is the contact surface and S b is the target surface (as shown in Fig. 16b), the boundary on object B is fixed, external load P is applied on object A, and the two points on which surface S a and S b contact each other are called contact point pairs. According to literature 21      www.nature.com/scientificreports/ In the analysis of FEM in this paper, bonded type is adopted for the contact model of multi-material and multi-body structure, that is, relative sliding or separation between contact surfaces is not allowed. The contact element uses the 4-node CONTACT174 element, and the corresponding target element uses the 8-node TAR-GET170 element. A single "contact pair" consists of a target element and a contact element. The program uses a shared real constant to identify a "contact pair". When establishing a "contact pair", the target element and the contact element must specify the same real constant 22 .

Convergence criterion of multi-material and multi-body structure. Whether a multi-material
and multi-body structure achieves convergence criteria is also determined by the type of design rules used by the constantly updated design variables. Assuming that the multi-material and multi-body structure contains N objects, the convergence criterion can be taken as the change �M Tot (t) of the total mass M Tot (t) , which can also be determined according to the average change of two consecutive iterations where, Q = 1, 2,…, N.

Topology optimization of three-dimensional multi-material and multi-body structures
Topology optimization design of cylinder-torus plate structure. The cylinder-torus plate assembly structure as shown in Fig. 17 is a multi-material and multi-body structure, which was consisted of a cylinder plate A with a radius r = 3 mm and an torus plate B with an outer diameter R = 10 mm. The outer perimeter of A fixed constraint is applied on the Surface A for the cylinder plate A. And a concentrated load P = 30 N is applied on the outer circumferential Surface I-IV of torus plate B. In this case study, tetrahedral elements are used as the cell shape of HCAM, and tetrahedral von Neumann type is selected as the neighboring cells. The cylindrical plate A is regarded as the rigid body, and the torus plate B as the flexible body. When calculating with FEM, select the mesh element as SOLID186 and contact type as Bonded. The 3D 8-node Target170 element is used at the target surface, and the 3D 4-node Contect174 element is used at the contact surface, and constrain the contact surface between them to produce static friction only. There are 107,330 tetrahedral cells (grids) after dividing the cylinder-torus plate structure. The torus plate B is divided into 96,923 tetrahedral cells (grids), and the cell size is consistent with the Finite Element size. The torus plate B is set as the design domain, and the optimization objective is set to remove 60% of the material. The initial design variable x i (0) = 1 , the convergence error ε = 10 −3 , the weight coefficient aif = 65,000, the local control rule is the linear control rule, the linear weight coefficient C p = 0.45.   Fig. 18a, which is very similar to the result obtained by ANSYS in Fig. 18b. The effectiveness and reliability of the proposed HCAM based tetrahedral irregular cell topology optimization algorithm for multi-material and multi-body are verified. www.nature.com/scientificreports/ In this paper, the obtained multi-material and multi-body topology optimization configuration is output to CAD software for geometric reconstructed of the model. The reconstructed model basically kept the original topology structure of the optimization result, as shown in Fig. 19.
In order to verify whether the geometrically reconstructed model can meet the requirements of working conditions, static simulation calculation is carried out under the same boundary conditions. Thus the stress and displacement comparison results are obtained, as shown in Figs. 20 and 21. As can be seen from the stress contour figures in Fig. 20, the maximum stress of both the original model and the reconstructed model appears at the contact position of the two plates. The maximum stress of the original model is about 89.72 MPa and that of the reconstructed model is about 94.11 MPa. When 60% materials are removed, the maximum stress increases by 4.9%. But the material yield limit of the torus plate B is 250 MPa, and the safety factor n is 2, which can still meet the strength requirements. At the same time, the maximum stress of cylinder plate A is 47.49 MPa, which is far lower than the material yield limit of cylinder plate A of 280 MPa. Figure 21 is the displacement contour figure of the structure before and after optimization, it can be seen from the figure that the maximum displacement positions of the original model and the reconstructed model are basically the same, both appearing near the place (I-IV) where the loads are applied, which are 1.32 × 10 −3 mm and 2.04 × 10 −3 mm respectively.
Topology optimization design of support-connecting rod-support structure. As shown in   www.nature.com/scientificreports/ In this case study, tetrahedral elements are used as the cell shape of HCAM, and tetrahedral von Neumann type is selected as the neighboring cells. The support A and B are taken as rigid body, connecting rod C is taken as a flexible body. The mesh element SOLID186 and contact type bonded are selected to calculate for FEM. The 8-node TARGET170 element is used at the target surface, and the 4-node Contect174 element is used at the contact surface. And constrain the contact surfaces between them to produce only static friction. A total of 116,111 tetrahedral cells are obtained by dividing the support-connecting rod-support boy structure into irregular cells. Among them, connecting rod C is divided into 80,342 tetrahedral cells (grids), and the cell size is consistent with the Finite Element size. Connecting rod C is set as the design domain and the optimization objective is set to remove 40% of the material. The initial design variable x i (0) = 1 , the convergence error ε = 10 −4 , the weight coefficient aif = 2.75 , the local control rule is the linear control rule, the linear weight coefficient C p = 0.225.
After HCAM calculation, the number of iterations to finally solve the result is 52, and the final topology optimization configuration obtained is shown in Fig. 23a, which is similar to the result obtained by ANSYS in Fig. 23b.
The obtained multi-material and multi-body topology optimization configuration is output to CAD software for geometric reconstructed of the model. The reconstructed model basically kept the original topology structure of the optimization result, as shown in Fig. 24.
In order to verify whether the geometrically reconstructed model can meet the requirements of working conditions, static simulation calculation is carried out under the same boundary conditions, and the stress and displacement comparison results are obtained as shown in Figs. 25 and 26. As can be seen from the stress contour figures in Fig. 25, the maximum stress in both the original model and the reconstructed model occurs below the contact position between connecting rod C and support A. The maximum stress of the original model is 79.95 MPa, while the maximum stress of the reconstructed model is 98.97 MPa. When 40% material is removed, the maximum stress increases by 23.80%, but the yield limit of connecting rod C material is 250 MPa, and the safety factor n is 2, which can still meet the strength requirements. At the same time, the maximum stress value of the non-design domain support A and B is 32.89 MPa, which is far lower than the material yield limit of www.nature.com/scientificreports/ the two supports 280 MPa, so the support-connecting rod-support body structure meets the strength requirements on the whole. Figure 26 is the displacement contour figure of the structure before and after optimization. It can be seen from the figure that the maximum displacement positions of before and after optimization are basically consistent, which all appear in the middle part of connecting rod C, which are 3.42 × 10 −3 mm and 5.21 × 10 −3 mm , respectively.

Concluding remarks
In this paper, a topology optimization design method for 3D multi-material and multi-body structures was proposed. The 3D tetrahedral cells with adaptive side length change are used to replace the traditional regular cube cells of HCAM, which can more flexibly adapt to complex geometric shapes. This avoids the problem that the cubic cells cannot uniformly cover the design area, so that more flexible solutions can be obtained when solving the topology configuration. At the same time, in order to solve the problem of contact relationship between objects in the multi-body structure, the contact problem processing method of FEM is applied to HCAM topology optimization, and the structure optimization configuration which is close to the existing commercial software module is obtained. Through the verification of typical numerical examples, the 3D multi-material and multi-body topology optimization method based on HCA can effectively remove the redundant materials