# To Refine or Not to Refine: Topology Optimization of Adaptively Refined Infill Structures for Additive Manufacturing

###### Abstract

We present a novel method for optimizing structures by adaptively refining the structural details. The adaptively refined structures topologically resemble the graph composed of edges in the quadtree mesh. However, while adaptive mesh refinement is employed in numerical analysis for reducing computational complexity, in this work we interpret the edges as structural elements carrying mechanical loads. The adaptivity and full coverage over the design domain make adaptively refined structures well-suited as infill for 3D printed parts, where uniform infill structures have been typically used in practice.

The topology optimization of adaptively refined structures is realized based on two novel ideas. First, we relax the binary design variables (i.e., to refine or not to refine) by using continuous variables, giving rise to sensitivity analysis for an efficient gradient-based optimization. Second, we propose a refinement filter to encode the dependence of design variables among multiple levels in the structural hierarchy. The refinement filter thus enables considering the design variables on all levels simultaneously in the optimization. Our numerical results demonstrate optimized structures spanning multiple levels in the quadtree, with cell sizes smoothly varying in the closed-walled design domain.

###### keywords:

Topology optimization, adaptive structures, additive manufacturing, infill^{†}

^{†}journal:

## 1 Introduction

### 1.1 Background: Topology Optimization for Additive Manufacturing

Topology optimization is a promising design method for additive manufacturing, as it fully exploits the manufacturing flexibility enabled by the additive process. It distributes a limited amount of material in the design domain with the objective of maximizing the structural performance Bendsoe04 (). Reviews of topology optimization techniques can be found in Sigmund13 (); Deaton14 (). Additive manufacturing (AM, also known as 3D printing) provides an effective means to fabricate such optimal structures, which are often too complex to be produced by conventional manufacturing such as machining. Nevertheless, AM also poses a few specific constraints and features which call for innovation in topology optimization Brackett2011 (); Zegard2016SMO ().

The focus of most research works in topology optimization targeting AM is on manufacturability constraints. Building upon density-based topology optimization, Langelaar Langelaar16 (), Qian Qian2017 (), and Gaynor and Guest Gaynor16 () proposed methods to ensure that optimized structures satisfy a constraint on the overhang angle, and thus the usage of additional support structures is reduced or completely eliminated. The overhang constraint has also been addressed in level set-based topology optimization Mirzendehdel2016CAD (); Allaire2017 (), and topology optimization using moving morphable components Guo2017CMAME () or using a virtual skeleton Amir2016ECCOMAS (). Of interest to AM is also the constraint on length scale Guest04 (); Guest09 (). A thorough examination of minimum and maximum length scale is presented in Lazarov16 (). Among others, Zhou et al. Zhou2015CMAME () presented geometric constraints to ensure the minimum length scale, while Zhang et al. Zhang2016CMAME () developed the minimum length scale based on moving morphable components. Lazarov and Wang Lazarov2017CMAME () proposed methods to impose maximum length scale in density-based topology optimization.

Closely related to our current work is research considering infill structures in topology optimization for AM. Rather than printing solid structures, engineering practices in additive manufacturing seem to favour porous structures Gibson10 (); Gao15Survey (). For instance, in fused deposition modelling (FDM), the interior of 3D models is often filled with repetitive triangular or rectangular grids known as infill. See Figure 1 (left) for an illustration of uniform, rectangular infill pattern within a 2D bracket-shaped domain. The concept of porous infill is initially introduced to roughly balance cost and mechanical properties: A larger infill volume percentage leads to a stronger print at the cost of more material and longer printing time. Recent research found out that infill structures can obtain significantly increased stability with respect to buckling Clausen16 () and unpredicted loading conditions Wu17-infill () at the expense of a minor decrease in stiffness. However, these works have assumed a constant infill volume percentage across the design domain: In Clausen16 (); Clausen15 () the base material is interpreted as uniform, isotropic infill, and in Wu17-infill (); Wu17-shell-infill () anisotropic infill is generated yet the effective local volume percentage is constant.

### 1.2 Overview of Adaptive Infill Optimization

In this paper we propose a topology optimization method to generate adaptively refined infill structures. An optimized quadtree for the bracket shape is shown in Figure 1 (right). Adaptively refined structures are a good option for 3D printed infill. On the one hand, the adaptivity can be exploited to enhance regions where denser structures are beneficial to support the mechanical loads. For instance, in the bracket example the optimized quadtree infill performs four times stiffer than the uniform infill. On the other hand, the structures on a uniform coarse grid filling the part’s interior make the printed part robust for uncertain and smaller loads, such as forces during transportation and manipulation that possibly happen to any surface region of the printed model.

Doing topology optimization by adaptive structure refinement allows for an intuitive control over some structural features. For instance, the coarsest and finest grids effectively determine the maximum and minimum void sizes, respectively. The interior structures have a uniform thickness, which simplifies tool-path generation in 3D printing. The uniform thickness and the minimum void size have positive implications on thermal effects in the additive manufacturing process Gibson10 (). Going beyond additive manufacturing, structures with a certain level of regularity facilitate the production of large engineering structures by assembling a few basic sub-structures. In this case, the final quadtree can be interpreted as being assembled from frames. The frames have only a few variations in size, and thus each variation can be massively produced.

In topology optimization of adaptive structures we consider structural refinement (i.e., to refine or not to refine) as design variables. This type of design variables was suggested in Wu16Rhombic (), where rhombic cells to selectively refined based on a heuristic criterion. Starting with a coarse grid, it refines the grid cells in a greedy manner: In each iteration, the cells with a high strain energy density are subdivided. As cells are refined new edges are added into the structure. This leads to an increase in the structural stiffness. The process is terminated once the material volume reaches the prescribed budget. This refinement strategy can be enhanced by combing selective coarsening, which revokes cells that are low in strain energy density. However, accurately solving discrete optimization problem is generally difficult. The greedy approach results in a local optimal solution, and it might be away from the global optimum.

In this work we address the adaptive structure refinement problem using gradient-based continuous optimization. To do so, we solve two challenging issues. First, for each cell, to refine it or not to refine it essentially represents a binary design variable. To enable efficient continuous optimization, we relax the binary design variables to continuous values between - not to refine and - to refine. Following projection approaches in density-based topology optimization, a Heaviside projection Guest04 (); Wang10 () is embedded within the optimization in order to improve the convergence of obtaining discrete 0-1 values.

The continuous representation of refinements also gives rise to a solution to the second challenge — Recursive refinement implies that the design variable of refinement for a fine cell is not activated until its parent coarse cell has been refined. A consequence is that during the optimization iterations the number of design variables is not constant: The coarse cells that have already been subdivided are removed from the list of subdivision candidates, and meanwhile some new subdivision candidates corresponding to newly created fine cells come into action. Using the continuous design variables, we propose a method to simultaneously take into account the refinements on all levels. The dependency of the refinement of a fine cell on refinements of cells on successively coarser levels is encoded by a novel operator, which we call a refinement filter. Consequently, this results in a constant number of design variables in the whole optimization process, and thus facilitates the implementation. More importantly, considering the design variables on all levels altogether enlarges the solution space, compared with the greedy approach where in each iteration only a part of the design variables is activated.

It shall be noted that obtaining spatially adaptive structures clearly distinguishes from previous use of adaptive mesh refinement in topology optimization, where the motivation is to reduce the intensive finite element analysis by reducing the number of elements. Spatial adaptive meshes (i.e., quadtree in 2D and octree in 3D) are often employed in numerical analysis to attain a required accuracy for a minimum amount of computation Berger84 (). In the context of topology optimization, Maute and Ramm Maute95 () are among the first to employ adaptive techniques to decrease the number of design variables and to generate smooth structures. Some recent developments along this direction include adaptive polygonal elements by Nguyen-Xuan Nguyen-Xuan2017 (), and error control in adaptive meshing by Wang et al. Wang14ATO () and Lambe and Czekanski Lambe2017 (). The adaptive mesh used in the literature represents the finite elements for elasticity analysis, and its refinement criterion is based on estimated errors. The obtained structures are similar, if not identical, to the structures optimized with uniform fine meshes. In contrast, here the adaptive mesh, its edges in particular, represents the physical structures to be obtained from optimization. The structural refinement respects global measures on the structural stiffness and the material volume. The structures our method is aiming at cover the entire closed-walled design domain (see Figure 1 right).

## 2 Formulation

This section starts with the discretization of the design domain, including the design variables representing mesh refinements, the finite elements for linear elasticity, and the mapping from design variables to finite elements (Section 2.1). We then proceed to the mathematical formulation of the optimization problem (Section 2.2). The proposed refinement filter is presented in Section 2.3, and the consequent sensitivity analysis is given in Section 2.4.

### 2.1 Hybrid Discretization

We maintain a hybrid discretization. As illustrated in Figure 2, it is composed of a quadtree grid (a) which encodes the topology information, and a discretization by uniform square elements (b) upon which elasticity analysis is performed. On the left, the red lines form the coarse grid where refinements start from. The coarse cells (red) are mapped to blocks of square elements in the finite element grid on the right. In this example, each coarse cell corresponds to a block of elements. Within each block, elements in the first boundary layer are assigned solid (i.e., a density value of ), indicated by black colour. Adjacent coarse cells share a common edge, and this edge in the finite element grid has a thickness of two elements — one element on each side.

For illustrative purpose two coarse cells in the bottom right of the domain are refined, creating two cross-shaped sub-structures (blue). One quadrant is further refined, indicated by the green cross. The crosses are mapped to solid elements on the right, with different grey intensities indicating the different levels of refinement. The mapping from the two-level refinements to finite elements is highlighted in (c) and (d). Here for clarity the mapping of the coarse grid (red) to finite elements is not included. Note that there is a gap between the two bold crosses in (d). This is because the refinement does not affect the elements that have already been assigned a value. This arrangement ensures the mapping from the density field back to the refinement crosses is one-to-one.

In this example two refinement levels are applied. A further refinement level would merge the edges and form fully solid blocks in the finite element grid. While such fully solid blocks cause no problem for the optimization, here we exclude this situation for the purpose of generating a uniform thickness inside the design domain. The number of supported refinement levels depends on the resolution within the coarse cell. For instance, increasing from the current within each coarse cell (red) to would allow two more refinement levels.

#### 2.1.1 Multi-level Design Variables ()

A design variable is assigned to each of the cells in the quadtree grid, with (resp. ) indicating a refinement (resp. non-refinement) for the cell with the 2D indices . Design variables are assigned on each level in the quadtree, denoted by , where indicates the -th refinement level. The initial coarse grid is denoted as the first level. From the first level downwards, the resolution of cells, which is also the resolution of the design field, is doubled on each dimension from the previous level. Denoting the resolution on the -th level by , and the resolution of the initial coarse grid by , the relation can be written as

(1) |

where is the maximum allowed refinement level. For a coarse block of elements, is calculated by

(2) |

When the refinement level reaches , it creates voids in the size of finite elements.

#### 2.1.2 Mapping ()

The mapping from the quadtree grid to the density field () has been illustrated in Figure 2. Solid finite elements can be mapped either from the design variables , , corresponding to cross-shaped sub-structures, or from the initial coarse cells (red). For compactness in the following formulation, let us denote the initial coarse cells by which has the same resolution as , i.e., . Note that has a constant value for all coarse cells. This corresponds to a coarse frame for ensuring the minimum stiffness of the optimized structure.

Reshaping the multi-level design fields , , and the density field into column vectors and , respectively, the mapping can be written compactly as a transformation by

(3) |

The transformation is a sparse matrix with a dimension of . The non-zero entries in have a constant value of , and concern the specific finite elements that will be assigned by refinements . Since the non-zero entries have a unit value, the transformation assigns the same value of the refinement to the affected finite elements. For instance, an intermediate refinement by leads to intermediate densities for the corresponding finite elements, .

### 2.2 Optimization Problem

With the design variables () and the physical density field () defined in the previous subsection, we can write the commonly-studied compliance minimization problem (e.g., in Sigmund01 (); Andreassen10 ()) as follows:

(4) | ||||

s.t. | (5) | |||

(6) | ||||

(7) |

The objective is to minimize the compliance represented by . The density vector is calculated from the design variable field by Eq. (3). The first constraint is the static equation of linear elasticity. , , and are stiffness matrix, displacement vector, and force vector, respectively. The global stiffness matrix is assembled from element stiffness matrix . The modified SIMP interpolation model Andreassen10 () is used to compute the element stiffness matrix:

(8) | ||||

(9) |

where is the element stiffness matrix for an element with unit Young’s modulus, is the stiffness of the material, is a very small stiffness assigned to void elements in order to prevent the stiffness matrix from becoming singular, and is a penalization factor (typically ) introduced to ensure black-and-white solutions. The second constraint measures the volume occupied by solid elements. is the unit volume for each element, and is the maximum allowed volume. The last constraint restricts the design variables to take a value between (non-refinement) and (refinement).

Solving this optimization problem for the cantilever beam presents the design shown in Figure 3. It generates many crossing sub-structures which correspond to the design variables. However, many of such sub-structures (pointed out by the red arrows) are suspended. This is due to the fact that the (imaginary) fine cells are refined, even though the fine cells themselves have not yet been created from refining their parent cells. The suspended sub-structures don’t effectively carry mechanical loads, and disobey the intention to create quadtree-like mechanical structures. In the following, a refinement filter is proposed to encode the recursive refinement rules, and thus to effectively filter out such suspended sub-structures.

### 2.3 The Refinement Filter ()

In creating a quadtree the refinement is only applicable to a cell that already exists. This implies dependency of the refinement of a fine cell on the refinement of its parent cell. In the binary setting, the refinement of a cell on the th level, , will not happen, as long as the refinement of its parent cell has not yet taken place, i.e., . Here refers to the indices of the parent cell. To encode this dependency in the setting of continuous design variables, we introduce a filtered design variable :

(10) |

In case of , i.e., the parent cell is not refined, the filtered value is

(11) |

This min function thus restricts the design variable on the fine level when the parent cell has not been refined. In case of , the filtered value is unaffected,

(12) |

Applying refinement recursively creates a quadtree with multiple levels. The dependency thus goes from the fine cell through all the intermediate levels to one of the coarse cells on the first level. For instance, consider the design variable illustrated in Figure 4. The filtered value is , where . At first glance this requests multiple filters executed sequentially, leading to much complex sensitive analysis afterwards. However, the recursive min functions can be consolidated into an equivalent min function by taking in as input all the involved arguments in the individual min functions, i.e.,

(13) |

This consolidation does not introduce approximate error, and significantly simplifies sensitivity analysis.

After the design variables have been filtered, the density vector, , is updated. This is done by replacing design vector by the filtered design vector in Eq. (3), i.e.,

(14) |

The effect of the recursive refinement filter is visualized in Figure 5. On the left, the multi-level design variables are initialized randomly with a value of or . A few suspended sub-structures can be observed on the left, and they disappear on the right. Figure 6 shows the result of integrating this filter into the optimization for the cantilever beam. The effect can be observed by comparing to Figure 3 which contains many suspended sub-structures.

#### 2.3.1 Smooth Approximation

Gradient-based numerical optimization necessitates differentiable functions. To this end, the non-differentiable min function (Eq. (13)) is approximated by a smooth -norm function, with a negative exponent,

(15) | ||||

(16) |

Here represents the -norm exponent, to be distinguished from the penalty in SIMP interpolation (Eq. (9)). As approaches negative infinity, the -norm equals the minimum value. The -norm is normalized by a factor considering the number of arguments, in this formulation. The normalization reduces the approximation error when a practical value is not infinite Wu17-infill (). In our tests, we use .

The indices of parent cells are needed in Eq. (16). When using a regular discretization, the indices can be found by a recursive function,

(17) |

where the function is defined as

(18) |

### 2.4 Sensitivity Analysis

By using the chain rule, the sensitivity of the objective function with respect to the design variable is given by:

(19) |

The first term is derived from adjoint analysis. The individual entry is

(20) |

The second term can be evaluated from Eq. (14),

(21) |

The third term comes from the smoothed refinement filter, Eq. (16). The non-zero entries (i.e., if ) are calculated by

(22) |

The sensitivity of the volume with respect to the design variable is derived in a similar way:

(23) |

Individual entries in the first term are given by

(24) |

where is the unit volume for each element.

## 3 Extensions

### 3.1 Balanced Quadtree

We provide an option to restrict the quadtree such that the level difference between neighbouring cells is at most one, known as balanced or restricted quadtree. A balanced quadtree maintains a smooth transition from fine cells to coarse cells. For 3D printed infill, a smooth transition is expected to make the part more stable for uncertain loads.

In an unbalanced quadtree, the refinement of a cell depends only on its parent cell (Eq. (10)). In the balanced quadtree, the refinement additionally depends on the neighbours of its parent cell. The refinement filter (Eq. (10)) is thus updated by

(25) |

Here, it is assumed that the indices belong to and belong to , otherwise the invalid arguments are excluded. A consolidation of this recursive function is obtained similar to Eq. (13). The recursive function can point to the same coarser cell multiple times, due to the introduction of neighbours. Only one copy of each dependent cell is included in the consolidated filter. Since the encoding of neighbourhood and hierarchy remains constant during the optimization, finding the dependency is only performed once in the initialization.

The effect of the balanced refinement filter is visualized in Figure 7, by comparing the unbalanced refinement (left) and the balanced refinement (right). An optimized balanced quadtree for the cantilever beam is shown in Figure 8.

### 3.2 Integration with other Filters

In density-based topology optimization using square elements, regulations such as density filter Guest04 () or sensitivity filter Sigmund01 () are often applied to get rid of checkerboard patterns (i.e., regions of alternating black and white elements). In our formulation, the refinement itself can be regarded as a regulation, and an additional density/sensitivity filter is no longer necessary.

In topology optimization with continuous design variables, projection filters are often applied to improve the convergence of approaching a binary white-or-black design. Our continuous reformulation of refinement also benefits from such a projection. In particular, we integrate the projection filter proposed in Wang10 ():

(26) |

Here controls the sharpness of projection. In the limit of the projection approaches a discontinuous step function. We employ a continuation of this parameter. is the threshold value, and we choose .

## 4 Numerical Examples

The method has been implemented based on the Matlab code provided by Andreassen et al. Andreassen10 (). The optimization problem is solved by using the method of moving asymptotes (MMA) Svanberg87 ().

In the following we report numerical results tested on regular design domains (Section 4.1) and curved design domains (section 4.2), and compare results from continuous and discrete optimization (Section 4.3). The following parameters are used in all tests: The SIMP penalization , the -norm , the threshold for projection . The sharpness of projection starts with and is increased up to 32, by doubling its value every iterations. The maximum number of iterations is .

### 4.1 Regular Domain

#### 4.1.1 Cantilever Beam

The first example is a cantilever beam that has appeared multiple times during the exposition of the method. The boundary condition is illustrated in Figure 3. The beam is fixed on the left side, while a force is applied on the right side in the middle. The prescribed material volume is of the design domain. The rectangular design domain is initialized with a coarse grid of . Each grid cell contains square elements, leading to a finite element resolution of for the entire domain. The maximum refinement level is .

Figure 10 shows a sequence of intermediate structures during the optimization of a balanced quadtree structure. As the optimization progresses, the quadtree structure emerges from the grey density distribution. To measure how close the continuous density field is to a binary field, the sharpness factor Sigmund07 () is defined as

(27) |

where is the number of finite elements. The sharpness factor becomes when the density values are converged to a strict 0-1 solution, and it becomes if all elements take a value of . The optimized quadtree (Figure 10 bottom right) has a sharpness factor of , which corresponds the field where density values on average take or .

The convergence plots are shown in Figure 10. The plot on the left shows that the compliance value decreases almost monotonically, except for a few jumps corresponding to the increase of in the Heaviside projection. The middle shows the volume fraction during the process. The design variables are initialized homogeneously with a value of such that the volume fraction of has been satisfied from the beginning. The plot on the right shows the convergence of sharpness. The optimization is terminated at the iteration of as the maximum change in the density field becomes smaller than .

#### 4.1.2 MBB Beam

The second example is the MBB design problem shown in Figure 11. Due to symmetry only one half of the design domain is optimized. The prescribed material volume is of the design domain. The half domain is initialized with a coarse grid of , where each grid cell is composed of elements, leading to a total of .

Figure 13 compares the quadtree structures obtained by using different refinement levels. The three rows from top to bottom show optimized quadtree using a maximum refinement level of three, four, and five, respectively. In each row, the left corresponds to the unbalanced refinement, while the right is obtained in the setting of balanced refinement. Comparing left and right, it can be found that the balanced quadtree (right) has a larger compliance and a larger sharpness value than in the unbalanced setting (left). This is due to the fact that the balanced refinement additionally involves the neighbours of the parent cells. Comparing the three rows, we observe smaller compliances when a larger refinement level is applied. This is because a larger refinement level increases the solution space. At a refinement level of five (bottom), small cavities from the previous refinement (middle) are filled with material. The optimizations with these different settings all converge well, as can be seen from the convergence plots in Figure 13.

### 4.2 Curved Design Domain

The adaptive structure refinement works as well on curved design domains. For a curved design domain, for simplicity in implementation we simulate the smallest axis-aligned bounding box enclosing the curved domain. The voids outside the curved domain are assigned as passive void elements, while a boundary layer with a thickness of 2 elements are passive solid elements. This requires updating the corresponding entries in the transformation matrix in Eq. (14). In the following we show two examples.

#### 4.2.1 Bracket

Figure 14 illustrates the dimension and boundary conditions of the bracket design problem. It contains curved boundary where a quadtree mesh is not conformal. The tight bounding box of the bracket is discretized using a coarse grid of , where each grid cell has elements, i.e., a total of elements for the bounding box. First we uniformly refine the coarse grid by two levels. The uniform infill is shown in Figure 1 (left). Its volume fraction is . Using this volume fraction as a constraint we perform adaptive structure optimization, allowing three levels of refinement. The optimized quadtree infill is shown in Figure 1 (right). Numerical analysis suggests that the adaptive infill is four times stiffer than the uniform infill.

Figure 15 further compares the structural performance using a different volume fraction. With one refinement level the uniform infill takes a volume fraction of and results in a compliance value of . Using the same volume fraction and allowing two levels of refinement the optimized infill has a compliance value of .