# Topology optimization of graded truss lattices based on on-the-fly
homogenization

Bastian Telgen  Mechanics & Materials Lab
Department of Mechanical and Process Engineering
ETH Zürich
8092 Zürich, Switzerland
Email<EMAIL_ADDRESS>Ole Sigmund  TopOpt group
Department of Mechanical Engineering
Technical University of Denmark
2800 Kongens Lyngby, Denmark
Email<EMAIL_ADDRESS>Dennis M. Kochmann Corresponding author.  Mechanics
& Materials Lab
Department of Mechanical and Process Engineering
ETH Zürich
8092 Zürich, Switzerland
Email<EMAIL_ADDRESS>
###### Abstract

We introduce a computational framework for the topology optimization of
cellular structures with spatially varying architecture, which is applied to
functionally graded truss lattices under quasistatic loading. We make use of a
first-order homogenization approach, which replaces the discrete truss by an
effective continuum description to be treated by finite elements in a
macroscale boundary value problem. By defining the local truss architecture
through a set of Bravais vectors, we formulate the optimization problem with
regards to the spatially varying basis vectors and demonstrate its feasibility
and performance through a series of benchmark problems in 2D (though the
method is sufficiently general to also apply in 3D, as discussed). Both the
displacement field and the topology are continuously varying unknown fields on
the macroscale, and a regularization is included for well-posedness. We argue
that prior solutions obtained from aligning trusses along the directions of
principal stresses are included as a special case. The outlined approach
results in heterogeneous truss architectures with a smoothly varying unit
cell, enabling easy fabrication with a tunable length scale (the latter
avoiding the ill-posedness stemming from classical nonconvex methods without
an intrinsic length scale).

##  1 Introduction

Cellular materials achieve excellent mechanical properties at a considerably
lower weight than conventional solids. Having been characterized extensively
with regards to the wide material property space they span [1, 2], cellular
metamaterials enable their functional optimization through a careful choice of
the base material(s) and, moreover, the topology and geometry of the cellular
design [3]. The latter particularly applies in hierarchical networks, whose
overall properties depend on the arrangement of structural members at each
level [4]. The rapid development of modern additive manufacturing techniques
has enabled the integration of cellular design principles across length scales
down to the nanoscale [5, 6], which has opened an unprecedented design space.
Cellular solids are also omnipresent in nature (e.g., in bone, wood, shell, or
bamboo), whose key distinguishing feature from most man-made metamaterials is
their non-periodicity stemming from self-assembly processes. The design of
metamaterials, by constrast, typically explores periodic architectures in
engineering design, owing to the simplicity of extracting their effective
properties by homogenization.

Truss lattices are one such subclass of periodic cellular materials, exploited
extensively over the past few decades [7, 8, 9]. A multitude of architected
truss (meta-)materials has been introduced theoretically and experimentally
with as-designed properties. Classical examples include extreme properties
such as designs with high-strength- and/or stiffness-to-weight ratios [10],
further unconventional properties such as auxeticity [11, 12] or near-infinite
bulk-to-shear modulus ratios [12, 13, 14, 15]. Advancing from linear to
nonlinear and dynamic material behavior, properties of interest have included
energy absorption [16], acoustic wave tuning [17, 18], controllable nonlinear
stress-strain behavior and shape recoverability under large deformations [19,
6], insensitivity to imperfections [20, 21, 22], and damage tolerance [23].

The almost infinite freedom in engineering the topology of cellular materials
and the resulting enormous potential reflected by the above examples
highlights the challenging inverse problem. Required are systematic approaches
to obtain structures with desired properties. When it comes to optimizing the
effective material behavior of periodic networks, a multitude of analytical
[3], numerical [12, 24, 25] and empirical methods [26] have been deployed,
successfully reaching, e.g., theoretical limits such as stiffness upper bounds
[27, 28] or demonstrating extreme band gaps in their dispersion relations [29,
30]. The integration of such periodic architectures into real-world
applications (from lightweight space structures to efficient infills or
patient-specific medical devices) is another challenge, since complex
macroscale load conditions and constraints must be accounted for [9, 31, 32],
so uniform periodic structures are not a solution, and structural-level
topology optimization becomes a method of choice.

While topology optimization has evolved into a mature systematic design
approach for structural optimization, the classical approach faces high
computational costs in case of large-scale optimization problems involving
several length scales. The recent benchmark presented in [33] proved the
impressive capability of a modern high-performance implementation to optimize
large structures in detail (based on more than one billion finite elements).
If unregularized, the optimization problem lacks existence of solutions,
meaning that increasing resolution leads to the appearance of finer and finer
features (for a detailed summary on regularization approaches and
manufacturability constraints, the reader is referred to [34]). In fact, for
linear elastic problems with isotropic phases and multiple loads, the optimal
local solution is given by laminations at multiple length-scales (3 in 2D and
6 in 3D) [35]. However, in a recent 2D microstructure optimization study [36]
it was demonstrated that single-scale realizations can be obtained at low loss
in performance. So-called homogenization-based approaches that use knowledge
of optimal microstructures offer great potential for time savings. Such
procedures require a separation of scales between the macroscale boundary
value problem and the fine-scale microstructure. However, experience [37]
shows that the optimal local material orientations resulting from such
approaches ensure stretching-dominated microstructures that are less prone to
scale separation issues.

A number of hierarchical multiscale topology optimization schemes have been
introduced, which optimize a macroscale problem while accounting for the
microscale topology by homogenization [38, 39, 40, 41] or computational FE2
techniques [42]. Here, the optimization problem on the macroscale depends on a
parametrization of the structural architecture on the lower scale, and the
latter is optimized in a spatially varying fashion. Such schemes have often
led to locally optimized (and spatially varying) unit cells (e.g., at the
quadrature-point level) without ensuring proper compatibility and connectivity
across the macroscale body, thus also lacking manufacturability. Different
solutions have been proposed to assure the connectivity of spatially varying
lattice topologies. In [43] and [44] a square-/tetrahedron-based topology
description linked to a constant quadrilaterial/tetrahedral mesh of the design
domain in 2D and was used, respectively, with varying unit cell densities. An
appealing generalization to allow for microstructural topology changes is the
introduction of microstructural families, in which all members are
interconnectable [45, 46]. In [47] non-designable connector cells to link
adjacent cells were defined. In [48] and [49] a connectivity index used as a
constraint during optimization was introduced. In [50] different
microstructures were stitched to have a better match, and in [51] transition
layers that are optimized for optimal connectivity were introduced. All of
these approaches produce elegant structures; however, the fixes come at the
expense of a restricted design space and, more severely, the way different
unit cells are connected may violate the homogenization assumptions of scale
separation and periodicity.

An elegant alternative is the introduction of microstructures that are
continuously graded in space – analogous to the transition from classical
multi-material topology optimization defining a unique material at each point
(e.g., [52, 53, 54]) to a continuous grading of materials [55]. In [56] a
framework for free material optimization with constrained grading was
presented and demonstrated its performance for rank-2 laminates;
unfortunately, this approach lacked a projection onto a finite length scale.
Recently the approaches shown in [37], [57], [58] and [59], in line with
earlier work in [60], were used to construct spatially variant lattices on a
finite length scale, based on a rectangular microstructure representation for
which a homogenization-based surrogate model was used. A Similar hybrid
approach has also been applied in 3D as presented in [61]. A recent review of
such multi-scale topology optimization approaches can be found in [62].

In summary, identifying optimal structural architectures with fine-scale
features for complex macroscale boundary value problems can be tackled by
either high-performance topology optimization schemes on the macroscale (but
those come with significant expenses) or by two-scale optimization approaches
(which often lack compatibility or a highly limited design space). We here
present a computational framework which aims at establishing compatible,
manufacturable truss networks with a tailored length scale suitable for
complex macroscale problems. Based on a chosen microstructure parametrization,
we employ a semi-analytical homogenized continuum model for truss lattices
accounting for finite rotations [63], which is used within a macroscale finite
element setup (using both translations and rotations as primary fields). Based
on the assumption of a continuously varying truss architecture on the
macroscale, our framework allows for the optimization of spatially variant
cellular lattice structures. A projection of the obtained set of truss unit
cells at the macroscale onto a chosen finite length scale [64] results in the
final truss design, thereby assuring manufacturability and respecting the
assumptions of periodicity and scale separation. The general framework is
applicable to, in principle, any periodic cellular material and hence enables
the structural realization, combination and symbiosis of a vast range of
microstructural topologies within a single body. We demonstrate the approach
by seeking functionally-graded truss architectures having minimal compliance.
Inspired by the atomic arrangement in crystals, we parametrize the truss unit
cell via a set of Bravais lattice vectors along with Voronoi tessellation.
Overall, our approach may be viewed as a generalization of those introduced
previously in [37], [65], and [57]. Other recent approaches in the field
include those in [66] and [67].

In the following Section 2, we summarize the optimization methodology in four
steps. First, Section 2.1 introduces the microstructure parametrizaion,
followed by a brief summary of the homogenization approach in Section 2.2 and
of the topology optimization problem in Section 2.3. Lastly, the projection
onto the finite length scale is presented in Section 2.4. To demonstrate the
performance of this framework, Section 4 presents a number of benchmarks for
functionally graded truss lattices with minimized compliance. Finally, Section
5 concludes this work and discusses its benefits and limitations.

##  2 Methodology

The present optimization methodology is based on a two-scale description of
the cellular body, which separates the macroscale boundary value problem from
the underlying discrete structural architecture. By assuming a separation of
scales, the small-scale truss network is condensed into an effective,
homogenized continuum model to be used at the larger scale. The connecting
link between both is the definition of the truss topology, which is defined on
the macroscale in a spatially varying fashion. Like the displacement field as
the primary unknown of the macroscale boundary value problem, the truss
topology is continuously varied, and the local unit cell at any point on the
macroscale is the basis for the homogenized mechanical response at that point.
Figure 1 schematically illustrates the two-scale setup to be explained in the
following.

Figure 1: The discrete spatially variant structure is described on the
macroscale as a mechanically equivalent continuous body, whose local effective
constitutive behavior derives from the local definition of the truss topology,
parameterized by the set $\bm{\psi}$ of all spatially varying topology
variables and based on the assumption of a (quasi-)periodic lattice. The
continuously varying topology variables $\bm{\psi}$ are optimized over the
macro-domain $\Omega$ towards a minimum compliance for given boundary
conditions (considering Neumann and Dirichlet boundaries $\delta\Omega_{N}$
and $\delta\Omega_{D}$, respectively). The final discrete truss representation
is obtained by introducing a mapping $\tau^{\prime}$, which projects the
continuous topology variable field $\bm{\psi}$ onto a discrete structure of
finite length scale.

###  2.1 Spatially variant truss topology parametrization

Figure 2: Construction of a unit-volume RUC in an $\varepsilon$-neighborhood
$B_{1}$ around point $\bm{X}\in\Omega$. Based on the Bravais lattice
parameters $\alpha$, $\beta$, and $\gamma=|\bm{a}_{2}|/|\bm{a}_{1}|$ the
lattice sites (red), and hence the Wigner-Seitz cell $\Omega_{\text{WSC}}$,
are uniquely defined. We choose to place nodes $\mathcal{B}_{n}^{\text{RUC}}$
and beams $\mathcal{B}_{b}^{\text{RUC}}$ (blue) on the RUC’s corners and
edges, respectively. All beam properties are collected in the set
$\bm{\lambda}$.

Figure 3: Effective directional Young’s modulus vs. orientation for three
different lattices at $\rho^{RUC}=0.1$: the tetragonal lattice with a square
RUC (dotted red, $\alpha=0$, $\beta=\pi/2$, $\gamma=1$), the orthorhombic
lattice with a rectangular RUC (solid green, $\alpha=0$, $\beta=\pi/2$,
$\gamma=2$), and the hexagonal lattice with a hexagonal RUC (dash-dotted blue,
$\alpha=0$, $\beta=\pi/3$, $\gamma=1$).

For the geometric description of spatially variant truss lattices, we borrow
principles from (spatially invariant) crystal structures in solid state
physics [68]. Confining ourselves to simple Bravais lattices [69], each site
in an infinite lattice in $d$ dimensions is located at position

$\bm{r}^{\prime}=\bm{r}+\sum_{i=1}^{d}n_{i}\bm{a}_{i}\qquad\textrm{with}\qquad
n_{i}\in\mathbb{Z},$ (1)

where $\bm{r}$ denotes a reference point, $\\{\bm{a}_{1},\ldots,\bm{a}_{d}\\}$
are the Bravais basis of linearly independent unit vectors, and
$\\{n_{1},\ldots,n_{d}\\}$ is a unique set of integers. Towards the
description of spatially variant structures, we now allow the primitive
translation vectors to vary smoothly in space on the macroscale, introducing
the (at least $\mathcal{C}^{1}$-continuous) primitive translation vector
fields $\bm{a}_{i}(\bm{X})$ ($i=1,\ldots,d$), where $\bm{X}\in\Omega$ denotes
a position on the macroscale body $\Omega$ in the undeformed configuration. In
an $\varepsilon$-neighborhood of a point $\bm{X}\in\Omega$, the local lattice
site locations are now approximated by

$\bm{r}^{\prime}\approx\bm{X}+\sum_{i=1}^{d}n_{i}\bm{a}_{i}(\bm{X})\qquad\textrm{with}\qquad
n_{i}\in B_{\varepsilon}[0]\subseteq\mathbb{Z},$ (2)

where $B_{\varepsilon}[x]=\\{y\in\mathbb{Z}:|x-y|\leq\varepsilon\\}$ defines a
closed $d$-ball of integer radius $\varepsilon>0$. Since the approximation (2)
is only suitable if the basis vectors and hence the lattice spacing must be
significantly smaller than any macroscopic fields of interest
($|\bm{a}_{i}|\ll|\Omega|$), we require a separation of scales between the
fields describing the structural topology (here, $\bm{a}_{i}(\bm{X})$) and all
macroscopic mechanical fields of interest (such as the displacement field).
The topology variations are consequently so smooth on the macroscale that at
each point $P\in\Omega$ the underlying structure within an
$\varepsilon$-neighborhood can be approximated by a periodic lattice.
Naturally, a spatially constant structure (i.e., a perfect lattice) is
recovered on the macrscoale if $\bm{a}_{i}(\bm{X})=\text{const}.$ for all
$\bm{X}\in\Omega$, in which case (2) is exact for $\varepsilon\to\infty$ and
coincides with (1).

The above Bravais construction only defines a point set. In order to arrive at
a truss topology, one could potentially connect those sites by beams – this
approach, however, is non-unique and rather restrictive in parametrizing the
structural topology. Therefore, we instead use the Wigner-Seitz
construction111We point out that this choice is only one possible option among
various others (e.g., one could extend the unit cell definition to include
several particles, as originally done for multi-lattices in [68]). The general
approach laid out here is independent of the specific choice of topology
parametrization. [70] – based on the set of lattice sites – to define a unique
primitive unit cell topology for a given set of Bravais basis vectors
$\\{\bm{a}_{1},\ldots,\bm{a}_{d}\\}$. The Wigner-Seitz cell
$\Omega_{\textrm{WSC}}$ is defined as the locus of points in space closer to
one lattice site than to any other lattice site (see Figure 3). Its geometric
construction can conveniently and efficiently be carried out by Voronoi
tessellation over the neighborhood of lattice sites defined in (2) around any
point $\bm{X}\in\Omega$. The Wigner-Seitz cell serves as representative unit
cell (RUC) $\Omega_{\textrm{RUC}}=\Omega_{\textrm{WSC}}$ (see Figure 3), whose
$n_{n}^{\textrm{RUC}}$ outer vertices and $n_{b}^{\textrm{RUC}}$ edges are
identified as nodes and beams of a truss, respectively (and we denote the sets
of all nodes and of all beams in the RUC by $\mathcal{B}_{n}^{\textrm{RUC}}$
and $\mathcal{B}_{b}^{\textrm{RUC}}$, respectively). This completes the unique
definition of the truss unit cell topology based on a given set of (local)
Bravais basis vectors.

For convenience, we represent the Bravais lattice vectors in fractional
coordinates, i.e., in 2D
$\bm{a}_{1}=(|\bm{a}_{1}|\cos\alpha,|\bm{a}_{1}|\sin\alpha)$ and
$\bm{a}_{2}=(|\bm{a}_{2}|\cos(\alpha+\beta),|\bm{a}_{2}|\sin(\alpha+\beta))$,
where angles $\alpha$ and $\beta$ define, respectively, the orientation and
the relative orientation of the two Bravais vectors (see Figure 3). The scale
separation allows us to define the homogenization problem on any scale. To
ensure RUC uniqueness, we enforce the representative unit cell to have a
normalized volume of
$|\Omega_{\text{RUC}}|=|\text{det}(\bm{a}_{1},\ldots,\bm{a}_{d})|=1$. The
volume constraint reduces the number of independent Bravais lattice variables
to $d^{2}-1$. In 2D, we choose $\alpha$, $\beta$, and
$\gamma=|\bm{a}_{2}|/|\bm{a}_{1}|$.

Since the mechanical response of a truss depends not only on topology but also
on its material parameters and strut geometry, we here assume that all beams
are of rectangular cross-section with constant out-of-plane thickness $t$ and
spatially varying width $w(\bm{X})$. (It is straightforward to relax this
constrained by introducing additional geometric design variables.) The beam
width $w$ is a direct consequence of the effective density
$\bar{\rho}(\bm{X})$ through $w=\bar{\rho}(\bm{X})|\Omega_{\text{RUC}}|/L$.
All beams are assumed to be composed of the same isotropic linear elastic
material with Young’s modulus $E$ and mass density $\rho$. The beam properties
are summarized in the set $\bm{\lambda}(\bm{X})=\\{E,\rho,t,w(\bm{X})\\}$,
which in the spatially variant context varies with position $\bm{X}\in\Omega$
on the macroscale, analogous to the basis vectors. From now on, we refer to
the set of all design variables, uniquely defining the associated RUC
$\Omega_{\text{RUC}}$ as well as $\mathcal{B}_{n}^{\textrm{RUC}}$ and
$\mathcal{B}_{b}^{\textrm{RUC}}$, as

$\bm{\psi}(\bm{X})=\\{\bm{a}_{1}(\bm{X}),\ldots,\bm{a}_{d}(\bm{X});\bm{\lambda}(\bm{X})\\}.$
(3)

Any further quantities required for the subsequent optimization problem, such
as the local mass density $\rho(\bm{X})$ or the effective mechanical response
(of an equivalent periodic structure), can be extracted from a given set
$\bm{\psi}(\bm{X})$ via homogenization, as described in the following.

###  2.2 Effective continuum description based on first-order homogenization

The topology optimization problem requires information about the local
mechanical constitutive behavior on the macroscale, which must depend on the
truss architecture on the microscale. Our natural solution is a first-order
homogenization approach, exploiting the assumed separation of scales and
extracting the effective response from the RUC. We here employ the semi-
analytical approach in [63], whose framework provides an efficient description
of periodic truss lattices of general linear elastic beam architectures,
accounting for finite rotations and hence geometric nonlinearity. In the
following, we summarize the homogenization technique to the extent necessary
for subsequent discussions and discuss its validity for the investigated type
of spatially variant trusses. For further details of the model and its
numerical implementation we refer to [63].

Acknowledging that the deformation of Timoshenko-Ehrenfest beams is described
by both rotations and translations, we employ a macroscale kinematic continuum
description of micropolar type, introducing two primary fields on the
macroscale: a deformation mapping $\bm{x}=\bm{\varphi}(\bm{X})$ and a rotation
field $\bm{\theta}=\bm{\theta}(\bm{X})$, where $\bm{x}\in\varphi(\Omega)$ and
$\bm{X}\in\Omega$ are the deformed and undeformed positions of a material
point. Assuming the deformation fields vary smoothly (at least
$\mathcal{C}^{1}$-continuous), the local deformation within an infinitesimal
neighborhood of a point $\bm{X}\in\Omega$ is approximated to first order222For
reduced computational complexity, we choose to truncate the Taylor series to
linear order. A detailed discussion and a comparison with second-order
approximations can be found in [63]. by

$\begin{split}x_{i}(\bm{X}+\;\\!\mathrm{d}\bm{X})=\varphi_{i}(\bm{X})+F_{iJ}(\bm{X})\,\;\\!\mathrm{d}X_{J}+\text{h.o.t.},\\\
\theta_{i}(\bm{X}+\;\\!\mathrm{d}\bm{X})=\theta_{i}(\bm{X})+\kappa_{iJ}(\bm{X})\;\\!\mathrm{d}X_{J}+\text{h.o.t.},\end{split}$
(4)

with deformation gradient $\bm{F}(\bm{X})=\nabla\bm{\varphi}(\bm{X})$ and
curvature tensor $\bm{\kappa}(\bm{X})=\nabla\bm{\theta}(\bm{X})$. Under the
assumption of a separation of scales, we apply (4) in the neighborhood of
$\bm{X}\in\Omega$ to the nodes in the truss RUC centered at $\bm{X}$ and
defined by $\bm{\psi}(\bm{X})$. Consequently, every node
$\alpha\in\mathcal{B}_{n}^{\textrm{RUC}}$ initially located at
$\bm{X}_{\alpha}^{m}\in\Omega_{\textrm{RUC}}$ is deformed according to

$\begin{split}\bm{x}_{\alpha}^{m}=\bm{\varphi}(\bm{X})+\bm{F}(\bm{X})\bm{X}_{\alpha}^{m}+\delta\bm{x}_{\alpha}+\text{h.o.t.},\\\
\bm{\theta}_{\alpha}^{m}=\bm{\theta}(\bm{X})+\bm{\kappa}(\bm{X})\bm{X}_{\alpha}^{m}+\delta\bm{\theta}_{\alpha}+\text{h.o.t.}.\end{split}$
(5)

Here and in the following, the superscript m denotes microscale quantities.
Importantly, unlike the affine deformation in (4), the approximations (5)
introduce shifts
$\tilde{\bm{\delta}}_{\alpha}=(\delta\bm{x}_{\alpha},\delta\bm{\theta}_{\alpha})$,
which allow nodes to deform relative to the RUC average in a periodic fashion.
By defining one set of shifts for each unique periodic node in the unit cell,
this effectively imposes periodic boundary conditions for both nodal
translations and rotations and hence ensures that the average primary fields
($\bm{\varphi}$ and $\bm{\theta}$) and their gradients ($\bm{F}$ and
$\bm{\kappa}$) are imposed as averages across the RUC. The relaxation from an
affine to a general periodic deformation of the RUC proves to be especially
important for bending-dominated truss architectures, in which soft modes are
accommodated primarily by bending in a manner far from affine deformation
[63].

Based on the co-rotational beam description in [71], the strain energy of each
beam in the RUC is a function of the nodal displacements and rotations, such
that an effective strain energy density can be defined as the RUC average

$W(\bm{\psi};\bm{F};\bm{\theta},\bm{\kappa};\tilde{\bm{\delta}})=\frac{1}{V_{\textrm{RUC}}}\sum_{\nu\in\mathcal{B}^{\textrm{RUC}}_{b}(\bm{\psi})}W_{\nu}(\tilde{\bm{x}}_{\nu_{1}}^{m},\tilde{\bm{x}}_{\nu_{2}}^{m}),$
(6)

where $V_{\textrm{RUC}}=|\Omega_{\textrm{RUC}}|$ is the volume (or area in 2D)
of the RUC, and $W_{v}$ is the strain energy of beam $\nu$ having the two
nodes $\nu_{1}$ and $\nu_{2}$ with degrees of freedom
$\tilde{\bm{x}}_{\alpha}=(\bm{x}^{m}_{\alpha},\bm{\theta}^{m}_{\alpha})$ from
(5) for $\alpha=\nu_{1},\nu_{2}$. We note that the resulting energy density
$W$ is effectively independent of $\bm{\varphi}$, since $\bm{\varphi}(\bm{X})$
in (5) merely constitutes rigid-body translations, which do not affect the
strain energy (whereas $\bm{\theta}(\bm{X})$ implies nodal rotations and beam
bending, hence not implying rigid-body rotations).

For linear elasticity, the shifts are expected to minimize the unit cell
energy (which is local in the shifts), so that we may condense them out on the
microscale. The resulting condensed energy density may be stated as

$W^{*}(\bm{\psi};\bm{F};\bm{\theta},\bm{\kappa})=\min_{\tilde{\bm{\delta}}}\big{\\{}W(\bm{\psi};\bm{F};\bm{\theta},\bm{\kappa};\tilde{\bm{\delta}})\
:\ \text{constr.}\big{\\}},$ (7)

with constraints ensuring that the shifts $\tilde{\bm{\theta}}$ do not violate
the imposition of the average kinematic variables from the macroscale (see
[63]). (7) defines a generally nonlinear optimization problem to be solved
numerically in an iterative fashion333Although this minimization must be
carried out numerically in general, all first and second derivatives with
respect to the kinematic variables can be computed in closed-form, thus
allowing for the efficient use of consistent tangents [63].. With the above
effective energy density $W^{*}$ (dependent on the local unit cell
architecture $\bm{\psi}$), we may employ principles of continuum mechanics to
derive, e.g., the first Piola-Kirchhoff stress tensor as $\bm{P}=\partial
W^{*}/\partial\bm{F}$ as well as analogous couple stress tensors associated
with the remaining macroscale fields, as well as consistent stiffness matrices
such as $\mathbb{C}=\partial^{2}W^{*}/\partial\bm{F}\,\partial\bm{F}$.

Finally, for a given field $\bm{\psi}(\bm{x})$ defining the local truss
architecture, we can solve the macroscale boundary value problem on $\Omega$,
viz.

$\begin{split}(\bm{\varphi},\bm{\theta})_{\bm{\psi}}&=\arg\inf\left\\{\int_{\Omega}W^{*}(\bm{\psi};\bm{F},\bm{\theta},\bm{\kappa})\,\;\\!\mathrm{d}V-\int_{\partial\Omega_{T}}\hat{\bm{T}}\cdot\bm{\varphi}\;\\!\mathrm{d}S\right.\\\
&\qquad\qquad\qquad\qquad\left.-\int_{\partial\Omega_{m}}\hat{\bm{m}}\cdot\bm{\theta}\;\\!\mathrm{d}S\
:\ \bm{\varphi}=\hat{\bm{\varphi}}\ \text{on}\ \partial\Omega_{u},\
\bm{\theta}=\hat{\bm{\theta}}\ \textrm{on}\
\partial\Omega_{\theta}\right\\}\end{split}$ (8)

with tractions $\hat{\bm{T}}$ acting on the boundary
$\partial\Omega_{T}\subset\partial\Omega$, and moments $\hat{\bm{m}}$ imposed
on $\partial\Omega_{m}\subset\partial\Omega$. For conciseness, we introduce
the generalized deformation mapping
$\tilde{\bm{\varphi}}(\bm{X})=(\bm{\varphi}(\bm{X}),\bm{\theta}(\bm{X}))$
along with $\tilde{\bm{F}}=\bm{\nabla}\tilde{\varphi}$ and
$\hat{\bm{t}}=(\hat{\bm{T}},\hat{\bm{m}})$, such that (8) simplifies to

$\tilde{\bm{\varphi}}=\arg\inf\left\\{\int_{\Omega}W^{*}(\bm{\psi};\tilde{\bm{\varphi}},\tilde{\bm{F}})\,\;\\!\mathrm{d}V-\int_{\partial\Omega_{N}}\hat{\bm{t}}\cdot\tilde{\bm{\varphi}}\;\\!\mathrm{d}S:\tilde{\bm{\varphi}}=\hat{\bm{\varphi}}\textrm{
on }\partial\Omega_{D}\right\\}$ (9)

with the combined Dirichlet boundary $\partial\Omega_{D}$.

###  2.3 Topology optimization

In topology optimization we classically seek to find the optimal material
distribution that minimizes an objective function under given constraints.
Here, the global design field is $\bm{\psi}(\bm{X})$, which we define to be
varying continuously across $\Omega$ (or its finite element discretization
$\Omega_{h}$). In our examples, we focus on classical benchmarks that minimize
the compliance for single-load-case problems while enforcing a total volume
(or mass) constraint. However, the framework presented here is sufficiently
general to be extended to more complex objective functions and, specifically,
it is already laid out for nonlinear problems (since the homogenization scheme
accounts for finite beam rotations).

The optimization problem is formalized by

$\begin{split}\min_{\bm{\psi}}:&\quad\mathcal{L}(\bm{\psi};\tilde{\bm{\varphi}})=\int_{\Omega}\ell(\bm{\psi};\tilde{\bm{\varphi}})\;\\!\mathrm{d}V=\int_{\Omega}\left[W^{*}(\bm{\psi};\tilde{\bm{\varphi}},\tilde{\bm{F}})+c_{p}||\bm{\nabla}\psi_{i}||_{p}\right]\;\\!\mathrm{d}V\\\
\textrm{s.t.:}&\quad\int_{\Omega}W^{*}(\bm{\psi};\tilde{\bm{\varphi}},\tilde{\bm{F}})\,\;\\!\mathrm{d}V=\int_{\partial\Omega_{S}}\hat{\bm{t}}\cdot\tilde{\bm{\varphi}}\;\\!\mathrm{d}S,\\\
:&\quad\int_{\Omega}\rho(\bm{\psi})\;\\!\mathrm{d}V-V_{\textrm{max}}\leq 0,\\\
:&\quad\text{RUC topology constraints},\end{split}$ (10)

where $||(*)||_{p}$ denotes the $\bm{L}^{p}$ norm (in the following, we use
$p=2$), and the $\bm{L}^{p}$-constraint [72] governs the smoothness of the
global design variable field $\bm{\psi}(\bm{X})$ through a regularization
constant $c_{p}\geq 0$. $V_{\textrm{max}}$ defines the maximum material volume
in the domain $\Omega$ (in subsequent examples we choose
$\rho^{\Omega}=V_{\textrm{max}}/|\Omega|=0.2$), and the “RUC topology
constraints” stand for any constraints imposed on the local truss topology and
architecture (such as constraints on the unit cell topology, beam geometry,
etc.). For the design updates of $\bm{\psi}$ we use the method of moving
asymptotes [73]. Due to its nonlinear nature, the above optimization problem
is solved iteratively, minimizing the objective with respect to the design
variables $\bm{\psi}$ in a step-wise fashion. At each iteration step, the
multiscale finite element problem defined in (2.2) is solved for
$\tilde{\varphi}$ in equilibrium. Computing the increments of the design field
$\bm{\psi}$ requires computing its gradient, which is classically termed the
sensitivity and given by

$\frac{\partial\ell}{\partial\psi_{i}}=\frac{\partial
W^{*}}{\partial\psi_{i}}+c_{p}\frac{\bm{\nabla}\psi_{i}|\bm{\nabla}\psi_{i}|^{p-2}}{||\bm{\nabla}\psi_{j}||^{p-1}_{p}},$
(11)

where the semi-analytical homogenization scheme introduced in Section 2.2
allows us to efficiently calculate the change of the condensed energy density
$\partial W^{*}$, locally in each quadrature point, with respect to a design
variable change $\partial\psi_{i}$ through a local forward-difference scheme
(introducing a small perturbation on each design variable of $10^{-9}$).

(a) Continuum cantilever beam problem as defined in [74]. All nodes on the
left edge are clamped. The load $P=10^{-6}$ is applied as distributed load
over the central $20\%$ of the right edge. Solid black indicates $50\%$ volume
fraction.

(b) Discrete cantilever beam. All nodes with distance to the left edge smaller
than $0.1h$ are clamped. The load $P$ is distributed to all nodes with
distance smaller $0.1h$ to the centerright point. This image serves for
visualization purposes only.

Figure 4: The cantilever beam is well studied and serves as benchmark. In this
study, the beam domain $\Omega$ has a length $L=2$ and a height $h=1$.
Furthermore, we choose the load $P=10^{-6}$.

For a performance comparison to literature benchmarks, we restrict our study
to the cantilever beam in 4(a) with an applied load of magnitude $P$ and
define the beam compliance, normalized by the effective density
$\rho^{\Omega}$, as

$\mathcal{C}_{C}(\tilde{\bm{x}})=\frac{\rho^{\Omega}}{P^{2}}\int_{\Omega}W^{*}(\bm{\psi};\tilde{\bm{\varphi}},\tilde{\bm{F}})\;\\!\mathrm{d}V.$
(12)

In subsequent examples, we study three different RUC families in $d=2$
dimensions (2D), whose truss architectures are constrained by the RUC topology
constraints in Table 1. The square lattice family resembles the study in [37]
and serves as a benchmark. The rectangular lattice family, and ultimately the
monoclinic Bravais lattice, exemplify two generalizations facilitated by this
approach. Note that the effective density $\bar{\rho}$ is analogous to
classical density-based topology optimization, while the Bravais lattice
variables ($\alpha$, $\beta$, $\gamma$) allow for tailoring of the local
anisotropy of the effective elastic surfaces. By letting the effective density
vary up to $50\%$, we are stretching the assumption of slender beams. Please
find a detailed discussion and reasoning of the upper density bound in A. As
initial guess we use a homogeneous design variable field with values in the
center of the design variable ranges.

Table 1: RUC topology constraints for tetragonal and orthorhombic crystal lattice families, and a “monoclinic” Bravais lattice parametrization that includes all primitive Bravais lattice families in 2D (tetragonal, hexagonal, orthorhombic, monoclinic). While the tetragonal lattice family is widely used in the literature, other lattice families are oftentimes neglected, which restricts the design space. tetragonal | orthorhombic | monoclinic
---|---|---
$-\pi\leq\alpha\leq\pi$ | $-\pi\leq\alpha\leq\pi$ | $-\pi\leq\alpha\leq\pi$
$\beta=\pi/2$ | $\beta=\pi/2$ | $0\leq\beta\leq\pi$
$\gamma=1$ | $0.5\leq\gamma\leq 2$ | $0.5\leq\gamma\leq 2$
$10^{-5}\leq\bar{\rho}\leq 0.5$ | $10^{-5}\leq\bar{\rho}\leq 0.5$ | $10^{-5}\leq\bar{\rho}\leq 0.5$

###  2.4 From the continuous solution to a discrete structure

Solving the optimization problem (10) results in a continuous distribution of
the truss design parameter field $\bm{\psi}$ across $\Omega$ on the
macroscale, whose underlying microstructure is, in principle, infinitely fine
based on the separation of scales. In order to construct a discrete truss at a
finite length scale (based on the continuous optimal field $\bm{\psi}$
obtained from (10)), we here extend the approach presented in [64].

We introduce a mapping $\epsilon(\bm{X}):\Omega\to[0,1]$, such that a value of
$0$ refers to one phase (e.g., void) and a value of $1$ to another phase
(e.g., solid). Most generally, we represent $\epsilon(\bm{X})$ of a periodic
unit cell in terms of a Fourier series (e.g., we find the complex amplitudes
for a discretely defined unit cell through a fast Fourier transform). When
assuming translational symmetry such as in a lattice, we use the symmetry
information and represent the infinite structure in terms of the Fourier
series over reciprocal lattice vectors, most generally written as

$\epsilon(\bm{X})=\sum_{\bm{k}}b_{\bm{k}}e^{i2\pi\mu\bm{k}\cdot\bm{X}}\quad\textrm{with}\quad\bm{k}=\sum_{i=1}^{d}m_{i}\bm{k}_{i}$
(13)

being an integer combination of the primitive reciprocal lattice vectors
$\bm{k}_{i}$, and $b_{\bm{k}}$ denoting the complex amplitude related to the
spatial harmonic with reciprocal lattice vector $\bm{k}$. $\mu$ represents a
dimensionless length scale introduced to govern the unit cell size.
Translational symmetry is ensured by the primitive reciprocal lattice vectors
which, in $\mathbb{R}^{3}$, are given by

$\bm{k}_{1}=\frac{\bm{a}_{2}\times\bm{a}_{3}}{\bm{a}_{1}\cdot(\bm{a}_{2}\times\bm{a}_{3})},\qquad\bm{k}_{2}=\frac{\bm{a}_{3}\times\bm{a}_{1}}{\bm{a}_{2}\cdot(\bm{a}_{3}\times\bm{a}_{1})},\qquad\bm{k}_{3}=\frac{\bm{a}_{1}\times\bm{a}_{2}}{\bm{a}_{3}\cdot(\bm{a}_{1}\times\bm{a}_{2})},$
(14)

based on the primitive Bravais lattice vectors $\bm{a}_{i}$. In 2D, the
reciprocal lattice vectors $\bm{k}_{1}$ and $\bm{k}_{2}$ follow from (14) with
$\bm{a}_{3}=(0,0,1)$.

When constructing a spatially variant truss, the spatial variation of
$\bm{k}(\bm{X})$ renders (13) problematic, since the latter requires a
constant history of the frequency vectors $\bm{k}$ in space. When violated,
highly distorted projections arise. As a remedy, in [64] a projection function
$\phi_{\bm{k}}(\bm{X})$ was introduced to obtain the mapping

$\epsilon^{\prime}(\bm{X})=\sum_{\bm{k}}b_{\bm{k}}e^{i2\pi\mu\phi_{\bm{k}}(\bm{X})}$
(15)

with lattice symmetry enforced locally through

$\nabla\phi_{\bm{k}}(\bm{X})=\bm{k}(\bm{X}).$ (16)

That is, we ensure that $\nabla\phi_{\bm{k}}$ is locally tangential to the
(local) $\bm{k}$-vector defined through (14) at every point on the macroscale.

Finding the field $\phi:\Omega\to\mathbb{R}^{d}$, in principle, requires the
solution of an over-determined system. When the lattice is perfectly periodic,
the system has the solution $\phi_{\bm{k}}(\bm{X})=\bm{k}\cdot\bm{X}+c$ with
arbitrary $c\in\mathbb{R}$. By contrast, for all spatially variant trusses, we
find $\phi_{\bm{k}}$ as a least-square minimizer:

$\phi_{\bm{k}}(\bm{X})=\arg\min\int_{\Omega}||\nabla\phi_{\bm{k}}(\bm{X})-\bm{k}||^{2}\;\\!\mathrm{d}\Omega.$
(17)

Given a separation of scales between the spatial variation of the lattice
vectors and the introduced RVE length scale $\mu$, this ensures quasi-
periodicity in the neighborhood of each RVE (on the microscale). The global
residual defined by (17) (which is generally non-zero) can thus be interpreted
as a measure of the violation of scale separation.

For the lattice description introduced in Section 2.1, we define a mapping
function $\tau(\bm{X})$ which attains a value of $1$ at each lattice site
while being $0$ elsewhere:

$\tau(\bm{X})=\prod_{i=1}^{d}\sum_{n=-\infty}^{\infty}\delta\left(\mu\frac{\bm{k}_{i}\cdot\bm{X}}{||\bm{k}_{i}||}-n||\bm{k}_{i}||\right)=\prod_{i=1}^{d}\Sha\left(\mu\bm{k}_{i}\cdot\bm{X}\right),$
(18)

where $\delta(x)$ and $\Sha(x)$ denote the Kronecker delta and Kronecker delta
impulse train, respectively, in $d$ dimensions. Decomposing (18) into spatial
harmonics, the lattice mapping function can be represented as the Fourier
series

$\tau(\bm{X})=\prod_{i=1}^{d}\frac{1}{1+2M}\sum_{m=-M}^{M}e^{i2\pi
m\mu\bm{k}_{i}\cdot\bm{X}}=\prod_{i=1}^{d}\frac{1}{1+M}\left[1+\sum_{m=1}^{M}\cos(2\pi
m\mu\bm{k}_{i}\cdot\bm{X})\right].$ (19)

Note that, in general, an infinite number ($M\to\infty$) of spatial harmonics
is required to represent the Dirac combs exactly. For implementation purposes,
we truncate the above Fourier series and, using the fact that the
translational symmetry is contained in all spatial harmonics, choose $M=1$
such that

$\tau^{\prime}(\bm{X})=\frac{1}{2^{d}}\prod_{i=1}^{d}\left[1+\cos(2\pi\mu\phi_{\bm{k}_{i}}\big{(}\bm{X})\big{)}\right]$
(20)

represents the truncated lattice mapping function with locally enforced
lattice symmetry according to (16). We thus converted the strictly periodic
and discontinuous lattice mapping (18) into a spatially varying and
$\mathcal{C}^{\infty}$-continuous truncated mapping function (20), whose local
maxima represent lattice site locations and have value $1$ for convenience of
implementation. We thus identify the lattice sites within $\Omega$ as the
local maxima of (20) (see 5(b)). Finally, the truss is constructed by
extracting all nodes $\mathcal{B}^{\Omega}_{n}$ and beams
$\mathcal{B}^{\Omega}_{b}$ of width $w(\bm{X})>w_{\text{min}}$ – analogous to
the construction of the RUC described in Section 2.1 – from the Voronoi
diagram over all lattice sites in the global domain $\Omega$ (see 5(c)). The
threshold $w_{\text{min}}$ reflects the lower bound on $\bar{\rho}$ set in
Table 1.

(a) lattice mapping function $\tau^{\prime}(\bm{X})$

(b) lattice sites

(c) Voronoi diagram

Figure 5: The discrete lattice is constructed from the lattice mapping
function by first identifying the lattice sites, e.g., at lattice mapping
function minima, and second constructing the Voronoi diagram over all lattice
sites.

###  2.5 Discrete description

With the effective continuum description introduced in Section 2.2 we
efficiently approximate the lattice behavior, as demonstrated previously for
periodic lattices in [63]. To validate its applicability to the spatially-
variant lattices introduced here, we compare to full discrete simulations.

Based on the global sets of nodes and beams, $\mathcal{B}^{\Omega}_{n}$ and
$\mathcal{B}^{\Omega}_{b}$, respectively (see Section 2.1), we define the
total potential energy of the beam network subjected to generalized external
nodal tractions $\hat{\bm{t}}$ acting on all nodes in
$\mathcal{B}^{\partial\Omega_{N}}_{n}$ as

$I(\tilde{\bm{x}})=\sum_{\nu\in\mathcal{B}^{\Omega}_{b}}W_{\nu}(\tilde{\bm{x}}_{\nu_{1}},\tilde{\bm{x}}_{\nu_{2}})-\sum_{\alpha\in\mathcal{B}^{\partial\Omega_{N}}_{n}}(\hat{\bm{t}}_{\alpha}\cdot\tilde{\bm{x}}_{\alpha}),$
(21)

where $\tilde{\bm{x}}=\\{\tilde{\bm{x}}_{1},\dots,\tilde{\bm{x}}_{n_{n}}\\}$
is the set of all nodal degrees of freedom. Unlike for the continuum
cantilever beam in 4(a), the load $P$ is distributed over a semi-circular area
$\partial\Omega_{N}$, as shown in 4(b), to assure an even load distribution
between discrete nodes. For analogous reasons, the fixed boundary condition at
the left boundary of the domain is spread over a thin slice
$\mathcal{B}^{\partial\Omega_{D}}_{n}$. Classically, the generalized force
vector of node $\alpha$, containing both forces and moments, and the
associated stiffness matrix follow as, respectively,

$\tilde{\bm{t}}_{\alpha}(\tilde{\bm{x}})=\frac{\partial
I}{\partial\tilde{\bm{x}}_{\alpha}}(\tilde{\bm{x}}),\qquad\tilde{\bm{K}}_{\alpha\beta}(\tilde{\bm{x}})=\frac{\partial
I}{\partial\tilde{\bm{x}}_{\alpha}\partial\tilde{\bm{x}}_{\beta}}(\tilde{\bm{x}}).$
(22)

We obtain the quasistatic response from energy minimization, i.e.,

$\tilde{\bm{x}}=\arg\min\\{I:\tilde{\bm{x}}_{\alpha}=\bf 0\ \text{for}\
\alpha\in\mathcal{B}^{\partial\Omega_{D}}_{n}\\},$ (23)

where the generalized displacements $\tilde{\bm{x}}_{\alpha}$ of all nodes in
$\mathcal{B}^{\partial\Omega_{D}}_{n}$ near the left boundary are enforced to
be zero. For comparability with literature and the compliance of the effective
continuum description (12), we define the normalized compliance for the
discrete beam network as

$\mathcal{C}_{D}(\tilde{\bm{x}})=\frac{\rho^{\Omega}}{P^{2}}\sum_{\nu\in\mathcal{B}^{\Omega}_{b}}W_{\nu}(\tilde{\bm{x}}_{\nu_{1}},\tilde{\bm{x}}_{\nu_{2}}),\qquad\text{where}\qquad\rho^{\Omega}=\frac{1}{|\Omega|}\sum_{\nu\in\mathcal{B}^{\Omega}_{b}}w_{\nu}L_{\nu}t\rho$
(24)

with $L_{\nu}$ denoting the length of beam $\nu$.

##  3 Numerical implementation

The mechanical boundary value problem on the macroscale, described by the
variational optimization problem (9) or (23), is solved via a finite-element
approximation which discretizes $\Omega$ into $\Omega_{h}$, using bi-linear
quadrilateral elements to result in bi-linear interpolations of both the
displacement field $\bm{\varphi}_{h}$ and the rotation field $\bm{\theta}_{h}$
across $\Omega_{h}$. To accommodate the micropolar nature of the
discretization, we use an in-house finite-element code, which solves the
stationarity conditions (i.e., the equations of linear momentum balance) via
Newton-Raphson iteration for quasistatic equilibrium. Further details about
the numerical treatment of the homogenized continuum model can be found in
[63].

The design parameter field $\psi(\bm{X})$, defined in (3), is interpolated
using the exact same finite element approximation based on bi-linear
quadrilateral elements, thus providing a continuous field of spatially varying
truss unit cell architectures. As pointed out in Section 2.3, finding the
solution for $\psi(\bm{X})$ involves the method of moving asymptotes. The
overall optimization problem is solved in an iterative fashion, minimizing the
objective function with respect to field $\psi(\bm{x})$, while solving the
macroscale boundary value problem.

To construct the discrete truss in a postprocessing step from the given
interpolated $\psi$-field, we first solve for the projection function (17),
using again the same finite element approximation, and then identify the
lattice sites as local maxima of (20). Only at this step we introduce the
finite length scale $\mu$. Finally, we build the global Voronoi diagram and
extract node locations and beam connections.

##  4 Numerical examples

The introduced methodology is a tool for the optimization of generic
structures based on a free choice of the underlying Bravais RUC topology and
constituents. In the following, we showcase the potential of this approach by
studying the optimization of tetragonal, orthorhombic and monoclinic truss
lattices, and we investigate the influence of the RUC constraints. Moreover,
we discuss the underlying assumption of a separation of scales and demonstrate
how it can be (approximately) fulfilled in the lattice projection as a trade-
off between continuum design variable field regularization and choice of
length scale.

###  4.1 Tetragonal lattice family

The tetragonal truss lattice family (see Table 1) has been extensively studied
in the literature and is the basis of various lattice optimization schemes
(e.g. [43, 48, 37]). Hence, the tetragonal lattice family serves as an
excellent introductory example and benchmark for our method. We first present
optimized continuum design variable results and discuss the influence of mesh
element size and the relaxation constant. Afterwards, we project the
continuous fields onto a finite scale and compare the performance of the
discrete truss lattice with that of the optimized continuum.

#### Continuum

(a) $c_{p}=10^{-13}$

(b) $c_{p}=10^{-12}$

(c) $c_{p}=10^{-11}$

Figure 6: Optimized cantilever beam based on a tetragonal truss lattice for different regularizations. Local Bravais bases are indicated by arrows, and the density is shown by the color code. The arrows are scaled and shown at selected locations for readability. Resolution: 400x200 elements. All plots are rotated by 90∘. Full design variable fields can be found in Figure 19. Table 2: Minimized continuum compliance $\mathcal{C}_{\text{C}}$ for tetragonal, orthorhombic, and monoclinic crystal lattice families with different regularization constants $c_{p}$ and initial guesses. The (H)omogeneous initial guess assumes initial angles $\alpha=\pi/4$, $\beta=\pi/2$, $\gamma=1$, and $\rho=\rho^{\Omega}=0.2$. The (O)rthotropic initial guess uses the optimized fields from the orthotropic lattice family as initial guess. | tetragonal |  | orthorhombic |  | monoclinic
---|---|---|---|---|---
$c_{p}$ | (H) |  | (H) |  | (H) | (O)
$10^{-13}$ | 66.6 |  | 62.2 |  | 76.3 | 62.2
$10^{-12}$ | 67.4 |  | 63.3 |  | 78.4 | 63.3
$10^{-11}$ | 84.8 |  | 79.9 |  | 88.2 | 79.9
Figure 7: Compliance versus mesh resolution for tetragonal lattice cantilever
problem. A mesh of 80x40 elements is within $2\%$ of the converged compliance
$66.6$. $c_{p}=10^{-13}$.

We solve the optimization problem stated in (10), subjected to the tetragonal
RUC constraints (see Table 1), on a continuum mesh of 400x200 Q4 elements. The
mesh resolution was chosen as a trade-off between convergence (see Figure 7),
fine resolution for projection step, and acceptable computation time ($<20$
minutes on a regular laptop until convergence). The infimizing design variable
fields $\bm{\psi}(\bm{X})$ for low ($c_{p}=10^{-13}$), medium
($c_{p}=10^{-12}$) and strong ($c_{p}=10^{-11}$) regularization are presented
in Figure 19. The resulting Bravais bases vectors and density field can be
found in Figure 6. The corresponding compliance values are shown in Table 2.

The regularization, as a constraint on the design variable field, has a
detrimental effect on the effective compliance. The regularization’s influence
on the compliance and density field is illustrated in 8(a) and 8(b),
respectively. While for low regularization, its influence on compliance and
density is negligible (see case (i)), the compliance is affected with
increasing $c_{p}$ (case (ii)). When further increasing, the compliance
degrades strongly (case (iii)) and finally reaches a second plateau (case
(iv)). This upper plateau refers to spatially constant values of $\alpha$ and
$\rho$. Given the negative effect on structural performance, one might wonder
why we use regularization. The reasons are twofold.

(a) Compliance versus regularization constant. For a low regularization
constant the compliance converges to a minimum, while for high regularization
constants the compliance plateaus at its maximum.

(b) Effective density for the optimized cantilever beam based on a tetragonal
lattice for different relaxation constants $c_{p}$, as indicated by the
markers in (a). Discretization: 80x40 elements.

Figure 8: Influence of the regularization constant $c_{p}$.

First, with no or too little regularization, numerical (interpolation-
dependent) artifact patterns arise (see case (i)). The effect is comparable to
the well-known checkerboarding in classical topology optimization with
discrete density fields [75]. The chosen FE discretization leads to
artificially stiff results, and the projection of such laminated lattices onto
a finite scale would result in disconnected lattices. Regularization
alleviates this problem. Alternatively, the problem could be avoided by
increasing the element interpolation order. However, we here choose the
regularization to control the length scale and restrict our finite elements to
linear interpolation for simplicity.

Second, the continuously described RUCs may feature abrupt changes of the
design variables (e.g., when crossing finite element boundaries on the
macroscale). When projecting those onto a finite discrete lattice, this may
either locally break the assumption of quasi-periodicity (see Equation 2) or
require high resolution of the discrete projection (in fact, infinite
resolution in the case of discrete design variable jumps, so $\mu\to 0$). The
regularization constant gives us control over the locally permitted change of
the unit cell in the lattice and hence guarantees a finite length scale
($\mu>0$) of the discrete lattice.

In comparison with optimization results of other topology optimization
schemes, the compliance for the tetragonal lattice with regularizations
$c_{p}=10^{-13}$ ($c_{p}=10^{-12}$, $c_{p}=10^{-11}$) is $34.9\%$ ($36.6\%$,
$71.8\%$) higher than the analytical solution ($49.35$) in [76] and $14.1\%$
($15.5\%$, $45.3\%$) higher than the benchmark ($58.35$) in [37]. This is
expected, as the topology space for the tetragonal unit cell description is
the most limited of this study. It can reproduce the same unit cell topology
as in [37], but it lacks the possibility to tune the directionality by
adjusting the widths of each beam in a unit cell. Moreover, in [77] it was
shown that the optimal orientation of an orthotropic material – like the
tetragonal lattice (see Figure 3) – aligns with the principal stress
directions. Here, $\alpha$ governs the unit cell orientation and indeed aligns
beams with the principle stress directions. One may argue that this knowledge
could be used in the optimization process. However, when the other RUC
parameters ($\beta$, $\gamma$, and $\bm{\lambda}$) are subject to
optimization, $\alpha=0$ does not necessarily point in the direction of
highest stiffness (see Figure 3). In the following sections we will study the
influence of changing orthotropy through changes of $\beta$ and $\gamma$.

#### Discrete lattice

Following the projection approach described in Section 2.4, we construct the
spatially-variant tetragonal lattice from the design variable fields in Figure
19 with different length scales $\mu$. Examples are presented in Figure 10.

To evaluate the performance of the projection method, we calculate the
compliance of the discrete lattice (cf. (24)) subjected to the boundary value
problem depicted in 4(b). In Figure 10 we compare the discrete compliance
$\mathcal{C}_{D}$ to the previously obtained continuum compliance
$\mathcal{C}_{C}$ at different regularization intensities. Note that at a
length scale of $\mu=1/320$ the lattice contains more than 200,000 elements.
Direct optimization of such high-resolution discrete structures would be very
time-demanding, while the desired resolution follows (almost) for free from
this multiscale approach (especially since changing the length scale $\mu$
does not require solving a new topology optimization problem but is merely a
change in postprocessing).

The discrete compliance converges towards the continuum approximation for all
cases. Depending on the regularization strength, the projected lattice
compliance is within $30\%$ of the optimization result for a coarse length
scale of $\mu\approx 1/80$ at $c_{p}=10^{-10}$, $\mu\approx 100$ at
$c_{p}=10^{-11}$, and within $20\%$ of the optimization result for a length
scale of $\mu\approx 1/150$ at $c_{p}=10^{-10}$, and $\mu\approx 1/170$ at
$c_{p}=10^{-11}$. For the lowest regularization ($c_{p}=10^{-13}$), the
discrete compliance is still $>7\%$ higher for the highest considered
resolution ($\mu=1/320$).

The convergence behavior indicates that – for arbitrary regularization – the
discrete compliance will ultimately near the continuum compliance when the
length scale approaches zero ($\mu\to 0$); a small offset is expected to stem
from the inevitable differences in both clamping and loading conditions
between the formulations of the continuum and the discrete problem (see 4(a)
and 4(b)).

Figure 9: Tetragonal lattice cantilever beam at different length scales for
$c_{p}=10^{-12}$. For large values of $\mu$, the structure resembles gradually
the density field in Figure 19. At $\mu=1/160$ the lattice consists of more
than 120,000 elements.

Figure 10: Convergence of the discrete compliance $\mathcal{C}_{D}$ to the
continuum approximation $\mathcal{C}_{C}$ for regularization constants
$c_{p}=10^{-13}$ (solid, blue), $c_{p}=10^{-12}$ (dash-dotted, green), and
$c_{p}=10^{-11}$ (dashed, red).

In summary, optimziation of the tetragonal-lattice cantilever beam reached
minimized compliance values close to those of comparable studies in the
literature – but still above the reported minima. This is attributed to a lack
of directional stiffness tunability in the tetragonal RUC definition and our
choice of $\bm{\lambda}$. Therefore, we proceed to enrich the truss design
space by considering the orthorhombic lattice family.

###  4.2 Orthorhombic lattice family

In addition to tailoring the local unit cell orientation through angle
$\alpha$ as for the tetragonal lattices in Section 4.1, the orthorhombic
lattice family (see Table 1) allows us to tune the anisotropy of the lattice
through design variable $\gamma$, as shown in Figure 3. As before, we first
discuss the continuum solution, before demonstrating discrete lattices.

#### Continuum

Analogous to Section 4.1, we solve the optimization problem (10), now
subjected to orthorhombic RUC constraints (see Table 1), using a continuum
mesh of 400x200 bilinear elements. The resulting design variables for
regularization constants $c_{p}=\\{10^{-11},10^{-12},10^{-13}\\}$ are
presented in Figure 20. Examples are shown in Figure 11. The corresponding
compliance values are included in Table 2.˙

In comparison to the tetragonal lattice, the compliance is up to $7\%$ lower.
Comparing to [37] (where no local volume constraint was used), the result at
$c_{p}=10^{-13}$ is within $6\%$ difference, which is notable given the
constraints on local volume fraction ($\bar{\rho}(\bm{X})\leq 0.5$) shown in
Table 1.

These results are promising, but they should be interpreted with caution, as
some of the beams in the design become relatively thick (as a consequence of
high local densities). This stretches the beam slenderness assumption. We
refer to A for a detailed study of suitable beam slenderness ratios in
compliance with the chosen Timoshenko-Ehrenfest beam model. If we allow for
relatively high local densities, then (as shown in A.2) we accept local beam
stiffness overpredictions of around $10\%$. This may result in potentially
underestimated structural compliance values. As we are primarily presenting an
optimization framework, where homogenization and RUC definition is
interchangeable, respectively, we here accept these uncertainties.

(a) $c_{p}=10^{-13}$

(b) $c_{p}=10^{-12}$

(c) $c_{p}=10^{-11}$

Figure 11: Optimized orthorhombic-lattice cantilever beam for different
regularization constants. Local Bravais bases are indicated by arrows, and the
density is shown by the color code. The arrows are scaled and shown at
selected locations for readability. Resolution: 400x200 elements. All plots
are rotated by 90∘. Full design variable fields can be found in Figure 20.

#### Discrete lattice

Analogous to the tetragonal lattice, we project the design variable fields
presented in Figure 20 onto different finite scales. Results are shown in
Figure 13. As the $\alpha$-fields in Figure 11 indicated, the projections are
similar in orientation to the tetrageonal case. The anisotropy of the
orthorhombic lattice compared to the tetragonal lattice is more pronounced but
equal in directionality (see Figure 3). In [77] it was showon that the optimal
orientation of an orthotropic composite coincides with the principal stress
directions. Therefore for similar orthotropic directionality the orientation
must be same.

In addition, the free design variable $\gamma$ allows to pronounce the
directionality which is strongly utilized towards the supports, where the
principle stresses are high. Results indicate that a larger range for
$\gamma$, allowing for even more strongly pronounced directionality, would
allow for further minimization of the compliance. However, a RUC that spans
more space in one direction than in the other also requires smaller length
scales to approach the continuum compliance. Furthermore, the approximation
concerns regarding the chosen homogenization approach based on slender beams
(discussed in A) increase with decreasing slenderness. As a trade-off we hence
choose to restrict $\gamma$ to the given interval ($0.5\leq\gamma\leq 2$).

Figure 12: Orthorhombic-lattice cantilever beam at different length scales
$\mu$ for $c_{p}=10^{-12}$.

Figure 13: Monoclinic-lattice cantilever beam at different length scales $\mu$
for $c_{p}=10^{-12}$.

###  4.3 Monoclinic lattice family

In an ultimate step, we enrich the set of design variables $\bm{\psi}$ by
additionally allowing the parameter $\beta$ to be optimized (see Table 1) and
study its influence on the truss design and on the effective compliance.
$\beta$ not only pronounces the orthotropy (in addition to $\gamma$), but it
also changes the topology of the elastic surface and introduces asymmetry (see
Figure 3).

#### Continuum

Analogous to the previous lattice families, we solve the optimization problem
(10), subjected to the monoclinic RUC constraints of Table 1, on a continuum
mesh of 400x200 bilinear elements. The optimized design variable fields for
moderate regularization constants $c_{p}=\\{10^{-11},10^{-12}\\}$ are
summarized in Figure 21 and Figure 22. Examples are shown in Figure 14. Both
configurations emerge from the orthorhombic initial guess. The resulting
compliance values are, again, included in Table 2.

In terms of the compliance, the optimized monoclinic lattice performs worse
than the optimized orthorhombic and tetragonal lattices when starting from
homogeneous initial conditions. By contrast, when using the results from the
orthorhombic lattice optimizations as an initial guess, all design variables,
including the monoclinic design parameter $\beta$, remain almost unchanged and
hence the same compliance is attained. The additional design parameter hence
worsens the challenge of finding a good initial guess, since more local minima
of the optimization problem arise and the gradient-based optimization approach
converges into the closest to the initial guess. Therefore, no solution with
better performance than the orthorhombic lattice family could be found. In
fact, some of the local minima can be higher than those of the more
constrained orthorhombic lattice (cf. Table 2). This finding is in line with
previous studies (e.g. [78] and [74]), which suggested that the orthorhombic
lattice family is performing close to theoretical limits. The monoclinic
lattice family includes these results and, in principle, better informed
initial guesses may lead to an improvement in terms of the global compliance
lower than the presented one.

(a) $c_{p}=10^{-11}$

(b) $c_{p}=10^{-12}$

Figure 14: Monoclinic-lattice cantilever beam for different regularization
values. Local Bravais bases are indicated by arrows, and the density is shown
by the color code. The arrows are scaled and shown at selected locations for
readability. Resolution: 400x200 elements. Full design variable fields can be
found in Figure 21.

#### Discrete

The projections of the design variable fields in Figure 21 onto different
finite scales are shown in Figure 13. The resulting lattice has a slight
asymmetry, due to the asymmetric effective properties of monoclinic RUCs (see
Figure 3) and shows different underlying RUC topologies.

The additional design parameter $\beta$ introduces asymmetry in the RUC
anisotropy. This anisotropy reflects in non-rectangular RUCs in the discrete
projection and different orientations of the RUCs in comparison to the
tetragonal and orthorhombic lattices. The monoclinic lattice family hence
shows the most complex lattice architecture yields higher compliance values
than previous results.

##  5 Conclusion

Cellular materials show enormous promise as light-weight materials spanning a
wide mechanical property space. While the large design freedom allows for
target-specific structural integration of spatially variant microstructures,
it also poses a challenging inverse-design problem with respect to
computational costs and manufacturability.

We presented a two-scale topology optimization scheme for complex macroscale
boundary value problems, which guarantees compatible and manufacturable
networks of spatially-variant cellular networks. The approach leverages
crystallographic principles for continuous and generic unit cell definitions
in arbitrary dimensions (here illustrated, without loss of generality, in 2D).
The approach is suitable for a wide range of microscale constitutive
descriptions, including, e.g., analytical, numerical, or ML-based constitutive
laws as well as tabulated material data. We here explored an FE2-type
numerical homogenization approach, which derives the effective constitutive
response from a discrete structural unit cell on the microscale.

We demonstrated the framework by optimizing spatially-variant trusses in 2D
for minimum compliance. A semi-analytical homogenized continuum model allowed
for computationally efficient structural optimization on a regular laptop,
while providing architectures of high resolution. The assumptions intrinsic in
the underlying slender-beam theory put constraints on the beam slenderness and
hence on the volume density of the RUC. While the Timoshenko-Ehrenfest beam
theory has proven to accurately represent beam behavior at low slenderness,
the nodal joints between beams introduce errors. As shown in the Appendix, we
investigated the expected resulting error and expect maximum errors of ca.
$10\%$ for local relative density values up to $50\%$. Even with the imposed
constraints, the minimized compliance was within $6\%$ of the comparable
(unconstrained) benchmark in [37].

In general, we observed a competition between the regularization of the
continuous design variable fields and the length scale of the finite-scale
discrete truss in line with the assumption of a separation of scales. For low
regularization values, where usually best optimization results are obtained, a
very fine length scale is required to ensure similar performance of the
discrete lattice versus the continuum prediction. Conversely, if a coarse
length scale is required (e.g., by manufacturing constraints), stronger
regularization is needed to ensure consistency between the final discrete
structure and the continuum approximation.

While the large design freedom in cellular materials offers an enormous
mechanical property space, its nonconvex nature requires special attention in
topology optimization. Especially, when utilizing a multi-variable
parametrization, the optimization is prone to identify only local minima based
on the initial guess, which can be far off from the global minimum. To avoid
unfavorable local minima, it can be beneficial to explore the design space by
optimizing distinct sets of design variables separately before optimizing the
whole set. In the case of the monoclinic truss lattice family, we observed
that optimization of the (more restricted) orthorhombic lattice family
produced better results. We could only reproduce these results with the
monoclinic lattice when using the output of the orthorhombic lattice
optimization as the initial guess.

Future possible directions include the exploration of different unit cell
definitions in 2D and 3D and their use for different target objectives. A
natural extension with the utilized nonlinear homogenization framework is the
exploration of finite deformations, e.g., for soft robotic applications.

O.S. was supported by the Villum Foundation through the Villum Investigator
Project “InnoTop”.

## References

  * [1] Gibson, L. J., and Ashby, M. F., 1997,  Cellular Solids: Structure and Properties, 2 ed.  Cambridge Solid State Science Series. Cambridge University Press.
  * [2] Ashby, M. F., 2016,  Materials Selection in Mechanical Design Butterworth-Heinemann, Sept.  Google-Books-ID: K4h4CgAAQBAJ.
  * [3] Evans, A. G., Hutchinson, J. W., Fleck, N. A., Ashby, M. F., and Wadley, H. N. G., 2001,  “The topological design of multifunctional cellular metals,”  Progress in Materials Science, 46(3), Jan., pp. 309–327.
  * [4] Lakes, R., 1993,  “Materials with structural hierarchy,”  Nature, 361(6412), Feb., pp. 511–515.
  * [5] Zheng, X., Lee, H., Weisgraber, T. H., Shusteff, M., DeOtte, J., Duoss, E. B., Kuntz, J. D., Biener, M. M., Ge, Q., Jackson, J. A., Kucheyev, S. O., Fang, N. X., and Spadaccini, C. M., 2014,  “Ultralight, ultrastiff mechanical metamaterials,”  Science, 344(6190), June, pp. 1373–1377.
  * [6] Meza, L. R., Zelhofer, A. J., Clarke, N., Mateos, A. J., Kochmann, D. M., and Greer, J. R., 2015,  “Resilient 3d hierarchical architected metamaterials,”  Proceedings of the National Academy of Sciences, 112(37), Sept., pp. 11502–11507.
  * [7] Fleck, N. A., Deshpande, V. S., and Ashby, M. F., 2010,  “Micro-architectured materials: past, present and future,”  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 466(2121), Sept., pp. 2495–2516.
  * [8] Schaedler, T. A., and Carter, W. B., 2016,  “Architected Cellular Materials,”  Annual Review of Materials Research, 46(1), July, pp. 187–210.
  * [9] Valdevit, L., Bertoldi, K., Guest, J., and Spadaccini, C., 2018,  “Architected materials: Synthesis, characterization, modeling, and optimal design,”  Journal of Materials Research, 33(3), Feb., pp. 241–246.
  * [10] Schaedler, T. A., Jacobsen, A. J., Torrents, A., Sorensen, A. E., Lian, J., Greer, J. R., Valdevit, L., and Carter, W. B., 2011,  “Ultralight Metallic Microlattices,”  Science, 334(6058), Nov., pp. 962–965.
  * [11] Lakes, R., 1987,  “Foam Structures with a Negative Poisson’s Ratio,”  Science (New York, N.Y.), 235(4792), Feb., pp. 1038–1040.
  * [12] Sigmund, O., 1994,  “Materials with prescribed constitutive parameters: An inverse homogenization problem,”  International Journal of Solids and Structures, 31(17), Sept., pp. 2313–2329.
  * [13] Milton, G. W., and Cherkaev, A. V., 1995,  “Which Elasticity Tensors are Realizable?,”  Journal of Engineering Materials and Technology, 117(4), 10, pp. 483–493.
  * [14] Huang, X., Radman, A., and Xie, Y. M., 2011,  “Topological design of microstructures of cellular materials for maximum bulk or shear modulus,”  Computational Materials Science, 50(6), Apr., pp. 1861–1870.
  * [15] Bückmann, T., Thiel, M., Kadic, M., Schittny, R., and Wegener, M., 2014,  “An elasto-mechanical unfeelability cloak made of pentamode metamaterials,”  Nature Communications, 5(1), p. 4130.
  * [16] Wadley, H., Dharmasena, K., Chen, Y., Dudt, P., Knight, D., Charette, R., and Kiddy, K., 2008,  “Compressive response of multilayered pyramidal lattices during underwater shock loading,”  International Journal of Impact Engineering, 35(9), Sept., pp. 1102–1114.
  * [17] Hussein, M. I., Leamy, M. J., and Ruzzene, M., 2014,  “Dynamics of Phononic Materials and Structures: Historical Origins, Recent Progress, and Future Outlook,”  Applied Mechanics Reviews, 66(4), 05  040802\.
  * [18] Zelhofer, A. J., and Kochmann, D. M., 2017,  “On acoustic wave beaming in two-dimensional structural lattices,”  International Journal of Solids and Structures, 115-116, pp. 248–269.
  * [19] Clausen, A., Wang, F., Jensen, J. S., Sigmund, O., and Lewis, J. A., 2015,  “Topology Optimized Architectures with Programmable Poisson’s Ratio over Large Deformations,”  Advanced Materials, 27(37), pp. 5523–5527  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.201502485.
  * [20] Symons, D. D., and Fleck, N. A., 2008,  “The Imperfection Sensitivity of Isotropic Two-Dimensional Elastic Lattices,”  Journal of Applied Mechanics, 75(5), Sept.
  * [21] Montemayor, L. C., Wong, W. H., Zhang, Y.-W., and Greer, J. R., 2016,  “Insensitivity to Flaws Leads to Damage Tolerance in Brittle Architected Meta-Materials,”  Scientific Reports, 6, Feb., p. 20570.
  * [22] Thomsen, C. R., Wang, F., and Sigmund, O., 2018,  “Buckling strength topology optimization of 2D periodic materials based on linearized bifurcation analysis,”  Computer Methods in Applied Mechanics and Engineering, 339, Sept., pp. 115–136.
  * [23] Pham, M.-S., Liu, C., Todd, I., and Lertthanasarn, J., 2019,  “Damage-tolerant architected materials inspired by crystal microstructure,”  Nature, 565(7739), Jan., pp. 305–311.
  * [24] Cadman, J. E., Zhou, S., Chen, Y., and Li, Q., 2013,  “On design of multi-functional microstructural materials,”  Journal of Materials Science, 48(1), Jan., pp. 51–66.
  * [25] Osanov, M., and Guest, J. K., 2016,  “Topology Optimization for Architected Materials Design,”  Annual Review of Materials Research, 46(1), pp. 211–233.
  * [26] Jain, A., Ong, S. P., Hautier, G., Chen, W., Richards, W. D., Dacek, S., Cholia, S., Gunter, D., Skinner, D., Ceder, G., and Persson, K. A., 2013,  “Commentary: The Materials Project: A materials genome approach to accelerating materials innovation,”  APL Materials, 1(1), July, p. 011002.
  * [27] Sigmund, O., 2000,  “A new class of extremal composites,”  Journal of the Mechanics and Physics of Solids, 48(2), Feb., pp. 397–428.
  * [28] Berger, J. B., Wadley, H. N. G., and McMeeking, R. M., 2017,  “Mechanical metamaterials at the theoretical limit of isotropic elastic stiffness,”  Nature, 543, Feb, pp. 533 EP –.
  * [29] Sigmund, O., and Jensen, J. S., 2003,  “Systematic design of phononic band–gap materials and structures by topology optimization,”  Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 361(1806), May, pp. 1001–1019.
  * [30] Bilal, O. R., and Hussein, M. I., 2011,  “Ultrawide phononic band gap for combined in-plane and out-of-plane waves,”  Physical Review E, 84(6), Dec.  arXiv: 1111.1457.
  * [31] Pasini, D., and Guest, J. K., 2019,  “Imperfect architected materials: Mechanics and topology optimization,”  MRS Bulletin, 44(10), p. 766–772.
  * [32] Kochmann, D. M., Hopkins, J. B., and Valdevit, L., 2019,  “Multiscale modeling and optimization of the mechanics of hierarchical metamaterials,”  MRS Bulletin, 44(10), p. 773–781.
  * [33] Aage, N., Andreassen, E., Lazarov, B. S., and Sigmund, O., 2017,  “Giga-voxel computational morphogenesis for structural design,”  Nature, 550(7674), Oct., pp. 84–86.
  * [34] Lazarov, B., Wang, F., and Sigmund, O., 2016,  “Length scale and manufacturability in density-based topology optimization,”  Archives of Applied Mechanics, 86, pp. 189–218.
  * [35] Milton, G. W., 1986,  “Modelling the Properties of Composites by Laminates,”  In Homogenization and Effective Moduli of Materials and Media, J. L. Ericksen, D. Kinderlehrer, R. Kohn, and J.-L. Lions, eds., The IMA Volumes in Mathematics and its Applications, Springer New York, pp. 150–174.
  * [36] Träff, E., Sigmund, O., and Groen, J., 2018,  “Simple single-scale microstructures based on optimal rank-3 laminates,”  Structural and Multidisciplinary Optimization, 59, pp. 1021–1031.
  * [37] Groen, J. P., and Sigmund, O., 2018,  “Homogenization-based topology optimization for high-resolution manufacturable microstructures,”  International Journal for Numerical Methods in Engineering, 113(8), pp. 1148–1163.
  * [38] Rodrigues, H., Guedes, J., and Bendsoe, M., 2002,  “Hierarchical optimization of material and structure,”  Structural and Multidisciplinary Optimization, 24(1), Aug., pp. 1–10.
  * [39] Coelho, P. G., Fernandes, P. R., Guedes, J. M., and Rodrigues, H. C., 2008,  “A hierarchical model for concurrent material and topology optimisation of three-dimensional structures,”  Structural and Multidisciplinary Optimization, 35(2), Feb., pp. 107–115.
  * [40] Sivapuram, R., Dunning, P. D., and Kim, H. A., 2016,  “Simultaneous material and structural optimization by multiscale topology optimization,”  Structural and Multidisciplinary Optimization, 54(5), Nov., pp. 1267–1281.
  * [41] Wang, Y., Xu, H., and Pasini, D., 2017,  “Multiscale isogeometric topology optimization for lattice materials,”  Computer Methods in Applied Mechanics and Engineering, 316, Apr., pp. 568–585.
  * [42] Xia, L., and Breitkopf, P., 2014,  “Concurrent topology optimization design of material and structure within FE2 nonlinear multiscale analysis framework,”  Computer Methods in Applied Mechanics and Engineering, 278, Aug., pp. 524–542.
  * [43] Khanoki, S. A., and Pasini, D., 2012,  “Multiscale design and multiobjective optimization of orthopedic hip implants with functionally graded cellular material.,”  Journal of biomechanical engineering, 134(3), p. 031004.
  * [44] Wang, Y., Arabnejad, S., Tanzer, M., and Pasini, D., 2018,  “Hip Implant Design With Three-Dimensional Porous Architecture of Optimized Graded Density,”  Journal of Mechanical Design, 140(11), Nov.
  * [45] Schumacher, C., Bickel, B., Rys, J., Marschner, S., Daraio, C., and Gross, M., 2015,  “Microstructures to control elasticity in 3d printing,”  ACM Transactions on Graphics, 34(4), July, pp. 136:1–136:13.
  * [46] Chen, D., Skouras, M., Zhu, B., and Matusik, W., 2018,  “Computational discovery of extremal microstructure families,”  Science Advances, 4(1), Jan., p. eaao7005.
  * [47] Radman, A., Huang, X., and Xie, Y. M., 2013,  “Topology optimization of functionally graded cellular materials,”  Journal of Materials Science, 48(4), Feb., pp. 1503–1510.
  * [48] Zhu, B., Skouras, M., Chen, D., and Matusik, W., 2017,  “Two-Scale Topology Optimization with Microstructures,”  arXiv:1706.03189 [cs], June  arXiv: 1706.03189.
  * [49] Du, Z., Zhou, X.-Y., Picelli, R., and Kim, H. A., 2018,  “Connecting Microstructures for Multiscale Topology Optimization With Connectivity Index Constraints,”  Journal of Mechanical Design, 140(11), Nov.
  * [50] Panetta, J., Zhou, Q., Malomo, L., Pietroni, N., Cignoni, P., and Zorin, D., 2015,  “Elastic textures for additive fabrication,”  ACM Transactions on Graphics, 34(4), July, pp. 135:1–135:12.
  * [51] Schury, F., Stingl, M., and Wein, F., 2012,  “Efficient Two-Scale Optimization of Manufacturable Graded Structures,”  SIAM Journal on Scientific Computing, 34(6), Jan., pp. B711–B733.
  * [52] Sanders, E. D., Aguiló, M. A., and Paulino, G. H., 2018,  “Multi-material continuum topology optimization with arbitrary volume and mass constraints,”  Computer Methods in Applied Mechanics and Engineering, 340, Oct., pp. 798–823.
  * [53] Zhang, X. S., Chi, H., and Paulino, G. H., 2020,  “Adaptive multi-material topology optimization with hyperelastic materials under large deformations: A virtual element approach,”  Computer Methods in Applied Mechanics and Engineering, 370, Oct., p. 112976.
  * [54] Sanders, E. D., Pereira, A., and Paulino, G. H., 2021,  “Optimal and continuous multilattice embedding,”  Science Advances, Apr.  Publisher: American Association for the Advancement of Science.
  * [55] Dunning, P. D., Brampton, C. J., and Kim, H. A., 2015,  “Simultaneous optimisation of structural topology and material grading using level set method,”  Materials Science and Technology, 31(8), June, pp. 884–894.
  * [56] Greifenstein, J., and Stingl, M., 2016,  “Simultaneous parametric material and topology optimization with constrained material grading,”  Structural and Multidisciplinary Optimization, 54(4), Oct., pp. 985–998.
  * [57] Allaire, G., Geoffroy-Donders, P., and Pantz, O., 2019,  “Topology optimization of modulated and oriented periodic microstructures by the homogenization method,”  Computers & Mathematics with Applications, 78(7), Oct., pp. 2197–2229.
  * [58] Xu, L., and Qian, Z., 2021,  “Topology optimization and de-homogenization of graded lattice structures based on asymptotic homogenization,”  Composite Structures, 277, Dec., p. 114633.
  * [59] Kumar, T., and Suresh, K., 2020,  “A density-and-strain-based K-clustering approach to microstructural topology optimization,”  Structural and Multidisciplinary Optimization, 61(4), Apr., pp. 1399–1415.
  * [60] Pantz, O., and Trabelsi, K., 2008,  “A Post-Treatment of the Homogenization Method for Shape Optimization,”  SIAM Journal on Control and Optimization, Mar.
  * [61] Wu, J., Wang, W., and Gao, X., 2019,  “Design and Optimization of Conforming Lattice Structures,”  IEEE Transactions on Visualization and Computer Graphics, pp. 1–1  arXiv: 1905.02902.
  * [62] Wu, J., Sigmund, O., and Groen, J. P., 2021,  “Topology optimization of multi-scale structures: a review,”  Structural and Multidisciplinary Optimization, 63(3), Mar., pp. 1455–1480.
  * [63] Glaesener, R. N., Lestringant, C., Telgen, B., and Kochmann, D. M., 2019,  “Continuum models for stretching- and bending-dominated periodic trusses undergoing finite deformations,”  International Journal of Solids and Structures, 171, Oct., pp. 117–134.
  * [64] Rumpf, R. C., and Pazos, J., 2012,  “Synthesis of spatially variant lattices,”  Optics Express, 20(14), July, pp. 15263–15274.
  * [65] Groen, J., Stutz, F., Aage, N., Bærentzen, J., and Sigmund, O., 2020,  “De-homogenization of optimal multi-scale 3d topologies,”  Computer Methods in Applied Mechanics and Engineering, 364, p. 112979.
  * [66] Xue, D., Zhu, Y., Li, S., Liu, C., Zhang, W., and Guo, X., 2020,  “On speeding up an asymptotic-analysis-based homogenisation scheme for designing gradient porous structured materials using a zoning strategy,”  Structural and Multidisciplinary Optimization, pp. 1–17.
  * [67] Geoffroy-Donders, P., Allaire, G., and Pantz, O., 2020,  “3-d topology optimization of modulated and oriented periodic microstructures by the homogenization method,”  Journal of Computational Physics, 401, p. 108994.
  * [68] Kittel, C., 1986,  Introduction to solid state physics, 8th ed.  Wiley.
  * [69] Bravais, A., 1850,  “Mémoire sur les systèmes formés par des points distribués régulièrement sur un plan ou dans l’espace,”  J. École polytechnique, 19, pp. 1––128.
  * [70] Wigner, E., and Seitz, F., 1933,  “On the Constitution of Metallic Sodium,”  Physical Review, 43(10), May, pp. 804–810.
  * [71] Crisfield, M. A., 1990,  “A consistent co-rotational formulation for non-linear, three-dimensional, beam-elements,”  Computer methods in applied mechanics and engineering, 81(2), pp. 131–150.
  * [72] Borrvall, T., 2001,  “Topology optimization of elastic continua using restriction,”  Archives of Computational Methods in Engineering, 8(4), Dec., pp. 351–385.
  * [73] Svanberg, K., 1987,  “The method of moving asymptotes—a new method for structural optimization,”  International Journal for Numerical Methods in Engineering, 24(2), pp. 359–373.
  * [74] Sigmund, O., Aage, N., and Andreassen, E., 2016,  “On the (non-)optimality of Michell structures,”  Structural and Multidisciplinary Optimization, 54(2), Aug., pp. 361–373.
  * [75] Sigmund, O., and Petersson, J., 1998,  “Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, mesh-dependencies and local minima,”  Structural Optimization, 16(1), pp. 68–75.
  * [76] Chan, A. S. L., 1962,  “The Design of Michell Optimum Structures,”  p. 43.
  * [77] Pedersen, P., 1989,  “On optimal orientation of orthotropic materials,”  Structural optimization, 1(2), June, pp. 101–106.
  * [78] Groen, J. P., Wu, J., and Sigmund, O., 2019,  “Homogenization-based stiffness optimization and projection of 2d coated structures with orthotropic infill,”  Computer Methods in Applied Mechanics and Engineering, 349, June, pp. 722–742.

##  Appendix A Homogenization based on beams

The utilized homogenization approach of [63] is based on modeling the RUC by
beam elements, which imply rigid (welded) connections at junctions. Beam
theory is generally deemed applicable whenever the length of a beam is much
larger than its width, i.e., when the slenderness $L/w\gg 1$. In addition,
assuming rigid connections between beams neglects nodal effects such as stress
concentrations at nodes or overlaps between connected beams. These
assumptions, however, are challenged when modeling RUCs with high volume
fractions $\bar{\rho}$, implying relatively thick beams. To evaluate the
applicability of the applied homogenization approach, we present in this
section a discussion of the influence of the effective RUC density on the
approximation of homogenized quantities and ways to extend the applicable
range. For all calculations we assume a rectangular beam cross-section with
out-of-plane thickness $t=1$, Young’s modulus $E=1$, Poisson’s ration
$\nu=0.3$, and the Timoshenko shear coefficient for rectangular cross-sections
$\kappa=1.2$ (all normalized for simplicity).

###  A.1 Beam slenderness

As the beam is the smallest building block of the RUC, the exactness of the
beam element has primary relevance. Generally well-known beam formulations
include the theories of Euler-Bernoulli and Timoshenko-Ehrenfest. To evaluate
their accuracy for low slenderness ratios $L/w$, we evaluate, using both
formulations, the strain energy of a 2D cantilever beam of length $L=1$ and
varying thickness $w$. The beam is clamped at one end and loaded by a constant
displacement at the other end. For comparison, we establish a reference by
resolving the cantilever beam ($20$ elements over the width $w$) with square
bilinear finite elements, assuming isotropic linear elasticity under the
plane-stress assumption.

Figure 15: Total strain energy of a cantilever beam (clamped at one end,
loaded by a constant deflection at the other end) for varying slenderness
$w/L$. The fully resolved result is shown in black (solid line), compared to
results obtained from Timoshenko-Ehrenfest beam theory (blue stars) and Euler-
Bernoulli beam theory (red crosses).

Figure 16: Different node treatments in a square RUC ($\bar{\rho}=0.5$):
variant A ignores the overlap of beams at a node (dark areas) and uses a beam
element from node center to node center. In variant B the nodes are assumed
rigid and beams have the effective length between nodes without overlap.
Variant C is a volume-preserving compromise between variants A and B.

A comparison of the strain energy inside the bent cantilever beam is shown in
Figure 16. Both beam implementations converge to the exact solution for
$L/w\gg 1$. With increasing $w/L$, Euler-Bernoulli beam theory exceedingly
overestimates the strain energy. In contrast, Timoshenko-Ehrenfest beam theory
shows significantly better agreement and generally tends to underestimate the
strain energy. For $w/L=0.5$, Euler-Bernoulli beam theory overstimates the
energy by roughly $27\%$, while Timoshenko-Ehrenfest stays within a $5\%$
error margin for the full range of slenderness ratios up to $w=L$.

We conclude that Euler-Bernoulli beam theory is not suited for modeling beams
with low slenderness ratio $L/w$. In contrast, when accepting a maximum error
of $5\%$, Timoshenko-Ehrenfest beam theory is applicable for a wide range of
slenderness ratios and even for extreme cases ($w=L$). Therefore, all
following discussions (and the truss optimization studies) are based on
Timoshenko-Ehrenfest beam theory.

###  A.2 Node influence

Truss networks consist of many beams that are connected to each other through
nodes. While the previous section confirmed the validity of beam theory for an
individual beam, we here investigate the influence of nodes on the
approximation of homogenized properties using beam elements. Generally, the
influence of nodes is neglected when using beam models, assuming a slenderness
ratio of $L/w\gg 1$. However, when beams are not slender, or equivalently when
the density $\bar{\rho}$ of the RUC is high, the influence of the nodes on the
effective mechanical properties cannot be neglected.

In Figure 16, we introduce three variants for approximating the geometry of a
RUC with non-slender beams. As in classical beam theory, all variants assume
that the node is rigid and hence described well by Timoshenko-Ehrenfest
theory. In variant A, beams connect from node centroid to node centroid; beams
overlap at nodes (this is the standard approach). In variant B, the
overlapping nodes are accounted for by assuming beams whose effective length
is not from node centroid to node centroid but the reduced length between
nodes without overlapping volumes. Variants A and B present the soft and stiff
extremes (lower and upper bounds), respectively. Finally, variant C is a
compromise, where the overlap of beams at each node is chosen such that the
volume of all beams in a RUC matches the volume fill fraction according to the
effective density $\bar{\rho}$, which defines the effective beam length
(between those of variants A and B). Note that, when $L/w\gg 1$, all three
variants must converge to the same solution.

(a) square RUC, $\bar{\rho}=0.1$

(b) square RUC, $\bar{\rho}=0.5$

(c) hexagon RUC, $\bar{\rho}=0.1$

(d) hexagon RUC, $\bar{\rho}=0.5$

(e) square RUC, $\bar{\rho}=0.1$

(f) square RUC, $\bar{\rho}=0.5$

(g) hexagonal RUC, $\bar{\rho}=0.1$

(h) hexagonal RUC, $\bar{\rho}=0.5$

Figure 17: Influence of the node treatment on the effective directional
Young’s modulus for different effective RUC densities and topologies. Figures
17(a) to 17(d) show the CST mesh of the square RUC ($\alpha=0$, $\beta=\pi/2$,
$\gamma=1$) and hexagonal RUC ($\alpha=0$, $\beta=\pi/3$, $\gamma=1$), each
for $\bar{\rho}=0.1$ and $\bar{\rho}=0.5$. Figures 17(e) to 17(h) show the
respective directional Young’s moduli based on homogenization of the shown
RUCs (black solid lines) and obtained from the modified beam representation in
variant A, ignoring the overlap at nodes (green dashed lines), assuming rigid
nodes in variant B (blue dotted lines), and volume-preserving node overlaps in
variant C (red dash-dotted lines). Expoiting symmetry, only half of the first
quadrant is shown in each elastic modulus plot.

The influence of the three variants on the effective stiffness is demonstrated
for the square and hexagon RUCs, as representatives of stretching- and
bending-dominated RUCs, in Figure 17. As a reference, we calculated the
effective directional Young’s modulus of fully discretized (ref) RUCs, also
shown in Figures 17(a) to 17(d), with classical homogenization assuming
periodic boundary conditions and linear elasticity under the plane-stress
assumption. These reference values are indicated as black solid lines. The
directional Young’s modulus based on truss homogenization [63] for the three
variants A, B, C are shown in green, blue, and red, respectively.

The results confirm that for low effective density (e.g., $\bar{\rho}=0.1$)
all three variants converge to the fully-resolved reference solution. Variant
A converges from below, variant B converges from above, and there is no clear
trend for variant C. For the large density $\bar{\rho}=0.5$, the trends are
the same, but the deviations are more severe. Variant A underestimates the
directional stiffness of the square RUC by maximally $5\%$, variants B and C
overestimate the directional stiffness by up to $34\%$ and $11\%$,
respectively. For the hexagonal RUC, variant B still overestimates the
directional stiffness by $16\%$, and variants A and C underestimate the
stiffness by $23\%$ and $7\%$, respectively.

The comparison is extended to the full homogenized orthotropic elasticity
tensor $\bm{C}$ of the RUC (having $3\times 3$ components, using Voigt
notation in 2D) . Results are shown in Figure 18 with the fully discretized
reference solution indicated by $C^{\text{ref}}$. As expected, the errors
shown in 18(a) and 18(d) converge to zero in the limit of low densities
$\bar{\rho}$. For $\bar{\rho}=0.5$ the error of variant B ($>30\%$) is highest
in the square, and variant A shows the largest error in the hexagon ($>20\%$).
For the $C_{11}$-entry, 18(b) and 18(e)) confirm the observations from the
directional Young’s moduli (cf. Figure 17): with decreasing $\bar{\rho}$,
variant A converges from below, variant B from above, and variant C is
scattered around the reference values. Lastly, we compare a measure for the
shear stiffness in form of invariant $(C_{11}+C_{22})/2-(C_{12}+2C_{66})$, as
presented by [77]. When the invariant assumes a negative value, the shear
stiffness is regarded as relatively high and could indicate artificial shear
stiffening. By contrast, a positive value indicates a relatively low shear
stiffness. In 18(f), we observe a constant zero invariant for the hexagon,
whereas the square RUC loses relative shear stiffness with increasing
$\bar{\rho}$, as shown in 18(c).

###  A.3 Applicability of beam-based homogenization

We showed in A.1 that the Timoshenko-Ehrenfest beam theory provides accurate
results up to low slenderness values (thick beams), so that we deem this model
applicable for the homogenization of RUCs with relatively high filling
fraction $\bar{\rho}$. In A.2, we observed that the nodes in a RUC can be
assumed rigid for relatively low $\bar{\rho}$. However, the error increases
drastically for increasing $\bar{\rho}$. Clear trends are observable: Variant
A consistently underestimates the reference, while variant B consistently
overestimates the effective properties. Variant C lies in between those two
bounds. No artificial effects, such as shear stiffening, could be detected.

In conclusion, one must carefully evaluate the acceptable error when choosing
to homogenize RUCs with high $\bar{\rho}$ while assuming rigid nodes (as done
in this study based on Timoshenko-Ehrenfest beam theory). For this
optimization study, we chose to use variant C as a compromise and accept
maximum local errors on the order of $~{}10\%$ by allowing for local relative
densities $\rho(\bm{X})$ up to $50\%$.

(a) square RUC

(b) square RUC

(c) square RUC

(d) hexagon RUC

(e) hexagon RUC

(f) hexagon RUC

Figure 18: Influence of effective density and node treatment on estimating the
elasticity tensor, comparing the results of beam elements of variants A, B,
and C to the fully-meshed reference solution (using Voigt notation in 2D for
the stiffness tensor components).

##  Appendix B Continuum design variable fields

All studied design variable fields are shown in Figures 19, 20 and 21.

1

Figure 19: Optimized design variable fields for the tetragonal-lattice
cantilever problem, using 400x200 elements. Top: $c_{p}=10^{-13}$; middle:
$10^{-12}$; bottom: $10^{-11}$. Figure 20: Optimized design variable fields
for the orthorhombic-lattice cantilever problem, using 400x200 elements. From
top to bottom, $c_{p}=(10^{-13},10^{-12},10^{-11})$. Figure 21: Optimized
design variable fields for the monoclinic-lattice cantilever problem, using
400x200 elements and the optimized fields of the for the orthorhombic lattice
family as starting guess. From top to bottom, $c_{p}=(10^{-12},10^{-11})$.
Figure 22: Optimized design variable fields for the monoclinic-lattice
cantilever problem, using 400x200 elements and the homogeneous starting guess.
From top to bottom, $c_{p}=(10^{-12},10^{-11})$.