Topology optimization for multiscale design of porous composites with multi-domain microstructures

This paper proposes a new multiscale topology optimization method for the design of porous composites composed of the multi-domain material microstructures considering three design elements: the topology of the macrostructure, the topologies of multiple material microstructures and their overall distribution in the macrostructure. The multiscale design involves two optimization stages: the free material distribution optimization and the concurrent topology optimization. Firstly, the variable thickness sheet (VTS) method with the regularization mechanism is used to generate multiple element density distributions in the macro design domain. Hence, different groups of elements with the identical densities can be uniformly arranged in their corresponding domains, and each domain in the space will be periodically configured by a unique representative microstructure. Secondly, with the discrete material distributions achieved in the macro domain, the topology of the macrostructure and topologies of multiple representative microstructures are concurrently optimized by a parametric level set method combined with the numerical homogenization method. Finally. Several 2D and 3D numerical examples are provided to demonstrate the effectiveness of the proposed multiscale topology optimization method.

Porous composites, comprising solids and voids, are artificially engineered to have the superior structural performances but lightweight, such as the higher specific stiffness and strength, better fatigue strength and improved corrosion-resistance [1,2] . Over the last decades, porous composites with periodically distributed microstructures have received great popularity in many engineering applications due to their easiness in manufacturing and cost effectiveness in mass production. As an example, the cellular honeycomb structures have been considerably employed in the automobile and aerospace industries [3,4].
Hence, an increasing number of analytical, numerical or experimental methods for the design of porous composites to achieve the advantage properties have been developed [1,2]. However, an efficient design method to systematically obtain the optimal layout of porous composites is still in demand.
Topology optimization [5,6] aims to seek for the best material layout in a given design domain subject to the loads and boundary conditions until the concerned performance is optimized. A wide range of methods have been developed in recent years, such as the homogenization method [7], the solid isotropic material with penalization (SIMP) [8,9], the evolutionary structural optimization (ESO) [10] and the level set method (LSM) [11][12][13], as well as the point wise-density interpolation (PDI) method [14,15]. Among of which, the LSM has attracted much attention in the field of structural design, due to their unique characteristics in evolving structural boundaries [12,13]. However, several complicated numerical issues are involved in the implementation of the conventional LSMs when applied to structural optimization, mainly due to the direct solution of the Hamilton-Jacobi partial differential equation (H-J PDE) [12,13], such as the extension of the boundary velocity, the re-initialization and the CFL condition.
To remove the numerical difficulties in the traditional LSMs, several variants of LSMs have been developed [16][17][18][19][20]. In particular, the parametric level set method (PLSM) [16,17] has been emerged as a powerful alternative approach for topology optimization, because it can not only inherit the unique characteristics of traditional LSMs, but also eliminate their numerical difficulties. The key concept of the PLSM is to interpolate the level set function by a system of the compactly supported RBFs (CSRBFs).
Topology optimization has been extensively used not only in the structural design, but also in the design of cellular composite materials together with the homogenization method [28]. Sigmund [29] proposed the inverse homogenization method to design the unit cell, periodically distributed in the micro-structured materials with the prescribed constitutive parameters, as the effective properties of microstructure depend on the topology of the material microstructure rather than the intrinsic composition. After that, many different optimization methods have been developed to create a range of micro-structured composites with the tailored or extreme properties, like the extreme material properties [30,31], the maximum stiffness and fluid permeability [32,33], negative Poisson's ratio [23,34], and exotic thermomechanical properties [26]. The numerical homogenization method has been mostly used to evaluate the effective properties.
It is noted that all the above works including structures and micro-structured composites focus on the design of monoscale structures. To achieve better structural performance, the multiscale design concept has been introduced into topology optimization [35]. Currently, the multiscale topological design methods can be roughly classified into two branches. The first is developed under the assumption that the identical material microstructures are uniformly distributed in the macrostructure [36][37][38][39][40][41]. Hence, the overall macro domain is characterized with the same effective properties. The earlier work is that the identical microstructure is uniformly configured within the macrostructure subject to the macro loads and boundary conditions, while the overall topology of the macrostructure kept unchanged [39,40]. Later, the macrostructure was also topologically optimized, aligning with the optimization of the representative microstructure [37,38]. This kind of multiscale designs for porous composites can save computational time, and no connectivity issue raised because all the material microstructures are identical in size, shape and topology. However, this assumption limits the design freedom for the further improvement of the performance. The second type of methods assume that different material microstructures can exist within the macrostructure [42]. Some studies have focused on heuristic multiscale topology optimization of functionally graded materials with different types of the representative microstructures, subject to a fixed topology of the macrostructure [43][44][45][46][47]. After that, the multiscale topology optimization has been developed for concurrent designs of both the macro and micro topologies [48][49][50][51][52][53]. However, the computational cost is still prohibitive and manufacture is challenging, because the topologies of all 4 microstructures may demand optimization. Some works also provide the criteria for distributing multiple material microstructures in the overall macrostructure, such as a trial-and-error criterion to divide the geometrical domain [54] or the principal stress distribution [55]. However, it is noted that how to devise a uniform criterion to distribute multiple microstructures is still not available in the multiscale topology optimization. Hence, an efficient multiscale topology optimization method for composites, to consider the overall topology of the macrostructure and the topologies of different microstructures, as well as their overall distribution in the macrostructure, is still in demand.
The main motivation of this paper is to develop a new multiscale topology optimization method for the design of porous composites considering three design pillars, in order to meet both ends for the structural performance and computational efficiency. Firstly, a free material distribution optimization is developed, where the VTS method [56] is applied to generate an element-wise varied distribution of the densities and then a regularization mechanism with the defined thresholds transforms the initial distribution into a domain-wise pattern. The discrete element densities in the domain-wise pattern can determine the overall   3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65   5

Parametric level set method for porous composites
As shown in Fig. 1, the macrostructure contains two kinds of microstructures and each of them is periodically distributed within different regions of the macrostructure, where the global coordinate system x shows the macrostructure and the local coordinate system y denotes the microstructures. In the following, the superscript indicates the macroscale quantities, and the is related to the microscale quantities.

Level set-based implicit boundary representation
In the LSM, the structural boundary is implicitly embedded into the zero-level set of a higher-dimensional level set function (LSF) [11], as shown in Fig. 2.  [16,17], which restrains the applicability and the efficiency of the H-J PDE-driven structural optimization method.

Parameterization of the level set function
In the PLSM, the CSRBFs with C2 continuity [57] is employed here: where is the number of the CSRBFs knots in the macrostructure, and all the representative microstructures have the same number of the CSRBFs knots .
is the radius of the support domain of the basis function at the current knot, and the detailed information at different scales is given as: where controls the scale of the support domain at the macrostructure, and is the corresponding term in the microstructures. Then, the time-dependent LSFs are interpolated by a given set of CSRBFs that are fixed in the space and their expansion coefficients by (6) The series of CSRBFs are formulated as a vector at every knot: The expansion coefficients are expressed in the same manner, as follows: Once time and space are decoupled by the above interpolation, the term of CSRBFs is only dependent on space, and the term of expansion coefficient is only dependent on time. Substituting Eq. (6) into Eq. (3), the H-J PDEs are transformed into the time-space independent forms as: Then, the normal velocities and at two scales can be given by: In Eq. (10), it can be seen that the normal velocity fields have been naturally extended to the whole design domain at two scales due to the definition of the CSRBFs over the whole design domain. The H-J PDEs 8 at two scales have been converted into two new forms, respectively, where the ordinary differential equations (ODEs) are defined with a set of unknown expansion coefficients serving as the design variables.

Multiscale topology optimization for porous composites
The core idea of the multiscale topology optimization method for porous composites has been clearly illustrated in Fig. 3. The whole procedure can be involved into two stages: the free material distribution optimization and the concurrent optimization considering three design pillars in porous composites. As an example of a cantilever beam, the initial structure is defined in Fig. 3 (a).

Fig. 3. Multiscale topology optimization of porous composites
At the stage 1, the VTS [56] method is firstly applied to generate the element densities which are distributed in an element-wise pattern, as illustrated in the first plot of Fig. 3 13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65   9 group. Then, the average of the element densities within each group is utilized to represent the density of this group. The re-distribution of the regularized element densities is displayed in the second plot of Fig.   3 (b). It can be seen that the regularized element densities are distributed in a domain-wise pattern. Each sub domain is only characterized by one equivalent density.
After that, the regularized density in each region will be chosen as the maximum volume fraction for the corresponding microstructure. Only one representative microstructure will be topologically optimized and periodically arranged within each region. Hence, it can be confirmed that the distribution of the In conclusion, free material distribution optimization determines the number and overall distribution pattern of the representative microstructures to be designed within the macrostructure, while the topologies of both the macrostructure and the macrostructure are optimized in the concurrent optimization stage.
is the structural compliance. and are the external load vector and the global displacement vector, respectively. is the global stiffness matrix, and is the stiffness matrix of the solid element. is the displacement vector in the e  element.
is the volume constraint, where is the allowable volume fraction and is the solid volume fraction of the solid element.

The regularization mechanism
Assuming that the densities in the initial distribution are classified into groups (also the number of distinct material microstructures), the regularized mechanism is defined as: (12) where is the element density in the group. The and indicate the lower and upper thresholds for the element densities in the group, respectively. is the total number of the elements in the group, and is the regularized density of the group which is defined by the mean value of all the element densities in this group. Hence, the whole macro design domain is divided into sub domains, and each of them is homogeneously occupied by one equivalent element density.

Concurrent topology optimization
The concurrent topology optimization for the macrostructure and multiple representative microstructures to minimize the structural static compliance is developed by the PLSM with the numerical homogenization method, subject to the overall distribution of the multiple microstructures, given as:  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 11 (13) where is the macro design variable bounded by and . and are the corresponding lower and upper bounds of the design variable for the material representative microstructure. is the objective function defined by the structural compliance. is the total material volume constraint, which is calculated considering the macrostructure and its microstructures.
is the allowable volume fraction for the global volume constraint. is the volume constraint for the representative material microstructure, which is subject to the maximum volume fraction defined by the regularized density in the free material distribution optimization. H is the Heaviside function [12,13]. is the displacement field of the macrostructure. is the macro virtual displacement field belonging to the kinematically admissible space .
The equilibrium equation is developed by the principle of virtual work, where the bilinear energy term a and the linear load term l are given by: (14) where is the macro body force, and is the macro boundary traction. is the partial derivative of the Heaviside function, namely the Dirac function [58]. The homogenized effective elastic tensor of the representative microstructure is evaluated by the numerical homogenization method [28], as: (15) where is the area in 2D or the volume in 3D of the representative microstructure, and is the elasticity property of the solid material. is the initial unit test strain field, which contains three 12 or six unit vectors in 2D or 3D [29]. is the locally varying strain field induced by . The unknown displacement is calculated by the micro equilibrium equation under the initial unit test strain .
The weak forms for both the energy a and load l terms for the representative material microstructure are defined as, respectively: where is the micro virtual displacement field of the representative microstructure belonging to the micro kinematically admissible displacement space .

Design sensitivity analysis
In the proposed multiscale topology optimization, both the VTS and PLSM employed in the problems can be solved by many well-established optimization algorithms [21,22]. Hence, the sensitivity information of the objective and constraint functions are required. In the above optimization formulation, three kinds of design variables are existed, including the element densities in the free material distribution optimization, the expansion coefficients for macro and micro structures in the concurrent optimization.

Macro sensitivity analysis in the concurrent optimization
given as: The shape derivative of the Lagrangian function is derived by the shape sensitivity analysis [58] as: (20) where (21) and is the mean curvature in two dimensions [12,16,17]. Now, recalling the normal velocity field in Eq. (10) and substituting it into Eq. (21), and the shape derivative of the Lagrangian function L can be written as: Eq. (22) can be expanded as: where and are defined as: where (25) On the other side, the shape derivative of the Lagrangian function L can be gained based on the chain rule.
14 (26) Thus, the design sensitivities of the objective and constraint functions with respect to the expansion coefficients at macro scale are obtained by comparing the corresponding terms in Eqs. (24) and (26): In order to improve the numerical efficiency of the optimization, the design sensitivities expressed by the boundary integration scheme are transformed into the following volume integration [16,24], given as.

Micro sensitivity analysis in the concurrent optimization
Here, the first-order derivatives of the objective and constraint functions with respect to the micro design variables (the micro expansion coefficients) are computed based on the chain rule, given by: It can be found that the key to calculate the sensitivities in Eq. (29) is the derivative of with respect to the micro design variable, which can be derived via the shape derivative [23,58], given as: Recalling the micro normal velocity field in Eq. (10), and substituting it into Eq. (30) will yield: Based on the chain rule, the first-order derivative of with respect to time t is given by:  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65   15 According to the Eqs. (30) and (31), it can be found that the first-order derivative of with respect to the micro expansion coefficient can be given as: Finally, the first-order derivatives of the objective function with respect to micro design variables can be attained by substituting Eq. (33) into Eq. (29). Similarly, the first-order derivative of the volume constraint with respect to the micro design variable is given as:

Numerical implementations
Here, the OC algorithm [21] is used to update the design variables in the multiscale design. The details for updating mechanism in the OC method can refer to [17]. The flowchart of the proposed multiscale topology optimization method is shown in Fig. 4. The first stage aims to seek for the appropriate distribution of the representative material microstructures with allowable volume fractions in the macrostructure by the VTS method combined with the regularization mechanism in Eqs. (11) and (12).

Numerical Examples
In this section, several numerical examples in 2D and 3D are provided to showcase the effectiveness of the multiscale topology optimization method for porous composites. The "ersatz material" model [13] is used to approximately calculate the related properties of those elements cut by the moving boundary in the finite element analysis at two scales. It is known that material microstructures have no specific sizes but should meet the periodic and multiscale conditions to enable the application of the homogenization method [28]. In all examples, the normal sizes of material microstructures are set as 1 mm. The Young's modulus is 210GPa, and the Poisson's ratio is 0.3. In the OC algorithm, the moving limit is m = 0.002 and the damping factor is set as 0.3. The termination criteria for the free material distribution optimization is that the difference of the objective values between two successive steps is less than 1e-4, and the concurrent optimization will terminate if the difference is less than 1e-4 or the maximum 200 steps are reached. 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65   17 In this example, the cantilever beam is used to study the efficiency of the proposed multiscale topology optimization method. As displayed in Fig. 5, the beam is fixed along the left side and a force (F=5e5 N) is loaded at the middle point of the right side, with the L=1.4m and H=0.7m. The macrostructure is discretized with 140×70 finite elements with four corner nodes, and 80×80 finite elements are used to discretize the microstructures. The dimensions of the macro finite element are equal to 1cm, ten times as the sizes of the microstructures, which can satisfy the conditions of the numerical homogenization to evaluate the effective property of the microstructures. The maximum volume fractions and are defined as 32% and 50%, respectively.

Multiscale design of the cantilever beam (a) Free material distribution optimization
The distribution of the element densities by the VTS method is displayed in Fig. 6. It can be seen that the   Table 1 lists a scheme (S1) of the regularization mechanism with the defined thresholds. At the same time, it is assumed that the regularized element densities will be 0 if they are located between 0.0 and 0.2, 3 4 5 6 7 8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 18 and the regularized element densities will be 1 if they are within 0.8-1.0, and the intermediate densities are regularized based on Eq. (12). As given in Fig. 7, the number of the distinct densities is reduced to 5, and the number of the intermediate densities after regularization is equal to 3. Then, the macro design domain is divided into five sub regions plotted with different colors (white, blue, green, red and black), and each of them is uniformly arranged by a number of identical element densities. In this way, each sub macro domain is occupied by a representative material microstructure. It should be noted that the regularized discrete densities are chosen as the material volume constraints for the corresponding microstructures to be devised in the latter concurrent topology optimization. Table 1. Scheme 1 (S1) of the regularization mechanism Scheme The defined thresholds in different groups The topology of the macrostructure and topologies of five representative microstructures are concurrently optimized, in terms of the overall distribution of these five microstructures in the macro domain in Fig. 7.
The initial design of the macrostructure is given in Fig. 8 (a). The optimizations for five representative microstructures employ the same initial design, defined in Fig. 8 (b). Fig. 9 gives the optimized topology of the macrostructure. Table 2 shows the details of four microstructures (excluding the microstructure with the void), including the optimized micro topologies and the 10×10 repetitive microstructures, as well as the corresponding effective elastic tensors after the homogenization.  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 19 Fig. 9. The optimized topology of the macrostructure In Fig. 10, the final multiscale design of the cantilever beam is provided. It can be found that the  This phenomenon also complies with the fact that the topology optimization for stiffness design should maintain the connectivity of loading propagation in the macro design domain [45], under the given supports and loading conditions. Finally, it can be found that the optimized topologies at two scales are all featured with the smooth boundaries. Hence, the proposed formulation can fully make use of the favorable features of the VTS method and the PLSM.  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65   21 The convergent histories are displayed in Fig. 11, in which the structural mean compliance and the total volume fraction are shown in Fig. 11 (a). The iterations of volume fractions for five representative material microstructures are illustrated in Fig. 11 (b). It can be seen that the iterations are smooth and fast, which shows that the proposed multiscale design method has a high optimization efficiency. The multiscale design optimization remarkably reduces the number of microstructures to be designed, and the PLSM removes several strict numerical implementations in most conventional LSMs.

Influence of the regularization mechanism
In order to address the influence of the free material distribution optimization on the structural performance, we define four regularization schemes (S2-S5) in Table 3. All the design parameters keep the same as Section 6.1.1. The defined four schemes are applied to regularize the element densities in Fig.   6, and we can achieve four different regularized distributions presented in Fig. 12. Similar to Fig. 7, the regularized element densities are distributed in a multi-domain pattern. Meanwhile, the macrostructure is divided into more sub domains, subject to a finer regularization of the continuous distribution of element densities. Hence, four overall distributions of the representative microstructures with allowable volume fractions to be designed are obtained for the latter concurrent topology optimization.
The concurrent design for the topology of the macrostructure and the topologies of multiple microstructures are performed in four cases under four different distributions of the representative 22 microstructures. Hence, four different multiscale designs are presented in Fig. 13. Similar to the design in  Moreover, the optimized objectives in four multiscale designs gradually decrease with the increasing of the number of the representative microstructures, as J 2 >J 3 >J 4 >J 5 . The main reason is that a finer regularization generates a finer distribution of multiple representative microstructures in the macrostructure, so that the design freedom of the structural performance is expanded [45,50]. In Fig. 13 (d), it is noted that the final multiscale design 5 does not comprise the microstructure 2 with the volume fraction 15%. It may be because the microstructure 2 contributes less to improve the overall stiffness for the macrostructure, and it is not included in the final design. In this way, we may conclude that when the types of the microstructures reach to a certain number, a further increase of their diversity cannot lift the performance of the macro structure. 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65   23 In this subsection, we will study the influence of the symmetry condition on the structural performance.

Influence of the symmetry condition
All the design parameters are consistent with Section 6.1.1. The S1 scheme of the regularization mechanism is used to reduce the number of the intermediate densities. The regularized distribution should be same as Fig. 7. However, in order to remove the effect of the symmetry of the macro loads and boundary conditions on the optimization of the microstructures, we partition the regularized distribution along the transverse central axis of the cantilever beam. The final distribution is shown in Fig. 14. It can be seen that the macrostructure is divided into eight sub domains.

Comparison with conventional multiscale design
To display the positive features of the proposed multiscale design, the cantilever beam is optimized by the conventional multiscale design [38,40] that the identical microstructure is uniformly and periodically arrayed within the macrostructure. It is noted that the concurrent topology optimization formulation in this paper will be degraded to the conventional multiscale design, if the distribution only contains a number of identical microstructures. The design parameters are consistent with those given in Section 6.1.1.
The final multiscale topology optimization design in shown in Fig. 16, and the objective value is J 0 =207.96. Hence, it can be seen that the optimized structural performance by the proposed multiscale design is much better than the conventional multiscale design (J 0 >J 2 >J 3 >J 1 >J 4 >J 5 ), which shows the benefit of the current multiscale topology optimization method. The free material distribution optimization can guide an adaptive configuration of the microstructures in the macrostructure for the enhancement of the performance.  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65   25 In this example, we further study the efficiency of the proposed multiscale topology optimization method.
The Michell structure (Fig. 17) is loaded with a concentrated vertical force (F=1e6 N)

Free material distribution optimization
In the distribution optimization, we firstly employ the VTS method to optimize the element densities, which is shown in Fig. 18 (a). The element densities are continuously varied from 0 to 1 in the macrostructure, which leads to a large amount of microstructures to design. We use the S1 scheme in Table 1 to regularize the continuously distributed element densities. Meanwhile, the symmetry of the macro loads and boundary conditions on the design of multiple microstructures is not considered.
Without the consideration of the symmetry may improve the stiffness of the structure. Hence, the overall distribution of the regularized element densities is partitioned along the vertical axis, as illustrated in Fig.   18(b). It can be seen that the element densities also distribute in a domain-wise pattern, also similar to the results in Section 6.1.

Concurrent topology optimization
The topology of the macrostructure and topologies of multiple microstructures are concurrently optimized, subject to the overall distribution gained in Fig. 18 (b). The initial design of the macrostructure is defined in Fig. 19 (a). The designs for multiple representative microstructures start from the same initial guesses displayed in Fig. 19 (b), due to the topologies of microstructures not known.
The optimized topology of the macrostructure is shown in Fig. 20. The optimized results of representative microstructures (excluding the microstructures with void) are clearly listed in Table 4, including the optimized topology, 10×10=100 repetitive microstructures, and the corresponding homogenized elastic tensors. In Table 4, the material microstructures 2, 4, and 6 are symmetric with respect to the microstructures 3, 5, 7, respectively. It can also be seen from the effective elastic tensors of the microstructures after homogenization.   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64   The optimized multiscale design of the Michell structure is provided in Fig. 21. It is noted that the microstructures with solid and void are excluded in the final design to simplify the description. We can see that the optimized topology of the macrostructure consists of eight sub domains, indicated by eight colors. Each sub domain is periodically distributed by the corresponding representative microstructure.
Hence, the effectiveness of the proposed multiscale method considering three design elements is further demonstrated. Additionally, it can be seen that the upper and lower edges of the Michell structure are occupied by the solid microstructures, in order to provide the sufficient stiffness and limit bending deformation in the Michell structure. It is also noticed that the connectivity of microstructures is basically maintained in the macro domain to ensure a reasonable macroscopic response [45].

3D supported structure
In this example, we study the effectiveness of the proposed multiscale optimization method on 3D porous composites. A 3D supported structure is given in Fig. 22, which is fixed at four corners in the bottom surface and loaded with a force (F=1e7 N)