Topology Optimization for FDM Parts Considering the Hybrid Deposition Path Pattern

Based on a solid orthotropic material with penalization (SOMP) and a double smoothing and projection (DSP) approach, this work proposes a methodology to find an optimal structure design which takes the hybrid deposition path (HDP) pattern and the anisotropic material properties into consideration. The optimized structure consists of a boundary layer and a substrate. The substrate domain is assumed to be filled with unidirectional zig-zag deposition paths and customized infill patterns, while the boundary is made by the contour offset deposition paths. This HDP is the most commonly employed path pattern for the fused deposition modeling (FDM) process. A critical derivative of the sensitivity analysis is presented in this paper, which ensures the optimality of the final design solutions. The effectiveness of the proposed method is validated through several 2D numerical examples.


Introduction
Additive manufacturing (AM) has gained fast development in research and industrial applications. The layer-by-layer material deposition nature of AM could enable graded material compositions and eliminate the design complexity constraints in conventional manufacturing methods [1]. AM breaks the restrictions between design and manufacturing and makes the greatest design freedom possible. Design for additive manufacturing (DfAM) has therefore attracted a great deal of attention [2,3].
However, there are new challenges introduced by AM [2,17], especially for the fused deposition modeling (FDM) process [18,19]. As is widely recognized, AM produces anisotropic material properties whose mechanical performance in the raster direction, the transverse direction, and the build direction are evidently different [20][21][22]. It means that the deposition direction of the material significantly affects the structural performance. Focusing on this point, build direction has recently been explored as an optimization variable to improve structural performance [23][24][25][26][27]. In reality, most fused deposition modeling (FDM) machines only support the hybrid deposition path (HDP) pattern, i.e., external contour attention has been drawn to designing a mechanical system while considering the effect of the material interface by topology optimization. Clausen et al. [33] proposed a double smoothing and projection (DSP) approach to design the coated structure with an enhanced solid shell and a weak base structure. Luo et al. [34] developed an erosion-based method to design shell-infill structures. Based on the level set method, a topology optimization method, which considers bi-material coated structures, was proposed in [35]. Besides, a new density filter was developed by Yoon et al. [36] to conduct topology optimization, considering the coating structure. More currently, using moving morphable sandwich bars, an explicit topology optimization method for coated structures was developed in [37]. Note that the stress constrained interface problem was also investigated by Yu et al. [38].
Inspired by the interface problem, a new algorithm based on solid orthotropic material with penalization (SOMP) and DSP approaches is developed to find out an optimal structure design which considers the HDP pattern and the anisotropic material properties. In this work, the external offset contour (Figure 1) is assumed to be a uniform boundary layer structure which has a different material from the substrate structure. In particular, the DSP method allows the identification of the boundary layer and achieves the length control for it; the anisotropic material topology optimization could be realized by the SOMP interposition model.  Currently, interface technology, like coating, is widely used in designing structures with enhanced functionalities and visual properties. To consider the interface issue from the view of structural optimization, other than simply designing the substrate structures, much research attention has been drawn to designing a mechanical system while considering the effect of the material interface by topology optimization. Clausen et al. [33] proposed a double smoothing and projection (DSP) approach to design the coated structure with an enhanced solid shell and a weak base structure. Luo et al. [34] developed an erosion-based method to design shell-infill structures. Based on the level set method, a topology optimization method, which considers bi-material coated structures, was proposed in [35]. Besides, a new density filter was developed by Yoon et al. [36] to conduct topology optimization, considering the coating structure. More currently, using moving morphable sandwich bars, an explicit topology optimization method for coated structures was developed in [37]. Note that the stress constrained interface problem was also investigated by Yu et al. [38].
Inspired by the interface problem, a new algorithm based on solid orthotropic material with penalization (SOMP) and DSP approaches is developed to find out an optimal structure design which considers the HDP pattern and the anisotropic material properties. In this work, the external offset contour ( Figure 1) is assumed to be a uniform boundary layer structure which has a different material from the substrate structure. In particular, the DSP method allows the identification of the boundary layer and achieves the length control for it; the anisotropic material topology optimization could be realized by the SOMP interposition model. The remainder of the paper is structured as follows: Section 2 presents the problem formulation. This includes the material model and the corresponding interpolation scheme, as well as the optimization problem and the sensitivity analysis. Section 3 presents the numerical implementation. Several numerical results are presented and discussed Section 4. Section 5 concludes the paper.

Problem Formulation
In this section, the optimization problem is defined. This includes defining an appropriate material model, formally defining the optimization problem, and deriving sensitivities. The material model and characteristic properties are derived analytically based on continuous versions of the design field and filters. The energy-based SOMP method is used to determine the structural topology; the DSP approach is adopted to distinguish the substrate structure and the boundary layer, whose implementation process is shown in Figure 2.
Micromachines 2020, 11, x FOR PEER REVIEW 3 of 23 The remainder of the paper is structured as follows: Section 2 presents the problem formulation. This includes the material model and the corresponding interpolation scheme, as well as the optimization problem and the sensitivity analysis. Section 3 presents the numerical implementation. Several numerical results are presented and discussed Section 4. Section 5 concludes the paper.

Problem Formulation
In this section, the optimization problem is defined. This includes defining an appropriate material model, formally defining the optimization problem, and deriving sensitivities. The material model and characteristic properties are derived analytically based on continuous versions of the design field and filters. The energy-based SOMP method is used to determine the structural topology; the DSP approach is adopted to distinguish the substrate structure and the boundary layer, whose implementation process is shown in Figure 2.

The Formulation for the Boundary Layer
According to the classical laminate theory, D 0 is the laminae unrotated compliance tensor [39,40]: and D(θ) is the 2D orthotropic elasticity tensor given any angle θ: with T(θ) being the transform matrix which is used to conduct the matrix coordinate transform: In this work, the local fiber orientation θ could be analytically expressed by Equation (3) and counted in the counter-clockwise direction, as shown in Figure 3.
∂ ϕ ∂y and ∂ ϕ ∂x are the spatial gradients of the filtered design field ϕ, and their derivations could refer to later content.

The Formulation for the Boundary Layer
According to the classical laminate theory, is the laminae unrotated compliance tensor [39,40]: and ( ) is the 2D orthotropic elasticity tensor given any angle : with ( ) being the transform matrix which is used to conduct the matrix coordinate transform: In this work, the local fiber orientation could be analytically expressed by Equation (3) and counted in the counter-clockwise direction, as shown in Figure 3.
and are the spatial gradients of the filtered design field , and their derivations could refer to later content. Thus, the element stiffness matrix could be written as: where ( * ) is the element stiffness matrix assembly operator.

The Optimization Model
In this paper, a standard compliance minimization problem, subject to mass constraint, is studied. The optimization problem is formulated as: where is the design variable and the Eth element density and is the total number of elements. , , and are the global stiffness matrix, displacement vector, and force vector, respectively. Thus, the element stiffness matrix could be written as: where δ( * ) is the element stiffness matrix assembly operator.

The Optimization Model
In this paper, a standard compliance minimization problem, subject to mass constraint, is studied. The optimization problem is formulated as: Minimize : C = U T KU Sub ject to :

Mass Constraint
In this paper, the mass constraint function could be written as: where G E is the mass of the element E: It is evident that the macroscale volume constraint considers both the substrate and coating materials. Like the notation form used in the last subsection, the following expression is introduced: where M 0 is the design element standard mass; S E , G E , and SG E are the three non-penalized pseudo-design fields, as follows:

Sensitivity Analysis
The updates of the design variables are performed based on sensitivity analysis using the Method of Moving Asymptote (MMA) algorithm, which requires first order sensitivity information of the constraints and the objective function. In this subsection, a critical derivative of the sensitivity analysis is presented.

Sensitivity Analysis for Objective Function
Recalling Equation (13), could be written in the following form: The first derivative term could be obtained using the chain rule: The term represents the standard modification of sensitivities due to projection and the detail expression can be found in [41].
is the standard modification of the smoothing filter and its derivative will be discussed in a later subsection.
The second derivative term of Equation (20) is given by: and the first derivative term could be written as the following form: in which the term indicates the standard modification of sensitivities due to the projection of the ∇ ϕ E α field. The derivative of the normalized gradient norm could be written as: in which: where N is a vector of the four shape functions relating nodal variable ϕ N with the elemental variable ϕ E ; B x and B y are the gradient computation matrices for N in the x and y directions, respectively, and are independent from the design variables. Therefore, we have: Micromachines 2020, 11, 709 8 of 22 Thus, and the following notation of Equation (28) is introduced: where For the term , it could have following form: where Substituting the material interpolation of Equation (2) into the term ∂D E ∂θ E will yield: The elastic tensor of matrix material D 0 is independent of µ E . Using the chain rule again, we could arrive at: where Then, the derivative of Similarly, the following notation of Equation (36) is introduced: All the derivative terms in Equation (41) could be obtained by the previous expressions.

Filtering Based on Helmholtz-Type Differential Equations
The smoothing filter adopted in this work could be implicitly represented by the solution of a Helmholtz-type partial differential equation (PDE) [42]: by imposing the homogeneous Neumann boundary conditions ( ∂ µ ∂n = 0) on the boundary of the design domain. The solution of Equation (42) can be written in a convolution integral form, which has a similar function to the classical filter [43]. In Equation (42), µ represents the unfiltered design field, and µ is the filtered field; the parameter R plays a similar role as the minimum filter radius (r min ) in the classical filter [43]. An approximate relation between the length scales for the classical filter and the PDE filter is given by [42]: Using the finite element method to discrete Equation (42): where K x is the standard stiffness matrix in finite element method for the scalar problem corresponding to R x and µ N is the representation of the filtered nodal field. Thus, the derivative of the filtered sensitives, with respect to the design variable, can be written as: Micromachines 2020, 11, 709 10 of 22 The element node representation of the field is obtained by: where T F is a matrix which maps the elemental values µ to a vector with nodal values µ N . Similarly, the element node representation of the filtered field could be expressed as: and According to the above derivation, the filtered sensitivities of the term ω S can be calculated as: the filtered sensitivities of the term ω G could be rewritten as: The terms ∂µ can be computed as: and Again, the filtered derivative for the term ω SG and the filtered sensitivity analysis for the mass constraint could be obtained in a similar way and are therefore omitted here.

Numerical Implementations
The proposed method is validated with several classical 2D benchmark cases in the next section. Four-node quadrilateral elements are adopted in all numerical examples. For the MMA optimizer, the default move limit is 0.3. Additionally, following the suggestion in Ref. [33], a continuation strategy for projection is adopted, where the sharpness factor for the substrate projection is set as β S = 1 at the beginning of optimization and gradually increased to 64 by doubling every 50 iterations (or at convergence), while the sharpness factor for the boundary layer projection is initialized with β G = 4 to ensure a sharp coating from the first iteration, and doubled every 50 iterations (or at convergence) until it is increased to 128. A projection threshold of 0.5 is used for µ S and µ G . The iterative process terminates when no further improvement in the objective function can be achieved, namely, when the difference in the objective values between two adjacent iterations is less than 0.01 or the maximum iterative number is exceeded. The whole process of the proposed method is shown in Figure 4.
doubling every 50 iterations (or at convergence), while the sharpness factor for the boundary layer projection is initialized with = 4 to ensure a sharp coating from the first iteration, and doubled every 50 iterations (or at convergence) until it is increased to 128. A projection threshold of 0.5 is used for and . The iterative process terminates when no further improvement in the objective function can be achieved, namely, when the difference in the objective values between two adjacent iterations is less than 0.01 or the maximum iterative number is exceeded. The whole process of the proposed method is shown in Figure 4.  The MBB beam problem is investigated to minimize the structural compliance under the maximum material volume ratio of 0.5, whose boundary condition is shown in Figure 5. The structural sizes are defined with L = 30 and H = 10. Only one half of the structure is optimized due to the symmetry condition. The MBB structure is loaded with a concentrated vertical force (F = 1) at the upper left corner; the bottom right corner is supported on a roller; and the asymmetrical boundary condition is applied to the left edge. The nodal displacement in the x-direction is restricted, while in the y-direction it is free. The MBB beam problem is investigated to minimize the structural compliance under the maximum material volume ratio of 0.5, whose boundary condition is shown in Figure 5. The structural sizes are defined with = 30 and = 10. Only one half of the structure is optimized due to the symmetry condition. The MBB structure is loaded with a concentrated vertical force ( = 1) at the upper left corner; the bottom right corner is supported on a roller; and the asymmetrical boundary condition is applied to the left edge. The nodal displacement in the x-direction is restricted, while in the y-direction it is free. The substrate material is assumed to be fully infilled ( = ) and the direction and the rotation angle of the raster direction is defined positively in the counter-clockwise direction (which is consistent with the depiction in Figure 3). A solid material with a Young's Modulus of 2.0 GPa in the raster direction and 0.5 GPa in the transverse direction is used. In addition, the Poisson's ratio is 0.4, and the shear modulus is 0.35 GPa. The substrate material is assumed to be fully infilled (ρ S = ρ B ) and the direction and the rotation angle θ of the raster direction is defined positively in the counter-clockwise direction (which is consistent with the depiction in Figure 3).
The boundary layer width T = 4 (R 1 = 14 and R 2 = 10) is investigated in the first test. The raster direction for the substrate material is assumed to be 0 • . In order to get a clear-cut solid structural design within the solid area, ρ E ≥ 0.95 indicates a clearly formed substrate domain, and ∇ ϕ E α ≥ 0.5 represents a clearly formed boundary layer. The final compliance is 88.3850, and the optimization terminates at the 320th iteration. The optimized result is shown in Figure 6. A solid material with a Young's Modulus of 2.0 GPa in the raster direction and 0.5 GPa in the transverse direction is used. In addition, the Poisson's ratio is 0.4, and the shear modulus is 0.35 GPa. The substrate material is assumed to be fully infilled ( = ) and the direction and the rotation angle of the raster direction is defined positively in the counter-clockwise direction (which is consistent with the depiction in Figure 3).
The boundary layer width = 4 ( = 14 and = 10) is investigated in the first test. The raster direction for the substrate material is assumed to be 0°. In order to get a clear-cut solid structural design within the solid area, ≥ 0.95 indicates a clearly formed substrate domain, and ‖ ‖ ≥ 0.5 represents a clearly formed boundary layer. The final compliance is 88.3850, and the optimization terminates at the 320th iteration. The optimized result is shown in Figure 6. From the convergence history graphs (Figure 7), it is seen that several sudden changes happen after the update of and , but it gradually becomes stable and finally converges. From the convergence history graphs (Figure 7), it is seen that several sudden changes happen after the update of β S and β G , but it gradually becomes stable and finally converges.
The detailed evolution of the topology of the substrate domain and boundary layer for the case in Figure 6 is given in Figure 8. As can be observed, the approximated structural topology is formed before the 200th iteration, and a clearer substate domain and boundary could be given in the last 120 iterations. Note that the boundary layer initially appears at the top left and bottom right corners, because the boundary layer (or solid substrate material) is required at all loads and helps to overrule the zero Dirichlet condition for the PDE filter [33]. The detailed evolution of the topology of the substrate domain and boundary layer for the case in Figure 6 is given in Figure 8. As can be observed, the approximated structural topology is formed before the 200th iteration, and a clearer substate domain and boundary could be given in the last 120 iterations. Note that the boundary layer initially appears at the top left and bottom right

Comparing with the Result from the Non-Boundary Layer Structure
In comparison, the structure without the boundary layer, under the same condition, is optimized in this test, and its optimized result is depicted in Figure 9. The raster direction of the substrate material is defined as 0°. To be specific, the filter radius ( = 15) is consistent with the case in Figure 6. Again, within the solid area, ≥ 0.95 indicates a clearly formed domain. It is clear to see that, under the same raster direction, the optimization result with the boundary layer ( = 88.3850) is better than the one without the boundary layer ( = 100.6204), and the compliance performance is 12.16% smaller than the latter one.

Comparing with the Result from the Non-Boundary Layer Structure
In comparison, the structure without the boundary layer, under the same condition, is optimized in this test, and its optimized result is depicted in Figure 9. The raster direction of the substrate material is defined as 0 • . To be specific, the filter radius (R 2 = 15) is consistent with the case in Figure 6. Again, within the solid area, ρ E ≥ 0.95 indicates a clearly formed domain. It is clear to see that, under the same raster direction, the optimization result with the boundary layer (c = 88.3850) is better than the one without the boundary layer (c = 100.6204), and the compliance performance is 12.16% smaller than the latter one.

The Influence of Different Raster Directions
In order to investigate the influence of different raster directions, this example explores the topology optimization with three designable raster directions (starting from 0°, 90°, and 45°) with the same boundary layer width ( = 4) and boundary condition. Correspondingly, the optimization results are demonstrated in Figure 10. The optimization results with the raster directions of 0°, 90°, and 45° have distinct structural topologies and shapes, and their structural performance is also different. The optimization result with the raster direction of 0° is better than the ones with the raster directions of 45° and 90°. This is reasonable from a mechanics point of view in sense that the principal stresses are distributed along the horizontal direction of the beam in the presented MBB problem.

The Influence of Different Raster Directions
In order to investigate the influence of different raster directions, this example explores the topology optimization with three designable raster directions (starting from 0 • , 90 • , and 45 • ) with the same boundary layer width (T = 4) and boundary condition. Correspondingly, the optimization results are demonstrated in Figure 10.

The Influence of Different Raster Directions
In order to investigate the influence of different raster directions, this example explores the topology optimization with three designable raster directions (starting from 0°, 90°, and 45°) with the same boundary layer width ( = 4) and boundary condition. Correspondingly, the optimization results are demonstrated in Figure 10. The optimization results with the raster directions of 0°, 90°, and 45° have distinct structural topologies and shapes, and their structural performance is also different. The optimization result with the raster direction of 0° is better than the ones with the raster directions of 45° and 90°. This is reasonable from a mechanics point of view in sense that the principal stresses are distributed along the horizontal direction of the beam in the presented MBB problem.

The Influence of Different Boundary Layer Widths
In Figure 11, the same design problem as above is solved with the raster direction θ = 0 • under the same conditions and varying boundary layer widths (T = 2, T = 4 and T = 6). The modeled boundary layer width is clearly controlled by modifying the filter radius R 2 , and the compliance improves when increasing the boundary layer width. In order to assure sufficiently wide features in the base structure, the first smoothing radius R 1 should be greater than or equal to the second smoothing radius R 2 [33]. Note that, in this work, the length control function could be implicitly achieved through the application of smoothing and projection filters. Increasing R 1 will lead to an increase in the minimum feature length and thereby eliminate some small geometry features.

The Influence of Different Boundary Layer Widths
In Figure 11, the same design problem as above is solved with the raster direction = 0° under the same conditions and varying boundary layer widths ( = 2, = 4 and = 6). The modeled boundary layer width is clearly controlled by modifying the filter radius , and the compliance improves when increasing the boundary layer width. In order to assure sufficiently wide features in the base structure, the first smoothing radius should be greater than or equal to the second smoothing radius [33]. Note that, in this work, the length control function could be implicitly achieved through the application of smoothing and projection filters. Increasing will lead to an increase in the minimum feature length and thereby eliminate some small geometry features.  Figure 12 shows the optimized structure discretized by three different element side lengths: 0.05, 0.1, and 0.2, respectively. The same topology and similar shapes could be found in the three final designs. The thickness of the boundary layer is almost independent of mesh size and highly uniform.  Figure 12 shows the optimized structure discretized by three different element side lengths: 0.05, 0.1, and 0.2, respectively. The same topology and similar shapes could be found in the three final designs. The thickness of the boundary layer is almost independent of mesh size and highly uniform.

The Fully Infilled Substrate Problem
Next, the optimization of the cantilever problem is conducted. The boundary conditions are presented in Figure 13; its left side is clamped and the middle point of the right side is loaded with a constant force (F = 1 ). The structural sizes are defined by L = 30 and H = 15. The mesh is discretized by 200 × 100 elements. The maximum material volume ratio is 0.5 in this case.

The Fully Infilled Substrate Problem
Next, the optimization of the cantilever problem is conducted. The boundary conditions are presented in Figure 13; its left side is clamped and the middle point of the right side is loaded with a constant force ( = 1 ). The structural sizes are defined by = 30 and = 15. The mesh is discretized by 200 × 100 elements. The maximum material volume ratio is 0.5 in this case. Three cases with fully infilled substrate materials under different raster directions (0°, 45°, and 90°) are given in Figure 14.

The Thick Boundary Problem
Finally, the structure with a thick boundary layer under the raster direction of 0° is investigated. The optimized objective value is 32.0089, which is lower than the result in Figure 14a. In order to guarantee a constant boundary layer thickness for the optimized structure, a relatively high minimum feature size is needed, and the boundary layer thickness should be much smaller than the feature size for the substrate structure. Besides, within our thick width results (Figure 15), we find the desired interface width dictated by the mesh resolution. For example, under the coarse mesh resolution (100 × 50) adopted in the same case in Figure 15, it is impossible to identify a boundary layer that is equal or above = 10.

The Thick Boundary Problem
Finally, the structure with a thick boundary layer under the raster direction of 0 • is investigated. The optimized objective value is 32.0089, which is lower than the result in Figure 14a. In order to guarantee a constant boundary layer thickness for the optimized structure, a relatively high minimum feature size is needed, and the boundary layer thickness should be much smaller than the feature size for the substrate structure. Besides, within our thick width results (Figure 15), we find the desired interface width dictated by the mesh resolution. For example, under the coarse mesh resolution (100 × 50) adopted in the same case in Figure 15, it is impossible to identify a boundary layer that is equal or above T = 10.
In order to guarantee a constant boundary layer thickness for the optimized structure, a relatively high minimum feature size is needed, and the boundary layer thickness should be much smaller than the feature size for the substrate structure. Besides, within our thick width results (Figure 15), we find the desired interface width dictated by the mesh resolution. For example, under the coarse mesh resolution (100 × 50) adopted in the same case in Figure 15, it is impossible to identify a boundary layer that is equal or above = 10.  Figure 16 shows the boundary layer raster direction distribution. The arrow indicates the raster direction for each boundary element. As can be seen in Figure 16, it is less prone to sudden orientation changes in the boundary layer domain, and it is matched well with the boundary layer structure.  Figure 16 shows the boundary layer raster direction distribution. The arrow indicates the raster direction for each boundary element. As can be seen in Figure 16, it is less prone to sudden orientation changes in the boundary layer domain, and it is matched well with the boundary layer structure.

The Fully Infilled Substrate Problem
The third numerical test is a short cantilever beam illustrated in Figure 17, where the left side is clamped and the middle point of the right side is loaded with a constant force ( = 1 ). The structural sizes are defined as = 15 and = 30. The whole design domain is meshed by 100 × 200 elements. Firstly, four fully infilled substrate materials with different raster directions (0°, 90°, and 45°) under a mass fraction constraint of 0.25 are considered in this subsection, and their optimized results are shown in Figure 18.

The Fully Infilled Substrate Problem
The third numerical test is a short cantilever beam illustrated in Figure 17, where the left side is clamped and the middle point of the right side is loaded with a constant force (F = 1 ). The structural sizes are defined as L = 15 and H = 30. The whole design domain is meshed by 100 × 200 elements.

The Fully Infilled Substrate Problem
The third numerical test is a short cantilever beam illustrated in Figure 17, where the left side is clamped and the middle point of the right side is loaded with a constant force ( = 1 ). The structural sizes are defined as = 15 and = 30. The whole design domain is meshed by 100 × 200 elements. Firstly, four fully infilled substrate materials with different raster directions (0°, 90°, and 45°) under a mass fraction constraint of 0.25 are considered in this subsection, and their optimized results are shown in Figure 18.

The Customized Infilled Pattern Problem
In this subsection, four customized infill patterns (wiggle, honeycomb, and two line infills) with 64.75%, 53.52%, 79.21%, and 79.10% density, respectively, are considered in this subsection. Their effective elastic tensors are predicted through the energy-based homogenization method [44,45], and the detail geometry structure, raster direction distribution, and effective elastic tensors are demonstrated in Figure 19. The problem configuration and the material properties are the same as the previous example, except the boundary layer width is = 3.

The Customized Infilled Pattern Problem
In this subsection, four customized infill patterns (wiggle, honeycomb, and two line infills) with 64.75%, 53.52%, 79.21%, and 79.10% density, respectively, are considered in this subsection. Their effective elastic tensors are predicted through the energy-based homogenization method [44,45], and the detail geometry structure, raster direction distribution, and effective elastic tensors are demonstrated in Figure 19. The problem configuration and the material properties are the same as the previous example, except the boundary layer width is T = 3.
The results are demonstrated in Figure 20   Finally, the computing time is briefly discussed. All the above cases were run on a desktop computer with Intel Xeon W-2145 CPU and 64GB RAM. Then, for the cases with the mesh dimension 300 × 100, an average of 18s is taken for each iteration: the FEM part takes 63.16%, the sensitivity analysis takes 18.42%, the MMA solver takes 15.79%, and the other parts take 2.63%. Meanwhile, for the cases with 200 × 100 elements, the algorithm takes 13s for each iteration on average: the FEM part takes 56.92%, the sensitivity analysis takes 17.12%, the MMA solver takes 19.23%, and the other parts take 6.73%. We could conclude that the FEM part takes more than half of the time in this method, and its computational cost grows rapidly with an increase in the dimension of the mesh. Therefore, a high efficiency FEM solver is demanded in this work.

Conclusions
The HDP pattern could be supported by most commercial tool path planning toolkits. Therefore, the HDP-based structure optimization could get closer to practice. Compared with the work proposed in [28] under the level set framework, based on solid SOMP and DSP approaches, this work proposes a methodology to find an optimal structure design, which takes the HDP pattern and the anisotropic material properties into consideration. The HDP pattern optimization here is assumed to be a structure optimization problem including coated structures, and the anisotropic material topology optimization is achieved by SOMP. The effectiveness of the proposed method are proven by several case studies, and the influence of different substrate raster directions under different boundary layer thicknesses is investigated. Note that the hybrid deposition paths produced in this work only provide the pattern where the zig-zag domain plays a significant role.
However, a unidirectional zig-zag deposition path is defined inside the substrate domain for the sake of simplicity. In fact, an optimized deposition path could achieve an even better design performance. The authors also intend to extend the proposed methods to address 3D problems. Besides, experimental validation is a must. These aspects will be explored in our future work as well.