Adaptive Density Optimization of Lattice Structures Sustaining the External Multi-load

Abstract

In recent years, additive manufacturing have attracted increasing attention and promoted the development of lightweight modeling methods. Some studies have been carried out using internal filling structure to optimize the 3D model, while reducing the weight of the model, it can also satisfy some physical properties, such as to withstand the external loads. This paper presents an adaptive infilling structure optimization method based on the triply periodic minimal surface (TPMS), which can be represented by a pure mathematical implicit function. The morphology of these lattice structures within different stress regions can be adaptively adjusted, so as to reduce the weight of 3D printed objects while sustaining the external multi-load constraints. Firstly, finite element method is used to analyze the stress distribution of the original model infilled with uniform lattice structure. According to its stress value, the internal lattice structure is divided into three regions consists of high region (HR), transition region (TR) and low region (LR). Then, the inner structure within different stress regions is adaptively adjusted to relieve the stress concentration. Finally, we demonstrated that the proposed algorithm can effectively reduce the weight of 3D model while sustaining its mechanical strength.

1 Introduction

Additive manufacturing (AM) refers to fabricate the part layer by layer, it enables the manufacture of novel, highly complex geometric models possible [1]. To optimize a model’s physical attributes, such as strength, stability and weight [2,3,4,5]. A typical practice in AM is to use a uniform infilling structure, which allows a balance between the part strength and the amount of material. However, most of the optimized internal structures are prone to create large overhangs, which make the structure impossible to be directly manufactured by AM without deformation [6,7,8], so it is inevitable to generate the support structures to assist manufacturing [9]. When the processing is finished, the support structures need to be removed, which is a hard and cumbersome process. It not only increases the workload but also causes waste of materials. So, there is a great demand for the study of the internal filling structures with both self-supporting functions and lightweight features [10].

This paper puts forward an adaptive density optimization method for the self-supporting lattice infilling structures sustaining the external multi–load. The lattice structures are distributed in relation to a density function defined by the stress map, allowing simple control to trade off between the strength and material consumption for AM. The contributions of this method are as follows:

• We adapt three types of self-supporting lattice model, which derived from triply periodic minimal surfaces [11], to generate the internal filling structures. The morphology of the lattice structure can be effectively controlled, such as its size, horizontal angle and pillar thickness, which makes the design of the internal structure more flexible. Meanwhile, the self-supporting characteristic could reduce the material consumption and improve the forming efficiency [12].

• We introduce a material optimization algorithm with respect to the exterior multi-load. Based on the finite element analysis, the stress of the object can be computed and the initial lattice infilling structures will be divided into three regions according to its stress value. Then the morphology of lattice structures within different stress regions can be adaptively adjusted, which makes the optimization of the internal filling structures reasonable and executable.

• Our method allows a simple control to effectively reduce the weight of 3D model while sustaining its mechanical strength.

The rest of this paper is organized as follows. Section 2 illustrates the main methodology. Section 3 introduces the lattice infilling structures derived from TPMS and explores its structural features. Section 4 indicates a stress analysis method for the 3D model and an adaptive optimization method formula. Section 5 provides three examples to express the optimization process. Section 6 summarizes the conclusions of this paper and lists some ideas for future work.

2 Methodology

The basic perspective embedded within the optimization is increasing the volume of high stress region, decreasing the volume of the low stress region and progressively changing the rest of the region. At first, original 3D model preprocessing should be done, which contains three steps, (1) Using the contour offset method to generate a hollowed shell structure from the input 3D model [13]. (2) Using the marching cube (MC) algorithm to generate a uniform internal structure, then filling the minimum Axis-Aligned Bounding Box (AABB) [14] that belongs to the inner shell. (3) Using Boolean operation to fill the inner shell with homogeneous lattice structure (Figs. 1 and 2).

Fig. 1.

Flow chart of optimization process

Fig. 2.

A bead head optimized example

Then, the finite element analysis need to implement to calculate the stress distribution of the original model infilled with uniform lattice structure. Some conditions should be initialized including external force, constraint conditions, material properties and so on. This paper computes an initial stress map using public OOFEM finite element library [15].

Finally, by means of the stress diffusion method proposed in Sect. 4.1, the internal structure is divided into three regions and the density of the lattice structure within different stress regions is adaptively adjusted.

3 TPMS Internal Structure

A TPMS structure is a surface that is locally area-minimizing, that is, a small piece has the smallest possible area for a surface spanning the boundary of that piece. Wu et al. [3] studied the self-supporting internal structure, which showed that it is no need to add the auxiliary support structure when the direction of the structure is horizontally at a certain angle ($\ge 45^{^\circ }$). Lattice structures derived from TPMS in this paper also have the self-supporting characteristics. Figure 3 gives three types of TPMS lattice structures we used. Figure 3(a) is the Schwarz’ P-lattice structure, which meets the self-support requirement between the growth direction and the horizontal angle. Figure 3(b) is the Schoen’s G-lattice structure, which has zero mean curvature, i.e. the sum of the principal curvatures at each point is zero. The angle between the direction of stretching and the horizontal axis can also be adjusted. Figure 3(c) is called Schwarz’ D-lattice structure, whose skeletons have a diamond lattice pattern and its porous structure is about 45 degrees with horizontal direction.

Fig. 3.

Three types of TPMS lattice structure, (a) P-lattice, (b) G-lattice and (c) D-lattice

From the viewpoint of Yoo [16], the implicit surface functions of the three kinds of TPMS in Fig. 3 can be expressed as Eqs. (1), (2) and (3), respectively, where $X = {{2\pi x} \mathord{\left/ {\vphantom {{2\pi x} {L , { }Y = {{2\pi y} \mathord{\left/ {\vphantom {{2\pi y} L}} \right. \kern-0pt} L} , {\text{ Z = }}{{ 2\pi {\text{z}}} \mathord{\left/ {\vphantom {{ 2\pi {\text{z}}} L}} \right. \kern-0pt} L}}}} \right. \kern-0pt} {L , { }Y = {{2\pi y} \mathord{\left/ {\vphantom {{2\pi y} L}} \right. \kern-0pt} L} , {\text{ Z = }}{{ 2\pi {\text{z}}} \mathord{\left/ {\vphantom {{ 2\pi {\text{z}}} L}} \right. \kern-0pt} L}}}$, x, y, z refer to the Cartesian coordinates. Parameter L determines the edge length of the cubic unit cell, which will set as a constant value in this paper, parameter t determines the volume fraction of the regions that are separated by the surface [5, 17]. Based on the two parameters, this paper studies the relationship between the edge length and the volume fraction of the three kinds of TPMS lattice structures. As shown in Fig. 4, where $t \in [ - 0.5, \, 0.9]$, the vertical axis represents the percentage of TPMS structure to the same scale in a solid structure. Besides, there is a linear relationship between the t-value and the volume fraction. Equations (4, 5 and 6) illustrated the linear functions of the three kinds of surfaces according to t-value respectively. In this paper, they are called as density function which indicates the percentage of the lattice relative to the same size of a solid volume.

Fig. 4.

The relationship between the t-value of three TPMS structure and its volume fraction (Color figure online)

$$ F(X,Y,Z) = \cos (X) + \cos (Y) + \cos (Z) - t $$

(1)

$$ F(X,Y,Z) = \sin (X)\cos (Y) + \sin (Y)\cos (Z) + \sin (Z)\cos (X) - t $$

(2)

$$ \begin{array}{*{20}l} {F(X,Y,Z) = \sin (X)\sin (Y)\sin (Z) + \sin (X)\cos (Y)\cos (Z)} \hfill \\ {\quad \quad \quad \quad \quad + \cos (X)\sin (Y)\cos (Z) + \cos (X)\cos (Y)\sin (Z) - t} \hfill \\ \end{array} $$

(3)

$$ \left\{ {\begin{array}{*{20}c} {\rho_{P} (t) = 0.5 - 0.288t} \\ {s.t. \, - 0.5 \le t \le 0.9} \\ \end{array} } \right. $$

(4)

$$ \left\{ {\begin{array}{*{20}c} {\rho_{D} (t) = 0.5 - 0.418t} \\ {s.t. \, - 0.5 \le t \le 0.9} \\ \end{array} } \right. $$

(5)

$$ \left\{ {\begin{array}{*{20}c} {\rho_{G} (t) = 0.5 - 0.327t} \\ {s.t. \, - 0.5 \le t \le 0.9} \\ \end{array} } \right. $$

(6)

4 Adaptive Density Optimization of Lattice Structure

To balance the mechanical strength of the lattice structure and reduce the material consumption effectively, this paper proposed an adaptive optimization method which can change the local morphology of the lattice structure while sustaining the given multiple forces. Starting from the model infilled with uniform lattice structure, the finite element method will be implemented so to analyze the stress value of the initial input model. Then the stress analysis can be associated with the internal lattice volume. As presented in Sect. 3, the parameter t deciding the lattice volume, when it changing, the diameter of the lattice pillars will be changed. Based on this phenomenon, we can adjust the parameter t to change the morphology of the internal lattice structure, so as to meet the specific strength and weight requirements.

4.1 Lattice Stress Analysis

To facilitate the execution of finite element calculations, we use the TetGen library to generate the tetrahedral meshes firstly, then by means of OOFEM library to carry out the finite element analysis. As described in Sect. 2, we will calculate the stress distribution under multiple loads. For simplicity, we assume the direction of each load is perpendicular to each other and use the Von-Mises stress to indicate the stress of the object, as shown in Eq. (7), $\sigma_{xy} , { }\sigma_{yz} , \, \sigma_{zx}$ are x-axis principal stress, y-axis principal stress and z-axis principal stress, respectively. Similarly, $\tau_{xy} , { }\tau_{yz} , \, \tau_{zx}$ are XOY surface shear, YOZ surface shear and ZOX surface shear, respectively.

$$ \sigma = \sqrt {\frac{1}{2}[(\sigma_{xx} - \sigma_{yy} )^{2} + (\sigma_{yy} - \sigma_{zz} )^{2} + (\sigma_{zz} - \sigma_{xx} )^{2} + 6(\tau_{xy}^{2} + \tau_{yz}^{2} + \tau_{zx}^{2} )]} $$

(7)

Stress values of each tetrahedron in the model are obtained by the above formula. After dispersing the tetrahedrons into points, the stress regions of the model are then divided into HR, TR and LR according to the stress value of the points, as shown in Fig. 5. Among which, the high stress region and the low stress region are divided by the stress diffusion method proposed in this paper, and the remaining part is automatically classified as the transition region. Stress diffusion steps are as follows(using low region as an example to explain).

Fig. 5.

Stress interval extracting process

• 1. Sort the stress points of the model in ascending order, then the lowest stress (LS) will be found;

• 2. Search the stress points inside the lattice and take the top 10 points as initial points set $S_{0}$;

• 3. Define the scale value as SV and then collect the points whose stress value $\sigma \in [LS, \, SV * LS]$ as candidate set $S_{1}$, if $p_{j}^{1} = \left\{ {\exists p_{i}^{0} \in S_{0} |\left\| {p_{i}^{0} - p_{j}^{1} } \right\|^{2} \le D} \right\}$, where $p_{j}^{1} \in S_{1}$, D is distance threshold, $p_{j}^{1}$ will be added to point set $S_{0}$ to form a new set $S_{0}^{'}$. Points in set $S_{1}$ won’t be visited until the size of set $S_{0}$ is not increasing any more;

• 4. Use the length L to express the voxel unit, then construct a minimum AABB bounding box containing the model by this unit, where L value is related to parameter t in lattice;

• 5. Calculate the point number $N_{i}$ inside the voxel unit $i$, when $N_{i} \ge CN$ (CN is the default threshold), the boundary of voxel is marked as low stress region. The process will continue until all the voxels have been visited.

The division of the high stress region can be completed by using the same method, where the AABB follows the bounding box in the low stress step, searching for high stress regions, which ensures that the voxel cell size is same. Inside the bounding box, voxels which are not marked as high stress regions or low stress regions are automatically divided into transition stress region.

4.2 Lattice Density Optimization

After Sect. 4.1, the stress distribution of the model can be generated. From Sect. 3, we know the parameter t directly influences the volume fraction and the physical property of the structure. The material cost of the lattice structure will be reduced with the decrease of the volume fraction. In this paper, the relationship between the stress distribution and the volume fraction also can be established by the parameter t, so that the lattice structure can change the morphology characteristics of the lattice in different stress regions according to its stress value. Therefore, if the model need to meet strength and weight requirement at the same time, the objective function can be defined as follows,

Minimize

$$ E = \frac{{\lambda_{1} }}{2}u^{T} K(\rho_{TPMS} )u + \lambda_{2} V(\rho_{TPMS} ) $$

(8)

Subject to

$$ K(\rho_{TPMS}^{{}} )u = f $$

(9)

$$ V = \sum\limits_{i} {\rho_{TPMS}^{i} V_{solid} } \le V_{threshold} $$

(10)

$$ \rho_{TPMS} (t) = \left\{ {\begin{array}{*{20}c} {\rho_{P} (t)} \\ {\rho_{D} (t)} \\ {\rho_{G} (t)} \\ \end{array} } \right.{ - 0} . 5\le {\text{t}} \le 0. 9 $$

(11)

The first term $\frac{1}{2}u^{T} K(\rho_{TPMS} )u$ represents the strain energy. The second term $V(\rho_{TPMS} )$ is the volume of the internal structure, where $\lambda_{1,} \, \lambda_{2}$ are the balance factors during each iteration. Equation (9) describes the static state of the object under the given external force $f$, where $u$ means the displacement vector. Equation (10) represents the termination condition of the model volume, where $V_{threshold}$ is volume threshold, $\rho_{TPMS}^{i}$ is the density value of TPMS lattice of voxel unit $i$, $V_{solid}$ is a volume of unit cube. Equation (11) expressed the lattice density function which has three choices. In practical, we use the power-law relationship to compute the stiffness matrices $K(\rho_{TPMS} ) = \sum\nolimits_{e} {K_{e} (\rho_{TPMS} ) = (\rho_{TPMS}^{i} )}^{ - p} K_{0}$, where $K_{0}$ is the stiffness matrix of a solid element and $p$ is a penalization parameter [18].

5 Results and Discussion

To demonstrate the feasibility of the proposed method, several models are studied in this paper. In the work, the finite element analysis result of the model before and after optimization under the same load are compared, and the variables of the strength and volume of the interior structure in the iterative process are analyzed.

Strength Optimization.

The first model is a circular sandwich structure with a skin thickness of 1 mm, an internally filled D lattice structure. The external loads are the top of the model with four constraint points shown in Fig. 6(a). It shows the initial model and the stress of the interior lattice in Fig. 6(b) and (c). Figure 6(d) and (e) are optimized model and stress map of the interior lattice. Comparing the initial stress map, it can be found that the optimized lattice structure under the same load can effectively reduce the high stress region. In addition, we show the change of the high stress region and volume with the increase of iterations in the optimization process, as shown in Fig. 9(a), where high stress point (HSP) percent represents the percentage of the optimized high stress point and the initial high stress point (black square), and the volume ratio represents the percentage of the optimized lattice volume and the initial volume (red circle). It can be found that the proportion of high stress points decreases with the increase of iterations, which illustrates that the optimization process can reduce the high stress region and achieve the model strength optimization. At the end of iterations, the final HSP only has one third of the initial value. As the strength of the main optimization, in the optimization process, the black curve is located in the red curve below the trend.

Fig. 6.

Optimization of circular sandwich structure filled with D lattice structure

Volume Optimization.

Here is a model of the letter ‘E’ with a skin thickness of 1 mm, which is filled with P lattice structure. Figure 7(a) shows the external forces and constraint points. It shows the transparent map of the primitive model in Fig. 7(b). Figure 7(c) is the stress result of the interior lattice of initial model. Figure 7(d) and (e) are optimized model and its stress map of the internal lattice. Before and after optimization, the model morphology changed distinctly. Meanwhile, we show the change of the volume and high stress region with the increase of iterations in the optimization process, as showed in Fig. 9(b). It can be learned that the volume percentage gradually decreases and tends to be stable with the increase in the number of iterations, HSP percentage remains stable with the number of iterations increases. At the end of iteration, it achieves nearly 40% reduction in the volume. As the volume of the main optimization, in the optimization process, the red curve is located in black curve the below the trend.

Fig. 7.

Optimization of the model of the letter ‘E’ filled with P lattice structure

Balance Optimization.

In Fig. 8, it is a cuboid model filled with G lattice structure which has a skin thickness of 2 mm. Figure 8(a) to (d) show the morphological change of the lattice within the model during the optimization process. Figure 8(e) to (h) indicate the change of the stress during the optimization process. Combining both of them, it can be noted that the volume of the internal lattice which locates in the low stress region is obviously reduced. At the same time, the volume located in the high stress region is increased. Furthermore, we have calculated the high stress point percentage and the volume ratio of the model in the optimization process. It can be seen from the Fig. 9(c) that with the increase of the number of iterations, the volume ratio of the internal lattice of the model is obviously reduced, and the proportion of the high stress point is also decreased. At the end of iteration, the volume of optimized lattice is less than the original one-half, the HSP reduction is nearly 35%. As the balance optimization, both of the curves gradually decrease and tend to be stable with the increase of the number of iterations.

Fig. 8.

Optimization of cuboid model filled with G lattice structure

Fig. 9.

HSP percent and volume ratio of the three models, (a) circular sandwich structure, (b) the letter ‘E’, (c) cuboid model (Color figure online)

6 Conclusion

In this paper, we propose an adaptive density optimization method for the self-supporting lattice infilling structures sustaining the external multi–load. It divides the internal region into three parts according to the stress values which are computed by finite element method. In the different stress region, the morphology of the lattice structure can be adjusted adaptively, so the stress of the infilling structures will be relieve and the volume of the 3D model can be reduced. We demonstrated the effectiveness of the proposed method by several examples.

However, there are still some shortcomings in the paper. The optimized result need to be verified by multi-load mechanics experiment if the measuring instrument can apply multiple forces in a plurality of different directions. Therefore, new measuring method needs to be considered for the further work. In addition, although the method can calculate the stress distribution more accurately, the finite element method for the larger model will cost too much time. In the subsequent work, we will simplify the analysis of objects and optimize the analysis process to accelerate the overall optimization process.