Available online at www.sciencedirect.

com
ScienceDirect

Comput. Methods Appl. Mech. Engrg. 415 (2023) 116209

www.elsevier.com/locate/cma

Multiscale concurrent topology optimization of hierarchal

multi-morphology lattice structures
Xiliang Liu, Liang Gao, Mi Xiao ∗
National Center of Technology Innovation for Intelligent Design and Numerical Control, State Key Laboratory of Intelligent Manufacturing
Equipment and Technology, Huazhong University of Science and Technology, Wuhan 430074, China
Received 9 April 2023; received in revised form 23 May 2023; accepted 23 June 2023
Available online 14 July 2023

Abstract
This paper proposes a multiscale concurrent topology optimization method for design of hierarchal multi-morphology
lattice structures (HMMLSs), which features in the Kriging metamodel assisted Uniform Multiphase Materials Interpolation
(KUMMI) model and the sigmoid function (SF) based hybrid transition strategy. Specifically, level set functions are adopted
to model hierarchal lattice unit cells (LUCs). Then LUCs with different morphology are treated as distinct lattice materials
(LMs) and the SF-based hybrid transition strategy is employed to realize smooth transition between multi-morphology LUCs.
Besides, Kriging metamodels are constructed to predict the equivalent properties of LUCs efficiently so that the topology
optimization of HMMLSs can be achieved with an affordable computational cost. Two kinds of design variables are involved
in the topology optimization of HMMLSs, where discrete variables indicate the categories of lattice materials and continuous
variables determine the relative densities of lattice microstructures. The proposed KUMMI model can couple the two kinds
of design variables dexterously to realize concurrent optimization of relative densities of LUCs, and distribution regions and
percentages of different LMs within HMMLSs. Several numerical examples are presented to demonstrate the effectiveness
and applicability of the proposed method. Results indicate the optimized HMMLSs exhibit rational distribution of hierarchal
multi-morphology LUCs, hence achieving superior structural performance.
© 2023 Elsevier B.V. All rights reserved.

Keywords: Hierarchal multi-morphology lattice structures; Multiscale topology optimization; Concurrent design; Connectivity; Kriging metamodel

1. Introduction
Lattice structures are becoming increasingly popular in aerospace, automotive, and biomedicine areas due to
their lightweight design schemes [1–4]. Generally, the properties of lattice structures can be altered by varying the
configurations and porosity of lattice microstructures. For instance, lattice structures with controlled anisotropy can
be obtained by combining different categories of lattice microstructures [5], the plastic energy absorption ability
of lattice structures will be strengthened by introducing multi-morphology lattice microstructures [6], improved
structural stiffness is observed on the lattice structures with spatially varying porosity [7], high-efficiency heat
dissipation structures are designed through changing the configurations and porosity of lattice microstructures
∗ Corresponding author.
E-mail address: [email protected] (M. Xiao).

https://doi.org/10.1016/j.cma.2023.116209
0045-7825/© 2023 Elsevier B.V. All rights reserved.
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

concurrently [8]. Therefore, introducing hierarchical multi-morphology lattice structures (HMMLSs) is essential
and applicable for improving the mechanical performance of lattice structures. Additionally, recent advancements in
additive manufacturing (AM) technology have greatly facilitated the manufacturing of HMMLSs, making it possible
to employ these structures with complex geometric morphology to achieve extraordinary performance.
Different from heuristic design and bionics design, topology optimization is a generative structural design
method that aims at achieving the optimal material distribution within the given design domain [9]. Earlier
attempts at designing hierarchical structures using topology optimization can be traced back to the work of
Bendsœet al. [10], who designed composites using the homogenization method (HM). Then the classical “Solid
Isotropic Material with Penalization” (SIMP) topology optimization method was proposed to avoid the complexity
and heavy computational burden of HM-based topology optimization [11,12]. As researchers explored the physical
significance of gray elements that appeared in the optimization process, multiscale topology optimization drew
increasing attention [13–16].
Multiscale topology optimization aspires to pursue extraordinary structural performance by optimizing the
configurations of microstructures and their macro periodic distribution pattern. Rodrigues et al. [17] accomplished
pioneering work in this area by introducing a hierarchical optimization model to concurrently optimize material
and structure. However, the heavy computational cost of generating point-by-point varied microstructures hinders
its widespread application, even though these designs are more theoretically optimal. To address this issue, several
improved multiscale topology optimization approaches have been proposed to reduce the computational burden. For
example, Liu et al. [18] proposed a concurrent topology optimization method for minimum structural compliance
problems, in which both the macrostructure and its uniformly distributed microstructures are optimized simulta-
neously. Yan et al. [19] developed Bi-directional evolutionary structural optimization (BESO) based concurrent
design frame for designing composite microstructures. Sivapuram et al. [20] achieved material and structural
optimization simultaneously by predefining the positions of all microstructures. Xu et al. [21] presented two-
scale concurrent topology optimization with multiple hierarchal microstructures by introducing misplaced material
volume constraints. Zhang et al. [22] adopted the Kriging metamodel to reduce the computational burden of
estimating microstructures in performing multiscale topology optimization. These improvement measures can be
roughly divided into two categories: material microstructures are assumed to be identical throughout the entire
macrostructure, or the macrostructure is composed of several different kinds of microstructures. Although these
methods have brought some relief to the computational burden of estimating microstructural properties in multiscale
topology optimization, there are still two important issues that need to be addressed. Firstly, the manufacturability
of multiscale structures should be taken into consideration during the optimization process. Secondly, there is a
need to strike a balance between the design space and computational efficiency to achieve satisfactory results.
Based on the above considerations, the multiscale design of hierarchical lattice structures has gained increasing
attention due to its manufacturability and high efficiency in properties estimation during the iteration procedure
in recent years. What’s more, hierarchical lattice microstructures are appropriate for executing parameterization,
interpolation, or function fitting for high-efficiency properties estimation during the iteration procedure. For example,
Cheng et al. [7,23] validated the structural manufacturability of hierarchical lattice structures through feature
evolution. Li et al. [24,25] employed the density mapping and interpolation approach for designing functionally
graded gyroid lattice structures. Wang et al. [26] realized a hierarchical lattice material infilled hip implant
design using proportional topology optimization. Wu et al. [27] presented the multiscale optimization method
to solve the size distribution problem of conformal lattice materials by the extended multiscale finite element
method. Murphy et al. [28] proposed a robust three-dimensional multiscale structural optimization framework
that reduced the computational expense of optimization by parametrizing the lattice geometry. Ngoc et al. [29]
developed an adaptive geometric components technique to explore coated structures with non-homogeneous material
microstructures under buckling criteria, where the design layout, coated layer, and the interior functionally graded
microstructures were concurrently optimized. Lu et al. [30] proposed a novel approach for designing material
microstructure in mechanical cloaking by formulating a concurrent multiscale optimization model of the structural
topology and material properties of the cloaking devices. Banh et al. [31,32] performed the density-driven topology
optimization of functionally graded structures considering structural static and dynamic performance by proposing
a generally refined interpolation scheme. Rastegarzadeh et al. [33] employed the neural network as the optimizer
in the topology optimization of triply periodic minimal surfaces (TPMS) based hierarchal lattice structures, and
the sharp transition problem of adjacent microstructures is addressed without increasing computational cost. Wang
2
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

et al. [34] accomplished the multiscale design of graded lattice sandwich structures by human-aided isogeometric
topology optimization, where editable geometric models of optimized results can be automatically generated.
Compared to traditional multiscale topology optimization, these methods have the advantage of significantly
reducing computational costs and better manufacturability. However, one of the trade-offs for this improvement
is the limited design freedom that these methods offer.
To further enhance structural performance, the multi-morphology lattice unit cells (LUCs) infill design strategy
has been introduced into the multiscale topology optimization of lattice structures. For example, Imediegwu et al.
[35] ensured the connectivity of neighboring microstructures by smoothly changing the geometric parameterization
of hybrid hierarchical lattice structures. Similarly, Wang et al. [36,37] developed a Parameterized Interpolation for
Lattice Material (PILM) model and introduced two new design variables to accomplish the topology optimization of
hierarchical structures filled with non-uniform parameterized lattice microstructures. However, the parameterization
in these works is only applicable to truss-based LUCs and the design space is limited since only one LUC is
parameterized before the optimization. In addition to the parameterization method, heuristic methods have been
attempted to select the morphology of LUCs in designing lattice structures. Kang et al. [38] employed two types of
truss-based lattice microstructures in designing sandwich structures where the relative density criterion and analyzing
geometry were adopted to determine the type of lattice within the macrostructure. Ren et al. [39] established the
correlation mapping between principal stress direction and lattice type to optimize multi-scale and Multi-TPMS
lattices with geometric continuity. Nevertheless, selecting several categories of LUCs in advance and then applying
the heuristic criteria to fill the structure with them will loss of computational efficiency and design potential. To
automatically choose multi-morphology LUCs, Liu et al. [40] proposed an Approximation of Reduced Substructure
with Penalization (ARSP) model for the design of 2D HMMLSs, which considered two different yet connected
scales concurrently by assuming each truss lattice microstructure as a substructure in the optimization process. Wang
et al. [41,42] developed latent-variable Gaussian process models for truss-based LUC libraries, which embedded
both qualitative and quantitative design variables into a latent space and enabled the design of multiscale structures
with multiclass microstructures infill. Unfortunately, sufficient microstructures are needed to construct the latent
space. Selecting multiclass microstructures and obtaining necessary sensitivity information are challenging.
In existing works, the applicable morphology number and categories of LUCs are restricted, which limits the
potential in achieving superior structural performance of lattice structures. To address this limitation, this paper
proposes a multiscale concurrent topology optimization method for the design of lattice structures with hierarchical
multi-morphology LUCs. In particular, different categories of hierarchical LUCs are treated as disparate lattice
materials (LMs), even though they are all composed of the same base material. Then the Kriging metamodel
assisted Uniform Multiphase Materials Interpolation (KUMMI) model is proposed to characterize the interpolated
mechanical properties of hierarchical multi-morphology LUCs. In the KUMMI model, the quantitative design
variables are employed to associate LUCs’ equivalent elasticity tensor with their relative densities, while the
qualitative design variables determine the distribution domain of different LMs. Based on the KUMMI model,
the multiscale concurrent topology optimization of HMMLSs is then conducted, involving concurrently optimizing
the distribution patterns of distinct categories of LMs, the relative density of LUCs, and the percentage of different
categories of LMs within the HMMLSs. Additionally, to guarantee the connectivity of neighbor multi-morphology
LUCs, the level set functions are employed for the configuration description of LUCs, and the sigmoid function
(SF) based hybrid transition strategy is proposed to realize smooth transition between adjacent LUCs. Design cases
with different combinations of candidate LMs and different boundary conditions are tested, and the optimized LUCs
within HMMLSs performed satisfactory matching with the stress distribution of the design domain, demonstrating
the generality and effectiveness of the proposed method. Although the Uniform Multiphase Materials Interpolation
(UMMI) model [43,44] is employed in this paper, it’s believed that other kinds of multiphase material interpolation
methods are also worth trying in designing HMMLSs [45–52].

2. Generation of hierarchal multi-morphology lattice structures

As shown in Fig. 1, the considered design domain is filled with several categories of hierarchical LUCs
with various morphology. However, poor connectivity between neighbor multi-morphology LUCs will harm the
structural bearing capacity. Though several treatment methods have been proposed to guarantee the connectivity
between multi-morphology LUCs [33,39,53,54], they are confined to several TPMS-based LUCs which lose
generality in designing HMMLSs. In this section, the configuration interpolating method is employed to generate
hierarchical LUCs, and the SF-based hybrid transition strategy is proposed to guarantee smooth connectivity between
incompatible neighboring multi-morphology LUCs.
3
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Fig. 1. Illustration of the hierarchal multi-morphology lattice structures

2.1. Configuration mapping of hierarchical LUCs

For each LUC, its structural boundaries can be defined by a higher-dimensional level set function φ (x), as
illustrated in Fig. 2. Then the geometry model of a LUC can be represented by:
⎨φ (x) > 0, (Solid)
⎧
∀x ∈ Ω \∂Ω
φ (x) = 0, ∀x ∈ ∂Ω (Boundary) (1)
φ (x) < 0, ∀x ∈ D\Ω (Void)
⎩

where x is any point in the fixed design domain D, and ∂Ω are the boundaries of the solid structural domain Ω .
Considering the relative densities of hierarchical LUCs are graded varying, a configuration interpolating method is
employed to generate hierarchical multi-morphology LUCs, which can be stated as follows:
φmn (x) = φmpr o (x) − ϕmn , (m = 1, 2, . . . , M; n = 1, 2, . . . , N ) (2)
where M is the number of categories for hierarchical LUCs (i.e. the number of LMs). N represents the number
of finite elements within the design domain of HMMLS. φmn (x) represents the interpolated level set function of
pr o
the mth category of hierarchical LUCs, and the solid region of a LUC is defined by {φmn (x) > 0, ∀x ∈ D}. φm (x)
denotes the level set function of the mth prototype LUC, and there are M prototype LUCs for generating HMMLS.
ϕmn is the interpolation coefficient which( prcan )be calculated by (theprbisection algorithm based on the category and
o o
relative density of the LUC, and min φm (x) ≤ ϕmn ≤ max φm (x) . Besides, for the nth LUC generated by
)
level set function, its relative density can be computed by:
∫
1
ρmn = H(φmn (x))dD (3)
|D| D
where H(φmn (x)) is a Heaviside function used to indicate the solid part of the LUC which can be stated as:
1 i f φmn (x) ≥ 0
{
H(φmn (x)) = (4)
0 i f φmn (x) < 0
D and |D| are the design domain of a LUC and the area or volume of the design domain, respectively. In
addition, dozens of LUCs are employed in this paper for the topology optimization of HMMLS, e.g., Body-centered
cubic, Face-centered cubic, Schoen-Gyroid Sheet-based, Schwarz-Primitive, et al. and their detailed explicit level
set functions can be obtained in [55,56].

2.2. Sigmoid function based hybrid transition strategy for multi-morphology LUCs

The core idea of the SF-based hybrid transition strategy is introducing narrow transition regions that act as bridges
connecting different categories of LUCs. Within these transition regions, structures are obtained by hybridizing
4
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Fig. 2. The generation of lattice unit cell: (a) the level set function of a lattice unit cell and its zero cutting plane; (b) solid structural
boundaries of a lattice unit cell on the zeros cutting plane.

Fig. 3. Schematic of the sigmoid function.

neighbor LUCs to ensure smooth transition between multi-morphology LUCs, which is crucial for maintaining the
load-carrying properties of HMMLSs. Considering both the truss-based and TPMS-based LUCs are employed in this
paper, different processing strategy is implemented for these LUCs. Generally, for HMMLSs filled with truss-based
LUCs, no hybridization is needed as they can connect naturally, thus the SF-based hybrid transition strategy is
adopted when the TPMS-based LUCs appeared in the optimized HMMLSs.
In order to determine the topology configuration of hybridized structures within the transition regions, the SF is
employed to hybridize different level set functions of neighbor multi-morphology LUCs which can be stated as:
φnhyb (x) = αn (x)φm 1 n (x) + (1 − αn (x)) φm 2 n (x), (m 1 ̸= m 2 , m 1 = 1, 2, . . . , M; m 2 = 1, 2, . . . , M; ) (5)
hyb
where φn (x) is the hybridized level set function of the nth LUC within HMMLS, and the superscript hyb stands
for the LUC is hybridized. m 1 and m 2 are category numbers used to distinguish different categories of LUCs. αn (x)
is a weight function that satisfies αn (x) ∈ [0, 1], and the SF (shown in Fig. 3) is chosen as αn (x) here to achieve a
smooth transition between different categories of LUCs, which can be stated as:
1
αn (x) = (6)
1 + e−kG n (x)
where parameter k (k > 0) determines the curvature of SF as shown in Fig. 3(a), and the widths of the transition
regions can be altered by choosing different k. G n (x) defines the location of the transition boundary among different
categories of LUCs, e.g. G n (x) ≡ x = 0 determines the transition boundary in the 2D case shown in Fig. 3(b).
5
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Fig. 4. Illustration of SF-based hybrid transition strategy (LM. 1: Schoen Gyroid Sheet-based lattice material, LM. 2: Schwarz Primitive
lattice material, LM. 3: Cubic truss lattice material).

Considering the connectivity of 2D LUCs is easy to guarantee, the SF-based hybrid transition strategy is only
implemented in the 3D HMMLSs in this paper. The first step of implementing the SF-based hybrid transition
strategy is to detect the shared boundaries of different categories of LUCs. Specifically, for 3D HMMLS, the LUCs
are detected from up-to-bottom, left-to-right, and front-to-back to find the shared boundaries. What’s more, for each
3D LUC within HMMLS, only its right face, bottom face, and front face are considered in detecting the shared
boundaries for simplicity. Note that the central LUC in Fig. 4(c) needs to be hybridized twice to better connect its
right and bottom LUCs, and twice hybridizations are executed as a result. Then the level set function of this LUC
can be written as:

φnhyb (x) = αn1 (x) αn2 (x)φm 1 n (x) + 1 − αn2 (x) φm 2 n (x) + 1 − αn1 (x) φm 3 n (x), (m 1 ̸= m 2 , m 1 ̸= m 3 )
( ( ) ) ( )
(7)
6
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

1
αn1 (x) = 1
(8)
1+ e−kG n (x)
1
αn2 (x) = 2
(9)
1 + e−kG n (x)
where G 1n (x) and G 2n (x) are the two different transition boundaries within the LUC, and αn1 (x) and αn2 (x) are the
corresponding SF respectively.
Fig. 4 illustrates the effect of the SF-based hybrid transition strategy for multi-morphology LUCs. By comparing
Fig. 4(a) and (b), (c) and (d), it can be observed that the proposed strategy significantly changes the LUCs’
topology configuration and alleviates the sharp transitions within the multi-morphology lattices. Therefore, the
proposed strategy is a convenient, effective means to guarantee smooth transition between multi-morphology LUCs
in designing HMMLSs. Note that the level set function of LUCs should be normalized before the hybridization to
achieve better beautiful hybridization effect. Furthermore, the SF-based hybrid transition strategy can also be applied
to realizing the hybridization of 3D LUCs in more complex scenarios via recursive operations, similar to Eqs. (5)
and (7). As for the ‘stair-stepped’ material interfaces between adjacent hierarchical LUCs, the Kriging-assisted
morphological post-process method is employed to realize smoothly-varying hierarchical LUCs [57,58].
In addition, the categories of LUCs play an important role in the hierarchal lattice structure design, and many
pieces of research about the performance of the TPMS-based LUCs and truss-based LUCs show that they have
their particular strengths and weaknesses in specific situations [25,59,60]. Considering the SF-based hybrid transition
strategy can bring favorable improvement in the multi-morphology lattice structural topological continuity as shown
in Fig. 4 (a) and (c), hence TPMS-based LUCs and truss-based LUCs are suggested to be employed together in the
topology optimization of HMMLSs to gives full play to their strengths, resulting in excellent structural performance.
To the best of our knowledge, the TPMS-based LUCs and truss-based LUCs have rarely appeared in HMMLSs
simultaneously in previous.

3. Mechanical properties estimation of hierarchal multi-morphology LUCs

3.1. Homogenization method for properties calculation of hierarchal multi-morphology LUCs

In multiscale concurrent topology optimization of HMMLSs, it’s vital to quantitatively analyze the mechanical
properties of hierarchical multi-morphology LUCs, which is the basis for constructing the topology optimization
model.
In this paper, it is assumed that the base material used to construct multi-morphology LUCs has linear elastic
constitutive behavior. Then the equivalent properties of LUCs can be estimated by the HM by imposing periodic
boundary conditions [61,62]. The equivalent elasticity tensor of LUCs can be calculated by:
∫
1
ε − ε∗pq u(in j) D pqr s εr0(kl) − εr∗s u(kl) H (φn (x)) dΩn
( H ) ( 0(i j) ( )) ( ( ))
Di jkl n = (10)
|Ωn | Ωn pq s n

( )
where DiHjkl is the equivalent elasticity tensor of the nth LUC within the HMMLSs. Ωn and |Ωn | denote the
n
(i j)
( )
structural domain and the area or volume of n-th LUC, respectively. ε∗pq un stands for the varying microscale
0(i j)
strain field corresponding to the given macroscopic strain field ε pq with subscripts i, j, k, and l are equal to 1, 2,
. . . , d, and d is the spatial dimension. D pqr s denotes the elasticity tensor of the base material of LUCs. H (φn (x))
(i j)
is a Heaviside function used to indicate the solid( structural domain of the n-th LUC. The displacement field un
(i j)
)
corresponds to the locally varied strain field ε∗pq un which can be computed by solving the following equation
∫
ε pq − ε∗pq u(in j) D pqr s εr∗s v(kl) H (φn (x)) dΩn = 0, ∀v(kl) ∈ U (Ωn )
( 0(i j) ( )) ( )
n n (11)
Ωn

where vn(kl) represents the virtual displacement field of the nth LUC, and U (Ωn ) denotes the kinematically admissible
displacement field under the assumption of periodic boundary conditions.
7
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

3.2. Efficient property prediction of hierarchal multi-morphology LUCs by Kriging metamodels

Estimating the effective mechanical properties of hierarchal multi-morphology LUCs within HMMLSs by HM
will consume massive computation time and resources during the optimization process, and the Kriging metamodels
are adopted in this paper to achieve properties prediction.
To be more precise, in the proposed method, M Kriging metamodels need to be constructed for the M categories
of hierarchal LUCs (i.e. LM). For the mth category of LM, several hierarchal LUCs are selected as the sample
LUCs with their relative densities denote by ρ m . Then the mth Kriging metamodel assumes the equivalent elasticity
tensors of the mth category LM are the combination of a global response and local deviations:
( )T
DmH ρ m = f ρ m β + Z ρ m
( ) ( )
(12)
T
where ρ m denotes the prediction points,
( )and Dm ρ m is the observed response. f ρ m β provides the global
H
( ) ( )

( ) estimation for the input ρ m . f ρ m and β are the regression function and regression coefficient, respectively.
response
Z ρ m is a realization function of a stochastic process with mean zero and variance σz2 , and its nonzero covariance
can be computed by
Cov Z ρ mi , Z ρ m j = σz2 R ρ mi , ρ m j ; θ
[ ( ) ( )] ( )
(13)
where ρ mi and ρ m j are two sample points. R is the stochastic process correlation function concerning unknown
correlation parameter vector θ , and the Gaussian correlation function is widely adopted as R. And according to the
Gaussian process regression theory, the response of ρ m is subject to a normal distribution [63,64]:
( ( ))
D̂mH ρ m ∼ N D̂mH ρ m , σ̂D̂2 H ρ m
( ) ( )
(14)
m

Therefore, for n-th LUC belonging to the m-th category of LM, its elasticity tensor DiHjkl (ρmn ) can be approximated
by the constructed Kriging metamodel based on its relative density ρmn :
DiHjkl (ρmn ) ≈ D̂mH (ρmn ) (15)

4. Multiscale concurrent topology optimization model of HMMLSs

To obtain the optimized design of HMMLSs, determining the distribution of distinct LMs and the hierarchical
pattern of LUCs’ relative densities are the keys to be solved. Therefore, two kinds of design variables are used in
this paper. The first kind of design variable, ρmn
cat
is a discrete variable indicating the category of the LUC within
the HMMLS. The second kind of design variable, ρmn is a continuous variable representing the relative density of
the LUC. They are defined as follows:
ρmn ρmn = 0 or 1, (m = 1, 2, . . . , M; n = 1, 2, . . . , N )
cat
∑M cat
Discrete variable : = 0 or 1, m=1
(16)
Continuous variable : ρmin ≤ ρmn ≤ ρmax , (m = 1, 2, . . . , M; n = 1, 2, . . . , N )
where ρmin and ρmax define the lower and upper bounds of ρmn , and their values are determined by considering
the manufacturability of the mth category of hierarchical LUCs. If two categories of LMs exist in the HMMLS
(i.e. M = 2), then each phase LM within the HMMLS can be determined by:
∑2
m=1 ρmn = 0
cat
Void :
Lattice material 1 : ρ1n = 1, ρ2n
cat cat
=0 (17)
Lattice material 2 : ρ1n = 0, ρ2n = 1
cat cat

Then the elasticity tensor of the nth element within the HMMLS Dn ρmn , ρmn can be obtained by the proposed
( cat )
KUMMI model, which can be stated as:
⎧ ⎡ ⎤ ⎫
M ⎨ M ( ( )p) ⎬
) ∑ ( cat ) p ∏
Dn ρmn , ρmn = (1 − δ) ⎣ ρmn 1 − ρmcat1 n ⎦ + δ DiHjkl (ρ̃mn )
( cat
(18)
⎩ ⎭
m=1 m 1 =1(̸=m)

where (ρ̃mn ) is the elasticity tensor of nth LUC belonging to the mth category of LM, and its value is predicted
DiHjkl
by the constructed mth Kriging metamodel. ρ̃mn is the density filtered value of ρmn to ensure the existence of
solutions and avoid the formation of checkerboard patterns [65], and the density filtering scheme can be stated as
Hni = max (0, rmin − ∆ (n, i)) (19)
8
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

1 ∑
ρ̃mn = ∑ Hni ρmn (20)
f
i∈Nn
Hni f
i∈Nn

where Hni is the weight factor for element n in the density filtering scheme. ∆ (n, i) denotes the center-to-center
f
distance between element n and elements i. Nn is the set of elements i for which ∆ (n, i) is smaller than the filter
radius rmin , and the superscript f denotes the filter scheme. p denotes the penalization factor with the typical value
of 3 in the topology optimization problem. Considering the void regions within the HMMLSs, δ is an extremely
small positive parameter to prevent the stiffness matrix singularity, e.g. δ = 10−9 . Obviously, solving the values of
ρmn
cat
directly is inefficient as it’s a binary variable, and the threshold projection is adopted in this paper to relax the
discrete design variable ρmn cat
to a continuous design variable smn (0 ≤ smn ≤ 1), which can be expressed as [66]:
tanh (βη) + tanh (β (s̃mn − η))
ρmn
cat
= (21)
tanh (βη) + tanh (β (1 − η))
1 ∑
s̃mn =∑ Hni smn (22)
f Hni
i∈N n f
i∈Nn

where s̃mn is the density filtered value of the continuous design variable smn . β denotes the projection parameter
which can be used to control the curvature of the projection function, and the right-hand side of Eq. (21) is equivalent
to a step function when β → ∞. All the filtered design variable s̃mn above the threshold η are projected to 1 and the
values below to 0 after several iteration steps. Therefore the obtained ρmn
cat
could achieve binary values to determine
the LUC’s category. Note that 0 ≤ η ≤ 1 is required to execute the projection, good convergence can be obtained
when the value of η is close to 0.5, and almost discrete designs can be obtained by large β [67].
Based on the discussions above, the multiscale concurrent topology optimization model for HMMLS, which
takes into account the traditional minimum-compliance problem, can be formulated as follows:
Find : smn( , ρ mn )
Minimize : C smn , ρ mn = FT U
Subjectto : G smn ,(ρ mn = ) m=1 n=1 ρmn ρ̃mn vn − Vmax ≤ 0
( ) ∑M ∑N cat

F = K smn , ρ mn U (23)
K smn , ρ mn = n=1 Ωn B Dn ρmn (smn ) , ρmn BdΩn
( ) ∑N ∫ T
( cat )

0 ≤ smn ≤ 1, (m = 1, 2, . . . , M; n = 1, 2, . . . , N )
ρmin ≤ ρmn ≤ ρmax , (m = 1, 2, . . . , M; n = 1, 2, . . . , N ) .
where smn and ρ mn are two vectors of continuous variables smn and ρmn . C smn , ρ mn denotes the structural
( )
design
compliance of the HMMLS. F, U and K (smn , ρ mn ) represent the external load vector, displacement vector, and
( )

structural stiffness matrix, respectively. G smn , ρ mn stands for the global volume constraint function with vn is
the volume of a finite element. Vmax is the prescribed maximum volume of solid structures. B and Ωn are the
strain–displacement matrix and the domain of the nth finite element, respectively. What’s more, to pace the way
for gradient-based iteration for the relative density of the n-th LUC, the elasticity tensor of n-th LUC belonging to
the m-th category of LM can be interpolated by:
DiHjkl (ρ̃mn ) = ρ̃mn D̃0mn (ρ̃mn ) (24)
where D̃0mn (ρ̃mn ) is a temporal variable varied with the optimization process, and its application is for better
numerical implementation in designing HMMLS. Then according to Eq. (15), D̃0mn (ρ̃mn ) can be calculated by:

D̃0mn (ρ̃mn ) ≈ D̂mH (ρ̃mn ) /ρ̃mn (25)

One thing must be pointed out that D̃0mn (ρ̃mn ) does not represent the natural property of LUC, and its existence
aims to realize the “penalization” similar to the SIMP interpolation model. Before each iteration, D̃0mn (ρ̃mn ) will
be computed based on Eq. (25) element-wisely where D̂mH (ρ̃mn ) and ρ̃mn are obtained from the previous iteration.
9
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

5. Sensitivity analysis
To resolve the multivariable multiscale concurrent topology optimization of HMMLSs efficiently, the method
of moving asymptotes (MMA) [68] is employed in this paper to update the two kinds of design variables
simultaneously. Therefore, it’s vital to deduce the sensitivity of the objective and constraint functions with respect
to the design variables.
According to the chain rule, the first-order derivatives of the objective function C smn , ρ mn can be derived as
( )

∂C smn , ρ mn ∂C smn , ρ mn ∂ρmn ∂s̃mn

( ) ( )
cat
= (26)
∂smn ∂ρmn
cat ∂s̃mn ∂smn
∂C smn , ρ mn ∂C smn , ρ mn ∂ρ̃mn
( ) ( )
= (27)
∂ρmn ∂ρ̃mn ∂ρmn
The derivative of the objective function C smn , ρ mn with respect to ρmn
( ) cat
can be deduced as:
∂ Ωn B Dn ρmn (smn ) , ρmn BdΩn
∫ ( )
∂C smn , ρ mn T ∂K smn , ρ mn
T cat
( ) ( )
T
= −U U = −Un Un (28)
∂ρmn
cat ∂ρmn
cat ∂ρmn
cat

Here the middle term of the right-hand side is expanded by introducing Eqs. (15) and (18) as
∂ Ωn BT Dn ρmn (smn ) , ρmn BdΩn
∫ ( cat )
T ∂Dn ρmn (smn ) , ρmn
∫ ( cat )
= B BdΩn (29)
∂ρmn
cat
Ωn ∂ρmn
cat

∂Dn ρmn (smn ) , ρmn

( cat )
= (1 − δ)
⎡ ∂ρmncat
⎛
M ( ( )p)
( cat ) p−1 ∏
⎣ p ρmn ⎝ 1 − ρmcat1 n D̂mH (ρ̃mn ) (30)
m 1 =1(̸=m)
]
( )p ∏ ( ( )p) ))
ρmcat2 n 1 − ρmcat3 n D̂mH2 ρ̃m 2 n
∑M M (
− m 2 =1(̸=m) m 3 =1(̸=m,̸=m 2 )

Then based on Eq. (21) the middle term of the right-hand side of Eq. (26) can be calculated as
∂ρmn β 1 − (tanh (β (s̃mn − η)))2
( )
cat
= (31)
∂s̃mn tanh (βη) + tanh (β (1 − η))
and the last term of the right-hand side of Eq. (26) can be substituted by
∂s̃mn H jn
=∑ (32)
∂smn i∈N
f Hni
n
f
where Hni , H jn , Nn are the indexing matrices used in the density filtering scheme [65].
Similarly, Eq. (27) can be calculated with the help of Eq. (18) as follows:
(∫ )
∂C smn , ρ mn T ∂K smn , ρ mn T ∂Dn ρmn (smn ) , ρmn
( ) ( ) ( cat )
T
= −U U = −Un B BdΩn Un (33)
∂ρ̃mn ∂ρ̃mn Ωn ∂ρ̃mn
⎧ ⎡ ⎤ ⎫
∂Dn ρmn
(
cat
(smn ) , ρmn
) ⎨ M ( ( )p) ⎬ ∂D H (ρ̃mn )
cat p i jkl
∏
= (1 − δ) ⎣ ρmn 1 − ρmcat1 n ⎦+δ
( )
(34)
∂ρ̃mn ⎩
m 1 =1(̸=m)
⎭ ∂ρ̃mn
∂ρ̃mn H jn
=∑ (35)
∂ρmn i∈N
f Hni
n

It’s noticeable that based on the constructed Kriging metamodels in Eq. (15), the rightmost term of Eq. (34) can
be rewritten as
∂DiHjkl (ρ̃mn ) ∂ D̂mH (ρ̃mn )
≈ (36)
∂ρ̃mn ∂ρ̃mn
10
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Note that the values of D̂mH (ρ̃mn ) are directly predicted by the Kriging metamodels based on the input vectors
ρ̃mn . However, obtaining the derivatives of the Kriging metamodels with respect to ρ̃mn is computationally inefficient
since the stochastic process is involved. In this work, the material interpolation scheme in Eq. (24) is adopted as
an alternative solution, and the solution of Eq. (36) can be obtained as
∂DiHjkl (ρ̃mn )
= D̃0mn (ρ̃mn ) (37)
∂ρ̃mn
What’s more, the sensitivities of constraint function G smn , ρ mn with respect to smn and ρmn are
( )

∂G smn , ρ mn ∂ρ cat ∂s̃mn

( )
= ρ̃mn mn vn (38)
∂smn ∂s̃mn ∂smn
∂G smn , ρ mn cat ∂ρ̃mn
( )
= ρmn vn (39)
∂ρmn ∂ρmn

6. Numerical examples
In this section, three numerical examples are provided to validate the effectiveness of the proposed multiscale
concurrent topology optimization method for designing HMMLSs. Especially, for the given design domain and
boundary conditions, different numbers of LMs are employed to validate the advantages of introducing multi-
morphology LUCs in designing hierarchical lattice structures. It’s assumed that the base material used to construct
LUCs is linearly elastic, and its Young’s modulus E = 2750 MPa, and Poisson’s ratio µ = 0.38. The value of the
threshold η is equal to 0.5, the β is set to 1 at the beginning and will be continuously increased after each iteration
step to avoid local optima and accelerate convergence.

6.1. 2D cantilever beam design

The first example investigates the optimization problem of a 2D cantilever beam structure with length L =
400 mm and height W = 250 mm. The beam is clamped at its left end and a vertical load F = 100 N is located at
the center of its right-hand edge as shown in Fig. 5. The entire design domain is discretized into 80 × 50 four-node
quadrilateral plane stress elements with the elemental size of 5 mm. The optimization objective is minimizing
the structural compliance of HMMLS under the global volume constraint of 35%. The optimization process will
terminate when the maximum element-wise change of design variables in two successive iteration steps is below
0.001, or the maximum 200 iteration steps is reached. What’s more, as shown in Fig. 6, fifteen categories of 2D
LUCs are employed in this subsection to construct candidate LMs for the topology optimization of HMMLSs, and
the values of LUCs’ relative densities are restricted to [ρmin , 1].
To investigate the benefits of introducing multiple LMs in designing hierarchical lattice structures, different
combinations of candidate LMs are employed in designing HMMLSs for 2D cantilever beam. The optimized design
of HMMLSs is presented in Fig. 7, where different colors denote disparate LMs. Firstly, as shown in Fig. 7(a), (b),
(d), (e), and (f), it’s evident that LMs with high x-direction tensile modulus (e.g. LM. 2) are mainly distributed in
the upper and lower layers of the HMMLSs, LM with good y-direction tensile modulus (e.g. LM. 9) is located near

Fig. 5. Sketch of the 2D cantilever beam with its load and boundary conditions.

11
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Fig. 6. Different categories of 2D candidate LUCs used to construct LMs.

the load region. As for the other LMs with preferable tensile modulus along the 45-degree axis and 135-degree axis
(e.g. LM. 7), they predominately occupy the interior regions of the beam. To further investigate this phenomenon,
the stress distribution of the cantilever beam with solid material is simulated using the finite element commercial
software COMSOL Multiphysics 5.5, and the results are presented in Fig. 8. The magnitude of von Mises stress for
the 2D cantilever beam is handled to better present the spatially varying stress, and the directions of principal stress
are denoted by red and blue cones together with their size determined by the magnitude of the principal stress.
It’s fascinating to discover the distribution pattern of LMs matches perfectly with the orientations of the major
principal stress within the cantilever beam, and LUCs with high relative densities are mainly located at the regions
with higher stress. It needs to be emphasized that no stress analysis is conducted in the compliance minimizing
procedure, and the stress distribution presented here is solely to demonstrate the proposed design method is capable
of achieving rational distribution of candidate LMs. Specially speaking, LUCs with different relative densities and
categories characteristics are guided by the topology optimization model to adaptively fill the cantilever beam for
better adaptation to various deformation conditions.
In addition, the distribution regions of some categories of LMs (e.g. LM. 1, LM. 7) are symmetric about the x-axis
of the cantilever. This is partly due to the boundary and load conditions of the design domain being symmetrical
so that the regions where the axial loads or tangential loads play the dominant role within the beam are symmetric.
Besides, the symmetry along the x-axis and y-axis of the LUCs also make contributions to this phenomenon.
Meanwhile, some LMs with a directional strengthened structure (e.g. LM. 3 and LM. 5) within the cantilever
beam are asymmetric. This can be attributed to the main load path in these regions preferring LMs with specifically
directional characteristics of elasticity tensor, which is observed in all design cases of Fig. 7. What’s more, the LUCs
with high relative densities are mainly located on the main load path of the cantilever, as shown in the left column
in Fig. 7. It’s easy to conclude that the combination modes of candidate LMs will also influence the distribution of
relative densities for LUCs. It is important to note that the optimization of LUCs’ relative densities is completed
simultaneously with the evolutionary process of the distribution of LMs, thereby achieving multiscale concurrent
topology optimization of HMMLSs.
The introduction of more candidate LMs can further enhance the structural stiffness (i.e. decrease the compli-
ance). For instance, the compliance of the cantilever beam in Fig. 7(a) is 232.55, whereas the value for the beam in
Fig. 7(c) is 4.15% lower than in Fig. 7(a). The reason for the performance improvement is that by introducing more
candidate LMs, better compatibility between the main load-bearing directions of the cantilever beam and hierarchical
12
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Fig. 7. The optimized design of 2D cantilever beams by using different combinations of candidate LMs under the global volume constraint of
35%. (In each design case, the left picture shows the optimized distribution of relative densities for LUCs where red and blue colors denote
LUCs with high and low relative density respectively, and the right picture shows the optimized HMMLSs together with the percentages
of different LMs within the structure). (For interpretation of the references to color in this figure legend, the reader is referred to the web
version of this article.)

multi-morphology LUCs is obtained. However, it must be pointed out that introducing more candidate LMs does
not always guarantee performance improvements since the introduced LMs may not appear in the final design,
which can be partly verified by comparing the optimized design of Fig. 7(a) and (b). The difference in compliances
for the cantilever beam in Fig. 7(a) and (b) is only 0.11%, despite the use of different candidate LMs. Further
13
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Fig. 8. Sketch of the distribution of stress for the 2D cantilever beam. (The directions of principal stress are denoted by red and blue cones).
(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

conclusions can be drawn that a sufficient quantity of candidate LMs with special configurations are worth trying
when designing HMMLSs, as this will provide more design space to pursue excellent mechanical performance.
For simplicity, the iterative process of designing the cantilever beam in Fig. 7(f) is illustrated in Fig. 9 and
Fig. 10, where the evolution history of the distribution regions and percentages of different candidate LMs are both
presented. Violent fluctuations of structural compliance can be observed in the first several iterative steps due to
the penalization factor p in the KUMMI model. Specifically, the category of LUC within a finite element is not
determined in the earlier stage, and this led to a serious penalty on the elasticity tensor of LUCs, huge compliance
is presented as a result. After a quick decrease of the structural compliance in the first 10 iterations, the evolution
becomes moderate until it converges to the final value C = 215.19 when the maximum iteration step is reached.
During the iteration steps, some LMs’ volume fractions approach zero, while others increase or decrease initially
and then remain stable. For better visualization of the optimized design for the 2D cantilever beam, Fig. 11 presents
the final structure of Fig. 8(f), where only 6 of the 15 candidate LMs reserved. The percentages, relative densities,
and distribution regions of all candidate LMs are redistributed more efficiently for better mechanical performance.
The CPU times of the proposed multiscale concurrent topology optimization of HMMLSs are shown in Fig. 12.
Especially, T D H denotes the CPU time to obtain the effective elasticity tensors of hierarchical multi-morphology
LUCs. T F E A is the CPU time for the finite element analysis (FEA) of HMMLSs. Tsen and Tita represent the
CPU times for computing the sensitivities of the objective function with respect to the design variables and
updating the design variables by the MMA optimizer, respectively. It can be seen that TF E A and Tita are far less
than T D H and Tsen , which means that the main computational burden involves obtaining the effective elasticity
tensors of hierarchal multi-morphology LUCs. However, the total time in an optimization iteration step is Ttotal =
T D H + T F E A + Tsen + Tita = 4.7694 s, an acceptable computational burden for realizing the multiscale topology
optimization of HMMLSs. By contrast, T D H in an iteration step will be 2046.73s without the proposed KUMMI
model, unacceptable in designing HMMLSs.

6.2. 3D Messerschmidt-Bölkow-Blohm (MBB) beam design

To investigate the feasibility of the proposed multiscale concurrent topology optimization method for designing
HMMLSs in 3D cases, different combinations of 3D candidate LMs are employed in this paper for design of 3D
MBB beams. As shown in Fig. 13, the dimensions of the MBB beam are length L = 300 mm, width W = 70 mm,
and height H = 60 mm. A vertically loaded distributed force F = 3000 N is located at the center of the top face
of the MBB beam. Besides, the bottom-left edge of the MBB beam is on fixed supports while the bottom-right
edge is supported on rollers. A mesh with 60 × 14 × 12 eight-node hexahedral elements is used to discretize the
design domain. 20 categories of 3D candidate LMs are adopted in this section as presented in Fig. 14, and the upper
14
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Fig. 9. Iteration histories of objective and constraint function for 2D cantilever beam with 15 categories of candidate LMs under the global
volume constraint of 35%.

Fig. 10. Topology optimization process for the 2D cantilever beam with 15 categories of candidate LMs under the global volume constraint
of 35%.

and lower limits of relative densities for LUCs are set to guarantee that there is no fine feature, fault structure, or
enclosed voids within the LUC. And the maximum number of design iterations for 3D cases is 150 unless the
greatest element-wise change falls below the threshold of 0.001.
In order to evaluate the effectiveness of the optimized layout of LMs within the HMMLSs, the stress distribution
of the MBB beam with solid material is shown in Fig. 15. By observing the distribution pattern of principal stress
orientations, it’s not hard to get the speculation that the regions near the top face and bottom face should be filled
by LM with high anisotropy. Moreover, LUCs with different direction reinforcements are required to fill the left
15
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Fig. 11. Optimized design of HMMLS for 2D cantilever beam with 15 categories of candidate LMs.

and right domains of the MBB beam, since the directions of principal stress in these domains are different, although
the magnitudes of stress are similar.
The optimized results of HMMLSs are given in Fig. 16, where various combinations of 3D candidate LMs are
employed in different design cases. A global volume constraint of 30% is applied in all design cases for better
comparison. Not surprisingly, LMs with high anisotropy (e.g. LM.10, LM.11) are mainly located on the left and
right ends of the MBB beam, as shown in Fig. 16(d), (f), and (g). As for the middle region of the beam, nearly
isotropic LMs (e.g. LM.16, LM.17) prefer to be settled here if they are selected as candidate LMs. Additionally,
as presented in Fig. 16(d), (f), and (g), although the distribution region of LMs is not symmetrical within the
beam (e.g. LM.8, LM.9), the whole multiscale structure is almost symmetric along the x-z plane and y-z plane to
accommodate the symmetric boundary conditions. Furthermore, when LMs with similar mechanical properties are
employed as the candidate materials, LMs with better stiffness-to-density will be adopted in the optimized HMMLSs
and the weaker is discarded. For example, LM.5 within the center region of the MBB beam will be substituted by
16
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Fig. 12. CPU times of an optimization iteration step for 2D cantilever beam with 15 categories of candidate LMs.

Fig. 13. Sketch of 3D MBB beam with its load and boundary conditions.

LM.16 when LM.16 is appended as candidate LMs, and LM.17 will replace LM.16 because the stiffness-to-density
ratio of LM.17 is superior to LM.16, as shown in Fig. 16(a), (e), and (f).
Despite having the same boundary conditions and constraint function, the optimized distribution of relative
densities for LUCs varies in different design cases, as illustrated in the left pictures of Fig. 16. In other words,
the distribution pattern of LUCs’ relative densities will be affected by the categories of candidate LMs. Besides,
multi-morphology LUCs with identical densities are allocated in different regions to better resist various deformation
behavior of HMMLSs. For example, LUCs with relative densities of ρmax = 0.7 are located in the central region
of the top face of the MBB beam, the LM. 8 and LM.9 here are mutually matched to better resist compressive
deformation along different directions near the loaded area.
As shown in Fig. 16, the structural compliance of the HMMLS in design case (g) is the smallest under the
global volume constraint of 30%, where six categories of LMs are retained in the final design. The best structural
compliance is 6064.82 in case (g), 8.53% lower than the worst 6630.19 in case (a). What’s more, compared with
HMMLS purely filled with truss-based LMs in case (d), the introduction of TPMS-based LMs could ulteriorly
enhance the structural performance, as shown in case (f) and case (g). Therefore, it is recommended to introduce
truss-based and TPMS-based candidate LMs while designing HMMLSs. The corresponding optimization process
of case (g) is presented in Fig. 17, together with the evolution of distribution regions for different candidate LMs
presented in Fig. 18. The KUMMI model proposed in this study has been found to catch suitable LMs for designing
HMMLSs, leading to fast and stable convergence procedures of objective and constraint functions. As the iteration
proceeds, clear distribution regions of different categories of 3D LMs are obtained, accompanied by a decrease in
17
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Fig. 14. Different categories of 3D LUCs used to construct LMs.

Fig. 15. Sketch of the distribution of stress for the 3D MBB beam (the directions of principal stress are denoted by red, blue, and green
cones). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

18
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Fig. 16. The optimized design of 3D MBB beams by using different combinations of candidate LMs under the global volume constraint of
30% (In each design case, the left picture shows the optimized distribution of relative densities for LUCs where red and blue colors denote
LUCs with high and low relative density respectively, and the right picture shows the optimized HMMLSs together with the percentages
of different LMs within the structure). (For interpretation of the references to color in this figure legend, the reader is referred to the web
version of this article.)

19
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Fig. 17. Iteration histories of objective and constraint function for 3D MBB beam with 20 categories of candidate LMs under the global
volume constraint of 30%.

Fig. 18. Topology optimization process for the 3D MBB beam with 20 categories of candidate LMs under the global volume constraint of
30%.

structural compliance. Fig. 19 presents the final structures after the iteration loop presented in Fig. 17 in detail.
The LUCs with high relative densities are distributed along the main load paths for bridging the loaded and
constraint boundaries. Moreover, thanks to the proposed SF-based hybrid transition strategy, the hierarchical LUCs
of different LMs are smoothly connected, which is beneficial to alleviate stress concentration phenomena resulting
from the incompatible neighboring LUCs. In other words, the structural connectivity between TPMS-based LMs
and truss-based LMs is improved.

6.3. 3D supported structure design

In this example, the proposed multiscale concurrent topology optimization method is employed to design a 3D
supported structure. The design domain with length L = 120 mm, width W = 120 mm, and height H = 70 mm
is illustrated in Fig. 20. Its four corners at the bottom surface are fixed, and five vertically loaded concentrated
20
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Fig. 19. Optimized design of HMMLS for 3D MBB beam with 20 categories of candidate LMs.

Fig. 20. Sketch of 3D supported structure with its load and boundary conditions.

force F = 600 N are applied on the top face of the MBB beam. A mesh with 24 × 24 × 14 eight-node hexahedral
elements is used to discretize the design domain. All 20 categories of 3D candidate LMs are adapted for the topology
optimization of the supported structure. The optimization objective is to minimize the overall structural compliance
under the global volume constraint of 30%.
Fig. 21 illustrates the stress distributions of the 3D supported structure, the directions of principal stress are
spatially varying which confirmed that adaptable hierarchical multi-morphology LUCs are needed to fill the design
domain. The optimization histories of the 3D supported structure is presented in Figs. 22 and 23, and the optimized
structural compliance is 2717.27. The iterative curves of the objective function and volume fractions of different
candidate LMs demonstrate the excellent capacity of the proposed topology optimization scheme in concurrently
21
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Fig. 21. Sketch of the distribution of stress for the 3D supported structure (the directions of principal stress are denoted by red, blue, and
green cones). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 22. Iteration histories of objective and constraint function for 3D supported structure with 20 categories of candidate LMs under the
global volume constraint of 30%.

determining the suitable combinations of LMs, the relative densities of LUCs, and the percentages of different LMs
within the optimized HMMLSs. The optimized design for 3D supported structure is presented in Fig. 24, with six
categories of candidate LMs reserved, and the load-bearing direction of LUCs exhibits excellent coordination to the
supported structure’s load and boundary conditions. The interior region of the design domain is mainly occupied by
near isotropic LMs (i.e. LM.5 and LM.17). Meanwhile, LMs with strong directional characteristics dominate the
outer border, and thicker rods with different orientations connect the fixed point and force-bearing point to enhance
force transfer paths. This also reveals the effectiveness and adaptive capacity of the proposed method for design of
3D HMMLSs. The above observations effectively clarify the necessity of introducing hierarchical multi-morphology
LUCs in the multiscale design of HMMLSs.
22
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Fig. 23. Topology optimization process for the 3D supported structure with 20 categories of candidate LMs under the global volume constraint
of 30%.

7. Conclusions
This paper introduces a novel multiscale concurrent topology optimization method for designing hierarchical
multi-morphology lattice structures, where multiple hierarchical multi-morphology LUCs are employed in a partic-
ular pattern to fill the design domain. The proposed method utilizes level set functions to depict the configuration of
multi-morphology LUCs so that configuration mapping can be employed to generate hierarchical multi-morphology
LUCs. Additionally, the method uses the SF-based hybrid transition strategy to ensure the smooth transition between
multi-morphology LUCs. Kriging metamodels are constructed based on the homogenized properties of sample LUCs
to realize efficient prediction of the elasticity tensor for hierarchical multi-morphology LUCs within the HMMLSs.
Based on the constructed Kriging metamodels, KUMMI model is put forward to characterize the interpolated
material properties during the topology optimization procedure. The optimization model of HMMLSs for compliance
minimization is constructed using two types of design variables, discrete variables indicating categories of LUCs
and continuous variables determining relative densities of LUCs. Sensitivity analysis is carried out to evolve the two
kinds of design variables concurrently, and the optimized design of HMMLSs with smooth connected hierarchical
multi-morphology LUCs infill can be obtained.
Three numerical examples are presented in this paper to demonstrate the effectiveness of the proposed multiscale
concurrent topology optimization method in enhancing lattice structural performance. Additionally, the proposed
SF-based hybrid transition strategy realizes the smooth connectivity of multi-morphology LUCs, which makes
it possible to combine TPMS-based LUCs and truss-based LUCs in optimizing HMMLSs. It also means that
more categories of LUCs with distinct morphology can be employed for the multiscale topology optimization
of HMMLSs to seek better structural properties. By comparing the static structural analysis and optimization
results, the optimized rational distribution patterns of hierarchical multi-morphology LUCs within HMMLSs clearly
demonstrate the effectiveness of the proposed multiscale concurrent topology optimization method. Apart from the
compliance minimization problem here, future work will focus on the implementation of the proposed method for
designing HMMLSs considering other physical fields, such as dynamics, thermology, and acoustics. Furthermore,
the categories of LUCs in the proposed method are not fixed, other kinds of LUCs can also be designed and then
employed to generate candidate LMs, realizing more satisfactory performance improvement in designing HMMLSs.
23
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

Fig. 24. Optimized design of HMMLS for 3D supported structure with 20 categories of candidate LMs.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could
have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.

Acknowledgments
This research was supported by the National Key Research and Development Program of China [grant number
2020YFB1708300] and the XPLORER PRIZE.

References
[1] H.N.G. Wadley, Cellular metals manufacturing, Adv. Eng. Mater. 4 (10) (2002) 726–733, http://dx.doi.org/10.1002/1527-2648(20021014)
4:10<726::AID-ADEM726>3.0.CO;2-Y.
[2] D.-J. Yoo, Heterogeneous porous scaffold design for tissue engineering using triply periodic minimal surfaces, Int. J. Precis. Eng.
Manuf. 13 (4) (2012) 527–537, http://dx.doi.org/10.1007/s12541-012-0068-5.
[3] A.O. Aremu, et al., A voxel-based method of constructing and skinning conformal and functionally graded lattice structures suitable
for additive manufacturing, Addit. Manuf. 13 (2017) 1–13, http://dx.doi.org/10.1016/j.addma.2016.10.006.
24
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

[4] A. Clausen, N. Aage, O. Sigmund, Exploiting additive manufacturing infill in topology optimization for improved buckling load,
Engineering 2 (2) (2016) 250–257, http://dx.doi.org/10.1016/J.ENG.2016.02.006.
[5] S. Xu, J. Shen, S. Zhou, X. Huang, Y.M. Xie, Design of lattice structures with controlled anisotropy, Mater. Des. 93 (2016) 443–447,
http://dx.doi.org/10.1016/j.matdes.2016.01.007.
[6] R. Alberdi, et al., Multi-morphology lattices lead to improved plastic energy absorption, Mater. Des. 194 (2020) 108883, http:
//dx.doi.org/10.1016/j.matdes.2020.108883.
[7] L. Cheng, P. Zhang, E. Biyikli, J. Bai, J. Robbins, A. To, Efficient design optimization of variable-density cellular structures for
additive manufacturing: Theory and experimental validation, Rapid Prototyp. J. 23 (4) (2017) 660–677, http://dx.doi.org/10.1108/RPJ-
04-2016-0069.
[8] B. Wang, G.D. Cheng, L. Jiang, Design of multi-tubular heat exchangers for optimum efficiency of heat dissipation, Eng. Optim. 40
(8) (2008) 767–788, http://dx.doi.org/10.1080/03052150802054027.
[9] H.A. Eschenauer, N. Olhoff, Topology optimization of continuum structures: A review*, Appl. Mech. Rev. 54 (4) (2001) 331–390,
http://dx.doi.org/10.1115/1.1388075.
[10] M.P. Bendsøe, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Comput. Methods Appl.
Mech. Engrg. 71 (2) (1988) 197–224, http://dx.doi.org/10.1016/0045-7825(88)90086-2.
[11] M.P. Bendsøe, O. Sigmund, Material interpolation schemes in topology optimization, Arch. Appl. Mech. 69 (9) (1999) 635–654,
http://dx.doi.org/10.1007/s004190050248.
[12] O. Sigmund, A 99 line topology optimization code written in matlab, Struct. Multidiscip. Optim. 21 (2) (2001) 120–127, http:
//dx.doi.org/10.1007/s001580050176.
[13] M.P. Bendsoe, O. Sigmund, Topology Optimization: Theory, Method and Applications, Springer, 2003.
[14] O. Sigmund, K. Maute, Topology optimization approaches: A comparative review, Struct. Multidiscip. Optim. 48 (6) (2013) 1031–1055,
http://dx.doi.org/10.1007/s00158-013-0978-6.
[15] H.A. Eschenauer, N. Olhoff, Topology optimization of continuum structures: A review*, Appl. Mech. Rev. 54 (4) (2001) 331–390,
http://dx.doi.org/10.1115/1.1388075.
[16] Y. Zhang, L. Zhang, Z. Ding, L. Gao, M. Xiao, W.-H. Liao, A multiscale topological design method of geometrically asymmetric
porous sandwich structures for minimizing dynamic compliance, Mater. Des. 214 (2022) 110404, http://dx.doi.org/10.1016/j.matdes.
2022.110404.
[17] H. Rodrigues, J.M. Guedes, M.P. Bendsoe, Hierarchical optimization of material and structure, Struct. Multidiscip. Optim. 24 (1) (2002)
1–10, http://dx.doi.org/10.1007/s00158-002-0209-z.
[18] L. Liu, J. Yan, G. Cheng, Optimum structure with homogeneous optimum truss-like material, Comput. Struct. 86 (13) (2008) 1417–1425,
http://dx.doi.org/10.1016/j.compstruc.2007.04.030.
[19] X. Yan, X. Huang, Y. Zha, Y.M. Xie, Concurrent topology optimization of structures and their composite microstructures, Comput.
Struct. 133 (2014) 103–110, http://dx.doi.org/10.1016/j.compstruc.2013.12.001.
[20] R. Sivapuram, P.D. Dunning, H.A. Kim, Simultaneous material and structural optimization by multiscale topology optimization, Struct.
Multidiscip. Optim. 54 (5) (2016) 1267–1281, http://dx.doi.org/10.1007/s00158-016-1519-x.
[21] L. Xu, G. Cheng, Two-scale concurrent topology optimization with multiple micro materials based on principal stress orientation,
Struct. Multidiscip. Optim. 57 (5) (2018) 2093–2107, http://dx.doi.org/10.1007/s00158-018-1916-4.
[22] Y. Zhang, M. Xiao, X. Zhang, L. Gao, Topological design of sandwich structures with graded cellular cores by multiscale optimization,
Comput. Methods Appl. Mech. Engrg. (2019) 112749, http://dx.doi.org/10.1016/j.cma.2019.112749.
[23] L. Cheng, J. Liu, A.C. To, Concurrent lattice infill with feature evolution optimization for additive manufactured heat conduction
design, Struct. Multidiscip. Optim. 58 (2) (2018) 511–535, http://dx.doi.org/10.1007/s00158-018-1905-7.
[24] D. Li, W. Liao, N. Dai, G. Dong, Y. Tang, Y.M. Xie, Optimal design and modeling of gyroid-based functionally graded cellular
structures for additive manufacturing, Comput.-Aided Des. 104 (2018) 87–99, http://dx.doi.org/10.1016/j.cad.2018.06.003.
[25] D. Li, N. Dai, Y. Tang, G. Dong, Y.F. Zhao, Design and optimization of graded cellular structures with triply periodic level surface-based
topological shapes, J. Mech. Des. 141 (7) (2019) 071402, http://dx.doi.org/10.1115/1.4042617.
[26] Y. Wang, S. Arabnejad, M. Tanzer, D. Pasini, Hip implant design with three-dimensional porous architecture of optimized graded
density, J. Mech. Des. 140 (11) (2018) 111406, http://dx.doi.org/10.1115/1.4041208.
[27] T. Wu, S. Li, An efficient multiscale optimization method for conformal lattice materials, Struct. Multidiscip. Optim. 63 (3) (2021)
1063–1083, http://dx.doi.org/10.1007/s00158-020-02739-5.
[28] R. Murphy, C. Imediegwu, R. Hewson, M. Santer, Multiscale structural optimization with concurrent coupling between scales, Struct.
Multidiscip. Optim. 63 (4) (2021) 1721–1741, http://dx.doi.org/10.1007/s00158-020-02773-3.
[29] N.M. Ngoc, V.-N. Hoang, D. Lee, Concurrent topology optimization of coated structure for non-homogeneous materials under buckling
criteria, Eng. Comput. 38 (6) (2022) 5635–5656, http://dx.doi.org/10.1007/s00366-022-01718-2.
[30] Y. Lu, L. Tong, Concurrent multiscale topology optimization of metamaterials for mechanical cloak, Comput. Methods Appl. Mech.
Engrg. 409 (2023) 115966, http://dx.doi.org/10.1016/j.cma.2023.115966.
[31] T.T. Banh, N.G. Luu, D. Lee, A non-homogeneous multi-material topology optimization approach for functionally graded structures
with cracks, Compos. Struct. (2021) 273, http://dx.doi.org/10.1016/j.compstruct.2021.114230.
[32] T.T. Banh, Q.X. Lieu, J. Lee, J. Kang, D. Lee, A robust dynamic unified multi-material topology optimization method for functionally
graded structures, Struct. Multidiscip. Optim. 66 (4) (2023) 75, http://dx.doi.org/10.1007/s00158-023-03501-3.
[33] S. Rastegarzadeh, J. Wang, J. Huang, Multi-scale topology optimization with neural network-assisted optimizer, presented at the ASME
2022, in: International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, American
Society of Mechanical Engineers Digital Collection, 2022, http://dx.doi.org/10.1115/DETC2022-89538.
25
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

[34] Y. Wang, M. Xiao, Z. Xia, P. Li, L. Gao, From computer-aided design (CAD) toward human-aided design (HAD): An isogeometric
topology optimization approach, Engineering (2022) http://dx.doi.org/10.1016/j.eng.2022.07.013.
[35] C. Imediegwu, R. Murphy, R. Hewson, M. Santer, Multiscale structural optimization towards three-dimensional printable structures,
Struct. Multidiscip. Optim. 60 (2) (2019) 513–525, http://dx.doi.org/10.1007/s00158-019-02220-y.
[36] C. Wang, J.H. Zhu, W.H. Zhang, S.Y. Li, J. Kong, Concurrent topology optimization design of structures and non-uniform parameterized
lattice microstructures, Struct. Multidiscip. Optim. 58 (1) (2018) 35–50, http://dx.doi.org/10.1007/s00158-018-2009-0.
[37] C. Wang, X. Gu, J. Zhu, H. Zhou, S. Li, W. Zhang, Concurrent design of hierarchical structures with three-dimensional parameterized
lattice microstructures for additive manufacturing, Struct. Multidiscip. Optim. 61 (3) (2020) 869–894, http://dx.doi.org/10.1007/s00158-
019-02408-2.
[38] D. Kang, S. Park, Y. Son, S. Yeon, S.H. Kim, I. Kim, Multi-lattice inner structures for high-strength and light-weight in metal selective
laser melting process, Mater. Des. 175 (2019) 107786, http://dx.doi.org/10.1016/j.matdes.2019.107786.
[39] F. Ren, et al., Transition boundaries and stiffness optimal design for multi-TPMS lattices, Mater. Des. 210 (2021) 110062, http:
//dx.doi.org/10.1016/j.matdes.2021.110062.
[40] Z. Liu, L. Xia, Q. Xia, T. Shi, Data-driven design approach to hierarchical hybrid structures with multiple lattice configurations, Struct.
Multidiscip. Optim. 61 (6) (2020) 2227–2235, http://dx.doi.org/10.1007/s00158-020-02497-4.
[41] L. Wang, S. Tao, P. Zhu, W. Chen, Data-driven topology optimization with multiclass microstructures using latent variable Gaussian
process, J. Mech. Des. 143 (3) (2020) http://dx.doi.org/10.1115/1.4048628.
[42] L. Wang, A. van Beek, D. Da, Y.-C. Chan, P. Zhu, W. Chen, Data-driven multiscale design of cellular composites with multiclass
microstructures for natural frequency maximization, Compos. Struct. 280 (2022) 114949, http://dx.doi.org/10.1016/j.compstruct.2021.
114949.
[43] J. Stegmann, E. Lund, Discrete material optimization of general composite shell structures, Internat. J. Numer. Methods Engrg. 62 (14)
(2005) 2009–2027, http://dx.doi.org/10.1002/nme.1259.
[44] T. Gao, W. Zhang, A mass constraint formulation for structural topology optimization with multiphase materials, Internat. J. Numer.
Methods Engrg. 88 (8) (2011) 774–796, http://dx.doi.org/10.1002/nme.3197.
[45] M.P. Bendsøe, O. Sigmund, Material interpolation schemes in topology optimization, Arch. Appl. Mech. Ing. Arch. 69 (9–10) (1999)
http://dx.doi.org/10.1007/s004190050248.
[46] M.Y. Wang, X. Wang, ‘Color’ level sets: A multi-phase method for structural topology optimization with multiple materials, Comput.
Methods Appl. Mech. Engrg. 193 (6–8) (2004) http://dx.doi.org/10.1016/j.cma.2003.10.008.
[47] X. Huang, Y.M. Xie, Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials, Comput.
Mech. 43 (3) (2009) 393–401, http://dx.doi.org/10.1007/s00466-008-0312-0.
[48] S. Zhou, M.Y. Wang, Multimaterial structural topology optimization with a generalized Cahn–Hilliard model of multiphase transition,
Struct. Multidiscip. Optim. 33 (2) (2006) 89–111, http://dx.doi.org/10.1007/s00158-006-0035-9.
[49] R. Tavakoli, S.M. Mohseni, Alternating active-phase algorithm for multimaterial topology optimization problems: A 115-line MATLAB
implementation, Struct. Multidiscip. Optim. 49 (4) (2014) 621–642, http://dx.doi.org/10.1007/s00158-013-0999-1.
[50] R. Tavakoli, Multimaterial topology optimization by volume constrained Allen–Cahn system and regularized projected steepest descent
method, Comput. Methods Appl. Mech. Engrg. 276 (2014) 534–565, http://dx.doi.org/10.1016/j.cma.2014.04.005.
[51] W. Sha, M. Xiao, L. Gao, Y. Zhang, A new level set based multi-material topology optimization method using alternating active-phase
algorithm, Comput. Methods Appl. Mech. Engrg. 377 (2021) 113674, http://dx.doi.org/10.1016/j.cma.2021.113674.
[52] Q.H. Doan, D. Lee, Optimum topology design of multi-material structures with non-spurious buckling constraints, Adv. Eng. Softw.
114 (2017) 110–120, http://dx.doi.org/10.1016/j.advengsoft.2017.06.002.
[53] N. Yang, Z. Quan, D. Zhang, Y. Tian, Multi-morphology transition hybridization CAD design of minimal surface porous structures
for use in tissue engineering, Comput.-Aided Des. 56 (2014) 11–21, http://dx.doi.org/10.1016/j.cad.2014.06.006.
[54] D.-J. Yoo, K.-H. Kim, An advanced multi-morphology porous scaffold design method using volumetric distance field and beta growth
function, Int. J. Precis. Eng. Manuf. 16 (9) (2015) 2021–2032, http://dx.doi.org/10.1007/s12541-015-0263-2.
[55] O. Al-Ketan, R.K. Abu Al-Rub, Multifunctional mechanical metamaterials based on triply periodic minimal surface lattices, Adv. Eng.
Mater. 21 (10) (2019) 1900524, http://dx.doi.org/10.1002/adem.201900524.
[56] M. Xiao, X. Liu, Y. Zhang, L. Gao, J. Gao, S. Chu, Design of graded lattice sandwich structures by multiscale topology optimization,
Comput. Methods Appl. Mech. Engrg. 384 (2021) 113949, http://dx.doi.org/10.1016/j.cma.2021.113949.
[57] X. Liu, L. Gao, M. Xiao, Y. Zhang, Kriging-assisted design of functionally graded cellular structures with smoothly-varying lattice
unit cells, Comput. Methods Appl. Mech. Engrg. 390 (2022) 114466, http://dx.doi.org/10.1016/j.cma.2021.114466.
[58] M. Xiao, W. Sha, Y. Zhang, X. Liu, P. Li, L. Gao, CMTO: Configurable-design-element multiscale topology optimization, Addit.
Manuf. 69 (2023) 103545, http://dx.doi.org/10.1016/j.addma.2023.103545.
[59] F. Tamburrino, S. Graziosi, M. Bordegoni, The design process of additively manufactured mesoscale lattice structures: A review, J.
Comput. Inf. Sci. Eng. 18 (4) (2018) 040801, http://dx.doi.org/10.1115/1.4040131.
[60] W. Chen, S. Watts, J.A. Jackson, W.L. Smith, D.A. Tortorelli, C.M. Spadaccini, Stiff isotropic lattices beyond the Maxwell criterion,
Sci. Adv. 5 (9) (2019) eaaw1937, http://dx.doi.org/10.1126/sciadv.aaw1937.
[61] E. Andreassen, C.S. Andreasen, How to determine composite material properties using numerical homogenization, Comput. Mater. Sci.
83 (2014) 488–495, http://dx.doi.org/10.1016/j.commatsci.2013.09.006.
[62] L. Xia, P. Breitkopf, Design of materials using topology optimization and energy-based homogenization approach in matlab, Struct.
Multidiscip. Optim. 52 (6) (2015) 1229–1241, http://dx.doi.org/10.1007/s00158-015-1294-0.
[63] J. Zhang, M. Xiao, L. Gao, J. Fu, A novel projection outline based active learning method and its combination with Kriging
metamodel for hybrid reliability analysis with random and interval variables, Comput. Methods Appl. Mech. Engrg. 341 (2018)
32–52, http://dx.doi.org/10.1016/j.cma.2018.06.032.
26
X. Liu, L. Gao and M. Xiao Computer Methods in Applied Mechanics and Engineering 415 (2023) 116209

[64] M. Xiao, J. Zhang, L. Gao, S. Lee, A.T. Eshghi, An efficient Kriging-based subset simulation method for hybrid reliability
analysis under random and interval variables with small failure probability, Struct. Multidiscip. Optim. 59 (6) (2019) 2077–2092,
http://dx.doi.org/10.1007/s00158-018-2176-z.
[65] E. Andreassen, A. Clausen, M. Schevenels, B.S. Lazarov, O. Sigmund, Efficient topology optimization in MATLAB using 88 lines of
code, Struct. Multidiscip. Optim. 43 (1) (2011) 1–16, http://dx.doi.org/10.1007/s00158-010-0594-7.
[66] F. Wang, B.S. Lazarov, O. Sigmund, On projection methods, convergence and robust formulations in topology optimization, Struct.
Multidiscip. Optim. 43 (6) (2011) 767–784, http://dx.doi.org/10.1007/s00158-010-0602-y.
[67] A. Kawamoto, T. Matsumori, S. Yamasaki, T. Nomura, T. Kondoh, S. Nishiwaki, Heaviside projection based topology optimization by
a PDE-filtered scalar function, Struct. Multidiscip. Optim. 44 (1) (2011) 19–24, http://dx.doi.org/10.1007/s00158-010-0562-2.
[68] K. Svanberg, The method of moving asymptotes—a new method for structural optimization, Internat. J. Numer. Methods Engrg. 24
(2) (1987) 359–373, http://dx.doi.org/10.1002/nme.1620240207.

27