On speeding up an asymptotic-analysis-based homogenisation scheme for designing gradient porous structured materials using a zoning strategy

09/25/2019 ∙ by Dingchuan Xue, et al. ∙ Dalian University of Technology 0

Gradient porous structured materials possess significant potential of being applied in many engineering fields. To help accelerate the design of infill graded microstructures, a novel asymptotic homogenisation topology optimisation method has been proposed by Zhu et al.(2019), aiming for 1) significantly enriching the pool of representable graded microstructures; 2) deriving an homogenised formulation underpinning its stress analysis in consistency with fine-scale results. But the work is severely confined from being widely applied mainly due to the following two reasons. Firstly, to circumvent the extremely time-consuming computation of the microscopic cell problems at every macroscopic point, its linearised form has to be used in numerical implementation, and this significantly reduces the design freedom. Secondly, lacking of an associated formulation of sensitive analysis, genetic algorithm was chosen for optimisation, which inevitably decreases the computational efficiency. To address these bottleneck challenging issues, a high-efficiency approach for the analysis and design of structures filled with quasi-periodic graded microstructures is proposed. A zoning scheme in favour of parallel computation is introduced, so as to make the best use of the new asymptotic homogenisation results. The proposed algorithm is validated through a comparison with existing two-dimensional benchmark cases and fine-scale brute force simulation results. The present approach is also applied for the topology optimisation of three-dimensional graded microstructures, which have barely been examined in literature, possibly because of the high computational cost, and the significantly high computational efficiency of the proposed scheme is clearly demonstrated by the provided examples.



There are no comments yet.


page 12

page 13

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Functional equipments and devices infilled with well-designed 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 efficiency-versus-accuracy dilemma, is to develop multiscale approaches in a viewpoint of scale separation. The basic idea lying behind the scale-separation 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 multi-functional 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 high-performance 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 scale-separation 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 th-order 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 re-build 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 time-consuming part of the AHTO plus calculation, only depends on macroscopic positions. In the present work, a divide-and-conquer 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 three-dimensional GMC optimisation problem which has rarely been examined in literature, possibly because of its high computational cost, it takes an eight-processor 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 manifold-like 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 so-called divide-and-conquer 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 gradient-based 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 manifold-like 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 quasi-periodic graded microstructures

Graded microstructural configurations generally possess a multiple-length-scale 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,


In this article, we term the length scale characterised by the “microscale”, and term the length scale characterised by the “macroscale”.

Figure 1: A configuration filled with graded microstructures and the corresponding mapping treatment.

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 bottom-right 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 quasi-periodic 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 bottom-right 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 bottom-left 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.,


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



denotes vector of the displacement field,

is vector of the body force density and

is the fourth-order 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:


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


where denotes the porous solid region within the unit generating cell and denotes its measure; (, , , ) is the Jacobean matrix associated with the mapping function


measuring the degree of spatial variance of the GMC; is a third-order tensor satisfying a microscale cell problem


equipped with periodic boundary conditions, where is given by


Once the homogenised displacement field is obtained, the homogenised stress field can be computed by


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.


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 above-described 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 third-order 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 spatially-periodic configuration. As a result, the design freedoms reduce significantly. For example, the linearised AHTO plus method simply outputs an optimised two-dimensional 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 resource-demanding 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 high-efficiency solution of the corresponding optimisation problems.

3 A parallel divide-and-conquer strategy

In this section, a parallel divide-and-conquer 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 time-consuming 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 divide-and-conquer strategy is to divide the macroscopic region into subdomains as shown in Fig. 2.

Figure 2: Domain decomposition.

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 third-order tensor , , , , , , can be written as


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 subdomain-wise form)


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


The above-listed 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 two-dimensional 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


where is a set of basis functions and is a set of weight parameters. Then the Jacobean matrix can be written by


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


for and is thus expressed by


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 “divide-and-conquer” optimisation scheme based on the AHTO plus formulation can be illustrated by the flowchart given in Fig. 3.

Figure 3: Flowchart of the divided-and-conquer scheme.

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 Gradient-based sensitivity analysis

To further speed up the optimisation process, the gradient-based 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


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


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







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


By adopting a similar treatment as in Liu et al. (2002), the sensitivity with respect to the microscale variables can be finally expressed by


where is the second-order identity tensor.

In this work, the so-called MMC-based 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 hyper-elliptical/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 divide-and-conquer scheme are presented. Firstly, several benchmarked two-dimensional examples will be considered, so as to demonstrate the effectiveness of the zoning strategy with parallel computing treatments. Then several three-dimensional 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 two-dimensional 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, four-node quadrilateral elements are used in all two-dimensional examples while eight-node hexahedral brick elements are adopted in all three-dimensional examples.

In addition, we adopt degree-3 polynomials formulated by Eq. (17) as the basis functions in all examples, and the method of moving asymptotes (MMA) Svanberg (1987) are employed as the optimiser in this work.

5.2 Two-dimensional 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. 


Figure 4: A short beam subjected to uniformly distributed force.

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%.

Figure 5: Convergence curve of the GMC compliance.

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 fine-scale 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 thin-boundary 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.

Figure 6: The optimised design obtained with X-shape material distribution in the unit generating cell.

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 .

Figure 7: A short beam loaded at the middle point on its left side. The upper-half region chosen as the design domain with symmetric condition imposed across the middle line, and the design domain is analysed and it is divided into 82 identical subdomains, each of which is a square.

The upper-half 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 building-block 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.

Figure 8: An optimised GMC design for the short beam example. The building block microstructure in the unit generating cell is of X-shape as shown in the left panel.

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.

Figure 9: Left: the optimised GMC obtained with the present method; Right: the optimised GMC using the method in Groen et al. (2019).

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.

Figure 10: A plot of the energy density distribution (in logarithm term) and the maximum principal stress component where the direction is denoted by an arrow.

The GMC compliance by the present homogenisation-based 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.

It is also worth noting that compared to the treatment in Groen et al. (2019) works for microstructures of laminate type, but that is not a restriction for the use of the present approach, and more details and examples can be seen in § 5.3.3.

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 square-shape 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.

Figure 11: The values of the objective function (i.e. compliance) and optimised configurations obtained with different numbers of subdomains.

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
Table 1: Solution time decomposition in one iteration step (2D short beam example)

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.

Figure 12: Non-dimensional solution time for each step with different numbers of workers.

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 three-dimensional cases and this will be shown in Sec 5.4.2.

Figure 13: Optimised design with different choices of unit generating cells.

5.3.3 Diversity of the describable GMC

By assuming that the microstructural cell takes an nth-order laminate configuration, the projection-based approach Groen and Sigmund (2018) also provides an effective way to design GMC (in two-dimensional space). The AHTO plus method, in contrast, does not pose any restriction in the form of involved microstructures. For the same short-beam 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 non-laminated shape. The AHTO plus method has also shown its potential of doing optimisation considering other physical performances of GMC.

5.4 Three-dimensional cases

In this subsection, the stiffness optimisation of three-dimensional GMC is studied. We will start with a quasi-two-dimensional example where the microstructural change only takes place in two directions, to validate the three-dimensional algorithm. Then full three-dimensional cases are studied, and the effectiveness of using parallel computing is also demonstrated.

5.4.1 A quasi-two-dimensional example

Figure 14: A quasi-two-dimensional short beam 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 two-dimensional 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 quasi-two-dimensional 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.

Figure 15: The optimised design obtained with a X-shape generating cell.

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 three-dimensional example

Now we consider a more general case where the building-block 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.

Figure 16: The optimised design obtained with an Orthogonal-shape generating cell and 4 subdomains.

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).

Figure 17: Deformation of unit generating cell without imposing the constraint on the determinant of under the mapping function . Note that under this circumstance, the size of the constituting cells may go beyond the validity range of the AHTO plus formulation.

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 third-oder 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
Table 2: Solution time decomposition in one iteration step (3D short beam example)

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 manifold-like 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.

Figure 18: Virtual mechanically driven complex 3D structural design.

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 manifold-like 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 divide-and-conquer 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 gradient-based 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 sub-function 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.

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 Chang-jiang 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.


  • N. Aage, E. Andreassen, B. S. Lazarov, and O. Sigmund (2017) Giga-voxel computational morphogenesis for structural design. Nature 550 (7674), pp. 84–86. Cited by: §1.
  • J. Alexandersen and B. S. Lazarov (2015) 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.
  • G. Allaire, F. Jouve, and A. M. Toader (2004) Structural optimization using sensitivity analysis and a level-set method. J. Comp. Phys. 194 (1), pp. 363–393. Cited by: §4.
  • M. P. Bendsoe and N. Kikuchi (1988) Generating optimal topologies in structural design using a homogenization method. Comput. Methods Appl. Mech. Engrg. 71 (2), pp. 197–224. Cited by: §1.
  • M. P. Bendsoe (1989) Optimal shape design as a material distribution problem. Struct. Optim. 1 (4), pp. 193–202. Cited by: §1, §4.
  • A. Bensoussan, J. L. Lions, and G. Papanicolaou (1978) Asymptotic analysis for periodic structures. Studies in mathematics and its applications, North Holland, Amsterdam. Cited by: §1.
  • L. Cheng, J. X. Bai, and A. C. To (2019) 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.
  • Z. Cheng, H. F. Zhou, Q. H. Lu, H. J. Gao, and L. Lu (2018) Extra strengthening and work hardening in gradient nanotwinned metals. Science 362 (6414), pp. 1925. Cited by: §1.
  • J. D. Deng and W. Chen (2017) 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.
  • J. D. Deng, J. Yan, and G. D. Cheng (2013) Multi-objective concurrent topology optimization of thermoelastic structures composed of homogeneous porous material. Struct. Multidiscip. Optim. 47 (4), pp. 583–597. Cited by: §1.
  • H. W. Dong, S. D. Zhao, Y. S. Wang, and C. Z. Zhang (2017) Topology optimization of anisotropic broadband double-negative elastic metamaterials. J. Mech. Phys. Solids 105, pp. 54–80. Cited by: §1.
  • H. Dumontet (1985) 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.
  • P. Fratzl and G. F. Barth (2009) Biomaterial systems for mechanosensing and actuation. Nature 462, pp. 442–448. Cited by: §1.
  • J. Gao, H. Li, L. Gao, and M. Xiao (2018) Topological shape optimization of 3d micro-structured materials using energy-based homogenization method. Adv. Engrg. Softw. 116, pp. 89–102. Cited by: §1.
  • J. P. Groen and O. Sigmund (2018) Homogenization-based topology optimization for high-resolution manufacturable microstructures. Internat. J. Numer. Methods Engrg. 113 (8), pp. 1148–1163. Cited by: §1, §2.2, §5.3.3.
  • J. P. Groen, J. Wu, and O. Sigmund (2019) Homogenization-based 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.
  • X. Guo, W. S. Zhang, J. Zhang, and J. Yuan (2016) 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.
  • X. Guo, W. S. Zhang, and W. L. Zhong (2014) 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.
  • O. Jorgensen, A. E. Giannakopoulos, and S. Suresh (1998) Spherical indentation of composite laminates with controlled gradients in elastic anisotropy. Int. J. Solids Struct. 35 (36), pp. 5097–5113. Cited by: §1.
  • M. S. Kushwaha, P. Halevi, L. Dobrzynski, and B. Djafari-Rouhani (1993) Acoustic band structure of periodic elastic composites. Phys. Rev. Lett. 71 (13), pp. 2022–2025. Cited by: §1.
  • R. Lakes (1993) Materials with structural hierarchy. Nature 361 (6412), pp. 511–515. Cited by: §1.
  • X. Lei, C. Liu, Z. L. Du, W. S. Zhang, and X. Guo (2018) Machine learning-driven real-time topology optimization under moving morphable component-based framework. ASME J. Appl. Mech. 86 (1), pp. 011004–1–011004–9. Cited by: §4.
  • C. Liu, Z. L. Du, Z. Sun, H. J. Gao, and X. Guo (2015) Frequency-preserved acoustic diode model with high forward-power-transmission rate. Phys. Rev. Appl. 3 (6), pp. 064014. Cited by: §1.
  • C. Liu, Z. L. Du, W. S. Zhang, Y. C. Zhu, and X. Guo (2017) Additive manufacturing-oriented 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.
  • L. Liu, J. Yan, and G. D. Cheng (2008) Optimum structure with homogeneous optimum truss-like material. Comput. Struct. 86 (13), pp. 1417–1425. Cited by: §1.
  • S. T. Liu, G. D. Cheng, Y. Gu, and X. G. Zheng (2002) Mapping method for sensitivity analysis of composite material property. Struct. Multidiscip. Optim. 24 (3), pp. 212–217. Cited by: §4.
  • M. A. Meyers, J. McKittrick1, and P. Y. Chen (2013) Structural biological materials: critical mechanics-materials connections. Science 339 (6121), pp. 773–779. Cited by: §1.
  • B. Niu, J. Yan, and G. D. Cheng (2009) Optimum structure with homogeneous optimum cellular material for maximum fundamental frequency. Struct. Multidiscip. Optim. 39 (2), pp. 115–132. Cited by: §1.
  • O. Pantz and K. Trabelsi (2008) A post-treatment of the homogenization method for shape optimization. SIAM J. Control Optim. 47 (3), pp. 1380–1398. Cited by: §1.
  • A. Radman, X. Huang, and Y. M. Xie (2013) Topology optimization of functionally graded cellular materials. J. Mater. Sci. 48 (4), pp. 1503–1510. Cited by: §1.
  • A. Radman, X. Huang, and Y. M. Xie (2014) Maximizing stiffness of functionally graded materials with prescribed variation of thermal conductivity. Comp. Mater. Sci. 82, pp. 457–463. Cited by: §1.
  • H. Rodrigues, J. M. Guedes, and M. P. Bendsoe (2002) Hierarchical optimization of material and structure. Struct. Multidiscip. Optim. 24 (1), pp. 1–10. Cited by: §1.
  • C. Sanchez, H. Arribart, and M. M. G. Guille (2005) Biomimetism and bioinspiration as tools for the design of innovative materials and systems. Nat. Mater. 4, pp. 277–288. Cited by: §1.
  • O. Sigmund and S. Torquato (1996) Composites with extremal thermal expansion coefficients. Appl. Phys. Lett. 69 (21), pp. 3203–3205. Cited by: §1.
  • O. Sigmund (1994) Materials with prescribed constitutive parameters: an inverse homogenization problem. Int. J. Solids Struct. 31 (17), pp. 2313–2329. Cited by: §1.
  • K. Svanberg (1987) The method of moving asymptotes—a new method for structural optimization. Internat. J. Numer. Methods Engrg. 24, pp. 359–373. Cited by: §5.1.
  • P. Vogiatzis, M. Ma, S. K. Chen, and X. F. D. Gu (2018) 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.
  • M. Y. Wang, X. M. Wang, and D. M. Guo (2003) A level set method for structural topology optimization. Comput. Methods Appl. Mech. Engrg. 192 (1), pp. 227–246. Cited by: §4.
  • Y. Q. Wang, F. F. Chen, and M. Y. Wang (2017) Concurrent design with connectable graded microstructures. Comput. Methods Appl. Mech. Engrg. 317, pp. 84–101. Cited by: §1.
  • J. Yan, X. Guo, and G. D. Cheng (2016) Multi-scale concurrent material and structural design under mechanical and thermal loads. Comput. Mech. 57 (3), pp. 437–446. Cited by: §1.
  • W. S. Zhang, D. Li, J. Yuan, J. F. Song, and X. Guo (2017) A new three-dimensional topology optimization method based on moving morphable components (mmcs). Comput. Mech. 59 (4), pp. 647–665. Cited by: §4.
  • W. S. Zhang, W. Y. Yang, J. H. Zhou, D. Li, and X. Guo (2016a) Structural topology optimization through explicit boundary evolution. ASME J. Appl. Mech. 84 (1), pp. 011011–1–011011–10. Cited by: §4.
  • W. S. Zhang, J. Yuan, J. Zhang, and X. Guo (2016b) 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.
  • Y. Zhang, H. Li, M. Xiao, L. Gao, S. Chu, and J. Zhang (2019) 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.
  • Y. Zhang, M. Xiao, H. Li, L. Gao, and S. Chu (2018)

    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.
  • M. Zhou and G. I. N. Rozvany (1991) The coc algorithm, part ii: topological, geometrical and generalized shape optimization. Comput. Methods Appl. Mech. Engrg. 89 (1), pp. 309–336. Cited by: §4.
  • S. W. Zhou and Q. Li (2008) Design of graded two-phase microstructures for tailored elasticity gradients. J. Mater. Sci. 43 (15), pp. 5157. Cited by: §1.
  • Y. C. Zhu, S. S. Li, Z. L. Du, C. Liu, X. Guo, and W. S. Zhang (2019) A novel asymptotic-analysis-based homogenisation approach towards fast design of infill graded microstructures. J. Mech. Phys. Solids 124, pp. 612–633. Cited by: On speeding up an asymptotic-analysis-based 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.
  • Y. C. Zhu, Y. H. Wei, and X. Guo (2017) Gurtin-murdoch surface elasticity theory revisit: an orbital-free density functional theory perspective. J. Mech. Phys. Solids 109, pp. 178–197. Cited by: §6.