Infill graded microstructural (IGM) configurations have shown their vast engineering potentials in achieving exceptional material/structural properties Lakes_Nature1993 ; DongHW_JMPS2017 ; Sigmund_APL1996 ; Kushwaha_PRL1993 ; LiuC_PRAppl2015 ; Aage_Nature2017 . For IGM design, early treatments encompass the use of asymptotic homogenisation (AH) approaches (Bensoussan1978, )
. The basic idea is to treat a spatially periodic configuration as a continuum, whose elasticity tensor at every point is calculated through solving a single set of partial differential equations (PDEs) defined over one constituting cell. The homogenised formulation was then adopted for spatially-varying infill configurations, where every macroscopic pixel is associated with an individual cell(Bendsoe_CMAME1988, ; Bendsoe_StructOpt1989, ; Sigmund_IJSS1994, ). The method was further extended to other cases, such as concurrent multi-scale design Rodrigues_SMO2002 ; LiuL_CompStruct2008 , multi-functional design (NiuB_SMO2009, ; DengJD_SMO2013, ; YanJ_CompMech2016, ) and problems concerning load uncertainties DengJD_SMO2017 . However, several limiting issues arise against the treatments based on conventional AH formulation. For instance, the smoothness across neighbouring cells may not be guaranteed. Moreover, the asymptotic analytical results may become inaccurate, since the presumption of microstructural periodicity no longer holds.
For improvement, various types of methods for IGM representation and stress analysis have been proposed. An intuitive way of thinking is to pile up building-block structures, whose topology description functions (TDFs) originate from a single set of continuous functions. This can be achieved either explicitly or implicitly. The explicit way involves cell configurations composed of parameterised members (ZhouSW_JMaterSci2008, ; Radman_JMaterSci2013, ; Radman_CompMaterSci2014, ; ChengL_CMAME2019, ; LiQH_CompStruct2019, ). Thus an IGM configuration is generated by letting the corresponding parameters vary gradually (but discretely) in space. The implicit treatment is usually achieved through the introduction of a level-set function, whose multiple (but discrete) level-set contours are used to characterise the boundaries of a family of prototype cells WangYQ_CMAME2017 ; Gao_AES2018 ; ZhangY_CMS2018 ; ZhangY_SMO2019 . By either mean, graded change in microstructure should be achieved, but several challenging issues still persist. First, slight non-smoothness across the boundaries of neighbouring cells is still an unresolved issue, since only discrete values of the controlling functions are concerned. Second, the variation of the microstructural infill patterns is only permitted to take place along directions parallel to cell edges. This greatly squeezes the resulting IGM variety.
Recently, a number of novel approaches underpinning fast IGM design are proposed, where smoothness connections within IGM configurations can be fully guaranteed. Interestingly, the restriction on oriented rectangle/cuboid-shape cells is also removed in these recently proposed approaches. Besides simulations directly conducted on fine scales Alexandersen_CMAME2015 ; LiuC_JAM2017
, the novelly introduced treatments can be roughly classified into three categories: 1) conformal-mapping-based approaches; 2) the de-homogenisation approaches; 3) the modified asymptotic homogenisation topology optimisation (AHTO plus) approaches.
The conformal-mapping-based approaches Vogiatzis_CMAME2018 mainly concern decorating microstructural configurations on a given surface in three-dimensional space. In this scenario, conforming mapping is used to ensure that the local orthogonality in microstructure is preserved as a planar microstructure is deployed on the prescribed surface. As a cost, the Ricci flow equation, a useful differential equation defined on manifolds, needs to be numerically solved to a steady state.
The de-homogenisation methods Groen_IJNME2018 ; Groen_CMAME2019 ; Allaire_CompMathAppl2019 ; WuJ_IEEE2019 ; GeoffroyDonders_JCP2020 focus on cells of laminate configurations, which theoretically guarantee the solution optimality for most cases in compliance optimisation. They first determine an optimised distribution of the laminate volume fraction and orientation (on a coarser grid). Then the actual microstructure is resolved by solving another optimisation problem projecting the laminate configurations in alignment with the desired directions of the local principal stresses Pantz_SIAMJCO2008 .
The AHTO plus approaches ZhuYC_JMPS2019 describe an IGM configuration with function composition. The generated IGM configuration can always be mapped to a spatially periodic configuration (in a fictitious space) by means of a continuous mapping function. Then the TDF of the IGM configuration is obtained through the function composition of the mapping function with the TDF describing the cell configuration (of one period in the fictitious space). Under this setting, asymptotic analysis is carried out for deriving a homogenised formulation governing the mechanical behaviour of the resulting IGM configuration. Compared with the above-mentioned two methods, no post-processing procedure, such as solving the Ricci flow equation or conducting projecting operations, is needed in the AHTO plus framework. Nevertheless, other challenging issues arise over the following three aspects. Firstly, since the constituting cells (of parallelogram/parallelepiped shape) originate from a single generating cell, thus they must be self-similar to each other. This dramatically reduces the IGM variety in concern. For instance, spatial change in microstructural volume fraction is still disabled. Secondly, the number of design variables, which parameterise the controlling function, is relatively small under the AHTO plus framework. Thus the identification of an appropriate set of design variables becomes a critical issue, which has yet been properly addressed. Finally, the derived homogenisation formulation involves solving various PDE systems at individual macroscopic points, generating a huge computational burden. Hence finding an efficient way of solving these PDEs XueDC_arxiv2019 is also a key factor affecting the applicability of the method for IGM design.
The present article is aimed for addressing the first aforementioned issue, that is, to essentially improve the IGM variety under the AHTO plus framework. Compared to the original treatment in ZhuYC_JMPS2019 where local microstructures have to stem from a single generating cell, here we consider a continuous menu of generating cell configurations defined in a fictitious space, which can be systematically identified by the continuous level sets of a single function, say, , termed as a “structural generating function” in this article. Then a macroscopic continuous function defined in the actual space, say , is introduced to indicate which level set it takes around point . Here the function is termed as a “level-set indicator” function. Incorporating with the original function composition treatment, a TDF associated with an IGM configuration is thus introduced, and the “graded microstructural” behaviour is carried out by the combinative effects of the following two sets of operations: manipulation of a single cell and transition among cells. The action of cell manipulation involves 1) cell rescaling; 2) cell rotation; 3) cell distortion in an iso-volumetric manner. Mathematically, these operations over a single cell can be carried out by manipulating the local Jacobian matrix, which consists of the spatial derivatives of the mapping function. Such an idea for IGM representation originates from the use of level-set contours for modelling dislocation loops in crystalline materials XiangY_JMPS2009 ; ZhuYC_JMPS2015 ; Chapman_SIAP2016 . It bears certain similarity with the treatment by WangYQ_CMAME2017 , but essential differences exist. For instance, level sets of discrete values are used in WangYQ_CMAME2017 , but continuous level sets are considered here. Moreover, local rotation and distortion of the generating cells are not taken into account by WangYQ_CMAME2017 , resulting in pixel-wise microstructures only. This is, however, no longer a limiting issue if the present scheme is implemented. In fact, we will show that the TDF form given by Eq. (4) effectively offers a platform for unifying most existing treatments for IGM generation, and detailed discussion is carried out in Sec. 3.4. The present algorithm also sees its generalised applications in two important aspects. Firstly, it offers a simple treatment for continuously patching arbitrarily given cell configurations, which are termed as “seminal cell” configurations in the present article. Secondly, it breeds a general algorithm for generating an IGM surface.
The IGM configurations obtained here inherit all the attractive features from that generated based on the original AHTO plus framework, such as microstructural smoothness and simple forms in TDFs. Moreover, the corresponding asymptotic formulation for stress analysis is also available, and its accuracy is further shown against fine-scale simulation results. Under the present framework, one can monitor an IGM configuration by varying the mapping function , the level-set indicator function and the structural generating function . Mathematically, this can be fulfilled by a topology optimisation formulation with the associated sensitivity analysis results given. All these should pave way for adapting the AHTO plus framework for fast digital design of IGM configurations.
The rest of the article is arranged as follows. In Sec. 2, the original AHTO plus method is briefly revisited with the challenging issues identified. In Sec. 3, our algorithm for IGM generation from a continuous menu of cell configurations is introduced, and two of its generalised applications are discussed: 1) IGM generation continuously patching different seminal cells; 2) IGM generation on a given surface. The section concludes with discussion over the unification of the present methods for IGM representation under the present framework. In Sec. 4, the asymptotic formulation for analysing IGM mechanical behaviour is firstly introduced. It is followed by the corresponding topology optimisation formulation, along with the associated sensitivity analysis results. In Sec. 5, two-dimensional examples will be firstly presented with several key issues of the algorithm addressed. Three-dimensional cases will also be shown in brief in the end of Sec. 5. The article concludes with further discussion in Sec. 6.
2 IGM description under the original AHTO plus framework and challenging issues
2.1 IGM description under the original AHTO plus framework
A graded microstructural configuration filling a certain domain, as shown in the bottom-left panel of Fig. 1, normally involves at least two length scales:
a microscopic length scale on which the structural details are fully resolved, and a macroscopic length scale on which the trend in microstructural variation is observed. Here we characterise the “microscale” by a length parameter , and the “macroscale” by the domain size . Practically, we have
The mathematical representation of a structural configuration is usually realised through the introduction of a topology description function, denoted by here. For a porous region with the solid parts occupying a domain of , its associated TDF can be defined by
As for an IGM configuration, evaluating its associated TDF piece by piece on the fine scale is sometimes too time-consuming. For improvement, a treatment using function composition has been introduced Pantz_SIAMJCO2008 ; LiuC_JAM2017 , and further developped ZhuYC_JMPS2019 . The key is to make use of the smoothly-varying property of an IGM configuration.
As shown in Fig. 1, a macroscopically smooth mapping function is firstly introduced, so as to map an IGM configuration into a spatially periodic one in a fictitious space, whose constituting cells are uniform and of a same length scale . Then one may pick any cell (of rectangle/cuboid shape in the fictitious space) as a generating cell, and non-dimensionalise it into a region of unit dimensions, say with denoting the dimensionality number (as shown in the bottom-right panel of Fig. 1). Now a TDF , can be assigned to represent the structural configuration within the (non-dimensional) fictitious unit cell. Finally, the TDF characterising the actual porous region should be obtained by a composition of the (macroscopically smooth) mapping function and the TDF associated with the (fictitious) porous region in :
The TDF generated through function composition based on Eq. (3) has shown many attractive features. Firstly, by Eq. (3) is a continous function, because is smooth and is periodic over its input. Thus the resulting IGM configuration always takes a smoothly-varying profile. Secondly, Eq. (3) has a simple form, and this facilitates its implementation in a computer aided design (CAD) platform. Thirdly, it overcomes the limitation that the actual building blocks have to be of rectangle/cuboid shape. Actually, a variety of IMG configurations can be described by using Eq. (3), as demonstrated in ZhuYC_JMPS2019 . Fourthly, Eq. (3) poses low demand on digital resources needed for memorising it. In theory, the mapping function can be defined on a coarser grid, and the TDF which is of high resolution spans only a single cell. Finally, a homogenisation formulation used for analysing the mechanical behaviour of an IGM configuration generated by Eq. (3) can be derived through mathematically rigorous asymptotic analysis ZhuYC_JMPS2019 , as tends to 0. Thus the stress analysis of an IGM configuration can still be performed at a coarse-grained level, and this renders significantly bring down the requirement on computational resources, when compared to simulations directly conducted on a fine scale.
2.2 Challenging issues
Albeit many attractive properties as mentioned above, the TDF generated by function composition through Eq. (3) still sees severe limitation. For instance, it struggles to describe an IGM configuration, where the microstructural volume fraction varies in space. This is because the mapping function simply rescales, rotates and distorts a single generating cell configuration. Thus the volume fraction stays unchanged after being transformed by . For a detailed illustration over this point, one may refer to Appendix B in ZhuYC_JMPS2019 . Thus the present study is motivated: how to expand the coverage of the AHTO plus framework for describing more general IGM configurations.
3 Representation of general IGM configurations
In this section, we will first discuss on how to extend the function composition treatment described in Eq. (3) for a much more general use. Then two of its generalised applications will be discussed in depth. We first introduce an algorithm for generating IGM configurations smoothly patching various seminal cells. Then the formulation is further extended for representing graded microstructures decorated on a given manifold. Finally, we will show that the present formulation effectively offers a theoretical framework that unifies various ways of IGM representation in literature.
3.1 General formulation
For the purpose of greatly improve the variety of describable IGM configurations using the function composition treatment, the TDF given by Eq. (3) is now modified as
where , termed as the structural generating function in this article, is periodic over its input; is another smooth function in . A similar treatment has also been used in GuoX_CMAME2013 to account for possible boundary perturbations in topology optimization considering shape uncertainties.
, a level-set function is defined within the whole design domain, where only its zero-level-set contour has a geometric meaning of representing the structural boundaries. This idea is now generalised through Eq. (4), where all level-set values are geometrically meaningful. Practically, one may introduce a function , which is periodic over the fictitious non-dimensional domain of unit dimension measured by . As shown in Fig. 2, its corresponding level-set contour of height , given by , specifies a certain structural boundary. Since is continuous, the change in boundary profiles as varies should be smooth, too. Hence a continuous menu of generating cell patterns are obtained in the fictitious (non-dimensional) space measured by . Then for a point in the actual space, one needs to specify and . Here we let . Thus Eq. (4) indicates that the actual microstructural configuration around is generated by manipulating the generating cell identified by the level-set contour of of height .
The IGM configurations generated from Eq. (4) naturally inherit all the attractive features from the original treatment with function composition by Eq.(3). For instance, the resulting IGM configuration varies smoothly in space, because and are both continuous over their inputs and is periodic over its input.
The transition between generating cells is thus enabled, and so does the volume change in microstructures. To formulate the microstructural volume fraction, we first identify the correlation between the contour of of height and the volume fraction associated with the corresponding generating cells. Mathematically, this is given by
where is the Heaviside function. It is noted that the correlation formulated by Eq. (5) is already set as the continuous menu of generating cells are present, which proceeds the representation of the final IGM configurations. Besides, cell manipulation, involving rescaling, rotation and iso-volumetric distortion, carried out through , does not alter the corresponding local volume fraction. Thus the microstructural volume fraction around a spatial point solely depends on the value of the level-set indicator function , but is independent of . Therefore, the overall volume fraction of an IGM configuration characterised by Eq. (4) is formulated by
where is given by Eq. (5); “” denotes the volume occupied by the corresponding domain.
In theory, the significantly expanded IGM describability by Eq. (4) is demonstrated. In the following, we will further explore its applicability in a more practically meaningful context.
3.2 IGM generation from seminal cell patterns
A number of microstructural configurations have been proposed for achieving exceptional functions (e.g. Sigmund_APL1996, ; Frenzel_Science2017, ). One extended application of Eq. (4) is to smoothly transit between cells outputting different high-end performances. For instance, two seminal cell structures are present: a smiling face and a weeping face, as shown in Fig. 3.
Then the question is how to design an IGM configuration where the facial expression continuously turns from smiling to weeping. Note that this is by no means a trivial task. For instance, the two seminal structures are topologically different. Hence it is of difficulty to use a single set of parameters to consistently express both configurations. Besides, the heights of “cheeks” of the two faces are different. Using spatially discrete parameters to patch them, e.g. WangYQ_CMAME2017 ; ChengL_CMAME2019 , may result in persistent mismatch between neighbouring cells.
Suppose the two facial configurations shown in Fig. 3 are identified by the zero-level-set contours of two functions, and
, respectively. Then smooth connection between the two configurations should be obtained through interpolation between them. Hence we first construct a structural generating function by
where is the interpolation parameter. Compared to the structural generating function constructed in the previous subsection, defined by Eq. (7) contains an extra parameter . This time, the structural boundaries of the generating cell patterns are still identified by the implicit relationship given by . It is checked that when , becomes , corresponding to the structural boundary of the seminal cell given by Fig. 3(a). So on for the case when , corresponding to the weeping facial configuration identified by .
Although coming from Eq. (4), the TDF expression defined by Eq. (8) has a clear mathematical meaning. It is generated through direct interpolation between the TDFs identifying the two seminal cell configurations, which is controlled by the level-set indicator function .
Note that tiny structural evolution takes place within every single cell. Thus the cells constituting the final IGM configuration are normally not exactly the same as any of the generating cells. This can be visualised in the third face to the left in Fig. 4, where tear drops only appear beneath its left eye. This again addresses the key different feature of the present method, as compared with that in WangYQ_CMAME2017 , where the constituting cells have to coincide with one of the prototype cells and tiny mismatch persists across cell boundaries.
When the effect from the mapping function is taken into account, a wall of IGM configuration with mixed expressions, for instance, can be generated, as shown in Fig. 5.
As a mapping treatment is now involved, several faces get rotated or distorted in the output IGM configuration.
The treatment can be extended to smoothly connect a number of seminal configurations. Suppose we have seminal cell configurations, represented by the zero level sets of , , , , , respectively. Now we introduce a number set , and our goal is to generate a TDF such that it becomes when the interpolation parameter equals . This can be achieved in a similar sense as Lagrangian interpolation, and a general TDF expression is thus defined by
where the coefficients are formulated by
It can be checked that when , the solid region is identified by , which actually corresponds to the -th seminal cell configuration. Therefore, the proposed method can be used for generating IGM configurations from any given set of seminal cell patterns.
In fact, Eq. (9) may be extended to a more general form for describing a continuous menu of parameterised cell patterns, say identified by the generating function , where is a parameter. The corresponding TDF can be given by
By this mean, the pool of generating cells can be enriched by “mixing” seminal cells at a certain ratio given by .
3.3 Generation of surfaces with smoothly-varying microstructures
The method can be naturally generalised for deploying microstructures on a given manifold. A surface in three dimensions can be parameterised by two free variables, say, and , and it can thus be represented by
are the vector containing the coordinates of the points lying on the manifold;is a vector-valued function from to . Here a superscript “” is attached to a variable indicating its affiliation with a manifold. Now our target is to decorate it with smoothly-varying microstructures. For implementing the proposed method, we need to determine the TDFs of the generating cells and the mapping function .
This time, the (continuous) menu of generating cells should be planar structures in two-dimensional space, and suppose they are given by a periodic structural generating function , where the level-set of value identifies a (two-dimensional) cell configuration. Then we consider how to construct the mapping function so as to decorate the menu of cells on the manifold given by Eq. (12). It takes the following steps.
First, since is actually a map of
, an extra degree of freedom (in conjunction withand ), say , should be introduced. A natural way of thinking is to let it be in alignment with the surface normal , which is mathematically given by
where , and , denote the derivative of each entry of with respect to its -th entry (e.g. ), and the symbol “” denotes the cross product of two vectors. Then a map between and all points around the manifold is established:
where the subscript indicates the two inputs of the corresponding function. Then we choose , and . A mapping function from the fictitious space of to the actual space of is defined:
and the desired mapping function appearing in Eq. (4) should be given by . Normally, an actual manifold is assigned with a thickness, say , and this can be controlled by constraining the range of by letting . Finally, one may use the treatment of function composition to define the TDF for a given manifold decorated with microstructures (generated from , ) by
where is an arbitrary positive number; denotes the -th entry of . As from Eq. (16), the TDF values for all points falling outside the thickness of the manifold is negative, and no solid materials are distributed there according to Eq. (2).
One illustrative microstructural surface is shown in Fig. 6.
The generating cell has a circular hole in its middle, which is given by with . The manifold is chosen as a parabolic surface given by . Following the aforementioned steps, the map between the fictitious and the actual spaces is found to be
It is noted that for a same manifold, different parameterisations lead to disparate microstructural surfaces. For instance, both
correspond to a same hyperboloid. But they give rise to different manifold decorations, as shown in Fig 7.
3.4 Comparison among existing methods for IGM representation
The present treatment based on Eq. (4) effectively offers a unified theoretical framework, which encompasses most existing approaches for IGM representation. Actually, Eq. (4) (and its extended form by Eq. (8)) enables two sets of operations: selection of generating cell configurations by the value of , and the manipulation of the selected cell (involving cell rescaling, cell rotation and cell distortion in an iso-volumetric manner) controlled by the mapping function , or alternatively, by the corresponding Jacobian matrix given by
In theory, controls the operation exerted on a single cell, and controls the transition between different cells. In this viewpoint, most existing treatments for IGM generation can be considered as being carried out through varying and , and a summary is made in Table 1.
denotes an orthogonal matrix.
Here the reported IGM generation approaches in literature are in general classified into five groups:
For the first three groups of methods, the constituting cells must take a rectangle/cuboid shape. In a viewpoint of Eq. (4), each is only a function , for , , , or alternatively, the associated Jacobian matrix must be diagonal, i.e., , with a diagonal matrix. In contrast, with the projection algorithm Pantz_SIAMJCO2008 used in de-homogenisation methods, cell rotation is also permitted. In the context of Eq. (4), with an orthogonal matrix. And for the original AHTO plus method, the choice of is more arbitrary provided that . But in theory, a pool with should suffice for optimising the IGM mechanical behaviour, and this is well demonstrated by Groen et al. Groen_IJNME2018 . It should also be noted that, the edges of the constituting cells, in categories 4) and 5), are slightly curved in the actual space so as to accommodate smooth connection across cell boundaries.
The IGM smoothness is also controlled by the continuity in function appearing in Eq. (4). For 1) classical homogenisation-based approaches, each pixel is individually assigned with a microstructure. Thus is effectively piece-wisely defined over individual pixels. Therefore, non-smoothness emerges in the resulting IGM configurations. For category 2) MPCs, the microstructure is generated through gradually varying the parameter values of its building blocks, but stays the same in every single cell. Hence should be a step function, where the gaps between neighbouring steps are tiny. Consequently, the resulting IGM configurations appear gradually-varying, but slightly non-smooth connection still emerges across cell boundaries. Similar case arises for the category 3) MLSMs, where the microstructure is generated through gradually controlling the outputs of a single level-set function. In contrast, the de-homogenisation methods enable smooth transition across cell boundaries in the actual space, which should correspond to a continuous function . As for the original AHTO plus method, no transition between cells is permitted, giving rise to a constant-value function of . But this limitation is removed in the present article by Eq. (4).
Therefore, most existing methods for IGM generation should find their corresponding formulation through Eq. (4) by selecting appropriate and .
4 Asymptotic formulation towards optimising IGM compliance
In this section, we will first consider how to formulate the mechanical behaviour of IGM configurations obtained based on Eq. (4). Then an optimisation problem for IGM compliance is established, with the corresponding sensitivity analysis being discussed at the end of the section.
4.1 IGM mechanical behaviour
Under the original AHTO plus framework ZhuYC_JMPS2019 , the multiscale mechanical response of a loaded IGM configuration can be predicted through asymptotic analysis. This is carried out by approximating the (multiscale) displacement field by a series expanded in terms of the small parameter defined by Eq. (1). Each term in the asymptotic expansion is a function containing two sets of spatial variables, say and . The variable captures the homogenised behaviour as if the IGM configuration is envisaged as a continuum. The (non-dimensional) variable captures the large displacement gradient accommodating the local microstructures. Through asymptotic analysis, the originally multiple-scale governing equations are solved on different scales individually.
In the following, the key equations constituting the original AHTO plus formulation ZhuYC_JMPS2019 are outlined. A weak form of the equilibrium equation is set up for the homogenised displacement field , given by
|where denotes the effective elasticity tensor of the homogenised continuum; is a virtual displacement field belonging to the functional space of ; denotes the boundary on where displacement boundary condition is imposed; is the boundary on where a surface traction of density is applied.|
Meanwhile, a set of characteristic functions, , , , , , are also needed, and they are obtained by solving the following set of PDEs:
subjected to periodic boundary conditions over , where is the Jacobian matrix defined by Eq. (19); is a piecewisely constant function: it equals the corresponding elasticity tensor at a solid point, and vanishes elsewhere in . It is noted that the spatial variables of are , while simply serves as a parameter. But this means one has to solve Eq. (20b) individually at every single macroscopic point in the actual space.
where is recalled to be the solid part in .
Such a scale-separation treatment roots in the fact that an IGM configuration from the original AHTO plus framework becomes periodic in the fictitious space measured by , or equivalently by . For the present case, however, the generated IGM configurations are not exactly periodic in the non-dimensional fictitious space measured by . Here we still adopt the same asymptotic formulae as given by Eqs. (20a) - (20c). The accuracy of the treatment here will be demonstrated through a comparison with the corresponding fine-scale simulation results in Sec. 5.
4.2 Topology optimisation
In this subsection, based on the results above, a mathematical formulation for optimising IGM compliance can be set up as follows
where denotes the functional space of continuous functions; is a certain functional space formed by periodic functions defined over the non-dimensional unit cell ; is the homogenised structural compliance; is the upper limit of the IGM solid volume. The last inequality in Eq. (21), as from Eq. (6), formulates the volume fraction constraint.
Eq. (21) establishes a scheme for minimising IGM compliance. Compared to the original AHTO plus framework, it exhibits several advantageous features. Firstly, the available IGM configurations are far more abundant. This is because one can not only conduct IGM design by rotating and distorting a generating cell through the mapping function , but also select among generating cells through the level-set indicator function . Secondly, under the original AHTO plus framework, one can not increase the volume of solid materials to strengthen the local microstructure, because the volume fraction of all constituting cells must be identical to each other. But this limitation is fully removed now. It will be shown in Sec. 5 that enabling transition among different transition cells should bring down the compliance value of an optimised structure by almost a half.
4.3 Gradient-based sensitivity analysis
For speeding up the optimisation process based on Eq. (21), its associated sensitivity analysis has been performed. Suppose that the vector quantifies all the design variables. Through using the idea of the adjoint method LiuST_SMO2002 , the derivative of the IGM compliance with respect to the -th design variable can be formulated by
Thus one simply needs to compute the derivative of the elasticity tensor with respect to all the design variables. Here the design variables are classified into two groups: the macroscopic design variables parameterising the mapping function and the level-set indicator function ; the microscopic design variables parameterising the structural generating function .
For sensitivity analysis with respect to typical macroscale variables, say, , we follow XueDC_arxiv2019 to write down
As for the sensitivity with respect to the microscale variables, say, , we also follow XueDC_arxiv2019 to write down
where , , , , , is the Kronecker delta.
5 Numerical results and analysis
In this section, numerical results based on the proposed scheme are presented. Two-dimensional examples will first be examined in depth, so as to reveal the key features of the present method. Then three-dimensional numerical examples will also be briefly studied at the end of the section.
5.1 General problem setting
For all the examples examined in this section, linearly elastic material response is assumed with the isotropic parameters chosen as follows: the (nondimensional) Young’s modulus and Poisson’s ratio . For two-dimensional examples, the specimen’s (non-dimensional) thickness in the third dimension is fixed to be 1, and scenarios of plane stress are considered. The upper limit of IGM volume fraction is set to be 0.3 unless otherwise specified.
Here the continuous menu of generating cell configurations, once chosen, are set fixed during the optimisation processes. The reason for this action is that IGM generation based on Eq. (4) should output abundant configurations. This is in particular reasonable, when IGM configurations are produced through continuously patching seminal cell structures. Thus the evolution of the design variables parameterising the structural generating function is temporarily suspended here. As for the parameterisation of the mapping function , we follow XueDC_arxiv2019 to assume
where symmetric conditions and are imposed.
As for the level-set indicator function , we let
where for , , , , is required.
Upon this set-up, the design variables are collected by
Hence we have 24 design variables for two-dimensional problems, and 67 design variables for three-dimensional cases. Here the length-scale characterising parameter is set to be 0.05, and the method of moving asymptotes (MMA) Svanberg_IJNME1987 are employed as the optimiser in this work.
Note that the homogenisation formulation for stress analysis involves individually solving a cell problem defined by Eq. (20b) at every single macroscopic point . This significantly brings up the computational cost. For overcoming this issue, a zoning scheme XueDC_arxiv2019 is employed, where the design domain is partitioned into several subdomains. In each subdomain, a point is selected and the effective elasticity tensor at that point obtained through Eq. (20b) is used to represent that of the whole subdomain. Here all design domains are divided into rectangles/cuboids, and the center of each rectangle/cuboild is chosen as the representative point. It has been shown XueDC_arxiv2019 that, when a subdomain contains no more than 25 finite elements (for two-dimensional cases), the deviation from the corresponding fine-scale results is less than 0.04. We will resume using this criterion for conducting simulations for this work. As for finite element analysis, four-node quadrilateral elements are used in all two-dimensional examples while eight-node hexahedral brick elements are adopted in all three-dimensional examples.
5.2 Two-dimensional cases
5.2.1 Laminar IGM configurations
We start with an example of a short beam fixed on its right side, as shown in Fig. 8.
A uniformly distributed compression of magnitude 0.1 is applied over its top side. The design domain is of size , and it is divided into identical square-shape subdomains. A mesh is used for solving the homogenised problem.
The continuous menu of generating cells are selected as a set of “X”-shape configurations, which are controlled by the thickness of the solid region. Mathematically, all generating cells can be represented by the multiple level-set contours of a single function:
and the structural boundaries of several generating cell configurations are shown in the left panel of Fig. 9.
Then according to Eq. (4), the TDF associated with the corresponding IGM configuration is thus given by
Note that the range of by Eq. (30) defines the bounds for , that is, . In practice, if an output of Eq. (28) exceeds , we require it equal , and so does an output falling below . Under the present setting, the overall IGM volume fraction, according to Eq. (6), is formulated by
The optimised result is then shown in the right panel of Fig. 9, and the IGM volume fraction is roughly 29.98%. The structural compliance based on the homogenisation formulation is computed to be 126.60, which delivers a 2.0% deviation from the corresponding fine-scale simulation result (124.24). Fig. 10 shows the evolution of both the compliance value and the corresponding volume fraction during optimisation. Compared to the initial structure composed of a spatially periodic microstructure, the structural compliance value drops more than 90%. It is noted that the start of evolution of the level-set indicator function is 30 steps later than that of the mapping function . This is because for minimising structural compliance, rotation and distortion of a single cell seems more effective than simply transiting among generating cells. This will be discussed with greater details later.
Compared to the spatially periodic configuration, the optimised one arranges the constituting fibres such that the external load on the top side of the design domain smoothly transited to the fixed boundary. Besides, the material distribution is denser at the top-right corner, where more materials are needed for accommodating the relatively high local bending moment.
Eq. (4) effectively generates an IGM configuration combinatively by two means. One is through distorting and rotating a single cell by the mapping function , and the other is through transition among generating cells by the level-set indicator . Their individual roles in compliance optimisation are examined with the results summarised in Fig. 11 and Table 2.
|Periodic||Cell manipulation||Cell transition||Both|
|structure||enabled only||enabled only||enabled|
|design||Fig. 11(a)||Fig. 11(b)||Fig. 11(c)||Fig. 11(d)|
It is shown that optimisation by simply permitting rotation and distortion of a single cell reduces the compliance value to 1/8 of that of the periodic structure shown in Fig. 11(a). This is in contrast to a drop of 1/2 if only transition between generating cells is enabled. This means arranging laminar orientations in align with the local directions of principal stress provides a more effective mean than simply redistributing materials by means of cell transition. However, compared to the optimised configuration in Fig. 11(b), further enabling transition among the generating cells results in an extra 50% drop in the optimised compliance. This demonstrates the necessities of involving transition among different cell configurations during IGM generation.
5.2.2 IGM generated through interpolating between seminal cells
In this case, the continuous menu of generating cells are obtained by interpolating between two seminal cells, which are chosen as a full solid square with a diminishing circular hole at its centre and a unit square with a hyperelliptic hole in its middle. Mathematically, a full-solid configuration with a diminishing hole can be identified by the zero level-set contour of
The cell with a hyperelliptic hole at its centre is indicated by the zero level-set contour of
Thus based on Eq. (7), one introduces a structural generating function by interpolating and controlled by a parameter :
such that the implicit relation corresponds to the structural boundary of one generating cell configuration in the fictitious space measured by . The structural boundaries of several generating cells are shown in the left panel of Fig. 12.
It can be checked that when is near 0, casts a dominant weight over . In this scenario, identifies a single point, and this corresponds to a square of full solid. As keeps growing, the hole shape looks more like a mixture of a circle and a hyper-ellipse. When , becomes , corresponding to a cell with a hyperelliptic hole. Note that the transition in the hole shape is carried out smoothly, since given by Eq. (35) is a continuous function (in ).
The volume fraction of the generating cells, obtained based on Eq. (35), spans from 100% () to 3.81% (). In contrast to the previous case, no explicit expressions linking a generating cell with its corresponding volume fraction are available. Nevertheless, one may pre-calculate a data set correlating the two quantities, and fit for a simple function approximating the volume fraction of generating cells against the parameter by
Finally, by introducing the mapping function by , and the level-set indicator function , we express the TDF of the IGM configurations for topology optimisation:
The corresponding optimised IGM design is shown in the right panel of Fig. 12 with the volume fraction being 29.94%. The homogenised result for structural compliance is 124.40, which delivers a deviation of 2.1% from the fine-scale computation (121.73). Fig. 13 shows the evolution in the values of both the structural compliance and the corresponding volume fraction during optimisation.
Again, the individual effects of cell manipulation (through ) and cell transition (through ) are examined for this case. The results are summarised in Fig. 14 and Table 3, and a conclusion similar as that in Sec. 5.2.1 can be drawn.
|Periodic||Cell manipulation||Cell transition||Both|
|structure||enabled only||enabled only||enabled|
|design||Fig. 14(a)||Fig. 14(b)||Fig. 14(c)||Fig. 14(d)|
Several aspects of similarities are read. Firstly, microstructural conduits transiting the load to the fixed boundaries are seen in both optimised configurations. Secondly, the trends in volume change are the same. Thirdly, both optimised configurations output similar compliance values. But the design given by Fig. 12 tends to reveal how a trade-off is made during optimisation. It has been widely agreed that Allaire_book2002 , for compliance optimisation, laminate microstructures should deliver a near-optimal performance, provided that the laminate orientation properly align with the desired principal stresses. Compared to circular holes, configurations with a hyper-elliptic hole are morphologically closer to laminate structures. Fig. 12 indicates laminate microstructures are more desirable away from the fixed boundary and the loaded edge, while denser material distribution close to the fixed and loaded boundaries seems to be more effective for compliance optimisation.
5.3 Three dimensional examples
In this subsection, topology optimisation results of three-dimensional cases are briefly presented, and further investigations are highly desired for its improvement. Here we consider a short beam of size with its right side fixed, as shown in Fig. 16.
On its upper surface, a uniformly distributed load is applied. The homogenised design domain is partitioned into identical cubic subdomains. For FE analysis, a mesh of is adopted.
The first example concerns a continuous menu of generating cells, whose representatives are shown in Fig. 17.
Every cell consists of three identical pillars, with each pillar crossing the other two in their middles. The cross-sectional shape continuously changes from a unit circle to an almost diminishing hyper-ellipse, as the volume fraction drops. Mathematically, the solid parts of these generating cells can be described using a parameter :
Starting from an initially periodic structure shown in Fig. 18(a), the optimised IGM configuration is given in Fig. 18(b). In a side view, the three-dimensional optimised result bears a strong similarity with the corresponding two-dimensional configuration shown in Fig. 12, that is, the microstructural cells arrange themselves for more effectively transiting the load to the fixed boundary.
Note that for better visualising the results, cell rotation is only permitted within the planes perpendicular to the -direction of the design domain.
Another set of generating cells are also concerned with some of their representatives shown in Fig. 19.
This time, every generating cell is obtained by digging holes from a solid cube. The void region can be envisaged as the union of three cylinders, one crossing the other two. The cross-section of each cylinder smoothly transits from an almost diminishing circle to a hyperellipse of size 1, as the volume fraction drops. Mathematically, the solid regions within the generating cells can be described by
Starting from an initially periodic configuration in Fig. 20(a), the optimised configuration is shown in Fig. 20(b), where the microstructural orientation obtained is roughly the same as that in Fig. 18(b).
6 Conclusion and discussion
The present article proposes an algorithm for generating infill graded microstructural configurations from a continuous menu of simple generating cells. With the present scheme, the IGM profile is controlled through the introduction of two functions: a level-set indicator function identifying which cell from the menu is of concern around a macroscopic point ; and a mapping function indicating how the cell is rotated and distorted locally. With these two macroscopically continuous functions, IGM configurations can be generated based on TDF expressed by Eq. (4). Eq. (4) effectively offers a unified expression where most existing methods for IGM generation find their traces, with detailed information summarised in Table 1. The present scheme also sees its generalised applications in two aspects: 1) IGM generation through interpolation among a set of desired seminal cell configurations; 2) IGM configurations decorated on a given manifold. The IGM mechanical behaviour can then be predicted using asymptotic analysis, and the accuracy of the homogenised formulation is further demonstrated through comparisons with the corresponding fine-scale simulation results. The sensitivity analysis is also studied, and a (homogenised) topology optimisation formulation for IGM design is also proposed. Numerical examples for two-dimensional cases are systematically analysed, and three-dimensional results are also briefly shown.
The present approach brings about considerable expansions in IGM describability under a topology optimisation framework underpinned by asymptotic analysis, which is termed as an AHTO plus framework here. Compared to its original version ZhuYC_JMPS2019 , the IGM configurations obtained from the present generation scheme inherit all the attractive IGM features, such as smooth connection and simple form in TDF expression. Moreover, IGM configurations with other practically useful features, such as smooth change of volume fraction in space, are also enabled by the present algorithm. As discussed in the beginning of this article, the original AHTO plus method ZhuYC_JMPS2019 needs essential improvements mainly over three issues: 1) to substantially raise the computational efficiency when implementing the homogenised formulation for stress analysis; 2) to overcome the limitation where all constituting cells are restricted to be self-similar to each other; 3) to find more appropriate basis functions governing the IGM generation through Eq. (4). The present article aims for addressing the second issue mentioned above, while the first issue is partly resolved by XueDC_arxiv2019 and the third issue still needs further investigations. This also points to one promising direction for future studies.
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), Program for Changjiang Scholars, Innovative Research Team in University (PCSIRT) and 111 project (B14013) are gratefully acknowledged.
- (1) Y. C. Zhu, S. S. Li, Z. L. Du, C. Liu, X. Guo, and W. S. Zhang. A novel asymptotic-analysis-based homogenisation approach towards fast design of infill graded microstructures. J. Mech. Phys. Solids, 124:612–633, 2019.
- (2) R. Lakes. Materials with structural hierarchy. Nature, 361(6412):511–515, 1993.
- (3) H. W. Dong, S. D. Zhao, Y. S. Wang, and C. Z. Zhang. Topology optimization of anisotropic broadband double-negative elastic metamaterials. J. Mech. Phys. Solids, 105:54–80, 2017.
- (4) O. Sigmund and S. Torquato. Composites with extremal thermal expansion coefficients. Appl. Phys. Lett., 69(21):3203–3205, 1996.
- (5) M. S. Kushwaha, P. Halevi, L. Dobrzynski, and B. Djafari-Rouhani. Acoustic band structure of periodic elastic composites. Phys. Rev. Lett., 71(13):2022–2025, 1993.
- (6) C. Liu, Z. L. Du, Z. Sun, H. J. Gao, and X. Guo. Frequency-preserved acoustic diode model with high forward-power-transmission rate. Phys. Rev. Appl., 3(6):064014, 2015.
- (7) N. Aage, E. Andreassen, B. S. Lazarov, and O. Sigmund. Giga-voxel computational morphogenesis for structural design. Nature, 550(7674):84–86, 2017.
- (8) A. Bensoussan, J. L. Lions, and G. Papanicolaou. Asymptotic analysis for periodic structures. Studies in mathematics and its applications. North Holland, Amsterdam, 1978.
- (9) M. P. Bendsoe and N. Kikuchi. Generating optimal topologies in structural design using a homogenization method. Comput. Methods Appl. Mech. Engrg., 71(2):197–224, 1988.
- (10) M. P. Bendsoe. Optimal shape design as a material distribution problem. Struct. Optim., 1(4):193–202, 1989.
- (11) O. Sigmund. Materials with prescribed constitutive parameters: An inverse homogenization problem. Int. J. Solids Struct., 31(17):2313–2329, 1994.
- (12) H. Rodrigues, J. M. Guedes, and M. P. Bendsoe. Hierarchical optimization of material and structure. Struct. Multidiscip. Optim., 24(1):1–10, 2002.
- (13) L. Liu, J. Yan, and G. D. Cheng. Optimum structure with homogeneous optimum truss-like material. Comput. Struct., 86(13):1417–1425, 2008.
- (14) B. Niu, J. Yan, and G. D. Cheng. Optimum structure with homogeneous optimum cellular material for maximum fundamental frequency. Struct. Multidiscip. Optim., 39(2):115–132, 2009.
- (15) J. D. Deng, J. Yan, and G. D. Cheng. Multi-objective concurrent topology optimization of thermoelastic structures composed of homogeneous porous material. Struct. Multidiscip. Optim., 47(4):583–597, 2013.
- (16) J. Yan, X. Guo, and G. D. Cheng. Multi-scale concurrent material and structural design under mechanical and thermal loads. Comput. Mech., 57(3):437–446, 2016.
- (17) J. D. Deng and W. Chen. Concurrent topology optimization of multiscale structures with multiple porous materials under random field loading uncertainty. Struct. Multidiscip. Optim., 56(1):1–19, 2017.
- (18) S. W. Zhou and Q. Li. Design of graded two-phase microstructures for tailored elasticity gradients. J. Mater. Sci., 43(15):5157, 2008.
- (19) A. Radman, X. Huang, and Y. M. Xie. Topology optimization of functionally graded cellular materials. J. Mater. Sci., 48(4):1503–1510, 2013.
- (20) A. Radman, X. Huang, and Y. M. Xie. Maximizing stiffness of functionally graded materials with prescribed variation of thermal conductivity. Comp. Mater. Sci., 82:457–463, 2014.
- (21) L. Cheng, J. X. Bai, and A. C. To. Functionally graded lattice structure topology optimization for the design of additive manufactured components with stress constraints. Comput. Methods Appl. Mech. Engrg., 344:334–359, 2019.
- (22) Q. H. Li, R. Xu, J. Liu, S. T. Liu, and S. Zhang. Topology optimization design of multi-scale structures with alterable microstructural length-width ratios. Comp. Struct., 230:111454, 2019.
- (23) Y. Q. Wang, F. F. Chen, and M. Y. Wang. Concurrent design with connectable graded microstructures. Comput. Methods Appl. Mech. Engrg., 317:84–101, 2017.
- (24) J. Gao, H. Li, L. Gao, and M. Xiao. Topological shape optimization of 3d micro-structured materials using energy-based homogenization method. Adv. Engrg. Softw., 116:89–102, 2018.
- (25) Y. Zhang, M. Xiao, H. Li, L. Gao, and S. Chu. Multiscale concurrent topology optimization for cellular structures with multiple microstructures based on ordered simp interpolation. Comput. Mater. Sci., 155(4):74–91, 2018.
- (26) Y. Zhang, H. Li, M. Xiao, L. Gao, S. Chu, and J. Zhang. Concurrent topology optimization for cellular structures with nonuniform microstructures based on the kriging metamodel. Struct. Multidiscip. Optim., 59(4):1273–1299, 2019.
- (27) J. Alexandersen and B. S. Lazarov. Topology optimisation of manufacturable microstructural details without length scale separation using a spectral coarse basis preconditioner. Comput. Methods Appl. Mech. Engrg., 290:156–182, 2015.
- (28) C. Liu, Z. L. Du, W. S. Zhang, Y. C. Zhu, and X. Guo. Additive manufacturing-oriented design of graded lattice structures through explicit topology optimization. ASME J. Appl. Mech., 84(8):081008–1—081008–12, 2017.
- (29) P. Vogiatzis, M. Ma, S. K. Chen, and X. F. D. Gu. Computational design and additive manufacturing of periodic conformal metasurfaces by synthesizing topology optimization with conformal mapping. Comput. Methods Appl. Mech. Engrg., 328:477–497, 2018.
- (30) J. P. Groen and O. Sigmund. Homogenization-based topology optimization for high-resolution manufacturable microstructures. Internat. J. Numer. Methods Engrg., 113(8):1148–1163, 2018.
- (31) J. P. Groen, J. Wu, and O. Sigmund. Homogenization-based stiffness optimization and projection of 2d coated structures with orthotropic infill. Comput. Methods Appl. Mech. Engrg., 349:722–742, 2019.
- (32) G. Allaire, P. Geoffroy-Donders, and O. Pantz. Topology optimization of modulated and oriented periodic microstructures by the homogenization method. Computers and Mathematics with Applications, 78(7):2197–2229, 2019.
- (33) J. Wu, W. M. Wang, and X. F. Gao. Design and optimization of conforming lattice structures. IEEE Transactions on Visualization and Computer Graphics, pages 1–1, 2019.
- (34) P. Geoffroy-Donders, G. Allaire, and O. Pantz. 3-d topology optimization of modulated and oriented periodic microstructures by the homogenization method. J. Comp. Phys., 401:108994, 2020.
- (35) O. Pantz and K. Trabelsi. A post-treatment of the homogenization method for shape optimization. SIAM J. Control Optim., 47(3):1380–1398, 2008.
- (36) D. C. Xue, Y. C. Zhu, S. S. Li, C. Liu, W. S. Zhang, and X. Guo. On speeding up an asymptotic-analysis-based homogenisation scheme for designing gradient porous structured materials using a zoning strategy. arXiv:1909.11182, 2019.
- (37) Y. Xiang. Continuum approximation of the peach-koehler force on dislocations in a slip plane. J. Mech. Phys. Solids, 57:728–743, 2009.
- (38) Y. C. Zhu and Y. Xiang. A continuum model for dislocation dynamics in three dimensions using the dislocation density potential functions and its application to micro-pillars. J. Mech. Phys. Solids, 84:230–253, 2015.
- (39) S. J. Chapman, Y. Xiang, and Y. C. Zhu. Homogenization of a row of dislocation dipoles from discrete dislocation dynamics. SIAM J. Appl. Math., 76:750–775, 2016.
- (40) X. Guo, W. S. Zhang, and L. Zhang. Robust structural topology optimization considering boundary uncertainties. Comput. Methods Appl. Mech. Engrg., 253:356–368, 2013.
- (41) M. Y. Wang, X. M. Wang, and D. M. Guo. A level set method for structural topology optimization. Comput. Methods Appl. Mech. Engrg., 192(1):227–246, 2003.
- (42) G. Allaire, F. Jouve, and A. M. Toader. Structural optimization using sensitivity analysis and a level-set method. J. Comp. Phys., 194(1):363–393, 2004.
- (43) Tobias Frenzel, Muamer Kadic, and Martin Wegener. Three-dimensional mechanical metamaterials with a twist. Science, 358(6366):1072, 2017.
- (44) P. G. Coelho, P. R. Fernandes, J. M. Guedes, and H. C. Rodrigues. A hierarchical model for concurrent material and topology optimisation of three-dimensional structures. Struct. Multidiscip. Optim., 35(2):107–115, 2008.
- (45) C. Wang, J. H. Zhu, W. H. Zhang, S. Y. Li, and J. Kong. Concurrent topology optimization design of structures and non-uniform parameterized lattice microstructures. Struct Multidiscip O, 58(1):35–50, 2018.
- (46) S. T. Liu, G. D. Cheng, Y. Gu, and X. G. Zheng. Mapping method for sensitivity analysis of composite material property. Struct. Multidiscip. Optim., 24(3):212–217, 2002.
- (47) K. Svanberg. The method of moving asymptotes—a new method for structural optimization. Internat. J. Numer. Methods Engrg., 24:359–373, 1987.
- (48) G. Allaire. Shape Optimization by the Homogenization Method, volume 146. Springer Science and Business Media, 2002.