WGMwgm QEqe EPep PMSpms BECbec DEde

[ orcid=0000-0002-4607-844X]

Conceptualization, Methodology, Software, Writing - Original Draft

1]organization=Technical University of Munich, addressline=Boltzmannstr. 3,
city=Garching, postcode=85748, country=Germany

[ orcid=0000-0002-3394-0731] Conceptualization, Methodology, Writing - Review
& Editing, Supervision, Funding acquisition

[ orcid=0000-0003-0829-7075] Methodology, Writing - Review & Editing

2]organization=Lightspeed Studios, city=Bellevue, state=WA, country=USA

[orcid=0000-0003-4237-1806] [1]

Conceptualization, Methodology, Writing - Review & Editing, Supervision

3]organization=Delft University of Technology, addressline=Landbergstraat 15,
city=Delft, postcode=2628 CE, country=The Netherlands

[cor1]Corresponding author

# Design and Optimization of Functionally-graded Triangular Lattices for
Multiple Loading Conditions

Junpeng Wang<EMAIL_ADDRESS>[  Rüdiger Westermann<EMAIL_ADDRESS>Xifeng Gao<EMAIL_ADDRESS>[  Jun Wu<EMAIL_ADDRESS>[

###### Abstract

Aligning lattices based on local stress distribution is crucial for achieving
exceptional structural stiffness. However, this aspect has primarily been
investigated under a single load condition, where stress in 2D can be
described by two orthogonal principal stress directions. In this paper, we
introduce a novel approach for designing and optimizing triangular lattice
structures to accommodate multiple loading conditions, which means multiple
stress fields. Our method comprises two main steps: homogenization-based
topology optimization and geometry-based de-homogenization. To ensure the
geometric regularity of triangular lattices, we propose a simplified version
of the general rank-$3$ laminate and parameterize the design domain using
equilateral triangles with unique thickness per edge. During optimization, the
thicknesses and orientation of each equilateral triangle are adjusted based on
the homogenized properties of triangular lattices. Our numerical findings
demonstrate that this proposed simplification results in only a slight
decrease in stiffness of less than $5\%$ compared to using the general rank-3
laminate, while achieving triangular lattice structures with a compelling
geometric regularity. In geometry-based de-homogenization, we adopt a field-
aligned triangulation approach to generate a globally consistent triangle
mesh, with each triangle oriented according to the optimized orientation
field. Our approach for handling multiple loading conditions, akin to de-
homogenization techniques for single loading conditions, yields highly
detailed, optimized, spatially varying lattice structures. The method is
computationally efficient, as simulations and optimizations are conducted at a
low-resolution discretization of the design domain. Furthermore, since our
approach is geometry-based, obtained structures are encoded into a compact
geometric format that facilitates downstream operations such as editing and
fabrication.

###### keywords:

Topology optimization De-homogenization Lattice structures Multiple loading
conditions

##  1 Introduction

Topology optimization and lattice infill are primary strategies for designing
lightweight structures for additive manufacturing. The integration of these
two strategies has been a topic of intensive research in the past decade [20].
Topology optimization transforms structural design into an optimization
problem of determining where to place material, while lattice infill can be
considered as a form of metamaterial. Topology optimization with lattice
infill results in multi-scale structures that are lightweight and robust [5,
34]. For an overview of topology optimization methods for designing multi-
scale structures, please refer to the recent review article by Wu et al. [36].

For designing load-bearing structures, tailoring the layout of lattice infill
based on the stress distribution is key to achieving exceptional structural
performance. Lattice infill can be categorized as isotropic or anisotropic. By
aligning the anisotropy of the lattice with the principal stress direction at
each location, the structural rigidity can be significantly increased [22].
This principle has been well explored and several approaches have been
developed [18, 6]. Among them, an effective and efficient approach is de-
homogenization [21, 9, 37]. It involves two major steps. The first step,
homogenization-based topology optimization, specifies the optimal lattice
infill configuration at each location in the design domain [3]. The second
step, de-homogenization, is crucial in this approach and it transforms the
optimized specifications into a globally consistent structure.

De-homogenization approaches can be classified, depending on the
representation of optimized designs, into image-based approaches and geometry-
based approaches. Image-based approaches represent the structural layout using
binary density fields, while geometry-based approaches represent it using a
compact geometric format. For example, a truss-like structure can be
represented by a mesh where each edge in the mesh represents a truss.
Geometry-based approaches offer direct control of geometric features, which is
beneficial in, e.g., addressing manufacturing constraints [14]. The geometric
format also offers a compact parameter space for subsequent editing operations
on the optimized structures [37, 30]. From a computational perspective, image-
based approaches center on computing a fine-grid scalar field, whose gradients
are aligned with the orientation of optimized lattice infill from
homogenization-based topology optimization, e.g., the projection-based methods
[21]. Geometric approaches, in contrast, involve directly computing a set of
field-aligned geometric primitives to represent the optimized lattice infill,
e.g., the conforming lattice structure [37]. De-homogenization for a single
loading condition is already a prosperous realm. The projection-based approach
proposed by Pantz and Trabelsi [21] was simplified by Groen and Sigmund [9]
and further developed by Allaire _et al._ [1]. It has since been extended to
design shell-infill structures [11], and to work in 3D [10, 8]. Different
approaches for the projection have been proposed, from related research fields
such as visualization and computer graphics, using stream surfaces [25],
convolutional neural networks [7], and phasor noises [33]. Geometric
approaches include those based on field-aligned hex-dominant meshing [37] and
streamline tracing [30]. The singularity issue in 2D has been investigated
respectively in the image- and geometry-based routines by Stutz _et al._ [24]
and Wang _et al._ [30].

The development of de-homogenization approaches considers almost exclusively a
single loading condition, i.e., constant external loads. In many engineering
applications, however, multiple loading conditions are omnipresent. For
instance, in airplane design, the major loads acting on the aircraft wing vary
at different stages during the flight, including lift and drag, gravitational
forces (particularly when the fuel tanks are fully loaded), propulsion forces
from the forward engines, and opposing thrust from thrust reversers. All these
loads are exerted on the airplane during the flight, thus requiring to
incorporate them into the structural design and optimization process.

De-homogenization under multiple loading conditions, however, is challenging.
Firstly, the optimized metamaterial properties can be realized by different
lattice configurations, i.e., the optimal lattice configurations are not
uniquely defined. Furthermore, the optimized lattice configuration exhibits
severe discontinuities. This further complicates the de-homogenization of a
consistent structural layout. A first solution to these challenges was
recently proposed by Jensen _et al._ [16]. They developed an image-based
method, extending the projection approach in [9] from single to multiple
loading conditions. The rank-3 laminate was used as the lattice infill for
homogenization. It describes a composite material composed of three
differently oriented material layers, where each material layer is described
by an orientation variable and a width variable. I.e., a rank-3 laminate
involves three independent width variables and three orientation variables.
The structural designs by this approach achieve comparable mechanical
performance (i.e., stiffness) with the ones obtained from conventional
density-based topology optimization, while the latter incurs a significantly
larger computational burden. The image-based projection method inherently
requires the integrability of the optimized specifications, making it
necessary to ensure the smoothness of the optimized orientation field. This,
however, becomes difficult due to the non-uniqueness of the optimal solutions
of the homogenization model. In addition, as the three orientation fields are
de-homogenized independently, the sub-structures aligning with different
orientation fields can come very close to each other, leaving small void
regions in the final designs. While this feature does not necessarily impact
mechanical performance, it may not be the preferred choice in certain design
or manufacturing scenarios due to the increased geometric complexity it
introduces. See the close-up view of fig:Overviewb for an illustration of
irregular intersections of three orientation fields.

In this paper, we present a novel method to generate functionally graded
triangular lattice structures that are resilient to multiple loading
conditions. Our method follows the established de-homogenization pipeline,
consisting two steps: (1) performing homogenization-based topology
optimization to obtain a set of optimized specifications; (2) converting the
optimized specifications into a graded triangular lattice structure.
fig:Overview shows an overview of the proposed method.

Figure 1: Method overview. (a) The optimization problem: Two loading
conditions $F_{1}$ (orange) and $F_{2}$ (green) are applied, and the bottom is
fixed. (b) The optimized density layout (background) and orientations
(streamlines) generated by homogenization-based topology optimization. (c) The
triangular parametrization with edges aligned to the tangents of the optimized
rank-3 layers. The entire design domain and the optimized density layout are
partitioned into sub-domains by this field-aligned mesh. (d) The triangular
lattice design obtained via de-homogenization.

Our technical contributions are two-fold. Firstly, we propose a simple and
effective infill model for multiple loading conditions — triangular lattices.
It can be viewed as a simplified version of the general rank-3 laminate. In
triangular lattices, the width of each edge is allowed to change
independently, creating a large range of attainable metamaterial properties.
More importantly, the triangles are rotatable, allowing local alignment of
their anisotropy with preferred directions. Secondly, we develop a geometric
approach to translate the optimized specifications into a globally conforming
triangular lattice structure, whose edges align with the three optimized
orientation fields and whose edge thicknesses convey the mechanical
anisotropy. Our de-homogenization approach is based on field-aligned
triangulation [15]. This method does not require the integrability of the
optimized specifications, as it is based on a local smoothing operator,
avoiding global parameterization. In addition, compared to the randomly
intersected substructures, our de-homogenization using equilateral triangles
restricts the geometric complexity by carefully placing the intersections,
resulting in triangular lattices with equal spacing, as shown in the close-up
views of fig:Overviewc and d.

The remainder of this paper is organized as follows. In Sec:Opti we first
present the triangular lattices by simplifying the general rank-$3$ laminate,
and we introduce the corresponding equation of homogenization-based topology
optimization. The construction of de-homogenization is presented in
Sec:DeHomoDiagram. In Sec:Result, we perform numerical experiments and use
them to discuss and compare the proposed method. Finally, Sec:Conclusion
concludes the paper with a summary and an outlook.

##  2 Topology Optimization with Triangular Lattices

###  2.1 Triangular Lattices

Our triangular infill model is a simplified version of the general rank-$3$
laminate. Stiffness-optimal structures can be described by a set of spatially
varying rank-$N$ laminates with local periodicity [3]. Rank-$N$ describes a
composite material composed of $N$ differently oriented layers. For a single
loading condition, the optimal laminate can be described as two orthogonally
oriented layers, i.e., rank-$2$ laminate, which was simplified as a square
cell with a rectangular hole in homogenization-based topology optimization
[3]. However, under multiple loading conditions, there are more than two major
stress directions in the structure, one thus needs to resort to the
anisotropic rank-$3$ to achieve optimal stiffness. The optimal stiffness under
multiple loading conditions may also be achieved by laminates with a higher
rank, yet this introduces unnecessary computational complexity and numerical
instability during optimization [16].

The rank-$3$ laminate is constructed by piling up three layers of materials
that orient differently, assuming perfect bound condition among different
layers (fig:rank3Diagb). On each layer, two phases of material are considered,
i.e., stiff material ($+$) and compliant material ($-$). The stiff material
indicates an isotropic solid material deposition, while the compliant material
varies among the different layers. A detailed description of the composition
of a rank-3 laminate can be found in existing works [17, 12, 16]. It is
briefly explained here for completeness.

Figure 2: Schematic composition of a rank-3 laminate and its geometric
interpretation. (a) The design domain and boundary conditions for
optimization, where each simulation cell is depicted by a set of specific
specifications of a rank-3 laminate. (b) The multi-scale rank-3 laminate in a
selected simulation cell. The material is composed of three differently
oriented layers, shown in (c), (d), and (e), respectively. $\alpha$ is the
relative width of each layer. $\bm{n}$ and $\bm{t}$ refer to the normal and
tangent direction of a layer, respectively. (f) The proposed single-scale
geometric interpretation of the rank-3 laminate as an equilateral triangle,
with a unique thickness per edge.

Let a macro-cell have unit size, with the base layer being layer $3$, i.e.,
the rank-$3$ laminate is constructed by starting from layer $3$
(fig:rank3Diagc). The solid material deposition ($\rho_{3}$) on layer $3$ is
then described by the width variable $\alpha_{3}$ and the compliant material
$1-\alpha_{3}$. Here the compliant material is rank-$2$. Sequentially on layer
$2$ with width variable $\alpha_{2}$ for the solid material (fig:rank3Diagd),
the remaining region to fill is $1-\alpha_{3}$. The actual consumption of
solid material ($\rho_{2}$) on layer $2$ is $(1-\alpha_{3})\alpha_{2}$.
Correspondingly, the compliant material on layer $2$ is rank-$1$, and its
volume fraction is $(1-\alpha_{3})(1-\alpha_{2})$. On layer $1$ with width
variable $\alpha_{1}$ (fig:rank3Diage), the solid material consumption
($\rho_{1}$) is $(1-\alpha_{3})(1-\alpha_{2})\alpha_{1}$, and the
corresponding compliant material is ‘void’, with a ratio of
$(1-\alpha_{3})(1-\alpha_{2})(1-\alpha_{1})$. The ‘void’ material is assigned
a small Young’s modulus $E^{-}$ to avoid numerical instabilities, yet which is
negligible compared to the Young’s modulus $E^{+}$ of the solid material
(i.e., $E^{-}=10^{-6}E^{+}$). The Poisson’s ratio $v_{0}$ of void material is
the same as that of the solid material. The orientation of layer $n$ is
described by its normal $\bm{n}_{n}$ and tangent $\bm{t}_{n}$, as depicted in
fig:rank3Diag. Through this process, the volume fraction of a rank-$3$ cell
with unit size can be described by

$\rho=1-\prod_{n=1}^{3}(1-\alpha_{n}).$ (1)

Alternatively, the volume fraction can be described by

$\rho=\rho_{1}+\rho_{2}+\rho_{3}.$ (2)

Here, $\rho_{1}=(1-\alpha_{3})(1-\alpha_{2})\alpha_{1}$,
$\rho_{2}=(1-\alpha_{3})\alpha_{2}$, and $\rho_{3}=\alpha_{3}$.

An interpretation of a rank-$3$ laminate as a triangular is illustrated in
fig:rank3Diagf. A general rank-$3$ laminate can be described by six
parameters, i.e., three width variables ($\alpha_{1},\alpha_{2},\alpha_{3}$)
and three orientation variables ($\theta_{1},\theta_{2},\theta_{3}$). However,
allowing the orientation variables to change independently of each other may
result in triangles with bad aspect ratios (i.e., sliver triangles). To
improve the geometric regularity of the de-homogenized triangular lattice, we
propose to reduce the design space by allowing only equilateral triangles,
i.e., by letting $\theta_{2}=\theta_{3}+\frac{\pi}{3}$ and
$\theta_{1}=\theta_{3}+\frac{2\pi}{3}$. Consequently, the independent design
variables become $\alpha_{1},\alpha_{2},\alpha_{3}$ and $\theta_{3}$. The
restriction of orientation valuables has a consequence on the mechanical
performance of optimized structures. Our experiments demonstrate that the
decrease in stiffness under this restriction is mild — less than $5\%$ in
various tests. This restriction allows for significant improvement in the
geometric regularity, and is also beneficial for numerical convergence, as it
reduces the discontinuities in the optimized laminate configurations.

###  2.2 Homogenization-based Topology Optimization

In topology optimization, the objective function quantifies the design target
with respect to the corresponding design variables. The optimization process
minimizes the objective function by iteratively adjusting the design
variables. For designing stiff structures, the objective function usually
considers the strain energy, also known as compliance. For multiple loading
cases the contributions from different loads need to be taken into
consideration concurrently. Thus, a weighted summation of the strain energy
subject to each loading condition is used as a compliance metric ($C$):

$C=\sum\limits_{q=1}^{M}w_{q}\bm{F}_{q}^{T}\bm{U}_{q}.$ (3)

Here, $M$ is the number of loading conditions, $\bm{F_{q}}$ is the nodal force
vector of the $q-$th loading condition and $\bm{U_{q}}$ is the corresponding
displacement vector obtained by solving the static equilibrium equation using
finite element analysis. $w_{q}$ is the weighting factor of the $q-$th loading
condition, with $w_{q}=\frac{1}{M},\>q=1:M$ in our implementation.

Let $O$ be the objective function, which is a function concerning the
compliance metric $C$. By taking triangular lattices as design parameters, the
homogenization-based topology optimization is formulated as

$\displaystyle\displaystyle\min\limits_{\bm{\alpha}_{1},\bm{\alpha}_{2},\bm{\alpha}_{3},\bm{\theta_{3}}}\quad$
$\displaystyle
O(\bm{\alpha}_{1},\bm{\alpha}_{2},\bm{\alpha}_{3},\bm{\theta_{3}}),$ (4)
$\displaystyle\displaystyle\mathrm{subject\,to}\quad$
$\displaystyle\bm{K}(\bm{\alpha}_{1},\bm{\alpha}_{2},\bm{\alpha}_{3},\bm{\theta_{3}})\bm{U}_{q}=\bm{F}_{q},\quad
q=1:M,$ (5)
$\displaystyle\frac{1}{N_{e}}\sum\limits_{e}\bm{\rho}_{e}-V_{F}\leq 0,$ (6)
$\displaystyle
L_{min}\leqslant\bm{\alpha}_{1},\bm{\alpha}_{2},\bm{\alpha}_{3}\leqslant
L_{max},$ (7) $\displaystyle-4\pi\leqslant\theta_{3}\leqslant 4\pi.$ (8)

Here, $\bm{K}$ is the global stiffness matrix assembled from the element
stiffness matrices $\bm{K}_{e}$.
$\bm{K}_{e}=\int_{\Omega_{e}}\bm{B}^{T}\bm{S}\bm{B}\,\mathrm{d}x$, with
$\bm{B}$ being the element strain matrix and $\bm{S}$ the elasticity tensor of
the rank-3 laminate. The elasticity tensor of the rank-n laminate can be
determined analytically [13, 17, 12]. In Sec:ElasticityTensor we summarize the
derivation of $\bm{S}$ for the rank-3 laminate. $N_{e}$ is the number of
simulation cells, $\bm{\rho}_{e}$ corresponds to the relative density value of
each cell (Eqn:VF). $V_{F}$ is the given threshold of volume fraction that
controls the material budget. $L_{min}$ and $L_{max}$ describe the lower and
upper bounds of the width values during optimization, respectively, with
$0<L_{min}<L_{max}<1$. The lower bound $L_{min}$ is introduced to avoid the
appearance of large regions that are fully empty in the design domain, and the
upper bound $L_{max}$ is used to avoid fully solid regions. We set the range
of $\bm{\theta_{3}}$ according to [17].

A hybrid scheme is utilized to solve the optimization problem. The three width
variables $\alpha_{1},\alpha_{2},\alpha_{3}$ are updated by a standard
optimality criteria [2], and the orientation design variable $\theta$ is
updated by the _Method of Moving Asymptotes_ (MMA) [26]. We found that a
moving step size of $0.01$ for the width variables and $\frac{\pi}{180}$ for
the orientation variable resulted in a good balance between convergence speed
and quality of the optimized structure. A standard convolution filter [4, 28]
is applied to the width and orientation design variables for improved
numerical robustness.

####  2.2.1 Orientation regularization

As mentioned in Sec:Rank-N, when using an anisotropic rank-$3$ laminate the
optimal solution is not unique. I.e., while the design variables of each cell
can locally converge to the optimum, they might not form a continuous design
specification globally. Since this prevents the generation of a consistent
structure in de-homogenization, regularizations are needed to counteract or at
least mitigate this non-uniqueness issue.

We adapt the orientation regularization method in [16] to alleviate the
potential discontinuity of the optimized orientation $\bm{\theta}_{3}$. This
is achieved by including a penalty term in the objective function.
Specifically, we use a per-edge penalty function, i.e., for two adjacent
elements $e_{i}$ and $e_{j}$ with shared edge $e_{ij}$

$P^{\theta}_{e_{ij}}=\frac{1}{2}-\frac{1}{2}\cos\left(h_{0}\cdot(\theta_{3}(e_{i})-\theta_{3}(e_{j}))\right).$
(9)

For each cell $e_{i}$, its adjacent cells are defined as the cells sharing an
element edge with $e_{i}$. $P^{\theta}_{e_{ij}}$ with values between 0 and 1
quantifies the orientation deviation between two adjacent cells. Considering
only the orientation variable, a rotation of an equilateral triangle by
$\frac{\pi}{3}$ is invariant. Thus we set $h_{0}=6$. The orientation deviation
between adjacent cells becomes zero if
$\theta_{3}(e_{i})-\theta_{3}(e_{j})=k\frac{\pi}{3}$ ($k$ is an integer).
Accordingly, the largest deviation occurs when
$\theta_{3}(e_{i})-\theta_{3}(e_{j})=k\frac{\pi}{3}+\frac{\pi}{6}$, leading to
$P^{\theta}_{e_{ij}}=1$.

Let $N^{i}_{e}$ be the number of cells adjacent to the considered cell
$e_{i}$, and $P^{\theta}_{e_{ij}}$ the orientation deviation of a cell $e_{i}$
to its $j-$th adjacent cell. Then, a global measure of the orientation
deviation across the design domain is given by

$P^{(\theta)}=\displaystyle\sum\limits_{i=1}^{N_{e}}\displaystyle\sum\limits_{j=1}^{N^{i}_{e}}P^{\theta}_{e_{ij}}.$
(10)

In the end, the objective function for optimization is given by a weighted
summation of the compliance and the orientation regularization

$O(\bm{\alpha}_{1},\bm{\alpha}_{2},\bm{\alpha}_{3},\bm{\theta_{3}})=W\frac{C}{C^{*}}+(1-W)\frac{P}{P^{*}}.$
(11)

The weighting factor $W$ takes on values in $\left(0,\>1\right]$. $C^{*}$ and
$P^{*}$, respectively, refer to the compliance and the orientation
regularization corresponding to the initial values of the design variables.
Orientation regularization can be disabled by setting $W=1$. In all of our
experiments, we have set $W=0.5$.

####  2.2.2 Initialization

The outcome of the non-convex optimization depends on the initial values of
the width and orientation variables $\Phi$. Different initializations may lead
to different density layouts and also affect the convergence behavior of the
optimization process. In previous works, these dependencies have been
demonstrated for single [31, 32] and multiple [16] loading conditions. The
initialization is particularly important since the optimal solution of a
rank-3 laminate is not unique, i.e., different values of the width and
orientation variables may lead to the same elasticity tensor [27].

We take inspiration from the work by Jensen et al.[16], where the initial
principal stress directions are used to construct the initialization. We adapt
it for our triangular infill model. This adaptation considers that only the
orientation of layer 3 is considered independently in our framework, i.e.,
only a single direction field is involved. Specifically, we first compute the
principal stress fields of the fully solid domain subject to the $M$ loading
conditions respectively, resulting in $2\times M$ principal stress values and
the corresponding principal stress directions for each simulation cell. For
each simulation cell, then, we choose the principal stress direction
$\bm{\theta}_{p}$ vector whose associated absolute principal stress value is
largest as the initialization of the orientation of layer 3.

Figure 3: (a) Streamlines generated from the initial orientation field. (b)
Optimization history of compliance and orientation regularization over the
entire optimization process.

It’s worth noting that the orientation variable in rank-3 is represented by
the normals of each layer, while the ‘dominant’ principal stress direction is
supposed to indicate the layering orientation, i.e., the layer tangent. Thus,
$\bm{\theta}_{p}+\frac{\pi}{2}$ is used as the initial value of
$\bm{\theta}_{3}$ in the optimization. For the initialization of the width
variable, we simply let $\bm{\alpha}_{n}=\frac{V_{F}}{3},\>n=1:3$.
fig:SolvingOptia shows the principal stress-based initialization
$\Phi_{stress}$ via a set of streamlines generated from the direction field
$\bm{\theta}_{p}+\frac{\pi}{2}$. The optimization histories of the compliance
$C$ and the orientation regularization $P_{\theta}$ are shown in
fig:SolvingOptib. The optimized design is shown in fig:Overview. For
comparison, fig:SolvingOptib also shows the optimization history of the
compliance when a uniform initialization $\Phi_{uniform}$ is used, where
$\theta_{3}=\frac{\pi}{2}$. Here, to release more design space, orientation
regularization is omitted. It is seen that the optimization process converges
to a smaller compliance value when using $\Phi_{stress}$ instead of
$\Phi_{uniform}$.

fig:SolvingOptia shows that the proposed initialization does not necessarily
yield a globally continuous direction field. This, however, does not pose a
problem in our approach. On the one hand, due to the use of orientation
regularization, the direction field will be smoothed to the maximum extent
during optimization. On the other hand, our proposed geometry-based de-
homogenization is less sensitive to the potential local discontinuity.

##  3 De-homogenization using Field-aligned Triangulation

The proposed de-homogenization converts the spatially varying optimized
specifications ($\alpha_{1},\alpha_{2},\alpha_{3}$, $\theta_{3}$) into an
explicitly represented structural design. This consists of three steps.
Firstly, we create an equilateral triangle mesh with its edges following the
layer orientations $\theta_{3}$. Secondly, we assign a thickness to each edge
of the mesh, determined based on the optimized layer widths
($\alpha_{1},\alpha_{2},\alpha_{3}$), creating a spatially graded triangular
lattice structure. Last, we create a global lattice structure by correcting
the inconsistency of edge thickness at the intersections, i.e., the nodes of
the triangle mesh.

###  3.1 Field-aligned Triangulation

In the first step, instead of employing streamlines to obtain the mesh [29],
which is hard to capture the common intersections for any three streamlines,
we use a field-aligned triangulation algorithm [15] to generate an
orientation-aligned equilateral triangle mesh. Here, we give an overview of
this triangulation method and fit it into the current context. For the
detailed mechanism of this method, refer to the original paper [15]. The
complete pipeline of this triangulation method is composed of three subsequent
stages:

  * •

_Orientation field optimization_. A direction field is computed. It consists
of a set of directions to which the edges of the output mesh should be
aligned. The smoothness of the orientation field is achieved by a Gauss-Seidel
process.

  * •

_Position field optimization_. A local parameterization (vertex positions) is
computed so that the edges between vertices are aligned with the optimized
orientation field.

  * •

_Mesh extraction_. The computed orientation and position fields are turned
into a graph structure from which a triangular mesh is constructed.

Figure 4: (a) The domain discretization (background) and streamlines in the
three different orientation fields of an optimized 6-RoSy field. (b) The
randomly distributed positions (and respective directions) used for position
optimization. The distribution is generated by considering a pre-defined edge
length of the output mesh and the resolution of the input grid. (c) The
position field after optimization and the representative positions (green
dots) used to construct the output mesh. (d) The output mesh whose edges align
with the optimized orientation field. Orientations are indicated by arrows,
and edge colors distinguish corresponding layer orientations. To reduce visual
clutter, only three of the six pairwise bidirectional 6-RoSy field directions
are illustrated.

In the original field-aligned triangulation, the direction field is a smooth
$6$-RoSy field111The $6$-RoSy field refers to a $6$-rotational symmetric
field, which intuitively represents phenomena that are invariant under
rotations of an integer multiple $\frac{2\pi}{6}$ [19, 23]. that aligns with
the features on the boundary domain, obtained by optimizing from a set of
randomly distributed vectors. In the context of de-homogenization, the
direction field represents the optimized layer orientations. Thus, the
original method is adjusted to take the optimized layer orientations as input,
and then perform position field optimization and mesh extraction sequentially.
This is illustrated in fig:meshingDiagram. Given a grid mesh discretization of
the domain, guided by our generated orientation field that is visualized by
streamlines, we apply the triangulation method to convert the orientation
field into a triangle mesh with its edges following the orientation field. In
the following, we elaborate on the conversion of this

###  3.2 Anisotropic triangle mesh

The de-homogenization is conducted in a per-triangle manner. As shown in
fig:DeHomoDiagrama, each triangular element covers a certain region in the
domain. The target deposition ratio $v_{k}^{*}$ of the $k$-th triangular
element is measured by $\frac{\rho_{k}}{M_{k}}$, where $M_{k}$ is the number
of cells located in the region covered by the triangular element, and
$\rho_{k}$ is the sum of the density values over all these cells. The
available material in each triangular region should be re-distributed so that
(i) a binary material layout is generated, (ii) a continuous transition at the
element boundaries is obtained, and (iii) the specific depositions (layer
widths optimized with respect to the object’s compliance) on each oriented
layer of the cells is reflected in the binary material layout. I.e., the
underlying anisotropy of the optimized cells should be preserved in the de-
homogenized lattice elements. In cases where the triangular mesh has a similar
or higher resolution than the rank-3 simulation cells, the specifications of
the cells are resampled using bilinear interpolation so that each triangular
element covers at least a certain number of cells, (e.g., $\geq 10$).

Figure 5: Diagrammatic view of per-lattice de-homogenization. (a) The field-
aligned triangular element with the covered cells. $\tilde{\theta}^{j}_{i}$
refer to the tangents of layers $i\in{1,2,3}$ in element $j$. Three arrow
glyphs on each cell indicate the tangents. Corresponding width variables
$\alpha^{j}_{i}$ are encoded into the arrow lengths. (b) The representative
values of widths and orientations of the cells in the triangular element. (c)
Relation of the representative values of orientations of three layers to the
corresponding edges of the triangular element by detecting their orientation
deviations. (d) The de-homogenized triangular lattice element whose unique
edge thicknesses reflect the representative values of widths. The solid
regions correspond to the deposition ratio indicated in (a).

De-homogenization, in particular, requires that the de-homogenized triangular
element maintains the deposition ratio after de-homogenization, i.e.,

$\frac{\rho_{k}}{M_{k}}=\frac{A-A_{void}(t)}{A}$ (12)

Here, $A$ is the area of the triangular element, and $A_{void}$ is the area of
the void region left in the triangular element after de-homogenization
(fig:DeHomoDiagramd). $t$ is the thickness parameter that needs to be
determined to satisfy the equation. The void area $A_{void}$ is up to the edge
thicknesses of the triangular element. To convey the anisotropy via the edge
thickness, $t$ is scaled for each edge by the scaling factors
$\alpha_{1}^{a}$, $\alpha_{2}^{a}$ and $\alpha_{3}^{a}$ (fig:DeHomoDiagramb),
which are determined by the representative values of the actual material
deposition on each layer.

For each cell $j$ in a triangular element, the material depositions of each
layer are evaluated according to Eqn:VFAlternative, i.e.,

$\begin{array}[]{l}\hat{\alpha}_{1}^{j}=(1-\alpha_{3}^{j})(1-\alpha_{2}^{j})\alpha_{1}^{j}\\\
\hat{\alpha}_{2}^{j}=(1-\alpha_{3}^{j})\alpha_{2}^{j}\\\
\hat{\alpha}_{3}^{j}=\alpha_{3}^{j}.\end{array}$ (13)

$\alpha_{1}^{j}$, $\alpha_{2}^{j}$, and $\alpha_{3}^{j}$ are the corresponding
width variables of layer 1, 2, and 3 of the $j$-th cell, respectively. From
Eqn:DefScalingFactors, the representative values ($\alpha_{1}^{a}$,
$\alpha_{2}^{a}$ and $\alpha_{3}^{a}$) of the actual material deposition along
each layering orientation in a triangular element are given by
$\alpha_{1}^{a}=\sum_{j=1}^{M_{k}}{\hat{\alpha}_{1}^{j}}$,
$\alpha_{2}^{a}=\sum_{j=1}^{M_{k}}{\hat{\alpha}_{2}^{j}}$, and
$\alpha_{3}^{a}=\sum_{j=1}^{M_{k}}{\hat{\alpha}_{3}^{j}}$, respectively. In
practice, $\alpha_{1}^{a}$, $\alpha_{2}^{a}$ and $\alpha_{3}^{a}$ are
normalized by $\max{(\alpha_{1}^{a},\alpha_{2}^{a},\alpha_{3}^{a})}$.

Relating each edge of the triangular element to the corresponding layer
orientation is required to convey the anisotropy of the underlying cells.
However, since the orientations vary across the cells, in general,
representative orientations first need to be computed. Here, we make use of a
weighted average strategy, i.e., the actual material deposition on each layer
is used to weight the corresponding orientation. With $\tilde{\theta}_{i}^{j}$
referring to the tangent of layer $i$ of the $j$-th cell within the triangular
element, the weighted average is computed as

$\tilde{\theta}_{i}^{a}=\frac{\displaystyle\sum\limits_{j=1}^{M_{k}}\tilde{\theta}_{i}^{j}\hat{\alpha}_{i}^{j}}{\displaystyle\sum\limits_{j=1}^{M_{k}}\hat{\alpha}_{i}^{j}},\quad
i=1:3.$ (14)

We let the mesh edges correspond to $\tilde{\theta}_{1}^{a}$,
$\tilde{\theta}_{2}^{a}$ or $\tilde{\theta}_{3}^{a}$, depending on which edge
they have the least directional deviation, see fig:DeHomoDiagramc.

###  3.3 Lattice structure

When considering triangular elements the per-lattice de-homogenization does
not necessarily produce a well-posed structural design. This is because gaps
are left at the mesh vertices when the thicknesses of mesh edges sharing a
vertex vary significantly, see fig:FillLatticeVtxGapa. To avoid this issue, we
introduce additional small triangular patches to fill these gaps.

Figure 6: The structural designs after de-homogenization. (a) Original result
by only directly performing the per-lattice de-homogenization. (b) Improved
result by the introduced triangular patches.

In particular, we first determine those target vertices that need to be fixed.
This requires to inspect — and repair upon validation — three vertices in each
lattice. At the vertex $p_{1}$ of a considered lattice, a gap-filling triangle
$S_{\hat{p}_{1}\hat{p}_{2}\hat{p}_{3}}$ is considered, see
fig:FillLatticeVtxGapb. $p_{2}$ and $p_{3}$ are the vertices of the
corresponding internal empty triangles of the two adjacent lattices that share
the same edges with the considered lattice. Then, it is checked whether the
corresponding vertex of the internal empty triangle of the considered lattice
$\hat{p}_{1}$ is within the triangle $S_{p_{1}p_{2}p_{3}}$. If $\hat{p}_{1}$
is outside of $S_{p_{1}p_{2}p_{3}}$, then no gap is to be filled. Otherwise, a
solid triangle $S_{\hat{p}_{1}\hat{p}_{2}\hat{p}_{3}}$ is placed, where
$\hat{p}_{2}$ and $\hat{p}_{2}$ are the intersection points of $p_{2}p_{3}$
and the corresponding edges of the internal empty triangle of the considered
lattice. In this way, it is ensured that there is no material overlap in the
final structural design. The improved structural design is shown in
fig:FillLatticeVtxGapb. The additional material consumption caused by gap-
filling triangles is compensated by scaling down the deposition ratio in
Eqn:equilibriumDeposition.

##  4 Results

In the following, we evaluate the proposed method with a variety of
experiments. Firstly, a parameter study is performed to showcase the behavior
of our method with respect to different settings. This includes the different
strategies for setting the initialization, with/without the orientation
regularization, and the different resolutions of the triangular
parametrization. Secondly, we compare our results to those of alternative
approaches, including the uniform triangular lattice structure, the porous
infill produced by density-based topology optimization under a local volume
constraint [34], and the alternative de-homogenization approach [16]. We then
consider two practical examples to further verify the usability of our method.
Finally, the adaptivity to situations with concurrently varying loading and
constraint conditions is demonstrated via a specific test case. To
differentiate our approach from alternative approaches, our generated
structural design is presented by a density layout in black together with the
corresponding triangular mesh in green. In this way, it is highlighted that
our approach results in both a structural design and an explicit geometric
format.

In all examples, the design domains are discretized by Cartesian grids with
unit size. Young’s modulus and Poisson’s ratio are set to 1.0 and 0.3,
respectively. Except for the field-aligned triangular meshing part, which is
performed in the open source software _Instant-meshes_
222https://github.com/wjakob/instant-meshes present by Jakob _et al._ [15],
all other parts of the approach, including homogenization-based topology
optimization and de-homogenization, is implemented in MatLab. We let the
optimization process terminate after 300 iterations. The lower and upper
bounds of the relative width of each cell are set to $L_{min}=0.1$ and
$L_{max}=0.5$ if there is no special explanation. Correspondingly, the
relative density value varies from 0.271 to 0.875 according to Eqn:VF. All
experiments have been carried out on a desktop PC with an Intel Xeon CPU at
3.60GHz.

Performance statistics for all examples can be found in
tab:qualitativeStatistics. In this table, we also provide the resulting
compliance $C^{*}$ of homogenization-based topology optimization without
restricting the layer orientations, i.e., all three orientation variables
($\theta_{1:3}$) can be adjusted arbitrarily during optimization. All other
settings are kept the same as the ones directing to $C_{0}$. By comparing
$C_{0}$ to the corresponding $C^{*}$, we see that our additional orientation
constraint does increase the compliance, whereas, it’s within $5\%$ (measured
by $\frac{C_{0}-C^{*}}{C^{*}}$) in all of the experiments using the proposed
$\Phi_{stress}$ and orientation regularization.

Table 1: Performance statistics for homogenization-based topology optimization
and field-aligned triangulation-guided de-homogenization. $C_{0}$ and $V_{0}$,
respectively, are the compliance and volume fraction of the optimal layout
resulting from homogenization-based topology optimization. Accordingly, $C$
and $V$ correspond to the de-homogenized result. The design deviation of the
de-homogenized result from the optimal layout is evaluated by
$\xi=\frac{C\cdot V-C_{0}\cdot V_{0}}{C_{0}\cdot V_{0}}$. $t_{0}$ and $t$ are
the time consumptions of homogenization-based topology optimization and de-
homogenization, respectively. Provided the same volume fraction, $C^{*}$ is
the resulting compliance of homogenization-based topology optimization without
restricting the layer orientations.

2*Examples |  | Homogenization-based Opti. | De-homogenization |  2*$\xi$
---|---|---|---|---
|  | $C^{*}$ | $C_{0}$ | $V_{0}$ |  $t_{0}$ (s) | $C$ | $V$ |  $t$ (s) |
4*Femur |  fig:Overviewh | 73.34 | 76.86 | 0.50 | 366 | 84.37 | 0.50 | 25 | $9.78\%$
|  fig:parameterStudyc | 75.26 | 83.25 | 0.50 | 356 | 92.17 | 0.50 | 27 | $10.93\%$
|  fig:parameterStudyf | 71.25 | 73.30 | 0.50 | 354 | 84.19 | 0.50 | 30 | $14.86\%$
|  fig:parameterStudyg | 73.34 | 76.86 | 0.50 | 366 | 85.01 | 0.50 | 29 | $10.60\%$
Wheel |  fig:wheelb | 7.78 | 8.03 | 0.50 | 280 | 8.30 | 0.50 | 27 | $3.29\%$
Beam |  fig:comp2Projb | 2.61 | 2.66 | 0.30 | 212 | 2.78 | 0.31 | 21 | $8.00\%$
Bracket 1 |  fig:additionalExpsc | 19.46 | 19.92 | 0.50 | 326 | 20.45 | 0.50 | 28 | $2.67\%$
Bracket 2 |  fig:additionalExpsf | 26.67 | 27.45 | 0.50 | 641 | 29.50 | 0.50 | 37 | $7.47\%$
Triangle |  fig:MultiBCsResultb | 11.05 | 11.07 | 0.50 | 539 | 12.12 | 0.50 | 33 | $10.15\%$

In the design of methodology, our method performs homogenization-based
topology optimization on a simulation mesh ($\Re_{sim}$) with low resolution
to obtain the optimized specification, and the de-homogenized structural
designs are described in a compact geometric format, i.e., a triangular mesh
with each edge being assigned a unique thickness. Herein, to facilitate
evaluating and comparing the obtained structural designs, we also introduce a
high-resolution validation mesh ($\Re_{fine}$). This is because the obtained
structural design can be composed of very fine geometric details, which can
only be captured by a sufficiently high resolution in compliance evaluation
using finite element analysis. To alleviate the potential inconsistency of
boundary conditions among different resolutions, we borrow the idea of mesh
hierarchy in the geometric multi-grid solver that is frequently employed in
high-resolution topology optimization [35]. Specifically, the design domain is
discretized into a fine Cartesian mesh as $\Re_{fine}$, from which, a
relatively coarse resolution grid ($\Re_{sim}$) and corresponding boundary
conditions are constructed for homogenization-based topology optimization and
the subsequent de-homogenization. This is achieved by the restriction
operation in the geometric multi-grid solver. When evaluating the compliance
of the obtained design, the triangular lattice structure is discretized
according to $\Re_{fine}$, by identifying the cells of $\Re_{fine}$ that are
located on the obtained triangular lattice edges.

To ensure consistent boundary conditions and discretization when comparing
$C_{0}$ and $C$, $C_{0}$ is re-evaluated after optimization by projecting the
optimized specifications onto $\Re_{fine}$ through bilinear interpolation. For
density-based topology optimization, $\Re_{fine}$ is directly used in the
simulation to provide sufficient design space. Thus, it is ensured that the
resolution and boundary conditions are consistent when comparing our method to
other alternatives. tab:PD provides information regarding the resolution and
number of elements of $\Re_{fine}$ and $\Re_{sim}$, as well as the number of
triangular lattices involved in de-homogenization. The number of triangular
lattices is determined as follows: First, a reference length for the
triangular lattice edges is specified. This reference considers the preferred
fineness of the structural design and a possible size limitation stemming from
a downstream manufacturing task. Since the triangular lattices are from a set
of approximately equilateral triangular elements, a reference area of the
triangular lattices can then be derived from the reference length. Finally,
the number of triangular lattices is computed by dividing the total area of
the design domain by the reference area.

Table 2: Statistics for the validation meshes ($\Re_{fine}$) used to evaluate
the final structural designs, the simulation meshes ($\Re_{sim}$) used for
homogenization-based topology optimization, and the triangular lattice
elements ($\\#$_Lattice_) used for de-homogenization.

2*Examples |  | $\Re_{fine}$ | $\Re_{sim}$ |  2*$\\#$_Lattice_
---|---|---|---|---
|  | $\\#Resolution$ | $\\#Elements$ | $\\#Resolution$ | $\\#Elements$ |
Femur |  fig:Overviewa | $920\times 1200$ | 632,468 | $116\times 150$ | 10,118 | 684 (fig:Overviewf)
Wheel |  fig:wheela | $1600\times 1600$ | 1,930,188 | $100\times 100$ | 7,740 | 1,637 (fig:wheelb)
Beam |  fig:comp2Proja | $1600\times 800$ | 1,280,000 | $100\times 50$ | 5,000 | 699 (fig:comp2Projb)
Bracket 1 |  fig:additionalExpsa | $2296\times 1328$ | 2,181,496 | $144\times 84$ | 8,749 | 2,101 (fig:additionalExpsc)
Bracket 2 |  fig:additionalExpsd | $2332\times 1000$ | 1,144,564 | $292\times 126$ | 18,271 | 1,167 (fig:additionalExpsf)
Triangle |  fig:MultiBCsResulta | $1360\times 1184$ | 790,132 | $170\times 148$ | 12,652 | 1,065 (fig:MultiBCsResultb)

##### Design verification.

With the ‘Femur’ model shown in fig:Overview as an example, de-homogenization
is performed using several test settings to verify our design decisions.
Specifically, the results when using a uniform initialization $\Phi_{uniform}$
and of original homogenization-based topology optimization, i.e., no
orientation regularization is used, are given fig:parameterStudya-c.
fig:parameterStudyd-f shows the results the principal stress-guided
initialization $\Phi_{stress}$ is used without orientation regularization. In
fig:parameterStudyg-i, the settings in homogenization-based topology
optimization remain the same as in fig:Overview, i.e., the initialization
$\Phi_{stress}$ and orientation regularization is used, but a higher
resolution of the triangular parametrization is chosen. fig:parameterStudy(c),
(f), and (i) use the consistent resolution for creating the field-aligned
triangular mesh.

Figure 7:  The homogenization-based results (optimized orientations and
density layout) and the corresponding de-homogenization results of three
different settings are demonstrated in each box. In the left box (a-c), the
homogenization-based topology optimization is conducted with the uniform
initialization ($\Phi_{uniform}$). In the middle one (d-f), $\Phi_{uniform}$
is replaced with the proposed initialization $\Phi_{stress}$, whereas other
settings are kept the same. In the right box, both the $\Phi_{stress}$ and
orientation regularization are considered in the homogenization-based topology
optimization. Among (c), (f), and (i), the resolution of the used triangular
meshes is consistent.

The test results in fig:parameterStudy and fig:Overview, and the statistics in
tab:qualitativeStatistics demonstrate that the design deviation between
homogenization-based topology optimization and de-homogenization results do
not exceed 15%, and most are around at 10%, regardless of whether the
optimized orientation field is continuous or not. This demonstrates that the
proposed de-homogenization has relaxed requirements on the optimized
specifications.

In terms of design quality measured by compliance and geometric regularity,
however, different design decisions give rise to different results. Compared
to homogenization-based topology optimization using $\Phi_{uniform}$
(fig:parameterStudya, b), we see that $\Phi_{stress}$ does improve the
convergence behavior (fig:parameterStudyd, e). The compliance is reduced from
$83.25$ to $73.30$ (fig:parameterStudyb, e). However, $\Phi_{stress}$ incurs
more discontinuities in the optimized orientation field. Though the de-
homogenized lattice structure corresponding to $\Phi_{uniform}$
(fig:parameterStudyc) is less stiff than the one obtained with $\Phi_{stress}$
(fig:parameterStudyf), the former design gives better geometric regularity.
The difference in the setting for fig:parameterStudyd, e and
fig:parameterStudyg, h is that the orientation regularization is introduced in
the latter. As a result, fig:parameterStudyg behaves more continuously than
fig:parameterStudyd), yet fig:parameterStudye) is slightly stiffer than
fig:parameterStudyh. This is because the orientation regularization acts like
an additional constraint in the optimization process. The same phenomenon is
kept in the de-homogenized results. Interestingly, the corresponding de-
homogenization result fig:parameterStudyf is only slightly stiffer than
fig:parameterStudyi compared to their difference in the homogenization-based
results. This is regardless of the introduced geometric regularity. This
behaviour is interesting due to the severe discontinuity of the optimized
orientation field (fig:parameterStudyd) that still hurts the precision of de-
homogenization, though it does not necessarily lead to failure in the de-
homogenization process.

From fig:parameterStudyi and fig:Overviewf, we further see that the proposed
de-homogenization is stable with respect to changes of the triangular mesh
resolution. By comparing fig:parameterStudyc, f, and i, one can observe that
the homogenization-based topology optimization combined with proposed
$\Phi_{stress}$ and orientation regularization directs to a more desirable
structural design, which achieves better stiffness and geometric regularity
concurrently. Thus, this setting is used in the following experiments.

##### Comparison analysis.

We first compare our results to the uniform triangular lattice structure and
porous infill structure. The latter is produced via density-based topology
optimization under a local volume constraint [34].
fig:comparison2existing_femur and fig:wheel show that the uniform triangular
lattice structure always gives rise to the highest compliance. The triangular
lattice structure is not significantly stiffer than the porous infill
structure. However, density-based topology optimization is a full-scale method
and requires a sufficiently high simulation resolution to obtain a fine
structural design. This significantly increases the computational cost. In
addition, the method does not always converge [31]. This is shown in
fig:comparison2existing_femurb, where a binary layout is still unavailable
even after 1000 iterations.

Figure 8: (a) The uniform triangular lattice structure. (b) The porous infill
structure from density-based topology optimization under the local volume
constraint. (c) The proposed de-homogenization result. All examples share the
same problem description. Figure 9: (a) A triangular design domain with 2
different loading cases, the normal force $F_{n}$ and shear force $F_{t}$. (b)
Result by our method. (c) Uniform triangular lattice structure. (d) Porous
infill structure after 2000 iterations.

While our method is primarily developed to generate infill structures that
spread across the entire design domain, it can be adapted to generate lattice
structures that cover only part of the design domain, by setting
$L_{min}=1.0\times 10^{-6}$, $L_{max}=1.0$. The obtained structural layout is
a subset of the initial design domain. With this adaptation, we evaluate the
trade-off between the structural performance and regularity of the lattice,
taking the results from the image-based de-homogenization method [16] as a
reference. The design domain and load conditions are illustrated in
fig:comp2Proja. The results from our method and the image-based approach are
shown in fig:comp2Projb and c, respectively. The structure from our method
consists of equilateral triangles with varying thicknesses, while that from
the image-based approach exhibits a large number of intersections of
differently oriented beam-like substructures. The improved regularity in
lattice structures is achieved with a compromise in the stiffness. The
difference in stiffness is 8.2% measured by
$\frac{C_{ours}-C_{alternative}}{C_{alternative}}$ when the volume fraction is
the same, with $C_{ours}=2.78$ and $C_{alternative}=2.57$.

Figure 10: (a) Beam with 5 different loads. (b) De-homogenized structural
design by our method. (c) The corresponding result from the image-based de-
homogenization, reproduced from [16].

##### Mechanical parts.

With ‘Bracket-1’ and ‘Bracket-2’ we consider two mechanical parts for
emphasizing the basic features of the proposed de-homogenization scheme (see
Fig. 11). In particular, the results are in line with those obtained for
‘femur’, concerning compliance and geometric regularity.

Figure 11: (a), (d) Problem descriptions of bracket 1 and bracket 2,
respectively. (b), (e) The corresponding density layouts generated by
homogenization-based topology optimization. (c), (f) The resulting de-
homogenized structural designs.

##### Multiple boundary conditions.

We further demonstrate the adaptivity of the proposed method to the
optimization problems with totally different boundary conditions, i.e., the
fixations and loads vary concurrently. Here we consider an equilateral
triangular design domain for infill optimization (fig:MultiBCsResulta).

In fig:MultiBCsResult, we show the optimization and de-homogenization results,
where the design deviation between the de-homogenized result and the
corresponding homogenization-based topology optimization result is at
$\xi=10.15\%$ and close to other examples. In addition, the triangular lattice
structure exhibits good regularity.

Figure 12: (a) The equilateral triangular design domain and the corresponding
boundary conditions. Here the thick lines in black indicate the fixation
positions. The loads are shown in different colors, each load is outward of
the design domain and perpendicular to the corresponding domain boundary. (b)
The optimized density layout and the corresponding layer orientations shown by
streamlines. (c) The de-homogenized triangular lattice structure.

##  5 Conclusion and Future Work

In this paper, we have introduced an innovative geometry-based de-
homogenization technique for designing and optimizing functionally graded
triangular lattices capable of withstanding various loading conditions. By
conducting simulations and optimizations at a low-resolution discretization of
the design domain, our de-homogenization approach, akin to other methods in
this domain, generates highly detailed structures with significantly reduced
computational demands compared to traditional density-based topology
optimization methods.

Unique to our approach is its ability to explicitly represent the resulting
structural design in a concise geometric format, i.e., a global triangular
mesh where each edge is assigned a distinct thickness. This explicit geometric
representation facilitates downstream processes such as editing and additive
manufacturing. In terms of mechanical performance, the de-homogenized
structure exhibits stiffness comparable to that of the optimal design achieved
through homogenization-based topology optimization. Additionally, when
compared to porous infill generated using local volume constraints, our de-
homogenized structure consistently demonstrates higher stiffness.

In future work, we aim to extend this geometry-based de-homogenization method
to address 3D problems. Just as we simplified a general rank-3 laminate to an
equilateral triangle, we speculate that a general rank-6 laminate in 3D could
potentially be approximated by an equilateral tetrahedron. We plan to verify
this hypothesis and devise geometric techniques for the de-homogenization
process.

## Declaration of competing interest

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

## Availability of code

The demonstration code for homogenization-based topology optimization can be
accessed through https://github.com/333will be made publicly available after
the review process..

## Acknowledgment

We thank Peter Dørffler Ladegaard Jensen from the Technical University of
Denmark for the helpful discussion. This work was supported in part by a grant
from German Research Foundation (DFG) under grant number WE 2754/10-1.

##  Appendix A Elasticity Tensor of Rank-$3$

In a variety of related works [17, 12, 16], the general description of the
elasticity tensor ($\bm{S}$) of rank-$3$ is given as below

$S_{ijkl}=S^{+}_{ijkl}-(1-\rho)\left[(S^{+}_{ijkl}-S^{-}_{ijkl})^{-1}-\frac{\rho(1-v^{2}_{0})}{E^{+}}\left(\displaystyle\sum\limits_{n=1}^{3}P_{n}\Lambda^{n}_{ijkl}\right)\right]^{-1}$
(15)

Here, same as the elasticity tensor of the plane problem, $S_{ijkl}$ is a
fourth-order tensor. $S^{+}_{ijkl}$ and $S^{-}_{ijkl}$ are the elasticity
tensors of the isotropic solid and ’void’ materials and are given by Hooke’s
law

$\bmatrix S_{1111}^{+/-}&S_{1122}^{+/-}S_{1112}^{+/-}\\\
S_{2222}^{+/-}S_{2212}^{+/-}\\\
symS_{1212}^{+/-}=\frac{E^{+/-}}{1-v^{2}_{0}}\bmatrix 1&v0\\\ 10\\\
sym2(1-v_{0})$ (16)

$E^{+}$ and $v_{0}$ are Young’s modulus and Poisson’s ratio of the isotropic
solid material. $P_{n}>0$ depicts the relative contribution of layer $n$ to
the stiffness, which subjects to $\sum_{n=1}^{3}P_{n}=1$, and
$P_{n}=\frac{\rho_{n}}{\rho}$. Combining $P_{n}$ and $\Lambda^{n}_{ijkl}$, the
item $\sum\limits_{n=1}^{3}P_{n}\Lambda^{n}_{ijkl}$ parametrizes all of the
possible properties of rank-$3$ laminate in the entire design space with a
fixed volume fraction. $\Lambda$ is defined as

$\Lambda^{n}_{ijkl}=n^{n}_{i}n^{n}_{j}n^{n}_{k}n^{n}_{l}+\frac{t^{n}_{i}n^{n}_{j}t^{n}_{k}n^{n}_{l}+n^{n}_{i}t^{n}_{j}t^{n}_{k}n^{n}_{l}+t^{n}_{i}n^{n}_{j}n^{n}_{k}t^{n}_{l}+n^{n}_{i}t^{n}_{j}n^{n}_{k}t^{n}_{l}}{2\left(1-v_{0}\right)}$
(17)

Here, $\bm{n}^{n}$ and $\bm{t}^{n}$ are the components of normal and tangent
of layer $n$. Let $\bm{n}^{n}=[\cos\theta_{n}\>\sin\theta_{n}]$ and
$\bm{t}^{n}=[-\sin\theta_{n}\>\cos\theta_{n}]$, here $\theta_{n}$ is the
orientation angle of the normal of layer $n$ to the $x$-axis in a standard
coordinate system, all the components of $\Lambda^{n}_{ijkl}$ are respectively
expressed by

$\begin{array}[]{l}\Lambda^{n}_{1111}=\frac{\cos(4\theta_{n})+4\cos(2\theta_{n})+3}{8}+\frac{1-\cos(4\theta_{n})}{4(1-v_{0})};\\\
\\\
\Lambda^{n}_{2222}=\frac{\cos(4\theta_{n})-4\cos(2\theta_{n})+3}{8}+\frac{1-\cos(4\theta_{n})}{4(1-v_{0})};\\\
\\\
\Lambda^{n}_{1122}=\Lambda^{n}_{2211}=\frac{1-\cos(4\theta_{n})}{8}-\frac{1-\cos(4\theta_{n})}{4(1-v_{0})};\\\
\\\
\Lambda^{n}_{1112}=\Lambda^{n}_{1121}=\Lambda^{n}_{2111}=\frac{2\sin(2\theta_{n})+\sin(4\theta_{n})}{8}-\frac{\sin(4\theta_{n})}{4(1-v_{0})};\\\
\\\
\Lambda^{n}_{1222}=\Lambda^{n}_{2221}=\Lambda^{n}_{2122}=\Lambda^{n}_{2212}=\frac{2\sin(2\theta_{n})-\sin(4\theta_{n})}{8}+\frac{\sin(4\theta_{n})}{4(1-v_{0})};\\\
\\\
\Lambda^{n}_{1212}=\Lambda^{n}_{1221}=\Lambda^{n}_{2112}=\Lambda^{n}_{2121}=\frac{1-\cos(4\theta_{n})}{8}+\frac{1+\cos(4\theta_{n})}{4(1-v_{0})}.\end{array}$
(18)

In practical computation, Eqn:Stensor is usually given in a convenient matrix
form ([17])

$\bm{S}^{[2]}=\bm{S}^{+}-(1-\rho)\left[(\bm{S}^{+}-\bm{S}^{-})^{-1}-\frac{\rho(1-v^{2}_{0})}{E^{+}}\mathcal{M}\right]^{-1}$
(19)

Here, $\mathcal{M}=\displaystyle\sum\limits_{n=1}^{3}P_{n}\bm{\Lambda}^{n}$,
and the matrices $\bm{S^{+}}$, $\bm{S^{-}}$ and $\bm{\Lambda}^{n}$ are
obtained from the corresponding fourth-order tensor components through the
transformation below

$\bm{A}=\bmatrix\frac{1}{2}(A_{1111}+A_{2222})-A_{1122}&A_{1112}-A_{2221}\frac{1}{2}(A_{1111}-A_{2222})\\\
2A_{1212}A_{1112}+A_{2221}\\\ sym\frac{1}{2}(A_{1111}+A_{2222})+A_{1122}$ (20)

Let

$\begin{array}[]{l}m_{1}=\displaystyle\sum_{n=1}^{3}p_{n}\cos(2\theta_{n}),\quad
m_{2}=\displaystyle\sum_{n=1}^{3}p_{n}\sin(2\theta_{n}),\\\
m_{3}=\displaystyle\sum_{n=1}^{3}p_{n}\cos(4\theta_{n}),\quad
m_{4}=\displaystyle\sum_{n=1}^{3}p_{n}\sin(4\theta_{n}).\end{array}$ (21)

Then, the matrix $\mathcal{M}$ is rearranged into a simpler form

$\mathcal{M}=\bmatrix\frac{3-v_{0}-(1+v_{0})m_{3}}{4(1-v_{0})}&-\frac{(1+v_{0})m_{4}}{4(1-v_{0})}\frac{m_{1}}{2}\\\
\frac{3-v_{0}+(1+v_{0})m_{3}}{4(1-v_{0})}\frac{m_{2}}{2}\\\ sym\frac{1}{2}$
(22)

Upon obtaining $\bm{S}$, one can get the components of $S_{ijkl}$ in the
following way

$\bm{S}=\bmatrix S_{1111}&S_{1122}S_{1112}\\\ S_{2222}S_{2221}\\\
symS_{1212}=\bmatrix\frac{1}{2}(S_{11}^{[2]}+S_{33}^{[2]})+S_{13}^{[2]}&-\frac{1}{2}(S_{11}^{[2]}-S_{33}^{[2]})\frac{1}{2}(S_{12}^{[2]}+S_{23}^{[2]})\\\
\frac{1}{2}(S_{11}^{[2]}+S_{33}^{[2]})-S_{13}^{[2]}-\frac{1}{2}(S_{12}^{[2]}-S_{23}^{[2]})\\\
sym\frac{1}{2}S_{22}^{[2]}$ (23)

For details of the derivation process above, we refer to [17].

## References

  * [1] Grégoire Allaire, Perle Geoffroy-Donders, and Olivier Pantz.  Topology optimization of modulated and oriented periodic microstructures by the homogenization method.  Computers & Mathematics with Applications, 78(7):2197–2229, 2019\.  https://doi.org/10.1016/j.camwa.2018.08.007 doi:10.1016/j.camwa.2018.08.007.
  * [2] Martin P Bendsøe.  Optimization of structural topology, shape, and material, volume 414.  Springer, 1995.  https://doi.org/10.1007/978-3-662-03115-5 doi:10.1007/978-3-662-03115-5.
  * [3] Martin Philip Bendsøe and Noboru Kikuchi.  Generating optimal topologies in structural design using a homogenization method.  Computer Methods in Applied Mechanics and Engineering, 71(2):197–224, 1988.  https://doi.org/10.1016/0045-7825(88)90086-2 doi:10.1016/0045-7825(88)90086-2.
  * [4] Blaise Bourdin.  Filters in topology optimization.  International journal for numerical methods in engineering, 50(9):2143–2158, 2001.  https://doi.org/10.1002/nme.116 doi:10.1002/nme.116.
  * [5] Anders Clausen, Niels Aage, and Ole Sigmund.  Exploiting additive manufacturing infill in topology optimization for improved buckling load.  Engineering, 2(2):250–257, 2016.  https://doi.org/10.1016/J.ENG.2016.02.006 doi:10.1016/J.ENG.2016.02.006.
  * [6] Stephen Daynes, Stefanie Feih, Wen Feng Lu, and Jun Wei.  Optimisation of functionally graded lattice structures using isostatic lines.  Materials & Design, 127:215–223, 2017.  https://doi.org/10.1016/j.matdes.2017.04.082 doi:10.1016/j.matdes.2017.04.082.
  * [7] Martin Ohrt Elingaard, Niels Aage, Jakob Andreas Bærentzen, and Ole Sigmund.  De-homogenization using convolutional neural networks.  Computer Methods in Applied Mechanics and Engineering, 388:114197, 2022.  https://doi.org/10.1016/j.cma.2021.114197 doi:10.1016/j.cma.2021.114197.
  * [8] Perle Geoffroy-Donders, Grégoire Allaire, and Olivier Pantz.  3-d topology optimization of modulated and oriented periodic microstructures by the homogenization method.  Journal of Computational Physics, 401:108994, jan 2020.  https://doi.org/10.1016/j.jcp.2019.108994 doi:10.1016/j.jcp.2019.108994.
  * [9] Jeroen P. Groen and Ole Sigmund.  Homogenization-based topology optimization for high-resolution manufacturable microstructures.  International Journal for Numerical Methods in Engineering, 113(8):1148–1163, 2018.  https://doi.org/10.1002/nme.5575 doi:10.1002/nme.5575.
  * [10] Jeroen P. Groen, Florian C. Stutz, Niels Aage, Jakob A. Bærentzen, and Ole Sigmund.  De-homogenization of optimal multi-scale 3D topologies.  Computer Methods in Applied Mechanics and Engineering, 364:112979, jun 2020.  https://doi.org/10.1016/j.cma.2020.112979 doi:10.1016/j.cma.2020.112979.
  * [11] Jeroen P Groen, Jun Wu, and Ole Sigmund.  Homogenization-based stiffness optimization and projection of 2D coated structures with orthotropic infill.  Computer Methods in Applied Mechanics and Engineering, 349:722–742, 2019.  https://doi.org/10.1016/j.cma.2019.02.031 doi:10.1016/j.cma.2019.02.031.
  * [12] Jeroen Peter Groen.  Multi-scale design methods for topology optimization.  DTU Mechanical Engineering, 2018.
  * [13] B. Hassani and E. Hinton.  A review of homogenization and topology optimization i@—homogenization theory for media with periodic structure.  Computers & Structures, 69(6):707–717, 1998.  https://doi.org/https://doi.org/10.1016/S0045-7949(98)00131-X doi:https://doi.org/10.1016/S0045-7949(98)00131-X.
  * [14] Yuming Huang, Guoxin Fang, Tianyu Zhang, and Charlie CL Wang.  Turning-angle optimized printing path of continuous carbon fiber for cellular structures.  Additive Manufacturing, 68:103501, 2023.  https://doi.org/10.1016/j.addma.2023.103501 doi:10.1016/j.addma.2023.103501.
  * [15] Wenzel Jakob, Marco Tarini, Daniele Panozzo, Olga Sorkine-Hornung, et al.  Instant field-aligned meshes.  ACM Trans. Graph., 34(6):1–15, 2015.  https://doi.org/10.1145/2816795.2818078 doi:10.1145/2816795.2818078.
  * [16] Peter Dørffler Ladegaard Jensen, Ole Sigmund, and Jeroen P Groen.  De-homogenization of optimal 2D topologies for multiple loading cases.  Computer Methods in Applied Mechanics and Engineering, 399:115426, 2022.  https://doi.org/10.1016/j.cma.2022.115426 doi:10.1016/j.cma.2022.115426.
  * [17] LA Krog and Niels Olhoff.  Topology and reinforcement layout optimization of disk, plate, and shell structures.  Springer, 1997.  https://doi.org/10.1007/978-3-7091-2566-3_9 doi:10.1007/978-3-7091-2566-3_9.
  * [18] Tsz-Ho Kwok, Yongqiang Li, and Yong Chen.  A structural topology design method based on principal stress line.  Computer-Aided Design, 80:19–31, 2016.  https://doi.org/10.1016/j.cad.2016.07.005 doi:10.1016/j.cad.2016.07.005.
  * [19] Jonathan Palacios and Eugene Zhang.  Rotational symmetry field design on surfaces.  ACM Transactions on Graphics (TOG), 26(3):55–es, 2007.  https://doi.org/10.1145/1276377.1276446 doi:10.1145/1276377.1276446.
  * [20] Ajit Panesar, Meisam Abdi, Duncan Hickman, and Ian Ashcroft.  Strategies for functionally graded lattice structures derived using topology optimisation for additive manufacturing.  Additive Manufacturing, 19:81–94, 2018.  https://doi.org/10.1016/j.addma.2017.11.008 doi:10.1016/j.addma.2017.11.008.
  * [21] O. Pantz and K. Trabelsi.  A post-treatment of the homogenization method for shape optimization.  SIAM Journal on Control and Optimization, 47(3):1380–1398, 2008\.  https://doi.org/10.1137/070688900 doi:10.1137/070688900.
  * [22] Pauli Pedersen.  On optimal orientation of orthotropic materials.  Structural optimization, 1(2):101–106, 1989.  https://doi.org/10.1007/BF01637666 doi:10.1007/BF01637666.
  * [23] Nicolas Ray, Bruno Vallet, Wan Chiu Li, and Bruno Lévy.  N-symmetry direction field design.  ACM Transactions on Graphics (TOG), 27(2):1–13, 2008.  https://doi.org/10.1145/1356682.1356683 doi:10.1145/1356682.1356683.
  * [24] F. C. Stutz, J. P. Groen, O. Sigmund, and J. A. Bærentzen.  Singularity aware de-homogenization for high-resolution topology optimized structures.  Structural and Multidisciplinary Optimization, aug 2020.  https://doi.org/10.1007/s00158-020-02681-6 doi:10.1007/s00158-020-02681-6.
  * [25] Florian Cyril Stutz, Tim Felle Olsen, Jeroen Peter Groen, Niels Aage, Ole Sigmund, Justin Solomon, and Jakob Andreas Bærentzen.  Synthesis of frame field-aligned multi-laminar structures.  ACM Transactions on Graphics, 2022.  https://doi.org/https://doi.org/10.1145/3516522 doi:https://doi.org/10.1145/3516522.
  * [26] Krister Svanberg.  The method of moving asymptotes—a new method for structural optimization.  International journal for numerical methods in engineering, 24(2):359–373, 1987.  https://doi.org/10.1002/nme.1620240207 doi:10.1002/nme.1620240207.
  * [27] Erik Träff, Ole Sigmund, and JP Groen.  Simple single-scale microstructures based on optimal rank-3 laminates.  Structural and Multidisciplinary Optimization, 59(4):1021–1031, 2019.  https://doi.org/10.1007/s00158-018-2180-3 doi:10.1007/s00158-018-2180-3.
  * [28] Fengwen Wang, Boyan Stefanov Lazarov, and Ole Sigmund.  On projection methods, convergence and robust formulations in topology optimization.  Structural and Multidisciplinary Optimization, 43(6):767–784, 2011\.  https://doi.org/10.1007/s00158-010-0602-y doi:10.1007/s00158-010-0602-y.
  * [29] Junpeng Wang, Christoph Neuhauser, Jun Wu, Xifeng Gao, and Rüdiger Westermann.  3D-TSV: The 3D trajectory-based stress visualizer.  Advances in Engineering Software, 170:103144, 2022.  https://doi.org/10.1016/j.advengsoft.2022.103144 doi:10.1016/j.advengsoft.2022.103144.
  * [30] Junpeng Wang, Rüdiger Westermann, and Jun Wu.  A streamline-guided dehomogenization approach for structural design.  Journal of Mechanical Design, 145(2):021702, 2023.  https://doi.org/10.1115/1.4056148 doi:10.1115/1.4056148.
  * [31] Junpeng Wang, Jun Wu, and Rüdiger Westermann.  Stress topology analysis for porous infill optimization.  Structural and Multidisciplinary Optimization, 65(3):1–13, 2022\.  https://doi.org/10.1007/s00158-022-03186-0 doi:10.1007/s00158-022-03186-0.
  * [32] Junpeng Wang, Jun Wu, and Rüdiger Westermann.  Stress trajectory guided structural design and topology optimization.  In International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, volume 1, page 1. American Society of Mechanical Engineers, 2022.  https://doi.org/10.1115/DETC2022-89030 doi:10.1115/DETC2022-89030.
  * [33] Rebekka V. Woldseth, J. Andreas Bærentzen, and Ole Sigmund.  Phasor noise for dehomogenisation in 2d multiscale topology optimisation.  Computer Methods in Applied Mechanics and Engineering, 418:116551, 2024.  https://doi.org/10.1016/j.cma.2023.116551 doi:10.1016/j.cma.2023.116551.
  * [34] Jun Wu, Niels Aage, Rüdiger Westermann, and Ole Sigmund.  Infill optimization for additive manufacturing – approaching bone-like porous structures.  IEEE Transactions on Visualization and Computer Graphics, 24(2):1127–1140, 2018.  https://doi.org/10.1109/TVCG.2017.2655523 doi:10.1109/TVCG.2017.2655523.
  * [35] Jun Wu, Christian Dick, and Rüdiger Westermann.  A system for high-resolution topology optimization.  IEEE transactions on visualization and computer graphics, 22(3):1195–1208, 2015.  https://doi.org/10.1109/TVCG.2015.2502588 doi:10.1109/TVCG.2015.2502588.
  * [36] Jun Wu, Ole Sigmund, and Jeroen P Groen.  Topology optimization of multi-scale structures: a review.  Structural and Multidisciplinary Optimization, pages 1–26, 2021\.  https://doi.org/10.1007/s00158-021-02881-8 doi:10.1007/s00158-021-02881-8.
  * [37] Jun Wu, Weiming Wang, and Xifeng Gao.  Design and optimization of conforming lattice structures.  IEEE Transactions on Visualization and Computer Graphics, 27(1):43–56, 2021.  https://doi.org/10.1109/TVCG.2019.2938946 doi:10.1109/TVCG.2019.2938946.