1 Introduction
Functional equipments and devices infilled with welldesigned microstructures have been widely used in mechanical, thermal, acoustics and optical aspects Lakes (1993); Dong et al. (2017); Sigmund and Torquato (1996); Kushwaha et al. (1993); Liu et al. (2015); Vogiatzis et al. (2018). Driven by the development of additive manufacturing techniques, the fabrication of complex microstructures has become increasingly accessible, and this naturally stimulates the development of the corresponding intelligent/automatic design algorithms. Among them, topology optimisation (TO) methods, which aim at optimally distributing a certain amount of materials within a design domain for achiveing optimised properties/functions, have demonstrated their effectiveness for microstructural design Aage et al. (2017) and received great increasing research attention in recent years.
Computational approaches for microstructural design normally fall into the category of multiscale topology optimisation, i.e., at least two length scales are involved: a macroscopic scale on which the overall structural response is concerned and a microscopic scale on which the microstructural members are resolved. Doing topology optimisation straightforward on the microscopic level, albeit a number of valuable works achieved Alexandersen and Lazarov (2015); Liu et al. (2017), usually poses high demand on computational resources. One way to manage this efficiencyversusaccuracy dilemma, is to develop multiscale approaches in a viewpoint of scale separation. The basic idea lying behind the scaleseparation methods is to treat the design domain as a macroscopic continuum, and the locally homogenised properties of the continuum at each macroscopic point are quantified by solving a set of boundary value problems defined on the microscopic unit cell.
The mathematical foundation beneath the TO approaches based on separation of scale is the asymptotic homogenisation (AH) results for periodic structures Bensoussan et al. (1978). They were pioneerly implemented in the optimal design of spatially periodic microstructures Bendsoe and Kikuchi (1988); Bendsoe (1989); Sigmund (1994). The idea was later implemented for the concurrent design on multiple scales Rodrigues et al. (2002); Liu et al. (2008), and was further generalised for multifunctional designs Niu et al. (2009); Deng et al. (2013); Yan et al. (2016) and for the consideration of loading uncertainty Deng and Chen (2017).
The aforementioned works mainly consider the optimal design of microstructures whose constituting cells are presumed spatially periodic. However, a huge amount of highperformance microstructures, both naturally formed Sanchez et al. (2005); Fratzl and Barth (2009); Meyers et al. (2013) and manufactured Jorgensen et al. (1998); Cheng et al. (2018), displays a spatially varying nature. As for the optimal design of such configurations filled with graded microstructures or “graded microstructural configurations”(GMCs), considerable efforts have been denoted to the modification of the existing AH formulations simply devised for spatially periodic configurations Zhou and Li (2008); Radman et al. (2013, 2014); Wang et al. (2017); Cheng et al. (2019); Gao et al. (2018); Zhang et al. (2018, 2019)
. Such treatments generally see their limitations on three aspects: 1) constituting cells are normally of rectangle or cuboid shape, thus the spatial variance is only permitted along directions parallel to one of the cell edges; 2) smooth connectivity across cell boundaries can not be ensured; 3) the computational accuracy of the calculated structural responses can not be controlled.
To circumvent these challenging issues in using AHTO for the minimisation of GMC compliance, a novel scaleseparation scheme was recently proposed Groen and Sigmund (2018); Groen et al. (2019). They first searched for an optimal continuum distribution of laminar orientations whose homogenised elastic properties are obtained either by AH or microscale precomputation. Then the details of GMC can be fully resolved referring to the projection method Pantz and Trabelsi (2008). The methods have shown their effectiveness in the compliance optimisation involving GMC composed of thorder laminate cells.
An alternative route to substantially improve the AHTO methods for effective GMC optimisation, is to systematically rederive the asymptotic homogenisation results suitable for configurations infilled with graded microstructures. This motivates the work by Zhu et al. Zhu et al. (2019). They performed asymptotic analysis to rebuild the mathematical foundation of the original framework, so that the description, compliance computation and topology optimisation of GMC are harmonically integrated in the developed framework. In order to distinguish it from the traditional AHTO formulation, the solution framework Zhu et al. (2019) is termed as “AHTO plus” in this article.
The method of AHTO plus established a solid theoretical foundation to address the bottleneck issues limiting the conventional AHTO approaches for effective GMC optimisation Zhu et al. (2019). From a practical viewpoint, however, improvement is still highly desired so that the AHTO plus method can be widely applied for the design of GMC, mainly due to the following two reasons. Firstly, no explicit formulation of sensitivity analysis has been derived in accordance with the AHTO plus approach. Secondly, to maintain high accuracy, one has to compute the microscopic equilibrium problem as many times as the number of the Guassian points in the finite element (FE) mesh used for computing the macroscopic structural response in each iteration step. Due to these difficulties, Zhu et al. Zhu et al. (2019) had to adopt a less efficient generic algorithm and only tested the linearised form of the AHTO plus formulation for numerical implementation.
In face of these challenging issues lying on the application side of the AHTO plus formulation, this article is presented. Note that solving the microscopic cell problem, which comprises of the most timeconsuming part of the AHTO plus calculation, only depends on macroscopic positions. In the present work, a divideandconquer strategy is proposed to make the best use of this feature of pointwise computation. By decomposing the macroscopic computational domain into several zones, the evaluation of the microscopic cell problems can be carried out in a fully parallel way. This significantly improves the computational efficiency. We will show that for a threedimensional GMC optimisation problem which has rarely been examined in literature, possibly because of its high computational cost, it takes an eightprocessor computer very few hours to output an optimised design with the use of the present zoning scheme. Detailed analysis is also performed to thoroughly assess the efficiency brought out by the use of the proposed zoning strategy. Moreover, the sensitivity analysis associated with the AHTO plus formulation is also performed in the present work, and this further speeds up the proposed algorithm. Finally it will be shown that the present method is not only useful for designing bulk structures filled with graded microstructures, but also has the potential in designing manifoldlike structures (e.g., surfaces) decorated with graded microstructures.
The article is arranged as follows. In Sec. 2, the AHTO plus formulation is firstly outlined, and the challenging issues limiting its potential application are discussed. Then a socalled divideandconquer strategy empowered by parallel computation, which aims at resolving the challenging issues associated with the original AHTO plus approach, is introduced in Sec. 3. This is followed by the derivation of the gradientbased sensitivity results associated with the AHTO plus formulation in Sec. 4. In Sec. 5, two dimensions and three dimensions examples are presented to illustrate the effectiveness of the proposed approach. In this section, the corresponding numerical issues are also examined in depth. The potential of the proposed approach to design manifoldlike structures with GMC is discussed in Sec. 6. Finally, some concluding remarks are provided in Sec. 7. In this article, indexation is needed in various occasions with the following rules applied. English letters in lower cases are used for the indexation associated with spatial dimensions; English letters in upper cases are used for the indexation associated with subdomains; Greek letters are used for the indexation associated with design variables.
2 Outline of AHTO plus formulation and challenging issues associated with its numerical implementation
In this section, we shall first review the AHTO plus formulation proposed in Zhu et al. (2019), which comprises of three parts, namely, the representation of GMC, homogenisation formulae of GMC compliance and its topology optimisation, respectively. Then the challenging issues associated with the efficient numerical implementation of the AHTO plus solution framework are addressed.
2.1 Outline of the AHTO plus formulation
2.1.1 Multiscale representation of a configuration filled with quasiperiodic graded microstructures
Graded microstructural configurations generally possess a multiplelengthscale feature. A configuration filled with graded microstructures is composed of members of characteristic length , while it also has an overall size characterised by another length scale parameter . In general,
(1) 
In this article, we term the length scale characterised by the “microscale”, and term the length scale characterised by the “macroscale”.
Within any topology optimisation framework, a GMC in can be described by a topology description function (TDF) . Suppose that a porous region filled with certain graded microstructures is denoted by , as shown in the bottomright panel of Fig. 1, and is used to denote the solid part of . Then a TDF can be assigned in the whole design domain , such that for all and for all .
Given the normally large contrast in the two length scales characterised by and , respectively, extremely fine mesh and a large amount of data storage space are usually needed for describing the GMC of numerically. To resolve this issue, a multiscale scheme for the representation of quasiperiodic graded microstructures was introduced Liu et al. (2017) and further modified Zhu et al. (2019). The key idea is illustrated in Fig. 1. For a given GMC as shown in the bottomright panel of Fig. 1, a macroscopically smooth mapping function is first defined to map it to a spatially configuration whose constituting cells are uniformly of length scale . Then for the spatially configuration in the space measured by (as shown in the middle panel of Fig. 1), one simply needs to assign a TDF, say , to its generating cell as shown in the bottomleft panel of Fig. 1. Here the microscale length parameter appearing in the denominator is to make the TDF ’s domain of definition be a unit cell denoted by (=2 or 3 unless otherwise stated). In this way, the TDF of can be represented by the composition of a macroscopically smooth mapping function and the TDF defined in the unit generating cell, i.e.,
(2) 
Although in a seemingly simple form, the multiscale representation given by Eq. (2) has seen its abundance in its describable GMCs, and one may refer to Zhu et al. (2019) for more illustrative examples. Besides, in a view point of numerical implementation, Eq. (2) enables the storage of the mapping function on a grid far coarser than the microscale parameter , while fine mesh is only necessary for representing the TDF in one unit cell. This treatment has the advantage of significantly squeezing the digital space used for GMC.
2.1.2 Asymptotic analysis of elasticity problems involving GMC
The governing equation of a typical elasticity problem involving GMC can be formulated as
(3) 
where
denotes vector of the displacement field,
is vector of the body force density andis the fourthorder elasticity tensor of the solid materials constituting the GMC. In this work, the Einstein summation convention rule is applied unless other stated. Nonetheless, given its multiscale nature, evaluating the stress response of a structure constituting by graded microstructures usually requires extremely fine FE meshes, and this is often computationally intractable.
To address this issue, the continuum limit of Eq. (3), as the ratio of the two length scale parameters, tends to zero, has been studied Zhu et al. (2019). They showed that the original multiscale problem is asymptotically equivalent to a homogenised problem where the region occupied by GMC is treated as a solid continuum covering the whole domain of . Here a superscript “” is affiliated with an variable to indicate that it is a homogenised quantity. Asymptotic analysis shows that the homogenised displacement field satisfies a macroscale equilibrium relation:
(4) 
with appropriate homogenised boundary conditions satisfied, where (, , , , , ) is the homogenised elasticity tensor of the equivalent continuum and (, , ) measures the homogenised body force density which is the (vector) values averaged over the generating cell .
The homogenised elasticity tensor is calculated by
(5) 
where denotes the porous solid region within the unit generating cell and denotes its measure; (, , , ) is the Jacobean matrix associated with the mapping function
(6) 
measuring the degree of spatial variance of the GMC; is a thirdorder tensor satisfying a microscale cell problem
(7) 
equipped with periodic boundary conditions, where is given by
(8) 
Once the homogenised displacement field is obtained, the homogenised stress field can be computed by
(9) 
for , ,, .
2.1.3 The optimisation formulation
In the present work, the homogenised compliance is minimised by considering the available volume of the solid material in a specified design domain filled in by graded microstructures. Under this circumstance, the optimisation problem can be formulated as follows.
(10)  
s.t.  
where is the homogenised structural compliance; is vector of virtual displacement field; and are the function spaces for the macroscopic mapping functions and the TDF within the unit generating cell , respectively; is the boundary part of on which displacement boundary condition is imposed; is the boundary part of on which a traction of is applied, and is outer normal vector on the domain boundary . Note that other constraints, such as volume constraints, may also be needed in the actual numerical implementation of the AHTO plus formulation.
In the abovedescribed AHTO plus formulation, the representation, stress analysis and compliance minimisation of GMC are systematically integrated. Since the formulation is derived using asymptotic analysis, computational efficiency and certain accuracy should be simultaneously attained.
2.2 Challenges in the implementation of the AHTO plus approach
The AHTO plus approach Zhu et al. (2019), however, still suffers from being effectively applied mainly for two reasons.
Firstly, the AHTO plus formulation involves solving two linear partial differential equation systems: the macroscale equilibrium equation (
4) for the homogenised displacement field , and the microscale equilibrium equation (7) for the thirdorder tensor , , and , , . Note that Eq. (7) includes , the matrix measuring the degree of GMC’s spatial variance, which is a function of the macroscopic position . This means that one needs to solve the microscopic equilibrium equation (7) as many times as the number of Guassian points in the adopted FE mesh used for solving the macroscopic equilibrium equation (4). This severely reduces the overall computational efficiency.To circumvent this challenging issue, only the linearised form of the AHTO plus formulation has used for simulation Zhu et al. (2019). Upon linearisation, the computed compliance is only accurate for GMC that are produced by slightly perturbing a spatiallyperiodic configuration. As a result, the design freedoms reduce significantly. For example, the linearised AHTO plus method simply outputs an optimised twodimensional cantilever filled with graded microstructures, whose compliance value is still two times larger than that computed by using the projection method proposed in Groen and Sigmund (2018). Therefore, speeding up the pointwise computation of the microscale equilibrium equation (7) is the key to realize the full potential of the AHTO plus approach.
Furthermore, it is also worth noting that no sensitive analysis has been considered with respect to the AHTO plus formulation thus far. As a result, one has to employ resourcedemanding algorithms, such as generic algorithm, for optimization Zhu et al. (2019). This, however, also inevitably prevents the efficient numerical implementation of the AHTO plus method.
Facing these challenges, this article is presented. In consideration of the generic consistency of the AHTO plus method with parallel computation. In the following section, a zoning strategy is proposed to essentially accelerate the pointwise computational of Eq. (7). Moreover, sensitive analysis is also carried out accordingly, so as to facilitate the highefficiency solution of the corresponding optimisation problems.
3 A parallel divideandconquer strategy
In this section, a parallel divideandconquer scheme (or alternatively a zoning scheme) is proposed to speed up the numerical implementation of the AHTO plus formulation.
3.1 Domain decomposition with parallel computing
The most timeconsuming part when implementing AHTO plus approach, as mentioned above, is to solve the pointwise microscale cell problem governed by Eq. (7). The reason is that Eq. (7) depends on , which varies within the homogenised continuum region . Hence usually one has to solve Eq. (7) at every Guassian point in the FE mesh.
The proposed divideandconquer strategy is to divide the macroscopic region into subdomains as shown in Fig. 2.
Then the macroscopic continuum is considered as constituting by subdomains made from homogenous materials. Thus the effective elasticity tensor associated with each subdomain, denoted by , takes constant values. To derive , a representative point with respect to each subdomain, say for the subdomain , is selected as shown in Fig. 2. Then we use the homogenised approximately elasticity tensor evaluated at to represent the elastic property of the whole subdomain. Mathematically, by introducing a constant Jacobean matrix at , the corresponding introducing cell problem for solving the thirdorder tensor , , , , , , can be written as
(11) 
where and . When solving Eq. (11), periodic boundary conditions imposed on should be utilised.
Finally, the homogenised elasticity tensor associated with the subdomain can be calculated as (in subdomainwise form)
(12) 
Compared with the original cell problem governed by Eq. (7), the coefficients in Eq. (11) which are contained in no longer depend on the macroscopic position, but take constant values within each subdomain. This means that one simply needs to solve Eq. (7) as times recalling that here is the total number of subdomains.
Now the elasticity tensor in the macroscopic equilibrium equation (4) becomes piecewisely constant with the homogenised displacement field continuous within the whole domain . Hence the homogenised compliance of GMC can be approximated by
(13) 
The abovelisted formulation provides an approximation to the original AHTO plus formulation. In theory, the compliance computed from this zoning scheme should converge to that obtained from the original formulation outline in Sec 2.1.3, as the subdomains become more coincident with the finite elements for macroscopic computation. However, it will be shown numerically that good accuracy can be attained with the use of a certain number of subdomains for twodimensional problems.
Note that when solving the microscopic problem described in Eq. (11), there is no data exchange between computations taking place in different subdomains. This indicates that the zoning scheme may become extremely efficient if parallel computing strategy is employed. If we define a single “task” by the process of calculating by solving the microscale Eq. (11) in one subdomain, then one can assign on each cluster equal number of tasks, and run them in parallel. Another feature of the zoning scheme favoured by computational parallelism is that, for each task, the input is and the output are , both requiring little storage space if compared with that needed for FE computation of Eq. (7). It will be numerically shown in Sec 5.2.2 that the simulations get accelerated dramatically with the use of the proposed parallel computation scheme.
3.2 Numerical implementation of the zoning scheme
To ensure the smoothness of the GMC, the mapping function is expressed by
(14) 
where is a set of basis functions and is a set of weight parameters. Then the Jacobean matrix can be written by
(15) 
in the whole macroscopic domain , where denotes the total number of basis functions.
In the present work, we adopt polynomial functions as the basic functions, although other choices, e.g. trigonometric functions, are also applicable. Here we assume
(16) 
for and is thus expressed by
(17) 
where symmetrical conditions are imposed , . It is noted that all subdomains share the same expression for and a single TDF associated with the unit generating cell. This is homogenised to ensure smooth transition in microstructure across the subdomain boundaries.
In order to optimise the GMC, we search for the optimised coefficients of the mapping function given by Eq. (17) and the microscale TDF , such that the compliance given by Eq. (13) is minimised. The overall procedure of the “divideandconquer” optimisation scheme based on the AHTO plus formulation can be illustrated by the flowchart given in Fig. 3.
Firstly, evaluating Eq. (16) at each representative point and assigning its results for Eq. (11) give the microscale governing equation. Then, one can compute in parallel the homogenised elasticity tensor with the use of Eqs. (11)(12). Then one solves the macroscopic governing equation, and calculate the homogenised compliance by Eq. (13). With the use of the sensitivity analysis results to be presented in Sec.4, the values of all design variables are updated, if the optimality criterion is not met. Finally, the optimised design variables characterising the GMC are output.
4 Gradientbased sensitivity analysis
To further speed up the optimisation process, the gradientbased sensitivity analysis in accordance with the AHTO plus formulation should also be carried out. Actually, by resolving to the classical adjoint method, the sensitivity of the homogenised compliance can thus be expressed by
(18) 
where the vector contains all the design variables. As for the term , we can first follow Zhu et al. (2019) to rewrite given by Eq. (12) by
(19)  
The design variables associated with the AHTO plus approach can be classified into two groups: the macroscopic design variables identified by subindex
which are contained in the Jacobean matrix ; the microscopic design variables identified by subindex , which are the parameters controlling the TDF of the unit generating cell influencing the microscale solution , , , , , , to Eq. (11).For sensitivity analysis with respect to the macroscale variables , taking derivatives on both sides of Eq. (19) with respect to (where is the total number of macroscale design variables) gives
(20) 
where
(21) 
and
(22)  
respectively.
As for the sensitivity with respect to the microscale variables, taking derivatives on both sides of Eq. (19) with respect to (where is the total number of microscale design variables) gives
(23) 
By adopting a similar treatment as in Liu et al. (2002), the sensitivity with respect to the microscale variables can be finally expressed by
(24) 
where is the secondorder identity tensor.
In this work, the socalled MMCbased framework Guo et al. (2014); Zhang et al. (2016a, 2017, b); Guo et al. (2016); Lei et al. (2018) is adopted to describe the structural topology within the unit representative cell. Under this framework, the layout of material distribution in a specified design domain is described by a set of explicit hyperelliptical/ellipsoidal components whose explicit geometrical descriptions forms are available. It is also worth noting that the structural layout can also be described within other types of TO frameworks such as solid isotropic materials penalty (SIMP) method Bendsoe (1989); Zhou and Rozvany (1991) and Level Set Method Wang et al. (2003); Allaire et al. (2004), etc..
5 Numerical results and analysis
In this section, numerical results based on the proposed divideandconquer scheme are presented. Firstly, several benchmarked twodimensional examples will be considered, so as to demonstrate the effectiveness of the zoning strategy with parallel computing treatments. Then several threedimensional numerical examples will also be presented to illustrate the capability of the proposed approach for dealing with relatively large scale problems.
5.1 Problem setting
For all examples examined in this section, linearly elastic response and isotropic material are considered with the parameters chosen as follows: the (nondimensional) Young’s modulus and Poisson s ratio . For twodimensional examples, the specimen’s nondimensional thickness in the third dimension is fixed to be 1, and plane stress model is used unless otherwise stated.
The presented algorithm should work in a design space spanned by design variables defined on both scales. Nonetheless, it is found that allowing too much degrees of change in microscale variables (which control the structural topology within the unit generating cell) may result in GMC that are too complicated to be manufactured. Therefor, in order to demonstrate the effectiveness of the present algorithm, we choose to keep the material distribution within the unit generating cell unchanging throughout the whole optimisation process. Moreover, for FE discretised, fournode quadrilateral elements are used in all twodimensional examples while eightnode hexahedral brick elements are adopted in all threedimensional examples.
5.2 Twodimensional cases
5.2.1 First example and comparison with fine scale results
We start with an example of a short beam subject to a uniformly distributed pressure as shown in Fig.
4.The right side of the beam is fixed and a uniformly distributed pressure of magnitude 1 is applied on its top. The design domain is of size , and it is divided into identical square subdomains. In our simulation, we fix the material distribution in the microscale unit generating cell to be of “X”shape, and a mesh is employed for the computation of the macroscale equation. The volume fraction of the solid material is set to be 39.12%.
The optimised configuration is shown in the right panel of Fig. 6. Fig. 5 shows the compliance value of the GMC against the iteration steps during optimisation. Compared to the (initial) structure composed of spatially periodic microstructure, the value of the objective function for the optimised GMC drops nearly 70%. This is sensible from a mechanical viewpoint, since the fibers in the optimised configuration are assigned such that the pressure on the top is smoothly transformed to the fixed boundary in a more efficient way.
With the present homogenisation formulation, the structural compliance is computed to be 4,780, while a finescale simulation (1600 800) outputs a compliance of 4,450. The present homogenisation formulation seems to slightly overestimate a GMC’s compliance by no more than 10 percents. The main reason to such a deviation is because the asymptotic homogenisation result is no longer valid in a thinboundary region where the assumption of local periodicity breaks down. In theory, the resulting deviation is roughly at the order of the cell size divided by the domain size Dumontet (1985); Zhu et al. (2019), which is about 1/10 in this case.
5.2.2 Comparison with existing benchmarked results
More systematic study is performed toward another typical example of a short beam fixed on one side as shown in Fig. 7. At the middle point on its left edge, a unit vertical load is applied. The macroscopic design domain is a rectangle of length and width .
The upperhalf region is chosen as the design domain given the symmetry in geometry, and the design domain is considered and discretised by a FE mesh. For numerical implementation of the proposed zoning scheme, the half domain is divided into 82 identical subdomains of square shape, and the center of each square is chosen as the representative point of each subdomain. The buildingblock microstructure in the unit generating cell for this example, is taken to be of “X”shape as shown in the left panel of Fig. 8.
Based on the presented approach, the optimised GMC is shown in the right panel of Fig. 8, and the optimised GMC bears a strong similarity with that obtained based on the projection method Groen et al. (2019) as shown in Fig. 9. For instance, both structures are composed of fibers connecting the loaded region to the fixed boundaries.
Moreover, the energy density distribution and the maximum principal stress component are shown in Fig. 10, identifying the path transiting the load to the fixed boundary.
The GMC compliance by the present homogenisationbased formulation is 127.05 with a volume fraction of 34.3%, while the value by the method proposed in Groen et al. (2019) is roughly 114.02. The deviation in the compliance values by both methods can be rationalised as follows. Practically, a small solid region is added in the neighbourhood of the loading point Groen et al. (2019), thus the geometrically singular point load is treated as a uniform force distribution within the solid region. Thus the singular behaviour near the loading site gets much more regularised, and the resulting optimised compliance value is naturally smaller than that from the present method, where the load singularity is kept (to its best) within a single FE element.
5.3 Analysis of the zoning scheme
The proposed zoning scheme has been shown outputting considerably optimised GMCs upon loading. Now we further analyse the factors that characterising the present algorithm.
5.3.1 Effect by subdomain number
As discussed above, the proposed zoning scheme is expected to deliver more accurate results of structural response with denser domain partition. However, increasing the number of subdomains poses higher demand on computational resources. In this subsection, we will examine the appropriate partition strategy in the the zoning scheme, i.e., what is the best choice for the number of subdomains.
The short beam example shown in Fig. 7 is considered again. Optimisations have been carried out with 1, 2, 4, 16 and 64 identical squareshape subdomains, respectively. The total compliance value along with the corresponding optimised configuration about the upper part for each case is complied in Fig. 11. The computed compliance drops as the number of the subdomains increases. When the number of subdomains exceeds 16, the change in GMC compliance becomes tiny for this example, and the corresponding optimised GMCs are visually identical. This suggests that for the considered problem, using 16 subdomains is sufficient for finding a considerably optimised GMC.
Further analysis of the computational time against the number of subdomain is shown in Table 1. The averaged time needed for proceeding one step increases roughly linearly with the subdomain number. We further divide the time needed for carrying out one computational step into three parts, a) computing the homogenised elasticity tensor by solving the microscale equation (11)(12) with FEA; b) FEA for solving the macroscale equation (4); c) updating the design variables based on the MMA algorithm and the sensitivity analytical results derived in Sec. 4. The fraction for each part is shown in Table 1 for various cases of domain decomposition. Note that for all cases, calculating the homogenised elasticity tensors associated with all subdomains takes up most of the computational time in a single step. Moreover, the fraction of the computational time associated with this part exceeds 99% as the number of subdomains increase to 16. As discussed in Sec 3, the microscale problem governed by Eq. (11) is naturally suitable for being solved in parallel, and the effect of employing computational parallelism will be examined in the next subsection.
No. of  total  microscale  macroscale  data 

subdomains  time(s)  FEA(%)  FEA(%)  updating(%) 
1  2.2162  90.24  9.46  0.30 
2  4.1658  94.74  5.01  0.25 
4  8.1095  97.20  2.62  0.18 
16  31.8499  99.21  0.67  0.12 
5.3.2 Effect of computational parallelism
In section 3.1, we have theoretically demonstrated that the proposed zoning scheme should be significantly speeded up with the use of parallel computing treatment. Here we will numerically examine the effect brought by computational parallelism based on the short beam example shown in Fig. 7.
We fix the number of subdomains to be 16 and the optimisation is carried out with 1, 2, 4, 8 and 16 workers. The average time needed for proceeding one time step is plotted against the number of workers in Fig. 12. Here time is measured in relative terms, i.e. the average time using only 1 worker is set to be 1. It is observed that simply in a MATLAB environment, employing 16 workers saves roughly 80% of computational time compared to the case of a single worker.
This demonstrates the effectiveness of implementing the proposed zoning scheme with the use of parallel computing treatment. In particular, when the number of workers is 8, the overall time needed for optimisation is 2.17 hours (with 500 iteration steps). This means that the GM design becomes computationally tractable with the use of a relatively advanced personal computer. The speed up effect of using parallel computing is more significant for threedimensional cases and this will be shown in Sec 5.4.2.
5.3.3 Diversity of the describable GMC
By assuming that the microstructural cell takes an nthorder laminate configuration, the projectionbased approach Groen and Sigmund (2018) also provides an effective way to design GMC (in twodimensional space). The AHTO plus method, in contrast, does not pose any restriction in the form of involved microstructures. For the same shortbeam problems as shown in Fig. 8, optimisation results with other choices of unit generating cells are shown Fig. 13.
Note that when the mechanical behaviour of a GMC is in concern, cells of laminated structures are usually the theoretically optimal choice. However, when other properties of GMC are also of interest, one may need unit generating cells of nonlaminated shape. The AHTO plus method has also shown its potential of doing optimisation considering other physical performances of GMC.
5.4 Threedimensional cases
In this subsection, the stiffness optimisation of threedimensional GMC is studied. We will start with a quasitwodimensional example where the microstructural change only takes place in two directions, to validate the threedimensional algorithm. Then full threedimensional cases are studied, and the effectiveness of using parallel computing is also demonstrated.
5.4.1 A quasitwodimensional example
As shown in Fig. 14, a short beam of length , hight and thickness is considered. The beam is fixed on its right side and its left side is loaded with a line force distributed uniformly vertical to its middle plane. Similar to the twodimensional short beam case, only upper half of this part of the computational domain is considered due to symmetry of the design domain. The homogenised half domain is discretized by a FE mesh, and it is further divided into 4 identical cubic subdomains.
To validate our 3D algorithm, we first consider a quasitwodimensional situation, while material distribution is uniform along the direction of thickness. Hence the unit generating cell is set as two solid components crossing each other to form an “X”shape if viewed from its front side, as shown in the left panel of Fig. 15.
From a side view, the optimised GMC shown in Fig. 15 bears a strong similarity with the corresponding optimised 2D structure shown in Fig. 8. The slight difference may be due to the fact that the 2D case is treated as a plane stress problem, while the third dimension is taken into account for the 3D case shown in Fig. 15.
5.4.2 A full threedimensional example
Now we consider a more general case where the buildingblock microstructure takes a fully 3D configuration, as shown in the left panel of Fig. 16. Now the problem setting is assumed the same as that shown in Fig. 14, and the upper half domain is now divided into 4 subdomains. Again, the unit generating cells deform themselves in space such that solid materials are distributed along the path connecting the load to the fixed end. It is noted that within the AHTO plus framework, the optimised GMC displays a good smoothness and inner connectivity.
Two critical issues that make 3D cases essentially different form 2D ones should be mentioned. The degree of spatial variance for 3D cases is far greater than the 2D cases. As a result, some of the 3D GMC may became too complicated to be fabricated. Thus in order to ensure manufacturability, we suggest that a limited number of basis functions consitituting the mapping function should be adopted. For the simulation results presented in Fig. 15, we exclude the spatial variation along the thickness direction at the continuum level. How to choose basis functions for consistent with the manufacturability condition, is an interesting issue worth being further investigated (elsewhere).
Besides, as the optimisation proceeds, the unit generating cell tends to be elongated to a size that is comparable to the continuum domain size, as shown in Fig. 17. This has gone beyond the validity range of the AHTO plus formulation. To avoid its emergence, one can set bounds to the determinant of defined by Eq. (16).
In Table 2, the time consumed to carry out one single optimisation step, and its fraction value for different parts of the algorithm are listed. The microscale FEA costs more than 95% of the time consumed in one step. This fraction value is even higher than that in 2D cases shown in Table 1 for a same number of subdomains. This is because the microscale analysis, governed by Eq. (11), deals with a thirdoder tensor , , , , , . When the spatial dimensionality increases, the number of equation systems to be solved rise from to . Besides, the time fraction needed for microscale FEA increases with the number of subdomains, and the value is more than 99% with 4 subdomains. Therefore, the use of parallel computing should bring about even higher efficiency for 3D computation.
No. of  total  microscale  macroscale  data 

subdomains  time(s)  FEA(%)  FEA(%)  updating(%) 
1  74.7950  96.88  2.59  0.53 
2  137.9675  98.57  1.30  0.13 
4  284.8819  99.26  0.67  0.07 
6 GMC design on manifolds in a virtually mechanistic way
In this subsection, we will discuss in brief another potential application of the present method: how to generate a manifoldlike decorated with GMC in a virtually mechanistic manner.
Suppose that we look for a closed surface decorated with GMC, which has a strong resistance to a torque load, as shown in the right panel of Fig. 18. In the following is a virtually mechanistic route that may achieve it.
We start with a 3D laminate configuration as shown in the left panel of Fig, 18. Each layer consists of spatially periodic microstructures whose basic cells take an “X”shape.
Upon a virtually applied torque as shown in the lower panel of the Fig. 18, thus each layer should virtually fold itself to resist the applied torque, and the whole structure deforms like a Matryoshka Doll as shown in the upper panel of Fig. 18. From the mathematical point of view, this deformation process should correspond to a gradual change in the mapping function , e.g., given by Eq. (17). Note that during the whole process of virtually deformation, neighbouring layers are kept disconnected, since the mapping function does not change the topology within the unit generating cell. Thus as the desire deformation status is reached, we simply extract the corresponding layer to generate a GMC surface as demanded by taking the corresponding contour surface with respect to the mapping function.
The concept demonstrated above indicates that one may use the proposed scheme for generating manifoldlike structures composed by GMC from a relatively simple configuration. The key is to apply a virtual load to a set of sheets decorated with simple graded microstructures, and then computationally deform it to a desired state. Although the homogenisation formulation may not be as accurate in the case of manifolds because the surface effect Zhu et al. (2017) and large deformation kick in, the proposed algorithm may still output a reasonable manifold filled with graded microstructures, as in the example schematically shown in Fig. 18.
7 Conclusion
In the present work, a parallel divideandconquer scheme is proposed to help accelerate the automatically design of configurations filled with graded microstructures concerning their mechanical behaviour. The algorithm is built on a solid theoretical ground, the AHTO plus framework where the representation, response analysis and topology optimisation of GMCs are systematically integrated. We have shown that the new asymptotic framework possesses a strong favour for computational parallelism, and this feature is made the best use of in the proposed zoning scheme. The proposed method gets validated through a comparison with 2D benchmark examples, and the high efficiency brought by the use of parallel computing is demonstrated, especially through 3D examples that have been barely investigated in existing literature. We have also derived the gradientbased sensitivity analysis formulation, which helps further speed up the proposed algorithm. Finally, we have also discussed another potential use of the proposed scheme in generating manifolds filled with graded microstructures.
Future works based on the present algorithm are expected in two directions: its further improvement and its utilisation in wider fields. For improvement, the role played by the choice of the basis functions constituting the mapping functions should be examined in depth, and a concurrent design strategy with the materials topology within the generating cell also changing is also worth being studied in consideration of manufacturability conditions. Besides, the present treatment mainly considers GMCs whose volume fraction is kept as a constant throughout the macroscopic design domain. A natural generalisation of the present work is to consider the optimal performance of GMC whose volume fraction is allowed to vary in space. For its wider utilisation, the presented algorithm can be naturally generalised to the optimisation of GMCs for not only mechanical functions, but also other purposes, such as for thermal and acoustics use, and this highlights the potential in the development of the proposed algorithm on its application aspect.
Replication of results
All the data underlying the argument in the article are generated by locally devised MATLAB codes consisting of a set of systematically arranged subfunction modules (such as homogenisation etc.). We are willing to satisfy the reasonable and responsible demand for the data and the source codes underpinning the present article.
Acknowledgements.
We thank Ole Sigmund and Jeroen Groen for providing the comparative simulation example shown in Fig. 9. The financial supports from National Key Research and Development Plan (2016YFB0201601) from the Ministry of Science and Technology of the People’s Republic of China, the National Natural Science Foundation of China (11772076, 11732004, 11821202), the Fundamental Research Funds for the Central Universities (DUT16RC(3)091), Program for Changjiang Scholars, Innovative Research Team in University (PCSIRT) are gratefully acknowledged.Conflict of Interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.
References
 Gigavoxel computational morphogenesis for structural design. Nature 550 (7674), pp. 84–86. Cited by: §1.
 Topology optimisation of manufacturable microstructural details without length scale separation using a spectral coarse basis preconditioner. Comput. Methods Appl. Mech. Engrg. 290, pp. 156–182. Cited by: §1.
 Structural optimization using sensitivity analysis and a levelset method. J. Comp. Phys. 194 (1), pp. 363–393. Cited by: §4.
 Generating optimal topologies in structural design using a homogenization method. Comput. Methods Appl. Mech. Engrg. 71 (2), pp. 197–224. Cited by: §1.
 Optimal shape design as a material distribution problem. Struct. Optim. 1 (4), pp. 193–202. Cited by: §1, §4.
 Asymptotic analysis for periodic structures. Studies in mathematics and its applications, North Holland, Amsterdam. Cited by: §1.
 Functionally graded lattice structure topology optimization for the design of additive manufactured components with stress constraints. Comput. Methods Appl. Mech. Engrg. 344, pp. 334–359. Cited by: §1.
 Extra strengthening and work hardening in gradient nanotwinned metals. Science 362 (6414), pp. 1925. Cited by: §1.
 Concurrent topology optimization of multiscale structures with multiple porous materials under random field loading uncertainty. Struct. Multidiscip. Optim. 56 (1), pp. 1–19. Cited by: §1.
 Multiobjective concurrent topology optimization of thermoelastic structures composed of homogeneous porous material. Struct. Multidiscip. Optim. 47 (4), pp. 583–597. Cited by: §1.
 Topology optimization of anisotropic broadband doublenegative elastic metamaterials. J. Mech. Phys. Solids 105, pp. 54–80. Cited by: §1.
 Boundary layers stresses in elastic composites. In Studies in Applied Mechanics, L. Pierre (Ed.), Studies in Applied Mechanics, Vol. 12, pp. 215–232. Cited by: §5.2.1.
 Biomaterial systems for mechanosensing and actuation. Nature 462, pp. 442–448. Cited by: §1.
 Topological shape optimization of 3d microstructured materials using energybased homogenization method. Adv. Engrg. Softw. 116, pp. 89–102. Cited by: §1.
 Homogenizationbased topology optimization for highresolution manufacturable microstructures. Internat. J. Numer. Methods Engrg. 113 (8), pp. 1148–1163. Cited by: §1, §2.2, §5.3.3.
 Homogenizationbased stiffness optimization and projection of 2d coated structures with orthotropic infill. Comput. Methods Appl. Mech. Engrg. 349, pp. 722–742. Cited by: §1, Figure 9, §5.2.2, §5.2.2, §5.2.2.
 Explicit structural topology optimization based on moving morphable components (mmc) with curved skeletons. Comput. Methods Appl. Mech. Engrg. 310, pp. 711–748. Cited by: §4.
 Doing topology optimization explicitly and geometrically  a new moving morphable components based framework. ASME J. Appl. Mech. 81 (8), pp. 081009–1–081009–12. Cited by: §4.
 Spherical indentation of composite laminates with controlled gradients in elastic anisotropy. Int. J. Solids Struct. 35 (36), pp. 5097–5113. Cited by: §1.
 Acoustic band structure of periodic elastic composites. Phys. Rev. Lett. 71 (13), pp. 2022–2025. Cited by: §1.
 Materials with structural hierarchy. Nature 361 (6412), pp. 511–515. Cited by: §1.
 Machine learningdriven realtime topology optimization under moving morphable componentbased framework. ASME J. Appl. Mech. 86 (1), pp. 011004–1–011004–9. Cited by: §4.
 Frequencypreserved acoustic diode model with high forwardpowertransmission rate. Phys. Rev. Appl. 3 (6), pp. 064014. Cited by: §1.
 Additive manufacturingoriented design of graded lattice structures through explicit topology optimization. ASME J. Appl. Mech. 84 (8), pp. 081008–1–081008–12. Cited by: §1, §2.1.1.
 Optimum structure with homogeneous optimum trusslike material. Comput. Struct. 86 (13), pp. 1417–1425. Cited by: §1.
 Mapping method for sensitivity analysis of composite material property. Struct. Multidiscip. Optim. 24 (3), pp. 212–217. Cited by: §4.
 Structural biological materials: critical mechanicsmaterials connections. Science 339 (6121), pp. 773–779. Cited by: §1.
 Optimum structure with homogeneous optimum cellular material for maximum fundamental frequency. Struct. Multidiscip. Optim. 39 (2), pp. 115–132. Cited by: §1.
 A posttreatment of the homogenization method for shape optimization. SIAM J. Control Optim. 47 (3), pp. 1380–1398. Cited by: §1.
 Topology optimization of functionally graded cellular materials. J. Mater. Sci. 48 (4), pp. 1503–1510. Cited by: §1.
 Maximizing stiffness of functionally graded materials with prescribed variation of thermal conductivity. Comp. Mater. Sci. 82, pp. 457–463. Cited by: §1.
 Hierarchical optimization of material and structure. Struct. Multidiscip. Optim. 24 (1), pp. 1–10. Cited by: §1.
 Biomimetism and bioinspiration as tools for the design of innovative materials and systems. Nat. Mater. 4, pp. 277–288. Cited by: §1.
 Composites with extremal thermal expansion coefficients. Appl. Phys. Lett. 69 (21), pp. 3203–3205. Cited by: §1.
 Materials with prescribed constitutive parameters: an inverse homogenization problem. Int. J. Solids Struct. 31 (17), pp. 2313–2329. Cited by: §1.
 The method of moving asymptotes—a new method for structural optimization. Internat. J. Numer. Methods Engrg. 24, pp. 359–373. Cited by: §5.1.
 Computational design and additive manufacturing of periodic conformal metasurfaces by synthesizing topology optimization with conformal mapping. Comput. Methods Appl. Mech. Engrg. 328, pp. 477–497. Cited by: §1.
 A level set method for structural topology optimization. Comput. Methods Appl. Mech. Engrg. 192 (1), pp. 227–246. Cited by: §4.
 Concurrent design with connectable graded microstructures. Comput. Methods Appl. Mech. Engrg. 317, pp. 84–101. Cited by: §1.
 Multiscale concurrent material and structural design under mechanical and thermal loads. Comput. Mech. 57 (3), pp. 437–446. Cited by: §1.
 A new threedimensional topology optimization method based on moving morphable components (mmcs). Comput. Mech. 59 (4), pp. 647–665. Cited by: §4.
 Structural topology optimization through explicit boundary evolution. ASME J. Appl. Mech. 84 (1), pp. 011011–1–011011–10. Cited by: §4.
 A new topology optimization approach based on moving morphable components (mmc) and the ersatz material model. Struct. Multidiscip. Optim. 53 (6), pp. 1243–1260. Cited by: §4.
 Concurrent topology optimization for cellular structures with nonuniform microstructures based on the kriging metamodel. Struct. Multidiscip. Optim. 59 (4), pp. 1273–1299. Cited by: §1.

Multiscale concurrent topology optimization for cellular structures with multiple microstructures based on ordered simp interpolation
. Comput. Mater. Sci. 155 (4), pp. 74–91. Cited by: §1.  The coc algorithm, part ii: topological, geometrical and generalized shape optimization. Comput. Methods Appl. Mech. Engrg. 89 (1), pp. 309–336. Cited by: §4.
 Design of graded twophase microstructures for tailored elasticity gradients. J. Mater. Sci. 43 (15), pp. 5157. Cited by: §1.
 A novel asymptoticanalysisbased homogenisation approach towards fast design of infill graded microstructures. J. Mech. Phys. Solids 124, pp. 612–633. Cited by: On speeding up an asymptoticanalysisbased homogenisation scheme for designing gradient porous structured materials using a zoning strategy, §1, §1, §2.1.1, §2.1.1, §2.1.2, §2.2, §2.2, §2.2, §2, §4, §5.2.1.
 Gurtinmurdoch surface elasticity theory revisit: an orbitalfree density functional theory perspective. J. Mech. Phys. Solids 109, pp. 178–197. Cited by: §6.
Comments
There are no comments yet.