Design and Optimization of Conforming Lattice Structures

TL;DR

This paper introduces a new method for designing conforming lattice structures that optimize shape and internal material distribution to maximize stiffness, applicable to both 2D and 3D domains, with demonstrated structural benefits.

Contribution

It presents a novel two-step approach for shape and lattice material optimization, including conforming lattice extraction, which is robust and applicable to complex geometries.

Findings

Enhanced structural stiffness in optimized lattice designs

Effective conforming lattice structures for 2D and 3D domains

Physical verification confirms numerical results

Abstract

Inspired by natural cellular materials such as trabecular bone, lattice structures have been developed as a new type of lightweight material. In this paper we present a novel method to design lattice structures that conform with both the principal stress directions and the boundary of the optimized shape. Our method consists of two major steps: the first optimizes concurrently the shape (including its topology) and the distribution of orthotropic lattice materials inside the shape to maximize stiffness under application-specific external loads; the second takes the optimized configuration (i.e. locally-defined orientation, porosity, and anisotropy) of lattice materials from the previous step, and extracts a globally consistent lattice structure by field-aligned parameterization. Our approach is robust and works for both 2D planar and 3D volumetric domains. Numerical results and physical verifications demonstrate remarkable structural properties of conforming lattice structures generated by our method.

Tables

Table 1. TABLE I : The difference in compliance predicted by homogenization and a full resolution analysis, for the lattice structures shown in Fig.  13 .

	a	b	c	d	e	f
Homo.	418.33	282.62	239.97	332.81	277.27	232.64
Full res.	444.78	300.15	255.48	323.18	292.66	241.94
Diff.	6.32%	6.20%	6.46%	2.89%	5.55%	4.00%

Table 2. TABLE II : Statistics of 3D model complexity and computational performance. The timing is reported in minutes.

Model	Resolution	#Ele.	Vol.	#It.	$J_{com}$	$T_{FEA}$	$T_{Opt}$	#vertex	#tet	#strut	$T_{pre}$	$T_{posy}$	$T_{extr}$	$T_{Total}$
3D cantilever (Fig. 8a)	$100\times 50\times 50$	250k	0.2	60	110.84	3.11	0.85	0.89m	5.0m	48k	1.80	5.28	0.33	11.36
3D cantilever (Fig. 8b)	$100\times 50\times 50$	250k	0.2	60	96.03	3.97	1.62	1.65m	9.62m	87k	1.75	9.98	0.70	18.03
3D cantilever (Fig. 8c)	$100\times 50\times 50$	250k	0.2	60	85.85	5.83	2.65	1.62m	9.41m	25k	1.83	6.50	0.57	17.38
Bridge (Fig. 8d)	$200\times 38\times 88$	644k	0.1	60	230.52	15.13	1.88	1.18m	6.54m	63k	1.41	7.09	0.47	25.97
Bridge (Fig. 8e)	$200\times 38\times 88$	644k	0.1	60	177.86	16.76	3.80	2.02m	11.57m	111k	2.84	13.58	1.14	38.12
Bridge (Fig. 8f)	$200\times 38\times 88$	644k	0.1	60	149.96	21.10	6.29	1.89m	10.77m	35k	2.40	8.29	0.72	38.79
chair (Fig. 1)	$140\times 100\times 200$	1.8m	0.1	60	193.5	30.92	5.03	3.32m	18.60m	178k	4.15	18.66	1.87	60.63
femur (Fig. 9)	$140\times 93\times 182$	696k	0.5	6	163.4	0.99	0	5.86m	14.26m	305k	12.36	35.50	5.94	54.79
dragon (Fig. 9)	$200\times 90\times 143$	461k	0.5	6	99.4	1.12	0	4.09m	23.31m	200k	5.09	24.84	2.88	33.92
3D cantilever (Fig. 14)	$100\times 50\times 50$	250k	0.2	60	110.84	3.11	0.85	6.65m	38.50m	351k	7.19	33.44	6.25	50.84
Bridge (Fig. 14)	$200\times 38\times 88$	644k	0.1	60	230.52	15.13	1.88	8.64m	49.63m	462k	12.50	56.56	14.35	100.42

Equations

$$v_e({\alpha }_e)=1-\frac{({\alpha }_{e,x}l-2t)({\alpha }_{e,y}l-2t)}{{\alpha }_{e,x}{\alpha }_{e,y}{l}^2}.$$

$${\rho }_e({\phi }_e,{\alpha }_e)={\phi }_e{v}_e({\alpha }_e).$$

$$\underset{\mathbf{\phi },\mathbf{\alpha },\mathbf{R}}{\min }$$

$$s.t.:$$

$${\phi }_e\in [0.0,1.0],\forall e$$

$${\alpha }_{e,k}\in [{\underline{\alpha }}_k,{\overline{\alpha }}_k],k\in \{x,y,z\},\forall e.$$

$$\mathbf{K}(\mathbf{\phi },\mathbf{\alpha },\mathbf{R})\mathbf{U}=\mathbf{F}.$$

$${\mathbf{K}}_e={\int }_{{\Omega }_e}{B}^{\mathsf{T}}{D}_{e}Bdx,$$

$${\mathbf{K}}_e(1,{\alpha }_e,R_e)={\int }_{{\Omega }_e}{B}^{\mathsf{T}}{D}_e({\alpha }_e,R_e)Bdx,$$

$$D_e({\alpha }_e,R_e)={\overline{R}}_e(R_e)D_e({\alpha }_e){\overline{R}}_e^{\mathsf{T}}(R_e),$$

$${\mathbf{K}}_e({\phi }_e,{\alpha }_e,R_e)={\phi }_e^p{\mathbf{K}}_e(1,{\alpha }_e,R_e),$$

$$E(\mathbf{P})=\sum _{i\in V}\sum _{j\in N(i)}||{\mathbf{p}}_i-({\mathbf{M}}_{ij}{\mathbf{t}}_{ij}+{\mathbf{p}}_j)|{|}^2,$$

$${\mathbf{M}}_{ij}=\frac{({\mathbf{R}}_i+{\mathbf{R}}_{j}r({\mathbf{R}}_i,{\mathbf{R}}_j))}{||{\mathbf{R}}_i+{\mathbf{R}}_{j}r({\mathbf{R}}_i,{\mathbf{R}}_j)||}\cdot \frac{({\mathbf{S}}_i+{\mathbf{S}}_j)}{2},$$

$${\mathbf{t}}_{ij}=round[{\mathbf{M}}_{ij}^{-1}({\mathbf{p}}_i-{\mathbf{p}}_j)].$$

$$R=\left(\begin{array}{ccc}{l}_1& l_2& l_3\\ m_1& m_2& m_3\\ n_1& n_2& n_3\end{array}\right),$$

$$\overline{R}=\left(\begin{array}{cc}A& B\\ C& D\end{array}\right),$$

$$A=\left(\begin{array}{ccc}{l}_1^2& m_1^2& n_1^2\\ l_2^2& m_2^2& n_2^2\\ l_3^2& m_3^2& n_3^2\end{array}\right),$$

$$B=\left(\begin{array}{ccc}2m_1{n}_1& 2n_1{l}_1& 2l_1{m}_1\\ 2m_2{n}_2& 2n_2{l}_2& 2l_2{m}_2\\ 2m_3{n}_3& 2n_3{l}_3& 2l_3{m}_3\end{array}\right),$$

$$C=\left(\begin{array}{ccc}{l}_2{l}_3& m_2{m}_3& n_2{n}_3\\ l_3{l}_1& m_3{m}_1& n_3{n}_1\\ l_1{l}_2& m_1{m}_2& n_1{n}_2\end{array}\right),$$

$$D=\left(\begin{array}{ccc}{m}_2{n}_3+m_3{n}_2& n_2{l}_3+n_3{l}_2& m_2{l}_3+m_3{l}_2\\ m_3{n}_1+m_1{n}_3& n_3{l}_1+n_1{l}_3& m_3{l}_1+m_1{l}_3\\ m_1{n}_2+m_2{n}_1& n_1{l}_2+n_2{l}_1& m_1{l}_2+m_2{l}_1\end{array}\right).$$

Full text

1 Introduction

The design of lightweight structures by optimization is a classical and still active topic in engineering. Stimulated by the increasingly high flexibility and resolution offered by 3D printing, there has been a growing interest in optimizing structures that are composed of delicate microstructures [1, 2]. These approaches assume that the microstructures are aligned with a prescribed regular grid. This simplifies modelling, simulation and optimization. It, however, limits the solution space, and thus the achievable structural performance. The microstructures are typically anisotropic (e.g. a hollowed cubic cell with uniform thickness is stiffer along its axes than along its diagonals). It is known that the orientation of anisotropic materials in stiffness-optimal structure coincides with the principal stress directions resulting from forces acting on these materials [3]. Furthermore, axis-aligned microstructures do not match the curved surfaces of 3D shapes.

To address the above issues, in this paper we propose an efficient and robust method to generate conforming lattice structures. A lattice is a connected array of struts. The lattice structure generated by our method is conforming in two aspects: the struts align with principal stress directions, maximizing structural stiffness; and, struts on the boundary capture the curved surface of the optimized shape. We note that the shape, according to design options accessible to the user, is allowed to evolve together with the optimization of lattice distribution, i.e. the optimized shape is a subset of the design domain.

Our method has two major steps. In the first step, both the shape and the spatially-varying orientation of lattices inside the shape evolve simultaneously according to stress analysis and numerical optimization. Rather than relying on extremely high-resolution finite elements to capture the evolving lattice geometric details, we develop a homogenization-based topology optimization method which allows to efficiently simulate and optimize the lattice material distribution on a relative coarse level. By introducing a novel parameterization of the unit cell, our method ensures a uniform thickness of struts while allowing a sufficient degree in lattice anisotropy. The second step, which we call lattice compilation, extracts a globally consistent lattice structure from the optimized, locally-defined lattice configuration, including orientation, porosity, and anisotropy. We address the challenging problem of extracting connected lattices across cells with spatially-varying orientation, by extending field-aligned meshing techniques. This extension allows a fast and robust lattice compilation where anisotropic geometric features are incorporated.

The specific contributions of our paper include:

• •

A novel workflow for designing, in both 2D and 3D, conforming lattice structures based on homogenization-based topology optimization and field-aligned parameterization.

• •

Insights into optimal lattice structures from a mechanical perspective, analyzed via a detailed parameter study.

• •

A simple and effective parameterization of the unit cell for allowing structural anisotropy while ensuring a uniform thickness of struts.

• •

A new formulation to allow simultaneous optimization of the shape and the lattice distribution.

• •

A novel approach for extracting globally consistent lattice structures that accommodate anisotropy and heterogeneity.

Our method generates highly detailed lattice structures. The optimized lattice chair in Figure 1, for instance, consists of $178,291$ struts, achieved on a simulation resolution of $140\times 100\times 200$.

The remainder of the paper is organized as follows. In the next section we review related work. In Section 3 we give an overview of the proposed method. In Sections 4 and 5 the two major steps of our method, lattice optimization and compilation, are presented. Results and analysis are presented in Section 6, before the conclusions are given in Section 7.

2 Related Work

2.1 Structural Optimization for 3D Printing

In the era of 3D printing (and more broadly, digital fabrication), structural optimization becomes increasingly relevant in computational design [4]. Skin-frame structures [5], honeycomb-like Voronoi structures [6], tree-like supporting structures [7], and bone-inspired porous structures [8] have been optimized as lightweight infill for prescribed 3D shapes. Guided by outputs from structural optimization, Martinez et al. proposed to use graded orthotropic foams as a parameterized metamaterial to fill a prescribed shape [9, 10]. In contrast to design and optimize internal structures for prescribed shapes, our method optimizes concurrently the shape and its internal microstructures for application-specific loads. Different from two-scale structural optimizations (e.g. [11, 1]) which assume axis-aligned microstructures, our method optimizes the orientation of microstructures, in particular, to align it with spatially-varying stress directions. We restrict our design method to lightweight microstructures that are composed of struts, i.e. lattice structures.

Lattice structures are typically aligned with a regular grid [12, 13]. Rosen and his co-authors proposed a method to design lattice structures that conform with the boundary surface of a prescribed 3D shape [14, 15]. Our method optimizes concurrently the shape and align the lattices with stress directions. The alignment of structures along principal stress directions improves structural performance [16, 3]. This principle has been applied to 2D planes (e.g. [17]) and curved surfaces [18, 19, 20]. The appealing 2D results are achieved by tracing stress directions or based on a ground structure approach [21, 22]. Due to their inherent challenges associated with the initialization of samples/nodes, a direction extension of these methods to 3D volumetric lattices is not applicable. Our method constructs stress-aligned 3D volumetric lattices, relying on homogenization-based topology optimization and field-aligned meshing.

Our method is among recent efforts to structural analysis and optimization for 3D printing. Stava et al. proposed a method to detect and correct structural defects [23]. Recent efforts include worst-case structural analysis [24, 25], and stochastic structural analysis [26]. Chen et al. proposed a solver for inverse elastic shape design [27]. Ulu et al. optimized structures under force location uncertainty [28]. Our method, targeting on stiffness maximization of lattice structures under certain static loads, is complementary to these efforts. Yet the integration goes beyond the scope of this paper.

2.2 Homogenization-based Topology Optimization

Topology optimization is an important design method for 3D printing, as it effectively leverages the fabrication flexibility to create structures with exceptional (mechanical) properties. Topology optimization transforms optimal shape design as a material distribution problem. In their seminal work, Bendsøe and Kikuchi proposed a homogenization method, which optimizes the distribution of square unit cells with variable rectangular holes [29]. Due to the lack of manufacturing means for such microstructures back then, the homogenization method was replaced by density-based approaches (e.g. SIMP [30]) which have since been widely used in industry and in many academic contributions (e.g. large scale optimization [31, 32, 33]).

In light of the capability of 3D printing to fabricate microstructures, the homogenization method was recently revisited to design structures with manufacturable microstructures [34, 35, 36], based on the rectangular hole model [29]. A challenge is to compile a continuous structure from hollowed cells that are defined on a regular grid, and that, after optimization, have different orientations. To this end, a projection approach proposed by Pantz and Trabelsi [37] was improved to connect the orthotropic microstructures [34, 38]. The output structure is represented by high-resolution pixels or voxels.

Our approach to design conforming lattice structures follows the homogenization approaches, but differs in three aspects. First, we propose a new parameterization of cells to ensure that the variable cells have a constant thickness, while allowing a large degree in anisotropy. Uniform thickness is a common practice in (metal) additive manufacturing of (axis-aligned) lattices [12, 13]. Such structures can also be fabricated by a direct extrusion in 3D [39], or by robotic fabrication [40, 41]. We note that variable thickness is not impossible with 3D printing. However, accurately achieving variable thickness requires precise control of the material flow rate in 3D printing hardware and robust algorithms in slicing [42]. Ensuring a uniform thickness alleviates these requirements. Second, our method simultaneously optimizes the lattice distribution and the shape, achieved using multiple design variables. Last but not least, while existing works exploited projection methods for generating high-resolution pixel or voxel models, we develop a novel approach based on field-aligned meshing to compile the lattice structure. The optimized structure is compactly represented by a graph. This direction shares a similar goal with the recent work by Arora et al. [43]. In contrast to the design approach [43], our method unlocks a large solution space by optimizing the porosity, anisotropy, and orientation of lattices.

2.3 Field-aligned Parameterization

We develop a lattice compilation method based on field-aligned parameterization which has been researched intensively in the past decade, especially for generating quadrilateral (quad-) mesh. We review briefly the more recent works on hexahedral (hex-) meshing, and for quad-meshing we refer an interested reader to the survey by Bommes et al. [44].

For a given 3D closed shape, field-aligned hex/hex-dominant meshing techniques typically consist of three steps [45, 46, 47, 48, 49, 50]. It starts by estimating the gradients of a volumetric parameterization using a directional field [51, 52], where the field is discretized per vertex or per tetrahedron and smoothly interpolated within the volume under a boundary alignment constraint. This is followed by computing a parameterization aligned with the estimated gradients by fitting. Finally it extracts the hex-mesh in the parametric space [53]. Robust hex-meshing remains a challenging problem. A promising direction is to topologically correct the directional fields [47, 48, 54, 55]. Lei et al. introduced a hex-mesh generation method based on surface foliation theory [56]. This approach, however, requires heavy topological pre-processing of the input.

The field-aligned parameterization pipeline is primarily used for generating semi-regular meshes. To ensure the validity of the mesh, complex geometric and topological computations are involved. In this paper we make use of field-aligned parameterization to generate lattice structures. This new application differs from mesh generation, as lattices are encoded by graphs rather than meshes. This goal sidesteps the numerical stability issue and geometrical and topological complexities typically occurred during mesh extraction from the parameterization.

To efficiently and robustly extract consistent lattice structures, we extend the robust meshing approach that was proposed by Jacob et al. [57] and further developed by Gao et al. [50]. The per vertex local parameterization from [50] fits our purpose well since the local parameterization aligns exactly with the direction field by construction and permits fast and scalable computations. The extension proposed in this paper allows to incorporate anisotropy and heterogeneity.

The recent work by Arora et al [43] shares a similar goal as ours, i.e. to extract field aligned struts from stress directions. Our approach takes the optimized stress fields as input, without a field smoothing operation that compromises the accuracy of input fields. During lattice compilation, while they extract the struts by tracing stress lines and simplifying the duplicated ones, our approach directly generates struts by simple and efficient graph operations. This makes our approach fast and scalable, taking a couple of minutes for an input with tens of millions of tetrahedral elements (see Table II).

3 Overview

Given a design domain and application-specific loads, our method generates a lattice structure that maximizes structural stiffness. The struts in the optimized lattice structure conform with principal stress directions. Moreover, the struts on the boundary span a smooth surface faithfully approximating the optimized shape.

As illustrated in Figure 2 for 2D and Figure 1 for 3D, our approach consists of two steps. The first optimizes the shape (including its topology) and the distribution of lattice material within the shape. The input includes a design domain and boundary conditions (Figure 1 left and Figure 2a), as well as design specifications such as the target fraction of solid material. The design domain in 3D is represented by a closed triangle surface mesh. This mesh is voxelized, generating finite elements for simulation and optimization. The output is a set of fields, indicating, per element, the occupancy of lattice material, and the orientation and anisotropy of lattice material (Figure 1 second left and Figure 2b). A surface mesh is then reconstructed, representing the optimized shape, i.e. the interface between elements that are filled with lattice material and that are empty. The shape enclosed by this surface mesh is tetrahedralized. The optimized fields are then interpolated on the vertices of the tetrahedral model. The vertices, including their connectivity and their associated field values, serve as input for the second step, which compiles a globally connected lattice structure that are composed of struts (Figure 1 second right and Figure 2c). The output lattice structure is encoded by a graph.

4 Lattice Optimization

The goal of our optimization is to find the optimal distribution of lattice material that maximizes structural stiffness, subject to a number of design constraints. To this end, the design space is discretized into a regular grid of bilinear square elements in 2D or trilinear cubic elements in 3D. As illustrated in Figure 3 for a 2D rectangular design domain, each element is to be filled by repeating a unique, rectangular-shaped cell. The cells are adapted from a unit cell by scaling and rotation. The scaling factors and rotation matrices are to be optimized. The scaling factors for the cell in element $e$ are denoted by $\alpha_{e}$, and in 2D by $(\alpha_{e,x},\alpha_{e,y})$ and in 3D by $(\alpha_{e,x},\alpha_{e,y},\alpha_{e,z})$. The rotation matrix for element $e$ is denoted by $R_{e}$.

The unit cell is a hollowed square with a side length of $l$ and a thickness of $t$, which are specified by the user. To ensure a uniform thickness of struts across the design domain, during scaling the side length of cells is elongated, while the thickness ($t$) is kept constant. This creates cells with gradation in the fraction of solid material ($v_{e}$),

[TABLE]

This gradation allows the optimization to place adapted cells with a smaller fraction of solid material in regions where the stress is relatively small. Furthermore, per axis elongation allows to increase the mechanical anisotropy of cells. This is beneficial since the stress tensors are typically anisotropic.

Besides a scaling factor per axis and a rotation matrix, each element is assigned a variable $\varphi_{e}$, to indicate whether the element is empty ($\varphi_{e}=0$) or filled ($\varphi_{e}=1$) with lattice material. The set of elements that are filled with lattice material defines the overall shape of the optimized structure. To allow for gradient-based numerical optimization, the variable $\varphi_{e}$ is relaxed to take intermediate values, i.e. $\varphi_{e}\in[0,1]$. This variable is akin to the density variable in classical density-based topology optimization, which in that context indicates the fraction of solid material. In the context of lattice optimization, this variable shall be interpreted as the fraction of lattice material. The fraction of solid material per element ($\rho_{e}$) depends on $\varphi_{e}$ and the fraction of solid material within an adapted cell ($v_{e}$), i.e.

[TABLE]

As the design space is parameterized by the fraction of lattice material ($\bm{\varphi}$), scaling factor ($\bm{\alpha}$), and orientation matrix ($\bm{R}$), the optimization problem is given as

[TABLE]

Here the objective is to minimize the work done by the external force, which is equivalent to minimize compliance (i.e. stiffness maximization). $\mathbf{F}$ denotes the force vector that is applied to the design domain. The force vector is constant. $\mathbf{U}$ denotes the displacement vector of the shape when it comes to its static equilibrium under the external force $\mathbf{F}$. The first constraint, Eq. 3b, restricts the amount of solid material, where $\overline{v}$ is the fraction of available solid material, and $N$ is the number of elements in the design domain. The second constraint, Eq. 3c, sets bounds for the fraction of lattice material ($\bm{\varphi}$). The third constraint, Eq. 3d, sets bounds for the scaling factors ($\bm{\alpha_{x}},\bm{\alpha_{y}},\bm{\alpha_{z}}$). The lower and upper bounds of the scaling factors are user-defined.

The novelty of this formulation is two-folds. First, by optimizing the scaling factors rather than the thickness of hollowed cells, it ensures that all struts in the optimized structure have the same thickness. As discussed in Section 2.2, this eases the control of the 3D printing process. Second, we assign an additional variable $\varphi$ to indicate the occupation of lattice material. This makes the formulation more general. Prescribing $\varphi=1$ leads to optimized lattices that fill the entire design domain. This is useful as infill for prescribed shapes. Allowing $\varphi$ to be decided by the optimization enables both the shape and the lattice to evolve simultaneously, achieving a higher stiffness.

4.1 Stiffness Matrix for Lattices

The objective function, Eq. 3a, involves the displacement vector ($\mathbf{U}$), which is related to the external force ($\mathbf{F}$). The unknown $\mathbf{U}$ is computed by solving the equilibrium equation with finite element analysis,

[TABLE]

Here the stiffness matrix, $\mathbf{K}$, is assembled from per element stiffness matrix, $\mathbf{K}_{e}(\varphi_{e},\alpha_{e},R_{e})$.

In standard finite element analysis of solids [58], the element stiffness matrix $\mathbf{K_{e}}$ is computed by integrating over the domain of an element, $\Omega_{e}$,

[TABLE]

where $B$ is the element strain-displacement matrix for linear basis functions [58]. $D_{e}$ is the fourth order elasticity tensor, computed based on the Young’s modulus and Poisson’s ration of the solid material.

For analyzing elements that are filled with lattice material, the elasticity tensor $D_{e}$ is not constant but rather depends on design variables $\alpha_{e}$ and $R_{e}$. Let us first consider an element that is filled lattice with $\varphi_{e}=1$. The stiffness matrix for lattices is calculated by

[TABLE]

The elasticity tensor of a rotated lattice cell, $D_{e}(\alpha_{e},R_{e})$, is computed by rotating the elasticity tensor of this cell in its local coordinate system, $D_{e}(\alpha_{e})$. In engineering notation, $D_{e}$ is represented as a ${3\times 3}$ matrix for 2D problems or a ${6\times 6}$ matrix for 3D. The rotation of tensor is realized by

[TABLE]

where the tensor rotation matrix $\overline{R}$ is given in the Appendix.

The effective elasticity tensor of an elongated cell, $D_{e}(\alpha_{e})$, is evaluated by numerical homogenization. We make use of the Matlab code provided in [59] and  [60] for homogenization in 2D and 3D, respectively. Given the scaling factors, the domain of the elongated unit cell is discretized by square finite elements with linear basis functions. To avoid performing homogenization for every $\alpha_{e}$ during the optimization process, we pre-compute $D_{e}$ for regularly sampled $\alpha$ values. In 2D, we fit a surface for every non-zero entry in $D$ over the 2D domain of $[\underline{\alpha}_{x},\overline{\alpha}_{x}]\times[\underline{\alpha}_{y},\overline{\alpha}_{y}]$. In 3D we use trilinear interpolation.

For elements with $\varphi_{e}$ between [math] and $1$, we use the power law from density-based approaches [30] to interpolate,

[TABLE]

where the parameter $p$ (typically $p=3$) is introduced to penalize intermediate values in $\varphi_{e}$, and consequently the optimization steers $\varphi_{e}$ towards either [math] or $1$.

4.2 Solving

The optimization problem (Eq. 3) is solved in an iterative process, as in [34]. In each iteration the following computational steps are performed, until the maximum change in design variables is smaller than a threshold (or the maximum number of iterations is reached).

First, the equilibrium equation (Eq. 4) is solved, obtaining the displacement vector, $\mathbf{U}$. From the element displacement vector ($\mathbf{U}_{e}$), strain ($\overline{\epsilon}_{e}$) and stress ($\overline{\sigma}_{e}$) per element, in engineering notation, are calculated by $\overline{\epsilon}_{e}={B}\mathbf{U}_{e}$ and $\overline{\sigma}_{e}=D_{e}(\alpha_{e},R_{e})\epsilon_{e}$, respectively.

Second, design variables $\varphi$ and $\alpha$ are updated using a gradient-based solver. We make use of the method of moving asymptotes (MMA) [61]. To avoid checkerboard patterns, the design variables are regulated into $\tilde{\varphi}$ and $\tilde{\alpha}$ using the so-called density filter. $\tilde{\varphi}$ is then projected into $\overline{\tilde{\varphi}}$ by a (smoothed) Heaviside operation, to approach a 0-1 solution. The filter and Heaviside operator are widely used in density-based approaches, e.g. in e.g. [62, 8].

Third, the orientation of each element ($R_{e}$) is updated based on the associated stress tensor ($\sigma_{e}$). The stress tensor is symmetric positive-definite. By eigendecomposition we obtain three mutually orthogonal principal stress directions ($v_{1},v_{2},v_{3}$). The eigenvectors are ordered by respective eigenvalues in ascending order, i.e. $\gamma_{1}\leq\gamma_{2}\leq\gamma_{3}$. As shown by Pedersen [3], the optimal orientation of an orthotropic material coincides with the principal stress directions, hence the element is rotated by $R_{e}=[v_{1}^{\mathsf{T}};v_{2}^{\mathsf{T}};v_{3}^{\mathsf{T}}]$.

Fourth, the stiffness matrices of lattices, $K_{e}$, are re-calculated based on the updated orientation ($R_{e}$) and regulated variables $\overline{\tilde{\varphi}}$ and $\tilde{\alpha}$, according to Section 4.1.

4.3 Example

The output of our optimization is a set of fields defined on the design domain. Figure 4 visualizes these fields for a rectangular domain, which is discretized by a grid of $80\times 40$ elements. The unit cell has $l=10t$. The maximum fraction of solid material is $0.15$. Figure 4a shows the optimized lattice fraction field. Even with a fraction of solid material as small as $0.15$, the lattice covers a large portion of the design space. This is due to the fact that the unit cell has a small fraction of solid material (i.e. $36\%$, with $l=10t$). Figure 4b visualizes the orientation of optimized cells. Here the rotated frame is elongated according to the respective scaling factor per axis. For clarity the frames are shown for regularly-spaced samples. On the right, the frame field is visualized for elements which have a fraction of lattice material ($\varphi_{e}$) that is larger than a threshold ($0.5$).

5 Lattice Compilation

Up to this step, we have equipped with a design domain with a set of fields including fraction of lattice infill, orientation, and scaling, that are optimized for the prescribed external loads (cf Figure 4). Since a region with a low fraction implies that few materials is required, we extract a sub-area (volume) from the design space by thresholding ($\geq 0.5$) out low infill regions. With the actual shape being extracted, we now focus on generating a lattice structure that conforms to both the boundary of the shape and the directional and scaling fields.

Our problem setting differs from the typical meshing problem in that both of our input and output are quite relaxed from the conditions of being a mesh. For the input, we put no constraints on its geometrical quality (i.e. angles, edge ratios, etc) nor its topological correctness (i.e. manifoldness, no holes, and free of intersections). This maximizes the scope of the problem but poses a great challenge to the design of a robust solution. For the lattice output, it does not require face (solid) elements, making complex topological operations in most of the meshing methods unnecessary for our purpose. Moreover, considering that it is not a hard requirement for our lattice structure to be all-hex cells for the designed structure to function, we choose the parameterization optimization in [50] that can be easily adapted to handle graphs and propose a simple extraction strategy to generate a lattice structure. The produced lattice structure contains mostly quad (hex)-like connections while allowing certain irregularity to adapt for rapid changes in the directions and scales.

Our method takes a graph with vertices of the optimized shape as the input: $G=(V,E)$, where every vertex $\mathbf{v}\in V$ has a position $\mathbf{x}\in\mathbb{R}^{k}$ ($k$ is 2 for 2D and 3 for 3D), an orientation matrix $\mathbf{R}\in\mathbb{R}^{k\times k}$ encoding the cross directions and also denoting a local coordinate system, and a scaling vector $\mathbf{s}\in\mathbb{R}^{k}$ (i.e. $\bm{\alpha}$ in the previous section) composed of scales for the $k$ axes of the local coordinate system. Our goal is to extract a lattice structure, which is another graph $G^{\prime}=(V^{\prime},E^{\prime})$ that (1) reproduces the input direction and the anisotropic property as much as possible, and (2) has a resolution that can be flexibly controlled by a target edge length $h$.

In the following, we first describe the parameterization optimization that incorporates anisotropic orientations, and then present the lattice structure extraction.

5.1 Parameterization

Given an orientation field $\mathbf{O}$ that includes the cross directions for all the vertices, we want to compute a parameterization $\mathbf{P}$ with the gradient aligned to $\mathbf{O}$. Methods that compute a global parameterization with the gradient aligning to the orientation field in a least-square sense (e.g. [63, 64, 43]), rely on non-linear optimization solvers which are not scalable to large dataset. We instead compute a parameterization for the input graph by representing it with a set of local parameterization and minimizing an energy between the local parameterizations of adjacent vertices [57, 50].

The local nature of the parameterization makes it easily parallelizable and scalable to large inputs.

As illustrated in the inset, the local parameterization of a vertex in 2D plane (or 3D volume) can be uniquely determined by its origin $\mathbf{p}$, the orientation matrix $\mathbf{R}$, and unit lengths $h\cdotp\mathbf{s}$, where $h$ is the user-defined global target edge length. The unit lengths are fixed through the entire process. Unlike the previous approaches [50, 57] that treat directions as a 4 rotational symmetric field in 2D or a 24 rotational symmetric field in 3D, since the unit length varies for different axis, our coordinate system is mutable only by flipping the signs of each axis while maintaining the right-hand rule. The origin with a random initialization, is the variable we need to optimize.

Given the above setting, the optimization energy of the parameterization $\mathbf{P}$ is defined as the summation of all the squared differences of local parameterizations for each edge:

[TABLE]

where $N(i)$ is a set of all the vertices sharing an edge with vertex $i$, $\mathbf{M}_{ij}$ is an interpolation of $\mathbf{M}_{i}$ and $\mathbf{M}_{j}$ where $\mathbf{M}=\mathbf{R}\mathbf{S}$ and $\mathbf{S}$ is the scaling matrix converted from $h\cdotp\mathbf{s}$, and $\mathbf{t}_{ij}\in\mathbb{Z}^{k}$ encodes the integer translations of $\mathbf{p}_{j}$. $\mathbf{M}_{ij}\mathbf{t}_{ij}+\mathbf{p}_{j}$ translates $\mathbf{p}_{j}$ by integer moves to the nearest position to $\mathbf{p}_{i}$, effectively avoiding the integer jumps between the two local parameterization and only the difference of their fractional parts is measured. The computation of $\mathbf{M}_{ij}$ is to interpolate the directions and scales separately,

[TABLE]

where $r(a,b)$ is the closest matching that gives the smallest difference between two coordinate systems which can be computed efficiently by enumerating all the cases. There are only two cases to compare in 2D and six cases for 3D.

The integer translation between two connecting vertices in the parameterization space, $\mathbf{t}_{ij}$, is computed by a rounding operation,

[TABLE]

By doing so, the energy between the two vertices will be minimized.

We minimize $E(\mathbf{P})$ in a Gauss-Seidel style by iteratively visiting every vertex and smoothing the origin of each vertex by computing an average of all its neighbors. The pseudo code of the optimization is provided in Algorithm 1.

The last step in line 8 rounds each origin of a local parameterization $\mathbf{p}_{i}$ to the integer position closest to the vertex position $\mathbf{x}_{i}$. Consequently, each component of $\mathbf{t}_{ij}$ becomes -1, 0, or 1. For example, as illustrated in Fig. 5, $\mathbf{t}_{ij}=(0,0)$ for Fig. 5 left, $\mathbf{t}_{ij}=(\pm 1,0)$ or $(0,\pm 1)$ for Fig. 5 middle, and $\mathbf{t}_{ij}=(\pm 1,\pm 1)$ for Fig. 5 right.

To speed up the optimization, similar to [50], we construct a hierarchical structure of the input graph by halving the number of vertices for each level and perform the optimization on each level of the hierarchy by 50 iterations for 2D and 200 iterations for 3D.

5.2 Graph Extraction

In the input graph $G$, each vertex $\mathbf{v}$ has a smoothed local parameterization. The origin $\mathbf{p}$ of $\mathbf{v}$ provides a guidance for the vertex position in the output graph $G^{\prime}=(V^{\prime},E^{\prime})$. Besides, the integer translation associated with each edge ($\mathbf{v}_{i}\mathbf{v}_{j}$) of $G$, $\mathbf{t}_{ij}\in\mathbb{Z}^{k}(k=2,3)$, categorizes this edge as a specific element in $G^{\prime}$, depending on the $L_{0}$ norm of $\mathbf{t}_{ij}$ which is the number of $\pm 1$s in $\mathbf{t}_{ij}$. In 3D ($k=3$), this number can be

• •

0 (i.e. $\mathbf{t}_{ij}=(0,0,0)$), indicating that the two vertices are very close in the parameterization space, and thus will be collapsed into to a point in $G^{\prime}$,

• •

1 ($\mathbf{t}_{ij}=(\pm 1,0,0),(0,\pm 1,0),\textrm{or}\,(0,0,\pm 1)$), meaning that the edge is parallel to one of the stress directions, and thus will be kept in $G^{\prime}$,

• •

2 ($\mathbf{t}_{ij}=(\pm 1,\pm 1,0),(\pm 1,0,\pm 1),\textrm{or}\,(0,\pm 1,\pm 1)$) or 3 ($\mathbf{t}_{ij}=(\pm 1,\pm 1,\pm 1)$), respectively corresponding to a rectangular or cuboid diagonal, which deviates from the stress directions and thus shall not appear in $G^{\prime}$.

For example, black and dashed green edges in Figure 6 left correspond to $||\mathbf{t}_{ij}||_{0}=$ 1 and 2, respectively.

By utilizing the positional guidance of $\mathbf{p}$ and the indication of $\mathbf{t}_{ij}$, the graph extraction is straightforward: collapse the edges with $||\mathbf{t}_{ij}||_{0}=0$ (dots in Figure 6 represent the averaged positions of collapsed edges), keep the edges with $||\mathbf{t}_{ij}||_{0}=1$, and remove the diagonal edges (i.e. $||\mathbf{t}_{ij}||_{0}=2\,\textrm{or}\,3$).

While the above procedure generates a graph with mostly right angles, we notice T-junctions in the final graph with near flat angles that are suboptimal for the stiffness of the lattice structure. Figure 6 left illustrates a vertex with T-junction in 2D. This can be attributed to the fact that the removal of the diagonal edges is aggressive. The T-junctions appear near singularities of the parameterization (similar to the positional singularities in [50]) which result in elements with non-right angles, for example, triangles and pentagons in 2D, and prisms and general polyhedra in 3D.

To address this issue, we propose to keep some diagonal edges in the final graph. Specifically, right after collapsing edges with $||\mathbf{t}_{ij}||_{0}=0$, we check the configuration of every vertex in the graph and identify critical diagonals. As illustrated in Figure 6, for a vertex in black, the process is done by first normalizing all of its adjacent edge vectors onto a unit circle (sphere in 3D), then computing their nearest directions over 4 rotational-symmetric ones in 2D (6 in 3D), e.g. red and dark green arrows, and finally relabelling a diagonal edge to be an edge with $||\mathbf{t}_{ij}||_{0}=1$ such that each of the 4 (6 in 3D) stress directions is represented by an edge that is close to the direction (Figure 6 right).

In summary, the process to extract the graph $G^{\prime}$, i.e. a lattice structure, consists of the following steps.

Categorize the edges in $G$ based on $||\mathbf{t}_{ij}||_{0}$. 2. 2.

Group vertices in $G$ according to $||\mathbf{t}_{ij}||_{0}$ such that groups are connected by edges with $||\mathbf{t}_{ij}||_{0}\neq 0$. Note that a group might contain only a single vertex. 3. 3.

Generate the initial $G^{\prime}$. For each group, a new vertex is positioned at the average of the origins of vertices in $G$. This vertex inherits the edges to new vertices that are converted from neighbouring groups. 4. 4.

Categorize the edges in $G^{\prime}$ based on $||\mathbf{t}_{ij}||_{0}$. 5. 5.

Identify and relabel critical diagonals in $G^{\prime}$ to avoid T-junctions, and remove other diagonal edges.

6 Results

6.1 Examples

Our method works for both 2D and 3D. Figure 8 shows three optimized 2D lattice structures. In (a), the lattice distributes across the prescribed curved shape, with spatial variations in orientation, porosity, and anisotropy. In (b) and (c), the optimized lattices cover a subset of a rectangular design domain, with variations in orientation in (c) and additionally variations in porosity and anisotropy in (b). The unit cell in 2D is specified with $l=10t$, $\underline{\alpha}=1$, and $\overline{\alpha}=4$.

Figure 8 show 3D lattice structures optimized by our method. Our method is also applicable to design lattices that spread over a prescribed 3D curved shape, as shown in Figure 9. 3D examples are optimized with a unit cell using $l=4t$, $\underline{\alpha}=1$, and $\overline{\alpha}=2$. Figure 10 shows fabricated bridge and cantilevers. The models are 3D printed with Form 2 which uses stereolithography (SLA). The dimension of models and the thickness of struts are scaled to comply with the volume and feature size of the printer. The printed femur (Figure 11) has a dimension of $112.4\times 77.9\times 133.1\,mm^{3}$, with a thickness of $0.5\,mm$. The chair (Figure 1) is $110.8\times 76.6\times 142.1\,mm^{3}$, with a thickness of $0.6\,mm$.

6.2 Evaluation

6.2.1 Design options

We evaluate the influence of various design options on the resulting lattice structures using the 2D cantilever ( Section 4.3, Figure 4), with the fraction of solid material $\overline{v}=0.15$, and bounds for the scaling factors $\underline{\alpha}=1$ and $\overline{\alpha}=4$. The optimized fields and compiled lattice structures are shown in Figures 13 and 13, respectively.

In the first row of Figure 13, the fraction of lattice is fixed, $\varphi=1$. Consequently the lattice distributes across the entire rectangular design domain. In (a) the scaling is also fixed, while in (b) the optimization of scaling is enabled. The enlarged solution space leads a decrease in compliance (i.e. improved stiffness), $418.33$ (a) vs. $282.62$ (b). In (c), the scaling factors along individual axes are decoupled, resulting in a further decrease in compliance to $239.97$.

In the second row of Figure 13, the fraction of lattice is optimized. Consequently, a shape evolves from the optimization, corresponding to $\varphi_{e}\geq 0.5$. Similar to the trend of compliance in the first row, it decreases from (d), to (e), and to (f), along with the increased flexibility in optimization. (f) has the smallest compliance among the six cases. It decreases from (a) by $44.39\%$. This confirms the significance of adaptive porosity and anisotropic features for stiffness maximization.

As a reference, an axis-aligned uniform lattice structure covering the entire domain (i.e. corresponding to the initialization of Figure 13) is evaluated. Its compliance is $852.30$, which is more than twice larger than the design in Figure 13a, and $3.66$ times larger than the design in Figure 13f. This comparison confirms the importance of aligning anisotropic microstructures with internal stress directions for stiffness maximization.

6.2.2 Accuracy

To evaluate the accuracy of our lattice compilation method, we perform a comparison of the compliance predicted by homogenization with the compliance of lattice structures by a full-resolution finite element analysis. To this end, the six lattice structures in Figure 13 are discretized by a finite element resolution of $4096\times 2048$, and analyzed using a geometric multigrid elasticity solver [65]. The comparison is summarized in Table I. The difference in compliance is between $2.89\%$ and $6.46\%$. This demonstrates that our lattice compilation introduces little error to the predicted performance from homogenization-based optimization. We note that homogenization theory assumes infinite periodicity of the cells, while for fabrication the compiled lattice has a finite physical size.

6.2.3 Computational performance

Table II presents statistics of 3D model complexity and computational performance. The experiments were run on a standard desktop PC equipped with an Intel Xeon E5-1650 v3 processor (12 cores) running at 3.5 GHz, 64 GB RAM, and an Nvidia GTX1080 graphics card with 8 GB memory. The optimization and compilation together take less than 1 hour even for complex models such as the chair and femur.

The group of columns 2-8 is related to lattice optimization. From the cantilever and bridge examples, it can be observed that with increasing design flexibility the compliance ($J_{com}$) decreases. This agrees with the 2D tests in Figure 13. The increased design flexibility is also reflected by an increase of time associated with updating stiffness matrices, which is counted in $T_{FEA}$. The optimization time of the gradient-based solver for $\bm{\varphi}$ and $\bm{\alpha}$, $T_{Opt}$, increases accordingly as well.

The resolution of optimized fields is refined by a regular subdivision (1 element $\to$ $2^{3}$ elements), followed by tri-linear interpolation of the fields. While our lattice compilation algorithm takes a general graph as the input, in our implementation, we use triangle meshes and tetrahedral meshes which are purely for the convenience of computing vertex normal. This step costs 45$\sim$70 seconds (cantilevers, Figure 8) to 4 minutes and 26 seconds (chair, Figure 1). The refinement generates a large number of vertices (#vertex) organized as tetrahedral elements (#tet), supplied to the lattice compilation. The compiled lattice has as many as 305k struts, for the femur model. Timings for pre-processing, i.e. building data structures ($T_{pre}$), local parameterization ($T_{posy}$), and graph extraction ($T_{extr}$) are reported.

In the last two rows, the optimized fields are refined twice (1 element $\to$ $4^{3}$ elements). This creates highly detailed lattice structures as shown in Figure 14.

6.3 Comparison and Discussion

Comparison with solid structures from density-based topology optimization [30]  A 2D simply supported beam is optimized using our method and the classic density-based approach – Solid Isotropic Material with Penalization (SIMP). The lattice and solid structure generated by our method and SIMP, are shown in Figure 15 a) and b), respectively. The physical sizes are $294.8\times 74.2\times 60\;mm^{3}$, and the struts have an in-plane thickness of $0.8\;mm$, which is twice the nozzle size. They were fabricated by a Ultimaker S5 printer using flexible TPU material. While the digital models were designed using the same fraction of solid material, with 3D printing the lattice structure is heavier ($52$ gram vs. $46$ gram) due to the delicate tool-path.

The load condition of the 3D printed specimen is shown in Figure 15c. It is supported at the two ends on the bottom, while a downward force is applied on the top middle. To avoid out-of-plane buckling of these thin specimens, two wooden plates (with open square windows for observation) are placed to clamp the specimen (Figure 15d). Clamping plates are placed with a gap of $62\;mm$, slightly larger than the thickness of the specimen.

The force-displacement plots for multiple tests are shown in Figure 15e. The forces on the solid structure increase steeper than on the porous structure, meaning that the solid structure from SIMP has a higher stiffness. However, the forces on the solid structure turn down after they reach a peak of about $62$ N. This is due to the (in-plane) buckling of the compressed bars. In contrast, the lattice structure can support a maximum force that is twice larger before it buckles. This is due to the increased effective cross-section area of the substrctures. This test, in agreement with previous physical tests on 3D printed isotropic infill [66], confirms the significance of lattice structures for buckling stability. We note that directly accounting for buckling stability in topology optimization is a much complicated problem, due to the less intuitive definition of the buckling mechanism and demanding eigenvalue problems [67]. Optimizing structures with lattice materials provides an efficient and effective solution to increase buckling stability.

Comparison with bone-like porous structures [8]  Wu et al. proposed a density-based approach to design bone-like porous structures using constraints on local material volume [8]. Figure 16 compares the porous structure and the lattice structure, generated with the same boundary conditions (see Fig. 2a) and the same fraction of solid material. The porous structure was optimized with a local volume fraction of $0.36$, leading to a total volume fraction of $0.288$. We then run lattice optimization with this total volume fraction, with the design options of rotation and scaling. The bone-inspired infill was optimized with a finite element resolution of $400\times 200$, while the lattice was obtained with a simulation resolution of $80\times 40$.

The convergence in compliance is plotted in Figure 17. The compliance of bone-like infill and conforming lattice is $184.64$ and $177.29$, respectively, meaning that the lattice structure is stiffer. Lattice optimization converges much faster, and since it runs on a lower resolution, this leads to a significant speed up. The optimization of lattice took 1 minute and 7 seconds (60 iterations), while the optimization of bone-like infill cost 40 minutes (1000 iterations).

Discussion of Arora et al. [43]  Arora et al. proposed a method to construct a graph with its edges aligned with stress directions from simulation of the solid shape [43]. This approach does not explicitly optimize the orientation, porosity, nor anisotropy. It relates to the option in our method with fixed $\varphi$ and $\alpha$ (cf. Figures 13a and 13a). The result in Figure 13f has demonstrated that with optimized porosity and anisotropy, the compliance reduces by $44.39\%$. We note that under the option of fixed $\varphi$ and $\alpha$, after aligning the lattice, our method re-calculates stress directions and update the lattice orientation. This leads a minor but noticeable decrease in compliance than aligning the lattice with stress directions from the solid shape (420.58 $\to$ 418.90).

Our lattice compilation approach is scalable, for example, the number of struts is more than two orders of magnitude compared with examples shown in [43]. This allows to design highly detailed lattice structures. Figure 14 shows optimized lattice structures with 462k struts (bridge) and 351k struts (cantilever).

7 Conclusions

In this paper we have presented a novel method to design conforming lattice structures by extending homogenization-based topology optimization and field-aligned parameterization. It can compute not only an optimized lattice structure that occupies certain subregions of regular design domains but also lattices that spread over prescribed (curved) shapes. The optimized lattice structures conform with principal stress directions and the boundary of the (optimized) shape. Our method is scalable and allows to optimize highly detailed lattice structures, which can be fabricated by 3D printing. Numerical analysis on different design options confirms the importance of aligning anisotropic lattice with internal stress directions and the importance of lattice gradation in porosity and anisotropy. The compiled lattice structure, by full resolution finite element analysis, has a compliance very close to the compliance predicted by homogenization-based optimization. By physical tests we demonstrate that the optimized lattice structure can support a buckling load twice as large as topology optimized solid structures, at the price of a slight decrease in stiffness. Besides quantified structural performance, the optimized conforming lattice structures look remarkably appealing.

Future work  Our method generates lattice structures particularly optimized for mechanical properties. It provides options to steer the optimization by configuring the design variables, and to adapt the output graph resolution in lattice compilation. It is found that the generated lattice, in certain areas where the stress is small, is less regular. This is attributed to the rapid spatial variation in the underlying stress directions. In some applications the designer might wish to trade the mechanical performance for regularity also in such rapidly change areas. An elegant solution for this balancing, ideally embedded in the lattice optimization step, remains an open question.

Acknowledgements

The authors gratefully acknowledge the support from the LEaDing Fellows Programme at the Delft University of Technology, which has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 707404. Weiming Wang wishes to thank the Natural Science Foundation of China (No. 61702079).

Appendix A Appendix

Denoting a $3\times 3$ rotation matrix by

[TABLE]

the $6\times 6$ rotation matrix for the elasticity tensor in engineering notation is written as

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]