Topology optimization of multi-scale structures: a review

Abstract

Multi-scale structures, as found in nature (e.g., bone and bamboo), hold the promise of achieving superior performance while being intrinsically lightweight, robust, and multi-functional. Recent years have seen a rapid development in topology optimization approaches for designing multi-scale structures, but the field actually dates back to the seminal paper by Bendsøe and Kikuchi from 1988 (Computer Methods in Applied Mechanics and Engineering 71(2): pp. 197–224). In this review, we intend to categorize existing approaches, explain the principles of each category, analyze their strengths and applicabilities, and discuss open research questions. The review and associated analyses will hopefully form a basis for future research and development in this exciting field.

1 Introduction

Topology optimization is a computational design method for automatically generating a structural layout with maximized performance under relevant design specifications. In other words, the structural design problem can be formulated as optimizing the distribution of material in a discretized design domain (Bendsøe 1989). The optimized layout is not restricted to its initial topology, opening for superior structural performance over manual designs based on engineers’ intuition and experience. This capability is especially attractive and has been successfully applied in high-tech industries, such as aerospace (Zhu et al. 2016), automotive (Yang and Chahande 1995), architecture (Beghini et al. 2014), and healthcare (Wang et al. 2016).

Early works that lead to the establishment of this field, and in particular the seminal work by Bendsøe and Kikuchi (1988), made use of a material model corresponding to infinitely small square cells with rectangular holes and evaluated equivalent mechanical properties of these cells by numerical homogenization. Theoretically, optimal structures comprise spatially varying geometric patterns spanning multiple length scales. Due to manufacturing difficulties of these multi-scale structures, focus in the late 1990s shifted from homogenization-based approaches to “mono-scale” approaches, optimizing the distribution of a homogeneous isotropic material (Bendsøe 1989; Zhou and Rozvany 1991; Mlejnek 1992; Bendsøe and Sigmund 1999). Popular methods such as those based on density (Sigmund 2001), level sets (Wang et al. 2003; Allaire et al. 2004), and evolutionary procedures (Xie and Steven 1993) all belong to this latter category.

Over the past few years, along with advances in additive manufacturing (AM, also known as 3D printing), there has been a resurgent interest in optimal design of multi-scale structures. In additive manufacturing, parts are produced layer upon layer by, e.g., extruding small flattened strings of molten material or melting and fusing powder material or wire. AM provides an effective means to fabricate complex mono-scale structures as well as delicate multi-scale structures. In fact, topology optimization of (mono-scale) structures and the use of lattice infill are recognized as two dominating strategies for designing next-generation lightweight structures (Brackett et al. 2011; Thompson et al. 2016; Plocher and Panesar 2019). The combination, i.e., topology optimization of multi-scale structures, thus holds the promise of superior performance in a general sense. These expectations are sometimes coupled with observation of multi-scale structures found in nature (e.g., bone and bamboo (Lakes 1993; Fratzl and Weinkamer 2007)), while the specific function and the driving mechanisms underlying those natural multi-scale structures are not always known precisely.

Figure 1 illustrates a multi-scale structure. Each point in the macroscale structure effectively represents a periodic repetition of a local microstructure. Here, as common in the topology optimization literature, microstructure shall be interpreted relatively—it does not specifically refer to a physical size under, e.g., 1 mm but rather to a scale much smaller than the macroscale. It is thus assumed that continuum mechanics is applicable to both macro- and microscales. Some examples of 3D multi-scale structures designed by topology optimization are shown in Fig. 2.

Fig. 1

Illustration of a multi-scale structure

Fig. 2

Some examples of 3D multi-scale structures designed by topology optimization. From left to right: a bone-inspired infill structures, b variable-density lattice optimization, and (c, d, e) three de-homogenization schemes in 3D. The first three are produced by 3D printing, while the last two are renderings. a and c Reprinted with permission from IEEE, from Wu et al. (2018) and Wu et al. (2021), respectively. b, d and e Reprinted with permission from Elsevier, from Zhang et al. (2017), Geoffroy-Donders et al. (2020)and Groen et al. (2020), respectively

In view of the rapid development of topology optimization of multi-scale structures, it seems timely to critically review the variety of methods that have been proposed. The term multi-scale has been used extensively in the literature (and also in this review) to describe structures and modelling techniques, as well as design approaches. Many design approaches make use of multi-scale modelling, i.e., they assume separation of length scales. However, mono-scale modelling can also be used to design multi-scale structures. In this review, we intend to categorize existing approaches, explain the principles of each category, analyze their strengths and applicabilities, and discuss open research questions. This review is supplemented with our implementation of representative approaches, based on which we quantitatively evaluate structural performance. This review and analysis hopefully will form a basis for future research and development in this exciting field.

2 Background

In this section, we categorize three fundamental formulations that provide the basis for most multi-scale approaches.

2.1 Homogenization-based structural optimization

A starting point in the history of structural optimization is the classical work by Michell (1904) where an optimal truss design is represented by a continuum description. Arguably, this is the first work on multi-scale optimization since the truss-like continuum description is actually a limit case that can only be realized by using an infinite number of discrete truss members; see Prager and Rozvany (1977) for more details. Similarly, Cheng and Olhoff (1981) found that for the case of optimizing the stiffness of a plate by varying the thickness, more stiffener-like reinforcement appeared when the design space was refined. They made the important conclusion that in the limit of an infinitely fine mesh an infinite number of stiffener-like reinforcing members would occur (Cheng 1981). Only by restricting the variation in shape can existence of a solution be proven (Niordson 1983). For compliance minimization in 2 or 3 dimensions, Kohn and Strang (1986) proved that the optimal material distribution can only be found by relaxing the design space. This means the use of a continuous composite material description using a characteristic unit-cell length $\epsilon \rightarrow 0$, allowing for much more detail than a single-scale discretized point-wise material or void description. Here, the theory of homogenization (Bensoussan et al. 1978) enters the field of structural optimization, since this is the perfect tool to bridge the scale between the microscopic periodic composite material and its effective (homogenized) properties on the macroscale. Using this theory, several research groups independently and more or less simultaneously realized that there exists a class of sequential laminates, the so-called rank-N laminates that can achieve the theoretical upper bounds for maximum strain energy density (Lurie and Cherkaev 1984; Norris 1985; Francfort and Murat 1986; Milton 1986). These composite materials consist of several scales in themselves, which is a necessary condition to achieve ultimate stiffness (Allaire and Aubry 1999). With this knowledge, Bendsøe and Kikuchi (1988) came up with the landmark paper that provided the computational framework to do homogenization-based topology optimization of continuum structures. Contrary to using rank-2 laminates that are optimal for a single load case in 2D, they came up with a single-scale interpretation, the square microstructure with a rectangular hole, shown in Fig. 3. By changing the widths μ₁ and μ₂ as well as the orientation 𝜃 throughout the design domain, near-optimal structures could be achieved. The use of rank-2 laminates as a microstructure was included in later works (Bendsøe 1989; Allaire and Kohn 1993).

Fig. 3

Layout of the unit cell with a rectangular hole, in local (y₁, y₂), and global (x₁, x₂) coordinate systems

2.2 Density-based topology optimization

Shortly after the homogenization-based topology optimization approach was introduced, an alternative known as the SIMP (Solid Isotropic Material with Penalization) or power-law approach was suggested (Bendsøe 1989; Zhou and Rozvany 1991; Mlejnek 1992). Here, the material distribution in a design domain is represented by a scalar field, with a relative density per element in the discretized domain. Here, ρ = 1 means solid and ρ = 0 means empty. This integer-programming problem is relaxed to allow for intermediate densities during optimization. Regardless of the value of the relative density, the material within each element is assumed to be isotropic and homogeneous. Material properties (e.g., Young’s modulus) are related to relative densities by a power-law interpolation. This seemingly artificial relationship was later proved physically permissible for a properly chosen penalty power, e.g., with a power p ≥ 3 for Poisson’s ratio ν = 1/3 (Bendsøe and Sigmund 1999). Based on this, a heterogeneous material distribution within the element domain, i.e., a microstructure of solid base material and void, can be found to match the expected material properties of an intermediate density. Hence, even the simple “mono-scale” density approaches may be interpreted as multi-scale approaches albeit with suboptimal microstructures. With a justified penalization power, the SIMP approach converges to near 0-1 solutions, representing mono-scale structures. In a number of multi-scale approaches to be discussed later, a smaller p is chosen for obtaining solutions with a large portion of intermediate densities, which provide a basis for filling porous microstructures. However, this may not be mechanically valid, c.f. if p < 3 for ν = 1/3, this power-law scheme violates the Hashin-Strikhman bounds (Hashin and Shtrikman 1963). The SIMP approach is easy to implement; see educational codes (Sigmund 2001; Andreassen et al. 2011; Ferrari and Sigmund 2020). Here, the three-field projection approach is recommended (Guest et al. 2004; Xu et al. 2010; Wang et al. 2011).

Mono-scale structures can also be optimized by methods based on level sets (Wang et al. 2003; Allaire et al. 2004), evolutionary procedures (Xie and Steven 1993), and geometric morphing and projection (Norato et al. 2004; Guo et al. 2014). For an overview of these methods, we refer to the reviews by Sigmund and Maute (2013) and Deaton and Grandhi (2014). The use of non-gradient approaches for solving topology optimization problems with many variables is not recommended (Sigmund 2011).

2.3 Inverse homogenization

The theory of numerical homogenization, which approximates effective properties of a periodic unit cell, can also be exploited to design a microstructure with desired properties. This method was introduced by Sigmund (1994) and is generally referred to as inverse homogenization. With compliance minimization in mind, it is common practice to formulate an optimization problem to maximize the microstructure stiffness with respect to the applied stresses or strains. For the difference in stiffness between such microstructures and the optimal energy bounds, see, e.g., Guedes et al. (2003) and Träff et al. (2019). Furthermore, microstructures with negative Poisson’s ratio (Sigmund 1994; Andreassen et al. 2014; Clausen et al. 2015b), materials with maximum shear and bulk moduli (Sigmund 2000), or materials with increased buckling strength (Thomsen et al. 2018; Wang and Sigmund 2020) can be designed using inverse homogenization.

Besides problems in elasticity, inverse homogenization can be applied to many other types of physics. Sigmund and Torquato (1996) considered thermal expansion problems, while Torquato et al. (2002) included both thermal and electrical conductivity. Sigmund (1999) and Challis et al. (2012) considered combined stiffness and fluid permeability. Inverse homogenization for the design of phononic and photonic materials has been considered by Sigmund and Jensen (2003) and Jensen and Sigmund (2011), respectively. Lately, inverse homogenization has been applied to the design of photonic isolaters by Christiansen et al. (2019). For more details on material design using inverse homogenization, we refer interested readers to the review papers by Cadman et al. (2013) and Osanov and Guest (2016) and the Ph.D. thesis by Andreassen (2015).

3 Full-scale approaches

As mentioned earlier, multi-scale structures can be designed using either multi-scale or mono-scale modelling. In this review, approaches based on multi-scale modelling (i.e., with the assumption of separation of length scales) are referred to as multi-scale approaches. Approaches that do not make this assumption, i.e., “mono-scale” approaches, optimize distribution of a homogeneous material. When the design domain is discretized by a finite-sized mesh, such “mono-scale” approaches typically result in mono-scale structures. However, as the discretization of the design domain goes high, it can directly be used to achieve multi-scale structures, since theoretical stiffness optimal structures span multiple scales. By employing careful continuation techniques and sufficiently fine meshes, and in the absence of regularization for mesh independence such as control of minimum length, perimeter, or slope, multi-scale structures should appear naturally. However, the appearance of fine scale structures may also be stimulated by controlling the layout locally. In these approaches, analysis and optimization of structures are performed in the full resolution of the domain, and we thus refer to them as full-scale approaches.

For local control approaches, two strategies can be distinguished: pattern repetition and local volume constraints. In the former, the design domain is partitioned into a number of subdomains that are further refined. The layout in each subdomain is enforced to be the same as that in the others, leading to periodic patterns in the full domain. In a variation, the subdomains are grouped, and an identical layout is enforced in subdomains per group, resulting in, for instance, periodic patterns along one axis with gradation along another. The enforcement of an identical layout in subdomains can be achieved by creating a template that links all subdomains. The sensitivity of each point in the template is the sum of sensitivities of corresponding points in all linked subdomains. This strategy reduces the solution space from the full design domain to one subdomain, or a few for a gradation or partitioning, but the analysis is still performed on the full mesh.

Pattern repetition is pursued in full-scale approaches as well as multi-scale approaches. We discuss the former here, and will return to the latter in subsequent sections. Zhang and Sun (2006) presented a first and comprehensive study on creating periodic patterns using topology optimization. It was shown that topology design results were greatly influenced by both the number and aspect ratio of the subdomains. Similar results were later reproduced using evolutionary procedures (Huang and Xie 2008) and level sets (Liu et al. 2018b). Almeida et al. (2010) presented a density-based approach to create pattern repetition as well as structural symmetry. It is worth noting that the subdomains do not have to be of the same size. Repetitive pattern with a variation in size was demonstrated by Stromberg et al. (2011) and Wu et al. (2016a) for architectural design, by establishing a mapping between points in size-varying subdomains. Lately, the principle of pattern repetition has been used to design mechanisms (Wu et al. 2020). Designing periodic patterns using full-scale analysis involves intensive computation. Zhang and Sun (2006) proposed to reduce computation time through condensation—unfortunately, the details were missing in the original paper. Alexandersen and Lazarov (2015) showed an enormous reduction in computational cost by using a spectral coarse basis pre-conditioner that exploits the repeated patterns, i.e., for each repeated pattern the same spectral basis can be used.

The second strategy to stimulate fine scale structures is to apply local volume constraints (Wu et al. 2018). The idea is to impose an upper bound on the fraction of solid elements in the neighborhood of every point in the full design domain. This involves two important parameters: the radius of the neighborhood which affects the pore size, and the prescribed upper bound on the local volume fraction which controls the porosity. Since this creates a large number of constraints (equal to the number of elements), aggregation schemes such as the p-norm are often applied to aggregate the per-element constraints into a global one, facilitating the optimization process. It is found that local volume constraints yield bone-like porous structures that are aligned with principal stress directions (see Fig. 2a). These structures are robust with respect to load variations, local material failure, and buckling, at the price of some decrease in stiffness. A typical way to ensure robustness is to model uncertainties such as material failure and load variations in topology optimization, resulting in distributed structures (e.g., Jansen et al. (2014)). This demands multiple finite element analyses per optimization iteration. In this regard, applying local volume constraints is an efficient alternative for robust design. The concept of local volume constraints is similar to maximum length scale control, which was introduced by Guest (2009) and later studied in Lazarov and Wang (2017) and Carstensen and Guest (2018).

A number of further developments of the local volume constraints have been proposed. Wu et al. (2017) combined it with a coating approach (Clausen et al. 2015a) to simultaneously evolve the (macroscale) structure and the microstructures therein, referred to as shell-infill composites. Gradation in the porosity and pore size was demonstrated by Schmidt et al. (2019). Dou (2020) reformulated the constraints by a projection method, while Cang et al. (2019) made use of machine learning algorithms to predict optimized structures with local volume constraints. This was recently extended to incorporate multiple materials (Li et al. 2020). Furthermore, the local volume constraints have been applied in the context of other physics problems, e.g., heat conduction (Yan et al. 2018) and structural mechanics coupled with heat conduction (Das and Sutradhar 2020).

Figure 4 compares compliance-minimized structures optimized using a conventional mono-scale formulation under a total volume constraint (left), with pattern repetition (middle) and with local volume constraints (right). From this comparison it can be concluded that:

• Both pattern repetition and local volume constraints restrict the solution space, and the structure is expected to be less optimal than obtained from a formulation on the same resolution without these constraints. While more tests are not included here, it can be found that the restriction becomes more pronounced if the size of the repetitive pattern or the filter radius for defining local volume is reduced. Oppositely, when the pattern size or filter radius is larger than the size of the domain, the constraints have no influence on the solution; the local volume constraints degenerate to the conventional global volume constraint.

• Porous structures from local volume constraints exhibit continuous variations in orientation, while periodic patterns have a constant orientation. In this sense, local volume constraints are less restrictive than pattern repetition in constraining the optimization problem. Note that pattern repetition also implicitly controls the volume fraction in subdomains.

Fig. 4

Full-scale optimized structures using the same amount of material. The design space is a rectangular domain discretized by 400 × 200 square elements. The boundary condition is illustrated on the reference design (left). Middle: Optimized design with enforced pattern repetition. The size of the pattern is 40 × 40, indicated by a blue transparent square. Right: Optimized design with local volume constraints. The radius of the filter for calculating local volumes is 20, indicated by a blue transparent disk. From left to right, the compliance values are 100.56, 282.82, and 152.28, respectively

4 Classification of multi-scale approaches

Many works on multi-scale optimization repeat the theory to do numerical homogenization using finite elements, which was described in detail by Guedes and Kikuchi (1990). Based on this theory, Andreassen and Andreasen (2014) provided an open-source educational Matlab code to do numerical homogenization for plane problems, which has later been extended to 3D by Dong et al. (2019). However, more important is how to accurately define the multi-scale optimization problem. Therefore, we will now review some theoretical aspects about the problem formulation for compliance minimization problems for n_load load cases. Using spatially varying microstructures, we can write the corresponding optimization problem as (e.g., Rodrigues et al. 2002),

$$ \underset{{\int}_{{\varOmega}} \rho \text{d}{\varOmega} \leq V_{\max}}{\underset{\rho}{\max}} \underset{\boldsymbol{E}^{H}(\rho) \in \boldsymbol{E}_{ad}}{\max} \underset{\boldsymbol{u}^{k} \in \boldsymbol{U}}{\min} \left\{ \sum\limits_{k=1}^{n_{load}} w^{k} {\varPi}_{tot}^{k} \right\}. $$

(1)

Here, ρ is the macroscale variable describing the porosity of the varying microstructures, which is subject to an upper bound on the available material V_max. E^H describes the homogenized elasticity tensor of the microstructures, which has to be in the physically admissible set of elasticity tensors E_ad. The inner problem describes the minimization of potential energy ${\varPi }_{tot}^{k}$, i.e., the equilibrium constraint that has to hold for each load case k, where we solve for displacement field u^k in the space of kinematically admissible displacement fields U. Finally, w^k is a weighting factor to scale the energy for each load case. The compliance $\mathcal {J}$ can be calculated as,

$$ \mathcal{J} = -2\left\{ \sum\limits_{k=1}^{n_{load}} w^{k} {\varPi}_{tot}^{k} \right\}, $$

(2)

hence, maximizing the total potential energy is equivalent to minimizing the compliance.

Optimization problem (1) can be solved in different ways, e.g., as one system of equations using the simultaneous analysis and design (SAND) approach (Jog and Haber 1996), or using a sequence of separable convex approximations using a so-called nested analysis and design (NAND) approach (Fleury 1993). The latter is the de facto standard for topology optimization problems. However, solving problem (1) would, depending on the way the material of the microstructures is modeled, result in a large KKT matrix for the maximization with respect to E^H. To circumvent this, we can interchange the maximization with respect to E^H and the minimization with respect to u,

$$ \underset{{\int}_{{\varOmega}} \rho \text{d}{\varOmega} \leq V_{max}}{\underset{\rho}{\max}} \underset{\boldsymbol{u}^{k} \in \boldsymbol{U}}{\min}\! \left\{\! \underset{\boldsymbol{E}^{H}(\rho) \in \boldsymbol{E}_{ad}}{\max}\! \left\{\! \mathcal{E}(\boldsymbol{u}^{k},\rho)\! \right\}\! -\! \sum\limits_{k=1}^{n_{load}}\! w^{k} {\varPi}_{ext}^{k}\! \right\}.\! $$

(3)

Here, π_ext is the potential energy of the applied loads for the k-th load case, and the point-wise optimal strain energy $\mathcal {E}$ can be written as,

$$ \mathcal{E} = \frac{1}{2} \sum\limits_{k=1}^{n_{load}} w^{k} \Big({\int}_{{\varOmega}} \boldsymbol{\varepsilon}(\boldsymbol{u}^{k}):\boldsymbol{E}^{H}(\boldsymbol{\rho}): \boldsymbol{\varepsilon}(\boldsymbol{u}^{k})\text{d}{\varOmega} \Big). $$

(4)

Hence, we now have an inner point-wise microstructure optimization problem. As proven by Lipton (1994b), problem (3) is equivalent to problem (1) if a microstructure parameterization is used such that the strain energy is concave with respect to E^H. This concavity is a property of optimal rank-N microstructures. Hence, if suboptimal microstructures are used, the optimality of the solution to problem (3) cannot be guaranteed. Nevertheless, problem (3) results in an optimization problem that is independent of the material description and thus forms the basis for the vast majority of multi-scale topology optimization approaches. The most common way to solve the optimization problem (3) is to use a hierarchical solution procedure (see, e.g., Bendsøe and Sigmund (2004) for a detailed description), which can be summarized in the flowchart below.

• 1. Initialize data: Geometry, mesh, boundary conditions, density ρ, microstructure description E^H, etc.

• 2. Analysis: Solve for macroscale displacements u^k.

• 3. Microstructure optimization: Solve point-wise for E^H by maximizing $\mathcal {E}$ for given strains and ρ.

• 4. Density update: Update the outer maximization problem, subject to volume constraint.

• 5. Check convergence: Check if the optimization problem is converged. If not, repeat steps 2–5.

It is not necessary to solve the problem exactly like this. For example, one can combine the global density update and all local problems in one large optimization step (see, e.g., Bendsøe and Kikuchi (1988)). However, for problems with a large number of design variables, the separation in local and global problems allows for reducing the computational time through the use of parallel computations (Coelho et al. 2008).

An alternative to the hierarchical formulation (problem (1)) is to assign two sets of design variables, i.e., on both macro- and microscales, and concurrently optimize these two sets. This concurrent formulation, known as PAMP (Porous Anisotropic Material with Penalization), was proposed by Liu et al. (2008) to design microstructures and their macroscopic layout.

Problem (3) and the abovementioned flowchart will be used to categorize multi-scale approaches later on. Although the equation concerns compliance minimization in the elasticity setting, similar looking equations and flowcharts can be made for different types of physics including additional constraints. Each of these multi-scale problems can be divided into a sequence of three related optimization problems, i.e., a state problem; a global parameter distribution problem; and one or many local microstructure optimization problems. The key point in classifying the different multi-scale topology optimization approaches is to determine what kinds of restrictions are applied to the solution space of each of the sub-problems. For now, we assume that at least the state problem is solved up to full accuracy in each of the multi-scale problems, although this may not apply to deep learning or iterative solution approaches. Hence, we classify the different approaches by (1) the restrictions that are applied to the density distribution and (2) the restrictions that are applied to the admissible set of properties E_ad that can be achieved by the parameterized microstructure.

It should be mentioned that even in the elasticity setting it is not for all cases exactly known which constitutive tensors are realizable (Milton et al. 2017). For example, microstructures that reach minimum shear modulus at the Walpole point have never been realized (Sigmund 2000). However, as discussed above, rank-N laminates cover the full space of optimal designs for compliance minimization problems. In other words, if one uses rank-N laminates for the microstructure description one can solve problem (3) to the true optimum. It is known that for a single load case in a 2D setting, rank-2 laminates are sufficient to describe the optimal solution, while for multiple load cases a rank-3 laminate is required (Avellaneda 1987). Hence, if one would only use a rank-2 parameterization for a multiple load case problem, the optimality could be severely reduced. Since the vast majority of works on multi-scale topology optimization use (strongly) restricted sets E_ad (for computational, manufacturing or other reasons), we categorize approaches based on the following restrictions, starting with the least restricted category:

• I Optimal set of elasticity tensors: E_ad is represented by a geometry parameterization that allows the local problem to be solved to optimality. This is for elasticity and conduction problems the set of rank-N laminates (Allaire 2002). However, note that, for many other types of problems, the optimal set of constitutive tensors is not exactly known yet.

• II Unrestricted unit-cell design: E_ad contains the set of unit cells that can be obtained using inverse homogenization, without restrictions on the material distribution, shape, connectivity, or orientation of the unit cell. This means that if a fine enough discretization is used, the microstructures should converge to what is theoretically possible. However, in practice, limited mesh resolution will cause suboptimal microstructures.

• III Restricted unit-cell design: E_ad contains the set of unit cells that can be obtained using inverse homogenization with restrictions on material distribution, unit-cell shape, connectivity, or orientation. For example, this can be a square unit cell or a design with pre-defined solid elements, both resulting in severely restricted design freedom.

• IV Parameterized unit cell with multiple parameters: E_ad contains a set of pre-computed parameterized unit cells such that the microstructure is parameterized using more than 1 parameter. For example, this can be the rectangular hole microstructure by Bendsøe and Kikuchi (1988). Due to the rotational freedom, the rectangular hole cell actually performs almost as well as rank-2 microstructures (Bendsøe and Sigmund 1999), if properly de-homogenized.

• V Parameterized unit cell based on density: E_ad contains a single constitutive tensor E^H for a given microstructure density ρ. This is the most restrictive case since it does not involve a local optimization problem. Isotropic microstructures satisfying the Hashin-Shtrikman bounds (Hashin and Shtrikman 1963) fall in this category since the isotropic elasticity tensor follows only from the density. Likewise, the SIMP approach with a penalty parameter p that satisfies the Hashin-Shtrikman bounds (Bendsøe and Sigmund 1999) falls into this category.

The difference between category III where the microstructure topology is optimized during the optimization process and category IV is that for category IV a database of pre-computed microstructural properties is generated. Hence, category IV will require a heavy pre-computation step, but afterwards the optimization problem can be solved more efficiently.

There are many works that consider a set of elasticity tensors E_ad outside of the bounds on what is theoretically possible. An example of this is the free material optimization (FMO) approach (Bendsøe et al. 1994), which has also been used in two-scale structural optimization (Schury et al. 2012; Hu et al. 2020). In this method, the entries of E^H are the design variables, and using either a constraint on a norm or on eigenvalues of E^H one can try to relate the tensor to a density distribution. Nevertheless, there is no direct relation between the bounds used on the tensor and the Hashin-Shtrikman bounds; and therefore, we will not discuss the FMO method further in this review. Another good example of violation of the bounds is the variable thickness sheet problem (Rossow and Taylor 1973), i.e., the SIMP method with p = 1. In 2 dimensions, we can justify these designs by interpreting density as a thickness in the out-of-plane direction. Hence, it should therefore be seen as a sizing problem. However, the model becomes invalid if used to pre-compute a global density distribution for a subsequent porous microstructure realization. One cannot make an in-plane realization of these variable thickness designs and therefore works that use SIMP with p = 1 for any porous or 3D design problem should be avoided.

Finally, we can categorize the outer problem into three different categories.

• A Unrestricted density: There are no restrictions on the density, i.e., ρ ∈ [0,1].

• B Restricted density: Only a few values of ρ are allowed. For example, this applies to the SIMP method (possibly combined with a projection method (Guest et al. 2004) to reduce elements with intermediate densities), and also the PAMP approach (Liu et al. 2008) falls into this category. This also applies to interface bounded approaches (Clausen et al. 2015a; Groen et al. 2019; Luo et al. 2019) with a fixed infill density and a solid outer shell.

• C Fixed density: The density field is fixed, i.e., there is no outer optimization problem. This is, for example, a uniform density field or a density distribution based on some prior optimization problem. Note that a fixed density does not necessarily restrict the unit cell to be periodic, though this is the case in Cherkaev et al. (1998) and Fujii et al. (2001). A fixed density permits unit cell adaptation, e.g., rotation (Wu et al. 2021).

In total, we can thus identify 5 × 3 different categories of multi-scale topology optimization problems, which are summarized in Table 1. Here, we also list a fundamental paper for each category (if available).

Table 1 Categories of multi-scale topology optimization problems and representative papers

5 Multi-scale approaches

This section is organized in accordance with the above classification based on restrictions on unit cells. We start with the least restricted category—the optimal set of constitutive properties (I), and then proceed directly to the most restricted one where unit cells are parameterized by a single parameter, i.e., density (V). This is followed by categories with increasing flexibility in the unit cell parameterization (II, III, IV). This order is found to align with the chronological order of the first article in each category (see Table 1, column A).

5.1 Optimal set of constitutive properties

In the context of elasticity, the minimum number of sequential layers that is required to parameterize the optimal solution using rank-N microstructures depends on both the number of load cases and the dimension of the problem (Avellaneda 1987; Francfort and Murat 1986):

• Rank-2 laminates with orthogonal layers, for plane problems subject to a single load case

• Rank-3 laminates, for plane problems subject to multiple load cases

• Rank-3 laminates with orthogonal layers, for problems in 3D subject to a single load case

• Rank-6 laminates for problems in 3D subject to multiple load cases

Hence, at most, 7 different length scales are needed to describe an optimal material parameterization in 3D, i.e., six length scales for the microstructure and one for the macroscopic parameterization. The elasticity tensors of these multi-scale laminates can be derived analytically, allowing the local microstructure optimization problems to be solved in a very efficient manner. For example, for a rank-2 laminate shown in Fig. 5 (right) the first layering is constructed on the infinitesimal length scale x/𝜖² shown in Fig. 5 (left). This layering can be described by the relative layer width μ₁ that describes the ratio of the isotropic stiff material (+) compared to the weak void material (-) and the interface normal n¹ and tangent t¹. Subsequently, the rank-2 laminate is constructed using a second layering on length scale x/𝜖 using relative layer width μ₂ that describes the ratio of the stiff material to the first layering. By ensuring orthogonal interface vectors, the elasticity tensor is thus defined by only three parameters, μ₁, μ₂, and 𝜃₂, besides the material properties of the stiff (+) and weak (-) phase respectively.

Fig. 5

Visualization of how a rank-2 microstructure is constructed from a stiff isotropic material (+) and a weak/compliant isotropic material (-). Please note the different length scales, on the left x/𝜖² and on the right x/𝜖, with the characteristic unit-cell length $\epsilon \rightarrow 0$

The first work in which topology optimization using rank-2 laminates was considered is the work by Bendsøe (1989). Subsequently, Allaire and Kohn (1993) proved that the optimal design for plane problem subject to a single load case can be described purely by the stress distribution in the macroscopic domain. Multiple load case problems, where the design is parameterized using optimal rank-3 laminates, have been considered for both plane problems (Allaire et al. 1996; Cherkaev et al. 1998) and plate/shell problems (Díaz et al. 1995; Hammer et al. 1997; Krog and Olhoff 1999). By choosing different plate layouts, e.g., using a fixed core thickness, Soto and Díaz (1993) and Krog and Olhoff (1999) showed the effect of using several density restrictions on the performance. The effect of restricting the amount of unique microstructures in the domain has been investigated by Cherkaev et al. (1998), who demonstrated that the performance can be at least twice as bad if only 1 or a few unique micostructures are used in the design domain.

The method has been applied to 3D topology optimization problems, where orthogonal rank-3 laminates are used for single load case problems (Cherkaev and Palais 1996; Allaire et al. 1997; Díaz and Lipton 1997; Olhoff et al. 1998; Czarnecki and Lewiński 2006). Multiple load cases are considered by Díaz and Lipton (2000). A nice feature of rank-N laminates is that they provide a lower bound on the theoretical compliance that can be reached. It is thus recommended to include this comparison when presenting a multi-scale approach. To help doing so, we supply a code for the optimization of rank-2 microstructures, which is discussed in Appendix.

A downside of structures optimized using rank-N laminates is that they consist of several length scales. This poses two challenges. Firstly, although technology is rapidly developing and lab work has demonstrated multi-scale manufacturing capability (e.g., Zheng et al. (2016)), a majority of manufacturing techniques do not yet support precise production of structures spanning multiple length scales. Hence, single scale interpretation of multi-scale rank-N laminates is required. Träff et al. (2019) have shown that the multi-scale rank-N microstructures can be approximated on a single scale with only a small loss in performance. Secondly, it is necessary to compile globally consistent structures from locally defined rank-N laminates (or their simplified versions), i.e., to de-homogenize the results. This problem was first addressed by Pantz and Trabelsi (2008), who used an implicit geometry description to interpret the optimized multi-scale designs on a single length scale with little loss in performance. More details on de-homogenization will be presented later in Section 5.3.1.

Another difficulty with the interpretation of rank-N microstructures is that, for multiple load cases, there exist an infinite number of different microstructure parameterizations that can reach the same elasticity tensor. As an example, consider the planar unit-cell optimization problem using four independent load cases shown in Fig. 6 (top). If we optimize for minimum complementary energy and find the corresponding rank-3 laminates using the method of Lipton (1994a), we obtain the four distinct microstructure realizations shown in Fig. 6 (bottom). These microstructures and an infinite number of combinations thereof can achieve the same minimum complementary energy, which is both an advantage and a difficulty when finding practical interpretations of rank-N designs.

Fig. 6

Top: Illustration of the four stress cases and their respective weights. Bottom: Visualizations of optimal rank-3 laminates with indicated hierarchy, optimized for the 4 stress cases using ρ = 0.5 (Groen 2019)

Finally, it should be mentioned again that for many problems beyond elasticity and heat conduction, a full parameterization of the optimal set of constitutive properties is yet to be found. Finding these bounds on the properties and geometries that reach these bounds is a research field on its own. The interested reader is referred to the books on the theory of composites by Milton (2002, 2016).

5.2 Parameterized unit cell based on density

Predefined cellular structures such as those shown in Fig. 7 have been used to efficiently design multi-scale components. Common in many of these multi-scale approaches is the use of homogenization to evaluate the effective material properties of predefined unit cells. Since the cells are predefined, homogenization can be performed off-line prior to optimization. This allows to generate structures with fine geometric details at a run-time computational complexity comparable to mono-scale approaches (Section 2.2). This is attractive, as evidenced by the large number of publications and industrial examples falling into this category. Predefined unit cells reduce the solution space of the microstructures to a dimension of 1 (or a few, see Section 5.3). This distinguishes this approach from hierarchical or concurrent approaches (Sections 5.4 and 5.5), where microstructures are represented by multi-variable density fields. This distinction is meaningful for implementation. Since the solution space is reduced, it is convenient to construct a differentiable function that maps the reduced parameters to homogenized properties, and use this function in macroscale optimization.

Fig. 7

Three types of 3D unit cells, from left to right, truss- or beam-like lattice structure, triply periodic minimal surface (TPMS), and topology optimized microstructure (Reprinted from Andreassen et al. 2014 with permission from Elsevier)

A typical application of these approaches is to design 3D printed components with uniform infill patterns. From the optimization perspective, this requires no or little adaption of mono-scale optimization approaches. The repeating unit cell can be interpreted as a material, the distribution of which is then optimized by mono-scale approaches. Note that the unit cells interpreted as a material can be isotropic or anisotropic. An example of the former is a 2D beam-like lattice structure constructed from a regular triangular tiling (Clausen et al. 2016). In contrast, the 2D lattice structure from a regular square tiling is anisotropic. Most mono-scale optimization approaches in their standard form deal with isotropic materials. To handle anisotropic materials, the stiffness matrices in finite element analysis must be adapted to reflect the homogenized elasticity tensor of the anisotropic material. Also note that if not optimally oriented along local principal stress directions, anisotropic microstructures may result in highly suboptimal results.

There has been a growing interest in designing components with graded cellular structures. Here, the input to optimization extends from a single unit cell to a set of cells with gradation in porosity and effective material properties (e.g., elasticity tensor). This set can thus be referred to as functionally graded cells. They are typically parameterized by the fraction of solid material within each cell (ρ), naturally serving as the design variable in optimization. The homogenized elasticity tensor (E^H) is a function of ρ, E^H(ρ), which can be constructed by interpolating elasticity tensors of sample cells with equally spaced material fractions. E^H(ρ) thus depends on the specific functionally graded cells, and in general it deviates from the power-law relation. This function for X-shaped lattice structures with varying thicknesses is visualized in Fig. 8a, where the power-law curves with p = 1 (i.e., unpenalized), 2 and 3 are also plotted. For accurately representing the corresponding mechanical properties, the cell-specific E^H(ρ) shall be used in lieu of the generic power-law in density-based topology optimization. The cell-specific E^H(ρ) is commonly constructed using numerical homogenization. Wu et al. (2019) proposed an alternative surrogate model for mapping the effective properties of density parameterized unit cells. It used static condensation to reduce the degrees of freedom of unit cells, followed by proper orthogonal decomposition and diffuse approximation for mapping the density to unit cell stiffness matrix.

Fig. 8

Comparison of different interpolations for topology optimization of a cantilever beam using a volume fraction of 25%. a Interpolation functions for square-shaped (E^s) and X-shaped (E^x) lattice cells with varying member thicknesses. The samples of square- and X-shaped cells are marked by $\square$ and ×, respectively. b, c Optimized density distributions using square- and X-shaped cells with cell-specific interpolation functions, respectively. d, e Optimized density distributions using a SIMP interpolation, with p = 1 and 3, respectively. In d and e, c indicates the compliance estimated by the generic (unphysical) power-law, while c^s and c^x refer to the compliance re-evaluated with cell-specific elasticity tensors. SIMP with p = 1 d significantly overestimates the stiffness of lattice cells. f A reference design by using rank-2 material (obtained with code in Appendix)

Cell-specific material interpolation models have been used in the optimization for components consisting of gyroid-based cellular structures (Li et al. 2018a) and beam-like lattice structures (Wang et al. 2018; Watts et al. 2019), as well as topology optimized microstructures (Garner et al. 2019). Graded cellular structures can also be obtained by post-processing an unmodified SIMP-based mono-scale topology optimization approach, e.g., by replacing the intermediate densities by unit cells of the corresponding material fractions (Brackett et al. 2011; Panesar et al. 2018). Figure 8 compares interpolation schemes for two different unit cells (c.f., (b) and (c)) with varying member thicknesses. The comparison is performed on the 2D cantilever problem, thoroughly studied by Sigmund et al. (2016). The domain is fixed at its left edge and loaded with a unit vertical load distributed over the central 10% of the right edge. Young’s moduli for the solid and void material are E₀ = 1 and $E_{\min \limits }=0.001$, respectively¹. Poisson’s ratio is ν₀ = 0.3 for both materials². The compliances of the designs using the cell-specific interpolation functions (b and c) are much smaller than the design using generic (unphysical) power-law with p = 1 (d) and post-processed by interpreting gray elements with square- or X-shaped cells, demonstrating the significance of using an accurate material interpolation model in optimization. It simply does not make mechanical sense to use a non-physical macroscale model like SIMP with p = 1 to control the macroscopic density in multi-scale approaches. In (d), the post-process of replacing gray elements by square-shaped cells leads to a large compliance of c^s = 1017.03, since low-density square-shaped cells are weak in shear. The designs using cell-specific interpolations (b and c) have higher compliance values than design (e) using the power law with p = 3 (p starts from 1 and is increased by 0.2 every 50 iterations). This suggests the optimization using cell-specific interpolations gets stuck in local minima. The number of iterations for obtaining both (b) and (c) is 400. Design (b) has a density distribution close to black and white, while design (c) has large gray regions. X-Shaped cells have large shear moduli, making it economical to place gray elements in the middle of the beam which is predominately under shear stress. Rank-2 material achieves optimal stiffness (f) due to the alignment of material anisotropy with principal stress directions, whereas the fixed orientation of the square- and X-shaped lattice cells restricts the adaptation of material anisotropy.

The curves for ${E}_{ij}^{H}(\rho )$ for the square- and X-shaped cells are lower than the theoretical limit according to the Hashin–Shtrikman bounds (Hashin and Shtrikman 1963). This means a suboptimal use of material using cellular structures with a simple geometry. To improve the properties, Zhou and Li (2008) applied inverse homogenization (see Section 2.3) to design a series of unit cells with varying material fractions. Figure 9 (top) shows unit cells independently optimized for an increasing material fraction, using the code provided by Xia and Breitkopf (2015). These optimized cells have distinct topology, owing to the large solution space in topology optimization. It is also noticeable that neighboring cells, despite being close in terms of material fractions, may have poor connectivity across the interface (see Fig. 9 (top), between the first and second and between the third and fourth cells). The connectivity issue is also recognizable in hierarchical approaches (Rodrigues et al. 2002), to be discussed in Section 5.4. Note that this is less of a problem in cellular structures parameterized by thickness since the topology there is constant, as long as the thickness does not vanish to 0.

Fig. 9

Top: Illustration of poorly connected microstructures. The unit cells are individually optimized for maximum bulk modulus under linearly varying volume fractions from 30 to 50%, from left to right. Bottom: Optimized connectivity, achieved by maximizing the bulk modulus of extended domains covering adjacent cells (Garner et al. 2019). Reprinted from Garner et al. (2019) with permission from Elsevier

Zhou and Li (2008) proposed three methods to address the connectivity issue between topology optimized microstructures: kinematic connectivity constraint, pseudo load, and unified formulation with non-linear diffusion. In the first two, unit cells are optimized individually, while solid non-design regions or pseudo loads are prescribed along the domain boundary to stimulate connectivity. In the last, unit cells are optimized altogether, while a nonlinear diffusion term defined on the domain covering all cells is integrated in the objective function to penalize disconnection and suppress checkerboard patterns. Radman (2013a, 2013b) improved the computational efficiency of the unified formulation by successively optimizing new unit cells, while considering connection to already optimized cells. Du et al. (2018) proposed a physics-independent connectivity index, which measures the amount of overlap in adjacent cells across the shared interface. Garner et al. (2019) proposed to optimize the connectivity, quantified by the physical properties of interest (e.g., bulk modulus) of an extended domain that covers adjacent cells. Poor connectivity between adjacent cells leads to an inferior bulk modulus for the extended domain as a compound cell, and thus is effectively suppressed in optimization. Figure 9 (bottom) shows results from this compound formulation. Mapping these compatibility optimized cells to a compliance-minimized macroscale density distribution, and evaluating structural compliance using both full-scale and homogenization-based analyses, it was found that the discrepancy is small (a relative error of 2%). This small error is also attributed to the limited and slowly varying nature of the microstructure which satisfies basic assumptions of the homogenization-based multi-scale modelling. This method has been demonstrated for optimizing up to 100 varying unit cells in 2D. Optimizing a large number of unit cells in 3D can be computationally expensive. Cramer et al. (2016) proposed geometric interpolation to obtain transitioning microstructures between individually optimized unit cells. This geometric approach works for microstructures of similar topology.

Cellular structures, both assemblies of geometric primitives and topology optimized microstructures, are typically defined within a square or cube domain. Such a fixed design domain may limit the achievable properties (c.f., Träff et al. (2019)). These cellular structures may be square symmetric (2D) or cubic symmetric (3D) but not (necessarily) isotropic. If orthotropic directions are oriented along regular finite element grid directions and not along principal stress directions as recommended by Pedersen (1989), inferior stiffness may be the result.

5.3 Parameterized unit cell with multiple parameters

Functionally graded cells can be obtained by uniformly varying the thickness of the geometric primitives in the cell domain, e.g., the square- and X-shaped 2D lattices in the previous subsection. Consequently, each element in the macroscale optimization is assigned a single design variable, the material volume fraction, ρ(d), where d is the thickness. To enlarge the solution space, the parameterization of unit cells can be extended. The unit cell in Fig. 10 (middle) has two superimposed geometric patterns, i.e., an X shape and a plus shape, each with an independent thickness, leading to more variations in the attainable elasticity tensor, E^H(d, t). The number of design variables per element can be further increased, e.g., by assigning an individual thickness per geometric primitive, indicated by different colors in Fig. 10 (right). Graded lattice structures in 2D and 3D optimized with multiple design variables per macro element have been demonstrated by Wang (2018, 2020). Imediegwu et al. (2019) used a lattice cell with seven independent parameters for 3D optimization. As the number of independent parameters for describing the unit cell increases, a large number of numerical homogenization operations is required. White et al. (2019) proposed a neural network surrogate model to provide the mapping from parameters to resulting metamaterial properties. Zhu et al. (2017) pre-computed the space of attainable Young’s modulus and Poisson’s ratio for certain types of microstructures, and used this space as a constraint for macroscale topology optimization with microstructures. Zhu et al. (2019) presented an approach for concurrent optimization of the microscale material distribution and a macroscale mapping function which transforms periodic microstructures into graded ones. This was extended and demonstrated for 3D problems (Xue et al. 2020).

Fig. 10

Illustration of unit cells with 1, 2, and 4 independent parameters, from left to right

As the number of parameters increases, the potential to generate a variety in the anisotropic directions of the microstructures increases, and this, to some extent, allows an adaption of microstructural anisotropy to the local stress directions during optimization. To ultimately align the anisotropy of microstructures to stress directions, it is critical to release the rotational freedom in the unit cell parameterization, as pioneered by Bendsøe and Kikuchi (1988).

5.3.1 De-homogenization

A benefit of the approaches using parameterized unit cells with multiple parameters including rotation is that one can get optimized designs that perform very close to the theoretical limit. For example, let us consider a classical optimization example: the Michell cantilever clamped on the left with a unit load on 20% of the right boundary. Using the square microstructure with a rectangular hole by Bendsøe and Kikuchi (1988) shown in Fig. 3, and parameterizing the design domain by 80 × 40 bi-linear elements, we can obtain a compliance value of 58.35 for an upper bound on the material volume fraction of 0.4. This is very close to the value of 56.73 obtained using optimal rank-2 microstructures reported in Sigmund et al. (2016). Although the microstructures are on a single length scale, it is still difficult to interpret the spatially varying multi-scale design as a well-connected manufacturable design. A naive approach would be to enlarge each microstructure to the size of a bi-linear element and apply the appropriate orientation. As can be seen in Fig. 11 (left), this results in a disconnected design. The process of constructing connected and physically realizable designs from homogenization-based optimization is referred to as de-homogenization—a term coined by G. Allaire and colleagues.

Fig. 11

Single-scale interpretations of a multi-scale design for the Michell cantilever using a volume fraction of 0.4 and the microstructures by Bendsøe and Kikuchi (1988). Left: Naive interpretation, right: de-homogenized design using the approach presented in Groen and Sigmund (2018). Reprinted from Groen et al. (2019) with permission from John Wiley & Sons

A very promising post-processing method to obtain well-connected mono-scale designs from a spatially varying multi-scale design is the de-homogenization method introduced by Pantz and Trabelsi (2008). Using this method, an implicit geometry description is created to represent/enlarge the multi-scale design to a fine but realizable single length scale. In recent years, the interest in this approach has renewed, resulting in simplifications and improvements of the approach (Groen and Sigmund 2018; Allaire et al. 2018). The approach of Groen and Sigmund (2018), applied to the Michell cantilever, resulted in the high-resolution mono-scale design (1600 × 800 elements) shown in Fig. 11 (right). A minimum feature size is applied and resulting compliance is 59.55. The core idea of the approach is to create a set of smooth mapping functions ϕ_i(x) that convert from the global frame of reference x to the microstructure frame of reference y. In other words, we create a conformal-like map in a similar fashion as texture mapping (Lévy et al. 2002). The most important requirement is that the spatially varying microstructure orientation is smooth throughout the domain such that smooth mapping functions ϕ_i(x) can be generated. A visual overview of how the approach works is given in Fig. 12. Consider a coated domain (Fig. 12a) with a porous interior using the microstructure in Fig. 3 with μ₁ = μ₂ = 0.1, the corresponding microstructure orientation 𝜃 is shown in Fig. 12b. One of the two mapping functions ϕ₁(x) is shown in Fig. 12c and the corresponding de-homogenized design is seen in Fig. 12d. A main advantage of this approach is that since it uses an implicit geometry description, the design can be evaluated at an infinitely fine resolution. Furthermore, the periodicity (i.e., the amount of microstructures) can be explicitly controlled, even to have a more uniform spacing as can be seen in Fig. 12e (Groen et al. 2019). It has to be mentioned that different kinds of microstructures (i.e., with non-orthogonal features) can be de-homogenized as well (Geoffroy-Donders 2018; Groen 2019; Kumar and Suresh 2020; Tamijani et al. 2020).

Fig. 12

Example of the de-homogenization procedure. From left to right, a coated domain, b microstructure orientation 𝜃, c mapping function ϕ₁(x), d de-homogenized design, e de-homogenized design using adaptive periodicity approach from (Groen et al. 2019). Reprinted from Groen et al. (2019) with permission from Elsevier

A challenge of the de-homogenization approach is to allow a smooth microstructure orientation throughout the domain. This is especially an issue since the elasticity tensor has a rotational symmetry of π, i.e., rotating the microstructure 180° does not affect the properties; however, it complicates the de-homogenization. To circumvent this, Pantz and Trabelsi (2008) proposed a two-field approach to solve for ϕ_i(x). In a different approach, Groen and Sigmund (2018) used an image-based sorting approach, while Allaire et al. (2018) used a discontinuous Galerkin approach to solve for ϕ_i(x). Furthermore, it is known that singularities can occur in the design domain, either caused by the stress field or caused by regularization methods. Hence, the algorithms should be able to handle these singularities. Different approaches to do so are discussed in Pantz and Trabelsi (2010), Geoffroy-Donders (2018), and Stutz et al. (2020). The extension of the approach to 3D requires a solution to the problem that the principal stress directions are not well-ordered in 3D. To circumvent this, Geoffroy-Donders et al. (2020) introduced a regularization functional, while Groen et al. (2020) combined a regularization and an image-based sorting approach.

A further challenge is the de-homogenization of rank-N designs optimized for multiple load cases. The non-uniqueness of the optimal microstructures, which allows for an infinite amount of designs, has to be taken into account, since the de-homogenization procedure requires smooth vector fields throughout the design domain. Hence, making sure that the microstructure orientation is continuous throughout Ω, without restricting the performance is a key challenge (Groen 2019).

Above approaches seek a global parameterization of each of the principal stress directions. In a different approach, de-homogenization is cast as finding a quadrilateral (2D) or hexahedral (3D) mesh with each edge being aligned with the optimized stress directions (Wu et al. 2021), borrowing ideas from field-aligned meshing (Jakob et al. 2015; Gao et al. 2017). Instead of decomposing a tensor field into three vector fields, the difference between the orientation of a (locally parameterized) quad/hex element and a stress tensor is measured by comparing all possible perturbations of decomposed directions. The cumulative difference over all elements is minimized by an iterative local updating scheme. This approach is demonstrated in combination with a modified version of the rectangular hole model. As illustrated in Fig. 13a for 2D cases, the unit cell is allowed to rotate (𝜃) and elongate independently along each axis (α_x and α_y). In addition, another design variable (not shown in (a)), similarly used in PAMP (Liu et al. 2008), is introduced to encode whether or not a finite element is filled with lattices. The optimized lattice distribution and the compiled continuous lattice structure are shown in Fig. 13b and c, respectively. The optimized structure possesses spatial variations in orientation, porosity, and anisotropy. A 3D example is shown in Fig. 2c.

Fig. 13

a The unit cell, a modified version of the rectangular hole model, adapts by rotating and elongation. b The optimized lattice distribution for the Michell cantilever using 15% solid material. The frames are rotated and elongated according to the optimized fields. c A continuous lattice structure compiled from b, demonstrating spatial variations in orientation, porosity, and anisotropy (Wu et al. 2021). Ⓒ2021 IEEE. Reprinted, with permission, from Wu et al. (2021)

Yet another mapping strategy is to modify crossing of contour lines of aforementioned ϕ fields, connect these by truss or frame elements, and then tune nodal positions and connectivities in a subsequent shape optimization process (Larsen et al. 2018).

5.4 Restricted unit-cell design

Predefined cellular structures are beneficial for computational efficiency. However, they may significantly restrict the solution space. To access the full solution space, while avoiding intensive full-scale analysis, an idea is to use a hierarchical formulation. Hierarchical optimization of material and structure dates back to Rodrigues et al. (2002) which was later extended to 3D (Coelho et al. 2008), both using SIMP. Comparable results were reproduced with level sets (Sivapuram et al. 2016). Hierarchical optimization combined with inverse homogenization is also referred to as concurrent (or simultaneous) optimization of structure and material (or microstructure). Like in approaches in other categories, here, a two-scale discretization of the domain is employed. The formulation involves one problem at the global (or macro) scale and many problems at the local (or micro) scale. The global problem determines the macroscopic spatial distribution of homogenized material, and local problems determine microscopic spatial distribution of solid and void phases by optimizing for homogenized properties. The structural equilibrium in the macroscale is in general nonlinear due to the microstructure adaptation (Jog and Haber 1996). A nonlinear resolution framework based on FE² scheme was developed to address this nonlinearity (Xia and Breitkopf 2014). This subsection discusses approaches that prescribe the domain shape of microstructures and/or their orientation.

The microstructural optimization problem (inverse homogenization) should result in single-scale approximations of theoretically optimal rank-N composites. In a recent work, it was found that inverse homogenization at low volume fractions is prone to producing local optima, and a simple mapping approach was suggested to approximate rank-3 laminates (Träff et al. 2019). These mapped microstructures perform relatively close to theoretical energy bounds, and can serve as starting guesses for inverse homogenization problems to achieve performances even closer to the bounds.

In each iteration of a hierarchical solution process, following a solved global problem, the local problems become independent from each other. On the positive side, the independent problems can be solved in parallel by sending sets of local problems to different processors (Coelho et al. 2011). This gains a computational speedup and thus allows solving two-scale problems in 3D. On the other hand, the independent nature of the local problems creates a critical challenge regarding the compatibility of microstructures across the shared boundary. We emphasize that the problem of concern is related to structural properties beyond the disconnected geometry, and thus choose to use compatibility in lieu of connectivity. The compatibility problem arises since disconnections between adjacent microstructures are not captured in the global analysis using homogenized properties (separation of scales). A mechanical indication of compatibility is the discrepancy between the objective (e.g., compliance) evaluated by a full-scale analysis and by an analysis using the homogenized properties (Garner et al. 2019). Some approaches summarized in Section 5.2 for improving compatibility in functionally graded microstructures are also applicable to the hierarchical optimization. For instance, by using extended domains that overlap in local optimizations, the compatibility can be significantly improved, reducing the discrepancy in compliance values between full-scale and homogenization analyses from six orders of magnitude to two (Garner et al. 2019) (see Fig. 14). While this is good progress, a discrepancy of two orders of magnitude is still alarming. To examine the optimality, we can visually compare hierarchically optimized structures with those from full-scale approaches with local volume constraints (Section 3) and de-homogenization approaches (Section 5.3). For a single load, orthogonal microstructures that are individually aligned with principal stress directions are known to be close to optimal. Such orthogonal microstructures are distinct in full-scale approaches with local volume constraints as well as de-homogenization approaches, but are difficult to discover in large regions of hierarchically optimized structures. A reason behind the poor compatibility is that the fixed, axis-aligned rectangular domain used in inverse homogenization is incapable of accommodating rotation of these orthogonal microstructures.

Fig. 14

Illustrative result of hierarchical optimization of material and structure of a simply supported beam, without (left) and with (right) a compound formulation to improve microstructural compatibility (Garner et al. 2019). Half of the optimized structure from each formulation is shown, for compactness. The compliance of the left part by a full-scale analysis and an analysis using homogenized properties is 8.32 × 10¹³ and 6.54 × 10⁷, respectively, while the compliance of the right part is 4.96 × 10⁹ and 8.42 × 10⁷, respectively. Reprinted from Garner et al. (2019) with permission from Elsevier

Other strategies have been proposed to enhance connectivity. Wang et al. (2017) proposed a shape metamorphosis method based on level-set representations to interpolate a prototype microstructure to generate a family of graded microstructures. The interpolated microstructures are connectable to each other in a natural way since they present similar topological features and material distribution patterns at their edges. Li et al. (2018b) adopted the kinematic constraint approach (Zhou and Li 2008) for level-set–based topology optimization of functionally graded cellular composites hosting auxetic metamaterials. Zhou et al. (2019) proposed a geometric connectivity index upon which constraints were defined and included in the macroscale optimization to improve connectivity. Liu et al. (2020) recently proposed to ensure connectivity between any two types of the microstructures by introducing pre-defined connective regions in the microstructural unit cells and keeping these regions sharing the same topology. In these approaches, the connectivity was often visually assessed, and a mechanical assessment was unfortunately missing.

The compatibility issue is circumvented if the optimization problem is reformulated to design structures consisting of repetitive microstructures, at the cost of reduced structural performance. Such a formulation was presented by Fujii et al. (2001) for designing repetitive microstructures in the entire design space. Liu et al. (2008) incorporated a macroscale variable to concurrently design the microstructure and its distribution. Here, the solution space is reduced to a single microstructure, a strategy similar to pattern repetition in full-scale approaches (Section 3). The global analysis is performed using homogenized properties rather than on the full scale, assisted with an interpolation scheme called Porous Anisotropic Material with Penalization (PAMP). It was first demonstrated for compliance minimization, and extended for maximizing fundamental frequency (Niu et al. 2009), for considering load uncertainties (Guo et al. 2015) and more. Similar formulations and extensions using evolutionary procedures can be found in Huang et al. (2013) and more.

In between the spectrum from a single microstructure to the full solution space, approaches have been developed to design the structural layout of a few unique, unprescribed microstructures (Deng and Chen 2017; Zhang et al. 2018; Liu et al. 2020). Pizzolato et al. (2019) built on the PAMP framework (Liu et al. 2008) and developed a level-set approach to optimize the distribution of multiple concurrently optimized microstructures. It was studied in the context of heat transfer problems. When the location of unique microstructures is not prescribed, it is often cast as a multi-material optimization problem, for which discrete material optimization (Stegmann and Lund 2005) and ordered SIMP interpolation (Zuo and Saitou 2017) are applicable.

As discussed in Section 3, periodic and graded microstructures can also be designed using full-scale approaches. When the structural analysis is performed on the full resolution, a poor connectivity is reflected by a suboptimal objective. Thus, full-scale approaches naturally ensure good connectivity between microstructures or subdomains, at the price of intensive full-scale analyses. Therefore, results from relevant full-scale approaches may serve as a reference for multi-scale approaches of graded microstructures. From the comparison in Fig. 4, it is clear that a logical (but often underreported) consequence of using periodic microstructures with a fixed orientation is a large reduction in stiffness.

5.5 Unrestricted unit-cell design

All approaches in the previous section have one thing in common; the shape of the unit-cell domain is rectangular. However, it is very important to acknowledge the effect of the unit-cell domain on the performance. In a large amount of works that consider a single load case, the numerical examples show that the optimized unit cells resemble a rotated version of the microstructure by Bendsøe and Kikuchi (1988) (Fig. 3). However, since the unit-cell domain cannot be rotated, the shape of the microstructures is modified to account for a periodic shape, in turn reducing performance. In other words, if the unit-cell design domain was allowed to rotate, a simpler and close to optimal microstructure would have been found.

To investigate the effect of the unit-cell shape and orientation on the performance, Träff et al. (2019) did an extensive comparison of the effect of the starting guess, when inverse homogenization is used to optimize the microstructure performance. They observed that the material design problem is non-unique and highly non-convex, i.e., different starting guesses result in completely different unit-cell designs and performances. Furthermore, a starting guess of the unit-cell shape and topology based on an optimal rank-3 microstructure resulted in a significantly better performance compared to the use of a rectangular domain. The parallelogram shape of the unit cell, which was optimized as well, allowed for periodicity patterns that cannot be described by a rectangular domain as seen in Fig. 15. Especially for lower volume fractions, Träff et al. (2019) observed that a unit cell with a starting guess and shape based on an optimal rank-3 laminate significantly outperformed (e.g., up to 30% more efficient) the designs using a traditional rectangular domain. Hence, the choice of the unit-cell domain significantly influences the performance and researchers have to be aware of this when choosing a multi-scale topology optimization algorithm.

Fig. 15

Left: Rank-3 microstructure with indicated hierarchy. Right: approximated single-scale microstructure using the method by Träff et al. (2019) to indicate different periodicity patterns that cannot be achieved by a rectangular unit cell. Reprinted from Träff et al. (2019) with permission from Springer Nature

Recently, Wang et al. (2019) showed that near-optimal and periodic truss lattice structures could be obtained for multiple load cases by distorting simple Bravais-like lattice structures to a parallelepiped. Besides using a parallelogram in 2D or a parallelepiped in 3D, one can use many more different types of polygons to solve the homogenization equations (Barbarosie et al. 2017; Podestá et al. 2019). For example, in 2D, a hexagon can be used to describe a periodic isotropic hexagonal microstructure (Sigmund 2000).

From all studies discussed above, we can conclude that performing multi-scale optimization with unit cells that are optimized using inverse homogenization on a rectangular domain reduces the optimality. This effect has been acknowledged by Barbarosie and Toader (2014) who combined the optimization of the microstructure topology and shape. This approach can achieve structural designs that are close to optimal; however, an extensive comparison unfortunately has not been performed. To address the large computational cost included in solving the problem, the approach has been parallelized. To the authors’ knowledge (and surprise), this is the only work on multi-scale topology optimization that simultaneously addresses the macroscopic design and both unit-cell topology and shape to get as close to the optimal design that can be achieved on a single length scale. Nevertheless, this method does not address the connectivity of the spatially varying unit cells over the design domain. More research has to be done to efficiently blend unit cells of different shapes together without a loss in performance to further advance this method.

6 Discussion

6.1 Motivation

Two important questions that people working on multi-scale structures should ask themselves are:

• What are the benefits of multi-scale structures?

• What are the expected benefits of multi-scale approaches?

These two questions are related but apparently independent from each other. While multi-scale structures are most often optimized, not surprisingly, by multi-scale approaches, they can also be designed using full-scale approaches with some additional constraints (Section 3). One important motivation of multi-scale topology optimization is to accelerate the computation for optimizing structures at high resolution. Here, it shall be reminded that theoretical stiffness-optimal structures span multiple scales. A higher resolution discretization enables the appearance of fine geometric details, which may bring the performance of optimized structures closer to theoretical limit. This, however, implies significant computational cost. Thus, multi-scale approaches are introduced for means of acceleration. These include the original hierarchical approach by Rodrigues et al. (2002), de-homogenization approaches (Pantz and Trabelsi 2008; Groen and Sigmund 2018), and works along these directions.

Many multi-scale approaches, while being formulated as an optimization problem (e.g., for maximizing stiffness), restrict the solution space by enforcing one or a few repetitive microstructures that are either predefined or concurrently optimized, with a fixed orientation and/or limited density range. This restriction may be beneficial under various considerations, from structural such as buckling strength (Clausen et al. 2016), robustness (Wu et al. 2018), and multi-functionality to operational such as inspection and repair, and from manufacturability over aesthetics to sustainability. In this regard, the microstructures play the role of accounting for these considerations. From a mathematical perspective, it would be interesting to explicitly model these requirements and integrate them in an optimization problem with a high-resolution discretization. This is challenging, since a mathematical model of some of the requirements is not available yet or comes with computational complications. Thus, restricting the solution space by a reduced parameterization can be understood as a strategy to balance the defined objective and various other considerations. A general recommendation is to always motivate use of stiffness suboptimal microstructures. Too many works seem to forget this aspect and show optimized designs that are clearly not optimal with respect to stiffness.

Analogous to technical optimization, nature has been an advocate for multi-scale structures. Hierarchically organized, functionally graded structures can be found in plant and animal bodies such as bamboo and bone (Lakes 1993; Fratzl and Weinkamer 2007). These bio-structures are remarkable from a mechanical perspective while supporting biological functionalities.

6.2 Evaluation

Most (if not all) multi-scale approaches make use of homogenization, which assumes separation of scales, i.e., microstructure should be much smaller than the macrostructure. This assumption often becomes invalid when considering the finite resolution of manufacturing processes. A general rule of thumb is that cells should be repeated 5 to 10 times before effective properties can be trusted. This certainly means that commonly seen approaches where one macroscale finite element corresponds to one microstructure cell, where cells may vary between each element, should be used with extreme caution. Interestingly, however, multi-scale optimization, where microstructure is appropriately adapted to local stress fields, unimpeded by cell geometry or orientation (c.f., de-homogenization approaches) provides extremely good performances even for quite large periodicities, i.e., with lack of scale separation (see Groen and Sigmund (2018) and Wu et al. (2021)). The reason is that the optimal microstructures ensure purely tension-compression-dominated deformations at both micro- and macroscale, which homogenization-wise need fewer cell repetitions for accuracy than required for bending or shear dominated deformations. Under all circumstances, it is strongly recommended to verify any assumptions of scale separation with subsequent full-scale analyses. Full-scale analysis of course implies heavy computation. Fortunately, however, recent adoption of advanced linear solvers and parallel computing has partially alleviated this problem (Amir et al. 2014; Aage et al. 2015; Wu et al. 2016a). Validation by full-scale analysis in 2D starts to appear in a handful of recent papers (Groen and Sigmund 2018; Garner et al. 2019; Wu et al. 2021), and should be an integral part of the validation of all multi-scale works.

Optimization in general and multi-scale approaches in particular are about making a delicate balance between conflicting objectives and constraints, and between the quality of results and computational efficiency. In this regard, gains in one aspect are often, understandably, accompanied by losses in other aspects. Both sides are valuable for inspiring future development. These are especially important for researchers who may be less familiar with the topic, e.g., students, application engineers, and colleagues from other disciplines. Therefore, we strongly recommend to include a quantitative comparison (or a discussion if a quantitative analysis is difficult to perform, e.g., regarding aesthetics) when a new approach is introduced. Comparisons can be made on different levels, with non-optimized designs, with optimized mono-scale structures, with designs from alternative multi-scale approaches, etc. In the case of a 2D compliance design subject to a single load case, we recommend to perform a quantitative comparison with what is theoretically possible using rank-2 microstructures. To do this, we provide a Matlab script in Appendix.

6.3 Extensions, alternative formulations, open questions

Our review has been focusing on design parameterizations, which are typically demonstrated in compliance minimization under the assumption of small deformations. Accommodating stress constraints, buckling constraints, and geometric and material nonlinearities represents an important and non-trivial next step. Some recent works have started to tackle these challenges.

• Stress constraints. Collet et al. (2018) proposed a formulation for optimizing periodic microstructures with stress constraints, assuming that the classical von Mises stress criterion remains valid at the microscale. Ferrer et al. (2020) proposed a square cell with a super-ellipsoidal hole, avoiding right angles and thus stress concentration in the classic rectangular hole model. In concurrent optimization of shell and parameterized lattice infill, Yu et al. (2020) proposed using two stress constraints on the macroscale, a von Mises stress-based constraint for the solid shell layer, and a Tsai-Hill yield criteria-based constraint for the homogenized infill.

• Buckling, geometric and material nonlinearities, and plasticity. Yuge and Kikuchi (1995) developed elasto-plastic finite element analysis based on homogenization to optimize frame structures subjected to plastic deformation. Neves et al. (2002) addressed the problem of determining highly localized buckling modes in perfectly periodic cellular microstructures of infinite extent. Thomsen et al. (2018) proposed a optimization model to predict local and global microstructural buckling, based on which periodic cellular materials are optimized for maximized strength under compressive load. A systematic investigation of the performances of simple and optimized periodic infill structures in terms of finite scale stiffness and buckling was presented by Wang and Sigmund (2020). Strength of common 3D microstructures is treated in Andersen et al. (2021). Furthermore, Bluhm et al. (2020) proposed a framework for benchmarking the ability of periodic microstructures to maintain stiffness under large deformations, accounting in a unified manner both for buckling and softening due to geometric and material nonlinearities.

Going beyond structural mechanics, design of multi-scale structures involving other physics and even multi-physics is another important venue to explore. Some of the works along this direction have been mentioned in previous sections when specific multi-scale approaches were introduced. For instance, Das and Sutradhar (2020) presented an extension of the full-scale approach with the local volume constraints (Wu et al. 2018) to optimize heat-dissipating structures considering structural and thermal performance. Fluid flow through the pores in lattice structures is another interesting topic, relevant for design of actuators (Andreasen and Sigmund 2011) and biomedical implants (Challis et al. 2012). Furthermore, as discussed in Section 2.3, inverse homogenization has been used to design metamaterials with extreme or counter-intuitive physical properties such as a negative Poisson’s ratio and negative thermal expansion (Sigmund and Torquato 1996). Designing multi-scale structures composed of such metamaterials is fascinating. It can potentially open innovative application areas in shape-morphing products, soft robots, and 4D printing (i.e., with an extra dimension of transformation over time).

Optimized structures from multi-scale approaches typically have two scales. It is worth mentioning that buckling-enhanced microstructures themselves exhibit two (or more) hierarchical scales (Thomsen et al. 2018). Integrating them in macroscale structural optimization would effectively lead to three-scale structures. The microstructures in two-scale approaches are normally defined in a uniform discretization of the macro domain, and thus have the same spatial extent. Wu et al. (2016b) proposed a full-scale approach to design octree-tree like structures, where the cells exhibit spatially varying sizes. This allows to achieve a large range of variations in porosity and pore sizes. The discrete problem of hierarchical subdivision was later reformulated by continuous variables to facilitate gradient-based optimization (Wu 2018).

Microstructures are typically defined and analyzed based on a Cartesian grid discretization with a fixed orientation. Approximating free-form macroscale shape by square or cubic microstructures leads to staircasing on the boundaries. This artifact may be less of a concern when the macroscale domain is orders of magnitude larger than the microstructure size. However, with the limited resolution of manufacturing processes, the effects of this artifact may be non-negligible, both visually and mechanically. Conforming quad/hex meshing can be an alternative for constructing a boundary-aligned grid (Wu et al. 2021). However, mapping square or cubic microstructures into irregular grids may introduce error in mechanical properties. Note that homogenization assumes a repeatable domain, and thus direct homogenization over irregular cells seems not possible.

As stated in the Introduction section, much of the recent interest in designing multi-scale structures is triggered by contemporary advances in additive manufacturing (AM). While being able to produce highly complex shapes and even structural features spanning seven orders of magnitude (Zheng et al. 2016), existing AM processes are not free from manufacturability and post-processing requirements. Structural design thus needs to consider, e.g., a minimum feature size, self-supportingness (i.e., free of critical overhang), accessibility for removing unsolidified powder or resin, and auxiliary support structure (Liu et al. 2018a). In the context of multi-scale design, these requirements (or some of them) can be satisfied by carefully selecting unit cells, e.g., beam-like lattice cells (see prints in Fig. 2) or self-supporting rhombic cells (Wu et al. 2016b). Stiffness optimal rank-N laminates and closed-walled cells with stiffness close to theoretical bounds (Sigmund et al. 2016) may, depending on process, be less favorable than open-walled cells regarding the removal of unsolidified powder or resin. Along with developments of advanced manufacturing technologies, design for manufacturability will continue to be an important topic in (multi-scale) topology optimization.

The design of multi-scale structures is a topic of interest in multiple other disciplines apart from applied mathematics, computational mechanics, and mechanical engineering. We envisage that cross-pollination with material science (e.g., multi-scale materials modelling (van der Giessen et al. 2020)) as well as geometry modelling and processing (e.g., texture synthesis (Dumas et al. 2015)) will further advance optimization and inverse design of multi-scale structures. The cross-pollination is actually already benefiting the field. For instance, the conformal-like mapping and field-aligned meshing have been used for de-homogenization.

7 Conclusion

We have reviewed the development of multi-scale topology optimization, from its inception to the current state-of-the-art. Existing approaches are classified into full-scale and multi-scale approaches, according to whether or not the separation of length scales is assumed in the modelling. Full-scale approaches control structural details by imposing, e.g., pattern repetition and local volume constraints, while performing structural analysis on the fine scale, which is computation-intensive. Multi-scale approaches reduce computation by using analytical or numerical homogenization, under the assumption of separation of length scales, which may give rise to compatibility issues. These approaches are categorized in this review based on the parameterization of the micro- and macroscales.

The development in this field is exciting. Innovative approaches and applications continue to appear. To construct an objective understanding of multi-scale structures and multi-scale optimization approaches, we make the following recommendations:

• Since the assumption of separation of length scales may not be respected in the optimized multi-scale structures, structural performance evaluated by homogenization may not faithfully represent reality. It is thus strongly recommended to compare with full-scale analysis when homogenization is used in the optimization approach.

• Many parameterizations have been introduced for designing multi-scale structures while reducing the gap between full-scale and homogenization-based analyses. These parameterizations may (severely) limit the achievable structural objective. It is thus also recommended to compare with standard mono-scale approaches under a comparable computation time, when investigating homogenization-based methods.

• Multi-scale structures hold the promise of achieving superior performance while being intrinsically lightweight, robust, and multi-functional. The true benefits of novel multi-scale structures need to be validated, numerically, and/or experimentally.

Appendix: MATLAB code topRank2

Together with this review, we present a MATLAB code for the topology optimization of 2D structures subject to a single load case using optimal rank-2 microstructures consisting of solid material (with stiffness E⁺) and void. The code is based on the MATLAB code by Andreassen et al. (2011), and we will therefore only discuss the differences. The most obvious difference is the use of a rank-2 material model (see Fig. 5 for the visualization of a unit cell). The effective material property can be analytically derived using the homogenization equation (see, e.g., Allaire (2002) and Bendsøe and Sigmund (2004)), and can be summarized in Voigt notation as,

$$ \begin{array}{@{}rcl@{}} \ {E}^{H}(\mu_{1},\mu_{2},\theta) & =& \frac{E^{-}}{1-\nu^{2}} \left[\begin{array}{ccc} 1& \nu & 0 \\ \nu& 1& 0 \\ 0 & 0 & \frac{1-\nu} {2} \end{array}\right] + \frac{E^{+}}{1-\mu_{2}+\mu_{1}\mu_{2}(1-\nu)} \\ &&\textbf{R}^{T}(\theta) \left[\begin{array}{ccc} \mu_{1} & \mu_{1}\mu_{2}\nu & 0 \\ \mu_{1}\mu_{2}\nu& \mu_{2}(1-\mu_{2}+\mu_{1}\mu_{2})& 0 \\ 0 & 0 & 0 \end{array}\right] \textbf{R}(\theta) \end{array} $$

(5)

Since a rank-2 microstructure consisting of solid material and void contains no stiffness against shearing, a small isotropic background stiffness E⁻ is added to make this material model work stable on a finite mesh. To minimize the effect of this background stiffness and avoid getting stuck in local minima, we start with E⁻ = 0.1E⁺ and gradually reduce every 50 iterations until E⁻ = 10^− 6E⁺. For efficient assembly of the stiffness matrix, we pre-integrated 6 element matrices, i.e., one for each unique index of the rotated elasticity tensor. As is discussed by Pedersen (1989), a microstructure is optimally aligned with the principal stresses/strains. Therefore, we update the microstructure orientation (𝜃) during each design iteration based on the principal stress directions. Subsequently, we update the relative layer widths (μ_i) based on the gradients using the optimality criterion approach. Finally, it should be noted that we use the density filter to avoid the formation of checkerboard-like patterns that are artificially stiff. Similar to the 88-line MATLAB code, the default design problem is the half MBB-beam, where the compliance is minimized subject to an upper bound on the volume fraction of the stiff material. The code can be called as follows:

toprank2(nelX,nelY,rMin,volFrac)

where nelX and nelY are the number of bi-linear elements in x- and y-direction, respectively, rMin is the filter radius in element length h used for the density filter, and volFrac is the volume fraction that the stiff material is allowed to use.