Micro- and Macro-Scale Topology Optimization of Multi-Material Functionally Graded Lattice Structures

Abstract

Lattice structures are becoming an increasingly attractive design approach for the most diverse engineering applications. This increase in popularity is mainly due to their high specific strength and stiffness, considerable heat dissipation, and relatively light weight, among many other advantages. Additive manufacturing techniques have made it possible to achieve greater flexibility and resolution, enabling more complex and better-performing lattice structures. Unrestricted material unit cell designs are often associated with high computational power and connectivity problems, and highly restricted lattice unit cell designs may not reach the optimal desired properties despite their lower computational cost. This work focuses on increasing the flexibility of a restricted unit cell design while achieving a lower computational cost. It is based on a two-scale concurrent optimization of the lattice structure, which involves simultaneously optimizing the topology at both the macro- and micro-scales to achieve an optimal topology. To ensure a continuous optimization approach, surrogate models are used to define material and geometrical properties. The elasticity tensors for a lattice unit cell are obtained using an energy-based homogenization method combined with voxelization. A multi-variable parameterization of the material unit cell is defined to allow for the synthesis of functionally graded lattice structures.

1. Introduction

Engineers across a variety of industries, including aerospace, automotive, and biomedical, have long sought lightweight yet high-efficiency structures, such as those with high stiffness, high specific strength, or good heat dissipation properties [1,2]. One approach to achieving these properties is through numerical topology optimization, which gained widespread attention following the publication of the seminal work by Bendsøe and Kikuchi in 1988 [3].

It is important to investigate the potential benefits and drawbacks of using different geometries at the micro-scale and evaluate the applicability of the methodology as proposed by Wang et al. [4]. In order to minimize computational complexity and reduce processing time, a voxelization method is suggested as a viable alternative to a highly discretized conformal mesh for constructing the continuous surrogate model that represents the graded evolution of the unit cells.

A lattice structure consists of an arrangement of struts or faces connected by nodes, which when assembled builds up to a geometry within a determined unit cell. Connecting these unit cells results in a complete lattice structure. Lattice structures are comparable to other materials, such as foam or honeycomb types of material. The research community has disseminated works which indicate that lattice structures have better mechanical properties in compression and shear compared to foams and honeycombs [5,6,7]. The enhanced mechanical properties are mainly due to the fact that lattice structures can be tailor-made for each desired application. The unit cells or the struts/faces within the unit cell can have different design variables assigned to them, enabling numerous combinations of the geometry in use and optimization of them to satisfy specific requirements.

Topology optimization is a process used to find the optimal material distribution that satisfies one or more objective functions while adhering to a set of constraints and boundary conditions. The design domain is discretized into several finite elements in a finite element approach and the material distribution is controlled and described by one or more design variables, where the simplest case assumes the density of the finite element to be the only design variable.

Topology optimization can achieve optimal structures by considering geometric patterns that transform spatially over the defined design domain and through multiple length-scales. When a single scale is considered for optimization, it is classified as mono-scale, while when more than one scale is considered it is categorized as multi-scale. Multi-scale approaches can optimize both material micro-structure and structural macro-scale concurrently, providing better results than approaches performed separately and/or successively.

Homogenization theories form the basis of multi-scale approaches, allowing for the computation of effective material properties by considering the material to be made up of repeating patterns of similar unit cells. The review paper by Wu et al. [8] categorizes unit cell topology design approaches into five major categories. The more flexible the design is, the more problems may arise concerning connectivity, integrity, and manufacturing issues.

Multi-scale optimization has been the focus of many research works, including those that employ a single-design variable parameterization. Wang et al. [9] and Watts et al. [10] successfully employed surrogate models to represent the elasticity tensor characterizing the material properties in 2D and 3D unit cell optimizations, respectively. Ypsilantis et al. [11] extended this approach to multimaterial 3D optimization of the micro-scale. Nguyen and Choi explored the effect of increasing design flexibility by optimizing the unit cell build orientation and addressing the concept of equivalent unit cells for conformal lattice designs [12].

In more recent works, multi-variable multi-scale optimization has also been studied. Imediegwu et al. [13] analyzed a multi-variable parameterization of the unit cell using 7 design variables, which has been expanded to 10 design variables in a recent work by the same group [14]. However, the high number of design variables in such models leads to highly complex computational models that may require a dynamic database that updates during the optimization process to reduce computational effort and time.

A previous work by Costa et al. [15] presented a method that enabled (in short) the discretization of unit cell geometries at the micro-scale and the study of its mechanical capabilities at a macro-scale level. Taking the aforementioned works into consideration, the present work focuses on the discretizion and study of different geometries, while also presenting four distinct case studies from the literature to be used in the macro-scale part of the process.

2. Theoretical Background

2.1. Topology Optimization

Topology optimization (TO) is a process that aims to find the optimal distribution of one or more materials within a given design domain while considering a set of constraints, such as volume fraction. An optimization problem can be formulated to minimize the compliance of a material subject to a volume constraint, which can be represented by Equation (1).

In this equation, c refers to compliance, |$\Omega$|$\left(\rho \right)$ represents the volume of material, $\mathit{U}$ is the displacement vector, and F is the global force vector [16].

$$\begin{array}{cc}\hfill \mathrm{Find}& :\rho =\{{\rho }_1,\dots ,{\rho }_N\},\hfill \\ \hfill \mathrm{Minimize}& :c\left(\rho \right)={\mathit{U}}^{\mathit{T}}K\mathit{U}=\sum _{e=1}^N{u}_e^T{k}_e{u}_e,\hfill \\ \hfill \mathrm{subject}\mathrm{to}& :|\Omega |\left(\rho \right)/|\Omega |\le \Lambda ,\hfill \\ \hfill & :K\mathit{U}=\mathit{F},\hfill \\ \hfill & :0\le {\rho }_e\le 1,e=1,\dots ,N.\hfill \end{array}$$

(1)

The element stiffness matrix, $k_e$, can be obtained by using the elasticity tensor, D, which can be determined for isotropic materials.

$$k_e={\int }_{{\Omega }_e}{\mathit{B}}^{\mathit{T}}\mathit{DB}d{\Omega }_e.$$

(2)

To compute the strain-displacement matrix, B, an isoparametric finite element formulation is used, where the physical coordinate system is transformed into a reference coordinate system centered in the element.

The expression in Equation (2) can be reformulated in terms of natural coordinates, as shown in Equation (3), where $|J|$ denotes the determinant of the Jacobian matrix. To accurately calculate the element stiffness matrix, the computation of matrices B and J is critical. Gao et al. [17] explain the process for calculating these matrices in detail. Gaussian quadrature can be used to evaluate the integral in Equation (3), and it is considered more effective than the Newton–Cotes quadrature because it optimizes the position and weights of the sampling points at which the integral is evaluated.

$$k_e={\int }_{-1}^1{\int }_{-1}^1{\int }_{-1}^1{\mathit{B}}^{\mathit{T}}\mathit{DB}|J|d\xi d\eta d\zeta$$

(3)

2.1.1. Sensitivity Analysis

To solve a topology optimization problem, it is necessary to perform a sensitivity analysis, which involves calculating the derivatives of both the objective function, $c\left(\rho \right)$, and the volume constraint, |$\Omega$|$\left(\rho \right)$ (as defined in Equation (1)), with respect to the design variables.

By representing the global stiffness matrix as a sum of N element stiffness matrices, $k_e^G$, the formulation can be written in the global reference frame (Equation (4)), where $u_e$ is the element vector of nodal displacements.

$$\frac{tialc\left(\rho \right)}{tial{\rho }_e}=-u_e{\left(\rho \right)}^T\left(\frac{tial}{tial{\rho }_e}{k}_e\left({\rho }_e\right)\right)u_e\left(\rho \right)$$

(4)

The derivative of the volume constraint, |$\Omega$|$\left(\rho \right)$, with respect to the density variable is given by Equation (5).

$$\frac{tial|\Omega |\left(\rho \right)}{tial{\rho }_e}=|{\Omega }_e|$$

(5)

2.1.2. Filtering

Filtering techniques are commonly employed in the topology optimization (TO) process to mitigate numerical issues, such as checkerboard patterns and mesh dependency. Mesh-independent filtering methods are the most widely used despite other potential solutions, such as perimeter and global gradient control [18].

The modified element density, $x_e$, is based on a neighboring region of that element, determined by the filter radius, $r_{min}$ [19]. The density of the finite element, e, is computed using Equation (6), where $N_e$ represents the set of neighboring elements to which the distance from centre to centre of the elements, $\Delta (e,i)$, is smaller than the filter radius dimension, and $H_{ei}$ is a weight factor defined in Equation (7) [16].

$${\tilde{x}}_e=\frac{1}{{\displaystyle \sum _{i\in N_e}}{H}_{ei}}\sum _{i\in N_e}{H}_{ei}{x}_i$$

(6)

$$H_{ei}=max(0,r_{min}-\Delta (e,i))$$

(7)

However, the application of a density filter results in the loss of physical meaning for the original densities, $x_i$, and the solution must rely on filtered densities. This approach produces filtered sensitivities, $\left(\frac{\widetilde{tial\psi }}{tial{x}_e}\right)$, which can be mathematically expressed using Equation (8). A small value is incorporated into the equation to avoid division by zero when the element densities ${\rho }_e$ are zero.

$$\frac{\widetilde{tial\psi }}{tial{x}_e}=\frac{1}{max(\gamma ,x_e){\displaystyle \sum _{i\in N_e}}{H}_{ei}}\sum _{i\in N_e}{H}_{ei}{x}_i\frac{tialc}{tial{x}_i}$$

(8)

2.2. Material Homogenization

To simulate the mechanical behavior of a lattice material, a homogenization technique may be employed. This approach enables the determination of the effective elastic properties of a material, assuming two conditions hold: first, the macro-scale characteristics are much greater than the micro-scale dimensions (large separation of scales), and second, the micro-structure exhibits a periodic distribution on the macro-scale, composed of periodically repeating unit cells [20].

One can use a finite element approach which solves the equations using the finite element method. The periodic unit cell can be divided into N finite elements to write Equation (9). Superscript i pertains to the six independent test strains in a 3D scenario, while subscript e denotes the element domain. $B_e$ represents the element strain-displacement matrix and $D_e$ is the elasticity tensor for the typically isotropic material that makes up the corresponding finite element.

$$K{\chi }^i=f^i,i=1,\dots ,6$$

(9)

$$K=\sum _{e=1}^N{\int }_{V_e}{\mathit{B}}_{\mathit{e}}^{\mathit{T}}{\mathit{D}}_{\mathit{e}}{\mathit{B}}_{\mathit{e}}d{V}_e$$

(10)

$$f^i=\sum _{e=1}^N{\int }_{V_e}{\mathit{B}}_{\mathit{e}}^{\mathit{T}}{\mathit{D}}_{\mathit{e}}{ϵ}^{0\left(i\right)}d{V}_e$$

(11)

Once the solution for the displacement field, ${\chi }^i$, is established, the term ${\chi }_e^i$ can be determined, and this can be achieved by solving Equation (9).

$$D_{ij}^H=\frac{1}{|V|}\sum _{e=1}^N{\left(u_e^{A\left(i\right)}\right)}^T{k}_e{u}_e^{A\left(j\right)}$$

(12)

If the unit cell is subjected to an imposed strain field, ${ϵ}_{ij}^0$, the displacement field can be expressed, as shown in Equation (13) [21,22]. The fluctuation field term corresponds to a deviation from the linear displacement field that arises due to the heterogeneous structure of the composites.

For a continuous physical body assumption to be valid, two continuity conditions must be met at the boundaries of the repeating unit cells. The first condition is that the displacements must be continuous and the unit cell boundaries cannot overlap or separate from others of adjoined unit cells. Equation (13) ensures this condition.

$${\chi }_i={ϵ}_{ij}^0{y}_j+{\chi }_i^{*}$$

(13)

The second condition is that the traction distributions at the opposite parallel boundaries of a periodic unit cell must be the same. To address this, Equation (13) can be transformed into a set of explicit expressions that consider the displacements on pairs of parallel opposite surfaces of the unit cell. Equations (14) and (15) provide the general formulation for the displacements on a pair of parallel opposite surfaces, with the superscripts $k+$ and $k-$ denoting the opposite parallel surfaces perpendicular to the kth = 1, 2, 3 principal direction (i.e., X, Y, or Z axis).

$${\chi }_i^{k+}={ϵ}_{ij}^0{y}_j^{k+}+{\chi }_i^{*}$$

(14)

$${\chi }_i^{k-}={ϵ}_{ij}^0{y}_j^{k-}+{\chi }_i^{*}$$

(15)

To eliminate the fluctuation term, ${\chi }_i^{*}$, Equation (15) can be subtracted from Equation (14), resulting in Equation (16).

Even though Equation (16) guarantees the first continuity condition, known as the traction continuity, where the normal and shear stresses must be equal at parallel boundaries, it may not always be met. However, for a displacement-based finite element,

$${\chi }_i^{k+}-{\chi }_i^{k-}={ϵ}_{ij}^0(y_j^{k+}-y_j^{k-})={ϵ}_{ij}^0\Delta y_j^k.$$

(16)

2.3. Lattice Structures

Depending on their individual properties, lattice structures can be categorized into two groups: stretch-dominated and bending-dominated lattices.

The choice between these two types of lattices depends on the application’s primary objective. In stretch-dominated lattices, $M\ge 0$, external loads are balanced by axial tension and compression in struts, meaning that no bending occurs at nodes, while in bending-dominated lattices there are too few struts to balance the external forces without balancing the moments induced at the nodes, causing bending stresses to develop in struts [23,24].

$$M=s-3\times j+6.$$

(17)

The Maxwell number measures the stiffness of the lattice structure by taking into account the number of struts, s, and joints, j, that make up the unique geometry of the unit cell (Equation (17)).

2.4. Sensitivity Analysis

In order for the optimizer to update the design variables and achieve convergence for the final solution, it requires access to the derivatives of both the objective function and the volume constraints with respect to those design variables. These sensitivities are crucial for the optimizer to make informed decisions about how to adjust the design variables. Specifically, the sensitivities regarding the objective function can be expressed in mathematical terms to provide a more precise understanding of the relationship between the design variables and the objective function, as expressed in Equation (18):

$$\frac{tialc\left(X\right)}{tial{\tilde{\rho }}_e}=-u_e{\left(X\right)}^T\left(\frac{tial{k}_e({\tilde{\rho }}_e,{\tilde{x}}_{ma{t}_e})}{tial{\tilde{\rho }}_e}\right)u_e\left(X\right).$$

(18)

3. Methodology

3.1. Surrogate Model

The first step involves building a sampling-point database that includes various combinations of design variables, lattice branch radius, and the resulting homogenized elasticity tensor for each combination. The second step involves connecting all sampling points to create the continuous surrogate model, which requires selecting a suitable method to fit the data into a continuous model using regression analysis.

A voxel-based model is chosen to generate the discretization of the micro-scale. The unit cube geometry is divided into structured hexahedral finite elements with eight nodes per element, resulting in a highly regular mesh. This choice is based on the fact that the voxel-based discretization eases the definition of periodic boundary conditions that are pivotal to the adopted homogenization method [13].

To address this issue, Imediegwu et al. [13] proposed a scheme that considers the linear dependency between the material assigned to each voxel and the number of nodes inside the domain delimited by the lattice structure. If this distance is less than or equal to the radius, the node is considered active and counts to set the density to the elements that share it. If it is more than the radius, the node is considered inactive. Research papers such as [10,12,25] are examples of how voxelization has been used to compute effective material properties for different lattice structures.

3.1.1. Unit Cell Parameterization

The two-variable parameterization of the unit cell is accomplished using the relative density and aspect ratio as design variables. The relative density represents the percentage of solid material volume in each unit cell and is defined in Equation (19). The aspect ratio, ${\rho }_{ext}$, is defined as a ratio between the volume of external lattice branches, $V_{ext}$, and the total lattice material volume, $V_{strut}$, as presented in Equation (20).

$$\rho =\frac{V_{strut}}{V_{uc}}$$

(19)

$${\rho }_{ext}=\frac{V_{ext}}{V_{strut}}$$

(20)

Overall, the two-variable parameterization approach follows the rationale presented by Wang et al. [4], with only the variable names changing. The approach enables the concurrent optimization of both the macro- and micro-scale in a multi-scale topology optimization.

3.1.2. Unit Cell Geometries

The geometries used at a micro-scale level in each unit cell play a fundamental part when compiling and creating the study cases in the TO process, in the later stages of the work. As stated previously, six geometries are to be presented and studied at both micro- and macro-scale.

In the literature [5,26,27,28], there might be specific titles for such geometries/branches presented next, although (and to better identify the final geometry combination) in the context of this work these were the titles given. Each geometry is identified as the arbitrary name of the outer branches plus arbitrary name of the inner branches.

The geometry in Figure 1 has a Maxwell number of $-12$ and is therefore bending-dominated.

The geometry in Figure 2 has a Maxwell number of $-1$ and is bending-dominated.

The geometry in Figure 3 has a Maxwell number of $-1$ and is bending-dominated.

The geometry in Figure 4 has a Maxwell number of 0; therefore, it is stretch-dominated.

The geometry illustrated in Figure 5 has a Maxwell number is $-9$, meaning this geometry is bending-dominated.

Lastly, the geometry in Figure 6 has a Maxwell number of $-13$; this geometry is bending-dominated.

3.1.3. Sampling Point Database

The assembling of a complete database is a crucial step in the development of a surrogate model. In this case, the range of relative densities of the sampling points varies from 0 to 1, initially with steps of 0.025 up to rel. density of 0.1, then steps of 0.1 until the value of 0.9 is reached, and finally steps of 0.025 until a value of 1 of relative density. The aspect ratio varies from 0 to 1 as well, with a constant step interval of 0.1. Thus, the database has a total of 165 sampling points that each have a unique combination of the two design variables, relative density and aspect ratio, leading to different lattice branch diameters for each sampling point. Once the sampling-point database script has collected all the required data, the points are connected together during the data-fitting process.

3.1.4. Regression Model

This approach aims to connect all the sampling points and can be visualized in a three-dimensional graph as a surface, considering a two-design variable parameterization of the lattice material unit cell. One of the simplest ways to perform regression analysis is by using polynomial regression. The least squares method (LSM) is a popular approach that follows the concept of parameterized interpolation for lattice material (PILM) introduced in [4]. To reduce the relative error between the sampling data and the regression model, the percentage robust least squares method was chosen.

$$D=D^H(\rho ,{\rho }_{ext})$$

(21)

Using the PILM model, it is possible to represent each of the independent elements of the elasticity tensor (6 for the general anisotropic case and only 3 if the unit cell displays cubic symmetry) as a polynomial function of the design variables.

$$\begin{array}{cc}\hfill D_{ij}& =a_0+a_1{\rho }_{ext}+a_2\rho +\hfill \\ & +a_3{\rho }_{ext}^2+a_4{\rho }_{ext}\rho +a_5{\rho }^2+\hfill \\ & +a_6{\rho }_{ext}^3+a_7{\rho }_{ext}^2\rho +a_8{\rho }_{ext}{\rho }^2+a_9{\rho }^3+\hfill \\ & +a_{10}{\rho }_{ext}^4+a_{11}{\rho }_{ext}^3\rho +a_{12}{\rho }_{ext}^2{\rho }^2+a_{13}{\rho }_{ext}{\rho }^3+a_{14}{\rho }^4+\hfill \\ & +a_{15}{\rho }_{ext}^5+a_{16}{\rho }_{ext}^4\rho +a_{17}{\rho }_{ext}^3{\rho }^2+a_{18}{\rho }_{ext}^2{\rho }^3+a_{19}{\rho }_{ext}{\rho }^4+a_{20}{\rho }^5+\hfill \\ & +a_{21}{\rho }_{ext}^6+a_{22}{\rho }_{ext}^5\rho +a_{23}{\rho }_{ext}^4{\rho }^2+a_{24}{\rho }_{ext}^3{\rho }^3+a_{25}{\rho }_{ext}^2{\rho }^4+a_{26}{\rho }_{ext}{\rho }^5+a_{27}{\rho }^6+\hfill \\ & +a_{28}{\rho }_{ext}^7+a_{29}{\rho }_{ext}^6\rho +a_{30}{\rho }_{ext}^5{\rho }^2+a_{31}{\rho }_{ext}^4{\rho }^3+a_{32}{\rho }_{ext}^3{\rho }^4+a_{33}{\rho }_{ext}^2{\rho }^5+a_{34}{\rho }_{ext}^1{\rho }^6+\hfill \\ & +a_{35}{\rho }^7\hfill \end{array}$$

(22)

By using the polynomial function displayed in Equation (22) to represent the independent elements of the elasticity tensor, the surrogate model can be built by expressing the elasticity tensor of the anisotropic lattice material as being a function of the defined design variables. This allows for an immediate representation of the relationship between the design variables and the material properties of the lattice.

3.2. Multi-Material TO

A modified version of the SIMP method has been developed by Venugopal et al. [29], which connects it with a surrogate model describing the elasticity tensor of each finite element based on its density value.

The product of a SIMP-based equation (defining which material is being allocated) with the surrogate model’s equation results in an equation that translates the stiffness matrix of element e ($k_e$). The material variable is denoted as ${\tilde{x}}_{mat}$ and the density variable as $\tilde{\rho }$, with the tilt indicating that these variables are filtered. Assuming a material with a Young’s modulus of unity, we can express the element stiffness matrix as $k^0$; q represents a penalization factor.

$$k_e({\tilde{\rho }}_e,{\tilde{x}}_{ma{t}_e})=\left[E_1{(1-{\tilde{x}}_{ma{t}_e})}^q+E_2{\tilde{x}}_{ma{t}_e}^q\right]k_e^0\left({\tilde{\rho }}_e\right)$$

(23)

The equations for the volume constraints are presented in Equation (24) (representing the volume constraint of material 1) and Equation (25) (representing the volume constraint of material 2). Here, |${\Omega }_e$| represents the element volume, |$\Omega$| represents the design space volume, and the volume fraction of material i is denoted as $\Lambda$. As a result, the volume constraints are dependent on both design variables.

$$G_1({\tilde{\rho }}_e,{\tilde{x}}_{ma{t}_e})=\sum _{e=1}^N|{\Omega }_e|{\tilde{\rho }}_e(1-{\tilde{x}}_{ma{t}_e})-{\Lambda }_1|\Omega |\le 0$$

(24)

$$G_2({\tilde{\rho }}_e,{\tilde{x}}_{ma{t}_e})=\sum _{e=1}^N|{\Omega }_e|{\tilde{\rho }}_e\left({\tilde{x}}_{ma{t}_e}\right)-{\Lambda }_2|\Omega |\le 0$$

(25)

To minimize the compliance objective functional, the formulation of the multi-material multi-scale TO problem is used:

$$\begin{array}{cc}\hfill \mathrm{Find}& :X=\{{\rho }_1,\dots ,{\rho }_N;{\rho }_{ex{t}_1},\dots ,{\rho }_{ex{t}_N};x_{ma{t}_1},\dots ,x_{ma{t}_N}\},\hfill \\ \hfill \mathrm{Minimize}& :c\left(\tilde{X}\right)=U^{T}KU=\sum _{e=1}^N{u}_e^T{k}_e{u}_e,\hfill \\ \hfill \mathrm{subject}\mathrm{to}& :\sum _{e=1}^N|{\Omega }_e|{\tilde{\rho }}_e(1-{\tilde{x}}_{ma{t}_e})-{\Lambda }_1|\Omega |\le 0,\hfill \\ \hfill & :\sum _{e=1}^N|{\Omega }_e|{\tilde{\rho }}_e\left({\tilde{x}}_{ma{t}_e}\right)-{\Lambda }_2|\Omega |\le 0,\hfill \\ \hfill & :KU=F,\hfill \\ \hfill & :0\le {\rho }_{min}\le {\rho }_e\le {\rho }_{max}\le 1,e=1,\dots ,N,\hfill \\ \hfill & :0\le x_{ma{t}_e}\le 1,e=1,\dots ,N,\hfill \\ \hfill & :0\le {\rho }_{ex{t}_e}\le 1,e=1,\dots ,N.\hfill \end{array}$$

(26)

Furthermore, and considering that for the multi-material multi-scale TO three design variables (in this case) are needed, we are able to write the elasticity tensor, D, as it is presented in Equation (27). The term $D^0({\tilde{\rho }}_{ext},\tilde{\rho })$ represents the surrogate model created for the unit cell’s geometry to be considered for a unitary Young’s modulus.

$$D({\tilde{x}}_{mat},{\tilde{\rho }}_{ext},\tilde{\rho })=[E_1{(1-{\tilde{x}}_{mat})}^q+E_2{\tilde{x}}_{mat}^q]D^0({\tilde{\rho }}_{ext},\tilde{\rho })$$

(27)

Sensitivity Analysis

The sensitivity analysis to be performed shares an identical path to the one presented previously in Section 2.4, although now a third design variable has to be taken into consideration: the aspect ratio, ${\rho }_{ext}$. As stated, Equations (28) and (29) express the additional derivatives of the volume constraints with regards to ${\rho }_{ext}$, to properly conclude this step.

$$\frac{tial{G}_1({\tilde{x}}_{mat},{\tilde{\rho }}_{ext},\tilde{\rho })}{tial{\tilde{\rho }}_{ext}}=0$$

(28)

$$\frac{tial{G}_2({\tilde{x}}_{mat},{\tilde{\rho }}_{ext},\tilde{\rho })}{tial{\tilde{\rho }}_{ext}}=0$$

(29)

Additionally, the volume constraint’s sensitivities are as presented in the following equations:

$$\frac{tial{G}_1({\tilde{\rho }}_e,{\tilde{x}}_{ma{t}_e})}{tial{\tilde{\rho }}_e}=(1-{\tilde{x}}_{ma{t}_e})|{\Omega }_e|,$$

(30)

$$\frac{tial{G}_2({\tilde{\rho }}_e,{\tilde{x}}_{ma{t}_e})}{tial{\tilde{\rho }}_e}={\tilde{x}}_{ma{t}_e}|{\Omega }_e|.$$

(31)

4. Results

In this chapter, the six geometries presented in Section 3.1.2 will be used in four case studies. This totals 24 case studies. The surrogate models created for the tensor elements ($D_{11}$, $D_{12}$, $D_{44}$) of each geometry are illustrated in 3D graphics, showing the relative density ($\rho$) and the aspect ratio (${\rho }_{ext}$) design variables.

4.1. Multi-Material Multi-Scale TO

For all topology optimization case studies, the value of the force applied was of 60 kN. The overall volume fraction is also common to all cases and it was set at 50% of the total design-domain volume, 25% being allocated to material 1 and the other 25% to material 2. Both design variable intervals were set between $[0,1]$.

For this work, material 1 was set with a Young’s modulus of $E_1=145$ MPa and material 2 was set with a Young modulus of $E_2=210$ MPa.

4.1.1. Cantilever Case Study

All the TO results of the six geometries for the Cantilever case study are presented. This structure is fully constrained at one end (to the left of the figure) and free to move at the opposite end (to the right of the figure). It has a 12 × 12 finite element-sized square section with 30 unit cells of length.

The Cantilever case study (Figure 7a) starts as a solid bar with a square section. One end of the bar is fully constrained, with neither linear nor rotational degrees of freedom, whereas the other end is free to translate and rotate. A vertical force, pointing downwards, is applied to the middle of the horizontal bottom edge of the free end of the bar. These properties should translate to a TO result that is symmetric in relation to a plane parallel to $xy$.

The renders of the following results of this case study are presented in Figure A1, Figure A2, Figure A3, Figure A4, Figure A5 and Figure A6.

Figure 8a, Figure 9a, Figure 10a, Figure 11a, Figure 12a and Figure 13a illustrate the density as well as the material distribution within the volume domain. On the constrained area (to the left of the illustration), there is a significant concentration of material 2 at the top and the bottom faces along the whole width of the structure. The material tends to shape similarly to a nozzle when we move away from the constraint end towards the free end and culminates on the area where the force is applied. Towards the free end of the cantilever structure, material 1 was allocated (less stiff) near the area where the force is applied. The qualitative density as well as the material distribution is identical to all geometries in this case study.

4.1.2. Messerchimitt–Bolkow–Blohm (MBB) Case Study

The renders of the following results of this case study are presented in Figure A7, Figure A8, Figure A9, Figure A10, Figure A11 and Figure A12.

In this case study (Figure 7b), a solid bar with a square section is horizontally supported in the four bottom edges of the solid. On these four points, the structure is fully constrained, meaning there is no linear or rotational movement. A downwards vertical force is applied to the centre of the upper face of the bar. This case study can perhaps be regarded as a bridge, and the TO should result in a final structure that is symmetrical to a plane parallel to $xy$ as well as to $yz$.

From the figures that illustrate the density distribution within the design domain (Figure 14a, Figure 15a, Figure 16a, Figure 17a, Figure 18a and Figure 19a), we can observe that there is a distribution of material 2 (in red) in the shape of an arc, connecting the constrained elements, on the bottom face, to the middle of the top face of the structure. This arc-shaped density distribution of material 2 (stiffer material) is identical (from a qualitative point of view) in all the MBB case studies of the different geometries.

The material and density distribution of material 1 (in blue) in the structure, on the other hand, presents slight differences in term of “filling” the inner part of the arc previously described with material 2, which is worth highlighting.

4.1.3. Hook Case Study

The renders of the following results of this case study are presented in Figure A13, Figure A14, Figure A15, Figure A16, Figure A17 and Figure A18.

The Hook case study starts with a structure in the L shape with a square section. The top square section is fully constrained, whereas the other square section (other end of the “L” shaped structure) is free to move. A downward vertical force is to be applied on the centre of the free-to-move section. This case study TO should result in a final structure that is symmetrical in relation to a plane parallel to $xy$.

The case studies that are presented below are quite similar in terms of overall relative density distribution and present qualitatively identical areas where material 1 and material 2 was allocated. In the figures that relate to relative density distribution (Figure 20a, Figure 21a, Figure 22a, Figure 23a, Figure 24a and Figure 25a), the salient observation that can be made is that there is a concentration of material 1 (in blue, less stiff) on the corner area of the L-shaped structure, which then switches into material 2 (in red, stiffer) as we move from the corner towards the constrained and free edges. The exact positioning of material does not appear to be the same for the different geometries presented, although the qualitative position of it within the structure is, in a way, identical.

4.1.4. Wheel (“Table Top”) Case Study

The renders of the following results of this case study are presented in Figure 26, Figure 27, Figure 28, Figure 29, Figure 30 and Figure 31.

The Wheel case study features a parallel piped structure that is supported in the four bottom corners. These corners are fully constrained. A downward vertical force is applied to the center of the bottom face of the structure. This case study might resemble a table top. The TO final structure should be symmetrical to a plane parallel to $xy$ and $yz$.

In the Wheel case studies, the qualitative similarities of density as well as material distribution amongst the various case studies with different geometries can be observed. On the bottom corners of the design domain, all figures show a concentration of material 2 (in red) starting from the support points, where the structure is fully constrained. From this point upwards and inwards, the “legs” of the final structure develop towards the middle of the design domain and eventually change into material 1 (in blue).

4.2. Results Comparison

The compliance values at which the TO has converged is displayed in bar-type graphs. In order to check the efficiency of combining the outer and inner branches of the geometries proposed, the converged compliance values for three different scenarios are presented for each geometry: when only the outer branches are used (${\rho }_{ext}=1$, in orange), for when only the inner branches are used (${\rho }_{ext}=0$, in blue), and finally for when the aspect ratio varies freely within the 0 to 1 scope (${\rho }_{ext}\in [0,1]$, in grey), resulting in a combination of the outer and inner branches.

An immediate observation that can be made from the information presented in Figure 32 and Figure 33 is that the use of the full geometry presented, meaning having a combination of the outer and inner branches of the geometry, tends to give a slightly or even a substantially better result (compliance value wise) when compared to the cases where the aspect ratio was set to 0 or 1, either only inner branches or only outer branches.

Figure 32a,b illustrate the compliance results obtained in the Cantilever case study for the six geometries. From the values of the grey bars (representing the 0–1 aspect ratio interval), the geometry with the lowest value of compliance was the Diagonal Cross + Inner Cross geometry. On the other hand, the geometry with the greatest compliance value was the Diagonal Cross + Inner Cage (for the scenario where the aspect ratio varies from 0 to 1 only).

Regarding the results obtained for the MBB case study (Figure 32c,d), the geometry which resulted in the lowest value of compliance was the Diagonal Cross + Inner Cage. On the other hand, the geometry with the highest value of compliance was for both the Cube + Inner Diagonal Cross and the Diagonal Cross + Inner Diagonal Cross geometries. The final compliance results are within $3.5\%$ from each other (for the 0–1 aspect ratio scenario).

In the figures presented next (Figure 33a,b), it can be verified that from the Cube + Inner Diagonal Cross geometry resulted the lowest value of compliance for the Hook case study. The Diagonal Cross + Inner Diagonal Cross geometry had the highest value of compliance (for the 0–1 aspect ratio scenario). As a side note, both of these geometries share the same inner branches architecture, although they have different outer branches.

Regarding the Wheel case study, given the information from Figure 33c,d the geometry which resulted in the lowest value of compliance was the Cross + Diamond geometry. The geometry with the highest value of compliance was the Diagonal Cross + Inner Cross geometry (for the 0–1 aspect ratio scenario). The final compliance values for the 0–1 aspect ratio scenario of all 6 geometries vary roughly $1\%$ from one another.

5. Concluding Remarks

The novelties of this work stem from the creation of six surrogate models that effectively represent the elasticity tensor of each of the six geometries concerning two design variables: relative density and aspect ratio. These design variables determine both the volume occupied by the geometry within the designated finite element volume and the distribution of that volume between the inner and outer branches of the geometry.

For the topology optimization process, surrogate models were created for six hierarchical lattice materials during the first solution phase. The creation of these models involved constructing a database based on material homogenization theory. Here, it should be noted that considerable computational time was needed for each run. It was assumed that the number of finite elements would adequately represent the strut-based lattice shape and that the lattice-independent branches would not significantly intersect for the entire density range. While this assumption did not prevent the creation of a surrogate model, it could result in failure to converge for a multi-scale topology optimization. Therefore, it is essential to consider branch interception when using the PILM principle to develop new surrogate models. A set of sampling points for each elasticity tensor were calculated and later interpolated using LSM with a 7th order polynomial, which resulted in the coefficients of the best fitting surface to those points.

Next, multi-scale, multi-variable topology optimization was performed. The overall average improvement considering the compliance result when the $[0,1]$ interval of aspect ratio was considered, and the second best result (either with the result from having aspect ratio set as 0 or 1) was $3.32\%$. In more detail, in the Cantilever case study the average improvement in compliance value, in the same circumstances as previously described, was $1.9\%$. In the MBB case study, the improvement was $1.67\%$. In the Wheel case study, the improvement was $1.7\%$. Lastly, the improvement in the Hook case study was $8.02\%$. These improvements bring to light the advantages of the combination of the outer and inner branches presented into more complex although better-performing geometries.

Furthermore, the amount of finite elements seems to have been sufficient, given that none of the TO runs resulted in a divergent behaviour. The number of iterations enabled the TO process to run until an acceptable level of convergence had been reached.